Identification and validation of a glutamine metabolism-related gene signature as a diagnostic and therapeutic target in osteoarthritis through integrated multi-omics analysis

Data acquisition and analysis strategy

We implemented a multi-stage integrative strategy, comprising four main phases: (1) identification of glutamine metabolism-related differentially expressed genes; (2) functional analysis of candidate genes and diagnostic model construction; (3) mechanistic exploration of diagnostic genes; and (4) experimental validation using independent clinical samples. The detailed study design workflow is illustrated in Fig. 1.

Dataset selection was based on the following criteria to ensure quality and relevance: (i) Inclusion of transcriptomic or epigenomic data from both OA patients and healthy controls. While multiple joint types were considered, priority was given to datasets focusing on knee joint pathology; (ii) Availability of complete raw data from key joint-related tissues, specifically articular cartilage, synovium, subchondral bone, meniscus, or peripheral blood components; (iii) Availability of clear clinical annotations, including but not limited to Kellgren-Lawrence (K-L) radiographic grading for all included knee OA samples; and (iv) sufficient sample size for statistical analyses.

Exclusion criteria included: (i) absence of healthy or non-OA control group; (ii) insufficient methodological descriptions necessary for data reprocessing; and (iii) use of non-standardized platforms or tissues irrelevant to joint pathology. Public data from GEO (https://www.ncbi.nlm.nih.gov/geo/) are summarized in Table S1, including the corresponding tissue type, joint anatomical location, and the analytical purpose assigned to each dataset within this study.

Clinical sample collection

To experimentally validate our bioinformatics findings, we prospectively recruited an independent clinical cohort for RT-qPCR and metabolomics analyses. The study protocol received formal approval from the Institutional Review Board at Chongqing Hospital of Traditional Chinese Medicine (CQTCM, Approval No. 2022-KY-YJS-LY) and adhered to the Declaration of Helsinki principles. All participants provided written informed consent before their inclusion in the study.

A priori power analysis (G*Power) indicated that a minimum of 30 samples per group was required to detect an effect size of 0.6 with 80% power at α = 0.05, based on prior OA biomarker studies demonstrating the utility of power analysis in this field (Attur et al., 2020). Informed by this analysis, we prospectively recruited a cohort comprising 31 patients with primary KOA and 31 demographically matched healthy volunteers between November 2022 and March 2024. Inclusion was based on American College of Rheumatology (ACR) criteria (Altman et al., 1986) with a Kellgren-Lawrence (K-L) radiographic grade ≥2. Exclusion criteria comprised secondary OA, inflammatory arthropathies, intra-articular injections within the past 6 months, and major systemic/immune disorders.

Following an overnight fast (≥8 h), 10 mL of venous blood was collected from each participant into EDTA-containing vacutainer tubes for subsequent RT-qPCR and metabolomic analyses. Peripheral blood mononuclear cells (PBMCs) were isolated immediately via density gradient centrifugation using Ficoll-Paque PLUS (GE Healthcare, Chicago, IL, USA), washed twice with phosphate-buffered saline (PBS), and either processed directly or cryopreserved at −80 °C for ≤6 months prior to RNA extraction. Key clinical parameters, including age, sex, BMI, symptom duration, Visual Analog Scale (VAS) pain scores, Western Ontario and McMaster Universities Osteoarthritis Index (WOMAC) scores, and K-L radiographic grades, were systematically recorded for all subjects. Any samples with missing clinical or metabolomic data (>5%) were excluded from the final analysis.

Data preprocessing and differential expression analysis

Raw expression data were normalized using R (version 4.1.2). Microarray datasets were processed via Robust Multi-array Average (RMA) algorithm (Fang, Li & Xu, 2021) for background correction, quantile normalization, and log₂ transformation. RNA-seq data were normalized using the TMM method. To ensure data integrity, missing values were inspected and handled using gene-wise filtering, whereby probes or genes with excessive missingness were removed, and sporadic missing entries were imputed using k-nearest-neighbor (kNN) imputation, a widely adopted approach for high-dimensional transcriptomic data (Troyanskaya et al., 2001). Outlier samples were identified through boxplot inspection and PCA-based multivariate distance analysis. Pre-analysis quality control also included assessment of batch effects using principal component analysis (PCA). To mitigate batch-associated variation across datasets, batch correction was performed using the ComBat algorithm from the sva package (v3.42.0) (Tran et al., 2020). The effectiveness of ComBat adjustment was confirmed by comparing PCA distributions before and after correction. Differential expression analysis was conducted using limma (v3.48.3) (Smyth, 2004) with thresholds of raw P < 0.05 and |log2FC| > 0.5, adopting raw P-values to maximize sensitivity. Data visualization was performed using the pheatmap and ggplot2 packages.

Identification of glutamine metabolism-related genes

Glutamine metabolism-related genes were collected through an integrated analysis of gene sets from the MSigDB database (KEGG_GLUTAMINE_METABOLISM); (Liberzon et al., 2015) and the GeneCards database (keyword: “glutamine metabolism”, relevance score ≥5). GeneCards relevance scores were used as described in the database’s documentation. The collected GMRGs were subsequently cross-referenced with differentially expressed genes (DEGs) to derive a candidate list of glutamine metabolism-associated differentially expressed genes (DE-GMRGs). A multi-step sensitivity analysis was then performed on this candidate gene list. This analysis involved: (i) cross-referencing with core gene sets from REACTOME (REACTOME_GLUTAMATE_AND_GLUTAMINE_METABOLISM.v2025.1.Hs.gmt) and KEGG (hsa00250) databases; (ii) varying the GeneCards relevance scores (from 3 to 10) and using alternative source databases; and (iii) applying different |log₂FC| thresholds (from 0.3 to 1.0) for the identification of DEGs.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA R package (v1.70-3) (Langfelder & Horvath, 2008) on the GSE57218 dataset. Sample quality was first assessed by average linkage hierarchical clustering, and one outlier sample (GSM1380887) was removed to improve network stability. To determine the optimal network construction parameters, a soft-thresholding power was selected by evaluating scale-free topology fit across powers 1–20. A power of 10 was chosen, as it achieved a scale-free topology model fit (R² ≈ 0.87) while maintaining adequate mean connectivity. A signed co-expression network was then constructed, and gene modules were identified using dynamic tree cutting (deepSplit = 2, minModuleSize = 30). Module eigengenes (MEs) were calculated, and their correlations with the OA phenotype were assessed using Pearson correlation to identify disease-associated modules.

Functional enrichment and network analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using clusterProfiler (v4.0.5) (Wu et al., 2021). Multiple testing correction was applied using the Benjamini-Hochberg method, and terms or pathways with adjusted P-value < 0.05 were considered significant. For completeness, enrichment terms with 0.05 ≤ adjusted P-value < 0.10 were also reported as trend-level findings. Within the Gene Ontology framework, three primary domains were analyzed: biological processes (BP), cellular components (CC), and molecular functions (MF).

Gene functional similarity was calculated using GOSemSim (v2.20.0) (Yu et al., 2010) with the Wang metric and Best-Match Average (BMA) aggregation. Pairwise GO-based semantic similarity scores were then computed for the 13-WGCNA-derived DE-GMRGs to generate a similarity matrix, and each gene’s overall similarity was obtained by averaging its scores across all others. Genes with higher cumulative similarity were considered more functionally central (“Friends”), whereas those with lower similarity were regarded as functionally distinct.

Diagnostic model construction and evaluation

We developed and validated a diagnostic signature using a two-stage workflow integrating WGCNA-derived candidate genes, LASSO-based feature selection, and multivariable logistic regression (LR). First, hub genes identified by WGCNA from glutamine metabolism-associated differentially expressed genes were subjected to LASSO logistic regression (glmnet R package) (Tibshirani, 1996) for feature selection, using an independent GEO dataset for discovery (GSE57218: n = 40). A fixed random seed (set.seed (12345)) was used for reproducibility. The regularization parameter (λ) was selected via 10-fold cross-validation by choosing λ_min that minimized the mean cross-validated binomial deviance; genes with non-zero coefficients at λ_min were retained as candidate features. Second, the selected genes were used to train a multivariable LR model on an independent cohort (GSE48556: n = 139). Internal performance was evaluated on GSE48556 by: (i) discrimination using ROC/AUC with 95% CIs estimated by the DeLong method; (ii) calibration using bootstrap calibration curves with 1,000 resamples; (iii) clinical utility using Decision Curve Analysis (DCA); and (iv) clinical visualization via a nomogram (rms package). A decision threshold was derived from the training cohort using Youden’s index and kept fixed thereafter. For external validation, the trained LR model and the predetermined threshold were applied—without parameter retraining—to multiple independent GEO cohorts (GSE55457, GSE114007, GSE51588, and GSE29746), using the same preprocessing and normalization/scaling procedures as in the training dataset.

Integrated analysis of diagnostic gene mechanisms

DNA methylation data from the GSE73626 (n =18) dataset were analyzed using β-values downloaded from GEO. Differential methylation analysis between OA and control samples was performed using the limma package (v1.40.0) in R. For each CpG site, linear modeling and empirical Bayes moderation were applied, and differentially methylated positions (DMPs) were identified using thresholds of P < 0.05 and |Δβ| > 0.2 (Campagna et al., 2021). CpG sites were annotated to corresponding genes using the UCSC hg19 reference, allowing multiple CpG sites to map to the same gene. DMPs located within diagnostic genes were extracted, and their β-values were visualized to assess methylation alterations at the gene level.

To identify biological pathways associated with each diagnostic gene, we performed correlation-based GSEA using clusterProfiler (v4.0.5). For each diagnostic gene, Spearman correlation coefficients between its expression and all other genes in the GSE57218 dataset were calculated to generate ranked gene lists. This approach identifies genes whose expression patterns co-vary with the diagnostic gene, and may reflect shared regulatory mechanisms. The ranked lists were subjected to GSEA against MSigDB gene sets (v7.5.1), including GO terms (c5.all) and KEGG pathways (c2.cp.kegg). Pathways with an FDR <0.05 and |normalized enrichment score (NES)| >1 were considered significantly enriched (Subramanian et al., 2005).

Potential therapeutic compounds targeting the diagnostic genes were identified using the Comparative Toxicogenomics Database (CTD, accessed November 2023) and the Drug-Gene Interaction Database (DGIdb v4.2.0). From CTD, we extracted curated chemical-gene interactions with interaction scores >5 limited to Homo sapiens. From DGIdb, we retrieved interactions classified as “inhibitor,” “agonist,” “antagonist,” or “modulator” with evidence from at least two independent sources. Overlapping drugs identified by both databases were prioritized, and the drug-gene interaction network was visualized using Cytoscape (v3.8.2).

ceRNA regulatory network construction

We constructed a ceRNA regulatory network integrating the three diagnostic genes (F13A1, IRS2, RELA) with miRNAs and circRNAs based on the competing endogenous RNA hypothesis (Salmena et al., 2011). Potential miRNA binding sites on diagnostic gene 3′UTRs were predicted using miRWalk v3.0 (Sticht et al., 2018). To balance prediction sensitivity and specificity, miRNA–mRNA pairs were retained when predicted by ≥8/12 algorithms and exhibiting minimum free binding energy ΔG <−25 kcal/mol, consistent with established rules for miRNA–target interactions (Kertesz et al., 2007).

To ensure biological relevance, predicted miRNAs were filtered against the GSE105027 miRNA-expression dataset using limma with thresholds of raw P < 0.05 and |log₂FC| >0.5, retaining only dysregulated miRNAs likely to participate in pathological ceRNA interactions.

For each differentially expressed miRNA, upstream circRNA sponges were predicted using StarBase v3.0 (Li et al., 2014). To increase reliability, circRNAs supported by chipExpNum >30 were retained, in accordance with the experimental evidence filtering strategy of StarBase. This criterion matches the computational pipeline used in the original analysis (Bosson, Zamudio & Sharp, 2014).

Predicted circRNAs were further compared with circRNA expression profiles in GSE175959. Because the dataset contains only 3 OA and 3 control samples, raw P < 0.05 and |log₂FC| >0.5 were applied as hypothesis-generating thresholds, consistent with exploratory ceRNA analysis.

The final integrated circRNA–miRNA–mRNA network was visualized with Cytoscape (v3.8.2) (Shannon et al., 2003), where node shapes represent molecular classes and node size reflects degree centrality. We note that this ceRNA network represents computational predictions requiring experimental validation.

RNA extraction and RT-qPCR

Total RNA was extracted from approximately 2 × 10⁶ PBMCs using a TRIzol-based reagent (SparkJade, Shandong, China) and subsequently treated with DNase I. RNA quality was verified by spectrophotometry (A₂₆₀/A₂₈₀ ratio 1.8–2.0; A₂₆₀/A₂₃₀ >1.8) and by observing intact 18S/28S rRNA bands on a 1.5% agarose gel. One microgram of qualified RNA was reverse-transcribed into cDNA in a 20 μL reaction using the 2×SPARKscript II All-in-One RT SuperMix (SparkJade, Shandong, China), following the thermal profile: 25 °C for 5 min, 42 °C for 15 min, and 85°C for 5 s. Quantitative PCR was performed on an ABI 7500 System using SPARKscript II Master Mix (SYBR Green), with each 20 μL reaction containing two μL of 1:10 diluted cDNA and 400 nM of each primer (Shanghai Bioengineering, Shanghai, China; Table S2). The amplification program consisted of 95 °C for 2 min, followed by 40 cycles of 95 °C for 15 s, 60 °C for 30 s, and 72 °C for 30 s. All reactions were run in triplicate, and product specificity was confirmed by melt curve analysis; no-template and no-reverse-transcriptase controls showed no amplification (Cq > 38). Relative gene expression was calculated using the 2^−ΔΔCT method, for which primer efficiencies were confirmed to be 95–105% via standard curve analysis. To identify suitable reference genes, a panel of six candidates (GAPDH, B2M, YWHAZ, ACTB, RPL13A, and TBP) was selected based on literature precedent (Pombo-Suarez et al., 2008). These genes were chosen to represent diverse biological pathways to minimize the risk of co-regulation. Their expression stability was evaluated across all experimental groups using the geNorm algorithm (Vandesompele et al., 2002), which identified B2M and YWHAZ as the most stable pair (M-value < 0.5) (Table S3). Consequently, normalization was performed against the geometric mean of B2M and YWHAZ. All RT-qPCR experiments were MIQE-compliant, as detailed in Supplemental File MIQE Checklist.

Metabolomics analysis

Plasma samples (100 μL) were thawed on ice and proteins were precipitated with 400 μL of pre-chilled methanol/water (80:20, v/v). The mixture was vortexed (2 min) and centrifuged (12,000 rpm, 10 min, 4 °C). The resulting supernatant was diluted with ultrapure water to a final methanol concentration of 53% to match the initial mobile phase, and then transferred to LC-MS vials. A pooled quality control (QC) sample, prepared by mixing equal aliquots from each sample, and blank samples (extraction solvent) were analyzed periodically throughout the sequence to monitor system stability (Broadhurst et al., 2018).

Metabolomic analysis was performed on a Vanquish UHPLC system coupled to a Q Exactive™ HF-X mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA) with a heated electrospray ionization (HESI) source. A 2 µL aliquot was injected onto a Hypersil Gold C18 column (2.1 × 100 mm, 1.9 μm) maintained at 40 °C. Chromatographic separation was achieved using a mobile phase of (A) 0.1% formic acid in water and (B) 0.1% formic acid in acetonitrile at a flow rate of 0.3 mL/min. The gradient was: 0–2 min, 5% B; 2–12 min, 5–95% B; 12–14 min, 95% B; 14–14.1 min, 95–5% B; 14.1–16 min, 5% B. The mass spectrometer operated in data-dependent acquisition (DDA) mode with polarity switching. Full MS scans were acquired from m/z 100–1,500 at 60,000 resolution (AGC target: 3 × 10⁶, max IT: 100 ms). MS/MS scans for the top 10 most abundant ions were acquired at 15,000 resolution (AGC target: 1e5, max IT: 50 ms) with a 1.0 m/z isolation window and stepped Normalized Collision Energy (NCE) of 20, 30, and 40 eV. Source parameters included a spray voltage of 3.5 kV (positive)/−2.8 kV (negative) and a capillary temperature of 320 °C.

