Library size-stabilized metacells construction enhances co-expression network analysis in single-cell data

TL;DR

LSMetacell improves co-expression network analysis in single-cell RNA data by stabilizing library sizes, reducing false correlations and revealing biologically meaningful interactions.

Contribution

LSMetacell introduces a novel computational framework that stabilizes library sizes during metacell construction to correct compositional biases in co-expression analysis.

Findings

LSMetacell reduces false-positive correlations caused by library size variance in single-cell RNA data.

Applied to Alzheimer’s disease data, LSMetacell identified microglia-specific co-expression modules linked to immune dysregulation and neurodegeneration.

LSMetacell outperforms conventional methods in preserving biological heterogeneity while mitigating technical noise.

Abstract

Single-cell RNA sequencing (scRNA-seq) deciphers cell type-specific co-expression networks to resolve biological functions but remains constrained by data sparsity and compositional biases. Conventional metacells construction strategies mitigate sparsity by aggregating transcriptionally similar cells but often neglect systematic biases introduced by compositional data. This problem leads to spurious co-expression correlations and obscuring biologically meaningful interactions. Through mathematical modeling and simulations, we demonstrate that uncontrolled library size variance in traditional metacells inflates false-positive correlations and distorts co-expression networks. Here, we present LSMetacell (Library Size-stabilized Metacells), a computational framework that explicitly stabilizes library sizes across metacells to reduce compositional noise while preserving cellular heterogeneity. LSMetacell addresses this by stabilizing library sizes during metacells aggregation, thereby enhancing the accuracy of downstream analyses such as Weighted Gene Co-expression Network Analysis (WGCNA). Applied to a postmortem Alzheimer’s disease brain scRNA-seq dataset, LSMetacell revealed robust, cell type-specific co-expression modules enriched for disease-relevant pathways, outperforming the conventional metacells approach. Our work establishes a principled strategy for resolving compositional biases in scRNA-seq data, advancing the reliability of co-expression network inference in studying complex biological systems. This framework provides a generalizable solution for improving transcriptional analyses in single-cell studies.

Gene co-expression analysis is a widely used method to infer functional relationships between genes by measuring correlations in their normalized gene expression level. However, in this paper, through mathematical modeling and simulations, we demonstrate that these correlations are systematically skewed—particularly due to biases caused by variability in sequencing depth (library size). This issue distorts co-expression analysis results, inflating false correlations and masking true biological interactions. Traditional methods fail to address library size biases in single-cell studies where data sparsity compounds these challenges. We introduce LSMetacell, a computational framework that simultaneously tackles single-cell data sparsity and corrects for library size-induced correlation biases. By constructing metacells with stabilized sequencing depths, our method reduces technical noise while preserving biological heterogeneity. Applied to Alzheimer’s disease brain data, LSMetacell uncovered microglia-specific co-expression networks linking immune dysregulation to neurodegeneration. Our work provides a dual solution: enhancing single-cell resolution through cell aggregation and mitigating systemic biases that plague co-expression studies. LSMetacell integrates technical approaches with biological analysis, enabling researchers to extract precise and reproducible findings from compositional data.

Figures

Relationship between correlation bias and library-size variation.a) Library-size variance-dependent correlation inflation in null data. Scatter plot demonstrates systematic overestimation of co-expression (Pearson correlation coefficients between gene pairs) (y-axis) with increasing mean expression levels (x-axis) under varying library size variance (Variance size factors). Variance size factors scale the variance of library size of the simulated data. The higher the variance size factor, the greater the library size variation in the simulated data. b) Library-size variance-dependent correlation absolute bias. Scatter plot showing the scale factor (variance of library size/mean of library size) and average absolute bias of the estimated correlation. The higher the variance size factor, the greater the library size variation in the simulated data.

Preservation of co-expression modules constructed by different methods across major cell types.Co-expression modules were constructed separately for excitatory neurons (EN), inhibitory neurons (IN), oligodendrocytes (OC), astrocytes (AC), and microglia (MC). The preservation of these modules was subsequently evaluated in a fully independent single-cell sequencing dataset using the Z-Summary statistic. A Z-Summary > 10 indicates strong evidence of preservation, while a Z-Summary > 2 suggests moderate preservation.

Microglia gene co-expression modules.Module names are prefixed with “MC” (abbreviation for microglia) followed by the color label. a) UMAP representation of the co-expression network visualizes individual genes (kMEs > 0.2) with eigengene of each module indicated by larger circles colored by the module color. Positive module eigengene correlations > 0.5 are indicated by grey lines connecting the eigengenes. b) Modules-phenotype association analysis. Correlation between the module eigengene and phenotypes. Numbers indicate correlation coefficients and p-values. Braaksc: Braak Stage is the semiquantitative measure (from no effect on neocortex to severity effect) of severity of neurofibrillary tangle (NFT) pathology. Ceradsc: CERAD score is the semiquantitative measure of neuritic plague (from severity to no effect). Cogdx: Clinical consensus diagnosis of cognitive status at time of death (from no cognitive impairment to dementia). Dcfdx_lv: Clinical diagnosis of cognitive status at last visit (from no cognitive impairment to dementia). c) The biological process gene ontology enrichment analysis of genes. GeneRatio is the ratio between genes of interest in the gene set and total genes of interest. Dot color represents the adjusted p-value (Benjamini–Hochberg method) of enrichment analysis.

