Recovering gene regulatory networks in single-cell multi-omics data with PRISM-GRN

The framework of PRISM-GRN

As illustrated in Figure 1A, TFs bind to their specific target genes, thereby forming cell type–specific GRNs. In conjunction with TF expression levels and the chromatin accessibility of target genes, these GRNs drive the emergence of diverse cellular phenotypes through complex regulatory mechanisms. Motivated by such established principles of gene expression, we proposed PRISM-GRN, a probabilistic Bayesian model for reconstructing cell type–specific GRNs based on single-cell multi-omics. Based on the aforementioned gene expression principle, PRISM-GRN incorporates a biologically interpretable architecture into the probabilistic framework, which is shown in Figure 1B. By integrating gene expression profiles, chromatin accessibility data, and prior GRNs, PRISM-GRN decomposes the observable expression of each gene into latent variables that capture the unobserved contributions of its regulatory factors. This is achieved through a regulation mechanism–informed variational inference process, enabling the precise reconstruction of cell type–specific GRNs.

Figure 1.

The framework of PRISM-GRN. (A) The motivation of PRISM-GRN. Within the same cell type, a gene's expression is influenced by the expression levels of the TF targeting it and the accessibility of its chromatin through the cell type–specific GRN. This regulatory mechanism shapes distinct gene expression patterns, ultimately defining the cellular phenotype with different states and functional identities. (B) The variational inference process. PRISM-GRN models the generation of gene expression as a Dirichlet-multinomial distribution, integrating three key latent variables: expression-related factor z_Exp, chromatin accessibility-related factor z_Acc, and regulatory interaction–related factor z_GRN. These latent variables are inferred from observed data and structured by biological mechanisms. The generative process reconstructs expression levels and GRNs based on these factors, ensuring PRISM-GRN adheres to regulatory principles. (C) Performance evaluations of PRISM-GRN. PRISM-GRN's performance on recovering GRNs and capturing causality was comprehensively evaluated and compared with the baseline methods. The robustness of PRISM-GRN was further validated by its strong performance with limited prior knowledge and unpaired omics data. Additionally, downstream biological analyses were conducted to highlight PRISM-GRN's real-world efficacy in improving the understanding of gene functions and finding GRNs.

Open in new tabDownload PowerPointLink to figure

More specifically, for a given cell type–specific data set, the input of PRISM-GRN includes the gene expression profile E ∈ ℝ^{M × N} from scRNA-seq; the chromatin accessibility score profile R ∈ ℝ^{M × N} from scATAC-seq, where M is the number of genes and N is the number of cells; and a binary adjacency matrix A ∈ {0, 1}^{M × M} representing a cell type–specific prior GRN G sourced from a public database. The adjacency matrix of the prior GRN, where A_ij = 1 indicates the presence of a known regulatory interaction between gene i and gene j, and A_ij = 0 indicates an unknown regulatory interaction, also acts as the prior knowledge to inform the prior distributions of GRN connectivity across genes during the inference process. More details of preprocessing procedures are provided in Supplemental Content S1. The output of PRISM-GRN is a reconstructed GRN G′, for which the unknown regulatory interactions between the given genes are predicted.

The mechanism-informed variational process ensures the model's flexibility in capturing complex regulatory interactions while adhering to biological principles. Constrained and guided by the biological mechanism–informed generative process, PRISM-GRN can achieve greater precision and interpretability in reconstructing GRNs. As shown in Figure 1C, comprehensive benchmark experiments were conducted to demonstrate the efficacy of PRISM-GRN in reconstructing GRNs and predicting causality. Experiments on data sets with unpaired omics data and limited prior GRN information were performed to demonstrate PRISM-GRN's adaptability across diverse real-world biological contexts. Finally, practical applications and biological analyses reveal that PRISM-GRN facilitates the identification of cell type–specific or context-specific regulatory relationships and core regulators under different biological research scenarios.

Baselines

In this paper, we proposed PRISM-GRN, a probabilistic Bayesian model designed to recover cell type–specific GRNs. To evaluate the effectiveness of PRISM-GRN in GRN recovery, we benchmarked it against seven baseline methods, including two unsupervised methods (GENIE3 [Huynh-Thu et al. 2010] and GRNBoost2 [Moerman et al. 2019]), two supervised methods (GeneLink [Chen and Liu 2022] and Grace [Wang et al. 2024a]), and three multi-omics methods (scAI [Jin et al. 2020], CellOracle [Kamimoto et al. 2023], and SupirFactor [Tjärnberg et al. 2024]). Among these methods, GENIE3 and GRNBoost2 rely solely on scRNA-seq data, as do GeneLink and Grace, although the latter two require prior GRNs as part of their supervised training process. In contrast, scAI, CellOracle, SupirFactor, and PRISM-GRN integrate both scRNA-seq and scATAC-seq data for GRN inference. Notably, SupirFactor, as an unsupervised method, also requires the same prior GRNs as architectural constraints to guide model structure to get the explained relative variance. The details of the baseline methods are shown in Supplemental Table S1. All baseline methods were implemented using default configurations to ensure consistency and reproducibility, whereas the default hyperparameters of PRISM-GRN were selected using a grid search and are presented in Supplemental Table S2.

Data set description and experimental designs

We used four data sets with paired scRNA-seq and scATAC-seq data to benchmark PRISM-GRN against seven baseline methods, evaluating their performance in GRN reconstruction and causality detection. These data sets include the cultured human cell lines A549 (Cao et al. 2018), GM12878 (Ma et al. 2020), and K562 (Ma et al. 2020), as well as the mouse cell line 3T3 (Ma et al. 2020) with homogeneous cells. Four data sets with unpaired scRNA-seq and scATAC-seq data, including two benchmark data sets, A549 and GM12878, and two additional data sets, human embryonic stem cells (hESCs) (Chu et al. 2016; Samara et al. 2022) and mouse embryonic stem cells (mESCs) (Hayashi et al. 2018; Zamanighomi et al. 2018), were used to assess PRISM-GRN's capability to recover GRNs with incomplete data modalities. The detailed information of the benchmark data sets is summarized in Table 1, including the cell number, TF number, gene number, ground-truth GRN number, and reference genome, as well as the data source. Notably, although the hg19 (GRCh37) reference genome was used in this study, the conclusions are not significantly affected by the choice of assembly, such as hg38 (GRCh38) or a more recent genome assembly, because our analyses depend on gene-level features rather than assembly-specific coordinates, and the benchmarking data sets include genes that are consistently annotated across both assemblies.

