RA: Rheumatoid Arthritis; SpA: Spondyloarthritis; DEG: Differentially Expressed Gene; LGBM: Light Gradient Boosting Machine; SFS: Sequential Feature Selection; RFECV: Recursive Feature Elimination with Cross-Validation; CNN: Convolutional Neural Network; SVM: Support Vector Machine; AUC: Area Under the Curve; F1: F1 Score; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes ; DO: Disease Ontology; UPR: Unfolded Protein Response;HDF5: Hierarchical Data Format; RNA-seq: RNA sequencing; SMOTE: Synthetic Minority Over-sampling Technique; SHAP: SHapley Additive exPlanations; ROI: Region of Interest; FC: Fold Change

Rheumatoid arthritis (RA) and Spondyloarthritis (SpA) are chronic inflammatory diseases characterized by joint inflammation and varying degrees of systemic involvement. RA primarily affects the synovium, often leading to joint destruction and severe disability, while SpA includes diseases such as ankylosing spondylitis and psoriatic arthritis, impacting both axial and peripheral joints [1]. These conditions represent significant public health challenges, affecting millions worldwide and contributing substantially to the global burden of disability [2-4].

The pathophysiology of RA and SpA is complex, involving mechanisms such as autoimmunity, cytokine dysregulation, and aberrant immune responses [5,6]. Recent definitions have highlighted the need for advanced diagnostic tools and more targeted therapeutic approaches. Conventional methods for studying RA and SpA have been limited by insufficient specificity and sensitivity, often leading to delayed or inaccurate diagnosis [7-9].

Recent advances in transcriptomic technologies and machine learning have significantly transformed our understanding of these diseases. Techniques such as RNA sequencing (RNA-seq) and spatial transcriptomics have offered deep insights into the gene expression patterns within affected tissues, allowing for the identification of disease-specific cellular and molecular signatures [10,11]. Particularly, spatial transcriptomics has provided a nuanced view of gene expression within the context of tissue architecture, which is critical for understanding the complex interactions at the site of inflammation [12-14].

RA and SpA are both heterogeneous inflammatory diseases and there are continued endeavors to decipher the molecular mechanisms underlying these diseases using advanced complex approaches such as spatial transcriptomics and machine learning. These approaches allow for assembling of the map of disease-related gene expression profiles and also for investigation of cellular heterogeneity particularly in tissues, which opens news avenue for biomarker discovery and disease modeling. These approaches enable the identification of disease-relevant gene expression profiles at single-cell resolution, especially within tissues leading to new opportunities for biomarker identification and disease understanding. Developments in recent years have focused on the combination of spatial Transcriptomics with supervisors, to identify potential biomarkers for autoimmune diseases. For instance, investigating gene expression profiles of synovial tissues from RA patients using spatial transcriptomics revealed the upregulated immune-related pathways within the affected tissues [15]. Similarly, Zheng et al. [16] assessed differential gene expression (DEGs) in SpA specimens using spatial transcriptomics to establish specific pathophysiological processes at a cellular level.

Support Vector Machines (SVM) and Convolutional Neural Networks (CNN) have been used to improve the classification and prediction of autoimmune disease status using the transcriptome. The application of machine learning models in RA has also helped in the use of omics data to predict the disease prognosis [17]. However, there are several drawbacks like small sample sizes that can affect its findings. Nascimben [18] expanded about this problem as pathological for machine learning, and discussed possible remedies, such as cross-validation and synthetic data generation, particularly in relation to autoimmune diseases. In addition, Zhou et al. [19] also performed a systematic review on the challenges and possibilities of employing machine learning approaches to biomarker confirmation step in autoimmune diseases and the need for large and precise validation datasets. The review mentioned that since the major limitation was the small sample sizes, the use of multi-omics data integration and statistical methods can minimize these challenges and improve the prediction accuracy of RA and SpA [20].

Feng et al. [21] identified the gene expression profile using the spatial transcriptomics method of synovial tissues in SpA and observed that immune cell infiltration has comparatively worsened the disease which is agreed with our study. Similarly, Shulman et al. [22] used machine learning models trained from spatial transcriptomic data to predict disease activity in SpA, despite the noted drawbacks and limited sample size in their clinical studies [22].

Moreover, the integration of machine learning with transcriptomic has brought forward new approaches for diagnosis and finding therapeutic targets. Machine learning algorithms are also capable of processing large amounts of data and finding patterns that are not easily discovered through other means. It has been highly effective in distinguishing between different types of arthritis and exploring both individual and general biomarkers to create individual treatment plans [23-26].

This study aims to leverage cutting-edge transcriptomic technologies and machine learning to explore novel diagnostic and prognostic biomarkers, as well as potential therapeutic targets for RA and SpA, ultimately aiming to enhance patient outcomes through more precise diagnostics and targeted therapies.

Data acquisition

