Model specification

We consider an mvLMM [5] to model two phenotypes of interest jointly. These two phenotypes could be collected from a single study or from two studies with either non-overlapping or partially overlapping individuals. We denote y₁ and y₂ as the two phenotype vectors, measured on n₁ and n₂ individuals, respectively. We denote n_s as the number of overlapping individuals who have both phenotypes measured. In the case of two separate studies, the sample size n₁ generally does not equal to n₂ and n_s = 0. In the case of a single study with completely overlapping individuals, the sample size n₁ = n₂ = n_s. Besides phenotypic measurements, all individuals have their genotypes measured on a common set of m SNPs. We denote the n₁ by m matrix X₁ as the genotype matrix for the n₁ individuals measured with the first phenotype and denote the n₂ by m matrix X₂ as the genotype matrix for the n₂ individuals measured with the second phenotype. To facilitate computation, we center and standardize each phenotype vector as well as each column of the genotype matrices to have a mean of zero and a standard deviation of one. We consider the following regression model to related phenotypes to genotypes (1) where β₁ and β₂ are m-vectors of genotype effect sizes for the two phenotypes, respectively; and ϵ₁ and ϵ₂ are vectors of environmental effects/residual errors, with dimension n₁ and n₂, respectively.

We follow the standard mvLMM assumption and assume that the genetic effect sizes of j’th SNP on the two phenotypes follow a multivariate normal (MVN) distribution a priori (2) where represents the heritability of the first phenotype; represents the heritability of the second phenotype; and ρ_g represents the genetic covariance, which characterizes the phenotypic covariance explained by genetic effects. We further define the genetic correlation as .

For the n_s individuals who have both phenotypes measured, we assume that their environmental effects for the i’th individual follow (3) where ρ_e represents the environmental covariance, which characterizes the phenotypic covariance explained by environmental effects. The environmental correlation is further defined as For the remaining individuals who only has one phenotype measured, we assume that ϵ_1i for the i'th individual follows a normal distribution N(0, ) while ϵ_2j for the j’th individual follows another normal distribution N(0, ).

We develop a composite likelihood based algorithm to perform estimation and inference for the mvLMM model defined in Eqs (1)–(3). The composite likelihood based algorithm requires only summary statistics as input and is both computationally and statistically efficient. For input summary statistics, we obtain the marginal z-score for each SNP-phenotype pair by fitting a univariate linear regression model. Specifically, the marginal z-score of j’th SNP for d’th phenotype (d = 1, 2) is , where X_dj represents j’th column of the genotype matrix X_d. We derive the marginal distribution for the two z-scores of each j’th SNP based on mvLMM. Such marginal distribution is in the following form (details in the S1 Text): (4)

Here ρ = ρ_e+ρ_g; is the LD score for the j’th SNP computed based on a reference panel, with r_ij being the correlation coefficient between j’th SNP and i’th SNP.

We denote as the likelihood obtained based on Eq (4). The marginal z-scores for different SNPs are correlated with each other due to linkage disequilibrium (LD). Consequently, the joint likelihood for genome-wide marginal z-scores are in a complicated form, on which parameter estimation is computationally challenging to carry out. To enable estimation with summary statistics, we approximate the complicated joint likelihood with a relatively simple composite likelihood, which is represented as a weighted product of individual marginal likelihood across genome-wide SNPs [13]: (5) where each individual likelihood is powered to w_j, a j’th SNP specific weight. The weight w_j is critical for ensuring estimation accuracy of the composite likelihood based algorithm in the presence of LD [14]. Intuitively, if j’th SNP is in LD with many other SNPs, then the marginal likelihood of (z_1j, z_2j) contains similar information as the marginal likelihood of the other SNPs that are in LD with the j’th SNP. Consequently, we can choose a small w_j to down-weight the j’th marginal likelihood. In contrast, if the j’th SNP is in LD with only a limited number of SNPs, then the marginal likelihood of (z_1j, z_2j) contains important information that cannot be replaced by that of the other SNPs. Consequently, we want to choose a large w_j to up-weight the j’th marginal likelihood. To captivate such intuition, we set w_j = 1/l_j with l_j being the LD score of the j’th SNP. Such choice of w_j is optimum for achieving maximum estimation accuracy of composite likelihood for the specific setting where the SNP correlation matrix is block diagonal [14]: SNPs within blocks have identical genotypes with correlations among them being one, while SNPs between blocks have zero correlations. In this setting, it has been shown that the optimal weighting choice is w_j = 1/(m_i−1), where m_i is the number of SNPs in i’th block to which the j’th SNP belongs [14]. This optimal weight choice w_j = 1/(m_i−1) is equivalent to w_j = 1/l_j, the weighting choice we make for our method. In addition, the composite likelihood under our weighting choice of w_j = 1/l_j is equivalent to the full likelihood when SNPs are uncorrelated with each other. The weight choice of w_j = 1/l_j also relates our method to LDSC and MQS in the single phenotype setting (details in S1 Text) [12]. Therefore, our method can also be viewed as a natural extension of LDSC and MQS towards modeling multiple phenotypes.

