A single-cell and spatially resolved atlas of human osteosarcomas

Abstract

Osteosarcomas are intricate cellular ecosystems, where heterotypic interactions significantly influence disease progression and therapeutic outcomes. Despite their importance, a detailed understanding of their cellular composition and organizational structure remains elusive. In this study, we provide a comprehensive single-cell and spatially resolved transcriptomics analysis of human osteosarcomas. We construct a cellular meta-map to dissect spatial transcriptomic data, unveiling a detailed atlas of osteosarcoma compositional subgroups. We meticulously characterize the unique gene signatures and functional states of each subgroup and investigate the impact of chemotherapy on these cellular subpopulations. Additionally, our spatial transcriptomics analysis identifies a distinct spatial niche, located at the forefront of tumor necrotic zones, potentially associated with chemotherapy resistance. We also delve into the crosstalk between different cellular subgroups. This study furnishes a comprehensive transcriptional atlas of osteosarcoma's cellular architecture, enriching our comprehension of its complexity and laying the groundwork for more targeted therapeutic approaches.

To the editor

Osteosarcoma, a primary bone tumor mainly affecting children and young adults, has seen limited improvement in patient prognosis despite chemotherapy advancements in the 1980s [1]. Single-cell RNA sequencing (scRNA-seq) has shed light on the complex tumor microenvironment, but the rarity of osteosarcoma has limited comprehensive studies, with only 17 samples across current datasets. Despite ongoing meaningful in-depth analyses of public data [2], there remains a significant gap in our understanding of the spatial distribution and interactions among various cell types within tumors. To address this, we integrated our newly generated dataset of 10 osteosarcoma samples with existing datasets (GSE162454 [3] and GSE152048 [4]), forming the largest scRNA-seq cohort of 27 samples. We also pioneered the establishment and analysis of spatial transcriptomic data for osteosarcoma. This allowed us to map a detailed scRNA-seq atlas of osteosarcoma and establish spatial transcriptomic data, similar to other scRNA-seq atlases [5], providing insights across different patient profiles (Fig. 1A).

Fig. 1

A single-cell and spatially resolved atlas of human osteosarcomas. A Workflow of this study. B UMAP plot displaying 49,356 single cells from osteosarcoma tissues, annotated by major cell types in NCC sc-RNA seq data. (osteosarcoma cells (n = 10,467), myeloid cells (n = 13,966), mesenchymal stromal cells (MSCs; n = 3,240), endothelial cells (ECs; n = 2,084), osteoclast cells (OCs; n = 2,555), B cells (n = 462), plasma cells (n = 206), mast cells (n = 184), myocytes (n = 1,790) and TNK cells (n = 14,402)). C UMAP plot displaying 194,083 single cells from osteosarcoma tissues, annotated by major cell types in integrated sc-RNA seq data. (osteosarcoma cells (n = 87,256), myeloid cells (n = 51,894), MSCs (n = 9,867), ECs (n = 6,960), OCs (n = 12,910), B cells (n = 1,094), plasma cells (n = 874), mast cells (n = 592), myocytes (n = 2,050) and TNK cells (n = 21,306)). D Pie chart depicting the proportion of major cell types within the sample in NCC sc-RNA seq data. E Pie chart depicting the proportion of major cell types within the sample in integrated sc-RNA seq data. F Spearman correlation analysis showing relationships between proportions of major cell types; MSC and endothelium demonstrate the correlation in NCC sc-RNA seq data. G. Spearman correlation analysis showing relationships between proportions of major cell types; MSC and endothelium demonstrate the correlation in integrated sc-RNA seq data. H MRI images of the BS3 sample before and after three chemotherapy cycles, indicating no significant tumor reduction, suggestive of potential chemotherapy resistance. I H&E staining of a BS3 tissue slice reveals less than 50% necrosis, indicating potential chemotherapy resistance. Spatial transcriptomics show a low number of detected genes in the necrotic areas of BS3. J Cell2location deconvolution analysis estimates the proportions and spatial distribution of primary cell types (OS, MSC, endothelium, myeloid, osteoclasts, and TNK cells) in SP_BS3. K Analysis by Scanpy identifies six spatial niches within SP_BS3, with cluster_6 located at the tumor necrosis front, highly expressing genes such as COL4A1, suggesting a potential link to chemotherapy resistance

Analysis of new osteosarcoma sc-RNA seq dataset from National Cancer Center (NCC) and the largest osteosarcoma sc-RNA seq dataset

After conducting quality control and data processing, We analyzed the NCC osteosarcoma sc-RNA seq dataset, identifying 49,356 cells across 10 clusters (Fig. 1B; Fig. S1A, B). By integrating this with existing data, we created the largest osteosarcoma sc-RNA seq dataset (Fig. 1C; Fig. S1C–F). The heterogeneity in cell distribution, intercellular correlations, and changes in gene expression of major cell types before and after chemotherapy were all described, providing insights into the cellular response to treatment. (Fig. 1D–G and Fig. S1 G–K).

