Cubic regression-based degree of correction predicts the performance of whole bisulfitome amplified DNA methylation analysis.

Epigenetic mechanisms, including DNA methylation, are important determinants in development and disease. There is a need for technologies capable of detecting small variations in methylation levels in an accurate and reproducible manner, even if only limited amounts of DNA are available (which is the case in many studies in humans). Quantitative methylation analysis of minute DNA amounts after whole bisulfitome amplification (qMAMBA) has been proposed as an alternative, but this technique has not been adequately standardized and no comparative study against conventional methods has been performed, that includes a wide range of methylation percentages and different target assays. We designed an experiment to compare the performance of qMAMBA and bisulfite-treated genomic (non-amplified) DNA pyrosequencing. Reactions were performed in duplicate for each technique in eight different target genes, using nine artificially constructed DNA samples with methylation levels ranging between 0% and 100% with intervals of 12.5%. Cubic polynomial curves were plotted from the experimental results and the real methylation values and the resulting equation was used to estimate new corrected data points. The use of the cubic regression-based correction benefits the accuracy and the power of discrimination in methylation studies. Additionally, dispersion of the new estimated data around a y = x line (R ( 2) ) served to fix a cutoff that can discriminate, with a single 9-point curve experiment, whether whole bisulfitome amplification and subsequent qMAMBA can produce accurate methylation results. Finally, even with an optimized reagent kit, DNA samples subjected to whole bisulfitome amplification enhance the preferential amplification of unmethylated alleles, and subtle changes in methylation levels cannot be detected confidently.

Introduction

DNA methylation plays a crucial role in normal development and has relevant effects on gene regulation, X chromosome inactivation and genomic stabilization through transcriptional silencing of repetitive elements. In humans, DNA methylation occurs mostly in cytosines within CpG dinucleotides, which are underrepresented in the genome because of their high mutation rate, but are especially abundant in 5′ regulatory regions of genes.

Common complex diseases with strong genetic and environmental risk determinants are excellent targets for the study of gene methylation, which is known to link environmental factors and the genome. Moreover, most common complex diseases have a progressive and quantitative nature that could be partly explained by accumulating epigenetic variation. For instance, differences in DNA methylation levels have been identified in the frontal cortex of patients diagnosed with major psychotic symptoms, in genomic regions that had been previously linked to disease etiology.

Studies in monozygotic twins have demonstrated purely epigenetic variation underlying complex traits, as is the case of the decrease in methylation in several genes associated with systemic lupus erythematosus. SNP associations detected by genome wide association (GWA) studies have been revisited in order to identify parent of origin effects or differentially methylated transcription factor binding sites and have been successful in the case of different diseases such as type 1 and type 2 diabetes. It is now believed that the integration of epigenomic and genomic analyses will help reveal novel functional determinants within associated loci and several associations between epigenomic perturbations and human diseases have been reported.

On the other hand, several studies have provided evidence on the existence of hundreds of genes with differential methylation of their promoters in various types of cancers, and such genes have been shown to be consistent and reproducible across different data sets and publications. Global hypomethylation of the genome and hypermethylation at specific gene promoters are known to be partly responsible for important events in carcinogenesis, including chromosomal instability and transcriptional silencing of tumor suppressor genes., Although robustly replicated across studies, these methylation changes are often quantitatively subtle, as in breast cancer, where the absolute differences between tumoral and normal tissues are below 5%. Thus, it is necessary to develop technologies capable of determining slight variations in methylation levels in an accurate and reproducible manner.

Clinical tissue specimens from particular stages of a disease are limited and often irreplaceable, making those DNA samples very valuable material for research. Quantitative methylation analysis of minute DNA amounts after whole bisulfitome amplification (qMAMBA) has been proposed as a solution for limited sample availability (for example, in the case of DNA extracted from body fluids)., Taking into account that acceptable methods for epigenetic studies should ensure that subtle methylation differences are robustly detected, we designed a study to compare qMAMBA and bisulfite-treated DNA pyrosequencing in terms of accuracy and reproducibility.

In this work, we validate the previously proposed cubic polynomial regression correction method and use it to fix a cutoff, based on a single 9-point dilution curve experiment that can determine whether whole bisulfitome amplification and subsequent qMAMBA technology can be performed without producing an artifactual result.

