Gene expression datasets and clinical phenotype data matching the search criteria “Alzheimer’s disease,” “human,” and “gene expression” were acquired through microarray dataset analysis. All gene expression profiles and corresponding platform probe annotations are publicly accessible for download from the Gene Expression Omnibus (GEO) database¹ (Zhu et al., 2020).

Using R software (version 4.3.2), we performed dataset-specific preprocessing for GSE5281, GSE29378, GSE37263, and GSE138260, which involved data reading and initial normalization using gene expression matrices and annotation files downloaded from the GEO database. After individual preprocessing, the datasets were merged to combine 134 normal samples and 142 Alzheimer’s disease (AD) samples, followed by batch effect correction and variance-stabilizing transformation. Differential gene screening was conducted using the “limma” package with empirical Bayesian analysis, applying significance thresholds of P < 0.05 and absolute log2 fold change (logFC) > 0.585 (Yang et al., 2020). The “pheatmap” package was utilized to generate visualizations, including volcano plots for differential expression analysis and heatmaps for clustering patterns of significant genes. Principal component analysis (PCA) was performed via the “prcomp” function to evaluate sample clustering, assess batch effect mitigation, and visualize key gene expression signatures distinguishing AD cases from healthy controls (Li et al., 2019; Vacchio et al., 2019). This integrated analytical pipeline ensured robust data normalization, rigorous statistical testing, and comprehensive visualization of molecular markers associated with Alzheimer’s disease.

Differential genes were analyzed by GO functional annotation and KEGG pathway enrichment using the “clusterProfiler” R software package (Yang-Chun et al., 2020), and the filtering criterion was set at P < 0.05 to understand the potential functional pathways and pathogenesis (Lu et al., 2020).

To identify genetic variants associated with gene expression, we conducted eQTL analysis using transcriptomic and genotypic data from multiple cohorts. Specifically, peripheral blood eQTL data comprising 5,311 European individuals were incorporated (Westra et al., 2013). The aggregated eQTL dataset utilized in this study was obtained from the GWAS Catalog website² (Cao et al., 2022). Employing the R package “TwoSampleMR,” we identified single-nucleotide polymorphisms (SNPs) with strong statistical associations (P < 5 × 10⁸) to serve as instrumental variables. Stringent linkage disequilibrium (LD) parameters were applied, setting the LD threshold at r² < 0.001 and defining an aggregation distance of 10,000 kb (Wootton et al., 2020). SNPs exhibiting weak trait associations or insufficient explanatory power for phenotypic variance were excluded through filtering based on an F-test value > 10 (Rosoff et al., 2021), ensuring only robust genetic instruments were retained for subsequent analyses.

The outcome data were sourced from the Genetic Association Database (see text footnote 2) within the GWAS Summary Dataset (IEU) (Wu et al., 2020). The specific GWAS identifier utilized was ebi-a-GCST90027158, which included 39,106 case samples and 46,828 control samples from European pedigree populations, encompassing a total of 20,921,626 single-nucleotide polymorphisms (SNPs). All GWAS summary statistics employed in this study are publicly accessible and available for free download.

Mendelian randomization (MR) analysis was conducted using the TwoSampleMR software package. To explore causal associations between Alzheimer’s disease and differentially expressed genes, we employed inverse variance weighting (IVW), MR-Egger, simple mode, weighted median, and weighted mode methods, complemented by sensitivity analyses (Chen et al., 2020). Co-expressed genes—including both upregulated and downregulated transcripts—were identified by intersecting disease-associated gene sets with differentially expressed gene lists. Subsequently, individual MR analyses were performed for each gene in this intersection to determine its causal relationship with Alzheimer’s disease. These analyses incorporated heterogeneity testing, multiple validity assessments, and leave-one-out sensitivity analysis to evaluate result robustness and reliability (Nie et al., 2020).

Immune cell infiltration profiles in the AD and normal tissue samples from the GEO AD dataset were quantified using the “LM22” signature matrix and the “CIBERSORT” algorithm in R software (Zhu et al., 2019). Statistical significance of differences in immune cell proportions between groups was evaluated with 1,000 permutations, and a P-value < 0.05 was set as the threshold for meaningful results. Visualization of immune cell infiltration patterns was accomplished by generating violin plots and heatmaps using the “pheatmap” and “ggplot2” packages, which clearly displayed the distribution of 22 immune cell subsets across samples. For immune correlation analysis, the Spearman correlation coefficient was calculated to assess the associations between infiltrating immune cell subsets using the “corrplot” package (Li J. et al., 2025). Meanwhile, the relationship between immune cell infiltration levels and the expression of immune checkpoint genes was explored through scatter plots and linear regression analysis with the “ggpubr” package, and statistical significance was determined by adjusting for multiple comparisons using the Benjamini-Hochberg method (FDR < 0.05) (Wang, 2025).