Raw data were processed using Compound Discoverer™ 3.1 software and custom R scripts. The workflow included peak picking, retention time alignment, and normalization. Features with a coefficient of variation (CV) >30% in QC samples were excluded (Schiffman et al., 2019). Metabolite annotation was based on matching accurate mass (<5 ppm tolerance) and MS/MS fragmentation patterns against public (KEGG, HMDB) and in-house spectral libraries (Wishart et al., 2007; Zhou, Wang & Ressom, 2012).

Statistical analysis

Data processing utilized IBM SPSS Statistics (Version 25.0) and R (V4.1.2), initiating with Shapiro–Wilk normality testing to determine analytic pathways. Normally distributed continuous variables were expressed as mean ± SD (independent t-tests), non-normal data as median (IQR) (Mann–Whitney U tests), and categorical variables analyzed by χ² test or Fisher’s exact test. Correlation analyses employed Pearson’s r (parametric) or Spearman’s ρ (non-parametric), categorized as weak (|r/ρ| <0.3), moderate (0.3–0.7), or strong (≥0.7).

For metabolomics data, we initially employed multivariate pattern recognition methods, including Partial Least Squares Discriminant Analysis (PLS-DA) and Orthogonal PLS-DA (OPLS-DA) (Worley & Powers, 2012), to maximize the separation between groups. Differential metabolites were identified using a combination of criteria: a Variable Importance in Projection (VIP) score >1 and a P-value < 0.05 from an independent samples t-test (Wold, Sjöström & Eriksson, 2001). Multiple group comparisons used ANOVA (Bonferroni post-hoc) or Kruskal–Wallis tests. Multivariate linear regression modeled metabolite-clinical relationships with age/sex/BMI adjustment.