Equations

Full text

Introduction

Gene co-expression network analysis serves as a powerful and efficient framework for deciphering transcriptional regulatory mechanisms across diverse biological contexts, offering an intuitive approach to unravel complex transcriptional interactions. [1,2]. While bulk RNA-seq enabled early discoveries, its fundamental limitation—confounding by cell-type heterogeneity—persists, masking cell-specific interactions. Single-cell RNA sequencing (scRNA-seq) holds promise to address this limitation by resolving transcriptional profiles at cellular resolution. However, the challenge of data sparsity, arising from the constraints of current single-cell sequencing workflows and technologies, continues to pose significant hurdles for downstream analysis. [3,4]. A conventional workflow accounts for these considerations by collapsing highly similar cells into “metacells” to reduce sparsity while retaining cellular heterogeneity [5,6]. However, the compositional nature of normalized scRNA-seq data (both before and after metacells aggregation) introduces critical yet underappreciated biases that undermine co-expression network inference [7].

In compositional data, where measurements represent relative abundances (e.g., gene counts normalized to library size, i.e., the total number of gene counts sequenced from a sample), spurious correlations arise due to the “closed-sum” constraint: an increase in one gene’s expression inherently distorts the apparent proportions of others [5,8]. While this issue plagues bulk RNA-seq analyses, it is exacerbated in scRNA-seq by two factors: (1) extreme library size variability across cells (e.g., often spanning 10-fold differences in unique molecular identifier [UMI] counts, where UMIs are short nucleotide tags used to accurately quantify transcript molecules and distinguish them from PCR duplicates), amplifying normalization artifacts; and (2) data sparsity, which compounds compositional noise by inflating stochastic zeros [9–11]. Conventional metacells construction strategies, designed primarily to mitigate sparsity, often neglect library size variance, generating metacells where normalized expression values (e.g., CPM/TPM, Log-CPM/TPM) remain compositionally biased, propagating technical noise into downstream network analyses [7,12,13].

Here, we rigorously demonstrate through mathematical modeling and in silico simulations that uncontrolled library size variance in metacells introduces systematic distortions in co-expression network inference. Specifically, it leads to false-positive correlations and obscures biologically meaningful interactions. Furthermore, to address this challenge, we introduce LSMetacell (Library Size-stabilized Metacells), a novel framework that explicitly stabilizes library sizes across metacells while preserving transcriptional heterogeneity. By explicitly accounting for library size variability, LSMetacell reduces compositionally induced biases, enabling more accurate and reliable co-expression network analyses, particularly for methods like WGCNA (Weighted Gene Co-expression Network Analysis) that rely on precise correlation estimates. Our approach not only enhances the robustness of co-expression network inference but also provides a generalizable solution for improving the accuracy of downstream transcriptional analyses in scRNA-seq studies. Finally, we evaluated the utility of LSMetacell by applying it to postmortem brain samples from Alzheimer’s disease patients and controls scRNA-seq data set and identified biological meaningful cell type-specific co-expression network.

Results

The magnitudes of correlation are overestimated by sequencing depth variations

To verify the impact of confounding effects of sequencing depth variations on the correlation of independent genes, we generated null datasets where genes are not co-expressed (null data) by permuting the single-nucleus RNA-seq (snRNA-seq) data from reference [14], while introducing different library size variation across cells (as detailed in the Methods section). Fig 1a shows that in null data with no library size variations, there were minimal biases. However, as the variations increased, both the mean and variance of correlation estimation increased. Especially, when the variance approached that of the original data (Variance Size Factor = 1), the biases were almost comparable to those observed with a tenfold increase in the original data. Furthermore, we prove mathematically that the growing variations of library size induce higher type I error, tending to over-reject independence between gene co-expression (S1 File). Our theoretical derivation aligns with the results from simulation.

Relationship between correlation bias and library-size variation.a) Library-size variance-dependent correlation inflation in null data. Scatter plot demonstrates systematic overestimation of co-expression (Pearson correlation coefficients between gene pairs) (y-axis) with increasing mean expression levels (x-axis) under varying library size variance (Variance size factors). Variance size factors scale the variance of library size of the simulated data. The higher the variance size factor, the greater the library size variation in the simulated data. b) Library-size variance-dependent correlation absolute bias. Scatter plot showing the scale factor (variance of library size/mean of library size) and average absolute bias of the estimated correlation. The higher the variance size factor, the greater the library size variation in the simulated data.

We further create a simulated snRNA-seq dataset, incorporating predefined potential correlations that are known in advance (as detailed in the Methods section). We observed that as the library size variability increased, the deviations in the estimated correlations also grew larger (Figs 1b and S1). Specifically, when the variance in library sizes was high, the estimated correlations exhibited significant biases, deviating substantially from the true predefined correlations. This finding underscores the critical impact of library size variability on the accuracy of correlation estimates in snRNA-seq data, highlighting the need for robust normalization and variance stabilization techniques to mitigate such biases.

The library stabilized metacell method achieved robust co-expression modules