The dataset comprises high-dimensional spatial transcriptomics data and high-resolution tissue images from 6 synovial biopsy samples: three from RA patients and three from SpA patients. These biopsies were obtained from knee or hip joints during orthopedic total replacement surgery. The data was originally collected by [10] and subsequently processed and made available by [27]. For each biopsy sample, the dataset provides normalized gene expression count data stored in Hierarchical Data Format (HDF5) files with the ‘.h5ad’ extension, containing a matrix of normalized expression values for each gene across multiple spatial locations within the tissue section. Additionally, high-resolution tissue images are provided, corresponding to the areas where spatial transcriptomics was performed. The dataset includes spatial coordinates for each spot where gene expression was measured, allowing the transcriptomic data to be mapped back to specific locations on the tissue image.

Data preprocessing

The preprocessing pipeline was designed to filter, annotate, and visualize spatial transcriptomics data against the backdrop of corresponding tissue sections. The primary purpose of this preprocessing was to label the samples as either infiltrated or non-infiltrated by immune cells, as this information was not originally present in the genetic dataset. This labeling is crucial for subsequent analyses of differential gene expression between infiltrated and non-infiltrated regions.

The input for this preprocessing step consists of normalized gene expression data in .h5ad format, high-resolution tissue images corresponding to the gene expression data, and spatial coordinates for each gene expression measurement spot. The goal is to label each sampling spot in the dataset as either ”infiltrated” or ”non-infiltrated” based on its location relative to regions of immune cell infiltration visible in the tissue images. The output is an annotated dataset where each gene expression measurement is labeled with its infiltration status, along with visualizations of these annotated regions.

The preprocessing steps were as follows:

1. Region of Interest (ROI) Detection: Tissue images were subjected to grayscale conversion and thresholding techniques to isolate specific regions annotated for analysis. Contours delineating these regions were identified, scaled, and subsequently used to generate masks that isolated ROIs within the tissue images. Pixel coordinates within these ROIs were extracted for annotation to samples.
2. Spatial Coordinate Transformation: Spatial coordinates from the transcriptomic data were converted to match the pixel coordinates of the tissue images, allowing for a direct comparison and overlay of gene expression data onto the tissue images. This transformation facilitated the identification of transcriptomic signals within the specified ROIs.
3. Immune Cell Infiltration Sample Identification: Utilizing the ROI-specific pixel coordinates, samples within the black line marked regions were identified as infiltrated. These samples were flagged within the dataset for subsequent analysis.
4. Data Annotation and Visualization: The final step involved the annotation of the processed spatial transcriptomics data with the infiltration status of each sample. Visualization tools were employed to generate plots displaying the spatial distribution of gene expression with overlays indicating ROIs and infiltrated regions.

Identification of differentially expressed genes

Differentially expressed genes (DEGs) between RA and SpA samples and control samples were identified using the ‘limma’ package in R [28,29]. The selection criteria were adjusted P value <0.05 and |log2-fold change (FC)| >0.6. Volcano plots and heatmaps of DEGs were created using the ‘ggplot2’ and ‘pheatmap’ packages in R, respectively.

GO, KEGG, and DO enrichment analysis

The R packages ‘clusterProfiler’ [30] and ‘DOSE’ [31] facilitated the analyses of Gene Ontology (GO), including Biological Process, Molecular Function, and Cellular Component, as well as KEGG pathway and Disease Ontology (DO) assessments. For GO analysis, both unadjusted and adjusted p-value thresholds were established at less than 0.05. For the KEGG and DO analyses, the significance threshold was set at a p-value of less than 0.1.

Machine learning methodology

Feature selection: Robust feature selection methods were employed before model training to harness the full potential of our high-dimensional spatial transcriptomics data. We utilized Sequential Feature Selection (SFS), Lasso Regression, and Recursive Feature Elimination with Cross-Validation (RFECV) within Python to identify significant predictors. We incorporated the LogisticRegressionCV and RFECV modules from Python’s scikit-learn library for implementing Lasso Regression and RFECV. SFS was implemented using a custom-written algorithm in Python. Additionally, to explore the overlapping genes among SFS, Lasso Regression, and RFECV, we employed the venn3 function from Python’s matplotlib library to generate Venn diagrams. To further understand the impact of various features on model predictions, we applied the SHAP (SHapley Additive exPlanations) library in Python to conduct SHAP value analysis, revealing the relative importance of different features within the model [25].

• LASSO: This method imposes a penalty proportional to the absolute value of the regression coefficients, effectively reducing less significant feature coefficients to zero. It is particularly beneficial in scenarios with high multicollinearity or when simplifying the model’s interpretation by minimizing the number of predictor variables.
• RFECV: This technique iteratively removes features while building models and assessing accuracy via cross-validation. It aims to isolate a subset of features that most significantly predict the target variable, thereby enhancing computational efficiency and predictive accuracy. Specifically, RFECV utilizes the LGBM (Light Gradient Boosting Machine) model to ensure robust feature selection and model training.
• SFS: Unlike RFECV, which removes features, SFS adds them sequentially based on a specific criterion until no further improvements are seen in model performance. This approach is advantageous for models sensitive to feature interactions, as it carefully considers the effects of each new feature added to those previously selected. SFS is also implemented using the LGBM model, allowing for precise and effective feature integration.

