Airqtl dissects cell state-specific causal gene regulatory networks with efficient single-cell eQTL mapping

Domain: neuroscience
Published: 2025-12-10
Editorial outcome: minor_revision

Abstract
Single-cell expression quantitative trait loci (sceQTL) mapping offers a powerful approach for understanding gene regulation and its heterogeneity across cell types and states. It has profound applications in genetics and genomics, particularly causal gene regulatory network (cGRN) inference to unravel the molecular circuits governing cell identity and function. However, computational scalability remains a critical bottleneck for sceQTL mapping, prohibiting thorough benchmarking and optimization of statistical accuracy. We present airqtl, an efficient method to overcome these challenges through algorithmic advances and efficient implementations of linear mixed models. Airqtl achieves superior time complexity and over 10^8 times of acceleration, enabling objective method benchmarking and optimization. Airqtl offers de novo inference of robust, experimentally validated cell state-specific cGRNs that reflect perturbation outcomes. Our results dissect the drivers of cGRN heterogeneity and underscore the value of natural genetic variations in primary human cell types for biologically relevant single-cell cGRN inference.

Introduction
Single-cell technology offers a unique opportunity to infer gene regulatory networks (GRNs) at an unprecedented scale. Recent computational advances1–4 have enabled the inference of cell type-specific and dynamic GRNs from a single sample. These approaches significantly reduce experimental costs and effort while enhancing the specificity of GRNs for distinct cell types, developmental stages, and extracellular conditions (collectively referred to as “states”). However, many existing methods focus on optimizing predictive performance between genes, often confusing causation with reverse causation or confounding on its own. These noncausal GRNs often fail to predict perturbation outcomes, limiting their reliability in applications such as disease understanding and drug discovery. While some approaches introduced causal evidence with transcription factors (TFs) and their DNA-binding events, these partial GRNs overlook other important regulatory mechanisms, such as signaling pathways, chromatin modifications, and non-coding RNAs.

Causal inference, particularly with perturbations, provides a robust framework for distinguishing causation and accounting for diverse regulatory mechanisms. This paradigm has been successful across various natural and social sciences, including single-cell experiments such as Perturb-seq5–8. However, these studies primarily involve engineered cell lines subjected to strong artificial perturbations. Since GRNs are highly dynamic and cell state-specific1–4, the relevance of these findings to primary cell types in natural human populations with both strong and mild perturbations remains suboptimal. In contrast, Mendelian randomization (MR) leverages natural genetic perturbations as instrumental variables to infer causality between organismal and molecular traits9–13. MR has been effectively employed to infer tissue-specific causal GRNs (cGRNs) from population-scale bulk RNA sequencing (RNA-seq), using expression quantitative trait loci (eQTLs) that directly regulate cis-gene expression and indirectly influence trans-genes14.

Applying MR to population-scale single-cell RNA-seq (scRNA-seq) promises greater cell state specificity and cost-effectiveness than bulk RNA-seq for cGRN inference. While dozens of population-scale scRNA-seq studies have been conducted (e.g.,15–19), the challenges in single-cell eQTL (sceQTL) mapping15,20 have limited its utility for downstream cGRN inference. On the statistical front, single-cell data sparsity necessitates specialized normalization techniques5,21,22, while the shared genome of cells from the same donor requires (generalized) linear mixed models16,23. However, the primary bottleneck lies in computational scalability. Without scalable methods, it is impractical to benchmark and improve statistical performance at a reasonable cost and time. Consequently, most population-scale scRNA-seq studies rely on ad hoc sceQTL mapping computer programs, with few capable of mapping all cis- and trans-sceQTLs across all single nucleotide polymorphism (SNP)-gene pairs. Some studies employ pseudo-bulk strategies for eQTL mapping, but they face challenges with continuous cell states, technical confounders, and limited statistical power15,20,24. Accurate and efficient eQTL mapping is also pivotal in human genetics25 and genomic foundation models26,27. The absence of a scalable sceQTL mapping method has hindered progress in these domains.

In this paper, we present airqtl, an efficient method to overcome these challenges in sceQTL mapping and infer cell state-specific cGRNs. Airqtl introduces a scalable exact algorithm that accelerates basic matrix operations by leveraging the donor-cell hierarchy in population-scale single-cell datasets. This advancement improves the time complexity of linear mixed models beyond state-of-the-art methods such as FaST-LMM28 and GEMMA29. Together with GPU utilization and other implementation optimizations, airqtl achieves over eight orders of magnitude in computational acceleration at a comparable or higher statistical accuracy. This scalability enables objective benchmarking and refinement of cell state-specific sceQTL mapping and admixed population modeling. Finally, airqtl’s efficient and accurate sceQTL mapping allows the de novo inference of robust, experimentally validated cell state-specific cGRNs. Comparative analyses reveal intrinsic and extrinsic drivers of cGRN heterogeneity, underlining the importance of cell state-specific experimental data and computational methods.

Results
Airqtl overview
The core of airqtl is an Array of Interleaved Repeats (AIR), a compressed data structure designed to efficiently handle repeated values, a pattern common in sceQTL mapping (Fig. 1a). SceQTL mapping requires genotype matrix expansion from donor- to cell-level, leading to repeated values that incur substantial computational and storage overhead. AIR circumvents this inefficiency by storing only their original values and repeat patterns, enabling scalable exact algorithms to accelerate basic matrix operations essential for sceQTL mapping, such as multiplication and singular value decomposition (Supplementary Note).

For instance, AIR factorizes cell-level genetic relationship matrix (GRM) in O(n_donor^3) time complexity, compared to O(n_cell * n_donor^2) required by state-of-the-art methods like FaST-LMM28 and GEMMA29. This principle underpins similar accelerations across all linear mixed model (LMM) steps in sceQTL mapping, yielding a total time complexity of O(n_SNP * n_gene * n_donor + (n_donor + n_cell) * (n_donor + n_cell) * n_donor) (Methods, Supplementary Note). When compared against the current state-of-the-art time complexity28,29 of O((n_SNP + n_gene) * (n_cell + n_donor) * n_cell), AIR’s algorithm alone would deliver a ~1000-fold speedup for sceQTL mapping (Methods).

Airqtl performs sceQTL mapping for each cell state separately using a streamlined two-step process (Fig. 1b, Methods). First, quality-controlled scRNA-seq read counts and covariates are normalized for each cell state using Normalisr5. Normalisr is optimized for scRNA-seq normalization and differential expression analysis with linear models, delivering both high accuracy and exceptional speed compared to generalized linear model-based methods. Normalisr’s accuracy and efficiency has been uniquely demonstrated in mapping the cis- and trans-effects of CRISPRi gRNAs and inferring cGRNs from large-scale Perturb-seq studies with high multiplicity of infection, a scenario that closely resembles sceQTL mapping.

Next, airqtl applies LMMs to identify sceQTLs from the normalized data (Fig. 1c, Methods). LMMs are particularly suited for scRNA-seq data as it accounts for the genetic similarity among cells from the same donor using the genetic relationship matrix (GRM). This approach is analogous to its use in genome-wide association studies (GWAS) for handling population structures, such as twins, to improve statistical power and reduce false positives30. Besides the algorithmic accelerations provided by AIR (Fig. 1a), airqtl implements additional optimizations such as GPU utilization and out-of-core computing. Airqtl also efficiently maps cell type/state-specific sceQTLs with an extended LMM including all linear and quadratic interaction terms that receives similar accelerations (Methods).

Efficient and accurate sceQTL mapping with airqtl
We conducted comprehensive benchmarking to evaluate the performance improvements offered by airqtl. In GWAS and eQTL mapping, simulated datasets have been pivotal for method benchmarking due to the lack of ground-truth information28,31–33. Following this precedent, we simulated population-scale scRNA-seq datasets of various sizes, incorporating ground-truth characteristics closely resembling real-world studies (Methods). These simulations were based on Randolph et al.19, a published scRNA-seq dataset comprising influenza A virus (flu) or mock (NI)-exposed human peripheral blood mononuclear cells (PBMCs) from male donors of African and European descent. This dataset was subsequently used for sceQTL mapping and cGRN inference in this study.

