基于单细胞测序的骨肉瘤异质性解析及预后风险预测模型的构建

Heterogeneity analysis of osteosarcoma based on single-cell sequencing and construction of a prognostic prediction model

摘要

目的 基于骨肉瘤单细胞测序数据集鉴定不同细胞亚群并分析其功能，筛选各细胞群中核心基因构建风险评估模型。 方法 基于GEO数据库中骨肉瘤单细胞测序数据集GSE152048，根据已知细胞标志物基因的表达对各细胞簇进行注释，识别其中肿瘤细胞进行亚群分析，采用GO和KEGG分析各亚群的主要细胞功能，鉴定不同亚群间表达差异最显著的100个基因，进一步分析基因功能并构建预后风险预测模型。 结果 基于GEO骨肉瘤单细胞测序数据集，共筛选得到9组肿瘤细胞亚群，各亚群在细胞增殖、有丝分裂、免疫及代谢调控等方面发挥不同功能。进一步根据各亚群中表达差异最大的100个基因，构建了由S100A11，CFH，RPL37A以及HIST1H4C四基因组成的模型，基于该模型可以将骨肉瘤病例分为高风险组和低风险组，其中，高风险组患者的总生存期显著短于低风险组，且患者肿瘤增殖能力显著强于低风险组（P<0.001）。 结论 本研究在单细胞水平上系统描绘了骨肉瘤细胞的异质性，根据不同细胞间基因表达差异，构建了一组预后风险评估模型，为骨肉瘤的精准分型、预后评估及潜在治疗靶点的开发提供了新的理论依据。

Abstract

Objective To identify distinct cellular clusters and analyze their functional characteristics based on osteosarcoma single-cell RNA sequencing（scRNA-seq） data， and construct a risk assessment model based on the core genes derived from these subpopulations. Methods Osteosarcoma scRNA-seq data GSE152048 from GEO database were firstly applied to annotate the cell types based on the expression of known marker genes， and the tumor cells were focused for subpopulation analysis. Afterwards， the primary biological functions of each subpopulation were investigated using Gene Ontology（GO） and Kyoto Encyclopedia of Genes and Genomes（KEGG） pathway enrichment analyses. Then， the top 100 most significantly differentially expressed genes（DEGs） across subpopulations were identified and analyzed to develop a prognostic risk prediction model. Results Based on GEO osteoarcoma scRNA-seq data， a total of nine tumor cell subpopulations were identified， and different subpopulations exerted distinct functions in cell proliferation， mitosis， immune regulation， and metabolic regulation. Moreover， a four-gene prognostic signature model was successfully constructed based on the most significant DEGs among these subpopulations. Based on the model， osteosarcoma patients were divided into high-risk group and low-risk group， and the survival of patients was significantly lower in high-risk group than in low risk group， and the proliferation ability of tumor was stronger in high-risk group（P<0.001）. Conclusion The study systematically delineates the heterogeneity of osteosarcoma at single-cell level， and establishes a prognostic risk assessment model based on intercellular gene expression differences， which provides new insights for precise subtyping， prognostic evaluation and identification of potential therapeutic targets in osteosarcoma.

Graphical abstract

关键词

骨肉瘤 / 异质性 / 单细胞测序 / 基因特征 / 风险评估

Key words

osteosarcoma / heterogeneity / single-cell sequencing / gene signature / risk assessment

骨肉瘤是最常见的原发性恶性骨肿瘤，好发于儿童和青少年^［1］。尽管随着新辅助化疗和手术治疗的进展，骨肉瘤患者5年生存率整体已提升至60%~70%^［2］，但对于发生转移或复发的患者，预后仍然极差，5年生存率不足30%^［3］。肿瘤异质性是导致治疗失败和预后不良的核心挑战之一，其不仅表现为不同患者之间基因组的差异性，而且体现在同一患者肿瘤组织内部不同区域肿瘤细胞间遗传学差异，以及转录组和细胞功能状态的多样性，这种异质性促进了肿瘤的进化、适应性生存以及对化疗的耐受^［4，5］。