With the above composite likelihood, we develop an expectation maximization (EM) algorithm [15] paired with an Newton-Raphson (NR) algorithm [16,17] to perform parameter estimation. In addition, we use the jackknife algorithm [10,18] to estimate standard errors for the parameter estimates. The composite likelihood-based algorithm allows us to perform unbiased parameter estimation [13] and make inference using only summary statistics in a computationally scalable fashion. The composite likelihood-based algorithm can also be easily extended to handle multiple genetic components in the presence of SNP annotations (details in the S1 Text). We refer to our method as genetic and environmental covariance estimate by composite-likelihood optimization (GECKO). GECKO is implemented as an R package, freely available at www.xzlab.org/software.html.

Simulation settings

We used genotype data from the Welcome Trust Case Control Consortium (WTCCC) study and simulated phenotypes. Specifically, we obtained 384,733 SNPs from 2,938 control samples in the WTCCC. We filtered out SNPs with a minor allele frequency less than 5%, with a Hardy–Weinberg equilibrium p-value less than 0.001, or that are strand-ambiguous. We focused on the remaining 289,994 SNPs and simulated pairs of traits. In each simulation, we first simulated the genetic effect sizes for all SNPs based on Eq (2) with different genetic covariance values (more details below). We simulated the environmental effects based on Eq (3) with different environmental covariance values. We then summed the genetic effects and the environmental effects to obtain simulated phenotypes based on Eq (1). Afterwards, we computed marginal z-scores for each trait SNP pair using linear regression. We also computed LD scores using all individuals, with a window size set to be 1MB to minimize the influence of LD decay following the recommendation of LDSC [19]. Finally, we applied different methods to estimate the genetic and environmental covariance based on the resulting z-scores and LD scores. In the simulations, we examined different combinations of genetic and environmental covariance values, examined the influence of sample overlap on estimation accuracy, and examined the setting where genetic effects vary across different genetic annotations. In total, we examined four main simulation settings as described below:

Compared methods

We mainly compared our method with two existing summary statistic based methods: LDSC and GNOVA. Both these two methods estimate the genetic covariance through MoM and both of them rely on marginal z-scores and LD scores as input. In LDSC, we used the default setting for genetic covariance estimation. In GNOVA, we also used the default setting in the one study design and the two partially overlapped study design; and we followed the recommendation of GNOVA in the two separate study design to improve estimation performance. Neither GNOVA nor LDSC estimate can directly estimate or test environmental covariance. LDSC outputs an estimate for the intercept covariance, which can be approximately considered as the summation of the genetic covariance and the environmental covariance. Therefore, we post hoc processed the LDSC output, obtained the difference between the intercept covariance estimate and the genetic covariance estimate and treated it as an ad hoc estimate for the environmental covariance. However, we are unable to test for the environmental covariance using LDSC: LDSC does not output the standard error for the environmental covariance and it is not straightforward to extract such standard error in a post hoc fashion. For GNOVA, we also obtained the LDSC estimate of the LDSC intercept covariance, from which we subtracted the GNOVA genetic covariance estimate to obtain an ad hoc estimate for the environmental covariance. In stratified analysis, we obtained the GNOVA genetic covariance estimates through the above procedure in each annotation region separately. Afterwards, we obtained the environmental covariance estimate by taking the difference between the LDSC intercept estimate and the summation of GNOVA genetic covariance estimates across annotation regions. Besides the above summary statistics based methods, we also compared to the individual level data based method mvLMM. Specifically, we fitted mvLMM using GCTA and examined the simulation setting III with different sample overlap compositions. Finally, we note that we used 1/l_j as the weight for the composite likelihood in GECKO throughout the text as explained above, where l_j is the LD score of the j’th SNP.

