通过结合批量测序、单细胞和空间转录组学分析,以及基于机器学习的生存建模和实验验证,本研究提出PPARG是骨肉瘤中一个候选的与衰老相关的预后基因,并将其表达水平降低与TARGET-OS队列中较差的生存率以及血管和微环境特征相关联。
研究文章
* These authors contributed equally
通过结合批量测序、单细胞和空间转录组学分析,以及基于机器学习的生存建模和实验验证,本研究提出PPARG是骨肉瘤中一个候选的与衰老相关的预后基因,并将其表达水平降低与TARGET-OS队列中较差的生存率以及血管和微环境特征相关联。
骨肉瘤在转移性、复发性或治疗耐药性疾病中仍具有挑战性。本研究旨在鉴定与衰老相关的预后基因,并表征其空间背景。在GSE99671数据集中使用DESeq2进行配对差异表达分析,随后与CellAge衰老基因集取交集。从UCSC Xena获取TARGET-OS队列的转录组数据。采用单变量Cox回归、Kaplan-Meier分析、时间依赖性受试者工作特征分析、LASSO Cox回归、重复LASSO分析及随机生存森林模型评估候选基因,并将临床协变量纳入校正后的Cox模型中。通过功能富集分析、免疫微环境分析、单细胞转录组学、SP_BS3空间转录组学、GSE36001表达验证以及在骨肉瘤143B细胞和成骨细胞中进行qRT-PCR和Western blot验证,对候选基因进行表征。在GSE99671中,共鉴定出2,248个差异表达基因(校正后P < 0.05),其与866个CellAge基因的交集得到105个与衰老相关的差异表达基因。在TARGET-OS队列中,PPARG表达水平较低与较高的死亡风险相关(单变量HR = 0.603,95% CI = 0.454–0.802,P = 0.000494;校正后HR = 0.224,95% CI = 0.085–0.589,P = 0.00241)。将PPARG加入临床模型后,C指数从0.707提升至0.829。PPARG在GSE99671和GSE36001中均呈下调,qRT-PCR和Western blot结果证实,与成骨细胞相比,骨肉瘤143B细胞中PPARG的mRNA和蛋白表达水平均较低。单细胞分析显示PPARG主要定位于内皮细胞、周细胞、巨噬细胞/单核细胞及肿瘤相关基质细胞。空间转录组分析显示,PPARG表达水平与CellAge衰老评分、内皮相关评分及周细胞相关评分之间存在微弱但显著的正相关关系。这些结果提示PPARG是一个源自CellAge的候选预后生物标志物,与骨肉瘤不良生存预后及血管微环境特征相关,支持进一步探讨其在风险分层及衰老相关肿瘤微环境中的潜在意义。GSE36001仅提供了外部表达验证,未进行独立的生存验证。
骨肉瘤是儿童、青少年和年轻成人中最常见的原发性恶性骨肿瘤1。尽管联合手术与多药化疗已改善了局限性疾病的预后2,3,但转移性、复发性或治疗耐药性骨肉瘤患者的长期生存率仍然较差3,4。新兴证据表明,氧化应激诱导的表观遗传重塑可促进转移性适应和肿瘤进展,凸显了侵袭性癌症表型背后的复杂分子可塑性5。目前仍缺乏兼具临床可解读性和生物学信息价值的可靠生物标志物。因此,识别能够跨多个数据层面反映骨肉瘤异质性及预后风险的分子特征仍具有重要意义。
细胞衰老是一种由端粒功能障碍、DNA损伤、氧化应激、癌基因激活以及治疗压力所诱导的稳定的细胞周期阻滞程序6。衰老可抑制异常增殖;然而,衰老细胞也可通过炎症因子、趋化因子、生长因子以及细胞外基质重塑程序来重塑肿瘤微环境7,8。在骨肉瘤中,衰老相关基因可能同时反映肿瘤细胞内在的应激状态和非恶性微环境组分,但其预后意义及空间分布特征尚未得到系统评估。
PPARG 编码过氧化物酶体增殖物激活受体γ(peroxisome proliferator-activated receptor gamma),这是一种配体激活的核受体,参与脂质代谢、炎症调控、细胞分化以及免疫调节9。PPARG 在癌症中的作用具有环境依赖性10。在某些情况下,PPARG 与细胞分化和抗炎状态相关,而在其他情况下,它可能支持肿瘤或间质细胞的适应性程序。然而,PPARG 在骨肉瘤中的表达模式、预后价值及其细胞和空间定位尚未被充分阐明。
在本研究中,鉴定了GSE99671中差异表达的基因,并将其与CellAge衰老基因集取交集,得到105个与衰老相关的差异表达基因。随后结合TARGET-OS生存数据,采用多种机器学习和生存模型分析方法,并进行临床因素校正,筛选出PPARG为核心基因。进一步通过批量功能分析和免疫微环境分析、单细胞转录组学、空间转录组学,以及来自GSE36001的外部表达验证,在骨肉瘤细胞系143B和人成骨细胞中进行qRT-PCR和Western blot实验,对PPARG进行了系统表征。
本研究使用公共数据集和细胞系进行分析和验证,未涉及人类参与者或临床组织样本,因此无需伦理审批。
GSE99671 与 CellAge 数据交集的差异表达分析
GSE99671 的原始计数数据和样本分组信息来自 GEO11,12。共分析了 18 对骨肉瘤组织及其配对的正常组织样本。使用 DESeq2 进行配对差异表达分析,设计公式为 ~ pair_id + condition,其中 pair_id 用于校正个体配对效应,condition 用于比较肿瘤组织与正常组织13。仅保留至少在三个样本中读段计数不低于 10 的基因。差异表达定义为经校正的 P < 0.05;在可视化时采用更严格的阈值,即经校正的 P < 0.05 且 |log2FC| ≥ 1。在将基因符号转换为大写后,将差异表达基因与 866 个 CellAge 衰老相关基因取交集14。差异表达分析在 R(版本 4.3.2)中使用 DESeq2(版本 1.40.2)完成,校正后的 P 值采用 Benjamini-Hochberg 方法计算。
TARGET-OS 队列与预后模型构建
TARGET-OS 的转录组学和临床数据来自 UCSC Xena15。提取了 105 个衰老相关差异表达基因的表达值。共纳入 85 例具有完整生存时间、生存状态及候选基因表达数据的患者,其中发生死亡事件 27 例。采用单变量 Cox 回归16、Kaplan-Meier 生存分析和时间依赖性受试者工作特征(ROC)分析17对标准化后的基因表达值进行分析。通过 LASSO Cox 回归18、重复 LASSO 稳定性分析以及随机生存森林模型评估变量选择的稳定性与重要性19。整合枢纽评分(integrated hub score)和临床整合排序(clinical-integrated ranking)依据下文详述的明确二分类标准进行计算。所有生存分析均在 R 软件中进行,使用的软件包包括 survival(版本 3.5-7)、timeROC(版本 0.4)、glmnet(版本 4.1-8)和 randomForestSRC(版本 3.2.2)。针对 105 个候选基因的单变量 Cox 筛选,采用 Benjamini-Hochberg 方法进行假发现率(FDR)校正,FDR < 0.05 的基因被视为具有统计学显著性。
模型预处理、PPARG 分组及时间依赖性 ROC 分析
在 85 例患者中,共观察到 27 例死亡事件,排除了方差为零的基因,对缺失的候选基因表达值进行中位数填补,并对每个候选基因的表达值进行 z 分数标准化。在 Kaplan-Meier 分析中,以队列中位数对表达水平进行二分类:严格高于中位数的值归为高表达组,等于或低于中位数的值归为低表达组(PPARG:高表达组 n = 42,低表达组 n = 43)。Log-rank 检验采用双侧检验。时间依赖性 ROC 分析使用 timeROC 方法,事件原因设为 1,采用边际逆概率删失加权法,评估时间点为 365、1,095 和 1,825 天,iid 参数设为 FALSE。为确保较大的标志物值始终表示更高的风险,对于 Cox 回归系数为正的基因,直接使用标准化表达值;而对于具有负系数的保护性基因,其标准化表达值乘以 -1。
LASSO 与重复 LASSO
使用 glmnet 拟合 Cox LASSO 模型,参数设置为 family = "cox"、alpha = 1、预先对 z 分数进行标准化(因此 standardize = FALSE)、五折交叉验证、type.measure = "deviance",并设定随机种子为 123。主要的系数解采用 lambda.min。稳定性分析将相同的五折交叉验证重复 300 次;第 b 次重复使用的随机种子为 1000 + b(b = 1,...,300)。对于每个基因,选择频率是指在 lambda.min 处具有非零系数的重复次数所占的比例;同时记录在 lambda.1se 处具有非零系数的情况。
随机生存森林
使用 randomForestSRC(版本 3.2.2)对全部 105 个标准化候选基因拟合生存森林,设定随机种子为 123,ntree = 1,000,importance = TRUE,na.action = "na.impute"。对于生存数据,采用软件包默认设置:对数秩分割(log-rank splitting),mtry = 11(105 个预测变量平方根的上取整),最小终端节点大小 = 15,nsplit = 10 个随机分割点,无放回抽样且抽样比例为 0.632,以及反分割变量重要性(anti-split variable importance)。
综合枢纽评分
每个与CellAge相关的105个差异表达候选基因根据六项二元标准各获得一分:(1)差异表达与CellAge交集成员资格(所有候选基因均获得此项分数,因为经调整后) P < 0.05 才进行评分);(2)名义单变量 Cox P < 0.05;(3)Kaplan-Meier 对数秩检验 P < 0.05;(4)3年和5年时依时性AUC的平均值 ≥ 0.65;(5)在所有候选变量中,lambda.min重复LASSO选择频率位于第70百分位数或以上 > 0;以及(6)在所有候选变量中,随机生存森林重要性位于第70百分位数及以上 > 0. 所有标准均具有相等的单位权重,得到0到6之间的枢纽评分;评分≥4的基因进入临床调整阶段(共18个基因)。同时报告单变量Cox FDR值,以及FDR < 0.05 用于表示多重检验显著性,但预先指定的评分指标采用名义显著性水平 P < 0.05.
临床整合评分
在合并表达数据与临床记录后,校正分析共纳入40例具有完整协变量记录的患者,其中13例死亡。最终得分为初始中心性得分加上满足以下七项标准的每一项各1分:校正Cox P < 0.05、校正Cox P < 0.10、剔除根治性手术后的敏感性Cox P < 0.05、敏感性Cox P < 0.10、3年和5年平均AUC ≥ 0.65、临床+基因模型相较于仅临床模型的似然比检验 P < 0.10,以及delta AIC < 0。由于0.05和0.10的阈值存在嵌套关系,当 P < 0.05 时计为2分,从而对具有明确显著性的校正和敏感性Cox证据赋予更高权重。总分范围为0至13分;若得分相同,则按较小的校正Cox P 值排序,其次按较大的3年/5年平均AUC排序。PPARG获得了全部6项初始得分和全部7项临床整合得分(13/13),排名第一。C指数的提升仅作描述性报告,未计入评分。
临床校正
将候选枢纽基因与 TARGET-OS 临床变量整合,包括性别、年龄、诊断时疾病状态、原发肿瘤部位、特定肿瘤区域以及是否接受根治性手术。纳入40例具有完整表达谱和临床记录的患者(包括13例死亡病例)进行临床校正分析。将仅包含临床变量的Cox模型与包含临床变量及基因表达的模型进行比较。采用C指数、赤池信息准则(AIC)和似然比检验的P值评估模型的改善情况。在剔除手术变量后进行敏感性分析。
功能富集与免疫微环境分析
根据 PPARG 表达水平对 TARGET-OS 样本进行分层。利用 PPARG 高表达组与低表达组之间的差异表达结果生成排序基因列表,用于基因集富集分析(GSEA)20。展示的通路包括:Nemeth Inflammatory Response LPS Up、Burton Adipogenesis 5、Burton Adipogenesis 6、Krieg KDM3A Targets Not Hypoxia、Reactome: Transcriptional Regulation By TP53、Fulcher Inflammatory Response Lectin Vs LPS Dn、Hollmann Apoptosis Via CD40 Dn、Zhou Inflammatory Response Live Dn、WP: Fatty Acids And Lipoproteins Transport In Hepatocytes、Sweet Lung Cancer KRAS Up、KEGG Medicus Pathogen HIV Tat To TLR2/4 NF-kB Signaling Pathway 和 Reactome: Fatty Acids。计算 bulk CellAge 衰老评分,并采用 Spearman 相关性分析评估 PPARG 与衰老相关基因或免疫微环境特征之间的关联21。使用非参数检验结合多重检验校正方法,评估 PPARG 高表达组与低表达组之间微环境评分的差异。
单细胞转录组学分析
采用已发表的人类骨肉瘤单细胞转录组数据集进行分析,所用数据对象已完成预处理,包括质控、降维、聚类及人工注释22,23,24。该数据集共包含 68,336 个细胞和 32,297 个基因。为便于正文解读,细胞注释被简化为 13 种主要细胞类型:B 细胞、癌相关成纤维细胞(CAFs)、增殖细胞、内皮细胞、红系细胞、巨噬细胞/单核细胞、恶性骨肉瘤细胞、肌源性细胞、中性粒细胞、类破骨细胞、周细胞、T/NK 细胞以及肿瘤相关基质细胞。通过降维图、特征表达图、点图和小提琴图可视化 PPARG 的分布情况。PPARG 阳性细胞定义为表达水平大于零的细胞。采用 Kruskal-Wallis 检验和 Wilcoxon 秩和检验(经 Benjamini-Hochberg 校正)评估不同细胞类型间的差异。单细胞分析在 R 语言环境中使用 Seurat 软件(版本 5.0.1)完成。
空间转录组学分析
使用 SP_BS3 空间转录组样本构建空间表达对象25,26。质量控制阈值设定为 nFeature_Spatial ≥ 200 且 percent.mt ≤ 30,最终保留 4,572 个有效点用于分析。对数据进行标准化处理,筛选出 3,000 个高变基因,并进行数据缩放、主成分分析、邻域图构建、空间点聚类及降维分析。在去除基因集中的 PPARG 以避免循环相关后,计算了空间 CellAge 衰老评分。构建并评分了内皮细胞、周细胞、巨噬细胞/单核细胞、肿瘤相关基质细胞、恶性骨肉瘤及类破骨细胞的特征基因表达谱。采用斯皮尔曼相关性分析评估 PPARG 表达与各空间评分之间的关联性。以单细胞数据集为参考,空间数据集为查询,进行标签转移分析,以推断每个空间点的预测细胞类型评分23,24。空间转录组分析在 R 中使用 Seurat(版本 5.0.1)完成,所有空间相关性 P 值均采用 Benjamini-Hochberg 方法进行 FDR 校正。
在 GSE36001 中进行外部表达验证
GEO 数据集 GSE36001 仅用作独立的表达验证队列;由于缺乏生存结局数据,未用于预后验证11,27。该数据集包含 19 个骨肉瘤样本和 6 个正常对照样本。使用 GPL6102 平台注释将探针标识符转换为基因符号。当多个探针映射到同一基因时,保留平均表达水平最高的探针。采用 limma28 评估肿瘤组与正常组之间的差异表达。所有分析在 R 中使用 limma(版本 3.56.2)完成,并通过 Benjamini-Hochberg 方法计算校正后的 P 值。
qRT-PCR 与 Western blot 验证
实验验证使用人骨肉瘤细胞系 143B 和人成骨细胞进行。骨肉瘤细胞在含 10% 胎牛血清和 1% 青霉素-链霉素的杜氏改良伊格尔培养基(Dulbecco's modified Eagle's medium)中培养,于 37 °C、含 5% CO₂ 的湿润环境中孵育,并在细胞融合度达 80%–90% 时用 0.25% 胰蛋白酶-EDTA 进行传代。人成骨细胞按照推荐的培养条件进行维持。所有细胞系均经确认无支原体污染。对于 qRT-PCR,使用基于苯酚-胍的 RNA 提取试剂提取总 RNA,并通过分光光度法测定 RNA 浓度和纯度。取 1 微克总 RNA,按照推荐方案使用逆转录试剂进行逆转录。qRT-PCR 采用基于荧光 DNA 结合染料的化学体系,扩增条件如下:95 °C 预变性 30 秒,随后进行 40 个循环(95 °C 变性 5 秒,60 °C 退火延伸 30 秒),并进行熔解曲线分析以确认扩增特异性。每个反应设技术重复三次,共进行三次独立的生物学重复实验。以 GAPDH 作为内参对照,采用 2-ΔΔCt 法计算 PPARG 的相对表达水平29。PPARG 的上游引物序列为 5'-CGAAGACATTCCATTCACAAGAACAG-3',下游引物序列为 5'-AGATGCAGGCTCCACTTTGATTG-3'。
通过蛋白质印迹法检测PPARG蛋白的表达水平。采用添加蛋白酶抑制剂的放射免疫沉淀法缓冲液裂解细胞,并使用二辛可宁酸法测定蛋白浓度。取等量蛋白(每泳道30 µg)经10%十二烷基硫酸钠-聚丙烯酰胺凝胶电泳分离后,转移至聚偏二氟乙烯膜上。室温下用5%脱脂牛奶封闭1小时后,膜在4 °C下与抗PPARG(1:1,000)和抗GAPDH(1:5,000)的一抗孵育过夜,随后在室温下与辣根过氧化物酶标记的二抗(1:5,000)孵育1小时。使用化学发光检测法显影蛋白条带,并重复进行三次独立实验。采用图像分析软件对条带强度进行定量分析30。组间差异采用双尾非配对Student's t检验进行分析。数据以三次独立实验的均值±标准差(SD)表示,P < 0.05被认为具有统计学显著性。实验数据的统计分析使用统计分析软件(版本9.0)完成。
统计分析
除非另有说明,所有生物信息学分析均在 R(版本 4.3.2)中进行。双侧 P 数值 < 0.05 被认为具有统计学意义。相关性分析采用斯皮尔曼等级相关系数(Spearman's rank correlation coefficient)进行评估。ρ)。在适用情况下,采用 Benjamini-Hochberg FDR 方法进行多重检验校正。实验数据以均值表示 ± SD,并采用双尾非配对 t 检验进行比较 t 实验性统计分析使用统计分析软件(版本 9.0)进行。
GSE99671 鉴定出 105 个与细胞衰老相关的差异表达基因(CellAge)
GSE99671 包含来自 18 对配对组织的 36 个样本。经过低计数过滤后,保留了 16,683 个基因。在校正后 P < 0.05 条件下,共有 2,248 个基因差异表达。在更严格的阈值(校正后 P < 0.05 且 |log2FC| ≥ 1)下,有 594 个基因显著差异表达,其中包括 102 个在肿瘤中上调的基因和 492 个下调的基因(图 1A、B)。将 2,248 个差异表达基因与 866 个 CellAge 基因取交集,得到 105 个与衰老相关的差异表达基因(图 1C)。PPARG 在 GSE99671 中呈下调,其 log2FC = -0.644,P = 0.00451,校正后 P 值 = 0.0309。在 18 对样本中的 13 对中,PPARG 在正常组织中的表达高于肿瘤组织,配对 Wilcoxon 检验 P = 0.0294(图 1D)。
多模型预后筛选确定 PPARG 为关键候选基因
在 TARGET-OS 队列中,共纳入 85 例患者,其中发生死亡事件 27 例。通过对 105 个衰老相关差异表达基因进行单变量 Cox 回归、Kaplan-Meier 分析、生存 ROC 分析、LASSO、重复 LASSO 及随机生存森林建模的综合筛选,在未校正临床因素前共获得 18 个候选枢纽基因(图 2A)。单变量 Cox 分析显示,PPARG 与总生存期显著相关(HR = 0.603,95% CI = 0.454–0.802,P = 0.000494,FDR = 0.0447),提示 PPARG 表达水平越高,死亡风险越低。高表达组与低表达组的 Kaplan-Meier 分析结果为 P = 0.00784(图 2B)。1 年、3 年和 5 年的时变 AUC 分别为 0.603、0.760 和 0.776(图 2C)。PPARG 的重复 LASSO 选择频率为 0.920,随机生存森林重要性得分为 0.0398(图 2D–F)。PPARG 的初始枢纽评分达到 6/6,因其满足全部六项预设筛选标准。
临床校正支持PPARG的预后关联
在纳入临床协变量后,PPARG仍与总生存期显著相关(校正后HR = 0.224,95% CI = 0.085–0.589,P = 0.00241;图2G)。仅包含临床变量的模型C指数为0.707,AIC为89.921(补充表1)。加入PPARG后,C指数提高至0.829,AIC降低至78.466,并且根据似然比检验,模型拟合显著改善(P = 0.000244;补充表2,图2H、I)。在敏感性分析中剔除接受根治性手术的病例后,PPARG的保护性关联仍然存在(HR = 0.249,P = 0.00185;补充表3,图2J)。PPARG获得最高的临床整合评分13分(初始枢纽评分6分,加上7个临床整合得分点),其3年和5年AUC最终分别为0.770和0.813(图2K)。
PPARG 相关的功能与免疫微环境特征
PPARG 高表达组与低表达组的 GSEA 分析显示,Nemeth 炎症反应 LPS 上调、Burton 脂肪生成 5、Burton 脂肪生成 6、Krieg KDM3A 非缺氧靶基因、Reactome:TP53 转录调控、Fulcher 炎症反应凝集素对比 LPS 下调、Hollmann CD40 介导的凋亡下调、Zhou 活体炎症反应下调、WP:肝细胞中脂肪酸与脂蛋白转运、Sweet 肺癌 KRAS 上调、KEGG Medicus 病原体 HIV Tat 至 TLR2/4 NF-kB 信号通路,以及 Reactome:脂肪酸富集显著(图 3A)。在整体 TARGET-OS 数据中,PPARG 与总体 CellAge 衰老评分无显著相关性(Spearman ρ = 0.022,P = 0.837;图 3B),但与多个单独的 CellAge 基因存在相关性(图 3C)。免疫微环境分析显示,PPARG 与巨噬细胞(ρ = 0.485,FDR = 2.7 × 10-5)、CD8 T 细胞(ρ = 0.410,FDR = 5.88 × 10-4)、破骨细胞样特征(ρ = 0.383,FDR = 0.00120)、中性粒细胞(ρ = 0.376,FDR = 0.00120)、树突状细胞(ρ = 0.370,FDR = 0.00123)呈正相关。在经过 FDR 校正后,PPARG 高表达肿瘤表现出更高的破骨细胞样、巨噬细胞、CD8 T 细胞、树突状细胞、单核细胞、中性粒细胞、NK 细胞和内皮细胞特征(补充表 4,图 3D、E)。
单细胞转录组学将 PPARG 定位于血管及微环境区室
该单细胞数据集包含 68,336 个细胞和 32,297 个基因。PPARG 表达在不同细胞类型间存在显著差异(图 4A)。内皮细胞(平均表达量 = 0.540;阳性比例 = 44.33%)、周细胞(平均表达量 = 0.439;阳性比例 = 41.61%)、巨噬细胞/单核细胞(平均表达量 = 0.363;阳性比例 = 31.62%)以及肿瘤相关基质细胞(平均表达量 = 0.361;阳性比例 = 44.10%)中观察到最高的平均表达水平(图 4B–E)。一部分恶性骨肉瘤细胞表达 PPARG(平均表达量 = 0.163;阳性比例 = 16.70%),但其在恶性骨肉瘤细胞中的表达水平并未显著高于其他细胞(FDR = 0.151)。这些结果提示,骨肉瘤中 PPARG 的表达主要反映的是血管、髓系及基质微环境状态,而非局限于恶性细胞本身(图 4F)。
空间转录组学将 PPARG 与衰老相关空间状态及血管微环境联系起来
经过质量控制后,保留了 4,572 个 SP_BS3 空间斑点,并将其聚类为七个空间簇(图5A). 图5B 显示 nFeature_Spatial(每个位点检测到的基因)的空间分布。在 866 个 CellAge 基因中,有 845 个在空间表达矩阵中匹配(占 97.58%)。PPARG 呈现局灶性空间表达(图5C)。在去除PPARG后计算的CellAge空间衰老评分显示,其与PPARG表达水平呈微弱但具有统计学显著性的正相关关系(ρ = 0.0692, P = 3.0 × 10-6,FDR = 1.9 × 10-5; 图5D). 空间生态位评分显示PPARG与内皮细胞评分呈正相关(ρ = 0.0433,FDR = 0.00592)和周细胞评分(ρ = 0.0367,FDR = 0.0181),而PPARG与恶性骨肉瘤评分呈负相关(ρ = -0.0592,FDR = 0.000219)和肿瘤间质评分(ρ = -0.0531,FDR = 0.000774; 补充表5, 图 5E)。标签转移分析同样显示与内皮细胞预测评分呈正相关(ρ = 0.0507,FDR = 0.00120)和周细胞预测评分(ρ = 0.0394,FDR = 0.0123),同时与恶性骨肉瘤细胞预测评分呈负相关(ρ = -0.0699,FDR = 1.1 × 10-5; 图5F–H).
外部表达数据及实验验证支持 PPARG 的下调
GSE36001 包含 19 个骨肉瘤样本和 6 个正常对照样本。在骨肉瘤中,PPARG 显著下调(logFC = -1.429,P = 0.00730,校正后 P 值 = 0.0435;图 6A)。在基于细胞的验证实验中,通过 qRT-PCR 检测发现,骨肉瘤 143B 细胞中 PPARG 的 mRNA 表达水平显著低于人成骨细胞(P < 0.001;图 6B)。Western blotting 结果也显示,143B 细胞中 PPARG 蛋白表达显著降低(P < 0.01;图 6C,D)。这些来自外部队列、mRNA 和蛋白水平的结果一致支持 PPARG 在骨肉瘤中表达下调。GSE36001 未包含生存结局数据,因此仅能提供外部表达验证,无法进行独立的预后验证。
数据可用性:
本研究中使用的所有数据集均可公开获取。GSE99671 和 GSE36001 数据来自基因表达综合数据库(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi/acc=GSE99671;https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi/acc=GSE36001)。TARGET-OS 转录组和临床数据从 UCSC Xena 下载(https://xena.ucsc.edu/)。衰老相关基因来自 CellAge:细胞衰老基因数据库,属于人类衰老基因组资源的一部分(https://genomics.senescence.info/cells/)。人类骨肉瘤单细胞和空间转录组数据集来自已发表的图谱及其关联的 GitHub 仓库(https://github.com/zhengxj1/A-Single-Cell-and-Spatially-Resolved-Atlas-of-Human-Osteosarcomas)。富集分析所用的基因集来自 MSigDB(https://www.gsea-msigdb.org/gsea/msigdb/)。本研究生成的处理后数据以及用于复现报告结果的分析脚本已整理并提交为 补充文件 1。

图1:骨肉瘤中差异表达基因及CellAge来源的衰老相关候选基因的鉴定。 (A)火山图显示GSE99671数据集中骨肉瘤组织与配对非肿瘤对照组织之间的差异表达基因。根据预定义的阈值标准,显著上调和下调的基因被突出显示。(B)热图显示GSE99671数据集中代表性差异表达基因在骨肉瘤样本及配对对照样本中的表达模式。(C)维恩图显示GSE99671差异表达基因与CellAge衰老相关基因的交集。(D)GSE99671数据集中PPARG在骨肉瘤组织与匹配的非肿瘤对照组织中的配对表达比较。请点击此处查看该图的放大版本。

图2机器学习与临床校正的生存分析鉴定出PPARG是骨肉瘤中一个核心的衰老相关预后枢纽基因。 (A森林图显示 TARGET-OS 队列中候选衰老相关基因的单变量 Cox 回归结果。B) 高 PPARG 与低 PPARG 患者总体生存率的 Kaplan-Meier 生存曲线比较。C) 基于PPARG预测总体生存率的时变ROC曲线。DLASSO Cox 回归筛选预后候选基因的交叉验证曲线。E) 重复LASSO稳定性分析,显示在300次五折重复中lambda.min的选择频率。F) 由1,000棵树构成的随机生存森林分析,显示变量重要性得分。G森林图显示经临床因素校正的候选枢纽基因Cox回归分析结果。H) 添加单个枢纽基因至临床模型后AIC的变化。 (I) 加入单个枢纽基因后临床模型的C指数改善情况。J) 候选枢纽基因的最终临床整合评分(范围0–13)排序。K) 基于40例患者的临床分析子集绘制的PPARG时间依赖性ROC曲线,显示3年和5年的AUC;1年AUC无法估计。 请点击此处查看此图的放大版本。