长期以来，对骨肉瘤的研究主要依赖于批量转录组测序（Bulk-RNA）技术，该方法将组织样本中的所有细胞进行混合检测，所得到的结果是各种细胞类型的平均信号^［6］。这种方法虽然能够识别与疾病整体相关的分子特征，但无法解析肿瘤微环境中不同细胞组分（如恶性细胞、免疫细胞、基质细胞）各自的信号，更难以揭示肿瘤细胞本身不同亚群的存在及其功能特性。单细胞RNA测序（scRNA-seq）技术的出现，为在单细胞分辨率下系统性解析肿瘤异质性提供了革命性的工具^［7，8］，通过scRNA-seq可以精确识别肿瘤细胞，并对这些恶性细胞进行更精细的亚群划分，从而深入探究驱动肿瘤恶性进展的关键细胞状态^［9］。

本研究利用scRNA-seq技术深入探索骨肉瘤肿瘤细胞的内在异质性。我们首先从公共数据库中获取骨肉瘤的scRNA-seq数据，并通过细胞注释识别肿瘤细胞群体。在此基础上，对肿瘤细胞进一步进行亚群分析，并深入探索每个亚群独特的生物学功能。进而，通过识别不同亚群间表达差异最显著的基因，筛选出驱动骨肉瘤异质性和恶性进展的核心基因集合。最后，基于这些关键基因构建一个稳健的预后风险预测模型，将骨肉瘤患者划分为具有显著临床预后差异的高、低风险组，并初步探索两组患者在病理学特征上的差异，从而为骨肉瘤的精准医疗提供新的分子分型依据和潜在的干预靶点。

1 材料与方法

1.1 GEO骨肉瘤单细胞测序数据的获取与处理

本研究首先从GEO公共数据库筛选骨肉瘤单细胞测序数据集，筛选标准：①基于人体组织进行测序，而非犬类等其他物种；②测序标本为实体组织，而非细胞系；③测序标本数量大于3例。

最终我们筛选得到数据集：GSE152048，该数据集包含了4例原发骨肉瘤标本及2例骨肉瘤肺转移标本的单细胞转录组数据。下载原始测序数据后，使用R语言Seurat软件包进行数据质控与标准化，剔除基因数小于500或大于2 500、线粒体基因比例超过5%的低质量细胞，以排除死细胞或受损细胞的影响。

1.2 细胞注释与肿瘤细胞鉴定

将单细胞数据经过质控和标准化后，对所有细胞进行主成分分析（PCA）和基于图论的聚类分析。随后，利用典型细胞标志物对细胞群进行注释：①T细胞：CD3D，CD3E，CD8A，CD4；②B细胞：CD79A，CD20；③髓系细胞：CD68，CD14，CD163；④内皮细胞：CD31，VWF；⑤成纤维细胞：COL1A1，DCN；⑥破骨细胞：CTSK。

最后，在将所有免疫细胞、内皮细胞等非恶性细胞注释并排除后，基于成骨标志物（CD4、RUNX2）的高表达及增殖活性，将剩余细胞定义为肿瘤细胞群，并纳入后续分析。

1.3 肿瘤细胞亚群分析

在识别得到肿瘤细胞群后，使用PCA对肿瘤细胞群进行降维分析，采用FindClusters函数（分辨率设置为0.8）进行无监督聚类，将肿瘤细胞划分为多个亚群，并使用UMAP进行可视化。利用FindAllMarkers函数识别每个亚群相较于其他群的特异性表达基因（log₂FC>0.5，调整后P<0.05），对这些特异性基因进行基因本体（GO）、京都基因与基因组百科全书（KEGG）通路、单样本基因富集分析（ssGSEA）以及基因集模块评分（AddModuleScore）功能富集分析以揭示各亚群的核心生物学功能。

1.4 肿瘤细胞亚群拟时序分析及关键基因的鉴定

在分析各肿瘤细胞亚群的初步功能后，使用Monocle3软件包对各细胞群进行拟时序分析，首先探索不同细胞亚群之间的前后演化顺序；然后筛选在细胞群拟时序演化过程中表达差异变化最大的前100个基因，并且对基因进行GO/KEGG功能富集分析，以探索其主要参与的生物学过程与信号通路。