Real data analysis

We applied our method to analyze GWAS summary statistics for 22 traits collected from five GWASs. These 22 traits include nine non-impedance traits from the UK biobank [21] that have SNP heritability greater than 0.05, four lipid traits from the Global Lipids Genetics Consortium (GLGC) [22], five birth traits from Early Growth Genetics (EGG) Consortium [23–28], two anthropometric measures from the Genetic Investigation of ANthropometric Traits (GIANT) Consortium [29,30], and two bone density measurements from the GEnetic Factors for OSteoporosis Consortium (GEFOS) [31]. Details of these traits are provided in the S1 Table. The 22 traits can be classified into four categories: anthropometric traits (e.g., height and BMI), hematological traits (e.g., WBC and RBC), birth traits (e.g., birth weight and birth length), and metabolic traits (e.g., TG and HDL). For analysis, we first obtained a common set of SNPs shared across all traits. We filtered out SNPs with minor allele frequency (MAF) less than 0.05, with Hardy-Weinberg equilibrium (HWE) p-value less than 0.001, or that are strand-ambiguous. We used the R package liftOver to match SNP base positions onto NCBI 37 build whenever necessary. We further overlapped SNPs with those available in the 1,000 Genomes project and focused on a final set of 611,444 SNPs for analysis. On the final set of SNPs, we calculated LD scores based on 503 individuals with European ancestry from the 1,000 Genomes project [32] with a window size of 1 MB as recommended by LDSC [19]. We applied GECKO together with LDSC and GNOVA to analyze all trait pairs among the 22 traits. We declared significance in covariance estimates based on Bonferroni corrected thresholds. For each trait pair in turn, we also performed gene set enrichment analysis using GSA-SNP2 [33] to identify enriched gene sets. We further compared the proportion of enriched gene sets in trait pairs with either significant or non-significant genetic covariances.

In addition, we examined the performance of different methods on denser genotypes in the real data. Specifically, we focused on the 13 phenotypes from GLGC (4) and UKBB (9) that have a much higher number of genotypes as compared to the remaining phenotypes. We obtained a common set of SNPs among these 13 phenotypes and retained SNPs that are also available in the 1,000 Genomes project reference panel. We then filtered out SNPs with minor allele frequency (MAF) less than 0.05, with Hardy-Weinberg equilibrium (HWE) p-value less than 0.001, or that are strand-ambiguous. After filtering, we obtained a final set of 1,731,189 SNPs for analysis.

Genetic and environmental covariance estimation

We performed simulations to compare the estimation accuracy of GECKO with two existing approaches, LDSC and GNOVA. Simulation and comparison details are provided in the Methods. Briefly, we used genotypes from WTCCC to simulate pairs of traits. We examined a total of 238 simulation scenarios in four main simulation settings. We first examined the one study design where both traits are measured on the same set of individuals (n = 2,938). We examined the estimation accuracy of the three methods for a range of genetic covariance values in the setting where the environmental covariance is either zero (Fig 1A and 1B), positive (S1A and S1B Fig), or negative (S1C and S1D Fig). All three methods provide approximately unbiased estimates of the genetic covariance regardless of the environmental covariance values (Figs 1A and S1A and S1C). To quantify estimation accuracy, we computed the mean squared error (MSE) for estimates obtained from different methods. To facilitate visualization, we also contrasted the MSE from each of the other two methods to GECKO by computing an MSE ratio (Figs 1B and S1B and S1D). The MSE ratio measures the inverse of the relative statistical efficiency of a given method with respect to GECKO: an MSE ratio above one suggests that GECKO provides more accurate estimates while an MSE ratio below one suggests the opposite. In addition, the MSE ratio directly quantifies how effective the two methods are in using the samples, as the estimation variance is a function of sample size due to the central limit theorem. For example, an MSE ratio of 1.11 would imply that GECKO effectively increases the sample size by 1.11−1 = 11% as compared to the other method. The MSE ratio results show that GECKO provides more accurate genetic covariance estimates than the other two methods across a wide range of scenarios. The performance of GECKO is often followed by LDSC. For example, when the environmental covariance is zero, the MSE ratio of LDSC ranges from 0.92 to 1.39 with a mean of 1.11 (with the ratio in 70.6% of simulation replicates being above 1), while the MSE ratio of GNOVA ranges from 1.01 to 1.87 with a mean of 1.22. The overall results do not vary much with respect to the environmental covariance values. For example, when the environmental covariance is 0.1, the MSE ratio of LDSC ranges from 0.86 to 1.48 with a mean of 1.16 (with 82.3% above 1), while the MSE ratio of GNOVA ranges from 1.04 to 1.69 with a mean of 1.26. Similarly, when the environmental covariance is -0.1, the MSE ratio of LDSC ranges from 0.95 to 1.59 with a mean of 1.13 (with 76.5% above 1), while the MSE ratio of GNOVA ranges from 1.05 to 1.65 with a mean of 1.25. The estimation accuracy of the genetic correlation with varying environmental correlation showed a similar overall pattern as the results on genetic covariance estimation (S2A, S2B and S3A–S3D Figs). Besides genetic covariance estimation, we also examined the estimation accuracy of environmental covariance with varying genetic covariances. Because LDSC and GNOVA do not provide environmental covariance estimates, we had to post-process LDSC and GNOVA output to obtain these estimates in an ad hoc fashion (details in Methods). The results on environmental covariance estimation are largely consistent with genetic covariance estimation results. For example, all three methods provide reasonably unbiased environmental covariance estimates regardless of the genetic covariance values (Figs 1G and S1E and S1G). GECKO produces more accurate environmental covariance estimates than both LDSC and GNOVA in most (78.4%) simulation scenarios (Figs 1H and S1F and S1H).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Comparison of genetic and environmental covariance estimation of different methods in simulations.

