Identification and validation of hub genes and potential drugs involved in osteoarthritis through bioinformatics analysis

Purpose: Osteoarthritis (OA) is a common degenerative disease, which still lacks specific therapeutic drugs. Synovitis is one of the most important pathological process in OA. Therefore, we aim to identify and analyze the hub genes and their related networks of OA synovium with bioinformatics tools to provide theoretical basis for potential drugs.

Materials and methods: Two datasets were obtained from GEO. DEGs and hub genes of OA synovial tissue were screened through Gene Ontology (GO) annotation, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment as well as protein—protein interaction (PPI) network analysis. Subsequently, the correlation between expression of hub genes and ferroptosis or pyroptosis was analyzed. CeRNA regulatory network was constructed after predicting the upstream miRNAs and lncRNAs. The validation of hub genes was undertook through RT-qPCR and ELISA. Finally, potential drugs targeting pathways and hub genes were identified, followed by the validation of the effect of two potential drugs on OA.

Results: A total of 161 commom DEGs were obtained, of which 8 genes were finally identified as hub genes through GO and KEGG enrichment analysis as well as PPI network analysis. Eight genes related to ferroptosis and pyroptosis respectively were significantly correlated to the expression of hub genes. 24 miRNAs and 69 lncRNAs were identified to construct the ceRNA regulatory network. The validation of EGR1, JUN, MYC, FOSL1, and FOSL2 met the trend of bioinformatics analysis. Etanercept and Iguratimod reduced the secretion of MMP-13 and ADAMTS5 of fibroblast-like synoviocyte.

Conclusion: EGR1, JUN, MYC, FOSL1, and FOSL2 were identified as hub genes in the development of OA after series of bioinformatics analysis and validation. Etanercept and Iguratimod seemed to have opportunities to be novel drugs for OA.

1 Background

Osteoarthritis (OA) is a common degenerative disease in middle-aged and elderly people all over the world, of which the prevalence has gradually increased due to the aging population and the trend of the overweight (Loeser et al., 2016; Roemer et al., 2022). According to the data from World Health Organization (WHO) in 2019, about 250 million people worldwide suffered from OA (Hunter and Bierma-Zeinstra, 2019). Synovitis is one of the most important pathological manifestations in the occurrence and development of OA, running through the whole process (Sharma, 2021). Studies have shown that synovial lesions usually occur earlier than cartilage, which can be detected by MRI in several patients with small joint injury on X-ray. In addition, low-grade synovitis contributes to radiographic and pain progression (Sanchez-Lopez et al., 2022). The pathological changes of synovitis are complex and diverse, leading to the lack of specific treatment for OA. Therefore, it is necessary to explore the pathogenesis and diagnostic markers of OA from the perspective of synovium to find the therapeutic targets of OA, alleviate symptoms and promote the prognosis.

Ferroptosis is a programmed cell death mode caused by abnormal oxidation, which is regulated by glutathione peroxidase 4 (GPX4), with characteristics of iron-dependent accumulation of lipid peroxide. Kennish et al. (2014) showed that the iron level in serum of OA patients was positively correlated with the severity of OA, suggesting the existence of abnormal iron homeostasis, but its effect in synovium has not been reported. Pyroptosis is another mode of programmed cell death, which occurs when pattern recognition receptors (PRRS) induce the activation of cystine aspartic protein 1 (caspase-1) or caspase-11, which is characterized by the destruction of cell membrane and the release of cytokines. Pyroptosis has been shown to be involved in synovits, that is, Interleukin-1β (IL-1β) in the inflammatory environment derives from synovium rather than cartilage (Borgonio et al., 2014).

Bioinformatics tools have been widely used to process microarray data to determine differentially expressed genes (DEGs) and conduct various analysis. In this study, combining with bioinformatics analysis and verification, we explored and screened the hub genes of OA synovial tissue, and discussed the correlation between them and ferroptosis or pyroptosisy, then constructed the upstream ceRNA regulatory network to evaluate the reliability of these genes as the prediction or treatment targets of OA. Finally, we identified and validated potential drugs targeting pathways and hub genes through Comparative Toxicogenomics Database (CTD), Drugbank and Drugs-Genes interaction (DGI) database.

2 Materials and methods

2.1 Data collection

Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) was used as our data source (Barrett et al., 2013), where “Osteoarthritis” was entered as the keyword in the search box for detection. The selected data set included the gene expression array of osteoarthritis and normal synovial tissue of human samples. Finally, GSE55235 (GPL96 platform), GSE55457 (GPL96 platform) were determined as the data sets of this study.

2.2 Data preprocessing and identification of DEGs

The data of GSE55235, GSE55457 were downloaded with the format of MINiML. The mRNA expression data of OA and normal synovial tissue in each data set were analyzed by limma software package of R software (Yu et al., 2012). P-value was analyzed in GEO to correct false-positive results, which <0.05 and |log2 (fold change)| > 1 was defined as the threshold. DEGs were then obtained, and were visualized by volcano map and heat map. common DEGs were identified after the intersection of the two datasets.

2.3 Enrichment analysis of common DEGs and PPI network construction

Gene Ontology (GO) annotation of the common DEGs was performed by the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 2021) (Sherman et al., 2022), including biological process (BP), cellular component (CC), and molecular function (MF). p < 0.05 was determined as a significant margin for all analysis. Column chart was plotted by “http://www.bioinformatics.com.cn,” a free online platform for data analysis and visualization. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the common DEGs was carried out by R package “clusterprofiler” (Version 4). Hypergeometric test was performed to evaluate the significance of pathway–pathway association: phyper (k-1, M, N-M, n, lower. tail = F). Metascape (http://metascape.org) was used for protein - protein interaction (PPI) network analysis of common DEGs (Zhou et al., 2019). PPI analysis was performed using the following databases: STRING, BioGrid (Kondo et al., 2019), OmniPath and InWeb_IM (Li et al., 2017). Physical score >0.132 was selected to be the standard of gene screening. In addition, molecular complexity detection (MCODE) was used to identify densely connected network components and obtain the hub genes of OA (Bader and Hogue, 2003). The network was visualized using Cytoscape.

