Skip to main content

A Bayesian mixed modeling approach for estimating heritability

Abstract

Background

A Bayesian mixed model approach using integrated nested Laplace approximations (INLA) allows us to construct flexible models that can account for pedigree structure. Using these models, we estimate genome-wide patterns of DNA methylation heritability (h2), which are currently not well understood, as well as h2 of blood lipid measurements.

Methods

We included individuals from the Genetics of Lipid Lowering Drugs and Diet Network (GOLDN) study with Infinium 450 K cytosine-phosphate-guanine (CpG) methylation and blood lipid data pre- and posttreatment with fenofibrate in families with up to three-generation pedigrees. For genome-wide patterns, we constructed 1 model per CpG with methylation as the response variable, with a random effect to model kinship, and age and gender as fixed effects.

Results

In total, 425,791 CpG sites pre-, but only 199,027 CpG sites posttreatment were found to have nonzero heritability. Across these CpG sites, the distributions of h2 estimates are similar in pre- and posttreatment (pre: median = 0.31, interquartile range [IQR] = 0.16; post: median = 0.34, IQR = 0.20). Blood lipid h2 estimates were similar pre- and posttreatment with overlapping 95% credibility intervals. Heritability was nonzero for treatment effect, that is, the difference between pre- and posttreatment blood lipids. Estimates for triglycerides h2 are 0.48 (pre), 0.42 (post), and 0.21 (difference); likewise for high-density lipoprotein cholesterol h2 the estimates are 0.61, 0.68, and 0.10.

Conclusions

We show that with INLA, a fully Bayesian approach to estimate DNA methylation h2 is possible on a genome-wide scale. This provides uncertainty assessment of the estimates, and allows us to perform model selection via deviance information criterion (DIC) to identify CpGs with strong evidence for nonzero heritability.

Background

Narrow-sense heritability (h2), traditionally estimated using twins or other constrained family relationships, can also be estimated in wider pedigrees using a linear mixed-model approach [1, 2]. This approach is well-established for traits having moderate to high h2 [2, 3]. What is less clear is how well these models perform for traits with low h2, which is the case for some proportion of cytosine-phosphate-guanine (CpG) methylation sites genome wide. This is of particular interest for h2 because it is a ratio of variance components, that is, the proportion of phenotypic variance explained by additive genetic variance. Compared to estimates of the mean, estimates of variance tend to have large uncertainty.

DNA methylation (DNAm) is an epigenetic mark that is implicated in many heritable diseases and traits. DNAm patterns are influenced by the environment, change over lifetime, and exhibit mitotic heritability [4]. There are multiple reports of DNAm h2 with CpG-specific estimates ranging from 0 to 1, and with genome-wide mean or median estimates ranging from approximately 15 to 30% [5,6,7].

Any point estimates of h2 should ideally be accompanied by an assessment of the uncertainty. A Bayesian mixed model with inference by integrated nested Laplace approximations (INLA) [2, 8] presents an improvement on the traditionally used simple, balanced designs (eg, parent–offspring or twin-based correlations) as it can include the full pedigree structure. Additionally, full posterior inference provides uncertainty for the estimated h2, which is not obtainable using frequentist methods, and allows an assessment of whether the true h2 is nonzero via comparison of models with and without a genetic component.

In this article, we use three-generation pedigrees from the Genetics of Lipid Lowering Drug and Diet Network (GOLDN) [9] to estimate h2 in a fully Bayesian framework. Blood lipid and DNAm measurements were available pre- and posttreatment with a lipid-lowering drug (fenofibrate). We first investigate continuous blood lipid traits with known moderate to high h2 and compare the Bayesian estimates with previously published frequentist estimates [10]. Second, we estimate the genome-wide distribution of DNAm h2.

The Bayesian mixed-model approach provides posterior distributions and therefore uncertainty for h2, the additive genetic variance, and the environmental variance. We further investigate if it is reasonable to estimate h2 for each CpG site given the data by performing model selection, to assess if the genetic component should be included at all. When the additive genetic variance is truly zero, a linear mixed model including the family structure will misspecify the model, making any resultant h2 estimates untrustworthy. When h2 is truly nonzero but smaller than the minimum h2 detectable in a given data set (ie, when the study is underpowered), we would expect extremely large standard errors for the h2 estimates.

