We retrieved the GEO platform (http://www.ncbi.nlm.nih.gov/geo/) for eligible RA data from whole blood. We used “Rheumatoid arthritis”, “Homo sapiens”, and “Expression profiling by array” as search keywords, and suitable datasets were screened based on the following criteria: the dataset includes a minimum of 80 samples comprised of downloadable raw data and series matrix files, with at least 40 samples each in the control and RA groups. After careful screening, we selected the dataset GSE93272 (14), which fully meet our criteria with 232 RA cases and 43 controls.

We adopted annotation package (R4.1.2) from Bioconductor (http://bioconductor.org/) to transform microarray probes into gene symbols. Then, the data was standardized through quantile normalization, including 232 RA samples and 43 controls. The FAM-related genes were retrieved and identified using the GeneCards database (https://www.genecards.org/) with “fatty acid metabolism” as a keyword. Totally 104 FAM-related genes (shown in Supplementary Table 1) were screened with a relevance score ≥50 (15) for the subsequent analysis. We used Limma package to identify differentially expressed FAM regulators between controls and RA patients. We screened the significant FAM regulators according to screening thresholds of |log₂ fold change (FC)|>0 and P-Value<0.05 (16). Then, the R package “clusterProfiler” was used to perform GO and KEGG enrichment analysis to explore the underlying mechanism of the FAM regulators implicated in RA. Moreover, we constructed the protein-protein interaction (PPI) network of these FAM regulators by the STRING database (https://string-db.org/), conducted network topology analysis and screened the top ten targets as FAM hub genes in terms of degree through Cytoscape software (v3.8.0).

Two machine learning algorithms including random forest (RF) and support vector machine (SVM) models were adopted to predict the occurrence of RA. The vital FAM modulators were screened in virtue of the R package “RandomForest” when their significance scores (Mean Decrease Gini) were greater than 2. In the SVM model, the variable n signifies the count of FAM hub genes, with each data point depicted as a dot within an n-dimensional space. We then selected an optimal hyperplane that distinctly separated the control and RA groups (17). Subsequently, the “rms” R package was employed to develop a nomogram model for predicting the prevalence of RA patients based on the identified candidate FAM regulators. Calibration curves assessed the accuracy of the prediction values against actual outcomes. Decision curve analysis (DCA) was conducted to generate a clinical impact curve and evaluate whether model-based decisions were advantageous for patients (18).

Through consensus clustering with resampling, each member and its corresponding subcluster number were identified, demonstrating the validity of the clusters (18). Using the “ConsensusClusterPlus” R package, different FAM patterns were identified based on FAM hub genes (19).

Differentially expressed genes (DEGs) between different FAM clusters were identified using the Limma package, applying a threshold of adjusted P-Value <0.05 and |log₂FC| >0.585. GO analysis was then conducted with the “clusterProfiler” R package to explore the involvement of DEGs in the process of RA (20).

To quantify the FAM clusters, principal component analysis (PCA) was adopted to assess the FAM score for each sample. This score was calculated using the following formula: FAM score = PC1_i, where PC1 represents principal component 1, and i indicates distinct FAM gene expression (21).

Single sample gene set enrichment analysis (ssGSEA) was used to quantify the levels of immune infiltration in RA group. Initially, the gene expression levels in the samples were ranked through sequencing using ssGSEA. Subsequently, we examined the input dataset for FAM hub genes and compiled their expression levels. From this analysis, we determined the quantity of immune cells present in each sample (22).

The Ethics Committee of Laboratory Animals in Guangdong Provincial Hospital of Chinese Medicine approved all studies. Female Sprague–Dawley(SD) rats, aged 8 weeks and weighing 200–220g, were purchased from the Experimental Animal Center of Guangzhou University of Chinese Medicine (Guangzhou, China). They were maintained under standard environmental conditions (22 ± 2°C, 50% humidity, and a 12-h light/dark cycle) with unrestricted access to food and water. The rats were euthanized under isoflurane anesthesia.

All animal experiments were approved by the Ethics Committee of Laboratory Animals in Guangdong Provincial Hospital of Chinese Medicine (No. 2023081) and conducted in accordance with the relevant guidelines. The study was carried out in compliance with the ARRIVE guidelines.

To isolate BMMs, we flushed long bones from 8-week-old rats using warm, serum-free alpha-minimum essential medium (α-MEM). The isolated BMMs were cultured with M-CSF (100 ng/mL) for 2 days to recruit macrophages, followed by the addition of RANKL (50 ng/mL) to induce osteoclast differentiation. RNA-seq analysis was then performed to examine the differential expression of FAM-related genes between groups with and without RANKL induction during osteoclast differentiation. Libraries from different samples were pooled according to quantitative assessments, and the final data were used for sequencing. DEGs were identified by comparing control and RANKL-induced samples using the Limma R package. FAM modulators were subsequently identified, and their expression profiles were established based on the data. The criteria for detecting FAM DEGs were set at P < 0.05.

To evaluate the relationships among significant FAM genes, linear regression analyses were used. Group comparisons in the bioinformatics analysis were conducted with Kruskal-Wallis tests, and corrected t-tests were applied to assess RNA-seq data. All parametric tests were two-tailed, with P<0.05 deemed statistically significant. Results are shown as mean ± standard deviation.

We totally screened 53 distinct FAM regulators through difference analysis of gene expression profiles between RA group and the controls (Figure 2A). Our analysis revealed that GO enrichment predominantly identified entries related to biological processes (notably fatty acid metabolic process), cellular components (specifically peroxisomal matrix), and molecular function (including lipid transporter activity) (Figure 2B). Moreover, KEGG pathway enrichment analysis uncovered that PPAR signaling pathway and fatty acid metabolism were notably significant pathways (Figure 2C). The PPI network of 53 distinct FAM regulators was plotted in Figure 2D. We ultimately screened 10 FAM hub genes (SREBF1, SCD, PPARG, PPARA, INSR, FASN, CD36, ACADVL, ACADM, ACACA), which were shown in Figure 2E. We observed that the expressions of ACADM, CD36, PPARG were upregulated in RA samples in comparison with controls, but the other FAM hub genes showed opposite outcomes (Figures 2F–O).

To elucidate the potential correlations among significant FAM genes in RA patients, Pearson correlation analysis was conducted utilizing R statistical software. FAM hub genes in RA exhibited different relationships with each other (Figure 3A). Thereafter, the remarkable correlations with R>|0.25| were selected for visualization. Significantly positive relationships were observed between the gene expressions of ACACA-PPARA, ACADVL-INSR, ACADVL-PPARA in RA cases (Figures 3B–D), but gene expressions of ACADM-ACADVL, ACADM-INSR, ACADM-SREBF1, CD36-FASN in RA cases showed significantly negative relationships (Figures 3E–H).

The RF model was validated to have the smaller residual according to reverse cumulative distribution of residual (Figure 4A) and boxplots of residual (Figure 4B). Most of the model samples have relatively small residuals, which indicates that the RF model is superior to the SVM model. Moreover, we utilized ROC curves to evaluate the models, and according to their AUC values, we discovered that the RF model exhibited higher accuracy than the SVM model (Figure 4C). As a result, we came to the conclusion that the RF model is the best one for predicting the occurrence of RA. Finally, we presented these 10 FAM hub genes based on their importance score (mean decrease Gini) and selected candidate genes with importance score>2, including SREBF1, SCD, PPARG, PPARA, INSR, FASN, CD36, ACADVL, ACADM, ACACA (Figure 4D).

To predict the prevalence of RA patients, a nomogram model was constructed using the “rms” package in R based on 10 candidate FAM regulators (Figure 5A). The calibration curves indicated high prediction accuracy of the nomogram model (Figure 5B), and the DCA curve suggested potential benefits for RA patient judgments using this model (Figure 5C). Furthermore, the clinical impact curve demonstrated remarkable predictive capacity of the nomogram model (Figure 5D).

Two FAM clusters (clusterA and clusterB) were identified on the basis of the ten FAM hub genes in virtue of the R package “ConsensusClusterPlus” (Figures 6A–D). Cluster A consisted of 159 samples, while cluster B included 73 samples. Subsequently, the heat map and boxplot clearly showed the differential expression levels of the 10 important FAM regulators between the two clusters. We observed that clusterA exhibited higher expression levels of CD36 and ACADM compared to clusterB, whereas SREBF1, PPARA, FASN, and ACADVL showed higher expression levels in clusterB than in clusterA. The expression levels of SCD, PPARG, INSR, and ACACA did not exhibit any noticeable variances between the two clusters (Figures 6E, F). The 10 FAM regulators were able to distinguish between the two FAM clusters based on the PCA results (Figure 6G). We identified 74 DEGs associated with FAM between the two FAM patterns. To gain further insight into the role of these DEGs in RA, we conducted GO enrichment analysis (Figure 6H). We observed that GO:0002181 (cytoplasmic translation), GO:0003735 (structural constituent of ribosome) and GO:0005840 (ribosome) were the mainly enriched entries.

We then explored the relationship between immune cells and 10 important FAM regulators by using ssGSEA to assess the abundance of immune cells in RA samples. We observed a positive association between INSR and multiple immune cells (Figure 7A). We compared the differences in immune cell infiltration between patients with high and low INSR expressions. Our results showed that patients with high INSR expression had significantly increased immune cell infiltration compared to those with low INSR expression (Figure 7B). Furthermore, we found that clusterA cases were associated with activated CD4 T cell, activated CD8 T cell, eosinophil, Gamma delta T cell, immature dendritic cell, MDSC, macrophage, regulatory T cell, and Type 2 T helper cell immunity; while clusterB was linked to CD56dim natural killer cell, Natural killer T cell, T follicular helper cell, Type 1 T helper cell immunity (Figure 7C).

Based on the 74 FAM-associated DEGs, we used a consensus clustering technique to classify RA cases into different genomic subtypes in order to understand FAM patterns. We identified two distinct FAM gene clusters (gene clusters A and B) that aligned with the sectionalization of FAM patterns (Figures 8A–D). The expression levels of the 74 FAM-related DEGs in gene cluster A and gene cluster B were depicted in Figure 8E. Similarities in immune cell infiltration levels and expressions of 10 significant FAM modulators between gene clusterA and gene clusterB also mirrored those in the FAM clusters (Figures 8F, G). These results once again confirmed the accuracy of our sectionalization using the consensus clustering method.

The Sankey diagram (Figure 9A) illustrated the connection between FAM scores, FAM clusters, and FAM gene clusters. PCA methods were employed to measure the FAM clusters by determining the FAM scores for each sample across the two distinct FAM clusters. In comparison to clusterA, we observed that clusterB showed a higher FAM score (Figure 9B). To explore the associations between FAM clusters and RA, we assessed the relationships between FAM clusters and three specific genes including IL17RA, TBXA2R, and RXRA, which have close association with osteoclast differentiation. We found that clusterB exhibited higher levels of IL17RA, TBXA2R, and RXRA expression than clusterA, indicating that clusterB may be strongly connected with RA defined by osteoclast differentiation (Figure 9C).

The expression heat map (Figure 10A) showed the differential expression profiles during osteoclast differentiation. Specifically, the FAM modulator CD36 exhibited increased expression levels in RANKL-induced group compared with controls, while the FAM modulators SREBF1, FASN, SCD1 and SCD2 exhibited decreased expression levels in RANKL-induced group compared with controls (Figure 10B), which verified the bioinformatics results.