Decoding apoptosis, ferroptosis, and inflammatory cell death in adenomyosis at single-cell resolution

To systematically evaluate the ability of GraphPAS to resolve cellular heterogeneity at the pathway level, we performed benchmarking on three publicly available scRNA-seq datasets, including the human pancreas dataset from Muraro et al . [ 21 ], the human testis dataset from Zhao et al . [ 22 ], and the human islet dataset from Segerstolpe et al . [ 23 ]. For each dataset, pathway activity matrices were generated using GraphPAS, AUCell, and scapGNN, respectively, and were used as input features for downstream cell clustering. We then applied 10 commonly used single-cell clustering algorithms to the pathway-level feature matrices produced by each method. Clustering performance was evaluated using ARI, NMI, and SW, enabling systematic comparison across methods. Compared with AUCell and scapGNN, GraphPAS consistently achieved the best performance across all three datasets and clustering algorithms ( Fig. 2 ). These results indicate that pathway activity representations generated by GraphPAS more effectively distinguish cell types and capture functional differences associated with cellular states at the pathway level. Relative to AUCell and scapGNN, the pathway-level feature space learned by GraphPAS exhibits clearer cell-type separation, thereby improving clustering accuracy based on pathway features. Moreover, GraphPAS demonstrates more stable performance across datasets, suggesting that its pathway representations generalize well across tissues with distinct cellular compositions. Collectively, these findings show that GraphPAS produces more discriminative and robust pathway activity representations, enabling more accurate characterization of cellular heterogeneity in single-cell data at the pathway level. Clustering performance comparison of GraphPAS, scapGNN, and AUCell; pathway activity matrices inferred by each method were used for downstream clustering with 10 clustering algorithms; clustering performance was evaluated using ARI (left), NMI (middle), and SW (right) across three benchmark scRNA-seq datasets, including the Segerstolpe, Muraro, and Zhao datasets. To evaluate pathway identification performance in the absence of gold-standard pathway annotations for single-cell data, we curated cell type-specific marker gene sets from the CellMarker, BioGPS, and Harmonizome databases as reference functional programs [ 24–27 ]. These curated marker gene sets, representing well-defined cell type-associated functional signatures, were used to assess whether different pathway activity scoring methods preferentially prioritize biologically relevant programs within the corresponding cell types. Detection accuracy of cell-type-specific marker gene sets using GraphPAS, scapGNN, and AUCell; for each individual cell, pathway activity scores of all reference marker gene sets were ranked, and a recovery event was defined when the marker gene set corresponding to the annotated cell type appeared among the top- k highest-scoring marker programs; Top1 represents the most stringent criterion, requiring the biologically matched marker program to receive the highest pathway activity score within a cell, whereas Top5 indicates recovery within the five highest-scoring marker programs; bars represent the proportion of cells satisfying the Top1–Top5 recovery criteria for each method; GraphPAS consistently achieved higher recovery rates across all top- k thresholds, indicating improved identification of cell type-associated functional programs. Robustness evaluation of GraphPAS under simulated technical noise; dropout noise was introduced by randomly masking nonzero expression values, and Gaussian noise was added by perturbing expression values with Gaussian-distributed noise; clustering performance derived from pathway activity matrices was evaluated across noise levels, and robustness was summarized using the area under the noise-performance curve (AUC); top panels show AUC values under dropout noise, and bottom panels show AUC values under Gaussian noise; left panels correspond to ARI, and right panels correspond to NMI. Using the embryonic stem cell scRNA-seq dataset from Yan et al . [ 28 ] as a benchmark, we computed pathway activity scores with GraphPAS, AUCell, and scapGNN, respectively. For each individual cell, pathway activity scores were ranked, and we quantified the proportion of cells in which the biologically matched cell type-associated marker program appeared among the top- k highest-scoring functional programs. Thus, Top1 represented the most stringent criterion, requiring the biologically matched functional program to receive the highest pathway activity score within a cell, whereas Top5 evaluated whether the corresponding cell type-associated program can be prioritized among the five highest-scoring pathways. We then calculated the proportion of cells satisfying this criterion across Top1–Top5 thresholds. GraphPAS consistently identified a higher proportion of cells with correctly prioritized cell type-associated functional programs across all Top- k thresholds ( Fig. 3 ), indicating improved recovery of biologically relevant pathway signals. Importantly, this metric does not measure absolute pathway accuracy, but instead evaluates whether a pathway scoring method preferentially assigns high activity scores to functional programs consistent with known cellular identity. The superior performance of GraphPAS likely stems from its ability to integrate transcriptional similarity structure across both genes and cells through hierarchical graph representation learning, resulting in more stable pathway activity estimation under sparse single-cell conditions. Collectively, these results demonstrate that GraphPAS more effectively captures cell identity-associated functional programs at single-cell resolution. Macrophages display prominent PCD-related activity in the human uterine single-cell transcriptome; (a) t-SNE visualization of single-cell transcriptomes from human uterine tissue colored by annotated cell types (left) and projection of the overall PCD score computed from apoptosis, ferroptosis, necroptosis, NETosis, and pyroptosis signatures (right); (b) heatmap showing average PCD signature scores across major uterine cell populations, including endothelial cells, epithelial cells, macrophages, mural cells, stromal cells, and T cells, highlighting the enrichment of multiple cell death-related programs in macrophages; (c) comparison of PCD signature scores in macrophages between the WithPain and WithoutPain groups; violin plots show enrichment scores for five PCD programs, including apoptosis, ferroptosis, necroptosis, NETosis, and pyroptosis; statistical significance was assessed using a two-sided Wilcoxon rank-sum test. To assess the robustness of GraphPAS under noisy conditions, we introduced artificial perturbations into the single-cell expression matrix and used the area under the noise-level clustering performance curve (AUC) as a comprehensive evaluation metric. These perturbation analyses were designed to mimic two major sources of technical variability commonly observed in scRNA-seq data: sparsity caused by dropout events and continuous measurement noise arising from sequencing and quantification variability. Specifically, dropout noise was simulated by randomly setting a proportion of nonzero expression values to zero, with noise levels of 5%, 10%, 15%, and 20%. At each noise level, pathway activity scores were recomputed, followed by cell clustering based on pathway-level features. Clustering performance across noise conditions was then summarized using AUC ( Fig. 4 ), where higher AUC values indicate greater preservation of cell-type structure under increasing noise conditions. As dropout noise increased, GraphPAS consistently achieved higher AUC values for both ARI and NMI compared with AUCell and scapGNN, indicating improved robustness to sparsity-induced noise. We further introduced Gaussian noise into the expression matrix to simulate continuous measurement errors and repeated the same analysis ( Fig. 4 ). Consistent with the dropout experiments, GraphPAS again showed higher AUC values for both ARI and NMI than AUCell and scapGNN. These results demonstrate that GraphPAS maintains strong robustness under both sparse dropout noise and continuous measurement noise, enabling stable recovery of biologically meaningful pathway activity patterns in noisy single-cell data. We then applied GraphPAS to examine programmed cell death (PCD)-related functional states in adenomyosis. The analysis used the single-cell transcriptomic dataset reported by Chen et al . [ 19 ], comprising two adenomyosis samples with pain and two control samples. After quality control and preprocessing, 33 617 cells were retained for downstream analysis (Methods). After dimensionality reduction and cell type annotation, six major cell populations were identified: endothelial cells, epithelial cells, macrophages, mural cells, stromal cells, and T cells ( Fig. 5a ). GraphPAS was used to quantify pathway activity for five PCD-related programs: apoptosis, ferroptosis, necroptosis, NETosis, and pyroptosis ( Supplementary Table S1 ). Projection of overall PCD activity onto the t-distributed stochastic neighbor embedding (t-SNE) space revealed that cell death-related programs were broadly distributed across multiple cell types, while exhibiting pronounced heterogeneity both across and within cell populations ( Fig. 5a ), suggesting substantial remodeling of cell death states within the uterine microenvironment. Comparison across cell types further showed that macrophages had the strongest enrichment of death-related pathways ( Fig. 5b ), suggesting that they are a major cellular compartment associated with altered PCD activity in adenomyosis. We therefore focused on macrophages and compared PCD activity between the pain-associated (WithPain) and control (WithoutPain) groups. Macrophages from the WithPain group displayed significantly elevated apoptosis, ferroptosis, and necroptosis activity, whereas NETosis showed no significant difference between groups ( Fig. 5c ). These findings indicate that macrophages represent the most prominent cellular compartment with aberrant activation of cell death programs in pain-associated adenomyosis, particularly characterized by increased apoptosis, ferroptosis, and inflammatory PCD-related programs. To investigate the molecular context of this activated cell death state, we performed KEGG and GO enrichment analyses of differentially expressed genes in macrophages. KEGG analysis revealed significant enrichment in multiple inflammation- and immune-related pathways, including the interleukin-17 (IL-17) signaling pathway, chemokine signaling pathway, phagosome, and other immune response-associated pathways ( Fig. 6a ). Consistently, GO analysis indicated that these genes were primarily involved in biological processes such as leukocyte proliferation, lymphocyte proliferation and differentiation, immune cell adhesion, and regulation of T cell activation. In addition, enrichment was observed in cellular component and molecular function categories related to membrane-associated structures, ribosome-related terms, transcriptional regulation, and immune-binding activities ( Fig. 6b ). Collectively, GraphPAS reveals pronounced cell type-specific heterogeneity in PCD programs in adenomyosis and identifies macrophages as the central cellular population underlying pain-associated remodeling of cell death states. In the pain-associated condition, macrophages simultaneously exhibit increased cell death-related activity and immune-inflammatory transcriptional signatures, suggesting that macrophage-associated reprogramming of PCD may play a key role in pathological remodeling of the uterine microenvironment and the development of pain. Functional characterization of macrophage-associated cell death-related transcriptional programs; (a) KEGG pathway enrichment analysis of differentially expressed genes in macrophages between the WithPain and WithoutPain groups; the circular network plot shows enriched pathways and their associated genes; node size indicates the number of genes within each pathway, and node color represents the log 2 fold change of gene expression; (b) GO enrichment analysis of macrophage-related differential genes; bar plots display significantly enriched biological process (BP), cellular component (CC), and molecular function (MF) categories; enriched terms highlight immune cell proliferation, lymphocyte activation and differentiation, leukocyte adhesion, and ribosome- and membrane-associated components, indicating an immune-inflammatory transcriptional state in macrophages.

