PathNetDRP: a novel biomarker discovery framework using pathway and protein–protein interaction networks for immune checkpoint inhibitor response prediction

Abstract

Background

Predicting immune checkpoint inhibitor (ICI) response remains a significant challenge in cancer immunotherapy. Many existing approaches rely on differential gene expression analysis or predefined immune signatures, which may fail to capture the complex regulatory mechanisms underlying immune response. Network-based models attempt to integrate biological interactions, but they often lack a quantitative framework to assess how individual genes contribute within pathways, limiting the specificity and interpretability of biomarkers. Given these limitations, we developed PathNetDRP, a framework that integrates biological pathways, protein-protein interaction networks, and machine learning to identify functionally relevant biomarkers for ICI response prediction.

Results

We introduce PathNetDRP, a novel biomarker discovery approach that applies the PageRank algorithm to prioritize ICI-associated genes, maps them to relevant biological pathways, and calculates PathNetGene scores to quantify their contribution to immune response. Unlike conventional methods that focus solely on gene expression differences, PathNetDRP systematically incorporates biological context to improve biomarker selection. Validation across multiple independent cancer cohorts showed that PathNetDRP achieved strong predictive performance, with cross-validation the area under the receiver operating characteristic curves increasing from 0.780 to 0.940. Interestingly, PathNetDRP did not merely improve predictive accuracy; it also provided insights into key immune-related pathways, reinforcing its potential for identifying clinically relevant biomarkers.

Conclusion

The biomarkers identified by PathNetDRP demonstrated robust predictive performance across cross-validation and independent validation datasets, suggesting their potential utility in clinical applications. Furthermore, enrichment analysis highlighted key immune-related pathways, providing a deeper understanding of their role in ICI response regulation. While these findings underscore the promise of PathNetDRP, future work will explore the integration of additional predictive features, such as tumor mutational burden and microsatellite instability, to further refine its applicability.

Background

Immune checkpoint inhibitors (ICIs) have revolutionized the landscape of cancer therapy, offering new hope to patients with advanced malignancies [1, 2]. Unlike traditional therapies that directly target cancer cells, ICIs, such as nivolumab and pembrolizumab, promote the immune-mediated elimination of tumor cells. These therapies have demonstrated remarkable success in clinical trials, significantly improving survival rates across various cancers. For instance, a phase 3 study showed that nivolumab increased the overall survival rate in patients with metastatic melanoma without a BRAF mutation to 72.9%, compared to 42.1% with dacarbazine [3]. Despite these promising results, only a minority of patients respond favorably to ICIs [1, 4], underscoring the urgent need for predictive biomarkers to identify which patients will benefit from these treatments.

Numerous studies have identified several biomarkers associated with ICI response, shedding light on the complex interactions within the tumor microenvironment (TME). Key predictors include the expression of programmed death-ligand 1 (PD-L1) and mismatch repair deficiency [5,6,7]. PD-L1, extensively studied and typically measured by immunohistochemistry, varies in its expression cutoffs across different cancer types such as advanced non-small-cell lung cancer (NSCLC), urothelial cancer, and triple-negative breast cancer [5, 8, 9]. Mismatch repair status is another well-known predictor, as deficiency in mismatch repair significantly promotes the expression of several immune checkpoint ligands, including PD-1, PD-L1, CTLA-4, LAG-3, and IDO [6, 10].

Moreover, gene expression profiles, such as the T cell-inflamed gene expression profile (GEP), and tumor mutation burden along with specific genetic alterations, have been linked to better outcomes with ICIs [11, 12]. Components of the TME, like cancer-associated fibroblasts (CAFs) and tumor-associated macrophages (TAMs), also play significant roles in modulating immune responses [13, 14]. CAFs influence cancer progression and treatment efficacy through their secretion of cytokines and growth factors, impacting the recruitment and regulation of both innate and adaptive immune cells [13].

With the advent of bioinformatics technologies and the accumulation of vast omics datasets, omics-based gene prioritization methods have shown promise in predicting ICI responses [15, 16]. The TIDE model, for example, predicts ICI response by modeling T cell dysfunction and exclusion, offering more accuracy than assessments based on PD-L1 expression or mutation load alone [15]. It has also identified new candidate regulators of ICI resistance. Similarly, the IMPRES predictor, which utilizes pairwise transcriptomics relations between immune checkpoint genes, has identified 15 predictive features for ICI response in melanoma, achieving high accuracy across multiple published datasets [16].

The integration of machine learning (ML) techniques has significantly enhanced the predictive accuracy and understanding of immunotherapy responses [17]. ML models can process and analyze large, complex datasets, including genomic, transcriptomic, and proteomic data, to uncover patterns and relationships that might not be evident through traditional methods. For instance, DeepGeneX utilizes a deep neural network model and feature elimination to reduce single-cell RNA-seq data from 26,000 genes to six key genes, identifying potential therapeutic targets in high LGALS1 and WARS-expressing macrophage populations [17]. These methodologies highlight the potential of ML and deep learning in advancing precision oncology, offering robust tools for predicting immunotherapy responses and optimizing treatment strategies.

