
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Can I use DESeq and edgeR for mixed ANOVA  abebe  Illumina/Solexa  18  07132017 04:45 AM 
Can edgeR/DESeq have more than one covariate?  arrchi  Bioinformatics  8  10282013 02:37 PM 
DESeq and edgeR up/down regulation  murphycj  Bioinformatics  7  09212011 07:09 AM 
BaySeq vs GLM in EdgeR and DESeq  Hilary April Smith  Bioinformatics  0  08032011 10:14 AM 
edgeR vs DESeq vs bayseq  Azazel  Bioinformatics  1  10072010 07:11 AM 

Thread Tools 
05212010, 06:39 AM  #1 
Member
Location: New York Join Date: Dec 2008
Posts: 44

Can DESeq and edgeR deal with inbalanced RNAseq data?
We did three biological replicates for our treatment and control using RNAseq to find out which transcripts have differential expression. To make sure we are obtaining genuine changes, we did another batch of experiments several months later. Now, we have:
Batch1: 3 treatments vs. 3 controls Batch2: 3 treatments vs. 3 controls The two batches were done under the same conditions(hopefully). However, there is a significant difference in total read count. The first batch contains ~10 million reads for each replicate but the second batch contains ~30 million reads for each. It is because Illumina has improved chemicals and software. I applied several tools (including DESeq, edgeR and limma) to identify differential genes from the two batches of data. The 1st batch yields ~500 genes and the 2nd batch yields ~200 genes. To our disappoint, the two lists contain very small overlaps. We suspect one set of treatments or controls was screwed so decided to switch the treatment and control of the two batches to identify the bad ones. To our surprise, the two batches yield 10 fold more genes after switching! That means, each batch now contains ~5000 differential genes and they overlap by 70%!! This cannot be biologically true and I suspect it is related with the unbalanced inputs of treatment vs. control. To my knowledge, both DESeq and edgeR try to normalize the library sizes internally before performing statistical tests. However, the question is how well is that done? Any input or suggestions? 
05212010, 08:34 AM  #2 
Senior Member
Location: Heidelberg, Germany Join Date: Feb 2010
Posts: 994

I'm not quite sure what you mean by switching. Are you now comparing treatment from batch 1 with control from batch 2?
But two answer your question: Both DESeq and edgeR adjust for library size. While edgeR uses the library sizes that you tell it, DESeq tries to estimate them from the data. To see whether this worked well, I'd suggest that you choose pairs of samples and divide all the counts from one sample by the size factor for this sample (for DESeq; for edgeR, take the total read count) and do likewise for the other. Then plot one against the other in a loglog scatter plot and mark the diagonal (with abline(a=0,b=1) ). Check that the points scatter symmetrically around the diagonal. Do this for a couple of sample pairs. In my experience, however, the library size normalisation works well and is unlikely to be the culprit. A good idea might be to check sample distances: With DESeq, make a CountDataSet containing all 12 of your samples. The perform a variance stabilizing transformation, get a distance matrix for the variance transformed matrix and plot it as a heatmap. I have described this procedure in the DESeq vignette. If all is well, the replicates should cluster together. If a sample does not cluster with its replicates, you might want to exclude it from the analysis. Lastly, have a look at the scvPlots in your four batchcondition combinations. What is the raw SCV value in the region of highest count density, i.e., at the peak of the black density curve? Is it maybe much larger in some cases than in others? Cheers Simon 
05232010, 04:01 PM  #3 
Junior Member
Location: Zurich Join Date: Mar 2010
Posts: 7

A couple comments from the edgeR camp ...
I agree with Simon that just a pairs() plots of read counts is a useful initial diagnostic, especially if you think you might have sample switching (I didn't fully understand what was switched from your description). Also, MvsA plots (edgeR does 'smear' plots) would be quite useful. One clarification of what Simon said with respect to edgeR. While its true that edgeR uses the library sizes "that you tell it", there is a function in there for calculating normalization factors from the data  calcNormFactors()  and a description in the manual of how to build that into your library sizes. I haven't compared directly, but its roughly similar to the DEseq calculation for this. The normalization (which is beyond just accounting for library size) is described at: http://genomebiology.com/2010/11/3/R25/abstract Another alternative to explore sample relations is the plotMDS.dge() function in edgeR. This is essentially a principal components plot, but specific to count data. Hope that helps. Cheers, Mark 
05262010, 04:45 AM  #4 
Senior Member
Location: Heidelberg, Germany Join Date: Feb 2010
Posts: 994

Hi Mark
It seems I haven looked into the edgeR vignette for a while and missed that you added a size estimation by now. You are right, your and our scheme are very much the same. We all had the same idea of looking at the quotient between individual gene counts and taking some robust location estimator of their distribution. The only difference is that you used a trimmed mean and we went all the way to maximal trimming, i.e., used the median. This definitely shouldn't make much of a difference. Simon 
Thread Tools  