GraphPAS is an unsupervised framework based on hierarchical graph representation learning, designed to infer stable gene–cell association structures from single-cell expression matrices and subsequently quantify pathway activity at single-cell resolution. Unlike approaches that directly aggregate gene sets from raw expression values, GraphPAS formulates the single-cell expression matrix as a graph learning problem comprising gene nodes, cell nodes, and multi-layer relational structures. Specifically, observed gene–cell relationships define cross-type associations, while local similarities among genes and among cells provide topological constraints for representation learning. Within this framework, GraphPAS jointly learns low-dimensional embeddings for genes and cells, reconstructs a more robust gene–cell score matrix, and leverages this refined representation to estimate cell state-associated pathway activity. GraphPAS consists of three major steps ( Fig. 1 ): (i) first, two independent autoencoders are employed to learn low-dimensional representations of genes and cells from the single-cell expression matrix. (ii) Based on the learned gene and cell embeddings, global gene–gene and cell–cell graphs are constructed, and homogeneous graph attention networks are applied to refine the local topological structure of each node type. Building upon this, a HGT performs cross-type message passing to jointly update gene and cell embeddings, yielding globally consistent gene–cell association representations. (iii) Finally, a decoder reconstructs the gene–cell score matrix, which is then combined with predefined seed gene sets to quantify pathway activity, enabling identification of cell state heterogeneity associated with specific biological processes. Overview of GraphPAS; GraphPAS infers pathway activity through hierarchical graph representation learning; dual autoencoders learn representations of genes and cells from single-cell expression profiles; gene–gene and cell–cell similarity graphs are constructed and refined using graph attention networks to capture local structure; an HGT then performs cross-type message passing to jointly update gene and cell representations; the resulting representations are used to reconstruct a gene–cell association matrix, from which pathway activity is computed by aggregating predefined pathway gene sets. GraphPAS employs two independent autoencoders to learn low-dimensional embeddings of cells and genes from the single-cell expression matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $X=\left\{{x}_{ij}\right\}\in{\mathrm{R}}^{G\times C}$\end{document} . From the cell-centric perspective, each cell is represented by its expression vector across all genes, whereas from the gene-centric perspective, each gene is characterized by its activity profile across all cells. The initial cell embedding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Z}_c^{(0)}$\end{document} and gene embedding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Z}_g^{(0)}$\end{document} are learned by minimizing the mean squared error between the input and reconstructed representations, defined as follows: Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${x}_j^{(c)}$\end{document} denotes the input vector of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} -th cell, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${x}_i^{(g)}$\end{document} denotes the input vector of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} -th gene. GraphPAS further constructs a gene–gene graph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_{gg}=\left({V}_g,{E}_{gg}\right)$\end{document} and a cell–cell graph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_{cc}=\left({V}_c,{E}_{cc}\right)$\end{document} to explicitly encode local similarity structures among nodes of the same type. Edges in both graphs are established based on neighborhood similarity in the latent embedding space. To incorporate these local topological structures into node representations, GraphPAS applies graph attention networks separately on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_{gg}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_{cc}$\end{document} . For an arbitrary node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $u$\end{document} , the updated representation is computed by aggregating its own features with those of its neighboring nodes using attention-based weighting: Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${h}_v$\end{document} denotes the input feature of node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $v$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $W$\end{document} is a learnable linear transformation matrix, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\alpha}_{uv}$\end{document} represents the attention coefficient measuring the contribution of node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $v$\end{document} to the update of node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $u$\end{document} . The function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\sigma$\end{document} denotes a nonlinear activation. The attention mechanism adaptively weights neighboring nodes according to their relative importance, allowing the model to differentially aggregate informative neighbors when updating node representations. GraphPAS organizes genes and cells into a bipartite heterogeneous graph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $G=\left(V,E,A,R\right)$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $V$\end{document} denotes the set of nodes, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $E$\end{document} the set of edges, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document} the set of node types, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $R$\end{document} the set of relation types. For each element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${x}_{ij}$\end{document} in the expression matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $X$\end{document} , an edge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${e}_{ij}$\end{document} is introduced between the gene node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_i^g$\end{document} and the cell node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_j^c$\end{document} if gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} exhibits observed activity in cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} . Accordingly, the edge set is defined as: To jointly learn representations of genes and cells on the gene–cell heterogeneous graph, GraphPAS employs a HGT as the core encoder [ 16 ]. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${H}^{(l)}\left[v\right]$\end{document} denote the representation of node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $v$\end{document} at layer \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $l$\end{document} . In each layer, a node is treated as a target node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_t$\end{document} , while its neighbors are regarded as source nodes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_s$\end{document} . Information is aggregated from source nodes to update the target representation through node type-specific and relation type-specific parameterized transformations. For the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $h$\end{document} -th attention head, the target and source nodes are projected into query, key, and value spaces, respectively: Gene nodes and cell nodes correspond to two distinct representational entities: the former capture gene activity patterns across the cellular population, whereas the latter characterize cellular molecular states in the gene expression space. GraphPAS learns separate projection matrices for different node types, enabling type-specific adaptation prior to message passing and facilitating more effective modeling of cross-type associations in the heterogeneous graph. For an edge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${e}_{s,t}$\end{document} , the relation-specific attention score is defined as follows: Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${W}_{\phi \left({e}_{s,t}\right)}^{ATT}$\end{document} denotes the relation-specific attention transformation matrix, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${d}_h$\end{document} represents the dimensionality of each attention head. The attention coefficients are normalized across all neighbors of the same target node using a softmax function, yielding the attention weights: During the message passing stage, the source node features are first transformed through a relation-specific message projection: Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${W}_{\phi\ \left({e}_{s,t}\right)}^{MSG}$\end{document} denotes the relation-specific message transformation matrix. The final representation of the target node is obtained by aggregating messages from all neighbors across all attention heads: Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\Big\Vert$\end{document} denotes multi-head concatenation. The aggregated representation is then combined with the previous layer’s embedding through a residual connection to produce the updated node representation: Through layer-wise heterogeneous propagation, GraphPAS ultimately obtains refined gene embeddings \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Z}_g$\end{document} and cell embeddings \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Z}_c$\end{document} , which jointly encode cross-type associations, local homogeneous structures, and higher-order topological dependencies within a shared latent space. GraphPAS adopts a mini-batch training strategy based on local subgraphs rather than optimizing over the full graph. Specifically, the model starts from a set of seed genes and constructs local heterogeneous subgraphs through alternating gene-to-cell and cell-to-gene neighborhood expansion. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${V}_g^{(b)}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${V}_c^{(b)}$\end{document} denote the sets of gene nodes and cell nodes in the current subgraph, respectively. The corresponding local gene–cell adjacency matrix is defined as: For each sampled subgraph, the global K -nearest-neighbor gene–gene and cell–cell graphs are cropped to retain only nodes present in the current batch. This strategy allows each training batch to preserve both cross-type gene–cell relationships and local same-type neighborhoods. The HGT encoder maps nodes in the sampled subgraph to low-dimensional representations, and an inner-product decoder reconstructs the local gene–cell association matrix: Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\hat{A}}^{(b)}$\end{document} is not interpreted as a direct reconstruction of the original expression matrix, but rather as a locally refined gene–cell association score matrix learned under graph structural constraints. A higher reconstruction score indicates a stronger relative association of a gene within the current cellular context, suggesting a greater likelihood that the gene contributes to defining the corresponding cell state. To enforce consistency between the observed adjacency matrix and the reconstructed score matrix, GraphPAS employs Kullback–Leibler (KL) divergence to align their normalized distributions: GraphPAS constructs a global gene–cell score matrix using the final gene and cell embeddings: GraphPAS then standardizes gene scores within each cell. Specifically, for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c$\end{document} -th cell, the normalized gene–cell score is defined as: where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mu}_c$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\sigma}_c$\end{document} denote the mean and standard deviation of scores across all genes in cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c$\end{document} , respectively, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\varepsilon$\end{document} is a small constant for numerical stability. Given a pathway seed gene set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P=\left\{{g}_1,{g}_2,\dots, {g}_m\right\}$\end{document} , the pathway activity in cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c$\end{document} is defined as where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathrm{Agg}\left(\bullet \right)$\end{document} denotes an aggregation operation applied to the normalized scores of pathway genes. In this formulation, the resulting pathway activity score no longer relies on direct aggregation of raw expression values, but instead is derived from the learned gene–cell association structure, thereby emphasizing the relative functional importance of pathway genes within specific cellular states. We benchmarked pathway activity matrices generated by AUCell and scapGNN within the sciPath framework [ 5 ]. To systematically evaluate GraphPAS under comparable analytical conditions, all benchmark datasets were processed using a unified preprocessing and evaluation framework. Raw count matrices were obtained from the original publications. Pathway gene sets were obtained from the KEGG, with genes in each pathway treated as pathway terms. For each dataset and method, pathway activity matrices were generated with pathways as features and cells as observations, and the resulting matrices were supplied to sciPath for downstream clustering analysis. This design ensured that the compared methods were evaluated using the same downstream clustering framework. Within sciPath, 10 state-of-the-art clustering algorithms were applied to infer cell clusters based on pathway activity representations, using the implementation and parameter settings provided by sciPath. Specifically, the Leiden algorithm implemented in Seurat was evaluated at three resolutions (0.5, 1.0, and 1.5). AUCell scores were computed using the AUCell_calcAUC function from the R package AUCell (version 1.8.0) with default parameters. scapGNN-derived pathway activity matrices were generated using the default model configuration. All pathways in each activity matrix were used as input features for clustering. Clustering performance was evaluated using three widely used metrics: adjusted Rand index (ARI), normalized mutual information (NMI), and silhouette width (SW). ARI and NMI quantify the agreement between inferred clusters and reference cell labels, with higher values indicating better clustering accuracy and consistency [ 17 , 18 ]. SW measures the separation and compactness of clusters in the feature space, with higher values reflecting improved intra-cluster similarity and inter-cluster separation [ 7 ]. To further evaluate GraphPAS on real biological data, we analyzed the 10× Genomics scRNA-seq dataset reported by Chen et al . [ 19 ], which comprises four uterine tissue samples, including two adenomyosis cases with dysmenorrhea and two control samples. Raw data were reprocessed using the Seurat package (v5) in R. For each sample, gene expression matrices were independently imported and converted into Seurat objects. Quality control was performed to remove low-quality and potential doublet cells. Cells with total unique molecular identifier (UMI) counts within the top 2% of each sample were removed to mitigate potential doublets or transcriptionally aberrant cells. Cells with fewer than 200 detected UMIs or with mitochondrial transcript proportions exceeding 20% were also excluded. After filtering, all samples were merged for downstream analysis. Expression values were normalized prior to feature selection and dimensionality reduction. Highly variable genes were identified to capture major sources of transcriptional heterogeneity, and scaled expression values were used for principal component analysis (PCA). To correct batch effects across samples, Harmony (v0.1.0) [ 20 ] was applied in PCA space. The Harmony-corrected embeddings were subsequently used to construct the cell neighborhood graph and perform unsupervised clustering. Multiple clustering resolutions ranging from 0.1 to 1.2 were evaluated, with cluster stability assessed using clustree (v0.5.0), and a resolution of 0.8 was selected for downstream analysis. For visualization, uniform manifold approximation and projection was computed based on the first 40 Harmony dimensions. Cluster-specific marker genes were identified using FindAllMarkers, retaining only positively enriched genes detected in at least 25% of cells with a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\log}_2$\end{document} fold change ≥0.25. Cell identities were assigned based on canonical marker genes and validated by marker expression patterns.