For each benchmark data set, cell type–specific prior GRNs were derived from ChIP-seq data archived in the ChIP-Atlas database (Oki et al. 2018), which is a large-scale integrative platform that aggregates and curates publicly available ChIP-seq and related epigenomic data sets. ChIP-Atlas annotates experiments with detailed metadata, including cell type or cell line, and aggregates peak-calling results to offer cell type–specific binding profiles, thereby providing a reliable ground truth for evaluating GRN inference methods. The ground-truth TF–gene regulatory interactions were divided into training and test sets as the positive samples. Specifically, 90% of the interactions were used for training to provide prior knowledge, whereas the remaining 10% were held out as positive gold standards for evaluation. To incorporate high-confidence negative samples, following the method of Zou et al. (Wang et al. 2024a), for each transcription factor TF_i and all its known target genes Gene_j, unobserved regulatory pairs (TF_i, Gene_NT), where Gene_NT denotes genes not targeted by TF_i, were sampled as negative samples and split between training and test sets at the same 9:1 ratio. Additionally, to mitigate overfitting and optimize model performance, one-ninth of the training set was randomly selected as a validation set.

PRISM-GRN, along with supervised baseline methods including GeneLink and Grace, was trained on the same training set and evaluated on the corresponding test set. In contrast, unsupervised methods were directly evaluated on the test set for each data set. As shown in Table 1, with the relatively sparse known regulatory interactions, GRN reconstruction is inherently a highly imbalanced classification task; therefore, multiple evaluation metrics were employed. The area under the receiver operating characteristic curve (AUROC) and the area under the precision-recall curve (AUPRC) were used as primary metrics to assess model performance. Additionally, accuracy, precision, recall, and F1-score were included to provide a more comprehensive evaluation. Detailed descriptions of these metrics are provided in Supplemental Content S2.

PRISM-GRN outperforms baseline methods in GRN reconstruction

We first evaluated the performance of PRISM-GRN in reconstructing undirected GRNs, a task naturally formulated as a highly imbalanced binary classification problem, in which the number of noninteracting gene pairs (negative class) is roughly 100 times greater than that of true regulatory interactions (positive class). Specifically, we compared PRISM-GRN against seven baseline methods on four benchmark data sets with paired scRNA-seq and scATAC-seq data. Following the aforementioned training strategy, the supervised methods, including PRISM-GRN, Grace, and GeneLink, were trained on the same training set, whereas the unsupervised methods, including GENIE3, GRNBoost2, scAI, CellOracle, and SupirFactor, were directly evaluated on the same test set without a training phase.

Given the extreme class imbalance, we prioritized the AUPRC as the primary evaluation metric, as it more accurately reflects performance in identifying rare but meaningful regulatory interactions. As illustrated in Figure 2, PRISM-GRN consistently outperforms all baseline methods across all four benchmark data sets, achieving a mean AUPRC of 0.476. This is a substantial improvement over the supervised methods, including the second-best method, GeneLink, which achieves a mean AUPRC of 0.322, and the third-best, Grace, at 0.263. In contrast, the unsupervised methods, especially those relying on multi-omics input, are more vulnerable to noise and low sequencing quality. For instance, both scAI and CellOracle achieve an AUPRC of approximately 0.004 on the 3T3 data set. The superiority on AUPRC highlights PRISM-GRN's capability to recover true regulatory relationships even when they are vastly outnumbered by negatives. This capability is critical in real-world GRN inference tasks, in which identifying a small number of biologically meaningful interactions from a vast candidate space is the norm.

Figure 2.

The performance comparison of PRISM-GRN with the baseline methods on the benchmark data sets in GRN reconstruction. The AUROC and AUPRC metrics are used to evaluate the performance, and the results demonstrate that PRISM-GRN outperforms the baseline methods across all data sets.

Open in new tabDownload PowerPointLink to figure

For completeness, we also report AUROC results. All supervised methods achieve high AUROC scores above 0.9, with PRISM-GRN again leading at a mean AUROC of 0.984, followed by GeneLink, which achieves 0.970, and Grace, which achieves 0.964. Similar to the results shown by AUPRC, the unsupervised methods show relatively poor performance. This can be attributed to their sensitivity to weak or noisy expression and chromatin accessibility values, which leads to difficulty in separating true regulatory relationships from noise.

Detailed metrics, including accuracy, precision, recall, and F1, were used to more comprehensively evaluate the performance of PRISM-GRN and the baseline methods. As listed in Supplemental Table S3, the results further highlight the effectiveness of PRISM-GRN in dealing with class imbalance. PRISM-GRN's advantage in AUPRC is largely driven by its high precision, with a mean precision of 0.634, which reflects a low false-positive rate. This is critical in biological validation, in which experimental resources are limited, and only high-confidence predictions are actionable. Although GeneLink and Grace maintain competitive AUROC scores and achieve higher recall scores, their precision on rare positive events is significantly lower, leading to a drop in AUPRC and suggesting a tendency to generate more false positives.

We adopted a commonly used distance threshold of ±2000 bp from the transcription start site (TSS) to construct the ground-truth GRNs in cell type–specific benchmark data sets. To assess the sensitivity of PRISM-GRN to this heuristic, we additionally constructed alternative ground-truth GRNs using stricter proximity thresholds of ±1500 bp, ±1000 bp, and ±500 bp. We then evaluated PRISM-GRN alongside seven baseline methods on the GRN reconstruction task, using each of these ground-truth variations. As shown in Supplemental Figure S1, PRISM-GRN consistently outperforms all baseline methods under both AUPRC and AUROC metrics across all data sets and threshold settings.