Single-gene GSEA enrichment analysis is a common method used to assess the enrichment of individual genes in a dataset (Bourdely et al., 2020). Instead of relying on differential genes, this method takes an enrichment perspective of the dataset by considering each gene in the expression matrix, ranking the genes according to a specific metric, and then checking whether the genes in the dataset are enriched at the top or bottom of the ranked list (Cai et al., 2019). Single-gene GSEA enrichment analysis is a powerful analytical tool that provides a more comprehensive assessment of the enrichment of all genes in a dataset, thus providing a deeper understanding of gene expression data (Sande-Melón et al., 2019). In this study, in order to more comprehensively explore the potential regulatory mechanism of each co-expressed gene in AD, we employed GSEA (Gene Set Enrichment Analysis) enrichment analysis and visualization in R, and selected “C2: KEGG gene sets” as the database (Li and Guo, 2020), and then performed single gene GSEA enrichment analysis for each co-expressed gene. P < 0.05 was considered as significant enrichment.

This study employed three machine learning algorithms—random forests (RF), least absolute shrinkage and selection operator (LASSO) logistic regression, and support vector machine-recursive feature elimination (SVM-RFE)—to screen for the characteristic genes of AD (Li G. et al., 2025). Specifically, the RF algorithm was implemented using the “randomForest” package in R software, LASSO logistic regression analysis was conducted via the “glmnet” package in R software, and the SVM-RFE algorithm was executed with the “e107” package in R software (Engebretsen and Bohlin, 2019). AD-related differentially expressed genes (DEGs) were obtained by taking the intersection of the characteristic genes identified by the RF, LASSO logistic regression, and SVM-RFE algorithms. Furthermore, the efficacy of these common AD-related DEGs in diagnosing AD was evaluated using the receiver operating characteristic (ROC) curve.

SK-N-SH cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM), which was supplemented with 10% fetal bovine serum (FBS) and 1% each of penicillin and streptomycin. All cell cultures were maintained in a humidified incubator at 37°C with 5% CO₂ (Goerges et al., 2024). To establish in vitro models that mimic distinct pathological features of AD, the SK-N-SH cells underwent two separate treatment regimens. First, the cells were incubated with 20 nmol/L okadaic acid for 48 h to simulate AD-related tau pathology (Boban et al., 2019). The second regimen involved exposure to 10 μmol/L amyloid-β 1-42 (Aβ_1–42) oligomers for 24 h to simulate AD-associated Aβ pathology (Han et al., 2017).

Total RNA was extracted from treated SK-N-SH cells using the TRIzol^® Plus RNA Purification Kit (Thermo Fisher, United States, REF:12183-555). The purity and concentration of the RNA was determined and then reverse transcribed to cDNA using the RNase-Free DNase Set (Qiagen, Shanghai, China, REF:79254). It was then processed using the Start-up reagent: SuperScriptIII First-Strand Synthesis SuperMix (Thermo Fisher, United States, REF:11752-050) and Power SYBR^® Green PCR Master Mix (Applied Biosystems, United States, REF:4367659). Finally, PCR was conducted on the CFX384 instrument (Bio-Rad, United States). The β-actin primer pairs was used as the internal control (Wang et al., 2020). The primer sequences used are shown in Table 1.

Four Alzheimer’s disease microarray datasets were retrieved from the GEO database as experimental datasets. The four datasets comprised 142 Alzheimer’s disease patients and 134 healthy controls in total. Details of the included datasets are provided in Table 2. Using R version 4.3.2, we performed normalization and integration of gene expression values across respective datasets and mitigated batch effects via principal component analysis (PCA). As illustrated in Figure 2A, pronounced batch effects were evident among the four Alzheimer’s disease gene datasets prior to correction. Following normalization and PCA-based batch effect adjustment, all samples within the integrated dataset exhibited satisfactory homogeneity, as demonstrated in the post-correction PCA analysis shown in Figure 2B.

In the analytical results, smaller P-values indicated stronger statistical significance for both gene sequencing consistency and differential gene expression. Overall, we identified 294 significantly upregulated and 330 significantly downregulated differentially expressed genes (DEGs). Supplementary Table 1 lists detailed annotations for these DEGs, including gene symbols, entrez IDs, and adjusted P-values. Figures 2C, D display the top 50 upregulated and top 50 downregulated DEGs, respectively, ranked by absolute fold-change values.

