- Open Access
A Bayesian mixed modeling approach for estimating heritability
- Haakon E. Nustad†1, 2, 3,
- Christian M. Page†4, 5,
- Andrew H. Reiner4,
- Manuela Zucknick6 and
- Marissa LeBlanc4Email author
© The Author(s). 2018
- Published: 17 September 2018
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.
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.
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.
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.
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 . 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–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)  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 . 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.
Individuals in the GOLDN study having DNAm and blood lipid measurements, as well as information on gender and age, were included in the analysis . 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 . We calculated the M values for each probe and used these in the analysis.
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%.
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 . 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 .
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..
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.
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.
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.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis 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.
- Almasy L, Blangero J. Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998;62(5):1198–211.View ArticleGoogle Scholar
- Holand AM, Steinsland I, Martino S, Jensen H. Animal models and integrated nested Laplace approximations. G3 (Bethesda). 2013;3(8):1241–51.View ArticleGoogle Scholar
- Ralston SH, Uitterlinden AG. Genetics of osteoporosis. Endocr Rev. 2010;31(5):629–62.View ArticleGoogle Scholar
- Slatkin M. Epigenetic inheritance and the missing heritability problem. Genetics. 2009;182(3):845–50.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.
- 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.View ArticleGoogle Scholar
- Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian Data Analysis, vol. 2. Boca Raton: Chapman & Hall/CRC; 2014.Google Scholar
- 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.View ArticleGoogle Scholar
- 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.