Compared methods include GECKO (purple), GNOVA (grey), and LDSC (green). Results are shown for the one study design (first row: A, B, G, H), two partially overlapped study design (second row: C, D, I, J), and two separate study design (third row: E, F). Boxplots display estimated genetic covariances (A, C, E) and environmental covariances (G, I) on y-axis versus the true covariances on x-axis across simulation replicates. Estimation accuracy is measured by the ratio of mean square errors (MSE), which contrast the MSE from GNOVA or LDSC with respect to GECKO, across various true covariances on x-axis, for genetic (B, D, F) and environmental covariances (H, J). An MSE ratio below one suggests that GECKO performs worse than the other method; above one otherwise.

https://doi.org/10.1371/journal.pgen.1009293.g001

Examining the detailed scenarios in which different methods perform well yields further insights. Specifically, LDSC is particularly effective in the scenario where the genetic and environmental covariance sums to zero. For example, when the environmental covariance is zero, the estimate by LDSC is accurate when the true genetic covariance is close to zero, even slightly more so than GECKO. However, the performance of LDSC degrades quickly when the summation of genetic and environmental covariance differs far away from zero. The performance dependence of LDSC with respect to the summation of the true covariances closely resembles its performance dependence on the true heritability for heritability estimation–LDSC heritability estimate is also more accurate when the true heritability is close to zero [12]. In contrast, the performance of GNOVA is almost in the opposite direction of LDSC: the estimation accuracy of GNOVA is the worst when the genetic and environmental covariance sum to zero and improves when the sum deviates away from zero. The different performance dependence on the true covariance values in GNOVA and LDSC is likely resulted from different MoM type algorithms used in the two methods [12]. Different from LDSC and GNOVA, the performance of GECKO is reasonably stable across a range of genetic and environmental covariance values (Figs 1B and S1B and S1D), supporting the benefits of performing inference in a composite likelihood framework.

Sample overlap, type I error control, and stratified analysis