Through cross-tabulation analysis, we identified co-expressed genes from the intersection of related genes and differentially expressed genes, comprising five up-regulated genes (METTL7A, SERPINB6, VASP, ENTPD2, and CXCL1) and five down-regulated genes (FIBP, FUCA1, TARBP1, SORCS3, and DMXL2), as illustrated in Figures 3A, B. To further characterize the chromosomal localization of these genes, we generated a visualization of the co-expressed gene distribution across the genome (Figure 3C). Subsequently, we conducted a MR analysis on the 10 genes co-expressed with AD to evaluate the causal effects of each gene on the disease. The results indicated that all five upregulated co-expressed genes exhibited a significant positive causal association with AD in the MR analysis using the inverse-variance weighting method. Specifically, five upregulated co-expressed genes exhibited significant positive associations with Alzheimer’s disease: METTL7A (OR = 1.067; 95% CI: 1.026–1.110; P = 0.001), SERPINB6 (OR = 1.022; 95% CI: 1.002–1.043; P = 0.033), VASP (OR = 1.046; 95% CI: 1.002–1.092; P = 0.040), ENTPD2 (OR = 1.015; 95% CI: 1.015–1.099; P = 0.007), and CXCL1 (OR = 1.060; 95% CI: 1.019–1.104; P = 0.004). Conversely, all five downregulated co-expressed genes showed significant negative causal associations with the disease: FIBP (OR = 0.934; 95% CI: 0.897–0.973; P = 0.001), FUCA1 (OR = 0.943; 95% CI: 0.904–0.984; P = 0.007), TARBP1 (OR = 0.920; 95% CI: 0.853–0.993; P = 0.033), SORCS3 (OR = 0.909; 95% CI: 0.840–0.982; P = 0.016), and DMXL2 (OR = 0.950; 95% CI: 0.913–0.988; P = 0.011). Beyond the MR-Egger approach, additional validation analyses were conducted employing simple mode, weighted median, and weighted mode methodologies. For the five upregulated genes, all analytical methods consistently revealed an elevated risk of Alzheimer’s disease, as evidenced by odds ratios (ORs) greater than 1. Conversely, across all applied methods, the five downregulated genes consistently indicated a reduced risk of Alzheimer’s disease, with ORs consistently below 1 (Figure 4). The heterogeneity and pleiotropy tests for co-expressed genes yielded non-significant results (all P-values > 0.05), indicating no statistical evidence of heterogeneity or pleiotropic effects that would necessitate adjustment for these biases. Results from the leave-one-out sensitivity analysis demonstrated consistency between the effect estimates when each instrumental variable was excluded individually and the overall combined effect size, confirming the robustness of the analytical framework.

Sensitivity analyses were conducted on 10 key AD genes using MR-Egger regression and Cochran’s test. The results indicated no heterogeneity or pleiotropy, thus confirming the reliability of the findings (Table 3). The analysis of the funnel plot indicated that no individual single-nucleotide polymorphism (SNP) affected the outcome, implying the absence of directional pleiotropy for individual SNP non-violation and bias estimation. The leave-one-out analysis further confirmed the absence of horizontal pleiotropy, thereby demonstrating the robustness and reliability of the analytical methods and results (Figure 5). Furthermore, the present study examined the variations in expression of 10 pivotal genes in AD by utilizing the validation set GSE48350 dataset. The findings indicated notable distinctions in the levels of expression of these critical genes in AD, thereby validating their differential expression (Figure 6).

After the screening process, we successfully identified 624 genes associated with AD. To delve deeper into the potential functions of these differentially expressed genes, we performed GO and KEGG enrichment analyses. The GO enrichment analysis revealed that these genes were significantly enriched in biological processes, cellular components, and molecular function, including neuronal cell body organization, regulation of membrane potential, neuronal cell body, and passive transmembrane transporter protein activity (Figure 7A). In the KEGG pathway analysis, the differentially expressed genes were primarily enriched in Pathways of neurodegeneration-multiple diseases and the signaling pathways of Alzheimer’s disease (Figure 7B).

The CIBERSORT algorithm was employed to characterize immune cell profiles and investigate the association between Alzheimer’s disease co-expressed genes and immune cell infiltration. Figure 8A illustrates the distribution of 22 immune cell types across individual samples, depicting their proportional composition in each sample. We identified significant differences in specific immune cell subsets, specifically naive CD4⁺ T cells, between AD and healthy controls. Notably, the proportion of naive CD4⁺ T cell phenotypes was significantly elevated in AD samples relative to healthy controls (Figure 8B). Correlation analyses with 22 immune cell types (Figure 8C) revealed distinct associations for co-expressed genes: METTL7A exhibited positive correlations with naive B cells, resting memory CD4⁺ T cells, and M2 macrophages, while negatively correlating with memory B cells, plasma cells, CD8⁺ T cells, and follicular helper T cells. SERPINB6 showed negative associations with plasma cells and eosinophils. VASP was positively linked to naive B cells, resting natural killer (NK) cells, and M1 macrophages, but negatively associated with plasma cells and follicular helper T cells. ENTPD2 displayed a negative correlation with naive CD4⁺ T cell phenotypes. CXCL1 correlated positively with regulatory T cells (Tregs) and M1 macrophages, and negatively with resting mast cells. FUCA1 was positively associated with M2 macrophages and neutrophils, but negatively correlated with M0 macrophages. TARBP1 showed positive associations with memory B cells and plasma cells, while negatively correlating with Tregs and M1 macrophages. SORCS3 and DMXL2 both demonstrated positive correlations with plasma cells. Additionally, DMXL2 was positively associated with eosinophils and negatively correlated with Tregs and M1 macrophages.

