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.
 A: 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.