We benchmarked airqtl against CellRegMap23, the only existing method explicitly designed for sceQTL mapping with frequentist inference. We excluded other methods such as GASPACHO34 and SAIGE-QTL24 for three reasons: i) their Bayesian or heuristic frameworks reduce compatibility with downstream analyses or population-scale scRNA-seq datasets (Methods); ii) they lack formal statistical tests for sceQTL context specificity; iii) we encountered critical challenges while trying to run these methods on official or custom datasets.

Due to CellRegMap’s inability to scale to full cis- and trans-sceQTL mapping in large population-scale scRNA-seq datasets, we first tested both methods on small-scale datasets (Fig. 2a), where airqtl already demonstrated substantial advantage. Because running time often scales nonlinearly with problem dimensions, we employed a log-log model to extrapolate CellRegMap’s running time for large-scale datasets and compared it with the actual runtime of airqtl (Methods). The results revealed that airqtl achieved over eight orders of magnitude of acceleration and six orders of magnitude of cost reduction (Fig. 2b).

Accurate effect size estimation is crucial for eQTL mapping and downstream analyses, such as cGRN inference. However, it is often overlooked in single-cell differential expression studies and is particularly challenging due to expression-dependent biases5. To assess this issue in sceQTL mapping, we followed a bias-variance decomposition approach5, stratifying errors by expression level (Methods). Airqtl greatly mitigated the expression-dependent underestimation bias of sceQTL effect sizes observed in CellRegMap (Fig. 2c–f). Additional analyses attributed this improvement to the use of Normalisr for scRNA-seq normalization (Supplementary Fig. 1, and Supplementary Note), which has been shown to reduce such biases in differential expression analyses5. This advantage is expected and confirmed in Supplementary Note to be especially critical for cell type/state-specific sceQTL mapping, where expression-dependent biases can erroneously attribute changes in expression levels as cell type/state-specific sceQTLs.

Besides its advancements in efficiency and effect size estimation, airqtl demonstrated comparable or superior sensitivity and specificity (Supplementary Fig. 2, and Supplementary Data 1, Methods). Taken together, airqtl provides a robust and comprehensive solution for efficient and accurate sceQTL mapping.

Efficient and optimized cell type-specific sceQTL mapping
The interaction between genotype and other factors is well established in genetics35,36. Numerous bulk and single-cell eQTL studies have revealed tissue or cell type-specific eQTL effects15–19,37. While sceQTLs have been mapped separately for each cell type, statistical tests are essential to differentiate true cell type specificity from randomness in data but remain rare and underexplored. Airqtl overcomes these limitations by enabling efficient cell type-specific sceQTL mapping with an extended linear mixed model that incorporates linear and quadratic interactions between genotype, cell type, and other covariates (Methods). To systematically evaluate performance, we simulated population-scale scRNA-seq datasets with cell type-specific eQTLs based on CD4+ and CD8+ T cell subpopulations from the Randolph et al. dataset19 (Methods).

In terms of computational efficiency, airqtl was over five orders of magnitude faster and three orders of magnitude cheaper than CellRegMap (Fig. 3a, b, Methods). This scalability enables genome-wide mapping of cell type-specific sceQTLs. Although CellRegMap’s limited scalability prevented statistical benchmarking on large datasets, airqtl showed robust performance in estimating cell type-specific effect sizes (Fig. 3c, d, and Supplementary Fig. 3, Supplementary Fig. 4, Methods). Interestingly, cell type-specific effects were easier to estimate than sceQTL effects (Fig. 2ef), displaying lower bias and variance.

False positives in cell type specificity, in the form of biased null P value distributions, can lead to unnecessary and costly follow-up experiments. Initial results from airqtl showed an overabundance of small null P values for both non-eQTLs and non-specific eQTLs (Kolmogorov-Smirnov/KS P < 10^-99 and P < 10^-50 respectively, Fig. 3e, and Supplementary Data 2, Methods). Given the weak but nonzero expression-dependent biases in sceQTL effect size estimation (Fig. 2e), this result was anticipated and motivated an optimization of the model. To address this, we calibrated P values using Beta distributions following existing bulk eQTL studies37,38. To overcome expression-dependent biases in scRNA-seq, we further incorporated stratification with additional parameters. We evaluated 10 candidate parameters, including known single-cell differential expression confounders, their derived parameters, and key population genetics parameters (Methods). Each parameter was used individually to parameterize the Beta distribution, which then models the null P values obtained from separately simulated datasets (Supplementary Fig. 5). The learned parameterization was then used to calibrate P values in the original dataset. Calibration using the parameter “n0min” (proportion of zero-read cells for one gene, minimum across cell types compared) yielded the best null P value performance, with KS P = 0.009 and P = 0.12 for non-eQTLs and non-specific eQTLs, respectively (Supplementary Data 2, Supplementary Fig. 6). This optimization substantially improved the uniformity of null P value distributions (Fig. 3e) and also the separate metric of Storey’s π1 (raw: 0.08, calibrated: 0.006)39.

In summary, airqtl provides an efficient and scalable solution for cell type-specific sceQTL mapping, uniquely enabling comprehensive benchmarking and objective optimization of statistical performance.

Efficient sceQTL mapping enabled cGRN inference
Genetic variations are effective instrumental variables for inferring cGRNs via Mendelian randomization (MR), given their predominantly cis-regulatory effects. By identifying instrumental (variable) SNPs and their significant cis- and trans-genes, causal regulations from cis- to trans-genes can be inferred (Fig. 4a). Population-scale scRNA-seq studies offer unprecedented power for cell state-specific cGRN inference but also pose substantial scalability challenges. Airqtl’s unique efficiency provides a solution.

Population admixture is a well-known source of false positives in GWAS40, yet its impact on eQTL mapping and cGRN inference remains less understood (Fig. 4b). The Randolph et al. dataset19 includes admixed populations with strong inter-chromosomal genotype correlations (Supplementary Fig. 7). While principal component (PC) covariates of the genotype matrix can capture much of this correlation (Supplementary Fig. 7), they may also introduce spurious associations via collider bias41. An alternative approach involves estimating a genetic relationship matrix (GRM) from genotype data42. To assess these methods, we compared two donor-level GRMs (identity matrix versus genotype-based GRM) with and without PC covariates. Using SNPs with highly insignificant cis-associations (all raw P > 0.5) as negative controls for trans-sceQTL mapping on different chromosomes (Methods), the vanilla approach (identity GRM without PC covariates) outperformed other methods in 7/12 cell states (cell type-flu exposure combination, Fig. 4c, and Supplementary Fig. 8). Notably, genotype-based GRMs had almost identical performances with/without PCs in the QQ plots. This indicated little impact from adding PC covariates, possibly because GRMs already effectively account for admixed population structure. Airqtl’s speed enabled mapping of the trans-effects of insignificant SNPs, warranting further investigations into population admixture modeling.

Therefore, we used the vanilla model to map sceQTLs for each cell state (Supplementary Data 3). Strong cis-sceQTL SNPs (BH Q < 0.1) were more likely to exhibit trans-associations and had more significant trans-genes (Fig. 4d, e, BH Q < 0.5), supporting their utility for cell state-specific cGRN inference. To minimize pleiotropy-induced false positives, possible pleiotropic SNPs that may affect multiple cis-genes were excluded (Methods). Each remaining SNP, namely “instrumental SNP”, could then infer causal gene regulations from its cis-gene to associated trans-genes via two-stage least squares (2SLS, Methods).

Multiple instrumental SNPs allow to evaluate and improve the robustness of inferred cGRNs. Combining gene regulatory pairs (namely gene regulations) identified for each cell state, we identified 28,171 gene regulations supported by ≥ 2 unique instrumental SNPs (Fig. 4f). Among these, only 1 (0.004%) displayed conflicting directionality predictions (activation vs. repression, Methods). Additionally, gene regulation strength estimates were highly consistent across different instrumental SNPs, with 75% exhibiting very low relative variability (< 10%, Fig. 4g, Methods). These findings reaffirm the robustness of cGRNs inferred from population-scale scRNA-seq data. Our final cell state-specific cGRNs consisted of gene regulations supported by at least two unique instrumental SNPs, whose regulation strength estimates were aggregated with median (Supplementary Data 5).

