this modification adds a methyl group ( ch3 ) to the 5 position of cytosine in a dna sequence and is inheritable through cell division.1,2 for mammalian cells , it occurs only at cpg sites ( cytosines paired with guanines ) .
dna methylation plays an essential role for both normal and cancerous cell development36 and is closely related to significant processes , including x - chromosome inactivation , genomic imprinting , and tumor growth.710 in a genome , there are different types of methylation patterns , including hypermethylation , hypomethylation , and hemimethylation .
hypermethylation occurs when samples in one group ( eg , cancer patients ) have more methylation than the samples in another group ( eg , normal individuals ) .
hypomethylation occurs when samples in one group have less methylation than the samples in another group .
hemimethylation means that at a cpg site , only one strand of the dna is methylated ( denoted m in fig .
recent research studies1113 show that , on a number of genes , hemimethylation may occur as a same - strand cluster ( fig .
1b ) with only two cpg sites , or a different - strand cluster with more than two cpg sites ( fig .
the cluster pattern means that two or more consecutive cpg sites are methylated only on one dna strand , and not on the other strand , as shown in figure 1a .
the polarity ( or reverse ) hemimethylation clusters imply that , at the two consecutive cpg sites , the methylation patterns on the positive and negative strands are  mu and um ,  respectively , as shown in figure 1b .  mu and um  means that , at two adjacent cpg sites , the first site is hemimethylated as mu on the positive and negative strands , while the next adjacent site is hemimethylated as um on two strands , which has a reversed pattern .
the hemimethylation patterns shown in figure 1 are like the footprints of dna demethylation ( ie , methyl groups are removed ) in cancer.12 this  footprint  role exists because hemimethylation is a transitional state between being methylated and having no methylation at specific cpg sites .
therefore , the identification of hemimethylation is important for understanding both methylation events and the establishment of different methylation patterns . a major experimental limitation in hemimethylation studies has been the difficulty in obtaining methylation signals from the two complementary strands of dna molecules for all cpg sites in an entire genome .
previous hemimethylation studies can only obtain the hemimethylation data for a few genes using the traditional sanger sequencing and hairpin sequencing methods.1113 even though microarray technologies can obtain methylation levels genome - wide , they can not produce methylation signals on two dna strands separately . however ,
the next - generation sequencing ( ngs ) technology,14,15 combined with the bisulfite conversion technique ( ie , c is converted to u and then becomes t ) , makes it possible to obtain methylation signals at the cpg site level on both dna strands in an entire genome.1619 during the last several years , a number of pioneering research groups have successfully used the bisulfite - converted methylation sequencing method on either arabidopsis thaliana or human samples.16,1825 using bisulfite - converted methylation sequencing data , we can detect both the gain and loss of methylation by investigating hemimethylation patterns on two complementary strands .
nevertheless , the ngs technology produces a large amount of data.26 the quality of bisulfite sequencing data may be poor because of incomplete bisulfite conversion , genome variation , and sequencing errors.27 all of these features make processing and analyzing data challenging when identifying hemimethylation patterns . to address this challenge ,
this pipeline can identify hemimethylated cpg sites and characterize different patterns in an entire genome .
hmpl locates hemimethylated sites in each of the two different samples and then compares them . in the next section
part i ( preprocessing ) utilizes available software packages in three steps : sequencing data quality assessment , trimming , and alignment ( ie , steps 13 of the workflow ) .
this part includes data parsing and summary reports ( ie , steps 4 and 5 of the work - flow ) .
hmpl code and resource files can be downloaded from the following web link : http://hal.case.edu/~sun/hmpl/hmpl.zip .
step 1 : assess sequencing qualities using fastqc.28 fastqc is a software package for assessing sequencing qualities by generating basic and informative diagnostic plots for sequencing data .
this package provides a modular set of analyses for users to obtain a quick impression as to whether or not there are any obvious and serious problems before they start any downstream data analysis .
fastqc produces basic statistics plots and summary reports for ( 1 ) per base sequence quality , ( 2 ) sequence quality scores , ( 3 ) summary of per base sequence content , ( 4 ) per sequence gc content , ( 5 ) per base n content , ( 6 ) sequence length distribution , ( 7 ) duplicate sequences , ( 8) overrepresented sequences , ( 9 ) adapter content , ( 10 ) kmer ( or k - base ) content , and ( 11 ) per tile sequence quality .
quite often , sequencing quality is very low at the 3 end in sequencing reads , and raw reads may include adapter sequences . therefore , hmpl has included the quality trimming and adapter - trimming step . in particular , dynamic trimming ( the trim function provided in the software package brat29 ) and fixed - number - base trimming options are provided for quality trimming .
step 3 : align reads using brat - bw31 and obtain methylation ratios at all cytosine sites . after trimming
, alignment is done using brat - bw.31 after alignment , the methylation level ( or ratio ) at each cytosine ( or c ) site is obtained using the  acgt - count  function of the brat - bw package .
the  acgt - count  function provides two options for generating output files : ( 1 ) the counts of  a ,   c ,  
g ,  and  t  at each cytosine site and ( 2 ) the methylation level for each cytosine site . in the hmpl package
, we have chosen the second option of the  acgt - count  function . at each cytosine site
, the methylation level is calculated as the ratio of the count of  c  ( or the number of sequencing reads with methylated cytosine ) to the count of  c  and  t  ( or the total number of reads covering that site ) .
the  acgt - count  function produces methylation levels for positive and negative strands in two separate output files .
each output file includes the following columns for each cytosine site ( ie , each row ) : chromosome , start position , end position , total number of sequencing reads , methylation level , and dna strand (  +  or    ) .
step 4 : identify hemimethylation patterns and compare these patterns in two samples . in this step ,
a hemimethylated cpg site is defined as a singleton if its adjacent cpg sites within d bases ( eg , d = 100 , a user - specified distance ) are not hemimethylated .
a hemimethylation cluster consists of at least two cpg sites that are all hemimethylated , and any two adjacent cpg sites in this cluster are within d bases .
starting from the first cpg site on a chromosome , for each hemimethylated cpg site , hmpl checks if the next cpg site is hemimethylated .
if it is hemimethylated and is within a d - base region of the previous cpg site , these two are grouped together as a cluster and we continue to check the next ( or third ) cpg site ; otherwise , this hemimethylated site is defined as a singleton .
hemimethylation clusters may have the following patterns : ( 1 ) consecutive cpg sites hemimethylated in the same dna strand ( eg , the three cpg sites in fig .
1a ) , ( 2 ) the polarity or reverse hemimethylation cluster with only two cpg sites ( as shown in fig .
1b ) , and ( 3 ) consecutive cpg sites methylated on different strands and with more than two cpg sites , eg , the methylation of three cpg sites are muu and umm on the positive and negative strands , respectively ( as shown in fig .
hmpl can also compare hemimethylated singleton sites and clusters of two samples . because high - throughput sequencing data may include sequencing and alignment errors
first , the user may determine a coverage cutoff value b ( b > 0 , eg , b = 5 ) depending on the sequencing quality and coverage level . on each strand
, there must be at least b sequencing reads to cover a specific cpg site in order for hmpl to check whether or not this site is hemimethylated .
second , if the methylation level of a specific cpg site at one strand ( eg , the positive strand ) is larger than the cutoff value h0 ( eg , h0 = 0.9 ) , it is identified as m ( methylated ) ; if the methylation level is lower than the cutoff value l0 ( eg , l0 = 0.1 ) , it is identified as u ( unmethylated ) . using the above criteria
, hmpl defines a cpg site as mu ( ie , methylated on the positive strand and unmethylated on the negative strand ) , um , mm , or uu . in the case of identifying hemimethylation clusters
. if a dataset has poor sequencing quality and low coverage , the results of using different cutoff values may be very different .
therefore , we recommend that users select very stringent cutoff values for h0 and l0 to reduce the false positive discovery rate because of poor sequencing quality and low coverage .
for example , the user may use h0 = 0.9 or 0.95 instead of 0.8 and l0 = 0.1 , or 0.05 instead of 0.2 .
if the user s dataset has good sequencing quality and high coverage , changing the cutoff values may not affect the results significantly .
step 5 : provide genetic annotation and report findings for each hemimethylated cpg site , hmpl provides the following genetic information . 

gene : if a hemimethylated cpg site is located on a gene , hmpl will report the name of that gene.pomoter : if a hemimethylated cpg site is within a d - base long region ( eg , d = 1000 ) of the promoter of a gene ( ie , d - base before the transcription starting site of a gene ) , hmpl will report the name of that gene .
gene : if a hemimethylated cpg site is located on a gene , hmpl will report the name of that gene .
pomoter : if a hemimethylated cpg site is within a d - base long region ( eg , d = 1000 ) of the promoter of a gene ( ie , d - base before the transcription starting site of a gene ) , hmpl will report the name of that gene .
the hmpl uses raw sequencing reads ( in fastq format ) as input in step 1 and step 2 . in steps 3 , 4 , and 5 ,
more detailed information about the input and output files of hmpl can be found in the user manual , which can be downloaded from http://hal.case.edu/~sun/hmpl/hmpl.user.manual.pdf .
the preprocessing step of hmpl ( part i ) can be implemented with the following command ( the command options of pre.hmpl.pl are explained in table 1 ) .
< reference_name > [ options ] the preprocessing pipeline ties the software and source code together with the appropriate dataflow to ensure that the correct output is achieved .
users need to have perl , python 2.6 , r , fastqc , and brat - bw software installed on their system .
if users have finished the hemimethylation pre - processing pipeline and have obtained the files of combined cpg sites , they may only run the parsing analysis using the part ii of hmpl ( ie ,  parse.hmpl.pl  ) .
the usage of part ii is given below ( the command options of parse.hmpl.pl are explained in table 2 ) .
we have prepared a small dataset with 5 million sequencing reads ( a fastq * .fastq file ) , the human chromosome 22 reference sequence ( a fasta * .fa file ) , and the example scripts .
users may align these 5-million reads to chromosome 22 and then identify hemimethylated sites using both part i and part ii of hmpl .
users simply need to download the data , install the hmpl package , and then change the data path / directory in the example script accordingly to run the hmpl .
it takes about 11 minutes to run this small dataset using a linux computer with 4 gb ram .
in addition , in order for users to explore the different options and arguments of hmpl , we have prepared a file named  readme . 
this file includes example scripts of running the hmpl using different command options and also explains to the user what to expect in the screen output . the small dataset , example scripts , and the  readme  file can be downloaded from the following web link : http://hal.case.edu/~sun/hmpl/hmpl.zip . in order to show the running time of hmpl in practice , we use two human example datasets .
each dataset has ~50 million raw sequencing reads , which will be introduced in detail in the results section . using the linux server with dual quad - core 2.66 ghz xeon e5430 processor that has 4 gb ram for each core
, it takes ~4 hours to run part i ( the pre - processing part ) of the hmpl , pre.hmpl.pl , if the reference index is provided .
if the index file were not provided , it would first take about three to four additional hours to build a reference index for the whole human genome , which has about 3 billion bases ( ~3 gb data ) .
therefore , it is more efficient to first build a reference index for the alignment tool brat - bw before running the hmpl .
it requires ~19 minutes to run part ii of hmpl ( the parsing pipeline , parse.hmpl.pl ) when using uncombined input with the default coverage setting .
the uncombined input means that positive and negative strand methylation levels of all cpg sites are provided as two separate files .
if the positive and negative strand methylation data are combined , it will only take 15 minutes for hmpl to generate the results .
if users have a faster linux server or high - performance computing clusters that have more memory and computing power , it will take much less time ( eg , a couple of hours ) to run the hmpl preprocessing pipeline ( pre.hmpl.pl ) , and it will take just a few minutes to get the results of parsing and comparing two samples using the hmpl part ii ( parse .
part i ( preprocessing ) utilizes available software packages in three steps : sequencing data quality assessment , trimming , and alignment ( ie , steps 13 of the workflow ) .
this part includes data parsing and summary reports ( ie , steps 4 and 5 of the work - flow ) .
hmpl code and resource files can be downloaded from the following web link : http://hal.case.edu/~sun/hmpl/hmpl.zip .
step 1 : assess sequencing qualities using fastqc.28 fastqc is a software package for assessing sequencing qualities by generating basic and informative diagnostic plots for sequencing data .
this package provides a modular set of analyses for users to obtain a quick impression as to whether or not there are any obvious and serious problems before they start any downstream data analysis .
fastqc produces basic statistics plots and summary reports for ( 1 ) per base sequence quality , ( 2 ) sequence quality scores , ( 3 ) summary of per base sequence content , ( 4 ) per sequence gc content , ( 5 ) per base n content , ( 6 ) sequence length distribution , ( 7 ) duplicate sequences , ( 8) overrepresented sequences , ( 9 ) adapter content , ( 10 ) kmer ( or k - base ) content , and ( 11 ) per tile sequence quality .
quite often , sequencing quality is very low at the 3 end in sequencing reads , and raw reads may include adapter sequences . therefore , hmpl has included the quality trimming and adapter - trimming step . in particular , dynamic trimming ( the trim function provided in the software package brat29 ) and fixed - number - base trimming options are provided for quality trimming .
step 3 : align reads using brat - bw31 and obtain methylation ratios at all cytosine sites . after trimming
, alignment is done using brat - bw.31 after alignment , the methylation level ( or ratio ) at each cytosine ( or c ) site is obtained using the  acgt - count  function of the brat - bw package .
the  acgt - count  function provides two options for generating output files : ( 1 ) the counts of  a ,   c ,  
g ,  and  t  at each cytosine site and ( 2 ) the methylation level for each cytosine site . in the hmpl package
, we have chosen the second option of the  acgt - count  function . at each cytosine site
, the methylation level is calculated as the ratio of the count of  c  ( or the number of sequencing reads with methylated cytosine ) to the count of  c  and  t  ( or the total number of reads covering that site ) .
the  acgt - count  function produces methylation levels for positive and negative strands in two separate output files .
each output file includes the following columns for each cytosine site ( ie , each row ) : chromosome , start position , end position , total number of sequencing reads , methylation level , and dna strand (  +  or    ) .
step 4 : identify hemimethylation patterns and compare these patterns in two samples . in this step ,
a hemimethylated cpg site is defined as a singleton if its adjacent cpg sites within d bases ( eg , d = 100 , a user - specified distance ) are not hemimethylated .
a hemimethylation cluster consists of at least two cpg sites that are all hemimethylated , and any two adjacent cpg sites in this cluster are within d bases .
starting from the first cpg site on a chromosome , for each hemimethylated cpg site , hmpl checks if the next cpg site is hemimethylated .
if it is hemimethylated and is within a d - base region of the previous cpg site , these two are grouped together as a cluster and we continue to check the next ( or third ) cpg site ; otherwise , this hemimethylated site is defined as a singleton .
hemimethylation clusters may have the following patterns : ( 1 ) consecutive cpg sites hemimethylated in the same dna strand ( eg , the three cpg sites in fig .
1a ) , ( 2 ) the polarity or reverse hemimethylation cluster with only two cpg sites ( as shown in fig .
1b ) , and ( 3 ) consecutive cpg sites methylated on different strands and with more than two cpg sites , eg , the methylation of three cpg sites are muu and umm on the positive and negative strands , respectively ( as shown in fig .
hmpl can also compare hemimethylated singleton sites and clusters of two samples . because high - throughput sequencing data may include sequencing and alignment errors
first , the user may determine a coverage cutoff value b ( b > 0 , eg , b = 5 ) depending on the sequencing quality and coverage level . on each strand
, there must be at least b sequencing reads to cover a specific cpg site in order for hmpl to check whether or not this site is hemimethylated .
second , if the methylation level of a specific cpg site at one strand ( eg , the positive strand ) is larger than the cutoff value h0 ( eg , h0 = 0.9 ) , it is identified as m ( methylated ) ; if the methylation level is lower than the cutoff value l0 ( eg , l0 = 0.1 ) , it is identified as u ( unmethylated ) . using the above criteria
, hmpl defines a cpg site as mu ( ie , methylated on the positive strand and unmethylated on the negative strand ) , um , mm , or uu . in the case of identifying hemimethylation clusters
. if a dataset has poor sequencing quality and low coverage , the results of using different cutoff values may be very different .
therefore , we recommend that users select very stringent cutoff values for h0 and l0 to reduce the false positive discovery rate because of poor sequencing quality and low coverage .
for example , the user may use h0 = 0.9 or 0.95 instead of 0.8 and l0 = 0.1 , or 0.05 instead of 0.2 .
if the user s dataset has good sequencing quality and high coverage , changing the cutoff values may not affect the results significantly .
step 5 : provide genetic annotation and report findings for each hemimethylated cpg site , hmpl provides the following genetic information . 

gene : if a hemimethylated cpg site is located on a gene , hmpl will report the name of that gene.pomoter : if a hemimethylated cpg site is within a d - base long region ( eg , d = 1000 ) of the promoter of a gene ( ie , d - base before the transcription starting site of a gene ) , hmpl will report the name of that gene .
gene : if a hemimethylated cpg site is located on a gene , hmpl will report the name of that gene .
pomoter : if a hemimethylated cpg site is within a d - base long region ( eg , d = 1000 ) of the promoter of a gene ( ie , d - base before the transcription starting site of a gene ) , hmpl will report the name of that gene .
the hmpl uses raw sequencing reads ( in fastq format ) as input in step 1 and step 2 . in steps 3 , 4 , and 5 ,
more detailed information about the input and output files of hmpl can be found in the user manual , which can be downloaded from http://hal.case.edu/~sun/hmpl/hmpl.user.manual.pdf .
the preprocessing step of hmpl ( part i ) can be implemented with the following command ( the command options of pre.hmpl.pl are explained in table 1 ) .
1 < fastq_input_file > -p < prefix > -r < reference_name > [ options ] the preprocessing pipeline ties the software and source code together with the appropriate dataflow to ensure that the correct output is achieved
. users need to have perl , python 2.6 , r , fastqc , and brat - bw software installed on their system .
if users have finished the hemimethylation pre - processing pipeline and have obtained the files of combined cpg sites , they may only run the parsing analysis using the part ii of hmpl ( ie ,  parse.hmpl.pl  ) .
the usage of part ii is given below ( the command options of parse.hmpl.pl are explained in table 2 ) .
1 < input 1 > [ options ] . in order for users to test the hmpl ,
we have prepared a small dataset with 5 million sequencing reads ( a fastq * .fastq file ) , the human chromosome 22 reference sequence ( a fasta * .fa file ) , and the example scripts .
users may align these 5-million reads to chromosome 22 and then identify hemimethylated sites using both part i and part ii of hmpl .
users simply need to download the data , install the hmpl package , and then change the data path / directory in the example script accordingly to run the hmpl .
it takes about 11 minutes to run this small dataset using a linux computer with 4 gb ram .
in addition , in order for users to explore the different options and arguments of hmpl , we have prepared a file named  readme . 
this file includes example scripts of running the hmpl using different command options and also explains to the user what to expect in the screen output .
the small dataset , example scripts , and the  readme  file can be downloaded from the following web link : http://hal.case.edu/~sun/hmpl/hmpl.zip . in order to show the running time of hmpl in practice , we use two human example datasets .
each dataset has ~50 million raw sequencing reads , which will be introduced in detail in the results section . using the linux server with dual quad - core 2.66 ghz xeon e5430 processor that has 4 gb ram for each core
, it takes ~4 hours to run part i ( the pre - processing part ) of the hmpl , pre.hmpl.pl , if the reference index is provided .
if the index file were not provided , it would first take about three to four additional hours to build a reference index for the whole human genome , which has about 3 billion bases ( ~3 gb data ) .
therefore , it is more efficient to first build a reference index for the alignment tool brat - bw before running the hmpl .
it requires ~19 minutes to run part ii of hmpl ( the parsing pipeline , parse.hmpl.pl ) when using uncombined input with the default coverage setting .
the uncombined input means that positive and negative strand methylation levels of all cpg sites are provided as two separate files .
if the positive and negative strand methylation data are combined , it will only take 15 minutes for hmpl to generate the results .
if users have a faster linux server or high - performance computing clusters that have more memory and computing power , it will take much less time ( eg , a couple of hours ) to run the hmpl preprocessing pipeline ( pre.hmpl.pl ) , and it will take just a few minutes to get the results of parsing and comparing two samples using the hmpl part ii ( parse .
the hmpl pipeline can be used to compare any two samples of bisulfite sequencing data . in this paper , we demonstrate the use of hmpl using publicly available bisulfite - treated methylation sequencing datasets for cell lines mcf10a and mcf7.25 these two samples are breast cancer cell lines . because a number of hypermethylated genes have been reported for breast cancer cells35 and there are hemimethylation patterns reported in individual genes,1113 it is very likely that many cpg sites have been hemimethylated in these two cell lines .
the bisulfite methylation sequencing data of mcf10a and mcf7 and more information about these two cell lines can be found from the corresponding references.25,36 the sequencing reads of mcf10a and mcf7 are generated using the reduced representative bisulfite sequencing ( rrbs ) protocol.16 there are 54,295,326 and 50,054,248 sequencing reads for mcf10a and mcf7 , respectively , and the read length is 50-base for each dataset . in order to identify hemimethylation singletons and clusters in both mcf10a and mcf7 and compare these two samples
, we have run both part i ( pre.hmpl.pl ) and part ii ( parse.hmpl.pl ) of the hmpl . in this section
, we mainly focus on showing the hemimethylation result , which is the summary of the hmpl part ii ( parse.hmpl.pl ) output . for a detailed description of the output files ,
see table 3 . the hemimethylated singleton and cluster patterns are compared and summarized in figure 3 . in this figure , 
i. singleton  means comparing the singleton hemimethylated cpg sites in mcf10a and mcf7 . 
consecutive polarity cluster  shows the results of comparing the polarity ( or reverse ) clusters that include two consecutive cpg sites ; no other cpg sites are located between these two sites .
non - consecutive polarity cluster  means comparing the polarity ( or reverse ) clusters that include two cpg sites that are not consecutive
. there is at least one cpg site located between these two sites , but there are either no sequencing reads ( or data ) or no hemimethylation sites between them .  iv . non - polarity cluster
 refers to the clusters that do not have the polarity ( or reverse ) pattern ( eg , the patterns shown in fig .
figure 3 and table 4 show the results of comparing mcf10a with mcf7 , which are summarized based on the output files named  * compare  as shown in the last row of table 3 .
the comparison results indicate that the hemimethylation patterns between the non - tumorigenic sample mcf10a and tumorigenic sample mcf7 are different at some cpg sites and/or genomic regions .
our pipeline hmpl has provided gene annotation files for all cpg sites by giving names of genes in which hemimethylated cpg sites are located , as well as the names of genes in whose promoter regions the hemimethylated cpg sites are located .
the results of comparing mcf10a and mcf7 show that there are more polarity ( or reverse ) clusters ( eg , fig .
1b ) than the single - strand hemimethylation clusters ( eg , fig . 1a ) .
this may be as a result of the fact that the methylation sequencing data we have used are generated by the rrbs protocol,16 which only sequences a small percentage of cpg sites in a human genome . in fact , there are ~5% of the cpg sites with at least 3 coverage in the rrbs data we have analyzed .
if we use the whole genome bisulfite sequencing ( wgbs ) data , it is very likely that more single - strand hemimethylation clusters would be identified .
currently , we are not aware of any wgbs data for either mcf10a or mcf7 . therefore , we have used the available rrbs data to demonstrate the usage of hmpl .
even though rrbs data are not ideal for identifying all hemimethylation sites in an entire genome , we have found many hemimethylation singletons and clusters , which show the capability of our hmpl package .
in fact , hmpl can be used to compare any two samples with data generated using either the rrbs or wgbs protocol .
in order to see if the genes with hemimethylated cpg sites are biologically important or meaningful , we have further investigated the 532 genes by comparing them with oncogenes , breast cancer methylated genes , and transcription factors .
the comparison results show that seven of these genes are methylated , 17 are oncogenes , and 62 are transcription factors .
we have also conducted the gene set enrichment analysis ( gsea ) for these 532 genes using the gsea software package and the molecular signature database provided by the broad institute.37 the analysis results show that 87 genes are significantly represented in ( or overlapped with ) 10 cancer modules ( with p - value < 0.05 ) , which are gene sets that are significantly changed in a variety of cancer conditions . these 87 genes and the 10 cancer modules they belong to are provided in table 5 . a detailed description of these 10 modules can be found online.38 the 532-gene list is included in the hmpl.zip file that can be downloaded from the following web link : http://hal.case.edu/~sun/hmpl/hmpl.zip .
there are a number of alignment tools for bisulfite methylation sequencing data , such as brat,29 brat - bw,31 bsmap,39 bs seeker,40 bismark,41 methylcoder,42 rmapbs,43 pash,44 and batmeth.45 a comprehensive list of these tools can be found at omictools.com.46 among all these available tools , we have tested brat - bw , brat , and bsmap .
we have found that they provide similar results , but brat - bw is faster .
brat - bw is user - friendly and has the following useful features : ( 1 ) it can align both single - end and paired - end reads , ( 2 ) it can produce the acgt count for all cytosines in a genome , ( 3 ) it can account for overlapping paired - end reads , and ( 4 ) it can check strands . if users prefer another alignment tool , they can obtain the alignment results and methylation levels using their preferred alignment tool , and then run part ii of the hmpl to obtain hemimethylated singletons and clusters . if necessary , this can be done with some minor format changes of their alignment output files .
reformatting is easy because the parsing pipeline of hmpl ( ie , part ii ) only requires data with the following columns that most alignment tools provide for two dna strands : chromosome , start position , end position , total number of sequencing reads , methylation level , and dna strand .
the cutoff values used in hmpl , especially the ones provided in table 2 for the parsing pipeline ( ie , part ii ) , should be determined depending on the sequencing quality and coverage .
if the user has a sequencing dataset with very good quality and high coverage , changing the cutoff value may not significantly affect the results .
however , if the sequencing dataset has low coverage and poor quality , changing the cutoff values may lead to very different results . for this case
, we recommend that users set up very stringent cutoff values to reduce false discovery rates .
we set up the default values based on findings in previous publications , some basic and common knowledge of methylation sequencing data , and our experience with bisulfite sequencing data ; however , every dataset is different .
we suggest that the users first determine the quality and coverage of their data , and then try different values to see if their results are dramatically different .
if results are dramatically different , users may choose results that are obtained based on stringent cutoff values , especially when they plan to do experimental validation .
using results based on stringent cutoff values can ensure a high validation rate and then the user may expand their validation list by adding more hemimethylated singletons or clusters from the list obtained with less stringent cutoff values .
ideally , it is best to have a list of known hemimethylated and nonhemimethylated sites to study the true and false discovery rates ( or sensitivity and specificity ) of hmpl . for the nonhemimethylated sites
, we can choose the cpg sites that are located in the housekeeping genes , which are relatively stable and not likely to be methylated on either strand .
in fact , housekeeping genes have been used for the purpose of  negative control  in previous methylation studies.4749 for the example dataset mcf7 , using the coverage cutoff of 5 ( that is , at least five reads to cover each strand ) and other default settings , there are 532 genes with at least three hemimethylated sites .
we compare these 532 genes with the 205 known housekeeping genes used in a previous study,49 and there is no overlap ; therefore , none of these 532 genes are housekeeping genes . as for the known hemimethylated sites , to the best of our knowledge
, no available list of genes or sites can be used as a  positive control . 
this is because genome - wide hemimethylation study is rarely done , and a few hemimethylated sites have been experimentally validated .
however , we have decided not to use this approach because there is little knowledge about genome - wide hemimethylation patterns . with little known information
, a simulation would be very arbitrary and the simulated data would not reflect true unknown patterns .
the purpose of our hmpl development is to provide an exploratory tool for this research topic . for genes identified with a number of hemimethylated cpg sites , users may do further investigation by studying their biological functions and relationships with other genes using pathway analyses .
first , the hmpl is designed for identifying hemimethylation only at cpg sites , but not at any non - cpg sites ( eg , chg and chh sites , where h represents a , c , or t in a dna sequence ) .
if users are interested , our algorithms may be modified to study the hemimethylation of non - cpg sites .
second , our pipeline is developed for identifying hemimethylated sites ( or clusters ) by comparing two samples , but it is not designed for comparing multiple samples in two or more groups . because the hemimethylation study is an area with little research , hmpl is good for preliminary studies . as for the topic of identifying hemimethylation patterns in multiple samples and in multiple groups , more sophisticated statistical methods may be used , and our group is working on projects related to this approach .
for the first step of the hmpl workflow , we have used the software package fastqc . even though fastqc is not designed for bisulfite - treated methylation sequencing data
if users find some serious sequencing quality issues in their data , we recommend that they check data more thoroughly using other available software packages , such as saap - rrbs,50 bseqc,51 and methyqa,27 before they interpret their hmpl results .
therefore , it is important to develop a software package to identify such patterns . to address this need ,
we have developed a new software package , hmpl , which includes both preprocessing and data parsing . for two samples , each with 50 million reads , it takes a few hours for hmpl to align the sequencing reads , and it only takes a few minutes to process the methylation level data .
if users have obtained their coverage and methylation ratio data , part ii of hmpl can identify hemimethylation patterns in minutes .
 dna methylation ( the addition of a methyl group to a cytosine ) is an important epigenetic event in mammalian cells because it plays a key role in regulating gene expression . 
 most previous methylation studies assume that dna methylation occurs on both positive and negative strands . 
 however , a few studies have reported that in some genes , methylation occurs only on one strand ( ie , hemimethylation ) and has clustering patterns . 
 these studies report that hemimethylation occurs on individual genes . 
 it is unclear whether hemimethylation occurs genome - wide and whether there are hemimethylation differences between cancerous and noncancerous cells . to address these questions 
 , we have developed the first - ever pipeline , named hemimethylation pipeline ( hmpl ) , to identify hemimethylation patterns . utilizing the available software and the newly developed perl and r scripts 
 , hmpl can identify hemimethylation patterns for a single sample and can also compare two different samples . 
Introduction
Methods
The workflow of HMPL
Input and output
Usage, command options, and running time
Results
Discussion
Conclusion