Cufflinks/Cuffdiff significant differential expression

[memo]

Hi!

I'm new to RNAseq, and I'm trying to find genes/isoforms that are differentially expressed in two samples (one from wild type mouse and one from knock out). I have no replicates.

To do this I used SpliceMap, then Cufflinks/Cuffcompare/Cuffdiff.

The problem is that about 40% of the genes, 30% of the isoforms and practically all of the entries in 0_1_tss_group_exp.diff are considered significantly differentially expressed. This is obviously not logical, so I have concluded that something must be wrong. I just don't know what.

Is there anyone out there who can help me shed some ligth on where my mistake is?

These are the options used:

cufflinks -o data/KO_cufflinks -m 78 -s 35 -p 8 good_hits_KO.sam
cufflinks -o data/WT_cufflinks -m 78 -s 35 -p 8 good_hits_WT.sam

cuffcompare -o ccKOogWT -r mm9ch.gtf -R data/transcriptsKO.gtf data/transcriptsWT.gtf


cuffdiff -m 80 -p 8 -o data/cuffdiff ccKOogWT.combined.gtf good_hits_KO.sam good_hits_WT.sam

[Uwe Appelt]

Hi!

I'm new to RNAseq, and I'm trying to find genes/isoforms that are differentially expressed in two samples (one from wild type mouse and one from knock out). I have no replicates.

To do this I used SpliceMap, then Cufflinks/Cuffcompare/Cuffdiff.

The problem is that about 40% of the genes, 30% of the isoforms and practically all of the entries in 0_1_tss_group_exp.diff are considered significantly differentially expressed. This is obviously not logical, so I have concluded that something must be wrong. I just don't know what.

Is there anyone out there who can help me shed some ligth on where my mistake is?

These are the options used:

cufflinks -o data/KO_cufflinks -m 78 -s 35 -p 8 good_hits_KO.sam
cufflinks -o data/WT_cufflinks -m 78 -s 35 -p 8 good_hits_WT.sam

cuffcompare -o ccKOogWT -r mm9ch.gtf -R data/transcriptsKO.gtf data/transcriptsWT.gtf


cuffdiff -m 80 -p 8 -o data/cuffdiff ccKOogWT.combined.gtf good_hits_KO.sam good_hits_WT.sam

Hi Memo,

this actually does appear logical to me, because Cole described in his paper (http://dx.doi.org/10.1038/nbt.1621, see online-Methods, section "Analysis of gene expression and regulation dynamics.") that significance is based on a t-test and Benjamini-Hochberg correction for multiple testing.
Thus, the true accuracies of the various transriptome sequencing approaches aren't taken into account, despite accuracies obviously differ from experiment to experiment (e.g. depend on parameters like actual sequencing coverage produced). After all, there's no "simple" way to provide a empirical FDR, but you could certainly establish one of your own, by for example deriving it from (technical or biological) replicates.

Best, Uwe

ps: especially for isoform and gene abundance estimations it would also make sense to notice that genes or isoforms are constantly reported differentially expressed, if expressed in just either sample (XOR) no matter of the estimated expression magnitude. Cole has already clarified this to be an inconsistency (i unfortunately don't remember the thread). Abundance thresholding might hold as a quick fix here!?

[Simon Anders]

Hi memo,

what Uwe said was also my first guess. Just to double-check, have a look at the intensity distribution of your hits. The (relative) accuracy of the estimation of transcript molecule concentration improves with count rate, i.e., for strongly expressed genes, even small differences will appear significant. And, of course, there will always be small differences between two samples. Hence, I assume that above a certain expression threshold, nearly all genes are in your hit list.

Our software, DESeq, and Robinson et al.'s edgeR, are meant to address this issue. They estimate the biological variability from the differences between replicates and so can tell you whether an observed difference is significantly stronger than what you would expect as variation even within samples with the same genotype.

Of course, as you don't have replicates, you are screwed. There is no sound way to guess the biological variability without replicates.

Simon

[Cole Trapnell]

Hi memo,

You might want to also take a look at this paper by Bullard et al (http://www.biomedcentral.com/1471-2105/11/94/abstract). The authors make the point that genuine differences in the highly expressed genes between samples will shift the values of genes and transcripts lower in the expression profile, creating differences that appear to be significant but are really just artifacts of normalization. In my experience, these are attributable to a very small handful of loci in the genome - typically ribosomal or mitochondrial regions. Because rRNA is so abundant, and because library prep typically involves a polyA enrichment or an rRNA depletion step which itself has a pretty variable efficiency, I often mask out alignments from rRNA repeats in the genome as well as any reads that align to chrM, before running Cufflinks. You might also see this in certain tissues such as liver, which have a handful of genes that are drastically higher than nearly everything else.

[lahoman]

Hi, Uwe Appelt,

May I ask you a question? How did you choose the parameter of -p 8 for cufflink/cuffdiff? In other words, why did you use -p 8?

cufflinks -o data/KO_cufflinks -m 78 -s 35 -p 8 good_hits_KO.sam

cuffdiff -m 80 -p 8 -o data/cuffdiff ccKOogWT.combined.gtf good_hits_KO.sam good_hits_WT.sam

Thank you so much,

Lahoman

[jb2]

Hi Lahoman,

To my knowledge the -p argument specifies the number of threads to use. This is based on using computers/clusters with multiple cores/processors. If you are running a dual core computer, you could run Cufflinks in parallel with 2 threads instead of just one, speeding up the time it takes Cufflinks to run. You would select this based on the hardware you are using and the number of threads/processors available when you submit the Cufflinks job.