BACKGROUND
methylated dna immunoprecipitation sequencing, or medip-seq, is a common methodology used to study dna methylation profiles within plant and animal genomes  <cit> .
dna methylation is increasingly being recognised as playing a central role in the regulation of gene expression and chromatin structure.
cytosine methylation is a particular type of dna methylation that most often refers to the reversible addition of a methyl group to the carbon- <dig> position of the cytosine pyrimidine ring, resulting in the formation of 5-methylcytosine  modifications within double stranded dna  <cit> .
these epigenetic dna modifications can be further classified based on the nucleobase context in which 5mc is observed.
double stranded dna methylation primarily occurs if the methylated cytosine occurs within cytosine-guanosine dinucleotide sequences i.e.
cg context .
this is due to nucleobase symmetry between sense  and antisense  dna strands.
however, methylation in a cg context is not always double-stranded; it can also occur on just a single strand.
single stranded  dna methylation can additionally be detected if the 5mc is in a non-cg methylation context, i.e.
a ch dinucleotide context where h = a, t or c .
asymmetric dna methylation is a common phenomenon of plant genomes  <cit> .
for example it has been associated with methylation changes on retrotransposons during the development of pollen cells  <cit> .
however, for many years it was thought to be absent in mammals where much of the initial focus was on the detection of symmetric cg methylation.
nevertheless, the presence of non-cg methylation in mammalian embryonic stem and induced pluripotent stem cells  <cit> , as well as somatic tissues, e.g.
adult brain  is becoming more commonly observed.
although asymmetric dna methylation might not have the same function as symmetric double-stranded methylation , there is preliminary evidence implying a link between non-cg methylation and transcriptional regulation  <cit> .
in contrast to symmetric cg methylation, cytosine methylation in a non-cg context is prevalent within gene bodies, particularly exons, rather than surrounding the transcriptional start site  <cit> .
there is substantial evidence suggesting that non-cg methylation serves a distinct function depending on the cell type.
for instance, non-cg methylated regions are prevalent in undifferentiated cells such as embryonic stem cells  <cit> , as well as differentiated cells such as brain and skeletal muscle  <cit> .
additionally, non-cg methylation is generally associated with a positive correlation to gene expression: this modification is thought to act as a repressive mark in the adult mammalian brain and appears to be crucial for neural function  <cit> .
dna methylation studies of other animals  initially were hampered by assumptions translated from mammalian systems, such as the predomination of methylation occurring at cg dinucleotide sequences  <cit> .
recently, however, this view is being challenged with the honey bee, apis mellifera, becoming an increasingly important invertebrate model for studying how dna methylation  affects metazoan development, behaviour and gene splicing .
for example, cingolani et al.
<cit>  have described how biases in both software and experimental design towards studying cg methylation have led to dna methylation in other dinucleotide contexts being poorly identified and functionally underappreciated.
after correcting for such biases, they were surprised to discover 5-fold more chh methylation than cg methylation in a. mellifera.
they report that of all the cytosines present in cg, chg, and chh contexts in the genome,  <dig>  % were methylated; of these, only 21 % of cg cytosines were symmetrically methylated with 0 % of chh and  <dig>  % of chg cytosines symmetrically methylated.
in total, symmetric methylation occurred in only  <dig>  % of the methylated cytosine positions.
thus, asymmetric methylation predominates in this invertebrate species.
cingolani et al.
further suggest that non-cg modifications might have a significant influence on the regulation of alternative splicing, and this highlights the need to detect and quantify strand-specific methylation in other metazoan organisms  <cit> .
the current medip-seq library preparation procedure preserves strand identity, meaning that medip-seq reads contain both symmetric and asymmetric strand methylation signals that are detectable  <cit> .
however, this information is not exploited by any of the currently available software packages .
it therefore follows that the current downstream analysis methodology for medip-seq data is not optimal as asymmetric dna methylation is not detectable.
the ability to study strand-specific dna methylation signals would significantly increase the utility of medip-seq and lead to further insight into the function of asymmetric dna methylation.
we, therefore, have developed a new software  to decompose medip-seq derived dna methylation information into individual, strand-specific signals .
we have used dismiss to confidently assign strand-associated dna methylation signals in medip-seq datasets derived from a. mellifera.
we have also quantified how well dismiss predicts strand-associated dna methylation in a. mellifera when compared to data obtained from a technology that is able to resolve methylation signals at the nucleotide-level, namely bs-seq .
by demonstrating how dismiss can increase the resolution and sensitivity of medip-seq analysis of gdna , we present an enabling methodology that can be used to more accurately study dna methylomes.
implementation
strand conservation in medip-seq libraries
the retrieval of strand identity during medip-seq analyses of genomes is possible due to sequencing adapter directionality  <cit> .
the methodological step-by-step process of retrieving strand-specific information is illustrated in fig.
1a.
in part  ), genomic 5mc can be represented by three possible scenarios: firstly, the cytosine on both dna strands is methylated , secondly the 5mc is found only on the plus dna strand , and lastly when the 5mc occurs only on the minus dna strand .
during the first phase of library preparation , gdna is fragmented, end-repaired and the forked adapters are ligated directionally.
subsequently, the gdna is denatured  prior to the medip 5mc immuno-capture step , which leads to the selection of single stranded 5mc containing dna fragments with the 5′ adaptor shown as a bold solid black line – the presence of this adaptor is important for resolving the strand specific methylation signal.
following the 5mc immuno-capture step, library preparation continues by performing a pcr amplification step , resulting in double stranded dna fragments.
no matter which of the three scenarios the original 5mc is found in, the sequencing of first mate reads is initiated from the 5′ adaptor  – and due to the previous selection of asymmetrically methylated fragments  in step  this creates a disparity between the numbers of c nucleobases and g nucleobases in the first fastq file.
the strand origin of the original fragments may be recovered by aligning the sequenced reads to the genome sequence: reads containing 5mc from the minus strand will align in reverse-compliment mode.fig.
<dig> dismiss utilises strand biases generated during medip-seq library preparation to detect strand-specific dna methylation.
a outline of medip-seq library preparation.
genomic 5mc can be present in three possible scenarios:  with 5mc found on both dna strands,  with 5mc found only on the plus dna strand  and  with 5mc found only on the minus dna strand .
dna is fragmented, end-repaired and the forked adapters ligated directionally to the dna fragments – solid black  and dashed black ;  dna is denatured.
medip enrichment  selects single-stranded 5mc containing fragments; a sampling bias arises due to the selection of strands originating from dna fragments containing 5mc on only a single strand  – note the black adaptor is retained, identifying the strand origin of the selected fragment.
pcr library enrichment is performed .
black arrows show sequencing of first mate reads  occurring from the black adaptor.
dna fragments from all three scenarios are sequenced from both ends, with those originating from black and grey adaptors deposited in a single pair of first and second fastq files, respectively.
reads from the minus strand will tend to align to the genome sequence in reverse-compliment mode.
b c to g ratio bias detected in a medip-seq fastq file.
base count per cycle plot for the first fastq file.
the x-axis shows the base position in the read, while the y-axis shows the count of each nucleotide at that position.
the a to t ratios do not show a bias, as is expected for a random library not selecting for either of these nucleobases.
however, the c to g ratio demonstrates a clear bias due to the selection of 5mc containing fragments
the disparity in the g to c ratio of the medip reads can be readily observed using quality control software for high-throughput sequencing data .
for instance, software packages such as fastqc  are commonly used to assess a number of read metrics, including the ‘base frequency per cycle’.
in a typical random gdna library, the ratio between the counts of a and t nucleobases should be approximately equal, as should the g to c ratio.
in a medip-enriched gdna library, however, the plot will clearly show different counts for c and g, as shown in the ‘base frequency per cycle plot’ of fig.
1b.
this qc abnormality is more pronounced in cases where dna methylation is primarily found in an asymmetric context.
medip-seq strand separation
if the medip-seq library has been prepared following the methodological outline described in fig.
<dig> , then this information can be exploited by dismiss to assign 5mc strand specificity to each peak called by a peak caller .
thus, each dna methylation peak can be classified into one of the three classes: 1) minus stranded, 2) plus stranded or 3) peaks on both strands.
the methods required to perform this classification are implemented in the dismiss software and are described in the following paragraphs.
the count of medip reads mapping to a genomic location can be modelled as a poisson distribution  <cit> .
let y be the number of independent reads mapping to a genomic location that has previously been identified as a peak region by macs <dig>  the distribution of y is assumed to be poisson with the parameter lambda : y~poiλ
the reads mapped to the macs <dig> peak region are divided into two sub events representing the number of independent reads mapping to each of the two strands, yplus and yminus.
these sub events are also distributed according to a binomial distribution: yplus~poiλθandyminus~poiλ1−θ, where θ is a fraction given by yplus/y.
the goal of the analysis is to find out if there are integer values of λ that are plausible for both y <dig> and y <dig>  the likelihood function gives the probability of seeing a poisson distributed value y, given a rate λ:  <dig> lλ=pryλ
for all the plausible values of λ there is a unique value of λ that maximises the probability given by eq.
<dig>  called the maximum likelihood estimate of λ:  <dig> Λ=argmaxλlλ
in addition to finding the maximum likelihood estimate, we also want to know what other values of λ have a reasonably high likelihood in order to capture the randomness of the biological process.
a likelihood set is a set of λ values that explain the data in a statistically significant manner:  <dig> lsα≡λ:lλ/lΛ≥α
where α ∈  is the cut-off value and we frequently use α ≈  <dig>  for convenience and custom.
further details of the poisson distribution and the comparison of distributions using likelihood ratios are thoroughly explained elsewhere  <cit> .
for each peak region with y reads, the decision to assign a class is performed in the following manner.
given that yplus are the number of first mate reads aligned to the plus strand and yminus the number of first mate reads aligned to the minus strand, then a set of integer values from and including yplus to yminus that are plausible values of λ to explain the data are chosen.
using this set of λ values and eq.
<dig>  two likelihood sets are generated.
if the intersection of these two sets is not empty, then there are values of λ that are plausible for both yplus and yminus and the peak is assigned to class ‘peaks on both strands’.
otherwise the intersection of the sets is empty and the peak will be assigned as single-stranded: if yplus is greater than yminus then the peak is assigned to class ‘plus stranded’, else the peak is assigned to class ‘minus stranded’.
for example, if the number of reads mapping to the peak identified by macs <dig> is given by y =  <dig>  with the number of reads mapping to the two strands comprising this peak given by yplus
=  <dig> and yminus
=  <dig>  then the set of plausible values for λ  is given by { <dig>   <dig> …,  <dig>  60} and θ will be 60/100 =  <dig> .
using the poisson distributions for yplus and yminus and eq.
<dig>  the likelihood for each value of λ in this set is calculated and eq.
<dig> is used to select two values of λ that maximise the likelihood for each plus and minus strand.
furthermore, using eq.
<dig>  two likelihood sets  containing plausible values of λ  are generated: lsplus = { <dig>   <dig>  …,  <dig>  60} and lsminus = { <dig>   <dig>  …,  <dig>  55}.
the intersection of these two sets is not empty, hence there are common values of λ with a high likelihood for the observed number of reads yplus and yminus and the peak is assigned to both strands.
data sets used and analysis procedures
medip-seq data derived from a. mellifera  <cit>  consisted of adult worker nurse with sequence read archive  number srr <dig> , adult worker forager with sra number srr <dig>  and adult worker reverted-nurse with sra number srr <dig> .
the bs-seq data originated from a. mellifera worker larvae with sra number srx <dig>  <cit>  .
all data sets used in this study were retrieved using the sra software toolkit provided by the ncbi .fig.
<dig> dismiss predictions of the dna strand on which the 5mc occurs are not influenced by read depth.
a macs <dig> peaks.
the frequency of occurrence of peaks, with a certain depth of first mate reads, on the plus and minus strands.
the minimum read depth is  <dig> first mate reads.
b dismiss peaks.
the depth of first mate reads on the plus strand is plotted against the depth of first mate reads on the minus strand.
the predictions made by dismiss are indicated in grey scale and marked on the graph; they show well defined linear boundaries.
the minimum read depth is  <dig> first mate reads
fig.
<dig> dismiss interrogation of medip-seq datasets can accurately predict the dna strand where 5mc occurs.
a comparative analysis of 5mc detected by dismiss interrogation of medip-seq data with 5mc detected by bs-seq is illustrated  a despite originating from different developmental forms, 5mc detected by dismiss analysis of medip-seq datasets and bs-seq datasets are concordant.
each horizontal bar represents a dismiss quantified peak class  of medip-seq data.
from all of the dismiss predicted peak regions in a class, the percentage of 5mcs on the plus and minus strands was calculated from the bs-seq dataset.
each bar shows the percentage of 5mc identified by bs-seq on either plus strands  or minus strands .
the dismiss predicted ‘minus stranded’ and ‘plus stranded’ peaks predominantly contain bs-seq identified 5mcs on either minus or plus strands, respectively.
b loci-specific examples of 5mc detected by dismiss analysis of medip-seq datasets, with the corresponding bs-seq predictions, over three genomic regions.
three randomly selected regions from the a. mellifera genome , corresponding to dismiss predicted peak classes  are indicated.
the x-axis shows the nc_ <dig>  base pair coordinates of each selected region.
the first mate read coverage for each base using dismiss-assessed medip-seq alignment data is shown along the y-axis.
the y-axis is divided into two parts by a central line where the dark-grey coverage above the line represents reads predicted to map to the plus strand and the light-grey coverage shows reads predicted to map to the minus strand .
the black dots above the line indicate bs-seq identified 5mcs found on the plus strand, while those below the line show bs-seq identified 5mcs found on the minus strand
fig.
<dig> dna methylation is accurately predicted by dismiss over gene splice sites in a. mellifera.
the x-axis shows two genomic regions representing intron splice donor and acceptor sites.
the y-axis shows the percentage of single stranded medip-seq peaks , predicted by dismiss, on the sense  and antisense  sides of these regions.
dismiss detected more dna methylation on the antisense sides of the splice acceptor and donor sites, compared to the sense side, which is statistically significant  with p values less than  <dig>
fig.
<dig> dismiss detection of dna methylation in medip-seq datasets approaches the resolution obtained by bs-seq analysis.
a medip-seq dataset was analysed by macs <dig> and dismiss.
the proportion of dna methylation found in four different genomic regions  was deduced and compared to the proportion of dna methylation signal  detected in a bs-seq dataset
medip-seq and bs-seq data analysis
the medip-seq fastq files were aligned to the a. mellifera genome  <cit>  using bowtie <dig>  <cit> .
samtools  <cit>  was subsequently used to remove aligned reads with a mapping quality below  <dig> as well as to eliminate pcr and optical duplicates, following the protocol suggested by taiwo et al.
<cit> .
the peak calling to identify peak regions  was performed using macs <dig>   <cit> .
the a. mellifera bs-seq data set used in fig.
<dig> was processed by bismark  <cit>  using default options and bowtie  to identify 5mcs; only 5mcs covered by at least three methylated reads were included in the analysis.
additionally, to enable the a. mellifera datasets to be viewed in the integrative genomics viewer   <cit>  various bedgraph files were generated .
these files include stranded output for the bs-seq 5mc calls in all methylation contexts, the coverage of the medip reads including a separate track for the first mates and the stranded calls made by dismiss for each of the three a. mellifera developmental stages analysed.
all of these files are available  at http://uhkniazi.github.io/dismiss.
annotation
the genome and annotation information for a. mellifera   <cit>  were obtained from the ncbi.
counting overlaps of genomic features
a location in the genome can be represented as a set of start and end coordinates, a strand  and sequence name .
two locations or features are overlapping if they have the same sequence name, compatible strand pairing and overlapping coordinates.
the bioconductor package genomicranges provides this data structure to handle genomic features along with various utility functions to perform operations on these features.
in our analysis, the counting of overlapping features was performed using the bioconductor countoverlaps function:  <dig> c=countoverlapsqs
where q is the query parameter ; s is the subject parameter ; c is a set of integer values of length equal to q, giving counts of overlaps of q and s.
RESULTS
here, for the first time, we describe a tool to improve the resolution of 5mc analysis using medip-seq, which allows for the accurate determination of symmetric vs asymmetric dna methylation.
previously, the accuracy of medip-seq has been investigated via comparison with bs-seq and a good concordance in cg methylation between these technologies for symmetric methylation has been demonstrated  <cit> .
following this strategy, we demonstrate the accuracy of dismiss by performing a comparison of 5mc strand-specific predictions made on an existing medip-seq data set to 5mc sites verified in an approximately equivalent bs-seq data set, both derived from an invertebrate genome with a relatively high level of asymmetric methylation, the honeybee .
asymmetric methylation prediction: investigating potential sources of bias
dismiss was used to identify strand specific methylated dna regions in existing medip datasets derived from three adult worker life stages: nurse, forager, and reverted nurse .
depending on the life stage, macs <dig> identified approximately  <dig>  to  <dig>  methylated regions within the whole honey bee genome; dismiss identified  <dig> to 85 % of those as being single-stranded regions with the remaining  <dig> to 20 % being double-strand regions.
plots showing the distribution of strand read depth from macs <dig> peaks are given in fig.
2a.
in order to demonstrate that strand decomposition of the macs <dig> peaks by dismiss is not influenced in a negative way by read depth, the distributions of read depths for plus, minus and double stranded peaks are shown in fig.
2b.
the figure shows a clear demarcation of the stranded categories.
to investigate the possibility of bias arising in peak regions that are wide enough to include both single and double stranded methylation we performed the following tests.
based on the quantiles of the peak widths, we split the peak regions into three groups with the following widths: group 1: from  <dig> to 439 bp; group 2: from  <dig> to 770 bp; and group 3: from  <dig> to 2500 bp.
groups  <dig>   <dig> and  <dig> contain  <dig>   <dig> and 5 % of the peaks, respectively.
group  <dig> has the widest peaks and would therefore be most likely to contain a mix of double and single-stranded peaks.
in each of these groups we counted the numbers of plus stranded, minus stranded and double stranded peaks.
these counts were used in a contingency table to test the null hypothesis that the distribution of plus, minus and double stranded peaks is independent of peak width: a chi-squared test gave a p-value of  <dig> .
thus, there is a slight bias in the widest peak regions, but the trend is not statistically significant.
asymmetric methylation prediction: concordance between dismiss and bs-seq
using bs-seq data from an additional life stage , bismark ascertained that  <dig>   <dig> and 42 % of the 5mcs in cg, chg, and chh contexts were methylated, respectively.
thereafter, only regions identified by dismiss were considered in this analysis.
these regions were centred on the peaks identified by macs <dig> from the medip data and ranged from about  <dig> to 2500 bp in length with an average size of about 430 bp.
the 5mcs identified within these regions by bismark  were separated into the fractions  present on each of the plus and minus strands , and counted using eq.
<dig>  to independently validate the dismiss predictions, a p-value was calculated using a paired two-tailed t-test and was performed three times: one for each peak class from the dismiss prediction  to test the null hypothesis that the number of 5mcs on the plus and minus strands are equal.
figure 3a summarises the concordance between the two approaches within the strand-associated dna methylation regions identified by dismiss, across the whole honeybee genome.
within the dismiss predicted ‘minus stranded’ dna methylation regions, more bs-seq identified 5mcs were classified on the minus strand .
similarly, in the dismiss predicted ‘plus stranded’ dna methylation regions, more bs-seq identified 5mcs were detected on the plus strand .
within the regions predicted by dismiss to contain methylation on both strands , there were approximately equal numbers of bs-seq identified 5mcs detected on both plus and minus strands .
exemplar genomic regions selected to show instances of good concordance between bs-seq and medip for each of the three cases of strand-associated dna methylation are illustrated in fig.
3b.
the y-axis depicts the first mate read coverage from the medip-seq alignment data and the x-axis represents the 5mcs observed in bs-seq data as black dots .
for a more detailed overview of the concordance between bs-seq, medip, and dismiss at specific loci we have prepared the appropriate files  as tracks for viewing in a genome viewer such as igv   <cit>  .
single stranded dna methylation is asymmetrically distributed on sense and antisense strands
using honeybee bs-seq data, cingolani, p. et al.
have reported asymmetrical distribution of non-cg methylation at gene splice sites and have speculated about the role that this type of dna methylation plays in alternative mrna splicing  <cit> .
although we cannot match the nucleotide-level resolution of bs-seq, here we test the functionality of dismiss in identifying the presence of 5mc signals in either splice donor or splice acceptor sites in a. mellifera medip-seq data.
a genomic feature such as a splice donor site can be present on either strand of the dna duplex: the sense side is the strand on which the feature is present, which can be either the plus or minus strand and the antisense side is always the opposite strand.
the number of single stranded dna methylation signals  on the sense and antisense sides of each splice site were counted using eq.
<dig>  to independently validate the dismiss predictions, a p-value was calculated using a two tailed proportion test  <cit>  , to test the null hypothesis that the proportion of single stranded dna methylation on the sense and antisense sides of the feature are equal.
in fig.
<dig>  we show that dismiss-mediated dna methylation predictions derived from medip-seq datasets indicate significantly greater  levels of dna methylation are found on the antisense strand compared to the sense strand at genomic splice junctions in adult worker reverted-nurse samples.
this observation is consistent with those made by cingolani, p. et al.
using bs-seq  <cit> .
importantly, this is the first time that splice junction dna methylation signals have been resolved from medip-seq datasets and demonstrates that dismiss can play an important role in understanding the function of these genomic modifications during mrna splicing.
comparison of dna methylation in gene regions
the detection of dna methylation over four distinct genomic features in the honeybee genome was subsequently compared using three different analysis strategies: medip-seq data processed by macs2; medip-seq data processed by macs <dig> and strands assigned by dismiss; and bs-seq data identified 5mcs using bismark.
the four genomic features analysed were upstream loci , exons, introns, and downstream loci .
genomic features occur in the genome with different degrees of abundance, for instance there is about three times more intronic sequence than exonic sequence in the a. mellifera genome  <cit> .
in order to adjust for this variation we used eq.
<dig> to count, per  <dig> instances of a particular genomic feature, the number of overlaps between the methylation signal and that feature: the posterior distribution of dna methylation signal over the four genomic features was modelled as four gamma-distributed variables using a conjugate poisson sampling model.
in fig.
<dig> the counts per  <dig> instances have been converted to percentages, and error bars show 95 % confidence intervals that were calculated using simulation  <cit> .
fig.
<dig> shows that the distribution of detectable medip-seq signal converges towards the bs-seq signal after assigning strands to the methylated regions using dismiss.
these results demonstrate that dismiss analysis of medip-seq datasets improves the resolution of detecting stand-associated dna methylation signals  and offers a competitive alternative to bs-seq for global methylome studies.
discussion
while a variety of sequencing-based methods for genome-wide dna methylation profiling have been developed  <cit> , bisulfite sequencing  and methylated dna immunoprecipitation sequencing  strategies are the most commonly used.
bs-seq results in single nucleobase classification of symmetric  or asymmetric  dna methylation  <cit> .
nevertheless, it suffers from limitations that include incomplete bisulphite conversion of cytosines into uracils leading to false positives, possible pcr bias due to preferential amplification of methylated  or unmethylated  versions of the template, bisulfite conversion of methylated cytosines to uracils  leading to false negatives and, finally, the inability of sodium bisulfite to convert cytosines into uracils within repeat sequences forming snap-back structures  <cit> .
furthermore, the sensitivity in detecting dna methylation when using bs-seq is generally a function of sequencing depth, which makes it relatively expensive to achieve sufficient resolution for detecting 5mc in poorly methylated genomes  <cit> .
all of these reasons could explain why previously bs-seq experiments failed to detect dna methylation in drosophila melanogaster  <cit>  when liquid chromatography-tandem mass spectrometry  clearly demonstrated its presence – estimating that  <dig>  % of cytosines were methylated  <cit> .
in contrast, medip-seq has several advantages over bs-seq in characterising dna methylomes.
as medip-seq is an enrichment-based technique, methylated dna fragments are preferentially sequenced, resulting in greater sequencing depth of regions containing 5mc at an over-all lower total cost per genome  <cit> .
additionally, medip-seq can accurately detect 5mc in repeat regions and does not require the same amount of input gdna for sample analysis in comparison to bs-seq  <cit> .
therefore, if single-base 5mc resolution is not required, or the organism under study contains low dna methylation levels, then medip-seq provides an excellent choice for global dna methylome studies  <cit> .
in terms of detecting methylation signals in genomes, medip-seq and bs-seq generally show good concordance.
harris et al.
compared two bisulphite-based techniques: methylc-seq and reduced representation bisulphite sequencing, rrbs  to two enrichment-based techniques .
they demonstrated that the assessment of dna methylation in human embryonic stem cells by all four methods was 97 % concordant, using binary dna methylation calls  <cit> .
however, these techniques produced divergent results in terms of cg coverage, resolution, and quantitative accuracy.
in particular, enrichment based techniques such as medip-seq lacked precision when quantifying dna methylation levels due to the need to analyse genomic regions by averaging the number of cgs covered by variable numbers of reads in windows hundreds of base-pairs long.
whilst exploring the greater resolution afforded to medip-seq studies by dismiss, we had to focus on loci of methylation rather than the nucleotide-level resolution that is possible with bs-seq.
additionally, we were unable to find any whole genome data sets that directly compare bs-seq with medip-seq.
therefore, potentially confounding factors to our comparison between bs-seq and medip-seq include the fact that the datasets used were derived from different laboratories, populations, and honey-bee life stages .
nonetheless, despite originating from different developmental forms, overall analysis across the whole genome shows that the stranded predictions of 5mc made by dismiss on medip-seq datasets are highly correlated with 5mcs detected within the bs-seq datasets.
our results demonstrate that by utilising the strand information present in medip-seq data sets, dismiss can improve the resolution of the dna methylation signal retrieved by other software packages .
specifically, medip-seq data analysed by macs <dig> and dismiss achieves better concordance with bs-seq than data analysed by macs <dig> alone .
furthermore, by applying dismiss to medip-seq datasets, improved resolution of strand-associated dna methylation signals can be achieved revealing biological features hidden in previous analyses .
finally, as dismiss takes advantage of a standard laboratory protocol for generating medip-seq libraries and acquiring data – there is no need to modify these for implementation.
CONCLUSIONS
the standard laboratory protocol for medip-seq, together with dismiss, offer a unique way to explore dna methylation in either lowly methylated genomes where bs-seq may fail to detect methylation, or in genomes that contain relatively high amounts of asymmetric methylation.
we believe dismiss is a significant contribution to the field of epigenetic data analysis and to facilitate its use by the community, we have provided both the software, and complete galaxy  workflows for both dismiss as a stand-alone application and for dismiss embedded within a complete medip-seq data analysis.
abbreviations
bs-seq, sodium bisulfite conversion of dna followed by sequencing; medip-seq, methylated dna immunoprecipitation sequencing

 methylated dna immunoprecipitation sequencing, or medip-seq, is a common methodology used to study dna methylation profiles within plant and animal genomes. dna methylation is increasingly being recognised as playing a central role in the regulation of gene expression and chromatin structure. cytosine methylation is a particular type of dna methylation that most often refers to the reversible addition of a methyl group to the carbon- position of the cytosine pyrimidine ring, resulting in the formation of 5-methylcytosine modifications within double stranded dna. these epigenetic dna modifications can be further classified based on the nucleobase context in which 5mc is observed. double stranded dna methylation primarily occurs if the methylated cytosine occurs within cytosine-guanosine dinucleotide sequences i.e. cg context. this is due to nucleobase symmetry between sense and antisense dna strands. however, methylation in a cg context is not always double-stranded; it can also occur on just a single strand. single stranded dna methylation can additionally be detected if the 5mc is in a non-cg methylation context, i.e. a ch dinucleotide context where h = a, t or c. asymmetric dna methylation is a common phenomenon of plant genomes. for example it has been associated with methylation changes on retrotransposons during the development of pollen cells. however, for many years it was thought to be absent in mammals where much of the initial focus was on the detection of symmetric cg methylation. nevertheless, the presence of non-cg methylation in mammalian embryonic stem and induced pluripotent stem cells, as well as somatic tissues, e.g. adult brain is becoming more commonly observed. although asymmetric dna methylation might not have the same function as symmetric double-stranded methylation, there is preliminary evidence implying a link between non-cg methylation and transcriptional regulation. in contrast to symmetric cg methylation, cytosine methylation in a non-cg context is prevalent within gene bodies, particularly exons, rather than surrounding the transcriptional start site. there is substantial evidence suggesting that non-cg methylation serves a distinct function depending on the cell type. for instance, non-cg methylated regions are prevalent in undifferentiated cells such as embryonic stem cells, as well as differentiated cells such as brain and skeletal muscle. additionally, non-cg methylation is generally associated with a positive correlation to gene expression: this modification is thought to act as a repressive mark in the adult mammalian brain and appears to be crucial for neural function. dna methylation studies of other animals initially were hampered by assumptions translated from mammalian systems, such as the predomination of methylation occurring at cg dinucleotide sequences. recently, however, this view is being challenged with the honey bee, apis mellifera, becoming an increasingly important invertebrate model for studying how dna methylation affects metazoan development, behaviour and gene splicing. for example, cingolani et al. have described how biases in both software and experimental design towards studying cg methylation have led to dna methylation in other dinucleotide contexts being poorly identified and functionally underappreciated. after correcting for such biases, they were surprised to discover 5-fold more chh methylation than cg methylation in a. mellifera. they report that of all the cytosines present in cg, chg, and chh contexts in the genome,  % were methylated; of these, only 21 % of cg cytosines were symmetrically methylated with 0 % of chh and  % of chg cytosines symmetrically methylated. in total, symmetric methylation occurred in only  % of the methylated cytosine positions. thus, asymmetric methylation predominates in this invertebrate species. cingolani et al. further suggest that non-cg modifications might have a significant influence on the regulation of alternative splicing, and this highlights the need to detect and quantify strand-specific methylation in other metazoan organisms. the current medip-seq library preparation procedure preserves strand identity, meaning that medip-seq reads contain both symmetric and asymmetric strand methylation signals that are detectable. however, this information is not exploited by any of the currently available software packages. it therefore follows that the current downstream analysis methodology for medip-seq data is not optimal as asymmetric dna methylation is not detectable. the ability to study strand-specific dna methylation signals would significantly increase the utility of medip-seq and lead to further insight into the function of asymmetric dna methylation. we, therefore, have developed a new software to decompose medip-seq derived dna methylation information into individual, strand-specific signals. we have used dismiss to confidently assign strand-associated dna methylation signals in medip-seq datasets derived from a. mellifera. we have also quantified how well dismiss predicts strand-associated dna methylation in a. mellifera when compared to data obtained from a technology that is able to resolve methylation signals at the nucleotide-level, namely bs-seq. by demonstrating how dismiss can increase the resolution and sensitivity of medip-seq analysis of gdna, we present an enabling methodology that can be used to more accurately study dna methylomes. implementation strand conservation in medip-seq libraries the retrieval of strand identity during medip-seq analyses of genomes is possible due to sequencing adapter directionality. the methodological step-by-step process of retrieving strand-specific information is illustrated in fig. 1a. in part ), genomic 5mc can be represented by three possible scenarios: firstly, the cytosine on both dna strands is methylated, secondly the 5mc is found only on the plus dna strand, and lastly when the 5mc occurs only on the minus dna strand. during the first phase of library preparation, gdna is fragmented, end-repaired and the forked adapters are ligated directionally. subsequently, the gdna is denatured prior to the medip 5mc immuno-capture step, which leads to the selection of single stranded 5mc containing dna fragments with the 5′ adaptor shown as a bold solid black line – the presence of this adaptor is important for resolving the strand specific methylation signal. following the 5mc immuno-capture step, library preparation continues by performing a pcr amplification step, resulting in double stranded dna fragments. no matter which of the three scenarios the original 5mc is found in, the sequencing of first mate reads is initiated from the 5′ adaptor – and due to the previous selection of asymmetrically methylated fragments in step this creates a disparity between the numbers of c nucleobases and g nucleobases in the first fastq file. the strand origin of the original fragments may be recovered by aligning the sequenced reads to the genome sequence: reads containing 5mc from the minus strand will align in reverse-compliment mode.fig. dismiss utilises strand biases generated during medip-seq library preparation to detect strand-specific dna methylation. a outline of medip-seq library preparation. genomic 5mc can be present in three possible scenarios: with 5mc found on both dna strands, with 5mc found only on the plus dna strand and with 5mc found only on the minus dna strand. dna is fragmented, end-repaired and the forked adapters ligated directionally to the dna fragments – solid black and dashed black; dna is denatured. medip enrichment selects single-stranded 5mc containing fragments; a sampling bias arises due to the selection of strands originating from dna fragments containing 5mc on only a single strand – note the black adaptor is retained, identifying the strand origin of the selected fragment. pcr library enrichment is performed. black arrows show sequencing of first mate reads occurring from the black adaptor. dna fragments from all three scenarios are sequenced from both ends, with those originating from black and grey adaptors deposited in a single pair of first and second fastq files, respectively. reads from the minus strand will tend to align to the genome sequence in reverse-compliment mode. b c to g ratio bias detected in a medip-seq fastq file. base count per cycle plot for the first fastq file. the x-axis shows the base position in the read, while the y-axis shows the count of each nucleotide at that position. the a to t ratios do not show a bias, as is expected for a random library not selecting for either of these nucleobases. however, the c to g ratio demonstrates a clear bias due to the selection of 5mc containing fragments the disparity in the g to c ratio of the medip reads can be readily observed using quality control software for high-throughput sequencing data. for instance, software packages such as fastqc are commonly used to assess a number of read metrics, including the ‘base frequency per cycle’. in a typical random gdna library, the ratio between the counts of a and t nucleobases should be approximately equal, as should the g to c ratio. in a medip-enriched gdna library, however, the plot will clearly show different counts for c and g, as shown in the ‘base frequency per cycle plot’ of fig. 1b. this qc abnormality is more pronounced in cases where dna methylation is primarily found in an asymmetric context. medip-seq strand separation if the medip-seq library has been prepared following the methodological outline described in fig. , then this information can be exploited by dismiss to assign 5mc strand specificity to each peak called by a peak caller. thus, each dna methylation peak can be classified into one of the three classes: 1) minus stranded, 2) plus stranded or 3) peaks on both strands. the methods required to perform this classification are implemented in the dismiss software and are described in the following paragraphs. the count of medip reads mapping to a genomic location can be modelled as a poisson distribution. let y be the number of independent reads mapping to a genomic location that has previously been identified as a peak region by macs the distribution of y is assumed to be poisson with the parameter lambda: y~poiλ the reads mapped to the macs peak region are divided into two sub events representing the number of independent reads mapping to each of the two strands, yplus and yminus. these sub events are also distributed according to a binomial distribution: yplus~poiλθandyminus~poiλ1−θ, where θ is a fraction given by yplus/y. the goal of the analysis is to find out if there are integer values of λ that are plausible for both y and y the likelihood function gives the probability of seeing a poisson distributed value y, given a rate λ: lλ=pryλ for all the plausible values of λ there is a unique value of λ that maximises the probability given by eq.  called the maximum likelihood estimate of λ: Λ=argmaxλlλ in addition to finding the maximum likelihood estimate, we also want to know what other values of λ have a reasonably high likelihood in order to capture the randomness of the biological process. a likelihood set is a set of λ values that explain the data in a statistically significant manner: lsα≡λ:lλ/lΛ≥α where α ∈ is the cut-off value and we frequently use α ≈  for convenience and custom. further details of the poisson distribution and the comparison of distributions using likelihood ratios are thoroughly explained elsewhere. for each peak region with y reads, the decision to assign a class is performed in the following manner. given that yplus are the number of first mate reads aligned to the plus strand and yminus the number of first mate reads aligned to the minus strand, then a set of integer values from and including yplus to yminus that are plausible values of λ to explain the data are chosen. using this set of λ values and eq.  two likelihood sets are generated. if the intersection of these two sets is not empty, then there are values of λ that are plausible for both yplus and yminus and the peak is assigned to class ‘peaks on both strands’. otherwise the intersection of the sets is empty and the peak will be assigned as single-stranded: if yplus is greater than yminus then the peak is assigned to class ‘plus stranded’, else the peak is assigned to class ‘minus stranded’. for example, if the number of reads mapping to the peak identified by macs is given by y =  with the number of reads mapping to the two strands comprising this peak given by yplus  =  and yminus  =  then the set of plausible values for λ is given by { …, 60} and θ will be 60/100 = . using the poisson distributions for yplus and yminus and eq.  the likelihood for each value of λ in this set is calculated and eq.  is used to select two values of λ that maximise the likelihood for each plus and minus strand. furthermore, using eq.  two likelihood sets containing plausible values of λ are generated: lsplus = { …, 60} and lsminus = { …, 55}. the intersection of these two sets is not empty, hence there are common values of λ with a high likelihood for the observed number of reads yplus and yminus and the peak is assigned to both strands. data sets used and analysis procedures medip-seq data derived from a. mellifera consisted of adult worker nurse with sequence read archive number srr, adult worker forager with sra number srr and adult worker reverted-nurse with sra number srr. the bs-seq data originated from a. mellifera worker larvae with sra number srx. all data sets used in this study were retrieved using the sra software toolkit provided by the ncbi.fig. dismiss predictions of the dna strand on which the 5mc occurs are not influenced by read depth. a macs peaks. the frequency of occurrence of peaks, with a certain depth of first mate reads, on the plus and minus strands. the minimum read depth is first mate reads. b dismiss peaks. the depth of first mate reads on the plus strand is plotted against the depth of first mate reads on the minus strand. the predictions made by dismiss are indicated in grey scale and marked on the graph; they show well defined linear boundaries. the minimum read depth is first mate reads fig. dismiss interrogation of medip-seq datasets can accurately predict the dna strand where 5mc occurs. a comparative analysis of 5mc detected by dismiss interrogation of medip-seq data with 5mc detected by bs-seq is illustrated a despite originating from different developmental forms, 5mc detected by dismiss analysis of medip-seq datasets and bs-seq datasets are concordant. each horizontal bar represents a dismiss quantified peak class of medip-seq data. from all of the dismiss predicted peak regions in a class, the percentage of 5mcs on the plus and minus strands was calculated from the bs-seq dataset. each bar shows the percentage of 5mc identified by bs-seq on either plus strands or minus strands. the dismiss predicted ‘minus stranded’ and ‘plus stranded’ peaks predominantly contain bs-seq identified 5mcs on either minus or plus strands, respectively. b loci-specific examples of 5mc detected by dismiss analysis of medip-seq datasets, with the corresponding bs-seq predictions, over three genomic regions. three randomly selected regions from the a. mellifera genome, corresponding to dismiss predicted peak classes are indicated. the x-axis shows the nc_ base pair coordinates of each selected region. the first mate read coverage for each base using dismiss-assessed medip-seq alignment data is shown along the y-axis. the y-axis is divided into two parts by a central line where the dark-grey coverage above the line represents reads predicted to map to the plus strand and the light-grey coverage shows reads predicted to map to the minus strand. the black dots above the line indicate bs-seq identified 5mcs found on the plus strand, while those below the line show bs-seq identified 5mcs found on the minus strand fig. dna methylation is accurately predicted by dismiss over gene splice sites in a. mellifera. the x-axis shows two genomic regions representing intron splice donor and acceptor sites. the y-axis shows the percentage of single stranded medip-seq peaks, predicted by dismiss, on the sense and antisense sides of these regions. dismiss detected more dna methylation on the antisense sides of the splice acceptor and donor sites, compared to the sense side, which is statistically significant with p values less than fig. dismiss detection of dna methylation in medip-seq datasets approaches the resolution obtained by bs-seq analysis. a medip-seq dataset was analysed by macs and dismiss. the proportion of dna methylation found in four different genomic regions was deduced and compared to the proportion of dna methylation signal detected in a bs-seq dataset medip-seq and bs-seq data analysis the medip-seq fastq files were aligned to the a. mellifera genome using bowtie. samtools was subsequently used to remove aligned reads with a mapping quality below as well as to eliminate pcr and optical duplicates, following the protocol suggested by taiwo et al.. the peak calling to identify peak regions was performed using macs. the a. mellifera bs-seq data set used in fig.  was processed by bismark using default options and bowtie to identify 5mcs; only 5mcs covered by at least three methylated reads were included in the analysis. additionally, to enable the a. mellifera datasets to be viewed in the integrative genomics viewer various bedgraph files were generated. these files include stranded output for the bs-seq 5mc calls in all methylation contexts, the coverage of the medip reads including a separate track for the first mates and the stranded calls made by dismiss for each of the three a. mellifera developmental stages analysed. all of these files are available at http://uhkniazi.github.io/dismiss. annotation the genome and annotation information for a. mellifera were obtained from the ncbi. counting overlaps of genomic features a location in the genome can be represented as a set of start and end coordinates, a strand and sequence name. two locations or features are overlapping if they have the same sequence name, compatible strand pairing and overlapping coordinates. the bioconductor package genomicranges provides this data structure to handle genomic features along with various utility functions to perform operations on these features. in our analysis, the counting of overlapping features was performed using the bioconductor countoverlaps function: c=countoverlapsqs where q is the query parameter; s is the subject parameter; c is a set of integer values of length equal to q, giving counts of overlaps of q and s.  here, for the first time, we describe a tool to improve the resolution of 5mc analysis using medip-seq, which allows for the accurate determination of symmetric vs asymmetric dna methylation. previously, the accuracy of medip-seq has been investigated via comparison with bs-seq and a good concordance in cg methylation between these technologies for symmetric methylation has been demonstrated. following this strategy, we demonstrate the accuracy of dismiss by performing a comparison of 5mc strand-specific predictions made on an existing medip-seq data set to 5mc sites verified in an approximately equivalent bs-seq data set, both derived from an invertebrate genome with a relatively high level of asymmetric methylation, the honeybee. asymmetric methylation prediction: investigating potential sources of bias dismiss was used to identify strand specific methylated dna regions in existing medip datasets derived from three adult worker life stages: nurse, forager, and reverted nurse. depending on the life stage, macs identified approximately to methylated regions within the whole honey bee genome; dismiss identified to 85 % of those as being single-stranded regions with the remaining to 20 % being double-strand regions. plots showing the distribution of strand read depth from macs peaks are given in fig. 2a. in order to demonstrate that strand decomposition of the macs peaks by dismiss is not influenced in a negative way by read depth, the distributions of read depths for plus, minus and double stranded peaks are shown in fig. 2b. the figure shows a clear demarcation of the stranded categories. to investigate the possibility of bias arising in peak regions that are wide enough to include both single and double stranded methylation we performed the following tests. based on the quantiles of the peak widths, we split the peak regions into three groups with the following widths: group 1: from to 439 bp; group 2: from to 770 bp; and group 3: from to 2500 bp. groups and contain and 5 % of the peaks, respectively. group has the widest peaks and would therefore be most likely to contain a mix of double and single-stranded peaks. in each of these groups we counted the numbers of plus stranded, minus stranded and double stranded peaks. these counts were used in a contingency table to test the null hypothesis that the distribution of plus, minus and double stranded peaks is independent of peak width: a chi-squared test gave a p-value of. thus, there is a slight bias in the widest peak regions, but the trend is not statistically significant. asymmetric methylation prediction: concordance between dismiss and bs-seq using bs-seq data from an additional life stage, bismark ascertained that and 42 % of the 5mcs in cg, chg, and chh contexts were methylated, respectively. thereafter, only regions identified by dismiss were considered in this analysis. these regions were centred on the peaks identified by macs from the medip data and ranged from about to 2500 bp in length with an average size of about 430 bp. the 5mcs identified within these regions by bismark were separated into the fractions present on each of the plus and minus strands, and counted using eq.  to independently validate the dismiss predictions, a p-value was calculated using a paired two-tailed t-test and was performed three times: one for each peak class from the dismiss prediction to test the null hypothesis that the number of 5mcs on the plus and minus strands are equal. figure 3a summarises the concordance between the two approaches within the strand-associated dna methylation regions identified by dismiss, across the whole honeybee genome. within the dismiss predicted ‘minus stranded’ dna methylation regions, more bs-seq identified 5mcs were classified on the minus strand. similarly, in the dismiss predicted ‘plus stranded’ dna methylation regions, more bs-seq identified 5mcs were detected on the plus strand. within the regions predicted by dismiss to contain methylation on both strands, there were approximately equal numbers of bs-seq identified 5mcs detected on both plus and minus strands. exemplar genomic regions selected to show instances of good concordance between bs-seq and medip for each of the three cases of strand-associated dna methylation are illustrated in fig. 3b. the y-axis depicts the first mate read coverage from the medip-seq alignment data and the x-axis represents the 5mcs observed in bs-seq data as black dots. for a more detailed overview of the concordance between bs-seq, medip, and dismiss at specific loci we have prepared the appropriate files as tracks for viewing in a genome viewer such as igv. single stranded dna methylation is asymmetrically distributed on sense and antisense strands using honeybee bs-seq data, cingolani, p. et al. have reported asymmetrical distribution of non-cg methylation at gene splice sites and have speculated about the role that this type of dna methylation plays in alternative mrna splicing. although we cannot match the nucleotide-level resolution of bs-seq, here we test the functionality of dismiss in identifying the presence of 5mc signals in either splice donor or splice acceptor sites in a. mellifera medip-seq data. a genomic feature such as a splice donor site can be present on either strand of the dna duplex: the sense side is the strand on which the feature is present, which can be either the plus or minus strand and the antisense side is always the opposite strand. the number of single stranded dna methylation signals on the sense and antisense sides of each splice site were counted using eq.  to independently validate the dismiss predictions, a p-value was calculated using a two tailed proportion test, to test the null hypothesis that the proportion of single stranded dna methylation on the sense and antisense sides of the feature are equal. in fig.  we show that dismiss-mediated dna methylation predictions derived from medip-seq datasets indicate significantly greater levels of dna methylation are found on the antisense strand compared to the sense strand at genomic splice junctions in adult worker reverted-nurse samples. this observation is consistent with those made by cingolani, p. et al. using bs-seq. importantly, this is the first time that splice junction dna methylation signals have been resolved from medip-seq datasets and demonstrates that dismiss can play an important role in understanding the function of these genomic modifications during mrna splicing. comparison of dna methylation in gene regions the detection of dna methylation over four distinct genomic features in the honeybee genome was subsequently compared using three different analysis strategies: medip-seq data processed by macs2; medip-seq data processed by macs and strands assigned by dismiss; and bs-seq data identified 5mcs using bismark. the four genomic features analysed were upstream loci, exons, introns, and downstream loci. genomic features occur in the genome with different degrees of abundance, for instance there is about three times more intronic sequence than exonic sequence in the a. mellifera genome. in order to adjust for this variation we used eq.  to count, per instances of a particular genomic feature, the number of overlaps between the methylation signal and that feature: the posterior distribution of dna methylation signal over the four genomic features was modelled as four gamma-distributed variables using a conjugate poisson sampling model. in fig.  the counts per instances have been converted to percentages, and error bars show 95 % confidence intervals that were calculated using simulation. fig.  shows that the distribution of detectable medip-seq signal converges towards the bs-seq signal after assigning strands to the methylated regions using dismiss. these results demonstrate that dismiss analysis of medip-seq datasets improves the resolution of detecting stand-associated dna methylation signals and offers a competitive alternative to bs-seq for global methylome studies. discussion while a variety of sequencing-based methods for genome-wide dna methylation profiling have been developed, bisulfite sequencing and methylated dna immunoprecipitation sequencing strategies are the most commonly used. bs-seq results in single nucleobase classification of symmetric or asymmetric dna methylation. nevertheless, it suffers from limitations that include incomplete bisulphite conversion of cytosines into uracils leading to false positives, possible pcr bias due to preferential amplification of methylated or unmethylated versions of the template, bisulfite conversion of methylated cytosines to uracils leading to false negatives and, finally, the inability of sodium bisulfite to convert cytosines into uracils within repeat sequences forming snap-back structures. furthermore, the sensitivity in detecting dna methylation when using bs-seq is generally a function of sequencing depth, which makes it relatively expensive to achieve sufficient resolution for detecting 5mc in poorly methylated genomes. all of these reasons could explain why previously bs-seq experiments failed to detect dna methylation in drosophila melanogaster when liquid chromatography-tandem mass spectrometry clearly demonstrated its presence – estimating that  % of cytosines were methylated. in contrast, medip-seq has several advantages over bs-seq in characterising dna methylomes. as medip-seq is an enrichment-based technique, methylated dna fragments are preferentially sequenced, resulting in greater sequencing depth of regions containing 5mc at an over-all lower total cost per genome. additionally, medip-seq can accurately detect 5mc in repeat regions and does not require the same amount of input gdna for sample analysis in comparison to bs-seq. therefore, if single-base 5mc resolution is not required, or the organism under study contains low dna methylation levels, then medip-seq provides an excellent choice for global dna methylome studies. in terms of detecting methylation signals in genomes, medip-seq and bs-seq generally show good concordance. harris et al. compared two bisulphite-based techniques: methylc-seq and reduced representation bisulphite sequencing, rrbs to two enrichment-based techniques. they demonstrated that the assessment of dna methylation in human embryonic stem cells by all four methods was 97 % concordant, using binary dna methylation calls. however, these techniques produced divergent results in terms of cg coverage, resolution, and quantitative accuracy. in particular, enrichment based techniques such as medip-seq lacked precision when quantifying dna methylation levels due to the need to analyse genomic regions by averaging the number of cgs covered by variable numbers of reads in windows hundreds of base-pairs long. whilst exploring the greater resolution afforded to medip-seq studies by dismiss, we had to focus on loci of methylation rather than the nucleotide-level resolution that is possible with bs-seq. additionally, we were unable to find any whole genome data sets that directly compare bs-seq with medip-seq. therefore, potentially confounding factors to our comparison between bs-seq and medip-seq include the fact that the datasets used were derived from different laboratories, populations, and honey-bee life stages. nonetheless, despite originating from different developmental forms, overall analysis across the whole genome shows that the stranded predictions of 5mc made by dismiss on medip-seq datasets are highly correlated with 5mcs detected within the bs-seq datasets. our results demonstrate that by utilising the strand information present in medip-seq data sets, dismiss can improve the resolution of the dna methylation signal retrieved by other software packages. specifically, medip-seq data analysed by macs and dismiss achieves better concordance with bs-seq than data analysed by macs alone. furthermore, by applying dismiss to medip-seq datasets, improved resolution of strand-associated dna methylation signals can be achieved revealing biological features hidden in previous analyses. finally, as dismiss takes advantage of a standard laboratory protocol for generating medip-seq libraries and acquiring data – there is no need to modify these for implementation.  the standard laboratory protocol for medip-seq, together with dismiss, offer a unique way to explore dna methylation in either lowly methylated genomes where bs-seq may fail to detect methylation, or in genomes that contain relatively high amounts of asymmetric methylation. we believe dismiss is a significant contribution to the field of epigenetic data analysis and to facilitate its use by the community, we have provided both the software, and complete galaxy workflows for both dismiss as a stand-alone application and for dismiss embedded within a complete medip-seq data analysis. abbreviations bs-seq, sodium bisulfite conversion of dna followed by sequencing; medip-seq, methylated dna immunoprecipitation sequencing 
