Single-cell sequencing data and clinical information for 11 OS patients were downloaded from GSE152048. The 10X genomics data for each individual sample were loaded by the Seurat package (v4.0.0) in R software (v4.0.2) (21). Low-quality cells with ≤ 300 detected genes or ≥ 10% mitochondrial genes were removed. Feature counts for each cell were divided by the total and multiplied by 10,000. The log(x+1) method was used to perform natural-log transformation. The top 2,000 highly variable genes (HVGs) in the normalized expression matrix were identified, centered, and scaled before we performed principal component analysis (PCA) based on these HVGs. Batch effects were removed with the Harmony package (v1.0) of R based on the top 50 PCA components identified (22).

Eighteen tumor-normal pairs of osteosarcoma RNA-sequencing data were downloaded from GSE99671. The edgeR package (v 3.30.3) was used to identify up- and down-regulated genes in osteosarcoma (log2FoldChange >1 and adjusted P-value < 0.1). Bulk RNA-seq and clinical data for 80 OS patients from the TARGET database and 53 from GSE21257 were collected. Patients from TARGET were randomly divided into two groups: a training group (40), and an internal verification group (40). The 53 patients from GSE21257 were used as an external training group. All the raw data are available in the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Based on harmony-corrected data, k-nearest neighbors were calculated, and a shared nearest neighbor (SNN) graph was constructed. Then, the modular function was optimized to realize cluster recognition based on the clustering algorithm. The identified clusters were visualized on a 2D map by employing the uniform manifold approximation and projection (UMAP) dimensional reduction technique. We first used the “SingleR” package (R software) as an auxiliary tool to identify cells, and then used the Wilcoxon rank-sum test with Bonferroni correction to identify marker genes for double annotation of cell clusters.

Crosstalk between OS cells and other cells within the tumor microenvironment was calculated with the iTALK package (v0.1.0) in R (23). The top 50% of highly expressed genes were uploaded to the ligand-receptor database to determine their locations within an intercellular signaling network.

Gene set enrichment analysis (GSEA) was used to detect whether the PI3K-AKT signaling pathway was significantly enriched in the subpopulations of OS cells. The enrichment score (ES) and statistical significance (nominal P value) of the ES were calculated; a positive ES suggested that PI3K-AKT signaling was activated in a subpopulation, whereas a negative ES suggested that PI3K-AKT signaling was not activated.

Monocle (v2.16.0) aims to resolve cellular transitions during differentiation through pseudotemporal profiling of scRNA-seq data (24). After inputting the cell-gene matrix into the “newCellDataSet” function with its clustering information, it was computed into a lower dimensional space based on the discriminative dimensionality reduction with trees (DDRTree) method, a more recent manifold learning algorithm, and then OS cells were ordered according to pseudotime.

To establish a generalized linear model and minimize the risk of overfitting, we applied the least absolute shrinkage and selection operator (LASSO) penalty in the COX proportional hazards regression model (“glmnet” package, R software). In the modelling process, “lambda” is the variable and causes the gene coefficients to converge. As lambda increases, the coefficients for individual genes gradually approach zero. When the coefficient becomes zero, it means that the expression of the gene no longer affects the risk score and the gene can be removed from the model. We used 10-fold cross-validation to calculate partial-likelihood, a key metric for evaluating model performance. Value of lambda that gives minimum partial-likelihood is optimal, at which the model constructed of candidate genes with nonzero regression coefficients has the best performance. Risk scores for each patient were calculated using the following formula:

where RS is the risk score for OS patients; N is the number of genes in the model; and the ith gene expression and coefficient are denoted by e _i and c _i respectively. The independent variables are the expression of the modelling genes for each patient, and the dependent variable is the risk score.

According to the median scores for each group, the patients were divided into high-risk and low-risk subgroups.

