Does MetaCell2 perform normalization of raw UMI counts?

[benchsar]

Hello,

Thank you for developing MetaCell2 it's a very promising tool that I'm currently exploring for CD34+ single-cell datasets.

While reviewing the documentation and pipeline, I couldn’t find a clear explanation regarding normalization of raw UMI counts (e.g., CPM, log-normalization, etc.) prior to or during metacell construction.

Could you please clarify whether MetaCell2 performs any implicit normalization, for example through probabilistic downsampling in compute_direct_metacells, to correct for differences in sequencing depth between cells?

If so, I would appreciate it if you could confirm:

At which step this normalization is applied,

And whether the values in adata.layers["downsampled"] can be considered normalized and suitable for comparing expression levels across cells or metacells.

Thank you very much for your help, and again for sharing this great tool!

Sara

[orenbenkiki]

We do downsample cells during metacells computation, due to the wild difference in depth of sampling of cells within a single dataset. This is however not a good way to get estimates of gene expression levels in the results. Downsampling is done only for the selected ("feature") genes and is done separately in each pile (when the data is large enough that divide-and-conquer is needed). It is discarded once we have built the pile's KNN graph (the paper goes into more details). Yes, we also do a log internally on the downsampled UMIs when we compute the metacells (we correlate the logs and not the fractions, to compute distances to build the KNN graphs), but again, that's not something that you can or should use for examining the resulting metacells.

Instead, we provide the total UMIs for each metacell and for each gene in each metacell, which you can use to compute the linear fraction of each gene out of the total. Such an estimate is wildly unreliable when applied to single cells, but the whole point of metacells is that they are just large enough to make this robust (and not too large so you can still see details of cell behavior within the "same" cell type).

We compute a fraction per metacell (by default using geomean, but there's a flag to force this back to linear fractions - see the documentation of the collect_metacells function). Geomean is less sensitive to one or two cells with very high expression dominating the results. This looks better on gene-gene plots, but because we renormalize these fractions to sum to 1, it is is essentially impossible to compare these normalized geomean fractions with UMI counts in single cells - for that you'd be better off with using the linear fractions. TBH we are re-evaluating the geomean approach...

Final note - because we (like everyone else doing scRNA-seq) are forced to use relative fractions of UMIs out of the total, the choice of "total" is crucial. In particular you can't meaningfully compare fractions between two data sets unless you ensure you use the same set of genes as the denominator in both (typically, choosing the common set of genes). This of course requires recomputing the fractions... For example, we do such a renormalization as one of the 1st steps in the projection algorithm (to make the atlas and query fractions comparable).

As for log - yes, most of our graphs use log2(fraction + 1e-5), we find these to be convenient units. It also allows us to consider "fold factor" (difference of such logs) to express distance between expression (fold factor 1 is 2x, 2 is 4x, etc.). However that's entirely up to you.

Final side note - while the vignette recommends excluding the mitochondrial genes from analysis (and therefore the denominator), we also discovered it is beneficial to also exclude the ribosomal genes, as in some (most?) data sets they vary wildly between cell (types) - having only 5% ribosomal genes in one cell and 35% in another means that all the other genes appear to differ by ~30% just because the denominator is different. We've been working on an updated vignette but are holding off on publishing it because we want to stabilize the new methodology. In the meanwhile, just adding "RP[SL]." to the "MT-." pattern of excluded genes is a quick and dirty solution.

Hope this helps - you can also see a list of all the properties with an explanation of each in the main package README.

[benchsar]

Hi,

Thank you for your detailed explanation — it helped me better understand how expression fractions are derived in MetaCell2.

After computing metacells using collect_metacells(), I performed UMAP and then clustered the metacells (e.g., using Leiden) based on the expression profiles. I used linear gene fractions (derived from the total_umis matrix normalized per metacell) to compare expression between clusters.

Just to clarify two things:

Is it reasonable to cluster metacells themselves after collect_metacells() to further explore subpopulations?

And if so, can we use the normalized gene fractions for basic comparisons (e.g., fold change or t-test) between those metacell clusters, within the same dataset?

Thanks again for your time and for developing this great tool.

Best,
Sara

[orenbenkiki]

Sure, it is standard practice to group metacells together and assign each such group a type name. When we do this we don't give up on the individual metacells though; the fact that they have the same "type" doesn't mean they don't vary (e.g. gradients). Also sometimes what is a "type" is subjective - depending on the study, some type labels may be broad and some may be more specific. We also don't typically use Leiden to define the groups - that is, we do use it as an initial estiomte, but we quickly switch to looking at scatter plots showing the (log) fraction of all the metacells with one gene in the X and another gene in the Y axis - what we call "gene-gene plots". We use these and look at genes of interest to more finely define the metacell groups (and ensure that they are "real" and that we don't have lateral genes messing up the picture). We look at genes based on both prior biological knowledge, and also by looking at "marker" genes (genes which have variance between metacells and high expression in at least some metacells). This isn't a trivial process and often we discover more genes we needed to mark as lateral, so we have to re-compute the metacells - the vignette goes into detail with examples of such a process.

Yes, we use the normalized (log) fractions of the genes in this process. It is tempting to use t-test to estimate whether a difference between metacells is "real", but, here be dragons. The actual distribution of the genes is very much not normal (or Poisson, or log-normal, or negative binomial, or ...). That is, some genes are well behaved, but most aren't (in particular we see way more zero UMIs than you'd expect for reasonable distribution models). So you have to be very careful applying such statistics. We use various techniques and at some points even resort to saying "unless there are at least a total of 40 UMIs in both compared metacells for the same gene, don't believe the fold factor for that gene". In general we don't care about the fold factor between weakly-expressed genes so it works out, more or less....

Another approach is to group genes into modules (this can be specific to a type), and sum the UMIs of all the genes in a module, to get an even more robust signal of some biological behavior. We are working on ways to automate detection of these gene modules, right now it is essentially a manual process. A good tool for this is showing a heatmap of gene expression across the metacells (this also works when looking at the marker genes to get an idea of which genes to use to define the types, and which types exist in the data), and computing correlations between genes. If done "well" you get gene modules that make sense biologically and that provide very robust scores for biological behaviors of interest (differentiation etc.). These modules are often specific to some type(s) and don't play a significant role in other type(s).

Proper analysis of the data is, alas, as much an art as it is a science at this point. It is an active area of research for us to try and automate these steps in a way that preserves biological signals of interest and ignores the very large noise, especially from lateral genes (including genes that weren't marked as lateral, but actually are, or are highly correlated to lateral genes). We are making some progress but right now this is a manual process that depends on the skill and experience of the analyst. I wish I had a better answer. Ask me again next year :-)

[orenbenkiki]

Final thought - we have the MCView tool, which generates the graphs and assists in the task of type annotations, in the browser. It is written in R - see https://github.com/tanaylab/MCView - you can ask Aviezer for help if you have trouble running it.