Data sources
From the GEO database (https://www.ncbi.nlm.nih.gov/geo/) to download osteoarthritis and normal cartilage transcriptome data set (GSE12021, GSE55235 and GSE55457). The GSE12021 dataset contains 9 normal and 22 OA samples, and the GSE55235 dataset contains 10 normal and 20 OA samples. The GSE55457 dataset contains 10 normal and 23 OA samples. To reduce error bias for individual samples, we combined and corrected the GSE12021, GSE55235 and GSE55457 datasets using the "LIMMA" and "SVA" R packages. A dataset of chromatin regulatory factors was obtained from Lu et al. [10] (Table S1). Then, by integrating the above three GEO datasets and chromatin regulatory factor datasets, OA gene expression datasets related to chromatin regulatory factors were finally obtained.
Identification of CRs-related DEGs
In order to reveal the difference in CRs-related gene expression between normal subjects and osteoarthritis patients, differential expression analysis of CRs-related genes in OA dataset was performed using the "Limma" language package of R software, and adj.p < 0.05 and |log2FC|>0.5 were used as the threshold values for screening DEGs [19]. At the same time, the "ggplot17" and "pheatmap" R language packages were used to create volcano maps and heat maps, respectively, to visualize differentially expressed genes.
Functional and pathway enrichment analysis
In order to understand the biological mechanisms and signaling pathways of CRs-related genes in OA, the "ClusterProfiler" language package of R software was used to conduct function enrichment analysis of CRs-related differentially expressed genes based on gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [20]. As part of the GO analysis, we explored the biological functions of CRS-related DEGs, including biological processes (BP), cellular components (CC), and molecular functions (MF). The GO term or KEGG path with adj.p <0.05 was considered statistically significant.
Construction of PPI network and identification of hub genes
The Search Tool for the Retrieval of Interacting Genes database (STRING, https://string-db.org/) was used for PPI analysis of CRs-related differentially expressed genes [21], and Cytosscape software was used to visualize the STRING analysis data [22]. The MCC algorithm of Cytohubba plug-in in Cytoscape was used to identify hub genes, and genes with top 10 scores were considered as potential hub genes [23].
Immune cell infiltration and correlation analysis
SsGSEA method was used to uncoil the expression matrix of 16 human immune cell subtypes and 13 immune functions, so as to explore the composition of immune cells and immune functions in OA patients and normal controls [24]. We then assessed the correlation between immune cells and immune function in patients with OA and healthy controls, and finally developed a heat map. Further, we explored the differences in infiltration levels of various immune cells and immune functions in OA patients and healthy controls, and performed rank sum tests through the "vioplot" R language package, which was finally demonstrated as a violin diagram. Finally, in order to understand the association of these immune cells and immune function with signature gene expression in OA patients and healthy controls, we used the Spearman method to analyze the correlation with the first 5 hub genes.
Construction and evaluation of signature genes risk model for OA
To further analyze the diagnostic value of these 5 hub genes in OA, we established a logistic regression risk model using CRs expression matrix and sample grouping as input files. Then, using the "rms" R language package, we plotted the OA risk model prediction nomogram, calculated the risk score for 5 hub genes, and validated it with ROC curve and calibration curve to measure the nomogram model's recognition ability. The calibration curve is displayed graphically, and the 45° line represents the best predicted value.
Cell culture and the establishment of osteoarthritis model in vitro
ATDC5 chondrocytes at a density of 1×106 cells/well were maintained in regular DMEM/F12 supplemented with 10% FBS and 1% penicillin/streptomycin and incubated in a humidified atmosphere of 5% CO2 at 37 °C. To construct an osteoarthritis model in vitro, ADTC5 cells were stimulated with IL-1β (10ng/ml) for 24h.
Total RNA extraction and real-time quantitative polymerase chain teaction (RT-qPCR)
Total RNA was respectively extracted from chondrocytes and cartilage tissues by TRIzol reagent (Takara, Japan) following the manufacturer's instruction. The purity and concentration of extracted RNA was measured by a spectrophotometer and 1000ng RNA was used to synthesize cDNA. Then, real-time PCR supermix kit (Mei5bio) was used for amplification. Finally, relative mRNA expression was calculated using the 2−ΔΔCT method. ACTB was used to normalize as an internal control. The primer sequences were: AURKB sense, 5ʹ-GATCCCAGAACAAGCAGCCT -3ʹ; antisense, 5ʹ-TCGATTTCGATCTCTCGGCG-3ʹ.
Western blot analysis
Total proteins harvested from ATDC5 cells and articular cartilage using RIPA lysis buffer (Applygen, China) containing a mixture of phosphatase and protease inhibitor (Beyotime, China). The supernatants were collected after lysates were centrifuged at 12,000 × g at 4℃ for 10 min. Subsequently, the bicinchoninic acid assay (BCA) was used to qualify the protein concentration. Protein samples were separated by 10% SDS-PAGE and transferred to polyvinylidene (PVDF) membranes (Millipore, USA). After blocking with 5% skim milk for two hours, membranes were incubated with AURKB (1:1000, Boster, China) or actin (1:5000, proteintech, China) overnight at 4°C. These membranes were washed with Tris-buffered saline-Tween 20 (TBST)and incubated with the corresponding secondary antibodies for 2 h at room temperature on the next day. Eventually, the protein bands were detected by the Bio-Rad Imaging System. Actin was used as an internal control.
Immunofluorescence analysis
After treatment, ADTC5 were washed with PBS for three times and then fixed with 4% paraformaldehyde for 15 min at the room temperature. After permeabilization with 0.1% Triton X-100 for 15 min, the cells were incubated with bovine serum albumin for 30 min.. Subsequently, the cells were incubated with rabbit anti-AURKB antibody (1:200) overnight at 4°C, followed by incubation with the revelant secondary antibody conjugated with Alexa Fluor 488 for 1 h and stained with 4, 6-diamidino-2- phenylindole (DAPI) for 10 min at room temperature in dark. Finally, the stained cells were observed by Olympus microscope.
DMM-induced OA mouse model
All animal experiment procedures were strictly performed following the instruction of the Animal Ethics Committee of the First Affiliated Hospital of Nanchang University. All methods in the experiment were carried out in accordance with relevant guidelines and regulations. Osteoarthritis model was established by surgical destabilization of the medial meniscus (DMM) [25].In brief, C57BL/6 male mice (8-week-old, n=6) were anesthetized with intraperitoneal injection of 10% chloral hydrate. Subsequently, a microsurgical scissor was used to incise the joint capsule of the right knee to the medial meniscotibial ligament. An arthrotomy without the transaction of medial meniscotibial ligament in the right knee was considered as sham group. After surgery, these mice were divided into two groups: Sham and DMM groups. All mice were given the same amount of water and food every day. Finally, mice were sacrificed at 8 weeks post-OA surgery, and the samples of knee joints were dissected and collected for further evaluation.
Potential therapeutic drug prediction
Using DSigDB protein database (http://tanlab.ucdenver.edu/DSigDB) - drug interaction data to predict the potential drugs in the treatment of OA, with adj.p < 0.05 and the comprehensive score > 5000 as the threshold.
Statistical analysis
Statistical analysis was performed using GraphPad Prism 8 and R software 4.2.2. Differences between groups were compared using student t tests. Spearman test is used to analyze the correlation between two variables. To examine the predictive value of key genes for OA and healthy specimens, the ROC curve and AUC were calculated using the R packet "proc". All statistical tests were two-tailed and the statistical significance level was set to 0.05.