Volume 5 Supplement 9

## Genetic Analysis Workshop 17: Unraveling Human Exome Data

# Application of Bayesian regression with singular value decomposition method in association studies for sequence data

- Soonil Kwon
^{1, 2}, - Xiaofei Yan
^{1, 2}, - Jinrui Cui
^{1, 2}, - Jie Yao
^{1, 2}, - Kai Yang
^{1}, - Donald Tsiang
^{1}, - Xiaohui Li
^{1, 2}, - Jerome I Rotter
^{1, 3}and - Xiuqing Guo
^{1, 2, 3}Email author

**5(Suppl 9)**:S57

**DOI: **10.1186/1753-6561-5-S9-S57

© Kwon et al; licensee BioMed Central Ltd. 2011

**Published: **29 November 2011

## Abstract

Genetic association studies usually involve a large number of single-nucleotide polymorphisms (SNPs) (*k*) and a relative small sample size (*n*), which produces the situation that *k* is much greater than *n*. Because conventional statistical approaches are unable to deal with multiple SNPs simultaneously when *k* is much greater than *n*, single-SNP association studies have been used to identify genes involved in a disease’s pathophysiology, which causes a multiple testing problem. To evaluate the contribution of multiple SNPs simultaneously to disease traits when *k* is much greater than *n*, we developed the Bayesian regression with singular value decomposition (BRSVD) method. The method reduces the dimension of the design matrix from *k* to *n* by applying singular value decomposition to the design matrix. We evaluated the model using a Markov chain Monte Carlo simulation with Gibbs sampler constructed from the posterior densities driven by conjugate prior densities. Permutation was incorporated to generate empirical *p*-values. We applied the BRSVD method to the sequence data provided by Genetic Analysis Workshop 17 and found that the BRSVD method is a practical method that can be used to analyze sequence data in comparison to the single-SNP association test and the penalized regression method.

## Background

Association studies usually involve a large number of single-nucleotide polymorphisms (SNPs) (*k*) and a relatively small number of samples (*n*). To avoid multiple testing problems and to consider the effect of multiple SNPs simultaneously, investigators need statistical models that will test multiple SNPs simultaneously. Because standard statistical methods are unable to analyze multiple SNPs simultaneously when *k* is much greater than *n*, Tibshirani [1] introduced the penalized regression (PR) method as an alternative. The method reduces the size of SNP coefficients by treating the coefficients with little effect as zero. In other words, only those SNPs that significantly improve prediction are kept in the model. A potential drawback of this method is that a SNP with a strong marginal effect might be removed from the model if some other SNPs can explain the effect. A second drawback is that the number of SNPs evaluated in the model is controlled by the chosen penalization parameter. Even though the PR method does evaluate multiple SNPs simultaneously when *k* is much greater than *n*, the maximum number of SNPs that can be evaluated in the model is limited by sample size; that is, the method usually cannot test all SNPs simultaneously in large-scale genetic association studies, such as genome-wide association studies.

To evaluate all SNPs simultaneously in one statistical model, we introduced the Bayesian classification with singular value decomposition (BCSVD) method [2]. The BCSVD method can be applied to a dichotomous response variable when *k* is much greater than *n*. The method achieves a massive dimension reduction by applying singular value decomposition to the design matrix in a binary probit model; it estimates the effect of SNPs through the reduced model. Selection of significant SNPs can be achieved by using the empirical *p*-values obtained from permutation. The BCSVD method handles small sample sizes quite well.

To analyze quantitative traits when k is much greater than n, we further developed the Bayesian regression with singular value decomposition (BRSVD) method. We applied the BRSVD method to the sequence data provided by Genetic Analysis Workshop 17 (GAW17). We show that the BRSVD method is a practical method that can be used to analyze sequence data by comparison to the single-SNP association test and PR methods.

## Methods

### BRSVD method

*y*

_{ n }

_{×1}is a vector of quantitative dependent variables,

*X*

_{ n }

_{×}

_{ k }is the design matrix,

*β*

_{ k }

_{×1}is a vector of parameters to be estimated,

*I*

_{ n }is an

*n*×

*n*identity matrix, and

*σ*

^{2}is an unknown variance; as before,

*k*and

