RNA-seq data of primary gliomas (TCGA-LGG and TCGA-GBM) from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/) and normal cerebral tissues from the Genotype-Tissue Expression (GTEX) project (https://www.gtexportal.org/) were downloaded from the UCSC Xena (https://xena.ucsc.edu/) database as the TPM data type. The Sequencing and Array datasets (CGGA-693, CGGA-325, and CGGA-301) were downloaded from the CGGA database (http://www.cgga.org.cn/) as validated cohorts. The sequencing data from CGGA were converted into the TPM data type to correct the effect of transcription sequencing depth and gene length. All the data were screened with criteria of completed survival information and survival times >30 (overall survival <30 was considered that glioma is not the critical factor for the short survival). Finally, we used the “sva” R package to adjust the batch between different cohorts and transformed “log2 (TPM+1)” for the follow-up analysis. The clinical information of three cohorts were showed in the Table 1 . Additionally, the protein expression levels of necroptosis-related genes between normal cerebral and glioblastoma tissues were obtained from the University of Alabama at Birmingham, cancer data analysis is portal database (UALCAN, Ualcan.path.uab.edu/analysis).

Necroptosis-related genes were obtained from the GSEA website (https://www.gsea-msigdb.org/) using the keyword necroptosis. The necroptosis genes were successively filtered using univariate Cox regression analysis, least absolute shrinkage, selection operator (LASSO) Cox regression analysis, and multivariate Cox regression analysis to determine the risk genes and build the risk model. The risk signature was constructed based on the expression levels of critical necroptosis genes and their coefficients. The risk score was calculated using the following formulas:

In this formula, Exp and β represent the expression levels and coefficients of each critical gene, respectively. We used the Kaplan–Meier survival curve in the TCGA cohort to define the optimal cutoff risk score. Based on the optimal cutoff value, glioma patients from the TCGA cohort were split into low and high-risk score subgroups. Kaplan–Meier survival curves and receiver operating characteristic (ROC) curves were used to detect the efficiency of the risk signature for predicting the survival of primary glioma patients. Two CGGA cohorts were used to verify the feasibility of the necroptosis risk signature.

Univariate and multivariate analyses were performed to investigate whether the risk signature can serve as a biomarker to evaluate the prognosis of glioma patients. Subsequently, we explored the relation of the risk signature to various clinical traits using stratified analysis. We also developed a nomogram based on the risk score and identified clinical characteristics with prognostic significance based on univariate Cox analysis. The nomogram performance was evaluated by the calibration curve and ROC curves at 1, 3, and 5 years.

To illuminate the difference of enrichment in the high-risk group compared with the low-risk group, we used GSEA to determine the results of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment. Based on the enrichment results associated with the immune response, we evaluated the tumor microenvironment (TME) of glioma patients in the TCGA cohort using the “estimate” R package. Subsequently, the absolute infiltration fraction of 22 immune cells of glioma patients from the TCGA dataset was calculated using the online website CIBERSORTx (https://cibersortx.stanford.edu/). From the article “The Immune Landscape of Cancer” (19), we obtained information on the immune subtype of glioma patients in the TCGA database and investigated the relationship with the risk score. Finally, we downloaded the results of the correlation of the risk genes RIPK1, RIPK3, FAS, and FADD with immune cells in LGG and GBM from the TIMER database(https://cistrome.shinyapps.io/timer/).

The infiltration of immunosuppressive cancer-associated fibroblasts (CAFs) and myeloid-derived suppressor cells (MDSCs) in glioma patients from the TCGA cohort was calculated using the Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard.edu/) (20). We also obtained the predicted response of glioma patients for predicting immune checkpoint inhibitor (ICI) therapy.

The bioinformatic analyses were completed in the R (version 4.1.1) programming environment. The difference between the two groups was determined using the Wilcoxon test. Spearman’s test was adopted for correlation analyses. The KM analysis and Cox regression analysis were performed using the R packages “Survival” and “Survminer”.

We constructed a protein–protein interaction (PPI) network to reveal the interactions among the 8 necroptosis-related genes ( Figure 1A ). The expression levels of the 8 necroptosis genes were compared between the 589 primary glioma patients and 207 normal cerebral cortex tissues from the TCGA and GTEX datasets. As shown by the bean plot, 7 genes (RIPK1, RIPK3, FAS, FADD, FASLG, TLR3, and TNF) were upregulated in the glioma tissues with p <0.001 ( Figure 1B ). Additionally, we compared the protein expression levels of the genes between normal and glioma tissues. The results downloaded from the UALCAN database demonstrated that most necroptosis genes had higher protein expression levels in gliomas ( Figure 1C ).

To identify the critical risk genes for constructing the risk model, 7 prognostic necroptosis-related genes were screened by univariate Cox regression analysis with a criterion of p <0.05 and used as candidate genes for the next step of the analysis ( Figure 1D ). Then, the candidate genes were filtered using LASSO regression analysis, an algorithm that minimized the risk of overfitting. Finally, 4 necroptosis genes, RIPK1, RIPK3, FAS, and FADD, were determined to be critical genes for building the risk model by multivariate regression Cox analysis. All four genes were associated with a poor prognosis (hazard ratio [HR] >1, Figure 1E ). The risk scores of glioma patients were calculated based on the coefficient and expression levels of each gene and the formula was as follows: risk score = 0.431 * RIPK1 +0.227 * RIPK3 + 0.302 * FAS +0.373 *FADD. To obtain the best effective grouping, we used the “surv_cutpoint” algorithm to calculate the optimal cutoff value of glioma patients in the training set and split the glioma patients from the TCGA cohort into low and high-risk score subgroups. To inspect the consequences of grouping, the principal component analysis (PCA) was completed and the results showed that the glioma population was clearly separated into two subgroups ( Supplementary Figure 1A ). The heatmap indicated that patients with high-risk scores had higher expression levels of the risk genes ( Figure 1F ). The Kaplan–Meier survival curve and the survival status scatter plot indicated that populations with high-risk scores had a poor overall survival outcome compared with that of the low-risk group in the TCGA cohort ( Figure 1G and Supplementary Figure 1D ). Additionally, the sensitivity and specificity of the risk scores to predict the overall survival (OS) of patients were assessed by ROC curves. The areas under the curve (AUCs) of the risk score at 1, 3, and 5 years were 0.851, 0.836, and 0.788, respectively ( Figure 1H ). Based on the training cohort optimal cutoff risk score, patients from the two CGGA cohorts were divided into two subgroups to verify the accuracy of the risk signature for predicting survival. Similar to the training set, significant differences in the distribution ( Supplementary Figures 1B, C ) and survival status ( Supplementary Figures 1E, F ) of patients were noted between the low- and high-risk subgroups in the test datasets. Similar consequences were obtained while the expression heatmap ( Figures 1I, L ), KM survival curve ( Figures 1J, M ), and time-dependent ROC analyses ( Figures 1K, N ) were obtained using the two CGGA cohorts. In conclusion, the four-necroptosis gene signature can accurately predict the overall survival of glioma patients.

The association of risk and clinicopathological parameters with OS was evaluated using univariate Cox regression and multivariate Cox analysis. As shown by the forest plots of univariate analysis, similar to various clinical features including age, tumor grade, IDH mutation status, 1p19q codeletion, and MGMT methylation, the risk score could serve as an independent prognostic factor (p<0.001, HR:2.715, 1.698, 2.132, Figures 2A–C ). Subsequently, we further performed multivariate Cox analysis to counterweigh the effects of other factors and found that the risk score remained statistically significant (p<0.05, HR: 1.295, 1.311, 2.370, Figures 2D–F ). We also compared the risk score between glioma patients with different clinical characteristics in the TCGA set and found that the population classified as the high-risk group had higher risk scores ( Figure 2G ).

Given the importance of the risk score and various clinical features, a nomogram was built that combined the risk score and prognostic clinical traits screened by univariate Cox analysis to maximize the efficiency for predicting individual prognosis ( Figure 3A ). As shown by the plot ( Figure 3B ), the 1-, 3-, and 5-year calibration curves for the nomogram demonstrated that the predicted nomogram was capable of forecasting the survival time of patients from the training set accurately and obtained similar results in the two test sets. In addition, we utilized ROC to evaluate the sensitivity of the nomogram model. The AUCs of the nomogram at 1, 3, and 5 years were 0.740, 0.737, and 0.685, respectively, in the TCGA cohort and validated by the two CGGA cohorts with all AUCs values greater than 0.65 ( Figure 3C ), indicating that the nomogram could distinguish patients with good or poor prognoses.

To determine the difference in the biological processes and pathways between the subgroups categorized by the necroptosis gene model, we utilized the R package “limma” to identify the differentially expressed genes and defined the log2-fold change (logFC) as the median gene expression level of the high-risk subgroup minus the value of the opposite subgroup in the TCGA cohort. Subsequently, GO and KEGG enrichment analyses between the low- and high-risk groups were performed using the GSEA algorithm. The results showed that these differentially expressed genes were enriched in the biological process of multiple immune cell-mediated immune response, regulation of cytokine production, multiple immune cells activation, angiogenesis ( Figure 4A ), and pathways related to cancer or immune response, including cytokine−cytokine receptor interaction, apoptosis, cell cycle, PI3K−Akt signaling pathway, NOD−like receptor signaling pathway, and cellular senescence( Figure 4B ).

To explore the association of necroptosis with immune activities, we calculated the tumor microenvironment scores using the R package “ESTIMATE”. The results of correlation scatter plots revealed that the risk score was positively correlated with the immune score, stromal score, and ESTIMATE score (R=0.61, 0.63, and 0.63, respectively), but negatively correlated with the tumor purity (R= -0.63, Figure 5A ). Additionally, we analyzed the infiltration fraction of 22 immune cells using the CIBERSORTx algorithm. The violin plot indicated the difference in 22 immune cells between the low- and high-risk groups. Multiple immune cells, including Macrophages (M0, M1, M2), regulatory T cells (Tregs), neutrophils, monocytes, resting NK cells, and CD8+ T cells were highly infiltrated in the high-risk group, whereas only plasma B cells were enriched in the low-risk group ( Figure 5B ). In addition, we visualized the absolute infiltration fraction of 22 immune cells using a percentage histogram, which suggested that the glioma population with high expression of necroptosis genes exhibited a greater level of total immune cell infiltration, which is similar to the results of TME analysis. We also combined the immune subtype information of glioma patients from the TCGA cohort to further evaluate the difference in immune activity between the low- and high-risk subgroups. The lymphocyte-depleted immune subtype (C4) had a higher infiltration fraction than the immunologically quiet subtype (C5) ( Figure 5C ). Because the proportion of the inflammatory immune subtype (C3) in the TCGA glioma cohort was very low, we excluded these C3 subtype populations for further analysis. In addition, we found that the low-risk group was mainly dominated by the C5 subtype (74%), whereas the high-risk group was dominated by the C4 subtype (95%, Figure 5C ). A significant correlation was noted between the risk grouping and the immunotype as assessed using Fisher’s test (p=0.001, Figure 5D ). The Kaplan–Meier survival curve revealed that patients with “low-risk + C5” subtypes had the best prognosis, while those with “high-risk + C4” subtypes had the worst prognosis ( Figure 5E ). Additionally, we downloaded the results of the correlation analysis of the risk genes RIPK1, RIPK3, FAS, and FADD with multiple immune cells using the online website Timer database to validate the accuracy of our immune analysis. Based on this analysis, the necroptosis genes were highly positively correlated with the infiltration of lymphocytes, neutrophils, macrophages, and dendritic cells (DC) ( Supplementary Figures 2A–D ). The correlation suggested the accuracy of our immune analysis.

To investigate the infiltration of CAFs and MDSCs in glioma patients in the TCGA cohort, we used the TIDE database. Furthermore, we constructed a heatmap to summarize and reveal the distribution of various traits in the risk signature. As shown in the plot, we found that almost all patients from the high-risk subgroup had the C4 immune subtype, tumor grade 3 or 4, and high mortality, with a higher proportion of poor prognostic factors such as IDH wild type, 1p19q non-codeletion, and MGMT unmethylation. The opposite trends were noted in the populations with a low risk. ( Figure 6A ). CAFs were significantly positively related to the risk score (R=0.65, Figure 6B ). Comparing the infiltration score of CAFs between different risk groupings and immune subtypes, we found that CAFs were more enriched in the populations with high-risk scores and C4 immune subtypes ( Figure 6C ). We also visualized the correlation of MDSCs with risk scores and different MDSCs between different groupings and tumor grades. Unexpectedly, the results indicated that MDSCs were weakly correlated with risk scores ( Figure 6D ) but they had a higher infiltration fraction in the high-risk subgroup ( Figure 6E ) and tumor grade 4 ( Figure 6F ). Additionally, we assessed the difference of CAF and MDSC infiltration fractions in patients with different clinical characteristics. Unexpectedly, the infiltration score of MDSCs in the 1p19q codeletion group was higher than in the non-codeletion group and CAFs were highly enriched in the 1p19q non-codeletion population. The remaining results showed that MDSCs and CAFs were less infiltrated in patients with favorable clinical characteristics, such as young age, IDH mutant, and MGMT methylation ( Supplementary Figures 3A, B ). Since TCGA patients have not received ICI treatment, we aimed to calculate the response rate of glioma patients to immune checkpoint inhibitors by the TIDE algorithm. The predicted results suggested that individuals with higher risk scores ( Figure 6G ) or tumor grade 4 ( Figure 6H ) obtained a better response compared to other group patients. Finally, we compared the expression levels of immunosuppressive cytokine genes in the low and high-risk groups, revealing that these genes were highly enriched in the high-risk group ( Figure 6I ). Additionally, tumor-associated chemokines ( Figure 6J ) and immune checkpoint genes ( Figure 6K ) were enriched to the same levels as immunosuppressive cytokines.