Method of detecting cancer through generalized loss of stability of epigenetic domains, and compositions thereof

Generalized Loss of Stability of Epigenetic Domains Across Cancer Types

The following experimental protocols and materials were utilized.

Custom Array Design.

Custom bead arrays were created based on the 151 most statistically significant cDMRs identified in Irizarry et al. 4 for which successful design probes using the Illumina protocol 7 could be designed. Cryogenically stored freshly frozen samples were obtained from the Cooperative Human Tissue Network (NCI, Bethesda, Md.), the National Wilms Tumor Study tissue bank (Edmonton, Alberta, Canada) and the Johns Hopkins Hospital, under an IRB-approved waiver of consent. A pathologist validated classification of all samples independently and blindly. Samples were used for DNA purification, the resulting DNA bisulfite treated, and the custom GoldenGate assay performed. The raw intensity data was processed using quantile normalization 28 and ruled out batch effects following the procedures described by Leek et al. 29 . Differences in variance were tested using a standard F-test.

Custom Methylation Array Processing and Analysis.

The raw intensity data from the Cy5 and Cy3 channels representing methylated and umethylated DNA was quantile normalized 2 separately, and methylation level was calculated as the ratio of the Cy5 intensity over the sum of the intensities from both channels. To control for array quality, arrays for which the average of the median log intensities from the two channels was small (<7), or for which the median absolute deviation of the overall methylation signal was small (<1.9) were removed from the dataset. Differences in methylation variability were measured and tested using an F-test. Differences in mean methylation levels were measured and tested using a t-test. Significance was taken as 0.01.

Illumina HumanMethylation27k™ Array Analysis.

A publicly available dataset was downloaded of methylation levels of 22 matched colon normal/tumor samples assayed using Illumina's HumanMethylation27 K™ array (Gene Expression Omnibus accession number GSE17648). The methylation measurements were used with no further preprocessing. Differences in methylation variability were measured and tested using an F-test. Differences in mean methylation levels were measured and tested using a t-test. Significance was taken as 0.01.

Whole Genome Bisulfite Sequencing.

Bisulfite sequencing libraries were prepared using the approach previously described by Bormann Chung et al 30 , with some modifications so the entire protocol is provided in detail here. 5 ug of genomic DNA was sheared using a Covaris E2 sonicator. 40 ng of Alul-digested unmethylated A DNA was added to each sample in order to monitor the efficiency of bisulfite conversion. Sample DNA ends were then repaired using 1× End Polish Buffer, 400 nM each of dATP, dGTP and dTTP (leaving out dCTP), 40 U of End Polishing Enzyme 1 (Life Tech) and 80 U of End Polishing Enzyme 2 (Life Tech). Adaptor sequences, as detailed by Bormann Chung et al. 3 , were then ligated onto the samples, using 1×T4 ligase buffer, 4.5 uM of methyl-protect PI adaptor, 4.5 uM of P2 adaptor and 50 U of T4 ligase. Nick translation was performed in IX Exo-Klenow buffer, 0.5 mM dNTPs containing methyl dCTP and 20 units of Exo-Klenow-Fragment (Ambion) for 1 hour at 16° C. 500 ng aliquots of the resulting product was then bisulfite converted using formamide as an additional denaturant 4 . 24 uL of formamide was added to an equal volume of DNA and incubated at 95° C. for 5 minutes. Subsequently, 100 uL of Zymo Gold bisulfite conversion reagent (Zymo) was added, and the mixture was incubated for 8 hours in 50° C. Samples were then desulphonated and purified using spin columns following the EZ-DNA Zymo Methylation-Gold protocol. 5 uL of the bisulfite converted library was amplified in 1×PCR buffer, 0.2 mM dNTP, 1 mM each of the standard SOLiD™ fragment library primers, 5 U of Taq (Denville), and 0.25 U of Pfu Turbo Taq (Stratagene). The DNA was subjected to 8 cycles of PCR and the resulting product was purified using AMpure™ SPRI beads (Beckman Genomics). The libraries were then sequenced on the SOLiD™ 3+ platform yielding 50 base pair reads. Corresponding colon cancer and normal mucosa samples were sequenced simultaneously in adjacent flow cells.

Alignment of SOLiD™ Sequencing Reads from Bisulfite-Treated DNA.

A custom alignment tool for SOLiD™ reads was developed from bisulfite-treated DNA. Reads were aligned with the aid of a spaced-seed index of the genome while biasing neither toward nor against methylated cytosines in CpGs described previously 31 to work with reads from SOLiD™ 3 Plus instruments. Note that aligners can introduce a bias when an unmethylated C (which becomes a T) is penalized for aligning to a C in the genome, but a methylated C (which remains a C) is not. The opposite bias can also occur, e.g., if all CpGs are converted to TpGs in the reference prior to alignment. Other projects address this in part by additionally converting Cs to Ts in the reads. But this approach is not applicable to the colorspace reads generated by the SOLiD™, for which nucleotide positions encoding Cs cannot be accurately determined prior to alignment. The aligner used here leaves each read as-is but penalizes neither C-to-C norT-to-C partial alignments in CpGs.

For alignment the approach taken by the BSMAP tool was extended. The approach used supports a broad range of spaced-seed designs and extends the BSMAP approach to additionally allow alignment of SOLiD™ colorspace reads. C/T bias is avoided by creating and storing multiple copies (potentially) of each reference subsequence indexed, one copy for each distinct assignment of Cs and Ts to genomic Cs or CpGs present in the subsequence. For colorspace reads, the algorithm extracts subsequences of colors (rather than nucleotides), and a copy it created for every distinct assignment of Cs and Ts to genomic Cs present in any nucleotide overlapped by any color in the extracted subsequence.

For these experiments, the alignment algorithm is configured to remove the penalty associated with either a C or a T aligning to a C in a CpG, and to treat non-CpG Cs in the genome as Ts. This policy removes bias from CpG methylation measurements, but assumes that there is little or no non-CpG cytosine methylation. Where non-CpG cytosine methylation occurs, this approach is more likely to fail to find alignments overlapping the methylated cytosine, and the resulting consensus will contain an anti-methylation bias at the methylated cytosine.

The algorithm was used to align a total of 7.79 billion reads obtained from 8 runs of a SOLiD™ 3 Plus instrument against a reference sequence collection consisting of the GRCh37 human genome assembly (including mitochondrial DNA and “unplaced” contigs) plus the sequence of the spiked-in phage genome. Each read obtained from the SOLiD™ 3 Plus instrument consists of a primer nucleotide followed by a string of 50 “colors,” where each color encodes a class of dinucleotides according to the SOLiD™ colorspace encoding scheme. Prior to alignment, the initial primer base and 5′-most color were trimmed from all input reads, yielding a string of 49 colors. The alignment policy was selected to guarantee that all alignments with up to 3 color mismatches would be found, and some but not all alignments with 4-6 color mismatches would be found. The alignment of a T or C to a C in a CpG does not incur a mismatch penalty (except in some cases where a sequencing error is also present). The policy was also set to distinguish between reads that align uniquely and those that align non-uniquely.

Color-to-color alignments are decoded into nucleotide alignments with a Viterbi-like algorithm. The final alignment, when expressed in nucleotides, is one character shorter than the input read, i.e. 48 nucleotides long. Alignments for reads aligning non-uniquely are ignored in subsequent stages. Alignment results are summarized in Table 1 and FIG. 1 .