In this article, our aims are to estimate h2 for (a) blood lipid traits and (b) genome-wide DNAm using a Bayesian approach, including a model selection step to identify CpG sites with strong evidence for nonzero heritability. We highlight the importance of the model selection step and of quantifying the uncertainty in h2 estimates.

Methods

Data

Individuals in the GOLDN study having DNAm and blood lipid measurements, as well as information on gender and age, were included in the analysis [11]. Genome-wide DNAm was measured using the 450 K Infinium array at 463,995 CpG sites and was available for 995 and 530 individuals at pre- and postfenofibrate treatment, respectively. From the original 463,995 sites available, we removed a list of known single nucleotide polymorphism–related CpG sites. Blood lipid levels were available both pre- and posttreatment for 818 and 861 individuals for triglycerides (TGs) and high-density lipoproteins (HDLs), respectively. In our analyses, we use log-transformed TG and HDL (averaged over 2 measurements taken 1 day apart), and simply refer to these transformed average values as TG and HDL from here on.

Standard quality control inspection of the DNAm data showed systematic differences with respect to probe-type chemistry. We therefore normalized the data with respect to the probe type using beta-mixture quantile normalization using the R package wateRmelon [12]. We calculated the M values for each probe and used these in the analysis.

Statistical model

We estimate h2 in a Bayesian framework using the INLA package in Rue et al. [8]. Using the following model, we investigate h2 for TG, HDL, and DNAm in the GOLDN study pre- and posttreatment:

$$ {y}_i\mid {\eta}_i,{\sigma}_e^2\sim N\left({\eta}_i,{\sigma}_e^2\right)\forall i=1,\dots, n\ \mathrm{with}\ {\eta}_i={\beta}_0+{\boldsymbol{x}}_i^T\boldsymbol{\beta} +{u}_i $$
$$ \boldsymbol{u}\mid {\sigma}_g^2\sim N\left(\mathbf{0},2{\sigma}_g^2\boldsymbol{K}\right)\ \mathrm{with}\ \boldsymbol{u}=\left({u}_1,\dots, {u}_n\right) $$
(1)

where y = (y1, …, yn) is the outcome vector, that is, DNAm for a given CpG, TG, or HDL, for n subjects. \( N\left({\eta}_i,{\sigma}_e^2\right) \) is the likelihood function for subject i, where ηi is the linear predictor defining the latent variables of interest and \( {\sigma}_e^2 \) denotes the residual unexplained variance in the model. Fixed effects are the intercept β0 and β = (β1, β2) corresponding to the effects of age (x1i) and gender (x2i). Family structure was modeled with u as a random effect, with distribution\( N\left(\mathbf{0},2{\sigma}_g^2\boldsymbol{K}\right) \), where K is the kinship matrix calculated based on pedigree information using the {kinship2} package in R [13], and \( {\sigma}_g^2 \) is the additive genetic variance. We assume default noninformative priors for all hyperparameters for the genome-wide investigation of h2 (\( \mathrm{in}\ \mathrm{particular}\ {\sigma}_{\ast}^{-1}\sim Gamma\left(\mathrm{1,0.00005}\right) \)), but other priors were also investigated. From this model,\( {h}^2=\frac{\sigma_g^2}{\sigma_g^2+{\sigma}_e^2} \).

With the deviance information criterion (DIC) [14, 15], we can investigate the fit of a model to the data and do model selection. We can compare the DIC of the 2 latent models:

$$ {\eta}_i={\beta}_0+{\boldsymbol{x}}_i^T\boldsymbol{\beta} $$
(2)
$$ {\eta}_i={\beta}_0+{\boldsymbol{x}}_i^T\boldsymbol{\beta} +{u}_i $$
(3)

and use the recommended ΔDIC = DIC(2) − DIC(3) > 10 to choose if we should include the genetic component (ui) or not.

As a result of the setup of INLA, we can only obtain point estimates of h2 with this model formulation, but not the marginal posterior distribution. If this is of interest, the following equivalent model formulation must be used instead [2]. The nuisance parameter \( {\sigma}_0^2 \) is introduced to ensure that \( {\sigma}_e^2 \) remains identifiable. It is fixed at a small value, so that it does not absorb much of the variation in the data and thus does not interfere with \( {\sigma}_g^2 \) or \( {\sigma}_e^2 \). We used \( {\sigma}_0^2=4.54\mathrm{e}-5 \), which corresponds to log(σ20) = − 10 and was suggested in an earlier publication [2].