While numerous biomarker identification techniques have been developed to predict ICI responses, there remains a significant need for further analysis of biomarkers and mechanisms related to drug responses, as well as enhancements in prediction accuracy. For instance, TIDE [15] can identify biomarkers based on genes associated with tumor immune dysfunction and exclusion using Pearson’s correlation analysis. However, the predictive performance of TIDE is limited by the immune system’s complexity, which results from the interaction of various factors. DeepGeneX [17] applies deep learning to select features associated with ICI response, but its effectiveness is hindered by both the limited size of available ICI-related datasets, which are often insufficient for training deep learning models effectively, and the “black box” nature of deep learning, which makes interpretation challenging. IMPRES [16] analyzes combinations of known genes related to ICI response, which restricts its utility for discovering new biomarkers. To address these challenges in understanding drug response mechanisms and improving predictive accuracy, an approach that leverages biological network information could be applied.

The network propagation algorithm has been combined with protein-protein interaction (PPI) networks to identify biomarkers for predicting drug responses [18, 19]. The NetBio framework, which integrates a network propagation algorithm with pathway enrichment analysis, identifies biomarkers associated with ICI response using clinical outcomes and transcriptomic data from over 700 ICI-treated patients [18]. This approach has demonstrated superior predictive accuracy in melanoma, gastric cancer, and bladder cancer compared to conventional biomarkers. ICINet [19] applies PageRank to PPI networks to identify selected nodes and then constructs a gene regulatory integration graph using 14 prior knowledge databases, including gene ontology. This regulatory graph is then used with a graph neural network (GNN) to predict immune therapy response in cancer. These studies underscore the potential of combining network propagation algorithms with PPI networks to improve biomarker identification for predicting therapeutic responses. However, they also present some inherent limitations: while NetBio effectively analyzes biological interactions relevant to immune therapy, it lacks the capability to investigate biomarkers at the gene level, which limits its utility in elucidating ICI response mechanisms. Similarly, ICINet integrates diverse biological data to predict immune response using a GNN; however, it lacks transparency in identifying specific biological markers that drive these predictions.

To overcome the limitations of existing methods, we developed PathNetDRP, a novel biomarker identification framework that integrates PPIs and biological pathway data for improved ICI response prediction. Unlike conventional approaches, PathNetDRP applies the PageRank algorithm to identify ICI-response-associated pathways and incorporates statistical tests to select predictive biomarkers, enabling more effective classification of gene expression profiles into responder and non-responder groups. Specifically, PathNetDRP constructs pathway-specific subnetworks and applies PageRank to each subnetwork, allowing it to capture complex relationships between gene expression patterns and ICI response. This approach enables the identification of functionally relevant pathways and co-expression patterns of genes that improve ICI response prediction compared to existing biomarker selection methods. The major contributions of this study are as follows:

• Integration of Pathways, PPIs, and ICI Target Information: Unlike conventional ICI biomarker models, which rely primarily on gene expression differences or predefined immune signatures, PathNetDRP incorporates PPIs and ICI target information to identify biologically meaningful and functionally relevant biomarkers. By applying PageRank to individual pathways, our method provides a more precise, context-aware analysis of gene contributions to ICI response prediction.

• Comprehensive Performance Evaluation Across Multiple Cohorts: We conducted a rigorous evaluation of PathNetDRP across eight independent ICI-treated patient cohorts, demonstrating high predictive accuracy in drug response prediction compared to state-of-the-art (SOTA) methods. Our leave-one-out cross-validation (LOOCV) and independent validation tests further confirm the robustness and generalizability of our approach.

• Discovery of Novel ICI Response Biomarkers: PathNetDRP identified novel biomarker candidates that exhibited strong predictive performance compared to differentially expressed genes. Additionally, our analysis uncovered previously unrecognized biomarkers with high potential for distinguishing ICI responders from non-responders, contributing to a deeper understanding of ICI response mechanisms.

Methods

We designed a novel biomarker identification algorithm, PathNetDRP, for predicting ICI response using ICI-related pathways, PPI networks, and gene expression data. Figure 1 provides an overview of PathNetDRP, which initially identifies biomarker candidates using ICI target information and a network propagation algorithm. It then identifies ICI-response-related pathways using the candidates and a hypergeometric test, and finally, calculates PathNetGene scores and selects biomarkers using the identified pathways and PageRank algorithm. The pseudo-code for the PathNetDRP is provided in the supplementary Fig. S1.

Fig. 1

PathNetDRP Biomarker Identification Process. First, PageRank is applied to a PPI network to identify ICI-related genes. Next, ICI-related pathways are determined using these selected genes. Finally, the genes are evaluated within the identified pathways, PathNetGene scores are computed, and biomarkers are selected via Kolmogorov-Smirnov significance testing

PathNetDRP

The calculation of PathNetGene scores involves three steps: (1) ICI-target-based biomarker candidate selection using a PPI network and PageRank algorithm; (2) identification of ICI-related biological pathways through hypergeometric testing of the candidates; and (3) calculation of gene importance via network analysis on the identified pathways. Figure 2 illustrates the process of calculating the PathNetGene Score within PathNetDRP.

Fig. 2

Overview of the Three Key Processes in PathNetDRP: A Application of the PageRank algorithm on the interaction network to identify the top 2000 genes closely associated with ICI targets. B Identification of biological pathways significantly associated with the selected 2000 gene candidates. C Configuration of subnetworks for each identified pathway and calculation of PathNetGene scores, quantifying the relevance of each gene within the context of its pathway

ICI-related gene selection via PageRank