For the data of Table 1, alignment was performed against a collection of reference sequences consisting of the GRCh37 human genome assembly, including mitochondrial DNA and “unplaced” contigs, plus the sequence of the spiked-in A DNA. A read is said to align “uniquely” if it has exactly one valid alignment to the reference according to the alignment policy. A read is said to align “non-uniquely” if it has more than one valid alignment according to the alignment policy. A read “fails” to align if it has zero valid alignments. Cytosine conversion percentage is estimated as the fraction of high-quality evidence from unique A phage alignments indicating lack of methylation.

FIG. 1 shows the number of reads and alignments obtained from SOLiD™ 3+ bisulfite sequencing. The stacked bar chart illustrates the number of reads sequenced per flowcell, with the colors indicating uniquely reads, unaligned reads and non-uniquely aligned reads. A total of 7.79 billion reads were obtained from 8 runs (16 flowcells) of a SOLiD™ 3+ instrument.

After alignment, a series of scripts extract and summarize CpG methylation evidence present in the unique alignments. The evidence is compiled into a set of per-sample, per-chromosome evidence tables. Alignments to the A phage genome are also compiled into a separate table. A piece of CpG “evidence” is created when an alignment overlaps the cytosine position of a CpG in the reference and the overlapping nucleotide in the alignment is either a T (indicating a lack of methylation) or a C (indicating presence of methylation). Once a piece of evidence is extracted from a unique alignment, it is subjected to a filter. The filter removes evidence that is either refuted by one or both of the overlapping colors from the original read, or is within 4 positions of either end of the nucleotide alignment. All evidence that passes the filter is added to the CpG summary table. A record in the table summarizes, for a given CpG: the filtered evidence nucleotides that aligned to it, the filtered quality values (i.e. of the two colors overlapping the evidence nucleotide) that aligned to it, the number of distinct alignment positions from which filtered evidence was taken, the “mapability” of the CpG and surrounding bases (i.e. the number of 50-mers overlapping the CpG that are unique up to 3 mismatches), (e) the local CG content of the bases surrounded the CpG. The mapability measure for each genome position is pre-calculated using Bowtie 8 .

Table 2 summarizes the amount and type of evidence extracted at each stage. For Table 2, a piece of CpG evidence occurs when an alignment overlaps the cytosine position of a CpG in the reference sequence and the overlapping alignment nucleotide is either a T (indicating a lack of methylation) or a C (indicating presence of methylation). A filter is applied to remove nucleotide evidence that (a) is refuted by one or both of the overlapping colors from the original read, or (b) is within 4 positions of either end of the nucleotide alignment.

Table 3 summarizes the resulting CpG evidence coverage with respect to the GRCh37 human genome assembly for each sample. For Table 3, each sample is sequenced on two flowcells of a SOLiD™ 3+ instrument; here, results are calculated after pooling evidence from the two flowcells.

Finally, Table 4 summarizes per-sample average coverage both genome-wide and for CpG cytosines. For Table 4, average coverage of mappable CpG cytosines (column 3) is calculated by dividing the total number of pieces of C and T nucleotide evidence aligning to CpG cytosines after filtering (column 2) by the total number of CpG dinucleotides in GRCh37 that are “mappable” in the experiment: 20.9 million. A CpG is considered “mappable” if it is overlapped by at least one non-ambiguous reference 50 bp substring; non-ambiguous substrings are those for which no other 50 bp reference substring exists within 3 mismatches.

Evidence from reads that aligned uniquely to the λ genome were used to estimate the bisulfite conversion rate for unmethylated cytosines. The conversion rate is estimated as the fraction of high-quality evidence from reads aligning uniquely to the λ phage genome that indicates lack of methylation. FIG. 2 and the final column of Table 2 show the estimates, which all lie between 99.7% and 99.8%.

FIG. 2 shows the estimated unmethylated cytosine conversion rate per sample. Bisulfite conversion efficiency is plotted per SOLiD™ 3+ flowcell. Conversion efficiency is estimated as the fraction of high-quality evidence aligning to CpG cytosines in the unmethylated λ phage genome that indicates lack of methylation.

To measure global prevalence of non-CpG cytosine methylation, all filtered nucleotide evidence overlapping non-CpG cytosine positions in the human reference genome were examined. Filtered nucleotide evidence consists of evidence (a) from reads that aligned uniquely, (b) where both overlapping colors from the original read agree with the decoded nucleotide and, (c) where nucleotides within 4 positions of either end of the alignment are excluded. For each subject the overall fraction of evidence at CpG cytosine positions where the overlapping nucleotide is a T or a C. The same was done for non-CpG cytosine positions were measured. Table 5 summarizes the results, comparing them with the rate of cytosine non-conversion estimated from filtered evidence aligning to the λ phage genome.

For Table 5, the differentially expressed genes comparing colon cancer samples to normal samples were divided into two groups: those associated hypomethyiated boundary shifts and the rest. The table shows the five categories that with statistically significant enrichment (FDRO.01). It was observed that for all subjects, the fraction of Cs observed overlapping non-CpG cytosines does not rise above the approximate fraction expected from unconverted cytosines.

FIG. 3 shows the results of a simple diagnostic to assess the basic accuracy of the filtered evidence from the reads that aligned uniquely to the human genome. FIG. 3 shows SOLiD™ 3+ read position bias in evidence for methylation. The horizontal axis represents an offset into the nucleotide alignment from the 5′ end. The vertical axis represents the fraction of filtered CpG methylation evidence from that offset that indicates that methylation is present. Only reads aligning uniquely to the GRCh37 human genome assembly are considered.

The diagnostic is to calculate, for each alignment position (1 through 48), the proportion of filtered evidence taken from that position that indicates that methylation is present. In a perfect assay, this proportion should be independent of position, and the plot should trace a straight horizontal line. In practice, the lines are not flat due to sequencing error and other noise arising from sample preparation and alignment. The relative flatness of the lines obtained for the samples gives us some assurance that the signal obtained is not substantially affected by noise such as sequencing error. Note that, while this plot shows the proportion in all alignment positions, evidence from the first four and last four positions of the alignment was additionally filtered.

After alignment, reads with a quality score of 10 or less filtered out before smoothing the methylation data (next section).

Reads were mapped to CpGs and determined if they were evidence for methylated or unmethylated cytosines. These data were then processed with the local likelihood approach, which resulted in pair-wise correlations between methylation estimates for the three normal samples of 0.97, 0.96, and 0.96.

DMR Detection.

A smoothing approach was developed, based on local likelihood estimation 33 , which leveraged the fact that the methylation status of neighboring CpG is highly correlated 34 to greatly improve precision of methylation measurements. Two related approaches; one for blocks and one for small DMRs was developed.

The alignment algorithm described in the previous section provided two counts for each CpG: number of pieces of filtered evidence indicating presence of methylation (M) and number of pieces of filtered evidence indicating unmethylation (U). The sum of these two counts was the coverage N. It was assumed that for each CpG, M followed a binomial distribution with success probability p, equal to the true methylation level, and N trials. Thus, M/N provided a naive estimate of p with standard error (M/N)×(1−M/N)/VN. Note that N, the coverage, ranged from 0 to 18,090 with a sample mean of 5.8-6.2 (after excluding CpGs with no coverage in all samples). Precision was greatly improved by leveraging the fact that proximal CpG have similar methylation levels and using a smoothing technique similar to that used in CHARM.

Because the data was binomially distributed, local likelihood estimation was used. This approach assumes that the p(L), the methylation level at genomic location L, is a smooth function of L; in other words, that CpGs that are close have similar methylation levels. The local likelihood approach uses data within windows of predefined sizes to estimates p(L) and weighing data based on distance to L (based on a tricube kernel). In addition, the binomial model ensures that data points with high coverage receive greater weight. Two window sizes were defined to detect the two different types of DMRs; for the blocks, a large window to detect low frequency differences, and a smaller window to detect high frequency differences, the small DMRs.

