Pan-cancer analysis of necroptosis-related gene signature for the identification of prognosis and immune significance

Background Necroptosis is a novel programmed cell death mode independent on caspase. A number of studies have revealed that the induction of necroptosis could act as an alternative therapeutic strategy for drug-resistant tumors as well as affect tumor immune microenvironment. Methods Gene expression profiles and clinical data were downloaded from XENA-UCSC (including The Cancer Genome Atlas and Genotype-Tissue Expression), Gene Expression Omnibus, International Cancer Genome Consortium and Chinese Glioma Genome Atlas. We used non-negative matrix factorization method to conduct tumor classification. The least absolute shrinkage and selection operator regression was applied to establish risk models, whose prognostic effectiveness was examined in both training and testing sets with Kaplan–Meier analysis, time-dependent receiver operating characteristic curves as well as uni- and multi-variate survival analysis. Principal Component Analysis, t-distributed Stochastic Neighbor Embedding and Uniform Manifold Approximation and Projection were conducted to check the risk group distribution. Gene Set Enrichment Analyses, immune infiltration analysis based on CIBERSORT, EPIC, MCPcounter, ssGSEA and ESTIMATE, gene mutation and drug sensitivity between the risk groups were also taken into consideration. Results There were eight types of cancer with at least ten differentially expressed necroptosis-related genes which could influence patients’ prognosis, namely, adrenocortical carcinoma (ACC), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), pancreatic adenocarcinoma (PAAD), liver hepatocellular carcinoma (LIHC), skin cutaneous melanoma (SKCM) and thymoma (THYM). Patients could be divided into different clusters with distinct overall survival in all cancers above except for LIHC. The risk models could efficiently predict prognosis of ACC, LAML, LGG, LIHC, SKCM and THYM patients. LGG patients from high-risk group had a higher infiltration level of M2 macrophages and cancer-associated fibroblasts. There were more CD8+ T cells, Th1 cells and M1 macrophages in low-risk SKCM patients’ tumor microenvironment. Gene mutation status and drug sensitivity are also different between low- and high-risk groups in the six cancers. Conclusions Necroptosis-related genes can predict clinical outcomes of ACC, LAML, LGG, LIHC, SKCM and THYM patients and help to distinguish immune infiltration status for LGG and SKCM. Supplementary Information The online version contains supplementary material available at 10.1007/s12672-022-00477-2.


