The colon adenocarcinoma-associated mRNA sequences and clinical data used in this study were from TCGA ¹ and GEO (GSE39582) databases. Patients who met the following selection criteria were included: (A) the primary site of the tumor must be the colon, not colon metastasis; (B) The pathological type of colon cancer must be adenocarcinoma; (C) Available gene expression data; (D) Available survival information. The clinical data details of two train datasets used in this study are shown in Supplementary Table S1. The list of apoptosis-related genes was extracted from molecular signatures database ² (Liberzon et al., 2015). The search strategy includes the following keywords: “apoptosis” and “Homo sapiens.” Finally, 2,737 apoptosis-related genes were listed. PFI (progression-free interval), the most reliable prognostic index recommended in the TCGA database, was selected as the end point of disease progression. PFI is the period from the date of diagnosis until the date of the first occurrence of a new tumor event, which includes progression of the disease, locoregional recurrence, distant metastasis, new primary tumor, or death with tumor (Liu J. et al., 2018). Supplementary Figure S1 showed the flowchart for the present research.

The differentially expressed genes (DEGs) between colon cancer samples (379 cases) and normal samples (39 cases) in the TCGA dataset were identified using the limma package in R (p-value < 0.05, |log2FC| >1) (Ritchie et al., 2015). The apoptosis-related genes obtained from the molecular signatures database were intersected with the DEGs of colon cancer to obtain apoptosis-related DEGs. The DEGs related to apoptosis were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using Metascape ³ to determine the potential molecular mechanism of gene expression profile (Zhou et al., 2019).

Weighted gene co-expression network analysis (WGCNA) is a bioinformatics analysis method used to describe the gene association patterns between different samples (Langfelder and Horvath, 2008). Clustering of genes with similar expression patterns promotes web-based hub gene screening, which can be used to identify candidate biomarkers (Tang J. et al., 2018; Chen et al., 2021). The DEGs related to apoptosis were analyzed by the WGCNA package in R. First, the Pearson correlation was calculated from the expression data, and the relationship matrix was established. Then, the correlation matrix was transformed into an adjacency matrix, and the power exponential weighting was introduced to construct the scale-free network. Next, based on the adjacency matrix, the Tom matrix was established to calculate the Tom difference degree (1-TOM) between genes, and then the gene module was established. Finally, the results of the module were visualized, and the target module was selected for analysis. In this study, 127 genes of the most significantly related module (Turquoise module) were selected for further exploration. Twenty apoptosis-related hub genes significantly related to progression were screened by the R package (“survival” and “surviviner”) for further analysis.

From 20 apoptosis-related hub genes, the genes independently related to the prognosis of colon cancer were screened by multivariate Cox regression analysis. The prognostic index of apoptosis related genes was determined using the following formula: ARGPI = Expression gene 1 * Coefficient + Expression gene 2 * Coefficient + + Expression gene * coefficient. The value of the ARGPI was defined as the risk score of each patient. According to the median risk score, the patients in the two databases were divided into the high-risk group and low-risk group for follow-up evaluation. The prognostic evaluation ability of the ARGPI was evaluated by plotting the K-M survival curve of TCGA and GEO cohort (p < 0.05, log rank test).

Gene set enrichment analysis (GSEA) is a method to evaluate microarray data at the level of gene set (Subramanian et al., 2005). In the GSEA software, the gene set (c2. cp.kegg.v7.4. symbols.gmt) was used as the reference gene set to explore the important signal pathways enriched in high-risk subgroups and low-risk subgroups (p < 0.05 and FDR <0.25). In the analysis of gene mutation, information on mutations related to colon cancer samples was obtained from the TCGA database. Somatic variants in tumor patients in the high-risk and low-risk subgroups were comprehensively analyzed using the Maftools package in R (Mayakonda et al., 2018).

CIBERSORT ⁴ is a tool for deconvoluting the expression matrix of immune cell subtypes based on the principle of linear support vector regression. This method utilized expression profile data of colon cancer to estimate the infiltration of 22 types of immune cells in different risk subgroups (Newman et al., 2015). ssGSEA (single-sample gene set enrichment analysis) is used to study the enrichment level of immune-related functions in tissue expression profiles (Barbie et al., 2009). The enrichment levels of 29 immune-related functional features in different risk subgroups were analyzed using the GSVA package in R.

To improve the prognostic risk stratification of colon cancer and assist in clinical diagnosis and treatment, a nomogram model integrating the ARGPI and clinicopathological features was constructed. It is used as a quantitative tool to predict the prognosis of patients with colon cancer. The effectiveness of the model was evaluated using a calibration curve.

Differences among variables were analyzed using an independent t-test and a chi-squared test. Univariate survival was assessed using Kaplan–Meier survival analysis and log rank test. Univariate Cox regression and multivariate Cox regression were used to analyze the effects of the ARGPI and clinical characteristics on prognosis. The data were statistically analyzed by R (version 4.1.1) and SPSS software (version 25.0). Bilateral p < 0.05 was considered statistically significant.

A total of 7,859 DEGs were obtained by comparing colon adenocarcinoma samples with normal tissue samples, of which 5,432 genes were upregulated and 2,427 genes were downregulated (Figure 1A). Among these DEGs, 655 genes were related to apoptosis. Compared with normal samples, 392 genes were upregulated and 263 genes were downregulated in tumor samples (Figure 1B). Functional enrichment analysis identified GO terms and KEGG pathways that were enriched by 655 apoptosis-related DEGs (Figures 1C,D).