PathNetDRP begins by leveraging ICI target gene information and applying the PageRank algorithm to a PPI network to identify candidate genes associated with drug response (Fig. 2A). The underlying assumption is that genes neighboring ICI targets within the PPI network are likely to exhibit strong functional interactions and contribute to immune response mechanisms. To capture this biological relevance, ICI target genes are used to initialize gene scores, and their influence is then propagated across the PPI network via PageRank. This algorithm effectively prioritizes genes based on network connectivity and influence, enabling the identification of biologically significant genes that may play a crucial role in drug response. The PageRank algorithm iteratively updates gene scores based on the topology of the PPI network. For a given gene $g_i$, the gene score at iteration t is computed as follows:

$$\begin{aligned} PR(g_i; t) = \frac{1-d}{N} + d \sum _{g_j \in B(g_i)}{\frac{PR(g_j; t-1)}{L(g_j)}}, \end{aligned}$$

(1)

where N is the total number of genes in the network, d is the damping factor controlling the influence of neighboring nodes, $B(g_i)$ represents the set of genes linking to $g_i$, and $L(g_j)$ represents the number of genes linking from $g_j$. This computation was implemented using Python v3.8.18 and NetworkX v3.1. Initially, ICI target genes were assigned a score of 1, while all other genes were set to 0. After the PageRank algorithm converged, the top 2,000 genes were selected as biomarker candidates for further analysis.

ICI-related pathway identification via hypergeometric test

Since drug response mechanisms can be complex, with many genes cooperating across various biological pathways, PathNetDRP performs pathway enrichment analysis using statistical tests to find ICI-related pathways (Fig. 2B). This step ensures that PathNetDRP identifies biomarkers that are directly linked to immune-related biological processes, rather than relying solely on differential gene expression patterns, which may overlook functionally important interactions. To identify ICI-related pathways, we applied a hypergeometric test to assess the enrichment of selected genes within biological pathways. The analysis was performed on the 2,000 ICI-related genes identified in the previous step, and pathways with statistically significant enrichment ($p<0.01$) were retained for further gene scoring. For each pathway, the hypergeometric p-value was computed using the following formula:

$$\begin{aligned} P(X=k) = \frac{(\frac{K}{k})(\frac{N-K}{n-k})}{(\frac{N}{n})}, \end{aligned}$$

(2)

where K is the number of candidate genes, N is the total number of genes, n is the number of genes in the pathway, and k is the number of candidate genes within the pathway. To account for multiple hypothesis testing, p-values were adjusted using the Holm-Sidak method, and pathways with adjusted $p<0.01$ were selected as ICI-related pathways. The pathway database used for this analysis was MSigDB [20], a curated collection of molecular signatures representing biological processes. All statistical computations were performed using the Scipy v1.10.1 [21] and statsmodels v0.14.0 [22].

PathNetGene scoring based on ICI-related pathways

PathNetDRP introduces PathNetGene, a gene importance score that integrates biological pathway relevance with ICI response-associated interactions. This approach provides a quantitative assessment of each gene’s contribution within key immune-related pathways, leading to more precise biomarker identification (Fig. 2C). Specifically, we construct PPI subnetworks that include only genes within the previously identified ICI-related biological pathways. Within each subnetwork, we recalculate gene importance using the PageRank algorithm, where the initial PageRank score is set using each patient’s normalized gene expression levels (min-max scaling). This allows us to derive personalized gene importance scores, ensuring that gene contributions are evaluated based on both network topology and patient-specific gene expression data. Finally, we aggregate pathway-specific scores to compute the final PathNetGene score. Since each gene may be involved in multiple pathways, we calculate a weighted average of its PageRank scores across all associated pathways, adjusting weights based on pathway membership. Consequently, genes that hold high importance across multiple key ICI response pathways receive a higher final score. Through this process, the PathNetGene score serves as a biologically meaningful biomarker evaluation metric, capturing not only individual gene expression levels but also the functional role of genes within immune-related pathways. This integrative approach enhances the biological relevance of ICI response predictions, making PathNetDRP a more robust and interpretable framework for biomarker discovery. The PathNetGene scores of a gene $G_i$ for a patient $S_j$ is calculated by:

$$\begin{aligned} s_{ij} = w_{i} \sum _{P \in \mathcal {P}}{ \mathbb {1}_{G_i \in P} s_{ij}^{P}} \end{aligned}$$

(3)

where $\mathcal {P}$ is a set of biological pathways, $w_i$ is the weight of $G_i$, $\mathbb {1}_{G_i \in P}$ is an indicator function that is 1 if a pathway P contains $G_i$ otherwise 0, and $s_{ij}^{P}$ is a pathway-specific PageRank score of $G_i$ for $S_j$ on P. In this study, we defined the weight $w_i$ as the inverse of the number of pathways to which $G_i$ belongs, ensuring that its influence is adjusted accordingly. Specifically, $w_i$ is calculated as:

$$\begin{aligned} w_{i}=1/(|{P\in \mathcal {P}|G_i\in P}|) \end{aligned}$$

(4)

Identification of biomarkers for ICI response prediction

