Patent Application
20150356243
The present invention relates to processes, systems and methods for estimating the effects of genetic polymorphisms associated with traits and diseases, based on distributions of observed effects across multiple loci. In particular, the present invention provides systems and methods for analyzing genetic variant data including estimating the proportion of polymorphisms truly associated with the phenotypes of interest, the probability that a given polymorphism has a true association with the phenotypes of interest, and the predicted effect size of a given genetic variant in independent de novo samples given effect size distributions in observed samples. The present invention also relates to using the described systems and methods and use of genetic polymorphisms across a plurality of loci and a plurality of phenotypes to diagnose, characterize, optimize treatment and predict diseases and traits.
Many devastating human diseases are heritable, including many of the largest health care burden today, including cardiovascular diseases, brain disorders, rheumatologic and immunological disorders. However, only a small fraction of genetic variance has been identified, even after using large genome-wide association studies (GWAS). Several lines of evidence support the existence of numerous small genetic effects that cannot be detected with traditional GWAS analyses.
Converging evidence suggest that complex human phenotypes are influenced by numerous genes each with small effects. Though thousands of single nucleotide polymorphisms (SNPs) have been identified by genome-wide association studies (GWAS), these SNPs fail to explain a large proportion of the heritability of most complex phenotypes studied, often referred to as the “missing heritability” problem. Recent findings indicate that GWAS have the potential to explain a greater proportion of the heritability of common complex phenotypes, and more SNPs are likely to be identified in larger samples. Due to the polygenic architecture of most complex traits and disorders, a large number of SNPs are likely to have associations too weak to be identified with the currently available sample sizes.
New analytical methods are needed to reliably identify a larger proportion of SNPs associated with complex diseases and phenotypes, since recruitment and genotyping of new samples are expensive.
The present invention relates to processes, systems and methods for estimating the effects of genetic polymorphisms associated with traits and diseases, based on distributions of observed effects across multiple loci. In particular, the present invention provides systems and methods for analyzing genetic variant data including estimating the proportion of polymorphisms truly associated with the phenotypes of interest, the probability that a given polymorphism has a true association with the phenotypes of interest, and the predicted effect size of a given genetic variant in independent de novo samples given effect, size distributions in observed samples. The present invention also relates to using the described systems and methods and use of genetic polymorphisms across a plurality of loci and a plurality of phenotypes to diagnose, characterize, optimize treatment and predict diseases and traits.
For example, in some embodiments the present invention provides a computer implemented process of identifying polymorphisms associated with a specific condition, comprising at least one of: a) inputting polymorphism information for a plurality of gene variants (e.g., single nucleotide polymorphisms (SNP)); b) assigning a linkage disequilibrium (LD) score to each SNP; c) testing each gene variant for enrichment using scores derived from conditional distribution analysis (e.g., Q-Q plots); d) assigning a ranking (e.g., false discovery rate (FDR) or local false discovery rate) to each gene variant using unconditional and conditional distributions; e) performing a Bayesian, resampling, or likelihood-based analysis on a combination of all or some enriching factors; f) applying a regression model to combine information; and g) identifying or quantifying the probability that the gene variants are associated with the condition. In some embodiments, identifying comprises listing identified gene variants in a priority order. In some embodiments, the LD assigns each of the gene variants to a functional category. In some embodiments, the Q-Q score provides a true discovery rate and a FDR for each SNP. In some embodiments, the FDR for a specific gene variant is defined as the nominal p-value divided by the empirical quantile. In some embodiments, gene variants with FDRs less than a threshold value (e.g., 0.01) are defined as associated with the condition. In some embodiments, empirical quantiles are plotted as Q-Q plots. In some embodiments, Q-Q plots identify pleiotropic enrichment. In some embodiments, polymorphism information is obtained from at least 2 subjects. In some embodiments, polymorphism information comprises at least 1000, 5000, or 10000 or more individual gene variants. In some embodiments, gene variants are intergenic. In some embodiments, the method further comprises the step of plotting FDRs within an LD block in relation to their chromosomal location. In some embodiments, the condition is, for example, a disease, a trait, a response to a particular therapeutic agent, or a prognosis, although other conditions are specifically contemplated.
In some embodiments, distributions of gene variant, effect sizes for a given trait or disease are used to determine Bayesian posterior effect sizes across a plurality of polymorphisms. In some embodiments, Bayesian posterior effect sizes are computed across a plurality of diseases or traits simultaneously. In some embodiments, prior information regarding genes, functional roles of SNPs, LD scores, or other covariates is used to improve estimates of Bayesian posterior effect sizes. In some embodiments, distributions of Bayesian posterior effect size for one or more diseases or traits is used to identify genetic loci associated with a disease or trait. In some embodiments, Bayesian posterior effect sizes in one or more diseases or traits is used to explain observed variance in a disease or trait. In some embodiments, Bayesian posterior effect size distributions for one or more diseases or traits is used to compute a polygenic risk score for the a disease or trait. In some embodiments, the polygenic risk score for a disease or trait is used to predict the risk of an individual having a disease or trait. In further embodiments, the predicted risk of an individual have the disease or trait includes confidence intervals indicating the degree of precision of the estimated risk. In some embodiments, distributions of Bayesian posterior effect sizes is used to produce estimates of power for identifying polymorphisms associated with a disease or trait in genetic studies for a given study sample size.
In further embodiments, the present provides a plurality of gene variants identified by the process described herein, wherein the plurality of gene variants are associated with a specific condition.
In yet other embodiments, the present invention provides a method, comprising: a) identifying a plurality of gene variants from a subject associated with a given condition using the process described herein; and b) characterizing one or more conditions in the subject based on the plurality of gene variants. In some embodiments, the method further comprises the step of providing a diagnosis or a prognosis to the subject. In some embodiments, the method further comprises the step of determining a treatment course of action based on the characterizing (e.g., choosing a therapeutic agent and/or choosing a dosage of a therapeutic agent.
In some embodiments, the present invention provides computer implemented processes and methods calculating polygenic personalized risk scores associated with a specific condition, comprising: computing gene variant, (e.g., single nucleotide polymorphisms (SNP)) posterior effect sizes (e.g. by randomly dividing subjects from a given group into disjoint training and replication subsamples); calculating sample mean replication effect sizes conditional on training effect sizes; and determining a polygenic risk score based on the effect sizes. In some embodiments, the polygenic risk score is computed as a linear or nonlinear function of the estimated statistical parameters. In some embodiments, the linear or nonlinear function of the estimated statistical parameters includes per gene variant allele effect size mean and/or estimates of variability. In some embodiments, computing comprises linear weighting of each gene variant by its estimated posterior effect size divided by its estimated posterior variance. In some embodiments, the process further comprises the step of obtaining maximal correlation of genetic risk scores with phenotypes in de novo subject samples by obtaining posterior effect size estimates for each SNP modulated by genie annotations and/or strength of association with pleiotropic phenotypes. In some embodiments, the posterior effect sizes for each gene variant are multiplied by the corresponding gene variant values for a de novo subject and added together to calculate an overall risk score for the condition or the posterior effect sizes for each SNP are scaled by dividing by a measure of its variability before computing the polygenic risk score. In some embodiments, gene variant effect sizes below a given threshold are deleted before computing polygenic risk scores. In some embodiments, the comprises subjects from a single study or collection of studies. In some embodiments, the polygenic personalized risk scores summarize patient-level genomic variation as a single score per subject, summed over assayed gene variants. In some embodiments, the polygenic personalized risk score includes other biomarkers of the condition, for example, including but not limited to, age, gender, family history, or results of diagnostic testing. In some embodiments, the process further comprises the step of predicting the likelihood of an offprising of two parents developing the condition. In some embodiments, predicting comprises the step of randomly simulating multiple offspring and estimating polygenic risk scores for each simulated offspring and using the scores across offspring to predict the likelihood of said offspring developing the condition.
Additional embodiments will be apparent to persons skilled in the relevant art based on the teachings contained herein.
To facilitate an understanding of the present invention, a number of terms and phrases are defined below:
As used herein, the term “sensitivity” is defined as a statistical measure of performance of an assay (e.g., method, test), calculated by dividing the number of true positives by the sum of the true positives and the false negatives.
As used herein, the term “specificity” is defined as a statistical measure of performance of an assay (e.g., method, test), calculated by dividing the number of true negatives by the sum of true negatives and false positives.
As used herein, the term “informative” or “informativeness” refers to a quality of a marker or panel of markers, and specifically to the likelihood of finding a marker (or panel of markers) in a positive sample.
As used herein, the term “amplicon” refers to a nucleic acid generated using one or more primers (e.g., two primers). The amplicon is typically single-stranded DNA (e.g., the result of asymmetric amplification), however, it may be RNA or dsDNA.
The term “amplifying” or “amplification” in the context of nucleic acids refers to the production of multiple copies of a polynucleotide, or a portion of the polynucleotide, typically starting from a small amount of the polynucleotide (e.g., a single polynucleotide molecule), where the amplification products or amplicons are generally detectable.
As used herein, the term “primer” refers to an oligonucleotide, whether occurring naturally as in a purified restriction digest or produced synthetically, that is capable of acting as a point of initiation of synthesis when placed under conditions in which synthesis of a primer extension product that is complementary to a nucleic acid strand is induced (e.g., in the presence of nucleotides and an inducing agent such as a biocatalyst (e.g., a DNA polymerase or the like) and at a suitable temperature and pH). The primer is typically single stranded for maximum efficiency in amplification, but may alternatively be double stranded. If double stranded, the primer is generally first treated to separate its strands before being used to prepare extension products, in some embodiments, the primer is an oligodeoxyribonucleotide. The primer is sufficiently long to prime the synthesis of extension products in the presence of the inducing agent. The exact lengths of the primers will depend on many factors, including temperature, source of primer and the use of the method. In certain embodiments, the primer is a capture primer.
A “sequence” of a biopolymer refers to the order and identity of monomer units (e.g., nucleotides, etc.) in the biopolymer. The sequence (e.g., base sequence) of a nucleic acid is typically read in the 5′ to 3′ direction.
As used herein, the term “subject” refers to any animal (e.g., a mammal), including, but not limited to, humans, non-human primates, rodents, and the like, which is to be the recipient of a particular treatment. Typically, the terms “subject” and “patient” are used interchangeably herein in reference to a human subject.
As used herein, the term “non-human animals” refers to all non-human animals including, but are not limited to, vertebrates such as rodents, non-human primates, ovines, bovines, ruminants, lagomorphs, porcines, caprines, equines, canines, felines, aves, etc.
The term “locus” as used herein refers to a nucleic acid sequence on a chromosome or on a linkage map and includes the coding sequence as well as 5′ and 3′ sequences involved in regulation of the gene.
In the present context the term “psychiatric disease” refers to brain disorders with a psychological or behavioral pattern that occurs in an individual and cause distress or disability that is not expected as part of normal development or culture, including symptoms related to behavior, emotion, cognition, perception, thought disorder. Non-limiting examples of psychiatric diseases are schizophrenia, other psychotic disorders, depression, bipolar disorder, depression, anxiety, OCD, Personality disorders, PTSD, Alzheimer's disease, eating disorders, child psychiatry disorders.
In the present context the term “neurological disease” refers to brain disorders involving the central, peripheral, and autonomic nervous systems, including their coverings, blood vessels, and all effector tissue, such as muscle, with primarily symptoms related to movement, but often other symptoms in addition, such as memory impairment, fatigue, pain, sensitivity abnormalities. Non-limiting examples of neurological diseases are stroke, epilepsy, neurodegenerative disorders, headache, multiple sclerosis.
As used herein, the term “gene variant” refers to any change in nucleotide sequence or dosage within a gene relative to the native or wild type sequences or copy number. Examples include, but are not limited to, mutations, single nucleotide polymorphisms (SNPs), copy number variants, deletions, inversions, duplications, splice variants, or haplotypes.
In the present, context the term “genotype information” refers information which can be obtained from the genome of an individual. Thus, genotype information may only be information from, part of the whole genome of the person. Non-limiting examples of genotype information which can be used in the present methods include SNPs (single-nucleotide polymorphisms), copy number variants (CNV), deletions, inversions, duplications, sequence variants, haplotypes. Preferably the genotype information obtained from a person are SNP's. Thus, in the present description, genotype information is used as a generic term for various genetic polymorphisms.
In the present context the phrase “SNP dose” refers to the number of times a specific SNP is present. Thus, for an individual the SNP dose can be 0, 1 or 2, meaning that a SNP dose of 0 means the specific SNP is not present in any of the two alleles, whereas a SNP dose of 1 means the SNP is present in one of the two alleles and a SNP dose of 2 means that the SNP is present on both alleles.
The present invention relates to processes, systems and methods for estimating the effects of genetic polymorphisms associated with traits and diseases, based on distributions of observed effects across multiple loci. In particular, the present invention provides systems and methods for analyzing genetic variant data including estimating the proportion of polymorphisms truly associated with the phenotypes of interest, the probability that a given polymorphism has a true association with the phenotypes of interest, and the predicted effect size of a given genetic variant in independent de novo samples given effect size distributions in observed samples. The present invention also relates to using the described systems and methods and use of genetic polymorphisms across a plurality of loci and a plurality of phenotypes to diagnose, characterize, optimize treatment and predict diseases and traits.
Embodiments of the present invention provide processes, systems, and methods (e.g., computer implemented) for analysis of gene variant data and characterization of conditions. The below description is exemplified with SNPs. However, the systems and methods described herein find use in the analysis of any type of gene variant. Examples of gene variants include, but are not limited to, mutations, single nucleotide polymorphisms (SNPs), copy number variants, deletions, inversions, duplications, splice variants, or haplotypes.
In the present study the power of GWAS data was leveraged to demonstrate how GWAS from disorders can improve discovery of novel susceptibility loci. Using standard GWAS analytical methods, only one significant locus was identified. By applying the stratified FDR method (Yoo et al, (2009) BMC Proc 3 Suppl 7: S103; Sun et al., (2006) Genet Epidemiol 30:519-530), an additional 7 loci (2 in bipolar disorder, 5 in schizophrenia) were found. Combining the independent schizophrenia and bipolar disorder GWAS samples, a total of 58 loci were identified in schizophrenia and 35 in bipolar disorders, with FDR<0.05 as a threshold. These results demonstrate the feasibility of using a cost-effective, pleiotropy-informed stratified FDR approach to discover common variants in schizophrenia and bipolar disorders.
The current statistical framework is based on the fact that SNPs are not interchangeable. Rather, a SNP with effects in two associated phenotypes has a higher probability of being true nonnulls, and hence also a higher probability of being replicated in independent studies. A conditional FDR approach was developed for GWAS summary statistics, adapting stratification methods originally used for linkage analysis and microarray expression data (Yoo et al, (2009) BMC Proc 3 Suppl 7: S103; Sun et al., (2006) Genet Epidemiol 30:519-530). Decreased conditional FDR (equivalently, increased conditional TDR) for a given nominal p-value increases power to detect true non-null effects. Increased conditional TDR is directly related to increased replication effect sizes and replication rates in de novo samples. Using this stratified approach, it was possible to increase power to detect true non-null signals in independent studies for given nominal p-values cut-offs. Equivalently, in the stratified approach the FDR can be used to control FDR at a given level while increasing power to discover non-null SNPs over approaches that treat all SNPs as interchangeable (Craiu R V, Sun L (2008) Statistica Sinica 18: 861-879). A conjunction FDR approach was developed to investigate which SNPs are pleiotropic. SNPs that exceed a stringent, conjunction FDR threshold are highly likely to be non-null in two phenotypes simultaneously.
The current findings of polygenic enrichment indicate that genetic pleiotropy is important in severe mental disorders. However, the datasets utilized herein are exemplary. The present disclosure is not limited to a particular condition or disorder. By using a stratified FDR approach, it was possible to leverage the overlapping polygenetic architecture to identify more of the specific SNPs involved. The current approach identified 58 loci in schizophrenia compared to 7 in the original publication. In bipolar disorder, the added power from schizophrenia GWAS identified 35 loci compared to two loci in the original study. It is important to note that this improvement in gene discovery was obtained despite the much smaller number of controls in the current analyses because the original analyses of the two disorders used largely overlapping control samples. Since 1KGP data was used to calculate LD structure, the number of loci can vary somewhat compared to the original analysis. For both disorders, most of the current findings were borderline significant in the original GWAS mega-analysis, or identified in other GWAS of partly overlapping samples, such as TRANK1 and SYNE1.
The current findings provide genes and polymorphisms related to bipolar disorder and schizophrenia. However, the processes, systems, and methods described herein find use in the characterization of a variety of disorder and conditions.
In some embodiments, the present invention provides processes, systems, and methods for analyzing gene variant data, identifying gene variants useful for characterizing and diagnosing conditions and diseases. In some embodiments, the process comprises, a computer implemented process, system, or method of identifying polymorphisms associated with a specific condition, comprising at least one of: a) inputting polymorphism information for a plurality of gene variants (e.g., single nucleotide polymorphisms (SNPs)0: b) assigning a linkage disequilibrium (LD) score to each gene variant; c) testing each SNP for enrichment using a Q-Q score; d) assigning a FDR to each gene variant using a look up table; e) performing a baysesian analysis on a combination all enriching factors; f) applying a regression model to combine information; and g) identifying gene variants associated with the condition. In some embodiments, identifying comprises listing identified SNPs in a priority order. In some embodiments, the LD assigns each of the gene variants to a functional category. In some embodiments, the Q-Q score provides a true discovery rate and a FDR for each gene variant. In some embodiments, the FDR for a specific gene variant is defined as the nominal p-value divided by the empirical quantile. In some embodiments, gene variants with false discovery rates less than 0.01 are defined as associated with the condition. In some embodiments, Q-Q scores are plotted as Q-Q plots. In some embodiments, Q-Q plots identify pleiotropic enrichment. In some embodiments, polymorphism information is obtained from at least 2 subjects. In some embodiments, polymorphism information comprises at least 1000, 5000, or 10,000 or more individual SNPs. In some embodiments, gene variants are intergenic. In some embodiments, the method further comprises the step of plotting false discovery rates within a LD block in relation of their chromosomal location. In some embodiments, the condition is, for example, a disease, a trait, a response to a particular therapeutic agent, or a prognosis, although other conditions are specifically contemplated.
In some embodiments, systems and methods utilize the following steps as illustrated in
1) The first step is to input the GWAS data of a particular train or disease as one data file or individual chip/sequence data. The data file includes the p-values (the significance of association with disease) for each SNPs from the GWAS (this can be original chipped SNPs or imputed SNPs). In some embodiments, raw data (e.g., unthresholded SNP list) is used.
2) Each SNPs is then annotated to the most recent catalogue of the human genome, such as 1000 genomes project (1KGP) for the ethnic group in question—so far most data are from Caucasians. In some embodiments, more detailed human genome variation maps for specific populations are used. In some embodiments, Linkage disequilibrium based annotation is used.
3) Obtain information about the enrichment factor (prior) from the literature or public databases, such as location of the SNP within a region of the genome. Several enrichment factors, such as, for example, regulatory regions of a gene, exons (coding region of the gene), microRNA binding sites and evolutionary measures, are used, although others may be utilized. Some of these are general for most phenotypes, while some vary between phenotypes. Another enrichment factor is associated or co-morbid phenotypes. For example, it was shown how SNPs associated with bipolar disorder greatly increase the signal in schizophrenia.
4) The statistical package includes tools according to the utility. In some embodiments, model-free methods or model-based analysis is used. The model-based tool is useful for quantification. In short, Q-Q plots were used to visualize enrichment, and to aid in obtaining TDR values for the SNPs and increase replication rate. One can then calculate a FDR value for each SNP, after using a look-up table. The FDR value for each SNP is the output of the package, and a much improved tool for gene discovery is provided (very strong improvement in schizophrenia, 4-5 times more genes), discovery of overlapping genes (pleiotropy, e.g., between CVD risk and schizophrenia) etc.
5) In some embodiments, the model-based tools are used for improving technical calculations of the GWAS, such as correcting for inflation (Genomic Control), for calculating power, and for quantification of overlap between phenotypes (and identification of the SNPs involved in the overlap), and for estimating the polygenicity of a trait (how many genes have an effect, 1000-10000).
6) In some embodiments, a regression tool it used to combine all the enrichment factors including pleiotropic enrichment. This tool produces a FDR value for each SNP for the phenotype in question. In some embodiments, this forms the basis of the tool used for generalization performance (e.g., prediction of individuals based on their GWAS or deep sequencing profile). It was shown that the generalization performance increase 3-4 times compared to standard tools (See e.g.,
7) In some embodiments, systems and methods include updates on gene function (e.g., enrichment factors, system for continuous updates when new information becomes available), and all available GWAS studies (e.g., human traits of disorders, anonymous summary statistics, new GWAS as they become available), and a script for each utility. For example, some exemplary applications include: i) providing FDR values to new GWAS to improve discovery, and all the technical information needed (e.g., GC correction, power, etc) and providing pleiotropy information with all available phenotypes; ii) taking two new GWAS from two phenotypes and providing information about pleiotropy measures between the new phenotypes in addition; iii) taking deep sequencing data and providing information; and iv) providing an estimate of risk for specific phenotypes using a GWAS from an individual person.
The present invention also provides a variety of computer-related embodiments. Specifically, in some embodiments the invention provides computer programming for analyzing and comparing polymorphism to identify and characterize conditions.
The methods and systems described herein can be implemented in numerous ways. In one embodiment, the methods involve use of a communications infrastructure, for example the internet. Several embodiments of the invention are discussed below. It is also to be understood that the present invention may be implemented in various forms of hardware, software, firmware, processors, distributed servers (e.g., as used in cloud computing) or a combination thereof. The methods and systems described herein can be implemented as a combination of hardware and software. The software can be implemented as an application program tangibly embodied on a program storage device, or different portions of the software implemented in the user's computing environment (e.g., as an applet) and on the reviewer's computing environment, where the reviewer may be located at a remote site (e.g., at a service provider's facility).
For example, during or after data input by the user, portions of the data processing can be performed in the user-side computing environment. For example, the user-side computing environment can be programmed to provide for defined test codes to denote platform, carrier/diagnostic test, or both; processing of data using defined flags, and/or generation of flag configurations, where the responses are transmitted as processed or partially processed responses to the reviewer's computing environment in the form of test code and flag configurations for subsequent execution of one or more algorithms to provide a results and/or generate a report in the reviewer's computing environment.
The application program for executing the algorithms described herein may be uploaded to, and executed by, a machine comprising any suitable architecture. In general, the machine involves a computer platform having hardware such as one or more central processing units (CPU), a random access memory (RAM), and input/output (I/O) interface(s). The computer platform also includes an operating system and microinstruction code. The various processes and functions described herein may either be part of the microinstruction code or part of the application program (or a combination thereof) which is executed via the operating system. In addition, various other peripheral devices may be connected to the computer platform such as an additional data storage device and a printing device.
As a computer system, the system generally includes a processor unit. The processor unit operates to receive information, which generally includes test data (e.g., specific gene products assayed), and test result data, (e.g., the pattern of gastrointestinal neoplasm-specific marker detection results from a sample). This information received can be stored at least temporarily in a database, and data analyzed in comparison to a library of marker patterns known to be indicative of the presence or absence of a condition.
Part or all of the input and output data can also be sent electronically; certain output data (e.g., reports) can be sent electronically or telephonically (e.g., by facsimile, e.g., using devices such as fax back). Exemplary output receiving devices can include a display element, a printer, a facsimile device and the like. Electronic forms of transmission and/or display can include email, interactive television, and the like. In some embodiments, all or a portion of the input data and/or all or a portion of the output data (e.g., diagnosis or characterization of a condition) are maintained on a server for access, e.g., confidential access. The results may be accessed or sent to professionals as desired.
A system for use in the methods described herein generally includes at least one computer processor (e.g., where the method is carried out in its entirety at a single site) or at least two networked computer processors (e.g., where detected marker data for a sample obtained from a subject is to be input by a user (e.g., a technician or someone performing the assays)) and transmitted to a remote site to a second computer processor for analysis detection results is compared to a library of patterns known to be indicative of the presence or absence of a disease or condition, where the first and second computer processors are connected by a network, e.g., via an intranet or internet). The system can also include a user component(s) for input; and a reviewer component(s) for review of data, and generation of reports. Additional components of the system can include a server component(s); and a database(s) for storing data (e.g., as in a database or report), or a relational database (RDB) which can include data input by the user and data output. The computer processors can be processors that are typically found in personal desktop computers (e.g., IBM, Dell, Macintosh), portable computers, mainframes, minicomputers, tablet computer, smart phone, or other computing devices.
The input components can be complete, stand-alone personal computers offering a full range of power and features to ran applications. The user component usually operates under any desired operating system and includes a communication element (e.g., a modem or other hardware for connecting to a network using a cellular phone network, Wi-Fi, Bluetooth, Ethernet, etc.), one or more input devices (e.g., a keyboard, mouse, keypad, or other device used to transfer information or commands), a storage element (e.g., a hard drive or other computer-readable, computer-writable storage medium), and a display element (e.g., a monitor, television, LCD, LED, or other display device that conveys information to the user). The user enters input commands into the computer processor through an input device. Generally, the user interface is a graphical user interface (GUI) written for web browser applications.
The server component(s) can be a personal computer, a minicomputer, or a mainframe, or distributed across multiple servers (e.g., as in cloud computing applications) and offers data management, information sharing between clients, network administration and security. The application and any databases used can be on the same or different servers. Other computing arrangements for the user and server(s), including processing on a single machine such as a mainframe, a collection of machines, or other suitable configuration are contemplated. In general, the user and server machines work together to accomplish the processing of the present invention.
Where used, the database(s) is usually connected to the database server component and can be any device which will hold data. For example, the database can be any magnetic or optical storing device for a computer (e.g., CDROM, internal hard drive, tape drive). The database can be located remote to the server component (with access via a network, modem, etc.) or locally to the server component.
Where used in the system and methods, the database can be a relational database that is organized and accessed according to relationships between data items. The relational database is generally composed of a plurality of tables (entities). The rows of a table represent records (collections of information about separate items) and the columns represent fields (particular attributes of a record). In its simplest conception, the relational database is a collection of data entries that “relate” to each other through at least one common field.
Additional workstations equipped with computers and printers may be used at point of service to enter data and, in some embodiments, generate appropriate reports, if desired. The computers) can have a shortcut (e.g., on the desktop) to launch the application to facilitate initiation of data entry, transmission, analysis, report receipt, etc. as desired.
Embodiments of the present invention provide diagnostic, prognostic, and screening compositions, kits, and methods. In some embodiments, compositions, kits, and methods characterize and diagnose diseases and traits using one or more polymorphisms identified using the systems and methods described herein.
Embodiments of the present invention provide compositions and methods for detecting polymorphisms in one or more genes (e.g., to identity or diagnose diseases and traits). The present invention is not limited to particular variants. Exemplary variants for several traits are described in Examples 1-3, although the systems and methods described herein find use in the identification of polymorphisms in additional diseases and traits.
In some embodiments, 1 or more (e.g., 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40, 50, 100, 1000, 5000, or more) gene variants associated with a given disease or trait are utilized to diagnose or characterize a condition. The specific number of necessary, useful, or sufficient to diagnose or characterize a given trait can vary based on posterior effect sizes of the gene variants or the pleiotropy of the condition being diagnosed and characterized. The system and methods described herein find use in identifying the number of polymorphisms necessary, useful, or sufficient for diagnosing or characterizing a given condition.
In some embodiments, the systems and method described herein identify particular combinations of markers that show optimal function with different ethnic groups or sex, different geographic distributions, different stages of disease, different degrees of specificity or different degrees of sensitivity. Particular combinations may also be developed which are particularly sensitive to the effect of therapeutic regimens on disease progression (e.g., to customize treatment). Subjects may be monitored after a therapy and/or course of action to determine the effectiveness of that specific therapy and/or course of action.
In some embodiments, the present, invention provides information that indicates if a particular individual is predisposed to a particular disease or trait. In some embodiments, the present invention provides information useful in determining a treatment course of action (e.g., determining a particular drug or treatment regimen that is customized to the individual).
In some embodiments, the systems and methods described herein find use in research applications (e.g., in the analysis of polymorphism information to identify markers or identify pleiotropy information).
In some embodiments, the present invention provides systems and method for computation of polygenic personalized risk scores leveraging linkage disequilibrium (LD) genie annotation scores employing the statistical methodology described herein. In some embodiments, gene variant (e.g., single nucleotide polymorphisms (SNP)) posterior effect sizes are computed by repeatedly and randomly dividing subjects from a given study or collection of studies into disjoint training and replication subsamples and computing sample mean replication effect sizes conditional on training effect sizes. In some embodiments, computation of polygenic risk scores leverages pleiotropic effects with other traits. In some embodiments, computation of polygenic risk scores leverages LD genie annotation scores and pleiotropy simultaneously. In some embodiments, computation of polygenic risk scores leverages other types of prior information.
In some embodiments, genetic personalized risk scores summarize patient-level genomic variation as a single score per subject, summed over assayed gene variants. The polygenic risk score is computed as a linear or nonlinear function of the estimated statistical parameters, including per SNP allele effect size mean and/or estimates of variability. In some embodiments, linear weighting of each gene variant by its estimated posterior effect size optionally divided by its estimated posterior variance, given the observed association statistics with a given complex phenotype or disease diagnosis is utilized. In some embodiments, statistical methods are utilized to obtain maximal correlation of genetic risk scores with phenotypes in de novo subject samples, by obtaining posterior effect size estimates for each gene variant modulated by genie annotations and/or strength of association with pleiotropic phenotypes. In some embodiments, posterior effect sizes for each gene variant are multiplied by the corresponding gene variant values for a de novo subject and added together to calculate an overall risk score for a given trait or illness. In other embodiments, the posterior effect size for each gene variant are scaled by dividing by a measure of its variability before computing the polygenic risk score. In some embodiments, gene variant effect sizes below a given threshold are deleted before computing polygenic risk scores.
In some embodiments, polygenic risk scores also include other biomarkers of complex phenotypes or disease diagnosis. Other biomarkers of risk include, but are not limited to, age, gender, family history of illness, brain imaging phenotypes, etc.
In some embodiments, the statistical methodology leverages LD-weighted annotation scores and pleiotropic associations to compute polygenic normative variation scores, accounting for non-risk related genetic variation in complex phenotypes. Non-risk related variation in genotypes is genotypic variation correlated with (and hence predictive of) normal phenotypic variation in a complex phenotype. Variation in non-risk related genotypic variation is used to compute a single personalized non-risk genetic score per subject, summed over assayed non-risk gene variants. Each gene variant is weighted by its estimated posterior effect size and divided by its estimated posterior variance, given the observed association statistics with a given complex phenotype. In some embodiments, non-risk related genetic scores are used to determine phenotypic and/or developmental norms for subjects with specific genetic backgrounds.
In some embodiments, the statistical methodology is used to assist in the development of specialized genotyping chips that enable computation of genetic personalized risk scores and polygenic normative variation scores with maximal power to predict normative and non-normative variation in complex phenotypes and diseases in de novo samples. For example, in some embodiments, arrays that focus on a specific disease or population group are developed.
In some embodiments, the statistical methodology is used to predict complex phenotypes and disease diagnosis of offspring of two parents, given the parents' genotypes. In some embodiments, this is accomplished by randomly simulating multiple offspring and estimating polygenic risk scores for each simulated offspring. The distribution of polygenic risk scores across offspring is used to determine a distribution of polygenetic risk for a given complex phenotype or disease.
The following examples are provided in order to demonstrate and further illustrate certain preferred embodiments and aspects of the present invention and are not to be construed as limiting the scope thereof.
Ethics Statement
The relevant institutional review boards or ethics committees approved the research protocol of the individual GWAS used in the current analysis and all human participants gave written informed consent.
Participant Samples
GWAS results were obtained in the form of summary statistics p-values from the Psychiatric GWAS Consortium (PGC)—Schizophrenia and Bipolar Disorder Working Groups. The schizophrenia (SCZ) GWAS summary statistics results were obtained from the PGC Schizophrenia Work Group[12], which consisted of 9,394 cases with schizophrenia or schizoaffective disorder and 12,462 controls (52% screened) from a total of 17 samples from 11 countries. Semi-structured interviews were used by trained interviewers to collect clinical information, and operational criteria were used to establish diagnosis. The quality of phenotypic data was verified by a systematic review of data collection methods and procedures at each site, and only studies that fulfilled these criteria were included. Controls were selected from the same geographical and ethnic populations as cases. For further details on sample characteristics and quality control procedures applied, please see Ripke et al[12].
The bipolar disorder (BD) GWAS summary statistics results were obtained from, the PGC Bipolar Disorder Working Group[13], which consisted of n=16,731 including 7481 cases and 9250 controls, from 11 studies from 7 countries. Standardized semi-structured interviews were used by trained interviewers to collect clinical information about lifetime history of psychiatric illness and operational criteria applied to make lifetime diagnosis according to recognized classifications. All cases have experienced pathologically relevant episodes of elevated mood (mania or hypomania) and meet operational criteria for a BD diagnosis. The sample consisted of BD I (84%), BD II (11%), schizoaffective disorder bipolar type (4%), and BD NOS (1%). Controls were selected from the same geographical and ethnic populations as cases. For further details on sample characteristics and quality control procedures applied, please see Sklar et al[13].
Due to overlapping control samples in these studies, the common controls were split randomly, and divided between the two case-control analyses. All results presented here are based on these nonoverlapping control samples, with n=9379 cases and n=7736 samples in schizophrenia, and n=6990 cases and n=4820 controls in bipolar disorder analyses.
Statistical Analyses
Analyses implemented here were motivated by previously published stratified FDR methods[5,33]. However, it was found that stratified empirical cdfs exhibited a high degree of variability. Instead, empirical cdfs were obtained for the first phenotype conditional on nominal p-values of the second being at or below a given threshold. These conditional empirical cdfs vary more smoothly as a function of pvalue thresholds in the second (associated) phenotype than do empirical cdfs employing disjoint strata. Conditional FDR estimates derived from the conditional empirical cdfs are a simple extension of Efron's Empirical Bayes FDR methods[40].
One advantage of the model-free empirical cdf approach is the avoidance of bias in conditional FDR estimates from model misspecification. However, there are inherent, limitations to model-free approaches, especially with respect to inferring properties of the non-null distribution and, consequently, estimating power to detect non-null effects. Complementary model-based analyses are provided that estimate conditional and conjunctional local false discovery rate (fdr)[27].
Stratified Q-Q Plots
Q-Q plots compare a nominal probability distribution against an empirical distribution. In the presence of all null relationships, nominal p-values form a straight line on a Q-Q) plot when plotted against the empirical distribution. For each phenotype, for all SNPs and for each categorical subset (strata), −log10 nominal p-values were plotted against −log10 empirical p-values (stratified Q-Q plots). Leftward deflections of the observed distribution from the projected null line reflect increased tail probabilities in the distribution of test statistics (z-scores) and consequently an over-abundance of low p-values compared to that expected by chance, also termed “enrichment”.
Genomic Control
The empirical null distribution in GWAS is affected by global variance inflation due to population stratification and cryptic relatedness[39] and deflation due to over-correction of test statistics for polygenic traits by standard genomic control methods[40]. A control method leveraging only intergenic SNPs which are likely depleted for true associations (Schork et al., under review) was applied. First, the SNPs was annotated to genie (5″UTR, exon, intron, 3″UTR) and intergenic regions using information from the 1000 Genomes Project (1KGP). As illustrated in
Q-Q Plots for Pleiotropic Enrichment
To assess pleiotropic enrichment, a Q-Q plot conditioned by “pleiotropic” effects was used. For a given associated phenotype, enrichment for pleiotropic signals is present if the degree of deflection from the expected null line is dependent on SNP associations with the second phenotype. Conditional Q-Q plots were constructed of empirical quantiles of nominal −log 10(p) values for SNP association with schizophrenia for all SNPs, and for subsets (strata) of SNPs determined by the nominal p-values of their association with bipolar disorder. Specifically, the empirical cumulative distribution of nominal p-values for a given phenotype for all SNPs and for SNPs with significance levels below the indicated cut-offs for the other phenotype (−log10(p)≧0, −log10(p)≧1, −log10(p)≧2, −log10(p)≧3 corresponding to p<1, p<0.1, p<0.01, p<0.001, respectively) was computed. The nominal p-values (−log10(p)) are plotted on the y-axis, and the empirical quantiles (−log10(q), where q=1−cdf(p)) are plotted on the x-axis. To assess for polygenic effects below the standard GWAS significance threshold, the conditional Q-Q plots were focused on SNPs with nominal −log 10(p)<7.3 corresponding to p>5×10−8).
Conditional FDR
Enrichment seen in conditional Q-Q plots can be directly interpreted in terms of FDR [29]), The stratified FDR method[26], previously used for enrichment of GWAS based on linkage information[5] was applied. Specifically, for a given p-value cutoff, the FDR is defined as
FDR(p)=π0F0(p)/F(p), [1]
where π0 is the proportion of null SNPs, F0 is the null cumulative distribution function (cdf), and F is the cdf of all SNPs, both null and non-null; see below for details on this simple mixture model formulation[41]. Under the null hypothesis, F0 is the cdf of the uniform distribution on the unit interval [0,1], so that Eq. [1] reduces to
FDR(p)=π0p/F(p), [2]
The cdf F can be estimated by the empirical cdf q=Np/N, where Np is the number of SNPs with pvalues
less than or equal to p, and N is the total number of SNPs. Replacing F by q in Eq. [2], one gets
FDR(p)≈p/q, [3]
which is biased upwards as an estimate of the FDR[41]. Replacing π0 in Equation [3] with unity gives an estimated FDR that is further biased upward. If π0 is close to one, as is likely true for most GWAS, the increase in bias from Eq. [3] is minimal. The quantity 1−p/q, is therefore biased downward, and hence is a conservative estimate of the TDR. Note, Eq. [3] is the Empirical Bayes estimate of the Bayesian FDR described by Efron[40]. Referring to the formulation of the Q-Q plots, that Eq. [3] is equivalent to the nominal p-value divided by the empirical quantile, as defined earlier. Given the −log 10 of the Q-Q plots one obtains:
−log 10(FDR(p))≈log10(q)−log10(p) [4]
demonstrating that the (conservatively) estimated FDR is directly related to the horizontal shift of the curves in the conditional Q-Q plots from the expected line x=y, with a larger shift corresponding to a smaller FDR. This is illustrated in
Conditional Statistics—Probability of Association with One Disorder
The conditional FDR is defined as the posterior probability that a given SNP is null for the first phenotype given that the p-values for both phenotypes are as small or smaller as the observed p-values. Formally, this is given by
FDR(p1|p2)=π0(p2)p1/F(p1|p2), [5]
where p1 is the p-value for the first phenotype, p2 is the p-value for the second, and F(p1|p2) is the conditional cdf and π0(p2) the conditional proportion of null SNPs for the first phenotype given that pvalues for the second phenotype are p2 or smaller. Eq. [5] makes the assumption, reasonable for independent GWAS, that summary statistics are independent across phenotypes if they are null for at least one phenotype. A conservative estimate of FDR(p1|p2) is produced by setting π0(p2)=1 and using the empirical conditional cdf in place of F(p1|p2) in Eq. [5]. This is a straightforward generalization of the Empirical Bayes approach developed by Efron[40]. A conditional FDR value for schizophrenia given bipolar disorder p-values (denoted by FDR SCZ BD) is assigned to each SNP by computing conditional FDR estimates on a grid and interpolating these estimates into a twodimensional look-up table (
Conjunction Statistics—Test of Association with Both Phenotypes
In order to identify which of the SNPs associated with schizophrenia and bipolar disorder, a conjunction testing procedure as outlined for p-value statistics in Nichols et al.[42], adopted to FDR statistics based on the stratified FDR approach[5,26], was used. Conjunction FDR is defined as the posterior probability that a given SNP is null for both phenotypes simultaneously when the p-values for both phenotypes are as small or smaller than the observed p-values. Formally, conjunction FDR is given by
FDR(p1,p2)=π0(p1,p2)F0(p1,p2)/F(p1,p2), [6]
where π0(p1, p2) is the proportion of SNPs null for both phenotypes simultaneously, F0(p1, p2)=p1 p2 is the joint null cdf, and F(p1, p2) is the joint overall cdf.
Conditional empirical cdfs provide a model-free method to obtain conservative estimates of Eq. [6]. This can be seen as follows. Estimate the conjunction FDR by
FDRSCZ&BD=max{FDRSCZ|BDFDRBD|SCZ} [7]
where FDR SCZ|BD and FDR BD|SCZ (the estimated conditional FDRs described above) are conservative (upwardly biased) estimates of Eq. [5]. Thus, Eq. [7] is a conservative estimate of max {p1/F(p1|p2), p2/F(p2|p1)}=max{p1 F2(p2)/F(p1, p2), p2 F1(p1)/F(p1, p2)}. For enriched samples, pvalues will tend to be smaller than predicted from the uniform distribution, so that F1(p1)≧p1 and F2(p2)≧p2. Hence, max{p1 F2(p2)/F(p1, p2), p2 F2(p1)/F(p1, p2)}≧max{p1 p2/F(p1, p2), p2 p1/F(p1, p2)}=p1 p2/F(p1, p2)≧π0(p1, p2) p1 p2/F(p1, p2). The last quantity is precisely the conjunction FDR defined by Eq. [6]. Thus, Eq. [7] is a conservative model-free estimate of the conjunction FDR.
The conjunction FDR values were assigned by interpolation into a bi-directional two-dimensional look-up table (
Conditional Manhattan Plots
To illustrate the localization of the genetic markers associated with schizophrenia given bipolar disorder effect, and vice versa, a “Conditional Manhattan plot”, plotting all SNPs within an LD block in relation to their chromosomal location was used. As illustrated in
Conjunction Manhattan Plots
To illustrate the localization of the pleiotropic genetic markers association with both schizophrenia and bipolar disorder, a “Conjunction Manhattan plot”, plotting all SNPs with a significant conjunction FDR within an LD block in relation to their chromosomal location was used. As illustrated in
Four-Groups Mixture Model
Here, a model-based methodology for computing pleitropy-informed conditional and conjunction analyses, complementary to the model-free approach presented in the main text is described. Let z be the GWAS test statistic (z-score) with corresponding nominal significance p (two-tailed probability of observed z-score under the null hypothesis of no effect). A standard Bayesian two-groups mixture model [Efron B (2010) Large-scale inference: empirical Bayes methods for estimation, testing, and prediction. Cambridge; New York: Cambridge University Press. xii, 263] is given by
f(z)=π0f0(z)+(1−π0)f1(z) [S1]
where f0 is the null distribution (e.g., standard normal after appropriate genomic control), f1 is the non-null distribution (which may be estimated parametrically or non-parametrically, and π0 is the proportion of null SNPs. From model [S1] the Bayesian False Discovery Rate (denoted as FDR) and the local False Discovery Rate (denoted as fdr) for a given effect size z are
FDR(z)=π0F0(z)/F(z) [S2]
fdr(z)=π0f0(z)/f(z) [S3]
where F0(z) and F(z) are the cumulative distribution functions (cdfs) corresponding to f0(z) and f(z), respectively. Following is an extension to conditional and conjunctional fdr (Eq. [S3]); it is straightforward to extend this to include conditional and conjunction FDR (Eq. [S2]). Eq. [S1] is generalized to bivariate z-scores from two phenotypes (z1 for phenotype 1 and z2 for phenotype 2) using a bivariate density from a four-groups mixture model
f(z1,z2)=π0f0(z1,z2)+π1f1(z1,z2)+π2f2(z1,z2)+π3f3(z1,z2) [S4]
where π0 is the proportion of SNPs for which both phenotypes are null, π1 is the proportion of SNPs where phenotype 1 is non-null and phenotype 2 is null, π2 is the proportion of SNPs where phenotype 1 is null and phenotype 2 is non-null, and 3 is the proportion of SNPs where both phenotypes are non-null (i.e., the pleiotropic SNPs). The mixture densities in [S4] are given by
f
0(z1,z2)=φ(z1)φ(z2)
f
1(z1,z2)=g1(z1)φ(z2)
f
2(z1,z2)=φ(z1)g2(z2)
f
3(z1,z2)=g1(z1)g2(z2) [S5]
where φ( ) denotes the theoretical null density and g1 and g2 denote the non-null marginal densities of z1 and z2, respectively. Modeling the φ with the standard normal and g1 and g2 with Normal-Laplace densities fits the empirical z-scores well. Another parametric model providing a very good fit to the squared z-scores (z2) sets φ to a central chi-squared density with one degree of freedom (χ21) and g1 and g2 to Weibull densities with scale parameters α1 and α2 and shape parameters β1 and β2 for g1 and g2, respectively. More generally f3 is modeled with marginal densities as above but allowing for dependence between pleiotropic (jointly non-null) SNPS using, for example, a copula formulation [Joe H (1997) Multivariate models and multivariate dependence concepts: Chapman & Hall/CRC]. The proportions π=(π0,π1,π2,π3) and the parameters of the non-null distributions can be estimated using Bayesian methods such as Markov Chain Monte Carlo (MCMC) algorithms or maximum likelihood (ML) estimation.
The estimated vector of probabilities π=(π0,π1,π2,π3) from these fits can be used to test whether the degree of pleiotropy is significantly higher than expected by chance if both phenotypes were independent. Independence implies that the joint pdf of both phenotype summary scores is a product of two two-group mixture models (two independent versions of Eq. [S1]). It is easy to show that testing for excess pleiotropy over that predicted by independence is equivalent to showing that π3>π1π2/π0 in Eq. [S4] or equivalently that the log-odds ratio
LOR(Phen. 1. Phen. 2)=log {π3/1−π3}−log {(π1π2/π0)/(1−π1π2/π0)} [S6]
is greater than zero. Using a multivariate normal approximation to the ML estimates with covariance obtained from the inverse Fisher information matrix, estimates of LOR with 95% confidence intervals are: LOR(SCZ,BD)=10.3 [4.1, 16.4], LOR(SCZ,T2D)=1.3 [0.2, 2.5], and LOR(BD,T2D)=1.5 [0.6, 2.4]. In particular, the departure from independence of SCZ and BD is highly significant, with a 95% CI bounded well above zero. ML estimates and 95% CIs were produced using the SCZ/BD data z2 values estimated using non-overlapping controls, and include an adjustment to account for correlation of SNPs (e.g., LD) that assumes an effective degree of freedom of 500,000 independent SNPs.
The proportion of pleiotropic SNPs is estimated for each phenotype. For example, π3/(π1+π3) is the proportion of pleiotropic SNPs for phenotype 1 (e.g., the proportion of non-null SNPs for phenotype 1 that are also non-null for phenotype 2). Again using the ML estimates from the χ21-Weibull model, the proportion of pleiotropic SNPs for BD with SCZ was 0.56 (95% CI: [0.48, 0.64]), the proportion for SCZ with BD was 0.94 [0.37, 1.00], the proportion for SCZ with T2D was 0.04 [0.01, 0.10], the proportion for BD with T2D was 0.05 [0.02, 0.09]. ML estimates and 95% CIs were again produced using the SCZ/BD data z-score estimates with non-overlapping controls, and include an adjustment to account for correlation of SNPs. The huge increase in power for BD|SCZ noted below is due to high proportion of non-null SCZ SNPs that are also non-null BD SNPs. As a point of comparison, two split-half samples are produced using the SCZ data, showing a pleiotropic overlap of 0.992 [0.988, 0.996] of SCZ with itself.
Conditional and Conjunction Local False Discovery Rate
From the ML-estimates of the four-groups mixture pdf (Eq. [S4]) one can compute ML estimates of the conditional pdf of z1 given z2 and hence the conditional fdr of the first phenotype given the second
fdr(z1|z2)=f(z1|z1null,z2)Pr(z1null|z2)/f(z1|z2) [S6]
where f(z1|z1 null, z2) is the null density of z1 conditional on z2, Pr(z1 null|z2) is the probability that z1 is null given z2, and f(z1|z2) is the mixture density of z1 conditional on z2. With component densities as given in Eq. [S5], this becomes
fdr(z1|z2)=φ(z1)[π0φ(z2)+π2(z2)]/f(z1,z2), [S7]
where f(z1, z2) is the joint density given in Eq. [S4]. Look-up tables were produced using Eq. [S7], with ML estimates of unknown parameters, again assuming the χ21-Weibiull model for z2.
The conjunctional fdr of both phenotypes is computed as
fdr(z1,z2)=f(z1,z2|z1null,z2null)Pr(z1null,z2null)/f(z1,z2) [S8]
where f(z1, z2|z1 null, z2 null)=φ(z1) φ(z2) is the joint null density of z1 and z2, Pr(z1 null, z2 null) is the probability that both z1 and z2 are null, and f(z1, z2) is the joint pdf of z1 and z2. With densities given in Eq. [S5], this becomes
fdr(z1,z2)=π0φ(z1)φ(z2)/f(z1,z2) [9]
A joint fdr look-up table for SCZ & BD is presented in
Conditional Local False Discovery Rate and Power
Conditional local false discovery rates fdr(z1|z2) can lead to significant increases in power when two phenotypes are genuinely pleiotropic (i.e., when LOR(Phen. 1, Phen. 2) is significantly larger than zero). Here, power is defined in terms of the probability of rejecting the null hypothesis for SNPs that are in fact non-null for a given fdr threshold α. In this sense power corresponds to sensitivity to detect non-null SNPs and power diagnostics correspond can be presented as ROC-type curves as detailed in Efron [Efron B (2007) Size, power and false discovery rates. The Annals of Statistics 35: 1351-1377]. In
Note, estimates of power in the sense described above are sensitive to assumptions about the shape of the non-null distribution near zero. However, relative power (the ratio of sensitivity of conditional fdr with marginal fdr for a given threshold α) is well identified. For example, using the fdr cut-off α≦0.05, the ratio of power for conditional fdr of BD|SCZ vs. marginal fdr of BD is 44.4. The ratio of power for unconditional vs. conditional fdr for SCZ|BD is 2.4, indicating improvement of power but to a much lesser degree. In contrast, the ratio of power for unconditional vs. conditional fdr for SCZ|T2D is 1.00, indicating no improvement whatsoever.
Results
Q-Q plots of schizophrenia SNPs stratified by association with bipolar disorder and vice versa Under large-scale testing paradigms, such as GWAS, quantitative estimates of likely true associations can be estimated from the distributions of summary statistics[27,28]. A common method for visualizing the “enrichment” of statistical association relative to that expected under the global null hypothesis is through Q-Q plots of nominal p-values obtained from GWAS summary statistics. The usual Q-Q curve has as the y-ordinate the nominal p-value, denoted by “p”, and the x-ordinate the corresponding value of the empirical cdf, denoted by “q”. Under the global null hypothesis the theoretical distribution is uniform on the interval [0,1]. As is common in GWAS, one instead plots −log10 p against −log10 q to emphasize tail probabilities of the theoretical and empirical distributions. Thus, enrichment results in a leftward shift in the Q-Q curve, corresponding to a larger fraction of SNPs with nominal −log10 p-value greater than or equal to a given threshold. Conditional Q-Q plots are formed by creating subsets of SNPs based on levels of an auxiliary measure for each SNP, and computing Q-Q plots separately for each level. If SNP enrichment is captured by variation in the auxiliary measure, this is expressed as successive leftward deflections in a conditional Q-Q plot as levels of the auxiliary measure increase.
Conditional Q-Q plots for schizophrenia conditioned on nominal p-values of association with bipolar disorder (SCZ|BD;
Conditional True Discovery Rate (TDR) in schizophrenia is increased by bipolar disorder, and vice versa.
Since categories of SNPs with stronger pleiotropic enrichment are more likely to be associated with schizophrenia, to maximize power for discovery all tag SNPs should not be treated interchangeably. Specifically, variation in enrichment across pleiotropic categories is expected to be associated with corresponding variation in the TDR (equivalent to 1-FDR)[29] for association of SNPs with schizophrenia. A conservative estimate of the TDR for each nominal p-value is equivalent to 1−(p/q), easily read from the stratified Q-Q plots (see Material and Methods). This relationship is shown for schizophrenia conditioned on nominal bipolar disorder p-values (SCZ|BD;
Schizophrenia Gene Loci Identified with Conditional FDR
A “conditional” Manhattan plot for schizophrenia showing the FDR conditional on bipolar disorder (
Bipolar Disorder Gene Loci Identified with Conditional FDR
A “conditional” Manhattan plot for bipolar disorder showing the FDR conditional on schizophrenia (
Pleiotropic Gene Loci in Both Schizophrenia and Bipolar Disorder Identified with Conjunctional FDR
To identify pleiotropic loci in schizophrenia and bipolar disorder, a conjunctional FDR analysis was performed and used to construct a “conjunction” Manhattan plot (
The model-based analysis using a bivariate mixture model showed that a very high proportion of the non-null schizophrenia SNPs are also non-null for bipolar disorder, leading to large increases in power (
0.014
0.022
0.025
0.037
0.018
0.049
0.035
0.00115
0.01859
kinase 4) regulator
-specific transcript 1 (non-protein coding)
9
homolog 9
(Borderline)
KA2
deficient-like 1
significant p-values.
were used as indentity previous results.
indicates data missing or illegible when filed
Genome-Wide Association Study (GWAS) Data
Fourteen phenotypes, body mass index (BMI) [30], height, waist to hip ratio [31](WHR), Crohn's disease [32](CD), ulcerative colitis [33](UC), schizophrenia [34](SCZ), bipolar disorder [35](BD), smoking behavior as measured by cigarettes per day [36](CPD), systolic and diastolic blood pressure [37](SBP, DBP), and plasma lipids [38](triglycerides, TG, total cholesterol, TC, high density lipoprotein, HDL, low density lipoprotein, LDL), were considered. Genome-wide association study (GWAS) results were obtained as summary statistics (p-values or z-scores) from public access websites (BMI, Height, WHR, TC, TG, HDL, LDL: GIANT consortium data files; IBD Genetics; Psychiatric Genomics Consortium; Center for statistical genetics and the University of Michigan; Geneva University Hospital—Tulipe Center For Cardiovascular Research), published supplementary material (SBP, DBP; The International Consortium for Blood Pressure Genome—Wide Association Studies, Nature 478, 103-109 (6 Oct. 2011)), or through collaborations with investigators (CD, UC, SCZ, BD). For CD pre-meta-analysis, sub-study specific p-values and effect sizes (z-scores) were obtained from the study principal investigators. In total these studies considered more than 1.3 million phenotypic observations, but considerable sample overlap makes the number of unique individuals much less.
GWAS Summary Statistics Processing.
The summary statistics from the respective GWAS meta-analyses, derived according to best practices, were used as-is. No further processing was performed, with the exception of intergenic inflation control (described below). Results from SNPs with reference SNP (rs) numbers that did not map to the 1000 genomes project (1KGP) reference panel were excluded.
Positional Annotation Categories
Bi-allelic SNP genotypes from the European reference sample provided by the November 2010 release of Phase 1 of the 1KGP were obtained in pre-processed form. Using Plink version 1.07 [39,40] 1KGP SNPs with a minor allele frequency less than 1%, missing in more than 5% of individuals and/or violating Hardy-Weinberg equilibrium (p<1×10−6) were excluded from the reference panel. Individuals missing more than 10% of genotypes were excluded. Each remaining 1KGP SNP was assigned a single, mutually exclusive genic annotation category based on its genomic position (hg19). Genic annotation categories were: 1) 10,000 to 1,001 base pairs upstream (10 k Up); 2) 1,000 to 1 base pair upstream (1 k Up); 3) 5′ untranslated region (5′UTR); 4) exon; 5) intron; 6) 3′ untranslated region (3′UTR); 7) 1 to 1,000 base pairs downstream (1 k Down); 8) 1,001 to 10,000 base pairs downstream (10 k Down), all with reference to protein coding genes only. Annotations were assigned based on the first gene transcript listed in the UCSC known genes database [41]. In total 9,078,405 1KGP SNPs were assigned positional categories. All positional categories were scored 0 or 1.
Linkage Disequilibrium (LD) Weighted Scoring
For each GWAS tag SNP a pairwise correlation coefficient approximation to LD (r2) was calculated for all 1KGP SNPs within 1,000,000 base pairs (1 Mb) of the SNP using Plink version 1.07 [39,40]. LD scores were thresholded providing continuous valued estimates from 0.2 to 1.0; r2 values<0.2 were set to 0 and each SNP was assigned an r2 value of 1.0 with itself. LD-weighted annotation scores were computed as the sum of r2 LD between the tag SNP and all 1KGP SNPs positioned in a particular category. Each tag SNP was assigned to every LD-weighted annotation category for which its annotation score was greater than or equal to 1.0. The resulting LD-weighted annotation categories were not mutually exclusive such that each GWAS tag SNP could be annotated with multiple categories. All analyses were repeated using a second set of LD thresholding parameters and found to be robust.
Intergenic SNPs.
Intergenic SNPs were determined after LD-weighted scoring and defined as having LD-weighted annotations scores for each of the eight categories equal to zero. In addition they were defined to not be in LD with any SNPs in the 1KGP reference panel located within 100.000 base pairs of a protein coding gene, within a noncoding RNA, within a transcription factor binding site nor within a microRNA binding site. SNPs labeled intergenic were defined to be a specific collection of non-genic SNPs chosen to not represent any functional elements within the genome (including through LD). Because of how they are defined these SNPs are hypothesized to represent a collection of null associations. Other non-genic categories (1 k up, 10 k up, 1 k down and 10 k down) were included in the analyses to ensure SNPs not too far away from genes, but not within protein coding genes, were represented by non-genic categories and enrichment due to these SNPs was not solely attributed to LD with genie categories.
Stratified Q-Q Plots and Enrichment
Q-Q plots compare two probability distributions. For each phenotype, for all SNPs and for each categorical subset, −log10 nominal p-values were plotted against −log10 empirical p-values. Leftward deflections of the observed distribution from the projected null line reflect increased tail probabilities in the distribution of test statistics (z-scores) and consequently an over-abundance of low p-values compared to that expected by chance. This deflection is referred to as “enrichment (
The significance of the annotation enrichment was estimated using two sample Kolmogorov-Smirnov (KS) Tests to compare the distribution of test statistics in each genic annotation category to the distribution of the intergenic category, for each phenotype. SNPs were pruned randomly to approximate independence (r2<0.2) ten times.
Intergenic Inflation Control
The empirical null distribution in GWAS is affected by global variance inflation due to factors including population stratification and cryptic relatedness [17] and deflation due to over-correction of test statistics for polygenic traits by standard genomic control methods. A control method leveraging was applied only intergenic SNPs which are likely depleted for true associations. All p-values were converted into z-scores, and, for each phenotype, the genomic inflation factor [16], λGC, was estimated for intergenic SNPs. All test statistics were divided by λ GC.
The inflation factor, λGC was computed as the median z-score squared divided by the expected median of a chi-square distribution with one degree of freedom or all phenotypes except CPD, where the 0.95 quantile was used in place of the median. 4.
Quantification of Categorical Enrichment
For each phenotype, enrichment was measured as the mean(z-score2 −1) for each category and normalized by the largest value per phenotype. The mean(z-score2 −1) is a conservative estimate of the variance attributable to non-null SNPs, given a standard normal null distribution and a non-null distribution symmetric around zero.
Q-Q Plots and False Discovery Rate (FDR)
Enrichment seen in the conditional Q-Q plots can be directly interpreted in terms of the FDR. Specifically, for a given p-value cutoff, the Bayes FDR [17] is defined as
FDR(p)=π0F0(p)/F(p), [1]
where π0 is the proportion of null SNPs, F0 is the null cdf, and F is the cdf of all SNPs, both null and non-null. Under the null hypothesis, F0 is the cdf of the uniform distribution on the unit interval [0,1], so that Eq. [1] reduces to
FDR(p)=π0p/F(p). [2]
The cdf F can be estimated by the empirical cdf q=Np/N, where Np is the number of SNPs with p-values less than or equal to p, and N is the total number of SNPs. Replacing F by q and replacing π0 with unity in Eq. [2]
FDR(p)≈p/q, [3]
This is upwardly biased, and hence p/q is conservative estimate of the FDR, and 1−p/q is a conservative estimate of the Bayes TDR[17].
If π0 is close to one, as is likely true for most GWAS, the increase in bias from setting π0 to one in Eq. [3] is minimal. The quantity 1−p/q, is therefore biased downward, and hence a conservative estimate of the TDR.
Referring to the formulation of the Q-Q plots, FDR(p) is equivalent to the nominal p-value under the null hypothesis divided by the empirical quantile of the p-values. Given the −log10 transformation applied to the Q-Q plots,
−log10(FDR(p))≈log10(q)−log10(p) [4]
demonstrating that the (conservatively) estimated FDR is directly related to the horizontal shift of the curves in the stratified Q-Q plots from the expected line x=y, with a larger shift corresponding to a smaller FDR. For the TDR plots in
Eq. [3] is the Empirical Bayes point estimate of the Bayes FDR given in Efron (2010). Using Eq. [3] to control FDR (e.g., the expected proportion of falsely rejected null hypotheses) [21] is closely related to the “fixed rejection region” approach of Storey[47,48]. Specifically, Storey[47] showed, for a given FDR α, rejecting all null hypotheses such that p/q<α is equivalent to the Benjamini-Hochberg procedure and provides asymptotic control of the FDR to α if the true null p-values are independent and uniformly distributed. Storey[47] also noted that asymptotic control is preserved under positive blockwise dependence, whereas Schwartzman and Lin [49] showed that Eq. [3] is a consistent estimator of FDR for asymptotically sparse dependence (e.g., the proportion of correlated pairs of p-values goes to zero as the number of hypothesis tests becomes large). Sparse dependence is a good description of the dependence present in GWAS data; for example, based on a threshold of R2>0.05 within 1,000,000 basepairs, one can estimate the ratio of correlated pairs to total pairs of p-values at 0.000128.
Replication Rate
For each of eight sub-studies contributing to the final meta-analysis in the CD report z-scores were independently adjusted using intergenic inflation control. For each of 70 (8 choose 4) possible combinations of four-study discovery and four-study replication sets, the four-study combined discovery z-score and four-study combined replication z-score for each SNP were calculated as the average z-score across the four studies, multiplied by two (the square root of the number of studies). For discovery samples the z-scores were converted to two-tailed p-values, while replication samples were converted to one-tailed p-values preserving the direction of effect in the discovery sample. For each of the 70 discovery-replication pairs cumulative rates of replication were calculated over 1000 equally-spaced bins spanning the range of negative log10(p-values) observed in the discovery samples. The cumulative replication rate for any bin was calculated as the proportion of SNPs with a −log 10(discovery p-value) greater than the lower bound of the bin with a replication p-value<0.05. Cumulative replication rates were calculated independently for each of the eight genic annotation categories as well as intergenic SNPs and all SNPs. For each category, the cumulative replication rate for each bin was averaged across the 70 discovery-replication pairs and the results are reported in
Stratified False Discovery Rates:
A multiple linear regression was used to predict the tagged variance (z2) for each SNP in the height GWAS from the unthresholded LD-weighted annotation scores. Using the category weights determined from the variance regression on the height GWAS, the tagged variance for each SNP was predicted for each other phenotype. For each phenotype, SNPs were grouped into strata according to the rank of their predicted tagged variance. Enrichment for each stratum was demonstrated using QQ-plots as described above. Sun et al [9] described a stratified false discovery rate (sFDR) procedure which results in improved statistical power over traditional FDR methods [16] when a collection of statistical tests can be grouped into disjoint strata with different levels of enrichment. In order to demonstrate the utility of using genic annotation categories in combination with sFDR for increasing power, the number of SNPs deemed significant at a given FDR threshold using both traditional[21] and stratified FDR was computed, where the strata were determined by the predicted tagged variance for each SNP based on regression weights determined from the height GWAS summary statistics (
For all studies, Genome-wide association study (GWAS) results in the form of summary statistic p-values were obtained from public access websites (Speliotes E K, Willer C J, Berndt S I, Monda K L, Thorleifsson G, et al. (2010) Association analyses of 249,796 individuals reveal 18 new loci associated with body mass index. Nat Genet 42: 937-948; Lango Allen H, Estrada K, Lettre G, Berndt S I, Weedon M N, et al. (2010) Hundreds of variants clustered in genomic loci and biological pathways affect human height. Nature 467: 832-838; Heid I M, Jackson A U, Randall J C, Winkler T W, Qi L, et al. (2010) Meta-analysis identifies 13 new loci associated with waist-hip ratio and reveals sexual dimorphism in the genetic basis of fat distribution. Nat Genet 42: 949-960; Teslovich T M, Musunuru K, Smith A V, Edmondson A C, Stylianou I M, et al. (2010) Biological, clinical and population relevance of 95 loci for blood lipids. Nature 466: 707-713), (Ehret G B, Munroe P B, Rice K M, Bochud M, Johnson A D, et al. (2011) Genetic variants in novel pathways influence blood pressure and cardiovascular disease risk. Nature 478: 103-109) or through collaboration with investigators (Franke A, McGovern D P, Barrett J C, Wang K, Radford-Smith G L, et al. (2010) Genome-wide meta-analysis increases to 71 the number of confirmed Crohn's disease susceptibility loci. Nat Genet 42: 1118-1125; Anderson C A, Boucher G, Lees C W, Franke A, D'Amato M, et al. (2011) Meta-analysis identifies 29 additional ulcerative colitis risk loci, increasing the number of confirmed associations to 47. Nat Genet 43: 246-252; Consortium TSPG-WASG (2011) Genome-wide association study identifies five new schizophrenia loci. Nat Genet 43: 969-976; Group PGCBDW (2011) Large-scale genome-wide association analysis of bipolar disorder identifies a new susceptibility locus near ODZ4. Nat Genet 43: 977-983). For Crohn's disease (CD) (Franke et al., supra) pre-meta-analysis, sub-study specific p-values and effect sizes (z-scores) were obtained from the study principal investigators. See Table 11.
In total over 1.3 million phenotypic observations were considered; however, due to considerable overlap in samples, the number of unique individuals surveyed is significantly less. Blood pressure phenotypes (systolic blood pressure; SBP, diastolic blood pressure; DBP) were a part of one study sample (Ehret et al., supra) as were lipid traits (triglycerides; TG, total Cholesterol; TC, High density lipoprotein; HDL, Low density lipoprotein; LDL) (Teslovich et al., supra). In addition, Body Mass Index (BMI) (Speliotes et al., sura), Height (Lango et al., supra) and Waist-hip-ratio (WHR) (Heid et al., supra) all arose from the GIANT consortium and there is thus much sample redundancy.
The samples used in the lipids GWAS (Teslovich et al., supra) overlap considerably with the GIANT consortium samples, as do the samples used in the smoking GWAS (Consortium TaG (2010) Genome-wide meta-analyses identify multiple loci associated with smoking behavior. Nat Genet 42: 441-447). The Schizophrenia (Consortium, supra) and Bipolar Disorder GWAS (Group, supra) share some controls. The phenotypes, however, are diverse.
Bi-allelic SNP genotypes from the European reference sample provided by the November 2010 release of Phase 1 of the 1000 Genomes Project (1KGP) were obtained in pre-processed form. Additional quality control was performed on the 1KGP data using Plink version 1.07 (Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira M A, et al. (2007) PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81: 559-575). 1KGP genotypes were pruned according to standard GWAS procedures, removing all SNPs with a minor allele frequency less than 1%, missing in more than 5% of individuals or violating Hardy-Weinberg equilibrium (p<1×10−6). Individuals missing more than 10% of genotypes were excluded. Plink implementations of identity by state (IBS) and identity by descent (IBD) analysis were used to remove one individual from each related pair present and implementations of multidimensional scaling were used to ensure population homogeneity within the reference sample.
Each SNP in the 1KGP based reference sample was assigned a mutually exclusive category based on its position within the genome. A computational annotation pipeline (Torkamani A, Scott-Van Zeeland A A, Topol E J, Schork N J (2011) Annotating individual human genomes. Genomics 98: 233-241), which calls upon a variety of publicly available tools and databases to aggregate comprehensive functional and positional information for any one variant, was utilized. For variants in genes with multiple transcripts or at positions that correspond to multiple genes categories were assigned based only on the position within the first gene listed in the UCSC known genes database (Hsu F, Kent W J, Clawson H, Kuhn R M, Diekhans M, et al. (2006) The UCSC Known Genes. Bioinformatics 22: 1036-1046). In total 9,078,405 1KGP SNPs were annotated with positional categories. All positional categories were scored 0 or 1.
The following genic annotation categories were used:
10 k Up. This category consisted of all 1KGP SNPs that were between 10,000 and 1,001 base pairs upstream of the transcription start site for the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). For SNPs gene dense areas, priority was given to upstream category over downstream. Thus SNPs both 10,000 base pairs upstream and downstream from a protein coding gene were only annotated with the upstream category.
1 k Up. This category consisted of all 1KGP SNPs that were between 1,000 and 1 base pair(s) upstream of the transcription start site for the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). For SNPs gene dense areas, priority was given to upstream category over downstream. Thus SNPs both 1,000 base pairs upstream and downstream from a protein coding gene were only annotated with the upstream category.
5′UTR. This category consisted of all 1KGP SNPs that were located within the five prime untranslated region (5′UTR) of the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). All regions that are transcribed, but not translated, are assigned to UTR categories. If a polymorphism was within an exon or intron within a 5′UTR, it was annotated only as 5′UTR.
Exon. This category consisted of all 1KGP SNPs that were located within an exon of the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). If a polymorphism was within an exon that fell within the 5′UTR or 3′UTR of a gene, it was annotated only as 5′UTR or 3′UTR.
Intron. This category consisted of all 1KGP SNPs that were located within an intron of the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). If a polymorphism was within an intron that fell within the 5′UTR or 3′UTR of a gene, it was annotated only as 5′UTR or 3′UTR.
3′UTR. This category consisted of all 1KGP SNPs that were located within the three prime untranslated region (3′UTR) of the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). All regions that are transcribed, but not translated, are assigned to UTR categories. If a polymorphism was within an exon or intron within a 3′UTR, it was annotated only as 5′UTR.
1 k Down. This category consisted of all 1KGP SNPs that were between 1 and 1,000 base pair(s) downstream of the transcription start site for the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). For SNPs gene dense areas, priority was given to upstream category over downstream. Thus SNPs both 1,000 base pairs upstream and downstream from a protein coding gene were only annotated with the upstream category.
10 k Down. This category consisted of all 1KGP SNPs that were between 1,001 and 10,000 base pair(s) downstream of the transcription start site for the primary listing of protein coding genes in the UCSC known genes database (Hsu et al., supra). For SNPs gene dense areas, priority was given to upstream category over downstream. Thus SNPs both 10,000 base pairs upstream and downstream from a protein coding gene were only annotated with the upstream category.
Additional categories were recorded, including 10,001-100,000 BP up and downstream of protein coding genes, presence within a non-coding RNA, presence within a transcription factor binding site, and presence within a microRNA binding site. These categories were used to help select intergenic SNPs but were not analyzed in terms of differential enrichment (see discussion below).
The above positional annotations were leverages in the densely mapped 1KGP to characterize the types of variants that each GWAS studied SNP was a surrogate for, or tagged, as a result of Linkage Disequilibrium (LD). Each GWAS performed quality control according to best practices, as describes in detail in each of the original publications (See above). GWAS SNPs with reference SNP (rs) numbers that did not map to the 1KGP were excluded.
In order to assign LD-weighted annotation scores, a correlation coefficient approximation to r2 pairwise linkage disequilibrium (LD) was calculated using Plink version 1.07 (Purcell et al., supra). For each GWAS tag SNP present in the 1KGP pairwise LD was calculated to all other 1KGP SNPs within 1,000,000 base pairs (1 Mb) on either side of the SNP. This provided, for each SNP, a 2 Mb window in which LD scores were considered. LD scores were thresholded at r2≧0.2. LD scores were continuous valued from 0.2 to 1. Each SNP was assigned an LD value of 1 with itself (The robustness of the results to these parameter settings is discussed below in the section entitled Robustness of LD Weighted Scoring Procedure).
For each GWAS tag SNP, continuous, non-exclusive LD-weighted category scores were assigned as the LD weighted sum of the positional category scores for variants tagged in each of the eight categories mentioned above as annotated in the 1KGP reference panel. Summary statistics describing the distribution of scores in each category for the 2,558,411 SNPs representing the union of all GWAS considered are provided in Table 12.
Intergenic SNPs were determined after LD-weighted scoring. They were defined by weighted LD scores for each of the eight categories equal to zero. In addition these SNPs did not tag any SNPs in the 1KGP reference panel located within 100,000 base pairs of a protein coding gene, within a noncoding RNA, within a transcription factor binding site nor within a microRNA binding site.
For comparison and to assess the effect of leveraging LD weighted scoring in this way comparisons were made between LD-weighted scores (
The empirical null distribution in GWAS is affected by global variance inflation due to population stratification and cryptic relatedness (Devlin B, Roeder K (1999) Genomic control for association studies. Biometrics 55: 997-1004) and deflation due to over-correction of test statistics for polygenic traits (Yang J, Weedon M N, Purcell S, Lettre G, Estrada K, et al. (2011) Genomic inflation factors under polygenic inheritance. Eur J Hum Genet 19: 807-812) by standard genomic control methods. A control method leveraging only intergenic SNPs which are likely depleted for true associations was applied. All p-values were converted into z-scores, and, for each phenotype, the genomic inflation factor (Devlin et al., supra), λGC, was estimated for intergenic SNPs. All test statistics were divided by λGC. The inflation factor, λGC was computed as the median z-score squared divided by the expected median of a chi-square distribution with one degree of freedom for all phenotypes except CPD, where the 0.95 quantile was used in place of the median. For correction statistics see Table 14.
The intergenic SNPs were leveraged to estimate inflation because their relative depletion of associations suggests they provide a robust estimate of true null SNPs that is uncontaminated by polygenic effects. Using annotation categories in this fashion is important given concerns posed by recent GWAS about the over-correction of test statistics using standard genomic control. Statistics from this procedure are shown in Table 14. The traditional GC value for the summary statistics from each GWAS in their received state are reported. Original values less than 1.0 suggest an over correction by traditional GC metrics, while values greater than 1.0 suggest an under correction or no correction at all. The values that remain after intergenic inflation correction are likely to represent variance inflation due to true polygenic effects.
Q-Q plots are standard tools for assessing similarity or differences between two cumulative distribution functions (cdfs) (Schweder T, Spjotvoll E (1982) Plots of P-Values to Evaluate Many Tests Simultaneously. Biometrika 69: 493-502). When the probability distribution of GWAS summary statistic two-tailed p-values is of interest, under the global null hypothesis the theoretical distribution is uniform on the interval [0,1]. If nominal p-values are ordered from smallest to largest, so that p(1)<p(2)< . . . <p(N), the corresponding empirical cdf, denoted by “q”, is simply q(i)=i/N (in practice adjusted slightly to account for the discreteness of the empirical cdf), where N is the number of SNPs in the GWAS (or genic category). Thus, for a given index i, the x-coordinate of the Q-Q curve is simply q(i), since the theoretical inverse cdf is the identity function, and the y-coordinate is simply the nominal p-value p(i). As is common practice in GWAS, −log10 p is plotted against the −log10 q to emphasize tail probabilities of the theoretical and empirical distributions; these coordinates are labeled “nominal −log10 (p)” and “empirical −log10 (q)” in the Q-Q plots. For a given threshold of GC-controlled p-values, category ‘enrichment’ is seen as a horizontal (not vertical) deflection of the Q-Q curve from the identity line (or from one genic category to another) as described in detail next.
The ‘enrichment’ seen in the Q-Q plots can be directly interpreted in terms of False Discovery Rate (FDR)[18]. For a given p-value cutoff, the Bayes FDR (Efron B (2010) Large-scale inference: empirical Bayes methods for estimation, testing, and prediction. Cambridge; New York: Cambridge University Press. xii, 263 p. p) is defined as
FDR(p)=π0F0(p)/F(p), [S1]
where π0 is the proportion of null SNPs, F0 is the null cumulative distribution function (cdf), and F is the cdf of all SNPs, both null and non-null; see below for details on this simple mixture model formulation (Efron B (2007) Size, power and false discovery rates. The Annals of Statistics 35: 1351-1377). Under the null hypothesis, F0 is the cdf of the uniform distribution on the unit interval [0,1], so that Eq. [S1] reduces to
FDR(p)=π0p/F(p), [S2]
The cdf F can be estimated by the empirical cdf q=Np/N, where Np is the number of SNPs with p-values less than or equal to p, and N is the total number of SNPs. Replacing F by q in Eq. [S2]
FDR(p)≈π0p/q, [S3]
which is biased upwards as an estimate of the FDR[20]. Replacing no in Equation [S3] with unity gives an estimated FDR that is further biased upward;
FDR(p)≈p/q [S4]
If π0 is close to one, as is likely true for most GWAS, the increase in bias from Eq. [S3] is minimal. The quantity 1−p/q, is therefore biased downward, and hence a conservative estimate of the True Discovery Rate (TDR, equal to 1-FDR). Given the −log10 of the Q-Q plots
−log10(FDR(p))≈log10(q)−log10(p) [S5]
demonstrating that the (conservatively) estimated FDR is directly related to the horizontal shift of the curves in the Q-Q plots from the expected line x=y, with a larger shift corresponding to a smaller FDR. As before, the estimated true discovery rate can be obtained as one minus the estimated FDR. For each TDR plot in
After appropriate genomic control enrichment can be assessed by its genic category-specific TDR for a given z-score (equivalently, nominal p-value). Categories of SNPs that have a higher TDR for a given nominal p-value are more “enriched” than categories of SNPs with a lower TDR for the same nominal p-value. This measure of enrichment depends on choice of p-value threshold.
An overall single number summary of category-specific enrichment is the sample mean of z minus one, where the mean is taken over all SNP z-scores in the given category. Both the TDR and the mean (z2)−1 are justified as measures of enrichment based on a simple Bayesian mixture model framework. Specifically, let f(z) be the probability density for the SNP summary statistic z-scores. This is modeled as the mixture of a null probability density f0 and a non-null density f1
f(z)=π0f0(z)+π1f1(z), [S6]
where, as above, π0 is the proportion of SNPs with no association with the trait and π1=1−π0 the proportion of SNPs with a non-zero association with the trait. Assuming that the z-scores are symmetric about zero, the variance of this distribution is
∫z2f(z)dz=∫z2π0f0(z)dz+∫z2π1f1(z)dz=π0+π1∫z2f1(z)dz, [S7]
since the variance of the null distribution is one after appropriate genomic control. Under the assumption that the proportion of null SNPs (π0) is close to one, a mildly conservative estimate of the excess in variance attributable to non-null SNPs is given by ∫z2 f(z) dz−1. An unbiased estimate of this expression is the sample mean of z2 minus 1. Note, non-null z-scores are scaled by the square root of the sample size, and hence mean(z2)−1 is proportional to, not identical with, π1 times the tagged phenotypic variance of the non-null SNPs.
Consistency with Local False Discovery Rate Estimates
Under scenarios of multiple testing, such as GWAS, quantitative estimates of likely true associations can be estimated from the distributions of summary statistics. Efron (Efron B (2010) Large-scale inference: empirical Bayes methods for estimation, testing, and prediction. Cambridge; New York: Cambridge University Press. xii, 263 p. p) has developed a flexible framework for quantitatively estimating the null, non-null and mixture distributions from the resulting test statistics. Similar approaches have been applied in other fields, most relevantly to gene expression array data (Allison D B, Gadbury G L, Heo M S, Fernandez J R, Lee C K, et al. (2002) A mixture model approach for the analysis of microarray gene expression data. Computational Statistics & Data Analysis 39: 1-20) and linkage analysis (Ginns E I, St Jean P, Philibert R A, Galdzicka M, Damschroder-Williams P, et al. (1998) A genome-wide search for chromosomal loci linked to mental health wellness in relatives at high risk for bipolar affective disorder among the Old Order Amish. Proc Natl Acad Sci USA 95: 15531-15536). As a demonstration, the CD statistics were fit using this model (
The empirical Bayesian modeling approach described by Efron (2010; supra) is implemented in the freely available R package locfdr (Efron B, Turnbull B B, Narasimhan B (2011) locfdr: Computes local false discovery rates). The approach is to model the mixture density of effects in terms of z-scores as in Eq. [S6] above, or as a mixture density consisting of a weighted linear combination of a null density f0(z) for the z-scores of SNPs with no association, and a non-null density f1(z) for z-scores from trait-associated SNPs. The local false discovery rate (locfdr) is then given by
locfdr(z)=π0f0(z)/f(z), [S8]
where f(z) is given by Eq. [S6]. Using this model, the empirical null density (assumed to be normal, with mean 0 and data determined standard deviation) was estimated. The null for intergenic SNPs was estimated and all statistics were adjusted accordingly such that the intergenic test statistics conformed to the theoretical distribution (normal with mean 0 and standard deviation 1). This approach mirrors the intergenic inflation control described previously. The locfdr library was used to estimate the mixture density, fixing the null distribution to the theoretical standard normal and estimating the mixture density non-parametrically as a smoothed histogram. This model was fit to the overall data and per category (
This framework also allows us to estimate the a posteriori expected z-scores, as described in chapter 11, pp. 218 of (Efron, 2010; supra), based on the nonparametric estimates of the mixture density f(z) (Eq. [S6]) obtained with locfdr. For each of the 70 discovery sets used to calculate cumulative replication rates, the expected a posteriori effect size across the same 120 equally sized z-score bins ranging from −5.33 to 5.33 (corresponding to the GWAS p-value of 5×10−8) were calculated. The results were averaged across the 70 iterations and plotted as a function of discovery z-score independently for each genic annotation category. Because the direction of effect (z-score sign) is arbitrary with respect to the allele and strand chosen as causal, the data were duplicated with opposite sign to enforce symmetry. Again this procedure was carried out for the overall data and per category (
For comparison, empirical replication z-scores were calculated using the same 70 discovery-replication pairs and averaged across iterations. For visualization a cubic smoothing spline was fit relating the discovery z-score bin midpoints to the corresponding average replication z-scores. The empirical z-score replications (
In addition to the the non-parametric approach to estimating the mixture model (Eq. [S6]) implemented in the locfdr package, a parametric model was implimate, to facilitate simulations and extensions of the basic locfdr model to include covariates, described below. Specifically, w=−2 ln(p) was modeled as a mixture of a (null) χ2 density with two degrees of freedom and a (non-null) Weibull density with shape parameter a and scale parameter b. Note, under the null hypothesis the p-values are uniformly distributed and hence w has a χ2 density with two degrees of freedom (df), equivalent to a Weibull density with a=1 and b=2. Hence, the mixture density for w is given by
f(w)=π0f0(w)+π1f1(w), [S9]
where f0(w) is Weibull(a0=1, b0=2) and f1(w) is Weibull(a1, b1), where the parameters (π0, a1, b1) are estimated from the data. For identifiability, the model is fit under the assumption (in common with the locfdr package) that the non-null density is zero in a small interval around zero, accomplished here by shifting f1 to the right by a fixed margin, e.g., the median of the χ2 distribution with 2 df. This is equivalent to the assumption that the vast majority of SNPs with z-scores close to zero are true nulls[19]. For parameter estimation, a Bayesian Monte Carlo Markov Chain (MCMC) algorithm was used, placing vague priors on the parameters (π0, a1, b1). Q-Q plots and model fits for Height and CD for SNPs below the GWAS-level significance threshold of 5×10−8 are given in
The CD parameter estimates were used to determine the impact of sample size and polygenicity on Q-Q plots and enrichment indices in the context of mixture models.
The basic parametric mixture model [S9] was extended by allowing for covariates (e.g., genic annotations). Specifically, let x be a vector of annotations for a given SNP. The covariate-modulated mixture model is given by
f(w|x)=π0(x)f0(w)+π1(x)f1(w|x), [S10]
where π0(x)=1/(1+exp(x′ν)) is a logistic function of the covariates, and f1(w|x) is a Weibull distribution with shape parameter a=exp(x′α) and scale parameter b=exp(x′β). The model is estimated using an MCMC algorithm (Gibbs sampler with Metropolis-Hastings steps), placing non-informative priors on unknown parameters (ν, α, β). Estimates from this model, not presented here, could be used to replace the stratified FDR analyses in the main text by directly using Eq. [S10] to estimate the local fdr (Eq. [S8]). Control for potential confounds: LD and MAF
Significant categorical differences in terms of total LD and total number of SNPs captured by each GWAS SNP that mirrors the enrichment findings were observed (Tables 17 and 18). To rule out total LD as a potential confound, a multiple regression was performed on height GWAS summary values (log of z2 after intergenic inflation control) using SNP annotation category scores and total summed LD as predictors. Each category score is computed as described in the main text. The category score of each SNP is pre-multiplied by the genetic variance (MAF*(1−MAF)) of that SNP. Annotations categories were centered to have mean zero. The analysis reveals only a minor effect of total LD on predicting log(z2) and strong individual category effects which mirror the enrichment findings (Table 20).
Systematic differences in the average minor allele frequency (MAF) could confound enrichment analysis as MAF acts multiplicatively with effect size to give z-scores. The average minor allele frequency per category are shown in Table 19.
The estimated TDR can be thought of as the replication rate in an independent sample as the replication sample size goes to infinity. In practice, both the estimated TDR and the replication sample effect sizes will be measured with error, and hence the estimated TDR will not perfectly predict the independent sample replication rate. Nonetheless, there should be a close correspondence for reasonable discovery and replication sample sizes. Thus, to provide empirical support for the findings, category-specific rates of replication across eight truly independent GWAS samples studying CD were investigated. For each of eight sub-studies contributing to the final meta-analysis in the CD report, the reported z-scores were adjusted according to the intergenic inflation correction method described above. For each of the 70 (8 choose 4) possible combinations of four-study discovery and four-study replication sets, the four-study combined discovery z-score and four-study combined replication z-score for each SNP as the average z-score across the four studies was calculated, multiplied by the square root of the number of studies. For discovery samples the z-scores were converted to two-tailed p-values, while replication samples were converted to one-tailed p-values preserving the direction of effect in the discovery sample. Replication was defined as a one-tailed p-value less than 0.05 in the replication set. For each of the 70 discovery-replication pairs cumulative rates of replication were calculated over 1000 equally-spaced bins spanning the range of negative log10(p-values) observed in the discovery samples. The cumulative replication rate calculated for any bin was the total number of replicated SNPs (p<0.05, one-tailed test with direction of effect given by the discovery sample) with a negative log10(discovery p-value) greater than or equal to the lower bound of the bin divided by the total number of SNPs with a negative log10 (discovery p-value) greater than or equal to the lower bound of that bin. This analysis was repeated for each of the eight genic annotation categories as well as intergenic SNPs and all SNPs. The cumulative replication rates were averaged across the 70 discovery-replication pairs and the results are reported in
The original LD weighted annotation scoring approach (see: Linkage Disequilibrium (LD) Weighted Annotation Score above) only considered pairwise r2 LD greater than 0.2 and within 1 megabase of the target GWAS SNP. However, it is likely that true correlations exist at lower level than r2=0.2 and beyond 1 megabase. To test the dependence of the results upon the parameters of the scoring approach, each SNP was reclassified following the same procedure as before, but including estimated r2 LD greater than 0.05 and within 2 megabases. The pattern of enrichment described in the original stratified QQ-plots appears robust to these changes (
Results
LD Based Enrichment of Genic Elements in Height
Under multiple testing paradigms, such as GWAS, quantitative estimates of likely true associations can be estimated from the distributions of summary statistics [12,13]. A common method for visualizing the enrichment of statistical association relative to that expected under the global null hypothesis is through Q-Q plots of the nominal p-values resulting from GWAS. Under the global null hypothesis the theoretical distribution is uniform on the interval [0,1]. Thus, the usual Q-Q curve has as the y-coordinate the nominal p-value, denoted by “p”, and the x-coordinate the value of the empirical cdf at p, denoted by “q”. As is common in GWAS, −log 10 p is plotted against the −log10 q to emphasize tail probabilities of the theoretical and empirical distributions. In such plots, enrichment results in a leftward shift in the Q-Q curve, corresponding to a larger fraction of SNPs with nominal −log10 p-value greater than or equal to a given threshold.
The stratified Q-Q plot for height (
An earlier departure from the null line (leftward shift) suggests a greater proportion of true associations, for a given nominal p-value. The divergence of the curves for different categories thus suggests that the proportion of non-null effects varies considerably across annotation categories of genic elements. For example, the proportion of SNPs in the 5′UTR category reaching a given significance level (−log10(p)>10) is roughly 10 times greater than for all SNPs, and 50-100 times greater than for intergenic SNPs.
Polygenic Enrichment Across Diverse Phenotypes
Recently Yang et al [14] demonstrated that an abundance of low p-values beyond what is expected under null hypotheses in GWAS, but not necessarily reaching stringent multiple comparison thresholds, and often seen as ‘spurious inflation,’ can also be consistent with an enrichment of true ‘polygenic’ effects [14]. The prevalence of enrichment below the established genome-wide significance threshold of p<5×10−8 (−log10(p)>7.3;) in height (
The enrichment patterns among annotation categories are consistent across phenotypes, including schizophrenia (SCZ) and tobacco smoking (cigarettes per day; CPD;
Significance values were computed for the curves for each annotation category relative to those for intergenic SNPs, using a two-sample Kolmogorov-Smirnov Test. The enrichment for height was highly significant for all categories when compared with the intergenic category, with all p-values less than 2.2×10−16. Nearly all genic categories were also significantly enriched for all the other phenotypes (Table 15).
While the pattern of enrichment is consistent, the shape of the curves varies across phenotypes. In particular, the point at which the curves deviate from the expected null line occurs earliest for height, followed by SCZ, and finally CPD (
Intergenic Genomic Control
The relative absence of enrichment in intergenic SNPs indicates minimal inflation due to polygenic effects and a more robust estimate of the global null. This fact can be exploited for estimation of variance inflation due to stratification [15] that is minimally confounded by true polygenic effects [14], by confining the estimation of the genomic inflation factor [15], λGC, to only intergenic SNPs. Here, summary statistics were adjusted for all phenotypes according to this “intergenic inflation control” procedure.
Category Specific True Discovery Rate
Since specific genic tag SNP categories are significantly more likely to be associated with common phenotypes, while intergenic ones are less likely, all tag SNPs should not be treated as exchangeable. Variation in enrichment across diverse genic categories is expected to be associated with corresponding variation in true discovery rate TDR for a given nominal p-value threshold. A conservative estimate of the TDR for each nominal p-value is equivalent to 1−(p/q) as plotted on the Q-Q plots. This relationship is shown for height, SCZ and CPD (
Quantification of Enrichment
While the TDR provides a quantification of enrichment for a given nominal p-value threshold (equivalently, SNP z-score threshold), a single number quantification of enrichment for each LD-weighted annotation category within each phenotype, computed as the sample mean (z2)−1 is provided. The sample mean, taken over all SNPs in a given category, provides an estimate of the variance due to null and non-null SNPs; by subtracting one can obtain a conservative estimate of the variance in effect sizes attributable to non-null SNPs alone. Both TDR and mean (z2)−1 are justified based on a standard mixture model formulation. These enrichment scores, normalized by the maximum value across categories within each phenotype, are presented in
Categories where each SNP, on average, tags more SNPs or represents a larger total amount of LD could spuriously appear enriched. Categorical differences in the number of SNPs and total summed LD captured by each SNP were observed but multiple regression shows the effect is negligible and independent categorical effects persist despite the significant correlation among categories. Likewise, systematic deviations in minor allele frequency (MAF) across categories could bias annotation category effects as MAF acts multiplicatively with effect size to explain variance. Minimal categorical stratification was found for MAF not consistent with it driving the enrichment findings. To further address the possibility that some of the differential enrichment of categories could be due to category-specific genomic inflation from the above factors, null-GWAS simulations based on genotypes from the 1000 Genome Project were performed. The results indicate that such effects are non-existent or negligible.
Replication Rate
To further address the possibility that the observed pattern of differential enrichment results from spurious (e.g., non-generalizable) associations due to category-specific confounding effects or statistical modeling errors, the empirical replication rate across independent sub-studies for one phenotype (CD), for which the required sub-study summary statistics were available was studied.
Increased Power Using Stratified False Discovery Rates
In order to demonstrate the utility of the enriched category information for improved discovery, an established method for computing stratified False Discovery Rates [9] was utilized. The sFDR method extends the traditional methods for FDR control [21], improving power by taking advantage of pre-defined, differentially enriched strata among multiple hypothesis testing p-values. Here, an increase in power from using stratified (vs. unstratified) methods is defined as a decreased Non-Discovery Rate (NDR) for a given level of FDR control α, where NDR is the proportion of false negatives among all tests [22]. Specifically, the ratio of NDR from stratified FDR control vs. NDR was estimated from unstratified FDR control. A ratio above one is equivalent to sFDR rejecting more SNPs than unstratified FDR for a common level α.
For each phenotype, the SNPs are divided into independent strata according to their predicted tagged variance (z2) based on a linear regression predictor with regression weights for each annotation category trained using the height GWAS summary statistics. An increase in the number of discovered SNPs was observed. For example, for α=0.05 the increased proportion of declared non-null SNPs using sFDR ranges from 20% in height to 300% in schizophrenia. Leveraging the genic annotation categories in the sFDR framework provides one possible avenue for improving the output of likely non-null SNPs in GWAS by taking advantage of the non-exchangeability of SNPs demonstrated by the genic annotation category enrichment analyses.
1.98 × 10−323
5.77 × 10−131
6.71 × 10−240
enrichment values used to create FIG. 2 and FIG. 27. All values are expressed in relative proportions highest category for each phenotype. BD, Bipolar Disorder; BMI, Body Mass Index; CD, Crohn's disease;
Cigarettes per Day; DBP, Diastolic blood pressure; HDL, High density lipoprotein; LDL, Low d
lipoprotein; SBP, systolic blood pressure; SCZ, Schizophrenia; TC, total Cholesterol; TG, triglycerides
Ulcerative Colitis; WHR, Waist-hip-ratio.
indicates data missing or illegible when filed
Participant Samples
Complete GWAS results in the form of summary statistics p-values were obtained from public access websites or through collaboration with investigators (T2D cases and controls from the DIAGRAM Consortium and schizophrenia cases and controls from the Psychiatric GWAS Consortium (PGC)—Table 25). There was no overlap among participants in the CVD GWAS and the schizophrenia case-control sample (n=21,856), except for 2,974 of 12,462 controls (24%)137. The schizophrenia GWAS summary statistics results were obtained from the Psychiatric GWAS Consortium (PGC)13, which consisted of 9,394 cases with schizophrenia or schizoaffective disorder and 12,462 controls (52% screened) from a total of 17 samples from 11 countries. The quality of phenotypic data was verified by a systematic review of data collection methods and procedures at each site, and only studies that fulfilled these criteria were included. This involved nine key items: i) the use of a structured psychiatric interview, ii) systematic training of interviewers in the use of the instrument, iii) systematic quality control of diagnostic accuracy, iv) reliability trials, v) review of medical record information, vi) best-estimate procedure employed, vii) specific inclusion and exclusion criteria developed and utilized, viii) MDs or PhDs as making the final diagnostic determination, and ix) special additional training for the final Schizophrenia PGC. One sample from Sweden used another approach, but further empirical support for the validity of this approach was provided. Controls consisted of 12,462 samples of European ancestry collected from the same countries. As the prevalence of schizophrenia is low, a large control sample where some controls were not screened for schizophrenia was utilized. For further details on sample characteristics and quality control procedures applied, please see Ripke et al 13. There were 2974 controls in the schizophrenia UK case control sample from the Welcome Trust Case Control Consortium that were also included in several of the CVD risk factor GWAS. This constitutes 24% of the total number of controls (n=12,462) in the Schizophrenia PGC sample13. More information about inclusion criteria and phenotype characteristics of the Cardiovascular Disease (CVD) risk factors samples of the different GWAS are described in the original publications 29-33. The relevant institutional review boards or ethics committees approved the research protocol of the individual GWAS used in the current analysis and all human participants gave written informed consent.
Statistical Analyses
Stratified Q-Q Plots
Q-Q plots compare a nominal probability distribution against an empirical distribution. In the presence of all null relationships, nominal p-values form a straight line on a Q-Q plot when plotted against the empirical distribution. For each phenotype, for all SNPs and for each categorical subset (strata), −log 10 nominal p-values were plotted against −log 10 empirical p-values (stratified Q-Q plots). Leftward deflections of the observed distribution from the projected null line reflect increased tail probabilities in the distribution of test statistics (z-scores) and consequently an over-abundance of low p-values compared to that expected by chance, also termed “enrichment”. Under large-scale testing paradigms, such as GWAS, quantitative estimates of likely true associations can be estimated from the distributions of summary statistics36; 37. A common method for visualizing the enrichment of statistical association relative to that expected under the global null hypothesis is through Q-Q plots of nominal p-values obtained from GWAS summary statistics. The usual Q-Q curve has as the y-ordinate the nominal p-value, denoted by “p”, and as the x-ordinate the corresponding value of the empirical cdf, denoted by “q”. Under the global null hypothesis the theoretical distribution is uniform on the interval [0.1]. As is common in GWAS, −log 10p is plotted against −log 10 q to empha 1 size tail probabilities of the theoretical and empirical distributions. Therefore, genetic enrichment results in a leftward shift in the Q-Q curve, corresponding to a larger fraction of SNPs with nominal −log 10 p-value greater than or equal to a given threshold. Stratified Q-Q plots are constructed by creating subsets of SNPs based on levels of an auxiliary measure for each SNP, and computing Q-Q plots separately for each level. If SNP enrichment is captured by variation in the auxiliary measure, this is expressed as successive leftward deflections in a stratified Q-Q plot as levels of the auxiliary measure increase.
Genomic Control
The empirical null distribution in GWAS is affected by global variance inflation due to population stratification and cryptic relatedness38 and deflation due to over-correction of test statistics for polygenic traits by standard genomic control methods39. A control method leveraging only intergenic SNPs, which are likely depleted for true associations (Example 2), was applied. First, the SNPs were annotated to genic (5″UTR, exon, intron, 3″UTR) and intergenic regions using information from the 1000 Genomes Project (1KGP). As illustrated in
Stratified 1 Q-Q Plots for Pleiotropic Enrichment
To assess pleiotropic enrichment, Q-Q plot stratified by “pleiotropic” effects were used. For a given associated phenotype, enrichment for pleiotropic signals is present if the degree of deflection from the expected null line is dependent on SNP associations with the second phenotype. Stratified Q-Q plots of empirical quantiles of nominal −log10(p) values were constructed for SNP association with schizophrenia for all SNPs, and for subsets (strata) of SNPs determined by the nominal p-values of their association with a given CVD risk factor. Specifically, the empirical cumulative distribution of nominal p-values was computed for a given phenotype for all SNPs and for SNPs with significance levels below the indicated cut-offs for the other phenotype (−log10(p)≧0, −log10(p)≧1, −log10(p)≧2, −log10(p)≧3 corresponding to p<1, p<0.1, p<0.01, p<0.001, respectively). The nominal p-values (−log10(p)) are plotted on the y-axis, and the empirical quantiles (−log10(q), where q=1−cdf(p)) are plotted on the x-axis. To assess for polygenic effects below the standard GWAS significance threshold, the stratified Q-Q plots were focused on SNPs with nominal −log10(p)<7.3 (corresponding to p>5×10−8).
Stratified True Discovery Rate (TDR)
Enrichment seen in the stratified Q-Q plots can be directly interpreted in terms of TDR (equivalent to one minus the FDR40). The stratified FDR method35, previously used for enrichment of GWAS based on linkage information were applied 34. Specifically, for a given p-value cutoff, the FDR is defined as
FDR(p)=π0F0(p)/F(p), [1]
where π0 is the proportion of null SNPs, F0 is the null cdf, and F is the cdf of all SNPs, both null and non-null; see below for details on this simple mixture model formulation41. Under the null hypothesis, F0 is the cdf of the uniform distribution on the unit interval [0,1], so that Eq. [1] reduces to
FDR(p)=π0p/F(p), [2].
The cdf F can be estimated by the empirical cdf 1 q=Np/N, where Np is the number of SNPs with p2 values less than or equal to p, and N is the total number of SNPs. Replacing F by q in Eq. [2]
Estimated FDR(p)=π0p/q, [3],
which is biased upwards as an estimate of the FDR41. Replacing π0 4 in Equation [3] with unity gives an estimated FDR that is further biased upward; q*=p/q [4]. If no is close to one, as is likely true for most GWAS, the increase in bias from Eq. [3] is minimal. The quantity 1−p/q, is therefore biased downward, and hence is a conservative estimate of the TDR. Referring to the formulation of the Q-Q plots, q* is equivalent to the nominal p-value divided by the empirical quantile, as defined earlier. Given the −log 10 of the Q-Q plots
−log10(q)=log10(q)−log10(p) [5]
demonstrating that the (conservatively) estimated FDR is directly related to the horizontal shift of the curves in the stratified Q-Q plots from the expected line x=y, with a larger shift corresponding to a smaller FDR, as illustrated in
Stratified Replication Rate
For each of the 17 sub-studies contributing to the final meta-analysis in schizophrenia, z-scores were independently adjusted using intergenic inflation control. For 1000 of the possible combinations of eight-study discovery and nine study replication sets, the eight-study combined discovery z-score and eight or nine-study combined replication z-score for each SNP as the average z-score across the eight or nine 1 studies, multiplied by two (the square root of the number of 2 studies). For discovery samples the z-scores were converted to two-tailed p-values, while replication samples were converted to one-tailed p-values preserving the direction of effect in the discovery sample. For each of the 1000 discovery-replication pairs cumulative rates of replication were calculated over 1000 equally-spaced bins spanning the range of negative log 10(p-values) observed in the discovery samples. The cumulative replication rate for any bin was calculated as the proportion of SNPs with a −log 10(discovery p-value) greater than the lower bound of the bin with a replication p-value<0.05. Cumulative replication rates were calculated independently for each of the four pleiotropic enrichment categories as well as intergenic SNPs and all SNPs. For each category, the cumulative replication rate for each bin was averaged across the 1000 discovery-replication pairs and the results are reported in
Stratified Replication Effect Sizes
Stratified TDR is directly related to stratified replication effect sizes and hence replication rates. As before, for each of the 17 sub-studies contributing to the final meta-analysis in schizophrenia z-scores were independently adjusted using intergenic inflation control. For 1000 of the possible combinations of eight-study discovery and nine study replication sets, the eight-study combined discovery z-score and eight or nine-study combined replication z-score were calculated for each SNP. The effect sizes were stratified by levels of log 10(p-values) from the triglycerides GWAS. For visualization, a cubic smoothing spline was fit relating the discovery z-score bin midpoints to the corresponding average replication z-scores (see
Conditional Statistics—Test of Association with Schizophrenia
To improve detection of SNPs associated with schizophrenia, a stratified FDR approach was used, leveraging pleiotropic phenotypes using established stratified FDR methods34; 35. Specifically, SNPs were stratified based on p-values in the pleiotropic phenotype (e.g. Triglycerides; TG). A conditional FDR value (denoted as FDR SCZ|TG) for schizophrenia (SCZ) was assigned to each SNP, based on the combination of p-value for the SNP in schizophrenia and the pleiotropic trait, by interpolation into a 2-D look-up table (
Conditional Manhattan Plots
To illustrate the localization of the genetic markers associated with schizophrenia given the CVD risk factor effect, a “Conditional Manhattan plot” was used, plotting all SNPs within an LD block in relation to their chromosomal location. As illustrated in
Conjunction Statistics 1—Test of Association with Both Phenotypes
In order to identify which of the SNPs associated with schizophrenia given the CVD risk factor (SCZ|CVD, Table 23) were also associated with CVD risk factors given schizophrenia (opposite direction), the conditional FDR was calculated in the other direction (CVD|SCZ). This is reported in in Table 24. The corresponding z-scores are listed in Table 27. The z-scores were calculated from the p-values and the direction of effect was determined by the risk allele. In addition, to make a comprehensive, unselected map of pleiotropic signals, a conjunction testing procedure was used, as outlined for p-value statistics in Nichols et al.42 and adapted this method for FDR statistics based on the conditional FDR approach34; 35. The conjunction statistics (denoted as FDR SCZ & TG) were defined as the max of the conditional FDR in both directions, i.e. FDR SCZ & TG=max(FDR SCZ|TG, FDR TG|SCZ) based on the combination of p-value for the SNP in schizophrenia and the pleiotropic trait, by interpolation into a bidirectional 2-D look-up table (
Conjunction Manhattan Plots
To illustrate the localization of the pleiotropic genetic markers association with both schizophrenia and CVD risk factors, a “Conjunction Manhattan plot” was used, plotting all SNPs with a significant conjunction FDR within an LD block in relation to their chromosomal location. As illustrated in
Results
Q-Q Plots of Schizophrenia SNPs Stratified by Association with Pleiotropic CVD Risk Factors
Stratified Q-Q plots for schizophrenia conditioned on nominal p-values of association with triglycerides (TG) showed enrichment across different levels of significance for TG (
Conditional True Discovery Rate (TDR) in Schizophrenia is Increased by CVD Risk Factors
Since categories of SNPs with stronger pleiotropic enrichment are more likely to be associated with schizophrenia, to maximize power for discovery, all tag SNPs should not be treated exchangeably. Specifically, variation in enrichment across pleiotropic categories is expected to be associated with corresponding variation in the TDR (equivalent to 1-FDR)40 for association of SNPs with schizophrenia. A conservative estimate of the TDR for each nominal p-value is equivalent to 1−(p/q), obtained from the stratified Q-Q plots. This relationship is shown for schizophrenia conditioned on TG (
Replication Rate in Schizophrenia is Increased by Pleiotropic CVD Risk Factors
To demonstrate that the observed pattern of differential enrichment does not result from spurious (e.g., non-generalizable) associations due to category-specific stratification or errors in statistical modeling, the empirical replication rate across independent sub-studies for schizophrenia was studied.
Schizophrenia Gene Loci Identified with Conditional FDR
To identify SNPs associated with schizophrenia, a “conditional” Manhattan plot was constructed for schizophrenia showing the FDR conditional on each of the CVD risk factors (
Pleiotropic Gene Loci in Schizophrenia and CVD Risk Factors Identified with Conjunction FDR
As a secondary analysis, it was investigated if any of the SNPs associated with schizophrenia conditioned on CVD (SCZ|CVD) were also significantly associated with CVD risk factors conditioned on SCZ (CVD|SCZ), i.e. the 1 conditional FDR in the opposite direction. 10 independent loci (pruned based on LD>0.2) were identified with a significant association also with the CVD risk factor (conditional FDR<0.05), including 3 complex loci, and 7 single gene loci. Of these, the ITIH4 region (3p21.1), and the CNNM2/NT5C2 region (10q24.32), in addition to the HLA region (chr. 6) have been identified in previous schizophrenia studies after including large replication samples 13. The significant loci were found in the TG|SCZ (6 loci), LDL|SCZ (3 loci), HDL|SCZ (4 loci), SBP|SCZ (2 loci), BMI|SCZ (1 locus) and WHR|SCZ (4 loci), and 6 loci were jointly associated with schizophrenia and more than one CVD risk factor (Table 24). This indicates that overlapping genetic pathways are involved in schizophrenia and CVD risk factors. The direction of the different SNP associations (z-scores) is shown in Table 27. There was no clear evidence for systematic directions across all the SNPs in the different phenotypes, probably due to complex LD structures, especially on chromosome 6.
Further, to provide a comprehensive, unselected map of pleiotropic loci between schizophrenia and CVD risk factors in addition to those primarily associated with schizophrenia a conjunction FDR analysis was performed and a “conjunction” Manhattan plot was constructed. 26 independent pleiotropic loci were identified (pruned based on LD>0.2, black line around large circles) with a significance threshold of conjunctional FDR<0.05, located on a total of 14 chromosomes. See Table 28 for more details.
These methods have been described in detail in a series of studies investigating psychiatric 11-13 and nonpsychiatric disorders.13,14
Q-Q Plots and False Discovery Rates
Q-Q plots are standard tools for assessing similarity or differences between two cumulative distribution functions (CDFs). When the probability distribution of GWAS summary statistic two-tailed P values is of interest, under the global null hypothesis, the theoretical distribution is uniform on the interval [0,1]. If nominal P values are ordered from smallest to largest, so that P(1)<P(2)< . . . <P(N), the corresponding empirical CDF, denoted by “Q,” is simply Q(i)=i/N (in practice, adjusted slightly to account for the discreteness of the empirical CDF), where N is the number of SNPs in the GWAS (or genic category). Thus, for a given index i, the x-coordinate of the Q-Q curve is Q(i) (since the theoretical inverse CDF is the identity function) and the y-coordinate is the nominal P value P(i). It is a common practice in GWAS to instead plot −log 10 P against the −log 10 Q to emphasize tail probabilities of the theoretical and empirical distributions. For a given threshold of genomic control-corrected P values, “enrichment” is seen as a horizontal deflection of the Q-Q curves from the identity line.
Enrichment seen in the Q-Q plots can be directly interpreted in terms of false discovery rate (FDR). For a given P value cutoff, the Bayes FDR, defined as the posterior probability of a given SNP is null, given its observed P value, is given by:
FDR(P)=π0F0(P)/F(P), (1)
where π0 is the proportion of null SNPs, F0 is the CDF under the null hypothesis, and F is the CDF of all SNPs, both null and non-null. Here, F0 is the CDF of the uniform distribution on the unit interval [0,1], and F(P) can be estimated with the empirical CDF Q, so that an estimate of equation (1) is given by:
FDR(P)≈π0·P/Qt, (2)
which is biased upwards as an estimate of the FDR. Setting π0=1 in equation (2), an estimated FDR is further biased upward; if π0 is close to 1, as is likely true for most GWAS, the increase in bias from equation (2) is minimal. The quantity 1−P/Q is, therefore, biased downward, and hence a conservative estimate of the true discovery rate (equal to 1 FDR). Given the −log 10 of the Q-Q plots:
−log10(FDR(P))≈log10(Q)−log10(P), (3)
demonstrating that the (conservatively) estimated FDR is directly related to the horizontal shift of the curves in the Q-Q plots from the expected line x=y, with a larger leftward shift corresponding to a smaller FDR.
Conditional Q-Q Plots and FDR. The Conditional
FDR as the posterior probability that a SNP belonging to a category c is null for a phenotype, given a P value as small as the observed P value. Formally, this is given by:
FDR(P|c)=π0(c)·P/F(P|c), (4)
where P is the P value for the phenotype, c=1, . . . , C is one of C possible categories, F(P|c) is the conditional CDF, and π0(c) is the proportion of null SNPs in category c. A conservative estimate of FDR(P|c) is produced by setting π0(c)=1 and using the empirical conditional CDF in place of F(P1|c) in equation (4). This is a straightforward generalization of the empirical Bayes approach developed by Efron.10
In terms of Q-Q plots, enrichment of category c2 compared with category c1 is expressed as a leftward deflection of the Q-Q curve for category c2 compared with c1. Given equation (3), this is equivalent to showing that the conditional FDR is smaller for SNPs in category c2 compared with c1 for the same P value, ie, FDR(P|c2)<FDR(P|c1). Thus, by choosing a priori categories that result in differentially enriched samples, a larger proportion of SNPs can be discovered for a given FDR threshold than can be obtained from typical (unconditional) FDR or P value-based analyses.
Covariate-Modulated FDR
Using summary statistics derived from SNP associations of huge GWAS, it was shown that functional genic elements show differential contribution to phenotypic variance, with some categories (eg, regulatory elements and exons) showing strong enrichment (ie, more likely to have an effect) for phenotypic association.13 The enrichment of SNPs in genic elements of the genome (the 5′UTR and 3′UTR regions) was present across a wide spectrum of complex phenotypes and traits, including SCZ.13 This shows that SNPs in 5′UTR, in particular, but also in exons and 3′UTR regions are more likely to be involved in susceptibility to SCZ. This information can be used in Bayesian statistical models to enhance gene discovery by including information on the genic region in which each SNP is located, as this indicates how likely it is for each SNP to have an effect. By applying this approach to data from the Psychiatric Genomics Consortium (PGC) SCZ sample,16 the power for detecting small genetic effects was improved, leading to discovery of new susceptibility loci that did not reach threshold of significance in traditional GWAS analyses.13
Empirical independent replication remains the gold standard for confirming statistical findings. The replication rates, defined as proportion of SNPs declared significant in training samples with P values below a given threshold in the replication sample and with z-scores with the same sign in both discovery and replication samples were tested in independent SCZ substudies from the PGC17 and it was found that annotation categories with the greatest enrichment (5′UTR, exons, 3′UTR) showed the highest replication rate for a given nominal P value, confirming that the observed enrichment is due to true associations and not to inflation due to population stratification or other potential sources of spurious effects (
In order to illustrate the increased sensitivity and specificity for gene discovery, the publically available PGC SCZ sample was utilized.16 Applying the CMFDR method to the PGC SCZ sample, a total of 86 gene loci (CMFDR<0.05) were identified. By computing a posteriori effect sizes from the CMFDR model, it is expected that a very large proportion of these loci will replicate in a SCZ GWAS of similar size.
Gene Discovery Due to Pleiotropy Enrichment
The small number of genes relative to the vast number of human phenotypes necessitates pleiotropy—the influence of one gene or haplotype on two or more distinct phenotypes. The value of pleiotropy for improved understanding of disease pathogenesis and classification, identification of new molecular targets for drug development, and genetic risk profiling have been recognized.18 But few studies have systematically investigated pleiotropy in human complex traits and disorders, and those that have have looked for pleiotropy only among SNPs that reach a threshold level of significance in one or both phenotypes.18 This approach fails to capitalize on the power inherent in pleiotropy to robustly detect weak genetic effects.
The pleiotropy approach described herein was used to assess the contribution of all SNPs from two independent GWAS to determine their common association with two distinct phenotypes. SCZ and bipolar disorder share several clinical phenotypes, and there is growing evidence indicating overlapping gene variants.6,16 This approach was used to increase gene discovery in these disorders, using two large GWAS from the PGC,6,16 where overlapping controls had been removed with same procedure as in the recent cross-disorder analysis.19 A very high degree of polygenic overlap between SCZ and bipolar disorder was discovered.12 This information was used to increase the power of the GWAS, by including level of pleiotropy as a factor in the statistical models. This resulted in an improved yield (sensitivity) of genes discovered for SCZ and bipolar disorder compared to standard methods at a given significance level (specificity). 12 Thus, by applying the pleiotropy enrichment method and leveraging the bipolar disorder GWAS, gene discovery in the SCZ GWAS was increased. Note, while the power to detect nonpleiotropic loci is not increased using the pleiotropy enrichment method, neither is power lost.
Simulations showed that a larger increase in gene discovery would occur, using standard GWAS approaches, if the SCZ sample was as large as the combined SCZ bipolar disorder GWAS.12 However, it is very expensive to recruit and genotype new samples; applying the new statistical tools to existing samples is a cost-efficient way to improve gene discovery.
The results also showed that an estimated 1.2% of all SNPs analyzed are pleiotropic for SCZ and bipolar disorder. With approximately 1 million SNPs analyzed, this means that there are approximately 12 000 SNPs involved. This is very similar to the estimate from a recent large SCZ GWAS.7 This quantification of the polygenicity further emphasizes that most of these variants must have very small effects.
The new statistical tools can also be used to investigate genetic overlap between SCZ and nonpsychiatric diseases and traits to gain more knowledge about shared genetic mechanisms. There is a well-known comorbidity between SCZ and cardiovascular risk factors, including obesity, hypertension, and dyslipidemia.20 For each of these phenotypes, results are available from large GWAS. The pleiotropy methods were used to investigate polygenic pleiotropy. A genetic overlap between SCZ and several cardiovascular risk factors, particularly blood lipids (cholesterol, triglycerides) was found. This enrichment was leveraged to boost gene discovery and identify several gene loci associated with SCZ,11 strongly indicating that common molecular genetic mechanisms are underlying some of the epidemiological relationships between SCZ and cardiovascular risk factors.
Immune factors have been implicated in SCZ. By investigating pleiotropy with multiple sclerosis, a demyelination disorder with clear evidence for involvement of immune genes, the statistical tools were applied to determine polygenic overlap. A strong genetic overlap between SCZ and multiple sclerosis were found 21 and several independent loci associated with SCZ were identified. In contrast, no genetic overlap was found between bipolar disorder and multiple sclerosis. Imputation of the major histocompatibility complex (MHC) alleles indicated opposite direction of effect in multiple sclerosis and SCZ. As most of the overlap between multiple sclerosis and SCZ was located in the MHC region, and there is previous evidence for large genetic overlap between bipolar disorder and SCZ, the findings indicate that the MHC region could differentiate between bipolar disorder and SCZ.
Polygenic Architecture: Implications for Disease Mechanisms and Clinic
The underlying biology of complex brain disorders such as SCZ remains mostly unknown. Structural magnetic resonance imaging (MRI) brain phenotypes are highly heritable (80%-90%),22 and a new cluster analytical method has shown how pleiotropic brain phenotypes cluster together.17 Previous work has shown how a selected number of SNPs can be used to identify genetically determined brain structure variation.23,24 Recent large meta-analysis showed how brain structure volumes can be successfully used in a GWAS, and SNPs associated with hippocampal volume were identified.25 By extending a twin study-based approach to a large MRI sample across different behavioral phenotypes, combined with the statistical framework for analysis of GWAS data to identify polygenic effects, it is possible to identify genetically determined brain substrates related to SCZ and core disease phenotypes.
Review of fdr
Efron and Tibshirani (2002) Efron and Tibshirani (2002) made the assumption that the test statistic zi, 1≦i≦n, has a different distribution based on whether the null hypothesis H0,i is true or false, where n is the total number of tests (SNPs). The non-null distribution will tend to have more extreme values of the test statistic. Hence, zi follows a two-group
mixture model f(zi)=π0f0(zi)+π1f1(zi), (1) where π0 is the proportion of true null hypotheses, π1=1−π0 is the proportion of true non-null hypotheses, f0 is the probability density function if H0 is true, and f1 is the probability density function if H0 is false. Local fdr is the posterior probability that the ith test is null given zi, which by
Bayes rule is given by
The null density was assumed to be standard normal (theoretical null) or normal with mean and variance estimated from the data (empirical null). The mixture density π0f0(z)+π1f1(z) (z) was estimated by fitting a high degree polynomial to histogram counts (Efron, 2010). If a set of SNPs are selected with an estimated fdr≦α for some αε 2 (0; 1), then on average (1−α)×100% of these will be true non-null SNPs.
Covariate-Modulated fdr
A set of external covariates observed for each hypothesis test may influence the distribution of the test statistic (Sun et al., 2006; Efron, 2010). Under this scenario, incorporating the covariate effects into fdr estimation can dramatically increase power for gene discovery. For example, the distribution of GWAS z-scores may depend on SNP-level functional annotations (Schork et al., 2013), pleiotropic relationships with related phenotypes (Andreassen et al.a, 2013; Andreassen et al.b, 2013), gene expression levels in certain tissues, evolutionary conservation scores, and so forth. These external covariates can be used to break the exchangeability assumption implicit in Eq. (1) and potentially increase the power for gene discovery over using standard local fdr given in Eq. (2).
Let xi=(1, x1i, x2i, . . . , xmi)T, where xi denotes an (m+1)-dimensional vector of covariates (including intercept) for the ith SNP. The cmfdr is defined as
where π1(xi)=1−π0(xi) is the prior probability that the ith test is non-null given xi and fi(zi|xi) is the non-null density of zi given xi. By Bayes' rule cmfdr is the posterior probability that the ith test is null given both zi and xi. It was assumed that the density under the null hypothesis does not depend on covariates. Both the probability of null status and the non-null density are allowed to depend on covariates, as described below.
Central to the estimation of the null proportion is the assumption that π0 is large (say greater than 0.90) and that the vast majority of SNPs with test statistics close to zero are in fact null. These assumptions are reasonable for GWA data (Hon-Cheong et al., 2010).
A Bayesian Two-Group Model
Summary statistics from GWAS are often made publicly available only as two-tailed p-values, and hence the magnitude of the z score is recoverable but not the sign. Moreover, the sign of the z score is a result of arbitrary allele coding. Hence, the mixture model was formulated for the absolute z-scores. The extension of the method to signed z-scores is straightforward. Folded Normal-Gamma Mixture Model The distribution of z under H0 is assumed to have the folded normal distribution, with null density f0(z)=φσ0(z)Iz≧0, where φ(z) is the normal density with mean zero and standard deviation σ0 and Iz≧0 is an indicator function which takes the value 1 when z≧0 and 0 otherwise. The density of z under the alternative hypothesis H1 is assumed to have a gamma distribution with shape parameter a(x) and rate parameter β.
Additionally, a location parameter μ>0 was specified to bound the nonnull gamma densities away from zero. The “zero assumption” of Efron (2007) states that the central peak of the z-scores consists primarily of null cases. Such an assumption is necessary to make the non-null distribution identifiable and for the MCMC sampling algorithm to converge. The assumption that the vast majority of SNPs with z-scores close to zero are null is already commonly made in GWAS. Hence, the location parameter μ=0.68 is set in the gamma distribution, corresponding to the median of the null density f0. All SNPs with absolute z-scores less than 0.68 are thus a priori considered null.
The mixture model formulation was completed by positing a latent indicator δ=(δ1, . . . , δn), where δi=1 if the ith SNP is non-null and zero otherwise. Then π1(xi) is the prior probability that δi=1 given covariates xi. The dependence of —1 on x is modelled via a logistic regression
where =z=(z1, . . . , zn)T is a vector of test statistics and X is a vector of unknown parameters.
The augmented likelihood function is then given by
where z=(z1, . . . , zn)T is the vector of test statistics and X is the n×(m+1) design matrix. Integrating out the latent indicators δ gives the mixture model corresponding to Eq. (3).
Prior Distributions Weakly-informative priors were applied to unknown parameters {β, α, γ, σ02}:
α˜N(0,Σα),
γ˜N(0,Σγ),
β˜Gamma(a0,b0),
σ02˜Inverse Gamma(aσ0,bσ0), (5)
0 g
where Σα and Σγ have large values on the diagonal, a0 and b= are shape and rate parameters of gamma distribution, and a—0 and b—0 are shape and scale parameters of inverse gamma distribution. Hyperparameters are fixed by the user. In the applications below, the dispersion matrices Σα and Σγ are set to be diagonal with variance 10,000; (a0; b0) and (a—0; b—0) were both set to (0.001,0.001).
Sampling Scheme The parameters sampled were α, β, γ and σ02 in turn from their full conditional distributions via a Gibbs sampler using Metroplis-Hastings (M-H) steps. Combining (4) and (5), the full conditional distributions are given by:
where I(β=0) is an indicator function and f(| . . . ) denotes the probability density of a parameter conditional on all other parameters and the data. The full conditional posteriors for α and γ in (6) do not take standard forms and are sampled using a multiple-try M-H sampler (Givens and Hoeting, 2005) with a multivariate t-distribution candidate. The full conditional for β has a gamma distribution and for σ02 an inverse gamma distribution, so that both can be sampled directly. Each iteration of the Gibbs sampler also includes generation of δ, with a Bernoulli full conditional distribution. For
One can obtain an a posteriori estimate of cmfdr(zi) for each zi as follows.
Assume that {(β(i), α(i), γ(i), σ02(i)) <1≦i≦L} from the posterior distribution of the parameters. For each draw 1
Then, for example, the posterior median of cmfdr(zi) can be estimated by taking the median of cmfdr(1)(zi) across all L posterior draws. The algorithm has been implemented in the R statistical package.
Simulation
Phenotypes were simulated under different settings of generative parameters from real genotype data obtained in n=3,719 healthy individuals. For each permutation of simulation settings 100 unique phenotypes were generated. The simulations were restricted to chromosome 1 (N=191,128 SNPs) for computational efficiency, assuming it was representative of the whole genome. These simulations allow us to evaluate the performance of the method in scenarios that approximate realistic GWAS conditions, including correlated SNPs according to true linkage disequilibrium (LD) patterns.
Table 29 displays the number of SNPs rejected and the False Discovery Proportion (FDP), or the proportion of rejected SNPs not in LD with a causal SNP. The cmfdr performs reasonably well across enrichment settings for more highly polygenic phenotypes, rejected SNPs conservatively for 1=0:05, but becoming progressively worse at controlling the FDP for phenotypes with low 1. In comparison, the fdr of Efron (2007) is much more conservative over the entire range of 1, but also has less power. The 2 mixture model of Lewinger et al. (2007) is performs similarly to that of cmfdr, but does not control fdr throughout the range of 1 considered. In particular, their model is very unstable for null GWAS, and performs poorly in the presence of population stratification; if no genomic control (GC) is applied (Devlin and Roeder, 1999), the Lewinger et al. (2007) method rejects far too many SNPs. If standard GC is applied, their method becomes overly conservative, as seen in the real data analysis below.
Median number of SNPs rejected (Rejected) and False Discovery Proportion (FDP) for the proposed cmfdr methodology. Settings include level of covariate enrichment (Enrich.), level of population statification (Strat.), and level of polygenicity (π1). Numbers in brackets give middle 95% of distributions across 100 simulations for each setting.
Real Data Application
The data consist of n=942,772 SNP summary test statistics (SNP z-scores) from a GWAS meta-analysis of eight sub-studies of Crohn's Disease (CD) on a total of 51,109 subjects, obtained through a publicly accessible database Franke et al. (2010). CD is a type of inflammatory bowel disease that is caused by multiple factors in genetically susceptible individuals. For this example the five SNP annotations from Schork et al. (2013) displayed in
Plots of posterior draws showed convergence to stable posterior distributions for all parameters.
The proposed cmfdr methodology rejected far more SNPs than fdr (Efron, 2007). For example, for a 0.05 cut-off, cmfdr rejects 2,742 SNPs whereas fdr rejects only 592. The Lewinger et al. (2007) method rejected 782 SNPs with the same cut-off. The lower number of rejected SNPs compared to cmfdr is due in part to the combination of GC and the lack of empirical null option with their methodology (Lewinger et al., 2007).
The 2,742 SNPS consisted of 108 independent loci (leading SNP cmfdr≦0:05 and more than 1 Mb apart from each other). Of these 108 independent loci, 66 had been previously described in Franke et al. (2010). Franke et al. (2010) described an additional 5 loci that were not discovered using a 0:05 cut-off; however, in this analysis, each of these loci had a cmfdr<0:06. 42 novel loci where the leading SNP had a cmfdr≦0:05. To demonstrate that the method identifies candidate SNPs pleiotropy analysis was performed. Given that Crohn's disease is known to share etiology, including pleiotropic genetic factors (Cho and Brant, 2011) with Ulcerative Colitis (UC), it is likely that causal SNPs would show joint associations. Significant enrichment was found for nomial associations (p<0:05) with UC (Anderson et al., 2011) for both the 71 previously discovered loci (bonferroni adjusted hypergeometric p-value=1.33×10−36) and the 42 novel loci (bonferroni adjusted hypergeometric p-value=6.24×10−5).
Power to detect non-null SNPs using cmfdr vs. usual fdr is displayed in
Further analyses was performed on CD substudies to determine whether this observed increase in power translates to increased replication rates in de novo samples. The CD meta-analysis was composed of summary statistics from eight substudies (Franke et al., 2010). Z-scores were computed from each of the 70 possible combinations of four substudies, leaving the z-scores computed from the remaining four independent substudies as test samples. Fdr and cmfdr were then estimated for each training sample. For a given fdr cut-off, the number of SNPs that replicated in the test sample was determined. Replication was defined as p≦0:05 and with the same sign as the corresponding z score in the training sample.
Number of replicated SNPs was much higher using cmfdr compared to fdr. For example, for usual fdr there was an average of 192 replicated SNPs (44% of SNPs declared significant) with an fdr cut-off of 0:05 in the training sample. In contrast, with the same cut-off using cmfdr there was an average of 1,068 SNPs (47% of declared significant SNPs) that replicated according to this definition, or almost 5.6 times as many SNPs. Similar increases in number of replicated SNPs were observed for other cutoffs in the range. Note, replication rates (44% and 47%) were much lower than the nominal fdr level of 0:05 would suggest. This is due to a significant degree of heterogeneity in the substudies (Franke et al., 2010), as well as limited sample sizes. For comparison, the usual GWAS threshold of 5×10=8 resulted in an average of 89 replicated SNPs, comprising 54% of declared significant SNPs from the training samples. In general, fdr provides a conservative estimate of the non-replication rate in an infinitely sized replication sample from a population like that of the training sample. Application of the cmfdr methodology in other GWAS samples with more homogeneous training and test sets has lead to replication rates much closer to nominal levels while maintaining large advantages in number of replicated SNPs over usual fdr.
Anderson, C. A. and Boucher, G. and Lees, C. W. and et al. (2011). Meta-analysis identifies 29 additional ulcerative colitis risk loci, increasing the number of confirmed associations to 47. Nature Genetics, 43, 246-252.
Participant Samples
Summary statistics from a large MS GWAS study performed by IMSGC (15), n=27 148, and from two large GWAS studies from the Psychiatric GWAS Consortium (PGC), PGC Schizophrenia sample (7), n=21 856, PGC Bipolar disorder sample (12), n=16 731. P-values and minor allele frequencies from the discovery samples were included in the analyses. For follow up analysis, the PGC Major depressive disorder (MDD)(25), Autism Spectrum Disorder (AUT)(26) and Attention Deficit/Hyperactivity Disorder (ADHD) (27) GWAS summary statistics were utilized.
Statistical Analyses
Conditional Q-Q Plots for Pleiotropic Enrichment
To visually assess pleiotropic enrichment, Q-Q plots conditioned on ‘pleiotropic’ effects (13, 23) (
Conditional Replication Rate
For each of the 17 sub-studies contributing to the final meta-analysis in SCZ, the z-scores were independently adjusted using intergenic inflation control (29). 1000 combinations of eight and nine sub-study groupings were randomly sampled. The eight-or-nine-study combined discovery zscore and eight-or-nine-study combined replication z-score was calculated for each SNP as the average z-score across the sub-studies multiplied by the square root of the number of studies. For discovery samples the zscores were converted to two-tailed p-values, while replication samples were converted to one-tailed pvalues preserving the direction of effect in the discovery sample. For each of the 1000 discovery replication pairs cumulative rates of replication were computed over 1000 equally-spaced bins spanning the range of −log 10(p-values) observed in the discovery samples. The cumulative replication rate for any bin was the proportion of SNPs with a −log 10 (discovery p-value) greater than the lower bound of the bin with a replication p-value<0.05 and the same sign as the discovery sample. Cumulative replication rates were calculated independently for each of the four pleiotropic enrichment categories. For each category, the cumulative replication rate for each bin was averaged across the 1000 discovery-replication pairs and the results are reported in
Conditional Replication Effect Size
Using the same z-score adjustment scheme and sampling method used for estimating cumulative replication rates (see above), the relationship of replication effect size of the discovery sample versus replication samples (
Improving Discovery of SNPs in SCZ and BD Using Conditional FDR
To improve detection of SNPs associated with SCZ and BD, a genetic epidemiology approach was employed, leveraging the MS phenotype from an independent GWAS using conditional FDR as outlined in Andreassen (13, 23). Specifically, conditional FDR is defined as the posterior probability that a given SNP is null for the first phenotype given that the p-values for both phenotypes are as small as or smaller than their observed p-values. A conditional FDR value for each SNP in SCZ given the p-value in MS (denoted as FDRSCZ|MS). The same procedure was applied to compute FDRBD|MS for each SNP. To display the localization of the genetic markers associated with SCZ and BD given the MS effect, a ‘Conditional Manhattan plot’, plotting all SNPs within an LD block in relation to their chromosomal location was used. As illustrated for SCZ in
Annotation of Novel Loci
Based on 1KGP linkage disequilibrium (LD) structure, significant SNPs identified by conditional FDR were clustered into LD blocks at the LD−r2>0.2 level. These blocks are numbered (locus #) in Tables 31 and 32. Any block may contain more than one SNP. Genes close to each SNP were obtained from the NCBI gene database. Only blocks that did not contain previously reported SNPs or genes related to previously reported SNPs were deemed as novel loci in the current study (Tables 31 and 32). Loci that contained either SNPs or genes known to be associated with SCZ were considered as replication findings.
HLA Allele Analysis
The PGC1 genotype data from the 17 sub-studies were used for HLA imputation (a detailed description of the datasets, quality control procedures, imputation methods, and, principal components estimation, are given in reference 7). First, genotypes of SNPs in the extended MHC (Major Histocompatibility Complex) (chr6: 25652429-33368333) of each individual in all the samples were extracted. Then, the program HIBAG30 was used to impute genotypes of classical HLA alleles for each sample separately, using the parameters trained on the Scottish 1958 birth cohort data. HLA alleles with posterior probabilities≧0.5 and frequency>0.01 were used in subsequent analysis. The genotypes of the 63 HLA alleles meeting these criteria were encoded as binary variables for the following conditional analysis.
Samples with imputed HLA genotypes were combined before the analysis. First, the logistic regression method implemented in PLINK31 was employed to test HLA alleles for associations with SCZ, using the first 5 principal components and sample indicator variable as covariates. After Bonferroni correction, 5 alleles passed the genomic significance threshold (7.9×10-4). The dosages of SNPs in the MHC, imputed based on HapMap3 data, were tested using logistic regression. The analysis was first performed with only sample indicator variables and the first 5 principal components as covariates and then including, in turn, one of the significant HLA alleles from the previous step as an additional covariate. In addition to the SCZ associated HLA alleles, 4 other alleles reported to be associated with MS were also tested in this framework. A large increase in a SNP's association p-value upon conditioning on HLA alleles is considered to indicate overlap with that HLA allele (Supplementary
Conditional Q-Q Plots
Q-Q plots compare a nominal probability distribution against an empirical distribution. In the presence of all null relationships, nominal p-values form a straight line on a Q-Q plot when plotted against the empirical distribution. For each phenotype, for all SNPs and for each categorical subset (strata), −log 10 nominal p-values were plotted against −log 10 empirical p-values (conditional Q-Q plots). Leftward deflections of the observed distribution from the projected null line reflect increased tail probabilities in the distribution of test statistics (z-scores) and consequently an over-abundance of low p-values compared to that expected by chance, also named ‘enrichment’.
Conditional True Discovery Rate (TDR)
The ‘enrichment’ seen in the conditional Q-Q plots can be directly interpreted in terms of true discovery rate (TDR=1−FDR). More specifically, for a given p-value cutoff, the FDR is defined as
FDR(p)=π0F0(p)/F(p), [1]
where π0 is the proportion of null SNPs, F0 is the null cumulative distribution function (cdf), and F is the cdf of all SNPs, both null and non-null7. Under the null hypothesis, F0 is the cdf of the uniform distribution on the unit interval [0,1], so that Eq. [1] reduces to
FDR(p)=π0p/F(p), [2].
The cdf F can be estimated by the empirical cdf q=Np/N, where Np is the number of SNPs with pvalues less than or equal to p, and N is the total number of SNPs. Replacing F by q in Eq. [2],
Estimated FDR(p)=π0p/q, [3],
which is biased upwards as an estimate of the FDR32. Replacing π0 in Equation [3] with unity gives an estimated FDR that is further biased upward;
q=p/q [4].
If π0 is close to one, as is likely true for most GWASs, the increase in bias from Eq. [3] is minimal. The quantity 1−p/q, is therefore biased downward, and hence a conservative estimate of the TDR. Referring to the Q-Q plots, q* is equivalent to the nominal p-value divided by the empirical quantile, as defined earlier. The FDR estimate is ready directly off the Q-Q plot as
−log 10(q*)=log10(q)−log10(p), [5]
e.g., the horizontal shift of the curves in the Q-Q plots from the expected line x=y, with a larger shift corresponding to a smaller FDR. This is illustrated in
Further Analyses Performed
Significance of Conditional Enrichment
After pruning the SNPs by removing SNPs in linkage disequilibrium (r2≧0.2), 95% confidence intervals were calculated for the conditional Q-Q plots. From these confidence intervals standard errors were calculated and two sample t-tests were used to estimate the difference (degree of departure) of the empirical distribution of SNPs in SCZ or BD (phenotype 1) that are above a given association threshold (−log 10(p)≧1, −log 10(p)≧2, −log 10(p)≧3, −log 10(p)≧4; red lines) in MS (phenotype 2) compared to the −log 10(p)≧0 in phenotype 1 category (blue line). The same procedure was used for the “censored data” of MS conditional on SCZ.
Conditional Analysis of HLA Alleles
It was tested if the associated HLA signals were independent of each other by conditional analysis between them. Samples with imputed HLA allele genotypes were combined before the analysis. The logistic regression method implemented in PLINK8 was employed to test each significant HLA allele for associations with SCZ, including another significant HLA allele, the first 5 principal components and sample indicator variable as covariates. It is more probable that the observed associations were driven by a single haplotype-block, consisting of the 5 individual HLA alleles.
The Effect of HLA Region on Enrichment
The enrichment method was reapplied to the same dataset with SNPs either located within the HLA region or in LD (r2>0.2) with such SNPs (in total 9379 SNPs). These results indicate that the enrichment of SCZ conditional on MS is largely the consequence of the HLA region (Supplementary
Enrichment of SCZ SNPs Due to Association with MS—Conditional Q-Q Plots
Conditional Q-Q plots for SCZ given level of association with MS (
Association with MS Increases Conditional True Discovery Rate (TDR) in SCZ
Variation in enrichment in pleiotropic SNPs is associated with corresponding variation in conditional TDR, equivalent to one minus the conditional FDR (28). A conservative estimate of the conditional TDR for each nominal p-value is equivalent to 1−(p/q) as plotted on the conditional Q-Q plots (see Methods). This relationship is shown for SCZ conditioned on MS in a conditional TDR plot (
Replication Rate in SCZ is Increased by MS Association
To address the possibility that the observed pattern of differential enrichment results from spurious (e.g., non-generalizable) associations due to category-specific stratification or statistical modeling errors, the empirical replication rate was examined across independent sub-studies for SCZ.
Replication Effect Size Depends Upon MS Association
Consistent with the pattern observed for replication rates in SCZ sub-studies (see above), it was found that the effect sizes of SNPs in enriched categories (e.g. −log 10 (pMS)≧3) replicated better than effect sizes of SNPs in less enriched categories (e.g. −log 10(pMS)≧0;
SCZ Gene Loci Identified with Conditional FDR
Conditional FDR methods (13, 23) improve the ability to detect SNPs associated with SCZ due to the additional power generated by use of the MS GWAS data. Using the conditional FDR for each SNP, a ‘conditional FDR Manhattan plot’ for SCZ and MS (
Effect of the Size of Strata on Enrichment
The observed enrichment was further confirmed by performing the same analysis on additional categories (−log 10 (pMS)≧4, −log 10(pMS)≧5 and −log 10(pMS)≧6.
Distribution of Allele Frequencies in Strata
The distribution of the minor allele frequencies (MAF) of the corresponding SNPs of each stratum were identified from the 1KGP.
HLA Imputation and Association Analysis
Among the loci identified by conditional FDR methods, eight are located in the MHC (Table 32). It is possible that these signals may be driven by common HLA alleles affecting both SCZ and MS. To test this hypothesis, HLA class I and class II alleles were investigated using the PGC1 genotype data (see Methods). Association analysis between imputed HLA alleles and SCZ was performed. The alleles HLA-B*08:01, HLA-C*07:01, HLA-DRB1*03:01. HLADQA1*05:01 and HLA-DQB1*02:01 are negatively associated with SCZ (p<7.8×10−4). Among these, HLA-DRB1*03:01 and HLA-DQB1*02:01 have been reported to be positively associated with MS 15. However, no association was seen with SCZ for the strong MS predisposing HLA-DRB1*15:01 and HLA-DRB1*13:03 alleles, nor for the protective HLA-A*02:01 allele. It was further tested whether SNPs in the MHC with conditional FDR<0.05 were independent of the association signal with the classical HLA alleles (see Methods). SNPs rs9379780, rs3857546, rs7746199, rs853676 and rs2844776 are to be independent of the HLA allelic signal (
It was further tested if the associated HLA alleles were independent of each other by conditional analysis between them (see Methods). The results indicate that the observed associations are driven by a single haplotype-block (i.e. ancestral haplotype 8.1), consisting of the 5 individual HLA alleles.
The Effect of MHC SNPs on Enrichment
The effect of MHC-related SNPs (SNPs located within the MHC or SNPs within 1 Mb and in LD (r2>0.2) with such SNPs) on the observed enrichment for SCZ and BD conditional on MS was investigated (see
Enrichment Analysis of Other Psychiatry Disorders
Using the analysis approach described above, genetic enrichment in Major depressive disorder (MDD)25, Autism spectrum disorder (AUT)26 and Attention Deficit/Hyperactivity Disorder (ADHD)27 was analyzed. GWAS summary statistics from the PGC conditioned on MS. In contrast to SCZ, none of these phenotypes demonstrated significant enrichment (
1 Loci identified by GWASs without leveraging genetic pleiotropy structure between phenotypes.
2 Loci identified using conditional FDR method on SCZ with CVD.
3 Loci identified using conditional FDR method on SCZ with BD.
1Loci identified by GWASs without leveraging genetic pleiotropy structure between phenotypes.
2Loci identified using conditional FDR method on SCZ with CVD.
3Loci identified using conditional FDR method on SCZ with BD.
Participant Samples
Complete GWAS results in the form of summary statistics p-values were obtained from public access websites or through collaboration with investigators (Table 33). Details on the inclusion criteria and phenotype characteristics of the different GWAS are described in the original publications 4,25-28. There was some overlap among several of the participants in the CVD risk factor GWAS and the SBP GWAS sample4. The relevant institutional review boards or ethics committees approved the research protocol of the individual GWAS and all participants gave written informed consent. All studies adhered to the principles of the Declaration of Helsinki.
Statistical Analyses
Genomic Control
A control method was applied using only intergenic SNPs to compute the inflation factor, λGC and all test statistics were divided by λGC, as detailed in prior publications21,22.
Conditional Quantile-Quantile (Q-Q) Plots for Pleiotropic Enrichment
Enrichment of statistical association relative to that expected under the global null hypothesis can be visualized through Q-Q plots of nominal p-values obtained from GWAS summary statistics. Genetic enrichment results in a leftward shift in the Q-Q curve, corresponding to a larger fraction of SNPs with nominal −log 10 p-value greater than or equal to a given threshold. Conditional Q-Q plots are constructed by creating subsets of SNPs based on the significance of each SNP's association with a related phenotype, and computing Q-Q plots separately for each level of association (for further details, see references 21, 22). Conditional Q-Q plots of empirical quantiles of nominal −log 10(p) values were constructed for SNP association with SBP for all SNPs, and for subsets of SNPs determined by the nominal p-values of their association with each of the 12 related phenotypes (−log 10(p)≧0, −log 10(p)≧1, −2 log 10(p)≧2, and −log 10(p)≧3 corresponding to p≦1, p≦0.1, p≦0.01, and p≦0.001, respectively). The nominal p-values (−log 10(p)) are plotted on the y-axis, and the empirical quantiles (−log 10(q), where q=1−cdf(p)) are plotted on the x-axis. To assess polygenic effects, the conditional Q-Q plots were focused on SNPs with nominal −log 10(p)<7.3 (corresponding to p>5×10-8).
Conditional False Discovery Rate (FDR)
Enrichment seen in the conditional Q-Q plots can be directly interpreted in terms of False Discovery Rate (FDR)21,22 (equivalent to 1−True Discovery Rate (TDR)35). A conditional FDR method22,36,37 was applied, and TDR plots were constructed, as described earlier21,22.
Conditional Statistics—Test of Association with Systolic Blood Pressure
To improve detection of SNPs associated with SBP, SNPs were conditioned based on p-values in the related phenotype21.22. A conditional FDR value (denoted as FDRSBP|related-phenotype) was assigned for SBP to each SNP, for each related phenotype by interpolation, using a two-dimensional look-up table of conditional FDR values21,22 computed for each of the specific datasets used in the current study. All SNPs with FDRSBP|related-phenotype<0.01 (−log 10(FDRSBP|related-phenotype)>2) in SBP given association with any of the 12 related phenotypes are listed in Table 33 after ‘pruning’ (i.e., removing all SNPs with r2>0.2 based on 1000 Genomes Project linkage disequilibrium (LD) structure). A significance threshold of FDR<0.01 corresponds to 1 false positive per 100 reported associations. To illustrate the localization of the genetic markers associated with SBP given the related phenotype effect, a ‘Conditional FDR Manhattan plot’ was generated, plotting all SNPs within an LD block in relation to their chromosomal locations. The strongest signal in each LD block was identified by ranking all SNPs in increasing order, based on the conditional FDR value for SBP, and then removing SNPs in LD r2>0.2 with any higher ranked SNP. Thus, the selected locus was the most significantly associated with SBP in each LD block.
Results
Pleiotropic Enrichment—Polygenic Overlap.
Conditional Q-Q plots for SBP conditioned on nominal p3 values of association with LDL, BMI, BMD, TID, SCZ, and CeD showed enrichment across different levels of significance (
Gene Loci Associated with SBP.
A “conditional FDR” Manhattan plot showed the 62 independent gene loci significantly associated with SBP based on conditional FDR<0.01 obtained from associated phenotypes. The 30 complex loci and 32 single gene loci (after pruning) were located on 16 chromosomes (Table 34). Only 11 of these loci would have been discovered using standard statistical methods (Bonferroni correction; bold values in the “SBP p-value” column, Table 34). Using the FDR method, 25 loci were identified (bold values in the “SBP-FDR” column, Table 34). The remaining 37 loci would not have been identified in the current sample without using the pleiotropy informed conditional FDR method. Of the 62 loci identified, 42 were novel; 20 were reported in the primary analysis of the current sample4. Many of these new loci are located in regions with borderline significant association with SBP in previous studies4. Of interest, several loci had multiple pleiotropic SNPs from several associated phenotypes, indicating overlapping genetic factors among these phenotypes. Follow-up Ingenuity Pathways Analysis (IPA) identifying the traits in the categories “Cardiovascular disease” or “Cardiovascular System Development and Function”, respectively, that may be affected by the gene heterogeneities in the vicinity of the indicated SBP associated genes were identified. A large proportion of SBP associated genes are functionally related.
indicates data missing or illegible when filed
9.86E−11
0.00003
2.38E−10
0.00003
0.00665
4.56E−08
0.00112
0.00909
1.08E−09
0.00004
0.00271
0.00665
0.00162
0.00558
0.00909
0.00386
4.07E−08
0.00093
0.00271
3.40E−08
0.00077
0.00271
5.48E−09
0.00015
2.77E−10
0.00003
2.43E−10
0.00003
0.00135
1.25E−09
0.00005
2.98E−08
0.00077
5.35E−10
0.00003
5.14E−13
0.00003
0.00783
4.58E−10
0.00003
1.63E−10
0.00003
0.00323
0.00162
0.00783
0.00558
0.00112
0.00665
0.00135
All publications and patents mentioned in the above specification are herein incorporated by reference. Various modifications and variations of the described method and system of the invention will be apparent to those skilled in the art without departing from the scope and spirit of the invention. Although the invention has been described in connection with specific preferred embodiments, it should be understood that the invention as claimed should not be unduly limited to such specific embodiments. Indeed, various modifications of the described modes for carrying out the invention that are obvious to those skilled in the medical sciences are intended to be within the scope of the following claims.
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/US2014/011014 | 1/10/2014 | WO | 00 |
| Number | Date | Country | |
|---|---|---|---|
| 61751420 | Jan 2013 | US |
