Transcription factor expression repertoire basis for epigenetic and transcriptional subtypes of colorectal cancers

Aberrant cancer-specific CGI promoter DNA methylation alterations are an early event during tumorigenesis, occurring as early as in the aberrant crypt foci precursor lesions of CRCs and polyps (6, 11, 36, 37). The relationships of methylation changes at promoter and enhancer elements during CRC development are not known. We established the relationships of methylation changes at promoter and enhancer elements during CRC development (Materials and Methods) by combining gene expression with genome-wide DNA methylation at the comprehensive set of promoter and enhancer candidate cis-regulatory elements (cCREs). The cCREs represent regulatory DNA elements which are defined by integrating DNA accessibility (DNase I-hypersensitive sites) and chromatin modification data in the ENCODE project and represent a robust set of well-annotated regulatory elements across a wide array of tissues and cell types (38). The cCREs are classified into promoter-like signatures, proximal enhancer-like signatures that are within 2 kb of a TSS, and distal enhancer-like signatures that are more than 2 kb away from the nearest TSS (38) (Fig. 1A). In our analyses, we further parsed these regions by the presence or absence of CGI.

Combined analyses of methylation patterns in CRCs with that of tubulovillous adenomas and serrated adenomas (detailed in Dataset S1) (1, 14) show that most of the PLS-CGI DNA methylation alterations observed in CRCs are conserved in the polyps (Fig. 1B). For the most variable PLS-CGI probes in the TCGA dataset (COAD), colon polyps harbor the range of promoter DNA methylation patterns similar to that of CRCs and segregate well with the different CIMP classes (Fig. 1B). Further, as in previous reports (11, 39), the tubulovillous and serrated polyps clearly separate out within this spectrum of CIMP categories with the latter clustering with the CIMP-H and CIMP-L CRCs. The dysplastic serrated polyps appear very similar to that of the nondysplastic serrated polyps (SI Appendix, Table S1). These data recapitulate that genome-wide CIMP features at CGI promoters are established very early during CRC development and support the model that serrated and conventional tubular/tubulovillous adenomas represent precursors that lead to CRCs having distinct epigenetic alterations along the CIMP spectrum (11, 39).

We show that the cancer-specific gain of DNA methylation significantly extends to the pELS and dELS enhancer regulatory regions in a manner quantitatively similar to the spectrum of CIMP phenotypes. To relate the CIMP subtypes, which is based on promoter-CGI methylation, to the degree of enhancer methylation, we have introduced the term CpG-enhancer methylator phenotype or CEMP classification based on the methylation changes in the pELS and dELS regions. The enhancer methylation is similar in CRCs and polyps and occurs in both CGI (Fig. 1 C and D), and non-CGI regions (SI Appendix, Fig. S1 A and B). Further, this spectrum of methylation patterns is also observed for an independent set of putative enhancer regions identified as non-promoter-associated accessible chromatin using ATAC-seq data for a subset of colon cancers in the TCGA (40) (SI Appendix, Fig. S1C). Hierarchical clustering of methylation at these putative CRC-specific enhancer elements mirrors the same spectrum of CIMP-associated hypermethylation patterns observed in the CRCs and adenomas. Importantly, enhancer regions overall have a cancer-specific gradient in hypermethylator phenotypes (CEMP) classically shown for promoter CGI regions (CIMP). This is observed in the high correspondence of the CRC samples assigned to the epigenetic subtypes (CIMP-H, L, L3, and L4) based on methylation at the CGI promoter regions (classical promoter CIMP subtypes) and those based on methylation at the enhancer regions (pELS and dELS; CEMP-H, L, L3, and L4) (Fig. 1 E and F).