To tackle the confounding effects of library size variability in single-cell co-expression analysis, we developed the Library size-Stabilized Metacell (LSMetacell) algorithm. This method is based on the hypothesis that cells of the same type share similar gene expression patterns. The LSMetacell method simultaneously accomplishes two objectives: (1) aggregating transcriptionally similar cells into metacells (2) stabilizing library size variation across metacells. This dual approach, which aims to preserve biological signals and minimize technical variance, addresses the key challenge in single-cell co-expression analysis, where library size heterogeneity often distorts downstream network inference.

The LSMetacell algorithm processes three key inputs: a cell-cell similarity network computed from transcriptional profiles (Pearson correlation among cells by default), a gene expression count matrix, and the target number of metacells. It then iteratively builds metacells by: (1) initiating with seed cells and incorporating their most similar neighbors; (2) dynamically evaluating the cumulative library size relative to the global mean metacells library size $[eqn]$ ; and (3) excluding cell additions that would substantially alter the library size distribution. This dual refinement strategy guarantees transcriptional homogeneity while maintaining uniform library sizes across all metacells (See Methods for implementation details). When evaluated on null datasets (under the condition variance = 1 in Fig 1a), LSMetacell introduced minimal technical bias compared to other benchmark methods, demonstrating its enhanced accuracy in mitigating artifactual correlations (S2 Fig).

To rigorously evaluate the enhancement of LSMetacell on WGCNA for single-cell data, we carried out a comprehensive comparison with the conventional high-dimensional WGCNA (hdWGCNA) method, which leverages k-nearest neighbor aggregation without library size stabilization [7]. In addition to this primary comparison, we also incorporated four other cutting-edge algorithms capable of generating metacells, even though they were not specifically designed for gene co-expression analysis. The first is metacells2 [15,16], which adopts a divide-and-conquer strategy to partition large single-cell datasets and identifies metacells by constructing k-nearest neighbor graphs within subsets. Another is metaQ [17], which leverages a deep learning framework to quantize single cells into a limited codebook. We also included SEACells, a metacell identification method that uses a kernel-based approach to learn a low-dimensional embedding of single-cell data and then employs sparse combinatorial optimization to group cells into archetypal aggregates [18]. Additionally, SuperCell is a method that constructs metacells by merging the most transcriptionally similar cells with a k-nearest neighbor network [19]. Furthermore, we incorporated the primary version of LSMetacell into the comparison, where cells are aggregated as LSMetacell without strict library size control. Using human dorsolateral prefrontal cortex scRNA-seq data spanning five major cell types (excitatory neurons [EN], inhibitory neurons [IN], astrocytes [AC], oligodendrocytes [OC], and microglia [MC]), we evaluated two key aspects of network quality: biological relevance and technical robustness. Notably, LSMetacell achieved the lowest coefficient of variation, indicating superior library size stability of the generated meta-cells, which is a property essential for reducing noise in downstream analyses (S2 Table). For biological validation, we leveraged the principle that functionally related genes (as evidenced by protein-protein interactions) should show stronger co-expression. Analysis of the top 10,000 correlated gene pairs revealed superior Protein-Protein Interaction (PPI) enrichment with LSMetacell across most cell types (Table 1).

Table 1: PPI enrichment among the top 10,000 correlated gene pairs across cell types. Values denote the number of gene pairs exhibiting significant PPI enrichment for each method within each cell type. Paired Wilcoxon test comparing LSMetacell with all benchmark methods across cell types yielded P = 0.031 for every comparison.

Technical robustness was evaluated through module preservation analysis in an independent human prefrontal cortex scRNA-seq data (Fig 2). LSMetacell demonstrated higher and stable module preservation, indicating greater reproducibility of the identified co-expression patterns. The improved PPI enrichment and module preservation underscore its utility for identifying conserved regulatory programs. These findings demonstrate that stabilizing library size variation in metacells enhances co-expression analysis. It is worth noting that the method of aggregation used in the primary version of LSMetacell has the potential to generate metacells with reduced library size variation. On this basis, LSMetacell further decouples technical artifacts from biological variation through library size variation stabilization, thereby enhancing co-expression analysis.

Preservation of co-expression modules constructed by different methods across major cell types.Co-expression modules were constructed separately for excitatory neurons (EN), inhibitory neurons (IN), oligodendrocytes (OC), astrocytes (AC), and microglia (MC). The preservation of these modules was subsequently evaluated in a fully independent single-cell sequencing dataset using the Z-Summary statistic. A Z-Summary > 10 indicates strong evidence of preservation, while a Z-Summary > 2 suggests moderate preservation.

The microglia immune response in AD

Using LSMetacell method followed by weighted gene co-expression network analysis (WGCNA) on a comprehensive single-cell transcriptomic dataset, we identified 73 modules across five major central nervous system cell types (EN: 18; IN: 16; AC: 14; MC: 12; OC: 13). Strikingly, only four microglial modules (MC-Brown, MC-Green, MC-Purple, MC-Yellow) exhibited significant enrichment for immune-related pathways (GO analysis, FDR < 0.01), whereas neuronal modules predominantly mapped to synaptic function, metabolic processes, etc. This cell-type-specific functional segregation validates our algorithm’s precision in recovering biologically coherent signals. This highlights the unique position of microglia as the core immune effector of central nervous system in AD pathology (Fig 3a).

