Integration of graph neural networks and transcriptomics analysis identify key pathways and gene signature for immunotherapy response and prognosis of skin melanoma

Abstract

Objective

The assessment of immunotherapy plays a pivotal role in the clinical management of skin melanoma. Graph neural networks (GNNs), alongside other deep learning algorithms and bioinformatics approaches, have demonstrated substantial promise in advancing cancer diagnosis and treatment strategies.

Methods

GNNs models were developed to predict the response to immunotherapy and to pinpoint key pathways. Utilizing the genes from these key pathways, multi-omics bioinformatics methods were employed to refine the construction of a gene signature, termed responseScore, aimed at enhancing the precision of immunotherapy response predictions. Subsequently, responseScore was explored from the perspectives of prognosis, genetic variation, pathway enrichment, and the tumor microenvironment. Concurrently, the association among 13 genes contributing to responseScore and factors such as immunotherapy response, prognosis, and the tumor microenvironment was investigated. Among these genes, PSMB6 was subjected to an in-depth analysis of its biological effect through experimental approaches like transfection and co-culture.

Results

In the finalized model utilizing GNNs, it has revealed an AUC of 0.854 within the training dataset and 0.824 within the testing set, pinpointing key pathways such as R-HSA-70,268. The indicator named as responseScore excelled in its predictive accuracy regarding immunotherapy response and patient prognosis. Investigations into genetic variation, pathway enrichment, tumor microenvironment disclosed a profound association between responseScore and the enhancement of immune cell infiltration and anti-tumor immunity. A negative correlation was observed between the expression of PSMB6 and immune genes, with elevated PSMB6 expression correlating with poor prognosis. ELISA detection after co-cultivation experiments revealed significant reductions in the levels of cytokines IL-6 and IL-1β in specimens from the PCDH-PSMB6 group.

Conclusion

The GNNs prediction model and the responseScore developed in this research effectively indicate the immunotherapy response and prognosis for patients with skin melanoma. Additionally, responseScore provides insights into the tumor microenvironment and the characteristics of tumor immunity of melanoma. Thirteen genes identified in this study show promise as potential tumor markers or therapeutic targets. Notably, PSMB6 emerges as a potential therapeutic target for skin melanoma, where its elevated expression exhibits an inhibitory effect on the tumor immunity.

Introduction

Skin melanoma, originating from melanocytes within skin tissue, accounts for over 90% of skin tumor-related fatalities [1,2,3]. Presently, a variety of treatment strategies have been developed for skin melanoma, including surgical resection [4, 5], chemotherapy [6, 7], targeted therapy (BRAF and MEK inhibitors) [8,9,10], and immunotherapy [11,12,13]. The effectiveness of these treatments for skin melanoma is influenced by clinical and pathological factors such as tumor staging, which categorizes the severity of melanoma into four stages based on tumor thickness, lymph node involvement, and metastasis [14], serving as an important prognostic indicator. Yet, these clinical and pathological characteristics often fall short in accurately reflecting patient prognosis which show considerable differences in the same stage due to melanoma’s heterogeneity. Thus, there is a pressing need to identify more reliable biomarkers for risk stratification and clinical management of skin melanoma. Furthermore, the recent advancements in immunotherapy, particularly immune checkpoint inhibitors like CTLA-4 and PD-1 inhibitors, have ushered in a new era in skin melanoma treatment over the last decade. These inhibitors activate the adaptive immune system, enhance immune surveillance, and facilitate the recognition and destruction of tumors by the host, significantly elevating melanoma survival rates [15]. Nonetheless, immunotherapy does not prove effective for all patients, with varied responses to treatments like PD-1 inhibitors [16]. Considering the substantial costs and potential permanent life-changing side effects associated with immune checkpoint inhibitors [17], it becomes imperative to assess the therapeutic impact of immunotherapy on individual patients, weighing the risks and benefits to make optimal treatment decisions. Precise pre-treatment evaluation of immunotherapy’s therapeutic effect can offer a solid basis for clinicians’ decisions, enabling tailored treatment approaches. Identifying biomarkers capable of assessing the therapeutic effect of immunotherapy has emerged as a critical area of focus in skin melanoma research.

Currently, pathway-based gene expression profiling has demonstrated notable stability and excellent performance. Pathways, as topological features, encapsulate the inherent relationships among various genes to a certain extent. Consequently, pathway-based gene expression profiling is more adept at identifying biomarker genes associated with disease phenotypes compared to conventional gene expression profiling. Graph neural networks (GNNs) have garnered significant interest for their proficiency in discerning the hidden topological information in pathways [18]. As an advanced deep learning method, GNNs have increasingly found application in contemporary medical biology research. More and more cancer researchers are recognizing the value of GNNs in cancer research [19]. For example, Peng W and colleagues [20] developed a GNNs model for oncogene prediction, while Isaac Furtney et al. [21] used GNNs to subtype breast cancer. Rui Yin et al. [22] applied Graph Neural Networks to identify potential miRNA targets in colorectal cancer. The greatest advantage of applying GNNs in transcriptomics lies in their ability to reflect the overall impact of specific pathways on prediction outcomes, i.e., their importance. This interpretability in pathways makes GNNs an essential tool for pathway-based transcriptomic analysis. Based on the superior performance of GNNs in bioinformatics analysis and their interpretability in pathways, we have decided to explore the application value of GNNs in transcriptomics and predicting immunotherapy responses. In this study, GNNs was utilized to extract transcriptomic and pathway-related topological features to model and predict skin melanoma patients’ response to immunotherapy and to identify key pathways influencing immunotherapy response.

Based on the key pathways identified through GNNs, we advanced the construction of a gene signature employing multi-omics bioinformatics approaches for a more precise prediction of immunotherapy response. Furthermore, significant genes influencing immunotherapy response were identified, and their functional impact on tumor immunity was investigated. The objective of this research is to develop an improved gene signature for predicting patient immunotherapy response and prognosis, explore the mechanisms underlying this signature, and identify novel biomarkers to guide the diagnosis and targeted treatment of skin melanoma.

