The single-cell dataset GSE162454 and bulk-RNA transcriptome dataset GSE36001 for OS were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/gds) database.

Seurat 4.0 (Hao et al., 2021) was used for quality control, dimensionality reduction, clustering, and marker gene screening of the single-cell dataset, GSE162454. Cell types were annotated and differential gene analysis was performed using singleR (Aran et al., 2019). Subsequently, the CIBERSORT algorithm (Newman et al., 2015) was used to calculate the abundance of cell types in the OS and normal samples from GSE36001.

OS-related targets were collected in the Genecards (Safran et al., 2010) and DisGeNET (Piñero et al., 2020) databases by searching with the keyword “osteosarcoma.” The search results were then merged and duplicates were deleted. Subsequently, these target genes were intersected with marker genes in the differentially expressed cells. Pearson correlation analysis was performed to screen genes whose expression levels were significantly correlated with differential cell type abundance (p <.05).

The protein-protein interaction (PPI) network of key targets was analyzed using the STRING database (Szklarczyk et al., 2021) and visualized using the R packages ggraph (version 2.0.5) and igraph (version 1.3.1).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of key targets were performed using the R package clusterProfiler (version 4.4.2) (Yu et al., 2012).

In drug-disease association prediction, the deep learning-based algorithm deepDR (Zeng et al., 2019) learns the high-level features of drugs from a heterogeneous network through a multimodal deep autoencoder. Then, the learned low-dimensional representations of drugs and drug-disease pairs are encoded and decoded by the autoencoder to perform drug indication inference and filter the drug-disease pairs with high association based on the association score.

Drug-target interactions (DTI) are used to indicate the strength of the binding ability of a compound to a protein target. The deep-learning algorithm DeepPurpose (Huang et al., 2021) was used to perform DTI prediction using a simplified molecular-input line-entry system (SMILES) of compounds and amino acid sequence pairs of proteins as input data. Drug-target pairs with higher scores were screened based on the prediction scores.