Previous studies have shown that promoter-CGIs associated with the PcG mark (H3K27me3) in progenitor tissue stem cells and/or embryonic stem cells have a high likelihood of being hypermethylated in cancers (41–44). We show a similar likelihood for enhancer hypermethylation to be associated with PcG mark in polyps and CRCs by analyzing the degree of hyper and hypomethylation in enhancer regulatory elements categorized into various activity states based on H3K27me3, H3K4me1, and H3K27ac marks (ChromHMM patterns from ENCODE) (45). In both CRCs, and polyps, the degree of hypermethylation at promoter and enhancer elements harboring the various chromatin marks follow the gradient of hypermethylator phenotypes, with CIMP-H showing the highest fraction of probes methylated in each category of regulatory elements (Fig. 2 A and B). The degree of hypermethylation varies by the chromatin marks and presence of CGI, with the PcG-associated promoter elements in the context of CGI showing the highest proportion of hypermethylated probes: For example, in the CRC samples, about 52% of probes in PcG-TSS-CGI are hypermethylated in CIMP-H CRCs, and this drops to ~35% for PcG-TSS-non-CGI regions, i.e., PcG containing promoter elements not associated with CGI (SI Appendix, Fig. S2A). In contrast and as expected, active promoters exhibit a lower degree of hypermethylation in both CGI and non-CGI (the highest observed being ~7 and 8%, respectively, for the CIMP-H CRCs). Like observed for the PcG-marked promoters, PcG-marked enhancers, harboring the bivalent enhancer mark (H3K27me3 + H3K4me1), have a high frequency of hypermethylated probes in the context of CGI (~60%), which drops to ~40% in the non-CGI regions (SI Appendix, Fig. S2A). Although the proportions of methylated probes significantly decrease for the non-PcG marked poised and active enhancers, the fraction of probes hypermethylated in these regulatory elements follow the same increasing trend from low to high hypermethylator phenotypes, indicating that tumor-specific hypermethylation along the CIMP spectrum is a general phenomenon at these gene regulatory elements (Fig. 2A and SI Appendix, Fig. S2A). In relation to H3K27me3-marked regions, H3K9me3-marked constitutive heterochromatin showed weaker relation to DNA hypermethylation along the CIMP spectrum. Only H3K9me3-marked regions in the context of CGI showed an increasing trend to be methylated with increasing CIMP phenotype, but this was far less compared to that observed for H3K27me3-marked regions, while the H3K9me3-marked regions associated with non-CGI have a low tendency to be methylated across the CIMP spectrum. Similar hypermethylation patterns and association with the chromatin marks are observed in the adenoma samples for these ChromHMM promoter and enhancer regulatory element classes (Fig. 2 A and B and SI Appendix, Fig. S2A). Further, the genome-wide hypermethylation patterns at the different regulatory elements are tightly linked as the fraction of probes methylated in each of these regulatory element categories are highly correlated in CRC and adenoma samples (SI Appendix, Fig. S2 C and D). Taken together, the above analyses using multiple independent enhancer annotations (Figs. 1 and 2) show that methylation patterns in the enhancer elements across the CRCs and the polyp datasets have the same relationship with hypermethylator subtypes classically observed for the promoter regions and that these linked promoter (CIMP) and enhancer (CEMP) hypermethylator patterns occur early during colon tumorigenesis.

Similar analyses as above, but for the degree of hypomethylation across the CIMP subtypes, show that in contrast to the increases in hypermethylation along the spectrum of CIMP phenotypes in CRCs and polyps (Fig. 2 A and B), hypomethylation is relatively similar across all CIMP subtypes (Fig. 2 C and D). Genomic regions harboring active marks, especially in the context of non-CGI, tend to show a higher degree of hypomethylation compared to inactive marks or those associated with CGI. The degree of hypomethylation at the different regulatory elements is less correlated to each other in both the CRCs and adenoma samples (SI Appendix, Fig. S2 E and F) in contrast to the highly correlated hypermethylation occurring genome-wide at these regulatory elements (SI Appendix, Fig. S2 C and D). The similarities in hyper- and hypomethylation changes in polyps and CRCs indicate that both methylation gains and losses are established early during tumor development.

