Bayesian hierarchical modeling of means and covariances of gene expression data within families

Abstract

We describe a hierarchical Bayes model for the influence of constitutional genotypes from a linkage scan on the expression of a large number of genes. The model comprises linear regression models for the means in relation to genotypes and for the covariances between pairs of related individuals in relation to their identity-by-descent estimates. The matrices of regression coefficients for all possible pairs of single-nucleotide polymorphisms (SNPs) by all possible expressed genes are in turn modeled as a mixture of null values and a normal distribution of non-null values, with probabilities and means given by a third-level model of SNP and trait random effects and a spatial regression on the distance between the SNP and the expressed gene. The latter provides a way of testing for cis and trans effects. The method was applied to data on 116 SNPs and 189 genes on chromosome 11, for which Morley et al. (Nature 2004, 430: 743–747) had previously reported linkage. We were able to confirm the association of the expression of HSD17B12 with a SNP in the same region reported by Morley et al., and also detected a SNP that appeared to affect the expression of many genes on this chromosome. The approach appears to be a promising way to address the huge multiple comparisons problem for relating genome-wide genotype × expression data.

Background

Recent advances in genomic technology now allow genotyping of hundreds of thousands of single-nucleotide polymorphisms (SNPs) and measurement of the expression of tens of thousands of genes on single microarrays or chips at a manageable cost. Extensive literature on the analysis of gene expression data has evolved over the last five years, and since the advent of ultra-high-volume genotyping platforms, genome-wide association and linkage scans using SNPs have also become feasible. The multiple comparisons problem is central to the analysis of either type of high-volume data. In 2001, Jansen and Nap [1] proposed combining the analysis of the two technologies in what he called "genetical genomics" to provide insight into the genetic regulation of gene expression.

However, only quite recently have attempts been made to relate the two technologies, first by Morley et al. [2] in a linkage scan for 3554 expressed genes in relation to 2756 autosomal SNP markers, and subsequently by the same group [3] in a genome-wide association scan of 27 of the expressed genes with the highest linkage in the first study, in relation to >770,000 SNPs. (See also Schadt et al. [4] and Stranger et al. [5] for similar analyses.) Independently, Tsalenko et al. [6] proposed a biclustering method to visualize SNPs and the transcripts they regulate, using an approach that is more visual than statistical. The multiple comparisons problem in such analyses (2.7 billion comparisons in the association analysis) dwarfs those from either genome-wide linkage or association analyses of single traits or supervised cluster analyses of expression data in relation to single outcomes.

Therefore, there is a need to develop new statistical methods to analyze all transcripts and genotypes together. Here, we describe a novel hierarchical Bayesian approach to the analysis of all possible pairs of associations and linkages between expressed genes and SNP markers. We demonstrate the results for chromosome 11 and we argue that the method can be extended to cover the entire genome and transcriptome.

Methods

Statistical model

Let Y _ij ⁿdenote the expression of gene n in member j of family i and let G _ij ^mbe the corresponding SNP genotype at marker m at location x _m . For the means and covariances of the expression traits, we adopted a generalized estimating equations model of the form used by Thomas et al. [7]

E(Y _ij ⁿ) ≡ μ _ij ⁿ= α₀ⁿ+ Σ _m A^nmG _ij ^m (1)

E(C _ijk ⁿ) ≡ χ _ijk ⁿ= β₀ⁿ+ BⁿZ _ijk (Xⁿ), (2)

where C _ijk ⁿ= (Y _ij ⁿ- μ _ij ⁿ)(Y _ik ⁿ- μ _ik ⁿ) and Z _ijk (x) are the estimated E(IBD _ijk (x)|G _i ) at chromosomal location x for pairs (j, k) from nuclear family i, based on the complete multilocus marker data. Xⁿis a latent variable for location of the unobserved causal locus linked to expression trait n. For j = k, V(Y _ij ⁿ) = χⁿmodels the gene expression variance in Eq. (2).

