Characterizing gene perturbations in single cells via network divergence analysis

Domain: neuroscience
Published: 2026-04-11
Editorial outcome: major_revision

Abstract
Functional perturbations of genes do not always cause expression changes, but can manifest through network rewiring or context-specific shifts in regulatory activity. However, inferring gene functional shift and linking them to the specific cell populations remains challenging, as most current scRNA-seq data analysis focuses either on differential gene expression or on cell abundance/state changes, but rarely associate gene perturbations with particular cell populations. Here we present scDNS, a framework that quantifies gene-specific functional perturbations by measuring information-theoretic divergence between condition-specific gene interaction network configurations. In simulated stress tests and multiple experimental datasets, scDNS prioritizes key regulators and the perturbed cell populations-even when mRNA changes are minimal but network rewiring is pronounced. Applications to immunodeficiency mutations, stimulus responses, and viral infection reveal hidden regulatory programs and heterogeneous responder cell states. In pancreatic cancer, scDNS nominates TIMM44 as a mitochondrial sensitizer enhancing gemcitabine efficacy. Together, scDNS provides a powerful tool for inferring dynamic gene perturbations on single cells.

Introduction
Gene interaction networks (GINs) integrate multiple layers of molecular interactions, such as transcription factor-DNA binding, protein-protein interaction, and signaling pathways, to coordinate gene expression programs that define cellular identity, function, and adaptation to environmental stimuli. These networks exhibit inherent dynamism and modularity, wherein perturbations—including genetic mutations, infections, environmental stressors, or pharmacological agents—can alter regulatory relationships and modular structures, thereby reprogramming cellular states and phenotypes1,2. For example, the T95R mutation in IRF4 modifies its DNA-binding affinity, thereby restructuring its target gene network and resulting in defective B-cell maturation, reduced immunoglobulin class switching, and immunodeficiency3. Mutations in TP53 reorganize downstream transcriptional modules, yielding lineage-specific oncogenic outcomes that differ between epithelial and hematopoietic cells4. Pyruvate kinase M2 (PKM2), conventionally regarded as a glycolytic enzyme, undergoes striking remodeling of its interaction network upon activation of the EGFR–ERK pathway. Phosphorylation-dependent nuclear translocation enables PKM2 to relinquish its metabolic role and act as a transcriptional coactivator of β-catenin5. These cases highlight those perturbations extend beyond changes in gene expression, fundamentally restructuring the gene interaction networks that govern cellular processes. Thus, quantifying such network reorganizations will provide a valuable approach to deciphering how perturbations modify molecular circuits, identifying critical regulatory nodes, and elucidating adaptive mechanisms in disease progression and therapeutic responses.

Importantly, further mechanistic interpretation requires more than ranking “perturbed genes” alone: they also demand resolving in which cell states these regulators exert their effects—particularly when the perturbed program is confined to a specific cell subpopulation or manifests heterogeneously across cell states.

Advances in single-cell RNA sequencing (scRNA-seq) have revealed the heterogeneous and dynamic nature of gene regulation across diverse cell types and states. This inherent complexity poses a significant challenge for comprehensively identifying functionally perturbed genes by simply reconstructing a universal GINs and analyzing network reconfiguration6,7. Current strategies for scRNA-seq analysis still predominantly rely on differential expression analysis (DEA) methods coupled with pathway enrichment analysis. Methods, such as Wilcoxon rank-sum or DESeq2-based methods, assess differences in mean expression or variance but fail to capture alterations in regulatory connectivity8,9. As a result, genes with unchanged expression yet profoundly modified regulatory associations often evade detection. Emerging DEA approaches, such as waddR10 and covariance-oriented frameworks like Memento11, capture distributional changes and gene-gene correlations. However, they remain insufficient for quantifying functional perturbations that specifically modulate gene-gene interactions. Similarly, regulon-centric tools like VIPER12 and SCENIC13 depend on predefined, static regulatory models. These methods often underperform in dynamic contexts where regulators and targets exhibit condition-specific rewiring, as exemplified by the IRF4 pT95R mutation3. Critically, even when a candidate regulator is nominated, existing pipelines rarely provide a principled way to map that regulator’s functional perturbation back onto individual cells, leaving the downstream cellular “footprint” of rewiring ambiguous. This disconnection limits mechanistic follow-up, because it becomes unclear whether a putative driver reflects a broad state transition, a rare responder subset, or a context-restricted program within a specific lineage.

Detecting cell populations subject to condition-specific perturbations represents an additional analytical challenge. Current strategies, including Milo14 (which evaluates differential abundance via k-nearest neighbor graphs), scDist15 (which quantifies donor-specific variability using mixed-effects models), MELD16 (which computes perturbation probabilities from density variations in latent space), and Augur17 (which employs machine learning to prioritize cell types based on separability across conditions), offer insights at the population level but primarily focus on compositional or broad transcriptional shifts rather than the underlying rewiring of interaction networks driving functional alterations. As a result, cell- and gene-focused analyses are often performed in parallel but remain decoupled: cell-state changes may be identified without linking to upstream regulators, and regulators may be prioritized without resolving their cell-state–specific effects. Consequently, there is a pressing need for a comprehensive computational framework that quantifies perturbation-driven GIN reorganization at both gene and cellular resolutions, incorporating prior biological knowledge while accounting for the heterogeneity characteristic of scRNA-seq data.

Here, we present scDNS, a computational platform that frames GIN remodeling as an information-theoretic divergence between network architectures across contrasting conditions. scDNS first assesses a gene's functional perturbation by quantifying dynamic alterations in its gene-gene dependency patterns, subsequently inferring the resultant effects on individual cells for detecting the perturbed cell populations. By coupling regulator-level rewiring with cell-resolved perturbation readouts, scDNS bridges upstream molecular drivers and downstream cellular consequences within a unified framework. In both simulated gene and cell perturbation tests, scDNS demonstrated superior sensitivity and robustness in detecting predefined perturbed genes and cells compared to current methods. Furthermore, application across diverse biological contexts—including DSP knockdown, concurrent mutations in TERT, TP53, and PDGFRA, IRF4pT95R-associated immunodeficiency, IFN-β stimulation in lupus-derived T cells, COVID-19 lung single-nucleus RNA-seq, and gemcitabine treatment in pancreatic cancer—consistently revealed previously overlooked perturbed regulators and delineated their cell-type-specific effects. We discovered the multifaceted effects of IRF4pT95R on B cells, natural killer (NK) and natural killer T (NKT) lineage development, upstream mediators of IFN-β-induced T-cell tolerance, alveolar cell susceptibilities in COVID-19 with therapeutic implications, and TIMM44 as a mitochondrial factor enhancing gemcitabine cytotoxicity. Collectively, these findings establish scDNS as a robust and versatile platform for quantifying functional perturbations within GINs. This capability significantly enhances our power to unravel the molecular basis of cellular responses, reveal previously latent regulatory processes, and ultimately accelerate the identification of therapeutic candidates across a wide array of biological systems.

Results
Overview of scDNS
Gene interaction networks, comprising multifaceted molecular interactions such as protein-DNA binding, protein-protein interactions, and signaling cascades, underpin the orchestration of transcriptional programs in cellular systems. Perturbations, whether experimental or pathological, can profoundly remodel these networks by modulating gene-gene dependencies, inducing transcriptional reprogramming, and reorganizing modular connectivity. Quantifying such reconfiguration provides a rigorous approach to elucidating perturbation effects on cellular phenotypes and revealing core regulatory principles. We introduce scDNS, a computational framework that conceptualizes network reconfiguration as the information divergence between two conditions. It offers a robust metric to evaluate gene-level functional disruptions and to characterize the cellular consequences of perturbations under contrasting conditions. The methodology proceeds through three sequential stages.

The first step reconstructs a context-specific GIN by integrating prior biological knowledge with single-cell transcriptomic data to capture condition-dependent interactions (Fig. 1a). We employed a Graph Attention Network (GAT)18, a specialized Graph Neural Network (GNN), to refine an initial reference network assembled from curated regulatory interactions (TRRUST19 and Dorothea20), protein-protein interactions (IntAct21 and STRING22) and signaling pathway databases (KEGG23 and SuperPathway24) (see Methods). This integrated prior network covers >90% of genes detected in the human single-cell datasets analyzed in our study, and ablation analyses confirm that incorporating these curated databases improves scDNS’s ability to infer biologically meaningful perturbation signals (Supplementary Fig. 1a&b). The model represents the network as a graph G=(V, E), where vertices V correspond to genes and edges E denote regulatory interactions. Vertex features are encoded using node2vec25 to encode network topology, while edge features are represented by embeddings that integrate three attributes of neighborhood relationships. To train the GAT, we adopted a semi-supervised contrastive learning strategy by leveraging high-confidence positive and negative edge labels to guide model optimization. Positive edges were defined as those involving housekeeping genes (e.g., consistently expressed across cell types) or cell-type-specific markers (e.g., derived from scRNA-seq cluster annotations), enabling the model to learn both universal and cell type-specific interactions. Negative edges were assigned to pairs with zero geometric mean expression in scRNA-seq data or randomly sampled gene pairs lacking supporting evidence in prior knowledge databases. GAT is optimized using a composite loss function comprising: (1) an edge classification loss, which evaluates attention scores against labeled positives to prioritize biologically relevant interactions, and (2) a graph information loss based on Deep Graph Infomax (DGI)26, which maximizes the mutual information between the encoder representations of the original graph and a corrupted graph generated by randomly shuffling vertices and edges. To ensure robustness, the model is trained ten times, with attention scores aggregated to rank edge reliability. Edges consistently receiving high attention scores (top 80% by ranking aggregated score) were retained as strong, context-specific interactions, while the lowest 20% were pruned to eliminate noisy or spurious connections. This GAT-based network inference is stable, reproducible, and scalable, yielding robust attention patterns and reliable convergence across marker-gene choices, high-confidence edge construction, stochastic training runs, and upstream clustering/annotation variability (Supplementary Note 1). It effectively distinguished true interactions from random ones, achieving a high area under the Receiver Operating Characteristic (ROC) curve of 0.99 (Supplementary Fig. 1c&d). This process yielded a refined GIN that balances universal regulatory interactions with cell-type-specific and condition-dependent edges, providing a robust foundation for downstream analyses.

The second step involves quantifying perturbation-induced alterations in gene–gene relationships within the refined GIN to identify regulatory genes across two contrasting conditions (Fig. 1b). To address data sparsity and heterogeneity inherent in single-cell data, we model the relationship between gene pairs using two complementary probability distributions: (1) a joint probability distribution that captures universal expression relationships prevalent in abundant cell populations, and (2) a conditional probability distribution that enhances sensitivity to regulatory shifts in sparse or transitional cell subpopulations27,28. Together, this dual-model framework encapsulates both global and context-specific gene-gene relationships (Supplementary Fig. 1e-h). The Jensen–Shannon Divergence (JSD) measures the information-theoretic distance of these distributions between the two conditions, enabling the detection of both broad and localized perturbations in gene expression patterns, and also providing insights into changes in cell population abundance or cell states. For each gene pair, the JSD value is weighted by the geometric mean of their expression levels to emphasize changes in transcriptionally active regions. The weighted JSD values are aggregated for each gene across its interacting partners to derive a functional perturbation score. To mitigate bias from heterogeneous regulatory module sizes, raw scores are transformed into degree-normalized Z-scores based on a null distribution derived from randomized networks.

The final step distributes each gene pair's JSD value to individual cells based on its pointwise divergence contribution (Fig. 1c). By summing the pointwise divergence contributions of all associated gene pairs for each gene, scDNS generates cell-specific perturbation scores that measure the impact of gene-level functional changes within each cell. All gene perturbation scores are then combined to compute a composite cellular effect score, which reflects the overall impact of perturbations on cell phenotypes, explaining changes in either cellular states or population proportions. In summary, scDNS quantifies functional perturbations at both gene and cellular levels by measuring network divergence across contrasting biological conditions, identifying key regulators and their responsive cell populations.

Systematic assessment of scDNS’s performance
We evaluated the performance of scDNS using both simulated and real scRNA-seq datasets. For stress testing, three simulated datasets were generated from iPSC-derived type 2 alveolar epithelial cells (iAT2) 29, H1 human embryonic stem cells (H1-hESCs), and PANC-1 cells30, respectively, with cells randomly assigned to control or in silico–perturbed groups. We perturbed each simulated dataset by modifying a subset of 100 genes under three biologically motivated regimes: (i) differential gene expression only (DGE), (ii) gene regulatory rewiring (DGR), and (iii) their combination (DGE+DGR); importantly, the simulation engine introduces these changes as structured, context-dependent perturbations to better mimic diverse, biologically realistic scenarios (Fig. 2a and Supplementary Note 2). To comprehensively characterize the problem difficulty and assess the power of our method to detect the perturbation signal, we generated ground truth case-control datasets by varying two key aspects. First, to further assess the effect of perturbation strength, we defined perturbed cells as hybrids of control cells and perturbed cells of variable relative weight (Methods). Decreasing the strength of the transcriptional difference between perturbed and unperturbed cells increased the difficulty of identifying the perturbed cell cluster. Second, perturbations were applied to varying fractions of cells, ranging from 10% to 100%. In these challenging scenarios, scDNS achieved accuracy exceeding 80% for either slight perturbation strength (>0.1) or a low fraction of perturbed cells (>10%) (Fig. 2b&c and Supplementary Fig. 2a). These results highlight the robustness of scDNS in detecting functionally perturbed genes, whether arising from altered expression levels or perturbations in gene-gene relationships.

We next assessed scDNS’s ability to recover perturbed cell subpopulations driven by diverse gene perturbations. In this validation setting, we artificially dysregulated 30 genes across three distinct subpopulations derived from the major cell clusters of the perturbed group. As an illustrative example, the perturbed subpopulations were not visually separated in the original UMAP embedding; however, incorporating the scDNS-derived perturbation scores (P-Score_cell) into the PCA embeddings successfully and clearly distinguished them (Fig. 2d & Supplementary Fig. 2b). The quantified stress tests across all three simulated datasets showed that scDNS maintains high accuracy even under weak and sparse perturbations (α = 0.25; 25% of cells perturbed), achieving AUC > 0.9 for perturbed-cell identification (Fig. 2e & Supplementary Fig. 2c).

To evaluate scDNS in real biological systems, we analyzed two scRNA-seq datasets: iAT2 cells with or without Desmoplakin (DSP) knockdown29, and H1-hESCs either wild-type or carrying triple mutations (TERTC228T + TP53G743A + PDGFRAD842V)31 (Fig. 2f). In both datasets, DSP (log₂FC = -0.27, false discovery rate (FDR) = 2.6 × 10⁻⁵, Wilcoxon test) and PDGFRA (log₂FC = 0.49, FDR < 6.4 × 10⁻384, Wilcoxon test) showed significant expression changes, therefore could be identified by standard differential expression analysis (e.g., Wilcoxon test), distributional differential expression methods (i.e., waddR), and Memento-de for evaluating differences in gene expression levels (Fig. 2g&h). However, the TP53 mutation did not significantly alter its expression (log₂FC = 0.05, FDR = 1, Wilcoxon test) but reconfigured its network of interacting genes (Supplementary Fig. 2d), rendering it undetectable by most conventional differential expression analyses. While Memento’s differential-variability module identified TP53, it failed to detect DSP and PDGFRA. Remarkably, scDNS was the only method capable of identifying all perturbed genes, including both differentially expressed and non-differentially expressed but network-perturbed genes. Furthermore, scDNS achieved superior accuracy compared with VIPER’s inferred protein activity analysis (99.8% vs. 70.5% for PDGFRA; 95.1% vs. 78.5% for DSP; 91.4% vs. 18.1% for TP53) (Supplementary Fig. 2e). Additional perturbed genes identified by scDNS were strongly enriched for biological processes associated with the phenotypes affected by these perturbed driver genes (Supplementary Fig. 2f). Despite its reliance on the reconstructed gene relationship network, scDNS exhibited remarkable robustness to the false-positive interactions, with a notable decline in accuracy only observed when the network contained more than 50% spurious interactions (Supplementary Fig. 2g). Moreover, scDNS demonstrated stable performance across genes with diverse degrees of network connectivity (Supplementary Fig. 2h&i).

We extended our validation by benchmarking scDNS against existing algorithms for detecting perturbed cells influenced by DSP knockdown or triple-gene mutations. For analyzing a single gene's perturbation effects, scDNS demonstrated superior discrimination between contrasting groups compared to metrics derived from inferred protein activity (VIPER) and raw gene expression levels (Fig. 2i). Furthermore, when assessing the overall perturbation effects across all genes, scDNS achieved the highest classification performance among all comparable methods (AUC = 0.99 for triple mutations dataset and 0.95 for DSP knockdown dataset) (Fig. 2j). Notably, Milo emerged as the second-best method, showing only slightly lower accuracy and specificity than scDNS (AUC = 0.92 for triple mutations dataset and 0.94 for DSP knockdown dataset). Collectively, these results robustly establish scDNS as a precise and powerful approach for identifying functionally perturbed genes and their affected cell populations across both simulated and real-world scRNA-seq datasets.