The high correlation between the CIMP and CEMP phenotypes indicates that the genome-wide hypermethylation alterations at gene-regulatory elements constitute systemic epigenetic dysregulation phenomenon during early stages of tumorigenesis which potentially has common mechanistic bases. We identified that among genes whose expression is positively or negatively correlated with the hypermethylator phenotypes, there is a significant enrichment for the class of DNA binding TFs (229 TFs enriched in 1813 genes that are negatively correlated with hypermethylator phenotype) (SI Appendix, Fig. S3 A and B). As reduction in TF binding to target sites has been directly related to de novo methylation alterations (46, 47), we examined the role of TF expression changes in the hypermethylator phenotypes. TFs with regulatory roles in the hypermethylated and hypomethylated promoter and enhancers cCREs were identified in the TCGA colon cancer dataset using the ELMER approach, which scans for TFs with binding sites in these regions and for which the expression of corresponding TFs is negatively correlated with the cancer-specific differential methylation alterations in these regions (see workflow in Fig. 1A, model in Fig. 3A, and details in Materials and Methods) (33). These TF sets associated with hypermethylated or hypomethylated promoter (PLS) and enhancer (pELS and dELS) elements significantly overlap with each other (Fig. 3B, TFs listed in Dataset S2). As expected, overall expression of TFs associated with hypermethylated cCRE sites is down-regulated, while the TFs associated with hypomethylated cCRE sites are up-regulated, in the CRC samples compared to normal samples (Fig. 3C). The degree of downregulation of TFs with binding sites enriched in hypermethylated cCREs is significantly higher (~threefold) than the degree of upregulation of TFs with binding sites enriched in hypomethylated cCRE, as seen in the higher negative median log gene expression ratio (−1.54) for the down-regulated TFs compared to that for the up-regulated TFs (+0.80) in relation to normal samples (Fig. 3C). Importantly, these factors are also deregulated to a similar extent in the polyp dataset, indicating that the TFs and their expression trends identified in the CRC dataset represent early events associated with neoplastic transformation. Further, the hypermethylated cCRE–associated TFs that are down-regulated are primarily enriched for TFs relevant to cancer-related pathways, while the hypomethylated cCRE–associated TFs that are up-regulated are enriched in pathways regulating cell cycle, viral response pathways, among other pathways (Fig. 3D).

Given that TFs have a direct role in modulating methylation (46, 47), deregulated expression of these cCRE methylation–associated TFs may potentially influence the DNA methylation landscape at regulatory regions and the manifestation of the hypermethylator phenotypes. In support of this, the expression of the cCRE methylation–associated TFs occurs in a manner correlated with the quantitative CIMP phenotypes. We identified the TFs associated with hyper- and hypomethylated cCRE (PLS and dELS) for the two extremes of the hypermethylator phenotypes, highest (CIMP-H) and lowest (CIMP-L4), to determine the relation between TF expression changes and the CIMP phenotypes (Fig. 4 A and E; TFs listed in Dataset S2; data for PLS and dELS shown in the following figures; similar results obtained for pELS due to the high overlap of the TFs among these regulatory elements). While the TFs associated with hypermethylated PLS or dELS regions are expected to be down-regulated, the following key characteristics associated with the degree of their downregulation and correlation with the hypermethylator phenotypes suggest that TF expression profiles strongly influence the cancer-specific global hypermethylator patterns. Although TFs commonly associated with hypermethylated PLS (Fig. 4B) and dELS (Fig. 4F), regulatory elements in both the CIMP-H and L4 subtypes are down-regulated in both these subtypes, CIMP-H cancers tend to display higher levels of downregulation. The latter is exemplified in the median expression ratio of the CIMP-H- and CIMP-L4-specific cCRE hypermethylation–associated TFs in the CIMP-H and L4 subtypes with respect to normal samples (Fig. 4 L, Top). Noticeably, the CIMP-H-specific cCRE hypermethylation–associated TFs are markedly down-regulated in the CIMP-H cancers (Fig. 4 C and L), which is in stark contrast to the degree of downregulation of the CIMP-L4-specific cCRE hypermethylation–associated TFs in CIMP-L4 cancers (Fig. 4 D and L). The cCRE hypomethylation–associated TFs identified in CIMP-H or CIMP-L4 are up-regulated to similar levels across the CIMP subtypes (Fig. 4 B–D and F–H and SI Appendix, Fig. S4A), which is consistent with findings that the degree of hypomethylation is similar across the CIMP subtypes (shown in Fig. 2 C and D and SI Appendix, Fig. S2 B, E, and F). Importantly, the cCRE hypermethylation–associated TFs specific to CIMP-H, in contrast to those specific to CIMP-L4, display a robust trend for decreasing expression with increasing degree of the hypermethylator phenotype from CIMP-L4 to H (SI Appendix, Fig. S4 A, Left).

