Microarray datasets were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) (19), which hosts a wide range of high-throughput sequencing and gene expression data. The terms “systemic sclerosis” and “sarcopenia” were used as search keywords to identify relevant gene expression datasets, and non-human samples were excluded.

The GSE181549 dataset (20), based on the GPL13497 platform, includes 113 skin samples from SSc patients and 44 from healthy controls. The GSE167186 dataset (21), generated using the GPL20301 platform, comprises 72 skin samples, including 24 from patients with sarcopenia and 48 from healthy controls. All patients with SSc met the diagnostic criteria established by the American College of Rheumatology (ACR) (22), and sarcopenia diagnoses were based on the guidelines set by the European Working Group on Sarcopenia in Older People (EWGSOP) (23).

As these datasets are publicly available and freely accessible, no ethics committee approval was required. The complete workflow is illustrated in Figure 1 .

To ensure that the same gene was compared in both datasets, we first converted Agilent probe IDs (GSE181549) and Illumina Ensembl IDs (GSE167186) to official HGNC gene symbols using the most recent platform annotation files. Only genes unambiguously mapped to identical symbols in both datasets were retained. DEGs were identified separately from the GSE181549 and GSE167186 datasets using the “limma” R package. Benjamini and Hochberg’s approach was used to adjust P-values for controlling the false discovery rate (FDR). Genes with a adjusted P-value < 0.05 and an absolute log fold change (|logFC|) > 0.5 were considered statistically significant. The results were visualized through gene clustering heatmaps and volcano plots. A comparative analysis of DEGs from the two datasets was performed, and overlapping genes potentially involved in the pathogenesis of both SSc and sarcopenia were identified using a Venn diagram.

Further screening of the shared genes between SSc and sarcopenia was conducted using three widely accepted machine learning algorithms: Least Absolute Shrinkage and Selection Operator (LASSO), Random Forest (RF), and Support Vector Machine - Recursive Feature Elimination (SVM-RFE).

Initially, the 19 previously identified common genes were input into the LASSO algorithm within the SSc dataset. A regression model was constructed using the “glmnet” R package with 10-fold cross-validation and a seed value of 100. The “family” parameter was set to “binomial,” and the optimal value of the regularization parameter was determined using “lambda.1min.” The log(λ) profiles of the LASSO coefficients for the 19 features were plotted, followed by the generation of the binomial deviance curve and the log(λ) curve. The optimal value for 1 standard error (1-SE) from the minimum criterion was subsequently determined.

Next, the Random Forest algorithm was implemented using the R package randomForest (v4.7-1.1). The model was trained on the SSc cohort with the following parameters: ntree = 500 (number of trees grown), mtry = floor(sqrt(p)) = 4 (number of variables randomly sampled as candidates at each split, where p=19 is the total number of input genes), nodesize = 1 (minimum size of terminal nodes) and the Gini impurity criterion. To ensure reproducibility, a random seed was set to 100. Variable importance was assessed with the Gini impurity decrease (MeanDecreaseGini), and the top 10 predictors were visualized.

Finally, SVM-RFE recursively eliminated the least important features based on their SVM weight coefficients while conducting cross-validation to evaluate model accuracy and error rates (24). This algorithm was performed using the e1071 package (v1.7-14) with default parameters: radial basis function kernel, cost = 1, and gamma = 1/19. The algorithm recursively eliminated features with the smallest absolute weight coefficients. Ten-fold cross-validation monitored prediction accuracy across progressively smaller feature subsets (from 19 to 1 feature), with the subset achieving minimum cross-validation error selected as optimal.

Genes identified by the intersection of all three machine learning algorithms, visualized through a Venn diagram, were selected as the diagnostic gene targets for SSc-associated sarcopenia.

Additionally, a similarity analysis of the common DEGs was performed using the “GOSemSim” R package to compare these genes with reference genomes. The geometric mean of Gene Ontology (GO) semantic similarity scores—encompassing biological process, molecular function, and cellular component—was calculated. The results were illustrated using a cloud and rain plot generated by the “ggplot2” R package.

