Identification of Disulfidptosis-Related Genes and Molecular Subgroups in Rheumatoid Arthritis for Diagnostic Model and Patient Stratification

Introduction

Rheumatoid arthritis (RA) is a systemic autoimmune disease characterized by persistent synovitis and the destruction of joint cartilage and bone.¹ Additionally, it affects the function of various organs, such as the heart, lungs, and kidneys, leading to a decline in quality of life.^2,3 The prognosis of RA is influenced by early diagnosis and proper treatment.^4,5 The biological and targeted synthesis disease modifying anti-rheumatic drugs (DMARDs) have been demonstrated to yield superior outcomes for patients suffering RA compared with conventional DMARDs. However, approximately 30%–40% patients have no response to the current treatment. Additionally, no treatment can be demonstrated to fully cure the disease.⁶ In accordance with the well-known heterogeneity of RA, variable prognosis and response to RA treatment are linked to the variable underlying pathogenesis in individual patients.^7–9 The stratification of RA patients is paramount for the evaluation of prognosis and the selection of appropriate treatment options.^10,11 However, this is challenging due to the lack of clarity regarding the pathogenesis of RA.

It is generally accepted that RA is associated with genetic susceptibility^12–14 and environmental factors.^15,16 Recent studies indicated that a variety of immune cells, including B cells,^17,18 T cells,^19,20 macrophages,²¹ dendritic cells,^22,23 mast cells,²⁴ and neutrophils,²⁵ were involved in the development and progression of RA. Furthermore, cell death contributes to the pathogenesis of RA through various pathways. The diverse cell death modes, including ferroptosis,^26,27 cuproptosis,²⁸ and pyroptosis,^29,30 were involved in RA. For instance, the inhibition of ferroptosis may result in a reduction of RA synovium proliferation by means of the alleviation of oxidative stress response, and may also affect the polarization of macrophages to reduce the inflammatory response.^31,32 Liu et al identified disulfidptosis, as a novel form of cell death, was associated with abnormal accumulation of intracellular disulfide bonds, particularly under conditions of glucose deprivation. Cells with high expression of SLC7A11 increase intracellular disulfide bond formation under glucose-starved conditions by promoting cystine uptake. This accumulation depletes intracellular NADPH and then leads to disulfidptosis.³³ Unlike traditional apoptosis or necrosis, disulfidptosis promotes cell death by altering the disulfide bond status of cellular scaffold proteins, a process closely linked to cellular redox balance.³⁴ It is common for RA synovial tissue characterized by anoxia present of oxidative stress, namely redox imbalance.^35,36 However, to our knowledge, there are no studies on disulfidptosis in RA.

To identify disulfidptosis-related genes associated with RA, we first used transcriptomic data from public databases to confirm disulfidptosis-related genes by weighted gene co-expression network analysis (WGCNA) and differential analysis, and then applied three machine learning algorithms, LASSO regression, SVM-RFE, and Random Forest, to find hub genes. Second, we constructed a prediction model for RA diagnosis using the hub genes. To evaluate the diagnosis value of the model, we performed internal and external validation. Third, we determined two RA subtypes based on the hub genes and compared the disease activity scores, clinical markers, immune checkpoints and immune cells between 2 subgroups of RA patients, trying to indicate the hub genes for RA patient stratification. Finally, the expression of hub genes was experimentally validated between RA patients and healthy controls.

The objective of this study was to construct a promising diagnostic model and validate its diagnostic value, and then provide a patient stratification, effectively distinguishing two RA subtypes by distinct clinical and molecular constituents. This has the potential to identify new therapeutic targets and enhance personalized treatment strategies for RA patients.

Methods

Dataset Sources and Preprocessing

The dataset used in the current study was obtained from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) using “rheumatoid arthritis” as the search term. Data sets were selected according to the following criteria: gene expression profiles were from human samples; samples included both healthy controls and patients. Three datasets were selected from GSE110169³⁷ and GSE93272,³⁸ which contained whole-blood gene expression data using microarray analysis, and GSE89408^39,40 containing synovial biopsy data by high-throughput sequencing. The information of three datasets is shown in Table 1.

