Mixture SNPs effect on phenotype in genome-wide association studies

Background Recently mixed linear models are used to address the issue of “missing" heritability in traditional Genome-wide association studies (GWAS). The models assume that all single-nucleotide polymorphisms (SNPs) are associated with the phenotypes of interest. However, it is more common that only a small proportion of SNPs have significant effects on the phenotypes, while most SNPs have no or very small effects. To incorporate this feature, we propose an efficient Hierarchical Bayesian Model (HBM) that extends the existing mixed models to enforce automatic selection of significant SNPs. The HBM models the SNP effects using a mixture distribution of a point mass at zero and a normal distribution, where the point mass corresponds to those non-associative SNPs. Results We estimate the HBM using Gibbs sampling. The estimation performance of our method is first demonstrated through two simulation studies. We make the simulation setups realistic by using parameters fitted on the Framingham Heart Study (FHS) data. The simulation studies show that our method can accurately estimate the proportion of SNPs associated with the simulated phenotype and identify these SNPs, as well as adapt to certain model mis-specification than the standard mixed models. In addition, we analyze data from the FHS and the Health and Retirement Study (HRS) to study the association between Body Mass Index (BMI) and SNPs on Chromosome 16, and replicate the identified genetic associations. The analysis of the FHS data identifies 0.3% SNPs on Chromosome 16 that affect BMI, including rs9939609 and rs9939973 on the FTO gene. These two SNPs are in strong linkage disequilibrium with rs1558902 (Rsq =0.901 for rs9939609 and Rsq =0.905 for rs9939973), which has been reported to be linked with obesity in previous GWAS. We then replicate the findings using the HRS data: the analysis finds 0.4% of SNPs associated with BMI on Chromosome 16. Furthermore, around 25% of the genes that are identified to be associated with BMI are common between the two studies. Conclusions The results demonstrate that the HBM and the associated estimation algorithm offer a powerful tool for identifying significant genetic associations with phenotypes of interest, among a large number of SNPs that are common in modern genetics studies.

