scRNA-Seq pseudobulk for marker detection

Question

A standard approach for scRNA-Seq is to partition the single cells into individual clusters, then use a Wilcoxon test to find markers that characterize each cluster (or other statistical methods that consider single cells as replicates).

Recently, there has been a push towards using pseudobulk approaches for "differential state" (DS) analysis, which avoid the pseudoreplication bias of considering single cells as replicates. See e.g. Crowell (2020), Zimmerman (2021), Thurman (2021), Squair (2021), or the OSCA book.

However, all the above papers seem to use pseudobulk methods for DS of a single cluster between two conditions (e.g. drug treatment), or occasionally to compare a specific pair of clusters. I haven't seen any example of using pseudobulk to find markers distinguishing clusters (comparisons of one cluster against all others).

I can think of a few reasons, for example:

• these pseudobulk methods can be more computationally intensive,
• perhaps the Wilcoxon test is good enough to simply find markers,
• pseudobulk methods can only work when each cluster is represented in several samples, which is hard to guarantee for every single cluster in the dataset.

On the other hand it seems to me a pseudobulk approach should yield more relevant markers (genes that are robustly enriched across biological replicates). So, am I missing something? Is there a statistical reason not to use pseudobulk approaches for finding cluster-specific markers?

Answer 1

The only reason I see for not using pseudobulks is lack of biological replication. Even though single-cell sequencing is common these days and standardized it is still expensive. If you have several conditions plus biological replicates and need many cells, maybe due to heterogeneity, the bill can easily reach the 5-digit realm of money. Biological replication is always desired but sometimes simply not possible. We recently did a study in which we were interested in organ heterogeneity in different conditions and both mouse work and library prep was so extensive that even hash tagging the mice we had to pool was simply not possible. So we were bound to using each cell as a replicate. I prefer limma-trend for this, not the Wilcox test. It is not computationally very demanding, neither is pseudobulks. But even if it was, I would happily wait 30min rather than 10sec to get the best possible rather than the fastest results. Pseudobulks conceptually have the advantage that you basically are forced to have biological replication (because only then they make sense) so findings are intrinsically more robust. Theoretically, you can of course have biological replicates and still treat each cell as replicate, and then blocking or adding the biological replication information into the DE testing design, but the statistical rigor of methods used for pseudobulks (DESeq2, edgeR and other) is for me the reason to always use pseudobulks if I can.

I personally use pseudobulks also for markers. After all marker detection comes downstream of DE testing by aggregating the results into lists that contain unique or enriched genes per cluster. I see no reason to use multiple testing strategies.

Answer 2

If you want a statistical reason not to use pseudobulk to do DE for marker finding, I think it is the same reason not to use any DE method for marker finding: the double dipping problem (see Gao and Witten 2020).

Briefly, if you cluster a dataset based on genes (we can safely say that clustering in PCA space is still based on genes), then asking which genes are differentially expressed between clusters is a statistical fallacy. We could just as easily divide apples in 2 groups according to their shade of red and then ask which RGB channels exhibit a statistically significant difference.

This means that even in the absence of real structure, you will identify very low p-values for your DE test. The OSCA book section you mention links to another section in the advanced book that tells you the same thing with a small practical demonstration generating a 200-dimensional Gaussian and clustering it, then testing for differences in clusters.

This does not mean, however, that markers found this way are meaningless; it simply means that you cannot use a p-value cutoff to define whether a gene is a marker in this scenario, or say that since you have statistically significant markers two clusters are different.

You can still rank these markers by fold change, difference in percentage expression, AUC, entropy, etc as well as p-value. Moreover, markers identified using the traditional methods tend to make biological sense, and statistically rigorous methods are still in development (there was a beautiful talk by Anna Neufeld on an approach involving count splitting, you can read the preprint here).

The test you suggest is interesting, but I don't think you would solve the double dipping problem; this being said it looks more robust and appropriate than using every single cell as a replicate.