The expression levels of candidate diagnostic genes were evaluated separately in the SSc and sarcopenia cohorts. The diagnostic utility of these biomarkers was assessed using receiver operating characteristic (ROC) curve analysis, and the area under the curve (AUC) was calculated.

Multivariate logistic regression analysis was conducted to construct a risk score for SSc-associated sarcopenia based on the shared diagnostic genes. For each sample, the risk score was computed using the gene expression levels and the corresponding logistic regression coefficients. The formula for the risk score was defined as follows:

where Coefficient_i represents the regression coefficient for gene i, and Expression_i denotes the expression level of gene i in each patient.

A nomogram was constructed using the “rms” R package to enhance the clinical applicability of the risk model. The performance of the model was evaluated through ROC curve analysis and AUC values. Model calibration and predictive consistency were assessed using calibration curves and concordance indexes (C-index).

To ensure the internal validation of model performance and guard against the risk of overfitting, this study adopted a 5-fold cross-validation strategy in the systemic sclerosis dataset. We repeated the entire 5-fold cross-validation process 200 times. During each iteration of cross-validation, we meticulously calculated key performance metrics, including the Discriminant Index (Dxy), AUC value, Coefficient of Determination (R²), and Brier score. Finally, we aggregated the mean values of all performance metrics to obtain a comprehensive assessment of the model’s performance.

Peripheral blood samples were collected from 32 SSc patients and 20 healthy individuals at Fujian Provincial Hospital. All patients were diagnosed according to the ACR classification criteria for SSc (22). The study protocol was approved by the Ethics Committee of Fujian Provincial Hospital. All procedures complied with the Declaration of Helsinki, and written informed consent was obtained from all participants.

qPCR was performed to validate the expression levels of the hub mRNAs, NOX4 and NEK6. Total RNA was extracted using TRI REAGENT BD (MRCGENE) according to the manufacturer’s instructions. RNA was then reverse transcribed into complementary DNA (cDNA) using the GeneAmp PCR System 9700 (Applied Biosystems). qPCR was conducted on a QuantStudio™ 5 Real-Time PCR System (Applied Biosystems), using GAPDH as the internal control. The relative expression levels of target genes were calculated using the 2^−ΔΔCt method. Primer sequences for the hub genes are listed in Table 1 .

Using PCR data obtained from our center as the validation set, we assessed the diagnostic performance of NOX4, NEK6, and a dual-gene model for SSc by constructing ROC curves. Calibration curves were also generated to provide a comprehensive evaluation of model performance. To further explore the potential involvement of NOX4 and NEK6 in the pathogenesis of SSc-associated sarcopenia, we conducted comparative analyses across three distinct cohorts: 12 SSc patients with sarcopenia, 20 SSc patients without sarcopenia, and 20 healthy controls. This approach allowed for systematic evaluation of both diagnostic efficacy and mechanistic associations. Patients classified as having sarcopenia also met the diagnostic criteria established by the Asian Working Group for Sarcopenia 2014 (AWGS2014) (25).

To investigate immune infiltration patterns in the GSE181549 dataset, we utilized the CIBERSORT algorithm to perform deconvolution analysis, quantifying the relative proportions of 22 immune cell subpopulations (26). The distribution of these immune subsets was visualized using a stacked bar chart generated with the “ggplot2 (v3.3.0)” package. Comparative analysis between SSc patients and healthy controls was then performed. Spearman’s rank correlation coefficients were computed to evaluate associations between candidate diagnostic genes and the abundance of immune cells. We also have applied the Benjamini-Hochberg procedure to adjust the p-values for controlling the FDR across the multiple comparisons. Hierarchical clustering heatmaps were constructed using the “pheatmap” package to visualize the relationships between the two conserved genes and immune cell composition. The adjusted P-values < 0.05 were considered statistically significant.

All statistical analyses and visualizations were performed using R software (version 4.2.1). Differences in gene expression levels and immune cell fractions between clinical groups were assessed using a two-sided Wilcoxon test. Correlation analyses were conducted using Spearman’s rank test. A P-value < 0.05 was considered indicative of statistical significance.

