Dynamic DNA methylation across diverse human cell lines and tissues

Quantifying DNA methylation

We modified a previously published protocol for RRBS (Meissner et al. 2008) to create a streamlined workflow for this larger-scale implementation (Supplemental Fig. S1). For each reference cytosine sequenced, we computed the percentage of reads in which the cytosine was methylated (remained a C after bisulfite treatment) out of the total reads covering that position (Supplemental Fig. S2). This percent methylated (PM) value represents the percentage of molecules that were methylated at each cytosine.

In replicate growths of the human myeloid cell line K562, we found that lower read depth provided less reproducible measurements of PM between replicates. When restricted to CpGs with at least 10× coverage, the reproducibility of PM measurements across the full range of values improved between replicates (r = 0.987) and resulted in an average of 3.96 PM difference per CpG between replicates (Supplemental Fig. S3A). The resulting PM measurements were highly correlated with values obtained from an array-based methylation assay (Illumina Methyl450K, r = 0.954) (Supplemental Fig. S3B). We performed RRBS on 82 human cell lines and tissues in duplicate (sample information in Supplemental Table S1) and obtained at least 10× coverage for an average of 700,000 CpGs in each sample. The MspI restriction digest utilized for RRBS enriches for CpG-dense regions of the genome, including genes (twofold enrichment) and CpG islands (111-fold enrichment).

Global observations

We found that all samples, regardless of the disease state or tissue type, had similar distributions of methylation among the assayed CpGs (mean of pairwise R² = 0.96) (Supplemental Fig. S3C). In each sample, 5%–15% of assayed CpGs were completely methylated (PM ≥ 90), and 65%–80% of assayed CpGs were unmethylated (PM ≤ 10). This consistency could appear because the same CpGs are always methylated in all samples, or it could result from the same net amount of methylation placed on different loci between samples. Our data set supports the latter; the PM values of individual CpGs varied substantially across cell lines and tissues. Only 4% (27,053) of CpGs are unmethylated (≤10 PM) across all cell lines and tissues that we assayed, and these are located near the transcription start sites (χ², P < 0.0004) of genes with housekeeping functions (P < 1.35 × 10⁻¹⁰). The remaining 670,000 CpGs that we assayed exhibited differential methylation in this study, providing a rich data set for investigating epigenetic patterns.

To characterize cell-type specific methylation patterns, we analyzed 440,974 autosomal CpGs with at least 10× coverage in at least 90% of the samples. We performed unsupervised hierarchical clustering on the PM values for the top 5% of CpGs with the most variable methylation (N = 22,696, σ ≥ 32.6). The samples all paired with their replicates and clustered into clades with distinct methylation patterns, and these clades corresponded to distinct types of samples, namely cancer cell lines, primary cell lines, tissues, and blood leukocytes (Fig. 1A; detailed tree in Supplemental Fig. S4). We divided these most variable CpGs based on whether they are located in gene regulatory regions (<2000 bp from the transcription start site) or in the body of genes (>2000 bp from the transcription start site) and found that both subsets recapitulate the classification of samples in the four clades (Supplemental Fig. S5). To determine if the detection of cell-type specific differences would be confounded by epigenetic variability introduced when cell lines were grown in different ENCODE Consortium laboratories, we obtained growth replicates of the same cell lines from different laboratories. We found that they clustered together based on cell type, not laboratory. Furthermore, for a particular cell type, replicates from within a lab and replicates grown in distant labs were equally similar (Supplemental Fig. S6).

Aberrant methylation across cancer cell lines

