MGAS: a powerful tool for multivariate gene-based genome-wide association analysis

Motivation: Standard genome-wide association studies, testing the association between one phenotype and a large number of single nucleotide polymorphisms (SNPs), are limited in two ways: (i) traits are often multivariate, and analysis of composite scores entails loss in statistical power and (ii) gene-based analyses may be preferred, e.g. to decrease the multiple testing problem. Results: Here we present a new method, multivariate gene-based association test by extended Simes procedure (MGAS), that allows gene-based testing of multivariate phenotypes in unrelated individuals. Through extensive simulation, we show that under most trait-generating genotype–phenotype models MGAS has superior statistical power to detect associated genes compared with gene-based analyses of univariate phenotypic composite scores (i.e. GATES, multiple regression), and multivariate analysis of variance (MANOVA). Re-analysis of metabolic data revealed 32 False Discovery Rate controlled genome-wide significant genes, and 12 regions harboring multiple genes; of these 44 regions, 30 were not reported in the original analysis. Conclusion: MGAS allows researchers to conduct their multivariate gene-based analyses efficiently, and without the loss of power that is often associated with an incorrectly specified genotype–phenotype models. Availability and implementation: MGAS is freely available in KGG v3.0 (http://statgenpro.psychiatry.hku.hk/limx/kgg/download.php). Access to the metabolic dataset can be requested at dbGaP (https://dbgap.ncbi.nlm.nih.gov/). The R-simulation code is available from http://ctglab.nl/people/sophie_van_der_sluis. Contact: mxli@hku.hk Supplementary information: Supplementary data are available at Bioinformatics online.


Content
Supplemental Information Simulations p3 Figure S1: The relationship between r*ρ (i.e., correlations between genotypes r and correlations between phenotypes ρ: x-axis) and correlations between p-values (y-axis) p8 Figure S2: Example QQplots of 4 whole-genome multivariate gene-based Type I error rate simulations p9 in which real genome-wide genetic data were used and simulated phenotypes with a correlational structure similar to that of the NFBC1966 metabolic phenotypes (Table S8). Figure S3: QQplot of the orginal p-values (provided by the authors (Sabatti et al) and available in the dbGaP data base) per phenotype. p10 Figure S4: QQplot of the univariate SNP-based p-values (combined across all phenotypes), subdivided in all SNPs (green), SNPs in genes (orange), and SNPs outside genes (pink). The multivariate gene-based p-values calculated by MGAS are printed in blue. p11 Table S1a: Visual representation of the LD patterns used for the simulations including small genes (10 SNPs) p12 Table S1b: Visual representation of the LD patterns used for the simulations including small genes (10 SNPs) p13 Table S1c: Visual representation of the LD patterns used for the simulations including large genes (60 SNPs) p14 Table S1d: Visual representation of the LD patterns used for the simulations including large genes (60 SNPs) p15 Table S1e: Visual representation of the LD patterns used for the simulations including large genes (60 SNPs) p16 Table S2: 1 factor model, gene-effect on factor p17 Table S3: 1 factor model, gene-effect specific to one phenotype p18 Table S4: 4 factor model, gene-effect on one of the four factors p19 Table S5: 4 factor model, gene-effect specific to one phenotype p20 Table S6: Network model mimicking 1-factor model (i.e., all phenotypes equally strongly related), gene-effect specific to one phenotype p21 Table S7: Network model mimicking 4-factor model (i.e., clusters of phenotypes), gene-effect specific to one phenotype p22 Table S8: Phenotypic correlations between 9 metabolic traits, corrected for sex, oral contraceptive use and pregnancy p23 Table S9: List of FDR corrected genes/genetic regions identified by MGAS in the metabolic data Table S10: Annotation of the genes/genetic regions identified by MGAS by reference to the Catalog of Published Genome-wide Association Studies (http://www.genome.gov/gwastudies/) Table S11: Genome-wide significant SNPs identified by TATES using the 9 continuous metabolic phenotypes Table S12: Genome-wide significant genes identified by GATES on a) the univariate sum score of five metabolic syndrome risk factors and b) the individual five dichotomous metabolic risk factors Table S13: Genome-wide significant genes identified by GATES on the 9 continuous metabolic phenotypes separately Table S14: Comparison results MGAS, TATES, and GATES when all analyses are conducted on the original 9 continuous metabolic phenotypes For the large gene, 4 scenarios were simulated. In each of these scenarios, the gene covered 8 LD blocks: 4 blocks of 10 SNPs each, and 4 blocks of 5 SNPs each.
1. 0 DSL (to control Type I error rate); 2. 1 DSL in a large block; 3. 8 DSL, 1 in each LD block; 4. 8 DSL, 1 in each block, opposite effects**. See Supplemental Tables S1a-e for a visual representation of these LD settings.
Note that, assuming the same number of DSL with the same effect size, the power to detect a gene is negatively correlated to the size of the gene. Specifically, because the weight q e , denoting the effective number of p-values (i.e., tests and thus the effective number of phenotypes + SNPs), will, given the same LD structure, invariably be higher for larger genes, the weight q e /q ej will also be higher. As a consequence, p-values of DSL with the same effect size will be multiplied by a larger weight when these DSL are located in a larger gene. That is, a stronger genetic signal is required to implicate a large gene. To obtain illustrative power results, we therefore chose to simulate the data such that the DSL in small genes together explained .5% of the variance, and the DSL in large genes together explained 1% variance. As a consequence, power results as depicted in Figure 2 and reported in Tables S2-S7 are not directly comparable across genes of different size.
When a simulated gene included multiple DSL, the effect was evenly divided over the multiple DSL (e.g., in a large gene with 2 DSL, both DSL were simulated to explain 1/2=.5% variance).
Important to note is that because of the strong LD between the SNPs within an LD block, SNPs that were simulated to not have an effect are related to the phenotype if one regresses the phenotype on that SNP because of that SNP's strong LD with the actual DSL(s) in the block (i.e., the regression coefficient β in that regression will be equal to the β of the actual DSL weighted by the LD measure r). In addition, when multiple DSL are simulated within a block, the effect of every DSL is augmented in the simple regression of a phenotype on that DSL because of the DSL's correlation to other DSL in a block (i.e., the β in that regression is equal to the simulated β of that DSL plus the β of other DSL in the block weighted by the LD measure r). As a consequence, a) the variance explained by the entire gene (the R 2 as calculated in e.g. multiple regression) was, given our simulation settings, larger for genes harbouring multiple DSL, and b) in the simple regression of a phenotype on a single DSL, the variance explained by the DSL was almost equal across scenarios 2, 3 and 4 of the small gene (i.e., when 1, 2 or 4 DSL resided in the same LD-block).

