The available RNA-seq transcriptome data, gene mutation information and clinicopathological information were directly downloaded from the TCGA database (https://portal.gdc.cancer.gov/). There were 476 samples with gene transcriptome data (41 normal and 435 tumor), 347 patients with mutation information, and 400 patients with available survival data. TCGA-Assembler 2 (Wei et al., 2018) was used to obtain Level 3 methylation data (37 normal and 279 tumor) from the TCGA Methylation 450 k Bead chip. Then the function CalculateSingleValueMethylationData of TCGA-Assembler two was used to calculate the average methylation level of each gene. The methylation levels of genes were scored using β values ranging from 0 to 1 (unmethylated to totally methylated). The external validation cohort was from a GEO dataset of gene expression arrays [GSE39582 (N = 550)] (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582). “Deseq2” package of R was used to normalize the raw RNA-seq data and identify differentially expressed genes (DEGs) between normal and cancer groups (Anders and Huber, 2010).

Through the integration of gene expression and DNA methylation datasets, “MethylMix” package of R software was employed to recognize DNA methylation-driven genes (Gevaert, 2015). There are three steps. Firstly, a linear regression model was built to estimate the association between gene methylation and gene expression, and genes with a significant inverse relationship (p value < 0.01) were selected. Such genes were defined as transcriptionally predictive genes. Secondly, a beta mixture model was created to determine the methylation states of each gene. Thirdly, the Wilcoxon rank tests were computed to compare the methylation levels between each methylation state and normal tissue samples. Differential genes were collected. Lastly, genes that were both transcriptionally predictive and differential were selected as DNA methylation-driven genes.

Metascape (Zhou et al., 2019) was used to perform GO Biological Processes enrichment and KEGG pathway analysis of the genes identified by MethylMix. Only terms with p < 0.01, a minimum overlap of 3 and an enrichment factor >1.5 were considered significant.

Initially, the genes identified by MethylMix were applied to a univariable Cox regression and a LASSO Cox regression. In univariable Cox regression analysis, genes with p < 0.05 were selected. In LASSO regression, genes were screened 1,000 times, and if specific genes were detected more than 700 times, they were regarded as candidates. The intersection of univariable Cox regression and LASSO regression (nine genes) were brought into a multivariable Cox regression analysis. The linear combination of the regression coefficient derived from the multivariable Cox regression model (β) multiplied by its mRNA level generated a prognostic risk score with nine genes.

Using X-tile software (Camp et al., 2004) to determine the optimal cut-off value, we divided patients into low- and high-risk groups. Then the Kaplan–Meier curves and time-dependent receiver operating characteristic curve (tdROC) analysis were performed to evaluate the predictive accuracy of the nine-gene signature risk model for OS. The survivalROC function in “survivalROC” package of R was used to draw tdROC and the span argument was set as 0.0657. A subgroup analysis was performed by dividing the patients based on clinicopathological characteristics.

The significance of the risk score model and other traditional clinicopathological characteristics to predict OS was evaluated by univariable Cox regression analysis. Then, we used multivariable Cox regression analysis and stepwise regression method to distinguish significant predictive factors, from which we built a nomogram predictive model. The concordance index (C-index) was computed to evaluate the accuracy of the nomogram. The prognostic risk value of each patient was calculated using the nomogram and tdROC curve analysis was used to further validate the predictive performance of the nomogram. The span argument of survivalROC function was set as 0.0657. The survival estimation and ROC curve of patients above were analyzed with the Kaplan–Meier method by “survival”, “survivalROC” and “plotROC” packages of R.

In the validation phase, we verified the nomogram by using another CRC dataset, GSE39582 in the GEO. The discrimination of model was assessed by the AUC of tdROC and concordance index. The span argument of survivalROC function was set as 0.0616. The calibration of model was visualized in calibration plot and quantified by the slope of the calibration line.

All robustly expressed genes (nearly 18,000 genes) were used for WGCNA analysis. The analysis was performed as described previously (Langfelder and Horvath, 2008). Briefly, we first screened the best soft threshold by “WGCNA” package (Langfelder and Horvath, 2008) of R to ensure a scale-free network. In the co-expression network, genes with high absolute correlations were clustered into the same module. Furthermore, module eigengenes (MEs) were defined as the first principal component of each gene module and the expression of MEs was considered as a representative of all genes in a given module. The correlation between MEs and clinical trait was calculated to identify the clinical key module. GO and KEGG enrichment analysis of the key module resulted from WGCNA were performed by Metascape. ModuleMembership (MM) was defined as the Pearson’s correlation between gene expression and ME. Gene significance (GS) was defined as the Pearson’s correlation between gene expression and certain clinical trait. The cut-off criteria was set as |MM|>0.8, |GS|>0.3 to identify hub genes with high connectivity in the clinical key module. Gene Expression Profiling Interactive Analysis (GEPIA) database (Tang et al., 2017) (http://gepia.cancer-pku.cn/) is an interactive website application that includes the transcriptome data in TCGA and GTEx projects and integrate them in a widely accepted process. We inputted the hub genes into the GEPIA and validated these hub genes.

Differential expression analysis was first performed on all genes to analyze the samples with high and low nine-gene score using “DESeq2” package of R. Enrichment analysis to determine the signaling pathways in which the differentially expressed genes are involved was then carried out by using GSEA method based on the GO Biological Processes and Hallmark gene sets with the “clusterProfiler” (Yu et al., 2012) package of R. A false discovery rate (FDR) less than 0.25 and an absolute value of the normalized enrichment score (NES) greater than 1 were defined as the cutoff criteria.

In the gene mutation analysis, information on genetic alterations was directly obtained from the TCGA database portal, and the quantity and quality of gene mutations were analyzed by using the “Maftools” package (Mayakonda et al., 2018) of R.

The ssGSEA by “GSVA” package (Hänzelmann et al., 2013) of R was introduced to quantify the relative infiltration of 28 immune cell types in the tumor microenvironment. Feature gene sets for each immune cell type were obtained from a recent publication (Charoentong et al., 2017). The gene sets include 782 genes for predicting the abundance of 28 TIICs in individual tissue samples. The following 28 types of immune cells include: activated B cells (Ba), activated CD4⁺ T cells (CD4⁺ Ta), activated CD8⁺ T cells (CD8⁺ Ta), activated dendritic cells (DCa), CD56bright natural killer cells (CD56⁺ NK), CD56dim natural killer cells (CD56⁻ NK), central memory CD4⁺ T cells (CD4⁺ Tcm), central memory CD8⁺ T cells (CD8⁺ Tcm), effector memory CD4⁺ T cells (CD4⁺ Tem), effector memory CD8⁺ T cells (CD8⁺ Tem), eosinophils, gamma delta T cells (γδT), immature B cells (Bim), immature dendritic cells (DCim), mast cells, myeloid-derived suppressor cells (MDSC), memory B cells (Bm), monocytes, natural killer cells (NK), natural killer T cells (NK T), neutrophils, plasmacytoid dendritic cells (DCp), macrophages, regulatory T cells (Tregs), follicular helper T cells (Tfh), type-1 T helper cells (Th1), type-17 T helper cells (Th17), and type-2 T helper cells (Th2). The relative abundance of each immune cell type was represented by an enrichment score in ssGSEA analysis. The ssGSEA score was normalized to unity distribution, for which zero is the minimal and one is the maximal score for each immune cell type.

The bilateral log-rank test was used in the Kaplan–Meier survival examination to compare OS between different subgroups. Wilcoxon rank-sum test was used to evaluate the differences between two groups. All statistical analyses were conducted with R software (version 4.0.3). All statistical tests were two-sided, and p values less than 0.05 were considered statistically significant.

The study flowchart describing the process is shown in Figure 1. To identify DNA methylation-driven genes in CRC, we performed MethylMix analysis on 273 tumor samples and 37 normal tissue samples. A total of 705 genes were identified and a mixture model of each gene was constructed (Supplementary Figure S1). Their methylation levels were visualized by heat map (Figure 2A). Metascape analysis showed the first 20 clusters of enriched sets (Figure 2B). For biological process (BP), these genes showed enrichment in pattern specification process, lipid catabolic process, multicellular organismal homeostasis, cell fate commitment and so on. The KEGG pathway data were enriched in peroxisome and PPAR signaling pathway.

We included 400 patients with RNA-seq data and complete clinical information from the TCGA database for subsequent analysis. The patients’ characteristics are summarized in Table 1. Univariable Cox regression analysis was performed for 705 genes and 40 genes with p < 0.05 were selected. Besides, the same 705 genes were brought into LASSO regression analysis and 10 genes were selected (Figures 3A,B). The intersection of above two gene sets, including nine genes (LINC01555,GSTM1,HSPA1A, VWDE, MAGEA12,ARHGAP4,PTPRD, ABHD12B,TMEM88) were eventually screened out as prognosis-related genes. Supplementary Figure S2 showed the correlation between gene expression and corresponding methylation level of nine genes. Then, multivariable Cox regression analyses were performed, and a nine-gene model was constructed according to their expression levels and coefficients (Figure 3C). The formula was as follow: risk score = (−0.421*LINC01555)+(0.168*GSTM1)+(0.241*HSPA1A)+(0.581*VWDE)+(0.107*MAGEA12)+(0.153*ARHGAP4)+(−0.168*PTPRD)+(0.281*ABHD12B)+(0.157*TMEM88). Based on the formula, we calculated the risk score for each patient and used X-tile software to identify the optimal cut-off value (2.51). Those with a risk score over the cut-off value, 172 patients in total, were classified as the high-risk group, while the remaining 228 patients were classified as the low-risk group (Figure 3D). Intuitively, more people died in the high-risk group than in the low-risk group (Figure 3E). The expression profiles of nine genes in all patients and the corresponding risk group were presented in the form of heat map (Figure 3F). The Kaplan–Meier analysis of all patients demonstrated that the high-risk group had a significantly shorter OS (P = 2e-08) (Figure 3G). The AUC of this nine-gene risk assessment model was 0.745 (95% CI 0.634–0.856), 0.708 (95% CI 0.615–0.801) and 0.721 (95% CI 0.621–0.821) at 1-, 3- and 5- years, respectively (Figure 3H).

We further tested the effect of the model in different subgroups divided by some clinicopathologic factors (age, TNM stage, lymph node metastatic and distant metastatic). The K-M survival analysis showed the same prediction trend in all subgroups, proving that the model has certain reliability and practicability in evaluating prognosis (Figure 4).

Risk score, age, T stage, lymph node metastatic and distant metastatic were selected as significant predictive factors after univariable regression analysis (Table 2). Stepwise regression analysis showed that after removing the factor lymph node metastatic, AIC of multivariable Cox regression went down from 751.65 to 750.64, so it was further ruled out. Therefore, four independent risk factors (risk score, age, T stage and distant metastatic) were ultimately retained to build a visualized and applicable nomogram (Figures 5A,B). The C-index of the model is 0.81 (95% CI: 0.754–0.868). According to Nomogram, each variable was assigned a corresponding score, as a result we could calculate a total score for each patient. The AUC of nomogram was 0.815 (95% CI 0.712–0.917), 0.794 (95% CI 0.708–0.881) and 0.802 (95% CI 0.720–0.884) at 1-, 3- and 5- year, respectively (Figure 5C). The results revealed that the predicted survival possibility by the nomogram was close to the actual survival situation. Furthermore, our nomogram model performed better than the model only using significant clinicopathological factors (Supplementary Table S1).

The nomogram established above was further validated in the GEO dataset GSE39582 (Figures 6A–E). Totally 550 patients with complete basic clinical information were included for analysis (Table 1). The calibration curves for the predicted possibility of 1-, 3- and 5-years survival displayed obvious concordance between the predicted results and the actual observations in the GEO dataset (Figures 6B–D). The calibration slope was 1.56 (95% CI 1.31–1.82) at 1-year survival, 1.21 (95% CI 1.05–1.36) at 3-years survival and 1.00 (95% CI 0.74–1.26) at 5-years survival. Similar to the performance in TCGA cohort, the AUC of tdROC was 0.788 (95% CI 0.684–0.892) at 1-year survival, 0.743 (95% CI 0.679–0.807) at 3-years survival and 0.714 (95% CI 0.647–0.780) at 5-years survival (Figure 6E). The concordance index of the nomogram was 0.722 (95%CI 0.683–0.761). To sum up, the results showed that our model performed well in the validation cohort.

We used two approaches to comprehensively explore molecular characteristic of different risk groups. Firstly, we performed WGCNA analysis in the TCGA cohort. In current study, the power of β = 10 (scale-free R ² = 0.74) was selected as a soft threshold parameter to ensure a scale-free network (Figures 7A,B). The co-expression gene modules were then constructed and divided into 28 meaningful modules via the average linkage hierarchical clustering (Figure 7C). Blue module was found to have the highest association with 9-gene risk model (Cor = 0.28, P = 1e-08 for risk score; Cor = 0.2, P = 6e-05 for risk group) (Figure 7D). There were a total of 1,577 genes in the blue module. We then used Metascape to analyze these genes. Top 20 clusters of functional enriched sets were presented (Figure 8A). Based on the cut-off criteria (|MM|>0.8, |GS|>0.3), twelve hub genes (LZTS1, TIE1, STARD8, VEGFC, KCNE4, ADAMTS1, AFAP1L1, ITGA5, BCL6B, MMRN2, CAVIN1 and CCBE1) with high connectivity in the clinical significant module were identified (Figure 8B). The expression of all twelve genes showed significant correlation with risk score (Supplementary Figures S3A–L). The hub genes were then analyzed in the GEPIA database. For LZTS1, VEGFC, KCNE4, ITGA5, CCBE1 and CAVIN1, K-M survival analysis revealed that a higher expression was meaningfully associated with a worse prognosis (Figures 8C–H).

Besides, GSEA was performed in the two risk groups. A false discovery rate (FDR) less than 0.25 and an absolute value of the normalized enrichment score (NES) greater than 1 were defined as the cut-off criteria. The gene sets of the high-risk group were enriched in tumor progression and metastasis-related pathways and inflammatory response-related pathways (Figure 9A), while the gene sets of the low-risk group were enriched in DNA integration, epithelial structure maintenance related pathways (Figure 9B).

Next, we analyzed gene mutations to further understand different molecular subgroup associated with the nine-gene risk score (Figures 9C,D). The mutation rates of APC, TP53, TTN, KRAS, SYNE1, MUC16 and PIK3CA were higher than 20% in both groups. The mutation of the RYR2, ABCA13 and PCLO genes was more common in the high-risk group, while the mutation of FAT4, OBSCN and ZFHX4 genes was more common in the low-risk group.

Given that the features of tumor-infiltrating immune cells (TIICs) are correlated with the development and progression of cancer and may influence the prognosis, including CRC, we then used ssGSEA method to analyze the infiltration of various immune cells in tumor samples. We found that cells executing pro-tumor suppression (such as MDSCs, regulatory T cells, immature dendritic cells) were more abundant in the high-risk group (Figure 10A). Additionally, the correlation between the risk and common ICPs, i.e., PDCD1 (PD1), CD274 (PDL1), CTLA4, LAG3, HAVCR2 (TIM3) and TIGIT, was analyzed. Consistent with the higher infiltration of immunosuppressive cells in high-risk group, several common ICPs showed a significantly higher expression in the high-risk group (Figure 10B). The immune landscape and clinicopathological characteristics of different risk groups are presented as heat map in Figure 10C.