scDNS uncovers pleiotropic effects of IRF4 mutation on immune cells’ regulatory program and development
scDNS was applied to a single-cell RNA-seq profile of 18,640 B cells from healthy donors and individuals heterozygous for the IRF4pT95R mutation (Fig. 3a). The IRF4pT95R variant, which causes autosomal-dominant combined immunodeficiency (CID), preserves normal IRF4 mRNA and protein expression yet alters its canonical DNA-binding specificity, leading to a global rewiring of its regulatory network3. Such mutation-induced topological reconfigurations of gene–gene interactions exemplified the type of functional network perturbation that scDNS was designed to detect and quantify. Accordingly, scDNS pinpointed IRF4 as a significantly perturbed regulator (top 2%), characterized by significantly altered interaction relationships within its regulatory network (Supplementary Fig. 3a). This sensitivity contrasted with conventional single-cell analytical methods—such as Memento, waddR, and DEA—which either ranked IRF4 poorly or failed to detect its dysregulation, since the T95R mutation does not alter IRF4 mRNA abundance (log2 FC = -0.12, FDR = 1, Wilcoxon test) (Fig. 3b & Supplementary Fig. 3b). Even protein-activity inference methods like SCENIC only detected IRF4 dysregulation modestly (79% accuracy) (Supplementary Fig. 3c).

We observed that the perturbation score of IRF4 markedly elevated in plasmablast and plasma-cell subsets (Fig. 3c), consistent with the differentiation blockade reported in IRF4pT95R immunodeficiency3. However, Milo, MELD, and Augur were unable to identify this remarkable perturbation in plasma cell development except scDist (Supplementary Fig. 3d-f). To contextualize these perturbations, we reconstructed the B-cell differentiation trajectory (Fig. 3d) and mapped the scDNS-derived IRF4 perturbation scores across pseudotime (Fig. 3e). The scDNS method successfully captured functional perturbations even in early B cells, while metrics based on IRF4 gene expression and SCENIC activity were insensitive to these early changes, exhibiting weak changes (log2 FC < 0.1) (Fig. 3e). For quantitative validation, the IRF4 perturbation scores exhibited a strong correlation (R = 0.86, p = 2.2×10-16) with the aggregate expression of established IRF4-bound target genes identified by chip-seq data analysis, substantially exceeding the correlations obtained from gene expression (R = 0.49, p = 8.4×10-5) or SCENIC activity (R = 0.42, p = 9×10-4) (Supplementary Fig. 3g). Analysis of B cells before they commit to becoming plasmablasts and plasma cells revealed two distinct subpopulations (Fig. 3f), exhibiting differential expression of IRF4-bound target genes during B cell differentiation (Fig. 3g & Supplementary Fig. 3h). These two groups displayed an inverse relationship between IRF4 perturbation score and differentiation pseudotime, though they showed different perturbation levels (Fig. 3f). Specifically, one subpopulation underwent advanced differentiation despite minimal IRF4 perturbation (Fig. 3h). In contrast, the second group suffered from significantly stronger IRF4 dysfunction, leading to restricted differentiation. This suppression of differentiation was further supported by upregulation of key transcriptional repressors of B cell differentiation—including NOTCH232, BACH233, PAX534, and FOXO135 (Fig. 3i). This pattern ultimately indicates that B cells with mild perturbation proceed toward progressive differentiation, while those with strong perturbation remain arrested in immature states.

Beyond IRF4 itself, scDNS identified top 250 dysregulated genes, implicating pathways central to B-cell immunity and antiviral defense (Fig. 3j). Seventeen of these genes are annotated in Online Mendelian Inheritance in Man (OMIM) database36 as causative for immunodeficiency syndromes. Significantly, all but one (SEC61A1) showed a positive correlation with the IRF4 perturbation score and were components of a relevant protein-protein interaction network (Fig. 3k). These results suggest that interacting proteins undergo co-modulatory or synergistic functional changes, implying shared underlying pathogenic mechanisms. This is exemplified by IRF8, a known IRF4 functional partner during early B-cell development37, which showed coordinated functional perturbations with IRF4 in the progenitor stages (Fig. 3l).

Extending beyond the B-cell compartment, analysis of the whole peripheral blood mononuclear cells (PBMCs) dataset revealed that IRF4 dysfunction also perturbs the NK and NKT lineages, as evidenced by significantly elevated IRF4 perturbation scores (Fig. S3m&n). This finding contradicted the original clinical report, which found no alteration in NK-cell activation or abundance3. Our analysis further delineated the impact of IRF4 perturbation in the NK lineage, revealing a functional shift from short-lived effector cells (SLEC) toward memory-precursor cells (MPEC). This shift was transcriptionally characterized by the corresponding down-regulation of SLEC and up-regulation of MPEC gene signatures (Fig. 3o & Supplementary Fig. 3i). Collectively, these results demonstrate that scDNS is highly effective at pinpointing gene dysfunctions arising from pathogenic mutations and delineating their diverse impacts on cell differentiation across various immune lineages.

scDNS uncovers upstream regulators driving heterogeneous cellular responses to external stimuli in PBMCTs
To evaluate the value of scDNS in uncovering regulatory mechanisms of cell response to external stimuli, we applied it to a single-cell RNA-seq dataset of 13,999 PBMCs of 8 lupus patients treated or untreated with recombinant human IFN-β38 (Fig. 4a). IFN-β plays a central role in immune activation by interacting with the interferon alpha/beta receptor (IFNAR); however, its functional impacts across diverse cell types and developmental states remain incompletely understood. For perturbed genes identification, scDNS detected 2784 genes significantly responsive to IFN-β (adjusted p < 0.05) (Fig. 4b). Notably, IFNAR1—whose expression change was minimal (log2 FC = -0.13, Supplementary Fig. 4a)—was ranked much higher by scDNS (66th, 99.4%) than by conventional differential expression approaches. The best-performing DEA was memento’s differential-variability analysis (accuracy 88%), while differential analysis of VIPER-inferred protein activity reached 73% accuracy. These results underscore scDNS’s superior sensitivity for identifying direct treatment targets (Fig. 4c). Gene set enrichment analysis (GSEA) further showed that scDNS-identified genes were strongly enriched in interferon-related pathways, including “cellular response to type I interferon” and “type I interferon signaling pathway” with stronger significance levels compared to alternative methods (Fig. 4d). This improvement arose because scDNS not only captured canonical transcriptional responses to IFN-β but also revealed upstream regulators of the signaling cascade that were missed by conventional gene expression analysis (Fig. 4e). For example, JAK1 and WTAP—despite showing negligible transcriptional changes under IFN-β stimulation (log₂FC = 0.01, Wilcoxon rank = 7551; log₂FC = 0.06, Wilcoxon rank = 8011) (Supplementary Fig. 4b) —were prioritized by scDNS (JAK1: scDNS rank=2263, FDR = 0.02; WTAP: scDNS rank=1692, FDR = 0.007). This observation suggested that IFN-β modulates their activity via non-transcriptional mechanisms. Indeed, IFN-β stimulates JAK1 activation through receptor-mediated autophosphorylation within the IFNAR–JAK–STAT signaling axis39. For WTAP, a recent study demonstrated that IFN-β induces its phase separation into liquid-like condensates, thereby modulating its role in m6A RNA methylation40. Consistently, we observed significant enrichment of scDNS-identified genes in related pathways governed by JAK1 and WTAP activity (Supplementary Fig. 4c).

We next examined cell type–specific effects of IFN-β treatment across 8 annotated clusters. By mapping scDNS-derived integrated sc-Zscores onto the UMAP embedding, we observed heterogeneous responses, indicating distinct sensitivities of individual cell states to IFN-β stimulation (Fig. 4f). CD4+ T cells and monocytes exhibited relatively stronger responses than other cell types, with high variability (Fig. 4g). In monocytes, scDNS revealed a significant and stronger response to IFN-β in CD14+ (classical) subpopulations than CD16+ (non-classical intermediate) monocytes, aligning with prior studies41, whereas MELED and Milo did not detect this difference (Supplementary Fig. 4d-h). Within CD4+ T cells, memory and naïve subsets were the most responsive, while activated states showed less perturbation (Fig. 4h). Beyond shared responsive genes, scDNS identified distinct mechanisms across three subsets, with memory cells upregulating genes associated with antigen processing and presentation, while naïve cells responded primarily through T cell proliferation (Fig. 4i). Importantly, scDNS revealed that activated T cells responded through the induction of cell tolerance, a key mechanism in lupus treatment. The majority of these responsive genes were not detected by standard differential expression analysis. Collectively, these findings demonstrate that scDNS is a valuable tool for discovering key upstream regulators of cascade signaling pathways and heterogeneous cellular perturbations in response to external treatments.

scDNS identifies cellular vulnerabilities and therapeutic targets in COVID-19
To further assess scDNS’s accuracy in the presence of cohort effects (i.e., individual-to-individual variation), we applied it to lung single-nucleus RNA-sequencing (snRNA-seq) data from 19 COVID-19 patients who succumbed to severe infection and 7 healthy donors42(Fig. 5a). scDNS assigned higher ranks to 1738 host genes known to interact with SARS-CoV-2 RNA and viral proteins, with their perturbation scores exceeding those of other genes (Fig. 5b). Based on the aggregated perturbation score of these host genes, we found that SARS-CoV-2 exerted a significantly high perturbation effect on alveolar type 1 (AT1) and alveolar type 2 (AT2) cells (Fig. 5c & Supplementary Fig. 5a). This finding aligns with the well-established role of alveolar epithelial cells as the primary host cells for SARS-CoV-243. Notably, scDNS-prioritized genes showed significant enrichment in phosphoproteomic and proteomic signatures detected from SARS-CoV-2–treated Vero E6 cell lines in vitro44, with enrichment scores rising with treatment duration (Supplementary Fig. 5b). These findings demonstrate scDNS’s ability to identify genes with functional roles at the translational and post-translational levels.

At a finer resolution, scDNS revealed heterogeneous effects of SARS-CoV-2 on AT2 cells, characterized by a functional gradient ranging from low to high perturbation scores (Fig. 5d). Cells with high perturbation scores exhibited significant upregulation of 616 genes that can be considered as a transcriptional signature of severely infected cells (Fig. 5e). These signatures were enriched in pathways critical for SARS-CoV-2 survival and replication, including antigen presentation, protein translation, viral replication, and life cycle regulation (Supplementary Fig. 5c). Among these, we identified elevated expression of secretory leukoprotease inhibitor (SLPI) in cells with high perturbation scores (Fig. 5e). Spatial transcriptomics45 revealed that AT2 cells with these signatures were primarily located in bronchial-proximal regions (Fig. 5f), highlighting this niche as a key site of SARS-CoV-2 infection. This spatial-functional relationship was underscored by increased SLPI expression, an anti-protease predominantly produced in the airway46 (Fig. 5g). We further applied scDNS to immune cells, revealing pronounced heterogeneity in response to SARS-CoV-2 infection (Supplementary Fig. 5d). Signature genes from cells with high perturbation scores were mapped onto bulk PBMCs RNA-seq data from 274 COVID-19 patients47, showing that signatures of CD4⁺ T cells, CD8⁺ T cells, cycling NK/T cells, and dendritic cells increased markedly with disease severity—more prominently than genes identified by direct case–control comparisons (Supplementary Fig. 5d-h). Functional enrichment analysis further demonstrated that these upregulated signatures were significantly enriched in inflammatory pathways, COVID-19 infection–related signaling cascades, and previously reported severity-associated gene program (Supplementary Fig. 5i-k, Supplementary Data 1 and Supplementary Note 3). These signatures increased markedly with disease severity, establishing a strong link between transcriptional dysfunction of these signatures and COVID-19 progression.

To assess the therapeutic potential of scDNS-prioritized genes, we evaluated the antiviral activity of drugs and compounds targeting them. For this analysis, we utilized antiviral efficacy profiles of ~10,000 approved drugs and compounds from the NCATS COVID-19 portal and linked each compound to its targets using DrugBank annotations. Across NCATS assays, only 40% of drugs showed antiviral activity, and this proportion increased modestly to 45% when we restricted to virus-interacting host proteins (Fig. 5h). In contrast, prioritizing targets by scDNS based on high perturbation scores with increasingly stringent thresholds substantially improved success rates, reaching 56% among the top 50 genes and 62% when further limited to those host genes directly interacting with SARS-CoV-2. Furthermore, projection of scDNS-nominated genes onto drug–target interaction networks identified 76 genes linked to 257 candidate compounds, all of which exhibited anti–SARS-CoV-2 activity in at least two NCATS assays (Fig. 5i). These compounds included modulators of viral replication and host immune signaling pathways, providing a diverse repertoire of translationally promising therapeutic options. Collectively, these findings highlight scDNS as an effective framework for prioritizing therapeutic candidates in cohort dataset analysis, while reducing the risk of failure.

scDNS identifies TIMM44 as a regulator sensitizing pancreatic cancer cells to gemcitabine
To demonstrate the ability of scDNS to facilitate novel biological discovery, we applied it to a scRNA-seq dataset of 1,749 in vitro PANC-1 pancreatic cancer cells subjected to gemcitabine treatment or vehicle control30 (Fig. 6a). scDNS identified 2669 genes implicated in gemcitabine’s mechanism of action (MoA). The top-enriched pathway was the DNA damage response, consistent with gemcitabine’s well-established function in inhibiting DNA replication and inducing DNA damage through incorporation into DNA and depletion of deoxyribonucleotide pools (Fig. 6b). Notably, scDNS successfully prioritized the known gemcitabine targets RRM1 and RRM2 as the 44th and 27th most perturbed genes, respectively. In addition, 448 highly ranked genes were linked to cell-cycle regulation and DNA repair, 23 of which have been experimentally confirmed to contribute to gemcitabine’s MoA (Fig. 6c). For example, CHK1, a central mediator of gemcitabine-induced DNA damage signaling48, showed minimal differential expression (log2 FC = 0.14, FDR = 1, Wilcoxon test) yet a strong perturbation score (Z-score = 6.3, FDR = 2.1×10-8), highlighting scDNS’s ability to capture functional disruptions beyond expression changes. Beyond canonical DNA damage pathways, scDNS also revealed non-canonical effects of gemcitabine on mitochondrial function (Fig. 6d). Among these, SIRT1 (DEA: log2 FC = 0.18, FDR = 0.06, Wilcoxon test; scDNS: Z-score = 6.3, FDR = 2.5×10-14), an NAD⁺-dependent deacetylase that promotes mitochondrial biogenesis and redox homeostasis, was markedly perturbed and has been reported to mediate gemcitabine’s MoA49.

Based on cell-level perturbation scores, we identified six cellular subpopulations exhibiting distinct responses to gemcitabine, with cluster 3 showing the strongest perturbation (Fig. 6e). Projection of each subpopulation’s transcriptional signature onto TCGA pancreatic adenocarcinoma (PAAD) patient survival data revealed significantly different survival outcomes, where the signatures of cluster 3 are associated with the poorest prognosis (hazard ratio = 1.76, p = 3.3 × 10⁻5, two-sided Wald test) (Fig. 6f). In contrast, clusters defined solely by gene expression showed weaker prognostic stratification (Supplementary Fig. 6a&b). Cluster 3 was characterized by metastatic potential, impaired apoptotic regulation, and enhanced invasiveness (Fig. 6g). Genes with the highest perturbation scores in cluster 3 were enriched for cell-cycle and mitochondrial functions (Fig. 6h). Among them, SDHA (succinate dehydrogenase complex flavoprotein subunit A) was the most differentially perturbed gene (Fig. 6i). SDHA deficiency leads to succinate accumulation, which has been implicated in drug resistance through ubiquitin-cullin pathway dysregulation in several cancers, including acute myeloid leukemia50. The second-most perturbed gene, TIMM44, a core component of the presequence translocase-associated motor (PAM) complex of the mitochondrial inner membrane, exhibited no significant change in expression (log2 FC = 0.33, FDR = 1, Wilcoxon test) but a marked perturbation in network activity (Z-score = 9.0, FDR = 1.1×10-15) (Fig. 6i). Western blot analyses confirmed stable TIMM44 expression before and after gemcitabine treatment (Supplementary Fig. 6c). Complementary proteomics/phosphoproteomics further showed stable total TIMM44 protein but a significant increase in S180 phosphorylation after treatment (FC = 2.6, p = 0.006, two-sided t-test) (Supplementary Fig. 6d). S180 maps to a functionally important domain involved in TIMM44’s association with the translocon and recruitment of the mtHsp7051, and is supported as a bona fide regulatory site by recurrent detection in PhosphoSitePlus and strong evolutionary conservation (Supplementary Fig. 6e–g). These orthogonal data are consistent with a potential activity or interaction-state shift of TIMM44, providing supportive evidence for its prioritization by scDNS even in the absence of differential expression. Independent bulk RNA-seq datasets from TCGA and GTEx showed abnormally elevated TIMM44 expression in tumor versus normal tissues, indicating its potential as a therapeutic target (Supplementary Fig. 6h). Similarly, an independent single-cell dataset52 confirmed TIMM44 upregulation specifically in tumor cells (Supplementary Fig. 6i), supporting its potential as a tumor-selective target with low off-target-risk.

