All clinical information and microarray data in the two molecular levels were downloaded from TCGA Data Portal (https://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp). A total of 153 colon Adenocarcinoma cancer samples with gene expression microarray data and DNA methylation microarray data had subtype labels from previous study (3). The platforms of gene expression and DNA methylation microarray are Custom Agilent 244K Gene Expression Microarray (AMDID019760) and Illumina Infinium DNA methylation (HumanMethylation27 BeadChip), respectively. The data level 3 in the data portal was used in this study, which means that the gene expression data was Lowess normalized and the ratio of the Cy5 channel and Cy3 channel were log2 transformed to create gene expression values for 23199 probe-sets, resulting in 17814 genes available for further analysis, and the DNA methylation data contain beta-value calculations, HUGO gene symbol, chromosome number and genomic coordinate for each targeted CpG site on the array. Approximately 27578 CpG sites were located in proximity to the transcription start sites of 14475 consensus coding sequences.

We combined gene expression data of 153 samples into one file, and imputed the missing value using KNN Imputed (13). The informative genes for clustering analysis were selected using a threshold standard deviation SD>1 across all samples, and this resulted in 1393 genes. To perform consensus clustering (14) we used K-mean approach with average linkage to detect robust clusters, where the metric was 1 minus the Pearson's correlation coefficient. The procedure was run over 2000 iterations and with a sub-sampling ratio of 0.8. To evaluate the heterogeneity of the subtypes we applied silhouette width values to identify the most ‘core’ members of each subtype (15–17), and samples with Silhouette Score>0.5 were considered as core samples. Significance analysis of microarrays (18) (SAM) was applied to identify differentially expressed genes between subgroups, and the Prediction analysis of microarrays (19) (PAM) was used to obtain marker genes and establish classifiers. The training set for PAM is 70% of 153 samples selected randomly and the testing set is the other 30% of the samples. The Gene Ontology analysis was performed for each subtype using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (20,21), and GeneMANIA (22) was applied to find the co-expressed network of marker genes.

After combining data into one file, we removed the probes containing any ‘NA’ marked data points and the probes that were designed for the sequences on the X and Y chromosomes. We then conducted a filtering process to reach a final data matrix with 1491 probes, which exhibited sufficient variable methylation levels with a threshold standard deviation value (SD>0.2) across all samples. The DNA methylation microarray data were β-value, following β-distribution. To use the consensus clustering method, a data set must be transformed so that it follows a normal distribution. We used the Transfer Function (23,24) to transform the β-value into M-value which is normally distributed, which was similar with RPMM (25) used in β-value in Hinoue et al (11). Since some subtyping systems were reported in early studies on DNA methylation of colon cancer, we only performed the PAM on all samples and did not build testing sets. DAVID and GeneMANIA were also used on DNA methylation data.

All data analyses were done in R platform (Windows version 2.15.2) (26,27). For the categorical variables in clinical information table such as gender, tumor subtype (previous studies), oncogene mutation (Yes or No), the Fisher's exact test was used to assess the significance of their association to the subtype derived in this study. For age levels, we used ANOVA to assess differences among subtypes. The package ConsensusClusterPlus was used to perform unsupervised clustering analysis. Package SAMr and PAMr were applied to identify the differentially expressed genes, to build the classifier and to determine the marker genes, respectively.

Clinical and pathologic features of the patients and their tumors were summarized for further analysis. All 153 patients had information on age, gender, AJCC stage, vital status, tumor location and subtypes from earlier studies (Table I).

Unsupervised K-mean consensus clustering was used to uncover potential subgroups of colon cancer on the basis of the similarities of their gene expression values of 1393 informative genes. We let K=2 to 6 in core K-mean clustering, two subgroups could be identified when K=2 and the cluster consensus are 0.98 and 0.99 for each subgroup (Fig. 1A and D), thus the first subgroup was named as ECL1 with 104 samples (68%) and the second subgroup was named as ECL2 with 49 samples (32%). When K=2 to 4, the ECL2 subgroup showed steady and consistency (Fig. 1A–C). The relationship between two subgroups and their clinical characteristics were listed in Table II.

In ECL2, the age of onset (73.3±11.47) is significantly higher than ECL1 (P<0.049, ANOVA). We found that the majority samples of ELC2 are right sided tumors. All the MSI-H samples were found in the ECL2, and all the CIN samples in the ECL1. Furthermore, these two subgroups showed no significant difference in AJCC stage and history of polyps. Mutations of KRAS, BRAF and TP53 were investigated in many studies, we found that all samples with BRAF mutation were in ECL2 and most of samples with TP53 mutation were in ECL1 (Table III).

Nearly 62% of the samples in ECL1 were left sided tumors. Most of ECL1 samples were MSS status and the majority samples of Invasive subtype (3) were in ECL1. Compared with those reported in previous studies, we found that ECL1 contained both CIN and Invasive subtypes (3), and therefore we examined the heterogeneity of this subgroup (Fig. 2A).

We carried out unsupervised clustering analysis only on ECL1 samples with K=2 to 6. When K=3, we discovered three distinct subclasses with very clear boundaries (Fig. 3B). There are 87 samples with Silhouette Score >0.5 considered as core samples and retained, with 29 samples in subclass 1, 30 samples in subclass 2 and 28 samples in subclass 3. There are 18 CIN samples in subclass 1, 27 CIN samples in subclass 2, and 23 Invasive samples in subclass 3, and Fig. 2B demonstrates the relationship between subtypes reported earlier (3) and the subclasses derived from ECL1 (P<1.065e-10, Fisher's exact test).

There were two subclasses correlated to the CIN subtype, and due to the heterogeneity of ECL1 subgroup we investigated the difference of these two CIN groups. CIN samples extracted from two subclasses were compared using SAM with Wilcoxon rank sum test. There were 250 differentially expressed genes found with 2-fold change, and only 6 genes were upregulated in subclass 2, namely, SLC25A21, POPDC3, GREG2, HOTAIR, GYPB and SLC35F4. The proportion of either the metastatic samples or the death samples in subclass 2 was roughly two-fold of that in subclass 1 (Table IV).

On the top level we identified two subgroups in colon cancer, ECL1 had relatively high heterogeneity and it was associated with CIN and Invasive subtype derived from earlier studies, whereas ECL2 showed high homogeneity. On the secondary level, three subclasses were derived from ECL1, where the subclass 1 and 2 were associated with CIN subtype and the subclass 3 was associated with Invasive subtype.

PAM analysis was carried out to identify marker genes that could discriminate the two subgroups on the top level. When Δ =4.16 (overall error rate 0.019 at minimum), 256 genes were selected from the 107 training samples. The testing set was used for independent validation, and only 2 samples were classified into wrong groups with an overall error rate of 0.043.

There were 137 genes out of the 256 marker genes that were upregulated in ECL2, among them SPP1 and POSTN were associated with metastasis and poor prognosis in colorectal cancer, which were reported in earlier studies. DAVID analysis showed that these 137 genes were enriched in immune response, defense response, response to wounding, inflammatory response and carbohydrate binding GO terms. Furthermore, the GSEA (28) analysis of these genes showed that they were upregulated in advanced gastric cancer and basal subtype of breast cancer. There were 119 genes upregulated in ECL1, and they were enriched in ERBB receptor signaling network, and β-oxidation of pristanoyl-CoA pathways. Finally, we plotted a heating map with the 256 marker genes for all 153 cancer samples (Fig. 4) with sample resorted hierarchical clustering and only 5 samples were classified into incorrect groups. This suggested that these genes could serve as feature genes for the subtype classification.

To investigate the subtypes using DNA methylation data, we applied the same method to the transformed methylation array data. When K=3 or 4 (Fig. 5A and B), the clustering reached the highest consensus. When K=3, we named these subgroups as MCL1 with 57 samples (37%), MCL2 with 40 samples (26%) and MCL3 with 56 samples (37%). We found that the gender proportion among the three subgroups showed significant difference (P<0.029, Fisher's exact test). The age distribution among three subgroups also showed significant difference (P<2.24e-3, ANOVA, Table V).

Majority of the samples in MCL1 were left tumors (~79%), of the minimum mean age, MSS status and no BRAF mutation. More than 50% of the samples in MCL1 had TP53 mutation and a few samples had KRAS mutation. Almost all samples in MCL2 were male, right tumors (~93%), of the maximum mean age and more than 50% of the samples were MSI-H status; all samples with BRAF mutation were in MCL2 and a few samples in this subgroup had KRAS mutation and TP53 mutation. More than 50% of the samples in MCL3 were female, right tumors, MSS status, and there were no BRAF mutation and nearly 50% of the samples had KRAS mutation and TP53 mutation (Table VI).

Compared with the results of TCGA and Hinoue et al (3,11) (Fig. 6A), most samples in cluster 4 fell into MCL1, and all of CIMP-H samples were in MCL2; majority of the samples in CIMP-L and cluster 3 fell into MCL3 (P<2.2e-16, Fisher's exact test).

Characteristics of MCL3 were quite similar with CIMP2. The CIMP2 showed more heterogeneity than the other two (12), hence we further examined the subdivision of MCL3. When K=4, the four subgroups generated were largely overlapped with the previous classification (Fig. 6B), but the cluster consensus were lower than that when K=3 (Fig. 5D). To judge whether the CIMP-L and cluster 3 were distinct subtypes of colon cancer, we examined the data in Table I of Hinoue et al (11), and we found that tumor location and the frequence of TP53 mutation exhibited significant difference between the two clusters.

PAM analysis was applied in these three subgroups to identify DNA methylation gene maker panels which could discriminate the subgroups. Firstly, MCL2 (CIMP-H) was compared with the combination of the MCL1 and MCL3 (Non-CIMP-H), when Δ=11.4, 52 probes corresponding to 47 genes were selected as the first panel, and the overall error rate was 0.052. DSC3, LOX, RUNX3, SLC30A2 and TLR2 harbored two hypermethylation sites in the samples from MCL2 subgroup. Secondly, regardless of MCL2, MCL1 (cluster 4) was compared with MCL3, when Δ =6.99 and overall error rate was 0.079, 39 probes corresponding to 33 genes were selected as the second panel. ELMO1, JAKMIP1, NCAM1, NDRG4 harbored two hypermethylation sites in the samples from MCL3.

Combining two marker panels, there were 80 methylation genes. DAVID analysis on these genes showed that they were enriched in cell fate commitment, neuron differentiation, extracellular matrix, and sequence-specific DNA binding GO terms. We also used GeneMANIA to build the co-expression network of these 80 genes, and it turned out that the Wnt receptor signaling pathway and the digestive system development pathway were involved in the network.

We performed hierarchical clustering on all 153 samples with the genes in the two panels and were able to find three subtypes in DNA methylation data. Labels of ECL1 and ECL2 in each sample were also listed. Almost all samples in MCL2 were overlapped with those in ECL2; moreover, the ECL1 comprised MCL1 and MCL3 (Fig. 7).