Intrinsic and extrinsic drivers of cGRN state specificity
Single-cell cGRNs offer considerable advantages over noncausal GRNs for understanding molecular circuit rewiring, as they directly quantify perturbation outcomes and indirectly link to organismal trait outcomes. To delineate distinct modes of cGRN rewiring—including detectability, directionality, and strength changes—we first performed pairwise comparisons across cell states (Fig. 5a, Methods). Detectability varied widely between cell states, reflected in differing overlap rates. These variations stem from both genuine (in)activation of specific gene regulations and technical factors such as cell count. Among overlapping regulations, directionality was highly consistent, with differences confined to a small fraction (up to 3%), which could arise from chromatin-level regulation rewiring. We further quantified changes in regulation strength with regulation scale and Pearson R, reflecting differences in overall magnitude and individual gene regulations, respectively, akin to our benchmarks for sceQTL effect size estimation (Fig. 2e, f and Fig. 3c, d, Methods). Pearson R consistently exceeded 0.8 between all cell states, though regulation scale exhibited up to a two fold variation (Fig. 5a). Based on these metrics, cGRNs from the same cell type but different conditions were most similar, followed by those from identical conditions (Fig. 5a–c, Supplementary Data 4). These findings suggest limited cGRN rewiring in regulation strength and directionality, driven predominantly by cell type and followed by external condition.

To dissect genuine (in)activation of specific gene regulations from technical detectability variations, we next compared CD4+ T cells under flu vs. NI conditions, the most abundant cell states in the dataset with the highest numbers of trans-eQTLs (Fig. 4d, e). Under the flu condition, we identified 10,195 regulations involving 1,938 genes (Supplementary Fig. 9, Supplementary Data 5), of which 8,539 (84%) were cell state-specific and absent in NI. Notably, most regulators (641/681 or 94%) were not known TFs (based on known motifs4, Methods), highlighting airqtl’s unique capacity for de novo inference of cell state-specific cGRNs. A force-directed layout segregated genes in the cGRN into two distinct clusters that repressed each other (Methods), with upregulated genes enriched in one cluster (Supplementary Fig. 9).

We further focused on the top 10 master regulators with the highest number of target genes. These master regulators segregated into two distinct groups, each connected exclusively by activation edges within the group and repression edges between groups (Fig. 5d). Notably, 50/51 (98%) regulations among these master regulators were specific to the flu condition. In contrast, the top 10 master regulators under the NI condition exhibited only 21/89 (24%) regulations as cell state-specific (Supplementary Fig. 10). Furthermore, 5/6 master regulators within the upregulated group (TRIM38, STAT1, RBCK1, LNPEP, and TRIM22) had no detectable target under the NI condition. Since all master regulators and many target genes remain expressed in both conditions (Supplementary Data 5), the cell state specificity is consistent with changes in gene regulation (cGRN edge) instead of (in)activation of individual genes (cGRN node). A plausible explanation is that the NI condition predominantly maintains house-keeping regulations while lacking the key gene regulations activated under the flu condition. These findings underscore that active cGRNs are highly specific to the cell state, including both cell type and condition, thereby necessitating cell state-matched data for accurate characterization.

Most of the identified top master regulators43–53 and their target genes are well-established components of cellular response to viral infections. The activation or repression targets of most (7/10) master regulators were most significantly enriched for the gene ontology (GO) category “defense response to virus” (Bonferroni hypergeometric P values < 0.05, Fig. 5e, Supplementary Data 6, Methods). Importantly, the activation/repression target gene set exhibiting this enrichment aligned with the regulatory groups of each master regulator shown in Fig. 5d. In contrast, target genes regulated in the opposite direction showed no statistically significant enrichment in any GO category. These enrichments were driven by diverse target genes, with ~20% unique to a single master regulator (Supplementary Fig. 11). These findings confirm that the inferred cGRNs and their directional regulation information are biologically meaningful and highly specific to the cell state.

Finally, we validated the cell state-specific cGRNs using published perturbation experiments in CD4+ T cells. We first queried the whole Gene Expression Omnibus (GEO) for datasets targeting any master regulator and identified one dataset54 targeting STAT1, which also happened to match the viral exposure cell state. Among their overlapping targets, the cGRN’s estimated effect sizes for STAT1 were highly consistent with the perturbation experiment, with Pearson R = 0.82 and P = 5 × 10^-7 (Fig. 5f). We further found one STAT1 Chromatin Immunoprecipitation sequencing (ChIP-seq) dataset in human CD4+ T cells in the whole GEO55. Using this dataset, we identified STAT1’s putative direct targets through active binding sites (Methods). These exhibited substantial overlap with STAT1 targets in our cGRN (Supplementary Fig. 12, Supplementary Fig. 13, Supplementary Data 7), while the partial overlap confirmed airqtl’s ability to capture indirect targets undetectable by ChIP. Together, these results highlight the capability of our method to infer cell state-specific cGRNs that are experimentally validated and directly reflect perturbation outcomes.

Methods
Airqtl
Array of Interleaved Repeats (AIR). AIR is a compressed data structure optimized for matrices containing repeated rows or columns, a common pattern in hierarchical data structures. Rather than fully expanding matrices, AIR maintains: i) the compressed matrix containing only one copy of each repeated rows/columns and ii) the number of repeats. This design enables faster matrix operations (addition, multiplication, singular value decomposition or SVD) by operating directly on compressed matrices. The acceleration leverages several mathematical properties including computing compressed results, variance rescaling, and other operand preprocessing. Based on and compatible with torch.Tensor, AIR executes all matrix operations within pytorch, simultaneously benefiting from AIR’s scalable algorithms and pytorch’s efficient GPU acceleration.

See Supplementary Note for algorithmic details. See “SceQTL mapping—AIR accelerations” in Methods below for AIR’s application to sceQTL mapping acceleration.

SceQTL mapping—linear mixed model. For sceQTL mapping, we employed the LMM described in Fig. 1c using restricted maximum likelihood (reML)28,29. Briefly, for each gene, λ was first estimated as bλ via maximum likelihood under the null hypothesis (β = 0 for all SNPs). Genes with over 75% of unexplained variance attributable to donor-level variation (i.e., bλ > 3) were excluded from downstream analyses, as such high donor-level variation relative to scRNA-seq measurement noise likely indicates strong technical batch effects. Then, for each SNP, the reML was computed by fixing λ = bλ and optimizing all other parameters (α, β, σ) to maximize the likelihood. The parameter estimates were obtained through reML. Hypothesis testing was conducted using the likelihood ratio test (LRT) to compare the full model against the null model under the same restriction λ = bλ. In practice, we equivalently tested conditional R2 (coefficient of determination) to obtain P values after transforming the data with the GRM and removing covariates. Under the null hypothesis, the conditional R2 follows a beta distribution, which can be well-approximated by a χ2 distribution for computational efficiency (default in airqtl).

Cell state specific sceQTL mapping—linear mixed model. We used the following model to map the cell state specificity of each sceQTL effect size:
y = α0 + α1(C, s, x) + α2(C, s) ⊗ (C, s) + α3 C ⊗ x + α4 x ⊗ x + βs ⊗ x + ε
(1)
This model incorporates all linear and quadratic interaction terms between genotype, cell state, and other covariates. Here, x ∈ {0, 1, 2}^n_cell and y ∈ R^n_cell denote bi-allelic genotype and normalized expression level data for the SNP and gene considered, subjected to the same QC and normalization steps as in sceQTL mapping. C ∈ R^n_covariate × n_cell and s ∈ {0, 1}^n_cell include all covariates, with s specifically containing the cell state covariate (e.g., cell type) being tested for eQTL specificity. α0 is the intercept, and parentheses indicate row-wise matrix concatenation. Covariates C include scRNA-seq technical covariates from Normalisr5 to capture library size confounding and any other data-specific or user supplied covariates (e.g., see “Latent dimensions as covariates” below). Coefficients α0 ∈ R, α1 ∈ R^n_covariate+2, α2 ∈ R^(n_covariate+1)(n_covariate+2)/2, α3 ∈ R^n_covariate, α4 ∈ R are effect sizes from covariates and β ∈ R is the effect size from the tested term of cell state specific sceQTL.

