Cell and tissue type independent age-associated DNA methylation changes are not rare but common

Most age-associated DMPs are shared between blood cell subtypes

First, we aimed to demonstrate that age-associated DMPs (aDMPs) derived from different cell or tissue types do exhibit a highly statistically significant overlap. We note that a highly significant overlap would be a strong cue that aDMPs shared between relevant cell or tissue-types is the norm and not the exception. In order to avoid the confounding effect of cell-type heterogeneity we first focused on three DNAm datasets profiling purified blood cell subtypes. Specifically, we analysed Illumina 450k DNAm data of 1199 purified CD14+ monocyte and 214 purified CD4+ T-cell samples from Reynolds et al. [17], and of 100 purified CD8+ T-cells from Tserel et al. [18]. For each dataset, we derived age-associated DMPs (aDMPs), adjusting for potential technical confounders (Methods). We used two different significance thresholds: a false discovery rate (FDR) < 0.05, and a Bonferroni-adjusted P_bonf < 0.05 thresholds. The former admits 5% of false positives but ensures that the FNR is not abnormally high, whereas the Bonferroni threshold ensures in principle no false positive but at the expense of a much larger FNR. Using an FDR < 0.05 threshold, we observed a very strong and highly statistically significant overlap of aDMPs between all 3 cell-types (Fig. 1A). For instance, almost 4000 aDMPs were found to be shared between all 3 cell-types using an FDR < 0.05 threshold (Fig. 1A). The high statistical significance of the overlaps remained if Bonferroni thresholds were used (Fig. 1B). Thus, even though “only” 198 aDMPs were in common between all 3 cell-types when using Bonferroni thresholds, this number was highly significant, i.e. much higher than what would be expected if aDMPs were cell-type specific.

Figure 1. Most age-DMPs are shared between blood cell subtypes. (A) Landscape overlap diagram for age-DMPs defined using FDR<0.05 threshold in three separate purified blood cell subtype datasets (1199 Monocytes from Reynolds et al, 214 CD4+ T-cells from Reynolds et al and 100 CD8+ T-cells from Tserel et al). Barplots indicate the number of aDMPs found in each purified cell-type, or the corresponding aDMP overlap between cell-types. For the overlapping categories, the P-value of the overlap is indicated by the color of the bar, as shown. (B) As A), but now using a Bonferroni corrected P<0.05 threshold. (C) Smoothed scatterplots of t-statistics of age-DMPs called in a given blood cell subtype vs. the corresponding t-statistics in another cell subtype, as indicated for the 3 pairwise comparisons. In each panel, we only depict the 100 most outlier data points, provide the number of probes in each significant quadrant and the P-value assessing agreement is from a one-tailed Fisher-test. The vertical red lines indicate the line of FDR<0.05, whilst the horizontal lines depict the “validation threshold” of P<0.05.

An alternative to estimating significance of overlaps, is to evaluate the consistency of the t-statistics between aDMPs, which therefore also takes into account the directionality of DNAm change. This can be done by generating scatterplots of the t-statistics of selected aDMPs in one dataset against the corresponding t-statistics in another dataset, the rationale being that if an aDMP declared in one dataset is not an aDMP in another, then its t-statistic in this other set will be distributed with a mean value of 0. Thus, if aDMPs are cell-type specific and only valid in one dataset, their t-statistics in the other set profiling a different cell-type should form a data cloud centered around 0. Applying this strategy to the three purified blood cell type datasets above, revealed in each case that t-statistics were strongly skewed away from 0 and very strongly positively correlated (Fisher-test P-values <1e-100, Fig. 1C).

In order to validate these findings, we analyzed independent data from a multi cell-type EWAS from the BLUEPRINT consortium [19], encompassing Illumina 450k data from CD4+ T-cells, CD14+ monocytes and CD16+ neutrophils for a total of 139 individuals. We note that the matched nature of this dataset naturally adjusts for age-range and sample size, since all 3 cell-types are available for each of the 139 individuals. We performed exactly the same analysis as before, which confirmed that overlaps were highly statistically significant, using either a FDR or Bonferroni-based threshold (Fig. 2A-B), and importantly, that there was also a very strong correlation between the t-statistics of corresponding aDMPs (Fisher test P-values < 1e-100, Fig. 2C), further supporting the view that a large fraction of aDMPs in one blood cell-type constitute aDMPs in another.