All tests were two-tailed with a significance threshold of P < 0.05. For omics analyses involving multiple testing, the Benjamini–Hochberg method was applied to control the false discovery rate. Differential expression analyses used an adjusted P-value (adj. P) < 0.05 as the significance threshold, while pathway enrichment analyses employed a more lenient exploratory threshold of FDR < 0.1 for hypothesis generation.

Identification of glutamine-related differentially expressed genes

To ensure robust data integration, we first performed principal component analysis (PCA) to assess potential batch effects between the GSE57218 and GSE55235 datasets. The application of the ComBat algorithm effectively removed these variations, ensuring reliable comparability between cohorts (Figs. S1A and S1B). Following quality control, analysis of the GSE57218 dataset identified 1,597 DEGs (Fig. 2A), while the GSE55235 dataset yielded 2,279 DEGs (Fig. 2B). Hierarchical clustering demonstrated clear segregation of expression patterns in both datasets (Figs. 2C and 2D). Integration of MSigDB and GeneCards databases yielded 1,190 GMRGs. By intersecting these GMRGs with DEGs from both datasets, we identified a core set of 22 DE-GMRGs (Table S4).

To assess the robustness of this intersection strategy, we performed a supplementary differential analysis on the merged, batch-corrected expression matrix. This joint analysis identified 74 significant genes (Table S4). Notably, all 22 DE-GMRGs identified by our intersection approach were contained within this expanded set, confirming that our strategy successfully captured the stable core gene signature while minimizing false positives. Finally, a comprehensive sensitivity analysis was conducted. The candidate signature remained highly stable across various gene-set definitions and statistical thresholds (Figs. S1C and S1D), demonstrating the robustness of our selection process.

