Colon cancer-specific diagnostic and prognostic biomarkers based on genome-wide abnormal DNA methylation

Genomic distribution of hypermethylated and hypomethylated DMPs

To explore the abnormal methylation status of the entire genome, we conducted in-depth studies on the early diagnosis and prognostic evaluation of COAD patients (Figure 1). First, we performed CpG site expression profiling analysis between COAD tumor samples (N = 25) and paired normal samples (N = 25) from The Cancer Genome Atlas (TCGA) cohort. Then, 13716 Hyper-DMPs and 11403 Hypo-DMPs were obtained in these included cohorts. Specifically, unsupervised cluster analysis distinguished these Hyper-DMPs and Hypo-DMPs in 25 paired COAD and normal samples from TCGA (Figure 2A). When we assessed the locations of these Hyper-DMPs and Hypo-DMPs among genomic region types, we observed that Hyper-DMPs were most abundant in Island regions (39%), whereas Hypo-DMPs were mainly distributed in Open sea regions (42%) (Figure 2B). We also determined the enrichment of the DMPs by calculating the ratio of Hyper-DMPs to Hypo-DMPs in each region. The results indicated that Hyper-DMPs were enriched in Island regions (66%; Hyper/Hypo = 5421/2689), whereas Hypo-DMPs occurred more frequently in Open sea regions (58%; Hypo/Hyper = 4882/3588).

More importantly, Hyper-DMPs were mainly located near promoter regions, including TSS1500 (the region 200 to 1500 nucleotides upstream of the transcription start site), TSS200 (the region from the transcription start site to 200 nucleotides upstream of the transcription start site), the 5′ untranslated region (UTR) and the 1^st Exon (Figure 2C). However, Hypo-DMPs were mostly enriched in the Body and the 3′UTR, which occupied a large percentage of the regions, genome-wide. The DMP distribution ratio also indicated that proximal promoter regions were mainly hypermethylated (69%; hyper/hypo = 6245/2738), while the proportions of Hyper- and Hypo-DMPs in the Body and 3′UTR were almost equal (51%; hypo/hyper = 4551/4307). Notably, both Hyper-DMPs and Hypo-DMPs occupied a large proportion of the whole genome, about 3.42%.

Next, we calculated Pearson correlation coefficients to determine the correlation between the DNA methylation of the DMPs and the expression of their corresponding genes (Figure 2D). Among the 17112 DMPs for which both the DNA methylation levels and the corresponding mRNA expression profiles were available, the methylation levels of 6565 Hyper-DMPs and 4112 Hypo-DMPs were significantly associated with the mRNA levels of the corresponding genes (|r| > 0.1, false discovery rate [FDR] < 0.05). When we analyzed the distance between these DMPs and promoter regions, we found that DMPs in or near promoter regions (i.e., in the 1^st Exon, 5′ UTR, TSS200 or TSS1500) were negatively associated with mRNA expression, whereas those outside promoter regions (i.e., in the Body or 3′ UTR) were positively associated with gene expression. Moreover, the DMPs had a higher distribution frequency on chromosomes 7 and 1 than on the other autosomes (Figure 2E). Since the DMPs that significantly altered the expression of their corresponding genes were not limited to promoter regions, we screened the whole genome for potential biomarkers of COAD and constructed a diagnostic prediction model based on the genome-wide Hyper- and Hypo-DMPs.

Identification of COAD-specific methylation biomarkers and construction of a COAD-specific diagnostic model

Next, in order to construct a diagnostic model to distinguish COAD tumor tissues from normal intestinal epithelial tissues and the tumor tissues of nine other cancer types (BLCA, BRCA, CESC, glioblastoma [GBM], head and neck cancer [HNSC], liver cancer [LIHC], lung adenocarcinoma [LUAD], lung squamous cell carcinoma [LUSC] and endometrial cancer [UCEC]), we performed conditional screening and machine learning studies based on the genome-wide DMPs obtained from TCGA above (13716 Hyper-DMPs and 11403 Hypo-DMPs). For the conditional screening, we determined the average β values (a measure of CpG site methylation) of these Hyper-/Hypo-DMPs in all the samples for the nine cancer types in TCGA. Then, we selected DMPs based on an average methylation level difference of at least 0.2 units in COAD. After screening under the above conditions, we obtained 17 Hyper-DMPs and 8 Hypo-DMPs as candidate biomarkers. Further analysis revealed that these 17 Hyper-DMPs and 8 Hypo-DMPs were among the 6565 Hyper-DMPs and 4112 Hypo-DMPs that were significantly associated with the expression of their corresponding genes (log2 |fold change| > 1, FDR < 0.05).