These results collectively demonstrate that PRISM-GRN not only outperforms all baseline methods across both AUPRC and AUROC but also shows remarkable robustness to data imbalance, species variation, and gene scale differences. Its biologically informed variational inference framework allows PRISM-GRN to effectively capture rare but functionally critical regulatory interactions, making it a powerful and practical tool for real-world GRN reconstruction from multi-omics single-cell data. The relatively poor performance of the unsupervised methods observed here, especially those based on multi-omics data, may be attributed to their sensitivity to input data quality and the lack of task-specific supervision; we further analyze their limitations and discuss their suitable application scenarios in the Discussion section.

PRISM-GRN captures causality in gene regulation

The variational inference used in PRISM-GRN can capture the GRN associations informed by the gene regulation mechanism; therefore, PRISM-GRN has enhanced ability in causality inference among genes. Here, we further evaluated the causality prediction performance of PRISM-GRN with the two supervised methods, Grace and GeneLink, which can also predict causality through training, on the four benchmark data sets with directional GRN ground truth. Specifically, for genes i and j, there are three causal probabilities: y = 1 if gene i regulates gene j; y = 2 if gene j regulates gene i; and y = 0 if no relationship exists between the genes. By using variational inference, we can easily adapt PRISM-GRN for a triple classification task by fixing the categorical distribution output number in the GRN decoder. For Grace and GeneLink, following the method of Zou et al. (Wang et al. 2024a), we configure their output layers to consist of three neurons, with a Softmax layer to produce probabilistic predictions.

As illustrated in Figure 3A, PRISM-GRN consistently outperforms both Grace and GeneLink across all evaluation metrics except recall. Notably, although PRISM-GRN lags in recall, it exhibits superior performance in AUROC and AUPRC, both of which are crucial indicators of overall predictive accuracy. Specifically, PRISM-GRN achieves an average AUROC of 0.988 and an average AUPRC of 0.730 across the four benchmark data sets. In comparison, GeneLink obtains an average AUROC of 0.952 and AUPRC of 0.528, whereas Grace lags further behind with an AUROC of 0.856 and AUPRC of 0.467. PRISM-GRN improves AUPRC by 38.3% over GeneLink and by 56.3% over Grace, demonstrating its superior ability to rank true regulatory interactions more effectively. Furthermore, PRISM-GRN achieves an average precision of 0.738, significantly outperforming GeneLink (0.446) and Grace (0.445), with an improvement of 65.5% and 65.9%, respectively. This substantial increase in precision further highlights PRISM-GRN's capability to accurately identify true regulatory interactions. Meanwhile, the performance metrics that were calculated using the weighted-average are presented in Supplemental Table S4, PRISM-GRN still outperforms all baseline methods across the data sets.

Figure 3.

Performance evaluations on the causality prediction. (A) Performance comparisons of PRISM-GRN with the supervised methods, Grace and GeneLink, on the benchmark data sets, using comprehensive metrics including accuracy, precision, recall, F1, AUROC, and AUPRC. (B) The receiver operating characteristic (ROC) curves and precision-recall (PRC) curves for each class on the A549 data set, illustrating the detailed performance of PRISM-GRN and the baselines in causality prediction.

Open in new tabDownload PowerPointLink to figure

We further plotted the ROC and PRC curves for each data set and class individually to evaluate the ability to predict causality. As demonstrated in Figure 3B and Supplemental Figure S2, PRISM-GRN consistently outperforms the baseline methods in terms of AUC across all classes, demonstrating its robust performance in causal inference. Given the highly imbalanced nature of GRNs, in which the majority of gene pairs belong to the noninteracting 0 class, PRISM-GRN exhibits superior performance in precision-recall curves, particularly in smaller subclasses that contain regulatory interactions. Although Grace and GeneLink exhibit relatively low AUPRC values in these sparse gene classes, indicating their limited ability to identify true causal interactions, PRISM-GRN achieves AUPRC values of approximately 0.99 for larger, more noninteractions and 0.6 for smaller and sparser causal interactions.

We also benchmarked the causality prediction results of PRISM-GRN with the unsupervised method SCENIC+ (Bravo González-Blas et al. 2023). SCENIC+ constructs initial GRNs using coexpression-based algorithms such as GENIE3 or GRNBoost2 (which are also baselines) and then refines the network by integrating motif enrichment to enforce biologically informed constraints. However, instead of inferring regulatory directionality directly from data, it relies on a user-specified list of TFs to assign direction, which limits its utility for causal inference, especially when the TFs are unknown. In contrast, PRISM-GRN learns regulatory directionality in a data-driven manner without requiring prior specification of TFs.

To ensure a fair comparison between PRISM-GRN and SCENIC+, we constructed a subset of the test set consisting only of TF-to-target gene relationships. This subset matches the assumptions of SCENIC+, which treats all TFs as causal regulators, allowing us to assign directionality to its predictions by default. This benchmarking strategy isolates the direction-specific predictive power of each method and ensures that all models are evaluated under comparable conditions. As demonstrated in Supplemental Table S5, PRISM-GRN consistently outperforms both SCENIC+ and its base methods, GENIE3 and GRNBoost2, in recovering directed TF–gene interactions across all data sets and performance metrics, including AUROC, AUPRC, precision, recall, and F1-score. Although SCENIC+ achieves higher precision after motif-based filtering, the recall remains extremely limited, as it heavily relies on a predefined list of TFs and the presence of known motifs. In contrast, PRISM-GRN does not require known TF lists or motif annotations, highlighting its strong discriminative power and practical utility in causal GRN inference.

We further use CRISPR data sets to more comprehensively evaluate the predicted causality results of PRISM-GRN. Specifically, we utilized the Cancer Dependency Map (DepMap) data repository (Arafeh et al. 2025), which archives genome-wide CRISPR knockout screens across a wide range of human cancer cell lines. For each TF, DepMap identified its most significantly affected downstream genes based on gene expression changes upon TF knockout. Given that the CRISPR data sets and the benchmark data sets have different TFs and genes, we performed causality prediction on four data sets and obtained predicted TF–gene relationships on their respective test sets; the intersections of the TF–gene pairs present in both data sets were further selected for validation. As shown in Supplemental Table S6, among the overlapping gene pairs between the CRISPR and benchmark test sets, PRISM-GRN was able to recover most of the validated regulatory edges, further supporting the biological relevance and causal validity of our predictions. Also, the crossing validation using CRISPR data sets can further provide orthogonal evidence for PRISM-GRN's predicted causal edges.