Table 1 GEO Datasets Information

The study flowchart was shown (Figure 1). The raw datasets were normalized using the Limma package. In instances where multiple probes were mapped to an identical gene, the median of these values was assigned as the gene’s expression level. This process resulted in the formulation of a refined expression matrix by converting the probe identifiers to their respective gene symbols. The GSE89408 dataset was log2 transformed to normalize the variance within gene expression levels and prepared the data suitable for further analytical processes. Meanwhile, the GSE93272 dataset, which presented with gaps in clinical data, underwent imputation via the k-nearest neighbor (KNN) algorithm, leveraging the 10 nearest neighbors to surmise the missing entries.⁴¹

Figure 1 Flowchart of the study. Abbreviations: DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RA, rheumatoid arthritis; ROC, receiver operating characteristic; ssGSEA, Single-sample Gene Set Enrichment Analysis; SVM-RFE, Support Vector Machine Recursive Feature Elimination; WGCNA, Weighted Gene Co-expression Network Analysis.

Weighted Gene Co-Expression Network Analysis and Identification of Significant Module Genes

The Weighted Gene Co-expression Network Analysis (WGCNA) was performed using the WGCNA R package (version 1.70.3). First, sample clustering was conducted to identify potential outliers (Supplementary Figure 1A). Second, a scale-free co-expression network was constructed using the pickSoftThreshold function in R to compute the soft thresholding power. The soft power was determined to be 14 (scale independence of 0.9) (Supplementary Figure 1B). The adjacency matrix was then transformed into a topological overlap matrix (TOM), and the gene ratio along with its corresponding dissimilarity (1-TOM) were computed. Then, hierarchical clustering and dynamic tree-cutting methods were applied to cluster the genes. The minimum number of gene modules was set to 30, and an appropriate cutoff value was chosen to merge similar modules (Supplementary Figure 1C and D). Finally, gene significance (GS) and module membership (MM) were calculated to explore the relationship between modules and clinical traits. Meanwhile, we obtained a heatmap illustrating the correlation between 20 gene modules and clinical traits.

Identification of Differentially Expressed Genes (DEGs) Associated With Disulfidptosis

According to the previous research, ten genes were associated with disulfidptosis, including GYS1, LRPPRC, NCKAP1, NDUFS1, NDUFA11, NUBPL, OXSM, SLC3A2, RPN1, and SLC7A11.³³ To identify DEGs, we estimated disulfidptosis scores within the 232 RA samples of the GSE93272 dataset using the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm in the GSVA package (version 1.46.0). The RA samples were divided into two groups with high and low disulfidptosis scores, based on the median score. Consistent with the previous study,⁴² differential expression analysis between the two groups was performed using the limma package (version 3.54.0) with P value < 0.05 and |log2FC| > 0.25. The genes meeting these criteria were considered as DEGs. These genes were displayed using volcano plot and heatmap, providing a clear visual differentiation of gene expression to facilitate the interpretation of the results.

Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis

The pathway enrichment analysis provided by KEGG served as a comprehensive knowledge base for systematic gene function analysis.⁴³ GO analysis, covering Biological Processes (BP), Molecular Functions (MF), and Cellular Components (CC), yielded invaluable insights into gene functionality.⁴⁴ Enrichment analysis was performed by the R package “clusterProfiler (version 3.14.3)”,⁴⁵ and the top 10 GO terms in each category were visualized using the R package “ggplot2”. These analyses were based on screening criteria: P value < 0.05.

Machine Learning Algorithms and Hub Genes Identification

To identify hub genes, three machine learning algorithms were employed for variable selection in this research. LASSO regression constructed generalized linear models and was used to mitigate model overfitting, which incorporated penalty terms to shrink unimportant feature coefficients to zero, effectively facilitating variable selection.⁴⁶ Support Vector Machine Recursive Feature Elimination (SVM-RFE) was applicable to small datasets and was an algorithm that retained variables relevant to outcomes while effectively eliminating redundant factors.⁴⁷ Additionally, Random Forest was capable of handling high-dimensional data, allowing predictive models to be built and ranking the importance of feature variables.⁴⁸ The intersection of genes obtained from the three machine learning algorithms was hub genes.

