Genetic and copy number alterations of NRGs in colon cancer
Our research contains 67 NRGs. Biological analysis of the frequency of somatic mutations in these NRGs showed that the incidence of mutations in colon cancer was significantly increased (Figure 1a). 46.12% of the 399 colon cancer patient samples had mutations. Of these, BRAF was the most mutation-prone gene (15%), followed by ATRX (8%), while four NRGs (TNFRSF1B, CDKN2A, TNFSF10 and ID1) were not mutated.
We next examined copy number alterations of somatic in these NRGs and identified copy number variations (CNV) that were prevalent in all 67 NRGs. Of these, CNV was generally increased for MYC, ID1, ZBP1, CD40, SPATA2, and STAT3, while CNV was reduced for TARDBP, TNFRSF1B, TLR3 and FAS (Figure. 1b). Figure 1c shows the position of CNV alterations in the NRGs on the respective chromosomes. We further explored the genetic differences between colon cancer tissues and normal tissues and showed that the expression levels of most NRGs showed significant differences. NRGs with low copy number variation frequencies, such as FAS and TLR3, were lowly expressed in colon cancer samples compared to normal samples, whereas NRGs with high copy number variation frequencies, such as MYC and SPATA2, were highly expressed in colon cancer samples (Figure. 1d), and based on the results, it was hypothesized that CNV may affect the gene expression of NRGs. Nevertheless, there were also differences in the expression of some NRGs between colon cancer samples and normal samples that were not associated with CNV gain or deletion. Therefore, gene expression levels are not only regulated by CNV10. Other factors can also be involved in regulating gene expression, such as DNA methylation11. Our study suggests that NRGs have a potential function in the development of colon cancer.
Identification of necroptosis subtypes in colon cancer
For the purpose of analysis, the expression pattern of NRGs in the biology of colon cancer, we used 1147 samples from 3 downloaded colon cancer datasets (TCGA-COAD, GSE39582 and GSE17536) for subsequent analysis in our research. Table S3 demonstrates the prognostic value of the 67 NRGs in the univariate Cox regression and Kaplan-Meier analysis in the colon cancer sample, filtered at p < 0.05. Subsequently, the necrotic apoptosis network demonstrated the linkage between NRGs and their prognostic value for patients with colon cancer (Figure 2a and Table S4).
To further characterize the expression pattern of NRGs in colon cancer, according to the expression of 67 NRGs, we applied a clustering algorithm to classify the colon cancer samples. The results demonstrate that dividing all samples into two subtypes (k = 2), A (n = 521) and B (n = 624), produces the smallest experimental error (Figure 2b). PCA analysis showed striking differences in the expression of NRGs between the A and B subtypes (Figure 2c). The survival curve demonstrates that the survival rate of subtype A samples is lower than that of subtype B samples over a long period of time (Figure 2d). In addition, significant differences in the NRGs expression of the two subtypes and clinical features can be seen in the NRGs expression heat map (Figure 2e).
GSVA enrichment and Immune cell infiltration analysis in distinct subtypes
GSVA enrichment analysis revealed significant differences in immune pathway enrichment between subtypes A and B. Significantly different pathways included natural killer cell-mediated cytotoxicity, T and B cell receptor signaling pathway, Toll-like receptor signaling pathways, ECM receptor interaction, antigen processing and presentation and leukocyte trans endothelial migration (Figure 3a; Table S5). Subsequently, we applied the CIBERSORT algorithm to analyze the relationship between A and B subtypes in 22 immune cell types. The results showed that the expression of the two subtypes A and B was significantly different in the immune cell species (Figure 3b). Among the 19 immune cells that showed significant differences (p < 0.001), subtype B showed lower expression of immune infiltration compared to subtype A.
Identification of gene subtypes based on DEGs
To further investigate the biological behavior of necroptosis, we screened 599 necroptosis subtype-related DEGs and performed biological function and pathway enrichment analysis (Figure 4a-b). These necroptosis subtype-related DEGs are markedly enriched in modules of biological processes associated with leukocytes (Figure 4a). The main enrichment pathways of DEGs are also immune-related as shown by the results of KEGG analysis (Figure 4b). This suggests a significant role for necroptosis in immune regulation. We then applied a univariate Cox analysis to screen 599 necroptosis subtype-related genes for prognostic value (p < 0.05), and we finally identified 301 genes associated with survival for subsequent analysis (Table S6). To further explore the relationship between necroptosis and biological processes, consensus clustering analysis was applied to divide the samples into two gene subtypes, A and B, based on the expression of 301 prognostic genes (Figure 4c). Survival curve results show that patients with gene subtype B have lower OS than those with gene subtype A (p < 0.001; Figure 4d). Furthermore, most of the genes were highly expressed in the necrotrophic apoptotic gene subtype B. (Figure 4e). As expected, the two necroptosis-related gene subtypes appeared to differ significantly in the expression results of NRGs (Figure 4f).
Construction and assessment of prognostic models and nomogram
Based on the necroptosis subtype-related 301 prognostic genes, we constructed NRG_score. Figure 5a shows the distribution of samples across two subtypes of necroptosis, two gene subtypes and high and low risk groups. In the first place, we applied the “caret” package from R to divide the sample randomly and equally into training (n = 568) groups and test (n = 567) groups. LASSO analysis and multivariate Cox statistical analysis were used to select the best prognostic features for 301 necroptosis prognostic genes. We then performed a LASSO regression model analysis and finally screened for 13 genes associated with OS, including PNMA1, CXCL13, PDGFRL, GREM1, SLC2A3, FABP4, HOXC6, NT5E, IL7R, CXCL11, IL1R2, NOX1 and CKMT2 (Figure S2). Based on the results of the analysis, the NRG_score model equation is as follows:
Risk score = (0.2101* expression of PNMA1) + (−0.1008* expression of CXCL13) + (−0.4855* expression of PDGFRL) + (0.1179* expression of GREM1) + (0.2148* expression of SLC2A3) + (0.1958* expression of FABP4) + (0.1632* expression of HOXC6) + (0.1925* expression of NT5E) + (−0.2543* expression of IL7R) + (−0.1110* expression of CXCL11) + (−0.2111* expression of IL1R2) + (−0.1159* expression of NOX1) + (−0.1153* expression of CKMT2).
Based on the risk scores of the samples, we visualized the distribution of each sample. In typing according to necroptosis-related DEGs, gene subtype B had a significantly higher NRG_score than gene subtype A (Figure 5b). But in typing according to NRGs subtype A had a higher NRG_score than subtype B (Figure 5c). Next, we divided the patients in the training group into two groups based on the median NRG_score, with patients with a risk score below the median NRG_score being classified as low-risk group (n = 284) and those with a risk score above the median NRG_score being classified as high-risk group (n = 284). Figure 5d-e shows that as NRG_score increases, the number of patients in the high-risk group gradually increases and the chance of recurrence of colon cancer gradually rises. Survival curve results suggest that OS is significantly higher in patients with low NRG_score than in those with high NRG_score (p < 0.001; Figure 5f). In addition, we produced ROC curves for NRG_score and Figure 5g shows that the accuracy of predicting patients by NRG_score at 1, 3 and 5 years was 0.764, 0.759 and 0.735 respectively.
Subsequently, we selected data from the test (n = 567) group to validate the prognostic accuracy of the NRG_score. In the same way as in the previous study, we divided the patients in the training group into a high-risk group (n = 272) and a low-risk group (n = 295) based on the median value of NRG_scores. The results showed that the number of patients in the high-risk group tended to rise as the NRG_score increased (Figure 6a-b). The results of the survival curve reveal that the prognosis of patients in the low-risk group is significantly better than that of the high-risk group (p < .001; Figure 6c). The predictive accuracy of the ROC curve for the validation group was 0.707, 0.679 and 0.692 for 1-, 3- and 5 years respectively, indicating that the NRG_score is an excellent predictor of survival in patients with colon cancer (Figure 6d).
To improve the clinical usefulness of the NRG_score, we constructed a nomogram by using risk scores and clinical data (Figure 6e). The factors affecting prediction include, gender, age, risk score and patient stage. Subsequently, the calibration plots show that the nomogram we constructed are functionally similar to the ideal prognostic model (Figure 6f).
Monitoring of immune cell checkpoints between the high and low risk groups
We applied the CIBERSORT algorithm to observe the relationship between risk score and immune cell abundance. Testing of 22 human immune cells showed the risk score was positively correlated with M2 macrophages, M0 macrophages and neutrophils and negatively correlated with M1 macrophages, plasma cells, activated memory CD4 + T cells, testing memory CD4 + T cells and CD8 + T cells (Figure 7a). Low-risk group significantly correlated with immune score; high-risk group correlated with stromal score (Figure 7b). Subsequently, we assessed the association of 13 genes in the NRG_score with immune cell infiltration. It is clearly observed that most of the genes are significantly associated with immune cells (Figure 7c).
CSC Index, mutation and drug susceptibility analysis of NRG_score
We explored the potential correlation between NRG_score and CSC index in colon cancer. Based on the results of the study, we conclude that the CSC index (R =−0.24, p < 0.001) decreases as the NRG_score increases, indicating that colon cancer cells with lower risk scores have lower levels of cell differentiation (Figure 8a). We then explored mutations between the two risk groups in the TCGA-COAD mutation dataset. The gene with the highest mutation rate is APC, and all of the genes shown in the graph are mutated to varying degrees (Figure 8b-c). We next explored differences in the sensitivity of chemotherapeutic agents for colon cancer in high and low risk groups. We observed that patients in the high-risk group were more sensitive to shikimic, while those in the low-risk group were more sensitive to gemcitabine and Cyclopamine In conclusion, the findings point to a significant correlation between NRGs and drug sensitivity (Figure 8d-f).