1.5 预后风险预测模型的构建与验证

筛选得到骨肉瘤细胞亚群之间表达差异最大的前100个候选基因后，基于TARGET骨肉瘤数据集（一个专用于儿童肿瘤分析的数据集，包含患者的生存随访等预后信息以及肿瘤RNA测序及基因芯片数据），利用LASSO回归分析，构建风险评估模型。LASSO回归分析的参数设置为：设定standar⁃dize=TRUE对基因表达数据标准化，将收敛阈值设为1×10^-7，最大迭代次数maxit设为100 000，以保障高维基因矩阵拟合收敛；搭配cv.glmnet开展十折交叉验证，以C-index作为模型评价指标。

风险评分（risk score）的计算公式为：

Risk Score = Coef*Exp*Gene1+ Coef*Exp*Gene2+… Coef*Exp*GeneN

其中，Exp是基因的表达量，Coef是该基因在多变量Cox回归中的系数，Gene代表基因名称。根据风险评分的中位数，将患者分为高风险组和低风险组。使用Kaplan-Meier生存分析和log-rank检验比较两组之间的生存差异，通过受试者工作特征（ROC）曲线评估模型的预测准确性。

1.6 高、低风险组之间基因功能分析

在构建得到基因风险预测模型后，基于该模型将骨肉瘤患者标本分为高、低风险组，首先，对比两组标本之间主要的差异表达基因，制作差异基因火山图；然后分别利用GO/KEGG以及ssGSEA进行差异功能富集分析，阐释两组患者之间生物学功能的主要差异点。

2 结果

2.1 骨肉瘤细胞单细胞图谱与异质性亚群鉴定

基于GSE152048中6例骨肉瘤单细胞测序数据，经过细胞注释后，共得到8个细胞群，包括：肿瘤细胞、破骨细胞、T淋巴细胞、B淋巴细胞、内皮细胞、成纤维细胞、巨噬细胞及间充质细胞。筛选其中肿瘤细胞群进一步分析，通过无监督聚类分析将这些肿瘤细胞清晰地划分为9个细胞亚群（见图1A、B）。

观察细胞群在肿瘤样本中的占比，可以发现，9个细胞亚群在骨肉瘤原发及转移灶中同时存在，但是相对占比存在差异，例如在原发灶中cancer cell-0群占比明显高于转移灶，而在转移灶中cancer cell-7及cancer cell-8群占比明显高于原发灶（见图1C），提示在两种类型标本中肿瘤细胞亚群存在功能差异性。

2.2 肿瘤细胞亚群特征基因及生物学功能分析

针对9个肿瘤细胞亚群进一步分析，发现：各细胞亚群具有特征性表达基因，而且生物学功能存在差异性，其中：cancer cell-1群、cancer cell-2群及cancer cell-6群细胞，特征性高表达细胞周期相关基因（如MKI67，TOP2A，CENPE，KIF4A），生物学功能富集于“细胞周期”“DNA复制”“有丝分裂”“姐妹染色体分离”等通路，提示其为高增殖性肿瘤细胞（见图2）。而cancer cell-3群、cancer cell-4群特征性表达GSTA1、CRYAB等基因，而生物学功能富集于代谢相关通路，例如：“氧化磷酸化”“脂肪酸代谢”等，提示这些细胞属于高代谢亚群（见图2）。cancer cell-5群特征性表达CCL3、CCL3L1等细胞因子，而生物学功能主要富集于“炎性调控”“细胞因子分泌”“白细胞激活”等方面，提示这些细胞与免疫调控相关（见图2）。cancer cell-7群特征性表达细胞外基质（ECM）相关基因（如FGF7、CA2、KRBOX1），富集于“细胞外基质”“ECM-受体相互作用”等通路，表现出间质样表型（见图2）。此外，cancer cell-0群及cancer cell-8群高表达应激反应相关基因（如NEAT1、GAGE1、GAGE12F），富集于“蛋白质折叠”“MAPK信号通路”，提示其处于应激状态（见图2）。

