Skip to main content
Download PDF

Abnormal signal pathways and tumor heterogeneity in osteosarcoma

Journal of Translational Medicine volume 21, Article number: 99 (2023) Cite this article

Abstract

Background

Osteosarcoma (OS) is the most frequent and aggressive primary malignant sarcoma among adolescents and chemotherapy has not substantially progressed for decades. New insights into OS development and therapeutic strategies are urgently needed.

Methods

We analyzed integrated single-cell transcriptomes, bulk RNA-seq, and microarray data from Gene Expression Omnibus (GEO) datasets. We also used Weighted Gene Co-expression Network Analysis (WGCNA), Gene set enrichment analysis (GSEA), and Gene set variation analysis (GSVA), along with Simple ClinVar and Enrichr web servers.

Results

The findings of integrated single-cell analysis showed that OS arises from imperfect osteogenesis during development. Novel abnormalities comprised deficient TGFβ and P53 signal pathways, and cell cycle pathway activation, and a potentially new driver mutation in the interferon induced transmembrane protein 5 (IFITM5) that might function as a pathogenic factor in OS. Osteosarcoma is characterized by oncocyte heterogeneity, especially in immunogenic and adipocyte-like subtypes that respectively promote and hamper OS treatment. Etoposide is a promising chemotherapeutic that provides palliation by affecting the subtype of OS and correcting the abnormal pathways.

Conclusion

Various abnormal signal pathways play indispensable roles in OS development. We explored the heterogeneity and underlying mechanisms of OS and generated findings that will assist with OS assessment and selecting optimal therapies.

Introduction

Osteosarcoma (OS) is an extremely aggressive sarcoma for which chemotherapeutic strategies have not substantially progressed for decades [1, 2]. The mechanism of OS awaits explanation and new molecular targets are urgently needed to improve the prognosis of patients. The causes of OS remain unclear in terms of oncogenetic mechanisms, but it might be due to genetic abnormalities, changes in tumor suppressor pathways (such as p53 and Rb), telomerase, and alternative telomere lengthening (ALT) inside bone cells [3]. Osteosarcoma driver mutations have been found in TP53, RB, PTCH1, MYC, NOTCH1, BRCA2, APC, and PRKAR1A genes. Clinical evaluations of targeted agents have yielded disappointing results [4, 5]. However, some new driver genes are currently under development as potential therapies for OS.

Interrelated biological signaling pathways in OS comprise WNT/βcatenin, Hedgehog, mTOR, and RANKL/NF-κB [6]. Achievements in the underlying molecular biology of OS have contributed to transformative advances in understanding this type of malignancy. Osteosarcoma originates from bone mesenchymal stem cells (MSCs) with osteoblastic lineage commitment. However, understanding the molecular mechanisms of OS development remains insufficient [3, 7]. For example, the TGFβ signal pathway exerts dual effects in terms of tumor prevention and carcinogenesis, depending on the timing of tumor development [8,9,10].

Osteosarcoma is significantly heterogeneous at the genomic, transcriptomic, and epigenetic levels. Malignant OS cells have stemness properties that are closely involved in chemotherapy resistance, relapse, or refractory and metastatic processes [11,12,13]. Current approaches to intra-tumor heterogeneity mainly depend on cellular genomic or transcriptional sequencing. However, this lacks high resolution in authenticating the complex cellular subtypes and intra-tumoral heterogeneity in OS. Single-cell RNA-sequencing (scRNA-seq) is a promising tool that can analyze the composition of heterogeneous cell populations and lineage hierarchies [14, 15].

Here we mapped single-cell (sc) RNA sequences of human embryonic long bones, normal bone tissue, and OS at single-cell high resolution to address the development of OS and oncocyte heterogeneity. The results of our integrated bioinformatic analysis identified and validated the process of osteogenesis imperfecta, potential new driver genes, abnormal signaling pathways, transcriptomic heterogeneity, and targeted therapeutic agents for OS.

Methods

Single cell RNA sequencing (scRNA-seq) data processing

We analyzed scRNA-seq data using R statistical software (v. 4.1.2, The R Foundation for Statistical Computing, Vienna, Austria), the Seurat v4.0 toolkit and MPLAB Harmony v1.0. Public scRNA-seq data were acquired from the GSE143753, GSE196678, and GSE162454 datasets in the Gene Expression Omnibus (GEO). Low-quality cells (< 3; features, < 200; > 10% mitochondrial genes) were filtered, the data were initially normalized, then batch effects were removed using LogNormalize (features = 3000) and the Harmony function (max.iter.harmony = 20). Dimensions were reduced using Uniform Manifold Approximation and Projection (UMAP). We found marker genes for cell groups in datasets using the COSine similarity-based marker gene identification (COSG) package [16]. We quantified marker genes using Z-scores and the results are shown in a heatmap. We calculated copy number variations to identify malignant oncocytes using the inferCNV package.

Cell types were identified in each subpopulation based on their known lineage markers: mesenchymal cells, vimentin (VIM); immune cells, cluster of differentiation (CD)45; myeloid cells: CD83, CD14, CD68; lymphoid cells, CD3 and B-cell maturation antigen (BCMA); endothelial cells, von Willebrand factor, (VWF) and carcinoembryonic antigen-related cell adhesion molecule-1 (CEACAM1); myogenic cells, myogenic differentiation 1 (MYOD1) and myogenin (MYOG); bone-related cells, collagen type 1 alpha 1 (COL1A1) and alpha 2 (COL1A2) chains, lumican (LUM), SRY-Box Transcription Factor 9 (SOX9), alkaline phosphatase (ALP) and RUNX Family Transcription Factor 2 (RUNX2)).

Pseudotime and trajectory analysis

Single-cell pseudotime and trajectory analyses were constructed using the Monocle3 toolkit (v. 2.14.0). Evolutionary processes were organized into potentially discontinuous trajectories by the learn_graph function. Pseudotime was defined using the order_cells function with a selected node representing development. Genes that were differentially expressed over these trajectories were then identified [17, 18].

Weighted gene co-expression network analysis (WGCNA)

We also searched related public expression profiles in 12 MSCs, 3 osteoblasts, and 84 OS samples (GSE33383). A co-expression network was constructed to identify co-expressed modules using the WGCNA package in R. The expression matrix was restricted to only the top 25% of expressed genes according to variance analyses. Relationships between gene sets (modules) were explored using a hierarchical cluster dendrogram. A clustering tree based on the eigengenes of modules calculated the dissimilarity of the module eigengenes (MEs). Associations between gene sets (modules) and clinical features were assessed using Pearson correlations.

Clinical datasets and resources

We downloaded 53 samples including expression profiles and clinical outcomes of patients with OS from GSE19743 for further analysis. Except for the expression profile, the clinical information mainly included living status, survival duration, Huvos grading scores, and metastatic status. Patients were assigned to high- (n = 27) and low- (n = 26) risk groups based on the median numbers of samples. We applied Cox regression and Kaplan–Meier curves to estimate overall survival.

Microarray data processing

We used the gene set with the expression profile GSE84500 to obtain significantly differentially expressed genes (DEGs). A data matrix was downloaded to annotate the probe into gene symbol sets. Significance was analyzed using Limma, or the Deseq2 package. The most significant changes in gene expression were identified using false discovery rates (FDR < 0.01) and absolute fold change (FC > 1).

Gene set enrichment analysis (GSEA)

We assessed gene set enrichment using GSEA [19]. Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway gene sets were curated from the Molecular Signatures Database (Human MSigDB database v2022.1.Hs). Immunogenic and adipocyte-like subtype molecular signatures were created by identifying significant marker genes in scRNA-seq. Normalized enrichment scores (NES) > 1 and p < 0.05 reliably filtered significant pathways.

Gene set variation analysis (GSVA)

We established the molecular signatures of OS oncocyte subtypes depending on the significant marker genes, which were acquired by single cell analysis [20]. We also calculated final scores for further analyses using GSVA (v. 1.47.0).

Databases and datasets

Simple ClinVar (https://simple-clinvar.broadinstitute.org/) is a web-server application that summarize variant, gene, gene- and disease-wise statistics based on the entire ClinVar database in a dynamic and user-friendly web-interface [21]. Enrichr (https://maayanlab.cloud/Enrichr/) is a robust web-server that contains many types datasets, among which we used the MGI Mammalian Phenotype, Jensen DISEASES, and BioPlanet pathway databases, as well as the PanglaoDB Augmented, and DrugMatrix datasets [22, 23]. Statistical significance was set at < 0.05.

Results

Osteosarcoma arises from osteogenesis imperfecta during development

We integrated data from 8-week-old long bones (GSE143753), human normal bone tissue (GSE196678) and OS (GSE162454) to determine whether a developmental disorder is associated with OS. After quality control and batch correction, we obtained 60,553 cells for downstream analysis. After dimensionality reduction, UMAP-based cell clustering identified myeloid, lymphoid, endothelial, myogenic, and bone-related cell clusters (Fig. 1A). Based on the lineage markers (Fig. 1B), then split them depending on their tissue origins. We found that bone mesenchymal stem cells (BMCs), normal osteocytes, and OS oncocytes mapped together, indicating that oncocytes might be derived from BMSCs due to a developmental disorder (Fig. 1C). We compared transcriptional profiles to identify DEGs between normal osteocytes and OS oncocytes during development. We identified 349 DEGs (log2FC > 1; P < 0.01), of which 143 were elevated in OS oncocytes (Fig. 1D; Additional file 2: Table S1). We analyzed the upregulated DEGs using the MGI Mammalian Phenotype and Jensen DISEASES databases to gain insight into their functional relevance. We found enrichment of decreased compact bone thickness, decreased trabecular bone volume, decreased length of long bones, and osteogenesis imperfecta (Fig. 1E and F; black arrows). This further supported the notion of a developmental disorder in OS. We analyzed trajectories to gain insight into the origin and development of bone-related cells (Fig. 2A and B). The results revealed partial BMSCs at the start of the developmental trajectory. This indicated that BMSCs can develop into OS.

Fig. 1
figure 1

Osteosarcoma arises from osteogenesis imperfecta during development. A Single cell RNA-sequence data (scRNA-seq) analysis for 8-week-old long bones (GSE143753), human normal bone tissue (GSE196678) and OS (GSE162454). Cluster analyzed by using Uniform Manifold Approximation and Projection (UMAP) method and populations were identified by color. B Markers used to define subpopulations. C Separated UMAP depending on tissue origins. D Volcano plot of significant different genes (| log2FC | > 1; P < 0.01) between normal osteocytes and OS oncocytes. E, F Enrichment of MGI Mammalian Phenotype and Jensen DISEASES by using the upregulated genes in D (p.value<0.01)

Full size image
Fig. 2
figure 2

IFITM5 is a potential new driver gene in OS development. A, B Single-cell trajectory analysis using Monocle3, cells colored by tissue origins and pseudotime. C Venn plot show the shared genes among the associated genes in osteogenesis imperfecta, and decreased bone thickness, volume, and length. D Plots of expression for ITIFM and MT2A by pseudotime. E, F UMAP of cell data set branch subsets. Cells are colored by the expression of IFITM5 and MT2A. G Various mutation categories of IFITM5 in Simple ClinVar datasets. H Coding variant mapping and domain annotation for IFITM5

Full size image

Interferon induced transmembrane protein 5 (IFITM5) is a potential new driver gene in OS development

Among the functional ontologies, only IFITM5 played a crucial role in the process of OS oncogenesis, including osteogenesis imperfecta, and decreased bone thickness, volume, and length (Fig. 2C). Furthermore, IFITM5 and Metallothionein 2A (MT2A) were the most significantly upregulated and downregulated genes, respectively, among the DEGs between normal osteocytes and OS oncocytes (Fig. 1D). The trajectory analyses also offered a clearer view of gene dynamics along a single path. We found that IFITM5 significantly increased from BMSCs to OS, whereas MT2A did not (Fig. 2D). The pseudotime-dependent IFITM5 gene changed as cells progressed along the oncocyte trajectory. Therefore knowing the order in which IFITM5 changes can lead to the generation of new models of OS development (Fig. 2E and F). We confirmed the IFITM5 variant using Simple ClinVar. We summarized various mutation categories of IFITM5 and found that a missense mutation was the most prevalent type (Fig. 2G), and that three variations of ITITM5, including the point mutation (c.-14C > T), the missense c.119C > T (p.S40L), and c.143A > G (p.N48S), can lead to pathogenic states according to the Simple ClinVar dataset and current literature (Fig. 2H). These findings together indicated that IFITM is a potential new driver gene in OS development that warrants further investigation.

Osteosarcoma was attributed to complicated regulation of abnormal signal pathways during development

We further investigated OS development by analyzing the expression profiles of 12 MSCs, 3 osteoblasts, and 84 OS samples (GSE33383). We established correlation networks and identified co-expression modules using WGCNA. Nine co-expression modules with 4,637 genes were obtained using the Dynamic Tree Cut algorithm (Fig. 3A). The clustering tree mainly showed two branches based on the module eigengenes (Fig. 3B). Correlations between modules and samples of MSCs, osteoblasts, and OS are shown in (Fig. 3C). The heat map shows that MEgreen and MEblue modules correlated negatively with MSCs and osteoblasts, and positively with OS. The MEgery module correlated positively with MSCs and osteoblasts, but negatively with OS. We analyzed BioPlanet pathways using the co-expressed genes in the MEgreen, MEblue and MEgery modules to gain insight into the signal pathway relevance (Additional file 3: Table S2). The three most significant pathways in each module, namely, Cell cycle, Antigen processing and presentation, and TGFβ regulation of extracellular matrix were enriched (Fig. 3D‒F). The module correlation coefficients for MSCs were consistent with those of osteoblasts, in contrast to OS (Fig. 3G and H). That is, OS was attributed to constitutive activation of the Cell cycle and Antigen processing and presentation pathways, and deficient TGFβ regulation of extracellular matrix signals. The atypical signal pathways also had complex connections to regulate bone development (ConsensusPathDB; Fig. 3I).

Fig. 3
figure 3

Osteosarcoma was attributed to complicated regulation of abnormal signal pathways during development. A Hierarchical cluster dendrogram of identified co-expressed genes in modules in 12 MSCs, 3 osteoblasts, and 84 OS samples (GSE33383). B Clustering tree based on the module eigengenes of modules. C Heatmap of the correlation between hub genes and cell types. D, E, F Enrichment of BioPlanet signal pathway in MEgreen, MEblue, and MEgrey modules. G, H The correlation analysis between MSCs and osteoblast, as well as MSCs and osteoblast by Pearson method. I The connections of atypical signal pathways in ConsensusPathDB datasets

Full size image

Transcriptomic heterogeneity of oncocytes in OS

We analyzed the transcriptomic heterogeneity of OS oncocytes by re-clustering bone-related cells in OS (n = 8758 cells) and found six distinct transcriptomic sub-clusters (Fig. 4A). The positive expression of COL1A1, COL1A2, ALP and the absence of CD45 revealed that the six sub-clusters were bone-related oncocytes (Additional file 1: Fig. S1 B). The copy number variation (CNV) was estimated using inferCNV (with myogenic cells as a reference) to distinguish oncocytes and normal cells (Fig. 4B). The CNV from scRNA-seq revealed numerous alterations in all sub-clusters compared with the reference cells (Fig. 4C). We determined the top 100 differentially expressed transcripts among the six sub-populations and defined their cell types in the Enrichr datasets (Additional file 4: Table S3). A heat map shows the hub genes in each sub-cluster (Fig. 4D). Six sub-populations were annotated and visualized by UMAP as neuron-like, immunogenic, fibroblastic, chondroblastic, adipocyte-like, and osteoblastic oncocytes (Fig. 4A). The top 100 gene signatures of each sub-population were assessed using GSVA and the oncocyte subtype abundance was estimated from 53 human transcriptomic profiles (GSE21257). Kaplan–Meier estimates associated higher immunogenic and adipocyte-like scores with longer progression-free survival (p = 0.0202), whereas and unfavorable prognoses, respectively (p = 0.0065). Other subtypes had no significant prognostic value (Fig. 4E). IFITM5 is a potential new driver gene in OS development in our study. Correlation analysis was performed to further identify the relationship between tumor subtypes and IFITM5. We found that the expression of IFITM5 has a strongly significant positive correlation with adipocyte-like and chondroblastic subtypes and a significant negative correlation with immunogenic and fibroblastic subtypes. Therefore, the high expression of IFITM5 was associated with the shorter overall survival subtype which is in line with the previous study (Additional file 1: Fig. S1A) [24].

Fig. 4
figure 4

Transcriptomic heterogeneity of oncocytes in OS. A Re-clustering of OS oncocytes. B The copy number variation (CNV) was estimated using inferCNV (with myogenic cells as a reference) to distinguish oncocytes and normal cells. The red and blue colors represent the copy number variation. C The CNV score of each cell was calculated as a quadraic sum of CNV. D A heat map shows the hub genes in each sub-cluster. E Kaplan–Meier survival curve shows the survival of OS patients with high- and low risk-scores for each subcluster signature (Nhigh-risk = 27 vs. Nlow-risk = 26)

Full size image

Deficient TGFβ signaling leads to adipocyte-like subtype of OS during development

The adipocyte-like OS subtype was associated with a poor prognosis. We profiled gene expression in human MSCs (hMSCs) in the presence of adipogenic/osteogenic factors (GSE84500) to assess the development of adipocyte-like subtype of OS. A comparison of transcriptional profiles between panels of hMSC_BMP2 + IBMX and hMSC_BMP2 + IBMX + TGFβ. Principal Component Analysis (PCA) initially revealed correlations (Fig. 5A). We identified 469 upregulated transcripts (log2FC > 1, P < 0.01), of which 240 were elevated in the hMSC_BMP2 + IBMX group (Additional file 5: Table S4) and named Adipocytes in PanglaoDB Augmented Datasets and 229 that were abundantly expressed in the hMSC_BMP2 + IBMX + TGFβ group (Additional file 6: Table S5) were labeled Chondrocytes, Osteoblasts, and Fibroblasts in the same datasets (Fig. 5B). Meanwhile, BioPlanet pathways confirmed activated TGFβ regulation of extracellular matrix in the hMSC_BMP2 + IBMX + TGFβ group. Therefore, we judged that BMP2 and IBMX cause hMSCs to differentiate into adipocytes, TGFβ inhibits this process and redirects these cells to differentiate into chondrocytes, osteoblasts, and fibroblasts. The TGFβ signal was defective during OS development (Fig. 3C and F), and adipocytes could not differentiate into normal bone cells without the TGFβ signal. This stalled adipocyte development and led to the adipocyte-like subtype of OS (Fig. 5D). In clinical translation, patients with metastasis are more likely to have high-risk scores (Fig. 5E). We applied the Huvos grading system to judge the effectiveness of neoadjuvant chemotherapy on OS; higher scores reflect larger necrotic areas in the clinical data (GSE21257). However, the risk score of the adipocyte-like subtype did not differ among grades, indicating that the adipocyte-like subtype is resistant to neoadjuvant chemotherapy (Fig. 5F).

Fig. 5
figure 5

Deficient TGFβ signaling leads to adipocyte-like subtype of OS during development. A Principal component analysis (PCA) biplot of transcriptional profiles with 95% confidence ellipsoids in the panels of hMSC, hMSC_BMP2+IBMX and hMSC_BMP2+IBMX+TGFβ (GSE84500). B Heat map shows the significantly different genes in these three groups (log2FC > 1; P < 0.01), and enrichment of cell types in PanglaoDB Augmented Datasets. C Enrichment of BioPlanet signal pathway in hMSC_BMP2+IBMX+TGFβ group. D Sketch map of TGFβ signaling in OS development. E Box plot and corresponding stacked histogram demonstrate the correlation between metastasis and adipocyte-like subtype score. F Huvos grading system was used to judge the effectiveness of neoadjuvant chemotherapy on OS

Full size image

Etoposide exerts effective palliation by affecting the subtype of OS and correcting abnormal pathways

Patients with OS and higher-, than lower-risk scores for the immunogenic subtype survived longer (Fig. 4E) and metastasis was less likely to develop (Fig. 6A). Furthermore, the Huvos system associated the immunogenic subtype with higher grades and sensitivity to neoadjuvant chemotherapy (Fig. 6B). We applied the immunogenic subtype gene signature to the DrugMatrix dataset that contains comprehensive information about gene expression in rats under treatment with various drugs. Transcriptional changes caused by thioguanine and etoposide were enriched (Fig. 6C). An extensive literature review revealed that etoposide is a feasible candidate as salvage therapy for relapsed and metastatic OS. A high-throughput transcriptomic screen of OS cells (U-2 OS) exposed to seven concentrations of etoposide was evaluated by GSEA (GSE200845) to determine the underlying mechanism of etoposide. The GSEA enrichment plot associated DEGs with the activation of immunogenic subtype signatures (p = 0.000, NES = 1,80; Fig. 6D left). The heat map also shows that the top 10 genes with immunogenic subtype signatures gradually increased as the dose increased (Fig. 6D right). Meanwhile, the enrichment plot and a heat map showed that the adipocyte-like subtype signatures correlated negatively with the dose-dependent profile of etoposide treatment (p = 0.015; NES = − 1.49; Fig. 6E). We also used GSVA to determine a dose-dependent tendency between immunogenic and adipocyte-like subtype signatures. The results were the same as those of the GSEA analysis (p < 0.01, Fig. 6F). We identified the high-risk genes checkpoint kinase 2 (CHEK2) and retinoblastoma transcriptional corepressor 1 (RB1), which represented 46.15% and 15.39% of the OS driver genes respectively, and play crucial roles in the cell cycle pathway, which we confirmed was activated (Figs. 6G; 3D). The mutation rate of the important tumor suppressor gene tumor protein 53 (TP53) was also 15.39% in OS (Fig. 6G). Etoposide can inhibit the cell cycle pathway, and activate the P53 signaling pathway (Fig. 6H and I). The immunogenic OS subtype originated in MSCs, because it correlated closely with Antigen processing and presentation as described in Fig. 3E (Fig. 6J). Therefore, etoposide has promise as a palliative therapeutic by promoting immunogenic subtypes, inhibiting adipocyte-like subtypes, and correcting the abnormal cell cycle and p53 signaling pathways.

Fig. 6
figure 6

Etoposide exerts effective palliation by affecting the subtype of OS and correcting abnormal pathways. A Box plot and corresponding stacked histogram show immunogenic subtype score in metastasis and no metastasis OS samples. B Immunogenic subtype score in OS samples across various Huvos grades. C Enrichment of immunogenic subtype gene signature in DrugMatrix dataset. D, E The GSEA enrichment plot and dose-dependent heatmap showed the effect of etoposide on immunogenic and adipocyte-like subtypes. F GSVA analysis was used to determine the different tendencies between immunogenic and adipocyte-like subtypes with different concentrations of etoposide. G Osteosarcoma associated mutant genes in Simple ClinVar datasets. H, I The GSEA enrichment plot and dose-dependent heatmap revealed the effect of etoposide on the cell cycle and p53 signaling pathway. J The correlation analysis between immunogenic subtype and Antigen processing and presentation pathway. K Overview of the abnormal signal pathways and tumor heterogeneity in osteosarcoma

Full size image

Discussion

Osteosarcoma is the most frequent primary malignant sarcoma in adolescents. Given the treatment dilemma, new insights into the OS development and therapeutic strategies are urgently needed [25, 26]. Here, we discovered using integrated single-cell analysis, that OS arises from osteogenesis imperfecta during development. We identified deficient TGFβ and P53 signal pathways, an activated cell cycle pathway, and a potentially novel driver gene IFITM5 mutant as potential pathogenic factors of OS. We characterized OS heterogeneity and found that the immunogenic and adipocyte-like subtypes respectively exerted beneficial and detrimental effects on OS. Etoposide is a promising chemotherapeutic drug that achieves palliation by affecting the OS subtypes and correcting aberrant pathways.

Osteosarcoma is derived from MSCs and is characterized by osteogenesis imperfecta [27, 28]. The IFITM5 gene encodes an osteoblast-specific membrane protein that is an established positive regulatory factor during bone mineralization [29, 30]. Single base variants (c.119 C > T), c.143A > G and (c.-14 C > T) in the coding region of IFITM5 have been identified in patients with osteogenesis imperfecta type V [31, 32]. Furthermore, IFITM5 is overexpressed in abnormal bone hyperplasia in rat primary osteoblasts, UMR106 rat osteosarcoma cells, human primary osteoblasts and Saos-2 human osteosarcoma cells [33,34,35]. There are already proven cases published about osteosarcoma occurring in osteogenesis imperfecta due to IFITM5 mutation [36, 37]. The skeletal disorders caused by alterations in the IFITM5 gene, including c.-14C > T, c.119C > T (p.S40L), c.143A > G (p.N48S), while more case series are required to establish detailed genotype-phenotypes for these alterations in the IFITM5 gene [32, 38]. Here, we determined that IFITM is a potential new driver gene in OS development depending on pseudotime using trajectory analysis and the ClinVar database.

The various abnormal signal pathways that prompt MSCs towards abnormal bone growth also play indispensable roles in OS development [6, 25]. The dysregulated cell cycle (CHEK2, RB1 mutant) and p53 (TP53 deficiency) signaling pathways often result in pathogenesis and aberrant OS growth. The p53 signal pathway acts as a tumor suppressor and exerts a crucial role in safeguarding our body from developing OS [39]. There are many studies reported that human osteosarcomas can have a deletion, mutations, and/or rearrangements of the p53 gene, which may cause loss of normal constraints on cellular growth, cell cycle, senescence and metabolism [40,41,42].P53 also regulates osteogenic, chondrogenic, myogenic, adipogenic differentiation of MSCs [43, 44]. The present findings are also in line with the fact that Rb and p53 knockout in bone-marrow-derived MSCs (BM-MSCs) results in OS-like tumors that can be serially transplanted [45, 46]. The TGF-β signal favors normal bone formation in the mesenchymal osteoblastic lineage [47, 48]. However, TGF-β seems to mainly have a pro-tumoral effect on OS [49,50,51,52]. Several clinical strategies to block TGF-β signaling pathways, such as neutralizing antibodies (GC-1008), ligands or receptors inhibitors (AP12009), or chemical compounds (SB-431542, or LY2157299), have not been successful against OS [10, 53,54,55,56]. In contrast to previous findings, we have the new perspective that a TGF-β deficiency is a major cause of OS and leads to a detrimental adipocyte-like subtype. This is an alternative experimental approach to treating patients with OS. Briefly, the TGF-β signal pathway in MSCs guards bone formation in the early stage and then compromises their osteogenic differentiation in the late stage, therefore it dictates the conditions for bone development [47, 57]. In this context, TGF-β signal pathway mutations and/or alterations mutations to TGF-β cascade components have been associated with several bone disorders and many carcinomas [58, 59].Several high-throughput genomic and transcriptomic studies have delineated the intra-tumoral heterogeneity of OS [11, 60,61,62]. Zhang et.al identified seven OS tumor cell clusters with three differentiation branches by using single-cell RNA sequencing of conventional OS and cancellous bone (CB) samples, which have different prognoses and possible drug sensitivities [63]. Jiang et.al classified OS into four subtypes according to the genomic, epigenomic, and transcriptomic data, while our classification depended on the single cell sequencing. We used different methods to identify the heterogeneity of OS, but partially with similar molecular features and clinical prognosis, such as Immune activated (S-IA) and our immunogenic subtype both have immunological traits and better clinical prognosis. Their Immune suppressed (S-IS) subtype has activated adipogenesis and fatty acid metabolism-related pathways and encodes the fatty acid scavenger receptor CD36 which is in line with our adipocyte-like subtype [64]. We explored the heterogeneity and the underlying mechanism of OS that could help to provide better, customized therapy. We identified six oncocyte subtypes of OS, of which four did not affect prognosis (Additional file 1: Fig. S1C and D). The immunogenic and adipocyte-like subtypes were respectively associated with a better and worse prognosis. To our knowledge, this is a novel molecular signature classification of OS, which will facilitate OS assessment, and choosing the optimal therapy. For example, etoposide is associated with the activation and inhibition, respectively, of the immunogenic and adipocyte-like subtype signatures, and combining it with other strategies has shown promising anti-tumour activity in clinical trials [65,66,67,68].

However, it should be noted that the study we examined is limited to bioinformatics analysis, and a significant number of challenges still need to be overcome for a successful in vitro and vivo study. Large-scale samples including both single-cell sequencing and clinical outcomes should be taken into account in the future. In addition, Some different signaling networks in osteosarcoma, including RANKL/RANK, Wnt, Notch, PI3K/Akt/mTOR, and mechanotransduction pathways, contribute to osteosarcoma progression and metastasis and tumor heterogeneity, which certainly looks worthy of further investigation.

Conclusions

We attributed the occurrence of OS to an IFITM5 mutant, deficient TGFβ signaling, continuous activation of the cell cycle signal pathway, and P53 signal inhibition during development. Among six OS subtypes, the immunogenic subtype was beneficial, whereas the adipocyte-like subtype was detrimental to health. Etoposide could be a useful palliative because it altered the OS subtype, inhibited the cell cycle pathway, and improved the p53 signal pathway.

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Abbreviations

ALP:

Alkaline phosphatase

ALT:

Alternative telomere lengthening

APC:

Antigen presenting cell

BCMA:

B-cell maturation antigen (endothelial cells)

BM-MSCs:

Bone marrow mesenchymal stem cells

BMP2:

Bone morphogenetic protein 2

BRCA2:

Breast cancer gene 2

CD:

Cluster of differentiation

CEACAM1:

Carcinoembryonic antigen-related cell adhesion molecule-1

CHEK2:

Checkpoint kinase 2

CNV:

Copy number variation

COL1A1:

Collagen type 1 alpha 1

COL1A2:

Collagen type 1 alpha 2

COSG:

COSine similarity-based marker Gene identification

GEO:

Gene Expression Omnibus

GSEA:

Gene set enrichment analysis

GSVA:

Gene set variation analysis

hMSC:

Human mesenchymal stem cells

IBMX:

3-Isobutyl-1-methylxanthine

IFITM5:

Interferon induced transmembrane protein 5

LUM:

Lumican

MT2A:

Metallothionein 2A

mTOR:

Mammalian target of rapamycin

MYC:

Myelocytomatosis proto-oncogene

MYOD1:

Myogenic differentiation 1

MYOG:

Myogenin

NF-κB:

Nuclear factor kappa B

NOTCH1:

Neurogenic locus notch homolog protein 1

OS:

Osteosarcoma

PRKAR1A:

Protein Kinase CAMP-Dependent Type I Regulatory Subunit Alpha

PTCH1:

Protein patched homolog 1

RANKL:

Receptor activator of nuclear factor kappa beta

RB1:

Retinoblastoma transcriptional corepressor 1

RUNX2:

RUNX Family Transcription Factor 2

scRNA-seq:

Single-cell RNA-sequencing

SOX9:

SRY-Box Transcription Factor 9

TGFβ:

Transforming growth factor beta

TP53:

Transcriptional repressor 53

UMAP:

Uniform Manifold Approximation and Projection

VIM:

Vimentin

VWF:

Von Willebrand factor

WGCNA:

Weighted Gene Co-expression Network Analysis

WNT:

Wingless-related integration site

References

  1. Xie L, Ji T, Guo W. Anti-angiogenesis target therapy for advanced osteosarcoma. Oncol Rep. 2017;38:625–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Yu D, Zhang S, Feng A, Xu D, Zhu Q, Mao Y, Zhao Y, Lv Y, Han C, Liu R. Methotrexate, doxorubicin, and cisplatinum regimen is still the preferred option for osteosarcoma chemotherapy: a meta-analysis and clinical observation. Medicine. 2019;98:e15582.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Gorlick R, Anderson P, Andrulis I, Arndt C, Beardsley GP, Bernstein M, Bridge J, Cheung N-K, Dome JS, Ebb D. Biology of childhood osteogenic sarcoma and potential targets for therapeutic development: meeting summary. Clin Cancer Res. 2003;9:5442–53.

    PubMed  Google Scholar 

  4. Walkley CR, Qudsi R, Sankaran VG, Perry JA, Gostissa M, Roth SI, Rodda SJ, Snay E, Dunning P, Fahey FH. Conditional mouse osteosarcoma, dependent on p53 loss and potentiated by loss of Rb, mimics the human disease. Genes Dev. 2008;22:1662–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Czarnecka AM, Synoradzki K, Firlej W, Bartnik E, Sobczuk P, Fiedorowicz M, Grieb P, Rutkowski P. Molecular biology of osteosarcoma. Cancers. 2020;12:2130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Saraf AJ, Fenger JM, Roberts RD. Osteosarcoma: accelerating progress makes for a hopeful future. Front Oncol. 2018;8:4.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Cortini M, Avnet S, Baldini N. Mesenchymal stroma: role in osteosarcoma progression. Cancer Lett. 2017;405:90–9.

    Article  CAS  PubMed  Google Scholar 

  8. Roberts AB, Wakefield LM. The two faces of transforming growth factor β in carcinogenesis. Proc Natl Acad Sci. 2003;100:8621–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Principe DR, Doll JA, Bauer J, Jung B, Munshi HG, Bartholin L, Pasche B, Lee C, Grippo PJ. TGF-β: duality of function between tumor prevention and carcinogenesis. JNCI J Natl Cancer Inst. 2014;106:djt369.

    Article  PubMed  Google Scholar 

  10. Colak S, Ten Dijke P. Targeting TGF-β signaling in cancer. Trends Cancer. 2017;3:56–71.

    Article  CAS  PubMed  Google Scholar 

  11. Schiavone K, Garnier D, Heymann M-F, Heymann D. The heterogeneity of osteosarcoma: the role played by cancer stem cells. Stem Cells Heterogeneity Cancer. 2019;1139:187–200.

    Article  CAS  Google Scholar 

  12. Mutsaers AJ, Walkley CR. Cells of origin in osteosarcoma: mesenchymal stem cells or osteoblast committed cells? Bone. 2014;62:56–63.

    Article  PubMed  Google Scholar 

  13. Mohseny AB, Szuhai K, Romeo S, Buddingh EP, Briaire-de Bruijn I, de Jong D, van Pel M, Cleton-Jansen AM, Hogendoorn PC. Osteosarcoma originates from mesenchymal stem cells in consequence of aneuploidization and genomic loss of Cdkn2. J Pathol. 2009;219:294–305.

    Article  CAS  PubMed  Google Scholar 

  14. Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, Rotem A, Rodman C, Lian C, Murphy G. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352:189–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Mereu E, Lafzi A, Moutinho C, Ziegenhain C, McCarthy DJ, Álvarez-Varela A, Batlle E, Grün D, Lau JK, Boutet SC. Benchmarking single-cell RNA-sequencing protocols for cell atlas projects. Nat Biotechnol. 2020;38:747–55.

    Article  CAS  PubMed  Google Scholar 

  16. Dai M, Pei X, Wang X-J. Accurate and fast cell marker gene identification with COSG. Briefings Bioinform. 2022;23:bbab579.

    Article  Google Scholar 

  17. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn JL. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32:381–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, Trapnell C. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14:979–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102:15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Ferreira MR, Santos GA, Biagi CA, Silva Junior WA, Zambuzzi WF. GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J Biomed Mater Res, Part A. 2021;109:1004–14.

    Article  CAS  Google Scholar 

  21. Pérez-Palma E, Gramm M, Nürnberg P, May P, Lal D. Simple ClinVar: an interactive web server to explore and retrieve gene and disease variants aggregated in ClinVar database. Nucleic Acids Res. 2019;47:W99–105.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44:W90-97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Xie Z, Bailey A, Kuleshov MV, Clarke DJB, Evangelista JE, Jenkins SL, Lachmann A, Wojciechowicz ML, Kropiwnicki E, Jagodnik KM, et al. Gene set knowledge discovery with Enrichr. Curr Protoc. 2021;1: e90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Jiang H, Du H, Liu Y, Tian X, Xia J, Yang S. Identification of novel prognostic biomarkers for osteosarcoma: a bioinformatics analysis of differentially expressed genes in the mesenchymal stem cells from single-cell sequencing data set. Transl Cancer Res. 2022;11:3841–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Rathore R, Van Tine BA. Pathogenesis and current treatment of osteosarcoma: perspectives for future therapies. J Clin Med. 2021;10:1182.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Raymond AK, Jaffe N. Osteosarcoma multidisciplinary approach to the management from the pathologist’s perspective. Pediatr Adolesc Osteosarcoma. 2009;152:63–84.

    Article  Google Scholar 

  27. Broadhead ML, Sivaji S, Balogh Z, Choong PF. Osteosarcoma: from molecular biology to mesenchymal stem cells. Osteosarcoma Biol Behav Mech 2017.

  28. Abarrategi A, Tornin J, Martinez-Cruzado L, Hamilton A, Martinez-Campos E, Rodrigo JP, González M, Baldini N, Garcia-Castro J, Rodriguez R. Osteosarcoma: cells-of-origin, cancer stem cells, and targeted therapies. Stem Cells Int. 2016;2016:3631764.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Moffatt P, Salois P, Gaumond MH, St-Amant N, Godin É, Lanctôt C. Engineered viruses to select genes encoding secreted and membrane-bound proteins in mammalian cells. Nucleic Acids Res. 2002;30:4285–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Hanagata N, Takemura T, Monkawa A, Ikoma T, Tanaka J. Phenotype and gene expression pattern of osteoblast-like cells cultured on polystyrene and hydroxyapatite with pre-adsorbed type-I collagen. J Biomed Mater Res Part A. 2007;83:362–71.

    Article  Google Scholar 

  31. Hanagata N, Li X, Morita H, Takemura T, Li J, Minowa T. Characterization of the osteoblast-specific transmembrane protein IFITM5 and analysis of IFITM5-deficient mice. J Bone Miner Metab. 2011;29:279–90.

    Article  CAS  PubMed  Google Scholar 

  32. Pachajoa H, Giraldo-Ocampo S. A patient with bone fragility, multiple fractures, osteosarcoma, and the variant c. 143A> G in the IFITM5 gene: a case report. Orthop Res Rev. 2022;14:453–8.

    PubMed  PubMed Central  Google Scholar 

  33. Liu BY, Lu YQ, Han F, Wang Y, Mo XK, Han JX. Effects of the overexpression of IFITM5 and IFITM5 c.-14C> T mutation on human osteosarcoma cells. Oncol Lett. 2017;13:111–8.

    Article  CAS  PubMed  Google Scholar 

  34. Moffatt P, Gaumond MH, Salois P, Sellin K, Bessette MC, Godin É, de Oliveira PT, Atkins GJ, Nanci A, Thomas G. Bril: a novel bone-specific modulator of mineralization. J Bone Miner Res. 2008;23:1497–508.

    Article  CAS  PubMed  Google Scholar 

  35. Gorski JP, Wang A, Lovitch D, Law D, Powell K, Midura RJ. Extracellular bone acidic glycoprotein-75 defines condensed mesenchyme regions to be mineralized and localizes with bone sialoprotein during intramembranous bone formation. J Biol Chem. 2004;279:25455–63.

    Article  CAS  PubMed  Google Scholar 

  36. Mathew SE, Santhanam M, Madhuri V. Interferon-induced transmembrane protein 5 mutation causing type-V osteogenesis imperfecta: a case report. JBJS Case Connector. 2015;5: e15.

    Article  PubMed  Google Scholar 

  37. Maiya S, Grimer R, Ramaswamy R, Deshmukh N. Osteosarcoma occurring in osteogenesis imperfecta tarda. Int Orthop. 2002;26:126–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zhytnik L, Maasalu K, Duy BH, Pashenko A, Khmyzov S, Reimann E, Prans E, Kõks S, Märtson A. IFITM5 pathogenic variant causes osteogenesis imperfecta V with various phenotype severity in Ukrainian and Vietnamese patients. Hum Genomics. 2019;13:1–11.

    Article  CAS  Google Scholar 

  39. Velletri T, Xie N, Wang Y, Huang Y, Yang Q, Chen X, Chen Q, Shou P, Gan Y, Cao G. P53 functional abnormality in mesenchymal stem cells promotes osteosarcoma development. Cell Death Dis. 2016;7:e2015–e2015.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Fu H-L, Shao L, Wang Q, Jia T, Li M, Yang D-P. A systematic review of p53 as a biomarker of survival in patients with osteosarcoma. Tumor Biology. 2013;34:3817–21.

    Article  CAS  PubMed  Google Scholar 

  41. Gokgoz N, Wunder JS, Mousses S, Eskandarian S, Bell RS, Andrulis IL. Comparison of p53 mutations in patients with localized osteosarcoma and metastatic osteosarcoma. Cancer. 2001;92:2181–9.

    Article  CAS  PubMed  Google Scholar 

  42. Miller CW, Aslo A, Tsay C, Slamon D, Ishizaki K, Toguchida J, Yamamuro T, Lampkin B, Koeffler HP. Frequency and structure of p53 rearrangements in human osteosarcoma. Can Res. 1990;50:7950–4.

    CAS  Google Scholar 

  43. Molchadsky A, Shats I, Goldfinger N, Pevsner-Fischer M, Olson M, Rinon A, Tzahor E, Lozano G, Zipori D, Sarig R. p53 plays a role in mesenchymal differentiation programs, in a cell fate dependent manner. PLoS ONE. 2008;3: e3707.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Molchadsky A, Rivlin N, Brosh R, Rotter V, Sarig R. p53 is balancing development, differentiation and de-differentiation to assure cancer prevention. Carcinogenesis. 2010;31:1501–8.

    Article  CAS  PubMed  Google Scholar 

  45. Cheng L, Wang C, Jing J. Cell cycle kinases in osteosarcoma: potential for therapeutic intervention. Curr Pharm Des. 2016;22:4830–4.

    Article  CAS  PubMed  Google Scholar 

  46. Rubio R, Gutierrez-Aranda I, Sáez-Castillo A, Labarga A, Rosu-Myles M, González-García S, Toribio ML, Menendez P, Rodriguez R. The differentiation stage of p53-Rb-deficient bone marrow mesenchymal stem cells imposes the phenotype of in vivo sarcoma development. Oncogene. 2013;32:4970–80.

    Article  CAS  PubMed  Google Scholar 

  47. Janssens K, Ten Dijke P, Janssens S, Van Hul W. Transforming growth factor-β1 to the bone. Endocr Rev. 2005;26:743–74.

    Article  CAS  PubMed  Google Scholar 

  48. Juárez P, Guise TA. TGF-β in cancer and bone: implications for treatment of bone metastases. Bone. 2011;48:23–9.

    Article  PubMed  Google Scholar 

  49. Lamora A, Talbot J, Mullard M, Brounais-Le Royer B, Redini F, Verrecchia F. TGF-β signaling in bone remodeling and osteosarcoma progression. J Clin Med. 2016;5:96.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Lamora A, Talbot J, Bougras G, Amiaud J, Leduc M, Chesneau J, Taurelle J, Stresing V, Le Deley MC, Heymann MF. Overexpression of smad7 blocks primary tumor growth and lung metastasis development in osteosarcoma. Clin Cancer Res. 2014;20:5097–112.

    Article  CAS  PubMed  Google Scholar 

  51. Yang R-S, Wu C-T, Lin K-H, Hong R-L, Liu T-K, Lin K-S. Relation between histological intensity of transforming growth factor-β isoforms in human osteosarcoma and the rate of lung metastasis. Tohoku J Exp Med. 1998;184:133–42.

    Article  CAS  PubMed  Google Scholar 

  52. Xu S, Yang S, Sun G, Huang W, Zhang Y. Transforming growth factor-beta polymorphisms and serum level in the development of osteosarcoma. DNA Cell Biol. 2014;33:802–6.

    Article  CAS  PubMed  Google Scholar 

  53. Mohammad KS, Javelaud D, Fournier PG, Niewolna M, McKenna CR, Peng XH, Duong V, Dunn LK, Mauviel A, Guise TA. TGF-β-RI kinase inhibitor SD-208 reduces the development and progression of melanoma bone metastasesTGF-β inhibition to treat melanoma bone metastases. Can Res. 2011;71:175–84.

    Article  CAS  Google Scholar 

  54. Hau P, Jachimczak P, Schlingensiepen R, Schulmeyer F, Jauch T, Steinbrecher A, Brawanski A, Proescholdt M, Schlaier J, Buchroithner J. Inhibition of TGF-β 2 with ap 12009 in recurrent malignant gliomas: from preclinical to phase I/II studies. Oligonucleotides. 2007;17:201–12.

    Article  CAS  PubMed  Google Scholar 

  55. Morris JC, Tan AR, Olencki TE, Shapiro GI, Dezube BJ, Reiss M, Hsu FJ, Berzofsky JA, Lawrence DP. Phase I study of GC1008 (fresolimumab): a human anti-transforming growth factor-beta (TGFβ) monoclonal antibody in patients with advanced malignant melanoma or renal cell carcinoma. PLoS ONE. 2014;9: e90353.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Rodon Ahnert J, Baselga J, Calvo E, Seoane J, Brana I, Sicart E, Gueorguieva I, Cleverly A, Lahn M, Pillay S. First human dose (FHD) study of the oral transforming growth factor-beta receptor I kinase inhibitor LY2157299 in patients with treatment-refractory malignant glioma. J Clin Oncol. 2011;29:3011–3011.

    Article  Google Scholar 

  57. Maeda S, Hayashi M, Komiya S, Imamura T, Miyazono K. Endogenous TGF-β signaling suppresses maturation of osteoblastic mesenchymal cells. EMBO J. 2004;23:552–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Alliston T, Piek E, Derynck R. TGF-beta family signaling in skeletal development, maintenance, and disease. Cold Spring Harbor Monogr Ser. 2008;50:667.

    Google Scholar 

  59. Levy L, Hill CS. Alterations in components of the TGF-β superfamily signaling pathways in human cancer. Cytokine Growth Factor Rev. 2006;17:41–58.

    Article  CAS  PubMed  Google Scholar 

  60. Martin JW, Squire JA, Zielenska M. The genetics of osteosarcoma. Sarcoma. 2012;2012:627254.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Chen X, Bahrami A, Pappo A, Easton J, Dalton J, Hedlund E, Ellison D, Shurtleff S, Wu G, Wei L. Recurrent somatic structural variations contribute to tumorigenesis in pediatric osteosarcoma. Cell Rep. 2014;7:104–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Zhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, Wang Y, Zhang Z, Yuan T, Ding X. Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma. Nat Commun. 2020;11:1–17.

    Article  CAS  Google Scholar 

  63. Adamopoulos C, Gargalionis AN, Basdra EK, Papavassiliou AG. Deciphering signaling networks in osteosarcoma pathobiology. Exp Biol Med. 2016;241:1296–305.

    Article  CAS  Google Scholar 

  64. Jiang Y, Wang J, Sun M, Zuo D, Wang H, Shen J, Jiang W, Mu H, Ma X, Yin F. Multi-omics analysis identifies osteosarcoma subtypes with distinct prognosis indicating stratified treatment. Nat Commun. 2022;13:1–17.

    Article  Google Scholar 

  65. Gaspar N, Venkatramani R, Hecker-Nolting S, Melcon SG, Locatelli F, Bautista F, Longhi A, Lervat C, Entz-Werle N, Casanova M. Lenvatinib with etoposide plus ifosfamide in patients with refractory or relapsed osteosarcoma (ITCC-050): A multicentre, open-label, multicohort, phase 1/2 study. Lancet Oncol. 2021;22:1312–21.

    Article  CAS  PubMed  Google Scholar 

  66. Akazawa R, Umeda K, Saida S, Kato I, Hiramatsu H, Sakamoto A, Arakawa Y, Sumiyoshi S, Okamoto T, Moritake H. Temozolomide and etoposide combination for the treatment of relapsed osteosarcoma. Jpn J Clin Oncol. 2020;50:948–52.

    Article  PubMed  Google Scholar 

  67. Gaspar N, Occean B-V, Pacquement H, Bompas E, Bouvier C, Brisse HJ, Castex M-P, Cheurfa N, Corradini N, Delaye J. Results of methotrexate-etoposide-ifosfamide based regimen (M-EI) in osteosarcoma patients included in the French OS2006/sarcome-09 study. Eur J Cancer. 2018;88:57–66.

    Article  CAS  PubMed  Google Scholar 

  68. Perret A, Dômont J, Chamseddine AN, Dumont SN, Verret B, Briand S, Court C, Lazure T, Adam J, Ngo C. Efficacy and safety of oral metronomic etoposide in adult patients with metastatic osteosarcoma. Cancer Med. 2021;10:230–6.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We are very grateful to the participating universities and China Scholarship Council.

Funding

This work was supported by the National Natural Science Foundation of China (81702667), Natural Science Foundation of Shandong (ZR202111050005), China Postdoctoral Science Foundation (8206300728). This work also funded by China Scholarship Council.

Author information

Authors and Affiliations

  1. Department of Orthopedic Surgery, The First Affiliated Hospital of Shandong First Medical University & Shandong Provincial Qianfoshan Hospital, Shandong Key Laboratory of Rheumatic Disease and Translational Medicine, Jinan, 250014, Shandong, People’s Republic of China

    Yifeng Sun, Chunming Zhang, Wenqiang Zhang & Wei Liu

  2. Department of Surgery, Ulm University Hospital, Ulm University, Ulm, Germany

    Yifeng Sun

  3. MOE Key Laboratory of Cell Proliferation and Differentiation, School of Life Sciences, Peking University, Beijing, 100871, China

    Qiongxuan Fang

Authors
  1. Yifeng Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chunming Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Qiongxuan Fang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Wenqiang Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Wei Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

WL and CZ conceived the study and collected the data, YS analysed the data, generated figures and wrote the manuscript. QF performed a part of data analysis. WZ help to generated figures. All authors reviewed the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Wei Liu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig.S1.

A The correlation analysis between IFITM5 expression and each tumor subtype. B Markers used to define subpopulations. C Box plots show the correlation between metastasis and other subtype scores. D Huvos grades of neuron-like, fibroblastic, chondroblastic, and osteoblastic subtypes.

Additional file 2: Table S1.

Different genes between normal osteocytes and OS oncocytes.

Additional file 3: Table S2.

Signatures of each module in WGCNA analysis.

Additional file 4: Table S3.

Signatures of each OS subtype.

Additional file 5: Table S4.

Significantly different genes of the hMSC_BMP2+IBMX group.

Additional file 6: Table S5.

Significantly different genes of the hMSC_BMP2+IBMX+TGFβ group.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, Y., Zhang, C., Fang, Q. et al. Abnormal signal pathways and tumor heterogeneity in osteosarcoma. J Transl Med 21, 99 (2023). https://doi.org/10.1186/s12967-023-03961-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12967-023-03961-7

Keywords

Journal of Translational Medicine

ISSN: 1479-5876

Contact us