Simultaneous SNP identification in association studies with missing data

The response is assumed to be continuous, following a normal distribution. The data set has fully observed family covariates for all the observations. Missing values are imputed only among SNPs, although the method can be modified to impute missing values for phenotypes as well. We focus on testing the relationship between the response and the SNPs. We assume only additive effects among SNPs, although the method can be adapted to quantifying additive and dominance effects of SNPs.

We begin with the linear mixed model

where $\mathbf{Y}_{n\times 1}$ is the phenotypic trait, $\mathbf{X}_{n\times p}$ is the design matrix for family covariates, $\bolds{\beta}_{p\times 1}$ are the coefficients for the family effect, $\mathbf{Z}_{n\times s}$ is the design matrix for SNPs (genotypes), $\bolds{\gamma}_{s\times 1}$ are the coefficients of the additive effect for SNPs, and $\bolds{\varepsilon}_{n\times 1}\sim N(0,\sigma^{2}\mathbf{R})$. The matrix $\mathbf{R}$ is the numerator relationship matrix, describing the degree of kinship between different individuals. (Details on the calculation of $\mathbf{R}$ are given in Appendix A.)

For our application here (carbon isotope data), we have $n=1000$ and $s=450$. In another application of BAMD [Quesada et al. (2010)], they used $n=450$ and $s=400$. In both cases the number of covariates, $p$, was less than $6$. With a fully Bayesian implementation, it is possible to adapt BAMD to the $p\gg n$ case.

Each row of the matrix $\mathbf{Z}$, $\mathbf{Z}_{i},i=1,\ldots,n,$ corresponds to the SNP genotype information of one individual, which can be homozygous for either of the two nucleotides $(-1,1)$ or heterozygous $(0)$. Some of this information may be missing, and we write $\mathbf{Z}_{i}=(\mathbf{Z}_{i}^{\mathrm{obs}},\mathbf{Z}_{i}^{\mathrm{miss}})$, where $\mathbf{Z}_{i}^{\mathrm{obs}}$ are the observed genotypes for the $i$th individual, and $\mathbf{Z}_{i}^{\mathrm{miss}}$ are the missing genotypes. Note two aspects of this framework: {longlist}[(2)]

The values of $\mathbf{Z}_{i}^{\mathrm{miss}}$ are not observed. Thus, if $*$ denotes one missing SNP, a possible $\mathbf{Z}_{i}$ is $\mathbf{Z}_{i}=(1,*,0,0,*,*,1)$.

Individuals are likely to have missing data at different SNP loci. So for 2 different individuals, we might have

For a Bayesian model, we want to put prior distributions on the parameters. We put a noninformative uniform prior for $\bolds{\beta}$, which essentially leads us to least squares estimation. For $\bolds{\gamma}$, we use the normal prior $\bolds{\gamma}\sim N(\mathbf{0},\sigma^{2}\phi^{2}\mathbf{I}_{s})$. Here $\phi^{2}$ is a scale parameter for the variance and $\sigma^{2}$ is the variance parameter. For $\sigma^{2}$ and $\phi^{2}$, we use inverted Gamma priors: $\sigma^{2}\sim\operatorname{IG}(a,b)$ and $\phi^{2}\sim\operatorname{IG}(c,d)$, where IG stands for the inverted Gamma distribution, and $a,b,c,$ and $d$ are constants used in the priors. For specified $a,b,c,d$, the resulting posterior distribution is proper [see Hobert and Casella (1996)].

We consider the case of tag SNP markers in the loblolly pine genome with no significant linkage disequilibrium between them [González-Martínez et al. (2006)]. Therefore, noninformative priors are used for the missing SNPs, meaning that missing data have equal probability of any allelic state. As information increases due to higher marker density, or parental information, or allele frequency in the population, missing data imputation could be constrained accordingly.

