Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Optimized design and data analysis of tag-based cytosine methylation assays

Optimized design and data analysis of tag-based cytosine methylation assays Using the type III restriction-modification enzyme EcoP15I, we isolated sequences flanking sites digested by the methylation-sensitive HpaII enzyme or its methylation-insensitive MspI isoschizomer for massively parallel sequencing. A novel data transformation allows us to normalise HpaII by MspI counts, resulting in more accurate quantification of methylation at >1.8 million loci in the human genome. This HELP-tagging assay is not sensitive to sequence polymorphism or base composition and allows exploration of both CG-rich and depleted genomic contexts. Background starting material, due to the resistance of methylcytosine Epigenetic mechanisms of transcriptional regulation are to bisulfite conversion compared with unmethylated increasingly being studied for their potential influences in cytosines. This allows nucleotide resolution, strand-spe- human disease pathogenesis. Much of this interest is cific, quantitative assessment of cytosine methylation, based on the paradigm of neoplastic transformation, in with such studies performed in Arabidopsis [3-5] and which epigenetic changes appear to be universal, wide- human cells to date [6]. spread throughout the genome, causative of critical tran- While this approach represents the ideal means of test- scriptional changes and predictive of disease prognosis ing cytosine methylation, the amount of sequencing nec- (reviewed in [1]). Furthermore, these epigenetic changes essary (for the human genome, over 1 billion sequences of represent potential pharmacological targets for reversal ~75 bp each [6]) to generate quantitative information and amelioration of the disease process [2]. genome-wide remains prohibitive in terms of cost, limit- Of the large number of regulatory processes referred to ing these studies to the few referred to above. When as epigenetic, there exist numerous assays to study chro- studying human disease, the emphasis remains on cyto- matin component distribution, cytosine methylation and sine methylation assays, as it is generally easier to collect microRNA expression genome-wide. The chromatin clinical samples for DNA purification than for ChIP or components include a large number of post-translational even RNA assays. However, the cell populations har- modifications of histones, variant histones, DNA-binding vested are rarely of high purity, and we generally do not proteins and associated complexes, all tested by chroma- know the degree of change in cytosine methylation in the tin immunoprecipitation (ChIP) approaches coupled disease of interest and thus the quantitative discrimina- with microarray hybridization or massively parallel tion required for an assay, with some studies to date indi- sequencing (MPS). MicroRNAs can be identified and cating that the changes may be quite subtle [7]. These quantified by using microarrays and MPS, while cytosine concerns emphasize the need for cytosine methylation methylation can be definitively studied by converting the assays that can detect methylation levels intermediate in DNA of the genome using sodium bisulfite, shotgun value and changes in disease that are relatively modest in sequencing the product using MPS and mapping this magnitude. Certain microarray-based assays to study back to the genome to count how frequently cytosines cytosine methylation have addressed this issue, with the remain unconverted, indicating their methylation in the methylated DNA immunoprecipitation (meDIP) assay amenable to such quantification when used for CpG islands [8] and possibly also for less CG dinucleotide-rich * Correspondence: [email protected] regions [9]. Restriction enzyme-based assays used with Department of Genetics (Computational Genetics), Center for Epigenomics, Albert Einstein College of Medicine, 1301 Morris Park Avenue, Bronx, NY 10461, microarrays have also proven to be reasonably quantita- USA tive, whether based on methylation-sensitive (for exam- © 2010 Suzuki et al.; licensee BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in BioMed Central any medium, provided the original work is properly cited Suzuki et al. Genome Biology 2010, 11:R36 Page 2 of 11 http://genomebiology.com/2010/11/4/R36 ple, the HELP assay [10]) or methylation-dependent (for were also tested, using the H1 (WA01) human embryonic example, MethylMapper [11]) enzymes. A promising new stem (ES) cell line. The outcome is a modified assay that MPS-based assay is reduced representation bisulfite combines the strengths of MSCC and HELP-seq/Methyl- sequencing (RRBS), which is designed to study the CG- seq, and the supporting analytical workflow that maxi- dense regions defined by short MspI fragments, and pro- mizes the quantitative capabilities of the data generated. vides nucleotide resolution, quantitative data [12]. The use of MPS for what were previously microarray- Results and discussion based assays has been associated with improved perfor- Library preparation and sequencing mance [13], as we found when we modified our HELP We generated HELP tagging libraries with HpaII- or (HpaII tiny fragment Enrichment by Ligation-mediated MspI-digested DNA derived from human ES cells using PCR) assay [10] for MPS, creating an assay similar to the experimental approach shown in Figure 1. The assay Methyl-Seq [14]. The strength of the HELP assay involves differs from MSCC [15] by using EcoP15I instead of the comparison of the HpaII with the methylation-insen- MmeI, generating longer flanking sequences (27 as sitive MspI representation, allowing a normalization step opposed to 18 to 19 bp) and the addition of a T7 poly- that makes the assay semi-quantitative [10]. The HELP merase and reverse transcription step to allow the gener- representation approach was improved upon by Ball et al. ation of the library without contaminating single-adapter [15], who developed the Methyl-Sensitive Cut Counting products, while in addition obviating the need for gel (MSCC) assay, which involves digesting DNA with HpaII, extraction. After the library preparation, a single band of ligating an adapter to the cohesive end formed, using a 125 bp in length was generated, as expected. Libraries restriction enzyme site within the adapter to digest at a flanking sequence and thus capturing the sequence immediate adjacent to the HpaII site. By adding a second me CGG CGG C CCGG C MPS-compatible adapter, a library can be generated for MspI/HpaII digestion MPS, allowing the counting of reads at these sites to rep- CGG T7 resent the degree of methylation at the site. The authors First adapter ligation demonstrated the assay to be reasonably quantitative, T7 testing over 1.3 million sites in the human genome, repre- EcoP15I digestion senting not only HpaII sites clustered in CG-dense T7 regions of the genome (approximately 12% of all HpaII sites are located in annotated CpG islands in the human genome [16]) but also the remaining majority of the A genome in which CG dinucleotides are depleted, a Second adapter ligation genomic compartment not tested by RRBS as currently A RNA designed. A focus on the CG-dense minority of the In vitro transcription genome will fail to observe changes such as those at CG- cDNA depleted promoters (such as OCT4 [17]) and CpG island RT reaction shores [18], and within gene bodies where cytosine meth- ylation has been found to be positively correlated with gene transcription [15]. It is likely, therefore, that an assay PCR amplification system that can study both CG-dense and CG-depleted regions will acquire substantially more information about Massively parallel sequencing with Illumina GA epigenomic states than those directed at the CG-dense compartment alone. Figure 1 HELP-tagging assay design and library preparation. The genomic DNA is digested by HpaII or MspI, the former only cutting at In the current study, we tested whether the use of an CCGG sequences where the central CG dinucleotide is unmethylated. MspI control would improve MSCC assay performance, The first Illumina adapter (AE) is ligated to the compatible cohesive as we had found for microarray-based HELP, and whether end created, juxtaposing an EcoP15I site beside the HpaII/MspI diges- we could develop an analytical pipeline for routine use of tion site and allowing EcoP15I to digest within the flanking DNA se- this assay in epigenome-wide association studies. We also quence as shown. An A overhang is created, allowing the ligation of the second Illumina adapter (AS, green). This will create not only AE-in- explored the use of longer tags than those employed in sert-AS products but also AS-insert-AS molecules. By performing a T7 the MSCC, and added T7 RNA polymerase and reverse polymerase-mediated in vitro transcription from a promoter sequence transcription steps to allow the generation of libraries located on the AE adapter, we can selectively enrich for the AE-insert- without contaminating products, thus obviating the need AS product, following which limited PCR amplification is performed to for gel extraction. The influence of base composition and generate a single sized product for Illumina sequencing. RT, reverse transcription. fragment length parameters as potential sources of bias T7 T7 CGG C T7 CGG C CGG C Suzuki et al. Genome Biology 2010, 11:R36 Page 3 of 11 http://genomebiology.com/2010/11/4/R36 (a) 0 1020304050 MspI count (b) Less methylated More methylated (HpaII count>MspI count) (HpaII count<MspI count) (b,a) (b,a) b BC Angle B=degree(arctan2(b,a)) a a MspI count MspI count Length c= (c) (d) HpaII counts Angle 6.0 R = 0.502 R = 0.826 4.0 2.0 0 20 40 60 80 100 0 20 40 60 80 100 % methylation % methylation Figure 2 Data transformation and bisulfite validation. (a) Scatter plot showing the relationship between the number of HpaII and MspI reads at each locus. (b) The location of the data point on the scatter plot indicates whether it is likely to be less or more methylated with larger or smaller angles B subtended as shown, while the confidence of the measurement will be greater when more reads represent the data point, represented by the length of line c. (c) The HpaII count correlates negatively with the degree of methylation, with more counts occurring at loci with less methylation. (d) Transformation of the data to the B angle measure to normalize HpaII by MspI counts substantially improves the correlation with bisulfite MassAr- ray validation data. HpaII counts HpaII count HpaII count 0 1020304050 HpaII count Angle (degrees) Suzuki et al. Genome Biology 2010, 11:R36 Page 4 of 11 http://genomebiology.com/2010/11/4/R36 were sequenced using an Illumina Genome Analyzer (36 major groups of values in the plot (Figure 2a), separated bp single end reads) and the sequences were analyzed and into loci with high or with minimal HpaII counts. This aligned using Illumina pipeline software version 1.3 or plot helped us to develop a new method for normalizing 1.4. A summary of the Illumina analysis results for each HpaII counts in terms of variability of the MspI represen- replicate is shown in Table S1 in Additional file 1. tation. We recognize that hypomethylated loci are associ- ated with relatively greater HpaII counts and a larger Data quality and reproducibility angle B (Figure 2b, left) whereas methylated loci will be Based on our experimental design, successfully generated defined by smaller angle values (Figure 2b, middle). Fur- products would be expected to possess a 5'-CGCTGCTG thermore, some loci will tend to be sequenced more read- sequence at the 3' end of the read, the first two nucle- ily than others, and may have identical B values but otides (CG) representing the cohesive end for ligation of differing distances from the origin (c distance), allowing a HpaII/MspI digestion products, the remaining six nucle- confidence score for identical methylation values (B) in otides the EcoP15I restriction enzyme recognition site. In terms of the c distance values (Figure 2b, right). To test order to evaluate the yield of desired products, we this model, we used bisulfite MassArray to test quantita- counted the number of reads containing this sequence tively the cytosine methylation values for 61 HpaII sites and plotted the starting positions of this sequence within (Tables S4, S5 and S6 in Additional file 1), choosing loci the reads obtained. We observed that approximately two- representing all components of the B angle spectrum of thirds of the reads contained the expected sequence, and values. In Figure 2c, d we show the correlations between found that the majority was located at base positions 25 these gold standard cytosine methylation values and raw and 26, consistent with the known digestion properties of HpaII counts or B angle values. We find that there is the the restriction enzyme [19]. Removal of the approxi- same negative correlation (R = 0.502) between HpaII mately 30% of reads lacking the CG-EcoP15I sequence counts and cytosine methylation values as demonstrated was performed to eliminate spurious sequences. In order in the MSCC technique [15], and that the angular trans- to investigate sequence quality further, we also deter- formation of the data incorporating the MspI normaliza- mined the number and relative position of Ns (ambigu- tion substantially improves this correlation (R = 0.826), ous base calls) within the reads obtained. Overall, few defining the optimal approach for processing of these reads were found to contain Ns, and where they were data. We represent the data for University of California present, they were found to be evenly distributed by posi- Santa Cruz (UCSC) genome browser visualization as wig- tion within the sequence. To test data reproducibility, we gle tracks, with higher B angle values defining less methy- compared the results of three experimental replicates lated loci. Methylated loci with zero values that would be against each other using the Pearson correlation coeffi- otherwise difficult to visualize as having been tested are cient metric. The results of this study showed that all rep- represented as small negative values. We show the details licates were highly correlated (all the r values exceed 0.9), of the analytical workflow in Figure 3 and an example of a which confirmed that the technical reproducibility of this UCSC genome browser representation of HELP-tagging assay was excellent (Table S2 in Additional file 1). data in Figure S1 in Additional file 1. All data are available through the Gene Expression Omnibus database (acces- Distribution of MspI/HpaII sequence tags sion number [GEO:GSE19937]) and as UCSC genome We merged three lanes of MspI data and observed that browser tracks [20]. approximately 80% of the 2,292,198 annotated HpaII sites in the human genome (hg18) were represented by at least Potential sources of bias: base composition and fragment one read, for a total of over 1.8 million loci throughout length the genome. The mean numbers of reads per locus for As the number of reads at CCGG sites following MspI MspI and HpaII were 3.94 and 1.82, respectively, and digestion should not be influenced by methylation, the MspI counts were distributed evenly across all genomic representation obtained from MspI digestion allowed us compartments examined (Table S3 in Additional file 1). to look for systematic sources of bias inherent to the We hypothesize that a combination of incomplete assay. A major concern was that base composition could genomic coverage and polymorphisms within some be a source of such bias, as it has been reported that Illu- CCGG sites (as we have previously observed [10]) mina sequencing can be influenced by GC composition accounts for the 20% of HpaII sites that were not repre- [21], possibly because of the gel extraction step [22]. Our sented by any reads. protocol does not require gel extraction and only begins to show an under-representation of sequences when the Normalization of HpaII by MspI counts and data (G+C) content exceeds approximately 80% (Figure 4a). transformation We also tested to see whether the sizes of the MspI frag- When we plot the MspI count on the x-axis and HpaII ments generated influenced the counts obtained, as the count on the y-axis for each HpaII site, we can see two Suzuki et al. Genome Biology 2010, 11:R36 Page 5 of 11 http://genomebiology.com/2010/11/4/R36 ELAND (1-36 bp alignment) Scan for EcoP15I tag No Discard the reads Yes Mask tag sequence as “n” ELAND (2-28 bp alignment) Quality failed/not aligned Discard the reads Unique/multiple aligned Distance from neighboring HpaII site < 27 bp Discard the reads ≥27 bp Multiple Unique aligned hits hit aligned >10 Discard the reads aligned ≤10 Weighted based on aligned number Map to annotated HpaII site? No Putative polymorphic HpaII sites dbSNP check Yes Count hit number MspI=0 Putative polymorphic HpaII sites dbSNP check MspI>0 Normalize HpaII with MspI by angle calculation UCSC genome browser Figure 3 HELP-tagging analysis workflow. The analysis workflow for HELP-tagging data is illustrated. Only sequence reads that contain the adapter sequence and map to a single or ≤ 10 sites are retained, the latter repetitive sequences distributed by weighting among the matched loci. Potential polymorphic loci are annotated. Normalization of HpaII by MspI using the angle calculation described in the previous figure is performed and files are generated for genome browser visualization. UCSC, University of California Santa Cruz. Suzuki et al. Genome Biology 2010, 11:R36 Page 6 of 11 http://genomebiology.com/2010/11/4/R36 Identification of polymorphic CCGG sequences Whereas MSCC used MmeI and generates an 18- to 19- (a) 0.12 bp sequence flanking the HpaII site [15], our use of Actual EcoP15I generates a 27-bp flanking sequence. We asked 0.10 MspI count whether this size difference influenced our ability to align sequences to the reference genome. We truncated our 0.08 sequence reads to 19 bp to mimic the MSCC read length 0.06 and found that this caused a profound loss of ability to align reads unambiguously (Table S7 in Additional file 1). 0.04 To compensate for the low alignment rate, the MSCC report described an ingenious strategy of alignment to 0.02 the sequences immediately flanking the annotated HpaII 0 sites in the reference genome [15], an approach suffi- ciently powerful that it generated the well-validated data 0 20406080 100 that they described. However, it does not offer the possi- GC content (%) bility of identifying polymorphic HpaII sites at the high frequencies that we previously observed for our HELP- (b) seq assay [10]. We tested whether our longer sequences 0.025 allowed the identification of loci at which an HpaII site is Actual annotated in the reference genome but we obtain no MspI count 0.020 sequence reads, and the opposite situation where we observed at least four MspI reads (the average number 0.015 per annotated MspI/HpaII site) flanking a locus not annotated in the reference genome. In Table S8 in Addi- 0.010 tional file 1 we list approximately 6,600 candidate poly- morphic HpaII sites, of which examples are shown in 0.005 Figure 5, confirmed by targeted resequencing of those loci. The 6,600 loci were selected based on overlap with dbSNP entries, allowing us to evaluate the pattern of 0 100 200 300 400 500 sequence variability at these loci. Approximately 80% of the SNPs are C:G to T:A transversions, consistent with Length (bp) deamination-mediated decay of methylcytosine being the Figure 4 Base composition and fragment length influences on se- cause of the polymorphism [23]. Polymorphic CG dinu- quence counts. (a) The proportion of (G+C) nucleotides was calculat- cleotides are major potential sources of error not only for ed for the 50-bp sequence centered around each annotated CCGG in microarrays, which are designed to a consensus genomic the reference human genome. The base composition of all of the MspI sequence, but also for both bisulfite sequencing, which sequences generated from the human ES cell line studied was also cal- culated. The relative proportion for (G+C) content in 2% bins for each would read the C to T transversion as unmethylated, and set of data was calculated and plotted as shown. The black line shows mass spectrometry-based assays, requiring the develop- the proportions in the reference genome, while the red line illustrates ment of specific analytical approaches such as we have the distribution we observed in our MspI experiment. Two peaks rep- described [24]. resenting base composition in repetitive sequences are apparent. The MspI distribution closely matches the expected distribution except DNA methylation studies of human embryonic stem cells when the base composition exceeds approximately 80%, when it is slightly under-represented. (b) We calculated the relative frequencies To test whether the HELP-tagging assay was generating of MspI digestion product sizes in the human reference genome. In data that are biologically plausible, we tested the methyla- this case we found that the shorter fragments are more likely to be se- tion of different genomic sequence compartments as den- quenced than larger (≥300 bp) fragments. The three major peaks ob- sity plots of B angle values for the human ES cells used in served represent Alu short interspersed repetitive element (SINE) these studies. In Figure S2a in Additional file 1 we show sequences. how promoters (defined as -2 kb to 2 kb from the tran- scription start site of RefSeq genes), gene bodies (the digestion by type III endonucleases like EcoP15I is most remaining region within the RefSeq gene) and intergenic efficient when a pair of enzymes is present in convergent (all other) sequences compare, finding the expected orientation on the same DNA molecule [19]. We find that enrichment of hypomethylated loci with larger B angle there is indeed an over-representation for shorter (≤300 values in promoter regions. When we compared unique bp) and a corresponding modest under-representation with repetitive sequences, again we found the expected for larger MspI fragments (Figure 4b). Counts relative to total Counts relative to total Suzuki et al. Genome Biology 2010, 11:R36 Page 7 of 11 http://genomebiology.com/2010/11/4/R36 (a) chr15: 25725400 25725500 25725600 25725700 25725800 25725900 MspI hit HpaII OCA2 CpG island dbSNP chr15: 25725640 25725650 25725660 25725670 25725680 A G T C T C T T CA C T C T CA CA T T C T A G C C C G G G G CT CC T G CC C A C A T T C T G C A T G G C A T G G C C T Reference HpaII OCA2 dbSNP rs12916836 rs12905726 A G TC T C T TCA C T C T CA CA T T C T A G C C C A A A GG C T CC T G CCC A C A C T C T G C A T GG C A T G G C C T Observed Trace data (b) 3870000 3875000 3880000 3885000 3890000 3895000 chr2: MspI hit HpaII dbSNP 3882300 3882310 3882320 3882330 3882340 3882350 chr2: G T C T G G A G C A G A G G C T T C T A A G CA CA GC A T C TT G G C C A A C G A A GC CA G C AC CAC A G G C A G G C A C T Reference MspI hit dbSNP rs6748872 G T C T G G A G CA G A G G C T TC T A A G CA CA G CA T G G C CA A C G A A G C C A G CA C CA CA G G C A G G C A C T Observed Trace data Figure 5 Polymorphic HpaII sites identified by HELP-tagging. Examples of HpaII sites (a) annotated in the reference genome sequence but not represented by MspI reads or (b) not annotated in the reference human genome and represented by at least four MspI reads are shown. In each case there is a SNP defined by dbSNP that indicates the C:T to G:A transversion that eliminates or restores the CCGG HpaII site. pattern of increased methylation of repetitive DNA com- that while transposable elements are generally methy- pared with unique sequences (Figure S2b in Additional lated and are depleted near gene promoters, those that file 1). Combining these observations, we tested whether are proximal to promoters tend to be less methylated the transposable element component of annotated repeti- than those located more distally. While many types of tive DNA sequences showed any tendency to unusual transposable elements were represented in this pro- methylation near gene promoters. In Figure 6 we show moter-proximal hypomethylated group, we found a sub- Suzuki et al. Genome Biology 2010, 11:R36 Page 8 of 11 http://genomebiology.com/2010/11/4/R36 (a) (b) 100% Angle 0-30 31-60 80% 61-90 60% 400 40% 200 20% 0 0% Length from TSS (bp) Length from TSS (bp) Figure 6 Identification of a position effect on DNA methylation in transposable elements located close to gene promoters. The distance from RefSeq gene transcription start sites and DNA methylation status are shown. The x-axis displays the distance from transcription start sites (TSSs). HpaII sites were categorized into three groups by angle, 0 to 30 (blue), 31 to 60 (red) and 61 to 90 (green)). (a) Number of HpaII sites; (b) proportions of each angle category (%). set to be the most markedly over-represented, as shown neous sites is that they contribute to the proportion of in Figure 6c. loci at which we could not obtain sequence reads. The outcome of these studies was an improvement in Our exploration of the distribution of cytosine methy- the previously described MSCC [15] and HELP-seq [10] lation in the same human ES cell line studied by Lister et assays, not only by means of technical modifications such al. [6] showed consistent results, with hypomethylation as the use of EcoP15I but also because of the concurrent of transcription start sites and methylation of transpos- use of MspI for normalization. The effect of these modifi- able elements, as expected from long-standing observa- cations was not only to increase the accuracy of the assay tions in the field. We furthermore discovered a limited but also to enhance the ability to align sequences to the subset of transposable elements that is hypomethylated genome and thus identify polymorphic HpaII/MspI sites. when in close proximity to transcription start sites. When The means of normalization of HpaII by MspI using an this subset was studied to determine whether certain angular metric is an innovation that improved the data types of transposable elements were disproportionately accuracy substantially and may have applications in other over-represented, we found two broad classes, one of MPS assay normalization strategies. We were also able to transposable element fossils with no innate capacity to discard reads that did not contain the expected adapter replicate themselves (the ancient DNA, long interspersed sequences, and created a straightforward data analytical repetitive elements (LINEs) and short interspersed repet- pipeline that will facilitate processing of these HELP-tag- itive elements (SINEs) shown in Figure 6c) and younger ging data by others. ERV1 long terminal repeat retroelements. Loss of methy- The potential sources of systematic artifacts due to base lation of functionally inactive transposable elements is composition or digestion product size were evaluated. likely to be of no negative consequence to the host Apart from a modest decrease in representation in genome, consistent with the host defense hypothesis [25], regions above approximately 80% (G+C) content, base while the young ERV long terminal repeats represent a composition did not cause biases in representations, pos- group of transposons whose function has been harnessed sibly in part due to our avoidance of a gel purification step as promoters of endogeneous genes [26,27]. This obser- in library preparation [22]. Fragment length does influ- vation demonstrates the value of a high-resolution, ence the outcome, most likely due to effects on EcoP15I genome-wide assay like HELP-tagging to define potential digestion [19], although the effects should be similar for functional elements in an unbiased manner. both HpaII and MspI and should, therefore, largely cancel each other out in the normalization step. It is possible Conclusions that endogeneous EcoP15I sites could influence the rep- We propose that MPS-based assays such as RRBS [12], resentations, but to have an effect they would have to be MSCC [15] and HELP-tagging will prove to be the assays located within the 27 bp adjacent to HpaII/MspI sites and of choice for epigenome-wide association studies in would cause digestion of the ligated adapter, causing human disease, with the latter two preferable as we begin those loci to be under-represented in both HpaII and to explore the CG-depleted majority of the genome. It MspI datasets. The most likely effect of these endoge- should not be necessary to run MspI assays every time a HELP-tagging assay is performed, suggesting that a com- -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 HpaII counts Percentage Suzuki et al. Genome Biology 2010, 11:R36 Page 9 of 11 http://genomebiology.com/2010/11/4/R36 mon MspI dataset can serve as a universal reference for a ments using a New England Biolabs Quick Ligation Kit species, allowing a single lane of Illumina sequencing of (25 μl of 2× Quick ligase buffer, 2.5 μl of adapter AS (10 the HpaII library to provide the methylation data for that μM), 2.5 μl of Quick Ligase in a final volume 50 μl). After sample. The development of analytical pipelines to sup- ligation, products were purified using the MinElute PCR port analysis of these datasets will be critical to the suc- purification kit (Qiagen, Hilden, Germany) and in vitro- cess of these projects, while the careful ongoing transcribed using the Ambion MEGAshortscriptkit (Life assessment of potential sources of bias will also be essen- Technologies, Carlsbad, CA, USA). Following in vitro tial for improving assay performance. transcription, products were purified with an RNeasy clean-up kit (Qiagen) before reverse transcription was performed using the Invitrogen SuperScript III kit (Life Materials and methods Cell preparation and DNA purification Technologies). The first strand cDNA produced was used H1 human ES cells (NIH code WA01 from Wicell as a template for PCR using the following conditions: Research Institute, Madison, WI, USA) were cultured on 96°C for 2 minutes, then 18 cycles of 96°C for 15 seconds matrigel (BD Biosciences, San Diego, CA, USA), at 37°C, and 72°C for 15 seconds followed by 5 minutes at 72°C for 5% O and 5% CO . Amplified human ES cell pluripo- the final extension. After PCR, the library was purified 2 2 using a QIAQuick PCR clean-up kit (Qiagen). tency was assessed by flow cytometry with SSEA4, CD24 and Oct4 markers. To extract DNA, the cells were sus- Single-locus quantitative validation assays pended in 10 ml of a solution of 10 mM Tris-HCl (pH Bisulfite conversion and MassArray (Sequenom, San 8.0), 0.1 M EDTA and 1 ml of 10% SDS to which 10 μl of Diego, CA, USA) were performed using an aliquot of the RNase A (20 mg/ml) was added. After incubation for 1 same sample of DNA as was used for the high-through- hour at 37°C, 50 μl of proteinase K (20 mg/ml) was added put assays described above. Bisulfite conversion was per- and the solution was gently mixed and incubated in a formed with an EZ DNA Methylation kit (Zymo 50°C water bath overnight. To purify the lysate, it was Research, Orange, CA, USA). Bisulfite primers were extracted three times using saturated phenol, then twice designed using MethPrimer [28], specifying the desired with chloroform, and dialyzed for 16 hours at 4°C against product length (250 to 450 bp), primer length (23 to 29 three changes of 0.2× SSC. Following dialysis, the DNA bp) and primer Tm (56 to 62°C). PCR was performed was concentrated by coating the dialysis bags in polyeth- using FastStart High Fidelity Taq polymerase (Roche, ylene glycol (molecular weight 20,000). The purity and Basel, Switzerland) with the following conditions: 95°C final concentration of the purified DNA was checked by for 10 minutes, then 42 cycles of 95°C for 30 seconds, spectrometry (Nanodrop, Wilmington, DE, USA). primer-specific Tm for 30 seconds and 72°C for 1 minute, followed by 72°C for 10 minutes for the final extension. Illumina library preparation Primer-specific Tm and sequence information are pro- The sample preparation steps are illustrated in Figure 1. vided in Table S6 in Additional file 1. Bisulfite MassArray Two custom adapters were created for HELP-tagging, assays were performed by the institutional Genomics referred to as AE and AS. As well as an Illumina adapter Core Facility. The data were analyzed using the analytical sequence, adapter AE contains an EcoP15I recognition pipeline we have previously described [24]. site and a T7 promoter sequence. Adapter AS contains an Illumina sequencing primer sequence. The adapter and Bioinformatic analysis primer sequences for library preparation are listed in Four lanes of sequencing were performed using an Illu- Table S9 in Additional file 1. Genomic DNA (5 μg) was mina GA IIx Sequencer at the institutional Epigenomics digested with HpaII and MspI in separate 200 μl reactions Shared Facility. Three lanes were used for technical repli- and purified by phenol/chloroform extraction followed cates of MspI, for the methylation-insensitive reference by ethanol precipitation. The digested genomic DNA was dataset. Images generated by the Illumina sequencer were ligated to adapter AE using a New England Biolabs Quick analyzed by Illumina pipeline software (versions 1.3 to Ligation Kit (25 μl of 2× Quick ligase buffer, 3 μg of 1.4). Initial data processing was performed using the HpaII-digested DNA or 1 μg of MspI-digested DNA, 0.1 default read length of 36 bp, after which we isolated the μl of Adapter AE (1 μM), 3 μl of Quick Ligase in a final sequences in which we found adapter sequences on the volume of 50 μl). After AE ligation, the products were 3'-end, replaced the adapter sequence with a poly(N) purified using Agencourt AMpure beads (Beckman sequence of the same length, and re-ran the Illumina Coulter, Brea CA, USA), then digested with EcoP15I ELAND pipeline again on these sequences with the (New England Biolabs). The restriction fragments were sequence length set at 27 bp (the 2 to 28 bp subsequence). end-repaired to inhibit to dimerization of adapters, and The data within the ELAND_extended.txt files were used tailed with a single dA, at the 3' end. After the dA tailing for counting the number of aligned sequences adjacent to reaction, adapter AS was ligated to the dA-tailed frag- Suzuki et al. Genome Biology 2010, 11:R36 Page 10 of 11 http://genomebiology.com/2010/11/4/R36 reduced representation bisulfite sequencing; UCSC: University of California each CCGG (HpaII/MspI) site annotated in the hg18 Santa Cruz. freeze of the human genome at the UCSC genome browser. We permitted up to two mismatches in each Authors' contributions MS and JMG designed the assays and strategies for its analysis, MS performed sequence, and allowed a sequence to align to up to a max- all library preparation and characterization, MS, DL and MP performed bisulfite imum of 10 locations within the genome. For non-unique validation studies, while QJ and AMcL performed computational analyses. JMG alignments, a sequence was assigned a partial count for and MS prepared the manuscript. each alignment location amounting to 1/n, where n rep- Acknowledgements resents the total number of aligned positions. To normal- This work is supported by a grant from the National Institute of Health (NIH, ize the data between experiments, the number of R01 HG004401) to JMG. The authors thank Shahina Maqbool PhD, Raul Olea sequences associated with each HpaII site was divided by and Gael Westby of the Einstein Epigenomics Shared Facility for their contribu- tions, Drs Eric Bouhassira and Emmanuel Olivier (Einstein) for the WA01/H1 the total number of sequences (including partial counts) human ES cell line, and Einstein's Center for Epigenomics. aligning to all HpaII sites in the same sample. We refer to this figure as the fixed count below. Author Details Department of Genetics (Computational Genetics), Center for Epigenomics, To examine an influence of (G+C) mononucleotide Albert Einstein College of Medicine, 1301 Morris Park Avenue, Bronx, NY 10461, content on counts of sequences obtained, we extracted USA the (G+C) annotation from the hg18 freeze of the human Received: 8 January 2010 Revised: 16 March 2010 genome at UCSC and examined the distribution of Accepted: 1 April 2010 Published: 1 April 2010 T © T Ge h h 2010 S i i no s s arti ime s an cle B u o iz o i p u ls o e k ava g n i e y acce 2010, t al. ilable ;ss arti lice f 11 rn o:R36 s m cle ee : h Bi dttp:/ istri oMbu / eg de te C nd e on u m tn ral L ed be io r th tlo d.g e y t .co erm ms o /2010/ f the 11/ Cre 4/ aR ti3 ve 6 Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited sequence counts according to (G+C) content. Annotated percentages of (G+C) content were available for adjacent References 1. Esteller M: Epigenetics in cancer. N Engl J Med 2008, 358:1148-1159. 5-bp windows. For each annotated HpaII site, we calcu- 2. Kopelovich L, Crowell JA, Fay JR: The epigenome as a target for cancer lated the mean percentage (G+C) for a 50-bp region cen- chemoprevention. J Natl Cancer Inst 2003, 95:1747-1757. tered at the restriction site. Counts of sequences 3. Zilberman D, Gehring M, Tran RK, Ballinger T, Henikoff S: Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an associated with HpaII sites were obtained for 50 sequen- interdependence between methylation and transcription. Nat Genet tial non-overlapping windows of 2% (G+C) (the mini- 2007, 39:61-69. mum possible in a sample of 50-bp regions). These data 4. Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE: Shotgun bisulphite were then normalized as a proportion of the total number sequencing of the Arabidopsis genome reveals DNA methylation of fragments. Comparisons were made to the expected patterning. Nature 2008, 452:215-219. frequencies, which for each 2% (G+C) bin was repre- 5. Lister R, O'Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, Ecker JR: Highly integrated single-base resolution maps of the sented by the counts of HpaII sites falling within a range epigenome in Arabidopsis. Cell 2008, 133:523-536. relative to the total number of HpaII sites in the genome. 6. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery This analysis was performed on both HpaII and MspI- JR, Lee L, Ye Z, Ngo QM, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar AH, Thomson JA, Ren B, Ecker JR: Human DNA methylomes at digested DNA for comparison. base resolution show widespread epigenomic differences. Nature The potential effect of distance between HpaII sites on 2009, 462:315-322. sequences counts obtained at each HpaII site was mea- 7. Heijmans BT, Tobi EW, Stein AD, Putter H, Blauw GJ, Susser ES, Slagboom PE, Lumey LH: Persistent epigenetic differences associated with sured by summing the counts of sequences aligning prenatal exposure to famine in humans. Proc Natl Acad Sci USA 2008, within each restriction fragment, and normalizing the 105:17046-17049. result with respect to total sequence count. As with the 8. Yamashita S, Hosoya K, Gyobu K, Takeshima H, Ushijima T: Development of a novel output value for quantitative assessment in methylated DNA (G+C) analysis above, this was performed for both MspI immunoprecipitation-CpG island microarray analysis. DNA Res 2009, and HpaII digested restriction fragments. The data were 16:275-286. compared with the expected distribution determined by 9. Down TA, Rakyan VK, Turner DJ, Flicek P, Li H, Kulesha E, Graf S, Johnson N, Herrero J, Tomazou EM, Thorne NP, Backdahl L, Herberth M, Howe KL, performing virtual restriction digestion using genomic Jackson DK, Miretti MM, Marioni JC, Birney E, Hubbard TJ, Durbin R, Tavare HpaII site coordinates, and normalizing the number of S, Beck S: A Bayesian deconvolution strategy for immunoprecipitation- virtual fragments of each size with respect to the total based DNA methylome analysis. Nat Biotechnol 2008, 26:779-785. 10. Oda M, Glass JL, Thompson RF, Mo Y, Olivier EN, Figueroa ME, Selzer RR, number of these virtual fragments. Richmond TA, Zhang X, Dannenberg L, Green RD, Melnick A, Hatchwell E, Bouhassira EE, Verma A, Suzuki M, Greally JM: High-resolution genome- Additional material wide cytosine methylation profiling with simultaneous copy number analysis and optimization for limited cell numbers. Nucleic Acids Res 2009, 37:3829-3839. Additional file 1 Supplemental data containing two figures (Figures S1 11. Ordway JM, Bedell JA, Citek RW, Nunberg AN, Jeddeloh JA: and S2) and nine tables ( Tables S1 to S9). MethylMapper: a method for high-throughput, multilocus bisulfite sequence analysis and reporting. Biotechniques 2005, 39:. 464, 466, 468 Abbreviations passim. bp: base pair; CG/CpG: cytosine-guanine dinucleotide; ChIP: chromatin immu- 12. Meissner A, Mikkelsen TS, Gu H, Wernig M, Hanna J, Sivachenko A, Zhang noprecipitation; ES: embryonic stem; (G+C): guanine and cytosine mononucle- X, Bernstein BE, Nusbaum C, Jaffe DB, Gnirke A, Jaenisch R, Lander ES: otides; HELP: HpaII tiny fragment Enrichment by Ligation-mediated PCR; MPS: Genome-scale DNA methylation maps of pluripotent and massively-parallel sequencing; MSCC: methyl-sensitive cut counting; RRBS: differentiated cells. Nature 2008, 454:766-770. Suzuki et al. Genome Biology 2010, 11:R36 Page 11 of 11 http://genomebiology.com/2010/11/4/R36 13. t Hoen PA, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RH, de Menezes RX, Boer JM, van Ommen GJ, den Dunnen JT: Deep sequencing- based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res 2008, 36:e141. 14. Brunner AL, Johnson DS, Kim SW, Valouev A, Reddy TE, Neff NF, Anton E, Medina C, Nguyen L, Chiao E, Oyolu CB, Schroth GP, Absher DM, Baker JC, Myers RM: Distinct DNA methylation patterns characterize differentiated human embryonic stem cells and developing human fetal liver. Genome Res 2009, 19:1044-1056. 15. Ball MP, Li JB, Gao Y, Lee JH, LeProust EM, Park IH, Xie B, Daley GQ, Church GM: Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells. Nat Biotechnol 2009, 27:361-368. 16. Fazzari MJ, Greally JM: Epigenomics: beyond CpG islands. Nat Rev Genet 2004, 5:446-455. 17. Hattori N, Nishino K, Ko YG, Hattori N, Ohgane J, Tanaka S, Shiota K: Epigenetic control of mouse Oct-4 gene expression in embryonic stem cells and trophoblast stem cells. J Biol Chem 2004, 279:17063-17069. 18. Irizarry RA, Ladd-Acosta C, Wen B, Wu Z, Montano C, Onyango P, Cui H, Gabo K, Rongione M, Webster M, Ji H, Potash JB, Sabunciyan S, Feinberg AP: The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific CpG island shores. Nat Genet 2009, 41:178-186. 19. Moncke-Buchner E, Rothenberg M, Reich S, Wagenfuhr K, Matsumura H, Terauchi R, Kruger DH, Reuter M: Functional characterization and modulation of the DNA cleavage efficiency of type III restriction endonuclease EcoP15I in its interaction with two sites in the DNA target. J Mol Biol 2009, 387:1309-1319. 20. Greallylab [http://greallylab.aecom.yu.edu/HELP-tagging/] 21. Dohm JC, Lottaz C, Borodina T, Himmelbauer H: Substantial biases in ultra-short read data sets from high-throughput DNA sequencing. Nucleic Acids Res 2008, 36:e105. 22. Quail MA, Kozarewa I, Smith F, Scally A, Stephens PJ, Durbin R, Swerdlow H, Turner DJ: A large genome center's improvements to the Illumina sequencing system. Nat Methods 2008, 5:1005-1010. 23. Duncan BK, Miller JH: Mutagenic deamination of cytosine residues in DNA. Nature 1980, 287:560-561. 24. Thompson RF, Suzuki M, Lau KW, Greally JM: A pipeline for the quantitative analysis of CG dinucleotide methylation using mass spectrometry. Bioinformatics 2009, 25:2164-2170. 25. Yoder JA, Walsh CP, Bestor TH: Cytosine methylation and the ecology of intragenomic parasites. Trends Genet 1997, 13:335-340. 26. Conley AB, Piriyapongsa J, Jordan IK: Retroviral promoters in the human genome. Bioinformatics 2008, 24:1563-1567. 27. Cohen CJ, Lock WM, Mager DL: Endogenous retroviral LTRs as promoters for human genes: a critical assessment. Gene 2009, 448:105-114. 28. MethPrimer [http://www.urogene.org/methprimer/] doi: 10.1186/gb-2010-11-4-r36 Cite this article as: Suzuki et al., Optimized design and data analysis of tag- based cytosine methylation assays Genome Biology 2010, 11:R36 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology Springer Journals

Optimized design and data analysis of tag-based cytosine methylation assays

Download PDF
11 pages

Loading next page...
 
/lp/springer-journals/optimized-design-and-data-analysis-of-tag-based-cytosine-methylation-j5XkccJkNw

References (60)

Publisher
Springer Journals
Copyright
Copyright © 2010 by Suzuki et al.; licensee BioMed Central Ltd.
Subject
Life Sciences; Animal Genetics and Genomics; Human Genetics; Plant Genetics & Genomics; Microbial Genetics and Genomics; Bioinformatics; Evolutionary Biology
eISSN
1474-760X
DOI
10.1186/gb-2010-11-4-r36
pmid
20359321
Publisher site
See Article on Publisher Site

Abstract

Using the type III restriction-modification enzyme EcoP15I, we isolated sequences flanking sites digested by the methylation-sensitive HpaII enzyme or its methylation-insensitive MspI isoschizomer for massively parallel sequencing. A novel data transformation allows us to normalise HpaII by MspI counts, resulting in more accurate quantification of methylation at >1.8 million loci in the human genome. This HELP-tagging assay is not sensitive to sequence polymorphism or base composition and allows exploration of both CG-rich and depleted genomic contexts. Background starting material, due to the resistance of methylcytosine Epigenetic mechanisms of transcriptional regulation are to bisulfite conversion compared with unmethylated increasingly being studied for their potential influences in cytosines. This allows nucleotide resolution, strand-spe- human disease pathogenesis. Much of this interest is cific, quantitative assessment of cytosine methylation, based on the paradigm of neoplastic transformation, in with such studies performed in Arabidopsis [3-5] and which epigenetic changes appear to be universal, wide- human cells to date [6]. spread throughout the genome, causative of critical tran- While this approach represents the ideal means of test- scriptional changes and predictive of disease prognosis ing cytosine methylation, the amount of sequencing nec- (reviewed in [1]). Furthermore, these epigenetic changes essary (for the human genome, over 1 billion sequences of represent potential pharmacological targets for reversal ~75 bp each [6]) to generate quantitative information and amelioration of the disease process [2]. genome-wide remains prohibitive in terms of cost, limit- Of the large number of regulatory processes referred to ing these studies to the few referred to above. When as epigenetic, there exist numerous assays to study chro- studying human disease, the emphasis remains on cyto- matin component distribution, cytosine methylation and sine methylation assays, as it is generally easier to collect microRNA expression genome-wide. The chromatin clinical samples for DNA purification than for ChIP or components include a large number of post-translational even RNA assays. However, the cell populations har- modifications of histones, variant histones, DNA-binding vested are rarely of high purity, and we generally do not proteins and associated complexes, all tested by chroma- know the degree of change in cytosine methylation in the tin immunoprecipitation (ChIP) approaches coupled disease of interest and thus the quantitative discrimina- with microarray hybridization or massively parallel tion required for an assay, with some studies to date indi- sequencing (MPS). MicroRNAs can be identified and cating that the changes may be quite subtle [7]. These quantified by using microarrays and MPS, while cytosine concerns emphasize the need for cytosine methylation methylation can be definitively studied by converting the assays that can detect methylation levels intermediate in DNA of the genome using sodium bisulfite, shotgun value and changes in disease that are relatively modest in sequencing the product using MPS and mapping this magnitude. Certain microarray-based assays to study back to the genome to count how frequently cytosines cytosine methylation have addressed this issue, with the remain unconverted, indicating their methylation in the methylated DNA immunoprecipitation (meDIP) assay amenable to such quantification when used for CpG islands [8] and possibly also for less CG dinucleotide-rich * Correspondence: [email protected] regions [9]. Restriction enzyme-based assays used with Department of Genetics (Computational Genetics), Center for Epigenomics, Albert Einstein College of Medicine, 1301 Morris Park Avenue, Bronx, NY 10461, microarrays have also proven to be reasonably quantita- USA tive, whether based on methylation-sensitive (for exam- © 2010 Suzuki et al.; licensee BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in BioMed Central any medium, provided the original work is properly cited Suzuki et al. Genome Biology 2010, 11:R36 Page 2 of 11 http://genomebiology.com/2010/11/4/R36 ple, the HELP assay [10]) or methylation-dependent (for were also tested, using the H1 (WA01) human embryonic example, MethylMapper [11]) enzymes. A promising new stem (ES) cell line. The outcome is a modified assay that MPS-based assay is reduced representation bisulfite combines the strengths of MSCC and HELP-seq/Methyl- sequencing (RRBS), which is designed to study the CG- seq, and the supporting analytical workflow that maxi- dense regions defined by short MspI fragments, and pro- mizes the quantitative capabilities of the data generated. vides nucleotide resolution, quantitative data [12]. The use of MPS for what were previously microarray- Results and discussion based assays has been associated with improved perfor- Library preparation and sequencing mance [13], as we found when we modified our HELP We generated HELP tagging libraries with HpaII- or (HpaII tiny fragment Enrichment by Ligation-mediated MspI-digested DNA derived from human ES cells using PCR) assay [10] for MPS, creating an assay similar to the experimental approach shown in Figure 1. The assay Methyl-Seq [14]. The strength of the HELP assay involves differs from MSCC [15] by using EcoP15I instead of the comparison of the HpaII with the methylation-insen- MmeI, generating longer flanking sequences (27 as sitive MspI representation, allowing a normalization step opposed to 18 to 19 bp) and the addition of a T7 poly- that makes the assay semi-quantitative [10]. The HELP merase and reverse transcription step to allow the gener- representation approach was improved upon by Ball et al. ation of the library without contaminating single-adapter [15], who developed the Methyl-Sensitive Cut Counting products, while in addition obviating the need for gel (MSCC) assay, which involves digesting DNA with HpaII, extraction. After the library preparation, a single band of ligating an adapter to the cohesive end formed, using a 125 bp in length was generated, as expected. Libraries restriction enzyme site within the adapter to digest at a flanking sequence and thus capturing the sequence immediate adjacent to the HpaII site. By adding a second me CGG CGG C CCGG C MPS-compatible adapter, a library can be generated for MspI/HpaII digestion MPS, allowing the counting of reads at these sites to rep- CGG T7 resent the degree of methylation at the site. The authors First adapter ligation demonstrated the assay to be reasonably quantitative, T7 testing over 1.3 million sites in the human genome, repre- EcoP15I digestion senting not only HpaII sites clustered in CG-dense T7 regions of the genome (approximately 12% of all HpaII sites are located in annotated CpG islands in the human genome [16]) but also the remaining majority of the A genome in which CG dinucleotides are depleted, a Second adapter ligation genomic compartment not tested by RRBS as currently A RNA designed. A focus on the CG-dense minority of the In vitro transcription genome will fail to observe changes such as those at CG- cDNA depleted promoters (such as OCT4 [17]) and CpG island RT reaction shores [18], and within gene bodies where cytosine meth- ylation has been found to be positively correlated with gene transcription [15]. It is likely, therefore, that an assay PCR amplification system that can study both CG-dense and CG-depleted regions will acquire substantially more information about Massively parallel sequencing with Illumina GA epigenomic states than those directed at the CG-dense compartment alone. Figure 1 HELP-tagging assay design and library preparation. The genomic DNA is digested by HpaII or MspI, the former only cutting at In the current study, we tested whether the use of an CCGG sequences where the central CG dinucleotide is unmethylated. MspI control would improve MSCC assay performance, The first Illumina adapter (AE) is ligated to the compatible cohesive as we had found for microarray-based HELP, and whether end created, juxtaposing an EcoP15I site beside the HpaII/MspI diges- we could develop an analytical pipeline for routine use of tion site and allowing EcoP15I to digest within the flanking DNA se- this assay in epigenome-wide association studies. We also quence as shown. An A overhang is created, allowing the ligation of the second Illumina adapter (AS, green). This will create not only AE-in- explored the use of longer tags than those employed in sert-AS products but also AS-insert-AS molecules. By performing a T7 the MSCC, and added T7 RNA polymerase and reverse polymerase-mediated in vitro transcription from a promoter sequence transcription steps to allow the generation of libraries located on the AE adapter, we can selectively enrich for the AE-insert- without contaminating products, thus obviating the need AS product, following which limited PCR amplification is performed to for gel extraction. The influence of base composition and generate a single sized product for Illumina sequencing. RT, reverse transcription. fragment length parameters as potential sources of bias T7 T7 CGG C T7 CGG C CGG C Suzuki et al. Genome Biology 2010, 11:R36 Page 3 of 11 http://genomebiology.com/2010/11/4/R36 (a) 0 1020304050 MspI count (b) Less methylated More methylated (HpaII count>MspI count) (HpaII count<MspI count) (b,a) (b,a) b BC Angle B=degree(arctan2(b,a)) a a MspI count MspI count Length c= (c) (d) HpaII counts Angle 6.0 R = 0.502 R = 0.826 4.0 2.0 0 20 40 60 80 100 0 20 40 60 80 100 % methylation % methylation Figure 2 Data transformation and bisulfite validation. (a) Scatter plot showing the relationship between the number of HpaII and MspI reads at each locus. (b) The location of the data point on the scatter plot indicates whether it is likely to be less or more methylated with larger or smaller angles B subtended as shown, while the confidence of the measurement will be greater when more reads represent the data point, represented by the length of line c. (c) The HpaII count correlates negatively with the degree of methylation, with more counts occurring at loci with less methylation. (d) Transformation of the data to the B angle measure to normalize HpaII by MspI counts substantially improves the correlation with bisulfite MassAr- ray validation data. HpaII counts HpaII count HpaII count 0 1020304050 HpaII count Angle (degrees) Suzuki et al. Genome Biology 2010, 11:R36 Page 4 of 11 http://genomebiology.com/2010/11/4/R36 were sequenced using an Illumina Genome Analyzer (36 major groups of values in the plot (Figure 2a), separated bp single end reads) and the sequences were analyzed and into loci with high or with minimal HpaII counts. This aligned using Illumina pipeline software version 1.3 or plot helped us to develop a new method for normalizing 1.4. A summary of the Illumina analysis results for each HpaII counts in terms of variability of the MspI represen- replicate is shown in Table S1 in Additional file 1. tation. We recognize that hypomethylated loci are associ- ated with relatively greater HpaII counts and a larger Data quality and reproducibility angle B (Figure 2b, left) whereas methylated loci will be Based on our experimental design, successfully generated defined by smaller angle values (Figure 2b, middle). Fur- products would be expected to possess a 5'-CGCTGCTG thermore, some loci will tend to be sequenced more read- sequence at the 3' end of the read, the first two nucle- ily than others, and may have identical B values but otides (CG) representing the cohesive end for ligation of differing distances from the origin (c distance), allowing a HpaII/MspI digestion products, the remaining six nucle- confidence score for identical methylation values (B) in otides the EcoP15I restriction enzyme recognition site. In terms of the c distance values (Figure 2b, right). To test order to evaluate the yield of desired products, we this model, we used bisulfite MassArray to test quantita- counted the number of reads containing this sequence tively the cytosine methylation values for 61 HpaII sites and plotted the starting positions of this sequence within (Tables S4, S5 and S6 in Additional file 1), choosing loci the reads obtained. We observed that approximately two- representing all components of the B angle spectrum of thirds of the reads contained the expected sequence, and values. In Figure 2c, d we show the correlations between found that the majority was located at base positions 25 these gold standard cytosine methylation values and raw and 26, consistent with the known digestion properties of HpaII counts or B angle values. We find that there is the the restriction enzyme [19]. Removal of the approxi- same negative correlation (R = 0.502) between HpaII mately 30% of reads lacking the CG-EcoP15I sequence counts and cytosine methylation values as demonstrated was performed to eliminate spurious sequences. In order in the MSCC technique [15], and that the angular trans- to investigate sequence quality further, we also deter- formation of the data incorporating the MspI normaliza- mined the number and relative position of Ns (ambigu- tion substantially improves this correlation (R = 0.826), ous base calls) within the reads obtained. Overall, few defining the optimal approach for processing of these reads were found to contain Ns, and where they were data. We represent the data for University of California present, they were found to be evenly distributed by posi- Santa Cruz (UCSC) genome browser visualization as wig- tion within the sequence. To test data reproducibility, we gle tracks, with higher B angle values defining less methy- compared the results of three experimental replicates lated loci. Methylated loci with zero values that would be against each other using the Pearson correlation coeffi- otherwise difficult to visualize as having been tested are cient metric. The results of this study showed that all rep- represented as small negative values. We show the details licates were highly correlated (all the r values exceed 0.9), of the analytical workflow in Figure 3 and an example of a which confirmed that the technical reproducibility of this UCSC genome browser representation of HELP-tagging assay was excellent (Table S2 in Additional file 1). data in Figure S1 in Additional file 1. All data are available through the Gene Expression Omnibus database (acces- Distribution of MspI/HpaII sequence tags sion number [GEO:GSE19937]) and as UCSC genome We merged three lanes of MspI data and observed that browser tracks [20]. approximately 80% of the 2,292,198 annotated HpaII sites in the human genome (hg18) were represented by at least Potential sources of bias: base composition and fragment one read, for a total of over 1.8 million loci throughout length the genome. The mean numbers of reads per locus for As the number of reads at CCGG sites following MspI MspI and HpaII were 3.94 and 1.82, respectively, and digestion should not be influenced by methylation, the MspI counts were distributed evenly across all genomic representation obtained from MspI digestion allowed us compartments examined (Table S3 in Additional file 1). to look for systematic sources of bias inherent to the We hypothesize that a combination of incomplete assay. A major concern was that base composition could genomic coverage and polymorphisms within some be a source of such bias, as it has been reported that Illu- CCGG sites (as we have previously observed [10]) mina sequencing can be influenced by GC composition accounts for the 20% of HpaII sites that were not repre- [21], possibly because of the gel extraction step [22]. Our sented by any reads. protocol does not require gel extraction and only begins to show an under-representation of sequences when the Normalization of HpaII by MspI counts and data (G+C) content exceeds approximately 80% (Figure 4a). transformation We also tested to see whether the sizes of the MspI frag- When we plot the MspI count on the x-axis and HpaII ments generated influenced the counts obtained, as the count on the y-axis for each HpaII site, we can see two Suzuki et al. Genome Biology 2010, 11:R36 Page 5 of 11 http://genomebiology.com/2010/11/4/R36 ELAND (1-36 bp alignment) Scan for EcoP15I tag No Discard the reads Yes Mask tag sequence as “n” ELAND (2-28 bp alignment) Quality failed/not aligned Discard the reads Unique/multiple aligned Distance from neighboring HpaII site < 27 bp Discard the reads ≥27 bp Multiple Unique aligned hits hit aligned >10 Discard the reads aligned ≤10 Weighted based on aligned number Map to annotated HpaII site? No Putative polymorphic HpaII sites dbSNP check Yes Count hit number MspI=0 Putative polymorphic HpaII sites dbSNP check MspI>0 Normalize HpaII with MspI by angle calculation UCSC genome browser Figure 3 HELP-tagging analysis workflow. The analysis workflow for HELP-tagging data is illustrated. Only sequence reads that contain the adapter sequence and map to a single or ≤ 10 sites are retained, the latter repetitive sequences distributed by weighting among the matched loci. Potential polymorphic loci are annotated. Normalization of HpaII by MspI using the angle calculation described in the previous figure is performed and files are generated for genome browser visualization. UCSC, University of California Santa Cruz. Suzuki et al. Genome Biology 2010, 11:R36 Page 6 of 11 http://genomebiology.com/2010/11/4/R36 Identification of polymorphic CCGG sequences Whereas MSCC used MmeI and generates an 18- to 19- (a) 0.12 bp sequence flanking the HpaII site [15], our use of Actual EcoP15I generates a 27-bp flanking sequence. We asked 0.10 MspI count whether this size difference influenced our ability to align sequences to the reference genome. We truncated our 0.08 sequence reads to 19 bp to mimic the MSCC read length 0.06 and found that this caused a profound loss of ability to align reads unambiguously (Table S7 in Additional file 1). 0.04 To compensate for the low alignment rate, the MSCC report described an ingenious strategy of alignment to 0.02 the sequences immediately flanking the annotated HpaII 0 sites in the reference genome [15], an approach suffi- ciently powerful that it generated the well-validated data 0 20406080 100 that they described. However, it does not offer the possi- GC content (%) bility of identifying polymorphic HpaII sites at the high frequencies that we previously observed for our HELP- (b) seq assay [10]. We tested whether our longer sequences 0.025 allowed the identification of loci at which an HpaII site is Actual annotated in the reference genome but we obtain no MspI count 0.020 sequence reads, and the opposite situation where we observed at least four MspI reads (the average number 0.015 per annotated MspI/HpaII site) flanking a locus not annotated in the reference genome. In Table S8 in Addi- 0.010 tional file 1 we list approximately 6,600 candidate poly- morphic HpaII sites, of which examples are shown in 0.005 Figure 5, confirmed by targeted resequencing of those loci. The 6,600 loci were selected based on overlap with dbSNP entries, allowing us to evaluate the pattern of 0 100 200 300 400 500 sequence variability at these loci. Approximately 80% of the SNPs are C:G to T:A transversions, consistent with Length (bp) deamination-mediated decay of methylcytosine being the Figure 4 Base composition and fragment length influences on se- cause of the polymorphism [23]. Polymorphic CG dinu- quence counts. (a) The proportion of (G+C) nucleotides was calculat- cleotides are major potential sources of error not only for ed for the 50-bp sequence centered around each annotated CCGG in microarrays, which are designed to a consensus genomic the reference human genome. The base composition of all of the MspI sequence, but also for both bisulfite sequencing, which sequences generated from the human ES cell line studied was also cal- culated. The relative proportion for (G+C) content in 2% bins for each would read the C to T transversion as unmethylated, and set of data was calculated and plotted as shown. The black line shows mass spectrometry-based assays, requiring the develop- the proportions in the reference genome, while the red line illustrates ment of specific analytical approaches such as we have the distribution we observed in our MspI experiment. Two peaks rep- described [24]. resenting base composition in repetitive sequences are apparent. The MspI distribution closely matches the expected distribution except DNA methylation studies of human embryonic stem cells when the base composition exceeds approximately 80%, when it is slightly under-represented. (b) We calculated the relative frequencies To test whether the HELP-tagging assay was generating of MspI digestion product sizes in the human reference genome. In data that are biologically plausible, we tested the methyla- this case we found that the shorter fragments are more likely to be se- tion of different genomic sequence compartments as den- quenced than larger (≥300 bp) fragments. The three major peaks ob- sity plots of B angle values for the human ES cells used in served represent Alu short interspersed repetitive element (SINE) these studies. In Figure S2a in Additional file 1 we show sequences. how promoters (defined as -2 kb to 2 kb from the tran- scription start site of RefSeq genes), gene bodies (the digestion by type III endonucleases like EcoP15I is most remaining region within the RefSeq gene) and intergenic efficient when a pair of enzymes is present in convergent (all other) sequences compare, finding the expected orientation on the same DNA molecule [19]. We find that enrichment of hypomethylated loci with larger B angle there is indeed an over-representation for shorter (≤300 values in promoter regions. When we compared unique bp) and a corresponding modest under-representation with repetitive sequences, again we found the expected for larger MspI fragments (Figure 4b). Counts relative to total Counts relative to total Suzuki et al. Genome Biology 2010, 11:R36 Page 7 of 11 http://genomebiology.com/2010/11/4/R36 (a) chr15: 25725400 25725500 25725600 25725700 25725800 25725900 MspI hit HpaII OCA2 CpG island dbSNP chr15: 25725640 25725650 25725660 25725670 25725680 A G T C T C T T CA C T C T CA CA T T C T A G C C C G G G G CT CC T G CC C A C A T T C T G C A T G G C A T G G C C T Reference HpaII OCA2 dbSNP rs12916836 rs12905726 A G TC T C T TCA C T C T CA CA T T C T A G C C C A A A GG C T CC T G CCC A C A C T C T G C A T GG C A T G G C C T Observed Trace data (b) 3870000 3875000 3880000 3885000 3890000 3895000 chr2: MspI hit HpaII dbSNP 3882300 3882310 3882320 3882330 3882340 3882350 chr2: G T C T G G A G C A G A G G C T T C T A A G CA CA GC A T C TT G G C C A A C G A A GC CA G C AC CAC A G G C A G G C A C T Reference MspI hit dbSNP rs6748872 G T C T G G A G CA G A G G C T TC T A A G CA CA G CA T G G C CA A C G A A G C C A G CA C CA CA G G C A G G C A C T Observed Trace data Figure 5 Polymorphic HpaII sites identified by HELP-tagging. Examples of HpaII sites (a) annotated in the reference genome sequence but not represented by MspI reads or (b) not annotated in the reference human genome and represented by at least four MspI reads are shown. In each case there is a SNP defined by dbSNP that indicates the C:T to G:A transversion that eliminates or restores the CCGG HpaII site. pattern of increased methylation of repetitive DNA com- that while transposable elements are generally methy- pared with unique sequences (Figure S2b in Additional lated and are depleted near gene promoters, those that file 1). Combining these observations, we tested whether are proximal to promoters tend to be less methylated the transposable element component of annotated repeti- than those located more distally. While many types of tive DNA sequences showed any tendency to unusual transposable elements were represented in this pro- methylation near gene promoters. In Figure 6 we show moter-proximal hypomethylated group, we found a sub- Suzuki et al. Genome Biology 2010, 11:R36 Page 8 of 11 http://genomebiology.com/2010/11/4/R36 (a) (b) 100% Angle 0-30 31-60 80% 61-90 60% 400 40% 200 20% 0 0% Length from TSS (bp) Length from TSS (bp) Figure 6 Identification of a position effect on DNA methylation in transposable elements located close to gene promoters. The distance from RefSeq gene transcription start sites and DNA methylation status are shown. The x-axis displays the distance from transcription start sites (TSSs). HpaII sites were categorized into three groups by angle, 0 to 30 (blue), 31 to 60 (red) and 61 to 90 (green)). (a) Number of HpaII sites; (b) proportions of each angle category (%). set to be the most markedly over-represented, as shown neous sites is that they contribute to the proportion of in Figure 6c. loci at which we could not obtain sequence reads. The outcome of these studies was an improvement in Our exploration of the distribution of cytosine methy- the previously described MSCC [15] and HELP-seq [10] lation in the same human ES cell line studied by Lister et assays, not only by means of technical modifications such al. [6] showed consistent results, with hypomethylation as the use of EcoP15I but also because of the concurrent of transcription start sites and methylation of transpos- use of MspI for normalization. The effect of these modifi- able elements, as expected from long-standing observa- cations was not only to increase the accuracy of the assay tions in the field. We furthermore discovered a limited but also to enhance the ability to align sequences to the subset of transposable elements that is hypomethylated genome and thus identify polymorphic HpaII/MspI sites. when in close proximity to transcription start sites. When The means of normalization of HpaII by MspI using an this subset was studied to determine whether certain angular metric is an innovation that improved the data types of transposable elements were disproportionately accuracy substantially and may have applications in other over-represented, we found two broad classes, one of MPS assay normalization strategies. We were also able to transposable element fossils with no innate capacity to discard reads that did not contain the expected adapter replicate themselves (the ancient DNA, long interspersed sequences, and created a straightforward data analytical repetitive elements (LINEs) and short interspersed repet- pipeline that will facilitate processing of these HELP-tag- itive elements (SINEs) shown in Figure 6c) and younger ging data by others. ERV1 long terminal repeat retroelements. Loss of methy- The potential sources of systematic artifacts due to base lation of functionally inactive transposable elements is composition or digestion product size were evaluated. likely to be of no negative consequence to the host Apart from a modest decrease in representation in genome, consistent with the host defense hypothesis [25], regions above approximately 80% (G+C) content, base while the young ERV long terminal repeats represent a composition did not cause biases in representations, pos- group of transposons whose function has been harnessed sibly in part due to our avoidance of a gel purification step as promoters of endogeneous genes [26,27]. This obser- in library preparation [22]. Fragment length does influ- vation demonstrates the value of a high-resolution, ence the outcome, most likely due to effects on EcoP15I genome-wide assay like HELP-tagging to define potential digestion [19], although the effects should be similar for functional elements in an unbiased manner. both HpaII and MspI and should, therefore, largely cancel each other out in the normalization step. It is possible Conclusions that endogeneous EcoP15I sites could influence the rep- We propose that MPS-based assays such as RRBS [12], resentations, but to have an effect they would have to be MSCC [15] and HELP-tagging will prove to be the assays located within the 27 bp adjacent to HpaII/MspI sites and of choice for epigenome-wide association studies in would cause digestion of the ligated adapter, causing human disease, with the latter two preferable as we begin those loci to be under-represented in both HpaII and to explore the CG-depleted majority of the genome. It MspI datasets. The most likely effect of these endoge- should not be necessary to run MspI assays every time a HELP-tagging assay is performed, suggesting that a com- -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 HpaII counts Percentage Suzuki et al. Genome Biology 2010, 11:R36 Page 9 of 11 http://genomebiology.com/2010/11/4/R36 mon MspI dataset can serve as a universal reference for a ments using a New England Biolabs Quick Ligation Kit species, allowing a single lane of Illumina sequencing of (25 μl of 2× Quick ligase buffer, 2.5 μl of adapter AS (10 the HpaII library to provide the methylation data for that μM), 2.5 μl of Quick Ligase in a final volume 50 μl). After sample. The development of analytical pipelines to sup- ligation, products were purified using the MinElute PCR port analysis of these datasets will be critical to the suc- purification kit (Qiagen, Hilden, Germany) and in vitro- cess of these projects, while the careful ongoing transcribed using the Ambion MEGAshortscriptkit (Life assessment of potential sources of bias will also be essen- Technologies, Carlsbad, CA, USA). Following in vitro tial for improving assay performance. transcription, products were purified with an RNeasy clean-up kit (Qiagen) before reverse transcription was performed using the Invitrogen SuperScript III kit (Life Materials and methods Cell preparation and DNA purification Technologies). The first strand cDNA produced was used H1 human ES cells (NIH code WA01 from Wicell as a template for PCR using the following conditions: Research Institute, Madison, WI, USA) were cultured on 96°C for 2 minutes, then 18 cycles of 96°C for 15 seconds matrigel (BD Biosciences, San Diego, CA, USA), at 37°C, and 72°C for 15 seconds followed by 5 minutes at 72°C for 5% O and 5% CO . Amplified human ES cell pluripo- the final extension. After PCR, the library was purified 2 2 using a QIAQuick PCR clean-up kit (Qiagen). tency was assessed by flow cytometry with SSEA4, CD24 and Oct4 markers. To extract DNA, the cells were sus- Single-locus quantitative validation assays pended in 10 ml of a solution of 10 mM Tris-HCl (pH Bisulfite conversion and MassArray (Sequenom, San 8.0), 0.1 M EDTA and 1 ml of 10% SDS to which 10 μl of Diego, CA, USA) were performed using an aliquot of the RNase A (20 mg/ml) was added. After incubation for 1 same sample of DNA as was used for the high-through- hour at 37°C, 50 μl of proteinase K (20 mg/ml) was added put assays described above. Bisulfite conversion was per- and the solution was gently mixed and incubated in a formed with an EZ DNA Methylation kit (Zymo 50°C water bath overnight. To purify the lysate, it was Research, Orange, CA, USA). Bisulfite primers were extracted three times using saturated phenol, then twice designed using MethPrimer [28], specifying the desired with chloroform, and dialyzed for 16 hours at 4°C against product length (250 to 450 bp), primer length (23 to 29 three changes of 0.2× SSC. Following dialysis, the DNA bp) and primer Tm (56 to 62°C). PCR was performed was concentrated by coating the dialysis bags in polyeth- using FastStart High Fidelity Taq polymerase (Roche, ylene glycol (molecular weight 20,000). The purity and Basel, Switzerland) with the following conditions: 95°C final concentration of the purified DNA was checked by for 10 minutes, then 42 cycles of 95°C for 30 seconds, spectrometry (Nanodrop, Wilmington, DE, USA). primer-specific Tm for 30 seconds and 72°C for 1 minute, followed by 72°C for 10 minutes for the final extension. Illumina library preparation Primer-specific Tm and sequence information are pro- The sample preparation steps are illustrated in Figure 1. vided in Table S6 in Additional file 1. Bisulfite MassArray Two custom adapters were created for HELP-tagging, assays were performed by the institutional Genomics referred to as AE and AS. As well as an Illumina adapter Core Facility. The data were analyzed using the analytical sequence, adapter AE contains an EcoP15I recognition pipeline we have previously described [24]. site and a T7 promoter sequence. Adapter AS contains an Illumina sequencing primer sequence. The adapter and Bioinformatic analysis primer sequences for library preparation are listed in Four lanes of sequencing were performed using an Illu- Table S9 in Additional file 1. Genomic DNA (5 μg) was mina GA IIx Sequencer at the institutional Epigenomics digested with HpaII and MspI in separate 200 μl reactions Shared Facility. Three lanes were used for technical repli- and purified by phenol/chloroform extraction followed cates of MspI, for the methylation-insensitive reference by ethanol precipitation. The digested genomic DNA was dataset. Images generated by the Illumina sequencer were ligated to adapter AE using a New England Biolabs Quick analyzed by Illumina pipeline software (versions 1.3 to Ligation Kit (25 μl of 2× Quick ligase buffer, 3 μg of 1.4). Initial data processing was performed using the HpaII-digested DNA or 1 μg of MspI-digested DNA, 0.1 default read length of 36 bp, after which we isolated the μl of Adapter AE (1 μM), 3 μl of Quick Ligase in a final sequences in which we found adapter sequences on the volume of 50 μl). After AE ligation, the products were 3'-end, replaced the adapter sequence with a poly(N) purified using Agencourt AMpure beads (Beckman sequence of the same length, and re-ran the Illumina Coulter, Brea CA, USA), then digested with EcoP15I ELAND pipeline again on these sequences with the (New England Biolabs). The restriction fragments were sequence length set at 27 bp (the 2 to 28 bp subsequence). end-repaired to inhibit to dimerization of adapters, and The data within the ELAND_extended.txt files were used tailed with a single dA, at the 3' end. After the dA tailing for counting the number of aligned sequences adjacent to reaction, adapter AS was ligated to the dA-tailed frag- Suzuki et al. Genome Biology 2010, 11:R36 Page 10 of 11 http://genomebiology.com/2010/11/4/R36 reduced representation bisulfite sequencing; UCSC: University of California each CCGG (HpaII/MspI) site annotated in the hg18 Santa Cruz. freeze of the human genome at the UCSC genome browser. We permitted up to two mismatches in each Authors' contributions MS and JMG designed the assays and strategies for its analysis, MS performed sequence, and allowed a sequence to align to up to a max- all library preparation and characterization, MS, DL and MP performed bisulfite imum of 10 locations within the genome. For non-unique validation studies, while QJ and AMcL performed computational analyses. JMG alignments, a sequence was assigned a partial count for and MS prepared the manuscript. each alignment location amounting to 1/n, where n rep- Acknowledgements resents the total number of aligned positions. To normal- This work is supported by a grant from the National Institute of Health (NIH, ize the data between experiments, the number of R01 HG004401) to JMG. The authors thank Shahina Maqbool PhD, Raul Olea sequences associated with each HpaII site was divided by and Gael Westby of the Einstein Epigenomics Shared Facility for their contribu- tions, Drs Eric Bouhassira and Emmanuel Olivier (Einstein) for the WA01/H1 the total number of sequences (including partial counts) human ES cell line, and Einstein's Center for Epigenomics. aligning to all HpaII sites in the same sample. We refer to this figure as the fixed count below. Author Details Department of Genetics (Computational Genetics), Center for Epigenomics, To examine an influence of (G+C) mononucleotide Albert Einstein College of Medicine, 1301 Morris Park Avenue, Bronx, NY 10461, content on counts of sequences obtained, we extracted USA the (G+C) annotation from the hg18 freeze of the human Received: 8 January 2010 Revised: 16 March 2010 genome at UCSC and examined the distribution of Accepted: 1 April 2010 Published: 1 April 2010 T © T Ge h h 2010 S i i no s s arti ime s an cle B u o iz o i p u ls o e k ava g n i e y acce 2010, t al. ilable ;ss arti lice f 11 rn o:R36 s m cle ee : h Bi dttp:/ istri oMbu / eg de te C nd e on u m tn ral L ed be io r th tlo d.g e y t .co erm ms o /2010/ f the 11/ Cre 4/ aR ti3 ve 6 Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited sequence counts according to (G+C) content. Annotated percentages of (G+C) content were available for adjacent References 1. Esteller M: Epigenetics in cancer. N Engl J Med 2008, 358:1148-1159. 5-bp windows. For each annotated HpaII site, we calcu- 2. Kopelovich L, Crowell JA, Fay JR: The epigenome as a target for cancer lated the mean percentage (G+C) for a 50-bp region cen- chemoprevention. J Natl Cancer Inst 2003, 95:1747-1757. tered at the restriction site. Counts of sequences 3. Zilberman D, Gehring M, Tran RK, Ballinger T, Henikoff S: Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an associated with HpaII sites were obtained for 50 sequen- interdependence between methylation and transcription. Nat Genet tial non-overlapping windows of 2% (G+C) (the mini- 2007, 39:61-69. mum possible in a sample of 50-bp regions). These data 4. Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE: Shotgun bisulphite were then normalized as a proportion of the total number sequencing of the Arabidopsis genome reveals DNA methylation of fragments. Comparisons were made to the expected patterning. Nature 2008, 452:215-219. frequencies, which for each 2% (G+C) bin was repre- 5. Lister R, O'Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, Ecker JR: Highly integrated single-base resolution maps of the sented by the counts of HpaII sites falling within a range epigenome in Arabidopsis. Cell 2008, 133:523-536. relative to the total number of HpaII sites in the genome. 6. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery This analysis was performed on both HpaII and MspI- JR, Lee L, Ye Z, Ngo QM, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar AH, Thomson JA, Ren B, Ecker JR: Human DNA methylomes at digested DNA for comparison. base resolution show widespread epigenomic differences. Nature The potential effect of distance between HpaII sites on 2009, 462:315-322. sequences counts obtained at each HpaII site was mea- 7. Heijmans BT, Tobi EW, Stein AD, Putter H, Blauw GJ, Susser ES, Slagboom PE, Lumey LH: Persistent epigenetic differences associated with sured by summing the counts of sequences aligning prenatal exposure to famine in humans. Proc Natl Acad Sci USA 2008, within each restriction fragment, and normalizing the 105:17046-17049. result with respect to total sequence count. As with the 8. Yamashita S, Hosoya K, Gyobu K, Takeshima H, Ushijima T: Development of a novel output value for quantitative assessment in methylated DNA (G+C) analysis above, this was performed for both MspI immunoprecipitation-CpG island microarray analysis. DNA Res 2009, and HpaII digested restriction fragments. The data were 16:275-286. compared with the expected distribution determined by 9. Down TA, Rakyan VK, Turner DJ, Flicek P, Li H, Kulesha E, Graf S, Johnson N, Herrero J, Tomazou EM, Thorne NP, Backdahl L, Herberth M, Howe KL, performing virtual restriction digestion using genomic Jackson DK, Miretti MM, Marioni JC, Birney E, Hubbard TJ, Durbin R, Tavare HpaII site coordinates, and normalizing the number of S, Beck S: A Bayesian deconvolution strategy for immunoprecipitation- virtual fragments of each size with respect to the total based DNA methylome analysis. Nat Biotechnol 2008, 26:779-785. 10. Oda M, Glass JL, Thompson RF, Mo Y, Olivier EN, Figueroa ME, Selzer RR, number of these virtual fragments. Richmond TA, Zhang X, Dannenberg L, Green RD, Melnick A, Hatchwell E, Bouhassira EE, Verma A, Suzuki M, Greally JM: High-resolution genome- Additional material wide cytosine methylation profiling with simultaneous copy number analysis and optimization for limited cell numbers. Nucleic Acids Res 2009, 37:3829-3839. Additional file 1 Supplemental data containing two figures (Figures S1 11. Ordway JM, Bedell JA, Citek RW, Nunberg AN, Jeddeloh JA: and S2) and nine tables ( Tables S1 to S9). MethylMapper: a method for high-throughput, multilocus bisulfite sequence analysis and reporting. Biotechniques 2005, 39:. 464, 466, 468 Abbreviations passim. bp: base pair; CG/CpG: cytosine-guanine dinucleotide; ChIP: chromatin immu- 12. Meissner A, Mikkelsen TS, Gu H, Wernig M, Hanna J, Sivachenko A, Zhang noprecipitation; ES: embryonic stem; (G+C): guanine and cytosine mononucle- X, Bernstein BE, Nusbaum C, Jaffe DB, Gnirke A, Jaenisch R, Lander ES: otides; HELP: HpaII tiny fragment Enrichment by Ligation-mediated PCR; MPS: Genome-scale DNA methylation maps of pluripotent and massively-parallel sequencing; MSCC: methyl-sensitive cut counting; RRBS: differentiated cells. Nature 2008, 454:766-770. Suzuki et al. Genome Biology 2010, 11:R36 Page 11 of 11 http://genomebiology.com/2010/11/4/R36 13. t Hoen PA, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RH, de Menezes RX, Boer JM, van Ommen GJ, den Dunnen JT: Deep sequencing- based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res 2008, 36:e141. 14. Brunner AL, Johnson DS, Kim SW, Valouev A, Reddy TE, Neff NF, Anton E, Medina C, Nguyen L, Chiao E, Oyolu CB, Schroth GP, Absher DM, Baker JC, Myers RM: Distinct DNA methylation patterns characterize differentiated human embryonic stem cells and developing human fetal liver. Genome Res 2009, 19:1044-1056. 15. Ball MP, Li JB, Gao Y, Lee JH, LeProust EM, Park IH, Xie B, Daley GQ, Church GM: Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells. Nat Biotechnol 2009, 27:361-368. 16. Fazzari MJ, Greally JM: Epigenomics: beyond CpG islands. Nat Rev Genet 2004, 5:446-455. 17. Hattori N, Nishino K, Ko YG, Hattori N, Ohgane J, Tanaka S, Shiota K: Epigenetic control of mouse Oct-4 gene expression in embryonic stem cells and trophoblast stem cells. J Biol Chem 2004, 279:17063-17069. 18. Irizarry RA, Ladd-Acosta C, Wen B, Wu Z, Montano C, Onyango P, Cui H, Gabo K, Rongione M, Webster M, Ji H, Potash JB, Sabunciyan S, Feinberg AP: The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific CpG island shores. Nat Genet 2009, 41:178-186. 19. Moncke-Buchner E, Rothenberg M, Reich S, Wagenfuhr K, Matsumura H, Terauchi R, Kruger DH, Reuter M: Functional characterization and modulation of the DNA cleavage efficiency of type III restriction endonuclease EcoP15I in its interaction with two sites in the DNA target. J Mol Biol 2009, 387:1309-1319. 20. Greallylab [http://greallylab.aecom.yu.edu/HELP-tagging/] 21. Dohm JC, Lottaz C, Borodina T, Himmelbauer H: Substantial biases in ultra-short read data sets from high-throughput DNA sequencing. Nucleic Acids Res 2008, 36:e105. 22. Quail MA, Kozarewa I, Smith F, Scally A, Stephens PJ, Durbin R, Swerdlow H, Turner DJ: A large genome center's improvements to the Illumina sequencing system. Nat Methods 2008, 5:1005-1010. 23. Duncan BK, Miller JH: Mutagenic deamination of cytosine residues in DNA. Nature 1980, 287:560-561. 24. Thompson RF, Suzuki M, Lau KW, Greally JM: A pipeline for the quantitative analysis of CG dinucleotide methylation using mass spectrometry. Bioinformatics 2009, 25:2164-2170. 25. Yoder JA, Walsh CP, Bestor TH: Cytosine methylation and the ecology of intragenomic parasites. Trends Genet 1997, 13:335-340. 26. Conley AB, Piriyapongsa J, Jordan IK: Retroviral promoters in the human genome. Bioinformatics 2008, 24:1563-1567. 27. Cohen CJ, Lock WM, Mager DL: Endogenous retroviral LTRs as promoters for human genes: a critical assessment. Gene 2009, 448:105-114. 28. MethPrimer [http://www.urogene.org/methprimer/] doi: 10.1186/gb-2010-11-4-r36 Cite this article as: Suzuki et al., Optimized design and data analysis of tag- based cytosine methylation assays Genome Biology 2010, 11:R36

Journal

Genome Biology – Springer Journals

Published: Apr 1, 2010

There are no references for this article.