For machine learning, two-thirds of the total tumor and normal samples from the COAD cohort of TCGA (200 COAD and 25 normal samples) were randomly set as the training cohort, while the remaining one-third of the total samples (99 COAD and 13 normal samples) were used as the validation cohort. The β values of the 17 Hyper-DMPs and 8 Hypo-DMPs in the training cohort were input into WEKA, and the selected attributes function of WEKA was used to filter these candidate biomarkers. As potential diagnostic biomarkers, nine Hyper-DMPs (cg26036626, cg03882585, cg08130988, cg16733654, cg12587766, cg08808128, cg13004587, cg05038216 and cg09746736) and one Hypo-DMP (cg26718707) were selected to construct a COAD-specific diagnostic model (Table 1). Finally, based on the nine Hyper-DMPs and one Hypo-DMP, we constructed a COAD-specific diagnostic model with BayesNet [23].

The average β values of the nine Hyper-DMPs and one Hypo-DMP selected for our diagnostic model in all the COAD tissues, normal tissues and nine types of cancerous tissues from TCGA are visualized in Figure 3A. We performed an unsupervised cluster analysis to evaluate these β values (Figure 3B), and found that they were clearly divided into four clusters. The COAD tumor samples were significantly differentiated from all the normal samples and the tumor samples from the nine other cancer types based on the nine Hyper-DMPs and one Hypo-DMP.

Subsequently, we used the COAD-specific diagnostic model to train the training cohort (including 200 COAD and 25 normal samples from TCGA) in WEKA. Using BayesNet, we determined that the COAD-specific diagnostic model had a sensitivity of 100%, specificity of 99.5% and accuracy of 99.6% in the training cohort (Figure 3C). In the validation cohort (99 COAD and 13 normal samples from TCGA), the diagnostic model had a sensitivity of 100%, specificity of 96% and accuracy of 96.4% (Figure 3D). Therefore, our COAD-specific diagnostic model was confirmed to perfectly distinguish between COAD and normal samples in the training cohort (AUC = 1) (Figure 3E) and the validation cohort (AUC = 1) from TCGA (Figure 3F).

To demonstrate the versatility of our diagnostic model, we conducted a population heterogeneity analysis using the validation cohort from TCGA, which included samples from 4 Asian, 25 black or African American and 73 white patients. The sensitivity, specificity and accuracy are shown in Table 2. Our diagnostic model exhibited no significant population heterogeneity, suggesting that it can be applied to people of different races. In addition, we used the five independent GEO COAD cohorts mentioned above (GSE42752, GSE53051, GSE77718, GSE48684 and GSE77954) as validation cohorts. In receiver operating characteristic (ROC) analyses, the AUCs of these five cohorts were 0.991, 0.964, 0.979, 0.951 and 0.966, respectively (Figure 3G). These results further illustrated the reproducibility and stability of our COAD-specific diagnostic model.

We also analyzed the correlation between the DMP methylation level and the expression of the corresponding genes for the nine Hyper-DMPs and the one Hypo-DMP in our diagnostic model. The Hypo-DMP (cg26718707) corresponded to DIP2C, while the nine Hyper-DMPs corresponded to eight genes: FBLIM1 (cg26036626), SYNE1 (cg03882585), EFEMP1 (cg08130988), PTPRS (cg16733654), LIFR (cg12587766), CLIP4 (cg08808128 and cg05038216), SCGB3A1 (cg13004587) and SLC6A2 (cg09746736). The results of the correlation analysis are shown in Figure 3H and Supplementary Table 1. The expression of DIP2C correlated positively with the methylation level of the Hypo-DMP (r > 0.1, FDR < 0.05), and the expression of the other eight genes correlated negatively with the methylation levels of the corresponding nine Hyper-DMPs (r < -0.1, FDR < 0.05).