$$ {y}_i\mid {\eta}_i,{\sigma}_e^2\sim N\left({\eta}_i,{\sigma}_0^2\right)\forall i=1,\dots, n\ \mathrm{with}\ {\eta}_i={\beta}_0+{\boldsymbol{x}}_i^T\boldsymbol{\beta} +{u}_i+{\varepsilon}_i $$
$$ \boldsymbol{u}\mid {\sigma}_g^2\sim N\left(\mathbf{0},2{\sigma}_g^2\boldsymbol{K}\right)\mathrm{with}\ \boldsymbol{u}=\left({u}_1,\dots, {u}_n\right) $$
$$ \boldsymbol{\varepsilon} \mid {\sigma}_e^2\sim N\left(\mathbf{0},{\sigma}_e^2\boldsymbol{I}\right)\ \mathrm{with}\ \boldsymbol{\varepsilon} =\left({\varepsilon}_1,\dots, {\varepsilon}_n\right) $$

Results

Figure 1 shows the marginal posterior distribution of h2 for HDL and TG levels pre- and posttreatment, as well as for the change (post–pre). The change in blood lipid levels can be interpreted as the response to treatment. The model selection step supported the inclusion of ui in the model and thus estimation of h2 for each case. The mode estimates of h2 for TG are 0.48, 0.42, and 0.21 for baseline, posttreatment, and change, respectively. Likewise, the mode estimates for h2 for HDL are 0.61, 0.68, and 0.10. Based on point estimates alone, HDL appears to be more heritable than TG but their 95% credible intervals overlap. TG is more heritable pretreatment than posttreatment, but the opposite is true for HDL. However, the 95% credibility intervals for h2 pre- and posttreatment overlap each other substantially for both traits.

Fig. 1
figure 1

Heritability (h2) estimates for high-density lipoprotein (HDL) and triglycerides (TG). Posterior h2 distribution for HDL and TG pre- and posttreatment, as well as for the change from pre- to posttreatment (“Diff”). The points indicate the 95% credibility region for each h2 estimate

In Fig. 2, the histogram of point estimates of h2 of DNAm for 448,040 CpG sites is displayed, for both pre- and posttreatment. Model selection indicates if there is evidence for nonzero h2.

Fig. 2
figure 2

Genome-wide heritability (h2) estimates for DNA CpG methylation (DNAm). Histogram of point estimates for DNAm h2 genome wide. The top panel shows the h2 estimates pretreatment and the lower panel shows the h2 estimates posttreatment. The blue color is for CpG sites with strong evidence for heritability (“delta DIC >10”) and the red color is for the CpG sites with no evidence for heritability given the data (“delta DIC ≤10”), based on model selection comparing the models with and without pedigree information [eqs. (2) and (3)]. The two colors are stacked on top, so that both colors are visible on each bar. The actual height of the zero-stack bar for the posttreatment is 238,165 CpG sites

For pretreatment, we detect nonzero h2 for 425,791 of 448,040 CpG sites (95.0%). For the CpG sites with strong evidence for nonzero h2, the mean, median, and interquartile range (IQR) are 0.33, 0.31, and 0.16. In contrast, for posttreatment we detect nonzero h2 for 199,027 of 448,040 CpG sites (44.4%) with mean, median, and IQR equal to 0.36, 0.34, and 0.20. Both cases show evidence of zero inflation and a right skewed distribution of DNAm h2. Indirectly, we see that GOLDN is underpowered for detecting h2 of DNAm less than approximately 10 to 15%.

Figure 2 clearly shows a striking difference in the distribution of genome-wide h2 estimates between pre- and posttreatment. Figure 3 shows the direct pairwise comparison of the two h2 estimates for DNAm across those CpG sites with strong evidence for nonzero h2. For these CpG sites, we observe a moderate Pearson correlation between the baseline and posttreatment heritability estimates of 0.53.

Fig. 3
figure 3

Contour plot of correlation of DNA CpG methylation heritability (h2) pre- and posttreatment. The plot is based on the 194,741 CpG sites that have evidence for nonzero heritability in both pre- and posttreatment methylation measurements

Discussion and conclusions

In this article, we use a Bayesian approach to estimate h2 for blood lipid levels and DNAm in extended pedigrees, and show that it is feasible to do genome-wide Bayesian estimation of h2.

