- 1Science Research Center, China-Japan Union Hospital of Jilin University, Changchun, China
- 2Department of Geriatrics, China-Japan Union Hospital of Jilin University, Changchun, China
- 3Department of Gastrointestinal Colorectal and Anal Surgery, China-Japan Union Hospital of Jilin University, Changchun, China
Background: Osteosarcoma (OS) is a common primary tumor with extensive heterogeneity. In this study, we used single-cell RNA sequencing (scRNA-seq) and network pharmacology to analyze effective targets for Osteosarcoma treatment.
Methods: The cell heterogeneity of the Osteosarcoma single-cell dataset GSE162454 was analyzed using the Seurat package. The bulk-RNA transcriptome dataset GSE36001 was downloaded and analyzed using the CIBERSORT algorithm. The key targets for OS therapy were determined using Pearson’s correlation analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on key targets. The DeepDR algorithm was used to predict potential drugs for Osteosarcoma treatment. Molecular docking analysis was performed to verify the binding abilities of the predicted drugs and key targets. qRT-PCR assay was used to detect the expression of key targets in osteoblasts and OS cells.
Results: A total of 21 cell clusters were obtained based on the GSE162454 dataset, which were labeled as eight cell types by marker gene tagging. Four cell types (B cells, cancer-associated fibroblasts (CAFs), endothelial cells, and plasmocytes) were identified in Osteosarcoma and normal tissues, based on differences in cell abundance. In total, 17 key targets were identified by Pearson’s correlation analysis. GO and KEGG analysis showed that these 17 genes were associated with immune regulation pathways. Molecular docking analysis showed that RUNX2, OMD, and CD4 all bound well to vincristine, dexamethasone, and vinblastine. The expression of CD4, OMD, and JUN was decreased in Osteosarcoma cells compared with osteoblasts, whereas RUNX2 and COL9A3 expression was increased.
Conclusion: We identified five key targets (CD4, RUNX2, OMD, COL9A3, and JUN) that are associated with Osteosarcoma progression. Vincristine, dexamethasone, and vinblastine may form a promising drug–target pair with RUNX2, OMD, and CD4 for Osteosarcoma treatment.
1 Introduction
Osteosarcoma (OS) is a malignant tumor that occurs mostly in children and adolescents, and accounts for 20% of primary bone tumors worldwide (Zhu et al., 2019; Zhao et al., 2021). Currently, the methods of OS treatment mainly include surgery, chemotherapy, and radiotherapy (Rothzerg et al., 2022). However, current treatment methods for OS are unsatisfactory, with an overall 5 year survival rate of 65%–70% (Rothzerg et al., 2021). Therefore, there is an urgent need to develop novel treatment options for OS.
Gene mutations are considered the underlying cause of OS (Kuijjer et al., 2013). With the development of bioinformatics analyses, an increasing number of genes have been shown to be involved in the development of OS. For example, mutant p53 promotes cell proliferation, migration, and chemoresistance in OS (Tang et al., 2019a). LIM kinase one overexpression contributes to metastasis, invasion, and multidrug resistance in OS (Yang et al., 2018). Overexpression of Notch homolog protein three is associated with metastasis and poor prognosis in patients with OS (Tang et al., 2019b). All evidence suggests that bioinformatic analysis may provide valuable clues for the treatment of OS.
Single-cell sequencing (scRNA-seq) is a new bioinformatic analysis technique that fills the gap in other bioinformatic techniques for single-cell studies. Moreover, scRNA-seq technology has been prominent in exploring tumor heterogeneity and providing new therapeutic leads for the treatment of many cancers, including pancreatic cancer (Han et al., 2021), gastric cancer (Jiang et al., 2022), and Ewing sarcoma (Aynaud et al., 2020). Tumor heterogeneity plays an important role in cancer progression. Therefore, understanding the gene expression patterns of individual cells is particularly important in cancer treatment (Zhang et al., 2021a). OS exhibits great tumor heterogeneity; thus, we explored the potential targets for the clinical diagnosis and treatment of patients with OS using scRNA-seq and network pharmacology analysis.
2 Materials and methods
2.1 RNA sequencing (RNA-seq) data download and analysis
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.
2.2 Identification of key targets in OS
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).
2.3 Protein-protein interaction network construction
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).
2.4 Functional enrichment analysis
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).
2.5 OS therapeutic drug prediction
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.
2.6 Drug–target interaction prediction
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.
2.7 Molecular docking analysis
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).
2.8 Cell culture
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).
2.9 Quantitative real-time PCR (qRT-PCR)
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.
2.10 Statistics analysis
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.
3 Results
3.1 Identification of cell clusters and dimension reduction analysis
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).
FIGURE 1. Identification of cell clusters and dimension reduction analysis. (A) Twenty-one cell clusters of dimension reduction analysis; (B) Ten cell types were identified by marker genes; (C) The expression of marker genes in eight cell types; (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for marker genes of eight cell types.
3.2 Collection of differential cell types
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).
FIGURE 2. Collection of differential cell types. (A) The CIBERSORT algorithm was used to obtain the differential abundance of eight cell types in OS and normal tissues from the GSE36001 dataset; (B) Scatter plot of the correlation between the abundance of the four cell types (B cells, cancer-associated fibroblasts (CAFs), endothelial cells, and plasmocytes) and the expression of immune checkpoints TDO2, PDCD1, LGALS9, and PVR.
3.3 Screening out key targets for OS
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).
FIGURE 3. Screening out key targets for OS. (A) A Venn diagram of OS-related targets and four different cell types marker genes; (B) Heat map of 17 potential targets correlated with cell abundance in B cells, CAFs, endothelial cells, and plasmocytes. *p <.05, **p <.01.
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:
FIGURE 4. Prognostic analysis of five key targets (CD4, RUNX2, OMD, COL9A3, and JUN) in OS. (A) The survival analysis of these five genes. (B) The receiver operator characteristic (ROC) curve of these five genes. (C) The expression of immune checkpoints BTLA and PDL1 (CD274) were measured in the high- and low-risk groups. *p <.05.
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).
3.4 PPI network construction and functional enrichment analysis
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).
FIGURE 5. Protein-protein interaction (PPI) network of 17 key targets. The line from purple to red indicates a higher score of the reciprocal relationship, and the bigger size of the node indicates the higher degree.
FIGURE 6. Function enrichment analysis of 17 key targets. (A) Top five terms of Gene Ontology (GO) enrichment analysis. GO consists of biological processes (BP), cellular component (CC), molecular function (MF); (B) Top 15 pathways of Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis.
3.5 Drug–target interaction prediction
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).
FIGURE 7. Drug-target interaction score dispersion points. The horizontal coordinate is the action score, and the vertical coordinate is the drug-target pair. The green dashed line is the 75% quartile of the score.
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).
FIGURE 8. Molecular docking analysis of vincristine, dexamethasone, and vinblastine with RUNX2, OMD, and CD4.
3.6 Cell validation assay
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).
FIGURE 9. The relative expression of CD4, OMD, JUN, COL9A3, and RUNX2 was detected using qRT-PCR assay. *p <.05 and **p <.01 compared with the hFOB1.19 group.
4 Discussion
Over the past decades, traditional bioinformatics techniques have partially revealed the pathological mechanisms of OS; however, further research is still required. scRNA-seq technology can reveal the development of disease at the individual cell level. Therefore, in this study, we used scRNA-seq technology combined with network pharmacology analysis to explore effective therapeutic targets for OS treatment.
Tumors are heterogeneous cell populations that contain transformed cancer, supporting, and tumor-infiltrating cells (Prasetyanti and Medema, 2017). In this study, scRNA-seq with an unbiased approach was used to characterize cellular heterogeneity in OS. We identified the abundance of B cells, CAFs, endothelial cells, and plasmocytes has correction with the expression of immune checkpoints TDO2, PDCD1, LGALS9, and PVR. B cells are the main component of humoral immunity. It has been suggested that an increase in B cells is associated with good prognosis of OS (Li et al., 2021). CAFs are activated fibroblasts present within the tumor microenvironment that promote tumor cell growth, invasion, metastasis, and drug resistance (Zhang et al., 2021b). Based on these data, we confirmed that OS exhibits cell heterogeneity, and the differential abundance of cell types could affect the malignant progression of OS.
A total of 17 key targets were obtained from the intersection of marker genes of the four cell types and OS-related targets using Pearson correction analysis. Network pharmacology analysis was performed to determine the functions of the 17 key targets. GO analysis showed that these targets were involved in the regulation of ERK1 and ERK2 cascade, MAPK cascade, and chemotaxis. Reportedly, ERK signaling is involved in various cellular progress, such as proliferation, migration, and differentiation (Bonjardim, 2017). A study found that blocking ERK1/2 signaling pathway inhibits OS cell growth and metastasis (Yuan et al., 2020). Another study found that blocking the nuclear translocation of phosphor-ERK suppresses the migration and invasion of OS cells (Kim et al., 2021). MAPK signaling has a major effect on cell survival and apoptosis (Chen et al., 2016). It has been reported that activating the MAPK signaling pathway induces cell death in human OS (Lv et al., 2020). Zhou et al. (2020) found that inhibiting the function of TLR4-mediated MAPK-NF-ĸB signaling pathway against the oncogenesis of OS.
KEGG analysis showed that these targets act against OS via the MAPK signaling pathway, PI3K-Akt signaling pathway, human T-cell leukemia virus 1 infection, and so on. CD4-positive T cell are involved in the tumor immune environment in OS (Liu et al., 2020). T follicular helper (Tfh) cells are a subpopulation of CD4-positive T cell that may play an important role in the tumor microenvironment (Ochando and Braza, 2017). Study has shown that inhibition of the PI3K/Akt/mTOR pathway enhances the ability of OS Tfh cells to promote B cell maturation and immune function (Jiang et al., 2021). RUNX2 is an important transcription factor for bone development and osteoblast differentiation, and both metastasis and chemoresistance are associated with dysregulation of RUNX2 in OS (Vega et al., 2017). A review reported that mutual activation of the PI3K/Akt pathway and RUNX2 may be one of the main drivers of tumor progression or migration (Cohen-Solal et al., 2015). JUN is a factor of the JNK pathway, JNK pathway activation induces apoptosis and autophagy of OS cells (Li et al., 2020). In our study, qRT-PCR assay exhibited 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. Moreover, multi-factor Cox regression algorithm and ROC curve analyses showed that five key genes are effective in the diagnosis and prognosis of OS. Furthermore, molecular docking results showed that vincristine, dexamethasone, and vinblastine all bound well to the key targets RUNX2, OMD, and CD4. Vincristine, dexamethasone, and vinblastine are common drugs used to treat OS (Boscoboinik et al., 1994; Ahlström et al., 2005; Chen et al., 2019). These results suggested that these five key targets could be potential targets for OS treatment.
In conclusion, four cell types (B cells, CAFs, endothelial cells, and plasmocytes) were identified in OS and normal tissues, based on differences in cell abundance. We identified five key targets (CD4, RUNX2, OMD, COL9A3, and JUN) that are associated with OS progression. GO and KEGG analysis showed that these five genes were associated with immune regulation pathways. Vincristine, dexamethasone, and vinblastine may form a promising drug–target pair with RUNX2, OMD, and CD4 for OS treatment. However, there are still shortcomings in our study. For example, we have only carried out a simple experimental validation of the bioinformatics results and our results lack the support of clinical results. Our study identified potential targets for the treatment of OS.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
YW: Conceptualization, validation, investigation, writing—original draft. DQ: Methodology, software, validation. YG: Formal analysis, isualization. YZ: Resources, methodology, data curation. YL: Methodology, software. LH: Validation, writing—Review and editing, funding acquisition. All authors reviewed the results and approved the final version of the manuscript.
Funding
This work was supported by Natural Science Foundation of Science and Technology Department of Jilin Province (grant number 20210101266JC) and Development and Reform Commission of Jilin Province (grant number 2020C034-3).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.1098800/full#supplementary-material
References
Ahlström, M., Pekkinen, M., Huttunen, M., and Lamberg-Allardt, C. (2005). Dexamethasone down-regulates cAMP-phosphodiesterase in human osteosarcoma cells. Biochem. Pharmacol. 69 (2), 267–275. doi:10.1016/j.bcp.2004.09.012
Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol. 20 (2), 163–172. doi:10.1038/s41590-018-0276-y
Aynaud, M. M., Mirabeau, O., Gruel, N., Grossetête, S., Boeva, V., Durand, S., et al. (2020). Transcriptional programs define intratumoral heterogeneity of ewing sarcoma at single-cell resolution. Cell Rep. 30 (6), 1767–1779. e6. doi:10.1016/j.celrep.2020.01.049
Bonjardim, C. A. (2017). Viral exploitation of the MEK/ERK pathway - a tale of vaccinia virus and other viruses. Virology 507, 267–275. doi:10.1016/j.virol.2016.12.011
Boscoboinik, D., Galeotti, T., and Azzi, A. (1994). Vinblastine-dependent down-modulation of TNF receptors in human osteosarcoma cells is mediated by protein kinase C activity. Biochem. Biophys. Res. Commun. 199 (1), 374–379. doi:10.1006/bbrc.1994.1239
Burley, S. K., Berman, H. M., Kleywegt, G. J., Markley, J. L., Nakamura, H., and Velankar, S. (2017). Protein Data Bank (PDB): The single global macromolecular structure archive. Methods Mol. Biol. Clift. NJ) 1607, 627–641. doi:10.1007/978-1-4939-7000-1_26
Chen, J. C., Hsieh, M. J., Chen, C. J., Lin, J. T., Lo, Y. S., Chuang, Y. C., et al. (2016). Polyphyllin G induce apoptosis and autophagy in human nasopharyngeal cancer cells by modulation of AKT and mitogen-activated protein kinase pathways in vitro and in vivo. Oncotarget 7 (43), 70276–70289. doi:10.18632/oncotarget.11839
Chen, R., Huang, L. H., Gao, Y. Y., Yang, J. Z., and Wang, Y. (2019). Identification of differentially expressed genes in MG63 osteosarcoma cells with drug-resistance by microarray analysis. Mol. Med. Rep. 19 (3), 1571–1580. doi:10.3892/mmr.2018.9774
Cohen-Solal, K. A., Boregowda, R. K., and Lasfar, A. (2015). RUNX2 and the PI3K/AKT axis reciprocal activation as a driving force for tumor progression. Mol. cancer 14, 137. doi:10.1186/s12943-015-0404-3
Cui, J., Tian, Y., Liu, T., Lin, X., Li, L., Li, Z., et al. (2022). Pancancer analysis of revealed TDO2 as a biomarker of prognosis and immunotherapy. Dis. markers 2022, 5447017. doi:10.1155/2022/5447017
Fan, Y., Li, T., Xu, L., and Kuang, T. (2020). Comprehensive analysis of immunoinhibitors identifies LGALS9 and TGFBR1 as potential prognostic biomarkers for pancreatic cancer. Comput. Math. methods Med. 2020, 6138039. doi:10.1155/2020/6138039
Feng, C., Zhao, M., Jiang, L., Hu, Z., and Fan, X. (2021). Mechanism of modified danggui sini decoction for knee osteoarthritis based on network pharmacology and molecular docking. Evid. Based Complement. Altern. Med. 2021, 6680637. doi:10.1155/2021/6680637
Fisher, L. D., and Lin, D. Y. (1999). Time-dependent covariates in the Cox proportional-hazards regression model. Annu. Rev. public health 20, 145–157. doi:10.1146/annurev.publhealth.20.1.145
Han, J., DePinho, R. A., and Maitra, A. (2021). Single-cell RNA sequencing in pancreatic cancer. Nat. Rev. Gastroenterology hepatology 18 (7), 451–452. doi:10.1038/s41575-021-00471-z
Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., Zheng, S., Butler, A., et al. (2021). Integrated analysis of multimodal single-cell data. Cell 184 (13), 3573–3587.e29. e29. doi:10.1016/j.cell.2021.04.048
Huang, K., Fu, T., Glass, L. M., Zitnik, M., Xiao, C., and Sun, J. (2021). DeepPurpose: A deep learning library for drug-target interaction prediction. Bioinforma. Oxf. Engl. 36 (22-23), 5545–5547. doi:10.1093/bioinformatics/btaa1005
Jiang, B., Kang, X., Zhao, G., Lu, J., and Wang, Z. (2021). miR-138 reduces the dysfunction of T follicular helper cells in osteosarcoma via the PI3K/Akt/mTOR pathway by targeting PDK1. Comput. Math. methods Med. 2021, 2895893. doi:10.1155/2021/2895893
Jiang, H., Yu, D., Yang, P., Guo, R., Kong, M., Gao, Y., et al. (2022). Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing. Clin. Transl. Med. 12 (2), e730. doi:10.1002/ctm2.730
Kim, J. Y., Park, H. H., Yong, T. S., and Jeon, S. H. (2021). Lithium chloride inhibits the migration and invasion of osteosarcoma cells by blocking nuclear translocation of phospho-Erk. Biochem. Biophys. Res. Commun. 581, 74–80. doi:10.1016/j.bbrc.2021.10.025
Kuijjer, M. L., Hogendoorn, P. C., and Cleton-Jansen, A. M. (2013). Genome-wide analyses on high-grade osteosarcoma: Making sense of a genomically most unstable tumor. Int. J. cancer 133 (11), 2512–2521. doi:10.1002/ijc.28124
Li, B., Zhou, P., Xu, K., Chen, T., Jiao, J., Wei, H., et al. (2020). Metformin induces cell cycle arrest, apoptosis and autophagy through ROS/JNK signaling pathway in human osteosarcoma. Int. J. Biol. Sci. 16 (1), 74–84. doi:10.7150/ijbs.33787
Li, G. Q., Wang, Y. K., Zhou, H., Jin, L. G., Wang, C. Y., Albahde, M., et al. (2021). Application of immune infiltration signature and machine learning model in the differential diagnosis and prognosis of bone-related malignancies. Front. Cell Dev. Biol. 9, 630355. doi:10.3389/fcell.2021.630355
Liu, Z., Wen, J., Wu, C., Hu, C., Wang, J., Bao, Q., et al. (2020). MicroRNA-200a induces immunosuppression by promoting PTEN-mediated PD-L1 upregulation in osteosarcoma. Aging 12 (2), 1213–1236. doi:10.18632/aging.102679
Lv, H., Zhen, C., Liu, J., and Shang, P. (2020). β-Phenethyl isothiocyanate induces cell death in human osteosarcoma through altering iron metabolism, disturbing the redox balance, and activating the MAPK signaling pathway. Oxidative Med. Cell. Longev. 2020, 5021983. doi:10.1155/2020/5021983
Miao, Y., Wang, J., Li, Q., Quan, W., Wang, Y., Li, C., et al. (2020). Prognostic value and immunological role of PDCD1 gene in pan-cancer. Int. Immunopharmacol. 89, 107080. doi:10.1016/j.intimp.2020.107080
Morad, G., Helmink, B. A., Sharma, P., and Wargo, J. A. (2021). Hallmarks of response, resistance, and toxicity to immune checkpoint blockade. Cell 184 (21), 5309–5337. doi:10.1016/j.cell.2021.09.020
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. methods 12 (5), 453–457. doi:10.1038/nmeth.3337
Ochando, J., and Braza, M. S. (2017). T follicular helper cells: A potential therapeutic target in follicular lymphoma. Oncotarget 8 (67), 112116–112131. doi:10.18632/oncotarget.22788
Piñero, J., Ramírez-Anguita, J. M., Saüch-Pitarch, J., Ronzano, F., Centeno, E., Sanz, F., et al. (2020). The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res. 48 (D1), D845-d55. doi:10.1093/nar/gkz1021
Prasetyanti, P. R., and Medema, J. P. (2017). Intra-tumor heterogeneity from a cancer stem cell perspective. Mol. cancer 16 (1), 41. doi:10.1186/s12943-017-0600-4
Rothzerg, E., Ingley, E., Mullin, B., Xue, W., Wood, D., and Xu, J. (2021). The Hippo in the room: Targeting the Hippo signalling pathway for osteosarcoma therapies. J. Cell. physiology 236 (3), 1606–1615. doi:10.1002/jcp.29967
Rothzerg, E., Pfaff, A. L., and Koks, S. (2022). Innovative approaches for treatment of osteosarcoma. Exp. Biol. Med. (Maywood, NJ) 247 (4), 310–316. doi:10.1177/15353702211067718
Safran, M., Dalah, I., Alexander, J., Rosen, N., Iny Stein, T., Shmoish, M., et al. (2010). GeneCards version 3: The human gene integrator. Database J. Biol. databases curation 2010, baq020. doi:10.1093/database/baq020
Seeliger, D., and de Groot, B. L. (2010). Ligand docking and binding site analysis with PyMOL and Autodock/Vina. J. computer-aided Mol. Des. 24 (5), 417–422. doi:10.1007/s10822-010-9352-6
Stamm, H., Wellbrock, J., and Fiedler, W. (2018). Interaction of PVR/PVRL2 with TIGIT/DNAM-1 as a novel immune checkpoint axis and therapeutic target in cancer. Mammalian genome official J. Int. Mammalian Genome Soc. 29 (11-12), 694–702. doi:10.1007/s00335-018-9770-7
Sterling, T., and Irwin, J. J. (2015). ZINC 15-ligand discovery for everyone. J. Chem. Inf. Model 55 (11), 2324–2337. doi:10.1021/acs.jcim.5b00559
Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING database in 2021: Customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49 (D1), D605–D612. doi:10.1093/nar/gkaa1074
Tang, F., Min, L., Seebacher, N. A., Li, X., Zhou, Y., Hornicek, F. J., et al. (2019). Targeting mutant TP53 as a potential therapeutic strategy for the treatment of osteosarcoma. J. Orthop. Res. official Publ. Orthop. Res. Soc. 37 (3), 789–798. doi:10.1002/jor.24227
Tang, X. F., Cao, Y., Peng, D. B., Zhao, G. S., Zeng, Y., Gao, Z. R., et al. (2019). Overexpression of Notch3 is associated with metastasis and poor prognosis in osteosarcoma patients. Cancer Manag. Res. 11, 547–559. doi:10.2147/CMAR.S185495
Vega, O. A., Lucero, C. M. J., Araya, H. F., Jerez, S., Tapia, J. C., Antonelli, M., et al. (2017). Wnt/β-Catenin signaling activates expression of the bone-related transcription factor RUNX2 in select human osteosarcoma cell types. J. Cell Biochem. 118 (11), 3662–3674. doi:10.1002/jcb.26011
Yang, J-Z., Huang, L-H., Chen, R., Meng, L-J., Gao, Y-Y., Ji, Q-Y., et al. (2018). LIM kinase 1 serves an important role in the multidrug resistance of osteosarcoma cells. Oncol. Lett. 15 (1), 250–256. doi:10.3892/ol.2017.7317
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics a J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118
Yuan, X. H., Zhang, P., Yu, T. T., Huang, H. K., Zhang, L. L., Yang, C. M., et al. (2020). Lycorine inhibits tumor growth of human osteosarcoma cells by blocking Wnt/β-catenin, ERK1/2/MAPK and PI3K/AKT signaling pathway. Am. J. Transl. Res. 12 (9), 5381–5398.
Zeng, X., Zhu, S., Liu, X., Zhou, Y., Nussinov, R., and Cheng, F. (2019). deepDR: a network-based deep learning approach to in silico drug repositioning. Bioinforma. Oxf. Engl. 35 (24), 5191–5198. doi:10.1093/bioinformatics/btz418
Zhang, Y., Liu, Z., Yang, X., Lu, W., Chen, Y., Lin, Y., et al. (2021). H3K27 acetylation activated-COL6A1 promotes osteosarcoma lung metastasis by repressing STAT1 and activating pulmonary cancer-associated fibroblasts. Theranostics 11 (3), 1473–1492. doi:10.7150/thno.51245
Zhang, Y., Wang, D., Peng, M., Tang, L., Ouyang, J., Xiong, F., et al. (2021). Single-cell RNA sequencing in cancer research. J. Exp. Clin. cancer Res. CR 40 (1), 81. doi:10.1186/s13046-021-01874-1
Zhao, X., Wu, Q., Gong, X., Liu, J., and Ma, Y. (2021). Osteosarcoma: A review of current and future therapeutic approaches. Biomed. Eng. online 20 (1), 24. doi:10.1186/s12938-021-00860-0
Zhou, J., Liu, Q., Qian, R., Liu, S., Hu, W., and Liu, Z. (2020). Paeonol antagonizes oncogenesis of osteosarcoma by inhibiting the function of TLR4/MAPK/NF-κB pathway. Acta Histochem. 122 (1), 151455. doi:10.1016/j.acthis.2019.151455
Keywords: osteosarcoma, single-cell RNA sequencing, network pharmacology, molecular docking, therapeutic target
Citation: Wang Y, Qin D, Gao Y, Zhang Y, Liu Y and Huang L (2023) Identification of therapeutic targets for osteosarcoma by integrating single-cell RNA sequencing and network pharmacology. Front. Pharmacol. 13:1098800. doi: 10.3389/fphar.2022.1098800
Received: 15 November 2022; Accepted: 19 December 2022;
Published: 06 January 2023.
Edited by:
Lin Qi, Second Xiangya Hospital, Central South University, ChinaReviewed by:
Xiaoqing Guan, Zhejiang Cancer Hospital, ChinaYancheng Tang, Hong Kong Baptist University, Hong Kong SAR, China
Copyright © 2023 Wang, Qin, Gao, Zhang, Liu and Huang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Lihong Huang, bGhodWFuZ0BqbHUuZWR1LmNu