The proposed method identifies ICI-response-related biomarker genes by leveraging PathNetGene scores and applying the Kolmogorov-Smirnov (K-S) test [23]. A two-sample K-S test is conducted to determine whether there is a statistically significant difference in the distribution of PathNetGene scores between responder and non-responder groups. To control for multiple hypothesis testing, p-values are adjusted using the Holm-Sidak method, and the top 30 genes with $p$ $<$0.01 are identified as biomarkers of PathNetDRP for ICI response prediction. The number of biomarkers was heuristically determined in this study. All statistical analyses, including the K-S test and Holm-Sidak adjustment, were performed using SciPy v1.10.1.

ML-based predictive models

Various ML-based binary classifiers are employed to predict ICI response labels using biomarkers identified by PathNetDRP. We utilized four ML models: logistic regression (LR), SVM, RF, and multi-layer perceptron (MLP). These models take PathNetGene scores as input and predict whether a given sample belongs to the responder or non-responder group. MLP was implemented using PyTorch v2.2.1, while the other models were implemented using Scikit-learn v1.3.0.

• LR: A statistical method for binary classification that applies the logistic function to convert a linear combination of inputs into a probability, typically classifying the input based on a threshold.

• SVM: A robust classifier that separates data points in high-dimensional spaces by finding the hyperplane with the maximum margin, using support vectors to define this separation.

• RF: An ensemble classifier that improves predictive accuracy and controls overfitting by building multiple decision trees and aggregating their outcomes to determine the final class.

• MLP: A type of neural network with multiple layers of nodes, excelling in tasks requiring complex pattern recognition. It learns non-linear relationships through activation functions in each layer, making it particularly useful for image and speech processing.

• AdaBoost: An ensemble learning method that combines multiple weak learners to improve prediction accuracy. It iteratively adjusts the weights of misclassified samples, assigning them higher importance in subsequent models to focus the learning process on difficult cases.

Performance evaluation metrics

To evaluate the predictive performance of PathNetDRP and baseline models, we used three widely recognized metrics: Accuracy, F1 score, and area under the curve of a receiver operating characteristic (AUCROC). Each of these metrics provides unique insights into the model’s ability to classify ICI responders and non-responders effectively. Accuracy measures the proportion of correctly classified samples out of all predictions and is defined as:

$$\begin{aligned} Accuracy=\ \ (TP+TN)/(TP+TN+FP+FN) \end{aligned}$$

(5)

where TP, TN, FP, and FN represent true positives, true negatives, false positives, and false negatives, respectively. While accuracy provides an overall measure of correctness, it can be misleading when dealing with imbalanced datasets, where one class significantly outnumbers the other.

To mitigate this issue, we also employed the F1 score, which considers both precision and recall, making it particularly valuable for imbalanced datasets. The F1 score is the harmonic mean of precision and recall and is calculated as:

$$\begin{aligned} F1\ score=(2\times Precision\times Recall)/(Precision+Recall) \end{aligned}$$

(6)

where Precision represents the proportion of correctly predicted positive cases among all predicted positives, and Recall reflects the proportion of correctly predicted positive cases among all actual positives. A higher F1 score indicates a well-balanced model that minimizes both false positives and false negatives, which is crucial in clinical applications where misclassification can have significant consequences. In addition to these metrics, we used AUCROC, which evaluates the model’s ability to distinguish between responders and non-responders by analyzing the trade-off between sensitivity and specificity across different classification thresholds. The AUCROC measures the area under the receiver operating characteristic curve, where a score of 0.5 suggests the model performs no better than random guessing, while a score of 1.0 indicates perfect classification. This metric is particularly useful for comparing models with varying decision thresholds, providing a comprehensive measure of discrimination performance. All metric calculations were performed using Scikit-learn v1.3.0.

Results

Materials and experimental details

Eight ICI cohort datasets

We collected eight cohort datasets of patients treated with ICIs for melanoma, gastric cancer, and bladder cancer. These datasets include pre-treatment gene expression data, drug response labels, and patient survival data, making them highly relevant for studying ICI response prediction. The Gide, Liu, Kim, and Huang datasets were sourced from the data repository of Lee et al. [24] and can be accessed at Zenodo (https://zenodo.org/records/4661265). Data preprocessing for these datasets was conducted using R version 4.4.3. The IMvigor210 dataset was obtained through the IMvigor210CoreBiologies package [25]. This dataset was downloaded using R 3.4 and Bioconductor 3.6, following the guidelines and required libraries available at IMvigor210CoreBiologies. The Auslander and Riaz datasets, also known as GSE115821 and GSE91061, were downloaded from the GEO database [16, 26]. Additionally, the Prat dataset for melanoma was obtained from the supplementary materials of a prior study [27]. A detailed description of these eight datasets is provided in Table 1.

Table 1 Descriptions of the Eight ICI Cohort Datasets

ICI response labels

We labeled cancer patients based on the effectiveness of ICI response for each dataset. For the Auslander and Kim datasets, we used label data provided in the original papers [16, 28]. For the Huang dataset, we classified patients without cancer recurrence as ICI responders, while others were categorized as non-responders. For the remaining datasets, ICI response labels were determined based on the Response Evaluation Criteria in Solid Tumors (RECIST), which classifies response into four categories: complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD) [29]. The details of RECIST were provided in Supplementary Table S1. Following Kong et al. [18], patients labeled with CR or PR were categorized as the ICI responder group, while those with SD or PD were categorized as non-response.

PPI network and pathways

