The mRNA expression data and clinical details for osteosarcoma patients (TARGET-OS dataset) were obtained from the Cancer Genome Atlas Program (TCGA) database (https://portal.gdc.cancer.gov/). The RNA-seq raw read count from the TCGA database was converted to transcripts per kilobase million (TPM) and further log-2 transformed. The GSE21257 dataset from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) supplied as an external validation set mRNA expression data and clinical information for osteosarcoma specimens. The CIBERSORT algorithm was employed to quantify the presence of 22 infiltrating immune cell types within the TME using the TARGET-OS dataset (23). Pearson correlation analysis identified 281 genes exhibiting a correlation with M1 macrophage expression (|R2| > 0.4, p < 0.001). Samples with a survival duration of less than 30 days were excluded from the analysis.

A total of 84 samples, consisting of survival and expression data from the TARGET-OS dataset, were divided into a training set (n=59) and an internal validation set (n=25) at a 7:3 ratio. The GSE21257 dataset provided 53 samples for an external validation cohort. In the training set, univariate Cox regression analysis identified M1 macrophage genes with prognostic relevance. The least absolute shrinkage and selection operator (LASSO) algorithm, with 1000 iterations, was then employed to select the optimal subset of prognostic genes, leading to the development of an M1 macrophage gene signature (MRS). The multivariate Cox regression model determined the final genes after LASSO algorithm application. Risk scores were calculated using the linear combination of each chosen gene, following the formula: Risk score = ∑(coef (β) * EXP(β)), where β denotes the regression coefficient. Patients were categorized into high- and low-risk groups based on the median risk score as the threshold, and the clinical differences between these groups were investigated using Kaplan-Meier survival analysis. The predictive accuracy of the model was evaluated via the ROC curve and C-index. Moreover, stratified analysis was performed to assess the additional prognostic value of the MRS.

The patients of the internal and external validation set were subjected to the identical grouping methodology utilized in the training set, after which their survival was assessed through the application of Kaplan-Meier survival analysis and risk plot.

A nomogram was constructed to forecast the 1-, 3-, and 5-year survival rates of OS patients using the risk score in conjunction with the clinicopathological factors of age, gender, race, and metastasis. The accuracy of the nomogram’s predictions was then tested using a calibration curve to compare actual overall survival with predicted survival rates.

The cohort was partitioned into high- and low-risk groups using the predetermined risk score threshold. Subsequently, gene expression fold changes were analyzed using the R package “limma”. Pathway analysis was performed using the R package “clusterProfiler” and focused on identifying significantly enriched pathways within the reference gene set for the high- and low-risk groups. The reference gene set was defined as the hallmark gene sets described by Subramanian et al. (24).

Single-sample gene set enrichment analysis (ssGSEA) algorithm using R packages, specifically limma, GSVA, and GSEABase, was used to evaluate the disparity in immune function between high- and low-risk groups classified according to MRS (25). TME and immune cell infiltration were evaluated using the ESTIMATE and CIBERSORT algorithms to determine the proportions of its components (23, 26). Subsequently, the association between the expression levels of immune checkpoint genes and the two groups was investigated.

The Genomics of Drug Sensitivity in Cancer (GDSC) public repository offers valuable insights into cancer cell drug sensitivity and associated molecular markers for drug responses (27). By employing the oncoPredict package, gene expression profiles from GDSC2 and corresponding drug response data were obtained to build a ridge regression model tailored for osteosarcoma transcriptomic information. Following this, sensitivity scores were calculated to estimate the half-maximal inhibitory concentration (IC50) for various drugs in the context of OS patients.

Human OS cell lines (143B, MG63, and U2OS) and human normal osteoblast cell line (hFOB1.19) were purchased from Wuhan Servicebio Technology Co., Ltd. (Wuhan, China). Each cell line was cultured in its specific medium (Wuhan Servicebio Technology Co. Ltd., Wuhan, China). Human OS cell lines and hFOB1.19 cells were cultured, respectively, at 37°C in humidified incubator with 5% CO₂ and at 34°C in an incubator with 5% CO₂.

Total RNA was extracted from OS cell lines and the normal osteoblast cell line, hFOB1.19, employing TRIzol reagent (Invitrogen, USA). First-strand cDNAs were synthesized using the PrimeScript™ RT reagent kit (Takara, Japan) as per the manufacturer’s guidelines. MRG mRNA levels were assessed through qRT-PCR, utilizing SYBR Premix Ex Taq (TaKaRa, Japan) and specific primers ( Table 1 ). The relative expression levels were determined and adjusted to the reference control GAPDH using the 2^–ΔΔCt method.

Total protein from OS cell lines and hFOB1.19 was isolated using RIPA buffer (Beyotime, China). Protein samples (20 μg each) were then separated using SDS-PAGE and transferred onto PVDF membranes. Following this, the membranes were blocked using 5% non-fat milk in TBST for 1 hour at room temperature and subsequently incubated with specific primary antibodies at 4˚C overnight. The employed antibodies included: CD37 (Abcam, ab300400, 1:1,000), ARHGAP25 (Abcam, ab181202, 1:10,000), GABRD (Abcam, ab93619, 1:1,000), β-Tubulin (Cell Signaling Technology, #2146, 1:1,000), and GAPDH (Cell Signaling Technology, #5174, 1:1,000).

The statistical software R version 4.1.2 and Prism 8.0 (GraphPad) were utilized to conduct data analysis and generate visual representations of the findings. To assess the associations between variables, Spearman or Pearson correlation coefficients were calculated. Kaplan-Meier methodology was utilized to construct survival curves, and the log-rank test was performed to compare the curves. Univariable and multivariable Cox regression models were used to identify prognostic factors for overall survival. Experimental data were presented as mean ± standard error of the mean (SEM). One-way ANOVA followed by Tukey post hoc analysis was used to compare 3 or more groups. P values of <0.05 were considered to be statistically significant.

Figure 1 illustrates the study’s workflow. The CIBERSORT algorithm was employed to analyze the M1 macrophage subpopulation abundance in each sample. A total of 281 genes exhibiting a correlation with M1 macrophages were pinpointed through Pearson correlation analysis (|R2| > 0.4, p < 0.001) ( Figure 1 , Supplementary Table S1 ). Forty-nine genes associated with M1 macrophages were recognized as potential prognostic indicators via univariate Cox analysis ( Figure 1 ). To mitigate the risk of overfitting, LASSO Cox regression was subsequently performed ( Figures 1A–E ). After applying the LASSO algorithm, the multivariate Cox regression model was used to identify the final gene set, which consisted of 3 robust genes, forming a prognostic signature for overall survival.

Coefficients of the three M1 macrophage-associated genes were employed to determine scores for each patient. The risk score computation was as follows: Risk score = (-2.284 × CD37 expression) + (3.845 × GABRD expression) + (-3.632 × ARHGAP25 expression). Subsequently, participants were assigned to low- or high-risk groups based on the median value of the risk score. In the training set, high-risk patients demonstrated shorter overall survival than low-risk counterparts (p < 0.001, Figure 2 ). Similar trends were observed in both internal and external validation cohorts (p < 0.001, Figures 2A–C ).

To assess the signature’s effectiveness, a time-dependent ROC curve was utilized. AUC values for 1-year, 3-year, and 5-year survival periods were 0.746, 0.839, and 0.850, respectively ( Figure 2 ). The one-year survival rate AUC suggested that both risk score (0.850) and metastasis (0.694) had satisfactory predictive capabilities ( Figure 2 ). As depicted in Figure 2 , the risk score’s C-index surpassed that of clinicopathological factors.

Univariate and multivariate Cox regression analyses were performed to evaluate the prognostic value of risk score and other factors. The risk score and metastasis emerged as significant independent prognostic factors, as evidenced by the HR and CI values: 4.905 (95% CI = 2.677–8.989, p < 0.001) and 4.516 (95% CI = 2.044-9.979, p < 0.001) for univariate analysis, and 4.669 (95% CI = 2.455-8.879, p < 0.001) and 3.392 (95% CI = 1.478-7.782, p = 0.004) for multivariate analysis ( Figures 2A–H ).

In order to furnish a numerical means of clinical utilization, a nomogram was devised incorporating factors such as age, gender, race, metastasis, and risk score to forecast the overall survival of patients ( Figure 3 ). The calibration plot demonstrates a high degree of agreement between the observed and predicted rates of 1, 3, and 5-year overall survival ( Figure 3 ). Finally, the applicability of MRS was evaluated by grouping patients based on their age, gender, race and metastasis. The results showed that patients with high-risk scores within each subgroup had an unfavorable prognosis, demonstrating the efficacy of MRS in predicting outcomes for all patients ( Figures 3A–F ).

A comprehensive analysis of GO and KEGG methodologies was carried out to explore the potential mechanisms responsible for the differing prognoses observed between high- and low-risk groups. A total of 631 differentially expressed genes were identified between the two groups. GO enrichment analysis of biological processes (BP) indicated that the DEGs were predominantly involved in “leukocyte activation regulation,” “positive modulation of lymphocyte activation,” and “leukocyte-driven immunity.” In relation to cellular components (CC), the DEGs were chiefly linked to “immunoglobulin complexes,” “external aspect of plasma membrane,” and “circulating immunoglobulin complexes.” Likewise, molecular function (MF) analysis revealed that the DEGs primarily focused on “antigen binding” and “immunoglobulin receptor binding” ( Figure 4 ). KEGG analysis results showed that the DEGs were mainly enriched in several pathways, including “Th1 and Th2 cell differentiation”, “Cytokine-cytokine receptor interactions,” “Antigen processing and presentation,” “Osteoclast differentiation,” “PD-L1 expression and PD-1 checkpoint pathway in cancer”, and “NF-kappa B signaling pathway” ( Figure 4 ).

In order to further investigate the relationship between risk score and infiltration of immune cells into tumors, the CIBERSORT algorithm was employed to compare the proportions of 22 different types of immune cells between groups of individuals classified as either low or high risk. The findings revealed that the low-risk group had higher fractions of plasma cells, CD8+ T cells, regulatory T cells, and memory CD4+ T cells (p < 0.05), while the high-risk group had higher fractions of M0 macrophages, which is associated with immunosuppressive activity (p < 0.05) ( Figure 5 ). There was an indication of a greater prevalence of immune and stromal cells within the low-risk group according to the ESTIMATE algorithm ( Figures 5A–C ). The ssGSEA algorithm was utilized to deduce the immune function. Figure 5 reveals that there exists a significant disparity between the two groups in immune function. These findings suggest that the low-risk group exhibits a higher level of immune function activity. Then we investigated the potential association between risk scores and the expression levels of immune checkpoint genes. Patients categorized as low-risk exhibited significantly elevated expression levels of 26 immune checkpoint genes, namely CD274, HAVCR2, SELPLG, LAG3, CD27, ICOS, TIGHT, TNFRSF9, TNFRSF14, CD28, LGALS9, CD80, TNFRSF15, NRP1, CD40, TNFSF14, CD86, KIR3DL1, CD48, LAIR1, CD40LG, TMIGD2, CD200R1, CD44, CD96, SIGLEC7, while TNFSF9 was found elevated expression in the high-risk group ( Figure 5 ).

We investigated the correlation between the risk score and the effectiveness of targeted therapy and chemotherapy for OS patients. The findings indicated a positive correlation between the IC50 of Dasatinib and the risk score. In contrast, a negative correlation was observed between the IC50 of Daporinad, Linsitinib, Sabutoclax, and Dihydrorotenone and the risk score ( Supplementary Figure S1 ).

In order to investigate the disparities in gene expression patterns between tumor and normal cell lines, we chose three OS cell lines to determine their expression levels of mRNA and protein, with the normal osteoblast hFOB1.19 serving as the control group. Western blotting results indicate that CD37 protein expression was markedly higher in MG63 and U2OS cell line when compared with normal osteoblast hFOB1.19, and GABRD expression levels were significantly upregulated in U2OS cell line. In addition, the protein expression levels of ARHGAP25 were elevated in both 143B and U2OS cell lines. Figures 6A illustrate these findings.

The mRNA expression levels of CD37, GABRD and ARHGAP25 were analyzed by qRT-PCR ( Figure 6 ). The outcomes revealed that the levels of CD37, GABRD and ARHGAP25 were significantly elevated in the tumor cells than that in the normal cell line.