For each sample, the smoothed data was evaluated on the same grid of data points termed “covered CpGs”, consisting of those CpGs where at least two normals samples had a coverage of at least 2.

For the small DMRs, the high-frequency analysis used a window size of 70 CpG or 1,000 bps; whichever generated a larger region. For the blocks, the low-frequency analysis used a window size of 500 CpGs or 2,000 basepairs; again whichever was larger. Note that the use of a tricube kernel ensures that data points far from the center of the window receive a smaller weight. This approach provided highly precise estimates of CpG methylation levels p(L) for each sample. The standard errors ranged from 0-0.11 (mean of 0.04) for the high frequency smoothing and 0.01-0.04 (mean 0.02) for low frequency. An example of the results from the high-frequency smoothing is provided as FIG. 4 . FIG. 4 shows precise methylation estimates obtained by high-frequency smoothing. The circles represent the single CpG estimates of methylation, which are plotted against the CpG location. The areas of the circles are proportional to the coverage. The high-frequency smoothed values (described in detail in the Methods) are plotted as solid lines. Dashed lines represent 95% pointwise confidence intervals. The region shown in the bottom panel illustrates the statistical approach.

A method for finding differences based on t-statistics that take into account biological variability was then developed. The highly precise estimates of pi(L) for each sample i at each CpG location L was started with. The average difference between the three tumor samples and the three normal samples referred to as d(L) was obtained. To properly account for biological variability ( FIG. 5 ) the standard error of d(L) using the normal samples was estimated.

FIG. 5 illustrates the need for biological replication for detecting DMRs. FIG. 5( a ) , in the left panel, methylation patterns for three normal samples (top lines) and matched cancers (bottom) are shown. The detected DMR is shaded in pink. In the right panel, the average methylation values within the DMR for the three paired samples is shown (normal in blue, cancer in red, the matched sample pairs indicated by numbers). A t-test for the difference between normal and cancer and obtained a p value of 0.0056 was performed. FIG. 5( b ) shows the same analysis as (a) for a region in which if only normal-cancer pair 3 (thick lines) were analyzed, there would appear to be a methylation difference between cancer and normal. However, the p value when all three samples are compared is 0.52. Notice that the methylation estimates are very precise (standard error for points in right panel are <0.02). Therefore the differences seen are biological not technical, and will be seen regardless of the measurement technology.

Only the normal samples were used because as demonstrated, cancer samples are prone to high variability ( FIG. 6 ). FIG. 6 shows the increased methylation variance of common CpG sites across human cancer types. Methylation levels measured at 384 CpG sites using a custom Illumina array exhibit an increase in across-sample variability in (a) colon, (b) lung, (c) breast, (d) thyroid, and (e) kidney (Wilms tumor) cancers. Each panel shows the across-sample standard deviation of methylation level for each CpG in normal and matched cancer samples. The solid line is the identity line; CpGs above this line have greater variability in cancer. The dashed line indicates the threshold at which differences in methylation variance become significant (F-test at 99% level). In all five tissue types, the vast majority of CpGs are above the solid line, indicating that variability is larger in cancer samples than in normal. Colors indicate the location of each CpG with respect to canonical annotated CpG islands. (f) Using the CpGs that showed the largest increase in variability hierarchical clustering was performed on the normal samples. The heatmap of the methylation values for these CpGs clearly distinguishes the tissue types, indicating that these sites of increased methylation heterogeneity in cancer are tissue-specific DMRs.

The standard error se[d(L)] was therefore estimated as σ(L)*V(2/3) with σ(L) the standard deviation of the Pi(L) for the three normal samples. To improve standard error estimates, these were smoothed using a running mean with a window size of 101 observations. To avoid inflated t-statistics as a result of artificially low variance, a threshold for the standard deviation of its 75th percentile was set, before computing the smoothed result. With the standard deviation in place the t-statistic t(L)=d(L)/se[d(L)] was constructed.

For the high frequency analysis the t-statistic was further corrected for low frequency changes. This allowed discovery of local features, such as a hypermethylated small DMR, inside global features, i.e. hypo- or hypermethylated block. This correction factor was calculated by forming a fixed grid of positions 2,000 bp apart in the genome, linearly interpolating the neighboring t-statistics to obtain measurements at these positions and then smoothing this dataset with a robust smoother based on the Huber family 10 and a bandwidth of 25,000 bp. Small DMRs were then defined as contiguous CpGs within 300 bp of each other, with the t-statistics above 4.6 or below −4.6 (corresponding to the 95th quantile of the empirical distribution of the t-statistics) and all differences in the same direction. For the low frequency analysis the t-statistics cutoff was 2 and contiguous CpG were defined as within 10,000 bps from each other.

These sets of regions formed the small DMRs and blocks that were subsequently filtered and processed according to the following criteria. 1) A small DMR needed to contain at least 3 covered CpGs and have at least 1 covered CpG per 300 bp. Furthermore, the mean difference in methylation percentages between tumors and normals across the small DMR had to be greater than 0.1. 2) A block needed to be longer than 5 kb. Blocks containing CpG Islands with a mean methylation of less than 0.25 in the normal samples were separated into two. Putative blocks that were shorter than 5 kb were included as small DMRs provided they satisfied the small DMR filters above.

After filtering, pairs of small DMRs were merged if they were less than 1 kb apart, changed in the same direction (both hypermethylated or both hypomethylated), and had no covered CpGs in the area separating them. A partial list of blocks is shown as Table 6, and of small DMRs as Table 7. Additional blocks and small DMRs were determined which are not shown in Tables 6 and 7.

The data from the adenoma samples were smoothed in the same way.

DMR Classification.

Small DMRs were classified into categories based methylation profiles of the tumor and normal samples within the DMR and the two flanking regions (within 800 bp). Based on these results, the DMRs that were discovered from data exploration could be classified into three types termed loss of methylation boundaries, shifting of methylation boundaries, and novel hypomethylation ( FIG. 7 ). FIG. 7 shows loss of epigenetic stability at small DMRs. Methylation estimates are plotted against genomic location for normal samples and cancer samples. The small DMR locations are shaded. Grey bars indicate the location of blocks, CpG islands, and gene exons. Tick marks along the bottom axis indicate the location of CpGs. Pictured are examples of (a) a methylation boundary shift outward (normal-top lines in shaded area; cancer-bottom lines in shaded area), (b) a methylation boundary shift inward (cancer-top lines in shaded area; normal-bottom lines in shaded area), (c) a loss of methylation boundary (cancer-top lines in shaded area; normal-bottom lines in shaded area), and (d) a novel hypomethylation DMR (normal-top lines in shaded area; cancer-bottom lines in shaded area).

A mathematical algorithm was used to automatically classify DMRs. Briefly, mean methylation was computed for both tumor and normal samples within the DMR and in the flanking the DMR both upstream and downstream. This provided three numbers for each of the six samples. If all the normal samples showed high methylation values (>50%) in the flanking regions and low methylation values (<0.25%) and the tumor samples all showed intermediate values across DMR and flanking regions, the DMR was classified as loss of methylation boundary. If one of the flanking regions had low methylation values in both the normal and the tumor samples, the region was classified as a shift of methylation boundary. Finally, if all the normal samples showed high methylation values in the DMR and flanking regions while the tumor samples were lower in the DMR, the region was classified as novel hypomethylation. The details of the algorithm are best understood by viewing the computer code (made available upon request).