Single-Sample Gene Set Enrichment Analysis (ssGSEA)

We performed ssGSEA to identify the relationship between hub genes and hallmark gene sets.⁴⁹ First, we downloaded the h.all.v7.4.symbols.gmt subset from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/downloads.jsp) to identify relevant pathways and molecular mechanisms. Then, we calculated the enrichment scores of each sample in the gene sets using gene expression profiles. Finally, an enrichment score matrix was obtained, and the results were visualized as a correlation heatmap.

The Construction and Validation Accuracy of the Nomogram

To construct a diagnostic model, we employed a regression equation in the GSE110169 dataset. Subsequently, we visualized the nomogram using the R package “rms (version 6.2.0)”. The nomogram was based on the regression coefficients of all features. We established scoring criteria for each feature’s different expression levels. By calculating the total score for each sample, we were able to predict the probability of each sample developing RA. Furthermore, we evaluated the accuracy of the model predictions and validated the predictive value of hub genes using ROC curve, calibration curve, and decision curve.

External Verification of Nomogram and Hub Genes

We performed external validation of the nomogram and hub genes using the GSE93272 and GSE89408 datasets. The GSE93272 dataset included 275 whole blood samples, of which 232 were RA samples and 43 were controls. The GSE89408 dataset comprised 180 synovial biopsy samples, including 152 RA samples and 28 controls. ROC curves were plotted for each dataset to assess the predictive accuracy.

Immune Infiltration Analysis

The “Cibersort” algorithm was used in the GSE93272 dataset to convert the normalized gene expression matrix into the composition of infiltrating immune cells.⁵⁰ We removed immune cells with zero abundance and obtained the abundance of 20 immune cell types in RA samples and controls utilizing the R package “Cibersort”. The relative proportions of the 20 immune cell types in both groups were visualized using a barplot, and the expression differences of immune cells were compared between the two groups using boxplots. The correlations were displayed among the 20 immune cell types using the R package “corrplot”. Additionally, the correlation between hub genes and the 20 immune cell types was assessed.

Clinical Data Analysis

The clinical scores, including disease activity score in 28 joints (DAS28), clinical disease activity index (CDAI), and simplified disease activity index (SDAI), were commonly used to assess disease activity for RA patients.^51–54 In this analysis, we further clarified the association between hub genes and clinical scores, including DAS28-CRP, DAS28-ESR, CDAI, and SDAI, using data from the GSE93272 dataset. The RA patients were categorized into four groups based on different scoring systems as described in Table 2. We compared the expression differences of hub genes among different groups and performed correlation analyses between hub genes and clinical indicators. Scatter plots were generated to illustrate the correlation between hub genes and DAS28-CRP, DAS28-ESR, CDAI, SDAI, MMP3 and RF.

Table 2 The Assessment of RA Disease Activity by Different Evaluation Tools

Consensus Clustering Analysis

Stratifying patients based on the disease characteristics was essential for personalized treatment of RA.⁵⁵ The “ConsensusClusterPlus” package was used to generate a consensus matrix and cumulative distribution function (CDF) plot. The relative change in the area under the CDF curve was assessed, and a tracking plot was used to determine the optimal number of clusters. Then, principal component analysis (PCA) was used to assess the ability to distinguish patients with different subtypes. The clinical markers differences between the two subtypes were compared and the results were displayed using boxplots. The visualization of differential expression of hub genes among different subtypes was achieved using a heatmap and boxplot. Finally, we analyzed the infiltration of 28 immune cells by using ssGSEA. The immune checkpoints were compared among different subtypes.

Validation of Hub Gene Expression by qPCR Analysis

This study included blood samples from 11 RA patients and 12 healthy individuals. The RA patients met the classification criteria of the ACR/EULAR 2010.⁵⁶ All participants signed a written informed consent form, and the study was approved by the Ethics Committee of the First Affiliated Hospital of Xi’an Jiaotong University (Approval No. XJTU1AF2022LSK-168).

