Disease Markers
Volume 2021 (2021), Article ID 9911184, 17 pages
https://doi.org/10.1155/2021/9911184

Research Article

Identification of Disease-Specific Hub Biomarkers and Immune Infiltration in Osteoarthritis and Rheumatoid Arthritis Synovial Tissues by Bioinformatics Analysis

Fei Sun, Jian lin Zhou, Pu ji Peng, Chen Qiu, Jia rui Cao, and Hao Peng

Department of Orthopedics, Renmin Hospital of Wuhan University, Wuhan 430060, China

Correspondence should be addressed to Hao Peng; penghao8599@163.com

Received 5 March 2021; Accepted 26 April 2021; Published 17 May 2021

Academic Editor:  Shabana

Copyright © 2021 Fei Sun et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Background. Osteoarthritis (OA) and rheumatoid arthritis (RA) are well-known cause of joint disability. Although they have shown the analogous clinical features involving chronic synovitis that progresses to cartilage and bone destruction, the pathogenesis that initiates and perpetuates synovial lesions between RA and OA remains elusive. Objective. This study is aimed at identifying disease-specific hub genes, exploring immune cell infiltration, and elucidating the underlying mechanisms associated with RA and OA synovial lesion. Methods. Gene expression profiles (GSE55235, GSE55457, GSE55584, and GSE12021) were selected from Gene Expression Omnibus for analysis. Differentially expressed genes (DEGs) were identified by the “LIMMA” package in Bioconductor. The DEGs were identified by Gene Ontology (GO) and KEGG pathway analysis. A protein-protein interaction network was constructed to identify candidate hub genes by using STRING and Cytoscape. Hub genes were identified by validating from GSE12021. Furthermore, we employed the CIBERSORT website to assess immune cell infiltration between OA and RA. Finally, we explored the correlation between the levels of hub genes and relative proportion of immune cells in OA and RA. Results. We identified 68 DEGs which were mainly enriched in immune response and chemokine signaling pathway. Six hub genes with a cutoff of by ROC analysis and relative expression of were identified successfully. Compared with OA, the RA synovial tissues consisted of a higher proportion of 7 immune cells, whereas 4 immune cells were found in relatively lower proportion (). In addition, the levels of 6 hub genes were closely associated with relative proportion of 11 immune cells in OA and RA. Conclusions. We used bioinformatics analysis to identify hub genes and explored immune cell infiltration of immune microenvironment in synovial tissues. Our results should offer insights into the underlying molecular mechanisms of synovial lesion and provide potential target for immune-based therapies of OA and RA.

1. Introduction

Rheumatoid arthritis (RA) affects approximately 1% of the population, has a female : male ratio of 2.5 : 1, and mostly occurs among adults aged 40–70 years [1]. RA is characterized by synovitis, synovial inflammation and hyperplasia, inflammatory cell infiltration, cartilage and bone destruction, and other systemic features, which may lead to joint deformity, severe disability, and even premature mortality [2–5]. Osteoarthritis (OA) is the most prevalent chronic degenerative arthritis. The knee is the most frequently affected joint with 10% of men and 13% of women over the age of 60 years suffering from symptomatic knee OA and a quarter of OA patients suffering disabilities [6, 7]. Recent studies revealed that synovial lesions make an important role in pathological change of OA and RA [8, 9]. The immune microenvironment in the synovium of OA and RA is significantly different. There is no or little inflammatory cell infiltration in OA, while in RA, there is a variety of inflammatory cell infiltration and the release of abnormal cytokines [10]. These clinical features pose an important mechanistic question that needs to be urgently determined: Why does synovial inflammation perpetuate? What is the difference of synovitis between OA and RA? What causes these biological mechanisms in OA and RA?

With the development of bioinformatics technology, microarrays based on high-throughput platforms have emerged as an efficient tool to investigate gene expression profiles and explore molecular mechanisms of various diseases [11]. Synovial tissues are the important tissue in gene analysis [12]. The difference in the pathological changes and immune microenvironment between OA and RA synovial tissues may be highly related to key genes. The differentially expressed genes (DEGs) in synovial tissue will bring new insights into the pathogenesis and precise treatment target of OA and RA. However, the pathophysiologic mechanisms of arthritis in RA and OA have not been investigated thoroughly. The promising biomarkers for diagnosis, prognosis, and therapeutic response of OA and RA are lacking. Thus, the purpose of this study is to employ bioinformatics analysis to identify DEGs, screen out promising biomarkers and analysis pathways and immune infiltration patterns between OA and RA, which may shed light on the understanding of the pathogenesis, occurrence, development, and consequently provide new strategies for its diagnosis, prevention, and treatment.