The dominant signal in this data set is cancer-specific hypermethylation found at 66,570 CpGs in the cancer cell lines (Fig. 1A), including lines derived from breast, prostate, lung, ovarian, endometrial, liver, and pancreatic cancer, as well as neuroblastoma and several leukemias (Kolmogorov-Smirnov, P < 1 × 10⁻⁷). Of the loci that are hypermethylated across cancers, 48,787 (73%) reside in CGIs (Table 1) and represent a significant portion of the 148,465 island CpGs assayed (33% vs. expected 15%, Fisher's exact test, P < 0.05). The observation of hypermethylation at CGIs across the genome, regardless of their proximity to genes, in 18 diverse cancer cell lines is consistent with reports of a CpG island methylator phenotype (CIMP) that was first described in colorectal cancer (Toyota et al. 1999) and has since been documented in many cancer types (Teodoridis et al. 2008; Liu et al. 2011; Chen et al. 2012; Turcan et al. 2012). An additional 7377 (11%) nonisland CpGs are significantly hypermethylated in cancer, and these reside in promoters and bodies of genes encoding proteins with sequence-specific DNA binding transcription factor activity (hypergeometric enrichment for GO:0003700, P = 5.4 × 10⁻¹⁵) (Table 1). The presence of increased methylation in both the promoter and gene body of these transcription factor genes indicates dysregulation of methylation at these genes in cancer cell lines.

It has been observed that hypermethylation in cancer is enriched at loci that in embryonic stem cells (ESCs) are unmethylated, have bivalent chromatin marks, and are reversibly repressed by the Polycomb Repressive Complex. This observation has led to the model of a stem/progenitor cell signature in cancer (Ohm et al. 2007; Schlesinger et al. 2007; Widschwendter et al. 2007; Easwaran et al. 2012). We determined the overlap between the loci that are hypermethylated across the cancer cell lines in our data set and the location of binding sites for 149 transcription factors that were assayed by chromatin immunoprecipitation sequencing (ChIP-seq) experiments performed by our lab and others in The ENCODE Project Consortium. For each transcription factor that overlapped loci that are hypermethylated in cancer cell lines, we report the enrichment and statistics in Supplemental Table S2. Consistent with previous reports, we found that hypermethylated CpGs have a significant overlap with loci that are bound by SUZ12 in ESCs, a component of Polycomb Repressive Complex 2 (Fisher's exact test, Benjamini–Hochberg (BH)-adjusted, P = 9.27 × 10⁻¹⁴²). The corepressor CTBP2 was also enriched at these sites and significantly overlapped SUZ12-bound loci in ESCs (Fisher's exact test, BH-adjusted, P = 9.75 × 10⁻²⁰). Additionally, we discovered that hypermethylated loci were also enriched for sites where NANOG is bound in ESCs, which is a transcription regulator essential for maintaining pluripotency in ESCs (Fisher's exact test, BH-adjusted, P = 4.06 × 10⁻⁴). The hypermethylated NANOG binding sites do not overlap the SUZ12-bound loci, and the genes nearest the hypermethylated NANOG binding sites are enriched for genes that are overexpressed in human ESCs (P = 6.64 × 10⁻⁵) (Ben-Porath et al. 2008). These observations support the role of embryonic transcriptional regulators in directing aberrant methylation that leads to cancer but add complexity to the model. In contrast to the polycomb repressive complex recruiting hypermethylation in cancer to persistently silence differentiation-inducing genes, our results suggest that hypermethylation also occurs at NANOG-bound loci that are active in ESCs. While these are seemingly opposing effects, it is possible that both mechanisms lead to tumorigenesis by silencing genes that lead to differentiation, as well as silencing genes that control the stable, nonneoplastic division of a pluripotent cell. This is supported by the observation that NANOG binding site hypermethylation in cancer occurs at genes enriched for transcription factor activity (P = 9.78 × 10⁻⁵), including developmental regulators such as FOXD3, a negative regulator of cell cycle (Abel and Aplin 2010), and CDX2, whose silencing is associated with cancer progression (Huang et al. 2012; Knosel et al. 2012). Figure 1B depicts another example, the SFRP2 gene, a modulator of Wnt signaling, whose promoter is unmethylated and bound by NANOG in the H1 ESCs where the gene is expressed. In cancer, the NANOG binding site in the promoter is methylated, and the gene is silenced.

Hypomethylation across the genome has been reported in cancer (Irizarry et al. 2009; Hansen et al. 2011; Berman et al. 2012; Hon et al. 2012). Although RRBS enriches for CpG-rich regions of the genome that tend to be hypermethylated in cancer, we also query thousands of positions in low CpG-density regions. We detected 6691 positions that were significantly hypomethylated across the cancer types when compared to the primary cell lines and tissues (Kolmogorov-Smirnov, P < 1 × 10⁻⁷). We discovered that these hypomethylated loci were colocalized in the genome and identified 114 independent megabase windows that had significantly more hypomethylated loci than expected, taking into account the nonrandom genomic coverage of the assay (Fisher's exact test, P = 2 × 10⁻⁵) (genomic coordinates listed in Supplemental Table S3). These megabase-size hypomethylated domains were significantly depleted of genes (binomial, P < 3.48 × 10⁻⁴⁰) and were significantly enriched in the 10-Mbp ends of chromosomes and near the centromere, although they are not usually directly adjacent to the telomeres or centromeres (binomial, p = 4.41 × 10⁻¹⁹) (Fig. 2A). Recent reports have found that some tumors exhibit hypomethylation corresponding to lamina-associated domains (Berman et al. 2012), but we did not observe an enrichment for lamin B1 binding in the hypomethylated domains we identified. It has also been reported that hypomethylated regions in cancer correspond to H3K27me3 (Hon et al. 2012; Statham et al. 2012), and we found that the domains we identified overlap long tracks of H3K27me3, as well as corresponding stretches of EZH2 binding. We compared the presence of H3K27me3 in these regions between cancer cell lines and noncancer cell lines and discovered that the long tracks of H3K27me3 were specific to the cancer cell lines (χ² test, P = 2.4 × 10⁻²⁷). Brinkman et al. demonstrated that when DNA methyltransferases are knocked-out in ESCs, broad local enrichments (BLOCs) of H3K27me3 appear in place of high levels of methylation (Brinkman et al. 2012). It is plausible that the same process is directing H3K27me3 BLOC formation at these hypomethylated domains in cancer. Among the ChIP-seq experiments comprising 149 transcription factors, we did not identify any factors that were particularly enriched in these domains. Notably, these hypomethylated domains occasionally contain a gene that is expressed in specific cancer samples, and those samples exhibit gene-body methylation within the hypomethylated domain corresponding to the gene's expression. This indicates that gene expression and gene-body methylation are not occluded from these regions by the unmethylated tracks of H3K27me3 surrounding the gene. The prevalence of these domains across cancer types warrants further investigation of these regions as predictive biomarkers and to uncover the mechanisms driving these massive defects. Figure 2B depicts an example of a domain at the end of the q-arm of chromosome 22, where a 2-Mb gene-depleted region is specifically hypomethylated across cancer cell lines and exhibits long tracks of cancer-specific H3K27me3 and EZH2 binding. This hypomethylated domain is flanked by methylation corresponding to the gene-body methylation of expressed genes.

The single nucleotide resolution of bisulfite sequencing allows us to detect both methylation and DNA sequence variants in the same molecules, and we used this information to identify loci with allele-specific or allele-biased methylation (Gertz et al. 2011). We identified 1144 CpGs that were adjacent to an allelic variant and exhibited allele-biased methylation in the noncancer samples from different tissues and individuals. We found that 1027 (90%) of these CpGs, which are allelically methylated in noncancer samples, exhibit aberrant methylation in cancer cell lines (Kolmogorov-Smirnov, FDR-adjusted, P < 0.05) (Supplemental Fig. S7A). This aberrant methylation occurs as either hypermethylation (gain of methylation on the unmethylated allele) (example in Supplemental Fig. S7B) or hypomethylation (loss of methylation on the methylated allele) (example in Supplemental Fig. S7C). Loss-of-imprinting has been previously reported at particular imprinted loci in cancer (Feinberg et al. 2002; Bjornsson et al. 2007; Feinberg 2007; Monk 2010), and these observations demonstrate that the majority of allelically methylated loci are dysregulated in cancer, regardless of whether they are imprinted.

Cell-type specific methylation

When we isolated and subjected noncancer samples to unsupervised hierarchical clustering of the top 5% of CpGs with the most variable methylation (N = 22,152, σ ≥ 24.5), the samples again clustered into clades corresponding to distinct types of biological samples: tissues, primary cell lines, embryonic cell lines, and blood leukocytes (Fig. 3A; detailed tree in Supplemental Fig. S8). The blood leukocyte-derived samples, including both peripheral blood leukocytes and Epstein-Barr virus (EBV)-immortalized lymphoblastoid cell lines, clustered together and displayed a distinct pattern of methylation at 6511 CpGs (Kolmogorov-Smirnov, P < 1 × 10⁻⁷). Blood-specific hypermethylation occurs at CGIs (Table 1), but unlike the ubiquitous CGI methylation in cancer cell lines, only a select subset of CGIs (66%) exhibit hypermethylation when compared to tissue and primary cell lines. Epigenetic regulation of genes involved in blood leukocyte development was evident and included specific methylation at genes involved in regulation of body fluid levels (GO:0050878, P = 3.25 × 10⁻⁴), blood coagulation (GO:0007596, P = 5.33 × 10⁻⁴), hematopoietic or lymphoid organ development (GO:0048534, P = 8.52 × 10⁻⁴), and B cell activation (GO:0042113, P = 9.24 × 10⁻⁴).

We identified seven tissue types that were represented by both primary cell lines and primary tissues and performed ANOVA to identify CpGs whose methylation is significantly associated with the tissue of origin (regardless of whether it has been grown in culture). We identified 117,795 CpGs whose methylation was significantly associated with the tissue of origin and was consistent in both the primary cell lines and the primary tissues (FDR < 0.05) (subset visualized in Fig. 3B). These data support previous observations that there is a large number of tissue-specific differentially methylated regions (tDMRs) (Rakyan et al. 2008) and that primary cell lines can serve as models for understanding epigenetic tissue-specific gene regulation at a large number of loci. However, as described above, we found that unsupervised hierarchical clustering of the most variable CpGs across samples divided the primary cell lines from the tissues (Figs. 1A, 3A). We identified 2238 CpGs that significantly discriminate primary cell lines from tissue samples (Kolmogorov-Smirnov, P < 1 × 10⁻⁷) (Table 1; loci listed in Supplemental Table S4). These methylation differences associated with cell culture occur predominately as unmethylated CpGs in the bodies of genes involved in regulating cellular proliferation (GO:0042127, P = 5.17 × 10⁻⁴). Studies of DNA methylation that use cell lines as model systems could use this list to reduce false positives due to the epigenetic effects of cell culture.

Context-dependent DNA methylation signatures of gene expression

The current models describing the relationship between DNA methylation and gene expression indicate that promoter methylation is associated with gene silencing, and gene body methylation is associated with expression (Doerfler et al. 1989; Jones and Baylin 2002; Lorincz et al. 2004; Ball et al. 2009; Illingworth et al. 2010; Maunakea et al. 2010; Aran et al. 2011; Deaton et al. 2011). RNA-seq has been performed as part of The ENCODE Project on several of the samples included in our DNA methylation study, providing an opportunity to explore and quantify this relationship. We divided the CpGs near genes into four categories that were distinguished based on whether or not they reside in CGIs and whether they are near the TSS (<2 kb upstream or downstream) or far from the transcription start site in the gene body (≥2 kb). We computed the Pearson correlation coefficient between RPKM values (reads per kilobase of transcript per million reads) (Mortazavi et al. 2008) measured by RNA-seq and PM values measured by RRBS. In this data set, the vast majority of CpGs close to the TSS, regardless of whether they reside in CGIs, are negatively correlated with gene expression (median r = −0.3756), i.e., increased methylation is associated with lower levels of gene expression (Fig. 4A,B). In contrast, the nonisland CpGs far from the TSS in the body of the gene are positively correlated with expression (median r = 0.44450), i.e., increased methylation is associated with higher levels of gene expression (Fig. 4D). For these promoter and nonisland gene body CpGs, the current model of the relationship between methylation and expression holds across these cell lines. However, we also identified an exception to the expected pattern: island CpGs residing in gene bodies displayed a bimodal distribution of correlation with expression levels (Fig. 4C). This indicates that a substantial subset of CGIs in the gene body are negatively correlated with expression, similar to promoters, rather than positively correlated with expression like other gene body CpGs.

To investigate gene regulatory processes that could account for this finding, we integrated data sets from The ENCODE Project, including CAGE tag sequencing, histone modification, and transcription factor ChIP-seq. Recent studies have proposed that cryptic or alternative promoters, marked by H3K4me3 and CAGE tags, may appear as promoter-like methylation in gene bodies (Illingworth et al. 2010; Maunakea et al. 2010; Deaton et al. 2011). We found that H3K4me3 and CAGE tags account for 8.5% (479/6155) of the gene body CGI CpGs that were negatively correlated with expression, providing evidence that this subset of CGIs indeed reside in alternative promoters (Fisher's exact test, P-value = 1.06 × 10⁻⁷ [H3K4me3] and 0.009124 [CAGE tags]).

We sought to identify other function elements that could explain the remaining gene body CGIs that have epigenetic regulation similar to promoters. A recent study showed that DNA binding factors influence DNA methylation and that the methylation signature could be used to identify CpG-poor distal enhancers in the mouse genome (Stadler et al. 2011). The causes and consequences of gene-body methylation are not well understood, making it unclear whether active enhancers could be unmethylated when embedded in the densely methylated gene body of an expressed gene. We investigated whether the unmethylated CpG-rich regions that we observe within methylated gene bodies might be intragenic enhancers. We performed ChIP-seq for the transcriptional coactivator EP300, a factor known to bind to transcriptional enhancers (Visel et al. 2009). The vast majority of CpGs in EP300 binding sites were unmethylated (PM ≤ 10; 99.3% in HepG2, 98.2% in GM12878, 99.6% in hESC H1). When gene body CGIs contain EP300 binding sites, they are more strongly inversely correlated with gene expression than are CGIs not bound by EP300 (Kolmogorov-Smirnov test, P-value < 2.2 × 10⁻¹⁶) (Supplemental Fig. S9). These EP300-bound CGI intragenic enhancers account for an additional 8% (452/5676) of gene body CGIs that are negatively correlated with expression. The signature of CpG-rich active enhancers in gene bodies was not observed at CpG-poor regions, indicating that this escape from gene-body methylation is associated with CGIs. To identify other transcription factors associated with these unmethylated gene body CpG islands that may reveal the role of the remaining loci that are not associated with EP300 binding, H3K4me3, or CAGE tags, we overlapped them with the binding sites for the remaining 148 transcription factors in ChIP-seq data sets generated by The ENCODE Project Consortium. We did not identify any significant enrichment for specific transcription factors. The catalog of EP300-bound enhancers and other nonpromoter regulatory elements is not comprehensive. It is possible that these are enhancers bound by EP300 below the sensitivity of ChIP-seq, or they could be regulatory regions bound by transcription factors or noncoding RNAs not yet studied by The ENCODE Project Consortium.

Together, these results support a revised model of the expected DNA methylation state found at expressed and silenced genes, depicted in Figure 4, E and F. In expressed genes, active intragenic enhancers bound by EP300 appear as patches of unmethylated CGIs amid the dense methylation found in the body of the expressed genes, and these enhancers have methylation that is concordant with the unmethylated CpGs several thousand base pairs away in the promoter and near the TSS. The reverse patterns are associated with silenced genes. This revised model changes our expectations of the type of methylation we find in gene regions and helps to more accurately interpret gene body methylation.

Non-CpG cytosine methylation

DNA methylation predominately occurs at CpG dinucleotides in the human genome, but there have been recent reports that non-CpG cytosine methylation occurs at a lower, but appreciable, level in embryonic and pluripotent cells (Lister et al. 2009; Ziller et al. 2011). We identified 56,287 cytosines that were not located at CpG positions in the reference genome that exhibited methylation (PM > 10) in at least one sample. We eliminated 30,773 (55%) of these methylated cytosines from further analysis because they were adjacent to genetic variants in these samples that created a CpG dinucleotide that became methylated, a finding that suggests that there is a large number of polymorphic CpGs that show epigenetic diversity between individuals. To reduce false-positives due to stochastic errors in bisulfite conversion, we identified 2466 non-CpG cytosines that were methylated (PM > 10) in both replicates of any sample and determined how many of these loci were methylated in each sample (Supplemental Table S5; data in Supplemental Table S6). We found the largest number of methylated non-CpG cytosines (N = 1418) in the human embryonic stem cell line H1 (H1-hESC), consistent with previous reports in this cell type (Fig. 5A). Adult human brain tissue had the second highest number of methylated non-CpG cytosines, with 666 non-CpG cytosines methylated in both replicates, and was more than twice as high as the following sample (Supplemental Table S5). The non-CpG methylation we observed in the brain samples occurred at a different set of loci than the non-CpG methylation observed in ESCs (Fig. 5A). This was unexpected in light of the recent report of the near complete absence of non-CpG cytosine methylation in human somatic cell types, although adult human brain tissue was not studied (Ziller et al. 2011). The observation of non-CpG methylation in somatic tissues challenges current hypotheses that a pluripotent-specific regulation or noise is involved in establishing this mark. The other samples with more than 200 non-CpG cytosines methylated in both replicates included skeletal muscle tissue from two different individuals and placental tissue (Fig. 5A; Supplemental Table S5). To further investigate non-CpG cytosine methylation in the brain, we performed RRBS on brain samples from 24 additional individuals. These fresh-frozen samples were collected post-mortem from the dorsolateral prefrontal cortex (DLPFC) of healthy control donors as part of the Pritzker Neuropsychiatric Disorders Research Consortium. The non-CpG cytosines that were methylated in the first brain sample were also methylated (PM > 10) in these additional brain samples (Fig. 5B; data in Supplemental Table S7). The non-CpG loci that are methylated in adult brain tissue are distinct from the set of non-CpG cytosines that were methylated in ESCs (Fig. 5B). The non-CpG cytosine methylation that we observed in the brain predominately occurs at CAC trinucleotides (Fig. 5D), which is different from the CAG trinucleotide context that we and others have observed in ESCs (Fig. 5C; Lister et al. 2009; Ziller et al. 2011). Non-CpG methylation was recently reported in mouse frontal cortex (Xie et al. 2012) in a similar sequence context that we find in human brain, suggesting that mice may serve as a valuable experimental model for understanding this new methylation pattern. We used GREAT (Genomic Regions Enrichment of Annotations Tool) (McLean et al. 2010) to determine if these methylated non-CpG cytosines are associated with any functional enrichment and found that loci in brain tissue are found near genes enriched for “blood vessel development” (GO:0001568, GREAT hyper FDR q-val = 5.34 × 10⁻⁶), while the loci methylated in hESC are found near genes related to small GTPase regulator activity (GO:0005083 hyper FDR q-val =4.91 × 10⁻⁴), suggesting that these events are associated with different gene regulatory processes in each cell type.

Cell lines and tissues

Samples included in this study are listed in Supplemental Table S1. Detailed information about the samples can be obtained from the ENCODE Common Cell Types websites at (http://genome.ucsc.edu/ENCODE/cellTypes.html).

Reduced representation bisulfite sequencing experimental procedure

We modified the previously published protocol for RRBS (Meissner et al. 2008) to create a streamlined workflow for this larger-scale implementation. We designed the reactions to eliminate the phenol extraction and ethanol precipitation steps as well as one of the gel extraction steps. We also changed the PCR conditions to amplify fragments with diverse GC content and a broad range of sizes uniformly. A protocol overview is depicted in Supplemental Figure S1, and specifics are provided as follows. We used the Qiagen DNeasy Blood and Tissue Kit to extract DNA. We then digested 1 μg genomic DNA with 1 μL 20U/μL MspI restriction enzyme (New England Biolabs [NEB]) in 1× NEBuffer 2 in a total reaction volume of 50 μL. This reaction was incubated at 37°C for 30 min, followed by heat inactivation at 80°C for 20 min. We then filled in the overhangs and added a 3′ A tail by adding dNTP mix to 33 μM and 1 μL 5U/ μL Klenow Fragment (3′–5′ exo-) in a total reaction volume of 55 μL. This reaction was incubated at 37°C for 30 min, followed by heat inactivation at 75°C for 20 min. We then purified the DNA with a Qiagen MinElute column. We purchased two methylated DNA oligonucleotides from IDT (www.idtdna.com) as follows: ilAdap Methyl PE1 (ACACTCTTTCCCTACACGACGCTCTTCCGATC*T) and ilAdap Methyl PE2 (5′P-GATCGGAAGAGCGGTTCAGCAGGAATGCCGA*G), where all C's are 5-methyl cytosine DNA nucleotides, 5′P indicates a 5′ phosphate, and an asterisk indicates a phosphorothioate bond. We then annealed these oligos to form a stock of 40 μM duplex DNA adapters in a reaction containing 40 μM ilAdap Methyl PE1, 40 μM Methyl PE2, 1× T4 DNA Ligase Buffer (NEB) in a total volume of 50 μL. We incubated this reaction at 95°C for 5 min, then 70°C for 1 min, then 60°C for 1 min, then 50°C for 1 min, then 40°C for 1 min, then 30°C for 1 min. We stored these annealed adapters at −30°C for future use. We ligated the annealed methylated Illumina adapters in a reaction containing 10 μL purified DNA, 1× T4 DNA Ligase Buffer (NEB), 1 μL 400U/μL T4 DNA Ligase (NEB), and 1 μL 40 μM annealed methylated adapters in a total volume of 20 μL. This reaction was incubated at 20°C for 15 min, followed by heat inactivation at 65°C for 10 min. We electrophoresed the 20 μL ligation reaction in a 2.5% Seaplaque Agarose (Lonza) gel. The desired restriction fragments are between 40 and 120 bp, and the adapters add 33 bp onto each end of the restriction fragments, so we used a razor blade to isolate the agarose gel section containing DNA between 106 and 186 bp, while avoiding the adapter self-ligation products that appear <100 bp. We then purified the DNA using a Qiagen Qiaquick Gel Extraction kit as described in the manufacturer's instructions, except that we did not heat the gel fragment to dissolve it, and we eluted the purified DNA from the column using 22 μL buffer EB. We then used 20 μL of this purified DNA in the sodium bisulfite conversion, which was performed using the EZ DNA Methylation Gold Kit (Zymo Research). We purchased PCR primers that would amplify the adapter-ligated DNA and add the cluster generation sequences to the amplicons for Illumina sequencing. We purchased these PCR primers from IDT as follows: ilPCR PE1 (AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATC*T) and ilPCR PE2 (CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATC*T), where an asterisk indicates a phosphorothioate bond. The DNA that was purified from the bisulfite conversion kit was then PCR-amplified in a reaction containing 5 U Platinum Taq DNA Polymerase (Invitrogen), 10× PCR Buffer without MgCl₂ (Invitrogen), 2 mM MgCl₂, 0.5 μM ilPCR PE1 DNA oligonucleotide primer, 0.5 μM ilPCR PE2 DNA oligonucleotide primer, 0.5 mM each dNTP and 0.5 M Betaine (Sigma-Aldrich) in a total reaction volume of 50 μL. The reaction was incubated at 98°C for 1 min, followed by 20 cycles of (95°C for 30 sec and 62°C for 3 min). We confirmed the amplification and correct product size range by running one-fifth of the reaction on a 2% agarose gel. We then purified the remaining PCR product with a Qiagen Qiaquick column, eluting in 25 μL buffer EB, and quantified the purified product using the Quant-IT High Sensitivity dye kit (Invitrogen) on a Qubit fluorometer (Invitrogen). We typically obtained 250–750 ng of library material. We then diluted each library to 10 nM and proceeded to sequence each library in a single lane on the Illumina Genome Analyzer IIx sequencing machine according the manufacturers instructions. We typically achieved better quality scores and alignment with slightly lower cluster density compared to other libraries with more even base-representation, and we empirically determined that clustering the sample at 5 pM was optimal.

Sequence alignment and calculating percent methylated value for each cytosine

The sequence data for this project were acquired between June 2009 and December 2010. On our Illumina Genome Analyzer IIx sequencing machine, we sequenced one library per lane and obtained between 8 million and 31 million single-end 36-bp reads per library. We aligned these reads to a modified reference genome sequence that was created to reflect both the reduced representation of the genome due to the MspI restriction digest as well as the sodium bisulfite conversion which creates a T in the sequencing reads rather than a C at all unmethylated bases. To create this reference, we first parsed the hg19 reference genome to identify all of the MspI restriction fragments <500 bp. We then isolated the 36-bp ends of these fragments into a fasta file and converted every C in the reference to a T and recorded the position of these reference cytosines in the name of the reference sequence. To achieve optimal alignment that is not biased by the methylation state of a molecule, we also created a copy of our sequencing read files, converted every C in the read to a T, and recorded the position of these read cytosines in the name of the read (Supplemental Fig. S2). We then used bowtie (Langmead et al. 2009) to align these converted reads to the custom reference sequence and required that the alignment be optimal and unique in the reference and only align to the proper strand (bowtie options –best, -m 1,–norc). On average, we uniquely aligned 53.2% (200,000/375,603) of the genomic MspI digest restriction fragments in the selected size range (40–120 bp), which resulted in coverage of an average of 1.2 million CpGs in each sample. This is only 8.6% of the 14 million nonrepetitive CpGs in the human genome but represents a 1.9-fold enrichment for genic regions and a 111-fold enrichment for CGIs. We then parsed the alignment file and the encoded read and reference names to determine how many reads covered each reference cytosine position and what percentage of those reads contained a C at each reference cytosine position (Supplemental Fig. S2). This percent methylated value approximates the percentage of molecules in the sample that were methylated at each individual cytosine. We compute the bisulfite conversion rate of each sample by determining the percent of non-CpG cytosines that are methylated (PM ≥ 10), which is an underestimate of the conversion rate in samples with biological non-CpG methylation. Each sample must meet quality control criteria before data release, including a bisulfite conversion rate ≥98.5%, a complex library with more than 500,000 CpGs with at least 10× coverage, and a correlation coefficient of greater than 0.9 between replicates.

Methyl 450 array methods

Illumina Methylation450 arrays were run using standard Illumina protocols. Briefly, 500 ng of DNA from each cell line was bisulfite-converted with the Zymo Research EZ DNA Methylation kit, amplified, hybridized, and stained with standard Illumina reagents. The intensity data were imported into Illumina's GenomeStudio software, and standard beta scores were exported and used in the analysis.

Analysis of methylation across samples

Once we compiled the percent methylated values for all cytosines shared across samples, we then performed extensive analysis of the trends in these data. Statistical associations including mean calculations, standard deviation calculations, Pearson correlations, binomial tests, Fisher's exact tests, χ² tests, hypergeometric tests, and Kolmogorov-Smirnov tests were calculated using Matlab (The Mathworks, Inc.), and the statistical package R (The R Foundation for Statistical Computing; http://www.R-project.org). Clustergrams were created using the average linkage of Euclidean distance in Cluster3.0 (de Hoon et al. 2004) and visualized using Java TreeView 1.1.4r3 (Saldanha 2004). For the annotation of cytosine positions relative to gene features, we used the genomic coordinates for gene features from the hg19 refGene table of the UCSC Genome Browser (Fujita et al. 2011; Raney et al. 2011). Similarly, we used the genomic positions of the CpG islands track on the UCSC Genome Browser to annotate CGI occupancy. To measure the Pearson correlation between methylation and expression, we excluded CpGs whose methylation did not vary by 10 PM across the cell lines to avoid spurious correlations to noise in the methylation measurements. For identifying transcription factor binding sites associated with loci hypermethylated across cancer cell lines, we used the compiled supertrack data set containing binding sites for 149 transcription factors from ENCODE ChIP-seq experiments, which is available from the UCSC Genome Browser (http://www.genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=305048059&c=chr2&g=wgEncodeRegTfbsClusteredV2). For analyzing the binding sites of the specific transcription factors, we used individual ChIP-seq data sets. The data for SUZ12 binding sites in H1-hESC are available from the UCSC Genome Browser (http://www.genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=305048059&g=wgEncodeSydhTfb). The data for NANOG binding sites in H1-hESC are also available from the UCSC Genome Browser (http://www.genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=305048059&g=wgEncodeHaibTfbs). For characterizing the hypomethylated domains found across cancer cell lines, we used the H3K27me3 and EZH2 binding ChiP-seq data collected by the Broad Institute as part of The ENCODE Project, which were obtained from the UCSC Genome Browser (http://www.genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=286312585&c=chr4&g=wgEncodeBroadHistone). The nuclear lamina-associated domains data were obtained from the UCSC Genome Browser (http://www.genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=305048059&c=chr2&g=laminB1Super). The RNA-seq data for the HeLa, hESC H1, K562, HepG2, and GM12878 cell lines were collected as part of The ENCODE Project and can be found under “RNA-seq from ENCODE/Caltech” on the UCSC Genome Browser (http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=193248635&c=chr10&g=wgEncodeCaltechRnaSeq). The H3K4me3 ChIP-seq data for the HeLa, K562, HepG2, and GM12878 cell lines were collected as part of The ENCODE Project and can be found under “Histone Modifications by ChIP-seq from ENCODE/University of Washington” on the UCSC Genome Browser (http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=203697013&c=chr5&g=wgEncodeUwHistone). The CAGE tag data from whole cell polyA+ fractions of the HeLa, hESC H1, K562, HepG2, and GM12878 cell lines were collected as part of The ENCODE Project and can be found under (http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=210114571&c=chr21&g=wgEncodeRikenCage). The EP300 ChIP binding site information was generated by our group as part of The ENCODE Project and can be found under “ENCODE Transcription Factor Binding Sites by ChIP-seq from HudsonAlpha Institute” on the UCSC Genome Browser (http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=210114571&c=chr21&g=wgEncodeHaibTfbs). Binding sites for EP300 identified in GM12878, H1 hESC and HepG2 were combined to create a list of potential enhancers. Genetic polymorphism that created CpG dinucleotides were identified as those positions where the bowtie alignment identified the same mismatched base in at least 10% of the reads with a minimum read depth of five. Functional annotation and enrichment of genes was obtained using the gene ontology search program Gorilla (http://cbl-gorilla.cs.technion.ac.il/) (Eden et al. 2009), using the option that calculates enrichment in a target gene list over a background list of all genes covered in the RRBS libraries. The motif logos representing the sequence context of non-CpG cytosine methylation were created using WebLogo (http://weblogo.berkeley.edu) (Crooks et al. 2004).