Model training: After feature selection, we trained various machine learning models to classify disease status based on gene expression data. The models included:

• Convolutional Neural Networks (CNN): Implemented using TensorFlow and Keras libraries, the CNN model was trained on visualized gene expression images. Hyperparameters such as learning rate, batch size, and number of epochs were optimized using grid search. The input images were standardized, and data augmentation techniques were applied to enhance model robustness.

The CNN model was trained through a visualized image of the genes. The 5404 genes provided by the dataset were filled into a 74x74 image, each gene representing a pixel. The grey scale intensity of each pixel/gene, on a scale from 0-255, was mapped from the normalized expression count values provided by StomicDB [27]. Figure 2 is an example visualization of the image input into the CNN model. Each square on the figure represents the normalized count of a gene.

• Light Gradient Boosting Machine (LGBM): The LGBM model was implemented using the LightGBM library. We set the number of leaves to 31, the learning rate to 0.05, and used 100 boosting iterations. Early stopping was applied based on validation loss to prevent overfitting.
• Support Vector Machine (SVM): The SVM model was implemented using scikit-learn’s SVC class with a radial basis function (RBF) kernel. Hyperparameters such as the regularization parameter C and kernel coefficient γ were optimized using cross-validation.

Addressing Class Imbalance: Given the limited and potentially imbalanced dataset, we employed strategies to mitigate class imbalance effects. Techniques such as synthetic minority over-sampling (SMOTE) were considered, but due to the small sample size, we opted for stratified cross-validation to maintain the original class distribution in training and validation splits. Performance metrics were selected to account for imbalances, focusing on area under the curve (AUC) and F1 score, which provide a balanced assessment of model performance.

Model Evaluation: The final step involved evaluating each model’s performance, to ensure the robustness and generalizability of the predictive models. Metrics such as Area Under the Curve (AUC), accuracy, precision, sensitivity, and F1 score were computed to evaluate each model comprehensively.

Result of SPA

DEGs between diseased and normal samples in SPA were screened with a criteria of p <0.05 and |logFC| >0.6. As depicted in the volcano plot in (Figure 3A), there were 49 DEGs identified between the two groups, all of which were upregulated. The corresponding heatmap is shown in (Figure 3B). Intersecting the DEGs identified by the LASSO model, RFECV, and SFS resulted in the Venn diagram shown in (Figure 3C), including XBP1, UBC, SSR4 and HLA-A. The descriptions of these four DEGs are listed in (Table 1).

SHAP summary plots were utilized to elucidate the roles of DEGs identified through various feature selection methods in assessing SPA risk. The DEGs identified by LASSO, RFECV, and SFS were systematically ranked based on their mean absolute SHAP values, highlighting the genes most significantly associated with increased SPA risk (Figures 4A, 4C, and 4E). Specifically, the top ten DEGs identified by LASSO included VIM, PSAP, HLA-A, DDX5, PIM2, C4B, DDIT4, ATP5F1E, SERP1, and HSP90B1, as shown in (Figure 4B). The RFECV method identified PSAP, VIM, HLA-A, PIM2, CALR, XBP1, DDX5, DDIT4, HSP90B1, and B2M as the principal genes (Figure 4D), while SFS highlighted XBP1, UBC, SSR4, HLA-A, and BTG2 as the most influential (Figure 4F). Genes with higher SHAP values were found to exert a greater impact on SPA, indicating their potential as critical biomarkers or therapeutic targets.

Result of RA

DEGs between diseased and normal samples in RA were identified using a cutoff of p <0.05 and an absolute logFC greater than 0.6. As illustrated in the volcano plot (Figure 5A), a total of 30 DEGs were detected, all of which were upregulated. The corresponding heatmap is displayed in (Figure 5B). Additionally, the Venn diagram in (Figure 5C) shows the overlap of DEGs identified by the LASSO model, RFECV, and SFS, including VIM, TRAM1, FTL, ATP5F1E, TPT1, SARAF, FAU, B2M, and SSR4. The descriptions of these nine DEGs are listed in (Table 2).

The SHAP summary plots were employed to elucidate the impact of DEGs identified through various feature selection methods, contributing to our understanding of their roles in RA. The DEGs identified by LASSO, RFECV, and SFS were ranked based on their mean absolute SHAP values, highlighting the genes most associated with increased RA risk (Figures 6A, 6C, and 6E). The top-10 mean SHAP values for LASSO, ranked from highest to lowest, included CST3, HLA-B, MMP3, CD52, PTGDS, EVL, TMSB4X, PFN1, HLA-E, and ACTB (Figure 6B). For RFECV, the top-10 leading genes were CST3, CD52, PFN1, HLA-E, EVL, SSR4, UBA52, MMP3, HLA-B, and TMSB4X (Figure 6D). SFS highlighted CD52, ACTB, FAU, EVL, MMP3, TMSB4X, SSR4, B2M, PTGDS, and HLA-C as the top influencers (Figure 6F). Genes depicted with higher SHAP values suggest a greater influence on RA.

