In this study, we used GEO database (https://www.ncbi.nlm.nih.gov/geo/) to find the genes related to sarcopenia by differentially expressed genes (DEGs) and WGCNA analysis. Then we further screened out the core genes by machine learning and PPI network analysis. SHAP algorithm helps us to explain the diagnostic model of sarcopenia. At the same time, we also combine the GWAS data about sarcopenia with the eQTL data of specific genes to find out the possible causal relationship between key genes and sarcopenia. We also established a model of sarcopenia with C2C12 cells, so that we can verify the key genes through quantitative polymerase chain reaction (qPCR) experiments in vitro.The overall study design and analytical workflow are summarized in Figure 1, which was constructed as a graphical abstract to visually illustrate the integrated research framework, following the principles of effective graphical abstract design (Kamble et al., 2025).

We got the data set about sarcopenia from GEO database. There are four data sets that can be used publicly: GSE1428 (Giresi et al., 2005), GSE8479 (Melov et al., 2007), GSE136344 (Gueugneau et al., 2021) and GSE111016 (Migliavacca et al., 2019). In this study, we selected GSE1428 and GSE8479 as training data sets, which included 37 sarcopenia samples (12 from GSE1428 and 25 from GSE8479) and 36 control muscle samples (10 from GSE1428 and 26 from GSE8479). In order to deal with the batch effect and ensure the consistency of data, we first use the “SVA” toolkit of R language to correct the batch effect. After that, we combined the arrays of each group and used principal component analysis (PCA) to reduce the dimension and the possible batch effect. Then, we use R’s “Limma” toolkit to analyze the differential expression of the integrated data set to find those DEGs that meet the conditions of |logFC| > 0.5 and P.adjust <0.05. Finally, we use the “Pheatmap” toolkit to draw these DEGs in the form of heat maps. Because all the data sets we use are public, this study does not need to be approved by the ethics Committee.

WGCNA constructed a scale-free network by linking the gene expression level with clinical features (Langfelder and Horvath, 2008). In this study, we used the “WGCNA” of R package to find the gene module which is particularly closely related to sarcopenia. First of all, we choose the most suitable soft threshold to build a network structure that conforms to the scale-free distribution. Then, we build a weighted adjacency matrix, and then turn it into a topological overlap matrix (TOM). After that, we calculated the difference between genes according to TOM, and grouped the genes by dynamic tree cutting algorithm, so that we could find out all kinds of co-expression modules. Finally, we chose the gene module which is most related to the phenotype of sarcopenia for further analysis.

GO enrichment analysis is a commonly used bioinformatics method, which divides the functional description of genes into three different types: molecular function, biological process and cellular component (Aleksander et al., 2023). Similarly, KEGG pathway enrichment analysis is often used to study the biological mechanism and function of genes (Kanehisa and Goto, 2000). These descriptions cover many fields, such as human diseases, biological processes, and signal pathways, which can help us find out the possible connections between genes related to sarcopenia and let us better understand how they lead to the development of this disease. Based on this, we conducted a comprehensive GO functional analysis and KEGG pathway enrichment analysis on the common genes, using the two toolkits of “clusterProfiler” and “DOSE” in R language. We set an important standard for these analyses, that is, the p value should be less than 0.05 to be meaningful. Finally, we use the “GOplot” toolkit to draw a circular diagram to show these results.

In the R language environment, we use 12 different algorithms for variable selection and model building, including ElasticNet, Stepglm, Lasso, Generalized Linear Model Boosting (glmBoost), Ridge, Elastic Net (Enet), Support Vector Machine (SVM), Partial Least Squares Regression for Generalized Linear Models (plsRglm), Linear Discriminant Analysis (LDA), Random Forest (RF), NaiveBayes and eXtreme Gradient Boosting (XGBoost) (Liang et al., 2025). This method combines linear regression model with advanced nonlinear learning method, which can give us a lot of analytical inspiration. The whole calculation framework adopts double algorithm strategy: (1) Firstly, the prediction factors are selected by feature sorting through recursive feature elimination, and then the prediction model is established by stacking generalization technology. (2) This method runs under the framework of hierarchical 10-fold cross-validation, and 113 different configurations are generated by system optimization superparameter (Li T. et al., 2025). Then, we calculated the consistency index (C-Index) of all combinations with the external verification set, and visualized the results through the heat map for evaluation (Zhang W. et al., 2025). Based on GSE1428 and GSE8479 data sets, we use the combination of 113 algorithms mentioned above to build a diagnosis model. In addition, the model is externally verified on two different data sets, GSE136344 and GSE111016.Finally, the model achieving the highest average area under the curve (AUC) across both the training and test sets was selected as the optimal model for this study.

A PPI network was established utilizing information obtained from the STRING database to depict the interconnections among various proteins (Szklarczyk et al., 2019). This network specifically highlights the overlap between DEGs and the genes identified through the WGCNA analysis. Using the Cytoscape softwareplugin “Cytohubba” (Chin et al., 2014) the top 15 genes were selected based on their scores in various network centrality measures, including MCC, EPC, Degree, MNC, EcCentricity, DMNC, Radiality, Bottle Neck, Closeness, Betweenness, Clustering Coefficient, and Stress.The core genes of Cytohubba were then identified from these genes.

In order to assess the diagnostic potential and expression levels of critical genes associated with sarcopenia, we employed the R package “pROC” to generate receiver operating characteristic (ROC) curves. This method helps us to calculate the AUC, and we use this AUC to evaluate the classification ability of each gene in distinguishing the people with sarcopenia from the control group. In the training and verification data set, as long as the AUC value of genes exceeds 0.7, we think that they have a good diagnostic effect, indicating that they may become biomarkers of sarcopenia. In addition, the R package “PerformanceAnalytics” was applied to analyze and visualize the correlations among the key genes in the validation dataset.

To improve the transparency and interpretability of the model, the SHAP algorithm was applied to explain the sarcopenia diagnostic model. This method will score each feature and calculate the SHAP value of each feature, so that we can understand how each feature affects the prediction results of the model (Zhang M. et al., 2025).

The “CIBERSORT” package was employed to evaluate the degree of immune cell infiltration within the gene expression dataset associated with sarcopenia.In order to show the number and distribution of immune cells in different samples, we draw a histogram with the tool of “ggplot2”. Then, we use the R toolkit ggpubr to compare the proportion of 22 different immune cells in sarcopenia samples and normal control samples, and use Wilcoxon test for statistical analysis. The analysis results are displayed by the stack histogram made by ggplot2. In addition, we also use the “corrplot” toolkit to draw the correlation between 22 kinds of immune cells. When the p value is less than 0.05, we think that the correlation is significant.

MR is conducted to verify whether the direction of effects observed in the differential gene expression model aligns with the results from the MR approach.The exposure data and gene eQTL data we used were all obtained from GWAS summary data. As for the results related to sarcopenia-such as lean body weight of limbs, grip strength of hands, and walking speed-they are all from UK Biobank, and only use the information of European participants (Supplementary Table S1). In MR analysis, we regard the gene expression level related to sarcopenia as the exposure factor, while the gene variation, especially single nucleotide polymorphisms (SNPs) in eQTL data set, is used as a tool variable to study their possible causal relationship with the risk of sarcopenia. In order to ensure the reliability of the tool, we only consider those genes with at least three effective SNPs. We will evaluate the strength of these tool variables by calculating the F statistics of each SNP set. If the F statistic is greater than 10, it means that it is a strong tool variable, which can help reduce the possible bias in MR analysis.

When we do MR analysis, we mainly use the method of Inverse Variance Weighted (IVW). At the same time, we also use other methods to help, such as MR-Egger, weighted median, simple mode and weighted mode. In order to ensure the reliability of the results, we will carry out pleiotropic test, heterogeneity evaluation and sensitivity analysis. The method of sensitivity analysis is to remove one data at a time and then see if the results will change. The screening criteria we set are: the P value of IVW should be less than 0.05, the odds ratios (OR) obtained by different MR methods should be consistent, and the P value of pleiotropic test should be greater than 0.05. In order to improve the accuracy and reliability of MR analysis, we will use MR-PRESSO to find out the biased SNPs and then remove them. We evaluate the heterogeneity by IVW and MR-Egger test. If the p value is greater than 0.05, it means that there is no obvious heterogeneity in the data. In addition, we use MR-Egger intercept test to check pleiotropy. If the p value is greater than 0.05, it means that there is no pleiotropic effect. For those genes that meet all these criteria, we will use various graphs for further analysis and display, such as forest plots, scatter plots, funnel plots and leave-one-out analysis, so that we can see the causal relationship more clearly.

Mouse C2C12 myoblasts were originally obtained from the Cell Bank of Chinese Academy of Sciences (CBCC) and procured for this study from Jinyuan Biotechnology (Shanghai, China).Cells were divided into Control and Model groups. To establish an in vitro sarcopenia model, cells in the Model group were treated with 100 μM hydrogen peroxide (H₂O₂) for 24 h, while the Control group was cultured under identical conditions without H₂O₂ (Feng et al., 2022; Wu X. et al., 2025). After treatment, total RNA was extracted from the C2C12 cells using the standard TRIzol method. The purity and concentration of extracted RNA were measured by NanoDrop spectrophotometer. According to the method of reverse transcription kit, we reverse transcription a microgram of RNA and made a cDNA.qPCR was performed on an ABI 7500 Real-Time PCR System using SYBR Green Master Mix. The expression levels of the target genes were normalized to β-Actin as the internal reference gene. The primer sequences used in this study are listed in Table 1.

We use R software (version 4.5.1) and GraphPad Prism for statistics. Continuous variables are presented as mean ± standard deviation, while categorical variables are reported as frequencies and percentages. Using the method of T-test, we compared the mRNA expression levels between different groups.To evaluate the diagnostic accuracy of the candidate mRNA, a ROC curve was plotted, and the AUC was determined through the application of a logistic regression model. Additionally, sensitivity and specificity at the optimal cutoff value were determined. The ROC curve was constructed by plotting sensitivity against (100-specificity), with the results expressed as the AUC and its 95% confidence interval. All statistical evaluations were performed utilizing the R programming language.

The GSE1428 and GSE8479 microarray datasets were subjected to normalization and batch correction, and subsequently merged to form a larger cohort for both training and internal validation purposes. After the normalization and batch correction processes, PCA was performed, demonstrating successful integration of the samples from both datasets, which confirmed effective harmonization across platforms (Figures 2A,B). Using the limma package with criteria of FDR <0.05 and |Log2 FC| > 0.5, we identified a total of 318 DEGs (Supplementary Table S2) when comparing the sarcopenia patient cohort to the healthy control group (Figures 2C,D).

In order to systematically identify key genes associated with the sarcopenia phenotype,we undertook an assessment of the weighted gene co-expression network. As illustrated in Figure 3, the MEbrown module exhibited the most pronounced negative correlation with sarcopenia within the studied cohort (correlation = −0.32, P = 0.006, refer to Figures 3A,B). Additionally, we observed a statistically significant variation in the distribution of gene significance (GS) across various modules (P = 0.05). Subsequently, we conducted an intra-module analysis, which uncovered a significant positive correlation between module membership (MM) and GS specifically in the MEbrown module (correlation = 0.52, P = 3.7e−27) (Figure 3C). Finally, we intersected the 372 genes in the MEbrown module with the previously identified DEGs, resulting in the selection of 109 candidate genes for further machine learning analysis (Figure 3D).

GO and KEGG enrichment analyses were performed on the genes identified as common across the datasets. The GO enrichment results showed that these genes were mainly involved in gland morphogenesis, retinoic acid metabolic process, complement activation, and humoral immune response in the Biological Process (BP) category, suggesting their potential roles in tissue development, immune regulation, and retinoid metabolism (Figure 4A). In the Cellular Component (CC) category, the significantly enriched terms included collagen-containing extracellular matrix, contractile muscle fiber, actin filament bundle, and myosin filament, indicating that DEGs are closely associated with cytoskeletal organization and extracellular matrix remodeling. In the Molecular Function (MF) category, extracellular matrix structural constituent, actin binding, oxidoreductase activity, and NAD/NAD⁺ binding were prominently enriched, implying the involvement of these genes in structural maintenance, cytoskeletal dynamics, and redox-related metabolic processes. Furthermore, KEGG pathway analysis demonstrated that DEGs were significantly enriched in pathways such as cytoskeleton in muscle cells, transcriptional misregulation in cancer, complement and coagulation cascades, phagosome, and retinol metabolism (Figure 4B). These pathways suggest that the dysregulated genes may participate in cytoskeletal reorganization, immune and inflammatory responses, metabolic reprogramming, collectively highlighting the potential roles of DEGs in extracellular matrix remodeling, immune modulation, and metabolic alterations underlying the studied condition.

For the 109 potential key genes identified earlier, leading to the generation of 113 distinct combinations utilizing 12 varied machine learning algorithms aimed at determining the most dependable diagnostic model. To ensure the robustness of our findings, we extended this comprehensive analysis beyond the initial training model to include two independent external validation datasets (Figure 5A). We assessed the utility of each feature selection technique by strictly monitoring the classification accuracy of the resulting models on the validation set. For the internal training dataset, we applied a 10-fold cross-validation strategy across every algorithmic combination to compute the associated AUC values. As detailed in the rankings presented in Figure 5A, the pairing of RF with Ridge regression emerged as the clear standout. This combination demonstrated superior predictive capability across both the internal and external cohorts, achieving AUC values of 1.000 (95% CI: 1.000–1.000), 0.852 (95% CI: 0.720–0.953), and 1.000 (95% CI: 1.000–1.000), respectively (Figure 5B). The confusion matrix shown in Figure 5C further indicates that the model achieved an impressive overall accuracy of 98.6% on the training set, calculated as (35 + 37)/(35 + 0+1 + 37). The model exhibited strong classification performance across all datasets, demonstrating high predictive accuracy (Figure 5C). The RF + Ridge algorithm identified 25 hub genes (HOXB2, SLPI, TPPP3, PLAG1, COL4A5, PPIF, KLF13, LYVE1, FOXO1, ST8SIA5, TMEM158, KAZALD1, ZBTB16, FST, ASXL1, SIRT4, ACTC1, TAF5, SLC7A6, CYP26B1, CCDC69, GADD45G, C1S, ATP1B4, MGP).

Figure 6A illustrates the distinct expression profiles of 25 key genes in individuals with sarcopenia compared to control participants, as represented in a volcano plot.We used the Cytoscape plugin cytoHubba to perform centrality analysis on the PPI network composed of 25 hub genes (Figure 6B). The intersection of all twelve ranking algorithms integrated by Cytohubba revealed that FOXO1, ZBTB16, HOXB2, LYVE1, MGP, and CYP26B1 were identified as key genes for further analysis (Figure 6C).

Single-gene ROC analysis indicated that the AUCs for FOXO1, HOXB2, ZBTB16, LYVE1, MGP, and CYP26B1 were 0.860, 0.982, 0.843, 0.873, 0.808, and 0.835, respectively, with HOXB2 showing an AUC >0.98, demonstrating excellent discriminative power (Figure 7A). The analysis of gene expression demonstrated a notable upregulation of the FOXO1, ZBTB16, LYVE1, MGP, and CYP26B1 genes within the sarcopenia cohort. Conversely, the HOXB2 gene exhibited a significant downregulation in this same group, thereby underscoring their possible involvement in the pathology of the disease (Figure 7B). The SHAP bee swarm plot and bar plot illustrated the impact of six key genes on model predictions. HOXB2 exhibited the highest average absolute SHAP value (0.213) and showed negative regulation, indicating its dominant role in discrimination, but increased expression was inversely related to predicted risk. LYVE1, ZBTB16, and FOXO1 showed positive SHAP values (0.057–0.066), indicating that higher expression levels of these genes were linked to an elevation in the predicted outcomes. MGP and CYP26B1 had weaker effects (SHAP <0.045) (Figures 7C,D). SHAP dependence plots revealed distinct non-linear and dose-dependent relationships between gene expression levels and their contributions to model output. Specifically, the negative contribution of HOXB2 plateaued at higher expression levels, suggesting a saturation effect (Figure 7Ei), while FOXO1 exhibited a transition from a negative to a positive impact as its expression increased, indicating a threshold-dependent shift in its influence on model prediction (Figure 7iv). Similar non-linear patterns were also observed for LYVE1, ZBTB16, MGP, and CYP26B1 (Figure 7Eii, iii, v, vi), further supporting the robustness of SHAP-based interpretability in revealing complex gene–model interactions.

We employed the CIBERSORT algorithm to deconstruct the immunological landscape, quantifying the relative proportions of 22 infiltrating immune cell types across both sarcopenia and control samples (Figure 8A). Subsequent comparative analysis using the Wilcoxon test highlighted distinct immunological signatures, specifically revealing significant disparities in the abundance of M2 Macrophages and resting Dendritic cells between the two cohorts (Figure 8B). To further elucidate the interplay between molecular drivers and the immune microenvironment, we constructed correlation network heatmaps linking core genes to specific immune subsets (Figure 8C) and analyzed intra-immune interactions (Figure 8D). Notably, genes such as HOXB2, LYVE1, MGP, and CYP26B1 exhibited significant negative associations with resting NK cells, while ZBTB16 revealed a negative association with M1 macrophages (Figure 8E). These patterns strongly imply that macrophage and NK cell dynamics are central to the immune-regulatory pathways modulated by these hub genes.

To probe the causal underpinnings of sarcopenia, we established a two-sample MR framework, mining the IEU Open GWAS database for valid SNPs associated with our six candidate genes (FOXO1, HOXB2, ZBTB16, LYVE1, MGP, and CYP26B1).By applying stringent instrumental variable criteria (P < 5e-8), we successfully isolated robust genetic instruments.The MR analysis demonstrated a noteworthy positive causal relationship between the expression levels of CYP26B1 (Odds Ratio = 1.006, 95% Confidence Interval: 1.002–1.011, P = 0.004) and the likelihood of developing sarcopenia (Figure 9A). While MR-Egger regression and pattern-based analyses yielded consistent directional outcomes, they did not achieve statistical significance. The heterogeneity assessments (Inverse Variance Weighted Q = 5.46, P = 0.14; MR-Egger Q = 1.18, P = 0.552) and pleiotropy evaluations (Egger intercept = 0.007, P = 0.174) revealed no significant heterogeneity or pleiotropic effects. Furthermore, the leave-one-out sensitivity analysis (Residual Sum of Squares observed = 7.01, P = 0.451) substantiated that no individual SNP exerted a considerable impact on the overall results, thereby reinforcing the reliability of these findings (Figures 9B–D).

To validate the reliability of the bioinformatic screening and the causal relationships identified through Mendelian randomization, we established an in vitro model of sarcopenia using C2C12 mouse myoblasts. As illustrated in Figures 10A–F, the mRNA expression levels of the six hub genes were quantified via qPCR. Compared with the Control group, the expression levels of Cyp26b1 (Figure 10A, P < 0.05), Foxo1 (Figure 10B, P < 0.01), Lyve1 (Figure 10D, P < 0.0001), Mgp (Figure 10E, P < 0.01), and Zbtb16 (Figure 10F, P < 0.001) were significantly elevated in the Model group. Conversely, the expression of Hoxb2 (Figure 10C, P < 0.05) exhibited a significant decrease following induction. These experimental results are highly consistent with the mRNA expression profiles identified in the GEO datasets (GSE1428 and GSE8479) and the machine learning predictions, confirming that these genes are key regulatory nodes in the progression of sarcopenia and providing a solid molecular basis for further mechanistic exploration.