Integrative network analysis reveals novel moderators of Aβ-Tau interaction in Alzheimer's disease

Abstract

Background

Although interactions between amyloid-beta and tau proteins have been implicated in Alzheimer's disease (AD), the precise mechanisms by which these interactions contribute to disease progression are not yet fully understood. Moreover, despite the growing application of deep learning in various biomedical fields, its application in integrating networks to analyze disease mechanisms in AD research remains limited. In this study, we employed BIONIC, a deep learning-based network integration method, to integrate proteomics and protein–protein interaction data, with an aim to uncover factors that moderate the effects of the Aβ-tau interaction on mild cognitive impairment (MCI) and early-stage AD.

Methods

Proteomic data from the ROSMAP cohort were integrated with protein–protein interaction (PPI) data using a Deep Learning-based model. Linear regression analysis was applied to histopathological and gene expression data, and mutual information was used to detect moderating factors. Statistical significance was determined using the Benjamini–Hochberg correction (p < 0.05).

Results

Our results suggested that astrocytes and GPNMB + microglia moderate the Aβ-tau interaction. Based on linear regression with histopathological and gene expression data, GFAP and IBA1 levels and GPNMB gene expression positively contributed to the interaction of tau with Aβ in non-dementia cases, replicating the results of the network analysis.

Conclusions