Overall, these findings demonstrate that PRISM-GRN outperforms the supervised methods, Grace and GeneLink, and the unsupervised method, SCENIC+, in causal inference tasks. By leveraging variational inference, PRISM-GRN can capture hidden and less obvious causal relationships within sparse and complex GRNs. More importantly, in GRNs, many TFs simultaneously regulate other genes and are themselves regulated. PRISM-GRN's high performance in causal prediction enables it to identify potential directed regulatory relationships, even in the absence of prior specification of the TFs.

PRISM-GRN integrates unpaired omics data to reconstruct GRNs

PRISM-GRN utilizes a probabilistic framework for the reconstruction of cell type–specific GRNs, effectively integrating both scRNA-seq and scATAC-seq data, along with prior biological knowledge. By leveraging variational inference to approximate the posterior distribution of latent variables, PRISM-GRN effectively captures the underlying biological signals amid noise, ensuring robust integration of transcriptomic and epigenomic data. Notably, the mechanism-driven variational process decomposes gene expression into contributions from TFs that regulate it, the chromatin accessibility of its regulatory regions, and the TF–gene interactions. This structure allows PRISM-GRN to independently use scRNA-seq and scATAC-seq data, enhancing the utility of PRISM-GRN in a wide range of experimental settings, even in the absence of paired omics data, which is a common issue owing to experimental constraints.

To validate PRISM-GRN's ability to integrate unpaired omics data, we evaluated it on two benchmark data sets, GM12878 and A549, each containing independently obtained scATAC-seq data. Specifically, for GM12878, a replicate of the paired scATAC-seq data set from the same sequencing batch was used as the unpaired input. For A549, scATAC-seq data were derived from an entirely separate sequencing experiment. Additionally, we incorporated two widely studied cell lines, hESCs and mESCs, which are frequently used for GRN inference owing to the availability of curated prior networks. For these data sets, we utilized the processed scRNA-seq and prior GRNs from BEELINE (Pratapa et al. 2020), pairing them with independently sequenced scATAC-seq data from the same cell lines. In all cases, the scATAC-seq data were transformed into chromatin accessibility scores following a unified pipeline (see Methods) (Supplemental Content S2). Performance was evaluated using precision, AUROC, and AUPRC, in comparison with Grace and GeneLink, two supervised methods that have demonstrated secondary performance in benchmarks.

PRISM-GRN was evaluated in the tasks of GRN reconstruction (Fig. 4A) and causality prediction (Fig. 4B), and it consistently maintains stable and reliable performance, even in the absence of paired scRNA-seq and scATAC-seq data. In the A549 and GM12878 data sets, PRISM-GRN exhibits minimal performance degradation when using paired data sets for both GRN reconstruction and causality prediction tasks. Notably, in the GM12878 data set, in which the scATAC-seq data were derived from the same sequencing batch and included a greater number of cells, PRISM-GRN demonstrates superior performance when utilizing unpaired scATAC-seq data compared with paired omics data. Specifically, for the GRN reconstruction task, PRISM-GRN achieves an AUPRC of 0.495 with unpaired scATAC-seq data, outperforming the AUPRC of 0.465 obtained with paired data. Similarly, in causality prediction, PRISM-GRN achieves an AUROC of 0.780 with unpaired data, surpassing the AUROC of 0.758 achieved with paired data sets. More importantly, PRISM-GRN achieves a comparable performance with Grace and GeneLink, on AUROC and AUPRC across all four data sets, and still significantly outperforms them in precision when using the unpaired omics data, showcasing its superior capability in stably reconstructing casual GRNs.

Figure 4.

Radar plots of performance comparisons on unpaired omics data sets using precision, AUROC, and AUPRC. (A) For GRN reconstruction, the radar plot shows the performance comparisons of PRISM-GRN using unpaired data on the A549 and GM12878 data sets, compared with results from paired data, as well as with Grace and GeneLink. On the hESC and mESC data sets that with no paired data, PRISM-GRN is compared with Grace and GeneLink. (B) For causality prediction, the radar plot illustrates the performance comparison of PRISM-GRN, Grace, and GeneLink on the same data sets. PRISM-GRN can achieve comparable results to those obtained with paired data and those obtained by Grace and GeneLink.

Open in new tabDownload PowerPointLink to figure

When integrating unpaired data sets, a key challenge is the batch effect introduced by unpaired scATAC-seq data. PRISM-GRN does not directly operate on raw peak-resolution scATAC-seq data. As described in the preprocessing section, we preprocess the ATAC data using Signac, which converts the high-resolution peak matrix into gene-level gene accessibility scores. Although this transformation would sacrifice some fine-scale chromatin accessibility information, it serves two critical purposes: (1) it enables direct alignment with scRNA-seq gene expression data, and (2) it reduces batch-specific artifacts arising from peak-calling inconsistencies and sequencing depth variations. By aggregating signals at the gene level, this preprocessing acts as a normalization step that smooths out noise and mitigates variability across batches. We further illustrate the impact of batch effects in Supplemental Content S3 and Supplemental Figure S3.

Overall, the Bayesian probabilistic framework allows PRISM-GRN to effectively handle unpaired omics data, providing a more comprehensive understanding of GRNs. This capability, combined with its robustness and accuracy in reconstructing cell type–specific regulatory networks, posits PRISM-GRN as a powerful tool for advancing the study of gene regulation. Its flexibility in integrating unpaired data without compromising performance makes it particularly valuable in scenarios in which paired data sets are difficult to obtain, enabling broader applications in genomics and single-cell research.

PRISM-GRN achieves robust performance with limited prior knowledge