WGCNA was performed to identify co-expression modules associated with OA. A soft-thresholding power of 10 was selected based on scale-free topology fit (R² ≈ 0.87) and mean connectivity (Fig. 3A). Dynamic tree cutting identified 12 co-expression modules (Fig. 3B), among which three modules (magenta: r = −0.45, P = 0.004; blue: r = −0.79, P < 0.001; red: r = −0.58, P < 0.001) showed significant negative correlations with the OA phenotype (Fig. 3C). These modules contained 3,099 genes in total. Intersection of the 22 DE-GMRGs (Figs. 4A and 4B; Table S4) with the 3,099 WGCNA-derived genes yielded 13 candidate genes: F13A1, IL10, CTH, PDK4, IRS2, SLC3A2, ENO2, RELA, VAMP2, GLB1, CYBA, CHST6, and SDC1 (Fig. 4C). Chromosomal localization of these genes is shown in Fig. S1E. All 13 genes exhibited significant differential expression between OA and control samples (Fig. 4D), with F13A1 and RELA showing the most pronounced downregulation in OA.

Functional analysis of candidate genes

Spearman correlation analysis revealed significant co-expression relationships among the 13 candidate genes (Fig. 5A), with the strongest positive correlations observed between IRS2 and PDK4 (ρ = 0.72, P < 0.001), and a notable correlation between RELA and SLC3A2 (ρ = 0.54, P = 0.002). Functional enrichment analysis demonstrated that the 13 candidate genes converge on several biological processes relevant to OA.