We assume that missing SNPs in the data set are Missing at Random (MAR). In particular, let the value of the random variable $T$ denote whether $Z$ is observed, with $T=1$ if the value is observed and $T=0$ if it is missing. If $\xi$ is the parameter of the missing mechanism, then, under the model (1), the MAR assumption results in

So the distribution of the missing SNP could depend on the observed SNPs, and the observed phenotypes. Of course, this does not require such a dependence, it only allows for it.

Other assumptions about missing data mechanisms are less common than MAR. The strongest assumption, and the most difficult to justify, is Missing Completely at Random (MCAR). Under this assumption, the missing data distribution is independent of all observed data, the complete cases can be regarded as sub-samples from the population, and statistical inference with regard to the complete cases is totally valid. Under the model (1), the MCAR assumption can be expressed as

MCAR is regarded as unrealistic and, in most cases, it is not satisfied. It is typically not used to model missing data, and we do not use it here.

Conditional on the MAR assumption, we impute the missing SNPs based on the correlation between SNPs within individuals and between individuals, and use the phenotypic trait information to improve the power of imputation.

In this model, the covariance matrix, $\mathbf{R}$, models the covariance between individuals within the same family, and covariance between individuals across families. Phenotypic traits of related individuals are alike because they share some proportion of SNPs, and genotypes of relatives are similar because they share the same alleles passed on from parents. Various methods can be used to calculate the relationship matrix, such as using a co-ancestry matrix, a kinship matrix, etc. The basic idea is to calculate the probability that 2 individuals share SNPs that are identical by descent. Some methods use pairwise calculations and thus do not guarantee a positive definite relationship matrix, which is unsatisfactory when the relationship matrix is used as covariance matrix. We use the recursive calculation method of Henderson (1976), which gives a numerator relationship matrix that quantifies the probability of sharing a SNP from the same ancestry, based on known family pedigree and parent pedigree in the population. So by calculating this relationship matrix we obtain a probability of 0.5 for the case that two siblings are within the same control-pollinated family and therefore share the same copy of a SNP, or a $0.25$ probability if the two siblings only have one parent in common. For the complex pedigree that we analyze here, there are a total of $9$ categories of relatedness.

The model (1) along with the prior specification allows the use of a Gibbs sampler to estimate parameters. We can iteratively sample from the the full conditionals, given by

The SNPs are contained in the $\mathbf{Z}$ matrix, which includes both the observed SNPs and missing SNPs, and we use the Gibbs sampler to impute the missing SNPs. The Gibbs sampler for the missing data simulates the samples of $Z_{i}^{m}$ according to the distribution of each missing SNP conditional on the rest of observed SNPs and sampled missing SNPs. For a particular SNP $Z^{m}_{ij}$, the $j$th missing SNP in the $i$th individual, the conditional distribution given the rest of the vector $Z^{m}_{i(-j)}$ and all other parameters in the model is

The value $c$ is the genotype currently being considered for that missing SNP, and $c_{l}$ represents any one of the possible genotypes for the SNP. Notice there are only $3$ terms in the denominator sum for each SNP and this is a key point why Gibbs sampling is computationally feasible for our situation with many SNPs and many observations. We also note that the EM algorithm, which provides an alternative method of parameter estimation, can require a prohibitive amount of computation. See Appendix B.

In order to compare models, we shall use the Bayes factor comparing each candidate model to the full model, given by

where $p$ denotes the full model. We now can compare models $\delta$ and $\delta^{\prime}$ through their Bayes factors, as a larger Bayes factor corresponds to a model that explains more variability, when compared to the full model. These pairwise comparisons result in a consistent model selector [O’Hagan and Forster (2004)], and have an advantage over BIC, which is overly biased toward smaller models [Casella et al. (2009)].

We now set up a Metropolis–Hastings (MH) search that has target distribution proportional to the the Bayes factor, $\mathrm{BF}_{\delta}$. Given that we are at model $\delta$, we choose a candidate $\delta^{\prime}$ from a random walk (choose one component at random and switch $0\rightarrow 1$ or $1\rightarrow 0$) with probability $a$ and, with probability $1-a$, we do an independent jump. This is a symmetric candidate, and $\delta^{\prime}$ is accepted with probability $\min\{1,\mathrm{BF}_{\delta^{\prime}}/\mathrm{BF}_{\delta}\}$.

