Analysis of human brain RNA-seq data reveals combined effects of 4 types of RNA modifications and 18 types of programmed cell death on Alzheimer’s disease

Abstract

Background

RNA modification plays a critical role in Alzheimer’s disease (AD) by modulating the expression and function of AD-related genes, thereby affecting AD occurrence and progression. Programmed cell death is closely related to neuronal death and associated with neuronal loss and cognitive function changes in AD. However, the mechanism of their joint action on AD remains unknown and requires further exploration.

Methods

We used the MSBB RNA-seq dataset to analyze the correlation between RNA modification, programmed cell death, and AD. We used combined studies of RNA modification and programmed cell death to distinguish subgroups of patients, and the results highlight the strong correlation between RNA modification-related programmed cell death and AD. A weighted gene co-expression network was constructed, and the pivotal roles of programmed cell death genes in key modules were identified. Finally, by combining unsupervised consensus clustering, gene co-expression networks, and machine learning algorithms, an RNA modification-related programmed cell death network was constructed, and the pivotal roles of programmed cell death genes in key modules were identified. An RNA modification-related programmed cell death risk score was calculated to predict the occurrence of AD.

Results

RPCD-related genes classified patients into subgroups with distinct clinical characteristics. Nineteen key genes were identified and an RPCD risk score was constructed based on the key genes. This score can be used for the diagnosis of AD and the assessment of disease progression in patients. The diagnostic efficacy of the RPCD risk score and the key genes was validated in the ROSMAP, GEO, and ADNI datasets.

Conclusion

This study uncovered that RNA modification-related PCD is of significance for AD progression and early prediction, providing insights from a new perspective for the study of disease mechanisms in AD.

Introduction

Alzheimer’s disease (AD), an incurable and progressive neurodegenerative disorder primarily affecting the elderly or pre-aged, is the most prevalent form of dementia, accounting for 60–70% of cases. It is characterized by progressive cognitive impairment and behavioral deficits [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25].

RNA modifications play crucial regulatory roles in various biological processes. Methylation modifications, such as m6A, m1A, m7G, and m5C, are common types of RNA chemical modifications. In neurodegenerative diseases like AD, these modifications modulate gene expression by affecting RNA stability, splicing, and translation efficiency, thereby influencing cellular functions and biological activities [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45]. For instance, m6A modifications have been shown to be associated with AD synapse plasticity. They affect synaptic regulation in patients by influencing the expression of the key plasticity protein CAMKII, which in turn impacts patient cognition [30, 31]. Additionally, changes in the levels of Tau protein and p-Tau are highly correlated with changes in key RNA modifiers [26, 44]. m5C modifications have also been studied in AD patients in recent years, and the coexistence of m6A and m5C has been demonstrated, highlighting the need to explore the role of m5C in AD [39, 46]. Although the role of m7G in AD is not well understood, its perturbation can induce cellular senescence and aging, which are key factors in AD [35, 43]. Dysregulation of m1A methylation is highly correlated with mitochondrial dysfunction, as it leads to translational inhibition of ND5 by catalyzing elevated levels of the TRMT10C protein, ultimately resulting in mitochondrial dysfunction [33].

Cell death, a fundamental physiological process accompanying cell growth, development, aging, and pathology, can be classified into non-programmed cell death and programmed cell death (PCD). PCD, a regulated form of cell death mediated by specific molecular mechanisms and influenced by genetic changes and pharmacological factors, plays important roles in neurodegenerative diseases like AD [47, 48]. By 2018, 12 types of PCD had been identified, and more recent studies have discovered novel modes, bringing the total to 18 types, including autosis, cuproptosis, anoikis, disulfidptosis, alkaliptosis, oxeiptosis, and mitotic cell death [49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64]. In AD, PCD is involved in neuronal death and neuroinflammatory processes. For example, ferroptosis has been demonstrated to have an impact on AD, with significant iron accumulation and ferroptotic features observed in AD neurons [49, 65,66,67]. Autophagy, another form of PCD, has a bidirectional role in AD. It is involved in the removal of abnormally accumulated proteins, long-lived proteins, and damaged organelles to maintain intracellular homeostasis and cell survival, but its dysregulation can also contribute to the pathogenesis of AD [68,69,70,71]. The role of other PCD modalities in AD warrants further investigation.

Previous studies have shown a close relationship between RNA modifications and PCD. In cancer research, m6A, the most common RNA modification, has been shown to be involved in the regulation of various forms of PCD, such as apoptosis, cellular pyroptosis, ferroptosis, necroptosis, and cellular autophagy [51, 52, 54, 72, 73]. M6A methylation can ameliorate pathological damage in hepatic fibrosis by modifying the post-transcriptional regulatory mechanism of hepatic stellate cell ferroptosis [74]. In addition to m6A, other RNA modifications like m7G, m1A, and m5C have also received increasing attention in cancer cell death studies. In recent years, the joint analysis of RNA modifications and forms of PCD has been applied in the study of some diseases [74, 75]. However, no studies have comprehensively analyzed the joint effects of RNA modifications and PCD in AD.

In this study, we explored the independent effects of RNA modification and PCD on AD. Then, by combining the two for comprehensive analysis, we screened for RNA modification-related programmed cell death (RPCD) using unsupervised consistency clustering to identify disease subgroups with RPCD features. We investigated the effects of RPCD on the synapses and immunity in patients with AD. The key role of PCD genes in AD was further demonstrated using gene co-expression network analysis. A combination of unsupervised uniform clustering, gene co-expression network analysis, and machine learning identified 19 key PCD genes (CNN2, CFI, GLP2R, C4B, SLC25A14, ANO6, PLEKHF1, GFAP, TP53, CDKN2B, ZWILCH, ME3, BCAR3, ELK3, OCA2, SEMA3A, COL12A1, TCEA3, and RWDD2A). We then introduced an RPCD risk score to predict disease onset and identify characteristics of high-risk populations. This score effectively predicted AD onset across multiple sets of data from multiple brain regions, thus establishing a strong association between RPCD and AD. The 19 genes used for the calculation also provided new potential targets for AD research, offering insights from a fresh perspective.

Four types of RNA modifications and 18 types of PCD affect AD occurrence and progression

We obtained parahippocampal gyrus (PHG) RNA-seq data from the MSBB study through the AD Knowledge Portal. Employing the Clinical Dementia Rating (CDR) scale for AD grouping, we selected 89 AD patients and 68 healthy controls as our experimental cohort. The PHG is a key site for AD pathological alterations, with these changes intricately linked to patients’ memory and cognitive functions. AD patients typically display reduced PHG volume and neuronal loss detectable microscopically [76, 77]. Thus, we initiated our investigation here before extending to other brain regions, as outlined in Fig. 1A.

Fig. 1

Correlation of 4 RNA modifications with AD. A For the experimental flowchart of this study, we analyzed the effects of 4 types of RNA modifications and 18 types of PCDs on AD, respectively. Subsequently, a joint analysis was performed to screen PCD genes related to RNA modifications to explore their potential roles in AD. Key genes were screened by unsupervised consensus clustering and WGCNA and risk scores were calculated. Patients were categorized into high-risk and low-risk groups based on risk scores to assess the impact of risk scores on disease progression. Finally, we assessed the diagnostic effect of risk scores on AD using ROC curves, which were validated using multiple external datasets (By Figdraw). B Box plots showed significant changes in the expression levels of four RNA modification regulators between the AD and Control groups. T-test was used to detect differential expression of regulators between AD and control groups. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. C Correlation analysis between the expression of RNA regulatory factors and CDR, Braak and PlaqueMean Volcano Plots. Horizontal coordinates indicate the magnitude of the correlation and vertical coordinates are − log10 (p-value). Red dots indicate positively correlated genes and blue dots indicate negatively correlated genes. D, E LASSO regression visualization of the MSBB RNA modifier cohort. The optimal λ is obtained when the partial likelihood of the deviation is minimized. F RNAM score and ROC analysis curves for the RNA modifier used to calculate the RNAM score