GO analysis revealed significant enrichment in pathways related to amino acid metabolism, energy metabolism, peptide hormone response, inflammatory signaling, and oxidative stress regulation (all q < 0.05; Fig. 5B). Key enriched terms included response to insulin (P = 7.8 × 10⁻⁷, q = 3.6 × 10⁻⁴), response to peptide hormone (P = 7.0 × 10⁻⁶, q = 1.6 × 10⁻³), and response to oxidative stress (P = 1.8 × 10⁻⁴, q = 0.0041), highlighting the metabolic and inflammatory features of OA. KEGG pathway analysis revealed significant enrichment of the Leishmaniasis pathway (P = 2.17 × 10⁻⁴, q = 0.021), reflecting immune-related signaling processes relevant to OA inflammation (Fig. 5C). Several additional pathways demonstrated strong nominal enrichment and trend-level significance (0.05 < q < 0.10), suggesting potential biological relevance. Among these, fluid shear stress and atherosclerosis (P = 1.23 × 10⁻³, q = 0.060) and the adipocytokine signaling pathway (P = 5.19 × 10⁻³, q = 0.077) were closely associated with chondrocyte mechanotransduction and metabolic regulation in OA. Pathways related to hypoxia- and stress-response mechanisms, including HIF-1 signaling (P = 1.26 × 10⁻², q = 0.104) and FoxO signaling (P = 1.78 × 10⁻², q = 0.125), also showed suggestive enrichment. Collectively, these pathways highlight the involvement of immune activation, mechanical stress responses, metabolic imbalance, and hypoxia-associated signaling in OA pathogenesis.

