Gene expression and clinically relevant information were obtained from the Cancer Genome Atlas (TCGA) and Gene-Expression Omnibus (GEO) databases. Three independent CRC cohorts (TCGA-COAD, GSE17538, and GSE39582) were collected for this study. We acquired gene expression data in FPKM (Fragments per Kilobase Million) format from TCGA and converted them to TPM (Transcripts per million) format, downloaded normalized matrix files from the GEO database, and removed batch effects between these datasets using the R package SVA. We downloaded clinical information of the three datasets, and detailed clinical information was provided in Supplementary Tables 1-3 . We extracted survival status and survival time from clinical information, and excluded data with a follow-up time fewer than 31 days and duplicate data. Additionally, we downloaded the tumor immune subtype and stemness score files from the UCSC database for subsequent analysis.

By searching necroptosis-related literature, we identified 67 necroptosis-related genes (17). To investigate whether NRGs play a role in the development of CRC, we compared NRGs expression in normal and CRC tissues, and we subsequently performed Kaplan-Meier survival analysis and univariate COX analysis to determine the prognostic value of NRGs.

Unsupervised clustering analysis was performed to identify different NRG-related molecular subtypes based on the gene expression profiles of 67 NRGs, and patients were classified for follow-up studies. We used a consensus clustering algorithm to determine the number of clusters (K) and their stability, using the K value corresponding to the flattest cumulative distribution function (CDF) curve as the determined number of clusters. The R package consusclusterplus was utilized to perform the above analysis steps. Principal component analysis was used to determine whether there was good discrimination between subtypes. To investigate the biological process differences among NRG-related molecular subtypes, we performed gene set variation analysis (GSVA) based on gene set c2.cp.kegg.v7.4.symbols.

The above consensus clustering analysis divided patients into different NRG-related molecular subtypes and DEGs between NRG subtypes were identified using the R package limma with screening criteria of fold-change ≥ 1.5 and adjusted P-value < 0.05.

To explore the biological functions of differentially expressed genes between NRG-related molecular subtypes, we performed functional enrichment analysis using the clusterProfiler package, including GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis. The screening criteria were P < 0.05 and FDR < 0.05. Univariate Cox regression analysis was performed on DEGs, and genes associated with CRC prognosis were selected for unsupervised clustering of samples to identify gene subtypes. The number and stability of gene clusters were determined by the consensus clustering algorithm.

ssGSEA (single sample GSEA) is an implementation method mainly proposed for a single sample that cannot do Gene Set Enrichment Analysis (GSEA). We used ssGSEA to explore the infiltration of various types of immune cells and immune-related functions in CRC samples. The R packages limma, GSEABase and GSVA were used to perform the above analysis. CIBERSORT is a deconvolution algorithm. The R package CIBERSORT can calculate the proportion of different types of immune cells in a sample based on the Leukocyte signature matrix (LM22) containing 547 reference genes. The ESTIMATE algorithm infers tumor purity and cell density from the RNA sequencing data of the samples. Based on the R package ESTIMATE, we assessed the immune and stromal content (immune and stromal scores) in each sample. In addition, we also collected 47 immune checkpoint genes to compare their expression in different subtypes.

We first randomly divided all samples into two groups: the training group (n=570) and the testing group (n=570). The training group was used to learn the sample characteristics and construct the signature, and the testing group was used as a validation cohort to test the prediction performance of the signature.

Based on the DEGs screened between different NRG subtypes, we performed a univariate Cox analysis to further identify the DEGs associated with CRC prognosis. Next, we performed a 10-fold cross-validation lasso (least absolute shrinkage and selection operator) regression and Cox proportional hazards regression analysis. Lasso regression was run for 1000 cycles to obtain the gene combination with the smallest cross-validation error, followed by Cox analysis and model construction. The final risk score was calculated as follows:

where Exp is the expression value of each gene in the signature, and W is the coefficient of multiple cox regression analysis for each gene. The median risk score in the training set was used as the cut-off value, and patients were divided into a high-risk group and a low-risk group. Survival curves were plotted to evaluate the overall survival (OS) of two groups, and the predictive performance of this signature for 1-, 3-, and 5-year survival was evaluated using the receiver operating characteristic (ROC) curve. Subsequently, the predictive NR-signature was applied to the testing group and the whole samples for validation.

