BACKGROUND
methylation of cytosine  in genomic dna is a well-characterized epigenetic modification involved in many cellular processes, including differentiation, genomic imprinting, x-chromosome inactivation, transposon silencing, chromosome stability, and maintenance of homeostasis.
aberrant dna methylation has been reported in a growing number of human diseases, such as cancer, developmental diseases, and metabolic disorders  <cit> .
until a decade ago, dna methylation studies only focused on small regions of the genome because of technical limitations.
however, recent advances in dna sequencing technology has made it possible to construct single-base resolution maps of 5-methylcytosine  at the genome-wide scale  <cit> .
the technology is collectively called whole-genome bisulfite sequencing  or methylome analysis, and its practical methods include methylc-seq  <cit>  and post-bisulfite adaptor tagging sequencing   <cit> .
reduced representation bisulfite sequencing  is also used for single-base resolution 5mc mapping in cpg-rich regions of the genome  <cit> .
the first complete methylome maps were constructed by methylc-seq , and then the pbat method was developed for performing amplification-free wgbs of a nanogram quantity of dna  <cit> .
with this method, methylome maps were constructed for human and mouse cells  as well as plant and fungal cells  <cit> .
wgbs is increasingly important in biology and medicine and the international human epigenome consortium  recommends wgbs as the standard method for dna methylation analysis .
the illumina hiseq platform accounts for the majority of wgbs studies that are currently under way because this technology generates the largest amount of data per run at the lowest cost per base among the high-throughput sequencers  <cit> .
like many other laboratories, we have been using the hiseq platform for our pbat-based wgbs.
in the base calling system of hiseq, the hiseq control software  locates clusters, extracts intensity, and calculates color matrix before the real-time analysis  performs base calling and quality scoring .
accurate base calling requires sequence diversity because identification of individual clusters and determination of their coordinates by hcs relies on the diversity.
thus, low sequence diversity samples, including bisulfite-converted dnas, are obviously not the best substrates for hiseq sequencing  <cit> .
in the course of our wgbs study on mouse spermatogonia  <cit> , we realized that different versions of hcs and rta installed on the hiseq system provided different global cpg methylation levels  for the same libraries.
this problem was reproduced in our system using different wgbs libraries and also in hiseq systems of other laboratories.
we found that the first read of pbat-seq and the second read of methylc-seq were affected.
thus, it appeared that the problem resides in inaccurate calling of guanine , which appears at the position corresponding to 5mc in the complementary strand.
these and other observations suggest that software updates can affect the sequence outputs from low diversity libraries such as wgbs libraries.
here we describe the details of the problem, determine which versions of hcs and rta are more reliable, and discuss our recommendations to minimize the problem.
RESULTS
different hcs and rta versions provide different cpg methylation levels
wgbs relies on bisulfite conversion of unmethylated c, but not 5mc, to uracil.
because 5mc normally occupies only a small proportion of cs, bisulfite-treated dna shows depletion of c, resulting in a low diversity sequence.
pbat-seq is designed to generate sequence reads complementary to the bisulfite-converted strand, and thus 5mc appears as g in the first read  and as c in the second read    <cit> .
in the course of our wgbs study of early postnatal mouse spermatogonia  <cit> , we realized that the global cpg methylation level determined by pbat-seq of the same library significantly changed upon hcs and rta updates of the hiseq sequencer.
therefore, we set out to examine the generality of the problem and to explore the causes.
throughout this paper, global cpg methylation levels refer to weighted levels, which take sequencing depth into account  <cit> .
we performed a series of single-end runs using three pbat libraries .
each library was sequenced multiple times  and each replicate run was performed using a lane of flow cell on a hiseq  <dig> or hiseq  <dig> sequencers .
we mapped single-end sequence reads trimmed to  <dig> bases on the human  or the mouse  reference genome and obtained 59– <dig> million uniquely mapped reads per lane .
we confirmed that different combinations of hcs and rta versions provided different global cpg methylation levels  for the same libraries .
such differences were observed even when an identical hiseq sequencer was used .
the global cpg methylation difference of 5% was not negligible because similar differences were observed in several types of mouse cells during differentiation .fig.
<dig> different cpg methylation levels obtained from identical pbat libraries using different hcs and rta versions.
a cpg methylation levels determined by single-end pbat-seq.
cpg methylation levels in 100 kb windows are shown as a box plot .
different proportions of g among the four bases in r <dig> obtained using different hcs and rta versions are shown as a line plot .
b correlation between the cpg methylation levels determined using hcs v <dig> .
<dig>  v <dig> .
<dig>  and v <dig> .
<dig>  cpg methylation values of 100-kb non-overlapping sliding windows across the autosomes are plotted with a linear regression line .
c differences between the cpg methylation levels determined using hcs v <dig> .
<dig> and v <dig> .
<dig> against the cpg density.
cpg methylation values were calculated in 100-kb non-overlapping sliding windows across the autosomes.
all 100 kb windows were grouped into nine classes according to the number of contained cpg.
d cpg methylation levels determined by paired-end pbat-seq .
cpg methylation levels in 100 kb windows are shown as a box plot .
different proportions of g in r <dig> and c in r <dig> obtained using different hcs and rta versions are shown as a line plot
for the three libraries, hcs v <dig> .
<dig> always provided the highest cpg methylation level, and hcs v <dig> .
<dig> the lowest .
compared with the other hcs versions, hcs v <dig> .
<dig> provided less gs , indicating that a decreased g count is the cause of the lower methylation levels.
hcs v <dig> .
<dig> provided a methylation level approximately 5% higher than that obtained by hcs v <dig> .
<dig>  even though the same rta version  was used .
thus, hcs and not rta was the major determinant of the observed differences.
regions with higher cpg density tended to show larger differences between hcs v <dig> .
<dig> and v <dig> .
<dig> .
furthermore, when we tried to identify partially methylated domains  in epilcs, the two versions gave very different results .
the pmds are observed in several cell types including cancer cells  <cit>  and associated with intermediate levels of methylation, specific histone modifications, nuclear lamina, and gene silencing  <cit> , showing that different hcs versions can impact biological outcomes.
we next prepared a paired-end pbat library from imr- <dig> human fibroblasts and compared the results obtained using the three hcs versions .
as mentioned above, 5mc appears as g in r <dig> and as c in r <dig> in paired-end pbat-seq .
r <dig> data showed cpg methylation differences with different hcs versions , similar to those observed by single-end pbat-seq .
interestingly, r <dig> data showed smaller differences  .
thus, r <dig> data derived by hcs v <dig> .
<dig> produced a significantly lower  methylation level than the corresponding r <dig> data .
ideally, r <dig> and r <dig> data should provide identical methylation levels.
among the hcs versions, v <dig> .
<dig> produced the closest r <dig> and r <dig> methylation levels .
we previously performed paired-end methylc-seq with human genomic dna  using hcs v <dig> .
<dig> and rta v <dig> .
<dig>    <cit> .
methylc-seq is designed to read 5mc as c in r <dig> and as g in r <dig>  <cit> , which is the opposite of pbat-seq.
the global cpg methylation level determined using r <dig> data  was lower than that determined using r <dig> data  as expected.
however, the difference between r <dig> and r <dig> methylation levels was relatively small  compared with the other paired-end wgbs cases.
we then realized that this particular paired-end methylc-seq library had contained unconverted phix dna at 50% w/w.
thus, the 50% w/w phix dna spike-in may have alleviated the problem, perhaps through increasing the sequence diversity, for this run using hcs v <dig> .
<dig>  these results suggest that r <dig> of pbat-seq and r <dig> of methylc-seq provide lower methylation levels than the other reads and that the methylation level is lower when 5mc is read as g or is higher when 5mc is read as c.
the hcs version suitable for wgbs
we then attempted to determine which version of hcs is most suitable for wgbs.
to generate substrates with known cpg methylation levels, lambda phage dna was methylated in vitro to near completion by treatment with sssi methyltransferase.
we confirmed the overall resistance of the treated dna to methylation-sensitive restriction enzyme hpaii .
furthermore, we performed bisulfite sequencing at three loci , which demonstrated  <dig> % cpg methylation .
we prepared a series of lambda dna mixtures with increasing proportions of the methylated dna  and performed paired-end pbat-seq using the three hcs versions.
we mapped the reads onto the lambda phage genome  and obtained  <dig> – <dig>  million uniquely mapped reads , with the calculated average depths of  <dig> – <dig>  per strand.
r <dig> data revealed striking differences in cpg methylation between the hcs versions, with hcs v <dig> .
<dig> showing the best performance  and v <dig> .
<dig> the poorest .
in contrast, r <dig> data showed smaller differences, and all data were close to the predetermined level .
hcs v <dig> .
<dig> provided the least differences between r <dig> and r <dig> methylation levels , consistent with the findings in human and mouse genomic dnas .
hcs v <dig> .
<dig> always provided the lowest methylation levels in r <dig> among the different versions .
these results indicate that hcs v <dig> .
<dig> is most suitable for wgbs among the three hcs versions and that hcs v <dig> .
<dig> provides methylation levels lower than the real values when 5mc is read as g. since r <dig> data obtained with different hcs versions were all close to the real values, single-end methylc-seq, where 5mc is read as c, should not be affected by the versions.fig.
<dig> observed versus predetermined cpg methylation levels of a series of mixture of unmethylated and in vitro methylated lambda dnas.
a the differences between the observed and predetermined cpg methylation levels are plotted against the predetermined cpg methylation levels for each hcs version.
r <dig> and r <dig> data from paired-end pbat-seq runs were separately analyzed.
b differences between the r <dig> and r <dig> cpg methylation levels are shown for each hcs version
quality scores assigned to 5mcs
next, we examined the quality scores assigned to the respective bases of the pbat-seq reads.
in the single-end pbat-seq data obtained from imr- <dig> human fibroblasts, mouse epilcs, and mouse spermatogonia, quality scores over  <dig>  were assigned to over 85% of the bases other than g .
however, the quality scores assigned to g greatly changed depending on the hcs versions.
in particular, hcs v <dig> .
<dig> assigned low quality scores to g  .
steep drops  in quality score were observed at gs in  <dig> % of the sequence reads containing at least one g  .
in contrast, high quality scores were consistently observed at gs in the unconverted phix phage control lane of the same flow cell .
hcs v <dig> .
<dig> assigned better quality scores to gs than hcs v <dig> .
<dig> with the same rta version  .
we did not find low score assignments to gs in our pbat-seq data generated using earlier hcs versions, including hcs v <dig> .
<dig>  v <dig> .
<dig>  and v <dig> .
<dig> , suggesting that hcs v <dig> may not have the problem in g calling.fig.
<dig> quality scores assigned to base g. a quality scores assigned to gs in the raw reads obtained using different hcs versions.
all gs were grouped into four classes according to the assigned quality score.
b examples of drops in quality score at gs.
representative sequence reads from the pbat-seq  and control  data generated using hcs v <dig> .
<dig> are shown.
gs in the imr- <dig> read are shown in red.
c quality scores assigned to the base representing 5mc  in r <dig> and r <dig> of the paired-end pbat-seq and methylc-seq using hcs v <dig> .
<dig>  in the methylc-seq, 50% w/w phix dna was spiked in   <cit>
in the paired-end pbat-seq using hcs v <dig> .
<dig> , gs in r <dig> showed lower quality scores than cs in r <dig> .
in contrast, in the paired-end methylc-seq using hcs v <dig> .
<dig>   <cit> , gs in r <dig> had lower quality scores than cs in r <dig>  even though phix dna was added at 50% w/w to confer sufficient sequence diversity .
these results showed that hcs v <dig> .
<dig> has problems in scoring the fewest g bases in low diversity samples showing depletion of gs.
the low quality scores and fewer g outputs may be linked to each other, and both are likely due to the fact that base g has the lowest fluorescence intensity among the four bases  <cit> .
effect of cluster density on wgbs
the identification of individual clusters and determination of their coordinates by hcs relies on sequence diversity because discrimination of clusters in close proximity requires different fluorescence signals in the initial cycles.
thus, it has been reported that low diversity sequencing is affected in a high cluster density range  <cit> .
to investigate the relationships between cluster density, hcs version, and cpg methylation level in wgbs, we prepared a series of dilutions of the paired-end pbat library , loaded them onto flow cell lanes, and created different cluster densities  .
the three hcs versions provided similar cpg methylation levels  at different cluster densities , indicating that lower cluster densities have little impact in this case.
however, hcs v <dig> .
<dig> assigned lower quality scores to cs in r <dig>  as the cluster density increases .
such a density-dependent decrease in quality score for c was not observed with the other hcs versions  .
also, gs of r <dig> of the same sequencing run did not show such drops in quality score .
because the phix control at a high cluster density  showed good quality scores at cs in r <dig> , we speculate that hcs v <dig> .
<dig> provides lower quality scores to the fewest bases in low-diversity r <dig> data.
furthermore, r <dig> data generated using hcs v <dig> .
<dig> at 629 k per mm <dig> provided a high global non-cpg  methylation level  for imr- <dig> human fibroblasts, which is clearly different from other data   <cit> , suggesting that r <dig> data may be less accurate.fig.
<dig> effect of cluster density on paired-end pbat-seq .
a quality scores assigned to the four bases in r <dig> generated using hcs v <dig> .
<dig> at different cluster densities.
the quality scores in r <dig> of phix control on the same flow cell are also shown.
b quality scores assigned to cs in r <dig> generated using different hcs versions at different cluster densities.
c quality scores assigned to the four bases in r <dig> and r <dig> generated using the latest hcs v <dig> .
<dig> at a modest cluster density .
the quality scores in r <dig> of phix control on the same flow cell are also shown
recently, hcs v <dig> .
<dig> and rta v <dig> .
<dig> were released from illumina.
to examine the performance of the latest versions, we performed paired-end pbat-seq on an imr- <dig> library constructed from the same dna as the above studies .
these versions provided a global cpg methylation level similar to that obtained by hcs v <dig> .
<dig>  in r <dig>  however, they assigned very low quality scores to overall r <dig> data at a modest cluster density   and quality scores over  <dig> were assigned to only  <dig> % of all bases.
wgbs data generated by new illumina systems
finally, we analyzed published paired-end wgbs data generated by new illumina systems, hiseq  <dig> and nextseq  <dig>  hcs v <dig> .x is installed on hiseq  <dig>  which uses patterned flow cell technology.
we analyzed methylc-seq data  generated by hiseq  <dig> and found that the difference between r <dig> and r <dig> global cpg methylation levels was  <dig> % in both data sets.
as expected, r <dig> produced a higher methylation level than corresponding r <dig> .
quality scores over  <dig> were assigned to approximately 50% of the gs in r <dig> .
thus, the performance of hiseq  <dig> seemed better than hcs v. <dig> .
<dig>  but it was not clear whether its performance was better than hcs v <dig> .
<dig>
we also analyzed methylc-seq data  generated by nextseq  <dig>  which uses a two-color chemistry.
calls for g are made where there is actually no signal on a flow cell.
in this run, 30% phix dna was spiked in  <cit> .
we found that the difference between r <dig> and r <dig> global cpg methylation levels was relatively small  .
quality scores over  <dig> were assigned to 73% of the gs in r <dig> .
taken together, wgbs outputs by hiseq  <dig> and nextseq  <dig> produced better results than hcs v <dig> .
<dig>  since we had no information on the software versions and cluster densities, further validation requires replicates from the same wgbs libraries.
CONCLUSIONS
in this study, we found the following regarding the use of illumina hiseq sequencers for wgbs.
hcs v <dig> .
<dig> is currently the best hcs version among the hcs v <dig> for wgbs .
this version provides cpg methylation levels closest to the real values.
it is better to avoid using hcs v <dig> .
<dig> for wgbs.
this version provides methylation levels lower than the real values  and assigns very low quality scores to g bases.
r <dig> of paired-end sequencing of low diversity libraries using hcs v <dig> .
<dig> or the latest hcs v <dig> .
<dig> is greatly affected by cluster densities.
thus, when using hcs v <dig> .
<dig> or v <dig> .
<dig>  it is better to choose single-end sequencing.
based on these findings, we suggest the following for wgbs using the hiseq platform.
the same hcs version should be used for data comparisons, whenever possible.
avoid using protocols that read 5mc as g, whenever possible.
choose single-end  sequencing for methylc-seq and rrbs.
we are currently developing a single-end pbat protocol where 5mc can be represented by c. addition of 50% w/w phix dna may alleviate the problem, but not fully.
check the quality scores of 5mc in each read.
data showing lower quality scores at 5mc appear to be less accurate, even if its overall quality score is high.
describe the versions of hcs and rta when publishing wgbs results.
it is also helpful to provide the information as metadata in databases.
this will help the users to judge whether they can be used for comparison.

 methylation of cytosine in genomic dna is a well-characterized epigenetic modification involved in many cellular processes, including differentiation, genomic imprinting, x-chromosome inactivation, transposon silencing, chromosome stability, and maintenance of homeostasis. aberrant dna methylation has been reported in a growing number of human diseases, such as cancer, developmental diseases, and metabolic disorders. until a decade ago, dna methylation studies only focused on small regions of the genome because of technical limitations. however, recent advances in dna sequencing technology has made it possible to construct single-base resolution maps of 5-methylcytosine at the genome-wide scale. the technology is collectively called whole-genome bisulfite sequencing or methylome analysis, and its practical methods include methylc-seq and post-bisulfite adaptor tagging sequencing. reduced representation bisulfite sequencing is also used for single-base resolution 5mc mapping in cpg-rich regions of the genome. the first complete methylome maps were constructed by methylc-seq, and then the pbat method was developed for performing amplification-free wgbs of a nanogram quantity of dna. with this method, methylome maps were constructed for human and mouse cells as well as plant and fungal cells. wgbs is increasingly important in biology and medicine and the international human epigenome consortium recommends wgbs as the standard method for dna methylation analysis. the illumina hiseq platform accounts for the majority of wgbs studies that are currently under way because this technology generates the largest amount of data per run at the lowest cost per base among the high-throughput sequencers. like many other laboratories, we have been using the hiseq platform for our pbat-based wgbs. in the base calling system of hiseq, the hiseq control software locates clusters, extracts intensity, and calculates color matrix before the real-time analysis performs base calling and quality scoring. accurate base calling requires sequence diversity because identification of individual clusters and determination of their coordinates by hcs relies on the diversity. thus, low sequence diversity samples, including bisulfite-converted dnas, are obviously not the best substrates for hiseq sequencing. in the course of our wgbs study on mouse spermatogonia, we realized that different versions of hcs and rta installed on the hiseq system provided different global cpg methylation levels for the same libraries. this problem was reproduced in our system using different wgbs libraries and also in hiseq systems of other laboratories. we found that the first read of pbat-seq and the second read of methylc-seq were affected. thus, it appeared that the problem resides in inaccurate calling of guanine, which appears at the position corresponding to 5mc in the complementary strand. these and other observations suggest that software updates can affect the sequence outputs from low diversity libraries such as wgbs libraries. here we describe the details of the problem, determine which versions of hcs and rta are more reliable, and discuss our recommendations to minimize the problem.  different hcs and rta versions provide different cpg methylation levels wgbs relies on bisulfite conversion of unmethylated c, but not 5mc, to uracil. because 5mc normally occupies only a small proportion of cs, bisulfite-treated dna shows depletion of c, resulting in a low diversity sequence. pbat-seq is designed to generate sequence reads complementary to the bisulfite-converted strand, and thus 5mc appears as g in the first read and as c in the second read. in the course of our wgbs study of early postnatal mouse spermatogonia, we realized that the global cpg methylation level determined by pbat-seq of the same library significantly changed upon hcs and rta updates of the hiseq sequencer. therefore, we set out to examine the generality of the problem and to explore the causes. throughout this paper, global cpg methylation levels refer to weighted levels, which take sequencing depth into account. we performed a series of single-end runs using three pbat libraries. each library was sequenced multiple times and each replicate run was performed using a lane of flow cell on a hiseq or hiseq sequencers. we mapped single-end sequence reads trimmed to bases on the human or the mouse reference genome and obtained 59– million uniquely mapped reads per lane. we confirmed that different combinations of hcs and rta versions provided different global cpg methylation levels for the same libraries. such differences were observed even when an identical hiseq sequencer was used. the global cpg methylation difference of 5% was not negligible because similar differences were observed in several types of mouse cells during differentiation.fig. different cpg methylation levels obtained from identical pbat libraries using different hcs and rta versions. a cpg methylation levels determined by single-end pbat-seq. cpg methylation levels in 100 kb windows are shown as a box plot. different proportions of g among the four bases in r obtained using different hcs and rta versions are shown as a line plot. b correlation between the cpg methylation levels determined using hcs v. v. and v. cpg methylation values of 100-kb non-overlapping sliding windows across the autosomes are plotted with a linear regression line. c differences between the cpg methylation levels determined using hcs v. and v. against the cpg density. cpg methylation values were calculated in 100-kb non-overlapping sliding windows across the autosomes. all 100 kb windows were grouped into nine classes according to the number of contained cpg. d cpg methylation levels determined by paired-end pbat-seq. cpg methylation levels in 100 kb windows are shown as a box plot. different proportions of g in r and c in r obtained using different hcs and rta versions are shown as a line plot for the three libraries, hcs v. always provided the highest cpg methylation level, and hcs v. the lowest. compared with the other hcs versions, hcs v. provided less gs, indicating that a decreased g count is the cause of the lower methylation levels. hcs v. provided a methylation level approximately 5% higher than that obtained by hcs v. even though the same rta version was used. thus, hcs and not rta was the major determinant of the observed differences. regions with higher cpg density tended to show larger differences between hcs v. and v.. furthermore, when we tried to identify partially methylated domains in epilcs, the two versions gave very different results. the pmds are observed in several cell types including cancer cells and associated with intermediate levels of methylation, specific histone modifications, nuclear lamina, and gene silencing, showing that different hcs versions can impact biological outcomes. we next prepared a paired-end pbat library from imr- human fibroblasts and compared the results obtained using the three hcs versions. as mentioned above, 5mc appears as g in r and as c in r in paired-end pbat-seq. r data showed cpg methylation differences with different hcs versions, similar to those observed by single-end pbat-seq. interestingly, r data showed smaller differences. thus, r data derived by hcs v. produced a significantly lower methylation level than the corresponding r data. ideally, r and r data should provide identical methylation levels. among the hcs versions, v. produced the closest r and r methylation levels. we previously performed paired-end methylc-seq with human genomic dna using hcs v. and rta v.. methylc-seq is designed to read 5mc as c in r and as g in r, which is the opposite of pbat-seq. the global cpg methylation level determined using r data was lower than that determined using r data as expected. however, the difference between r and r methylation levels was relatively small compared with the other paired-end wgbs cases. we then realized that this particular paired-end methylc-seq library had contained unconverted phix dna at 50% w/w. thus, the 50% w/w phix dna spike-in may have alleviated the problem, perhaps through increasing the sequence diversity, for this run using hcs v. these results suggest that r of pbat-seq and r of methylc-seq provide lower methylation levels than the other reads and that the methylation level is lower when 5mc is read as g or is higher when 5mc is read as c. the hcs version suitable for wgbs we then attempted to determine which version of hcs is most suitable for wgbs. to generate substrates with known cpg methylation levels, lambda phage dna was methylated in vitro to near completion by treatment with sssi methyltransferase. we confirmed the overall resistance of the treated dna to methylation-sensitive restriction enzyme hpaii. furthermore, we performed bisulfite sequencing at three loci, which demonstrated % cpg methylation. we prepared a series of lambda dna mixtures with increasing proportions of the methylated dna and performed paired-end pbat-seq using the three hcs versions. we mapped the reads onto the lambda phage genome and obtained – million uniquely mapped reads, with the calculated average depths of – per strand. r data revealed striking differences in cpg methylation between the hcs versions, with hcs v. showing the best performance and v. the poorest. in contrast, r data showed smaller differences, and all data were close to the predetermined level. hcs v. provided the least differences between r and r methylation levels, consistent with the findings in human and mouse genomic dnas. hcs v. always provided the lowest methylation levels in r among the different versions. these results indicate that hcs v. is most suitable for wgbs among the three hcs versions and that hcs v. provides methylation levels lower than the real values when 5mc is read as g. since r data obtained with different hcs versions were all close to the real values, single-end methylc-seq, where 5mc is read as c, should not be affected by the versions.fig. observed versus predetermined cpg methylation levels of a series of mixture of unmethylated and in vitro methylated lambda dnas. a the differences between the observed and predetermined cpg methylation levels are plotted against the predetermined cpg methylation levels for each hcs version. r and r data from paired-end pbat-seq runs were separately analyzed. b differences between the r and r cpg methylation levels are shown for each hcs version quality scores assigned to 5mcs next, we examined the quality scores assigned to the respective bases of the pbat-seq reads. in the single-end pbat-seq data obtained from imr- human fibroblasts, mouse epilcs, and mouse spermatogonia, quality scores over were assigned to over 85% of the bases other than g. however, the quality scores assigned to g greatly changed depending on the hcs versions. in particular, hcs v. assigned low quality scores to g. steep drops in quality score were observed at gs in % of the sequence reads containing at least one g. in contrast, high quality scores were consistently observed at gs in the unconverted phix phage control lane of the same flow cell. hcs v. assigned better quality scores to gs than hcs v. with the same rta version. we did not find low score assignments to gs in our pbat-seq data generated using earlier hcs versions, including hcs v. v. and v., suggesting that hcs v may not have the problem in g calling.fig. quality scores assigned to base g. a quality scores assigned to gs in the raw reads obtained using different hcs versions. all gs were grouped into four classes according to the assigned quality score. b examples of drops in quality score at gs. representative sequence reads from the pbat-seq and control data generated using hcs v. are shown. gs in the imr- read are shown in red. c quality scores assigned to the base representing 5mc in r and r of the paired-end pbat-seq and methylc-seq using hcs v. in the methylc-seq, 50% w/w phix dna was spiked in in the paired-end pbat-seq using hcs v., gs in r showed lower quality scores than cs in r. in contrast, in the paired-end methylc-seq using hcs v., gs in r had lower quality scores than cs in r even though phix dna was added at 50% w/w to confer sufficient sequence diversity. these results showed that hcs v. has problems in scoring the fewest g bases in low diversity samples showing depletion of gs. the low quality scores and fewer g outputs may be linked to each other, and both are likely due to the fact that base g has the lowest fluorescence intensity among the four bases. effect of cluster density on wgbs the identification of individual clusters and determination of their coordinates by hcs relies on sequence diversity because discrimination of clusters in close proximity requires different fluorescence signals in the initial cycles. thus, it has been reported that low diversity sequencing is affected in a high cluster density range. to investigate the relationships between cluster density, hcs version, and cpg methylation level in wgbs, we prepared a series of dilutions of the paired-end pbat library, loaded them onto flow cell lanes, and created different cluster densities. the three hcs versions provided similar cpg methylation levels at different cluster densities, indicating that lower cluster densities have little impact in this case. however, hcs v. assigned lower quality scores to cs in r as the cluster density increases. such a density-dependent decrease in quality score for c was not observed with the other hcs versions. also, gs of r of the same sequencing run did not show such drops in quality score. because the phix control at a high cluster density showed good quality scores at cs in r, we speculate that hcs v. provides lower quality scores to the fewest bases in low-diversity r data. furthermore, r data generated using hcs v. at 629 k per mm provided a high global non-cpg methylation level for imr- human fibroblasts, which is clearly different from other data, suggesting that r data may be less accurate.fig. effect of cluster density on paired-end pbat-seq. a quality scores assigned to the four bases in r generated using hcs v. at different cluster densities. the quality scores in r of phix control on the same flow cell are also shown. b quality scores assigned to cs in r generated using different hcs versions at different cluster densities. c quality scores assigned to the four bases in r and r generated using the latest hcs v. at a modest cluster density. the quality scores in r of phix control on the same flow cell are also shown recently, hcs v. and rta v. were released from illumina. to examine the performance of the latest versions, we performed paired-end pbat-seq on an imr- library constructed from the same dna as the above studies. these versions provided a global cpg methylation level similar to that obtained by hcs v. in r however, they assigned very low quality scores to overall r data at a modest cluster density and quality scores over were assigned to only % of all bases. wgbs data generated by new illumina systems finally, we analyzed published paired-end wgbs data generated by new illumina systems, hiseq and nextseq hcs v.x is installed on hiseq which uses patterned flow cell technology. we analyzed methylc-seq data generated by hiseq and found that the difference between r and r global cpg methylation levels was % in both data sets. as expected, r produced a higher methylation level than corresponding r. quality scores over were assigned to approximately 50% of the gs in r. thus, the performance of hiseq seemed better than hcs v.. but it was not clear whether its performance was better than hcs v. we also analyzed methylc-seq data generated by nextseq which uses a two-color chemistry. calls for g are made where there is actually no signal on a flow cell. in this run, 30% phix dna was spiked in. we found that the difference between r and r global cpg methylation levels was relatively small. quality scores over were assigned to 73% of the gs in r. taken together, wgbs outputs by hiseq and nextseq produced better results than hcs v. since we had no information on the software versions and cluster densities, further validation requires replicates from the same wgbs libraries.  in this study, we found the following regarding the use of illumina hiseq sequencers for wgbs. hcs v. is currently the best hcs version among the hcs v for wgbs. this version provides cpg methylation levels closest to the real values. it is better to avoid using hcs v. for wgbs. this version provides methylation levels lower than the real values and assigns very low quality scores to g bases. r of paired-end sequencing of low diversity libraries using hcs v. or the latest hcs v. is greatly affected by cluster densities. thus, when using hcs v. or v. it is better to choose single-end sequencing. based on these findings, we suggest the following for wgbs using the hiseq platform. the same hcs version should be used for data comparisons, whenever possible. avoid using protocols that read 5mc as g, whenever possible. choose single-end sequencing for methylc-seq and rrbs. we are currently developing a single-end pbat protocol where 5mc can be represented by c. addition of 50% w/w phix dna may alleviate the problem, but not fully. check the quality scores of 5mc in each read. data showing lower quality scores at 5mc appear to be less accurate, even if its overall quality score is high. describe the versions of hcs and rta when publishing wgbs results. it is also helpful to provide the information as metadata in databases. this will help the users to judge whether they can be used for comparison. 
