Development and validation of a risk prediction model and nomogram for colon adenocarcinoma based on methylation-driven genes

Identification of aberrantly methylated and differentially expressed genes in COAD

An analysis flow chart of our bioinformatics workflow is shown in Figure 1A. A total of 1940 differentially expressed genes (DEGs) were detected by overlapping date from the TCGA database and GSE39582 (Figure 1B). Similarly, 6681 differentially methylated genes (DMGs) were identified by overlapping TCGA data and GSE48684 (Figure 1C). Subsequently, we overlapped these DEGs and DMGs, and identified 659 aberrantly methylated DEGs, including 129 genes with high expression and hypermethylation, 188 genes with low expression and hypermethylation, 192 genes with high expression and hypomethylation, and 150 genes with low expression and hypomethylation (Figure 1D).

Figure 1. Identification of differentially expressed genes (DEGs) and differentially methylated genes (DMGs) in colon adenocarcinoma (COAD). (A) Flowchart showing overall design and analytic procedure of this study. (B) Overlapping of DEGs from TCGA and GSE39582. (C) Overlapping of DMGs from TCGA and GSE48684. (D) Identification of aberrantly methylated DEGs.

Identification of MDGs in COAD

Promoter hypermethylation can trigger transcriptional silencing of cancer-related genes. Therefore, we selected genes with high methylation and low expression for further analysis. We evaluated the Pearson coefficients from gene expression and methylation values for aberrantly methylated DEGs. In total, 97 aberrantly methylated DEGs were identified as MDGs (Pearson coefficient < -0.3 and P < 0.05; Figure 2 and Supplementary Table 1).

Figure 2. Correlation between the expression value and methylation value of the methylation-driven genes (MDGs) in colon adenocarcinoma (COAD) tissues.

Development of a risk prediction model of COAD patients

There were 414 COAD patients with both expression data and complete clinical information in the TCGA database, thus, we used these datasets to identify prognostic genes for COAD. Univariate Cox regression analysis initially identified that among 97 MDGs, 6 MDGs (CBLN2, GSTM1, RBM47, SH3GL3, SLCO4C1, and TMEM220) were significantly correlated with OS of COAD patients (Table 1, P < 0.05). GSTM1 and SH3GL3 were excluded due to having a Hazard ratio (HR) >1. These four prognostic genes were then utilized to build a best-fit risk prediction model using least absolute shrinkage and selection operator (LASSO) Cox regression analysis. The risk prediction formula was as follows: Risk score = (-0.121 * Expression level of CBLN2) + (-0.377 * Expression level of RBM47) + (-0.065* Expression level of SLCO4C1) + (-0.136 * Expression level of TMEM220). We then calculated the risk scores of 414 COAD patients using the formula above. The distribution of risk scores and the patients’ survival status are shown in Figure 3A. A risk heatmap was used to visualize the expression profiles of these four prognostic genes (Figure 3A). The median risk (-5.979) was used as a cutoff point to divide COAD patients into a high-risk group (n = 207) and a low-risk group (n = 207). Kaplan-Meier (K-M) analysis showed the patients in the high-risk group had worse prognosis than those in the low-risk group. (Figure 3B, P = 0.004). The areas under the curves (AUCs) of the 1-, 2-, and 5-year OS predictions were 0.669, 0.651 and 0.652, respectively (Figure 3C). Meanwhile, compared with any single mRNA, the signature from all four genes had higher accuracy for predicting a patients' OS (Supplementary Figure 1). These results showed that this genetic signature was effective for OS prediction.