Hypomethylation in Blocks and Repeat Regions.

Repeat regions were identified based on the UCSC repeatMasker track 11 . Based on the repeats and/or blocks, the genome was segmented into regions both repeats and blocks, repeats but not blocks, not repeats but blocks, and neither repeats nor blocks. The methylation levels were computed as the average of the high-frequency smoothed methylation levels of all CpGs in the 4 different regions. Density estimates were computed from the same distribution.

Enrichment of Overlap Between Different Genomic Domains.

For each pair of different genomic domains (like blocks and LOCKs) a 2×2 table is formed containing the number of CpGs inside and outside the two genomic domains (like inside blocks and inside LOCKS, inside blocks and outside LOCKs, etc). Odds ratios and p-values were calculated using Fisher's exact text.

Copy Number Analysis.

Estimates of copy number were based on the per-base coverage obtained after alignment. The filters developed were not applied specifically for methylation measurements (described in the Bisulfite alignment Section). Note that the coverage considered here is not specific to CpGs: every genomic position is assigned a coverage value. For each sample, the average coverage in non-overlapping 10,000 bp windows was then computed, yielding two coverage vectors for each tumor-normal pair denoted cov(T) and cov(N). For each tumor-normal pair the corrected log-ratio: log 2(CN)=log 2(cov(T))−log 2(cov(N))+c was defined. Here c is a correction factor to account for different yields in each sequencing run; c is defined as the log of total sequencing yield of the normal sample divide by total yield of the tumor sample. The copy number log-ratios were segmented using circular binary segmentation (CBS) 12 . For illustrative purposes, copy number log-ratios and the associated segmentation on chromosome 20 were depicted ( FIG. 8 a ).

To determine if copy number had an effect on methylation estimates, each segment provided by CBS was divided into 100 kb regions. For each of these regions average copy number ratios as well as average methylation ratios was computed. These were then plotted ( FIG. 8 b ) and no relationship between CNV and methylation blocks was observed.

FIG. 8 shows that block and small DMR detection not affected by copy number variation. For FIG. 8( a ) for a 25 megabase region of chromosome 20 (1-25 MB), the methylation differences between all three normal-cancer pairs plotted along the chromosome, with red lines representing the average values in blocks are shown. Notice that the location of blocks is consistent across all three normal-cancer pairs. For illustrative purposed, highlighted (shade) are the seven largest blocks. For FIG. 8( b ) for each normal-cancer pair, copy number alterations is shown, quantified by log ratios (base 2) of coverage in cancer sample to coverage in normal sample, for the same region as in (a). Log-ratios of 0 are associated with lack of copy number alternations in cancer, while values larger or equal to log 2 (3/2) (dashed line) are associated with gain of copy number in cancer. The red lines represent segments obtained with the CBS algorithm (described in the Methods). Notice that each sample shows different copy number alterations. For FIG. 8( c ) differences in methylation are plotted against differences in copy number (log-ratios) for the entire genome. Specifically, for each of the segments detected by CBS, the average difference in methylation and the average log-ratio associated with copy number alteration was computed. This was done for each sample and all the points combined in one scatter-plot.

Capture Bisulfite Sequencing.

Bisulfite Padlock Probes (BSPP) were used to capture bisulfite converted gDNA using a previously described method 32 . Sequencing libraries were generated from the captured DNA, and sequenced on an Illumina™ GA II instrument, generating 79.3 million reads. Capture reads were aligned using the same methodology as the whole genome bisulfite reads, but adapted to work in nucleotide space instead of color space.

Assays for Large Organized Chromatin K9-Modifications (LOCKs).

Primary human pulmonary fibroblasts (HPF) were purchased from ScienCell Research Laboratories (San Diego, Calif.). Cell culture was conducted using the media and protocols recommended by ScienCell. Primary cells at the second passage were used for H3K9Me2 LOCK analysis. ChIP-on-chip experiments and microarray data analysis were performed as described earlier 12 .

Gene Expression Analysis.

Expression data was obtained from the gene expression barcode (rafalab.jhsph.edu/barcode). This resource combines all the expression data from the public repositories purportedly to standardize data in a way that allows one to call a gene expressed or not expressed. From this source, two independent colon cancer datasets (GSE8671 35 and GSE4183 36,37 ) were used. To define hypervariable genes an F-test using a across sample variance in tumor and normal samples computed from the original log expression was performed. A gene was defined as expressed if it had a gene expression barcode standardized value above 2.54 (p=0.01). For the fibroblast analysis datasets (GSE7890 38 , GSE11418 39 , GSE11919 40 ) were downloaded. These expression values from these datasets were also standardized using the gene expression barcode. The standardized values were used to determine if genes were expressed or not each sample.

To determine the correlation between small DMRs and expression, a gene and a small DMR were considered associated if the DMR was within 2,000 bps from the transcription start site of the gene; 6,869 genes mapped to a DMR in this way.

Data Annotation.

Annotation was obtained from the UCSC genome browser based on hgl9. In the cases where a data track was only available for hgl8 or hgl7, the UCSC liftOver tool was used to map between builds of the human genome. Specifically the repeatMasker track, the RefSeq mRNA track, and the UCSC known genes track were used.

Laminin Associated Domain (LAD) coordinates were obtained from the NKI LADs track from UCSC, generated from microarrays in fibroblast cells. PMDs were obtained from Lister et al., generated from bisulfite sequencing in fibroblast cells 11 . DNase I hypersensitive sites were obtained from the UCSC ENCODE track, using the Hies, Caco2rep1, and Caco2rep2, both the narrow and broad peak.

Results

Stochastic variation in DNA methylation as a defining feature across cancer types.

