Signature consist of nine differentially methylated gene pairs to predict the prognosis for stage II colorectal cancer patients

Abstract

Colorectal cancer (CRC) is the third most commonly malignance cancer worldwide. The quality of life and life expectancy of CRC survivors are seriously impacted by the recurrence after resection surgery. Therefore, there is a need to develop a reliable signature to predict the recurrence after resection surgery in CRC patients. Epigenetic alterations, especially gene methylation, have been found highly related to the metastasis of tumors. However, limited evidence is available for CRC patients. Moreover, there remains a knowledge gap of a comprehensive view of gene methylation in stage II CRC. In this study, we conducted and validated a differentially methylated 9-gene-pair prognosis signature to predict recurrence after resection surgery in stage II CRC patients. In addition, we use Univariable and Multivariable Cox regression to estimate the association between relapse free survival times in relation to each candidate gene pairs and the 9-gene-pair score. Moreover, the network analysis found 8 hub genes, including NDRG4 and BMP3, out of all candidate genes have large degree of interactions in the protein-protein interaction network. Seven hub genes have been found associated with CRC progression. Furthermore, the function enrichment analysis suggested that differentially methylated genes were mostly enriched in the cell cycle pathway. In summary, a 9-gene-pair signature was developed and validated in this study. The signatures identified in this study have the potential to be applied in the prognosis of post-surgery recurrence in the stage II CRC patents.

INTRODUCTION

Colorectal cancer (CRC) is one of the leading causes of cancer death worldwide. According to a recently published report, it is estimated that 693,900 deaths are caused by CRC annually worldwide [1]. The incidence and mortality rates of CRC have significantly decreased, due to improved screening, treatment and modification of environmental factors. However, the quality of life and life expectancy of CRC survivors are seriously impact by the recurrence after resection surgery. Approximately, 20% to 30% patients going through resection surgery will suffer from recurrence [2]. This recurrence leads to a necessarily of chemotherapy to prolong their survival time, especially in stage II CRC patient [3, 4]. Considering the reduced life quality caused by chemotherapy and the poor prognosis of recurrence of CRC, there is an urgent need of developing a reliable signature which can predict the prognosis of stage II CRC survivor after resection surgery.

Epigenetic alterations are clearly related to the cancer occurrence and development [5]. CRC is also a type of cancer which arises through the accumulation of epigenetic alterations [5, 6]. DNA methylation is the major type of epigenetic alteration that closely related to CRC [7]. In CRC, the abnormal DNA methylation can occur in around 75% to 80% of the cases [5]. The most recent review has summarized the findings concerning epigenetic biomarkers of CRC, especially differentially methylated gene-specific DNA [8]. Given that evidence is accumulating that there exist a crosstalk among genes, essentials have been pointed out to provide a comprehensive view of differential methylated genes [9, 10]. To our best knowledge, the prognostic signatures based on gene-specific DNA methylation for CRC patients are still lacking. Moreover, gene expression files are prone to systematic measurement biases, which normally caused by batch effects in experiment procedures [11-13]. Thus, a signature, which takes into consideration data normalization and batch effects, is highly required to be generally used to independent datasets [13]. Furthermore, it has been emphasized that patients in different stage of CRC should be treated differently [2]. This brings up a concern that whether it is possible a signature can be generating in patients in different stages. To illustrate this point, patients in stage II are under the highest risk of recurrence after resection surgery compare to other stages. It is understandable that stage 0 and stage I patients are at the early stage of the disease, and are easier to control. Patients under stage III are generally have adjuvant chemotherapy along with the resection surgery, which also decreased their risk of recurrence [2]. Therefore, it has utmost importance to establish a prognosis signature for stage II CRC patients.

To fill up the current knowledge gap, we conducted and validated a prognostic gene-pair based signature for stage II CRC patients to predicate the risk of relapse after resection surgery. We further performed the functional enrichment analysis to provide insights into the mechanism underlying the association between differentially methylated gene-specific DNA and prognosis of stage II CRC patients.

RESULTS

Identification of differential methylated gene and prognosis related gene

In total, information on differential methylated gene-specific DNA information was available for cancer tissues from 458 CRC patients and para-cancerous normal tissues from 75 patients. Altogether, 11,287 significantly hyper-methylated CpG sites and 990 significantly hypo-methylated CpG sites were detected by using SAM. After the exclusion of genes with inconsistent methylation aberrations, 7,873 significantly hyper-methylated genes and 616 hypo-methylated genes were identified. By using univariable Cox regression models, 460 hyper-methylated genes and 7 hypo-methylated genes were identified to be associated with the RFS, and together defined as DM genes. The general flowchart for the analysis was shown in Figure 1.