Figure 3. Validation and development of risk prediction model in colon adenocarcinoma (COAD) patients. (A) Risk score distribution in COAD patients, survival status of COAD patients, and expression heatmap of four methylation-driven genes (MDGs) in a Cancer Genome Atlas (TCGA) training cohort. (B) The K-M curve of overall survival (OS) for COAD patients between two different groups in our TCGA training cohort. (C) Time-dependent ROC curves at 1 year, 2 years, and 5 years in the TCGA training cohort. (D) Risk score distribution of COAD patients, survival status of COAD patients, and expression heatmap of four MDGs in a validation cohort. (E) The K-M curve of OS for COAD patients between two different groups in a validation cohort. (F) Time-dependent ROC curves at 1 year, 2 years, and 5 years in a validation cohort.

In order to clarify the importance of the four prognostic genes above in COAD patients, we used the GSE17536 array data as an independent validation set. We calculated the risk scores of all patients according to the risk score formula mentioned above. As expected, scatter plots, risk heatmaps and K-M curves were able to accurately distinguish high- and low-risk group patients using this dataset (Figure 3D, 3E). The AUCs of the 1-, 2-, and 5-year OS predictions were 0.576, 0.544, and 0.574 using the GSE17536 array data, respectively (Figure 3F). These results indicated that our model was effective for OS prediction in COAD patients.

Stratification analysis

A stratification analysis was then performed to investigate the applicable population using this risk prediction model and TCGA data. Patients were assigned to different subgroups based on their age (≤ 60 / > 60), sex (female/male), T stage (1 + 2/ 3 + 4) and TNM stage (I + II / III + IV). The results showed that our risk prediction model could delineate high- and low-risk patients in each age, sex, T staging, and TNM staging subgroup. Unfortunately, the risk prediction model lost its prognostic value in patients of grade T (1 + 2) and TNM (I + II), which we attributed to the small sample size of these subgroups in the TCGA COAD dataset (Figure 4).

Figure 4. K-M analysis of overall survival (OS) for patients stratified by age, gender, T stage, and TNM stage.

Establishment and evaluation of a predictive nomogram

Based on our stepwise and multivariate Cox regression analyses, we found that age, TNM stage and risk score were independent prognostic factors in 414 COAD patients (Figure 5A). Considering the different clinical characteristics of each patient, we built a nomogram that combined age, TNM stage, and risk score to individually predict the 1-, 2-, and 5-year OS of COAD patients (Figure 5B). The AUCs of the 1-, 2-, and 5-year of this nomogram were 0.776, 0.761 and 0.740, respectively (Figure 5C). The AUCs of 1-, 2-, and 5-year prognosis from traditional TNM stage were 0.728, 0.710, and 0.672, respectively (Figure 5C). At the same time, the concordance index (C-index) of the nomogram was significantly higher than traditional TNM stage (0.755 versus 0.706, P < 0.05). Therefore, in terms of predicting the OS of COAD patients, our nomogram was better than traditional TNM staging. Based on the median of the nomogram score as a cutoff value, patients were then divided into high-risk and low-risk groups. K-M analysis revealed that the high-risk group had significantly poorer OS (Figure 5D, P < 0.001). The calibration curves of our nomogram suggested that the predicted OS was consistent with the observed OS (Figure 5E).

Figure 5. Establishment of an overall survival (OS) nomogram for colon adenocarcinoma (COAD) patients. (A) Univariate and multivariate analyses of risk score and clinical variables. Red solid dots represent significant difference, and black solid dots mean no difference. (B) A nomogram individually predicting OS in COAD patients. (C) The time-dependent ROC of our nomogram and TNM stage in the prediction of prognosis at 1-, 2-, and 5-year time points. (D) The K-M curve of our nomogram. (E) Calibration plot of the nomogram. The predicted and the actual probabilities of OS are plotted using blue solid and black dotted lines, respectively.

In the validation phase, a new nomogram still showed a higher predictive efficacy in using GSE17536 array. Similar to its performance in the TCGA cohort, the AUCs of the 1-, 2-, and 5-year nomograms were greater than those from TNM stage, respectively (Figure 6A). The calibration curves of the nomograms from 1-, 2- and 5-year OS displayed obvious concordance between the predicted OS and the observed OS, respectively (Figure 6B). In addition, the C-index values of the nomogram and TNM stage were 0.778 and 0.774, respectively. Meanwhile, K-M curves could still distinguish high - and low-risk group patients (Figure 6C).