In light of the established upstream regulatory role of RNA modifications in neurodegeneration [78]. We scrutinized four common types: m6A, m7G, m1A, and m5C. Previous research has associated m6A with AD, while m7G, m1A, and m5C have also been explored in this context, suggesting potential involvements. We collated 67 regulators of these modifications from published literature and analyzed the MSBB PHG RNA-seq data (Table S1) [79,80,81].

We observed differential expression of different RNA modification factors between the AD group and the Control group (|LogFC|> 0, p < 0.05, Fig. 1B, Figure S1A). Many regulators exhibited significant correlations with AD-related clinical parameters, mostly negative ones (|r|≥ 0.2, p < 0.01, Fig. 1C). There were 27 CDR score-related regulators (10 positive and 17 negative correlations), suggesting an association of RNA modification with patients’ cognition and dementia level (Fig. 1C). A total of 34 (9 positive and 25 negative) correlations were associated with Braak stage-related regulators, suggesting a potential association with the severity of brain lesions in the patients (Fig. 1C). A total of 18 (six positive and 12 negative) correlations were associated with PlaqueMean (mean plaque density)-related Regulators, suggesting a relationship with inflammation in the pathological state of AD (Fig. 1C). m5C modifiers, such as TEM2, XB1, and USUN7, showed positive correlations with all three clinical correlates, whereas the other three RNA modifications were mostly negatively correlated, suggesting that they have potential roles in AD (Fig. 1C).

Subsequently, we assembled 39 potential AD-related modifiers and utilized LASSO—logistic regression with dichotomous variables to screen for those with the minimum λ and calculate coefficients (Fig. 1D and E). The RNA modification score (RNAM score) was then computed, and its predictive efficacy was evaluated. Boasting an AUC value of 0.802, the RNAM score demonstrated robust predictive potential for AD onset (Fig. 1F). However, when patients were stratified into high and low RNAM score groups based on the median, no significant disparities emerged in cognition, pathologic involvement, or inflammation (Figure S1C–S1J). Evidently, while RNA-modifying factors and the RNAM score revealed potential associations with AD, they failed to discriminate between patient groups based on disease characteristics.

Volume changes in the PHG are a striking pathological feature of patients with AD. Studies have shown that patients with AD have reduced gray matter volume in the hippocampus and PHG compared to healthy individuals, and we speculate that this is related to neuronal death. As various forms of PCD play a key role in neurodegenerative diseases such as AD, including ferroptosis, autophagy, apoptosis, and pyroptosis (Fig. 2A), we collated 11,206 PCD-related genes from 18 species from databases and literature collections (Table S2) [82,83,84]. After removing duplicates, 6366 genes were expressed in the experimental dataset. These factors were subsequently assessed. Similar to the RNAM score, a PCD score was constructed and had a good correlation and predictive effect on AD based on CDR score, Braak stage, and PlaqueMean (Figure S2A-C). Nevertheless, the CDR score alone was insufficient to effectively distinguish between high and low PCD score groups (Figure S2D and S2F).

Fig. 2

Patient clustering based on RPCD gene expression and co-expression network analysis. A The screening process for RNA-modification-associated PCD genes (By Figdraw). B Unsupervised consistent clustering based on RPCD Gene expression, when k = 2, MSBB parahippocampal gyrus cohort patients were divided into two Clusters. C Comparison of clinically relevant information between the two RPCD Cluster and Control groups, including age, CDR, Braak and PlaqueMean. T-test was used to test the significance of changes between two clusters. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. D Differential expression of PCD genes in two comparisons between RPCD Cluster A, RPCD Cluster B and Control groups. E Dot plot of top10 GO terms for differentially expressed genes between RPCD cluster A and B. F Enrichment analysis results of differential genes on SynGO between RPCD cluster A and B. G Heatmap of correlation between gene co-expression network modules and RPCD clusters with clinically relevant information. Numbers in parentheses are p-values and numbers outside parentheses indicate correlations. H Dot plot of top10 GO terms for green module hub gene. I Dot plot of top10 GO terms for the blue module hub gene

In summary, despite the potential impacts of both RNA modifications and PCD on AD, their inability to differentiate patient characteristics calls for further exploration of the effects of RNA modification-related PCD genes on AD, aiming to unveil possible mechanistic variations among patients.

RNA modification-related PCD genes are hub genes for changes in synaptic function, inflammatory state, and disease progression in patients with AD

To investigate the impacts of RNA modification-related PCD genes on AD, we employed Pearson correlations to screen PCD genes associated with RNA modification factors. A total of 3,173 genes meeting the criteria of |r|≥ 0.6 and p < 0.01 were identified as potential targets of RNA modification genes (Fig. 2A). Their regulatory relationships were then queried using the RM2Target database [85]. We then applied unsupervised consistent clustering, using PAM classification and the PAC algorithm to determine the optimal number of clusters. Patients were categorized into two subgroups: RPCD Cluster A and RPCD Cluster B (Fig. 2B). Comparing the clinically relevant information between the Control and patient groups, RPCD Cluster A had more severe dementia (CDR score), more severe pathological involvement (Braak stage), and a more severe inflammatory state of the disease relative to the control and RPCD Cluster B (PlaqueMean, Fig. 2C). GO analysis revealed a large number of differentially expressed PCD genes between RPCD Cluster A and the other two groups, while only seven such genes differed between RPCD Cluster B and the control group (Fig. 2D).

Differentially expressed genes in RPCD Clusters A and B were mainly enriched for the positive regulation of cell death, positive regulation of programmed cell death, and inflammation-related pathways (Fig. 2E, Table S3). Synaptic Gene Ontology (SynGO) analysis highlighted the relationship between PCD and synapses, with PCD genes significantly enriched in synapse- and presynapse-related functions (Fig. 2F). Functional enrichment between the RPCD Cluster A and Control groups also showed similar results (Table S3). In addition, we performed overall-level differential gene expression analysis and functional enrichment analysis, which also reflected the enrichment of programmed death-related pathways and synapse-related pathways in RPCD Cluster A compared to the other two groups (Figure S3A, Table S4).

Subsequently, we explored the role of RPCD genes in co-expression networks. We used weighted gene co-expression network analysis (WGCNA) to build co-expression networks, calculated soft thresholds, and subsequently partitioned the genes into different modules (Fig. 2G, Figure S3B). The green and blue modules showed the highest negative and positive correlations with expression in Cluster A, CDR score, and Braak stage, which we considered to be the key modules (Fig. 2G). We further screened the hub genes in the modules according to Gene Significance (GS) > 0.25 and Module Membership (MM) > 0.6 (Figure S3C and S3D). GO analysis showed that the hub genes in the blue module were enriched in pathways related to cellular morphogenesis, neuronal differentiation, regulation of autophagy, and synaptic modulation, whereas the hub genes in the green module were characteristically enriched in synaptic-associated pathways (Fig. 2H, I, Tables S5 and S6). RPCD genes were highly correlated with AD onset and progression, and multiple PCD pathways were crucial in differentiating patient populations. This indicates that RPCD is pivotal in AD progression, helping us better understand the disease. PCD genes' key positions in the co-expression network underscore their importance in the disease process.

Establishment and assessment of the RNA modification-related PCD risk score