*n*are the number of SNPs and the number of samples, respectively. By applying singular value decomposition (SVD) to the design matrix

*X*′ =

*ADF*′, the model in Eq. (1) with the SVD of

*X*can be written:

*L*=

*FD*and:

*γ*a superfactor vector because it is expressed as a linear combination of the original parameters

*β*. The statistical inference will be held on the superfactor vector instead of on

*β.*From Eq. (2), the likelihood function of

*y*given (

*γ*,

*σ*

^{2}) can be obtained as:

*γ.*Let us choose prior densities for (

*β*|

*σ*

^{2}) and

*σ*

^{2}as:

*IG*is the inverted gamma distribution and (

*β**,

*m*,

*a*,

*b*) are known hyperparameters. Because

*γ*=

*A*′

*β*, the conjugate prior density on

*β*implies the conjugate prior density on

*γ*so that:

*γ*,

*σ*

^{2}) can be expressed as:

*γ*,

*σ*

^{2}) can be obtained by multiplying the likelihood function in Eq. (4) to the prior density in Eq. (9):

*γ*and

*σ*

^{2}can be obtained by integrating Eq. (10) with respect to

*σ*

^{2}and

*γ*, respectively. Given the observed data, the marginal posterior density for

*γ*is a multivariate Student’s

*t*distribution in which each element is a Student’s

*t*distribution with (

*n*+

*a*) degrees of freedom and the marginal density for

*σ*

^{2}is:

*γ*can be estimated through a Markov chain Monte Carlo simulation with Gibbs sampler, which starts with the maximum-likelihood estimate. To transform the superfactor vector (

*γ*) in Eq. (2) back to

*β*, which is our original parameter of interest vector, we use the most general solution form for the linear equation (

*γ*=

*A*′

*β*) and achieve the unique solution for

*β*by choosing the generalized inverse of

*A*′ as

*A*[3]. We use a permutation test to estimate the significance of the SNP effects on the phenotype. Let be the estimate of the

*i*th SNP effect from the raw data, and let be the estimate of the

*i*th SNP effect from the jth shuffled data that were obtained by permuting the quantitative trait (y). Define as the difference between and . Then the test statistic can be defined as:

*H*

_{0}:

*β*

_{ i }= 0), the statistic Λ

_{ i }follows the standard normal distribution when

*J*is large:

### Study sample and association analysis

We used the unrelated individuals data distributed by GAW17, which includes 697 individuals, 24,487 SNPs, and 3 covariates (sex, age, and smoking status). We analyzed the first 10 replicates of phenotypes for quantitative risk factor Q1. We first performed the single-SNP association test using the simple linear regression model option in PLINK [4]. Second, we applied the PR method with L1 penalty introduced by Tibshirani [1] using the R package monomvn [5]. We evaluated SNP association with Q1 within the maximum number of SNPs allowed by the package in each step, which is min(*k*, *n* − intercept). Because the package does not provide *p*-values, we used the same permutation technique as in the BRSVD method to obtain empirical *p*-values. Third, we implemented the BRSVD method. To define significant SNPs for each method, we considered the following statistical models: quantitative risk factor Q1 versus the single SNP and the three covariates for the single-SNP association test; quantitative risk factor Q1 versus the maximum number of SNPs allowed by the package plus the three covariates for the PR method; and quantitative risk factor Q1 versus all SNPs (24,487) and the three covariates for the BRSVD method. All SNPs identified as significant for each model were compared to the 39 SNPs listed in the answer sheet distributed by GAW17. The analyses were run for each of the first 10 replicates, and the average of the 10 replicates was summarized (see Results section).

## Results and discussion

### Single-SNP association

*p*-value less than 10

^{−5}, which is an approximate value of 0.05 genome-wide level using Bonferroni correction, as the cut point, the single-SNP test identified age and 50 SNPs as the risk factors for Q1 (Figure 1). By comparison with the answer sheet distributed by GAW17, which listed 39 SNPs that were associated with Q1, only 2 SNPs (C13S522 and C13S523) out of the 50 were correctly identified.

### PR method

*p*= 0.05 (i.e., −log

_{10}(

*p*) = 1.3), age was again identified as a risk factor for Q1. In addition, 15 SNPs were also found to be significant. However, only 3 SNPs (C13S523, C13S522, and C4S1884) out of the 15 were on the 39 risk SNPs list.

