Weighted gene co-expression network analysis reveals specific modules and hub genes related to immune infiltration of osteoarthritis

Introduction

Osteoarthritis (OA) is a degenerative disease that is caused by numerous factors, including aging, obesity, strain, trauma, congenital abnormality, joint deformity, and so on (1,2). OA is common in middle-aged and elderly people, especially in weight-bearing joints (knee and hip). With the aging of the population, the incidence of OA is increasing every year, which becomes a leading cause of disability and confers a heavy burden on the both patient’s family and society (3). OA is recognized as having a complicated pathological process and a long course of disease. At present, diagnosis based on clinical manifestations is the standard to confirm OA, while early diagnosis and treatment intervention remain challenging (4,5). Some studies have shown that biomarkers provide useful diagnostic information by detecting cartilage degradation in OA, reflecting the biological activities related to the disease, and predicting the process of disease progression. Therefore, exploring diagnostic biomarkers of OA may have important clinical application value (6,7).

In recent years, an increasing number of studies have shown that immune cell infiltration plays an important role in the development of OA (8,9). Pro-inflammatory cytokines secreted by infiltrated immune cells can aggravate chondrocyte apoptosis and cartilage matrix proteolysis (10,11). However, the pathological mechanism of immune cell infiltration in the process of OA is still unclear. Although exploring the pathological mechanism of osteoarthritis can provide potential therapeutic targets for treatment (12,13), it is more significant to explore the biomarkers for early diagnosis (14). Therefore, evaluating the infiltration of immune cells and determining the different composition of infiltrating immune cells might aid in the development of new diagnostic biomarkers and immunotherapy targets.

With the evolution of high-throughput sequencing technology, microarray has been used in a large number of OA studies (15-17). Weighted gene co-expression network analysis (WGCNA) involves the construction of a network based on systematic gene expression levels and analyzes gene expression microarray profiling data, with the aim of determining the co-expression relationship between genes. Thus, genes with similar expression patterns may be co-regulated, functionally related, or in the same pathway (18,19). In this study, we used microarray data from OA patients and normal control subjects to analyze differentially expressed genes (DEGs) and immune cell infiltration by WGCNA, in order to identify effective biomarkers for the early diagnosis in order to control the progress of OA. We present the following article in accordance with the MDAR reporting checklist (available at https://dx.doi.org/10.21037/atm-21-4566).

Methods

Microarray data sources and processing

The series matrix file, GSE51588, was downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) public database (20). There were 40 groups of human OA data (including 20 LT groups and 20 MT groups) and 10 groups of human non-OA data (including five LT groups and five MT groups) for the WGCNA co-expression network construction. The series matrix file data file, GSE55235, was also downloaded from the NCBI GEO public database (21). There were data from 20 transcriptional groups, including human non-OA groups (n=10) and human OA groups (n=10), which were used for subsequent model validation. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

DEG identification and WGCNA

