BACKGROUND
locus-specific dna methylation analysis is used widely in many research fields.
traditionally, sanger sequencing was used as the standard technique to quantify the methylation state of a specific bisulfite-treated locus at single nucleotide resolution.
nowadays, next-generation sequencing techniques are used for high-throughput sequencing of bisulfite polymerase chain reaction  amplicons obtaining many thousands of sequences in a single sequencing run  <cit> .
in such a way, the methylation heterogeneity of a given locus can be studied at the single molecule level.
with high-throughput sequencing of bisulfite pcr amplicons, it is possible to investigate methylation diversity in a sample by looking directly at methylation profiles  of the individual cells in a population, rather than considering a single profile where cpg methylation is analyzed as a mixture of methylated and unmethylated cpgs  <cit> .
analysis of epihaplotype diversity is applicable to fields as diverse as carcinogenesis, developmental biology and plant biology .
using this high-throughput approach, the epihaplotypes of the pool of cells that comprise the study sample can be treated as a population of haploid organisms.
when considered in this way, notions and techniques derived from other fields, such as population genetics, ecology and metagenomics can be incorporated into protocols.
in particular, several metrics, statistical methods and tools developed to analyze population structure can be easily imported and adapted for the analysis of methylation profiles generated from deep targeted sequencing.
it is, therefore, important to develop tools that are able to extract locus-specific ngs methylation data in a format that can be easily imported into already available statistical tools, and that allow a user-friendly, basic statistical interpretation of this particular kind of data.
here, we present amplimethprofiler, a pipeline for the extraction and analysis of methylation profiles at the single molecule level from deep targeted bisulfite sequencing of multiple dna regions.
this tool provides functions to demultiplex, filter and extract methylation profiles directly from fasta files.
among the various output formats that are available for the representation of methylation profile composition, amplimethprofiler provides the biological observation matrix   <cit>  format, which allows the user to directly import methylation profiles into a wide range of meta-genomics analysis software tools.
also, several core analyses of the epihaplotype population structure of input samples can be automatically performed by the pipeline using a local installation of qiime software  <cit> .
implementation
input data
amplimethprofiler  requires three types of input files: a file containing the reads from the sequencer in fasta format, a comma-separated file containing information on the sequenced regions, and a fasta file containing the reference sequences of the analyzed regions.
optionally, a file containing metadata associated with each sample can be provided to enable the tool to perform a series of basic epihaplotype based analyses  on the pipeline outcome.fig.
<dig>
amplimethprofiler workflow.
functional modules are represented as trapezes.
input and output files are represented as dashed and solid rectangles, respectively
demultiplexing and filtering
reads from targeted bisulfite sequencing of multiple regions are demultiplexed by comparing their 5′ and 3′ ends with a list of provided pcr primers.
the demultiplexing procedure is based on a user-provided percentage of similarity between the 5′ or 3′ end of a read sequence and the corresponding pcr primer sequences.
reads are filtered out if no match is found between at least one of the read ends or if, given a user-provided threshold, their length does not match.
extraction of methylation profiles
first, amplicons from targeted bisulfite sequencing are aligned to the corresponding bisulfite-converted reference sequence using the locally installed blastn program  <cit> .
then, the tool inspects the c and cpg aligned positions for each input read.
bisulfite efficiency for each aligned read is computed as the percentage of conversion of non-cpg cytosine residues  to thymine residues .
if the percentage of non-cpg deaminated c residues  over the total number of non-cpg c residues is below the given threshold, the read is discarded.
in this latter case, positions for which residues other than c or t  or gaps are found are excluded from the assay .
a user provided threshold defines the minimum percentage of reference non-cpg cytosine residues to be assayed to consider the bisulfite efficiency estimate valid; if this percentage is below the given threshold the read is discarded.
the methylation profile for each aligned read is determined by evaluating the deamination of cpg sites as a result of the bisulfite treatment.
for each cpg position in the aligned reference sequence , the corresponding position in the aligned read sequence is inspected.
if a c is found in that position, then that site is considered methylated; if a t is found, then the site is considered unmethylated; and if alignment gaps or other bases  are found, the methylation state of the cpg site is reported as uncertain .
methylation percentages for each site are then computed as the number of non-deaminated bases mapped on that site over the total number of non-ambiguously mapped positions.
the same procedure is used to compute bisulfite efficiency for all c  sites.
then, the abundance of each distinct methylation pattern is evaluated for each sample.
such reports are created by counting, for each of the possible 2ncpg epihaplotypes , the number of passing filter reads containing the pattern.
epihaplotype based analysis
a series of exploratory ehas are performed on the sample profile abundances obtained in the previous steps.
these analyses are performed starting from the biom file containing methylation profile abundances and a metadata file reporting the characteristics for each analyzed sample.
a local installation of biom tool  <cit>  and qiime software suite are employed for this purpose.
three kinds of analyses are performed to summarize sample epihaplotype composition:i) a series of summary statistics on the number of passing filter profiles in each sample are performed using the “biom summarize-table” command;
ii) a summary of samples’ taxonomic composition, computed as the abundance of profiles stratified by the number of methylated cpgs, is performed through qiime’s summarize_taxa_through_plots.py module; and
iii) a heatmap, comparing relative abundances of methylation profiles between samples, where profiles  are clustered by upgma hierarchical clustering, is created with qiime’s make_otu_heatmap.py script.
within-sample diversity , for samples and groups of samples in the study, is evaluated using qiime’s alpha_rarefaction.py workflow, which performs the following steps:generate rarefied profile abundance tables for each sample ;
compute measures of alpha diversity for each rarefied otu table ;
collate alpha diversity results ; and
generate alpha rarefaction plots .
the between-sample diversity  between all pairs of samples in the study is evaluated using qiime’s beta_diversity_through_plots.py workflow, which performs the following steps:rarefy profile abundance tables to remove sampling depth heterogeneity ;
compute beta diversity metrics  using bray–curtis dissimilarity between methylation profile abundances of samples;
run principal coordinates analysis ;
generate 3d emperor pcoa plots  and 2d pcoa plots ; and
compare distances within and between groups of samples using boxplots .
RESULTS
amplimethprofiler pipeline
the amplimethprofiler pipeline is composed of three functional modules , implemented in three python modules: preprocessfasta.py, methprofiles.py, qiime_analysis.py.
the preprocessfasta.py module generates, for each sequenced region, a quality filtered fasta file containing the reads from that region that passed filtering.
importantly, it creates a new fasta file for each analyzed region, whose entries are annotated with the id of the region and of the sample.
the methprofiles.py module runs on each demultiplexed, filtered fasta file generated by preprocessfasta.py and computes cpg methylation profile matrices, profile counts and several summary and quality statistics.
for each analyzed region, methprofile.py returns the following output files.
summary and quality statistics file
this file contains information about the number of reads that pass the filtering, the methylation percentage of each c in cpg sites, and the bisulfite efficiency for each c in non-cpg sites .fig.
<dig>
amplimethprofiler output files.
a content example of a summary and quality statistics file.
b content example of a plain text alignment file.
c content example of a methylation profiles file.
d content example of a methylation profile abundances file
alignment file
these files contain blast-aligned sequences in the standard blast xml output format and in plain text format.
the plain text format  allows the user to easily inspect alignments.
each entry of this file contains a filter-passed aligned read, represented by five rows that provide the following information:read identification, read length, experiment identification, region identification;
bisulfite efficiency, calculated as the percentage of deaminated cs  over all cs ;
alignment of the read sequence against its bisulfite-converted reference sequence.
methylation profiles file
this file contains the cpg methylation profile matrix  in which columns and rows represent cpg sites and single reads, respectively.
the methylation status of each cpg site in a read is coded  <dig> if the site is recognized as unmethylated,  <dig> if the site is recognized as methylated, and  <dig> if the methylation state could not be assessed .
row entries are reported in the same order as in the “alignment file”, and column order represents the cpg positions as they appear in the reference sequence.
each row can be considered as the cpg methylation profile of a single read and defines an epihaplotype in subsequent analyses.
profile abundance reports
these files contain counts of the occurrence of each epihaplotype in the sample.
such reports are provided in two formats: tabular and biom.
each entry of the tabular file  represents a distinct methylation profile along with the following information:□id: sample identification;
□profile: string representation of the methylation profile;
□count: number of times the profile has been found in the sample; and
□n_meths: number of methylated cytosines characterizing the profile.
the biom format is a common general-use format for representing biological samples that uses observation contingency tables.
the format is designed for general use in broad areas of comparative-omics and is based on the json format  <cit> .
methylation profile abundances are coded in the rich and sparse biom format .
the methylation profiles are coded as taxonomic units and the number of methylated cytosines constituting each profile, hereafter denoted as methylation class, is used as their first-level grouping factor in an ideal phylogeny.
importantly, biom coded files from different samples can be merged together in a single biom file using suitable ad-hoc scripts.
finally, qiime_analysis.py returns a first level of exploratory ehas on the input sample.
for each analyzed region a folder is created containing the following reports:a text summary file, containing summary statistics about the number of profiles present in the set of input samples.
in particular, the file reports the number of samples, the total number of observations  in all analyzed samples, the total read count, the table density , the summary of read counts per sample  and a detailed list of read counts per sample.
the profilesummary folder contains text reports and plots reporting the distribution of methylation classes among samples.
the file heatmap.pdf contains a heatmap representing the distribution of each distinct epihaplotype among all the input samples.
the alpha folder contains information and plots based on alpha diversity metrics for each provided sample.
five alpha diversity metrics are computed for each sample: number of different methylation profiles in the sample, shannon entropy, simpson index, chao  <dig> index and number of singletons .
such metrics are computed through a rarefaction procedure to evaluate eventual biases deriving from different sequencing depths.
the beta folder contains information and plots based on beta diversity between the provided samples.
all beta diversity analyses are based on a distance function between samples.
to achieve this, bray–curtis dissimilarity among the epihaplotype compositions of samples has been employed.
the files bray_curtis_dm.txt and bray_curtis_pc.txt contain pairwise distances among samples and principal component analysis data .
the bray_curtis_emperor_pcoa_plot and the pca folders contain principal coordinate analysis  plots in html format.
the first plot shows the first three components of the pcoa through an interactive 3d html interface and relies on an emperor browser tool, the second plot shows pcoa plots in 2d using combinations of the first three components.
finally, the dist_boxplot folder contains a series of boxplots reporting the distribution of pairwise differences within and between user defined groups of samples.
amplimethprofiler pipeline
the whole set of analyses presented above can be executed by running each module alone on each analyzed sample or runs can be pipelined together.
the amplimethprofiler.py module implements the whole flowchart described above by sequential application .
the “demultiplexing and filtering” and “extraction of methylation profiles” steps are first applied to each analyzed region and each provided sample separately.
thus, for each region, the module creates a single methylation profile abundance file, in the two formats described above, containing epihaplotype abundances for the whole set of analyzed samples.
finally, ehas of each analyzed region are carried out by this module using the biom file containing computed abundances for each sample.
case study
as a proof of concept, we report amplimethprofiler pipeline analysis of targeted deep bisulfite sequencing of a genomic region in the promoter of the ddo gene from gut tissues of three newborn mice  and three adult mice .
we analyzed the region spanning from − <dig> to −63 bp upstream of the transcription start site of the ddo- <dig> transcript .
to evaluate the methylation levels of the target region, we used a double-step pcr strategy to generate an amplicon library of bisulfite dna that could be sequenced by an illumina miseq sequencer.
in the first pcr reaction, we designed primers to generate tiled amplicons.
the 5′ ends of these primers contained overhang adapter sequences  to be used in the second pcr step to add multiplexing indices and illumina sequencing adapters.
paired-end reads from illumina miseq sequencing were merged together using the pear tool  <cit>  using as threshold a minimum of  <dig> overlapping residues, then quality filtered using as threshold a mean phread score of at least  <dig>  and finally converted to fasta format using prinseq  <cit> .
we then used the resulting fasta files as input to the amplimethprofiler pipeline using the following parameters:length ±50% compared with the reference sequence length;
at least 80% sequence similarity with the primer of the corresponding target region; and
at least 98% read bisulfite efficiency.
the whole analysis took 23 m  <dig>  s on a 2 × 6-core intel xeon x5660@ <dig>  ghz with 64 gb ram, running the ubuntu  <dig> .
<dig> lts operating system.
table  <dig> reports the characterization for each input sample, along with the number of input and passing filter reads.table  <dig> sample characteristics
the obtained methylation profile compositions  were then analyzed by the qiime_analisys.py module to describe samples by methylation class  and epihaplotype frequencies .
as expected, both analyses showed clear differences between mice at stage p <dig> and mice at stage p <dig>  the analysis also revealed that profile composition is consistent at a developmental stage in different mice.fig.
<dig> profile abundances plots.
a profile composition summary charts.
bar charts representing relative abundances of profiles grouped by number of methylated cpgs.
b heatmap representing methylation profile abundances in each sample
within-sample diversity indices were then computed by the qiime_analisys.py module through rarefaction at the minimum depth found in the pool of input samples.
figure  <dig> shows rarefaction curves computed by amplimethprofiler for five different alpha diversity metrics: observed species, shannon entropy, simpson index, chao  <dig> index and number of singletons .
alpha diversity curves are provided for each sample, as well as averages for the two groups along with the corresponding confidence intervals.fig.
<dig> alpha diversity rarefaction plots at sample level  and developmental stage level
also in this case, the analysis was able to identify differences between the two groups of mice, and in particular which phenotype  was richer in terms of epihaplotype composition.
in this case, p <dig> mice showed a more heterogeneous composition than fully developed mice.
finally, between-sample diversity was computed for the two groups of samples.
we let the tool compute distances between epihaplotype composition of input samples using bray-curtis distance.
differences in epihaplotype composition between the two developmental stages are represented by means of pcoa plots.
figure 5a reports the layout of a 3d emperor plot of the first three principal components from pcoa with colors representing the two developmental stages.
samples from the two groups clearly separate in the 3d space and also tend to cluster together.fig.
<dig> beta diversity plots.
a 3d emperor plot snapshot representing the first three principal components of the pcoa.
b from left to right are reported: bray-curtis distance boxplots of pairwise distances computed between samples from the same developmental stage, pairwise distances computed between pairs of samples from different developmental stages, distances within p <dig> mice, distances within p <dig> mice and distances between pairs of p <dig> and p <dig> mice
distance boxplots of epihaplotype composition are a useful graphical tool to validate this last statement.
in particular, the qiime_analisys.py module analyzes and summarizes distances within and between user defined groups of samples, reporting distributions of distances through a series of boxplots.
the first two boxplots of fig.
5b show how distances between pairs of samples from the same group are appreciably lower than distances between pairs from different groups.
finally, the third and fourth boxplots show that methylation profiles of p <dig> samples are on average 2-fold closer to each other compared with those of p <dig> samples.
discussion
in this manuscript we present the amplimethprofiler pipeline, a tool aimed at epihaplotype based analysis of data from targeted deep bisulfite sequencing experiments.
classic quantitative methylation analyses only consider percentages of methylation by characterizing each cpg site in a region, thus flattening the information on local conformation heterogeneity carried in the pool of analyzed amplicons.
these kinds of approaches unavoidably mask the intrinsic complexity of the local methylation patterns in each cell of an analyzed sample.
epihaplotype based approaches, on the other hand, offer the possibility to study the methylation state of a sample from a complementary point of view, namely by considering the methylation conformation of each single molecule in the pool of analyzed cells.
to perform such analyses with sufficient power for a biological sample, it is essential to analyze the methylation profiles of a very large number of sequences.
this can now be accomplished through targeted deep sequencing of bisulfite-treated dna.
analyzing the methylation conformation of single reads in a multi-clonal population, such as cells from cancer tissues, offers the possibility to track the progression of distinct methylation patterns among different pathological forms/stages.
a similar approach has been adopted in the study of driver and passenger dna mutations in various form of cancer  <cit> .
likewise, the proposed epihaplotype based approach to study methylation patterns, if applied at relevant selected genomic regions, such as promoters of cancer-related genes, could lead to the discovery of driver and passenger epi-mutations.
the approach proposed here is based on the idea that the epihaplotype composition of a sample can be considered as a biological community, were each distinct methylation profile can be studied exactly as a distinct taxonomical unit is studied in a metagenomics analysis.
in this way, several notions and metrics used in ecology and population genetics can be exploited to describe the heterogeneous methylation patterns in a population of cells from a sample and to assess the compositional differences between different samples.
for example, the diversity and distribution of methylation profiles characterizing a sample can be described with alpha diversity metrics, such as the number of different taxonomic units or the shannon entropy index.
likewise, differences among epihaplotype compositions of samples can be measured through beta diversity metrics, such as bray-curtis distance or euclidean distance.
the recent diffusion of metagenomics analyses, linked to the advent of microbiome analysis from raw dna sequencing data, was accompanied by the production of multiple bioinformatics tools for the analysis of biological communities  <cit> , as well as the development of standards to represent biological communities.
one of the most widely used formats in this field is the biom, which is recognized by the vast majority of tools for the analysis of biological communities.
in this regard, it can be useful to represent epihaplotype compositions as biological observation matrices.
in fact, this format gives the possibility to carry out ehas on methylation data by taking advantage of the already available repository of tools available for ecology and metagenomics.
the amplimethprofiler tool provides a complete analysis pipeline that, starting from fasta files containing reads from targeted bisulfite sequencing experiments, extracts methylation profiles from the input samples along with a series of exploratory analyses of their profile compositions.
it provides functions to demultiplex, filter and quality check input reads along with the classic quantitative assessment of cpg methylation percent per site.
by taking advantage of the local installation of the qiime suite, amplimethprofiler enables a series of basic exploratory analyses of the methylation profiles in the given experimental samples.
the core set of the analyses provided by amplimethprofiler were chosen to be instrumental for all studies investigating methylation patterns.
if more specific analyses are needed, the biom files produced by the tool, in combination with the vast collection of qiime scripts, enable the user to easily perform more sophisticated tasks depending on the specific experimental design.
table  <dig> presents a comparison of amplimethprofiler with state of the art tools for methylation analysis of bisulfite sequencing experiments.
in particular, several tools have been described in the literature for the analysis of bisulfite sequencing data  <cit>  but the majority of them were designed to explicitly provide quantitative measurement of methylation for each analyzed cpg site.
few of these tools provide outputs containing a direct representation of methylation profiles for each analyzed read and none provide output formats and statistical tools that are specifically designed for eha of methylation heterogeneity.table  <dig> comparison of existing software programs for bisulfite sequencing analysis
the computation and listing of epihaplotype abundances are certainly important, but, especially when the number of samples  begins to grow, it’s essential to provide biologists with statistical tools able to quantify and summarize the individual sample composition and the differences between samples.
compared to existing tools, amplimethprofiler pipeline offers two main advantages:it automatically provides a large number of statistical analyses and representations of intra- and inter-sample diversity in term of their epihaplotype composition;
it provides epihaplotype abundances in several output formats, which, in turn, are easy to import in other statistical and/or population genetics tools that are borrowed from ecology.
CONCLUSIONS
in conclusion, our tool provides an easy and user friendly way to extract and analyze the epihaplotype composition of reads from targeted bisulfite sequencing experiments.
amplimethprofiler is written in python language and requires a local installation of blast and  qiime tools.
it can be run on linux and os x platforms.
the software is open source and freely available at http://amplimethprofiler.sourceforge.net.
availability of data and materials
project name: amplimethprofiler
project home page: https://sourceforge.net/projects/amplimethprofiler.
operating system: linux, macos x.
programming language: python.
other requirements: biom  <dig> .
<dig> or higher , qiime  <dig>  or higher , biopython  <dig>  or higher, blast  <dig> .
<dig> or higher .
license: gnu gpl.
additional file
additional file 1: table s <dig>  the table reports for each possible epihaplotype the number of methylated cpg and the number of filter-passed aligned reads containing the epihaplotype in each sample.
abbreviations
biombiological observation matrix
ehaepihaplotype based analysis
ncpgnumber of cpg sites in a region
ngsnext generation sequencing
otuoperational taxonomic unit
pcaprincipal components analysis
pcoaprincipal coordinate analysis
pcrpolymerase chain reaction

 locus-specific dna methylation analysis is used widely in many research fields. traditionally, sanger sequencing was used as the standard technique to quantify the methylation state of a specific bisulfite-treated locus at single nucleotide resolution. nowadays, next-generation sequencing techniques are used for high-throughput sequencing of bisulfite polymerase chain reaction amplicons obtaining many thousands of sequences in a single sequencing run. in such a way, the methylation heterogeneity of a given locus can be studied at the single molecule level. with high-throughput sequencing of bisulfite pcr amplicons, it is possible to investigate methylation diversity in a sample by looking directly at methylation profiles of the individual cells in a population, rather than considering a single profile where cpg methylation is analyzed as a mixture of methylated and unmethylated cpgs. analysis of epihaplotype diversity is applicable to fields as diverse as carcinogenesis, developmental biology and plant biology. using this high-throughput approach, the epihaplotypes of the pool of cells that comprise the study sample can be treated as a population of haploid organisms. when considered in this way, notions and techniques derived from other fields, such as population genetics, ecology and metagenomics can be incorporated into protocols. in particular, several metrics, statistical methods and tools developed to analyze population structure can be easily imported and adapted for the analysis of methylation profiles generated from deep targeted sequencing. it is, therefore, important to develop tools that are able to extract locus-specific ngs methylation data in a format that can be easily imported into already available statistical tools, and that allow a user-friendly, basic statistical interpretation of this particular kind of data. here, we present amplimethprofiler, a pipeline for the extraction and analysis of methylation profiles at the single molecule level from deep targeted bisulfite sequencing of multiple dna regions. this tool provides functions to demultiplex, filter and extract methylation profiles directly from fasta files. among the various output formats that are available for the representation of methylation profile composition, amplimethprofiler provides the biological observation matrix format, which allows the user to directly import methylation profiles into a wide range of meta-genomics analysis software tools. also, several core analyses of the epihaplotype population structure of input samples can be automatically performed by the pipeline using a local installation of qiime software. implementation input data amplimethprofiler requires three types of input files: a file containing the reads from the sequencer in fasta format, a comma-separated file containing information on the sequenced regions, and a fasta file containing the reference sequences of the analyzed regions. optionally, a file containing metadata associated with each sample can be provided to enable the tool to perform a series of basic epihaplotype based analyses on the pipeline outcome.fig. amplimethprofiler workflow. functional modules are represented as trapezes. input and output files are represented as dashed and solid rectangles, respectively demultiplexing and filtering reads from targeted bisulfite sequencing of multiple regions are demultiplexed by comparing their 5′ and 3′ ends with a list of provided pcr primers. the demultiplexing procedure is based on a user-provided percentage of similarity between the 5′ or 3′ end of a read sequence and the corresponding pcr primer sequences. reads are filtered out if no match is found between at least one of the read ends or if, given a user-provided threshold, their length does not match. extraction of methylation profiles first, amplicons from targeted bisulfite sequencing are aligned to the corresponding bisulfite-converted reference sequence using the locally installed blastn program. then, the tool inspects the c and cpg aligned positions for each input read. bisulfite efficiency for each aligned read is computed as the percentage of conversion of non-cpg cytosine residues to thymine residues. if the percentage of non-cpg deaminated c residues over the total number of non-cpg c residues is below the given threshold, the read is discarded. in this latter case, positions for which residues other than c or t or gaps are found are excluded from the assay. a user provided threshold defines the minimum percentage of reference non-cpg cytosine residues to be assayed to consider the bisulfite efficiency estimate valid; if this percentage is below the given threshold the read is discarded. the methylation profile for each aligned read is determined by evaluating the deamination of cpg sites as a result of the bisulfite treatment. for each cpg position in the aligned reference sequence, the corresponding position in the aligned read sequence is inspected. if a c is found in that position, then that site is considered methylated; if a t is found, then the site is considered unmethylated; and if alignment gaps or other bases are found, the methylation state of the cpg site is reported as uncertain. methylation percentages for each site are then computed as the number of non-deaminated bases mapped on that site over the total number of non-ambiguously mapped positions. the same procedure is used to compute bisulfite efficiency for all c sites. then, the abundance of each distinct methylation pattern is evaluated for each sample. such reports are created by counting, for each of the possible 2ncpg epihaplotypes, the number of passing filter reads containing the pattern. epihaplotype based analysis a series of exploratory ehas are performed on the sample profile abundances obtained in the previous steps. these analyses are performed starting from the biom file containing methylation profile abundances and a metadata file reporting the characteristics for each analyzed sample. a local installation of biom tool and qiime software suite are employed for this purpose. three kinds of analyses are performed to summarize sample epihaplotype composition:i) a series of summary statistics on the number of passing filter profiles in each sample are performed using the “biom summarize-table” command; ii) a summary of samples’ taxonomic composition, computed as the abundance of profiles stratified by the number of methylated cpgs, is performed through qiime’s summarize_taxa_through_plots.py module; and iii) a heatmap, comparing relative abundances of methylation profiles between samples, where profiles are clustered by upgma hierarchical clustering, is created with qiime’s make_otu_heatmap.py script. within-sample diversity, for samples and groups of samples in the study, is evaluated using qiime’s alpha_rarefaction.py workflow, which performs the following steps:generate rarefied profile abundance tables for each sample; compute measures of alpha diversity for each rarefied otu table; collate alpha diversity results; and generate alpha rarefaction plots. the between-sample diversity between all pairs of samples in the study is evaluated using qiime’s beta_diversity_through_plots.py workflow, which performs the following steps:rarefy profile abundance tables to remove sampling depth heterogeneity; compute beta diversity metrics using bray–curtis dissimilarity between methylation profile abundances of samples; run principal coordinates analysis; generate 3d emperor pcoa plots and 2d pcoa plots; and compare distances within and between groups of samples using boxplots.  amplimethprofiler pipeline the amplimethprofiler pipeline is composed of three functional modules, implemented in three python modules: preprocessfasta.py, methprofiles.py, qiime_analysis.py. the preprocessfasta.py module generates, for each sequenced region, a quality filtered fasta file containing the reads from that region that passed filtering. importantly, it creates a new fasta file for each analyzed region, whose entries are annotated with the id of the region and of the sample. the methprofiles.py module runs on each demultiplexed, filtered fasta file generated by preprocessfasta.py and computes cpg methylation profile matrices, profile counts and several summary and quality statistics. for each analyzed region, methprofile.py returns the following output files. summary and quality statistics file this file contains information about the number of reads that pass the filtering, the methylation percentage of each c in cpg sites, and the bisulfite efficiency for each c in non-cpg sites.fig. amplimethprofiler output files. a content example of a summary and quality statistics file. b content example of a plain text alignment file. c content example of a methylation profiles file. d content example of a methylation profile abundances file alignment file these files contain blast-aligned sequences in the standard blast xml output format and in plain text format. the plain text format allows the user to easily inspect alignments. each entry of this file contains a filter-passed aligned read, represented by five rows that provide the following information:read identification, read length, experiment identification, region identification; bisulfite efficiency, calculated as the percentage of deaminated cs over all cs; alignment of the read sequence against its bisulfite-converted reference sequence. methylation profiles file this file contains the cpg methylation profile matrix in which columns and rows represent cpg sites and single reads, respectively. the methylation status of each cpg site in a read is coded if the site is recognized as unmethylated, if the site is recognized as methylated, and if the methylation state could not be assessed. row entries are reported in the same order as in the “alignment file”, and column order represents the cpg positions as they appear in the reference sequence. each row can be considered as the cpg methylation profile of a single read and defines an epihaplotype in subsequent analyses. profile abundance reports these files contain counts of the occurrence of each epihaplotype in the sample. such reports are provided in two formats: tabular and biom. each entry of the tabular file represents a distinct methylation profile along with the following information:□id: sample identification; □profile: string representation of the methylation profile; □count: number of times the profile has been found in the sample; and □n_meths: number of methylated cytosines characterizing the profile. the biom format is a common general-use format for representing biological samples that uses observation contingency tables. the format is designed for general use in broad areas of comparative-omics and is based on the json format. methylation profile abundances are coded in the rich and sparse biom format. the methylation profiles are coded as taxonomic units and the number of methylated cytosines constituting each profile, hereafter denoted as methylation class, is used as their first-level grouping factor in an ideal phylogeny. importantly, biom coded files from different samples can be merged together in a single biom file using suitable ad-hoc scripts. finally, qiime_analysis.py returns a first level of exploratory ehas on the input sample. for each analyzed region a folder is created containing the following reports:a text summary file, containing summary statistics about the number of profiles present in the set of input samples. in particular, the file reports the number of samples, the total number of observations in all analyzed samples, the total read count, the table density, the summary of read counts per sample and a detailed list of read counts per sample. the profilesummary folder contains text reports and plots reporting the distribution of methylation classes among samples. the file heatmap.pdf contains a heatmap representing the distribution of each distinct epihaplotype among all the input samples. the alpha folder contains information and plots based on alpha diversity metrics for each provided sample. five alpha diversity metrics are computed for each sample: number of different methylation profiles in the sample, shannon entropy, simpson index, chao index and number of singletons. such metrics are computed through a rarefaction procedure to evaluate eventual biases deriving from different sequencing depths. the beta folder contains information and plots based on beta diversity between the provided samples. all beta diversity analyses are based on a distance function between samples. to achieve this, bray–curtis dissimilarity among the epihaplotype compositions of samples has been employed. the files bray_curtis_dm.txt and bray_curtis_pc.txt contain pairwise distances among samples and principal component analysis data. the bray_curtis_emperor_pcoa_plot and the pca folders contain principal coordinate analysis plots in html format. the first plot shows the first three components of the pcoa through an interactive 3d html interface and relies on an emperor browser tool, the second plot shows pcoa plots in 2d using combinations of the first three components. finally, the dist_boxplot folder contains a series of boxplots reporting the distribution of pairwise differences within and between user defined groups of samples. amplimethprofiler pipeline the whole set of analyses presented above can be executed by running each module alone on each analyzed sample or runs can be pipelined together. the amplimethprofiler.py module implements the whole flowchart described above by sequential application. the “demultiplexing and filtering” and “extraction of methylation profiles” steps are first applied to each analyzed region and each provided sample separately. thus, for each region, the module creates a single methylation profile abundance file, in the two formats described above, containing epihaplotype abundances for the whole set of analyzed samples. finally, ehas of each analyzed region are carried out by this module using the biom file containing computed abundances for each sample. case study as a proof of concept, we report amplimethprofiler pipeline analysis of targeted deep bisulfite sequencing of a genomic region in the promoter of the ddo gene from gut tissues of three newborn mice and three adult mice. we analyzed the region spanning from − to −63 bp upstream of the transcription start site of the ddo- transcript. to evaluate the methylation levels of the target region, we used a double-step pcr strategy to generate an amplicon library of bisulfite dna that could be sequenced by an illumina miseq sequencer. in the first pcr reaction, we designed primers to generate tiled amplicons. the 5′ ends of these primers contained overhang adapter sequences to be used in the second pcr step to add multiplexing indices and illumina sequencing adapters. paired-end reads from illumina miseq sequencing were merged together using the pear tool using as threshold a minimum of overlapping residues, then quality filtered using as threshold a mean phread score of at least and finally converted to fasta format using prinseq. we then used the resulting fasta files as input to the amplimethprofiler pipeline using the following parameters:length ±50% compared with the reference sequence length; at least 80% sequence similarity with the primer of the corresponding target region; and at least 98% read bisulfite efficiency. the whole analysis took 23 m  s on a 2 × 6-core intel xeon x5660@  ghz with 64 gb ram, running the ubuntu. lts operating system. table  reports the characterization for each input sample, along with the number of input and passing filter reads.table sample characteristics the obtained methylation profile compositions were then analyzed by the qiime_analisys.py module to describe samples by methylation class and epihaplotype frequencies. as expected, both analyses showed clear differences between mice at stage p and mice at stage p the analysis also revealed that profile composition is consistent at a developmental stage in different mice.fig. profile abundances plots. a profile composition summary charts. bar charts representing relative abundances of profiles grouped by number of methylated cpgs. b heatmap representing methylation profile abundances in each sample within-sample diversity indices were then computed by the qiime_analisys.py module through rarefaction at the minimum depth found in the pool of input samples. figure  shows rarefaction curves computed by amplimethprofiler for five different alpha diversity metrics: observed species, shannon entropy, simpson index, chao index and number of singletons. alpha diversity curves are provided for each sample, as well as averages for the two groups along with the corresponding confidence intervals.fig. alpha diversity rarefaction plots at sample level and developmental stage level also in this case, the analysis was able to identify differences between the two groups of mice, and in particular which phenotype was richer in terms of epihaplotype composition. in this case, p mice showed a more heterogeneous composition than fully developed mice. finally, between-sample diversity was computed for the two groups of samples. we let the tool compute distances between epihaplotype composition of input samples using bray-curtis distance. differences in epihaplotype composition between the two developmental stages are represented by means of pcoa plots. figure 5a reports the layout of a 3d emperor plot of the first three principal components from pcoa with colors representing the two developmental stages. samples from the two groups clearly separate in the 3d space and also tend to cluster together.fig. beta diversity plots. a 3d emperor plot snapshot representing the first three principal components of the pcoa. b from left to right are reported: bray-curtis distance boxplots of pairwise distances computed between samples from the same developmental stage, pairwise distances computed between pairs of samples from different developmental stages, distances within p mice, distances within p mice and distances between pairs of p and p mice distance boxplots of epihaplotype composition are a useful graphical tool to validate this last statement. in particular, the qiime_analisys.py module analyzes and summarizes distances within and between user defined groups of samples, reporting distributions of distances through a series of boxplots. the first two boxplots of fig. 5b show how distances between pairs of samples from the same group are appreciably lower than distances between pairs from different groups. finally, the third and fourth boxplots show that methylation profiles of p samples are on average 2-fold closer to each other compared with those of p samples. discussion in this manuscript we present the amplimethprofiler pipeline, a tool aimed at epihaplotype based analysis of data from targeted deep bisulfite sequencing experiments. classic quantitative methylation analyses only consider percentages of methylation by characterizing each cpg site in a region, thus flattening the information on local conformation heterogeneity carried in the pool of analyzed amplicons. these kinds of approaches unavoidably mask the intrinsic complexity of the local methylation patterns in each cell of an analyzed sample. epihaplotype based approaches, on the other hand, offer the possibility to study the methylation state of a sample from a complementary point of view, namely by considering the methylation conformation of each single molecule in the pool of analyzed cells. to perform such analyses with sufficient power for a biological sample, it is essential to analyze the methylation profiles of a very large number of sequences. this can now be accomplished through targeted deep sequencing of bisulfite-treated dna. analyzing the methylation conformation of single reads in a multi-clonal population, such as cells from cancer tissues, offers the possibility to track the progression of distinct methylation patterns among different pathological forms/stages. a similar approach has been adopted in the study of driver and passenger dna mutations in various form of cancer. likewise, the proposed epihaplotype based approach to study methylation patterns, if applied at relevant selected genomic regions, such as promoters of cancer-related genes, could lead to the discovery of driver and passenger epi-mutations. the approach proposed here is based on the idea that the epihaplotype composition of a sample can be considered as a biological community, were each distinct methylation profile can be studied exactly as a distinct taxonomical unit is studied in a metagenomics analysis. in this way, several notions and metrics used in ecology and population genetics can be exploited to describe the heterogeneous methylation patterns in a population of cells from a sample and to assess the compositional differences between different samples. for example, the diversity and distribution of methylation profiles characterizing a sample can be described with alpha diversity metrics, such as the number of different taxonomic units or the shannon entropy index. likewise, differences among epihaplotype compositions of samples can be measured through beta diversity metrics, such as bray-curtis distance or euclidean distance. the recent diffusion of metagenomics analyses, linked to the advent of microbiome analysis from raw dna sequencing data, was accompanied by the production of multiple bioinformatics tools for the analysis of biological communities, as well as the development of standards to represent biological communities. one of the most widely used formats in this field is the biom, which is recognized by the vast majority of tools for the analysis of biological communities. in this regard, it can be useful to represent epihaplotype compositions as biological observation matrices. in fact, this format gives the possibility to carry out ehas on methylation data by taking advantage of the already available repository of tools available for ecology and metagenomics. the amplimethprofiler tool provides a complete analysis pipeline that, starting from fasta files containing reads from targeted bisulfite sequencing experiments, extracts methylation profiles from the input samples along with a series of exploratory analyses of their profile compositions. it provides functions to demultiplex, filter and quality check input reads along with the classic quantitative assessment of cpg methylation percent per site. by taking advantage of the local installation of the qiime suite, amplimethprofiler enables a series of basic exploratory analyses of the methylation profiles in the given experimental samples. the core set of the analyses provided by amplimethprofiler were chosen to be instrumental for all studies investigating methylation patterns. if more specific analyses are needed, the biom files produced by the tool, in combination with the vast collection of qiime scripts, enable the user to easily perform more sophisticated tasks depending on the specific experimental design. table  presents a comparison of amplimethprofiler with state of the art tools for methylation analysis of bisulfite sequencing experiments. in particular, several tools have been described in the literature for the analysis of bisulfite sequencing data but the majority of them were designed to explicitly provide quantitative measurement of methylation for each analyzed cpg site. few of these tools provide outputs containing a direct representation of methylation profiles for each analyzed read and none provide output formats and statistical tools that are specifically designed for eha of methylation heterogeneity.table comparison of existing software programs for bisulfite sequencing analysis the computation and listing of epihaplotype abundances are certainly important, but, especially when the number of samples begins to grow, it’s essential to provide biologists with statistical tools able to quantify and summarize the individual sample composition and the differences between samples. compared to existing tools, amplimethprofiler pipeline offers two main advantages:it automatically provides a large number of statistical analyses and representations of intra- and inter-sample diversity in term of their epihaplotype composition; it provides epihaplotype abundances in several output formats, which, in turn, are easy to import in other statistical and/or population genetics tools that are borrowed from ecology.  in conclusion, our tool provides an easy and user friendly way to extract and analyze the epihaplotype composition of reads from targeted bisulfite sequencing experiments. amplimethprofiler is written in python language and requires a local installation of blast and qiime tools. it can be run on linux and os x platforms. the software is open source and freely available at http://amplimethprofiler.sourceforge.net. availability of data and materials project name: amplimethprofiler project home page: https://sourceforge.net/projects/amplimethprofiler. operating system: linux, macos x. programming language: python. other requirements: biom. or higher, qiime or higher, biopython or higher, blast. or higher. license: gnu gpl. additional file additional file 1: table s the table reports for each possible epihaplotype the number of methylated cpg and the number of filter-passed aligned reads containing the epihaplotype in each sample. abbreviations biombiological observation matrix ehaepihaplotype based analysis ncpgnumber of cpg sites in a region ngsnext generation sequencing otuoperational taxonomic unit pcaprincipal components analysis pcoaprincipal coordinate analysis pcrpolymerase chain reaction 