In Eq. (1), the regression coefficients A^nmare modeled as a mixtures of null values with probabilities 1-π^nmand a normal distribution of non-null values with means α^nmexpressed in terms of row and column effects:

A^nm~ (1 - π^nm) δ (0) + π^nmN(α^nm, σ²), (3)

where

α^nm= γ₀^A+ γ₁^AI(x _m ∈ R _n ) + e _m ^A+ h _n ^A (4a)

logit(π^nm) = γ₀^P+ γ₁^PI(x _m ∈ R _n ) + e _m ^P+ h _n ^P. (4b)

The parameter γ₁ distinguishes between cis and trans effects, a cis interaction occurs when the chromosomal location x _m of SNP m is within the interval R _n , the alignment region for the gene expression probe n. The random effects e and h are distributed as

(e _m ^A, e _m ^P) ~ N₂(0, T) (5a)

(h _n ^A, h _n ^P) ~ N₂(0, W) (5b)

and the γs, T, W have uninformative normal and Wishart priors.

The regression coefficients Bⁿin the covariance model in Eq. (2) are handled similarly, except that we assume each trait has at most one region linked to it. (This is not essential to the method, because Eq. (2) could be extended to a summation over multiple independent linkage regions, but it would not make sense to offer all marker locations simultaneously, since the IBD variables are highly correlated from one location to the next.) Thus, we assume

Bⁿ~ N [γ₀^B+ γ₁^BI(Xⁿ∈ R _n ), τ²] (6)

and pick a uniform prior on Xⁿ; to simplify the calculations, we restrict Xⁿto the observed marker locations x _m and compute IBD probabilities only at these locations. Thus, Xⁿhas a discrete distribution with prior masses inversely proportional to the local marker density, here estimated simply as |x_m+1 - x_m-1|. The full model is represented in the directed acyclic graph shown in Figure 1.

Figure 1

Directed acyclic graph for the analysis model. Squares represent observed data, circles represent parameters or latent variables, triangles represent deterministic nodes.

We fitted the model using a Markov-chain Monte Carlo (MCMC) approach, implemented in Matlab. Updates of all parameters except the location parameters Xⁿvalues involve standard Gibbs sampling from their respective full conditional distributions, e.g., [α₀ⁿ|Yⁿ, G, A], [β₀ⁿ|Cⁿ, Z, B], [A^nm|Yⁿ, G^m, α₀ⁿ, π^nm, α^nm, σ²], etc. The updates of the X values are based on a Metropolis-Hastings procedure with a random walk proposal. The sequence was started ten times from several initial points chosen from an overdispersed prior around rough estimates. Half of the initial samples are discarded and the second half is kept. The number of kept samples, L = 4000, is chosen to be large enough so that for all parameters of interest the variance between sequences V _B is comparable to that within sequence V _W , R < 1.10:

\widehat{R}=\sqrt{\frac{L-1}{L}+\frac{1}{L}\frac{V_B}{V_W}}.

(7)

The rationale behind this convergence monitoring procedure is described and justified by Gelman et al. [8].

Subjects, genotypes, and phenotypes

In order to keep the computation to a manageable level, we restricted this analysis to the SNP genotypes and expressed genes on chromosome 11, since previous analyses by Morley et al. [2] had found evidence of both cis and trans linkages at this chromosome. The final data set thus had 116 SNPs and 189 expressed genes. IBD status was estimated from the complete two-generation pedigrees (excluding grandparents) by a program written by JM based on the Lander-Green algorithm [9]. All 378 sib pairs (110 individuals) from the available 14 families were included in the phenotype analysis.

Results

After convergence has been reached, the number of regression coefficients with nonzero coefficients in Eq. (1) is very small. This is because in the mixture model employed in Eq. (3), a large number of the probabilities are close to 0 (Figure 2).

Figure 2

