Single-cell RNA-seq: deep & few or shallow & many?

My lab does a lot of single-cell RNA-seq. Our goals are to find out what cell types are present, characterize their gene activity, and (if possible) understand how the different subpopulations are related to each other (for example through pseudotemporal analysis).

We often wonder which is better: should we sequence many cells at a shallow level, or few cells deeply? There's a paper on this here, but it is completely empirical. I am posting this question because I want to develop intuition via statistical analysis: given a reasonable stochastic model, what would we expect from shallow versus deep sequencing?

To help formulate a stochastic model, suppose a typical mammalian cell has 200,000 RNA molecules, and we detect each molecule with probability 1% (for low coverage) or 2.5% (for medium coverage).

I posted an answer below. I would be interested in answers that improve on mine in the following ways.

• Account for correlations between difference genes, or generally follow up on EdM's suggestion to account for natural variability in mRNA levels.
• Correct my math, especially the downright incorrect assumption that $\sum_jz_{ij}$ is constant.
• Quantify tradeoffs for more interesting target parameters. For instance, if you want to know whether a given cell state transition is continuous or discrete, you might want to test the adequacy of a two-cluster model (say, with K-means) versus a continuous model (say, 1-dimensional PCA). This ought to make the shallow sequencing worse, because higher technical noise will blur the clusters together, even though their centers will still be accurately estimated.

I am not so interested in incorporating domain knowledge about single cell sequencing. For instance, I am aware that deeper sequencing leads to redundancy as the same barcode is detected twice, but that is beyond the scope of my question and it belongs in another discussion.

• Another set of considerations you should add to your question, to get a more useful answer: what is the distribution of expression among different RNA species, and how much do you care about detecting RNA species with low levels of expression?
– EdM
Dec 20, 2017 at 20:08
• The title should probably read "deep & few or shallow & many"? Dec 20, 2017 at 20:22
• @amoeba: indeed. Fixed. Dec 20, 2017 at 20:56
• @EdM: unfortunately, I don't know much about the distribution of mRNA counts, though I'd like to learn more. I care about low-expressed genes, such as transcription factors, enough to want decent cluster-level averages for them. Dec 20, 2017 at 21:06
• Not worth a full answer, but Gaussian mixture models (GMM's) have sample complexities that depend inversely on some power of the total-variation distance between the distributions. This total variation distance decreases with shallow sequencing, because it depends a lot on the variance of a single cell. For GMM's, the sample complexity should go up with shallower sequencing. See this: people.csail.mit.edu/moitra/docs/disentangling.pdf Dec 20, 2017 at 21:13

Here's an answer, but it makes naive assumptions about the data generating process, it is mathematically incorrect, and it only works with a very simple estimand (group average expression). Throughout this answer, $i$ indexes a cell and $j$ a gene.

Model

Suppose the cells come from a homogeneous population, with mRNA counts varying due to natural stochasticity. (Link is to a tutorial that derives a stochastic model for gene expression from first principles). Suppose the molecule counts for cell $i$, call them $x_i$, have approximately mean $\vec\mu$, variance $\delta \vec\mu$, and zero covariance (each like an overdispersed Poisson). We can sample $y_{ij} \approx Poisson(0.01*x_{ij})$ for $i \in 1 ... 10k$ via shallow sequencing or $z_{ij} \approx Poisson(0.025*x_{ij})$ for $i \in 1 ... 4k$.

Variance decomposition

We know that, for the measurement of mRNA species $j$,

\begin{align}Var(y_{ij}) =& Var[E[y_{ij}|x_{ij}]] + E[Var[y_{ij}|x_{ij}]] \\ =& Var[0.01x_{ij}] + E[0.01 x_{ij}] \\ =& 0.01^2\delta\mu_{j} + 0.01\mu_j\\ \end{align}.

For $z$, the analogous result is $0.025^2\delta\mu_{j} + 0.025\mu_j$. The first summand can be regarded as biological variance and the second as technical.

Estimation

When we go to estimate $\mu$, it's typically done using "normalized" data: we do the bulk of the exploratory analysis on $\tilde{z_j} = 10,000*z_{i}/\sum_jz_{ij}$ or $\tilde{y_i} = 10,000*y_{i}/\sum_jy_{ij}$. This is not optimal, because deeper-sequenced cells should be given more weight, but I am answering within this framework because this is how most people get things done with single cell data.

Derivation of consequences

Suppose we ignore the variability in $\sum_j y_{ij}$ and $\sum_jz_{ij}$. For all $i$, $\sum_jz_{ij} = 5000$ and $\sum_jy_{ij} = 2000$. (If you can do better than this, I'm interested! Please write your comment or answer below.) To simply estimate mean expression, the technical variances will be equal across strategies. If $\bar z\equiv \sum_i\tilde z_i / 4000$ and $\bar y\equiv \sum_i\tilde y_i / 10000$:

Z:

• $Var[\bar z_j] = (1/4000)*Var[\tilde z_{ij}]$
• $Var[\tilde z_{ij}] = Var[2*z_{ij}] = 4*Var[z_{ij}]$
• $Var[z_{ij}] = 0.025\mu_{j}$ (technical component only)
• $Var[\bar z_j] = (1/4000)*0.025*4\mu_j = 2.5\times 10^{-5}\mu_j$ (technical component only)

Y:

• $Var[\bar y_j] = (1/10000)*Var[\tilde y_{ij}]$
• $Var[\tilde y_{ij}] = Var[5*z_{ij}] = 25*Var[y_{ij}]$
• $Var[y_{ij}] = 0.01\mu_{j}$ (technical component only)
• $Var[\bar y_j] = (1/10000)*0.01*25\mu_j = 2.5\times 10^{-5}\mu_j$ (technical component only)

The biological component of the variability is lower in the shallow sample: the multiplier over the technical variance is $0.01\delta$, rather than $0.025\delta$.

Summary

This answer indicates that technical variation is equal across samples, but the biological component of the variability is lower in the shallow sample.

• Are you sure the link is correct? The text seems to be about "master equation" in chemistry and does not even contain the word "gene". Dec 20, 2017 at 21:48
• Yes, that is what I meant to post. It's very low-level and probably not useful, except that it gives a mechanistic justification for certain stochastic models of biochemical systems. Dec 21, 2017 at 14:04