https://orcid.org/0000-0003-3762-2746
Abstract
Alzheimer’s disease (AD) is a devastating neurological disorder characterized by changes in cell-type proportions and consequently marked alterations of the transcriptome. Here we use a data-driven systems biology meta-analytical approach across three human AD cohorts, encompassing six cortical brain regions, and integrate with multi-scale datasets comprising of DNA methylation, histone acetylation, transcriptome- and genome-wide association studies and quantitative trait loci to further characterize the genetic architecture of AD. We perform co-expression network analysis across more than 1200 human brain samples, identifying robust AD-associated dysregulation of the transcriptome, unaltered in normal human aging. We assess the cell-type specificity of AD gene co-expression changes and estimate cell-type proportion changes in human AD by integrating co-expression modules with single-cell transcriptome data generated from 27 321 nuclei from human postmortem prefrontal cortical tissue. We also show that genetic variants of AD are enriched in a microglial AD-associated module and identify key transcription factors regulating co-expressed modules. Additionally, we validate our results in multiple published human AD gene expression datasets, which can be easily accessed using our online resource (https://swaruplab.bio.uci.edu/consensusAD).
Introduction
Alzheimer’s disease (AD) is a progressive neurodegenerative disease, which is symptomatically characterized by impairment in cognitive and executive functions, including memory loss (1,2). AD is pathologically characterized by the presence of neurofibrillary tangles and senile plaques, which are composed of hyperphosphorylated microtubule associated protein tau (MAPT; tau) inclusions and amyloid beta (Aβ) deposits, respectively (2,3). While aging is the major non-genetic risk factor for AD, several genetic risk factors, including mutations in MAPT, have been found and associated with the disease (4–7). Despite the identification of causal genetic mutations and disease risk genes implicated by genome-wide association studies (GWAS) of AD (6,7), the molecular mechanisms of neurodegeneration are highly complex and are still poorly understood, impeding the design of effective therapeutic intervention.
Transcriptomics coupled with gene network analysis are powerful tools for investigating quantitative molecular phenotypes and pathways underlying disease progression in a genome-wide manner. These approaches have the potential to unravel previously unidentified disease-associated pathways, as well as highlighting disease regulators, which themselves serve as candidate targets for drug development. Within this context, the National Institute of Aging has developed a Target Identification and Preclinical Validation Project of the Accelerating Medicines Project–Alzheimer’s Disease (AMP-AD) consortium (8), whose goal is to integrate high-throughput genomic and molecular data from the diseased brain within a network driven structure (8,9). Under this program, several AD-specific large-scale RNA sequencing (RNA-seq) projects have been conducted. These projects have identified transcriptomic networks and splicing events that are altered in the AD cortex (9–15). However, these studies were limited to specific brain banks, and a holistic meta-analysis across hundreds of samples with integrative genomic approaches is pertinent to robustly identify disease-specific signatures in an unbiased manner.
Here we have taken a functional genomics and integrative systems biology approach to identify human AD-specific transcriptomic alterations conserved across multiple cohorts and multiple cortical regions. We applied consensus weighted gene co-expression network analysis (cWGCNA), a meta-analytical approach, to more than 1200 postmortem brain RNA-seq samples across brain banks and identified robust, disease relevant co-expression modules. With this approach, we were not only able to identify strong, disease-specific signatures across multiple cohorts but also were able to define brain-region specific alterations that are fundamental to disease pathophysiology. We found several neuronal and glial specific modules that are altered in AD but not in normal human aging, and one of the glial modules was enriched in AD candidate risk genes from AD GWAS datasets. We also extensively examined module enrichment of GWAS SNPs associated with other neurological disorders. We used generic and brain-specific ATAC-seq and ChIP-seq datasets to define transcription factors (TFs) regulating the neuronal and glial modules, identifying members of the NF-κB pathway as regulators of AD-associated co-expression changes. In addition, we have leveraged various orthogonal datasets, including DNA methylation, histone acetylation (H3K9ac) and expression and methylation quantitative loci (eQTLs and mQTLs), to provide a comprehensive picture of the genetic alterations associated with AD. While most co-expression modules were enriched in AD-associated epigenetic changes, one glial module was notably and consistently not enriched in these epigenetic marks.
Massive changes in cell composition are concurrent with the progression of AD, impairing accurate estimation of transcriptomic changes from bulk tissue data. To circumvent this problem, we generated single-nuclei transcriptomics data from aged, cognitively healthy postmortem human brain samples and used their transcriptomic signatures to accurately estimate cell-type proportion changes from bulk RNA-seq samples and to classify AD-specific co-expression modules by cell type. Finally, we confirm the robustness of our findings by computing the preservation of co-expression modules in multiple orthogonal published gene expression datasets of AD encompassing several brain regions. Such a comprehensive study across multi-scale datasets has not been performed in neurodegeneration and provides a framework for future meta-analytic integration studies. Together, our findings define a core set of transcriptomic changes altered during the progression of the disease and highlight several regulators of these changes, which may be used as therapeutic targets.
Results
In order to identify robust, key transcriptional changes in AD, we analyzed RNA-seq data from three different studies—Mayo Clinic Brain Bank (Mayo) (10), the Religious Orders Study and Memory and Aging Project (ROSMAP) (14) and Mount Sinai School of Medicine (MSSM) (16). The Mayo dataset contained temporal cortex samples (n = 160) (10), while the ROSMAP dataset was composed of prefrontal cortex samples (n = 632) (14). The MSSM cohort included samples from four different regions—parahippocampal gyrus, inferior frontal gyrus, superior temporal gyrus and the frontal pole (n = 476, Fig. 1) (16). Control samples were classified as Braak stage 0-I. Early-stage pathology samples were defined as Braak stage II-IV and CERAD score of possible AD, while late-stage pathology samples were Braak stage V-VI and CERAD score of probable and definite AD (see Supplementary Material, Table S1 for all case information). All datasets were uniformly processed and regressed for cell-type proportion changes and other technical and biological covariates (see Materials and Methods).
Schematic representation of the various analyses performed in the study. Consensus weighted gene co-expression network analysis (cWGCNA) with datasets from three different studies—Mayo Clinic Brain Bank (Mayo) temporal cortex (TC); Religious Orders Study and Memory and Aging Project (ROSMAP) prefrontal cortex (PFC); Mount Sinai School of Medicine (MSSM) para-hippocampal gyrus (PHG), inferior frontal gyrus (IFG), superior temporal gyrus (STG) and frontal pole (FP)—identified Alzheimer’s disease (AD)-associated modules. Distinct cell-types were detected using single-nuclei RNA-sequencing (snRNA-seq) from postmortem human brain tissue revealing cell-type proportion changes in AD. Disease-associated modules and cell-type clusters were used for functional pathway analysis in addition to GWAS enrichment, transcription factor binding site enrichment and validation in published datasets.
cWGCNA (17,18) is a widely used analysis method to discover biologically relevant gene modules across multiple gene expression datasets. We performed cWGCNA on the three bulk tissue RNA-seq datasets—Mayo, ROSMAP and MSSM—to examine disease-related changes conserved across more than a thousand individuals from different cortical regions and different brain tissue repositories (Fig. 2A). We identified 10 consensus modules significantly correlated with AD diagnosis—four negatively (CM7, CM1, CM10 and CM12) and six positively correlated (CM16, CM5, CM4, CM8, CM9 and CM23; FDR < 0.05; Fig. 2B, Supplementary Material, Table S2). To better understand these modules, we examined the relationships between co-expression modules using a bi-weighted mid-correlation (Fig. 2C). We found modules CM9 and CM23 were strongly positively correlated, but CM23 was strongly anti-correlated with CM1, providing us an idea of how these modules function together. In addition, cWGCNA modules were found to be concordant with gene modules constructed using several other co-expression and gene clustering algorithms (19–22) (see Materials and Methods, Supplementary Material, Fig. S1).
Consensus transcriptomic analysis identifies disease-relevant modules. (A) Consensus WGCNA dendrogram showing consensus modules from Mayo, ROSMAP and MSSM datasets. Red indicates positive correlation of covariate (diagnosis, age, RIN) with gene expression, and blue is negative correlation. (B) Signed correlation of mRNA module eigengenes with AD diagnosis. (C) Multidimensional scaling plot demonstrates relationships between consensus modules significantly correlated with diagnosis. Color signifies biweighted mid-correlation (R2) value. Modules are also annotated with summarizing terms.
Identification of cell-type clusters in the aged human brain
Due to the cell-type composition changes during the progression of the disease, profiling gene expression in only bulk tissue samples may obscure biologically relevant cell-type specific changes. While single-cell RNA-seq allows us to evaluate transcriptional changes within cell-types, it is prohibitively costly to execute on large cohorts (i.e. hundreds of individuals). To circumvent this issue, we developed a framework that leverages single-cell RNA-seq data from a smaller cohort in conjunction with co-expression network analysis in order to estimate cell-type specific transcriptomic changes in large, bulk tissue RNA-seq datasets.
We isolated nuclei and performed single-nuclei RNA-seq (snRNA-seq, n = 27 321 nuclei) on postmortem human brain tissue from aged, neurologically healthy controls (n = 5, 67 to 90+ years old, PFC, Supplementary Material, Table S1) to clarify cell-type proportions and the corresponding transcriptional profiles of the aged, adult human cortex. Through unsupervised Louvain community detection, we identified 18 transcriptionally distinct cell populations that we visualized using t-distributed stochastic neighbor embedding (t-SNE; Fig. 3A, Supplementary Material, Table S3). Canonical cell-type marker genes, such as SLC17A7 for neurons, were examined in order to assign a cell type to each of the 18 clusters (Fig. 3B, Supplementary Material, Fig. S2A). The majority of the profiled nuclei was classified as excitatory neurons (Fig. 3C). This strategy identified both neuronal and non-neuronal subpopulations: seven excitatory neuron, four inhibitory neuron, two astrocyte and two oligodendrocyte subpopulations. We additionally found clusters consistent with endothelial cells, microglia and oligodendrocyte precursor cells (OPCs), and there was one cluster we were unable to confirm its cell-type origin.
Single nuclei RNA-seq on the aged human brain elucidates bulk tissue RNA-seq analyses. (A) Cell-type clusters visualized with t-SNE plot. Cell types were labeled post hoc based on marker gene expression. (B) Marker gene expression in cell-type clusters visualized with same t-SNE plot as in (A). (C) Bar graph indicating fraction of nuclei identified as each cell type. (D–E) Dot plots showing differentially expressed cell-type markers (D) and excitatory neuronal markers (E) in cell-type clusters. Color denotes average gene expression, while size denotes percent of nuclei expressing the gene. (F) Boxplots depicting estimated cell-type abundance with pathological state for four different cell-type clusters—excitatory neuron cluster 1, astrocyte cluster 1, microglia cluster and oligodendrocyte cluster 1 (top to bottom)—in the Mayo TC dataset. (G) Enrichment of consensus modules in cell-type clusters identified by snRNA-seq.
Moreover, differential gene expression analysis identified gene markers specific to each cell-type cluster (Fig. 3D, Supplementary Material, Fig. S2B). These included previously established cell-type markers, such as PVALB, NEUROD6 and CX3CR1, and layer-specific markers like SV2C and HTR2C (23,24). We also utilized the Allen Human Brain Atlas to verify the expression of these marker genes in the human cortex (Supplementary Material, Fig. S2C). Since we found multiple neuronal subpopulations, we aimed to characterize them further, distinguishing three selectively and highly co-expressed genes in each neuronal cluster (Fig. 3E, Supplementary Material, Fig. S2D). Based on these genes we were able to identify our excitatory neuron clusters 7 and 8 as deep layer neurons, while our inhibitory neuron cluster 3 likely originated from layer 3 neurons. In addition, we found inhibitory neuron subpopulations similar to Darmanis et al. (25)—our inhibitory neuron clusters 1 and 4, CALB2-VIP and CRH-PVALB neurons, respectively. We also examined genes differentially expressed across the cell-type clusters in respect to sex, identifying X chromosome inactivation gene XIST upregulated in astrocytes from females (Supplementary Material, Fig. S3, Supplementary Material, Table S4). Additionally, NRXN1 was upregulated in several clusters from males and has been previously associated with male-specific behavioral changes (26,27).
Clarifying AD bulk tissue RNA-seq changes with snRNA-seq
We reasoned that our cell-type clusters from snRNA-seq in the aged human brain could be leveraged to estimate AD-related changes in cell-type proportions across different brain regions and by disease stage. We devised a statistical framework to estimate the cell-type composition in bulk tissue RNA-seq data, not regressed for cell-type proportion changes, using our cell-type clusters (see Materials and Methods). Using the specificity score obtained from single-nuclei data, we independently estimated cell-type composition in each of the bulk RNA-seq datasets used. In the Mayo TC Set, we found that all excitatory and inhibitory neuronal subpopulations significantly decreased with late stage pathology (P < 0.001), as expected; however, there was little change with early stage pathology (P > 0.05; Fig. 3F, Supplementary Material, Fig. S4A). On the other hand, microglia and astrocytes increased with late stage pathology (P < 0.001). Microglia only tended to increase at early stage pathology (P > 0.05), while both astrocyte subpopulations significantly increased (P < 0.05). We also saw that at late stage pathology, OPCs greatly decreased (P < 0.001), but both oligodendrocyte subpopulations tended to increase (P > 0.05), suggesting the differentiation of OPCs into oligodendrocytes with increasing pathology. In addition, we used a published deconvolution algorithm, CIBERSORT (28), and found that our results were comparable to CIBERSORT (Storey’s π1 = 0.93).
Notably, AD progresses in an anatomical manner (1,29), and gene expression changes may be causal or an effect of these region-specific cellular changes. While we found similar patterns in the other datasets, we noticed that RNA-seq datasets from temporal cortical regions demonstrated more pronounced changes in cell-type proportions than frontal cortical regions (Supplementary Material, Fig. S4B, Supplementary Material, Figs S5–S6). This regional difference resembles the progression of pathology in which the temporal cortex is affected before the frontal cortex in AD (1,29).
Accordingly, we also utilized the cell-type clusters from snRNA-seq in order to further characterize the consensus modules we identified through cWGCNA. We classified the consensus modules as neuronal or non-neuronal based on enrichment of cell-type markers and additionally confirmed our classification by integrating our snRNA-seq results (Fig. 3G). For example, upregulated consensus module CM8 was selectively enriched in the microglia cell-type cluster, while CM23 was selectively enriched in both astrocyte clusters.
Neuron-specific co-expression changes in AD
We next wanted to examine the AD-associated modules more closely and used experimentally derived databases of human protein–protein interactions (PPIs) from Inweb (30) and Biogrid (31) to guide our attention to important genes. We created integrated co-expression-PPI networks based on the expectation that a subset of the edges defined by co-expression are due to PPI (32–35); such interactions functionally annotate the networks’ edges and provide orthogonal module validation with independent data (35). While gene expression and protein levels do not always correlate, we used a recently published dataset of proteomics in human AD (n > 2000 samples) to further investigate each cWGCNA module (36). We projected cWGCNA modules into a human AD proteomics dataset and observed agreement between module eigenprotein and module eigengene trajectories throughout AD progression (see Materials and Methods, Supplementary Material, Fig. S7), demonstrating that AD-associated transcriptomic cWGCNA modules show similar expression patterns at the protein level. In these integrated networks, the edges between genes (nodes) represent both gene co-expression and PPI (Fig. 4A, N), permitting us to focus on hub genes observed at both the RNA and protein level.
Consensus transcriptomics analysis identifies neuron-specific changes in disease. (A–B) Co-expression plot (A) and gene ontology term enrichment (B) for neuronal consensus module CM1. (C–D) Module eigengene trajectory with pathological state (C) and with Braak stage (D) for CM1 in Mayo TC dataset. (E–H) Module eigengene trajectory with pathological state (E), Braak stage (F), CERAD score (G) and MMSE score (H) for CM1 in ROSMAP PFC dataset. (I–L) Module eigengene trajectory with pathological state for CM1 in MSSM PHG (I), MSSM STG (J), MSSM IFG (K) and MSSM FP (L). (M) Module eigengene trajectory with age for CM1 in NABEC dataset. (N–O) Co-expression plot (N) and gene ontology term enrichment (O) for neuronal consensus module CM10. (P–Q) Module eigengene trajectory with pathological state (P) and with Braak stage (Q) for CM10 in Mayo TC dataset. (R–T) Module eigengene trajectory with pathological state (R), Braak stage (S) and MMSE score (T) for CM10 in ROSMAP PFC dataset. (U) Module eigengene trajectory with age for CM10 in NABEC dataset.
Consistent with the neuronal dysfunction and neuronal loss seen in AD, we determined the downregulated consensus module CM1 was highly enriched in neuronal markers and our neuronal cell-type clusters (Fig. 3G). CM1’s hub genes (CALM1, CALM2, CALM3, CAMKV, CLTC and ATP2B1) are indicative of calcium signaling dysregulation, providing evidence for the calcium hypothesis of AD (37) (Fig. 4A, Supplementary Material, Table S2). GO analysis of CM1 demonstrated enrichment of terms like synaptic transmission and axonogenesis, further supporting its neuronal specificity (Fig. 4B). One of CM1’s members was also PTK2B, which has been previously identified as an AD risk gene (7). In the Mayo TC dataset we found that CM1’s module eigengene significantly decreased with late stage pathology (P < 0.001) but was unchanged with early stage pathology (P > 0.05; Fig. 4C), and the module’s trajectory declined with increasing Braak stage (R = −0.47, P = 1.4e-06; Fig. 4D). Additionally, we examined CM1’s trajectory in the ROSMAP PFC dataset and found a similar but weaker decline with increasing pathology (Fig. 4E–G). We also saw that CM1’s module eigengene increased with MMSE score (R = 0.13, P = 1.1e-03; Fig. 4H). However, similar to the regional differences we found with cell-type proportions, we noticed that the decrease in CM1 module eigengene was greater in the MSSM datasets from temporal cortical regions (parahippocampal gyrus and superior temporal gyrus) compared to those from frontal cortical regions (inferior frontal gyrus and frontal plate; Fig. 4I–L). Further, we wanted to confirm that CM1 was associated with disease state rather than normal human aging, so we calculated a synthetic module eigengene in a North American Brain Expression Consortium (NABEC) dataset, comprising brain gene expression data from hundreds of normal human patients between the age of 16 and 100 years old (38). We found that CM1 was not correlated with age (adjusted R2 = 0.031, Fig. 4M), supporting that the co-expression changes are specific to AD.
We then performed similar analyses on another downregulated consensus module, CM10, and we likewise identified this module as neuronal. One of CM10’s hub genes is APP, classically studied in AD because it transcribes the precursor of Aβ (Fig. 4N). Notably, CM10 is enriched in GO terms related to proteasomal degradation and ubiquitination (Fig. 4O), and CM10’s hub gene CAND1 regulates neddylation (39). Neddylation has been found dysregulated in AD and may impact ubiquitination (39). Like CM1, we found the module eigengene decreased with late stage pathology (P < 0.01), but we also surprisingly discovered CM10’s module eigengene tended to increase with early stage pathology (P > 0.05; Fig. 4P). This may suggest that in the early stages of disease, neurons attempt to degrade pathological protein aggregates. Again, we saw that the module’s declining trajectory was more pronounced in the temporal cortex in comparison to the frontal cortex (Fig. 4Q–T, Supplementary Material, Fig. S8) and CM10 was disease specific, rather than an effect of aging (adjusted R2 = 0.093, Fig. 4U).
In addition, we examined consensus module CM7, which we also found neuronal and downregulated in AD. CM7’s hub genes include ATP synthase subunits (ATP5G1, ATP5D, ATP5L, etc.), consistent with CM7’s enrichment of GO terms related to mitochondrial function and ATP synthesis (Fig. 5A–B). While mitochondrial dysfunction and oxidative stress have been studied in neurodegeneration in the context of both neuronal and non-neuronal cell types (40–43), mitochondria in the synapse, specifically, have been found affected in AD (44), and the gene members of CM7 may contribute to this phenomenon. Similar to CM1 and CM10, we saw CM7’s module eigengene declined with increasing pathology, and this decrease was greater in datasets from the temporal cortex compared to those from the frontal cortex (Fig. 5C–G, Supplementary Material, Fig. S8). Moreover, CM7 also showed no correlation with age, indicating the module as disease specific (adjusted R2 = 0.0072, Fig. 5H).
Consensus transcriptomics analysis identifies neuron-specific mitochondrial changes and glial-specific immune-related changes in disease. (A–B) Co-expression plot (A) and gene ontology term enrichment (B) for neuronal consensus module CM7. (C–D) Module eigengene trajectory with pathological state (C) and with Braak stage (D) for CM7 in Mayo TC dataset. (E–G) Module eigengene trajectory with pathological state (E), Braak stage (F) and MMSE score (G) for CM7 in ROSMAP PFC dataset. (H) Module eigengene trajectory with age for CM7 in NABEC dataset. (I–J) Co-expression plot (I) and gene ontology term enrichment (J) for non-neuronal consensus module CM8. (K–L) Module eigengene trajectory with pathological state (K) and with Braak stage (L) for CM8 in Mayo TC dataset. (M–O) Module eigengene trajectory with pathological state (M), Braak stage (N) and MMSE score (O) for CM8 in ROSMAP PFC dataset. (P) Module eigengene trajectory with age for CM8 in NABEC dataset. (Q–R) Co-expression plot (Q) and gene ontology term enrichment (R) for non-neuronal consensus module CM9. (S–T) Module eigengene trajectory with pathological state (S) and with Braak stage (T) for CM9 in Mayo TC dataset. (U–W) Module eigengene trajectory with pathological state (U), Braak stage (V) and MMSE score (W) for CM9 in ROSMAP PFC dataset. (X) Module eigengene trajectory with age for CM9 in NABEC dataset.
Glial-specific co-expression changes in AD
We next inspected modules that were classified as non-neuronal, specifically consensus modules, CM8, CM9, CM4 and CM5. These modules were positively correlated with AD diagnosis, and we identified CM8 as highly enriched in the microglia cell-type cluster from our snRNA-seq (Fig. 3G). Its hub genes (C1QA, C3, CASP1, CD14, etc.) and GO term enrichment of immune-related mechanisms support a critical role of neuroinflammation in AD (Fig. 5I–J). Microglia have been shown to secrete cytokines that induce a deleterious astrocytic state (A1)—one of these cytokines is C1Q (45). Similarly, members of consensus module CM9 included genes involved in some of the immune-related pathways recently implicated in AD, such as TGFβ/BMP (TGFBR2, BMPR1B), JAK/STAT (STAT3) and MAPK signaling (EGFR, FGFR1, FGF2; Fig. 5Q) (46,47). GO analysis indicated enrichment of terms related to positive regulation of signaling pathways (Fig. 5R). CM9 was highly enriched in astrocyte and endothelial cell-type clusters (Fig. 3G), and the JAK/STAT3 pathway previously has been shown to mediate astrocyte reactivity in AD (46).
In the Mayo TC dataset, both CM8 and CM9’s module eigengenes significantly increased with late stage pathology (P < 0.001), but while CM9 significantly increased in early stage pathology (P < 0.05), there was little change for CM8 (P > 0.05, Fig. 5K, S). In addition, we again saw that the changes were more prominent in the temporal cortex than the frontal cortex (Fig. 5, Supplementary Material, Fig. S8), and both CM8 and CM9 demonstrated little correlation with age, supporting that these immune-related changes are specific to disease (adjusted R2 = 0.0014, adjusted R2 = 0.056; Fig. 5P, X).
Further, consensus module CM4 contained gene members involved in various signaling pathways (NOTCH1, NFKB, ROCK1, RHOA, SMAD2, etc.; Fig. 6A), and GO analysis demonstrated enrichment of terms related to transcriptional regulation (Fig. 6B). With our snRNA-seq data, we determined this module to be enriched in OPC and both oligodendrocyte cell-type clusters (Fig. 3G), and this provided us a better functional understanding of CM4. It has been shown that the RhoA-ROCK pathway and Notch1 inhibit OPC differentiation (48–51); however, another member of CM4, ITGB1, has been demonstrated to promote myelination (52). Smad signaling also has been implicated in myelination (49,51). We additionally found that CM4’s module eigengene increased with pathology, and this increase again was greater in datasets from the temporal cortex than those from the frontal cortex (Fig. 6C–G, Supplementary Material, Fig. S8). Moreover, we found CM4 was not dependent on aging but was specific to disease (adjusted R2 = 0.018; Fig. 6h).
On the other hand, CM5’s GO term enrichment and genes (HNRNPA2B1, PAN2, CLK2, etc.) were associated with microRNA (miRNA) processing (Fig. 6I–J). miRNAs have been largely studied in the context of development (53). For example, two of miRNA-125’s predicted targets are SMAD4 and RHOA (members of CM4) (51), but miRNAs have been more recently demonstrated to play roles also in the immune system and neurodegeneration (34,54,55). Alterations in miRNA processing may mediate differential expression of miRNAs in AD. Similar to the previously mentioned non-neuronal modules, CM5’s module eigengene increased with late stage pathology (P < 0.01); however, the module eigengene also tended to decrease with early stage pathology (P > 0.05, Fig. 6K). We again saw that the increase in CM5’s module eigengene with late stage pathology was more pronounced in the temporal cortex than frontal cortex and that CM5 showed little correlation with age (adjusted R2 = 0.11; Fig. 6L–P, Supplementary Material, Fig. S8). The remaining modules can be found in Supplementary Material, Figure S9.
Assessment of genetic risk within consensus modules
The previously described transcriptomic alterations can be a cause or consequence of the disease. To address this question, we used GWAS SNPs as causal anchors and integrated our module-level information to find which modules harbor these causal anchors. We examined the enrichment of GWAS genes in our modules using the Multi-marker Analysis of GenoMic Annotation (MAGMA) pipeline (56), which controls for confounds like gene length. Each gene was assigned a score based on the best P-value of a single nucleotide polymorphism from an AD GWAS study (5), and MAGMA P-value <0.05 was used as the cut-off to include a gene as an AD candidate. We found that AD candidate genes were highly enriched in the microglial module CM8, further asserting a key role of the immune system in AD (Fig. 7A, Supplementary Material, Table S5). Several candidate genes, including CSF2RB, MS4A6A, CD33, SPI1 and CD14, were mapped in CM8, and these have been previously implicated in AD (57–60). In addition, our findings were consistent with a previous study (61) in which they reported AD risk genes in a microglial module. We wanted to further confirm our results by examining the largest number of GWAS enrichment analysis to date, including other neurodegenerative and neurological disorders (7,62–99). We found that CM8 was enriched only in GWAS hits from an AD study and a multiple sclerosis study (Supplementary Material, Fig. S10). Interestingly, CM8 was not enriched in risk genes associated with other tauopathies (frontotemporal dementia and progressive supranuclear palsy).
In addition, we validated our results in data from genome-wide gene-based association analysis (GWGAS), employed by Jansen et al., using MAGMA on a much larger cohort comprising 71 880 AD cases and 383 378 controls (6). Similar to our GWAS enrichment scores, we found that GWGAS genes were enriched in CM8 (Fig. 7B). Moreover, we inspected brain-specific eQTLs and mQTLs overlapping with AD GWAS SNPs (AD-associated brain eQTLs and mQTLs, see Materials and Methods) and found that CM8 was enriched in AD-associated brain eQTLs. On the other hand, astrocytic module, CM9, was enriched in AD-associated brain mQTLs; these results again highlight glial immune-related processes in AD.
Integration with orthogonal datasets and drug databases
We sought to understand the regulatory landscape of these AD-associated modules and integrated multiple layers of data to gain novel insights. We performed a transcriptome-wide association study (TWAS) to identify genes with cis-regulated expression associated with pathological and clinical traits of AD (100). We used plaque and tangle density as quantitative measures of the pathological changes in AD along with clinical dementia rating (CDR) scores to find significant transcriptome-wide associations. Using hypergeometric overlap tests, we found that downregulated cWGCNA modules—CM1, CM7, CM10 and CM12—were enriched in TWAS genes anti-correlated with plaque and tangle density, and these modules also were enriched in TWAS genes positively correlated with cognitive decline (Fig. 7C). Conversely, upregulated modules—CM4, CM16 and CM23, for example—were enriched in TWAS genes positively correlated with plaque and tangle density and TWAS genes negatively correlated with cognitive decline. These findings confirm our module trajectories with pathology and cognition previously discussed, but interestingly, module enrichment of TWAS genes associated with tangle density was not only less but also sometimes incongruent to module enrichment of TWAS genes associated with plaque density. Interestingly, microglial module CM8 was enriched in TWAS genes positively correlated with plaque density but was not enriched in TWAS genes associated with tangle density. This suggests that the upregulation of cis-regulated microglial immune-related genes may be either driving or affected by Aβ pathology but not by tau pathology. On the other hand, CM7, which we identified as neuron-specific and mitochondrial-related, was highly enriched in TWAS genes anti-correlated with plaque density but is also enriched in TWAS genes positively correlated with tau density, indicating contrasting relationships between tau and Aβ pathology with neuronal mitochondrial function.
As a further step to characterize our modules and to identify co-expression changes associated with epigenetic changes, we leveraged histone H3 lysine 9 acetylation (H3K9ac) marks, which are found near active transcription start sites and enhancers, associated with plaque or tangle density (101). We found that all AD-associated modules except microglial module CM8 were enriched in H3K9ac marks significantly associated with tangles (Fig. 7D). In addition, only neuronal module CM1 and astrocytic module CM9 were enriched in plaque-associated H3K9ac marks at promoter regions. We also examined our consensus modules with a dataset of AD-associated DNA methylated regions (102,103) and found only upregulated modules CM5, CM16 and CM23 were enriched in AD-associated hypo-methylated DNA regions. However, downregulated modules CM10, CM12, CM1 and CM7 were enriched in hyper-methylated DNA regions. The notable lack of CM8 enrichment in both AD-associated H3K9ac marks and AD-associated DNA methylation hints at an intriguing possibility that microglial immune-related changes in disease may not be mediated by epigenetic changes.
Additionally, we queried the Library of Integrated Network-Based Cellular Signatures (LINCS) database to identify potential drugs to target each consensus module (see Materials and Methods, Supplementary Material, Fig. S11). For modules upregulated in AD, we queried for drugs that would inhibit module expression, whereas for modules downregulated in AD, we queried for drugs that activate module expression. Notably, we found several drugs previously reported to be beneficial in AD and other neurodegenerative diseases. For example, histone deacetylase inhibitor Vorinostat, found to modulate CM9, is currently in human clinical trials (NCT03056495). Both saracatinib, an inhibitor of the Src/abl family of kinases and a modulator of CM8, and resveratrol, a modulator of CM5, have been found to have beneficial effects in human clinical trials (104,105).
Identification of TFs regulating consensus modules
Network analysis identifies co-expressed modules that are likely co-regulated (33,106,107). With this in mind, we performed transcription factor binding site (TFBS) enrichment analysis (see Materials and Methods) to identify TFs that may be regulating disease-associated modules. This method has been successfully used previously to accurately predict TFs (32). TFBS revealed that master TF—RELA, otherwise known as p65, appears to control both neuronal and non-neuronal modules (Fig. 7E–F). It becomes further evident that NFκB signaling may be mediating transcriptional changes in AD since we found TFs NKFB1 and RELB associated with non-neuronal modules. Previous studies have reported increased NFκB activation in AD and pointed to a relationship with Aβ pathology (108–111). Moreover, Withaferin A has been identified as an inhibitor of NFκB signaling (112,113) and has shown therapeutic potential to protect against Aβ toxicity (114), and we found Withaferin A is a modulator of CM9, CM4 and CM1 (Supplementary Material, Fig. S11).
Validation of consensus modules in external datasets
Since modules can be affected by non-biological variables like differences in sample preparation, we needed to authenticate the existence of our modules of interest in other datasets.
In order to do this, we utilized module preservation analysis, which helps us quantify the features of a particular module that are preserved across multiple networks (see Materials and Methods). We found that our consensus modules were highly preserved across published human RNA-seq and microarray datasets (61,115,116), indicated by Z-summary preservation (Fig. 8A, Supplementary Material, Table S6). Moreover, we examined the module eigengene trajectories of modules—CM10, CM1 and CM8—in these datasets and found comparable changes with diagnosis (Fig. 8B-D). As expected and also in line with our findings, we observed that changes were notably greater in the temporal cortex than the frontal cortex, again reflecting AD’s progressive atrophy of different brain regions (1) (Supplementary Material, Fig. S12). These results altogether support the robustness of our findings.
Consensus transcriptomics analysis identifies non-neuronal changes in disease. (A–B) Co-expression plot (A) and gene ontology term enrichment (B) for non-neuronal consensus module CM4. (C–D) Module eigengene trajectory with pathological state (C) and with Braak stage (D) for CM4 in Mayo TC dataset. (E–G) Module eigengene trajectory with pathological state (E), Braak stage (F) and MMSE score (G) for CM4 in ROSMAP PFC dataset. (H) Module eigengene trajectory with age for CM4 in NABEC dataset. (I–J) Co-expression plot (I) and gene ontology term enrichment (J) for non-neuronal consensus module CM5. (K–L) Module eigengene trajectory with pathological state (K) and with Braak stage (L) for CM5 in Mayo TC dataset. (M–O) Module eigengene trajectory with pathological state (M), Braak stage (N) and MMSE score (O) for CM5 in ROSMAP PFC dataset. (P) Module eigengene trajectory with age for CM5 in NABEC dataset.
Integration of orthogonal datasets to identify associated gene regulatory programs. (A) Mean scaled enrichment of AD genome-wide association study (GWAS) hits from Kunkle et al. (5) in consensus modules. (B) Heat map showing enrichment of consensus modules in genome-wide gene-based association analysis (GWGAS) (6) hits and AD-associated brain expression and methylation quantitative trait loci (eQTLs and mQTLs). (C) Heat map indicating enrichment of consensus modules with genes identified through transcriptome-wide association analysis (TWAS) (100). TWAS genes were measured in relation to plaque and tangle density, as well as clinical dementia rating (cognitive decline). (D) Heat map showing enrichment of consensus modules with H3K9 acetylation marks, correlated with tangle and plaque burden, in promoter, enhancer or other genomic regions and AD-associated hypermethylated and hypomethylated DNA regions. (E) Co-expression plot demonstrating demonstrates relationships between neuronal consensus modules and transcription factors. Red edges reflect low correlation and green reflect high correlation. (F) Co-expression plot demonstrating demonstrates relationships between non-neuronal consensus modules and transcription factors.
Validation of findings in external datasets. (A) Heat map indicating module preservation in other human RNA-seq and microarray datasets. (E–I) Boxplot showing module eigengene trajectory with diagnosis of consensus modules CM10, CM8 and CM1 in datasets Zhang et al. PFC (61) (E), Berchtold et al. superior frontal gyrus (SFG, F) and post-central gyrus (115) (PCG, G), and Webster et al. FC (H) and TC (116) (I). The upper and lower lines represent the 75th and 25th percentile, respectively, and the central line represents the median.
Discussion
Previous AD transcriptomic studies have been limited to tissue from a single brain repository or a single brain region (11,14,15); we performed cWGCNA (17,18) on three bulk tissue RNA-seq datasets (ROSMAP, Mayo and MSSM) in order to identify robust, conserved and biologically important transcriptomic alterations in AD across six cortical regions. We then integrated various genetic and epigenetic datasets, along with our own snRNA-seq data, in order to provide a more comprehensive picture of the genome-wide transcriptomic alterations occurring in AD and to begin to examine the associated gene regulatory programs.
A limitation of bulk tissue RNA-seq is that it does not address cell-type heterogeneity-associated changes that occur in disease, thus in order to probe the cell-type-specific transcriptional changes occurring in AD, we generated single-nuclei transcriptomic profiles on more than 27 000 nuclei from the aged human prefrontal cortex. We used a novel algorithm in order to estimate the cell-type abundance in bulk tissue RNA-seq data with our single-nuclei transcriptomic dataset. With this approach, we were able to show that the abundance of several excitatory and inhibitory neuronal clusters decreased with late stage pathology, whereas the abundance of microglial and astrocyte clusters increased. In addition, these changes in cell-type proportion were dependent on the brain region, modeling the progression of AD pathology in the brain (1,29).
Furthermore, we integrated our snRNA-seq data in order to help classify the cWGCNA modules based on cell type, allowing us to determine that AD co-expression changes related to calcium signaling (CM1), ubiquitination (CM10) and mitochondria function (CM7) are neuron specific. Gene set enrichment analysis showed that CM4 was associated with transcriptional regulation, but integration of our snRNA-seq data helped to clarify that this module was related to myelination and oligodendrocyte maturation. While white matter changes have been previously implicated in AD, the associated molecular changes occurring have been largely unstudied (11,117–121), and a recent snRNA-seq study of AD also has highlighted myelination-related genes (122), supporting the importance of future work into understanding the gene members of CM4.
We also were able to confirm the role of astrocytes and microglia (CM9 and CM8) in immune-related processes in AD. Notably, we found that the microglial module CM8 was highly enriched in genetic risk signals associated with AD but not other tauopathies, suggesting that the concurrence of Aβ and tau pathology may be important. However, recently, a transcriptomic study of microglia isolated from APP and tau transgenic mice investigated the differential effects of Aβ and tau pathology, finding that Aβ pathology results in an exacerbated inflammatory response from microglia, unseen in the tau transgenic-derived microglia (123). Similarly, we found that microglial module CM8 was enriched in TWAS genes associated with plaque density but not TWAS genes associated with tangle density, reinforcing the relationship between Aβ pathology and microglial immune-related transcriptional changes.
Moreover, we wanted to clarify the epigenetic changes associated with AD co-expression changes, examining AD-associated DNA methylation regions and mQTLs, along with H3K9ac marks associated with tangle and plaque density. Interestingly, despite finding other modules significantly enriched in H3K9ac marks associated with pathology, we found that microglial module CM8 was not enriched in H3K9ac marks associated with either tangle or plaque density. Accordingly, CM8 was not enriched in AD-associated hypermethylated or hypomethylated DNA regions or AD-associated mQTLs, suggesting a lack of epigenetic regulation of CM8. Histone deacetylase inhibitors, however, have been shown to suppress microglial activation, which is associated with transcriptional changes (124,125), indicating that further work to elucidate the epigenetic mechanisms underlying AD is needed.
Since aging is the most important non-genetic risk factor for AD, we wanted to investigate if the AD-associated consensus modules are altered during normal human aging. Normal aging has been shown to alter glial and neuronal gene-expression changes in the brain (38); however, we hypothesized that these changes are distinct from pathological changes occurring in neurodegenerative diseases. By calculating synthetic module eigengenes from AD-associated consensus modules in a NABEC dataset (38), we were able to confirm that the AD-associated modules were not significantly correlated with normal aging.
In addition, one of the biggest challenges in co-expression analysis is to find robust co-expressed modules that are altered in different cohorts and datasets. In order to tackle this problem, we used module preservation analysis, which is a powerful tool to find conserved modules in different datasets. Through module preservation analysis, we found that our AD-associated consensus modules are robustly preserved in multiple published AD datasets from different brain regions (61,115,116). This finding in combination with the modules’ lack of correlation with normal human aging verified that we have identified disease-specific transcriptomic signatures.
Altogether, our study emphasizes the value of integrative analyses in order to fully understand the molecular changes occurring in disease. While single-cell RNA-seq provides the advantage of profiling gene expression changes at the cellular level, it is still not only expensive but also limited to mRNA. Thus, at this time, analyses incorporating both single-cell and bulk tissue transcriptomic data are important for thorough investigation of AD transcriptomic alterations. Additionally, we have examined epigenetic and TF regulation by integrating multiple orthogonal datasets, providing a framework for disease-associated gene regulatory networks in AD and further highlighting a need to employ integrative analyses of multiple types of data. We also hope that our findings will help pave way for accelerated drug development, a critical need in AD, since our study has identified robust co-expression modules and their hubs, in addition to TFs, which may be regulating these modules, all of which can be targets for novel therapeutics.
Materials and Methods
Data and code availability
RNA-seq data and details can be found on https://www.synapse.org using corresponding synapse (syn) IDs- Mayo RNA-seq (syn5550404), Mount Sinai Brain Bank (MSBB) (syn3159438) and ROSMAP (syn3219045). Single-cell data can be found on synapse website using synID# syn18915937. Additionally, we have created shiny apps to easily visualize gene, network-level and single-cell data, which can be accessed through our lab website at https://swaruplab.bio.uci.edu/consensusAD. Code for consensus WGCNA can be found here: https://github.com/vswarup/ConsensusWGCNA.
Description of human samples and RNA-seq data
For single-cell analysis control human prefrontal cortex brain samples were obtained from UCI MIND’s Alzheimer’s Disease Research Center brain bank. These samples were obtained under UCI’s Institutional Review Board. Post-mortem human brain tissues were obtained from three cohorts—Mayo study, MSBB study and ROSMAP study. All the human samples were obtained under their respective institutional review board. For details on human samples used in this study and their availability, please see Supplementary Material, Table S1. RNA-seq data and details can be found on synapse.org website using corresponding synapse (syn) IDs—Mayo RNA-seq (syn5550404), MSBB (syn3159438) and ROSMAP (syn3219045).
Other published datasets used in the study
For normal human aging, we used NABEC data ranging from 16 to 101 years old from frontal cortex (38,126,127). The gene-expression arrays are available from Gene Expression Omnibus (GEO)—GSE36192. Published AD microarray datasets were also downloaded from GEO including Zhang et al. prefrontal cortex data (61) (GSE44770), Berchtold et al. frontal and temporal cortex data (115) (GSE48350) and Webster et al. (116) (GSE15222).
RNA-seq data processing
Consensus WGCNA
cWGCNA was performed with ROSMAP, Mayo and MSSM datasets using WGCNA package in R (17). We employed a signed cWGCNA approach by calculating component-wise values for topological overlap for individual brain banks. First, bi-weighted mid-correlations were calculated for all pairs of genes, and then a signed similarity matrix was created. In the signed network, the similarity between genes reflects the sign of the correlation of their expression profiles. The signed similarity matrix was then raised to power β (β = 12) to emphasize strong correlations and reduce the emphasis of weak correlations on an exponential scale. The resulting adjacency matrix was then transformed into a topological overlap matrix as described elsewhere (33). Modules were defined using specific module cutting parameters, which included minimum module size of 100 genes, deepSplit = 4 and threshold of correlation = 0.2. Modules with correlation greater than 0.8 were merged together. We used first principal component of the module, called module eigengene, to correlate with diagnosis and other variables. Hub genes were defined using intra-modular connectivity (kME) parameter of the WGCNA package. Gene-set enrichment analysis was done using enrichR package (128). Potential drugs targeting consensus modules were identified using enrichR to query the LINCS (www.clue.io) chemical perturbations database.
Comparison to other gene-clustering algorithms
In addition to cWGCNA, we constructed gene modules using several alternate methods to demonstrate that these modules are detected using different techniques. For this analysis, we used samples from the Mayo cohort to construct gene modules using DICER, DiffCoEx, MEGENA and self-organizing maps (SOMs). We compared the overlap between the modules resulting from these algorithms to the top 200 genes in each cWGCNA module by intra-modular connectivity (kME) using a hypergeometric test (phyper function in R). Default parameters were used for DICER (v0.6.0) and MEGENA (v1.3.7) in their respective R packages. To perform DiffCoEx, we used the diffcoex function from the R package MODifieRDev (https://github.com/ddeweerd/MODifieRDev) with default parameters. We constructed a SOM of the scaled gene expression data using the som function from the R package kohonen (v 3.0.10) using a grid of 40 by 40 nodes, whereby each gene in the expression dataset is assigned to one of these nodes throughout the SOM training routine. Gene modules are then constructed by hierarchically clustering the resulting SOM nodes using the hclust function in R. We then visualized bi-weighted mid correlations between genes and relevant biological variables such as AD diagnosis, age at death and sex on the SOM to ensure that the SOM gene modules were capturing relevant biological features.
Proteomic Module Eigengenes—Processed proteomic AD dataset (36) was downloaded from synapse website (syn20933797). Consensus AD transcriptomic module definitions were used to calculate synthetic module eigengenes using the moduleEigengenes function in WGCNA package and plotted as boxplots stratified by control, asymptomatic AD and AD samples.
Co-expression PPI networks—Co-expression PPI networks were created as described previously (34). Briefly, we used two PPI resources, InWeb (30,129) and BioGRID (31) and using the union of the two resources, the subset of compiled PPIs between genes in a module was extracted and all edges were counted. The PPI dataset was matrix-multiplied with co-expression data from the RNA-seq data, the edges that were present both in PPI and co-expression datasets were eventually kept. This approach allowed us to infer tissue and species specificity in the PPI network. The PPI network was then visualized using the igraph package in R.
Module Preservation Analysis
To understand the functional relationship of modules and its network properties in orthogonal datasets, we performed module preservation analysis. Module definitions from the cWGCNA analysis was used and Z-statistics were calculated using the module preservation function in WGCNA package in R (130).
Genetic enrichment with eQTLs, mQTLs and GWGAS
Brain specific eQTLs were generated from databases including the CommonMind Consortium Portal, AMP-AD, BRAINEAC and xQTL Serve. A false discovery rate of 0.05 was used to define significant eQTL associations. eQTLs that overlapped with AD associated SNPs from either Lambert et al. or Jansen et al. GWAS datasets (6,7) were annotated as AD-associated eQTLs. Similarly, methylation QTLs (mQTLs) defined in AMP-AD, and CommonMind Consortium Portal were annotated as AD-associated mQTLS if they overlapped with AD GWAS data. We also used GWGAS data using MAGMA (56). GWAS enrichment analysis in modules was done using MAGMA pipeline. Briefly, GWAS summary statistics was obtained for AD from IGAP website (7). Using MAGMA approach, we assigned each gene a score based on the best P-value of an SNP in a given GWAS study within 20 kB of the gene, and then set a P-value cut-off at 0.05 to define the gene as included in the common variant set related to AD. We performed enrichment analysis with logistic regression controlling for gene length, LD blocks and other biases (33).
TWAS Enrichment Analysis
We performed AD TWAS using the FUSION package (http://gusevlab.org/projects/fusion/) (100) with gene-expression measured in ROSMAP PFC data and GWAS summary statistics from AD GWAS (7). We chose the best model of accuracy in the FUSION package (best cis-eQTL, best linear unbiased predictor, Bayesian sparse linear mixed model (131), Elastic-net regression, LASSO regression) and used for downstream association analyses. TWAS association statistics were Bonferroni corrected per GWAS, gene and intron separately. Association for TWAS was assessed for plaque and tangle densities and genome-wide significant genes were used for enrichment analysis. Similarly, CDR reported in ROSMAP dataset was used for association with cognitive decline.
Single-cell Library Preparation
Five control human brain samples were processed for nuclei isolation following modified version of the published protocol (132). All procedures were carried on ice. Approximately 50 mg of the brain sample was homogenized in 500 μl chilled Nuclei EZ lysis buffer (Sigma-Aldrich), using a motorized tissue grinder and disposable pellet pestles (Fisher Scientific) by stroking ~10–20 times. The homogenate was then incubated for 5 min after adding 1 ml of the Nuclei EZ lysis buffer. Homogenate was then filtered using a 70 μm MACS strainer (Miltenyi Biotec). The flow through was centrifuged at 500 g for 5 min at 4°C and the supernatant was carefully removed. The nuclei pellet was then re-suspended in another 1.5 ml of Nuclei EZ lysis buffer and incubated for 5 min. The cell suspension was again centrifuged (500 g, 5 min, 4°C) and the supernatant was removed. A total of 500 μl of Nuclei wash and re-suspension buffer (NWR buffer: 1× PBS (Fisher Scientific), 1% BSA (Sigma-Aldrich) and 0.2 U/μl RNase inhibitor (New England Bio Labs)) was then added and incubated without re-suspending, to allow buffer interchange. After the incubation, 1 ml of NWR buffer was added and the nuclei were re-suspended. The resuspension was centrifuged (500 g, 5 min, 4°C) and the supernatant was removed. The nuclei pellet was then gently resuspended in 1.4 ml of NWR buffer, centrifuged (500 g, 5 min, 4°C), removed supernatant and was gently re-suspended in another 500 μl of NWR buffer. All nuclei were collected by washing the walls of the centrifuge tube. Nuclei re-suspension was then filtered using 40-μm Flowmi cell strainer (Bel-Art). The number of nuclei was then counted using an automated cell counter (Bio-Rad). Samples were directly sorted by FACS using DAPI positive gates at the UCI Immunology FACS core in the 10X buffer and libraries were prepared using Single-cell RNA-seq library kit (v3, 10X Genomics). Samples were pooled and sequenced at an average read-depth of 100 000 reads/cell.
Single-cell Clustering and analysis
Raw sequencing Illumina bcl2 files were demultiplexed using 10X Genomics Cell Ranger software. Further analysis was performed in Seurat V3 package in R (133). We used sctransform function in Seurat to normalize the data and top 75 PCs were used for clustering and visualization using Fourier transformation t-SNE (FItsne (134)). Differential expression between clusters were performed using MAST package in Seurat (135). Additionally, differential expression was performed to compare cells originating from male and female samples within each cluster using MAST.
Cell-type composition analysis
For cell-type composition analysis, unregressed FPKM data which did not account for cell-type principal components, were used. Using the most variable genes in a cluster, we performed WGCNA analysis and determined the intramodular connectivity (kME) of each gene. kME values are correlation of first principal components of a module/cluster to the expression of that gene. Thus, higher kME values can often explain most of the variance in the data because of their high correlation to PC1. Using Fisher’s transformation, we converted kME values to specificity scores (Z-scores). We then used the top one percentile of genes by their kME values and calculated their module eigengene in respective bulk RNA-seq dataset. These module eigengene values were then plotted to understand the cell-composition changes in control and AD samples.
Transcription factor binding site enrichment
TFBS enrichment was done using pipeline published elsewhere (106). Briefly, the computational pipeline identifies the TFs (in the promoter region) in the given gene list using experimentally validated TFs. The pipeline used TF position-weighted matric from JASPAR and TRANSFAC databases as well as included experimental associated datasets obtained from ENCODE (136) and ChEA (137). The pipeline uses clover algorithm and corrects for background using human chromosome 22.
Acknowledgements
This work was supported by University of California Irvine startup funds and American Federation of Aging Research (AFAR) young investigator award to V.S. For ROSMAP dataset, the study data were provided by the Rush Alzheimer’s Disease Center, Rush University Medical Center, Chicago. Data collection was supported by National Institute on Aging [P30AG10161, R01AG15819, R01AG17917, R01AG30146, R01AG36836, U01AG32984, U01AG46152]; the Illinois Department of Public Health; and the Translational Genomics Research Institute. The results published here are in whole or in part based on data obtained from the AMP-AD Knowledge Portal (doi:10.7303/syn2580853). For Mayo dataset, the study data were provided by the following sources: The Mayo Clinic Alzheimer’s Disease Genetic Studies, led by Dr Nilufer Taner and Dr Steven G. Younkin, Mayo Clinic, Jacksonville, FL, using samples from the Mayo Clinic Study of Aging, the Mayo Clinic Alzheimer’s Disease Research Center and the Mayo Clinic Brain Bank. Data collection was supported by National Institute on Aging [AG016574, R01 AG032990, U01 AG046139, R01 AG018023, U01 AG006576, U01 AG006786, R01 AG025711, R01 AG017216, R01 AG003949]; National Institute of Neurological Disorders and Stroke [R01 NS080820]; CurePSP Foundation; and support from Mayo Foundation. Study data include samples collected through the Sun Health Research Institute Brain and Body Donation Program of Sun City, Arizona. The Brain and Body Donation Program is supported by the National Institute of Neurological Disorders and Stroke [U24 NS072026 to National Brain and Tissue Resource for Parkinson’s Disease and Related Disorders]; the National Institute on Aging [P30 AG19610 to Arizona Alzheimer’s Disease Core Center]; the Arizona Department of Health Services (contract 211002, Arizona Alzheimer’s Research Center); the Arizona Biomedical Research Commission (contracts 4001, 0011, 05-901 and 1001 to the Arizona Parkinson’s Disease Consortium); and the Michael J. Fox Foundation for Parkinson’s Research. The results published here are in whole or in part based on data obtained from the AMP-AD Knowledge Portal (doi:10.7303/syn2580853). For MSSM dataset, the data were generated from postmortem brain tissue collected through the Mount Sinai VA Medical Center Brain Bank and were provided by Dr Eric Schadt from Mount Sinai School of Medicine. The results published here are in whole or in part based on data obtained from the AMP-AD Knowledge Portal (doi:10.7303/syn2580853). In addition, we are grateful to Thanneer Perumal and Ben Logsdon at Sage bionetworks for processing the fastq and bam files used in the analysis.
Conflict of Interest statement. The authors declare no conflict of interests.