Figure 2. Most age-DMPs are shared between blood cell subtypes: validation in BLUEPRINT. (A) Landscape overlap diagram for age-DMPs defined using FDR<0.05 threshold in the matched multi cell-types DNAm dataset from BLUEPRINT (139 monocyte samples, 139 matched CD4+ T-cell samples and 139 matched neutrophil samples. Barplots indicate the number of aDMPs found in each purified cell-type, or the corresponding aDMP overlap between cell-types. For the overlapping categories, the P-value of the overlap is indicated by the color of the bar, as shown. (B) As (A), but now using a Bonferroni corrected P<0.05 threshold. (C) Smoothed scatterplots of t-statistics of age-DMPs called in a given blood cell subtype vs. the corresponding t-statistics in another cell subtype, as indicated for the 3 pairwise comparisons. In each panel, we only depict the 100 most outlier data points, provide the number of probes in each significant quadrant and the P-value assessing agreement is from a one-tailed Fisher-test. The vertical red lines indicate the line of FDR<0.05, whilst the horizontal lines depict the “validation threshold” of P<0.05.

Many age-associated DMPs are shared between tissue types

Next, we decided to compare aDMPs across different tissues. These comparisons are particularly problematic if tissues derive from different datasets that profiled variable numbers of samples, with different age-ranges and subject to potentially different confounding factors, all of which could greatly impact on statistical significance estimates [20]. Thus, ideally, a cross-tissue comparison should include multiple tissue samples from the same set of individuals, all profiled as part of the same study. Therefore, we analyzed Illumina 850k data from an EWAS profiling blood, buccal and cervical samples from a common set of 263 women (Methods). Because blood is a complex mixture of many immune-cell (IC) subtypes, and buccal and cervical samples are highly contaminated by immune cells [15], we identified aDMPs in each tissue after adjustment for batch effects and cell-type heterogeneity using EpiDISH [21] (in the case of blood) and HEpiDISH [15] (in the case of buccal and cervix). Although aDMPs identified with or without cell-type correction were highly correlated (Fig. 3A), we observed that adjusting or not for cell-type heterogeneity did have a marked impact on the number of aDMPs, and that the number of aDMPs also varied substantially between tissue-types (Fig. 3B). Of note, using either an FDR < 0.05 or Bonferroni adjusted P-value < 0.05 thresholds, the overlap of aDMPs between the 3 tissues was highly significant (Fig. 3C, P<1e-100), mimicking the result obtained on blood cell subtypes. For instance, we observed a total of 2200 aDMPs in common between blood, buccal and cervix, an overlap which cannot be explained by random chance (Fig. 3C, P<1e-100). Scatterplots of t-statistics of aDMPs between tissues further supported an extremely strong correlation, suggesting that shared aDMPs between blood, buccal and cervix is the norm and not the exception (Fig. 3D-F).

Figure 3. Most age-DMPs are shared between tissue-types. (A) Smoothed scatterplots of age-DMP (aDMP) t-statistics obtained after adjustment for cell-type heterogeneity (x-axis) against the corresponding t-statistics without adjustment (y-axis), for 3 different tissue-types (Blood, Buccal and Cervix) separately. In each scatterplot, we only depict the 100 most outlier data points, provide the number of probes in each significant quadrant and the P-value is from a one-tailed Fisher-test. The vertical red lines indicate the line of FDR<0.05, whilst the horizontal lines depict the “validation threshold” of P<0.05. (B) Barplot of the number of aDMPs (FDR<0.05) in each tissue-type before and after adjustment for cell-type heterogeneity. (C) Venn-diagrams representing the number of overlapping aDMPs between the 3 tissue-types using an FDR<0.05 threshold for calling aDMPs (top panel) and using a Bonferroni threshold (lower panel). P-value as estimated using a nested Hypergeometric test. (D) Smoothed scatterplots of aDMP t-statistics (adjusted for cell-type heterogeneity) in blood (x-axis) against their corresponding t-statistics in buccal (y-axis) for top panel, with lower panel depicting the reverse analysis, as indicated. (E-F) As (D), but for the combinations blood-cervix and buccal-cervix, respectively.

Studies profiling only a few hundred samples are underpowered to detect most age-DMPs

The previous analyses clearly support the view that at least thousands of aDMPs are shared between two and three tissue/cell-types. However, is this lower bound a significant underestimate on the true number of aDMPs that are shared between any two given cell or tissue types? To address this question requires careful consideration of the expected power of the studies. To estimate empirically the expected power as a function of sample size, we devised a subsampling strategy using the Reynolds monocyte dataset (Methods). We reasoned that this set, due to its large size (n=1199) and purified nature of the cell population, would allow us to objectively define a gold-standard set of aDMPs in monocytes. Using a stringent Bonferroni-adjusted P<0.05 threshold and using all 1199 monocyte samples, we thus defined a gold-standard set of 18596 monocyte aDMPs. We note that although this is certainly only a small subset of all true monocyte aDMPs, that it would nevertheless allow us to assess the impact of sample size on power. We next subsampled 100 samples from the 1199, and re-derived a new set of aDMPs at the same Bonferroni threshold. This subsampling strategy was performed for increasing subsampling size (from 100 to 1000, in units of 100), and a total of 10 times for each subsampling size. Sensitivity to detect the 18596 gold-standard aDMPs was computed and plotted against subsampling size, revealing that for sample sizes on the order of 100 or 200, the sensitivity was very low (Fig. 4A). Indeed, for a subsample of 100, we estimated a mean sensitivity of only 0.00001, for a subsample of size 200 the mean sensitivity was 0.005, for 300 the value was 0.04, and at a value of 600 (i.e. about half of the full set) the sensitivity was still only 0.35. Thus, in light of this, if we were to now compare aDMPs across different tissue-types with some of the corresponding datasets in the order of 100-200 samples, as done in [12], then even if all aDMPs were shared between tissue-types, we would never be able to detect much overlap and we would wrongly conclude that most aDMPs are “dataset-specific”, i.e. tissue-specific in the context of the analysis presented in [12]. For instance, using a Bonferroni threshold we observed an overlap of “only” 213 aDMPs between the 18596 gold-standard monocyte aDMPs and the 291 aDMPs derived from the 214 CD4+ T-cell samples (Fig. 1B), and so we would be inclined to conclude that most aDMPs are cell-type specific. However, assuming that all aDMPs are shared between monocytes and CD4+ T-cells, our subsampling analysis (Fig. 4A) would suggest that the expected sensitivity to detect the 18596 gold-standard monocyte aDMPs using the 214 T-cell samples would be a value between 0.005 (n=200) and 0.04 (n=300) (Fig. 4A), with the value closer to 0.005. Indeed, under this shared aDMP scenario, the detected overlap of 213 aDMPs represents a sensitivity fraction estimate of 213/18596 ≈ 0.01, in line with our subsampling estimate. Thus, the observed overlap of aDMPs at the given sample sizes of the two studies is not inconsistent with the great majority of aDMPs being shared between monocytes and T-cells. We note that this result is also highly congruent with the statistical significance estimates obtained previously via the Fisher-test (Fig. 1B-C).

Figure 4. Empirical Power Analysis. (A) Boxplot of the sensitivity (y-axis) to detect gold-standard aDMPs, defined at the full purified monocyte sample size n=1199, at random subsampling sizes (x-axis), as indicated. Each boxplot displays the sensitivity over 10 separate runs, in each run aDMPs were defined at the Bonferroni 0.05 level. (B) As A), but now defining aDMPs in each run as those with an FDR<0.05. Because the FDR estimates are more unstable, we performed 100 runs at each subsampling size. (C) As A), but now defining the gold-standard set of aDMPs by imposing a threshold on the effect size (2% DNAm change over 10 years), in addition to a Bonferroni adjusted P-value < 0.05. At each subsampling run, aDMPs were also defined using the same criterion and 10 runs were performed at each subsampling size. (D) As A), but now only using the threshold on the effect size to define gold-standard aDMPs and to define aDMPs each subsampling size. A total of 100 runs at each subsampling size were performed.

