Apoptosis-associated genetic mechanisms in the transition from rheumatoid arthritis to osteoporosis: A bioinformatics and functional analysis approach

Osteoporosis (OP) and rheumatoid arthritis (RA) are two prevalent, interrelated conditions with growing global health burdens. Osteoporosis, primarily defined by a loss of bone mass and density, is associated with increased fracture risk, particularly in aging populations.^1–4 The increasing longevity of populations worldwide has significantly heightened the medical and economic impacts of OP.^1,3 On the other hand, RA is an autoimmune condition characterized by chronic inflammation of the joints, leading to pain, swelling, and progressive joint destruction.^5–7 Affecting about 1% of the global population, RA not only reduces patients' quality of life but also increases morbidity and mortality.^6,7

In recent decades, advancements in the treatment of RA, such as biologic therapies and disease-modifying antirheumatic drugs (DMARDs), have significantly improved disease outcomes, slowing joint destruction and inflammation.^8,9 However, RA patients, especially those undergoing glucocorticoid therapy, are at a heightened risk for developing osteoporosis.¹⁰ This dual burden of disease is driven by shared mechanisms, including chronic inflammation and the impact of RA treatments on bone metabolism.¹¹ The accelerated bone loss in RA, particularly through mechanisms involving immune system activity, makes the disease a key contributor to secondary osteoporosis.¹² 

One crucial process contributing to both RA and OP progression is apoptosis or programmed cell death.¹³ Apoptosis plays a major role in maintaining the balance between bone formation and resorption, and its dysregulation can lead to both excessive bone degradation (as seen in osteoporosis) and joint damage (as observed in RA).^14–16 Specifically, immune cells infiltrating joint tissues in RA release pro-inflammatory cytokines that promote apoptosis of bone-forming cells, leading to reduced bone density and strength.¹⁶ Despite some understanding of these processes, the precise genetic and molecular pathways linking apoptosis to the progression from RA to OP remain poorly defined.

The goal of this study is to address this gap by identifying apoptosis-related genes that contribute to both RA and OP using comprehensive bioinformatics analysis. We aim to explore the underlying genetic mechanisms involved in apoptosis, immune infiltration, and the progression of RA to OP. This research may unveil potential therapeutic targets that could help to mitigate bone loss in RA patients.¹⁷ 

In this study, the identification of key genes was conducted using a multi-step approach to ensure that the genes selected were highly relevant to both rheumatoid arthritis (RA) and osteoporosis (OP). We began by collecting differentially expressed genes (DEGs) from the GEO datasets GSE12021 (RA) and GSE56814 (OP). Differential analysis of dataset GSE12021 identified 2660 genes with p-values less than 0.05. Furthermore, volcano plots [Fig. 1(a)] and heatmaps [Fig. 1(b)] of differentially expressed genes were constructed. In GSE12021, through the WGCNA analysis, modules with p-values less than 0.05 were screened, yielding 1571 module genes related to RA [Figs. 1(c) and 1(d)]. This network analysis allowed us to focus on gene co-expression patterns, highlighting genes that are likely to be functionally related and potentially significant in RA. The GSE56814 dataset for OP revealed 418 key genes through differential expression analysis.

The next step involved intersecting these gene sets to identify genes common to both RA and OP using a differential analysis of the OP dataset, and by using a selection criterion of p < 0.05, we identified 418 key genes. We then created a heatmap [Fig. 2(a)] and a volcano plot [Fig. 2(b)] of these differentially expressed genes. Finally, we intersected the differentially expressed genes from RA, genes related to RA's WGCNA modules, and the differentially expressed genes from OP [Fig. 2(c)], ultimately identifying 23 key genes for further analysis. This intersection allowed us to narrow down genes that are potentially implicated in both diseases. To further refine this list, we cross-referenced the intersected genes with an apoptosis gene set, focusing specifically on genes involved in programmed cell death, a process known to play a significant role in the pathogenesis of both RA and OP.