Calculating the Bayes factor in (5) requires knowing the $\mathbf{Z}$ matrix, which is not the case with missing data. Thus, to calculate the Bayes factor, we need to use the imputed $\mathbf{Z}$ matrix from the Gibbs sampler. Thus, we run two Markov chains simultaneously: {longlist}[(2)]

A Gibbs sampler on the full model, to impute the missing data in $\mathbf{Z}$ and estimate all parameters.

A Metropolis–Hastings algorithm on $\delta$, in model space, to find the best model. This algorithm uses an estimated Bayes factor based on the current values in the Gibbs chain. The aim is to search for $\delta^{*}$ such that $\delta^{*}=\arg\max_{\delta}\mathrm{BF}_{\delta}$, but since we are not able to compute $\mathrm{BF}_{\delta}$ exactly for any given $\delta$, we estimate it using samples from the Gibbs sampler, which yields a strongly consistent estimator. We then use the estimated Bayes factor as the target in a stochastic search driven by a Metropolis–Hastings algorithm.

A typical method of estimating a quantity such as (5) would be to use bridge sampling [Meng and Wong (1996)]. However, since the numerator and denominator have different dimensions (but the numerator model is always nested in the denominator model), ordinary bridge sampling will not work. A variation [Chen and Shao (1997)] which accounts for this introduces a weight function to handle the dimension difference. We summarize this strategy in the following proposition.

One particular $g$ function is defined as follows. Let $P_{\delta^{c}}:=\mathbf{Z}_{\delta^{c}}(\mathbf{Z}^{\prime}_{\delta^{c}}\mathbf{Z}_{\delta^{c}})^{-1}\mathbf{Z}^{\prime}_{\delta^{c}}$, $\mathbf{C}_{\delta}:=(\mathbf{Y}-\mathbf{X}\bolds{\beta}-\mathbf{Z}_{\delta}\bolds{\gamma}_{\delta})$, and

which leads to the strongly consistent Bayes factor estimator

Details and proofs of the results given here are in Supplemental Information, Section D [Li et al. (2011)].

For data sets with large numbers of SNPs and phenotypes, the slow computation speed of the Gibbs sampler can be a major problem. We have identified two bottlenecks. First, if the number of SNPs is increased, then for each iteration, the number of missing SNPs to be updated will also increase. Second, in the iterations of the Gibbs sampler, the generation of $\gamma$ involves inverting the matrix $\mathbf{Z}^{\prime}R^{-1}\mathbf{Z}+(1/\phi^{2})I$ each time, as the $\mathbf{Z}$ matrix changes at each iteration. We address these in the following sections.

In this subsection we applyBAMD to simulated data in order to assess the procedure’s performance as we increase the percentage of missing data in the $\mathbf{Z}$ matrix. We simulated a data set with six families, 20 observations in each family and 5 SNPs per observation. The five SNPs are independent of each other. The six families are also independent, so that the parents of the six families are not related and individuals across families are independent. On the other hand, the individuals within each family share the same parents; this relationship is captured via the numerator relationship matrix. From this data set, four data sets with different percentages of missing values, 5%, 10%, 15%, and 20%, were randomly derived. The family effects, $\beta$, which were used to simulate the data, are listed in Table 1. The true SNP effects (additive and dominant effects) used to generate the data are listed in Table 2. When simulating, we let the variance parameter $\sigma^{2}=1$. Our proposed methodology was applied to analyze the data without missing values and also to the new data containing missing values.

Note that for this small simulation, we used a parameterization different from the $\{-1,0,1\}$ coding that we use for larger numbers of SNPs. In this example, each SNP effect is represented as $(\gamma_{a},\gamma_{d})$—the additive and dominant effects of the SNP genotypes.