图 3:PPARG 相关的功能富集与免疫微环境分析。(A)GSEA 气泡图比较 PPARG 高表达组与 PPARG 低表达组,显示了 Nemeth 炎症反应 LPS 上调、Burton 脂肪生成 5、Burton 脂肪生成 6、Krieg KDM3A 靶基因非缺氧相关、Reactome:TP53 转录调控、Fulcher 炎症反应凝集素对比 LPS 下调、Hollmann CD40 介导的凋亡下调、Zhou 活菌诱导炎症反应下调、WP:肝细胞中脂肪酸与脂蛋白的转运、Sweet 肺癌 KRAS 上调、KEGG Medicus 病原体 HIV Tat 至 TLR2/4 NF-kB 信号通路,以及 Reactome:脂肪酸。(B)PPARG 与整体 CellAge 衰老评分在 bulk TARGET-OS 数据中的相关性。(C)PPARG 与各个 CellAge 基因之间的相关性。(D)PPARG 与免疫微环境特征之间的相关性。(E)PPARG 高表达组与低表达组之间微环境评分的差异。相关性采用 Spearman 秩相关系数(ρ)进行评估,校正后的 P 值使用 Benjamini-Hochberg 方法计算。缩写:GSEA = 基因集富集分析;NF-κB = 核因子 kappa-B;JAK-STAT = Janus 激酶-信号转导与转录激活因子;IL-12 = 白细胞介素-12;FDR = 错误发现率。请点击此处查看该图的放大版本。

