Identification of key genes associated with rheumatoid arthritis with bioinformatics approach.

We aimed to identify key genes associated with rheumatoid arthritis (RA).The microarray datasets of GSE1919, GSE12021, and GSE21959 (35 RA samples and 32 normal controls) were downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) in RA samples were identified using the t test in limma package. Functional enrichment analysis was performed using clusterProfiler package. A protein-protein interaction (PPI) network of selected DEGs was constructed based on the Human Protein Reference Database. Active modules were explored using the jActiveModules plug-in in the Cytoscape Network Modeling package.In total, 537 DEGs in RA samples were identified, including 241 upregulated and 296 downregulated genes. A total of 24,451 PPI pairs were collected, and 5 active modules were screened. Furthermore, 19 submodules were acquired from the 5 active modules. Discs large homolog 1 (DLG1) and related DEGs such as guanylate cyclase 1, soluble, alpha 2 (GUCY1A2), N-methyl d-aspartate receptor 2A subunit (GRIN2A), and potassium voltage-gated channel member 1 (KCNA1) were identified in 8 submodules. Plasminogen (PLG) and related DEGs such as chemokine (C-X-C motif) ligand 2 (CXCL2), laminin, alpha 3 (LAMA3), complement component 7 (C7), and coagulation factor X (F10) were identified in 4 submodules.Our results indicate that DLG1, GUCY1A2, GRIN2A, KCNA1, PLG, CXCL2, LAMA3, C7, and F10 may play key roles in the progression of RA and may serve as putative therapeutic targets for treating RA.

Introduction

