The Gene Expression Omnibus (GEO) is a functional genomics data repository, maintained by the National Center for Biotechnology Information (NCBI). To identify DEGs associated with OS, we selected three datasets: GSE12865 (17), GSE14359 (18), and GSE36001 (19). The inclusion criteria for the selected datasets were as follows: i) Tumor samples included primary or metastatic OS tissue. ii) Control samples included normal human bone tissue or osteoblasts. iii) Statistical thresholds for DEGs were p < 0.05 and |log2 fold change| > 1 iv) At least 10 overlapping DEGs in multiple datasets were required.

DEGs were identified by GEO2R, an interactive web tool that compares two or more sample sets in GEO and ranks genes based on statistical significance. We used |log2FC| ≥ 1 and p < 0.05 as thresholds. Venn plots were generated using the Venn tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) to visualize overlapping DEGs.

Gene Ontology (GO) was analyzed to explore the enrichment of DEGs in three major categories: biological process (BP), cellular component (CC) and molecular function (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used to investigate functional interactions at the molecular level (20). Both analyses were performed by the Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.ncifcrf.gov/), a widely used bioinformatics tool for systematic analysis of large-scale gene and protein datasets.

A PPI network was created using the STRING (Search Tool for Retrieving Gene/Protein Interactions) database and visualized through Cytoscape software. We then identified the most critical hub genes using the Molecular Complexity Detection (MCODE), an algorithm that allows for highly interrelated genes for scoring and screening. Moreover, the CytoHubba algorithm was used to further identify hub genes based on different topological properties, which can detect key hub genes in the network, thus providing a comprehensive view of molecular interactions.