We further note that using a more relaxed FDR < 0.05 threshold, sensitivity would be substantially higher: at about n=600, sensitivity is already close to 1, and for 300 samples, sensitivity is over 0.4 (Fig. 4B). Thus, when comparing aDMPs between multiple cell or tissue-types it is even more critical to use FDR-based thresholds, since otherwise using Bonferroni-based adjustment, the expected overlap of aDMPs derived from say 4 separate studies will be zero, even if all aDMPs are common to all 4 cell/tissue-types. Our analysis suggests that many hundreds if not thousands of samples would be needed to ensure that overlaps over 3 or more studies would have the appropriate sensitivity to detect the majority of shared aDMPs (Fig. 4A-B). We verified that all these results are independent of whether an additional threshold on the effect size is used to select aDMPs (Fig. 4C). Indeed, using an additional and identical threshold on the effect size to define aDMPs as used in [12], i.e demanding at least a 2% change in DNAm over 10 years in addition to Bonferroni significance, yielded higher sensitivities, but at sample sizes of 100 and 200, the expected sensitivity was still only 0.0002 and 0.05, respectively (Fig. 4C).

Using only a threshold on effect size for feature selection suffers from strong selection bias

If we were to ignore statistical significance estimates (which depend on sample size) altogether, and instead rank features by effect size using the above mentioned threshold (a 2% DNAm change over 10 years) to declare aDMPs, we can see that sensitivities increase substantially, exhibiting a much lower dependency on sample size (Fig. 4D). For instance, at n=100, the sensitivity would be as high as 0.6 (Fig. 4D). At first, this seems to support the argument by Slieker et al. that an observed lack of overlapping aDMPs defined via an effect size threshold would imply that aDMPs are largely tissue-specific. However, this argument is problematic, primarily for two reasons. First, selecting aDMPs based on effect size still suffers from selection bias, i.e. the fact that in the dataset where aDMPs are selected effect sizes will naturally be higher than in independent studies. This selection bias arises in real data because of numerous study-specific confounders which can significantly inflate or deflate effect sizes. Indeed, our subsampling analysis shows that the sensitivity to detect aDMPs at a sample size of 100 to be approximately 40% lower than at the full sample size (n=1199) (Fig. 4D), which is a substantial reduction given that the subsample derives from the same dataset. While this also demonstrates that using an effect size to rank and select aDMPs does not guarantee that the ranking and selection is independent of sample size, as claimed in [12], we stress that the selection bias will be even more pronounced when comparing across independent studies. Thus, the observation made by Slieker et al. [12] that the effect sizes of aDMPs selected from one dataset appear reduced in another set profiling a different tissue could easily be the result of selection bias, and nothing to do with the nature of the different tissue being profiled.

