Construction of a diagnostic model for osteoarthritis based on transcriptomic immune-related genes

Abstract

Background: Osteoarthritis (OA) is a leading cause of disability globally, affecting over 500 million individuals worldwide. However, accurate and early diagnosis of OA is challenging to achieve. Immune-related genes play an essential role in OA development. Therefore, the objective of this study was to develop a diagnostic model for OA based on immune-related genes identified in synovial membrane.

Methods: The gene expression profile of OA were downloaded based on four datasets. The significantly differentially expressed genes (DEGs) between OA and control groups were selected. The differential immune cells were analyzed, followed by immune-related DEGs screening. WGCNA was used to screen module genes and these genes were further selected through optimization algorithm. Then, nomogram model was constructed. Chemical drug small molecule related to OA was predicted. Finally, expression levels of several key genes were validated by qRT-PCR through construction of OA rat models.

Results: The total 656 DEGs were obtained. Eight immune cells were significantly differential between two groups, and 317 immune-related DEGs were obtained. WGCNA identified three modules. The genes in modules were significantly involved in 15 pathways, involving in 65 genes. Then 12 DEGs were screened as the final optimal combination of DEGs, such as CEBPB, CXCL1, JUND, GABARAPL2 and PDGFC. The Nomogram model was also constructed. Furthermore, the chemical small molecules, such as acetaminophen, aspirin, and caffeine were predicted. The expression levels of CEBPB, CXCL1, GABARAPL2 and PDGFC were validated in OA rat models.

Conclusion: A diagnostic model based on twelve immune related genes was constructed. These model genes, such as CEBPB, CXCL1, GABARAPL2, and PDGFC, may serve as diagnostic biomarkers and immunotherapeutic targets.

Keywords: Biomarker; Gene; Model; Osteoarthritis; immune.

Fig. 1

Screening of significantly differentially expressed genes (DEGs). A and B: The sample relationships before (A) and after (B)batch effect removal, respectively. C and D: The (C) volcano map and (D) heatmap of DEGs. The blue and red dots in the volcano map indicate significantly downregulated and upregulated genes, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2

Percentage of immune cells. A: The composition distribution of each immune cell in each sample. B: Eight immune cells were significantly differential between OA and control groups.

Fig. 3

Screening of DEGs significantly associated with important immune cells. A: Network diagram of immune cells and DEGs. Circles indicate DEGs and colors indicate degree of difference; Red squares indicate eight immune cell types. B: GO functions and KEGG signaling pathways enriched by immune-related genes. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4

The constructed interaction network. The color indicates the degree of significant difference and the size of the node indicates the degree of connectivity. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5

Screening of modules related to disease progression and immunity. A: Left: The selection process for the adjacency matrix weight parameter power. The red line signifies the standard line where the squared value of the correlation coefficient reaches 0.9. On the right, the average connectivity of genes is displayed under different power parameters. B: The tree diagram visually presents the division of modules, with each color representing a distinct module. C: The MDS diagram showcases the genes contained within each module. D: A heatmap demonstrates the correlation between the proportion of immune cells in the sample and the modules obtained through partitioning. E: The enriched GO functions and KEGG signaling pathways attributed to module genes. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6

Construction and verification of a Nomogram diagnostic model. A–C: Parameter diagram of feature gene screening based on (A) LASSO, (B) RFE and (C) RF algorithms. D: Nomogram model based on the expression levels of 12 feature genes in the combined training dataset. E: The calibration curves of Nomogram model of the combined training dataset. The horizontal and vertical axes represented the calculated and actual probabilities, respectively. The solid line depicted the performance of an ideal diagnostic model, while the dotted line illustrated the performance of the developed nomogram model. F: Nomogram model based on the expression levels of 12 feature genes in the validation dataset. G: The calibration curves of Nomogram of the validation dataset.

Fig. 7

Verification of key genes in the nomogram model. The expression levels of 12 genes in the training dataset (A) and GSE82107 validation dataset (B). (C) The expression levels of CEBPB, CXCL1, JUND, GABARAPL2 and PDGFC in OA group compared with control detected by qRT-PCR. *P < 0.05, **P < 0.01, ***P < 0.001.

Fig. 8

The network of the 12 feature immune-related genes and small molecular chemical drugs.