The functional annotation of 4 DEGs for SPA was conducted using GO, KEGG, and DO enrichment analyses. The GO enrichment analysis demonstrated that these DEGs are primarily enriched in immune-related functions. Specifically, within the subset of Biological Processes (BP), the genes are mainly involved in T cell differentiation (Figure 7A). In the Cellular Component (CC) subset, enrichment was observed in the vesicle lumen (Figure 7B), and in the Molecular Function (MF) subset, the genes showed significant activity in immune receptor functions (Figure 7C). These findings highlight the critical role these genes play within the immune system, potentially influencing the development, communication, and response of immune cells. Additionally, the KEGG pathway analysis revealed significant enrichment in processes such as protein processing within the endoplasmic reticulum (ER), interactions with Kaposi sarcoma-associated herpesvirus, and the pathophysiology of Parkinson’s disease, among other neurodegenerative disorders. The analysis also highlighted pathways related to immune-related disorders including allograft rejection, Type I diabetes mellitus, graft-versus-host disease, and autoimmune thyroid disease. Furthermore, the DO enrichment analysis demonstrated that the differentially expressed genes were associated with a range of diseases such as multiple myeloma, vulva cancer, myeloid neoplasms, borna disease, bone marrow cancer, fibrous histiocytoma, and histiocytoma, underscoring their potential implications in diverse pathological conditions.

GO, KEGG, and DO enrichment analysis

The 9 DEGs for RA underwent functional annotation through GO, KEGG, and DO enrichment analyses. The results revealed that these DEGs are predominantly enriched in immune-related functions. Specifically, they are involved in the negative regulation of binding within the BP category (Figure 8A), localized to the rough ER in the CC category (Figure 8B), and associated with profilin binding in the MF category (Figure 8C). In addition, the significantly enriched KEGG pathways included Rap1 signaling pathway, regulation of actin cytoskeleton, coronavirus disease - COVID-19, Vibrio cholerae infection, arachidonic acid metabolism, viral myocarditis, gastric acid secretion, and bacterial invasion of epithelial cells (Figure 8D). Finally, the DEGs were closely related to various diseases as per the DO analysis, including carotid artery disease, Moyamoya disease, cerebrovascular disease, amyloidosis, intracranial arterial disease, cerebral arterial disease, osteoarthritis, and bone sarcoma (Figure 8E).

Model evaluation

Table 3 list performance Metrics of CNN, LGBM, and SVM Models. For the diseases SPA and RA, the LGBM model demonstrated superior performance across most metrics, particularly excelling in AUC (0.897) and Accuracy (0.933) for SPA. The SVM model led in Sensitivity (0.754) and F1 Score (0.694) for SPA, while the CNN model showed strong performance in AUC (0.906) and Accuracy (0.899) for RA. Each model displayed varied results in Precision and Sensitivity, reflecting their applicability and balance in different clinical scenarios.

The identification of differentially expressed genes (DEGs) inThe identification of differentially expressed genes (DEGs) in this study provides valuable insights into potential biomarkers for RA and SpA. To contextualize these findings, we compare the DEGs identified here with known biomarkers from related autoimmune conditions, focusing on common and unique biomarkers that could enhance disease understanding and therapeutic strategies.

Both RA and SpA share several biomarkers associated with immune activation and inflammation. For instance, HLA antigens like HLA-A and HLA-B have In fact, been proven in both the conditions. HLA-B27 antigen is a confirmed risk factor gene in SpA with significant focus on AS and it has been shown that this antigen has significant role to play in antigen presentation process of immune reaction. Thus, the results of this study offer a rationale for the use of HLA-A as candidate biomarker of immune dysregulation and disease activity, and RA and SpA patient cohorts [32,33].

In addition to HLA markers, other shared biomarkers include MMP3, CD52, and B2M. These molecules play a role in immune moderation, tissue repair, and antigen display, and are related to both RA and SpA [34]. For instance, MMP3 is understood to mark joint degradation in RA and is sensitive to the progressive severity of the pathological condition [35]. Likewise, B2M has been found to have significant positive associations with disease activity in both RA and SpA due to immune response [36].

While many biomarkers are shared between RA and SpA, there are also disease-specific markers that may offer new avenues for diagnosis and treatment. For instance, XBP1 was found in our study as a DEG that is specific to SpA. This transcription factor is involved in unfolded protein response (UPR) and cellular stress, which is highly implicated in SpA, where immune cell infiltration and tissue remodeling occur [37]. Likewise, UBC (Ubiquitin C), which controls protein degradation, was also overexpressed in the SpA samples and has been associated with immune cell activation to inflammation in autoimmune disorders [38]. According to the results of RA genes, ACTB (Actin Beta) and PTGDS (Prostaglandin D2 Synthase) could be selected as significant biomarkers in the current work. ACTB contributes to cytoskeletal rearrangement which is important for immune cell migration in RA and therefore crucial in the disease progression [39]. PTGDS synthesized and released prostaglandin D2, which has both inflammatory and anti-inflammatory effects in RA [40].

