Skip to main content

METHimpute: imputation-guided construction of complete methylomes from WGBS data



Whole-genome bisulfite sequencing (WGBS) has become the standard method for interrogating plant methylomes at base resolution. However, deep WGBS measurements remain cost prohibitive for large, complex genomes and for population-level studies. As a result, most published plant methylomes are sequenced far below saturation, with a large proportion of cytosines having either missing data or insufficient coverage.


Here we present METHimpute, a Hidden Markov Model (HMM) based imputation algorithm for the analysis of WGBS data. Unlike existing methods, METHimpute enables the construction of complete methylomes by inferring the methylation status and level of all cytosines in the genome regardless of coverage. Application of METHimpute to maize, rice and Arabidopsis shows that the algorithm infers cytosine-resolution methylomes with high accuracy from data as low as 6X, compared to data with 60X, thus making it a cost-effective solution for large-scale studies.


METHimpute provides methylation status calls and levels for all cytosines in the genome regardless of coverage, thus yielding complete methylomes even with low-coverage WGBS datasets. The method has been extensively tested in plants, but should also be applicable to other species. An implementation is available on Bioconductor.


Cytosine methylation (5mC) is a widely conserved epigenetic mark [14] with important roles in the regulation of gene expression and the silencing of transposable elements (TEs) and repeats [5, 6]. Experimentally-induced changes in 5mC patterns have been shown to affect plant phenotypes [79], rates of meiotic recombination [1013], genome stability [1418] and alter plant-environment interactions [1922]. Similar to genetic mutations, changes in 5mC patterns can also occur spontaneously as a result of errors in DNA methylation maintenance [2326]. There is substantial evidence in plants that experimentally-induced as well as spontaneously occurring 5mC changes can be stably inherited across multiple generations, independently of genetic changes [27]. Cytosine methylation has therefore emerged as a potentially important factor in plant evolution [2830] and as a possible molecular target for the improvement of commercial crops [31, 32].

Plant methylomes are now routinely studied using whole-genome bisulfite sequencing (WGBS), a next generation sequencing (NGS) method that can interrogate the methylation status of individual cytosines at the genome-wide scale. The application of this technology has been instrumental in dissecting the molecular pathways that establish and maintain 5mC patterns in plant genomes. Unlike in animals, plants methylate cytosines in context CG, but also extensively in contexts CHG and CHH, where H= A, T, C [5]. Methylation at CG dinucleotides (mCG) is maintained by methyltransferase 1 (MET1), which is recruited to hemi-methylated CG sites in order to methylate the complimentary strand in a template-dependent manner during DNA replication [33]. By contrast, mCHG is maintained dynamically by the plant specific chromomethylase 3 (CMT3) [34], and requires continuous interactions with H3K9me2 (dimethylation of lysine 9 on histone 3) [35]. Asymmetrical methylation of CHH sites (mCHH) is established and maintained by another member of the CMT family, CMT2 [2, 36]. Similar to CMT3, CMT2 dynamically methylates CHH in H3K9me2-associated regions. In addition to these context-specific maintenance mechanisms, all three sequence contexts can also be methylated de novo via RNA-directed DNA methylation (RdDM) [5], which involves short-interfering 24 nucleotide small RNAs (siRNA) that guide the de novo methyltransferase domains rearranged methyltransferase 2 (DRM2) to homologous target sites throughout the genome [37, 38].

Although these methylation pathways appear to be broadly conserved across plant species, recent data indicates that there is extensive variation in 5mC patterns both between but also within species [3, 39]. Efforts to explore the origin of this variation and its implications for plant evolution, ecology and agriculture will require large inter- and intraspecific methylome datasets. Such datasets are currently emerging. To date, the methylomes of over 50 plant species have been analyzed using WGBS [3, 4], including representative species of major taxonomic groups such as angiosperms (flowering plants), gymnosperms, ferns, and non-vascular plants. In addition, the methylomes of over 1000 natural A. thaliana accessions are now available [40], as well as those of several experimentally derived populations [41]. However, deep inter- and intraspecific WGBS measurements remain cost-prohibitive, particularly for species with large genomes. Most published plant methylomes have therefore been sequenced far below saturation (i.e. large number of cytosines in the genome are not covered). Indeed, even simple genomes, like that of the model plant A. thaliana (Col-0 accession), are typically only sequenced to about 10-30X. At this depth, about 5-10% of cytosines have missing data (i.e. zero read coverage) and about 15-20% have nearly uninformative read coverage (<3 reads), and this problem is exacerbated in more complex genomes, like those of rice and maize (see Fig. 1).

Fig. 1
figure 1

Coverage distributions. a-c Percentage of cytosines with X coverage (strand-specific). d-f Percentage of cytosines with missing data (red) and “uninformative” coverage (green), defined as less than three reads

Low to moderate sequencing depths in individual samples have cumulative consequences for analyzing population-level data. For instance, in the recently released 1000 A. thaliana methylome data [40] (measured at 5X coverage per strand on average), 92% of cytosines have missing data in at least one sample when 100 accessions are compared (Additional file 1: Figure S1). These incomplete measurements will reduce statistical power in genome-wide methylation QTL (meQTL) mapping studies, in estimates of epimutation rates, or in ecological studies that aim to correlate site-specific methylation levels with environmental/climatic variables. Moreover, incomplete measurements also complicate and potentially bias methylome scans for signatures of epigenetic selection using methylation site frequency spectrum (mSFS) analytic approaches [28]. One way to circumvent the missing data problem is to calculate methylation levels over larger regions, ranging anywhere from several hundred to several thousand basepairs and to use these methylation levels for downstream population-level analyses. In the above-mentioned A. thaliana population data, only 36% of 100bp regions in the genome are missing in at least one sample of the 100 accessions, compared with 92% of individual cytosines, and this percentage further decreases with larger region sizes. However, while region-based methylation levels are useful measures for descriptive and correlative analyses, these measures obscure detailed insights into the base-resolution methylation status calls, and thus arguably undermine the key advantages of WGBS over other lower resolution technologies such as MeDIP-seq. Base-resolution status calls are needed to be able to apply existing population (epi)genetic models to population methylome data, and to be able to test explicit hypotheses about the evolutionary forces that shape methylome variation patterns within and among species [28].

In order to maximize the information contained in WGBS data and to facilitate cost-effective sequencing decisions for future studies, we developed METHimpute, a Hidden Markov Model (HMM) based imputation algorithm for the construction of base-resolution methylomes from WGBS data. The unique feature of this algorithm is its ability to impute the methylation status and level of cytosines with missing or uninformative coverage, thus yielding complete methylomes even with low-coverage WGBS datasets. Indeed, using published WGBS data from Arabidopsis thaliana (thale cress), Oryza sativa (rice) and Zea mays (maize), we demonstrate that METHimpute accurately constructs base-resolution methylomes from data with an average coverage as low as 6X, suggesting that typical sequencing costs could be cut without a significant loss of information.


Conceptual overview

WGBS is an NGS-based method in which DNA is treated with sodium bisulfite prior to sequencing in order to convert unmethylated cytosines into uracils and finally into thymines during PCR amplification. Hence, a cytosine in a bisulfite treated read that maps to a cytosine in the reference genome provides evidence for methylation, while a thymine that maps to a cytosine does not. Many specialized short read mapping programs make use of this information and output so-called methylation levels [4244]; that is, the proportion of aligned reads that support that a cytosine is methylated out of all the reads covering the site. Methylation levels are inherently noisy due to inefficiencies in the sodium bisulfite conversion step. Moreover, tissue heterogeneity and the highly dynamic maintenance methylation at CHH and CHG, which requires feedback loops with histone modifications and small RNAs [5, 6], lead to intermediate methylation levels which are very susceptible to experimental variation. Finally, in WGBS data a large proportion of cytosines are often either not covered by any sequencing read or are covered only by a few number of reads (Fig. 1), meaning that methylation levels at these positions cannot be determined.