2. Materials and Methods

2.1. Microarray Data Collection and Workflow

The workflow of this study is shown in Figure 1. Gene expression profiles (see S1 in Supplementary Materials) (GSE55235, GSE55457, GSE55584, and GSE12021) were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) [13, 14]. The datasets consisted of a total of 81 synovial membrane samples, including 36 OA tissue samples and 45 RA samples (Table 1). Platform annotation information was downloaded to match the gene probes with gene names, and the “sva” package of Bioconductor was applied to remove batch effects. Further, the raw datasets were preprocessed by performing background correction and normalization using R (version 4.0.2) to delete unrecognized and duplicate genes.

Figure 1: Workflow of this study.

Table 1: Information of datasets acquired from GEO.

Datasets
Platform
Total sample number
OA sample
RA sample

GSE55235
GPL96
20
10
10

GSE55457
GPL96
23
10
13

GSE55584
GPL96
16
6
10

GSE12021
GPL96
22
10
12

2.2. Identification of Differentially Expressed Genes (DEGs)

The R package linear models for microarray data (limma) have been a popular tool for determining DEGs from microarray and RNA sequence datasets [15]. Here, the R package “limma” was applied to the GEO datasets (GSE55235, GSE55457, and GSE55584) to screen genes that satisfied the following cutoff criteria: Benjamini-Hochberg adjusted value (adj. ) < 0.05 and . The value was adjusted to control the false discovery rate [16]. Genes that satisfied the criteria were designated as DEGs and were visualized as a heatmap by using the “pheatmap” package and as a volcano plot by using the “ggplot2” package of R.

2.3. Functional Annotation of the Intersection Genes

To perform Gene Ontology annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of the intersection genes, we utilized R packages “clusterProfiler” and “org.hs.eg.db” [17–19]. GO annotation enrichment terms were determined as per biological processes (BP), cellular components (CC), and molecular functions (MF) to identify the biological properties, and KEGG enrichment analysis was employed to identify the critical signaling pathways associated with the intersection genes. An adj. was set as the threshold significance level.

2.4. Construction of the Protein-Protein Interaction (PPI) Network and Identification of the Candidate Hub Genes