To confirm the DEGs found in this study, the following experimental assays should be performed; qPCR, IHC, Western blotting, and FACS. These methods will not only validate gene expression but also offer further understanding of the functional characteristics of these biomarkers in the development of RA and SpA.

Real-time PCR is one of the most frequently utilized methods for the validation of gene expression. It shows alteration in the gene expression that is either up or down-regulated and quantifying mRNA ensures that there is a real and accurate increase of expression levels. It is particularly useful for ensuring the upregulation of DEGs in both RA and SpA. For instance, the HLA-A and CD52 with their roles in immune system activation and antigen presentation could be further examined using qPCR to determine the expression of the corresponding mRNA in affected tissues [41].

Another valuable technique of affirming gene expression at protein level is Immunohistochemistry (IHC). IHC also allows one to study the distribution of proteins within tissues; this embeds gene expression data into correct three-dimensional context. This method would be very useful for the validation of the presence and the location of protein such as HLA-A and CD52 in synovial tissue of samples of patients with RA and SpA. Hence, by using respective antibodies for these proteins, IHC would be applied to measure the inflamed tissue for amount of HLA-A, or CD52 proteins if the HLA-A or CD52 gene is transcribed as equal corresponding amounts of proteins. It also contains valuable data regarding the distribution of such biomarkers in tissue, which in turn helps in understanding the processes of immune cell infiltration as well as tissue remodeling [42].

In addition to qPCR and IHC, Western blotting is also very effective in protein validation. This technique could be used to confirm proteins like MMP3 and PTGDS that play the role of tissue remodeling and inflammation respectively. The direct quantity of protein can be measured through procedures such as the western blotting method which will also provide a strong evidence of protein expression at the post-translational level. Thus, confirmation of these proteins’ upregulation will aid in identifying their roles in joint degradation and inflammation in both RA and SpA [43].

Finally, the expression of CD 52 or any other immune cell stain on the specific cell type can be confirmed by flow cytometry. This method comes in handy as it permits viewing of protein expression on individual cells while especially applicable for analysis of immune cell activation and phenotype. For instance, flow cytometry may be used to determine the relative content of CD52 in T cells or B cells that infiltrate the synovial tissue in RA and SpA and demonstrate roles of this protein in immunopathogenesis and disease progression [44].

Regarding the DEGs uncovered in the current study, including HLA-A, CD52, MMP3, and UBC, the experimental validation utilizing these methods will offer additional evidence as to the biomarker candidacy of RA and SpA. Therefore, HLA-A and CD52 are singled out for further validation because of their roles in the immune system and as potential therapeutic targets. qPCR and IHC will be particularly applicable to these markers, as the technique will further indicate whether the genes and proteins of interest are expressed in diseased tissues. Second, due to the strong association of MMP3 and PTGDS with tissue remodeling and inflammation, these two proteins should be subjected to western blot analysis to verify their protein level and explore their biological function in disease development.

However, the small sample size in this study is the biggest restriction that has to be taken into account in terms of both statistical and generalizability effects. Statistical power in turn is the capacity of a given study to identify a genuine effect if indeed there is one. A small number of participants tends to lower the sensitivity of the study—the probability of a Type II error rises, meaning that the researcher may fail to see a difference or an effect when they do exist.

In studies involving small cohorts, the reliability of findings becomes more uncertain. Small samples are generally more variable than large samples since their representation of the whole population is much smaller. This means that it is very difficult to generalize results based on a small sample, which is often the case when conducting such studies. Besides, narrower confidence intervals are obtained from studies with larger sample sizes; therefore, the studies included in this analysis had a small sample size leading to broader confidence intervals, and therefore the estimates are less precise. However, due to the small sample size, the CIs were wider and repeating the study with a bigger sample would give inconsistent results [45]. A primary challenge in studies with small sample sizes is low statistical power, which reduces the ability to detect significant differences, even when they exist. Small sample sizes increase the risk of not detecting an effect that exists due to inadequate subjects. Also, in any given study, differences in gene expression between two different people may not be truly significant, and variability obscures real differences. Small cohorts further have more artefacts and cannot be extrapolated to the general RA or SpA population with the consequent risk of making conclusions that only apply to the study sample. In addition, many subjects may restrict the choice of statistical techniques, and subsequently, the reliability of the results can be less strong as compared to those samples constituting large and diverse populations [46].

In light of these challenges, future investigations need to expand the number of participants to boost statistical conclusiveness. A sample that is involved with a large number will offer a good estimate and also enable the detection of small effects that could biologically be important. Also, having more participants requires the investigators to be certain that all patients are selected, hence the findings will be generalizable to the population.