GO-based functional similarity analysis revealed distinct patterns of functional relatedness among the 13 WGCNA-derived DE-GMRGs (Fig. 5D). The similarity matrix showed substantial heterogeneity, with several genes clustering tightly together while others remained more isolated. By averaging similarity scores across all pairwise comparisons, IRS2 emerged as the most functionally central gene, displaying the highest cumulative similarity to its peers. In contrast, CHST6 exhibited the lowest similarity, indicating a more functionally distinct role within the gene set. These findings highlight IRS2 as a key functional hub within the shared GO landscape.

Diagnostic model construction and validation

In the GSE57218 discovery cohort, LASSO logistic regression was applied to 13 WGCNA-derived candidate DE-GMRGs. Ten-fold cross-validation identified λ_min = 2.1 × 10⁻⁴ (minimum binomial deviance) and λ_1se = 2.9 × 10⁻³ (most regularized model) (Fig. 6A). At lambda.min, only three genes, F13A1, IRS2, and RELA, retained non-zero coefficients (Fig. 6B) and were selected as the final diagnostic features.

A multivariable logistic regression (LR) model based on F13A1, IRS2, and RELA was then trained in an independent cohort (GSE48556, n = 139). A nomogram was constructed to visualize the model and facilitate individual risk estimation (Fig. 6E). In this nomogram, each gene contributes a point score according to its expression level, and the total points are mapped to a linear predictor and a predicted probability of OA, which ranges from approximately 0.05 to 0.95 across the observed score range. Internal validation in GSE48556 showed good calibration and discrimination. Bootstrap-corrected calibration curves (1,000 resamples) demonstrated close agreement between predicted and observed OA probabilities, with a mean absolute calibration error of 0.028 (Fig. 6C). Decision curve analysis indicated that the three-gene LR model provided a higher net benefit than any single gene, as well as the “treat-all” and “treat-none” strategies, across almost the entire range of threshold probabilities (0–1) (Fig. 6D).