The crystal structures of key target proteins were downloaded from the RCSB Protein Data Bank (http://www.pdb.org/) (Burley et al., 2017). Protein conformations were modified using PyMOL and AutoDock 1.5.6, including the removal of the original ligands and water, addition of hydrogen, optimization of amino acids, and calculation of charges (Seeliger and de Groot, 2010). The structure file of the drug in “mol2” format was downloaded via ZINC (https://zinc.docking.org/) (Sterling and Irwin, 2015). The downloaded protein and drug files were then converted to PDBQT format using Open Babel GUI software. Finally, molecular docking was performed using AutoDock 1.5.6, and the results were visualized using PyMOL. The screening criteria were binding energy less than −5.0 kcal/mol and the formation of hydrogen bonds between ligand receptors (Feng et al., 2021).

OS cell lines (U-2 OS, MG-63, and Saos-2) and human osteoblast cell line (hFOB1.19) were purchased from ATCC (VA, United States). OS cells were cultured in Dulbecco’s modified eagle medium (DMEM, Gibco, CA, United States) containing 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (Invitrogen, CA, United States). hFOB1.19 cells were maintained in DMEM/F-12 (Gibco) supplemented with 10% FBS, 2.5 mM L-glutamine, and .3 mg/mL geneticin (Invitrogen).

Total RNAs were extracted with Trizol regent (Tiangen, China) and reverse-transcribed to cDNA using the Reverse Transcription Kit (Promega, China). Then, qRT-PCR was conducted on CFX96 Touch (Bio-Rad, Hercules, CA, United States) system with thermocycling conditions followed as 95°C for 3 min, 40 cycles for 95°C for 12 s and 62°C for 40 s. The relative expression of CD4, RUNX2, OMD, COL9A3, and JUN was calculated with 2^−ΔΔCT method and GAPDH as a housekeeping gene. The primer sequences of these genes are listed in Supplementary Table S1.

Data were expressed as mean ± standard deviation and analyzed by Prism 8.0 software with t-test or one-way ANOVA. p <.05 has a significant difference.

Using scRNA-seq analysis, 21 cell clusters and eight cell types were identified (Figures 1A,B). Marker gene expression in the eight cell types was different between the different cell types (Figure 1C). KEGG enrichment analysis revealed significant heterogeneity in the enriched pathways for marker genes in the eight cell types (Figure 1D).

To further screen for differential cell types in OS, the abundance of eight cell types in OS and normal tissues from the GSE36001 dataset was analyzed using the CIBERSORT algorithm (Figure 2A). Studies show that immune checkpoints TDO2, PDCD1, LGALS9, and PVR play an important role in cancer treatment and prognosis (Stamm et al., 2018; Fan et al., 2020; Miao et al., 2020; Cui et al., 2022). Based on Pearson correlation analysis, cell abundance was screened to be significantly correlated with immune checkpoint expression levels, and the results showed that the abundance of B cells, cancer-associated fibroblasts (CAFs), endothelial cells, and plasmocytes were significantly correlated with immune checkpoints TDO2 (p = .043), PDCD1 (p = .017), LGALS9 (p = .017), and PVR (p = .026), respectively (Figure 2B).

A total of 4,236 OS-related targets were retrieved from the Genecards and DisGeNET databases. Next, overlapping with the marker genes of the four cell types, we obtained 289 common targets (Figure 3A). Finally, further screening by Pearson correlation analysis yielded 17 key targets (GZMB, IL1A, IGFBP4, JUN, KDR, CD2, FLT1, CCR7, CD4, COL9A3, SPRY4, OMD, RUNX2, PTPRG, HSPA1A, SOX18, and PDGFRB) that were significantly associated with the abundance of four differential cells (Figure 3B).

Subsequently, we obtained expression data for 17 key targets corresponding to OS survival time information through the TARGET database (https://ocg.cancer.gov/programs/target). The five targets (CD4, RUNX2, OMD, COL9A3, JUN) may be the important prognostic factor via the multivariate Cox regression analysis. A combined prognostic marker model consisting of these five genes was constructed based on a multivariate Cox regression algorithm (Fisher and Lin, 1999). The risk calculation formula was as follows: R i s k   s c o r e = − 0.31634 × C D 4 + 0.68318 × R U N X 2 − 0.26491 × O M D + 0.16464 × C O L 9 A 3 − 0.25361 × J U N . The results showed that the five key genes were able to accurately grade the risk of OS (p = .0019, Figure 4A). Moreover, the area under the curve (AUC) values of the receiver operator characteristic (ROC) curve were greater than .72 in the 3-, and 5-year survival analyses of OS (Figure 4B).

Additionally, cancer cells can undergo immune escape by dysregulating immune checkpoint proteins (Morad et al., 2021). To verify the accuracy of the high- and low-risk groupings of the five key genes, the expression of the immune checkpoints BTLA and PDL1 was detected. The expression levels of BTLA and PDL1 (CD274) were higher in the low-risk group than in the high-risk group (Figure 4C).

A total of 161 pairs of reciprocal relationships were obtained from the PPI network analysis of 17 key targets using the STRING database (Figure 5). GO and KEGG enrichment analyses were performed to explore the functions of the 17 key targets. GO enrichment analysis revealed 787 enriched terms (p <.05), and the top five terms of biological processes (BP), cellular component (CC), and molecular function (MF) are shown in Figure 6A. BP terms were primarily related to the regulation of ERK1 and ERK2 cascade, MAPK cascade, and chemotaxis. CC terms are located in the transcription regulator complex, RNA polymerase II transcription regulator complex, and external side of plasma membrane. The MF terms are associated with DNA-binding transcription activator activity, RNA polymerase II-specific, growth factor binding, and protein tyrosine kinase activity. Simultaneously, KEGG enrichment analysis was enriched in 77 pathways (p <.05), primarily in the MAPK signaling pathway, PI3K-Akt signaling pathway, and Human T-cell leukemia virus one infection (Figure 6B).

Using the deepDR algorithm, ten drugs with a high score association with OS were obtained (Table 1). The interaction relationship between 17 key targets and 10 drugs was predicted using the DeepPurpose algorithm (Figure 7). The results of screening scores with greater than 75% quartiles yielded six drugs (DB04572, DB01005, DB01234, DB00541, DB00570, and DB00309) that may act on 17 key targets (Table 2).

To confirm that these six drugs were suitable for the treatment of OS, molecular docking analysis of six drugs and CD4, RUNX2, OMD, COL9A3, and JUN were performed using AutoDock Vina. The docking results showed RUNX2 has good binding affinity with vincristine (DB00541), vinblastine (DB00570), and dexamethasone (DB01234); OMD has good binding affinity with vincristine (DB00541), vinblastine (DB00570), and dexamethasone (DB01234); CD4 has good binding affinity with vincristine (DB00541), vinblastine (DB00570), and dexamethasone (DB01234) (Table 3; Figure 8).

To evaluate the prognostic effects of key targets on OS, qRT-PCR assay was used to detect the expression of CD4, RUNX2, OMD, COL9A3, and JUN in OS and osteoblast cells. The results showed that the expression of CD4, OMD, and JUN was decreased in OS cells compared with hFOB1.19 cells, whereas RUNX2 and COL9A3 expression was increased (p < .01; Figure 9).