Analysis of the first osteosarcoma spatial transcriptomics data

Sample BS3, after three cycles of chemotherapy, showed minimal lesion shrinkage on MRI (Fig. 1H). Spatial RNA sequencing revealed characteristics of chemoresistance, with post-surgical histopathology (H&E staining) showing a necrosis rate below 50% and a low gene count in necrotic areas (Fig. 1I). We focused on analyzing these necrotic regions, exploring the spatial distribution of major cell types and gene expression to uncover potential links between gene expression and chemotherapy resistance in osteosarcoma (Fig. 1J-K).

Subcluster analysis of osteosarcoma cells from integrated sc-RNA seq data

We performed unbiased clustering of malignant osteosarcoma cells and identified a total of seven subclusters. (Fig. 2A and Fig. S2A–D). We identified the characteristics and functions of each subcluster based on their gene expression profiles and pathway enrichment analysis. Specifically, the OS_c3_TAGLN subcluster appears to highlight a fibroblast-like trajectory. The OS_c5_TOP2A subcluster exhibits a proliferative phenotype, while the OS_c1_JUN subcluster is involved in stress response mechanisms and is enriched in TNFα signaling through the NFKB, WNT/β-catenin, and TGF-β signaling pathways., among others. (Fig. S2E, F).

Fig. 2

Subcluster analysis of cells from integrated sc-RNA seq data. A Dimensionality reduction plot (Dimplot) for main osteosarcoma cell types: OS_C1_JUN, OS_C2_MT1X, OS_C3_TAGLN, OS_C4_SPP1, OS_C5_TOP2A, OS_C6_IBSP, and OS_C7_COL2A1. B Dimplot for five MSC subtypes: matrix CAF (mCAF), inflammatory CAF (iCAF), myogenesis CAF (myoCAF), vascular CAF (vCAF), and pericytes. C Cellchat analysis of cell communication between ten main cell types, with MSCs showing the strongest interaction with endothelium. D Cellchat analysis of cell communication between mCAF, myoCAF, iCAF, pericytes, and endothelium, with mCAF showing the strongest interaction with endothelium. E Key receptors and ligands involved in cell communication between mCAF, myoCAF, iCAF, pericytes, and endothelium. F Dimplot for nine TNK subtypes: three CD4 + subtypes, four CD8 + subtypes, and two NK cell subtypes. G Dimplot for nine myeloid subtypes: Macro_C1QC, Macro_LYVE1, Macro_SPP1, Macro_NLRP3, Mast_KIT, Mono_CD14, cDC1_CLEC9A, cDC2_CD1C, and Cdc3_LAMP3. H Dimplot of five main osteoclast (OC) subclusters with pseudotime analysis showing progression from OC_immature to OC_early response, and finally to OC_mature. I Monocle3 analysis identifies differentially expressed genes (DEGs) during pseudotime, with ACP5 and CTSK (known osteoclast markers) increasingly expressed from OC_immature to OC_mature

Subcluster analysis of mesenchymal stromal cells (MSCs) from integrated sc-RNA seq data

MSCs were divided into distinct clusters based on gene expression, revealing diverse origins and functions (Fig. 2B and Fig. S2G–L). Using Gene Set Variation Analysis (GSVA) and examining changes in cell abundance before and after chemotherapy, we identified the functional roles of these clusters (Fig. S2M, N). For instance, inflammatory fibroblasts (iCAF) are involved in antitumor activities, myofibroblasts (myoCAFs) play a role in wound healing, and vascular CAFs (vCAFs) are linked to tumor resistance. Cell communication analysis showed MSCs closely interact with endothelial cells, highlighting FN1 − (ITGA5 + ITGB1) as a key communication signal between pericytes and endothelial cells, as well as between matrix-associated fibroblasts (mCAFs) and endothelial cells, suggesting a potential synergistic interaction (Fig. 2C–E).

Subcluster analysis of TNK cells and myeloid cells from integrated sc-RNA seq data

We conducted subcluster analyses of TNK and myeloid cells from integrated scRNA-seq data. TNK cells were divided into nine subpopulations, including distinct CD4 + and CD8 + T cell types, and NK cells, highlighting roles in immune activation, suppression, and cytotoxicity (Fig. 2F, Fig. S3A–E). Myeloid cells formed nine clusters, with osteoclasts showing a maturation trajectory, crucial for osteosarcoma progression (Fig. 2G, H, Fig. S3F–M). Pseudotime analysis traced osteoclast evolution from immature to mature stages, with increasing expression of ACP5 and CTSK (Fig. 2H, I).

In summary, we provided a detailed description of the osteosarcoma cellular landscape at both single-cell and spatial transcriptomic levels, explored spatial niches, partially uncovered mechanisms of chemotherapy resistance, and identified potential therapeutic targets.