Next, using the five independent GEO cohorts, we compared our COAD-specific diagnostic model with three previously reported methylation-based diagnostic models: a Bayesian model including four CpG sites from Azuara et al. [24], a logistic regression model including five CpG sites from Beggs et al. [25] and a logistic regression model including 12 CpG sites from Naumov et al. [26] (Figure 4A). As expected, our model exhibited better sensitivity, specificity and accuracy than the three previously reported diagnostic models in most cases for the five GEO COAD cohorts. We also compared our diagnostic model with these three models in terms of its ability to distinguish COAD from normal tissues and nine other types of cancerous tissues. For this purpose, we divided all the samples from the five GEO COAD cohorts and nine various TCGA cancerous cohorts into a tumor group and a normal group, and we calculated the proportion of samples that were predicted to be COAD among all the samples (Figure 4B, Table 3). In the cohorts from the nine different types of cancers, the ideal proportion would have been 0. When we tested our COAD-specific diagnostic model, almost none of the normal intestinal epithelial samples or the tumor tissues from the nine other cancer types were predicted as COAD (0-5%, median: 0%). However, when we tested the three previously reported diagnostic models, 0-97.7% of the normal tissues (median: 0%, 33.3% and 0%, respectively) and 20.2-98.7% of the tumor tissues from the nine other cancer types (median: 45.4%, 93.5% and 75.2%, respectively) were predicted as COAD. Therefore, our COAD-specific diagnostic model based on nine Hyper-DMPs and one Hypo-DMP not only distinguished COAD from normal samples, but also compensated for the deficiencies of previous COAD diagnostic models that could not differentiate COAD from nine other cancer types.

Then, we used the STRING database to construct a protein-protein interaction network for the nine genes corresponding to the nine Hyper-DMPs and the one Hypo-DMP (Figure 4C). We also performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses on these genes (Figure 4D and 4E). We found that the nine genes were involved in important signaling pathways of tumorigenesis and development, such as Salmonella infection, Janus kinase (JAK)/signal transducer and activator of transcription (STAT) signaling, focal adhesion, proteoglycans in cancer, cytokine-cytokine receptor interactions, etc. All the results of the KEGG pathway analysis and the top 10 results of the GO analysis are shown in Supplementary Tables 2 and 3.

The above results demonstrated that our COAD-specific diagnostic model could accurately and precisely distinguish COAD tissues from normal intestinal epithelial samples and tumor samples from nine cancer types, and that the nine Hyper-DMPs and one Hypo-DMP included in this model may be potential biomarkers for the early prediction and specific diagnosis of COAD.

Identification of prognostic biomarkers of COAD and construction of a combined COAD prognostic model

A total of 287 COAD tissue samples in the cohort from TCGA had both methylated β values and corresponding prognostic information. The distribution and corresponding demographic characteristics of these patients are summarized in Table 4. The patients were divided into a training cohort (N = 143) and a validation cohort (N = 144). The training cohort was used to obtain prognostic biomarkers and construct a COAD prognostic model, while the validation cohort was used to test the COAD prognostic model. Univariate Cox hazard regression analysis of the training cohort revealed 64 CpG sites that correlated significantly with the OS of COAD patients (FDR < 0.05); thus, these CpG sites were identified as candidate prognostic biomarkers. Multivariate Cox stepwise regression analysis was applied to these 64 CpG sites, and six sites (cg00177496, cg01963906, cg05165940, cg12921795, cg19414598 and cg25783173) were included in our final hazard ratio model, which was constructed as a combined COAD prognostic model for OS prediction (Table 5). The six CpG sites from our COAD prognostic model were found to correspond to BDH1 (cg00177496), SYTL1 (cg01963906), SATB2 (cg05165940), WDR20 (cg12921795), DMC1 (cg19414598) and ZNF35 (cg25783173) (Figure 5A and Supplementary Table 4). The expression of SYTL1 correlated positively with the methylation level of cg01963906 (r > 0.1, FDR < 0.05), and the expression of the other five genes correlated negatively with the methylation levels of the remaining CpG sites (r < -0.1, FDR < 0.05).

