TOPHAT2 vs STAR and CUFFMERGE vs CUFFCOMPARE: CASE STUDY

[tonup69]

I've been busy. I'm still analyzing and re-analyzing some fly data that we are getting ready for publication and I found some interesting caveats I thought I would share with the SEQanswers community. Please keep in mind that all of these analyses were done on our local install of GALAXY - I don't have time for command line foolery and these are the same programs with the same options.

Basic Study Design:
UAS-Gene X (3 tech replicate)
Repo-GAL4>UAS-Gene X (3 replicates)
C155-GAL4>UAS-Gene X (3 replicates)

Samples were run on an Ion Proton machine here at the Molecular Resource Center.

Pipeline for Analysis: Tuxedo (Cuff-etc...)

Mapping:
I used TOPHAT2 for mapping with a reference sequence and GTF file for flies. Based on discussions here, I also used RNA-STAR. I got approximately 50% mapped reads (average 30M reads to start) using TOPHAT2 and 80% mapped reads with RNA-STAR. Sounds great, right? Well, along the way I found some things that are disturbing.

Cuffmerge vs Cuffcompare:
I realize these two programs are not exactly the same, but they should be similar, right? They both produce a reference GTF file for Cuffdiff to use to determine differential gene expression. Here are JUST the RNA-STAR mapped results using the classic-fpkm option.





Looks like Cuffmerge produced more significantly different genes per group across all groups, but the numbers are not crazy different.

Looking closer at just genes that were unique to Repo expression or C155 expression, I was surprised to find that there are rather large differences. in fact only a few genes are the same for both Cuffmerge and Cuffcompare.





There were bumps along the road here and things I would like to improve with this pipeline, but what is most disturbing is this huge difference in both mapping and GTF file production depending on which option you choose. I am going to try to use DeSeq and EdgR (which are not working on our install right now). May show the output for that as well just for the record.

Its a good example of "be careful what you believe" in bioinformatic analysis. This is why we do experiments to VALIDATE even RNAseq data.

[mastal]

You should probably use DESeq2 instead of DESeq.

I think the statisticians would probably say that the reason you get results like that is because even 3 replicates is not really a large sample.

[tonup69]

You should probably use DESeq2 instead of DESeq.

I think the statisticians would probably say that the reason you get results like that is because even 3 replicates is not really a large sample.

I would agree if these were human or mouse samples, but these are iosgenic fly stocks being used and all of the transgenes were placed into the same w1118 background. We usually get very tight statistics even for behavior studies in flies.

[tonup69]

Also, this is the SAME mapping data using Cuffdiff with the same settings in both cases. It is only Cuffcompare vs Cuffmerge that is different here.

[Brian Bushnell]

Well gosh, this sounds like a perfect case to bring in a third aligner as a tie-breaker! Not all RNA-seq aligners play nice with reads containing short indels, you know. But I happen to know that BBMap does

As far as the post-mapping analysis part goes... when tool X and tool Y give different results, I would place much more faith in the tool that does not start with "Cuff".

[gringer]

As far as the post-mapping analysis part goes... when tool X and tool Y give different results, I would place much more faith in the tool that does not start with "Cuff".

Indeed. This is especially so when the group that developed Cufflinks has released another program that supersedes it:

Nonetheless, StringTie consistently outperforms Cufflinks by a substantial amount, as shown below on four real data sets: GSM981256, GSM981244, GSM984609, and SRP041943

http://ccb.jhu.edu/software/stringtie/index.shtml?t=example#comp

The same applies to Tophat2:

HISAT2 is a successor to both HISAT and TopHat2. We recommend that HISAT and TopHat2 users switch to HISAT2.

http://ccb.jhu.edu/software/hisat2/index.shtml

It surprises me that you are aware of STAR, but not of HISAT2 and StringTie.

FWIW, the authors of DESeq have also indicated that DESeq2 is better:

Well, yes, it is intended. DESeq had a rather ugly hack to ensure type-I error control, which cost a lot of power, and fixing this was the core aim of the development of DESeq2.

http://seqanswers.com/forums/showthread.php?t=30907

[tonup69]

Other programs....

Well, as usual, analysis of sequence data is a moving target. Thanks for the other program options - I will look into these with our local Bioinformatics person who maintains our server.

That being said, I am not looking for perfection. I would like to find a set of significantly dysregulated transcripts which can be validated biologically. It is my feeling that even the Tuxedo pipeline should be able to do to that!

I'll try to use what I have now and see if any of this is true (i.e. these transcripts really change significantly by qRT-PCR). We will go from there.

LTR