Besides the above settings where both traits are measured on the same set of individuals, we also examined settings where the two traits are measured on two studies, with either no overlapping or partially overlapping individuals between the two. Specifically, we set the first study sample size n₁ = 1,500 and the second study sample size n₂ = 1,438 in the two separate study design. We also set n₁ = 2,000, n₂ = 1,938, with overlap sample size n_s = 1,000 in the two partially overlapped study design. The results are largely consistent with the one study design described in the previous section. For example, all methods provide approximately unbiased estimates for both the genetic covariance (Fig 1C and 1E) and the environmental covariance (Fig 1G and 1I). Compared to the other two methods, GECKO produces more accurate estimates of the genetic covariance in most (82.3%) scenarios (Fig 1D and 1F) and more accurate estimates of the environmental covariance in most (70.6%) scenarios (Fig 1H and 1J). As expected, the proportion of individual overlap between studies positively influences the estimation accuracy of all methods. For example, the average MSE of the genetic covariance estimates by GECKO increases from 0.012 in the one study design to 0.026 in the two partially overlapped study design, and further to 0.043 in the two separate study design. Certainly, all these compared methods are based on summary statistics and use either method of moments or composite likelihood. Consequently, these methods produce less accurate genetic and environmental covariance estimates as compared to the full likelihood based method mvLMM, which uses individual level data (S4 Fig). For example, the MSE ratio of mvLMM over GECKO for genetic covariance estimation is much smaller than one in the presence of sample overlap though it becomes closer to one in the two separate study design. Overall, the results suggest that the performance of GECKO is robust with respect to sample compositions.

Besides estimation, we found that both LDSC and GECKO control for type I error well for genetic covariance testing and produce conservative p-values across study designs (Table 1). Consistent with [11], GNOVA also produces reasonably calibrated type I error control in the two separate study design. However, GNOVA produces inflated type I error in both one study design and two overlapping study design, sometimes several times larger than expected (Table 1); both these two study designs were not explored in [11]. In addition, GECKO is the only method capable of testing the environment covariance, producing controlled type I error and conservative p-values across settings (Fig 2D and 2E and Table 1). Because different methods have different control of type I error, we compared the power of different methods at a fixed type I error rate instead of a nominal p-value threshold. In the analysis, we found that GECKO is also more powerful than the other two methods in testing genetic covariance, despite the relatively small sample sizes used in simulations. For example, in the one study design, GECKO improves power by 4.05% and 32.4% on top of LDSC and GNOVA, respectively (Fig 2A). In the partially overlapped study design, GECKO improves power by 3.72% and 9.44% on top of LDSC and GNOVA, respectively (Fig 2C). GNOVA is even more powerful than the other two methods in the two separate study design (Fig 2E).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Power comparison of different methods in detecting non-zero genetic and environmental covariances in simulations.

Compared methods include GECKO (purple, solid line), GNOVA (grey, dotted line), and LDSC (green, dashed line). Results are shown for the one study design (first row: A, D), two partially overlapped study design (second row: B, E) and two separate study design (third row: C). Power (y-axis) of different methods are shown with respect to the true covariances (x-axis) for detecting non-zero genetic covariances (A, B, C) and environmental covariances (D, E). Power is shown based on a type I error of 0.05 but not a nominal p-value of 0.05. The power of GNOVA and LDSC for detecting non-zero environmental covariance are not shown because GNOVA and LDSC cannot test for environmental covariance.

https://doi.org/10.1371/journal.pgen.1009293.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Type I error control of different methods under null simulations.

https://doi.org/10.1371/journal.pgen.1009293.t001

Next, we examined the performance of GECKO in estimating annotation-stratified genetic covariance and environmental covariance. To do so, we partitioned the whole genome into two non-overlapping annotation categories and applied both GECKO and GNOVA to estimate the two annotation-stratified genetic covariances. We did not include LDSC here as it cannot estimate annotation-stratified covariances. The results in the annotation-stratified analysis are largely similar to the results for non-stratified analysis. For example, both GECKO and GNOVA provided approximately unbiased estimates of the partitioned genetic covariances regardless of the study design (Fig 3A, 3C and 3E). GECKO also produces more accurate genetic covariance estimates than GNOVA. For example, in the one study design, the MSE ratio of GNOVA ranges from 1.01 to 1.08 with a mean of 1.05 (Fig 3B). The relative accuracy of GECKO with respect to GNOVA improves with reduced sample overlap. For example, the MSE ratio of GNOVA to GECKO ranges from 1.10 to 1.23 with a mean of 1.17 in the two separate study design (Fig 3D), and ranges from 1.16 to 1.44 with a mean of 1.36 in the two partially overlapping study design (Fig 3F). In terms of environmental covariance, GNOVA cannot directly perform estimation in the presence of stratified annotations so we had to use a post hoc procedure to obtain environmental covariance estimates from GNOVA. Consistent with results on the unstratified environmental covariance estimation, GECKO provided approximately unbiased estimate of environmental covariance regardless of study design (Fig 3G and 3I) and is more accurate than GNOVA (Fig 3H and 3J).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Comparison of genetic and environmental covariance estimation in presence of stratified genetic components for different methods in simulations.