Figure 1: Flowchart of construction and validation of the prognostic signature. The flowchart showed three major analysis steps: the development (step 1) and validation (step 2) of the gene pairs signature, and the function analysis for the CRC prognostic genes (setp 3).

Identification of prognostic gene pair signature

We paired the 467 DM genes into 108,811 gene pairs, from which 4,617 gene pairs were significantly associated with RFS (FDR<0.01). Using C-index and forward selection procedure, a 9-gene-pair signature was identified as prognosis related to RFS for stage II CRC patients (Table 1). The C-index of this 9-gene-pair signature is 0.816. We regard a patient at high risk of relapse when more than 5 out of the 9 gene pairs significantly verse methylated.

Table 1: The 9-gene-pair prognostic signature ^a

Signature	RMOs (Ma > Mb)		C-Index	Hazard Ratio	P-value
	Gene a	Gene b
gene pair 1	LAMA1	TMEM43	0.6953	5.0815	0.0001
gene pair 2	NUDT21	PHKG2	0.6847	3.4707	0.0011
gene pair 3	LAMA1	CCL28	0.6751	4.1125	0.0001
gene pair 4	NVL	PHKG2	0.6681	3.3206	0.0011
gene pair 5	CHST12	EIF4G3	0.6567	3.2913	0.0012
gene pair 6	TOB2	FAM13B	0.6478	3.1399	0.0019
gene pair 7	PPARD	TSPO	0.6456	3.9223	0.0006
gene pair 8	ARL14	VNN1	0.5982	2.9013	0.0088
gene pair 9	LYPD6	HCFC1R1	0.5976	7.9021	0.0000

^a Ma and Mb represent the methylation value of gene a and b. RMOs (Ma > Mb) represents the relative methylation ordering of gene pair (a, b). Hazard Ratio the parameters were calculated by the Cox regression model. Where Hazard Ratio > 1 indicates that Ma > Mb is a risk factor.

Survival analysis

In total, 145 stage II CRC patients with RFS information were available and comprised the test sample. The clinical information for this test sample was summarized in Table 2. We first performed a univariable Cox regression analysis for the association between RFS and each DM genes (data not shown). We then categorized these patients into high-risk and low-risk group based on each identified DM gene-pair. Figure 2 showed the univariable Cox regression analysis of the association between RFS and different risk groups. We further categorized patients into high-risk group and low-risk group according to the 9-gene-pair signature. Univariate and multivariate Cox regression analysis showed that the RFS was significantly different between high-risk and low-risk group (P-value all < 0.0001). Adjustment for age or gender did not change this observation (Figure 3).

Figure 2: The Kaplan-Meier curves of relapse free survival of prognostic groups predicated by the each gene-pair in the 145 stage II CRC patients. Kaplan-Meier curves of relapse free survival for gene pair 1 (A); gene pair 2 (B); gene pair 3 (C); gene pair 4 (D); gene pair 5 (E); gene pair 6 (F); gene pair 7 (G); gene pair 8 (H); gene pair 9 (I). The sample was classified into the high-risk group (red line) if and only if at least 5 of 9 gene pairs voted for high-risk; otherwise, the low-risk group (blue line). RFS, relapse-free survival.

Figure 3: The Kaplan-Meier curves of relapse free survival of prognostic groups predicated by the 9-gene-pair in the 145 stage II CRC patients. Kaplan-Meier curves conducted from (A) multivariable cox regression analysis and (B) univariable cox regression analysis of prognostic groups predicated by the 9-gene-pair in the 145 stage II CRC patients. (C) is additionally adjusted for gender, and (D) is additionally adjusted for age. The sample was classified into the high-risk group (red line) if and only if at least 5 of 9 gene pairs voted for high-risk; otherwise, the low-risk group (blue line). HR, hazard ratio.

Table 2: Characteristics of 145 stage II CRC patients

Covariates	Category	Total
Age, years	<60	34
	>=60	111
Gender	Male	79
	Female	66
RFS status	replase	30
	non-replase	115

Network analysis of prognosis related 9-gene-pair signature

We identified 381 DM genes between high-risk and low-risk stage II CRC patients, which were categorized by the 9-gene-pair signature. Based on these 381 DM genes, we conducted a one-step PPI network (Figure 4). Eight hub genes (ATM, DLX5, ESF1, AKT1, BMP3, NDRG4, CHEK2, ABL1) with the largest degrees (all ≥ 8) were found in this PPI network. We tested the ability of NDRG4 and BMP3 methylation as marker for detection of colon cancer by testing its ability to detect aberrant vimentin exon 1 methylation in fecal DNA samples. Fecal DNAs from 22 and 25 of these 32 patients tested positive for NDRG4 and BMP3 methylation, respectively, in this assay, yielding a 68% and 78% clinical sensitivity for detecting the presence of a colon cancer.

