Identification of diagnostic genes and vital microRNAs involved in rheumatoid arthritis: based on data mining and experimental verification

Data download and processing

Gene expression matrices and the clinical information of GSE55457 based on platform GPL96 [HG-U133A], were downloaded from the GEO database. Twenty-three samples were selected for our study. We performed normalization and log2 transformation on original expression matrices and used a platform annotation file to convert probes into gene symbols. If a gene was detected by more than one probe, the average expression value was used for subsequent analysis.

Differentially expressed genes (DEGs) between RA patients and healthy controls were analyzed by limma, which is a tailored R package for expression profiling studies (Ritchie et al., 2015). The DEGs filtering thresholds were set to |logFC|>1 and the adjusted p-value was < 0.05.

Weighted correlation network analysis

Gene co-expression analysis of dataset GSE55457 was performed in R software (version: 3.6.3) based on the WGCNA R package reference manual. We calculated the median absolute deviation (MAD) for each gene, which reflected the degree of variation in gene expression across samples, and selected the top 3,000 highly variable genes for network construction and modules detection. First, we clustered samples to detect the presence of outliers. Next, optimal soft-thresholding power (β) was selected, which was a pivotal step in the entire analysis process. In order to reduce interference from spurious correlations, the adjacency matrix was replaced by a topological overlap matrix (TOM), and then the hclust function was used to generate a hierarchical clustering tree (dendrogram) of genes. The detailed parameters used for modules detection were as follows: minModuleSize was set to 30 to obtain large modules; modules with a correlation of 0.75 were merged using mergeCloseModules function. We also quantified module–trait associations and their p-values, which allowed us to pick out modules that were closely related to the trait of interest for next step analysis. The absolute value of correlation coefficient > 0.6 and p-value < 0.05 were cutoffs for module screening. Finally, we measured the gene significance (GS) and module membership (MM) of individual genes and scatterplots were drawn.

Functional annotation and GSEA

R package clusterProfiler (Yu et al., 2012) and its dependency annotation package org.Hs.eg.db, sourced from Bioconductor, were used to analyze and visualize the functional profiles (GO and KEGG) of key module genes. To obtain highly credible analysis results, only enrichment items satisfying both p-value < 0.05 and q-value < 0.05 were outputted. Moreover, gene set enrichment analysis (GSEA) (Mootha et al., 2003; Subramanian et al., 2005), which can be used to determine whether an already-defined gene set exhibits a statistically significant difference in two different traits (e.g., RA patients and healthy individuals), was conducted to explore biological processes and pathways that may be involved in the pathogenesis of RA. The complete gene expression matrix and phenotype information of dataset GSE55457 were uploaded to GSEA and analyzed with hallmark gene sets that are stored in molecular signatures database (MSigDB) (Liberzon et al., 2015; Liberzon et al., 2011). Gene sets enriched at FDR < 0.25 and p-value < 0.01 were retained.

Building PPI network and mining hub genes

The intersection of DEGs and key module genes was entered into the STRING database to find interrelationships (Szklarczyk et al., 2015; Szklarczyk et al., 2019; Szklarczyk et al., 2017). Due to the large size of the initial protein-protein interaction (PPI) network, we hid isolated nodes and limited the minimum interaction score to 0.7. Then, the network was exported and further optimized in Cytoscape software (Cline et al., 2007; Otasek et al., 2019; Shannon et al., 2003). The MCODE plugin was installed in Cytoscape to find clusters made up of highly related genes in the network (Bader & Hogue, 2003). The parameters used in our analysis were default values. Genes in the highest scoring cluster were considered to be hub genes of RA and were re-performed for KEGG analysis to explore vital biological pathways.

Estimating proportion of immune cells in synovial tissue

Understanding the immune microenvironment in RA may provide new insights into the immunotherapy of the disease. The CIBERSORT tool used for estimating the immune cell contents of different samples in a gene expression admixture has a wide range of applicability and high sensitivity, which is important for disease diagnosis and therapeutic target predictions (Newman et al., 2015). A mixture file was generated according to the input format requirements, and we selected default LM22 as signature gene file to distinguish the hematopoietic populations and activation states. The analysis results were filtered based on a threshold of p-value < 0.05 to remove samples with a low confidence level. We conducted correlation analysis between the hub genes and immune cell subtypes significantly enriched in RA samples according to Pearson’s algorithm.

Prediction of gene-miRNA relationships