Phenotypic data Phenotypic factor models
With m phenotypes and k common factors, phenotypic data were simulated according to the following model: where Σ denotes the mxm phenotypic variance-covariance matrix, Λ denotes the mxk matrix of factor loadings (with t denoting matrix transpose), Ψ denotes the kxk covariance matrix between the k latent factors, and Θ denotes the diagonal mxm matrix of residual variances (i.e., the part of the phenotypic variances that is not explained by any of the k factors). In models with k>1, simple structure was maintained (i.e., all phenotypes loaded on only 1 latent factor).
Data were generated for models featuring either 1 or 4 factors. In the 1-factor models, factor loadings were set to .75, such that the latent factor explained .75 2 =.56 % of the variance in each phenotype, and residual variances were equal to 1-.75 2 =.44. Consequently, all intercorrelations between the 20 standard normally distributed phenotypes were .56. In the 1-factor model, all SNPs within a gene were either modelled to affect the latent factor (and through the latent factor all 20 phenotypes underlying the factor), or to affect one of the 20 phenotypes directly (see Figures 1a and 1b, respectively).
Note that when 1) all phenotypic intercorrelations are explained by 1 latent factor, 2) all phenotypes have equal factor loadings, and 3) all phenotypes have equal residual variances, the model is called a Rasch model (Rasch, 1980). Important to note is that only when the phenotypic data-generating model is a Rasch model and the effect of gene X is on the latent factor (i.e., Figure  1a), is the sum score calculated across all individual phenotypes a sufficient statistic for the genetic analysis of gene X, i.e., only then does the sum of scores on all phenotypes exhaustively summarize all phenotypic information that is relevant for the analysis testing the association with gene X. The sum score is, however, not a sufficient statistic for the genetic analysis when the genetic effect is directly on one of the phenotypes, or when the phenotypic Rasch model does not hold.
In the 4-factor models, 5 phenotypes were modelled per factors. All factor loadings were set to .75 within factors, all residual variances equalled 1-.75 2 =.44, and correlations between factors were set to .23. Consequently, correlations between phenotypes loading on the same factor equalled .75 2 =.56, while correlations between phenotypes loading on different factors equalled .75*.23*.75=.13. Effects of all SNPs within a gene were either modelled on one of the four latent factors (and through that factor on all 5 phenotypes defining that factor), or directly on one of the phenotypes (see Figures 1c and 1d, respectively). Note that the phenotypic sum score is not a sufficient statistic for the genetic analysis when the phenotypic model is a 4-factor model.