example, Vissher [1] discussed 54 loci associated with height which only explained 5% heritability; [2] described 32 loci associated with Body Mass Index (BMI) which explained 1.45% of the variance in BMI. More recently, [3] used mixed linear models (MLM) to simultaneously take into account all the SNPs, which is shown to alleviate the missing-heritability issue.
In this study, we extend the work of [3] to identify the subset of SNPs that are significantly associated with the phenotype of interest, instead of assuming all the SNPs are associative, through a Hierarchical Bayesian model http://www.biomedcentral.com/1471-2164/16/3 (HBM). Similar to [3], all SNPs are considered simultaneously to estimate the heritability, instead of one by one as in the traditional GWAS, hence our HBM also helps to capture missing heritability. Different from the authors in [3], we assume that the SNP effects are distributed as the mixture of a point mass at zero, for those non-effective SNPs, 6 and a normal distribution for those associative SNPs.
Our proposed Hierarchical Bayesian model (HBM) can be represented using the following set of Eqs. 1 and (2). Eq. 1 follows the same set-up as [3]: Y is the n × 1 response vector which corresponds to the individuals' phenotype in our study, X is the design matrix for the fixed effects, W is the standardized genotype matrix, and the vector b contains the N SNP random effects, where the jth element is the random effect corresponding to the jth SNP and is assumed to follow the mixture distribution as in (2), depending on the latent indicator I j , j = 1, . . . , N: (2) One key contribution of our HBM is its capability of automatically selecting significant SNPs while simultaneously incorporating all the SNPs. Eq. 2 is the technical reason behind the selection feature, which can be intuitively understood as follows. Imagine that each SNP is coupled with one Bernoulli indicator I j with success probability p, and all the N Bernoulli indicators are independent. The SNPs then fall into two categories, where the first category contains those with I j = 1, which are the 100 × p% associative SNPs with effects following a normal distribution, while the second group includes the remaining SNPs with I j = 0, who have zero effects. The selection of the associative SNPs is achieved through identifying the SNPs with I j = 1, which are chosen to be those with the largest posterior probability of being 1, through the HBM algorithm as described below: The HBM Algorithm Initialize Choose starting values of β (0) , b (0) , 1) .
a on the main diagonal and 0 everywhere else, with k being the dimension of β.
7. Repeat from Step 1 to Step 6 until convergence.
Several Bayesian variable selection algorithms have been proposed through hierarchical modeling, with applications in genomic studies [4]. Considered a variational Bayes algorithm for GWAS. This method approximates http://www.biomedcentral.com/1471-2164/16/3 the joint posterior density of the hierarchical regression model with a factorized form and minimizes the Kullback-Liebler distance between the factorized form and the full posterior distribution. Although this method is fast to compute, the accuracy of prediction depends on how well the factorized form approximates the posterior distribution of the hierarchical model [5]. Developed a Bayesian variable selection regression algorithm to solve the hierarchical model. They adopted several strategies to improve computational performance, for example, they used marginal associations of the SNPs on the traits as the initial screen step for the latent indicator I j in (2) [6]. This indicates that the distribution of the random effect b j is similar to the marginal estimates of the SNP effects on the traits.
In this study, we modify the standard MCMC algorithm based on the stochastic search algorithm proposed by [7]. The algorithm directly samples the parameters from their posterior distributions and obtain the inferences for the parameters. Because the number of SNPs is large, each iteration of the algorithm involves matrix inversion with the dimension being the number of SNPs. To reduce computation time, we modify the algorithm by sampling the random effects b j conditional on the indicator I j . The modified algorithm significantly reduces computation time, especially when the number of SNPs is large and the mixture probability p is small, while is still able to identify the significant predictors accurately. Detailed description of the algorithm will be stated in Section "Methods": Method. We also implement several computing tricks so that the algorithm can be used to estimate models with the number of SNPs in the order of 100,000 (Section "Example 2").
Our HBM is first applied to analyze simulated data sets in Section "Simulation studies" to show that the proposed algorithm is able to identify the SNPs that are significantly associated with the phenotype and correctly estimate the model parameters as well as PVE, which is defined as the proportion of total genetic variance over total phenotypic variance: where σ 2 g is the total genetic variance which equals σ 2 b in (2) times the number of SNPs. The total phenotypic variance is the sum of the genetic variance σ 2 g and the variance of the error terms of in (2), denoted as σ 2 .
We also compare HBM with the Genome-wide Complex Trait Analysis (GCTA) proposed by [3]. The basic concept of GCTA is to fit the effects of all the SNPs as random effects using a mixed linear model (MLM). Note that the MLM is a special case of our HBM when p = 1. It is shown in our studies that if a large number of SNPs have small/noisy effects on the phenotype, the MLM tends to over-estimate the PVE while the HBM is still able to correctly estimate it. We present in Section "Real data set results" two real data applications through the Framingham Heart Study [8] and the Health and Retirement Study [9], where we study the association between the SNPs on Chromosome 16 and the phenotype body mass index (BMI). We are able to identify associative SNPs on the FTO gene which are consistent with earlier findings in the literature and replicate the results in the two studies.

Simulation studies
The performance of the HBM and MLM is illustrated using two simulated examples with the identical simulation settings but different number of random effects. Example 1 (Section "Example 1") considers 10,000 random effects, while Example 2 (Section "Example 2") has 100,000 random effects and is closer to the scale of real GWAS. Each example also consists of two simulation cases: in Case 1 the random effects follow a mixture distribution of a point mass at zero and a normal distribution, while in Case 2, the random effects follow a mixture of two normal distribution with one of the two has a very small variance, trying to mimic scenarios with a large number of small/noisy effects on the phenotype.
For both simulated examples, genotype information of the individuals from the Framingham Heart Study (FHS) is used as input matrix. Detailed description of the FHS data is provided in Section "The Framingham heart study".