The Search Tool for the Retrieval of Interacting Genes (STRING; https://string-db.org) is an online database for predicting direct and indirect relationships between proteins [20]. The intersection genes were chosen to build a PPI network with a minimum interaction and hiding disconnected nodes [21]. Maximal Clique Centrality (MCC), degree, and closeness algorithms, regarded as the most effective methods to calculate hub nodes, were used to identify the candidate hub genes using the “cytoHubba” plugin in Cytoscape (version 3.7.2) [22, 23]. The intersection genes of these 3 algorithms were regarded as the candidate hub genes, which were visualized as an UpSetR diagram using the “UpSetR” package of R.

2.5. Identification of Hub Genes by Validation from GSE12021

To evaluate the accuracy of candidate genes, receiver operating characteristic (ROC) curve analysis was applied by using the R package “pROC” in GSE12021 [24]. Hub genes were identified with a validating criteria of and relatively expressed levels of with the wilcoxTest.

2.6. Determination of Immune Infiltration and Correlation of the Hub Genes and Immune Cells in OA and RA Synovial Tissues

The online analysis tool CIBERSORT (https://cibersort.stanford.edu/), a commonly used tool to assess the relative content of 22 types of immune cells, was applied to assess the immune microenvironment in the OA and RA tissues [25–28]. The proportion of these immune cells, calculated with significance criteria of value < 0.05, was visualized as a violin plot using the “vioplot” package in R. Finally, the correlation between the hub gene expression and differential proportion of immune cells in OA and RA synovial tissues was conducted by using R packages “ggExtra” and “ggpubr” (with a cutoff of ).

2.7. Statistical Analysis

The collation of original data, statistical analysis, and visualization of the statistical results was done by using the Report Language (Perl, version 5.32.1) and R software (version 4.0.2).

3. Results

3.1. Identification of DEGs

A total of 68 DEGs between OA and RA were screened out from GSE55235, GSE55457, and GSE55584 datasets, including 45 upregulated and 23 downregulated genes (see S2 in Supplementary Materials). The volcano plot and heatmap that reveal the clustering relationship of the DEGs between OA and RA samples are shown in Figures 2(a) and Figure 2(b), respectively.

(a)

(b)

Figure 2: (a) DEGs are visualized by volcano plot filtering from GSE55235, GSE55457, and GSE55584 datasets (red means upregulated genes and green means downregulated genes; and adjusted ). (b) A heatmap of partial DEGs clustering relationship between OA and RA synovial tissue.

3.2. Functional Enrichment Analyses of the Intersection Genes

In the GO analysis of the DEGs, the BP-enriched terms were humoral immune response, lymphocyte-mediated immunity, and B cell activation; the CC-enriched terms were immunoglobulin complex, external side of plasma membrane, and MHC protein complex; and the MF-enriched terms were antigen binding, immunoglobulin receptor binding, and chemokine activity (Figure 3(a)). KEGG enrichment analysis indicated that these genes were significantly enriched in the cytokine-cytokine receptor interaction, rheumatoid arthritis, and chemokine signaling pathway (Figure 3(b)).

(a)

(b)

Figure 3: GO (a) and KEGG (b) enrichment analysis of DEGs.

3.3. Construction of PPI and Identification of the Candidate Hub Genes

The protein-protein interaction network of 68 DEGs was constructed as shown in Figure 4(a). Top 10 genes were identified by Maximal Clique Centrality (MCC), degree, and closeness algorithms, respectively (see S3 in Supplementary Materials). Seven candidate hub genes from these three common algorithms were screened out including LEP, CXCL3, CXCL9, CXCL10, CXCL13, NPY1R, and POU2AF1 (Figure 4(b)).

(a)

(b)

Figure 4: (a) The protein-protein interaction (PPI) network of DEGs with hiding disconnected nodes in the network (red represents upregulated genes and green represents downregulated genes). (b) The intersection genes were screened out from 3 common algorithms (MCC, degree, and closeness).

3.4. Identification of the Hub Genes

GSE12021 was employed to validate these candidate hub genes. The result of 7 candidate genes by ROC analysis were visualized (Figure 5). Six genes (CXCL3, CXCL9, CXCL10, CXCL13, NPY1R, and POU2AF1), with AUC score more than 0.80 and relative expression levels of by the wilcoxTest, were regarded as hub genes (Figure 6). In addition, the expression levels of the six hub genes are shown in S4 in Supplementary Materials.

Figure 5: The ROC curve of 7 candidate hub genes validated by GSE12021.

Figure 6: The relative expression of 6 hub genes validated by GSE12021.

3.5. Analysis of Immune Infiltration

Using the CIBERSORT website, we calculated the relative proportion of subpopulations of different immune cells in the OA and RA synovial tissues (Figure 7(a)). Compared with the OA samples, the RA synovial tissues consisted of a higher proportion of 7 immune cells (plasma cells, CD8+ T cells, CD4+ naive T cells, activated memory CD4+ T cells, follicular helper T cells, T cells gamma delta, and M1 macrophages), whereas 4 immune cells (T cells regulatory, M2 macrophages, dendritic cells resting, and mast cells resting) were found in relatively lower proportion (Figure 7). In addition, the levels of 6 hub genes were closely associated with relative proportion of 11 immune cells in OA and RA (Figures 8 and 9). The expression of CXCL3 was positively related to resting dendritic cells in OA, whereas it was negatively correlated with resting dendritic cells in RA (Figures 8(a) and 9(a)).

(a)

(b)

Figure 7: Analysis of immune infiltration between OA and RA synovial tissues. (a) Relative proportion of 22 subpopulations of immune cells in 59 samples performed by using the CIBERSORT algorithm. (b) The difference of 22 subpopulations of immune cells among OA and RA synovial tissues (blue represents OA and red represents RA).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8: The correlation between expression levels of 6 hub genes and relative proportion of 11 immune cells in OA synovial tissues ().

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 9: The correlation between expression levels of 6 hub genes and relative proportion of 11 immune cells in RA synovial tissues ().

4. Discussion

Microarray technology and high-throughput technology, the important approaches to explore the expression levels of genes, have enabled a deeper understanding of the intrinsic molecular mechanism of complex disorders [29–31]. Increasing evidence has indicated that pathological changes in the synovium play a crucial role in the pathogenesis and progression of OA and RA [32, 33]; thus, it becomes essential to identify key genes that are differentially regulated and their target pathways. In this study, six hub genes (CXCL3, CXCL9, CXCL10, CXCL13, NPY1R, and POU2AF1), with excellent specificity and sensitivity, were successfully identified by bioinformatics analysis. CXCL3, one of the hub genes in our study and a neutrophil-activating chemokine that belongs to the CXC-type chemokine family, was found highly expressed in uterine cervical cancer(UCC) tissues. In addition, overexpression of CXCL3 could promote the tumorigenic potential of UCC. [34]. In our study, we found that the correlation of expression of CXCL3 and resting dendritic cells was diametrically opposite in OA and RA synovial tissues, which may give insights into immune microenvironment of OA and RA. C-X-C motif chemokine 10 (CXCL10), known as interferon γ-induced protein 10, was reported to increase the migration ability of inflammatory cells mediated by CXCR3. In Cxcl10 -/- and Cxcr3 -/- mice, F4/80 macrophages and CD4 T cells that infiltrated into the synovium are significantly reduced compared to those in WT mice [35]. CXCL13 is a chemokine for many immune cells including regulatory T cells, follicular T cells, and B cells. Serum CXCL13 levels in systemic sclerosis patients were abnormally elevated, reflecting that CXCL13 has a role in aberrant activation of the immune system [36, 37]. In addition, the CXCL9, CXCL10, CXCL11/CXCR3 axis was reported to regulate immune cell activation, differentiation, and migration, leading migration of immune cells to focal targets [38, 39]. Neuropeptide Y receptor Y1 (NPY1R), identified as a novel peripheral blood biomarker, may predict for the prognosis and metastasis of breast cancer patients [40]. Another research reported that increasing expression of NPY1 may promote cartilage matrix degradation and chondrocyte hypertrophy, affecting the progression of OA [41]. POU2AF1 known as a B cell transcriptional coactivator was found to have a significant correlation with disease progression by assessing gene expression in whole blood from RA patients [42]. Notably, these previous studies increase the credibility of the hub genes that were identified in our study.

Further, the functional enrichment analysis of DEGs revealed that the hub genes were mainly enriched in the immune response signaling pathway. We explored immune cell infiltration and the correlation between hub genes and immune cells in OA and RA. Compared with the OA synovium, the proportion of 7 types of immune cells in RA synovial tissue was increased and 4 types of immune cells decreased. Previous studies have revealed that immune cells may produce various soluble mediators that affect disease progression [43]. Other researchers have found that there are differential severity of histological synovitis, proportion of immune cells, varied degrees of immune cell infiltration, and differential upregulation of genes involved in B and T cell activation/function during different stages of RA [44, 45]. Are infiltration of different immune cells and their degree of infiltration related to the clinical stage of OA and RA? Is the proportion of different immune cells in the synovial membrane related to the clinical prognosis? Do these differential parameters provide guidance for the clinical treatment of OA and RA? These questions need to be investigated in future studies.

Although our study provides new hub genes as potential biomarkers between OA and RA, the limitations of our study must be acknowledged. First, our study mainly focused on data mining and analysis, and the results are rather preliminary; although useful for an initial screening, they must be validated by analytical approaches and larger groups. In addition, considering the inconsistency of the depth and breadth of the datasets from different sequencing platforms, we selected datasets from only the GPL96 platform for analysis and validation. Therefore, more data from different platforms should be collected to obtain more solid evidence to support our results. Furthermore, screening of hub genes is subject to different sample sizes and calculation methods [46]. Finally, due to lack of relevant clinical information, the expression levels of the hub genes in different clinical stages and the degree of immune cell infiltration could not be explored. We believe that addressing the above limitations may provide more insightful results.

5. Conclusions

We utilized bioinformatics analyses to screen out 68 DEGs and 6 hub genes between OA and RA synovial tissues. We found different types of immune cell infiltration between OA and RA synovial tissues, and these hub genes are highly related to the infiltration of inflammatory cells in the immune microenvironment of synovial tissues. We believe that the results of this study increase our understanding of the underlying molecular pathogenesis between OA and RA synovial tissues and consequently provide the new detection and targeting for better therapeutic modalities.

Data Availability

The datasets analyzed during the current study are available in Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

Fei Sun and Jian lin Zhou contributed equally to this article.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (No. 81672154).

Supplementary Materials

Supplementary 1. S1 (Supplementary Material 1): datasets of GSE55235, GSE55457, GSE55584, and GSE12021 downloaded from GEO.

Supplementary 2. S2 (Supplementary Material 2): list of 68 DEGs. S3 (Supplementary Material 3): list of different algorithms to identify the candidate hub genes using the “cytoHubba” plugin in Cytoscape (version 3.7.2). S4 (Supplementary Material 4): the expression levels of these hub genes (CXCL3, CXCL9, CXCL10, CXCL13, NPY1R, and POU2AF1) in OA/RA from datasets GSE55235, GSE55457, and GSE55584.

References

J. S. Smolen, D. Aletaha, and I. B. McInnes, “Rheumatoid arthritis,” The Lancet, vol. 388, no. 10055, pp. 2023–2038, 2016.

T. L. Barker and T. L. Puckett, “Rheumatoid arthritis: coping with disability,” Rehabilitation Nursing, vol. 35, no. 2, pp. 75–79, 2010.

D. Aletaha, T. Neogi, A. J. Silman et al., “2010 Rheumatoid arthritis classification criteria: an American College of Rheumatology/European League Against Rheumatism collaborative initiative,” Arthritis & rheumatism, vol. 62, no. 9, pp. 2569–2581, 2010.

D. L. Scott, D. P. Symmons, B. L. Coulton, and A. J. Popert, “Long-term outcome of treating rheumatoid arthritis: results after 20 years,” Lancet, vol. 1, no. 8542, pp. 1108–1111, 1987.

D. M. Mitchell, P. W. Spitz, D. Y. Young, D. A. Bloch, D. J. McShane, and J. F. Fries, “Survival, prognosis, and causes of death in rheumatoid arthritis,” Arthritis and Rheumatism, vol. 29, no. 6, pp. 706–714, 1986.

Y. Zhang and J. M. Jordan, “Epidemiology of osteoarthritis,” Clinics in Geriatric Medicine, vol. 26, no. 3, pp. 355–369, 2010.

M. H. Li, R. Xiao, J. B. Li, and Q. Zhu, “Regenerative approaches for cartilage repair in the treatment of osteoarthritis,” Osteoarthritis and Cartilage, vol. 25, no. 10, pp. 1577–1587, 2017.

B. J. de Lange-Brokaar, A. Ioan-Facsinay, E. Yusuf et al., “Association of pain in knee osteoarthritis with distinct patterns of synovitis,” Arthritis & Rhematology, vol. 67, no. 3, pp. 733–740, 2015.

I. B. McInnes and G. Schett, “The pathogenesis of rheumatoid arthritis,” The New England Journal of Medicine, vol. 365, no. 23, pp. 2205–2219, 2011.

I. Prieto-Potin, R. Largo, J. A. Roman-Blas, G. Herrero-Beaumont, and D. A. Walsh, “Characterization of multinucleated giant cells in synovium and subchondral bone in knee osteoarthritis and rheumatoid arthritis,” BMC Musculoskeletal Disorders, vol. 16, no. 1, 2015.

G. R. Oliver, S. N. Hart, and E. W. Klee, “Bioinformatics for clinical next generation sequencing,” Clinical Chemistry, vol. 61, no. 1, pp. 124–135, 2015.

U. Ungethuem, T. Haeupl, H. Witt et al., “Molecular signatures and new candidates to target the pathogenesis of rheumatoid arthritis,” Physiological Genomics, vol. 42A, no. 4, pp. 267–282, 2010.

R. Edgar, M. Domrachev, and A. E. Lash, “Gene Expression Omnibus: NCBI gene expression and hybridization array data repository,” Nucleic Acids Research, vol. 30, no. 1, pp. 207–210, 2002.

T. Barrett, S. E. Wilhite, P. Ledoux et al., “NCBI GEO: archive for functional genomics data sets--update,” Nucleic Acids Research, vol. 41, pp. D991–D995, 2013.

M. E. Ritchie, B. Phipson, D. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, p. e47, 2015.

M. J. Zhang, F. Xia, and J. Zou, “Fast and covariate-adaptive method amplifies detection power in large-scale multiple hypothesis testing,” Nature Communications, vol. 10, no. 1, p. 3433, 2019.

M. Kanehisa, M. Furumichi, M. Tanabe, Y. Sato, and K. Morishima, “KEGG: new perspectives on genomes, pathways, diseases and drugs,” Nucleic Acids Research, vol. 45, no. D1, pp. D353–D361, 2017.

M. Kanehisa, Y. Sato, M. Kawashima, M. Furumichi, and M. Tanabe, “KEGG as a reference resource for gene and protein annotation,” Nucleic Acids Research, vol. 44, no. D1, pp. D457–D462, 2016.

B. Nota, “Gogadget: an R package for interpretation and visualization of GO enrichment results,” Molecular informatics, vol. 36, no. 5-6, 2017.

A. Franceschini, D. Szklarczyk, S. Frankild et al., “STRING v9.1: protein-protein interaction networks, with increased coverage and integration,” Nucleic Acids Research, vol. 41, pp. D808–D815, 2013.

H. V. Cook, N. T. Doncheva, D. Szklarczyk, C. von Mering, and L. J. Jensen, “Viruses.STRING: a virus-host protein-protein interaction database,” Viruses, vol. 10, no. 10, p. 519, 2018.

R. Chaudhary, M. Balhara, D. K. Jangir, M. Dangi, M. Dangi, and A. K. Chhillar, “In silico protein interaction network analysis of virulence proteins associated with invasive Aspergillosis for drug discovery,” Current Topics in Medicinal Chemistry, vol. 19, no. 2, pp. 146–155, 2019.

H. Yang, Y. Wang, Z. Zhang, and H. Li, “Identification of KIF18B as a hub candidate gene in the metastasis of clear cell renal cell carcinoma by weighted gene co-expression network analysis,” Frontiers in Genetics, vol. 11, p. 905, 2020.

X. Robin, N. Turck, A. Hainard et al., “pROC: an open-source package for R and S+ to analyze and compare ROC curves,” BMC Bioinformatics, vol. 12, no. 1, 2011.

W. Jin, Y. Zhang, Z. Liu, Z. Che, M. Gao, and H. Peng, “Exploration of the molecular characteristics of the tumor-immune interaction and the development of an individualized immune prognostic signature for neuroblastoma,” Journal of Cellular Physiology, vol. 236, no. 1, pp. 294–308, 2021.

J. Falconer, A. N. Murphy, S. P. Young et al., “Review: synovial cell metabolism and chronic inflammation in rheumatoid arthritis,” Arthritis & Rhematology, vol. 70, no. 7, pp. 984–999, 2018.

A. M. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015.

R. Zhou, J. Zhang, D. Zeng et al., “Immune cell infiltration as a biomarker for the diagnosis and prognosis of stage I-III colon cancer,” Cancer Immunology, Immunotherapy, vol. 68, no. 3, pp. 433–442, 2019.

H. Yan, G. Zheng, J. Qu et al., “Identification of key candidate genes and pathways in multiple myeloma by integrated bioinformatics analysis,” Journal of Cellular Physiology, vol. 234, no. 12, pp. 23785–23797, 2019.

S. Udhaya Kumar, D. Thirumal Kumar, R. Siva et al., “Dysregulation of signaling pathways due to differentially expressed genes from the B-cell transcriptomes of systemic lupus erythematosus patients - a bioinformatics approach,” Frontiers in Bioengineering and Biotechnology, vol. 8, 2020.

S. Mishra, M. I. Shah, S. Udhaya Kumar et al., “Network analysis of transcriptomics data for the prediction and prioritization of membrane-associated biomarkers for idiopathic pulmonary fibrosis (IPF) by bioinformatics approach,” Advances in Protein Chemistry and Structural Biology, vol. 123, pp. 241–273, 2021.

B. J. de Lange-Brokaar, A. Ioan-Facsinay, G. J. van Osch et al., “Synovial inflammation, immune cells and their cytokines in osteoarthritis: a review,” Osteoarthritis and Cartilage, vol. 20, no. 12, pp. 1484–1499, 2012.

S. Rosengren, H. M. Hoffman, W. Bugbee, and D. L. Boyle, “Expression and regulation of cryopyrin and related proteins in rheumatoid arthritis synovium,” Annals of the Rheumatic Diseases, vol. 64, no. 5, pp. 708–714, 2005.

Y. L. Qi, Y. Li, X. X. Man et al., “CXCL3 overexpression promotes the tumorigenic potential of uterine cervical cancer cells via the MAPK/ERK pathway,” Journal of Cellular Physiology, vol. 235, no. 5, pp. 4756–4765, 2020.

J. H. Lee, B. Kim, W. J. Jin, H. H. Kim, H. Ha, and Z. H. Lee, “Pathogenic roles of CXCL10 signaling through CXCR3 and TLR4 in macrophages and T cells: relevance for arthritis,” Arthritis Research & Therapy, vol. 19, no. 1, p. 163, 2017.

H. S. Carlsen, E. S. Baekkevold, H. C. Morton, G. Haraldsen, and P. Brandtzaeg, “Monocyte-like and mature macrophages produce CXCL13 (B cell-attracting chemokine 1) in inflammatory lesions with lymphoid neogenesis,” Blood, vol. 104, no. 10, pp. 3021–3027, 2004.

T. Taniguchi, T. Miyagawa, S. Toyama et al., “CXCL13 produced by macrophages due to Fli1 deficiency may contribute to the development of tissue fibrosis, vasculopathy and immune activation in systemic sclerosis,” Experimental Dermatology, vol. 27, no. 9, pp. 1030–1037, 2018.

M. F. Hsieh, S. L. Lai, J. P. Chen et al., “Both CXCR3 and CXCL10/IFN-inducible protein 10 are required for resistance to primary infection by dengue virus,” Journal of Immunology, vol. 177, no. 3, pp. 1855–1863, 2006.

R. Tokunaga, W. Zhang, M. Naseem et al., “CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation - a target for novel cancer therapy,” Cancer Treatment Reviews, vol. 63, pp. 40–47, 2018.

L. Liu, Q. Xu, L. Cheng et al., “NPY1R is a novel peripheral blood marker predictive of metastasis and prognosis in breast cancer patients,” Oncology Letters, vol. 9, no. 2, pp. 891–896, 2015.

X. Kang, Z. Qian, J. Liu et al., “Neuropeptide Y acts directly on cartilage homeostasis and exacerbates progression of osteoarthritis through NPY2R,” Journal of Bone and Mineral Research, vol. 35, no. 7, pp. 1375–1384, 2020.

M. F. Romo-Garcia, H. S. Nava-Ramirez, M. Zapata-Zuniga et al., “Evaluation of SUMO1 and POU2AF1 in whole blood from rheumatoid arthritis patients and at risk relatives,” International Journal of Immunogenetics, vol. 46, no. 2, pp. 59–66, 2019.

S. Viatte, D. Plant, and S. Raychaudhuri, “Genetics and epigenetics of rheumatoid arthritis,” Nature Reviews Rheumatology, vol. 9, no. 3, pp. 141–153, 2013.

G. Lliso-Ribera, F. Humby, M. Lewis et al., “Synovial tissue signatures enhance clinical classification and prognostic/treatment response algorithms in early inflammatory arthritis and predict requirement for subsequent biological therapy: results from the pathobiology of early arthritis cohort (PEAC),” Annals of the Rheumatic Diseases, vol. 78, no. 12, pp. 1642–1652, 2019.

F. Rivellese, F. Humby, S. Bugatti et al., “B cell synovitis and clinical phenotypes in rheumatoid arthritis: relationship to disease stages and drug exposure,” Arthritis & Rhematology, vol. 72, no. 5, pp. 714–725, 2020.

W. C. Li, L. Bai, Y. Xu et al., “Identification of differentially expressed genes in synovial tissue of rheumatoid arthritis and osteoarthritis in patients,” Journal of Cellular Biochemistry, vol. 120, no. 3, pp. 4533–4544, 2019.