The risk score formula for our COAD prognostic model was based on the regression coefficients and methylation levels of the six CpG sites, as follows: risk score = (38.52 × cg00177496 β value) – (4.13 × cg01963906 β value) + (2.574 × cg05165940 β value) – (79.32 × cg12921795 β value) + (2.31 × cg19414598 β value) + (6.061 × cg25783173 β value). In the risk score formula, a positive coefficient for a CpG site (cg00177496, cg05165940, cg19414598 and cg25783173) indicates that hypermethylation of that site was associated with shorter OS in COAD patients. In contrast, a negative coefficient for a CpG site (cg01963906 and cg12921795) indicates that greater methylation of that site was associated with longer OS. Our COAD prognostic model revealed that there were significant differences in DNA methylation levels between patients with long-term (> 5 years) and short-term (< 5 years) survival (FDR < 0.05) (Figure 5B). Consistent with the results of the multivariate Cox stepwise regression analysis, the CpG sites with positive coefficients (cg00177496, cg05165940, cg19414598 and cg25783173) exhibited lower methylation levels in patients who survived long-term, while the CpG sites with negative coefficients (cg01963906 and cg12921795) exhibited higher methylation levels in patients who survived long-term. Thus, our combined COAD prognostic model based on six CpG sites successfully distinguished long-term from short-term surviving patients in the training cohort of 143 COAD samples from TCGA.

Based on the Cox regression analyses, risk scores were used as continuous variables in the training and validation cohorts. The risk scores obtained from the combined COAD prognostic model correlated significantly with the OS of COAD patients (Training cohort: likelihood ratio test = 45, p < 0.0001; Wald test = 41.93, p < 0.0001; score [log-rank] test = 50.34, p < 0.0001. Validation cohort: likelihood ratio test = 27.63, p < 0.0001; Wald test = 33.48, p < 0.0001; score [log-rank] test = 39.33, p < 0.0001). Using the median risk score of the training cohort as a cut-off value, we generated Kaplan-Meier curves and performed log-rank tests on the training cohort (Figure 5C) (p = 0.00011) and the validation cohort (Figure 5D) (p = 2e-05). Through these analyses, we sought to compare the OS of patients in the high-risk and low-risk groups and thus determine the predictive value of the combined COAD prognostic model based on six CpG sites. The risk scores for the training and validation cohorts are shown in Supplementary Tables 5 and 6. The survival rate of patients was significantly greater in the low-risk group than in the high-risk group. These results confirmed that our combined prognostic model based on six CpG sites could classify patients into high-risk and low-risk groups, demonstrating its clinical practicability.

To further evaluate the specificity of our combined COAD prognostic model in predicting survival, we used the AUC values obtained from time-dependent ROC analyses as categorical variables. In both the training cohort and the validation cohort, the combined COAD prognostic model precisely predicted the survival of COAD patients, with AUC values of 0.826 and 0.792, respectively (Figure 5E and 5F). We also performed univariate Cox regression analyses of the six individual CpG sites included in the COAD prognostic model (Supplementary Figure 1). The calculated AUC values indicated that the six individual CpG sites could also distinguish high-risk from low-risk patients; however, the predictive effect of any one CpG site was not as good as the predictive effect of the combined prognostic model using all six CpG sites. These results demonstrated that the six CpG sites may be potential prognostic biomarkers of COAD, but the combined COAD prognostic model based on six CpG sites is more valuable than the individual sites for clinical validation and prognostic evaluation.

Independence of the combined COAD prognostic model in OS prediction, and comparison of our prognostic model with other known prognostic models

