A Non-Gaussian, Nonparametric Structure for Gene-Gene and Gene-Environment Interactions in Case-Control Studies Based on Hierarchies of Dirichlet Processes

TL;DR

This paper introduces a novel nonparametric Bayesian model using hierarchies of Dirichlet processes to better capture complex gene-gene and gene-environment interactions in case-control studies, improving analysis of genetic data.

Contribution

It develops a new nonparametric Bayesian framework with an efficient MCMC algorithm for modeling intricate dependence structures in genetic data, surpassing previous models.

Findings

Effective detection of gene-environment interactions

Successful application to simulated datasets

Insights into gene-gene interactions in myocardial infarction

Abstract

It is becoming increasingly clear that complex interactions among genes and environmental factors play crucial roles in triggering complex diseases. Thus, understanding such interactions is vital, which is possible only through statistical models that adequately account for such intricate, albeit unknown, dependence structures. Bhattacharya & Bhattacharya (2016b) attempt such modeling, relating finite mixtures composed of Dirichlet processes that represent unknown number of genetic sub-populations through a hierarchical matrix-normal structure that incorporates gene-gene interactions, and possible mutations, induced by environmental variables. However, the product dependence structure implied by their matrix-normal model seems to be too simple to be appropriate for general complex, realistic situations. In this article, we propose and develop a novel nonparametric Bayesian model for case-control genotype data using hierarchies of Dirichlet processes that offers a more realistic and nonparametric dependence structure between the genes, induced by the environmental variables. In this regard, we propose a novel and highly parallelisable MCMC algorithm that is rendered quite efficient by the combination of modern parallel computing technology, effective Gibbs sampling steps, retrospective sampling and Transformation based Markov Chain Monte Carlo (TMCMC). We use appropriate Bayesian hypothesis testing procedures to detect the roles of genes and environment in case-control studies. We apply our ideas to 5 biologically realistic case-control genotype datasets simulated under distinct set-ups, and obtain encouraging results in each case. We finally apply our ideas to a real, myocardial infarction dataset, and obtain interesting results on gene-gene and gene-environment interaction, while broadly agreeing with the results reported in the literature.

Figures

Equations

$$[{\mathbf{X}}_{ijk}]=\sum _{m=1}^M{\pi }_{mijk}\prod _{r=1}^{L}f\left({\mathbf{x}}_{ijkr}|p_{mijkr}\right),$$

$$f\left({\mathbf{x}}_{ijkr}|p_{mijkr}\right)={\left\{p_{mijkr}\right\}}^{x_{ijkr}^1+x_{ijkr}^2}{\left\{1-p_{mijkr}\right\}}^{2-(x_{ijkr}^1+x_{ijkr}^2)}.$$

$$[z_{ijk}=m]={\pi }_{mijk},$$

$$[{\mathbf{X}}_{ijk}|z_{ijk}]=\prod _{r=1}^{L}f\left({\mathbf{x}}_{ijkr}|p_{z_{ijk}ijkr}\right).$$

$${\mathbf{p}}_{1ijk},{\mathbf{p}}_{2ijk},\dots ,{\mathbf{p}}_{Mijk}$$

$${\mathbf{G}}_{ijk}$$

$$\log ({\alpha }_{G,ik})={\mu }_G+{\mathbf{\beta }}_G^T{\mathbf{E}}_{ik},$$

$${\mathbf{G}}_{0,jk}\stackrel{iid}{\sim }$$

$$\log ({\alpha }_{G_{0,k}})={\mu }_{G_0}+{\mathbf{\beta }}_{G_0}^T{\bar{\mathbf{E}}}_k,$$

$${\bar{\mathbf{E}}}_k=\frac{1}{N_k}\sum _{i=1}^{N_k}{\mathbf{E}}_{ik}.$$

$${\mathbf{H}}_k\stackrel{iid}{\sim }DP\left({\alpha }_H\tilde{\mathbf{H}}\right);k=0,1,$$

$$\log ({\alpha }_H)={\mu }_H+{\mathbf{\beta }}_H^T\overline{\stackrel{ˉ}{\mathbf{E}}},$$

$$\overline{\stackrel{ˉ}{\mathbf{E}}}=\frac{{\bar{\mathbf{E}}}_0+{\bar{\mathbf{E}}}_1}{2}.$$

$$p_{mijkr}\stackrel{iid}{\sim }\text{mbox}Beta\left({\nu }_1,{\nu }_2\right),$$

$$[{\mathbf{p}}_{mijk}|{\mathbf{p}}_{lijk};l<m]=\frac{{\alpha }_{G,ik}}{{\alpha }_{G,ik}+m-1}{\mathbf{G}}_{0,jk}\left({\mathbf{p}}_{mijk}\right)+\frac{1}{{\alpha }_{G,ik}+m-1}\sum _{t=1}^{T_{mijk}}{\tilde{n}}_{tmijk}{\delta }_{{\mathbf{\varphi }}_{tijk}}\left({\mathbf{p}}_{mijk}\right),$$

$$[{\mathbf{\varphi }}_{tijk}|{\mathbf{\varphi }}_{lijk};l<t]=\frac{{\alpha }_{G_0,k}}{{\alpha }_{G_0,k}+t-1}{\mathbf{H}}_k\left({\mathbf{\varphi }}_{tijk}\right)+\frac{1}{{\alpha }_{G_0,k}+t-1}\sum _{l=1}^{R_{tk}}{\bar{n}}_{ltik}{\delta }_{{\mathbf{\xi }}_{lk}}\left({\mathbf{\varphi }}_{tijk}\right),$$

$$[{\mathbf{\xi }}_{sk}|{\mathbf{\xi }}_{lk};l<s]=\frac{{\alpha }_H}{{\alpha }_H+s-1}\tilde{\mathbf{H}}\left({\mathbf{\xi }}_{sk}\right)+\frac{1}{{\alpha }_H+s-1}\sum _{l=1}^{S_{sk}}{\breve{n}}_{lsk}{\delta }_{{\mathbf{\zeta }}_l}\left({\mathbf{\xi }}_{sk}\right),$$

$$H_{01}:h_{0j}=h_{1j};j=1,\dots ,J,$$

$$H_{11}:\text{mbox}not{H}_0,$$

$$h_{0j}(\cdot )$$

$$h_{1j}(\cdot )$$

$$H_{02\ell }:{\beta }_{G,\ell }=0\text{mbox}versus{H}_{12\ell }:{\beta }_{G,\ell }\ne 0,$$

$$H_{03\ell }:{\beta }_{G_0,\ell }=0\text{mbox}versus{H}_{13\ell }:{\beta }_{G_0,\ell }\ne 0,$$

$$H_{04\ell }:{\beta }_{H,\ell }=0\text{mbox}versus{H}_{14\ell }:{\beta }_{H,\ell }\ne 0.$$

$${\bar{p}}_{mijk}=\frac{{\sum }_{r=1}^{L_j}{p}_{mijkr}}{L_j}.$$

$$C(i,j_1,j_2,k)=cov\left(\text{mbox}logit({\bar{p}}_{z_{i{j}_{1}k}i{j}_{1}k}),\text{mbox}logit({\bar{p}}_{z_{i{j}_{2}k}i{j}_{2}k})\right),$$

$$H_{05i{j}_1{j}_{2}k}:C(i,j_1,j_2,k)=0\text{mbox}versus{H}_{15i{j}_1{j}_{2}k}:C(i,j_1,j_2,k)\ne 0.$$

$${\mathbf{P}}_{M{i}_{0}jk=0}=\left\{{\mathbf{p}}_{1i_{0}jk=0},{\mathbf{p}}_{2i_{0}jk=0},\dots ,{\mathbf{p}}_{M{i}_{0}jk=0}\right\}$$

$${\mathbf{P}}_{M{i}_{1}jk=1}=\left\{{\mathbf{p}}_{1i_{1}jk=1},{\mathbf{p}}_{2i_{1}jk=1},\dots ,{\mathbf{p}}_{M{i}_{1}jk=1}\right\},$$

$$\text{mbox}logit\left({\bar{\mathbf{P}}}_{M{i}_{0}jk=0}\right)=\left\{\text{mbox}logit\left({\bar{p}}_{1i_{0}jk=0}\right),\text{mbox}logit\left({\bar{p}}_{2i_{0}jk=0}\right),\dots ,\text{mbox}logit\left({\bar{p}}_{M{i}_{0}jk=0}\right)\right\}$$

$$\text{mbox}logit\left({\bar{\mathbf{P}}}_{M{i}_{1}jk=1}\right)=\left\{\text{mbox}logit\left({\bar{p}}_{1i_{1}jk=1}\right),\text{mbox}logit\left({\bar{p}}_{2i_{1}jk=1}\right),\dots ,\text{mbox}logit\left({\bar{p}}_{M{i}_{1}jk=1}\right)\right\}.$$

Full text

In this article, we propose and develop a novel nonparametric Bayesian model for case-control genotype data using hierarchies of Dirichlet processes that offers a more realistic and nonparametric dependence structure between the genes, induced by the environmental variables. In this regard, we propose a novel and highly parallelisable MCMC algorithm that is rendered quite efficient by the combination of modern parallel computing technology, effective Gibbs sampling steps, retrospective sampling and Transformation based Markov Chain Monte Carlo (TMCMC). We use appropriate Bayesian hypothesis testing procedures to detect the roles of genes and environment in case-control studies.

Applying our ideas to 5 biologically realistic case-control genotype datasets simulated under distinct set-ups, we obtain encouraging results in each case. We finally apply our ideas to a real, myocardial infarction dataset, and obtain interesting results on gene-gene and gene-environment interaction, that broadly agree with the results reported in the literature.

Keywords: Case-control study; Hierarchical Dirichlet process; Gene-gene and gene-environment interaction; Myocardial Infarction; Parallel processing; Transformation based MCMC.

1 Introduction

Thanks to intense research on gene-gene interaction, including genome-wide association studies (GWAS), it has become increasingly clear that gene-gene interaction alone is insufficient for explaining most complex diseases. Investigating environmental factors independently of the genetic factors is not sufficient either – biomedical research points towards the importance of interactions between genes and the environment in explaining complex diseases. Indeed, according to ? (see aso ?), considering only the separate contributions of genes and environment to a disease, ignoring their interactions, will lead to incorrect estimation of the disease proportion (the “population attributable fraction”) that is explained by genes, the environment, and their joint effect. In particular, environmental exposures are expected to influence gene-gene interactions of the individuals. A comprehensive overview of gene-environment interaction with various examples is provided in ?.

Since no simple relationship exists between the genes and environment, it is clear that linear or additive models, as are mostly used so far, are inadequate for modeling gene-environment interactions. Also, the logistic model based approaches, (see for example ?, ? and ?) resting on Fisher’s definition of interaction result in the inclusion of a large number of interaction terms even with a moderate number of genetic and environmental factors. The existing methods construct various measures of main effects and interaction effects from the genotype data. For instance, ?, use kernel-based methods, ?, interpret interactions in terms of Kullback-Leibler divergence, ?, use entropy-based information gain to interpret interactions, ?, use common homozygote, heterozygote and rare homozygote to model interactions, to name only a few. Thus, there is no unique notion of such measures of main and interaction effects, and since the final results can be very much dependent on the measures of interactions employed, it does not seem to facilitate comparison with our method of genuinely modeling interactions in clear, well-established, statistical terms based on Hierarchical Dirichlet Process based nonparametric methodologies. The existing methods, including those cited above, try to model interactions through pairwise SNP-SNP interaction kernels and logistic regression, ignoring genes as functional units and are hence not scalable for higher order interactions. Two stage approaches like the BOOST software are also SNP based, which test interactions only on the SNPs found marginally significant in the first stage, leading to higher rate of false positives. While almost all the classical methods are based on Fisher’s idea of GxG, existing Bayesian techniques like BEAM, EpiBN study interaction by identifying the SNPs that influence the disease risk given particular allele combinations, ignoring the genes as functional units. In a nutshell, unlike our Bayesian approach, none of the existing methods, classical or Bayesian, attempts simultaneous modelling of the uncertainties associated with the genes as the functional units along with the interactions, both at SNP and gene level through unified statistical models. The fact that the genetic data may arise from a stratified population with an unknown number of subpopulations makes the problem all the more demanding. ?, in their attempt to study the genetic association with respect to genetic data arising from multiple potentially-heterogeneous subgroups, assume the number of subgroups to be known in advance. Also, the problem of quantifying the strength of heterogeneity, as acknowledged by ?, remains unanswered due to the above considerations and the need of an appropriate prior. The Bayesian semiparametric model proposed by ? takes care of the above mentioned problems by proposing a model based on Dirichlet Processes (DP) and a hierarchical matrix-normal distribution, which encapsulates the mechanism of dependence among genes under environmental effects with respect to genotype data arising out of a possibly stratified population.