Here we defined the modified row-wise Kronecker product ⊗ satisfying the following conditions. i) For A ∈ R^n_A × k and B ∈ R^n_B × k, A ⊗ B ∈ R^(n_A * n_B) × k computes the element-wise product between each row pair of A and B. ii) Vectors are treated as 1-row matrices: A ∈ R^n_A × k and x ∈ R^k, A ⊗ x ∈ R^n_A × k regards x as a matrix in R^(1 × k). iii) Self-product A ⊗ A ∈ R^(n_A(n_A+1)/2) × k includes only unique row pairs: A ⊗ A = (A1 ⊗ A1, ..., A1 ⊗ An_A, A2 ⊗ A2, ..., A2 ⊗ An_A, ..., An_A ⊗ An_A). iv) Note {0, 1} ⊗ {0, 1, 2} ⊗ R.

The term s ⊗ x captures genotype-by-cell state interactions, with β quantifying the cell state-specific eQTL effect. We incorporated all other quadratic interaction terms as covariates, which offer several advantages for statistical inference. i) While confounders should be included to reduce false positives, the benefits of including non-confounder covariates, such as “prognostic covariates” that explain variance in the response variable, are less known. Prognostic covariates reduce residual variance and improve statistical power, even when not confounding s ⊗ x57. One example prognostic covariate is library size, which is known to confound scRNA-seq measurements of gene expression at experiment- and cell state-dependent levels. ii) Our comprehensive inclusion of quadratic terms follows established recommendations for covariate pre-selection (even if over-selection) before any statistical inference is drawn58. iii) On contrary, covariate selection based on association with target variable using the same dataset risks the statistical complexities of post-selection inference59.

Other practical aspects include replacing x ⊗ x with the presence of a dominant allele, x(dom) = (x(dom)_1, ...), where x(dom)_i = 0 or 1, depending on whether x_i ≠ 0. Since humans are biallelic (x_i = 0, 1, 2), x ⊗ x and x(dom) span the same linear space together with x and the intercept.

The residual variance ε ∈ R^n_cell was modeled as:
ε ~ MVN(0, Σ), Σ = σ^2(I + λK)
(2)
where K ∈ R^n_cell × n_cell is the GRM, same as in Fig. 1c.

Hypothesis testing for the null hypothesis β = 0 was conducted using reML and the LRT, as in sceQTL mapping.

SceQTL mapping—AIR accelerations. In sceQTL mapping, AIR is used to store genotype data for individual cells in population-scale scRNA-seq studies. With its unique exact accelerations for basic matrix operations (Supplementary Note), AIR further accelerates LMMs on top of state-of-the-art implementations for each sceQTL mapping step:
• Expanded genotype matrix. Conventional approaches expand the donor-level G ∈ R^n_SNP × n_donor to cell-level G(E) ∈ R^n_SNP × n_cell explicitly at time complexity O(n_SNP * n_cell). AIR instead stores the compressed G with cell counts per donor, eliminating expansion costs at O(n_donor).
• GRM estimation. Conventional methods estimate the cell-level GRM via matrix multiplication on the expanded genotype (after mean and variance normalization) as K = (G(E))^T G(E) ∈ R^n_cell × n_cell. This incurs a time complexity of O(n_SNP * n_cell^2). AIR leverages efficient matrix multiplication (Supplementary Note) and reduces the time complexity to O(n_SNP * n_donor^2).
• GRM factorization. A major optimization introduced in FaST-LMM and GEMMA is the eigendecomposition (via SVD) of GRM and its subsequent transformation of other matrices such as genotype, gene expression, and covariates. Conventional methods decompose the GRM K ∈ R^n_cell × n_cell in O(n_cell^2 * n_donor) time, even after recognizing the GRM has a rank of n_donor due to row and column repeats. Because all cells from the same donors share the same germline genotypes, their corresponding rows and columns in K would be identical. This allowed AIR to merge these duplicate rows and columns with rescaled variance, further reducing eigendecomposition or SVD time complexity to O(n_donor^3) (Supplementary Note).
• GRM transformation. Transforming the genotype and expression matrices with GRM eigenvectors, as required in the above optimization, involves matrix multiplication. Conventional methods compute this with a time complexity of O((n_SNP + n_gene) * n_cell * n_donor) after accounting for the GRM’s low rank. Using AIR’s efficient matrix multiplication (Supplementary Note), this complexity further is reduced to O((n_SNP + n_gene) * n_donor^2 + n_gene * n_cell).
• SNP-gene covariance. The naive covariance calculation between genotype and gene expression (without considering GRM) takes O(n_SNP * n_gene * n_cell) time to compute by conventional methods. AIR accelerates this computation to O(n_gene * (n_SNP * n_donor + n_cell)) by leveraging its efficient matrix multiplication (Supplementary Note).

SceQTL mapping—normalization with Normalisr. Normalization of scRNA-seq read counts and covariates was performed following the GSE120861 tutorial for Normalisr (1.0.0). Briefly, the pipeline included the functions lcpm, normcov, scaling_factor, compute_var, and normvar from normalisr.normalisr. These steps transformed scRNA-seq read counts into the minimum mean square error (MMSE) estimator of log gene expression levels. Technical variations were extracted from the scRNA-seq read counts as additional covariates (in addition to the constant 1 intercept) for inclusion in the linear mixed model during sceQTL mapping and for variance normalization.

SceQTL mapping—P value calibration. We calibrated P values for cell state-specific sceQTL mapping using null P values derived from a separately simulated null dataset (Supplementary Fig. 5). The null P value distribution was modeled as Beta(e^α0+α1x, e^β0+β1x). Unlike existing approaches37,38, we introduced x, one of the candidate parameters detailed below, computed separately for each gene or SNP to account for potential P value distribution biases. Before calibration, x was normalized to have zero mean and unit variance, and this normalization was stored for subsequent steps. Parameters α0, α1, β0, β1 were estimated via maximum likelihood using stochastic gradient descent in PyTorch, with a learning rate of 10^-6 over 20,000 epochs. To prevent divergence during estimation, P values smaller than 10^-300 were clipped to 10^-300. For calibration, x was computed for the original dataset, normalized as stored, and the P value was adjusted by computing the cumulative distribution function (CDF) of the Beta distribution using the parameters learned above.

Candidate parameters for x:
• Single-cell differential expression technical covariates: lcpm (pseudo-bulk log(CPM+1) for each gene), n0 (proportion of cells with zero reads for each gene);
• Derived parameters (applied to each of the above covariates as x): xmin (minimum of x’s values computed separately for each cell state), xdiff (difference in x’s values between two cell states), xdiffabs (absolute difference in x’s values between two cell states);
• Population genetics parameters: maf (minor allele frequency for each SNP), l0 (bλ for each gene).

CGRN inference—sceQTL mapping and filtering. For cGRN inference based on sceQTL summary statistics, we first performed full cis- (within 1Mbp window before or after the transcription start site) and trans-sceQTL mapping, limiting trans-sceQTL candidate outputs to those with P < 10^-4 for computational efficiency. Gene names containing ‘-’, ‘.’, or ‘_’ were excluded. SNPs not identified as cis-sceQTLs for any gene (BH Q < 0.1, computed separately for cis-sceQTL candidate pairs) were removed. SNPs considered potentially pleiotropic were excluded from cGRN inference unless associated with a single cis-gene (BH Q < 0.1) while associations with all other cis-genes were highly insignificant (raw P > 0.1). Among the remaining SNP-gene pairs, BH Q values were calculated separately for trans-sceQTL mapping. To accomplish this, we implemented a custom function for calculating BH Q values for P values below a cutoff. Significant cis-sceQTLs (BH Q < 0.1) and trans-sceQTLs (BH Q < 0.5, indicating at least half of identified trans-sceQTLs were true positives) were retained as SNP-gene pairs for cGRN inference.

CGRN inference—effect size estimation. To estimate the effect size of each candidate causal gene regulation in each cell state, we looked for instrumental SNP → cis-gene (regulator) → trans-gene (target) relationships. Here the SNP serves as the instrumental variable, as described and identified as a cis-sceQTL for the regulator and a trans-sceQTL for the target above. Using each instrumental SNP, we estimated the gene regulation effect size as the ratio between trans- and cis-sceQTL effect sizes following 2SLS.

CGRN inference—inference robustness evaluation. To evaluate inference robustness, we tested how much the sign and value of estimated gene regulation effect sizes vary across instrumental SNPs. We defined the relative variability of gene regulation strength estimates as their difference between maximum and minimum values divided by the twice of maximum of absolute value. This metric gives 0 when all instrumental SNPs yield identical estimates, 0.5 when one SNP estimates no effect while the others preserve directionality, and 1 when the maximum and minimum estimates have equal absolute values but opposite signs.