Microglia gene co-expression modules.Module names are prefixed with “MC” (abbreviation for microglia) followed by the color label. a) UMAP representation of the co-expression network visualizes individual genes (kMEs > 0.2) with eigengene of each module indicated by larger circles colored by the module color. Positive module eigengene correlations > 0.5 are indicated by grey lines connecting the eigengenes. b) Modules-phenotype association analysis. Correlation between the module eigengene and phenotypes. Numbers indicate correlation coefficients and p-values. Braaksc: Braak Stage is the semiquantitative measure (from no effect on neocortex to severity effect) of severity of neurofibrillary tangle (NFT) pathology. Ceradsc: CERAD score is the semiquantitative measure of neuritic plague (from severity to no effect). Cogdx: Clinical consensus diagnosis of cognitive status at time of death (from no cognitive impairment to dementia). Dcfdx_lv: Clinical diagnosis of cognitive status at last visit (from no cognitive impairment to dementia). c) The biological process gene ontology enrichment analysis of genes. GeneRatio is the ratio between genes of interest in the gene set and total genes of interest. Dot color represents the adjusted p-value (Benjamini–Hochberg method) of enrichment analysis.

Module MC-Brown exhibited strong positive correlations with AD severity markers (Fig 3b). This gene module orchestrates immune cell signaling and regulation through GTPase activity (e.g., RAC1, ARHGAPs, DOCK proteins), driving processes like leukocyte activation, phagocytosis, and chemotaxis [20,21]. Furthermore, lipid metabolism perturbations (dysregulated glycerophospholipid and phosphatidylinositol biosynthesis) and vesicle trafficking (autophagy, endocytosis) amplify microglial membrane dysfunction, impairing Aβ clearance and amplifying inflammatory cascades [22–24] (Fig 3c).

The MC-Green module is more significantly enriched in pathways associated with adaptive immunity compared to the MC-Brown module (Fig 3c). It is worth noting the GO results of T cells, especially CD8-positive T cells. Previous studies have shown the double-edged role of T cells in the mouse model of AD [25]. Su et al. reported the protective role of CD8 ⁺ T cells in AD development [26], while Chen et al. also reported that microglia-mediated infiltration of T cells exacerbates neurodegeneration in the tauopathy-associated model [27]. Our results show that this module is negatively correlated with AD indicators. This may indicate the protective role of this module in AD. The observed decline in MC-Green activity may reflect exhaustion of adaptive-immune surveillance as the disease progresses, thereby accelerating the emergence of AD-related signatures [28]. In contrast, this module profiled with other benchmark approaches showed a notably closer association with innate-immune rather than adaptive-immune processes, along with a weaker correlation with AD indicators (S3–S8 Figs).

MC-Yellow shows significant pathway in chaperone-mediated protein folding, ATP synthesis coupled electron transport, and immune related terms (Fig 3c). MC-Yellow hub genes contain ribosome-related and some known AD-related genes, including APOE [29], HSP90B1 [30], TREM2 [31], etc. The MC-Red module occupies a central position in the network analysis, with its functional enrichment linked to RNA metabolism (Fig 3c). Notably, the hub genes within this module are primarily associated with mitochondrial pathways (MT-ND1, MT-ND2, MT-ND3, MT-ND4, MT-ND4L, MT-CO1, MT-CO2, MT-ATP6, MTCO3), a finding with particular pathological relevance given the established role of microglial mitochondrial dysfunction in driving Alzheimer’s pathogenesis through mechanisms involving neuroinflammation and metabolic failure. [32]. However, when applying certain benchmark methods, many genes now assigned to the MC-Red module were merged into the MC-Brown or MC-Yellow modules, highlighting the enhanced resolution offered by the present workflow (S3–S8 Figs).

Discussion

Gene co-expression network analysis is indispensable for uncovering transcriptional regulatory mechanisms in diverse biological contexts. Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool to overcome the limitations of bulk RNA-seq in dealing with cell-type heterogeneity. However, the sparsity of single-cell data and the compositional nature of normalized data pose significant challenges for accurate co-expression network inference.

Our study rigorously demonstrates the detrimental impact of uncontrolled library size variation in metacells on co-expression network analysis. Through mathematical modeling and in silico perturbations, we found that this variance introduces false-positive correlations and obscures biologically meaningful interactions. The biases in correlation estimates increase with higher variability in sequencing depths, as evidenced by both the analysis of null datasets and simulated snRNA-seq data. These findings highlight the urgent need for methods that can mitigate the confounding effects of library size variability in single-cell co-expression analysis.

A standard approach to address differences in library size is by scaling the counts for each cell using a specific factor, e.g., its total UMI count. While effective in controlling for technical variability, this method introduces a critical artifact: applying different scaling factors across cells artificially inflates gene-gene correlations, as the same multiplicative factor is applied to all genes within a cell.

To address this issue, we introduced the Library Size-Stabilized Metacell (LSMetacell) algorithm. By aggregating transcriptionally similar cells into metacells while simultaneously stabilizing their library sizes, LSMetacell effectively reduces technical noise and preserves biological signals. Our proposed pepline controls for library size effects through a two-stage process. In the first stage, we form metacells by pooling UMI counts from clusters of similar cells. This initial aggregation step naturally equilibrates library sizes across metacells to a large extent, but more importantly, it maintains gene-specific variance. For instance, the counts for each gene in a metacell represent the sum from its constituent cells, thereby preserving co-variation patterns. In the second stage, we perform a library size normalization on the count profiles of these metacells. The key innovation lies in applying normalization after aggregation. Since the normalization factors are computed across metacells, which already exhibit a more uniform library size distribution, the variance among these factors is significantly reduced. As a result, scaling each metacell to a common total count does not materially distort the gene-gene correlation matrix. This two-step procedure thus effectively controls for any residual differences in metacell size while minimizing the introduction of spurious correlations, leading to a more reliable representation of the data for downstream analysis.

The significant improvement in protein-protein interaction (PPI) enrichment and module preservation in an independent scRNA-seq dataset of human prefrontal cortex demonstrates the biological accuracy and technical robustness of LSMetacell. This method provides a reliable framework for single-cell co-expression network analysis, enabling the decoupling of technical artifacts from biological variation.

Applying the LSMetacell method to a comprehensive single-cell transcriptomic dataset in Alzheimer’s disease (AD), we identified distinct co-expression modules in major central nervous system cell types, with a particular focus on microglia. The functional segregation of these modules, such as immune-related pathways in microglial modules and synaptic/metabolic functions in neuronal modules, validates the precision of our algorithm in recovering biologically coherent signals. These findings confirm the unique role of microglia as the core immune effector in AD pathology. Given the critical role of immune mechanisms in AD, we focused on the co-expression network modules of microglia.

Among the microglial modules, MC-Brown is strongly correlated with AD and is enriched in immune activation pathways. The regulation of small GTPase-mediated signal transduction, innate immunity, and lipid metabolism perturbations in this module provide potential therapeutic targets. For example, targeting small GTPases could potentially disrupt microglial chemotaxis towards amyloid-β (Aβ) plaques. Similarly, restoring normal lipid metabolism in microglia may enhance Aβ clearance and alleviate the inflammatory cascades. The MC-Green module, which is more significantly enriched in adaptive immunity pathways and negatively correlated with AD indicators, suggests a potential protective role. However, its changing correlation patterns in control and AD groups imply possible exhaustion of the adaptive immune system in AD. Future studies could explore ways to boost the function of this module to harness its protective effects against AD progression. MC-Yellow and MC-Red also present valuable insights into AD pathogenesis. The involvement of chaperone-mediated protein folding, ATP synthesis coupled electron transport, and RNA metabolism in these modules provides new perspectives on the complex molecular mechanisms underlying AD. For instance, enhancing chaperone-mediated protein folding may help prevent the accumulation of misfolded proteins and inhibit pro-inflammatory process which are characteristic features of AD.

In the broader context, our work not only provides a novel methodological framework for single-cell co-expression network analysis but also offers in-depth insights into the molecular mechanisms of AD. The LSMetacell algorithm has the potential to be applied in other biological systems and diseases, facilitating the discovery of conserved regulatory programs and potential therapeutic targets. Future research could focus on validating the proposed targets in in vivo models and exploring the translational potential of these findings in clinical settings.

However, it is important to acknowledge that, like all count-based sequencing data, our meta-cell expression matrices remain inherently compositional. Although meta-cell aggregation significantly reduces technical variability and mitigates spurious correlations caused by variable library sizes across single cells, it does not eliminate the compositional nature of the data, as the data remain subject to the unit-sum constraint. This means that the elevated expression of one gene necessarily leads to the relative decrease of others, potentially introducing negative biases in correlation estimates. Another notable limitation of LSMetacell is its dependence on pre-annotated and transcriptionally homogeneous cell populations. While our method effectively controls for technical variation in library size during metacell construction, it is restrictive that cell type labels are uncertain, incomplete, or derived from continuous differentiation trajectories. This design reflects a deliberate trade-off in Algorithm 1: by focusing on library size stabilization within annotated groups, we prioritize suppressing technical variation for co-expression analysis at the potential expense of broader metacell versatility. Interestingly, our benchmarking results revealed that metacell2, a method that also achieves relatively stable library sizes, performed well in module preservation and other metrics, as might be expected. Furthermore, the co-expression modules derived by metacell2 exhibited greater similarity to those from LSMetacell than those from other methods. Surprisingly, SuperCell, despite showing relatively higher library size variation, also achieved competitive module preservation results, and its co-expression patterns were likewise more similar to those of LSMetacell compared to others. These observations suggest that, in addition to library size control, other factors such as low compactness and high separation of metacells may also play important roles in obtaining biologically meaningful co-expression networks. Future work could explore complementary approaches to extend the utility of library size stabilization to more complex and heterogeneous single-cell datasets, while incorporating features that enhance metacell compactness and transcriptional coherence. Furthermore, while we have significantly reduced the number of parameters compared to previous metacells algorithms (only the number of metacells needs to be provided), it is still crucial to select a cell-cell similarity algorithm. Fortunately, we achieved good results by simply using Pearson correlation in our study, but this aspect warrants further investigation. Finally, the time complexity of LSMetacell scales quadratically ( $[eqn]$ ) with the number of cells ( $[eqn]$ ), making computational efficiency an important consideration for large datasets (S3 Table). Constructing metacells separately for each individual and cell type, rather than pooling cells across sources, enhances biological specificity while simultaneously reducing computational expense.

Conclusions

Through mathematical modeling and in silico simulations, we rigorously demonstrate that excessive library size variance in single-cell data systematically inflates gene-gene correlations, fundamentally compromising co-expression network fidelity. We introduce LSMetacell, a rational metacell framework that strategically balances transcriptional similarity with library size stabilization. By leveraging scRNA-seq’s cellular granularity to aggregate phenotypically coherent cells while controlling technical variability, LSMetacell mitigates false-positive interactions. Applied to Alzheimer’s disease data, our method identified robust cell type-specific modules enriched for neurodegeneration pathways, including microglia-driven immune dysregulation. This work establishes that systematic technical stabilization during metacell design is essential for reliable single-cell network inference, advancing mechanistic discovery in complex biological systems.

Methods

Data collection and preprocessing

Single-cell RNA sequencing (scRNA-seq) datasets were obtained from public repositories: dorsolateral prefrontal cortex (DLPFC) data from Synapse (syn18485175; [14]); prefrontal cortex data from Synapse (syn21261143; [33]). Gene count matrices for single-cell data were generated using the original authors’ preprocessing pipelines [14,33]. Then, we performed log-normalization with a scale factor of 10,000.

Simulated dataset

To generate null data sets from an authentic scRNA-seq dataset with co-expression levels at or close to zero among all gene pairs while pre-serving gene expression levels, we adopt the following approach that combines permutation with Poisson sampling. First, we normalize the expression level of each gene $[eqn]$ in cell $[eqn]$ , denoted as $[eqn]$ , where $[eqn]$ represents the raw expression count and $[eqn]$ is the scaling factor (e.g., sequencing depth) for cell $[eqn]$ . Subsequently, for each gene $[eqn]$ , we randomly shuffle its normalized expression values $[eqn]$ across all $[eqn]$ cells. This shuffling step decorrelates the gene expressions, effectively eliminating any potential co-expression patterns among gene pairs. Finally, to generate the UMI (Unique Molecular Identifier) counts for the permuted data, we sample from a Poisson distribution with parameters derived from the normalized, permuted expression levels. Specifically, for gene $[eqn]$ in cell $[eqn]$ , the UMI count is sampled from $[eqn]$ , where $[eqn]$ is the desired library size in cell $[eqn]$ . The sequencing depth $[eqn]$ is determined from a truncated normal distribution, with $[eqn]$ almost surely, using variance of the original scaling factors $[eqn]$ multiplied by different variance size factor. This process ensures that the null data sets maintain the overall gene expression characteristics of the original data while lacking any meaningful co-expression patterns.

To generate synthetic single-cell RNA-seq data with controlled co-expression structures, we developed a simulation framework integrating a Gamma-Poisson marginal model with a Gaussian copula. First, the Gamma distribution parameters, $[eqn]$ and $[eqn]$ for the gene $[eqn]$ , were estimated from authentic scRNA-seq data to ensure realistic overdispersion. Next, given a predefined co-expression matrix $[eqn]$ , calibrated from the real data, we sampled latent variables $[eqn]$ from multivariate normal distribution, $[eqn]$ . These latent variables were transformed by

$$[eqn]$$

where $[eqn]$ is the cumulative distribution function (CDF) of a standard normal distribution and $[eqn]$ is the CDF of $[eqn]$ . This approach simulates gene expression level $[eqn]$ , maintaining their predefined correlation structure for benchmark. Finally, the simulated UMI count, $[eqn]$ , is sampled from $[eqn]$ , where $[eqn]$ scales expression to match desired library sizes. By generating $[eqn]$ from a truncated normal distribution (larger than 0) with different variance we simulated datasets with low, moderate, and high library size variability.

LSMetacell method

The LSMetacell algorithm is specifically designed to operate within pre-defined, transcriptionally homogeneous cell populations. Before applying the method, users must first annotate cell types using established clustering approaches to ensure biologically meaningful aggregation. At its core, the algorithm constructs metacells through an iterative process that balances transcriptional similarity (default: Pearson correlation between cells) with library size stabilization. The algorithm begins by calculating the target average library size $[eqn]$ across all cells. Starting with the smallest-library cell as a seed, it iteratively aggregates cells based on transcriptional similarity. At each step, a candidate cell is probabilistically selected from unassigned cells based on its aggregate affinity to the current metacell. Crucially, the algorithm enforces library size stabilization: a candidate is retained only if adding it brings the metacell’s cumulative library size closer to $[eqn]$ . This process continues until the metacell’s library size approaches $[eqn]$ or no suitable candidates remain. The workflow proceeds as Algorithm 1.

Algorithm 1 LSMetacell

Input: $[eqn]$ , set of single cells with raw gene expression profiles; $[eqn]$ , target number of metacells; $[eqn]$ , cell-cell similarity matrix for cells in $[eqn]$ .

Output: $[eqn]$ , a list of metacells with stabilized library sizes.

Initialization:

Target average library size: $[eqn]$ ;

Unassigned cells: $[eqn]$ ;

Metacell list: $[eqn]$

Main Loop:

while $[eqn]$ do

Step1: Seed Selection

Select the cell with the smallest library size:

$[eqn]$ ;

Initialize current metacell:

$[eqn]$ ;

$[eqn]$ ;

$[eqn]$ ;

Step2: Iterative Aggregation

While $[eqn]$ and $[eqn]$ do

Compute affinity scores between $[eqn]$ and all $[eqn]$ :

$[eqn]$ , $[eqn]$ ;

Normalized weights:

$[eqn]$ , $[eqn]$ ;

$[eqn]$ with probability $[eqn]$ ;

Tentatively update:

$[eqn]$

if $[eqn]$ then

Accept candidate:

$[eqn]$ ;

$[eqn]$ ;

$[eqn]$ ;

else

Reject $[eqn]$ .

end if

end while

end while

To extend the LSMetacell algorithm for datasets with heterogeneous groups (e.g., individuals, biological replicates), we developed a group-aware metacells construction framework (Algorithm 2).

Algorithm 2 LSMetacell by Group

Input: $[eqn]$ , set of single cells in group $[eqn]$ ; $[eqn]$ , cell-cell similarity matrix for cells in group $[eqn]$ ; $[eqn]$ , the target number of metacells for group $[eqn]$ .

Output: $[eqn]$ , a list of metacells with stabilized library sizes.

Main Loop:

Step1: Proportional Metacells Allocation

For each group $[eqn]$ , calculate its target metacell count:

$[eqn]$ ;

Step2: Group-Specific Metacell Construction

For each group $[eqn]$ :

Apply Algorithm 1 LSMetacell.

WGCNA and identification of significant modules

Gene co-expression network construction.

Following metacells construction, weighted co-expression network were constructed by the common pipeline. First, the gene-gene similarity matrix $[eqn]$ is computed by taking the pairwise signed correlation of genes.

$$[eqn]$$

To comply with scale-free topology criterion and the recommendations of WGCNA use, we chose appropriate soft-thresholding powers to convert the gene expression matrices to adjacency matrices. Then topology overlap matrices (TOM) were calculated by adjacency matrices [34]. We then use hierarchical clustering and dynamic tree cut method to identify gene clusters [35].

Protein-Protein Interaction (PPI) enrichment.

The protein-protein interaction annotation were download from STRING database (12.0) [36].

Module preservation.

We used Z-Summary, a network preservation statistic aggregating multiple preservation statistics (3 density-based statistics and 3 connectivity-based statistics), to quantify the conservation of the co-expression network in another dataset [37].

Module eigengenes.

Given that each module comprises genes with correlated expression patterns, it is logical to summarize each module using a single representative expression profile, known as the module eigengene. The module eigengene is defined as the first principal component of the standardized expression matrix of the genes within the module, which captures the majority of the variance within the module. Module eigengenes lead to a natural measure of similarity (membership) of all individual genes to all modules [38]. A continuous measure of module membership of gene $[eqn]$ in module $[eqn]$ is defined as

$$[eqn]$$

where $[eqn]$ is the expression measurement of gene $[eqn]$ and $[eqn]$ is the eigengene of module $[eqn]$ [39].

Gene enrichment analysis.

We implemented gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis by ClusterProfiler (4.15) [40].

Comparison of LSMetacell with other metacell algorithms.

To compare LSMetacell with other metacell construction algorithms, including hdWGCNA, metacell2, metaQ, SEACells and SuperCell. we generated an identical number of metacells across the same cell types using each of these methods, all under their respective default parameter settings. For primary LSMetacell (LSMetacell without strict library size control), we omitted lines (25–29) in LSMetacell (Algorithm 1) and ceased the metacell aggregation process once the required cell number (number of cells/number of designed metacells) was achieved. To ensure a fair and biologically relevant comparison of co-expression networks across different metacell methods, we applied a unified gene filtering criterion prior to WGCNA analysis: only genes expressed in at least 10% of cells within each cell type were retained. This filtering step was consistently applied to metacells generated by all methods, including LSMetacell, hdWGCNA, metacell2, metaQ, SEACells, and SuperCell. Although the absolute number of retained genes varied slightly across methods due to differences in metacell construction, the overlap among gene sets was substantial. To enable a direct and methodologically consistent comparison of functional modules, all subsequent WGCNA analyses were performed using the intersection of these filtered gene sets across methods.

Supporting information

S1 FileProof for the main results.Contains the theoretical proof of our main results, where we establish theoretically that compositional gene expression data suffers from inflated Type I errors in correlation analyses.(PDF)

S1 TableVariability in library size distribution of metacells generated by different methods across major cell types.Values represent the coefficient of variation (standard deviation divided by the mean), indicating the degree of variability in library sizes within each metacell type.(PDF)

S2 TableResults of one-sided Wilcoxon rank-sum tests (left > right) comparing Zsummary scores between pairs of metacell construction methods presented in Fig 2.Each cell displays the p-value from the statistical test evaluating whether the method in the row exhibits significantly higher Zsummary scores than the method in the column. Empty cells indicate comparisons that were either not applicable or not performed.(PDF)