图4:骨肉瘤单细胞转录组数据中PPARG在不同细胞区室中的定位。(A)经简化人工注释后的人骨肉瘤单细胞转录组数据集中主要细胞类型的UMAP可视化图。(B)FeaturePlot展示PPARG在单细胞中的整体表达分布。(C)DotPlot展示PPARG在主要细胞类型中的表达情况。(D)小提琴图展示不同细胞类型中PPARG的表达水平。(E)柱状图展示每种主要细胞类型中PPARG阳性细胞的比例。(F)UMAP可视化图展示恶性骨肉瘤细胞中PPARG的表达。缩写:UMAP = 均匀流形近似与投影。请点击此处查看该图的放大版本。

图5:骨肉瘤中PPARG及衰老相关空间特征的空间转录组定位。(A)SP_BS3骨肉瘤空间转录组切片中由转录组定义的簇的空间分布。(B)每个空间点检测到的基因数目的空间分布,以nFeature_Spatial表示。(C)PPARG在SP_BS3各空间点中的空间表达模式。(D)基于CellAge的衰老评分的空间分布。(E)PPARG表达与基于CellAge的衰老评分或细胞生态位评分之间的相关性分析。(F)标签转移预测图,显示每个空间点中占主导地位的单细胞来源的细胞类型。(G)PPARG表达、基于CellAge的衰老评分以及基于标签转移的细胞类型预测评分之间的相关性分析。(H)PPARG高表达与PPARG低表达空间点的空间分布。请点击此处查看该图的放大版本。

