SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 08-10-2010, 02:57 AM   #1
memo
Junior Member
 
Location: Oslo, Norway

Join Date: Aug 2010
Posts: 2
Question 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
memo is offline   Reply With Quote
Old 08-16-2010, 03:30 AM   #2
Uwe Appelt
Member
 
Location: Heidelberg, Germany

Join Date: Oct 2009
Posts: 27
Default

Quote:
Originally Posted by memo View Post
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!?
Uwe Appelt is offline   Reply With Quote
Old 08-16-2010, 03:48 AM   #3
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

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
Simon Anders is offline   Reply With Quote
Old 08-17-2010, 07:02 PM   #4
Cole Trapnell
Senior Member
 
Location: Boston, MA

Join Date: Nov 2008
Posts: 212
Default

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.
Cole Trapnell is offline   Reply With Quote
Old 01-25-2011, 09:32 AM   #5
lahoman
Member
 
Location: Houston

Join Date: Jan 2011
Posts: 12
Default

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
lahoman is offline   Reply With Quote
Old 01-25-2011, 09:49 AM   #6
jb2
Member
 
Location: Boston, MA

Join Date: Jun 2010
Posts: 25
Default

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.
jb2 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:35 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO