Computational discovery of novel colony‑stimulating factor‑1 receptor as a potential therapeutic biomarker in osteosarcoma and a novel inhibitor from herbal sources

Introduction

Osteosarcoma is a highly aggressive primary bone malignancy that predominantly impacts pediatric, adolescent and young adult populations, with peak incidence occurring in individuals aged between 10–25 years. Epidemiologically, osteosarcoma constitutes nearly 20% of global primary bone tumors and is among the most common primary bone malignancies affecting pediatric and young adult populations (1,2). The current standard of care for integrating surgical resection with chemotherapy and radiotherapy (3) fails to achieve optimal therapeutic outcomes; 5-year survival rates have remained stagnant at 65–70% for decades, and survival outcomes are even lower (≤30%) in patients with metastatic or recurrent disease (4,5). Furthermore, existing clinical biomarkers for osteosarcoma [such as alkaline phosphatase (ALP), lactate dehydrogenase and tumor necrosis factor-related apoptosis-inducing ligand] possess notable limitations: i) Markers lack sufficient specificity, for example elevated ALP, is also common in benign bone disorders such as rickets; ii) markers exhibit low sensitivity for early-stage tumors, often remaining normal until disease progression; (iii) and markers cannot effectively stratify patients into prognostic subgroups or predict responses to specific therapies (6–8). As a result, both preoperative prognostic stratification and postoperative disease monitoring remain clinically challenging, highlighting the urgent need to identify novel molecular biomarkers that can improve early diagnostic accuracy, enable precise prognostic assessment and guide the development of targeted therapeutic strategies.

Genes that regulate osteoclast differentiation and activation are recognized as promising therapeutic targets for bone resorption-related disorders, including osteosarcoma (9–11). Among these receptors, the colony-stimulating factor-1 receptor (CSF1R), a key receptor tyrosine kinase, serves a pivotal role in mediating the differentiation, proliferation and functional homeostasis of tissue-resident macrophages and osteoclasts (12–17). Dysregulated CSF1R signaling has been linked to a broad spectrum of pathologies, including spanning bone metabolic disorders (such as osteoporosis), hematological and solid malignancies (such as gastric cancer, ovarian cancer and sarcoma), and neurodegenerative diseases (12,15). Mechanistically, binding of its cognate ligands (CSF1 or IL-34) to the extracellular domain (ectodomain) of CSF1R induces conformational changes in the receptor, triggering autophosphorylation of specific tyrosine residues within its intracellular kinase domain. This phosphorylation event subsequently activates a cascade of downstream pro-survival and pro-proliferation signaling pathways, including PI3K/Akt, ERK1/2 and JNK (16–18). Immunologically, CSF1R exerts critical regulatory effects on tumor-associated macrophages (TAMs), which are major components of the tumor microenvironment (TME). Specifically, CSF1R signaling promotes TAM recruitment and polarization towards an M2-like (protumoral) phenotype; these TAMs then secrete an array of tumor-promoting cytokines (such as IL-6 and TGF-β) and growth factors, which facilitate tumor cell proliferation, angiogenesis and immune evasion (19).

Given its multifaceted roles in tumor progression and bone homeostasis, CSF1R has emerged as a potential therapeutic target for antitumor and anti-bone-resorption strategies (20–22). However, the expression pattern, functional relevance and potential diagnostic or prognostic value of CSF1R in osteosarcoma remain largely uncharacterized and thus require further comprehensive investigation. The present study aimed to screen for key genes with potential regulatory roles in osteosarcoma progression, by systematic analysis of the expression profiles of osteoclast differentiation-related genes (ODRGs) using multiple Gene Expression Omnibus (GEO) datasets. Focusing on CSF1R, a candidate gene identified from the ODRG panel, a comprehensive bioinformatics analyses was then performed using the TARGET-osteosarcoma dataset to evaluate its clinical and immunological relevance. Specifically, the expression patterns of CSF1R across different clinical subgroups were characterized, CSF1R correlations with immune checkpoint molecules were analyzed (such as programmed death-ligand 1 and CTLA-4) and prognostic value was validated in patients with osteosarcoma. Furthermore, to clarify the cellular source of CSF1R in the TME, its cellular distribution was characterized using single-cell RNA sequencing (scRNA-seq) data from osteosarcoma tissues; the analysis was also extended to explore the pan-cancer expression landscape of CSF1R across multiple tumor types from The Cancer Genome Atlas (TCGA) database, providing a broader perspective on its role in oncogenesis. Finally, functional enrichment analyses were conducted to predict CSF1R-related signaling pathways involved in osteosarcoma, and this was supplemented with in vitro experimental validation to comprehensively investigate the potential of CSF1R as both a diagnostic biomarker and therapeutic target for osteosarcoma.