Example 1
In this example, we randomly select 10,000 SNPs on Chromosome 16 of the FHS data and use them as the input genotype matrix, W. The trait Y is then simulated according to the following model: where W is the standardized genotype matrix and b is the allelic effect of the SNPs that will be simulated. The residual effect ( ) is generated from a normal distribution with http://www.biomedcentral.com/1471-2164/16/3 a mean of zero and variance of σ 2 . As discussed above, two simulation cases are generated as follows.
• Simulation Case 1: The random effect b follows a mixture distribution of a point mass at zero plus a normal distribution. In this situation, the SNPs are either associated with the phenotype (whose random effects are distributed as a normal distribution) or not associated with the phenotype (whose random effects will be zero); • Simulation Case 2: The random effect b follows a mixture of two normal distributions with one of the two distributions has a very small variance. In practice, many SNPs might have very small/noisy effects on the complex traits [10]; hence, we are simulating those scenarios with letting some of the SNPs have noisy effects on the phenotype that are normally distributed with a very small variance.
For Simulation Case 1, we randomly select 100 × p% of the SNPs as the ones associated with the phenotype (namely, the association SNPs), and draw their random effects b from the distribution N 0, σ 2 b , and treat the remaining SNPs as non-association SNPs with zero effects. We then fix the PVE at the predetermined value, and simulate the residual from the distribution N 0, σ 2 where σ 2 = j Var(b j )(1/PVE − 1). Phenotype y is generated using W , b and according to Eq. 4. For Simulation Case 2, the data set is generated in a similar way as in Case 1, with the only difference being that the random effects for the non-association SNPs are simulated from N 0, σ 2 where σ is a very small number (e.g. σ =0.01) instead of zero. Table 1 shows the estimation results from the simulated data sets using the HBM and MLM along with the true model parameters. The estimated mixture probabilityp and the random effect varianceσ b by the HBM are close to their corresponding true values, 0.01 and 0.1, respectively. This demonstrates the good performance of our estimation method. In both simulation cases, the MLM severely underestimates σ 2 b , as it divides the total genetic variance onto all the SNPs, instead of just the 1% association SNPs (p = 0.01), which results in underestimation of the genetic effects. In addition, in Simulation Case 2, the estimated PVE from the MLM is much larger than the true value while the HBM gives a closer PVE estimate. The reason is that the MLM can not distinguish the "significant" SNP effects versus those "noisy" effects due to its assumption that all random effects follow the same distribution. Therefore,σ 2 g obtained by MLM would include both "significant" and "noisy" effects and thus lead to overestimation of PVE according to (3). We comment that the simulation model in this case is different from the underlying models assumed by our HBM and the MLM of GCTA. As the results indicate, the HBM is rather robust against such model misspecification.

Example 2
This simulation example is used to demonstrate the performance of the HBM algorithm when the number of SNPs is large (i.e. 100,000), in the order of real GWAS. We have to implement several computational optimizing strategies in order to speed up the computation on such a large number of SNPs as well as to efficiently use the computer memory.
First, in each iteration of the HBM algorithm, we need to invert a square matrix with the rank the same as the number of SNPs. Instead of inverting this matrix directly, we employ the Sherman-Morrison-Woodbury formula [12], to change the matrix inversion to one that only has the same rank as the number of observations, which usually is much smaller than the number of SNPs in genetic a Values in parenthesis are standard errors. b Genetic Varianceσ g 2 is defined in the same way in [11] which equals toσ b 2 × N. N is the number of SNPs whose effect b j follows the N (0,σ b 2 ) distribution. c PVE is calculated asσ studies. Secondly, computation using a large number of SNPs is intensive. Analyzing large datasets of SNPs seems to be impractical on uniprocessor machines. Thus, we carry out the analysis in parallel on UNC-CH's multi-core Linux-based cluster computing server. We write scripts to distribute the computation among multiple cores/CPUs and run multiple computing analyses simultaneously. Our study shows that parallel computing can speed up the computation by a factor of 20 on a 10-core computing node on the cluster. It takes 668.5 minutes and 158GB memory to finish the calculation for the simulated data set with 100,000 SNPs. To consider whole genome data with even more SNPs, the amount of memory and computation power of the server will be the main bottleneck. Similar to Example 1 in Section "Example 1", we consider the same two simulation cases. The estimation results are summarized in Table 2. Even for the larger number of SNPs, our HBM still performs well in both cases, while the same drawbacks exist for the MLM.
For Example 2 (with 100,000 SNPs), Table 3 reports the cross table results of the association SNPs identified by HBM against the truth. As one can see, the HBM can correctly detect 76% and 82% of the true association SNPs for the two simulation cases respectively, and more than 99.9% of the true non-association SNPs. This suggests that the HBM works very well at detecting association SNPs with the false positive rate as low as 0.062% and 0.041%, respectively.