We independently validated the above relationship of the augmented TF downregulation in the CIMP-H subtype by linking TF expression changes to the methylation frequency (measured as the relative risk for methylation) at TSS distal sites (>5 kb from TSS) occurring near the TF binding motifs of corresponding up- and down-regulated TFs in each tumor sample. In these analyses, the relative risk for hypermethylation at sites near TF motifs for up- and down-regulated TFs in each CRC sample was estimated without a priori requiring negative correlation of TF expression and target site DNA methylation. Relative risk for hypermethylation is higher at sites near motifs of down-regulated TFs, as opposed to up-regulated TFs, across all the CIMP categories, and these relationships are particularly highest for the CIMP-H subtype (SI Appendix, Fig. S4B). Further, the associated TFs are markedly down-regulated as a function of the increasing hypermethylator phenotype, with most downregulation in the CIMP-H subtype (SI Appendix, Fig. S4C).

As previous studies have shown that expression of epigenetic modifiers, such as UHRF1, is directly linked to the CIMP profiles (48), we compared the impact of expression of genes encoding key enzymes involved in DNA methylation modification (DNMT1, DNMT3A, DNMT3B, TET1, TET2, TET3, and UHRF1), in relation to the median expression of the TFs associated with the hypermethylated cCRE elements, on the variability in hypermethylation at promoters across the CRC CIMP spectrum. DNMT1 and UHRF1 both are highly expressed in the CIMP-H samples, which is in direct relation to their potential increased activity being associated with increased hypermethylation in the CIMP-H CRC and adenoma samples (SI Appendix, Fig. S4 D and E). However, the standardized regression coefficients associated with the expression of hypermethylation-associated TFs as well as these epigenetic modifiers computed using multiple linear regression model show that TF expression has the strongest and most significant association with the degree of hypermethylation across the CIMP spectrum (SI Appendix, Fig. S4F).

Furthermore, as immune cell infiltration varies based on the CRC subtypes, especially with the MSI samples enriched in the CIMP-H subtype showing favorable immune cell infiltration (49, 50), we tested whether the CRC tumor purity differences among the CIMP subtypes due to differences in the immune microenvironment may be directly related to the cCRE methylation–associated TF expression differences in the CIMP subtypes. We applied the CIBERSORT algorithm to identify the proportions of immune cell infiltration in the TCGA CRC cohort, which shows an increased immune activity in the CIMP-H subtype in agreement with previous reports (49, 50) (SI Appendix, Fig. S4G). Multiple regression analyses of the variability in hypermethylation at promoters across the CIMP categories as a function of the proportion of various immune cell type infiltration along with the median expression of TFs associated with the hypermethylated cCRE elements show that TF expression has the highest and significant standardized regression coefficient explaining the variability in hypermethylation across the CIMP spectrum (SI Appendix, Fig. S4H).