Results and Discussion

Following the experimental design outlined in Figure 1, we compared the performance of conventional genomic DNA pyrosequencing and qMAMBA in a methylated DNA dilution curve (0%, 12.5%, 25%, 37.5%, 50%, 62.5%, 75%, 87.5% and 100% methylated human DNA) for eight Pyromark methylation assays located in different genes (MALT1, MAP3K7, MAP3K7IP1, MAP3K14, NFKBIA, RELA, TNFAIP3 and TRADD), that were selected because of their central regulatory functions in NFκB and Toll-like receptor responses, and their interest in studies on inflammation and autoimmunity. A 9-point cubic polynomial curve was plotted with the results from each of the target genes, with real methylation levels in the x axis and experimentally obtained, Pyromark Q24 software derived values for each point of the curve in the y axis. Curves for all genes fitted very well a polynomial cubic equation (R2 > 0,975) (). For each experimental result, x was solved to correct for PCR bias, and estimated data presented b values very close to 1, confirming the high level of correction of PCR bias that is obtained with the cubic equation based correction, compared with the raw pyrosequencing results (Table 1). When these estimated values were plotted against the real methylation levels, linear regression showed very high coincidence with a y = x straight line (Fig. 2) and data-points fitted very well to the function, as shown by correlation coefficients (Table 1), especially for the conventional pyrosequencing experiments.

Figure 1. Schematic representation of the experimental design.

Table 1. Raw and corrected b values, sums of squared errors of the differences between estimated and real methylation levels, and coefficients of variation between replicates are shown for each gene and methodology, expressed as average values of the 9 points of the curve and standard deviations

Gene	Experimentalapproach	Raw b	b after correction	Sum of squared error	Mean CV	R2
MALT1	qMAMBA	0.25 ± 0.12	1.01 ± 0.61	603.68	0.15 ± 0.16	0.9410
	PYROSEQUENCING	0.55 ± 0.12	0.97 ± 0.22	100.35	0.10 ± 0.15	0.9886
MAP3K7	qMAMBA	0.48 ± 0.22	1.00 ± 0.30	224.42	0.24 ± 0.13	0.9719
	PYROSEQUENCING	1.71 ± 0.41	1.00 ± 0.21	78.12	0.07 ± 0.05	0.9955
MAP3K7IP1	qMAMBA	0.42 ± 0.12	1.02 ± 0.20	26.35	0.11 ± 0.14	0.9960
	PYROSEQUENCING	0.78 ± 0.23	0.96 ± 0.19	44.71	0.08 ± 0.10	0.9948
MAP3K14	qMAMBA	0.63 ± 0.25	1.17 ± 0.42	278.96	0.07 ± 0.10	0.9690
	PYROSEQUENCING	0.67 ± 0.19	0.96 ± 0.22	34.17	0.09 ± 0.11	0.9961
NFKBIA	qMAMBA	0.25 ± 0.07	0.95 ± 0.32	327.06	0.08 ± 0.12	0.9654
	PYROSEQUENCING	0.67 ± 0.14	1.02 ± 0.14	30.96	0.03 ± 0.04	0.9964
RELA	qMAMBA	0.18 ± 0.06	0.98 ± 0.28	196.92	0.07 ± 0.04	0.9747
	PYROSEQUENCING	0.62 ± 0.06	1.02 ± 0.10	18.49	0.04 ± 0.03	0.9979
TNFAIP3	qMAMBA	0.20 ± 0.09	0.98 ± 0.27	483.03	0.10 ± 0.05	0.9555
	PYROSEQUENCING	0.70 ± 0.12	1.00 ± 0.10	21.94	0.04 ± 0.05	0.9975
TRADD	qMAMBA	0.44 ± 0.12	0.95 ± 0.19	62.28	0.08 ± 0.05	0.9933
	PYROSEQUENCING	0.77 ± 0.11	1.01 ± 0.12	17.09	0.06 ± 0.04	0.9971

