
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
infinite fold change RNAseq  biofreak  Bioinformatics  6  07022012 10:02 AM 
DESeq DiffExpress Vs. Fold Change  KellerMac  Bioinformatics  3  06102011 10:49 AM 
fold change valuecuffdiff  madsaan  Bioinformatics  4  02102011 06:51 AM 
computing FOLD CHANGES  bogdan  RNA Sequencing  1  01312011 01:50 PM 
fold change for genes with 0 reads  rnaseq  General  2  11042010 07:44 PM 

Thread Tools 
11032010, 11:36 AM  #1 
Junior Member
Location: washington Join Date: Dec 2009
Posts: 7

fold change for genes with 0 reads
I am comparing Illumina RNAseq data sets for differentially expressed genes. A significant number of genes have 0 reads in one of the two conditions. What is the best way to include these genes? Is it acceptable to subsitute the 0 in the numerator or denominator with another number to get an approximate idea of the fold change?
Thanks! 
11032010, 01:13 PM  #2 
Member
Location: New Haven Join Date: Jul 2008
Posts: 36

you should be careful of these genes. In my points, you do not need calculate the fold change. You can split these cases into two situations: one condition is larger or smaller than threshold, e.g. gene RPKM>=5 (one Nature paper uses this scale). For the smaller, it is nothing, while the larger is significant different.

11032010, 01:46 PM  #3 
Junior Member
Location: washington Join Date: Dec 2009
Posts: 7

Thanks for the reply lmf_bill!
I am doing a timecourse comparison and typically the raw reads under one condition are 0 and in another condition are about 300 or higher so I know that there is at least a 300 fold increase in expression which is significant. I was wondering if there was a way to report this (approximate) numerical value rather than binning it as larger or smaller. I appreciate any additional feedback! 
11032010, 03:51 PM  #4 
Member
Location: Davis, CA Join Date: Sep 2010
Posts: 45

Hello all,
I am also very interested in everyone thinks about this as well. I've been debating on changing the zeros to ones to get some usable output, especially in cases where the denominator is zero (but not the numerator) because the output is undefined. Any suggestions? Cheers, Johnathon 
11032010, 04:06 PM  #5 
Member
Location: sweden Join Date: Aug 2010
Posts: 23

What about zeroinflated count models?
I once checked the histograms of many variables and many did not look like follow Poisson distribution or negative binomial distribution, as there were excessive number of zoers. 
11032010, 04:27 PM  #6 
Member
Location: Davis, CA Join Date: Sep 2010
Posts: 45

Hello Zhengz,
Thank you for your reply. I think part of the problem with my approach to this matter is my general lack of statistical savvy. Could you perhaps provide a good starting point or primer (or a useful open source tool) that could help explain this type of modeling to those of us who are naive in this area? "Zeroinflated count models" sounds very interesting. Cheers, Johnathon 
11032010, 11:58 PM  #7 
Senior Member
Location: Heidelberg, Germany Join Date: Feb 2010
Posts: 994

You don't need a "zeroinflated count model". The negative binomial (NB) model used, e.g., by edgeR, DESeq, BaySeq, is just fine. "Zeroinflated" means that you see zeroes more often than an NB (or Poisson) distribution would predict. This is an issue in some applications of statistics but I don't think we need it here.
For the question of reporting it: Well ,a fold change estimate of 300:0 is infinite, but if you don't like this, you can still say, it is ">300". In the end, what you want is some kind of confidence interval: If you have few counts, an observed count ratio of, say, 5:1, can be caused by a real expression ratio of anything from, say, 2:1 to 10:1, while, if the count ratio was 500:100, you can be way more sure that the estimate of 5:1 is close to the real value. Hence, a foldchange estimate should be supplemented by an indication of its precision (which, in RNASeq, would stringly depend on the number of reads the estimate is based on). There have been a few methods to get such confidence intervals for microarrays, but, to my knowledge, nothing is available yet for RNASeq data. Simon 
11042010, 12:53 AM  #8 
Member
Location: sweden Join Date: Aug 2010
Posts: 23

