DMRforPairs: identifying Differentially Methylated Regions between unique samples using array based methylation profiles

Background Array based methylation profiling is a cost-effective solution to study the association between genome methylation and human disease & development. Available tools to analyze the Illumina Infinium HumanMethylation450 BeadChip focus on comparing methylation levels per locus. Other tools combine multiple probes into a range, identifying differential methylated regions (DMRs). These tools all require groups of samples to compare. However, comparison of unique, individual samples is essential in situations where larger sample sizes are not possible. Results DMRforPairs was designed to compare regional methylation status between unique samples. It identifies probe dense genomic regions and quantifies/tests their (difference in) methylation level between the samples. As a proof of concept, DMRforPairs is applied to public data from four human cell lines: two lymphoblastoid cell lines from healthy individuals and the cancer cell lines A431 and MCF7 (including 2 technical replicates each). DMRforPairs identified an increasing number of DMRs related to the sample phenotype when biological similarity of the samples decreased. DMRs identified by DMRforPairs were related to the biological origin of the cell lines. Conclusion To our knowledge, DMRforPairs is the first tool to identify and visualize relevant and significant differentially methylated regions between unique samples.


Background
Epigenetic (de)regulation, including DNA (CpG) methylation, is associated with development, differentiation and many human diseases [1][2][3] including the initiation and progression of various cancers [3][4][5][6][7][8]. While the primary DNA sequence is mostly stable during the lifetime of an individual, the epigenome is highly dynamic and responsive. Because of this, it provides valuable information about (past) (micro-)environmental conditions in the context of human disease and development [9,10].
DNA CpG methylation is routinely investigated on a genome wide scale [2,3]. The methylation profile can be assessed using micro-arrays or sequencing by applying (1) methylation-sensitive restriction enzymes or immunoprecipitation (anti-5mC) or (2) bisulfite-based treatment, which converts unmethylated cytosines into uracils [11]. The Illumina Infinium HumanMethylation450 BeadChip (450 K) is a bisulfite-based, cost-effective, two-color array querying over 480,000 independent genomic positions (99% Refseq genes, 96% CpG islands) [12][13][14]. Various tools are available to pre-process and analyze the 450 K data, but differential methylation is primarily detected per locus or by comparing differential patterns across regions using groups of samples [15]. The latter is implemented in IMA and bumphunter. Indeed, IMA offers region based analysis [16], but it does not work when using unique samples. Bumphunter identifies regional changes in the regression coefficient between methylation status and phenotype. Therefore, bumphunter (like IMA) requires groups of samples of sufficient size to estimate this coefficient for each probe [17]. However, when analyzing small numbers of samples with unique characteristics (e.g. normal and affected tissue of a clinically unique patient or a manipulated cell model), large series of samples are not available and current methods cannot be applied. Although larger series of samples are preferred (biological replicates or more patients), comparison of unique samples is desired in such a situation. DMRforPairs was designed to address this problem by comparing regional methylation status between unique samples.

Implementation
The algorithm consists of a number of phases ( Figure 1A) with fully customizable parameters which will be discussed below: 1) Recoding of the probe classes 2) Identification of regions with sufficient probe density 3) Quantification and testing of (difference in) methylation status.
Data import and pre-processing As input DMRforPairs requires the methylation percentage of each CpG site in each sample. It was originally designed for the 450 K array, but is applicable to any platform that generates a methylation percentage per CpG site and has sufficient coverage. For example, Additional file 1 illustrates the algorithm's applicability to data generated using Nimblegen microarrays and the McrBC protocol (CHARM). DMRforPairs does not provide functions to import, filter (cross-hybrization, SNPs in probe sequencing) or preprocess 450 K data because of the existence of a number of excellent, well maintained pre-processing R-pipelines [11,15,16,[18][19][20][21][22]  pipelines. The output of these pipelines serves as input for DMRforPairs.

Recoding of the probe classes
Illumina assigns the majority of probes to eleven specific classes according to their association to one or more functional regions (relation to gene: Body, 5'UTR, 3'UTR, 1 st exon, TSS1500, TSS200; relation to CpG island: Island, Northern/Southern Shelf & Shore [12]). Highly detailed classification may result in too low probe density per class as DMRforPairs investigates probes in close proximity to each other within each class individually. DMRforPairs therefore allows custom grouping and/or selection of classes. Three commonly used schemes are hard-coded in the software: (0) retain all 11 classes, (1) group on relation to gene/transcription start site/CpG island or (2) put all probes in one class. The last option ignores the assigned classes as it might be desirable to just let DMRforPairs identify DMRs without providing information about probes that belong to the same functional class. This option can also be used in case this functional classification is unknown.