As CIMP subtypes have been reported in cancers of various tissues of origin, we queried whether downregulation of TF expression is associated with increasing methylator phenotypes in other cancers. We analyzed TCGA samples for which sufficient numbers of matched methylation and expression data are available and where CIMP categories have been well described. Accordingly, gastric, esophageal, and glioma cancer datasets (TCGA) were identified to satisfy these criteria and were analyzed to identify TFs significantly associated with hypermethylated and hypomethylated dELS regulatory elements. For esophagus and glioma cancer, the CIMP-H and L (intermediate methylation) subtypes were combined into a high methylation group (CIMP-H’), and subtypes L3 and L4 were combined into a low methylation group (CIMP-L’) due to low sample number availability. The TFs commonly associated with the hypermethylated cCREs in both the CIMP-H/H’ and L4\L’ subtypes are down-regulated to a similar extent in gastric, esophageal, and glioma cancers (Fig. 4 I–K). Like the observations in the CRCs, most significant downregulation is consistently observed for TFs associated with hypermethylated cCREs in the CIMP-H subtype of these cancers (CIMP-H and L panels in Fig. 4 I–K and Fig. 4 L, Top). These data thus suggest that the increased downregulation of subsets of TFs is associated with the CIMP-High phenotypes across cancer types.

Analyses of gene expression relationships of putative enhancer-target genes reveal that the enhancer methylation alterations are functionally linked to expression of their associated target genes early during CRC development. Enhancer gene targets were identified in the CRC data using the ELMER algorithm, which identifies target genes based on their proximity to enhancer probes and significant cancer-specific negative correlation in expression of the gene and enhancer methylation. To evaluate the functional relevance and early roles, these enhancer methylation and target gene expression relationships were then explored in the independent polyp dataset. As expected, based on the ELMER approach for identification of enhancer–gene pairs, the target genes of the hypermethylated dELS enhancer regions are down-regulated in CRCs compared to normal samples (Fig. 5A). Importantly, a similar trend for decreased expression of these target genes is also observed in the independent polyp dataset. In the case of hypomethylated dELS regions, the overall expression of target genes in CRCs is higher than that in normal samples, and again, this relationship is mirrored in the polyps (Fig. 5B). The negative correlations between dELS methylation and expression of its target gene, identified in the CRC samples, are also present in polyps for substantial numbers of enhancer–target gene pairs, with 49% (367/750) and 75% (544/725) of the hyper- and hypomethylated enhancer probe–gene pairs showing significant negative correlation (P < 0.05) in adenomas (Fig. 5 C and D). Similarly, for the hypermethylated CGI promoters in CRCs whose direct target genes have lower expression in CRCs, these genes are down-regulated in polyps (SI Appendix, Fig. S5A) with a significant negative correlation between promoter methylation and expression in the polyps (SI Appendix, Fig. S5B). Examination of individual genes highlights the relationship between enhancer methylation changes and target gene expression in CRCs and polyps (SI Appendix, Fig. S5 C–F). These example genes include those important in cancer pathways, like SOX17 (51, 52) and LMX1A (53, 54) that are down-regulated in conjunction with hypermethylation of their cognate enhancer probes, and MYC (55, 56) and SLC6A6 (57) that are up-regulated in conjunction with hypomethylation (SI Appendix, Fig. S5 C–F).

Consistent with the functional nature of individual genes described above, the genes affected by enhancer methylation are enriched for biological pathways relevant to tumorigenesis. Our previous studies have shown that promoter-methylated genes in CRCs are important regulators of tumor suppression and pathways important in cancers, development, and differentiation (44, 58). This is again the case for the target genes of hypermethylated enhancers which are enriched for pathways such as cancer development and Wnt signaling (Fig. 5 E, Bottom). In contrast, the hypomethylated enhancer target genes are involved in basic cell functions such as cell cycle, RNA transport, and DNA replication, important for increased cell turnover (Fig. 5 E, Top).

