Machine learning screening for disulfidptosis genes-associated immunosuppression status in osteosarcoma and rhabdomyosarcoma

Introduction

Sarcomas originate from cancerous transformations in mesenchymal tissues. In children, the two prevalent types of sarcoma are osteosarcoma (OS), which develops in bones, and soft tissue sarcoma, which arises in soft tissues, with rhabdomyosarcoma (RMS) being the most common, accounting for 50% of all soft tissue sarcomas (1). Distinguishing between OS and RMS in the limbs of children can be challenging. The poor prognosis associated with these sarcomas is often attributed to their resistance to chemotherapy and radiotherapy (2,3). To overcome this challenge, immunotherapy is being explored as a novel therapeutic approach (4). Recent advancements in understanding the tumor immune microenvironment and immune checkpoints have led to the increasing use of anti-programmed cell death protein-1 (PD-1) therapy for pediatric solid tumors. Although not yet widely adopted for OS and RMS, studies have shown the potential of this approach (5-7).

Cell death involves irreversible changes such as metabolic arrest, structural damage, and loss of function. Understanding its mechanisms enhances our knowledge of cellular homeostasis and provides insights into treating diseases like cancer. Certain cancer cells resistant to conventional therapies are particularly susceptible to induced ferroptosis (8). Cell death can also trigger changes in the immune response; for instance, immunogenic death, a specific form of cell death, is induced by conventional therapies (9). Recently, a new cell death mechanism, disulfidptosis, has been identified. This mechanism involves glucose starvation-induced disulfide stress and rapid cell death due to massive NADPH depletion and abnormal disulfide accumulation, such as cystine, in SLC7A11-overexpressing cancer cells (10). Previous research has shown that inhibiting NADPH oxidase 2 induces apoptosis in OS cells (11). However, the relationship between this new death mode and the immune microenvironment of OS and RMS remains unclear.

In this study, we integrated bulk RNA-sequencing (RNA-seq) data from public databases for OS and RMS. We first analyzed the immune microenvironment of these diseases using non-negative matrix factorization (NMF) based on the bisulfide death-associated gene set. This analysis identified and validated an immune-response-insensitive immunosuppressive subtype. We then screened the five most representative genes in the disulfidptosis gene collection using 100 combinations of machine learning algorithms. Prognostic models were developed for OS and RMS, revealing that high expression of IQGTPase-activating protein 1 (IQGAP1) indicates a poor prognosis in children with RMS and could serve as a new therapeutic target and independent prognostic marker. Our findings suggest that macrophages play a crucial role in the development and progression of OS and RMS and may be strongly associated with disulfidptosis, impacting patient prognosis. We hope that this study would encourage further research into the relationship between disulfidptosis and these diseases, advancing immunotherapy efforts. We present this article in accordance with the STREGA and MDAR reporting checklists (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1439/rc).

Methods

Data sources used for analysis

We obtained bulk RNA-seq data and clinical information for OS from Therapeutically Applicable Research to Generate Effective Treatments (TARGET) (n=88) and Gene Expression Omnibus (GEO) (GSE21257, n=53, GPL10295), and single-cell RNA-seq data for OS from GEO (GSE162454, n=6, GPL24676). Also, we obtained from TARGET (n=70) and GEO (GSE66533, n=58, GPL570) bulk RNA-seq data and clinical information. All data are in fragments per kilobase of transcript per million mapped fragments (FPKM) format. We used K-Nearest Neighbors (KNN) completion to process missing values and Z-score to normalize the data. We then obtained the data and clinical information of RMS from the preprint article Consensus Clustering Analysis for Disulfidptosis-related Genes, and descriptive pan-cancer genetic analysis of disulfidptosis-related gene set collected 24 disulfidptosis-related genes (DGs). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Analysis of genes of disulfidptosis

