Alzheimer’s Disease Biological Domains
The primary goal of defining a structured set of biological domains is to standardize areas of disease-associated biology that can serve as a common reference point for the analysis of large data sets. The CADRO ontology, discussed in the methods section, provided a useful starting place for the development of the biological domains, having been designed for a related purpose, and regularly employed to categorize clinical investigations6,109. Leveraging CADRO as a starting place and augmenting this ontology based on literature curation, we established a set of nineteen biological domains that cover the AD endophenotypic space to objectify the processes involved in AD pathogenesis (Supp Table 1). In total we use 7,127 unique GO terms (16.4% of all terms in the GO) to annotate the biological domains, and the number of GO terms annotated to each biological domain varies enormously (Fig. 1A, Supp Table 2). The ‘Synapse’ domain requires the largest number of GO terms to define, with 1,379 constitutive GO terms, while ‘Tau Homeostasis’ is the smallest with only 10 GO terms. ‘APP Metabolism’ is the second smallest biological domain; the two smallest biological domains focus on gene-centric processes, requiring fewer terms to annotate than larger domains with broader biological focus, such as ‘Lipid Metabolism’ or ‘Proteostasis’. Each biological domain was designed to be discrete from the others. This is reflected in the sparse overlap of shared GO terms between domains (Fig. 1A). The few GO terms that are present in more than one domain are truly inextricable (e.g. the term “mitophagy” defined as “The selective autophagy process in which a mitochondrion is degraded by macroautophagy'' legitimately resides in both the Mitochondrial Metabolism and Autophagy domains), in which case the repetition was allowed, as it represents a meaningful intersection of biological areas.
The genes annotated to the biological domains reveal a different pattern, with high levels of promiscuity between the domains (Fig. 1B). The number of genes annotated to each domain is roughly proportional to the number of GO terms per domain, with almost two orders of magnitude separating the largest domain ‘Proteostasis’ and the smallest ‘Tau Homeostasis’. Many genes are annotated to multiple GO terms subordinate to multiple biological domains. This may represent a convergence of biologically related processes, or the pleiotropy associated with a given gene’s function. While a plurality of annotated genes (30%) is only annotated to a single biological domain, many participate in multiple domains each (Fig. 1C), and genes annotated to more biological domains tend to have higher overall Target Risk Scores (see below) (Pearson r = 0.158, p = 2.8x10− 80).
Target Risk Score Overview
The TREAT-AD Target Risk Score (TRS) is a metric designed to rank potential disease involvement of specific genes based on multiple independent lines of evidence to objectify the prioritization of potential targets and disease-linked biology. The TRS has two component elements: genetic risk and multi-omic risk, each derived from a meta-analysis harmonizing multiple datasets. In the sections below, we discuss the results for each component, with genetic risk weighted more heavily, receiving up to 3 points, while multi-omic risk has a maximum of 2 points. The rationale for providing more weight to the genetics score component reflects the greater success in clinical trials for targets with genetic support110–112.
Genetic Risk Score Component
The genetic risk score component is a summary of genetic evidence supporting the target gene's association with late-onset Alzheimer's disease from multiple genetic studies. The score is based on genetic evidence retrieved from both genome wide association (GWAS) and GWAS by proxy (GWAX) along with QTL studies. In total, 27 different Alzheimer’s association data sources were queried (Supp Table 3). These datasets are not fully independent, as the patient cohorts used partially overlap between studies, and as such many genes are scored with associations across multiple studies (Supp Fig. 1A). Given the differences between the populations sampled and methodologies applied, we chose a windowed approach to assign variants to genes rather than accounting for the various linkage disequilibrium blocks in the source data. Furthermore, as the goal is to score as many targets as possible for genetic association with AD, we selected a permissive significance threshold for variant inclusion in order to stratify targets that are nominally significant in association with AD but do not meet a more standard criteria for genome-wide significance (e.g. p < 5x10− 8). Had we restricted our analyses to genes that meet this threshold, only 1,706 of the 60,664 assessed genes (2.81%) would be included (Supp Fig. 1B). The vast majority of the genome would be omitted from consideration, using this standard, including many genes that have notable association with AD, but have not yet reached genome-wide significance for various reasons, such as a lack of statistical power and/or the limited diversity in cohort ancestry. Performing gene set enrichment analyses (GSEA) using the combined GWAS score to rank genes (Supp Fig. 1C) enriches GO terms from the Immune Response, Synapse, and Lipid Metabolism biological domains (Supp Fig. 1D).
We analyzed genes with identified eQTL or pQTL from studies where genetic variation is associated with expression differences in postmortem AD brains. In contrast with the GWA studies, only 492 genes are identified with a significant QTL in all three study populations used (Supp Fig. 1E-G), which is consistent with previous reports of limited overlap between brain eQTL and pQTL58. Accordingly, GSEA using the integrated QTL signal across the studies to rank genes enriches GO terms from only 5 biological domains - namely Cell Cycle, Proteostasis, Oxidative Stress, Mitochondrial Metabolism, and Immune Response (Supp Fig. 1H).
In addition to accounting for the association of a gene with AD traits, we also sought to assess the predicted severity of identified variants. We collected 7,198,607 unique variants with at least a nominal association to AD or related phenotypes, including 50,833 variants predicted to affect the coding sequence of a target gene (Supp Fig. 2A) and 1,365,759 variants associated with altered gene or protein expression (eQTL or pQTL) in postmortem AD brain (Supp Fig. 2D). Variant severity is scored differently for coding and noncoding variants; while the preponderance of identified variants are noncoding (99.29%), the prediction of functional effects of variation is more straightforward for coding variants. Coding variant severity was assessed for 17,892 distinct genes, 13,088 of which (73.2%) are identified with at least one deleterious coding variant. In general, genes with larger CDS lengths tend to have a higher number of deleterious coding variants (Supp Fig. 2B). In addition, there are 939 genes that are among the least tolerant of loss-of-function variation (bottom 10% gnomAD LOEUF) where we identify at least one deleterious coding variant associated with AD. The coding variant summary is enriched for genes in Structural Stabilization, Proteostasis, Immune Response, Endolysosome, Synapse, Epigenetic, and Mitochondrial Metabolism biological domains (Supp Fig. 2C). Noncoding variant severity was assessed for variants identified through QTL studies (namely, those associated with altered expression of the target gene or protein) for a total of 15,812 distinct genes. The noncoding variant summary score integrates information from the RegulomeDB probability score and DeepSEA mean -log e-value (Supp Fig. 2E). The noncoding variant summary is enriched for terms from the Immune Response and Proteostasis biological domains (Supp Fig. 2F). For the 7,394 genes assessed for both coding and noncoding variant severity, there is a very weak positive correlation between the coding and noncoding variant summary scores (Pearson r = 0.05, p = 1.19x10− 5, Supp Fig. 2G).
The score incorporates phenotypic evidence supporting a given target’s potential to impact AD relevant phenotypes from both human and model organism sources. The human phenotype score leverages the Human Phenotype Ontology to query the phenotypic abnormalities annotated to each gene. 4,131 genes have at least one phenotype in common with AD (MONDO:0004975) or dementia (MONDO:0001627) (Supp Fig. 3A), and these are enriched for terms annotated to 15 biological domains, but predominantly to the Synapse, Mitochondrial Metabolism, Lipid Metabolism, and Proteostasis biological domains (Supp Fig. 3B). For relevant ortholog phenotypes we queried the uPheno ontology97 and scored 8,901 genes with orthologs that have at least one phenotype in common with AD or dementia (Supp Fig. 3C). Genes with a strong ortholog phenotype score are enriched in GO terms from 17 biological domains, with Synapse, Immune Response, and Lipid Metabolism being the three with the largest number of enriched GO terms (Supp Fig. 3D). For the 3,063 genes with relevant phenotypes detected in both human and model organism resources there is a relatively weak positive correlation in phenotype scores (pearson r = 0.20, p = 6.7x10− 30, Supp Fig. 3E). We provide an additional weight to genes that have been positively identified within MODEL-AD studies to be associated with LOAD113. MODEL-AD is an NIA funded effort to develop new animal models that better represent features of human AD pathophysiology based on signatures from AD GWAS. We include in our score whether genes have a model in development by MODEL-AD as well as the maximum correlation between human transcriptomic module gene expression from AMP-AD and mouse model gene expression, where available.
The genetic risk score was calculated for a total of 60,664 genes. After tallying all genetic evidence sources, we found support for 52,375 genes, though the majority have only weak or sparse evidence connecting the locus to Alzheimer’s risk. To limit spurious associations, the list was reduced to only the top 25 percent of ranked genes (15,166 genes) or those scored in the multi-omics analysis (see below). This resulted in a final list of 24,278 scored genes (Supp Table 4). The genes re-introduced to genetic scoring based on their inclusion in the multi-omic analyses are represented as the lower mode of the distribution (Fig. 2A, between 0–1). Genes contained within known AD GWAS loci are enriched among the top scores (Fig. 2A). Genes identified as high in AD risk by Open Targets (https://www.opentargets.org/)/)105–107,114, a large-scale effort to rank genes based on genetic support for translational relevance, were also among the top scoring genes (Fig. 2A). Comparing the TREAT-AD genetics score with the Open Targets genetic association score for AD (Fig. 2B) reveals a very weak positive correlation (Pearson r = 0.185, p = 4.4x10− 4) which is stronger when only considering known AD GWAS genes (Pearson r = 0.321, p = 6.7x10− 3). In general, most genes receive a relatively higher score from the TREAT-AD genetics score. GSEA using the genetics score enriches GO terms from 17 of 19 biological domains; the biological domains with the largest number of enriched GO terms are Synapse, Lipid Metabolism, and Structural Stabilization (Fig. 2C-D, Supp Table 5). The Open Targets genetic association score enriches GO terms from 10 biological domains, with terms in the APP Metabolism domain by far the most significantly enriched and Synapse, Immune Response, and Lipid Metabolism being the domains with most terms enriched (Supp Fig. 4C). The relative emphasis of APP Metabolism in the Open Targets score likely reflects the inclusion of evidence from early onset autosomal dominant forms of AD in which causal variants are clustered in the Presenilin genes (PSEN1, PSEN2) involved in proteolytic processing of APP, whereas the TREAT-AD genetics score draws primarily from genetic associations of late onset sporadic forms of the disease. Notably Synapse, Lipid Metabolism, Structural Stabilization, and Immune Response are among the biological domains with the largest number of enriched GO terms for both scores.
Multi-omic Risk Score Component
The multi-omic risk score is a summary metric encapsulating available evidence supporting that target gene expression is altered in the brains of AD patients. The score makes use of proteomic and transcriptomic datasets generated as part of the consortium efforts of AMP-AD (Supp Tables 6 and 7). For each data modality (i.e. transcriptomic or proteomic), a meta-analysis of samples is used to generate weights for significantly differentially expressed genes based on the observed fold changes (Fig. 3A & 3B). Using ratio of means meta-analysis with a random effects model for each data modality, we analyzed the directionality of expression change due to AD using GSEA (Fig. 3C & 3D). This analysis demonstrates that Synapse and Mitochondrial Metabolism are the biological domains with the largest number of down-regulated GO terms and Immune Response and Structural Stabilization are the biological domains with the largest number of up-regulated GO terms – across both transcriptomics (Fig. 3C) and proteomics (Fig. 3D).
The calculated weights for each modality are combined using a scoring harness (Supp Table 8) that yields a higher score for targets with (A) evidence of differential expression at both the protein level and the RNA level, followed by (B) those targets that are only significantly differentially expressed at the protein level, followed lastly by (C) those that are only significantly differentially expressed at the RNA level (Fig. 3E). The rationale for this harness is two-fold: first, concordant evidence from multiple data modalities leads to higher confidence that a target gene’s expression is altered in AD brains, and second, that protein levels more accurately relate to the state of the biology within an in vivo system35,115,116. Importantly the effect of the harness is tuned to prioritize this imposed hierarchy without disregarding results specific to a single profiling modality. This reflects a balance between higher relevance of proteomic evidence to disease state versus the increased sensitivity of detection for transcriptomics.
We compare the distribution of all scored targets (Supp Table 9) to those scored by the Open Targets Platform105–107, 114 and those nominated for follow-up by the AMP-AD consortium59,117. The 605 targets nominated by AMP-AD investigators tend to have higher multi-omics scores than the population of all targets (Fig. 3F), and there is a very weak correlation between the multi-omics score and the number of nominations received by a given target (r = 0.153; Supp Fig. 5E). Numerous targets receive only one nomination yet are ranked among the highest multi-omics scores. Conversely, there are many targets that received several nominations with a low multi-omics score. This likely reflects the fact that AMP-AD investigators use diverse methods to identify targets, beyond differential expression analysis, and that some modalities (e.g. metabolomics) are not included within the multi-omic score. Future efforts will include work to integrate additional data modalities into the score. Using the multi-omics score to rank genes to perform GSEA results in the enrichment of GO terms from 18 biological domains, with DNA Repair the only domain with no terms enriched. The biological domains with the largest number of enriched terms include Synapse, Immune Response, Structural Stabilization, and Mitochondrial Metabolism (Fig. 3G-H, Supp Table 10). Terms from the Mitochondrial Metabolism domain are both among the most significantly enriched (Fig. 3H) and have the highest normalized enrichment score (NES) (Fig. 3G). Interestingly, the top terms within Mitochondrial Metabolism focus upon mitochondrial translation and complex I of the electron transport chain (Fig. 3G), showing these terms are the most significantly down-regulated biological processes associated with LOAD.
Composite Target Risk Score
The TRS is a metric derived by summing the component risk scores (i.e. genetic and multi-omic scores). The maximum observed score for any target is 4.74 out of a total of 5 (Supp Table 11). As with the component scores, the top TRS scores are enriched for GWAS loci, AMP-AD nominated targets, and targets considered by the Open Targets platform (Fig. 4A). Considering both the genetic and multi-omic scores for each target (Fig. 4B), the targets with the top TRS tend to have relatively higher genetic scores. When we compare the TREAT-AD overall TRS with the Open Targets target score (Supp Fig. 4A) we see that many targets receive a relatively higher score from the TRS, likely due to the unique inclusion of disease-relevant transcriptomic and proteomic datasets from the AMP-AD consortium in the TRS. AD GWAS loci are generally scored highly by both the TREAT-AD TRS and the OpenTargets target score metrics. GSEA using the overall TRS enriches 3,142 GO terms, including 1,358 (43.5%) annotated to at least one of each of the 19 biological domains (Supp Table 12). The biological domains with the largest number of enriched GO terms are Synapse, Immune Response, and Lipid Metabolism (Fig. 4D). In comparison, the Open Targets target score enriches terms from 16 of 19 biological domains, and the domains with the largest number of enriched terms are also Synapse, Lipid Metabolism, and Immune Response (Supp Fig. 4B).
Considering all GO terms within any of the biological domains that are significantly enriched with the TRS, some have more significant enrichment from GSEA using the genetics score versus the omics score to rank genes (Fig. 5A, Supp Tables 5, 10, and 12). For example, from the Lipid Metabolism biological domain the “fatty acid metabolic process” term is enriched using both scores, but for the genetics score based tests the adjusted p-value for the term is 5.9x10− 6 while for the omics score-based test the adjusted p-value is 7.9x10− 3. The inverse is true for the “membrane raft” term which is more significant using the multi-omics scores (adjusted p = 1.3x10− 5) compared to the genetics scores (adjusted p = 1.8x10− 2). A similar dichotomy is observed for terms associated with developmental processes from the Immune Response domain, especially “myeloid leukocyte differentiation”, having relatively more significant enrichment from the genetics score whereas corresponding terms for mature structures and associated processes (e.g. “leukocyte degranulation”) have a relatively more significant enrichment from the multi-omics score (Fig. 5B). Other biological domain resident dichotomies are more consistent across risk modalities. For example, considering all presynaptic and postsynaptic terms from the Synapse biological domain, there are more postsynaptic terms enriched (40 terms) relative to presynaptic terms (9 terms) and the postsynaptic terms are enriched with more significance across the TRS as well as the genetic and multi-omic score components (Fig. 5C). The optimal targets reflect a parity of ranking between the two contributing modalities, with accumulation of both genetic and multi-omic risk. For example, enrichment of genetic risk without multi-omic risk can indicate risk in developmental processes that are not present in later stages or disease or postmortem tissue, whereas enrichment of multi-omic risk without genetic risk can indicate processes that are changing as a response to disease pathology but are not causal to underlying disease etiology. Moreover, processes that are risk enriched across two distinct measures increases our confidence of disease association. For example, while the multi-omics score enriches terms related to mitochondrial electron transport chain complexes I and III, the genetics score only enriches terms related to complex I (Fig. 5D), yet both measures point to the centrality of electron transport chain related events in our scored AD risk. GO terms close to the diagonal dashed line (Fig. 5A), as seen in Structural Stabilization with the term “cell-substrate adhesion”, suggest that the risk associated with genes in these terms are embedded equivalently in each data modality and therefore are worth consideration for further resource development. The GO terms that are significantly enriched but do not associate with any biological domain (i.e. Figure 5A, ‘none’) seem to reflect either biological process categorization that is too general to be mapped into endophenotypic space, such as “ATP Metabolic Process”, or terms that map to high order positions, or do not imply any specific associated biological process, in the molecular function or cell component ontologies within GO, such as “Cell Leading Edge” (Fig. 5A) also rendering them uncategorizable in endophenotypic space.
Example Use Case: AMP-AD Co-expression Modules Risk Scores and Biological Domains
To demonstrate how this information can be useful more broadly, we apply these approaches to published co-expression modules for transcripts20,22and proteins35 produced by AMP-AD consortia members. For each module we calculated the median TRS and genetics and multi-omics component scores (Supp Table 13), performed GO term enrichment analyses using the identities of the genes in each module, and categorized enriched GO terms into biological domains. The STGblue module has both the highest median TRS and the highest median multi-omic scores of the 30 transcript co-expression modules from Wan et al, while the DLPFCblue module has the highest median genetics score (Fig. 6A). Among the 45 protein co-expression modules from Johnson et al, the M11 Cell-ECM Interaction module has the highest median TRS and the highest median multi-omics score while M27 Extracellular Matrix has the highest median genetics score (Fig. 6B). The GO term enrichments for the three proteomics modules with the top median TRS are primarily from the Structural Stabilization (Fig. 6C & 6D) and Mitochondrial Metabolism (Fig. 6E) biological domains. Comparing the GO term enrichments for the modules from each study with the highest median genetics scores and highest median multi-omics scores, we note several common enrichments across modules from all three studies. The GO terms enriched for proteins in modules with the highest median multi-omics (M11, Supp Fig. 6B) and highest median genetics scores (M27, Supp Fig. 6C) are predominantly from the Structural Stabilization biological domain. This is also similar for the genes in the transcriptomic sub-modules from Milind et al, where the enriched GO terms from the sub-module with the highest median multi-omics score (IFGturquoise_2, Supp Fig. 6F) are primarily from the Structural Stabilization biological domain. The transcriptomic sub-module with the highest median genetics score (PHGturquoise_2, Supp Fig. 6E) is enriched for many terms from the Immune Response biological domain. The biological domains enriched for each of the top-scoring transcriptomic modules from Wan et al are very similar - primarily Immune Response, Structural Stabilization, Vasculature, and Lipid Metabolism, in varying orders (Supp Fig. 6H and I). In Wan et al these modules were associated with microglial, pericyte and endothelial cell types relevant to these biological domains. The large number of biological domains enriched is likely because these modules are large (between 504 and 4,673 genes), compared with the proteomic modules (between 28 and 2,231 proteins) and the transcriptomic sub-modules (between 77 and 2,090 genes), so more terms are enriched for each module and the overall pattern of enrichment is less specific. The top risk-enriched co-expression modules from each study implicate terms from the Structural Stabilization and Immune Response biological domains. Moreover, there are GO terms that are shared among these co-expression modules, including “focal adhesion” and “extracellular matrix” for the Structural Stabilization biological domain and “neutrophil degranulation” and “cytokine-mediated signaling pathway” for the Immune Response biological domain.