In terms of discrimination, the model achieved an AUC of 0.777 (95% CI [0.691–0.863]) in the training set (Fig. S2A). To assess the stability of model coefficients, we performed 1,000 bootstrap resamples. The bootstrap-based forest plot showed that the odds ratios and 95% confidence intervals for F13A1, IRS2, and RELA consistently deviated from the null line, indicating relatively stable contributions of these three genes to the model (Fig. S2B).

External validation across multiple independent cohorts confirmed the generalizability of the model. In GSE55457, GSE114007, GSE51588, and GSE29746, the LR model yielded AUCs of 0.810 (95% CI [0.612–1.000]), 0.925 (95% CI [0.836–1.000]), 0.848 (95% CI [0.738–0.957]), and 0.628 (95% CI [0.381–0.875]), respectively (Fig. 7A). A performance heatmap further indicated that, at a fixed decision threshold, the model maintained generally high sensitivity, specificity, accuracy, positive predictive value (PPV), and negative predictive value (NPV) in each dataset (Fig. 7B). A forest plot of corrected AUCs showed a summary AUC of 0.827 across these validation cohorts (Fig. 7C). Given the varying case–control ratios across validation cohorts, we additionally evaluated precision–recall performance. Faceted precision–recall curves showed that the areas under the precision–recall curve (AUPRC) in GSE51588, GSE55457, GSE114007, and GSE29746 were 0.944, 0.965, 0.650, and 0.856, respectively, all exceeding the baseline expected under random classification (Fig. S2C). Detailed classification performance at the pre-specified decision threshold, including sensitivity, specificity, accuracy, PPV, and NPV for each cohort, is summarized in Table S5.

Overall, these results indicate that the three-gene glutamine metabolism-based LR model has favorable discrimination, calibration, and clinical net benefit, and shows robust diagnostic performance across multiple independent OA cohorts.

Mechanistic insights into diagnostic genes

Single gene GSEA enrichment analysis results revealed that F13A1 was mainly enriched in extracellular matrix-related pathways, such as Focal Adhesion (NES = −1.73, FDR = 1.21 × 10⁻³), and showed strong negative enrichment in the ECM-receptor interaction pathway (NES = −2.20, FDR = 2.56 × 10⁻⁵) (Fig. 8A). IRS2 was associated with RNA metabolic processes, as evidenced by its positive enrichment in the Spliceosome pathway (NES = 1.75, FDR = 3.96 × 10⁻³) (Fig. 8B). RELA was negatively enriched in mitochondrial energy metabolism, particularly in the oxidative phosphorylation pathway (NES = −2.68, FDR = 9.20 × 10⁻⁹), and was also negatively associated with pathways like Parkinson’s disease (NES = −2.54, FDR = 9.20 × 10⁻⁹) (Fig. 8C).

Analysis of DNA methylation data from the GSE73626 dataset identified 18,297 differentially methylated positions (DMPs), including 9,147 hypermethylated and 9,150 hypomethylated sites. Among the diagnostic genes, F13A1 harbored two DMPs: cg03440267 (hypomethylated, Δβ = −0.57, P = 0.004) and cg16530177 (hypermethylated, Δβ = 0.47, P = 0.027); IRS2 contained one DMP: cg12085119 (hypermethylated, Δβ = 0.23, P = 0.044) (Fig. 9A).

To explore the regulatory mechanisms of F13A1, IRS2, and RELA, we predicted miRNAs targeting these three key genes through integrated analysis. For F13A1, 46 targeting miRNAs were predicted, among which hsa-miR-4731-5p showed significant regulatory interaction with notable expression correlation; for IRS2, 103 potential targeting miRNAs were predicted, displaying complex regulatory patterns with varied expression profiles; for RELA, 128 miRNAs were predicted, including the key hsa-miR-1908-5p, which demonstrated significant impact on RELA expression regulation. Based on miRNA-mRNA interaction predictions, a ceRNA regulatory network containing three mRNAs, 42 key miRNAs, and 158 regulatory nodes was constructed (Fig. 9B). This network comprised 203 nodes and 417 regulatory edges, revealing potential non-coding RNA regulatory mechanisms.