To enhance the reliability and generalizability of findings, future studies should aim to include larger and more diverse cohorts of patients. Increasing the sample size will increase the statistical power and thus conduct more robust tests on the data as well as increase the chances of forming a sound conclusion as to the biological variability in disease groups. Larger cohorts also minimize some of the factors that arise when sample sizes are small; such factors include sampling bias and variation in the patterns of gene expression in other individuals.

Increasing the size of the cohort in autoimmune diseases, especially the RA and SpA, can be made easier by partnering with other centers. In this way, undertakings like these collaborations where many research centers accumulate data from their patient samples can yield a more representative database of disease patients included in their studies. Regional factors can also be addressed in multi-center trials, and this enhances the external validity of the results in different geographical areas and other characteristics age, sex, ethnicity, and severity of disease amongst the patients. Such an approach may address the difficulties of a single-center study in terms of the accessibility to patients.

For instance, the International Genetics of Rheumatoid Arthritis Consortium (GRA) and the Spondyloarthritis Research Consortium of Canada (SPARCC) are large, centralized genetic studies that have previously merged data from multiple centers to compare the genetic factors and clinical attitudes of RA and SpA. Such large consortia have allowed for more extensive analyses and greater statistical confidence for the identification of disease etiology.

In this study, we leveraged advanced spatial transcriptomic technologies and machine learning algorithms to investigate the molecular underpinnings of RA and SpA. By analyzing high-dimensional spatial transcriptomics data from synovial biopsy samples, we identified DEGs that may serve as potential biomarkers or therapeutic targets. Our integrative approach allowed us to discern specific genes and pathways implicated in the pathogenesis of RA and SpA, providing novel insights that could enhance diagnostic precision and inform targeted therapies.

Differentially expressed genes and their implications

Our analysis revealed 49 DEGs in SpA and 30 DEGs in RA with significant differential expression compared to control samples. Through rigorous feature selection methods, including LASSO, RFECV, and SFS, we refined these to 4 key DEGs for SpA (HLA-A, UBC, XBP1, SSR4) and 9 for RA (CD52, ACTB, PTGDS, SSR4, EVL, MMP3, FAU, B2M, TMSB4X). The identification of these genes underscores their potential roles in the pathogenesis of these diseases.

In SpA, the DEGs identified are predominantly involved in immune regulation and cellular stress responses. For instance, HLA-A is a major histocompatibility complex class I gene crucial for antigen presentation to cytotoxic T cells, implicating it in immune surveillance and response [47,48]. The strong association of SpA, particularly ankylosing spondylitis, with HLA-B27 has been well-documented, and our findings suggest that other HLA class I molecules like HLA-A may also contribute to disease pathogenesis.

The ubiquitin C gene (UBC) plays a critical role in protein degradation through the ubiquitin-proteasome system, which is essential for regulating immune responses and maintaining cellular homeostasis [49]. Dysregulation of ubiquitination processes can lead to aberrant immune activation, a hallmark of SpA. XBP1 is a key transcription factor in the UPR, which manages ER stress [50,51]. Chronic activation of the UPR has been implicated in inflammatory diseases, suggesting that XBP1 may contribute to the persistent inflammation observed in SpA. SSR4, encoding the translocon-associated protein δ, is involved in protein translocation across the ER and has been linked to immune cell infiltration in various tissues [52]. The upregulation of SSR4 may reflect enhanced protein synthesis demands in activated immune cells within inflamed synovial tissue.

In RA, the identified DEGs are integral to immune modulation, cytoskeletal organization, and tissue remodeling. CD52 encodes a glycoprotein expressed on the surface of lymphocytes and other immune cells, playing a vital role in immune regulation [53]. Therapeutic targeting of CD52 with monoclonal antibodies like alemtuzumab has shown efficacy in depleting lymphocytes and modulating autoimmune responses [54,55], highlighting its potential as a therapeutic target in RA.

ACTB, encoding β-actin, is essential for cytoskeletal integrity and cell motility. Altered expression of β-actin can affect immune cell migration into the synovium, exacerbating inflammation and joint destruction [56,57]. PTGDS, which encodes prostaglandin D2 synthase, is involved in the production of prostaglandin D2 (PGD2), a lipid mediator with complex roles in inflammation [58]. PGD2 can exhibit both pro-inflammatory and anti-inflammatory effects depending on receptor interactions [59,60], suggesting that PTGDS expression levels may influence inflammatory pathways in RA.

The consistent identification of SSR4 in both RA and SpA underscores its potential importance in inflammatory joint diseases. Its role in protein translocation and the UPR may contribute to synovial hyperplasia and immune cell activation [61]. EVL, encoding the Enah/VASP-like protein, is involved in actin cytoskeleton remodeling, affecting immune cell migration and potentially the invasive behavior of synovial fibroblasts in RA [62-64].