MiRNAs are a class of non-coding RNAs that regulate gene expression by making mRNA unstable and inhibiting translation, a process by which mRNA is converted into protein. Here, MCODE-identified hub genes were uploaded to the online analysis tool miRDB; we searched for miRNAs that have regulatory roles in the development of RA (Chen & Wang, 2020; Liu & Wang, 2019). The parameters were set as follows: only functional miRNAs were considered; gene targets with a score below 60 and miRNAs with more than 2,000 targets were excluded. Cytoscape was then used to construct a gene-miRNA network in which miRNAs targeting two or more genes were considered to be vital regulators involved in RA pathogenesis.

Validating hub genes and miRNAs with independent microarray datasets

The expression matrix GSE12021 was selected as the validation dataset for hub genes. We plotted ROC curves in R software based on the pROC package to test whether screened hub genes have good diagnostic value for RA (Robin et al., 2011). The area under the curve (AUC) was computed according to the trapezoidal rule and genes with AUC > 85% were further validated by RT-qPCR. To understand the differences in transcriptional regulators between RA patients and healthy individuals, we validated vital miRNAs with GEO dataset GSE37425, which focused on the expression of miRNAs in synovial tissue using the TaqMan probe approach. The entire gene expression matrix was normalized and then transformed with log2.

Validating hub genes and miRNAs with RT-qPCR

Five RA patients and three healthy individuals who attended the health examination center of Xiaoshan Traditional Chinese Medicine Hospital were included and all subjects signed consent forms. This study was approved by Ethics Committees of Xiaoshan Traditional Chinese Medicine Hospital (NO. 2020012). We collected blood samples from each test subject and isolated the peripheral blood mononuclear cells (PBMCs) using human lymphocyte separation medium (Solarbio, China). Next, the total RNA of PBMCs was extracted using the TRIzol reagent (Invitrogen, USA) and its purity was assessed by NanoDrop (Thermo Scientific, USA). We used a cDNA synthesis kit (Thermo Scientific, USA) to generate the first strand cDNA and cDNA amplification was performed on a 7500 real-time PCR instrument (Applied Biosystems, USA). For the quantitative analysis of miRNA expression, Hairpin-it miRNAs qPCR quantitation kit (Genepharma, China) was applied. We chose GAPDH and U6 as internal reference genes and compared relative mRNA or miRNA expression based on the 2^−ΔΔCT method. The primer sequences are available in Table 1.

Western blot