Figure 4: PPI network constructed using the DMG genes between the high and low risk stage II CRC patients. A grey line indicates the interaction between two proteins that are linked by the grey line. Red round nodes, hub-genes with high amount of interactions. PPI, protein-protein interactions

Function analysis of the prognosis relate DM genes

We performed Spearman correlation analysis between 16 candidate genes of the 9-gene-pair signature and all 8,489 DM genes. This resulted in 5,532 high correlated DM genes, with FDR <0.01 and r > 0.7. KEGG pathway enrichment analysis showed that all these DM genes significantly enriched in the cell cycle pathways (Figure 5).

Figure 5: KEGG pathways analysis of 5532 significant differentially methylated genes.

DISCUSSION

In this study, we developed a set of prognostic signatures based on RMOs of 9-gene-pair for stage II CRC patients. We found 381 DM genes which are significantly differentially expressed between high-risk and low-risk groups categorized by the 9-gene-pair signature. Among these DM genes, 8 hub genes were observed to have large degrees in the PPI networking. Enrichment analysis suggested these DM genes, which differentially methylated between high-risk and low-risk stage II patients, were significantly enriched in the cell cycle pathways.

The identification of differentially methylated gene pair signatures related to prognosis of stage II CRC patients would benefit the personalized treatment. Here, we applied RMOs procedure in the construction of the prognosis methylation signature. The prognostic signature base on RMOs is robust against experimental batch effects and data normalization, and can be easily applied at the individual level to sample which are measured in different laboratories [13]. We performed Cox regression variable analyses to further investigate the association between RFS with candidate gene pairs and the 9-gene-pair signature. The univariable Cox regression analysis indicated all 9 candidate DM gene pairs are significantly associated with the prognosis of CRC. This further supported our observed association between the 9-gene-pair signature and RFS, which suggested stage II patients with altered signature have a higher risk of relapsing. Besides, the prognosis ability of this score is better than any individual candidate gene pair. This further emphasized the importance of developing a comprehensive signature for cancers [13]. We also performed multivariable Cox regression to estimate the association between RFS and the signature. The results from Cox regression suggested that this signature is independent from both gender and age in the prediction of the prognosis of the stage II patients. This is vital in the practice usage of this signature, which can be generalized into a wider patients group independent from gender or age.

The PPI network analysis of the candidate and DM genes, which were highly related to prognosis of CRC patients, provided insights into the underlying mechanisms of CRC progression. We found 8 hub genes which have high degree of interaction in the PPI network, namely ATM, DLX5, E2F1, AKT1, BMP3, NDRG4, CHEK2 and ABL1. Except for BMP3, the other 7 hub genes have been found related to neoplastic transplantation. Previous studies have found that ATM and AKT1 are associated with the progress of multiple types of cancers [14-16]. For instance, Sun et al. found that ATM can involve in the tumor metastasis thorough regulation of Snail in breast cancer [17]. Rachelle et al. found that AKT1 promotes mammary tumor induction but not metastasis [16]. In addition, ATM is regard as a key mediator in the DNA damage response, which plays a vital role in the stages of tumor development and also related to the effectiveness of standard cancer treatment. [18]. Moreover, ATM and AKT1, as tumor suppressors, has been found frequently notate in colorectal cancer, this might explain their relationship with relapsing in post-surgery patients [19]. Similarly to the finding of ATM, CHEK2 also participates in the DNA damage response progress, and further related to the tumor metastasis [20]. An existing study found the DLX5 showed DNA methylation in >50% CRC tissue, which is in line with our observation [21]. DLX5 also has been found important as a cancer-hypermethylated gene with demonstrated tumor-suppressor functions. Expression level of NDRG4 also has been found strongly associated with the prediction and prognosis of colorectal cancer. Evidence suggested E2F1, BMP3 and ABL1 are associated with the progression of CRC, but mostly induced by non-coding RNAs [22-24]. These findings further supported the reliability of the prognosis ability of our conducted 9-gene-pair signature. To further explore the underlying mechanism, we did the functional enrichment analysis of the candidate genes and genes which are highly correlated with these candidate genes. Our observed enrichment in the cell cycle is in line with previous studies, which suggested the cell cycle pathway is a determinant in the efficiency of CRC therapy [25, 26]. The limitation of this study was the lack of the detailed information about the detail distance between tumor tissues and adjacent normal tissues.

In conclusion, we developed a 9-gene-pair differentially methylated prognostic signature to predict the prognosis of stage II CRC patients. Furthermore, the network analysis and functional enrichment analysis provided insights into the underlying mechanism of the progression of the stage II CRC patients.

MATERIALS AND METHODS

Data sources and data preprocessing

