Clinical features and gene expression of TCGA, ICGC, and PCAWG patients were obtained in UCSC Xena (http://xena.ucsc.edu). The validated transcriptomic data and clinical characteristics from glioma were fetched from CGGA (http://www.cgga.org.cn). The external gene expression and prognosis datasets of LUAD, UVM, and HNSC (GSE30219, GSE31210, GSE37745, GSE50081, GSE22138, GSE41613) were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/). DRGs (ACTB, TLN1, CAPZB, STN, FLNB, IQGAP1, ACTN4, MYL6, FLNA, MYH9, MYH10, PDLIM1, CD2AP, and INF2) were extracted from Gan et al.’ disulfidptosis paper (13). Different immune cell infiltration markers were obtained from the cancer immunome atlas (TCIA) (35), Genecard (https://www.genecards.org/), GEPIA (http://gepia2.cancer-pku.cn/#index), Cibersot (https://cibersortx.stanford.edu/). The prognosis of different c-MET level glioblastoma patients treated with anti-PD1 therapy was obtained from Kaplan Meier-plotteR (http://kmplot.com/analysis/index.php?p=background).

Statistical analyses were performed with R (4.2.0) and GraphPad Prism (version 8.0.1). Discontinuous data were expressed as numbers/percentages, and continuous data were expressed as mean ± standard deviation (SD). P < 0.05 was considered a statistically significant difference.

In TCGA, the 14 validated disulfidptosis-related genes (DRGs) - ACTB, TLN1, CAPZB, STN, FLNB, IQGAP1, ACTN4, MYL6, FLNA, MYH9, MYH10, PDLIM1, CD2AP, and INF2 - were generally expressed in all 33 types of cancer (Figure 1A). The correlation analysis between the DRGs indicated that MYH9 and ACTN4 were the most positively related gene pair, while MYH10 and PDL1M1 were the most negatively related (Figure 1B). And the DRGs’ expression comparison between tumor and nontumor tissues implied a significant difference in most types of them (Figure 1C). MYH10 showed the highest 2.34-fold change between glioma and normal brain tissues among all the DRGs (Figure 1D). Moreover, the univariate Cox regression of the DRGs showed that almost all 14 DRGs could predict prognosis well in patients with glioma, kidney carcinoma (KCA), kidney renal clear cell carcinoma (KIRC), etc. (Figure 1E). Interestingly, DRGs were the completely hazardous factors in glioma (Figure 1F).

Following the ssGSEA analysis of different immune cell infiltration and programmed cell death, the correlation analysis indicated a strong association between disulfidptosis and most immune cells. For the most significant glioma, the R-value between disulfidptosis and exhausted T cells (TEX_Genecard), central memory CD8 T cell, effector memory CD8 T cell, gamma delta T cell, regulatory T cell, macrophage was over 0.5 (Figure 2A). Interestingly, the correlation between disulfidptosis and other modalities of cell death like ferroptosis (R-value = 0.651), necroptosis (R-value = 0.612), pyroptosis (R-value = 0.609), immunogenic cell death (ICD) (R-value = 0.559) are also very high in glioma compared with other types of cancer (Figure 2A). The univariate Cox regression indicated that T cell exhaustion (Tex), immature B cell infiltration, etc., were the dangerous factors in glioma patients. In contrast, the activated NK cells’ infiltration was a beneficial factor for survival (Figure 2B). More importantly, a higher T cell exhaustion (TEX_GEPIA or TEX_Genecard) could predict a lousy prognosis in the glioma cohort from TCGA (Figure 2C).

The univariate Cox regression, least absolute shrinkage and selection operator (LASSO) regression, and multivariate Cox regression were used to construct a gene signature for each type of cancer. Except for thyroid cancer (THCA) and uveal melanoma (UVM), the gene signatures that could predict the prognosis for patients with all other types of cancer, respectively, were successfully made (Figure 3). For the top 6 gene signatures ranked by c-index, i.e., the gene signature in adrenocortical carcinoma (ACC), pheochromocytoma and paraganglioma (PCPG), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), prostate adenocarcinoma (PRAD), kidney chromophobe (KICH), and thymoma (THYM), the receiver operating characteristic (ROC) curves showed a very high area under the curve (AUC) for 1-year, 2-year, 3-year, 4-year, and 5-year survival (Figure 3). And in glioma that showed the most outstanding relation between disulfidptosis and immune cell infiltration (Figures 2A, B), its 7-gene signature (risk score = 1.56709174 * APOBEC3C + (-3.2556028) * GLUD1 + (-2.0800874) * KIAA1671 + 1.08729963 * KIF4A + (-7.9141641) * RPL3 + 1.83720741 * TAGLN2 + 1.89252831 * TSPAN31) (Figures 4A–C) was further validated by dividing the TCGA cohort into a training group and a testing group. And both the Kaplan–Meier (KM) analysis and ROC curve (0.5-year, 1-year, 3-year, 5-year, and 10-year) indicated significant results in the training group, testing group, and the whole group (Figures 4D, E). Then, the multivariant Cox analysis of the gene signature and the clinical characteristics implied that the gene signature was an independent hazard factor for the prognosis of patients with glioma (Figure 4F). The nomogram indicated the relation of age, gender, DRGs gene signature, and the survival probability (0.5-year, 1-year, 3-year, 5-year, 7-year, and 10-year) for glioma patients (Figure 4G). Furthermore, the model based on age, gender, and DRGs gene signature was validated in the Chinese Glioma Genome Atlas (CGGA) with AUC over 0.72 (Figure 4I). In both glioma patients from TCGA and CGGA, there was a consistency between the predictive model and survival rate in the real world (Figures 4J, K).