MMP3 is a matrix metalloproteinase that degrades extracellular matrix components, facilitating tissue remodeling and contributing to cartilage destruction in RA [65,66]. Elevated levels of MMP3 correlate with disease severity, making it a valuable biomarker and potential therapeutic target. FAU, involved in protein degradation and ribosome synthesis, may influence immune responses and cell proliferation through the ubiquitin-proteasome system [67,68]. B2M, a component of MHC class I molecules, is essential for antigen presentation, and its elevated levels have been associated with increased disease activity in RA [69-71]. TMSB4X, encoding thymosin beta-4, plays a role in actin binding and may affect immune cell function and tissue repair mechanisms [72,73].

Interpretation of upregulated DEGs

An unusual finding in our analysis was that all identified DEGs in both RA and SpA were upregulated. Typically, differential expression studies reveal a mix of upregulated and downregulated genes. This pattern may suggest a global activation of certain pathways or cellular processes in diseased tissues compared to controls. One possible explanation is that the diseased synovial tissues exhibit heightened inflammatory activity, leading to increased expression of genes involved in immune responses and cellular stress. However, this could also result from limitations in our dataset or analytical approach. The small sample size and lack of technical replicates may have affected the sensitivity of detecting downregulated genes. Additionally, our criteria for DEG selection might have favored the identification of upregulated genes. This limitation highlights the need for cautious interpretation of the results and underscores the importance of validating these findings in larger, independent cohorts.

Biological explanations for upregulated DEGs

The finding that all DEGs in our study were upregulated in both RA and SpA is an atypical result. Generally, DEG analysis brings out both upregulated and downregulated genes, which is a sign of study on gene expression in diseases. Nonetheless, increased expression of the DEGs solely in these autoimmune diseases may be due to several reasons that are biological and interrelated, namely immune activation, inflammation, and stress.

Immune system activation in RA and SpA: RA and SpA are both autoimmune diseases that involve inflammation and immune overactivity, which might explain the global increases in gene activation that were detected. In these conditions, the immune system of the body becomes self-destructive and attacks the healthy tissues but most affected are the joint tissues as they always exhibit continuous inflammation. This is due to the continuous stimulation of several immune cells like T cells and macrophages leading to the enhancement of experimental expression of genes related to immune signaling. For example, HLA molecules play a role in antigen presentation in immune responses; depending on HLA types, these molecules may be upregulated in autoimmune diseases [74]. Such overactivation could well account for the worldwide increase in the expression of genes associated with immune response, cell migration, and tissue remodeling [75].

Cellular stress responses in chronic inflammation: The chronic inflammation observed in RA and SpA induces cellular stress, particularly within immune cells and synovial tissue. Another possible reason for the increased DEGs is the activation of the UPR. This response occurs when cells are under stress brought about by the presence of misfolded proteins, which are characteristic of chronic inflammation [76]. Both in RA and SpA immune cells invade the affected tissues requiring increased biosynthesis of proteins and utilization of quality control mechanisms. Some of the specific targets of ER stress include genes like XBP1 and SSR4, which can be upregulated due to the attempt made by the cell to respond to this specific stress [77]. In this regard, the global increase in gene expression could be indicative of a cellular attempt at restoration of normalcy when under stress, a process which may be generalized in inflamed tissues where protein synthesis is increased.

Tissue remodeling and immune cell infiltration: Another potential cause for the upregulation of all DEGs is the extensive tissue remodeling occurring in RA and SpA. These diseases are characterized by synovial inflammation where the synovial intima layer known as the synovial membrane becomes thicker and is made up of new fibroblast cells as well as immune cells. T cells, B cells, and macrophages are invoked in the synovium and cartilage and regulate genes that are involved in inflammation and remodeling of the matrix. For example, MMP3 and ACTB which are involved in tissue remodeling may be upregulated for cellular migration and matrix degradation in the inflamed synovium [78].

Chronic inflammatory state and immune dysregulation: Immune dysregulation in RA and SpA, additionally, entails a chronic elevation of inflammation in which cytokines such as TNF-α, IL-6, and IL-17 are implicated. These pro-inflammatory cytokines reinforce the signals of inflammation to continue the inflammatory process and promote the increased expression of genes concerning immune response, cell attachment, and cytokine signaling [79]. This enhancement of immune responses could lead to the exclusive increase of DEGs in our study, indicating an immune system that has been heightened to fight threats, threats which may not necessarily be posed by pathogens.

Possible mechanisms behind global upregulation: Several biological processes may be likely to attribute to the global upregulation of genes of RA and SpA; as both the disorders are characterized by persistent inflammation, one can also note that immune-related gene and cytokine signaling pathways are constantly active in these conditions and promote immune responses. However, ER stress and the UPR upregulated genes related to protein quality control, including immune-activated cells under stress conditions. Finally, the path of cytokine dysregulation factor causes an imbalance of pro-inflammatory cytokines which in turn activate other immune cells and result in constant gene expression programs. Altogether, these mechanisms are attributed to the upregulation of several genes that persist due to immune activation and cellular stresses in the inflamed joint.