Phenotypic network models
In all network simulations, we assumed the model to be stationary, i.e., the mutual interactions between the phenotypes have over time resulted in a stable variance-covariance matrix. Also, we assumed that all phenotypes were mutually related, i.e., in all regressions of the phenotypes on each other β ij ≠0. We did, however, not allow instantaneous self-activation, i.e., the phenotypes were not allowed to affect themselves on the next time point (i.e., β ii =0). (Note, however, that indirect selfactivation on subsequent time points was allowed, i.e., phenotype A affects phenotype B on time point t=1, and as all phenotypes were simulated to be mutually related, it affects itself via phenotype B on t=2). Network data were simulated according to the model: where Σ denotes the mxm phenotypic variance-covariance matrix, I is an mxm identity matrix, B is a full mxm matrix including the regression weights β of all phenotypes on each other (i.e., element B[i,j] includes the regression weight of phenotype i on phenotype j). The diagonal of B was fixed to zero to model absence of instantaneous self-activation. Ψ is the mxm diagonal matrix containing the variances of all phenotypes conditional on the other phenotypes.
Note that in network models, in contrast to the factor models, phenotype-specific genetic effects propagate through the network to indirectly affect other connected phenotypes in the network. In phenotypic network models, genes thus affect all connected phenotypes, with the size of the effect depending on the strength of the phenotypic relations (i.e., edges) in the network. In all network simulations, all SNPs within a gene affected only the first phenotype in the network, but because all phenotypes are connected, the genetic effects do spread through the network. Two types of networks were simulated. First, all regression weights in the B matrix were set to .04197, resulting in phenotypic intercorrelations among all 20 phenotypes of .56. That is, the phenotypic variance-covariance matrix of these network data is similar to the phenotypic variance-covariance matrices resulting from the 1-factor Rasch models discussed above. Second, a network was simulated with 4 clusters of phenotypes that correlated .56 within clusters, and .13 between clusters, resulting in a phenotypic variance-covariance matrix that is similar to that resulting from the 4-factor models discussed above. The network models are illustrated in Figures 1e and 1f, respectively).
It is important to note that network models and factor models can yield phenotypic data with very similar variance-covariance structure, in spite of the different underlying data-generating processes. As a consequence, even if the variance-covariance structure of data can be well described by a 1-factor model, this does not necessarily imply that the 1-factor model is the actual datagenerating model. This is relevant for genetic studies, because factor analytic results are often used to justify the reduction of multivariate data to univariate composite scores (e.g., sum scores, factor scores). Such reduction is, however, only valid if the data-generating process is a 1-factor model, but not if the data are generated by a network process.
Rasch G (1980) Probabilistic models for some intelligence and attainment tests. Chicago: The University of Chicago Press.

Figure S1
The relationship between r*ρ (i.e., correlations between genotypes r and correlations between phenotypes ρ: x-axis) and correlations between p-values (y-axis) obtained in the regression of a phenotype on a genetic variant. Simulations (see original manuscript §2.2 for details) show that this relationship can be accurately described by a 6 th order polynomial (coefficient of determination R 2 =.995). Note that the data points produced for the quantitative and qualitative traits are mixed together to generate a unified approximate function. Note also that using the same 6th order polynomial, similarly high R 2 values were obtained for different simulation settings. For instance, when the simulated sample size was only N=2000, the R 2 was 0.998, and when 4 phenotypes and 4 SNPs were simulated in a sample of 4000 subjects, the R 2 was 0.997.