2.3 肿瘤细胞亚群拟时序分析及关键差异基因鉴定

针对9个肿瘤细胞亚群进一步拟时序分析，结果显示：不同细胞亚群位于细胞演化的不同阶段，提示其生物学功能的差异性。例如，即使cancer cell-1群、cancer cell-2群及cancer cell-6群细胞均特征性高表达细胞周期及有丝分裂相关基因，但是在拟时序分析中，位于细胞演化的不同阶段，提示其细胞功能存在差异性（见图3）。

2.4 肿瘤细胞亚群间差异基因的生物学功能分析

进一步分析在肿瘤细胞拟时序演化过程中表达差异最明显的前100个基因，包括：随着拟时序进展，表达逐渐升高，或者逐渐降低，或者表达在中间段最高的基因，其表达变化趋势见图4。

对表达差异最明显前100个基因的功能分析显示，它们显著富集于与肿瘤进展密切相关的通路，包括“细胞氧化磷酸化”“细胞外基质重塑”“肿瘤免疫调控”及“MAPK”与“PI3K-Akt”等信号通路，提示这些生物学功能与骨肉瘤的异质性表型相关（见图5）。

2.5 基于肿瘤细胞亚群差异基因构建预后风险预测模型

利用TARGET数据库中89例骨肉瘤患者的转录组及临床数据，对上述100个基因进行单因素Cox生存分析，然后利用LASSO回归构建了一个包含S100A11、CFH、RPL37A以及HIST1H4C基因的预后风险预测模型，风险预测模型评分公式如下：

Risk Score =（-0.642）*S100A11+（-0.392）*CFH+0.213 9*RPL37A+（-0.543 6）*HIST1H4C

根据4个基因在骨肉瘤中的相对表达以及基因模型风险评分中位数，患者被分为高风险组和低风险组（见图6A-C）。Kaplan-Meier生存分析显示，高风险组患者的总生存期显著短于低风险组（P<0.001，见图6D）。ROC曲线显示，该模型1年、3年和5年生存率的曲线下面积（AUC）分别为0.91，0.87和0.89（见图6E），表明其具有良好的预后预测性能。

2.6 纳入模型的4个基因的单细胞表达定位

为了进一步分析LASSO分析所得的4个基因在骨肉瘤细胞中的表达定位，我们分析了其在9个肿瘤细胞亚群中的相对表达，结果显示，4个基因在肿瘤细胞亚群中表达存在一定差异，其中S100A11、RPL37A基因主要表达于cancer cell-0、cancer cell-2、cancer cell-4群，CFH基因在cancer cell-7、cancer cell-8群中最高，而HIST1H4C基因在cancer cell-1、cancer cell-2、cancer cell-6群中表达高于其他群（见图7）；4个基因分布于不同肿瘤细胞亚群，对提示肿瘤生物学特征及预后信息具有协同作用。

2.7 骨肉瘤高、低风险组间细胞功能差异性分析

高、低风险组患者中差异表达基因主要聚焦在“染色体有丝分裂”“细胞核分裂”“DNA复制”等细胞周期相关功能（见图8），提示两组患者最显著的生物学差异在于细胞增殖能力，这与临床上对肿瘤的认知一致，提示风险预测模型的可靠性。

3 讨论

骨肉瘤是儿童及青少年最常见的原发性恶性骨肿瘤，其主要发生于长骨干骺端，表现出特征性的溶骨性破坏和成骨性增生，具有高度侵袭性，极易突破骨皮质侵入周围软组织，初诊时即有15%~20%的患者已发生远处转移^{［10，11］}，然而，近年来骨肉瘤的临床治疗无明显进展，手术及新辅助化疗仍是主要手段，患者对化疗反应不一，内在性耐药和获得性耐药普遍存在，这与骨肉瘤的高度异质性直接相关^［1，12］。由于缺乏高频、可靶向的驱动基因突变，骨肉瘤的靶向治疗至今未取得实质性突破^{［13，14］}。尽管靶向治疗在黑色素瘤、乳腺癌、肺癌等恶性肿瘤中已取得显著进展^［15-17］，但骨肉瘤因基因组高度不稳定、驱动基因突变频率低且异质性极强，其靶向药物的研发仍面临巨大挑战^{［18，19］}。解析骨肉瘤的异质性、识别不同亚群肿瘤细胞的独特功能，对探索新型基因标志物具有重要意义^{［20，21］}。本研究利用单细胞测序技术，系统鉴定了骨肉瘤9个功能各异的肿瘤细胞亚群，并构建了包含S100A11、CFH、RPL37A和HIST1H4C四个基因的预后风险模型，为骨肉瘤的精准分型和潜在靶点发现提供了新的依据。

