Cufflinks2 analysis without significant genes

[rboettcher]

Dear all,

I'm currently analysing RNAseq data of two distinct cancer conditions using the tuxedo pipeline and I highly doubt my results so far. According to cuffdiff, there are absolutely NO genes found to be differentially expressed with a significant p-value which I cannot truly believe. When evaluating my results with cummeRbund, more than 60% of the transcripts are flagged with test status other than "OK".

Therefore, I would like to ask whether anyone has encountered similar problems and how to circumvent them. Since I already heard that cufflinks seems to have serious flaws, I run scripture in parallel although I still have to find tools to further analyse the assembled transcripts's abundance.

So here some details on my analysis:

• paired end RNAseq reads with 101 bases length

• mapping done via bowtie2/tophat2 against Hg19

• using latest releases of bowtie2, tophat2 and cufflinks2

here's the commands I use for tuxedo:

• tophat --b2-sensitive -p 8 -o mapping1 bowtie2-2.0.0-beta6/hg19/hg19 R1_001.fastq R2_001.fastq

• cufflinks --GTF UCSC_annotation_hg19.gtf -p 8 -o annotation1 mapping1/accepted_hits.bam

• cuffmerge -g UCSC_annotation_hg19.gtf -s bowtie2-2.0.0-beta6/hg19/hg19.fa -p 8 assemblies.txt

• cuffdiff -m 199 -s 38 -o diff_out -b bowtie2-2.0.0-beta6/hg19/hg19.fa -p 9 -L G1,G2 -u merged_asm/merged.gtf mapping1/accepted_hits.bam,mapping2/accepted_hits.bam,mapping3/accepted_hits.bam mapping4/accepted_hits.bam,mapping5/accepted_hits.bam,mapping6/accepted_hits.bam

note: I currently use the --GTF flag for cufflinks, though I will change that to --GTF-guide as soon as the results make sense.

and for scripture:

• for i in {1..22};
  do java -jar scripture-beta2.jar -alignment merged_align.bam -out chr$i.scriptureESTest.segments -sizeFile Hg19/chr$i\_size.txt -chr chr$i -chrSequence Hg19/Homo_sapiens.GRCh37.67.dna_rm.chromosome.$i.fa;
  done

Any help and suggestions are greatly appreciated!
Best regards

[sdriscoll]

you aren't alone in thinking the cufflinks/cuffdiff output doesn't make sense. i've never been happy with it either with novel isoform discovery and quantification or -G quantifaction only mode. if you're not after isoform level DE then you can go with raw read counts at the gene level and then use an R package like DESeq or EBSeq to analyze for differential expression. you'll need something to count your reads (you can try HTSeq, authored by the author of DESeq). HTSeq is a python library that comes with a script called htseq-count.

even if you ARE after isoform level DE it might be helpful to at least start at the gene/locus level using this alternative approach.

another possible pipeline is RSEM / EBseq. RSEM works by quantifying isoform level expression against a known annotation (any GTF). EBSeq is a new R package that takes isoform level uncertainty into account when calculating differential expression. I have not yet run all the way through this pipeline yet so I can't say how the output looks. of course with this pipeline you lose novel isoform discovery. i'd think you could take the output of Scripture and build a GTF which you could then use with RSEM for quantification.

[hlwright]

Yes, we have been experiencing similar problems (along with other problems) with the latest version of Cufflinks and are considering other methods for our DE analysis. I have some variance in my biological replicates which may be an issue for cuffdiff, but I have data for example where the RPKMs for Disease v Control are 250 vs 10 and this is not significant.

Here are some of the issues we have been trying to deal with. I have also emailed the cufflinks helpdesk but have not had a response.

http://seqanswers.com/forums/showthread.php?t=23962
http://seqanswers.com/forums/showthread.php?t=23910
http://seqanswers.com/forums/showthread.php?t=21708

[rboettcher]

Nice to see that I am not the only one experiencing this (though I did already find a lot of related threads).

I was hoping to get some replies by the cufflinks team as well, but so far it seems that they are busy preparing the next release version.

Despite that, my current plan is to switch to easyRNASeq for read counting and edgeR as well as DEXSeq for downstream analysis of my samples.
Thus far, transcript assembly is only a future prospective, but it will definitely end up on my table and scripture might be an option once I get it working properly. Unfortunately, the last time I tried to run it, scripture did not proof to be very user friendly.