PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of biostsLink to Publisher's site
 
Biostatistics. 2010 January; 11(1): 139–150.
Published online 2009 October 12. doi:  10.1093/biostatistics/kxp043
PMCID: PMC2800164

A hidden Markov random field model for genome-wide association studies

Hongzhe Li*
Department of Biostatistics and Epidemiology, University of Pennsylvania School of Medicine, Philadelphia, PA 19104, USA ; ude.nnepu@ehzgnoh
Zhi Wei
Department of Computer Science, New Jersey Institute of Technology, Newark, NJ 07102, USA

Abstract

Genome-wide association studies (GWAS) are increasingly utilized for identifying novel susceptible genetic variants for complex traits, but there is little consensus on analysis methods for such data. Most commonly used methods include single single nucleotide polymorphism (SNP) analysis or haplotype analysis with Bonferroni correction for multiple comparisons. Since the SNPs in typical GWAS are often in linkage disequilibrium (LD), at least locally, Bonferroni correction of multiple comparisons often leads to conservative error control and therefore lower statistical power. In this paper, we propose a hidden Markov random field model (HMRF) for GWAS analysis based on a weighted LD graph built from the prior LD information among the SNPs and an efficient iterative conditional mode algorithm for estimating the model parameters. This model effectively utilizes the LD information in calculating the posterior probability that an SNP is associated with the disease. These posterior probabilities can then be used to define a false discovery controlling procedure in order to select the disease-associated SNPs. Simulation studies demonstrated the potential gain in power over single SNP analysis. The proposed method is especially effective in identifying SNPs with borderline significance at the single-marker level that nonetheless are in high LD with significant SNPs. In addition, by simultaneously considering the SNPs in LD, the proposed method can also help to reduce the number of false identifications of disease-associated SNPs. We demonstrate the application of the proposed HMRF model using data from a case–control GWAS of neuroblastoma and identify 1 new SNP that is potentially associated with neuroblastoma.

Keywords: Empirical Bayes, False discovery, Iterative conditional model, Linkage disequilibrium

1. INTRODUCTION

Genome-wide association studies (GWAS) are increasingly utilized for identifying novel susceptible genetic variants for complex traits, but there is little consensus on optimal study design and analysis. GWAS are designed to scan the entire genome for the identification of genetic variations associated with phenotypic traits, such as a disease condition, blood pressure, or body mass index. These studies usually examine several hundred thousand single nucleotide polymorphisms (SNPs) that explain most of the genetic variation across the genome as well as SNPs in a large array of candidate gene regions. Most GWAS involve some form of multistage design, which includes an initial scan of hundreds of thousands of SNPs on a sample of cases and controls, followed by testing a subset of the most promising markers on independent samples. Additional markers may also be included at later stages to better characterize the full spectrum of genetic variation in the targeted regions. GWAS have been demonstrated to be a powerful approach for the detection of genetic variants related to complex traits, such as age-related macular degenerative diseases (Klein and others, 2005), prostate and breast cancers (Hunter and others, 2007), and type 2 diabetes (Scott and others, 2007). The Welcome Trust Case Control Consortium (2007) has recently published a GWAS of 7 diseases using 14000 cases and 3000 shared controls . The success of these studies has provided solid evidence that GWAS represent a powerful approach to the identification of genes involved in common human diseases.

To date, the analytical methods of GWAS have largely been limited to the single SNP or SNP–SNP pair analysis, coupled with statistical techniques such as the Bonferroni procedure for controlling multiple comparisons. However, since SNPs in typical GWAS are in linkage disequilibrium(LD) locally, simple Bonferroni correction can potentially be conservative and therefore lead to a loss of statistical power. In addition, if multiple SNPs are all in LD with the true disease variants, effectively utilizing the information from multiple SNPs in LD can potentially increase the power of detecting the SNPs associated with the disease. Local LD information has been mainly utilized in haplotype analysis based on sliding windows (Huang and others, 2007). However, such haplotype analysis may lead to loss of power due to high degrees of freedom. In addition, it is often arbitrary to decide how many SNPs to consider in typical haplotype analysis. Browning and Browning (2007) developed an efficient multilocus association test for whole-genome association studies using localized haplotype clustering, effectively using the local LD information, and demonstrated its application to the Wellcome Trust Case Control Consortium data (Browning and Browning, 2008). Eskin (2008) proposed a data-adaptive multiple comparison procedure to incorporate prior information of LD structure and molecular function in GWAS in order to increase power.