首先，本研究通过单细胞测序分析识别出9个肿瘤细胞亚群，各亚群呈现显著的功能分化：高增殖亚群（cancer cell-1、cancer cell-2、cancer cell-6）富集于细胞周期和有丝分裂通路，高代谢亚群（cancer cell-3、cancer cell-4）富集于氧化磷酸化和脂肪酸代谢通路，免疫调控亚群（cancer cell-5）高表达CCL3等趋化因子，间质样亚群（cancer cell-7）高表达细胞外基质相关基因，应激响应亚群（cancer cell-0、cancer cell-8）富集于蛋白质折叠和MAPK信号通路。已有研究通过单细胞测序在骨肉瘤临床样本中同样鉴定出多个功能各异的肿瘤细胞亚群，发现这些亚群具有截然不同的基因表达特征，涵盖了富含干性特征的细胞、处于不同细胞周期阶段的细胞以及表现出成骨、成软骨或成纤维分化的细胞等^［22-24］。本研究结果与上述文献相互印证，共同表明即使在同一肿瘤内，也存在着从增殖性祖细胞到相对分化状态的细胞连续谱系。然而，本研究揭示了各亚群在代谢调控和免疫调节维度的功能差异，拓展了既往主要聚焦于分化状态的认识。综合文献和本研究结果可以得出：骨肉瘤的恶性进展并非由单一细胞群体驱动，而是多个功能各异的亚群协同作用的结果，其中增殖亚群驱动肿瘤快速生长，间质样亚群通过细胞外基质重塑促进侵袭和转移，代谢及应激相关亚群则赋予肿瘤在恶劣微环境下的生存优势，这可能与化疗抵抗和不良预后直接相关。

在此基础上，我们基于肿瘤细胞亚群间差异表达基因，进一步构建了包含S100A11、CFH、RPL37A和HIST1H4C的预后风险预测模型，高风险组患者总生存期显著缩短，模型预测1，3，5年生存率的AUC分别达0.91，0.87和0.89。该模型的4个基因并非随机组合，而是分别从不同维度反映了肿瘤的恶性特征。已有研究报道，S100A11可通过肿瘤细胞来源的细胞外囊泡介导肺转移前微环境形成，其高表达与骨肉瘤患者的肺转移和不良预后显著相关^［25］；CFH已被鉴定为骨肉瘤中与预后相关的关键基因之一，并与免疫微环境密切相关^［26］；RPL37A在骨肉瘤细胞中高表达，敲低其表达可显著抑制骨肉瘤细胞的增殖、迁移和侵袭能力^［27］；HIST1H4C作为组蛋白家族成员，负责核小体结构的维持，在多种肿瘤中表达异常与不良预后相关^［28］。本研究的单细胞定位发现，这4个基因具有明确的亚群分布特征：S100A11和RPL37A主要表达于cancer cell-0、cancer cell-2、cancer cell-4群（应激、增殖和代谢亚群），CFH高表达于cancer cell-7、cancer cell-8群（间质样和应激亚群），HIST1H4C高表达于cancer cell-1、cancer cell-2、cancer cell-6群（高增殖亚群）。将文献中4个基因的已知功能与本研究中它们在不同亚群的特异性表达相结合可以推断，该模型实际上整合了多个功能亚群的“集体信号”：增殖信号（HIST1H4C、RPL37A）、转移与免疫逃逸信号（S100A11、CFH），而非单一亚群的孤立特征，这解释了为何该模型能够更全面地反映肿瘤的整体恶性潜能，并具有良好的预后预测性能。