It was desired to increase the precision of DNA over the previous tiling array-based approach, termed CHARM (Comprehensive High-throughput Array for Relative Methylation) 6 , analyzing 151 colon cDMRs from Irizarry et al. 4 . A custom nucleotide-specific bead array (see Methods) was designed on the type of Illumina GoldenGate™ platform previously shown to be quantitative within 6% of DNA methylation 7 . The resulting 384 probes covered 139 regions, with 1-7 probes per region. To test the generality of DNA methylation alterations across tumor types, 290 samples were studied, including cancers from colon (10), lung (24), breast (27), thyroid (36), and kidney (Wilms') (25), with matched normal tissues to 111 of these 122 cancers, as well as 30 colon premalignant adenomas, and an additional 18 normal colon and 9 normal breast samples. To minimize the risk of genetic heterogeneity arising from sampling multiple clones DNA from small (0.5 cm×0.2 cm) sections verified by histopathologic examination were purified.

Cluster analysis of the DNA methylation values revealed that the colon cancer cDMRs largely distinguished cancer from normal across tumor types ( FIG. 9 ).

FIG. 9 shows that methylation levels of differentially methylated CpGs in colon cancer largely differentiate cancer from normal in colon(a), lung(b), breast(c), thyroid(d) and kidney(e)(Wilms) tissues. Columns and rows in each panel are ordered by a hierearchical clustering of methylation profiles using Euclidean distance. The heights of dendrogram branches, larger between tumor samples than between normal samples, illustrate the increased across-sample variability in cancer seen in a majority of CpGs in all tissues.

Even more striking than the differences in mean methylation between cancer and normal was the increased across-sample methylation variability in the cancer samples compared to normal, within each type, but occurring generally across all cancer types ( FIG. 9 ). Therefore the across-sample variance within normal and cancer samples in all five tumor/normal tissue types at each CpG site was computed. For each tissue type, the great majority of CpGs exhibited greater variance in cancer than normal ( FIG. 6( a )-( e ) ), even though these sites had been selected for differences in mean values in colon cancer. This increase was statistically significant (p<0.01, using an F-test) for 81%, 92%, 81%, 70%, and 80% of the CpG sites in colon, lung, breast, thyroid, and Wilms' tumor, respectively. Furthermore, 157 CpG sites had statistically significant increased variability in all cancer types tested. This increased stochastic variation was found in CpG islands, CpG island shores, and regions distant from islands ( FIG. 6( a )-( e ) ). This increased methylation variability in cancer are not themselves due to mean differences because the cancers show marked departure in variation from the expected binomial distribution that is followed in the normal samples. These data suggest a potential mechanism of tumor heterogeneity, namely increased stochastic variation of DNA methylation in cancers compared to normal, within each tumor type tested (see Discussion).

To rule out artifactual methylation heterogeneity caused by a greater degree of cell type heterogeneity in cancers compared to normal, a careful histological analysis for cellular subtypes within all of the colon and kidney tumor and normal samples was performed by a board-certified oncology pathologist. The normal samples were on average much more heterogeneous than the tumor samples in cellular composition. Thus, the data showing greater methylation variability in cancers cannot be due to increased cellular heterogeneity in the cancers compared to normal. This increased variation was not due to patient age, as the great majority of the samples were normal/tumor matched, and an age-corrected analysis did not change the result.

To address a potential connection between genetic and epigenetic heterogeneity that might underlie the data two experiments were performed. First, five Wilms tumor samples were selected which showed little to no intertumoral genetic variation by microarray-based copy number analysis and compared them to five colon cancers which showed relatively high amounts of copy number variation. Similar levels of methylation hypervariability were present in both types of cancers despite the differences in copy number variation. Second, immunohistochemical staining was performed for aberrant p53 expression on sections taken from Wilms tumor and colon cancers (7 samples of each), since p53 expression is a marker of chromosomal instability in both of these types of cancer. All seven colon samples tested positive for aberrant p53 whereas none of the Wilms tumors were positive. Analysis of these same samples by a custom GoldenGate™ arrays showed similar levels of methylation hypervariability in both of these cancer types, despite differential expression of p53.

These same loci which show increased variability in cancer are also able to distinguish the normal tissues from each other, but this is a mean shift rather than a variation shift, apparent from cluster analysis. Interestingly, even using a reduced set of the 25 sites which are most variable in cancer the five normal tissue types are clearly discriminated from each other ( FIG. 6( f ) ). This result reinforces the concept of a biological relationship between normal tissue differentiation and stochastic variation in cancer DNA methylation.

To determine if the increased variability is a general property of cytosine methylation in cancer or a specific property of the CpGs selected for the custom array, as a control a publicly available methylation dataset comparing colorectal cancer to matched normal mucosa on the Illumina HumanMethylation™ 27 k beadchip array was used. In this dataset it was found that only 42% of the sites showed a statistically significant increase in methylation variability, compared to 81% in the custom array (p<0.01), confirming the specificity of the cancer DMRs included in the custom array. Examining variability with respect to location relative to CpG islands, it was noted that increased stochastic variation was more common in CpGs far from islands (57%) than in shores (44%) or islands (31%), contrasting the relative representation of these locations on the 27 k array which breaks down as: distal to islands (26.4%), shores (31.6%) and islands (42%) (see Methods). This was an intriguing result as it suggested that something other than relationship to CpG islands might be defining the largest fraction of sites of altered DNA methylation in cancer.

Hypomethylation of Large DNA Methylation Blocks in Colon Cancer.

The methylation stochasticity described above appears to be a general property of cancer, affecting cDMRs in both island and non-island regions, in all five cancer types tested. To further investigate this apparent universal loss of DNA methylation pattern integrity in cancer, and to analyze the lower CpG abundance regions not previously examined by array-based methods, shotgun bisulfite genome sequencing on 3 colorectal cancers and the matched normal colonic mucosa from these three patients using the ABI SOLiD™ platform was performed. Methylation estimates with enough precision to detect differences as small as 10% methylation were desired. Because a local likelihood approach was used, which aggregated information from neighboring CpGs, and combined data from 3 biological replicates, it was determined that 4× coverage would suffice to estimate methylation values at the desired precision with a standard error of at most 3% (See Methods for details). Between 12.5 and 13.5 gigabases of data for each sample was therefore obtained, which after quality control filtering (see Methods) provided ˜5× coverage for each CpG after aligning the bisulfite converted read (Table 4). After statistical processing, pair-wise correlations were obtained between methylation estimates for the three normal samples of 0.97, 0.96, and 0.96, confirming the adequacy of sequencing coverage. Details on the mapping and statistical algorithms are available in the Methods Section.

To verify the accuracy of the methylation values obtained via the local likelihood approach, capture bisulfite sequencing was performed on the same 6 samples that were sequenced, for 39,262 regions which yielded 39.3 k-125.6 k CpG with >30× coverage. These results showed a correlation of 0.82-0.91 between the local likelihood approach and the capture bisulfite sequencing, a remarkable agreement given that the experiments were performed in different laboratories using different sequencing platforms and different protocols. Close examination of individual loci demonstrated that the methylation estimates closely track the high-coverage capture data. While capture bisulfite sequencing is far more exhaustive and CpG inclusive, traditional bisulfite pyrosequencing was further performed, which also confirmed the accuracy of the approach.

The bisulfite sequencing analysis revealed the surprising presence of large blocks of contiguous relative hypomethylation in cancer compared to normal ( FIGS. 10( a )-( b ) ). FIG. 10 shows the existence of large hypomethylated genomic blocks in human colon cancer. Shown in (a) and (b) are smoothed methylation values from bisulfite sequencing data for cancer samples (red) and normal samples (blue) in two genomic regions. The hypomethylated blocks are shown with pink shading. Grey bars indicate the location of PMDs, LOCKs, LADs, CpG Islands, and gene exons. Note that the blocks coincide with the PMD, LOCKS, and LADs in panel (a) but not in (b). (c) The distribution of methylation values for the normal samples (blue) versus the cancer samples (red) demonstrates global hypomethylation of cancer compared to normal. (d) The distribution of methylation values in the blocks (solid lines) and outside the block (dashed lines) for normal samples (blue) and cancer samples (red). Note that while the normal and cancer distributions are similar outside the blocks, within the blocks methylation values for cancer exhibit a general shift. (e) Distribution of methylation differences between cancer and normal samples stratified by inclusion in repetitive DNA and blocks. Inside the blocks, the average difference was ˜20% in both in repeat and non-repeat areas. Outside the blocks, the average difference was ˜0% in repeat and non-repeat areas, indicating that blocks rather than repeats account for the observed differences in DNA methylation.

13,540 such regions of median size 39 kb, varying from 5 kb-10 MB were identified (partial listing in Tables 6 and 8). The across-cancer average loss of mean DNA methylation throughout all the blocks was 12%-23% (partial listing shown in Table 6). Remarkably, these hypomethylated blocks in cancer corresponded to more than half of the genome, even when accounting for the number of CpG sites within the blocks (Table 8). Analysis shows that these blocks are common across all three cancers; an analysis of the tumors individually versus a normal profile shows consistent block boundary locations (see FIG. 10 , FIG. 2 , and Methods). The existence of a small fraction of hypermethylated blocks (3% of all differentially methylated blocks) in cancer was also noted, with examples visible on the far right edge of FIG. 10( a ) and FIG. 10( b ) (Table 8). A histogram of smoothed methylation values shows the shift in distribution of global DNA methylation ( FIG. 10( c ) ). The relatively small compartments of fully unmethylated or fully methylated changed little in the aggregate between normal and cancer. Instead, the predominant change in block methylation in cancer was a loss in the more abundant compartment of intermediate methylation levels (mean 73% for all samples) to a significantly less methylated level (50-61%). Examining methylation density distributions stratified by block and non-block regions supported this finding ( FIG. 10( d ) ). It is important to note that copy number variation does not drive the results presented here. First, the locations of copy number variation changes in cancer, observed in the data, were not consistent across subjects, in contrast to the consistent locations of methylation changes in the blocks identified above. Second, the methylation difference estimates provided by the statistical approach did not correlate with copy number values.

Global quantitative hypomethylation in cancer 8 is widely attributed to the presence of normally methylated repetitive elements 9 and may be relevant to colon cancer as UNE-1 element hypomethylation is associated with worse prognosis in colon cancer 10 . It was observed that in normal tissues, repetitive elements were more methylated than non-repetitive regions (76% vs. 66%). To determine whether such repetitive elements were responsible for the block hypomethylation, differences in methylation levels inside and outside repeat elements (see Methods), both inside and outside the blocks were compared. It was clear from this analysis that most of the global hypomethylation was due to hypomethylated blocks ( FIG. 10( e ) ) and not to the presence of repetitive elements. As repetitive elements are slightly enriched in blocks (odds ratio 1.4), much of the apparent repeat-associated methylation may in fact be due to blocks. This result does not exclude repeat-associated hypomethylation, since not all repeats were mappable. However, 57% of L1 elements, 94% of L2 elements, 95% of MIR sequences, and 18% of Alu elements were covered by the data ( FIG. 11 ) and did not show repeat-specific hypomethylation ( FIG. 4 ). It is important to note it is possible that Alu sequences not covered by the data are somehow more hypomethylated than covered Alu sequences, and thus would contribute to global hypomethylation.

FIG. 11 shows the methylation changes of selected repetitive DNA families. Distribution of methylation differences between cancer and normal samples stratified by repeat family and inclusion in blocks. CpGs outside all repetitive elements were used as controls. Each graph is the percentage of mappable CpGs in each repeat family.

Lister et al. recently performed bisulfite sequencing analysis of the H1 human embryonic stem cell line compared to the IMR90 fibroblast line, identifying large regions of the genome that are less methylated in fibroblast cells than ES cells, referred to as partially methylated domains (PMDs) 11 . The intermediate-methylation level regions identified above largely coincided with the PMDs, containing 85% of CpGs inside PMDs (odds ratio 6.5, R<2×10 −16 , Table 9). Large organized chromatin lysine (K) modifications, or LOCKs, genome-wide in normal mouse cells, that are associated with both constitutive and tissue-specific gene silencing were previously described that 12 . Therefore LOCKs in primary human cells (see Methods) mapped and whether the LOCKs also correspond to the hypomethylated blocks was determined. Remarkably, 89% of the LOCKs were contained within the blocks (odds ratio 6.8, P<2×10 −16 ). LOCKs are also known to overlap with nuclear lamina-associated domains or LADs 12 . Approximately 83% of the LADs were also contained within the blocks, similar to the LOCKs (odds ratio 4.9, P<2×10 −16 ). In addition, DNase I hypersensitive sites, a structural signal for regulatory regions 13 were enriched within 1 kb of block boundaries and small DMRs (p<2×10 −16 for both). Thus the large hypomethylated blocks identified in cancer correspond to a genomic organization identified in normal cells by several complementary methods. Note that although the PMDs and the hypomethylated blocks largely overlap, significant differences in gene expression in cancer between non-overlapping blocks and PMDs is demonstrated later.

Furthermore, a relationship between the hypervariable CpGs across all cancer types, from the custom Illumina™ array described above, and the hypomethylated blocks identified by whole genome bisulfite sequencing was observed. Because the custom array was designed to contain both hypomethylated and hypermethylated CpG sites tested for enrichment of these hypervariable 157 CpG sites in hypomethylated blocks and in the small number of hypermethylated blocks. It was found that 63% of the hypomethylated hypervariable CpGs were indeed within hypomethylated blocks, and 37% of the hypermethylated hypervariable CpGs were within the hypermethylated blocks. To test for statistical significance, these percentages were compared to those found for the CpGs in the HumanMethylation27K array CpGs that did not vary differently in cancer and normal. These percentages were 13% and 1.5% for the hypomethylated and hypermethylated blocks, respectively, demonstrating high statistical significance for enrichment of hypervariably methylated CpGs in blocks (p<2×10 −16 ; Table 10).

Small DMRs in Cancer Involve Loss of Stability of DNA Methylation Boundaries.

A statistical algorithm (see Methods) was then developed for detecting consistent DNA methylation changes in regions smaller than the blocks (≤5 kb). Critical for this algorithm was the analysis of biological replicates, as it was found that regions showing across-subject variability in normal samples would otherwise be easily confused with DMRs if only one cancer-normal pair were available ( FIG. 5 ). Methylation measurements in these smaller regions exhibited good agreement with methylation measurements in the previous CHARM-based microarray analysis 4 ( FIG. 12 ).

FIG. 12 shows a comparison of CHARM microarray data and bisulfite sequencing for measuring methylation in (a) normal and (b) cancer samples. Average methylation level from previously published CHARM microarray data (Irizarry et al. (2009)) (y-axis) is plotted versus the average methylation obtained from high-frequency smoothed bisulfite sequencing data. Each point represents one of the cDMR regions originally identified in (Irizarry et al. (2009)). Note the high degree of correlation between CHARM and sequencing.

These are referred to as small DMRs to distinguish them from the large (>5 kb) differentially methylated blocks described above. The increased comprehensiveness of sequencing over CHARM and other published array-based analyses allowed us to detect many more small DMRs than previously reported, 5,810 hypermethylated and 4,315 hypomethylated small DMRs (partial listing of DMRs shown in Table 7). Previous findings 4 that hypermethylated cDMRs are enriched in CpG islands while hypomethylated cDMRs are enriched in CpG island shores were also confirmed (Table 8).

The increased comprehensiveness of sequencing also allowed us to quantify the different types of methylation changes at CpG Islands specifically (Table 11). It was found that: (1) 61% of islands were unmethylated in normal colon, and of these 17% showed a >10% absolute increase in methylation in cancer; (2) 20% of islands were found to be fully methylated in the normal colon, and of these 24% showed a >10% decrease in absolute methylation in cancer; and (3) 18% of islands show partial methylation in normal colon, and of these 17% showed increased and 36% decreased methylation in cancer. Of the 3907 hypermethylated small DMRs overlapping CpG Islands, 63.7% involved changes in either shift or loss boundary. Additionally, to determine if the same pattern held true for promoter, genic, intergenic, and repeat CpG islands analysis in Table 11 was stratified. It was found that in the normal samples promoter CpG Islands are largely unmethylated but that the findings above still hold. In contrast a balanced proportion of methylated and unmethylated (35% versus 26%) intergenic CpG islands with 45% hypomethylated in cancer was observed (Table 1).

When the architecture of the small DMRs was examined, the most common characteristic was a shift in one or both of the DNA methylation boundaries of a CpG island out of the island into the adjacent region ( FIG. 7( a ) ) or into the interior of the island ( FIG. 7( b ) ). Where the boundary shifts are observed, the same regions would appear on array-based data as hypermethylated islands, if the boundary is shifting into the island, or to hypomethylated shores, if the boundary is shifting outside the island.

The second most frequent category of small DMRs involved loss of methylation boundaries at CpG islands. For example, the architecture of many hypermethylated cDMRs were defined in normal samples by unmethylated regions (<10% methylation) surrounded by highly methylated regions (75-95% methylation), with boundaries ˜100 bp wide. In cancer, these regions exhibited an approximately constant methylation level of approximately 40-60% throughout (Table 8, FIG. 7( c ) ). These regions with loss of methylation boundaries largely correspond to what are typically classified as hypermethylated islands in cancer.

Hypomethylated cDMRs were also found that arose de novo in highly methylated regions outside of blocks. In these cases, regions in which normal colon tissue was 75-95% methylated dropped to lower levels (20-40%) in cancer ( FIG. 7( d ) ). These regions are referred to as novel hypomethylated DMRs, and they usually corresponded to CpG-rich regions that were not conventional islands (Table 8). In summary, in addition to the hypomethylated blocks, 10,125 small DMRs were found, 5,494 of which clearly fell in three categories: shifts of methylation boundaries, loss of methylation boundaries, and novel hypomethylation. Note that not all small DMRs followed a consistent pattern across all three sample pairs and therefore did not permit their classification (Table 8).

Methylation-Based Euclidean Distances Show Colon Adenomas Intermediate Between Normals and Cancers.

By performing a multidimensional scaling of the methylation values measured via the custom array in each colon sample, it was observed that the normal samples clustered tightly together in contrast to the dispersed methylation values measured in the cancer samples ( FIG. 13( a ) ). This is consistent with the observed increase in methylation variability found in cancer that was described earlier. 30 colon adenomas were also analyzed for methylation on the custom array, and found that the adenomas were intermediate in both variability within samples and distance to the cluster of normal samples ( FIG. 13( a ) ).

FIG. 13 shows that adenomas show intermediate methylation variability. (a) Multidimensional scaling of pairwise distances derived from methylation levels assayed on a custom array. Note that cancer samples (red) are largely far from the tight cluster of normal samples (blue), while adenoma samples (black) exhibit a range of distances: some are as close as other normal samples, others are as far as cancer samples, and many are at intermediate distances. (b) Multidimensional scaling of pairwise distances derived from average methylation values in blocks identified via bisulfite sequencing. Matching sequenced adenoma samples (labeled 1 and 2) appear in the same locations relative to the cluster of normal samples in both (a) and (b). (c) Methylation values for normal (blue), cancer (red) and two adenoma samples (black). Adenoma 1, which appeared closer to normal samples in the multidimensional scaling analysis (a), follows a similar methylation pattern to the normal samples. However, in some regions (shaded with pink) differences between Adenoma 1 and the normal samples are observed. Adenoma 2 shows a similar pattern to cancers.

Whole genome bisulfite sequencing was then performed on two of these adenomas, a premalignant colon adenoma with relatively small methylation-based distance to the normal colons and an adenoma with a large methylation-based distance to the normal colons, similar to the cancer samples. For each adenoma, the average methylation levels in each block and used these values to compute pairwise Euclidean distances were then computed. These results, based on measurements from the hypomethylated blocks, confirm the characteristic described above for the bead array data: genome-wide increased variability in cancers compared to normals with adenomas exhibiting intermediate values ( FIG. 13 ).

Expression of Cell Cycle Genes Associated with Hypomethylated Shores in Cancer.

Whole genome methylation analysis has shown an inverse relationship between gene expression and methylation, especially at transcriptional start sites 14 . To study the relationship between methylation and gene expression in small DMRs, public microarray gene expression data was obtained from cancer and normal colon samples (see Methods Section) and compared to results from the sequencing data. Each gene represented in the expression microarrays was mapped to the closest DMR. A gene and DMR were considered associated if the DMR was within 2 kb of the transcription start site of the gene; 6,869 genes mapped to a DMR in this way. The expected inverse relationship between DNA methylation and gene expression was observed (r=−0.27, p<2×10 −16 , FIG. 14 ).

FIG. 14 shows an inverse correlation of gene expression with methylation at small DMRs. Average log gene expression values, obtained from GEO dataset GSE8671, plotted versus the average difference in methylation of a nearby small DMR. A gene and a small DMR were considered associated if the DMR was within 2,000 bp of the transcription start site of the gene; 6,869 genes mapped to a DMR in this way. Different types of small DMRs are indicated by color, boundary shift (green), boundary loss (orange), novel hypomethylation (purple) and other (pink). The dashed lines represent a fold change of 2 in the gene expression comparison.

By examining the inverse relationship between methylation and gene expression for each category of small DMRs separately, the strongest relationship for hypomethylated shores due to methylation boundary shifts (Table 12) was noticed. Gene ontology enrichment analysis was performed for differentially expressed genes (FDR<0.05), comparing those associated with hypomethylated boundary shifts to the other categories. These categories (Table 13) were strongly enriched for cell-cycle related genes CEP55, CCNB1, CDCA2, PRC1, CDC2, FBXO5, AURKA, CDK1, CDKN3, CDK7, and CDC20B, among others (Table 14).

Increased Variation in Gene Expression in Hypomethylated Blocks and DMRs.

Given the increased DNA methylation variability in cancer found with the custom bead arrays ( FIG. 6( a )-( e ) ), the methylation patterns in the blocks was examined and found increased variability in cancer as well. Specifically, across-subject methylation variability levels within the blocks in cancer and normal were compared, and found a striking similarity to the results found with the Illumina custom bead array as described above ( FIG. 6( a )-( e ) compared to FIG. 15 ).

FIG. 15 shows increased variation in methylation between normal and cancer samples in blocks. Across-sample standard deviation of methylation level for each block of normal versus cancer samples. Average methylation levels were computed for each block using high frequency smoothed SOLiD™ bisulfite sequencing data. The solid line is the identity line; CpGs above this line have greater variability in cancer. As in FIG. 6 , the vast majority of blocks show an increased variation in cancer compared to normal samples.

To study the relationship to gene expression in colon cancer, public gene expression data was obtained from cancer and normal samples (see Methods Section). The expression data was processed with the gene expression barcode 16 , a statistical algorithm that leverages information from thousands of publicly available datasets to standardized expression measurements in a way that expression values for silenced genes follow a standard normal distribution (mean 0 and standard deviation 1). Note that while expression values reported in the public databases cannot be used to call genes silenced or expressed, the barcode approach permits the construction of statistical tests for silencing. Genes in the blocks were generally silenced (80% genes silenced in all samples) both in normal and cancer samples. Of the genes consistently transcribed in normal tissue, albeit at low levels, 36% are silenced in blocks in cancers, compared to 15% expected by chance. This is consistent with other reports in the literature, e.g. Frigola et al 17 .

More strikingly than subtle differences in gene silencing, substantial enrichment of genes exhibiting increased expression variability in cancer compared to normal samples in the hypomethylated blocks was found. First, the possibility that this observed increased variability was due to tissue-specific genes being over-represented in blocks and the potential high cellular heterogeneity of cancer by noting the lack of association between tissue-specific genes and cancer variability outside of the blocks was ruled out. A clear association between increased variability in expression of a gene and its location within a hypomethylated block was then noticed. For example, 26 of the 50 genes exhibiting the largest increase in expression variability were inside the blocks; 52% compared to the 17% expected by chance (p=3×10 −9 ); and there is significant enrichment at P=0.01 up to as many as 1,000 genes of the 4,940 genes on the array and within the blocks. Expression levels for 25 of these 26 most hypervariable genes exhibited an interesting pattern: while never expressed in normal samples, they exhibited stochastic expression in cancer ( FIGS. 16 and 17 ). For example the genes MMP3, MMP7, MMP10, SIM2, CHI3L1, STC1, and WISP (described in the Discussion) were expressed in 96%, 100%, 67%, 8%, 79%, 50%, and 17% of the cancer samples, respectively, but never expressed in normal samples (Table 14).

FIG. 16 shows high variability of gene expression associated with blocks. (a) An example of hypervariably expressed genes contained within a block, MMP7, MMP10, and MMP3 (red). Methylation values for cancer samples (red) and normal samples (blue) with hypomethylated block locations highlighted (pink shading) are plotted against genomic location. Grey bars are as in FIG. 8 . (b) Standardized log expression values for 26 hypervariable genes in cancer located within hypomethylated block regions (normal samples in blue, cancer samples in red). Standardization was performed using the gene expression barcode. Genes with standardized expression values below 2.54, or the 99.5th percentile of a normal distribution (horizontal dashed line) are determined to be silenced by the barcode method 15 . Vertical dashed lines separate the values for the different genes. Note there is consistent expression silencing in normal samples compared to hypervariable expression in cancer samples.

FIG. 17 shows hypervariable gene expression in cancer in hypomethylated blocks. Standardized (using gene expression barcode) expression values for the 26 of the 50 most hypervariable genes in cancer which are within hypomethylated block regions. Genes with standardized expression values below 2.54 (dotted horizontal line) are determined to be silenced by the barcode method (Zilliox and Irizarry, 2007). Expression values for each gene, separated by vertical dotted lines, from dataset GSE4183 are plotted for normal (blue) and cancer (red) samples. Note there is consistent expression silencing in normal samples compared to variable expression in cancer samples.

Functional Differences Between Hypomethylated Blocks and PMDs.

As noted above, the hypomethylated blocks observed, substantially overlapped the PMDs previously reported in a fibroblast cell line by Lister et al. 11 . The genomic regions that did not overlap between blocks and PMDs were therefore examined, in order to identify potential functional differences between them. They were grouped into two sets: 1) regions within the hypomethylated blocks not in the PMDs (B+P−) and 2) the counterpart regions within the PMDs but not in the hypomethylated blocks (B−P+). To examine the expression of genes within these regions, public microarray gene expression data from fibroblast samples was obtained (see Methods Section). As expected, the genes in the fibroblast PMDs were relatively silenced in the fibroblast samples; 79% of genes were silenced in all samples within the PMD compared to 49% outside, p<2×10 −16 ). Furthermore, genes that were silenced in fibroblast samples and consistently expressed in colon were enriched in the B−P+ regions (odds ratio of 3.2, p<2×10 −16 ), while genes consistently silenced in colon and consistently expressed in fibroblast samples were enriched in the B+P− regions (odds ratio 2.8, p=0.0004). B+P− and B−P+ genes were then examined with respect to expression variability in cancer, focusing on the 50 hypervariable genes described above. These genes were markedly enriched in the B+P− regions (p=0.00013), yet showed no enrichment in the B−P+ regions. These results suggest that hypervariable gene expression in colon cancer is related specifically to their presence in hypomethylated blocks.