Rheumatoid arthritis (RA) is a chronic autoimmune disorder that mainly affects the synovial membrane, cartilage, and bone.[ It affects 1% of the population and is related to significant morbidity and increased mortality.[ Because RA is incurable, it causes high economic burdens that severely reduce the quality of life.[ The cause of RA is still unknown.[ Reliable predictive biomarkers for RA prognosis and therapeutic response are limited.[ Therefore, elucidation of molecular mechanisms underlying RA development would help in conceiving more effective therapeutic strategies than before for treating the disease.

RA is characterized by chronic inflammation, synovial hyperplasia, cartilage and bone destruction, and formation of pannus.[ It involves a complex interplay among genotype, environmental triggers, and chance.[ Numerous studies have confirmed that many activated cell types such as T and B cells, monocytes/macrophages, endothelial cells, and synovial fibroblasts play a role in the development and progression of RA.[ Moreover, genome-wide analyses have verified that immune regulatory factors underlie the disease.[ Cytokines such as interleukin (IL)-1α, IL-8, and tumor necrosis factors (TNFs) have been directly implicated in a broad range of immune and inflammatory processes, thereby being responsible for joint damage and RA progression.[ Extended analyses of gene expression profiles have shown that several gene groups such as activator protein-1, nuclear factor kappa B, and Smad family members show binding activity at their cognate recognition sites in the promoters of cytokines, and are, therefore, key players in RA development and progression.[ Consequently, blocking TNFs and IL-1α has proved successful in treating RA.[ In addition, increased levels of sphingosine kinase 1 in the synovium membrane constitute RA biomarkers.[ However, identifying newer molecular targets remains necessary despite such advances in treating RA to foster the development of new therapeutics with improved outcomes and accuracy of diagnosis.

Previously, datasets GSE1919 and GSE12021 were used to decode differential biomarkers of RA and osteoarthritis for disease-specific therapeutics using statistical analysis.[ These datasets were also used to identify a panel of potential molecular targets for diagnosis and treatment of RA.[ The microarray dataset GSE21959 was used to show that gene expression is considerably altered in RA synovial fibroblasts compared with healthy synovial fibroblasts after hypoxia induction and to identify new factors and pathways relevant to the pathogenesis of chronic arthritis.[ In comparison with these, we integrated the above 3 gene expression datasets in the present study and used comprehensive bioinformatics methods to screen differentially expressed genes (DEGs) associated with RA. Functional enrichment analysis, protein–protein interaction (PPI) network construction, and active module mining were then performed to identify key genes involved in the progression of RA. Our findings will help in the discovery of potential targets for therapeutic intervention in patients with RA.

Methods

Data resources

Data on the expression profile of GSE1919,[ GSE12021,[ and GSE21959[ related to RA were acquired from the Gene Expression Omnibus repository by the National Center of Biotechnology Information based on the platform of GPL91 Affymetrix Human Genome U95A Array, GPL96 Affymetrix Human Genome U133A/B Array, and GPL4133 Agilent-014850 Whole Human Genome Microarray 4 × 44K G4112F (Feature Number version), respectively. GSE1919 datasets contained synovial membrane samples derived from 5 patients with RA and 5 normal controls; GSE12021 contained 9 RA samples and 12 normal controls; and GSE21959 contained 18 RA samples and 18 normal controls. Therefore, 67 datasets of synovial fibroblasts in total (35 RA samples and 32 normal controls) were used for subsequently analyses. Ethical approval was not necessary in our study because we downloaded the expression profiles from the public database and do not perform any experiments in patients or animals.

Data preprocessing

We downloaded the series matrix files of the 3 datasets, and the probe name of each series matrix was transformed into a gene symbol based on the Affy probe annotation files. If multiple probes corresponded to the same gene symbol, the mean value was calculated using the aggregate function[ in R as the expression value of that particular gene. If the expression value of the probe was absent, the nearest neighbor averaging (KNN)[ of the impute package[ in R was used to supplement. Finally, quartile normalization was performed using the preprocessCore[ in R, and the gene expression matrix of each sample was acquired.

All expression estimates were log2 transformed and merged using cross-platform normalization, which was performed using the CONOR[ package in R. If different studies comprised similar or common gene symbols, 2 expression data of the same gene symbols were normalized for producing a new dataset. Then, the newly produced data were renormalized with the next data.

Identification of DEGs

DEGs in RA samples compared with normal controls were screened using the t test[ method in the limma[ package in R. Then, the P value was adjusted using the Benjamini and Hochberg method.[ The | log2 fold change (FC) | >0.585 and an adjusted P value of <0.05 were set as the threshold.

Functional enrichment analysis of DEGs

Gene ontology (GO) analysis is increasingly applied for functional studies of large-scale genomic or transcriptomic data,[ which comprises 3 independent ontologies, including biological process (BP), molecular function (MF), and cellular component (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is the major public database containing the information of biochemistry pathways.[

We used the clusterProfiler[ package in R to annotate and visualize GO terms (associated with BP, MF, and CC) and KEGG pathway of upregulated and downregulated genes, respectively.

Construction of the PPnetwork and screening of active modules

Human Protein Reference Database (HPRD)[ is a database containing curated proteomic information pertaining to human proteins. The human PPI data on HPRD (Release 9) consists of 39,240 interactions among 9617 genes. We used HPRD (Release 9) to identify PPI pairs of DEGs and construct the PPI network.

Active module hypothesis was proposed by Ideker et al[ for clarifying underlying mechanisms governing observed changes in gene expression, which integrated expression profiling, large-scale proteomics, and PPI network construction along with using a rigorous statistical measure for identifying active modules with a search algorithm. Therefore, active modules were screened within the PPI network considering the active module hypothesis.

We also used the jActiveModules plug-in[ available for the Cytoscape Network Modeling package[ to identify active modules. At each step of the search, all adjacent proteins are considered for inclusion in the resultant network. Active modules were identified with a local search of “depth” = 2 and “max depth” = 2. Then, the connectivity degree of each node in modules was estimated to identify the hub nodes.

If many nodes were present in active modules, smaller submodules needed to be acquired. To produce smaller submodules, we repeated the active modules search within each original submodule with a local search of “depth” = 1 and “max depth” = 1 to identify singleton nodes with a significant number of neighbors.

Results

Data preprocessing and identification of DEGs

Based on the merged information, 57 samples corresponding to 8453 gene expression levels were acquired. After preprocessing, data were successfully normalized and could be used for further study.

Using the limma package with | log2 FC |>0.585 and an adjusted P value of <0.05 set as a threshold, we ultimately obtained 537 DEGs in RA samples compared with normal samples, including 241 upregulated and 296 downregulated genes.

GO terms and KEGG pathways were considerably enriched by upregulated and downregulated DEGs, respectively. Upregulated genes were significantly associated with immune response, response to stimulus, and immune system process (Table 1A). Downregulated genes were significantly associated with multicellular organism process, single-organism process, single-multicellular organism process (Table 1B). In addition, upregulated genes were significantly enriched in cell adhesion molecules, primary immunodeficiency, hematopoietic cell lineage, cytokine–cytokine receptor interaction, and Staphylococcus aureus infections (Table 2A). Downregulated genes were mainly enriched in tyrosine metabolism, pathways in cancer, insulin signaling pathway, cytokine–cytokine receptor interaction, and amoebiasis cancer (Table 2B).

Table 1

The enriched GO terms of differentially expressed genes.

Table 2

The top 5 enriched KEGG pathways of differentially expressed genes.

PPnetwork construction and active module screening

Based on the information on HPRD, 24,451 PPI pairs corresponding to 8453 genes were acquired by mapping into HPRD-. Then, 5 active modules and 19 submodules were screened (Fig. 1). In active module (A), the hub genes with a higher degree were androgen receptor (AR), v-src avian sarcoma (Schmidt-Ruppin A-2) viral oncogene homolog (SRC), and CREB-binding protein (CREBBP); in active module (B), they were SRC, caveolin 1, caveolae protein, 22kDa (CAV1), and amyloid beta (A4) precursor protein (APP); in active module (C), they were c-abl oncogene 1, nonreceptor tyrosine kinase (ABL1), CREBBP, and SRC; and in active module (D), they were AR, SRC, and epidermal growth factor receptor (EGFR). The hub genes in active module (E) were protein kinase C, alpha (PRKCA), SRC, and EGFR. Our results showed that most of DEGs present in active modules were downregulated genes rather than upregulated genes.

Figure 1

Five active modules and 19 submodules. The yellow nodes indicate normally expressed genes, red nodes indicate upregulated genes, and green nodes indicate downregulated genes. Color shades of nodes are proportional to the | log2 fold change (FC) |, and node size is inversely proportional to and adjusted P value. The edges between nodes indicate interaction between these genes.

The connectivity degrees of many hub nodes of submodules were not <3, and they were not DEGs. These hub nodes of submodules were as follows: plasminogen (PLG), phosphatase and tensin homolog (PTEN), discs large homolog 1 (DLG1), v-erb-b2 avian erythroblastic leukemia viral oncogene homolog 4 (ERBB4), growth factor receptor-bound protein 7 (GRB7), v-erb-b2 avian erythroblastic leukemia viral oncogene homolog 2 (ERBB2), epidermal growth factor (EGF), discs large homolog 3 (DLG3), kallikrein-related peptidase 2 (KLK2), and SRC.

DLG1 and related DEGs such as guanylate cyclase 1 soluble alpha 2 (GUCY1A2), N-methyl d-aspartate receptor 2A subunit (GRIN2A), and potassium voltage-gated channel member 1 (KCNA1) were present in 8 submodules. PLG and associated DEGs such as chemokine (C-X-C motif) ligand 2 (CXCL2), laminin, alpha 3 (LAMA3), complement component 7 (C7), and coagulation factor X (F10) were present in 4 submodules.

Discussion and conclusions

We used a bioinformatics approach to identify key genes associated with RA progression. In total, 537 DEGs in RA samples were identified, including 241 upregulated and 296 downregulated genes. Five active modules and 19 submodules were also acquired. DLG1 and related DEGs such as GUCY1A2, GRIN2A, and KCNA1, were present in 8 submodules. PLG and related DEGs such as CXCL2, LAMA3, C7, and F10 were present in 4 submodules. Although DLG1 and PLG were not DEGs, they were strongly associated with DEGs involved in RA progression.

DLG1 encodes a multidomain scaffolding protein that is required for T cell receptor (TCR)-induced activation of regulatory T cell function.[ A previous study showed that T cells play an important role in the development of RA.[ Factors involved in T cell activation, such as CD28 and T cell activation RhoGTPase-activating protein, are strongly linked to RA and, therefore, highlight the role of T cells in RA.[ DLG1 is involved in the generation of memory T cells and regulation of T cell proliferation, migration, and Ag-receptor signaling.[ Moreover, DLG1 functions as a negative regulator of TCR-induced proliferative responses, thereby decreasing the proliferation of T lymphocytes.[ Thus, our results are consistent with those of previous studies and suggest that the downregulated DLG1 may enhance the function of regulatory T cells, thereby exacerbating RA.

We also found that downregulated DEGs such as GUCY1A2, GRIN2A, and KCNA directly interacted with DGL1 in several submodules. GUCY1A2 is involved in EDNRB signaling, which may be a part of RA development.[GRIN2A (also known as NMDAR2A) helps increase the expression of proinflammatory molecules such as IL-1b and TNF-α following N-methyl D-aspartate-induced excitotoxicity in the postnatal brain.[ These cytokines have been confirmed to play important roles in the pathogenesis of RA.[ Furthermore, the expression of GUCY1A2 and GRIN2A is downregulated in RA synovial fibroblasts compared with that in healthy synovial fibroblasts, specifically in hypoxia.[KCNA1 is a member of potassium voltage-gated channel proteins, which regulate the membrane potential of T cells and thereby control the calcium signaling response necessary for lymphocyte activation.[ The function of T cells is strongly associated with RA[ and calcium signaling, which are few of the signaling pathways disrupted in RA.[ Therefore, GUCY1A2, GRIN2A, and KCNA1 may also be key molecules involved in the pathogenesis of RA. Our results also suggest that DLG1 may be indirectly involved in RA progression owing to its interaction with the above 3 DEGs.

PLG was also found to be a crucial protein in several submodules. PLG induces the production of IL-1β and TNF-α in response to reactive oxygen species,[ which are also important molecules in the development of RA. Moreover, plasmin is formed upon the cleavage of PLG by specific PLG activators, and further induces the expression of inflammatory cytokines.[ Enhanced urokinase plasminogen activator activity is a key component of both the inflammatory and tissue remodeling processes occurring in the joints of patients with RA.[ In other words, PLG has been shown to have both a positive and a negative influence with respect to RA, depending on joint location, in mice.[ Therefore, downregulated PLG may be a joint-specific determinant playing a role in the progression of RA.

CXCL2, LAMA3, C7, and F10 directly interacted with PLG. CXCL2 is considered a proinflammatory and matrix-destructive factor in RA.[ CXCR2, the receptor of CXCL2, inhibits acute IL-8/CXCL8- or LPS-induced arthritis in rabbits.[LAMA3, which encodes an angiogenic protein, is induced by hypoxia in macrophages and effectively regulates angiogenesis.[ Synovial hyperplasia in RA is accompanied by active angiogenesis and expansion of blood vessels,[ highlighting the role of LAMA3 in RA progression. C7 is a key component of the complement system, which contributes to the pathogenesis of several autoimmune and inflammatory conditions, including RA.[ C7 deficiency is strongly associated with the occurrence of RA.[ F10, one of the blood coagulation factors, is associated with the activation of the coagulation cascade that widely occurs in RA, and reduced levels of coagulation factors have been found in the synovial fluids of patients with RA.[ These DEGs are involved in different mechanisms underlying the pathogenesis of RA, and we speculate that they contribute to RA progression, albeit by playing different roles. Although current evidence on the direct association between PLG and DEGs is limited, we believe that PLG may be involved in the pathogenesis of RA via interacting with these DEGs based on our present findings.

Nevertheless, there was no experimental validation, such as quantitative real-time PCR and western blot analysis to determine the expression levels of these key genes we identified. Moreover, the significant role of these candidate genes in the development of RA was not investigated. More studies with high throughput data and experiment validation are still needed to verify our observation and speculation.

To summarize, DLG1, GUCY1A2, GRIN2A, KCNA1, PLG, CXCL2, LAMA3, C7, and F10 may play a key role in the progression and development of RA. These genes may be promising as therapeutic targets for treating RA.