### BRSVD method

*p*= 0.05 (i.e., −log

_{10}(

*p*) = 1.3) because all SNPs were evaluated in one test. Age and 45 SNPs were found to be significant. Comparing to the 39 SNPs listed in the answer sheet, 9 SNPs (C1S6533, C4S1877, C4S1889, C4S4935, C6S2981, C13S431, C13S522, C13S523, and C13S524) out of the 45 SNPs were correctly identified.

Summary of validation of the three methods

Empirical outcome | |||||||||
---|---|---|---|---|---|---|---|---|---|

Single-SNP association | PR method | BRSVD method | |||||||

Actual outcome | E′ (= 50) | IE′ (= 24,437) | E′ (= 16) | IE′ (= 24,471) | E′ (= 45) | IE′ (= 24,442) | |||

| TP = 2 | FN = 37 | Sen = 0.051 | TP = 3 | FN = 36 | Sen = 0.077 | TP = 9 | FN = 30 | Sen = 0.231 |

IE (= 24,448) | FP = 48 | TN = 24,400 | Spe = 0.998 | FP = 13 | TN = 24,435 | Spe = 0.9995 | FP = 36 | TN = 24,412 | Spe = 0.9985 |

PPV = 0.04 | NPV = 0.9984 | PPV = 0.187 | NPV = 0.9985 | PPV = 0.2 | NPV = 0.9988 |

## Conclusions

We used three different analysis methods (the single-SNP association analysis method implemented in PLINK, as widely used in genome-wide association studies; the PR method; and the BRSVD method) to identify SNPs that significantly influence the quantitative trait Q1 using the unrelated-individuals sample provided by GAW17. Both the PR and BRSVD methods out-performed the single-SNP association analysis method, suggesting that evaluating multiple SNPs simultaneously not only reduced the problems of multiple testing but also provided more power than single-SNP association in genetic association studies. The BRSVD method had a sensitivity almost three times as high as that of the PR method, suggesting that the BRSVD method is more optimal than the PR method. Another advantage of the BRSVD method is that it requires no specification of parameters compared to the PR method, which requires specification of the penalization parameter that controls the number of variables selected. Moreover, the BRSVD method takes much less computing time than the PR method does. For the association analysis of Q1 in the GAW17 unrelated individuals data, the PR methods used about 1.5 times as much run-time as the BRSVD method. With all factors considered, we believe that the BRSVD method is a good choice for large-scale genetic association study for quantitative traits.

## Declarations

### Acknowledgments

This research was partly supported by the Inflammatory Bowel Disease Program Project Grant DK046763, General Clinical Research Center (GCRC) grant RR00425-30, and Diabetes Endocrinology Research Center (DERC) grant NIDDK DK 063459. The Genetic Analysis Workshop is supported by National Institutes of Health grant R01 GM031575.

This article has been published as part of *BMC Proceedings* Volume 5 Supplement 9, 2011: Genetic Analysis Workshop 17. The full contents of the supplement are available online at http://www.biomedcentral.com/1753-6561/5?issue=S9.

## Authors’ Affiliations

## References

- Tibshirani R: Regression shrinkage via the Lasso. J Roy Stat Soc. 1996, 58: 267-288.
- Kwon S, Cui J, Rhodes S, Tsiang D, Rotter J, Guo X: Application of Bayesian classification with singular value decomposition method in genome-wide association studies. BMC Proc. 2009, 3 (suppl 7): S9-10.1186/1753-6561-3-s7-s9.PubMed CentralView ArticlePubMed
- Graybill F: Theory and application of the linear model. 1976, Belmont, CA, Duxbury Press
- Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, Maller J, Sklar P, de Bakker PIW, Daly MJ, et al: PLINK: a tool set for whole-genome association and population-based linkage analysis. Am J Hum Genet. 2007, 81: 559-575. 10.1086/519795.PubMed CentralView ArticlePubMed
- Gramacy RB: Monomvn: estimation for multivariate normal and student-
*t*data with monotone missingness. R package version 1.8-3. 2010, [http://CRAN.R-project.org/package=monomvn]

## Copyright

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.