Compared methods include GECKO (purple) and GNOVA (grey). Results are obtained in the presence of stratified genetic components. Results are shown for the one study design (first row: A, B, G, H), two partially overlapped study design (second row: C, D, I, J) and two separate study design (third row: E, F). Boxplots display the estimated total genetic covariances (A, C, E) and environmental covariances (G, I) on y-axis versus the true covariances on x-axis across simulation replicates. Estimation accuracy is measured by the ratio of mean square errors (MSE), which contrast the MSE from GNOVA or LDSC with respect to GECKO, across various true covariances on x-axis, for genetic (B, D, F) and environmental covariances (H, J). An MSE ratio below one suggests that GECKO performs worse than the other method; above one otherwise.

https://doi.org/10.1371/journal.pgen.1009293.g003

Besides estimation, we found that GECKO is more powerful than GNOVA in testing genetic covariance in the one study design and the two partially overlapped study design. For example, in the one study design, GECKO improves power by 17.94% on top of GNOVA (S5A Fig). In the partially overlapped study design, GECKO improves power by 2.59% on top of GNOVA (S5B Fig). However, GNOVA is more powerful than GECKO in the two separate study design (S5C Fig), presumably due to its implicit modeling assumption that the environmental covariance is zero, which happens to be true in the two separate study design.

Alternative simulations and model misspecifications

Besides the main simulations, we examined the performance of GECKO and the other methods with additional simulation settings including those under different model misspecifications. We first examined the performance of different methods under simulations with denser genotypes. Specifically, we used imputed genotypes and performed simulations under setting III. We found that the results under denser genotypes are largely consistent with the above main results (S6 Fig). We then examined the performance of different methods based on a moderate heritability of 0.15 instead of 0.5. A heritability of 0.15 is close to the mean heritability estimate obtained in our real data application. The results with moderate heritability are also largely consistent with the main results (S7 Fig). We further examined a LD mismatch setting, where we calculated LD scores and subsequent composite likelihood weights based on the 1,000 Genomes project reference panel instead of directly from the data at hand. The results are largely consistent with the main simulations: all methods provide approximately unbiased estimates for both the genetic (S8A, S8C and S8E Fig) and environmental covariances (S8G and S8I Fig) while GECKO produces more accurate estimates of the genetic covariance (S8B, S8D and S8F Fig) and the environmental covariance in most (82.3% and 76.5%) scenarios (S8H and S8J Fig). We also examined another model misspecification setting where the number of overlapping individuals is mis-specified. Because the genetic covariance and the number of overlapping individuals are two parameters that are separable from each other in the model, the estimation of genetic covariance from different methods would not be affected by the misspecification of the number of overlapping individuals (S9A, S9C, S9G, and S9I Fig). Indeed, the MSE ratio comparing different methods remains the same (S9B, S9D, S9H and S9J vs S9F Fig). However, the environmental covariance and the number of overlapped individuals are two parameters that are coupled together in the model and are thus not identifiable from each other. Consequently, when the number of overlapping individuals is mis-specified, the environmental covariance estimate would become biased for all three methods (S10A, S10C, S10G and S10I Fig), although the MSE ratio comparing different methods remains largely the same (S10B, S10D, S10H and S10J vs S10F Fig).

Finally, we note that the computational complexity of GECKO scales linearly with the number of individuals and with the number of SNPs. While GECKO requires iterative optimization, the linear computational complexity of GECKO makes it reasonably efficient as compared to MoM based approaches. Indeed, across simulation settings, the average computing time for covariance estimation is 95.43s, 7.41s, and 7.67s for GECKO, LDSC and GNOVA, respectively.

Real data applications