To overcome these limitations we developed METHimpute, a Hidden Markov Model (HMM) for the construction of base-resolution methylomes from WGBS data. METHimpute takes methylated and unmethylated read counts at every cytosine as input, and outputs discrete methylation status calls (unmethylated or methylated), together with recalibrated methylation levels between 0 and 1 for every cytosine in the genome, regardless of coverage (Fig. 2).

Fig. 2
figure 2

Conceptual overview of METHimpute. a Cytosines on the sequenced genome are assumed to be either unmethylated or methylated. b Bisulphite-sequencing and alignment yields methylation levels for each cytosine, i.e. the number of reads showing methylation divided by the total number of reads. c Emission densities for each state are obtained with a binomial test with state-specific parameters. Note that "imputed" cytosines, i.e. cytosines without any reads, are treated identically as all other cytosines. However, since the emission densities for all states are 1 for imputed cytosines, the methylation status call is purely driven by the neighborhood of cytosines. d Model fitting yields posterior probabilities for methylation status calls. e Inferred methylation status calls and methylation levels

The METHimpute algorithm fits a two-state HMM to the observed methylation counts. The two hidden states correspond to the unmethylated (U) and methylated (M) components, with component-specific binomial emission densities. The estimates of the binomial parameters (pU and pM) and the HMM transition matrix (i.e. the collection of probabilities to transition from one hidden state to another) are estimated freely during model training for different sequence contexts, thus requiring no empirical knowledge of the conversion rate. In the present analysis we have used contexts CG, CCG, CWG, CAA, CTA and CCA |CHY (where H={A,C,T},W={A,T} and Y={C,T}), following evidence of their different methylation characteristics [45]. If necessary the model could be extended to account for different emission and transition parameters for every context and annotation category (genes, TEs, CpG density, etc.), instead of context alone, although this would only be possible for well-annotated genomes.

Based on the model fits, the probability that a given cytosine belongs to one of the hidden states is given by the posterior probabilities γU and γM (Fig. 2d and “Methods” section). A cytosine’s maximum posterior probability represents its most likely methylation status (Fig. 2d, e), and the magnitude of this probability can be used as a measure of confidence in the underlying status call. In addition to methylation status calls, METHimpute outputs recalibrated methylation levels per cytosine, calculated as m=pU·γU+pM·γM (Fig. 2e). A key feature of METHimpute is its ability to infer the methylation level and status for cytosines with missing data (i.e. zero read coverage) or for those with poor read coverage (i.e. less than 3 reads). It achieves this inference iteratively during HMM training by borrowing information from neighboring sites. The algorithm therefore outputs complete, base-resolution methylomes, that can otherwise only be obtained through very high-depth sequencing experiments.

Imputation-guided construction of complete Arabidopsis, rice and maize methylomes

To demonstrate the performance of METHimpute we analyzed representative WGBS datasets from A. thaliana (Col-0) (replicate 1: 8.6X; rep.2: 15.7X coverage per cytosine per strand) [41], O. sativa (japonica nipponbare) (rep.1: 7.4X, rep.2.: 6.9X, rep.3: 4.6X) [46], and Z. mays (B73) (rep.1: 1.6X, rep.2: 3.3X, rep.3: 2.4X) [47]. These three species cover a wide spectrum of plant genomes in terms of length and complexity: the A. thaliana, O. sativa and Z. mays genomes are 120 Mb, 374 Mb and 2.1 Gb in length, respectively, and have an estimated repeat content of 10, 28-35 and 85% [4851]. The mapping statistics for each dataset are detailed in Additional file 2: Table S1. Alignment and pre-processing of the data was carried out using a single pipeline as described in the “Methods” section. Runtimes and memory requirements for METHimpute are listed in Additional file 2: Table S4.

Despite average coverage being relatively high, a substantial proportion of cytosines had either missing data or low coverage. For instance, in the A. thaliana (rep.1: 8.6X), O. sativa (rep.3: 4.6X) and Z. mays (rep.2: 3.3X) datasets, about 9% (3.71M), 24% (39.54M) and 26% (36.77M) of all cytosines had missing data (i.e. zero read coverage) and 24% (10.27M), 49% (79.38M) and 60% (85.5M) were nearly uninformative (here defined as coverage <3 reads) (Fig. 1d-f and Additional file 1: Figure S2 for the other replicates). Interestingly, the genome-wide proportions of missing or uninformative sites were highly context dependent, being highest for CCA |CHY, probably as a result of less unique short read alignments in this context as it is more abundant in repetitive regions of the genome (Additional file 1: Figures S3 and S4).

We applied METHimpute to the above-described datasets and evaluated the quality of the resulting methylation calls. For A. thaliana, O. sativa and Z. mays, the algorithm imputed the methylation status of all 3.71M, 39.54M and 36.77M missing data cytosines, respectively, and inferred the methylation status of all 10.27M, 79.38M and 85.5M uninformative cytosines.

Inferred methylation calls capture known biology

To evaluate the quality of the inferred methylation status calls and levels we examined the per-cytosine posterior probability of being either unmethylated (U) or methylated (M). This probability represents a measure of statistical confidence in the underlying methylation call, with a value of 1 being the most confident. We found that the distribution of maximum posterior probability values for imputed cytosines shows a clear peak around 1 and a tail of lower confidence values (Fig. 3 and Additional file 1: Figure S5 for the other replicates), suggesting that the algorithm produces high-confidence methylation calls for a large proportion of cytosines with missing data. Indeed, 58% (1.50M), 54% (3.96M) and 83% (6.43M) of imputed cytosines in A. thaliana, O. sativa and Z. mays were called with high confidence (defined as posterior probability ≥0.9), and these numbers increased to 91% (4.16M), 90% (6.64M) and 93% (9.56M) for cytosines covered by only one or two reads.

Fig. 3
figure 3

Maximum posterior distributions for imputed cytosines (coverage = 0), uninformative cytosines (coverage = 1 or 2) and informative cytosines (coverage ≥3). For Arabidopsis (a), rice (b) and maize (c), for each context. The figure shows the distributions of the maximum posterior probabilities with density on the y-axis and the maximum posterior probability on x-axis. The maximum posterior probability, i.e. the confidence in the methylation status calls, is generally lower for sites with less coverage

To assess whether the inferred methylation levels are consistent with known biology, we constructed meta-methylation profiles for annotated repeats and genes using cytosines separated in three different categories: informative (coverage ≥3), uninformative (coverage = 1 or 2) and imputed cytosines (coverage =0). Regardless of coverage category, METHimpute confirms that A. thaliana TE sequences are heavily methylated in all sequence contexts, with a marked decrease in methylation levels at their 5’ and 3’ ends (Fig. 4b and Additional file 1: Figure S6b for the other replicate). The CCA |CHY context shows the lowest methylation levels and CG shows the highest, consistent with Gouil and Baulcombe 2016 [45], and the ordering is conserved for imputed and uninformative cytosines. Similar profiles were detected for repeat elements in O. sativa and Z. mays, with high CG, CCG and CWG methylation, and very low levels of CAA, CTA, and particularly CCA |CHY methylation, consistent with known results (Fig. 4d, f and Additional file 1: Figure S6 for the other replicates) [52].