We screened out the DEGs using the “limma” R package (https://bioconductor.org/packages/release/bioc/html/limma.html) with P<0.05 & |log2 fold change (FC)| >2 (22). By building a weighted gene co-expression network, we could find the gene modules of co-expression, and explore the relationship between the gene network and phenotype, as well as the core genes in the network. We used the R package ‘WGCNA’ (http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA) to build the co-expression network of all genes in the GSE51588 dataset (19). The first 5,000 genes were screened by this algorithm for further analysis; the soft threshold was set to 12, and the soft threshold of GSE55235 was set to 4. The weighted adjacency matrix was converted into the topological overlap matrix (TOM) in order to estimate the network connectivity, and a hierarchical clustering method was used to generate a cluster tree structure of the TOM matrix (23). Different branches of the cluster tree represented the different gene modules, and different colors represented the different modules. In view of the weighted correlation coefficient of genes, the genes were classified based on the expression patterns. Genes showing similar patterns were divided into one same module, and thousands of genes were divided into multiple modules via the gene expression patterns. Finally, the overlapping hub module genes from GSE51588 and GSE55235 were identified by the Venn method.

Immune gene correlation analysis

The effect of genes on immune infiltration was evaluated. The level of immune cell infiltration in each sample was quantified by single sample gene set enrichment analysis (ssGSEA) (24). Spearman correlation analysis was performed on gene expression and immune cell content.

Gene set variation analysis (GSVA) analysis

GSVA is a nonparametric and unsupervised method for assessing transcriptome gene set enrichment (25). GSVA changes the gene level into the pathway level by comprehensively scoring the gene set of interest, and then evaluates the biological function of the sample. In this study, to evaluate the potential biological function changes of different samples, we downloaded gene sets from the Molecular signatures database (v7.4) (http://www.gsea-msigdb.org/gsea/msigdb/), and used the GSVA algorithm to comprehensively score each gene set.

GeneMANIA analysis

GeneMANIA (http://www.genemania.org) is a flexible and user-friendly protein-protein interaction (PPI) network construction database, which is used to visualize the functional network between genes and analyze gene functions and interactions (26). The website can set the data sources of gene nodes, and has a variety of bioinformatics analysis methods, such as physical interaction, gene co-expression, gene co-location, gene enrichment analysis, and website prediction. In this study, geneMANIA generated the core gene network to explore its possible mechanism in patients with OA.

Statistical analysis

All statistical analysis was conducted with R language (version 3.6.3) (https://cran.r-project.org/bin/windows/base/old/3.6.3/). Statistical analysis was performed using the two-tailed Student’s t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

Important DEGs in GSE51588 and WGCNA

The GSE51588 disease-related dataset was downloaded from the NCBI GEO public database to identify DEGs between the normal (n=10) and disease (n=40) groups. We further calculated the DEGs between the two groups using the “limma” R package. DEGs were screened out on the basis of the threshold of P<0.05 & |log2 fold change (FC)| >2. A total of 95 DEGs were screened out, including 24 up-regulated genes and 71 down-regulated genes (Figure 1A). According to the patients’ clinical characteristics, the WGCNA network was further constructed to detect the key modules in OA. The soft threshold β was determined by the “sft$powerEstimate" function and set to 12 (Figure 1B-1F). Based on the TOM, 12 gene modules were detected in this analysis, which were blue [1,166], cyan [186], grey [621], light cyan [692], light green [76], light yellow [59], magenta [262], midnight blue [149], pink [508], purple [553], red [283], and turquoise [445] (Figure 1G). Through module-trait correlations analysis, we found that the pink module had the highest correlation with clinical features (cor =0.69, P=2e-08), and thus, the pink module was selected for subsequent verification (Figure 1H).

Figure 1 Identification of DEGs in the GSE51588 dataset and construction of the WGCNA. (A) Volcano plot of DEGs in GSE51588; (B) clustering dendrogram of 50 samples; (C,D) soft-threshold power for WGCNA; (E,F) dendrogram of DEGs clustered based on a dissimilarity measure (1-TOM); (G) Heatmap of correlation between the ME and gene modules; (H) A scatterplot of gene significance vs. MM in the pink module. DEGs, differentially expressed genes; WGCNA, weighted gene co-expression network analysis; ME, module eigengenes; MM, module membership.

Important DEGs in GSE55235 and WGCNA

We downloaded the GSE55235 disease-related dataset from the NCBI GEO public database to identify DEGs between OA patients (n=10) and controls (n=10) using the “limma” R package. Based on the threshold of P<0.05 & |log2 FC|>2, 164 DEGs were screened out, including 84 up-regulated genes and 80 down-regulated genes (Figure 2A). Similarly, the WGCNA network was constructed to detect the key modules in OA based on the clinical characteristics of patients. The soft threshold β was determined by “sft$powerEstimate” function and set to 4 (Figure 2B-2F). Based on the TOM, eight gene modules were detected, including black [1,173], brown [776], green [363], green yellow [70], magenta [107], purple [79], red [2039], and yellow [393] (Figure 2G). Furthermore, through correlation analysis of the modules and traits, we found that red module had the highest correlation with clinical features (cor =0.94, P=1e-09), and thus, the red module was selected for subsequent verification (Figure 2H).

Figure 2 Identification of DEGs in the GSE55235 dataset and construction of the WGCNA. (A) Volcano plot of DEGs in GSE55235; (B) clustering dendrogram of 20 samples; (C,D) soft-threshold power for WGCNA; (E,F) dendrogram of DEGs clustered based on a dissimilarity measure (1-TOM); (G) Heatmap of correlation between the ME and gene modules; (H) a scatterplot of gene significance vs. MM in the red module. DEGs, differentially expressed genes; WGCNA, weighted gene co-expression network analysis; ME, module eigengenes; MM, module membership.

Identification of hub genes and diagnostic efficacy verification

A Venn diagram was used to detect the core genes in both datasets. Three hub DEGs were identified, including nuclear factor interleukin-3 (NFIL3), adrenomedullin (ADM), and osteoglycin (OGN) (Figure 3A). Receiver operating characteristic (ROC) curve assessment and the area under curve (AUC) value were utilized to verify the diagnostic efficacy. The results showed that the AUC values of the three core genes were greater than 0.8 (NFIL3: AUC =0.957; ADM: AUC =0.900; OGN: AUC =0.815), indicating that these three hub genes could be effective indicators for clear diagnosis (Figure 3B-3D).

Figure 3 Identification of hub genes and diagnostic efficacy verification. (A) The Venn diagram shows the intersecting genes of GSE51588 and GSE55235; (B-D) ROC curve analysis and AUC statistics were applied to verify diagnostic efficacy of NFIL3, ADM, and OGN separately. ROC, receiver operating characteristic; AUC, area under the curve.

GSVA of core genes and key pathways

Next, we explored the potential molecular mechanism and investigated the specific signaling pathways enriched in the three core genes by GSVA. The GSVA results showed that highly-expressed NFIL3 was related to tumor necrosis factor alpha (TNF-α) signaling via nuclear factor kappa-B (NF-κB), the reactive oxygen species (ROS) pathway, and myelocytomatosis (MYC) targets v2 (Figure 4A). Highly-expressed ADM pathways included TNF-α signaling via NF-κB, the ROS pathway, and ultraviolet (UV) response up (Figure 4B). OGN-enriched pathways included epithelial mesenchymal transition, coagulation, and peroxisome (Figure 4C).

Figure 4 GSVA of core genes and key pathways. (A) NFIL3-enriched pathways; (B) ADM-enriched pathways; (C) OGN-enriched pathways. The blue band represents a positive correlation and the green band represents a negative correlation. GSVA, gene set variation analysis; ADM, adrenomedullin; OGN, osteoglycin.

Differences in immune infiltration and PPI analysis

The inflammatory microenvironment of OA is mainly composed of immune cells, extracellular matrix, various growth factors, inflammatory factors, and special physical and chemical characteristics, which significantly affect the sensitivity of disease diagnosis and treatment. By analyzing the relationship between core genes and immune infiltration, we investigated the potential mechanism of core genes affecting disease progression. The results showed that there were strong correlations between the three hub genes and immune cells (Figure 5A), which was in line with expectations. The PPI network involved in the three core genes is shown in Figure 5B.

Figure 5 Immune infiltration and PPI analysis. (A) Correlation between immune cells and hub genes. Red, positive; blue, negative; (B) the PPI network of hub genes. PPI, protein-protein interaction.

Discussion

OA is a common chronic joint disease in the elderly, which is characterized by the degeneration of articular cartilage, thickening of synovial lining, subchondral sclerosis, and formation of osteophyte at the edges of joints (27). There is no effective clinical treatment for OA and the pathological mechanism remains unclear. Therefore, understanding the molecular mechanism and pathological process of OA is critical. A comprehensive investigation of OA, including the pathological process, early clinical diagnosis, clinical manifestations, clinical prevention, and drug treatment, requires systematic analysis. Since normal cartilage samples, as a control group, are limited, microarray data in the GEO database is sparse and the samples are often insufficient. In this study, we obtained microarray data from GSE51588 and GSE55235 to expand sample size and focused on identifying biological markers related to early diagnosis, for the first time.

Compared with other bioinformatics analyses, the advantage of WGCNA lies in the comprehensive study of the relationship between co-expression modules and clinical features, which provides highly reliable and biologically significant results (28). We screened out three hub genes based on the WGCNA and DEGs in this study, including NFIL3, ADM, and OGN. Meanwhile, based on functional enrichment via GSVA in this study, it was suggested that TNF-α signaling via NF-κB and the ROS pathway play important roles in the progression of OA, which has been confirmed in other studies. Chen et al. found that spermine can reduce the progression of OA by inhibiting the TNF-α/NF-κB pathway (29). It has also been reported that TNF-α plays a major pathological role in rheumatoid arthritis (RA) via NF-κB and the Janus kinase (JAK)/signal transducer and activator of transcription (STAT) pathway (30). Numerous studies have shown that ROS or oxidative stress regulates the intracellular signal transduction process, chondrocyte aging and apoptosis, synthesis and degradation of extracellular matrix, synovitis, and subchondral bone dysfunction (31-33).

Macrophages have been reported playing a crucial role in the progression of OA (34). Especially, M1 macrophages, known as pro-inflammatory phenotype, have an adverse effect on chondrocyte apoptosis which leading to cartilage degeneration (10). Besides, T cell, B cell and NK cell, which belong to innate immune system, are also defined as integral cells involved in the process of OA (35-37). Mesenchymal stem cells (MSCs) have pluripotent mesoderm differentiation potential and can differentiate into a variety of cell types, including bone cells, chondrocytes, muscle cells, and adipocytes (38). Studies have shown that MSCs can exert their ability to promote tissue repair by releasing paracrine factors, mainly a variety of growth factors, immunoregulatory cytokines and other nutrient mediators (39,40). It has been reported that MSCs could reduce chondrocyte apoptosis, inhibit the inflammatory response and provide a suitable microenvironment for cartilage repair (41,42). Besides, extracellular vesicles derived from MSCs also show the crucial effect on attenuating cartilage degeneration (43,44).

As for the three DEGs, only one study reported the correlation between NFIL3 and OA, which may be helpful to reveal the similarities and differences in the pathogenesis of femoral head necrosis and cartilage injury in hip OA (45). A series of studies revealed the role of OGN in OA, which suggested that highly-expressed OGN could promote the ossification and reconstruction of articular cartilage (46-48). Other studies suggested that ADM was regulated by TNF-α, and ADM down-regulation could inhibit adipogenesis and synthesis of osteocytes, which may provide new treatment strategies for RA and OA (49-51).

Taken together, these three identified hub genes provide important insights into the potential molecular mechanisms of OA, which may reveal new biomarkers for early diagnosis.

Acknowledgments

Funding: This work was supported by the Clinical Research Fund of Wu Jieping Medical Foundation (Project No. 320.6750.2020-03-13) National Natural Science foundation of China (82072439) and the International (Regional) Cooperation and Exchange Program of the National Natural Science Foundation of China (81620108018).

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at https://dx.doi.org/10.21037/atm-21-4566

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/atm-21-4566). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Martel-Pelletier J, Barr AJ, Cicuttini FM, et al. Osteoarthritis. Nat Rev Dis Primers 2016;2:16072. [Crossref] [PubMed]
2. Hunter DJ, Bierma-Zeinstra S. Osteoarthritis. Lancet 2019;393:1745-59. [PubMed]
3. Sacitharan PK. Ageing and Osteoarthritis. Subcell Biochem 2019;91:123-59. [Crossref] [PubMed]
4. Abramoff B, Caldera FE. Osteoarthritis: Pathology, Diagnosis, and Treatment Options. Med Clin North Am 2020;104:293-311. [Crossref] [PubMed]
5. Maruotti N, Corrado A, Cantatore FP. Osteoblast role in osteoarthritis pathogenesis. J Cell Physiol 2017;232:2957-63. [Crossref] [PubMed]
6. Kraus VB, Karsdal MA. Osteoarthritis: Current Molecular Biomarkers and the Way Forward. Calcif Tissue Int 2021;109:329-38. [Crossref] [PubMed]
7. van Spil WE, Szilagyi IA. Osteoarthritis year in review 2019: biomarkers (biochemical markers). Osteoarthritis Cartilage 2020;28:296-315. [Crossref] [PubMed]
8. Nees TA, Rosshirt N, Zhang JA, et al. T Helper Cell Infiltration in Osteoarthritis-Related Knee Pain and Disability. J Clin Med 2020;9:2423. [Crossref] [PubMed]
9. Cai W, Li H, Zhang Y, et al. Identification of key biomarkers and immune infiltration in the synovial tissue of osteoarthritis by bioinformatics analysis. PeerJ 2020;8:e8390 [Crossref] [PubMed]
10. Utomo L, Bastiaansen-Jenniskens YM, Verhaar JA, et al. Cartilage inflammation and degeneration is enhanced by pro-inflammatory (M1) macrophages in vitro, but not inhibited directly by anti-inflammatory (M2) macrophages. Osteoarthritis Cartilage 2016;24:2162-70. [Crossref] [PubMed]
11. Mobasheri A, Rayman MP, Gualillo O, et al. The role of metabolism in the pathogenesis of osteoarthritis. Nat Rev Rheumatol 2017;13:302-11. [Crossref] [PubMed]
12. Gao X, Sun Y, Li X. Identification of key gene modules and transcription factors for human osteoarthritis by weighted gene co-expression network analysis. Exp Ther Med 2019;18:2479-90. [Crossref] [PubMed]
13. Gu HY, Yang M, Guo J, et al. Identification of the Biomarkers and Pathological Process of Osteoarthritis: Weighted Gene Co-expression Network Analysis. Front Physiol 2019;10:275. [Crossref] [PubMed]
14. Kandahari AM, Yang X, Dighe AS, et al. Recognition of Immune Response for the Early Diagnosis and Treatment of Osteoarthritis. J Immunol Res 2015;2015:192415 [Crossref] [PubMed]
15. Cai P, Jiang T, Li B, et al. Comparison of rheumatoid arthritis (RA) and osteoarthritis (OA) based on microarray profiles of human joint fibroblast-like synoviocytes. Cell Biochem Funct 2019;37:31-41. [Crossref] [PubMed]
16. Liu W, Jiao Y, Tian C, et al. Gene Expression Profiling Studies Using Microarray in Osteoarthritis: Genes in Common and Different Conditions. Arch Immunol Ther Exp (Warsz) 2020;68:28. [Crossref] [PubMed]
17. Rai MF, Sandell LJ, Barrack TN, et al. A Microarray Study of Articular Cartilage in Relation to Obesity and Severity of Knee Osteoarthritis. Cartilage 2020;11:458-72. [Crossref] [PubMed]
18. Stuart JM, Segal E, Koller D, et al. A gene-coexpression network for global discovery of conserved genetic modules. Science 2003;302:249-55. [Crossref] [PubMed]
19. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
20. Yang Z, Ni J, Kuang L, et al. Identification of genes and pathways associated with subchondral bone in osteoarthritis via bioinformatic analysis. Medicine (Baltimore) 2020;99:e22142 [Crossref] [PubMed]
21. Huang Z, He Z, Kong Y, et al. Insight into osteoarthritis through integrative analysis of metabolomics and transcriptomics. Clin Chim Acta 2020;510:323-9. [Crossref] [PubMed]
22. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47 [Crossref] [PubMed]
23. Song L, Langfelder P, Horvath S. Comparison of co-expression measures: mutual information, correlation, and model based indices. BMC Bioinformatics 2012;13:328. [Crossref] [PubMed]
24. Xu J, Zeng Y, Si H, et al. Integrating transcriptome-wide association study and mRNA expression profile identified candidate genes related to hand osteoarthritis. Arthritis Res Ther 2021;23:81. [Crossref] [PubMed]
25. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
26. Warde-Farley D, Donaldson SL, Comes O, et al. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res 2010;38:W214-20 [Crossref] [PubMed]
27. Johnson VL, Hunter DJ. The epidemiology of osteoarthritis. Best Pract Res Clin Rheumatol 2014;28:5-15. [Crossref] [PubMed]
28. Pei G, Chen L, Zhang W. WGCNA Application to Proteomic and Metabolomic Data Analysis. Methods Enzymol 2017;585:135-58. [Crossref] [PubMed]
29. Chen Z, Lin CX, Song B, et al. Spermidine activates RIP1 deubiquitination to inhibit TNF-α-induced NF-κB/p65 signaling pathway in osteoarthritis. Cell Death Dis 2020;11:503. [Crossref] [PubMed]
30. Nishimura R, Hata K, Takahata Y, et al. Role of Signal Transduction Pathways and Transcription Factors in Cartilage and Joint Diseases. Int J Mol Sci 2020;21:1340. [Crossref] [PubMed]
31. Lepetsos P, Papavassiliou AG. ROS/oxidative stress signaling in osteoarthritis. Biochim Biophys Acta 2016;1862:576-91. [Crossref] [PubMed]
32. Drevet S, Gavazzi G, Grange L, et al. Reactive oxygen species and NADPH oxidase 4 involvement in osteoarthritis. Exp Gerontol 2018;111:107-17. [Crossref] [PubMed]
33. Bolduc JA, Collins JA, Loeser RF. Reactive oxygen species, aging and articular cartilage homeostasis. Free Radic Biol Med 2019;132:73-82. [Crossref] [PubMed]
34. Fernandes TL, Gomoll AH, Lattermann C, et al. Macrophage: A Potential Target on Cartilage Regeneration. Front Immunol 2020;11:111. [Crossref] [PubMed]
35. Haseeb A, Haqqi TM. Immunopathogenesis of osteoarthritis. Clin Immunol 2013;146:185-96. [Crossref] [PubMed]
36. de Lange-Brokaar BJ, Ioan-Facsinay A, van Osch GJ, et al. Synovial inflammation, immune cells and their cytokines in osteoarthritis: a review. Osteoarthritis Cartilage 2012;20:1484-99. [Crossref] [PubMed]
37. Huss RS, Huddleston JI, Goodman SB, et al. Synovial tissue-infiltrating natural killer cells in osteoarthritis and periprosthetic inflammation. Arthritis Rheum 2010;62:3799-805. [Crossref] [PubMed]
38. Hu L, Yin C, Zhao F, et al. Mesenchymal Stem Cells: Cell Fate Decision to Osteoblast or Adipocyte and Application in Osteoporosis Treatment. Int J Mol Sci 2018;19:360. [Crossref] [PubMed]
39. Mamidi MK, Das AK, Zakaria Z, et al. Mesenchymal stromal cells for cartilage repair in osteoarthritis. Osteoarthritis Cartilage 2016;24:1307-16. [Crossref] [PubMed]
40. Barry F, Murphy M. Mesenchymal stem cells in joint disease and repair. Nat Rev Rheumatol 2013;9:584-94. [Crossref] [PubMed]
41. Zhou J, Wang Y, Liu Y, et al. Adipose derived mesenchymal stem cells alleviated osteoarthritis and chondrocyte apoptosis through autophagy inducing. J Cell Biochem 2018; [Epub ahead of print]. [PubMed]
42. Ryu JS, Jung YH, Cho MY, et al. Co-culture with human synovium-derived mesenchymal stem cells inhibits inflammatory activity and increases cell proliferation of sodium nitroprusside-stimulated chondrocytes. Biochem Biophys Res Commun 2014;447:715-20. [Crossref] [PubMed]
43. Woo CH, Kim HK, Jung GY, et al. Small extracellular vesicles from human adipose-derived stem cells attenuate cartilage degeneration. J Extracell Vesicles 2020;9:1735249 [Crossref] [PubMed]
44. Wang Y, Yu D, Liu Z, et al. Exosomes from embryonic mesenchymal stem cells alleviate osteoarthritis through balancing synthesis and degradation of cartilage extracellular matrix. Stem Cell Res Ther 2017;8:189. [Crossref] [PubMed]
45. Wang W, Liu Y, Hao J, et al. Comparative analysis of gene expression profiles of hip articular cartilage between non-traumatic necrosis and osteoarthritis. Gene 2016;591:43-7. [Crossref] [PubMed]
46. Ruiz AR, Tuerlings M, Das A, et al. The role of TNFRSF11B in development of osteoarthritic cartilage. Rheumatology (Oxford) 2021. [Epub ahead of print].
47. Yan CL, An FY, Wang YF, et al. Intervention effect of Youguiwan on rats with knee osteoarthritis and its mechanism. Zhongguo Ying Yong Sheng Li Xue Za Zhi 2020;36:511-6. [PubMed]
48. Balakrishnan L, Nirujogi RS, Ahmad S, et al. Proteomic analysis of human osteoarthritis synovial fluid. Clin Proteomics 2014;11:6. [Crossref] [PubMed]
49. Hopwood B, Tsykin A, Findlay DM, et al. Gene expression profile of the bone microenvironment in human fragility fracture bone. Bone 2009;44:87-101. [Crossref] [PubMed]
50. Uzan B, Ea HK, Launay JM, et al. A critical role for adrenomedullin-calcitonin receptor-like receptor in regulating rheumatoid fibroblast-like synoviocyte apoptosis. J Immunol 2006;176:5548-58. [Crossref] [PubMed]
51. Qi X, Yu F, Wen Y, et al. Integration of transcriptome-wide association study and messenger RNA expression profile to identify genes associated with osteoarthritis. Bone Joint Res 2020;9:130-8. [Crossref] [PubMed]