We now elaborate on a possible drawback of the dependence structure induced by the modeling strategy of ?, which motivated us to develop our present work based on Hierarchical Dirichlet Processes.

In their model, the relevant gene-gene covariance matrix for individual $i$ is $\tilde{\sigma}_{ii}\boldsymbol{A}$, where $\boldsymbol{A}$ is the gene-gene interaction matrix common to all the individuals in the absence of environmental variables, and $\tilde{\sigma}_{ii}=\sigma_{ii}+\phi$, with $\sigma_{ii}$ being the $i$-th diagonal element of a symmetric, positive definite matrix not associated with the environmental variable, and $\phi$ is a non-negative parameter, to be interpreted as the effect of the environmental variable $\boldsymbol{E}$ on gene-gene interaction. Note that ? assumed that the covariance matrices for all the individuals are affected in the same way by the environmental variable, which seems to be a limitation of the covariance structure. The enviromental variables may affect the gene-gene interactions of individuals differently depending on the extent and type of their exposure to the environmental factors.

In this article, we introduce a novel Bayesian nonparametric model for gene-gene and gene-environment interactions for case-control genotype data that solves the issues detailed above. Our model represents the individual genotype data as finite mixtures based on Dirichlet processes as before, but instead of the hierarchical matrix normal distribution, we introduce a hierarchy of Dirichlet processes that create appropriate nonparametric dependence among the genes induced by the environment, case-control dependence, and dependence among the individuals. As we show, our modeling strategy satisfies all the desirable properties, bypassing the drawbacks of the matrix-normal based model of ?. Although our hierarchical Dirichlet process (HDP) model has parallels with the HDP introduced by ?, our HDP has one more level of hierarchy compared to ?. Moreover, we develop a novel and highly parallelisable Markov Chain Monte Carlo (MCMC) methodology that combines the efficiencies of modern parallel computing infrastructure, Gibbs steps, retrospective sampling methods, and Transformation based Markov Chain Monte Carlo (TMCMC). For the hypothesis testing procedures, we essentially adopt and extend the ideas provided in ?. Application of our model and methods to five simulation experiments for the validation purpose yielded quite encouraging results, and application to a real myocardial infarction (MI) case-control type dataset yielded results that are broadly in agreement with the results reported in the literature, but provided new and interesting insights into the mechanisms of gene-gene and gene-environment interactions.

The rest of our paper is structured as follows. We introduce our HDP-based Bayesian nonparametric gene-gene and gene-environment interaction model in Section 2, and in Section 3 discuss the relevant dependence structures induced by our model. In Section 4 we extend the Bayesian hypothesis testing procedures proposed in ? to learn about the roles of genes, environmental variables and their interactions in case-control studies, with respect to our current HDP model. In Section 5 we briefly discuss the results of application of our model and methodologies to $5$ biologically realistic simulated data sets, the details of which are provided in the supplement, described below. In Section 6 we analyze the real MI dataset using our ideas, demonstrating quite interesting and insightful outcome. Finally, we summarize our work with concluding remarks in Section 7.

Additional details are provided in the supplement, whose sections and figures have the prefix “S-” when referred to in this paper.

2 A new Bayesian nonparametric model for gene-gene and gene-environment interactions

2.1 Case-control genotype data

For $s=1,2$ denoting the two chromosomes, let $y^{s}_{ijkr}=1$ and $y^{s}_{ijkr}=0$ indicate the presence and absence of the minor allele of the $i$-th individual belonging to the $k$-th group (either control or case), for $k=0,1$, with $k=1$ denoting case; at the $r$-th locus of $j$-th gene, where $i=1,\ldots,N_{k}$; $r=1,\ldots,L_{j}$ and $j=1,\ldots,J$; let $N=N_{0}+N_{1}$. Let $\boldsymbol{E}_{i}$ denote a set of environmental variables associated with the $i$-th individual. In what follows, we model this case-control genotype and the environmental data using our Bayesian semiparametric model, described in the next few sections.

2.2 Mixture models based on Dirichlet processes

Let $\boldsymbol{y}_{ijkr}=(y^{1}_{ijkr},y^{2}_{ijkr})$, and if $L=\max\{L_{1},\ldots,L_{J}\}$, let $\boldsymbol{Y}_{ijk}=(\boldsymbol{y}_{ijk1},\boldsymbol{y}_{ijk2},\ldots,\boldsymbol{y}_{ijkL_{j}})$ and $\tilde{\boldsymbol{Y}}_{ijk}=(\tilde{\boldsymbol{y}}_{ijk,{L_{j}+1}},\ldots,\tilde{\boldsymbol{y}}_{ijkL})$, where $\tilde{\boldsymbol{Y}}_{ijk}$ are unobserved and assumed to be missing. We introduce these unobserved variables to match the number of loci for all the genes, which is required so that the vectors of minor allele frequencies come from the distribution having the same dimension. This “dimension-matching” is required for the theoretical development of our modeling ideas; see (2.5) and (2.6).

We assume that for every triplet $(i,j,k)$, $\boldsymbol{X}_{ijk}=(\mathbf{x}_{ijk1},\ldots,\mathbf{x}_{ijkL)}=(\boldsymbol{Y}_{ijk},\tilde{\boldsymbol{Y}}_{ijk})$ have the mixture distribution

[TABLE]

where $f\left(\cdot|p_{mijkr}\right)$ is the Bernoulli mass function given by

[TABLE]

In the above, $M$ denotes the maximum number of mixture components and $p_{mijkr}$ stands for the minor allele frequency at the $r$-th locus of the $j$-th gene for the $i$-th individual of the $k$-th case/control group. Note that minor allele frequency is the frequency at which the second most common allele occurs in a given population.

Allocation variables $z_{ijk}$, with probability distribution

[TABLE]

for $i=1,\ldots,N_{k}$ and $m=1,\ldots,M$, allow representation of (2.1) as

[TABLE]

Following ?, ?, we set $\pi_{mijk}=1/M$, for $m=1,\ldots,M$, and for all $(j,k)$.

Letting $\boldsymbol{p}_{mijk}=\left(p_{mijk1},p_{mijk2},\ldots,p_{mijkL}\right)$, we next assume that

[TABLE]

where $\mbox{DP}\left(\alpha_{G,ik}\boldsymbol{G}_{0,jk}\right)$ stands for Dirichlet process with expected probability measure $\boldsymbol{G}_{0,jk}$ having precision parameter $\alpha_{G,ik}$, with

[TABLE]

where $\boldsymbol{E}_{ik}$ is a $d$-dimensional vector of continuous environmental variable for the $i$-th individual in the $k$-th group, $\boldsymbol{\beta}_{G}$ is a $d$-dimensional vector of regression coefficients, and $\mu_{G}$ is the intercept term. The model can be easily extended to include categorical environmental variables along with the continuous ones.

We clarify a further property of our $M$-component mixture model (2.1). It follows directly from ? that conditional on $G_{ijk}$, (2.1) has the same form as the traditional infinite-dimensional Dirichlet process mixture model independent of $M$, but unlike the latter case, where the data are $iid$, in our case the data have a joint dependence structure involving $M$.

2.3 **Hierarchical Dirichlet processes to introduce dependence between the genes

and case-control status**

We further assume that for $k=0,1$,

[TABLE]

where

[TABLE]

with

[TABLE]

We postulate the last level of hierarchy as

[TABLE]

where

[TABLE]

with

[TABLE]

We specify the base probability measure $\tilde{\boldsymbol{H}}$ as follows: for $m=1,\ldots,M$, $i=1,\ldots,N_{k}$, $k=0,1$, and $r=1,\ldots,L$,

[TABLE]

under $\tilde{\boldsymbol{H}}$, where $\nu_{1},\nu_{2}>0$.

This completes the specification of a hierarchy of Dirichlet processes to build dependence between the genes and the distributions of genotypes of cases-controls given data. Note that our model consists of one more level of hierarchy of Dirichlet processes than considered in the applications of ?, who introduce hierarchical Dirichlet processes (HDP). Moreover, our likelihood based on Dirichlet processes ensuring at most $M$ mixture components, is significantly different from those considered in the applications of ?, which are based on the traditional DP mixture; see ?, ?, ? for details on the conceptual, computational and asymptotic advantages of our modeling style over the traditional DP mixture.

2.4 The Chinese restaurant analogy

An extended version of the Chinese restaurant metaphor used by ? may be considered to illustrate our model. For $k=0,1$, the set of random probability measures $\left\{\boldsymbol{G}_{0,jk};j=1,\ldots,J\right\}$ can be associated with $J$ restaurants. Letting $\tau_{ijk}$ denote the number of tables at the $j$-th restaurant associated with the $i$-th individual, we denote by $\boldsymbol{\phi}_{lijk}$ the dish being served at table $l$ of the $j$-th restaurant for the $i$-th individual. Note that $\left\{\boldsymbol{\phi}_{lijk};l=1,\ldots,\tau_{ijk};i=1,\ldots,N_{k}\right\}$ is a set of $iid$ realizations from $\boldsymbol{G}_{0,jk}$. Thus, we have different sets of realizations from $\boldsymbol{G}_{0,jk}$ for different individuals $i$.

For $k=0,1$, we also let $\boldsymbol{\Xi}_{R_{k}k}=\left\{\boldsymbol{\xi}_{1k},\ldots,\boldsymbol{\xi}_{R_{k}k}\right\}$ denote a set of $R_{k}$ $iid$ realizations from $\boldsymbol{H}_{k}$. Then it follows that for $l=1,\ldots,\tau_{ijk}$, $i=1,\ldots,N_{k}$, and for $j=1,\ldots,J$, $\boldsymbol{\phi}_{lijk}\in\boldsymbol{\Xi}_{R_{k}k}$. In other words, $\boldsymbol{\Xi}_{R_{k}k}$ is the set of distinct elements in the set $\{\boldsymbol{\phi}_{lijk};l=1,\ldots,\tau_{ijk};i=1,\ldots,N_{k};j=1,\ldots,J\}$, and, from the Chinese restaurant perspective, is the set of global dishes among all the restaurants, given $k$.

Finally, let $\boldsymbol{\zeta}_{S}=\left\{\boldsymbol{\eta}_{1},\ldots,\boldsymbol{\eta}_{S}\right\}$ denote a set of $S$ $iid$ realizations from $\tilde{\boldsymbol{H}}$. Then it follows that $\boldsymbol{\zeta}_{S}$ is the set of distinct elements in $\left\{\boldsymbol{\Xi}_{R_{k}k}:k=0,1\right\}$. In other words, $\boldsymbol{\zeta}_{S}$ is the set of global dishes served in all the restaurants, irrespective of $k=0$ or $k=1$.

3 Discussion of the dependence structure induced by our HDP-based model

3.1 Dependence among individuals

It follows from the discussion in Section 2.4 that $\left\{\boldsymbol{\phi}_{lijk};l=1,\ldots,T_{mijk};i=1,\ldots,N_{k}\right\}\in\left\{\boldsymbol{\xi}_{1k},\ldots,\boldsymbol{\xi}_{R_{mk}k}\right\}$, where $\boldsymbol{\xi}_{1k},\ldots,\boldsymbol{\xi}_{R_{mk}k}\stackrel{{\scriptstyle iid}}{{\sim}}\boldsymbol{H}_{k}$. This shows that $\left\{\boldsymbol{\phi}_{lijk};l=1,\ldots,T_{mijk};i=1,\ldots,N_{k}\right\}$ in (3.1) are shared among the individuals, thus creating dependence among the subjects.