Tables 1 and 2 summarize the parameter estimation capabilities of BAMD for family and SNP effects. All calculations were based on samples obtained after an initial burn-in of 20,000 iterations of BAMD. The results show that when the percentage of missing values is less than 15%, the proposed methodology yields good estimates for the parameters of direct interest. When the percentage of missing values is greater than 15%, we should be wary of interpreting the results. For example, the true dominant effect for SNP 3 is 0, but the estimate is 1.32 when the percentage of missing values is 20%. Note that the estimate in this case is accurate when the percentage of missing values is less than 10%. We believe the discrepancy arises because one category of genotype for SNP 3 has substantially higher probability and it overpowers the other two categories. When the percentage of missing values increases, the dominated genotype category has only a small chance to be well represented and thus may have unreliable estimates.

Our ultimate goal is to identify significant SNPs from the candidate SNPs. Since we believe that imputation is a tool to obtain better estimates of the parameters, we are not particularly interested in recovering the actual imputed values for the missing SNPs. Nonetheless, the simulation results in Table 3 show that when the probability of one genotype for a certain SNP is dominantly high, the imputed SNPs are correctly identified with probability ranging from $0.55$–$0.85$, being correctly imputed more frequently than the other genotypes distributions (see SNPs 3 and 4).

Here we compare our multiple-imputation missing data algorithm with a program called BIMBAM [Servin and Stephens (2007), http://stephenslab.uchicago.edu/software.html],which is a popular program among geneticists for association genetics and variable selection with missing data (using single imputation).

BAMD and BIMBAM both propose a two-stage procedure that involves first finding a set of significant SNPs, and then running these significant SNPs through a variable selection procedure that finds the best subset of the significant SNPs that describes the variation in the phenotype. Hence, in this study, BIMBAM and BAMD are assessed through the SNPs they find in the first stage and through the final model they put forth.

For the evaluation, we simulated data from the model given in equation (1). The dimensions of the model were fixed to be $n=50$, $p=3$, and $s=25$ throughout. The three families comprised 16, 17, and 17 individuals, respectively. In addition, the $\mathbf{X}$ and $\bolds{\beta}$ matrices were fixed. Entries in the $\mathbf{Z}$ matrix took three possible values, mirroring the real-life situation, when they would represent genotypes. Interest lies in discovering the significant coordinates of the $\bolds{\gamma}$ vector (which corresponds to SNP effects), in the presence of missing values in the $\mathbf{Z}$ matrix.

In the simulation study, three factors—percentage of missing values in $\mathbf{Z}$, magnitude of $\bolds{\gamma}$ effects, and the degree of correlation within a family—were varied across different levels. When a particular factor was being investigated, the others were held constant. Here in the main paper, we only present two specific comparisons. The remainder of the results from the simulation study can be found in Supplemental Information [Li et al. (2011)], Section F. In running the study, we simulated several data sets for each case, and observed very consistent results. Hence, in presenting our results, we focus on a single representative data set in each case.

In both of the studies presented here, the $\bolds{\gamma}$ vector was generated from a multivariate normal, and any values less than 3 in absolute value were set to 0. After generating the $\mathbf{Y}$ responses, 20% of the entries in the $\mathbf{Z}$ matrix were set to missing before being passed to BAMD and BIMBAM.

The first comparison measures the performance of the procedures when an equicorrelation structure ($\rho$ was set to be 0.8) exists within each of the three families. The second comparison presented here aims to see if BAMD turns up many false positives. The $\bolds{\gamma}$ vector was generated in the same way as earlier, but only the coordinates with the five largest values were kept. The rest were set to 0. In addition, the individuals were assumed to be uncorrelated, that is, $\mathbf{R}=\mathbf{I}_{n}$.

The results for each comparison are summarized through two diagrams. The first (the upper panels in Figures 1 and 2) display the SNPs that BAMD and BIMBAM found to be significant in the first stage. The lower panels in Figures 1 and 2 display the output from the variable selection procedure.

Each figure shows that BAMD significantly outperforms BIMBAM. In the first example, BIMBAM found only one of the $14$ significant SNPs, while BAMD found six. In the second example, there were five significant SNPs, and BIMBAM only found one again, while BAMD found four of them.

Multiple imputation of missing SNP data is the best way to ensure unbiased parameter estimates, which is an important consideration given that SNP effects tend to be small for complex traits of greatest biological interest, and given that results of association studies typically motivate more detailed and labor-intensive investigations of how and why associations were detected.

We used simulation to compare BAMD and BIMBAM for their detection of “correct” vs. “incorrect” SNPs, and found that BAMD performed better than BIMBAM. In practice, this advantage of BAMD over BIMBAM would likely be greatest when missing SNPs are not in LD with nearby SNPs (or adjacency cannot be determined). This is the case in many species, including loblolly pine, in which LD is low and genomic resources such as high-resolution genomic maps and high-density SNP chips for genome scanning are not as well developed as they are for the human genome. The higher computational intensity required for formal multiple imputation in BAMD is a trade-off, however, this has not restricted its practical use in most data sets. For very large data sets, parallel processing seems a logical next step in further increasing the computational efficiency of BAMD.

Our method can be applied to family-based association populations, populations of unrelated genotypes, or combination populations. It can incorporate prior information if known. [The application of BAMD in Quesada et al. (2010) was to a population of unrelated genotypes, where significant SNPs related to disease resistance were found.]

Although here we assume a continuous response variable, the method can be adapted to discrete phenotypes using a probit link. For example, in a case control study, the response would be either case or control status, and with a probit model we add a latent variable in the Gibbs sampler.

Although BAMD successfully detected the same significant SNPs as were previously detected using the family-based method QTDT [González-Martínez et al. (2008)], as well as an additional significant SNP, the BAMD variable selector indicated that a subset of the significant SNPs was sufficient to explain variation in the phenotype carbon isotope discrimination. This is a useful tool for biologists because a simultaneous solution for SNP effects enables detection of numerous SNPs that collectively explain phenotypes, which in turn enables further biological experiments to investigate their underlying basis.

However, the candidate SNPs found by BAMD and QTDT cannot necessarily be deemed “correct” or “incorrect” without additional biological experiments. As such, little more can be stated about the correctness of SNPs 3, 5, 29, and 31 without validation experiments. In the broader context of association testing, it is relevant to note that the use of QTDT is limited to families, whereas BAMD and BIMBAM can be used to detect associations in families as well as populations of unrelated individuals. The ability to use BAMD and BIMBAM in many different types of populations is appealing.

Genome-wide association studies are, typically, marginal processors of the data. That is, each SNP is assessed individually for its association, so simultaneous assessment, or epistasis, cannot be detected. A model such as (1) is assessing the SNPs simultaneously—that is its strength. But how many SNPs should we expect to be able to handle in one model? Computational issues aside, if the number of SNPs is greatly increased, we are then susceptible to the usual regression pitfalls—multicollinearity being the most prevalent. Thus, we recommend using BAMD on smaller sets of SNPs that have had some preprocessing. Thus far, BAMD has been used successfully on a model with $400$ SNPs [Quesada et al. (2010)], and we have tested it on as many as $800$ SNPs.

The missing data problem is common across all genomics data sets, so there is broad potential utility of this method. The assumption of MAR (missing at random), which is reasonable in these contexts, may bear additional research attention. If there are quality concerns about SNP data, there are some statistical steps forward, as noted by Wilson et al. (2010), such as using indicator variables of missingness as predictors. This approach can even be extended to test if missingness is a heritable trait and, if so, the MAR assumption is invalid. Next generation sequencing platforms may generate sufficient data to enable this assumption to be tested and, if borne out, may motivate placement of priors on SNP calls in certain sequence contexts.

Last, software to run the Gibbs sampler and variable selector is available in the R package BAMD.