Total RNA was extracted from the whole blood using an RNA isolation kit (5201050, Simgen, China). The RNA was then reverse transcribed into cDNA using the Reverse Transcription Kit (RR037A, Takara, Japan). The qPCR was performed using TB Green Premix Ex Taq II (RR820A, Takara, Japan). The sequences of the primers used were as follows: KLHL2 (forward primer: 5ʹ-AGAGCCGAGCAAAGAGAGTTA-3ʹ; reverse primer: 5ʹ-CCTAAGCAGTTGACAGGGTGA-3ʹ); POLK (forward primer: 5ʹ-TGAGGGACAATCCAGAATTGAAG-3ʹ; reverse primer: 5ʹ-CTGCACGAACACCAAATCTCC-3ʹ); CLEC4D (forward primer: 5ʹ-CTGATACCTTCGGTTATTGCTGT-3ʹ; reverse primer: 5ʹ-GCACTCCTGTGCCTCTCTTAC-3ʹ); NXT2 (forward primer: 5ʹ-GGACAAGGCCACCTTAATATGG-3ʹ; reverse primer: 5ʹ-TGGAACTCACTAGAAGGCAATGT-3ʹ); and 18S (forward primer: 5ʹ-GCAATTATTCCCCATGAACG-3ʹ; reverse primer: 5ʹ-GGCCTCACTAAACCATCCAA-3ʹ). The 18S gene served as an internal control, and the relative expression of genes was determined using the 2-ΔΔCT method.

Data Analysis

The bioinformatics analyses and data visualization were performed using R software (version 4.2.0). Statistical analysis was conducted using GraphPad Prism (version 9.0). Group comparisons were made using Student’s t-test for two groups, and one-way ANOVA for three or more groups (with Holm-Sidak’s multiple comparison), with statistical significance set at P < 0.05.

Results

Identification of Disulfidptosis-Related Genes in RA and Functional Enrichment Analysis

To identify significant module genes in RA, we performed WGCNA analysis using the GSE110169 dataset. The correlation heatmap revealed that the palevioletred3 module was the most significantly positively correlated with RA (r = 0.42, P = 3×10⁻⁸), while the darkorange2 module was the most significantly negatively correlated (r = −0.43, P = 2×10⁻⁸) (Figure 2A). Furthermore, the Module Membership and Gene Significance for the palevioletred3 and darkorange2 modules showed significant correlation (Supplementary Figure 2A and B). Genes were selected based on criteria of geneTraitSignificance >0.40 and geneModuleMembership >0.60 to further analyze these modules. Genes meeting these conditions were considered important for the modules. A total of 472 genes, identified as important genes from both modules, were merged for subsequent analysis.

Figure 2 Identification of disulfidptosis-related genes in RA. (A) The heatmap showed the relationship between 2 traits and 20 modules. (B) The volcano plot displayed DEGs associated with disulfidptosis. (C) The Venn plot showed the intersection of DEGs associated with disulfidptosis and significant module genes in RA. Abbreviations: DEGs, differentially expressed genes; RA, rheumatoid arthritis.

To explore the relationship between RA and disulfidptosis, we calculated the disulfidptosis score for each sample of RA. A total of 1576 disulfidptosis-related DEGs were identified comparing the high-scoring group with the low-scoring group. The volcano plot and heatmap displayed all DEGs, with 871 genes upregulated and 705 genes downregulated (Figure 2B, Supplementary Figure 3). Finally, we identified 49 important genes by intersecting the DEGs associated with disulfidptosis and the significant module genes in RA (Figure 2C).

The KEGG pathway enrichment analysis revealed the association between the intersecting genes and signaling pathways. The analysis results indicated that the intersecting genes were primarily enriched in pathways such as “Ribosome”, “Cardiac muscle contraction”, “Protein export”, “Thermogenesis”, and “Oxidative phosphorylation” (Supplementary Figure 4A). GO analysis, covering Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF), was employed to assess the functional enrichment of genes. In the BP term, significant enrichments included “amide biosynthetic process”, “protein localization to membrane”, and “establishment of protein localization to membrane” (Supplementary Figure 4B). Notable differences in CC were observed in “organelle envelope”, “envelope”, and “ribosomal subunit” (Supplementary Figure 4C). Within the MF category, prominent enrichments were seen in “structural constituent of ribosome”, “cytochrome-c oxidase activity”, and “heme-copper terminal oxidase activity” (Supplementary Figure 4D).