Figure 6. Validation of nomogram in a validation cohort. (A) Shown is the time-dependent ROC curves for 1-, 2-, and 5-year overall survival (OS) predictions from our nomogram compared with TNM stage. (B) Calibration curve for our nomogram in a validation cohort. The predicted and the actual probabilities of OS are plotted using blue solid and black dotted lines, respectively. (C) OS of our nomogram in a validation cohort.

Encyclopedia of genes and genomes (KEGG) enrichment of four candidate genes

We performed gene set enrichment analysis (GSEA) with our four candidate genes to investigate the potential biological mechanisms via these genes in COAD progression. Patients were divided into high-expression and low-expression groups based on the median expression value of these candidate genes. The results showed that the four candidate genes were involved in multiple tumor-associated pathways, such as the apoptosis, the calcium signaling pathway, the colorectal cancer, the Hedgehog signaling pathway, the JAK-STAT signaling pathway, and the TGF-β signaling pathway (Figure 7, Supplementary Table 2–5).

Figure 7. Representative enriched pathways in four candidate genes from gene set enrichment analysis (GSEA) software.

Validation of differential expression of CBLN2 and TMEM220 due to promoter methylation

First, to further verify our results based on data from the TCGA and GEO databases, we used quantitative real-time PCR (qPCR) to determine candidate gene expression levels in NCM460 cells and SW480 cells, respectively. We found that the expression of CBLN2 and TMEM220 was low in SW480 cells, but very high in NCM460 cells (Figure 8A). Second, in order to determine whether abnormally methylated promoter regions directly caused transcriptional silencing of CBLN2 and TMEM220, SW480 cells were treated with the DNA methyltransferase inhibitor 5-Aza-2′-deoxycytidine (5-aza), and the expression of CBLN2 and TMEM220 was determined via qPCR. This study found that the expression of CBLN2 and TMEM220 was restored in SW480 cells after treatment with 5-aza (Figure 8B). Third, methylation-specific PCR (MSP) was applied to identify the methylation status of the CBLN2 and TMEM220 promoter regions. Studies have previously shown that these regions are partially methylated in SW480 and SW620 cells (Figure 8C). CpG islands situated in the CBLN2 and TMEM220 promoter regions and the designed MSP primers are shown in Figure 8D. In summary, we confirmed that the expression of CBLN2 and TMEM220 was silenced by the methylation of these promoter regions in a COAD cell line.

Figure 8. Experimental verification in colon cells and tissues. (A) qPCR was performed to identify the relative expression of CBLN2 and TMEM220 in NCM460 and SW480 cells. (B) qPCR was carried out to assess CBLN2 and TMEM220 expression levels in SW480 cells before and after treatment with 5-aza. (C) Methylation status of CBLN2 and TMEM220 was determined by MSP in SW480 and SW620 cells. (D) Schematic diagrams of CpG islands in the promoter regions of CBLN2 and TMEM220. (E) Representative images of immunohistochemistry staining of colon sections from colon adenocarcinoma (COAD) tissues (n = 46). Original magnification, ×100. (F) Prognostic significance of CBLN2 and TMEM220 expression in COAD patients. (G) Knockdown of RBM47 gene in sw480 cells.

Relationship between CBLN2 and TMEM220 expression and OS of COAD patients

The expression of CBLN2 and TMEM220 in 46 COAD tissues was then examined by immunohistochemistry. CBLN2 and TMEM220 protein expression levels were significantly different in tumor tissues as compared to controls (Figure 8E). A marker was considered positive when 20% or more cells were stained [17]. The prognostic effects of CBLN2 and TMEM220 on the OS of COAD patients were next evaluated through K-M analysis. As shown in Figure 8F, both CBLN2 and TMEM220 led to significant survival differences in 46 COAD samples (P = 0.033 and P=0.047, respectively). The clinical information for these patients is listed in Table 2.