Owing to their structural diversity and favorable safety profiles, natural products, particularly those isolated from traditional medicinal plants, have long served as cornerstones for anticancer drug discovery (23). In the present study, a rational drug discovery pipeline was adopted to identify novel CSF1R inhibitors from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (https://old.tcmsp-e.com/tcmsp.php), a comprehensive repository of bioactive compounds derived from herbal medicines. Leveraging the high-resolution crystal structure of CSF1R, structure-based molecular docking simulations were performed using AutoDock Vina, a widely used tool for virtual ligand screening. The docking workflow was executed in two key steps: i) Computational predictions and validation of the ligand-binding sites of CSF1R (focusing on the intracellular kinase domain, the primary target for small-molecule inhibitors) to define high-confidence interaction pockets; and ii) virtual screening of the TCMSP compound library against these predicted binding pockets, ranking candidates based on docking scores (such as binding affinity values) to prioritize potential CSF1R inhibitors. These findings collectively validated the potential of herbal compound as a promising lead compound for CSF1R-targeted osteosarcoma therapy.

Materials and methods

Data sources, preprocessing and functional enrichment analysis

Osteosarcoma datasets GSE33382 (24), GSE14359 (25) and GSE218035 (26) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Gene expression and clinical data from patients with osteosarcoma was obtained from the TARGET-Osteosarcoma database (https://www.cancer.gov/ccg/research/genome-sequencing/target/studied-cancers/osteosarcoma), a publicly available resource for pediatric cancer, including microarray data for 89 patients and single-cell RNA-sequencing data for 98 patients. The expression distribution of CSF1R in pan-cancer cell lines was explored from the CCLE database (https://portals.broadinstitute.org/ccle).

Differentially expressed genes (DEGs) between normal and osteosarcoma tissues were identified using the R package ‘limma’ (25) (version 3.40.6) across the three GEO datasets, with thresholds set at |log2FC| >2 and P<0.05. Subsequently, the enriched DEGs were subjected to pathway analysis using the WikiPathways (https://www.wikipathways.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (https://www.kegg.jp).

Consensus clustering analysis for ODRGs

ODRGs, including acid phosphatase 5, tartrate resistant (ACP5), CSF1R, CSF1, CTSK, four and a half LIM domains 2 (FHL2), FOS-like 1, AP-1 transcription factor subunit (FOSL1), LILRB4, ITGB3, TYROBP and TNFRSF11B, were identified. The ODRGs identified from the enrichment analysis were clustered and scored based on their expression profiles using the ‘ConsensusClusterPlus’ R package (version 1.50.0) for consistency analysis. The ‘survival’ R package (26) (version 4.1) was utilized to perform survival analysis, specifically to compare overall survival (OS) differences between distinct patient clusters. To visualize disparities in clinical characteristics (such as age, tumor stage and pathological grade) across these clusters, a heatmap was constructed using the ‘pheatmap’ R package (26) (version 1.0.13), with rows representing clinical features and columns representing patient clusters for intuitive comparison.

Survival analysis, expression and immune checkpoint sensitivity analysis of CSF1R

Osteosarcoma patient samples were divided into two groups in equal numbers. For Kaplan-Meier (KM) survival curves, statistical significance (P-values) and hazard ratios (HR) with 95% confidence intervals were computed using the log-rank test and univariate Cox proportional hazards regression, respectively. Additionally, CSF1R expression patterns and immune checkpoint sensitivity across the TARGET and GEO datasets were analyzed.

Single-cell transcriptome analysis

The osteosarcoma single-cell RNA sequencing (scRNA-seq) dataset (GSE162454) was acquired from the Tumor Immune Single-Cell Hub 2 database (http://tisch.comp-genomics.org/) by implementing a text-mining-based data parsing workflow (27). Prior to downstream analysis, raw single-cell data were subjected to quality control filtering with the following criteria: i) Genes expressed in ≥3 cells were retained; ii) cells with ≥200 detected genes per cell were included; iii) the number of unique molecular identifiers (UMIs) per cell were filtered to the range of 500–6,500, adjusted according to the UMI distribution of the dataset; and iv) only cells with mitochondrial read percentages <80% were kept to exclude severely degraded cells.

The filtered single-cell data, comprising UMI transcript matrices (recording gene expression levels) and cellular metadata matrices (containing cell-level information), were merged into Seurat objects, the standard input format for scRNA-seq analysis. Data normalization was performed using the LogNormalize method (scaling gene expression by total UMI counts per cell, followed by log transformation). To assess potential batch effects between samples, uniform manifold approximation and projection dimensionality reduction was applied; no significant batch effects were observed across samples, confirming the reliability of subsequent analyses.

Comprehensive bioinformatics analysis of CSF1R in pan-cancer datasets

The expression levels of CSF1R were examined in 33 types of cancer from the TCGA cohort via Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn; version 3.0) and Genotype-Tissue Expression (GTEx; https://www.genome.gov/Funded-Programs-Projects/Genotype-Tissue-Expression-Project; version 8.0) analysis. Cox regression analysis was carried out to investigate if CSF1R could be an independent predictive factor in patients with multiple types of cancer. Lastly, the relationship between CSF1R and the invasion of 22 different types of immune cells in pan-cancer was evaluated using the CIBERSORT algorithm (https://cibersortx.stanford.edu; version 1.03).

Molecular docking and molecular dynamics (MD) simulation studies

Molecular docking was used to study the interactions between the compounds and the proteins. The 3D structure of sotuletinib was retrieved from PubChem (https://pubmed.ncbi.nlm.nih.gov/?term=sotuletinib+&filter=datesearch.y_5), while the protein structures of CSF1R were obtained from the RCSB Protein Data Bank (PDB) database (PDB no. 6WXJ; http://www.rcsb.org/structure/6WXJ). Molecular docking was performed using the AutoDock Vina (https://vina.scripps.edu; version 1.1.2) program, with binding site identification preprocessed via AutoDock Tools. Ligands were subjected to energy minimization using the Avogadro (https://avogadro.cc; version 2.0) program. To ensure docking reliability, the exhaustiveness parameter was set to 100. Following docking, putative hydrogen bonds were analyzed, and complex structures were visually inspected using PyMOL (https://www.pymol.org/pymol.html; version 3.1).

All docked complexes (CSF1R-sotuletinib and other candidate complexes) were subjected to 100 nsec all-atom MD simulations using the GROMACS (https://www.gromacs.org; version 2025) software. Partial charges for sotuletinib and CSF1R were generated using the Antechamber module in AmberTools23 (https://www.osc.edu/resources/technical_support/hpc_changelog/2023/ambertools_23_is_available; version 2023). The AMBER99SB-ILDN force field was employed to construct the protein topology. Each system was solvated in a dodecahedron unit cell using the TIP3P water model, and charge neutrality was achieved by adding sodium (Na+) and chloride (Cl-) counterions.

Prior to production simulations, energy minimization was conducted until the system reached a predefined force constant threshold to eliminate steric clashes. The minimized systems were then equilibrated in two stages: i) A constant particle number, volume and temperature equilibration at 300 K; and ii) a constant particle number, pressure and temperature equilibration at 1 atm pressure using the Parrinello-Rahman barostat, with a leap-frog integrator used for both stages.

For the production MD phase, a time step of 2 femtosec was applied, and trajectory frames were saved every 10 picosec. Post-simulation trajectory analysis was performed to evaluate system stability, including calculations of root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), solvent-accessible surface area (SASA) and hydrogen bond interactions.

In silico screening of herbal compounds

A structured in silico virtual screening workflow was established to identify potential CSF1R-targeting small molecules from natural herbal sources. First, the 3D structures of herbal compounds in SDF format were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). For ligand preparation prior to docking, OpenBabel (https://computing.ch.cam.ac.uk/software/openbabel; version 2.4.1) was used to process these bioactive compounds: Hydrogen atoms were added to simulate physiological pH (pH 7), and 3D structures were generated in PDB format. The resulting OpenBabel outputs served as the initial structures for AutoDockFR (https://ccsb.scripps.edu/adfr/; version 1.0), which converted them into PDBQT format (a required format for docking).

For protein preparation, the CSF1R structure (PDB ID. 6WXJ) was retrieved from the PDB. AutoDockFR was then used to prepare CSF1R for virtual screening and molecular docking: Gasteiger charges were assigned to the protein, and its structure was converted to PDBQT format. To define the binding region, a grid box (dimensions, 19.5×25.5×27.0 Å) was generated using AutoDockFR to fully enclose CSF1R's active binding pocket. Finally, AutoDock Vina was employed for docking simulations, using the prepared CSF1R PDBQT file and a library of 13,227 herbal compound ligands from TCMSP database (https://old.tcmsp-e.com/tcmsp.php).

Cell counting kit-8 (CCK-8) assay

Non-immortalized human umbilical vein endothelial cells (HUVECs) at passage 2 were purchased from the Cell Bank of the Chinese Academy of Sciences (cat. no. #PCS-01; http://www.cellbank.org.cn/search-detail.php?id=1375). MG63 osteosarcoma cells were purchased from Procell Life Science Technology Co., Ltd. (cat. no. #CL-0157; http://www.procell.com.cn/p/mg-63-cl-0157-72575). Saos-2 osteosarcoma cells were purchased from the Cell Bank of the Chinese Academy of Sciences (cat. no. #SCSP-5057; http://www.cellbank.org.cn/search-detail.php?id=729). HUVECs and osteosarcoma cells (MG63 and Saos-2) were cultured in RPMI-1640 (NEST®; Wuhan NEST Biotechnology Co., Ltd.) medium supplemented with 10% fetal bovine serum (NEST®; Wuhan NEST Biotechnology Co., Ltd.) under standard cell culture conditions (37°C and 5% CO2). Cell viability was evaluated using a CCK-8 assay (Solarbio®; Beijing Solarbio Life Sciences Co., Ltd.). Briefly, MG63 and Saos-2 cells were individually treated with gradient concentrations (0, 2, 4, 8, 20, and 40 µmol/l) of three test agents (pexidartinib, sotuletinib and sarsasapogenin) prior to being seeded into 96-well plates. Following 48 h of treatment, cells were further incubated for 40 min by adding CCK-8 reagent, and a multimode microplate reader was employed to measure the optical density of each well at a wavelength of 450 nm. Each experimental group was set up in triplicate to ensure reproducibility.

Western blotting

MG63 and Saos-2 cells were first washed twice with ice-cold phosphate-buffered saline to remove residual culture medium, then lysed in RIPA protein lysis buffer (cat. no. P0013VS; Beyotime Biotechnology Co., Ltd.) supplemented with a protease inhibitor cocktail (to prevent protein degradation). Total protein concentration in the cell lysates was quantified using a BCA Protein Assay Kit (cat. no. P0011; Beyotime Biotechnology), ensuring equal protein (30 ug) loading across samples. Equal amounts of protein were separated on 8% gels using sodium dodecyl sulfate-polyacrylamide gel electrophoresis, then transferred via electrophoresis onto a polyvinylidene difluoride membrane. The membrane was blocked with 5% non-fat milk in Tris-buffered saline with 0.05% Tween-20 (TBST) for 1 h at room temperature to reduce non-specific binding, followed by incubation with the indicated primary antibodies against CSF1R (1:500; cat. no. YT0881; ImmunoWay Biotechnology Company) and β-actin (1:5,000; cat. no 20536-1-AP; Proteintech Group, Inc.) at 4°C overnight. After three 10-min washes with TBST to remove unbound primary antibody, the membrane was incubated with a horseradish peroxidase-conjugated Goat Anti-Rabbit IgG (H+L) secondary antibody (1:20,000; cat. no. A21020; Abbkine Scientific Co., Ltd.) for 1 h at room temperature. Finally, protein bands were visualized using an ECL prime western blotting detection kit (cat. no. BMP3010; Abbkine Scientific Co., Ltd.). Protein band intensities were quantified using the ImageJ software (version 1.8.0; National Institutes of Health). Western blot analysis was performed in three independent biological replicates. For each experimental replicate, protein expression levels were normalized to those of the corresponding DMSO-treated control group. As a result, all control values were set to 1, with no variability shown for the control groups.

Statistical analysis

Bioinformatic analyses were conducted using the R software (version 4.0.3; Posit Software, PBC). The Kruskal-Wallis test was used for count sample. Wilcoxon rank-sum test was used to immune checkpoint sensitivity analysis. Data from in vitro experiments were analyzed using GraphPad Prism (version 9.0; Dotmatics). Correlations were evaluated using Spearman's rank correlation test. Survival risk and HR were derived from Cox proportional-hazards models. KM survival curves were compared using the log-rank test. Group differences were assessed by one-way ANOVA followed by Dunnett's multiple-comparison test. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of ODRGs in osteosarcoma

Differential gene expression analysis was conducted to compare gene profiles between normal bone tissues and osteosarcoma tissues across three independent GEO datasets, namely GSE33382, GSE14359 and GSE218035, using the limma R package, a widely used tool for microarray and RNA-seq data analysis.

From the GSE33382 dataset, 781 DEGs were identified, of which 399 were significantly upregulated and 382 were downregulated (Fig. 1A). From analysis of the GSE14359 dataset, a larger set of 2,312 DEGs was detected, with 1,299 genes showing increased expression and 1,013 showing decreased expression (Fig. 1B). Analysis of the GSE218035 dataset yielded 4,852 DEGs, including 2,389 upregulated and 2,463 downregulated genes (Fig. 1C). Venn diagram analysis of the DEGs from the three datasets was performed to identify genes whose expression levels were consistently dysregulated across the cohorts, which revealed 156 overlapping DEGs (Fig. 1D); these genes were considered core candidates for further functional exploration.

Figure 1. Identification of ODRGs in osteosarcoma. Volcano plots depicting the DEGs in the (A) GEO-GSE33382, (B) GEO-GSE14359 and (C) GEO-GSE218035 datasets (log2 FC >1; P<0.05). (D) A Venn diagram identifying 156 DEGs in the overlap of the three datasets retrieved from the GEO database. (E) KEGG enrichment analysis of the 156 DEGs. (F) WIKI pathways enrichment analysis of the 156 DEGs. (G) Analysis of expression levels of osteoclast differentiation genes (FHL2, ACP5, FOSL1, CSF1R and TNFRSF11B) in the GSE33382 datasets using an unpaired two-tailed Student's t-test. (H) The expression levels of osteoclast differentiation genes were assessed between osteosarcoma tissues and normal bone tissues (n=2) in the GSE14359 dataset. (I) Analysis of expression levels of osteoclast differentiation genes (FHL2, ACP5, FOSL1, CSF1R and TNFRSF11B) in the GSE218035 datasets using an unpaired two-tailed Student's t-test. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; CSF1R, colony-stimulating factor-1 receptor; FHL2, four and a half LIM domains 2; ACP5, acid phosphatase 5, tartrate resistant; FOSL1, FOS-like 1, AP-1 transcription factor subunit; TNFRSF11B, TNF receptor superfamily member 11b.

Functional enrichment analysis of 156 shared DEGs

Osteosarcoma cells secrete osteoclastogenic mediators, including RANKL and IL-6, to fuel osteoclast differentiation and bone resorption (9). This process releases tumor-promoting factors from the degraded bone matrix, establishing a self-perpetuating ‘tumor-osteoclast vicious cycle’ that exacerbates tumor invasion and metastasis (10). KEGG pathway analysis revealed significant enrichment in the ‘osteoclast differentiation’ pathway (adjusted P<0.01), which indicated that the core dysregulated genes were closely associated with osteoclast maturation processes (Fig. 1E). Wiki pathways analysis further confirmed association of the 156 DEGs with ‘osteoclast signaling’ pathways, which regulate key downstream events in osteoclast activation (Fig. 1F).

The expression patterns of five key ODRGs known to regulate the osteoclast differentiation pathway (6,28–31) were then analyzed, including FHL2, ACP5, FOSL1, CSF1R and TNFRSF11B, across the GEO datasets. FHL2, FOSL1 and TNFRSF11B expression levels were expressed increased in normal bone tissues compared with those in osteosarcoma tissues, whereas ACP5 and CSF1R were significantly upregulated in osteosarcoma tissues compared with those in normal tissues in the GSE33382 (P<0.001; Fig. 1G and I). In the GSE14359 dataset, the expression levels of FHL2 and FOSL1 were lower in osteosarcoma tissues compared with those in normal bone tissues (n=2), whereas ACP5 and CSF1R exhibited higher expression in osteosarcoma samples than in normal bone tissues. The expression of TNFRSF11B did not differ between the two groups (statistical comparison not possible as n=2; Fig. 1H). In the GSE218035 datasets, ACP5, FOSL1, CSF1R and TNFRSF11B were significantly upregulated in osteosarcoma tissues compared with those in normal tissues (P<0.05; Fig. 1I). Univariate Cox regression analysis was performed for FHL2, ACP5, FOSL1, CSF1R and TNFRSF11B expression levels and OS using the TARGET dataset. CSF1R expression levels were associated with the OS of patients with osteosarcoma based on the microarray data (P<0.05; Fig. S1A). KM survival analysis was then conducted, which demonstrated that patients with osteosarcoma with high expression of CSF1R had increased OS rates compared with that of the low-expression group (Fig. S1B). Time-dependent ROC analysis demonstrated that CSF1R had predictive value for 1-, 3- and 5-year OS with areas under the curve (AUC) values of 0.693, 0.616 and 0.620, respectively (Fig. S1C). CSF1R (P<0.01) and TNFRSF11B (P<0.05) were associated with OS in patients with osteosarcoma in the univariate Cox regression analysis of RNA-seq data (Fig. S1D). KM survival analysis revealed that patients with osteosarcoma with high expression of CSF1R had increased OS rates compared with the CSF1R low-expression group (Fig. S1E). The AUC values for the OS rate for CSF1R were 0.572, 0.610 and 0.596 for 1-, 3- and 5-year survival, respectively (Fig. S1F). By contrast, patients with low expression of TNFRSF11B had increased OS rates, compared with those with high expression, as determined by RNA-seq data (Fig. S1G). Time-dependent ROC analysis revealed that the TNFRSF11B had predictive value for 1-, 3- and 5-year OS rates with AUCs of 0.714, 0.637 and 0.689, respectively (Fig. S1H). These results highlighted the pronounced heterogeneity in the expression of osteoclast-related genes between healthy individuals and patients with osteosarcoma, implying that dysregulation of these genes may serve a critical role in driving osteosarcoma initiation and progression.

Immune checkpoint sensitivity analysis

To evaluate differences in immunotherapy sensitivity among patients with osteosarcoma, patients were stratified into two groups based on CSF1R expression levels, the high-expression (G1) and low-expression (G2) groups. Their immune checkpoint blockade (ICB) therapy response patterns were then compared using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm, a well-validated computational tool specifically designed to predict clinical responses to ICB. The analysis revealed significant differences in ICB sensitivity across the two groups; the G1 TIDE score was significantly lower compared with that of G2; a lower TIDE score indicates higher ICB responsiveness. This pattern was consistently validated across both the gene chip array (Fig. 2A) and RNA-seq (Fig. 2B) cohorts, which ruled out cohort-specific biases.

Figure 2. Assessment of the clinical response to immune checkpoint inhibitor therapy. (A) In the gene chip array cohort of osteosarcoma, differences in TIDE algorithm scores of CSF1R between the high- and low-expression groups in the TARGET dataset were assessed (Wilcoxon rank sum test). (B) In the RNA-seq data of osteosarcoma, differences in TIDE algorithm scores of CSF1R between the high- and low-expression groups in the TARGET dataset were assessed (Wilcoxon rank sum test). (C) The sensitivity of immune checkpoint between the high- and low-expression groups of CSF1R in the TARGET dataset (gene chip array data in osteosarcoma). (D) The sensitivity of immune checkpoint between the high- and low-expression groups of CSF1R in the TARGET dataset (RNA-seq data in osteosarcoma). **P<0.01 and ***P<0.001. ns, not significant; CSF1R, colony-stimulating factor-1 receptor; G1, high-expression group; G2, low-expression group; RNA-seq, Ribonucleic Acid-sequencing; TARGET, Therapeutically Applicable Research To Generate Effective Treatments; TIDE, Tumor Immune Dysfunction and Exclusion.

To further examine the mechanism underlying this difference, the expression of key immune checkpoint molecules were analyzed in the two groups: In the gene chip array cohort, compared with G2, G1 showed significantly increased expression levels of immune checkpoint genes, including hepatitis A virus cellular receptor 2 (HAVCR2), programmed cell death 1 (PDCD1), PDCD1 ligand 2 (PDCD1LG2) and sialic acid-binding ig-like lectin-15 (SIGLEC15) (P<0.01; Fig. 2C). RNA-seq data corroborated these findings; G1 showed significantly higher expression of the immune checkpoints HAVCR2, lymphocyte activation gene-3 and SIGLEC15 than G2 (P<0.001; Fig. 2D).

Notably, high expression levels of these immune checkpoints are often associated with a ‘hotter’ tumor immune microenvironment, one that is more responsive to ICB therapy (32). Collectively, the present results suggested that elevated CSF1R expression in osteosarcoma may be associated with enhanced sensitivity to immunotherapy.

Single-cell sequencing analysis of CSF1R

Using the human osteosarcoma scRNA-seq dataset GSE162454, the expression pattern of CSF1R within the osteosarcoma TME was systematically characterized. First, cell type annotation for the dataset was performed, which identified eight major cell populations: CD4+ conventional T cells (CD4Tconx), endothelial cells, cancer-associated fibroblasts, CD8+ exhausted T cells, malignant osteosarcoma cells, monocytes/macrophages, osteoblasts and plasmocytes (Fig. 3A). This classification provided a clear framework for mapping CSF1R expression across distinct TME components.

Figure 3. Single-cell sequencing analysis of CSF1R. (A) UMAP plot showing annotated cell categories after dimension reduction in the GSE162454 dataset. (B) UMAP plot showing the expression distribution of CSF1R in different cell categories. (C) Violin plot showing the results of CSF1R expression distribution in different cell categories from UMAP data analysis. (D) Violin plot representing relative expression levels of CSF1R for different cell types in 6 patients with osteosarcoma. CD4Tconx, CD4+ conventional T cells; CD8Tex, CD8+ exhausted T cells; mono/macro, monocytes/macrophages; CSF1R, colony-stimulating factor 1 receptor; N.S., not significant.

Subsequent expression analysis revealed that CSF1R exhibited striking cell type-specific expression, with notably high transcript abundance in two cell populations: Osteoblasts and monocytes/macrophages (Fig. 3B). Violin plots further visualized this differential expression, which showed that CSF1R expression was restricted to these two cell types and barely detectable in others, such as CD4Tconx and plasmocytes (Fig. 3C).

To validate this pattern across patients, CSF1R expression was analyzed in six osteosarcoma cases (GSE162454 dataset). Statistical analysis confirmed that CSF1R expression was significantly elevated in osteoblasts and monocytes/macrophages compared with all other cell populations (P<0.001; Fig. 3D), ruling out patient-specific variability and corroborating the consistent cell-type enrichment of CSF1R.

Pan-cancer analysis of CSF1R expression

To systematically characterize the pan-cancer expression levels and functional relevance of CSF1R, CSF1R expression patterns were analyzed across 33 tumor types using data from the TCGA and GEPIA2 databases, revealing distinct tumor-specific dysregulation patterns (Fig. S2A and B). Specifically, CSF1R was consistently upregulated in eight malignancies: Diffuse large B-cell lymphoma (DLBC), glioblastoma multiforme (GBM), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), pancreatic adenocarcinoma (PAAD) and testicular germ cell tumors (TGCT). By contrast, CSF1R expression was significantly downregulated in adrenocortical carcinoma (ACC).

Spearman's correlation analysis was next performed to explore the associations between CSF1R expression and key genomic instability markers. With respect to microsatellite instability (a marker of DNA mismatch repair defects), CSF1R expression was correlated in 18 types of cancer: Thymoma, colon adenocarcinoma (COAD), sarcoma (SARC), lung squamous cell carcinoma (LUSC), kidney chromophobe (KICH), breast invasive carcinoma (BRCA), lung adenocarcinoma (LUAD), pheochromocytoma and paraganglioma (PCPG), stomach adenocarcinoma (STAD), head and neck cancer (HNSC), uveal melanoma (UVM), ACC, PAAD, TGCT, liver hepatocellular carcinoma (LIHC), GBM, esophageal carcinoma (ESCA) and mesothelioma (Fig. S2C).

With respect to tumor mutational burden (a marker of somatic mutation accumulation), CSF1R expression was correlated with 18 malignancies, including COAD, SARC, KIRC, UVM, LAML DLBC, cholangiocarcinoma (CHOL), STAD, uterine carcinosarcoma, TGCT, LUSC, ACC, LUAD, ESCA, PAAD, KIRC, LIHC and THYM (Fig. S2D).

To further evaluate the clinical and immunological significance of CSF1R, univariate Cox regression identified CSF1R as a significant prognostic factor in skin cutaneous melanoma (P<0.0001; Fig. S3A), which indicated its potential value in predicting patient outcomes for this cancer type.

Immune checkpoint correlation was assessed, and CSF1R expression was positively correlated with the expression levels of most immune checkpoint molecules, with the most prominent correlations observed with HAVCR2 and PDCD1LG2 (P<0.01; Fig. S3B). These findings suggested that CSF1R expression may serve as an indirect indicator of immunotherapy response, given the role of these checkpoints in regulating ICB efficacy.

Immune infiltration analysis was performed using CIBERSORT to quantify immune cell subsets, and it was found that CSF1R expression was significantly positively correlated with macrophage populations (both M1 and M2 phenotypes) but negatively correlated with naïve B cells, plasma B cells and resting natural killer cells across various types of cancer (P<0.01; Fig. S3C). In SARC (a tumor type closely related to osteosarcoma), CSF1R expression was negatively correlated with M0 macrophages, activated mast cells and CD4+ naïve T cells, but strongly positively correlated with M1/M2 macrophages (P<0.001).

Collectively, this pan-cancer analysis underscored the conserved role of CSF1R in modulating the tumor immune microenvironment, particularly through regulating macrophage infiltration, and highlight its context-dependent relevance to genomic instability and patient prognosis.

Sensitization of the MG63 and Saos-2 cell lines to a CSF1R inhibitor

To evaluate the relationship of CSF1R inhibitor with CSF1R, small-molecule drug molecular docking was conducted. It was demonstrated that at a distance of 2.1 Å, sotuletinib established one hydrogen bond with the amino acid site SER807 in CSF1R (Fig. S4). Sotuletinib interacted with the CSF1R protein with a binding energy of −9.5 kcal/mol, which indicated the potential effectiveness of sotuletinib for osteosarcoma treatment due to its strong binding affinity for the target.

MD simulation investigations were executed to ascertain the stability and convergence of the CSF1R-sotuletinib complexes. These simulations lasted as long as 100 nsec, revealing a persistent and stable conformation, which was marked in the RMSD analysis. Specifically, regarding the stability of CSF1R-sotuletinib, deviations could still be identified up to ~60 nsec, after which it reached a higher stability level until the end of the 80 nsec period (Fig. S5A). The RMSD for CSF1R-sotuletinib was found to be 0.4 Å. The structural flexibility of CSF1R-sotuletinib was interpreted using the RMSF. The RMSF value of the CSF1R-sotuletinib complex was <1 Å, indicating that it had good stability (Fig. S5B). The Rg was a metric denoting protein compactness. The Rg for CSF1R-sotuletinib was 2.01 Å (Fig. S5C). The observed Rg of protein-ligand complexes was in the range of 1–3 Å, which is consistent with established acceptable standards (33). Two hydrogen bonds persisted throughout the entire 100 nsec simulation, indicating that the interaction of sotuletinib with CSF1R was both marked and lasting (Fig. S5D). Energy analysis of all amino acids of sototinib with CSF1R is shown in the Fig. S5E. The results of energy degradation revealed that sotuletinib strongly interacted with the key amino acids ASN808, LYS820 and TR821 of the CSF1R protein (Fig. S5F), which were thus identified as the key residues of the CSF1R-sotuletinib association.

To functionally validate the therapeutic potential of targeting CSF1R in osteosarcoma, the effects of two well-characterized CSF1R inhibitors, pexidartinib (PLX3397) and sotuletinib (22), on osteosarcoma cell viability were assessed. HUVECs (used as a normal cell control) and two osteosarcoma cell lines, MG63 and Saos-2 (both with low endogenous CSF1R expression), were treated with a gradient of inhibitor concentrations for 48 h. Cell viability was quantified using the CCK-8 assay, a standard method for evaluating drug-induced cytotoxicity.

Notable differences in inhibitor potency were demonstrated: Compared with sotuletinib, pexidartinib exhibited superior cytotoxic activity, with lower half-maximal inhibitory concentration (IC50) values in both MG63 and Saos-2 cells. Notably, compared with normal HUVECs, the CSF1R-low osteosarcoma cells (MG63 and Saos-2) displayed greater sensitivity to both inhibitors, as evidenced by steeper viability reduction curves and lower IC50 values in the MG63 and Saos-2 cancer cells (Fig. 4). These findings confirmed that the CSF1R inhibitors exhibited selective antitumor efficacy in osteosarcoma cell models, with minimal toxicity to normal vascular endothelial cells, supporting the potential of CSF1R as a safe and actionable therapeutic target for osteosarcoma.

Figure 4. Cell viability (%) of MG63 and Saos-2 cells treated with pexidartinib (via PLX3397) and sotuletinib.

Identification of sarsasopogenin as an inhibitor of CSF1R via virtual screening

To identify novel CSF1R inhibitors with high binding affinity, virtual screening of the TCMSP database, which contains 13,227 herbal-derived compounds, was performed. Using a binding score threshold of <-7 kcal/mol (which indicated significant affinity), 6,840 candidate compounds were initially identified. Among these, three compounds, linarin, saikosaponin C and sarsasapogenin, exhibited particularly strong interactions with CSF1R, with binding energies of −10.2, −10.0 and −10.6 kcal/mol, respectively. The molecular docking interactions between CSF1R and small molecules are shown in Fig. S6. Analysis of protein-ligand interactions revealed that the docked complex of CSF1R-linarin displayed 10 hydrophobic interactions and 17 non-bonded interactions. Additionally, the ligand formed three hydrogen bonds with amino acid residues GLU633, VAL647 and CYS774 (Fig. S6A) within the binding pocket. Furthermore, linarin interacted with the CSF1R amino acid residues GLY669, TYR665, THR633, GLU664, LEU649, ILE794, ILE646, ASN648, LEU640, LEU769, HIS776, ILE775, ILE636, ASP796, VAL661, LEU588, MET637 and PHE797 via Van der Waals contacts. Analysis of protein-ligand interactions revealed that the docked complex of CSF1R-saikosaponin c displayed 8 hydrophobic interactions and 10 non-bonded interactions. The findings showed that at a distance of 2.0 Å, saikosaponin c formed one hydrogen bond with the amino acid site ILE794, and bonded seven amino acid residues LEU640, ILE646, LEU769, HIS776, MET637, ILE636 and CYS774 of the CSF1R protein through hydrophobic interactions (Fig. S6B).

The effects of these three candidate drugs, linarin, saikosaponin C and sarsasapogenin, on osteosarcoma cell proliferation, were next evaluated using CCK-8 assays in MG63 and Saos-2 cells (Fig. 5A-5C). Notably, sarsasapogenin treatment resulted in the most pronounced reduction in cell viability in both cell lines (Fig. 5C). To prioritize candidates for further study, virtual screening results with drug-likeness predictions using Lipinski's and Veber's rules were integrated (Swiss ADME; http://swissadme.ch), which are pharmacological criteria that assess a compound's potential for clinical development. Based on this comprehensive evaluation (complied with Lipinski's and Veber's filters), sarsasapogenin was selected for in-depth biological characterization.

Figure 5. Sensitization of the MG63 and Saos-2 cell lines to a CSF1R inhibitor. (A) Dose-dependent cytotoxic effects of linarin on MG63 and Saos-2 osteosarcoma cell viability. (B) Dose-dependent cytotoxic effects of saikosaponin C on MG63 and Saos-2 osteosarcoma cell viability. (C) Dose-dependent cytotoxic effects of sarsasapogenin on MG63 and Saos-2 osteosarcoma cell viability. (D) Molecular docking of the binding between sarsasapogenin and the CSF1R protein. (E) Visual molecular docking of the binding acid residues between sarsasapogenin and CSF1R protein. (F) The expression levels of CSF1R as detected by western blot analysis in MG63 and Saos-2 cells after treatment with sarsasapogenin. *P<0.05. CSF1R, colony-stimulating factor-1 receptor.

To elucidate the molecular basis of the interaction between sarsasapogenin and CSF1R, molecular docking simulations were performed. Visual molecular docking analysis results about the interaction mode of sarsasapogenin with CSF1R protein 3D structure is observed in Fig. 5D, which showed a strong interaction of −10.6 Kcal/mol binding energy. These analyses revealed that sarsasapogenin binds to CSF1R through a combination of hydrogen bonds, van der Waals forces and hydrophobic interactions (Fig. 5E). Specifically, sarsasapogenin forms a key hydrogen bond with the ALA629 residue and engages in hydrophobic interactions with six additional residues: MET637, LEU640, LEU769, VAL647, ILE636 and HIS776.

Consistent with these in silico findings, western blot analysis demonstrated that CSF1R protein expression levels were significantly lower in sarsasapogenin-treated MG63 and Saos-2 cells compared with those in the untreated control cells (Fig. 5F). Collectively, these results demonstrated that sarsasapogenin inhibited osteosarcoma cell viability and downregulated CSF1R expression, supporting its potential as a CSF1R-targeted therapeutic agent.

MD simulation and analyses of CSF1R-sarsasopogenin complexes

RMSD analysis of the energy-minimized CSF1R protein revealed a flexible trajectory during the initial phase of the MD simulation, reflecting transient structural fluctuations. However, RMSD values stabilized after ~40 nsec, indicating that the CSF1R protein reached a structurally stable state for the remainder of the simulation. The RMSD of the CSF1R-sarsasapogenin complex was 2.3 Å (Fig. 6A), which confirmed the stability of the protein-ligand complex. The Rg of CSF1R bound to sarsasapogenin remained consistently <20 Å throughout the simulation (Fig. 6B), suggesting no significant changes in protein size or global shape and confirmed the preservation of its overall conformation. To evaluate local residue dynamics, the RMSF for the CSF1R-sarsasapogenin complex was calculated (Fig. 6C). The RMSF profile showed minimal overall residue fluctuation, with no pronounced peaks in the binding pocket region, further supporting the stability of the protein-ligand interaction. The SASA results revealed that the activity of the CSF1R-sarsasopogenin complex fluctuated slightly (Fig. 6D). The binding of small molecules is reported to affect the binding microenvironment and can lead to changes in SASA to a certain extent (33); two hydrogen bonds persisted throughout the entire 100-nsec simulation, which indicated that the interaction of sarsasopogenin with the protein, CSF1R, was both marked and lasting (Fig. 6E). The computed binding free energy for the CSF1R-sarsasopogenin complex was −34.38 Kcal/mol (Fig. 7A), with key contributions from residues ASP796, CYS774, HIS776, ILE636, ILE646, ILE775, LEU769 and MET637 (Fig. 7B).

Figure 6. Molecular dynamics simulation analyses. (A) RMSD curves of sarsasapogenin with the CSF1R protein. (B) The Rg analysis of sarsasapogenin with the CSF1R protein. (C) RMSF curves of sarsasapogenin with the CSF1R protein. (D) The SASA of sarsasapogenin with CSF1R protein. (E) H-bond number of sarsasapogenin with the CSF1R protein. CSF1R, colony-stimulating factor-1 receptor; RMSD, root mean square deviation; Rg, radius of gyration; RMSF, root mean square fluctuation; SASA, solvent-accessible surface area.

Figure 7. Relative binding free energies of sarsasapogenin-CSF1R complexes. (A) Binding free energy component analysis of sarsasapogenin with the CSF1R protein. (B) Key amino acid energy analysis of sarsasapogenin with the CSF1R protein. CSF1R, colony-stimulating factor-1 receptor. VDWAALs, Van der Waals; MM/PBSA, molecular mechanics Poisson-Boltzmann surface area.

Discussion

To systematically characterize ODRGs, the expression patterns of CSF1R, a key ODRG, was first analyzed in osteosarcoma and various types of cancer using multiomics datasets; specifically, the GEO database for osteosarcoma-focused analysis and the TCGA/GTEx database for pan-cancer analysis. CSF1R was significantly upregulated in osteosarcoma tissues compared with that in normal bone tissues. Pan-cancer analysis further revealed that CSF1R expression was aberrantly upregulated in multiple malignancies, with the highest fold-changes observed in DLBC, GBM, KIRC, KIRP, LAML, LGG, PAAD and TGCT, while moderate downregulation was noted in ACC. This widespread dysregulation across tumor types highlights CSF1R as a conserved oncogenic driver.

Further scRNA-seq analysis of the GSE162454 dataset confirmed that CSF1R was predominantly expressed in osteoblasts and monocytes/macrophages within the osteosarcoma microenvironment, which is consistent with its known expression in myeloid lineages, and its absence in hematopoietic stem cells (12,16). As a critical regulator of embryogenesis, tissue repair and immune cell homeostasis, dysregulated CSF1R signaling has been previously implicated in the initiation and progression of multiple tumor types (gastrointestinal stromal tumor, hepatocellular carcinoma, ovarian cancer, tenosynovial giant-cell tumor and triple-negative breast cancer) (16).

The present pharmacological experiments found that compared with the HUVEC control cell line, MG63 and Saos-2 osteosarcoma cell lines (which exhibit low endogenous CSF1R expression) were more sensitive to CSF1R inhibitors. Among the tested inhibitors, pexidartinib displayed the highest potency, as evidenced by its lower IC50 values compared with those of sotuletinib. Collectively, these findings support the potential of CSF1R as a diagnostic biomarker and therapeutic target for osteosarcoma.

In recent years, considerable progress has been made in the development of CSF1R inhibitors, with clinical efficacy already demonstrated in the management of tenosynovial giant cell tumors (34,35). However, safety concerns persist, primarily regarding off-target effects, most notably hepatotoxicity. For instance, current clinical guidelines mandate rigorous liver function monitoring when the CSF1R-selective inhibitor PLX3397 is administered due to documented risks of fetal liver injury in a preclinical study (35). Given the well-established therapeutic value of Chinese herbal medicines as a source of bioactive compounds for drug discovery, an in silico drug repurposing experiement was conducted to identify novel CSF1R inhibitors for osteosarcoma treatment. Through computational virtual screening of 13,227 herbal-derived compounds against the CSF1R kinase domain, three promising candidates were identified: Linarin, saikosaponin C and sarsasapogenin. Compared with the reference inhibitor sotuletinib, all three compounds exhibited stronger binding affinities to CSF1R.

Sarsasapogenin, a steroidal sapogenin primarily isolated from the rhizome of Anemarrhena asphodeloides Bunge (36), has emerged as a promising phytosteroid due to its broad spectrum of pharmacological activities, including anti-inflammatory, anticancer, antidiabetic, anti-osteoclastogenic and neuroprotective effects (37–43). Despite these well-characterized properties, its potential therapeutic effects and underlying mechanisms in osteosarcoma remain unexplored. In the present study, sarsasapogenin treatment dose-dependently inhibited the viability of both the MG63 and Saos-2 osteosarcoma cell lines, as assessed by CCK-8 assay. Western blot analysis further revealed that CSF1R protein expression was significantly lower in sarsasapogenin-treated cells than in vehicle-treated control cells, which suggested a potential mechanism underlying the antitumor activity of sarsasapogenin in osteosarcoma.

These findings demonstrated the efficacy of CSF1R inhibitors in reducing the viability of MG63 and Saos-2 cells; however, the present study had several limitations. First, the MAPK/ERK signaling pathways that mediate the therapeutic role of CSF1R in osteosarcoma require further investigation. Second, although molecular docking predicted strong sarsasapogenin-CSF1R binding, biochemical validation (such as cellular thermal shift assay or surface plasmon resonance assay) was not performed due to resource constraints. Future studies should employ protein- and cell-based binding assays to directly confirm this interaction.

In conclusion, the present findings systematically validated CSF1R as a clinically relevant diagnostic biomarker and actionable therapeutic target for osteosarcoma. Among the identified candidates, sarsasapogenin was identified as particularly promising lead compound, as it exerted potent growth-inhibitory effects on osteosarcoma cells through the specific downregulation of CSF1R. These results strongly supported its potential for further development as a natural-product-derived chemotherapeutic agent for osteosarcoma.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was funded by the Natural Science Foundation of Hunan Province (grant no. 2024JJ9532), the Chinese Medicine Research Project of Hunan Province (grant no. B2023048), the Natural Science Foundation of Changsha City (grant no. kq2403178) and the Hunan Medical Association Medical Research Program (grant nos. HMA202502017 and HMA202503013).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Author's contributions

HZ conceived and designed the study, performed experiments and analyzed the data. HZ and SW assisted with data analysis and interpretation. DL conducted the cell viability and western blot assays. JG supervised the project, participated in data analysis, and wrote the manuscript with contributions from all authors. HZ and DL confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References