In this study, we developed GraphPAS, a hierarchical graph representation learning framework for pathway activity scoring in scRNA-seq data. By jointly modeling gene–gene, cell–cell, and gene–cell relationships, GraphPAS learns a robust gene–cell association landscape and uses it to infer pathway activity at single-cell resolution. Across multiple benchmark datasets, GraphPAS improved clustering-based evaluation, marker gene set recovery, and robustness under noise perturbation relative to existing methods. These findings suggest that structured relational information can improve pathway-level representation of cellular states beyond direct aggregation of raw gene expression. The advantage of GraphPAS is linked to the structure of single-cell pathway signals. Biologically relevant programs may involve coordinated but modest changes across multiple genes, whereas scRNA-seq measurements are sparse and noisy. Direct gene set aggregation can therefore be sensitive to unstable gene-level values. GraphPAS instead estimates pathway activity from a graph-refined gene–cell score matrix, allowing pathway inference to use observed expression together with local neighborhood structure and cross-type dependencies. This design explains why the framework improved cell type discrimination, biologically expected pathway recovery, and robustness to dropout and Gaussian noise. Compared with previous graph-based methods, GraphPAS directly couples hierarchical representation learning with pathway quantification, rather than treating pathway analysis as a downstream step after graph embedding. Application of GraphPAS to adenomyosis further demonstrated its utility in resolving disease-associated functional states. By quantifying PCD-related pathway activity at single-cell resolution, we identified marked cell type-specific heterogeneity in the uterine microenvironment and found that macrophages showed the strongest enrichment of cell death-related programs. In pain-associated adenomyosis, macrophages exhibited increased apoptosis, ferroptosis, and necroptosis activity, together with enrichment of inflammatory and immune-related pathways, including IL-17 signaling, chemokine signaling, and phagosome-related processes. These findings suggest that macrophages may represent a central cellular compartment linking aberrant cell death programs with inflammatory activation in adenomyosis and highlight the value of pathway-oriented analysis for uncovering functionally relevant disease mechanisms. Several limitations should be noted. First, GraphPAS infers pathway-associated transcriptional states rather than direct biochemical pathway activity and may therefore not fully capture regulation at the protein, post-translational, or metabolic level. Second, as no gold-standard pathway activity labels are available for single-cell data, our evaluation relies mainly on indirect criteria, including clustering performance, marker set prioritization, and biological consistency in a disease setting. Future extensions could incorporate pathway topology, regulatory directionality, and multimodal constraints from chromatin accessibility, spatial transcriptomics, or proteomic data to further improve functional interpretability and pathway activity estimation. Overall, GraphPAS provides an effective framework for structure-aware pathway analysis in single-cell data and offers a practical strategy for linking graph representation learning with functional characterization of complex cellular states. Graph-based Pathway Activity Scoring (GraphPAS) provides a graph-based framework for inferring pathway activity from sparse and noisy single-cell RNA-sequencing data. GraphPAS integrates gene–gene, cell–cell, and gene–cell relationships to recover robust pathway-level cellular representations. GraphPAS outperforms AUCell and scapGNN in clustering accuracy, marker program recovery, and noise robustness. GraphPAS identifies macrophages as the major programmed cell death-enriched population in adenomyosis, with pain-associated lesions exhibiting enhanced apoptosis, ferroptosis, necroptosis, and immune-inflammatory activation.