Furthermore, through cross-validation using CTD and DGIdb databases, 14 potential drug intervention targets were predicted based on the three diagnostic genes (Fig. S3). These included several major drug categories: anti-inflammatory pathway drugs (such as troglitazone and dexamethasone targeting RELA), metabolic modulators (such as metformin targeting IRS2), and extracellular matrix regulators (such as thrombin targeting F13A1).

Clinical validation and metabolomic profiling of diagnostic genes

Expression validation in public datasets confirmed that F13A1, IRS2, and RELA were consistently downregulated in OA samples across the GSE57218, GSE55235, and GSE51588 datasets (all P < 0.01) (Figs. S4A, S4B and S4C). A prospective cohort of 31 OA patients and 31 age- and sex-matched healthy controls was then recruited (Table 1). Baseline characteristics (age, sex, BMI) were comparable between groups, whereas OA patients had significantly elevated inflammatory markers (ESR and CRP; both P < 0.01). RT-qPCR of PBMCs confirmed significant downregulation of the three genes in OA patients (all P < 0.01; Fig. 10A). Notably, IRS2 expression in OA patients was negatively correlated with VAS pain scores (ρ = −0.575, P < 0.001) and WOMAC functional scores (ρ = −0.603, P < 0.001) (Figs. 10B and 10C).

LC-MS/MS-based metabolomics of plasma samples identified multiple metabolites with significantly altered levels in OA patients (Fig. 10D). Levels of 1-methylguanosine (−75.5%, P < 0.001), hypoxanthine (−48.1%, P < 0.001), L-aspartic acid (−35.8%, P < 0.001), acetylcarnitine (−34.8%, P < 0.001), xanthine (−23.5%, P = 0.003), and ergothioneine (−21.2%, P = 0.012) were decreased, whereas pantothenic acid (+30.9%, P < 0.001), creatine (+27.8%, P = 0.006), and proline-hydroxyproline (+69.9%, P = 0.021), as well as guanosine (+77.3%, P = 0.032) and cytidine (+29.2%, P = 0.041) were increased.

To validate the mechanistic link between the identified key genes and glutamine metabolism in a clinical context, we analyzed the correlation between gene expression (RT-qPCR) and metabolite abundance (LC-MS) in the validation cohort. This analysis revealed distinct, disease-specific alterations in the gene-metabolite regulatory network. First, RELA expression exhibited a strictly disease-specific association with 1-methylguanosine. While no correlation was observed in healthy controls (Spearman’s ρ = 0.08, P > 0.05), a significant positive correlation emerged exclusively in OA patients (Spearman’s ρ = 0.472, P = 0.007; Fig. 9C), suggesting that RELA-mediated inflammatory signaling may drive specific metabolic stress responses under pathological conditions. Second, we observed a metabolic rewiring involving IRS2 and pyridoxal. In healthy controls, IRS2 expression showed a positive trend with pyridoxal levels. Strikingly, this relationship was inverted in the OA cohort, which exhibited a significant negative correlation (Spearman’s ρ = −0.474, P = 0.007; Fig. 9C). An interaction analysis confirmed that the regulatory link between IRS2 and pyridoxal was significantly modified by the disease state (P < 0.05), suggesting a potentially disrupted metabolic axis in OA pathology.

Beyond these regulatory alterations, multivariable linear regression adjusting for age, sex, and BMI identified pyridoxal (standardized β = 0.29, 95% CI [0.06–0.51]; P = 0.014) and plasma creatine concentration (standardized β = 0.27, 95% CI [0.05–0.50]; P = 0.019) as independent predictors of OA severity. In contrast, no significant association was observed for cytidine (standardized β = 0.21, 95% CI [−0.03 to 0.44]; P = 0.081; Table S6).