We utilized human protein network data downloaded from the STRING database (v11.0) with a confidence score greater than 700 [30]. This yielded 841,068 interactions and 17,012 gene symbols. We also sourced Reactome pathway data from MSigDB, comprising 1692 pathways and 11,155 gene symbols [20, 31].

ICI target genes

The benchmark datasets used in this study include four ICIs: nivolumab, pembrolizumab, ipilimumab, and atezolizumab. The target gene information for these drugs was obtained from DrugBank (Supplementary Table S2).

Model evaluation

Various cross-validation methods were employed to optimize hyperparameters and evaluate the performance of ML classifiers based on PathNetDRP biomarkers. We used 5-fold cross-validation and a grid search algorithm for model parameter optimization. Supplementary Table S3 lists the model parameters for LR, SVM, RF, MLP and AdaBoost. The regularization strength of LR and SVM, the number of trees in RF, the number of hidden layers in MLP, and the number of epochs for MLP were selected as hyperparameters and optimized using the area under the curve of receiver operating characteristics (AUCROC).

Cross-validation

LOOCV is used to evaluate accuracy and F1 score on test data (Supplementary Fig. S2), where accuracy measures the ratio of correctly predicted labels (i.e., true positive and true negative), and F1 score, the harmonic mean of precision and recall, measures the accuracy of true positive predictions. In these evaluations, responder and non-responder groups are designated as positive and negative, respectively.

Independent validation

To validate identified biomarkers, we conducted independent testing where training and test datasets were constructed from different datasets. We used five melanoma datasets including Gide, Auslander, Prat, Riaz, and Huang. The Gide dataset served as the training dataset, while the others were used as test datasets. Supplementary Fig. S3 illustrates the process of independent testing.

Experimental results

Model selection on PathNetDRP biomarkers

To effectively predict ICI response using the selected biomarkers, we compared four ML models: LR, RF, SVM, MLP, and AdaBoost. Supplementary Table S4 displays the accuracy and F1 scores of these models across four benchmark datasets: Gide, Liu, Kim, and IMvigor210. All predictive performances were evaluated using LOOCV, and the models were optimized through grid search. SVM and RF demonstrated suboptimal F1 scores on the Kim and IMvigor210 datasets, respectively, while LR consistently achieved better predictive performance than the other models. Consequently, we selected LR as the most suitable classifier and subsequently used it to compare the biomarkers of PathNetDRP with various baseline models.

We also examined the optimal number of biomarkers identified by PathNetDRP using varying gene set sizes. Genes were ranked by p-values, and F1 scores were evaluated using the top 10, 30, 50, 70, 100, and all genes. F1 scores were computed using LOOCV and LR. Our analysis determined that 30 genes provided the optimal biomarker set (Supplementary Table S5). Consequently, all experimental results in this study were obtained using 30 genes.

Comparison with known ICI response biomarkers

To assess the predictive performance of ICI response using PathNetDRP-identified biomarkers, we conducted comparisons with established ICI-related biomarkers. Employing the methodology described by Kong et al. [18], we utilized PD-1, PD-L1, CTLA4, CD8 T cells (CD8T) [32], T-cell exhaustion (Tex) [15], TAM [15], and CAF [33] as control group biomarkers (Supplementary Table S6). LR, trained on the gene expression levels of these biomarkers, was used as the prediction model. Table 2 shows the accuracy and F1 scores for four benchmark datasets: Gide, Liu, Kim, and IMvigor210. The biomarkers identified by PathNetDRP achieved the best accuracy scores of 0.747, 0.748, 0.956, and 0.735, respectively, and the highest F1 scores of 0.768, 0.700, 0.917, and 0.521, respectively, on these datasets. The fluctuations in performance can be attributed to differences in treatment regimens and cancer types. Despite these variations, PathNetDRP consistently showed better predictive performances than known ICI response biomarkers in terms of accuracy and F1 scores, with the exception of the F1 score in the Gide dataset. These results suggest that PathNetDRP effectively captures immune response mechanisms across different types of cancer while being influenced by treatment heterogeneity.

Table 2 Comparison of LOOCV Results on Four Benchmark Datasets

Survival analysis was conducted on three datasets providing overall survival data: Gide, Liu, and IMvigor210. Patients were divided into two groups based on the predicted response labels using biomarkers identified by PathNetDRP and control group biomarkers, and Kaplan-Meier curves and log-rank tests were calculated (Fig. 3 and Supplementary Fig. S4). Unlike the control group biomarkers in the Liu dataset, which showed higher survival rates in the non-responder group, PathNetDRP’s results aligned normally with expected outcomes. Across the Gide, Liu, and IMvigor210 datasets, our experimental results confirmed that the predicted responder group exhibited a low hazard risk. Specifically, the hazard ratio (HR) scores were 0.32 for Gide, 0.35 for Liu, and 0.48 for IMvigor210, demonstrating a consistent trend across different datasets.

Fig. 3

Kaplan–Meier curves comparison between PathNetDRP and ICI drug targets. A Gide, B Liu, and C IMvigor210

Benchmarking against state-of-the-art biomarker identification methods