We further conducted GO [Figs. 3(a) and 3(b)] and KEGG [Figs. 3(c) and 3(d)] enrichment analyses on these 23 key genes: C20orf27, EIF1, CDCA4, ARID5A, BRD4, JUND, ADRBK1, MMP14, NR4A1, ATXN2L, IER2, NFKB2, MED25, DND1, MAVS, GOLGA2, PLK3, ZC3H7B, JUNB, TSC22D4, ACIN1, FSTL3, and CCL25. The GO enrichment analysis results indicated that these key genes are mainly enriched in pathways, such as organelle inheritance, Golgi inheritance, and ossification. The KEGG enrichment results showed that these key genes are primarily concentrated in pathways like osteoclast differentiation and the MAPK signaling pathway.

To further explore the mechanisms from RA to OP, we intersected these key genes with genes related to apoptosis to investigate the role of apoptosis-related genes in these two diseases. As shown in Fig. 4(a), we identified 10 key genes: BRD4, MMP14, NR4A1, ATXN2L, NFKB2, MAVS, GOLGA2, PLK3, JUNB, and ACIN1. We conducted a protein–protein interaction analysis to study these genes and their related proteins, as depicted in Fig. 4(b).

Initially, we used LASSO on the OP dataset to screen the aforementioned 10 key genes, from which we identified six crucial genes [Figs. 5(a) and 5(b)], namely,BRD4, MMP14, NR4A1, ATXN2L, NFKB2, and ACIN1. Subsequently, a random forest analysis was conducted, resulting in the identification of 10 key genes [Figs. 5(c) and 5(d)].

Furthermore, we refined the screening of these 10 key genes in the RA dataset. Through the LASSO analysis [Figs. 6(a) and 6(b)], we identified four key genes: MMP14, ATXN2L, MAVS, and JUNB. The random forest analysis yielded the top 10 key genes [Figs. 6(d) and 6(e)]. By intersecting these sets of genes, we identified two crucial genes: ATXN2L and MMP14 [Fig. 6(c)].

Furthermore, we conducted the GSEA analysis on these two key genes to explore their potential pathways. The GO analysis results for ATXN2L showed significant enrichment in GOBP_PLATELET_DEGRANULATION and GOCC_AZUROPHIL_GRANULE. The KEGG analysis results for ATXN2L indicated significant enrichment in KEGG_CITRATE_CYCLE_TCA_CYCLE and KEGG_GLUTATHIONE_METABOLISM [Figs. 7(a) and 7(b)]. For MMP14, the GO analysis results revealed significant enrichment in GOBP_DETECTION_OF_CHEMICAL_STIMULUS and GOBP_MYELOID_LEUKOCYTE_MEDIATED_IMMUNITY. The KEGG analysis results for MMP14 showed significant enrichment in KEGG_LYSOSOME and KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION [Figs. 7(c) and 7(d)].

We conducted a correlation analysis between the selected key genes and the immune infiltration results from GSE56814 (Fig. 8). We found that ATXN2L is positively correlated with B cells naïve, T cells CD4 memory activated, and T cells regulatory Tregs [Fig. 8(a)], while MMP14 is positively correlated with B cells memory and negatively correlated with T cells CD4 memory resting [Fig. 8(b)]. All p-values were less than 0.05. Additionally, there were significant differences in T cells CD4 memory resting and Mast cells resting between the disease and normal groups [Fig. 8(c)].

Finally, we conducted ROC performance evaluation and found that ATXN2L has an AUC value of 0.667, and MMP14 has an AUC value of 0.649, indicating good diagnostic performance [Fig. 9(a)]. We also constructed a protein interaction network [Fig. 9(b)] and predicted the miRNA interaction networks associated with these two key genes [Figs. 9(c) and 9(d)].

Historically, managing osteoporosis (OP) has involved medication and lifestyle adjustments, with increasing public awareness about fracture prevention being crucial. Research indicates a significant correlation between osteoporosis and rheumatoid arthritis (RA), where RA-induced deterioration of bone quality accelerates the development of periarticular osteoporosis.¹⁸ Existing studies suggest that RA can exacerbate OP by accelerating the deterioration of cortical bone geometry and reducing bone mass, which in turn accelerates periarticular osteoporosis.¹⁹ Preventive measures for RA include a balanced diet, weight control, appropriate exercise, and correcting improper postures. Consuming foods rich in calcium and vitamin D, engaging in gentle exercises like walking and calisthenics, and correcting improper daily postures, such as prolonged squatting, are all effective prevention strategies. These measures, combined with targeted treatments informed by transcriptomic analyses, could offer more precise interventions for patients.²⁰ 