These findings suggest that GPNMB + microglia moderate the Aβ-tau interaction in early AD and therefore are a novel therapeutic target. To facilitate further research, we have made the integrated network available as a visualization tool for the scientific community (URL: https://igcore.cloud/GerOmics/AlzPPMap).

Background

Alzheimer's disease (AD) is the most common form of dementia, characterized by extracellular amyloid-beta (Aβ) plaques and intracellular tau protein tangles, which are believed to be central to its pathology [1,2,3]. Recent studies have reported an interaction between Aβ and tau proteins in the pathogenesis of AD, contributing to disease progression [4,5,6,7,8,9]. For instance, Lee et al. [7] revealed that Aβ-tau interactions are associated with the propagation of tau. In practice, lecanemab, an anti-Aβ antibody, which has been shown to slow cognitive decline in patients with early-stage AD, reduces total tau protein and P-Tau181 levels in cerebrospinal fluid [10]. This suggests that Aβ contributes to the effect of tau in AD. Moreover, several phenotypes are associated with the Aβ-tau interaction, such as astrocyte reactivity, blood pressure, vascular burden, and microglia [11,12,13,14,15]; however, the detailed mechanisms are not fully understood. Various cellular perturbations have been identified through single-nucleus RNA sequencing (snRNA-seq) in AD [16,17,18,19,20,21,22]. For example, a recent snRNA-seq study [17] revealed the presence of Trem2-dependent disease-associated microglia (DAM) and a unique Serpina3n + C4b + reactive oligodendrocyte population in 5XFAD mice. In human AD, distinct glial phenotypes including IRF8-driven reactive microglia and oligodendrocytes with impaired myelination and metabolic adaptation to neuronal degeneration have been observed. However, the specific cell types and effects of Aβ-tau interactions are yet to be elucidated.

Network analyses based on gene or protein co-expression have revealed gene groups and associated pathways correlated with pathological and clinical factors in AD [23,24,25]. Additionally, network integration methods based on multi-omics and other data types have been developed [26,27,28,29]. By integrating both physical and functional interactions simultaneously, network integration analyses can identify detailed molecular relationships. However, this approach has largely been applied to cancer. AD is a complex disease involving various pathologies, and network integration is a promising approach for determining comprehensive physical and functional interactions involving Aβ and tau in AD.

In this study, we investigated moderating factors in the Aβ-tau interaction using a systematic network integration approach. The aim of this study was to elucidate the cell types and moderating factors in the interaction between Aβ and tau proteins in early-stage AD. We show the overview of this study (Supplementary Fig. S1). First, we integrated the proteomics and protein–protein interaction (PPI) networks by using BIONIC [30], a network integration algorithm. In the downstream analysis, glial cell-related proteins, including those associated with microglia and astrocytes, were enriched around APP and MAPT subnetworks, which encode the precursors of Aβ and tau, respectively. Further analyses of microglial subtypes revealed enrichment in GPNMB + microglia. As a validation, linear regression analysis based on histopathological and gene expression data showed that interactions between Aβ and IBA1, GFAP, and GPNMB, were positively related to tau levels. Our findings highlight that deep learning-based network integration can serve as a novel framework for elucidating complex disease mechanisms. The integrated network, made available as a visualization tool for the scientific community, not only enhances our understanding of the Aβ-tau interaction but also provides a basis for generating experimental validation ideas and developing more effective treatments for AD. By leveraging this tool, researchers can gain insights into the intricate molecular interactions at play in AD, fostering the discovery of new therapeutic targets.

Methods

Brain tissue proteomics data from the ROSMAP cohort

For the proteomic analysis, ROSMAP tandem mass tag (TMT) proteomics data were obtained. TMT proteomics is a technique for the simultaneous quantitative analysis of protein expression in multiple samples using LC–MS/MS and is commonly employed in AD research [24, 25]. The ROSMAP data set included TMT proteomics data for 400 cases with differences in cognitive status (Alzheimer's dementia, other dementia, MCI, or no impairment), derived from brain tissue samples. Data normalization and batch correction were conducted according to previous described procedures [24]. Four outlier samples were excluded from subsequent analyses based on a principal component analysis. In this study, ROSMAP TMT data for MCI and early AD cases were included. Early-stage AD was defined according to previously reported criteria: meeting the NIA-AA core clinical criteria for probable AD dementia and having a global CDR score of 0.5 to 1.0 and a CDR memory box score of 0.5 or greater at screening and baseline [10]. Additionally, the MMSE effectively discriminates between CDR stages 0.5, 1, 2, and 3 but performs poorly in the separation of CDR stages zero and 0.5. The previously established MMSE ranges were 30 for no, 26–29 for questionable, 21–25 for mild, 11–20 for moderate, and 0–10 for severe dementia [31]. We defined early AD as having an MMSE score of ≥ 21.

Proteomics and PPI data integration using a graph attention network-based method

To integrate the proteomics and PPI data, total values were calculated for each gene symbol in the proteomics dataset. Pearson correlation coefficients were derived from the proteomics data, retaining only data with Benjamini-Hochberg (BH)-corrected p-values of < 0.05. To validate the chosen significance thresholds (p < 0.05 with BH correction), we conducted a statistical power analysis. The mean absolute correlation coefficient among the significant results (BH-adjusted p-value < 0.05) was 0.338 and was used as the effect size. The highest unadjusted p-value (0.00301) for these results was set as the significance level. Using a sample size of 122, the calculated power was 0.811, confirming that our study design is adequately powered to detect correlations of this magnitude. Normalization was performed to set maximum values to 1 and minimum values to 0. For PPI, STRING [32] filtering was used to obtain physical interactions within Homo sapiens. To ensure the reliability of the interactions, we applied a confidence score threshold of 0.7, classifying interactions as "high confidence" based on STRING’s scoring system. This threshold, which has been widely used in previous studies [33, 34], indicates a 70% or higher probability that a predicted interaction exists within the same metabolic map in the KEGG database. By adopting this threshold, we aimed to maintain consistency with established methodologies in PPI network analysis. Normalization was performed to set the maximum values to 1 and minimum values to 0. Proteomic and PPI data were integrated using a recently reported method, BIONIC [30], which integrates multiple biological networks while preserving information using graph attention networks [35]. Unlike linear integration methods that assign fixed weights to input networks, BIONIC dynamically balances information across networks by reconstructing their topological features through a data-driven approach. Default values were used for the BIONIC hyperparameters. For both proteomics and PPI data, 4,787 proteins common to both datasets were retained. After integrating the two networks, 512-dimensional features were output for each protein.

Clustering analysis and dimension reduction using features derived from BIONIC

Considering that the data derived from neural networks are nonlinear, clustering was performed following a method similar to that used in snRNA-Seq analyses [36]. After feature calculation using BIONIC and subsequent scaling to Z-scores, Euclidean distances were computed to create a distance matrix. Then, setting k = 25, the k-nearest neighbor method was used to calculate an adjacency matrix and form a graph structure. This graph structure was then subjected to clustering using the Louvain method [37] (resolution = 1) to classify each protein into clusters.

Cluster centroid analysis using BIONIC-derived features

After implementing BIONIC to evaluate the networks derived from proteomic and PPI data, a clustering analysis of proteins was performed. To analyze the relationships between clusters, the average of the feature values for each cluster was calculated as its centroid. Pearson correlation coefficients and Euclidean distances between the clusters were calculated. Furthermore, the network of clusters was visualized using the reciprocal values of the Euclidean distances as the weights of the edges.

Cell type enrichment analysis based on human brain and microglia snRNA-seq data

For a cell type enrichment analysis, publicly available human brain [38] and microglia [21, 22] snRNA-seq data were used. Cell annotations for each dataset were based on publicly available data. The AddModuleScore function was performed using the Seurat package [39] to calculate the average expression levels in each protein community. The combined score was calculated by multiplying the expression level by the proportion of positive cells in each community.

Cell type composition analysis based on human brain and microglia snRNA-seq data

For a cell type composition analysis, publicly available human microglia [21, 22] snRNA-seq data were used. ANOVA was used to compare the frequencies of different microglial subtypes across various conditions (Control, Aβ + , and Aβ + Tau +). Next, for the microglial subtypes "Microglia_GPNMB_LPL" and "Microglia_GPNMB_PLAT," a Tukey HSD test was used to determine specific differences between groups (for cases where ANOVA showed p < 0.05).

Gene Ontology functional enrichment analysis

A Gene Ontology (GO) analysis was used to explore the biological functions of our targeted gene set. The GO biological process (GObp) dataset version 2023.1 was downloaded from the Molecular Signatures Database (MsigDB) (https://www.gsea-msigdb.org/gsea/msigdb). For the analysis, the enrichment function in the clusterProfiler package was used, and values of p < 0.05, after BH correction, were considered significant. The Rrvgo [40] software was used to summarize the GO terms.

Modulators of the Aβ-tau interaction

To explore modulators of the Aβ-tau interaction based on proteomics data, a slightly modified version of the MINDy algorithm [41] was employed. Briefly, (1) Euclidean distances for APP and MAPT were calculated using BIONIC features, and common elements within the top 5% with respect to proximity were selected for analyses. (2) For each candidate modulator, the top and bottom 35% of proteins with respect to expression levels were segregated, and mutual information for APP and Tau was calculated for each group. Mutual information, which measures the interdependence between two random variables X and Y, was defined using Eq. (1) as follows:

$$\begin{array}{c}MI\left(X, Y\right)= \sum x\in X \sum y\in Yp\left(x,y\right)\text{log}\left(p\left(x\right)p\left(y\right)p\left(x,y\right)\right)\end{array}$$

(1)

Here, p(x,y) represents the joint probability distribution of X and Y. p(x) and p(y) are the marginal probability distributions of X and Y, respectively. The logarithm base 2 was used. Subsequently, (3) the difference in mutual information between the top and bottom 35% groups was calculated using Eq. (2):

$$\begin{array}{c}\Delta MI = \left|MI\left({MAPT}_{Mtop35}, { APP}_{Mtop35}\right)- MI\left({MAPT}_{Mbottom35}, { APP}_{Mbottom35}\right)\right|)\end{array}$$

(2)

Here, M_top35 and M_bottom35 indicate the APP and MAPT groups in the top 35% and bottom 35% of protein expression levels, respectively. This differentiation is crucial for understanding the difference in mutual information between the two groups, allowing for a clear analysis of the Aβ-Tau interaction modulators. Bootstrapping with 1,000 iterations was then used to compare mutual information, with values of p < 0.05 deemed significant.

Next, factors moderating Aβ-tau interactions were extracted using data obtained from BIONIC and a network analysis was conducted. Among candidate modulators identified in the MINDy analysis, HSPA5, which is closely associated with Aβ [42] and Tau [43], was included in the analysis. The inner product of the BIONIC-derived features was calculated to determine protein similarity, retaining the top 5% for network construction. From this network, the top 5% of factors based on Euclidean distances to APP, MAPT, and HSPA5 were retained to form a subnetwork. Community detection within this subnetwork was performed using the leading eigenvector function of the iGraph package [44].

Neuropathological and RNA-seq datasets

The neuropathological and RNA-Seq datasets used in this study were obtained from previously published works [45]. Neuropathological data include quantitative immunostaining results for the microglial marker IBA1, astrocyte marker GFAP, Aβ, tau, and phosphorylated tau. The datasets included 377 samples collected from the cortical grey (parietal and temporal) and white matter (parietal) as well as the hippocampus. Additionally, the frequency distributions of age, sex, CERAD scores, and Braak scores among the participants were evaluated, differentiating between the no-dementia and dementia groups. Age was compared between the dementia and no-dementia groups using a t-test, whereas sex, Braak score, and CERAD score were compared between the groups using Fisher's exact test.

Linear regression analysis for exploring interaction terms

To explore moderating factors in the Aβ-tau interaction, a linear regression analysis was used to evaluate interaction effects, as described previously [13]. Measurements were performed on formalin-fixed, paraffin-embedded tissue Sects. (5 μm). For Aβ, the percentage of area covered by Aβ immunoreactivity was used, identifying both Aβ40 and Aβ42. The analysis focused on the percentage of the area covered by Tau2 immunoreactivity, identifying mature neurofibrillary tangles and dystrophic neurites. AT8, an antibody that recognizes phosphorylated tau (Ser202, Thr205), was also utilized. Interaction effects with Aβ were investigated using immunoreactivities of IBA1 and GFAP and gene expression levels of GPNMB. IBA1 was measured as the percentage of area covered by IBA1 immunoreactivity, GFAP as the percentage of area covered by GFAP immunoreactivity (identifying activated astrocytes), and GPNMB gene expression levels were quantified as fragments per kilobase of exon per million reads mapped (FPKM). The interaction effects were evaluated using the t-test, with values of p < 0.05 considered significant.

Trajectory analysis of human microglia snRNA-seq data

A trajectory analysis was conducted using human microglial snRNA-seq data and the Slingshot package. The default parameters for the slingshot function were utilized, and the getCurves function was used with the following parameters: approx_points = 300, threshold = 0.01, stretch = 0.8, allow.breaks = FALSE, and shrink = 0.99. From the dataset reported by Sun et al. [21], only microglia were extracted from immune cells for clustering. Clustering was performed using the Seurat package, following the tutorial provided by the Satija Lab (https://satijalab.org/seurat/articles/pbmc3k_tutorial). The settings were FindVariableFeatures with nfeatures = 4000, FindNeighbors, RunUMAP with dims = 1:10, and FindClusters with resolution = 0.5. To visualize gene expression over time, ggplot2 was used and regression curves were generated using generalized additive models (GAM).

Comparative analysis of GPNMB + Microglia with DAM and MGnD Subtypes

To compare GPNMB + microglia with DAM and MGnD subtypes, we retrieved their respective marker genes. The marker genes for GPNMB + microglia were identified from existing microglial subtype datasets using the FindAllMarkers function in the Seurat package, while those for DAM and MGnD were obtained from the relevant literature. The overlap between these gene sets was visualized using the UpSetR library. The Jaccard index, a metric of similarity between sets, was calculated using the following formula (3):

$$\begin{array}{c}J(A, B) = \frac{|A\cap B|}{|A\cup B|}\end{array}$$

(3)

where |A ∩ B| represents the number of common elements between sets A and B, and |A ∪ B| represents the total number of unique elements in both sets.

Heatmaps of CERAD and Braak scores in the ROSMAP and MSBB cohort data

To investigate whether GPNMB protein expression moderates the Aβ-tau interaction in patients with MCI and early AD in the ROSMAP cohort, correlations between CERAD and Braak scores were evaluated. The CERAD score in the ROSMAP cohort is a semi-quantitative measure of neuritic plaque density, classified into four scores ranging from 4 to 1, corresponding to no AD, possible AD, probable AD, and definite AD, respectively. The Braak score is a semi-quantitative measure of the severity of neurofibrillary tangle pathology, ranging from 0 to 6, corresponding to Braak stages I–VI. For the top 35% and bottom 35% of patients based on GPNMB protein expression, the frequency of each CERAD and Braak score was visualized using a heatmap. The same analysis was conducted using the proteomics data from MSBB cohort [25], which includes TMT proteomics data from 198 individuals diagnosed with AD and non-AD condition. Notably, the MSBB cohort did not include MCI patients and consisted only of AD patients and control groups. Therefore, our analysis focused on patients with CDR scores of 0.5 and 1 in both groups.

ADNI cohort analysis

Data used in the preparation of this analysis were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). The ADNI was launched in 2003 as a public–private partnership, led by Principal Investigator Michael W. Weiner, MD. The original goal of ADNI was to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). The current goals include validating biomarkers for clinical trials, improving the generalizability of ADNI data by increasing diversity in the participant cohort, and to provide data concerning the diagnosis and progression of Alzheimer’s disease to the scientific community. For up-to-dateinformation, see adni.loni.usc.edu.

CSF proteomics data were obtained using the SOMAscan platform, a proteomics technology that employs SOMAmer (Slow Off-rate Modified Aptamer) reagents. SOMAmers are synthetic DNA-based aptamers that bind to specific protein targets with high affinity and specificity, enabling the simultaneous quantification of thousands of proteins in a single assay. We utilized previously reported SOMAscan data that were normalized and quality-checked according to established protocols [46]. For data on Aβ and tau protein levels in CSF, we used Roche Elecsys Immunoassay data [47]. Among the samples that included both SOMAscan proteomics and Immunoassay data, we selected patients diagnosed with MCI (n = 112) and patients diagnosed with mild dementia (n = 132) for the analysis.

To investigate the relationship between GPNMB expression and Aβ42, tau, and phosphorylated tau levels, patients were stratified into two groups: the top 35% and bottom 35% based on GPNMB expression levels from SOMAscan proteomics. Statistical significance for differences in Aβ42, tau, and phosphorylated tau levels between the two groups was evaluated using a two-sample t-test.

Establishment of AlzPPMap

Using the R Shiny app, we developed a tool named AlzPPMap to analyze and visualize the integrated network. By inputting one or more gene names, the tool allows users to select a percentage of neighboring genes within the integrated network. The selected percentage of data is then used to construct and display a subnetwork. For instance, in this study, selecting HSPA5, APP, and MAPT, extracting the top 5% of genes closest to each within the integrated network, and executing the tool allows reproduction of the results shown in Fig. 4. In the “Table” tab, users can output the community and degree of each gene in the network. This information can be leveraged for enrichment analysis for GObp or other databases, facilitating identification of pathways associated with genes in each community. Additionally, since the degree value represents a gene’s connectivity level (hubness) within the network, genes with higher degree are more centrally positioned. Such genes may serve as potential drug targets for network-wide regulation, offering insights for experimental validation.

Results

Study design

The analysis workflow, including an overview of the datasets, methods, and results, is presented in Fig. 1. We extracted proteomics data for patients with MCI or early-stage AD, constructed a co-expression network, and integrated it with a PPI network using BIONIC, a deep learning-based method. To characterize the integrated data, we utilized enrichment analysis and whole snRNA-seq data (Fig. 2). Using proteomics data, we searched for factors associated with tau and Aβ (Fig. 3), constructed a subnetwork of tau, Aβ, and moderators, and performed cell type and functional enrichment analyses (Fig. 4). We validated these findings using immunohistochemical data from a separate cohort (Fig. 5). Additionally, for the subnetworks, we obtained microglia-specific snRNA-seq data to perform a microglial cell subtype enrichment analysis (Fig. 6). Finally, we validated these findings using immunohistochemical and gene expression data (Fig. 7).

Fig. 1

Study design. A schematic workflow of how the dataset was filtered and utilized in each figure is depicted. Data used for analyses are outlined with black borders, while analyses are framed in blue borders and filled in blue. Specific data applied for sub-objectives are highlighted by black frames within the blue sections. The analysis was conducted in two phases: exploration and validation

Fig. 2

Characterization of integrated protein network modules. a Schematic representation of the workflow for integrating proteomics and protein–protein interaction (PPI) data and evaluation of integrated features. b Stochastic neighbor embedding (tSNE) plot displaying protein modules. c Upper panels: Uniform Manifold Approximation and Projection (UMAP) visualizations of single-nucleus RNA sequencing (snRNA-seq) data for patients with Alzheimer's disease (AD) (GSE174367) [40] clustered according to cell types (left) and conditions (right) Lower panels: expression levels of proteins within the identified modules. d Representative Gene Ontology (GO) terms enriched in each protein module. e Network depicting the closeness among modules. Closeness was defined by averaging BIONIC-derived features within modules, calculating the Euclidean distances among modules, and taking reciprocals

Fig. 3

Inference of APP-MAPT modulators using a mutual information-based method. a Graphical representation of the Modulator Inference by Network Dynamics (MINDy) analysis. b Scatter plot illustrating the relationships between –log10(p-values) and Scores. Points representing p-values less than 0.05 are highlighted in red and are labeled with their corresponding gene symbols

Fig. 4

Glial cells moderate APP and MAPT interactions. a Workflow for network reconstruction from BIONIC analysis results, subnetwork construction using proteins highly associated with intermediary factors, APP and MAPT, and community detection within the network, followed by cell type identification and enrichment analyses for each community. b GO enrichment analysis of the subnetwork, categorized by parent terms using the rrvgo [45] software. c Representative GO terms enriched in each community. d Dotplot displaying expression levels of proteins within each community

Fig. 5

Histopathological analysis validated the roles of microglia/astrocytes as Aβ-tau moderators in early Alzheimer’s disease. a Schematic of the analytical workflow for linear regression analyses. Samples were sourced from four brain regions: frontal white matter (FWM), hippocampus (HIP), parietal neocortex (PCx), and temporal neocortex (TCx). RNA sequencing and immunohistochemistry data were used [21, 22]. A linear regression was performed, repeated for GFAP and IBA1, and Aβ as predictors of tau interactions. b Significant results of linear regression analyses for each of the four distinct brain regions, with subjects categorized into no-dementia and dementia groups

Fig. 6

snRNA-seq Analysis Identifying Enrichment for GPNMB + Microglia in Community 1. a, c Dot plot showing the expression levels of genes aggregated by community in each microglial subtype. [21, 22]. b, d Min–max normalized combined score for community 1, calculated by multiplying the expression level by the proportion of expressing cells

Fig. 7

GPNMB + microglia as a possible moderator of the Aβ-tau interaction. a, b Heatmaps displaying the correlations between CERAD and Braak scores in samples categorized by the top (a) and bottom (b) 35% of GPNMB expression. c Significant results of linear regression analyses for each of the four distinct brain regions, with subjects categorized into no-dementia and dementia groups. d Schematic representation summarizing the findings supporting the roles of both astrocytes and GPNMB+ microglia in the Aβ-tau interaction

Integration of protein co-expression and physical PPI networks

We constructed an integrated framework as a basis for analyzing the effects of PPIs and their moderators in the AD postmortem brain. To this end, we integrated proteomic data and a physical PPI network to construct a functional-physical PPI model. We employed BIONIC [30], a deep learning method designed to integrate multiple biological networks while preserving their integrity (Fig. 2a). BIONIC uses a graph attention network [35] to integrate multiple networks, maintain the distinctive characteristics of each, and generate new features. A graph attention network is a neural network designed for graph-structured data. It assigns varying levels of importance to different nodes and edges, thereby improving the capacity to learn meaningful patterns from complex biological networks. It shows good performance [30] and has been used in analyses of gene interactions in cancer [48].

First, a functional protein interaction network was constructed based on protein abundance data. We compiled postmortem proteome profiles for patients with MCI and early-stage AD from the dorsolateral prefrontal cortex, which were available from the Religious Order Study (ROS) and Rush Memory and Aging Project (MAP) cohort (https://adknowledgeportal.synapse.org/). Early-stage AD was defined based on the clinical criteria outlined in the NIA-AA guidelines and lecanemab clinical trials [10, 31]. Subsequently, we generated a PPI network using the STRING database [32] (https://string-db.org/) and constructed subnetworks by prioritizing “high confidence” interactions. Using BIONIC, the integrated framework was embedded into a lower-dimensional vector space, and the embedded coordinates were obtained in a 512-dimensional space. Our downstream analysis used 512 integrated network features.

Characterizing integrated protein network modules

Next, we characterized the biological context captured by the integrated network by evaluating the cell types and biological functions associated with the protein modules. To investigate these modules, we built k-nearest neighbor graphs of proteins based on the integrated network to obtain clusters (see Methods and Figs. 2a, b). A k-nearest neighbor graph links each protein to its k most similar neighbors based on network features, thereby effectively clustering proteins with similar biological properties. The cell types and biological processes corresponding to each module were inferred using publicly available snRNA-seq data [38] and Gene Ontology terms. The modules in the integrated network contained cell type signatures (Fig. 2c) and were enriched in biological processes related to these cell types (Supplementary Table S1 and Fig. 2d). For example, module 1 was characterized by endothelial cells, microglia, and astrocytes, modules 5, 17, 22, and 27 were associated with oligodendrocytes, modules 9 and 13 with astrocytes, and module 11 with neurons (Fig. 2c). Consistent with the cell type results, in module 5, terms related to oligodendrocytes, such as “Ensheathment of neuron” and “Oligodendrocyte differentiation,” were identified, whereas module 11 was predominantly enriched in synapse-related terms (Fig. 2d).

Identification of APP and MAPT as adjacent modules in the integrated network

We explored the inter-modular relationships within the integrated network by specifically examining modules that included APP and MAPT, which serve as markers for Aβ and tau, respectively. The APP gene was located in module 11 and MAPT was located in module 12. Additionally, APOE, which plays a significant role in the accumulation and aggregation of Aβ, was assigned to module 1. We then constructed a network of module–module interactions (Fig. 2e, and Supplementary Fig. S2). Notably, these three modules were in close proximity in the network, with a strong correlation between the APP and MAPT modules (Pearson correlation coefficient = 0.7667, p < 0.001). Additionally, strong correlations were observed between the APOE module and both the APP (r = 0.6479, p < 0.001) and MAPT modules (r = 0.6289, p < 0.001). Furthermore, within the network, the MAPT module was closest to the APP module based on Euclidean distances; the APP module was second closest to the MAPT module. These module-to-module correlations in early-stage AD support the interactions between APOE, APP, and MAPT [49, 50].

Next, we evaluated whether similar findings could be obtained without network integration. Employing the STRING database with settings for high confidence and physical interactions alone, no interactions were detected between the APP and MAPT modules. Similarly, APOE did not interact with the other genes (Supplementary Fig. S3). When examining protein expression data separately, APP and MAPT did not show significant correlations (Pearson’s correlation coefficient = −0.0465, p = 0.61). These results underscore the importance of our integrated network approach in capturing the biological interactions between Aβ and tau, which are crucial for understanding early-stage AD development.

Stress response factors moderate Aβ-tau interactions

Next, we aimed to identify factors that could moderate the Aβ-tau interaction based on the integrated network. For this purpose, we prioritized proteins with strong moderating effects on APP-MAPT interactions against proteins close to the APP and MAPT modules. We utilized the computational framework of Modulator Inference by Network Dynamics (MINDy) [41] (Fig. 3a). Briefly, MINDy calculates mutual information for interaction pairs (such as APP and MAPT), conditioned by the expression of the modulator, which reflects the modulating effect. Mutual information quantifies the information shared between two variables―specifically, protein interactions in this context―and measures how much the expression of one protein provides information about the other under the influence of a modulator. In the present analysis, we extracted the neighboring proteins of the APP and MAPT modules in the integrated network (see Methods, "Modulators of the Aβ-tau interaction"). We then stratified the samples from the proteomics dataset into groups representing the top and bottom 35% of expression for each neighboring protein and evaluated differences in the magnitude of mutual information between APP and MAPT (Fig. 3a). We analyzed 312 neighboring proteins and identified the 11 most significant modulators (Fig. 3b). These candidates included factors involved in the cellular stress response (HSPD1 and HSPA5) [51], ubiquitin–proteasome system (USP48, KBTBD6, UBE2D4, and STAMBPL1) [52,53,54], and apoptosis and cellular responses (MAP3K5) [54]. These results suggest that proteostasis involving stress response genes is related to the Aβ-tau interaction.

To further identify proteins contributing to the Aβ-tau interaction from the moderator candidates, we investigated proteins known to physically interact with Aβ and/or tau sourced from the Amyloidome [42] and Tau Interactome [43]. Consequently, Heat Shock Protein Family A (Hsp70) Member 5 (HSPA5) was identified. HSPA5, a marker of ER stress associated with AD [55, 56], has been shown to interact in vitro with Aβ and tau, mitigating their toxicity [57,58,59,60,61]. Although in vivo studies have identified it as a candidate therapeutic target [62, 63], the detailed mechanisms of action and its relevance in the Aβ-tau interaction are not fully understood [64].

Glial cells are closely related to the APP-MAPT interaction in the integrated network

We further investigated whether stress responses involving HSPA5 are involved in Aβ-tau interactions and identified related cell types and biological pathways. To address this, we exploited the integrated network extensively and reconstructed HSPA5 and its neighboring subnetworks as well as APP and MAPT (Fig. 4a). Note that the subnetwork included other significant moderators, STAMBPL1 and MAP3K5, suggesting that these factors, in addition to HSPA5, could be associated with the Aβ-tau interaction. This subnetwork was characterized by functions in neuronal death, positive regulation of kinase activity, cellular response to chemical stress, synapse organization, and regulation of autophagy, confirming its relation to stress response pathways (Fig. 4b, see also Supplementary Table S2).

We classified the subnetwork proteins into six communities based on topology and then estimated the cell types and biological functions associated with each community (Figs. 4c, d, and Supplementary Table S3). We found that glial cells were the primary cell type related to the stress response subnetwork (Fig. 4c). Additionally, the combined score for each community, calculated by multiplying the expression level by the proportion of cells, showed enrichment for glial cells (Supplementary Fig. S4a). Similar results were obtained using snRNA-seq data from two different cohorts [21, 22] and different regions (Supplementary Figs. S4b, c). However, communities 2 and 6 were excluded from this analysis owing to the small number of genes included. In the APP-related Community 3, oligodendrocytes were enriched in protein folding, oxidative stress responses, neuronal death, and autophagy. Community 4, which contained MAPT and HSPA5, was characterized by microglia and excitatory neurons with enrichment in kinase activity. Other communities such as communities 1 and 5 were also enriched in glial cells (Figs. 4c, d). Microglia were detected as a feature of three of the four significantly enriched submodules (Fig. 4c).

Histopathological data analysis validated the roles of microglia and astrocytes in the Aβ-tau interaction in early AD

To further investigate the effect of Aβ-tau interaction, we examined publicly available histopathological (ACT study) data [45]. This dataset contains extensive multimodal histopathological information, including Aβ and tau pathology. The distributions of age, sex, CERAD score, and Braak score in groups with and without dementia are shown in Supplementary Fig. S5. Among these parameters, CERAD and Braak scores showed significant differences between the groups (Fisher’s exact test, p = 0.00237 and p = 0.000468, respectively). In the no-dementia group, the average CERAD score was 1.476 and the average Braak score was 3.421. Braak stages I and II indicate that neurofibrillary tangles are confined mainly to the entorhinal region of the brain, whereas stages III and IV involve the limbic regions, such as the hippocampus.

We analyzed the effect of the interaction between tau and Aβ using a linear modeling approach. Tau2, an antibody for a mature neurofibrillary tangles and dystrophic neurites, was used as the dependent variable, whereas Aβ, along with GFAP or IBA1 (activated astrocyte or microglia markers, respectively), was used as the independent variable. We investigated the interaction between GFAP or IBA1 and Aβ, stratifying the analysis by the dementia status and across four brain regions (Fig. 5a). For GFAP, a significant positive contribution to the interaction was observed in the hippocampus of both the dementia and no-dementia groups (HIP, β = 43.59, p = 0.022 and β = 24.852 and p = 0.001, respectively) and in the temporal neocortex of the no-dementia group (TCx, β = 13.873 and p < 0.001). In the case of IBA1, a positive contribution trend was observed in the HIP of the no-dementia group, though the relationship was not significant (β = 4.678 and p = 0.338), whereas a significant positive contribution was observed in the TCx (β = 4.352 and p < 0.001). Conversely, in the dementia group, a significant negative contribution was noted in the TCx (β = −11.124, p = 0.04; Fig. 5b and Supplementary Table S4). Analyses of AT8, an antibody that recognizes phospho-tau (Ser202 and Thr205) showed the same trend in the no-dementia group: GFAP (β = 33.912 and p = 0.235 in HIP; β = 62.732 and p < 0.001 in TCx) and IBA1 (β = 14.442 and p = 0.36 in HIP; β = 16.344 and p < 0.001 in TCx) (Supplementary Table S5). Of note, interactions were observed only in the HIP and TCx regions, which are vulnerable in patients with early-stage AD [65].

The average CERAD and Braak scores in the no-dementia group mentioned earlier suggest that, although not exhibiting dementia, this group exhibited mild progression of AD pathology (Supplementary Figs. S5c, d). Collectively, these results support the hypothesis that astrocytes and microglia moderate the Aβ-tau interaction in early-stage AD.

GPNMB + microglia function as moderators of the Aβ-tau interaction

We further examined microglial subtypes that contribute to Aβ-tau interactions. Microglia-overexpressing neuroinflammation-related processes can be assigned to subtypes that exhibit different cellular states under different AD conditions [66]. We evaluated whether a certain subtype of microglia moderates the Aβ-tau interaction in early-stage AD. We accessed recently reported snRNA-seq data [22] from frontal cortex biopsies of 52 patients with hydrocephalus (NPH). This cohort included samples with confirmed amyloid and tau pathologies, in which AD was in an asymptomatic stage. Of the 52 living biopsies, Aβ plaques were identified in 19 subjects, Aβ plaques and phosphorylated tau pathology in eight subjects, and no pathology in the remaining 25 subjects; for subjects with Aβ pathology, the subsequent onset of AD was determined in longitudinal studies. In total, 892,828 high-quality single-nucleus profiles were included in analyses. Here, we focused on the microglial cell populations.

We performed an enrichment analysis of each microglial subtype for the stress response-related subnetwork communities (communities 1–6) based on the integrated network. GPNMB + microglia were significantly enriched in community 1 (Fig. 6). The expression level and proportion of cells identified as the GPNMB + microglial subtype was high for the community 1 signature score (Fig. 6a). In addition, we compared the compositions of cell subtypes across different statuses in the same cohort. The frequency of GPNMB + microglia was significantly higher in Aβ + and Aβ + Tau + samples than in the control group (Supplementary Fig. S6 and Supplementary Table S6). The combined score for Community 1, calculated by multiplying the expression level and the proportion of cells with positive expression, further confirmed that the GPNMB + microglial subtype was the highest ranked among the examined subtypes (Fig. 6b). Importantly, the GPNMB + microglial subtype is expressed in the early stages of AD and expands as the disease progresses [22].

To further validate the involvement of GPNMB + microglia in the Aβ-tau interaction in early-stage AD, we used another large-scale snRNA-seq data set for microglia subtypes in AD [21], including transcriptome data for 194,000 single-nuclear microglia and multiple AD pathological phenotypes in 443 human subjects. The cell subtypes were annotated based on a previous study [22]. We found the same trend towards enrichment for GPNMB + in Module 1 (Figs. 6c, d). These results suggest that the stress response module moderating Aβ-tau interactions in early AD corresponds to GPNMB + microglia.

Aβ-tau stress response mediators are activated prior to GPMNB + microglial activation

We further examined the relationship between stress response-related mediators of the Aβ-tau interaction and GPNMB + microglia. Previous studies have reported that ER stress induces GPNMB expression in neurons but not in microglia [67]. We hypothesized that stress response mediators are involved in the transition of microglia from a homeostatic state to an activated GPNMB + state in the human brain. To test this hypothesis, trajectory analysis of the two snRNA-seq datasets used in the previous analyses was performed using Slingshot [68]. Briefly, in this analysis, cells are embedded in the manifold space based on expression similarity between cells. This trajectory reflects the cellular state and can be treated as a “pseudo-time” axis. Multiple branched trajectories were identified. For each trajectory, we identified the GPNMB + microglial coordinates on the pseudo-time trajectories (Supplementary Figs. S7a, c).

Next, we analyzed the points at which the candidate Aβ-tau-moderating factors related to the stress response predicted using MINDy were expressed in the pseudo-time analysis. In both datasets, the moderator proteins were highly expressed prior to GPNMB + microglia in all pseudo-time trajectories (Supplementary Figs. S7b, d). Similar results were obtained for HSPA5, a moderating factor based on the amyloidome (i.e., the collection of amyloid peptides and related proteins, peptides, and other molecules) [42] and tau interactome [43]. These results suggest that stress factors moderating the Aβ-tau interaction may act through the activation of GPNMB + microglia.

GPNMB + microglia share limited similarity to disease-associated microglia observed in mouse models

Previous studies have shown that GPNMB is expressed in both disease-associated microglia (DAM) and microglial neurodegenerative phenotype (MGnD) [66, 69]. To further investigate the relationship between these microglial subtypes and GPNMB + microglia, we conducted a comparative analysis. Our results revealed that the number of shared marker genes between these subtypes was limited (Fig. S8a). Additionally, we calculated the Jaccard index, a metric for quantifying similarity between sets, which supported these findings (Fig. S8b). Furthermore, a detailed comparison of marker gene distributions reported in the literature showed that while GPNMB, LPL, CD9, and LGALS3 exhibited similar expression patterns, LILRB4 expression was significantly lower in GPNMB + microglia (Fig. S8c).

GPMNB protein expression is correlated with the clinical severity of AD

We evaluated the hypothesis that GPMNB + microglia moderate Aβ-tau interactions in early AD. In particular, we examined whether GPMNB could serve as a biomarker for the clinical severity of AD.

For this purpose, we assessed brain tissue protein expression and clinical parameters from the ROSMAP cohort data using CERAD and Braak scores as indicators of the clinical severity of AD. We stratified the patients into two groups based on GPNMB protein expression: the top 35% and bottom 35%. The GPNMB high- and low-expression groups corresponded to the GPNMB + microglia-activated and non-activated groups, respectively. Within the ROSMAP cohort, the CERAD score is a semi-quantitative measure of neuritic plaques, where 1 indicates definite AD, 2 is probable AD, 3 is possible AD, and 4 represents no AD. The top 35% showed greater disease progression based on both CERAD and Braak scores than that of the bottom 35% (Figs. 7a, b). A replication study using the Mount Sinai Brain Bank (MSBB) cohort yielded similar results in the parahippocampal gyrus proteomics data (Figs. S9a, b). Furthermore, in the Alzheimer’s Disease Neuroimaging Initiative (ADNI) study cohort, we evaluated the relationship between GPNMB expression levels in cerebrospinal fluid (CSF), as measured by SOMAscan proteomics [46], and the protein levels of CSF Aβ42, Tau, and pTau181 [47]. While no significant differences were observed for Aβ42, Tau and pTau181 levels were significantly elevated in the group with higher GPNMB expression (Fig. S9c).

GPNMB shows significant Aβ-tau interaction effects in the preclinical stage

We further examined whether GPNMB moderates the effect of the Aβ-tau interaction in early-stage AD. To this end, we assessed the histopathological and AD diagnostic information from the ACT study as well as gene expression data. We tested interaction effects of GPNMB gene expression on Aβ and tau by linear regression analyses (Fig. 7c, Supplementary Tables S7 and S8).

There were significant positive interaction effects of GPNMB gene expression on Aβ-tau in the non-dementia group in the HIP and TCx brain regions (β = 0.089, p = 0.02 and β = 0.026, p < 0.001 in HIP and TCx, respectively). Notably, no significant interaction effects were observed at the same sites in the dementia group. Collectively, these findings suggest that GPNMB + microglia and astrocyte moderate the Aβ-tau interaction in the early stages of AD (Fig. 7d).

Establishment of an integrated network analysis tool, AlzPPMap

Finally, we developed a tool named AlzPPMap (Alzheimer's disease Proteomics PPI Integrated Network Map) for the analysis and visualization of integrated networks. This tool enables users to analyze the integrated functional network derived from proteomics data alongside physical networks derived from PPIs. Users can select one or more genes of interest and, as demonstrated in this study, choose a specific percentage of closely related genes within the network based on their proximity to the selected gene(s). The selected genes are then used to construct, visualize, and detect communities within a subnetwork. AlzPPMap facilitates the identification of molecules and molecular pathways that are both functionally and physically associated with the target genes. In this study, we employed AlzPPMap as a proof-of-concept, constructing a subnetwork centered on HSPA5, APP, and MAPT. This approach revealed that astrocytes and GPNMB + microglia moderate the Aβ-tau interaction through stress responses. The tool can be accessed at: https://igcore.cloud/GerOmics/AlzPPMap.

Discussion

The number of individuals with AD is expected to continue increasing [70], making the development of effective treatments a pressing issue. As mentioned in the Introduction, whereas Aβ and tau proteins are known to contribute to AD pathology, the details of their involvement remain unclear [1, 2]. The specific interactions between them have been reported [4,5,6,7], but not the underlying mechanisms. In this study, we employed a deep learning-based network integration analysis to investigate the mechanisms of the Aβ-tau interaction in patients with early-stage AD. The study was conducted in four steps: (1) network integration, (2) identification of moderators and subnetwork construction, (3) validation analysis using histopathological data, and (4) development of a user-friendly web application.

In the first phase, we constructed an integrated functional-physical interaction network where two different data types—protein co-expression and physical PPIs—were combined into a single network using BIONIC [30]. Interestingly, enrichment analysis and snRNA-seq data showed that the integrated network captured relevant cell types and pathways. We also confirmed that Aβ and tau were present in neighboring modules in the integrated network. Notably, neither co-expression information nor physical interaction data alone could capture the Aβ-tau interaction, emphasizing the need to evaluate functional and physical interactions simultaneously. Co-expression network analyses are widely used in AD research [23,24,25], but it is difficult to represent physical interactions, such as Aβ-tau interactions or the amyloidome. Therefore, an approach that integrates networks based on different types of biological information may be effective in elucidating the complex molecular pathogenesis of AD.

In the second phase, we used the MINDy [41] algorithm to search for factors that moderate the Aβ-tau interaction using proteomics expression data. We identified ER stress markers such as HSPA5, along with other stress response-related factors, including HSPD1, USP48, KBTBD6, UBE2D4, STAMBPL1, and MAP3K5, as potential moderators. Although ER stress responses have been frequently reported to be associated with AD [55,56,57,58,59,60,61,62,63], their detailed mechanisms remain unclear. In this study, our findings suggested that stress response pathways were associated with the Aβ-tau interaction, indicating that elucidating how these pathways influence the progression of AD could be a focus for future research. In particular, understanding how ER stress responses are involved in the Aβ-tau interaction and transmitted through glial cells could contribute to the development of novel therapeutic strategies. We next identified molecules, cell types, and pathways involved in the Aβ-tau interaction by constructing subnetworks from the factors identified in the second phase, extracting genes closely related to the two proteins within the integrated network. The results suggested that glial cells, particularly GPNMB + microglia, are involved in the Aβ-tau interaction.

In the third phase, we validated the findings using an independent dataset with histopathological data from the ACT study [45]. Through linear regression analysis, we examined the interaction between GFAP, IBA1, and GPNMB with Aβ in the hippocampus and temporal neocortex in the non-dementia group. Furthermore, in the ROSMAP, MSBB, and ADNI cohorts, variations in Aβ and tau pathology progression based on GPNMB expression levels not only confirmed the reproducibility of the GPNMB-mediated Aβ-tau interaction but also provided novel insights. In the ADNI cohort, analysis of CSF protein levels revealed that patients with high GPNMB expression also exhibited elevated CSF tau levels, suggesting the GPNMB’s potential as a biomarker for disease progression. On the other hand, results for Aβ levels did not align between CSF and brain tissue. This discrepancy may stem from the well-documented inverse correlation between CSF and brain Aβ levels in AD patients [71]. Significant positive interactions were detected in multiple cohorts, further suggesting that astrocytes and GPNMB + microglia are moderators of the Aβ-tau interaction in early-stage AD. GPNMB + microglia, a recently identified subtype, are elevated in the frontal white matter and cortex of patients with AD [72]. Additionally, GPNMB is expressed in activated microglia [73] and is associated with AD pathology [66]. Furthermore, the expression of GPNMB in microglia is correlated with recently reported subtypes, such as disease-associated microglia (DAM) [69] and microglial neurodegenerative phenotype (MGnD) [22], which are closely associated with plaques and neurodegeneration. However, their specific functions related to Aβ and tau have not been fully elucidated. Our comparative study between GPNMB + microglia, DAM, and MGnD demonstrated that the number of shared marker genes between these subtypes and GPNMB + microglia is limited. One potential reason for this finding is that DAM and MGnD were originally identified in studies using the 5XFAD and APP-PS1 mouse models, respectively. These mouse models fail to replicate tau pathology observed in human AD patients, which may account for the differences in gene expression profiles. Our results suggest that the distinct marker gene profile of GPNMB + microglia may indicate a unique role in AD pathogenesis, providing new insights into their potential as a therapeutic target. Modulating GPNMB + microglia may help suppress tau aggregation and, consequently, mitigate neuronal death and cognitive decline. However, the precise mechanisms by which GPNMB + microglia contribute to the Aβ-tau interaction remain unknown. To translate these findings into clinical applications, further studies are required to elucidate the underlying molecular and cellular pathways. The present study along with recent work by Vahid Gazestani et al. [26,27,28,29,30] ulilizing meta-analysis of snRNA-seq data has reported that the transcriptional states of GPNMB + microglia differ from those of microglial subtypes identified in mouse models. This discrepancy highlights the necessity of establishing experimental systems that accurately model the unique characteristics of GPNMB + microglia in humans. Additionally, our trajectory analysis using snRNA-seq data from microglia suggests that GPNMB + microglial activation may be influenced by stress response pathways. This finding raises the possibility that the HSPA5-mediated stress response contributes to GPNMB + microglial activation and the moderation of the Aβ-tau interaction. By developing such experimental systems and performing detailed mechanistic analyses, it may become possible to further elucidate how stress responses influence the activation of GPNMB + microglia and their involvement in the Aβ-tau interaction, potentially revealing actionable pathways and molecular targets for therapeutic intervention.

In summary, network integration analysis has been extensively reported in cancer research [30], but its use in AD research remains limited. Through this study, we were able to identify potential moderators of the Aβ-tau interaction by integrating proteomics and PPI data using BIONIC [22]. Future analyses using a framework like the one presented in this study, integrating other omics data, could lead to new insights into the mechanisms of AD and other disease areas. Additionally, in the last phase, we developed a tool named AlzPPMap to evaluate the integrated network constructed in this study. This tool enables the analysis and visualization of relationships between other genes, potentially providing new insights and ideas for experimental validation. However, as mentioned in the Limitations section below, AlzPPMap does not distinguish between different gene isoforms or post-translational modifications of proteins, so incorporating these factors would be necessary for a more detailed identification of molecular interactions. Moreover, although AlzPPMap considers physical molecular interactions, it does not focus exclusively on AD-associated PPIs, which may result in false positives.

This study had several limitations. First, except for the GPNMB + microglial data [73] and ADNI analyses, all data were derived from postmortem brains. This introduces potential confounding factors between postmortem effects and AD pathology. To validate the findings of GPNMB + microglia of this study in vivo, the development of PET probes is necessary. Notably, P2X7R and P2Y12R have been highlighted as tracers for pro-inflammatory and anti-inflammatory phenotypes, respectively [74], and GPNMB also has the potential to serve as a tracer for microglial subtypes. Additionally, due to the postmortem nature of the brain samples, it is possible that the degradation of specific proteins was accelerated, potentially influencing the structure of the Aβ-tau interaction network. In particular, this factor could have impacted the detection of neuron-associated interactions within the Aβ-tau network. Second, in the analyses summarized in Fig. 4, protein expression was aggregated by genes, not accounting for multiple isoforms or post-translational modifications, especially in APP and MAPT. For example, Aβ40 and Aβ42 are produced from APP through different processes, and Tau in MAPT undergoes various phosphorylation modifications [75], impacting the disease differently. Histopathological evaluations of Aβ and tau were used to validate the integrated protein network analysis results. Third, the GPNMB gene expression analysis was based on bulk RNA-seq data, which may not specifically reflect GPNMB + microglia. A detailed analysis of the role of GPNMB + microglia in AD requires further experimental studies. Fourth, although we applied a significance threshold of P < 0.05 with the Benjamini–Hochberg correction, there is ongoing discussion in the field regarding the potential benefits of using a more stringent threshold. Future studies with larger sample sizes could help to further support our findings. Nonetheless, our findings offer valuable new insights into the Aβ-tau interaction and highlight potential directions for the development of AD therapeutics.