<p>Hi my RNAseq workflow is ending up with wrongly formatted gene IDs (and separately transcript _IDs) after a Hisat2 -&gt;samtools sort -&gt; stringtie -&gt; DEprep.py workflow.
DEprep.py outputs a gene_counts.csv which includes this...</p>
<pre><code>gene:ENSRNOG00000000017         47
ENSRNOG00000000021              220
</code></pre>
<p>Spelling this out &quot;gene:&quot; should not be in the gene_ID.  First thought of course was to just delete the &quot;gene:&quot; string with search and replace.
However, sometimes, the same gene appears as both indicating the error could be really damaging:</p>
<pre><code>Gene_ID                   Counts 
ENSRNOG00000004937        5234
gene:ENSRNOG00000004937   0
</code></pre>
<p>Digging around I have found that the issue is present in the initial stringtie /ballgown/ GTF output prior to running DEprep.py</p>
<pre><code>chr1    ensembl exon    la la   gene_id 
 &quot;gene:ENSRNOG00000046319&quot; la la la
</code></pre>
<p>Curiously my workflow is copied from my own local Galaxy installation workflow that uses the same steps together with the SAME fastq's the same, Hisat2 index, the same Guide.GFF, but there is no issue there.</p>
<p>I have spent a day googling this and no success yet, lots of peep have missing ids. One other person had this symptom, but sadly no resolution.  Any suggestions gratefully received!!
Code details below and an attached folder in the comments with dummy data and r script that reproduces this error.
HiSat2:</p>
<blockquote>
<p>hisat2 --dta --add-chrname --time --summary-file YH2minihisat.txt
--rna-strandness FR  --pen-cansplice 0 --pen-noncansplice 12 --pen-canintronlen G,-8.0,1.0 --pen-noncanintronlen G,-8.0,1.0 --min-intronlen 20 --max-intronlen 500000 --threads 12 -x rn6/genome -q -1 YH2mini_R1.fastq.gz -2 YH2mini_R2.fastq.gz |  samtools sort -m20G -l 0 -T /Users/me/Documents/TEMP/b -O bam | samtools view -O bam -@ 12 -o YH2mini.bam</p>
</blockquote>
<p>Stringtie:</p>
<blockquote>
<p>stringtie -fr -p 12 -G rat6.gtf -e -b ballgown/YH2mini -f 0.15 -m 200
-a 10 -j 1 -c 2 -g 50 -M 0.95 -o ballgown/YH2mini/YH2mini.gtf YH2mini.bam</p>
</blockquote>
<p>PrepDE:</p>
<blockquote>
<p>prepDE.py -i ballgown -g YH2mini_gene_counts.csv -t
YH2mini_transcript_counts.csv -l 100</p>
</blockquote>
<p>So peep can see I'm not lazy, I have a work around where I can run this, extra prior to prepDE.py, but seems very wrong!!</p>
<pre><code>if(CORRECTGTF){
    
        command &lt;- paste0(&quot;sed 's/gene\\://g' ballgown/&quot;,PREFIX,&quot;/&quot;,PREFIX,&quot;.gtf &gt; temp.txt&quot;) 
        system(command)
        command &lt;- paste0(&quot;sed 's/transcript\\://g' temp.txt &gt; ballgown/&quot;,PREFIX,&quot;/&quot;,PREFIX,&quot;.gtf&quot;) 
        system(command)
        command &lt;- paste0(&quot;sed 's/gene\\://g' ballgown/&quot;,PREFIX,&quot;/&quot;,&quot;t_data.ctab &gt; t_datafix.ctab.tmp&quot;)
        system(command)
        command &lt;- paste0(&quot;sed 's/transcript\\://g ' t_datafix.ctab.tmp &gt; ballgown/&quot;,PREFIX,&quot;/&quot;,&quot;t_data.ctab&quot;)
        system(command)
    }
</code></pre>

CC BY-SA 4.0
2022-10-13T10:22:36.387Z
The `gene:ENSRNOG00000046319` in the GTF file is pretty standard. I don't think we can help without access to the actual data. Can you give us some examples? for instance, you seem to be searching through something for the gene names. What are you searching through? Presumably the output of your pipeline, but what kind of output is that? Then, you say that you sometimes see both but you show a count of `0` for `gene:ENSRNOG00000004937`. My guess is that you are getting the gene names from different fields, but I would need to understand what you're doing in more detail to be sure.
terdon
CC BY-SA 4.0
2022-10-13T11:07:29.323Z
@ Terdon Thanks, I will just need to investigate now how to share data and upload anything needed (after a meeting now!).
RichardBJ
CC BY-SA 4.0
2022-10-13T11:26:12.427Z
We don't need real data, per se, just some examples to understand what you mean.
terdon
CC BY-SA 4.0
2022-10-13T17:56:27.767Z
...OK, sorry... don't really know how to give examples beyond the one above and so the actual code and data are here... https://www.dropbox.com/sh/lqk0usgt8n2suvx/AADqZDNpKKPKFGWZQ9sn9ix6a?dl=0
RichardBJ
<rna-seq><gtf><gff3><gff><stringtie>
Stringtie/DEprep.py gene/transcript IDs are wrongly formatted
# Stringtie/DEprep.py gene/transcript IDs are wrongly formatted

