Skip to main content

methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles

Abstract

DNA methylation is a chemical modification of cytosine bases that is pivotal for gene regulation,cellular specification and cancer development. Here, we describe an R package, methylKit, thatrapidly analyzes genome-wide cytosine epigenetic profiles from high-throughput methylation andhydroxymethylation sequencing experiments. methylKit includes functions for clustering, samplequality visualization, differential methylation analysis and annotation features, thus automatingand simplifying many of the steps for discerning statistically significant bases or regions of DNAmethylation. Finally, we demonstrate methylKit on breast cancer data, in which we find statisticallysignificant regions of differential methylation and stratify tumor subtypes. methylKit is availableat http://code.google.com/p/methylkit.

Rationale

DNA methylation is a critical epigenetic modification that guides development, cellulardifferentiation and the manifestation of some cancers [1, 2]. Specifically, cytosine methylation is a widespread modification in the genome, and itmost often occurs in CpG dinucleotides, although non-CpG cytosines are also methylated in certaintissues such as embryonic stem cells [3]. DNA methylation is one of the many epigenetic control mechanisms associated with generegulation. Specifically, cytosine methylation can directly hinder binding of transcription factorsand methylated bases can also be bound by methyl-binding-domain proteins that recruitchromatin-remodeling factors [4, 5]. In addition, aberrant DNA methylation patterns have been observed in many humanmalignancies and can also be used to define the severity of leukemia subtypes [6]. In malignant tissues, DNA is either hypo-methylated or hyper-methylated compared to thenormal tissue. The location of hyper- and hypo-methylated sites gives distinct signatures withinmany diseases [7]. Often, hypomethylation is associated with gene activation and hypermethylation isassociated with gene repression, although there are many exceptions to this trend [7]. DNA methylation is also involved in genomic imprinting, where the methylation state of agene is inherited from the parents, but de novo methylation also can occur in the earlystages of development [8, 9].

A common technique for measuring DNA methylation is bisulfite sequencing, which has the advantageof providing single-base, quantitative cytosine methylation levels. In this technique, DNA istreated with sodium bisulfite, which deaminates cytosine residues to uracil, but leaves5-methylcytosine residues unaffected. Single-base resolution, %methylation levels are thencalculated by counting the ratio of C/(C+T) at each base. There are multiple techniques thatleverage high-throughput bisulfite sequencing such as: reduced representation bisulfite sequencing (RRBS)[10] and its variants [11], whole-genome shotgun bisulfite sequencing (BS-seq) [12], methylC-Seq [13], and target capture bisulfite sequencing [14]. In addition, 5-hydroxymethylcytosine (5hmC) levels can be measured through amodification of bisulfite sequencing techniques [15].

Yet, as bisulfite sequencing techniques have expanded, there are few computational toolsavailable to analyze the data. Moreover, there is a need for an end-to-end analysis package withcomprehensive features and ease of use. To address this, we have created methylKit, amulti-threaded R package that can rapidly analyze and characterize data from many methylationexperiments at once. methylKit can read DNA methylation information from a text file andalso from alignment files (for example, SAM files) and carry out operations such as differentialmethylation analysis, sample clustering and annotation, and visualization of DNA methylation events(See Figure 1 for a diagram of possible operations). methylKit hasopen-source code and is available at [16] and as Additional file 1 (see also Additional file 2 for the user guide and Additional file 3 for thepackage documentation ). Our data framework is also extensible to emerging methods in quantizationof other base modifications, such as 5hmC [14], or sites discovered through single molecule sequencing [17, 18]. For clarity, we describe only examples with DNA methylation data.

Figure 1
figure 1

Flowchart of possible operations by methylKit. A summary of the most importantmethylKit features is shown in a flow chart. It depicts the main features of methylKitand the sequential relationship between them. The functions that could be used for thosefeatures are also printed in the boxes.

Flexible data integration and regional analysis

High-throughput bisulfite sequencing experiments typically yield millions of reads with reducedcomplexity due to cytosine conversion, and there are several different aligners suited for mappingthese reads to the genome (see Frith et al. [19] and Krueger et al. [20] for a review and comparison between aligners). Since methylKit only requires amethylation score per base for all analyses, it is a modular package that can be applied independentof any aligner. Currently, there are two ways that information can be supplied tomethylKit:: 1) methylKit can read per base methylation scores from a text file(see Table 1 for an example of such a file); and, 2) methylKit canread SAM format [21] alignments files obtained from Bismark aligner [22]. If a SAM file is supplied, methylkit first processes the alignment file to get%methylation scores and then reads that information into memory.

Table 1 Sample text file that can be read by methylKit.