For blood lipids, our pretreatment point estimates for h2 using the Bayesian approach (h2TG = 0.48 and h2HDL = 0.61) are comparable to previously reported estimates using other methods [16]. Notably, the 95% credibility intervals are wide (see Fig. 1). Interestingly, the response to treatment (change in blood lipid levels from pre- to posttreatment; see Fig. 1) also appears to be heritable. This may have important implications for personalized medicine, especially if one were able to identify explanatory variables (genetic or environmental) that predict the individual genetic value (ui).

The genome-wide distribution of DNAm h2 (as ascertained by the Illumina Infinium 450 K chip) appears to follow a 2-group mixture model with a h2 = 0 component and a nonzero component following a smooth, unimodal right skewed distribution. The mixture proportion is strikingly different in the pre- and posttreatment (see Fig. 2). It is unlikely that treatment with fenofibrate could cause substantial genome-wide changes in DNAm h2. It is much more likely that technical differences between the DNAm time points explain the difference [17].

The h2 estimates we present in this article are generally similar to earlier genome-wide results in different populations [5, 7]. These previously published studies also show zero inflation, although the reported proportions vary substantially, from 5 to 17%. All show right-skewed distributions, but with slightly different modality. We report the median, h2pre = 0.31, h2post = 0.34, in the GOLDN study only for CpG sites with strong evidence for nonzero h2 (ie, ΔDIC > 10). This is higher than previously reported, although this is to be expected if previous reports included the zero component in the calculation of the median. With the zero component included, our median estimates decrease to h2pre = 0.30, h2post = 0.0036.

A Bayesian model holds some advantages over a frequentist model when estimating h2. The most notable difference is that a Bayesian model can estimate the full posterior distribution of the parameters in question. This includes the random effect variance parameters, which allows more direct assessment of the uncertainty of the random-effect estimates—and thus of the h2 estimates—than could be achieved with classical mixed-effects models. For example, by comparing posterior credibility intervals for h2 pre- versus posttreatment, we can determine if the h2 changes are a result of treatment (see Fig. 1). In addition, model selection can be performed in a principled way by DIC or Bayes factors to assess if the data provide sufficient evidence for nonzero h2. When the h2 is truly zero, the genetic values ui should not be included in the model, and therefore a model without the genetic values should be favored. Indeed, DIC-based model selection favored the model without random effect in the vast proportion of CpG sites for which h2 was estimated to be close to zero.

An important limitation in our study is the lack of information on shared households to estimate the effect of shared environment. This likely results in an overestimation of h2. The CpG sites on the 450 K array are a nonrandom selection of CpG sites with a bias toward CpG islands that often have regulatory function. This is important to consider, as an oversampling of biologically important CpG sites could result in a biased estimation for nonspecific genome-wide DNAm h2.

In this article, we use a Bayesian approach to estimate h2 in extended pedigrees and show the importance of using model selection to determine if there is strong evidence of nonzero h2 given the data. Future work should focus on simulation studies comparing the frequentist and Bayesian approaches, with particular focus on traits having low h2..