To determine whether RNA modification-related PCD genes are effective in predicting AD and to further screen for potential biomarkers of AD, we used the intersection of 1311 differentially expressed genes and 642 key modular genes in the RPCD Cluster to obtain 437 alternative genes for input into LASSO (Fig. 3A). The key genes were screened using LASSO, and 19 key genes were identified (Fig. 3B). The network diagram showed that they were involved in the regulation of PCD (Fig. 3C). In addition, we conducted preliminary assessments of the feature importance of these genes using both the random forest algorithm and the LightGBM algorithm. For the random forest, we calculated the cross-validation error rates for various gene combinations through ten rounds of cross-validation, selecting the combination with the lowest average, which comprised 119 genes. Subsequently, we ranked the feature importance of these genes based on MeanDecreaseAccuracy (Figures S4A and S4B). For LightGBM, we similarly conducted an overall ranking of 437 genes and then demonstrated the feature importance of the top 100 genes by calculating SHAP values (Figure S4C). Our 19 key genes contribute to the disease to varying extents in both algorithms. The calibration curve also indicates that the model exhibits good predictive performance (Figure S4D). As described in the Methods, we extracted the coefficients of each gene, calculated the RPCD risk score based on these genes, and performed a series of assessments to determine the role of this score in AD. First, the risk score was significantly higher in AD patients than in the Control group, and significant differences in risk scores were also observed between men and women. A correlation was also observed between the risk score and APOE genotype; patients carrying genotype 4 had a higher risk score relative to non-carriers (Fig. 3D). In addition, the risk scores showed a high correlation with inflammatory plaque density, CDR score, and Braak stage (Fig. 3E). Linear regression analysis showed significant correlations among inflammatory plaque density, age, CDR score, Braak stage, and risk score (Figure S5A-C). The risk score tended to increase as the patients' cognition and pathological progression worsened, which was reflected in a significant difference in risk scores between patients with different levels of dementia and those in the Braak stage (Fig. 3D, G). In addition, we compared the importance of the risk score and 19 key genes for AD using a random forest model, and the risk score was significantly more effective than the other genes, followed by ANO6, OCA2, CDKN2B, and GFAP (Figure S5D). In conclusion, our risk score constructed on the basis of RPCD showed a relatively close relationship with AD-related clinical information and has a potential non-negligible association with AD.

Fig. 3

Establishes and evaluates RPCD risk scores. A Veen plot: intersection of differentially expressed PCD genes and key module genes in RPCD clusters. B Visualization of LASSO regression model based on the intersected genes, the optimal λ is obtained when the partial likelihood of the deviation reaches the minimum value. C Delineation of the 19 key genes in PCD species after LASSO regression model screening. D Variation of RPCD risk score in different kinds of populations, including Control and AD, male and female, and four different kinds of APOE genotype populations. T-test was used to test the significance between Control and AD. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. E Linear regression analysis of the correlation between RPCD Risk score and PlaqueMean. F Changes in RPCD Risk score at different Braak stages. T-tests were used to test for changes in risk scores between patients with different Braak stages. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. G Variation of RPCD Risk score at different CDR scores. T-tests were used to test for changes in risk scores between patients with different CDR scores. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. H RPCD risk score and ROC analysis curves for 19 key genes. I Column line graph showing the importance of RPCD Risk score, Age, CDR and Braak on AD risk

RNA modification-related PCD risk genes have important diagnostic significance in AD

We sought to determine whether the risk scores and the RPCD genes comprising these scores could predict AD onset. To this end, we constructed ROC curves using logistic regression to evaluate the predictive capabilities of the risk score and the 19 RPCD genes in relation to AD development. The ROC curves indicated a strong predictive performance of both the risk score and the 19 key RPCD genes for AD, with an AUC value reaching 0.887. When tested individually, the 19 RPCD genes also yielded favorable outcomes, typically exhibiting AUC values greater than 0.7 (Fig. 3H). These findings suggest that the risk score and its constituent genes play a potentially crucial role in AD.

Subsequently, by leveraging column-line plots, we evaluated the AD risk prediction incorporating both the risk score and patient clinical information. The risk assessment based on the risk score surpassed that of age and patient dementia, trailing only the Braak stage in predictive significance (Fig. 3I). Employing the LightGBM algorithm, we assessed the importance of the 19 genes in AD, identifying OCA2, ANO6, CNN2, and PLEKHF1 as the four most influential ones (Fig. 4A). SHAP values were computed to discern whether these genes act as drivers or suppressors of the disease. Notably, most of the RPCD genes were found to drive AD progression, with OCA2, ANO6, GFAP, and CPL12A1 exerting the most substantial driving effects (Fig. 4B and C). These results underscore the potential impact of PCD on AD progression, highlighting the necessity of considering RNA modification-related PCD when investigating AD processes.

Fig. 4

The importance of RPCD risk scores in the assessment of AD and risk grouping. A Visualization of LightGBM analysis: assessment of the importance of 19 key genes in AD. B, C SHAP Value calculation and model evaluation. D Patients were ranked according to the RPCD risk score and categorized into high and low risk groups using the median (− 6.100465). E Comparison of relevant clinical information (Risk score, CDR, Braak and PlapueMean) between high and low-risk groups. T-tests were used to test for changes in clinical information between patients between high and low-risk groups. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. F Proportion of population distribution between High and Low Risk groups, including gender, CDR, Braak and APOE genotypes. G Sankey plots of high and low risk groups, RPCD Cluster, CDR, Braak, and gender, organized to observe patient distribution. H Box plots of differences in expression of 19 key genes between high and low risk groups. T-test was used to detect differential expression of 19 key genes high and low risk groups. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001

The RNA modification-related PCD risk score divides patients into two populations with different clinical and biological functional characteristics

Previous studies have shown that RNA modification, PCD, and RNA modification-related PCD are crucial in AD. Both RNA modification and PCD have demonstrated potential roles in AD pathogenesis and possess good predictive capabilities for AD onset. However, scores constructed based on RNA modifications and PCD do not accurately characterize patients, and a better way to differentiate patient characteristics is needed to better understand the progression of AD. Therefore, we explored whether RNA modification-related PCD could better differentiate patient subgroups. Using unsupervised consistent clustering based on RPCD genes, we successfully divided patients into two subgroups with distinct biological traits. Subsequently, we screened for key RPCD genes and constructed an RPCD risk score. Patients were differentiated into high- and low-risk groups based on the median risk score (Fig. 4D). Significant differences in clinical characteristics were observed between the high- and low-risk groups. The degree of dementia (CDR score), Braak stage, and mean inflammatory plaque density were significantly higher in patients in the high-risk group than in the low-risk group (Fig. 4E). Linear regression analysis indicated that the risk score was not related to the age of the AD patients, but had a significant positive correlation with CDR score, Braak stage, and PlaqueMean (Figure S6A). The proportion bars show that the proportion of female patients was higher in the high-risk group (71.11%) than in the low-risk group (56.82%). The proportion of patients with high CDR scores and Braak stages was also significantly higher in the high-risk group than in the low-risk group, indicating that patients in the high-risk group had higher levels of pathological progression and dementia status. The APOE genotypes of patients in the high- and low-risk groups also showed significant differences, with a significantly higher proportion of patients in the high-risk group than in the low-risk group carrying the type IV APOE gene (Fig. 4F). Our results illustrate that the risk score based on RPCD can divide patients into two subgroups with significant differences in clinical characteristics and that the risk score also plays an important role in predicting AD disease progression. The alluvial plot as a whole demonstrates the overall relationship between the risk group, RPCD Cluster, and patients' clinical information. Patients in RPCD Cluster A are mostly classified into the high-risk group, and they have higher CDR scores and Braak stages. Female patients were mostly classified as high risk (Fig. 4G). These results suggest that the risk score can be used as a potential AD-related biological indicator, which has an important association with the disease onset and progression of AD. This also illustrates the non-negligible role of RPCD in AD.