For more precise insights regarding the dependence structure, let us first marginalize over $\boldsymbol{G}_{ijk}$ to obtain the joint distribution of $\boldsymbol{P}_{Mijk}=\left\{\boldsymbol{p}_{1ijk},\ldots,\boldsymbol{p}_{Mijk}\right\}$ using the following Polya urn distributions: given $\boldsymbol{G}_{0,jk}$, $\boldsymbol{p}_{1ijk}\sim\boldsymbol{G}_{0,jk}$, and for $m=2,\ldots,M$,

[TABLE]

where $\sum_{t=1}^{T_{mijk}}\tilde{n}_{tmijk}=m-1$. Here $\tilde{n}_{tmijk}=\#\left\{l<m:\boldsymbol{p}_{lijk}=\boldsymbol{\phi}_{tijk}\right\}$.

Since conditionally on $\boldsymbol{G}_{0,jk}$, the marginal distribution of $\boldsymbol{p}_{mijk}$, for $m=1,\ldots,M$ and $i=1,\ldots,N_{k}$, is $\boldsymbol{G}_{0,jk}$, the marginal is unaffected by the environmental variable, but the joint distribution of $\boldsymbol{P}_{Mijk}$ implied by the Polya urn distributions (3.1) shows that the dependence structure of $\boldsymbol{P}_{Mijk}$ is influenced by the regression on $\boldsymbol{E}_{ik}$ through $\alpha_{G,ik}$. This is a very desirable property of our modeling approach, since, in reality, the population minor allele frequencies for the case-control group are not expected to be affected by environmental variables, although environmental exposure is expected to influence dependence among individuals and gene-gene interactions in individuals. Note that marginal distributions depending upon environmental variables may be envisaged only under mutation, but since it is an extremely rare phenomenon and the type of case control type genotype data we are dealing with is not appropriate for such studies, we do not include mutational effects in our model.

3.2 Dependence among the genes

We now show that the gene-gene interactions of the $i$-th individual are affected by $\boldsymbol{E}_{ik}$, but not the marginal effects of the genes.

Dependence among the genes for the $i$-th individual is induced by $\left\{\boldsymbol{\phi}_{tijk};t=1,\ldots,\tau_{ijk};j=1,\ldots,J\right\}$, where, for $t=1,\ldots,\tau_{ijk}$, $\boldsymbol{\phi}_{tijk}\stackrel{{\scriptstyle iid}}{{\sim}}\boldsymbol{G}_{0,jk}$, with $\boldsymbol{G}_{0,jk}\sim DP\left(\alpha_{G_{0},k}\boldsymbol{H}_{k}\right)$. In fact, marginalizing over $\boldsymbol{G}_{0,jk}$ yields the following Polya urn scheme for $\left\{\boldsymbol{\phi}_{tijk};t=1,\ldots,\tau_{ijk}\right\}$:

[TABLE]

where $\bar{n}_{ltik}=\#\left\{l^{\prime}<t:\boldsymbol{\phi}_{l^{\prime}ijk}=\boldsymbol{\xi}_{lk}\right\}$. Note that $\sum_{l=1}^{R_{tk}}\bar{n}_{ltik}=t-1$.

It is clear from (3.2) that $\left\{\boldsymbol{\phi}_{tijk};j=1,\ldots,J\right\}$ share $\left\{\boldsymbol{\xi}_{lk};l=1,\ldots,R_{k}\right\}$, so that the latter set creates dependence among the genes. Moreover, it is also clear from (3.2) that the dependence structure does not depend directly upon $\boldsymbol{E}_{ik}$, but upon $\bar{\boldsymbol{E}}_{k}$. In other words, the gene-gene dependence structure of any individual is not directly influenced by the corresponding environmental variable. However, the dependence structure is also influenced by $\bar{n}_{ltik}$, which depends upon the $i$-th individual in the $k$-th case-control group through $\tau_{ijk}$, which is directly influenced by $\boldsymbol{E}_{ik}$ through $\alpha_{G,ik}$. Thus, as is desirable, our modeling style induces gene-gene interactions that are specific to the individuals and are influenced by the corresponding environmental variables and the averages of the environmental variables within the case-control groups that the individuals belong to.

It is also interesting to observe that in spite of the individual-specific gene-gene interactions, the marginal distributions of $\boldsymbol{\phi}_{tijk}$ remains $\boldsymbol{G}_{0,jk}$ for the non-marginalized version and $\boldsymbol{H}_{k}$ for the marginalized version characterized by (3.2), signifying that the individual genes are not affected by $\boldsymbol{E}_{ik}$.

3.3 Case-control dependence

Finally, we note that

[TABLE]

where $\breve{n}_{lsk}=\#\left\{l^{\prime}<s:\boldsymbol{\xi}_{l^{\prime}k}=\boldsymbol{\zeta}_{l}\right\}$ and $\sum_{l=1}^{S_{sk}}\breve{n}_{lsk}=s-1$. So, $\left\{\boldsymbol{\xi}_{sk};s=1,\ldots,R_{k};k=0,1\right\}$ share $\left\{\boldsymbol{\zeta}_{l};l=1,\ldots,S\right\}$, creating dependence between case and control status. Dependence between case and control status are likely to be caused by various implicit factors and environmental variables that are not accounted for in the study. These factors and environmental variables may be insignificant individually, but together may exert non-negligible influence on cases and controls.

A schematic diagram of our HDP-based model and the dependence structure is depicted in Figure 3.1.

We remark that in a much simpler set-up, the original HDP proposed in ? has also been used by ? for inferring population admixture, allowing for correlations between loci due to linkage disequilibrium.

In Section S-1 we propose an MCMC procedure for the inferential purpose, and in Section S-2 we provide a parallel algorithm for implementing the MCMC method.

4 Detection of the roles of environment, genes and their interactions with respect to our HDP based model

4.1 Formulation of the tests and interpretation of their results

4.1.1 Bayesian test for the impact of the genes on case-control

To test if genes have any effect on case-control, we formulate as in ? and ?, the following hypotheses:

[TABLE]

versus

[TABLE]

where

[TABLE]

In the above, for $k=0,1$, $i_{k}$ is the index such that $\boldsymbol{P}_{Mi_{k}jk}=\left\{\boldsymbol{p}_{1i_{k}jk},\boldsymbol{p}_{2i_{k}jk},\ldots,\boldsymbol{p}_{Mi_{k}jk}\right\}$ is some measure of central tendency of $\left\{\boldsymbol{P}_{Mijk}=\left\{\boldsymbol{p}_{1ijk},\boldsymbol{p}_{2ijk},\ldots,\boldsymbol{p}_{Mijk}\right\};i=1,\ldots,N_{k}\right\}$. Appropriate measures of central tendency, based on clusterings, is discussed in Section 4.2.1.

4.1.2 Bayesian test for the significance of the environmental variables

To check if the environmental variables are significant, we shall test the following: for $\ell=1,\ldots,d$,

[TABLE]

[TABLE]

and

[TABLE]

4.1.3 Bayesian test for significance of gene-gene interaction

In our HDP based nonparametric model there is no readily available quantification of gene-gene interaction unlike the models of ? and ?. Thus, in order to test for gene-gene interaction, it is necessary to first reasonably define such a measurement.

A measure of gene-gene interaction influenced by environmental variables

For our purpose, we first define

[TABLE]

With the above definition, for subject $i$ belonging to case-control group $k$, we consider the following covariance

[TABLE]

as quantification of subject-wise gene-gene dependence that accounts for population memberships of subject $i$ with respect to genes $j_{1}$ and $j_{2}$, through $z_{ij_{1}k}$ and $z_{ij_{2}k}$, where for any $p\in(0,1)$, $\mbox{logit}(p)=\log\left(\frac{p}{1-p}\right)$. Thus, gene-gene interaction associated with our model is subject-specific.

While implementing our model using our parallelised MCMC methodology, we simulate $C(i,j_{1},j_{2},k)$ at each iteration by generating $\left\{p_{mijkr}:r=1,\ldots,L_{j}\right\}$ as many times as required from the respective full conditionals holding the remaining parameters fixed, and then compute the empirical covariance corresponding to (4.9) using the generated $iid$ samples conditionally on the remaining parameters to approximate (4.9).

Formulation of the Bayesian tests for gene-gene interactions

To test for subject-wise gene-gene interaction, we consider the following tests: for $i=1,\ldots,N_{k}$, $k=0,1$, and for $j_{1},j_{2}\in\{1,\ldots,J\}$,

[TABLE]

4.1.4 Interpretations of the results of the above tests

The cases that can possibly arise and the respective conclusions are the following:

• •

For some appropriate divergence measure $d$ between two distributions, if $\underset{1\leq j\leq J}{\max}~{}d(h_{0j},h_{1j})$, is significantly small with high posterior probability, then $H_{01}$ is to be accepted. If $h_{0j}$ and $h_{1j}$ are not significantly different, then it is plausible to conclude that the $j$-th gene is not marginally significant in the case-control study.

• •

Suppose that $H_{01}$ is accepted (so that genes have no significant role) and that at least one of $\beta_{G,\ell}$ or $\beta_{G_{0},\ell}$ or $\beta_{H,\ell}$ is significant, at least for some $\ell$. This may be interpreted as the environmental variable $\boldsymbol{E}$ having some altering effect on all the genes in a way that doesn’t affect the disease status. If $C(i,j_{1},j_{2},k)$ turns out to be significant, then this would additionally imply that the environmental variable $\boldsymbol{E}$ influences interaction between genes $j_{1}$ and $j_{2}$ for the $i$-th individual, but not in a way that is responsible for the case/control status.

• •

If $H_{01}$ is rejected, indicating that the genes are significant, but none of the $\beta_{G,\ell}$, $\beta_{G_{0},\ell}$, $\beta_{H,\ell}$ or $C(i,j_{1},j_{2},k)$ are significant, then only the genes, not $\boldsymbol{E}$, are responsible for the disease. In that case, one may conclude that the disease is of purely genetic nature.

• •

Suppose that $H_{01}$ is rejected, none of $\beta_{G,\ell}$, $\beta_{G_{0},\ell}$, $\beta_{H,\ell}$ is significant, but $C(i,j_{1},j_{2},k)$ is significant for at least some $i,j_{1},j_{2},k$. Then the environmental variable is not significant, and the case/control status of the individuals associated with significant gene-gene interactions can be attributed to purely genetic causes triggered by gene-gene interactions of the individuals.

• •

Now suppose that $H_{01}$ is rejected, and at least one of $\beta_{G,\ell}$, $\beta_{G_{0},\ell}$, $\beta_{H,\ell}$ is significant, but none of the subject-wise gene-gene interactions is significant. Then the environmental variable $\boldsymbol{E}$ does not significantly affect the interactions to determine the case/control status, and marginal effects of the individual genes are responsible for the case/control status of an individual.

• •

If, on the other hand, $H_{01}$ is rejected, at least one of $\beta_{G,\ell}$, $\beta_{G_{0},\ell}$, $\beta_{H,\ell}$ is significant, and $C(i,j_{1},j_{2},k)$ is significant for at least some $i,j_{1},j_{2},k$, then the environmental variable is significant and is responsible for influencing gene-gene interactions within the individuals with significant $C(i,j_{1},j_{2},k)$, which, in turn, affects the case/control status of the individuals.

4.2 Methodologies for implementing the Bayesian tests

4.2.1 Hypothesis testing based on clustering modes

As in ? and ?, here we exploit the concept of “central” clustering introduced by ?. Briefly, central clustering may be interpreted as a suitable measure of central tendency of a set of clusterings. ? particularly consider the mode(s) of the set of clusterings, and provide methods for appropriately obtaining the mode(s) using a suitable metric that they propose to quantify distances between any two clusterings. Their proposed metric is also computationally inexpensive, which makes the concept based on central clusterings extremely useful in practice.