A second major problem associated with using an effect size threshold to select and validate aDMPs is related to confounding factors such as cell-type heterogeneity, which may artificially deflate effect sizes in spite of associations with age remaining highly significant. To demonstrate this, we posited that the fraction of aDMPs derived in the monocyte set that would validate in a large whole blood set [10] (which contains monocytes) would be much reduced if an effect size criterion is used throughout, as compared to using a statistic and P-value. Confirming this, out of the 844 aDMPs with an effect size larger than 2% over 10 years in the Reynolds monocyte set, only about 40% passed this same threshold across the 656 whole blood samples from Hannum et al. [10] (Fig. 5A-B). While the fraction validating in whole blood was similar if a Bonferroni threshold is used (fraction was 41%), the fraction doubled to 84% if an FDR<0.05 threshold was used instead (Fig. 5B). Thus, using only an effect size threshold to select and evaluate overlap of aDMPs between studies could lower sensitivities by as much as another 40% in relation to using an FDR-based threshold. Of note, using effect size thresholds in the original beta-value basis, which is highly heteroscedastic, may also strongly bias aDMPs towards those with intermediate DNAm values. Indeed, we verified that selecting aDMPs using the 2% change over 10-year threshold resulted in very few or no aDMPs with mean DNAm values close to 0 or 1 (Fig. 5C). Since the mean DNAm level of a CpG site in a complex tissue may be strongly influenced by cell-type heterogeneity, using an effect size threshold to select or validate aDMPs could therefore easily miss a large fraction of true aDMPs.

Figure 5. Pitfalls of using a threshold on effect-size only to select aDMPs. (A) Boxplots comparing the effect size distribution for the 844 aDMPs defined in the Reynolds monocyte set against their effect sizes in the whole blood dataset from Hannum et al. (B) Barplots comparing the fraction of aDMPs, defined either by the effect size threshold (blue) or P-value threshold (red & magenta, Bonferroni-adjusted), that validate in the whole blood set from Hannum et al. In Hannum et al, validated aDMPs were defined either as those passing the same effect size threshold (blue), or the same Bonferroni-threshold (red), or a more relaxed FDR<0.05 threshold (magenta). The numbers below the bars indicate the absolute number of aDMPs validating in Hannum et al. This panel demonstrates that using the same effect size threshold to define aDMPs in a dataset of complex tissue samples could miss up to 40% of true aDMPs. (C) Comparison of the density profiles of the average DNAm for the 844 aDMPs defined by having an effect size larger (in absolute magnitude) than 0.002 (equivalent to a 2% DNAm change over 10 years) (blue line) across the 1199 monocyte samples of Reynolds et al, versus the corresponding density profile of the 18596 gold-standard aDMPs with Bonferroni adjusted P-values < 0.05.