图6: 骨肉瘤中PPARG下调的外部表达及实验验证。 (A)箱线图显示GSE36001数据集中骨肉瘤样本(n = 19)与正常对照样本(n = 6)的PPARG表达水平。(B)人骨肉瘤143B细胞和人成骨细胞对照细胞中PPARG mRNA表达的qRT-PCR分析。(C)代表性Western blot图显示人成骨细胞对照细胞和骨肉瘤143B细胞中PPARG和GAPDH蛋白的表达。GAPDH用作上样对照。(D)Western blot条带的密度测定定量结果,显示相对于GAPDH标准化后的PPARG蛋白相对水平。在(B)和(D)中,数据表示为三次独立实验的均值 ± SD。P < 0.01 且 P < 0.001,与人成骨细胞对照组相比,采用双尾非配对Student’s t-检验确定。缩写:qRT-PCR = 定量反转录聚合酶链式反应;GAPDH = 甘油醛-3-磷酸脱氢酶;SD = 标准差。请点击此处查看该图的放大版本。
补充表 1:仅临床变量的 Cox 模型在 TARGET-OS 队列中的性能。 使用性别、年龄、诊断时疾病状态、原发肿瘤部位、特定肿瘤区域和根治性手术状态等临床变量单独构建的 Cox 模型的 C 指数、AIC 及模型摘要。缩写:AIC = 赤池信息准则。 请点击此处下载该文件。
补充表 2:仅临床因素模型与临床因素加基因模型的 Cox 模型比较。 在临床模型中逐个加入候选枢纽基因后的模型比较结果,包括 C 指数、AIC、似然比检验统计量及模型改进指标。缩写:AIC = Akaike 信息准则。 请点击此处下载该文件。
补充表3:移除根治性手术变量后的敏感性分析。 敏感性Cox回归结果,用于评估在调整后的临床模型中排除根治性手术变量后,候选枢纽基因(特别是PPARG)的预后关联是否仍然稳定。 请点击此处下载该文件。
补充表 4:TARGET-OS 中与 PPARG 相关的免疫和基质微环境特征。 PPARG 表达水平与免疫、基质、血管、炎症及与衰老相关分泌表型(SASP)相关的单样本基因集富集分析(ssGSEA)特征之间的相关性及组间比较结果,包括斯皮尔曼相关系数、P 值、校正后的 P 值,以及 PPARG 高表达组与低表达组的比较。缩写:SASP = 衰老相关分泌表型;ssGSEA = 单样本基因集富集分析。 请点击此处下载该文件。
补充表 5:SP_BS3 中 PPARG 的空间转录组相关性分析。 PPARG 表达与 SP_BS3 骨肉瘤空间转录组切片中基于空间 CellAge 的衰老评分、细胞生态位评分以及基于标签转移的细胞类型预测评分之间的相关性分析结果。 请点击此处下载该文件。
本研究通过整合差异表达分析、CellAge基因交集、TARGET-OS生存模型构建、临床校正、多组学定位及实验验证,提名PPARG作为骨肉瘤中一个潜在的与衰老相关的预后基因。该分析框架与当前强调基于分子信息的骨肉瘤研究以及用于解读衰老相关生物学的已整理衰老基因资源保持一致14,31。相较于正常组织,PPARG在骨肉瘤中表达下调,且在TARGET-OS队列中,较低的PPARG表达与较差的总生存率相关。这些发现提示,PPARG不仅在骨肉瘤中存在转录水平的改变,还可能携带具有临床意义的预后信息。然而,由于候选基因筛选和模型评估均在同一TARGET-OS队列(85例患者,27例死亡事件)中进行,所观察到的模型性能提升可能存在乐观偏倚和过拟合风险;因此,PPARG的预后价值在经独立的骨肉瘤生存队列验证之前,应被视为假设生成性的。尽管如此,PPARG的生物学功能应结合具体背景进行解读,因为已有实验研究表明,PPAR-γ/核受体调控具有抗肿瘤作用,同时也有报道指出PPARG相关的破骨细胞程序可能促进疾病进展32,33,34。
一个重要的细节是,不应将PPARG解释为整体CellAge评分的简单替代指标。在批量的TARGET-OS数据中,PPARG与全局CellAge衰老评分无显著相关性;而在空间转录组学数据中,PPARG与排除PPARG后计算的CellAge评分呈微弱但显著的相关性。这种差异可能反映了批量数据中的细胞组成效应、衰老基因集的多功能性,以及空间位点中微环境生态位的局部富集。共识性研究和转录组学研究均强调,细胞衰老具有异质性、动态性,并依赖于细胞类型、应激源和组织背景,而SASP程序在癌症进展过程中可能产生相反的作用25,35,36,37。因此,PPARG被保守地定义为一种源自CellAge的衰老相关预后基因,而非已被证实的衰老驱动因子。相应地,PPARG的衰老相关性归类是基于其属于CellAge基因集,而非其在衰老过程中具有已证实的机制性作用,因此不应将PPARG用作整体衰老活性的定量替代指标。
总体而言,这些空间分析结果阐明了统计学显著性与生物学相关性之间的区别。在包含4,572个空间点位的情况下,即使是非常微弱的相关性也可能超过传统的显著性阈值;例如,PPARG与空间CellAge衰老评分之间的相关性(ρ = 0.0692)仅解释了约0.48%的方差,但由于点位数量庞大,提供了充足的统计效能,其P值仍达到3.0 × 10-6。因此,此类点位水平的关联应被视为在统计学上可检测但生物学效应较弱、具有假设生成意义的信号,其生物学解释应基于效应量大小,而不仅仅是P值。在大规模点位水平数据集中,统计学显著性不应等同于强烈的生物学效应。
单细胞和空间分析为PPARG提供了生物学背景。PPARG在内皮细胞、周细胞、巨噬细胞/单核细胞以及肿瘤相关间质细胞中富集,并且在空间上与内皮细胞和周细胞评分相关。这些结果提示,PPARG所携带的预后信息可能与骨肉瘤中的血管和微环境组分有关。该解释与近期的骨肉瘤单细胞/空间图谱以及更广泛的证据一致,即免疫、血管和间质生态位共同塑造了肿瘤异质性、治疗反应及疾病进展25,38,39。由于恶性骨肉瘤细胞仅在部分细胞中表达PPARG,且其表达水平并未显著高于其他细胞类型,因此仅从肿瘤细胞角度进行解释将是不完整的。相反,在GSE99671和GSE36001数据集中观察到的PPARG在整体水平上的下调,可能部分反映了肿瘤组织与正常组织之间间质、血管、骨髓、脂肪生成或免疫细胞组成差异,而非恶性骨肉瘤细胞中真正的下调。由于在整体分析中未明确校正肿瘤纯度和细胞类型丰度,这一可能性无法排除,需进一步开展针对性研究加以验证。
功能与免疫微环境分析结果与此解释一致。PPARG 与巨噬细胞、CD8 T细胞、树突状细胞、中性粒细胞、NK细胞及内皮细胞特征基因表达谱相关,GSEA分析揭示了与炎症、NF-κB、JAK-STAT/IL-12信号通路、TP53调控以及DNA损伤检查点相关的通路。综合来看,这些结果提示PPARG可能标记了一种复合的微环境状态,涉及衰老相关应激、免疫细胞浸润以及血管和基质组分。这一解释在生物学上是合理的,因为PPAR-γ在抑制巨噬细胞/单核细胞的炎症活化方面具有明确作用,包括对AP-1、STAT和NF-κB相关转录程序的影响;此外,骨肉瘤的免疫微环境中包含具有促瘤和抑瘤双重功能的髓系、淋巴系及血管成分38,40,41。这些单细胞和空间层面的观察结果具有描述性和假设生成性;其本身并未确立血管微环境机制、衰老程序或预后相关通路。
本研究存在若干局限性,需予以说明。首先,主要的预后分析基于回顾性的 TARGET-OS 公共队列,样本量和事件数量有限;因此,PPARG 的预后价值应依据公认的肿瘤标志物报告与验证原则,在独立队列中进一步验证42。具体而言,候选基因筛选在包含 27 例死亡病例的 85 例 TARGET-OS 患者中进行,而临床校正模型的评估则使用了其中重叠的 40 例患者亚组(含 13 例死亡)。由于两项分析均源自同一队列,所报告的 C 指数提升(从 0.707 至 0.829)和 AIC 下降可能偏于乐观。目前尚无独立的生存队列可用于外部验证 PPARG 的预后价值;GSE36001 仅用于肿瘤与正常组织表达的验证。此外,整合枢纽评分(integrated hub score)是一种探索性的内部排序启发式方法,而非经过验证的预后工具。第二,空间转录组学分析仅基于单个 SP_BS3 样本;尽管相关性具有统计学意义,但效应量较小,需在更多空间样本中加以验证。特别是,增加来自独立患者的空間轉錄組樣本數量,對獲得這些微弱關聯的更穩健估計至關重要,未來有必要開展基於更大空間隊列的研究。空間轉錄組學提供了寶貴的原位(in situ)分子背景,但其解讀仍受平台分辨率、取樣策略、組織質量及計算整合方法的影響43。此外,空間分析依賴於一個預處理的公共單細胞數據對象,其註釋較為簡化,並基於計算評分與標籤轉移方法;鑒於效應量極小,這些數據僅支持描述性定位結論,無法支持關於血管微環境或衰老程序的機制性推論。第三,實驗驗證基於骨肉瘤細胞系 143B 和人骨細胞;仍需更多骨肉瘤細胞系和臨床樣本進行驗證,單一細胞系的結果無法確立細胞類型特異性、臨床預後相關性或衰老生物學意義。第四,本研究僅展示了關聯性,而非因果關係。需在骨肉瘤細胞及微環境模型中進行功能干預實驗,以確定 PPARG 是否直接調控衰老相關程序、血管微環境或腫瘤進展。第五,衰老相關性是根據與 CellAge 基因集的重疊來定義的,而 PPARG 在 TARGET-OS 中與整體 CellAge 衰老評分無顯著相關性。此外,腫瘤與正常組織的整體比較未對腫瘤純度或細胞類型組成進行校正,因此觀察到的下調可能部分反映的是微環境組成差異,而非惡性細胞內在的改變。
PPARG 是一个源自 CellAge 的基因,其表达水平降低与 TARGET-OS 队列中较差的总体生存率相关。在骨肉瘤整体组织比较中,PPARG 表达下调,同时其在血管、髓系和基质组分中富集,提示整体组织中 PPARG 水平可能部分反映了微环境中的细胞组成,而非肿瘤细胞固有的表达。单细胞和空间分析结果为描述性发现,且该预后关联尚未在生存数据上进行独立验证,需在具有生存结局的外部队列中进一步验证。这些结果提示 PPARG 可作为骨肉瘤预后评估及衰老相关微环境研究的候选生物标志物,但有待外部验证和功能研究确认。在空间转录组学中,PPARG 与衰老相关或血管生态位评分之间的相关性效应量较弱,尽管由于空间点数量较大而具有统计学显著性,因此应谨慎解读。
作者声明无竞争利益。
作者贡献:
Yongwen Li 和 Wentao Qin 构思并设计了本研究。Yongwen Li 进行了生物信息学和计算分析。Tuo Liang 完成了实验验证。Rubiao Qiu 和 Zide Zhang 参与了图表制作。Rubiao Qiu 和 Zide Zhang 对研究进行了指导,并对稿件提出了关键性修改意见。所有作者均审阅并批准了最终稿件。
作者衷心感谢 GEO、TARGET-OS、UCSC Xena、CellAge、MSigDB 以及人类骨肉瘤单细胞与空间转录组图谱项目的研究人员和贡献者,感谢他们提供了公开可用的数据集和资源,使得本研究得以开展。本研究得到了广西自然科学基金(编号:2023GXNSFAA026111)的资助。
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| 抗GAPDH一抗 | Proteintech Group, Wuhan, China | 10494-1-AP | 用于Western blot中检测GAPDH作为上样对照的一抗;稀释比例为1:5,000。 |
| 抗PPARG一抗 | Proteintech Group, Wuhan, China | 16643-1-AP | 用于Western blot检测PPARG蛋白的一抗;稀释比例为1:1,000。 |
| BCA蛋白检测试剂盒 | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | PC0020 | 比色法检测,用于电泳前测定总蛋白浓度。 |
| DESeq2 | Bioconductor | Version 1.40.2 | 用于基于计数的转录组数据差异基因表达分析的R软件包。 |
| 杜尔贝科改良伊格尔培养基(DMEM) | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | 11995 | 基础培养基,用于维持143B骨肉瘤细胞的培养。 |
| ECL检测试剂 | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | PE0010 | 化学发光底物,用于Western blot中HRP标记二抗的检测。 |
| 胎牛血清(FBS) | Zhejiang Tianhang Biotechnology Co., Ltd. (Sijiqing), Huzhou, Zhejiang, China | 11011-8611 | 添加至培养基中的血清补充物,用于支持细胞生长和存活。 |
| glmnet | CRAN | Version 4.1-8 | 用于惩罚回归分析(包括LASSO和弹性网络建模)的R软件包。 |
| GraphPad Prism | GraphPad Software, San Diego, CA, USA | Version 9.0 | 用于统计分析、图表生成和实验数据可视化的软件。 |
| HRP标记二抗 | Proteintech Group, Wuhan, China | SA00001-2 | 辣根过氧化物酶标记的二抗,用于Western blot检测;稀释比例为1:5,000。 |
| 人成骨细胞 | Cell Applications, Inc., San Diego, CA, USA | 406-05A | 原代人成骨细胞,用作非恶性对照/比较细胞类型。 |
| 人骨肉瘤143B细胞系 | American Type Culture Collection (ATCC), Manassas, VA, USA | CRL-8303 | 人骨肉瘤细胞系,用于体外验证实验和分子检测。 |
| ImageJ | National Institutes of Health (NIH), USA | Version 1.53 | 图像分析软件,用于实验图像的定量分析。 |
| limma | Bioconductor | Version 3.56.2 | 用于差异表达分析和基于线性模型的统计检验的R软件包。 |
| 青霉素-链霉素 | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | P1400 | 细胞培养基中使用的抗生素补充剂,用于减少细菌污染。 |
| PPARG和GAPDH引物 | Sangon Biotech (Shanghai) Co., Ltd., Shanghai, China | 定制合成;序列见方法部分 | 用于qPCR分析PPARG和GAPDH表达的定制寡核苷酸引物。 |
| PVDF膜,0.45 µm | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | YA1701 | 用于Western blot中蛋白转膜的膜材料。 |
| R统计软件 | R Foundation for Statistical Computing, Vienna, Austria | Version 4.3.2 | 用于生物信息学分析、模型构建和可视化的统计计算环境。 |
| randomForestSRC | CRAN | Version 3.2.2 | 用于随机生存森林建模和特征重要性分析的R软件包。 |
| 逆转录试剂盒 | Beyotime Biotech Inc., Shanghai, China | D7168M | 用于在定量PCR前由提取的RNA合成互补DNA(cDNA)。 |
| RIPA裂解缓冲液 | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | R0010 | 用于裂解细胞以进行Western blot分析的蛋白提取缓冲液。 |
| Seurat | Satija Laboratory | Version 5.0.1 | 用于单细胞RNA测序数据处理、整合、聚类和可视化的R软件包。 |
| survival | CRAN | Version 3.5-7 | 用于生存分析(包括Cox比例风险模型)的R软件包。 |
| SYBR Green qPCR预混液 | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | SR1110 | 用于定量实时PCR扩增的荧光预混液。 |
| timeROC | CRAN | Version 0.4 | 用于生成时间依赖性受试者工作特征曲线并计算随时间变化的预测性能的R软件包。 |
| TRIzol试剂 | Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA | 15596026CN | 用于从培养细胞中提取总RNA的试剂。 |
| 胰蛋白酶-EDTA,0.25% | Beijing Solarbio Science & Technology Co., Ltd., Beijing, China | T1300 | 用于传代和收获贴壁细胞的细胞解离试剂。 |