CGRN inference—filters with instrumental SNPs. Only gene regulations with at least 2 unique instrumental SNPs were included in the cGRN. Here uniqueness is defined as having non-identical genotype data vectors between the two instrumental SNPs among the donors used to map cell state-specific sceQTLs. Differences merely due to opposite definitions of reference/alternative alleles were regarded as identical.

Data
Published population-scale scRNA-seq datasets. The Randolph et al. dataset19 contains scRNA-seq data of 255,731 peripheral blood mononuclear cells (PBMCs) sampled from 89 donors of African and European ancestry. The PBMCs were subjected to two distinct conditions: mock treatment (NI) and infection with the H1N1 Cal/04/09 influenza A virus (flu). Donors were genotyped with whole genome sequencing. Cell type and donor annotations provided in the original dataset were retained for this study. The Seurat object was provided in the original study and directly exported into the input files for airqtl without any processing. The code for reformatting is available on Zenodo under accession code 15925067.

Quality control. We primarily followed the quality control established for high-MOI single-cell CRISPR screen in the Normalisr study5, with additional criteria for donor-based statistics. Briefly, we removed cells with less than 100 expressed genes or 500 reads, genes expressed in less than 10% of donors, 100 cells, or 2% cells, (cells from) donors with 5 cells or fewer, and SNPs that could not be confidently identified for any donor in the original study.

The above removal was performed separately for each cell subset (e.g., by cell type according to existing annotations) in published (and possibly already QCed) population-scale scRNA-seq data and iteratively until no gene, cell, SNP, or donor was removed. Then, cell subsets with less than 500 cells, 40 donors, 10,000 SNPs, or 500 genes remaining were excluded from sceQTL mapping and downstream analyses. Here specifically, a gene was defined as “expressed” if it exhibited at least 1 read in a particular cell or any cell from a particular donor.

Simulation method for sceQTL mapping data. Our simulation strategy employed permutation and resampling to maximize the similarity with real datasets for benchmarking purposes. To generate a simulated dataset with n_SNP SNPs, n_gene genes, n_donor donors, and n_cell cells that mimics a real dataset, we followed the approach established in Normalisr5 and introduced the contributions of cis-eQTLs. First, to resample these dimensions, we removed non-autosomal genes and randomly downsampled to n_gene genes from the remaining ones. SNPs outside 1Mbp window of the transcriptional start site (TSS) of any remaining gene were removed and then downsampled to n_SNP randomly to maximize cis-eQTL simulation efficiency. Then, donors and cells were randomly downsampled or upsampled to n_donor and n_cell, respectively.

For cis-eQTL simulation, each SNP within a 1Mbp distance of any gene’s TSS was randomly selected as a cis-eQTL for that gene with a 50% probability. Cis-eQTL effect sizes were then simulated by performing multivariate differential expression using Normalisr on the genotypes of all simulated cis-eQTLs in the resampled dataset. For the log expression of gene i normalized via Normalisr as y_i (see “SceQTL mapping—normalization with Normalisr”), the maximum likelihood estimator for the effect size of SNP j was computed as β_ji within the linear model:
y_i = Σ_k α_ki c_k + β_ji x_j + ε(0)_i
(3)
where y_i represents the normalized expression of gene i, c_k denotes the normalized covariate k (including the intercept), α_ki is the effect size of covariate k on gene i’s expression, x_j is the genotype of SNP j, and ε(0)_i ~ N(0, σ(0)^2) is the cell-level error term. To account for effect size underestimation caused by genotype correlations from linkage disequilibrium (LD) and to extend the range of effect sizes covered in benchmarking, β_ji was scaled by a factor of 1.5. The α_ki term for the intercept was retained to preserve the mean expression level of gene i.

To simulate genotypes, the existing genotype values (0, 1, 2) in the resampled dataset for each SNP were randomly permuted across donors to remove LD. This allowed to benchmark false positives in cis-eQTL mapping while preserving the minor allele frequency. In this study, we regard LD a question to be addressed in fine-mapping, beyond the scope of the sceQTL mapping step. Additionally, since neither this study nor CellRegMap aimed to enhance the statistical accuracy of GRM estimation, the use of permuted genotypes provided an ideal test scenario for assessing their core functionalities without introducing unequal genetic relatedness among donors.

To simulate the (unnormalized) true log expression levels for gene i in each cell, y_i, we extended the log expression model from Normalisr:
y_i = α_i + Σ_j β_ji x_j + ε(0)_i + ε(1)_i
(4)
where α_i represents the intercept and mean expression level for gene i, x_j denotes the genotype of SNP j, β_ji is the effect size of SNP j on gene i, ε(0)_i ~ N(0, (σ(0))^2) is the cell-level error term, and ε(1)_i ~ MVN(0, (σ(1))^2 K) is the donor-level error term. Here, K is a block diagonal matrix of all-one square matrices, each sized according to the number of cells from a given donor. The intercept α_i was estimated as the pseudo-bulk mean expression level, and β_ji was estimated during the cis-eQTL simulation step. We set σ(0) = 0.5 and σ(1) = 1 to ensure a substantial contribution from each source of variation.

To simulate scRNA-seq read counts, we first calculated the normalized true expression level as the softmax across all genes for each cell. This step provided gene-level normalization and incorporated the exponential transformation from the unnormalized log expression level y_i. Then, scRNA-seq read counts for each cell were sampled randomly across all genes using a binomial distribution, with probabilities proportional to the normalized true expression levels. The total number of reads per cell was randomly sampled from its distribution derived from the resampled dataset.

Details about our simulation method can be found in airqtl.sim at (https://github.com/grnlab/airqtl).

Simulation method for cell type-specific sceQTL mapping data. The simulation of cell type-specific sceQTLs was based on the sceQTL simulation method with the following modifications: i) Cell type (CD4 vs. CD8 T) was included as an additional one-hot encoded covariate throughout the simulation process. Its effect size on gene expression levels was estimated in α_ki in Eq. (3) and retained for subsequent simulation. Technical covariates introduced by Normalisr were additionally normalized in advance to be orthogonal to cell type covariates, ensuring accurate effect size estimation. ii) Each cis-eQTL was assigned a 50% probability of having cell type-specific effects. For genes with cell type-specific cis-eQTLs, their expression level model (Eq. (3)) incorporated additional linear terms derived from the interaction between cell type and genotype. The corresponding effect sizes were estimated and stored for use in the simulation. iii) The simulation of true log expression levels (Eq. (4)) was updated to account for cell type-specific contributions, consistent with the modified effect size estimation model. Specifically, the two extra terms above were added along with contributions from technical covariates from Normalisr and the intercept. Their effect sizes were introduced as they were estimated above. Since these terms already characterized the mean expression levels, the α_i term was removed.

Simulation method for calibration. The simulation method for null calibration datasets of cell type-specific sceQTL mapping followed the same steps as cell type-specific sceQTL mapping benchmarking datasets but with the modifications: i) Simulations were based on the dataset for which cell type-specific sceQTL mapping P values were to be calibrated, with identical dimensions (n_SNP, n_gene, n_donor, and n_cell). Because simulated datasets were used for benchmarking, they were the input datasets for this simulation to mimic. ii) Genotypes were not permuted, as LD removal was only required for benchmarking. iii) Cis-eQTLs were present but none of them was cell type-specific. iv) The β_ji values were scaled by a factor of 2 (instead of 1.5) to ensure adequate coverage of effect size ranges for calibration purposes.

Simulated datasets for benchmarking. To benchmark the efficiency of sceQTL mapping and cell type-specific sceQTL mapping, we generated a series of simulated datasets designed to mimic the CD4+ T cells (the largest cell subset for sceQTL mapping) or the CD4+ and CD8+ T cells (the largest and most biologically related cell subsets, for cell type-specific sceQTL mapping) from the quality-controlled real dataset by Randolph et al.19. Small-scale datasets, designed to accommodate the limited scalability of CellRegMap, had default dimensions of n_SNP = 500, n_gene = 500, n_donor = 100, and n_cell = 1,000 for sceQTL mapping, and n_SNP = 200, n_gene = 1,000, n_donor = 50, and n_cell = 2,000 for cell type-specific sceQTL mapping. Large-scale datasets had default dimensions of n_SNP = 4,000,000, n_gene = 5,000, n_donor = 100, and n_cell = 50,000 for sceQTL mapping, and n_SNP = 200,000, n_gene = 2,000, n_donor = 100, and n_cell = 50,000 for cell type-specific sceQTL mapping. Using these default configurations, we also generated related datasets by varying each of these dimensions to assess the scalability of each method under different scenarios.