For $k=0,1$, let $i_{k}$ denote the index of the central clusterings of $\boldsymbol{P}_{Mijk}=\left\{\boldsymbol{p}_{1ijk},\boldsymbol{p}_{2ijk},\ldots,\boldsymbol{p}_{Mijk}\right\}$, $i=1,\ldots,N_{k}$. We then study the divergence between the two clusterings of

[TABLE]

and

[TABLE]

for $j=1,\ldots,J$. A schematic diagram illustrating the idea can be found in ?.

Significantly large divergence between the two clusterings clearly indicates that the $j$-th gene is marginally significant.

4.2.2 Enhancement of clustering metric based inference using Euclidean distance

As argued in ?, significantly large clustering distance between $\boldsymbol{P}_{Mjk=0}$ and $\boldsymbol{P}_{Mjk=1}$ indicates rejection of $H_{0}$, but insignificant clustering distance does not necessarily provide strong evidence in favour of the null. In this regard, ? (see also ?) argue that the Euclidean distance is an appropriate candidate to be tested for significance before arriving at the final conclusion. Briefly, we first compute the averages $\bar{p}_{mijk}=\sum_{r=1}^{L_{j}}p_{m,ijkr}/L_{j}$, then consider their logit transformations $\mbox{logit}\left(\bar{p}_{mijk}\right)=\log\left\{\bar{p}_{mijk}/(1-\bar{p}_{mijk})\right\}$. Then, we compute the Euclidean distance between the vectors

[TABLE]

and

[TABLE]

We denote the Euclidean distance associated with the $j$-th gene by

[TABLE]

and denote $\underset{1\leq j\leq J}{\max}~{}d_{E,j}$ by $d^{*}_{E}$.

4.2.3 Formal Bayesian hypothesis testing procedure integrating the above developments

In our problem, we need to test the following for reasonably small choices of $\varepsilon$’s:

[TABLE]

[TABLE]

for $\ell=1,\ldots,d$,

[TABLE]

[TABLE]

[TABLE]

and, for $i=1,\ldots,N_{k}$, $k=0,1$, $j_{1},j_{2}\in\{1,\ldots,J\}$,

[TABLE]

If $H_{0}$ is rejected in (4.11) or in (4.12), we could also test if the $j$-th gene is influential by testing, for $j=1,\ldots,J$, $H_{0,\hat{d}_{j}}:~{}\hat{d}_{j}<\varepsilon_{\hat{d}_{j}}\hskip 5.69054pt\mbox{versus}\hskip 5.69054ptH_{1,\hat{d}_{j}}:~{}\hat{d}_{j}\geq\varepsilon_{\hat{d}_{j}}$, where $\hat{d}_{j}=\hat{d}\left(\boldsymbol{P}_{Mi_{0}jk=0},\boldsymbol{P}_{Mi_{1}jk=0}\right)$; we could also test $H_{0,d_{E,j}}:~{}d_{E,j}<\varepsilon_{d_{E,j}}\hskip 5.69054pt\mbox{versus}\hskip 5.69054ptH_{1,d_{E,j}}:~{}d_{E,j}\geq\varepsilon_{d_{E,j}}$.

4.2.4 Null model and choice of $\varepsilon$

To obtain the null posterior distribution, we fit our HDP-based Bayesian model to the dataset generated from the HDP-based model where the genes are independent and not influenced by the environmental variable, and where there is no difference between the probabilities associated with case and control. For the null data we chose the same number of genes, the same number of loci for each gene, and the same number of cases and controls as the non-null data. We also choose the same value $M$ as in the non-null model, but set $\beta_{G}=\beta_{G_{0}}=\beta_{H}=0$. To generate the data from the null model, we first simulate, independently for $j=1,\ldots,J$, the set $\{\boldsymbol{p}_{m1j0}:m=1,\ldots,M\}$, using the Polya urn scheme involving $\tilde{\boldsymbol{H}}$ and $\alpha_{H}$, and set $\{\boldsymbol{p}_{m1j1}:m=1,\ldots,M\}=\{\boldsymbol{p}_{m1j0}:m=1,\ldots,M\}$, so that there is no difference between the probabilities associated with case and control, and that the genes are independent. Since the simulation method is independent of the environmental variable, it is clear that the genes are not influenced by the environment. Given the probabilities $\{\boldsymbol{p}_{m1j1}:m=1,\ldots,M\}$ and $\{\boldsymbol{p}_{m1j0}:m=1,\ldots,M\}$, we then simulate the data using our Bernoulli model. To the data thus generated, we fit our full HDP-based Bayesian model, to obtain the null posterior.

As in ? here also we specify $\varepsilon$’s as $F^{-1}\left(0.55\right)$, where $F$ is the distribution function of the relevant benchmark null posterior distribution. Recall that the choice $F^{-1}\left(0.55\right)$, rather than the median, ensures that the correct null hypothesis is accepted under the “$0-1$” loss. Note that, for the median, the posterior probability of the true null is $0.5$, while under the “$0-1$” loss, the true null will be accepted if its posterior probability is greater than $1/2$.

5 Simulation studies

For simulation studies, we first generate realistic biological data for stratified population with known gene-environment interaction from the GENS2 software of ?. To this data, we then apply our model and methodologies in an effort to detect gene-environment interaction effects that are present in the data. We consider simulation studies in $5$ different true model set-ups: (a) presence of gene-gene and gene-environment interaction, (b) absence of genetic or gene-environmental interaction effect, (c) absence of genetic and gene-gene interaction effects but presence of environmental effect, (d) presence of genetic and gene-gene interaction effects but absence of environmental effect, and (e) independent and additive genetic and environmental effects. The details of our simulation experiments are provided in Section S-3 of the supplement. Here we briefly summarize the results of our experiments.

In case (a), we correctly obtained clear significance of the influence of genetic effects. Also, $\beta_{H}$ turned out to be very significant, demonstrating significant overall impact of the environmental variable on tsignificant overall impact of the environmental variable on the genes. Interestingly, as one may expect, there are more instances of significant gene-gene interactions in the case group compared to the control group. The posteriors of the number of sub-populations gave high probabilities to the correct number of sub-populations in all the 5 simulation experiments. Quite importantly, we demonstrate in cases (a), (d) and (e) where the genes are relevant, that our HDP model can detect disease predisposing loci (DPL) with more precision compared to the matrix-normal-inverse-Wishart model for gene-environment interactions proposed in ?. In case (b) using our ideas in conjunction with significance testing in a simple logistic regression framework, we are correctly able to conclude that the genetic or gene-environmental effects are insignificant. As in ?, the right conclusion is arrived at in case (c) by utilizing our ideas in conjunction with the Akaike Information Criterion (AIC) in the context of simple logistic regression. Using our Bayesian testing procedure along with the aid of logistic regression, we have been able to correctly obtain insignificance of the environmental variable and significance of the genes. In this experiment, we found no gene-gene interaction in the control group and found two (marginal) instances of gene-gene interaction among the cases. As regards case (e), we note as in ? that additivity of genetic and environmental effects is not supported even by our current HDP-based Bayesian model. In spite of this, we correctly obtained significance of the environmental variable and precisely obtained the DPLs. But the lack of the additivity criterion in our model seems to have forced gene-environment interactions. ? report similar results, who obtained, after eventually resorting to logistic regression, the AIC-based best model consisting of the additive marginal effects of the first gene and the environmental variable, along with an additive intercept, which is broadly consistent with the data-generating mechanism.

6 Application of our HDP based ideas to a real, case-control dataset on Myocardial Infarction

We now consider application of our model and methods to a case-control dataset on early-onset of myocardial infarction (MI) from MI Gen study, obtained from the dbGaP database

http://www.ncbi.nlm.nih.gov/gap. The same dataset has been analyzed by ? without considering the sex variable as the covariate, and by ?, who incorporate the sex variable in their gene-environment interaction model. Although ? obtained significant genetic and gene-gene interaction effects, their later study after considering sex as the environmental variable, revealed strong effects of the sex variable but no significant gene-gene interaction, although many of the genes turned out to be individually significant. In our current HDP based analysis, we once again obtain strong effects of the sex variable, but in contrast with ?, although we obtain significant genetic effects, none of the genes turned out to be significant individually. Moreover, the subject-wise gene-gene interactions, although of small magnitude, turned out to be significant in some cases, and interestingly (and apparently counter-intuitively) seem to be instrumental in counter-acting the disease rather than provoking it.

6.1 Data description