the software is open source and freely available at http://amplimethprofiler.sourceforge.net. amplimethprofiler tool provides an easy and user friendly way to extract and analyze the epihaplotype composition of reads from targeted bisulfite sequencing experiments. cpg sites in an individual molecule may exist in a binary state and each individual dna molecule, containing a certain number of cpgs, is a combination of these states defining an epihaplotype. here we present amplimethprofiler, a python-based pipeline for the extraction and statistical epihaplotype analysis of amplicons from targeted deep bisulfite sequencing of multiple dna regions. epihaplotype based approaches, on the other hand, allow methylation profiles of cell populations to be studied at the single molecule level. currently available tools for methylation analysis lack output formats that explicitly report cpg methylation profiles at the single molecule level and that have suited statistical tools for their interpretation. dna methylationbisulfite sequencingepihaplotypeepihaplotype based analysismethylation profilesissue-copyright-statement© the author 2016 it can be run on linux and os x platforms. amplimethprofiler is written in python language and requires a local installation of blast and qiime tools.
 cpg sites in an individual molecule may exist in a binary state and each individual dna molecule, containing a certain number of cpgs, is a combination of these states defining an epihaplotype. classic quantification based approaches to study dna methylation are intrinsically unable to fully represent the complexity of the underlying methylation substrate. epihaplotype based approaches, on the other hand, allow methylation profiles of cell populations to be studied at the single molecule level. for such investigations, next-generation sequencing techniques can be used, both for quantitative and for epihaplotype analysis. currently available tools for methylation analysis lack output formats that explicitly report cpg methylation profiles at the single molecule level and that have suited statistical tools for their interpretation.  here we present amplimethprofiler, a python-based pipeline for the extraction and statistical epihaplotype analysis of amplicons from targeted deep bisulfite sequencing of multiple dna regions.  amplimethprofiler tool provides an easy and user friendly way to extract and analyze the epihaplotype composition of reads from targeted bisulfite sequencing experiments. amplimethprofiler is written in python language and requires a local installation of blast and qiime tools. it can be run on linux and os x platforms. the software is open source and freely available at http://amplimethprofiler.sourceforge.net. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords dna methylationbisulfite sequencingepihaplotypeepihaplotype based analysismethylation profilesissue-copyright-statement© the author 2016 