Identification of Hub Genes via Machine Learning Algorithms and ssGSEA

To select hub genes for subsequent model construction, we employed three machine learning algorithms for feature selection. Lasso regression results revealed 5 genes (CLEC4D, KLHL2, MRPL50, NXT2, POLK) with the lowest binomial deviation (Figure 3A). SVM-RFE results indicated 9 genes (ACAP2, ACSL4, CD46, C5orf28, CLEC4D, KLHL2, NXT2, POLK, TMEM55A) with the highest accuracy and the lowest error rate (Figure 3B). Random Forest analysis displayed the top 15 important genes, including ACAP2, CETN3, CLEC4D, COPS2, ESYT1, GIMAP2, GOPC, KLHL2, NXT2, POLK, RSL24D1, SEC62, SRP9, TMEM55A, and YEATS4 (Figure 3C). Table 3 summarizes all results. As shown in Figure 3D, 4 genes, including KLHL2, POLK, CLEC4D and NXT2, were identified by the overlapping areas of the three machine learning algorithms. Furthermore, we investigated the relationship between hub genes and hallmark gene sets. A correlation heatmap depicted that the associations of the four genes had a significant correlation with “IL6_JAK_STAT3_SIGNALING” and “WNT_BETA_CATENIN_SIGNALING” (Supplementary Figure 5). These hallmark gene sets were regarded as the signaling pathways associated with RA.

Table 3 Genes Selected by Three Machine Learning Algorithms

Figure 3 Identification of hub genes via machine learning algorithms. (A) The five genes were selected as the features with the lowest binominal deviation by the LASSO regression. (B) The nine genes were identified as the features with the lowest error rate and highest accuracy via the SVM-RFE. (C) The top fifteen genes were ranked according to their importance score based on the Random Forest. (D) The four hub genes were identified through the intersection of three machine learning algorithms. Abbreviations: CV, Cross Validation; SVM-RFE, Support Vector Machine Recursive Feature Elimination.

Figure 4 Validation of hub genes in RA. (A) The different expression of the four hub genes in RA and HC; ***: P < 0.001. (B) The hub genes were utilized in the construction of the nomogram to predict the risk of RA. (C) The predicted value of the nomogram and the hub genes using the ROC curve. (D) The calibration curve showed the accuracy of the nomogram. (E) The decision curve analysis indicated that the net benefit of the constructed model surpassed that of the default method. Abbreviations: AUC, area under the curve; HC, Healthy control; RA, rheumatoid arthritis; ROC, receiver operating characteristic.

The Model Construction Based on Hub Genes and the External Verification

To assess the value of the hub genes for RA diagnosis, we used data from three datasets: GSE110169 as the training set, and GSE89408 and GSE93272 as external validation sets, for the model construction and validation. Initially, we observed significantly higher expression of hub genes in RA compared with Healthy control (HC) across all three datasets (Figure 4A, Supplementary Figure 6A and B). Subsequently, the model was constructed and visualized using the nomogram in GSE110169. The nomogram displayed scoring criteria for each feature’s different expression levels. By calculating the total score for each sample, we predicted the probability of each sample developing RA (Figure 4B). ROC curves demonstrated that the model exhibited worthy predictive performance and the model’s predictive performance surpassed that of each hub genes (Figure 4C). Similarly, in an external validation set that included data from both whole blood and synovial biopsy samples, the model also showed promising predictive value (Supplementary Figure 6C and D). In the GSE110169 dataset, the AUC values and 95% confidence intervals for NOMOSCORE (AUC: 0.833, 95% CI: 0.768–0.899), KLHL2 (AUC: 0.813, 95% CI: 0.745–0.880), CLEC4D (AUC: 0.792, 95% CI: 0.722–0.862), NXT2 (AUC: 0.77, 95% CI: 0.696–0.843), and POLK (AUC: 0.78, 95% CI: 0.708–0.853) were shown (Supplementary Figure 7A). Supplementary Figure 7B shows that in the GSE93272 dataset, the NOMOSCORE model (AUC: 0.895, 95% CI: 0.847–0.943) demonstrated superior predictive power compared to each hub gene (Supplementary Figure 7B). Furthermore, we were somewhat surprised to find that in the GSE89408 dataset, which contained synovial biopsy samples, NOMOSCORE exhibited markedly higher predictive power (AUC: 0.977, 95% CI: 0.952–1.000) compared to each hub gene (Supplementary Figure 7C). The calibration curve indicated relatively high accuracy of the nomogram (Figure 4D). Decision curve analysis showed that the net benefit of the constructed model exceeded that of the default method (Figure 4E).