References

  1. Almasy L, Blangero J. Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998;62(5):1198–211.

    Article  CAS  Google Scholar 

  2. Holand AM, Steinsland I, Martino S, Jensen H. Animal models and integrated nested Laplace approximations. G3 (Bethesda). 2013;3(8):1241–51.

    Article  Google Scholar 

  3. Ralston SH, Uitterlinden AG. Genetics of osteoporosis. Endocr Rev. 2010;31(5):629–62.

    Article  CAS  Google Scholar 

  4. Slatkin M. Epigenetic inheritance and the missing heritability problem. Genetics. 2009;182(3):845–50.

    Article  Google Scholar 

  5. McRae AF, Powell JE, Henders AK, Bowdler L, Hemani G, Shah S, Painter JN, Martin NG, Visscher PM, Montgomery GW. Contribution of genetic variation to transgenerational inheritance of DNA methylation. Genome Biol. 2014;15(5):R73.

    Article  Google Scholar 

  6. Quon G, Lippert C, Heckerman D, Listgarten J. Patterns of methylation heritability in a genome-wide analysis of four brain regions. Nucleic Acids Res. 2013;41(4):2095–104.

    Article  CAS  Google Scholar 

  7. Kulkarni H, Kos MZ, Neary J, Dyer TD, Kent JW Jr, Goring HH, Cole SA, Comuzzie AG, Almasy L, Mahaney MC, et al. Novel epigenetic determinants of type 2 diabetes in Mexican-American families. Hum Mol Genet. 2015;24(18):5330–44.

    Article  CAS  Google Scholar 

  8. Rue H, Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J R Stat Soc Series B Stat Methodol. 2009;71:319–92.

    Article  Google Scholar 

  9. Irvin MR, Kabagambe EK, Tiwari HK, Parnell LD, Straka RJ, Tsai M, Ordovas JM, Arnett DK. Apolipoprotein E polymorphisms and postprandial triglyceridemia before and after fenofibrate treatment in the genetics of lipid lowering and diet network (GOLDN) study. Circ Cardiovasc Genet. 2010;3(5):462–7.

    Article  CAS  Google Scholar 

  10. Weiss LA, Pan L, Abney M, Ober C. The sex-specific genetic architecture of quantitative traits in humans. Nat Genet. 2006;38(2):218–22.

    Article  CAS  Google Scholar 

  11. Das M, Irvin MR, Sha J, Aslibekyan S, Hidalgo B, Perry RT, Zhi D, Tiwari HK, Absher D, Ordovas JM, et al. Lipid changes due to fenofibrate treatment are not associated with changes in DNA methylation patterns in the GOLDN study. Front Genet. 2015;6:304.

    Article  Google Scholar 

  12. Pidsley R, Wong CCY, Volta M, Lunnon K, Mill J, Schalkwyk LC. A data-driven approach to preprocessing Illumina 450K methylation array data. BMC Genomics. 2013;14:293.

    Article  CAS  Google Scholar 

  13. Therneau T, Atkinson E, Sinnwell J, Matsumoto M, Schaid D, Mcdonnell S: Package “kinship2” [Online] 2012. Available at: http://cran.r-project.org/web/packages/kinship2/kinship2.pdf.

  14. Spiegelhalter DJ, Best NG, Carlin BR, van der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc Series B Stat Methodol. 2002;64:583–616.

    Article  Google Scholar 

  15. Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian Data Analysis, vol. 2. Boca Raton: Chapman & Hall/CRC; 2014.

    Google Scholar 

  16. Kathiresan S, Manning AK, Demissie S, D’Agostino RB, Surti A, Guiducci C, Gianniny L, Burtt NP, Melander O, Orho-Melander M, et al. A genome-wide association study for blood lipid phenotypes in the Framingham heart study. BMC Med Genet. 2007;8(Suppl 1):S17.

    Article  Google Scholar 

  17. LeBlanc M, Nustad HE, Zucknick M, Page CM. Quality control for Illumina 450K methylation data in the absence of iDat files using correlation structure in pedigrees and repeated measures. BMC Genet. 2018;19(Suppl 1) https://doi.org/10.1186/s12863-018-0636-5.

Download references

Acknowledgements

Part of this work was performed on the Abel Cluster, owned by the University of Oslo. The authors also acknowledge input from Robert Lyle and Kristina Gervin.

Funding

Publication of this article was supported by NIH R01 GM031575. HEN received funding from the Norwegian Research Counsel. CMP and ML received internal funding from Oslo University Hospital.

Availability of data and materials

The data that support the findings of this study are available from the Genetic Analysis Workshop (GAW) but restrictions apply to the availability of these data, which were used under license for the current study. Qualified researchers may request these data directly from GAW.

About this supplement

This article has been published as part of BMC Proceedings Volume 12 Supplement 9, 2018: Genetic Analysis Workshop 20: envisioning the future of statistical genetics by exploring methods for epigenetic and pharmacogenomic data. The full contents of the supplement are available online at https://bmcproc.biomedcentral.com/articles/supplements/volume-12-supplement-9.

Author information

Authors and Affiliations

Authors

Contributions

HEN, ML, MZ, and CMP conceived and designed the experiment; HEN, ML, AHR, and CMP generated the results; and HEN, ML, MZ, and CMP interpreted the results. All authors wrote the manuscript and approved the final version.

Corresponding author

Correspondence to Marissa LeBlanc.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Nustad, H.E., Page, C.M., Reiner, A.H. et al. A Bayesian mixed modeling approach for estimating heritability. BMC Proc 12 (Suppl 9), 31 (2018). https://doi.org/10.1186/s12919-018-0131-z

Download citation

  • Published:

  • DOI: https://doi.org/10.1186/s12919-018-0131-z

Keywords