Introduction
Although recent authoritative statistics showed that the death rate of cancer declined over the past 30 years, cancer remains one of the primary causes of death worldwide no matter in developed or developing countries, which greatly increases economic burden and seriously affects life quality [1]. The occurrence and development of tumor involves a series of extremely complex biological processes, and the treatment effect of many tumors is still not satisfactory even under the combination of multiple therapies. It is urgent and of great importance to find novel insights and effective agents for cancer. The resistance to cell death has been identified as one of the most important characters of malignant tumors [2]. Classical theory divided cell death forms into apoptosis and necrosis, according to the whether it's under the programmed regulation of genetic materials [3]. However, in the 1990s, a new pattern of necrosis-like cell death featured by non-caspase dependency gradually emerged. Researchers found that, under the inhibition of key proteins in apoptosis pathway [such as Caspase-8 or Fas-associated protein with death domain (FADD)] and the stimulation of tumor necrosis factor α (TNF-α), the cell morphology was similar to the necrotic cell [4,5]. Then, at the beginning of the twenty-first century, the concept and process of programmed necrosis or necroptosis was gradually proposed and elaborated [6][7][8]. In 2018, the Nomenclature Committee on Cell Death (NCCD) officially defined this special form of cell death as necroptosis [9]. Unlike apoptosis which involves kinds of morphological changes, such as cell shrinkage and detachment from the surrounding cells, nucleoplasm concentration, fragmentation of nuclear membrane and nucleolus as well as the appearance of apoptotic bodies, several special biological events occur in cells undergoing necroptosis, including the damage of membranes, disorder of metabolism and the extravasation of inflammatory substances [8]. Necroptosis and apoptosis share the same initiating stage. When tumor necrosis factor receptor 1 (TNFR1) on the cell membrane surface is activated by TNF-α, TNF receptor 1-associated death domain protein (TRADD) and tumor necrosis factor and receptor related factor 2 (TRAF2) will be recruited by its death domain at C-terminal. Subsequently, TRADD and TRAF2 separately recognizes and binds to receptor-interacting protein kinase 1 (RIPK1) and cellular inhibitors of apoptosis 1 and 2 (CIAP1/2), and protein complex scaffold is formed by linear ubiquitin Chain assembly complex (LUBAC). Then, with the combination of these molecules and transforming growth factor-beta (TGF-β) activated kinase 1/TGF-β activated kinase 1 binding protein (TAK1/TAB) complex, the supramolecular structure (TNFR1 Complex I) come into being [10]. The deubiquitination of RIPK1 by the cylindromatosis (CYLD) can result in the cleavage of Complex I and the dissociation of RIPK1 as well as TRADD, where different endings of the cell happen. Complex IIa constituted of TRADD, FADD as well as Caspase-8 and Complex IIb composed of PIPK1, receptor-interacting protein kinase 3 (RIPK3), FADD and Caspase-8 would lead cell to apoptosis. The catalytic activity inhibition of caspase-8 would allow RIPK1 to phosphorylate RIPK3, which recruits mixed lineage kinase domain-like (MLKL) to form necroptosome [11,12]. MLKL migrates to cell membrane to result in necroptosis.
Necroptosis played an indispensable role in the maintenance of internal environment homeostasis and the progression of several inflammation-related diseases, such as neurodegenerative disease, ischemia-reperfusion injury and pathogen infection [10,13]. A number of studies have also revealed the significance of necroptosis induction at cancer treatment in recent years, which especially worked for the apoptosis-resistant tumors [14]. Meanwhile, with the rise of immunotherapy, the relationship between different forms of cell death and tumor immunity has gradually attracted extensive attention [15]. There was no effective anti-tumor immune response observed in the tumor area where apoptosis or necrosis occurred. However, increasing number of studies have revealed the influence of necroptosis on tumor immune microenvironment, where the results were opposite in different tumor models. Damage associated molecular patterns (DAMPs) and various cytokines and chemokines which leaked out of necroptotic cells of colon carcinoma and melanoma could strengthen cytotoxic function of CD8+ T cells and the activity of antigen-presenting cells [16][17][18]. However, the necroptotic cells of pancreatic ductal adenocarcinoma enhanced the immunosuppressive function of tumor-associated macrophage (TAM) by C-X-C motif chemokine ligand 1 (CXCL1) and Mincle signaling [19]. The studies also showed that the synergistic effect of necroptosis-promoting agents and immune checkpoint inhibitors (ICIs) could trigger long-term tumor-suppression effect in mouse models [17,18], indicating that the necroptosis induction of tumor cell was probably an effective complement to immunotherapy.
In this study, we comprehensively analyzed the necroptosis-related genes in different kinds of cancers based on data from The Cancer Genome Atlas (TCGA), Genotype-Tissue Expression (GTEx), Gene Express Omnibus (GEO), International Cancer Genome Consortium (ICGC) and Chinese Glioma Genome Atlas (CGGA). We developed novel tumor classification and constructed risk models based on necroptosis-related genes to predict patients' clinical outcomes. Immune infiltration, gene mutation and drug sensitivity were also taken into consideration.

