The search term “Rheumatoid Arthritis” was queried in the GEO to obtain expression profile data of genes related to RA from three datasets GSE55235, GSE55457, and GSE77298 based on 10 normal samples and 10 RA samples, 10 normal and 13 RA samples, and 7 normal and 16 RA samples, respectively.

The raw gene expression profile data of all three RA GEO datasets were downloaded and subjected to pre-processing, such as standardization, correction, and gene name annotation using the limma package (version 3.42.2). The data were combined, and the sva package (version 3.36.0) was used to perform batch-removal correction on all chip data. DEGs in RA were screened according to $|$log${}_2$ fold change (FC)$|$$>$1.0 and adjusted p value 0.05. The ggplot2 package (version 3.3.1) and pheatmap package (version 1.0.12) were used to plot the heat map and volcano map corresponding to the screened DEGs. All the data analysis was done in R language (version 4.0.2, https://www.r-project.org/).

The clusterProfiler package (version 3.14.3) was used for Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of all DEGs according to p 0.05, q 0.05, and Cytoscape (version 3.7.2, https://cytoscape.org/) software was used to construct an RA gene-pathway network.

All DEGs identified in RA were imported into the STRING (version 11.0, https://string-db.org/) database to construct a PPI network of RA targets, with the parameters set as the following: organism, Homo sapiens; combined score threshold, 0.7. The PPI network for RA was plotted using Gephi (version 0.9.2, https://gephi.org/).

To identify key RA genes in the PPI network, we used two methods for considering the network node degree, which differ from previous methods in that they focus on examining the efficiency of nodes in the disease network.

The maximum connected reduction (MCR${}_i$) refers to the relative maximum connected reduction of the network after deletion of node v${}_i$ in the network, which can be used as an index for measuring the importance of the deleted node, calculated as follows:

(1)$$MC{R}_i=\frac{N-N_i}{N}$$

where N represents the total number of nodes connected to the network before node deletion, and N${}_i$ represents the number of nodes in the maximum connected branch after deletion of node v${}_i$.

Similarly, the efficiency relative reduction (ERR${}_i$) of node v${}_i$ in the network is calculated as:

(2)

(3)$$ER{R}_i=\left|\frac{E-E_i}{E}\right|$$

where E is the average value of the reciprocal of the shortest distance between network nodes, d${}_{i\mathit{}j}$ is the shortest distance between two nodes in the network, E${}_i$ is the efficiency of the remaining networks after deleting node v${}_i$ and k${}_i$ edges connected to v${}_i$. Thus, ERR${}_i$ represents the relative reduction of network efficiency after deleting node v${}_i$; the larger the ERR${}_i$, the greater the disturbance degree of node v${}_i$ to the PPI network.

We included the top genes in terms of the above two indicators into the RA key gene set based on PPI network efficiency analysis.

According to the sample classification of the normal group and RA group, we used the RFECV function in the sklearn$$package (version 0.24.2) and adopted three different machine learning algorithms, namely least absolute shrinkage and selection operator (lasso) regression, random forest (RF) classifier, and linear support vector machine (SVM) classification, to perform feature extraction on the gene expression data of RA samples to screen out the key genes based on classification features. The gene sets screened out by these three different machine learning classifiers were intersected, and the characteristic genes common to more than two classifiers were selected as the key RA gene set based on the feature extraction method.

The key RA gene sets based on PPI network analysis and feature extraction were intersected again to obtain the core gene set of RA. To verify the disease classification effect of these core genes, the three different machine learning classifiers mentioned above were used to predict RA classification, which was evaluated by the area under the receiver operating characteristic (ROC) curve (AUC) value. All classification prediction analyses were performed on the Anaconda3 (https://www.anaconda.com/).

Based on GEO gene expression profile data, CIBERSORT Toolkit (version 1.03) was used to analyze immune cell infiltration in RA, and the abundances of 22 immune cell types in RA were calculated. The abundance difference heat map, correlation heat map, and violin plot of immune cells in each sample were plotted using the pheatmap (version 1.0.12), corrplot (version 0.84), and vioplot (version 0.3.5) packages, respectively. Spearman correlation analysis was used to analyze the correlations between 22 types of infiltrating immune cells and key genes to identify genes closely correlated with immune cell infiltration.

The chemical information and data of each component in our previously compiled active ingredient set of natural traditional Chinese medicines proven to have clear pharmacological effects, including 432 small molecules [11], were obtained through the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), including the English name, CAS, PubChem CID, molecular formula, canonical Simplified Molecular-Input Line-Entry System (SMILES), and SDF file.

Before molecular docking of TCM small molecules according to known target receptors, the receptor file containing the original ligands was downloaded from the PDB database (https://www.rcsb.org/); molecules with a similar structure to the original ligand were considered as the potential and optimal candidates for further screening using molecular fingerprint similarity.

According to canonical SMILES, the molecular fingerprint of each compound was calculated using DRAGON (version 7.0, Kode, Pisa, Italy) software, and the similarity between the original ligands and TCM small molecules was calculated by the Tanimoto coefficient:

(4)$$T(A,B)=\frac{c}{a+b-c}$$

Among them, a is the number of descriptors of the basic fragment in compound A, b is the number of descriptors of the basic fragment in compound B, and c is the number of descriptors of the basic fragment shared by compounds A and B.

Further, the molecular fingerprint similarity coefficient $S(A,B)$ between any pair of original ligand and TCM small molecule was obtained. Taking the overall Tanimoto coefficient value $T(A,B)$ as the background similarity value, the Z value was calculated according to Eqn. 5 to measure and screen whether the molecular fingerprint similarity between the original ligand and TCM small molecule was significantly greater than the mean value of overall background similarity, as follows:

(5)$$Z=\frac{T(A,B)-\overline{S}(A,B)}{sd(S(A,B))}$$

where $\overline{S}(A,B)$ is the overall mean and $sd(S(A,B))$ is the overall standard deviation.

The ligand SDF file of all collected TCM small molecules was downloaded from the PubChem, wherein the three-dimensional (3D) structure file was directly converted into a PDB file using PyMOL (version 1.6, https://pymol.org/2/) and the two-dimensional file was converted into a PDB file in 3D format using Chem3D (version 19.0, PerkinElmer, Waltham, Massachusetts, USA). The PDB file of target receptors was downloaded from the PDB database, and key target proteins were pre-processed using PyMOL to remove impurities such as water molecules. The PDB file of receptors and ligands was then converted to a PDBQT file using AutoDockTools (version 1.5.6, https://autodock.scripps.edu/).

The coordinates of the central position for molecular docking were determined by referring to the binding site of the protein receptor and the original ligand, and the docking center parameter was defined as a box size of 30 $\times$ 30 $\times$ 30. AutoDock Vina (version 1.1.2, https://vina.scripps.edu/) software was used for semi-flexible molecular docking, the affinity values (kcal/mol) of all small molecules with the key targets were calculated, and judgment and analysis were based on the affinity values and docking conformations.

After the combination and batch removal correction of the GSE55235, GSE55457, and GSE77298 datasets, a total of 527 DEGs in RA were screened out, including 348 up-regulated genes and 179 down-regulated genes. The heat map and volcano plot for the expression levels of the DEGs are shown in Fig. 2.

Fig. 2.

DEGs analysis results of RA. Heat map (A) and volcano plot (B) of expression analysis of 527 differentially expressed genes in rheumatoid arthritis (RA) compared with normal control (NC) samples after combination and batch removal of three groups of GEO chip datasets (GSE55235, GSE55457, and GSE77298).

Functional annotation and enrichment analysis suggested that the DEGs in RA are mainly involved in various immune-related functions and pathways. Among the 1233 associated GO biological functions, the DEGs were mainly concentrated in 24 GO functions, as shown in Fig. 3A., namely leukocyte migration (GO:0050900), T cell activation (GO:0042110), immune response-activating cell surface receptor signaling pathway (GO:0002429), immune response-activating signal transduction (GO:0002757), etc.

Fig. 3.

Functional annotation and enrichment analysis results of DEGs in RA. (A) GO functional enrichment analysis of differentially expressed genes in RA. (B) KEGG pathway enrichment analysis of differentially expressed genes in RA. (C) Gene-pathway enrichment network; the dots represent genes and the squares represent pathways.

KEGG pathway enrichment analysis suggested that the DEGs are mainly involved in 39 signaling pathways, including cytokine-cytokine receptor interaction (hsa04060), chemokine signaling pathway (hsa04062), Epstein-Barr virus infection (hsa05169), viral protein interaction with cytokine and cytokine receptor (hsa04061), and cell adhesion molecules (hsa04514); the top 30 main signaling pathways with enrichment of the highest number of genes are shown in Fig. 3B and the gene-enrichment pathway network is shown in Fig. 3C.

We imported the 527 DEGs into the STRING to obtain the PPI network of RA, as shown in Fig. 4A. Two types of importance degree values of all genes in the PPI network were obtained using the MCR and ERR algorithms, respectively, as shown in Fig. 4B.

Fig. 4.

Identification key RA genes by PPI network. (A) Protein-protein interaction network of RA. (B) Screening of key RA genes based on the MCR (top) and ERR (bottom) algorithms.

According to the calculation results, genes that met the threshold conditions (MCR$>$0.0065 and ERR$>$0.005) were selected, the intersection of which was taken, thus obtaining 31 key genes based on PPI network analysis, as shown in Fig. 5A: PCK1, RHOB, AGTR1, MYH11, CFP, VAMP8, PLA2G7, PTPRC, SEMA4D, STAT1, LY96, ATF3, SCD, ADCY2, IGJ, PLA2G2D, MCM5, EGFR, ITGA4, MYC, LPL, IRS2, ZAP70, LILRB2, CD48, GBP1, APOB, RAC2, EGR1, CSK, and BIRC3.

Fig. 5.

Identification key RA genes by machine learning. (A) Results of key RA gene screening based on the lasso regression model, random forest classifier (RFC), and support vector classifier (SVC) feature extraction, and RA core gene discovery based on PPI network analysis and machine learning feature extraction. (B) Receiver operating characteristic (ROC) curves of disease prediction classification of the eight core genes.

Using the RFECV tool to extract features and combining with each of the three machine learning algorithms, 486 genes related to classification were screened out by the lasso regression model, 111 were obtained by the RF classification model, and 35 were extracted by the linear SVM model. After taking the intersection of the genes extracted by the three classification models, a total of 134 genes with classification characteristics were obtained (Fig. 5A). After further intersecting with the core gene set obtained by PPI network analysis, eight common genes were found: MYH11, CFP, LY96, IGJ, LPL, CD48, RAC2, and CSK (Table 1).

Further, we took these eight genes as the classification feature indicators of data samples, used the three machine learning classification models (Lasso, RF, and SVM) in the sklearn package again for 50% cross-validation, and calculated and plotted AUC values and ROC curves of the three classifier models, as shown in Fig. 5B. The average classification prediction AUC value of these eight genes in the three classifiers was 0.63, and the highest value was 0.69. Even though the sample size of RA data included in this analysis was relatively small, the AUC values of the three classification models were all greater than 0.6, indicating that these eight core genes do have certain classification potential for the existing RA samples, suggesting that these genes have a certain correlation with RA.

The results of immune cell infiltration analysis are shown in Fig. 6. Among the 22 types of immune cells in RA, correlations were found between neutrophils and activated dendritic cells ($\rho$ = 0.58), neutrophils and resting natural killer (NK) cells ($\rho$ = 0.59), and resting NK cells and CD4+ naïve T cells ($\rho$ = 0.50); however, the correlation coefficients of the infiltration degrees between most other immune cells were generally low.

Fig. 6.

Results of immune cell infiltration evaluation and gene correlation analysis. (A) Infiltration and abundance of 22 types of immune cells in all RA samples. (B) Heat map of the infiltration and abundance of 22 types of immune cells in RA. (C) Correlation heat map of 22 types of immune cells in RA. (D) Differences in the infiltration of 22 types of immune cells between RA and normal samples.

Among the 22 types of immune cells, the infiltration of memory B cells, plasma cells, CD8+ T cells, CD4+ resting memory T cells, CD4+ activated memory T cells, follicular helper T cells, gamma delta T cells, activated NK cells, M1 macrophages, resting dendritic cells, activated mast cells, and eosinophils was significantly different between the normal group and RA group (p 0.05).

Spearman correlation analysis showed that six out of the eight key RA genes were significantly correlated with the infiltration and abundance of immune cells (p 0.05), as shown in Table 2. Among them, CD48 was correlated with infiltration of five types of immune cells, namely plasma cells ($\rho$ = 0.74), CD4+ naïve T cells ($\rho$ = 0.50), CD4+ activated memory T cells ($\rho$ = 0.64), gamma delta T cells ($\rho$ = 0.59), and activated NK cells ($\rho$ = –0.57). IGJ was correlated with infiltration of four types of immune cells, namely memory B cells ($\rho$ = 0.51), plasma cells ($\rho$ = 0.91), CD4+ activated memory T cells ($\rho$ = 0.62), and activated NK cells ($\rho$ = –0.53). RAC2 was correlated with infiltration of five types of immune cells, namely plasma cells ($\rho$ = 0.78), CD4+ resting memory T cells ($\rho$ = –0.54), CD4+ activated memory T cells ($\rho$ = 0.56), follicular helper T cells ($\rho$ = 0.59), and activated NK cells ($\rho$ = –0.57). LY96 and MYH11 were also correlated with plasma cells and CD4+ activated memory T cells, respectively, and CSK was correlated with plasma cells. Therefore, these six genes, namely CD48, CSK, IGJ, LY96, MYH11, and RAC2, are potential therapeutic markers related to immune cell infiltration in RA.

Moreover, the five types of immune cells correlated with RAC2 showed significant differences between the normal and RA samples. Therefore, we selected RAC2 protein for the subsequent screening and verification of TCM small molecules.

A PDB file (PDB ID: 2W2V) [12] of RAC2 protein was downloaded from the PDB database. All four chains of RAC2 protein had the known ligand guanosine-5’-triphosphate (GTP). Therefore, we screened for TCM small molecules with high similarity to GTP in terms of the chemical structure.

We collected one SMILES string of GTP (PubChem CID: 135398633) from the PubChem. The TCM small molecules and the SMILES strings of GTP were imported into DRAGON software to calculate the 1024-dimensional extended-connectivity molecular fingerprints of all small molecules.

A total of 31 TCM small molecules with similar structures to GTP were screened out according to a threshold of Z$>$ 1.0 for subsequent molecular docking screening.

Since RAC2 protein has four chains (A, B, C, and D) and each chain has a corresponding ligand-binding site, we performed batch molecular docking on four sites in the four chains of RAC2. Among them, the A-chain binding site coordinate is (5.711, –4.855, 41.735), the B-chain binding site coordinate is (–10.673, –5.506, 23.999), and the C-chain binding site coordinate is (10.654, 9.806, –11.671), the D-chain binding site coordinate is (14.897, 58.007, 29.835).

According to the scores of docking the 31 small molecules with the molecules on the different chains of RAC2, the top 14 TCM small molecules with average affinity values less than –7.0 kcal/mol and with relatively good scores were as follows: jatrophon, sanguinarine, miroestrol, withaferin A, pristimerin, chelidonine, picroside I, sesamin, gnidimacrin, cephalotaxine, honokiol, oleuropein, atractylenolide I, and dracorhodin.

Finally, based on the molecular docking results and the reported literature, we selected four TCM small molecular compounds for subsequent experimental verification: sanguinarine, chelidonine, sesamin, and honokiol. It should be noted that since Jatrophon was out of stock, we chosen sanguinarine with the second-ranked molecular docking score. The other three components were all reported small molecules of TCM with potential to treat RA (Table 3, Ref. [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]).

The 2D binding conformations of these 4 small molecules and RAC2 were displayed using LigPlot (version 2.2, https://www.ebi.ac.uk/thornton-srv/software/LIGPLOT/), as shown in Fig. 7.

Fig. 7.

The molecular docking results of these 4 small molecules and RAC2. The 2D conformation diagram of (A) Sanguinarine, (B) Chelidonine, (C) Sesamin, and (D) Honokiol in the 4 chains of RAC2, from left to right shown the ligand binding poses in chain A, chain B, chain C, and chain D.

The in vitro CCK-8 assays showed that none of the screened out TCM compounds had a significant effect on the proliferation of normal rat synovial cells; however, the proliferation of RA rat synovial cells to sanguinarine, chelidonine, sesamin, and honokiol was significantly decreased, with IC${}_{50}$ values of 344.3 $\pm$ 3.2 $\mu$g/mL, 511.7 $\pm$ 4.7 $\mu$g/mL, 756.8 $\pm$ 7.0 $\mu$g/mL, and 345.1 $\pm$ 3.2 $\mu$g/mL, respectively (Supplementary Material: The curves used to determine IC${}_{50}$ of these 4 small molecules).

Compared with the normal control group, the levels of TNF-$\alpha$ and IL-1$\beta$ produced by rat RA synovial cells significantly increased, whereas the levels of TNF-$\alpha$ and IL-1$\beta$ in cell the culture supernatants of the sanguinarine-, sesamin-, and honokiol-treated cells significantly decreased, suggesting that these compounds could significantly inhibit the hypersecretion of TNF-$\alpha$ and IL-1$\beta$ in rat RA synovial cells (Fig. 8A).

Fig. 8.

Verification with RA rat model synovial cells. (A) Effects of different small-molecule compounds of TCM on the secretion levels of TNF-$\alpha$ and IL-1$\beta$ in synovial cells of collagen-induced arthritis rats. *p 0.05, **p 0.01, ***p 0.001 compared with the model group. (B) Western blots of RAC2 protein in different groups: Lane 1, sanguinarine treatment; Lane 2, chelidonine treatment; Lane 3, sesamin treatment; Lane 4, honokiol treatment; Lane 5, Tripterygium glycosides treatment (control); Lane 6, RA model (untreated control); Lane 7, normal rat cells. (C) Gray value analysis of WB expression of RAC2 in synovial cells of collagen-induced arthritis rats between RA and different compounds (#p 0.05, ##p 0.01, ###p 0.001).

Western blotting confirmed that the expression of RAC2 protein in RA rat synovial cells was significantly higher than that in normal synovial cells; however, treatment with sanguinarine, sesamin, and honokiol could significantly suppress the expression of RAC2 protein (Fig. 8B), and the significant differences in the WB gray values of the RA group and the three small molecules were shown in Fig. 8C.