Identification of regions with sufficient probe density
A region of interest meets the following criteria: Default settings of d min are based on decreasing correlation between methylation status of adjacent loci when evaluated at inter-locus distances between 0 and 1 kb (200 bp is reported to correlate well) [11,23]. The default value for n min is based on the theoretical minimal number of 2×4 observations required for statistical testing using Mann-Whitney U test. Probes annotated to more than one class are included in multiple regions and fully identical regions from different classes are merged into one region with a combined class. Figure 2 illustrates the number of regions identified for various settings of d min and n min and the fraction of all probes included in the regions. A function is available in DMRforPairs to generate these benchmarking results for specific data sets and tune the settings of n min and d min .

Quantification and testing of methylation status
As recommended, the methylation percentage β and the M-values (logit 2 (β)) were used for visualization and statistical computations respectively [24]. Descriptive statistics are computed by DMRforPairs for all regions and samples (optional parallelization). These consist of median methylation levels (M and β values) and pairwise differences in median methylation level between all samples. If the median difference in M value between any pair of samples is sufficiently large in a specific region (> |ΔM|), the difference is formally tested using the Mann Whitney U or Kruskal-Wallis test. Pairwise testing is performed if indicated (n > 2 & p Kruskal-Wallis ≤ 0.05). An α of 0.05 after adjustment for multiple testing (Bejamini & Hochberg (FDR) [25]) is used to select significant regions (default settings). α and the method to correct for multiple testing can be specified by the user.
Several issues need to be kept in mind when choosing the algorithm's parameters and interpreting (test) results. In general, setting the algorithm's parameter more stringently (|ΔM|↑,n min ↑,d min ↓)) reduces the amount of regions to be tested, but also discards potential DMRs that are less optimally covered by the probes on the array. Concerning the |ΔM| threshold it is important to be  [24]. Also, correlation of methylation levels of CpG sites located closely together on the genome should be kept in mind. The potential presence of correlation warrants careful evaluation of statistical test results related to the independency assumption even though methylation levels at specific sites are technically (different probes) and biologically (different genomic positions) independent. Finally, comparisons with a higher number of probes per region have a higher power and are more likely to survive multiple testing. Therefore, the list of significant DMRs is theoretically biased towards regions with more probes (i. e. larger sample size). This bias was limited in a comparison of samples which are derived from a strongly biologically different origin ( Figure 1B,C). When comparing the more similar samples there was some overrepresentation of regions with a high number of probes (28 DMRs, Figure 1B,D).

Visualization, export and exploration
HTML tables listing all, only relevant (median difference ≥ ΔM) and only significant regions are generated with links to genome browsers ( Figure 3A, application of the R2HTML package [26]). Links are also provided to images depicting the observed methylation profiles and a text file with additional descriptive statistics (Figure 3). Pairwise plots are generated in case of more than two biological samples. For relevant and significant DMRs an extended output can be generated including thumbnails in the HTML tables and visualizations that also depict transcripts annotated (close) to the region ( Figure 3C, application of the Gviz & GenomicRanges packages [27,28]). In addition, DMRforPairs includes a number of functions to further inquire the data. Methylation status of genes of interest, regions identified by DMRforPairs and custom genomic intervals can be visualized, annotated and quantified/tested.

Results
In the Illumina manifest, 12% of the probes were not assigned to any of the 11 categories (discarded in this analysis with recode parameter set to 1). 35% of the remaining probes was included in one or more regions, leading to 29,404 potential regions of interest. Samples were compared pairwise in descending order of biological similarity: technical duplicates, lymphoblastoid cell lines and cancer cell lines (average of duplicates) (Figure 1, Additional file 2).
Between A431 and MCF7 2,626 DMRs were identified ( Figure 1B,C). On top of the list was FAM195A a gene with known low expression [33] and complete methylation in MCF7. In A431, the region showed complete demethylation, but no public expression data was available for this cell line. The rest of the top-5 consisted of homeobox genes which are frequently methylated in breast cancer and active in squamous cell carcinoma [34,35]. Cancer was by far the strongest overrepresented disease category in the enrichment/network analysis (989 genes, p = 1.31E-19 -2.71E-4). Enriched subcategories included breast cancer (n = 234, p = 2.06E-10), head and neck (squamous cell) carcinoma (n = 131, p = 1.30E-7) and genital tumor (n = 192, p = 1.94E-7).

Conclusions
DMRforPairs defines genomic regions using local probe density and optionally functional homogeneity. It quantifies, tests, annotates and visualizes (differential) methylation patterns between unique samples including pairwise comparison of samples if n > 2. Here, it is shown that in two lymphoblastoid cell lines from healthy individuals and cancer cell lines A431 and MCF7 (including 2 technical replicates each), DMRforPairs was able to identify an increasing number of DMRs related to the sample phenotype when biological similarity of the samples decreased. DMRs identified by DMRforPairs were related to the biological origin of the cell lines. In addition, DMRforPairs has been applied successfully in the analysis of integrated genome-wide epigenetic and expression profiles of germ cell cancer cell lines [36].