We applied GECKO and LDSC to analyze 22 quantitative traits from five different GWASs. These traits belong to five distinct phenotype groups that include birth traits, lipid traits, bone density traits, blood traits, and anthropometric measures. The heritability estimates of 22 traits range from 0.047 to 0.368 with a mean heritability estimate being 0.155. We focused on all 231 pairs of traits and a common set of 611,444 SNPs for analysis. We did not include GNOVA in comparison as it does not provide a well-controlled type I error for testing. Across pairs of traits, the genetic correlation estimates from GECKO range from -0.731 to 0.945, and the environmental correlation estimates from GECKO range from -0.363 to 0.875. In the analysis, GECKO identified 50 significant genetic correlations based on Bonferroni correction (p-value < 0.00022 (0.05/231); Fig 4A). Among them, a higher proportion of significant genetic correlation are obtained among trait pairs collected in the same study (19 out of 54; 35.19%) than that collected in separate studies (31 out of 177; 17.51%), consistent with the fact that traits are collected in the same study if they likely share a common genetic background. Consistent with the lower power of LDSC in simulations, LDSC identified 24 significant genetic correlations, the majority of which (19) were also identified by GECKO (S11 Fig) and 19 of which were obtained among trait pairs collected in the same study. Among the 19 pairs of traits with significant genetic correlation identified by both methods, the majority of them have previous literature support (S2 Table). For example, both methods identified positive genetic correlation between LSBMD and FNBMD, and between BMI and TG, as well as negative genetic correlation between BMI and HDL, and between TG and HDL [10,11]. In terms of the 31 significant genetic correlation obtained only by GECKO, many of them also have previous literature support. For example, the positive genetic correlation between BMI and WBC (correlation = 0.091; p-value = 2.24×10⁻¹²) is consistent with the previous study that common inflammatory pathways may underlie both traits: high BMI and obesity is considered as a chronic inflammatory condition with which elevated WBC is widely recognized to be associated with [34,35]. As another example, the positive genetic correlation between BMI and SBP (correlation = 0.06; p-value = 1.46×10⁻⁶) suggests that either common genetic pathways may underlie both BMI and SBP or that a potentially causal positive effect of BMI on hypertension exist as suggested previously [36]. Indeed, gene set enrichment analysis (GSEA) identified thyroid hormone receptor binding pathway (WBC q-value = 0.00079; BMI q-value = 0.033) and cytokine-mediated signaling pathway (WBC q-value = 0.0014; BMI q-value = 0.013) as enriched in genes associated with BMI or SBP [37–39]. Similarly, the positive genetic correlation between BMI and LDL (correlation = 0.23; p-value = 4.14×10⁻⁹) suggests that either common genetic pathways underlie both BMI and LDL or BMI has a positive and potentially causal effect on LDL as suggested previously [40,41]. Indeed, GSEA identified lipid storage pathway (BMI q-value = 0.022; LDL q-value = 1.29×10⁻⁷) and receptor biosynthetic process pathway (BMI q-value = 0.0067; LDL q-value = 0.0001) as enriched in genes associated with BMI or LDL, supporting common genetic regulation underlying the two traits. Besides positive genetic correlation, the identified negative genetic correlation between pairs of traits are also consistent with previous literature. For example, the negative genetic correlation identified between height and LDL (correlation = -0.16; p-value = 2.96×10⁻¹³) is consistent with previous studies that height associated SNPs display negative associations with LDL [42]. In addition, GSEA analysis revealed that, consistent with previous studies [43], the steroid biosynthetic process is enriched in genes associated with both height and LDL (LDL q-value = 0.0002; height q-value = 0.038). Similarly, the negative genetic correlation between WBC and height (correlation = -0.065, p-value = 1.60×10⁻¹²) is also supported by GSEA analysis, where we found that the receptor biosynthetic process (WBC q-value = 0.01, height q-value = 0.01) and the multicellular organism growth pathways (WBC q = 0.004, height q-value = 0.009) are both enriched in genes associated with either of the two traits [44]. Overall, the proportion of enriched pathways identified in trait pairs with a significant genetic correlation is 4.4 fold higher than that in trait pairs with non-significant genetic correlation (7.0% vs 1.6%; p-value < 1×10⁻⁴; Fig 4B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Genetic and environmental correlation estimation by GECKO across pairs of 22 GWAS traits.

A: heatmap displays genetic (lower triangle) and environmental (upper triangle) correlation estimates for pairs of traits from GECKO. A cross in the square represents statistically significant correlation estimates for the trait pair after Bonferroni correction. Environmental correlation estimation is carried out only for pairs of traits that are collected in the same study (non-grey boxes). B: violin plot displays the proportion of significantly enriched gene sets detected for pairs of traits with non-significant genetic correlation estimates (purple) and for pairs of traits with significant genetic correlation estimates (green). C: scatterplot contrasts the genetic correlation estimates (y-axis) versus the environmental correlation estimates (x-axis) for pairs of traits that are collected in the same study. Traits are organized into five different phenotype categories (figure legend).

https://doi.org/10.1371/journal.pgen.1009293.g004

Besides genetic correlation, we also applied GECKO to obtain environmental correlation estimates among 54 trait pairs in which both traits were measured on a common set of individuals. We did not apply LDSC for environmental correlation estimation becaues it cannot test the environmental correlation. The estimated environmental correlation estimates are positively correlated with the genetic correlation estimates across trait pairs (Spearman correlation = 0.707; Fig 4C), with the two estimates sharing the same sign for majority of trait pairs (45 out of 54; 83.3%). The positive correlation between the environmental correlation estimates and genetic correlation estimates is highest for pairs of traits in the lipid phenotype group (correlation = 0.98) where 100% of the trait pairs share the same sign, and is the lowest for pairs of traits in the blood phenotype group (correlation = 0.25) where 66.7% of the trait pairs share the same sign. 9 trait pairs (out of 54 pairs) showed opposite signs of the genetic and environmental covariance estimates, though none of these pairs were statistically significant. Among the 54 analyzed trait pairs, GECKO identifies 24 significant environmental correlation based on Bonferroni correction (p-value < 0.0009 (0.05/54); Fig 4A). Among them, 11 pairs have significant genetic correlation estimates with the same sign. The significant environmental correlation estimates are identified among pairs of traits within each of the three phenotypic categories that include birth traits, metabolic traits, and anthropometric traits (S3 Table). For example, significant positive environmental correlation is observed between birth length and birth weight (correlation = 0.31; p-value = 4.24×10⁻³⁷), between birth weight and infant head circumference (correlation = 0.18; p-value = 1.25×10⁻⁵), and between childhood BMI and pubertal growth (correlation = 0.19; p-value = 3.49×10⁻³). Such positive environmental correlation estimates suggest that common environmental factors, such as maternal health status and maternal dietary, may influence birth traits in a coherent fashion [45–48]. As a second example, the significantly positive environmental corelation estimates between bone mineral density at femoral neck (FNBMD) and lumbar spine (LSBMD; correlation = 0.52; p-value = 3.02×10⁻⁶³) suggest that environmental factors such as physical activity may influence bone mineral density at two different anatomical locations in a similar fashion [49]. As a third example, the significantly negative environmental correlation estimate between height and BMI (correlation = -0.09, p-value = 4.93×10⁻²⁵) suggests that environmental factors such as television watching, physical activity, and smoking habits may influence height and BMI in different directions [50–52]. Finally, paring the genetic correlation estimates with the environmental correlation estimates helps us better understand the main contribution of phenotypic correlation between pairs of traits. For example, previous studies have identified a positive phenotypic association between WBC and pulse rate [53]. Such positive phenotypic association appeared to be accounted for by a positive environmental correlation between the traits (genetic correlation = 0.097, p-value = 2.28×10⁻³; while environmental correlation = 0.16, p-value = 6.90×10⁻⁸), suggesting that common environmental factors for inducing inflammatory response could underlie heart rate variation as previous studies suggest [34,54].

Finally, we examined the influence of SNP density on covariance estimation. Specifically, we focused on 13 phentoypes that have a larger number of SNPs than the remaining phenotypes and analyzed the 78 trait pairs among them using the dense set of SNPs (details in Methods). In the analysis, GECKO identified 37 significant genetic correlations based on Bonferroni correction (p-value < 0.00064) while LDSC identified 24, 22 of which were also identified by GECKO. Importantly, the genetic correlation estimates obtained with the dense set of SNPs is highly consistent with the estimates obtained earlier (Spearman correlation = 0.865). In addition, 28 significant estimates are shared in common between the significant genetic correlations detected using the dense set of SNPs and the using early non-dense SNPs (S4 Table). Besides genetic corrections, the environmental correlation estimates obtained with the dense set of SNPs is highly correlated with the estimates obtained with the earlier non-dense set of SNPs, with Spearman correlation being 0.98 between the two. In addition, 15 significant environmental correlation estimates are shared in common between the significant environmental correlations detected using the dense set of SNPs and the early non-dense set of SNPs (S4 Table). Overall, SNP density does not appear to influence the GECKO results.