The human OS cell line HOS was kindly provided by Procell Life Science & Technology Co., Ltd (Wuhan, China) and cultured in Dulbecco’s Modified Eagle Medium/Nutrient Mixture F-12 (DMEM/F-12; Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (FBS; HyClone, South Logan, UT, USA) and 1% penicillin-streptomycin (Gibco, Grand Island, NY, USA) at 37°C with 5% CO₂.

OS cells were grown to 60-80% confluence. In all, 30 pmol EFEMP2 siRNA, GALNT14 siRNA or negative control siRNA (GenePharma, Shanghai, China) was diluted to 150 µl Opti-MEM^® Reduced-Serum Medium (Gibco, Grand Island, NY, USA) respectively and mixed with the diluted Lipofectamine^® RNAiMAX Reagent (1:1 ratio; Thermo Fisher Scientific, Waltham, MA, USA). The siRNA-lipid complex was incubated for 5 minutes at room temperature and then added to the cells.

Total RNA was extracted using the SimplyP RNA Kit (Bioer, Hangzhou, China), and then reversed transcription was carried out using the PrimeScript RT reagent Kit (TaKaRa, Tokyo, Japan). Quantitative RT-PCR was performed with SsoAdvanced Universal SYBR Green Supermix reagent (Bio-Rad, Hercules, CA, USA) on a PCR system (Bio-Rad, Hercules, CA, USA) per the manufacturer’s instructions. The primers were designed and synthesized by Sangon Biotech (Shanghai, China), and the sense and antisense primers are listed in Table S2 .

HOS cells (5x10⁵) were seeded in 6-well plates and transfected by siRNA described above. When the cells reached 90-95% confluence as a monolayer, a scratch wound was made as a marked line in the cultures using a sterile pipette tip. Wells were washed three times with phosphate-buffered saline (PBS; 0.01 M) to remove loose or dead cells, and fresh serum-free medium was added. The cells were photographed at 0-, 12- and 24-hour time points after treatment. The wound area was quantitated using ImageJ software (25).

Bilateral tests were performed for statistical tests. A P value less than 0.05 was considered statistically significant. R software version 4.0.0 (https://www.r-project.org/) was used for analysis. Some R packages, including “dplyr”, “Seurat”, “cowplot”, “Matrix”, “ggplot2”, “reshape2”, “iTALK”, “monocle”, “harmony”, and “survival”, were used in this study.

A total of 11 OS patients with scRNA-seq data were involved in this study (seven primary patients: P01-07, two recurrent patients: P08-09, and two lung-metastatic patients: P10-11) ( Figure 1A ). After quality control and removal of the batch effect between samples, we used the UMAP technique to reduce dimensionality, and graphed the output on a 2D scatter plot ( Figure 1B ). A total of 109,415 single cells were unbiasedly clustered into 14 major identities, and cluster-specific markers were used to annotate cell types: OS cells (COL1A1, CXCL12, MEPE and COL2A1, including osteoblastic OS cells and chondroblastic OS cells); macrophages (HLA-DRA, CD68 and CD14); mesenchymal stem cells (MSCs) (NES, HMGB2 and CCNB2); T cells (CCL5 and CD69); endothelial cells (CD34 and PECAM1); pericytes (ACTA2 and RGS5); B cells (JCHAIN and MZB1) and myoblasts (MYOG) ( Figures 1B–D ). Detailed annotation information and references are given in the Table S1 .

To investigate and visualize the cell-cell interaction network in the OS microenvironment, the iTALK package of R was used to analyze the ligand-receptor interactions in different cell types (23). Based on the primary function of the ligand, iTALK further classifies intercellular crosstalk into 4 categories: checkpoint, cytokine, growth factor and other. In terms of the checkpoint category, OS cells of clusters 1 and 6 expressed higher levels of CD24, which could interact with the SIGLEC10 receptor of cluster 0 macrophages to regulate the antitumor immune response of macrophages and promote immune evasion (26). The cytokine category showed internal connections among macrophages. CCL3, CCL4 and CCL3L1 secreted by cluster 0 and 12 macrophages could interact with CCR1, which was highly expressed on the surface of cluster 4 macrophages, promoted their migration, and regulated the inflammatory response (27). Regarding the growth factor category, we also found that cluster 6 OS cells likely produce VEGFA, which binds to ITGB1 receptors on endothelial cells (ECs) to regulate their migration and promote new vessel formation and tumor metastasis (28–30). ( Figure 2 )

To study the heterogeneity among OS cells and find new cell subpopulations closely related to OS metastasis, data for 53,441 OS cells were extracted, and the cells were clustered into nine subclusters (named Osteosarcoma_0 ~ Osteosarcoma_8), and Osteosarcoma_3 and 8 were chondroblastic and the rest were osteoblastic cells ( Figure 3A and Figure S1 ). As illustrated in Figures 3B, C , the genes CADM1, LINC00662, and RUVBL1, which have been shown to inhibit OS metastasis (31–33), were highly expressed by the cells in subclusters 6 and 7. In addition, Osteosarcoma_8 cells showed high expression levels of COL6A1, COL6A3 and MIF, and these genes have been confirmed to be associated with metastasis (34–37). In addition, many genes that play important roles in osteosarcoma display distinct expression patterns between Osteosarcoma_6, 7 and 8 ( Figure S2A ).

PI3K-AKT signaling is one of the key pathways promoting the metastasis in OS, and inhibition of PI3K-AKT signaling yields enhanced antitumor activity, impaired tumor growth, and reduced migratory and metastatic potential (38, 39). The GSEA results revealed that the PI3K-AKT pathway was activated in Osteosarcoma_8 (ES=0.35, P-value=0.004) but was inhibited in Osteosarcoma_6 (ES=-0.40, P-value=0.004) and Osteosarcoma_7 (ES=-0.39, P-value=0.001) ( Figure 3D ). We further constructed an osteosarcoma-related gene set using the results of the edgeR package and performed GSEA enrichment analysis. Similarly, the results also confirmed that there were significant differences in cell states between Osteosarcoma _6, 7 and 8 ( Figure S2B ). Thus, we demonstrated that subcluster 8 contained aggressive cells, while the cells in subclusters 6 and 7 showed limited associations with metastasis.

Next, we explored the dynamics and regulators of cell fate decisions in OS cells with Monocle 2. The pseudotime trajectory was visualized and is shown in Figure 3E . Strikingly, we noticed that Osteosarcoma_4 was mainly located at the beginning of the trajectory path, with high expression of many stem and progenitor cell marker genes (CD34, Oct-4, NGFR, SSEA3 and Nanog) (40–44). A set of stem cell-related genes were collected from the “CellMarker” database (http://biocc.hrbmu.edu.cn/CellMarker/index.jsp), and their expression gradually decreased with the cell trajectory ( Figure 3F ). Therefore, Osteosarcoma_4 was identified as the principal progenitor. As cells progressed along a differentiation trajectory, they diverged along two separate paths, which represented two distinctive cell fates ( Figures 3E–G ). Interestingly, Osteosarcoma_8, which represented the more aggressive subpopulation, mainly followed the right differentiation trajectory after the branch point. In contrast, the metastatically indolent OS cells in Osteosarcoma_6 and Osteosarcoma_7 predominantly followed the left trajectory ( Figure 3H ). In summary, Osteosarcoma_8 OS cells showed different differentiation states from the cells of Osteosarcoma_6 and Osteosarcoma_7, which may be the reason for their different levels of invasiveness.

Differential analysis was performed between weakly metastatic subclusters (Osteosarcoma_6 and Osteosarcoma_7) and the highly invasive subcluster (Osteosarcoma_8) of OS cells, and 813 metastasis-related genes were identified. As described in Materials and Methods, BCDEG, a model considering five metastasis-related genes (GALNT14, BNIP3, EFEMP2, CPE and DKK1), was developed ( Figures 4A, B ). Previous studies have confirmed that silencing BNIP3 and DKK1 can effectively inhibit OS proliferation and promote apoptosis (45–47) and that CPE overexpression enhances OS growth, migration, and invasion (48, 49). Since the roles of EFEMP2 and GALNT14 in OS have not been explored, which is worthy of further study. We used siRNA to significantly suppress the expression of EFEMP2 and GALNT14 ( Figures 4C, D ). Wound healing analysis showed that, within 24 hours, the scratches in the negative control group were close to healing, while those in the siRNA-EFEMP2 and the siRNA-GALNT14 groups were still relatively obvious ( Figures 4E ). These findings suggest that inhibiting EFEMP2 and GALNT14 could decrease the invasion and metastasis ability of OS.

The BCDEG model was applied to predict survival in the training group, and the results of Kaplan-Meier analysis showed that the risk score was significantly and negatively associated with the prognosis of OS patients ( Figure 4F ). When patients were divided into groups according to risk score, the mortality rate was significantly higher in the high-risk group than the low-risk group, and almost all patients with the highest risk scores died with a very short survival period ( Figure 4G , P=0.001).

The remaining patients were used to verify the BCDEG model, and similar results were observed in the internal ( Figure 5A ) and external validation groups ( Figure 5B ). Time-dependent ROC curves were produced to assess the prognostic ability of BCDEG, and the AUCs at 5 years were 0.838, 0.724, and 0.729 in the training, internal validation, and external validation groups, respectively ( Figure 5C ), indicating that the BCDEG model had a great ability to predict the survival of OS patients. We also wondered whether BCDEG could predict OS metastasis. We generated a box plot, which showed that the risk scores of metastatic patients were significantly higher than those of localized OS patients ( Figure 5D ), and the sensitivity and specificity of prediction were 0.81 and 0.74, respectively.

We also applied the BCDEG model for patients with liver hepatocellular carcinoma and colon adenocarcinoma from The Cancer Genome Atlas (TCGA). The results of TCGA-LIHC and TCGA-COAD stand in sharp contrast to OS patients, as that there were no significant differences between the survival of high-risk and low-risk groups ( Figure 5E ), indicating that this model is specific for OS.

To enable clinicians to have convenient access to this model, an open-source web app was generated. Users can calculate multiple patient risk scores online to evaluate prognosis and guide clinical treatment. The web app can be accessed at http://39.108.85.1:8072/common/home ( Figure 5F ).