S3 TableThe table shows the LSMetacell algorithm’s runtime on a standard laptop (64 GB RAM, Windows 11 OS) using datasets of varying sizes.The results demonstrate a near-quadratic scaling of runtime with cell count, which aligns with its theoretical O( $[eqn]$ ) complexity.(PDF)

S1 FigDistribution of correlation absolute bias under a specific scaling condition.This histogram illustrates the distribution of correlation absolute bias values across all gene pairs for different scale factor. Each scale factor is corresponding to one representative scaling condition shown in Fig 1b. The x-axis represents the correlation absolute bias, while the y-axis shows the frequency of gene pairs at each bias interval. Lower and more concentrated distributions indicate better control of technical artifacts induced by library size variation.(PDF)

S2 FigBenchmarking metacell algorithms under a synthetic null-correlation dataset.We generated a synthetic single-cell dataset in which genes possess no intrinsic pairwise correlation (see Fig 1 Methods) and applied seven metacell-building algorithms: hdWGCNA, Metacell2, SEACells, SuperCell, MetaQ, Primary, and LSMetacell. After normalization within each metacell set, all pairwise Pearson correlations between genes were recalculated; their frequency distributions are displayed above. To quantify differences in noise suppression, two-sided Wilcoxon signed-rank tests were conducted between every distribution and the LSMetacell-derived distribution (our reference). P-values are reported directly on the plots.(PDF)

S3 FigMicroglia gene co-expression modules identified by hdWGCNA employing the same analysis and presentation panels as in Fig 3 (see this figure for a detailed caption).To ensure comparability across methods, the modules identified by each algorithm were aligned with those from LSMetacell using Fisher’s exact test to assess significant overlaps. Thus, modules labeled as MCgreen, MCred, MCyellow, MCblue, MCturquoise, MCblack, MCbrown, MCpink, MCpurple, MCyellow, and MCmagenta in each method correspond significantly to their counterparts in LSMetacell. Any additional modules identified by a method that are not listed above do not exhibit consistent correspondence across the other algorithms.(PDF)

S4 FigMicroglia gene co-expression modules identified by Metacell2 employing the same analysis and presentation panels as in Fig 3 (see this figure for a detailed caption).To ensure comparability across methods, the modules identified by each algorithm were aligned with those from LSMetacell using Fisher’s exact test to assess significant overlaps. Thus, modules labeled as MCgreen, MCred, MCyellow, MCblue, MCturquoise, MCblack, MCbrown, MCpink, MCpurple, MCyellow, and MCmagenta in each method correspond significantly to their counterparts in LSMetacell. Any additional modules identified by a method that are not listed above do not exhibit consistent correspondence across the other algorithms.(PDF)

S5 FigMicroglia gene co-expression modules identified by MetaQ employing the same analysis and presentation panels as in Fig 3 (see this figure for a detailed caption).To ensure comparability across methods, the modules identified by each algorithm were aligned with those from LSMetacell using Fisher’s exact test to assess significant overlaps. Thus, modules labeled as MCgreen, MCred, MCyellow, MCblue, MCturquoise, MCblack, MCbrown, MCpink, MCpurple, MCyellow, and MCmagenta in each method correspond significantly to their counterparts in LSMetacell. Any additional modules identified by a method that are not listed above do not exhibit consistent correspondence across the other algorithms.(PDF)

S6 FigMicroglia gene co-expression modules identified by SEACells employing the same analysis and presentation panels as in Fig 3 (see this figure for a detailed caption).To ensure comparability across methods, the modules identified by each algorithm were aligned with those from LSMetacell using Fisher’s exact test to assess significant overlaps. Thus, modules labeled as MCgreen, MCred, MCyellow, MCblue, MCturquoise, MCblack, MCbrown, MCpink, MCpurple, MCyellow, and MCmagenta in each method correspond significantly to their counterparts in LSMetacell. Any additional modules identified by a method that are not listed above do not exhibit consistent correspondence across the other algorithms.(PDF)

S7 FigMicroglia gene co-expression modules identified by SuperCell employing the same analysis and presentation panels as in Fig 3 (see this figure for a detailed caption).To ensure comparability across methods, the modules identified by each algorithm were aligned with those from LSMetacell using Fisher’s exact test to assess significant overlaps. Thus, modules labeled as MCgreen, MCred, MCyellow, MCblue, MCturquoise, MCblack, MCbrown, MCpink, MCpurple, MCyellow, and MCmagenta in each method correspond significantly to their counterparts in LSMetacell. Any additional modules identified by a method that are not listed above do not exhibit consistent correspondence across the other algorithms.(PDF)

S8 FigMicroglia gene co-expression modules identified by Primary employing the same analysis and presentation panels as in Fig 3 (see this figure for a detailed caption).To ensure comparability across methods, the modules identified by each algorithm were aligned with those from LSMetacell using Fisher’s exact test to assess significant overlaps. Thus, modules labeled as MCgreen, MCred, MCyellow, MCblue, MCturquoise, MCblack, MCbrown, MCpink, MCpurple, MCyellow, and MCmagenta in each method correspond significantly to their counterparts in LSMetacell. Any additional modules identified by a method that are not listed above do not exhibit consistent correspondence across the other algorithms.(PDF)