为验证该模型的预测效能，基于四基因模型将患者分为高、低风险组，本研究发现两组之间的差异表达基因主要富集于“染色体有丝分裂”“细胞核分裂”“DNA复制”等细胞周期相关功能，高风险组患者肿瘤增殖能力显著强于低风险组。文献已明确指出骨肉瘤的侵袭性和不良预后与其快速增殖能力密切相关^{［10，11］}，而高风险组与低风险组之间最显著的生物学差异同样聚焦于增殖能力。这一发现与临床上预后较差肿瘤通常生长更快的观察高度一致。更重要的是，该结果与本研究所识别的肿瘤细胞亚群功能特征相呼应：高增殖亚群（cancer cell-1、cancer cell-2、cancer cell-6）所富集的生物学通路正是“细胞周期”“DNA复制”“有丝分裂”等，提示高风险组患者可能富集了更多高增殖亚群的细胞。这种从“单细胞亚群功能”到“临床预后分组”的关联链条，为理解骨肉瘤预后差异的细胞根源提供了直接证据，也侧面印证了该预后模型的可靠性。

然而，本研究存在一定局限性。首先，数据来源于公共数据库（GSE152048），标本量有限（仅6例），且缺乏患者详细的治疗信息，限制了预后模型与不同化疗方案反应之间关系的深入分析。此外，目前骨肉瘤的单细胞测序数据本身就较为稀缺，已有研究指出，骨组织质地坚硬、含大量钙化骨基质，导致高质量单细胞悬液制备困难，且测序成本较高、数据分析复杂，这使该领域的单细胞技术分析仍处于较早期阶段^［29］，这也在客观上限制了本研究利用单细胞层面数据进行更深入验证的可能。

综上所述，本研究在单细胞水平上系统描绘了骨肉瘤9个功能各异的肿瘤细胞亚群，揭示了骨肉瘤异质性的细胞和分子基础；在此基础上成功构建了由S100A11、CFH、RPL37A和HIST1H4C四基因组成的预后风险预测模型，该模型从增殖、转移、免疫逃逸和表观调控等多个维度综合反映了肿瘤的恶性潜能。这些发现为骨肉瘤的精准分型、预后评估及潜在治疗靶点的开发提供了新的理论依据。

参考文献

[1] Meltzer PS， Helman LJ. New horizons in the treatment of osteosarcoma［J］. N Engl J Med， 2021， 385（22）： 2066-2076.

[2] Robinson MJ， Davis EJ. Neoadjuvant chemotherapy for adults with osteogenic sarcoma［J］. Curr Treat Options Oncol， 2024， 25（11）： 1366-1373.

[3] Rubio⁃San⁃Simón A， Wilson W， Sironi G， et al. Prognostic factors in patients with relapsed high⁃grade osteosarcoma： a systematic review［J］. Br J Cancer， 2025， 133（7）： 1020-1028.

[4] Rickel K， Fang F， Tao J. Molecular genetics of osteosarcoma［J］. Bone， 2017， 102： 69-79.

[5] Baumhoer D. Pathogenesis and genetics of osteosarcoma： current concepts and developments［J］. Pathologe， 2018， 39（2）： 139-145.

[6] Sun Y， Zhang C， Fang Q， et al. Abnormal signal pathways and tumor heterogeneity in osteosarcoma［J］. J Transl Med， 2023， 21（1）： 99.

[7] Su M， Pan T， Chen QZ， et al. Data analysis guidelines for single⁃cell RNA⁃seq in biomedical studies and clinical applications［J］. Mil Med Res， 2022， 9（1）： 68.

[8] Huang D， Ma N， Li X， et al. Advances in single⁃cell RNA sequencing and its applications in cancer research［J］. J Hematol Oncol， 2023， 16（1）： 98.

[9] Truong DD， Weistuch C， Murgas KA， et al. Mapping the single⁃cell differentiation landscape of osteosarcoma［J］. Clin Cancer Res， 2024， 30（15）： 3259-3272.