In summary, the implicit assumption made by Slieker et al. [12] that the effect size of an aDMP should not depend on sample size and on other factors such as cell-type heterogeneity of the tissue or other confounders, appears to not hold and would lead to the false conclusion that a lack of overlapping aDMPs, all selected using effect size thresholds, is due to a lack of shared aDMPs. In fact, our analysis above clearly indicates that using only effect sizes to select aDMPs could result in a severe overall selection bias, with sensitivities reduced by as much as 80% if not more.

FDR analysis suggests that most of the DNA methylome is altered with age

A clue as to how many aDMPs are cell/tissue specific can also be inferred from the FDR characteristics in the largest datasets. For the 1199 monocyte sample set from Reynolds et al. we used the estimated FDR values (q-values) to further estimate that on average only 52419 of the 482,091 probes (i.e. 10%) are not associated with age, suggesting therefore that approximately 90% of the DNA methylome is altered with age. Because FDR estimates are notoriously sensitive to confounding factors, it is important to check the consistency of the FDR estimates in Reynolds et al. against those obtained in smaller studies, ideally profiling identical or related cell-types. For the 201 monocyte samples profiled with Illumina 450k as part of BLUEPRINT, we estimated 41% of the probes to be aDMPs, whereas the fraction of aDMPs was similar, around 45%, for the 104 monocyte samples (52 twin pairs) from Paul et al. [22]. In the case of whole blood, in Hannum et al. which encompassed 656 samples and therefore approximately half of the numbers in the Reynolds monocyte set, FDR values yielded an estimate of approximately 66% aDMPs. This value is very close to the one we estimated for the 689 whole blood samples profiled in Liu et al. [5]: the fraction of probes estimated there to be aDMPs was 68%. For a smaller set such as the 263 whole blood samples from FORECEE, FDR values yielded a correspondingly lower estimate of only 12% aDMPs.

Thus, assuming that the fraction of aDMPs (~90%) in the Reynolds monocyte set is inflated due to some confounder, it is unlikely to be inflated by more than 25%, since the fraction of predicted aDMPs in two large whole blood studies was consistently over 65%. This supports the view that a very high fraction (likely to be well over 65%) of the DNA methylome of blood cell subtypes is altered with age, and therefore this also means that there would be a large overlap of aDMPs between any two blood cell subtypes, consistent with our observations. Moreover, the fact that we observe equally strong overlaps of aDMPs between tissues like blood and cervix, further suggests that relatively high fractions of the DNA methylome of other cell-types (e.g. epithelial and fibroblasts) are also altered with age.

Most of Horvath’s clock CpGs are pan-tissue aDMPs

Finally, to demonstrate that there are indeed many examples of aDMPs that are shared between tissues, we analyzed in detail the 353 CpGs that make up Horvath’s clock [9]. Specifically, we computed their t-statistics and P-values using linear models against age, including potential covariates as confounding factors, across a total of 10 different cell/tissue types, which included Monocytes, CD4+ and CD8+ T-cells, B-cells, Neutrophils, Buccal, Cervix, Liver, Brain and Fibroblasts (Methods). In the case of buccal, cervix and liver we adjusted for variations in the epithelial, fibroblast and immune cell fractions, whereas in the case of brain we adjusted for variations in neuronal and non-neuronal fractions (Methods, Table S1). Although only 3% (i.e 10 CpGs) of the 353 clock CpGs were aDMPs in all 10 tissues, approximately 78% were aDMPs across at least 3 different cell/tissue types, with only 7% appearing to be “tissue or cell-type specific” (Fig. 6, Table S2). The 10 CpGs defining aDMPs in all 10 tissue/cell types mapped to genes that included VGF, GRIA2, FZD9, KLF14, RHBDD1, KCNC2, NHLRC1, P2RX6 and CECR6. Of note, FZD9 is a transmembrane receptor for Wnt signaling proteins, whilst GRIA2 is a glutamate neurotransmitter receptor, both of which have previously been demonstrated to be part of interactome “hotspots” of age-associated DNAm which occur independently of cell or tissue-type [23]. Thus, this confirms that ELOVL2 is not unique and that many of Horvath’s clock CpGs constitute aDMPs across several cell/tissue types.