The Framingham heart study
We further apply HBM and MLM to data from the Framingham Heart Study [8] to study genetic associations with the body mass index (BMI). The FHS is a communitybased, prospective, longitudinal study following three generations of participants.
Genotyping for FHS participants was performed using the Affymetrix 500K GeneChip array. Genotypes on the Y chromosome are not included in our analysis. A standard quality control filter is applied to the genotype data. Individuals with 5% or more missing genotype data were excluded from analysis. SNPs that are on the X chromosomes and have a call rate ≤ 99% or a minor allele frequency ≤ 0.01 were also eliminated from the analysis. The application of the quality control procedures resulted in 8,738 individuals with 287,525 SNPs from the 500,000 genotype data. Genotype data were converted to minor allele frequencies for the analysis. One individual of a pair is deleted if the genetic relationship is greater than 0.025. Note that the genetic relationship between individual j and individual k is defined as in [3]: where x ij /x ik is the number of copies of the reference allele for the i th SNP of the j th /k th individual and p i is the frequency of the reference allele. After the above preprocessing, there are 1,915 unrelated individuals in the analysis.
Because the total number of SNPs in the FHS data is close to 300,000, computation is limited by the memory of the UNC server if we include all SNPs in the analysis. Therefore, as a proof of concept, the 13,764 SNPs on Chromosome 16 are used in the analysis. Another reason for considering this chromosome is that it contains an enzyme fat mass and obesity associated protein also known as FTO. We would like to see whether the HBM can identify the SNPs that are significantly correlated with BMI on Chromosome 16, especially those SNPs on the FTO gene. We include the first seven Principle  Components (PCs) for BMI as fixed effects in the model to eliminate genotype correlation induced by biological ancestry.
The estimation results are shown in the left panel of Table 4. We see that the estimated mixture probabilityp from HBM is around 0.3%, which indicates only 0.3 percent of the SNPs on Chromosome 16 are associated with BMI.
The top panel of Table 5 lists the 43 SNPs that are identified by HBM as associated with BMI, which are ordered according to their names, along with the corresponding genes if available. Among these identified SNPs, the SNP rs9939609 variant has been found to be associated with obesity risk among children and adolescents of Beijing, China by [13], and with BMI and waist circumference among European-and African-American youth by [14]. The SNP rs9939973 on the FTO gene has also been found to be related with overweight of children in Korean by [15]. These two SNPs are in strong linkage disequilibrium (LD) with SNP rs1558902 (Rsq=0.901 for rs9939609 and Rsq=0.905 for rs9939973), which had been previously reported in a well-known GWAS [2]. The detection results can also be replicated to some extent: the three SNPs highlighted in red and the genes indicated in blue are also detected in the HRS analysis to be reported below in Section "The health and retirement study".
The predicted allele effects on BMI (kg/m 2 per allele) by HBM and MLM are compared in Table 5, which are calculated as the posterior mean of the random effects under each model. The allele effect predicted by HBM is closer to the findings in the previous GWAS. As an example, we compare SNP rs9939973's effect on BMI with rs1558902's, both of which are on the FTO gene and are highly correlated [2], found that the per allele change in   BMI for SNP rs1558902 is 0.39 (kg/m 2 ) based on a total of 249,796 individuals of European ancestry using a GWAS method. It is much closer to the estimate obtained by HBM (0.224 kg/m 2 ), rather than the much-lower estimate given by MLM (0.014 kg/m 2 ). This comparison indicates that the MLM, assuming that every SNP has an effect on the phenotype, underestimates the SNP effects. One can also see from Table 5 that the estimated PVEs are different (5.6% vs. 13.3%). In Section "Simulation studies", we have shown that the MLM tends to overestimate PVE if there exist many SNPs with small/noisy effects on the phenotype, which we think is also the case for the FHS data here. To demonstrate that there exist SNPs with small effects on BMI, we perform the following multi-scale analysis by varying the amount of SNPs on Chromosome 16 to be included in MLM and showing how the corresponding estimated PVE changes. We first regress BMI on every single SNP and obtain the corresponding p-value. Then we consider a range of varying thresholds on the p-values, and only include those SNPs with a p-value below the threshold in the MLM when estimating the PVE. We systematically increase the p-value threshold so that more and more SNPs that are "less" significant will be included. The idea is that as the p-value threshold increases, more SNPs with small effects on BMI will be included when estimating PVE, which will result in higher PVEs. The estimation results are presented in Table 6. The estimated PVE decreases from 18% to 1% as a decreasing number of SNPs with smaller p-values (below the thresholds from 10 −1 to 10 −7 ) are included in the analysis. The results indicate that when estimating PVE using MLM, the more SNPs with small effects on BMI are included, the higher the estimated PVE is.
In summary, the analysis of the Framingham data reveals several important empirical findings: (1) Among all the SNPs on Chromosome 16, only 0.3% of them are significantly associated with BMI according to HBM; (2) Several association SNPs identified by HBM have also been reported to be significantly related with BMI in previous studies; (3) The MLM tends to underestimate the allele effect on the phenotype while the HBM estimates much closer to previous GWA study results; (4) Because the MLM includes SNPs with small effects on BMI, the estimated PVE by MLM is much higher than the estimate from HBM.