Gene expression × Genotype associations and residual linkage summary. Left, Image describing the mean value of the association parameters π^nm between the gene expression phenotypes (rows) and the SNP genotypes (columns). The matrix shows that the interactions are very sparse (dark spots), meaning that phenotypes are controlled by small number of SNPs, with no apparent concentration along the cis region delimited by blue lines. However, there exist some SNPs (columns) that seem to be correlated with a large set of phenotypes, potentially indicating a master regulatory region. Right, Image describing the posterior probability of the linkage locus after removing the association effect from the covariance.

Figure 2 also shows that each gene expression phenotype is explained by relatively few genotypes that have a role in regulating their expression. Table 1 lists, for the best predicted phenotypes, the SNPs included most frequently in the model. Significantly, the top ranking phenotype, HSD17B12 (217869_at), associated with SNP rs1453389, is the same as the one reported by Cheung et al. as associated with another SNP in the same region (not included in the GAW data set). Figure 3 shows that some SNPs in chromosome 11, especially rs916482, are significantly associated with more phenotypes than others. These SNP may be within a master regulatory region of gene expression. The list of gene ontology terms that were over-represented in the list of its associated genes involved mostly metabolic functions (Figure 4).

Figure 3

Potential master regulatory region around rs916482 SNP. Bottom plot is the cross-section of column 84 of Figure 1, describing the association between all phenotypes in chromosome 11 and SNP m = rs916482. The top plot shows the sign of dependence on the genotype. This SNP has a large number of associated genotypes, providing a strong indication of a master regulatory region.

Figure 4

Gene ontology on potential master regulatory region. Overrepresented GO terms by the phenotypes associated to the SNP rs916482 analyzed using FatiGO http://www.fatigo.org/.

Figure 5

Linkage of residual gene expression variation after association. Phenotypes ranked by most significant coefficient of determination in the covariance model, the posterior distribution of their locus position X _n , and its mode. The coefficient of determination and its significance are calculated from samples drawn around the mode (10%).

Table 1 Top-ranking associations

The covariance model [Eq. (2)] results are summarized in the right panel of Figure 2, and the strongest linkage peaks are listed in Table 2. This linkage is for the remaining variation not explained by the association/means model [Eq. (1)], and the peaks would correspond to unseen genotypes that are in LD with a marker that was not used in the association model. Thus, this explains in part why linkage results are less compelling than the association ones. However, for those phenotypes for which significant linkage was found, the expression covariance increased with the IBD status, especially in 208964_s_at.

Table 2 Linkage of residual gene expression variation after association

Discussion

We have introduced a novel hierarchical Bayes model for genetic control of gene expression. Our approach to dealing with the multiple comparisons problem is to represent the matrices of all possible SNP × expressed gene association or linkage coefficients in terms of row and column random effects, along with a spatial regression on the distance between the two. Although this allows inference on specific pairs, we have greater interest in the variances of the row and column effects, which reflect systematic tendencies for SNPs to affect variable numbers of phenotypes and for phenotypes to be differentially expressed. Our mixture model also supports the possibility that the vast majority of such associations or linkages would be truly null, and allows separate estimation of both the probability and magnitude of non-null tests. So far we have not imposed any relationship between the parameters of the association (means) and linkage (covariance) models, but one might contemplate using the broad regions where linkage is seen for a particular phenotype as a prior for testing single-SNP associations with that phenotype.

The strongest gene-expression × SNP association reported by Cheung et al. [3] on chromosome 11 also appeared in our results as the most significant association, but with a SNP close to theirs (their reported SNP was not included in the data set). We also found evidence of at least one SNP that appears to be linked to a large number of expressed genes, suggesting the existence of master regulatory genes in that region.

We chose to restrict these analyses to a subset of genes and SNPs on a single chromosome to test the feasibility of the method. In principle the approach could be applied on a genome-wide scale, since the computation time increases linearly with m, n, sample size, and number of MCMC samples. Generating 4000 MCMC samples required 6 hours on a 2.2 GHz single-processor machine. However, one outstanding methodological challenge that would have to be addressed before the approach could be applied to dense SNP associations would be how to deal with the multicollinearity problem; for this reason, we chose to restrict this analysis to only a subset of SNPs that were not in strong LD with each other.