Box plots were created from RNA sequencing (RNA-seq) data from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases using the Gene Expression Profiling Interactive Analysis 2.0 (GEPIA2) platform (http://gepia2.cancer-pku.cn/). Gene expression values were normalized to transcripts per million (TPM), log2 transformed, and compared between tumor and normal tissues. Statistical significance was assessed via the Wilcoxon rank sum test, with significant differences indicated by an asterisk.

The prognostic significance of central genes was assessed using the Kaplan-Meier (KM) plotter (https://kmplot.com/analysis/), which integrates survival data from GEO and TCGA. OS patient samples were divided into high and low expression groups according to the expression levels of hub genes. Kaplan-Meyer survival curves, hazard ratios (HR), 95% confidence intervals (CI), and log-rank P-values were calculated to assess the potential prognostic value of centrality.

The abundance of TIIC in OS samples was inferred from the TIMER (https://cistrome.shinyapps.io/timer/) database, an integrated platform for immune cell infiltration analysis. TIMER used an inverse plethysmography approach to infer tumor-infiltrating immune cells the gene expression profiles of samples from various cancer types in the TCGA (B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils and dendritic cells) abundance.

Single cell suspensions were processed using the 10× Chromium Single Cell Platform (10× Genomics, 30 v3 chemistry). Individual cell mRNAs were captured, cDNAs synthesized, and libraries constructed according to the manufacturer’s protocol. Sequencing was performed on the NovaSeq 6000 platform (Illumina) with 150 bp paired-end reads. Data were normalized using logarithmic transformation. Raw data were processed using the Cell Ranger 3.0.1 pipeline for alignment, data processing, and initial clustering. Single-cell data were obtained from GEO (GSE162454). scMetabolism (v0.2.1) pipeline was used to assess the metabolic diversity of each cluster. It scores each cluster and calculates activity scores for metabolic pathways based on single-cell matrix files.

Total RNA was extracted using TRIzol reagent (Invitrogen, CA, USA), and then reverse-transcribed into cDNA using the Evo M-MLV RT Premix Kit (Accurate Biology, Hunan, China). Real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed using the LightCycler^®480 (Roche, USA) system and SYBR^®Green Premix Pro Taq HS qPCR Kit (Accurate Biology, Hunan, China). The procedures above were performed according to the manufacturer’s instructions.

Protein samples were prepared from the collected tissues, subjected to SDS-polyacrylamide gel electrophoresis, and transferred to nitrocellulose membranes. The membranes were blocked with 5% non-fat milk at room temperature for 2 hours, followed by incubation with the purchased primary antibody overnight. Subsequently, the membranes were washed three times with Tris-buffered saline (TBS) containing Tween 20, and then incubated with the secondary antibody at room temperature for 1–2 hours. After washing the membrane three times, the protein expression levels were detected using chemiluminescence reagents. The antibodies used were as follows: anti-CYBB antibody (1:1000; ABclonal, A19701); anti-PLEK antibody (1:1000; ABclonal, A6305); anti-Vinculin (1:10000; Abways, CY5164).

Cell proliferation was evaluated using CCK-8 assays (Yeasen Biotechnology). Briefly, OS cells (1×10⁵/well) were plated in 96-well plates and cultured for 24–96 h. After adding 10 μL CCK-8 solution per well, plates were incubated at 37°C for 1 h. Absorbance at 450 nm was measured using a microplate reader (SpectraMax i3x).

Cells were seeded at a density of 600 cells/well in 6-well plates and cultured for 10–14 days, with medium changes performed as needed. Following colony formation, cells were washed twice with PBS, fixed with 4% paraformaldehyde (15min), and stained with 1% crystal violet (30min). Colonies containing >50 cells were counted under a light microscope. All experiments were conducted in triplicate.

Cell proliferation was analyzed using the Cell Light™ EdU Apollo567 In Vitro Kit (RiboBio). Briefly, OS cells were co-cultured seeded in 96-well plates for 48 h, followed by incubation with 50 μM EdU for 2 h at 37°C. Cells were then fixed with 4% formaldehyde (30 min), treated with 2 mg/mL glycine (5min), and permeabilized with 0.5% Triton X-100 (10min). After PBS washing, samples were incubated with Apollo reaction cocktail (30min), washed twice with 0.5% Triton X-100, and counterstained with Hoechst 33342 (30min). Fluorescence images were acquired using an Olympus microscope.

Cell migration and invasion were assessed using Transwell chambers (8 µm pore size; Corning). For invasion assays, chambers were pre-coated with 50 µL of 1:8 diluted Matrigel (Corning), whereas uncoated chambers were used for migration assays. Cells (1×10⁵) suspended in 200 µL serum-free medium were seeded in the upper chamber, while the lower chamber contained medium supplemented with 10% FBS as a chemoattractant. After 24 h incubation at 37°C, cells that migrated/invaded to the lower chamber were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet (1h, RT), and quantified by counting five random fields per well under a microscope. All experiments were performed in triplicate.

The objective of this study is to elucidate the key molecular features of OS and the flowchart is shown in Figure 1 . We selected three publicly available datasets from the GEO database (GSE12865, GSE14359, and GSE36001) and the statistical thresholds for DEGs were |log2FC|≥1 and p < 0.05. The results are as follows: 4,243 DEGs in GSE12865 (1,712 up-regulated, 2,531 down-regulated), 3,322 DEGs in GSE14359 (1,883 up-regulated, 1,439 down-regulated) and 890 DEGs in GSE36001 (271 up-regulated, 619 down-regulated) ( Figures 2A–C ). We then identified 159 DEGs overlapping in the three datasets through Venn diagrams ( Figure 2D ), which are potentially significant in the pathogenesis of OS. These genes may represent core molecular drivers of OS progression and require further validation.

To explore the biological context of the identified DEGs, we conducted GO and KEGG pathway enrichment analyses. These analyses revealed that DEGs were predominantly enriched in pathways related to tumor progression and immune regulation, including cell adhesion, apoptotic processes, extracellular matrix organization, and actin cytoskeleton remodeling ( Supplementary Figure S1 ). These results provided useful biological insight and served as a foundation for downstream analysis. While such enrichment patterns are frequently observed in various cancer types, the findings supported the relevance of the selected DEGs to core tumor-associated processes and helped to guide subsequent hub gene identification.

To explore protein interactions among identified DEGs, we built a PPI network with the STRING database and visualized it by Cytoscape software ( Figure 3A ). The network consisted of 134 nodes and 806 edges. This analysis then evaluated the interactions between these genes and identified 14 hub genes: C1QA, CD36, CD86, CXCL12, CXCR4, CYBB, ENG, FCER1G, GZMA, ITGAM, LAPTM5, PECAM1, PLEK and SELL ( Figure 3B ). The MCODE algorithm was applied in Cytoscape to detect densely connected modules, and five key modules were finally identified, of which the highest scoring module is shown in Figure 3C , consisting of 21 nodes and 153 edges. The other modules are shown in Figures 3D–G . The intersection results of MCODE and CytoHubba algorithm identified nine key hub genes: CD36, CD86, CXCL12, CXCR4, CYBB, FCER1G, ITGAM, LAPTM5, and PLEK ( Figure 3H ).

We assessed the expression levels of these nine hub genes in normal and OS tissues by data from TCGA and GTEx databases to identify their potential as prognostic genes. Eight genes (CD36, CD86, CXCR4, CYBB, FCER1G, ITGAM, LAPTM5, and PLEK) were observed upregulated in OS tissues, while CXCL12 was downregulated. Notably, the CXCL12-CXCR4 axis may affect tumor cell migration and the immune microenvironment (21). Elevated expression of CD36 and CYBB indicates their significant involvement in regulating lipid metabolism and oxidative stress during OS development (22, 23). The upregulation of immune-related genes such as CD86, FCER1G, and ITGAM suggests involvement of immune activation and inflammatory signaling within the osteosarcoma microenvironment. PLEK, which encodes a signaling adaptor protein, also showed a notable upregulation, meriting further investigation ( Supplementary Figure S2 ). However, given the extremely limited number of normal control samples, the expression comparison based on public datasets should be regarded as exploratory in nature.

To further support the differential expression of hub genes between OS and normal tissues, we performed an additional expression analysis using our recently published OS transcriptomic dataset (24). Consistent with our observations from TCGA and GTEx, most hub genes—including PLEK, CYBB, CXCR4, and ITGAM—were significantly upregulated in osteosarcoma samples compared to normal controls ( Supplementary Figure S3 ). These independent data reinforce the robustness of the expression trends and further support the relevance of these hub genes in the OS context.

To evaluate the prognostic significance of these hub genes, we performed a Kaplan-Meier survival analysis. Six genes (CD36, CXCL12, CXCR4, CYBB, ITGAM, and PLEK) were significantly associated with overall survival in OS patients. Specifically, CD36 [HR = 0.61 (0.41-0.91), p = 0.015], CXCL12 [HR = 0.56 (0.37-0.85), p = 0.0054], CXCR4 [HR = 0.66 (0.44-0.99), p = 0.043], CYBB [HR = 0. 64 (0.43-0.96), p = 0.031], ITGAM [HR = 0.67 (0.45-1.00), p = 0.049], and PLEK [HR = 0.58 (0.39-0.86), p = 0.0063] were associated with significantly better prognosis ( Figure 4 ). This highlights their potential as prognostic biomarkers and therapeutic targets for OS.

The bar charts ( Figure 5A ) show the total alteration frequency of six hub genes (CD36, CXCL12, CXCR4, CYBB, ITGAM, PLEK) in pan-cancer. The most frequently altered genes were CD36 (7.88%) and CYBB (8.1%), mainly due to copy number amplification. Moderate alteration frequencies were observed in ITGAM and PLEK (3.63% and 2.53%), whereas lower frequencies were observed in CXCR4 and CXCL12 (1.87% and 2.14%). These findings suggest that variation in high-frequency genes is dominated by copy number amplification and that their overexpression may drive tumor-related functions. Bar charts ( Figure 5B ) showed the frequency and distribution of genetic alterations of these six genes in different cancer types. These hub genes showed different patterns of alterations in different tumors, suggesting their potential roles in a variety of cancer types except OS. Notably, CYBB showed a high frequency of deletions in various malignancies such as liposarcoma, bladder squamous cell carcinoma and severe ovarian cancer.

According to TIIC analysis ( Figure 6 ), all six hub genes were negatively correlated with tumor purity (p < 0.001), which may indicate that these genes are closely associated with non-tumor cells like immune cells within TME. Among them, CYBB, ITGAM, and PLEK showed significantly negative correlations with tumor purity (ρ = -0.555, -0.514, and -0.44), suggesting their potential as markers of the immune microenvironment. CD36 exhibited a positive association with macrophages (ρ = 0.329) and neutrophils (ρ = 0.354), which may correspond to its involvement in lipid metabolism and immune modulation (25).

CXCL12 and CXCR4 are positively correlated with neutrophils (ρ = 0.129, 0.362), which may be related to their role in promoting neutrophil recruitment through chemotaxis (26). This indicates that the CXCL12-CXCR4 axis plays a role in promoting immune cell migration (27). The CXCL12-CXCR4 axis can recruit immune cells such as T cells and dendritic cells to TME, thereby enhancing immune surveillance and influencing tumor progression through immune-mediated pathways, consistent with the results of TIIC analysis in this study (28–31).

CYBB showed a positive correlation with macrophages (ρ = 0.778) and dendritic cells (ρ = 0.696), highlighting its central role in reactive oxygen species (ROS) generation in myeloid cells (32). It also showed the strongest negative correlation with tumor purity (ρ = -0.555), reflecting its potential central role in TME and immune regulation (33, 34).

ITGAM was strongly correlated with macrophages (ρ = 0.752), dendritic cells (ρ = 0.735), and CD4+ T-cells (ρ = 0.579), suggesting its involvement in immune cell adhesion and antigen immune processes of the presenting cells (35, 36).

Finally, PLEK was positively correlated with all immune cell types, highlighting its potential role in shaping the tumor immune landscape. Functionally, PLEK regulates actin cytoskeleton remodeling via integrins and Rac GTPases downstream of PKC. It is involved in leukocyte activation and integrin-mediated signaling, suggesting a role in macrophage recruitment and polarization. Its expression aligns with enhanced phagocytic and migratory activity in macrophages, neutrophils, and Tregs, supporting its contribution to immune cell activation and infiltration (37–39).

We validated the expression of six hub genes (CXCR4, CYBB, ITGAM, PLEK, CD36, and CXCL12) in OS and normal tissues respectively. The results of qRT-PCR showed that the expression of CYBB and PLEK was up-regulated in OS tissues compared with normal tissues, whereas there was no significant trend in the expression of CXCR4, ITGAM, CD36 and CXCL12 in both tissue types (P > 0.05) ( Figure 7A ).Then, we used Western blot analysis to detect the differences in protein expression levels of these two key genes in OS tissue and normal tissue. We found that the expression level of PLEK protein in OS tissue was higher than that in normal tissue, which was consistent with the PCR results, while CYBB protein expression showed no significant trend ( Figure 7B ).

Next, different cell populations were identified by single-cell RNA sequencing (scRNA-seq) in OS tissues, including cancer-associated fibroblasts (CAFs), OS cells, NK/T cells, macrophages, plasma cells, endothelial cells, B cells, and mast cells ( Figure 7C ). The expression of hub genes for each cell population was shown in Figure 7D . Violin plots showed that PLEK is predominantly expressed in macrophage subsets, with minimal expression in non-immune compartments, suggesting that this gene may be involved in the immune microenvironment of tumors, particularly in relation to immune function, immunosuppression, or pro-tumorigenic effects of macrophages (40) ( Figure 7E ).

Analysis of immune cell distribution in the six groups of OS tissues revealed a high percentage of macrophages in the tumor microenvironment, which was consistent with the identified hub genes expression patterns ( Figure 7F ). In addition, metabolic pathway enrichment analysis revealed that key pathways like glycolysis, oxidative phosphorylation and lipid metabolism were observed upregulated in macrophages ( Figure 7G ). And this metabolic reprogramming is also an important feature of cancer (41), for instance, reversing glycolysis to OXPHOS in cancer cells has been shown to induce cell death (42).

Among the nine hub genes identified through integrative bioinformatic analysis, PLEK was strategically prioritized for in-depth investigation based on a combination of biological and translational criteria. First, among all candidates, PLEK was the only gene consistently upregulated in osteosarcoma tissues at both the mRNA and protein levels, as confirmed by qRT-PCR and Western blotting. Second, single-cell RNA sequencing revealed that PLEK exhibited the highest cell-type specificity within macrophage subsets in the OS microenvironment, whereas other hub genes, such as CXCR4, were more broadly expressed across diverse immune cell populations. Third, PLEK showed one of the strongest associations with overall survival in OS patients, underscoring its prognostic relevance. Fourth, PLEK was enriched within metabolically active macrophage clusters characterized by elevated glycolytic and oxidative phosphorylation signatures, linking it directly to the metabolic axis of the tumor microenvironment—a key focus of this study. Collectively, these features place PLEK at the intersection of immune specificity, prognostic significance, and metabolic regulation, providing a robust rationale for its selection as a candidate for mechanistic exploration.

To characterize the cellular communication network within the TME of OS, we performed a systematic ligand-receptor interaction analysis across OS cells and various stromal and immune populations. The global interaction landscape revealed that macrophages exhibited the highest number and strength of interactions with other cell types, most notably with OS cells ( Figure 8A ). This suggests that macrophages serve as central communication hubs within the TME and may play a pivotal role in shaping tumor behavior.

To dissect the signaling specificity driving these interactions, we next examined the outgoing signaling patterns of each cell type across key paracrine pathways ( Figure 8B ). Macrophages and OS cells emerged as dominant sources of multiple tumor-relevant signals, including SPP1, TNF, MIF, GALECTIN, ANNEXIN, and VEGF. These outgoing signals showed pathway-dependent specificity, reflecting diverse functional roles in modulating immune, stromal, and vascular responses. A more detailed analysis of pathway-specific interactions revealed distinct yet interconnected roles for each signaling axis ( Figure 8C ).

In the SPP1 pathway, macrophage-derived osteopontin (OPN) targeted OS cells, endothelial cells, and NK/T cells. This interaction aligns with OPN’s established roles in extracellular matrix remodeling, immune modulation, and metastatic niche formation. Encoded by SPP1, OPN reshapes the tumor immune microenvironment by inducing M2-like TAM polarization, suppressing dendritic cell migration, and promoting Th1 responses via IL-12 and TNF-α. OPN⁺ TAMs engage FAP⁺ CAFs through IL-1 and TGF-β signaling to restrict T cell infiltration and remodel the ECM, thereby limiting CTL access. In parallel, OPN enhances tumor growth, metastasis, and angiogenesis by activating survival and autophagy pathways, inducing EMT, and upregulating VEGF via integrin signaling (43).

TNF signaling, primarily from macrophages to OS cells and CAFs, points to a central role for inflammation in tumor-stroma interaction and possibly in stromal activation and immune suppression. TNFα and IL-1 inhibit osteogenic differentiation of OS cells by activating the ERK signaling pathway, thereby maintaining an undifferentiated state that facilitates tumor progression and metastasis. In the tumor microenvironment, macrophages are the primary source of TNFα (44).

MIF signaling was particularly robust between OS cells and macrophages, indicating a feedback loop wherein tumor-derived and macrophage-derived MIF mutually reinforce inflammatory signaling, angiogenesis, and macrophage recruitment. MIF is highly expressed in OS and promotes tumor proliferation and metastasis by activating the RAS/MAPK pathway, and is associated with poor prognosis. Targeting MIF can inhibit OS progression, enhance chemotherapy sensitivity, and exert anticancer effects by regulating the tumor microenvironment (45).

The GALECTIN pathway mediates bidirectional crosstalk among OS cells, macrophages, and CAFs. Galectins, a conserved family of glycan-binding proteins, regulate tumor progression by integrating intracellular and extracellular signals in both cancer and stromal compartments (46). Galectin-1 serves as a diagnostic marker distinguishing chondroblastic OS from conventional chondrosarcoma (47), and has also been linked to OS progression and metastasis (48). Galectin-3 promotes OS malignancy via a triple mechanism: driving bone resorption through RANKL/M-CSF signaling, enhancing metastasis via FAK/Src/β-catenin activation, and conferring chemoresistance by inhibiting apoptosis via its NWGR motif (49). Galectin-9, through interaction with Tim-3, modulates tumor immunity; blockade of the Gal-9/Tim-3 axis elicits anti-tumor immune responses and suppresses tumor growth (50).

In the ANNEXIN pathway, macrophages initiate key signals to OS cells, CAFs, and endothelial cells, highlighting annexin-mediated roles in inflammation resolution and membrane repair within the tumor microenvironment. AnxA1, abundantly expressed in macrophages, suppresses inflammation via formyl peptide receptors (FPRs), shaping the immune milieu. In OS cells, enhanced AnxA2 translocation to the plasma membrane releases TFEB from the AnxA2-TFEB complex, promoting autophagy and regulating differentiation. Additionally, AnxA2 enhances membrane repair and facilitates OS metastasis. Elevated ANXA3 expression further promotes OS cell proliferation, migration, and invasion (51–53).

Interestingly, VEGF signaling originated mainly from OS cells and was directed toward endothelial and immune cells, aligning with its classical role in promoting angiogenesis and vascular remodeling. Angiogenesis critically influences tumor growth and metastatic potential, and is recognized as one of the six hallmarks of cancer (54). Vascular endothelial growth factor (VEGF), a central mediator of angiogenesis and vasculogenesis, plays a pivotal role in the molecular pathogenesis of tumor progression and metastasis (54). Targeting VEGF signaling has been shown to suppress osteosarcoma cell proliferation and induce apoptosis (55).

Collectively, these findings reveal a complex, multi-pathway signaling network orchestrated by OS cells and macrophages, characterized by both unidirectional and reciprocal communication. This network underpins key hallmarks of OS progression, including chronic inflammation, stromal activation, immunomodulation, and angiogenesis. Macrophages, in particular, emerge as key amplifiers and mediators of tumor-promoting signals, highlighting their potential as targets for therapeutic intervention.

Based on previous studies demonstrating extensive intercellular signaling between macrophages and OS cells, particularly through pro-tumor signaling pathways such as MIF, SPP1, and TNF, we aimed to functionally assess the contribution of the cytoskeletal adaptor protein PLEK (associated with immune cell signaling) to the malignant behavior of OS cells. To this end, we established a Transwell-based co-culture system to mimic the interactions between macrophages and OS cells under paracrine conditions in TME. OS cells were co-cultured with macrophages transfected with control siRNA (NC) or PLEK-specific siRNA (si-PLEK), as shown in Figure 9A . This allowed us to conduct a controlled assessment of PLEK function in a tumor-supportive cytokine environment. First, we confirmed the successful knockdown of PLEK by transfecting macrophages with three independent siRNAs. qPCR results showed that PLEK mRNA expression levels were significantly downregulated in the siRNA groups, particularly in the si2 and si3 groups ( Figure 9B ), and corresponding PLEK protein expression levels were also significantly reduced ( Figure 9C ).

To assess cell proliferation, we performed a 96-hour CCK-8 assay. Compared with the NC group, the survival rate of OS cells in the si-PLEK group was significantly higher, with differences becoming evident after 48 hours and further amplified at 72 and 96 hours ( Figure 9D ). Subsequently, we assessed the clonogenic capacity of OS cells via colony formation assays. The number and volume of colonies formed by si-PLEK group cells were higher than those in the NC group ( Figure 9E ), consistent with the enhanced long-term proliferative capacity. These results were validated by EdU uptake experiments, showing that the percentage of EdU-positive (proliferating) nuclei in PLEK-silenced OS cells was significantly higher ( Figure 9F ), indicating more active DNA synthesis and cell cycle progression. To investigate whether PLEK influences the migration and invasion capabilities of OS cells (key determinants of metastatic potential), we conducted Transwell migration assays. PLEK knockdown significantly promoted cell migration and invasion ( Figure 9G ), with a significant increase in the number of cells that migrated through the membrane in the si-PLEK group. These results collectively indicate that PLEK inhibits the malignant transformation of OS cells by regulating multiple aspects, including proliferation, clonogenic survival, and invasiveness. Mechanistically, PLEK may exert its effects through macrophage-derived cytokine signaling pathways, such as the MIF and TNF pathways previously identified as active in macrophage-osteosarcoma cell communication. As a cytoskeletal regulator, PLEK may integrate upstream inflammatory signals to promote cytoskeletal remodeling and cell migration. Importantly, this study provides functional evidence that disrupting PLEK expression in OS cells enhances the tumor-promoting effects of macrophage-derived paracrine signals. These findings suggest that PLEK may represent a potential therapeutic target, potentially linking signals from the immune microenvironment to the malignant progression of OS.