Most bisulfite experiments have a set of test and control samples or samples across multipleconditions, and methylKit can read and store (in memory) methylation data simultaneouslyfor N-experiments, limited only by memory of the node or computer. The default setting of theprocessing algorithm requires that there be least 10 reads covering a base and each of the basescovering the genomic base position have at least 20 PHRED quality score. Also, since DNA methylationcan occur in CpG, CHG and CHH contexts (H = A, T, or C) [3], users of methylKit have the option to provide methylation information for allthese contexts: CpG, CHG and CHH from SAM files.

Summarizing DNA methylation information over pre-defined regions or tiling windows

Although base-pair resolution DNA methylation information is obtained through most bisulfitesequencing experiments, it might be desirable to summarize methylation information over tilingwindows or over a set of predefined regions (promoters, CpG islands, introns, and so on). Forexample, Smith et al. [9] investigated methylation profiles with RRBS experiments on gametes and zygote andsummarized methylation information on 100bp tiles across the genome. Their analysis revealed aunique set of differentially methylated regions maintained in early embryo. Using tiling windows orpredefined regions, such as promoters or CpG islands, is desirable when there is not enoughcoverage, when bases in close proximity will have similar methylation profiles, or where methylationproperties of a region as a whole determines its function. In accordance with these potentialanalytic foci, methylKit provides functionality to do either analysis on tiling windowsacross the genome or predefined regions of the genome. After reading the base pair methylationinformation, users can summarize the methylation information on pre-defined regions they select oron tiling windows covering the genome (parameter for tiles are user provided). Then, subsequentanalyses, such as clustering or differential methylation analysis, can be carried out with the samefunctions that are used for base pair resolution analysis.

Example methylation data set: breast cancer cell lines

We demonstrated the capabilities of methylKit using an example data set from sevenbreast cancer cell lines from Sun et al. [23]. Four of the cell lines express estrogen receptor-alpha (MCF7, T47D, BT474, ZR75-1), andfrom here on are referred to as ER+. The other three cell lines (BT20, MDA-MB-231, MDA-MB-468) donot express estrogen receptor-alpha, and from here on are referred to as ER-. It has been previouslyshown that ER+ and ER- tumor samples have divergent gene expression profiles and that those profilesare associated with disease outcome [24, 25]. Methylation profiles of these cell lines were measured using reduced RRBS [10]. The R objects contained the methylation information for breast cancer cell lines andfunctions that produce plots and other results that are shown in the remainder of this manuscriptare in Additional file 4.

Whole methylome characterization: descriptive statistics, sample correlation and clustering

Descriptive statistics on DNA methylation profiles

Read coverage per base and % methylation per base are the basic information contained in themethylKit data structures. methylKit has functions for easy visualization of suchinformation (Figure 2a and 2b for % methylation and readcoverage distributions, respectively - for code see Additional file 4). Innormal cells, % methylation will have a bimodal distribution, which denotes that the majority ofbases have either high or low methylation. The read coverage distribution is also an importantmetric that will help reveal if experiments suffer from PCR duplication bias (clonal reads). If suchbias occurs, some reads will be asymmetrically amplified and this will impair accurate determinationof % methylation scores for those regions. If there is a high degree of PCR duplication bias, readcoverage distribution will have a secondary peak on the right side. To correct for this issue,methylKit has the option to filter bases with very high read coverage.

Figure 2
figure 2

Descriptive statistics per sample. (a) Histogram of %methylation per cytosine forER+ T47D sample. Most of the bases have either high or low methylation. (b) Histogram of readcoverage per cytosine for ER+ T47D sample. ER+, estrogen receptor-alpha expressing.

Measuring and visualizing similarity between samples

We have also included methods to assess sample similarity. Users can calculate pairwisecorrelation coefficients (Pearson, Kendall or Spearman) between the %methylation profiles across allsamples. However, to ensure comparable statistics, a new data structure is formed before thesecalculations, wherein only cytosines covered in all samples are stored. Subsequently, pairwisecorrelations are calculated, to produce a correlation matrix. This matrix allows the user to easilycompare correlation coefficients between pairs of samples and can also be used to performhierarchical clustering using 1- correlation distance. methylKit can also further visualizesimilarities between all pairs of samples by creating scatterplots of the %methylation scores(Figure 3). These functions are essential for detecting sample outliers or forfunctional clustering of samples based on their molecular signatures.

Figure 3
figure 3

Scatter plots for sample pairs. Scatter plots of %methylation values for each pair inseven breast cancer cell lines. Numbers on upper right corner denote pair-wise Pearson's correlationscores. The histograms on the diagonal are %methylation histograms similar to Figure 2a for eachsample.