[10] Beird HC， Bielack SS， Flanagan AM， et al. Osteosarcoma［J］. Nat Rev Dis Primers， 2022， 8（1）： 77.

[11] Shoaib Z， Fan TM， Irudayaraj JMK. Osteosarcoma mechanobiology and therapeutic targets［J］. Br J Pharmacol， 2022， 179（2）： 201-217.

[12] Gill J， Gorlick R. Advancing therapy for osteosarcoma［J］. Nat Rev Clin Oncol， 2021， 18（10）： 609-624.

[13] Corre I， Verrecchia F， Crenn V， et al. The osteosarcoma microenvironment： a complex but targetable ecosystem［J］. Cells， 2020， 9（4）： 976.

[14] Rathore R， Van Tine BA. Pathogenesis and current treatment of osteosarcoma： perspectives for future therapies［J］. J Clin Med， 2021， 10（6）： 1182.

[15] Boutros A， Croce E， Ferrari M， et al. The treatment of advanced melanoma： current approaches and new challenges［J］. Crit Rev Oncol Hematol， 2024， 196： 104276.

[16] Hong R， Xu B. Breast cancer： an up⁃to⁃date review and future perspectives［J］. Cancer Commun（Lond）， 2022， 42（10）： 913-936.

[17] Tan AC， Tan DSW. Targeted therapies for lung cancer patients with oncogenic driver molecular alterations［J］. J Clin Oncol， 2022， 40（6）： 611-625.

[18] Sayles LC， Breese MR， Koehne AL， et al. Genome⁃informed targeted therapy for osteosarcoma［J］. Cancer Discov， 2019， 9（1）： 46-63.

[19] Shaikh AB， Li F， Li M， et al. Present advances and future perspectives of molecular targeted therapy for osteosarcoma［J］. Int J Mol Sci， 2016， 17（4）： 506.

[20] Yu S， Yao X. Advances on immunotherapy for osteosarcoma［J］. Mol Cancer， 2024， 23（1）： 192.

[21] Wu C， Tan J， Shen H， et al. Exploring the relationship between metabolism and immune microenvironment in osteosarcoma based on metabolic pathways［J］. J Biomed Sci， 2024， 31（1）： 4.

[22] Li J， Bai Y， Zhang H， et al. Single⁃cell RNA sequencing reveals the communications between tumor microenvironment components and tumor metastasis in osteosarcoma［J］. Front Immunol， 2024， 15： 1445555.

[23] Zhou Y， Yang D， Yang Q， et al. Single⁃cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma［J］. Nat Commun， 2020， 11（1）： 6322.

[24] Liu Y， Feng W， Dai Y， et al. Single⁃cell transcriptomics reveals the complexity of the tumor microenvironment of treatment⁃naive osteosarcoma［J］. Front Oncol， 2021， 11： 709210.

[25] Deng C， Chen H， Yang Y， et al. Extracellular⁃vesicle⁃packaged S100A11 from osteosarcoma cells mediates lung premetastatic niche formation by recruiting gMDSCs［J］.Cell Rep， 2024， 43（2）： 113751.

[26] Chu D， Huang R， Shi J， et al. NETs⁃related genes predict prognosis and are correlated with the immune microenvironment in osteosarcoma［J］. Front Oncol， 2025， 15： 1551074.

[27] Liu Z， Zhong Y， Meng S， et al. Identification of a seven⁃gene prognostic signature using the gene expression profile of osteosarcoma［J］. Ann Transl Med， 2022， 10（2）： 53.

[28] Han JA， Kim JY， Kim JI. Analysis of gene expression in cyclooxygenase⁃2⁃overexpressed human osteosarcoma cell lines［J］. Genomics Inform， 2014， 12（4）： 247-253.

[29] Wang T， Wang L， Zhang L， et al. Single⁃cell RNA sequencing in orthopedic research［J］. Bone Res， 2023 11（1）：10.

基金资助

国家自然科学基金项目(82572996)

山西省科技厅应用基础研究项目(202303021222333)

山西省科技厅应用基础研究项目(202403021211135)

山西省科技厅应用基础研究项目(202403021222467)