Bayesian hierarchical modeling of means and covariances of gene expression data within families
© Pique-Regi et al; licensee BioMed Central Ltd. 2007
Published: 18 December 2007
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.
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  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.  in a linkage scan for 3554 expressed genes in relation to 2756 autosomal SNP markers, and subsequently by the same group  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.  and Stranger et al.  for similar analyses.) Independently, Tsalenko et al.  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.
Let Y ij n 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. 
E(Y ij n ) ≡ μ ij n = α0 n + Σ m A nm G ij m (1)
E(C ijk n ) ≡ χ ijk n = β0 n + B n Z ijk (X n ), (2)
where C ijk n = (Y ij n - μ ij n )(Y ik n - μ ik n ) 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 n is a latent variable for location of the unobserved causal locus linked to expression trait n. For j = k, V(Y ij n ) = χ n 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 , σ2), (3)
α nm = γ0 A + γ1 A I(x m ∈ R n ) + e m A + h n A (4a)
logit(π nm ) = γ0 P + γ1 P I(x m ∈ R n ) + e m P + h n P . (4b)
The parameter γ1 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 ) ~ N2(0, T) (5a)
(h n A , h n P ) ~ N2(0, W) (5b)
and the γs, T, W have uninformative normal and Wishart priors.
The regression coefficients B n 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 ~ N [γ0 B + γ1 B I(X n ∈ R n ), τ2] (6)
The rationale behind this convergence monitoring procedure is described and justified by Gelman et al. .
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.  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 . All 378 sib pairs (110 individuals) from the available 14 families were included in the phenotype analysis.
Top SNPs used in the prediction
P(R2 > 0)
Linkage of residual gene expression variation after association
P(R2 > 0)
Mode [X n ]
Pr [X n = m]
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.  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.
Supported in part by NIH grant R01-GM58897 and U01-ES 015090. Implemented Matlab software is available upon request.
This article has been published as part of BMC Proceedings Volume 1 Supplement 1, 2007: Genetic Analysis Workshop 15: Gene Expression Analysis and Approaches to Detecting Multiple Functional Loci. The full contents of the supplement are available online at http://www.biomedcentral.com/1753-6561/1?issue=S1.
- Jansen RC, Nap JP: Genetical genomics: the added value from segregation. Trends Genet. 2001, 17: 388-391. 10.1016/S0168-9525(01)02310-1.View ArticlePubMedGoogle Scholar
- Morley M, Molony CM, Weber TM, Devlin JL, Ewens KG, Spielman RS, Cheung VG: Genetic analysis of genome-wide variation in human gene expression. Nature. 2004, 430: 743-747. 10.1038/nature02797.View ArticlePubMed CentralPubMedGoogle Scholar
- Cheung VG, Spielman RS, Ewens KG, Weber TM, Morley M, Burdick JT: Mapping determinants of human gene expression by regional and genome-wide association. Nature. 2005, 437: 1365-1369. 10.1038/nature04244.View ArticlePubMed CentralPubMedGoogle Scholar
- Schadt EE, Monks SA, Drake TA, Lusis AJ, Che N, Colinayo V, Ruff TG, Milligan SB, Lamb JR, Cavet G, Linsley PS, Mao M, Stoughton RB, Friend SH: Genetics of gene expression surveyed in maize, mouse and man. Nature. 2003, 422: 297-302. 10.1038/nature01434.View ArticlePubMedGoogle Scholar
- Stranger BE, Forrest MS, Clark AG, Minichiello MJ, Deutsch S, Lyle R, Hunt S, Kahl B, Antonarakis SE, Tavaré S, Deloukas P, Dermitzakis ET: Genome-wide associations of gene expression variation in humans. PLoS Genet. 2005, 1: e78-10.1371/journal.pgen.0010078.View ArticlePubMed CentralPubMedGoogle Scholar
- Tsalenko A, Sharan R, Edvardsen H, Kristensen V, Børresen-Dale AL, Ben-Dor A, Yakhini Z: Analysis of SNP-expression association matrices. Proc IEEE Comput Syst Bioinform Conf. 2005, 135-143.Google Scholar
- Thomas DC, Qian D, Gauderman WJ, Siegmund K, Morrison JL: A generalized estimating equations approach to linkage analysis in sibships in relation to multiple markers and exposure factors. Genet Epidemiol. 1999, 17 (Suppl 1): S733-S734.Google Scholar
- Gelman A, Carlin JB, Stern HS, Rubin DB, (Eds): Bayesian Data Analysis. 2004, New York: Chapman & Hall/CRC, 2Google Scholar
- Lander ES, Green P: Construction of multilocus genetic linkage maps in humans. Proc Natl Acad Sci USA. 1987, 84: 2363-2367. 10.1073/pnas.84.8.2363.View ArticlePubMed CentralPubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.