For each cell type–specific data set, PRISM-GRN uses the prior knowledge, that is, the ChIP-seq-derived prior GRNs, as the GAT skeleton and positive samples to inform the prior distributions of GRN connectivity across genes during the variational inference. Therefore, we additionally evaluated the impact of prior knowledge on the performance of PRISM-GRN. In the previous benchmarking experiments, the known GRNs were split into training and test sets with a 9:1 ratio, with the dominant portion being used for training. In this section, we investigated how PRISM-GRN performs when the proportion of prior knowledge in the training set is gradually reduced, while increasing the proportion of the test set. We progressively decreased the fraction of prior knowledge by 5% intervals, starting from 90% down to 50%. Given the results from previous experiments, in which PRISM-GRN demonstrates strong performance in terms of precision, we focused on evaluating its sensitivity to the amount of prior knowledge using precision, AUROC, and AUPRC as key metrics.

As shown in Figure 5, A and B, PRISM-GRN demonstrates stable performance across the metrics in the task of GRN reconstruction and causality detection, despite the reduction in prior GRN knowledge and the increase in test set size. Specifically, as the proportion of prior knowledge decreases, PRISM-GRN maintains consistently outstanding performance, with AUROC remaining above 0.99 in both tasks. More importantly, when the imbalance in the data set becomes more pronounced, PRISM-GRN still achieves robust predictive performance. In GRN reconstruction, PRISM-GRN achieves an average AUPRC of 0.469 across all the proportions and data sets, whereas in causality prediction, it reaches an AUPRC of 0.764. Notably, the precision consistently exceeds 0.6 across different data sets, further emphasizing PRISM-GRN's stability in maintaining a low false-positive rate in highly imbalanced sets.

Figure 5.

Performance evaluations with decreased prior knowledge. (A) Line plots showing the changes in precision, AUROC, and AUPRC across different data sets as the proportion of prior knowledge decreases from 0.9 to 0.5 with intervals of 0.05. PRISM-GRN maintains stable performance across all metrics during GRN reconstruction. (B) Line plots illustrating the changes in performance metrics for causality prediction. Similarly, PRISM-GRN maintains its performance even as prior knowledge decreases.

Open in new tabDownload PowerPointLink to figure

These results suggest that PRISM-GRN is not overly reliant on prior knowledge and can effectively identify regulatory interactions, making it a reliable tool for GRN reconstruction and causality inference, particularly in scenarios in which limited prior knowledge is available or in novel biological contexts in which prior data may be sparse. The ability of PRISM-GRN to effectively adapt to varying levels of prior knowledge ensures its applicability to a wide range of data sets and biological research tasks, further strengthening its potential for exploring gene regulatory mechanisms.

PRISM-GRN improves the understanding of gene functions

In PRISM-GRN, the latent variable z_GRN is utilized to generate gene regulatory relationships through a categorical distribution, thereby refining the causal effects among genes. To further assess PRISM-GRN's capabilities, we analyzed the A549 data set by clustering genes using the Leiden algorithm (Traag et al. 2019) and visualizing the clusters with UMAP, based on both gene expression data and the embeddings of z_GRN derived from the causal prediction task. As shown in Figure 6A, clustering based solely on gene expression results in limited distinction among genes. In contrast, clustering results using the trained z_GRN embeddings reveals clearer heterogeneity among genes, indicating that these embeddings encapsulate more biological information (Fig. 6B).

Figure 6.

Illustration of PRISM-GRN to improve understanding of gene functions. (A) UMAP plot of genes in the A549 data set based on the raw gene expression, colored by the Leiden clustering results. (B) UMAP plot of genes in the A549 data set based on the trained latent embeddings (z_GRN), colored by the Leiden clustering results. (C) KEGG enrichment analysis results using all genes in the A549 data set. (D) KEGG enrichment analysis results using genes in cluster 0, highlighting more significant pathway terms compared with the analysis based on all genes.

Open in new tabDownload PowerPointLink to figure

Considering that the analyzed genes are HVGs (Supplemental Content S1), we proceeded to perform KEGG and GO enrichment analyses on the entire gene set and on genes from cluster 0. As demonstrated in Figure 6, C and D, PRISM-GRN facilitated the identification of more significant KEGG enrichment results compared with analyses conducted directly on the highly variable genes. Notably, in the analysis of cluster 0, the most significantly enriched pathway identified was “central carbon metabolism in cancer,” which was not as prominent in the analysis of the entire gene set. A549 cells, as non-small-cell lung cancer cells, have been found to exhibit typical cancer metabolic features, including enhanced glycolysis and lactate fermentation, maintaining high levels of lactate production even under aerobic conditions (Liu et al. 2014; Sha et al. 2021). For the second-rank significantly enriched pathway, “circadian rhythm,” a pathway that is insignificant in the enrichment results based on the entire gene set, A549 cells have been proved to exhibit dysregulated expression of circadian rhythm–related genes such as TYMS, LGR4. These genes play crucial roles in regulating cellular processes like proliferation, migration, and invasion, suggesting that the “circadian rhythm” pathway is a key signaling pathway highly relevant to the gene functions of A549 cells (Xiang et al. 2018, Sun et al. 2024). Similarly, as shown in Supplemental Figure S4, the GO enrichment analysis of the specific gene sets in cluster 0 reveals additional functional terms, such as “regulation of pri-miRNA transcription by RNA polymerase II” and “epithelial cell development,” both of which are critically involved in the regulatory functions of A549 cells (Foster et al. 1998; Lin and Gregory 2015). The GO enrichment analysis results for genes from cluster 1 to cluster 9 are also illustrated in Supplemental Figure S5.

Beyond the analysis of latent variable, we also cross-validated the predicted directed interactions with those recorded in the TRRUST database (Han et al. 2018). To ensure consistency, we limited our analysis to the intersection of genes present in both data sets. As shown in Supplemental Table S7, PRISM-GRN successfully recovers nearly all TF–gene regulatory relationships overlapping with TRRUST (none of which were seen during training). Only two regulatory interactions were missed, highlighting PRISM-GRN's strong generalization ability and its robustness in recovering biologically meaningful, novel interactions.