Figure 2. Degree of bias and correction in each experiment. In all graphs, the x axis represents the real methylation levels, while the experimental (black diamonds and squares) and estimated (gray triangles) values are on the y axis. The derived cubic polynomial curve is shown for one of the sets of experimental results (diamonds), presenting in most cases a preferential amplification of unmethylated alleles (b < 1, except for MAP3K7 with conventional pyrosequencing). Gray triangles plot the real values vs. the levels estimated after the cubic equation-based correction and are fitted to a y = x linear regression (discontinuous line). Linear correlation coefficients are shown for each experiment.

Assays with a correlation coefficient higher than 0.99 in the estimated vs. real value linear regression (Table 1) were able to discriminate clearly between all adjacent points of the 9-point methylation curves. This was the case in 7 out of 8 target genes analyzed with the conventional pyrosequencing method but only in 2 of the genes assayed with qMAMBA. In the plots with lower correlation values, several adjacent points of the dilution curve (differing in 12.5% methylated DNA) were scored with similar estimated values and could not be discriminated, as highlighted in Figure 2.

Sum of squared errors were calculated in order to quantify the accumulated difference between real and estimated methylation values with different target genes and techniques. In general, sum of squared error values were notably higher in the qMAMBA results. On the other hand, variation coefficients and their standard deviations across the different curve points were also calculated and showed more variation between replicates in the case of qMAMBA results (Table 1).

As previously mentioned, in order to be able to assess the contribution of variation in DNA methylation to human disease, subtle changes in methylation levels must be robustly detected in the wider range of disorders as possible using small amounts of DNA. Technical approaches such as qMAMBA could be very useful in methylation analyses, but they should prove to be reliable enough. Previous analyses of whole bisulfitome amplified samples addressed only extreme differences in methylation levels, which indeed are correctly discriminated., The present work is the first study to compare qMAMBA and bisulfite-treated genomic DNA pyrosequencing in terms of accuracy and reproducibility in a systematic way, covering a wide range of methylation percentages and different target assays.

In general, our study shows that the performance of assays is worse when bisulfitome amplification precedes PCR and pyrosequencing reactions, suggesting an inflation of PCR bias in this technique, due to a very strong preferential amplification of unmethylated alleles in qMAMBA. The only exception is MAP3K7IP1, which shows a very similar performance by both conventional pyrosequencing and qMAMBA. In this work, we demonstrate that qMAMBA method can only be satisfactorily applied if the R2 value derived from the linear regression for the corrected vs. the real methylation values is very close to 1 (> 0.99) and thus, dispersion of estimations from the straight line is small enough to avoid stretches of overlapping values in which methylation differences cannot be accurately calculated (Fig. 2). Although less frequent, the same precautions should also be taken with conventional pyrosequencing.

In conclusion, we use the previously proposed cubic polynomial regression correction method and provide evidence on its efficiency in correcting the PCR bias. Moreover, we show that b values depart much more from 1 when calculated with uncorrected observed methylation values compared with corrected estimations, so that systematic use of the cubic regression-based correction enormously benefits both the accuracy and the power of discrimination in methylation studies.

Finally, we fix a cutoff that allows distinguishing between those assays that will have a satisfactory performance from those unable to discriminate subtle changes in methylation levels, by performing a single 9-point curve experiment. This threshold establishes that in order to avoid artifactual results, the correlation coefficient (R2) should be higher than 0.99 when a linear regression between the estimated and the real methylation values is plotted, following cubic polynomial equation correction of the raw results.

Methods

DNA samples

Methylation curves were constructed combining different amounts of commercially available HeLa Genomic DNA (used as unmethylated DNA) and CpG Methylated HeLa Genomic DNA (representing totally methylated human genomic DNA), both from New England Biolabs (cat. nos. N4006S and N4007S, respectively) (Fig. 1). The resulting methylation levels of the curve where the following: 0%, 12.5%, 25%, 37.5%, 50%, 62.5%, 75%, 87.5% and 100%, all in a final volume of 20 µl and a DNA concentration of 100 ng/µl.