As introduced earlier, recent transcriptome-based profiling of CRCs into CMS subtypes has revealed clinically significant subgroups linked to important molecular pathways and the CIMP subtypes (2). We thus analyzed the relation between expression of the cCRE methylation–associated TFs with transcriptome-based CMS subtypes to determine whether expression of these TFs may form a basis for the CIMP–CMS relationships. To determine the CIMP–CMS relationship, we applied TF expression–based clustering on CRC datasets (TCGA and GSE39582 cohort, which were used in previous studies for developing the CMS subtyping) (2). Clustering CRCs based on expression of CIMP-specific cCRE-associated TFs shows that these TFs not only distinguish the CIMP-H subtype but also separate the CMS categories into their respective subgroups (Fig. 6A). Using the K-means clustering approach, we identified three major CRC clusters (K1, K2, and K3) that showed significant overall association with both CIMP and CMS subtypes (SI Appendix, Fig. S6F, Fisher’s test, P < 0.001). In this combined cohort, CIMP status is designated as CIMP-H, CIMP-Low, and CIMP-Negative (CIMP-Neg, constituting CIMP-L3 and L4), based on methods used to identify the CIMP subtypes in the non-TCGA cohort (2). The majority (~77%) of CIMP-H samples separated in the distinct K1 cluster, which is enriched for CMS1 and three subtypes (Fig. 6B), indicating that expression of CIMP-specific cCRE methylation–associated TFs may provide a strong basis linking the CIMP-H and CMS1 and three subtypes. The CIMP-L and negative subtypes are mainly enriched in K2 and K3 clusters, which are discretely enriched for CMS2 and CMS4 subtypes, respectively. In order to test the enrichment of each CIMP/CMS subtype in the three K-means clusters, we have performed Fisher’s tests of pairs CIMP-H vs. non-CIMP-H, CIMP-Low vs. non-CIMP-Low, CIMP-Neg vs. non-CIMP-Neg, (CMS1+CMS3) vs. (CMS2+CMS4), CMS2 vs. (CMS1+CMS3+CMS4), and CMS4 vs. (CMS1+CMS2+CMS3) against the three K-means clusters (SI Appendix, Fig. S6F). Except for CIMP-Low vs. nonCIMP-Low comparison, every CIMP and CMS subtype is significantly enriched in specific K-means clusters. Taken together, expression of cCRE methylation–associated TFs not only separates the CIMP-H group but also sufficiently recapitulates the CMS subtypes. Thus, expression of these TFs, which potentially gives rise to the methylation spectrum defining the CIMP-H subtype, may provide a strong basis linking the CIMP and CMS subtypes.

CIMP-H CRCs have been variably associated with survival outcomes (17, 19, 21, 59, 60). A major indicator of the survival outcomes in the context of CIMP-H is the MSI-status, with the CIMP-H/MSS phenotype being an unfavorable indicator of overall survival (OS) while MSI-H being a favorable indicator (18, 21, 61–63). We show that stratification of CIMP-H samples based on expression of CIMP-H-specific hypermethylated cCRE-associated TFs allows better discrimination of patients by OS than that using MSI-based stratification (Fig. 6 C and D). Stratification of CIMP-H CRCs based on expression of TFs associated with hypomethylated regions specific to both CIMP-H as well as CIMP-L subtypes showed no significant differences in OS (SI Appendix, Fig. S6 C and D). In these analyses, the CIMP-H cases in general did not show overall differences in survival compared to CIMP-L and negative (SI Appendix, Fig. S6A). Stratifying the CIMP-H cases by MSI status (MSI-high and MSS) shows only marginal trends for better OS in MSI-high cases (P-value = 0.22, HR = 0.65) (SI Appendix, Fig. S6B). In contrast, CIMP-H patients with higher expression of cCRE-associated TFs are significantly linked to worst survival (Fig. 6C). While CIMP-H cases are known to have a high tendency to harbor mutations in BRAF compared to KRAS, no significant association of BRAF or KRAS mutation status was observed in the CIMP-H cases with high and low expression of CIMP-H specific hypermethylated cCRE-associated TFs, indicating that these TF expression profiles are independent of oncogenic mutations (Fig. 6D).