Subsequently, we assessed the RNA levels between the high- and low-risk groups. We examined the differential expression of 19 key genes between the Control and high- and low-risk groups; 11 genes showed significant upregulation between the Control, low-risk, and high-risk groups, and 8 genes were significantly downregulated (Fig. 4H). Overall differential expression analysis identified the ten most significant genes, with GFAP among them (Fig. 5A). The heatmap demonstrated the 20 genes with the highest LogFC values, including OCA2, LTF, and other genes identified earlier, as well as some of the long-stranded non-coding RNAs. GO analysis showed the enrichment of these genes in synapse-related biological pathways, such as regulation of chemosynaptic transmission, modulation of trans-synaptic signaling, vesicle-mediated transport in the synapse, and synaptic vesicle recycling, and significant enrichment of these genes in cellular components such as the synaptic membrane, the postsynaptic membrane, and the synaptic membrane (Table S7). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the genes were enriched in pathways such as neuroactive ligand-receptor interactions, calcium signaling pathways, GABAergic synapses, and synaptic vesicle cycling (Figure S6B, Table S8). Gene set enrichment analysis (GSEA) showed similar results, with synapse-associated biological process (BP) terms being upregulated in the high-risk group while inflammation-related pathways upregulated in the high-risk group. Neurodegenerative pathways, multiple diseases, oxidative phosphorylation, and glutamatergic synapses were significantly enriched in KEGG terms (Figure S6C and S6D, Table S9).

Fig. 5

Biological function differences between patients in the high and low risk groups. A Differential gene analysis volcano plot between high and low risk groups. B Heatmap of differential gene analysis between high and low risk groups. C Bubble plots of GO terms or differential genes between the high and low risk groups. D Heatmap of differential expression of cellular senescence-related genes between high-risk and low-risk groups. E Bubble plots of GO terms for differential genes related to Cell Senescence in the High and Low Risk Groups

Given the importance of cell senescence in AD, an aging-related disease, we examined its impact in our study. We compiled cell senescence-related genes from MsigDb and analyzed their differential expression between the high- and low-risk groups. A total of 164 genes showed differential expression, with 96 upregulated and 68 downregulated in the high-risk group. GO analysis demonstrated enrichment in functions such as cellular senescence regulation and neurogenesis (Fig. 5E, Table S10). GSEA, on the other hand, showed downregulation of synapse-related functions in the high-risk group (Figure S6E and S6F, Table S11). These results demonstrate a relationship between our RPCD risk group and cellular senescence, suggesting a potential association between the two. In conclusion, our results suggest a potential role for RPCD in AD, and the risk score can distinguish patient subgroups well, with patients in the high-risk group having more severe pathologies, cognitive and inflammatory states, and more risky genotypes. Therefore, we conclude that RPCD is a good indicator of disease progression in AD. In addition, we identified a potential association between RPCD and cellular senescence, providing new ideas for our study.

The RNA modification-associated PCD risk score has good diagnostic validity across multiple brain regions and multi-omics data

We undertook a comprehensive assessment to validate the efficacy of our constructed risk score. Initially, we evaluated the risk score in various brain regions of patients within the MSBB dataset. The risk score was computed in the frontal pole, inferior frontal gyrus, and temporomandibular regions based on key genes, followed by grouping patients according to the risk score. In these three brain regions, no significant difference in Braak stage was observed between the low- and high-risk groups (Fig. 6A). In the frontal pole and inferior frontal gyrus, the CDR score was significantly higher in the high-risk group than in the low-risk group, whereas there was no significant difference in cognition between patients relative to changes in the superior temporal gyrus (Fig. 6B). In the frontal pole, the high-risk group had a significantly higher mean density of inflammatory plaques than the low-risk group, with no significant differences between the other two brain regions (Fig. 6C). In all three brain regions, the proportion of patients carrying the APOE4 genotype was significantly higher in the high-risk group than in the low-risk group (Fig. 6D).

Fig. 6

Effect of RPCD Risk score on other brain regions in AD patients. A Box plot of differences in Braak stages in the frontal pole, inferior frontal gyrus, and superior temporal gyrus between high and low risk groups. T-tests were used to test for changes in Braak between patients between high and low risk groups. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. B Box plot of differences in CDR score in the frontal pole, inferior frontal gyrus, and superior temporal gyrus between high and Low risk groups. T-tests were used to test for changes in CDR between patients between high and low risk groups. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. C Box plot of differences in PlaqueMean in the frontal pole, inferior frontal gyrus, and superior temporal gyrus between high and Low risk groups. T-tests were used to test for changes in PlaqueMean between patients between high and low risk groups. p-value: * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001. D Proportion of people with different APOE genotypes in each group in the frontal pole region, inferior frontal gyrus, and superior temporal gyrus

Furthermore, we validated the predictive effect of our key genes and their constitutive risk scores for AD using multiple datasets. Validation was performed using the MSBB PHG proteomic data, where 10 corresponding proteins of 19 key genes were expressed in the proteomic data, and their individual proteins and the combination of the ten genes had high AUC values, showing good predictive effects (Fig. 7A). The risk score in RNA-seq data from other brain regions of the MSBB also showed good prediction (Fig. 7B). We also used external datasets for validation and collected integrated datasets from the GEO hippocampus, internal olfactory cortex, temporal cortex, and frontal cortex collated from the ALZDATA database; our risk score showed similarly good predictive results in these data (Fig. 7C). Similarly, good prediction results were achieved using the ROSMAP RNA array data (Fig. 7D). In summary, our results reveal a potential connection between RNA modification, PCD, and RPCD in AD, all of which show good predictive potential for AD onset.

Fig. 7

External dataset validation of the effect of RPCD Risk score on AD prediction. A ROC curves of RPCD Risk score and key proteins in MSBB parahippocampal gyrus proteomic data. B ROC curves of RPCD Risk score in MSBB frontal pole, inferior frontal gyrus and superior temporal gyrus RNA-Seq data. C RPCD Risk score ROC curves in GEO hippocampus, internal olfactory cortex, temporal cortex, and frontal cortex RNA data. D RPCD Risk score ROC curve in ROSMAP human brain RNA-Seq data. E RPCD Risk score ROC curve in ADNI cerebrospinal fluid proteomic data and in ADNI whole blood genomic data. F RPCD risk score and ROC curves of common AD biomarkers in MSBB hippocampal parahippocampal proteomic data. G, H, I Correlation analysis between the expression of 19 key genes and CDR, Braak and PlaqueMean Volcano Plots. Horizontal coordinates indicate the magnitude of the correlation and vertical coordinates are − log10 (p-value). Red dots indicate positively correlated genes and blue dots indicate negatively correlated genes

Subsequently, we assessed whether the RPCD risk score and the genes used to construct it have potential clinical applications. Examination of cerebrospinal fluid and blood can aid in the diagnosis of AD; therefore, we collated cerebrospinal fluid proteomic data and hematological genome-wide data from the ADNI database to assess the diagnostic efficacy of RNA-modification-associated PCD genes. In the cerebrospinal fluid proteomic data of patients at baseline, the expression of proteins corresponding to 10 key genes was detected, and the constructed risk score had an AUC value of 0.692 (Fig. 7E). The corresponding risk score in the hematology genome-wide data reached an AUC value of 0.73, with a good diagnostic effect (Fig. 7E). We have compiled ten biomarkers associated with AD from the literature and compared them with our model [86, 87]. The ROC curve demonstrates that our RPCD risk score maintains good predictive performance (Fig. 7F). Lastly, we assessed the correlation between the 19 key genes in our model and pertinent clinical indicators. These genes exhibited strong correlations with CDR, Braak, and plaqueMean. Notably, ANO6 and GFAP align with our previous machine learning findings, indicating a potential association with the onset and progression of the disease (Fig. 7G, I). Therefore, we conclude that RNA modification-related PCD genes hold potential clinical applications in AD, offering a new avenue for AD treatment and research. The above results suggest that RPCD positively influences AD progression and the classification of patient subgroups, serving as a potential biological indicator of AD.

Method

Data sources and processing