PathNetDRP demonstrated improved predictive performance in predicting ICI response compared to SOTA methods for identifying immune response biomarkers. The same four benchmark datasets—Gide, Liu, Kim, and IMvigor210—were used to perform LOOCV. Control group biomarkers were identified using four SOTA methods: NetBio [18], TIDE [15], IMPRES [16], and DeepGeneX [17]. DeepGeneX was excluded from the gastric cancer and bladder cancer datasets, Kim and IMvigor210 respectively, as it is specifically tailored for melanoma [17]. PathNetDRP showed competitive performance when compared to SOTA methods, based on receiver operating characteristic (ROC) curves and precision-recall (PR) curves (Fig. 4, Supplementary Fig. S5, and Supplementary Table S7). PathNetDRP achieved AUCROC scores of 0.820, 0.801, 0.934, and 0.784 on the Gide, Liu, Kim, and IMvigor210 datasets, respectively, demonstrating better predictive discrimination capability between ICI responder and non-responder groups compared to SOTA methods. Additionally, for the area under curve of PR (AUCPR) scores (a.k.a. average precision score), PathNetDRP recorded 0.841, 0.721, 0.936, and 0.532, respectively, demonstrating that our biomarkers were more effective in identifying true responder groups than the baseline methods.

Fig. 4

Receiver operating characteristic curves and areas under curves. A Gide, B Liu, C Kim, and D IMvigor210

We performed Kaplan-Meier survival analysis and log-rank tests to assess whether patients with higher drug response exhibited greater survival rates (Fig. 5 and Supplementary Fig. S6). The log-rank test for PathNetDRP revealed statistically significant differences in survival across the Gide, Liu, and IMvigor210 cohorts ($p<1e-3$; Fig. 5A–C). Additionally, patients classified as responders demonstrated higher survival rates compared to non-responders. In contrast, some cohorts analyzed using traditional methods did not yield significant p-values, highlighting limitations in their predictive accuracy. Furthermore, our method achieved the highest c-index values, indicating improved performance relative to conventional methods (Fig. 5D–F). These findings suggest that PathNetDRP provides a more reliable survival stratification, reinforcing its potential clinical utility in predicting ICI response.

We also evaluated the predictive performance of PathNetDRP using disease control rate (DCR), where PR/CR/SD patients were classified as responders, and confirmed that PathNetDRP maintained strong predictive performance (Supplementary Table S8). Across multiple cohorts, PathNetDRP consistently outperformed SOTA methods in accuracy and F1 score, further supporting the robustness of PathNetDRP-derived biomarkers. These findings highlight the reliability of PathNetDRP in predicting ICI response using DCR-based classification.

Fig. 5

A–C Log-rank test comparison with baseline biomarker identification methods. A Gide, B Liu, and C IMvigor210 datasets. D–F c-index Results Comparison with Baseline Biomarker Identification Methods. D Gide, E Liu, and F IMvigor210 datasets

Independent validation

To evaluate the generality of biomarkers identified by PathNetDRP, we conducted an independent test, which assesses whether a classifier constructed on given biomarkers is well performed on other datasets that were not part of the model construction process. Specifically, both biomarker identification and classifier training are conducted using whole data in the Gide dataset, and then four datasets-Auslander, Prat, Riaz, and Huang-were used as test datasets to validate the robustness and reliability of the biomarkers. In this experiment, we utilized AdaBoost as a classifier because it had good predictive power for independent test (Supplementary Table S9).

To ensure robust results, experiments were repeated 10 times, and final outcomes were reported as the average of these repetitions. We used AUCROC as the evaluation metric to compare the predictive performance of PathNetDRP with baseline methods, including TIDE, IMPRES, NetBio, DeepGeneX, and IRNet. DeepGeneX and IRNet are deep learning frameworks designed to predict ICI treatment response based on gene expression patterns, while NetBio is a network-based machine learning method that identifies ICI-associated biomarkers and predicts immunotherapy response. TIDE and IMPRES are statistical biomarker selection methods that evaluate gene-ICI response relationships. For PathNetDRP, TIDE, and IMPRES, we used an AdaBoost classifier to predict ICI response labels based on the identified biomarkers. Since raw count data required for IRNet were unavailable in the Huang dataset, its results were not included.

As shown in Table 3 and Supplementary Fig. S7, PathNetDRP achieved AUCROC scores of 0.779±0.023, 0.701±0.009, 0.655±0.016, and 0.893±0.023 on the Auslander, Prat, Riaz, and Huang datasets, respectively, showing competitive performance compared to baseline methods. PathNetDRP achieved the highest AUCROC scores in the Auslander, Prat, and Huang datasets. Notably, even in the drug relapse-related Huang dataset, PathNetDRP achieved an AUCROC of 0.893, suggesting that our framework can effectively identify biomarkers associated with immunotherapy response. Although PathNetDRP did not have the best score on the Riaz dataset, the score of 0.655±0.016 was significantly higher than the baseline average of 0.625, indicating that PathNetDRP demonstrates moderate or better accuracy in ICI response prediction.

Table 3 Comparison of Independent test results to all baseline methods

Differential expression analysis of PathNetDRP biomarkers

To explore the association between biomarkers identified by PathNetDRP and ICI response, we conducted differential gene expression (DGE) analysis on all the biomarker genes obtained from the Gide, Kim, Liu, and IMvigor210 datasets (Fig. 6 and Supplementary Tables S10–S13), investigating their biological functions related to ICI response. We applied the Mann–Whitney U test for DGE analysis. In the Gide dataset, genes such as HLA-A, HLA-B, and HLA-F were identified as biomarkers (Fig. 6A–C). These genes encode human leukocyte antigen class I (HLA class I) proteins, which present antigens to T cells and serve as ligands for various immune receptors on natural killer (NK) cells, T cells, and bone marrow cells [34, 35]. The gene expression for HLA-A, HLA-B, and HLA-F were statistically significantly higher in the responder group ($p<0.05$), suggesting that the expression patterns of these HLA class I proteins are significant biomarkers for predicting the response to ICI inhibitors.