We recall that the MI Gen data obtained from dbGaP consists of observations on presence/absence of minor alleles at $727478$ SNP markers associated with 22 autosomes and the sex chromosomes of $2967$ cases of early-onset myocardial infarction, $3075$ age and sex matched controls. The average age at the time of MI was 41 years among the male cases and 47 years among the female cases. The data broadly represents a mixture of four sub-populations: Caucasian, Han Chinese, Japanese and Yoruban. Using the Ensembl human genome database (http://www.ensembl.org/) we could categorize $446765$ markers out of $727478$ with respect to $37233$ genes.

As in ? we considered 32 genes covering 1251 loci, for 200 individuals. These 1251 loci include 33 SNPs that are believed to be associated with MI and also those that are believed to be associated with different cardiovascular end points like LDL cholesterol, smoking, blood pressure, body mass, etc. Other than the 33 SNPs, the remaining 1218 SNPs are not known to be associated with the disease. See ? for the details and the relevant references.

Since the four broad sub-populations are not unlikely to admit further genetic sub-divisions, it makes sense to set the maximum number of mixture components, $M$, to a value much larger than 4. As before, we set $M=30$; we also set $\nu_{1}=\nu_{2}=1$, so that $\tilde{\boldsymbol{H}}$ is the uniform distribution on $[0,1]$. As in the simulation experiments, here also the structures $\alpha_{G,ik}=0.1\times\exp\left(100+\mu_{G}+\beta_{G}E_{ik}\right)$, $\alpha_{G_{0},k}=0.1\times\exp\left(100+\mu_{G_{0}}+\beta_{G_{0}}\bar{E}_{k}\right)$ and $\alpha_{H}=0.1\times\exp\left(100+\mu_{H}+\beta_{H}\bar{\bar{E}}\right)$, where $\mu_{G},\mu_{G_{0}},\mu_{H}\stackrel{{\scriptstyle iid}}{{\sim}}U(0,1)$ and $\beta_{G},\beta_{G_{0}},\beta_{H}\stackrel{{\scriptstyle iid}}{{\sim}}U(-1,1)$, ensured adequate number of sub-populations and satisfactory mixing of MCMC. For the null data and model, we follow the same procedure as discussed in Section 4.2.4.

6.2 Remarks on model implementation

Our parallel MCMC algorithm detailed in Section S-2 takes about 7 days to generate 30,000 iterations on our VMware consisting of 1 TB RAM, 60 double-threaded, 64-bit physical cores, each running at 2.5 GHz; 50 such cores were available to us. We discard the first $10,000$ iterations as burn-in, using the subsequent 20,000 iterations for our Bayesian inference. Informal convergence diagnostics such as trace plots, although did not demosntarte excellent mixing properties, did not indicate evidence of non-convergence. Some instances are provided in Section S-3 of the supplement.

6.3 Results of the real data analysis

6.3.1 Effect of the sex variable

We obtain $P(|\beta_{G}|<\varepsilon_{\beta_{G}}|\mbox{Data})\approx 0$, $P(|\beta_{G_{0}}|<\varepsilon_{\beta_{G_{0}}}|\mbox{Data})\approx 0$ and $P(|\beta_{H}|<\varepsilon_{\beta_{H}}|\mbox{Data})\approx 1$. In other words, although $\bar{\bar{E}}$ (here $E$ being the sex variable) is insignificant, both $E_{ik}$ and $\bar{E}_{k}$ are very significant. Thus, in this study, sex seems to play an important role in influencing the genes.

6.3.2 Roles of individual genes

With the clustering metric we obtained $P\left(d^{*}<\epsilon_{1}|\mbox{Data}\right)\approx 0.030$ while that with the Euclidean distance we obtained $P\left(d^{*}_{E}<\epsilon_{2}|\mbox{Data}\right)\approx 0.540$. That is, the maximum of the gene-wise clustering metrics turns out to be significant, while the maximum of the gene-wise Euclidean metrics is seen to be insignificant. The same ambiguity was also obtained by ?. The tests of the marginal genes are expected to shed some light regarding this dilemma. The posterior probabilities of the null hypotheses (of no significant genetic influence) are shown in Figure 6.1.

As is observed, none of the individual genes turned out to be significant, for either the clustering metric or the Euclidean metric. Our result is not much different from that of ? who also note that their marginal probabilities, at least for the clustering metric, are not significantly small to provide strong enough evidences against the nulls.

Now, at least from the clustering metric perspective, it is necessary to explain the issue that all the genes are insignificant individually but still the maximum of the gene-wise clustering metrics is significant. The key to this issue seem to be the correlations between the distances, which are induced by gene-gene interactions. We explain this phenomenon using a bivariate normal example. Let $(X_{1},X_{2})$ have a bivariate normal distribution with means 0, variances 1, and correlation $\rho$.

Figure 6.2 depicts the median of $\max\{X_{1},X_{2}\}$ as a function of $\rho$, which is seen to be increasing as $\rho$ decreases from 1 to -1. On the other hand, the medians of the marginal distributions of $X_{1}$ and $X_{2}$ remain zero, irrespective of the value of $\rho$. Thus, it seems that gene-gene interaction does have some role to play in this study.

6.3.3 Gene-gene interactions

Unlike ?, where there is a single gene-gene correlation structure for all the individuals, our current model has provision for subject-specific gene-gene correlations. Figures 6.3 and 6.4 show the typical gene-gene correlations representative of cases and controls in all males and females respectively. Essentially, the pictures represent the gene-gene correlation patterns for all the subjects. The color intensities correspond to the absolute values of the correlations. Although the correlations are small in all the subjects, the tests of hypotheses reveal some interesting structures. Figures 6.5 and 6.6 represent the all possible interacting patterns found in the study. Panel (a) of Figure 6.5 represents 9 male cases where no gene-gene interaction is significant. Panel (b) shows the genetic interaction pattern in some male cases where $AP006216.10$ and $C6orf106$, interact with all the other genes. Panel (c) shows the results of significance tests of gene-gene interactions for some male cases, for whom only $AP006216.10$ interacts with all the other genes in the study. A representative interaction pattern for the male controls shown in panel (d), indicates that only $C6orf106$ or only $AP006216.10$ interacts with every gene, but in a few subjects both $AP006216.10$ and $C6orf106$ interact with all the genes.

Even for the females, the two genes, $AP006216.10$ and $C6orf106$, play significant role in gene-gene interactions. Indeed, in our data, unlike the 9 male cases, there is no female for whom all gene-gene interactions are insignificant. The relevant representative plots for the females, given by Figure 6.6, shows that for all the female cases, only $AP006216.10$ interacts with the other genes. For the female controls, either only $AP006216.10$ or only $C6orf106$ interacts with the other genes, or both $AP006216.10$ and $C6orf106$ interact significantly with the other genes included in the study.

The messages gained from our analysis seem to be somewhat counter-intuitive but perhaps quite insightful. Our tests indicate that the genes have insignificant marginal effect. Thus, some external or non-genetic factors might have some significant role to play. But for most of the subjects, at least one of the genes $AP006216.10$ and $C6orf106$ interact with every other gene. The subjects, for whom no significant genetic interactions involving $AP006216.10$ and $C6orf106$ were detected, turned out to be male cases, indicating that the lack of genetic interaction in these males failed to get them any preventive measure against MI. On the other hand, the interactions of the genes $AP006216.10$ and $C6orf106$ with all the genes seemed to reduce the risk of the disease for the other subjects. Thus, in this study, the gene-gene interactions seem to have a beneficial effect on the subjects. It also seems that only a small proportion of males are prone to the risk of having no beneficial gene-gene interactions.

Note that our results are broadly consistent with those obtained by ? but are more precise and informative. Indeed, they also noted relatively small impact of the individual genes and small gene-gene correlations. Our current ideas and analyses also support their conclusion that external factors (in particular, sex) are perhaps playing a bigger role in explaining case-control with respect to MI. We recall (see ?) that with respect to the data that we used, the empirical conditional probability of a male given case is about $0.38$, and that of a male given control is about $0.50$, so that females seem to be more at risk, given our data. The inherent coherence of the Bayesian paradigm upholds the sex factor by attaching little importance to the individual genes. However, in contrast with ? who found no interacting genes, here it turns out that the genes $AP006216.10$ and $C6orf106$ in interaction with other genes generally lower the risk of the individuals with respect to MI. Importantly, each of the few males having no such interactions turned out to be a case. This seems to be roughly in accordance with the popular belief that males are more susceptible to MI than females. Our Bayesian model coherently weaves together the prior and the data and brings out this information in spite of the data-driven information that females are more prone to MI than males. We also note that ?, who analyzed the same MI dataset using logistic regression, reached the conclusion that there is no significant gene-gene interaction. Thus, their result completely supports that of ? and are also very much in keeping with our current results.

6.3.4 Posteriors of the number of sub-populations

Figures 6.7 and 6.8 show the posteriors of the number of sub-populations for the same males and females associated with Figures 6.5 and 6.6, respectively. Observe that the posteriors are quite similar, with the mode at $3$ and $4$ components receiving the next highest probability. Thus, the 4 sub-populations, irrespective of sex, are well-supported by our model, showing that these can not be further sub-divided genetically. This is not unexpected, since the roles of the individual genes are not significant in our study. Our result broadly agrees with ? who obtained for different genes, the modes at $5$ components, with $4$ components receiving the next highest posterior mass.

7 Summary and conclusion

In this paper, we have proposed a novel Bayesian nonparametric gene-gene and gene-environment interaction model based on hierarchies of Dirichlet processes. This model is a significant improvement over that of ? in the sense of much clear interpretability and accounting for subject-specific gene-gene interactions. Moreover, the interactions arise as natural by-products of our nonparametric structure based on HDP, and are not based on matrix normal distributions, as in ? and ?, and hence, are more realistic. We propose a novel parallel MCMC algorithm to implement our model, that combines powerful technology with conditionally independent structures inherent within our HDP based model and efficient TMCMC methods. The Bayesian tests of hypotheses that we employ in this paper are appropriately modified versions of those proposed in ?.

Applications of our ideas to biologically realistic datasets generated under $5$ different set-ups characterized by different combinations and structures associated with gene-gene and gene-environment interactions demonstrated encouraging performance of our model and methods. Our analysis of the MI dataset showed strong impact of the sex variable, which is consistent with the results of ?. Our tests showed no effect of the individual genes, which is also in keeping with ? who obtained relatively weak marginal effects. But most interestingly, even though we obtained very weak gene-gene correlations in accordance with ? and ?, our tests on gene-gene interaction showed that two genes, $AP006216.10$ and $C6orf106$, generally interact with all the other genes in a beneficial way so as to fight the disease. Moreover, the only situations where all the gene-gene interactions turned out to be insignificant, were the male cases, showing that the usual belief that males are more prone to heart attack than females may hold some value from this perspective.

Although many standard methods are commonly used in GWAS to identify the genetic and the environmental effects, there are several reasons that point towards the fact that our approach is not comparable with the existing methods. Firstly, our main objective is to provide a very flexible Bayesian nonparametric approach that adequately accounts for various levels of uncertainties, whereas the objective of the existing approaches is to provide computationally feasible algorithms for analyzing very large datasets. Since the goals are different, different performance criteria must be devised to measure the performances of the approaches, and hence, they are not comparable. Another point which renders our Bayesian approach incomparable with the existing approaches such as PLINK, Bayesian variable selection regression (BVSR) model and the Bayesian epistasis association mapping (BEAM) etc, is that both the simulated and the real data sets are associated with multiple sub-populations and unlike our Bayesian approach none of the existing methods coherently accounts for the unknown number of sub-populations. As we argued in our paper, for different sub-populations, the genes and the loci may interact differently, and so it is of utmost importance to account for such uncertainty. Indeed, Bhattacharjee et al. (2010) empirically demon- strate that methods ignoring population sub-structures can incur severe bias leading to large-scale false positives.

So far, due to insufficient computational resources, we are compelled to restrict focus on a relatively small portion of the data. For improving our computing infrastructure, we have already taken the initiative of procuring supercomputing facilities from the Govt. of India, a project led, on behalf of Indian Statistical Institute, by the second author of this paper. With such a facility, we will be able to analyze the entire MI dataset with much ease.

Supplementary Material

S-1 An MCMC method using Gibbs sampling and TMCMC

S-1.1 Full conditionals

S-1.1.1 Full conditional of $\boldsymbol{H}_{k}$

First observe that for $k=0,1$, the full conditional of $\boldsymbol{H}_{k}$ is given by

[TABLE]

where $n_{sk}=\#\{r\in\{1,\ldots,R_{k}\}:\boldsymbol{\xi}_{rk}=\boldsymbol{\eta}_{s}\}$ and $n_{\cdot k}=\sum_{s=1}^{S}n_{sk}$.

S-1.1.2 Full conditional of $\boldsymbol{G}_{0,jk}$

Similarly, the full conditional of $\boldsymbol{G}_{0,jk}$ is given, for $j=1,\ldots,J$ and $k=0,1$, by

[TABLE]

where $n_{ljk}=\#\{(t,i)\in\{1,\ldots,\tau_{ijk}\}\times\{1,\ldots,N_{k}\}:\boldsymbol{\phi}_{tijk}=\boldsymbol{\xi}_{lk}\}$ and $n_{\cdot jk}=\sum_{l=1}^{R_{k}}n_{ljk}$.

The full conditionals of $\boldsymbol{H}_{k}$ and $\boldsymbol{G}_{0,jk}$ given by (S-1.1) and (S-1.2) indicate generating the infinite-dimensional random probability measures using Sethuraman’s characterization of Dirichlet processes (see ?). However, in our case, forming the infinite-dimensional Sethuraman’s construction is not necessary; instead, it will be required to simulate from the random probability measures having distributions (S-1.1) and (S-1.2). Such simulations are possible using the retrospective method (see ?) which avoids dealing with infinitely many objects.

S-1.1.3 Full conditional of $\boldsymbol{p}_{mijk}$

The associated Polya urn distribution of $\boldsymbol{p}_{mijk}$ given $\boldsymbol{P}_{Mijk}\backslash\{\boldsymbol{p}_{mijk}\}$, derived by marginalizing over $\boldsymbol{G}_{ijk}$, is the following:

[TABLE]

where $M_{tijk}=\#\{m^{\prime}\in\{1,\ldots,M\}\backslash\{m\}:\boldsymbol{p}_{m^{\prime}ijk}=\boldsymbol{\phi}_{tijk}\}$ and $\delta_{\boldsymbol{\phi}_{tijk}}(\cdot)$ denotes point mass at $\boldsymbol{\phi}_{tijk}$.

Given $z_{ijk}=m$, on combining the Polya urn distribution with the likelihood $\prod_{r=1}^{L}f(x_{ijkr}|p_{mijkr})$ we obtain the following full conditional of $\boldsymbol{p}_{mijk}$:

[TABLE]

Note that in (S-1.5), $\boldsymbol{G}_{0,jk}$, drawn from (S-1.2), is not available in closed form and only admits the form dictated by Sethuraman’s construction, given, almost surely, by

[TABLE]

where $\tilde{p}_{1}=V_{1}$, $\tilde{p}_{l}=V_{l}\prod_{s<l}(1-V_{s})$, for $l\geq 2$, with $V_{1},V_{2},\ldots\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Beta}\left(\alpha_{G_{0},k}+n_{\cdot jk},1\right)$, and for $l=1,2,\ldots$, $\tilde{\boldsymbol{\xi}}_{ljk}\stackrel{{\scriptstyle iid}}{{\sim}}\frac{\alpha_{G_{0},k}\boldsymbol{H}_{k}+\sum_{l=1}^{R_{k}}n_{ljk}\delta_{\boldsymbol{\xi}_{lk}}}{\alpha_{G_{0},k}+n_{\cdot jk}}$.