Effects of RBM47 gene knock-out on CBLN2, SLCO4C1 or TMEM220 expression levels

Considering that RBM47 had the greatest influence on our risk prediction model, we further explored the correlation between RBM47 and CBLN2, SLCO4C1 and TMEM220. The expression values of CBLN2 and TMEM220 were decreased and the expression value of SLCO4C1 was increased in SW480 cells when RBM47 was knocked out (Figure 8G).

Materials acquisition and preprocessing

DNA methylation data, transcriptome data and corresponding clinical data about COAD tissues were obtained from the TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases in April, 2019. Gene methylation data from the TCGA dataset was generated using the Illumina Infinium HumanMethylation450 microarray, which included 310 COAD and 37 adjacent non-tumor samples. If any gene had multiple cg sites, the empty sites were removed and the mean value of β was used to represent its methylation level [38]. Gene transcriptome data from the TCGA database (Level 3) was normalized and log2 scaled using the functions DEGList and calcNormFactors in the edgeR package for R [39], which included 473 COAD and 41 adjacent non-tumor samples. Gene methylation data from the GSE48684 arrays was generated using the GPL13534 platform (Illumina HumanMethylation450 BeadChip). Gene transcriptome data from the GSE39582 and GSE17536 arrays was generated using the GPL570 platform (Affymetrix Human Genome U133 plus 2.0 Array). The GSE48684 array consisted of 106 COAD and 41 adjacent non-tumor samples. The GSE39582 array consisted of 566 COAD and 19 adjacent non-tumor samples. The GSE17536 array consisted of 177 COAD samples. We also retrospectively collected 46 cases of COAD tissues from patients who underwent surgical resection in the Fourth Hospital of Hebei Medical University, China (from December 2010 to December 2013). All patients had resectable COAD, and none of them had received preoperative anticancer treatments. They were followed until June 2018. Ethical permission of this study protocol was granted by the ethical committee of Hebei Medical University. All patients were informed and signed informed consent forms prior to enrollment in the study.

Identification of aberrantly methylated DEGs

For TCGA transcriptome data, the EdgeR package was used to identify the DEGs between COAD and non-tumor samples, and an absolute value of the log2 fold change (|log₂FC|) >1 and false discovery rate (FDR) < 0.05 were considered statistically significant. For GEO transcriptome data, the limma package was used to identify DEGs between COAD and non-tumor samples, with the thresholds of FDR < 0.01 and |log₂FC| > 0.5. All methylation data was analyzed with the limma package. Herein, genes with FDR < 0.05 were considered as DMGs. Finally, aberrantly methylated DEGs were detected by overlapping DEGs and DMGs in Venny software 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/).

Correlation analysis of aberrantly methylated DEGs

To study the transcriptional regulation of DNA methylation, we evaluated the Pearson coefficient between gene expression and the methylation data for aberrantly methylated DEGs. A total of 322 COAD samples with matching methylation data and expression data were used for correlation analysis. Aberrantly methylated DEGs with a Pearson coefficient < -0.3 and P < 0.05 were defined as MDGs [40]. Scatter plot of these MDGs was plotted using ggplot2 in R.

Development of a risk prediction model

Initially, univariate Cox regression analysis was used to evaluate the association between MDGs and the OS of COAD patients, and MDGs with a P < 0.05 were selected for further analysis. Based on the expression value of MDGs, the LASSO Cox regression models were used to develop a best-fit risk prediction model with the R package “glmnet”. The risk score for each COAD patient was calculated as follows:

$$\text{Risk}\text{score}={\displaystyle \sum _{i=1}^n{\exp }_{\text{i}}{}^{\ast }\beta \text{i},}$$