(Asked by: username_0 on 2022-10-13)

Hi my RNAseq workflow is ending up with wrongly formatted gene IDs (and separately transcript _IDs) after a Hisat2 ->samtools sort -> stringtie -> DEprep.py workflow.
DEprep.py outputs a gene_counts.csv which includes this...
gene:ENSRNOG00000000017         47
ENSRNOG00000000021              220

Spelling this out "gene:" should not be in the gene_ID.  First thought of course was to just delete the "gene:" string with search and replace.
However, sometimes, the same gene appears as both indicating the error could be really damaging:
Gene_ID                   Counts 
ENSRNOG00000004937        5234
gene:ENSRNOG00000004937   0

Digging around I have found that the issue is present in the initial stringtie /ballgown/ GTF output prior to running DEprep.py
chr1    ensembl exon    la la   gene_id 
 "gene:ENSRNOG00000046319" la la la

Curiously my workflow is copied from my own local Galaxy installation workflow that uses the same steps together with the SAME fastq's the same, Hisat2 index, the same Guide.GFF, but there is no issue there.
I have spent a day googling this and no success yet, lots of peep have missing ids. One other person had this symptom, but sadly no resolution.  Any suggestions gratefully received!!
Code details below and an attached folder in the comments with dummy data and r script that reproduces this error.
HiSat2:

hisat2 --dta --add-chrname --time --summary-file YH2minihisat.txt
--rna-strandness FR  --pen-cansplice 0 --pen-noncansplice 12 --pen-canintronlen G,-8.0,1.0 --pen-noncanintronlen G,-8.0,1.0 --min-intronlen 20 --max-intronlen 500000 --threads 12 -x rn6/genome -q -1 YH2mini_R1.fastq.gz -2 YH2mini_R2.fastq.gz |  samtools sort -m20G -l 0 -T /Users/me/Documents/TEMP/b -O bam | samtools view -O bam -@ 12 -o YH2mini.bam

Stringtie:

stringtie -fr -p 12 -G rat6.gtf -e -b ballgown/YH2mini -f 0.15 -m 200
-a 10 -j 1 -c 2 -g 50 -M 0.95 -o ballgown/YH2mini/YH2mini.gtf YH2mini.bam

PrepDE:

prepDE.py -i ballgown -g YH2mini_gene_counts.csv -t
YH2mini_transcript_counts.csv -l 100

So peep can see I'm not lazy, I have a work around where I can run this, extra prior to prepDE.py, but seems very wrong!!
if(CORRECTGTF){
    
        command <- paste0("sed 's/gene\\://g' ballgown/",PREFIX,"/",PREFIX,".gtf > temp.txt") 
        system(command)
        command <- paste0("sed 's/transcript\\://g' temp.txt > ballgown/",PREFIX,"/",PREFIX,".gtf") 
        system(command)
        command <- paste0("sed 's/gene\\://g' ballgown/",PREFIX,"/","t_data.ctab > t_datafix.ctab.tmp")
        system(command)
        command <- paste0("sed 's/transcript\\://g ' t_datafix.ctab.tmp > ballgown/",PREFIX,"/","t_data.ctab")
        system(command)
    }

**Tags:** rna-seq, gtf, gff3, gff, stringtie
> **Comments:**
> **username_1**: The `gene:ENSRNOG00000046319` in the GTF file is pretty standard. I don't think we can help without access to the actual data. Can you give us some examples? for instance, you seem to be searching through something for the gene names. What are you searching through? Presumably the output of your pipeline, but what kind of output is that? Then, you say that you sometimes see both but you show a count of `0` for `gene:ENSRNOG00000004937`. My guess is that you are getting the gene names from different fields, but I would need to understand what you're doing in more detail to be sure.  
> **username_0**: @ Terdon Thanks, I will just need to investigate now how to share data and upload anything needed (after a meeting now!).  
> **username_1**: We don't need real data, per se, just some examples to understand what you mean.  
> **username_0**: ...OK, sorry... don't really know how to give examples beyond the one above and so the actual code and data are here... https://www.dropbox.com/sh/lqk0usgt8n2suvx/AADqZDNpKKPKFGWZQ9sn9ix6a?dl=0  

---

## 0 answers