Discussion

In summary, it is shown that colon cancer cDMRs are generally involved in the common solid tumors of adulthood, lung, breast, thyroid, and colon cancer, and the most common solid tumor of childhood, Wilms tumor, with tight clustering of methylation levels in normal tissues, and marked stochastic variation in cancers. Current efforts to exploit DNA methylation for cancer screening are focused on identifying a cancer-specific profile 18 , but the data shown here suggest that epigenetic variability per se may be pathognomonic for cancer. Future efforts might be directed at defining the cancer epigenome as the departure from a narrowly defined normal profile rather identifying a specific narrowly defined cancer profile.

Surprisingly, two-thirds of all methylation changes in colon cancer involve hypomethylation of large blocks, with consistent locations across samples, comprising more than half of the genome. The most variably expressed genes in cancer are extraordinarily enriched in the blocks, and they involve genes associated with tumor heterogeneity and progression, including three matrix metalloproteinase genes, MMP3, MMP7, and MMP10 19 , and a fourth, SIM2, which acts through metalloproteinases to promote tumor invasion 20 . Another, STC1, helps mediate the Warburg effect of reprogramming tumor metabolism 21 . CHI3L1 encodes a secreted glycoprotein associated with inflammatory responses and poor prognosis in multiple tumor types including the colon 22 . WISP genes are targets of Wnt-1 and are thought to contribute to tissue invasion in breast and colon cancer 23 . Gene ontology enrichment analysis 15 was performed for genes associated with hypervariable expression in blocks (FDR<0.05), which showed substantial enrichment for categories including extracellular matrix remodeling genes (Table 13). One cautionary note raised by these findings is that treatment of cancer patients with nonspecific DNA methylation inhibitors could have unintended consequences in the activation of tumor-promoting genes in hypomethylated blocks. It is also important to note that previous studies 24,25 have shown large-region hypermethylation or no regional methylation change, although the present study is based on whole-genome bisulfite sequencing.