To experimentally validate TIMM44-mediated sensitization to gemcitabine, we performed in vitro combination assays using gemcitabine and the mitochondrial import inhibitor MitoBloCK-10 (MB-10). MB-10 potentiated the cytotoxic effect of gemcitabine (Fig. 6j & Supplementary Fig. 6j). At low gemcitabine concentrations (≤8 μM), the combination of MB-10 and gemcitabine exhibited a synergistic effect with an average Chou–Talalay combination index (CI) of 0.44. Transcriptomic profiling revealed that gemcitabine induced major transcriptional changes, further amplified in the gemcitabine + MB-10 combination (Fig. 6k). Differential expression analysis showed that downregulated genes were enriched for cell-cycle and DNA repair pathways. Notably, several of these genes—including SIRT149, DUSP153, EP30054, and MCL155—have been previously implicated in gemcitabine resistance Fig. 6l&m). Collectively, these results demonstrate that scDNS effectively delineates heterogeneous drug responses, identifies non-canonical perturbation mechanisms, and uncovers novel sensitivity-associated genes such as TIMM44, highlighting its potential as a therapeutic target for improving gemcitabine efficacy.

Methods
Context-adaptive GINs reconstruction
A semi-supervised analysis framework based on graph attention networks was designed to identify context-adaptive regulatory networks by integrating single-cell transcriptomic data with a comprehensive prior knowledge network. The overall analysis pipeline was comprised of the two steps:
(1) Prior knowledge network construction: A comprehensive biological networks encompassing gene regulatory interactions (GRIs), protein-protein interactions (PPIs) and signaling pathways was firstly constructed. The GRIs were sourced from the TRRUST database19 and the R Dorothea package (v1.7.4)20. PPIs were obtained from the IntAct21, and STRING databases22. Signaling pathways were collected from the KEGG23 and SuperPathway24 databases. In total, 332,740 unique interactions involving 21,006 genes were collected.
(2) Pre-set the high-confidence positive and negative edges: Standard Seurat preprocessing was applied to the single-cell RNA sequencing data to identify major cell types. Positive edges were defined as those connecting highly expressed and biologically relevant gene pairs within each cell type. Specifically, the COSG57 method was used to identify the top 100 marker genes for each cell type. In addition, housekeeping genes with an expression frequency greater than 75% and average expression above 1 were incorporated. These genes, along with their interactions in the prior knowledge network, were classified as positive interactions for each cell type. Negative edges were defined as (i) gene pairs exhibiting a geometric mean expression of zero across all cell types, and (ii) those randomly sampled pairs lacking prior biological evidence of interaction.

Node and edge encoding
Node embeddings for the prior knowledge network were obtained by applying the Node2Vec algorithm to map each vertex 𝑖 to a 512-dimensional vector 𝐡⃗⃗⃗_𝑖. This process captured both local neighborhood structure and global topological properties. For a candidate edge between genes 𝐺_𝑖 and 𝐺_𝑗 in cell population 𝐶_𝑘, the edge feature vector was constructed by integrating: (i) cell-specific co-expression 𝑓_exp(𝐺_𝑖,𝐺_𝑗,𝐶_𝑘), defined as the geometric mean expression of gene pair within 𝐶_𝑘; (ii) global co-expression statistics 𝑓_MI(𝐺_1,𝐺_2) and 𝑓_DREMI(𝐺_1,𝐺_2), computed as mutual information (MI) and density-resampled mutual information (DREMI) between genes across the whole dataset, log-transformed and sigmoid-normalized; and (iii) contextual features (𝑓_context(𝐺_1,𝐺_2)), a one-hot encoding of biological context (e.g., cell identity and treatment). The final edge feature vector is therefore 𝐡⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗_𝐺_𝑖,𝐺_𝑗,𝐶_𝑘 =[𝑓_exp,𝑓_DREMI,𝑓_MI,𝑓_context] (1).

Graph attention mechanism for context-adaptive regulatory network prediction
A GATv2-based graph attention encoder was employed to couple with a joint contrastive learning module. This module incorporates graph-level representation learning using Deep Graph Infomax (DGI) and an attention score-based edge supervision mechanism.

GAT encoder module: For a cell-type-specific gene network 𝐺=(𝑉,𝐸), where 𝑉 is the node set and 𝐸 is the edge set representing potential interactions 𝑒_𝑖,𝑗,𝐶_𝑘 between genes 𝑔_𝑖 and 𝑔_𝑗 in cell type 𝐶_𝑘, these edges were used employed for joint encoding. Given the node input feature 𝐡⃗⃗⃗_𝑖 and edge feature 𝐡⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗_𝐺_𝑖,𝐺_𝑗,𝐶_𝑘, the features are linearly transformed using the learnable weight matrices 𝑊_𝑣 ∈ℝ^{|𝑉|×𝐹′} and 𝑊_𝑒 ∈ℝ^{|𝐸_all|×𝐹′}, where 𝐹′ is the embedding dimension, and |𝐸_all| is the total number of edges across all cell type networks. Then, the attention score for the message propagation from node 𝑗 to node 𝑖 in cell type 𝐶_𝑘 was calculated as:
𝛼_𝑖,𝑗,𝐶_𝑘 =Tanh(𝑎⃑^𝑇[𝑊_𝑣𝐡⃗⃗⃗_𝑖+𝑊_𝑣𝐡⃗⃗⃗_𝑗+𝑊_𝑒𝐡⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗_𝑖,𝑗,𝐶_𝑘]) (2)
where 𝑎⃑ is a learnable attention vector, and Tanh is the hyperbolic tangent activation function. This score reflects the importance of the neighbor 𝑗 and its connecting edge to the representation of node 𝑖. A Softmax normalization is applied over all incoming edges to node 𝑖 to obtain the normalized attention weight
𝛼′_𝑖,𝑗,𝐶_𝑘 =Softmax_𝑗∈𝑁_𝑖 (𝛼_𝑖,𝑗,𝐶_𝑘) (3).
This normalized score is used for weighted aggregation during message passing and serves as a proxy for the edge existence probability in the next supervised task.

Joint contrastive learning module: The model's GAT encoder was jointly optimized using both unsupervised contrastive learning and supervised edge classification. Deep Graph Infomax (DGI) is utilizedto learn structure-aware node embeddings 𝑍=𝑓_GATv2(𝑉,𝐸). This is achieved by maximizing the mutual information between the node representations 𝑍 and the graph-level summary vector 𝑠=sigmoild(mean(𝑍)), while minimizing the mutual information between the corrupted graph's node representations 𝑍̃ (generated by randomly shuffling input node features) and the summary 𝑠. The DGI loss is defined as:
ℒ_DGI =−𝐸_𝑍,𝑍̃ [log𝜎(𝑧^𝑇𝑠) + log(1 − 𝜎(𝑧̃^𝑇𝑠))] (4).
This loss ensures the model can distinguish the true graph structure from structural noise
Additionally, the normalized attention weights serve as the prediction 𝑦̂_𝑖,𝑗,𝐶_𝑘=𝛼′_𝑖,𝑗,𝐶_𝑘 for the ground-truth labeled edges 𝑦 (1 for positive interactions and 0 for negative ones). Focal Loss is employed to mitigate the severe positive/negative sample imbalance:
ℒ_edge classification = −𝛼(1-ŷ)^𝛾log(𝑦̂) − (1−𝛼)(1-y)𝑦̂^𝛾log(1−𝑦̂) (5)
Hyperparameters α and γ were set to 0.75 and 1.5, respectively, ensuring that the model focused on hard-to-classify samples. The final overall training objective is the weighted sum of the supervised and self-supervised losses:
ℒ_total =ℒ_edge classification + 𝜆∙ℒ_DGI (6)
The weight coefficient 𝜆 was set to 0.1 to balance the two objectives.

Training and edge prediction: Prior networks for each cell type across different conditions were integrated into a multiplex graph as the GAT input, allowing each gene pair to have multiple edges representing diverse cellular contexts. This design enables the model to simultaneously learn interaction patterns across these contexts. To enhance the model robustness and reproducibility, we independently trained the model across 10 random seeds. For each cell type, the attention scores for each edge across the 10 training runs were aggregated, and their Rank Product (RP) was then computed. The final ranking of gene interactions for each cell type was derived by independently sorting the edges based on their RP. To ensure high-confidence gene interactions, the minimum rank across all cell types was employed as a final evaluation criterion. Gene interactions that consistently ranked in the lowest 20% across all cell types were considered to have a low likelihood of occurring in the single-cell data system.

Estimate the perturbation scores for genes and cells
scDNS characterized the dynamic relationships of gene pairs in the network (G_i↔G_j) using joint and conditional gene expression probability densities. Each method had distinct advantages in capturing interactions within both high-density and low-density cell populations, allowing for a comprehensive analysis. The distribution densities were estimated based on a k-nearest neighbor (KNN) approach on a two-dimensional grid, following the widely used procedure of van Dijk et al28. A small random expression value was incorporated to improve robustness against extreme distributions, such as the zero-inflated distribution frequently observed in single-cell data.

Jensen-Shannon divergence (JSD) was utilized to quantify changes in probability densities across conditions. It provided valuable insights into both local and global information divergences between two probability distributions, providing a robust metric for assessing gene interaction dysregulation. To emphasize biologically relevant interactions, point-wise JSD was weighted by gene expression levels. The statistical significance of JSD was determined using a null distribution generated from 2×10^4 JSD values derived from randomized gene pairs not included in known interaction networks, with additional randomization of gene expression within each condition to account for background variability.

Finally, gene prioritization was performed at both global and single-cell levels. The significance of JSDs associated with a given gene was aggregated while accounting for network degree to avoid bias. Z-scores from joint and conditional probability densities were then combined to obtain an overall gene-level Z-score. Finally, single-cell Z-scores (sc-Zscores) were computed to quantify each gene’s effector cells based on weighted point-wise JSD. A detailed description of the algorithm was provided in the Supplemental information.

Simulation setup for benchmarking datasets
Three simulated datasets were generated from iAT2, H1-hESCs, and PANC-1 cells, with cells randomly assigned to control or in silico–perturbed groups. Perturbations were applied to a subset of 100 target genes under three biologically motivated regimes: (i) expression-only changes (DGE), (ii) rewiring-only changes affecting regulatory relationships or expression patterns (DGR), and (iii) coupled expression and rewiring changes (DGE+DGR). For rewiring changes, the relationships between each target gene and its partner genes were systematically modified to simulate context-dependent regulatory changes, including activation/inactivation and switches between positive and negative regulation. Only the top 60% of partner genes, ranked by their average Spearman correlation with other partners, were selected for rewiring. Perturbations were constrained within local neighborhoods and functional modules to preserve structured and biologically plausible alterations, enabling benchmarking of scDNS under realistic regulatory complexity. Expression-level perturbations were simulated by multiplying the expression of each target gene in perturbed cells by a factor of two. Combined perturbations were generated by applying both expression-level and rewiring changes simultaneously.

To systematically assess the effect of perturbation strength, perturbed cells were defined as hybrids between control cells and fully perturbed cells. We introduced a perturbation strength parameter, α (ranging from 0 to 1), to control the similarity between perturbed and control expression profiles. The resulting expression profile of each perturbed cell was computed as:
𝐸_𝛼 =(1−𝛼)·𝐸_before perturbation +𝛼·𝐸_after perturbation (7)
where 𝐸_before perturbation and 𝐸_after perturbation represent the expression profiles of the control and perturbed cells, respectively. Perturbation strength α was varied to create a gradient of transcriptional differences, simulating conditions ranging from subtle to strong perturbations. We additionally varied the fraction of cells subjected to perturbation from 10% to 100% to assess scDNS’s sensitivity across different levels of heterogeneity.

To further evaluate scDNS’s ability to detect subpopulation-specific perturbations, we first performed clustering of the single-cell data to define independent cell subpopulations. For three selected subpopulations, we then applied DGR, DGE, and combined perturbations to 30 high-expression marker genes per subpopulation, simulating biologically relevant perturbation scenarios across distinct cell clusters.

Comparison of perturbation detection methods
Gene-level perturbation analysis was performed using waddR (v1.15.0)10, memento (v0.1.1)11 in two modes (memento_de and memento_dv), and Seurat’s FindMarkers function (v4.3.0)58, including the following statistical frameworks: Wilcoxon rank-sum test, bimod likelihood-ratio test, ROC-based marker identification, Student’s t-test, negative binomial and Poisson generalized linear models, logistic regression, and MAST hurdle model. Cell-level perturbation detection was conducted using MELD (v1.0.2)16, Milo (v2.5.1)14, Augur (v1.0.3)17, and scDist (v1.1.4)15.

scRNA-seq data pre-processing
The Seurat R package (v4.3.0)58 was utilized to preprocess the scRNA-seq data. Following standard quality control procedures in single-cell data analysis, low-expression genes were filtered out and then normalized using the NormalizeData function with default parameters. Highly variable genes were identified using the FindVariableFeatures function with the default “vst” method. Subsequently, the principal component analysis (PCA) was applied to the scaled gene expression of highly variable genes. The datasets exhibiting batch-related differences were integrated using the Harmony R package (v0.1.0)59 with default parameters. The first 30 PCs were used to construct a shared nearest neighbor graph for cells and to identify cell clusters using the FindNeighbors and FindClusters functions, respectively. The clustering results were visualized using UMAP by the RunUMAP function. Cell annotation was performed using marker genes derived either from the original dataset or from a curated reference. Marker genes for each cluster were identified with the FindAllMarkers function with the “wilcox” method.

Single-cell protein activity inference
Single-cell transcriptomic data were first clustered to define coherent cell subpopulations. For each cluster, pseudo-bulk gene expression profiles were generated by averaging the expression of constituent single cells. Gene regulatory networks were reconstructed from these pseudo-bulk profiles using ARACNe60 for identifying the downstream targets (regulons) of each protein. Single-cell protein activity was subsequently inferred using the VIPER R package61. Proteins with small regulon size (<20 target genes) were excluded. Differences in protein activity between groups were assessed using the Wilcoxon rank-sum test.

Analysis of the expression levels of IRF4 target genes
The IRF4T95R-specific regulatory network was constructed using differential ChIP-seq peaks and expressed genes between IRF4T95R and wild-type samples (data available in reference3). The target genes of IRF4T95R were defined as those exhibiting an absolute log fold change greater than 1 in both gene expression levels and normalized peak counts. We classified the target genes into positive and negative regulators based on their upregulation and downregulation in expression levels. The regulatory network of IRF4pT95R was then employed to infer its specific abnormal activity in each cell through enriched regulon analysis implemented by the viper function in the viper R package (v1.28.0)61. IRF4 transcriptional activity was inferred using the pySCENIC Python package (v0.12.1)13. Pseudotime and stemness analyses of IRF4 p.T95R B cells were performed using Monocle3 (v1.3.1)62 and CytoTRACE (v 0.3.3)63 with default parameters.

Survival analysis
Cluster-specific marker genes (|log FC| > 0.25, adjusted p < 0.01) were identified from gemcitabine-treated single-cell data using the FindAllMarkers function in R Seurat. To evaluate clinical relevance, we retrieved TCGA pancreatic adenocarcinoma (TCGA-PAAD) raw count matrices and associated clinical metadata through the R TCGAbiolinks package (v2.35.3)64. Each cluster’s highly expressed genes were treated as a transcriptional signature and projected onto bulk tumor profiles by computing sample-level signature enrichment/activity scores using the VIPER algorithm (R viper package (v1.32.0)61), producing one quantitative score per signature per sample. Associations between signature scores and patient outcome were assessed using univariate Cox proportional-hazards regression with the R survival package (v3.5-3)65; signature score was modelled as a continuous covariate and results are reported as hazard ratios (HRs) with 95% confidence intervals.