Functional pathway analysis

Gene Ontology (GO), KEGG pathway, and Disease Ontology (DO) enrichment analyses provided further insights into the biological functions and pathways associated with the identified DEGs. In SpA, the DEGs were enriched in processes related to T cell differentiation, vesicle lumen localization, and immune receptor activity (Figure 7). KEGG pathways included protein processing in the ER and immune-related disorders, emphasizing the role of protein homeostasis and immune responses in SpA pathogenesis. The DO analysis linked these genes to diseases such as multiple myeloma and bone marrow cancer, suggesting common underlying mechanisms involving immune dysregulation.

For RA, GO analysis highlighted enrichment in negative regulation of binding and localization to the rough ER (Figure 8), indicating disruptions in protein processing and cellular trafficking. KEGG pathways implicated included the Rap1 signaling pathway, which influences cell adhesion and migration [80], and regulation of the actin cytoskeleton, essential for immune cell movement and synovial fibroblast invasion. The association with viral myocarditis and other diseases in the DO analysis may reflect shared inflammatory mechanisms.

These findings align with the current understanding that both RA and SpA involve complex interactions between immune cells, cytokines, and cellular stress responses. The identification of common pathways and genes provides a foundation for exploring shared therapeutic targets.

Integration of machine learning models

The application of machine learning models enhanced our ability to predict disease status based on gene expression data. The use of SHAP values allowed us to interpret the contribution of individual genes to the model predictions, adding transparency to the machine learning process [81]. The LGBM model demonstrated superior performance in AUC and accuracy for both RA and SpA (Table 3), suggesting its suitability for high-dimensional genomic data.

The CNN model, applied to visualized gene expression images, also showed strong performance, particularly in RA. This indicates that spatial patterns of gene expression captured through imaging can be effectively utilized for disease classification. The SVM model excelled in sensitivity and F1 score for SpA, highlighting its potential in identifying true positive cases.

However, it is important to acknowledge that the moderate precision and sensitivity values indicate room for improvement. The limited sample size may have affected the generalizability of the models, and future studies with larger cohorts are necessary to validate and refine these predictive models.

Limitations and future directions

While our study provides valuable insights, several limitations must be considered. The sample size was relatively small, which may limit the statistical power and generalizability of the findings. All data were derived from publicly available datasets and validated in silico; therefore, experimental validation using independent cohorts is essential to confirm the relevance of the identified DEGs.

However, despite functional annotations based on these results, continued and more intricate molecular characterization is needed for understanding the roles of these DEGs in RA and SpA pathogenesis. A study of these genes’ dynamics as well as the functions of the genes within distinct pathways may shed light on new targets for treatment.

Therefore, built on the results of this research, future research should concentrate on the expanding of the sample size to increase the statistical strength to the findings and generalizability. Incorporation of clinical researchers could also bring in more synovial biopsy samples from the other groups of patients as well. For the purposes of this review, experimental confirmation of the identified DEGs with techniques such as quantitative PCR or immunohistochemistry is necessary. The functional studies focused on examining the functions of these genes in disease development may give further understanding of their therapeutic function. Furthermore, if other related omics data experiences proteomics and metabolomics were included, there would be more accurate understanding of RA and SpA. With the clearer class sets and enhanced ability of the larger datasets to model more specific slices of reality, the models are expected to produce enhanced accuracy and reliability. Finally, it will be perspective to translate these findings into clinical practice and to evaluate the scenarios that defined biomarkers are really diagnostic and prognostic factors.

Recommendations for future studies

Research on RA and SpA in the future must contemplate sample size to have adequate power. This can be done by having multiple centers where recruitment is done to incorporate various areas of the world and patients with different characteristics and distribution of disease forms. The next priority is patient characteristics, for example, age, gender, race, and disease severity, to compare manifestations of RA and SpA in their different subgroups. Furthermore, since the new biomarkers will involve CD52 and HLA-A, longitudinal research that will follow the history of diseases will be useful in evaluating how these biomarkers change with time and the impact on the clinical results of the treatments administered.

Furthermore, the studies at the molecular level should be extended in the context of RA and SpA by identifying the multi-omics data, including genomics, proteomics, and metabolomics, associated with clinical data in the future. Multiple countries and healthcare settings will also be useful in confirming the results and minimizing sampling errors, making results generalizable. These approaches will also enhance the accuracy and validity of tested outcomes while increasing the sample statistic, therefore, shaping a strategy to implement suitable therapeutic intercessions for a variety of RA and SpA patients to increase their effective treatment and management.

Our integrative approach combining spatial transcriptomics and machine learning has identified key DEGs and pathways involved in RA and SpA. The findings enhance our understanding of the molecular mechanisms underlying these diseases and highlight potential biomarkers and therapeutic targets. While preliminary, these results lay the groundwork for future studies aimed at improving diagnostic accuracy and developing targeted therapies, ultimately contributing to better patient outcomes.