Materials and methods

Data acquisition

The melanoma patient cohort was obtained from the TIGER database (http://tiger.canceromics.org/, accessed on April 20, 2023). The dataset includes clinical information and gene expression data for 342 metastatic melanoma patients who received immunotherapy, including 136 responders and 206 non-responders. All tumor samples from patients were collected before receiving anti-PD1 immune checkpoint therapy (nivolumab or pembrolizumab). According to RECIST1.1 criteria [23], patients with complete or partial response were classified as responders, while those with stable disease or disease progression were classified as non-responders. The corresponding data can be found at https://www.jianguoyun.com/p/DU70LLUQu86UDRiL1uUFIAA. Data on gene expression and related clinical information for TCGA-SKCM were downloaded from The Cancer Genome Atlas Program (TCGA) database (accessed on April 25, 2023, containing information on 460 patients). RNA sequencing data were normalized to transcripts per thousand base million (TPM). Gene mutation data were acquired from the Genomic Data Commons (GDC; https://portal.gdc.cancer.gov/) database and analyzed using the ‘maftools’ package (PMID: 30341162).

Construction of GNNs models and identification of key pathways

We began by extracting pathway and gene network data from the Reactome database via the ‘graphite’ package (version 1.24), excluding pathways that had either less than 15 or more than 400 genes, as well as those that were duplicates, leaving us with 927 pathways for our first GNN model construction. Our study employed a GNNs framework inspired by PathGNN, as developed by Bilin Liang et al. [24], focusing on predicting the survival rates of cancer patients over a long term. The design of our model is detailed in Supplementary Fig. 1. The GNNs model for prediction is composed of two parts, Subnetwork1 and Subnetwork2. Subnetwork1 captures gene information and features hidden in pathways, while Subnetwork2 divides patients into Response and non-Response groups. All pathway graphs are input into Subnetwork1. Subnetwork1 consists of three blocks, a readout layer (using Set2Set al.gorithm) after each block, a graph normalization layer between blocks, and a multilayer perceptron (MLP). Each block includes a graph convolution (GraphSAGE) layer and a graph pooling (SAGPool) layer to learn the feature representation of pathways. Subnetwork2 is based on MLP and consists of an input layer upstream to two hidden layers with 50 and 10 nodes respectively, a dropout layer with a dropout rate of 40%, and an output layer downstream. A unit (ReLU) activation function follows each hidden layer. Finally, cross-entropy loss was used to train this model. Detailed explanations of the model’s mathematical foundations are available in Supplementary Table 1.

Pytorch (version 1.8) and Pytorch Geometric (version 2.0.3) were used for GNN model construction and construction of the training and testing framework. The melanoma patient cohort from the TIGER database was randomly divided into 80% for training and 20% for testing. This model was trained to predict responses to immunotherapy over 150 epochs, using the Adam optimization algorithm with a learning rate of 0.001 for the first 100 epochs, then reducing it to 0.0005 for the remainder. The importance of each pathway in prediction outcomes was analyzed using the Integrated Gradients (IG) algorithm from the ‘Captum’ module (https://github.com/pytorch/captum.git). Pathways with a z-transforming IG score above 4, were selected for constructing a more streamlined version of GNNs model. This model underwent training for 250 epochs, with learning rates adjusted to 0.001 for the initial 150 epochs and 0.0005 for the last 100. The model was uploaded to GitHub (https://github.com/maodongYe/Immunotherapy-Response-of-Skin-Melanoma). Model’s performance was assessed with the area under the curve (AUC) metric.

Screening of significant genes

The process began with the extraction of genes within pathways whose z-transforming IG score above 4. Following this, genes with a significant correlation to immunotherapy response were discerned through univariate logistic regression, adhering to a significance level of p < 0.05. This approach resulted in the selection of 346 genes for inclusion in subsequent analysis.

Clustering analysis

The “ConsensusClusterPlus” R package (PMID: 20427518) facilitated the implementation of consensus clustering, dividing patients into different subgroups according to gene expression matrix. Additionally, principal component analysis (PCA) was performed to evaluate the clustering effect.

Development and validation of the gene signature for a more concise immunotherapy response prediction

The set of 346 significant genes was used in the least absolute shrinkage and selection operator (LASSO) regression analysis. LASSO regression can further identify feature genes related to the target variable and provide regression coefficients for these feature genes relative to the target variable under different error conditions. Utilizing the melanoma patient cohort dataset, this research used the LASSO regression algorithm from the ‘glmnet’ R package (version 4.1-7) to identify genes crucial for immunotherapy response. By weighting expression values of these significant genes with their LASSO coefficients, the gene signature named responseScore was developed for a more concise immunotherapy response prediction. The TCGA-SKCM dataset served as an external validation dataset to ascertain the performance of responseScore.

Survival analysis and pathway enrichment

The melanoma patient cohort was divided into high and low groups based on responseScore using the median as the boundary. The ‘survminer’ R package (version 0.4.9) was employed to plot Kaplan-Meier (K-M) survival curves, offering insights into survival trends over specified periods. A multivariate Cox proportional hazards analysis was conducted, with responseScore and a range of clinical factors as variables. Based on results of this Cox analysis, the nomogram predicting overall survival (OS) time for 1, 3, and 5 years was constructed using R package ‘rms’ (version 6.9-0). Based on differentially expressed genes between high and low responseScore group, Gene Ontology (GO) and Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed using ‘clusterProfiler’ (version 4.4.4) [25] package.

Tumor immune-related analysis

To calculate the Immune Score and Stromal Score of patient samples, the ‘estimate’ package (PMID: 24113773) in R was utilized. Similarly, the ‘xCell’ package (PMID: 29141660) facilitated the assessment of immune cell infiltration within patient samples. Bagaev et al. proposed 29 functional gene expression signatures (Fges) to encapsulate the cells and functional properties of the tumor microenvironment [26]. The ‘PreMSIm’ R package [27] was employed to determine microsatellite instability (MSI) status for each specimen, categorizing them into MSI-high (MSI-H) and MSI-low/microsatellite stability (MSS). Tumor Immune Dysfunction and Exclusion [26] (TIDE, accessible at http://tide.dfci.harvard.edu/login/) provided metrics such as Exclusion score, etc.

Functional verification experiment of PSMB6

Two melanoma cell lines (A375 and A875) were cultured. Empty vector and PSMB6 overexpression vector were constructed and transfected into melanoma cell lines (A375, A875) to obtain PCDH-NC group and PCDH-PSMB6 group, respectively. PSMB6 expression was detected by qPCR 24 h after transfection and by Western Blot (WB) 48 h after transfection. CD8 + T cells were sorted from human blood samples (10 ml whole blood) using magnetic bead sorting kit, and CD8 was identified by flow cytometry. After 24 h of transfection, the sorted CD8 + T cells were mixed with PCDH-NC group and PCDH-PSMB6 group of two melanoma cell lines (A375, A875) for co-culture (ratio: adherent: suspension = 1:5), and the cell supernatant was collected after 48 h. ELISA was used to detect the cytokines released by T cells in the supernatant (tumor necrosis factor α (TNF-α), interferon γ (IFN-γ), IL-6, IL-1β). The ELISA experiment was repeated three times. For detailed experimental procedures, please refer to Supplementary Table 2.

Statistical analysis

The codes of the bioinformatics analysis were uploaded to Jianguoyun repository (https://www.jianguoyun.com/p/DRYrwBAQlc6UDRiK1uUFIAA). For the purpose of comparing continuous variable between various groups, researchers employ the Wilcoxon test. Meanwhile, the analysis of categorical variable across different groups is conducted using the Chi-square test. To examine inter-group difference of survival rate, the Log-Rank test is frequently utilized. Delong’s test is used for AUC comparison. The analysis of correlation between variables is made using either the Pearson or Spearman correlation coefficient, depending on the nature of the data distribution. A p-value below 0.05 is regarded as indicative of statistical significance, with significance levels further delineated by the number of asterisks (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001).

Results

Construction of GNNs models and development of the gene signature response score

The main flowchart of this study is depicted in Supplementary Fig. 2. In this study, a GNNs model was initially developed, taking 927 pathways as its input. Following 150 epochs of training, this model demonstrated an AUC of 0.961 within the training set and an AUC of 0.841 within the testing set. Concurrently, the IG scores for these 927 pathways were calculated, with comprehensive details provided in Supplementary Table 3. Notably, 15 pathways exhibited z-scores exceeding 4, as listed in Supplementary Table 4 in order of descending IG scores: R-HSA-70,268, R-HSA-1,606,322, R-HSA-159,231, R-HSA-5,620,922, R-HSA-159,230, R-HSA-2,871,809, R-HSA-400,685, R-HSA-5,693,532, R-HSA-2,871,796, R-HSA-163,200, R-HSA-162,906, R-HSA-70,263, R-HSA-977,606, R-HSA-5,674,135 and R-HSA-909,733. Subsequently, a new GNNs model was constructed using these 15 pathways as inputs. This model, after undergoing 250 epochs of training, achieved an AUC of 0.854 on the training set, depicted in Fig. 1A, and an AUC of 0.824 on the testing set, illustrated in Fig. 1B.

Fig. 1

Performance of GNNs model for predicting responses to immunotherapy in the melanoma patient cohort from the TIGER database: (A) ROC curve illustrating the model’s accuracy on the training set, comprising 80% of the cohort; (B) ROC curve for the testing set, representing the remaining 20%. The division between training and testing datasets adheres to an 80:20 ratio within the melanoma patient cohort sourced from the TIGER database. ROC --- The Receiver Operating Characteristic

Subsequently, we extracted a total of 762 genes within key pathways. Based on the melanoma patient cohort from the TIGER database, a univariate logistic regression analysis was performed to assess the correlation between gene expression and immunotherapy response. We identified 346 significant genes associated with immunotherapy response (detailed in Supplementary Table 5). These significant genes were then subjected to subsequent clustering analysis. The optimal number of clusters (k) was determined to be 3, based on the relative changes in the area under consensus clustering cumulative distribution function, as depicted in Fig. 2B. The clustering results are illustrated in Fig. 2A, and the PCA distribution of the three distinct clusters is presented in Fig. 2C. Significant difference in immunotherapy response among these clusters were evident, as demonstrated in Fig. 2D (p < 1.0 × 10^− 3). This analysis suggests that the identified genes have a significant capacity to differentiate between patients’ responses to immunotherapy.

Fig. 2

(A) Depicts the consensus clustering results of the melanoma patient cohort from the TIGER database; (B) Shows the variation in the area under the curve for the cumulative distribution function of consensus clustering with respect to k; (C) Illustrates the spatial distribution of clusters identified within the melanoma patient cohort through a two-dimensional principal component analysis; (D) Details the immunotherapy response across different clusters; (E) Describes the association between the lambda value in LASSO regression and the binomial deviance, leading to the inclusion of 13 genes for the development of an immunotherapy response score, termed responseScore; (F) Exhibits the ROC curve for the predictive accuracy of responseScore in evaluating immunotherapy response within the melanoma patient cohort. For comparison, the ROC curves of other published models and biomarkers (including Inflammatory response score, CXCL9, SETD, DSC2, TNFSF10, and CRKL) are also presented; (G) Presents the correlation graph highlighting the relationship between responseScore and the expression levels of 12 immune checkpoint genes in the TCGA-SKCM dataset, indicating statistical significance where *p < 0.05

Following this, a LASSO regression was applied to the set of 346 genes. The relationship between the lambda value in LASSO regression and the binomial deviance is depicted in Fig. 2E. Taking into account the balance between binomial deviance and gene count, we ultimately determined to construct a gene signature comprising 13 genes, which we have named after responseScore. The formula for calculating the responseScore employs the LASSO regression coefficients, each amplified tenfold to serve as weights. The detailed calculation formula for responseScore can be found in Supplementary Table 6.

Upon examining the efficacy of responseScore in predicting immunotherapy response in the melanoma patient cohort, it was observed that the AUC stood at 0.745 (95%CI: [0.693–0.797], as illustrated in Fig. 2F). To comprehensively assess the effectiveness of the model, we also validated the performance of other published models and biomarkers (including Inflammatory response score [28], CXCL9 [29], SETD [30], DSC2 [31], TNFSF10 [32], and CRKL [33]) in the melanoma patient cohort. The results (as illustrated in Fig. 2F) demonstrated that the responseScore exhibited superior performance in predicting immunotherapy response. Moreover, within the TCGA-SKCM dataset, a positive correlation was identified between responseScore and the expression of immune checkpoint genes, including PDCD1, CD274, and CTLA4, as depicted in Fig. 2G. Therefore, responseScore has proven to possess a significant capability in differentiating responses to immunotherapy.

The relationship between responseScore and prognosis

The dataset was divided into groups with high and low responseScore by setting the median as the threshold. Within the TCGA-SKCM dataset, the distribution of clinical variables and the expression of related genes in both the high and low groups is depicted in Fig. 3A. Kaplan-Meier analysis of the TCGA-SKCM dataset revealed a marked survival rate difference between the groups (p < 1.0 × 10^− 3, as shown in Fig. 3B). This result also highlighted responseScore’s robust capability in distinguishing prognosis (Fig. 3B). As shown in Fig. 3B, the 3-year and 5-year overall survival rates for the high-score group are 0.611 and 0.519, respectively, while those for the low-score group are 0.367 and 0.293. Results from multivariate Cox regression analysis confirmed that responseScore acted as an independent prognostic factor that did not depend on clinical factors including gender, race, age, and tumor stage (Fig. 3C). Further refinement of the multivariate Cox model by excluding non-significant variables such as gender solidified the independent prognostic value of responseScore (Fig. 3D). According to the outcomes of the multivariate Cox analysis (Fig. 3D), a prognostic nomogram was developed (Fig. 3E), which predicted 1-, 3-, and 5-year survival with AUCs of 0.703, 0.693, and 0.671, respectively (Fig. 3F).

Fig. 3

Insights from the TCGA-SKCM dataset: (A) Heatmap visualization delineates gene expression alongside clinical variable distribution across groups with high and low responseScore; (B) Kaplan-Meier curves coupled with stratified analysis distinguish survival outcomes between these high and low responseScore groups; (C) A forest plot illustrates results from a multiple-factor Cox regression analysis, assessing the significance of responseScore alongside diverse clinical variables; (D) Another forest plot elaborates on the significance of responseScore when examined with selected clinical variables through multiple-factor Cox regression; (E) A constructed nomogram integrates responseScore with various clinical variables to predict survival at 1, 3, and 5 years; (F) The ROC curve of the nomogram evaluates its predictive accuracy for survival periods of 1, 3, and 5 years. ROC --- The Receiver Operating Characteristic

The immune landscape of high and low groups

By analyzing with ESTIMATE algorithm, it was observed that the Immune Score in individuals grouped by higher responseScore significantly exceeded those with lower scores (p < 0.001, as shown in Fig. 4A), indicating enhanced immune cell infiltration in the former. Concurrently, a notable increase in the Stromal Score was also documented for the high-score group (p < 0.001, represented in Fig. 4B). Evaluation of immune cell infiltration through the XCELL algorithm revealed that the abundance of several immune cell types, including memory B cells, macrophages (M0 and M1 types), monocytes, activated NK cells, and CD4 + memory as well as CD8 + T cells, was substantially greater in the high-score group (Fig. 4C). Furthermore, the analysis correlating responseScore with the level of infiltration by 21 immune cell types, such as activated CD4 + T cells, activated CD8 + T cells, Gamma delta T cells, and NK cells, demonstrated a positive association (Fig. 4D). When comparing the scores of 26 Functional Gene Expression Signatures (Fges) between the two groups, the result indicated that both pro-tumor and anti-tumor immune signatures were significantly higher in the high-score group (Fig. 4E), with the latter encompassing Checkpoint molecules, Effector cells, T cells, Th1 signature, MHC I and other anti-tumor immune markers. Some studies [34] have divided melanoma into immune phenotypes ranging from 1 to 6, with phenotype 1 and 2 linked to a poor cytotoxic immune scenario, phenotype 3 and 4 to a mildly cytotoxic immune scenario, and phenotype 5 and 6 to a highly cytotoxic immune scenario. Our analysis (Fig. 4F) revealed that phenotypes 1 and 2 were less common in the high-score group (28%) compared to the low-score group (86%), whereas the frequencies of phenotypes 3, 4, and especially 5 and 6, were significantly higher (35% and 37%, respectively) compared to the low-score group (10% and 4%). The difference was significant, as determined by chi-square test. This suggests a highly cytotoxic immune environment in the high-score group, indicative of enhanced immune activity against tumor cells.

Fig. 4

Comparative analysis within the TCGA-SKCM dataset reveals: (A) Bar plots illustrating differences in Immune Score between the high responseScore group and low responseScore group, as determined by the Estimate algorithm; (B) Differences in Stromal Score between the two groups, also evaluated through the Estimate algorithm, depicted in bar plots; (C) Bar plots showcasing variance in immune cell infiltration assessed by the XCELL algorithm between groups with high and low responseScore; (D) A graph of correlation coefficients illustrating the relationship between responseScore and levels of immune cell infiltration; (E) A heatmap displaying differential scores of 26 functional gene expression signatures (Fges) between the two groups; (F) Analysis of 6 immune subtypes, with phenotypes 1 and 2 linked to a low cytotoxic immune scenario, phenotypes 3 and 4 to a moderate cytotoxic scenario, and phenotypes 5 and 6 to a high cytotoxic scenario, comparing high and low groups; (G) Bar plots comparing microsatellite instability (MSI), distinguishing between MSS/MSI-L and MSI-H subtypes, across high and low groups; (H) Comparison of responseScore across different MSI subtypes shown in bar plots; (I) Bar plots revealing differences in Exclusion scores between the groups. Significance levels are denoted as ns p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001

Investigating the relationship between responseScore and microsatellite instability revealed a distinct distribution between high and low responseScore groups. The high-score group exhibited a composition of 66% MSS/MSI-L and 34% MSI-H, in contrast to the low-score group, which had a 74% MSS/MSI-L and 26% MSI-H composition, as illustrated in Fig. 4G. Significantly, the MSI-H subtype showed a higher responseScore compared to the MSS/MSI-L subtype (p = 0.011, presented in Fig. 4H), underscoring a significant association between responseScore and microsatellite instability. Furthermore, analysis indicated that the low responseScore group exhibited elevated levels of Exclusion score (p < 0.001, Fig. 4I). This suggests that the tumor microenvironment in the low responseScore group is more conducive to immunosuppression.

Genetic variation

In the TCGA-SKCM dataset, our analysis extended to the genetic variations between the high and low responseScore groups, encompassing both somatic mutations and the tumor mutation burden. The mutation profiles for each group are depicted through waterfall plots in Fig. 5A and B. The prevalence of somatic mutations was 92.21% in the high group and 89.79% in the low group. Notably, the top five genes mutated at the highest frequencies were identical in both groups: TTN, MUC16, BRAF, DNAH5, and PCLO. Figure 5C highlights the difference in somatic gene mutations between the groups, with mutations in MXRA5 (p < 1.0 × 10^− 3), COL5A1 (p < 1.0 × 10^− 3), ADAMTS6 (p < 1.0 × 10^− 2) and SPTA1 (p < 1.0 × 10^− 2) being significantly more prevalent in the high group, whereas mutations in SLIT3 and other genes were more common in the low group (P < 5.0 × 10^− 2). The enrichment of tumor-related signal pathways based on the mutated genes of both groups is shown in Fig. 5D and E, revealing similar pathways enriched in both groups, including RTK-RAS, WNT, NOTCH, Hippo, PI3K, Cell_Cycle, MYC, TGF-Beta, TP53, and NRF2. Moreover, Fig. 5F compares the tumor mutation burden between two groups, indicating a higher mutation burden in the high group, though without statistical significance.

Fig. 5

Genetic variations in the TCGA-SKCM dataset: (A) Oncoplots depicting gene mutations within the low-score group; (B) Oncoplots for the high-score group illustrating distinct patterns of gene mutations; (C) A forest plot highlights the variances in gene mutations between the high and low-score groups; (D) Pathway analysis revealing significantly enriched pathways of mutated genes in the high-score group; (E) Similar analysis for the low-score group, showing distinct pathways significantly enriched due to mutations; (F) A boxplot comparison of tumor mutation burden (TMB) between the high and low-score groups, indicating statistical significance with ns p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001

Functional analysis

To delve deeper into the mechanisms associated with responseScore, pathway analyses utilizing the GO and KEGG databases were conducted. The GO analysis result highlighted pathways such as organelle fission, development, embryonic organ development, nuclear division, and axonogenesis, as depicted in Fig. 6A). The main pathways involved in KEGG pathway enrichment analysis encompass the cAMP signaling pathway, Ras signaling pathway, MAPK signaling pathway, and PI3K-Akt signaling pathway, among others (as shown in Figs. 6B and 7).

Fig. 6

Functional enrichment analysis in TCGA-SKCM dataset: (A) Barplot illustrating the results of GO enrichment for groups with high and low responseScore, showcasing the differential biological processes, cellular components, and molecular functions; (B) A bubble chart represents the outcomes of KEGG enrichment analysis, contrasting the pathway involvement for high versus low responseScore groups. GO --- Gene Ontology. KEGG --- Kyoto Encyclopedia of Genes and Genomes

Relationship between each included gene and prognosis

Utilizing the melanoma patient cohort from the TIGER database, a clustering analysis was conducted on the gene expression matrix comprising 13 genes to identify distinct clusters. The analysis determined the optimal number of clusters (k) to be 3, based on the shift in the area under the curve of the consensus clustering cumulative distribution function, as shown in Fig. 7A. The distribution of clusters is depicted in Fig. 7B, and a notable difference in immunotherapy response among the three clusters was observed, with statistical significance (p < 1.0 × 10^− 2) highlighted in Fig. 7C. Similarly, using the TCGA-SKCM dataset, the expression matrix of the same 13 genes was subjected to clustering analysis, corroborating the initial findings by setting the k value to 3, as evidenced in Fig. 7D. The clustering outcome is illustrated in Fig. 7E. Noteworthy, cluster A demonstrated a significantly higher OS than cluster C, with the difference being statistically significant (p < 5.0 × 10^− 2), as showcased in Fig. 7F.

Investigating the prognostic value of genes included in the LASSO model within the TCGA-SKCM dataset, univariate Cox analysis was conducted for each gene, whose results are illustrated in Fig. 7G. This analysis identified 8 genes as significant prognostic factors. While high expression levels of PSMB6, SLC25A11, and VPS28 were linked to a poorer prognosis, the remaining genes generally showed a positive association with favorable outcomes. Using the median gene expression as a cutoff point for grouping, K-M analysis of each gene further verified the ability of each gene to perform prognostic stratification. Specifically, the survival rates for patients with high expression of PSMB6 and SLC25A11 were significantly lower compared to those with lower expression levels (p < 0.05 and < 0.001 respectively, demonstrated in Fig. 7H).

Exploring the expression levels of 13 genes in tumor versus normal tissues within the TCGA-SKCM dataset revealed significant findings, as depicted in Fig. 7I. Notably, genes including PSMB6 and SLC25A11 were expressed at significantly higher levels in tumor tissues compared to normal tissues (both p < 1.0 × 10^− 3, Fig. 7I). Furthermore, an investigation into the co-expression patterns of these genes showed that PSMB6 expression positively correlated with that of SLC25A11, VPS28, among others (p < 1.0 × 10^− 3), while exhibiting a negative correlation with CD247, IRF8, ZBP1, and other genes (p < 1.0 × 10^− 3, as shown in Fig. 7J).

Fig. 7

Clustering and gene expression analysis across datasets: (A) Displays the variation in area under the curve for consensus clustering cumulative distribution function with respect to k in the melanoma patient cohort from the TIGER database; (B) Shows consensus clustering of the expression matrix of 13 genes included in responseScore within the melanoma patient cohort; (C) Details the immunotherapy response for different clusters in the melanoma patient cohort; (D) Illustrates the change in area under the curve for consensus clustering cumulative distribution function with k in the TCGA-SKCM dataset; (E) Depicts consensus clustering of the gene expression matrix for the 13 selected genes in the TCGA-SKCM dataset; (F) Compares overall survival periods across different clusters in the TCGA-SKCM dataset; (G) A forest plot presenting the results of single-factor Cox analysis for the 13 selected genes in the TCGA dataset; (H) Kaplan-Meier analysis results for each of the 13 included genes in the TCGA-SKCM dataset; (I) A barplot contrasting the expression levels of each gene between normal and cancer tissues in the TCGA-SKCM dataset; (J) Correlation coefficient graph for the expression levels of each included gene in the TCGA-SKCM dataset. (ns p > 0.05; * p < 0.05)

Each included gene and immune-related genes

Based on the TCGA-SKCM dataset, an analysis was conducted to explore how each gene under study interacts with the immune microenvironment, particularly focusing on their relationships with genes associated with chemokines, immune checkpoints, MHC molecules, immune stimulation, and receptors. It was found that genes such as PSMB6, SLC25A11, and VPS28 generally exhibited a negative correlation with a majority of chemokine-related genes (e.g., CCL5, CCL19, CCL25, CXCL13, as indicated in Fig. 8A), immune checkpoint genes (e.g., CD244, CD96, LAG3, TIGIT, as shown in Fig. 8B), MHC-related genes (e.g., HLA-DOB, HLA-DQA1, HLA-DMB, as seen in Fig. 8C), genes related to immune stimulation (e.g., TNFRSF17, CD40LG, as depicted in Fig. 8D), and receptor-related genes (e.g., CCR4, CCR8, as illustrated in Fig. 8E). Conversely, a positive correlation was observed with a select few chemokine-related genes (e.g., CCL28, shown in Fig. 8A), MHC-related genes (e.g., TAPBP, presented in Fig. 8C), genes involved in immune stimulation (e.g., CD276, TNFSF9, PVR, depicted in Fig. D), and receptor-related genes (e.g., CXCR1, indicated in Fig. 8E) (p < 1.0 × 10^− 2). The correlation patterns between immune-related genes and CD247, IRF8, ZBP1, and others were typically inverse to those observed with PSMB6, SLC25A11, and VPS28 (as detailed in Fig. 8A-E).

Fig. 8

Gene expression correlations in the TCGA-SKCM dataset: (A) Graphs depicting correlations between the expression of genes included in responseScore and chemokine factor-related genes; (B) Correlation analyses between the expression of the 13 selected genes and immune inhibitor-related genes; (C) Correlation graphs for the expression of the 13 genes and MHC-related genes; (D) Analyses illustrating correlations between the expression of the 13 genes and immune stimulation-related genes; (E) Graphs showing correlations between the expression of each gene included in responseScore and receptor-related genes. Significance levels are indicated with *p < 0.05; **p < 0.01; ***p < 0.001

Functional verification experiment results of PSMB6

To delve deeper into the role of genes associated with responseScore and identify potential targets for melanoma treatment, biological experiments were conducted focusing on PSMB6. We detected the expression of PSMB6 in PCDH-NC group and PCDH-PSMB6 group by qPCR and Western Blot. The qPCR results in Fig. 9A showed that PSMB6 was significantly overexpressed in the CDH-PSMB6 group. The WB results in Fig. 9B also showed significant overexpression of PSMB6. We sorted CD8 + T cells from human blood and then flow cytometry identified CD8. The flow cytometry results in Fig. 9C showed that CD8 + T cells had been effectively sorted. CD8 + T cells were mixed and co-cultured with PCDH-NC group and PCDH-PSMB6 group, and then ELISA was used to detect cytokines released by CD8 + T cells in the supernatant. The ELISA results in Fig. 9D showed that cytokine IL-6 was significantly reduced in PCDH-PSMB6 groups of both melanoma cell lines (A375, A875), while TNF-α and IFN-γ were not significantly different, and IL-1β was significantly reduced in A375 cell line but not significantly different in A875 cell line.

Fig. 9

Experimental results in melanoma cell lines A375 and A875: (A) qPCR assessed PSMB6 mRNA levels in both the PCDH-PSMB6 and PCDH-NC groups, revealing distinct expression patterns; (B) PSMB6 protein abundance was detected via Western Blot in the two aforementioned groups, demonstrating differential protein expression; (C) Flow cytometry analysis identified CD8 + T cells. The results showed that CD8 + T cells were effectively sorted; (D) Enzyme-linked immunosorbent assay (ELISA) measured the secretion of cytokines IL-6, IL-1β, TNF-α, and IFN-γ in both groups across the melanoma cell lines (A375, A875), indicating cytokine release profiles. The error bars in the figure represent the standard deviation

Discussion

In this study, the first focus was on the pathways that were associated with immunotherapy response. Pathway refers to the interactions among genes and other substances, that demonstrate distinct biological functionalities. Such interactions encompass numerous genes and molecules, with their regulation occurring through diverse mechanisms. Conventionally, these interactions are represented via graphical data structures. GNNs, a class of deep learning methodologies, excel in capturing the topological characteristics of these graphical data structures. This superiority of GNNs over convolutional neural networks and other machine learning methods is particularly evident when processing graphical data, owing to their enhanced performance and interpretability. In the context of this study, GNNs were used to predict the response of melanoma to immunotherapy, specifically immune checkpoint inhibitors like CTLA-4 and PD-1 inhibitors. The results were promising, with a model incorporating 927 pathways achieving an AUC of 0.961 in the training set. To generalize the predictive ability of the model, we selected 15 key pathways based on IG score and constructed another GNNs model with 15 pathways as input. The GNNs model leveraging these 15 pathways demonstrated a robust predictive capacity, evidenced by an AUC of 0.854 on the training set and 0.824 on the testing set. The model’s utilization on gene set information and the inherent topological information hidden in pathways contributes to its superior predictive accuracy. These findings suggest that pathway-based analysis and GNNs hold significant promise for predicting responses to immunotherapy. On the other hand, we used the first GNNs model to quantify the importance of different pathways in predicting immunotherapy responses, identifying 15 key pathways. The excellent predictive performance of the second GNNs model further validated the significance of these 15 pathways, suggesting the presence of key genes related to immunotherapy responses within these pathways. This provides a foundation for subsequent bioinformatics analyses.

In the process of this research, the application of GNNs interpretive algorithms facilitated the identification of key pathways related to the response of skin melanoma to immunotherapy. Among these, the pathway named Pyruvate metabolism, labeled R-HSA-70,268, which governs pyruvate metabolism, emerges as a critical node within the network of energy metabolism pathways. Pyruvate, emerging as a pivotal molecule, serves as a juncture where it can undergo conversion into acetyl-CoA through the action of the pyruvate dehydrogenase complex, thereby feeding into the tricarboxylic acid (TCA) cycle, or alternatively, it can be transformed into lactate by the lactate dehydrogenase. Lactate plays a significant role in modulating the behavior of various immune cells, such as lymphocytes, regulatory T cells, natural killer cells, monocytes, dendritic cells, and macrophages, all of which are important in the tumor microenvironment [35]. Reports further elucidate that an increased concentration of pyruvate in the mitochondria of macrophages may initiate the activation of the PD-1/PD-L1 immune checkpoint, thereby negatively impacting the immune efficacy of T cells [36]. This evidence the profound connection of the pyruvate metabolism pathway to the efficacy of immunotherapy. Another pathway, R-HSA-1,606,322, is named ZBP1(DAI) mediated induction of type I IFNs. ZBP1, or Z-DNA binding protein-1, has been identified as a regulator of necroptosis in cancer-associated fibroblasts, which are implicated in creating an immunosuppressive environment conducive to tumor metastasis [37]. The compound CBL0137 has shown promise in inducing ZBP1-dependent necroptosis of cancer-associated fibroblasts, a finding corroborated by studies in melanoma mouse models [37]. Additional pathways of interest identified in this study include R-HSA-159,231, named Transport of Mature mRNA Derived from an Intronless Transcript; R-HSA-5,620,922, named BBSome-mediated cargo-targeting to cilium; and R-HSA-159,230, named Transport of the SLBP Dependant Mature mRNA. The identification and analysis of these pathways underscore their potential in modulating the response to immunotherapy in melanoma, thereby offering novel avenues for enhancing the efficacy of immune checkpoint blockade and other immunotherapy interventions.

In this study, a gene signature was formulated, termed responseScore, which is both linear and computationally straightforward, derived from the GNNs model. The responseScore presented a significant positive correlation with the expression levels of immune checkpoint genes including PDCD1, CD274, CTLA4, among others within the TCGA-SKCM dataset. What’s more, it exhibited an AUC of 0.745 in predicting the outcome of immunotherapy in the melanoma patient cohort. We also compared the responseScore with currently common models and biomarkers (including Inflammatory response score, CXCL9, SETD, DSC2, TNFSF10, and CRKL), and found that the performance of responseScore far exceeded the others. These all underscores its robust capability in differentiating responses to immunotherapy. This means that the responseScore can be applied to the stratification of melanoma patients, distinguishing groups with good responses to immunotherapy. It can assist clinicians in making decisions regarding immunotherapy, thereby promoting precision treatment. In terms of decision-making for melanoma immunotherapy, the responseScore is expected to surpass other gene signatures and biomarkers, becoming the primary standard of measurement.

Further analysis was conducted to understand the underlying mechanisms of responseScore, examining genetic variations, pathway enrichment, and the tumor microenvironment. Genetic variation analysis revealed that the high responseScore group was associated with enhanced tumor immunity, potentially due to immunogenic events like high tumor mutation burden (Fig. 5F), thereby indicating a more favorable response to immunotherapy. Moreover, mutations in genes such as COL5A1 [38] and SPTA1 [39] within the high responseScore group demonstrated a strong linkage to tumor immunity and the effectiveness of immunotherapy. In the results of pathway enrichment analysis, PI3K-Akt signaling pathway was proven to affect the PD-L1 expression level of CD8 + T cells, thereby affecting tumor immunity [40]. The PI3K-Akt signaling pathway has been proven to regulate the process of metabolic reprogramming centered on pyruvate metabolism in various cancers [41,42,43,44]. This is highly correlated with the R-HSA-70,268 pathway identified in GNNs pathway analysis, suggesting that the interaction between PI3K-Akt and R-HSA-70,268 pathways could be a potential entry point for future cancer treatment strategies. Explorations into the MAPK signaling pathway’s impact on tumor immunity have revealed significant findings. It has been discovered that head and neck squamous cell carcinoma patients with MAPK mutations exhibit a tumor microenvironment rich in CD8 + T cells and active tumor immunity, leading to prolonged survival time after anti-PD1/PD-L1 immunotherapy compared to patients without such mutations [45]. The high correlation between the MAPK signaling pathway and the R-HSA-70,268 pathway [46, 47] suggests its potential application in tumor immunotherapy. Research by Chenggang Gao and colleagues [48] demonstrated that the Ras signaling pathway could up-regulate PD-L1 expression, thus suppressing the cytotoxic functionality of T cells within the tumor environment, while Jesse Boumelha et al. [49] elaborated from various aspects that Ras-controlled signaling network has immunosuppressive properties. These studies collectively underline the intricate link between the Ras pathway and immunotherapy. Extensive discussions have underscored the significance of the cAMP signaling pathway in modulating tumor immunity and enhancing the efficacy of immunotherapy [50, 51].

Investigations into the infiltration of immune cells of the two groups have identified a pronounced increase in the levels of immune cells, in individuals exhibiting higher responseScore relative to their counterparts with lower scores, as depicted in Fig. 4C. Furthermore, an correlation has been found between responseScore and the abundance of 21 immune cell types, including but not limited to, activated CD4 + T cells, activated CD8 + T cells, Gamma delta T cells, NK cells, Th1 cells, Th2 cells, each playing a crucial role in fortifying anti-tumor responses [52]. Analysis depicted in Fig. 4A reveals a marked elevation in the Immune Score within the high responseScore group compared to the low group. The scores of Fges marked by anti-tumor immunity such as Antitumor cytokines, Co-activation molecules, NK cells, Checkpoint molecules, Effector cells, T cells were also significantly higher in the high-score group (Fig. 4E). Regarding microsatellite instability, the high responseScore group had a higher incidence of the MSI-H subtype than the low-score group (Fig. 4G). Furthermore, the high-score group exhibited notably lower Exclusion Score (Fig. 4I) compared to the low-score group. These all confirmed from different aspects that the anti-tumor immunity of the high responseScore group was stronger than that of the low group, which lead to the better effect of immunotherapy.

This study investigated the prognostic value of responseScore. Through Kaplan-Meier analysis and multivariate Cox regression applied to the TCGA-SKCM dataset, responseScore was validated as an independent factor for prognosis. Additionally, a prognostic nomogram was developed to forecast survival rates in melanoma patients, demonstrating its predictive accuracy with AUC values of 0.703, 0.693, and 0.671 for 1-year, 3-year, and 5-year survival, respectively (Fig. 3F). These findings highlight its effectiveness in prognosis prediction. The ability of responseScore in prognostic prediction can be applied to the risk stratification of melanoma patients. It will assist clinicians in formulating personalized treatment plans based on the specific conditions of melanoma patients, thereby improving treatment outcomes. Additionally, it will help healthcare institutions allocate medical resources more efficiently, ensuring that high-risk patients receive timely and adequate treatment.

In this investigation, 13 genes were identified as integral components of the responseScore, exhibiting a strong association with melanoma prognosis and response to immunotherapy, meriting further exploration. Clustering analysis partitioned melanoma patients into distinct subgroups based on these genes, showing varied responses to immunotherapy (refer to Fig. 8C) and differing prognosis within the TCGA-SKCM dataset (refer to Fig. 8E). Thus, this gene set is pivotal for categorizing melanoma patients to differentiate between their immunotherapy responses and prognosis. All 13 genes are differentially expressed between melanoma and normal tissues. Their prognostic value was underscored by Kaplan-Meier analysis results (refer to Fig. 8H) and univariate Cox regression results (refer to Fig. 8G) within the TCGA-SKCM dataset. Furthermore, correlation analyses with chemokine factors (refer to Fig. 9A), immune inhibitor-related genes (refer to Fig. 9B), MHC genes (Fig. 9C), immune stimulation-related genes (refer to Fig. 9D), and receptor genes (refer to Fig. 9E) elucidated their deep ties to tumor immunity. These findings collectively underscore the utility of these 13 genes as potential biomarkers or therapeutic targets.

This study focused on PSMB6, a constituent of the 20 S core particle of the proteasome. The proteasome, a sizable complex with multiple subunits, is pivotal in degrading proteins through the ubiquitin-proteasome pathway and plays a role in various biological processes including apoptosis, angiogenesis, cell adhesion, and transcription [53]. The 20 S proteasome has been suggested as a biomarker for diverse diseases [54]. Studies have shown PSMB6’s association with chronic hypoxic pulmonary hypertension progression and its involvement in the hypoxia-induced pulmonary vascular remodeling in rats [55]. Additionally, PSMB6 expression has been observed to increase in lung and mesenchymal thyroid cancers [56], and it has been identified as an independent prognostic indicator in glioma [57]. This investigation revealed PSMB6 as a significant prognostic marker for melanoma, where its elevated expression in melanoma tissues correlates with poor prognosis. Examination of PSMB6 expression in relation to immune-related genes indicated a negative correlation with the majority of chemokine genes, immune inhibitor genes, MHC genes, immune stimulation genes, and receptor genes. To delve deeper into PSMB6’s impact on tumor immunity, experimental studies were conducted on melanoma cell lines transfected to create PCDH-NC and PCDH-PSMB6 groups. The two groups were co-cultured with CD8 + T cells, followed by the measurement of cytokines released by T cells using ELISA. The results demonstrated a significant reduction in cytokines IL-6 and IL-1β in the PCDH-PSMB6 group, confirming PSMB6’s high expression as an inhibitor of tumor cell immunity. Hence, PSMB6 emerges as a potential therapeutic target for melanoma, offering promising prospects for clinical treatment.

In conclusion, this study developed a GNNs model to predict the response of skin melanoma to immunotherapy, pinpointing key pathways and establishing the gene signature named responseScore. The study explored the application value and underlying mechanisms of responseScore from multiple perspectives. Concurrently, it pinpointed significant genes that play roles in the immunotherapy response, revealing their potential as markers for tumor identification or as targets for therapeutic intervention. Additionally, the study provided a focus on PSMB6’s significance in the immunotherapy response and prognosis for skin melanoma, conducting experimental inquiries into the gene’s impact on tumor immunity.

Nonetheless, the study is not without its limitations. On one hand, the validation of the model is limited to public datasets such as melanoma patient cohorts and TCGA-SKCM dataset, lacking broader clinical validation. The application value of the model requires further confirmation through multi-center, large-sample prospective studies. On the other hand, differences in ethnicity, region, and treatment strategies of various medical institutions may pose challenges to the generalizability of the model. When applying the model in clinical settings at different medical institutions, adjustments may be necessary to suit their clinical conditions. Despite these limitations, the research presented herein offers substantial value.