In this paper, we propose to formally incorporate LD information among the SNPs derived from the data set itself or from the HapMap data (Altshuler and others, 2005) into identifying the disease-associated SNPs. In particular, we will first build a weighted LD graph based on pairwise LD measures among the SNPs. Such measures can be derived either from the HapMap project (Altshuler and others, 2005) or from the data set itself. We then propose a hidden Markov random field model (HMRF) on such an LD graph in order to compute the posterior probability that an SNP is associated with the disease. The key is that the posterior probability that a given SNP is associated with the disease depends not only on the genotype data observed at this SNP but also on the genotypes of the SNPs that are in strong LD with this particular SNP. In addition, we propose a simple empirical Bayes (EB) method to borrow information across all the markers in estimating the model parameters. We provide an efficient iterative conditional mode (ICM) algorithm (Besag, 1986) to estimate the parameters and Gibbs sampling approach for estimating the posterior probabilities. These posterior probabilities can be used to define a false discovery rate (FDR) controlling procedure in order to select the relevant SNPs. Different from the standard FDR control procedure of Benjamini and Hochberg (1995), this posterior probability–based FDR control accounts for SNP dependency explicitly and is expected to gain power in detecting the disease-associated SNPs. The proposed HMRF model for accounting for the LD dependency of the SNPs is somewhat similar to some recently developed HMRF models for the analysis of microarray gene expression data in order to account for known regulatory network information (Wei and Li, 2007, 2008; Wei and Pan, 2008).

The rest of the paper is organized as follows: in Section 2, we present our proposed HMRF model for SNPs data and the ICM algorithm for estimating the parameters. We also present a Gibbs sampling approach to estimate the posterior probabilities and use these probabilities to define the FDR. In Section 3, we present results from the analysis of a case–control neuroblastoma (NB) data set and simulation results to demonstrate the methods and to compare the power with single-marker analysis. Finally, we present a brief summary and discussion of our methods and results.

2. AN HMRF MODEL FOR GWAS

2.1. Weighted LD graph and a Markov random field model

Suppose we have m cases and n controls that are genotyped over a set of p SNPs. Let S = {1,…, p} denote the SNP index. We want to determine which SNPs in S are associated with the disease. Let Y = (Y1,…, Ys,…, Yp) be the observed genotype data for the p SNPs, where Ys itself is a vector Ys = (ys1,…, ysm; ys(m+1),…, ys(m + n)), where ysi is the observed genotype for the ith individual at the sth SNP. The typical single SNP analysis often ignores the LD among these SNPs. Here we propose to develop an HMRF model to take into account the LD information in identifying the disease-associated SNPs. We first construct a weighted undirected LD graph G based on pairwise LD information derived from the data or from the HapMap project. Specifically, an edge between SNPs s and s′ is drawn with weight

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx1_ht.jpg

if wss ≠ 0, where I(·) is the indicator function, rss′2 is the r2 measurement of LD between SNPs s and s, and τ is a predetermined cutoff value. We use τ = 0.4 in our simulations and analysis. An example of such an LD graph is given in Figure 1 for the 998 SNPs we use for simulation in Section 3.2.

Fig. 1.

Example of an LD graph, which is used for the simulation study, where each node represents an SNP and a link between 2 SNPs indicates that r2 is greater than 0.40 between them.

For a given SNP s, we then define a random indicator variable as

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx2_ht.jpg

For 2 SNPs s and s′ that are linked on the LD graph, that is, if the r2 between these 2 SNPs are greater than τ, we expect that Xs and Xs are dependent. We propose to model such dependency using a simple discrete Markov random field (MRF) model (Besag, 1974, 1986) with the following joint probability function for X = (X1,…, Xp):

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx3_ht.jpg
(2.1)