In order to gain a deeper understanding of DGs, we used the “RCircos” in the R package to visualise the chromosomal locations of these genes in conjunction with the TARGET gene expression matrices of OS and RMS and the interactions of the corresponding proteins of these genes were explored using the STRING website (https://string-db.org/) and Cytoscape software. Afterwards, to further distinguish the different subtypes associated with disulfidptosis, we used the NMF algorithm in the R software to decompose the matrices

Analysis of subgroups associated with disulfidptosis

In order to explore the immune microenvironment of different subtypes of disulfidptosis, we used the ESTIMATE and CIBERSOFT algorithms to calculate the scores of the TARGET expression matrices of OS and RMS, respectively, to obtain the corresponding immune scores of the matrices and the infiltration of immune cells, and then classified them into subtypes according to the scores to perform gene set enrichment analysis (GSEA). Gene set and phenotypic marker expression files from the Kyoto Encyclopedia of Genes and Genomes (KEGG) were loaded into the GSEA software and run 1,000 times to ensure functional consistency. For data from the GEO databases (GSE21257, GSE66533), they were processed in the same way. In addition, we identified the differential expression of inhibitory receptors in different subtypes.

Five regression algorithms

In order to screen the most representative differential genes for the two subtypes, we selected the TARGET database gene expression matrix, which is more comprehensive in terms of sequenced genes, for a more in-depth analysis. For both diseases, we picked differentially expressed genes between the immune subtype and other subtypes according to |log₂ fold change| >2 and P<0.01, and then we took the intersecting genes obtained from the two diseases and generated a network of functional groupings of KEGG pathways using the plug-in ClueGO in Cytoscape software to explore cell signalling pathways; After that, we applied five machine learning algorithms, Lasso/Elasticnet/Ridge/RF/Rfe (all default parameters refer to the “sklean” Python library), to downscale these common differential genes in order to obtain the most representative differential genes.

Twenty classification algorithms

We performed Spearman correlation analysis between the sum of the above representative differential genes and DGs, and finally screened out the most representative disulfidptosis genes. In order to verify the classification of these disulfidptosis genes on immunosuppressive subtypes and other subtypes, we used 20 mainstream machine learning algorithms [Linear Regression, Ridge Regression, RidgeCV, Linear Lasso, Lasso, ElasticNet, BayesianRidge, Logistic Regression, Stochastic Gradient Descent (SGD), Support Vector Machine (SVM), KNN, Naive Bayes, Decision Tree, Bagging, Random Forest, Extra Tree, AdaBoost, GradientBoosting, Voting, Artificial Neural Network (ANN); all default parameters refer to the “sklean” Python library] to calculate the classification accuracy of the two DGs subtypes, which constitutes a combination of 100 machine learning algorithms, and finally, the risk score formulas consisting of the final screened DGs were obtained by using multifactorial Cox regression analyses among the two disease TARGET databases, and the risk scores were calculated using the corresponding receiver operating characteristic (ROC) curve values were calculated using the “timeROC” and “survivalROC” R packages.

SHAP analysis

For model interpretation, we computed SHapley Additive exPlanations (SHAP) values to quantify the contribution of the ACTN4, MYH9, FLNA, MYH10, and IQGAP1 to the model’s predictions. The mean absolute SHAP value was used to rank global feature importance. Analyses were conducted using the Python shappackage.

Cell lines, antibodies, and reagents

The human OS U2OS cells and RMS RD cell lines were obtained from Wuhan Pricella Biotechnology Co., Ltd. All media were supplemented with 10% (v/v) fetal bovine serum (FBS; HyClone, USA), and cells were incubated at 37 ℃ in a humidified atmosphere containing 5% CO₂. The following primary antibodies were used: IQGAP1 (#YA2388), Tubulin (#14555-1-AP) and Tubulin (#YA3401) from MedChemExpress (Shanghai, China).

Western blot (WB) analysis

Cells and tissue samples were lysed using RIPA buffer containing phenylmethanesulfonyl fluoride (PMSF) (Beyotime Biotechnology, Shanghai, China). Protein extracts were separated on 8–12% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) gels and transferred onto polyvinylidene difluoride (PVDF) membranes. After blocking with 5% bovine serum albumin (BSA), membranes were incubated with primary antibodies at 4 ℃ overnight, followed by incubation with horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 h at room temperature. Protein bands were visualized using a VIBER FUSION FX6 EDGE imaging system (Vilber, France).

Analysis in risk score group

In order to verify the accuracy of the risk score formulas, we plotted the Kaplan-Meier (KM) survival curves based on the high- and low-risk groups distinguished by the risk scores, combined with the overall survival time in days and event-free survival time in days of OS and RMS in the TARGET database using the “survminer” and “survival” R packages. “Survminer” and “survival” R packages were used to plot KM survival curves (with a log-rank P value <0.05, indicates a statistically significant difference in survival outcomes). After that, in order to investigate the differences in cell signalling pathways between the high- and low-risk groups, we used the KEGG gene set (c2.cp.kegg.v7.4) to carry out the Gene Set Variance Analysis (GSVA), and finally obtained the GSVA enrichment score matrix.

Collection of published prognostic models of OS

To investigate whether OS and RMS can use the same genes to construct prognostic models, we searched the currently published literature on OS but did RMS due to its less study. So for OS, we used the keywords “osteosarcoma”, “model”, “GSE21257”, “prognostic” and other keywords to count the published prognostic models in PubMed website (https://pubmed.ncbi.nlm.nih.gov/). In order to prevent the computational complexity caused by the large number of genes constituting the models, we only selected model genes with a gene number less than 10 for subsequent analyses. Then, we performed multifactorial Cox regression analyses of the corresponding model genes in the RMS TARGET database, and calculated the ROC curve values using the “survivalROC” R package. The results were compared with our screened DGs model. In addition, in order to explore the correspondence between different subgroups (divide into high-risk group and low-risk group based on the level of risk score) and clinical variables of the two diseases, we also used the “networkD3” R package to draw the Sankey diagrams of the TARGET databases of OS and RMS.

Analysis of independent prognostic genes

Since we only found skeletal muscle data from the Human Protein Atlas (https://www.proteinatlas.org/), we did not find any information on bone data. Therefore, in order to explore the model genes in RMS more deeply, we combined the clinical information corresponding to the RMS TARGET database with the “forestploter” R package to draw a forest plot, search for molecules with independent prognostic effects, and draw KM survival curves, which were then displayed in the INTERACTION module of the website. The INTERACTION module and finally the SINGLE CELL TYPE module were used to search for differences in the expression of this molecule in skeletal muscle in different cells.

Single-cell data analysis

Since we could not find any single-cell data of RMS in the GEO public database, we used “Seurat” and “NormalizeData” R packages to pre-process GSE162454 in our in-depth study of OS. Following quality control of the GSE162454 dataset using the Seurat R package (retaining cells with 200–4,500 unique genes, 1,000–35,000 unique molecular identifiers (UMIs), and a mitochondrial gene percentage below 10%), the data were normalized. Subsequently, the expression of DGs was projected onto the cell populations annotated by the SingleR package. Population annotated in the “SingleR” R package. In order to make a more accurate judgement on the trajectory of cell differentiation and reduce the chance as much as possible, we randomly divided the samples of the database into three groups with the same number of samples, and used Monocle3 for the proposed time series analysis of cell subpopulations, and cellchat and cellphoneDB for the analysis of subpopulation cell communication.

Statistical analysis

All statistical analyses were performed using R software (version 4.1.2) and its associated packages. For data preprocessing, missing values were imputed using the KNN algorithm, and continuous variables were normalized using Z-score transformation. Differential expression analysis between groups was conducted using the Limma package, with a significance threshold set at |log₂ fold change| >2 and a two-sided adjusted P value <0.05, employing the Benjamini-Hochberg method for multiple testing correction. For survival analysis, KM curves were generated and compared using the log-rank test via the “survival” and “survminer” packages. Univariate and multivariate Cox proportional hazards regression models were applied to identify independent prognostic factors, with results visualized using forest plots. The predictive accuracy of the prognostic model was evaluated by time-dependent ROC curves and area under the curve (AUC) values calculated with the “timeROC” package. Machine learning algorithms were implemented using the “caret” and “glmnet” packages, with model performance assessed through 10-fold cross-validation. For all analyses, a two-sided P value <0.05 was considered statistically significant, unless otherwise specified for multiple comparisons.

Results

Immunosuppressive states associated with disulfidptosis

We discovered that disulfidptosis-related gene-encoded proteins exhibit strong interactions with each other (Figure 1A). These genes are also amplified and co-expressed in both OS and RMS (Figure 1B,1C). This evidence suggests that disulfidptosis is a complex death process involving multiple gene regulations and protein interactions. Using the NMF algorithm, we identified three subtypes related to disulfidptosis in both diseases (Figure 1D,1E). Immune scoring of these subtypes revealed that subtype “1” in OS has higher tumor purity compared to the other two subtypes. The Tumor Immune Score (TIS) is a quantitative metric that evaluates the overall infiltration and functional activity of immune cells—such as T lymphocytes, B lymphocytes, and natural killer (NK) cells—within tumor tissues. An elevated TIS generally correlates with a robust anti-tumor immune response, while a reduced TIS indicates impaired immune cell infiltration or functional suppression. The Tumor Microenvironment Score (TME) provides a more comprehensive assessment of the tumor microenvironment composition, incorporating immune cells, stromal cells, angiogenic components, and extracellular matrix elements. A diminished TME score reflects a lower abundance of immune constituents, suggesting potential dominance by tumor cells or immunosuppressive stromal elements. Thus, the concurrent observation of low TIS and TME values provides a rationale for inferring the presence of an immunosuppressive subtype. With relatively lower TIS and TME scores, we defined this as the immunosuppressive subtype (Figure 1F,1G). In RMS, the immunosuppressive subtype comprises subtypes “2” and “3”. These findings indicate that the immunosuppressive subtype is not highly sensitive to immunotherapy (Figure 1H).

Figure 1 Immunosuppressive states associated with disulfidptosis in OS and RMS. (A) The PPI network among the DGs. (B,C) Circus plots of chromosome distributions of DGs from TARGET cohort. (D,E) The heatmap of TARGET samples by unsupervised NMF illustrates 3 distinct expression patterns. (F) Heatmap shows the cluster 1 exhibited low score of TumorPurity, Stromalscore, Immunescore and Estimatescore in TARGET-OS cohort. (G,H) Box plots show the differences of TumorPurity, Stromalscore, Immunescore and Estimatescore in different cluster of TARGET cohort. (I,J) GSEA outcome of different cluster. DGs, differentially expressed genes; GSEA, gene set enrichment analysis; NMF, non-negative matrix factor; OS, osteosarcoma; PPI, protein-protein interaction; RMS, rhabdomyosarcoma; TARGET, therapeutic applications research to generate effective treatments; TIS, Tumor Immune Score.

We then categorized the immunosuppressive subtypes and other subtypes into two groups: the “1 group” and the “other group”. GSEA results showed that in OS, the “other group” is significantly enriched in the lysosomal signaling pathway, apoptosis signaling pathway, and VEGF signaling pathway (Figure 1I). In RMS, the “other group” is significantly enriched in the purine metabolism-related signaling pathway, ribosome signaling pathway, and oxidative phosphorylation signaling pathway (Figure 1J).

Verification of the presence of immunosuppressive states

To further validate our hypothesis regarding the association between DGs and the “immunosuppressive” state, we analyzed data from the GEO database (GSE21257, GSE66533). We identified two subtypes of OS (Figure 2A), with subtype 1 characterized as the immunosuppressive state (Figure 2B). Similarly, among the four subtypes of RMS (Figure 2C), subtypes 2 and 3 were identified as immunosuppressive (Figure 2D).

Figure 2 Verification of immunosuppressive states. (A,C) The heatmap of TARGET samples by unsupervised NMF illustrates 2 and 3 distinct expression patterns. (B,D) Box plots show the differences of Stromalscore, Immunescore and Estimatescore in different clusters of the GEO cohort. (E,F) GSEA outcome of different cluster. (G) Box plots show the different expression levels of multiple inhibitory receptors between different clusters. Group 1, immunization group 1; Group 2, immunization group 2; Group 3, immunization group 3; Group 4, immunization group 4; Group other, other immunization group. GEO, Gene Expression Omnibus; GSEA, gene set enrichment analysis; NMF, non-negative matrix factor; OS, osteosarcoma; PD-1, programmed cell death protein-1; RMS, rhabdomyosarcoma; TARGET, therapeutic applications research to generate effective treatments.

We subsequently categorized the immunosuppressive subtypes and other subtypes into two groups: “1 group” for immunosuppressive subtypes and “other group” for the remaining subtypes. GSEA revealed that the immunosuppressive subtype (1 group) of OS was significantly enriched in the pyruvate metabolism pathway and the ubiquitin-mediated protein hydrolysis pathway. In contrast, the sphingolipid biosynthesis global series pathway was significantly enriched in the other subtypes (other group) (Figure 2E). For RMS, the other subtypes were significantly enriched in pathways associated with systemic lupus erythematosus, adipocytokine signaling, and Fc epsilon RI signaling (Figure 2F).

In addition, we analyzed the key markers of T-cell exhaustion, PD-1 (PDCD1) and CTLA-4, in the immunosuppressive subtypes of RMS and OS. We found that PD-1 and CTLA-4 are highly expressed in the immunosuppressive subtypes of both RMS and OS. When T cells are continuously exposed to antigen stimulation (such as in the tumor microenvironment), they upregulate inhibitory receptors like PD-1 and CTLA-4, which inhibits their function and prevents them from effectively killing tumor cells, thereby leading to immune escape. Therefore, this result further confirms the existence of the immunosuppressive subtypes of RMS and OS as we defined (Figure 2G). Interestingly, we found significant differences in T cell infiltration between the two subtypes of OS. However, no significant differences in immune infiltrating cells were observed between the two subtypes of RMS, possibly due to limited data (Figure 3A).

Figure 3 Combinations of machine learning algorithms. (A) The functionally grouped network of KEGG pathways by ClueGO for the interpretation of intersecting genes’ biological roles. (B) 100 combinations of machine learning algorithms; the numbers are the values of ROC curve and represents the accuracy of subtype classification of screened disulfidptosis-related genes. TARGET represents the training cohort, GEO represents the validation cohort, mean represents the average of the training cohort and validation cohort for each disease, Multicox-OS/RMS represents the ROC curve value corresponding to the disulfidptosis gene screened by the combination algorithm when constructing a multivariate Cox regression model in the TARGET cohort of the two diseases, the mean-risk score represents the average of the ROC curves of the two diseases, and mean-all is the average of all ROC curves. (C) The Ridge algorithm selects the top 20 genes for acute trait selection of intersection genes. ANN, Artificial Neural Network; GEO, Gene Expression Omnibus; KEGG, Kyoto Encyclopedia of Genes and Genomes; KNN, K-Nearest Neighbors; OS, osteosarcoma; ROC, receiver operating characteristic; RMS, rhabdomyosarcoma; SGD, Stochastic Gradient Descent; SVM, Support Vector Machine; TARGET, therapeutic applications research to generate effective treatments.

Furthermore, differentially expressed genes between immunosuppressive subtypes and other subtypes were enriched in signaling pathways such as nervous system development, actin filament-based processes, and VEGFA-VEGFR2 signaling (Figure 3A).

100 combinations of machine learning algorithms

Through the application of 100 machine learning algorithms (Figure 3B), we determined that the Ridge algorithm exhibited significantly higher comprehensive accuracy compared to the other four algorithms. Utilizing the Ridge algorithm, we screened feature genes for two disease types (Figure 3C) and identified five representative DGs: ACTN4, MYH9, FLNA, MYH10, and IQGAP1 (Figure 4A,4B). Based on these genes, we developed a risk score formula for OS:

Figure 4 Genes and model. (A) Disulfidptosis-related genes were analyzed in relation to the first 20 genes of OS and RMS, and the final 5 representative disulfidptosis genes were screened. (B) Calculated coefficients for five genes. (C) The PPI network among the 5 representative disulfidptosis genes. (D) Contribution values of five genes calculated by SHAP. (E) Differences in expression levels of IQGAP1 protein and PD-L1 protein between immunosuppressive subtypes and normal subtypes in OS and RMS. (F) OS TARGET cohort KM curve of overall survival time in days and event-free survival time in days. (G) RMS TARGET cohort KM curve of overall survival time in days and event free survival time in days. (H) GSVA of biological pathways between two distinct subgroups. *, P<0.05; **, P<0.01. CI, confidence interval; EFS, event-free survival; GSVA, Gene Set Variance Analysis; H, high; HR, hazard ratio; KM, Kaplan-Meier; L, low; OS, osteosarcoma; PD-1, programmed cell death protein-1; PPI, protein-protein interaction; RMS, rhabdomyosarcoma; SHAP, SHapley Additive exPlanations; TARGET, therapeutic applications research to generate effective treatments.

Risk score = (0.003651229) (ACTN4) − (0.005301413) (MYH9) − (0.003609198) (FLNA) − (0.004973973) (MYH10) + (0.003825834) (IQGAP1) (Figure 4C).

For RMS, the risk score formula is:

Risk score = (0.016037485) (ACTN4) − (0.027786077) (MYH9) + (0.008547981) (FLNA) − (0.003369) (MYH10) + (0.064156282) (IQGAP1) (Figure 4C).

Through SHAP model analysis, the feature gene with the greatest contribution, IQGAP1, was identified (Figure 4D). We identified a subset of cells with high IQGAP1 expression in OS and RMS cells, and simultaneously found that this subset also had high PD-L1 expression (Figure 4E). Due to the high expression of PD-L1, we preliminarily hypothesize that this subset is an immunosuppressive subset, and reasonably speculate that the immunosuppressive state of this subset may be associated with IQGAP1 in some way. The accuracy of these risk scores was confirmed by different KM survival curves (Figure 4F,4G). Interestingly, GSVA results (Figure 4H) indicated no significant difference in cellular signaling pathways between high and low-risk groups. This lack of difference may be attributed to bias in the selected database.

Comparison of different model genes

We collected 66 models that have been published so far (Figure 5A) and ultimately screened 54 models for multivariate Cox regression analysis in the TARGET database of RMS. The results indicated that approximately half of the models had a standard ROC curve >0.60, suggesting a similarity in molecular mechanisms between OS and RMS to some extent (Figure 5B). The risk formula composed of the five DGs we identified showed a standard ROC curve of 0.75 and 0.64 for the two diseases, with the 1-, 3-, and 5-year time-dependent ROC curve of 0.73, 0.75, 0.76 and 0.71, 0.73, 0.73, respectively (Figure 5B). Although this model did not exhibit outstanding advantages compared to other models, it demonstrated the feasibility of using DGs to predict the prognosis of the two diseases. Additionally, the Sankey plot results indicated no significant membership relationship between the different risk score groups and certain clinical characteristic variables corresponding to the two diseases (Figure 5C). Combined with the GSVA results, it appears that there is no strong correlation between the subtypes of immunosuppressive status and the high and low-risk groups distinguished by the risk scoring formula. Therefore, further in-depth analysis is required for the five genes ACTN4, MYH9, FLNA, MYH10, and IQGAP1.

Figure 5 Comparison of different model genes. (A) 66 published osteosarcoma models and the databases they used. (B) In the TARGET cohort of rhabdomyosarcoma 54 model genes, multivariate Cox regression analysis was performed, corresponding to standard ROC curve value ordering; and the time-dependent ROC curve values obtained after DGs performed multivariate Cox regression analysis. (C) Alluvial diagram of subgroup distributions in group with Immunization category, Risk score and Different clinical features. (D) Multivariate Cox regression analysis on ten variables (DGs, gender, age, stage, risk score and Immunization category) for RMS-TARGET cohort. (E) Single-gene Kaplan-Meier curve of IQGAP1. (F) IQGAP1 gene interactions. (G) Expression of the IQGAP1 gene in different cells in skeletal muscle. *, P<0.05; **, P<0.01; ***, P<0.001. AIC, Akaike Information Criterion; CI, confidence interval; DGs, differentially expressed genes; H, high; HR, hazard ratio; L, low; LncRNA, long non-coding RNA; mRNA, messenger RNA; NA, not available; nTPM, normalized Transcripts Per Million mapped reads; OS, osteosarcoma; ROC, receiver operating characteristic; RMS, rhabdomyosarcoma; TARGET, therapeutic applications research to generate effective treatments.

IQGAP1 molecular mechanism of action in RMS

IQGAP1 is a scaffold protein that plays a crucial role in regulating the actin cytoskeleton and cell migration (12-14). Our analysis of forest plots and KM survival curves indicates that IQGAP1 can serve as an independent prognostic factor for RMS (Figure 5D,5E). Notably, there are currently no relevant research reports on this finding. In our study, we identified the DGs ACTN4 and CD2AP in the interaction network of IQGAP1 (Figure 5F). This discovery further confirms the influence of IQGAP1 on the occurrence and development of RMS through the disulfidptosis mechanism. Additionally, we found that IQGAP1 is predominantly expressed in macrophages (Figure 5G). Previous studies have reported that IQGAP1 regulates macrophage infiltration and plays a crucial role in post-ischemic neovascularization (15). Furthermore, in the context of immune therapy resistance in hepatocellular carcinoma, extracellular vesicles secreted by tumor-associated macrophages can regulate PD-L1 expression through IQGAP1 (16). These findings provide theoretical support for the impact of IQGAP1 on the prognosis of RMS by affecting macrophages.

Role of macrophages in OS

In our study, we analyzed the cellular composition of OS, focusing on chondrocytes, endothelial cells, monocytes, macrophages, NK cells, and tissue stem cells. Among these, monocytes and macrophages were found to be the predominant cell types (Figure 6A). Notably, IQGAP1 and MYH9 were highly expressed in these cells. Monocle3 trajectory analysis revealed that monocytes and macrophages play a crucial role in the differentiation process, contributing to the formation of tissue stem cells and T cells (Figure 6B-6D). Interestingly, one group lacked a sufficient number of macrophages, likely due to sample selection bias. Consequently, we focused on the other two groups for cell communication analysis. Although the interaction between macrophages and monocytes was not the strongest (Figure 6B,6C), their ligand-receptor pairs, Siglec-15 and TNFRSF1B, were identified as potential targets for immunotherapy (Figure 6E,6F) (17,18). Additionally, LAIR1 was found to regulate the dynamics of monocytes and macrophages (19).

Figure 6 Single-cell transcriptomic analysis reveals cell-type-specific expression. (A) Distribution of six types of cells in the sample; the expression of DGs in each type of cell. (B-D) Differentiation trajectories of cell subsets, numbers with white circles indicate different endings; The numbers with black circles represent bifurcation points, from which cells can have multiple outcome numbers indicating the process of differentiation, the lower the number, the earlier the time of differentiation. (E,F) Ligand-receptor pairs for cell-cell interactions; the larger the point, the thicker the edge, the lighter the color in the circle, the greater the meaning. DGs, differentially expressed genes; NK, natural killer; tSNE, t-Distributed Stochastic Neighbor Embedding; UMAP, Uniform Manifold Approximation and Projection.

Discussion

Cell death is crucial for maintaining the dynamic balance of cell proliferation and metabolism in the body (20,21). In the tumor microenvironment, immune cell regulatory mechanisms enhance tumor cell apoptosis and necrosis, releasing chemical factors that influence the immune response. Disulfidptosis, a newly discovered mode of cell death, has not been studied in relation to OS and RMS. We conducted omics research on the gene set of disulfidptosis, rather than individual genes, to better understand its role in the immune microenvironment of these cancers. This approach provides valuable insights into the prognosis, molecular, and immune characteristics of OS and RMS, and proposes new strategies for optimizing clinical treatment.

To explore the relationship between disulfidptosis and the tumor immune microenvironment of RMS and OS, we used NMF to extract biological correlation coefficients from gene expression matrix data in patient cohorts. By organizing DGs and samples, we identified internal structural characteristics and grouped the samples (22). We discovered a special subtype in both OS and RMS characterized by low immune cell infiltration and low immune response efficiency. Previous studies have reported an immunosuppressive state in these cancers (23,24), leading us to define this subtype as “immunosuppressive”. In OS, the immunosuppressive subtype identified by disulfidptosis showed fewer T cells compared to other subtypes. However, the GSE21257 cohort revealed a significant increase in resting memory CD4⁺ T cells in this subtype. We speculate that the complex differences in memory T cells across environments, where cysteine accumulates during disulfidptosis (25), may create an acidic environment that disrupts T cell immunotherapy for OS (26), leading to the emergence of immunosuppressive subtypes.

We further analyzed the different enrichment levels of cellular signaling pathways between immunosuppressive and other subtypes to explore the relationship between disulfidptosis and immunosuppressive subtypes. In OS, we found that VEGF-related signaling pathways were enriched. VEGF family members are secreted dimeric glycoproteins involved in regulating angiogenesis (27,28). VEGFA and VEGF regulate angiogenesis and proliferation in OS and RMS (29-31), and previous studies have reported a correlation between VEGF and macrophages (32,33).

Single-cell database analysis of OS also revealed macrophages as a major cellular subgroup. Macrophages play a crucial role in the immune response. In RMS, cells are sensitive to macrophage-mediated cytotoxicity and secrete macrophage migration inhibitory factor (MIF) under cytotoxic drug induction, regulating tumor cell metastasis and inhibiting tumor-associated fibroblast recruitment (34-36). In OS, macrophages differentiate from monocytes, a process closely related to regulatory T cell (Treg) differentiation (37). Tumor-associated macrophages promote tumor growth and angiogenesis, upregulate tumor stem cell-like phenotypes (38), and their diversity makes macrophage repolarization a hot topic in immunotherapy (39). SLC7A11 is involved in macrophage ferroptosis (40), and macrophages release exosomes targeting SLC7A11 during heart injury, inducing ferroptosis by inhibiting glutathione synthesis (41,42). Inhibiting glutathione synthesis increases cysteine and cystine levels, meeting the conditions for disulfidptosis during glucose starvation (10). Our bioinformatic analysis identified a significant correlation between the expression of IQGAP1—an independent prognostic factor in RMS—and markers of macrophage infiltration. Furthermore, this association suggests a potential role for IQGAP1 in macrophage recruitment or polarization; however, this speculative link requires experimental validation.

There are limitations in this study. The collection of DGs is not comprehensive, leading to potential deviations in results. There may be selectivity bias in our database selection, and the sample size is insufficient, causing randomness and affecting result credibility. Additionally, the lack of a single-cell common database for RMS limits further in-depth analysis. Detailed in vivo and in vitro experiments are needed to supplement our findings on independent prognostic factors for RMS and the complex role of macrophages in OS and RMS.

Conclusions

Based on integrated multi-omics analysis and machine learning algorithms, this study identifies and validates a novel disulfidptosis-related immunosuppressive subtype in both OS and RMS, characterized by diminished immune infiltration and dysfunctional microenvironment. We further established a robust five-gene prognostic signature (ACTN4, MYH9, FLNA, MYH10, IQGAP1), which effectively stratifies patients into distinct risk groups with significant survival differences. Notably, IQGAP1 serves as an independent prognostic factor in RMS, suggesting its potential as both a therapeutic target and clinical biomarker.

In conclusion, our results provide a foundational framework for leveraging DGs to refine prognostic prediction and inform immunotherapeutic strategies for pediatric sarcomas, ultimately advocating for more personalized treatment approaches.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the STREGA and MDAR reporting checklists. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1439/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1439/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1439/prf

Funding: None.

Conflicts of Interest: Both authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1439/coif). 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. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

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. Hettmer S, Linardic CM, Kelsey A, et al. Molecular testing of rhabdomyosarcoma in clinical trials to improve risk stratification and outcome: A consensus view from European paediatric Soft tissue sarcoma Study Group, Children's Oncology Group and Cooperative Weichteilsarkom-Studiengruppe. Eur J Cancer 2022;172:367-86. [Crossref] [PubMed]
2. Hu Z, Wen S, Huo Z, et al. Current Status and Prospects of Targeted Therapy for Osteosarcoma. Cells 2022;11:3507. [Crossref] [PubMed]
3. Singh S, Abu-Zaid A, Jin H, et al. Targeting KDM4 for treating PAX3-FOXO1-driven alveolar rhabdomyosarcoma. Sci Transl Med 2022;14:eabq2096. [Crossref] [PubMed]
4. Breakthrough of the year 2013. How we did in 2013 and. Science 2013;342:1442. [Crossref] [PubMed]
5. Chen CY, Wang PY, Hutzen B, et al. Cooperation of Oncolytic Herpes Virotherapy and PD-1 Blockade in Murine Rhabdomyosarcoma Models. Sci Rep 2017;7:2396. [Crossref] [PubMed]
6. Zheng B, Ren T, Huang Y, et al. PD-1 axis expression in musculoskeletal tumors and antitumor effect of nivolumab in osteosarcoma model of humanized mouse. J Hematol Oncol 2018;11:16. [Crossref] [PubMed]
7. Hegde M, Joseph SK, Pashankar F, et al. Tumor response and endogenous immune reactivity after administration of HER2 CAR T cells in a child with metastatic rhabdomyosarcoma. Nat Commun 2020;11:3549. [Crossref] [PubMed]
8. Jiang X, Stockwell BR, Conrad M. Ferroptosis: mechanisms, biology and role in disease. Nat Rev Mol Cell Biol 2021;22:266-82. [Crossref] [PubMed]
9. Zhao R, Wang W, Pan L, et al. The prognostic value and response to immunotherapy of immunogenic cell death-associated genes in breast cancer. Front Oncol 2023;13:1047973. [Crossref] [PubMed]
10. Liu X, Nie L, Zhang Y, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol 2023;25:404-14. [Crossref] [PubMed]
11. Kitamoto K, Miura Y, Karnan S, et al. Inhibition of NADPH oxidase 2 induces apoptosis in osteosarcoma: The role of reactive oxygen species in cell proliferation. Oncol Lett 2018;15:7955-62. [Crossref] [PubMed]
12. Benseñor LB, Kan HM, Wang N, et al. IQGAP1 regulates cell motility by linking growth factor signaling to actin assembly. J Cell Sci 2007;120:658-69. [Crossref] [PubMed]
13. Brown MD, Sacks DB. IQGAP1 in cellular signaling: bridging the GAP. Trends Cell Biol 2006;16:242-9. [Crossref] [PubMed]
14. Le Clainche C, Schlaepfer D, Ferrari A, et al. IQGAP1 stimulates actin assembly through the N-WASP-Arp2/3 pathway. J Biol Chem 2007;282:426-35. [Crossref] [PubMed]
15. Urao N, Razvi M, Oshikawa J, et al. IQGAP1 is involved in post-ischemic neovascularization by regulating angiogenesis and macrophage infiltration. PLoS One 2010;5:e13440. [Crossref] [PubMed]
16. Wang X, Ye X, Chen Y, et al. Mechanism of M2 type macrophage-derived extracellular vesicles regulating PD-L1 expression via the MISP/IQGAP1 axis in hepatocellular carcinoma immunotherapy resistance. Int Immunopharmacol 2023;124:110848. [Crossref] [PubMed]
17. Wang J, Sun J, Liu LN, et al. Siglec-15 as an immune suppressor and potential target for normalization cancer immunotherapy. Nat Med 2019;25:656-66. [Crossref] [PubMed]
18. Medler J, Kucka K, Wajant H. Tumor Necrosis Factor Receptor 2 (TNFR2): An Emerging Target in Cancer Therapy. Cancers (Basel) 2022;14:2603. [Crossref] [PubMed]
19. Keerthivasan S, Şenbabaoğlu Y, Martinez-Martin N, et al. Homeostatic functions of monocytes and interstitial lung macrophages are regulated via collagen domain-binding receptor LAIR1. Immunity 2021;54:1511-1526.e8. [Crossref] [PubMed]
20. Villa E, Ricci JE. How does metabolism affect cell death in cancer? FEBS J 2016;283:2653-60. [Crossref] [PubMed]
21. Bao JH, Lu WC, Duan H, et al. Identification of a novel cuproptosis-related gene signature and integrative analyses in patients with lower-grade gliomas. Front Immunol 2022;13:933973. [Crossref] [PubMed]
22. Lee DD, Seung HS. Learning the parts of objects by non-negative matrix factorization. Nature 1999;401:788-91. [Crossref] [PubMed]
23. Liu W, Hu H, Shao Z, et al. Characterizing the tumor microenvironment at the single-cell level reveals a novel immune evasion mechanism in osteosarcoma. Bone Res 2023;11:4. [Crossref] [PubMed]
24. DeMartino J, Meister MT, Visser LL, et al. Single-cell transcriptomics reveals immune suppression and cell states predictive of patient outcomes in rhabdomyosarcoma. Nat Commun 2023;14:3074. [Crossref] [PubMed]
25. Radbruch A, McGrath MA, Siracusa F, et al. Homeostasis and Durability of T-Cell Memory-The Resting and the Restless T-Cell Memory. Cold Spring Harb Perspect Biol 2021;13:a038083. [Crossref] [PubMed]
26. Park JA, Cheung NV. Promise and Challenges of T Cell Immunotherapy for Osteosarcoma. Int J Mol Sci 2023;24:12520. [Crossref] [PubMed]
27. Olsson AK, Dimberg A, Kreuger J, et al. VEGF receptor signalling - in control of vascular function. Nat Rev Mol Cell Biol 2006;7:359-71. [Crossref] [PubMed]
28. Melincovici CS, Boşca AB, Şuşman S, et al. Vascular endothelial growth factor (VEGF) - key factor in normal and pathological angiogenesis. Rom J Morphol Embryol 2018;59:455-67.
29. Vimalraj S, Saravanan S, Raghunandhakumar S, et al. Melatonin regulates tumor angiogenesis via miR-424-5p/VEGFA signaling pathway in osteosarcoma. Life Sci 2020;256:118011. [Crossref] [PubMed]
30. Zhang L, Lv Z, Xu J, et al. MicroRNA-134 inhibits osteosarcoma angiogenesis and proliferation by targeting the VEGFA/VEGFR1 pathway. FEBS J 2018;285:1359-71. [Crossref] [PubMed]
31. Wan X, Shen N, Mendoza A, et al. CCI-779 inhibits rhabdomyosarcoma xenograft growth by an antiangiogenic mechanism linked to the targeting of mTOR/Hif-1alpha/VEGF signaling. Neoplasia 2006;8:394-401. [Crossref] [PubMed]
32. Ganta VC, Choi M, Farber CR, et al. Antiangiogenic VEGF(165)b Regulates Macrophage Polarization via S100A8/S100A9 in Peripheral Artery Disease. Circulation 2019;139:226-42. [Crossref] [PubMed]
33. Kim KH, Kim YH, Son JE, et al. Intermittent fasting promotes adipose thermogenesis and metabolic homeostasis via VEGF-mediated alternative activation of macrophage. Cell Res 2017;27:1309-26. [Crossref] [PubMed]
34. Tarnowski M, Grymula K, Liu R, et al. Macrophage migration inhibitory factor is secreted by rhabdomyosarcoma cells, modulates tumor metastasis by binding to CXCR4 and CXCR7 receptors and inhibits recruitment of cancer-associated fibroblasts. Mol Cancer Res 2010;8:1328-43. [Crossref] [PubMed]
35. Johler SM, Fuchs J, Seitz G, et al. Macrophage migration inhibitory factor (MIF) is induced by cytotoxic drugs and is involved in immune escape and migration in childhood rhabdomyosarcoma. Cancer Immunol Immunother 2016;65:1465-76. [Crossref] [PubMed]
36. Herrmann D, Seitz G, Fuchs J, et al. Susceptibility of rhabdomyosarcoma cells to macrophage-mediated cytotoxicity. Oncoimmunology 2012;1:279-86. [Crossref] [PubMed]
37. Mujal AM, Combes AJ, Rao AA, et al. Holistic Characterization of Tumor Monocyte-to-Macrophage Differentiation Integrates Distinct Immune Phenotypes in Kidney Cancer. Cancer Immunol Res 2022;10:403-19. [Crossref] [PubMed]
38. Zhao Y, Zhang B, Zhang Q, et al. Tumor-associated macrophages in osteosarcoma. J Zhejiang Univ Sci B 2021;22:885-92. [Crossref] [PubMed]
39. Anand N, Peh KH, Kolesar JM. Macrophage Repolarization as a Therapeutic Strategy for Osteosarcoma. Int J Mol Sci 2023;24:2858. [Crossref] [PubMed]
40. Bao X, Luo X, Bai X, et al. Cigarette tar mediates macrophage ferroptosis in atherosclerosis through the hepcidin/FPN/SLC7A11 signaling pathway. Free Radic Biol Med 2023;201:76-88. [Crossref] [PubMed]
41. Zhao X, Si L, Bian J, et al. Adipose tissue macrophage-derived exosomes induce ferroptosis via glutathione synthesis inhibition by targeting SLC7A11 in obesity-induced cardiac injury. Free Radic Biol Med 2022;182:232-45. [Crossref] [PubMed]
42. Wu Q, Liang Z, Jiang J, et al. Macrophages originated IL-33/ST2 inhibits ferroptosis in endometriosis via the ATF3/SLC7A11 axis. Cell Death Dis 2023;14:668. [Crossref] [PubMed]