For statistical performance benchmarking in sceQTL mapping or cell type-specific sceQTL mapping, we similarly generated one simulated dataset each with n_SNP = 20,000, n_donor = 100, n_cell = 50,000, and all post-QC genes.

In all cases, datasets were generated by regarding each donor-condition pair from the real dataset as a separate donor. This allows to expand the range of variability in the number of donors. Genotype permutations for LD removal also eliminate any residual population structure, ensuring that each simulated donor is distinct.

Benchmarking
Other methods. For CellRegMap, we followed its website and tutorials for single-cell preprocessing and association testing starting from raw read counts (https://github.com/annacuomo/CellRegMap_analyses/tree/4f26437d20e824e704385deeb2bec1d9f1f61f4e/endodiff/preprocessing, https://github.com/annacuomo/CellRegMap_analyses/blob/4f26437d20e824e704385deeb2bec1d9f1f61f4e/simulations/utils/run_tests.py, https://limix.github.io/CellRegMap/usage.html). To allow a consistent benchmarking, latent dimension covariates were not included (see below). For statistical performance benchmarking (below), only SNP-gene relations in the cis region (1Mbp up or downstream) were computed for CellRegMap’s efficiency.

We did not include GASPACHO34 and other Bayesian methods that are less compatible with downstream analyses such as fine mapping, transcriptome-wide association studies, and Mendelian randomization. We did not include BOLT-LMM32 or other heuristic methods that rely on the same approximation such as SAIGE-QTL24 for time complexity comparison or benchmarking. Their approximation is designed for GWAS and are prone to large errors without a large number of individuals (BOLT-LMM recommends 5000 samples), which is rarely met by current population-scale scRNA-seq studies. In addition, we could not run GASPACHO on custom datasets due to the lack of instruction on preparing input parameters. Our test run of SAIGE-QTL’s tutorial also met with critical errors. Nevertheless, airqtl is estimated to run approximately 10^3 times faster than SAIGE-QTL based on its projected running time in ref. 24. Our AIR data structure can also provide similar accelerations to these methods.

Estimation of realistic acceleration from AIR’s algorithm. For estimation purposes, a realistic sceQTL mapping task was assumed to have 100,000 cells, 1,000 donors, 8,000 genes, and 4,000,000 SNPs post-QC for each most abundant cell type because they take up the majority of computation time. These numbers were put into the computational complexities, namely (n_SNP + n_gene) * (n_donor + n_cell) * n_cell for FaST-LMM and GEMMA and n_SNP * n_gene * n_donor + (n_donor + n_cell) * (n_donor + n_cell) * n_donor for airqtl. Therefore, their ratio yield the estimated acceleration from AIR’s algorithm.

Latent dimensions as covariates. In our benchmarking, we did not account for latent dimensions (such as those derived from PCA) from genome and/or transcriptome data as covariates for several reasons. i) This study is focused on improving computational efficiency without sacrificing statistical performance. Identifying the optimal strategy for incorporating latent factors would be the subject of a separate study. ii) While latent dimension covariates may improve false discovery rate control in sceQTL mapping, they can also obscure true signals, reduce statistical power, and lead to spurious associations41. Although the impact of these drawbacks varies across dimension reduction methods and data types, we would ideally prefer to identify a method that improves false discovery rate control without introducing these trade-offs. iii) Our benchmark for sceQTL mapping employed independently permuted genotypes, which lacked population structure, and yielded satisfactory results without the need for latent dimension covariates. iv) Airqtl is based on a standard linear mixed model and offers the flexibility to incorporate any type of covariate, including latent dimensions, if necessary, as demonstrated using the real dataset from Randolph et al.

Efficiency benchmarking. We compared the efficiency of CellRegMap and airqtl on Google Cloud Compute Engine, using the machine type for which each method is designed, and recorded the running times with Snakemake. Specifically, airqtl was run on a 2-highgpu-1g with 1 NVIDIA A100 40GB GPU. Since each CellRegMap process was single-threaded, we ran up to 32 CellRegMap processes in parallel on n2d-highmem-32 with 32 vCPUs. No resource competition for vCPUs or memory was observed among the 32 parallel processes. The computing cost was estimated by multiplying the running time by the spot price of the Google Cloud Compute Engine instance ($1.469/hour for a2-highgpu-1g and $0.668/hour for n2d-highmem-32). The cost for CellRegMap was further divided by 32, as each process utilized only 1 vCPU, and the price was linearly proportional to the number of vCPUs. The datasets described in “Simulated datasets for benchmarking” were used for these evaluations.

Since CellRegMap could not scale to the large datasets, we projected its running time and associated cost based on its performance on smaller datasets, using the following linear model:
log t = α + β_SNP log n_SNP + β_gene log n_gene + β_donor log n_donor + β_cell log n_cell + ε
(5)
where ε ~ N(0, σ^2). For sceQTL mapping, all data points for CellRegMap in 2a were used to train this model and project the running times for CellRegMap in 2b. For cell type-specific sceQTL mapping, all data points for CellRegMap in Fig. 3a were used to train the model and project the running times for CellRegMap in 3b.

Statistical performance benchmarking. We used the datasets described in “Simulated datasets for benchmarking” for statistical performance benchmarking. We benchmarked the estimation accuracy of sceQTL effect sizes and cell type-specific sceQTL effect sizes with bias and variance, following 5. Briefly, bias characterizes the overall over- or under-estimation of effect sizes and can be quantified in the linear model:
y = βx + ε, ε ~ N(0, σ^2)
(6)
where y represents the estimated effect sizes, and x denotes the ground-truth effect sizes for cis-eQTLs in the simulated dataset. Therefore, β reflects the estimated effect size scale. The maximum likelihood estimator for β was computed as βb, and its deviation from 1 indicates systematic bias. Variance was estimated as the deviation of Pearson R from 1 between the estimated and ground-truth effect sizes. Given that expression-dependent underestimation of effect sizes is commonly observed in single-cell differential expression5, a simpler version of sceQTL mapping, we performed this benchmarking both for all genes together and for genes stratified by expression level (pseudo-bulk log CPM) to assess the severity of this problem in sceQTL mapping and cell type-specific sceQTL mapping.

To benchmark the specificity of each method, we separately compared the sceQTL P value distribution for non-eQTLs or non-specific sceQTLs against the uniform distribution. To assess the sensitivity of each sceQTL mapping method, we compared the P value distribution of eQTLs against those of non-eQTLs. Similarly, to evaluate the sensitivity of each cell type-specific sceQTL mapping method, we compared the P value distribution of cell type-specific eQTLs against those of non-eQTLs and non-specific eQTLs.

Type I error evaluations. We drew P value histograms for sceQTL mapping or cell type-specific sceQTL mapping P values for SNP-gene pairs categorized as non-eQTLs (SNP has no effect on gene expression), eQTLs (SNP has a nonzero effect on gene expression), non-specific eQTLs (SNP has a nonzero effect on gene expression independent of cell type), or cell type-specific eQTLs (SNP has a nonzero effect on gene expression dependent on cell type). The ground-truth eQTL status of these pairs was determined in the simulation data. Each histogram was divided into 20 equal bins from 0 to 1. Error bars were estimated as 3 * sqrt(N), where N is the number of entries in each bin. Dashed lines represent the uniform distribution.

We used two-sided Kolmogorov-Smirnov test against the standard uniform distribution to evaluate the specificity of each method. We also used Storey’s π1 estimated by fdrtool60 (1.2.18) to characterize the proportion of false positives identified by each method or configuration in fully null datasets.

State-specific cGRN inference and analysis
Population admixture detection and removal. We performed a Pearson correlation test between SNP pairs located on different chromosomes to avoid correlation due to linkage disequilibrium. To accelerate computation, only 100,000 post-QC SNPs were randomly sampled for this test. As a baseline, we also performed this test after permuting the donors separately for each SNP.