where γ and β ≥ 0 are the 2 model parameters and β measures dependencies of Xs for SNPs in LD. In this model, the parameter β > 0 encourages the SNPs that are in LD to have similar values of Xs. We assume that the true association state X is a realization of a locally dependent discrete MRF with a specified distribution {p(X)}. This is in contrast to the hidden Markov model where some time or spatial order of the SNPs has to be assumed. The conditional association state for SNP s, given the states of all neighboring SNPs, is

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx4_ht.jpg

where Ns represents the neighbors of the SNP s on the LD graph.

2.2. An EB model for genotype data

In order to relate the latent vector X to the observed genotypes, we further assume that given any particular realization of X, the random variables Y = (Y1,…, Yp) are conditionally independent with the following conditional density:

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx5_ht.jpg
(2.2)

where P(Ys|Xs) is the joint probability of the observed genotypes over m + n individuals at the SNP s given the latent state Xs.

In order to specify f(Ys|Xs), let θs = (θs1,θs2,θs3) be the genotype frequencies at the sth SNP in the case population and ρs = (ρs1,ρs2,ρs3) be the genotype frequencies at the sth SNP in the control population, for genotype values of 0, 1, and 2, respectively. We assume that both of these frequencies across all the SNPs have a Dirichlet prior with parameter α = (α1,α2,α3):

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx6_ht.jpg

The same prior is also assumed for ρs. For SNP s, let ys + = (ys + ,1,ys + ,2,ys + ,3) denote the observed genotype counts data in the m cases and ys = (ys − ,1,ys − ,2,ys − ,3) denote the observed genotype counts data in the n controls. So if SNP s is not associated with the disease, cases should have the same genotype frequencies as the controls. The combined genotype counts data ys0 = ys + + ys are generated from a trinomial distribution with the genotype frequencies of θs = (θs1,θs2,θs3). Thus, given Xs = 0 the probability of the combined genotype count data is

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx7_ht.jpg
(2.3)

On the other hand, if SNP s is associated with the disease, that is, when Xs = 1, cases and controls should have different genotype frequencies, in which case we have

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx8_ht.jpg
(2.4)

Together, the probability models (2.1), (2.3), and (2.4) define an HMRF model, where X is the vector of the hidden states that follows a discrete MRF, (2.1) and (2.3) and (2.4) define the emission probabilities.

2.3. Parameter estimation using the ICM algorithm and an FDR controlling procedure

In order to estimate the model parameters Φ = (γ, β) in the MRF model and α = (α1,α2,α3) in the emission probabilities, we propose to use the following ICM algorithm, which was originally proposed by Besag (1986) for statistical image analysis. The algorithm involves iterative updates of model parameters and latent states and has the following steps:

  1. Obtain an initial estimate An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx13_ht.jpg based on the single-marker trend test using a p-value of 0.0001.
  2. Estimate α with the value of An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx14_ht.jpg by maximizing the probability of the observed data given by (2.2), l(Y|An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx13_ht.jpg).
  3. Estimate Φ with the value of An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx15_ht.jpg by maximizing the following pseudo-likelihood function:
    An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx9_ht.jpg
    The reason for using the pseudo-likelihood function here is that it is computationally difficult to compute the full likelihood function due to an unknown normalizing constant in the MRF model (2.1).
  4. Carry out a single cycle of ICM based on the current An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx13_ht.jpg, An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx14_ht.jpg, and An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx15_ht.jpg to obtain a new An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx13_ht.jpg. Specifically, for s = 1 to p, update Xs based on
    An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx10_ht.jpg
    (2.5)
  5. Go to step 2 until An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx16_ht.jpg

After the convergence of the algorithm, we propose to sample the latent vector X M times using Gibbs sampling based on the conditional probability (2.5). Based on these samples, we can estimate the posterior probability of qs = Pr(Xs = 0|Y). Let q(s) be the sorted values of qs in descending order. For SNP s, the null hypothesis of interest is

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx11_ht.jpg

Based on these posterior probabilities, let An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx17_ht.jpg, then reject all H(s), s = 1,…, k. This posterior probability–based definition of FDR has been widely used in the analysis of microarray gene expression data (Newton and others, 2001) and has been shown to control the FDR at the appropriate level α and to obtain optimal false nondiscovery rates in the setting of hidden Markov models (Sun and Cai, 2009).

