Running differential expression analyses on count matrices with many zeroes

Question

Note: I also posted this issue (with less context) in the bioconductor support site: https://support.bioconductor.org/p/97424/

I'm working on a snakemake workflow that identifies various small RNA species in C. elegans small RNA-seq libraries. Some are endogenous siRNAs supposedly generated from RNA templates (through RNA-dependant RNA polymerases (RdRP)) that can be variously classified (protein-coding genes, transposons, and other repeat types).

I count such small RNA reads in a set of libraries and try to identify differentially producing sources. For this, I use DESeq2 (that I run using rpy2 from within snakemake).

I'm not sure DESeq2 is always appropriate for this kind of data, but so far, the analyses would at least complete. However, I recently added new potential types of small RNA source (simple repeats and satellites), and these happen to have low counts. I'm not 100% sure, but I suspect these low counts are the reason for failures during DESeq2 analyses (for debugging purposes, I ran this manually in R):

> dds <- DESeq(dds, betaPrior=T)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  : 
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
2: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
3: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
4: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
5: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
6: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
7: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
8: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
9: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
10: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
11: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
12: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight

The count matrix indeed contains a high proportion of zeroes, the rest are mostly ones, but there are some sources that seem a little more significantly producing small RNAs:

> mean(counts_data == 0)
[1] 0.7488889
> mean(counts_data == 1)
[1] 0.1507407
> max(counts_data)
[1] 34

How would you recommend handling such kind of data?

Should I for instance discard rows without enough counts? If so, what would be a reasonable threshold?

Edits: trying to get a dispersion plot:

The same error occurs when trying to estimate dispersion:

> dds <- estimateSizeFactors(dds)
> dds <- estimateDispersions(dds)
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  : 
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)

This seems to prevent the generation of dispersion plots:

> plotDispEsts(dds)
Error in min(py[py > 0], na.rm = TRUE) : 
  invalid 'type' (list) of argument
In addition: Warning message:
In structure(x, class = unique(c("AsIs", oldClass(x)))) :
  Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

Using fitType="local" also fails:

> dds <- DESeq(dds, betaPrior=T, fitType="local")
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  : 
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)

But fitType="mean" works:

> dds <- DESeq(dds, betaPrior=T, fitType="mean")
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

The dispersion plot then looks as follows:

What do these fitType options mean, and are there good reasons why they would have to be changed in order to accomodate for low counts?

Answer 1

I don't think that the issue is the low counts, but rather the number of features without any real variance (the black dots at the bottom).

So what the heck is the dispersion plot and why does one need to fit it anyway?

In a typical RNAseq experiment, one measures many thousands of genes with only a few replicates per biological group. This then leads to issues when performing testing, since your variance measurement for any given gene will be larger or smaller than reality, simply due to the lack of a large number of replicates. Consequently, at least since the limma package came out, most RNAseq packages have tried to fit the mean-dispersion with some sort of line and then used that as a background distribution toward which the variance assigned to each gene should be "shrunken" (i.e., moved toward). Getting this at least approximately correct is important, since if you have few replicates and the final "shrunken variance" is significantly wrong, then your p-values will also be wrong.

What are the various fitType options and what are they doing?

The three options that one can use to fit this relationship are parametric, local, and mean. parametric is the default and basically results in fitting to a bent line. This usually works quite well, but you can see that you have two clouds of data, one shaped like a bent line (at the top) and another (the black dots at the bottom) with no variance. That's causing this fit type to fail. Likewise, local is doing a local fit and failing for the same reason. mean is just taking the mean (average) of all of the estimates. This is fool-proof, but you're then regressing to something that really doesn't fit your data. mean is really a last choice if nothing else works.

For your particular dataset, my main concerns would be figuring out (1) why the maximum mean counts are around 20 (that's really low) and (2) why so many of these features have incredibly low variance.

Answer 2

Remove low-count features in advance. This is standard for most tools including DESeq2 and edgeR (see section 2.6).

This will keep you from testing a lot of features that cannot be differentially expressed (using NB model you need at least four reads - below that the uncertainty is too high to ever be DE) and you will make the false discovery correction unnecessarily hard.