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

Statistical model

Let Y _ij ⁿ denote the expression of gene n in member j of family i and let G _ij ^m be 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 ^nm G _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 ^nm are modeled as a mixtures of null values with probabilities 1-π ^nm and a normal distribution of non-null values with means α ^nm expressed in terms of row and column effects:

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

where

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

logit(π ^nm ) = γ₀ ^P + γ₁ ^P I(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 + γ₁ ^B I(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.

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.