Inferring genetic architecture of complex traits using Bayesian integrative analysis of genome and transcriptome data

Large amounts of genomic information generated from Single Nucleotide Polymorphism (SNP) microarrays have become available in recent years for many species[1–3]. This genomic information is used to detect polymorphisms that contribute to variation in economically important traits, such as production traits in farm animals[3]. Microarray technology is also used to screen the expression levels of thousands of genes, i.e., the transcriptome[4, 5]. Studies have shown that genetic background can have a large impact on differential expression[6]. Integrating genome and transcriptome information can help to elucidate the underlying biology of the genotype-phenotype map, using expression Quantitative Trait Locus (eQTL) mapping[7].

However, in the eQTL approach, associations between SNPs, transcript level, and phenotypes are analyzed individually. This is likely to lead to “missing heritability”[8], because corrections for multiple testing lead to a high false negative rate and multiple SNPs and transcript level that jointly explain the phenotype are ignored[9, 10]. Here we propose and demonstrate Bayesian models that model all SNPs and transcript level simultaneously to obtain explained variances by the whole genome and whole transcriptome. In these models, we identify eQTLs as those SNPs whose effects disappear when transcript level are added to the model. Genomic- and transcriptomic-explained variances are further partitioned by chromosome and genome sections to offer a view of the genetic architecture on different aggregation levels.

The choice of Bayesian variable selection (BVS) models was due to its features to separate markers with large/moderate or small effects, and to locate the important regions in the genome or transcriptome which serves a better QTL mapping method because it produces clearer signals for QTL[11]. Furthermore the prediction based on genomic variables using BVS is more accurate even when the prior is not correct[11–14]. It is important to say that simpler methods suffer from “missing heritability” too[15, 16].

The aim of this study was to explore the contributions of various sources of variation, such as population structure, SNP variants, and gene expression levels, to a set of growth related traits (body weight, feed intake, and feed efficiency) in mice. These traits are very important, both in terms of agricultural production and for obesity in humans. Bayesian mixed models and Bayesian variable selection models were applied to model pedigree, SNPs and/or gene expressions and to derive explained variances for these components. In addition, they were used to partition of SNPs and gene expression by chromosome and genome sections. To validate the estimates of explained variances, the predictive ability of these models was studied using cross validation.

where X is the design matrix for batch and sex effects, Z is a design matrix that links polygenic effects to the observed records, W is a matrix with 1806 SNP covariates, and Q is a matrix with 23,698 gene expression covariates. The SNP and gene expression covariates were centered and scaled to unit variance.

Based on work of[19–22], the Bayesian mixed model version assigns normal prior to the vectors u, a, g, and e in (1), i.e.,$u\sim N\left(0,A{{\displaystyle \sigma }}_u^2\right),a\sim N\left(0,I{{\displaystyle \sigma }}_s^2\right),g\sim N\left(0,I{{\displaystyle \sigma }}_g^2\right),e\sim N\left(0,I{{\displaystyle \sigma }}_s^2\right)$, where${{\displaystyle \sigma }}_e^2$ is the polygenic variance and A is the numerator relationship matrix based on pedigree information,${{\displaystyle \sigma }}_s^2$ is the per-SNP explained variance,${{\displaystyle \sigma }}_g^2$ is the per-gene expression explained variance, and${{\displaystyle \sigma }}_e^2$ is the residual or environmental variance. These four variances are estimated in the model using flat prior distributions, i.e.,${{\displaystyle \sigma }}_u^2,{{\displaystyle \sigma }}_s^2,{{\displaystyle \sigma }}_g^2,{{\displaystyle \sigma }}_e^2\sim \mathit{\text{Bern}}$. The remaining parameters in (1), μ and b, are assigned flat prior distributions, which is the Bayesian analog of fitting “fixed effects” (unshrunken) estimates. A Markov chain Monte Carlo (MCMC) algorithm was applied in the software bayz[23] to obtain samples from the posterior distribution of the model parameters$f\left(\mu ,b,u,a,g,{{\displaystyle \sigma }}_u^2,{{\displaystyle \sigma }}_s^2,{{\displaystyle \sigma }}_g^2,{{\displaystyle \sigma }}_e^2|y\right)$. MCMC algorithms for sampling effects and variances in mixed models have been extensively described, for a general overview see[24]. The Monte Carlo accuracy of the MCMC algorithm was evaluated by correlating repeated estimates for the parameter vectors u, a and g, requiring a correlation >0.999 from repeated MCMC runs, and by computing the effective sample sizes for the variance components using the R Coda package[25].

The explained variance in y from (1) is var(Zu) + var(Wa) + var(Qg) + var(e). To obtain posterior means (PMs) and posterior standard deviations (PSDs) on the explained variances for SNPs and gene expressions, var(Wa) and var(Qg) were evaluated based on the posterior samples for a and g from the MCMC, i.e., as the PM and PSD of var(Wa^t) values over MCMC cycles, where a^t is the posterior sample for a from MCMC cycle t. This procedure is not required for the polygenic variance, because Z is a design matrix, unlike W and Q, which are covariate matrices.

