Screening and identification of potential hub genes and immune cell infiltration in the synovial tissue of rheumatoid arthritis by bioinformatic approach

Abstract

Background: Rheumatoid arthritis (RA) is an autoimmune disease that affects individuals of all ages. The basic pathological manifestations are synovial inflammation, pannus formation, and erosion of articular cartilage, bone destruction will eventually lead to joint deformities and loss of function. However, the specific molecular mechanisms of synovitis tissue in RA are still unclear. Therefore, this study aimed to screen and explore the potential hub genes and immune cell infiltration in RA.

Methods: Three microarray datasets (GSE12021, GSE55457, and GSE55235), from the Gene Expression Omnibus (GEO) database, have been analyzed to explore the potential hub genes and immune cell infiltration in RA. First, the LIMMA package was used to screen the differentially expression genes (DEGs) after removing the batch effect. Then the clusterProfiler package was used to perform functional enrichment analyses. Second, through weighted coexpression network analysis (WGCNA), the key module was identified in the coexpression network of the gene set. Third, the protein-protein interaction (PPI) network was constructed through STRING website and the module analysis was performed using Cytoscape software. Fourth, the CIBERSORT and ssGSEA algorithm were used to analyze the immune status of RA and healthy synovial tissue, and the associations between immune cell infiltration and RA-related diagnostic biomarkers were evaluated. Fifth, we used the quantitative reverse transcription-polymerase chain reaction (qRT-PCR) to validate the expression levels of the hub genes, and ROC curve analysis of hub genes for discriminating between RA and healthy tissue. Finally, the gene-drug interaction network was constructed using DrugCentral database, and identification of drug molecules based on hub genes using the Drug Signature Database (DSigDB) by Enrichr.

Results: A total of 679 DEGs were identified, containing 270 downregulated genes and 409 upregulated genes. DEGs were primarily enriched in immune response and chemokine signaling pathways, according to functional enrichment analysis of DEGs. WGCNA explored the co-expression network of the gene set and identified key modules, the blue module was selected as the key module associated with RA. Seven hub genes are identified when PPI network and WGCNA core modules are intersected. Immune infiltration analysis using CIBERSORT and ssGSEA algorithms revealed that multiple types of immune infiltration were found to be upregulated in RA tissue compared to normal tissue. Furthermore, the levels of 7 hub genes were closely related to the relative proportions of multiple immune cells in RA. The results of the qRT-PCR demonstrated that the relative expression levels of 6 hub genes (CD27, LCK, CD2, GZMB, IL7R, and IL2RG) were up-regulated in RA synovial tissue, compared with normal tissue. Simultaneously, ROC curves indicated that the above 6 hub genes had strong biomarker potential for RA (AUC >0.8).

Conclusions: Through bioinformatics analysis and qRT-PCR experiment, our study ultimately discovered 6 hub genes (CD27, LCK, CD2, GZMB, IL7R, and IL2RG) that closely related to RA. These findings may provide valuable direction for future RA clinical diagnosis, treatment, and associated research.

Keywords: Bioinformatics approach; Hub genes; Immune infiltration; Rheumatoid arthritis; WGCNA.

Fig. 1

A flowchart of the studies that were considered for inclusion in the analysis.

Fig. 2

(A) PCA plot before removing batch effects. (B) PCA plot after removing batch effects. Clusters of sample points with different distances are from different batches and sequencing platforms. Whereas plot B shows a reduced difference in distance between batches after removing the batch effect. (C) Fold-change values and an adjusted P-value were used to create volcano graphs. The over-expressed mRNAs are shown by the red point in the figure, while the down-expressed mRNAs are represented by the blue point. (D) According to logFC, the top 20 downregulated genes and upregulated genes that were differentially expressed between RA and normal tissues are represented in the heatmap. The over-expressed mRNAs are shown by the red point in the figure, while the down-expressed mRNAs are represented by the blue point. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3

(A) The enriched KEGG signaling pathways were chosen to show the important biological activities of significant potential mRNA. The abscissa represents the gene ratio, whereas the ordinate depicts the enriched pathways. (B) Analysis of probable mRNA targets using the Gene Ontology (GO) system. The ClusterProfiler utility in R software was used to cluster prospective targets by biological pathways (BP), molecular function (MF), and cellular component (CC). q value < 0.05 was judged to have statistical significance in the enrichment result. (C) The top five GSEA enrichment analysis of DEGs. DEGs are for differentially expressed genes; GO stands for Gene Ontology; KEGG is for Kyoto Encyclopedia of Genes and Genomes; GSEA stands for gene set enrichment analysis.

Fig. 4

From merged datasets, coexpression modules are identified by WGCNA. (A) Calculation of the soft thresholding value for scale-free coexpression networks. (B) Cluster dendrogram of the modules identified. (C) Gene interactions in coexpression modules. (D) Module correlations with RA patients and controls.

Fig. 5

Construction of the PPI network and module analysis of DEGs between Rheumatoid arthritis and normal controls. (A) The most significant module was extracted from the PPI network through the MCODE plugin of Cytoscape. Upregulated genes are marked in red; downregulated genes are marked in green. (B)The hub genes were extracted from the PPI network through Degree algorithms of cytoHubba. (C) MCODE plugin, Degree algorithms of cytoHubba and hub genes of MMbrown module obtain a Venn diagram of common hub genes. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6

The immune infiltration landscape in rheumatoid arthritis and healthy tissue. (A) The proportion of 22 immune cell subpopulations in 64 samples from the GSE12021, GSE55235, and GSE55457 datasets. (B) All samples were subjected to principal component analysis. The first two major components, which account for the majority of data fluctuation, are displayed. (C) The difference between rheumatoid arthritis and healthy controls in terms of immune infiltration. (The normal controls group was color-coded blue, whereas the rheumatoid arthritis group was color-coded red. Statistical significance was defined as a P value < 0.05). (D) Positive and negative correlations among 22 immune cell types, with red and blue indicating positive and negative correlations, respectively. The absence of any association between the mentioned immune cell types is represented by the color white. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7

Immune gene sets identified by ssGSEA. Based on the ssGSEA algorithm, 16 immune cells and 13 immune-related functions were calculated in rheumatoid arthritis and healthy tissue. (∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ns, non-significant).

Fig. 8

The correlation between hub genes and the immune status in RA (A) The proportion of infiltrated immune cells was obtained by CIBERSORT algorithm analysis. (B) The proportion of infiltrated immune cells was obtained by ssGSEA algorithm analysis. Positive or negative correlations between the hub genes and 22 immune cell types, with red and blue indicating positive and negative correlations, respectively. The absence of any association between the mentioned immune cell types is represented by the color white. (∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 9

Validation of the hub genes in RA and normal tissue using qRT-PCR. All studies were carried out in triplicate, with the findings provided as M ± SD. (∗ P < 0.05, ∗∗ P < 0.01,∗∗∗P < 0.001, ns, non-significant).

Fig. 10

ROC curve for the 6 hub genes that are specifically expressed.

Fig. 11

Gene-drug interaction network constructed from hub genes and drugs based on DrugCentral. Available drugs that increase or decrease mRNA or protein expression of hub genes are shown in panels A and B. Hub genes are marked in red, drugs increase the expression of hub genes are marked in blue; drugs decrease the expression of hub genes are marked in green. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 12

The role of targeted drugs to antagonize hub genes expression in the pathogenesis of RA. Targeted drugs antagonize hub genes expression to inhibit immune cell infiltration, thereby inhibiting synovial activation and ultimately inhibiting the occurrence of RA.