Single-cell RNA landscape of the intratumoral heterogeneity and expression of angiogenesis-related genes in osteosarcoma

Single-cell RNA landscape of the intratumoral heterogeneity and expression of angiogenesis-related genes in osteosarcoma | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Single-cell RNA landscape of the intratumoral heterogeneity and expression of angiogenesis-related genes in osteosarcoma Hao Li, Tian Ma, Zhiqian Yi, Fang Gao, Xiaojuan Li, Mi Li This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-5305987/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background: Osteosarcoma is an aggressive malignancy of bone that poses significant treatment challenges and has been a focal point of extensive research due to its complex pathogenesis. Despite advances in traditional therapeutic approaches, the intricate genetic and cellular landscape of osteosarcoma remains inadequately understood, emphasizing the need for innovative research methodologies to unravel its underlying mechanisms. Objective: This study aims to leverage the power of single-cell transcriptome sequencing technology to elucidate the cellular heterogeneity, gene expression patterns, intercellular communication networks, and critical genetic pathways implicated in osteosarcoma. By doing so, we intend to contribute valuable insights into the biogenesis of this malignancy, which may ultimately inform precision treatment strategies. Methods: Utilizing single-cell sequencing, we conducted a comprehensive analysis of osteosarcoma samples to identify diverse cellular subpopulations within the tumor microenvironment. Our focus on gene expression profiles revealed significant differences across these subpopulations. Moreover, we employed bioinformatics approaches to explore the intercellular communication networks and identified key ligand-receptor pairings, substantiating the role of angiogenesis-related genes prominently expressed in osteoblasts and their proliferative counterparts. Results: Our findings underscore the critical involvement of angiogenesis in osteosarcoma pathogenesis, with notable pathway activity variations among distinct cellular subpopulations. Additionally, protein interaction network mapping has unveiled significant discrepancies in pathway activities and highlighted the potential functional roles of key genes involved in tumor progression. Conclusion: This study offers a comprehensive exploration of the biological characteristics of osteosarcoma through single-cell sequencing technology, thereby establishing a robust theoretical foundation that may facilitate the development of targeted and effective therapeutic strategies. However, it is essential to recognize that these findings are preliminary, necessitating further validation through expanded sample sizes and integration of multi-omics data. Future research will delve deeper into the mechanisms of the identified key pathways and genes, with the aspiration of enhancing the prognostic outcomes and quality of life for patients with osteosarcoma. Osteosarcoma Single-cell sequencing Cellular heterogeneity Gene expression Intercellular communication Pathway analysis Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Introduction Osteosarcoma, one of the most commonly occurring and notoriously malignant bone tumors, has long been a focal point of biomedical research, particularly in elucidating its mechanisms of occurrence, progression, and metastasis[ 1 ]. Osteosarcoma, arising from primitive mesenchymal-derived osteoblasts, commonly manifests in bones undergoing rapid growth[ 2 ]. The treatment protocol for osteosarcoma involves surgical resection and chemotherapy, whereas radiotherapy is recommended in cases of unresectable osteosarcoma[ 3 ]. It exhibits a high proclivity for local invasion and early metastasis, which poses a threat to patient survival. Despite the development of anti-cancer therapeutics, the overall survival rate of osteosarcoma patients has improved, however, the prognosis remains poor for those with metastatic or recurrent osteosarcoma[ 4 ]. Therefore, researching new treatment strategies and developing new targeted drugs for osteosarcoma are paramount. However, the intricate mechanism of occurrence, development, and metastasis of osteosarcoma remain elusve. In recent years, the evolution of single-cell sequencing technology has remarkable breakthroughs in understanding of tumor tissues heterogeneity and the interaction between tumor cells and their microenvironment[ 5 ]. Significant progress has also been made in the study of osteosarcoma cells using single-cell sequencing technology[ 6 ]. This article aims to delve into the cellular heterogeneity among osteosarcoma subpopulations, angiogenesis-related gene expression patterns, intercellular communication networks, and crucial genes and pathways. Through a comprehensive analysis of osteosarcoma's single-cell data, we seek to provide a noval theoretical framework for understanding its pathogenesis and developing precision treatment strategies. Utilizing single-cell transcriptome sequencing technology, this study conducted a comprehensive analysis of tumor samples from osteosarcoma patients. Through rigorous quality control measures and cell type annotation, we successfully identified multiple cell subpopulations and unveiled gene expression disparities among them. Notably, we uncovered the expression profiles of angiogenesis-related genes in these subpopulations and observed their elevated expression in specific cell types. This finding underscored the crucial role of angiogenesis in osteosarcoma's tumorigenesis and progression. Furthermore, leveraging bioinformatics techniques, we delved into the intercellular communication network and identified key genes. By quantifying the strength of intercellular communication, we discovered that certain cell types occupy pivotal positions within the communication network, serving as a decisive role in intercellular signaling. These insights not only enhanced our understanding of cellular interactions within the osteosarcoma microenvironment, but also presented potential therapeutic targets. Lastly, our analysis revealed significant differences in pathway activity among osteosarcoma cell subpopulations and the pivotal roles of certain genes. This comprehensive analysis offered a deeper understanding of osteosarcoma's biological complexities, paving the way for more targeted and effective therapeutic approaches. Materials and Methods Data download The Cancer Genome Atlas (TCGA)[7] collects various human cancers and tumor subtype data, including clinical information, genomic variations, mRNA expression profiles, miRNA expression data, methylation, and more, which is an important data source for cancer researchers. TCGA also contains some data from TARGET. For our study, we retrieved reliable mRNA expression data in FPKM format of osteosarcoma and corresponding clinical information data, survival data, and copy number variation data from TCGA database. The samples in the data are all from Homo Sapiens, and the platform is based on Illumina. Additionally, single-cell transcriptome sequencing data of normal transcriptome and osteosarcoma samples were downloaded from GEO database. To ensure data consistency and comparability, we standardized the expression data using the limma package in R, employing a Log2 transformation[8], and the normalized expression data was visualized using box diagrams. There are 88 cases of TARGET osteosarcoma in the TCGA dataset, of which 85 cases with integrated clinical information. After integrating copy number variations, 26 cases were included in this study. Moreover, we also utilized the GSE16088 dataset from GEO database, which comprised 6 normal tissue samples and 14 osteosarcoma samples. Additionally, we accessed the GSE152048 dataset, which comprised single-cell transcriptome sequencing data of osteosarcoma. And 6 primary osteosarcoma samples of GSE152048 were selected for this study. HUGO Gene Nomenclature Committee (HGNC) [9] is responsible for providing a unique, standard, and widely distributed symbol on protein-coding genes of the human genome. The mRNA expression profiles were obtained using HGNC mRNA gene annotation file. Single-cell data processed We processed the single-cell data for each sample used R package Seurat 4.3.0. Cells expressing fewer than 300 genes or with mitochondrial genes comprising more than 10% were filtered out from further analysis[10]. Subsequently, we utilized the R package DoubletFinder 2.0.3 to eliminate potential doublet cells from the sample[11]. Following this filtration, a total of 58,241 cells were retained for subsequent analysis. To process Seurat objects for each individual sample. we employed the Read10× function. After normalization, the gene expression score was multiplied by 10,000 after adding 1 to prevent the logarithm from being 0, and then it was converted to natural logarithm values. Next, 3000 highly variable genes (HVGs) in the dataset were identified by using the "SelectIntegrationFeatures" function. We then scaled the data using "ScaleData" to control for the effects of sequencing depth and mitochondrial genes expression. To merge the samples and mitigate the batch effect, the R package Harmony 2.0 was utilized[12]. Subsequently, principal component analysis (PCA) was applied to identify significant principal components (PCs)[13]. The Elbowplot function was used to visualize the distribution of p-value. Finally, 30 PCs were selected for t-distributed stochastic neighbor embedding (tSNE) analysis. Using the default "FindNeighbors" parameter and 30 PCs dimension parameters, we constructed a k-nearest neighborhood based on Euclidean distances in PCA space. The Louvain algorithm, implemented through the "FindClusters" function, was then applied to optimize cell clustering. Using a resolution of 0.1, the cells were divided into 9 distinct clusters. Finally, we employed the "RunTSNE" function to perform dimensionality reducti, enabling us to visualize and analyze the dataset. Cell type identification Cell types can be identified through the utilization of cell type marker genes. For Osteoblastic osteosarcoma, the marker genes consist of COL1A1, CDH11, and RUNX2. On the other hand, Osteoblastic_proli, which refers to proliferating osteoblastic osteosarcoma, is characterized by the marker genes PCNA and MKI67. Osteoclast cells are distinguished by the expression of CTSK and MMP9. Tumor Infiltrating Lymphocytes (TILs) are identifiable through the presence of IL7R, CD3D, and NKG7. Myeloid cells are marked by CD74, CD14, and FCGR3A. Fibroblasts are identified by the genes COL1A1, LUM, and DCN. Pericytes are characterized by the expression of ACTA2 and RGS5. Mesenchymal stem cells (MSCs) are recognized by the presence of CXCL12, SFRP2, and MME. Endothelial cells are distinguishable by the expression of PECAM1 and VWF[6]. To determine the differentially expressed genes (DEGs) among these cell types, we used FindAllMarkers and Wilcoxon rank-sum test to compare the gene expression profiles across different subpopulations. Differential expression of angiogenesis-related genes between cells We retrieved 36 angiogenesis-related genes from the HALLMARK_ANGIOGENESIS entry of "Hallmark Gene Set" in the MSigDB database[14]. Subsequently, we intersected them with the DEGs between cell types to identify angiogenesis-related genes that exhibited differential expression. To visualize the expression patterns of these genes, we employed the "DoHeatmap" function to generate a heat map. Additionally, we calculated the correlation among these genes using the "COR" function and depicted it as a correlation heat map. Analysis of intercellular communication By leveraging the Python package CellphoneDB[15], we integrated single-cell expression profiles to calculate the number of interactions cell-to-cell communication and identify the interacting receptors and ligands. To visualize the interaction intensity among nine distinct cell types (Osteoblastic, Osteoblastic proli, Osteoclast, TIL, Myeloid cells, Fibroblasts, Pericyte, MSC, Endothelial cells), we generated a heat map. We used R packet CellChat[16] to calculate the communication patterns between cell subpopulations. Using CellChat to identify the significant interaction of ligand-receptor pairs through ligand-receptor interaction probability and perturbation test. The resulting cell-cell communication network was then constructed by integrating the number or strength of significantly interact ligand-receptor pairs across cell types. The interaction strength of osteoblastic, osteoblastic proli, osteoclast, TIL, myeloid cells, fibroblasts, pericyte, MSC and endothelial cells was shown through a circular graph. Subsequently, cellular communication mediated by intercellular ligand receptors was visualized through bubble plots. We also demonstrated several kinds of intercellular ligand-receptor interaction networks using the R-packet iTalk (https://github.com/Coolgenome/iTALK). GSVA analysis Gene Set Variation Analysis (GSVA)[17] is a nonparametric and unsupervised algorithm that transforms gene expression data. We downloaded the C2.cp.kegg.v7.5.1.symbols.gmt dataset from Molecular Signatures Database (MSigDB). Using the "gsva" package in R, we analyzed the single-cell data of osteosarcoma to obtain the GSVA enrichment score for each cell corresponding to each pathway. Subsequently, we employed the R package limma 3.50.0 [8] to identify pathways with significant differences ( p value < 0.05). The scores of pathway activity of cells in each group and other remaining groups were evaluated for difference test analysis. Finally, the Top3 pathways with t values arranged from highest to lowest in each group were plotted. Protein-protein interaction analysis We inputed the differentially expressed TOP100 genes from each subpopulation into the String database[18] to conduct protein interaction network analysis. During this pocess, we preserved the interaction relationships that had been experimentally validated. Following this, degree algorithm in Cytoscape was used to identify key genes among protein interactions. The top20 genes were then selected for visualization. Abundance analysis of immune cell infiltration CIBERSORT[19] is a tool that employs linear support vector regression to deconvolute the expression matrix of human immune cell subtypes. By referencing the LM22 database, CIBERSORT quantify the relative expression abundance of 22 immune cell types, including B cells naive, B cells memory, plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells helper, T cells regulatory (Tregs), T cells gamma delta, NK cells resting, NK cells activated, monocyte, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, eosinophils and neutrophils. We analyzed osteosarcoma datasets from TARGET database to obtain the proportions of different immune cell types using CIBERSORT package in R. WGCNA analysis Weighted Gene Correlation Network Analysis (WGCNA)[20] aims to identify clusters of co-expressed gene, explore the association between gene networks and phenotypes, and investigate the core genes within the network. Firstly, we selected the scores of 22 immune cell type as the trait data for WGCNA analysis. Using the WGCNA package (version 1.71) in R, we calculated the soft threshold using pickSoftTreshold function. The optimal soft threshold was determined to be 3, ensuring the construction of a scale-free network. Subsequently, based on the soft threshold, we constructed according a scale-free network. The networkwas then used to generate a topology matrix and perform hierarchical clustering. By setting the minimum number of genes in the module to 30, we dynamically cutted the modules and calculated the eigengenes. Eigengenes represent the overall expression profile of a module and are crucial for assessing the correlation between modules. After calculating the eigengenes, we constructed a correlation matrix between the modules and performed hierarchical clustering. The modules with a correlation above 0.25 were merged again, and finally 16 modules were obtained. Next, we calculated the correlation between gene modules and phenotypes using pearson to identify the modules that were associated with the phenotypes. The modules with the highest correlation coefficient were intersected with angiogenesis genes to obtain immune-related angiogenesis genes. Differential expression analysis of immune-related angiogenesis genes Obtained the expression matrix of OSTEOSARCOMA from the GSE16088 dataset, the limma package in R was employed to identify the DEGs between the tumor group and the normal group. The DEGs fulfilled the requirements of adj. p value1. To visualize the differential expression patterns of these DEGs, the ggplot2 package was leveraged to generate a volcano plot, in which immune related angiogenesis genes were highlighted. Furthermore, the pheatmap package was utilized to create heat maps to visualize the expression variations of 10 immune-related angiogenesis genes between tumor and normal tissues in OSTEOSARCOMA patients. GO analysis of immune-associated angiogenesis genes Gene Ontology (GO) analysis is a prevalent technique for large-scale functional annotation and enrichment[21], encompassing biological process (BP), molecular function (MF) and cellular component (CC). We used the R package clusterProfiler[22] to perform GO annotation analysis of immune-related angiogenic genes. The screening standard q value was < 0.05, which was deemed statistically significant. Additionally, the Benjamini-Hochberg method (BH) for p-value correction was used. Consistent clustering for classification of osteosarcoma R packet ConsensusClusterPlus[23] was used to cluster the gene expression profile of osteosarcoma in TARGET database, focusing on a panel of 471 immune genes and angiogenesis genes. Spearman method was emploed to calculate the distance between genes, and the Partitioning Around Medoids (PAM) clustering algorithm was chosen for its robustness in handling outlier genes. Through a rigorous analysis of matrix heat map, consistency cumulative distribution function maps, and delta area plots, the osteosarcoma cells subtypes: Cluster1, Cluster2, Cluster3, and Cluster4, were identified. A box plot was generated to compare the expression levels of immune-related angiogenesis genes across these four subtypes, and the t-test was used for statistical significance. Significant differences were considered when P value < 0.05. To gain insights into the immune-related functions of these osteosarcoma subtypes, the immune.gmt dataset was downloaded from MSigDB. And then, the osteosarcoma data were analyzed using the "SSGSEA" method of GSVA package to obtain scores of immune-related functions of different samples. A box plot was constructed to visualize the differences in these immune-related functions among the four subtypes. Statistical significance was again determined using the t-test, with a P value < 0.05 deemed significant. Prognostic modeling We conducted a prognostic analysis utilizing gene expression data and survival information from 85 osteosarcoma samples source from TARGET database. Initially, a univariate Cox regression analysis was conducted on 471 immune genes and angiogenesis genes, enabling us to preliminarly identify those genes that exhibited significant associations with overall survival (P value < 0.05). Subsequently, the osteosarcoma samples were stratified into two distinct sets: comprising 53 cases, used to establish a prognostic model, and a validation set consisting of 32 cases. Employing the prognostic genes selected through the univariate Cox analysis, we utilized lasso-Cox regression analysis to construct the prognostic model. The risk score was determined using the following calculation formula: Coef (GeneI) represents lasso-Cox regression coefficient; Expression (Gene i ) represents the expression value of each gene, and n represents the number of genes[24]. 85 osteosarcoma samples were categorized into high-risk and low-risk groups based on the median risk score derived from the training set. To assess the prognostic value and predictive accuracy of the model, we employed Kaplan-Meier survival curve analysis to analyze overall survival and generated time-dependent receiver operating characteristic curves (ROC). Construction and correlation analysis of Nomogram model After removing clinical features with null from the TARGET osteosarcoma database using R, we obtained three characteristics: age, gender and M stage. To further investigate factors related to patient prognosis, we conducted both univariate and multivariate Cox risk regression analysis, incorporating the prognostic risk score alongside these clinical characteristics. For visualization, we employed the R package 'forestplot'. Leveraging the nomogram function within R package RMS, we constructed a nomogram model specific to the prognostic factors significantly associated with outcomes in the TARGET osteosarcoma database. For enhanced visualization, we utilized the 'ggplot' package. Time-dependent ROC analysis was performed on the predicted scores from nomogram model, considering the overall survival status and survival time of patients at one, three and five years, respectively, using R package timeROC. Multivariate ROC analysis of nomogram model was performed incorporating the predicted score, overall survival status and survival time of patients. Finally, we evaluated the accuracy and resolution of the nomogram using a calibration curve. We employed the R package RMS to devise the nomogram and calibration curve. We leveraged R package ggDCA to evaluate the one-year, three-year, and five-year survival rates of patients using the nomogram model. Immune infiltration analysis The abundance of 22 immune cell types was calculated based on the immune invasion matrix of osteosarcoma samples from TARGET. We analyzed the association between these immune cells and prognostic genes. In addition, t-test was used to compare the abundances of the 22 immune cell types between the high-risk and low-risk groups, with p value < 0.05 considered as statistically significant. The correlation between the expression of prognostic genes and the infiltration of immune cell was also calculated. Copy number analysis Copy number data of osteosarcoma were downloaded from TARGET database. A copy number of 0 was defined as double deletion. A copy number of 1 was designated as single deletion. A copy number of 2 was considered normal. A copy number of 3 was defined as single gain. And a copy number of 4 or more was classified as amplification. Then a boxplot was drawn to compare the differences in copy numbers and gene expression. Statistical analysis With the exception of CellphoneDB, which employed Python, all remaining data calculations and statistical analysis were carried out utilizing R language. When comparing continuous variables across two groups, the statistical significance of normally distributed variables was determined using the independent Student’s t test. For variables that did not follow a normally distributed, differences were assessed through the Mann-Whitney U test. To evaluate the statistical significance of categorical variables between the two groups, we employed either the Chi-square test or Fisher's exact test, depending on the circumstances. Furthermore, the correlation coefficients among different genes were derived through Pearson correlation analysis. Result 1. Cell heterogeneity in osteosarcoma We conducted a single-cell RNA sequencing (scNA-seq) analysis on tumor samples from 6 osteosarcoma patients within the GSE152048 dataset to delve into the cellular composition. After rigorous quality control, we filtered out cells with mitochondrial gene content >10%, those with feature counts 6000, and eliminated duplicates, resulting in a final dataset of 58,241 cells for further analysis. Nine cell types were identified by artificial annotation based on gene expression and marker genes (Figure 2A). The marker genes utilized for this classification were detailed in table S1. Cell subpopulations revealed the following distribution: osteoblastic cells (18014,30. 93%), osteoblastic_proli (3172, 5.44%), osteoclasts (6194, 10.64%), TILs (318, 5.53%), myeloid cells (15309, 26.29%), fibroblasts (6390, 10.97%), pericytes (1844, 3.17%), MSCs (115, 1.91%), endothelial cells (2984, 5.12%). The differential expression of 23 marker genes that distinguish these subpopulations was visualized in figure 2B (violin diagram) and figure 2C (bubble diagram), demonstrating the preferential expression of these marker genes in respective cell subpopulations. We calculated the DEGs between the cell subpopulations using the "FindAllMarkers" function. Figure 2D showed the top three DEGs for each of the nine main cell subpopulations. Additionally, we tabulated the number and proportion of each cells subpopulation across the 6 samples. Notably, osteoblastic cells comprised a higher proportion in BC5 and BC6, myeloid cells were more abundant in BC2 and BC16, and BC22 exhibite a higher content of fibroblast cells (Figure 2E and Figure 2F). 2. Differentially expressed angiogenesis-related genes among cell subsets We identified 19 differentially expressed angiogenes-associated genes (AAGs) by intersecting DEGs among cell subsets with angiogenesis-related genes. A heat map was used to visualize the expression of 19 angiogenesis-related genes across various cell subsets (Figure 3A). The figure illustrated that angiogenesis-related genes were significantly upregulated in osteoblastic osteosarcoma and proliferative osteoblastic OSTEOSARCOMA, whereas their expression was relatively low in myeloid cells. Furthermore, MSCs also exhibited high expression levels of angiogenesis genes, including COL3A1, COL5A2 and FSTL1. Additionally, certain AAGs, such as KCNJ8, NRP1, POSTN, TIMP1 and VCAN were expressed in pericytes. Notably, S100A4 and SPP1 were highly expressed in osteoclast. Our correlation analysis of these 19 angiogenesis-related genes revealed that THBD, NRP1, SLCO2A1, JAG2 and STC1 were positively correlated, whereas VCAN, COL5A2, FSTL1, FGFR1, TIMP1, POSTN, COL3A1 and LUM also demonstrated positive correlations (Figure 3B). These genes were further visualized in the tSNE dimension reduction diagram for osteoblastic and osteoblastic1 proli cells (Figure 3C). 3. Cell communication and hub genes We inferred and quantified the communication among 9 cell subtypes using CellphoneDB and CellChat respectively, and then visualized cell communication intensities through a circle diagram and a heat map (Figure 4 a-b). A cursory examination revealed that osteoblastic cells, fibroblasts and MSCs exhibited a relatively higher cumulative signal intensity, indicating a relatively higher cumulative signal intensity, indicating a greater level of communication activity among these cell types. Similarly, fibroblasts, osteoblastic proli cells, osteoclasts and MSCs also demonstrated a significant cumulative signal intensity. Subsequently, we calculated all important ligand-receptor pairs that mediated from TIL to various cell types, including fibroblasts, osteoblastic proli, myeloid cells, osteoclasts and MSCs (Figure 4C). Notably, CD47-related pathway emerged as a pivotal communication channel between TIL and myeloid cells. The communication intensity between TIL and osteosarcoma cells appeared to be relatively low, with CD44 standing out as the most important ligand-receptor pair. ITGB1 and SPP1 related ligand-receptor pathways played an important role in the communication between fibroblasts, osteoblastic proli, osteoclasts and MSCs (Figure 4D). Subsequently, GSVA analysis was performed on single-cell data of osteosarcoma to derive the GSVA enrichment score for each cell across various pathways. Then, the pathway activity among diverse cell subsets was evaluated to identify those exhibiting significant differences (P value < 0.05). Furthermore, the differential analysis was performed to ascertain the pathways that demonstrated notable variations in the relative pathway activity across different cell subsets. A heat map was generated to visualize the top 27 pathways from each group, arranged in descending order of each group (Figure 4E). Additionally, protein interaction analysis was also conducted focusing on the top100 DEGs across various subpopulations, and the top 20 hub genes were obtained: FN1, PTPRC, VEGFA, IL1B, COL1A1, MMP9, CD8A, JUN, ITGB1, PECAM1, CXCL8, COL1A2, MMP2, CCL2, COL3A1, VWF, SPP1, FCGR3A, ITGB2, and TYROBP (Figure 4F). 4. WGCNA analysis We performed WGCNA analysis of osteosarcoma data obtained from the TARGET database. We used immune-related genes and 36 angiogenesis-related genes as the gene set for WGCNA analysis. Immune-related genes were obtained from Immport database. Initially, we determined that a soft threshold β was 3 (R-square= 0.94) would optimize the consistency of the network with the characteristics of scale-free network (Figure 5A). Next, we gerneated a hierarchical clustering tree using the correlation coefficient between genes, leading to the identification of 16 similar gene modules (Figure 5B). A heat map was constructed to visualize the correlation of modules. We employed CIBERSORT to quantify the abundance of 22 immune cell types in osteosarcoma, considering them as phenotypic traits. Among the 10 modules, the MEred module exhibited the highest correlation with CD8+ T cells, while the MEturquoise module showed the strongest association with resting dendritic cells (Figure 5C). These two modules collectively encompassed 471 genes, including 10 angiogenesis-related genes (FGFR1, LPL, LRPAP1, KCNJ8, COL3A1, MSX1, JAG2, VAV2, PGLYRP1, APOH), all of which had immune-related angiogenesis functions. 5. Differential expression of hub genes in angiogenesis The GSE16088 dataset was used to conduct a differential expression analysis aomparing the tumor group with the normal group. We employed R package limma to analyze the transcriptome data of OSTEOSARCOMA patients, utilizing |log2FC| > 1 and adj.p value < 0.05 as threshold criteria for screening. We identified 2270 upregulated genes and 1981 downregulated gene. Volcano and heat maps depicting these DEGs were presented in Fig 6A and Fig 6B, respectively, with 10 immune-related angiogenesis genes highlighted in yellow in Fig 6A. The finding revealed that FGFR1, LPL, COL3A1 and MSX1 were highly expressed in tumor tissues, whereas APOH was low expressed in tumor tissues. The detailed expressions of these 10 angiogenesis-related genes were shown in Table S4. Subsequently, we conducted a GO enrichment analysis on these 10 genes, revealed that their functions were primarily concentrated in BP, such as positive regulation of lipase activity, regulation of plasma lipoprotein particle levels, blood coagulation, skeletal system morphogenesis, hemostasis. At the CC level, they are associated with chylomicron, very high-density lipoprotein particle, triglyceride plasma lipoprotein particle, plasma lipoprotein particle, lipoprotein particle, protein-lipid complex. In terms of MF, these genes were involved in glycosaminoglycan binding, heparin binding, sulfur compound binding, growth factor binding, 1-acyl-2-Lysophosphatidylserine acylhydrolase activity, lipoprotein lipase activity. 6. Classification of osteosarcoma The TARGET osteosarcoma samples were categorized through consistent clustering. Examination of the consistent clustering heat map (Figure 7A), consistent cumulative distribution function (CDF) curve (Figure 7B) and delta area curve (Figure 7C) revealed that the optimal clustering number was 4. Consequently, osteosarcoma was stratified into four subtypes: cluster1, cluster2, cluster3 and cluster4. A comparative analysis of the expression of 10 immune-related angiogenesis genes across these subtypes (Figure 7D) revealed significant differential expression among them (p value < 0.05). Notably, APOH exhibited differential expression not only between Cluster1 and Cluster2 but also between Cluster1 and Cluster3. Futhermore, we assessed the abundance of 22 immune cell types across the four osteosarcoma subtypes. The cell types that exhibited notable abundance differences among the four subtypes activated CD4 T cells, activated CD8 T cells, gamma delta T cells, MDSCs, macrophage, mast cells, monocytes and regulatory T cells. Among them, cluster3 osteosarcoma, which was characterized by higher abundance of activated CD8 T cell, gamma delta T cells, MDSCs and regulatory T cells, was identified as immune infiltrating osteosarcoma. 7. Prognostic modeling Utilizing gene expression data from TARGET database for osteosarcoma, 59 genes associated with overall survival were identified by univariate Cox regression analysis from 471 genes related to immunity and angiogenesis within MEred and MEturquoise modules. The 85 samples were stratified into a training set (consisting of 53 cases) and a validation set (consisiting of 32 cases). Employing Lasso-Cox regression analysis on the training set samples, we calculated the Cox regression coefficients for the 59 genes (Figure 8a-b). Notably, eight genes - CD48, ESRRA, GNAI1, GNRH1, JAG2, KRAS, SECTM1, and TRPC4AP - emerged as significant predictors with non-zero coefficients. The risk score for each sample was calculated according to the following formula: risk score = -2.034× Exp(CD48) + (2.969) × Exp(ESRRA) + 1.177× Exp(GNAI1) + (2.492) × Exp(GNRH1) + (0.907) × Exp(JAG2) + (-1.600) × Exp(KRAS) + (-0.779) × Exp(SECTM1) + (-3.855) × Exp(TRPC4AP). Based on the median risk score derived from the training set, the osteosarcoma samples in TARGET database were categorized into the high-risk and low-risk groups (training set: 26 high-risk, 27 low-risk; validation set: 20 high-risk and 12 low-risk). As evident from the distribution of risk scores and corresponding survival status for both the training set (Figure 8C) and validation set (Figure 8D), an increase in risk score correlated with a higher risk of mortality and shorter survival time. The prognostic model's risk score demonstrated excellent predictive performance, with 1-year, 3-year, and 5-year area under the curve (AUC)values exceeding 0.65 for both the training set (Figure 8E) and validation set (Figure 8F). Survival analysis confirmed a significant difference in overall survival between the high-risk and low-risk groups in both the training sets (Figure 8G) and validation sets (Figure 8H) (P value < 0.05). 8. Verify the prognostic independence of the model and the prognostic significance of other risk factors The results of the univariate Cox regression analysis from the TARGET database indicated that age (HR = 1.009, 95% CI = 0.894-1.139, P = 0.880), risk score (HR = 1.004, 95% CI = 1.000-1.008, P = 0.068), M stage (HR = 1.719, 95% CI = 0.444-6.654, P = 0.433), and gender (HR = 1.120, 95% CI = 0.316 - 3.973, P = 0.880) did not show significant associations with prognosis. However, it is noteworthy that the P value for risk score was close to significance (P = 0.068). (Figure 9A) (Table S5). A multivariate Cox regression analysis was conducted, employing these factors as covariates to further assess their prognostic significance. The analysis revealed that age (HR = 1.020, 95% CI = 0.882-1.180, P = 0.788), gender (HR = 1.427, 95% CI = 0.301-6.780, P = 0.654), and M stage (HR = 2.186, 95% CI = 0.502-9.517, P = 0.297) were not independent prognostic factors for osteosarcoma. However, the risk score emerged as a significant independent prognostic factor for overall survival among osteosarcoma patients (HR = 1.005, 95% CI = 1.000 - 1.010, P < 0.05) (Figure 9B) (Table S5). Subsequently, we further constructed a nomogram model to predict 1-year, 3-year and 5-year overall survival, which incorporated four factor: age, gender, M stage and prognostic risk score (Figure 9C). By leveraging the R package timeROC, we evaluated the predictive performance of the nomogram model by calculating AUC values for 1-year, 3-year, and 5-year survival predictions. The all AUC values exceeded 0.9, and outperformed the prognostic risk score alone. In addition, we conducted prognostic calibration analysis for 1-year, 2-year, and 3-year survival probabilities using variables from both univariate and multivariate Cox regression modles. And calibration curves were plotted (Figure 9D). 9. Immunocorrelation analysis The results of immunocorrelation analysis of the TARGET osteosarcoma dataset indicated a rich immune cells infiltration within the samples. We calculated the correlation between immune scores and immune cells infiltration, discovering that the infiltration of macrophages M1, T cells follicular helper, T cells CD8, and Tregs exhibited correlation. A correlation was observed between the infiltration of B cells naive and plasma cells was correlated. T cells CD4 memory resting and NK cells resting was correlated. We also analyzed the differences in the abundance of immune cells between the highrisk and low-risk groups. As shown in Figure 10A and Figure 10C, there were significant differences in the abundance of T cells CD8, plasma cells, T cells CD4 naive, T cells Helper, T cells regulatory (Tregs), monocyte, macrophages M1, and dendritic cells resting between the two group (P value < 0.05). The abundance of T cells CD8 was lower in the high-risk group. Next, we investigated the relationship between prognostic gene expression and T cells CD8 infiltration. The expression of CD48 and SECTM1 was positively correlated with T cells CD8 infiltration (Figure 10D and Figure 10G), while the expression of GNAI1 and JAG2 was negatively correlated with T cells CD8 infiltration (Figure 10E and Figure 10F). 10. Differential genes and CNV analysis in highrisk and low-risk groups The model was employed to assess the prognostic risk of each sample in the TARGET osteosarcoma database. The cutoff value was set as the median prognostic risk score across all samples. Using this cutoff value, the samples were stratified into high-risk and low-risk groups, and the DEGs between these two risk groups was calculated employing the R package limma. A total of 113 genes were overexpressed in the low-risk patients and 36 genes were overexpressed in the high-risk patients. We used a heat map to show the expression of TOP30 differential expressed genes (Figure 11A). Subsequently, KEGG pathway enrichment analysis was performed on the TOP30 genes. The results indicated that these DEGs were enriched in various BPs, including complement and coagulation cascades, staphylococcus aureus infection, pertussis, rheumatoid arthritis, asthma (Figure 11C). We also correlated the expression of TOP30 genes with the results of the KEGG enrichment analysis, preenting them in a string diagram (Figure 11B) and a network diagram (Figure 11D) respectively. Finally, we explored the correlation between TOP30 gene expression and copy number variation. Our analysis revealed that the expressions of CGREF1, NGEF, PDGFD, RHBDL2, TCN2 and TRAC were associated with copy number amplification, and most of the genes with significant associations exhibited copy number amplification of 4 or higher, namely, indicating a strong link between gene expression and amplification. Discussion Osteosarcoma, a notoriously aggressive bone tumor, has long posed significant therapeutic challenges[ 4 , 25 , 26 ]. However, the remarkable advancements in single-cell sequencing technology have revolutionized our understanding of the cellular composition, gene expression patterns, and intercellular dynamics of osteosarcoma. This intricate level of analysis has provided profound insights into its pathogenesis and offers hope for the discovery of novel therapeutic strategies. By analyzing the cellular heterogeneity, we have been able to delineate diverse 9 cell subpopulations from osteosarcoma samples within the GSE152048 dataset. These subpopulations, which include osteoblastic cells, osteoblastic_proli, osteoclasts, TILs, myeloid cells, fibroblasts, pericytes, MSCs and endothelial cells. Osteoblastic cells comprised a higher proportion in BC5 and BC6, myeloid cells were more abundant in BC2 and BC16, and BC22 exhibite a higher content of fibroblast cells. This heterogeneity not only reflects the diverse nature of the cell types involved but also indicates varying levels of gene expression across each subpopulation. This cellular diversity was presumably a key factor underlying osteosarcoma's resilience to traditional therapeutic approaches. Then, top three DEGs for each of the nine main cell subpopulations were analyzed. Consequently, the development of targeted therapeutic strategies, specifically tailored to address these distinct cell subpopulations, could potentially lead to more favorable treatment outcomes. Our meticulous examination of gene expression patterns has been particularly crucial in elucidating the role of angiogenesis, a crucial process in tumor growth and metastasis. We identified 19 differentially expressed AAGs by intersecting DEGs among cell subsets with angiogenesis-related genes. AAGs were highly expressed in osteoblastic osteosarcoma and proliferative osteoblastic osteosarcoma, whereas their expression was relatively low in myeloid cells. Osteosarcoma, as a highly vascularized malignancy, heavily relies on angiogenesis for its aggressive nature [ 27 ]. Our findings have revealed that genes related to angiogenesis were abundantly expressed in osteoblasts and their proliferative counterparts. MSCs exhibited high expression levels of FSTL1. Blocking FSTL1 could treat osteosarcoma effectively by inhibiting cell proliferation, invasion, and various other crucial cellular functions[ 28 ]. NRP1, POSTN, TIMP1 and VCAN were expressed in pericytes. MicroRNA-1247 inhibited the viability and metastasis of osteosarcoma cells via targeting NRP1 and mediating Wnt/β-catenin pathway[ 29 ]. Our correlation analysis of these 19 angiogenesis-related genes revealed that THBD, NRP1, SLCO2A1, JAG2 and STC1 were positively correlated, whereas VCAN, COL5A2, FSTL1, FGFR1, TIMP1, POSTN, COL3A1 and LUM also demonstrated positive correlations. This discovery aligns with prior research and suggests that suppressing the expression of these angiogenesis-associated genes could potentially disrupt the angiogenesis process in osteosarcoma, thereby inhibiting tumor growth and metastasis. Furthermore, our study employed bioinformatics techniques to delve deeper into the intricate intercellular interactions within osteosarcoma. Osteoblastic cells, fibroblasts and MSCs exhibited a relatively higher cumulative signal intensity, indicating a relatively higher cumulative signal intensity, indicating a greater level of communication activity among these cell types. Similarly, fibroblasts, osteoblastic proli cells, osteoclasts and MSCs also demonstrated a significant cumulative signal intensity. Intercellular communication plays a pivotal role in the tumor microenvironment, serving as a crucial conduit for signaling among diverse cell types and orchestrating tumor growth, invasion, and metastasis. Through a meticulous analysis of the communication network and the identification of key ligand-receptor pairs, we gained profound insights into the complex interplay between osteosarcoma cells and their surrounding microenvironment. This understanding could pave the way for the development of innovative therapeutic strategies that specifically target these interactions, thus halting tumor progression. GSVA analysis was performed to derive the GSVA enrichment score for each cell across various pathways. Additionally, protein interaction analysis was also performed on the top100 DEGs across various subgroups. WGCNA analysis of osteosarcoma data was performed to obtain the immune-related genes. The MEred module exhibited the highest correlation with CD8 + T cells, while the MEturquoise module showed the strongest association with resting dendritic cells. Research has demonstrated that IL-35 could lead to dysfunction of CD8 + T cells and thus constrain the anti-tumor immune response in osteosarcoma[ 30 ]. This study indicated that FGFR1, LPL, COL3A1 and MSX1 were highly expressed in osteosarcoma tumor tissues, whereas APOH was low expressed in osteosarcoma tumor tissues. The activation of nuclear FGFR1 induced histone H3 phosphorylation at Ser 10 and c-jun/c-fos expression to contribute osteosarcoma cell survival rendering radiation resistance[ 31 ]. Research has unequivocally confirmed that the miR-29 family exerts a tumor suppressive function in modulating MTX resistance and osteosarcoma cell apoptosis, accomplished through the regulation of COL3A1[ 32 ]. In this study, osteosarcoma was stratified into four subtypes: cluster1, cluster2, cluster3 and cluster4, and the abundance of 22 immune cells types across the four osteosarcoma subtypes was evaluated. Cluster3 osteosarcoma, which was characterized by higher abundance of activated CD8 T cell, gamma delta T cells, MDSCs and regulatory T cells, was identified as immune infiltrating osteosarcoma. Utilizing gene expression data from TARGET database for osteosarcoma, eight genes - CD48, ESRRA, GNAI1, GNRH1, JAG2, KRAS, SECTM1, and TRPC4AP - emerged as significant predictors with non-zero coefficients. The expression of high levels of JAG2 in osteosarcoma patients could potentially be linked to a more favorable response to immune checkpoint blockade therapy [ 33 ]. Hsa-miR-557 effectively inhibited osteosarcoma growth both in vivo and in vitro by modulating the expression of KRAS[ 34 ]. The univariate Cox regression analysis did not detect any significant associations between age, risk score, M stage, and gender with prognosis. However, upon further investigation using multivariate Cox regression analysis, it became evident that age, gender, and M stage were not independent prognostic indicators for overall survival. Notably, the risk score emerged as a significant and independent prognostic indicator for overall survival among osteosarcoma patients. This finding suggests that the risk score, likely encompassing multiple clinicopathological variables, may hold key information for predicting patient outcomes and guiding treatment strategies in osteosarcoma. In our analysis, we calculated the correlation between immune scores and the infiltration of immune cells. Our findings revealed a significant correlation between the infiltration of macrophages M1, follicular helper T cells, CD8 + T cells, and regulatory T cells. This suggests that these immune cell subsets play a pivotal role in modulating the immune response in the context of our study. Further investigation is needed to understand the functional significance of these correlations and their potential implications in disease pathogenesis and therapeutic responses. The results of our analysis indicate that the expressions of CGREF1, NGEF, PDGFD, RHBDL2, TCN2, and TRAC are closely associated with copy number amplification. Furthermore, it is noteworthy that the majority of genes displaying significant associations exhibited copy number amplification of 4 or higher. This observation strongly suggested a robust correlation between gene expression levels and copy number amplification in these specific genes. Such a link could have significant implications for understanding the molecular mechanisms underlying osteosarcoma and other related malignancies, as well as for developing potential therapeutic targets. Future studies should aim to further explore the functional roles of these genes and their association with disease progression and prognosis. Indeed, there are inherent limitations in this study that need to be acknowledged. Firstly, the absence of validation using a larger sample size is a significant limitation. This could have affected the reliability and generalizability of our findings, as a larger sample would provide a more robust basis for validating the observed associations. Secondly, the lack of further evidence from basic experiments is another shortcoming. Basic experiments would have provided deeper insights into the biological mechanisms underlying the observed gene expression patterns and copy number amplifications. Such evidence is crucial for validating and complementing the findings obtained from bioinformatic analysis. In future studies, it is recommended to address these limitations by including a larger sample size for validation and conducting basic experiments to further explore the functional roles of the identified genes. In summary, the utilization of single-cell sequencing technology in osteosarcoma research has significantly enhanced our comprehension of the disease's cellular landscape, gene expression patterns, and intercellular dynamics. These revelations provide precious leads for devising more effective and tailored therapeutic approaches against this aggressive bone tumor. Future endeavors in this field promise to further improve treatment outcomes for osteosarcoma patients. Declarations Ethics approval and consent to participate: Not applicable. Consent for publication ： Not applicable. Availability of data and materials: The data that support the findings of this study were obtained from TCGA and GEO. Derived data supporting the findings of this study are available from the corresponding author on reasonable request. Competing interests: The authors declare that we have no commercial or financial relationships that could potentially be construed as a conflict of interest. Funding: This study was supported by National Natural Science Foundation(82273436). Authors' contributions: HL analyzed the data and wrote the paper. TM and ZQY provided the help of the R language. ML, FG and XJL designed the project. ML selected the analyzed results. All authors read and approved the final manuscript. Acknowledgments: Not applicable. References Chen C, Xie L, Ren T, Huang Y, Xu J, Guo W: Immunotherapy for osteosarcoma: Fundamental mechanism, rationale, and recent breakthroughs . CANCER LETT 2021, 500 :1-10. Mutsaers AJ, Walkley CR: Cells of origin in osteosarcoma: mesenchymal stem cells or osteoblast committed cells? BONE 2014, 62 :56-63. Ciernik IF, Niemierko A, Harmon DC, Kobayashi W, Chen YL, Yock TI, Ebb DH, Choy E, Raskin KA, Liebsch N et al : Proton-based radiotherapy for unresectable or incompletely resected osteosarcoma . CANCER-AM CANCER SOC 2011, 117 (19):4522-4530. Wedekind MF, Wagner LM, Cripe TP: Immunotherapy for osteosarcoma: Where do we go from here? PEDIATR BLOOD CANCER 2018, 65 (9):e27227. Lee HW, Chung W, Lee HO, Jeong DE, Jo A, Lim JE, Hong JH, Nam DH, Jeong BC, Park SH et al : Single-cell RNA sequencing reveals the tumor microenvironment and facilitates strategic choices to circumvent treatment failure in a chemorefractory bladder cancer patient . GENOME MED 2020, 12 (1):47. Zhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, Wang Y, Zhang Z, Yuan T, Ding X et al : Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma . NAT COMMUN 2020, 11 (1):6322. Wang Z, Jensen MA, Zenklusen JC: A Practical Guide to The Cancer Genome Atlas (TCGA) . Methods Mol Biol 2016, 1418 :111-141. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies . NUCLEIC ACIDS RES 2015, 43 (7):e47. Bruford EA, Antonescu CR, Carroll AJ, Chinnaiyan A, Cree IA, Cross NCP, Dalgleish R, Gale RP, Harrison CJ, Hastings RJ et al : HUGO Gene Nomenclature Committee (HGNC) recommendations for the designation of gene fusions . LEUKEMIA 2021, 35 (11):3040-3043. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R: Integrating single-cell transcriptomic data across different conditions, technologies, and species . NAT BIOTECHNOL 2018, 36 (5):411-420. McGinnis CS, Murrow LM, Gartner ZJ: DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors . Cell Syst 2019, 8 (4):329-337 e324. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, Brenner M, Loh PR, Raychaudhuri S: Fast, sensitive and accurate integration of single-cell data with Harmony . NAT METHODS 2019, 16 (12):1289-1296. Kim S, Kang D, Huo Z, Park Y, Tseng GC: Meta-analytic principal component analysis in integrative omics application . BIOINFORMATICS 2018, 34 (8):1321-1328. Qing X, Xu W, Liu S, Chen Z, Ye C, Zhang Y: Molecular Characteristics, Clinical Significance, and Cancer Immune Interactions of Angiogenesis-Associated Genes in Gastric Cancer . FRONT IMMUNOL 2022, 13 :843077. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R: CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes . NAT PROTOC 2020, 15 (4):1484-1506. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q: Inference and analysis of cell-cell communication using CellChat . NAT COMMUN 2021, 12 (1):1088. Hanzelmann S, Castelo R, Guinney J: GSVA: gene set variation analysis for microarray and RNA-seq data . BMC BIOINFORMATICS 2013, 14 :7. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P et al : STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets . Nucleic Acids Res 2019, 47 (D1):D607-D613. Steen CB, Liu CL, Alizadeh AA, Newman AM: Profiling Cell Type Abundance and Expression in Bulk Tissues with CIBERSORTx . Methods Mol Biol 2020, 2117 :135-157. Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis . BMC BIOINFORMATICS 2008, 9 :559. Yu G: Gene Ontology Semantic Similarity Analysis Using GOSemSim . Methods Mol Biol 2020, 2117 :207-215. Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters . OMICS 2012, 16 (5):284-287. Wilkerson MD, Hayes DN: ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking . BIOINFORMATICS 2010, 26 (12):1572-1573. Zhang C, He H, Hu X, Liu A, Huang D, Xu Y, Chen L, Xu D: Development and validation of a metastasis-associated prognostic signature based on single-cell RNA-seq in clear cell renal cell carcinoma . Aging (Albany NY) 2019, 11 (22):10183-10202. Siegel RL, Miller KD, Wagle NS, Jemal A: Cancer statistics, 2023 . CA Cancer J Clin 2023, 73 (1):17-48. Thanindratarn P, Dean DC, Nelson SD, Hornicek FJ, Duan Z: Advances in immune checkpoint inhibitors for bone sarcoma therapy . J BONE ONCOL 2019, 15 :100221. Li YS, Liu Q, Tian J, He HB, Luo W: Angiogenesis Process in Osteosarcoma: An Updated Perspective of Pathophysiology and Therapeutics . AM J MED SCI 2019, 357 (4):280-288. Ogiwara Y, Nakagawa M, Nakatani F, Uemura Y, Zhang R, Kudo-Saito C: Blocking FSTL1 boosts NK immunity in treatment of osteosarcoma . CANCER LETT 2022, 537 :215690. Wei QF, Yao JS, Yang YT: MicroRNA-1247 inhibits the viability and metastasis of osteosarcoma cells via targeting NRP1 and mediating Wnt/beta-catenin pathway . Eur Rev Med Pharmacol Sci 2019, 23 (17):7266-7274. Liu MX, Liu QY, Liu Y, Cheng ZM, Liu L, Zhang L, Sun DH: Interleukin-35 suppresses antitumor activity of circulating CD8(+) T cells in osteosarcoma patients . CONNECT TISSUE RES 2019, 60 (4):367-375. Kim JA, Berlow NE, Lathara M, Bharathy N, Martin LR, Purohit R, Cleary MM, Liu Q, Michalek JE, Srinivasa G et al : Sensitization of osteosarcoma to irradiation by targeting nuclear FGFR1 . Biochem Biophys Res Commun 2022, 621 :101-108. Xu W, Li Z, Zhu X, Xu R, Xu Y: miR-29 Family Inhibits Resistance to Methotrexate and Promotes Cell Apoptosis by Targeting COL3A1 and MCL1 in Osteosarcoma . Med Sci Monit 2018, 24 :8812-8821. Yang L, Long Y, Xiao S: Osteosarcoma-Associated Immune Genes as Potential Immunotherapy and Prognosis Biomarkers . BIOCHEM GENET 2024, 62 (2):798-813. Qiao Z, Li J, Kou H, Chen X, Bao D, Shang G, Chen S, Ji Y, Cheng T, Wang Y et al : Hsa-miR-557 Inhibits Osteosarcoma Growth Through Targeting KRAS . Front Genet 2021, 12 :789823. Supplementary Files Table1.docx Table S1 Table2.docx Table S2 Table3.docx Table S3 Table4.docx Table S4 Table5.docx Table S5 Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-5305987","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":373929269,"identity":"b429c620-0c62-412e-86cb-1d028e7e9e67","order_by":0,"name":"Hao Li","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA4klEQVRIiWNgGAWjYBACPmYeIFnBwGwApAzAQgcIaGEDazljANTCTKwWBqAWxjYDBpAWBuK0sPMe/PBx3h92cwb+A0U32xjk+G4kMH4uwOswvmTJmdsMmC0bmBmMc9sYjCVvJDBLz8DvFwNpXqAWgwMQLYkbbiRAPIhHi/Fv3jkILfXEaDGT5m1AaEkwIEaL5YxjxswGh5kNjHPOSRjOPPOwWRqfFn7+M8Y3PtTIJRscb3xmnFNmI893PPngZ3xaYCAZGC1swKiUALIZG4jQwMBgB8TMD4hSOgpGwSgYBSMOAAD8aDxYeZ6XFAAAAABJRU5ErkJggg==","orcid":"https://orcid.org/0000-0002-1135-1944","institution":"Huazhong University of Science and Technology Tongji Medical College Tongji Hospital","correspondingAuthor":true,"prefix":"","firstName":"Hao","middleName":"","lastName":"Li","suffix":""},{"id":373929270,"identity":"5b76d507-64ea-476a-8992-081d682574b2","order_by":1,"name":"Tian Ma","email":"","orcid":"","institution":"Huazhong University of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Tian","middleName":"","lastName":"Ma","suffix":""},{"id":373929271,"identity":"a9d9248e-102d-4e89-8373-2a86da134a79","order_by":2,"name":"Zhiqian Yi","email":"","orcid":"","institution":"Huazhong University of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Zhiqian","middleName":"","lastName":"Yi","suffix":""},{"id":373929272,"identity":"a6bda71b-f803-4857-8435-57c081f0ff19","order_by":3,"name":"Fang Gao","email":"","orcid":"","institution":"Huazhong University of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Fang","middleName":"","lastName":"Gao","suffix":""},{"id":373929273,"identity":"7fbcf492-bb2c-47c1-a6d6-809b165abf4a","order_by":4,"name":"Xiaojuan Li","email":"","orcid":"","institution":"Huazhong University of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Xiaojuan","middleName":"","lastName":"Li","suffix":""},{"id":373929274,"identity":"b818867a-84c5-4636-9ea6-72a35dd99475","order_by":5,"name":"Mi Li","email":"","orcid":"","institution":"Huazhong University of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Mi","middleName":"","lastName":"Li","suffix":""}],"badges":[],"createdAt":"2024-10-21 16:22:45","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-5305987/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-5305987/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":70922354,"identity":"6bd993b9-74fc-4db6-a9d9-617387bf62e8","added_by":"auto","created_at":"2024-12-09 08:48:36","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":121329,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/e486847ad95d641a34320b6f.jpg"},{"id":70923922,"identity":"0a512cd8-ebcd-4652-b342-bd1ba911241a","added_by":"auto","created_at":"2024-12-09 08:56:36","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":570048,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/17a856126f34b290a2f76088.jpg"},{"id":70922357,"identity":"bcb1f2c6-be2b-425b-9674-c44c8f92c0af","added_by":"auto","created_at":"2024-12-09 08:48:36","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":459853,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/878683a19bb78e4af3a31382.jpg"},{"id":70922364,"identity":"6e36ae50-d316-4f1c-8eaa-9e3a6dce6b15","added_by":"auto","created_at":"2024-12-09 08:48:36","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":248605,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/dedf5511d6365d33ce5cda28.jpg"},{"id":70921584,"identity":"7e04a779-09c6-45cb-b5f9-f9a049a9f0cf","added_by":"auto","created_at":"2024-12-09 08:40:37","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":579172,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/5b55e4a76366625fac0bc83c.jpg"},{"id":70921579,"identity":"717aad8c-610c-4c04-a66b-90ecbd175cb8","added_by":"auto","created_at":"2024-12-09 08:40:36","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":356610,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/523a8601e1b7fa2a6f97d44a.jpg"},{"id":70921582,"identity":"7aa92185-e449-4b02-a6fe-dd757f9525c1","added_by":"auto","created_at":"2024-12-09 08:40:36","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":283105,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure7.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/5390b26ee2f3910000694924.jpg"},{"id":70922361,"identity":"8eff810a-e565-4119-a6f4-cdd600007e05","added_by":"auto","created_at":"2024-12-09 08:48:36","extension":"jpg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":202558,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure8.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/dcde79167905141bce158684.jpg"},{"id":70923923,"identity":"b2bc0c55-92d5-4859-b33a-3bf72127d468","added_by":"auto","created_at":"2024-12-09 08:56:36","extension":"jpg","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":108045,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure9.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/42c965661939f09d2fd62c6b.jpg"},{"id":70921576,"identity":"0b0b4f94-c2d7-42f3-a312-3ff96a859d0e","added_by":"auto","created_at":"2024-12-09 08:40:36","extension":"jpg","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":485293,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure10.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/136dce8355e0e89613cf318b.jpg"},{"id":70922366,"identity":"dc31e6ea-4bd9-43d8-8c4c-46d77c3b71e2","added_by":"auto","created_at":"2024-12-09 08:48:36","extension":"jpg","order_by":11,"title":"Figure 11","display":"","copyAsset":false,"role":"figure","size":264730,"visible":true,"origin":"","legend":"\u003cp\u003eLegend not included with this version.\u003c/p\u003e","description":"","filename":"figure11.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/0467120c66efe5211f035915.jpg"},{"id":70926530,"identity":"c928fda4-3719-488d-aa38-f13943edd755","added_by":"auto","created_at":"2024-12-09 09:12:40","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":5155083,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/8bd6c6b7-df06-434e-9a0e-648021f922c1.pdf"},{"id":70921570,"identity":"23efb966-5aa6-4dcc-878b-8f8b03abe55c","added_by":"auto","created_at":"2024-12-09 08:40:36","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":18449,"visible":true,"origin":"","legend":"\u003cp\u003eTable S1\u003c/p\u003e","description":"","filename":"Table1.docx","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/dd85e607e36dfefffb0caa8d.docx"},{"id":70923925,"identity":"002d8c12-b23a-4b6a-9223-b669faeeeabe","added_by":"auto","created_at":"2024-12-09 08:56:36","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":23261,"visible":true,"origin":"","legend":"\u003cp\u003eTable S2\u003c/p\u003e","description":"","filename":"Table2.docx","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/66deedaf0621876aadc37223.docx"},{"id":70921585,"identity":"3f4c9fe4-3eea-4490-af0b-4dbaad12e5e1","added_by":"auto","created_at":"2024-12-09 08:40:37","extension":"docx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":23869,"visible":true,"origin":"","legend":"\u003cp\u003eTable S3\u003c/p\u003e","description":"","filename":"Table3.docx","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/5ff05658c0796f7c96ba135e.docx"},{"id":70921571,"identity":"382ab059-211f-4ea0-97cf-e0872336ec4b","added_by":"auto","created_at":"2024-12-09 08:40:36","extension":"docx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":28537,"visible":true,"origin":"","legend":"\u003cp\u003eTable S4\u003c/p\u003e","description":"","filename":"Table4.docx","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/c443ec704b7223f0e66ee91e.docx"},{"id":70921573,"identity":"daf4d920-facd-49db-99e1-ff1f75c314e3","added_by":"auto","created_at":"2024-12-09 08:40:36","extension":"docx","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":21997,"visible":true,"origin":"","legend":"\u003cp\u003eTable S5\u003c/p\u003e","description":"","filename":"Table5.docx","url":"https://assets-eu.researchsquare.com/files/rs-5305987/v1/d35f4c74ccf9caeca405942f.docx"}],"financialInterests":"","formattedTitle":"Single-cell RNA landscape of the intratumoral heterogeneity and expression of angiogenesis-related genes in osteosarcoma","fulltext":[{"header":"Introduction","content":"\u003cp\u003eOsteosarcoma, one of the most commonly occurring and notoriously malignant bone tumors, has long been a focal point of biomedical research, particularly in elucidating its mechanisms of occurrence, progression, and metastasis[\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e]. Osteosarcoma, arising from primitive mesenchymal-derived osteoblasts, commonly manifests in bones undergoing rapid growth[\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. The treatment protocol for osteosarcoma involves surgical resection and chemotherapy, whereas radiotherapy is recommended in cases of unresectable osteosarcoma[\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. It exhibits a high proclivity for local invasion and early metastasis, which poses a threat to patient survival. Despite the development of anti-cancer therapeutics, the overall survival rate of osteosarcoma patients has improved, however, the prognosis remains poor for those with metastatic or recurrent osteosarcoma[\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. Therefore, researching new treatment strategies and developing new targeted drugs for osteosarcoma are paramount. However, the intricate mechanism of occurrence, development, and metastasis of osteosarcoma remain elusve. In recent years, the evolution of single-cell sequencing technology has remarkable breakthroughs in understanding of tumor tissues heterogeneity and the interaction between tumor cells and their microenvironment[\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. Significant progress has also been made in the study of osteosarcoma cells using single-cell sequencing technology[\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. This article aims to delve into the cellular heterogeneity among osteosarcoma subpopulations, angiogenesis-related gene expression patterns, intercellular communication networks, and crucial genes and pathways. Through a comprehensive analysis of osteosarcoma's single-cell data, we seek to provide a noval theoretical framework for understanding its pathogenesis and developing precision treatment strategies.\u003c/p\u003e \u003cp\u003eUtilizing single-cell transcriptome sequencing technology, this study conducted a comprehensive analysis of tumor samples from osteosarcoma patients. Through rigorous quality control measures and cell type annotation, we successfully identified multiple cell subpopulations and unveiled gene expression disparities among them. Notably, we uncovered the expression profiles of angiogenesis-related genes in these subpopulations and observed their elevated expression in specific cell types. This finding underscored the crucial role of angiogenesis in osteosarcoma's tumorigenesis and progression. Furthermore, leveraging bioinformatics techniques, we delved into the intercellular communication network and identified key genes. By quantifying the strength of intercellular communication, we discovered that certain cell types occupy pivotal positions within the communication network, serving as a decisive role in intercellular signaling. These insights not only enhanced our understanding of cellular interactions within the osteosarcoma microenvironment, but also presented potential therapeutic targets. Lastly, our analysis revealed significant differences in pathway activity among osteosarcoma cell subpopulations and the pivotal roles of certain genes. This comprehensive analysis offered a deeper understanding of osteosarcoma's biological complexities, paving the way for more targeted and effective therapeutic approaches.\u003c/p\u003e"},{"header":"Materials and Methods","content":"\u003cp\u003e\u003cstrong\u003eData download\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe Cancer Genome Atlas (TCGA)[7]\u0026nbsp;collects various human cancers and tumor subtype data, including clinical information, genomic variations, mRNA expression profiles, miRNA expression data, methylation, and more, which is an important data source for cancer researchers. TCGA also contains some data from TARGET. For our study, we retrieved reliable mRNA expression data in FPKM format of osteosarcoma and corresponding clinical information data, survival data, and copy number variation data from TCGA database. The samples in the data are all from Homo Sapiens, and the platform is based on Illumina. Additionally, single-cell transcriptome sequencing data of normal transcriptome and osteosarcoma samples were downloaded from GEO database. To ensure data consistency and comparability, we standardized the expression data using the limma package in R, employing a Log2 transformation[8], and the normalized expression data was visualized using box diagrams. There are 88 cases of TARGET osteosarcoma in the TCGA dataset, of which 85 cases with integrated clinical information. After integrating copy number variations, 26 cases were included in this study. Moreover, we also utilized the GSE16088 dataset from GEO database, which comprised 6 normal tissue samples and 14 osteosarcoma samples. Additionally, we accessed the GSE152048 dataset, which comprised single-cell transcriptome sequencing data of osteosarcoma. And 6 primary osteosarcoma samples of GSE152048 were selected for this study.\u003c/p\u003e\n\u003cp\u003eHUGO Gene Nomenclature Committee (HGNC)\u0026nbsp;[9]\u0026nbsp;is responsible for providing a unique, standard, and widely distributed symbol on protein-coding genes of the human genome. The mRNA expression profiles were obtained using HGNC mRNA gene annotation file.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSingle-cell data processed\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe processed the single-cell data for each sample used R package Seurat 4.3.0. Cells expressing fewer than 300 genes or with mitochondrial genes comprising more than 10% were filtered out from further analysis[10]. Subsequently, we utilized the R package DoubletFinder 2.0.3 to eliminate potential doublet cells from the sample[11]. Following this filtration, a total of 58,241 cells were retained for subsequent analysis. To process Seurat objects for each individual sample. we employed the Read10\u0026times; function. After normalization, the gene expression score was multiplied by 10,000 after adding 1 to prevent the logarithm from being 0, and then it was converted to natural logarithm values. Next, 3000 highly variable genes (HVGs) in the dataset were identified by using the \u0026quot;SelectIntegrationFeatures\u0026quot; function. We then scaled the data using \u0026quot;ScaleData\u0026quot; to control for the effects of sequencing depth and mitochondrial genes expression. To merge the samples and mitigate the batch effect, the R package Harmony 2.0 was utilized[12]. Subsequently, principal component analysis (PCA) was applied to identify significant principal components (PCs)[13]. The Elbowplot function was used to visualize the distribution of p-value. Finally, 30 PCs were selected for t-distributed stochastic neighbor embedding (tSNE) analysis. Using the default \u0026quot;FindNeighbors\u0026quot; parameter and 30 PCs dimension parameters, we constructed a k-nearest neighborhood based on Euclidean distances in PCA space. The Louvain algorithm, implemented through the \u0026quot;FindClusters\u0026quot; function, was then applied to optimize cell clustering. Using a resolution of 0.1, the cells were divided into 9 distinct clusters. Finally, we employed the \u0026quot;RunTSNE\u0026quot; function to perform dimensionality reducti, enabling us to visualize and analyze the dataset.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCell type identification\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eCell types can be identified through the utilization of cell type marker genes. For Osteoblastic\u0026nbsp;osteosarcoma, the marker genes consist of COL1A1, CDH11, and RUNX2. On the other hand, Osteoblastic_proli, which refers to proliferating osteoblastic osteosarcoma, is characterized by the marker genes PCNA and MKI67. Osteoclast cells are distinguished by the expression of CTSK and MMP9. Tumor Infiltrating Lymphocytes (TILs) are identifiable through the presence of IL7R, CD3D, and NKG7. Myeloid cells are marked by CD74, CD14, and FCGR3A. Fibroblasts are identified by the genes COL1A1, LUM, and DCN. Pericytes are characterized by the expression of ACTA2 and RGS5. Mesenchymal stem cells (MSCs) are recognized by the presence of CXCL12, SFRP2, and MME. Endothelial cells are distinguishable by the expression of PECAM1 and VWF[6]. To determine the differentially expressed genes (DEGs) among these cell types, we used FindAllMarkers and Wilcoxon rank-sum test to compare the gene expression profiles across different subpopulations.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDifferential expression of angiogenesis-related genes between cells\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe retrieved 36 angiogenesis-related genes from the HALLMARK_ANGIOGENESIS entry of \u0026quot;Hallmark Gene Set\u0026quot; in the MSigDB database[14]. Subsequently, we intersected them with the DEGs between cell types to identify angiogenesis-related genes that exhibited differential expression. To visualize the expression patterns of these genes, we employed the \u0026quot;DoHeatmap\u0026quot; function to generate a heat map. Additionally, we calculated the correlation among these genes using the \u0026quot;COR\u0026quot; function and depicted it as a correlation heat map.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAnalysis of intercellular communication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eBy leveraging the Python package CellphoneDB[15], we integrated single-cell expression profiles to calculate the number of interactions cell-to-cell communication and identify the interacting receptors and ligands. To visualize the interaction intensity among nine distinct cell types (Osteoblastic, Osteoblastic proli, Osteoclast, TIL, Myeloid cells, Fibroblasts, Pericyte, MSC, Endothelial cells), we generated a heat map. We used R packet CellChat[16] to calculate the communication patterns between cell subpopulations. Using CellChat to identify the significant interaction of ligand-receptor pairs through ligand-receptor interaction probability and perturbation test. The resulting cell-cell communication network was then constructed by integrating the number or strength of significantly interact ligand-receptor pairs across cell types. The interaction strength of osteoblastic, osteoblastic proli, osteoclast, TIL, myeloid cells, fibroblasts, pericyte, MSC and endothelial cells was shown through a circular graph. Subsequently, cellular communication mediated by intercellular ligand receptors was visualized through bubble plots. We also demonstrated several kinds of intercellular ligand-receptor interaction networks using the R-packet iTalk (https://github.com/Coolgenome/iTALK).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGSVA analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eGene Set Variation Analysis (GSVA)[17]\u0026nbsp;is a nonparametric and unsupervised algorithm that transforms gene expression data. We downloaded the C2.cp.kegg.v7.5.1.symbols.gmt dataset from Molecular Signatures Database (MSigDB). Using the \u0026quot;gsva\u0026quot; package in R, we analyzed the single-cell data of osteosarcoma to obtain the GSVA enrichment score for each cell corresponding to each pathway. Subsequently, we employed the R package limma 3.50.0\u0026nbsp;[8]\u0026nbsp;to identify pathways with significant differences (\u003cem\u003ep\u003c/em\u003e value \u0026lt; 0.05). The scores of pathway activity of cells in each group and other remaining groups were evaluated for difference test analysis. Finally, the Top3 pathways with t values arranged from highest to lowest in each group were plotted.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eProtein-protein interaction analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe inputed the differentially expressed TOP100 genes from each subpopulation into the String database[18]\u0026nbsp;to conduct protein interaction network analysis. During this pocess, we preserved the interaction relationships that had been experimentally validated. Following this, degree algorithm in Cytoscape was used to identify key genes among protein interactions. The top20 genes were then selected for visualization.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAbundance analysis of immune cell infiltration\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eCIBERSORT[19]\u0026nbsp;is a tool that employs linear support vector regression to deconvolute the expression matrix of human immune cell subtypes. By referencing the LM22 database, CIBERSORT quantify the relative expression abundance of 22 immune cell types, including B cells naive, B cells memory, plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells helper, T cells regulatory (Tregs), T cells gamma delta, NK cells resting, NK cells activated, monocyte, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, eosinophils and neutrophils. We analyzed osteosarcoma datasets from TARGET database to obtain the proportions of different immune cell types using CIBERSORT package in R.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eWGCNA analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWeighted Gene Correlation Network Analysis (WGCNA)[20]\u0026nbsp;aims to identify clusters of co-expressed gene, explore the association between gene networks and phenotypes, and investigate the core genes within the network. Firstly, we selected the scores of 22 immune cell type as the trait data for WGCNA analysis. Using the WGCNA package (version 1.71) in R, we calculated the soft threshold using pickSoftTreshold function. The optimal soft threshold was determined to be 3,\u0026nbsp;ensuring the construction of a scale-free network. Subsequently, based on the soft threshold, we constructed according a scale-free network. The networkwas then used to generate a topology matrix and perform hierarchical clustering. By setting the minimum number of genes in the module to 30, we dynamically cutted the modules and calculated the eigengenes. Eigengenes represent the overall expression profile of a module and are crucial for assessing the correlation between modules. After calculating the eigengenes, we constructed a correlation matrix between the modules and performed hierarchical clustering. The modules with a correlation above 0.25 were merged again, and finally 16 modules were obtained. Next, we calculated the correlation between gene modules and phenotypes using pearson to identify the modules that were associated with the phenotypes. The modules with the highest correlation coefficient were intersected with angiogenesis genes to obtain immune-related angiogenesis genes.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDifferential expression analysis of immune-related angiogenesis genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eObtained the expression matrix of OSTEOSARCOMA from the GSE16088 dataset, the limma package in R was employed to identify the DEGs between the tumor group and the normal group. The DEGs fulfilled the requirements of adj. p value\u0026lt;0.05 and | log2FC |\u0026gt;1. To visualize the differential expression patterns of these DEGs, the ggplot2 package was leveraged to generate a volcano plot, in which immune related angiogenesis genes were highlighted. Furthermore, the pheatmap package was utilized to create heat maps to visualize the expression variations of 10 immune-related angiogenesis genes between tumor and normal tissues in OSTEOSARCOMA patients.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGO analysis of immune-associated angiogenesis genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eGene Ontology (GO) analysis is a prevalent technique for large-scale functional annotation and enrichment[21], encompassing biological process (BP), molecular function (MF) and cellular component (CC). We used the R package clusterProfiler[22]\u0026nbsp;to perform GO annotation analysis of immune-related angiogenic genes. The screening standard q value was \u0026lt; 0.05, which was deemed statistically significant. Additionally, the Benjamini-Hochberg method (BH) for p-value correction was used.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsistent clustering for classification of osteosarcoma\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eR packet ConsensusClusterPlus[23]\u0026nbsp;was used to cluster the gene expression profile of osteosarcoma in TARGET database, focusing on a panel of 471 immune genes and angiogenesis genes. Spearman method was emploed to calculate the distance between genes, and the Partitioning Around Medoids (PAM) clustering algorithm was chosen for its robustness in handling outlier genes. Through a rigorous analysis of matrix heat map, consistency cumulative distribution function maps, and delta area plots, the osteosarcoma cells subtypes: Cluster1, Cluster2, Cluster3, and Cluster4, were identified. A box plot was generated to compare the expression levels of immune-related angiogenesis genes across these four subtypes, and the t-test was used for statistical significance. Significant differences were considered when P value \u0026lt; 0.05. To gain insights into the immune-related functions of these osteosarcoma subtypes, the immune.gmt dataset was downloaded from MSigDB. And then, the osteosarcoma data were analyzed using the \u0026quot;SSGSEA\u0026quot; method of GSVA package to obtain scores of immune-related functions of different samples. A box plot was constructed to visualize the differences in these immune-related functions among the four subtypes. Statistical significance was again determined using the t-test, with a P value \u0026lt; 0.05 deemed significant.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePrognostic modeling\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe conducted a prognostic analysis utilizing gene expression data and survival information from 85 osteosarcoma samples source from TARGET database. Initially, a univariate Cox regression analysis was conducted on 471 immune genes and angiogenesis genes, enabling us to preliminarly identify those genes that exhibited significant associations with overall survival (P value \u0026lt; 0.05). Subsequently, the osteosarcoma samples were stratified into two distinct sets: comprising 53 cases, used to establish a prognostic model, and a validation set consisting of 32 cases. Employing the prognostic genes selected through the univariate Cox analysis, we utilized lasso-Cox regression analysis to construct the prognostic model. The risk score was determined using the following calculation formula:\u003c/p\u003e\n\u003cp\u003e\u003cimg src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYIAAABLCAYAAAB9aoBBAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAAwlSURBVHhe7Z1PjBTFF8eL38ULogkcVMJBEyAeNBqiiQkYo5JA4EBAXCBx0QNi5KCIRi+CkQMc+HeBhKwXNbJk5UACK/EAAfWkUQnhoiZwYvUgB4Qjyf76U/Zba2u7e3p2epadqe8nqe3t6u6qV69e1Zt+1dMzZzzDCSGESJb/5VshhBCJIkcghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0fQB/z6669uYGDAbdiwwc2ZM8dt3bo1PyKEEK2RI+hx7ty54z744AM3MjLi3njjDTc2NubTn3/+mZ8hhBDVyBH0OHPnznVffvml27Ztm3vhhRfc/fff7x599NH8qBBCtEaOQAghEke/UNYHsCbwxRdfuAMHDvj9999/3w0ODrrPP//c7wshRBVyBEIIkTgKDQkhROLIEQghROLIEQghROLIEQghROLIEQghROLIEQghROLIEQghROLIEfQQBw8e9C+Vq5sWLFjgX0gnhBBV3BNHwDdhX375Zf/CtCr++usvt2vXLv9mzVbnit7kjz/+cM8+++yE4zp58mR+5D/I++yzz/K92QsyFslfBC8FpN2pO+q6c0GT9Io9QTs21Qmz+o7gww8/dIcOHXI3b97Mc8pBWUwkTCgMsG+++cYrsZ8G2osvvujmz5/v/+d1EnwpvCr9/fff7umnn/bn14FJ+c0335zQI1scMQ557969jQ9WJsOPP/7YffXVV+727dtu/fr1vs/CergL+umnn9ymTZvynNkLMiI/MrcCnd66dSvfmwo6YIK0u7uyVKcu8R+9ZE/Qjk11RDZhzGoGBwfHX3rppfFsoshzpnLmzJnxbIIcHx0d9fuc+9577/m8X375xef1C5kD4JUgPtHuphgeHvb62rZt2/jvv/+e5/6X36oPpgNtqSqX4/R/L0FbaFNZ39CeZ555xttlNtGPf/fdd36/rJ2ch/7RRQx5RfmimF60J2hlU03Q847AlBR3MPmvvvpq3zkCoK04gscee2x8bGwsz50+NtlUTUZMVk3UFVLVt9RlE2avUaWv7C5oPLvLHc/uvMYHBgb8ln3yi6AM+rlowufYJ598ku+JKnrZnqBbY9Dom8XiH374YdKPsfCe/h07duR7/cXRo0ddNoG6a9euuddee62jkA3X8sM28O677/ptDOGlNWvW5HvNQL03btzI96Zy4sQJ9/jjj7cV2potIDOy04YqHnjggfy/9jl79qz7+eef3Z49e/IcUUUv2xPUtanp0jVHEC70fv/99z5uT8yT/MOHD09aAMaoiUeTOLcsHk1eGDdloYkJf/Xq1X5SfOKJJyYtrDz//PNTOp54W7iWEK8hhGsNbDnfKGuTydqqbEBmk78scU4V9mM02SdFd/78+Y4cHusCly9fdk899ZRbvHhxnjsVJpyHH34436vWk1GmD/qbH9BBdhL/h3pke+7cOffkk0/6/RCOsY5h9VJWTFm9Yf/xv5XDOeGHCKjTl1UgO22wNhmse42Ojrq33nrL2yzb4eFhn98OX3/9df7fvz9VarKSiCeHdsZ+ke1a27AB6MS+q8Zw2TFSPBdA3THYqg8Nyu7EnqBTm+rUnqDMphohvzNoFAvXUDy3M8TuiTVzizs0NOTzLSTArU5mbBO3bNwCh+GCOHzA9dwOx+EE1gQoN6wzhnOsLBL/h+EVq9ti5MhMechg5xe1iespO6yTsrq9RkHZ1IFMRaGDOhB3tDbWpUpPRit9mD5JcV9au4pioqE9WL0kyyurl1j8+vXr/blLlizx4RjOpw30YTuy1wHZq64hHzlalYltIZ+1M0yhfmgHOrA6rXzaQTut7eHYsLZzHaEptuE5de27agxXHTPbs/3weNkYrNuHIdSNrNOxJ+jUppqwJ2hlU53Q1TUCFGGGFBIqn0aZwQJ5n3766UQnhOeiUDqrDDoAIy7qTKsnVCKKtTwbcLGxUH+o/KI2cYxzrN4w0endxAYTqcjQW2HX06461NFTHX3QL/RP2EcGZRcZvNUdl2H7deoN7QmsDGt/U31p5ZT1CcdjeywibjMgM/LGZdu5jIHsE+okG4W47WD9b2VxTrv2bW0pGsNVxyCUyeSP28U5oT3E7WAb9mEM5YXXG7FurZxQ163aDlXy1Lm+LlZWrJ8m6PoaAb+fy61/GYQjFi1a5GPQS5cu9bdmO3fu9OEP459//vG3uo888kjlY1+U9eOPP7pMYS4bDJPCJsSj79696x566CG/D2vXrnW//fabDx8Rb+VWPWbjxo3+8dUwnh23iWPEezPD8o9thonbxpAmQkMhtCEzOP//O++8M3GbX5eFCxf6R1JpQ51bzjp6akcfTXDlyhW/7aRea/9MyY7Nme21C2MDfccQujt16pTvn/vuu29SKK8M639kMdq176oxXGd8G+2MwSKsDzvF7Ak6sQeuRa8zORamyz1fLMYgTp8+PfF9gc2bN7t169ZN6tB58+Z5Y923b9+UGB7nxWsKDK4LFy74BdXQODAyjK0b8Ew4McJW8PORsUHEqZ2fmCQOyUI5LF++vDLOXwTnsz5w/fp1/yx/U9TVRzswsb399tvu2LFj3uHRduQOY79N1NsN2ZuGDwCkGOydCZefLi2Ld4fwwYiJimuqqNJJ1RiuM77vFXXsCTq1h16wp3vuCFA+i1l8SuALUMPDw37x8uLFi/kZ/4Ih8ani9ddfn7LQcvXq1SmfhDFAPu0Y9skn/tTMYCHZ8XARzshuH92yZcvyvalYPStXrvRfZDOQM1zo6gYfffSRd3A4PZ4mahf0xNNClEFZZYyMjPi+qqOnTvVh1zOpxWzZssUtWbLEPffcc37LpLJ9+3Z/rIl+qFMGemi14GeyW3ndwuwXkIeHCL799lt/l1g0VmL4YMREVSVnK51UjeG64xuoZ7pjsAqTv117glZtb0Wd6+vYE3TVprJPoF0jjp0ZYX52y+TjdBb7Jw5GjJMthOda7C0zlonjYV64fkBMkrwwnkZZNDlM1I0MQMyOPCuHNQfKDmN5VW2Kyw7l7AYmb9iG6WJlxV8o43/yQj3W1ZPpwVJRvxXp0o6F5RmsE1XptFW9cf8VydFEXyJ7UdvahTqpu0gX2Dhjhb4nsS5gfcc+x7ANy6NdlBX3G/lGrB+jSifUVTaGq45BXF9d2wqvYct+kdxgx6djT1DVdjteJU+r6+vSlE0V0RVHYIqwRhdNtiQURH7mGf3TQCjHDLWoDDOKOI/Fp8yr+lV7rueYlRNCmXS8XR8apGFyxGVUtQk4ztNMdm1R2U3CxGwyNlUPE8uqVasm2kiKHYNRpiejSh8me5iwhRDsJM6zSSW+lrKtL8rqLeo/8sPyrE+b6EtkL5p46hLLW5ZMR2yL9i2hc/Iokz4lj/bZE3hF+qlr35xXNIZbHYvnAqOdMVjWhzHTtScoa3tdeRg/ndoTIH8nNlVFV+8IRHfAiMyoumUY9xoGYtGACR15mGZSD8jGI5FFEw4gczyZzAbMETCBpUYv2xN026b65pvFqUA88ZVXXvELb9nAnlVPHjQJC3m7d+92+/fvn1hYZPvggw/6Re3MdidSNkj8kykzha3LlHHkyBH/BaM6T+2ImaGX7Qm6bVNyBD0EhssrJTpZHO4leCqGRTR7yuT48ePu0qVL3hka5DNwV6xYked0H3vFRwyy4Jj5Bmm/OuheptfsCWbMpjIPKHoEblfpsrI4aL9CjJi4MbFWi29b4tZ+pkMd1IcccR8gY7xWMhtA3jiWnZL9xPSKPcFM2ZQcQQ9BjDc02jopfNpHNIM55HjxUYjpMBvsaQ5/MiFED8A3jvmiUDtkjqDwi0dCCGHIEQghROJosVgIIRJHjkAIIRJHjkAIIRJHjkAIIRJHjqBPqPsGQyGEiNFTQ0IIkTi6IxBCiMSRI+gTCA3x4xfhe1OEEKIOcgR9Qp03GAohRBFyBH1C1RsMhRCiCjkCIYRIHDmCPoF3qw8NDVX+AL0QQhShx0eFECJxdEcghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBCJI0cghBBJ49z/AeXYeqjTOQ55AAAAAElFTkSuQmCC\"\u003e\u003c/p\u003e\n\u003cp\u003eCoef (GeneI) represents lasso-Cox regression coefficient; Expression (Gene\u003csub\u003ei\u003c/sub\u003e) represents the expression value of each gene, and n represents the number of genes[24].\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e85 osteosarcoma samples were categorized into high-risk and low-risk groups based on the median risk score derived from the training set. To assess the prognostic value and predictive accuracy of the model, we employed Kaplan-Meier survival curve analysis to analyze overall survival and generated time-dependent receiver operating characteristic curves (ROC).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConstruction and correlation analysis of Nomogram model\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAfter removing clinical features with null from the TARGET osteosarcoma database using R, we obtained three characteristics: age, gender and M stage. To further investigate factors related to patient prognosis, we conducted both univariate and multivariate Cox risk regression analysis, incorporating the prognostic risk score alongside these clinical characteristics. For visualization, we employed the R package \u0026apos;forestplot\u0026apos;. Leveraging the nomogram function within R package RMS, we constructed a nomogram model specific to the prognostic factors significantly associated with outcomes in the TARGET osteosarcoma database. For enhanced visualization, we utilized the \u0026apos;ggplot\u0026apos; package. Time-dependent ROC analysis was performed on the predicted scores from nomogram model, considering the overall survival status and survival time of patients at one, three and five years, respectively, using R package timeROC. Multivariate ROC analysis of nomogram model was performed incorporating the predicted score, overall survival status and survival time of patients.\u003c/p\u003e\n\u003cp\u003eFinally, we evaluated the accuracy and resolution of the nomogram using a calibration curve. We employed the R package RMS to devise the nomogram and calibration curve. We leveraged R package ggDCA to evaluate the one-year, three-year, and five-year survival rates of patients using the nomogram model.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eImmune infiltration analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe abundance of 22 immune cell types was calculated based on the immune invasion matrix of osteosarcoma samples from TARGET. We analyzed the association between these immune cells and prognostic genes. In addition, t-test was used to compare the abundances of the 22 immune cell types between the high-risk and low-risk groups, with p value \u0026lt; 0.05 considered as statistically significant. The correlation between the expression of prognostic genes and the infiltration of immune cell was also calculated.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCopy number analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eCopy number data of osteosarcoma were downloaded from TARGET database. A copy number of 0 was defined as double deletion. A copy number of 1 was designated as single deletion. A copy number of 2 was considered normal. A copy number of 3 was defined as single gain. And a copy number of 4 or more was classified as amplification.\u0026nbsp;Then a boxplot was drawn to compare the differences in copy numbers and gene expression.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eStatistical analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWith the exception of CellphoneDB, which employed Python, all remaining data calculations and statistical analysis were carried out utilizing R language. When comparing continuous variables across two groups, the statistical significance of normally distributed variables was determined using the independent Student\u0026rsquo;s t test. For variables that did not follow a normally distributed, differences were assessed through the Mann-Whitney U test. To evaluate the statistical significance of categorical variables between the two groups, we employed either the Chi-square test or Fisher\u0026apos;s exact test, depending on the circumstances. Furthermore, the correlation coefficients among different genes were derived through Pearson correlation analysis.\u003c/p\u003e"},{"header":"Result","content":"\u003cp\u003e\u003cstrong\u003e1. Cell heterogeneity in osteosarcoma\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe conducted a single-cell RNA sequencing (scNA-seq) analysis on tumor samples from 6 osteosarcoma patients within the GSE152048 dataset to delve into the cellular composition. After rigorous quality control, we filtered out cells with mitochondrial gene content \u0026gt;10%, those with feature counts \u0026lt;300 or \u0026gt;6000, and eliminated duplicates, resulting in a final dataset of 58,241 cells for further analysis. Nine cell types were identified by artificial annotation based on gene expression and marker genes (Figure 2A). The marker genes utilized for this classification were detailed in table S1. Cell subpopulations revealed the following distribution: osteoblastic cells (18014,30. 93%), osteoblastic_proli (3172, 5.44%), osteoclasts (6194, 10.64%), TILs (318, 5.53%), myeloid cells (15309, 26.29%), fibroblasts (6390, 10.97%), pericytes (1844, 3.17%), MSCs (115, 1.91%), endothelial cells (2984, 5.12%). The differential expression of 23 marker genes that distinguish these subpopulations was visualized in figure 2B (violin diagram) and figure 2C (bubble diagram), demonstrating the preferential expression of these marker genes in respective cell subpopulations. We calculated the DEGs between the cell subpopulations using the \"FindAllMarkers\" function. Figure 2D showed the top three DEGs for each of the nine main cell subpopulations. Additionally, we tabulated the number and proportion of each cells subpopulation across the 6 samples. Notably, osteoblastic cells comprised a higher proportion in BC5 and BC6, myeloid cells were more abundant in BC2 and BC16, and BC22 exhibite a higher content of fibroblast cells (Figure 2E and Figure 2F).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e2. Differentially expressed angiogenesis-related genes among cell subsets\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe identified 19 differentially expressed angiogenes-associated genes (AAGs) by intersecting DEGs among cell subsets with angiogenesis-related genes. A heat map was used to visualize the expression of 19 angiogenesis-related genes across various cell subsets (Figure 3A). The figure illustrated that angiogenesis-related genes were significantly upregulated in osteoblastic osteosarcoma and proliferative osteoblastic OSTEOSARCOMA, whereas their expression was relatively low in myeloid cells. Furthermore, MSCs also exhibited high expression levels of angiogenesis genes, including COL3A1, COL5A2 and FSTL1. Additionally, certain AAGs, such as KCNJ8, NRP1, POSTN, TIMP1 and VCAN were expressed in pericytes. Notably, S100A4 and SPP1 were highly expressed in osteoclast. Our correlation analysis of these 19 angiogenesis-related genes revealed that THBD, NRP1, SLCO2A1, JAG2 and STC1 were positively correlated, whereas VCAN, COL5A2, FSTL1, FGFR1, TIMP1, POSTN, COL3A1 and LUM also demonstrated positive correlations (Figure 3B). These genes were further visualized in the tSNE dimension reduction diagram for osteoblastic and osteoblastic1 proli cells (Figure 3C).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e3. Cell communication and hub genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe inferred and quantified the communication among 9 cell subtypes using CellphoneDB and CellChat respectively, and then visualized cell communication intensities through a circle diagram and a heat map (Figure 4 a-b). A cursory examination revealed that osteoblastic cells, fibroblasts and MSCs exhibited a relatively higher cumulative signal intensity, indicating a relatively higher cumulative signal intensity, indicating a greater level of communication activity among these cell types. Similarly, fibroblasts, osteoblastic proli cells, osteoclasts and MSCs also demonstrated a significant cumulative signal intensity. Subsequently, we calculated all important ligand-receptor pairs that mediated from TIL to various cell types, including fibroblasts, osteoblastic proli, myeloid cells, osteoclasts and MSCs (Figure 4C). Notably, CD47-related pathway emerged as a pivotal communication channel between TIL and myeloid cells. The communication intensity between TIL and osteosarcoma cells appeared to be relatively low, with CD44 standing out as the most important ligand-receptor pair. ITGB1 and SPP1 related ligand-receptor pathways played an important role in the communication between fibroblasts, osteoblastic proli, osteoclasts and MSCs (Figure 4D).\u003c/p\u003e\n\u003cp\u003eSubsequently, GSVA analysis was performed on single-cell data of osteosarcoma to derive the GSVA enrichment score for each cell across various pathways. Then, the pathway activity among diverse cell subsets was evaluated to identify those exhibiting significant differences (P value \u0026lt; 0.05). Furthermore, the differential analysis was performed to ascertain the pathways that demonstrated notable variations in the relative pathway activity across different cell subsets. A heat map was generated to visualize the top 27 pathways from each group, arranged in descending order of each group (Figure 4E). Additionally, protein interaction analysis was also conducted focusing on the top100 DEGs across various subpopulations, and the top 20 hub genes were obtained: FN1, PTPRC, VEGFA, IL1B, COL1A1, MMP9, CD8A, JUN, ITGB1, PECAM1, CXCL8, COL1A2, MMP2, CCL2, COL3A1, VWF, SPP1, FCGR3A, ITGB2, and TYROBP (Figure 4F).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e4. WGCNA analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe performed WGCNA analysis of osteosarcoma data obtained from the TARGET database. We used immune-related genes and 36 angiogenesis-related genes as the gene set for WGCNA analysis. Immune-related genes were obtained from Immport database. Initially, we determined that a soft threshold β was 3 (R-square= 0.94) would optimize the consistency of the network with the characteristics of scale-free network (Figure 5A). Next, we gerneated a hierarchical clustering tree using the correlation coefficient between genes, leading to the identification of 16 similar gene modules (Figure 5B). A heat map was constructed to visualize the correlation of modules. We employed CIBERSORT to quantify the abundance of 22 immune cell types in osteosarcoma, considering them as phenotypic traits. Among the 10 modules, the MEred module exhibited the highest correlation with CD8+ T cells, while the MEturquoise module showed the strongest association with resting dendritic cells (Figure 5C). These two modules collectively encompassed 471 genes, including 10 angiogenesis-related genes (FGFR1, LPL, LRPAP1, KCNJ8, COL3A1, MSX1, JAG2, VAV2, PGLYRP1, APOH), all of which had immune-related angiogenesis functions.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e5. Differential expression of hub genes in angiogenesis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe GSE16088 dataset was used to conduct a differential expression analysis aomparing the tumor group with the normal group. We employed R package limma to analyze the transcriptome data of OSTEOSARCOMA patients, utilizing |log2FC| \u0026gt; 1 and adj.p value \u0026lt; 0.05 as threshold criteria for screening. We identified 2270 upregulated genes and 1981 downregulated gene. Volcano and heat maps depicting these DEGs were presented in Fig 6A and Fig 6B, respectively, with 10 immune-related angiogenesis genes highlighted in yellow in Fig 6A. The finding revealed that FGFR1, LPL, COL3A1 and MSX1 were highly expressed in tumor tissues, whereas APOH was low expressed in tumor tissues. The detailed expressions of these 10 angiogenesis-related genes were shown in Table S4. Subsequently, we conducted a GO enrichment analysis on these 10 genes, revealed that their functions were primarily concentrated in BP, such as positive regulation of lipase activity, regulation of plasma lipoprotein particle levels, blood coagulation, skeletal system morphogenesis, hemostasis. At the CC level, they are associated with chylomicron, very high-density lipoprotein particle, triglyceride plasma lipoprotein particle, plasma lipoprotein particle, lipoprotein particle, protein-lipid complex. In terms of MF, these genes were involved in glycosaminoglycan binding, heparin binding, sulfur compound binding, growth factor binding, 1-acyl-2-Lysophosphatidylserine acylhydrolase activity, lipoprotein lipase activity.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e6. Classification of osteosarcoma\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe TARGET osteosarcoma samples were categorized through consistent clustering. Examination of the consistent clustering heat map (Figure 7A), consistent cumulative distribution function (CDF) curve (Figure 7B) and delta area curve (Figure 7C) revealed that the optimal clustering number was 4. Consequently, osteosarcoma was stratified into four subtypes: cluster1, cluster2, cluster3 and cluster4. A comparative analysis of the expression of 10 immune-related angiogenesis genes across these subtypes (Figure 7D) revealed significant differential expression among them (p value \u0026lt; 0.05). Notably, APOH exhibited differential expression not only between Cluster1 and Cluster2 but also between Cluster1 and Cluster3. Futhermore, we assessed the abundance of 22 immune cell types across the four osteosarcoma subtypes. The cell types that exhibited notable abundance differences among the four subtypes activated CD4 T cells, activated CD8 T cells, gamma delta T cells, MDSCs, macrophage, mast cells, monocytes and regulatory T cells. Among them, cluster3 osteosarcoma, which was characterized by higher abundance of activated CD8 T cell, gamma delta T cells, MDSCs and regulatory T cells, was identified as immune infiltrating osteosarcoma.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e7. Prognostic modeling\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eUtilizing gene expression data from TARGET database for osteosarcoma, 59 genes associated with overall survival were identified by univariate Cox regression analysis from 471 genes related to immunity and angiogenesis within MEred and MEturquoise modules. The 85 samples were stratified into a training set (consisting of 53 cases) and a validation set (consisiting of 32 cases). Employing Lasso-Cox regression analysis on the training set samples, we calculated the Cox regression coefficients for the 59 genes (Figure 8a-b). Notably, eight genes - CD48, ESRRA, GNAI1, GNRH1, JAG2, KRAS, SECTM1, and TRPC4AP - emerged as significant predictors with non-zero coefficients. The risk score for each sample was calculated according to the following formula: risk score = -2.034× Exp(CD48) + (2.969) × Exp(ESRRA) + 1.177× Exp(GNAI1) + (2.492) × Exp(GNRH1) + (0.907) × Exp(JAG2) + (-1.600) × Exp(KRAS) + (-0.779) × Exp(SECTM1) + (-3.855) × Exp(TRPC4AP). Based on the median risk score derived from the training set, the osteosarcoma samples in TARGET database were categorized into the high-risk and low-risk groups (training set: 26 high-risk, 27 low-risk; validation set: 20 high-risk and 12 low-risk). As evident from the distribution of risk scores and corresponding survival status for both the training set (Figure 8C) and validation set (Figure 8D), an increase in risk score correlated with a higher risk of mortality and shorter survival time. The prognostic model's risk score demonstrated excellent predictive performance, with 1-year, 3-year, and 5-year area under the curve (AUC)values exceeding 0.65 for both the training set (Figure 8E) and validation set (Figure 8F). Survival analysis confirmed a significant difference in overall survival between the high-risk and low-risk groups in both the training sets (Figure 8G) and validation sets (Figure 8H) (P value \u0026lt; 0.05).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e8. Verify the prognostic independence of the model and the prognostic significance of other risk factors\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe results of the univariate Cox regression analysis from the TARGET database indicated that age (HR = 1.009, 95% CI = 0.894-1.139, P = 0.880), risk score (HR = 1.004, 95% CI = 1.000-1.008, P = 0.068), M stage (HR = 1.719, 95% CI = 0.444-6.654, P = 0.433), and gender (HR = 1.120, 95% CI = 0.316 - 3.973, P = 0.880) did not show significant associations with prognosis. However, it is noteworthy that the P value for risk score was close to significance (P = 0.068). (Figure 9A) (Table S5). A multivariate Cox regression analysis was conducted, employing these factors as covariates to further assess their prognostic significance. The analysis revealed that age (HR = 1.020, 95% CI = 0.882-1.180, P = 0.788), gender (HR = 1.427, 95% CI = 0.301-6.780, P = 0.654), and M stage (HR = 2.186, 95% CI = 0.502-9.517, P = 0.297) were not independent prognostic factors for osteosarcoma. However, the risk score emerged as a significant independent prognostic factor for overall survival among osteosarcoma patients (HR = 1.005, 95% CI = 1.000 - 1.010, P \u0026lt; 0.05) (Figure 9B) (Table S5).\u003c/p\u003e\n\u003cp\u003eSubsequently, we further constructed a nomogram model to predict 1-year, 3-year and 5-year overall survival, which incorporated four factor: age, gender, M stage and prognostic risk score (Figure 9C). By leveraging the R package timeROC, we evaluated the predictive performance of the nomogram model by calculating AUC values for 1-year, 3-year, and 5-year survival predictions. The all AUC values exceeded 0.9, and outperformed the prognostic risk score alone. In addition, we conducted prognostic calibration analysis for 1-year, 2-year, and 3-year survival probabilities using variables from both univariate and multivariate Cox regression modles. And calibration curves were plotted (Figure 9D).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e9. Immunocorrelation analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe results of immunocorrelation analysis of the TARGET osteosarcoma dataset indicated a rich immune cells infiltration within the samples. We calculated the correlation between immune scores and immune cells infiltration, discovering that the infiltration of macrophages M1, T cells follicular helper, T cells CD8, and Tregs exhibited correlation. A correlation was observed between the infiltration of B cells naive and plasma cells was correlated. T cells CD4 memory resting and NK cells resting was correlated. We also analyzed the differences in the abundance of immune cells between the highrisk and low-risk groups. As shown in Figure 10A and Figure 10C, there were significant differences in the abundance of T cells CD8, plasma cells, T cells CD4 naive, T cells Helper, T cells regulatory (Tregs), monocyte, macrophages M1, and dendritic cells resting between the two group (P value \u0026lt; 0.05). The abundance of T cells CD8 was lower in the high-risk group. Next, we investigated the relationship between prognostic gene expression and T cells CD8 infiltration. The expression of CD48 and SECTM1 was positively correlated with T cells CD8 infiltration (Figure 10D and Figure 10G), while the expression of GNAI1 and JAG2 was negatively correlated with T cells CD8 infiltration (Figure 10E and Figure 10F).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e10. Differential genes and CNV analysis in highrisk and low-risk groups\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe model was employed to assess the prognostic risk of each sample in the TARGET osteosarcoma database. The cutoff value was set as the median prognostic risk score across all samples. Using this cutoff value, the samples were stratified into high-risk and low-risk groups, and the DEGs between these two risk groups was calculated employing the R package limma. A total of 113 genes were overexpressed in the low-risk patients and 36 genes were overexpressed in the high-risk patients. We used a heat map to show the expression of TOP30 differential expressed genes (Figure 11A). Subsequently, KEGG pathway enrichment analysis was performed on the TOP30 genes. The results indicated that these DEGs were enriched in various BPs, including complement and coagulation cascades, staphylococcus aureus infection, pertussis, rheumatoid arthritis, asthma (Figure 11C). We also correlated the expression of TOP30 genes with the results of the KEGG enrichment analysis, preenting them in a string diagram (Figure 11B) and a network diagram (Figure 11D) respectively. Finally, we explored the correlation between TOP30 gene expression and copy number variation. Our analysis revealed that the expressions of CGREF1, NGEF, PDGFD, RHBDL2, TCN2 and TRAC were associated with copy number amplification, and most of the genes with significant associations exhibited copy number amplification of 4 or higher, namely, indicating a strong link between gene expression and amplification.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eOsteosarcoma, a notoriously aggressive bone tumor, has long posed significant therapeutic challenges[\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e, \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e, \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]. However, the remarkable advancements in single-cell sequencing technology have revolutionized our understanding of the cellular composition, gene expression patterns, and intercellular dynamics of osteosarcoma. This intricate level of analysis has provided profound insights into its pathogenesis and offers hope for the discovery of novel therapeutic strategies.\u003c/p\u003e \u003cp\u003eBy analyzing the cellular heterogeneity, we have been able to delineate diverse 9 cell subpopulations from osteosarcoma samples within the GSE152048 dataset. These subpopulations, which include osteoblastic cells, osteoblastic_proli, osteoclasts, TILs, myeloid cells, fibroblasts, pericytes, MSCs and endothelial cells. Osteoblastic cells comprised a higher proportion in BC5 and BC6, myeloid cells were more abundant in BC2 and BC16, and BC22 exhibite a higher content of fibroblast cells. This heterogeneity not only reflects the diverse nature of the cell types involved but also indicates varying levels of gene expression across each subpopulation. This cellular diversity was presumably a key factor underlying osteosarcoma's resilience to traditional therapeutic approaches. Then, top three DEGs for each of the nine main cell subpopulations were analyzed. Consequently, the development of targeted therapeutic strategies, specifically tailored to address these distinct cell subpopulations, could potentially lead to more favorable treatment outcomes.\u003c/p\u003e \u003cp\u003eOur meticulous examination of gene expression patterns has been particularly crucial in elucidating the role of angiogenesis, a crucial process in tumor growth and metastasis. We identified 19 differentially expressed AAGs by intersecting DEGs among cell subsets with angiogenesis-related genes. AAGs were highly expressed in osteoblastic osteosarcoma and proliferative osteoblastic osteosarcoma, whereas their expression was relatively low in myeloid cells. Osteosarcoma, as a highly vascularized malignancy, heavily relies on angiogenesis for its aggressive nature [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. Our findings have revealed that genes related to angiogenesis were abundantly expressed in osteoblasts and their proliferative counterparts. MSCs exhibited high expression levels of FSTL1. Blocking FSTL1 could treat osteosarcoma effectively by inhibiting cell proliferation, invasion, and various other crucial cellular functions[\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]. NRP1, POSTN, TIMP1 and VCAN were expressed in pericytes. MicroRNA-1247 inhibited the viability and metastasis of osteosarcoma cells via targeting NRP1 and mediating Wnt/β-catenin pathway[\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]. Our correlation analysis of these 19 angiogenesis-related genes revealed that THBD, NRP1, SLCO2A1, JAG2 and STC1 were positively correlated, whereas VCAN, COL5A2, FSTL1, FGFR1, TIMP1, POSTN, COL3A1 and LUM also demonstrated positive correlations. This discovery aligns with prior research and suggests that suppressing the expression of these angiogenesis-associated genes could potentially disrupt the angiogenesis process in osteosarcoma, thereby inhibiting tumor growth and metastasis.\u003c/p\u003e \u003cp\u003eFurthermore, our study employed bioinformatics techniques to delve deeper into the intricate intercellular interactions within osteosarcoma. Osteoblastic cells, fibroblasts and MSCs exhibited a relatively higher cumulative signal intensity, indicating a relatively higher cumulative signal intensity, indicating a greater level of communication activity among these cell types. Similarly, fibroblasts, osteoblastic proli cells, osteoclasts and MSCs also demonstrated a significant cumulative signal intensity. Intercellular communication plays a pivotal role in the tumor microenvironment, serving as a crucial conduit for signaling among diverse cell types and orchestrating tumor growth, invasion, and metastasis. Through a meticulous analysis of the communication network and the identification of key ligand-receptor pairs, we gained profound insights into the complex interplay between osteosarcoma cells and their surrounding microenvironment. This understanding could pave the way for the development of innovative therapeutic strategies that specifically target these interactions, thus halting tumor progression. GSVA analysis was performed to derive the GSVA enrichment score for each cell across various pathways. Additionally, protein interaction analysis was also performed on the top100 DEGs across various subgroups. WGCNA analysis of osteosarcoma data was performed to obtain the immune-related genes. The MEred module exhibited the highest correlation with CD8\u0026thinsp;+\u0026thinsp;T cells, while the MEturquoise module showed the strongest association with resting dendritic cells. Research has demonstrated that IL-35 could lead to dysfunction of CD8\u0026thinsp;+\u0026thinsp;T cells and thus constrain the anti-tumor immune response in osteosarcoma[\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]. This study indicated that FGFR1, LPL, COL3A1 and MSX1 were highly expressed in osteosarcoma tumor tissues, whereas APOH was low expressed in osteosarcoma tumor tissues. The activation of nuclear FGFR1 induced histone H3 phosphorylation at Ser 10 and c-jun/c-fos expression to contribute osteosarcoma cell survival rendering radiation resistance[\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. Research has unequivocally confirmed that the miR-29 family exerts a tumor suppressive function in modulating MTX resistance and osteosarcoma cell apoptosis, accomplished through the regulation of COL3A1[\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eIn this study, osteosarcoma was stratified into four subtypes: cluster1, cluster2, cluster3 and cluster4, and the abundance of 22 immune cells types across the four osteosarcoma subtypes was evaluated. Cluster3 osteosarcoma, which was characterized by higher abundance of activated CD8 T cell, gamma delta T cells, MDSCs and regulatory T cells, was identified as immune infiltrating osteosarcoma. Utilizing gene expression data from TARGET database for osteosarcoma, eight genes - CD48, ESRRA, GNAI1, GNRH1, JAG2, KRAS, SECTM1, and TRPC4AP - emerged as significant predictors with non-zero coefficients. The expression of high levels of JAG2 in osteosarcoma patients could potentially be linked to a more favorable response to immune checkpoint blockade therapy [\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e]. Hsa-miR-557 effectively inhibited osteosarcoma growth both in vivo and in vitro by modulating the expression of KRAS[\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eThe univariate Cox regression analysis did not detect any significant associations between age, risk score, M stage, and gender with prognosis. However, upon further investigation using multivariate Cox regression analysis, it became evident that age, gender, and M stage were not independent prognostic indicators for overall survival. Notably, the risk score emerged as a significant and independent prognostic indicator for overall survival among osteosarcoma patients. This finding suggests that the risk score, likely encompassing multiple clinicopathological variables, may hold key information for predicting patient outcomes and guiding treatment strategies in osteosarcoma.\u003c/p\u003e \u003cp\u003eIn our analysis, we calculated the correlation between immune scores and the infiltration of immune cells. Our findings revealed a significant correlation between the infiltration of macrophages M1, follicular helper T cells, CD8\u0026thinsp;+\u0026thinsp;T cells, and regulatory T cells. This suggests that these immune cell subsets play a pivotal role in modulating the immune response in the context of our study. Further investigation is needed to understand the functional significance of these correlations and their potential implications in disease pathogenesis and therapeutic responses.\u003c/p\u003e \u003cp\u003eThe results of our analysis indicate that the expressions of CGREF1, NGEF, PDGFD, RHBDL2, TCN2, and TRAC are closely associated with copy number amplification. Furthermore, it is noteworthy that the majority of genes displaying significant associations exhibited copy number amplification of 4 or higher. This observation strongly suggested a robust correlation between gene expression levels and copy number amplification in these specific genes. Such a link could have significant implications for understanding the molecular mechanisms underlying osteosarcoma and other related malignancies, as well as for developing potential therapeutic targets. Future studies should aim to further explore the functional roles of these genes and their association with disease progression and prognosis.\u003c/p\u003e \u003cp\u003eIndeed, there are inherent limitations in this study that need to be acknowledged. Firstly, the absence of validation using a larger sample size is a significant limitation. This could have affected the reliability and generalizability of our findings, as a larger sample would provide a more robust basis for validating the observed associations. Secondly, the lack of further evidence from basic experiments is another shortcoming. Basic experiments would have provided deeper insights into the biological mechanisms underlying the observed gene expression patterns and copy number amplifications. Such evidence is crucial for validating and complementing the findings obtained from bioinformatic analysis. In future studies, it is recommended to address these limitations by including a larger sample size for validation and conducting basic experiments to further explore the functional roles of the identified genes.\u003c/p\u003e \u003cp\u003eIn summary, the utilization of single-cell sequencing technology in osteosarcoma research has significantly enhanced our comprehension of the disease's cellular landscape, gene expression patterns, and intercellular dynamics. These revelations provide precious leads for devising more effective and tailored therapeutic approaches against this aggressive bone tumor. Future endeavors in this field promise to further improve treatment outcomes for osteosarcoma patients.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate:\u003c/strong\u003e Not applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003cstrong\u003e：\u003c/strong\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of data and materials:\u003c/strong\u003e The data that support the findings of this study were obtained from TCGA and GEO. Derived data supporting the findings of this study are available from the corresponding author on reasonable request.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests:\u003c/strong\u003e The authors declare that we have no commercial or financial relationships that could potentially be construed as a conflict of interest.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding:\u003c/strong\u003e This study was supported by National Natural Science Foundation(82273436).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors' contributions:\u003c/strong\u003e HL analyzed the data and wrote the paper. TM and ZQY provided the help of the R language. ML, FG and XJL designed the project. ML selected the analyzed results. All authors read and approved the final manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments:\u003c/strong\u003e Not applicable.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eChen C, Xie L, Ren T, Huang Y, Xu J, Guo W: \u003cstrong\u003eImmunotherapy for osteosarcoma: Fundamental mechanism, rationale, and recent breakthroughs\u003c/strong\u003e. \u003cem\u003eCANCER LETT \u003c/em\u003e2021, \u003cstrong\u003e500\u003c/strong\u003e:1-10.\u003c/li\u003e\n\u003cli\u003eMutsaers AJ, Walkley CR: \u003cstrong\u003eCells of origin in osteosarcoma: mesenchymal stem cells or osteoblast committed cells?\u003c/strong\u003e \u003cem\u003eBONE \u003c/em\u003e2014, \u003cstrong\u003e62\u003c/strong\u003e:56-63.\u003c/li\u003e\n\u003cli\u003eCiernik IF, Niemierko A, Harmon DC, Kobayashi W, Chen YL, Yock TI, Ebb DH, Choy E, Raskin KA, Liebsch N\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eProton-based radiotherapy for unresectable or incompletely resected osteosarcoma\u003c/strong\u003e. \u003cem\u003eCANCER-AM CANCER SOC \u003c/em\u003e2011, \u003cstrong\u003e117\u003c/strong\u003e(19):4522-4530.\u003c/li\u003e\n\u003cli\u003eWedekind MF, Wagner LM, Cripe TP: \u003cstrong\u003eImmunotherapy for osteosarcoma: Where do we go from here?\u003c/strong\u003e \u003cem\u003ePEDIATR BLOOD CANCER \u003c/em\u003e2018, \u003cstrong\u003e65\u003c/strong\u003e(9):e27227.\u003c/li\u003e\n\u003cli\u003eLee HW, Chung W, Lee HO, Jeong DE, Jo A, Lim JE, Hong JH, Nam DH, Jeong BC, Park SH\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eSingle-cell RNA sequencing reveals the tumor microenvironment and facilitates strategic choices to circumvent treatment failure in a chemorefractory bladder cancer patient\u003c/strong\u003e. \u003cem\u003eGENOME MED \u003c/em\u003e2020, \u003cstrong\u003e12\u003c/strong\u003e(1):47.\u003c/li\u003e\n\u003cli\u003eZhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, Wang Y, Zhang Z, Yuan T, Ding X\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eSingle-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma\u003c/strong\u003e. \u003cem\u003eNAT COMMUN \u003c/em\u003e2020, \u003cstrong\u003e11\u003c/strong\u003e(1):6322.\u003c/li\u003e\n\u003cli\u003eWang Z, Jensen MA, Zenklusen JC: \u003cstrong\u003eA Practical Guide to The Cancer Genome Atlas (TCGA)\u003c/strong\u003e. \u003cem\u003eMethods Mol Biol \u003c/em\u003e2016, \u003cstrong\u003e1418\u003c/strong\u003e:111-141.\u003c/li\u003e\n\u003cli\u003eRitchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: \u003cstrong\u003elimma powers differential expression analyses for RNA-sequencing and microarray studies\u003c/strong\u003e. \u003cem\u003eNUCLEIC ACIDS RES \u003c/em\u003e2015, \u003cstrong\u003e43\u003c/strong\u003e(7):e47.\u003c/li\u003e\n\u003cli\u003eBruford EA, Antonescu CR, Carroll AJ, Chinnaiyan A, Cree IA, Cross NCP, Dalgleish R, Gale RP, Harrison CJ, Hastings RJ\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eHUGO Gene Nomenclature Committee (HGNC) recommendations for the designation of gene fusions\u003c/strong\u003e. \u003cem\u003eLEUKEMIA \u003c/em\u003e2021, \u003cstrong\u003e35\u003c/strong\u003e(11):3040-3043.\u003c/li\u003e\n\u003cli\u003eButler A, Hoffman P, Smibert P, Papalexi E, Satija R: \u003cstrong\u003eIntegrating single-cell transcriptomic data across different conditions, technologies, and species\u003c/strong\u003e. \u003cem\u003eNAT BIOTECHNOL \u003c/em\u003e2018, \u003cstrong\u003e36\u003c/strong\u003e(5):411-420.\u003c/li\u003e\n\u003cli\u003eMcGinnis CS, Murrow LM, Gartner ZJ: \u003cstrong\u003eDoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors\u003c/strong\u003e. \u003cem\u003eCell Syst \u003c/em\u003e2019, \u003cstrong\u003e8\u003c/strong\u003e(4):329-337 e324.\u003c/li\u003e\n\u003cli\u003eKorsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, Brenner M, Loh PR, Raychaudhuri S: \u003cstrong\u003eFast, sensitive and accurate integration of single-cell data with Harmony\u003c/strong\u003e. \u003cem\u003eNAT METHODS \u003c/em\u003e2019, \u003cstrong\u003e16\u003c/strong\u003e(12):1289-1296.\u003c/li\u003e\n\u003cli\u003eKim S, Kang D, Huo Z, Park Y, Tseng GC: \u003cstrong\u003eMeta-analytic principal component analysis in integrative omics application\u003c/strong\u003e. \u003cem\u003eBIOINFORMATICS \u003c/em\u003e2018, \u003cstrong\u003e34\u003c/strong\u003e(8):1321-1328.\u003c/li\u003e\n\u003cli\u003eQing X, Xu W, Liu S, Chen Z, Ye C, Zhang Y: \u003cstrong\u003eMolecular Characteristics, Clinical Significance, and Cancer Immune Interactions of Angiogenesis-Associated Genes in Gastric Cancer\u003c/strong\u003e. \u003cem\u003eFRONT IMMUNOL \u003c/em\u003e2022, \u003cstrong\u003e13\u003c/strong\u003e:843077.\u003c/li\u003e\n\u003cli\u003eEfremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R: \u003cstrong\u003eCellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes\u003c/strong\u003e. \u003cem\u003eNAT PROTOC \u003c/em\u003e2020, \u003cstrong\u003e15\u003c/strong\u003e(4):1484-1506.\u003c/li\u003e\n\u003cli\u003eJin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q: \u003cstrong\u003eInference and analysis of cell-cell communication using CellChat\u003c/strong\u003e. \u003cem\u003eNAT COMMUN \u003c/em\u003e2021, \u003cstrong\u003e12\u003c/strong\u003e(1):1088.\u003c/li\u003e\n\u003cli\u003eHanzelmann S, Castelo R, Guinney J: \u003cstrong\u003eGSVA: gene set variation analysis for microarray and RNA-seq data\u003c/strong\u003e. \u003cem\u003eBMC BIOINFORMATICS \u003c/em\u003e2013, \u003cstrong\u003e14\u003c/strong\u003e:7.\u003c/li\u003e\n\u003cli\u003eSzklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eSTRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets\u003c/strong\u003e. \u003cem\u003eNucleic Acids Res \u003c/em\u003e2019, \u003cstrong\u003e47\u003c/strong\u003e(D1):D607-D613.\u003c/li\u003e\n\u003cli\u003eSteen CB, Liu CL, Alizadeh AA, Newman AM: \u003cstrong\u003eProfiling Cell Type Abundance and Expression in Bulk Tissues with CIBERSORTx\u003c/strong\u003e. \u003cem\u003eMethods Mol Biol \u003c/em\u003e2020, \u003cstrong\u003e2117\u003c/strong\u003e:135-157.\u003c/li\u003e\n\u003cli\u003eLangfelder P, Horvath S: \u003cstrong\u003eWGCNA: an R package for weighted correlation network analysis\u003c/strong\u003e. \u003cem\u003eBMC BIOINFORMATICS \u003c/em\u003e2008, \u003cstrong\u003e9\u003c/strong\u003e:559.\u003c/li\u003e\n\u003cli\u003eYu G: \u003cstrong\u003eGene Ontology Semantic Similarity Analysis Using GOSemSim\u003c/strong\u003e. \u003cem\u003eMethods Mol Biol \u003c/em\u003e2020, \u003cstrong\u003e2117\u003c/strong\u003e:207-215.\u003c/li\u003e\n\u003cli\u003eYu G, Wang LG, Han Y, He QY: \u003cstrong\u003eclusterProfiler: an R package for comparing biological themes among gene clusters\u003c/strong\u003e. \u003cem\u003eOMICS \u003c/em\u003e2012, \u003cstrong\u003e16\u003c/strong\u003e(5):284-287.\u003c/li\u003e\n\u003cli\u003eWilkerson MD, Hayes DN: \u003cstrong\u003eConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking\u003c/strong\u003e. \u003cem\u003eBIOINFORMATICS \u003c/em\u003e2010, \u003cstrong\u003e26\u003c/strong\u003e(12):1572-1573.\u003c/li\u003e\n\u003cli\u003eZhang C, He H, Hu X, Liu A, Huang D, Xu Y, Chen L, Xu D: \u003cstrong\u003eDevelopment and validation of a metastasis-associated prognostic signature based on single-cell RNA-seq in clear cell renal cell carcinoma\u003c/strong\u003e. \u003cem\u003eAging (Albany NY) \u003c/em\u003e2019, \u003cstrong\u003e11\u003c/strong\u003e(22):10183-10202.\u003c/li\u003e\n\u003cli\u003eSiegel RL, Miller KD, Wagle NS, Jemal A: \u003cstrong\u003eCancer statistics, 2023\u003c/strong\u003e. \u003cem\u003eCA Cancer J Clin \u003c/em\u003e2023, \u003cstrong\u003e73\u003c/strong\u003e(1):17-48.\u003c/li\u003e\n\u003cli\u003eThanindratarn P, Dean DC, Nelson SD, Hornicek FJ, Duan Z: \u003cstrong\u003eAdvances in immune checkpoint inhibitors for bone sarcoma therapy\u003c/strong\u003e. \u003cem\u003eJ BONE ONCOL \u003c/em\u003e2019, \u003cstrong\u003e15\u003c/strong\u003e:100221.\u003c/li\u003e\n\u003cli\u003eLi YS, Liu Q, Tian J, He HB, Luo W: \u003cstrong\u003eAngiogenesis Process in Osteosarcoma: An Updated Perspective of Pathophysiology and Therapeutics\u003c/strong\u003e. \u003cem\u003eAM J MED SCI \u003c/em\u003e2019, \u003cstrong\u003e357\u003c/strong\u003e(4):280-288.\u003c/li\u003e\n\u003cli\u003eOgiwara Y, Nakagawa M, Nakatani F, Uemura Y, Zhang R, Kudo-Saito C: \u003cstrong\u003eBlocking FSTL1 boosts NK immunity in treatment of osteosarcoma\u003c/strong\u003e. \u003cem\u003eCANCER LETT \u003c/em\u003e2022, \u003cstrong\u003e537\u003c/strong\u003e:215690.\u003c/li\u003e\n\u003cli\u003eWei QF, Yao JS, Yang YT: \u003cstrong\u003eMicroRNA-1247 inhibits the viability and metastasis of osteosarcoma cells via targeting NRP1 and mediating Wnt/beta-catenin pathway\u003c/strong\u003e. \u003cem\u003eEur Rev Med Pharmacol Sci \u003c/em\u003e2019, \u003cstrong\u003e23\u003c/strong\u003e(17):7266-7274.\u003c/li\u003e\n\u003cli\u003eLiu MX, Liu QY, Liu Y, Cheng ZM, Liu L, Zhang L, Sun DH: \u003cstrong\u003eInterleukin-35 suppresses antitumor activity of circulating CD8(+) T cells in osteosarcoma patients\u003c/strong\u003e. \u003cem\u003eCONNECT TISSUE RES \u003c/em\u003e2019, \u003cstrong\u003e60\u003c/strong\u003e(4):367-375.\u003c/li\u003e\n\u003cli\u003eKim JA, Berlow NE, Lathara M, Bharathy N, Martin LR, Purohit R, Cleary MM, Liu Q, Michalek JE, Srinivasa G\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eSensitization of osteosarcoma to irradiation by targeting nuclear FGFR1\u003c/strong\u003e. \u003cem\u003eBiochem Biophys Res Commun \u003c/em\u003e2022, \u003cstrong\u003e621\u003c/strong\u003e:101-108.\u003c/li\u003e\n\u003cli\u003eXu W, Li Z, Zhu X, Xu R, Xu Y: \u003cstrong\u003emiR-29 Family Inhibits Resistance to Methotrexate and Promotes Cell Apoptosis by Targeting COL3A1 and MCL1 in Osteosarcoma\u003c/strong\u003e. \u003cem\u003eMed Sci Monit \u003c/em\u003e2018, \u003cstrong\u003e24\u003c/strong\u003e:8812-8821.\u003c/li\u003e\n\u003cli\u003eYang L, Long Y, Xiao S: \u003cstrong\u003eOsteosarcoma-Associated Immune Genes as Potential Immunotherapy and Prognosis Biomarkers\u003c/strong\u003e. \u003cem\u003eBIOCHEM GENET \u003c/em\u003e2024, \u003cstrong\u003e62\u003c/strong\u003e(2):798-813.\u003c/li\u003e\n\u003cli\u003eQiao Z, Li J, Kou H, Chen X, Bao D, Shang G, Chen S, Ji Y, Cheng T, Wang Y\u003cem\u003e et al\u003c/em\u003e: \u003cstrong\u003eHsa-miR-557 Inhibits Osteosarcoma Growth Through Targeting KRAS\u003c/strong\u003e. \u003cem\u003eFront Genet \u003c/em\u003e2021, \u003cstrong\u003e12\u003c/strong\u003e:789823.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Osteosarcoma, Single-cell sequencing, Cellular heterogeneity, Gene expression, Intercellular communication, Pathway analysis","lastPublishedDoi":"10.21203/rs.3.rs-5305987/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-5305987/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003e\u003cstrong\u003eBackground:\u003c/strong\u003e Osteosarcoma is an aggressive malignancy of bone that poses significant treatment challenges and has been a focal point of extensive research due to its complex pathogenesis. Despite advances in traditional therapeutic approaches, the intricate genetic and cellular landscape of osteosarcoma remains inadequately understood, emphasizing the need for innovative research methodologies to unravel its underlying mechanisms.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eObjective:\u003c/strong\u003e This study aims to leverage the power of single-cell transcriptome sequencing technology to elucidate the cellular heterogeneity, gene expression patterns, intercellular communication networks, and critical genetic pathways implicated in osteosarcoma. By doing so, we intend to contribute valuable insights into the biogenesis of this malignancy, which may ultimately inform precision treatment strategies.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eMethods:\u003c/strong\u003e Utilizing single-cell sequencing, we conducted a comprehensive analysis of osteosarcoma samples to identify diverse cellular subpopulations within the tumor microenvironment. Our focus on gene expression profiles revealed significant differences across these subpopulations. Moreover, we employed bioinformatics approaches to explore the intercellular communication networks and identified key ligand-receptor pairings, substantiating the role of angiogenesis-related genes prominently expressed in osteoblasts and their proliferative counterparts.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eResults:\u003c/strong\u003e Our findings underscore the critical involvement of angiogenesis in osteosarcoma pathogenesis, with notable pathway activity variations among distinct cellular subpopulations. Additionally, protein interaction network mapping has unveiled significant discrepancies in pathway activities and highlighted the potential functional roles of key genes involved in tumor progression.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConclusion:\u003c/strong\u003e This study offers a comprehensive exploration of the biological characteristics of osteosarcoma through single-cell sequencing technology, thereby establishing a robust theoretical foundation that may facilitate the development of targeted and effective therapeutic strategies. However, it is essential to recognize that these findings are preliminary, necessitating further validation through expanded sample sizes and integration of multi-omics data. Future research will delve deeper into the mechanisms of the identified key pathways and genes, with the aspiration of enhancing the prognostic outcomes and quality of life for patients with osteosarcoma.\u003c/p\u003e","manuscriptTitle":"Single-cell RNA landscape of the intratumoral heterogeneity and expression of angiogenesis-related genes in osteosarcoma","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-12-09 08:40:31","doi":"10.21203/rs.3.rs-5305987/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"a4b77a88-3e8b-404d-a50f-5eb5bfd67fcc","owner":[],"postedDate":"December 9th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-12-09T08:40:34+00:00","versionOfRecord":[],"versionCreatedAt":"2024-12-09 08:40:31","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-5305987","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-5305987","identity":"rs-5305987","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}