Comprehensive Bioinformatics Analysis Reveals Hub Genes and Inflammation State of Rheumatoid Arthritis

Abstract

Rheumatoid arthritis (RA) is an autoimmune disease characterized by erosive arthritis, which has not been thoroughly cured yet, and standardized treatment is helpful for alleviating clinical symptoms. Here, various bioinformatics analysis tools were comprehensively utilized, aiming to identify critical biomarkers and possible pathogenesis of RA. Three gene expression datasets profiled by microarray were obtained from GEO database. Dataset GSE55235 and GSE55457 were merged for subsequent analyses. We identified differentially expressed genes (DEGs) in RStudio with limma package, performing functional enrichment analysis based on GSEA software and clusterProfiler package. Next, protein-protein interaction (PPI) network was set up through STRING database and Cytoscape. Moreover, CIBERSORT website was used to assess the inflammatory state of RA. Finally, we validated the candidate hub genes with dataset GSE77298. As a result, we identified 106 DEGs (72 upregulated and 34 downregulated genes). Through GO, KEGG, and GSEA analysis, we found that DEGs were mainly involved in immune response and inflammatory signaling pathway. With the help of Cytoscape software and MCODE plug-in, the most prominent subnetwork was screened out, containing 14 genes and 45 edges. For ROC curve analysis, eight genes with AUC >0.80 were considered as hub genes of RA. In conclusion, compared with healthy controls, the DEGs and their closely related biological functions were analyzed, and we held that chemokines and immune cells infiltration promote the progression of rheumatoid arthritis. Targeting the eight biomarkers we identified may be useful for the diagnosis and treatment of rheumatoid arthritis.

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.

1. Introduction

Rheumatoid arthritis (RA) is an autoimmune disease, mainly destroying synovium and joints, characterized by autoantibodies that target immunoglobulin G (known as rheumatoid factor, RF) and citrullinated proteins (called anticitrullinated protein antibodies, ACPAs) [1]. Some epidemiological studies conducted in western countries showed that the prevalence of rheumatoid arthritis is about 0.5-1.0% [2, 3]. Rheumatoid arthritis is a complicated disease due to the changeable clinical manifestations and complications in different patients or disease stages, which brings difficulties to the clinical work of doctors.

The serological detection of autoantibodies is a crucial indicator in the diagnosis and prognosis of rheumatoid arthritis, but about 25% of patients are seronegative and thus may experience a delay in diagnosis as well as initiation of drug therapy [4]. Moreover, it was estimated that 50% of seropositive patients had negative serum test results at the beginning of the disease [5]. Previous studies have shown that proinflammatory cytokines in inflammatory synovium, such as interleukin-8, can stimulate osteoclasts proliferation and then result in bone resorption of RA patients [6–8]. However, some scholars found that bone destruction may also occur in ACPA-positive individuals without detectable inflammation conditions [9]. A recent study supporting the latter result demonstrated that monoclonal ACPAs derived from B cells in the synovial fluid of RA patients have obvious epitope specificity, which promotes the differentiation of osteoclasts in cell cultures [10]. Although the pathogenetic insights, classification criteria, and therapeutic strategies of RA have been updated in the past 20 years, some patients are still unable to achieve satisfactory clinical remission or have serious adverse reactions to antirheumatoid therapy, so more efforts are required to address these unmet needs.

The microarray technology has emerged for more than 20 years, which makes it possible to analyze the complete transcriptional information of various cell types and tissues [11]. Studies based on gene expression analysis have obtained new findings, elucidating how the transcriptome varies among distinct phenotypes and stages of disease [12, 13]. The Gene Expression Omnibus (GEO) is a user-friendly repository, in which stores microarray, next-generation sequencing, and other forms of genomics data for users to query and download. Here, we aimed to dissect biomarkers and inflammation state of rheumatoid arthritis by comprehensively applying multiple bioinformatics analysis tools including R packages from Bioconductor, STRING database, CIBERSORT website, Cytoscape, and GSEA software. The findings in our study may contribute to novel ideas for better diagnosis and treatment of rheumatoid arthritis.

2. Materials and Methods

2.1. Data Download and Processing

Three microarray datasets GSE55235, GSE55457 [14], and GSE77298 [15] were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo). A total of 20 normal synovial tissues and 23 diseased specimens were enrolled from dataset GSE55235 and GSE55457, whose detection platforms were identical (GPL96, HG-U133A). The dataset GSE77298 was based on GPL570 platform (HG-U133 Plus 2), containing 7 synovium samples from healthy controls (HC) and 16 from RA patients. According to the research plan, the former two datasets were merged as training dataset to explore hub genes, and mRNA profiles of GSE77298 were used to assess whether the discovered hub genes have excellent diagnostic value for RA.