where n is the number of prognostic genes, exp_i is the expression value of each gene i, and βi is the weighted regression coefficient in gene i from multivariate Cox regression analysis. Then, time-dependent ROC curve and K-M analyses were used to evaluate the predictive ability of our model. In the validation phase, we verified the risk prediction model using the GSE17536 dataset, another COAD cohort.

Construction and assessment of nomograms

Stepwise and multivariate Cox proportional hazard regression models were used to distinguish independent prognostic parameters of COAD patients, based on which we developed a nomogram. K-M analysis, time-dependent ROC curve analysis, calibration plot and C-indices were used to evaluate the discriminative ability of our nomogram. A C-index was calculated to assess nomogram discrimination by means of the bootstrap method with 1000 resamples. We assessed the performance of our nomogram on predicting OS for COAD patients using traditional TNM stage as a control. Meanwhile, the GSE17536 cohort was used to verify the new nomogram.

Pathway enrichment analysis of MDGs

GSEA was carried out to explore the underlying biological mechanisms of each marker. GSEA software was downloaded from the GSEA home (http://software.broadinstitute.org/gsea/index.jsp). “c2.cp.kegg.v7.2.symbols.gmt gene sets” was used as a reference gene set to enrich KEGG pathways for candidate genes in TCGA. Last, |NES| > 1 and P < 0.05 were set as thresholds.

Validation experiments in colonic cells

To verify the transcriptional silencing relationship of prognostic genes, we used three human colon cell lines (NCM460, SW480, and SW620) for validation. NCM460 and SW480 cells were cultured in RPMI 1640 (Gibco, Carlsbad, CA, USA). SW620 cells were cultured in DMEM (Gibco, Shanghai, China). All cell culture medium was supplemented with 10% fetal bovine serum (Invitrogen, Carlsbad, CA, USA) and 1% penicillin/streptomycin. To investigate the effect of 5-aza (Sigma, St. Louis, MO, USA) treatment, SW480 cells were treated with 3 μM for 72 h [41]. Meanwhile, SW480 cells were transfected with control and RBM47 siRNAs (Thermo Fisher, USA) according to the manufacturer's protocol. Cells were siRNAs were treated for 48 h and then switched to media lacking siRNA. Total RNA was then isolated from cells utilizing the Trizol method (Invitrogen, Shanghai, China). qPCR was performed on an ABI 7500 real-time PCR System (Applied Biosystems, Carlsbad, CA) using SYBR Green (Takara, Japan). GAPDH was used as an internal reference, and the relative expression level of each gene of interest was calculated with the formula 2^-ΔΔCt [42]. The methylation status of MDGs was tested in SW480 and SW620 cells by MSP. We predicted CpG islands and designed MSP primers with Methyl Primer Express software v1.0 (Thermo Fisher Scientific, Waltham, MA) based on the genomic sequence around the transcriptional start site (TSS) of each gene. qPCR and MSP primers are illustrated in Supplementary Table 6.

Immunohistochemistry

Paraffin-embedded specimens from colon tissues were sectioned to a 5 μm thickness. The sections were then deparaffinized in xylene and rehydrated through graded alcohol solutions. Antigen extraction was performed using citrate buffer (pH 6.0), and sections were stored in Tris buffered saline (TBS). Endogenous peroxidase activity was blocked by incubation in 3% hydrogen peroxide. The sections were incubated with the anti-CBLN2 antibody (1:100, Abcam, Shanghai, China) or anti-TMEM220 antibody (1:100, Abcam, Shanghai, China) overnight at 4° C. The reaction products were visualized with diaminobenzidine (Vector labs, Burlingame, CA, USA) as the chromogen and counterstained with hematoxylin. Finally, images were acquired with immunofluorescence microscopy.

Statistical analysis

All statistical analysis was performed in R 3.5.0 and GraphPad Prism 6.0 (GraphPad Software, La Jolla, CA, USA). A two-sample t-test was used to compare gene expression levels in colon cell lines. P < 0.05 was considered statistically significant.