Gene expression and clinical data collection
We obtained gene profiles, clinical features and survival information of 33 TCGA cancers from XENA-UCSC (https:// xena. ucsc. edu/). For thirteen types of cancer with no or very limited number of corresponding normal tissue samples (< 10), we obtained gene expression data of normal samples from GTEx at XENA-UCSC, namely, adrenocortical carcinoma (ACC), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), glioblastoma multiforme (GBM), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), rectum adenocarcinoma (READ), skin cutaneous melanoma (SKCM), testicular germ cell tumors (TGCT), thymoma (THYM) and uterine carcinosarcoma (UCS). Because of no relevant samples for pheochromocytoma and paraganglioma (PCPG) and sarcoma (SARC) found in GTEx, we only used TCGA data for the analysis. Mesothelioma (MESO) and uveal melanoma (UVM) were excluded from this study, for there were no normal samples in neither TCGA nor GTEx. Necroptosis-related gene list (hsa04217) was found in Kyoto Encyclopedia of Genes and Genomes (KEGG). The details of necroptosis-related genes were shown in Supplementary file 1.

Identification of differentially expressed necroptosis-related genes (DENGs), survival analysis and tumor classification
To identify DENGs between tumors and the corresponding normal tissues, the "limma" R package was applied, with |log2 (fold change)| > 1 and false discovery rate (FDR) < 0.05 as the thresholds. Then, we conducted survival analysis of DENGs in each particular type of cancer. The cancer types with at least 10 DENGs that significantly influence patients' overall survival (OS) were selected. Next, we constructed chord diagrams of the prognostic DENGs in the chosen cancers by using "circlize" and "corrplot" R packages, where Pearson correlation analysis was performed. The correlation at protein level was visualized by STRING (Version: 11.5, https:// cn. string-db. org/) through "Multiple protein" module with the "Homo sapiens" and "low confidence (0.150)" as the main parameters. Finally, based on prognostic DENGs, we used the non-negative matrix factorization (NMF) to conducted cancer classification. "NMF" R package was used, with "brunet", "10 iterations" and "clusters k ranks from 2 to 10" as the main parameters. Kaplan-Meier analysis was performed between patients' survival and the different clusters, where four survival endpoints were taken into consideration, namely, OS, disease specific survival (DSS), progression free survival (PFS) and disease free survival (DFS).

Construction and validation of DENGs-based risk model
First, batch corrections were performed between TCGA cohorts and the corresponding additional cohorts of the selected cancers by "sva" R package. Then TCGA and additional cohorts were appointed as the training sets and testing sets separately. For each cancer the training set was used to establish necroptosis-related risk model by the least absolute shrinkage and selection operator (LASSO) regression, employing "glmnet" R package, with fivefold cross-validation applied to optimize the model. Patients were classified into low-and high-risk groups according to the median risk score of training set. Kaplan-Meier analysis of OS and the risk groups was conducted. To assess the predictive efficiency of the risk model, time-dependent receiver operating characteristic (ROC) curves of 1, 3, 5-years were made using "survivalROC" R package.
Uni-and multi-variate survival analyses were employed to examine whether the risk score could independently affect patients' prognosis. Model genes expression heat maps were constructed with the increase of risk score by "pheatmap" R package, and some clinical factors between patients from low-and high-risk groups were also compared by the use of Fisher's exact test. Principal Component Analysis (PCA), t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) were carried out to verify the risk-group assignments according to the model genes expression data, where "stats", "Rtsne" and "umap" R packages were used. Distribution of patients' risk score and survival state was also analyzed. The same procedures were performed in the testing sets.