Hierarchical clustering of samples

methylKit can also be used to cluster samples hierarchically in a variety of ways. Theuser can specify the distance metric between samples ('1 - correlation' 'Euclidean', 'maximum','manhattan', 'canberra', 'binary' or 'minkowski') as well as the agglomeration method to be used inthe hierarchical clustering algorithm (for example, 'Ward's method', or 'single/complete linkage',and so on). Results can either be returned as a dendrogram object or a plot. Dendrogram plots willbe color coded based on user defined groupings of samples. For example, we found that most ER+ andER- samples clustered together except MDMB231 (Figure 4a). Moreover, the usermay be interested in employing other more model-intensive clustering algorithms to their data. Userscan easily obtain the %methylation data from methylKit object and perform their ownanalysis with the multitude of R-packages already available for clustering. An example of such aprocedure (k-means clustering) is shown in Additional file 4.

Figure 4
figure 4

Sample clustering. (a) Hierarchical clustering of seven breast cancer methylationprofiles using 1-Pearson's correlation distance. (b) Principal Component Analysis (PCA) of sevenbreast cancer methylation profiles, plot shows principal component 1 and principal component 2 foreach sample. Samples closer to each other in principal component space are similar in theirmethylation profiles.

Principal component analysis of samples

methylKit can be used to perform Principal Component Analysis (PCA) on the samples'%-methylation profiles (see for example [26]). PCA can reduce the high dimensionality of a data set by transforming the large numberof regions to a few principal components. The principal components are ordered so that the first fewretain most of the variation present in the original data and are often used to emphasize groupingstructure in the data. For example, a plot of the first two or three principal components couldpotentially reveal a biologically meaningful clustering of the samples. Before the PCA is performed,a new data matrix is formed, containing the samples and only those cytosines that are covered in allsamples. After PCA, methylKit then returns to the user a 'prcomp' object, which can be usedto extract and plot the principal components. We found that in the breast cancer data set, PCAreveals a similar clustering to the hierarchical clustering where MDMB231 is an outlier.

Differential methylation calculation

Parallelized methods for detecting significant methylation changes

Differential methylation patterns have been previously described in malignancies [27–29] and can be used to differentiate cancer and normal cells [30]. In addition, normal human tissues harbor unique DNA methylation profiles [7]. Differential DNA methylation is usually calculated by comparing methylation levelsbetween multiple conditions, which can reveal important locations of divergent changes between atest and a control set. We have designed methylKit to implement two main methods fordetermining differential methylation across all regions: logistic regression and Fisher's exacttest. However, the data frames in methylKit can easily be used with other statistical testsand an example is shown in Additional file 4 (using a moderated t-test,although we maintain that most natural tests for this kind of data are Fisher's exact and logisticregression based tests). For our example data set we compared ER+ to ER- samples, with our 'controlgroup' being the ER- set.

Method #1: logistic regression

In logistic regression, information from each sample is specified (the number of methylated Csand number of unmethylated Cs at a given region), and a logistic regression test will be applied tocompare fraction of methylated Cs across the test and the control groups. More specifically, at agiven base/region we model the methylation proportion Pi, for sample i= 1,...,n (where nis the number of biological samples) through the logistic regression model:

log ( P i /(1 -  P i )) =  β 0 + β 1 * T i
(1)

where Ti denotes the treatment indicator for sample i, Ti = 1 if sample iis in the treatment group and Ti = 0 if sample i is in control group. The parameterβ0 denotes the log odds of the control group and β1 the logoddsratio between the treatment and control group. Therefore, independent tests for all thebases/regions of interest are against the null hypothesis H0: β1= 0. Ifthe null hypothesis is rejected it implies that the logodds (and hence the methylation proportions)are different between the treatment and the control group and the base/region would subsequently beclassified as a differentially methylated cytosine (DMC) or region (DMR). However, if the nullhypothesis is not rejected it implies no statistically significant difference in methylation betweenthe two groups. One important consideration in logistic regression is the sample size and in manybiological experiments the number of biological samples in each group can be quite small. However,it is important to keep in mind that the relevant sample sizes in logistic regression are not merelythe number of biological samples but rather the total read coverages summed over all samples in eachgroup separately. For our example dataset, we used bases with at least 10 reads coverage for eachbiological sample and we advise (at least) the same for other users to improve power to detectDMCs/DMRs.

In addition, we have designed methylKit such that the logistic regression framework canbe generalized to handle more than two experimental groups or data types. In such a case, theinclusion of additional treatment indicators is analogous to multiple regression when there arecategorical variables with multiple groups. Additional covariates can be incorporated into model (1)by adding to the right side of the model:

α 1 * Covariat e 1,i + . . . + α K * Covariat e K,i

where Covariate1,i, ..., CovariateK,i denote K measured covariates(continuous or categorical) for sample i = 1,...,n and α1,..., αkdenote the corresponding parameters.

Method #2: Fisher's exact test

The Fisher's exact test compares the fraction of methylated Cs in test and control samples in theabsence of replicates. The main advantage of logistic regression over Fisher's exact test is that itallows for the inclusion of sample specific covariates (continuous or categorical) and the abilityto adjust for confounding variables. In practice, the number of samples per group will determinewhich of the two methods will be used (logistic regression or Fisher's exact test). If there aremultiple samples per group, methylKit will employ the logistic regression test. Otherwise,when there is one sample per group, Fisher's exact test will be used.

Following the differential methylation test and calculation of P-values, methylKitwill use the sliding linear model (SLIM) method to correct P-values to q-values [31], which corrects for the problem of multiple hypothesis testing [32, 33]. However, we also implemented the standard false discovery rate (FDR)-based method(Benjamini-Hochberg) as an option for P-value correction, which is faster but moreconservative. Finally, methylKit can use multi-threading so that differential methylationcalculations can be parallelized over multiple cores and be completed faster.

Extraction and visualization of differential methylation events

We have designed methylKit to allow a user to specify the parameters that define theDMCs/DMRs based on: q-value, %methylation difference, and type of differential methylation(hypo-/hyper-). By default, it will extract bases/regions with a q-value <0.01 and %methylationdifference >25%. These defaults can easily be changed when calling get.methylDiff()function. In addition, users can specify if they want hyper-methylated bases/regions(bases/regions with higher methylation compared to control samples) or hypo-methylated bases/regions(bases/regions with lower methylation compared to control samples). In the literature, hyper- orhypo-methylated DMCs/DMRs are usually defined relative to a control group. In our examples, and inmethylKit in general, a control group is defined when creating the objects through suppliedtreatment vector, and hyper-/hypomethylation definitions are based on that control group.

Furthermore, DMCs/DMRs can be visualized as horizontal barplots showing percentage of hyper- andhypo-methylated bases/regions out of covered cytosines over all chromosomes (Figure 5a). We observed higher levels of hypomethylation than hypermethylation in the breast cancercell lines, which indicates that ER+ cells have lower levels of methylation. Since another commonway to visualize differential methylation events is with a genome browser, methylKit canoutput bedgraph tracks (Figure 5b) for use with the UCSC Genome Browser orIntegrated Genome Viewer.

Figure 5
figure 5

Visualizing differential methylation events. (a) Horizontal bar plots show thenumber of hyper- and hypomethylation events per chromosome, as a percent of the sites with theminimum coverage and differential. By default this is a 25% change in methylation and all sampleswith 10X coverage. (b) Example of bedgraph file uploaded to UCSC browser. The bedraph file isfor differentially methylated CpGs with at least a 25% difference and q-value <0.01. Hyper- andhypo-methylated bases are color coded. The bar heights correspond to % methylation differencebetween ER+ and ER- sets. ER+, estrogen receptor-alpha expressing; ER-, estrogen receptor-alphanon-expressing. UCSC, University of California Santa Cruz.

Annotating differential methylation events

Annotation with gene models and CpG islands

To discern the biological impact of differential methylation events, each event must be put intoits genomic context for subsequent analysis. Indeed, Hansen et al. [34] showed that most variable regions in terms of methylation in the human genome are CpGisland shores, rather than CpG islands themselves. Thus, it is interesting to know the location ofdifferential methylation events with regard to CpG islands, their shores, and also the proximity tothe nearest transcription start site (TSS) and gene components. Accordingly, methylKit canannotate differential methylation events with regard to the nearest TSS (Figure 6a) and it also can annotate regions based on their overlap with CpG islands/shores andregions within genes (Figures 6b and 6c are output frommethylKit).

Figure 6
figure 6

Annotation of differentially methylated CpGs. (a) Distance to TSS fordifferentially methylated CpGs are plotted from ER+ versus ER- analysis. (b) Pie chartshowing percentages of differentially methylated CpGs on promoters, exons, introns and intergenicregions. (c) Pie chart showing percentages of differentially methylated CpGs on CpG islands,CpG island shores (defined as 2kb flanks of CpG islands) and other regions outside of shores and CpGislands. (d) Pie chart showing percentages of differentially methylated CpGs on enhancers andother regions. ER+, estrogen receptor-alpha expressing; ER-, estrogen receptor-alpha non-expressing,TSS, transcription start site.

Annotation with custom regions

As with most genome-wide assays, the regions of interest for DNA methylation analysis may bequite numerous. For example, several reports show that Alu elements are aberrantly methylated incancers [35, 36] and enhancers are also differentially methylated [37, 38]. Since users may need to focus on specific genomic regions and require customizedannotation for capturing differential DNA methylation events, methylKit can annotatedifferential methylation events using user-supplied regions. As an example, we identifieddifferentially methylated bases of ER+ and ER- cells that overlap with ENCODE enhancer regions [39], and we found a large proportion of differentially methylated CpGs overlapping with theenhancer marks, and then plotted them with methylKit (Figure 6d).

Analyzing 5-hydroxymethylcytosine data with methylKit

5-Hydroxymethylcytosine is a base modification associated with pluropotency, hematopoiesis andcertain brain tissues (reviewed in [40]). It is possible to measure base-pair resolution 5hmC levels using variations oftraditional bisulfite sequencing. Recently, Yu et al. [41] and Booth et al. [15] published similar methods for detecting 5hmC levels in base-pair resolution. Both methodsrequire measuring 5hmC and 5mC levels simultaneously and use 5hmC levels as a substrate to deducereal 5mC levels, since traditional bisulfite sequencing cannot distinguish between the two [42]. However, both the 5hmC and 5mC data generated by these protocols are bisulfitesequencing based, and the alignments and text files of 5hmC levels can be used directly inmethylKit. Furthermore, methylKit has an adjust.methylC() function toadjust 5mC levels based on 5hmC levels as described in Booth et al. [15].

Customizing analysis with convenience functions

methylKit is dependent on Bioconductor [43] packages such as GenomicRanges and its objects are coercible toGenomicRanges objects and regular R data structures such as data frames via providedconvenience functions. That means users can integrate methylKit objects to otherBioconductor and R packages and customize the analysis according to their needs or extend theanalysis further by using other packages available in R.

Conclusions

Methods for detecting methylation across the genome are widely used in research laboratories, andthey are also a substantial component of the National Institutes of Health's (NIH's) EpiGenomeroadmap and upcoming projects such as BLUEPRINT [44]. Thus, tools and techniques that enable researchers to process and utilize genome-widemethylation data in an easy and fast manner will be of critical utility.