where${\tau }_{a1}^2$ and${\tau }_{a0}^2$ are the “large” and “small” variances in the mixture distribution for a,${\tau }_{g1}^2$ and${\tau }_{g0}^2$ are the “large” and “small” variances in the mixture distribution for g, and${\gamma }_a$ and${\gamma }_g$ are vectors of 0/1 indicator variables for a and g, respectively, indicating whether the i th element in a or g, respectively, comes from the distribution with large or small variance. The variances${\tau }_{a1}^2,{\tau }_{a0}^2,{\tau }_{g1}^2,{\tau }_{g0}^2$ were all estimated from the data using unbounded flat prior distributions. The constraints${\tau }_{a1}^2>{\tau }_{a0}^2$ and${\tau }_{g1}^2>{\tau }_{g0}^2$ were applied using a rejection sampler, so that “large” and “small” effects remained identifiable. The priors for the indicator variables were taken as${\gamma }_{\mathit{ai}}\sim \mathit{\text{Bern}}\left({\pi }_a\right)$ and${\gamma }_{\mathit{gi}}\sim \mathit{\text{Bern}}\left({\pi }_g\right)$, where$\mathit{\text{Bern}}\left(\pi \right)$means a Bernoulli distribution for a 0/1 indicator with a probability π for a 1. The parameters${\pi }_a,{\pi }_g$were taken as known. The MCMC implementation of this model is relatively straightforward, because conditional on the indicator variables the model remains a mixed model. The updating of the mixture indicators is described in[26]. This model is also run in the software bayz[23], and the Monte Carlo accuracy was evaluated in the same way as the mixed model version.

From the posterior samples for a and g in the variable selection model, explained variances were computed and partitioned by chromosome and by genome section. The variable selection model is more suited to make such a partitioning, because unlike the mixed model version, it allows for different variance contributions per SNP. The explained variances were evaluated in the same way as for the mixed model, by evaluating var(Wa^t) and var(Qg^t) over MCMC cycles t, except that the a and g samples are obtained under the mixture model prior assumptions. The same expressions can be straightforwardly evaluated for parts of the SNPs or gene expressions to obtain explained variances per chromosome and for small windows of SNPs within chromosomes. Variance within a chromosome was computed using a 5-SNP sliding window to obtain a genomic variance profile.

It is difficult to choose an optimal windows size as it depends on extend of LD, marker density and an arbitrary cut-off for what is considered important LD. In the data analyzed here, average R² between adjacent SNPs was 0.55, and average R² between SNPs two apart was 0.39, which we considered sufficiently high to warrant computation of variances in a 5-SNP window. To study the relative importance of family structure, SNPs, and gene expressions, six sub models and the complete model (1) were used. These were models that use only pedigree information (PED), only SNP data (SNP), only gene expression data (GEX), SNP + GEX, PED + GEX, PED + SNP, and the complete model PED + SNP + GEX. These models always included sex and batch effects.

The predictive ability of the models was evaluated using an 11-fold cross-validation. For body weight, 440 records were divided randomly in 11 groups, each with 40 individuals. Feed intake and feed efficiency, with 337 records in total, were randomly divided in 10 groups of 30 records and one group of 37 records. The complete model, including all variance parameters, was re-estimated on each set of 10 folds and predictions were computed for the phenotypes in the remaining 11^th fold. All predictions from the 11-fold cross validation were collected to compute correlations between predicted and actual phenotypes, and regressions of predicted phenotypes on actual phenotypes, using the whole data set. The slope of the regression lines of predicted phenotypes on actual phenotypes are expected to be 1 if the model produces unbiased predictions, which would validate the estimates of explained variances. The University of Nebraska Institutional Animal Care and Use Committee approved all procedures and protocols.

Overall, explained variances increase by adding gene expression information (GEX; data from liver), i.e., in the most complete model (PED + SNP + GEX) explained variances were 88%, 75%, and 71% for BW, FI, and FE respectively. This confirms the assumption that gene expressions can explain a larger part of phenotypic variance than genetic or genomic information, by capturing environmental, and possibly non-additive, genetic effects through the gene expressions[5, 30]. Information on the genetic architecture of these traits is best judged from the relative contributions of genomic and transcriptomic data in the SNP + GEX model.

This model shows that, for these traits, the liver transcriptome contributes a larger portion of explained variance. This is most pronounced for BW, with 18% of explained variance from the genome and 63% from the liver transcriptome. Thus, in this case, the predominant model is that SNPs regulate gene expressions to exert their effect on the phenotype.