In this study, we used multiple datasets from MSBB, ROSMAP and GEO. The MSBB (Mount Sinai/JJ Peters VA Medical Center Brain Bankcohort) cohort consisted of RNA-seq data from four brain regions: parahippocampal gyrus, frontal pole, inferior frontal gyrus, and superior temporal gyrus. The MSBB was screened according to the CERAD subgroups, and we chose AD and normal individuals as comparisons to be studied. The parahippocampal gyrus data were studied as a discovery set, including a total of 157 samples (AD = 89, Control = 68). Quantitative TMT protein data from the corresponding brain regions included 136 samples (AD = 76, Control = 60) as a proteomic level validation. Frontal pole data included a total of 187 samples (AD = 111, Control = 76); inferior frontal gyrus data included 154 samples (AD = 90, Control = 64); and superior temporal gyrus data included 167 samples (AD = 102, Control = 65) were used for multi-brain region data to conduct the study. GEO data were obtained from the ALZDATA database integrating multiple brain region RNAseq data. ROSMAP (The Religious Orders Study and Memory and Aging Project) data were used for RNA array data (AD = 124, Control = 144) from CERAD grouping of 268. MSBB and ROSMAP data were downloaded from Synapse (https://www.synapse.org). ADNI cerebrospinal fluid data were selected from the CruchagaLab SomaScan7K baseline-phase proteomic data, with the stipulation that MMSE > 26, CDR = 0 in normal patients was used as Control group, and MMSE < 26, CDR > 0.5 in patients was used as AD control group. Finally, 141 columns of samples were identified as Control and 117 AD sample controls. ADNI haematological genomic data were downloaded from the ADNI page, clinical information cerebrospinal fluid data were matched for clinical information and samples without corresponding clinical data were excluded. For data processing, use the scale function and log2 function to standardize the data used.

Collection and collation of RNA modifier genes and PCD-related genes

RNA modification genes were collected and collated from literature collation and RM2Target database, which included a total of 65 modification factors for the four RNA modifications. We obtained 11,206 PCD-related genes from 18 species through databases and literature collections. After removing duplicates, 6366 genes were expressed in the experimental dataset.

Screening of RNA modification-related PCD genes

PCD genes associated with RNA modification were screened using Pearson correlation. PCD genes with a P < 0.01 and |r|≥ 0.6 were considered associated with RNA modification factors. The identified genes regulated by RNA modification factors using the RM2target database were finally identified as PCD genes associated with RNA modification for subsequent analyses. Pearson correlation analyses were calculated using the ‘cor’ function in the ‘stats’ R package, and the significance was calculated using the ‘corr.test’ function in the ‘psych’ R package.

Identification of differentially expressed genes

Identification of differentially expressed genes was analyzed using the ‘Limma’ R package, and genes were considered to be significantly differentially expressed with a p < 0.05. Genes were considered to be up-regulated using LogFC > 0 and down-regulated using LogFC < 0.

Correlation analysis between genes and CDR, Braak and plaqueMean

Spearman correlation analyses were conducted using the ‘cor’ function in the ‘Statistics’ R package to explore the associations between genes and MMSE scores, as well as Braak staging. For correlations between genes and plaqueMean, Pearson correlation coefficients were calculated using the same ‘cor’ function in the ‘Statistics’ R package. The significance of these correlation results was assessed using the ‘corr.test’ function from the ‘psych’ R package. A correlation with |r|≥ 0.2 was considered low, |r|≥ 0.6 was regarded as moderate, and a P value of < 0.01 was deemed statistically significant, providing a quantitative measure of the strength and significance of the associations.

Screening and evaluation of key genes

To screen for potential predictors of Alzheimer's disease and construct gene scores, we implemented three state-of-the-art machine learning algorithms. The Least Absolute Shrinkage and Selection Operator (LASSO) adds a penalty term to the logistic regression objective. This term is the sum of the absolute values of the regression coefficients, multiplied by a tuning parameter λ. As λ increases, more coefficients are shrunk towards zero. This helps in reducing overfitting by selecting only the most relevant features. We determined the optimal λ by minimizing the binomial bias.

The gene score was calculated using the formula:

$\text{Gene score}=\sum_{\text{i}=1}^{\text{n}}\text{Coef}i*\text{X}i$ (Coefi represents the coefficient of LASSO, Xi is the expression of the corresponding gene or protein).

Random Forest Algorithm: The algorithm builds multiple decision trees using random data and feature subsets. This reduces tree correlation and overfitting. Predictions are aggregated from all trees (majority vote for classification, average for regression). We used the ‘randomForest’ R package to construct the model and cross-validation to identify genes important for AD prediction. LightGBM Algorithm: LightGBM is a gradient-boosting framework. It uses histogram-based algorithms and builds trees leaf-wise for faster training and less memory use. We ran the algorithm with the ‘lightgbm’ R package, calculated SHAP values to evaluate gene importance, and constructed ROC curves from linear logistic regression to assess predictive performance.

Linear regression analysis

Linear regression analysis was performed using the ‘glmnet’ function in the ‘glmnet’ R package to assess the relationship between gene scores and clinical information. The “ggplot2” R package was used to visualize the results of the linear regression analysis.

Identification of RNA modification-related PCD gene expression models in AD patients

We performed unsupervised consistent clustering of Alzheimer's patients based on RNA modification-associated PCD gene expression using the ‘ConsensusClusterPlus’ R package. The primary objective of this unsupervised clustering approach is to categorize patients into distinct groups by leveraging the similarity of their gene expression patterns, without any prior assumptions regarding the group structure. For this clustering task, we employed the PAM (Partitioning Around Medoids) classification, which is renowned for its robustness against outliers, thereby enhancing the reliability of the clustering outcomes. The PAC (Principal Axes Clustering) algorithm was used to find the optimal number of clusters (k). PAC calculates the proportion of sample pairs with consensus indices in (u1, u2) ∈ [0, 1] to infer the optimal k, with a lower PAC indicating a more stable clustering.

Weighted correlation network analysis (WGCNA)

A weighted correlation network was constructed using the ‘WGCNA’ R package, and the disease-related modules within the network were analyzed, with a particular focus on the two modules exhibiting the most significant phenotypic correlations. WGCNA constructs a weighted gene co-expression network by calculating the pairwise correlation between genes and then raising the correlation values to a power (soft threshold) to emphasize strong correlations. We determined the soft threshold based on the scale-free topology criterion, aiming to achieve an approximate scale-free network. A scale-free network is characterized by a power-law degree distribution, where a small number of genes (hubs) have a large number of connections, while most genes have relatively few connections. We selected a soft threshold that achieved an R² > 0.8 for the scale-free topology fit index. Subsequently, we partitioned the genes into different modules based on their topological overlap measure. Genes within a module share similar expression patterns and are likely to be functionally related. We screened for hub genes in key modules, defining them as genes with a gene significance (GS) greater than 0.25 and a module membership (MM) greater than 0.6. GS measures the correlation between a gene's expression and a particular trait, while MM quantifies the similarity between a gene's expression profile and the module eigengene. Hub genes identified in this way are likely to play crucial roles in the biological processes associated with the module.

GO and KEGG analysis

GO analysis was performed using the ‘clusterProfiler’ R package and the ‘richGO’ function. The ‘richKEGG’ function performs KEGG analysis.

Gene set enrichment analysis

Gene Set Enrichment Analysis (GSEA) is used to assess trends in gene distribution in gene sets in a gene table, which is ordered by phenotypic relevance, to determine the contribution of genes to the phenotype. The basic principle is to use a predetermined set of genes to rank genes according to their degree of differential expression in two samples, and then test whether a pre-selected set of genes is enriched at the top or at the bottom of the ranked table, where GO BP data is a pre-selected set of genes downloaded from the MisgDb database, which includes the GO BP terms and their corresponding sets of genes for GSEA analysis. GSEA enrichment analysis was performed using the GSEA function in the ‘clusterProfiler’ R package to assess changes in biological function.

Statistical analyses

All statistical analyses were performed using R version 4.4.1. Pearson correlation was used to assess the co-expression relationship between genes, while Spearman correlation was used to assess the correlation between genomic and pathological data, and the statistical significance of correlation analyses were as follows: *p < 0.01; **p < 0.001; ***p < 0.0001; ****p < 0.00001. 0.00001. Differences in gene expression between clusters were assessed using t-tests with the following statistical significance: *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.

Discussion

In this study, we first investigated the impacts of four RNA modifications and PCD in AD, revealing their potential associations with AD, as evidenced by significant correlations with patient dementia levels, pathological conditions, and inflammation. Both RNA modifications and PCD were potentially associated with AD, with significant correlations between them and the degree of dementia, pathological involvement, and inflammatory plaques in patients. We then combined RNA modification and PCD to explore the potential mechanisms underlying RNA modification-related PCD in AD. Both showed good predictive effects; however, they did not have the independent ability to effectively differentiate patients. Considering the close association between RNA modification and PCD that has been discovered in recent years, we considered whether combining the two could be used to better explore the mechanisms of AD onset and progression, and to this end, we screened the RPCD genes for subsequent exploration. RPCD genes were also highly correlated with the onset and progression of AD. We divided the patients into two subgroups based on RPCD expression to explore the RPCD modification patterns in patients with AD and found a significant difference between the subgroups, with RPCD Cluster A showing a more severe disease state. Comparing RPCD Cluster A and RPCD Cluster B separately against the Control, RPCD Cluster A had higher PCD gene expression. There was also a significant difference in the biological function between RPCD Clusters A and B. RPCD Clusters A and B were also significantly different from each other in terms of their biological functions. The construction of the gene co-expression network module helped us screen the key node genes of the network, which are mostly PCD genes that play an indispensable role in AD.

Finally, we combined the genes screened in the co-expression network module and the RPCD Cluster into LASSO to screen key genes and calculate the RPCD risk score. The risk score demonstrated an association with the disease and a good predictive effect. Unlike the RNA and PCD scores, the risk score can also be used to differentiate between patient phenotypes, with the high-risk group of patients highly overlapping with RPCD Cluster A patients with more severe disease states. We also validated the classification effect of the risk score using other MSBB brain region data and achieved better results. In addition, good disease prediction was observed using the MSBB proteomic data and other external datasets. Our results, for the first time, integrate the effects of multiple RNA modifications and PCD on AD. We provide a comprehensive assessment from multiple perspectives, which offers new insights into the exploration of AD.

Previous studies predominantly focused on either RNA modifications or PCD in isolation within the context of AD. For example, research on RNA modifications mainly centered around their influence on individual gene expressions or specific pathways related to AD. Studies on PCD, on the other hand, typically examined the roles of one or a few types of PCD, such as ferroptosis or autophagy, in the disease process. In contrast, our study innovatively integrated four types of RNA modifications and 18 forms of PCD. By combining unsupervised consensus clustering, gene co-expression network analysis, and machine learning algorithms, we constructed an RNA modification-related PCD network. This approach allowed us to identify novel connections and regulatory mechanisms between RNA modifications and PCD in AD, which were previously overlooked. Our integrated analysis not only provided a more comprehensive understanding of the complex pathogenesis of AD but also uncovered new potential therapeutic targets and diagnostic markers, thereby opening up new avenues for future research in this field.

RNA modifications play crucial roles in regulating cellular functions and gene expression. RNA modifications can significantly affect the stability and function of RNA molecules, and their interactions with other cellular components can be influenced by precise regulation of the modification site to control protein production. These modifications are essential for enhancing the stability and functionality of RNA molecules [28, 34, 38, 88]. In addition, RNA modification affects gene expression through cis- or trans-mechanisms, which are critical for cell function and cell fate decisions [27, 29, 30, 37, 43].

RNA modifications play important roles in the AD pathogenesis. In recent years, changes in the methylation level of m6A have been closely associated with AD synapse plasticity, which affects synaptic regulation in patients by influencing the expression of the key plasticity protein CAMKII, which in turn affects patient cognition [89,90,91]. In addition, the expression of some key modifiers has been shown to be highly correlated with AD pathology, and changes in the levels of Tau protein and p-Tau are highly correlated with changes in key modifiers [26, 31, 44]. m5C modifications have also been studied in patients with AD in recent years. In addition, the coexistence of m6A and m5C has been demonstrated; therefore, it is necessary to explore the role of m5C in AD [46]. The role of m7G in AD has not been well studied; however, its perturbation induces cellular senescence and aging. As aging-related diseases, senescence, and cellular senescence play key roles in AD, exploring m7G was also meaningful for our study [92,93,94]. Dysregulation of m1A methylation is highly correlated with mitochondrial dysfunction, and m1A methylation leads to translational inhibition of ND5 by catalyzing elevated levels of the TRMT10C protein, which ultimately leads to mitochondrial dysfunction. m1A methylation is highly correlated with mitochondrial dysfunction [33].

RNA modifications have several important effects on AD, including immunity, inflammation, mitochondrial function, and synaptic function. The changes in mitochondrial function, immune status, and neuronal cell death caused by RNA modification remind us of types of PCD such as ferroptosis, immunity, and autophagy-dependent cell death. This phenomenon gives us some inspiration that there is a close connection between RNA modification and PCD. The strong correlation between GPX4 and AD suggests that ferroptosis plays an important role in AD [95,96,97]. Autophagy plays a key role in a series of neurodegenerative diseases including AD. Autophagy is not only involved in the removal of abnormally accumulated proteins, long-lived proteins, and damaged organelles to maintain intracellular homeostasis and cell survival, but it is also closely related to the pathogenesis of AD. It has been shown that autophagy is closely related to the production and metabolism of both pathogenic β-amyloid (Aβ) and tau proteins and has a specific role in neuronal death [69, 98, 99]. There is also a close association between autophagy and immunity, which is prominent in microglia, where autophagy directly or indirectly influences the inflammatory response. Appropriate autophagic processes promote recovery from the disease, whereas autophagic disorders may lead to deterioration. Studies have shown that inhibition of autophagic processes in microglia and macrophages contributes to excessive neuroinflammation after brain injury [69, 100, 101]. The immune status of the patient is the point at which we must focus; the immune and acquired immune responses are the central drivers of neuronal death and brain atrophy. Therefore, immunoprototypic death is a disease-inducing mechanism that we must consider. One of the distinguishing features of AD is the progressive degeneration and death of neurons in the brain, leading to memory and cognitive decline. Considering PCD diversity of and the link between different types of PCD, we considered the combined effects of 18 types of PCD in AD and, for the first time, assessed the impact of multiple forms in AD.

Our study explored the effects of RNA modifications, PCD, and RNA modification-related PCD on AD separately, and assessed AD for the first time from multiple perspectives. We first observed the differential expression of the four RNA modifications between the Control and AD groups. m6A regulators such as METTL3, FTO, and YTHDF2 have significant roles in AD, which have been partially demonstrated. m6A regulators such as METTL3, as a pro-AD progression regulator, are involved in the modification of β-amyloid and affect tau protein pathology as well as cognitive level changes. Tau proteopathy and changes in cognitive levels were consistent with our results. Knockdown of METTL3 can improve the disease status of patients by decreasing the level of methylation [44, 89]. Our results showed significant downregulation of FTO in AD, and its downregulation also had a significant negative correlation with AD progression, which was significantly higher than that of most other RNA modification factors. FTO downregulation is a key factor in the progression of type 2 diabetes mellitus. In addition, FTO is a risk gene for type 2 diabetes, obesity, and other diseases, and these disease populations are also at a high risk for AD; hence, the role of FTO in AD deserves further investigation [102,103,104,105]. TET2 has a significantly higher positive phenotypic correlation with AD than any other RNA modifier, has been shown to be involved in immune-regulatory processes of AD, and is highly correlated with neuroinflammation [106]. It is involved in both DNA and RNA methylation, induces a microglial inflammatory response that promotes AD progression, and is a key regulator of AD.

Despite the lack of studies, some other RNA modifiers had more effects on AD than some of the proven modifiers in our study, especially the modifiers of m7G and m1A; thus, studies on the mechanism of RNA modification in AD have great potential. We further screened and integrated key RNA modification factors to construct the RNAM score, which helped us assess the combined effects of multiple RNA modification factors on AD. It showed good disease prediction and a significant correlation with the patients' pathologic phenotypes and genotype, suggesting that the combined role of RNA modifications in AD is critical and significant. However, we were unable to use this score in the AD subgroups. We found similar results during our study of the effects of PCD, with more than 50% of PCD genes differentially expressed between the Control and AD groups, and more than 200 PCD gene variants out of higher correlation with clinical phenotypes, all of which suggest a critical role of PCD in AD. Some of the key regulators, including co-regulators of ferroptosis (SLC7A11 and GPX4), co-regulators of cuproptosis (CDKN2A, SLC31A1, and DLD), and co-regulators of PCD in a variety of other cell types, were closely correlated with AD. These genes suggest a complex and unknown role for multiple PCD pathways in AD. With an increasing number of in-depth studies on forms of PCD, such as ferroptosis, cuproptosis, autophagy, apoptosis, pyroptosis, and immunogenic death, the link between cell death mechanisms in AD and AD pathology will gradually help us further our understanding of AD, and the study of programmed death is indispensable in the study of AD pathology. The PCD score calculated based on the PCD also showed a high correlation with the pathologic phenotype of AD, and the ROC curve also demonstrated its excellent predictive effect. However, differentiation between patient characteristics remains unclear.

We aimed to characterize the relationship between RNA modifications and PCD in the AD patient population. Therefore, a better method was needed to further our exploration. Several studies have shown a close relationship between RNA modification and different types of PCD, and the study of the mechanism of action of PCD under the effect of RNA modification has also achieved good results in cancer-related studies, whereas few studies have been carried out on AD. Therefore, we combined RNA modification and PCD in a follow-up study to achieve the expected results. We found a general correlation between RNA modifiers and PCD genes, with more than 50% of the PCD genes highly correlated with one or more RNA modifiers. The targeted relationship between them was also confirmed in the MR2Target database, which supports our hypothesis. Based on the expression of the RPCD genes, we divided the patients into two subgroups. A significant clinical phenotypic difference was observed between RPCD Cluster A and B, which achieved the desired effect. By further exploring the biological functional differences, we determined that Cluster A was characterized by significant PCD. Under differential gene comparisons, more than 1000 PCD genes were differentially expressed in the RPCD Cluster A and Control groups, whereas only seven PCD genes were expressed between the RPCD Cluster B and Control groups. Significant differences were also observed between Clusters A and B, and the difference in differentially expressed genes between these two groups was similar to the difference between the differentially expressed genes between Cluster A and Control, with more than 80% overlap. This highlights the PCD profile of Cluster A. Cluster A, as a subgroup characterized by PCD, exhibited more severe levels of deteriorated patient cognition, pathology, and inflammation, and the GO analyses provided a more intuitive picture of the PCD status of the two subgroups. The GO analyses showed us a variety of functional profiles of Cluster A compared to the other two subgroups, the most important of which included the response to injury, positive regulation of cell death, positive regulation of apoptotic processes, positive regulation of programmed cell death, and other biological pathways highly relevant to PCD. Next was the enrichment of multiple immune- and inflammation-related pathways, especially neuroinflammation-related pathways, which suggests a potential link between PCD and neuroinflammation, which may be highly related to immunogenic death. Functional enrichment analysis at the overall level helped us to understand the characteristics between different groups more comprehensively; the functional enrichment analysis between Clusters A and B, Cluster A, and Control all showed the enrichment of synapse-related pathways, and the regulation-related pathway of metal ions was also one of the important pathways, which again suggests a close connection between AD and PCD.

Co-expression network analysis further identified two modules that were positively and negatively associated with AD; more than 60% of the hub node genes were PCD genes, which played a key role in the co-expression network. Hub genes were enriched in numerous synapse-related functions, indicating that PCD plays a key role in synaptic modulation. Abnormal synaptic function causes cognitive decline and pathological progression in patients with AD. There are few studies on the role of PCD in synapses, but our study demonstrates the critical role of PCD in synaptic regulation and that dysregulation of PCD function may be closely related to the loss of synapses in AD. Studies have also shown the critical role of some forms of PCD in synapses, among which ferroptosis is currently the most widely studied. In the pathogenesis of Parkinson's disease (PD), increased iron in the substantia nigra (SN) region produces substances such as reactive oxygen species and reactive nitrogen species, which promote the formation of intracellular α-synuclein, thus causing degenerative ferroptosis in dopaminergic neurons [107]. The effects of autophagy on synaptic plasticity are well documented and can influence learning and memory formation in the same way [108, 109]. Other types of PCD and their mechanisms of action in AD also require further investigation.

Finally, we combined the RPCD cluster and gene co-expression network modules and screened 19 key genes as potential biological targets for AD using a machine learning algorithm. We further explored the relationships between the identified genes and clinical subgroups. In the RPCD Cluster A, patients exhibited more severe dementia (higher CDR scores), greater pathological involvement (higher Braak stages), and a more intense inflammatory state (higher PlaqueMean) compared to those in RPCD Cluster B and the Control group. Among the 19 key genes, several demonstrated notable associations with these clinical characteristics.

For instance, ANO6 and GFAP showed significant positive correlations with CDR, Braak stage, and PlaqueMean. In patients with more advanced dementia (higher CDR scores), the expression of ANO6 and GFAP was elevated. This suggests that these genes might contribute to the progression of cognitive decline and the associated pathological processes. In terms of Braak staging, higher levels of ANO6 and GFAP were found in patients with more severe neurofibrillary tangle burden, indicating their potential roles in the development of this key pathological feature. Additionally, the positive correlation with PlaqueMean implies that ANO6 and GFAP could be involved in the inflammatory response around plaques, potentially exacerbating the disease state. Conversely, genes like ME3 and RWDD2A showed negative correlations with these clinical indicators. In patients with better cognitive function (lower CDR scores) and less severe pathology (lower Braak stages), the expression of ME3 and RWDD2A tended to be higher. This indicates that these genes might play protective roles against the progression of AD, perhaps by inhibiting harmful processes such as neuroinflammation and neuronal death. When comparing the high-risk and low-risk groups divided by the RPCD risk score, significant differences in the expression of these 19 key genes were observed. In the high-risk group, which was characterized by more severe clinical manifestations, genes related to cell death and inflammation were generally upregulated. This further supports the notion that the identified genes are closely associated with the clinical progression of AD.

These 19 genes regulate a variety of PCD processes, of which immunogenic death is dominant; ferroptosis, autophagy, and autosis are also involved. TP53 is even involved in the regulation of nine types of PCD. TP53 is a common key regulatory gene in cancer research. Although less research has been carried out on its role in AD, but there are studies that have demonstrated the effect of its mutation on AD, and it is a potential therapeutic target for AD [110]. Based on the random forest and LightGBM assessments, we also focused on several genes. One of the genes we focused on was ANO6, which regulates astrocyte polarization by mediating the Notch signaling pathway. This process is closely related to the pathological process of AD, including a series of adverse effects such as phosphorylation of Tau protein caused by over-accumulation of Aβ protein. ANO6 is not only a potential biomarker, but also a new therapeutic target, which is important for a deeper understanding of AD pathogenesis and the development of new therapeutic approaches [111, 112]. GFAP may serve as a biomarker for the early stages of AD. GFAP may reflect AD-induced brain changes that occur before tau protein accumulation and measurable neuronal damage. Serum GFAP not only discriminates between AD and frontotemporal dementia (FTD) but also predicts the transition from mild cognitive impairment (MCI) to dementia. These findings emphasize the importance of GFAP in the diagnosis of neurodegenerative dementia [113,114,115,116]. PLEKHF1, a glucocorticoid receptor pathway gene, is closely associated with AD pathogenesis. Changes in the expression of PLEKHF1 may be associated with the neuroinflammatory and neurodegenerative processes of AD. As a risk factor for AD, PLEKHF1 can also affect the pathogenesis and developmental process of AD by affecting neuroinflammatory and neuroprotective mechanisms [117]. The influential roles of OCA2 and CDKN2B remain unknown. Subsequent explorations are needed to clarify whether C4B is also a gene to focus on, and research on it has been increasing in recent years. Studies have shown that the mechanism by which the APOE2 genotype reduces the risk of AD may be related to C4B and that PPP2CB has the second most significant single nucleotide polymorphism associated with carriers of the APOE2 genotype. PPP2CB expression correlates with the ratio of phosphorylated tau231/total tau. CD4 expression is associated with brain PPP2CB expression and indirectly influences AD disease progression [118,119,120,121]. Each of these genes is involved in the regulation of various forms of PCD, although the mechanisms of action of some of these genes in AD have not been addressed in the literature. However, their mechanisms of action in cells may indirectly affect the progression of AD and are potential biological targets of AD, which is important for a comprehensive understanding of AD pathology. Finally, we constructed the RPCD risk score that comprehensively evaluates the overall impact of key PCD genes on AD.

Considering the strong correlation between the RPCD risk score and AD-related clinical features, integrating this score into the clinical setting holds great promise for early AD diagnosis and patient sub-group stratification [122, 123]. For early AD diagnosis, the RPCD risk score could be incorporated into routine biomarker panels. Given its good predictive performance demonstrated by the high AUC values across multiple datasets, it could serve as an additional screening tool. In clinical practice, when patients present with mild cognitive impairments or other early symptoms suggestive of AD, healthcare providers could measure the expression levels of the genes contributing to the RPCD risk score in patient samples, such as cerebrospinal fluid or blood [123].

By calculating the risk score based on these measurements, clinicians can estimate the likelihood of a patient developing AD. For instance, a high RPCD risk score could prompt further in-depth examinations, including more comprehensive cognitive tests, imaging studies like PET scans to detect amyloid plaques and tau tangles, and longitudinal monitoring of the patient's condition. This approach could potentially lead to earlier detection of AD, allowing for timely intervention and treatment, which may slow down disease progression. In terms of patient sub-group stratification, the RPCD risk score can be used to categorize patients into different risk groups. This stratification could guide personalized treatment strategies. Patients in the high-risk group, characterized by more severe clinical manifestations and a higher likelihood of rapid disease progression, could be prioritized for more aggressive therapeutic interventions. These may include experimental treatments targeting specific PCD pathways or RNA modification enzymes that are implicated in AD pathogenesis.

On the other hand, patients in the low-risk group could be managed with less intensive monitoring and treatment regimens, allowing for a more tailored approach to patient care. Additionally, stratifying patients based on the RPCD risk score could aid in clinical trial design. By ensuring a balanced distribution of patients with different risk levels in clinical trials, researchers can better evaluate the efficacy of new treatments in different patient populations, leading to more accurate and reliable results.

Unlike the RNAM and PCD scores, the RPCD risk score can better differentiate the patients' clinically relevant phenotypes and APOE genotypes. GSEA visually demonstrated the functional changes between the high-risk and low-risk groups, with the immune- and inflammation-related GO BP and KEGG pathways significantly upregulated in high-risk patients and the synapse-related GO BP pathway significantly downregulated in high-risk patients. Therefore, we hypothesize that RNA modification-associated PCD has an overall impact on the immune status and synaptic function in AD and that activation of inflammation and loss of synapses may be responsible for severe pathological progression in high-risk patients.

We also validated the categorization effect of the risk score in other brain regions of the patients. In all three brain regions we assessed, we found significant changes in the degree of dementia between the subgroups of risk score patients. A significant difference in Braak stage was found in the frontal pole brain region. Surprisingly, we found genotypic variations between the patient groups in all brain regions. The percentage of APOE type 4 gene carriers was significantly higher in high-risk patients than in the low-risk and Control groups. Therefore, we hypothesize that the risk of AD caused by RPCD genes is highly correlated with the patient genotype.

Subsequently, we used multiple datasets to validate the predictive effect of the risk score in the MSBB protein data, MSBB multiregional data, GEO multiregional data, and ROSMAP whole-brain RNA microarray data, all of which achieved good predictive effects. Cerebrospinal fluid and peripheral blood are auxiliary indicators of AD that can help screen for early AD and provide better treatment to mitigate the progression of AD. To verify whether RPCD has a potential clinical application, we used cerebrospinal fluid and peripheral blood data from ADNI for validation, and the AUC value (0.7) reached the international diagnostic standard, indicating that our study showed potential clinical application value. We demonstrated a high correlation between RPCD and AD, which significantly affects the occurrence and disease progression of AD, and screened excellent biomarkers, presenting another perspective for the study of AD disease mechanisms. In addition, our study provides new insights for clinical treatment by predicting the occurrence and early treatment of AD through the regulation of programmed death. The demonstrated roles of RNA modifications and programmed cell death (PCD) in Alzheimer's disease (AD) pathogenesis suggest promising directions for therapeutic development.

For RNA modifications, targeting the enzymes involved in these processes could be a viable strategy. For instance, given the connection between m6A modifications and AD synapse plasticity, modulating the activity of m6A writers like METTL3, erasers such as FTO, and readers including YTHDF2 could potentially normalize gene expression and enhance synaptic function in AD patients. As for m5C modifications, identifying and regulating its modifiers might also influence AD progression, although their roles are yet to be fully understood. Regarding PCD, targeting specific pathways holds potential. Ferroptosis, which is implicated in AD, could be inhibited by agents like ferrostatin-1 to prevent neuronal death caused by excessive iron accumulation and lipid peroxidation. Autophagy, with its dual role in AD, requires precise modulation. Developing agents that promote autophagy in a controlled manner could help clear abnormal proteins and damaged organelles without causing harmful overactivation [122]. However, translating these ideas into clinical applications faces challenges. The complex interplay between RNA modifications and PCD pathways needs further elucidation, and developing effective drug delivery systems to the brain remains a hurdle. Future research should focus on clarifying these mechanisms and overcoming delivery issues.

Conclusions

In this study, we analyzed the brain data of patients with AD for the first time from three separate perspectives and elucidated the effects of RNA modification, PCD, and RNA modification-associated PCD on AD. A strong correlation between PCD and AD pathology was established by comparing these three perspectives and introducing the RPCD score to assess the patients' RPCD. The correlation between PCD and AD pathology was established and validated in several brain regions. Nineteen potential target genes were screened, and good predictive results were achieved using multiple datasets. Despite our findings, several limitations should be noted. Our study mainly relies on bioinformatics analyses, lacking direct experimental validation. The associations between RNA modifications, PCD, and AD were deduced from correlations and database searches. Without in vitro and in vivo experiments, establishing causal relationships is difficult. For example, while we identified correlations between certain factors and AD, direct experimental evidence is needed to confirm causation. The datasets used might introduce biases. The MSBB, ROSMAP, and GEO datasets, though widely used, might not represent the entire AD patient population. Samples were collected from specific regions or cohorts, potentially leading to selection bias. Variations in data collection and sample preparation across datasets could also affect result accuracy and comparability. Future studies should include more diverse samples to improve generalizability. The generalizability of our findings across different populations is uncertain. AD is heterogeneous, and genetic, environmental, and lifestyle factors vary among populations. The identified key genes and risk scores might have different predictive values in different ethnic groups. Thus, validation in diverse cohorts is essential. Finally, many datasets lack complete clinical information, limiting our analyses. In this study, we couldn’t perform certain analyses due to insufficient data. Comprehensive clinical data, including medical history, lifestyle factors, and biomarkers, is crucial for a better understanding of the disease. Future studies should prioritize collecting such data.