2.4 Correlation between OA hub genes and ferroptosis and pyroptosis related genes

19 genes related to ferroptosis and 21 genes related to pyroptosis contained in microarray data were selected for correlation analysis with the hub genes (Liu et al., 2020a; Wu et al., 2021). Due to the small sample size of the three groups of data related to ferroptosis, the samples were combined and normalized (Zhang et al., 2019). The data was standardized using “normalize.quantiles” function in the preprocessCore package of R, and was evaluated through the box diagram. The batch effect of data was evaluated by comparing the visual PCA diagram before and after removal. The polygenic correlation map was displayed by pheatmap package of R. Spearman’s correlation analysis was used to describe the correlation between quantitative variables without normal distribution. p < 0.05 was considered statistically significant.

2.5 Analysis of ceRNA regulatory network of hub genes

ENCORI (Li et al., 2014) (https://starbase.sysu.edu.cn) and TargetScan (McGeary et al., 2019) (https://www. targetscan. org/vert_80/) were used to predict miRNAs which regulated hub genes. The first three reliable miRNAs were selected after the intersection of the prediction. LncBase database v3.0 (Paraskevopoulou et al., 2016) (https://diana.e-ce.uth.gr/lncbasev3/interactions) was used to predict lncRNAs which regulated miRNAs above. The first three reliable lncRNAs were selected for each miRNA. Finally, the ceRNA regulatory network of hub genes was constructed.

2.6 Extraction of human fibroblast-like synoviocyte (FLS)

The study was designed according to the Declaration of Helsinki, and was approved by Ethic Committee of the Second Affiliated Hospital of Harbin Medical University (KY 2021-256). Informed consent was obtained from each donor. Synovium of 3 OA patients (age 54–70 years, Kellgren-Lawrence grade 4) that underwent total joint arthroplasty (TKA) and 3 patients (age 56–68 years) that underwent meniscectomy without OA were obtained at the time after surgery. All patients were confirmed without Rheumatoid Arthritis, acute trauma, tumor or infection of knee joint. Briefly, synovium was cut into pieces at the final size about 0.5 mm*0.5 mm, and put into 0.1% type I collagenase (Biosharp, China, BS163). α-MEM medium (Cytiva, United States, SH30265.01) was added with 10% foetal bovine serum (ExCell Bio, China, FSD500) after 2 h. Primary cells could be seen climbing out after about 3–5 days. FLS of passage 6-8 (P6-8) were used in this study. 10 ng/ml of IL-1β (PEPROTECH, United States, 200-01B) were used to stimulate FLS of OA groups for 48 h in order to imitate the environment of OA, while complete medium was added into the FLS of control group. 10 μg/ml of Etanercept and Iguratimod (MedChem Express, China, HY-108847, HY-17009) were added along with IL-1β to FLS for the following test.

2.7 Screening of potential drugs for OA

CTD database (version 16766M) was used to search drugs targeting the KEGG pathways above. Drugs targeting hub genes were searched and compared using Drugbank (version 5.1.8) database and DGI database (version 4.2.0).

2.8 RT q-PCR

Trizol (Beyotime, China, R0016) was used to extract total RNA from FLS. After determining the concentration, 2 µg of total RNA was used to synthesize cDNA through cDNA synthesis kit. SYBR Green (ES Science, China, QP002) was used for qRT-PCR according to the instructions. The primer sequence of 8 hub genes and GAPDH were listed in Table 1. The mRNA level of a specific gene was calculated as 2-^ΔΔ Ct and normalized to GAPDH.

TABLE 1

TABLE 1. Primers used for RT-qPCR amplification.

2.9 Elisa analysis

The supernatant was collected and centrifuged at 1,000 × g at 4°C for 10 min, then was added to 96-well plates coverd by the antibody of each protein (Elabscience, United States). After incubated for 90 min at 37°C, Biotinylated Detection Ab was added to each well. After 60 min, the plates were washed for 3 times. HRP Conjugate was then added to each well and incubated for 30 min at 37°C in the dark. 90 μL of Substrate Reagent and 50 μL Stop Solution were then added. The optical density (OD value) of each well were determined at once by microplate reader at 450 nm.

2.10 Detection of intracellular ROS

FLSs of each group were seeded at a density about 2000 cells per well in 96-well plates. The medium was replaced by H2DCFDA (Biosharp, China, BL714A) working solution with the concentration of 10 μM. FLSs were then incubated at 37°C for 2 h in the dark and washed twice to fully remove the probes outside the cells. OD value of each well was determined at once by microplate reader at 525 nm.

2.11 Detection of lipid peroxidation

FLSs of each group were broken by ultrasonic cell crusher and centrifuged at 3,500 xg for 10 min at 4°C. MDA standard solution (Solarbio, China, BC0025) and the supernatant to be tested with chromogenic agent were boiled, centrifuged, and added to a 96-well plate respectively. The absorbance was measured at 450 nm, 520 nm, and 600 nm for detection by microplate reader. Then the concentration of MDA was calculated according to the instruction.

2.12 Cell viability assay

FLSs were seeded in 96-well plates at a density at about 2000 cells per well. After adding 10 ul of CCK-8 reagent (Beyotime, China, C0037) at 24 h or 48 h respectively, the plates were put into the incubator for 30 min away from light. The absorbance of each well was measured at 450 nm by microplate reader.

2.13 Statistical analysis

All the data of validation was presented as the means ± standard deviation (SD). Statistical analysis was performed using GraphPad Prism version 6.02. Differences in numerical data between two groups were determined by t-test, while four groups were determined by One-way ANOVA followed by a Bonferroni post-hoc test. p < 0.05 was defined as statistically significant.

3 Results

3.1 Data collection and identification of DEGs

The flow chart of the study design was showed in Figure 1. A total of 20 OA synovium samples and 20 normal synovium samples were obtained from two GEO database, including 10 OA samples and 10 normal samples from GSE55235 dataset and 10 OA samples and 10 normal samples from GSE55457 dataset. Relevant clinical data were not provided in the two datasets above. After screening the differential expression of mRNA between OA and normal synovium samples, 595 up-regulated DEGs and 346 down-regulated DEGs from GSE55235 as well as 175 up-regulated DEGs and 385 down-regulated DEGs from GSE55457 were obtained (Figures 2A, B). After the intersection, a total of 161 common DEGs were obtained (Figure 2C; Table 2).

FIGURE 1

FIGURE 1. Flow chart of the study design. GSE55235, GSE55457 were determined as the data sets of this study, from where we screened DEGs by the standard of p < 0.05 and log2|FC| ≥ 1. After the intersection, 71 up-regulated DEGs and 88 down-regulated DEGs were obtained, followed by GO analysis, PPI network construction and KEGG pathway analysis. The correlation between hub genes and ferroptosis or pyroptosis and the construction of ceRNA were undertaken. Moreover, the expression of 8 hub genes was validated. Finally, potential drugs targeting pathways and hub genes were identified and validated. DEGs = differentially expressed genes; FC = fold change; GO = gene ontology; PPI = protein - protein interaction; KEGG = Kyoto Encyclopedia of Genes and Genomes; ceRNA = competing endogenous RNA.

FIGURE 2

FIGURE 2. Data collection and identification of DEGs. (A) Volcano plot of DEGs in two datasets (p < 0.05 and log2|FC| ≥ 1). The red nodes represent up-regulated genes, and the blue nodes represent down-regulated genes. (B) Heat map of DEGs in two datasets, in which different colors represent different expression trend. Due to the large number of DEGs, the top 50 up-regulated and down-regulated genes with the largest FC were shown. (C) Venn diagram of two datasets which showed up-regulated and down-regulated genes respectively. DEGs = differentially expressed genes; FC = fold change.

TABLE 2

TABLE 2. 161 common DEGs of GSE 55235 and GSE 55457.

3.2 Enrichment analysis and protein-protein interaction analysis of DEGs in OA

In order to further explore the function of DEGs in OA, GO annotation and KEGG enrichment analysis were performed subsequently. As the results, 156 GO terms with p < 0.05 were screened out (S1). Top 5 BP, 5 CC and 4 MF terms were showed in Figure 3A. It showed that the most significant BP of 161 commom DEGs was “positive regulation of transcription from RNA polymerase II promoter”, the most significant CC was “extracellular space”, and the most significant MF was “transcriptional activator activity, RNA polymerase II transcription regulatory region sequence-specific binding”. As the result of KEGG enrichment analysis, 38 pathways with p < 0.05 were found (S2), and top 15 pathways were showed in Figure 3B. It showed that these genes were involved in pathways including tumor necrosis factor (TNF) signaling pathway, osteoclast differentiation, IL-17 signaling pathway, nuclear factor kappa-B (NF-κB) signaling pathway and so on. In addition, as the result of PPI analysis combined with MCODE component algorithm, top 3 modules with highest MCODE score were showed in Figures 4A, B. Finally, 8 hub genes in the biggest module were identified as hub genes, which were ATF3, EGR1, FOSB, FOSL1, FOSL2, JUN, JUNB, and MYC (Table 3).

FIGURE 3

FIGURE 3. Enrichment of GO and KEGG of 161 common DEGs. (A) Go annotation categories of 161 common DEGs including top 5 of BP, top 5 of CC and top 4 of MF. The horizontal axis represents counts of enriched genes. (B) Top 15 KEGG pathway categories after enrichment analysis of 161 common DEGs. The horizontal axis represents counts of enriched genes. GO = gene ontology; BP = biological processes; CC = cellular components; MF = molecular function; KEGG = Kyoto Encyclopedia of Genes and Genomes. (*represents “RNA polymerase II transcription regulatory region sequence-specific binding”, which was hided due to the space).

FIGURE 4

FIGURE 4. PPI analysis of 161 DEGs and screening of key modules. (A) Overview of the PPI network of 161 DEGs through metascape. The larger size of the points, the higher degree of the genes. (B) Genes in three key modules with the highest MCODE value.

TABLE 3

TABLE 3. Features and functions of 8 hub genes in OA screened from DEGs.

3.3 Correlation analysis between OA hub genes and ferroptosis related genes

We combined and normalized the data of 2 groups as described above, then evaluated the data standardization and batch effect through the box diagram (Figure 5A) and visual PCA diagram (Figure 5B), respectively. 19 ferroptosis related genes and 21 pyroptosis related genes were selected for correlation analysis and displayed by heat map. It showed that in the analysis of ferroptosis, ALOX15, CISD1, SAT1, and TFRC were positively correlated to hub genes, while ATP5MC3, GPX4, HSPB1, and MT1G were negatively correlated (Figure 5C). In the analysis of pyroptosis, Caspase-6, ELANE, GSDMB, IL-6 and NLRP1 were positively correlated to hub genes, while GPX4, NOD1 and PYCARD were negatively correlated (Figure 5D).

FIGURE 5

FIGURE 5. Correlation between OA hub genes and ferroptosis or pyroptosis related genes. (A) The data standardization evaluated by box diagram after the combination and normalization of 2 groups. (B, C) The visual PCA diagram before (B) and after (C) removal of the batch. (D, E) The heat map of the correlation between hub genes and genes related to ferroptosis and pyroptosis. The horizontal and vertical coordinates represent genes. Different colors represent the correlation coefficient (red represents positive correlation, blue represents negative correlation). PC = Principal Component.

3.4 Analysis of ceRNA regulatory network of OA hub genes

A total of 24 reliable miRNAs that regulated mRNA of 8 hub genes were obtained, after the intersection of the prediction, using ENCORI and TargetScan database. Then, a total of 69 reliable lncRNAs were obtained after the analysis under LncBase v3.0 database and selection, which regulated miRNAs above. Finally, we constructed the ceRNA regulatory network of hub genes (Figure 6).

FIGURE 6

FIGURE 6. The ceRNA regulatory network of hub genes. The red circles represent hub genes, the blue squares represent miRNAs that targeted hub genes, and the green V-shapes represent lncRNAs that targeted miRNAs.

3.5 Validation of hub genes

qRT-PCR and ELISA were used to evaluate the mRNA and protein expression level of hub genes in FLS of control group and OA group. The results showed that the mRNA expression levels of JUN, MYC, EGR1, FOSL1, and FOSL2 of FLS in OA group was significantly lower than that in normal FLS (p < 0.001, p < 0.01, p < 0.01, p < 0.05, p < 0.05, respectively), meeting the trend of bioinformatics analysis (Figures 7A–E). The expression of JUNB and FOSB had no significant difference with control group (Figures 7F, G), while the expression of ATF3 was higher than that in normal FLS (p < 0.001, Figure 7H), which was opposite to bioinformatics analysis. The result of ELISA showed that the expression of EGR1, JUN, MYC were all significantly lower in the supernatant of FLS in OA group, same to the result of qRT-PCR (p < 0.0001, Figure 7I).

FIGURE 7

FIGURE 7. Validation of hub genes through qRT-PCR and ELISA of FLS. (A–H) mRNA expression of JUN, MYC, EGR1, FOSL2, FOSL1, JUNB, FOSB and ATF3. JUN, MYC, EGR1, FOSL2 and FOSL1 showed the same trend as the results of bioinformatics analysis. JUNB and FOSB showed no significant difference between two groups, while MYC showed the opposite trend. (I) The result of ELISA analysis of EGR1, JUN and MYC in the supernatant of FLS. OA = osteoarthritis. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, ^****p < 0.0001. Error bars represent SD.

3.6 Identification of potential drugs for OA

Drugs targeting 8 hub genes and top 4 KEGG pathways were identified by CTD, DrugBank and DGI database. As the result, 85 drugs targeting 4 of hub genes were obtained, among which Quinapril targeted MYC and JUN at the same time (Figure 8A). 285 drugs were found targeting KEGG pathways, among which Elsubrutinib, Lenalidomide, Lenercept, Mifamurtide, Opinercept, Rebimastat, Tasonermin, Thalidomide and Tibulizumab targeted two pathways. Etanercept and Iguratimod targeted three pathways, which were selected for following validation (Table 4; Figure 8B).

FIGURE 8

FIGURE 8. Identification of potential drugs for OA. (A) Drugs targeting hub genes. Red circles represent genes, blue circles represent drugs. (B) Drugs targeting top four KEGG pathways. Red circles represent pathways, blue circles represent drugs. Green circles represent drugs which have three targets.

TABLE 4

TABLE 4. Drugs targeting KEGG pathways ≥2.

3.7 Etanercept and Iguratimod protect FLS of OA from inflammation and cartilage degeneration

To verify the effect of Etanercept and Iguratimod on FLS, we co-treated FLS with IL-1β. By the test of cell viability, we found that Etanercept can significantly inhibit the growth rate of FLS compared with OA group (p < 0.01, Figure 9A). The results of intracellular ROS and MDA detection showed that the levels of these two species in Etanercept and Iguratimod groups were significantly lower than those of OA group (p < 0.01, p < 0.05, Figures 9B, C). At the same time, the mRNA expression of EGR1, JUN and MYC in Etanercept and Iguratimod group were significantly higher than those in OA group (p < 0.001, p < 0.01 for EGR1, p < 0.001, p < 0.05 for JUN and MYC, Figures 9D–F). Moreover, the level of MMP-13 and ADAMTS5 in the supernatant of FLS was significantly lower in Etanercept and Iguratimod group than those in OA group (p < 0.01 for MMP-13, p < 0.001, p < 0.01 for ADAMTS5, Figures 9G, H).

FIGURE 9

FIGURE 9. Validation of the effect of Etanercept and Iguratimod on FLS. (A) Cell viability of FLS at 24 and 48 h. (B) Relative Intracellular ROS concentration of FLS. (C) MDA concentration of FLS of different groups. (D–F) mRNA expression of EGR1, JUN and MYC of the four groups. (G, H) The level of MMP-13 and ADAMTS5 in the supernatant of FLS in different groups. OA = osteoarthritis; Etan = Etanercept; Igu = Iguratimod; ROS = reactive oxygen species; MDA = Malonic dialdehyde; MMP-13 = matrix metalloproteinase 13; ADAMTS5 = a disintegrin and metalloproteinase with thrombospondin motifs-5. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, ^****p < 0.0001. Error bars represent SD.

4 Discussion

As one of the major chronic diseases endangering the middle-aged and elderly, the affected population of OA is also showing a trend of younger age (Mahmoudian et al., 2021). The compliance of non-drug treatment such as kinesitherapy and physiotherapy is usually inexact. The operations for KOA mainly includes arthroscopic debridement, high tibial osteotomy (HTO), unicondylar knee arthroplasty (UKA) and TKA (Zhang et al., 2020), which still have the problems of high risk and cost. The gathered evidence suggests that mononuclear infiltration and over expression of inflammatory mediators in synovium are seen in early OA and predate radiographic damage in OA (Sokolove and Lepus, 2013). Synovium and synovial fluid are the main contributors to inflammation that can secrete key cytokines, most of which are the main regulators of matrix metalloproteinases (MMP), resulting in the loss of cartilage (Kulkarni et al., 2021). Therefore, in-depth study of the molecular mechanism of synovitis is of great significance to the prevention and treatment for OA.

In this study, we screened out 8 hub genes in synovim that may be involved in the progression of OA through a series of bioinformatics methods. EGR1 is potentially involved in postnatal bone biology and implicated in the regulation of osteoblastic cell growth and differentiation. Studies have shown that when EGR1 was inhibited, the bone mass of the limbs decreased in mice, with lower bone volume fraction and mineral density, same to our study (Reumann et al., 2011). However, over-expression of EGR1 regulated the expression of KLF5 and β-catenin signaling pathway, leading to the acceleration of cartilage hypertrophy and degeneration (Sun et al., 2019). MYC encodes transcription factors to regulate transcriptional activity and cell proliferation, growth, and had also been found to regulate OA process through multiple pathways. Wu et al. (O'Donnell et al., 2005) found that the effects of miR-24 on OA chondrocytes may be achieved by targeting MYC and further regulating the MAPK signaling pathway. The silence of c-MYC could promote proliferation of rat chondrocyte (Zou et al., 2018). ATF3 deficiency in chondrocytes had been reported to alleviates OA development (Iezaki et al., 2016). This effect may be achieved by ATF3 directly affecting the expression of MMP-13 thus reducing cartilage loss (Chan et al., 2017). However, There is little study on the mechanism of ATF3 in OA synoviocyte, except for the bioinformatics analysis.

The rest five hub genes, which were FOSB, FOSL1, FOSL2, JUN and JUNB, belonged to the AP-1 family. AP-1 was thought to be critically involved in the pathogenesis of arthritis due to the binding activity for its cognate recognition sites in the promoters of inflammatory cytokines and matrix-degrading enzymes (Huber et al., 2019). Study showed that mitochondrial dysfunction triggered a catabolic response in chondrocytes via activation of the JNK/AP1 pathway (Ansari et al., 2020). On the other hand, JUNB-FBXO21-ERK axis promoted cartilage degeneration in OA by inhibiting autophagy (Lin et al., 2021).

As independent forms of programmed death, ferroptosis and pyroptosis have been shown to play essential roles in the pathological processes of tumor, Alzheimer’s disease, cerebral hemorrhage, ischemia-reperfusion injury, OA and so on (Kenny et al., 2019; Liang et al., 2019; Liu et al., 2020b). In this study, we selected several ferroptosis and pyroptosis related genes. Although these genes were not differentially expressed in the two datasets, analyzing their relationship with hub genes will be still beneficial to explore the mechanism of FLS participating in OA. We found that FOSL1 and JUNB were associated with 3 of the 19 ferroptosis related genes. Overexpression of c-JUN inhibits ferroptosis induced by erastin in Schwann cells to promotes the rehabilitation of facial nerve function (Gao et al., 2022). Vertically, the strongest two genes related to hub genes were CISD1 and HSPB1. CISD1, an iron-containing outer mitochondrial membrane protein, inhibits ferroptosis by protection against mitochondrial lipid peroxidation (Yuan et al., 2016). The phosphorylation of HSPB1 at the Ser-15 site induced by Erastin is the key of the protective response to ferroptosis stress (Wang et al., 2022). In addition, over-expression of HSPB1 attenuates ferroptosis in rats through promoting G6PD expression (Dai and Hu, 2022). Therefore, we believe that CISD1 and HSPB1 are worth studying in the future to explore the mechanism of ferroptosis in OA synovium. In addition, we found that JUN was related to the most pyroptosis related genes. Study had showed that the Jun N-terminal kinases (JNK) pathway is largely upstream of the NLRP3 inflammasome, which exerts a crucial regulatory impact on microglia pyroptosis and inflammatory responses (He et al., 2021). However, the relationship between JUN and pyroptosis of FLS still needs further study. As a key regulator of innate immunity and inflammasome activation, caspase-6 promotes the activation of pyroptosis (Zheng et al., 2020).

Studies have found the significance of miRNAs and lncRNAs on OA. Xu et al. (2021) showed that SNHG7 ameliorated the development of OA by suppressing apoptosis through miR-214-5p-PPARGC1B-PPARγ axis. Wang et al. (2019) found that lncRNA FOXD2 Adjacent Opposite Strand RNA 1 (FOXD2-AS1) served as the protector for OA patients by inducing chondrocyte proliferation. The downregulation of lncRNA LOC101928134, which acts as a promoter of OA, can block the process of OA (Yang et al., 2019). Here, we built a ceRNA network after predicting the upstream miRNAs and lncRNAs, in order to provide more possible early biomarkers for diagnosis, or targeting drugs for OA patients.

Among the 11 drugs targeting two or more than two KEGG pathways, 10 acted on NF-κB signaling pathway, 9 acted on TNF signaling pathway, 3 acted on Osteoclast differentiation pathway, and 2 acted on IL-17 signaling pathway, suggesting that NF-κB was still the key point for OA treatment. Moreover, Etanercept and Iguratimod related to three pathways. Etanercept, a soluble fusion protein that binds TNF-α, has been proved to be an effective choice targeting several inflammatory diseases, especially rheumatoid arthritis (Graudal et al., 2015). Recent study showed that Etanercept lead to the repairment of cartilage with a special scaffold (Campos et al., 2022), and the reduction of MMP-3 (Kroon et al., 2020), suggesting the potential effect on OA. The effect of Iguratimod on OA still needs further study. Our study proved that both Etanercept and Iguratimod could reduce level of intracellular ROS and MDA of FLS as well as level of MMP-13 and ADAMTS5 in the supernatant. In addition, Etanercept could significantly alleviate the abnormal proliferation of FLS under OA environment. Drugs targeted to the hub genes can be divided into antineoplastic, antipsychotic, anti-inflammatory, anti-hyperglycemia, anti-hyperlipidemia, etc. Their roles in OA still needs to be explored.

Our study still had some limitations. First, we obtained gene arrays in synovial tissues of OA patients from GEO database. However, due to the lack of information of patients, it was hard to correlate the DEGs and hub genes obtained in this study with the severity or grade of OA. Second, the sample size of each group was not so large, which was needed for further research.

In summary, our study aimed to identify key genes involved in the pathophysiology of OA. 161 common DEGs and 8 hub genes were screened through GO analysis, KEGG pathway analysis, PPI network construction as well as MCODE, which may become potential targeting clinical diagnosis and treatment of OA. Furthermore, 5 of 8 genes met the similar expression trend with our result through the validation. Subsequently, we discussed the relationship between hub genes and key genes of ferroptosis and pyroptosis. Moreover, miRNAs and lncRNAs were identified to construct the ceRNA regulatory network of hub genes. Finally, we found that Etanercept and Iguratimod, as top two of potential drugs targeting KEGG pathways, had the protective effect on FLS in the OA environment.

5 Conclusion

EGR1, JUN, MYC, FOSL1, and FOSL2 were identified and validated as hub genes in the development of OA after series of bioinformatics analysis. They may have effect on OA development through different kinds of pathways, including the process of ferroptosis and pyroptosis. Etanercept and Iguratimod seemed to have stronger opportunities to be novel drugs for OA.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by Ethic Committee of the Second Affiliated Hospital of Harbin Medical University (KY2021-256). The patients/participants provided their written informed consent to participate in this study.

Author contributions

DL, PK, and JY: conception and design. WX, XL, and XW: data collection. DL, CX, and WX: data analysis and interpretation. BW, WX, and DL: manuscript writing and validation. All authors have approved the final manuscript.

Funding

This work was supported by the Research & Practice Innovation Program of Harbin Medical University (No. 2020-KYYWF-1447) and National Natural Science Foundation of China (No. 82072472).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1117713/full#supplementary-material

References

Ansari, M. Y., Ahmad, N., Voleti, S., Wase, S. J., Novak, K., Haqqi, T. M., et al. (2020). Mitochondrial dysfunction triggers a catabolic response in chondrocytes via ROS-mediated activation of the JNK/AP1 pathway. J. Of Cell Sci. 133 (22), jcs247353. doi:10.1242/jcs.247353

CrossRef Full Text | Google Scholar

Bader, G. D., and Hogue, C. V. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinforma. 4, 2. doi:10.1186/1471-2105-4-2

CrossRef Full Text | Google Scholar

Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: Archive for functional genomics data sets-update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Borgonio, C., González-Huerta, N. C., Romero-Córdoba, S., Hidalgo-Miranda, A., and Miranda-Duarte, A. (2014). Altered expression of circulating microRNA in plasma of patients with primary osteoarthritis and in silico analysis of their pathways. PloS one 9 (6), e97690. doi:10.1371/journal.pone.0097690

PubMed Abstract | CrossRef Full Text | Google Scholar

Campos, Y., Fuentes, G., Almirall, A., Que, I., Schomann, T., Chung, C. K., et al. (2022). The incorporation of Etanercept into a porous tri-layer scaffold for restoring and repairing cartilage tissue. Pharmaceutics 14 (2), 282. doi:10.3390/pharmaceutics14020282

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, C., Macdonald, C., Litherland, G., Wilkinson, D. J., Skelton, A., Europe-Finner, G. N., et al. (2017). Cytokine-induced MMP13 expression in human chondrocytes is dependent on activating transcription factor 3 (ATF3) regulation. J. Of Biol. Chem. 292 (5), 1625–1636. doi:10.1074/jbc.M116.756601

CrossRef Full Text | Google Scholar

Dai, Y., and Hu, L. (2022). HSPB1 overexpression improves hypoxic-ischemic brain damage by attenuating ferroptosis in rats through promoting G6PD expression. J. Of Neurophysiology 128 (6), 1507–1517. doi:10.1152/jn.00306.2022

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhillon, A. S., and Tulchinsky, E. (2015). FRA-1 as a driver of tumour heterogeneity: A nexus between oncogenes and embryonic signalling pathways in cancer. Oncogene 34, 4421–4428. doi:10.1038/onc.2014.374

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, D., Huang, Y., Sun, X., Yang, J., Chen, J., and He, J. (2022). Overexpression of c-Jun inhibits erastin-induced ferroptosis in Schwann cells and promotes repair of facial nerve function. J. Of Cell. And Mol. Med. 26 (8), 2191–2204. doi:10.1111/jcmm.17241

CrossRef Full Text | Google Scholar

Graudal, N., Jürgens, G., and Combe, B. (2015). Etanercept tapering in rheumatoid arthritis. N. Engl. J. Of Med. 372 (5), 489–490. doi:10.1056/NEJMc1414787

CrossRef Full Text | Google Scholar

Guo, J., Shen, S., Liu, X., Ruan, X., Zheng, J., Liu, Y., et al. (2019). Role of linc00174/miR-138-5p (miR-150-5p)/FOSL2 feedback loop on regulating the blood-tumor barrier permeability. Mol. Ther. Nucleic Acids. 18, 1072–1090. doi:10.1016/j.omtn.2019.10.031

PubMed Abstract | CrossRef Full Text | Google Scholar

He, J., Liu, T., Li, Y., Mi, X., Han, D., Yang, N., et al. (2021). JNK inhibition alleviates delayed neurocognitive recovery after surgery by limiting microglia pyroptosis. Int. Immunopharmacol. 99, 107962. doi:10.1016/j.intimp.2021.107962

PubMed Abstract | CrossRef Full Text | Google Scholar

Huber, R., Kirsten, H., Näkki, A., Pohlers, D., Thude, H., Eidner, T., et al. (2019). Association of human FOS promoter variants with the occurrence of knee-osteoarthritis in a case control association study. Int. J. Mol. Sci. 20 (6), 1382. doi:10.3390/ijms20061382

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, D. J., and Bierma-Zeinstra, S. (2019). Osteoarthr. Lancet 393 (10182), 1745–1759. doi:10.1016/S0140-6736(19)30417-9

CrossRef Full Text

Iezaki, T., Ozaki, K., Fukasawa, K., Inoue, M., Kitajima, S., Muneta, T., et al. (2016). ATF3 deficiency in chondrocytes alleviates osteoarthritis development. J. Of Pathology 239 (4), 426–437. doi:10.1002/path.4739

CrossRef Full Text | Google Scholar

James, C. G., Woods, A., Underhill, T. M., and Beier, F. (2006). The transcription factor ATF3 is upregulated during chondrocyte differentiation and represses cyclin D1 and A gene transcription. BMC Mol. Biol. 7, 30. doi:10.1186/1471-2199-7-30

PubMed Abstract | CrossRef Full Text | Google Scholar

Kennish, L., Attur, M., Oh, C., Krasnokutsky, S., Samuels, J., Greenberg, J. D., et al. (2014). Age-dependent ferritin elevations and HFE C282Y mutation as risk factors for symptomatic knee osteoarthritis in males: A longitudinal cohort study. BMC Musculoskelet. Disord. 15, 8. doi:10.1186/1471-2474-15-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kenny, E. M., Fidan, E., Yang, Q., Anthonymuthu, T. S., New, L. A., Meyer, E. A., et al. (2019). Ferroptosis contributes to neuronal death and functional outcome after traumatic brain injury. Crit. Care Med. 47 (3), 410–418. doi:10.1097/CCM.0000000000003555

PubMed Abstract | CrossRef Full Text | Google Scholar

Kondo, K., Ohigashi, I., and Takahama, Y. (2019). Thymus machinery for T-cell selection. Int. Immunol. 31 (3), 119–125. doi:10.1093/intimm/dxy081

PubMed Abstract | CrossRef Full Text | Google Scholar

Kroon, F., Bay-Jensen, A. C., Wittoek, R., VerbruGGen, G., Smolen, J. S., Kloppenburg, M., et al. (2020). Etanercept therapy leads to reductions in matrix metalloproteinase-3 in patients with erosive hand osteoarthritis. Scand. J. Of Rheumatology 49 (2), 167–168. doi:10.1080/03009742.2019.1657493

CrossRef Full Text | Google Scholar

Kulkarni, P., Martson, A., Vidya, R., Chitnavis, S., and Harsulkar, A. (2021). Pathophysiological landscape of osteoarthritis. Adv. Clin. Chem. 100, 37–90. doi:10.1016/bs.acc.2020.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97. doi:10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Wernersson, R., Hansen, R. B., Horn, H., Mercer, J., Slodkowicz, G., et al. (2017). A scored human protein-protein interaction network to catalyze genomic interpretation. Nat. METHODS 14 (1), 61–64. doi:10.1038/nmeth.4083

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, C., Zhang, X., Yang, M., and Dong, X. (2019). Recent progress in ferroptosis inducers for cancer therapy. Adv. Mater. 31 (51), e1904197. doi:10.1002/adma.201904197

PubMed Abstract | CrossRef Full Text | Google Scholar

Licht, A. H., Pein, O. T., Florin, L., Hartenstein, B., Reuter, H., Arnold, B., et al. (2006). JunB is required for endothelial cell morphogenesis by regulating core-binding factor beta. J. Of Cell Biol. 175 (6), 981–991. doi:10.1083/jcb.200605149

CrossRef Full Text | Google Scholar

Lin, Z., Miao, J., Zhang, T., He, M., Wang, Z., Feng, X., et al. (2021). JUNB-FBXO21-ERK axis promotes cartilage degeneration in osteoarthritis by inhibiting autophagy. Aging Cell 20 (2), e13306. doi:10.1111/acel.13306

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, P., Feng, Y., Li, H., Chen, X., Wang, G., Xu, S., et al. (2020). Ferrostatin-1 alleviates lipopolysaccharide-induced acute lung injury via inhibiting ferroptosis. Cell Mol. Biol. Lett. 25 (10), 10. doi:10.1186/s11658-020-00205-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Zhao, Q., Zuo, Z. X., Yuan, S. Q., Yu, K., Zhang, Q., et al. (2020). Systematic analysis of the aberrances and functional implications of ferroptosis in cancer. iScience 23 (7), 101302. doi:10.1016/j.isci.2020.101302

PubMed Abstract | CrossRef Full Text | Google Scholar

Loeser, R. F., Collins, J. A., and Diekman, B. O. (2016). Ageing and the pathogenesis of osteoarthritis. Nat. Rev. Rheumatol. 12 (7), 412–420. doi:10.1038/nrrheum.2016.65

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, H., Hou, G., Zhang, Y., Dai, Y., and Zhao, H. (2014). c-Jun transactivates Puma gene expression to promote osteoarthritis. Mol. Med. Rep. 9 (5), 1606–1612. doi:10.3892/mmr.2014.1981

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahmoudian, A., Lohmander, L. S., Mobasheri, A., Englund, M., and Luyten, F. P. (2021). Early-stage symptomatic osteoarthritis of the knee - time for action. Nat. Rev. Rheumatol. 17 (10), 621–632. doi:10.1038/s41584-021-00673-4

PubMed Abstract | CrossRef Full Text | Google Scholar

McGeary, S. E., Lin, K. S., Shi, C. Y., Pham, T. M., Bisaria, N., Kelley, G. M., et al. (2019). The biochemical basis of microRNA targeting efficacy. Science 366 (6472), eaav1741. doi:10.1126/science.aav1741

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Donnell, K. A., Wentzel, E. A., Zeller, K. I., Dang, C. V., and Mendell, J. T. (2005). c-Myc-regulated microRNAs modulate E2F1 expression. Nature 435 (7043), 839–843. doi:10.1038/nature03677

PubMed Abstract | CrossRef Full Text | Google Scholar

Papoudou-Bai, A., Hatzimichael, E., Barbouti, A., and Kanavaros, P. (2017). Expression patterns of the activator protein-1 (AP-1) family members in lymphoid neoplasms. Clin. And Exp. Med. 17 (3), 291–304. doi:10.1007/s10238-016-0436-z

CrossRef Full Text | Google Scholar

Paraskevopoulou, M. D., Vlachos, I. S., Karagkouni, D., Georgakilas, G., Kanellos, I., Vergoulis, T., et al. (2016). DIANA-LncBase v2: Indexing microRNA targets on non-coding transcripts. Nucleic Acids Res. 44 (D1), D231–D238. doi:10.1093/nar/gkv1270

PubMed Abstract | CrossRef Full Text | Google Scholar

Reumann, M. K., Strachna, O., Yagerman, S., Torrecilla, D., Kim, J., Doty, S. B., et al. (2011). Loss of transcription factor early growth response gene 1 results in impaired endochondral bone repair. BONE 49 (4), 743–752. doi:10.1016/j.bone.2011.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Roemer, F. W., Guermazi, A., Hannon, M. J., Fujii, T., Omoumi, P., Hunter, D. J., et al. (2022). Presence of magnetic resonance imaging- defined inflammation particularly in overweight and obese women increases risk of radiographic knee osteoarthritis: The POMA study. Arthritis Care & Res. 74 (8), 1391–1398. doi:10.1002/acr.24568

CrossRef Full Text | Google Scholar

Sanchez-Lopez, E., Coras, R., Torres, A., Lane, N. E., and Guma, M. (2022). Synovial inflammation in osteoarthritis progression. Nat. Rev. Rheumatol. 18 (5), 258–275. doi:10.1038/s41584-022-00749-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, L. (2021). Osteoarthritis of the knee. N. Engl. J. Of Med. 384 (1), 51–59. doi:10.1056/NEJMcp1903768

CrossRef Full Text | Google Scholar

Sherman, B. T., Hao, M., Qiu, J., Jiao, X., Baseler, M. W., Lane, H. C., et al. (2022). David: A web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 50, W216–W221. doi:10.1093/nar/gkac194

PubMed Abstract | CrossRef Full Text | Google Scholar

Sokolove, J., and Lepus, C. M. (2013). Role of inflammation in the pathogenesis of osteoarthritis: Latest findings and interpretations. Ther. Adv. Musculoskelet. Dis. 5 (2), 77–94. doi:10.1177/1759720X12467868

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Huang, H., Pan, X., Li, S., Xie, Z., Ma, Y., et al. (2019). EGR1 promotes the cartilage degeneration and hypertrophy by activating the Krüppel-like factor 5 and β-catenin signaling. Biochimica Biophysica Acta-Molecular Basis Of Dis. 1865 (9), 2490–2503. doi:10.1016/j.bbadis.2019.06.010

CrossRef Full Text | Google Scholar

Tanel, A., Fonseca, S. G., Yassine-Diab, B., Bordi, R., Zeidan, J., Shi, Y., et al. (2009). Cellular and molecular mechanisms of memory T-cell survival. Expert Rev. Vaccines 8 (3), 299–312. doi:10.1586/14760584.8.3.299

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Wu, S., He, H., Ai, K., Xu, R., Zhang, L., et al. (2022). CircRNA-ST6GALNAC6 increases the sensitivity of bladder cancer cells to erastin-induced ferroptosis by regulating the HSPB1/P38 axis. Lab. Investig. 102 (12), 1323–1334. doi:10.1038/s41374-022-00826-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Cao, L., Wang, Q., Huang, J., and Xu, S. (2019). LncRNA FOXD2-AS1 induces chondrocyte proliferation through sponging miR-27a-3p in osteoarthritis. Artif. Cells Nanomedicine Biotechnol. 47 (1), 1241–1247. doi:10.1080/21691401.2019.1596940

CrossRef Full Text | Google Scholar

Wu, J. B., Zhu, Y., Luo, M. M., and Li, L. (2021). Comprehensive analysis of pyroptosis-related genes and tumor microenvironment infiltration characterization in breast cancer. Front. Immunol. 12, 748221. doi:10.3389/fimmu.2021.748221

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., Pei, Y., Lu, J., Liang, X., Li, Y., Wang, J., et al. (2021). LncRNA SNHG7 alleviates IL-1β-induced osteoarthritis by inhibiting miR-214-5p-mediated PPARGC1B signaling pathways. Int. Immunopharmacol. 90, 107150. doi:10.1016/j.intimp.2020.107150

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, D. W., Zhang, X., Qian, G. B., Jiang, M. J., Wang, P., Wang, K. Z., et al. (2019). Downregulation of long noncoding RNA LOC101928134 inhibits the synovial hyperplasia and cartilage destruction of osteoarthritis rats through the activation of the Janus kinase/signal transducers and activators of transcription signaling pathway by upregulating IFNA1. J. Of Cell. Physiology 234 (7), 10523–10534. doi:10.1002/jcp.27730

CrossRef Full Text | Google Scholar

Yu, G., Wang, L., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics-A J. Of Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

Yuan, H., Li, X., Zhang, X., Kang, R., and Tang, D. (2016). CISD1 inhibits ferroptosis by protection against mitochondrial lipid peroxidation. Biochem. And Biophysical Res. Commun. 478 (2), 838–844. doi:10.1016/j.bbrc.2016.08.034

CrossRef Full Text | Google Scholar

Zhang, X., Zhang, W., Jiang, Y., Liu, K., Ran, L., and Song, F. (2019). Identification of functional lncRNAs in gastric cancer by integrative analysis of GEO and TCGA data. J. Of Cell. Biochem. 120 (10), 17898–17911. doi:10.1002/jcb.29058

CrossRef Full Text | Google Scholar

Zhang, Z., Huang, C., Jiang, Q., Zheng, Y., Liu, Y., Liu, S., et al. (2020). Guidelines for the diagnosis and treatment of osteoarthritis in China (2019 edition). Ann. Transl. Med. 8 (19), 1213. doi:10.21037/atm-20-4665

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, M., Karki, R., Vogel, P., and Kanneganti, T. D. (2020). Caspase-6 is a key regulator of innate immunity, inflammasome activation, and host defense. Cell 181 (3), 674–687. doi:10.1016/j.cell.2020.03.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, J., Li, X., Shi, Z., and Xue, J. F. (2018). Effects of C-myc gene silencing on interleukin-1β-induced rat chondrocyte cell proliferation, apoptosis and cytokine expression. J. Of Bone And Mineral Metabolism 36 (3), 286–296. doi:10.1007/s00774-017-0845-4

CrossRef Full Text | Google Scholar