Fig. 6

Differential gene expression analysis on PathNetDRP miomarkers. A–C Biomarkers identified on Gide dataset, D Liu, C Kim, and D IMvigor210

In the Liu cohort, IFNLR1 was identified as a biomarker. Lasfar et al. [36] reported that IFNLR1 regulates IFN-$\lambda$ and its antitumor effect in melanoma is dose-dependent. Sato et al. [37] also found that overexpression of IFN-$\lambda$ significantly reduces lung metastasis in melanoma. The gene expression for IFNLR1 was higher in the responder group ($p=5.26e$-02) (Fig. 6D), indicating its potential as a valuable biomarker for predicting ICI response.

In the Kim dataset, HLA-G was identified as a biomarker. Recently studied as a novel immune checkpoint, HLA-G is known to interact with immune receptors and inhibit the biological functions of immune cells [38, 39]. According to Bartolome et al. [40], higher expression of HLA-G correlates with lower survival rates in solid tumors. The gene expression for HLA-G was statistically higher in the non-responder group ($p=1.38e$-04) (Fig. 6E), underscoring the robust biomarker identification capabilities of PathNetDRP.

In the IMvigor210 dataset, TGF-$\beta$1 was identified as a biomarker. TGF-$\beta$ proteins, which are highly expressed in advanced cancers, are associated with poor prognosis [41]. Ni et al. [42] reported that patients with low TGF-$\beta$ scores in gynecological cancer had significantly improved proliferation-free survival. The gene expression for TGF-$\beta$1 was lower in the responder group in the bladder cancer cohort of IMvigor210 ($p=6.31e$-05) (Fig. 6F), demonstrating that TGF-$\beta$1 could also serve as a crucial ICI response prediction biomarker for bladder cancer.

Beyond differential gene expression: advantages of PathNetDRP

To elucidate the biological significance of biomarkers identified by PathNetDRP, we compared them with those identified using Limma-Trend [43], a widely used tool for differential gene expression analysis. We selected 34 biomarker genes that were significantly differentially expressed in the Gide cohort using Trend (Supplementary Table S14) and evaluated their predictive performance via LOOCV. Trend achieved an AUCROC of 0.808 (Supplementary Fig. S8), which is lower than the PathNetDRP performance reported in Fig. 4A. Furthermore, independent validation on melanoma datasets confirmed that Trend biomarkers exhibited lower generalizability compared to those identified by PathNetDRP (Supplementary Fig. S9). Additionally, enrichment analysis using g:Profiler [44] (Supplementary Tables S15 and S16) revealed that PathNetDRP biomarkers are more closely associated with regulation of immune response (GO:0050776), whereas Trend biomarkers are associated with immune response (GO:0006955). Although Trend effectively identifies differentially expressed genes between distinct immune response groups, PathNetDRP integrates pathway and PPI information, allowing it to identify both genes specific to each immune response group and their biologically associated neighbors. This enables PathNetDRP to achieve more accurate ICI response predictions by capturing genes involved in immune response regulation.

Enrichment analysis on PathNetDRP biomarkers

To investigate the biological functionality of biomarkers derived from PathNetDRP, we performed enrichment analysis on biomarkers identified from the Gide dataset using both DAVID [45] and g:Profiler[44] (Supplementary Table S17). The analysis identified 10 commonly enriched pathways, including four immunity-associated pathways: interferon alpha/beta signaling (R-HSA-909733), interferon gamma signaling (R-HSA-877300), adaptive immune system (R-HSA-1280218), and interleukin-10 signaling (R-HSA-6783783). These pathways were visualized using Cytoscape (Fig. 7A). While most PathNetDRP biomarker genes are included in these four Reactome pathways, known ICI response-related biomarkers and ICI target genes are not part of these pathways. This observation suggests that PathNetDRP biomarker genes not included in these pathways might play crucial roles in predicting ICI response. Among such genes, IL2RB demonstrated interactions with all four identified pathways, prompting further analysis.

Fig. 7

Functional interaction network with enriched reactome pathways on PathNetDRP biomarkers from the gide dataset. A Visualization of enriched Reactome pathways. B Comparison of IL2RB expression levels between responder (blue) and non-responder groups using a t-test. C Chi-square test assessing the association between IL2RB and HLA-A co-expression patterns and ICI response groups. High and low expression levels were defined based on mean values. D Kaplan-Meier survival analysis based on IL2RB and HLA-A expression levels

IL2RB is activated through interactions with various pathway-associated genes, such as CD8A, CD8B, and HLA-A. Serving as the beta subunit of the Interleukin-2 (IL-2) receptor, IL2RB mediates IL-2 signaling to promote the activation, proliferation, and survival of T cells and NK cells. Identified as a growth factor for T cells, IL-2 is a key cytokine affecting the immune system and was an early candidate for cancer immunotherapy, approved in 1998 for treating metastatic melanoma [46]. Recently, complementary therapeutic approaches combining ICIs and IL-2 immunotherapy are being explored, with ongoing active research [47]. These studies suggest that IL2RB, closely related to IL-2, may be involved in the efficacy of IL-2 based immunotherapies.