Cell culture and drug treatment
Human pancreatic adenocarcinoma cell line, PANC-1 cells (CL-0184, Wuhan Pricella Biotechnology Co., Ltd.), was cultured in DMETM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin at 37°C in a humidified incubator with 5% CO2. Cells were seeded in 96-well plates (3 × 10³ cells/well) and allowed to adhere overnight. Then, cells were treated with MitoBlock-10 (HY-115467, MedChemExpress) and gemcitabine (HY-17026, MedChemExpress), separately or in combination, for 72 h. Cell viability was assessed using the CCK-8 assay (C0037, Beyotime Biotechnology. Briefly, 10 μL of CCK-8 reagent was added to each well and incubated for 2 h at 37°C. Absorbance was measured at 450 nm using a microplate reader. Background absorbance from blank wells was subtracted. Experiments were performed with three independent biological replicates.

Drug interactions were evaluated by calculating the combination index (CI) using the Chou-Talalay method (https://www.combosyn.com/). CI values were calculated using Compusyn software based on dose–response curves generated from multiple drug concentrations. CI < 1 indicates synergy, CI = 1 indicates an additive effect, and CI > 1 indicates antagonism.

For Western blot analysis, cells were treated with gemcitabine at indicated concentrations. Cells were lysed in RIPA buffer supplemented with protease inhibitors. Protein concentration was determined using a BCA protein assay. Equal amounts of protein (20 μg per lane) were separated by 10% SDS-PAGE and transferred onto PVDF membranes. Membranes were blocked with 5% non-fat milk for 1 h at room temperature and incubated overnight at 4°C with anti-TIMM44 antibody (A21141, Abclonal; 1:1,000 dilution). After washing, membranes were incubated with HRP-conjugated secondary antibody (1:5,000 dilution) for 1 h at room temperature. Signals were visualized using enhanced chemiluminescence (ECL). TIMM44 protein abundance was detected using anti-TIMM44 antibody (A21141, Abclonal) and HRP-conjugated secondary antibody, with visualization by ECL. Protein levels were normalized to GAPDH (A10868, Abclonal).

RNA sequencing for pancreatic cancer cells
For transcriptomic analysis, cells were treated with MB10 (0.64 μM), gemcitabine (0.50 μM), their combination, or DMSO as a vehicle control. Total RNA was extracted from PANC-1 cells using the Total RNA Extractor (Trizol) (B511311, Sangon Biotech) following the manufacturer's protocol. Libraries were prepared using the Illumina TruSeq RNA Sample Preparation Kit, following standard protocols. The constructed libraries were sequenced on Ran Illumina NovaSeq 6000 platform. Raw sequencing data were processed using FastQC for quality control and trimmed with Fastp to remove low-quality reads and adaptors. The resulting clean reads were aligned to the human genome (GRCh38) using STAR aligner. Gene expression profiles were quantified with feature Counts. Differential expression analysis was performed using DESeq2 (v1.34.0)8 with default parameters, and genes with absolute log fold change > 1 and adjusted p < 0.01 were defined as differentially expressed.

Phosphoproteomic profiling by LC–MS/MS
PANC cells were seeded 24 h before treatment and, at 60–70% confluency, treated with equal volumes of DMSO (control), 1 μM gemcitabine for 36 h. Cells were washed with ice-cold PBS, harvested by scraping and snap-frozen in liquid nitrogen. Frozen cell pellets were lysed in RIPA buffer (P0013B, Beyotime) supplemented with 1 mM PMSF (ST507, Beyotime) and 1× phosphatase inhibitor cocktail A (P1082, Beyotime). Lysates were sonicated on ice and centrifuged at 15,793 × g for 10 min at 4 °C. Supernatants were collected as total protein extracts, and protein concentration was determined using a BCA assay (23225, Thermo Scientific). Equal amounts of protein were processed using the SP3 (single-pot solid-phase-enhanced sample preparation) method, including acetonitrile-induced binding to magnetic beads, ethanol and acetonitrile washes, reduction with dithiothreitol, alkylation with chloroacetamide, and digestion with Trypsin/LysC. Phosphopeptides were enriched using an IMAC phosphopeptide enrichment kit (A32992, Thermo Scientific) according to the manufacturer’s instructions. Enriched peptides were desalted using SOLA™ SPE 96-well plates and dried in a vacuum concentrator prior to LC–MS/MS analysis. Peptides were separated on a Vanquish Neo UHPLC system coupled to an Orbitrap Astral mass spectrometer (Thermo Fisher Scientific). Samples were loaded onto a trap column and separated on a PepMap™ Neo analytical column (15 cm × 150 μm, 2 μm). Mobile phase A consisted of 0.1% formic acid in water and mobile phase B consisted of 0.1% formic acid in 80% acetonitrile.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Current strategies for scRNA-seq analysis focus on differential expression analysis methods coupled with pathway enrichment analysis. However, genes with unchanged expression yet profoundly modified regulatory associations remain hard to detect. To address this knowledge gap, in this manuscript, Huang et al. introduces scDNS, a computational framework designed to quantify \"network rewiring\" in single-cell data. The authors combine Graph Attention Networks (GAT) with Jensen-Shannon Divergence (JSD) to identify regulators whose functional impact stems from altered interaction patterns rather than changes in transcript abundance.\nThe concept is quite novel, which aims to identify key genes playing important roles during a certain altered condition without showing the commonly assumed expression changes. The study is well-structured, and provided benchmark against existing tools (VIPER, Milo, MELD, etc.). The tool also demonstrate broad utility across diverse research settings, including primary immunodeficiency (IRF4 mutation), viral infection (IFN and COVID-19), and oncology (Gemcitabine resistance in PDAC). The experimental validation of TIMM44 as a therapeutic target is also a highlight demonstrating the power and robustness of this computational tool. However, several major technical and computational points require clarification before publication.\nMajor Comments\n1. Simulated perturbations\nIn methods “Simulation setup for benchmarking datasets”, the scDNS tool introduces perturbations in three ways (shuffle, multiply, mix). While this is ok, this form of perturbation results in dramatic statistical divergence that may not accurately reflect the structured, non-random rewiring of regulatory networks observed in actual pathological conditions. Can the authors add more ways of perturbation to mimic the diverse scenarios of gene changes in biological settings? In addition, only one single cell dataset is used, it would beneficial if the authors could include more scRNA-seq datasets.\n2. Algorithmic overfit\nIn figure2e, scDNS achieve an AUC value of near 1.0 in simulated data while other softwares are around 0.5, which is quite poor. This result alerts that there is an overfitting problem here. The core point of scDNS is to quantify the differences in probability distributions using JSD divergence. The perturbation methods used, as also discussed above, directly maximizes the JSD value, creating an inherent \"home-court advantage\" for scDNS compared with other published methods.\n3. Benchmark\nFollowing the previous question, in figure 2e-g benchmarking, Milo and scDist showed near-zero sensitivity in identifying perturbed subpopulations in specific DGR (Differential Gene Relationship) scenarios, while scDNS perform very good. The baseline methods chosen for comparison (e.g., Milo, Augur, scDist, VIPER) were not originally engineered to detect statistical divergence caused by random shuffling.\nMilo and scDist are optimized for detecting shifts in neighborhood abundance or inter-condition distances. VIPER relies on established regulon structures to infer protein activity. Evaluated against a perturbation scheme that mirrors the specific mathematical logic of scDNS, their performance decline is expected.\n4. Potential bias of prior knowledge networks\nThe scDNS framework relies on integrating established databases (TRUST, STRING, KEGG, etc.) to construct the initial reference network. There is a concern regarding \"research bias\": the model might be predisposed to find perturbations in well-studied genes, while \"orphan genes\" with fewer known interactions remain undetectable.\nHow about the authors perform an experiment by randomly removing varying percentages (e.g., 20%-50%) of prior edges to assess the robustness of identifying key regulators (e.g., TP53)?\n5. DGRs validation\nThe scDNS can detect both genes with altered expression levels (DGE) and, interestingly, genes with unchanged expression levels but altered regulatory networks (DGR). In figure4, JAK1 and WTAP is identified as DGRs, while in figure6, TIMM44 is also considered as a DGR. These gene show a little mRNA change but high network perturbation. Some experimental validation in a certain setting will significantly strengthen the analytical power of this tool.\n6. In the lupus IFN-β and COVID-19 cohorts, analyses appear to be conducted primarily at the cell level with minimal explicit modelling of donor variation, clinical covariates, or batch structure, raising concerns about confounding and dominance by high-cell-count samples; donor- or sample-level summaries (pseudo-bulk/mixed-effects, leave-one-donor-out checks) are needed to demonstrate cohort-level robustness of key signals.\nMinor Comments\n1. The GAT model requires 10 iterations of training to ensure robustness. Please provide time and memory complexity benchmarks for datasets exceeding 10,000 cells to ensure the tool is scalable for current large cohort.\n2. Indicate in the Methods how hyperparameters were selected (grid search, prior experience, or ad hoc choices).\n3. The reported github link and code is not available.", "reviewer_2": "This study introduces scDNS, a novel computational framework for identifying dysregulated gene regulatory interactions and prioritizing effector genes and cells in single-cell data. The authors demonstrate that scDNS can identify perturbed genes and cell populations across diverse biological contexts. However, several aspects of the methodology and interpretation require further clarification.\n1. scDNS incorporates multiple interaction databases to construct the prior GIN, but the rationale for selecting these specific databases is unclear. Are these standard resources in the field, or were they chosen based on empirical testing? More importantly, the authors should demonstrate that incorporating the prior GIN actually improves performance through ablation analysis. Specifically, compare scDNS performance with and without the prior GIN (i.e., using only single-cell expression data) to justify the added complexity of integrating these databases.\n2. When defining edges, scDNS uses the top 100 marker genes and housekeeping genes with high expression frequency. It is unclear whether 100 marker genes are sufficient or not. Sensitivity analysis is required here. Also, the negative edges are constructed using zero expression or random sampling, but I wonder if they could include some undetected interactions, and whether that would affect the performance of scDNS.\n3. The authors train the GAT model 10 times with different random seeds and aggregate attention scores. However, it is unclear whether 10 runs are sufficient for stable estimation. The authors should report the variability of attention scores across runs (e.g., coefficient of variation or standard deviation) for both simulated and real data. This analysis would help users understand the robustness of the method and guide appropriate parameter settings.\n4. scDNS computes both joint probability densities p(Gᵢ, Gⱼ) and conditional densities p(Gᵢ | Gⱼ), claiming each has distinct advantages. However, since conditional densities are derived from joint densities, they contain overlapping information. The authors should perform ablation studies comparing performance using (a) joint density only, (b) conditional density only, and (c) both combined. They should provide examples where each density type captures different biological signals, and clarify the claim that each performs differently in \"high-density vs. low-density cell populations\".\n5. scDNS constructs cell-type-specific networks, making its results inherently dependent on cell type definitions. However, the robustness to cell type annotation choices is not evaluated. The paper uses Seurat standard clustering method, but the clustering could be very different based on different parameters or methods being used. The authors should demonstrate the findings from scDNS are reproducible across reasonable variations.\n6. The authors should report the computational runtime and memory requirements of scDNS on datasets of varying sizes. As single-cell datasets increasingly exceed hundreds of thousands cells, users need to understand how runtime scales with the number of cells and genes, and whether the method is feasible for large-scale datasets.\n7. The interpretation of gene rankings is inconsistent throughout the manuscript. For example,\n• Line 206: IRF4 is called the \"top perturbed regulator\" but ranks 273\n• Lines 266-268: JAK1 and WTAP are highlighted but have moderate ranks\nThe authors should clearly define what rank/p-value cutoff constitutes a \"top\" or \"significant\" signal. How should users interpret rankings (is rank 273 considered high)?\n8. For the COVID-19 analysis, scDNS identifies disease-driven perturbation signals, but the consistency of these signals across patients is not clear. If scDNS is applied to individual patients separately, are the identified signals replicated across patients? What is the degree of inter-patient variability in perturbation signals? Is such variability technical, or driven by real biological variability? It would be nice if the authors can show that scDNS can identify patient-specific dysregulation patterns that relate to disease severity or outcome.\n9. Both the COVID-19 and IFN-β PBMC analyses utilized datasets with batch effects, which may introduce technical variation. While the authors apply Harmony correction for UMAP visualization, scDNS operates on uncorrected count data. Is scDNA robust under such circumstances?"}, "author_response": "We sincerely thank all Reviewers for the careful evaluation of our manuscript and for the\nconstructive suggestions. We have substantially revised the manuscript to address the concerns\nand to improve clarity, robustness, and reproducibility. For ease of navigation, responses to\nReviewer 2 start on page 15. Reviewer’s comments are marked in black while our replies are\nindicatedinblue.\nOverviewofthemajorrevisions:\nIn this revision, we (i) expanded the simulation engine with biologically structured perturbation\nmodes and repositioned simulations as stress tests (rather than head-to-head rankings), so that\nthey serve to validate scDNS under diverse rewiring scenarios; accordingly, we removed\ncross-method comparisons in these simulations, since the compared tools were developed for\ndifferentobjectives.(ii) strengthened cohort-level robustness analyses in both the COVID-19 and\nIFN-β PBMCs cohorts; (iii) added additional validation for representative DGR candidates\nincluding TIMM44; (iv) performed sensitivity/ablation analyses for key design choices (prior GIN\nintegration, marker-gene selection, joint vs. conditional density components, and\nclustering/annotation robustness); and (v) reported computational benchmarks and\nimplementationdetailstoimproveusabilityandscalability.\nImportantly,theseadditionsalsoallowedustoarticulateacoreadvantageofscDNSthatwasnot\nsufficiently highlighted in the original submission: unlike conventional methods that primarily\nprioritizeeither dysregulated regulatorsor differential/abundantcell populations, scDNS bridges\nthis gap by enabling an integrated interrogation of regulatory rewiring and its cellular\nconsequences. By linking gene-level interaction rewiring to cell-state–resolved perturbation\ndynamics,scDNShelpsrevealhowspecificregulatorsandprogramsmaydrivecomplexbiological\nprocessessuchasimmuneactivationanddiseaseprogression.OurCOVID-19analysisexemplifies\nthis capability: scDNS not only identifies high-perturbation gene programs, but also localizes\nthese signals to the specific cell populations and donors in which they are enriched, thereby\nconnectingmoleculardysregulationtointer-patientheterogeneityandclinicalseverity.\nWe believe these revisions substantially strengthen the manuscript and clarify where scDNS\nprovides value beyond expression- or abundance-driven analyses. Below, we provide a detailed\npoint-by-pointresponsetoeachcommentandindicatethecorrespondingchangesintherevised\nmanuscript.\nREVIEWERCOMMENTS\nReviewer#1(RemarkstotheAuthor):\nCurrent strategies for scRNA-seq analysis focus on differential expression analysis methods\ncoupled with pathway enrichment analysis. However, genes with unchanged expression yet\nprofoundly modified regulatory associations remain hard to detect. To address this knowledge\ngap, in this manuscript, Huang et al. introduces scDNS, a computational framework designed to\nquantify \"network rewiring\" in single-cell data. The authors combine Graph Attention Networks\n(GAT) with Jensen-Shannon Divergence (JSD) to identify regulators whose functional impact\nstems from altered interaction patterns rather than changes in transcript abundance. The\nconceptisquitenovel,which aimstoidentify keygenes playingimportant rolesduring acertain\n1\naltered condition without showing the commonly assumed expression changes. The study is\nwell-structured, and provided benchmark against existing tools (VIPER, Milo, MELD, etc.). The\ntool also demonstrate broad utility across diverse research settings, including primary\nimmunodeficiency (IRF4 mutation), viral infection (IFN and COVID-19), and oncology\n(GemcitabineresistanceinPDAC).TheexperimentalvalidationofTIMM44asatherapeutictarget\nisalsoahighlightdemonstratingthepowerandrobustnessofthiscomputational tool.However,\nseveralmajortechnicalandcomputationalpointsrequireclarificationbeforepublication.\nMajorComments\n1.Simulatedperturbations\nIn methods “Simulation setup for benchmarking datasets”, the scDNS tool introduces\nperturbations in three ways (shuffle, multiply, mix). While this is ok, this form of perturbation\nresults in dramatic statistical divergence that may not accurately reflect the structured,\nnon-randomrewiring ofregulatorynetworksobserved inactual pathologicalconditions.Can the\nauthors add more ways of perturbation to mimic the diverse scenarios of gene changes in\nbiological settings? In addition, only one single cell dataset is used, it would beneficial if the\nauthorscouldincludemorescRNA-seqdatasets.\nResponse: We thank the reviewer for this insightful critique. We agree that, under real\npathological or stimulation conditions, gene-level changes are rarely “random perturbations.”\nInstead, they are typically driven by alterations in gene/protein function (e.g., transcription\nfactors, signaling mediators, and post-translational regulators) and manifest as structured,\nnon-random rewiring of regulatory programs with clear modularity and context dependence.\nAccordingly, without changing the core goal of scDNS, we now systematically expand our\nsimulation engine and redesigned the perturbations into a set of biologically motivated\nperturbationmodes.\nAll perturbation modes follow the same design principles: (i) perturbations are confined to\nselected regulators and edges within their local neighborhoods and/or between functional\nmodules; and (ii) when adding new targets/edges, candidates are not sampled uniformly from\nthegenome,butareinsteaddrawnfromthesamefunctionalmodule,pathway,orregulon(i.e.,a\nrelated regulatory program), thereby better reflecting structured reprogramming observed in\nrealisticbiologicalsettings.\nBased on common biological scenarios, we integrated the perturbations into the following\ncategories:\n1) Expression-levelperturbationwithouttopologicalrewiring\nWe simulate up- or down-regulation of key factors (e.g., copy-number–driven dosage\neffectsortranscriptionalchanges).Withoutexplicitlyalteringnetworktopology,weonly\nmodify the expression distribution of the perturbed nodes, which enables us to assess\nscDNSsensitivityandrobustnesstoexpression-drivendifferences.\n2) Condition-dependentswitchingofregulatoractivityandmode\nTo capture condition-dependent regulatory function, we simulate regulators whose\n2\nactivity is weak or absent at baseline but becomes active under stimulation/pathology\n(inactive→active; gain-of-activity), or conversely switches from active to inactive\n(active→inactive; loss-of-activity), consistent with post-translational or\nsignaling-mediated activation (e.g., phosphorylation) that does not necessarily require\nchanges in mRNA abundance. In addition, we model logic rewiring driven by altered\ncofactorsorcomplexcomposition,suchthatselectedregulatoryrelationshipswitchfrom\nactivation to repression (or vice versa), reflecting changes in regulatory sign rather than\nmerelymagnitude.\n3) Regulatorytargetrewiringandinteraction-strengthrescaling\nWe model scenarios in which a regulator’s expression remains largely unchanged, but\nmutationsoraltered binding specificity/activity leads to partial lossof regulation over a\nsubset of targets and neomorphic gain over new targets (loss / neomorphic gain). For\nexample, mutations (e.g., IRF4) may shift DNA-binding preferences and reprogram the\ntargetrepertoire, producing simultaneouslossof original targetsand acquisition ofnew\ntargets.\nSeparately, while preserving the primary connectivity, we continuously modulate the\nstrengths of selected edges or local subnetworks (increase/decrease), mimicking\n“strength-type” perturbations arising from changes in binding affinity, cooperative\ninteractions,orchromatinaccessibility.\n4) Modularreconfigurationwithcell-state-specificpartialpenetrance\nWe simulate pathway-level reconfiguration in disease, in which previously coordinated\nfunctional modules become decoupled and/or develop aberrant couplings to other\nmodules,yieldingamorerealisticmodularrewiringpattern.Inimplementation,wealter\nonly a limited number of inter-module edges while keeping intra-module structure as\nstableaspossible,therebyavoidingunrealisticglobalrandomization.\nTo reflect clonal heterogeneity and varying stimulus responsiveness, perturbations can\nfurther be restricted to specific cell states/subpopulations and/or applied to only a\nfraction of cells (partial penetrance). We systematically scan penetrance together with\nperturbation proportion and strength to evaluate scDNS performance under\nheterogeneousconditions.\nFinally, we explicitly construct complex settings by combining expression-level perturbations\n(category1)withtherewiringmodesincategories2–4,therebycoveringthreekeycategories:(i)\nexpression-onlychanges,(ii)rewiring-onlychanges(regulatoryrelationships/expressionpatterns),\nand(iii) coupled expression +rewiring changes. The local-neighborhood/module constraintsand\nbounded perturbation design ensure that the resulting changes remain structured and\nbiologicallyplausible,enablingbenchmarkingofscDNSunderrealisticcomplexity.\nIn revision, we use a unified parameterization to demonstrate robustness across diverse\nperturbation regimes. Concretely, for each mode we define a perturbed regulator set R,\nperturbation proportion r (fraction of affected edges/targets within a local neighborhood or\nfunctional modules), perturbation strength α (e.g., expression fold-change, edge-weight scaling,\nor correlation shift), and cell-level penetrance p (fraction of cells expressing the perturbed\nphenotype). Module constraints restrict candidate sets for new targets or new edges, ensuring\nthatperturbationsstructurallyconstrainedratherthangloballyrandomized.\nassessthegeneralizabilityofscDNS.AsshowninResponseFig.1,scDNSconsistentlydetectsthe\nimposed perturbations with high sensitivity and strong robustness across different perturbation\nmodesanddatasets.Wehaverevisedthemanuscriptaccordingly,providingdetaileddescriptions\nof the newly implemented perturbation schemes and reporting the corresponding evaluation\nresultsonallthreedatasets.\nResponse Fig. 1 | Benchmarking scDNS performance across diverse perturbation regimes and\ndatasets. a-c, Accuracy ofscDNS indetecting simulated perturbations acrossthree independent\nscRNA-seqdatasets:iPSC-derived alveolarepithelial typeII-likecells(a), humanembryonic stem\ncells (b), and PANC-1 cells (c). Performance is evaluated across three perturbation categories:\nexpression-level changes only (Expression level), regulatory rewiring/pattern changes only\n(Expression pattern), and combined effects (Mixed). Top panels: Accuracy as a function of\nperturbation strength (α), representing the magnitude of fold-change or correlation shifts.\nBottompanels:Accuracyrelativetothefractionofperturbation(%),representingtheproportion\nofaffectedcells.\n2.Algorithmicoverfit\nIn figure2e, DNS achieve an AUC value of near 1.0 in simulated data while other softwares are\naround 0.5, which is quite poor. This result alerts that there is an overfitting problem here. The\ncorepointofscDNSistoquantifythedifferencesinprobabilitydistributionsusingJSDdivergence.\nThe perturbation methods used, as also discussed above, directly maximizes the JSD value,\n4\ncreating an inherent \"home-court advantage\" for scDNS compared with other published\nmethods.\nResponse:Weappreciatethereviewer’sconstructivecritique.InthesimulationshowninFig.2e,\nthe perturbation modes were designed to generate three regimes—expression-only changes,\nrewiring-onlychanges,andcoupledexpression+rewiringchanges—whicharepreciselythetypes\nof distributional/network shifts that scDNS is designed to quantify via network divergence\nmeasures. As the reviewer correctly points out, this design can inadvertently create a\n“home-courtadvantage”forscDNS,andthereforetheresultingAUCinthissimulationshouldnot\nbeusedasahead-to-headrankingagainstmethodsthatweredevelopedfordifferentobjectives.\nImportantly, our goal for the simulation experiments was not to claim superiority over existing\ntools, but to demonstrate that scDNS can sensitively recover perturbed cells across a range of\nbiologically motivated perturbation scenarios. In contrast, the compared methods address\nfundamentallydifferenttasksand rely on different assumptions—forexample, VIPER focuses on\nregulon-based protein activity inference; MELD estimates cell-wise condition-associated\nlikelihood/enrichment on a manifold; and Augur prioritizes cell types by how well they\ndiscriminateconditions—soadirectAUCcomparisonunderthisscDNS-orientedsimulationisnot\nmethodologicallyfair.\nTo address this concern, in the revised manuscript we have removed the cross-method AUC\ncomparison on the simulated benchmark and instead use the simulation only as a stress test to\nillustrate how scDNS responds to diverse, structured perturbations. We now report scDNS\nperformance trends across perturbation strength and penetrance, and we clarify that these\nsimulations are intended to evaluate scDNS’s ability to detect divergence-driven regulatory\nabnormalitiesratherthantobenchmarktoolsoptimizedforotherobjectives.\nCrucially,ourmaincomparativeevidenceisbasedonrealdatasetswithexternalgroundtruth(Fig.\n2f–j), where scDNS is not advantaged by simulation design. In the revision, we emphasize that\ncross-methodcomparisonsandperformanceclaimsaredrawnfromthesereal-databenchmarks,\nandweexplicitlydirectreaderstothoseevaluationsforfairmethodcomparison.\n3.Benchmark\nFollowingthepreviousquestion,infigure2e-gbenchmarking,Milo andscDistshowednear-zero\nsensitivityinidentifyingperturbedsubpopulationsinspecificDGR(DifferentialGeneRelationship)\nscenarios, while scDNS perform very good. The baseline methods chosen for comparison (e.g.,\nMilo,Augur,scDist,VIPER)werenotoriginallyengineered todetectstatisticaldivergencecaused\nby random shuffling. Milo and scDist are optimized for detecting shifts in neighborhood\nabundance or inter-condition distances. VIPER relies on established regulon structures to infer\nprotein activity. Evaluated against a perturbation scheme thatmirrors the specific mathematical\nlogicofscDNS,theirperformancedeclineisexpected.\nResponse: We thank the reviewer for this comment and agree with the interpretation.\nBenchmarking these DGR-oriented simulations against methods with fundamentally different\n5\nobjectivescaninadvertentlyfavorscDNSandyieldapparentlylowsensitivityforothertools—not\nbecause they “fail,” but because they are not designed to detect divergence in gene–gene\nrelationships.Accordingly,theirlowsensitivityincertainDGRsettingsisexpectedandshouldnot\nbeinterpretedasafaircross-methodranking.\nConsistentwith ourresponseabove,wehavethereforeremovedthecross-methodcomparisons\nin Fig. 2e–g in the revised manuscript and now use these simulations only as “stress tests” to\ndemonstrate scDNS sensitivity and robustness across perturbation modes, strengths, and\npenetrance. To assess generalizability, we further extended the simulation benchmarks to\nmultiple scRNA-seq datasets rather than relying on a single dataset. Cross-method comparisons\nand performance claims are instead anchored in real-data benchmarks with external ground\ntruth (Fig. 2f–j), where the simulation design does not advantage scDNS. We have revised the\ntext accordingly to guide readers toward these ground-truth-based evaluations for objective\ncomparisons.\n4.Potentialbiasofpriorknowledgenetworks\nThescDNS frameworkreliesonintegratingestablished databases(TRUST,STRING,KEGG,etc.)to\nconstruct the initial reference network. There is a concern regarding \"research bias\": the model\nmight be predisposed to find perturbations in well-studied genes, while \"orphan genes\" with\nfewer known interactions remain undetectable. How about the authors perform an experiment\nby randomly removing varying percentages (e.g., 20%-50%) of prior edges to assess the\nrobustnessofidentifyingkeyregulators(e.g.,TP53)?\nResponse: We thank the reviewer for this insightful comment. We agree that the density of\nknowninteractionsvariesacrossdatabases.Wehaveaddressedthepotentialfor\"researchbias\"\nthrough algorithmic design, existing benchmarks, a proactive discussion of limitations, and the\nspecificrobustnessexperimentsuggested.\n1) Algorithmicmitigationofdegreebias\nFromtheoutset,scDNSwasdesignedtoaccountforheterogeneityinregulatorymodulesize(i.e.,\nnodedegreeinthepriornetwork).Topreventhigh-degreegenesfromdominatingtheresults,we\napply degree normalization: rawdivergence scores are converted to degree-normalized Z-scores\nusinganulldistributionderivedfromrandomizednetworks.Thisdesignensuresthatsignificance\nis assessed relative to the expected behavior of nodes with the same degree, rather than being\ndrivenbytheabsolutenumberofedges.\n2) Sensitivityanalysisacrossgenedegrees\nIn our original submission, we directly evaluated the impact of network connectivity on\nperformance. Asshown in Fig. S2i, scDNS maintainsstableperformanceacrossabroadrangeof\ngene degrees. A noticeable drop occurs only when a gene has very few interaction partners\n(degree < 5), where the lack of local topological information makes reliable rewiring inference\nstatistically challenging. This indicates that scDNS is robust for the vast majority of genes and is\nnotdrivensolelyby“high-degree”nodes.\n3) Mitigatingincompletenessviamulti-omicsanddata-drivenpriors\nWhile human gene interaction networks are extensively studied, we acknowledge that some\ngenes and pathways remain under-characterized. We already noted this limitation in the\ndata-driven, context-specific evidence, such as ATAC-seq co-accessibility, ChIP-seq binding, or\nothermulti-omicspriors,toimprovenetworkreconstruction. Inaddition,usersmayintegratede\nnovo regulatory relationships inferred by single-cell network inference methods as\ncomplementary priors. These integration strategies could extend scDNS to genes currently lack\ndocumentedinteractionsinpublicdatabases.\n4) Robustnesstoincompletepriorknowledge(newexperiment)\nAs suggested by the reviewer, we performed an additional robustness test to assess whether\nscDNS can still recover key regulators when prior knowledge is substantially missing. We\nrandomly removed 25% - 6%0 of edges incident to a representative regulator (e.g., TP53) from\nthe prior network and re-ran scDNS. As shown in Response Fig. 2, scDNS consistently identified\nthe regulator with a significantZ-score, remaining robusteven under severe edge removal (e.g.,\nZ-score > 3.5 at 60% removal). These results indicate that scDNS tolerates substantial\nincompletenessinpriornetworksandcanstilldetectunder-annotated(“orphan”)genes.\nResponse Fig. 2 | Impact of prior network incompleteness on the identification of TP53\ndysregulation.Thex-axisrepresentsthepercentageofremovededges,rangingfrom25%to60%\nwitha5%increment.BoxplotsillustratethedistributionofZ-scoresacross50runs.\n5.DGRsvalidation\nThe scDNS can detect both genes with altered expression levels (DGE) and, interestingly, genes\nwith unchanged expression levels but altered regulatory networks (DGR). In figure4, JAK1 and\nWTAP is identified as DGRs, while in figure6, TIMM44 is also considered as a DGR. These gene\nshow a little mRNA change but high network perturbation. Some experimental validation in a\ncertainsettingwillsignificantlystrengthentheanalyticalpowerofthistool.\nResponse:Wethankthereviewerforthisconstructivesuggestion.AcentralcapabilityofscDNSis\nto detect not only differential gene expression (DGE), but also differential gene relationships\n(DGRs)—cases where mRNA abundance remains largely unchanged while regulatory programs\narerewired.\nFirst,thebiologicalplausibilityofDGRsissupportedbymultipleground-truthexamplesanalyzed\nin our manuscript, including regulators whose functional activity can change without marked\n7\ntranscriptional shifts (e.g., TP53, IRF4, IFNAR). In addition, the specific DGR candidates\nhighlightedbythe reviewer—JAK1 and WTAP—further exemplifythis. Forinstance, JAK1 activity\ncan be modulated by phosphorylation/autophosphorylation, enabling substantial signaling\nchanges without requiring large changes in transcript level. Likewise, WTAP can undergo\nstimulus-dependentcondensateremodeling/phase behavior changes thatimpact itsm6A writer\ncomplex function, providing a mechanistic route to strong regulatory effects with limited mRNA\nchange.\nSecond,forournewlyprioritizedDGRcandidateTIMM44,ouroriginalmanuscriptconfirmedthat\nitsmRNAandproteinabundanceremainstableundergemcitabinetreatment.Toprobewhether\nTIMM44 undergoes an activity/interaction-state change, we performed complementary\nproteomics and phosphoproteomics profiling, as phosphorylation is one of the most prevalent\nand rapidly inducible mechanisms for modulating protein function and protein–protein\ninteractions. The new data (Response Fig. 3a) confirm that total TIMM44 protein levels remain\nstable,whereasphosphorylationatS180issignificantlyincreased aftertreatment.Thisfindingis\nparticularly striking given that TIMM44 serves as a critical tether for the mitochondrial Hsp70\n(mtHsp70) motorcomponents. Existing biochemicalevidence (e.g., Schiller etal., 2008. Mol Cell\nBiol.PMID:18426906;Tingetal.,2017.Elife.PMID:28440746)confirmsthatS180resideswithin\na functional domain essential for associating TIMM44 with the translocon and recruiting\nmtHsp70.\nMoreover, data from PhosphoSitePlus identify S180 as a recurrently detected phosphorylation\nsiteacross27independentstudies(ResponseFig.3b),underscoringitsstatusasabonafideand\nphysiological regulatory site. To evaluate its functional significance, we performed multi-species\nsequencealignmentandevolutionaryconservationanalysis.TheresultsrevealthatS180ishighly\nconserved across diverse taxa, suggesting it is subject to strong evolutionary constraints\n(Response Fig. 3c). Consistent with this conservation, AlphaMissense analysis predicts that\nsubstitutions at S180 are \"likely pathogenic\" (Response Fig. 3d). This high score reinforces the\nview that S180 is a critical functional hotspot, where biochemical alterations are likely to have\nsignificantbiologicalconsequences.\nWhiletheseobservationsdonotdirectlyestablishcausality,theyprovidecompelling,convergent\nevidence that S180 is a regulated and functionally indispensable residue. Such findings strongly\nsupport the hypothesis that S180 phosphorylation modulates the functional or interaction state\nof TIMM44, thereby providing a mechanistic basis for the DGR signal captured by scDNS.\nAlthough other regulatory mechanisms—such as alternative post-translational modifications or\nprotein complex remodeling—may also be at play, our phosphoproteomic profiling offers\northogonal validation that scDNS can effectively prioritize functionally relevant regulators that\nremain\"hidden\"instandardexpression-onlyanalyses.\nresponsive to gemcitabine. a, Protein-normalized TIMM44 pS180 abundance and total protein\nlevels. Left: Relative abundance of TIMM44 S180 phosphorylation after gemcitabine (GEM)\ntreatment, normalized to total TIMM44 protein abundance. Right: Relative total protein\nabundance of TIMM44 under control (DMSO) and gemcitabine (GEM.) (1 μM) conditions. Data\nare presented as mean ± SEM (n = 3); **p < 0.01; ns, not significant (Student’s t-test). b,\nPhysiological evidence of S180 phosphorylation from PhosphoSitePlus. c, Evolutionary\nconservationandsecondarystructureoftheTIMM44S180region. Multiplesequencealignment\nof TIMM44 orthologs from diverse species (Human to Yeast) was performed using MUSCLE.\nSecondary structure elements (top) are mapped from the human TIMM44 AlphaFold structure\n(O43615). d, In silico pathogenicity prediction of TIMM44 S180 mutations. Table summarized\nAlphaMissenseandREVELscoresforallpossibleaminoacidsubstitutionsatSer180.\n9\n6.InthelupusIFN-βandCOVID-19cohorts,analysesappeartobeconductedprimarilyatthecell\nlevel with minimal explicit modelling of donor variation, clinical covariates, or batch structure,\nraising concerns about confounding and dominance by high-cell-count samples; donor- or\nsample-levelsummaries(pseudo-bulk/mixed-effects,leave-one-donor-outchecks)areneededto\ndemonstratecohort-levelrobustnessofkeysignals.\nResponse: We appreciate the reviewers’ comments and would like to clarify the theoretical\nrationale behind scDNS’s relative robustness to batch effects and donor heterogeneity. Unlike\nabundance-based methods thatprimarilyrelyon meanexpression orneighborhood abundance,\nscDNS focuses on changes in gene–gene dependency structure and summarizes divergence\npatternsoverregulatorymodules,whichprovidesseveralinherentadvantagesforheterogeneous\ncohorts:\n1) Relativerobustnesstotechnicalartifacts\nscDNS quantifies differential gene relationships (DGRs) by capturing mutual information–based\ndependency shifts and network-topological divergence. From a theoretical standpoint, common\ntechnicalconfounders—suchasvariationinlibrarydepth,globalexpressionshifts,andstochastic\nmeasurement noise—primarily affect absolute expression levels (first-order statistics). While\nsuchfactorscanaddnoisetotheestimationofgene–genedependencies,theyaregenerallyless\nlikely to induce coherent, directionally consistent rewiring patterns across a regulatory network.\nTherefore,scDNSiscomparativelylesspronetofalsepositivesdrivenpurelybybatchstructureor\nsequencingdepth.\n2) Signalreinforcementthroughcross-donorconcordance\nThescDNSframeworkprioritizessignalsthatareconsistentlysupportedacrossthecohort.When\naregulatorychangeispresentinmanydonors,theresultingstatistics(e.g.,Z-scoresandp-values)\nbecome more significant because concordant divergence patterns reinforce one another. In\ncontrast,deviationsconfinedtoonlyasmallsubsetofdonorswillgenerallycontributelesstothe\ncohort-wide ranking. Importantly,this does not implythatscDNS ignoresinter-patientbiological\nheterogeneity:scDNSproducescell-state–resolvedperturbationreadoutsthatmakeinter-patient\nvariability explicit and enable downstream stratified analyses linking patient-specific patterns to\nclinical covariates (e.g., disease severity), as we demonstrate in our response to Reviewer #2\n(Comment#8).\n3) Interpretabilityandvisualizationofpatientheterogeneity.\nImportantly, scDNS produces cell-resolved perturbation scores that can be aggregated at the\nsample level and projected onto UMAP embeddings. This provides an interpretable way to\nexaminehowperturbationsignalsdistributeacrossdonorsandclinicalstrata.Byinspectingthese\nscore maps, users can readily determine whether a perturbation is broadly shared across the\ndisease state or concentrated in specific donors or covariate-defined subgroups, enabling\ntransparentposthocassessmentofpotentialconfoundingorclinicallymeaningfulheterogeneity\n(e.g.,patternsrelatedtodiseaseseverityoroutcome).\nconcerns, we performed additional donor-aware analysesin both the lupusIFN-βand COVID-19\ncohorts. These datasets provide complementary ground-truth anchors: in the lupus cohort, key\nIFN-responsive genes and pathways are well established, allowing evaluation of whether donor\nheterogeneity or covariates distort expected signals; in the COVID-19 cohort, disease-relevant\nepithelial populations (AT1/AT2) serve as a known biological reference, enabling assessment of\nwhetherinferredperturbed cellpopulationsareconsistentunderdonorvariabilityandcovariate\nstructure.Weassessedrobustnessfromtwoperspectives:\n(1)Cross-donorconsistency.\nIn the lupus dataset, we aggregated scDNS scores at the donor level by averaging scZ-scores\nwithineachdonorandcelltype,focusingoncanonicalinterferon-responsegenes.Weobserveda\nconsistent elevation of the IFN-response signal (e.g., IFNAR1, a direct target of IFN-β) across\ndonors, aligned with the expected direction of change between control and treatment groups\n(ResponseFig.4a),indicatingthatthekeyperturbationsignalsarenotdrivenbyasmallsubsetof\nhigh-cell-countsamples.\nInthe COVID-19dataset, wesimilarly computed donor-level aggregated scores(mean scZ-score)\nacross cell types. Although one control sample (C53ctr) showed an outlying score relative to\nothercontrols,theoverallcohort-leveltrendwaspreserved:AT1andAT2populationsexhibiteda\nconsistentreductionin aggregatedscZ-scoresacrossthemajority ofpatients(ResponseFig.4b).\nThese donor-level summaries provide a direct view of inter-patient variability and support that\nthemainscDNSsignalsreflectbroadlyconsistentcohortpatternsratherthansingle-donoreffects.\nMoreover, the inter-patient variability captured by scDNS can also prompt further clinical or\ntechnicalfollow-up. Forexample,asdescribed inourresponseto Reviewer 2(comment#8), the\npatient-to-patient variation in the fraction of high-scDNS-score cells is associated with disease\nseverity, supporting the biological and clinical relevance of these heterogeneous perturbation\npatterns.\n(2)Covariaterobustness.\nIn both cohorts, we examined whether scDNS scores were associated with common technical\ncovariates, including total UMI counts (nCounts) and number of detected genes (nFeature).\nSimple linear fits showed no significant correlation between these technical covariates and\nscZ-scores (Response Fig. 4c&d). We further employed a multivariable regression framework to\nassess potential confounding by technical and sample-level metadata. Specifically, we modeled\ndonor-aggregated scDNS scores as a function of technical covariates (nCounts, nFeature) and\ndemographic covariates (age, sex). After adjustment, none of these covariates was a significant\npredictorofscDNSscoresineithertheCOVID-19orlupuscohorts(small∣t∣ statistics; >0.1),\nsuggestingthatthescDNSsignalisnotmateriallydrivenbysequencingdepth,featuredetection,\n\nor demographic differences (Response Fig. 4e&f; t-statistics computed from the fitted model\nafteradjustment).\nInsummary,theseadditionalanalysessupportthatscDNScanrobustlyrecovercoreperturbation\nsignals in heterogeneous human cohorts, while providing an interpretable, cell-resolved\n11\nframeworktoinspectdonor-specificeffectsandpotentialclinicalassociations.\nResponse Fig. 4 | Robustness and donor-level consistency of scDNS scores. a, Donor-level\naggregationofscDNSscoresinthelupusdataset.EachdatapointrepresentsthemeanscZ-score\nof canonical interferon-response gene, IFNAR1, for an individual donor (n = 8donors). Different\nsymbolsdenoteuniquedonorIDs.b,Comparisonofdonor-levelmeanscZ-scoresacrosslungcell\ntypes in the COVID-19 dataset. Red and blue dots represent control and COVID-19 donors,\nrespectively. Dark red and blue linesindicate the cohort average. Cell typeswere ranked by the\nmaximum absolute value of their average scZ-score across Control and COVID-19 groups. c,d,\nCorrelationanalysisbetweentheabsolutescZ-scoreandtechnicalcovariates:totalUMIcounts(c)\nandthenumberofdetectedgenes(d).Bluelinesrepresentsimplelinearfits;Pearsoncorrelation\ncoefficients (R) and P values are shown (two-sided test). e,f, Multivariable linear regression\nanalysisofdonor-aggregatedscZ-scoresintheCOVID-19(e)andlupus(f)cohorts.Barplotsshow\nt-statistics for the indicated clinical and technical variables. The color gradient indicates\nlog -transformedPvalues.\n1. The GAT model requires 10 iterations of training to ensure robustness. Please provide time\nandmemorycomplexitybenchmarksfordatasetsexceeding10,000cellstoensurethetoolis\nscalableforcurrentlargecohort.\nResponse: We have included comprehensive benchmarking data in the revised manuscript to\n12\naddressconcernsregardingthemodel'sstabilityandscalability.\nTo ensure the robustness of the GAT component, we performed 10 independent training runs\nand quantified the reproducibility of the attention scores. As shown in Response Fig. 5, the\npairwise Spearman correlations of attention weights across independent runs are consistently >\n0.95. This high degree of concordance indicates that the attention mechanism reliably captures\nstableregulatorysignals,andtheinferredpatternsarenotsensitivetostochasticinitialization.\nWebenchmarkedscDNS using theIFN-β-treatedPBMCs dataset(10,194 cellsand 13,431 genes)\nonastandardhigh-performancecomputingnode(IntelXeonGold6226RCPU,768GBRAM,and\na single NVIDIA RTX A6000 GPU). The wall-clock time and peak memory consumption for the\nthreeprimarystagesaresummarizedbelow:\nTable1.RuntimeandmemoryusageofscDNSacrosskeycomputationalstages\nComputationalStage Runtime PeakMemoryUsage Platform\nI.Preprocessing&NetworkConstruction ~8.5min 21.4GB(RAM) CPU(R)\nII.NodeEmbeddingPreparation ~6.5min 1.17GB(VRAM) GPU(Python)\nIII.GAT+DGITraining(10iterations) ~87min 22.5GB(VRAM) GPU(Python)\nIt shows that the preprocessing and embedding steps(I & II) are comparatively lightweight. The\nmostcomputationally intensive stage of scDNS is GAT+ DGI training (III), which operates on the\ngene-level graph and embeddings. Consequently, once the initial embeddings are computed,\nruntimeandmemoryaredrivenmainlybythenumberofgenesandgraphdensity(connectivity)\nrather than the number of cells, enabling scDNS to scale efficiently to large single-cell datasets\nexceeding10,000cells.\nResponseFig.5|Reproducibilityofattentionweightsacrossindependentruns.Heatmaps\nshowingSpearmancorrelationcoefficientsofattentionweightsbetween10independentmodel\nruns(labeledR1–R10)forfourdistinctcelltypes:IFN-β-treatedPBMCs,PANC-1cells,Bcells,and\niAT2cells.Thecolorscaleindicatesthestrengthofthecorrelation,rangingfrom0.9(white)to\n1.0(darkred).\n2. Indicate in the Methods how hyperparameterswere selected (grid search, prior experience,\noradhocchoices).\nResponse: Hyperparameters were selected using domain-standard defaults rather than\nexhaustive grid search. Specifically, for data preprocessing and prior-network construction, we\nfollowestablishedbestpracticesinsingle-cellanalysistoensurereproducibility.FortheGAT+DGI\n13\ncomponent, the architecture is intentionally kept lightweight and aligned with widely used\nconfigurations in graph representation learning, rather than being extensively tuned on a\nper-datasetbasis.\nWe adopted this strategy for three reasons. First, our reconstruction uses weak supervision via\nprior-derived proxy labels (positive/negative edges derived from prior knowledge and\nmarker-based anchors), butthere isno definitivegold-standard regulatory network availablefor\nmost single-cell datasets to serve as an objective target for large-scale hyperparameter\noptimization. In this setting, extensive grid search can encourage dataset-specific tuning that\nharms reproducibility and generalization. Second, heavy hyperparameter optimization can\ninadvertently increase the risk of overfitting to cohort- or batch-specific idiosyncrasies, whereas\nour downstream goal is to quantify robust divergence patterns. Third, using a stable default\nconfigurationimprovesusabilityandmakesresultsmorecomparableacrossdatasetsandstudies.\nIn practice, we verified that the default settings yield stable training dynamics and complete\nconvergence. Specifically, as shown in Response Fig. 6, both the supervised edge-classification\nlossandtheDGIself-supervisedlossdecreaserapidlyandplateauafter~300epochs;importantly,\nthe combined (weighted) total loss shows the same convergence behavior. Additional epochs\nmainlyprovideasafetymarginanddonotmateriallyalterthelearned representation.InscDNS,\nthe default number of epochs was set to 1,000 per run to ensure reliable convergence across\ndiversedatasetswithoutrequiringdataset-specifictuning.\nFor transparency and reproducibility, we have added an explicit description in the Methods\ndetailing this selection strategy, and we have documented all hyperparameters—including\ndefinitions and default values—in the software documentation (GitHub: AttnEdgeMiner). All\nparametersremainfullyconfigurable,allowinguserstofine-tunesettingsforspecificdatasetsor\nbiologicalcontextsifdesired.\nResponse Fig. 6| Convergence of model training. Representative training loss curves for the\nGAT+DGI framework over1,000 epochs. The supervised taskloss (blue), self-supervised DGI loss\n(orange),andthecombinedtotalloss(green)alldemonstraterapidminimization.\n3.Thereportedgithublinkandcodeisnotavailable.\nResponse: We apologize for the inconvenience. The GitHub repository has now been made\npublicly accessible (https://github.com/xiaolab-xjtu/scDNS), and the manuscript has been\n14\nupdatedwiththecorrectedhyperlink.\nReviewer#1(Remarksoncodeavailability):\nThereportedgithublinkandcodeisnotavailable.\nResponse:thisissuehasbeenresolved.\nReviewer#2(RemarkstotheAuthor):\nThisstudyintroducesscDNS,anovelcomputationalframeworkforidentifyingdysregulatedgene\nregulatory interactions and prioritizing effector genes and cells in single-cell data. The authors\ndemonstrate that scDNS can identify perturbed genes and cell populations across diverse\nbiologicalcontexts.\nHowever,severalaspectsofthemethodologyandinterpretationrequirefurtherclarification.\n1.scDNSincorporatesmultipleinteractiondatabasestoconstructthepriorGIN,buttherationale\nforselectingthesespecificdatabasesisunclear.Arethesestandardresourcesinthefield,orwere\ntheychosenbasedonempiricaltesting?Moreimportantly,theauthorsshoulddemonstratethat\nincorporatingthepriorGINactuallyimprovesperformancethroughablationanalysis.Specifically,\ncompare scDNS performance with and without the prior GIN (i.e., using only single-cell\nexpressiondata)tojustifytheaddedcomplexityofintegratingthesedatabases.\nResponse: Wethank the reviewer for this important suggestion. We have clarified the rationale\nfor selecting the prior interaction resources and performed an ablation analysis to quantify the\ncontributionofthepriorGIN.\nRationalefordatabaseselection:TheuseofapriorGINisintendedtoprovideacomprehensive,\nwell-curated interaction scaffold that guides subsequent network refinement and construction.\nThisprior not onlyreduces the searchspace and computationalburden compared with denovo\nnetworkinference,butalsosuppliesahigh-qualityreferencetoimprovethebiologicalplausibility\nof the reconstructed network. Accordingly, we integrated (i) transcriptional regulatory\ninteractions (e.g., DoRothEA and TRUST), (ii) protein–protein interactions (e.g., IntAct and\nSTRING), and (iii) curated signaling/pathway knowledge (e.g., KEGG and SuperPathway). These\nare widely used, community-standard resources with established curation and quality-control\nprocedures.Similarregulatory/interactionpriorsarealsoroutinelyleveragedbyexistingmethods\n(e.g. VIPER). Among these, STRING and IntAct provide particularly broad and high-coverage\ninteraction evidence. The integrated prior covers >90% of total genes (Response Fig. 7a),\nensuringthatscDNSisnotdrivenbyasmallsubsetofwell-annotatedgenes.\nAblation analysis with vs. without prior GIN: As requested, we performed an ablation study on\nthreerealdatasetsbycomparing scDNS resultswith thepriorGINversusavariantthatrelieson\nsingle-cell expression information only (i.e., withoutintegrating the prior interaction graph). The\nresultsshowthatincorporating theprior GINsubstantiallystrengthenstheinferredperturbation\nsignals—most clearly reflected by higher and more decisive scDNS Z-scores for biologically\nexpected regulators/genes (Response Fig. 7b). This demonstrates that the database-derived\npriors provide meaningful, non-redundant information that improves sensitivity and confidence\nofscDNSoutputs,therebyjustifyingtheaddedintegrationstep.\ncategory,(ii) report coverage statistics, and (iii) include the ablation results to explicitly quantify\nthebenefitofthepriorGIN.\nResponseFig.7|GenecoverageanalysisandablationstudyofthepriorGINs.a,Genecoverage\nacrossdifferentdatabasesinthreecelltypes:Barplotsshowingthenumberofgenescoveredby\nindividualpriorknowledgeresources(Trrust,SuperPathway,KEGG,DoRothEA,IntAct,andSTRING)\nand the integrated prior for iAT2, H1-hESCs, and PANC-1 datasets. Dark blue bars represent\ncoveredgeneswithhighconnectivity(degree>10),lightbluebarsrepresentcoveredgeneswith\nlower connectivity (degree < 10), and grey bars represent uncovered genes. b, Ablation analysis\nofperturbationsignalstrength:BoxplotscomparingthescDNSZ-scoresobtainedwithoutaprior\nnetwork versus those obtained by integrating specific prior databases and the final integrated\nnetwork.\n2.Whendefiningedges,scDNSusesthetop100markergenesandhousekeepinggeneswithhigh\nexpression frequency. It is unclear whether 100 marker genes are sufficient or not. Sensitivity\nanalysis is required here. Also, the negative edges are constructed using zero expression or\nrandomsampling,butIwonderiftheycouldincludesomeundetectedinteractions,andwhether\nthatwouldaffecttheperformanceofscDNS.\nResponse: We thank the reviewer for raising these important points. The purpose of selecting\nmarker genes is to define a compact set of cell-type–informative anchor genes for constructing\npredefined,cell-type–specificedges.Inthisstep,thegoalisnottoexhaustivelyincludeallgenes,\nbut rather to include a sufficient number of highly characteristic genes that reliably capture the\nidentity and key programs of each cell type. To directly assess whether selecting 100 marker\ngenes is sufficient, we performed a sensitivity analysis on three datasets by varying the number\nof top marker genes used to define the predefined edges from 5 to 500 and quantifying the\nagreement of the learned GAT+DGI attention scores across settings. We found that once the\nmarker set reaches ~100 genes, the resulting attention scores become highly consistent with\nthose obtained using larger marker sets (Spearman’s ρ ≥ 0.8 for comparisons between 100\nand>100markers;ResponseFig.8a).Thisindicatesthatincreasing thenumberofmarkergenes\nbeyond 100 has little additional effect on the inferred attention patterns—especially under\nfull-scalenetworkreconstruction—andthushasminimalimpactondownstreamscDNSresults.\nRegarding negative-edge construction and the possibility of undetected interactions, scDNS\nconstructsnegativeedgesusingtwocomplementarystrategies.Zero-expression–basednegatives\narederivedfromgeneswith zeroexpression intherelevantcellcontext.Suchgenes areunlikely\nto participate in strong, cell-type–defining regulatory relationships in that context, and are\ntypically not key drivers of cell identity or function. Consistent with this, we analyzed the\nfunctionalcompositionoftheexcludedzero-expressiongenesacrossfourdatasets(IFN-βtreated\nPBMCs,pancreaticcancer(PANC-1cells),Bcells,andiAT2cells)andfoundthattheyareenriched\nforgenesthatarenotspecificallyrelevanttothecorrespondingcell-typeprograms(ResponseFig.\n8b). For random-sampling negatives, we draw gene pairs that are absent from the large\nbackground prior network. This is a standard negative-sampling strategy in graph machine\nlearningtoapproximatenonexistentinteractions.\nWhileitispossiblethatfewsuchpairsmightrepresentundiscoveredinteractions,thisisunlikely\nto materially affect performance in practice for two reasons: (i) the probability of repeatedly\nsampling a truly functional but entirely undocumented interaction is low; and (ii) both the\ntrainingobjectiveanddownstreamdivergenceestimationarerobusttoamodestamountoflabel\nnoiseinthenegativeset.\nMoreover, even if a small number of false negative edges were included (or if some predicted\nedgeswereimperfect), they areunlikelytoproducestrong, systematicdivergencesignals across\nconditions. scDNS also operates as an aggregation framework, further reducing sensitivity to a\nsmall fraction of mislabeled or spurious edges. Empirically, our robustness experiments support\nthis: in Supplementary Fig. 2g, we replaced a substantial fraction of true TP53 interactions with\nrandom incorrect edges and observed that scDNS maintained high accuracy even when 50% of\nedges were corrupted. In our response to Reviewer #1 (comment #4), we further show that\nscDNSremainsrobustwhenaportionoftrueedgesisremovedfromthepriornetwork.\nTakentogether,thesenewsensitivityandrobustnessanalysessupportthat(i)using~100marker\ngenes are sufficient and stable, and (ii) the negative-edge construction has minimal impact on\nscDNSperformanceanddoesnotunderminethereliabilityoftheinferredperturbationsignals.\nof excluded genes. a, Heatmap of pairwise Spearman correlations of attention scores obtained\nusing varying numbers of marker genes (top 5 to top 500) in IFN-β-treated PBMCs. Correlations\nstabilize(Spearman’s ρ > 0.80) once morethan 100 markersare included, supporting the use of\nthe top 100 marker genes for robust downstream analyses. b, Gene Ontology enrichment\nanalysis of zero-expression genes used to construct negative edges in four representative\ndatasets (IFN-β-treated PBMCs, PANC-1, B cells and iAT2 cells). The x-axis denotes enrichment\nsignificance(−log₁₀Pvalue),andthey-axislistsenrichedbiologicalprocesses.\n3. The authors train the GAT model 10 times with different random seeds and aggregate\nattentionscores. However,it isunclear whether 10 runsaresufficientforstableestimation. The\nauthorsshould reportthevariabilityof attention scoresacrossruns(e.g.,coefficientofvariation\nor standard deviation) for both simulated and real data. This analysis would help users\nunderstandtherobustnessofthemethodandguideappropriateparametersettings.\nResponse:Wethankthereviewerforthishelpfulsuggestionandagreethatreportingrun-to-run\nvariability is important forassessing robustness. As also raised by Reviewer #1 (Minor comment\n#1), we have now explicitly quantified the stability of the GAT+DGI attention scores across\nrepeatedtrainings.\nSpecifically, we trained the GAT+DGI model 10 times with different random seeds on three\nindependent datasets and compared the resulting attention scores across runs. As shown in\nResponse Fig. 5, the attention scores are highly reproducible, with pairwise correlations\nconsistently>0.95acrossruns.Inaddition,wequantifieddispersionusingthestandarddeviation\n(SD) and observed median SD < 0.02 (Response Fig. 9), indicating low variability. These results\ndemonstrate that aggregating attention scores over 10 runs provides a stable estimate in both\nsimulatedandreal-datasettings.\nWe have added these robustness statistics (correlation and CV) to the revised manuscript to\nguide users on appropriate parameter settings and to clarify that the default choice of 10 runs\nyieldsreliableandreproducibleattentionestimates.\nillustratingthedispersionofattentionscores,quantifiedbystandarddeviation(SD),asafunction\nof the number of independent model runs (from 2 to 10 runs). Results are shown for four\nrepresentativedatasets:IFN-β-treatedPBMCs,PANC-1cells,Bcells,andiAT2cells.\n4. scDNS computes both jointprobability densities p(Gᵢ, Gⱼ) and conditional densities p(Gᵢ| Gⱼ),\nclaimingeachhasdistinctadvantages.However,sinceconditionaldensitiesarederivedfromjoint\ndensities, they contain overlapping information. The authors should perform ablation studies\ncomparing performance using (a) joint density only, (b) conditional density only, and (c) both\ncombined. They should provide examples where each density type captures different biological\nsignals, and clarify the claim that each performs differently in \"high-density vs. low-density cell\npopulations\".\nResponse: We thank the reviewer for this thoughtful comment. The joint and conditional\nprobability densities capture different aspects of the bivariate distribution. When coupled with\nJSD, they can emphasize different types of distributional changes, and are therefore\ncomplementaryinpractice.\nConceptualrationale:\nThe joint density captures the overall bivariate distribution and global shifts of a gene pair. In\npractice, it is typically more tolerant to small local fluctuations and therefore tends to be more\nrobusttonoise,especiallyinsparseorlow-coveragesettings.Incontrast,theconditionaldensity\nemphasizesdirectional, localdependencystructure—i.e., howthedistribution ofG changesasa\ni\nfunction of G. This can be more sensitive to localized variation that occurs in a subset of cells\nj\n(e.g., specific subpopulations or state-restricted interactions). However, this increased locality\nandsensitivitycanalsomakeconditionaldensitiesmoresusceptibletolocalsamplingnoisewhen\nthe effective cell support is small. Thus, conditional densities can highlight dependency changes\nthat are weak in the global joint distribution, whereas the joint representation provides a\nstabilizing global view. We therefore combine them to mitigate each other’s limitations and\nimproverobustness.\nAblationanalysis:\nAsrequested,weperformedanablationstudyusingsimulatedperturbationsconstructedontop\nof a real scRNA-seq dataset. Specifically, we tested the predictive accuracy of three\nvariants—jointdensityonly,conditionaldensityonly,andtheirintegratedcombination—acrossa\ngradient of varying cell perturbation fractions. Our analysis revealed that when perturbations\naffect a large proportion of the cell population, all three variants recover the signal accurately.\nHowever,performancedifferencesemergewhenperturbationsareconfinedtosmallpopulations\n19\n(e.g., 10% of cells) (Response Fig. 10a). In such case, the conditional-density-only approach\nexhibits higher accuracy in detecting localized dependency changes but may also be susceptible\nto local noise. Crucially, the combined score exhibits superior overall predictive performance\nwhen aggregated across perturbation scales, while maintaining a favorable trade-off between\nsensitivityandrobustness(ResponseFig.10b).\nIllustrativeexamplesdemonstratingdistinctsignals:\nTo further clarify why the two density types can behave differently, we provide two\nrepresentativecases:Inoneexample,theglobaljointdistributionremainslargelystable,butthe\nrelationshipchangesinalocalizedregion(asubsetofcells),consistentwith interactionrewiring.\nHere, the joint-density signal is relatively muted, whereas the conditional density more clearly\ncapturesthedependencychange(ResponseFig.10c).Inasecondexample,thejointdistribution\nshifts globally, while the dependency structure between the two genes remains similar. In this\nscenario, conditional densities—being more focused on dependency—show a weaker change,\nwhereas the joint density appropriately reflects the global distributional shift (Response Fig.\n10d).\nOverall, these ablations and examples support our claim that joint- and\nconditional-density–based divergences have distinct practical strengths, and that combining\nthem yields a more accurate and robust divergence estimate across both abundant and rare\ncellularpopulations\nResponse Fig. 10 | Ablation study and illustrative examples of joint and conditional density\nscoring. a, Performance across varying perturbation fractions: Accuracy of perturbation gene\nprediction using joint density only (purple), conditional density only (yellow), and their\n20\ncombination (red) across a gradient of cell perturbation fractions (10% to 100%). b, Overall\npredictive performance when aggregated across perturbation scales. c, Illustrative example of\nlocalized dependency change: Scatter plots and density heatmaps showing a case where the\nglobaljointdistribution isstable(lowdivergence:0.04), buttheconditionalrelationshipchanges\nsignificantly in a subset of cells (higher divergence: 0.15), highlighting the sensitivity of\nconditional density to interaction rewiring. d, Illustrative example of global distribution shift: A\ncase where the joint distribution shifts globally (high divergence: 0.40), while the dependency\nstructure remains similar (lower conditional divergence: 0.05), demonstrating the robustness of\njointdensityincapturingglobaltranscriptionalshifts.\n5. scDNS constructs cell-type-specific networks, making its results inherently dependent on cell\ntype definitions. However, the robustness to cell type annotation choices is not evaluated. The\npaper uses Seurat standard clustering method, but the clustering could be very different based\non different parameters or methods being used. The authors should demonstrate the findings\nfromscDNSarereproducibleacrossreasonablevariations.\nResponse: We thank the reviewer for highlighting this important point. We agree that scDNS\nconstructs cell-type–specific networks and therefore depends on upstream cell-type definitions.\nIn our study, we used Seurat because it is a widely adopted and well-validated framework for\nsingle-cell clustering and annotation. While changing Seurat parameters—especially the\nclusteringresolution—canalterthenumberofclusters,stablemajorcellpopulationsaretypically\npreserved across a reasonable parameter range, and their marker genes remain largely\nconsistent.\nThatsaid,weacknowledgethatclustervariationcanaffecttheselectionofmarkergenesandthe\nconstruction of predefined positive edges. Importantly, as we demonstrated in our earlier\nresponses,boththeGAT+DGItrainingprocedureandthedownstreamscDNSscoringaretolerant\nto moderate perturbations in the predefined edge set. To directly evaluate robustness to\ncell-type annotation choices, we conducted an additional analysis on real data by systematically\nvaryingtheSeuratclusteringresolutiontogeneratealternativeclusteringwithdifferentnumbers\nof clusters. For each clustering, we re-derived marker genes, reconstructed the corresponding\ncell-type-specific networks with GAT+DGI, and quantified network similarity across resolutions.\nUsingaadjustedRand-based/graphsimilaritymeasure,weobservedconsistentlyhighagreement\nbetween networks inferred under different clustering resolutions (Response Fig. 11), indicating\nthat the inferred interaction structure and downstream scDNS signals are reproducible across\nreasonablevariationsincell-typeannotation.\nWe have added these robustness results to the revised manuscript to clarify that scDNS\nconclusions are not sensitive to modest changes in clustering parameters and cell-type\ndefinitions.\nof pairwise Spearman correlations of attention weights inferred across Seurat clustering\nresolutionsrangingfrom0.1to0.9intheIFN-β-treatedPBMCsdataset.Thenumberofclustersat\neachresolutionisshownontheleft.\n6.TheauthorsshouldreportthecomputationalruntimeandmemoryrequirementsofscDNSon\ndatasetsofvarying sizes. Assingle-celldatasetsincreasinglyexceedhundredsofthousandscells,\nusersneed to understand howruntime scales with the number of cells and genes, and whether\nthemethodisfeasibleforlarge-scaledatasets.\nResponse:Wethankthereviewerforthisrequest.InresponsetoReviewer#1’srelatedcomment\n(Minorcomment#1),wealreadyreported detailedruntimeandpeakmemory benchmarkson a\nlargedataset(>50,000cells),demonstratingthatscDNSisfeasibleforlarge-scaleanalyses.Briefly,\nwe profiled the major computational stages (preprocessing/network construction, embedding\npreparation, and GAT+DGI training) and observed that: (i) preprocessing and graph construction\nscale primarily with the expression matrix size, whereas (ii) the core GAT+DGI step operates on\nthe gene graph and therefore depends mainly on the number of genes and graph connectivity,\nand is relatively insensitive to cell count once embeddings are constructed. These benchmarks\nclarifyconfirmthatscDNSremainspracticalforlarge-scaledatasets.\nTofurtheraddressscalabilityacrossdatasetsizes,wehavenowbenchmarkedscDNSondatasets\nwithvaryingnumbersofcells/genesandreportbothwall-clockruntimeandpeakmemoryusage\n(Response Fig. 12 /Table 1). All experiments were performed on a standard high-performance\nworkstation(IntelXeonGold6226RCPU,768GBRAM,andasingleNVIDIARTXA6000GPU).\npipelineasafunctionofdatasetsize(numberofcells,log\n10\n-transformed).b,Peakmemoryusage\nduringscDNSanalysisasafunctionofdatasetsize.\n7.Theinterpretationofgenerankingsisinconsistentthroughoutthemanuscript.Forexample,\n•Line206:IRF4iscalledthe\"topperturbedregulator\"butranks273\n•Lines266-268:JAK1andWTAParehighlightedbuthavemoderateranks\nThe authors should clearly define what rank/p-value cutoff constitutes a \"top\" or \"significant\"\nsignal.Howshouldusersinterpretrankings(isrank273consideredhigh)?\nResponse: We thank the reviewer for noting this ambiguity. In scDNS, genes are ranked by the\ndegree-normalized perturbation statistic (scZ-score, with an associated p-value/FDR). Because\nthe rank is computed over all tested genes, its absolute value depends on the size of the\ncandidategeneuniverseandthestringencyoffiltering.Inaddition,scDNSisdesignedtocapture\nboth(i)geneswithsignificantexpressionchangesand(ii)geneswithmodestmRNAchangesbut\nstrong activity/relationship rewiring (DGR). As a result, a gene with a rank of ~200–300 can still\nrepresent a strong and statistically significant signal when the total number of tested genes is\nlargeandthegenepassesthesignificancethreshold.\nTo remove ambiguity, we have revised the manuscript to (i) avoid using the term “top” unless\nreferringtoanexplicitlydefinedsubset(e.g.,topN),and(ii)reportscZ-scoresandFDRsalongside\nrankwhenhighlightinggenes.\nRegarding the examples noted by the reviewer, we use volcano-style plots to jointly display\nexpression change and the scDNS perturbation score, and we specifically highlight genes with\nlimited expression change but strong scDNS-predicted perturbation—one of the key advantages\nof scDNS over expression-only analyses. For instance, IRF4 is highlighted because its inferred\nperturbation is strong despite minimal mRNA change; although its global rank is 273, it falls\nwithin the statistically significant set given the large gene universe considered. We therefore\nrevised thetexttorefertoIRF4asa“significantlyperturbed regulator”ratherthan ambiguously\n“top.” Similarly,JAK1 and WTAP are highlighted as illustrative DGR examples because they show\nstrongnetworkperturbationwithlimitedmRNAchange;wenowexplicitlyreporttheirscZ-scores\nand FDRs and clarify that they are discussed as biologically interpretable exemplars rather than\n23\nnecessarilythehighest-rankedgenes.\n8. For the COVID-19 analysis, scDNS identifies disease-driven perturbation signals, but the\nconsistencyof these signals acrosspatientsis notclear.If scDNS isapplied toindividual patients\nseparately, are the identified signals replicated across patients? What is the degree of\ninter-patient variability in perturbation signals? Is such variability technical, or driven by real\nbiological variability? It would be nice if the authors can show that scDNS can identify\npatient-specificdysregulationpatternsthatrelatetodiseaseseverityoroutcome.\nResponse:Wethankthereviewerforthisimportantquestion.Aswedetailedinourresponseto\nReviewer #1’s related comment (Comment #6), we performed additional analyses in the\nCOVID-19 cohort to assess cross-patient consistency and inter-patient variability. Briefly,\ndonor-level aggregation of scDNS scores shows that the major disease-driven perturbation\nsignals are replicated across most patients, rather than being driven by a small number of\nhigh-cell-count samples (Response Fig. 4b), and the scores are not materially explained by\ncommontechnicalcovariates(nCounts/nFeature;ResponseFig.4c–f).\nAt the same time, scDNS also reveals meaningful inter-patient heterogeneity. Notably, we\nobserved one control sample with an outlying perturbation profile (Response Fig. 4b). To\ndeterminewhether thisdeviationistechnicalorbiological,wecompareditsexpression program\nagainst other controls and found marked upregulation of genes involved in antigen processing\nand presentation, MHC complex assembly, and viral entry-related pathways, suggesting a\nplausible biological explanation (e.g., an unrecognized infection-related immune state) rather\nthanapurelytechnicalartifact(ResponseFig.13a-d).\nTo further assess the clinical relevance of inter-patient variability in perturbation signals, we\nintegratedscDNSscZ-scoreswithsingle-cellexpressionprofilesandperformedclusteringinaPCA\nembedding space. We then stratified cells into three response groups (low / medium / high)\nbasedonscZscorequantiles,representingincreasinglevelsoffunctionaldysregulationrelativeto\ncontrols(ResponseFig.13e&f).Notably,theproportionofhigh-scZscorecellsvariessubstantially\nacross COVID-19 patients, indicating heterogeneous enrichment of this perturbed cell state\nacrossthecohort(ResponseFig.13g).Toevaluatewhetherthisheterogeneityisbiologicallyand\nclinically meaningful, we derived a high-response gene signature by differential expression\nanalysis of the high-scZscore group and applied this signature to an independent bulk RNA-seq\ncohortofCOVID-19patientsandhealthycontrols.Thebulksignaturescorecanbeinterpretedas\nthe enrichment of the gene expression program of high-scZscore cells, and thus reflects the\nabundance and/or within-cell-type activation of high-response cells in each sample. Strikingly,\nthis score showsa strong positive association with clinical severity (mild → moderate → severe;\nR2 ≥ 0.92), supporting that scDNS-prioritized dysregulated cell programs capture clinically\nrelevant disease progression even though the underlying high-response cell fraction is not\nuniformly shared across all patients (Response Fig. 13h). In contrast, signatures derived from a\ndirectcomparison between all disease and control cells without scZscore-based stratification do\nnot reproduce this association with disease severity (Response Fig. 13i). Collectively, these\nresults support the biological and clinical relevance of scDNS-derived signals and argue against\npurelytechnicalartifacts.\nprogression. a-d, Characterization of inter-patient heterogeneity in control samples. Volcano\nplotsshowingdifferentialgeneexpressionanalysisbetweenoneoutlyingcontrolsample(C53ctr)\nand other control samples in AT2 cells (a) and AT1 cells (c). Key genes involved in immune\nresponse and viral entry are highlighted. Combined bar-dot plots illustrating Gene Ontology\nenrichmentanalysisoftheup-regulatedgenesinC53ctrAT2(b)andAT1(d)cells.e,UMAPplots\noffourrepresentativecelltypes(CD4+Tcells,CD8+Tcells,CyclingNK/Tcells,andDendriticcells)\ncolored byresponsegroups. Cellswere stratified into control (blue), low(light blue),mid (pink),\nand high (dark red) response states based on scZscore quantiles. f, Violin plots showing the\ndistribution of scZscores across the four response groups for each cell subpopulation. Black\nhorizontal linesindicate themean value. g,Bar plotsillustrating the proportion ofhigh-scZscore\ncells (red) versus low/mid-scZscore cells (blue) across individual COVID-19 patients, revealing\nsignificant inter-patient variability in the abundance of dysregulated cell states. h, Correlation\nbetween clinical severity (Control, Mild, Moderate, Severe) and the enrichment score of a\n\"high-response\" gene signature derived from high-scZscore cells. Data points represent\nindividual samples from an independent bulk RNA-seq cohort. i, Correlation between clinical\nseverityanda\"standard\"signaturederivedfromadirectcomparisonofallCOVID-19cellsversus\ncontrols without scZscore-based stratification. R2 (coefficient of determination) and P-values\n(two-sided)areshownforlinearregression.\n25\n9. Both the COVID-19 and IFN-β PBMC analyses utilized datasets with batch effects, which may\nintroduce technical variation. While the authors apply Harmony correction for UMAP\nvisualization, scDNS operates on uncorrected count data. Is scDNS robust under such\ncircumstances?\nResponse:Weagreethatbatcheffectsareanimportantconsiderationincohort-scalescRNA-seq.\nIn our analysis, Harmony is used only for embedding/UMAP visualization, consistent with\ncommon practice, because integration methods are designed to align batches in a latent space;\nusing batch-corrected representations as input for network/relationship inference can\nunintentionally alter gene–gene dependency structure. Since scDNS is designed to quantify\nrewiring of gene–gene dependencies, we therefore operate on uncorrected (but standardly\nnormalized)expressiontopreserveaninterpretablemutual-informationstructure.\nAs detailed in our response to Reviewer #1 (Comment #6), scDNS is robust to batch-related\ntechnicalvariationforboththeoreticalandempiricalreasons.First,commonbatchartifacts(e.g.,\nlibrary-depth shifts, global expression offsets, and increased measurement noise) primarily\nperturbmarginalexpressionlevelsandtendtodilutedependencyestimationratherthaninduce\ncoherent, directionally consistent rewiring across many genes. Second, empirically, in both the\nCOVID-19 and IFN-β PBMC cohorts, scDNS scores show no significant association with key\ntechnicalcovariatessuchastotalUMIcounts(nCounts)andnumberofdetectedgenes(nFeature),\nand the main disease/treatment signals remain consistent under donor-aware aggregation\n(ResponseFig.4c–f).Third,asanadditionalcontrolled stresstest,wesimulatedarepresentative\nbatcheffectdrivenbylibrary-sizescalingandfoundthatthecohort-levelmeandivergencescores\nshowed no substantial change, indicating that scDNS is not overly sensitive to\nlibrary-size–inducedtechnicalshifts(ResponseFig.14).\nFinally, scDNS produces cell-resolved perturbation readouts that can be aggregated by\nsample/batch and projected onto embeddings, providing an interpretable diagnostic view to\nexamine whether any signal is disproportionately driven by a specific batch. Collectively, these\nresults support that scDNS signals are not driven by batch structure despite operating on\nuncorrectedcounts.\nResponse Fig. 14| Robustness of network divergence estimates to batch-associated library size\neffects.Boxplotsshownetwork divergencescoresunder unperturbedand perturbed conditions\nacross increasing strengths of batch-associated library size effects, modeled by varying the\nstandarddeviationofalog-normaldistribution(0–0.2).Thehighlyconsistentdistributionsacross\nlibrarysizeshifts.\nI have successfully installed the scDNS R-package. I have tested the package using the \"ifnb\"\ndataset from SeuratData (the same ifnb-PBMC data used in this study), and the main functions\nworked well. The demo data provided by the authors on their GitHub page is not working\n(requirespermissiontoaccess)\nResponse:Weapologizeforthe inconvenienceregardingthe demo datalink. Thiswasduetoan\naccess-permission setting on the hosting location. We have now updated the repository and\nconfirmedthatthedemodatasetispubliclyaccessible.\n27"}, {"round": 2, "reviewer_comments": {"reviewer_1": "All my concerns have been addressed. I have no further comments", "reviewer_2": "I thank the authors for their thorough revision. I am pleased with the substantial additional analyses that address my concerns. The revision is overall strong, and I only have a few minor points for the authors to address.\nI have some follow-up questions related to comment 8:\n(1). The high-response gene signature in COVID-19 is very cool, but they are not listed anywhere. What are these genes? Any enriched pathways? Do they overlap with previously reported COVID severity signatures? Also it is a bit ambiguous to call them \"high-response\" genes cause a strong immune response is not equivalent to Covid19 disease severity.\n(2). The authors stated that an independent bulk RNA-seq cohort was used for validating the gene signature. But this cohort is not cited (either in rebuttal or Main text). Please specify the source of this cohort.\n(3). The interpretation of outlier C53ctr as having an 'unrecognized infection-related immune state' is plausible but speculative. Is clinical metadata available to support this?"}, "author_response": "Reviewer #2 (Remarks to the Author):\nI thank the authors for their thorough revision. I am pleased with the substantial\nadditional analyses that address my concerns. The revision is overall strong, and I\nonly have a few minor points for the authors to address.\nI have some follow-up questions related to comment 8:\n(1). The high-response gene signature in COVID-19 is very cool, but they are not\nlisted anywhere. What are these genes? Any enriched pathways? Do they overlap\nwith previously reported COVID severity signatures? Also it is a bit ambiguous to\ncall them \"high-response\" genes cause a strong immune response is not equivalent to\nCovid19 disease severity.\nResponse:\nWe thank the reviewer for this constructive feedback. We agree that the term “high-\nresponse” could be misinterpreted as implying a direct equivalence between the\nmagnitude of immune activation and COVID-19 severity. To avoid confusion, we have\nrevised the manuscript and now refer to these genes as “disease-driven perturbation\ngenes,” defined as genes showing the strongest scDNS-inferred network dysregulation\nwithin each immune lineage.\n(i) Gene lists.\nTo ensure full transparency, we now provide the complete lists of significant disease-\ndriven perturbation genes for each immune lineage (CD4⁺ T, CD8⁺ T, cycling NK/T-like,\nand dendritic cells) in Supplementary Table 2.\n(ii) Pathway enrichment.\nWe performed GO:BP and KEGG enrichment analyses and found that these lineage-\nspecific disease-driven genes are significantly enriched in immune processes highly\nrelevant to COVID-19 immunopathology (Response Fig. 1a; representative genes\nannotated). Across lineages, the top enriched themes consistently point to immune\nactivation and effector programs, including cytokine- and interferon-related signaling,\nregulation of T-cell activation, cytokine production (including TNF-related programs),\nand regulation of inflammatory responses. Collectively, these enrichments indicate a\ncoordinated immune activation/dysregulation program characterized by heightened\ninflammatory signaling and effector function.\nTo delineate shared core biology, we further derived a cross-lineage recurrent set\n(genes present in ≥2 lineages) and repeated the enrichment analysis, which yielded\ndominant terms related to leukocyte activation and positive regulation of\ninflammatory response. These results support the conclusion that scDNS captures\nseverity-associated state transitions in immune cells.\n(iii) Overlap with published COVID-19 severity signatures.\nWe note that a simple pooled, gene-by-gene overlap with previously reported\n“severity signatures” can be misleading, because many published signatures are\nderived from bulk/whole-blood analyses and can be strongly influenced by cell-\ncomposition shifts (often including epithelial- or neutrophil-associated transcripts),\nwhereas our signatures are immune-lineage- and cell-state-resolved. Therefore, we do\nnot necessarily expect large one-to-one overlap across studies; instead, we emphasize\nconcordance at the level of biological processes/modules, where we observe\nconsistent enrichment for immune/leukocyte activation (Response Fig. 1b)—an\noverarching theme widely linked to COVID-19 severity (PMID: 32788292, 32717743).\nNevertheless, as requested, we performed a direct overlap analysis using the union\ngene set across the four immune lineages and observed partial intersections with an\nestablished severe-versus-mild differential gene sets (31/1001, 3.1%) (PMID:\n33096026)(Response Fig. 1c). Notably, the overlapping genes prominently included\ninterferon-response transcripts such as CD38, IFI44L, HERC5, and SIGLEC1 (CD169),\nall\nof which are well-established type I interferon–inducible genes and have been\nconsistently associated with heightened inflammation and severe COVID-19 (PMID:\n33635833; 38772950; 41551039; 33825125; 36522333).\nThese revisions have been incorporated into the COVID-19 Results subsection (lines\n331–334) and Supplementary Note 3 (lines 613–614).\nResponse Fig. 1 | Enrichment analysis and validation of high-scZscore cell\nsignatures. a, Functional enrichment of upregulated genes in high-scZscore cells\nacross major immune compartments. Bar plots show significantly enriched Gene\nOntology Biological Process terms and KEGG pathways for CD4+ T cells, CD8+ T cells,\nDendritic cells, and Cycling NK/T cells. The x-axis represents the enrichment\nsignificance calculated as -log (adjusted P-values). Representative genes associated\n10\nwith each term are listed below the bars. b, GOBP and KEGG enrichment analysis of\ngenes commonly upregulated in high-scZscore cells across at least two cell types. c,\nVenn diagram illustrating the overlap between upregulated genes identified in high-\nscZscore cells (scDNS-derived) and those associated with COVID-19 severity in a\nbenchmark multi-omic study (Overmyer et al., 2021, PMID: 33096026). The\nintersection contains 31 common genes, representing a statistically significant overlap\n(P = 0.002, hypergeometric test), highlighting the biological relevance of the scDNS-\nidentified cell populations to clinical disease severity.\n(2). The authors stated that an independent bulk RNA-seq cohort was used for\nvalidating the gene signature. But this cohort is not cited (either in rebuttal or Main\ntext). Please specify the source of this cohort.\nResponse：\nWe thank the reviewer for highlighting this important point. The independent bulk\nRNA-seq cohort used to validate the high-response gene signature is GSE217948\n(Wang, Y. et al., Front. Immunol. 13:1043219, 2023). Although this dataset was already\ncited in the main text, we agree it was not sufficiently prominent. We have now added\nthe citation at the first mention of the bulk validation in the revised manuscript and\nalso included it in the relevant Data Availability section.\n(3). The interpretation of outlier C53ctr as having an 'unrecognized infection-related\nimmune state' is plausible but speculative. Is clinical metadata available to support\nthis?\nResponse:\nWe thank the reviewer for this important point. We agree that our previous\ninterpretation of C53ctr as having an “unrecognized infection-related immune state”\nwas overly specific given the absence of direct infection evidence. Upon re-\nexamination of the clinical annotations, donor C53ctr shows multiple comorbidities\nand respiratory manifestations, including pre-existing lung disease (e.g., asthma/COPD,\nrequiring home O₂, any chronic lung condition), pre-existing diabetes, prior smoking,\nand respiratory symptoms (see Supplementary Table S1 of the original study,\nhttps://static-content.springer.com/esm/art%3A10.1038%2Fs41586-021-03569-\n1/MediaObjects/41586_2021_3569_MOESM4_ESM.xlsx).\nWe have therefore revised the manuscript to describe C53ctr more conservatively as\nan immune-activated/inflammatory outlier. This is consistent with the elevated\nexpression of genes involved in antigen processing, MHC assembly, and related\nimmune programs in C53ctr compared with the other healthy controls\n(Supplementary Note Fig. 2g-j). Our donor-aware analyses confirm that this signal is\nnot explained by technical covariates (e.g., UMI counts, detected genes), supporting a\nbiological basis. As we demonstrated in last revision, this outlier illustrates donor-level\nheterogeneity captured by scDNS but does not alter our main cohort-level conclusions."}, {"round": 3, "reviewer_comments": {"reviewer_2": "I thank the authors for their additional revision. All my concerns have been addressed."}, "author_response": ""}]

Decision letter