Fig. 4
figure 4

Enrichment profiles for genes (left panels) and transposable elements or repeats (right panels) for Arabidopsis (a, b), rice (c, d) and maize (e, f), for each context. Sub-panels show the enrichment profiles for imputed (coverage = 0), uninformative (coverage = 1 or 2) and informative cytosines (coverage ≥3). See the “Methods” section for definition of the recalibrated methylation level

In line with numerous methylome studies in Arabidopsis (e.g. [45, 53, 54]), METHimpute finds that A. thaliana genes are intermediately methylated in CG context, and essentially unmethylated at all other contexts. (Fig. 4a and Additional file 1: Figure S6a for the other replicate). Genic meta-methylation profiles for O. sativa and Z. mays were generally similar to those of A. thaliana (Fig. 4c, e and Additional file 1: Figure S6 for the other replicates), with the exception that both crop species are known to also methylate genic CHG context, probably owing to the fact that genes in these complex genomes often overlap or contain heavily methylated TE or repeat copies.

Taken together the above analyses illustrate two points: first, METHimpute infers annotation-specific methylation profiles that are consistent with published reports; and second, the methylation profiles inferred from imputed or uninformative cytosines recapitulate the patterns seen for highly-informative cytosines, indicating that - regardless of coverage - the inferred methylation calls are robust and biologically meaningful.

Saturation analysis for the performance assessment of imputed methylomes

METHimpute achieves high quality imputations by leveraging information from neighboring cytosines via the estimated distance-dependent transition probabilities (see “Methods” section). Therefore, confidence in the imputed calls is higher for cytosines that are closer to informative sites (Additional file 1: Figure S7). This spatial dependency remains high over distances of 10-40 bp and then decays to background levels. We reasoned that our imputation method may therefore be relatively robust even in shallow WGBS experiments.

To test this directly, we implemented a saturation analysis similar to Libertini et al. 2016 [55], where we compared high-coverage datasets with low-coverage subsets of these datasets. Bam files with mapped reads for the Arabidopsis, rice and maize replicates were merged to obtain samples with 23.2X, 18.6X and 7.2X coverage per cytosine per strand, respectively (Additional file 2: Table S1). These merged files were downsampled to generate a series of reduced datasets, ranging from 90 to 10% of the original data (Additional file 2: Table S3). Upon downsampling, the proportion of cytosines with zero read coverage increased from 5% (23.2X) to 31% (13.47M, 2.6X) in A. thaliana, from 11% (18.6X) to 40% (65.41M, 1.8X) in O. sativa and from 14% (7.2X) to 37% (52.07M, 2.2X) in the Z. mays data (Fig. 5d-f). We ran METHimpute on each reduced dataset and calculated the F1-score in the status calls relative to those obtained with the full data. The F1-score is defined as the harmonic mean of precision and recall, and the status calls of the full dataset were assumed as ground truth.

Fig. 5
figure 5

Saturation analysis. a-c F1-score for METHimpute and the binomial test, compared to the full sample, respectively. The F1-score is the harmonic mean of precision and recall. d-f Proportion of imputed cytosines. g-i Proportion of the genome in each state. The x-axis shows the average strand-specific coverage per cytosine

Our analysis shows that performance remains remarkably high despite drastic decreases in sequencing depth (Fig. 5a-c, Additional file 1: Figure S8 with precision and recall, Additional file 1: Figure S9 F1-score per context). With data as low as 5X coverage per cytosine (strand-specific), the F1-score was as high as 95% (U: 95%, M: 74%) in Arabidopsis, 97% (U: 97%, M: 88%) in rice and 99% (U: 99% M: 98%) in maize. In general, annotations with a large percentage of missing cytosines in the high coverage datasets were less accurately called upon downsampling (Additional file 1: Figure S4). These include in particular transposable elements and repeats. The exception to this trend were 5’ UTRs, which in all three species showed a large percentage of cytosines with missing data but a low amount of miscalled sites upon downsampling.

To put the above accuracy analysis into perspective, we compared the HMM-based imputation method with a much simpler method based on the commonly used binomial test: Methylation states for informative cytosines (>=3 reads) were called with a binomial test, and methylation states for missing and uninformative cytosines (<3 reads) were imputed by assigning the majority methylation state of covered cytosines of the same context in the 200bp neighborhood of the missing or uninformative cytosine. Cytosines without any informative neighbors in a 200bp neighborhood were not imputed and treated as “undefined”, and therefore counted as false negatives in the downsampled data if the full dataset was informative in these positions. We find that the accuracy obtained with this approach is less robust to average sequencing depth. With only 5X data, the F1-score drops down to 93% (U: 93%, M: 74%) in Arabidopsis, 94% (U: 93%, M: 84%) in rice and 95% (U: 93%, M: 91%) in maize (Fig. 5a-c).

Finally, we also considered the fidelity of the recalibrated methylation levels upon downsampling. Recalibrated methylation levels can be interpreted as the probability of observing a methylated read at a given position, and they are highly correlated with original methylation levels: For Arabidopsis, rice and maize, the correlation (linear fit) was 0.91, 0.94 and 0.93, respectively (p-value ≤2e−16). To assess their fidelity upon downsampling, we calculated the correlation between recalibrated methylation levels per cytosine and per 100 bp window to the full coverage dataset, and compared that to the results obtained from the original methylation level (Additional file 1: Figure S10). Per-cytosine recalibrated methylation levels show slightly higher correlations than original methylation levels, and with 10% of the original data the correlations for Arabidopsis, rice and maize are 0.89, 0.90 and 0.93, respectively. Window-based recalibrated methylation levels showed the same correlation performance as the original ones, with remarkably high correlations even when only 10% of the original data was retained (0.95, 0.95, 0.83 for Arabidopsis, rice, maize). These results suggest that recalibrated methylation levels can be used for downstream methylation analysis, since they are correlated to original methylation levels and are robust upon downsampling, while providing base-resolution information even at low sequencing depth.

Overall, both for status calls and for recalibrated methylation levels, METHimpute produces robust results even at very low sequencing depth, suggesting that the algorithm offers a cost-effective solution for methylome studies of large genomes and for population-level studies involving a large number of samples.

Re-calibrated estimates of genome-wide and context-specific methylation levels

Plant species differ greatly in their genome-wide methylation levels (GMLs, i.e. the proportion of cytosines that are methylated) [3, 4]. In a recent survey of about 30 angiosperms, GMLs were found to be as low as 5% in Theobroma cacao to as high as 43% in Beta vulgaris, with a mean of about 16% [3, 39]. Much of this diversity appears to be the result of differences in genome size and repeat content, as well as differences in the efficiency of DNA methylation maintenance pathways [28]. Precise estimates of GMLs are important for studying the evolutionary forces that shape plant methylomes over short and long time-scales, and for understanding genome-epigenome co-evolution. However, obtaining GML estimates based on WGBS data is not trivial, as they are highly dependent on the method used for methylation status calling and on the depth of the sequencing experiment. In A. thaliana, for instance, reported GML estimates vary widely between studies. This dependency is even larger when considering context-specific GMLs (i.e. the proportion of methylated cytosines in a given context; CG-GMLs, CHG-GMLs, CHH-GMLs), with CHH-GMLs being by far the most variable between studies, with reported values ranging from as low as 1.51% [1] to as high as 3.91% [3].