Sodium bisulfite conversion of unmethylated cytosines was performed in each sample of the curve using the commercially available Epitect Bisulfite Kit (Qiagen cat. no. 59104) as indicated by the manufacturer, and column-purified into 40 µl. An aliquot (5 µl) of each sample was used to amplify the bisulfitome with the Epitect Whole Bisulfitome Amplification Kit (Qiagen cat. no. 59205) following the manufacturer’s instructions. Briefly, around 250 ng bisulfite treated DNA were amplified using 1 µl of REPLI-g Midi Polymerase in a final volume of 40 µl. Reactions were incubated at 28°C for 8h, followed by a 95°C step of 5 min and then kept overnight at 4°C. Finally, 30 µl each of bisulfite treated genomic and amplified samples were diluted 1:2 in ddH2O and used for PCR reactions.

Selection and pyrosequencing of the CpG islands

Eight predesigned, commercially available Pyromark CpG Assays (Qiagen) located in eight different genes (MALT1, MAP3K7, MAP3K7IP1, MAP3K14, NFKBIA, RELA, TNFAIP3 and TRADD) were selected. Each of the assays covered 4–6 CpG islands located in the 5′ non coding end and/or in the first exon of each target, regions known to be involved in gene expression regulation and thus, interesting candidates for subsequent disease related studies (Table 2).

Table 2. Pyromark CpG Assay list

Pyromark CpG assay target name	Cat. no.
MALT1	PM00185143
MAP3K14	PM00177569
MAP3K7	PM00122850
MAP3K7IP1	PM00199521
NFKBIA	PM00056287
RELA	PM00048895
TNFAIP3	PM00122129
TRADD	PM00061369

Each of the pyrosequencing experiments was performed in duplicate, so that two sets of results were available for each target gene and approach. Amplifications of the loci of interest were performed in a C1000 Thermal Cycler (BioRad) using 1.5 µl of DNA in 25-µl PCR reactions with the Pyromark PCR Kit (Qiagen cat. no. 978703) and the forward and the biotinylated reverse specific primers provided with each Pyromark CpG Assay, following the manufacturer’s protocol. Thermocycling conditions were as follows: 15 min at 95°C for enzyme activation followed by 45 cycles of denaturation for 30 sec at 95°C, annealing for 30 sec at 56°C and extension for 30 sec at 72°C, with a final extension of 10 min at 72°C. PCR reactions were tested by agarose gel electrophoresis.

After PCR amplification, 20 µl of each biotinylated PCR product were bound to streptavidin coated sepharose beads (GE Healthcare cat. no. 17-5113-01) in 24-well plates. PCR products were denatured and non-biotinylated strands were removed using the Pyromark Q24 vacuum workstation. The beads were then resuspended in 25 μl of annealing buffer containing 0.3 μM of the specific sequencing primer from each CpG Assay. Pyrosequencing with the Pyromark Q24 system was performed using the Pyro Gold Q24 reagents. Nucleotide, enzyme and substrate dispensation order and volumes, as well as the duration of each step, were planned using the Pyromark Q24 software v2.0 (all from Qiagen) based on the target sequence. Following the runs, raw data files were imported into the software for standard quality controls and preliminary analyses.

Data analysis and Statistics

For each target gene and experimental approach (direct pyrosequencing of genomic DNA or qMAMBA), results were fitted into a polynomial cubic equation (y = ax3 + cx2 + dx + e), with the real methylation levels plotted on the x axis and the observed, experimental values (extracted from de Pyromark Q24 software) on the y axis, using Microsoft Excel v.10.0. With this equation, methylation levels were estimated for each experimental result by solving the x with Cardano’s method, as previously described. Adjustment and accuracy of these estimations with respect to the real methylation levels were measured by linear correlation coefficients (R2).

The value of b, a measure of amplification bias that reflects the efficiency of primer binding and polymerase elongation in amplicons with differential composition of methylated and unmethylated DNA, was calculated for each curve point with the following equation: b = [y × (100 - x)]/[x × (100 - y)], where y is the uncorrected experimental or estimated and x is the real value. As described previously, b = 1 indicates unappreciable PCR bias, while preferential amplification of unmethylated alleles results in b < 1.

Average b values (representing an overall b value of each particular gene), sums of squared errors of the differences between estimated methylation values and real levels, and coefficients of variation between replicates and their standard deviations were calculated for each gene and approach using Prism v5.0 (GraphPad Software Inc.).

Click here for additional data file.