In (S-1.5), the posterior proportional to $\prod_{r=1}^{L}f(x_{ijkr}|p_{mijkr})\boldsymbol{G}_{0,jk}\left(\boldsymbol{p}_{mijk}\right)$, which we denote by $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$, is the discrete distribution that puts mass $C_{ijk}\tilde{p}_{t}\prod_{r=1}^{L}f(x_{ijkr}|\tilde{\xi}_{tjkr})$ to the point $\tilde{\boldsymbol{\xi}}_{tjk}$, for $t=1,2,\ldots$, where

[TABLE]

is the normalizing constant. Combining these with (S-1.5) it follows that

[TABLE]

where

[TABLE]

is the normalizing constant of $\left[\boldsymbol{p}_{mijk}|\cdots\right]$.

S-1.2 Retrospective method for simulating from $\left[\boldsymbol{p}_{mijk}|\cdots\right]$

From (S-1.8) it follows that, to draw from $\left[\boldsymbol{p}_{mijk}|\cdots\right]$, it is required to simulate from $[\boldsymbol{G}_{0,jk}\left(\boldsymbol{p}_{mijk}\right)|\boldsymbol{X}_{ijk}]$ with probability proportional to $C^{-1}_{ijk}$. However, since $C_{ijk}$ involves an infinite series, its calculation is infeasible. The same issue also prevents the traditional simulation methods to draw from the discrete distribution $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$. In this case, the retrospective sampling method proposed in Section 3.5 of ? is the appropriate method for our purpose. We first briefly discuss the role of such method in simulating from $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$, and then argue that a by-product of the method can be used to estimate $C_{ijk}$ arbitrarily accurately.

S-1.2.1 Retrospective method to draw from $[\boldsymbol{G}_{0,jk}\left(\boldsymbol{p}_{mijk}\right)|\boldsymbol{X}_{ijk}]$

Note that the retrospective method requires $\prod_{r=1}^{L}f(x_{ijkr}|\phi_{tijkr})$ in our case to be uniformly bounded for all $\boldsymbol{\phi}_{tijk}$, which holds in our case, as $f(x_{ijkr}|\phi_{tjkr})$ represents the Bernoulli distribution, which is bounded above by 1. We briefly describe the method as follows. Let

[TABLE]

and

[TABLE]

Let us also define $\breve{p}_{\ell,a}(K)=\tilde{p}_{a}\prod_{r=1}^{L}f(x_{ijkr}|\tilde{\xi}_{ajkr})/c_{\ell}(K)$ and $\breve{p}_{u,a}(K)=\tilde{p}_{a}\prod_{r=1}^{L}f(x_{ijkr}|\tilde{\xi}_{ajkr})/c_{u}(K)$. To simulate from $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$ we first generate $U\sim\mbox{Uniform}(0,1)$, and given $U$, choose $\tilde{\boldsymbol{\xi}}_{tjk}$ when

[TABLE]

In fact, $K$ needs to be increased and $\tilde{p}_{t}$ and $\prod_{r=1}^{L}f(x_{ijkr}|\tilde{\xi}_{ajkr})$ simulated retrospectively, till (S-1.11) is satisfied for some $t\leq K$.

S-1.2.2 Retrospective method for estimating $C_{ijk}$ arbitrarily accurately

By choosing $K$ to be large enough, the quantities $c_{\ell}(K)$ and $c_{u}(K)$ given by (S-1.9) and (S-1.10), respectively, can be made arbitrarily close. In other words, for any $\epsilon>0$, there exists $K_{0}\geq 1$ such that $|c_{\ell}(K)-c_{u}(K)|<\epsilon$, for $K\geq K_{0}$. Thus, for any such $K\geq K_{0}$, one may approximate $C_{ijk}$ with $\left[c_{\ell}(K)\right]^{-1}$. In practice, it is only required to simulate $\tilde{U}\sim\mbox{Uniform}(0,1)$ and simulate from $[\boldsymbol{G}_{0,jk}\left(\boldsymbol{p}_{mijk}\right)|\boldsymbol{X}_{ijk}]$ if $\tilde{U}\leq\bar{C}C^{-1}_{ijk}$. For sufficiently small $\epsilon$ and for finite number of simulations, it will generally hold that $\tilde{U}\leq\bar{C}C^{-1}_{ijk}$ if and only if $\tilde{U}\leq\bar{C}_{\epsilon}c_{\ell}(K)$, for $K\geq K_{0}$, where

[TABLE]

S-1.2.3 **Retrospective method to simulate from

$\frac{\alpha_{G_{0},k}\boldsymbol{H}_{k}+\sum_{l=1}^{R_{k}}n_{ljk}\delta_{\boldsymbol{\xi}_{lk}}}{\alpha_{G_{0},k}+n_{\cdot jk}}$**

Note that the retrospective simulation method requires simulation of $\tilde{\boldsymbol{\xi}}_{ljk}\stackrel{{\scriptstyle iid}}{{\sim}}\frac{\alpha_{G_{0},k}\boldsymbol{H}_{k}+\sum_{l=1}^{R_{k}}n_{ljk}\delta_{\boldsymbol{\xi}_{lk}}}{\alpha_{G_{0},k}+n_{\cdot jk}}$, for $l=1,2,\ldots$. This requires simulation from $\boldsymbol{H}_{k}$ with probability proportional to $\alpha_{G_{0},k}$. For this, we first simulate $U\sim\mbox{Uniform}(0,1)$. We then simulate a realization from $\boldsymbol{H}_{k}$ after generating $\boldsymbol{H}_{k}$ from the Dirichlet process given by (S-1.1). Note that we do not have to generate the entire random probability measure $\boldsymbol{H}_{k}$ for this; we only need to generate as many realizations $\boldsymbol{\eta}^{*}_{lk}$’s from $\frac{\alpha_{H}\tilde{\boldsymbol{H}}+\sum_{s=1}^{S}n_{sk}\delta_{\boldsymbol{\eta}_{s}}}{\alpha_{H}+n_{\cdot k}}$ and as many $p^{*}_{lk}=V^{*}_{lk}\prod_{s<l}(1-V^{*}_{lk})$; $l=1,2\ldots$, with $p^{*}_{1k}=V^{*}_{1k}$, with $V^{*}_{lk}\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Beta}\left(\alpha_{H}+n_{\cdot k},1\right)$, as required to satisfy $\sum_{l=1}^{t-1}p^{*}_{lk}<U\leq\sum_{l=1}^{t}p^{*}_{lk}$, for some $t\geq 1$ (with $p^{*}_{0}=0$). We then report $\tilde{\boldsymbol{\xi}}_{1jk}=\boldsymbol{\eta}^{*}_{tk}$ with probability proportional to $\alpha_{G_{0},k}$ and $\tilde{\boldsymbol{\xi}}_{1jk}=\boldsymbol{\xi}_{\tilde{l}jk}$ with probability proportional to $n_{\tilde{l}jk}$, for $\tilde{l}\in\{1,\ldots,R_{k}\}$. We repeat this procedure for generating $\boldsymbol{\xi}_{ljk}$; $l\geq 2$, by sequentially augmenting the existing simulations of $\boldsymbol{\eta}^{*}_{lk}$’s and $p^{*}_{lk}$’s with new draws from $\tilde{\boldsymbol{H}}$ and $\mbox{Beta}\left(\alpha_{H}+n_{\cdot k},1\right)$, if needed. Indeed, note that for augmentation of $p^{*}_{lk}$’s, only extra $V^{*}_{lk}$’s need to be generated from $\mbox{Beta}\left(\alpha_{H}+n_{\cdot k},1\right)$.

S-1.3 Updating procedure of $z_{ijk}$ and $\boldsymbol{p}_{mijk}$

The full conditional of $z_{ijk}$ is given by the following:

[TABLE]

for $m=1,\ldots,M$.

In Section S-1.2 we have devised a method of simulating from the full conditional of $\boldsymbol{p}_{mijk}$ given the data and the remaining variables. For our convenience, we re-formulate the full conditional in terms of the dishes $\boldsymbol{\phi}_{tjk}$ and the indicators of the dishes, which we denote by $t_{mijk}$, where $t_{mijk}=t$ if and only if $p_{mijk}=\phi_{tijk}$; $t=1,\ldots,\tau_{ijk}$.

Now let $\tau^{(m)}_{ijk}$ denote the number of elements in $\boldsymbol{P}_{Mijk}\backslash\left\{\boldsymbol{p}_{mijk}\right\}$ that arose from $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$. Also let ${\boldsymbol{\phi}^{m}}^{*}_{tijk}=\left\{{\phi^{m}}^{*}_{tijkr};~{}r=1,\ldots,L\right\};~{}t=1,\ldots,\tau^{(m)}_{ijk}$ denote the parameter vectors arising from $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$. Further, let ${\boldsymbol{\phi}^{m}}^{*}_{tijk}$ occur $M_{mtijk}$ times.

Then we update $t_{mijk}$ using Gibbs steps, where the full conditional distribution of $t_{mijk}$ is given by

[TABLE]

where

[TABLE]

In (S-1.14) and (S-1.15), $n_{1mijkr}$ and $n_{2mijkr}$ denote the number of $``a"$ and $``A"$ alleles, respectively, at the $r$-th locus of the $j$-th gene of the $i$-th individual, associated with the $m$-th mixture component. In other words, $n_{1mijkr}=x^{1}_{ijkr}+x^{2}_{ijkr}$ and $n_{2mijkr}=2-\left(x^{1}_{ijkr}+x^{2}_{ijkr}\right)$.

Let ${n}^{*}_{1tijkr}=\sum_{m:t_{mijk}=t}n_{1mijkr}$ and ${n}^{*}_{2tijkr}=\sum_{m:t_{mijk}=t}n_{2mijkr}$. Then, for $t=1,\ldots,\tau_{ijk}$; $r=1,\ldots,L_{j}$; $j=1,\ldots,J$ and $k=0,1$, update ${\boldsymbol{\phi}}^{*}_{tijk}$ by simulating from its full conditional distribution, given by

[TABLE]

The above simulation from $[{\boldsymbol{\phi}}^{*}_{tijk}|\cdots]$ is to be carried out by the retrospective method discussed in Sections S-1.2.1 and S-1.2.3.

S-1.4 Updating the missing data $\tilde{\boldsymbol{Y}}_{ijk}$

From (2.4) it follows that

[TABLE]

Hence, given the other unknowns, $\tilde{\boldsymbol{Y}}_{ijk}$ can be updated by simply simulating from the Bernoulli distributions given by (S-1.17).

S-1.5 **Updating $\mu_{G}$, $\boldsymbol{\beta}_{G}$, $\mu_{G_{0}}$, $\boldsymbol{\beta}_{G_{0}}$, $\mu_{H}$

and $\boldsymbol{\beta}_{H}$ using TMCMC**

S-1.5.1 Relevant factors for updating $\mu_{G}$ and $\boldsymbol{\beta}_{G}$

Let

[TABLE]

where $[\boldsymbol{p}_{mijk}|\boldsymbol{p}_{lijk};l<m]$ is given by (3.1). Let $\pi_{G}(\mu_{G},\boldsymbol{\beta}_{G})$ denote the prior on $(\mu_{G},\boldsymbol{\beta}_{G})$. Note that $\pi_{G}(\mu_{G},\boldsymbol{\beta}_{G})\mathcal{L}_{G}(\mu_{G},\boldsymbol{\beta}_{G})$ is the product of the only factors in the joint model consisting of $\mu_{G}$ and $\boldsymbol{\beta}_{G}$.

S-1.5.2 Relevant factors for updating $\mu_{G_{0}}$ and $\boldsymbol{\beta}_{G_{0}}$

Now let

[TABLE]