The health and retirement study
In this section, we try to replicate the results in Section "The Framingham heart study" using data from the Health and Retirement Study [9]. The HRS is a longitudinal study of Americans over age 50, conducted every two years from 1992 to 2012; it collects information on economic, health, social, and other factors relevant to aging and retirement. DNA samples were collected in 2006 and 2008. Out of the collected samples, 13,129 individuals were put into genotyping production and 12,507 passed the University of Washington Genetics Coordinating Center's standardized quality control process.
The HRS analysis was performed on 12,237 unrelated individuals and the 11,925 SNPs on Chromosome 16 that are common to those SNPs used in the FHS analysis of Section "The Framingham heart study".
The estimation results are shown in the right panel of Table 3 and Table 4. We first note that the HBM estimates of the proportion of association SNPs are very close in the two studies: 0.34% and 0.42% for FHS and HRS respectively. Both data sets identified the same set of six genes for BMI including the well-known FTO gene. These genes account for about 25% of the genes identified in our analysis.
Forty SNPs are identified to be associated with BMI by the HBM using HRS data set, which are listed in the bottom panel of Table 6. Between the two studies, the HBM identifies three common SNPs to be associated with BMI: rs4782578, rs4784621 and rs9939606 (shown in red), as well as a few common genes (shown in blue). Furthermore, SNP rs9940128 identified using the HRS data is also on the FTO gene, and has been found before to be correlated with BMI by [16,17] and [18].

Conclusion
In this paper, we propose a Hierarchical Bayesian Model (HBM) that extends the MLM of [3]. Our model allows SNP effects on phenotypes of interest to follow a mixture distribution of a point mass at zero and a normal distribution. Our approach addresses the challenge of high-dimensionality in GWAS data by incorporating simultaneous selection of genetics variables that are jointly significant in predicting the phenotype. We employ several computing tricks that enable us to analyze a large number of SNPs (in the order of 100,000). http://www.biomedcentral.com/1471-2164/16/3 We demonstrate the applicability of our approach using both simulated and real data. The simulations are first used to show the accuracy and robustness of the estimation algorithm. We then analyze real data from the FHS and the HRS to identify SNPs on Chromosome 16 that are associated with the body mass index (BMI). The identified SNPs are consistent with earlier findings in the literature, and the results can be replicated across the two studies. The results from both the simulations and the real applications suggest that the MLM tends to over-estimate the proportion of total genetic variance over total phenotypic variance, i.e. PVE. The reason is that the MLM assumes that all the SNPs have effect on the phenotype, including those SNPs with small or noisy effects.
Our work offers a flexible framework that can be extended in several aspects. We now offer some discussion regarding potential future work directions. To analyze the whole-genome data, we can follow [19] and [20] to analyze each chromosome separately. We believe that more work is needed to rigorously study how to aggregate the results, and leave that for future work. The current assumption on the mixture distribution, i.e. a point mass at zero plus a normal distribution, may not be flexible enough to capture genetic effects in certain situations. We intend to relax the distributional assumption to a mixture of a point mass at zero plus a nonparametric distribution as in [21]. One challenge is that the computational short cut we used in this study for Gibbs Sampling might not remain effective for more flexible distributions; hence alternative algorithm have to be considered. Another direction of extension is to relax the independence assumption to allow potential dependence among SNPs within LD blocks. One difficulty then is the estimation of (potentially arbitrary) correlation structure among the SNPs. We are experimenting with adapting the principal factor approximation idea of [22] into our current framework.