In this study, we employed a comprehensive bioinformatics approach to analyze microarray data from RA and OP using the GEO database datasets GSE12021 and GSE56814. Our differential expression analysis identified 418 key genes in the OP dataset and 1571 module genes related to RA. Heatmaps and volcano plots demonstrated the expression patterns of these genes in disease samples vs normal samples, underscoring the importance of cell apoptosis in OP and RA. By cross-analyzing the differentially expressed genes from RA and OP, we identified 23 key genes that potentially play critical roles in both diseases. GO and KEGG enrichment analyses of these genes revealed their primary involvement in biological processes, such as organelle inheritance, Golgi inheritance, and ossification, as well as pathways including osteoclast differentiation and the MAPK signaling pathway. These results provide deeper insights into the roles of these genes in both OP and RA, aligning with known mechanisms of bone metabolism and inflammation.

The application of machine learning methods, including LASSO and random forest analysis, refined the identification of key genes and confirmed their potential diagnostic value through ROC analysis. Immune infiltration analysis revealed significant correlations between key genes and various immune cell types, especially B cells, T cells, and regulatory T cells (Tregs). This finding highlights the importance of immune regulation in both OP and RA. Notably, ATXN2L (Ataxin-2-Like) and MMP14 (Matrix Metalloproteinase 14) emerged as genes of significant interest. ATXN2L, associated with RNA metabolism and neurodegenerative diseases, may contribute to the pathogenesis of OP through its involvement in cellular processes related to bone health.^21,22 MMP14, a key player in tissue remodeling and extracellular matrix degradation, is crucial for understanding the development of RA and its associated joint degradation.^23,24 Dysfunction in ATXN2L could potentially contribute to osteoporosis through disruptions in cellular processes, while altered MMP14 activity might impact the degradation of articular cartilage, influencing RA progression.²⁵ The identification of ATXN2L and MMP14 as key players in both RA and OP is clinically significant. ATXN2L, involved in RNA metabolism and linked to neurodegeneration, may also influence bone health through its role in post-transcriptional regulation, particularly in processes like apoptosis, which is integral to bone cell turnover. This makes ATXN2L a promising therapeutic target, especially in managing osteoporotic bone loss where apoptosis is dysregulated. MMP14, a member of the matrix metalloproteinase family, is essential for tissue remodeling and is implicated in extracellular matrix degradation. Its role in RA is particularly significant, as the breakdown of articular cartilage and joint tissue is a hallmark of RA progression. MMP14's role in cartilage degradation and tissue remodeling makes it a strong candidate for therapeutic intervention, with potential applications in both RA and OP. Modulating MMP14 activity could slow joint degradation in RA patients and reduce bone resorption in OP patients.

Our study provides a deeper understanding of how cell apoptosis influences RA and OP progression by integrating transcriptomic data with bioinformatics analysis. Our findings corroborate existing research showing that RA-induced inflammation accelerates bone loss and increases susceptibility to OP due to periarticular bone degradation. The identified pathways, such as osteoclast differentiation and MAPK signaling, are consistent with established mechanisms of bone metabolism and inflammation. The machine learning techniques that highlighted ATXN2L and MMP14 as key genes suggest that these genes could be promising targets for therapeutic interventions. Targeting these genes may help in developing more effective treatments for RA and OP by addressing the underlying molecular mechanisms. However, the reliance on public datasets introduces potential biases, and future research should validate these findings through experimental studies and expanded patient datasets. Overall, our study offers crucial insights into the molecular mechanisms of RA and OP and provides new targets for therapeutic strategies, potentially improving disease management and patient outcomes.

The ATXN2L and MMP14 genes may play important roles in the pathogenesis of osteoporosis and RA, offering new perspectives and potential therapeutic targets for future research. Further studies are needed to elucidate the specific roles of these genes and validate our findings through experimental methods. Our research provides valuable insights into the role of cell apoptosis in RA and osteoporosis and highlights potential new targets for treatment. By enhancing our understanding of these mechanisms, we can look forward to discovering new therapeutic approaches to more effectively manage these challenging conditions.