To identify apoptosis-related hub genes, WGCNA was used for candidate genes (n = 655). As shown in the left figure of Figure 2A, the horizontal axis is the soft threshold and the vertical axis is the evaluation parameters of scale-free network. The higher the value of evaluation parameters, the more consistent the network is with the characteristics of scale-free network. In the left graph, the horizontal line indicates that the threshold value is 0.90. As shown in the right figure of Figure 2A, the relationship between soft threshold (power) and mean connectivity shows that the optimal soft-thresholding power was seven based on the scale-free network. Finally, six modules are identified by hierarchical clustering and optimal soft threshold capability (Figure 2B). As shown in Figure 2C, the turquoise module with the highest absolute correlation with colon cancer was selected as the candidate module for further analysis, which contains 127 genes. Significantly enriched GO terms and KEGG pathways of turquoise modular genes are shown in Supplementary Figure S2. The expression profiles and clinical prognosis information of 127 genes in the turquoise module were obtained. Through K-M survival analysis, the expression of 20 apoptosis-related hub genes was closely related to the progression of colon cancer patients (Supplementary Figure S3; p < 0.05, log rank test).

Multivariate Cox regression analysis was performed on 20 apoptosis-related hub genes to determine independent prognostic genes. As shown in Figures 3A,B, five genes (FAS, VWA5A, SPTBN2, PCK1, and TIMP1) significantly affect the prognosis of patients with colon cancer. Then, we constructed the ARGPI for all colon cancer patients. Through the Risk formula = Expression level of FAS * (−0.2326) + Expression level of VWA5A * (−0.3434) + Expression level of SPTBN2 * 0.3991 + Expression level of PCK1 * (−0.1833) + Expression level of TIMP1 * 0.3263. Taking the median risk value as the critical value, low-risk patients had better prognosis than high-risk patients (p < 0.001, log rank test) (Figure 3C). Then, the role of the prognostic model was verified by using the GSE39582 (n = 517) colon cancer data set. As shown in Figure 3D, the prognosis of patients in the low-risk group was significantly better than that in the high-risk group, which was concordant to the results of the TCGA dataset (p = 0.018, log rank test).

The distribution of clinicopathological features of 379 colon cancer patients in the TCGA cohort in different risk subgroups is shown in Figure 4A. In Figure 4B shows significant differences in the distribution of tumor stages between the two risk subgroups. Univariate Cox regression analysis showed that stage and risk score were significantly correlated to colon cancer prognosis (Figure 4C). Multivariate Cox regression analysis confirmed that risk score was an independent prognostic factor after adjusting for other clinicopathological factors (Figure 4D).

GSEA analysis was performed to identify pathways enriched in different risk subgroups. High-risk subgroup enriched CELL_CYCLE and FOCAL_ADHESION and other related pathways (Figure 5A). The low-risk subgroup was enriched in DRUG_ METABOLISM_ CYTOCHROME_ P450, RETINOL_ Metabolism and other related pathways (Figure 5B). Supplementary Table S2 lists the detailed results of GSEA, in which apoptosis related pathways are also enriched.

Next, we analyzed gene mutations to gain further biological insight into the nature of subgroup differences. We found that the mutation count in the low-risk subgroup was higher than that in the high-risk subgroup. Missense mutations were the most common mutation type, followed by nonsense mutations and frameshift deletions. We then identified the top 20 genes with the highest mutation rate in the risk subgroup (Figures 5C,D). The mutation rates of APC, TP53, TTN, KRAS, SYNE1, PIK3CA, MUC16 and FAT4 in the two groups were >20%.

To analyze the tumor microenvironment status in different risk subgroups, we first used CIBERSORT algorithm to estimate the proportion of 22 types of invasive immune cells in different risk subgroups. We found that macrophages M0, macrophages M1 and T cells regulatory (Tregs) were higher in the high-risk subgroup, while B cells naïve, plasma cells, T cells CD4 memory resting and eosinophils were more abundant in the low-risk subgroup (Figure 6A). We then applied the ssGSEA score to quantify the enrichment level of 29 immune features in each colon cancer sample. The results showed that there were more macrophages in the high-risk subgroup and more IDCs and NK_cells correlation signals in the low-risk subgroup (Figure 6B). We further studied the effects of immune and molecular functional differences between different risk subgroups on survival. As shown in Figures 6C,D, patients with more IDCs and NK_cells have better prognosis.

To explore the relationship between apoptosis characteristics and treatment, we analyzed the chemotherapy data and molecular characteristics of colon cancer in the GEO (GSE39582) databases. As shown in Figure 7A, the high-risk group is more likely to progress after adjuvant chemotherapy than the low-risk group. Therefore, we observed the effect of the expression of five core genes on the progression after chemotherapy. As shown in Figures 7B,C, when the VWA5A gene is highly expressed, tumor patients make less progress after chemotherapy, suggesting that the upregulation of this gene promote chemotherapy sensitivity. When the TIMP1 gene is highly expressed, tumor patients progress more after chemotherapy, suggesting that the upregulation of TIMP1 gene reduce chemotherapy sensitivity. Figure 7D again shows that when the tumor shows upregulation of the VWA5A gene and downregulation of the TIMP1 gene, the tumor has better chemosensitivity. Stratified analysis showed that the apoptosis-related prognostic model was still a clinically and statistically significant prognostic model in patients with pMMR status, BRAF wild type, and KRAS mutation (Figures 7E–H).

Receiver operating characteristic curve predicted by two-year progression showed that the ARGPI had higher prognostic accuracy than other clinicopathological features (Figure 8A). Then, to provide clinicians with a quantitative method to predict the possible risk of cancer progression, we constructed a nomogram based PFI, which integrates the ARGPI and four clinicopathological risk factors (Figure 8B). The prediction accuracy of the nomogram is as follows. The calibration curve of the nomogram shows that it does not deviate from the reference line and does not need to be recalibrated [Figure 8C (1 year), Figure 8D (3 years), Figure 8E (5 years)].