Note that if we fix β = 0 in the MRF model (2.1), that is, if we do not consider the LD dependency of the SNPs, the procedure reduces to a simple EB procedure for the analysis of genetic association data in order to borrow information across all the SNPs.

3. RESULTS

3.1. Application to an NB case–control data set

NB is a common and lethal pediatric malignancy, but despite significant effort the genetic events that initiate tumorigenesis were until recently unknown (Mosse and others, 2008). We had hypothesized that NB is a complex disease that results from the interaction of mutant alleles with relatively low-to-moderate effect on tumor initiation. To identify these genetic variants, Maris and others (2008) reported a GWAS of NB where 1032 NB cases and 2043 controls of European descent were genotyped using the Illumina 550K SNP chips and they observed a significant association between NB and the common minor alleles of 3 consecutive SNPs at chromosome band 6p22 and containing the predicted genes FLJ22536 and FLJ44180 (p-value = 1.71×10 − 9 to 7.01×10 − 10; allelic odds ratio, 1.39–1.40) using single SNP trend tests. Homozygosity for the at-risk G allele of the most significantly associated SNP, rs6939340, resulted in an increased likelihood of NB development (odds ratio, 1.97; 95% confidence interval, 1.58–2.45).

To demonstrate our proposed HMRF modeling approach, we reanalyzed the 30 216 SNPs on chromosome 6 that passed the standard quality controls (see Maris and others, 2008, for details). We grouped these SNPs into groups of 1000 SNPs along the genome and fitted our proposed HMRF for each of these groups in order to facilitate the computation. We ran 10 000 Gibbs samples to obtain the posterior probability of an SNP being associated with the disease. Figure 2 shows the results from both single SNP analysis based on the chi-square tests and our proposed HMRF model. We observed that many SNPs have almost zero estimated probability of being associated with NB. We show in Table 1 the results for 9 SNPs with chi-square test p-value less than 10 − 5. In general, the results agree with each other well, where the 3 SNPs in gene FLJ22536 on chromosome 6, rs6939340, rs4712653 and rs9295536, showed both the most significant association and also the highest posterior probabilities (close to 1.00) of being associated with NB. However, 2 other SNPs, rs11759745 and rs9466269, with chi-square test p-values of 3.33×10 − 5 and 8.37×10 − 5 have a posterior probability of being associated with NB of 0.9933 and 0.9881, clearly indicating that both are also associated with NB. These 2 SNPs are in high LD with the 3 SNPs that have a very strong association with NB, with r2 values ranging from 0.42 to 0.69. This demonstrated that the proposed method is effective in identifying the SNPs with borderline significance at the single-marker level that nonetheless are in high LD with significant SNPs. If we chose these top 5 SNPs based on their estimated posterior probabilities, we expected an FDR of 0.0037.

Table 1.

Results of analysis of NB data set for top 9 SNPs with single SNP chi-square test p-value less than 10−5. Also shown are the estimated posterior probability of being associated with NB based on the HMRF model (Post.prob) and the single SNP p-value ...

Fig. 2.