In order to bypass many of the statistical issues in calling methylation states, especially in shallow WGBS data, recent studies have proposed so-called weighted genome-wide methylation levels (wGMLs) as a proxy for GMLs. A wGML is a non-statistical measure which is obtained by counting the number of methylated reads over the total number of reads at the genome-wide scale. Figure 5g-i shows clearly that wGMLs are robust upon down-sampling in any sequence context in the A. thaliana, O. sativa and Z. mays data, thus justifying its use. By contrast, GMLs calculated from base-resolution binomial status calls (i.e. #mC/all Cs) are highly unstable, particularly in non-CG contexts and when sequencing depth is low (Fig. 5g-i).

In order to assess whether the re-calibrated methylation levels provided by METHimpute can also be used to obtain robust estimates of GMLs, we calculated wGMLs by summing the per-cytosine re-calibrated methylation level genome-wide, weighted by coverage. Using this measure we find that METHimpute-derived wGMLs perform nearly identical to naive wGMLs, both in terms of robustness and magnitude (Fig. 5g-i, Additional file 1: Figure S11 with replicates). This demonstrates that METHimpute recalibrated levels are consistent with original methylation levels and known biology not only at the individual cytosine level, but also aggregated over 100bp windows and genome-wide, with the added advantage that they are available for all positions in the genome.

METHimpute facilitates insights into bisulfite conversion rates

One source of measurement noise in WGBS data is the bisulfite conversion procedure prior to sequencing. Bisulfite treatment of DNA is typically performed long enough so that all unmethylated cytosines are converted to uracils. The conversion success (or rate) is typically high. Most studies report conversion rates of about 0.99, implying that only about 1% of all unmethylated cytosines failed to convert. Knowledge of this rate is important not only to verify that bisulfite reaction was efficient but also to be able to separate biological signal from noise in downstream analyses of the data. Empirical estimates of the conversion rate are often obtained by including unmethylated chloroplast and virus genomes as controls in the WGBS workflow, and counting the number of non-converted cytosines from the mapped reads.

A helpful byproduct of the METHimpute fitting procedure is that the conversion rate can be directly estimated from the sequenced material without requiring auxiliary information from chloroplast or virus genomes. METHimpute achieves this in the HMM framework by estimating the probability, pU, of finding a methylated read given that the underlying cytosine is unmethylated (see “Methods” section), which can be used to derive the conversion rate. To obtain these rates we focus on estimates of pU in context CG to exclude potential biases arising from the “fuzzy” maintenance of methylation at CHG and CHH sites. For A. thaliana and Z. mays our estimated conversion rates were 0.989 and 0.961, respectively, which is remarkably close to chloroplast-based estimates of 0.993 and 0.970.

Although bisulfite conversion kits and protocols have been optimized to achieve the highest conversion rate possible the specificity of the reaction is not perfect. A well-known trade-off is that some methylated cytosines can be accidentally converted to uracils, and are later falsely detected as unmethylated. Some controls (commercial or artificially methylated DNA fragments) are available to estimate this inappropriate conversion rate, but, to our knowledge, they are not systematically used in WGBS experiments. Some studies using such controls have shown that the inappropriate conversion rate (% of methylated cytosines converted to uracils) ranges from 0.09 to 6.1% depending on the kit and protocol used [5658].

METHimpute approximates this value by estimating the parameter pM for the M component (see Methods), which can be used to calculate the probability of finding an unmethylated cytosine given that the underlying cytosine is truly methylated. Again, focusing on CG sites, we estimate the methylated cytosines conversion rate at 6.3, 11.5 and 16% in O. sativa, Z. mays and A. thaliana, respectively. Although these estimates are close to the empirical rates reported in the literature, they are biased upward most likely owing to the fact that the parameter pM is partly confounded with methylation variation arising from cellular heterogenity in the sampled tissues. We therefore suspect that our estimates become more accurate in situations where tissue heterogeneity is minimized.

Nonetheless, the ability of METHimpute to provide an accurate estimate of the conversion rate for unmethylated cytosines and an upper-bound estimate for methylated cytosines could be utilized to calibrate WGBS experiments in the laboratory when no controls are available.


A key advantage of WGBS over alternative measurement technologies is its ability to provide base-resolution measurements from bulk and - more recently - also from single cell data. Since its first application in the model plant A. thaliana in 2008 [53, 54], WGBS has become an integral tool for studying the methylomes of increasingly large plant genomes and for surveying patterns of natural methylome variation within and among plant species. However, the relatively high costs associated with this technology pose limits on the sequencing depths that can be achieved within most experimental budgets. A typical solution is to sequence methylomes far below saturation, which results in substantial measurement noise and missing data at the level of individual cytosines.

Here we presented METHimpute, a Hidden Markov Model for the construction and imputation of complete methylomes from shallow or deep WGBS data. We show that our approach imputes high-confidence methylation calls for uncovered cytosines that are sufficiently close to informative cytosines (<40 bp). For cytosines in widely uncovered regions, our approach imputes low-confidence calls which might be filtered out in downstream analyses, since they do not contain more information than background frequencies of methylation. A threshold for filtering can be determined from the asymptotic behavior of the maximum posterior probability as in Additional file 1: Figure S7. In the presented analysis we have used six different sequence contexts, but our implementation is general enough to allow also other context specifications. For example, the model could be extended to account for different emission and transition parameters for every context and annotation category (genes, TEs, CpG density, etc.), instead of context alone, although this would only be possible for well-annotated genomes.

We recommend the use of METHimpute instead of the binomial test for the analysis of WGBS data whenever methylation status calls are required. Furthermore, METHimpute solves the problem of missing data in population epigenetic studies, which will facilitate the estimation of epigenetic mutation rates and methylation site frequency spectrum analyses.


Here we introduced METHimpute, an imputation-based HMM for the construction of complete methylomes from shallow or deep WGBS data. Our analyses show that the algorithm can impute the methylation status of cytosines with missing (i.e. zero read coverage) or uninformative coverage (i.e. coverage of less than 3 reads), as well as their recalibrated methylation levels. We demonstrate that these imputations are not only statistically robust, but also biologically meaningful. Our estimates suggest that routine use of this algorithm could reduce sequencing costs of typically sized methylome experiments without a substantial loss of biological information. The method works with small, streamlined genomes like that of Arabidopsis but also with large, repeat-rich genomes like those of most commercial crops, thus making it a flexible software tool for the analysis of DNA methylomes of a wide spectrum of species. METHimpute is implemented as an R-package and seamlessly integrates with the extensive bioinformatic tool sets available through Bioconductor. The algorithm has been extensively tested in plants, but it should also be applicable in non-plant species.


Hidden Markov Model for methylation calling


We define an N=2 state Hidden Markov Model (HMM), where the states i represent unmethylated (U) and methylated (M) cytosines. The emission densities for each state are binomial distributions, which can be interpreted as a binomial test on the number of methylated counts m over total counts r. The probability parameter pi of the binomial test can be interpreted as the probability of finding m methylated counts out of r total counts, given the state i. Note that in this definition 1−pU is the conversion rate, i.e. the probability of a read showing non-methylation when the cytosine is indeed non-methylated. Cytosines are not equally spaced in the genome, and we therefore chose a distance dependent transition matrix A, where the distance dependent change in transition probabilities is modeled by an exponential function. Furthermore, to account for different sequence contexts, we implemented context-specificity for both the binomial test and the transition probabilities.

Mathematical description

The probability P of observing methylated mt and total rt read count at a particular cytosine t in context ct can be written as

$$ P_{t}\left(m_{t}, r_{t}, \mathbf{p}_{c_{t}}\right) = \sum_{i \in \{U,M\}} \gamma_{it} B_{{ic}_{t}}\left(m_{t}, r_{t}, p_{{ic}_{t}}\right), $$

where γi are the posteriors (mixing weights) and Bi are binomial distributions with context-specific parameter pic. The binomial distribution is defined as

$$ B(m, r, p) = {{r}\choose{m}} p^{m} (1-p)^{r-m}. $$

All probability parameters of the binomial tests (i.e. the probabilities of a success) are estimated freely during model training (next section). For C=6 contexts and N=2 states, N·C=12 independent parameters pic need to be fitted.

The distance dependent transition probabilities from cytosine t in state i to cytosine t+1 in state j, separated by distance dt,t+1 and in transition context ct,t+1, can be described as

$$ {\begin{aligned} A_{ij,c_{t,t+1}}\left(A^{o}_{ij,c_{t,t+1}}, d_{t,t+1}, D_{c_{t,t+1}}, N\right) &= A^{o}_{ij,c_{t,t+1}} \mathrm{e}^{-d_{t,t+1}/D_{c_{t,t+1}}} \\ & \quad + \frac{1}{N} \left(1-\mathrm{e}^{-d_{t,t+1}/D_{c_{t,t+1}}}\right). \end{aligned}} $$

Here, \(A^{o}_{ij,c_{t,t+1}}\) are the transition probabilities without distance dependency (or for adjacent cytosines with dt,t+1=0). \(D_{c_{t,t+1}}\) is a constant that reflects how fast neighboring cytosines lose correlation. The distance dependency is constructed in such a way that all transitions \(A_{ij,c_{t,t+1}}\) are equally likely for an infinite distance dt,t+1=. For C=6 contexts the model has C·C=36 transition contexts and thus 36 different transition matrices with dimensions N×N.

The constants Dc are determined by a non-linear least-squares (nls) fit to the correlation decay between cytosines in transition context ct,t+1 (see Additional file 1: Figure S12 for all used transition contexts). The formula for the fit is \(y_{c}(d) = a0 * \mathrm {e}^{-d/D_{c}}\), where yc is the correlation between neighboring cytosines at distance d in transition context c. The parameters a0 and Dc are fitted by the nls-fit.

The correlation is calculated between adjacent cytosines, with no other cytosines in between. This reflects the definition of the transition probabilities in the Hidden Markov Model, where transitions are defined from one cytosine to the next in the sequence.

Model fitting

Model parameters are fitted with the Baum-Welch algorithm [59]. The distance-dependent transition probabilities require modified updating formulas compared to a standard Baum-Welch algorithm without distance dependency. The derivation of the modified updating formulas is detailed below, and uses notation introduced in [60].

The conditional expectation Q that needs to be maximized can be written as

$$ {\begin{aligned} Q &= \sum_{i}^{N} \gamma_{i,t=0} ~ log\left(\pi_{i}\right) + \sum_{i,j,t}^{N,N,T-1} \xi_{ijt} ~ log\left(A_{ij,c_{t,t+1}}\right)\\ & \quad + \sum_{i,t}^{N,T} \gamma_{it} ~ log\left(f_{i}\right). \end{aligned}} $$

The updated transition probabilities \(A^{\prime {o}}_{ijc}\) can be obtained by solving \(\frac {\partial Q}{\partial A^{o}_{ijc}} = 0\) using the method of Lagrange multipliers to deal with the constraint \(\sum ^{N}_{j} A^{o}_{ijc} = 1\).

$$ {\begin{aligned} A^{\prime{o}}_{ijc} &= \left(\sum_{t}^{T-1} \delta_{c,c_{t,t+1}} ~ \xi_{ijt} ~ \frac{A^{o}_{ijc}}{A_{ij,c_{t,t+1}}} \frac{\partial A_{ij,c_{t,t+1}}}{\partial A^{o}_{ij,c_{t,t+1}}} \right) \\ & \quad \left/ \left(\sum_{t,j}^{T-1,N} \delta_{c,c_{t,t+1}} ~ \xi_{ijt} ~ \frac{A^{o}_{ijc}}{A_{ij,c_{t,t+1}}} \frac{\partial A_{ij,c_{t,t+1}}}{\partial A^{o}_{ij,c_{t,t+1}}} \right)\right.. \end{aligned}} $$

Here, \(\delta _{c,c_{t,t+1}}\) is the Kronecker delta function, which ensures that only terms in the correct transition context c are included into the sum.

Similarly, the updated parameters for the binomial test can be obtained by solving \(\frac {\partial Q}{\partial p_{ic}} = 0\). For independent binomial tests, this yields

$$ p^{\prime}_{ic} = \left(\sum^{T}_{t} \delta_{c,c_{t}} ~ \gamma_{it} ~ m_{t} \right) \left/ \left(\sum^{T}_{t} \delta_{c,c_{t}} ~ \gamma_{it} ~ r_{t} \right)\right.{.} $$

The methylation status it is determined by maximizing over the posterior probabilities it=argmaxi(γit).

Finally, we can use the posterior probabilities γU|M,t and estimated parameters pic to define a recalibrated methylation level \(m^{\prime }_{t}\) that is defined on every cytosine t in the genome and can serve as input for other applications:

$$ m^{\prime}_{t} = p_{U,c_{t}} \cdot \gamma_{U,t} + p_{M,c_{t}} \cdot \gamma_{M,t} $$

Plants DNA methylation data

We used published data (fastq files containing bisulfite sequencing reads) from three model plant species: Arabidopsis thaliana, rice (Oryza sativa Japonica cv. Nipponbare) and maize (Zea mays B73). We used three replicates for rice and maize, and two replicates for Arabidopsis. Each sample was mapped to the latest available version of the reference genome for this species. Details and references on these datasets, reference genomes and annotations files, as well as additional alignment metrics can be accessed in Additional file 2: Table S2.

Mapping of bisulphite sequenced (BS-seq) reads and construction of DNA methylomes

Read sequences (Additional file 2: Table S2) were quality trimmed and adapter sequences were removed with Cutadapt (version 1.9; python version 2.7.9; [61]). Trimming was performed on both ends using the single-end mode and the quality threshold was set to a phred score of 20 (q=20). We applied the default error rate of 10% for the removal of the adapter sequences. Afterwards, we discarded reads shorter than 40 base pairs. Reads were subsequently mapped to an indexed genome. The maximum allowed proportion of mismatches was set to 0.05 (m=0.05, 5 mismatches per 100 bp) and the maximum insert size was set to 1000 bp (X=1000). BS-Seeker2 (v2.0.10; [44]) using Bowtie2 (version 2.2.2; [62]) was chosen for the alignment of the reads. Samtools (version 1.3.1; using htslib 1.2.1; [63]) was used to remove duplicates (samtools rmdup -s) and to sort bam files (samtools sort). Methylomes were subsequently constructed through the module from BS-Seeker2 (v2.0.10; [44]). CGmap files containing methylome information were used as an input for METHimpute.



5-methyl cytosine


Cytosine followed by a guanine


Genome wide methylation level


Hidden Markov Model


Methylated state of the HMM


Methylated CG


Methylation quantitative trait locus


Methylation site frequency spectrum


Next generation sequencing


Transposable element


Unmethylated state of the HMM


Whole-genome bisulfite sequencing


Weighted genome wide methylation level


  1. Feng S, Cokus SJ, Zhang X, Chen P-Y, Bostick M, Goll MG, Hetzel J, Jain J, Strauss SH, Halpern ME, Ukomadu C, Sadler KC, Pradhan S, Pellegrini M, Jacobsen SE. Conservation and divergence of methylation patterning in plants and animals. Proc Natl Acad Sci. 2010; 107(19):8689–94.

    Article  PubMed  Google Scholar 

  2. Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-Wide Evolutionary Analysis of Eukaryotic DNA Methylation. Science. 2010; 328(5980):916–9.

    Article  PubMed  CAS  Google Scholar 

  3. Niederhuth CE, Bewick AJ, Ji L, Alabady M, Kim KD, Page JT, Li Q, Rohr NA, Rambani A, Burke JM, Udall JA, Egesi C, Schmutz J, Grimwood J, Jackson SA, Springer NM, Schmitz RJ. Widespread natural variation of DNA methylation within angiosperms. Genome Biol. 2016;17(194).

  4. Takuno S, Ran J-H, Gaut BS. Evolutionary patterns of genic DNA methylation vary across land plants. Nat Plants. 2016; 2(January):15222.

    Article  PubMed  CAS  Google Scholar 

  5. Law JA, Jacobsen SE. Establising, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010; 11(3):204–20.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. Matzke Ma, Kanno T, Matzke AJM. RNA-Directed DNA Methylation: The Evolution of a Complex Epigenetic Pathway in Flowering Plants. Annu Rev Plant Biol. 2014; (December 2014):1–25.

    Google Scholar 

  7. Cortijo S, Wardenaar R, Colome-Tatche M, Gilly A, Etcheverry M, Labadie K, Caillieux E, Hospital F, Aury J-M, Wincker P, Roudier F, Jansen RC, Colot V, Johannes F. Mapping the Epigenetic Basis of Complex Traits. Science. 2014; 343(6175):1145–8.

    Article  PubMed  CAS  Google Scholar 

  8. Johannes F, Porcher E, Teixeira FK, Saliba-Colombani V, Simon M, Agier N, Bulski A, Albuisson J, Heredia F, Audigier P, Bouchez D, Dillmann C, Guerche P, Hospital F, Colot V. Assessing the impact of transgenerational epigenetic variation on complex traits. PLoS Genet. 2009;5(6).

  9. Reinders J, Wulff BBH, Mirouze M, Marí-Ordóñez A, Dapp M, Rozhon W, Bucher E, Theiler G, Paszkowski J. Compromised stability of DNA methylation and transposon immobilization in mosaic Arabidopsis epigenomes. Genes and Development. 2009; 23(8):939–50.

    Article  PubMed  CAS  Google Scholar 

  10. Mirouze M, Lieberman-Lazarovich M, Aversano R, Bucher E, Nicolet J, Reinders J, Paszkowski J. Proc Natl Acad Sci USA. 2012; 109(15):5880–5.

  11. Yelina NE, Lambing C, Hardcastle TJ, Zhao X, Santos B, Henderson IR. DNA methylation epigenetically silences crossover hot spots and controls chromosomal domains of meiotic recombination in Arabidopsis. Genes Dev. 2015; 29(20):2183–202.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Colome-Tatche M, Cortijo S, Wardenaar R, Morgado L, Lahouze B, Sarazin A, Etcheverry M, Martin A, Feng S, Duvernois-Berthet E, Labadie K, Wincker P, Jacobsen SE, Jansen RC, Colot V, Johannes F. Features of the Arabidopsis recombination landscape resulting from the combined loss of sequence variation and DNA methylation. Proc Natl Acad Sci. 2012; 109(40):16240–5.

    Article  PubMed  Google Scholar 

  13. Melamed-Bessudo C, Levy aa. PNAS Plus: Deficiency in DNA methylation increases meiotic crossover rates in euchromatic but not in heterochromatic regions in Arabidopsis. Proc Natl Acad Sci. 2012; 109(16):981–8.

    Article  Google Scholar 

  14. Tsukahara S, Kobayashi A, Kawabe A, Mathieu O, Miura A, Kakutani T. Bursts of retrotransposition reproduced in Arabidopsis. Nature. 2009; 461(7262):423–6.

    Article  PubMed  CAS  Google Scholar 

  15. Mirouze M, Reinders J, Bucher E, Nishimura T, Schneeberger K, Ossowski S, Cao J, Weigel D, Paszkowski J, Mathieu O. Selective epigenetic control of retrotransposition in Arabidopsis. Nature. 2009; 461(September):1–5.

    Google Scholar 

  16. Miura A, Yonebayashi S, Watanabe K, Toyama T, Shimada H, Kakutani T. Mobilization of transposons by a mutation abolishing full DNA methylation in Arabidopsis. Nature. 2001; 411(May):212–4.

    Article  PubMed  CAS  Google Scholar 

  17. Singer T, Yordan C, Martienssen RA. Robertson’s Mutator transposons in A. thaliana are regulated by the chromatin-remodeling gene Decrease in DNA Methylation (DDM1). Genes Dev. 2001; 15(5):591–602.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. Cheng C, Tarutani Y, Miyao A, Ito T, Yamazaki M, Sakai H, Fukai E, Hirochika H. Loss of function mutations in the rice chromomethylase OsCMT3a cause a burst of transposition. Plant J. 2015; 83(6):1069–81.

    Article  PubMed  CAS  Google Scholar 

  19. Secco D, Wang C, Shou H, Schultz MD, Chiarenza S, Nussaume L, Ecker JR, Whelan J, Lister R. Stress induced gene expression drives transient DNA methylation changes at adjacent repetitive elements. eLife. 2015; 4(July):09343.

    Google Scholar 

  20. Zhang X. Dynamic differential methylation facilitates pathogen stress response in Arabidopsis. Proc Natl Acad Sci. 2012; 109(32):12842–3.

    Article  PubMed  Google Scholar 

  21. Yu A, Lepère G, Jay F, Wang J, Bapaume L, Wang Y, Abraham A-L, Penterman J, Fischer RL, Voinnet O, Navarro L. Proc Natl Acad Sci USA. 2013; 110(6):2389–94.

  22. López Sánchez A, Stassen JHM, Furci L, Smith LM, Ton J. The role of DNA (de)methylation in immune responsiveness of Arabidopsis. Plant J. 2016; 88(3):361–74.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Schmitz RJ, Schultz MD, Lewsey MG, O’Malley RC, Urich MA, Libiger O, Schork NJ, Ecker JR. Transgenerational Epigenetic Instability Is a Source of Novel Methylation Variants. Science. 2011; 334(6054):369–73.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Becker C, Hagmann J, Müller J, Koenig D, Stegle O, Borgwardt K, Weigel D. Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. Nature. 2011; 480(7376):245–9.

    Article  PubMed  CAS  Google Scholar 

  25. Jiang C, Mithani A, Belfield EJ, Mott R, Hurst LD, Harberd NP. Environmentally responsive genome-wide accumulation of de novo Arabidopsis thaliana mutations and epimutations. Genome Res. 2014; 24(11):1821–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. van der Graaf A, Wardenaar R, Neumann DA, Taudt A, Shaw RG, Jansen RC, Schmitz RJ, Colomé-Tatché M, Johannes F. Rate, spectrum, and evolutionary dynamics of spontaneous epimutations. Proc Natl Acad Sci USA. 2015; 112(21):6676–81.

    Article  PubMed  CAS  Google Scholar 

  27. Quadrana L, Colot V. Plant Transgenerational Epigenetics. Annu Rev Genet. 2016; 50(1):467–91.

    Article  PubMed  CAS  Google Scholar 

  28. Vidalis A, živković D, Wardenaar R, Roquis D, Tellier A, Johannes F. Methylome evolution in plants. Genome Biol. 2016; 17(1):264.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. Diez CM, Roessler K, Gaut BS. Epigenetics and plant genome evolution. Curr Opin Plant Biol. 2014; 18(1):1–8.

    Article  PubMed  CAS  Google Scholar 

  30. Weigel D, Colot V. Epialleles in plant evolution. Genome Biol. 2012; 13(10):249.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Springer NM. Epigenetics and crop improvement. Trends Genet. 2013; 29(4):241–7.

    Article  PubMed  CAS  Google Scholar 

  32. Ji L, Neumann DA, Schmitz RJ. Crop Epigenomics: Identifying, Unlocking, and Harnessing Cryptic Variation in Crop Genomes. Mol Plant. 2015; 8(6):860–70.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Finnegan EJ, Peacock WJ, Dennis ES. Reduced DNA methylation in Arabidopsis thaliana results in abnormal plant development. Proc Natl Acad Sci USA. 1996; 93(16):8449–54.

    Article  PubMed  CAS  Google Scholar 

  34. Lindroth AM, Cao X, Jackson JP, Zilberman D, McCallum CM, Henikoff S, Jacobsen SE. Requirement of CHROMOMETHYLASE3 for Maintenance of CpXpG Methylation. Science. 2001; 292(5524):2077–80.

    Article  PubMed  CAS  Google Scholar 

  35. Du J, Zhong X, Bernatavichute Y, Stroud H, Feng S, Caro E, Vashisht A, Terragni J, Chin H, Tu A, Hetzel J, Wohlschlegel J, Pradhan S, Patel D, Jacobsen S. Dual Binding of Chromomethylase Domains to H3K9me2-Containing Nucleosomes Directs DNA Methylation in Plants. Cell. 2012; 151(1):167–80.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Stroud H, Do T, Du J, Zhong X, Feng S, Johnson L, Patel DJ, Jacobsen SE. Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis. Nat Struct Mol Biol. 2013; 21(1):64–72.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Cao X, Jacobsen SE. Locus-specific control of asymmetric and CpNpG methylation by the DRM and CMT3 methyltransferase genes. Proc Natl Acad Sci. 2002; 99(Supplement 4):16491–8.

    Article  PubMed  CAS  Google Scholar 

  38. Cao X, Jacobsen SE. Role of the arabidopsis DRM methyltransferases in de novo DNA methylation and gene silencing. Curr Biol. 2002; 12(13):1138–44.

    Article  PubMed  CAS  Google Scholar 

  39. Alonso C, PÃrez R, Bazaga P, Herrera CM. Global DNA cytosine methylation as an evolving trait: phylogenetic signal and correlated evolution with genome size in angiosperms. Front Genet. 2015; 6:4.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Kawakatsu T, Huang S-sC, Jupe F, Sasaki E, Schmitz RJ, Urich MA, Castanon R, Nery JR, Barragan C, He Y, Chen H, Dubin M, Lee CR, Wang C, Bemm F, Becker C, O’Neil R, O’Malley RC, Quarless DX, The 1001 Genomes Consortium, Weigel D, Nordborg M, Ecker JR. Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions. Cell. 2016; 166(2):492–506.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Stroud H, Greenberg MVC, Feng S, Bernatavichute YV, Jacobsen SE. Comprehensive Analysis of Silencing Mutants Reveals Complex Regulation of the Arabidopsis Methylome. Cell. 2013; 152(17):352–64.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. Krueger F, Andrews SR. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011; 27(11):1571–2.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, Mason CE. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012; 13(10):87.

    Article  Google Scholar 

  44. Guo W, Fiziev P, Yan W, Cokus S, Sun X, Zhang MQ, Chen P-Y, Pellegrini M, Cokus S, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild C, Pradhan S, Nelson S, Pellegrini M, Jacobsen S, Lister R, Pelizzola M, Dowen R, Hawkins R, Hon G, Tonti-Filippini J, Nery J, Lee L, Ye Z, Ngo Q-M, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar A, Thomson J, Ren B, Ecker J, Meissner A, Mikkelsen T, Gu H, Wernig M, Hanna J, Sivachenko A, Zhang X, Bernstein B, Nusbaum C, Jaffe D, Gnirke A, Jaenisch R, Lander E, Wang J, Xia Y, Li L, Gong D, Yao Y, Luo H, Lu H, Yi N, Wu H, Zhang X, Tao Q, Gao F, Chen P, Cokus S, Pellegrini M, Langmead B, Trapnell C, Pop M, Salzberg S, Krueger F, Andrews S, Harris E, Ponts N, Roch KL, Lonardi S, Pedersen B, Hsieh T-F, Ibarra C, Fischer R, Xi Y, Li W, Smith A, Chung W-Y, Hodges E, Kendall J, Hannon G, Hicks J, Xuan Z, Zhang M, Wu T, Nacu S, Xi Y, Bock C, Müller F, Sun D, Meissner A, Li W, Langmead B, Salzberg S, Li R, Li Y, Kristiansen K, Wang J, Smith A, Xuan Z, Zhang M, Giardine B, Riemer C, Hardison R, Burhans R, Elnitski L, Shah P, Zhang Y, Blankenberg D, Albert I, Taylor J, Miller W, Kent W, Nekrutenko A, Thorvaldsdóttir H, Robinson J, Mesirov J, Molaro A, Hodges E, Fang F, Song Q, McCombie W, Hannon G, Smith A. BS-Seeker2: a versatile aligning pipeline for bisulfite sequencing data. BMC Genomics. 2013; 14(1):774.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Gouil Q, Baulcombe DC. DNA Methylation Signatures of the Plant Chromomethyltransferases. PLoS Genet. 2016; 12(12):1006526.

    Article  CAS  Google Scholar 

  46. Stroud H, Ding B, Simon SA, Feng S, Bellizzi M, Pellegrini M, Wang GL, Meyers BC, Jacobsen SE. Plants regenerated from tissue culture contain stable epigenome changes in rice. eLife. 2013; 2013(2):1–14.

    Google Scholar 

  47. Regulski M, Lu Z, Kendall J, Donoghue MTA, Reinders J, Llaca V, Deschamps S, Smith A, Levy D, McCombie WR, Tingey S, Rafalski A, Hicks J, Ware D, Martienssen RA. The maize methylome influences mRNA splice sites and reveals widespread paramutation-like switches guided by small RNA. Genome Res. 2013; 23(10):1651–62.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. The Arabidopsis Genome Initiative. Analysis of the genome sequence of the flowering plant Arabidopsis thaliana. Nature. 2000; 408(6814):796–815.

    Article  Google Scholar 

  49. Sequencing Project IRG. The map-based sequence of the rice genome. Nature. 2005; 436(7052):793–800.

    Article  CAS  Google Scholar 

  50. Rice Annotation Project T, Itoh T, Tanaka T, Barrero RA, Yamasaki C, Fujii Y, Hilton PB, Antonio BA, Aono H, Apweiler R, Bruskiewich R, Bureau T, Burr F, Costa de Oliveira A, Fuks G, Habara T, Haberer G, Han B, Harada E, Hiraki AT, Hirochika H, Hoen D, Hokari H, Hosokawa S, Hsing YI, Ikawa H, Ikeo K, Imanishi T, Ito Y, Jaiswal P, Kanno M, Kawahara Y, Kawamura T, Kawashima H, Khurana JP, Kikuchi S, Komatsu S, Koyanagi KO, Kubooka H, Lieberherr D, Lin Y-C, Lonsdale D, Matsumoto T, Matsuya A, McCombie WR, Messing J, Miyao A, Mulder N, Nagamura Y, Nam J, Namiki N, Numa H, Nurimoto S, O’Donovan C, Ohyanagi H, Okido T, Oota S, Osato N, Palmer LE, Quetier F, Raghuvanshi S, Saichi N, Sakai H, Sakai Y, Sakata K, Sakurai T, Sato F, Sato Y, Schoof H, Seki M, Shibata M, Shimizu Y, Shinozaki K, Shinso Y, Singh NK, Smith-White B, Takeda J-i, Tanino M, Tatusova T, Thongjuea S, Todokoro F, Tsugane M, Tyagi AK, Vanavichit A, Wang A, Wing RA, Yamaguchi K, Yamamoto M, Yamamoto N, Yu Y, Zhang H, Zhao Q, Higo K, Burr B, Gojobori T, Sasaki T. Curated genome annotation of Oryza sativa ssp. japonica and comparative genome analysis with Arabidopsis thaliana. Genome Res. 2007; 17(2):175–83.

    Article  Google Scholar 

  51. Schnable PS, Ware D, Fulton RS, Stein JC, Wei F, Pasternak S, Liang C, Zhang J, Fulton L, Graves TA, Minx P, Reily AD, Courtney L, Kruchowski SS, Tomlinson C, Strong C, Delehaunty K, Fronick C, Courtney B, Rock SM, Belter E, Du F, Kim K, Abbott RM, Cotton M, Levy A, Marchetto P, Ochoa K, Jackson SM, Gillam B, Chen W, Yan L, Higginbotham J, Cardenas M, Waligorski J, Applebaum E, Phelps L, Falcone J, Kanchi K, Thane T, Scimone A, Thane N, Henke J, Wang T, Ruppert J, Shah N, Rotter K, Hodges J, Ingenthron E, Cordes M, Kohlberg S, Sgro J, Delgado B, Mead K, Chinwalla A, Leonard S, Crouse K, Collura K, Kudrna D, Currie J, He R, Angelova A, Rajasekar S, Mueller T, Lomeli R, Scara G, Ko A, Delaney K, Wissotski M, Lopez G, Campos D, Braidotti M, Ashley E, Golser W, Kim H, Lee S, Lin J, Dujmic Z, Kim W, Talag J, Zuccolo A, Fan C, Sebastian A, Kramer M, Spiegel L, Nascimento L, Zutavern T, Miller B, Ambroise C, Muller S, Spooner W, Narechania A, Ren L, Wei S, Kumari S, Faga B, Levy MJ, McMahan L, Van Buren P, Vaughn MW, Ying K, Yeh C-T, Emrich SJ, Jia Y, Kalyanaraman A, Hsia A-P, Barbazuk WB, Baucom RS, Brutnell TP, Carpita NC, Chaparro C, Chia J-M, Deragon J-M, Estill JC, Fu Y, Jeddeloh JA, Han Y, Lee H, Li P, Lisch DR, Liu S, Liu Z, Nagel DH, McCann MC, SanMiguel P, Myers AM, Nettleton D, Nguyen J, Penning BW, Ponnala L, Schneider KL, Schwartz DC, Sharma A, Soderlund C, Springer NM, Sun Q, Wang H, Waterman M, Westerman R, Wolfgruber TK, Yang L, Yu Y, Zhang L, Zhou S, Zhu Q, Bennetzen JL, Dawe RK, Jiang J, Jiang N, Presting GG, Wessler SR, Aluru S, Martienssen RA, Clifton SW, McCombie WR, Wing RA, Wilson RK. The B73 Maize Genome: Complexity, Diversity, and Dynamics. Science. 2009; 326(5956):1112–5.

    Article  PubMed  CAS  Google Scholar 

  52. West PT, Li Q, Ji L, Eichten SR, Song J, Vaughn MW, Schmitz RJ, Springer NM. Genomic distribution of H3K9me2 and DNA methylation in a maize genome. PLoS ONE. 2014; 9(8):1–10.

    Article  CAS  Google Scholar 

  53. Lister R, Malley RCO, Tonti-filippini J, Gregory BD, Berry CC, Millar a. H, Ecker JR. Highly Integrated Single-Base Resolution Maps of the Epigenome in Arabidopsis. Cell. 2008; 133(3):523–36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  54. 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(7184):215–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  55. Libertini E, Heath SC, Hamoudi RA, Gut M, Ziller MJ, Herrero J, Czyz A, Ruotti V, Stunnenberg HG, Frontini M, Ouwehand WH, Meissner A, Gut IG, Beck S. Saturation analysis for wholegenome bisulfite sequencing data. Nat Publ Group. 2016:11–13.

  56. Holmes EE, Jung M, Meller S, Leisse A, Sailer V, Zech J, Mengdehl M, Garbe L-A, Uhl B, Kristiansen G, Dietrich D. Performance Evaluation of Kits for Bisulfite-Conversion of DNA from Tissues, Cell Lines, FFPE Tissues, Aspirates, Lavages, Effusions, Plasma, Serum, and Urine. PLoS ONE. 2014; 9(4):93933.

    Article  CAS  Google Scholar 

  57. Leontiou CA, Hadjidaniel MD, Mina P, Antoniou P, Ioannides M, Patsalis PC. Bisulfite Conversion of DNA: Performance Comparison of Different Kits and Methylation Quantitation of Epigenetic Biomarkers that Have the Potential to Be Used in Non-Invasive Prenatal Testing. PLoS ONE. 2015; 10(8):0135058.

    Article  CAS  Google Scholar 

  58. Genereux DP, Johnson WC, Burden AF, Stoger R, Laird CD. Errors in the bisulfite conversion of DNA: modulating inappropriate- and failed-conversion frequencies. Nucleic Acids Res. 2008; 36(22):150.

    Article  CAS  Google Scholar 

  59. Baum LE, Petrie T, Soules G, Weiss N. A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains. Ann Math Stat. 1970; 41(1):164–71.

    Article  Google Scholar 

  60. Rabiner L. A tutorial on hidden Markov models and selected applications in speech recognition. Proc IEEE. 1989; 77(2):257–286.

    Article  Google Scholar 

  61. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011; 17(1):10–12.

    Article  Google Scholar 

  62. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012; 9(4):357–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. 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(16):2078–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references


We thank R. Schmitz and C. Niederhuth with their help in accessing the maize and rice annotation files and providing data, and N. Springer and R. Schmitz for their quick feedback on this project.


FJ and DR acknowledge support from the Technical University of Munich-Institute for Advanced Study funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement #291763. MCT acknowledges support from the Helmholtz Association’s Initiative and Networking Fund and from the University of Groningen (Rosalind Franklin Fellowship).

Availability of data and materials

METHimpute can be downloaded from The data sources used in this publication are listed in Additional file 2: Table S1.

Author information

Authors and Affiliations



AT, MC-T and FJ conceived this project; AT implemented the model with input from MC-T, FJ and DR; AT, DR, RW and AV analyzed the data; AT, MC-T and FJ wrote the paper. All authors have read and approved the manuscript.

Corresponding authors

Correspondence to Frank Johannes or Maria Colomé-Tatché.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1

Figures S1-S12. (PDF 1832 kb)

Additional file 2

Tables S1-S5. (XLSX 23 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Taudt, A., Roquis, D., Vidalis, A. et al. METHimpute: imputation-guided construction of complete methylomes from WGBS data. BMC Genomics 19, 444 (2018).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Methylation
  • Whole-genome bisulfite sequencing
  • Imputation
  • Hidden Markov Model