To further explore the potential regulatory mechanisms of co-expression in AD, we performed single-gene GSEA enrichment analysis for each of the ten co-expression genes in the merged dataset of GSE5281, GSE29378, GSE37263, and GSE138260. We found that the expression of the ten co-expression genes was closely associated with multiple biological pathways. Examples: Cell adhesion molecules signaling pathway, Alanine, aspartate and glutamate metabolism signaling pathway, Alzheimer disease signaling pathway, Citrate cycle (TCA cycle) signaling pathway and so on. This again demonstrates that AD progression is a complex biological process and that the 10 co-expression genes may influence AD development by regulating different pathways. Among them, we noticed that several immune-related signaling pathways were significantly enriched (Figure 9). Therefore, we hypothesized that the expression of the co-expression gene may be closely associated with the immune response in AD.

The forest plot depicting the 10 AD-related DEGs is presented in Figure 3. Using the support vector machine (SVM) algorithm, we established that the model attained optimal accuracy with nine genes included (Figures 10A, B). We subsequently deployed the random forest (RF) algorithm to pinpoint potential diagnostic biomarkers (Figures 10C, D). Lastly, implementation of the least absolute shrinkage and selection operator (LASSO) regression algorithm generated nine candidate biomarkers, as depicted in Figures 10E, F. The nomogram indicated the importance of each gene in the diagnostic model (Figure 10G). The accuracy of the diagnostic model was evaluated using the calibration analysis, which showed high accuracy in diagnosing diseases, as demonstrated in Figures 10H, I. Furthermore, the area under the receiver operating characteristic curve (AUC) for the merged dataset (GSE5281, GSE29378, GSE37263, and GSE138260) was 0.860, indicative of robust diagnostic performance of the model for AD (Figure 10J). Finally, the intersection of genes identified by the SVM, RF, and LASSO regression analyses was visualized using a Venn diagram (Figure 10K). Nine common critical genes—METTL7A, SERPINB6, VASP, ENTPD2, FIBP, FUCA1, TARBP1, SORCS3, and DMXL2—were selected for final validation.

The specific expression levels of these nine common critical genes were compared between AD and control groups using the Wilcoxon rank sum test, with analyses performed on the merged dataset (GSE5281, GSE29378, GSE37263, and GSE138260) (Figure 11A). Nine critical genes exhibited statistically significant differences in the merged datasets. Receiver operating characteristic curves were then constructed to assess the diagnostic specificity and sensitivity of each gene in these datasets. In the merged dataset (Figure 11B), METTL7A (AUC = 0.740), SERPINB6 (AUC = 0.723), VASP (AUC = 0.723), ENTPD2 (AUC = 0.714), FIBP (AUC = 0.814), FUCA1 (AUC = 0.765), TARBP1 (AUC = 0.732), SORCS3 (AUC = 0.712), and DMXL2 (AUC = 0.696) all showed significant diagnostic value. In the GSE36980 dataset (Figure 11C), METTL7A (AUC = 0.622), SERPINB6 (AUC = 0.779), VASP (AUC = 0.653), ENTPD2 (AUC = 0.785), FIBP (AUC = 0.567), FUCA1 (AUC = 0.574), TARBP1 (AUC = 0.613), SORCS3 (AUC = 0.789), and DMXL2 (AUC = 0.700) exhibited diagnostic value.

Moreover, we validated the mRNA expression of METTL7A, SERPINB6, VASP, ENTPD2, FIBP, FUCA1, TARBP1, SORCS3, and DMXL2 in AD-associated tau and Aβ pathology model. The results revealed significantly increased mRNA levels of METTL7A, SERPINB6, VASP, and ENTPD2, whereas FIBP, FUCA1, TARBP1, SORCS3, and DMXL2 exhibited reduced mRNA expression in the AD-associated tau and Aβ pathology model (Figures 12, 13). Collectively, these findings indicate that all nine candidate genes could serve as potential diagnostic markers for AD, and may be involved in AD-associated tau and Aβ pathogenesis.