Analysis of the SSc dataset GSE181549 revealed 382 DEGs in the SSc group ( Figure 2A ). To identify DEGs associated with sarcopenia, skeletal muscle samples from 24 sarcopenia patients and 48 healthy controls were analyzed using dataset GSE167186, yielding 508 DEGs ( Figure 2B ). Intersecting the DEGs from both datasets identified 19 shared genes ( Figure 2C ), providing valuable insights into the common molecular mechanisms underlying SSc and sarcopenia. Heatmaps illustrating these 19 intersecting genes in both datasets ( Figures 2D, E ) demonstrated clear differences in gene expression patterns between disease and control groups. The 19 shared genes were: FCN1, NOX4, STC2, NEK6, LTBP2, H19, HAPLN3, SLC16A3, IGFBP3, IGSF10, VCAM1, PCOLCE2, NNMT, EMX2, CHST1, LPL, RPRML, ADIPOQ, ZNF469.

To further identify candidate diagnostic biomarkers capable of discriminating between disease and control states, three machine learning algorithms were applied to the 19 shared genes: LASSO regression, RF, and SVM-RFE. In the SSc group, LASSO coefficient profiling ( Figure 3A ) and optimal lambda selection ( Figure 3B ) indicated that the tuning parameter λ was 0.0146, which balanced model complexity and predictive performance. The LASSO model identified 10 candidate genes. These 19 genes were also evaluated using the RF classifier, and the top 10 genes were ranked by importance ( Figure 3C ). To further improve diagnostic accuracy, the SVM-RFE method was employed. Model accuracy peaked at 0.879 with 13 genes ( Figure 3D ), and the minimum error rate was 0.121 ( Figure 3E ). Ultimately, five genes—NOX4, STC2, NEK6, IGSF10, and EMX2—were identified through the intersection of all three algorithms ( Figure 3F ). Functional prioritization by the “Friends” analysis confirmed that NOX4, STC2, NEK6, and EMX2 were among the top 10 functionally relevant genes ( Figure 3G ).

For the five candidate biomarkers, single-gene expression analyses were conducted. In the GSE181549 dataset, the expression levels of NOX4, STC2, and NEK6 were significantly upregulated in the SSc group relative to controls, whereas IGSF10 and EMX2 were downregulated ( Figure 4A ). In the sarcopenia dataset GSE167186, NOX4 expression was elevated, and NEK6 expression was reduced in sarcopenia patients compared with controls ( Figure 4B ). The expression of STC2, IGSF10, and EMX2 showed no significant differences.

The objective of this study was to identify biomarkers with relevance to both SSc and sarcopenia. Since NOX4 and NEK6 were differentially expressed in both datasets, they were selected for further validation. ROC curve analysis was used to evaluate their diagnostic performance. In the SSc dataset, the AUC was 0.893 (95% CI: 0.845 – 0.941) for NOX4 and 0.795 (95% CI: 0.725 – 0.865) for NEK6 ( Figure 4C ), indicating robust diagnostic potential. In the sarcopenia dataset, the AUC for NOX4 was 0.662 (95% CI: 0.524 – 0.801), and for NEK6, it was 0.661 (95% CI: 0.533 – 0.789) ( Figure 4D ), suggesting limited individual predictive efficacy in sarcopenia.

To identify more effective diagnostic indicators, this study first confirmed that NOX4 and NEK6 serve as independent risk factors for both SSc and sarcopenia through multivariate regression analysis. A risk score model was subsequently constructed by weighting the normalized expression levels of NOX4 and NEK6 using the regression coefficients derived from multivariate logistic regression analysis. As shown in Figure 4E , the SSc risk score was calculated as: normalized expression of NOX4 × 1.624 + normalized expression of NEK6 × 0.704. For sarcopenia ( Figure 4F ), the risk score was: normalized expression of NOX4 × 2.229 + normalized expression of NEK6 × (–2.369). The ROC curves indicated high predictive accuracy for both diseases ( Figures 4G, I ). The C-index of the two-gene risk model exceeded that of the single-gene models, indicating superior predictive performance. The bias-corrected calibration curves were closely aligned with the ideal curve ( Figures 4H, J ). Additionally, to mitigate the risk of overfitting, we conducted internal cross-validation within the GSE181549 dataset. The results, as shown in Supplementary Table 1 , demonstrated an AUC of 0.890, which further confirms the model’s robust classification performance and stability.