The unsupervised clustering analysis based on the 14 DRGs’ expression was used to categorize the TCGA cohort into disulfidptosis (DSP)1, DSP2, and DSP3 groups (Figures 5A–E). The KM analysis suggested the DSP groups had significantly different survival in the disease‐free interval (DFI), disease‐specific survival (DSS), overall survival (OS), and progression‐free interval (PFI) (Figure 5F). In line with the KM analysis of pan-cancer, the KM analysis or univariate Cox regression in individual cancer type indicated that the 3 DSP clusters could serve as a significant survival-related factor in colon adenocarcinoma (COAD), CRCA [COAD + rectum adenocarcinoma (READ)], glioblastoma multiforme (GBM), glioma, head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), lung adenocarcinoma (LUAD), lung carcinoma (LCA), stomach adenocarcinoma (STAD), uterine corpus endometrial carcinoma (UCEC), and uveal melanoma (UVM) (Figures 5G–I). Next, the top 10 mutated genes (TP53, TTN, MUC16, etc.) were listed and compared among DSP1, DSP2, and DSP3 groups (Figures 6A, C). Besides, the disulfidptosis, stromal score, immune score, tumor purity, Tex, and tumor mutation burden (TMB) were significantly different among the 3 DSP groups (Figure 6B). Since the previous 7-gene model included APOEBC3C, the TMB between APOBEC-enriched and APOBEC-unenriched groups was also compared in each DSP group (Figure 6D). Immune cell infiltration and immune molecules differed greatly among the 3 DSP groups (Figures 6E, F). Each cancer type’s total T-cell infiltration ratio was also listed to give a whole landscape (Figure 6G). In particular, the glioma, in which DRGs models showed the most significant relationship with survival and immunity, implicated a significant difference in disulfidptosis, Tex_GEPIA, Tex_genecard, CD8 (+) T cell subtypes, immune score, and tumor purity between the two DSP subgroups (Figures 6H–J).

The weighted correlation network analysis (WGCNA) was used to extract the gene module most associated with disulfidptosis, immune cell infiltration, etc. (Figures 7A–C). Next, the ten hub genes (PRSS8, CRB3, ILDR1, ELF3, TMEM184A, AP1M2, TMC4, TJP3, CLDN7, HOXB7) within this cyan module were further abstracted by the STRING database and cytoHubba (Figures 7D, E). The refined DSP models based on the ten hub genes were then constructed by employing the best method of machine learning-randomForest, in which the training and testing cohorts have the highest AUC (Figure 7F). Moreover, compared with the original DSP groups, it could better predict the prognosis in pan-cancer patients (Figure 7G). The refined DSP models could differentiate the prognosis more evidently in patients with glioma (Figure 7H). After that, the new DSP model was also validated in pan-cancer cohorts from PCAWG and ICGC, glioma from CGGA, LUAD from GEO (GSE30219, GSE31210, GSE37745, GSE50081), and UVM from GEO (GSE22138) with significant p-value (Figures 8A–F).