Results of analysis on 30 216 SNPs on chromosome 6 based on 1251 NB cases and 2236 controls. Legend on the left margin (– log of the p-values) and the circles are results from single SNP analysis and the legend on the right margin (posterior probability ...

Another SNP, rs4487594 on the 6q23.2 band of the chromosome 6, has a p-value of 2.73×10 − 6 but with a relatively high posterior probability of 0.74 to be associated with NB. There are in fact no other SNPs that are in high LD with this SNP, and therefore the posterior probability was only determined by the observed genotypes at this SNP. This is in contrast to some other SNPs with similar single SNP p-values but low posterior probabilities; these SNPs all have neighboring SNPs that are in high LD with them but do not show any association with NB. If we also chose the SNP rs4487594 to be NB associated, the estimated FDR became 0.046. However, if we chose the top 10 SNPs, the FDR increased dramatically to 0.29. So at the FDR level of 0.046, we concluded that the SNP rs4487594 is also associated with NB risk. We tested the association between this SNP and NB in a replication set of 401 cases and 1178 controls and obtained a p-value of 0.054. This further indicated the possible association between the SNP rs4487594 and NB. In contrast, none of the 5 SNPs with similar single SNP p-values but with lower posterior posterior probabilities was significant in the replication set (see Table 1).

In order to assess the false-positive rate of the proposed procedure, we performed the same analysis on 100 data sets created by randomly permuting the disease phenotypes. For the 2 regions that include the SNPs with posterior probability of being associated with NB greater than 0.5 in the original data set, only in 1 permuted data set did we observe 3 SNPs with posterior probability greater than 0.5. The posterior probabilities were 0.55, 0.69, and 0.62, respectively. We did not identify any SNPs with posterior probability of being associated with NB greater than 0.50 all other 99 permuted data sets. These results indicate that the false-positive rate of our proposed procedure is very low.

3.2. Simulation studies

To demonstrate the proposed methods, we conducted a simulation study. In order to obtain more realistic LD patterns among the SNPs, our simulation was based on sampling from the 3075 individuals in the case–control study of NB analyzed in Section 3.2. Particularly, we consider a region of 998 SNPs around the chromosome band 6p22 and select 5 SNPs as the true causal SNPs with relative risk of − 1.4 for 2 such SNPs and 1.4 for another 3 SNPs. These disease-associated SNPs have a total of 33 neighboring SNPs with r2 ≥ 0.4. Figure 1 shows the constructed LD graph for 857 SNPs with at least one neighbor based on the data from 3075 samples from the NB GWAS study. Based on the true genotypes of these 3075 individuals, we simulated the disease probability using the following logistic regression model to give a disease rate of 0.10:

An external file that holds a picture, illustration, etc.
Object name is biostskxp043fx12_ht.jpg

where β0 = − 2.25, βk = − log(1.4) for k = 1,…, 2 and βk = log(1.4) for k = 3,4,5. Based on this model, we repeatedly generated the disease status for the 3075 individuals until we obtained 2000 cases and 2000 controls. The 5 true disease variants were then removed from our analysis. In order to assess the performance of selecting the relevant SNPs, the 33 SNPs that are in LD with the true variants with r2 ≥ 0.4 are defined as true positives, and this definition of true positives is used in defining the sensitivities (SEN) and specificities (SPE).

We repeated the simulation 100 times. For each simulation, we can calculate the SENs, SPEs, and FDRs using different cutoff values of the posterior probabilities. Figure 3 shows the receiver operating characteristic (ROC) curves for 3 different analysis methods, including the single SNP analysis, analysis based on an EB method without utilizing the LD information, and the proposed HMRF method using the LD information, averaged over 100 simulations. Due to the fact that the estimated posterior probability of being associated with disease for some SNPs can be 1, the ROC curve for the HMRF model does not start from the origin (0,0). It is clear that the methods without accounting for the LD resulted in lower SEN in identifying the disease-associated SNPs compared to our proposed HMRF method with LD information. However, we observed that the EB approach did not improve SENs over the single SNP analysis, resulting in almost identical ROC curves. Similar results were also observed when SNPs that are in LD with the true variants with r2 ≥ 0.6 or r2 ≥ 0.8 are defined as true positives (results not shown).

Fig. 3.

The ROC curves (SEN versus 1 – SEP) for the proposed HMRF model (MRF), EB method, and the single SNP Cochran–Armitage trend test (Logistic). Simulation results based on 100 replications.

In order to asses the false-positive rate of the proposed procedure, we performed the same analysis on 100 data sets created by randomly permuting the disease phenotypes for each of the 100 replications generated as above. We observed that 1 permuted data set resulted in 1 SNP with posterior probability of being associated with the disease of 0.77, 1 permuted data set resulted in 2 SNPs with posterior probabilities of 0.52 and 0.90, and the rest of the 98 data sets did not identify any SNPs with the posterior probability of being associated with the disease greater than 0.50. Like the permutation analysis results from the real data set, these results further indicate that the false identification rate of our proposed procedure is very low.

4. DISCUSSION

Since many SNPs in GWAS are in LD, it is important to efficiently utilize such LD information in order to increase the power of detecting disease-associated SNPs. In this paper, we introduce an HMRF model to account for LD in analysis of the SNP data from GWAS, where the pairwise LD measurements are used to define a weighted LD graph, and an MRF model based on the constructed LD graph is used to model the dependency among the SNPs. In addition, a simple EB model is proposed in order to borrow information across all the SNPs. Our simulation indicates that the proposed method can lead to an increase in SEN in identifying the SNPs that are associated with the disease. We used the GWAS of NB to demonstrate the computational feasibility of the methods and the advantage that one can gain over the simple commonly used single SNP analysis. The proposed method is especially effective in identifying the SNPs with borderline significance at the single-marker level that nonetheless are in high LD with significant SNPs, as demonstrated by our analysis of the NB data set. In addition, by simultaneously considering the SNPs in LD, the proposed method can also help to reduce the number of false identifications of disease-associated SNPs.

Although Bonferroni correction has been widely applied in GWAS, analytical and simulation studies by Sabatti and others (2003) have shown that the FDR procedure of Benjamini and Hochberg (1995) can effectively control the FDR for the dependent tests encountered in case–control association studies and increase power over more traditional methods. However, the direct application of the FDR procedure can lead to loss of power due to the dependency among tests, although the FDR can still be controlled (Benjamini and Yekutieli, 2001). The most effective way of correcting this problem relies on developing a precise model for the dependency among the SNPs and incorporating it in the definition of an FDR controlling procedure (Sabatti and others, 2003). The HMRF model proposed provides one such model, and we expect some gain in power over the standard FDR controlling procedure ignoring the dependency. For the hidden Markov model, Sun and Cai (2009) proved the optimal power of a posterior probability–based FDR procedure while controlling the FDR. It is easy to show that if the model we used is the true model and the true parameters are known, our proposed posterior probability–based FDR controlling procedure can indeed control the FDR with a minimal false nondiscovery rate and therefore the optimal power. However, it is not clear that such a result holds when the parameters are estimated using the ICM algorithm. This deserves further investigation.

The method in this paper was developed for case–control GWAS data; however, it can be easily extended to continuous quantitative traits where one only needs to redefine the emission probabilities. One such probability can be based on the EB linear models as in the Limma method for microarray gene expression data analysis (Smyth, 2004; Li and others, 2008). The method in this paper can also be extended to incorporate the prior genetic network such as the protein–protein interaction network information into the analysis of GWAS data, where one can create a weighted LD graph for the SNPs within a given gene. The LD graphs are then linked based on the a priori gene– or protein–protein interaction network to form a large SNP network. An HMRF model can then be defined on this combined network that can help to identify the subnetworks of SNPs that might be related to disease risk. This network-based analysis of genomic data has shown some promise in the analysis of microarray gene expression data (Wei and Li, 2007, 2008; Wei and Pan, 2008). We would also expect some potential gain in SEN in identifying the disease-associated genetic variants in the analysis of GWAS data. It would be interesting to investigate more on how much one can gain in the analysis of GWAS data when the prior biological information is effectively utilized in the analysis.

Finally, we comment on the computational complexity and feasibility of the proposed HMRF approach for analyzing the GWAS data with 500 000 to 1 million SNPs. First, we can analyze the SNPs data chromosome-by-chromosome without losing much information since the SNPs on different chromosomes are not expected to be linked on the LD graph (i.e. they are in linkage equilibrium). This will facilitate parallel computing. Second, for a given chromosome, we can analyze all the SNPs simultaneously or split them into segments of 1000 SNPs as we did in our analysis of the NB data set of 30 216 SNPs. The ICM algorithm converges very fast, usually within 10 iterations for our simulated data and the NB data set. Our implementation of the methods and analysis of data sets was done in R, which has memory limitation to consider all the SNPs. We plan to implement the methods in C/C++, which can further facilitate the analysis of GWAS data.

FUNDING

National Institutes of Health (R01-ES009911, R01-CA127334, and R01-CA78454).

Acknowledgments

We thank Tony Cai and Jichun Xie for discussing FDR controls and Mr Edmund Weisberg, MS at Penn Center for Clinical Epidemiology and Biostatistics, for editorial assistance. Conflict of Interest: None declared.

References

  • Altshuler D, Brooks L, Chakravarti A, Collins F, Daly M, Donnelly P, Consortium IH. A haplotype map of the human genome. Nature. 2005;437:1299–1320. [PMC free article] [PubMed]
  • Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B. 1995;57:289–300.
  • Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under independence. Annals of Statistics. 2001;29:1165–1188.
  • Besag J. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society, Series B. 1974;36:192–225.
  • Besag J. On the statistical analysis of dirty pictures. Journal of the Royal Statistical Society, Series B. 1986;48:259–302.
  • Browning BL, Browning SR. Efficient multilocus association testing for whole genome association studies using localized haplotype clustering. Genetic Epidemiology. 2007;31:365–375. [PubMed]
  • Browning BL, Browning SR. Haplotypic analysis of Wellcome Trust Case Control Consortium data. Human Genetics. 2008;123:273–280. [PMC free article] [PubMed]
  • Eskin E. Increasing power in association studies by using linkage disequilibrium structure and molecular function as prior information. Genome Research. 2008;18:653–660. [PubMed]
  • Huang BE, Amos CI, Lin DY. Detecting haplotype effects in genomewide association studies. Genetic Epidemiology. 2007;31:803–812. [PubMed]
  • Hunter DJ, Kraft P, Jacobs KB, Cox DG, Yeager M, Hankinson SE, Wacholder S, Wang Z, Welch R, Hutchinson A. and others. A genome-wide association study identifies alleles in FGFR2 associated with risk of sporadic postmenopausal breast cancer. Nature Genetics. 2007;39:870–874. [PMC free article] [PubMed]
  • Klein RJ, Zeiss C, Chew EY, Tsai JY, Sackler RS, Haynes C, Henning AK, SanGiovanni JP, Mane SM, Mayne ST. and others. Complement factor H polymorphism in age-related macular degeneration. Science. 2005;3085:385–389. [PMC free article] [PubMed]
  • Li C, Wei Z, Li H. UPenn Biostatistics Working Paper. University of Pennsylvania; 2008. A network-constrained empirical Bayes method for analysis of genomic data.
  • Maris JM, Yael PM, Bradfield JP, Hou C, Monni S, Scott RH, Asgharzadeh S, Attiveh EF, Diskin SJ, Laudenslager M. and others. A genome-wide association study identifies a susceptibility locus to clinically aggressive neuroblastoma at 6p22. The New England Journal of Medicine. 2008;358:2585–2593. [PMC free article] [PubMed]
  • Mosse YP, Laudenslager M, Longo L, Cole KA, Wood A, Attiyeh EF, Laquaglia MJ, Sennett R, Lynch JE, Perri P. and others. Identification of ALK as a major familial neuroblastoma predisposition gene. Nature. 2008;455:930–935. [PMC free article] [PubMed]
  • Newton MA, Kendziorski CM, Richmond CS, Blattner FR, Tsui KW. On differential variability of expression ratios: improving statistical inference about gene expression changes from microarray data. Journal of Computational Biology. 2001;8:37–52. [PubMed]
  • Sabatti C, Service S, Freimer N. False discovery rates in linkage and association linkage genome screens for complex disorders. Genetics. 2003;164:829–833. [PubMed]
  • Scott LJ, Mohlke KL, Bonnycastle LL, Willer CJ, Li Y, Duren WL, Erdos MR, Stringham HM, Chines PS, Jackson AU. and others. A genome-wide association study of type 2 diabetes in Finns detects multiple susceptibility variants. Science. 2007;316:1341–1345. [PMC free article] [PubMed]
  • Smyth GK. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology. 2004;3 Article 3. [PubMed]
  • Sun W, Cai T. Large-scale multiple testing under dependency. Journal of the Royal Statistical Society, Series B. 2009;71:393–424.
  • Wei P, Pan W. Incorporating gene networks into statistical tests for genomic data via a spatially correlated mixture model. Bioinformatics. 2008;24:404–411. [PubMed]
  • Wei Z, Li H. A Markov random field model for network-based analysis of genomic data. Bioinformatics. 2007;23:1537–1544. [PubMed]
  • Wei Z, Li H. A hidden spatial-temporal Markov random field model for network-based analysis of time course gene expression data. The Annals of Applied Statistics. 2008;2:408–429.
  • Welcome Trust Case Control Consortium. Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature. 2007;447:661–678. [PMC free article] [PubMed]

Articles from Biostatistics (Oxford, England) are provided here courtesy of Oxford University Press