The small DMRs, which while representing a relatively small fraction of the genome (0.3%), are numerous (10,125), and frequently involve loss of boundaries of DNA methylation at the edge of CpG islands, shifting of the DNA methylation boundaries, or the creation of novel hypomethylated regions in CG-dense regions that are not canonical islands. These data underscore the importance of hypomethylated CpG island shores in cancer. In particular, the large number of shores associated with hypomethylation and gene overexpression in cancer are enriched for cell cycle related genes, suggesting a role in the unregulated growth that characterizes cancer.

A model it proposed that might help to explain the relationship between tissue-specific DMRs and the sites of extreme methylation variability in cancer. Pluripotency itself might require a certain degree of stochastic gene expression at some loci, allowing for differentiation along alternative pathways in response to external stimuli, or even in part randomly. The epigenome could collaborate by changing its physical configuration to relax the stringency of epigenetic marks away from extremes near 100% and 0%, since variance increases away from the extremes. One potential way the epigenome may do this is by reducing the methylation of large blocks in cancer, which overlap with known functionally important domains, LOCKs and LADs, and thus could involve a change in the chromatin packing density or proximity to the nuclear lamina. Similarly, subtle shifts in DNA methylation boundaries near CpG islands may drive chromatin organization and normal tissue-specific and aberrant cancer-specific gene expression. Consistent with this idea, a shift in methylation boundary was shown to differentially regulate genes escaping X-inactivation 26 . It will therefore be important in the future to relate the data presented here to large-scale chromatin organization, and to normal developmental changes in DNA methylation and chromatin, as well as to expand future epigenetic investigations on the boundaries of blocks and CpG island shores.

It is not argued that epigenetic changes in cancer arise independently of genetic changes; both are clearly important and they are highly inter-related. However, the increased methylation and expression variability in each cancer type, compared to the normal samples, is consistent with the potential selective value of increased epigenetic plasticity in a varying environment first suggested for evolution but equally applicable to the strong but variable selective forces under which a cancer grows, such as varying oxygen tension or metastasis to a distant site 27 . Thus, increased epigenetic heterogeneity in cancer at cDMRs (which are shown also as tDMRs) could underlie the ability of cancer cells to adapt rapidly to changing environments, such as increased oxygen with neovascularization, then decreased oxygen with necrosis; or metastasis to a new intercellular milieu.