To test the effect of population admixture removal on genotype correlation, we introduced different numbers of top PCs of the genotype matrix for the tested SNPs as covariates for the (conditional) Pearson test.

To test the effect of population admixture removal on trans-sceQTL mapping, we used GRM estimated from the genotype matrix and/or introduced its the top principal component from the genotype matrix as a covariate without LD pruning.

Comparison of cell state-specific cGRNs. The overlap rate was first computed as the ratio between the number of overlapping gene regulations and the minimum number of gene regulations between the two cGRNs being compared. Overlapping gene regulations between two cell states were defined as those present in both cGRNs, irrespective of effect size or directionality estimation. Any deviation of the overlap rate from one may reflect contributions from differences in the cGRNs, differences in the detectability of individual gene regulations, and data randomness due to the finite number of cells.

Among the overlapping gene regulations, three metrics were computed to characterize the actual differences in cGRNs. The directionality difference was defined as the proportion of gene regulations that exhibited opposite signs in the estimated effect sizes. Pearson R was calculated for the estimated effect sizes between the two cell states to assess the overall similarity. Regulation scale was computed as the linear regression coefficient between the estimated effect sizes of the two cell states (inverted if < 1). Since only overlapping gene regulations were considered, these metrics are less influenced by factors other than differences in the cGRNs.

Visualization of cell state-specific cGRNs. We used Cytoscape (3.10.2) for cGRN visualization. Prefused force directed layout was applied on the largest connected component of the cGRN for CD4+ T cells under flu condition using regulation strength as edge weight.

Gene ontology enrichment analysis. For each master regulator, we used GOATOOLS (1.3.11) for GO enrichment analysis on its target genes using all other genes that have a regulator in the cGRN as background. This reduces the bias of background gene selection compared to using all genes or genes expressed in the corresponding cell state because they do not account for cell state specificity of cGRNs. Activated targets (with regulation strength > 0) and repressed targets (with regulation strength < 0) were analyzed for GO enrichment separately. To avoid evaluation biases for our expression-based analyses, for all GO analyses we excluded GO evidences that may also have an expression-based origin (IEP, HEP, RCA, TAS, NAS, IC, ND, and IEA). All three GO categories were used. Both over- and under-enrichments were calculated. Gene names were converted with mygene (3.2.2) when needed. Basic ontology (2024-09-08) and human GO association (2024-09-10) files were used.

Query of perturbation experiments. We queried GEO for perturbation experiments involving any of the 10 master regulators in CD4+ T cells under the flu condition. The experiments were required to involve primary CD4+ T cells or their derived cell lines. Perturbations should include overexpression, knockdown, or knockout of the master regulator, along with a control condition, followed by either bulk or single-cell RNA-seq or microarray measurements. Viral infection or exposure was not a requirement. Only target genes that were significantly identified in the perturbation experiment and also found in the inferred cGRN were included for comparison. Perturbation experiments with ≤ 5 overlapped targets were excluded from the comparison.

Analyses of ChIP-seq experiments. We queried the whole GEO for ChIP-seq experiments of STAT1 in CD4+ T cells in human and found one dataset with two conditions: INFg treated and control.

For each condition, we performed liftover of their annotated STAT1 peaks from hg18 (reference genome in the original study) to hg38 (reference genome in Randolph et al.). We used Homer (version 5.1, scanMotifGenomeWide function) to filter high-confidence STAT1 binding sites by filtering only peaks containing canonical STAT1 motifs. Each gene was then annotated with its nearest STAT1 binding site based on transcriptional start site (TSS) proximity.

For the Venn diagrams, we used a wide range of distance thresholds from 0.02 to 500 kb as a parameter to determine whether each gene was regarded a STAT1 direct target. A hypergeometric test was performed to assess the statistical significance of the overlap, using all genes with at least one regulator gene in the cGRN as the background.

For the box plots, we used 5 kb and 100 kb for thresholds representing proximal and distal regulation activity. To account for potential target gene sharing among flu-condition master regulators, we performed our comparisons using NI-condition master regulators as the baseline reference when assessing overlap with ChIP-seq identified targets.

