Let's imagine a single cell experiment in which we have 3 biological replicates, treated (TR) and untreated (UNT).
After all the necessary filtering and integration steps, we isolate a cluster of interest (cluster X), for which we want to test differential gene expression (DE) between TR and UNT.
Ideally one would filter lowly expressed genes, aggregate counts within each replicate x protocol to generate pseudobulks, and use tried and tested NB methods such as the ones implemented in DESeq2 or edgeR to get DE between conditions.
At the same time, we have also tested for differential abundance (DA) between TR and UNT (using a statistical framework such as Milo), and we know that TR is more abundant than UNT in X - say 60% TR and 40% UNT, for the sake of argument.
My question is: how do you account for DA when doing DE? Is it even necessary?
I had thought of 2 ways but I am not convinced by any of them:
- downsampling (by random sampling without replacement) each replicate so the number of cells matches the smallest number, in this case so that TR and UNT have the same number of cells. One could also use more refined sampling strategies such as geometrical sketching, k-means + KNN etc;
- adding DA as a numeric covariate in the design.
I don't like 1. because I would be literally throwing counts away (especially if DA is strong) and I would introduce an otherwise unnecessary random element, and I don't like 2. either, because it looks like an overly strong correction.
An argument against accounting for DA is that the inclusion of depth size factors in the model would take care of any systematic, global difference given by compositional effects. However, I worry that genes that are not consistently quantified across cells would be artificially inflated by this unequal coverage, i.e. by virtue of having many 0's they would have a higher chance of being DE after aggregation, just because they were supported by more cells. I don't know if the count depth normalization would be enough to overcome this issue.
Moreover, even if we didn't care about single genes and focused more global results such as GSEA, there is no guarantee that the same pathways will be significantly enriched before and after downsampling, or that the enrichment score will follow the same direction, because rankings may change if we downsample cells and "lose" counts for several genes. Honestly, I don't even know how to go about benchmarking this particular aspect - I believe it is difficult to simulate pathway enrichments with currently available tools. I tested GSEA on DE results before and after downsampling on a dataset and it looks like some NES values for significant (in at least one of the two enrichment tests) pathways have an opposite sign.
I am familiar with the OSCA book chapter on handling DE and DA, but I don't feel like it addresses this particular question as it just suggests to run both analyses (e.g. here).
Do you have an informed opinion, or any suggestion?