JessicaHayward
commented
4 years ago Hi, I had the exact same problem - different genes being tagged with the same unique MSTRG identifier from stringtie --merge command. My work-around was to use this python script (https://gist.github.com/e2b3e89fc0220e4363765fadc8c1dd73), which modifies the merged.gtf output file by adding a suffix onto MSTRG identifiers for each different gene_name, and if there is no gene_name, then it adds a suffix onto MSTRG identifiers for each different transcript_id. Then I used this modified gtf file in stringtie -e command. This isn't a solution to the problem for novel isoforms but thought I'd put it here in case it is helpful to others.
m-waqas
I am interested to perform the differential gene expression (DE) and alternative splicing analysis of my RNA-seq dataset using Texudo2 approach as recommended in Pertea et al., 2016 (doi:10.1038/nprot.2016.095). I also want to use Deseq2 for DE analysis.
we have built a genome index using the command:
hisat2-build \ --ss Brassica_napus.AST_PRJEB5043_v1.43.ss \ --exon Brassica_napus.AST_PRJEB5043_v1.43.exons \ -f Brassica_napus.AST_PRJEB5043_v1.dna.toplevel.fa \ -p 4 \ napus
Then we aligned the reads using the following command: hisat2 -x genome/napus -1 fastq/${R1[$i]} -2 fastq/${R2[$i]} -S S${sample[$i]}.sam -p 4
Here I forget to use --dta parameter
Next, we sort and index the files:
samtools sort -@ 4 -o S${sample[$i]}.bam S${sample[$i]}.sam samtools index S${sample[$i]}.bam S${sample[$i]}.bai
The bam files were assembled using command: stringtie S${sample[$i]}.bam -G genome/Brassica_napus.AST_PRJEB5043_v1.43.gtf -o S${sample[$i]}.bam -l S${sample[$i]} -p 4
Next, We merged the assembles using stringtie --merge:
stringtie --merge -G genome/Brassica_napus.AST_PRJEB5043_v1.43.gtf -o Bnapus.gtf -p 4 merge_Bnapus_gtf.names
Merged assembly was compared using gffcomapre:
gffcompare -r genome/Brassica_napus.AST_PRJEB5043_v1.43.gtf -o Bnapus_compare Bnapus.gtf
This merged assembly was used to quantify expression using the command:
stringtie -e -B -G Bnapus.gtf -o ballgown/$sampleName/$sampleName.gtf -p 4 bam/$sampleName.bam
Now the problem is that I am getting the same gene id for different genes, I have checked different threads and noticed that this issue is raised several times (#159, #118, #153, #190, #217). It addressed well in some threads and is quite clear from this statement:
StringTie makes assembly decisions based primarily on the read alignment data. Reference annotation is often imperfect and lacking, and in order to allow for the discovery of novel isoforms, StringTie always uses the read alignments as the basis of transcript assembly. Unfortunately read alignments can also be wrong/imperfect and may actually "bridge" neighboring genes, as it seems to happen in our case where we are getting the same gene id for different genes.
I believe this issue can be sorted by three main approaches as mentioned in previous threads:
1) Using a better or more stringent read alignment strategy may help with this problem. Or some post-alignment filtering can be applied to the alignment data in order to eliminate large, low scoring alignments which seem to spuriously "connect" neighboring genes. This could be very challenging and time-consuming.
2) We can try to use stringent parameters for assembly. E.g. increasing value for -f (fraction of isoforms) -j (junction coverage) -c (coverage allowed for the predicted transcripts). This approach will help us up to some extent as some merged ref will be separated but there will still remain some merged together.
3) An alternate, faster differential expression analysis workflow can be pursued if there is no interest in novel isoforms (i.e. assembled transcripts present in the samples but missing from the reference annotation), or if only a well known set of transcripts of interest are targeted by the analysis. This simplified protocol has only 3 steps as it bypasses the individual assembly of each RNA-Seq sample and the "transcript merge" step to avoid the case that One internally generated locus IDs like MSTRG with several gene names.
Keeping in mind all these points, I still need your suggestions for some points to proceed further with downstream analysis.
1) I didn't use --dta/--downstream-transcriptome-assembly (Report alignments tailored for transcript assemblers including StringTie. With this option, HISAT2 requires longer anchor lengths for de novo discovery of splice sites. This leads to fewer alignments with short-anchors, which helps transcript assemblers improve significantly in computation and memory usage.) which can be helpful for the downstream analysis. So my question is that do I really need to use it and have to do the alignment again or am I good to proceed for downstream analysis?
2) Is it a good idea to proceed with just uniquely mapped reads, maybe this help to reduce the false-positive and assist to sort out the problem of same gene id for multiple genes
3) If I used the alternate approach (which has only 3 steps) can I still detect the alternative splicing using ASTALAVISTA?
4) If I keep my focus on just known genes and perform DE analysis based on it, can we still publish good quality stuff?
5) How about performing analysis using following two complementary approaches to quantify expression and then perform DE analysis based on both approaches:
a) Stringtie b) Salmon
Looking forward to hearing from you. please let me know if you need any further details.