These results demonstrate PRISM-GRN's effectiveness in capturing and utilizing latent biological information to enhance the interpretation of gene functions and regulatory mechanisms. Using the probabilistic framework and integrating prior knowledge, PRISM-GRN provides a more comprehensive understanding of gene interactions and their functions in complex biological processes.

PRISM-GRN finds cell type–specific gene regulators

We applied PRISM-GRN to reconstruct cell type–specific GRNs from the data set of human peripheral blood mononuclear cells (PBMCs) with paired scRNA-seq and scATAC-seq data, derived from 10x Genomics, in which the cells were annotated by Cao and Gao (2022). To ensure robust network reconstruction, we focused on the five cell types with the highest number of cells (CD4 naive T cells, CD8 naive T cells, CD14 monocytes, natural killer cells, and naive B cells) and incorporated cell type–specific prior knowledge (Fig. 7A). Notably, to assess the practical ability of PRISM-GRN in detecting cell type–specific regulators, we utilized all known GRNs derived from ChIP-seq data as prior knowledge and trained PRISM-GRN in the GRN reconstruction task. PRISM-GRN was then applied to complete GRNs by predicting the unknown regulatory interactions among genes.

Figure 7.

Illustration of PRISM-GRN for finding gene regulators in cell type–specific GRNs. (A) UMAP plot of the top five cell types with the highest number in the PBMC data set based on the raw gene expressions. (B) Predicted GRN for CD4 naive cells, displaying the top 25 genes with the highest degree. The node color represents the betweenness centrality of genes in the predicted GRN, and node size reflects the degree of each gene. The color and size of the edges correspond to the probability assigned by PRISM-GRN in the predicted GRN. (C) GO enrichment analysis results for the 25 genes in CD4 naive cells. The FDR indicates the significance of enrichment, with P-values corrected for multiple testing using the Benjamini–Hochberg procedure. (Signal) Weighted harmonic mean of the observed/expected ratio and −log (FDR), (gene count) the number of proteins annotated with a particular term. The results highlight the genes as significant regulators associated with CD4 naive cells.

Open in new tabDownload PowerPointLink to figure

In the case of CD4 naive T cells, a subset of T lymphocytes that play a crucial role in the adaptive immune response by serving as precursors to various effector and memory T cells upon activation (Zhu and Paul 2008), the top 25 genes with the highest degree in the predicted GRN are shown in Figure 7B, alongside the predicted regulatory network. To further identify key regulatory genes and understand their functional significance, we additionally calculated betweenness centrality, a critical metric that quantifies a node's importance in a network by measuring how frequently it appears on the shortest paths between other nodes. The top 25 genes, all of which exhibit significant betweenness centrality according to the density plot (with a one-sided hypothesis t-test at the 0.05 level) (Supplemental Fig. S6), suggest their potential roles as key regulators in CD4 naive cells. Furthermore, we calculated both degree and betweenness centrality scores for these genes in the prior GRN, as listed in Supplemental Table S8. Notably, their rankings in the predicted GRN differ from those in the prior GRN, suggesting that PRISM-GRN does not solely depend on prior knowledge for GRN inference. Instead, PRISM-GRN effectively integrates multi-omics data and employs probabilistic modeling to infer regulatory interactions beyond the prior knowledge. As illustrated in Figure 7C, the Gene Ontology enrichment analysis based on the top 25 genes reveals a significant functional correlation with the CD4 cells.

More specifically, among these top 25 genes, 20 genes have been experimentally validated as TFs involved in regulatory functions within CD4 naive cells, even though they were not explicitly designated as TFs during the training phase. Among the remaining six genes, ETV3L, a member of the ETS TF family, has not yet been explicitly characterized in the context of CD4 naive cell regulation. However, given that the ETS family has been well established to play critical roles in immune function (Hsing et al. 2020; Qu et al. 2020), ETV3L may possess corresponding regulator potential in CD4 naive cells. Although HLA-DRA is not a TF, it plays a pivotal role in the immune system. It encodes the alpha chain of the HLA-DR molecule, a component of the major histocompatibility complex class II, which is essential for presenting extracellular antigenic peptides to CD4⁺ T-helper cells, thereby initiating and regulating immune responses (Palojärvi et al. 2013; Matern et al. 2020).

As shown in Supplemental Table S9, the top 10 gene regulatory relationships with the highest probability achieved by PRISM-GRN are listed, along with their literature supports, which can partially validate their interactions. For example, for the first-ranking gene relationship between CTCF and BCL6, Paul et al. (Lai et al. 2010) have found that CTCF regulates BCL6 expression by binding to specific intronic regions of BCL6 in a DNA methylation-sensitive manner, in which DNA methylation can prevent CTCF binding, leading to sustained BCL6 expression. Batlle-López et al. (2015) also claimed that CTCF binds to a methylation-sensitive site within BCL6 exon1A, maintaining BCL6 expression by preventing its negative autoregulation. TRIM28 has been shown to interact with transcription CEBPB, potentially modulating its transcriptional activity through the formation of multimeric complexes (Shibata et al. 2011).

These findings show the potential downstream applications of PRISM-GRN and provide further validation for the regulatory roles of the genes in cellular processes. The integration of prior biological knowledge and multi-omics data ensures the accuracy and biological relevance of the inferred GRN, highlighting the practicability of PRISM-GRN for uncovering novel regulatory mechanisms in complex biological systems.

PRISM-GRN recovers context-specific GRNs

PRISM-GRN integrates prior GRNs to provide a structured initialization for gene–gene interactions within a specific cell type. However, rather than being strictly dictated by these priors, PRISM-GRN remains inherently data-driven. It leverages the observed omics data to refine and adjust the regulatory network, allowing for the identification of novel or context-specific interactions that may not be captured by prior knowledge alone. Therefore, we implemented PRISM-GRN on another PBMC data set, denoted as PBMC-ID, to assess its ability to find context-specific GRNs.