Investigation of tumor immune microenvironment
Five algorithms were applied to assess immune infiltration status of each patient in both training and testing sets, namely, CIBERSORT, EPIC, MCPcounter, ssGSEA and ESTIMATE. Then, the immune infiltration level was compared between patients from low-and high-risk groups with Wilcoxon signed-rank test. The Spearman's correlation analysis of risk score and immune score, stromal score as well as ESTIMATE score was also conducted. Then, we compared tumor mutational burden (TMB) and microsatellite instability (MSI) between the patients from the two risk groups with Wilcoxon signed-rank test, and investigated the relationship of risk score and TMB as well as MSI using Spearman's correlation analysis. In addition, we explored whether there existed a correlation of risk score and immune related genes expression with Pearson correlation analysis, including immunoinhibitor genes, immunostimulator genes, Major Histocompatibility Complex (MHC) genes, chemokine genes and chemokine receptor genes. The corresponding genes were acquired from TISIDB (http:// cis. hku. hk/ TISIDB/ index. php).

Analysis of gene mutation
Somatic mutation data based on "VarScan2" software was acquired for TCGA samples. Then, we made oncoplots to show the mutation status of the top 20 most frequently mutated genes in low-and high-risk groups, with "maftools" R package. The mutation rate of the top 20 genes was compared by Fisher's exact test.

Drug sensitivity analysis
We downloaded the gene expression and z-score matrix from CellMiner (https:// disco ver. nci. nih. gov/ cellm iner/ home. do) and calculated the risk score of each sample according to the genes and corresponding coefficient of the different cancers' risk model. Then, we investigated whether there existed any correlation of risk score and the sensitivity of Food and Drug Administration (FDA)-approved drugs with Pearson correlation analysis.

Identification of prognostic DENGs in TCGA-cancers
As shown in Fig. 1, there were eight types of cancer with at least ten prognostic DENGs, namely, ACC, CESC, LAML, LGG, LIHC, PAAD, SKCM and THYM. The situation of other cancers was shown in Fig. S1, and no prognostic DENGs was found in colon adenocarcinoma (COAD) (d), stomach adenocarcinoma (STAD) (t), thyroid carcinoma (THCA) (u) and uterine corpus endometrial carcinoma (UCEC) (v). Notably, there were no DENGs observed in SARC. We also revealed the correlation between the prognostic DENGs in the eight cancers at both transcription and translation level (Fig. 2).

Tumor classification
We used NMF to classify cancer patients into different subgroups according to the expression profiles of the prognostic DENGs. NMF rank survey with multiple parameters and the consensus matrix heat maps were displayed at K value from 2 to 10 for ACC, CESC, LAML, LGG, LIHC, PAAD, SKCM and THYM (Fig. S2). The optimal K value was chosen for each cancer and the corresponding classification was shown (Fig. 3a, c, e, g, i, k, m, o). Notably, there existed significant difference of OS among the subgroups in all cancers except for LIHC (Fig. 3b, d, f, h, j, l, n, p).

LASSO regression risk models
The LASSO coefficient spectrum of the selected necroptosis-related genes for ACC, CESC, LAML, LGG, LIHC, PAAD, SKCM and THYM were shown in Figs. 4a, g, m, s and 5a, g, m, s. Figures 4b, h, n, t and 5b, h, n, t showed the fivefold cross-validation. The risk score calculation formulas of the eight cancers were shown in Supplementary file 2. In ACC, LAML, LGG, LIHC and SKCM, low-risk patients had obviously better OS compared with patients from high-risk group (Figs. 4c, o, u, 5c, o), and the time-dependent ROC curves of 1, 3 and 5 years in training and testing sets revealed the good efficiency of our risk models at predicting patients' prognosis (Figs. 4d, p, v, 5d, p). The risk score could independently influence patients' prognosis in both training and testing sets (Figs. 4f, r, x, 5f, r). However, In CESC and PAAD, we failed to observe the statistically significant difference of patients' OS between low-and high-risk groups in the testing sets (Figs. 4i, 5i). We didn't find a THYM cohort with sufficient prognostic information, so the analyses were only conducted in TCGA cohort (Fig. 5s-x). For ACC, LAML, LGG, LIHC, SKCM and THYM, the variation trend of model genes expression with the increase of risk score was shown, along with the comparison of some clinical factors between low-and high-risk groups (Fig. 6a, d, g, j, m, p). Dimensionality reduction analysis showed that the risk groups were largely in accordance with the two dimensional pattern of PCA, t-SNE and UMAP distribution, while in the testing set of LGG (CGGA cohort), the results were less satisfactory (Fig. 6b, e, h, k, n, q). With the increase of risk score, patients' survival period was shortened and the number of deaths increased (Fig. 6c, f, i, l, o, r).

GSEA result
Gene Ontology (GO) and KEGG pathways related to the cell cycle were enriched in the high-risk group of ACC (Fig. 7a, c) and LIHC (Fig. 7e, g) no matter at training or testing sets, such as cell cycle checkpoint, cell cycle G1-S phase transition, cell cycle G2-M phase transition, chromosome segregation, DNA dependent DNA replication and splicesome, with similar situation observed in low-risk group of THYM (Fig. 7j). In addition, innate and adaptive immune-related pathways were enriched in LGG high-risk group (Fig. 8e, g) and SKCM low-risk group (Fig. 8j, l) no matter at training or testing sets, including activation of immune response, adaptive immune response, antigen presenting and presentation as well as complement and coagulation cascades. Surprisingly, in the analysis of LAML, we found visible enrichment discrepancies in high-risk group at training and testing sets, with immune-related or cell-circle-related pathways separately enriched in the two sets (Fig. 8a, c).

Immune infiltration analysis of LGG and SKCM
Based on the GSEA results above, we further explored whether there existed any immune infiltration difference between low-and high-risk groups in LGG and SKCM. According to five immune infiltration assessment algorithms, high-risk LGG patients and low-risk SKCM patients had higher level of immune infiltration and function at both training and testing sets, which accorded with the GSEA enrichment results. For LGG patients, the infiltration level of B cells, plasma cells, CD8+ T cells, macrophages, endothelial cells, cancer-associated fibroblasts (CAFs) and dendritic cells was higher in high-risk group (Fig. 9a-d), while the situation of NK cells (Fig. 9a-d) and regulatory T (Treg) cells (Fig. 9a, d) (Fig. 9f-i). As shown in Fig. 9e, immune score, stromal score and ESTIMATE score were higher in LGG patients from high-risk group at both training and testing sets, which also positively correlated with the patients' risk score. For SKCM patients, the results were opposite (Fig. 9j). Then, we took TMB and MSI into consideration and found that high-risk LGG patients possessed higher TMB level (Fig. 10a), and TMB increased with risk score (Fig. 10b). Next, we explored the relationship of risk score and the gene expression of immunoinhibitors, immunostimulators, MHCs, chemokines and chemokine receptors. As shown in Fig. 10i-m, the expression of most immune-related genes positively correlated with risk score of LGG patients in both training and testing sets, while the results were opposite for SKCM patients (Fig. 10n-r).

Gene mutation status
We explored gene mutation status between low-and high-risk groups in TCGA cohorts of ACC, LAML, LGG, LIHC, SKCM and THYM, and screened out the top 20 genes with the highest mutation frequency. Higher mutation rate of tumor protein p53 (TP53) occurred in ACC and LIHC patients from high-risk group (Fig. 11a, d). For LAML and SKCM  genes (a, g, m, s) and the fivefold cross-validation (b, h, n, t) for variable selection were shown. Kaplan-Meier plots (c, i, o, u) showed the OS difference between patients from low-and high-risk groups sorted by median risk score of the training set, with logrank p value marked in the graphs. Time-dependent receiver operating characteristic (ROC) curves of 1, 3, 5-years (d, j, p, v) showed the predictive efficiency of the risk model, with area under curve (AUC) values noted in the graphs. The forest plots showed the results of univariate (e, k, q, w) and multivariate (f, l, r, x) survival analyses  genes (a, g, m, s) and the fivefold cross-validation (b, h, n, t) for variable selection were shown. Kaplan-Meier plots (c, i, o, u), time-dependent ROC curves of 1, 3, 5-years (d, j, p, v) and forest plots (e, f, k, l, q, r, w, x) showed the prognostic effectiveness of the risk models  maps (a, d, g, j, m, p) showed the variation trend of model genes expression with the increase of risk score and the comparison of several clinical factors between lowand high-risk groups in the six selected cancers. Fisher's exact test was used. *p < 0.05; **p < 0.01; ***p < 0.001. Principal Component Analysis (PCA), t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) (b, e, h, k, n, q) confirmed the stratification of patients into low-and high-risk clusters. The scatter diagrams (c, f, i, l, o, r) showed the condition of patients' risk score and distribution of their survival time and state, with dotted line separating patients into low-and high-risk groups

Discussion
Necroptosis is a novel programmed cell death mode independent on caspase, with increasing evidence of anti-tumor effects discovered in recent years. As we know, traditional chemotherapeutic agents usually induced cell apoptosis to exert anti-tumor effects [26]. However, tumor cells are inherently anti-apoptotic. In spite of the prevalence of heterogeneity in various tumors, there's a high possibility that the subpopulation of tumor cells with greater antiapoptotic selection superiority will gradually clone and govern the entire tumor as the treatment proceeds. Therefore, drug resistance has become a common fact during clinical practice, and tumors which relapse or progress after treatment are extremely difficult to deal with [26]. Thus, it became a natural idea to induce other types of cell death for drug-resistant tumors, and alternative choices mainly included ferroptosis, pyroptosis as well as necroptosis [27]. Numerous studies have proven that the transition of apoptosis to necroptosis or the direct induction of necroptosis could make for overcoming drug resistance and inhibiting tumor development for various cancers, such as acute myeloid leukemia [28,29], breast cancer [30], osteosarcoma [31], nasopharyngeal carcinoma [32], prostate cancer [33,34] and colon cancer [35,36].
In this study, based on TCGA and GTEx data, we identified eight types of cancer with the highest number of prognostic DENGs and for the first time sorted ACC, CESC, LAML, LGG, LIHC, PAAD, SKCM and THYM patients into different subgroups based on necroptosis-related genes. Kaplan-Meier analysis of four follow-up endpoints showed that the classification was excellent in distinguishing patients' OS in all cancers above except for LIHC. Then, the risk models were set up. Unfortunately, the risk models didn't work at testing sets of CESC and PAAD, but we do find a method to efficiently distinguish patients' OS in ACC, LAML, LGG, LIHC and SKCM. The testing set of LAML (GSE37642) lacked M3-subtype patients and the testing set of LGG (CGGA) only consisted of Asian patients, so there existed some intrinsic discrepancies between TCGA cohorts (used as training set) and these testing sets. This might cause the inconsistency of AUC values between training and testing set. Notably, among these five cancer types, ACC is relatively less studied. As a rare malignancy with great complexity, the 5-year DFS rate of ACC was only about 30%, and there still existed many therapeutic challenges [37,38]. Due to the heterogeneity of ACC, the prognostic efficiency of the most widely accepted TNM staging was inevitably limited [39]. Thus, it is necessary to seek new risk factors for ACC patients. Our ACC risk model based on necroptosis-related genes has good predictive ability for patient' survival, which might provide meaningful references for patients' prognosis in the future clinical practice.
Although kinds of immunotherapies have achieved remarkable success in cancer treatment, only limited number of patients could exhibit long-lasting anti-tumor response, where tumor immune infiltration status played a significant role [40]. Identification of cancer patients with abundant infiltration of immune cells is of great importance to screen out potential candidates for immunotherapy. Our GSEA results of SKCM and LGG cohorts highlighted immune-related GO and KEGG pathways in low-and high-risk groups, which along with results of the estimated immune infiltration level based on five algorithms could contribute to the distinction of "cold" and "hot" tumors.
As we know, immunotherapies have not acquired satisfactory results in glioma patients in recent years, including adoptive lymphocyte transfer, tumor associated vaccine, viral-based therapy and ICIs, where T-cell exhaustion played a dominant role, and tumor heterogeneity, blood brain barrier as well as lack of immune organs in central nerve systems also shared the blame [41]. Although there is a higher CD8+ T cells infiltration level in high-risk LGG patients, we failed to observe the difference of cytotoxic lymphocytes between the two risk groups according to MCPcounter. Noteworthy is the infiltration level of M2 macrophages and CAFs is higher in high-risk LGG patients. Recent studies have revealed the fact that M2 macrophages played a vital part in the development of glioma by promoting tumor invasion and metastasis, facilitating tumor stemness as well as suppressing immunity of the tumor area and the whole body [42,43]. CAFs were involved in tumor cell replication, angiogenesis, chemotherapy insensitivity and the suppression of CD8+ T cell function [44,45]. M2 macrophages and CAFs have been considered as promising therapeutic targets by number of studies [44][45][46], and high-risk LGG patients perhaps benefit from the agents which inhibit M2 macrophages or CAFs.
Unlike the situation in LGG, the infiltration level of immune cells widely known for suppressing tumor development is higher in low-risk SKCM patients, including CD8+ T cells, Th1 cells and M1 macrophages. According to the correlation analysis of risk score and immune-related gene expression, SKCM patients from low-risk group also possessed a higher gene expression level of plenty of immunosuppressive molecules, some of which were identified as immune checkpoints and their therapeutic potential has been proven by numerous studies. ICIs were initially studied and applied for the clinical application in melanoma, and Ipilimumab, targeting cytotoxic T-lymphocyte-associated protein 4 (CTLA4), is the first drug in history to significantly prolong the survival period of patients with this highly malignant tumor [47]. Programmed cell death protein 1 (PD-1) antibody was also approved for the treatment of advanced melanoma by FDA in the year of 2014 and phase 3 clinical trial of Relatlimab, targeting lymphocyte-activation gene 3 (LAG-3), has met its primary endpoint of PFS, which may offer new hope for SKCM patients in the future. It needs to be mentioned that there existed a higher mutation rate of MUC16 in low-risk SKCM patients. MUC16, also known as carbohydrate antigen 125 (CA125), ranks third in the list of gene mutation frequency of cancers, whose mutation occurs most frequently in SKCM [48]. The study also showed that MUC16-mutated melanoma patients treated with ICIs had significantly longer OS. Given that our study could help to recognize SKCM patients with higher level of immune infiltration and immune-checkpoint genes expression as well as higher MUC16 mutation rate, it is reasonable to believe that low-risk SKCM patients are more likely to benefit from ICIs treatment.
Although we failed to find a cohort to check the predictive ability of prognosis in the THYM risk model, there were still some results which could arouse our attention. First, the nine-genes risk model successfully assigned all death cases into high-risk group, and the following time-dependent ROC analysis exhibited an excellent predictive ability of the model with 1, 3, 5-year OS area under the ROC curve up to 0.854, 0.936 and 0.966. Regardless of the application of which dimensionality reduction method, the cases could be obviously divided into low-and high-risk clusters. Thymoma has a low incidence and favorable prognosis, so the associated studies are relatively limited compared with other common or highly malignant tumors. For patients classified as high-risk, their review period perhaps needs to be shortened so that the tumor progression can be detected and treated in time.
For THYM, it is still controversial whether adjuvant radiotherapy or chemotherapy should be applied after surgery. According to our result, some of the patients classified as high-risk might be the potential candidates for postoperative adjuvant therapies. We noticed a decline in the sensitivity of tumor cells to vincaleukoblastinum drugs with the increase of risk score based on THYM risk model. However, irofulven exhibited anti-tumor activity in cells with high risk score, which is a kind of cytotoxic drug proven to be an effective agent for tumors with DNA repair deficiency by several studies [49,50]. This finding may provide some useful information for the clinical chemotherapy of THYM. In addition, we noticed the mutation rate of GTF2I in the high-risk patients was about twice as high as that in low-risk patients. Researchers have found that there existed a high mutation rate of GTF2I in indolent thymomas, which was extremely rare in aggressive thymomas and thymic carcinomas [51]. Mutant GTF2I, identified as a novel tumorigenic driver, can promote growth, proliferation and transformation of epithelial cell as well as alter glucose and lipid metabolism [51,52], and whether it could work as a therapeutic target requires further research.

Conclusions
In summary, this is the first study to comprehensively investigate the genes of necroptosis pathway in all TCGA cancers. We conducted NMF to classify ACC, CESC, LAML, LGG, PAAD, SKCM and THYM patients into subgroups with different prognosis. The risk model based on necroptosis-related genes can effectively predict the prognosis of ACC, LAML, LGG, LIHC, SKCM and THYM patients. The risk score contributes to the identification of immune infiltration level for LGG and SKCM patients, which could help to screen out the potential candidates who might benefit from immunotherapy. Genetic mutation status and drug sensitivity were also different for patients from different risk groups, which may offer meaningful information for the future clinical practice. Fig. 10 Tumor mutational burden (TMB), microsatellite instability (MSI) and immune-related genes expression analysis. Bar graphs showed the comparison of TMB (a, e) and MSI (c, g) between low-and high-risk groups and scatter diagrams showed the correlation between TMB (b, f) or MSI (d, h) and the risk score of LGG and SKCM patients. Wilcoxon signed-rank test p value and Spearman's correlation coefficient R value as well p value were marked in the graphs. The correlations between risk score and the expression of immunoinhibitor genes (i, n), immunostimulator genes (j, o), MHC genes (k, p), chemokine genes (l, q) as well as chemokine receptor genes (m, r) were shown, with "*" representing Pearson correlation p value < 0.05