CIMP-H cases with low expression of hypermethylated cCRE-associated TFs, which have better survival, showed a tendency to be enriched for the MSI-H cases with borderline significance (P-value = 0.08, OR = 2.09) (Fig. 6D). Thus, we further tested whether the differences in the OS in CIMP-H cases with high or low expression of the CIMP-H-specific hypermethylation-associated cCRE TFs are related to MSI status (Fig. 6E). The favorable prognosis for CIMP-H patients with low expression of hypermethylated cCRE-associated TFs remained even after adjusting for MSI status (Fig. 6E, Cox regression analyses, HR = 0.45, P-value = 0.047). Furthermore, the association of OS with expression differences of the hypermethylated cCRE-associated TFs is specific to the CIMP-H patients as OS was not associated with CIMP-L and negative patients stratified based on expression differences of these TFs. Taken together, the above analyses show that the CIMP-specific hypermethylation-associated cCRE TFs allow prognostication of CIMP-H subtypes which performs better than MSI-based subtyping.

DNA was extracted from two to three slices (10 µm) from the FFPE samples (Qiagen), quantified using PicoGreen, and checked for quality using the Illumina FFPE QC Kit. The DNA was then subjected to bisulfite conversion (EZ DNA methylation kit, Zymo Research) and controlled by two in-house qPCRs. Next, 250 ng DNA was subjected to Illumina’s FFPE restoration protocol. Data were processed in two batches: the first batch was profiled on Infinium HumanMethylation 450 BeadChip, and the second batch was run on Infinium MethylationEPIC BeadChip. The methylation data retained for analysis consisted in the β-values (ratio of the methylated probe intensity and the overall intensity) of 390,277 probes that are common to both array types.

RNA was extracted from one slice (10 µm) from the FFPE samples (Qiagen RNEasy FFPE kit). Amplification and labeling of total RNA were performed using the GeneChip® PICO Reagent Kit following the manufacturer’s protocol (ThermoFisher 2016, P/N 902,790). Biotin-labeled target samples were hybridized to the Affymetrix HT HG-U133+ PM Array. Target hybridization was processed on the GeneTitan® Instrument according to the manufacturer’s instructions provided for Expression Array Plates (P/N 702,933). Images were analyzed using the GeneChip® Command Console Software (GCC) (ThermoFisher).

Following are the public datasets used in this study: 1) TCGA CRC (colon tumors 279 vs. normal 19) (Methylation and Expression); 2) sessile serrated adenomas (n = 78) (ArrayExpress, E-MTAB-7854) (Methylation); and 3) colon cancer dataset GSE39582 (n = 516) (Expression).

We integrated and normalized the DNA methylation array data from TCGA colon tumors (n = 317), TCGA colon normal (n = 38), tubular adenoma polyps (accrued by Janssen Pharmaceuticals, n = 127), and sessile serrated adenomas (n = 78) (ArrayExpress, E-MTAB-7854) (74). Package “Minfi” was used to process the methylation data (75). See SI Appendix, Materials and Methods, for a detailed description on analysis of methylation data.

Expression dataset for the Janssen polyps was derived using the Affymetrix Array platform (HT HG-U133 Plus PM) and was processed using the affy package (76) with Robust Multichip Average normalization. The Janssen polyp cohort contained 61 samples that had matched methylation and expression profiles. Expression data for both TCGA colon and Janssen polyps were integrated and normalized for removing batch effects by using the “ComBat” subroutine from the SVA package (77) (SI Appendix, Fig. S5G). Analyses of enhancer methylation and associated TFs and target gene expression were done only for the tubular/tubulovillous samples as the expression profiles are available for only these samples.