To evaluate the stability and independence of our combined COAD prognostic model based on six CpG sites, we stratified patients according to clinicopathological characteristics such as age, gender, race, AJCC stage, examined lymph node count and lymphatic invasion. Remarkably, the Kaplan-Meier plots displayed significant extension of OS in the low-risk groups for all these characteristics in the 287 COAD samples from TCGA. Nevertheless, the combined COAD prognostic model predicted the survival of COAD patients more precisely than these factors, with an AUC value of 0.687 (Figure 6A–6C, Figure 7A–7C and Supplementary Figure 2). These results confirmed that the combined COAD prognostic model based on six CpG sites provided an excellent reference for different populations and was an independent predictor of patient survival.

In recent years, DNA methylation biomarkers have been increasingly recognized as important prognostic predictors in COAD. Previously identified markers of COAD prognosis have included hypermethylation of FAM134B [27], higher expression of MMP-11 [28], abnormal methylation and expression of DIRAS1 [29], upregulation of the long non-coding RNA BLACAT1 (a cell cycle regulator) [30] and abnormal expression of 10 microRNAs (including hsa-mir-891a, hsa-mir-6854, etc.) [31]. Dai et al. [32] demonstrated that combined biomarkers of DNA methylation were more sensitive and specific than individual DNA methylation markers. To compare the survival prediction ability of our combined prognostic predictive model with those of previously reported biomarkers, we performed ROC analyses of the previous mRNA, long non-coding RNA and microRNA biomarkers in the validation cohort. Our combined COAD prognostic model based on six CpG sites had a much higher AUC value than the other prognostic biomarkers assayed by ROC analysis in the COAD validation cohort from TCGA (N = 144) (Figure 8A). These results suggested that our combined COAD prognostic model provided more reliable and sensitive predictions of OS than other biomarkers in COAD patients.

Next, we performed a gene set enrichment analysis (GSEA) on the high- and low-risk groups that had been classified according to the median risk score. The pathways that correlated significantly with risk are illustrated in Figure 8B, Supplementary Figure 3 and Supplementary Table 7 (enrichment score [ES] > 0.4, |normalized enrichment score [NES]| > 1, p-value < 0.05 and FDR q-value < 0.25). We selected the intestinal immune network for IgA production (IINIP) for further analysis, since this pathway is known to be involved in COAD. The core enrichment genes of the IINIP pathway were obtained via GSEA (Supplementary Table 8). To determine whether the genes corresponding to the combined COAD prognostic biomarkers functioned in IINIP pathways, we conducted a one-to-one correlation analysis on the expression of the core enrichment genes from the IINIP pathway, the combined methylation level of our prognostic model and the expression of the genes corresponding to the individual CpG sites of the COAD prognostic biomarkers (Figure 8C and Supplementary Figure 4). As expected, the expression of BDH1 and SATB2 and the combined methylation level of our COAD prognostic model based on six CpG sites correlated significantly with the IINIP signaling pathway (p < 0.05). The above results indicated that our combined COAD prognostic model not only accurately predicted prognosis, but also included CpG sites that may directly or indirectly influence the IINIP pathway.

Data source

