|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: ASF GJM UD MPR. Performed the experiments: ASF GJM UD JFF. Analyzed the data: ASF GJM UD JFF. Contributed reagents/materials/analysis tools: ASF GJM UD RL. Wrote the paper: ASF MPR.
Informing missing heritability for complex disease will likely require leveraging information across multiple SNPs within a gene region simultaneously to characterize gene and locus-level contributions to disease phenotypes. To this aim, we introduce a novel strategy, termed Mixed modeling of Meta-Analysis P-values (MixMAP), that draws on a principled statistical modeling framework and the vast array of summary data now available from genetic association studies, to test formally for locus level association. The primary inputs to this approach are: (a) single SNP level p-values for tests of association; and (b) the mapping of SNPs to genomic regions. The output of MixMAP is comprised of locus level estimates and tests of association. In application of MixMAP to summary data from the Global Lipids Gene Consortium, we suggest twelve new loci (PKN, FN1, UGT1A1, PPARG, DMDGH, PPARD, CDK6, VPS13B, GAD2, GAB2, APOH and NPC1) for low-density lipoprotein cholesterol (LDL-C), a causal risk factor for cardiovascular disease and we also demonstrate the potential utility of MixMAP in small data settings. Overall, MixMAP offers novel and complementary information as compared to SNP-based analysis approaches and is straightforward to implement with existing open-source statistical software tools.
Serum lipid levels are established determinants of cardiovascular disease morbidity with well-described heritability. Indeed, meta-analysis of data arising from recent genome-wide association studies has identified common genetic variants in at least 95 loci associated with low-density lipoprotein cholesterol (LDL-C), high-density lipoprotein cholesterol (HDL-C), triglycerides (TG) and total cholesterol (TC) . In practice, selection of statistically relevant genes or loci, subsequently referred to collectively as loci, is often based on simple linear regression of single-nucleotide polymorphisms (SNPs); that is, loci within which there is at least one SNP that reaches genome-wide significance, defined according to a Bonferroni level correction for multiple testing, are regarded as significantly associated with the trait under study. While this approach is valid, we conjecture that substantial, complementary knowledge about association can be acquired by considering available information on all SNPs within a locus simultaneously in characterizing association.
In order to address this, we apply a mixed effects modeling paradigm that uses SNP-level meta-analysis p-values to arrive at formal analytic characterization of underlying locus-level effects on complex disease phenotypes. Through application of principled, well-vetted statistical concepts for valid and reliable inference, this approach, termed Mixed modeling of Meta-Analysis P-values (MixMAP), draws strength from available information across all SNPs within a pre-defined region without requiring first-stage data reduction. That is, it is not necessary to eliminate redundancy among SNPs in high linkage disequilibrium (LD) within a locus. By incorporating all available information about SNPs within a locus, MixMAP results in increased sensitivity for identifying loci with multiple SNPs of moderate significance, as evidenced in the applications described below.
Applications of hierarchical linear and generalized linear models for analysis of data arising from genetic association studies have been described previously. Specifically, use of mixed effects models for family studies is common to account for correlation arising from familial level clustering where a random effect for family is included in the model (see for example , ). In the context of population-based investigations of unrelated individuals, applications of mixed models are described with random effects for SNPs, genes and sets of genes (see for example –). A primary conceptual distinction between MixMAP and alternative approaches is that MixMAP uses summary level -values as the dependent variable in the hierarchical model rather than raw data, allowing investigators to leverage existing, publicly-accessible data resources. MixMAP also provides locus level tests and these are not limited by an increasing number of SNPs within a gene as described in . Technically most similar to MixMAP, the approach of Wang et al.  involves fitting a hierarchical model to summary level data, specifically -statistics from single-SNP tests of association; however, in addition to differences in model specification, hypothesis testing within the Wang et al. framework is based on fixed effects parameters (representing pathway effects) while MixMAP involves generating prediction intervals for latent (locus-level) effects.
The purpose of this manuscript is to describe the MixMAP algorithm and to illustrate its application under a range of conditions. Application of MixMAP to two independent datasets of SNP associations with LDL-C, a causal risk factor for atherosclerotic cardiovascular disease (CVD), reveals that MixMAP offers novel and complementary information as compared to traditional analysis approaches. Notably, as an analytic tool that selects loci with multiple moderate to strong signals, and one that uses single SNP-level summary data (-values), MixMAP is intended to serve as an additional analysis framework that complements, rather than replaces, single-SNP based investigations. This is further supported by simulation studies suggesting MixMAP has a higher true positive rate than single-SNP based analysis in the context of more moderate gene level signals. MixMAP is straightforward to implement with the open-source MixMAP package in R (freely-available for download at http://cran.r-project.org/ and http://people.umass.edu/foulkes/software.html.).
This section reports the results of applying MixMAP to data arising from two independent sources: the Global Lipids Gene Consortium (GLGC) and the Penn Coronary Artery Calcification (PennCAC). These revealed that application of MixMAP: (1) supports published LDL-C loci; (2) suggests novel LD-CL genes and (3) complements a single-SNP testing approach, as described in detail below. Additionally, the results of simulation studies designed to characterize the performance of MixMAP under a range of conditions, also described below, support the concept that MixMAP serves as a complementary analysis strategy. To begin, we briefly outline the MixMAP algorithm and the available GLGC and PennCAC data. Further details are provided in the Material and Methods.
MixMAP is a statistical framework that uses the results of single SNP-level analysis to test formally for locus-level association. The MixMAP approach is designed to identify loci involving multiple SNPs each with moderate effects on the trait that may not be detected by single SNP analysis. Loci that are detected by a single SNP approach may also be detected by MixMAP if multiple SNPs have modest association across the locus. The primary inputs to the MixMAP algorithm are: (a) p-values corresponding to single SNP tests of association with a trait (e.g. LDL-C); (b) a mapping of SNPs to genomic regions. Additional SNP and gene-level covariate information, such as gene size, number of SNPs per Kb, number of recombinant hotspots per Kb and average linkage disequilibrium can be incorporated.
The MixMAP algorithm is summarized as follows, with additional detail provided below: [Step 1:] Fit a mixed effects model to inverse normally transformed (ranked) -values, with random locus-specific intercepts; [Step 2:] Predict random locus-specific effects using empirical Bayes estimation; [Step 3:] Calculate corresponding prediction intervals using a Bonferroni corrected threshold; and [Step 4:] Report a locus as statistically meaningful if the upper limit of the prediction interval corresponding to the locus effect is less than . In the present manuscript, we fit the mixed model with random gene-level effects and provide a post-hoc characterization of loci (genes or groups of genes) based on the genes that emerge through application of MixMAP. As more comprehensive information on the relationship between genes within loci becomes available, the MixMAP approach is flexible in that an alternative loci-level annotation can be used as input to the algorithm.
The current applied investigation focuses on loci for LDL-C, an important causal factor for CVD. Data arising from two independent data-sets are considered: (1) reported and publicly available meta-analysis SNP level -values for association with LDL-C derived from multiple independent association of approximately individuals in the Global Lipids Gene Consortium (GLGC) study (http://www.broadinstitute.org/mpg/pubs/lipids2010/); and (2) SNP level -values for association with LDL-C derived from analysis of ITMAT-Broad-CARe (IBC) 50K SNP array data in European ancestry individuals within the Penn Coronary Artery Calcification (PennCAC) sample –, a relatively small study ( Caucasians) that is considered underpowered by itself to identify the genetic determinants of a complex disease phenotype. Additional details on this cohort are provided in Supporting Information S1.
A total of 31827 SNPs in 2960 genes that are common to both GLGC and PennCAC are used in the current investigation. As an illustrative example, we chose to focus on the set of SNPs included on the IBC array which was specifically designed to provide SNP coverage in putative candidate CVD genes as well as emerging loci at the time of design . This facilitates (a) focus on a defined set of SNP within candidate loci and (b) direct comparison of findings across the two datasets that were examined. Representing the largest published lipids meta-analysis, GLGC significant SNPs are treated as the “gold-standard” to which application of MixMAP to PennCAC is compared. Further validation of the findings resulting from applying MixMAP to GLGC is not practical presently given the comprehensive inclusion of all prior genome wide association studies in GLGC, although this should be feasible in the near future with published and available larger datasets from on-going projects such as the GLGC Metabochip project . Instead, however, we use published literature on human data, animal/mouse models and cell biology to support our suggested novel discoveries.
Here we describe our interpretation of MixMAP findings in assigning genes to already established or novel loci in GLGC analysis. Briefly, we assigned MixMAP identified genes to an established GLGC locus (marked by SNP in GLGC Table 1 ) if LD values were for one or more SNPs in a MixMAP identified gene and a GLGC top SNP. Additionally, MixMAP identified genes that were in close physical proximity () to a GLGC top SNP but did not have strong LD () or were within a wider region with multiple candidate genes (e.g. HLA) were also assigned to the GLGC established LDL-C locus. The SNP Annotation and Proxy Search (SNAP)  web-based tool (http://www.broadinstitute.org/mpg/snap/) (for SNPs Kb apart) and the Genome-wide LInkage DisEquilibrium Repository and Search engine (GLIDERS)  (for SNPs apart) were used to determine pairwise LD between interrogated SNPs. In SNAP, we used the SNP dataset for the 1000 Genomes Pilot 1 [or HapMAP 3 (release 2) build 36 for SNPs not available in 1000 Genomes] and population panel “CEU” while in GLIDERS HapMAP 3 (release 2) build 36 was used.
Analysis of 31827 SNPs in 2960 genes in GLGC data identified genes in 26 loci based on a single SNP signal threshold of (Table 2). MixMAP detects genes within 21 of these 26 loci based on a Bonferroni corrected gene-level threshold of (Table 2; Figure 1; Table S1). Thus, of 26 LDL-C loci detected by the single SNP approach in GLGC GWAS data, 21 were also detected by MixMAP while 5 were not detected. A comparison of MixMAP detected and not detected GLGC genes/loci is of some interest to illustrate the strengths and weaknesses of the approach. First, coverage of the genes at loci detected by single SNP analysis but not detected by MixMAP is generally lower than those detected by both approaches (5 vs. 13 median number of SNPs respectively). Second, genes detected by a single SNP signal but not MixMAP tend to have a higher median SNP p-value (median = 0.006 for genes with greater than 2 SNPs) compared to genes detected by both approaches (median) suggesting that the distribution of p-values of SNPs within genes detected by both approaches is shifted downwards compared to that of genes detected by a single SNP only signal ( for two-sided Wilcoxon rank sum test of difference in medians). Importantly, from a biological perspective, almost all clinically important LDL-C genes/loci were detected (e.g., LDLR, APOB, APOE, HMGCR, PCSK9, LPA, SORT1, ABCG5/8, TRIB1, ABCA1, APOA5-A4-C3-A1 and CETP) while genes at loci not detected (LDLRAP1, ANGPTL3, HFE, HPR, TOP1) tended to be less well characterized functionally and clinically.
Twelve (12) additional loci are supported by MixMAP that are not detected using the single SNP signal threshold in the GLGC data (Table 3 and Figure 1). While the minimum -values within index genes at these loci are more moderate than those detected by the single SNP threshold, in that they do not reach genome-wide significance, the overall distribution of -values is lower than expected under no association (one-sided Wilcoxon signed rank test that the median of the within gene -values is greater , against the alternative that it is less than 0.50, .) The median -value of SNPs within these genes ranges from for a gene with SNPs to for a gene with relatively high coverage of SNPs. Several of the genes in this group (e.g. PPARG, PPARD and NPC1) are strongly implicated in lipid and lipoprotein metabolism in human and animal model systems, as indicated and referenced in the final column of Table 4 (See Supporting Information S2 for more detail).
Here we apply MixMAP to a small study (PennCAC) that generated IBC array SNP data. In terms of independent support, we treat published GLGC data as the gold standard to which PennCAC is compared. This analysis includes 31585 (of the original 31827 SNPs) in 2944 (of 2960) genes that remain after filtering out those with MAF and HWE. As might be anticipated for this small sample, no SNPs meet genome-wide significance () or even IBC array-wide significance (; an estimated threshold for independent SNP tests based on simulations ) for association with LDL-C using the single SNP signal approach. However, genes in loci contain SNPs that are significant at the level, a previously applied threshold for suggestive stage 1 evidence of association in IBC studies . Notably, in PennCAC, none of these loci coincide with significant MixMAP locus tests at the Bonferroni corrected threshold of ; however, interrogated genes (BUD13 and APOA5) in of these loci (APOA5-A4-C3-A1) do have single SNP signals for LDL-C at the Bonferroni corrected level in the GLGC study data . The same genes within this locus had significant MixMAP findings in GLGC data.
In PennCAC, an additional genes in distinct loci are supported by MixMAP based on the Bonferroni corrected threshold of (Table 4 and Figure 2). None of these loci have single SNP tests that reach even the suggestive threshold of in PennCAC; however, of these loci reach genome wide significance based on single SNP signals in GLGC (SORT1 and LPA). The same loci, as well as VPS13B had significant MixMAP findings in GLGC data. Furthermore, loci (IL1R2, VPS13B) have some support for modulation of lipids in animal models, as indicated in the final column of Table 3. Thus, in this PennCAC data-set, an illustrative example for small sample settings, MixMAP may add value to single SNP based testing in identifying loci for LDL-C and other complex traits.
In order to evaluate the performance of MixMAP relative to single SNP analysis, we investigate the ability to detect informative loci and the likelihood of false findings as functions of gene level effect size (measured by the shift parameter in a two-component Gaussian mixture distribution), the number of informative genes (measured by the number of genes with random effects arising from a non-zero mean normal prior) and coverage (measured by the proportion of observed SNPs included in the analysis). Details of the simulation approach are described in Materials and Methods below. Simulation results are reported in terms of: the true positive rate (TPR), defined as the proportion of true signal genes that are correctly identified; the false discovery rate (FDR), given by the proportion of selected genes that are not associated with the trait; and the false positive rate (FPR), defined as the proportion of truly uninformative genes that are incorrectly selected as significant. Explicit definitions of TPR, FDR and FPR are given in Table 1. In all simulation scenarios, the results of applying MixMAP as well as the single SNP approach are presented to illustrate the potential gains associated with using MixMAP as a complementary strategy.
Figures 3 illustrates the estimated TPR (Sensitivity), FDR and FPR (1-Specificity) for shift parameters ranging from to , where estimates are based on 500 simulations per condition. The left panel of this figure reports results when the true number of informative gene is , the number detected in the GLGC data, while the right panel considers twice as many () informative genes. The range of to informative genes was selected because it is generally consistent with the number of significantly associated loci in prior reports of complex disease traits. Notably for the GLGC data, the mean empirical Bayes (EB) (shrinkage) estimate for the detected genes was . This is consistent with a shift parameter of approximately 2.3, as reported in Figure 4; however, a broader range of shift parameters is presented in the simulations studies for illustration and generalizability.
As expected, the TPR increases with increasing shift parameter values, with MixMAP consistently higher than the single-SNP based approach for more moderate shift parameters (). For both the single-SNP and MixMAP approaches the rate of increase in the TPR (as a function of the shift parameter) is slower when the number of informative genes is 96 versus 48, though the difference is more pronounced in the context of MixMAP. This result is consistent with underlying statistical theory, as the estimated variance parameter for the single normal prior distribution of the random effects will be larger as more genes arise from an alternative distribution with non-zero mean. As a result, the corresponding gene-level prediction intervals will be wider and it will be more difficult to detect informative genes. This result suggests that as the number of truly informative genes increases, the added contribution of MixMAP over single SNP analysis is smaller; however, with as many as informative genes, MixMAP is detecting a greater percentage of truly informative genes than the single SNP approach for shift parameters less than . While the FDR is consistently low for the single-SNP based approach regardless of the shift parameter, the FDR for MixMAP is relatively high () for shift parameters less than . Notably, the average number of detected genes (across the 500 simulations) that are truly uninformative (FPs) is relatively constant across shift parameter values, ranging between to with a median of . At the same time, the FDR tends to decrease because the total number of detected genes (TPs+FPs) is increasing as the shift parameter gets larger. Finally, the FPR is consistently small for both the single-SNP based and MixMAP approaches across all shift parameter values with averages (assuming 48 truly informative genes) of and , respectively.
Figure 5 reports the estimated TPR, FDR and FPR when only a subset of the observed SNPs within each gene is included in the analysis. The proportion of included SNPs (given on the x-axis) ranges from to and each column in this figure represents a different value of the shift parameter (2.0, 3.0 or 4.0 for columns 1, 2 and 3 respectively). As expected for all values of the shift parameter, the TPR increases dramatically for MixMAP as the proportion of SNPs within a gene approaches . This follows from the dependency of the prediction interval for a gene level random effect on the number of SNPs within the corresponding gene. Further the MixMAP FDR is relatively high () for low proportions () of included SNPs within each gene, and this is more pronounced for small shift parameters. Notably, while the FDR is high, the average number of detected genes (across the 500 simulations) that are truly uninformative (FPs) is relatively constant across the proportion of SNPs included, ranging between to with a median of for shift parameter of . Finally, the FPR is consistently small for both the MixMAP and single-SNP approaches, with estimated values similar to those reported above with complete SNP data. These results suggests that the contribution of MixMAP will be more pronounced when genes are adequately covered, and application of MixMAP in the context of substantial missing information on SNPs can result in a higher proportion of false discoveries (relative to the total number of discoveries), though the absolute number of false discoveries remains constant.
Despite advances in the genomics of complex traits, only a portion of heritability for common human diseases has been elucidated. To date, most common variant discovery approaches have relied on tests of disease association using one SNP at a time. Methods that leverage existing datasets and exploit information across multiple SNPs within a gene region are likely to yield additional information regarding the locus association with traits. We developed a novel strategy, MixMAP, that relies on well-vetted statistical principles and draws from the vast array of summary data now available from genetic association studies, to test formally for locus-level association. The primary inputs required for this approach are single SNP level p-values for tests of trait association and mapping of SNPs to locus regions while the output is locus level estimates and tests of association. Application of MixMAP to SNP summary data for a pre-defined set of genes within the GLGC meta-analysis of association with LDL-C suggest that MixMAP can provide substantial value in discovery that is complementary to single SNP testing approaches in identifying novel loci for LDL-C. In addition, MixMAP analysis of PennCAC IBC array and LDL-C data support its application in combination with traditional SNP testing to enhance the power of discovery in small dataset settings. Thus, MixMAP provides a novel strategy, based on established statistical principles, for exploiting existing and emerging genomic data to provide advances in our understanding of complex human diseases.
Over the past decade sequencing of the human genome, definition of common SNP variation in human population and advances in genotyping technology have provided the possibility to discover common genetic contributions to complex traits in human. Indeed, very large scale applications of genome SNP scans in humans combined with rigorous statistical correction for multiple testing has led to an explosion of novel validated genomic discoveries for human diseases with exciting progress in functional genomics as well as promise for novel therapeutics and disease prediction. Despite this the majority of heritability for most complex traits remains to be discovered. Current statistical approaches for testing single SNP associations with disease are designed to protect against excess false positives but may be excessively conservative. Further, single SNP approaches to analysis do not draw strength from information gained by assessing simultaneously trends of association across a locus. These observations suggest that false negatives are a significant feature of existing association analysis and that additional genomic discovery should be possible in existing data if appropriate statistical methodologies are applied. Indeed, recent research suggests that common variants with individual level effects that are too small to be considered statistically significant using stringent significance thresholds account for a substantial proportion of this missing heritability for complex traits . However, differentiating the true signals within the vast amount of SNP data with moderate p-values remains an unsolved problem.
We chose to analyze genetic contributors to LDL-C for several reasons. First, LDL-C is an important complex trait that is causal for a substantial portion of CVD death and morbidity in our society. Second, LDL-C has a well described heritability and large rigorously performed meta-analyses have been performed and summary data are publicly available (GLGC). Third, although many loci for LDL-C have been identified through association studies at , only a modest portion of its heritability (approximately of genetic variability ) has been defined. Fourth, the basic biology of plasma lipids and LDL-C has been extensively studied in animal models and cell systems providing some additional mechanistic reference for any novel discoveries we might make. We hypothesized that we would identify novel loci for LDL-C, beyond the existing single SNP-based discoveries, through application of MixMAP in the large GLGC meta-analysis summary data. As an informative example, we chose to focus on the set of SNPs for CVD candidate genes included on the ITMAT-Broad-CARe (IBC) SNP array which was designed to provide dense SNP coverage in putative candidate CVD genes as well as some coverage of emerging loci at the time of its design . This approach allowed us to focus on a defined set of SNPs within candidate loci and to perform direct comparison of findings for this subset of SNPs within the GLGC data-set to those in the smaller PennCAC sample application.
In GLGC data, MixMAP confirmed association for over of the loci identified through single SNP testing of the 31827 SNPs in 2960 genes examined. Failure to detect more of the loci established by single SNP testing should not be surprising because MixMAP loses information for extremely low individual SNP p-values and is not designed for finding association when SNP coverage of a gene region is poor, as is the case for some loci that reached significance in the GLGC. Further, from a biological perspective, almost all clinically important LDL-C genes/loci were detected by MixMAP (e.g., LDLR, APOB, APOE, HMGCR, PCSK9, LPA, SORT1, ABCG5/8, TRIB1, ABCA1, APOA5-A4-C3-A1 and CETP). MixMAP, however, did provide evidence for new loci (corresponding to interrogated genes) for LDL-C in GLGC data that did not reach genome wide significance in single SNP testing. This may be an under estimate because we applied conservative criteria for our selection of novel loci (greater than 500kb from known GLGC locus, pairwise with top SNP at GLGC established LDL-C locus, and outside of region with multiple candidate genes). For example, interrogated SNPs in C2, which MixMAP identifies as a gene associated with LDL-C, have with the top GLGC SNP at the HLA locus in Teslovich et al. .
A more detailed description of the 12 genes/loci detected by MixMAP alone is provided in Supplementary Materials. For many of these 12 loci (FN1, UGT1A1, PPARG, GAB2 and APOH) the GLGC single SNP test p-value provided suggestive evidence of association () and published data in mice and human support specific biological processes and plausible mechanisms of association with LDL-C for some of the index genes (UGT1A1, PPARG and APOH) at these loci –. Notably, a recent meta-analysis of IBC array data for plasma lipids across 66,240 individuals also supports an association of APOH with LDL-C and suggests that UGT1A1 is a locus for total cholesterol levels . For some MixMAP significant loci with suggestive GLGC single SNP tests, there is no or limited published biology or mechanism for association with LDL-C (e.g. FN1 and GAB2). On the other hand, some loci that are significant by MixMAP have quite modest statistical support in GLGC single SNP analysis, but have strong published data supporting mechanisms by which genes at the locus may modulate LDL-C (e.g NPC1 and PPARD) –. Finally, a few MixMAP loci have neither suggestive single SNP support from GLGC nor reported biological plausibility for gene-lipid associations (e.g. PKN2 and CDK) and such loci require further focus and validation. Overall, these data support the utility of MixMAP, when used in combination with traditional single SNP testing, in discovery of true loci for LDL-C and other complex traits particularly.
A specific challenge in the genomics of complex traits is identifying loci for such a trait when power is low due to limited availability of human data. We chose to illustrate this issue in a small sample (PennCAC, n=2096) using LDL-C as an example in part because the large GLGC dataset for LDL-C provides an external reference for any MixMAP findings. In PennCAC, no individual SNPs meet criteria for association with LDL-C using the conservative genome-wide Bonferroni correction () or the less conservative IBC array-wide Bonferroni correction (). At a less stringent, suggestive single SNP criteria (), loci (represented by genes) are identified. At one of these loci, interrogated genes (APOA5 and BUD13) contain SNPs with genome-wide significant signals in the independent GLGC dataset. As expected at this less conservative threshold, however, most SNPs lack supporting signals in GLGC data and lack supporting biology for genes at the locus (Table 3) suggesting that several may be false positives. Using MixMAP genes, representing independent loci, were suggested for LDL-C. Of these, (SORT1 and LPA) have genome-wide significant signals in the independent GLGC data and one (VPS13B) had a significant MixMAP signal in GLGC data. Furthermore, loci (IL1R2 and VP13B) have some support for modulation of lipids in animal models –. Overall, these PennCAC LDL-C analyses suggest that application of MixMAP in small sample settings may provide complementary value to single SNP tests and other strategies to maximize genetic inference in settings where sample size is limited. Although in these small sample settings false positives will remain a challenge, MixMAP should enhance findings for prioritization and further follow-up.
We recognize that independent replication of findings is essential for complete validation of novel findings in genetic studies. Because the GLGC data represent the largest published lipids GWAS meta-analysis to date, we believe a comprehensive replication for LDL-C beyond these GLGC data is not possible at the current time. However, we will pursue this for lipid genes/loci in additional GLGC data when these data become available (e.g., Metabochip project data expected 2013 ). We also acknowledge that for common SNP variation, a single gene often can not specifically be assigned to the disease-associated variant. Further, simple proximity to a variant and even incorporation of expression QTL knowledge are not always correct in selecting causal genes. This problem can lead to incorrect assumptions of causal genes and raise concerns for validity of gene-based inference. However, this limitation is not unique to our illustration of MixMAP and is common to current gene and pathway analyses leveraging common SNP datasets (e.g. ). The challenge can be addressed in part by leveraging the maximum amount of linkage disequilibrium, eQTL, fine mapping and biological data when assigning genes to the associated SNPs. In the present investigation, we use gene as the cluster to which SNPs belong, though MixMAP is not limited by this specification. Importantly, the user can employ alternative and newly evolved classifications, as the primary input to the MixMAP algorithm.
The results of the simulation study further support the application of MixMAP as a complementary strategy to single-SNP based testing, particularly in the context of moderate gene level effects and adequate SNP coverage. Our on-going research is exploring calibrating the variance coefficient in the prediction interval, as an alternative to using , to obtain desired control of the FDR in specific well-defined settings. Because a first stage ranking of -values is applied prior to inverse normally transforming the data for model fitting, the implications of using p-values from a single cohort study (PennCAC) versus a meta-analysis (GLGC) are limited to the varying degrees of precision in each setting. That is, the full range of the quantitative data, and specifically the fact that p-values from a meta-analysis tend to be substantially smaller than those from a single cohort study, is not being incorporated into the analysis presented herein. We expect additional knowledge can be gained through a mixture modeling extension of MixMAP that can accommodate the quantitative nature of the summary data, and this is currently under investigation. The present investigation is based on common variants, and while incorporating the results of rare variant analysis poses an additional challenge as these variants tend to be grouped a priori for analysis, such an extension would also likely be informative.
Further extensions of MixMAP would also allow application to gene set and/or pathway-based analysis of association data. Specifically, through inclusion of multiple nested random effects, the MixMAP framework could be applied using both locus level and pathway information simultaneously. Through fully parametric modeling, this may offer advantages over gene set enrichment analysis, which similarly involves a first stage rank ordering . This extension of MixMAP would be notably distinct from the hierarchical modeling approach of  that similarly includes random gene specific effects, but separately models each gene set and focuses testing on fixed intercepts representing pathway effects rather than latent variables. Additional future work includes a specific evaluation of the influence of linkage disequilibrium, minor allele frequencies, gene size and numbers of recombinant hotspots as potential covariates in the models, as well as comprehensive evaluation of the complex statistical power considerations across a range of applications and conditions, including candidate gene studies, GWAS, pathway analysis and partial or whole-exome sequencing studies. Additional characterization of MixMAP may facilitate applications to summary findings from Metabochip and exome sequencing, as well in interrogation of gene sets and pathways utilizing such data. In conclusion, the approach we have described is intended to complement single SNP analysis and should provide a useful tool to potentiate existing summary data and reveal important novel loci, pathways and causal factors for complex diseases at little additional cost.
As a consequence of the LD structure within genetic loci, we expect SNP level -values, corresponding to single tests of association, to be potentially more similar within a gene than across genes, regardless of the level of association. Thus, a common statistical modeling framework for correlated data, the mixed effects model –, is a natural analytic framework to consider for this setting. The application of MixMAP presented herein is at a gene level, and thus the term “gene” is used throughout this section; however, we note that as additional, locus-level annotations become available, these can replace or enhance the gene level classifications.
We begin by transforming the SNP level analysis or meta-analysis -values, which will serve as the outcomes in our model, to normal variates in order to meet model assumptions. This is achieved by: (a) applying a simple rank transformation, to ensure uniformity over the interval from to ; and (b) applying an inverse normal (probit) transformation to normalize the data. For (a), the rank of the th SNP is set equal to where is the ordered ranking across all SNPs and is the total number of SNPs under study. For (b), we let where is the cumulative density function of a standard normally distributed random variable.
The first step of MixMAP is to fit a mixed effects model to appropriately transformed -values with gene-specific random intercept terms and fixed effects for any relevant covariates. Formally this model is given by:
where , is the transformed -value for the th SNP within gene , is a matrix with th row equal to a vector of SNP or gene level covariates, and , where is the total number of genes and is the number of SNPs in gene . Further, we let be an vector of 's, is the random effect of gene , and is independent of . Finally, is the corresponding vector of fixed effects. In the example presented in this manuscript, the term additionally reduces to an overall intercept, but we retain this for generalizability. In this model, represents a latent (unobservable) effect of gene on the corresponding transformed -values. Notably, as a result of the transformation described above, small -values correspond to large negative values of . Thus, values of that are less than would indicate a gene level effect.
Based on this model formulation, the best linear unbiased predictor of the random effect for gene is given by:
where where . The empirical Bayes estimate of , denoted , is calculated by replacing and with corresponding REML estimates. A measure of dispersion for this quantity in the intercept only model is given by:
and returned by the lmer() function in the R lme4 package (http://cran.r-project.org/web/packages/lme4/index.html). This is related in expectation to as described in , Chapter 7. Notably, the novel gene and locus findings reported in this manuscript were not sensitive to choice of prediction variance; however, if were applied in place of the measure returned by lmer(), then three genes (namely HAVCR2, HLA-DRA and LPAL2 in the TIMD4, HLA and LPA loci, respectively) would not be MixMAP significant.
In general, we are interested in testing the null hypothesis of no association between a given gene and the trait. To this aim, we construct a one-sided prediction interval for the true gene-level effect, given by for gene , with the upper limit defined as:
where is the ()-quantile from a standard normal distribution, is equal to evaluated at RML estimates of and , and is a pre-defined significance threshold. A Bonferroni level threshold given by where is the number of genes under study is suggested and applied in this manuscript. If the upper limit of the gene-level prediction interval is less than , then we conclude that the corresponding locus is significantly associated with the trait.
-values corresponding to tests of association with LDL-C were generated as follows, and according to the approach described in : (1) Regress LDL on age, age, gender and the first principal components derived using all available SNPs; (2) Calculate the residuals from this model fit; (3) Fit a separate simple linear regression for each SNP (coded as ordinal 0, 1, 2 variables) with the residuals as the outcome; and (4) Record the t-test statistic and corresponding p-value within each model for the test that the coefficient of the SNP in the linear regression was equal to . For the PennCAC data, SNP level -values were generated within Caucasians according to this same algorithm. For GLGC data, reported meta-analysis -values that were generated in the same manner and then meta-analyzed were used in analysis. Transformations of these -values and subsequent application of MixMAP proceeded as described above.
For all simulation studies, random gene level effects, for , are simulated according to a two-component Gaussian mixture distribution with elements arising from a and the remaining elements arising from a distribution, where is the total number of genes under study and is the number of informative genes. SNP level z-scores are then generated according to the model where . To begin for all simulations, the numbers of SNPs within each gene are set equal to the observed counts for the 31825 SNPs within the 2960 genes in the GLGC and PennCAC IBC subset. The median number of SNPs per gene is , the mean is and the range is to . The first and third quartiles are equal to and , respectively. Notably, due to the LD structure within genes, we expect -values to be correlated within these regions even under the complete null of no association between all genes and the trait under study. As a result, the gene level random effects, given by in the model above, are not identically equal to under this null. That is, even uninformative genes, whose latent effects are assumed to arise from a mean distribution, will have corresponding non-zero effects.
Summary of suggested novel genes/loci identified in GLGC data by MixMAP.
Empirical Bayes (EB) estimates and corresponding prediction intervals for MixMAP supported genes in GLGC. The Empirical Bayes estimate and corresponding prediction variance for each gene are used in the construction of an associated one-sided prediction interval. 1The upper limit of this interval is reported using a Bonferroni corrected α=0.05/(2960). An upper limit that is less than 0 implies the gene's random effect (on the inverse normal transformed ranked p-values) is significantly less than 0. Since negative transformed p-values that are large in absolute value correspond to small p-values on the original scale, this implies the gene as a whole has a significant effect on LDL-C. *Highlighted rows correspond to novel loci not identified using single SNP analysis.
Empirical Bayes (EB) estimates and corresponding prediction intervals for MixMAP supported genes in PennCAC. The Empirical Bayes estimate and corresponding prediction variance for each gene are used in the construction of an associated one-sided prediction interval. 1The upper limit of this interval is reported using a Bonferroni corrected α=0.05/(2944). An upper limit that is less than 0 implies the gene's random effect (on the inverse normal transformed ranked p-values) is significantly less than 0. Since negative transformed p-values that are large in absolute value correspond to small p-values on the original scale, this implies the gene as a whole has a significant effect on LDL-C. None of these loci were identified using single SNP analysis.