Information of methylated gene-specific DNA for colorectal cancer was collated from The Cancer Genome Atlas (TCGA) (downloaded on 2017-06-09). A total of 458 CRC patients and 75 para-cancerous normal tissues from patients were available. All gene-specific DNA methylation data was measured by Illumina Human Methylation 450 Beadchip (450K array). Stage II CRC patients with RFS information available were extracted to form a test sample. In this study, we included 25,978 CPG sites at the promoter regions of genes for analyzing.

Data pre-processing

For gene-specific DNA methylation profiles, a CPG site was excluded when it is missing in more than 50% of the total samples (n=533). We performed K Nearest Neighbors procedure to predicate the value for other CPG site when it is missing in less than half of the samples.

Identification of differential methylated gene-specific DNA

We carried out Significance Analysis of Microarrays (SAM) [27, 28] to identify the differential methylated (DM) CpG sites in both cancer tissues and para-cancerous normal tissues. To control the false discovery rate (FDR) [29], we adjusted the P-values with the Benjamin and Hochberg (BH) procedure in the multiple binomial tests. The preliminary identified DM CpG sites were annotated to genes according to the annotation table of 450K platform file. We excluded the genes with inconsistent methylation aberrations (hyper-methylation or hypo-methylation) states within their promoter regions from the analysis.

Survival analysis of the differentially methylated gene-specific DNA

The 50^th percentile of the abnormal methylation among all identified differentially methylated gene-specific DNA was defined as a cut-off point. Stage II CRC patients with RFS information, i.e. the test sample, were classified into a high-risk group or low-risk group according this cut-off point. Univariable Cox regression analysis [30] was performed for each identified differentially methylated gene-specific DNA in relation to the RFS of the patients.

Identification of prognostic gene-pair signature

From the all identified differentially methylated genes, we picked out gene A and gene B and paired them into gene pairs. The methylation levels of gene A and gene B were indicated by Ma and Mb respectively. Patients from the test sample were classified into two subgroups according to the relative methylation orderings (RMOs) (Ma > Mb or Ma < Mb) of the gene pair in each sample. Univariate Cox regression was used to evaluate whether there is a significant difference in RFS between patients in each subgroups. A prognosis related gene pair was defined as gene pair with 1% FDR-value control. Afterwards, C-index of each gene pair was evaluated. Gene pairs with highest C-index were picked up as the seed prognosis related gene pairs. These seed gene pairs were further added to a signature one at a time until the C-index of the signature was not improved by any additional gene pairs. The gene pairs included in this signature were defined as candidate gene pairs. Then, we performed a forward selection procedure to search the optimal subset of the prognosis related gene pairs that reached the highest C-index. To be specific, a patient was assigned to the high-risk group when at least 50% of the gene pairs voting for high risk, otherwise patient were assigned to the low-risk group. The optimal subset of gene pairs with the highest C-index was identified as the final prognostic gene pair signature.

Survival analysis for gene pairs and the prognosis gene-pair signature

The univariable Cox progression model [31] was used to evaluate the association between RFS with DM gene pairs or the prognosis gene-pair signature. Survival curves were estimated using the Kaplan-Meier method and were compared using the log-rank test [32]. Multivariable Cox regression models were further used to evaluate the age-independent and gender-independent prognostic value of the signature.

Testing of expression from Fecal DNA

Stools were collected from a total of normal individuals with no prior history of colon cancers or polyps and from a population (n = 32) of colorectal cancer patients. Samples were treated as previously described [33]. Target human vimentin DNA fragments were purified from total DNA preparations by acrylamide gel–based affinity capture as previously described [33]. Total captured DNA from each sample was then applied to bisulfite-modification and MS-PCR, and the results were analyzed in a manner blinded to patient’s disease status.

Protein-protein interaction network

To recognize potential prognosis related gene candidates, we conducted a protein-protein interaction network (PPI network). Cytoscape was used in this study to construct a PPI network [34]. Information on the interactions between proteins was collected from SIGnaling Network Open Resource (SIGNOR) database [35].

Functional enrichment analysis

Enrichment analysis for KEGG Pathway analysis was performed for gene categories which were downloaded from KEGG (downloaded on 2017/05). This analysis aimed to identify the biological functions that were most significantly to genes that differentially methylated in the stage II CRC patients. To limit the chance finding, the p-values were adjusted by the BH procedure.

CONFLICTS OF INTEREST

All authors declare that he/she has no conflicts of interest.

FUNDING

This study was supported by Science and Technology Planning Project of Guangdong Province, China (NO. 2014A020212045); the Science and Technology Project of Shenzhen, China (NO. JCYJ20150330102720178); Shenzhen Health and Family Planning Commission Scientific Research Project (NO. 201604130300); Clinical Doctor-Basic Scientist Combination Foundation of Shenzhen Second People’s Hospital (NO. 2015024).