Cufflinks/Cuffdiff significant differential expression

memo · 08-10-2010, 03:57 AM

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 · 08-16-2010, 04:30 AM

Quote:

Originally Posted by 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

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 · 08-16-2010, 04:48 AM

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 · 08-17-2010, 08:02 PM

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 · 01-25-2011, 10:32 AM

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 · 01-25-2011, 10:49 AM

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.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
CuffDiff 0 FPKM Output/Incorrectly Identified Differential Expression?	cw11	Bioinformatics	4	04-27-2015 04:59 AM
New differential testing of cuffdiff/cufflinks since 1.3.0	DerSeb	Bioinformatics	46	09-12-2013 06:46 AM
Cufflinks differential expression problem	Rachelly	Bioinformatics	2	05-12-2011 12:08 AM
Does Cuffdiff normalize FPKMS before it does differential expression?	greener	Bioinformatics	0	02-25-2011 03:37 PM
Cuffdiff differential expression two groups	jamminbeh	Bioinformatics	0	12-09-2010 03:13 PM

08-10-2010, 03:57 AM	#1
memo Junior Member Location: Oslo, Norway Join Date: Aug 2010 Posts: 2	Cufflinks/Cuffdiff significant differential expression 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

08-16-2010, 04:48 AM	#3
Simon Anders Senior Member Location: Heidelberg, Germany Join Date: Feb 2010 Posts: 994	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

08-17-2010, 08:02 PM	#4
Cole Trapnell Senior Member Location: Boston, MA Join Date: Nov 2008 Posts: 212	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.

01-25-2011, 10:32 AM	#5
lahoman Member Location: Houston Join Date: Jan 2011 Posts: 12	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

01-25-2011, 10:49 AM	#6
jb2 Member Location: Boston, MA Join Date: Jun 2010 Posts: 25	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.