Immune Infiltration Analysis and the Relevance With Hub Genes

The imbalance of immune cells played a pivotal place in the pathogenesis of RA, which necessitated a comprehensive analysis of immune infiltration in each sample. Utilizing the Cibersort analysis, we obtained profiles of 20 immune cell types after excluding two cells with zero abundance. The barplot illustrated the relative percentages of these 20 immune cells (Figure 5A). Then, significant differences were demonstrated among 7 immune cell types, including T cells CD8, T cells CD4 naïve, T cells CD4 memory activated, Tregs, T cells gamma delta, dendritic cells activated and mast cells resting among different groups (Figure 5B). Furthermore, correlation analysis revealed significant associations between hub genes and the 20 immune cell types. For example, CLEC4D demonstrated a significant positive correlation with T cells gamma delta and a significant negative correlation with T cells CD4 naïve. Similarly, KLHL2 exhibited a significant positive correlation with neutrophils and a significant negative correlation with T cells CD8. Moreover, NXT2 displayed a significant positive correlation with T cells gamma delta and a significant negative correlation with neutrophils. Lastly, POLK showed a significant positive correlation with T cells CD4 memory activated and a significant negative correlation with neutrophils (Figure 5C). The correlation heatmap displayed notable associations (Supplementary Figure 8), including dendritic cells resting exhibiting the highest positive correlation with macrophages M1 (r = 0.74), while B cells memory showed the greatest negative correlation with B cells naïve (r = −0.5).

Figure 5 Immune infiltration analysis from RA to Control and correlation between hub genes and immune cells. (A) The relative percent of 20 immune cells in each sample. (B) The different expression of 20 immune cells between RA and control; ns: no significance; *: P < 0.05; **: P < 0.01; ***: P < 0.001. (C) Correlation analysis of immune cells with hub genes. Abbreviation: RA, rheumatoid arthritis.

The Correlation Between Hub Genes Expression and RA Disease Activity

The scoring systems, including DAS28, CDAI and SDAI, were the most used to assess disease activity in RA. The RA patients were divided for 4 groups according to disease activity. Then, the expression of hub genes was compared among the 4 groups. Overall, the expression of CLEC4D, KLHL2 and NXT2 showed an elevated trend with increasing disease activity, conversely, the expression of POLK exhibited a decreasing trend in the patients stratified by DAS28-CRP (Figure 6A) as well as DAS28-ESR, CDAI, and SDAI (Supplementary Figures 9A and B, Supplementary Figures 10). The correlation heatmaps demonstrated the association between hub genes and clinical indicators (Figure 6B, Supplementary Figure 11). We observed significant positive correlations between CLEC4D, KLHL2, NXT2 and DAS28-CRP scores (r = 0.29, P = 6.6 × 10⁻⁶; r = 0.23, P = 3.51 × 10⁻⁴; r = 0.19, P = 3.04 × 10⁻³, respectively), while POLK showed a significant negative correlation with DAS28-CRP scores (r = −0.19, P = 3.04 × 10⁻³, Figure 6C). The scatter plots also illustrated the same correlations between hub genes and DAS28-ESR, CDAI, SDAI, MMP3, and RF (Supplementary Figure 12A–E).