Specifically, the PBMC-ID data set contains only scRNA-seq data and comprises both a control group of PBMCs and an experimental group subjected to IFNB1 stimulation, further denoted as PBMC-CTRL and PBMC-IFNB1. Within this data set, the control CD4 T cells can be classified as CD4 naive T cells, whereas the IFNB1-stimulated cells exhibit an activated phenotype, despite being broadly annotated as CD4 T cells, as illustrated in Figure 8A. To maintain consistency with the PBMC data set (Fig. 7A), we utilized the same prior GRN for CD4 T cells and the corresponding scATAC-seq data. Notably, although these scATAC-seq data were paired with scRNA-seq data in the PBMC data set, it serves as nonpaired scATAC-seq data in this experiment as the PBMC-ID data set initially contains only scRNA-seq data.

Figure 8.

Illustration of PRISM-GRN for recovering context-specific GRNs. (A) t-SNE plots of the PBMC-ID data sets, with cells colored by type and stimulation status annotations. Focusing on CD4 T cells, they reveal distinct phenotypic changes induced by IFNB1 stimulation. (B) UpSet plot of the three reconstructed GRNs for CD4 T cells in the PBMC and PBMC-IFNB data sets. The reconstructed GRNs of CD4 naive T cells in both the PBMC and PBMC-CTRL data sets exhibit high consistency, whereas the GRN of IFNB1-stimulated CD4 T cells shows substantial divergence, suggesting a regulatory network shift. (C) GO enrichment analysis on the top 200 genes with the highest degree in the reconstructed GRN for IFNB1-stimulated CD4 T cells. The enriched terms are closely associated with immune response and cell development, potentially reflecting the impact of IFNB1 stimulation.

Open in new tabDownload PowerPointLink to figure

As we used the same set of prior GRNs for CD4 T cells in both the PBMC and PBMC-ID data sets, the GRNs were reconstructed based on the same test set, in which the underlying gene regulatory relationships remained unknown. The intersection structure of the reconstructed GRNs across the three conditions, PBMC, PBMC-CTRL, and PBMC-IFNB1, is depicted in Figure 8B. The intersection matrix at the bottom indicates which data sets contribute to each intersection, with black dots denoting presence in a specific condition. The vertical bars above quantify the number of shared regulatory relationships for each combination, highlighting the extent of overlap or condition-specific interactions.

With 3435 gene regulatory relationships as prior GRNs, PRISM-GRN reconstructed 9088 regulatory relationships for the CD4 T cells in the PBMC data set, 9391 in the PBMC-CTRL data set, and 9898 in the PBMC-IFNB1 data set, respectively. Notably, the reconstructed GRNs for the CD4 T cells in the PBMC and PBMC-CTRL data sets exhibit a high degree of overlap (5470 overlapped relationships). This suggests that, despite being derived from sequencing data generated in different batches, the similar phenotypic profiles of the two conditions contribute to the consistency of the omics features used. As a result, the reconstructed GRNs maintain a high level of concordance. In contrast, in the PBMC-IFNB1 data set, CD4 T cells undergo significant phenotypic changes upon IFNB1 stimulation. Consequently, although the same prior GRN was applied, the reconstructed GRN was primarily driven by the distinct omics features, leading to a markedly different regulatory network structure. In the reconstructed GRN for CD4 T cells in PBMC-IFNB1, only 77 and 97 regulatory relationships overlap with those in PBMC and PBMC-CTRL, respectively. This highlights the dynamic nature of gene regulation under immune stimulation and underscores the sensitivity of GRN inference of PRISM-GRN to underlying biological changes.

As shown in Supplemental Table S10, we present the top 10 genes with the highest degree in the inferred GRN for CD4 T cells in both the PBMC-IFNB1 and PBMC-CTRL data sets, along with their betweenness centrality. We can observe that if the context of the current data set aligns well with the conditions described by the prior GRN, the key regulators tend to retain a high degree and betweenness centrality in the reconstructed GRN. However, when the data set deviates from the prior GRN's context, for example, the IFNB1-stimulated CD4 T cells, the inferred key regulators exhibit lower degree and betweenness centrality. This further demonstrates PRISM-GRN's data-driven capability in reconstructing context-specific GRNs. Importantly, among the 10 genes for the PBMC-IFNB1 data set, JAK2 has been shown to play a crucial role in mediating IFNB1 signaling in CD4 T cells, influencing their proliferation, cytokine production, and immune function by activating STAT proteins and modulating gene expression, with evidence showing its involvement in T cell responses under inflammatory conditions (Sugimoto et al. 2003; Conti et al. 2012; Wang et al. 2024b). ETV6 is recognized as an important transcriptional regulator in hematopoietic cell development and immune function (Hock et al. 2004), whereas MAPKAPK2 acts as a key downstream effector of the p38 mitogen-activated protein kinase pathway, mediating a variety of cellular processes such as inflammatory responses, cytokine production, and cell proliferation (Jeffrey et al. 2007). Additionally, we further performed the Gene Ontology enrichment analysis based on the top 200 genes in PBMC-IFNB1. As illustrated in Figure 8C, the enriched terms are strongly associated with immune response and cell development, potentially reflecting the impact of IFNB1 stimulation.

In the PBMC-IFNB1 data set, the significant divergence in the reconstructed GRN compared with the PBMC and PBMC-CTRL data sets highlights the impact of IFNB1 stimulation on CD4 T cell phenotype and gene regulation. The results demonstrate that PRISM-GRN is highly sensitive to immune-related changes, such as those induced by IFNB1, and can identify novel gene regulatory relationships that are not captured by prior knowledge alone. In summary, PRISM-GRN's ability to integrate prior GRNs while maintaining a data-driven approach allows it to accurately reconstruct context-specific GRNs that reflect the underlying biological conditions. This capability positions PRISM-GRN as a powerful tool for investigating gene regulatory mechanisms in diverse biological contexts and offers potential applications in interventional therapy or precision medicine.

The mechanism-informed generation process

As illustrated in Figure 1B, the generation process of gene expression is formulated on the established gene regulatory mechanism, which asserts that the expression is influenced by the contributions from TF expressions, chromatin accessibility states, and the TF–gene interactions, namely, the GRNs. Therefore, E′ is assumed to follow a Dirichlet-multinomial distribution (Elkan 2006) integrating the latent variables of z_Exp from the gene expression profile, z_Acc from the chromatin accessibility score, and z_GRN from prior TF–gene interactions, which is formulated as

