The GSE106241 data file was downloaded from the NCBI Gene Expression Omnibus public database (GEO, https://www.ncbi.nlm.nih.gov/geo/) annotated by GPL24170 as a Series Matrix File. The dataset included data on gene expression profiles from 60 human temporal cortical tissue samples with varying degrees of AD-related neurofibrillary pathology. The AD pathological features, such as Braak stages, α, β, and γ—secretase activity, and Aβ_1–42 levels of each sample were downloaded from the GEO database to perform Pearson’s correlation analysis with Braak stages-related immune hub genes (https://ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE106241). The AlzData database was a full collection of current high-throughput omics databases, such as genomics (GWAS and Whole Exome Sequencing), Proteomics, Functional genomics, and Transcriptomes data (Xu et al., 2018). In this study, we selected transcriptome expression data of temporal cortical from AlzData as a validation dataset, including 39 healthy controls and 52 patients with AD. In addition, the expression level of the Braak stages-related immune hub gene at the single-cell level in the brain was analyzed and visualized through the AlzData database.

Based on the genetic and clinical data in GSE106241, a weighted messenger RNA (mRNA) co-expression network was constructed using the WGCNA package in R. First of all, we used the gene expression spectrum to calculate the Median Absolute Deviation (MAD) of each gene, and removed the first 50% of the smallest MAD genes. Then, we used the goodSamplesGenes method of R software package WGCNA to remove outlier genes and samples. WGCNA was further used to build a scale-free co-expression network. After the acquisition of an appropriate power (β = 6), the adjacency matrix was transformed into the topological overlap matrix (TOM). Third, hierarchical clustering was performed to identify modules, and the eigengene was calculated. Finally, we calculated the correlation between Braak stages and each module through Pearson’s correlation analysis.

Functional enrichment was analyzed through the STRING online tool to investigate GO cellular components (CC), biological process (BP), molecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to potential immune-related pathogenesis of Braak stages in AD. The interaction score was 0.4. The false discovery rate (FDR) was < 0.05. The protein-protein interaction (PPI) network was constructed through the NetworkAnalyst based on the STRING interactome with a 900 confidence score (https://www.networkanalyst.ca/) (Xia et al., 2015).

Then, the total of 2,483 immune-related genes list was download from the Immunology Database and Analysis Portal (ImmPort) (https://www.immport.org/home) to screen out the DEIRGs (Bhattacharya et al., 2014). The cytoHubba plugin was adopted to screen out immune-related hub genes through three different algorithms [Edge Percolated Component (EPC), Maximum Neighborhood Component (MNC), and Degree]. In total, 91 temporal lobe transcriptomic data (39 healthy controls and 52 patients with AD) from the AlzData database were used as a validation dataset to analyze the difference in immune hub genes between the AD and HC groups.

In this study, CIBERSORT was used to assess the abundance of 22 types of immune cells in 60 samples with different Braak stages. CIBERSORT is an analytical algorithm that uses normalized gene expression profiles to assess the abundance of specific cells in complex tissues (Newman et al., 2015). After evaluating the abundance of 22 types of immune cells in each sample, we performed differential analysis and correlation analysis according to the Braak stages of the samples.

The Gene Set Enrichment Analysis (GSEA) was used to identify the different signal pathways between the high and low levels of immune hub genes in GSE106251. The annotated gene set c2.cp.kegg.v7.1.symbols.gmt was chosen as the reference gene list. The cut-off value for the GSEA was set as p < 0.05.

Statistical analysis and graphs were performed using Sangerbox online software (http://sangerbox.com/) and GraphPad Prism 5.0 software. A value of p less than 0.05 was considered statistically significant. Multiple testing corrections were made using the Bonferroni correction and Duncan’s multiple range test.

The WGCNA method was used to identify genes associated with Braak stages. First of all, we screened the top 50% highest variance of the expression profile (a total of 9,386 genes) from 60 samples for WGCNA analysis. Then, the scale-free network was constructed with a β value equal to 6 (R² = 0.74) (Figures 1A,B). Finally, a total of 14 co-expression modules were identified (Figure 1C). The connectivity was calculated and cluster analysis was performed among the 14 modules (Figure 1D). To further analyze the association between the models and phenotype, we calculated the correlation coefficients of each model with Braak stages. The results showed that the blue module (r = −0.32, p = 0.01) was the most negatively and the dark gray module (r = 0.31, p = 0.02) was the most positively associated with the Braak stage (Figures 1E–G). In total, 5,374 genes from these two modules were selected for the next analysis.

We intersected the above Braak stages-related immune genes screened by WGCNA with immune gene in the ImmPort database to screen out Braak stages-related immune genes. Among them, there were 260 Braak stages-related immune genes (Figure 2A and Supplementary Table 1). The enrichment analysis of GO cellular components revealed that these genes were mainly located at the extracellular region, cell surface, proteasome complex, MHC protein complex, etc. (Figure 2B). The biological processes of each of them were associated with signal transduction, cytokine-mediated signaling pathway, response to cytokine, cellular response to cytokine stimulus, etc. (Figure 2C). The enrichment analysis of GO molecular function showed that Braak stage-related immune genes were involved in signaling receptor binding, growth factor activity, cytokine activity, peptide antigen binding, etc. (Figure 2D). The KEGG enrichment analysis revealed that these genes were involved in AD, antigen processing and presentation, natural killer cell-mediated cytotoxicity, B-cell receptor signaling pathway, etc. (Figure 2E). The PPI network of the 260 Braak stages-related immune genes was constructed (Figures 2F,G). These results strongly suggest the role of immune-related genes in the pathological progression of Braak stages in AD.

The cytoHubba plugin was adopted to screen out Braak stages-related immune hub genes through three different algorithms (EPC, MNC, and Degree). The top 15 hub genes, filtered by the EPC algorithm, were CD86, HSPA4, FOS, GRB2, PTPN11, Kirsten rat sarcoma viral oncogene homolog (KRAS), PTPRC, CXCL12, TLR2, Phosphoinositide-3-Kinase Regulatory Subunit 1 (PIK3R1), IGF1, JUN, HSP90AA1, STAT3, and SOCS3 (Figure 3A). The DEGREE screened out TLR2, FOS, GRB2, HSPA8, IGF1, PTPRC, STAT3, HSP90AA1, JUN, CD86, PIK3R1, PTPN11, KRAS, HSPA4, and SOD1 (Figure 3B). TLR2, FOS, GRB2, HSPA8, IGF1, PTPRC, STAT3, HSP90AA1, JUN, CD86, PIK3R1, PTPN11, KRAS, HSPA4, and SOD1were found out by the MNC (Figure 3C). Finally, the co-existing gene of the three algorithms was selected as the hub gene, such as CD86, HSPA4, FOS, GRB2, KRAS, PTPN11, PTPRC, TLR2, PIK3R1, IGF1, JUN, HSP90AA1, and STAT3 (Figure 3D). After that, we then validated these Braak stages-related immune hub genes by using the temporal lobe transcriptome data of 39 healthy controls and 52 patients with AD from the Alzdata database. The results showed that the expression levels of GRB2, HSP90AA1, HSPA4, IGF1, KARS, PIK3R1, and PTPN11 were significantly decreased in the AD group compared with the HC group (Figures 3G–J,L–N), and the expression levels of CD86, FOS, JUN, PTPRC, STAT3, and TLR2 were not statistically different between the AD and HC groups (Figures 3E,F,K,O–Q).

We further analyzed the correlation between Braak stages-related immune core genes and AD pathological features, such as Braak stages, α, β, γ-secretase activity, and Aβ_1–42 levels. The expression levels of GRB2 (r = −0.321, p = 0.002), HSP90AA1 (r = −0.359, p < 0.001), IGF1 (r = −0.191, p = 0.07), KRAS (r = −0.344, p < 0.001), PIK3R1 (r = −0.467, p < 0.001), and PTPN11 (r = −0.291, p = 0.005) were negatively correlated with the grade of Braak stages (Figures 4A–G). The expression level of STAT3 was positively correlated with Braak stages (r = 0.264, p = 0.012) (Figure 4H). In amyloidogenic APP processing, we found that the γ-secretase activity was negatively correlated with the expression level of KARS (r = −0.40, p = 0.002) and PIK3R1 (r = −0.35, p = 0.009). The β-secretase activity was negatively correlated with the expression level of GRB2 (r = −0.36, p = 0.006), KRAS (r = −0.63, p < 0.001), PIK3R1 (r = −0.49, p < 0.001), and was positively correlated with PTPN11 expression level (r = 0.44, p < 0.001). In addition, the Aβ_1–42 levels were negatively correlated with KRAS (r = −0.29, p = 0.029) and PIK3R1 (r = −0.31, p = 0.019). However, the expression level of KRAS (r = −0.34, p = 0.011) was negatively correlated with the α-secretase activity in non-amyloidogenic APP processing (Figure 4I). Accordingly, we found that KRAS and PIK3R1 were not only involved in the Tau pathologically related Braak stages, but also closely related to the regulation of Aβ pathology.

Subsequently, we estimated the abundance of 22 kinds of immune cell infiltration in the GSE106241 samples by CIBERSORT to explore the role of the peripheral immune system in the progression of Braak stages in AD. We found that the ImmuneScore, which represents the total level of immune cells infiltration, increased with the increase of Braak stages (Figures 5B,C). In immune cell subtype analysis, four kinds of immune cells were significantly different in different Braak stages (Figure 5A). Among them, the Braak stages were negatively correlated with the abundance of follicular helper T cells (r = −0.337, p = 0.008), activated NK cells (r = −0.226, p = 0.082), and eosinophils (r = −0.348, p = 0.008) (Figures 5E,G,H). The abundance of M2 macrophages was positively correlated with Braak stages (Figure 5F). In addition, the abundance of plasma cells was negatively correlated with Braak stages (Figure 5D). Therefore, we speculated that peripheral immune cells play an important role in the pathological process of AD.

To explore whether these two core genes were involved in regulating abnormal infiltration of peripheral immune cells, we conducted a correlation analysis between immune hub genes and differentially infiltrated immune cells (Figures 6A,B). The ImmuneScore was negatively correlated with GRB2 (r = −0.31, p = 0.014), HSP90AA1 (r = −0.44, p < 0.001), HSPA4 (r = −0.36, p = 0.005), KRAS (r = −0.66, p < 0.001), and PIK3R1 (r = −0.73, p < 0.001) expression level, and was positively correlated with the PTPN11 (r = 0.29, p = 0.025). The abundance of plasma cells was positively correlated with GRB2 (r = 0.43, p < 0.001), HSP90AA1 (r = 0.45, p < 0.001), HSPA4 (r = 0.32, p = 0.013), KRAS (r = 0.54, p < 0.001), and PIK3R1 (r = 45, p < 0.001) expression levels. The abundance of follicular helper T cells was positively correlated with GRB2 (r = 0.32, p = 0.011), HSP90AA1 (r = 0.27, p = 0.034), KRAS (r = 0.37, p = 0.004), and PIK3R1 (r = 43, p < 0.001) expression levels. The abundance of activated NK cells was positively correlated with KRAS (r = 0.34, p = 0.008) and PIK3R1 (r = 0.43, p = 0.002). The abundance of M2 macrophage was negatively correlated with HSPA4 (r = −0.27, p = 0.038), KRAS (r = −0.26, p = 0.045), and PIK3R1 (r = −0.45, p < 0.001). Among them, we found that KRAS and PI3KR1 were the two most important hub genes in Braak stages-related immune regulation (Figure 6B).

To explore the function of Braak stages-related immune hub genes, we analyzed the expression level at the single-cell level of brain tissue. The results showed that the expression level of KRAS was the highest in neurons, but was limited in the other five cell types (endothelial, astrocytes, microglia, oligodendrocytes, and oligodendrocyte precursor cell) (Supplementary Figure 1A). The expression level of PIK3R1 was relatively high in six kinds of cells, among which the expression level was highest in neurons and astrocytes (Supplementary Figure 1B). After that, we explored the potential molecular mechanisms of KRAS and PI3KR1 associated with Braak stages in AD through GSEA. The results showed that axon guidance, long-term potentiation, inositol phosphate metabolism, and GnRH signaling pathway were significantly enriched in groups with high KRAS expression (Figure 7A). In addition, the high expression of PI3KR1 is involved in AD and Ubiquitin mediated proteolysis (Figure 7B). The PI3KR1 low expression group was related to cytokine–cytokine receptor interaction and RNA polymerase (Figure 7B). These results suggested that both KRAS and PI3KR1 were involved in the pathway of AD-related pathological mechanisms.