Microarray data for rheumatoid arthritis (RA) were obtained from the Gene Expression Omnibus (GEO) database, which is GSE12021 dataset,¹⁵ including 9 normal and 12 disease samples. For osteoporosis (OP), the GSE56814 dataset, containing 31 disease and 42 normal samples, was selected for analysis.

The Limma package in R was employed for differential gene expression analysis on the GSE12021 (RA) and GSE56814 (OP) datasets. The criterion for selecting differentially expressed genes (DEGs) was set to p < 0.05 to ensure significance.¹⁶ These DEGs were used for subsequent analysis to identify genes with altered expression between disease and normal samples.

The WGCNA package was employed for the WGCNA analysis of GSE12021, examining the correlation between modules and disease status. The parameters included power = 6, TOMType = “signed,” minModuleSize = 30, recreateThreshold = 0, and mergeCutHeight = 0.25. Genes selected for further analysis had a standard deviation greater than zero, excluding outliers. The most relevant modules to PMOP were identified, and the “pickSoftThreshold” function from the “WGCNA” package was used to set the optimal soft threshold, dividing the data into different modules, with b = 19 serving as the power for constructing an unscaled network. To merge similar modules in clusters, the threshold was set at 0.25, and the minimum number of modules was set to 30. Each module contained genes with similar co-expression characteristics.¹⁷ 

To identify key genes involved in both rheumatoid arthritis (RA) and osteoporosis (OP), we analyzed gene expression data from multiple sources. We intersected differentially expressed genes (DEGs) from the GSE12021 dataset (RA) with genes from RA-related Weighted Gene Co-expression Network Analysis (WGCNA) modules and compared these with DEGs from the GSE56814 dataset (OP). This intersection highlighted genes potentially involved in both diseases. To further refine our focus, we cross-referenced these genes with an apoptosis gene set to pinpoint those specifically related to apoptosis, thus identifying key genes implicated in the progression of both RA and OP through apoptotic pathways.

For functional enrichment analysis of the intersecting genes from rheumatoid arthritis (RA) and osteoporosis (OP), we utilized the “clusterProfiler” R package (v4.0). This analysis included Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment to identify significant biological processes, molecular functions, and pathways related to apoptosis and immune regulation. We clustered keywords based on similarity and selected the most representative terms with the highest enrichment. Statistical significance was determined with a threshold of p < 0.05.²⁶ 

Lasso regression, a machine learning method, was employed for apoptosis-related gene selection using the glmnet package, applied to both GSE12021 (RA) and GSE56814 (OP) datasets. Lasso helps to refine gene selection by applying a regularization penalty, effectively identifying key genes related to apoptosis, which minimizes overfitting and reduces irrelevant gene contributions.²⁷ The datasets were normalized, and cross-validation determined the optimal regularization parameter (lambda). The genes with non-zero coefficients were considered significant for further analysis. To validate their discriminatory power, Receiver Operating Characteristic (ROC) analysis was conducted, plotting the sensitivity and specificity to evaluate their predictive accuracy in distinguishing disease from normal states. The AUC (Area Under the Curve) values quantified the ability of each gene to differentiate between classes, with higher AUC values indicating better performance.²⁸ Following this, Single-Gene Gene Set Enrichment Analysis (GSEA) was performed to explore the functional relevance of these key genes. GSEA helps to reveal how each gene influences biological processes by identifying enriched pathways in the dataset. Statistically significant enrichment (p < 0.05) provided insight into the genes' involvement in pathways crucial to RA and OP pathophysiology, such as immune response, inflammation, and bone metabolism, further emphasizing their potential as diagnostic markers or therapeutic targets.²⁸ 

The CIBERSORT analysis technique was used to analyze the levels of immune cell infiltration between disease and normal samples in GSE56814, with the “PERM” parameter set to 1000 and the cutoff value at p < 0.05.²⁹ Additionally, the proportions of each immune cell type in the samples were calculated and displayed using bar graphs. The “pheatmap” package was used to create heatmaps of 22 immune cells, and the “vioplot” package was employed to display their abundance. Using the “corrplot” package, a correlation heatmap was created to visualize the relationships between 22 different infiltrating immune cells. Predictions were also made for miRNA regulatory genes related to key genes.³⁰ 

The authors have no conflicts to disclose.

Ethics approval is not required.

Hao-Ju Lo: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). Chun-Hao Tsai: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Tsan-Wen Huang: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.