Figure S2
Example QQplots of 4 of the 10 whole-genome multivariate gene-based simulations in which we used real genome-wide genetic data of N=4,763 subjects from the Northern Finland Birth Cohort (NFBC1966, Sabatti et al., 2009, and simulated 10 sets of 9 phenotypes, of which the correlational structure mimicked the observed correlations of the NFBC1966 metabolic phenotypes (see Supplemental Table S8), to check whether MGAS Type I error rates were correct regardless of LD pattern and gene size. Plotted are the univariate SNP-based p-values (i.e., pvalues of all 9 phenotypes in one plot), subdivided in all SNPs (green), SNPs in genes (orange), and SNPs outside genes (pink). The multivariate gene-based p-values (blue) are calculated through the MGAS procedure and are assembled from the SNPs located within genes only.

Figure S3
QQplot of the orginal p-values per phenotype obtained in covariate and PC corrected analyses, available in the dbGaP data base as provided by the original authors: Sabatti, C., Service, S.K., Hartikainen, A-L, Pouta, A., Ripatti, S., Brodsky, J., et al. Genome-wide association analysis of metabolic traits in a birth cohort from a founder population. Nat. Genet., 41, 35-46 (2009).

Figure S4
QQplot of the univariate SNP-based p-values (i.e., p-values of all 9 phenotypes in one plot), subdivided in all SNPs (green), SNPs in genes (orange), and SNPs outside genes (pink). The multivariate gene-based p-values (blue) are calculated through the MGAS procedure and are assembled from the SNPs located within genes only: the enrichment of association signals in the multivariate gene-based p-values is due to the selection of smallest weighted p-values by MGAS. Table S1a Visual representation of the LD patterns used for the simulations including small genes (10 SNPs) Gene 1 Gene 2 Gene 3 Scenario 1 0 DSL (to control Type I error rate) 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 .9 1 Scenario 2/3 1 DSL(in red) 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 .9 1 Scenario 4/5 2 DSL(in red) 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 .9 1 Note. Within blocks (grey shading) r=.9, between block (no shading), r=0. Disease Susceptibility Loci (DSL) are indicated in red. All small genes including DSL explain .5% of the phenotypic variance; when a gene includes multiple DSL, the .5% is evenly distributed across all DSL (e.g., in scenario 3, the 2 DSL each explain .25% of the variance). The scenario numbers refer to those in Table 1 of the original manuscript. Table S1b Visual representation of the LD patterns used for the simulations including small genes (10 SNPs) Gene 4 Gene 5 Gene 6 Scenario 6/7 4 DSL (in red) 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 1 .9 .9 .9 .9 .9 .9 .9 .9 .9 1 Scenario 8 1 DSL (in red) 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 Scenario 9/10 2 DSL(in red) 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 1 .9 1 .9 .9 1 .9 .9 .9 1 .9 .9 .9 .9 1 Note. Within blocks (grey shading) r=.9, between block (no shading), r=0. Disease Susceptibility Loci (DSL) are indicated in red. All small genes including DSL explain .5% of the phenotypic variance; when a gene includes multiple DSL, the .5% is evenly distributed across all DSL (e.g., in scenario 4, the 4 DSL each explain .125% of the variance). The scenario numbers refer to those in Table 1 of the original manuscript. Table S1c Visual representation of the LD patterns used for the simulations including larger genes (60 SNPs) Scenario 11: 0 Disease Susceptibility Loci (to control Type I error rate) Note. A gene covering 60 SNPs divided over 4 large (10 SNPs) and 4 small (5 SNPs) LD blocks. Within blocks (grey shading) r=.9, between block (no shading), r=0. The scenario numbers refer to those in Table 1 of the original manuscript. Table S1d Visual representation of the LD patterns used for the simulations including large genes (60 SNPs) Scenario 12: 1 Disease Susceptibility Locus (in red) explaining 1% of the variance Note. A gene covering 60 SNPs in 4 large (10 SNPs) and 4 small (5 SNPs) LD blocks. Within blocks (grey shading) r=.9, between block (no shading), r=0. Disease Susceptibility Loci (DSL) are indicated in red. The scenario numbers refer to those in Table 1 of the original manuscript. Table S1e Visual representation of the LD patterns used for the simulations including large genes (60 SNPs) Scenario 13/14: 8 Disease Susceptibility Loci (in red) together explaining 1% of the variance (i.e., .125% each) Note. A gene covering 60 SNPs in 4 large (10 SNPs) and 4 small (5 SNPs) LD blocks. Within blocks (grey shading) r=.9, between block (no shading), r=0. Disease Susceptibility Loci (DSL) are indicated in red. The scenario numbers refer to those in Table 1 of the original manuscript. Note. DSL=Disease Susceptibility Loci. All Nrep=2000 simulations included N=2000 subjects, Nphen=20 phenotypes, and small or large genes covering N SNP =10 or N SNP =60 SNPs, respectively. Simulated data were analyzed in 5 ways all yielding a gene-based p-value. Sum Gates: original Gates procedure conducted on the univariate sum score calculated across the 20 phenotypes; Sum multiple regression: univariate phenotypic sum score regressed on all SNPs within a gene; MANOVA: all 20 phenotypes as dependent variables and all SNPs within a gene as predictors; GATES-MANOVA: all 20 phenotypes as dependent variables and each SNP separately as predictor in MANOVA, the NSNP MANOVA p-values are then combined using GATES; MGAS: weighted selection of the mxn p-values obtained in regression of all m phenotypes on all n SNPs within a gene, while correcting for the expected correlations between these p-values resulting from the phenotypic correlations and the LD-structure within a gene. MAF was .2 for all simulated SNPs. Within and between LD-blocks, LD measure r was set to .9 and to 0, respectively. As the MGAS procedure has less power to detect larger genes (see Supplemental Information), small genes and large genes explained .5% and 1% of the variance, respectively. As a consequence, power results are not directly comparable between small and large genes. When a gene harboured multiple DSL, the total gene effect was equally distributed across all DSL. * The Roman numerals refer to Figure 2 in the paper. ** Half of the DSL in a gene conveyed a positive effect, the others a negative effect on the factors (Figure 1a,c) or phenotypes (Figure 1c,d,e,f). *** In latent factor models ( Figure 1a,c), the one DSL affected half of the phenotypic indicators of a factor positively and the other half negatively. In unclustered networks (Figure 1e), the one DSL affected 1 phenotype positively and another phenotype in the network negatively. In clustered networks (Figure 1f), the one DSL affected either 2 phenotypes in the same cluster, or 2 phenotypes in different clusters, in opposite directions. Note. DSL=Disease Susceptibility Loci. All Nrep=2000 simulations included N=2000 subjects, Nphen=20 phenotypes, and small or large genes covering N SNP =10 or N SNP =60 SNPs, respectively. Simulated data were analyzed in 5 ways all yielding a gene-based p-value. Sum Gates: original Gates procedure conducted on the univariate sum score calculated across the 20 phenotypes; Sum multiple regression: univariate phenotypic sum score regressed on all SNPs within a gene; MANOVA: all 20 phenotypes as dependent variables and all SNPs within a gene as predictors; GATES-MANOVA: all 20 phenotypes as dependent variables and each SNP separately as predictor in MANOVA, the NSNP MANOVA p-values are then combined using GATES; MGAS: weighted selection of the mxn p-values obtained in regression of all m phenotypes on all n SNPs within a gene, while correcting for the expected correlations between these p-values resulting from the phenotypic correlations and the LD-structure within a gene. MAF was .2 for all simulated SNPs. Within and between LD-blocks, LD measure r was set to .9 and to 0, respectively. As the MGAS procedure has less power to detect larger genes (see Supplemental Information), small genes and large genes explained .5% and 1% of the variance, respectively. As a consequence, power results are not directly comparable between small and large genes. When a gene harboured multiple DSL, the total gene effect was equally distributed across all DSL. * The Roman numerals refer to Figure 2 in the paper. ** Half of the DSL in a gene conveyed a positive effect, the others a negative effect on the factors (Figure 1a,c) or phenotypes (Figure 1c,d,e,f). *** In latent factor models ( Figure 1a,c), the one DSL affected half of the phenotypic indicators of a factor positively and the other half negatively. In unclustered networks (Figure 1e), the one DSL affected 1 phenotype positively and another phenotype in the network negatively. In clustered networks (Figure 1f), the one DSL affected either 2 phenotypes in the same cluster, or 2 phenotypes in different clusters, in opposite directions. Note. DSL=Disease Susceptibility Loci. All Nrep=2000 simulations included N=2000 subjects, Nphen=20 phenotypes, and small or large genes covering N SNP =10 or N SNP =60 SNPs, respectively. Simulated data were analyzed in 5 ways all yielding a gene-based p-value. Sum Gates: original Gates procedure conducted on the univariate sum score calculated across the 20 phenotypes; Sum multiple regression: univariate phenotypic sum score regressed on all SNPs within a gene; MANOVA: all 20 phenotypes as dependent variables and all SNPs within a gene as predictors; GATES-MANOVA: all 20 phenotypes as dependent variables and each SNP separately as predictor in MANOVA, the NSNP MANOVA p-values are then combined using GATES; MGAS: weighted selection of the mxn p-values obtained in regression of all m phenotypes on all n SNPs within a gene, while correcting for the expected correlations between these p-values resulting from the phenotypic correlations and the LD-structure within a gene. MAF was .2 for all simulated SNPs. Within and between LD-blocks, LD measure r was set to .9 and to 0, respectively. As the MGAS procedure has less power to detect larger genes (see Supplemental Information), small genes and large genes explained .5% and 1% of the variance, respectively. As a consequence, power results are not directly comparable between small and large genes. When a gene harboured multiple DSL, the total gene effect was equally distributed across all DSL. * The Roman numerals refer to Figure 2 in the paper. ** Half of the DSL in a gene conveyed a positive effect, the others a negative effect on the factors (Figure 1a,c) or phenotypes (Figure 1c,d,e,f). *** In latent factor models ( Figure 1a,c), the one DSL affected half of the phenotypic indicators of a factor positively and the other half negatively. In unclustered networks (Figure 1e), the one DSL affected 1 phenotype positively and another phenotype in the network negatively. In clustered networks (Figure 1f), the one DSL affected either 2 phenotypes in the same cluster, or 2 phenotypes in different clusters, in opposite directions. Note. DSL=Disease Susceptibility Loci. All Nrep=2000 simulations included N=2000 subjects, Nphen=20 phenotypes, and small or large genes covering N SNP =10 or N SNP =60 SNPs, respectively. Simulated data were analyzed in 5 ways all yielding a gene-based p-value. Sum Gates: original Gates procedure conducted on the univariate sum score calculated across the 20 phenotypes; Sum multiple regression: univariate phenotypic sum score regressed on all SNPs within a gene; MANOVA: all 20 phenotypes as dependent variables and all SNPs within a gene as predictors; GATES-MANOVA: all 20 phenotypes as dependent variables and each SNP separately as predictor in MANOVA, the N SNP MANOVA p-values are then combined using GATES; MGAS: weighted selection of the mxn p-values obtained in regression of all m phenotypes on all n SNPs within a gene, while correcting for the expected correlations between these p-values resulting from the phenotypic correlations and the LD-structure within a gene. MAF was .2 for all simulated SNPs. Within and between LD-blocks, LD measure r was set to .9 and to 0, respectively. As the MGAS procedure has less power to detect larger genes (see Supplemental Information), small genes and large genes explained .5% and 1% of the variance, respectively. As a consequence, power results are not directly comparable between small and large genes. When a gene harboured multiple DSL, the total gene effect was equally distributed across all DSL. * The Roman numerals refer to Figure   Note. DSL=Disease Susceptibility Loci. All Nrep=2000 simulations included N=2000 subjects, Nphen=20 phenotypes, and small or large genes covering N SNP =10 or N SNP =60 SNPs, respectively. Simulated data were analyzed in 5 ways all yielding a gene-based p-value. Sum Gates: original Gates procedure conducted on the univariate sum score calculated across the 20 phenotypes; Sum multiple regression: univariate phenotypic sum score regressed on all SNPs within a gene; MANOVA: all 20 phenotypes as dependent variables and all SNPs within a gene as predictors; GATES-MANOVA: all 20 phenotypes as dependent variables and each SNP separately as predictor in MANOVA, the NSNP MANOVA p-values are then combined using GATES; MGAS: weighted selection of the mxn p-values obtained in regression of all m phenotypes on all n SNPs within a gene, while correcting for the expected correlations between these p-values resulting from the phenotypic correlations and the LD-structure within a gene. MAF was .2 for all simulated SNPs. Within and between LD-blocks, LD measure r was set to .9 and to 0, respectively. As the MGAS procedure has less power to detect larger genes (see Supplemental Information), small genes and large genes explained .5% and 1% of the variance, respectively. As a consequence, power results are not directly comparable between small and large genes. When a gene harboured multiple DSL, the total gene effect was equally distributed across all DSL. * The Roman numerals refer to Figure 2 in the paper. ** Half of the DSL in a gene conveyed a positive effect, the others a negative effect on the factors (Figure 1a,c) or phenotypes (Figure 1c,d,e,f). *** In latent factor models ( Figure 1a,c), the one DSL affected half of the phenotypic indicators of a factor positively and the other half negatively. In unclustered networks (Figure 1e), the one DSL affected 1 phenotype positively and another phenotype in the network negatively. In clustered networks (Figure 1f), the one DSL affected either 2 phenotypes in the same cluster, or 2 phenotypes in different clusters, in opposite directions. Note. DSL=Disease Susceptibility Loci. All Nrep=2000 simulations included N=2000 subjects, Nphen=20 phenotypes, and small or large genes covering NSNP=10 or NSNP=60 SNPs, respectively. Simulated data were analyzed in 5 ways all yielding a gene-based p-value. Sum Gates: original Gates procedure conducted on the univariate sum score calculated across the 20 phenotypes; Sum multiple regression: univariate phenotypic sum score regressed on all SNPs within a gene; MANOVA: all 20 phenotypes as dependent variables and all SNPs within a gene as predictors; GATES-MANOVA: all 20 phenotypes as dependent variables and each SNP separately as predictor in MANOVA, the NSNP MANOVA p-values are then combined using GATES; MGAS: weighted selection of the mxn p-values obtained in regression of all m phenotypes on all n SNPs within a gene, while correcting for the expected correlations between these p-values resulting from the phenotypic correlations and the LD-structure within a gene. MAF was .2 for all simulated SNPs. Within and between LD-blocks, LD measure r was set to .9 and to 0, respectively. As the MGAS procedure has less power to detect larger genes (see Supplemental Information), small genes and large genes explained .5% and 1% of the variance, respectively. As a consequence, power results are not directly comparable between small and large genes. When a gene harboured multiple DSL, the total gene effect was equally distributed across all DSL. * The Roman numerals refer to Figure 2 in the paper. ** Half of the DSL in a gene conveyed a positive effect, the others a negative effect on the factors (Figure 1a,c) or phenotypes (Figure 1c,d,e,f). *** In latent factor models ( Figure 1a,c), the one DSL affected half of the phenotypic indicators of a factor positively and the other half negatively. In unclustered networks (Figure 1e), the one DSL affected 1 phenotype positively and another phenotype in the network negatively. In clustered networks (Figure 1f), the one DSL affected either 2 phenotypes in the same cluster (X), or 2 phenotypes in different clusters (XI), in opposite directions. **** Across multiple simulation runs, this scenario showed a slight tendency towards inflation (~.063). .177 1 Note: Abbreviations are: body mass index (BMI), C-reactive protein (CRP), diastolic blood pressure (DBP), glucose (GLU), high-density Lipoprotein (HDL), insulin (INS), low-density lipoprotein (LDL), systolic blood pressure (SBP), and triglycerides (TG). These 9 traits are part of the Northern Finland Birth Cohort data (NFBC1966, Sabatti et al. 2009), found at http://www.ncbi.nlm.nih.gov/gap, controlled through dbGaP study accession number phs000276.v1.p1. The NFBC1966 Study is conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with the Broad Institute, UCLA, University of Oulu, and the National Institute for Health and Welfare in Finland. This manuscript was not prepared in collaboration with investigators of the NFBC1966 Study and does not necessarily reflect the opinions or views of the NFBC1966 Study Investigators, Broad Institute, UCLA, University of Oulu, National Institute for Health and Welfare in Finland, and the NHLBI.