Can DESeq and edgeR deal with in-balanced RNA-seq data?

[asiangg]

We did three biological replicates for our treatment and control using RNA-seq 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?

[Simon Anders]

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 log-log 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 batch-condition 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

[markrobinsonca]

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, M-vs-A 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

[Simon Anders]

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