It would be nice to see some references addressing when to use zeroinflated NB and when ordinary NB.
UCLA Statistical Computing is a very nice stitistics source (http://www.ats.ucla.edu/stat/stata/dae/zinb.htm). "Some Strategies You Might Be Tempted To Try Before we show how you can analyze this with a zeroinflated negative binomial analysis, let's consider some other methods that you might use. OLS Regression  You could try to analyze these data using OLS regression. However, count data are highly nonnormal and are not well estimated by OLS regression. Zeroinflated Poisson Regression  Zeroinflated Poisson regression does better when the data is not overdispersed, i.e. variance much larger than the mean. Ordinary Count Models  Poisson or negative binomial models might be more appropriate if there are not excess zeros. " Hi Johnathon, unfortunately, I have not seen an R package for this. It may be "just fine" to use ordinary NB for some reviewers. I am definetely looking forward to seeing it be integreated in one of the open source packages. It would be even nicer to have the generlized estimated equation (GEE) or robust cluster variance function integreated with zeroinflated NB/Poisson models. 
11042010, 01:33 AM  #9  
Senior Member
Location: Heidelberg, Germany Join Date: Feb 2010
Posts: 994

Hi
Quote:
Just to make sure: If you take all the count values for the different genes in one sample and plot a histogram of these, there is no reason for it to look like NB, and many zeroes are totally expected. What I am talking about is the distribution of count values of one and the same gene in many different samples. This is what is expected to follow an NB. Unfortunately, this claim is made from theoretical reasoning, because no one has yet performed RNASeq experiments with hundreds of replicates, and with only a handful of numbers for each gene, you cannot estimate a distribution. Consequently, it would be hard to check whether there is an excess of zeroes, and hence, I suggest to assume there is not  not only to keep thing simple but also because it is hard to see a mechanistic cause for a zero excess. Simon 

11042010, 07:43 AM  #10 
Member
Location: Retirement  Not working with bioinformatics anymore. Join Date: Apr 2010
Posts: 63

What we've used for timecourse models is doing a pairwise comparison of each time point with the control (time 0) after quantile normalization. We've had a lot of success with doing a Fisher Exact Test for our pairwise comparisons (it's faster to do a chisquare approximation of the FET, but it isn't quite as good). For Illumina and SOLiD reads the hypergeometric distribution that the FET is based on describes the data really well after quantile normalization. Before normalization it doesn't work as well, and more length bias remains. In some samples we've seen a complete elimination of the length bias using this method, FWIW. The method has no problems determining if a zerocount sample is different from an ncount sample.
We also find this method needs some multiple hypothesis testing correction. We've had pretty good luck with the Bonferroni (Sidak) correction despite the fact it's really conservative. FDRs are commonly used in microarray data, and would probably work better than Bonferroni, but are a little more involved to calculate. This and other methods for pairwise comparison can be found in Bullard et al. (2010) in BMC Bioinformatics. I only tweaked one of their methods to get ours. 
11042010, 08:56 AM  #11  
Senior Member
Location: Heidelberg, Germany Join Date: Feb 2010
Posts: 994

Quote:
Simon 

11042010, 08:59 AM  #12 
Junior Member
Location: Canada Join Date: Jul 2010
Posts: 1

I am fascinated by the comments posted here. I would thank everybody who helps later when I have question about NGS, because you are the professionals.

11042010, 09:15 AM  #13 
Junior Member
Location: washington Join Date: Dec 2009
Posts: 7

Hi Simon Anders,
Thank you very much for the insight and for clarifying things. I think you've addressed all the follow up questions I had about the zeroinflated model. My inclination was to report these genes as >300 (0 reads vs 300 reads). I agree it might be better to indicate the precision as well. I appreciate all the feedback. 
11042010, 12:34 PM  #14 
Member
Location: sweden Join Date: Aug 2010
Posts: 23

It was microbiota composition data  histograms of each OTU (suppose this is the same as what you referred to as each gene) for groups of tens of samples. (OK, not RNAseq data) But I think for people doing 16S rRNA, an R package with zeroinflated model would be useful.
In cases of just a few replicates, certainly no way to check zero excess. But why use NB then (if this was what you meant to do)? Isn't there a simpler way? 
11042010, 07:20 PM  #15 
Senior Member
Location: WashU Join Date: Aug 2010
Posts: 117

Dealing with 0 counts in comparisons across conditions was a very common topic of debate in the good old days of SAGE analysis. Perhaps there are some insights in that literature. Thankfully I didn't have to analyze that data but perhaps a SAGE aficionado will jump in here...

11052010, 01:23 AM  #16  
Senior Member
Location: Heidelberg, Germany Join Date: Feb 2010
Posts: 994

Quote:
Simon 

Thread Tools  