Data processing was divided into four steps. First, the three probe expression matrix files ($\ast$series_matrix.txt) downloaded from GEO database were normalized and log2 transformed. Next, we matched the platform annotation file with each probe expression matrix and well-annotated probes were retained. For multiple probes corresponding to one gene, the average expression value was taken for further analysis. Third, we merged the expression matrix of GSE55235 and GSE55457 into one, and the order of samples in dataset was rearranged. Last, R package sva, installed from Bioconductor (https://bioconductor.org/), was applied to eliminate the heterogeneity caused by different experimental batches and platforms.

2.2. Identification of DEGs

Differentially expressed genes (DEGs) were screened out by limma package [16] based on the comparison of expression values between HC samples and RA samples. The screening criteria for DEGs were as follows: log2 fold change (FC) should be greater than 2 or less than -2 and adjusted $p$ value <0.05. The analysis results were presented by heatmap and volcano map drawn in RStudio software (version:1.2.1335).

2.3. Functional Enrichment Analysis

The Bioconductor package clusterProfiler [17] was applied to carry out Gene Ontology (GO) and KEGG pathway analysis for DEGs. Based on the threshold $p$ value <0.05, GO terms and signal pathways with significant enrichment were screened out. Gene Set Enrichment Analysis (GSEA) is a software codeveloped by UC San Diego and Broad Institute [18, 19], which assesses whether a predefined gene set shows statistically significant differences between two biological phenotypes (e.g., HC and RA). Expression dataset collapsed to gene symbol and phenotype information were uploaded to GSEA for enrichment analysis with default parameters. Hallmark gene sets (h.all.v7.1.symbols.gmt) selected in the present study were downloaded from Molecular Signatures Database (MSigDB) [20, 21]. Enrichment results with nominal $p$ value <0.01 as well as FDR <0.25 were considered statistically significant.

2.4. Immune Infiltration Analysis

Online analysis tool CIBERSORT(https://cibersort.stanford.edu/), a method for calculating the cell proportion of complex tissues based on gene expression profiles, is superior to other methods in terms of large-scale analysis of RNA mixtures [22]. In this study, we used CIBERSORT to characterize the inflammation state of RA and healthy joint tissue with default signature gene file (22 types of immune cells). The analysis result was filtered according to $p$ value <0.05, and the immune cell composition of each sample was shown in barplot.

2.5. PPI Network Construction and Module Analysis

In order to explore the mutual relationship between proteins encoded by different genes, DEGs were imported into STRING website (version:11.0) for further analysis [23]. The lowest interaction score should be greater than 0.4 and isolated nodes in the network were removed. Next, we output the analysis results to a TSV format file and used Cytoscape software (version:3.7.1) for details processing and module analysis. MCODE [24] is a plug-in downloaded from Cytoscape App Store, which can find closely connected nodes in a complex network based on topology. Therefore, we applied this plug-in to detect critical modules in PPI network with default parameters.

2.6. Verification of Hub Genes by ROC Analysis

Genes in the most significant module identified by MCODE plug-in were selected as candidate hub genes. To evaluate the role of candidate genes in the diagnosis of RA, receiver operating characteristic (ROC) curve analysis was conducted in RStudio with pROC package [25]. The genes with area under curve (AUC) >0.8 as well as $p$ value <0.05 were considered as hub genes of RA.

3. Results

3.1. Identification of DEGs

By using limma package to analyze the differential expression of the integrated dataset, we obtained 106 DEGs composed of 72 upregulated and 34 downregulated genes. DEGs screened by threshold were visualized by volcano map (Figure 1(a)). The expression of top 25 upregulated and downregulated genes ordered by adjusted $p$ value in each sample was shown by heatmap (Figure 1(b)).

[figures omitted; refer to PDF]

3.2. GO and KEGG Enrichment Analysis

Enrichment analysis conducted by clusterProfiler package revealed the biological functions and pathways related to DEGs. As is shown in Figure 2(a) and Table 1, GO annotation of DEGs consists three parts (BP, CC, MF) and top 8 terms of each category were listed. GO analysis results revealed that the biological functions of DEGs were mainly related to immune cells infiltration and inflammatory signaling pathway. The top 10 pathways of KEGG analysis (Figure 2(b) and Table 2) indicated that DEGs were involved in signal transduction and inflammatory response-related pathways, such as chemokine signaling pathway, cytokine-cytokine receptor interaction, Th cells differentiation, IL−17 signaling pathway, and Toll-like receptor signaling pathway.

[figures omitted; refer to PDF]

Table 1

GO analysis results of DEGs (top 8 terms of each category were listed).

Ontology	ID	Description	Adj. $p$ value	Count
BP	GO:0050900	Leukocyte migration	2.33E-08	18
BP	GO:0050851	Antigen receptor-mediated signaling pathway	2.33E-08	15
BP	GO:0002449	Lymphocyte mediated immunity	6.91E-08	15
BP	GO:0051249	Regulation of lymphocyte activation	6.91E-08	17
BP	GO:0050853	B cell receptor signaling pathway	2.54E-07	10
BP	GO:0002429	Immune response-activating cell surface receptor signaling pathway	3.05E-07	16
BP	GO:0030098	Lymphocyte differentiation	4.36E-07	14
BP	GO:0006959	Humoral immune response	4.36E-07	14
CC	GO:0009897	External side of plasma membrane	1.28E-08	17
CC	GO:0042571	Immunoglobulin complex, circulating	7.02E-06	7
CC	GO:0019814	Immunoglobulin complex	8.86E-06	9
CC	GO:0042101	T cell receptor complex	0.008265	5
CC	GO:0001772	Immunological synapse	0.012398	3
CC	GO:0072562	Blood microparticle	0.012678	5
CC	GO:0042613	MHC class II protein complex	0.030215	2
CC	GO:0098802	Plasma membrane receptor complex	0.033455	6
MF	GO:0034987	Immunoglobulin receptor binding	8.86E-06	7
MF	GO:0008009	Chemokine activity	1.34E-05	6
MF	GO:0045236	CXCR chemokine receptor binding	1.51E-05	4
MF	GO:0042379	Chemokine receptor binding	5.51E-05	6
MF	GO:0003823	Antigen binding	7.28E-05	8
MF	GO:0005125	Cytokine activity	7.80E-05	9
MF	GO:0048018	Receptor ligand activity	0.000201	12
MF	GO:0001664	G protein-coupled receptor binding	0.002229	8

Table 2

Top 10 pathways of KEGG analysis.

ID	Description	Adj. $p$ value	Count
hsa04060	Cytokine-cytokine receptor interaction	1.44E-05	13
hsa04061	Viral protein interaction with cytokine and cytokine receptor	0.000459	7
hsa04062	Chemokine signaling pathway	0.002564	8
hsa04659	Th17 cell differentiation	0.003911	6
hsa03320	PPAR signaling pathway	0.005839	5
hsa04658	Th1 and Th2 cell differentiation	0.011245	5
hsa04657	IL-17 signaling pathway	0.011245	5
hsa04640	Hematopoietic cell lineage	0.01248	5
hsa04620	Toll-like receptor signaling pathway	0.013885	5
hsa05020	Prion diseases	0.026379	3

3.3. GSEA Analysis

The overall expression data and hallmark gene sets were loaded into GSEA software for further analysis. Hallmark gene sets contain 50 gene sets, representing well-annotated biological functions or processes determined by integrating multiple MSigDB datasets. According to the filter threshold of analysis results, complement and interferon alpha response gene sets were significantly upregulated in RA samples. Normalized Enrichment Score (NES) and nominal $p$ value were presented in the upper right corner of the plot (Figures 3(a) and 3(b)).

[figures omitted; refer to PDF]

3.4. Immune Infiltration Analysis

CIBERSORT analytical tool can accurately calculate levels of 22 types of leukocyte subpopulations in synovial tissues profiled by microarray. Based on significance threshold $p$ value <0.05, 4 samples (GSM1332203, GSM1332205, GSM1332206, GSM1332208) were abandoned. Figure 4(a) shows the composition of immune cells in 23 RA samples and 16 normal samples. Wilcoxon test was conducted to determine whether there was a significant difference in immune cell infiltration between RA and HC samples. In violin plot (Figure 4(b)), it was revealed that memory B cells, plasma cells, CD8 T cells, activated CD4 memory T cells, T cells follicular helper, monocytes, and macrophages M1 were abundant in RA synovial membrane.

[figures omitted; refer to PDF]

3.5. PPI Network Construction and MCODE Plug-in Analysis

PPI network constructed by STRING database was adjusted and visualized by Cytoscape (Figure 5(a)). Upregulated genes were marked with red color, and downregulated genes were blue. The diameters of nodes were positively correlated with their connectivity degree. In order to figure out the core modules of complex network, we performed MCODE plug-in analysis and identified 5 modules. The most significant module with the highest score (module score:6.923) was shown in Figure 5(b), containing 14 genes and 45 edges.

[figures omitted; refer to PDF]

3.6. Validation of Hub Genes with GEO Database

The 14 genes screened by MCODE plug-in, which could be used to distinguish RA patients from healthy controls, were considered as candidate hub genes. To further validate the expression of these 14 genes in synovium of other patients, we selected GSE77298 as testing dataset and performed ROC analysis in RStudio. The analysis results were available in Figure 6. Among the 14 genes, eight genes (CCR5, CCL5, CXCL9, CXCL10, CXCL13, PNOC, TLR8, and CD52) with AUC more than 0.80 were considered as hub genes, indicating that they have the capability to diagnose RA patients with excellent specificity and sensitivity.

[figure omitted; refer to PDF]

4. Discussion

Rheumatoid arthritis is an autoimmune nature joint disease with irreversible cartilage destruction and bone erosion [1]. If not treated promptly and effectively, RA can seriously reduce the quality of life and even cause disability. There is no doubt that understanding diseases at the molecular level will help to improve their diagnosis and treatment [26, 27]. Up to now, various biomarkers have been identified to be associated with rheumatoid arthritis and may be selected as therapeutic targets, but the detailed mechanism of gene regulation leading to disease progression remains elusive [28, 29].

In our study, we aimed to identify biomarkers of RA and uncover their biological functions through bioinformatics analysis. Dataset GSE55235 and GSE55457 were selected as training dataset in our analysis. As a result, 72 upregulated and 34 downregulated DEGs, at least 4-fold change between RA and normal samples, were screened out. Next, the DEGs were annotated by performing functional enrichment analysis, and we observed that these genes were closely related to immune response and inflammatory signaling, such as humoral immune response, chemokine signaling pathway, antigen receptor-mediated signaling pathway, Th17 cell differentiation, and IL-17 signaling pathway. Previous studies have shown that the infiltration of Th17 cells in synovium is a typical pathological change of RA [30, 31], which releases abundant IL-17A, IL-17F, IL-6, and mediates the activation of neutrophils. An animal experiment indicated that blockade of IL-17 can delay the destruction of articular cartilage by inhibiting local inflammatory reaction in the collagen-induced arthritis (CIA) mice [32]. In addition to cytokines, chemokines were also suggested to be involved in systemic inflammatory disorders [33–35]. It has been reported that elevated plasma chemokine level was discovered in RA patients, and CXCL10 can be served as a diagnostic biomarker for active rheumatoid arthritis [36]. Similarly, the GSEA results based on all gene expression information revealed that 2 gene sets (interferon alpha response and complement) were significantly enriched in RA phenotype at nominal $p$ value <0.01. Considering the crucial role of inflammatory response and leukocyte infiltration in the pathogenesis of RA, we carried out CIBERSORT analysis to explore the differences between two groups in the distribution of immune cells. Compared with healthy controls, memory B cells, plasma cells, activated CD4 memory T cells, T cells follicular helper, monocytes, and M1 macrophage were abundant in RA synovial membrane, which was consistent with the published studies [37, 38].

A PPI network of DEGs was established using STRING website and Cytoscape software. With the help of MCODE plug-in, we screened out the most significant subnetwork, which consists of 14 nodes and 45 edges. Moreover, we validated the above 14 genes by performing ROC analysis with the testing dataset GSE77298. As a result, 8 genes with $p$ value <0.05 as well as AUC >0.80 showed excellent diagnostic value for rheumatoid arthritis, and thus were considered as hub genes of RA, including CCR5, CCL5, CXCL9, CXCL10, CXCL13, PNOC, TLR8, and CD52.

The chemokine system is a large and complicated family, containing more than 50 ligands and 25 receptors [39]. Chemokines can be divided into four subfamilies, known as CXC, CC, CX3C, and XC, according to the position difference of cysteine residues. Correspondingly, there are four types of receptors that can bind to their ligands and have the ability to trigger a Gαi-mediated signaling pathway [40, 41]. In inflammatory diseases like rheumatoid arthritis, chemokines mediate the migration of leukocytes into synovial membrane, participating in angiogenesis, endothelial activation, synovial hyperplasia, and regulation of cartilage metabolism [42–45]. Pandya [36] found that chemokine CXCL9, CXCL10, CXCL13, CCL4, and CCL22 were significantly higher in the blood plasma of RA patients compared to healthy people by multivariate discriminant analysis. Besides, various studies based on targeting of chemokines and their receptors have been widely taken out for therapeutic purposes. Tofacitinib is a Janus kinase inhibitor and has been proved to inhibit chemokines secretion in synovium, including CXCL10, CXCL13, and CCL2 [46]. In animal experiments, the use of receptors inhibitors also achieved fine therapeutic effect. J­113863, an antagonist of CCR1, has positive effects on the CIA model of murine [47], and adjuvant-induced arthritis in rats can be inhibited by Met-RANTES, a drug blocking CCR1 and CCR5 [48]. All these findings suggested that chemokines act as a key factor in the pathological process of RA and deserve more attention.

TLR8, a gene encoding Toll-like receptor 8, was predominantly expressed in peripheral blood leukocytes, playing a fundamental role in antimicrobial immune responses and autoimmune inflammation [49]. TLR7 and TLR8 are located at the membranes of endosomal compartment, involved in recognizing viral RNA [50, 51]. Toll-like receptor signaling pathway, stimulated by the adaptor protein MyD88, modulates NFкB, IRF-7, and MAPK activation, resulting in the release of proinflammatory cytokines and cell adhesion molecules [52]. CD52, also known as EDDM5, is a low molecular weight glycoprotein that is abundant in B and T cells [53]. Alemtuzumab, the first therapeutic anti-CD52 antibody, has been shown to be effective in the treatment of autoimmune diseases such as RA and inflammatory bowel disease. However, in recent studies, alemtuzumab has been found to result in long-term immunosuppression, particularly depleting CD4+ T cells, and increase the risk of opportunistic infection [54, 55]. Therefore, in the future clinical practice, it is essential to pay attention to the dynamics of immune reconstitution and the results of immunosuppression. Unlike the genes discussed above, PNOC is rarely mentioned in previous rheumatoid arthritis studies, and thus needs more exploration. This gene encodes a preproprotein, which is the precursor of nociceptin. A previous study demonstrated that the mRNA expression of PNOC and NOP was suppressed by lipopolysaccharide as well as inflammatory cytokines [56]. In the current study, PNOC gene is upregulated in RA patients, which may stimulate the secretion of cytokines, further aggravate the inflammatory response of synovium and become a vicious circle.

5. Conclusions

In summary, we integrated multiple bioinformatics tools and found that chemokines and immune cell infiltration were extremely critical factors in the progression of rheumatoid arthritis. The eight hub genes we identified may serve as potential therapeutic targets for RA and further investigations are required to support our conclusions.

Acknowledgments

The present study was supported by the National Natural Science Foundation of China (grant number 81904053); Plan Project of Hangzhou Health Science and Technology Department (grant number 2018B028); and Opening Project of Zhejiang Provincial Preponderant and Characteristic Subject of Key University, Zhejiang Chinese Medical University (grant number ZYX2018008).

References

[1] J. S. Smolen, D. Aletaha, A. Barton, G. R. Burmester, P. Emery, G. S. Firestein, A. Kavanaugh, I. B. McInnes, D. H. Solomon, V. Strand, K. Yamamoto, "Rheumatoid arthritis," Nature Reviews Disease Primers, vol. 4 no. 1,DOI: 10.1038/nrdp.2018.1, 2018.

[2] E. Myasoedova, C. S. Crowson, H. M. Kremers, T. M. Therneau, S. E. Gabriel, "Is the incidence of rheumatoid arthritis rising?: results from Olmsted County, Minnesota, 1955-2007," Arthritis and Rheumatism, vol. 62 no. 6, pp. 1576-1582, DOI: 10.1002/art.27425, 2010.

[3] G. J. Tobon, P. Youinou, A. Saraux, "The environment, geo-epidemiology, and autoimmune disease: rheumatoid arthritis," Journal of Autoimmunity, vol. 35 no. 1, pp. 10-14, DOI: 10.1016/j.jaut.2009.12.009, 2010.

[4] C. M. Coffey, C. S. Crowson, E. Myasoedova, E. L. Matteson, J. M. Davis, "Evidence of diagnostic and treatment delay in Seronegative rheumatoid arthritis: missing the window of opportunity," Mayo Clinic Proceedings, vol. 94 no. 11, pp. 2241-2248, DOI: 10.1016/j.mayocp.2019.05.023, 2019.

[5] P. F. Whiting, N. Smidt, J. A. Sterne, R. Harbord, A. Burton, M. Burke, R. Beynon, Y. Ben-Shlomo, J. Axford, P. Dieppe, "Systematic review: accuracy of anti-citrullinated peptide antibodies for diagnosing rheumatoid arthritis," Annals of Internal Medicine, vol. 152 no. 7, pp. 456-464, DOI: 10.7326/0003-4819-152-7-201004060-00010, 2010.

[6] A. Kleyer, G. Schett, "Arthritis and bone loss: a hen and egg story," Current Opinion in Rheumatology, vol. 26 no. 1, pp. 80-84, DOI: 10.1097/BOR.0000000000000007, 2014.

[7] P. Kopesky, K. Tiedemann, D. Alkekhia, C. Zechner, B. Millard, B. Schoeberl, S. V. Komarova, "Autocrine signaling is a key regulatory element during osteoclastogenesis," Biology Open, vol. 3 no. 8, pp. 767-776, DOI: 10.1242/bio.20148128, 2014.

[8] M. C. Kraan, D. D. Patel, J. J. Haringman, M. D. Smith, H. Weedon, M. J. Ahern, F. C. Breedveld, P. P. Tak, "The development of clinical signs of rheumatoid synovial inflammation is associated with increased synthesis of the chemokine CXCL8 (interleukin-8)," Arthritis Research, vol. 3 no. 1, pp. 65-71, DOI: 10.1186/ar141, 2001.

[9] A. Kleyer, S. Finzel, J. Rech, B. Manger, M. Krieter, F. Faustini, E. Araujo, A. J. Hueber, U. Harre, K. Engelke, G. Schett, "Bone loss before the clinical onset of rheumatoid arthritis in subjects with anticitrullinated protein antibodies," Annals of the Rheumatic Diseases, vol. 73 no. 5, pp. 854-860, DOI: 10.1136/annrheumdis-2012-202958, 2014.

[10] A. Krishnamurthy, V. Joshua, A. Haj Hensvold, T. Jin, M. Sun, N. Vivar, A. J. Ytterberg, M. Engström, C. Fernandes-Cerqueira, K. Amara, M. Magnusson, G. Wigerblad, J. Kato, J. M. Jiménez-Andrade, K. Tyson, S. Rapecki, K. Lundberg, S. B. Catrina, P. J. Jakobsson, C. Svensson, V. Malmström, L. Klareskog, H. Wähämaa, A. I. Catrina, "Identification of a novel chemokine-dependent molecular mechanism underlying rheumatoid arthritis-associated autoantibody-mediated bone loss," Annals of the Rheumatic Diseases, vol. 75 no. 4, pp. 721-729, DOI: 10.1136/annrheumdis-2015-208093, 2016.

[11] J. H. Malone, B. Oliver, "Microarrays, deep sequencing and the true measure of the transcriptome," BMC Biology, vol. 9 no. 1,DOI: 10.1186/1741-7007-9-34, 2011.

[12] A. A. Alizadeh, M. B. Eisen, R. E. Davis, C. Ma, I. S. Lossos, A. Rosenwald, J. C. Boldrick, H. Sabet, T. Tran, X. Yu, J. I. Powell, L. Yang, G. E. Marti, T. Moore, J. Hudson, L. Lu, D. B. Lewis, R. Tibshirani, G. Sherlock, W. C. Chan, T. C. Greiner, D. D. Weisenburger, J. O. Armitage, R. Warnke, R. Levy, W. Wilson, M. R. Grever, J. C. Byrd, D. Botstein, P. O. Brown, L. M. Staudt, "Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling," Nature, vol. 403 no. 6769, pp. 503-511, DOI: 10.1038/35000501, 2000.

[13] T. R. Golub, D. K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. P. Mesirov, H. Coller, M. L. Loh, J. R. Downing, M. A. Caligiuri, C. D. Bloomfield, E. S. Lander, "Molecular classification of cancer: class discovery and class prediction by gene expression monitoring," Science, vol. 286 no. 5439, pp. 531-537, DOI: 10.1126/science.286.5439.531, 1999.

[14] D. Woetzel, R. Huber, P. Kupfer, D. Pohlers, M. Pfaff, D. Driesch, T. Häupl, D. Koczan, P. Stiehl, R. Guthke, R. W. Kinne, "Identification of rheumatoid arthritis and osteoarthritis patients by transcriptome-based rule set generation," Arthritis Research & Therapy, vol. 16 no. 2,DOI: 10.1186/ar4526, 2014.

[15] M. G. A. Broeren, M. de Vries, M. B. Bennink, O. J. Arntz, A. B. Blom, M. I. Koenders, P. L. E. M. van Lent, P. M. van der Kraan, W. B. van den Berg, F. A. J. van de Loo, "Disease-regulated gene therapy with anti-inflammatory interleukin-10 under the control of the CXCL10 promoter for the treatment of rheumatoid arthritis," Human Gene Therapy, vol. 27 no. 3, pp. 244-254, DOI: 10.1089/hum.2015.127, 2016.

[16] M. E. Ritchie, B. Phipson, Y. H. Di Wu, C. W. Law, W. Shi, G. K. Smyth, "limma powers differential expression analyses for RNA-sequencing and microarray studies," Nucleic Acids Research, vol. 43 no. 7,DOI: 10.1093/nar/gkv007, 2015.

[17] G. Yu, L. G. Wang, Y. Han, Q. Y. He, "clusterProfiler: an R package for comparing biological themes among gene clusters," OMICS, vol. 16 no. 5, pp. 284-287, DOI: 10.1089/omi.2011.0118, 2012.

[18] A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, S. L. Pomeroy, T. R. Golub, E. S. Lander, J. P. Mesirov, "Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles," Proceedings of the National Academy of Sciences, vol. 102 no. 43, pp. 15545-15550, DOI: 10.1073/pnas.0506580102, 2005.

[19] V. K. Mootha, C. M. Lindgren, K. F. Eriksson, A. Subramanian, S. Sihag, J. Lehar, P. Puigserver, E. Carlsson, M. Ridderstråle, E. Laurila, N. Houstis, M. J. Daly, N. Patterson, J. P. Mesirov, T. R. Golub, P. Tamayo, B. Spiegelman, E. S. Lander, J. N. Hirschhorn, D. Altshuler, L. C. Groop, "PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes," Nature Genetics, vol. 34 no. 3, pp. 267-273, DOI: 10.1038/ng1180, 2003.

[20] A. Liberzon, A. Subramanian, R. Pinchback, H. Thorvaldsdottir, P. Tamayo, J. P. Mesirov, "Molecular signatures database (MSigDB) 3.0," Bioinformatics, vol. 27 no. 12, pp. 1739-1740, DOI: 10.1093/bioinformatics/btr260, 2011.

[21] A. Liberzon, C. Birger, H. Thorvaldsdóttir, M. Ghandi, J. P. Mesirov, P. Tamayo, "The molecular signatures database (MSigDB) hallmark gene set collection," Cell Systems, vol. 1 no. 6, pp. 417-425, DOI: 10.1016/j.cels.2015.12.004, 2015.

[22] A. M. Newman, C. L. Liu, M. R. Green, A. J. Gentles, W. Feng, Y. Xu, C. D. Hoang, M. Diehn, A. A. Alizadeh, "Robust enumeration of cell subsets from tissue expression profiles," Nature Methods, vol. 12 no. 5, pp. 453-457, DOI: 10.1038/nmeth.3337, 2015.

[23] D. Szklarczyk, A. L. Gable, D. Lyon, A. Junge, S. Wyder, J. Huerta-Cepas, M. Simonovic, N. T. Doncheva, J. H. Morris, P. Bork, L. J. Jensen, C. . Mering, "STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets," Nucleic Acids Research, vol. 47 no. D1, pp. D607-d613, DOI: 10.1093/nar/gky1131, 2019.

[24] G. D. Bader, C. W. Hogue, "An automated method for finding molecular complexes in large protein interaction networks," BMC Bioinformatics, vol. 4 no. 1,DOI: 10.1186/1471-2105-4-2, 2003.

[25] X. Robin, N. Turck, A. Hainard, N. Tiberti, F. Lisacek, J.-C. Sanchez, M. Müller, "pROC: an open-source package for R and S+ to analyze and compare ROC curves," BMC Bioinformatics, vol. 12 no. 1,DOI: 10.1186/1471-2105-12-77, 2011.

[26] Y. Liao, Y. Wang, M. Cheng, C. Huang, X. Fan, "Weighted gene coexpression network analysis of features that control cancer stem cells reveals prognostic biomarkers in lung adenocarcinoma," Frontiers in Genetics, vol. 11,DOI: 10.3389/fgene.2020.00311, 2020.

[27] Y. Liao, H. Xiao, M. Cheng, X. Fan, "Bioinformatics analysis reveals biomarkers with cancer stem cell characteristics in lung squamous cell carcinoma," Frontiers in Genetics, vol. 11,DOI: 10.3389/fgene.2020.00427, 2020.

[28] W. P. Maksymowych, M. A. Brown, "Genetics of ankylosing spondylitis and rheumatoid arthritis: where are we at currently, and how do they compare?," Clinical and Experimental Rheumatology, vol. 27, pp. S20-S25, 2009.

[29] Y. Huang, S. Zheng, R. Wang, C. Tang, J. Zhu, J. Li, "CCL5 and related genes might be the potential diagnostic biomarkers for the therapeutic strategies of rheumatoid arthritis," Clinical Rheumatology, vol. 38 no. 9, pp. 2629-2635, DOI: 10.1007/s10067-019-04533-1, 2019.

[30] A. P. Cope, H. Schulze-Koops, M. Aringer, "The central role of T cells in rheumatoid arthritis," Clinical and Experimental Rheumatology, vol. 25, pp. S4-11, 2007.

[31] T. Dhaouadi, M. Chahbi, Y. Haouami, I. Sfar, L. Abdelmoula, T. Ben Abdallah, Y. Gorgi, "IL-17A, IL-17RC polymorphisms and IL17 plasma levels in Tunisian patients with rheumatoid arthritis," PLoS One, vol. 13 no. 3, article e0194883,DOI: 10.1371/journal.pone.0194883, 2018.

[32] E. Lubberts, M. I. Koenders, B. Oppers-Walgreen, L. van den Bersselaar, C. J. J. Coenen-de Roo, L. A. B. Joosten, W. B. van den Berg, "Treatment with a neutralizing anti-murine interleukin-17 antibody after the onset of collagen-induced arthritis reduces joint inflammation, cartilage destruction, and bone erosion," Arthritis and Rheumatism, vol. 50 no. 2, pp. 650-659, DOI: 10.1002/art.20001, 2004.

[33] I. B. McInnes, G. Schett, "The pathogenesis of rheumatoid arthritis," New England Journal of Medicine, vol. 365 no. 23, pp. 2205-2219, DOI: 10.1056/NEJMra1004965, 2011.

[34] M. Mellado, L. Martínez-Muñoz, G. Cascio, P. Lucas, J. L. Pablos, J. M. Rodríguez-Frade, "T cell migration in rheumatoid arthritis," Frontiers in Immunology, vol. 6, 2015.

[35] Z. Szekanecz, A. E. Koch, "Successes and failures of chemokine-pathway targeting in rheumatoid arthritis," Nature Reviews Rheumatology, vol. 12 no. 1,DOI: 10.1038/nrrheum.2015.157, 2016.

[36] J. M. Pandya, A. C. Lundell, K. Andersson, I. Nordström, E. Theander, A. Rudin, "Blood chemokine profile in untreated early rheumatoid arthritis: CXCL10 as a disease activity marker," Arthritis Research & Therapy, vol. 19 no. 1,DOI: 10.1186/s13075-017-1224-1, 2017.

[37] F. Ginhoux, S. Jung, "Monocytes and macrophages: developmental pathways and tissue homeostasis," Nature Reviews. Immunology, vol. 14 no. 6, pp. 392-404, DOI: 10.1038/nri3671, 2014.

[38] P. Italiani, D. Boraschi, "From monocytes to M1/M2 macrophages: Phenotypical vs. functional differentiation," Frontiers in Immunology, vol. 5,DOI: 10.3389/fimmu.2014.00514, 2014.

[39] E. M. Borroni, B. Savino, R. Bonecchi, M. Locati, "Chemokines sound the alarmin: the role of atypical chemokine in inflammation and cancer," Seminars in Immunology, vol. 38, pp. 63-71, DOI: 10.1016/j.smim.2018.10.005, 2018.

[40] F. Bachelerie, A. Ben-Baruch, A. M. Burkhardt, C. Combadiere, J. M. Farber, G. J. Graham, R. Horuk, A. H. Sparre-Ulrich, M. Locati, A. D. Luster, A. Mantovani, K. Matsushima, P. M. Murphy, R. Nibbs, H. Nomiyama, C. A. Power, A. E. I. Proudfoot, M. M. Rosenkilde, A. Rot, S. Sozzani, M. Thelen, O. Yoshie, A. Zlotnik, "International Union of Basic and Clinical Pharmacology. [corrected]. LXXXIX. Update on the extended family of chemokine receptors and introducing a new nomenclature for atypical chemokine receptors," Pharmacological Reviews, vol. 66 no. 1,DOI: 10.1124/pr.113.007724, 2013.

[41] M. Thelen, "Dancing to the tune of chemokines," Nature Immunology, vol. 2 no. 2, pp. 129-134, DOI: 10.1038/84224, 2001.

[42] C. X. Maracle, S. W. Tas, "Inhibitors of angiogenesis: ready for prime time?," Best Practice & Research. Clinical Rheumatology, vol. 28 no. 4, pp. 637-649, DOI: 10.1016/j.berh.2014.10.012, 2014.

[43] Z. Szekanecz, A. E. Koch, "Chemokines and angiogenesis," Current Opinion in Rheumatology, vol. 13 no. 3, pp. 202-208, DOI: 10.1097/00002281-200105000-00009, 2001.

[44] Z. Szekanecz, A. E. Koch, P. P. Tak, "Chemokine and chemokine receptor blockade in arthritis, a prototype of immune-mediated inflammatory diseases," The Netherlands Journal of Medicine, vol. 69 no. 9, pp. 356-366, 2011.

[45] Z. Szekanecz, A. Vegvari, Z. Szabo, A. E. Koch, "Chemokines and chemokine receptors in arthritis," Frontiers in Bioscience (Scholar Edition), vol. 2, pp. 153-167, DOI: 10.2741/s53, 2010.

[46] D. L. Boyle, K. Soma, J. Hodge, A. Kavanaugh, D. Mandel, P. Mease, R. Shurmur, A. K. Singhal, N. Wei, S. Rosengren, I. Kaplan, S. Krishnaswami, Z. Luo, J. Bradley, G. S. Firestein, "The JAK inhibitor tofacitinib suppresses synovial JAK1-STAT signalling in rheumatoid arthritis," Annals of the Rheumatic Diseases, vol. 74 no. 6, pp. 1311-1316, DOI: 10.1136/annrheumdis-2014-206028, 2015.

[47] M. Amat, C. F. Benjamim, L. M. Williams, N. Prats, E. Terricabras, J. Beleta, S. L. Kunkel, N. Godessart, "Pharmacological blockade of CCR1 ameliorates murine arthritis and alters cytokine networks in vivo," British Journal of Pharmacology, vol. 149 no. 6, pp. 666-675, DOI: 10.1038/sj.bjp.0706912, 2006.

[48] S. Shahrara, A. E. I. Proudfoot, J. M. Woods, J. H. Ruth, M. A. Amin, C. C. Park, C. S. Haas, R. M. Pope, G. K. Haines, Y. Y. Zha, A. E. Koch, "Amelioration of rat adjuvant-induced arthritis by Met-RANTES," Arthritis and Rheumatism, vol. 52 no. 6, pp. 1907-1919, DOI: 10.1002/art.21033, 2005.

[49] N. Ishii, K. Funami, M. Tatematsu, T. Seya, M. Matsumoto, "Endosomal localization of TLR8 confers distinctive proteolytic processing on human myeloid cells," Journal of Immunology, vol. 193 no. 10, pp. 5118-5128, DOI: 10.4049/jimmunol.1401375, 2014.

[50] F. Heil, H. Hemmi, H. Hochrein, F. Ampenberger, C. Kirschning, S. Akira, G. Lipford, H. Wagner, S. Bauer, "Species-specific recognition of single-stranded RNA via toll-like receptor 7 and 8," Science, vol. 303 no. 5663, pp. 1526-1529, DOI: 10.1126/science.1093620, 2004.

[51] A. L. Blasius, B. Beutler, "Intracellular toll-like receptors," Immunity, vol. 32 no. 3, pp. 305-315, DOI: 10.1016/j.immuni.2010.03.012, 2010.

[52] U. Ohto, H. Tanji, T. Shimizu, "Structure and function of toll-like receptor 8," Microbes and Infection, vol. 16 no. 4, pp. 273-282, DOI: 10.1016/j.micinf.2014.01.007, 2014.

[53] G. Hale, "The CD52 antigen and development of the CAMPATH antibodies," Cytotherapy, vol. 3 no. 3, pp. 137-143, DOI: 10.1080/146532401753174098, 2001.

[54] L. Alinari, R. Lapalombella, L. Andritsos, R. A. Baiocchi, T. S. Lin, J. C. Byrd, "Alemtuzumab (Campath-1H) in the treatment of chronic lymphocytic leukemia," Oncogene, vol. 26 no. 25, pp. 3644-3653, DOI: 10.1038/sj.onc.1210380, 2007.

[55] C. E. Dearden, E. Matutes, "Alemtuzumab in T-cell lymphoproliferative disorders," Best Practice & Research. Clinical Haematology, vol. 19 no. 4, pp. 795-810, DOI: 10.1016/j.beha.2006.05.005, 2006.

[56] L. Zhang, F. Stuber, U. M. Stamer, "Inflammatory mediators influence the expression of nociceptin and its receptor in human whole blood cultures," PLoS One, vol. 8 no. 9, article e74138,DOI: 10.1371/journal.pone.0074138, 2013.