Since the refined DSP model performed exceptionally well in glioma among all the types of cancer, the unsupervised consensus clustering and non-negative matrix factorization (NMF) clustering were further utilized to categorize the DRGs into different groups (Figures 9A, C). Finally, the more practical and evident two-DSP-group classification by the NMF method was chosen for further construction of gene signature. Compared with a lack of significance between the survival of some subtypes by the consensus clustering (Figure 9B), the KM analysis indicated a significant difference (p < 0.0001) between DSP1 and DSP2 with Hazard Ratio (HR) equal to 5.47 (Figure 9D). Furthermore, the blue module, most correlated with DSP subtypes classification and immune cell infiltration, was extracted by WGCNA (Figures 9E–G). Ten hub genes (IL2RB, CD96, CD3D, HOXC9, HOXC5, SLAMF6, GZMH, CD3E, GZMK, and GZMA) from this module were screened by cytoHubba to construct an enhanced refined DSP clustering model by ML in glioma (Figure 9H). Surprisingly, the glioma-customized DSP model trained from TCGA could predict survival well in the glioma cohort from CGGA (Figures 9I, J). Moreover, The DSP1 has a 3-fold immune therapy response rate than the DSP2 group by oncoPredict package prediction (R.4.2.0).

The pathway enrichment of the blue gene module implied that these genes might be involved in PD1 regulation (Figures 10A–C). The c-MET inspired us to explore its function further since it was one of the top 2 genes in both the blue module and tumor driver genes (TDG) (36) (Figure 10D). High expression of c-MET was associated with poor survival among glioma patients from TCGA and CGGA (Figures 10D, E). More importantly, the survival tendency in glioblastoma patients receiving anti-PD1 therapy agreed with the previous two cohorts (Figure 10F). Interestingly, its expression differed significantly between tumor and nontumor samples in over 90% of cancer types (Figure 10G). Interestingly, most immune markers in glioma had an expression difference between the high-c-MET and low-c-MET groups (Figures 10H, I). The expression of c-MET was positively linked with PD-L1, PD2, IL-10, IRF1, JAK3, and STAT3 (Figure 10J). Furthermore, the in-vitro experiment results indicated that the knockdown of c-MET could decrease the survival (Figure 11A) and proliferation (Figure 11B) of glioblastoma cell line ln299, which could be further enhanced by the combination treatment with cabozantinib (2μM, a c-MET inhibitor) (Figures 11A, B). In line with our previous data, the decrease of c-MET could down-regulated the p-JAK3, p-STAT3, and PD-L1 (Figure 11C). Furthermore, the Jurkat T cell co-cultured with the ln299 of c-MET knockdown obtained a higher level of IL-2, IFN-γ, and PD-1 (Figure 11D).

To further verify the regulation of c-MET on PD1/PDL1, peripheral blood mononuclear cells (PBMC) were extracted from healthy females. Through the co-culture of PBMC and glioma cells, our data showed that down-regulation of c-MET in Ln299 significantly decreased the activation of STAT3 and the expression level of PDL1 in this cell (Figure 12A). In contrast, the expression level of IL2, IFN-γ, CD8 and CXCR9 were elevated in PBMC (Figure 12A). Furthermore, extracellular level of IL2, IFN-γ, and CXCL9 were also significantly increased in the culture media (Figure 12B). Next, FACS was applied to detect the c-MET-mediated CD8+ T cell immunity inhibition. In Figure 12C, we found that the proportion of CD8+ T cells was increased a little after co-culture with glioma cells while it could return to normal level (Figure 12C). However, this phenomenon was very marginal compared with the PD1 change in CD8+ T cells. The CD3+ CD8+ T cells with high PD-1 expression elevated from 8.8% to 16% after co-cultured with ln299 cells. In contrast, the knockdown of c-MET almost reversed the T-cell exhaustion completely (Figure 12D).