back to top

BACKGROUND

APPROACHES

Direct linear combination of correlated test statistics (dLC)

In general, let T= (T₁,T₂,…,T_k)^T denote a vector of K correlated test statistics obtained individually from each univariate analysis for a  specific trait against a genetic marker. Under most circumstances in current GWAS studies, T usually follows an asymptotical multivariate normal distribution with mean   )^T  and known or estimated covariance Σ, where Σ is a k × k symmetric matrix . Assume the null hypothesis we want to test in multivariate analysis is H₀: = 0. In other words, the genetic variant is not associated with any phenotype. In contrast, the general alternative hypothesis H₁ is at least one >0, k=1,…,k. Extending from O’Brien’s theory[33], we propose a new approach for combining correlated traits [40],[40] called direct linear combination of test statistics (dLC). The new test statistics of dLC can be written as:

 

Under the null hypothesis,  follows a  distribution with K degrees of freedom and can be effectively used to test the joint significance of dependent univariate test statistics.

Empirical linear combination of correlated test statistics (eLC)

As illustrated in Xu et al. [33], dLC may not have optimal power against specific alternatives resulting from the heavy tail of the  distribution. Therefore, we further proposed a data-driven empirical approach to combine correlated test statistics as:

]           (2)

where c is some given non-negative constant. The weight in this new test statistics will be optimally determined by the specific data structure. For instance, when c=0, the test statistics simply reduces into sum of squares of T_k. When c is relatively large, equal weight is assigned to each T_k. Ideally, we would like to find an optimal value of c, so the  performs as a linear combination of T_i when under the H₀; but, under the alternative H_A, more weight is given to the larger true T_i. The bona fide p-value for  then can be estimated by applying permutation or perturbation techniques (see Supplementary).

   Estimation of the Variance-Covariance Matrix  Σ

Several methods can be used to estimate the covariance matrix Σ of univariate test statistics. For simplicity, we demonstrate two estimation approaches in a bivariate scenario. Based on  the method used in Yang et al. [37], we first utilize the sample covariance matrix of the test statistics of all SNPs from univariate GWAS analyses as an approximation. Σ, the covariance matrix, thus can be estimated as:

where Z₁ consists of unbiased univariate test statistics of all the SNPs for the first trait on genome-wide scale, so does Z₂. On the other hand, Σ can be estimated by using generalized least squares from individual-level data, as suggested by O’Brien [34]. A similar approach is also demonstrated by Liu et al. [29].

In this paper, the combining tests utilizing a covariance matrix approximated from the summary GWAS test statistics will be denoted by OB_z, dLC_z, and eLC_z. Similarly, those with a covariance matrix calculated from individual-level data will be referred to as OB_i, dLC_i, and eLC_i.

Simulation Study

Monte-Carlo Simulations were employed to generate data for evaluating the validity and performance of all the multivariate methods, especially the proposed dLC and eLC approaches. The main scenarios and key parameters for simulating genotype data are shown in Table 1. Various conditions were considered, including equal genetic marker contributes to both traits (Scenario A) and unequal genetic effects (Scenario B). To evaluate the proposed pleiotropy test strategy in Stage 2, a special Scenario C was also introduced. In particular, a single continuous trait was first generated with an assigned effect size. The second trait was then simulated by artificially adding random noise on the first generated trait. Therefore, the genetic variant is directly associated with the first trait but indirectly linked with the second trait, as illustrated in R2 and R3 in Figure 2. In all simulations, only quantitative traits and unrelated subjects were generated. A sample size of 1000 subjects was simulated in each of 1000 replicates.

The effect size of the SNP on each trait is estimated from   , where p is the minor allele frequency (MAF) at the locus and  , heritability, is the phenotypic variance explained by the SNP. The heritability is simulated at 1% and 2% in the Simulation Series I, and 0.5% and 0.1% in Simulation Series II for each trait, respectively. The variance of the environmental effects was fixed as 1 in all simulation studies. The genotypes were then generated under Hardy-Weinberg equilibrium with a specified MAF at 10%.

Bivariate quantitative phenotypes were randomly drawn from a bivariate normal distribution to represent the pleiotropic relationship R1 in Figure 1:

Y= ~ N

where  is the additive genetic effect size for trait i, X is the additive score of the coded allele, and  is the residual correlation between  and ,respectively.

Various simulation scenarios with respect to  were further generated. Briefly,  was selected at -0.75, -0.25, 0.01, 0.25, and 0.75 to mimic the correlations we have observed in real GWAS data analyses in the Simulation Series I. A smaller range of  at -0.5, -0.25, 0.01, 0.25 and 0.5 was also employed in Simulation Series II with lower genetic heritability. Notably, a mixture of protective and deleterious genetic effects was introduced when negative  was used in generating data. For instance, we assigned β₁=- β₂ as the alternative hypothesis when  was -0.25.