Total proteins were extracted from PBMCs using RIPA lysis buffer (Beyotime Biotechnology, China), followed by concentration detection with the BCA protein assay kit (Sangon Biotech, China). Equal amounts of protein were separated by 8–20% SDS-PAGE and transferred onto polyvinylidene difluoride membranes. Afterwards, membranes were blocked with 5% skim milk for 1 h and then incubated with primary antibodies including β-actin (1:1000, Cell Signaling Technology (CST), #3700), and CXCL10 (1:1000, CST, #14969) at 4 °C overnight. Anti-mouse IgG (1:2000, CST, #7076) and anti-rabbit IgG (1:2000, CST, #7074) were used to bind the primary antibodies. Finally, membranes were scanned using FluorChem Q (Proteinsimple, USA). Image J software was used for the quantitative analysis of protein bands.

Statistical analysis

Statistical analysis was performed in GraphPad Prism 7.0. Data satisfying normal distribution were analyzed by unpaired t-test; the Mann–Whitney test was used for data that were not normalized. P-value < 0.05 was considered significant.

Screening for DEGs

Our differential expression assessment of dataset GSE55457 identified 367 DEGs, of which 186 genes were highly expressed and 181 genes had low expression levels in RA patients. We created a volcano map and heat map for DEGs as shown in Fig. 2.

Identification of key module closely related to RA

Two samples (GSM1337307 and GSM1337308) were identified as outliers and were excluded by the cutreeStatic function with parameter cutHeight = 90 (Fig. 3A). The soft-thresholding power (β) was set to 5, for which the scale-free topology fit index (R²) reached 0.90 (Fig. 3B). Co-expression analysis resulted in 14 modules, and the size of each module is shown in Figs. 3C and 3D. Figure 3E showed the relationships between modules and traits, with four modules (black, turquoise, blue, and purple) having strong positive or negative correlations with disease status (presence or absence of RA). Taking into account gene numbers, correlation coefficients, and p-values of the above four modules, we focused on the turquoise module containing 845 genes and performed an in-depth analysis. The GS and MM plot of the turquoise module (correlation coefficient = 0.79, p-value = 2.7e−181) demonstrated that the closer the gene is to the trait, the more important it is in the module (Fig. 3F).

Functional annotation of turquoise module genes and GSEA

The top 10 results for GO enrichment are shown in Fig. 4A and mainly involved leukocyte activation and cell adhesion. To associate a specific gene with its assigned GO terms, we adopted the GOChord function to generate a circle diagram (Fig. 4B). Figure 4C pointed out that inflammation-related signaling pathways, such as chemokine signaling pathway, NF-kappa B signaling pathway, Th17 cell differentiation and T cell receptor signaling pathway, were significantly linked to the pathogenesis of RA. Similarly, a circle plot was displayed to visualize the correspondence between genes and pathways (Fig. 4D). We identified five gene sets based on a reference database derived from MSigDB using the GSEA method. These were associated with the immune response or inflammatory pathways that were highly enriched in RA, which verified the GO and KEGG results (Figs. 5A–5E).

Constructing PPI network to find hub genes

A total of 171 intersection genes were entered into the STRING database to establish the network (Fig. 6A). The network was adjusted using Cytoscape software, resulting in 52 nodes and 169 relationship pairs. The MCODE plugin uncovered five clusters, with the highest scoring cluster containing nine nodes, which were potential hub genes for RA (Fig. 6B, Table 2). Furthermore, we re-performed KEGG analysis on the nine genes to investigate the signaling pathways in which they were involved. As shown in Fig. 6C, we found that nine hub genes played important roles in inflammatory or immune-related pathways, including the chemokine signaling pathway, toll-like receptor signaling pathway, and TNF signaling pathway.

Immune infiltration and Pearson correlation analysis

One sample (GSM1337307) was excluded from analysis results because its p-value did not meet the threshold we set. The contents of various immune cells in each sample were presented by stacked bar chart (Fig. 7A). M2 macrophages accounted for the highest proportion in both RA and normal tissue. We conducted a Wilcoxon test to further clarify whether the level of each immune cell type was statistically different between two groups. From the violin plot (Fig. 7B), we inferred that five cell types, including plasma cells, CD8 T cells, follicular helper T cells, γδ T cells, and M1 macrophages, were abundant in RA synovial tissue compared to the normal group. More importantly, Pearson correlation analysis confirmed that, in addition to CCL25, eight other hub genes identified by MCODE played pro-inflammatory roles. Their overexpression may recruit large numbers of inflammatory cells to lesion sites, which is a critical pathological basis for RA (Fig. 7C).

Predicting miRNAs that regulate nine hub genes

Normally, miRNAs can fine-tune gene expression to maintain the homeostasis of life, while abnormal miRNAs expression is thought to be associated with diseases, including benign or malignant tumors, and immune system dysfunction. Our prediction of miRNAs was made in the miRDB database and 72 miRNA-gene regulatory relationships were obtained. We then created an interactive network in Cytoscape and analyzed the number of genes targeted by each miRNA (Fig. 8). The miRNAs possessing two target genes are listed in Table 3 and were vital regulators in the pathogenesis of RA.

Validation of hub genes by ROC curves

The expression matrix of dataset GSE12021 was chosen to test the diagnostic efficacy of hub biomarkers for RA. Probability curves (ROC curves) were plotted and the AUC reflected the degree or measure of separability. A higher AUC indicated that the gene is better at discriminating between diseased and non-diseased individuals. As shown in Fig. 9, six genes had AUC values greater than 85%, including CCL25, CXCL10, CCL5, CXCL9, CXCL13, and CXCL11, indicating that these genes were potential biomarkers for the diagnosis of RA.

Validation of hub genes by RT-qPCR and Western blot

We collected peripheral blood from subjects for RT-qPCR validation to further evaluate the six genes with AUC > 85% identified by ROC curves. Five genes were differentially expressed between RA patients and normal subjects (Fig. 10A). Among them, CXCL9, CXCL10, CXCL11, and CXCL13 were up-regulated genes, while CCL25 was down-regulated. Moreover, to detect the expression of genes in protein level, we chose a representative gene (CXCL10) for Western blot validation. As shown in Figs. 10B and 10C, the protein and mRNA level of gene expression were consistent.

Vital miRNAs verification

The vital miRNAs possessing two target genes were validated with GSE37425. Among the six miRNAs investigated, has-miR-449b, has-miR-573, and has-miR-203 expression was suppressed in RA patients (Figs. 11A–11F). It is worth noting that only the differences in has-miR-573 and has-miR-34a were statistically significant. Previous studies have reported that miR-34a played an important role in the pathogenesis of RA, but the findings were controversial (Niederer et al., 2012; Xie et al., 2019). Hence, we analyzed the expression of miR-34a-5p in PBMCs using RT-qPCR and found that miR-34a-5p was significantly down-regulated in RA patients compared to health controls (Fig. 11G).