To further validate the expression levels of NOX4 and NEK6, whole blood samples were collected from 20 healthy controls and 32 SSc patients, including 12 with sarcopenia and 20 without. The results revealed that NOX4 expression was elevated in both the SSc sarcopenia and non-sarcopenia groups compared to healthy controls ( Figure 5A ). In contrast, NEK6 expression was reduced in both patient groups ( Figure 5B ). However, no statistically significant differences in NOX4 and NEK6 expression were observed between the SSc sarcopenia and non-sarcopenia groups.

Using PCR data from our cohort as the validation set, Figures 5C, D demonstrated satisfactory diagnostic performance of NOX4 (AUC = 0.798, 95% CI: 0.679 – 0.917) and NEK6 (AUC = 0.812, 95% CI: 0.690 – 0.935) for SSc. The two-gene risk score showed superior predictive capability, achieving an AUC of 0.894 ( Figure 5E ). Moreover, the bias-corrected calibration plot closely approximated the ideal curve ( Figure 5F ). To further assess the roles of NOX4 and NEK6 in the pathogenesis of SSc-associated sarcopenia, ROC curves were generated to evaluate their diagnostic performance. NOX4 (AUC = 0.798, 95% CI: 0.633 – 0.962) and NEK6 (AUC = 0.873, 95% CI: 0.749 – 0.996) exhibited high diagnostic accuracy for SSc-associated sarcopenia ( Figures 5G, H ). The two-gene model showed the highest discriminative performance, with an AUC of 0.941 (95% CI: 0.861 – 1.000), supporting its favorable predictive efficacy ( Figure 5I ). The calibration plot demonstrated excellent agreement between predicted probabilities and observed outcomes after bias correction, closely aligning with the ideal curve across the full probability range ( Figure 5J ).

To assess the diagnostic utility of NOX4 and NEK for sarcopenia in SSc, we performed a Pearson correlation analysis with ASMI. ASMI and NOX4 have a substantial negative connection (r=-0.31, p=0.041). Additionally, there is a positive connection (r = 0.45, p = 0.01) between ASMI and NEK6.

A comprehensive analysis was conducted to characterize the immune cell infiltration patterns and their functional relevance in SSc and control groups. Figure 6A illustrates the relative proportions of immune cell subsets, indicating significant alterations in the SSc group. The CIBERSORT algorithm was employed to calculate the proportions of 22 immune cell types in each sample from the GSE181549 dataset. Notably, several immune cells, including activated CD4⁺ memory T cells, activated NK cells, M1 macrophages, and both resting and activated dendritic cells and mast cells, exhibited significantly increased infiltration in the SSc group compared to controls ( Figure 6B ).

As shown in Figure 6C , NOX4 and NEK6 were significantly correlated with infiltration levels of various immune cells in the SSc cohort, including plasma cells, dendritic cells, naive CD4⁺ T cells, macrophages, NK cells, and neutrophils. Specifically, NOX4 displayed strong positive correlations with both M1 and M2 macrophages, and significant negative correlations with resting dendritic cells, activated dendritic cells, plasma cells, activated CD4⁺ memory T cells, and resting mast cells ( Figure 7A ). Correlation coefficients greater than 0.3 were observed between NOX4 and M1/M2 macrophages, plasma cells, and both resting and activated dendritic cells ( Figure 7B ). Similarly, NEK6 showed positive correlations with macrophages and monocytes, and negative correlations with plasma cells, resting and activated dendritic cells, and mast cells ( Figure 8A ). Significant correlations (R > 0.3) were noted between NEK6 and M1 macrophages as well as plasma cells ( Figure 8B ). These findings suggest that NOX4 and NEK6 may influence the pathophysiological progression of SSc by modulating immune cell infiltration, thereby offering novel perspectives for immunomodulatory therapeutic strategies.