Biological pathways comprise functionally coordinated regulatory units formed by the concerted action of multiple genes, reflecting the organized execution of intracellular molecular processes [ 1 , 2 ]. Resources such as the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) organize experimental evidence and prior knowledge into pathways or functional modules, providing a foundation for interpreting cellular states at the functional level [ 3 , 4 ]. In single-cell RNA sequencing (scRNA-seq) studies, pathway-level activity provides a higher-order representation of cellular programs compared with individual gene expression, which is often highly variable and susceptible to technical noise. Consequently, pathway activity has become an important strategy for dissecting cellular heterogeneity, identifying cell states, and uncovering disease-associated molecular mechanisms [ 5 , 6 ]. However, accurate pathway activity inference at single-cell resolution remains difficult because dropout events and measurement noise can obscure the coordinated signals that define functional programs [ 7 ]. Most pathway enrichment methods were developed for bulk transcriptomic data, where predefined gene sets are summarized using expression distributions, rank statistics, or related aggregate measures [ 8–11 ]. When adapted to single-cell data, these methods often treat pathways as independent gene collections and score them directly from sparse expression matrices. This strategy does not explicitly model gene–gene dependency or the organization of cellular states, and it can miss coordinated pathway activation driven by multiple weakly expressed genes [ 12 ]. As a result, aggregation-based scoring may overlook the cooperative network structure that underlies pathway activity in sparse single-cell data. Graph representation learning has emerged as a promising methodological paradigm for single-cell data analysis [ 13–15 ]. By introducing structured relationships between cells and genes and leveraging local topological information for representation learning, these approaches can derive more robust feature representations from highly sparse and noisy single-cell datasets. Graph-based models, exemplified by scapGNN, further demonstrate that incorporating multi-level relational structures facilitates the recovery of stable functional associations from single-cell expression matrices and supports downstream pathway analysis [ 7 ]. However, existing graph-based approaches are primarily designed to reconstruct association structures or learn low-dimensional embeddings, with pathway analysis typically performed as a downstream task rather than being directly optimized within the model. A unified framework is still needed to refine structural relationships, integrate cross-level information, and directly produce pathway-level representations for single-cell analysis. To address this need, we developed Graph-based Pathway Activity Scoring (GraphPAS), a graph-based framework for robust pathway activity inference at single-cell resolution. GraphPAS models genes and cells as connected entities, learns local gene–gene and cell–cell structures, and integrates cross-type associations through a heterogeneous graph transformer (HGT). The learned gene-cell association matrix is then transformed directly into pathway activity scores, rather than being used only as an intermediate embedding for downstream analysis. This design couples structural learning with pathway quantification in a single model, allowing GraphPAS to capture coordinated gene variation and similarity among cellular states. Systematic evaluation across multiple single-cell datasets demonstrates that GraphPAS consistently outperforms existing approaches in resolving cellular heterogeneity, identifying cell type-associated pathways, and maintaining robustness under noisy conditions.