this problem was reproduced multiple times with different wgbs libraries and likely to be associated with the low sequence diversity of bisulfite-converted dna. was greatly affected by cluster densities. software updates in the illumina hiseq platform can affect the outputs from low-diversity sequencing libraries such as wgbs libraries. methylation of cytosine in genomic dna is a well-characterized epigenetic modification involved in many cellular processes and diseases. whole-genome bisulfite sequencing, such as methylc-seq and post-bisulfite adaptor tagging sequencing, uses the power of high-throughput dna sequencers and provides genome-wide dna methylation profiles at single-base resolution. in addition, paired-end sequencing of low diversity libraries using hcs v. we have used the illumina hiseq platform for our pbat-based wgbs, and found that different versions of hiseq control software and real-time analysis installed on the system provided different global cpg methylation levels for the same libraries. to determine which version of hcs is most suitable for wgbs, we used substrates with predetermined cpg methylation levels, and found that hcs v. showed the poorest performance and provided artificially lower cpg methylation levels when 5-methylcytosine is read as guanine. thus, together with other experimental conditions, special care has to be taken on this point when cpg methylation levels are to be compared between different samples by wgbs. we found that hcs was the major determinant in the observed differences. whole-genome bisulfite sequencingdna methylationillumina hiseq platformhiseq control softwarehttp://dx.doi.org/ /501100003382core research for evolutional science and technology, japan science and technology agencygrant-in-aid for scientific research on innovative areas from the ministry of education, culture, sports, science, and technology of japan25112010sasaki hiroyuki http://dx.doi.org/ /501100001691japan society for the promotion of scienceissue-copyright-statement© the author 2017 is currently the best for wgbs among the examined hcs versions. however, the accuracy and consistency of wgbs outputs in relation to the operating conditions of high-throughput sequencers have not been explored.
 methylation of cytosine in genomic dna is a well-characterized epigenetic modification involved in many cellular processes and diseases. whole-genome bisulfite sequencing, such as methylc-seq and post-bisulfite adaptor tagging sequencing, uses the power of high-throughput dna sequencers and provides genome-wide dna methylation profiles at single-base resolution. however, the accuracy and consistency of wgbs outputs in relation to the operating conditions of high-throughput sequencers have not been explored.  we have used the illumina hiseq platform for our pbat-based wgbs, and found that different versions of hiseq control software and real-time analysis installed on the system provided different global cpg methylation levels for the same libraries. this problem was reproduced multiple times with different wgbs libraries and likely to be associated with the low sequence diversity of bisulfite-converted dna. we found that hcs was the major determinant in the observed differences. to determine which version of hcs is most suitable for wgbs, we used substrates with predetermined cpg methylation levels, and found that hcs v. is the best among the examined versions. hcs v. showed the poorest performance and provided artificially lower cpg methylation levels when 5-methylcytosine is read as guanine. in addition, paired-end sequencing of low diversity libraries using hcs v. or the latest hcs v. was greatly affected by cluster densities.  software updates in the illumina hiseq platform can affect the outputs from low-diversity sequencing libraries such as wgbs libraries. more recent versions are not necessarily the better, and hcs v. is currently the best for wgbs among the examined hcs versions. thus, together with other experimental conditions, special care has to be taken on this point when cpg methylation levels are to be compared between different samples by wgbs. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords whole-genome bisulfite sequencingdna methylationillumina hiseq platformhiseq control softwarehttp://dx.doi.org/ /501100003382core research for evolutional science and technology, japan science and technology agencygrant-in-aid for scientific research on innovative areas from the ministry of education, culture, sports, science, and technology of japan25112010sasaki hiroyuki http://dx.doi.org/ /501100001691japan society for the promotion of scienceissue-copyright-statement© the author 2017 