Here, we show a large set of tools and cross-sample analysis algorithms built intomethylKit, our open-source, multi-threaded R package that can be used for any base-leveldataset of DNA methylation or base modifications, including 5hmC. We demonstrate its utility withbreast cancer RRBS samples, provide test data sets, and also provide extensive documentation withthe release.

Abbreviations

5hmC:

5-hydroxymethylcytosine

5mC:

5-methylcytosine

bp:

base pair

BS-seq:

:bisulfitesequencing

DMC:

differentially methylated cytosine

DMR:

differentially methylated region

ER:

estrogen receptor alpha

FDR:

false discovery rate

PCA:

principal component analysis

PCR:

polymerase chain reaction

RRBS:

reduced representation bisulfite sequencing

SLIM:

sliding linearmodel

TSS:

transcription start site.

References

  1. Deaton AM, Bird A: CpG islands and the regulation of transcription. Genes Dev. 2011, 25: 1010-2210. 10.1101/gad.2037511.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  2. Suzuki MM, Bird A: DNA methylation landscapes: provocative insights from epigenomics. Nat Rev Genet. 2008, 9: 465-476.

    Article  PubMed  CAS  Google Scholar 

  3. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery JR, Lee L, Ye Z, Ngo Q-M, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar AH, Thomson JA, Ren B, Ecker JR: Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009, 462: 315-322. 10.1038/nature08514.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  4. Bird AP, Wolffe AP: Methylation-induced repression--belts, braces, and chromatin. Cell. 1999, 99: 451-454. 10.1016/S0092-8674(00)81532-9.

    Article  PubMed  CAS  Google Scholar 

  5. Hendrich B, Bird A: Identification and characterization of a family of mammalian methyl-CpG binding proteins. Mol Cell Biol. 1998, 18: 6538-6547.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  6. Figueroa ME, Abdel-Wahab O, Lu C, Ward PS, Patel J, Shih A, Li Y, Bhagwat N, Vasanthakumar A, Fernandez HF, Tallman MS, Sun Z, Wolniak K, Peeters JK, Liu W, Choe SE, Fantin VR, Paietta E, Löwenberg B, Licht JD, Godley LA, Delwel R, Valk PJM, Thompson CB, Levine RL, Melnick A: Leukemic IDH1 and IDH2 mutations result in a hypermethylation phenotype, disrupt TET2 function,and impair hematopoietic differentiation. Cancer Cell. 2010, 18: 553-567. 10.1016/j.ccr.2010.11.015.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  7. Fernandez AF, Assenov Y, Martin-Subero JI, Balint B, Siebert R, Taniguchi H, Yamamoto H, Hidalgo M, Tan A-C, Galm O, Ferrer I, Sanchez-Cespedes M, Villanueva A, Carmona J, Sanchez-Mut JV, Berdasco M, Moreno V, Capella G, Monk D, Ballestar E, Ropero S, Martinez R, Sanchez-Carbayo M, Prosper F, Agirre X, Fraga MF, Graña O, Perez-Jurado L, Mora J, Puig S, et al: A DNA methylation fingerprint of 1628 human samples. Genome Res. 2012, 22: 407-419. 10.1101/gr.119867.110.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  8. Li E, Beard C: Role for DNA methylation in genomic imprinting. Nature. 1993, 366: 362-365. 10.1038/366362a0.

    Article  PubMed  CAS  Google Scholar 

  9. Smith ZD, Chan MM, Mikkelsen TS, Gu H, Gnirke A, Regev A, Meissner A: A unique regulatory phase of DNA methylation in the early mammalian embryo. Nature. 2012, 484: 339-344. 10.1038/nature10960.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  10. Meissner A, Mikkelsen TS, Gu H, Wernig M, Hanna J, Sivachenko A, Zhang X, Bernstein BE, Nusbaum C, Jaffe DB, Gnirke A, Jaenisch R, Lander ES: Genome-scale DNA methylation maps of pluripotent and differentiated cells. Nature. 2008, 454: 766-770.

    PubMed  CAS  PubMed Central  Google Scholar 

  11. Akalin A, Garrett-Bakelman FE, Kormaksson M, Busuttil J, Zhang L, Khrebtukova I, Milne TA, Huang Y, Biswas D, Hess JL, Allis D, Roeder RG, Valk PJM, Lo B, Paietta E, Tallman MS, Schroth GP, Mason CE, Melnick A, Figueroa ME: Base-pair resolution DNA methylation sequencing reveals profoundly divergent epigeneticlandscapes in acute myeloid leukemia. PLoS Genet. 2012, 8: e1002781-10.1371/journal.pgen.1002781.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  12. Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE: Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008, 452: 215-219. 10.1038/nature06745.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  13. Lister R, O'Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, Ecker JR: Highly integrated single-base resolution maps of the epigenome in Arabidopsis. Cell. 2008, 133: 523-536. 10.1016/j.cell.2008.03.029.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  14. Ball MP, Li JB, Gao Y, Lee J-H, LeProust EM, Park I-H, Xie B, Daley GQ, Church GM: Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells. Nat Biotechnol. 2009, 27: 361-368. 10.1038/nbt.1533.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  15. Booth MJ, Branco MR, Ficz G, Oxley D, Krueger F, Reik W, Balasubramanian S: Quantitative sequencing of 5-methylcytosine and 5-hydroxymethylcytosine at single-baseresolution. Science. 2012, 336: 934-937. 10.1126/science.1220671.

    Article  PubMed  CAS  Google Scholar 

  16. methylKit. [http://code.google.com/p/methylkit]

  17. Flusberg BA, Webster DR, Lee JH, Travers KJ, Olivares EC, Clark TA, Korlach J, Turner SW: Direct detection of DNA methylation during single-molecule, real-time sequencing. Nat Methods. 2010, 7: 461-465. 10.1038/nmeth.1459.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  18. Cherf GM, Lieberman KR, Rashid H, Lam CE, Karplus K, Akeson M: Automated forward and reverse ratcheting of DNA in a nanopore at 5-Ã… precision. Nat Biotechnol. 2012, 30: 344-348. 10.1038/nbt.2147.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  19. Frith MC, Mori R, Asai K: A mostly traditional approach improves alignment of bisulfite-converted DNA. Nucleic Acids Res. 2012, 40: e100-10.1093/nar/gks275.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  20. Krueger F, Kreck B, Franke A, Andrews SR: DNA methylome analysis using short bisulfite sequencing data. Nat Methods. 2012, 9: 145-151. 10.1038/nmeth.1828.

    Article  PubMed  CAS  Google Scholar 

  21. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009, 25: 2078-2079. 10.1093/bioinformatics/btp352.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Krueger F, Andrews SR: Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011, 27: 1571-1572. 10.1093/bioinformatics/btr167.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  23. Sun Z, Asmann YW, Kalari KR, Bot B, Eckel-Passow JE, Baker TR, Carr JM, Khrebtukova I, Luo S, Zhang L, Schroth GP, Perez EA, Thompson EA: Integrated analysis of gene expression, CpG island methylation, and gene copy number in breastcancer cells by deep sequencing. PloS One. 2011, 6: e17490-10.1371/journal.pone.0017490.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  24. van 't Veer LJ, Dai H, van de Vijver MJ, He YD, Hart AAM, Mao M, Peterse HL, van der Kooy K, Marton MJ, Witteveen AT, Schreiber GJ, Kerkhoven RM, Roberts C, Linsley PS, Bernards R, Friend SH: Gene expression profiling predicts clinical outcome of breast cancer. Nature. 2002, 415: 530-536. 10.1038/415530a.

    Article  PubMed  Google Scholar 

  25. Sotiriou C, Neo S-Y, McShane LM, Korn EL, Long PM, Jazaeri A, Martiat P, Fox SB, Harris AL, Liu ET: Breast cancer classification and prognosis based on gene expression profiles from apopulation-based study. Proc Natl Acad Sci USA. 2003, 100: 10393-10398. 10.1073/pnas.1732912100.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  26. Joliffe I: Principal Component Analysis. 2002, New York, USA, Springer, 2

    Google Scholar 

  27. Esteller M, Corn PG, Baylin SB, Herman JG: A gene hypermethylation profile of human cancer. Cancer Res. 2001, 61: 3225-3229.

    PubMed  CAS  Google Scholar 

  28. Baylin SB, Herman JG: DNA hypermethylation in tumorigenesis: epigenetics joins genetics. Trends Genet. 2000, 16: 168-174. 10.1016/S0168-9525(99)01971-X.

    Article  PubMed  CAS  Google Scholar 

  29. Costello JF, Frühwald MC, Smiraglia DJ, Rush LJ, Robertson GP, Gao X, Wright FA, Feramisco JD, Peltomäki P, Lang JC, Schuller DE, Yu L, Bloomfield CD, Caligiuri MA, Yates A, Nishikawa R, Su Huang H, Petrelli NJ, Zhang X, O'Dorisio MS, Held WA, Cavenee WK, Plass C: Aberrant CpG-island methylation has non-random and tumour-type-specific patterns. Nat Genet. 2000, 24: 132-138. 10.1038/72785.

    Article  PubMed  CAS  Google Scholar 

  30. Doi A, Park I-H, Wen B, Murakami P, Aryee MJ, Irizarry R, Herb B, Ladd-Acosta C, Rho J, Loewer S, Miller J, Schlaeger T, Daley GQ, Feinberg AP: Differential methylation of tissue- and cancer-specific CpG island shores distinguishes humaninduced pluripotent stem cells, embryonic stem cells and fibroblasts. Nat Genet. 2009, 41: 1350-1353. 10.1038/ng.471.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  31. Wang H-Q, Tuominen LK, Tsai C-J: SLIM: a sliding linear model for estimating the proportion of true null hypotheses in datasetswith dependence structures. Bioinformatics. 2011, 27: 225-231. 10.1093/bioinformatics/btq650.

    Article  PubMed  Google Scholar 

  32. Storey J: A direct approach to false discovery rates. J R Stat Soc Series B Stat Methodol. 2002, 64: 479-498. 10.1111/1467-9868.00346.

    Article  Google Scholar 

  33. Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003, 100: 9440-9445. 10.1073/pnas.1530509100.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  34. Hansen KD, Timp W, Bravo HC, Sabunciyan S, Langmead B, McDonald OG, Wen B, Wu H, Liu Y, Diep D, Briem E, Zhang K, Irizarry R a, Feinberg AP: Increased methylation variation in epigenetic domains across cancer types. Nat Genet. 2011, 43: 768-775. 10.1038/ng.865.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  35. Ehrlich M: DNA hypomethylation in cancer cells. Epigenomics. 2009, 1: 239-259. 10.2217/epi.09.33.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  36. Rodriguez J, Vives L, Jordà M, Morales C, Muñoz M, Vendrell E, Peinado MA: Genome-wide tracking of unmethylated DNA Alu repeats in normal and cancer cells. Nucleic Acids Res. 2008, 36: 770-784.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  37. Stadler MB, Murr R, Burger L, Ivanek R, Lienert F, Schöler A, Wirbelauer C, Oakeley EJ, Gaidatzis D, Tiwari VK, Schübeler D: DNA-binding factors shape the mouse methylome at distal regulatory regions. Nature. 2011, 480: 490-495.

    PubMed  CAS  Google Scholar 

  38. Wiench M, John S, Baek S, Johnson TA, Sung M-H, Escobar T, Simmons CA, Pearce KH, Biddie SC, Sabo PJ, Thurman RE, Stamatoyannopoulos JA, Hager GL: DNA methylation status predicts cell type-specific enhancer activity. EMBO J. 2011, 30: 3028-3039. 10.1038/emboj.2011.210.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  39. Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, Zhang X, Wang L, Issner R, Coyne M, Ku M, Durham T, Kellis M, Bernstein BE: Mapping and analysis of chromatin state dynamics in nine human cell types. Nature. 2011, 473: 43-49. 10.1038/nature09906.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  40. Branco MR, Ficz G, Reik W: Uncovering the role of 5-hydroxymethylcytosine in the epigenome. Nat Rev Genet. 2011, 13: 7-13.

    Article  PubMed  Google Scholar 

  41. Yu M, Hon GC, Szulwach KE, Song C-X, Zhang L, Kim A, Li X, Dai Q, Shen Y, Park B, Min J-H, Jin P, Ren B, He C: Base-resolution analysis of 5-hydroxymethylcytosine in the mammalian genome. Cell. 2012, 149: 1368-1380. 10.1016/j.cell.2012.04.027.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  42. Huang Y, Pastor WA, Shen Y, Tahiliani M, Liu DR, Rao A: The behaviour of 5-hydroxymethylcytosine in bisulfite sequencing. PloS One. 2010, 5: e8888-10.1371/journal.pone.0008888.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Leisch F, Li C, Maechler M, Rossini AJ, Sawitzki G, Smith C, Smyth G, Tierney L, Yang JYH, Zhang J: Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004, 5: r80-10.1186/gb-2004-5-10-r80.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Adams D, Altucci L, Antonarakis SE, Ballesteros J, Beck S, Bird A, Bock C, Boehm B, Campo E, Caricasole A, Dahl F, Dermitzakis ET, Enver T, Esteller M, Estivill X, Ferguson-Smith A, Fitzgibbon J, Flicek P, Giehl C, Graf T, Grosveld F, Guigo R, Gut I, Helin K, Jarvius J, Küppers R, Lehrach H, Lengauer T, Lernmark Å, Leslie D, et al: BLUEPRINT to decode the epigenetic signature written in blood. Nat Biotechnol. 2012, 30: 224-226. 10.1038/nbt.2153.

    Article  PubMed  CAS  Google Scholar 

Download references

Acknowledgements

We wish to acknowledge the invaluable contribution of the WCMC Epigenomics Core Facility. MEF issupported by the Leukemia & Lymphoma Society Special Fellow Award and a Doris Duke ClinicalScientist Development Award. FGB is supported by a Sass Foundation Judah Folkman Fellowship. AM issupported by an LLS SCOR grant (7132-08) and a Burroughs Wellcome Clinical Translational ScientistAward. AM and CEM are supported by a Starr Cancer Consortium grant (I4-A442). CEM is supported bythe National Institutes of Health (I4-A411, I4-A442, and 1R01NS076465-01).

Author information

Authors and Affiliations

Authors

Corresponding authors

Correspondence to Altuna Akalin or Christopher E Mason.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

AA designed methylKit, developed the first codebase, and added most features. MK designed the logistic regression based statistical test for methylKit and worked onstatistical modeling and initial clustering features. SL wrote some of the features in methylKitand prepared plots for the manuscript. MEF, FGB and AM tested the code and provided initialdata for development of methylKit. CEM supervised the work, tested code, and coordinatedtest data for validation. All authors have read and approved the manuscript for publication.

Electronic supplementary material

13059_2012_2930_MOESM1_ESM.gz

Additional file 1: methylKit v0.5.3. This version of methylKit is included for archival purposes only. Pleasedownload the most recent version from [16]. (GZ 1 MB)

13059_2012_2930_MOESM2_ESM.pdf

Additional file 2: methylKit User Guide. A vignette file to accompany the methylKit software package; themost recent software and vignette can be downloaded at [16]. (PDF 677 KB)

13059_2012_2930_MOESM3_ESM.pdf

Additional file 3: methylKit documentation. Documentation for functions and classes in the methylKit softwarepackage; the most recent software and documentation can be downloaded at [16]. (PDF 131 KB)

13059_2012_2930_MOESM4_ESM.r

Additional file 4: R script for example analysis. The file contains R commands that are needed to do analysisand to produce graphs used in this manuscript. The file contains both the commands and detailedcomments on how those commands can be used. An up to date version of this script will beconsistently maintained at [16]. (R 5 KB)

Authors’ original submitted files for images

Rights and permissions

Reprints and permissions

About this article

Cite this article

Akalin, A., Kormaksson, M., Li, S. et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol 13, R87 (2012). https://doi.org/10.1186/gb-2012-13-10-r87

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/gb-2012-13-10-r87

Keywords