We downloaded DNA methylation, gene expression and COAD clinical data from TCGA using the University of California Santa Cruz Xena tool (http://xena.ucsc.edu). The DNA methylation data were generated using the Illumina Human Methylation 450 Bead Chip platform, with methylation levels ranging from 0 to 1. We collected the methylation levels of ten types of tumor tissues and normal tissues from TCGA: COAD (299 tumors, 38 normal), BLCA (413 tumors, 21 normal), BRCA (790 tumors, 98 normal), CESC (309 tumors, 3 normal), UCEC (432 tumors, 46 normal), GBM (153 tumors, 2 normal), HNSC (530 tumors, 50 normal), LIHC (379 tumors, 50 normal), LUAD (460 tumors, 32 normal) and LUSC (372 tumors, 43 normal). We calculated the average methylation level of multiple samples as the methylation level of a given CpG site. The level-3 gene expression data were downloaded from RNA-seq HiSeqV2 based on the Illumina HiSeq 2000 RNA sequencing platform, and were obtained as RNA-Seq by Expectation-Maximization (RSEM) [49] normalized counts converted by log2 (x + 1). RSEM software was used to normalize counts on the Xena website.

In addition, we downloaded five DNA methylation array cohorts from GEO (https://www.ncbi.nlm.nih.gov/geo/) as independent validation cohorts. (1) The GSE42752 cohort [26] included 22 COAD samples and 41 normal samples from Russia. Quality control was performed using the R package GenomeStudio (v. 2011.1). Samples were removed if p was > 0.05 and the CpG coverage was < 95%. Thereafter, GenomeStudio was used to normalize the DNA methylation data. The methylation levels of the CpG sites were calculated as β values, where β = intensity (methylated)/intensity (methylated + unmethylated). The data were further normalized using a peak correction algorithm embedded in the IMA R package. Finally, CpG sites on sex chromosomes were removed, and the remaining CpG sites were retained for further analysis. (2) The GSE53051 cohort [50] included 35 COAD samples and 18 normal samples from the US. The R package minfi was used to preprocess the methylation data. To analyze the methylation levels of CpG sites, we averaged the values of all individuals in normal colon and colon cancer. We compared the data from colon cancer and normal control tissues using t-tests. We defined CpG sites with q values < 0.05 and a difference magnitude > 0.1 as important CpG sites. (3) The GSE77718 cohort [51] included 96 COAD samples and 96 normal samples from New Zealand. The CpG sites of each sample were preprocessed and rescaled to ensure that the internal control CpG sites had a common mean across samples. CpG sites located on sex chromosomes or known to cross-react with other regions of the genome were excluded from further analyses. The methylation of the remaining CpG sites was corrected using the COMBAT algorithm to account for batch effects (between-array effects). (4) The GSE48684 cohort [52] included 64 COAD samples and 41 normal samples from the US. CpG sites were removed if p was > 0.05 in the Illumina Infinium DNA methylation data. Thereafter, the R package minfi was used for normalization, including Illumina background level correction, color adjustment and Subset-quantile Within Array Normalization. CpG sites beginning with “rs” on the array were excluded, along with non-CpG sites associated with the X chromosome. The COMBAT algorithm was used to evaluate and correct the batch processing effects of all array runs. (5) The GSE77954 cohort [53] included 20 COAD samples and 11 normal samples from the US. The microarray data were collected at Expression Analysis Inc. (Durham, NC, USA) and preprocessed using the R package methylamine. The array platform was the Human Methylation 450 Bead Chip (GPL13534). We normalized the data from each cohort using the R package limma.

Difference and correlation analysis of DNA methylation and corresponding gene expression

Paired samples are the most suitable for assessing differential methylation levels among individuals [54]. Therefore, we used 25 paired samples to obtain DMPs. CpG sites with > 10% missing values were excluded during the screening process. Missing values in the remaining CpG sites were replaced by the median of the cohort. Then, a paired t-test was used to obtain DMPs between COAD and normal tissues, and the FDR values were adjusted by the Bonferroni method. CpG sites on the sex chromosomes were removed. CpG sites with |Δβ| > 0.2 and FDR values < 0.05 were considered differentially methylated. The R package Pheatmap [55] was used for heat mapping and unsupervised cluster analysis.

When a CpG site corresponded to multiple genes, the optimal corresponding gene was obtained using the R package Champ [56]. We used TCGA data from 25 paired patients with both COAD and normal expression profiles for differential gene expression analysis. The R package limma was used to identify differentially expressed genes from the original data. Genes with a log2 |fold change| > 1 and FDR < 0.05 were considered differentially expressed. Pearson correlation coefficients were calculated to assess the association between DNA methylation and gene expression. Correlations were considered significant based on an |r (cor-value)| > 0.1 and FDR < 0.05. All FDR values were adjusted by the Bonferroni method.

Acquisition of candidate diagnostic biomarkers and construction of a diagnostic model

The 299 COAD and 38 normal samples from TCGA were randomly assigned to the training and validation cohorts at a ratio of 2:1. Five independent GEO cohorts (GSE42752, GSE53051, GSE77718, GSE48684 and GSE77954) were also used as validation cohorts. Firstly, 25 paired COAD and normal samples were compared to obtain DMPs. Secondly, TCGA data from 10 cancer types were used to further screen and narrow the range of candidate DMPs. Candidate Hyper-DMPs were required to have an average methylation level in COAD that was 0.2 units higher than the average methylation level in the normal samples and the nine other types of tumor and normal samples (BLCA, BRCA, CESC, GBM, HNSC, LIHC, LUAD, LUSC and UCEC). Candidate Hypo-DMPs were required to have an average methylation level in COAD that was 0.2 units lower than the average methylation level in the normal samples and the nine other types of tumor and normal samples. Thirdly, we evaluated whether the candidate Hyper-/Hypo-DMPs were significantly associated with the expression of their corresponding genes (|r| > 0.1 and FDR < 0.05), and whether the corresponding genes were differentially expressed genes (log2 |fold change| > 1 and FDR < 0.05). Then, the selected attributes function in the data mining tool WEKA [57] was used to obtain the final list of potential diagnostic biomarkers (nine Hypo-DMPs and one Hypo-DMP). Lastly, for machine learning in WEKA, we used BayesNet to construct the COAD-specific diagnostic model with the nine Hyper-DMPs and one Hypo-DMP in the training cohort from TCGA (including 200 COAD and 25 normal samples).

Then, the diagnostic model was verified with the validation cohort from TCGA (including 99 COAD and 13 normal samples) and five independent GEO COAD cohorts. Firstly, we imported the methylation values of the nine Hyper-DMPs and one Hypo-DMP from the training cohort using the Filter option of the WEKA Preprocess panel. Then, we selected BayesNet in the Classifier option of the Classify panel to build our model. Thereafter, we imported the methylation values of the nine Hyper-DMPs and one Hypo-DMP from the validation cohort (TCGA validation set or GSE42752 or GSE53051 or GSE77718 or GSE48684 or GSE77954) to the Supplied test set option of the Classify panel. Finally, selecting BayesNet in the Classifier option, we identified the sensitivity, specificity, and accuracy of our model in the evaluated cohort.

Acquisition of candidate prognostic biomarkers and construction of a prognostic model

A total of 287 COAD tissue samples in the cohort from TCGA had both methylated β values and corresponding prognostic information. These samples were randomly assigned to the training cohort and the test cohort at a ratio of 1:1. Firstly, univariate Cox proportional hazard regression analysis was performed in the training cohort from TCGA (143 COAD samples) to identify significant methylation markers associated with the OS of COAD patients (p < 0.05). Then, multivariate Cox stepwise regression analysis was performed on these CpG sites, and sites with p-values > 0.05 were removed from the feature cohort in each iteration. The R packages Survival and Mass were jointly used to complete the multivariate Cox stepwise regression analysis.

The Cox proportional risk model was used to determine patients’ hazard ratios and corresponding 95% confidence intervals. The linear combination of model predictors weighted by their regression coefficients was used as the formula to predict the survival risk of patients. The high-risk and low-risk groups were classified according to the median risk value. The R packages Survival and Plot were then used to plot Kaplan-Meier survival curves to visualize the cumulative survival of the patients at risk at some time point. A log-rank test was used to evaluate the difference in OS between the high- and low-risk groups. Finally, the area under the ROC curve was determined by ROC analysis using the R package SurvivalROC. The effectiveness of the risk score in predicting OS was also evaluated. All statistical calculations were performed using the R statistical environment (R version 3.5.4).

STRING database

The genes corresponding to the diagnostic model were analyzed using the STRING functional protein-protein interaction network (9.1) [58]. The same website was used to analyze the input GO biological processes and KEGG pathways. P-values < 0.05 were considered significant.

GSEA

After the prognostic prediction model was used to calculate patients’ risk scores, GSEA (JAVA version) [59] (http://software.broadinstitute.org/gsea/index.jsp) was performed for the high-risk and low-risk groups. GSEA includes four key statistics: the ES, NES, FDR q-value and p-value. An ES > 0.4, |NES| > 1, p-value < 0.05 and FDR q-value < 0.25 were used to filter the GSEA results. Based on these statistics, all the genes in the list of specific genes were scored and ranked.