Figure 6. Pan-tissue analysis of Horvath’s clock CpGs. Heatmap displays the signed P-values of the 353 Horvath Clock CpGs across 10 independent cell or tissue-types, where the P-values derive from a linear model of DNAm against age plus additional confounders as covariates. Blue denotes highly significantly age-associated hypermethylation (hyperM), red denotes highly significant hypomethylation (hypoM). The 353 clock CpGs have been ranked according to the number of cell/tissue types where they are age-DMPs (using FDR<0.05), indicated by horizontal black bars at top. The cumulative proportion of the 353 CpGs that are age-DMPs in 10, 9, 8, 7, 6, 5, 4, 3, 2 or 1 studies are shown as vertical bars.

DNAm datasets

Identification of age-DMPs

In each dataset we used linear models with the DNAm value as the response variable and with age of the sample as the predictor. Depending on the dataset, linear models were run with additional covariates to account for potential confounding factors. In the case of purified cell samples, covariates included those provided by the publications which included batch or ethnicity information. In the case of Reynolds et al., age-DMPs were derived by linear regression on 482127 (CD4T) and 482091 (Monocytes) probes, with gender and race as covariates (which dominated variation as determined by a PCA). In the case of the CD8+ T-cells, age-DMPs were derived by linear regression on 472484 probes, with gender and array number as covariates (which dominated variation as determined by a PCA). In the case of Blueprint data, age-DMPs were derived by linear regression on 473719 probes, with gender and batch number as covariates (which dominated variation as determined by a PCA).

In the case of complex tissues, besides adjusting for batch effects (if these were present), we also corrected for cell-type heterogeneity. Briefly, in the case of whole blood, we used our previously validated DNAm reference matrix for blood with our EpiDISH algorithm [21] to obtain cell-type fraction estimates for the 7 main immune cell subtypes: neutrophils, basophils, monocytes, B-cells, NK-cells, CD4+ and CD8+ T-cells. In the case of other tissues, like buccal and cervix, we used the corresponding DNAm reference matrix from our HEpiDISH algorithm [15] to obtain cell-type fractions for the total epithelial, total fibroblast and the 7 main immune cell subtypes. In the case of liver, we derived aDMPs using sex, body-mass index, cohort and cell-type fractions as covariates. In each case, the estimated cell-type fractions were used as covariates in the linear models. In the case of brain, we derived aDMPs using sex, plate, position and brain cell fraction (neuron and non-neuron) as covariates.

Age-DMPs (aDMPs) were generally defined at two distinct thresholds: at a false discovery rate (FDR) threshold less than 0.05, where FDR values were estimated using the q-value procedure [38], and using a much more conservative Bonferroni threshold (P < 0.05/n with n the number of probes for which the linear model was run).

Subsampling power analysis in Reynolds Monocyte set

We used the large (n=1199) purified monocyte sample set from Reynolds et al. to define a gold-standard list of 18596 monocyte aDMPs using a Bonferroni threshold. We then subsampled 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000 samples from the original 1199 and redefined aDMPs at each subsampling size using the same Bonferroni threshold. At each subsample size we estimated the sensitivity to detect the 18596 aDMPs from the full set. A total of 10 different runs were performed at each subsample size. We also derived aDMPs and sensitivities at each subsample size but now using an FDR<0.05 threshold. Because the FDR estimation is more unstable, we performed 100 different Monte-Carlo runs at each subsampling size.

The whole analysis above was repeated by defining aDMP by the additional requirement, that the effect size (i.e. slope) is larger than 2% over 10 years, i.e. a slope value of absolute magnitude larger than 0.002, which is the effect size threshold used in [12]. In a final analysis, we repeated the procedure but now defining aDMPs using only the threshold on the effect size, discarding statistics and P-values.

Data availability

All data analyzed in this manuscript are publicly available from EGA https://ega-archive.org/ accession number EGAS00001001456 and GEO (http: www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE117370, GSE56581, GSE56046 and GSE59065.