Conventions
All P values are two-sided raw unless otherwise specified. For box plots, boxes extend between 1st and 3rd quartiles. Lines indicate median. Whiskers extend from the box to the farthest datapoint within 1.5x of the inter-quartile range from the box, unless otherwise stated. Fliers are shown if beyond the whiskers.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Lingfei Wang presents a very efficient implementation for the detection of single-cell gene expression QTLs (sceQTL), which is called Airqtl. At its heart, Airqtl performs greatly accelerated matrix operations utilizing GPU functionalities. AIR stands for Array of Interleaved Repeats, which exploits the fact that all cells from the same donor have the same genotype. Further it uses torch.Tensor to perform the matrix operations in GPUs. Importantly, all of these computations are performed while correcting for population-substructure using linear mixed models (LMM), which is currently state-of-the-art for QTL mapping in human populations.\n\nThe submission presents a significant methodological advancement and it is mostly very clearly written. I only have a few, mostly minor comments.\n\n1. The core of the method should be better described. The paragraph on AIR in the Methods section is too brief to really understand what it is doing. The text in the Methods only says what AIR does, but provides hardly any details on how it does it. One has to refer to the supplement to understand the approach. I suggest to provide a more concrete description of how AIR works and in which way it exploits functionalities of torch.Tensors in the Methods section.\n2. Related to the first point: it should be better explained in which way AIR is different from pseudo-bulking. This has to become more clear already in the Results section (\"Airqtl overview\"). Note that also in Figure 1a, green box \"AIR acceleration\" it is not really clear why this is not simply pseudo-bulking per donor. The matrix on the right hand side (Gene x Donor) looks as if all cells of a donor were collapsed into a single vector. This point also has to be improved.\n3. Line 20: T is missing at the beginning of the sentence\n4. Figure 2 & 3, specify in the description/text (not only in methods) the values at which the constants n were kept.\n5. Figure 3e: First of all, the shades (representing confidence intervals) are hardly visible. Second, the confidence intervals seem very small given the variation (noise?) of the estimates along the x-axis. I can't believe that all of the small 'wiggles' would be 'true', which is suggested by the super small confidence intervals.\n6. Figure 4c: The differences between the models with and without PC1 are hardly visible; colors/dots are overlapping too much. Consider a scatter plot of p-values with versus without the 1st PC.\n7. Page 11: 'Regulation scale' needs to be explained here. Unlike Pearson's R, this is not a commonly known measure.\n8. Page 13, lines 214 - 219: Does the author interpret this as a cell-type dependent 're-wiring' of the network (alternative network topology) or is this rather a cell-type-dependent activation/repression of target genes? In other words: to what extent is the topology of the network actually changing as opposed to altering the activity of different parts of an otherwise static network topology?\n9. It would be very helpful to describe example 1.3 in the Supplementary Note in more detail and apply it to the definition above. Then it would also be easier to understand how the repeats are reduced to one representation.\n10. Suppl. Fig. 7: I am not sure this figure conveys any relevant information. In the main text is reads \"... with upregulated genes clearly enriched in one cluster\". I don't see this from this figure. Maybe a heatmap showing the clustered adjacency matrix might work better?", "reviewer_2": "This manuscript describes a new computational approach, AIRQTL, for single-cell-based eQTL mapping. By changing the strategy for eQTL mapping—usually based on SNP vs. cell and cell vs. gene—to SNP × donor and donor × gene matrices, AIR improves the computational efficiency of sc-eQTL mapping. While the computational improvements are clear, it is not evident how the eQTL estimates produced by the new algorithm advance the eQTL analysis problem. Additionally, the comparative analysis with alternative competing methods is limited and should be expanded and improved.\n\nMain Points:\nI could not find results supporting the statement: “Compared to state-of-the-art LMM methods such as FaST-LMM [27] and GEMMA [28], which operate at O((nSNP + ndonor)(ngene + ncell)ncell) time complexity, AIR delivers a ∼1000-fold algorithmic speedup for sc-eQTL mapping with 100,000 cells and 1,000 donors per cell type.” Currently, only a computational complexity analysis is provided in the supplement and an empirical evaluation with CellRegMap is included. The authors should also contrast their results with other similar approaches. How does the method compare with bulk-based approaches or cell-type bulking approaches? There are at least two other tools for sc-eQTL analysis in the literature (references 33, 34), which should also be discussed and compared in the manuscript.\n\nWhy is CellRegMap not evaluated in Fig. 3c–e? The authors should clarify this omission or include the comparison.\n\nThe authors discuss a conceptual improvement (line 114) potentially arising from the use of Normalizer. They should conduct ablation experiments to support this claim. Additionally, they could demonstrate how this improvement enhances cell-type-specific eQTL detection.\n\nOverall, it is difficult to determine whether AIRQTL offers a computational speed improvement, a conceptual improvement, or both. The manuscript should clarify how AIRQTL compares with current state-of-the-art methods in these respects.\n\nLine 206: Do the results imply that there are 641 proteins potentially acting as transcription factors (TFs), although their functions are currently unknown? The authors should investigate the potential functions of these proteins—e.g., whether they preferentially interact with known TFs or are associated with chromatin-related functions. Similarly, more insight into the functional roles of the proteins in the gene regulatory networks (GRNs) shown in Fig. 5d is needed.\n\nRegarding the GRN in Fig. 5, STAT1 is shown as a potential regulator. Since STAT1 is a TF, the authors are encouraged to investigate whether the predicted interactions are supported by TF binding data (e.g., ChIP-seq or TF binding motif analysis).\n\nThere is a lack of detail on how the input data for AIRQTL was generated from the Randolph et al. dataset (line 383). This section should be expanded and clarified.\n\nWe could run the software, but there is limited documentation and it mostly consists of scripts to re-run the experiments. This should be improved. Here are some additional comments.\nThe Airqtl code is clean and well-structured, with good comments on the parameters.\nThere is no reference for the python version needed to run.\nData format can only be understand via a script. Author is encouraged to make use of Data Containers (AnnData / scanpy). This would make the useability of the method possible for a broader audience.\nThe method uses PCA reduced dimensions of the transcriptome data as covariates in the linear mixed model to control hidden structures. The package requires a covariates matrix as input, which can be obtained from the expression matrix.\nThere is no option/output to visualize cGRNs.\n\nMinor Points:\nIn Fig. 1, it is unclear how AIRQTL obtains cell-type-specific eQTL estimates. This should be clarified in the figure or its legend.\nLine 196 mentions chromatin rewiring. What data supports this claim? Was chromatin directly analyzed? If not, the authors should rephrase this statement to more accurately reflect the observed results (e.g., “rewiring of the GRN”).\nLine 259 should be rephrased. While eQTL analysis offers a robust statistical framework for linking genetic variants to expression traits, the manuscript does not provide or cite supporting evidence for this claim.", "reviewer_3": "This paper presents Airqtl, a novel computational method for mapping single-cell expression quantitative trait loci (sceQTLs) to infer cell state-specific causal gene regulatory networks (cGRNs). The motivation stems from the computational limitation of existing methods that are unable to scale to population-level scRNA-seq studies. By employing a compressed data structure called Array of Interleaved Repeats (AIR) together with GPU and other implementation optimization, airqtl significantly improves computational efficiency, achieving orders of magnitude acceleration compared to existing methods. The scalability of airqtl facilitates cell state-specific sceQTL mapping and allows robust and experimentally validated cGRNs. The author benchmarked airqtl against existing methods using both simulated and real datasets, and the real data analyses inferred cGRNs consistent with established biological knowledge.\n\nThis work makes contributions to sceQTL mapping and cGRN inference. However, I still have some major and minor comments listed below that may be helpful to improve the study.\n\nMajor comments:\n1. Since computation acceleration is a major contribution of this work, the related computational complexity parts need further clarity and accuracy. In the 'SceQTL mapping — algorithmic accelerations' section, 1) for GRM estimation, G(E) is a matrix of n_snp by n_cell, so the complexity of G(E) multiplying its transpose should be O(n_snp^2*n_cell) instead of O(n_snp*n_cell^2). The same issue occurred to GRM for G of size n_snp by n_donor. 2) In general, decomposition of matrix n by n requires a complexity of O(n^3). Can the author elaborate on why GRM decomposition for the n_snp by n_snp matrix G(E)*G(E)^T requires O(n_cell^2 by n_donor) and ARI requires O(n_donor^3)? 3) The GRM transformation and SNP-gene coraiance parts also need further clarification regarding the computational complexity.\n2. As computational complexity was discussed in the paper, fig 2 (a) and (b) are somewhat misleading, as they appear to suggest that the computational cost scales linearly with all dimension parameters. A brief discussion of the results on the practical and theoretical scalability should be added here.\n3. In Fig. 1(a), if E(R) = sum_in_donor(E) denotes summation over all cells for a specific donor, is there information loss here (biologically or statistically), in comparison to using the cell-specific expression data? Can the author elaborate on this point?\n4. In fig. 1, since (d) uses data from (c), the readability can be improved by connecting the two parts with shared notations. For instance, using C to denote the normalized covariates in both (c) and (d).\n5. Comparing fig. 2(a) with 2(b), the second and third columns are a bit redundant as the numbers genes and donors do not change too much. Also, it would be good to clarify the fixed sized when one dimension size changes.\n6. The accuracy of section 'SceQTL mapping — linear mixed model for cell state specificity' needs further improved or at least additional clarification. If A is a matrix of n_A by k, B is a matrix of n_B by k, the Kronecker product is of n_A*n_B by k^2. The author also needs to clarify the self-product case.\n7. Can the author elaborate on why an LMM with all linear and quadratic interaction terms was adopted? Was such quadratic interaction reflected in the result? How to interpret both statistically and biologically?\n8. For equation (1), detailed description, including the dimension of each notation and their practical meaning, is needed to improve the readability of the manuscript.\n\nMinor comments:\n1. On page 2 line 20, there is a typo in 'hese approaches ...'.\n2. For better comparison, fig 2(c) and 2(d) should be combined with different colors denoting different methods."}, "author_response": "We sincerely thank all reviewers for their thoughtful and constructive feedback, which has significantly improved our manuscript. Below, we provide a high-level summary of key shared concerns and our responses before addressing individual reviewer comments in detail. Our response and proposed changes to the manuscript are marked in blue.\n\nShared Key Concerns and Responses:\n1. Clarification of AIR’s novelty and distinction from pseudo-bulking\nAll Reviewers raised questions about how AIR differs conceptually from pseudo-bulking and how it achieves computational acceleration without information loss. We have:\n• Expanded Results, Methods, and Supplementary Note to provide both high-level and technical descriptions on AIR’s exact algorithmic acceleration and its application to sceQTL mapping.\n• Explained four major distinctions between AIR and pseudo-bulking: exact acceleration, statistical power, statistical modelling, and implementation.\n• Revised Figure 1a to highlight AIR’s role in accelerating multiple LMM steps, not just SNP-gene covariance calculations.\n2. Benchmarking and validation of computational improvements\nReviewers #2 and #3 requested additional clarification and validation of AIR’s claimed speedup and comparisons to more methods. We have:\n• Added a new Methods section and revised Results to estimate realistic acceleration provided by airqtl.\n• Clarified limitations in benchmarking other tools (e.g., SAIGE-QTL, GASPACHO) due to technical challenges.\n• Included ablation experiments to compare with LogCPM as baseline normalization.\n\n[... (Author response continues for each reviewer point as provided in the source text) ...]"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have clarified most of my requests and clarified the advances of Airqtl lies on the computational speed up, which allows the proposed analysis for the first time. They also explained (mostly technical reasons) why a broader benchmarking was not feasible. This resolves most of my requests. I do have a pending request regardign the stat1 chip-seq analysis.\n\nSpecific point\nWhile the overlapp analysis (Sup. Fig. 12) suggests some relation between Stat1 targets and ChiP-seq targets; this analysis lacks statistical assesment. Authors could check the overlapp for randonly selected genes (same size as Stat1 targets) and estimate the enrichment with a hypergeometric test (or similar).", "reviewer_3": "The author has addressed all my comments. I do not have any further comments."}, "author_response": "We have included hypergeometric test summary statistics for each panel of Sup. Fig. 12 in Supplementary Data 7 as requested."}]

Decision letter