dna methylation is an important regulator of gene expression and chromatin structure. using medip-seq datasets derived from apis mellifera, an invertebrate species that contains more asymmetric- than symmetric- dna methylation, we demonstrate that dismiss can identify strand-specific dna methylation signals with similar accuracy as bisulfite sequencing. specifically, dismiss is able to confidently predict where dna methylation predominates in medip-seq datasets derived from a. mellifera samples. furthermore, dismiss can correctly detect the origin of dna methylation at splice site junctions in a. mellifera medip-seq datasets with a precision close to bs-seq analysis. finally, dismiss-mediated identification of dna methylation signals associated with upstream, exonic, intronic and downstream genomic loci from a. mellifera medip-seq datasets outperforms macs and closely approaches the results achieved by bs-seq. while asymmetric dna methylation is increasingly being found in growing numbers of eukaryotic species and is the predominant pattern observed in some invertebrate genomes, it has been difficult to detect in medip-seq datasets using existing software. dismiss now enables more sensitive examinations of medip-seq datasets and will be especially useful for the study of genomes containing either low levels of dna methylation or for genomes containing relatively high amounts of asymmetric methylation. while symmetric cg methylation has been relatively well-studied, asymmetric methylation in any dinucleotide context has received less attention. methylated dna immunoprecipitation sequencing is commonly used to identify regions of dna methylation in eukaryotic genomes. importantly, no currently available software for processing medip-seq reads is able to resolve these strand-specific dna methylation signals. when compared to dna methylation data derived from bs-seq analysis of a. mellifera worker larva, dismiss-mediated identification of strand-specific methylated cytosines is 80 % accurate. dna methylationmedip-seqbs-seqapis melliferadismissepigeneticsasymmetricgalaxyhttp://dx.doi.org/ /501100000268biotechnology and biological sciences research councilbb/k005448/1hoffmann karl f. issue-copyright-statement© the author 2016 here we introduce dismiss, a new software package that detects strand-associated dna methylation from existing medip-seq analyses. within medip-seq libraries, methylated cytosines can be found in both double-stranded and single-stranded genomic contexts.
 dna methylation is an important regulator of gene expression and chromatin structure. methylated dna immunoprecipitation sequencing is commonly used to identify regions of dna methylation in eukaryotic genomes. within medip-seq libraries, methylated cytosines can be found in both double-stranded and single-stranded genomic contexts. while symmetric cg methylation has been relatively well-studied, asymmetric methylation in any dinucleotide context has received less attention. importantly, no currently available software for processing medip-seq reads is able to resolve these strand-specific dna methylation signals. here we introduce dismiss, a new software package that detects strand-associated dna methylation from existing medip-seq analyses.  using medip-seq datasets derived from apis mellifera, an invertebrate species that contains more asymmetric- than symmetric- dna methylation, we demonstrate that dismiss can identify strand-specific dna methylation signals with similar accuracy as bisulfite sequencing. specifically, dismiss is able to confidently predict where dna methylation predominates in medip-seq datasets derived from a. mellifera samples. when compared to dna methylation data derived from bs-seq analysis of a. mellifera worker larva, dismiss-mediated identification of strand-specific methylated cytosines is 80 % accurate. furthermore, dismiss can correctly detect the origin of dna methylation at splice site junctions in a. mellifera medip-seq datasets with a precision close to bs-seq analysis. finally, dismiss-mediated identification of dna methylation signals associated with upstream, exonic, intronic and downstream genomic loci from a. mellifera medip-seq datasets outperforms macs and closely approaches the results achieved by bs-seq.  while asymmetric dna methylation is increasingly being found in growing numbers of eukaryotic species and is the predominant pattern observed in some invertebrate genomes, it has been difficult to detect in medip-seq datasets using existing software. dismiss now enables more sensitive examinations of medip-seq datasets and will be especially useful for the study of genomes containing either low levels of dna methylation or for genomes containing relatively high amounts of asymmetric methylation. keywords dna methylationmedip-seqbs-seqapis melliferadismissepigeneticsasymmetricgalaxyhttp://dx.doi.org/ /501100000268biotechnology and biological sciences research councilbb/k005448/1hoffmann karl f. issue-copyright-statement© the author 2016 