CoMM: A Collaborative Mixed Model That Integrates GWAS and eQTL Data Sets to Investigate the Genetic Architecture of Complex Traits

Abstract

Genome-wide association study (GWAS) analyses have identified thousands of associations between genetic variants and complex traits. However, it is still a challenge to uncover the mechanisms underlying the association. With the growing availability of transcriptome data sets, it has become possible to perform statistical analyses targeted at identifying influential genes whose expression levels correlate with the phenotype. Methods such as PrediXcan and transcriptome-wide association study (TWAS) use the transcriptome data set to fit a predictive model for gene expression, with genetic variants as covariates. The gene expression levels for the GWAS data set are then ‘imputed’ using the prediction model, and the imputed expression levels are tested for their association with the phenotype. These methods fail to account for the uncertainty in the GWAS imputation step, and we propose a collaborative mixed model (CoMM) that addresses this limitation by jointly modelling the multiple analysis steps. We illustrate CoMM’s ability to identify relevant genes in the Northern Finland Birth Cohort 1966 data set and extend the model to handle the more widely available GWAS summary statistics.

Comment on: Yang C, Wan X, Lin X, Chen M, Zhou X, Liu J. CoMM: a collaborative mixed model to dissecting genetic contributions to complex traits by leveraging regulatory information. Bioinformatics. 2018;35:1644-1652. doi:10.1093/bioinformatics/bty865. PubMed PMID: 30295737. https://www.ncbi.nlm.nih.gov/pubmed/30295737

Introduction

In the last decade, genome-wide association studies (GWASs) have identified thousands of genetic variants associated with complex traits. However, an understanding of the mechanisms linking a genetic variant to a complex trait is relatively limited. As many GWAS loci are located outside coding regions and regulatory variation plays an important role in shaping observed traits, gene expression has been proposed as an informative intermediate phenotype. ¹ Incorporating regulatory information into statistical analyses provides us with a principled approach to study the genetic contribution to complex traits through the regulation of gene expression.

One approach to incorporate functional information is to do so without explicitly modelling the relationship between gene expression and phenotype. Sequence Kernel Association Tests (SKATs), ² for instance, can be used to prioritize genetic variants based on functional annotation. It is based on the idea that genetic variants with known biological functions are more likely to be associated with a trait. Hence, when testing for association between genetic variants and a trait, the genetic variants are prioritized by placing larger weights on those with known functions. ²

Recently, the growing availability of transcriptome data has given rise to methods that evaluate genetically regulated gene expression using both GWAS and transcriptome data sets. Large-scale transcriptome data sets, which contain information on genotypes and gene expression levels, include the Genotype-Tissue Expression Consortium (GTEx), ³ the Genetic European in Health and Disease (GEUVADIS) Project, ⁴ Braineac, ⁵ Depression Genes and Networks ⁶ and eQTLGen. ⁷

Methods that have leveraged on both GWAS and transcriptome data include PrediXcan^8,9 and transcriptome-wide association study (TWAS) ¹⁰ and generally proceed in 3 steps. First, using transcriptome data, they fit predictive models for gene expression with genetic variants near a gene as covariates. PrediXcan proposed the use of elastic net regression or ridge regression to build a predictive model, while TWAS proposed the use of a linear mixed model. The fitted models are then used to predict the gene expression levels for individuals in the GWAS data set. Finally, a simple linear regression is used to examine the association between the predicted expression levels and the complex trait in the GWAS data set.

Methods that proceed in such a stage-wise manner do not account for the uncertainty that arises when imputing the gene expression levels in the GWAS data set, which may lead to a loss in statistical power. To address this limitation, we proposed the collaborative mixed model (CoMM), ¹¹ which accounts for the uncertainty in the ‘imputation’ model by jointly fitting the imputation and the association analysis models.

Method