Many prognostic models have been constructed and demonstrated good predictive performance in CRC. To investigate whether the signature based on necroptosis is effective and complementary to existing studies, we compared our model with other phenotype-based models, including those based on immune, autophagy, pyroptosis, ferroptosis, aging, and metabolism genes (14, 18–27). Based on the TCGA database, we calculated and compared the concordance index (C-index) of several signatures, and we plotted the ROC curve and survival curve of each signature to visualize the predictive effect.

We performed GSEA analysis on all samples to enrich the respective involved signaling pathways in high and low risk groups based on the KEGG gene set with clusterProfiler package. To explore the tumor mutational burden (TMB) in CRC, we calculated the total number of non-synonymous mutations in the samples from TCGA database. The R package maftools was used to draw the oncoprint of gene mutations in high and low risk groups. We also investigated the correlation of risk groups with clinical indicators to see if there were differences in the proportion of high- and low-risk samples in different tumor stages and immune subtypes. To explore the relationship between risk status and immune cells, we assessed the status of immune cell infiltration and immune function in samples using multiple methods including ssGSEA and the CIBERSORT algorithm. In addition, we also used various methods such as XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT to analyze the correlation between risk scores and infiltrating immune cells based on TCGA samples. We also calculated the TIDE score and Dysfunction score of the TCGA samples based on the online website TIDE (http://tide.dfci.harvard.edu/) to evaluate the immune escape and immunotherapy of the samples in the high and low risk groups.

We compared the differential expression of genes in normal and cancer tissues in pan-cancer. The correlation of gene expression with immune microenvironment scores (including immune and stromal scores, ESTIMATE score, and tumor purity), Stemness Score (including RNAss: RNA expression-based and DNAss: DNA methylation-based), and immune subtypes were calculated on a pan-cancer scale.

We obtained cancer and normal tissue samples from eight CRC patients who underwent curative resection at Tongji Hospital of Tongji Medical College (Wuhan, China). The MolPure^® Cell/Tissue Total RNA Kit (Yeasen) was used to extract total RNA, which was then reverse transcribed with the Hifair^® III 1st Strand cDNA Synthesis Kit (gDNA digester plus, Yeasen) in accordance with the manufacturer’s protocols. The target sequence was amplified with real-time PCR with the Hieff^® qPCR SYBR Green Master Mix(Low Rox Plus, Yeasen). The cycling parameters used were 95°C for 5 min, 95°C for 10 s, and 60°C for 30 s for 40 cycles. Melting curve analyses were performed, and Ct values were determined during the exponential amplification phase of real-time PCR. The 2^–ΔΔCt method was used to determine relative fold changes between tumor tissues and normal tissues as the following equation: 2 ^–ΔΔCt (ΔΔCt = ΔCt^tumor – ΔCt^normal). The primer sequences were listed in Supplementary Table 5 .

All statistical analysis and graphing were performed by R-4.1.2. Wilcoxon rank-sum test and Student’s t test were used for comparison between the two groups. Kruskal-Wallis test was used for comparison among more than two groups of samples. The Kaplan-Meier method was used to plot survival curves for prognostic analysis, and the log-rank test was used to determine the significance of differences. The correlation test was performed using Spearman correlation analysis and distance correlation analysis. Comparisons of composition ratios among groups were performed by chi-square test. All heatmaps were plotted by the R package pheatmap. All statistical P values were two-tailed, and P < 0.05 was used as the truncated value.

We first compared the expression of 67 NRGs in 473 CRC tissues and 41 normal tissues based on the TCGA database, of which 52 genes were differentially expressed ( Figure 1A ). Among them, MLKL was highly expressed in cancer tissues, while RIPK1 and RIPK3 were lowly expressed in cancer tissues. To provide a more comprehensive landscape of the expression of NRGs, we jointly analyzed TCGA and GTEx database through GEPIA to expand the sample size. NRGs with differential expression exceeding 1.5 times were shown in Supplementary Figure 1 . A total of 1247 tumor samples were obtained after combining TCGA and GEO data, and we further analyzed whether there was an impact of NRGs on the prognosis of CRC patients. We performed univariate Cox analysis ( Figure 1B ) and Kaplan-Meier survival analysis ( Supplementary Figure 2 ), respectively ( Supplementary Table 4 ). The analysis results showed that the expression of some NRGs was related to the survival of patients. The high expression of RIPK1 was associated with poor prognosis, while the high expression of RIPK3 represented a good prognosis, and the MLKL expression was not related to the survival time of patients.

To further explore the expression characteristics of NRGs in colorectal cancer, we qualitatively classified patients based on the expression profiles of 67 NRGs using the R package ConensusClusterPlus. By consensus clustering algorithm, a cluster number (K value) of 2 was the best choice to classify the whole samples into NRG cluster A (n=695) and NRG cluster B (n=552) ( Figures 2A, B ). The principal component analysis demonstrated significant differences of necroptosis gene profile between the two clusters ( Figure 2C ), and the heatmap showed the expression of NRGs in two NRG clusters ( Figure 2E ). Kaplan-Meier survival analysis indicated that patients with NRG cluster A had better overall survival (P = 0.010, Figure 2D ). To explore the potential biological change between distinct NRG clusters, we applied GSVA enrichment analysis, which showed that NRG cluster B was significantly enriched in immune-related pathways, including B_CELL_RECEPTOR_SIGNALING_PATHWAY, T_CELL_RECEPTOR_SIGNALING_PATHWAY, LEUKOCYTE_TRANSENDOTHELIAL_MIGRATION, CHEMOKINE_SIGNALING_PATHWAY, NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY, ANTIGEN_PROCESSING_AND_PRESENTATION, NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY, TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY and REGULATION_OF_ACTIN_CYTOSKELETON ( Figure 2F ). To further investigate the role of NRGs in TME of colorectal cancer, we applied the ssGSEA algorithm to evaluate the association between NRG subtypes and immune cell subpopulations, and the results showed significant differences in the infiltration of most immune cells between different clusters ( Figure 2G ). Immune cell infiltration was more abundant in NRG cluster B, including Activated.B.cell, Activated.CD4.T.cell, Activated.CD8.T.cell, Activated.dendritic.cell, Eosinophil, Gamma.delta.T.cell, MDSC (Myeloid-derived suppressor cell), Macrophage, Monocyte, Natural.killer.T.cell, Natural.killer.cell, Neutrophil, Regulatory.T.cell (Treg). In addition, the expression of immune checkpoint genes was generally upregulated in cluster B.

To reveal the potential biomolecular characteristics underlying different NRG subtypes, we identified 702 NRGcluster-related DEGs based on R package limma, followed by GO and KEGG functional enrichment analyses ( Figures 3A, B ). The results showed that these DEGs were mainly enriched in biological processes such as immune-related processes, cytokines-related processes, chemokines-related processes, and tumor signaling pathways, suggesting that necroptosis may play an important role in tumor development and immune regulation. We then performed a univariate Cox regression analysis to identify 361 NRGcluster-related DEGs with prognostic value, which were defined as necroptosis-related signature genes (NRSGs). To further verify the regulatory mechanism of necroptosis, we performed unsupervised clustering based on 361 NRSGs to divide the samples into two gene subtypes: gene cluster A and gene cluster B ( Figures 3C, D ). According to the Kaplan-Meier curve, gene cluster B showed a better prognosis (P < 0.001, Figure 3E ). The heatmap showed the expression distribution of 361 NRSGs in the two gene clusters ( Figure 3F ). Furthermore, significant differences in the expression of NRGs were observed between distinct gene clusters, which was consistent with the expected results of the necroptosis regulatory pattern ( Figure 3G ).

We assessed the potential role of NRPGs in the tumor microenvironment of CRC. Based on the ESTIMATE algorithm, gene cluster A exhibited higher estimate scores, stromal scores, immune scores, and lower tumor purity ( Figure 4A ). This suggested that the TME of gene cluster A had a higher content of stromal cells and immune cells. Through CIBERSORT algorithm, we further observed that plasma cells, CD8⁺ T cells, naive CD4⁺ T cells, resting and activated CD4 memory cells, Tregs, resting NK cells, and dendritic cells were more abundant in gene cluster B, while M0 macrophages, M1 macrophages, M2 macrophages, and neutrophils were significantly higher in gene cluster A ( Figure 4B ). In addition, we integrated immune cell infiltration information of TCGA samples from multiple existing databases. Through comparison, overall immune cell infiltration was more abundant in gene cluster A ( Figure 4D ). Finally, we compared the expression levels of immune checkpoint genes between the two clusters. The results showed that the vast majority of immune checkpoint genes were highly expressed in subtype A, including the most well-known genes like LAG3, CTLA4, ICOS, TIGIT, PDCD1, CD274 (PD-L1), BTLA ( Figure 4C ). Previous studies have shown that macrophage infiltration is associated with poor prognosis in CRC, whereas infiltration of NK cells, NKT cells, and γδT cells predicts good prognosis (28). Although cluster A generally had higher infiltration of immune cells and stromal cells, it had higher infiltration of macrophages and lower infiltration of CD8⁺ T cells. Besides, the high expression of immune checkpoint genes can cause the depletion of immune cells. These factors may ultimately lead to a worse prognosis of subtype A.

A prognostic signature was constructed based on NRSGs. We randomly divided the samples into a training group (n = 569) and a testing group (n = 569). LASSO and Cox analyses were subsequently performed on the 361 NRPGs to further construct the optimal prognostic model ( Supplementary Figure 3A ). Finally, we included 25 genes into the model, and the signature was constructed as follows: NRG riskscore=-0.364*PALLD+0.215*VSIG4+0.329*ANTXR2-0.552*PTPRM+0.269*IHH+0.595*SIGLEC1-0.186*CXCL13-0.509*CEBPA+0.185*PLK2+0.190*EGR2-0.339*IGFBP5-0.567*UBE2L6+0.275*AIFM3+0.303*ARMCX2+0.259*IGFBP3+0.192*ACE2+0.306*HOXC6-0.077*SPINK1+0.146*FABP4-0.296*THBS4-0.175*PHGR1-0.173*MMP12-0.161*CKMT2-0.167*PLCB4+0.278*VIP.

According to the median value of the risk score of the training group, we divided the patients into a high-risk group and a low-risk group. With the increase in risk scores, the proportion of death among patients also increased. The heatmap showed differential expression of 25 mRNAs in the high-risk and low-risk groups ( Figure 5A ). K-M survival analysis showed that the prognosis of the low-risk group was significantly better than that of the high-risk group ( Figure 5D ). To test the prediction performance of signature, we plotted ROC curves with AUC values of 0.813, 0.791, and 0.799 in the training group at 1, 3, and 5 years ( Figure 5E ). To further validate the signature, we applied it to the testing group and all samples, and the results showed a trend consistent with the training group ( Figures 5B–E ). And we also re-validated the prediction effect of the model in three independent datasets, respectively ( Supplementary Figure 3B, C ).

We compared the predictive performance of the risk model and other clinical indicators. In the prediction of 1-year OS, the effect of the stage was slightly better than that of the risk model, while in the prediction of 3- and 5-year OS, the prediction ability of the risk model was better than other indicators ( Figures 6A–C ). Similarly, we also re-compared each indicator in the three datasets, and the results showed good predictive performance of the risk model ( Supplementary Figure 4 ).

In addition, various types of prediction models have been constructed in CRC. To further test the effect of NR-signature, we compared the prediction accuracy of various models based on the TCGA database. By calculating, we found that the C-index of the NR-signature was higher than other types of models ( Figure 6D ). Besides, we also plotted the ROC curve and K-M survival curve of each model ( Figure 6E ; Supplementary Figure 5 ). The results demonstrated that the NR-signature had good prediction performance and can be used as a complement to existing types of prediction models.

To investigate the underlying biological processes in the high- and low-risk groups, we performed GSEA analysis based on the KEGG geneset. The results showed that tumor-related pathways were enriched in the high-risk group, including ECM RECEPTOR INTERACTION, FOCAL ADHESION, and PATHWAYS IN CANCER. In addition, CALCIUM SIGNALING PATHWAY, which regulates cell survival and death, was also enriched in the high-risk group ( Figure 7A ). In the low-risk group, DNA REPLICATION, OXIDATIVE PHOSPHORYLATION, PEROXISOME, and PROTEASOME were significantly enriched ( Figure 7B ). Next, we analyzed the distribution of somatic mutations in the high and low-risk groups based on TCGA-COAD data, and the top ten mutated genes in both groups were APC, TP53, TTN, KRAS, SYNE1, PIK3CA, MUC16, FAT4, ZFHX4, and RYR2 ( Figures 7C, D ). APC and TP53 are the most widespread and representative mutations in colorectal cancer, and detecting the mutation status of these genes can help improve the accuracy of diagnosis and guide individualized treatment. Moreover, several studies have focused on reactivating TP53 function to exert its anti-tumor effects (29, 30). Among the mutated genes, the mutation rate of KRAS was significantly higher in the high-risk group than that in the low-risk group (51% vs. 37%, P = 0.009). According to studies, the prevalence of KRAS mutation in colorectal cancer is about 40%, and KRAS mutation suggests poor prognosis (31–33). The 8th AJCC Cancer Staging Manual specifies that the level of evidence for KRAS mutation as a prognostic and predictive factor is grade I and II, respectively (33). In addition, KRAS gene mutation indicates poor response of colorectal cancer patients to anti-EGFR targeted therapy (32, 33). The higher frequency of KRAS mutations in the high-risk group in our analysis is consistent with these findings. Additionally, we assessed the correlation of risk scores with clinical indicators. The proportion of patients with high and low risk varied by stage, with fewer high-risk patients than low-risk patients within stage I-II and more high-risk patients within stage III-IV ( Figure 7E ). Moreover, the proportion of patients with high and low risk differed across immune subtypes ( Figure 7F ).

Based on the CIBERSORT algorithm, the heatmap showed a close correlation between the 25 genes in NR-signature and immune cells ( Figure 8A ). We explored the relationship between risk score and immune cells of TCGA samples according to immune cell infiltration information from multiple databases. Correlation analysis showed that the overall risk score was negatively correlated with immune infiltration ( Figure 8B ). We then assessed immune cell infiltration and immune function in all samples using different methods. ssGSEA analysis showed that the infiltration of CD8⁺ T cell, DC (dendritic cell), pDC (plasmacytoid dendritic cells), Tfh (follicular T-helper-cell), Th1, TIL (tumor-infiltrating lymphocytes) and Treg was lower in the high-risk group, and APC (antigen-presenting cell) co-inhibition, cytolytic activity, inflammation-promoting, T cell functions such as co-inhibition and co-stimulation were attenuated ( Figures 8C, D ). The CIBERSORT algorithm showed a similar trend, with lower immune infiltration and immune function overall in the high-risk group, and a significant decrease in CD8⁺ T cells, cytolytic activity, inflammation-promoting, and TIL ( Figures 8E, F ). We also performed a survival analysis of these immune cells and functions, and low-grade infiltration of the aforementioned important immune cells and functions was associated with poor prognosis ( Supplementary Figure 6 ). The expression levels of multiple immune checkpoint genes also differed between high and low-risk groups ( Figure 8G ). Finally, we calculated TIDE scores and Dysfunction scores based on samples from the TCGA database, and the results showed that patients in the high-risk group scored higher, predicting a higher potential for immune escape and immune dysregulation ( Figure 8H ).

To further validate the functions of the genes included in NR-signature, we explored these genes in pan-cancer based on the TCGA database. First, we compared the expression levels of 25 genes in cancer tissues and normal tissues, and the results showed that 25 genes showed significant expression differences in multiple cancer types ( Supplementary Figures 7, 8 ). In addition, we also investigated the association of genes with immune microenvironment score and Stemness Score. Among them, VSIG4, SIGLEC1, and CXCL13 were significantly positively correlated with estimate scores, stromal scores, and immune scores in pan-cancer, while negatively correlated with tumor purity ( Supplementary Figure 9 ). In general, genes were negatively correlated with RNAss. Subsequently, we analyzed the relationship between genes and immune subtypes in COAD, and boxplots indicated differences in gene expression among different immune subtypes ( Supplementary Figure 10 ). We also re-analyzed the association of 25 genes with immune microenvironment score and stemness score in COAD, and the expression of most genes showed a strong positive correlation with the immune microenvironment score ( Supplementary Figure 11 ).

To further validate the functions of the genes included in the model, we first performed KM survival analysis of these genes based on the TCGA database ( Supplementary Figure 12 ), and selected 10 survival-related genes for expression validation. RT-qPCR results based on eight pairs of colorectal cancer tissues showed that eight of the 10 genes were significantly differentially expressed between cancer and normal tissues. Among them, MMP12 was highly expressed in colorectal cancer tissues, while ARMCX2, CEBPA, CXCL13, FABP4, HOXC6, SIGLEC1 and VSIG4 were more expressed in normal tissues than in cancer tissues ( Figure 9 ).