Pan-cancer analysis of a novel indicator of necroptosis with its application in human cancer

Data sources

We downloaded the transcriptional data of 10496 samples from the UCSC Xena database (https://xenabrowser.net/datapages/) and normal expression data from GTEX project database. The list of cancer types is presented in Supplementary Table 1. The clinical parameters and follow-up information was also obtained. The gene-level copy number, DNA methylation, and somatic mutation data were also downloaded. The tumor mutation burden was calculated according to the single-nucleotide variation profiles of the cancers. Microsatellite instability data were obtained from previously published data. Necroptosis regulator data were obtained from previous studies and the Molecular Signatures Database. We ultimately included 67 necroptosis regulators from previous study [11]. The list of genes is provided in Additional file 1: Supplementary Table 2. Figure 1 presents the overall overview of our study.

Figure 1. The overall overview of the study.

Establishment of the necroptosis score (NPS) in cancers

Using the necroptosis regulator genes (NRGs), we first constructed a protein–protein interaction (PPI) network using the STRING database (STRING: functional protein association networks (string-db.org)), and we identified the hub genes by inputting the PPI data into Cytoscape software. Next, we explored the correlations among the necroptosis regulators in all cancers. Analyses of differentially expressed NRGs were also performed to compare the expression differences in tumor and normal samples (GTEX data). Then, we explored the function and pathway enrichment of these NRGs by performing GO and KEGG analyses. Subsequently, we evaluated the single-nucleotide variant (SNV), copy number variation (CNV), and methylation status in cancers. We also identified the prognostic roles played by these NRGs in cancers by log-rank test. Finally, we performed single-sample gene set variation analysis (ssGSEA) using the NRG set, and we assigned a necroptosis score (NPS) to each sample according to the enrichment levels. We completed the calculation using the R “GSEABase” packages. In ssGSEA, gene expression values for a given sample are rank-normalized, and an enrichment score is produced using the Empirical Cumulative Distribution Functions (ECDF) of the genes in the signature and the remaining genes. For a given signature G of size NG and single sample S, of the data set of N genes, the genes are replaced by their ranks according their absolute expression from high to low: An enrichment score ES (G, S) is obtained by a sum (integration) of the difference between a weighted ECDF of the genes in the signature and the ECDF of the remaining genes PNG [12]:

$$ES(G,S)={\displaystyle \sum _{i=1}^N[P_G^w(G,S, i)-P_{N_G}(G,S,i)]}$$

Expression profile based on the NPS and its correlation with clinical prognosis in cancers

After obtaining the NPS of each sample, we first evaluated the expression level in 33 cancers and the differential NPSs between tumor samples and normal samples in different cancers. Wilcoxon tests were performed to compare tumor and normal tissues in each cancer. The tumor samples were categorized into a high-NPS group and a low-NPS group on the basis of the optimal cutoff NPS value for each cancer. We evaluated the associations between NPS and prognosis in each cancer. The four types of prognosis were overall survival (OS), disease-specific survival (DSS), disease-progression interval (DFI), and progression-free interval (PFI). A Kaplan–Meier curve was assessed to compare survival curves, and univariate Cox regression was performed to calculate the hazard ratio (HR) and the corresponding confidence interval (CI). P<0.05 was considered to be significant.

Associations between NPS and related pathways and molecular features in cancers

We also performed gene set variation analysis (GSVA) using R packages (GSVA) and the latest HALLMARK pathways dataset from the MsigDB, including fifty common signaling pathways. As necroptosis is a type of programmed cell death, we explored the correlation between the NPS and certain cell death-related gene sets, including genes involved in pyroptosis, autophagy, and ferroptosis. The results were visualized in a heatmap.

The NPS and immune status in cancers

To explore the association between the NPS and immune status in each cancer, we followed three different methods. First, we used the ESTIMATE algorithm, which depicts a tumor microenvironment based on four parameters: immune, stromal, ESTIMATE scores, and tumor purity. The second method involved evaluating the correlations between NPS and the TME signature [13]. The third method was based on the CIBERSORT algorithm, which allows for the sensitive and specific discrimination of 22 human immune cell phenotypes. Finally, we explored the associations between the NPS and immune regulation genes, including MHC genes, immune-activated and inhibited genes, immune checkpoint genes, chemokines, and chemokine receptor genes, using Pearson correlation.

Effect of the NPS on immunotherapy

We performed a ssGSEA for the following datasets: GSE32894 (urothelial carcinoma), GSE135222 (advanced non-small-cell lung carcinoma patients), GSE176307 (metastatic urothelial cancer), and the Supplementary Materials of previous studies for KIRC (PMID: 32472114) [14]. Patients were categorized into a high-NPS group and a low-NPS group on the basis of the optimal cutoff value (minimum P value). Then, we performed Kaplan–Meier analysis to evaluate the effect of the NPS on prognosis in patients who received immunotherapy.

The NPS and chemotherapy sensitivity in cancers

Using the GDSC and CTRP databases, we performed Pearson correlation analysis to obtain the correlations between mRNA expression of NPGs and drug IC50. The P value was adjusted on the basis of the false discovery rate (FDR). The correlations were visualized with a bubble plot. Negative correlations indicated chemotherapy sensitivity, and positive correlations indicated chemotherapy resistance.

Data availability

All data can be download from UCSC Xena database (https://portal.gdc.cancer.gov/).

Landscape of necroptosis regulators across cancers

We ultimately obtained 67 necroptosis regulator genes. Figure 2A presents the PPI associated with these NRGs. We identified the 10 most interactive hub genes using maximal clique centrality methods: TNFRSF1A, RAF2, RIPK3, CFLAR, DIABLO, MLKL, CASP8, RIPK1, TNF, and FADD. A correlation heatmap of these genes is presented in Figure 2B, and it indicates that more positive associations were found among these NRGs. STUB1, KLF9, and BNPI3 tended to be negatively associated with other genes. The GO enrichment analysis indicated that these NRGs were mainly enriched in necroptotic progression, the extrinsic apoptotic signaling pathway, DNA-binding transcription factor activity, tumor necrosis factor receptor superfamily binding, and ubiquitin protein ligase binding and activity (Supplementary Figure 1A). The KEGG analysis indicated that these genes were involved in the TNF signaling pathway, RIG-I-like receptor, NF-kappa B, NOD-like receptor, FoxO, Toll-like receptor, IL-17, MAPK, and adipocytokine signaling pathways (Supplementary Figure 1B).

Figure 2. The landscape of necroptosis regulators. (A) Protein-protein interaction among necroptosis regulators (Blue: Top 10 hub genes; Red: non-hub genes). (B) Heatmap showed the expression profiling of necroptosis regulators between tumors and normal samples. (C) Heatmap indicated the correlations among necroptosis regulators in pan-cancer. (D) The prognosis role of necroptosis regulators in pan-cancer using log-rank method.

The single-nucleotide variation (SNV) analysis indicated that the variant classification was mainly missense, nonsense, frame shift, splice, and frame shift insertions. Most SNV classes included a C>T mutation. The 10 most mutated genes were ATRX BRAF, IDH1, EGFR, CDKN2A, ALK, HADC9, CASP8, GATA3, and FLT3 (Supplementary Figure 2A). BRAF showed higher mutation frequencies in SKCM and THCA. The mutation frequencies of ATRX and IDH1 were the highest in LGG. EGFR mutations were mainly found in GBM, UCEC, SKCM, LGG, LUAD, and STAD. The 10 most mutated genes were found in similar cancer types (Supplementary Figure 2B, 2C). The following genes showed the highest copy number variations: MTC, ID1, ZBP1, SPATA2, CD40, EGFR, HDAC9, FASLG, BRAF, TRIM11, TERT, TNFSF10, OTULIN, and TRNFRSF1A. These genes showed a high CNV percentage in each cancer (Supplementary Figure 3A). The Spearman correlation analysis indicated that CNV was correlated with mRNA expression (Supplementary Figure 3B). We then explored the differential expression levels of the NGRs between tumor and normal samples. TERT, CDKN2A, PLK1, FASLG, TNFRSF21, TNF, ZBP1, LEF1, HSPA4, and IDH2 were highly expressed in most tumor tissues, while TARDBP, KLF9, BACH2, CFLAR, TSC1, FLT3, and AXL were expressed at low levels in most cancers (Figure 2C). We further explored the correlation of these NGRs with clinical outcomes by performing log-rank tests. NGRs played protective roles in KIRC, SARC and SKCM and risk-enhancing roles in LGG, THYM, LIHC, LUAD, UVM, and ACC. NGRs did not seem to be associated with the prognosis of certain cancers, such as COAD, DLBC, UCS, and ESCA (Figure 2D). Finally, we evaluated the methylation differences in each cancer. NGRs showed different methylation levels in each cancer, particularly COAD, PRAD, BRCA, KIRP, ESCA, PAAD, HNSC, UCEC, LUAD, BLCA, LIHC, LUSC, KIRC and THCA (Supplementary Figure 4A). Overall, the methylation level was found to be negatively associated with mRNA expression level in each cancer (Supplementary Figure 4B).

Establishment of the NPS and its correlation with clinical outcomes

We calculated the NPS of each sample in each cancer and found that DLBC exhibited the highest NPS and that the NPS of UVM was the lowest (Figure 3A). We then compared the NPSs in paired tumor and normal samples. The results showed that the NPS was significantly higher in the ESCA, HNSC, KIRC, KIRP, LUAD, LUSC, SARC, and STAD samples than the normal samples (Figure 3B–3I). However, the NPS was significantly lower in LIHC samples than in normal samples (Figure 3J).

Figure 3. Establishment of NPS and its clinical relevance in pan-cancer. (A) Expression level of NPS in pan-cancer. (B–J) Significant expression levels of NPS between tumor and normal samples (HNSC, KIRC, KIRP, LUAD, LUSC, SARC, STAD, LIHC). (K) Association between NPS and overall survival in pan-caner using univariate cox regression (Left) and log-rank method (Right). (L) Association between NPS and disease-specific survival in pan-caner using univariate cox regression (Left) and log-rank method (Right).

Then, we investigated the correlation of NPS levels with four kinds of clinical outcomes.

For OS, the NPS was associated with a favorable prognosis in SKCM, while the NPS was found to be a risk factor for PAAD, UVM, THYM, TGCT and KIRP (Figure 3K). The Kaplan–Meier analysis indicated that a higher NPS indicated favorable OS in BRCA, DLBC, HNSC, KIRC, OV, READ, and STAD (Supplementary Figure 5). In contrast, the NPS was correlated with unfavorable survival in GBM, KICH, LGG, LUAD, and THYM (Supplementary Figure 5). For DSS (Figure 3L and Supplementary Figure 6), DFI (Supplementary Figures 7A, 8) and PFI (Supplementary Figures 7D–7L, 9), we found similar results. In summary, a high NPS was indicative of a good prognosis in SKCM, ESCA, OV, BRCA, and CHOL (Figure 3C), while a high NPS was an indicator of a poor prognosis in PAAD, KIRP, UVM, TGCT, BRCA and THYM.

Correlations of necroptosis with signaling pathways and cell death in cancers

To explore the correlation of the NPS with signaling pathways, we performed a GSVA of necroptosis in cancers (Figure 4A). We found that the NPS was significantly and positively associated with the interferon gamma response, the inflammatory response, the IL6/JAK/STAT3 signaling pathway, allograft rejection, the complement system, the IL2/STAT5 signaling pathway, the interferon alpha response, TNF-alpha signaling via NFK-beta, and the TGF beta signaling pathway, which were mainly involved in immune response-related signaling. Some pathways were related to cell growth and apoptosis, including the p53 signaling pathway, and certain transduction pathways (PI3K/AKT/mTOR signaling pathway and KRAS signaling up) were positively correlated with the NPS. The NPS was significantly and negatively associated with MYC target V2 in 14 cancers, fatty acid metabolism and bile acid metabolism in 18 cancers, peroxisome in 14 cancers, DNA repair in 11 cancers, oxidative phosphorylation in 20 cancers, reduced KRAS signaling in 15 cancers, and spermatogenesis in 26 cancers. The MYC target V2 was related to tumorigenesis, and fatty and bile acid metabolism and oxidative phosphorylation, which were related to energy supply in tumors. KRAS signaling downregulation was part of a tumor inhibition pathway, and anticancer inhibitors targeting the KRAS signaling pathway have been developed [15]. DNA repair suppressed the antitumor effect of treatment. These results indicated that a high NPS indicates the inactivation of oncogenes, a reduction in energy supply, and the inhibition of tumor progression. However, the NPS showed positive associations with these signaling pathways in OV, TGCT, THYM, and UVM. These results indicated that NPS may exhibit different patterns in different cancers.

Figure 4. GSVA of NPS and its correlations with cell death. (A) GSVA of NPS. (B, C) Correlations of NPS with pyroptosis and ferroptosis regulators.

Necroptosis is a type of programmed cell death. We explored the correlations of NPS with pyroptosis, ferroptosis, and autophagy in each cancer. For pyroptosis, the NPS showed was significantly and positively associated with certain key molecules of pyroptosis, including CASP1, CAPS4, CASP3, IL18, IL1B, and GSDME [16], which are components of classical and nonclassical pyroptosis pathways (Figure 4B). Some molecules, such CHMP2A, showed negative associations with the NPS in cancers. ACLS4 is a key molecule that promotes ferroptosis [17], and the NPS was positively associated with ACLS4. Moreover, the NPS showed a negative association with GPX4 (an anti-ferroptosis-related molecule) in most cancers except PRAD, TGCT, and LAML (Figure 4C). In autophagy, LC3 is the key molecule [18]. We found that the NPS was negatively associated with LC3 (MAP1LC3A) expression in most cancers. However, a positive association was found with OV. Almost all autophagy genes were positively associated with the NPS in OV and DLBC (Supplementary Figure 10).

Correlations of necroptosis with immune cell infiltration in cancers

GSVA indicated that NPS was related to many immune-related signaling pathways. We further evaluated the correlations of NPS with immune status. Our results indicated that the NPS was positively associated with immune and ESTIMATE scores in all cancers and stromal scores in most cancers (Figure 5A). No significant correlations were found between the NPS and stromal score in DLBC and THYM. In contrast, the NPS was negatively associated with tumor purity in call cancers. Next, we assessed the correlations of the NPS with immune-related pathways, stromal/metastasis-related signaling pathways, and DNA repair-related pathways. These pathways were used to assess the tumor microenvironment [19]. We found that the NPS was positively associated with immune checkpoints, effector CD8 T cells, antigen-processing machinery, EMT2, pan-F-TBRs, and EMT3 in cancers (Figure 5B). These results indicated that the NPS potentially reflects the tumor microenvironment status.

Figure 5. Associations with NPS with immune status. (A, B) Associations between NPS and tumor microenvironment in each cancer. (C, D) Associations between NPS and immune cell infiltrations in ImmuCell AI and TIMER2.

Then, we evaluated the correlations of the NPS with 22 human immune cells in ImmuCellAI and TIMER2. Our results indicated that the NPS was positively associated with many immune cells including T cells (iTreg, Th1, Tc, Tfh, and Tex cells), cell infiltration score, macrophages, DCs, nTegs, Tr1 cells, and CD8 T cells, NK cells, and CD4 T cells in various cancers, while the NPS was negatively correlated to naïve CD 4 T cells, Th17 cells, CD8 naïve T cells, and neutrophils in all cancers (Figure 5C). Finally, we explored the correlations between the NPS and immune cells using the EPIC, CIBERSORT, ABS, QUANTISEQ, and XCELL methods. We found that the NPS was positively associated with M1 and M2 macrophages, mast cell resting, monocytes, neutrophils, activated NK-cells, naïve CD8+ and CD4+ T cells, memory B cells, and myeloid dendritic cells (Figure 5D). These results further proved that NPS shows promise as an index that can predict the immune status.

Associations between necroptosis and molecular features and immune regulation

Microsatellite instability (MSI) and tumor mutation burden (TMB) are indicators of genome instability and mutation, respectively. We first evaluated the correlations of the NPS with MSI and TMB. The results indicated that the NPS was positively associated with MSI levels in COAD and GBM and was negatively associated with MSI levels in HNSC, PRAD, SKCM, TGCT, and DLBC (Figure 6A). NPS was positively associated with TBM in COAD but was negatively associated with TBM in TGCTs (Figure 6B).

Figure 6. Associations between NPS and molecular features and immune regulators. (A, B) Associations between NPS with MSI and TMB in each cancer. (C–F) Associations between NPS and immune regulators: MHC genes, immune-regulated and checkpoint-related genes, chemokines, and chemokine receptors.

Necroptosis is a process by which an organism attacks and kills the body’s cells by leveraging the immune system. We further explored the correlations of the NPS with certain important immune regulation gene sets. Major histocompatibility complex (MHC) molecules are involved in antigen recognition during immune responses. Human MHC is called human leukocyte antigen (HLA). The NPS was significantly and positively associated with HLA genes, including TAP1, HLA-E, HLA-DRA, B2M, TAP2, HLA-DPA1, and HLA-DOA in all cancers (Figure 6C). However, the NPS was negatively associated with HLA-G, an immune tolerance molecule that enables tumor cells to escape killing and lysis by NK cells and CTL cells, which is a previously unknown mechanism of tumor cell escape from immune surveillance. Next, we evaluated the correlations of the NPS with immune checkpoint and immune regulation genes in cancers. The results indicated that the NPS was highly correlated with immune checkpoint genes such as PD-L1 (CD274, Figure 6D), programmed cell death 1 ligand 2 (PDCD1LG2), tumor necrosis factor-related immune genes, and T lymphocyte-related immune genes (CD86, CD48, CD80, CD28, CD40, CD27, and CD70; Supplementary Figure 11). We also noticed that V-Set domain containing T-Cell activation inhibitor 1 (VTCN1) was negatively correlated with necroptosis levels in most cancers. These results suggested that the necroptosis level was closely related to most immune checkpoints. Chemokines and their receptors can induce the directional migration of immune cells, which is a necessary condition for the occurrence and completion of the immune response. We then evaluated the associations between necroptosis levels and chemokines and chemokine receptor genes in cancers (Figure 6E, 6F). Chemokines (CCR1, CCR2, CCR3, CCR4, CCR5, CCR6, CCR7, and CCR8) and receptor genes (CXCL1-17, CCL1-14, CCL16-26) were positively associated with necroptosis levels in various cancers (Figure 6). CCL27 and CCL15 showed negative correlations with necroptosis levels in certain cancers. These associations further indicated that the necroptosis level was highly correlated with immune regulation.

Necroptosis and immunotherapy

Our comprehensive analysis revealed that necroptosis was involved in immune infiltration and immune regulation in cancers. We then evaluated the effect of necroptosis on immunotherapy. In bladder cancer (GSE13507), we found that patients with a high NPS exhibited a lower progression ratio than those with a low NPS (Figure 7A), and the Kaplan–Meier curves indicated that patients with a high NPS showed favorable progression-free survival (P<0.04, Figure 7B). In metastatic urothelial cancer (GSE176307), the CR/PR response rate of immunotherapy was higher in the high-NPS group than in the low-NPS group (56% vs. 14%, P<0.001, Figure 7C). The Kaplan–Meier analysis indicated that the high-NPS group presented with a higher OS (P=0.033, Figure 7D) and PFS (P=0.031, Figure 7E). However, in KIRC (PMID: 32472114) and urothelium carcinoma (GSE32894), the high-NPS group showed a worsened PFS (KIRC: P=0.0042, Figure 7F) and disease-free survival (urothelium carcinoma: P=0.02, Figure 7G). These results indicated that NPS can potentially be a marker for predicting the prognosis of patients who receive immunotherapy.

Figure 7. Impacts of NPS on immunotherapy. (A, B) Associations between NPS and immune response and PFS in advanced non-small cell lung carcinoma patients who receiving immunotherapy. (C–E) Associations between NPS and immune response, OS and PFS in metastatic urothelial cancer patients who receiving immunotherapy. (F) Associations between NPS and PFS in KIRC patients who receiving immunotherapy. (G) Associations between NPS and PFS in urothelium carcinoma patients who receiving immunotherapy.

Necroptosis and chemotherapy sensitivity

A previous study has reported that chemotherapy can exert a treatment effect by inducing necroptosis. We therefore evaluated the correlation between necroptosis regulators and certain small-molecule compounds (30 kinds of drugs) in cancers using the GDSC and CTRP databases (Figure 8A, 8B). We observed that IDH1, TNFRSF21, SLC39A7, EGFR, APP, SQSTM1, PANX1, IDH1, AXL, TNFRSF1A, STAT3, STUB1, TLR3, and FADDBNIP3 induced chemotherapy resistance to common drugs. In contrast, ALK, BACH2, BCL2L11, BRAF, CYLT, DNMT1, IDH2, LEF1, MAP3K7, MYC, MTCN, SIRT1, TARDBP, TERT, TSC1, USP22, and FLT3 enhance the sensitivity of these drugs in cancers. More details can be found in Supplementary Tables 3, 4. These results provide potential therapeutic targets for chemotherapy in cancers.

Figure 8. Necroptosis and chemotherapy sensitivity. (A) Necroptosis regulators and GDSC drug sensitivity. (B) Necroptosis regulators and CTRP drug sensitivity.