Suppose that the gene expression levels for $G$ genes and the allele counts for $M$ single-nucleotide polymorphisms (SNPs) are measured for $n_1$ samples in the transcriptome data set and that the phenotype values and the allele counts for the same $M$ SNPs are measured for $n_2$ samples in the GWAS data set. The transcriptome data set consists of the $n_1$ by $G$ gene expression data matrix $\mathbf{Y}$ and the $n_1$ by $M$ genotype data matrix ${\mathbf{W}}_1$ . The GWAS data set consists of the phenotype vector $\mathbf{z}$ and the $n_2$ by $M$ genotype data matrix ${\mathbf{W}}_2$ . We predict the gene expression levels one gene at a time and denote the gene expression levels at the gth gene by ${\mathbf{y}}_g$ . As we are interested in the variation in gene expression attributable to variation in its cis-SNPs, we model the gene expression levels and phenotype value using only nearby SNPs. Let ${\mathbf{W}}_{1g}$ and ${\mathbf{W}}_{2g}$ denote the genotype matrix corresponding to the gene’s nearby SNPs, in the transcriptome data set and the GWAS data set, respectively. Let $M_g$ denote the number of SNPs corresponding to the gth gene. We assume that ${\mathbf{y}}_g$ is mean centred, and ${\mathbf{W}}_{1g}$ and ${\mathbf{W}}_{2g}$ are standardized (columns have zero mean and unit variance).

In the CoMM, we first model the relationship between gene expression ${\mathbf{y}}_g$ and genotype ${\mathbf{W}}_{1g}$ in the transcriptome data set

$$\begin{array}{c}{\mathbf{y}}_g={\mathbf{W}}_{1g}{\mathrm{\gamma }}_g+{\mathbf{e}}_1\end{array}$$

(1)

where ${\mathrm{\gamma }}_g$ is an $M_g$ vector of SNP effects on the gene expression level, and ${\mathbf{e}}_1\sim N(0,{\mathrm{\sigma }}_1^2{\mathbf{I}}_{n_1})$ is an $n_1$ vector representing the error associated with the expression level. Next, we model the relationship between phenotype $\mathbf{z}$ and genotype ${\mathbf{W}}_{2g}$ in the GWAS data set as

$$\begin{array}{c}z=X\mathrm{\beta }+{\mathrm{\alpha }}_g{\mathbf{W}}_{2g}{\mathrm{\gamma }}_g+{\mathbf{e}}_2\end{array}$$

(2)

where $\mathbf{X}$ contains covariates that control for population stratification and other confounding variables, $\mathrm{\beta }$ is a vector of fixed-effects coefficients, and ${\mathbf{e}}_2\sim N(0,{\mathrm{\sigma }}_2^2{\mathbf{I}}_{n_2})$ is an $n_2$ vector representing the error associated with the phenotype.

The quantity of interest is ${\mathrm{\alpha }}_g$ , the effect of gene g’s expression level on the phenotype. In addition, we assume the prior distribution on ${\mathrm{\gamma }}_g$

$$\begin{array}{c}{\mathrm{\gamma }}_g\sim N\left(0,\mathrm{\sigma }{}_g^2{\mathbf{I}}_{M_g}\right)\end{array}$$

(3)

effectively treating the effects of genotype on gene expression as random. An accelerated expectation–maximization (EM) algorithm using parameter expansion ¹² is used to estimate all parameters in the joint model given by equations (1) and (2). Figure 1 summarizes the data input, the joint model, and output for CoMM.

Figure 1.

Schematic of CoMM. The transcriptome and GWAS data sets are used to fit the parameters in the model given by equations (1) and (2). The parameter estimates are used to evaluate the likelihood ratio test statistic, which tests for association between the phenotype and genetically regulated gene expression. CoMM indicates collaborative mixed model; GWAS, genome-wide association studies.

As our objective is to evaluate whether a gene is associated with the phenotype via gene expression, we perform the hypothesis test

$$\begin{array}{c}{\mathscr{H}}_0:{\mathrm{\alpha }}_g=0,v.s.{\mathscr{H}}_1:{\mathrm{\alpha }}_g\ne 0\end{array}$$

(4)

The likelihood ratio test is used, and the test statistic is given by

$${\mathrm{\Lambda }}_g=2\left(\log \mathrm{Pr}\left({\mathbf{y}}_g,\mathbf{z},{\mathrm{\gamma }}_g|\widehat{\mathrm{\theta }}\right)-\log \mathrm{Pr}\left({\mathbf{y}}_g,\mathbf{z},{\mathrm{\gamma }}_g|{\widehat{\mathrm{\theta }}}_0\right)\right)$$

where $\widehat{\mathrm{\theta }}$ and ${\widehat{\mathrm{\theta }}}_0$ contain the parameter estimates obtained under the full model and ${\mathscr{H}}_0$ , respectively. The test statistic ${\mathrm{\Lambda }}_g$ is asymptotically distributed as ${\chi }_{df=1}^2$ under the null hypothesis. The key thing to note is that the likelihood reflects the uncertainty in both the imputation and association analysis models (equations (1) and (2)). As such, the CoMM test statistic for expression-trait association takes into account the uncertainty in the imputation model.