where $[\boldsymbol{\phi}_{tijk}|\boldsymbol{\phi}_{lijk};l<t]$ is given by (3.2).

Let $\pi_{G_{0}}(\mu_{G_{0}},\boldsymbol{\beta}_{G_{0}})$ denote the prior on $(\mu_{G_{0}},\boldsymbol{\beta}_{G_{0}})$. Then $\pi_{G_{0}}(\mu_{G_{0}},\boldsymbol{\beta}_{G_{0}})\mathcal{L}_{G_{0}}(\mu_{G_{0}},\boldsymbol{\beta}_{G_{0}})$ is the functional form associated with $\mu_{G_{0}}$ and $\boldsymbol{\beta}_{G_{0}}$.

S-1.5.3 Relevant factors for updating $\mu_{H}$ and $\boldsymbol{\beta}_{H}$

Finally, we let

[TABLE]

where $[\boldsymbol{\xi}_{sk}|\boldsymbol{\xi}_{lk};l<s]$ is given y (3.3).

Let $\pi_{H}(\mu_{H},\boldsymbol{\beta}_{H})$ be the prior on $(\mu_{H},\boldsymbol{\beta}_{H})$. Then $\pi_{H}(\mu_{H},\boldsymbol{\beta}_{H})\mathcal{L}_{H}(\mu_{H},\boldsymbol{\beta}_{H})$ is the functional form to be considered for updating $\mu_{H}$ and $\boldsymbol{\beta}_{H}$.

S-1.5.4 **Mixture of additive and multiplicative TMCMC for updating

$\mu_{G}$, $\boldsymbol{\beta}_{G}$, $\mu_{G_{0}}$, $\boldsymbol{\beta}_{G_{0}}$, $\mu_{H}$ and $\boldsymbol{\beta}_{H}$ in a single block**

We shall update all the parameters $\mu_{G}$, $\boldsymbol{\beta}_{G}$, $\mu_{G_{0}}$, $\boldsymbol{\beta}_{G_{0}}$, $\mu_{H}$ and $\boldsymbol{\beta}_{H}$ using a mixture of additive and multiplicative TMCMC, where all the aforementioned parameters are given either the additive move or the multiplicative move with equal probability, and where the acceptance ratio will be calculated by evaluating the functional form

[TABLE]

at the numerator and the denominator corresponding to the proposed and the current values of $\mu_{G}$, $\boldsymbol{\beta}_{G}$, $\mu_{G_{0}}$, $\boldsymbol{\beta}_{G_{0}}$, $\mu_{H}$ and $\boldsymbol{\beta}_{H}$, with all other unknowns held fixed at their current values, multiplied by an appropriate Jacobian whenever the multiplicative move is chosen. For details regarding mixture of additive and multiplicative TMCMC, see ?.

S-2 A parallel algorithm for implementing our MCMC procedure

Recall that the mixtures associated with gene $j\in\{1,\ldots,J\}$, and individual $i\in\{1,\ldots,N_{k}\}$ and case-control status $k\in\{0,1\}$, are conditionally independent of each other, given the interaction parameters. This allows us to update the mixture components in separate parallel processors, conditionally on the interaction parameters. Once the mixture components are updated, we update the interaction parameters using a specialized form of TMCMC, in a single processor. Furthermore, the parameters of the HDP are also amenable to efficient parallelization. The details are as follows.

• (1)

1. (a)

In processes numbered [math] and $1$, simultaneously obtain the set of distinct elements $\boldsymbol{\Xi}_{R_{k},k}$; $k=0,1$, from $\{\phi_{tijk};t=1,\ldots,\tau_{ijk};i=1,\ldots,N_{k};j=1,\ldots,J\}$; $k=0,1$. 2. (b)

Communicate $\boldsymbol{\Xi}_{R_{k},k}$; $k=0,1$, to all the processes.

• (2)

1. (a)

In process [math], obtain the set of distinct elements $\boldsymbol{\zeta}_{S}$ from $\left\{\boldsymbol{\Xi}_{R_{k},k};k=0,1\right\}$. 2. (b)

Communicate $\boldsymbol{\zeta}_{S}$ to all the processes.

• (3)

In processes numbered [math] and $1$, do the following in parallel for $k=0,1$:

1. (a)

Simulate, following the retrospective method detailed in Section S-1.2.3,

$\boldsymbol{\eta}^{*}_{lk}\stackrel{{\scriptstyle iid}}{{\sim}}\frac{\alpha_{H}\tilde{\boldsymbol{H}}+\sum_{s=1}^{S}n_{sk}\delta_{\boldsymbol{\eta}_{s}}}{\alpha_{H}+n_{\cdot k}}$; $l=1,2,\ldots,\mathcal{L}$, for sufficiently large $\mathcal{L}$. 2. (b)

Communicate the simulated values to all the processes.

• (3)

Split $\left\{(j,k):~{}j=1,\ldots,J;~{}k=0,1\right\}$ in the available parallel processes.

1. (a)

For each $(j,k)$, simulate, following the retrospective method detailed in Section S-1.2.3, $\tilde{\boldsymbol{\xi}}_{ljk}\stackrel{{\scriptstyle iid}}{{\sim}}\frac{\alpha_{G_{0},k}\boldsymbol{H}_{k}+\sum_{l=1}^{R_{k}}n_{ljk}\delta_{\boldsymbol{\xi}_{lk}}}{\alpha_{G_{0},k}+n_{\cdot jk}}$; $l=1,2,\ldots,\mathcal{L}$. 2. (b)

Communicate the simulated values to all the processes.

• (4)

1. (a)

Split the triplets $\left\{(i,j,k):~{}i=1,\ldots,N_{k};~{}j=1,\ldots,J;~{}k=0,1\right\}$ in the available parallel processes sequentially into

[TABLE]

and

[TABLE] 2. (b)

Then parallelise updation of the mixtures associated with $\mathcal{T}_{1}$, followed by those of $\mathcal{T}_{2}$. 3. (c)

If, for any $(i,j,k)$, retrospective simulation from $[\boldsymbol{G}_{0,jk}|\boldsymbol{X}_{ijk}]$ requires more than $\mathcal{L}$ simulations of $\tilde{\boldsymbol{\xi}}_{ljk}$ in step (3) (a), then increase $\mathcal{L}$ to $\mathcal{L}^{*}$, and

1. (i)

For $k=0,1$, augment the simulations of $\left\{\boldsymbol{\eta}^{*}_{lk};l=1,\ldots,\mathcal{L}\right\}$ with new simulations $\left\{\boldsymbol{\eta}^{*}_{lk};l=\mathcal{L}+1,\ldots,\mathcal{L}^{*}\right\}$. 2. (ii)

For $j=1,\ldots,J$ and for $k=0,1$, augment the simulations of $\left\{\tilde{\boldsymbol{\xi}}_{ljk};l=1,\ldots,\mathcal{L}\right\}$ with new simulations $\left\{\tilde{\boldsymbol{\xi}}_{ljk};l=\mathcal{L}+1,\ldots,\mathcal{L}^{*}\right\}$. 3. (iii)

Repeat (4) (a) and (4) (b).

• (5)

During each MCMC iteration, for each $(i,j,k)$ in each available parallel processor, update the allocation variables $z_{ijk}$, the proportions $\boldsymbol{p}_{mijk}$; $m=1,\ldots,M$, and the missing data $\tilde{\boldsymbol{Y}}_{ijk}$, using the methods proposed in Sections S-1.3 and S-1.4.

• (6)

Communicate the results of updating in (4) and (5) to all the processes.

• (7)

1. (a)

During each MCMC iteration, update the parameters $\mu_{G}$, $\boldsymbol{\beta}_{G}$, $\mu_{G_{0}}$, $\boldsymbol{\beta}_{G_{0}}$, $\mu_{H}$ and $\boldsymbol{\beta}_{H}$ using additive TMCMC in a single block, as proposed in Section S-1.5, in process number [math]. 2. (b)

Communicate the updated results to all the processes.

S-3 Simulation studies

For simulation studies, we first generate realistic biological data for stratified population with known gene-environment interaction from the GENS2 software of ?. To this data, we then apply our model and methodologies in an effort to detect gene-environment interaction effects that are present in the data. We consider simulation studies in $5$ different true model set-ups: (a) presence of gene-gene and gene-environment interaction, (b) absence of genetic or gene-environmental interaction effect, (c) absence of genetic and gene-gene interaction effects but presence of environmental effect, (d) presence of genetic and gene-gene interaction effects but absence of environmental effect, and (e) independent and additive genetic and environmental effects.

As we demonstrate, our model and methodologies successfully identify the effects of the individual genes, gene-gene and gene-environment interactions, and the number of sub-populations. In all our applications, we set $M=30$, $\nu_{1}=\nu_{2}=1$, so that $\tilde{\boldsymbol{H}}$ is the uniform distribution on $[0,1]$. We set $\alpha_{G,ik}=0.1\times\exp\left(100+\mu_{G}+\beta_{G}E_{ik}\right)$, $\alpha_{G_{0},k}=0.1\times\exp\left(100+\mu_{G_{0}}+\beta_{G_{0}}\bar{E}_{k}\right)$ and $\alpha_{H}=0.1\times\exp\left(100+\mu_{H}+\beta_{H}\bar{\bar{E}}\right)$, where we assumed $\mu_{G},\mu_{G_{0}},\mu_{H}\stackrel{{\scriptstyle iid}}{{\sim}}U(0,1)$ and $\beta_{G},\beta_{G_{0}},\beta_{H}\stackrel{{\scriptstyle iid}}{{\sim}}U(-1,1)$. This structure ensured adequate number of sub-populations and reasonable mixing of MCMC. Figure S-1 shows some typical trace plots associated with our HDP model. Although the mixing is not excellent, there is no evidence of non-convergence.

S-3.1 First simulation study: presence of gene-gene and gene-environment interaction

S-3.1.1 Data description

As in ? we consider two genetic factors as allowed by GENS2 and simulated 5 data sets with gene-gene and gene-environment interaction with a one-dimensional environmental variable, associated with 5 sub-populations. One of the genes consists of 1084 SNPs and another has 1206 SNPs, with one disease pre-disposing locus (DPL) at each gene. There are 113 individuals in each of the 5 data sets, from which we selected a total of 100 individuals without replacement with probabilities assigned to the 5 data sets being $(0.1,0.4,0.2,0.15,0.15)$. Our final dataset consists of 46 cases and 54 controls. Since, in our case, the environmental variable is one-dimensional, $d=1$.

S-3.1.2 Model implementation

We implemented our parallel MCMC algorithm on 50 cores in a 64-bit VMware with 64-bit physical cores, each running at 2793.269 MHz. Our code is written in C in conjunction with the Message Passing Interface (MPI) protocol for parallelisation.

The total time taken to implement $30,000$ MCMC iterations, where the first $10,000$ are discarded as burn-in, is approximately 20 hours. We assessed convergence informally with trace plots, which indicated adequate mixing properties of our algorithm.

S-3.1.3 Specifications of the thresholds $\varepsilon$’s using null distributions

Following the method outlined in Section 4.2.4 and setting $M$ to be 30, we obtain $\varepsilon_{d^{*}}=0.200$, $\varepsilon_{\hat{d}_{1}}=0.167$, $\varepsilon_{\hat{d}_{2}}=0.167$, $\varepsilon_{d^{*}_{E}}=0.250$, $\varepsilon_{d^{*}_{E,1}}=0.185$, $\varepsilon_{d^{*}_{E,2}}=0.173$, $\varepsilon_{\beta_{G}}=0.874$, $\varepsilon_{\beta_{G_{0}}}=0.128$, $\varepsilon_{\beta_{H}}=0.219$.

S-3.1.4 Results of fitting our model