For a deeper analysis, we examined the differential expression patterns of IL2RB between ICI response groups (Fig. 7B). A one-sample t-test confirmed that IL2RB was significantly upregulated in the responder group (p = 5.18e$-$07). Additionally, we found that the co-expression patterns of IL2RB and HLA-A may serve as potential biomarkers for ICI therapy prognosis. As shown in Fig. 7C and Supplementary Fig. S10, patients with simultaneous upregulation of IL2RB and HLA-A tended to exhibit better ICI responses compared to other groups. We statistically evaluated this observation using chi-square test, confirming its significance ($p<0.01$).

Furthermore, Kaplan-Meier survival analysis supported the prognostic relevance of IL2RB and HLA-A co-expression (Fig. 7D). The log-rank test comparing the high/high expression group against the remaining groups showed a statistically significant survival difference (p = 7.62e$-$03), suggesting the predictive potential of IL2RB and HLA-A as biomarkers for ICI response.

Uncovering novel ICI response biomarkers with PathNetDRP

Finally, we investigated novel biomarkers for practical ICI response prediction among the genes identified by PathNetDRP. To assess the novelty of these biomarkers, we compared biomarkers selected by PathNetDRP in the Gide cohort with known biomarkers (Supplementary Table S6) and categorized them into overlapping and non-overlapping genes (Supplementary Fig. S11). We then evaluated the predictive performance of these gene sets and confirmed that non-overlapping genes, which were exclusively identified by PathNetDRP, exhibited better predictive performance than overlapping genes, which consist of widely studied biomarkers (Supplementary Table S18). These findings suggest that PathNetDRP can identify novel gene signatures relevant to ICI response prediction.

Furthermore, we conducted functional analysis on the non-overlapping genes and identified potential biomarkers, including HLA-F, GBP1, and IFI16. While HLA-F has not yet been fully established as a predictive biomarker, it plays a role in immune modulation through interactions with NK cell receptors [48]. Similarly, GBP1 and IFI16 have been relatively understudied as ICI response biomarkers, yet they may serve as promising research targets for ICI response prediction (Supplementary Fig. S12). We found that both genes were significantly upregulated in responders. Moreover, when samples were stratified into two groups based on mean expression levels of these genes, the high-expression group demonstrated better survival outcomes in the Kaplan-Meier survival analysis.

Discussion

In the functional analysis of PathNetDRP biomarkers, DGE analysis was conducted across various datasets, highlighting significant genes such as IL2RB, HLA-A, HLA-B, HLA-F, and HLA-G. These genes play crucial roles in the immune response by serving as ligands for immune receptors. The PathNetGene scores for these genes were significantly higher in the responder groups, indicating their potential as predictive biomarkers for ICI response. Moreover, enrichment analysis revealed these biomarkers’ association with key Reactome pathways, such as interferon alpha/beta signaling and interleukin-10 signaling, which are crucial for immune response mechanisms. This suggests that PathNetDRP not only identifies relevant biomarkers but also connects them to biological pathways that could aid in understanding the underlying mechanisms of ICI response, enhancing the potential for targeted immunotherapy strategies.

The primary limitation of PathNetDRP lies in its dependence on large and diverse datasets for effective training and validation, which may not always be readily available. While network-based approaches provide valuable insights into gene interactions, PathNetDRP may not fully capture the dynamic nature of the tumor microenvironment and the complex interplay among cellular components. This limitation could affect the generalizability of its predictions across different cancer types and treatment conditions, particularly in cases where tumor heterogeneity plays a significant role. Additionally, although tumor mutational burden (TMB) and microsatellite instability (MSI) are well-established as important predictors of ICI response, our study is currently limited to gene expression-based analysis. The exclusion of these additional predictive factors represents a key constraint, as integrating multi-omics data could further enhance the robustness and accuracy of PathNetDRP.

Conclusion

PathNetDRP, our novel biomarker discovery framework that utilizes pathway and PPI networks, has demonstrated significant potential in predicting ICI responses. This study has effectively showcased the efficacy of integrating network-based methodologies with ML to enhance the precision of ICI response predictions across various cancer types. Through extensive validation and comparative analyses, PathNetDRP provided not only improved predictive accuracy but also insights into the biological processes influencing ICI response.

Beyond its application in ICI response prediction, PathNetDRP has the potential to be integrated into clinical workflows for personalized immunotherapy. By predicting drug response before treatment initiation, the framework could assist clinicians in patient stratification, optimizing treatment decisions, and enhancing therapeutic outcomes. Future studies will focus on the clinical validation of PathNetDRP, exploring how its biomarker predictions can be incorporated into decision-making processes for immunotherapy.

While this study focused on ICI response prediction, PathNetDRP is not limited to ICI applications. The framework can be extended to predict drug responses for other therapeutic agents, provided that drug target information, PPIs, pathways, and gene expression profiles are available. In the future work, we will focus on expanding PathNetDRP to include more comprehensive multi-omics data, such as metabolomics and proteomics, to enhance prediction accuracy. We also plan to develop methods that can accommodate smaller or less diverse datasets through advanced ML techniques, such as transfer learning or few-shot learning, to improve the applicability of PathNetDRP in less-studied cancers or personalized medicine scenarios.