The Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) website (https://portal.gdc.cancer.gov/) was used to obtain the expression spectrum data for colon cancer (colon adenocarcinoma, COAD) patients (n = 514), such as the count, and FPKM and TPM values; and patient clinical data (n=430), such as gender, age, TMN stage, and survival prognosis. “Masked somatic mutation” was selected as the somatic mutation data (n=420) and downloaded. The somatic mutations were visualized using maftools package (22) in R to obtain tumor mutation burden (TMB) for per patient. Additionally, the MSI data in the TCGA-COAD patients’ dataset was obtain the tumor mutation burden (TMB) per patient. Additionally, the MSI data in TCGA-COAD patient dataset were obtained from the cBioPortal database (https://www.cbioportal.org). The baseline information of TCGA-COAD patients is provided in Table 1 .

Gene expression data and the clinical characteristics of patients of GSE17536 (23) and GSE39582 (24) were downloaded from the GEO database. The data samples were obtained from Homo Sapiens. The chip platforms were grounded in the GPL570 [HG-U133_plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. After deleting patients lacking clinical information, 177 COAD tissue samples were included in GSE17536, all of which were used in the analysis. GSE39582 included 585 COAD tissue samples. As survival information was not available for 5 samples, 580 were included in the analysis. R’s limma package (25) was used to standardize the two data sets separately. Table S1 shows the information from GEO.

For all samples in the combined dataset and based on 36 necrotizing apoptosis-related genes from previously published literature (26), the necroptosis score (NPs) of every sample in the TCGA-COAD dataset was determined using the R package, GSVA (27), and the ssGSEA method according to the gene expression matrix of the respective sample.

Using the above method, the NPs of each sample was obtained, and the optimal cut-off value was selected in combination with the patient’s survival data. Based on the NPs score, the NPs group was separated into high and low NPs. To identify genes associated with NPs, the DEGs between high NPs and low NPs in TCGA-COAD samples were identified using the R package, limma (25). The screening threshold of the DEGs was set to |log2 fold change (FC)|>1 and adjusted P < 0.05. The results of the difference analysis are presented as a heatmap and a volcano plot.

WGCNA was implemented using the R package, WGCNA (28). First, the weighted values of the calculated correlation coefficients between any two genes were used to generate connections between genes in the network to assemble a scale-free network. A hierarchical clustering tree was then established according to the correlation coefficients. The branches of the cluster tree highlighted various genetic modules, and various colors signified different modules. The module saliency was then calculated. All mRNAs of the sample were input into WGCNA to measure the associations between the two NPs groups and different modules. All genes were recorded in their respective modules. The genes in the respective modules were considered as modular characteristic genes (MEs). The correlation between the NPs values and genes was determined based on the significance of genes. Module membership was determined according to the relevance between module genes and DEG expression profile. The modules of interest were selected using the MS score, and all genes in these modules were recognized to have a high correlation with NPs.

Based on the necroptosis characteristic genes and TCGA-COAD expression data, the k-means method in “ConsensusClusterPlus” R package (29) was used to perform unsupervised cluster analysis to identify the necroptosis subtypes. The concordant clustering algorithm was used to discern the cluster number. The analysis was repeated 1000 times to ensure stability of the category. Principal component analysis (PCA) was conducted for patients with grouped subtypes to determine the differences between samples. Survival analysis was implemented after grouping to determine the influence of various subtypes on prognosis.

The results of differential expression analysis combined with WGCNA analysis were used to acquire the NPs-related characteristic genes. The significantly differentially expressed NPs-related characteristic genes were involved in the model, and the prognostic genes were screened using univariate Cox regression analysis, with a cut-off P value of 0.1. Subsequently, the selected genes were regularized and dimensionally reduced using the least absolute shrinkage and selection operator regression (LASSO) algorithm to further identify prognostic-related genes. Thereafter, the weighted normalized gene expression value of the penalty coefficient acquired by multivariate Cox analysis (STEP method) was used to establish a risk score formula. Using the median risk score, patients were divided into high-risk and low-risk groups.

The above dataset based on TCGA-COAD served as a training set, and the internal test was conducted using the bootstrap method with 1000 re-sampling. Thereafter, the coefficient based on model variables was used to calculate the risk score for each sample in test sets, GSE17536 and GSE39582, using the predict function in the “survival” R package (https://CRAN.R-project.org/package=survival>). Finally, a time-dependent receiver operating characteristic (ROC) curve was plotted. The area under the curve (AUC) was used to reflect the performance of the model.

To acquire genes relevant to the NPs model, DEGs between the high-risk group and low-risk group of TCGA-COAD patients were analyzed using the R limma package. The screening threshold of the significantly different DEGs was defined as |LogFC| > 1 and adj. P value < 0.05. The DEGs were visualized using volcano maps and heat maps.

GO (30) analysis is an approach adopted for massive functional enrichment research, including biological processes (BP), molecular functions (MF), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) (31) is an extensively used database that stores data regarding genomes, biological pathways, diseases, and drugs. The ClusterProfiler package of R (32) was used for GO and KEGG analyses of significant DEGs. The critical value of FDR less than 0.05 indicated significant difference.

To explore the discrepancies in biological processes among different subgroups, GSEA was performed according to the gene expression profiling dataset of COAD patients. The gene set “c2.cp.v7.2.Symbols. gmt” obtained from the MSigDB (33) database was used for GSEA. An FDR < 0.25 indicated statistical significance.

To quantitatively analyze the relative tumor infiltration degree of various immunocytes in COAD, the ssGSEA algorithm was employed to differentiate between highly sensitive and specific phenotypes of various human immune cells in the tumor microenvironment (TME). The algorithm revealed 28 gene sets for labeling various tumor-infiltrating immunocyte types based on a study by Bindea et al. (34). The gene sets comprised various human immunocyte subtypes, such as macrophages, mast cells, etc. Enrichment scores obtained using ssGSEA in R’s GSVA package indicated the degree of infiltration of various immune cell types in every sample. Meanwhile, R’s ESTIMATE package (35) was employed to evaluate the immunological activity of the tumor. ESTIMATE quantitatively analyzes the immune activity of a tumor sample according to its gene expression profile to obtain an immune score per tumor sample. Herein, the discrepancies in immune infiltration features between the two groups of patients with COAD were compared.

To compare the copy number differences between the two groups of TCGA-COAD patients, the TCGAbiolinks package (36) of R was used to obtain the Masked Copy Number Segment information. The downloaded CNV fragments were subjected to GISTIC 2.0 analysis by GenePattern (https://cloud.genepattern.org); default parameters were used for the analysis.

The predictive power of the NPs risk score in combination with clinicopathological characteristics on OS based on univariate and multivariate Cox analysis was used to demonstrate that the NPs risk score in combination with clinicopathological features can be used to estimate patient prognosis. The risk scoring model was then combined with clinicopathological features to establish a nomogram, and the accuracy of the model was reflected by the AUC values under the time-ROC curve. The performance of the rosette was assessed using a calibration curve that compared the predicted values of the rosette with the observed actual values. Testing of the model was carried out using the bootstrap method, and internal re-sampling was performed 1000 times.

The STRING (37) online tool was applied to establish the PPI network. Genes with scores > 0.7, which indicates high credibility, were selected from the STRING database to construct the network model visualized using Cytoscape (version3.7.2) (38). The Maximal Clique Centrality (MCC) of each node was calculated using the cytoHubba plug-in (39) in Cytoscape software.

Information on the miRNA-mRNA interactions was collected from the miRTarBase database (40). The core mRNAs acquired from the PPI analysis were used to predict the miRNAs that might be regulated. The relevant lncRNAs were further predicted based on evidence from the luciferase reporter gene assay. The results of ceRNA analysis were visualized using Cytoscape software. The P values of all hypothetical tests were two-sided, and a p value of less than 0.05 was considered to indicate statistical significance.

Data analysis was performed using R software. version 4.0.2. Independent Student’s t test and Mann–Whitney U test (namely Wilcoxon rank-sum test) were used to estimate the differences between two groups of normally distributed and two groups of non-normally distributed continuous variables, respectively. The χ2 test or Fisher exact test was carried out to determine the difference between the two groups of categorical variables. Survival analysis was carried out using R’s survival package. The Kaplan–Meier survival curve was applied to display survival differences, and log-rank test was performed to compare the differences in survival. Univariate and multivariate Cox analyses were based on the R survival package. LASSO analysis was carried out using the glmnet R package (41).

The whole research design was illustrated in Figure 1 . First, 36 necroptotic genes were extracted from the RNA-seq data of TCGA-COAD and their expression differences were compared between the normal group and tumor group. A total of 30 necroptotic genes were found to be differentially expressed, and only 6 genes were not differentially expressed. Such results suggest that necroptosis may play crucial role in COAD ( Figures 2A, B ). According to the somatic mutation data of TCGA-COAD samples, mutation information was obtained for the 36 necroptotic genes using the maftools package. The necroptotic genes were not found to mutate significantly in COAD patients, except TP53, in which the mutation frequency of TNF, CASP6, and TNFSF10 was less than 1% ( Figure 2C ). In addition, the prognostic status of the 36 genes was analyzed. Only TRAF2, RIPK3, and IPMK genes were found to have significant prognostic differences, and may thus serve as potential prognostic markers ( Figure S1 ).

Based on the above differential necroptosis genes, the NPs of each COAD patient was obtained using the ssGSEA algorithm to represent the necroptosis level of the patient. The optimal cut-off value was determined using prognostic analysis based on the necroptosis score. Thereafter, TCGA-COAD samples were divided into groups with high NPs and low NPs. Survival analysis revealed that patients with low necrosis apoptosis scores had worse prognosis than patients with high necroptosis score (log-rank P < 0.024; Figure 3A ). Subsequently, 766 DEGs were acquired through rigorous analysis, including 380 significantly upregulated and 386 downregulated genes, respectively ( Figures 3B, C ). A gene co-expression network was also established to identify biologically significant gene modules through WGCNA and further identify genes closely related to COAD necroptosis. In this study, 4 modules (except grey module) were obtained for subsequent analysis ( Figures 3D, E ). As shown in Figure 3F , we integrated the difference analysis results with the MEturquoise, MEblue, and MEbrown modules to obtain a total of 209 necrotizing apoptotic characteristic genes.

Based on the above genes with necroptosis characteristics, consistent clustering was employed to cluster LIHC samples. Here, K=2 was selected and two subgroups, subgroup 1 and subgroup 2, were obtained ( Figure S2A ). Dimension reduction analysis was conducted via PCA and the PCA results of the two groups were plotted. Based on the results, the degree of differentiation between the two groups was not obvious, which may be due to the insignificant clustering gene characteristics, resulting in insignificant grouping differences ( Figure S2B ). The prognostic characteristics of subgroup 1 and subgroup 2 was subsequently analyzed using the KM curve; however, no distinct difference in prognosis was found between the two groups, suggesting that clustering subtypes based on all necroptotic characteristic genes could not distinguish differences in the prognosis of patients ( Figure S2C ). However, when the expression distributions of characteristic genes of type 1 and type 2 were compared, necroptotic characteristic genes were found to be significantly differentiated in the two subtypes, suggesting that the classification is of guiding significance for assessing the mechanism of necroptosis, but not suitable for identifying clinical prognostic markers ( Figure S2D ).

Based on the necroptosis characteristic genes, a necroptosis-related risk score system was constructed to quantitatively evaluate the prognostic information of each COAD patient by risk score. First, univariate Cox regression analysis revealed that 23 genes met the screening criteria (P < 0.05). Dimension reduction was analyzed using LASSO ( Figure 4A ). When 5 variables were present, the most stable model was obtained. Multivariate Cox regression analysis revealed that Calbindin 1 (CALB1), Carbohydrate sulfotransferases (CHST13), and solute carrier family 4 member 4 (SLC4A4) were independent prognostic factors ( Figure 4B ). Multivariate Cox analysis was also carried out to obtain the model coefficients of important characteristic genes. Thereafter, the gene expression was multiplied and summed with its coefficients to construct a risk score. The final risk score (necroptosis risk score related to prognosis) was calculated for each sample. In terms of the risk score and gene expression values of patients, a heat map of the risk factors was plotted to show the distribution of the risk score ( Figure 4C ). The time-dependent ROC curve analysis revealed AUC values of 0.684, 0.657, and 0.710 for the 1-, 3-, and 5-year OS, respectively, which indicated that risk score was an ideal predictor of OS in COAD patients ( Figure 4D ). For the external dataset test, the GSE39582 and GSE17536 datasets were employed. After data normalization, the model was tested. For GSE39582, the ROC curve revealed AUC values of 0.682, 0.623, and 0.708 for the 1-, 3-, and 5-year OS, respectively. For GSE17536, the AUC values were 0.665, 0.712, and 0.758 for the 1-, 3-, and 5-year OS, respectively, indicating that the model had been well tested in external datasets ( Figure S3 ).

To determine the role of the necroptosis-related risk model on the evolution of COAD samples, TCGA-COAD patients were divided into high-risk group and low-risk group based on the expression median value of the COAD patient risk model score in TCGA dataset. Subsequently, the DEGs in the two groups of patients were identified. Overall, 317 genes were significantly differentially expressed in COAD patients, among which 185 and 132 genes were significantly upregulated and downregulated, respectively ( Figures 5A, B ).

Functional enrichment analysis was performed using 317 genes identified as significantly DEGs. The GO analysis results revealed that the significant DEGs were related to GO:0044421 extracellular region part, GO:0005576 extracellular region, GO:0042588 zymogen granule, GO:0071752 secretory dimeric IgA immunoglobulin complex, and other functions ( Figure 5C ). KEGG functional analysis suggested that significant DEGs mainly had an impact on nitrogen metabolism, bile secretion, and rheumatoid arthritis, and other pathways ( Figure 4D ). Many pathways were found to be related to immunity, such as the chemokine signaling pathway, WNT signaling pathway, etc. Detailed GO and KEGG results are presented in Tables S2 and S3 .

Based on the results of expression analysis, we continued GSEA and summarized the related results of the pathway database based on C2. KEGG pathway results and GSEA results revealed distinct differences in the activity of the JAK-STAT signaling pathway, WNT signaling pathway, fructose and mannose metabolism, primary immunodeficiency, and nitrogen metabolism ( Figure 5 ). The detailed results of GSEA and the metabolism-related pathways are provided in Table S4 . These findings coincide with those obtained from the KEGG database, and suggest the activation and inhibition characteristics of the high- and low-risk groups.

In terms of the DEGs in the high- and low- risk groups, we aimed to identify the hub gene that played a key role, and its potential molecular interaction mechanism. First, the STRING database was used to analyze the protein interaction mechanism. As shown in Figure 6A , after screening with a confidence of 0.700, the number of PPI nodes (protein) was 233. Further, 81 edges were identified, with an average connection degree of 0.695 for each node. The enrichment statistic P value of the whole PPI network was less than 1.0e-16.

The interacting proteins were further identified as hub genes using the cytoHubba plugin in Cytoscape. After the calculations, MUC5AC, MUC5B, WNT16, WIF1, and other interacting proteins that had the top 8 scores were found ( Figure 6B ). Subsequently, miRNA molecules and lncRNAs that potentially regulate these hub genes were analyzed using miRTarBase database. Finally, a ceRNA regulatory network was established using Cytoscape ( Figure 6C ).

The influence of necroptosis-associated risk scores in patients with TCGA-COAD on their holistic immune characteristics and varying degrees of immunocyte infiltration was analyzed. Patients in the high-risk group were found to have significantly lower immune-related scores (P < 0.001); however, no distinct difference was found in the matrix score ( Figures 7A, B ). ssGSEA was used to appraise the changes and effects of immunological characteristics of COAD tissues during pathogenesis. Through ssGSEA, the relative enrichment scores of 28 different subtypes of immunocytes in the high- and low-risk groups of COAD patients was obtained. Heat maps were generated to illustrate their expression in different patients ( Figure 7C ). Based on the results, the expression abundance of immunocytes in the low-risk group was lower than that in the high-risk group. The correlation analysis results revealed that most immune cell infiltration levels were positively correlated ( Figure 7D ). Differential analysis also revealed distinctions in the infiltration levels of various immune cells between COAD samples in the high and low necrotizing apoptosis-related risk groups. Only effector memory CD8+ T cells, immature dendritic (iDC) cells, and other cells were not found to significantly differ between the groups ( Figure 7E ). A distinct difference was found between the HLA family expression levels and the various immunological targets of the high- and low-risk groups. Further, the immunoactive genes were almost all elevated in the low-risk group ( Figures 7F, G ).

The influence of the necroptosis-related risk score on changes in genetic variation levels, including single nucleotide polymorphism and CNV, in COAD patients was evaluated. The analysis of single-nucleotide mutations in common tumor-driven genes revealed that the high mutation levels were similar or close between patients with high and low scores in the necroptosis-related model ( Figure 8A ). Based on assessments of the frequency of CNV changes, CNV was found to be widely present in high- and low-risk samples. However, no distinct discrepancy in CNV was found between the two groups ( Figures 8B, C ). When TMB and MSI were compared between the two patient groups, no distinct discrepancy in TMB and MSI was found between the high and low risk groups. This result suggests that changes at the genomic level were not significant in the two groups ( Figures 8D, E ).

Based on the significant role of immunotherapy in tumors, the TIDE algorithm was employed to calculate the sensitivity of patients in the high- and low-risk groups to immunotherapy. The TIDE score in the low-risk group was higher than that in the high-risk group, suggesting that the immunotherapy response of the high-risk group might be better than that of the low-risk group ( Figure 8F ).

To further probe the clinical value of necroptosis-related risk score, the clinical characteristics related to the high-risk and low-risk groups, such as the discrepancy in age and TNM stage, were analyzed. Notably, no distinct discrepancy was found in the age of patients in the high-risk group ( Figure 9A ). For gender, the proportion of women in the high-risk group increased ( Figure 9B ). In terms of stage, a significantly higher proportion of advanced patients was identified in the high-risk group ( Figure 9C ). Subsequently, based on the risk scores associated with necroptosis and clinicopathologic features (age and TNM stage), we established a prognostic model for COAD patients ( Figure 9D ) and analyzed the model via 1000 resampling using the bootstrap method. Based on timeROC, the AUC values were 0.798, 0.772, and 0.741 for 1-, 3-, and 5- years, respectively ( Figure 9E ). Calibration curves were generated to present the consistency of the model. A good consistency was found between the model’s estimated 1-, 3-, and 5-year OS and the actual observed OS of patients ( Figure 9F ).