Methods
All research involving human subjects, human material, and human data in this paper has been performed in accordance with the Declaration of Helsinki, and with approval from the University of North Carolina-Chapel Hill Institutional Review Board.
The statistical setup of our model is closely related to that of [7,23], and [24]. Our estimation algorithm combines the good features of the three methods, and is the fastest to compute, which is crucial for analyzing GWAS data [7]. First proposed a stochastic search algorithm in order to identify the subset of "promising" subsets predictors through multiple regression. The key feature of the study assumes that the slope of each regressor comes from a mixture of two normal distributions with different variances. The set of slopes with the smaller variance can be considered as being equal to 0. By employing a mixture normal distribution, [7] avoided discontinuity of the mixture between point mass and a normal distribution. However, each step of the iteration would involve all the regressors which is time-consuming, especially when the number of the regressors is large. In another closely related study, [23] explicitly considered the situation in which the regressors' slopes are distributed as zero plus a continuous distribution [23], also allowed sign constraints on the continuous part of the distribution. http://www.biomedcentral.com/1471-2164/16/3 Although Geweke's approach incorporates more realistic assumptions compared with [7], one shortcoming is its slow computation, which makes it unrealistic for large scale genetic studies [24], also tackled the problem of gene selection using the Bayesian variable selection framework. Their algorithm is similar to ours in that both use computational shortcuts to derive the posterior distribution of the random effects conditional on the significant ones in each iteration. However, the proportion of the significant random effects p is pre-specified in their research, while we can estimate it in the process. We relax the known p assumption in our analysis by assuming a prior distribution for p and estimate it using its posterior distribution.
Automatic relevance determination (ARD) is a popular Bayesian variable selection approach [25]. ARD assumes that each random regressor slope follows a normal distribution with mean 0 and a (potentially) distinct variance. The hyperparameters, i.e. the variances, are estimated through maximizing the marginal likelihood, and the variables with zero variance estimates are pruned from the model. The flexibility of the hyperparameters and the estimation algorithm make it difficult to apply ARD to GWAS with a large number of SNPs, which is our primary interest. Our HBM can be viewed as a special case of ARD with only two choices for the hyper parameters: 0 for those non-associative SNPs and σ 2 b for those associative SNPs, and our model is estimated via Gibbs sampling instead of direct likelihood maximization. Similar to the setup in [7], we use a latent variable I j such that when I j = 1, the random effect of the j th SNP, b j , follows N(0, σ 2 b ) and when I j = 0, b j = 0. In addition, I j follows a Bernoulli distribution with Pr(I j = 1) = p, the mixture probability. We seek to estimate the parameters, p, β, σ 2 b and σ 2 e in (1) and (2), as well as predict the random effects b.
Because of the large number of random effects, which equals the number of SNPs in this study, a faster algorithm is employed in our approach based on (1) and (2). The algorithm first modifies the prior distribution of the random effects b to the following mixture normal distribution: and Pr(I j = 1) = p. (6) When σ is a really small number (e.g. σ =0.01), the above mixture normal distribution is approximately a mixture distribution of a normal distribution plus a point mass at zero. Secondly, rather than drawing from the posterior distribution of all the random effects b as a vector, we modify the algorithm to draw b j component wise conditional on the indicator I j . Specifically, if I j = 1, b j is drawn from the marginal conditional distribution f (b j |I j = 1), and for I j = 0, b j is set to zero in each iteration. Thus in each iteration of the Gibbs sampling, the conditional distribution, f (b j |I j = 1), would only involve the columns of W that correspond to I j = 1. In practice, this algorithm speeds up the computation considerably especially in the case when the random effects b have a high dimension and the true mixture probability p is small. For example, it takes 21,727 and 2,004 seconds respectively using the stochastic search algorithm of [7] and our algorithm on a simulated data set with 5,000 SNPs and 5,000 individuals and the true mixture probability p as 0.1.
The Gibbs sampling algorithm for estimating the HBM is provided in Table 1. After a burnin period of 5,000 iterations, the MCMC samples (t) , t = 5, 000, . . . , 7, 000, are obtained. Statistical inference and prediction can then be made based on the posterior distribution of these parameters.

Availability of supporting data
Our study used phenotype and genotype data from two national studies in the United States: the Framingham Heart Study (FHS) and the Health Retirement Survey (HRS), both were conducted with support from NIH. Our group applied for, received approval for using, and downloaded the data from the National Center for Biotechnology Information Genotypes and Phenotypes Database (NCBI dbGaP). The FHS and HRS data should be available to researchers who follow the same application procedures.