The posterior probabilities $P\left(d^{*}<\varepsilon_{d^{*}}|\mbox{Data}\right)$, $P\left(\hat{d}_{1}<\varepsilon_{\hat{d}_{1}}|\mbox{Data}\right)$ and $P\left(\hat{d}_{2}<\varepsilon_{\hat{d}_{2}}|\mbox{Data}\right)$ empirically obtained from $20,000$ MCMC samples, turned out to be $0.378$, $0.317$ and $0.324$, respectively. Hence, $H_{0,d^{*}}$, $H_{0,\hat{d}_{1}}$ and $H_{0,\hat{d}_{2}}$ are rejected, suggesting the influence of significant genetic effects in the case-control study.

However, $P\left(d^{*}_{E}<\varepsilon_{d^{*}_{E}}|\mbox{Data}\right)$, $P\left(\hat{d}_{E,1}<\varepsilon_{\hat{d}_{E,1}}|\mbox{Data}\right)$ and $P\left(\hat{d}_{E,2}<\varepsilon_{\hat{d}_{E,2}}|\mbox{Data}\right)$ are given, approximately, by $0.558$, $0.561$ and $0.550$, respectively, which seem to contradict the results of the clustering based hypothesis tests. This can be explained as follows. Since $\boldsymbol{G}_{0,jk}$ are discrete, the parameters $\boldsymbol{p}_{mijk}$, even if generated from $\boldsymbol{G}_{0,jk}$, coincide with positive probability, so that the effective dimensionalities of $\mbox{logit}\left(\bar{P}_{Mi_{0}jk=0}\right)$ and $\mbox{logit}\left(\bar{P}_{Mi_{1}jk=1}\right)$ are drastically reduced, so that the Euclidean distance between these two vectors is substantially small. As such, the Euclidean distance fails to reject the null even if it is false. As noted in ?, even the clustering metric in this scenario is not completely satisfactory since this involves clustering distance between two empirically obtained central clusterings which may not be very accurate unless the sample sizes for case and control are very large. However, compared to the Euclidean distance, the clustering metric turns out to be far more reliable.

To check the influence of the environmental variable on the genes we compute the posterior probabilities $P\left(|\beta_{G}|<\varepsilon_{\beta_{G}}|\mbox{Data}\right)$, $P\left(|\beta_{G_{0}}|<\varepsilon_{\beta_{G_{0}}}|\mbox{Data}\right)$ and $P\left(|\beta_{H}|<\varepsilon_{\beta_{H}}|\mbox{Data}\right)$. The probabilities turned out to be $0.544$, $0.550$ and $0.191$, respectively, showing that $\beta_{H}$ is very significant. That is, the environmental variable has a significant overall effect on the genes. Figure S-2 depicts the posterior probabilities of no gene-gene interactions for the controls and cases, showing the prominence of several gene-gene interactions in both control and case groups. As to be expected, in the case group, more instances of gene-gene interactions turned out to be significant compared to the control group.

The posteriors of the number of sub-populations, some of which are shown in Figure S-3, give high probabilities to $5$, the true number of sub-populations.

S-3.1.5 Detection of DPL

The correct positions of the DPL, provided by GENS2, are $rs13266634$ and $rs7903146$, for the first and second gene respectively. Due to the LD effects implied by the highly correlated structure of our current HDP based model, the actual DPL are difficult to locate. Notably, our model is considerably more structured than those of ? and ?, and any inappropriate dependence structure would render the task of DPL finding far more difficult than our previous models. Nevertheless, we demonstrate that our HDP model can detect DPLs with more precision compared to our previous matrix-normal-inverse-Wishart model for gene-environment interactions.

Following ? and ?, and writing $\boldsymbol{p}^{r}_{ijk}=\left\{p_{mijkr}:m=1,\ldots,M\right\}$, we declare the $r$-th locus of the $j$-th gene as disease pre-disposing if, for the $r$-th locus, the Euclidean distance $d^{r}_{j}\left(\mbox{logit}\left(\boldsymbol{p}^{r}_{i_{0}jk=0}\right),\mbox{logit}\left(\boldsymbol{p}^{r}_{i_{1}jk=1}\right)\right)$, between $\mbox{logit}\left(\boldsymbol{p}^{r}_{i_{0}jk=0}\right)$ and $\mbox{logit}\left(\boldsymbol{p}^{r}_{i_{1}jk=1}\right)$, is significantly larger than $d^{r^{\prime}}_{j}\left(\boldsymbol{p}^{r^{\prime}}_{i_{0}jk=0},\boldsymbol{p}^{r^{\prime}}_{i_{1}jk=1}\right)$, for $r^{\prime}\neq r$. We adopt the graphical method as in our previous works.

The red, horizontal lines in the panels of Figure S-4 represent the cut-off value such that the points above the horizontal line are those with the highest $2\%$ Euclidean distances. The actual DPLs of the two genes, as well as their nearest neighbours with Euclidean distances on or above the red, horizontal lines, are shown in the figures. That even such small sets of SNPs with highest 2% Euclidean distances consist of close neighbours of the true DPLs, is quite encouraging. Observe that the DPL detection is more precise for the second gene in the sense that the closest neighbour of the actual DPL above the red, horizontal line is closer to the true DPL than for the first gene.

The above results on DPL detection is also a significant improvement over ? where highest 10% Euclidean distances were considered, suggesting that our current HDP based model is more appropriate compared to our previous matrix-normal-inverse-Wishart model for gene-environment interaction.

S-3.2 Second simulation study: no genetic or environmental effect

Here we use the same case-control genotype data set as used by ? in their second simulation study where genetic effects are absent, consisting of 49 cases and 51 controls and 5 sub-populations with the mixing proportions $(0.1,0.4,0.2,0.15,0.15)$. We use the same environmental data set generated in our first simulation study described in Section S-3.1, which is unrelated to this genotype data.

Here we obtain $P\left(d^{*}<\varepsilon_{d^{*}}|\mbox{Data}\right)\approx 0.407$. Although this does not cross the $0.5$ benchmark, there is significant evidence in favour of the null, and falling short of $0.5$ can be attributed to the slight deficiency of the distance between the two approximate central clusterings associated with case and control, as already discussed in the context of the first simulation study.

Also, in this study, $P\left(|\beta_{G}|<\varepsilon_{\beta_{G}}|\mbox{Data}\right)$, $P\left(|\beta_{G_{0}}|<\varepsilon_{\beta_{G_{0}}}|\mbox{Data}\right)$ and $P\left(|\beta_{H}|<\varepsilon_{\beta_{H}}|\mbox{Data}\right)$ are given by $0.549$, $0.550$ and $0.649$, respectively, suggesting insignificance of the effect of the environmental variable on gene-gene interaction. As noted in ?, however, it is not straightforward to test whether or not the environment is responsible for the case-control status. This is because we have modeled the genotype data conditionally on case-control instead of modeling the case-control status conditionally on the environmental variable. ? use significance testing in a simple logistic regression framework to show insignificance of the environmental variable.

As before, our model assigned high posterior probability to 5 sub-populations.

Note that since there is no genetic effect in this study, the question of detecting DPLs does not arise here.

S-3.3 **Third simulation study: absence of genetic and gene-gene interaction effects

but presence of environmental effect**

In this study we consider a case-control genotype data set simulated from GENS2 where case-control status depends only upon the environmental data. The number of cases here is 47 and the number of controls is 53. This is the same case-control genotype data set as used by ? in their third simulation study.

In this case, we find that $P\left(d^{*}<\varepsilon_{d^{*}}|\mbox{Data}\right)\approx 0.400$, which provides reasonable evidence in favour of the null, even though the 0.5 benchmark is not crossed. Moreover, $P\left(|\beta_{G}|<\varepsilon_{\beta_{G}}|\mbox{Data}\right)\approx 0.536$, $P\left(|\beta_{G_{0}}|<\varepsilon_{\beta_{G_{0}}}|\mbox{Data}\right)\approx 0.518$ and $P\left(|\beta_{H}|<\varepsilon_{\beta_{H}}|\mbox{Data}\right)\approx 0.504$, suggesting that the environmental variable does not affect the genetic structure. ? show by means AIC, in the context of simple logistic regression, that the best model consists of the marginal effects of the second gene and the environment. In conjunction with our HDP-based model which produces reasonable evidence in favour of accepting the hypothesis of no genetic effect, it may be possible to conclude that the environmental variable is responsible for the case-control status.

As before, 5 subpopulations get significant weight by our posterior distribution, and again, the question of DPL detection is irrelevant here since there is no genetic effect.

S-3.4 **Fourth simulation study: presence of genetic and gene-gene interaction effects

but absence of environmental effect**

Here we use the same genotype data set as used by ? in their first simulation study associated with genetic and gene-gene interaction effects, consisting of 41 cases and 59 controls and 5 sub-populations with the mixing proportions $(0.1,0.4,0.2,0.15,0.15)$. We use the same environmental data set generated in our first simulation study described in Section S-3.1, which is unrelated to this case-control genotype data.

Here we obtain $P\left(|\beta_{G}|<\varepsilon_{\beta_{G}}|\mbox{Data}\right)\approx 0.549$, $P\left(|\beta_{G_{0}}|<\varepsilon_{\beta_{G_{0}}}|\mbox{Data}\right)\approx 0.542$ and $P\left(|\beta_{H}|<\varepsilon_{\beta_{H}}|\mbox{Data}\right)\approx 0.552$, correctly suggesting insignificance of the environmental variable with respect to its effect on the genetic structure. Using logistic regression, ? conclude that the environmental variable has no role to play in the case-control status. Furthermore, we obtain $P\left(d^{*}<\varepsilon_{d^{*}}|\mbox{Data}\right)\approx 0.390$, $P\left(\hat{d}_{1}<\varepsilon_{\hat{d}_{1}}|\mbox{Data}\right)\approx 0.336$ $P\left(\hat{d}_{2}<\varepsilon_{\hat{d}_{2}}|\mbox{Data}\right)\approx 0.324$. so that importance of genes is correctly indicated by our tests. Figure S-5 shows the posterior probabilities of no gene-gene interactions for controls and cases. Interestingly, there seems to be no gene-gene interaction in the control group and only two (marginal) instances of gene-gene interaction among the cases.

Figure S-6 shows the plots of Euclidean distances between cases and controls for the loci of the two genes. In this case, Gene-1 has been located quite precisely, and for Gene-2 the Euclidean distance for even the true DPL is very close to the red, horizontal line, indicating encouraging performance.

S-3.5 Fifth simulation study: independent and additive genetic and environmental effects

As in ?, we consider the situation where the genetic and environmental effects are independent of each other and additive; the data consists of 57 cases and 43 controls.

Note that, as in ?, in our current HDP-based Bayesian model also there is no provision for additivity of genetic and environmental effects. As such, it is not expected to capture the true data-generating mechanism accurately. Indeed, here we obtain $P\left(d^{*}<\varepsilon_{d^{*}}|\mbox{Data}\right)\approx 0.389$, $P\left(\hat{d}_{1}<\varepsilon_{\hat{d}_{1}}|\mbox{Data}\right)\approx 0.337$ and $P\left(\hat{d}_{2}<\varepsilon_{\hat{d}_{2}}|\mbox{Data}\right)\approx 0.331$, indicating significance of the genes. However, the test with $d^{*}_{E}$ does not yield overwhelming evidence against the null. Our tests of gene-gene interaction, as depicted in Figure S-7, indicate significant interactions for controls and particularly for cases.

Also, $P\left(|\beta_{G}|<\varepsilon_{\beta_{G}}|\mbox{Data}\right)$, $P\left(|\beta_{G_{0}}|<\varepsilon_{\beta_{G_{0}}}|\mbox{Data}\right)$ and $P\left(|\beta_{H}|<\varepsilon_{\beta_{H}}|\mbox{Data}\right)$ are given, approximately, by $0.547$, $0.550$ and $0.367$, the last value showing that the environmental variable does affect gene-gene interaction. The lack of the additivity provision in our model seems to have forced the gene-environment interaction in this case.

In spite of the lack of additivity of our model the Euclidean distances between cases and controls for the gene-wise SNPs are not adversely affected, and the actual DPLs are detected quite accurately; see Figure S-8. This brings forth the generality and usefulness of our nonparametric dependence structure.

As before, 5 sub-populations receive significant posterior probabilities.