Figure1 shows a decomposition of the explained variances at the chromosome and sub-chromosome level for the models using genomic (SNP) and genomic with transcriptomic (SNP + GEX) data for the trait BW. For the traits FI and FE, see Additional file2 and Additional file3 respectively. These results are based on the Bayesian variable selection model to better differentiate between genomic regions contributing more and less variance. The genomic variances at the sub-chromosome level are explained variances in a sliding 5-SNP window. At the chromosome level, chromosome 10 particularly stands out, with a relatively large contribution from the SNPs effects via transcriptome, but only a small contribution from the genome alone in explaining the phenotype. This does not mean there is no important QTL on this chromosome. In fact, there is a large QTL on chromosome 10; however, it is an eQTL whose effect can be captured by gene expressions. Figures 1b and c show the details at the sub-chromosome level, with Figure1b showing the explained genomic variances when fitting SNPs alone (SNP model), and Figure1c showing the explained genomic variances when adding gene expressions to the model (SNP + GEX model). The differences between these two graphs show locations of QTLs that regulate the liver transcripts and QTLs that exert their effect on the phenotype through another route. For BW, the main QTLs are found on chromosomes 1, 2, 9, 10, and 11, and the QTLs on 1, 9, and 10 appear to be eQTLs affecting gene expression in the liver. The QTL on chromosome 2 is an intermediate case whose effect is reduced, but does not completely disappear, when adding gene expressions to the model. Thus, this chromosome 2 QTL regulates liver transcripts, but must also have effects on BW through other routes, possibly by regulating genes outside the liver. The chromosome 11 QTL is a clear case of a QTL whose effect on BW does not work via the regulation of liver transcripts. The QTL locations are in agreement with QTLs detected for body weight in other studies[17, 31–33]. The same graphs for traits FI and FE are provided as supplementary material. These traits show relatively more cases where QTL effects remain after adding liver transcriptome data, which is in agreement with results in Table 1.

This method/approach is suitable for gene-level resolution. However, gene-level resolution is highly data dependent, i.e. it requires high marker density and a study population with LD blocks that span small genomic regions. In this work we have used F2 crosses from outbred lines, which has large LD blocks and this kind of data has limited resolution for fine-mapping of QTL.

One may argue that the most complete model is more interesting to investigate genetic architecture and chromosomal/sub-chromosomal variance but as we have shown SNPs and pedigree are largely confounded and they explain about the same variance. This confounded explained variance is getting worse in the case that both Pedigree and SNPs are in one model (PED + SNP model) which is shown in higher confidence intervals of explained variance by pedigree. The model with only omics information (SNP + GEX) is therefore simpler, more accurate and as effective as the model that also uses pedigree information. This is interesting for future applications of omics technologies, because we expect that pedigree information often will be absent.

Figures 2 and3 present detailed graphs of the genomic variances (left panels) and the distribution of chromosomal genomic values or breeding values[34] of the animals (right panels) for chromosomes 9 and 11, and for models fitting SNP only (top) or SNP + GEX (bottom). Breeding value is defined as the value of an individual as a parent based on sum of its genes effects[34]. Chromosome 9 is the case of an apparent eQTL, showing that genomic variance disappears, and that a tri-modal distribution of genomic values collapses, when gene expressions are added to the model. Chromosome 11 is the case of a QTL that does not regulate liver transcripts. The detailed picture of chromosome 11 shows that adding gene expressions to the model makes the effects of this QTL clearer: genomic variances outside the QTL region reduce, and a clear tri-modal distribution of chromosomal genomic values is seen in the SNP + GEX model, but not in the SNP-only model. Table 2 shows the rank correlations between the predicted values from using pedigree (PED), genomic (SNP), or transcriptomic (GEX) information. Pedigree and genomic values correlate better than pedigree/genomic values with transcriptomic values. This confirms that pedigree and genomic information overlap to a reasonable degree, but this is less true for transcriptomic information.

With increased availability of various -omics data, integrative approaches are promising tools for understanding the genetic architecture of complex traits. We have developed a complementary approach to the univariate “eQTL” mapping, by considering Bayesian models that fit all genome-wide SNPs and transcript abundances in one model, and that estimate and partition explained variances by chromosome and genome segments. Our results show that, using gene expressions, more of the phenotypic variance can be explained and phenotypes can be better predicted. Predictions were also shown to be unbiased, which validates the assessed explained variances. The improvement of phenotype predictions using gene expression data will be useful for several applications in agriculture and medicine, although it should be assessed on a case-by-case basis as to whether a suitable tissue can be sampled for the gene expression measurements. Partitioning of the explained genomic variance at the level of chromosomes and genome segments showed clear examples of eQTL locations as regions where genomic variance disappears when gene expressions are added to the model. Our study used only gene expressions from the liver, and an obvious further extension is to include expressions from other tissues. The QTLs that did not disappear when transcripts are added to the model may be eQTLs that affect gene expression in a tissue other than liver. The Bayesian model is quite efficient for handling large sets of covariates, and extensions to include multiple sets of expressions will be feasible. We have not provided formal statistical tests in this model, but the Bayesian approach lends itself naturally to obtaining confidence intervals for (differences between) parameter estimates. The estimates of total explained variances from the Bayesian mixed model can also be obtained by a residual maximum likelihood (REML) approach. We verified this, and the Bayesian and REML estimates generally agree. However, using REML it is not feasible to utilize mixture priors to better discriminate between SNPs which contribute more or less variance, and to partition the variances at the sub-chromosome level, which is all straightforward in a Bayesian approach.

Our approach can easily allow up scaling to higher-density arrays, even to whole-genome sequence data with the variance components analysis as it was for gene expression probes in this study.