Results in the Northern Finland Birth Cohort 1966

We analysed the GWAS data set from the Northern Finland Birth Cohort 1966 (NFBC1966) ¹³ with the aid of transcriptome data from GTEx (tissue: subcutaneous adipose). ³ The NFBC1966 data set records traits such as body mass index (BMI), low-density lipoprotein cholesterol (LDL), triglycerides (TGs), total cholesterol (TC), systolic blood pressure (SysBP), and diastolic blood pressure (DiaBP).

The CoMM returns a larger number of significant findings than PrediXcan and SKAT, as indicated by the QQ-plots of the P values (Figure 2). In particular, CoMM reported 12 significant genes associated with triglyceride (TG) levels, whereas PrediXcan:Enet, PrediXcan:Ridge, and SKAT reported 2, 1, and 0 significant genes, respectively. Among the 12 identified genes, 2 (OST4 and EIF2B4) have nonnegligible cellular heritability ( $h_C^2=2.03\%\mathrm{and}1.33\%$ , respectively) and have reported associations with TG in previous studies.^14,15 In this instance, CoMM performs better than SKAT due to the use of gene regulation information in the GTEx data, and outperforms PrediXcan by taking into account the uncertainty in the imputation model.

Figure 2.

The QQ-plots of P-values for the quantitative traits in NFBC1966.

Extension of CoMM to Analyse GWAS Summary Data

A limitation of CoMM is that it requires individual-level data and is unable to make use of large-scale GWAS that provide only summary statistics. To capitalize on these GWAS, we extend CoMM so that summary statistics, in the form of estimated SNP effect sizes and their variances, can replace the role of individual-level GWAS data.

We adapt the method of Zhu and Stephens, ¹⁶ which made use of summary statistics by introducing a regression with summary statistics (RSS) likelihood in a Bayesian framework. As gene expression levels are modelled using multiple SNPs, we additionally require information on the correlations among SNPs (linkage disequilibrium). Fortunately, such information is available in public data sets such as the 1000 Genomes Project Consortium. ¹⁷ In CoMMs for GWAS summary statistics (CoMM-S²), ¹⁸ the association between phenotype and genetically regulated gene expression is evaluated by combining the distribution for individual-level transcriptome data with an RSS distribution for GWAS summary statistics, while taking into account linkage disequilibrium as estimated from a reference panel.

Even though CoMM-S² utilizes GWAS summary statistics, it has comparable performance as CoMM. To illustrate the performance of CoMM-S² relative to CoMM, we use NFBC1966 as the GWAS data set, GTEx as the transcriptome data set and the 1000 Genomes Project as a reference panel to estimate linkage disequilibrium. In general, the test statistic values from CoMM-S² are close to their corresponding values from CoMM: the regression slope is around 1 and R² ranges from 0.91 to 0.99 (Figure 3). The close correspondence in test statistic values is most apparent in the null region. In the nonnull region, the test statistics for CoMM-S² may be inflated (Figure 3). One possible reason for this is linkage disequilibrium misspecification, ¹⁸ due to the genetic differences between the Finnish cohort in NFBC1966 and the European sample in the 1000 Genomes Project. ¹⁹ Nonetheless, as a strong inflation occurs only when the CoMM statistics is large, CoMM-S² maintains a reasonable false-positive rate in the presence of misspecified linkage disequilibrium.

Figure 3.

Scatter plot of test statistics from the likelihood ratio test for CoMM-S² versus CoMM, using GTEx (tissue: subcutaneous adipose) transcriptome data, NFBC1966 GWAS data, and 1000 Genomes Project reference panel data.

Finally, we note that both CoMM and CoMM-S² are designed for single-tissue analysis. When a multi-tissue transcriptome data set is available, an approach that takes into account genetic correlation across tissues may be preferable. Such an approach would be better equipped to identify biologically relevant tissues for each gene, and may also provide an increase in statistical power for tissues that are difficult to obtain. ²⁰ Recently, 2 multi-tissue approaches, UTMOST ²⁰ and MultiXcan, ²¹ have been proposed. They are more powerful than single-tissue approaches in expression-trait association analyses. However, they ignore the uncertainty due to the imputation step, and similar to what has been proposed for CoMM and CoMM-S², the ability to detect relevant genes can be further improved by combining the imputation model and association analysis model via a unified likelihood framework. This remains a promising avenue for further research.