Figure 6 The correlation between hub genes expression and RA disease activity. (A) The expression levels of four hub genes were compared among four different disease activity stages based on DAS28.CRP scores in RA; *: P < 0.05; **: P < 0.01; ***: P < 0.001. (B) The heatmap displayed the correlation and significant difference between the expression patterns of four genes and clinical markers. (C) The scatter plots demonstrated a significant correlation between the expression of the four hub genes and DAS28.CRP. Abbreviations: CDAI, clinical disease activity index; DAS28, disease activity score in 28 joints; HAQ, health assessment questionnaire; RA, rheumatoid arthritis; SDAI, simplified disease activity index.

Two RA Subtypes Based on Hub Genes Had Different Disease Activity

In RA samples, consensus clustering analysis was performed using the expression profiles of hub genes. Based on the consensus matrix, CDF plot, relative alterations in the area under the CDF curve, and tracking plot results, we determined the two subtypes (Supplementary Figure 13A–D). PCA results demonstrated significant differences between most RA samples between the two subtypes (Figure 7A). Boxplot and heatmap revealed significant differences in the expression of hub genes between the two subtypes. Specifically, cluster B had high expression of KLHL2, CLEC4D and NXT2 and low expression of POLK, on the contrary, cluster A had high expression of POLK as well as low expression of KLHL2, CLEC4D and NXT2 (Figure 7B, Supplementary Figure 14A). Furthermore, the markers associated with RA disease activity were compared between two groups, the patients of cluster B had significantly higher levels of DAS28-CRP, DAS28-ESR, CDAI, SDAI, RF, CRP, and MMP3 than the subjects of cluster A (Figure 7C). Additionally, we observed significant differences in 13 immune checkpoints and 15 immune cells between the two subtypes, mainly in activated dendritic cell, effector memory CD8 T cell, Eosinophil, Gamma delta T cell, natural killer cell and associated activation markers (Supplementary Figure 14B, Figure 7D).

Figure 7 Identification of two subtypes based on the hub genes in RA samples. (A) The PCA plot showed that RA samples were categorized as two subtypes. (B) The different expression of the hub genes between the two subtypes. (C) The boxplots displayed the differences of clinical markers between the two subtypes. (D) The different expression of 28 immune cells between the two subtypes; ns: no significance; *: P < 0.05; **: P < 0.01; ***: P < 0.001. Abbreviations: CDAI, clinical disease activity index; DAS28, disease activity score in 28 joints; PCA, principal component analysis; RA, rheumatoid arthritis; SDAI, simplified disease activity index.

The Expression Levels of Hub Genes in the in-House Validation Cohort

To validate the differential expression of hub genes between RA patients and healthy controls, we compared the expression levels of CLEC4D, KLHL2, NXT2, and POLK in 11 RA patients and 12 healthy controls. The results showed that CLEC4D, POLK, KLHL2, and NXT2 were more highly expressed in RA patients than those in the healthy controls (Figure 8), which was consistent with validation from public database analyses.

Figure 8 The expression differences of four hub genes (KLHL2, POLK, CLEC4D, NXT2) between RA and HC; *: P < 0.05; **: P < 0.01; ***: P < 0.001. Abbreviations: HC, Healthy control; RA, rheumatoid arthritis.

Discussion

The heterogeneity of RA has posed significant challenges in the early diagnosis and optimal treatment selection for RA.⁵⁷ A comprehensive understanding of the underlying pathogenesis of RA has been instrumental in developing patient stratification strategies. Recently, disulfidptosis, a novel form of cell death triggered by intracellular accumulation of excess cysteine, has been identified as a potential mechanism in the pathogenesis of RA.^33,34 Disulfidptosis modulated the equilibrium of redox reactions and cytoskeletal stability in a variety of cells, and may be involved in immunity activation and inflammatory response, which influenced the onset and progression of RA. However, there is currently a paucity of research on disulfidptosis in RA. To the best of our knowledge, this study represented the inaugural attempt to investigate the relationship between RA and disulfidptosis using bioinformatics techniques and machine learning algorithms. It encompassed comprehensive analysis of hub genes in both peripheral blood and synovial samples and evaluated the diagnostic value of the model constructed using the hub genes for RA. Based on the differential expression of the 4 hub genes, two subtypes of RA patients were determined, and different disease activity was compared between the two subgroups, which lay the groundwork for future approaches to patient stratification.