In addition to bivariate scenarios, we also generated three moderately correlated quantitative traits to evaluate the performance of various approaches when more than two traits were combined. For simplicity, the genetic heritability was equivalently set at 0.5% for each of the three traits. And the pairwise phenotypic residual correlation was chosen at 0.25, -0.15, and -0.20. The data were then generated under the null hypothesis H_o: β₁=β2=β₃=0 and two alternative hypotheses, H1: β₁=-β₂=-β₃ >0 and H₁: β₁=β₂=-β₃ >0, respectively.

Each replicate was analyzed by using OB, dLC, and eLC approaches individually. We estimated the covariance matrix Σ through approximation from summary data and individual-level data. Other multivariate methods were also compared, including LME with a random effect accounting for phenotypic correlations, GEE, MANOVA, and PCA, in which the first component was used as a dependent variable in the subsequent linear model. All analyses were conducted in R software (http://r-project.org/). Power and type I error rates of each approach were calculated as the proportion of replicates with a p-value less than a given significant threshold in the corresponding scenarios. Specifically, power was derived with 1000 replicates with the significance level at p-value equal to 10^-4 for each true scenario in Simulation Series I, and 10^-2 for those in Simulation Series II. Type I error rates were estimated in the settings that β₁=β₂=0 with 1000 replicates at nominal significance levels of 0.05 and 0.01. In the context of adjusting multiple testing for two univariate association tests of two phenotypes, standard Bonferroni corrections were applied with the significance level at 5*10^-5 and 5*10^-3 for Simulation Series I and II separately.

 

SIMULATIONS

Part I.  Performance of multivariate methods

Valid type I error rates, as expected, were consistently observed across all simulation scenarios at different nominal levels for all proposed and existing multivariate approaches using individual-level data (Table 2). In particular, our proposed approaches, direct linear combining (dLC) and empirical combining (eLC), demonstrate robustness under the null hypothesis, regardless of the directions of phenotypic residual correlation and approaches employed to estimate the covariance matrix Σ.

The estimated power using individual-level data in this series is presented in Table 3. Simulation Series I assessed the performance of various multivariate approaches under relevant hypotheses, in which relatively large genetic effects were simulated. The results are presented with respected to two alternative hypotheses, H₁:β₁=|β|₂>0 and H₁:β₁>|β₂|>0. In particular, we assigned a mixture of protective and deleterious genetic effects when negative phenotypic residual ρ was adopted in the simulation. All the methods performed well when test statistics were homogeneous, but their power varied considerably under the hypothesis of mixed genetic effects. For instance, our proposed approaches showed comparable power to MANOVA under the first alternative hypothesis; GEE, LME, and the O’Brien method (OB) had low power. In addition, principle component analysis (PCA) outperformed other methods when traits were highly correlated and showed less power than our proposed approaches when traits were barely related (ρ=0.01). The results under the second alternative hypothesis produced similar conclusions across various simulation scenarios. In general, our proposed approaches are substantially superior to univariate analysis approaches under all alternative hypotheses. Additionally, eLC outperformed dLC marginally in this simulation series. 

We next examined to what extent power is changed by applying our proposed approach to aggregated data under the same alternative hypotheses. MANOVA, PCA, GEE, and LME are unable to analyze aggregative results. Thus, we compared power between only OB, dLC, and eLC (Table 4). OB was inferior when the univariate test statistics were heterogeneous. dLC and eLC were not only superior to OB when directions of effects were opposite, but nearly as powerful as OB when genetic effects are homogenous. Moreover, similar power was consistently observed by applying two distinct methods to estimate covariance matrix Σ. Again, eLC and dLC showed significant advantages over univariate methods in detecting pleiotropic effects. Similar conclusions could be drawn on the results from Simulation Series II, where smaller genetic effects were simulated (see Supplementary Tables S1, S2). These results demonstrate that our proposed methods have comparable power to other approaches when all effects are in the same direction and much greater power when genetic effects are mixed. These can also be directly and effectively applied in aggregated data, thus efficiently increasing the sample size and power.

We extended our proposed approaches from combining two to combining three correlated traits. We generated data in the same manner as for the bivariate analysis under the null hypothesis, H_o: β₁=β2=β₃=0, and two alternative hypotheses separately, H₁: β₁=-β₂=-β₃ >0 and H₁: β₁=β₂=-β₃ >0. We excluded PCA in this series; because two principal components are typically needed to explain at least 80% of the total phenotypic variation, an increased number of subsequent regression tests and adjustments for multiple comparisons are necessary. For individual level-data, among classical methods MANOVA had better power compared to other conservative approaches, including GEE, LME, and OB (Table 5). Our proposed approaches (dLC and eLC) demonstrated comparable power to MANOVA under both alternative hypotheses. However, dLC had less power than eLC in this scenario, likely due to the increased degree of freedom of its underlying chi-square distribution. Additionally, our methods outperform OB using aggregated data when individual-level data are not available.

Part II.  Stage 2 Testing of Pleiotropy

We utilized the Simulation Series II dataset (low effect size described in Table 1) to act as positive and negative controls in evaluating the performance of our proposed stage 2 method of testing pleiotropy under relevant hypotheses.

A specific null hypothesis, “β₁>0,β₂=β₁+ε”, was simulated and used as the negative control. Table 6 presents the estimated type I error rates and power of our proposed approach (conditional testing of pleiotropy, cPLT) in comparison with the performance of the two conditional models strategy. Under the null hypothesis, cPLT demonstrated reasonable type I error rates, regardless of phenotypic residual correlations and effect directions.

Two alternative hypotheses, “β₁= |β₂|>0” and “β₁> |β₂|>0”, served as positive controls for detecting pleiotropy. cPLT demonstrated better performance than typical conditional strategies under these two alternative hypotheses (Table 6). However, the power was low for almost all methods when the degree of phenotypic correlation was high. Superior power was consistently found using cPLT under the first alternative hypothesis, in which the genetic variant contributed to phenotypes equally and the phenotypic correlation was modest. In general, statistical power was highest when phenotypic residual correlation was very small, e.g., ρ =0.01. Further, power of cPLT decreased as the degree of phenotypic residual correlation increased. Interestingly, power did not seem to be as strongly influenced by the mixture of protective and deleterious genetic effects in this stage as it was in the first stage.

These findings show that our proposed cPLT method performs well under the intended hypothesized pleiotropic relationship between genetic variants and correlated phenotypes, but its power can vary significantly as the correlation of phenotypes changes.

Joint stage 1 and 2 reducing type 1 error of identifying pleiotropic effects:

We further extensively examined the robustness of our proposed two-stage strategy under the null hypothesis, β₁>0,β₂=β₁+ε. Specifically, we presented type 1 error rates for separate and combined strategies with respect to various alpha levels (Table 7). Type 1 error rates were inflated for dLC applied alone in Stage 1, as expected. However, applying the proposed cPLT approach in Stage 2 would efficiently reduce the false positive results introduced in Stage 1. Notably, the overall false-positive rate was better controlled at the level of 0.05 when the screening threshold at Stage 1 was chosen at 0.025. Thus, this two-stage strategy would allow us to use less-stringent cut-off p-values in the screening stage to increase our power while effectively maintaining the overall false-positive rates through Stage 2.

Part III. An Application to GAW16 simulation data

To validate our novel two-stage approach, we applied the method to the first replicate of the GAW16 Problem 3 simulation dataset. We sought to identify pleiotropic variants (SNPs) associated with two blood lipids, high-density lipoprotein (HDL) and triglycerides (TG), known to influence cardiovascular disease and measured in the Framingham Heart Study, from which the dataset originates. The phenotypic correlation for HDL and TG was -0.28, and the estimated correlation of univariate test statistics of HDL and TG was -0.32. A unified genome-wide significance level was defined by false discovery rate (FDR) at p-value=3*10^-6, equivalent to q=0.05, for both bivariate and univariate analysis to screen potential pleiotropic SNPs. Similar performance of dLC and eLC and computational efficiency led us to use dLC to analyze combined correlated univariate GWAS results. The genomic control parameter was 1.01 for this bivariate GWAS analysis. Q-Q plot of the bivariate GWAS results (Supplementary Figure S1) showed no inflation beyond that expected by chance alone. A further comparison of bivariate results from dLC between individual-level data and aggregated data (Supplementary Figure S2) revealed no substantial differences.

The 25 leading SNPs in each genome-wide significant locus from the first stage were used in the second stage to infer their independent effects on HDL and TG by applying our proposed cPLT method with 10,000 permutations (Table 8). The significance level in this stage was 0.002 after Bonferroni correction. The two-stage analysis identified 2 SNPs, rs3200218 in the coding region of LPL and rs8192719 on the exon/intron boundary of CYP2B6. These loci served as positive controls in the GAW16 simulation dataset, and were validated as pleiotropic genes using our strategy. The simulated heritability of rs3200218 for HDL and TG was 0.3% and 0.4%, respectively. For rs8192719, the heritability for HDL and TG was reported as 0.3% for both [42]. Additionally, we detected a negative control SNP, rs7031748. Although it demonstrated a significant p-value from bivariate analysis in the first stage, it failed to reject the null hypothesis in our proposed method in the second stage. Therefore, the observed significant bivariate association of this genetic variant in the first stage likely resulted from indirect correlations, rather than a causal pleiotropic relationship.

Type 1 errors when filtering out SNPs with bivariate p-value > univariate p-value

Commonly, the multivariate p-value of a true pleiotropic marker is proposed to be smaller than that of both its univariate results. We investigated the plausibility of this hypothesis by comparing the distributions of p-values from both combining methods using our simulation data under the alternative hypotheses. The proportions of SNPs with smaller p-values in dLC under various scenarios are provided in Supplementary Table S3. In general, most p-values from dLC were smaller when the marker equally contributed to both traits and the phenotypic correlation was modest. On the other hand, p-values from dLC did not tend to be smaller when the phenotypic correlation was high or unequal genetic effect sizes existed. Adopting an additional significance threshold added little improvement. We further investigated the impact on statistical power by using this filter under relevant alternative hypotheses (Table 9). Applying the small multivariate p-value filtering criteria would introduce modest to significant power loss, largely depending on the magnitudes of phenotypic correlations and genetic effect sizes.