Estimation of the contribution of immune cell infiltration to CIMP phenotypes was carried out using expression data and CIBERSORT (78). See SI Appendix, Materials and Methods, for details.

See SI Appendix, Materials and Methods, for a detailed description of TF identification and the presence of chromatin marks. In brief, for TF analysis, TCGA level three data matched to both DNA methylation and expression were downloaded directly from the GDC Data Portal (website) by using the ELMER program (33). Limma package (79) was then used to identify differentially methylated probes associated with the various cCRE elements in the promoter and enhancer regions. Differentially hypermethylated cCRE sites were defined as sites that are unmethylated in normal samples, satisfying mean β-value (µβ_normal) <= 0.2 and µβ_cancer-µβ_normal >= 0.2, with FDR <= 0.05. Differentially hypomethylated cCRE sites were defined as sites that are methylated in normal samples, satisfying µβ_normal >= 0.6 and µβ_cancer-µβ_normal <= −0.2 with FDR <= 0.05. We used the TCGA colon tumor data for TF analyses as these samples represent the comprehensive spectrum of CIMP phenotypes while the polyp samples with expression and DNA methylation data mainly consisted of tubular adenomas. The ELMER algorithm was used to identify TFs whose expression is correlated to the methylation of their binding motifs. In order to obtain various chromatin states of interest, we used ChromHMM emission states derived from ChIP-seq of H3K4me3, H3K27ac, H3K4me1, H3K27me3, and H3K9me3 on normal human colon mucosa from the ENCODE project. Accordingly, Illumina methylation array probes overlapping following ChromHMM states summarized into the various activity states of the promoter, enhancer, and constitutive heterochromatin elements. All annotations were performed using hg19 genome version.

TCGA CRC data and GSE39582 (colon cancers, Affymetrix expression dataset downloaded from Synapse, Synapse ID syn2623706) datasets were normalized together using “ComBat” algorithm (77) for batch effect removal (SI Appendix, Fig. S6E). After filtering for samples with known CMS subtypes, there were 208 TCGA CRC samples and 516 GSE39582 samples for analysis. We used CIMP-specific cCRE-associated TF (N = 476) and performed an unsupervised clustering method called K-means clustering to categorize the samples into specific classes. The choice of the optimum number of K-means clusters (N = 3) was guided by using the elbow method and gap statistics.

Survival analyses were carried out on the combined CRC dataset (TCGA and GSE39582) by using the Kaplan–Meier method. In order to identify the low and high expression groups associated with the CIMP-H specific hypermethylated cCRE-associated TF (N = 145), the gene set variation analysis (GSVA) method was used (80). The GSVA method assigns an expression score to each sample by comparing the expression of genes inside and outside of the given gene set. The samples were grouped into high and low expression classes based on the median cutoff applied to the expression scores. Multivariate Cox regression analysis was used to identify the hazard ratio and adjust for MSI, age, gender, and cancer stage.

Y.R.B., J.R., L.V.W., K.E.B., and H.E. designed research; Y.R.B., J.R., L.V.W., K.E.B., S.B.B., and H.E. performed research; Y.R.B., V.K., S.-J.T., K.E.B., S.B.B., and H.E. analyzed data; Y.R.B., S.-J.T., and H.E. prepared figures; and Y.R.B., V.K., R.P., S.P., E.G., O.K., P.G., S.B.B., and H.E. wrote the paper.

V.K., P.G., J.R., L.V.W., and K.E.B. are employees of Janssen Research and Development Inc, and hold stock and stock options in this company. Methylation-specific PCR (MSP): S.B.B. consults for MDxHealth. MSP is licensed to MDxHealth in agreement with Johns Hopkins University (JHU). S.B.B. and JHU are entitled to royalty shares received from sales. Multiple co-authors are employed by Janssen Pharmaceuticals.