To investigate the potential correlation between RA and disulfidptosis, we identified 49 key genes by overlapping the key module genes in RA and disulfidptosis-related DEGs using WGCNA and differential analysis. Then, we identified 4 hub genes (KLHL2, POLK, CLEC4D and NXT2) among the 49 key genes by three machine learning algorithms. The 4 hub genes were associated with the autoimmune disease by involvement in the inflammatory response as well as differentiation and development of immune cells. Specifically, KLHL2, binding to With-No-Lysine Kinases (WNKs), modulated the activity of the NLRP3 inflammasome by sensing intracellular chloride and potassium ion levels, thereby leading to the production of IL-1β and other cytokines.^58,59 Moreover, KLHL2 was considered to be involved in B cell activation by regulating the B cell receptor (BCR) signaling pathway as well as inhibition of Treg-mediated tolerance by promoting FOXP3 deacetylation and destabilization.^60,61 CLEC4D was demonstrated to initiate signal cascades that induce robust innate immune and inflammatory responses by activating spleen tyrosine kinase (Syk) and the NF-κB signal pathway.^62–64 Furthermore, CLEC4D was linked to the differentiation of Th17 cells in certain inflammatory diseases.⁶⁵ The existing literature on NXT2 in autoimmune diseases is limited. In the current study, NXT2 was found to be enriched in multiple pathways, including the IL6-JAK-STAT3 pathway, which suggested a potential role for NXT2 in the pathogenesis of RA. Contrarily, POLK, a specialized DNA polymerase, may mitigate premature aging and functional decline of T cells by facilitating DNA damage response repair, thereby suppressing inflammation.^66–69 In summary, four hub genes played a potential role in immune responses and inflammation of RA.

In this study, we performed validation using an in-house cohort to confirm the significantly higher expression of 4 hub genes in RA compared with the healthy controls, which improved the reliability of the results. To assess the diagnostic value of four hub genes, we developed a model using machine learning algorithms and performed external verification to confirm the model power for RA diagnosis. Due to prediction performance improvement, error rates decline and dependability elevation of machine learning, the model was demonstrated superior diagnostic value compared with each hub gene individually by ROC analysis, calibration curve and decision curve analysis. Synovitis is the most basic pathological manifestation of RA, as is well known.¹ In the external validation set, the model exhibited a markedly enhanced predictive capacity in synovial biopsy samples, thereby reinforcing the pivotal role of the four hub genes in RA pathogenesis. To further determine the role of hub genes in RA, we compared the expression of the 4 gene among the RA patients with different disease activity. Consistent with the biological function of the genes mentioned before, the expression of CLEC4D, KLHL2, and NXT2 exhibited an elevated trend with increasing RA disease activity, while that of POLK exhibited a decreasing trend in patients stratified by four scoring systems of disease activity assessment in RA. Based on differential expression of the 4 hub genes, we further determined the two subtypes of RA patients by consensus clustering analysis, consensus matrix, CDF plot, relative alterations in the area under the CDF curve, and tracking plot results. Among the RA patients, the cluster A subgroup exhibited high expression of POLK and low expression of KLHL2, CLEC4D, and NXT2, whereas the cluster B subgroup exhibited high expression of KLHL2, CLEC4D, and NXT2 and low expression of POLK. As anticipated, the patients in the cluster B subgroup exhibited significantly higher RA disease activity compared with those in the cluster A subgroup, which provided a new approach to RA patients stratification.

This study inevitably has some limitations. Although we demonstrated the association of the 4 hub genes with clinical markers, cell or animal model was deficient. Integrating single-cell transcriptomic analysis may enhance the reliability of our findings and provide further insights into the role of the hub genes in RA.

Conclusion

This study systematically delineated the complex correlation between disulfidptosis and RA and identified four hub genes (KLHL2, CLEC4D, NXT2, POLK). Additionally, a promising diagnostic model and stratification strategy for RA was developed to assess the risk of RA onset and RA disease activity, providing unique insights into the role of disulfidptosis in RA, which may contribute to improving early diagnosis and personalized treatment strategies.