The distributions of the used latent variables are initialized and then constrained and guided by the mechanism-informed generative process. Specifically, the prior distribution of z_Exp and z_Acc both follow a standard normal distribution with a mean of zero and a variance of I, formulated as

Because the latent z_GRN is constrained by the prior GRN and guided by the expression generative process in a cell type–specific condition, it is supposed to indicate the hidden states for the potential GRN connections. Therefore, PRISM-GRN additionally uses a categorical distribution to recover the GRN based on z_GRN. Specifically, for two genes i and j, we first concatenate their respective z_GRN,i and z_GRN,j vectors to form a combined latent variable z_GRN,ij, which captures the regulatory relationship between them:

Finally, the regulatory relationship y_ij between gene i and gene j is generated by sampling from the categorical distribution:

The prior GRN-primed inference process

In the inference process of PRISM-GRN, the distribution parameters used in the generation process are inferred from the observable input data. A GAT (Veličković et al. 2017) is first used to obtain the initial gene embeddings, which contain the hidden information of gene regulatory interactions by incorporating a self-attention mechanism into the information propagation step. Taking as input the adjacency matrix A of the prior GRN and the gene expression profile E, GAT produces the initial gene feature embeddings z_GAT ∈ ℝ^{M × H} that are assigned the importance of their neighbors, namely, the regulatory-interacted genes, through the graph connectivity according to the prior GRN, where H is the size of the latent embedding.

More specifically, a shared weight matrix W ∈ ℝ^{H × N} is first applied to perform a linear transformation on the raw feature, the self-attention mechanism is subsequently implemented using a feedforward MLP layer parameterized by a weight vector α ∈ ℝ^2H and the LeakyReLU activation function is followed to introduce nonlinearity. Taking gene i and gene j as examples, the attention coefficient of gene j on gene i can be formulated as follows:

The attention coefficients of gene j of gene i are subsequently normalized across all the first-order neighbors using Softmax, which is defined as follows:

Finally, the initial embedding of gene i can be obtained by aggregating its neighbor genes through the single-head attention mechanism:

z_GAT produces the primary hidden state for each gene by incorporating the known regulatory information from prior GRN. To enable the causality prediction among the gene pairs, z_GRN is further inferred from z_GAT to learn the unobserved biological factor from GRN. The variational posterior distribution q(z_GRN|z_GAT) is used to approximate the real posterior distribution p(z_GRN|z_GAT) in the variational inference, which is parameterized as

As mentioned in the generation process (Fig. 1B), PRISM-GRN subsequently integrates the latent variable z_GAT with z_Exp and z_Acc into a mechanism-informed generation process. Similarly, given the observed data, we use an inference process to infer the posterior distribution parameters to be used in the generation process, including z_Exp from E, z_Acc from R.

Specifically, z_Exp captures the unobserved biological factors influencing gene expression of TFs; z_Acc captures the unobserved biological factors influencing the accessible state of target genes; and their posterior distributions are parameterized as

Notably, the reparameterization trick is implemented automatically by using Pyro (Bingham et al. 2019) to enable efficient and differentiable sampling of the latent variables z_Exp, z_Acc, and z_GRN, during the inference process. Based on the learned latent variables that represent the biological factors in gene expression, PRISM-GRN can accordingly perform the mechanism-informed generation process.

Loss function and parameter optimization

PRISM-GRN is a probabilistic model, in which the parameters are optimized through variational inference using the Adam optimizer. The loss function comprises three key components: (1) the KL divergence loss, (2) the reconstruction loss for gene expression data, and (3) the reconstruction loss for the GRN. Each of these components ensures that the model both approximates the true posterior distributions and effectively reconstructs the observed data.

The KL divergence is used to measure the divergence between the approximate posterior distributions of the latent variables and their corresponding prior distributions. Specifically, for the latent variables z_Exp, z_Acc, z_GRN, the KL divergence is computed between the approximate posterior distributions q(z|*) and the standard Gaussian prior distributions $p\left(z\right)\sim \mathcal{N}\left(0,I\right)$:

The reconstruction loss for gene expression data E is computed as the mean squared error (MSE) between the observed expression data and the reconstructed expression data E′, which is generated by the model from the latent variables z_Exp, z_Acc, and z_GRN:

The total loss function is the sum of KL divergence, gene expression reconstruction loss, and GRN reconstruction loss, as defined

Data source

All data sets used in this paper were downloaded from public databases. The paired scRNA-seq and scATAC-seq data for the A549 data set were obtained from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE117089 (Cao et al. 2018). The paired scRNA-seq and scATAC-seq data for data sets GM12878, K562, and 3T3 are from the SHRAE-seq data under GEO accession number GSE140203 (Ma et al. 2020), in which the replicant scATAC-seq data for GM12878 were also embraced. The unpaired scATAC-seq data for A549 were obtained from the study under GEO accession number GSE211617. For the additional two data sets, hESC and mESC, the scRNA-seq data were obtained from the study under GEO accession number GSE75748 (Chu et al. 2016) and GSE98664 (Hayashi et al. 2018), whereas the scATAC-seq data were independently obtained from the study under GEO accession number GSE192858 (Samara et al. 2022) and GSE107651 (Zamanighomi et al. 2018), respectively. All cell type–specific ChIP-seq data were obtained from the ChIP-Atlas database (Oki et al. 2018). The PBMC data sets can be accessed from 10x Genomics platform with the link https://www.10xgenomics.com/datasets/pbmc-from-a-healthy-donor-no-cell-sorting-10-k-1-standard-2-0-0. The PBMC-IFNB1 data sets were obtained from the study under GEO accession number GSE96583 (Kang et al. 2018).

Code availability

The source code of PRISM-GRN and the preprocessing scripts of this paper for reproducing the results are available at GitHub (http://github.com/Ying-Lab/PRISM-GRN), under the MIT License, and as Supplemental Code.