DNA甲基化+转录组+肿瘤微环境，深度分析识别区分癌症预后

Integrative analysis of immune-related multi-omics profiles identifies distinct prognosis and tumor microenvironment patterns in osteosarcoma

对免疫相关的多组学特征的综合分析确定了骨肉瘤的独特预后和肿瘤微环境模式

发表期刊：Mol Oncol

发表日期：2021 Dec 11

影响因子：6.603

DOI:  10.1002/1878-0261.13160

一、背景

        人类骨肉瘤（OS）是最常见的、侵袭性的原发性骨肿瘤之一，在青少年和青年中最为流行，主要影响长骨，如股骨远端、胫骨近端和肱骨。尽管明确的手术切除结合辅助化疗已经明显改善了局部OS患者的预后，但约20%的OS患者在最初诊断时就患有肺转移性疾病。

        DNA甲基化的改变在癌症的发展中起着至关重要的作用。在肿瘤发生过程中，高甲基化启动子区域可导致肿瘤抑制基因的表观遗传沉默，而非启动子元件的异常DNA甲基化是造成肿瘤内异质性的重要因素。

        免疫系统在抑制肿瘤发生和发展中的重要作用，使免疫治疗成为与手术、化疗和放疗并列的癌症治疗的第四大支柱。然而，并不是所有的患者都能从这些疗法中完全受益，其中许多患者，包括OS患者，未能获得完全的反应或经常复发。

二、材料与方法

1、数据来源

1）从TARGET获得OS患者的临床信息、RNA-seq（基因表达）、DNA甲基化阵列和DNA拷贝数变异数据：从总共98名患者中提取了计数和TPM（每百万转录本）格式的RNA-seq数据；从总共88名患者中提取了DNA拷贝数片段

2）公开的数据集进行独立验证：GSE21257数据集OS患者的基因表达谱；E-MTAB-9875数据集和E-MTAB-7263数据集软骨肉瘤患者的DNA甲基化谱

3）结合以前报道的两个数据库的基因组来构建免疫学基因列表：ImmuneSigDB（通过MSigDB集合）和ImmPort

4）IMvigor210队列

2、分析流程

1）对OS的免疫相关DNA甲基化模式进行聚类：单变量Cox回归分析；基于K-means算法的无监督共识聚类

2）识别差异表达基因和差异甲基化的CpG位点：应用deseq2 r软件包来处理RNA-seq count数据，然后识别两组之间的差异表达基因（DEG）

3）DNA甲基化与基因表达的相关性

4）评估TME细胞的浸润丰度：CIBERSORT、quanTIseq、xCell评估免疫浸润程度

5）功能富集分析：clusterprofiler和goplot r软件包被用于过度代表分析、预排序基因组富集分析（GSEA）和可视化；使用gsva r软件包进行多组间患者的非参数GSVA

6）IMP相关特征评分模型的构建和验证：单变量Cox回归分析；Lasso回归分析；使用Akaike信息准则（AIC）方法进行逐步回归；生存分析；构建nomogram

7）识别具有体细胞拷贝数改变的复发性区域：应用GISTIC 2.0来分析DNA拷贝数分割图

8）蛋白质-蛋白质相互作用的评价：基于STRING数据库构建了一个蛋白质-蛋白质相互作用（PPI）网络

9）对药物敏感性的IMP相关签名得分的评估：癌症药物敏感性基因组学（GDSC）数据库

10）泛癌TME标志物表达模式的聚类分析：通过GSVA计算富集度分数；应用无监督聚类算法来分析96名OS患者的TME特征的标准化GSVA分数

11）对免疫治疗反应的预测能力的综合评价：TIDE；TIS

三、实验结果

01 - 识别与OS患者预后相关的免疫相关DNA甲基化模式（IMPs）

        从TARGET下载了所有OS组织患者的临床信息和DNA甲基化图谱，并通过champ R包对甲基化β值矩阵进行了归一化。共有84例患者与DNA甲基化图谱和临床信息相匹配。作者通过单变量Cox回归分析发现，39335个CpG位点的甲基化水平与OS患者的总生存期相关。参照Illumina HumanMethylation450k注释文件，作者根据UCSC基因组注释获得了每个甲基化探针的靶基因，并进一步确定了39335个CpG位点中的25924个受免疫学特征基因集。提取了这些免疫学CpG位点的甲基化β值矩阵进行无监督的共识聚类。

        如图1A，将OS患者分为三个聚类，每个聚类中都有合理数量的患者（n = 32，26和26分别在亚群1，2和3）。进行了Kaplan-Meier生存分析，不同IMP的患者有不同的总生存概率。如图1B所示，亚群2的患者的预后令人满意，亚群3中的患者预后极为不利，2年随访的总生存概率下降到50%，5年随访的总生存概率下降到25%。此外，比较了OS患者的临床病理特征（图1C,D），三个IMP之间的年龄、性别和种族没有明显差异。值得注意的是，观察到亚群3的患者对化疗的组织学反应不理想，而且转移率较高。

图1 骨肉瘤中免疫相关甲基化模式（IMPs）的鉴定及不同IMPs的免疫浸润特征

        对两个独立的DNA甲基化队列（E-MTAB-9875和E-MTAB-7263）进行了验证，提取相同的免疫CpG位点并进行无监督共识聚类。对于E-MTAB-9875队列，为了避免非生物性的变异，排除了少量由450KBeadChip测试的OS患者，共有219名经EPIC BeadChip检测的OS患者被用于聚类。观察到三种IMP之间的肿瘤差异等级分布明显不同（图S1D-F）。对于E-MTAB-7263队列，所有102名软骨肉瘤患者被聚类为三个IMP，通过Kaplan-Meier生存分析，发现三个IMP之间存在明显的预后（图S1G-I），其中亚群3的患者的总生存期极差。

图S1

        作者提取了每个样本的转录组特征，并应用两种独立的去卷积算法：CIBERSORTx和quanTIseq，以计算TME中的特定免疫细胞组成。两种分析方法的结果在三个IMP的TME中呈现出明显不同的免疫浸润特征，表明不同IMP的患者具有不同的免疫细胞浸润景观。

        如图1E所示，与其他两个亚群相比，CD8+T细胞浸润在亚群2患者中明显较高，而T细胞CD4naive浸润在亚群3患者中则明显较高。此外，亚群2的患者被发现有明显较高的活化NK细胞水平，而静止NK细胞的水平较低，亚群3患者的γδT细胞的浸润水平相对较高。quanTIseq分析显示，在亚群3中巨噬细胞、中性粒细胞、CD4+T细胞的浸润明显较低，而B细胞、DCs的浸润较高（图S1J）。此外，应用xCell算法来评估不同集群中整体的TME组成，发现亚群2中的TME得分和免疫得分较高（图S1K）。

        作者进行了Kaplan-Meier生存分析，以考察渗透水平对OS患者预后的每种细胞类型的潜在有益或有害意义。如图S1L所示，观察到T细胞CD4记忆静止期、T细胞滤泡辅助期、T细胞CD8、单核细胞和巨噬细胞M2的高浸润水平对OS患者的总生存率有益，而T细胞CD4幼稚期、NK细胞静止期和树突状细胞静止期的浸润水平可能有害。综合来看，T细胞CD8、T细胞CD4 naive、NK细胞静止、单核细胞、树突状细胞静止和巨噬细胞M2浸润在IMPs中的分布和对生存的潜在影响是一致的。

        基因组变异分析（GSVA）表明，一些与癌症相关的标志和途径在IMPs之间存在差异，如GSVA富集分数的热图所示（图1F)。综上所述，在亚群2患者中发现了免疫浸润升高、细胞毒性潜能增强和抗肿瘤免疫反应激活，而亚群3则呈现免疫抑制性TME状态，因此，第2组的预后较好，而第3组的生存期较短。

        关于亚群2和3之间最明显的总生存率和TME浸润细胞组成，作者进行了差异甲基化探针（DMPs）和差异表达基因（DEGs）分析。与亚群2相比，亚群3中共有32805个DMPs（15410个低甲基化和17395个高甲基化）和2747个DEGs（1516个上调和1231个下调）被确定（图S2A,B）。其中，8176个DMPs和1910个DEGs形成了免疫学特征基因集（图2A,B）。

图2 IMPs的DNA甲基化和基因表达的特点

        为了研究DNA甲基化对1910个免疫学DEGs基因表达的调控作用，作者进行了Pearson相关分析，以检查顺式和反式调控。对于顺式调控效应，即基因表达通常与其启动子区域的DNA甲基化水平有关，发现对应于754个基因的2047个CpG位点的甲基化水平与基因表达水平明显相关（1309个负相关和738个正相关）。对IMP群组2和3中的DEGs和DMPs进行了跨调控效应的分析，即一个基因的启动子区域的甲基化水平与其他基因的表达水平相关联。发现大多数DEGs（1910个中的1638个）被免疫相关的DMPs转调。

02 - 构建IMP相关的特征评分模型

        作者试图建立一个特征评分模型来定量评估IMPs，并评估个别OS患者的预后情况。通过单变量Cox回归分析了IMPs 亚群2和3之间的1910个免疫学DEGs，发现662个基因的表达水平与总生存期相关。Lasso Cox回归分析确定了9个基因用于构建预后风险模型。基于使用Akaike信息准则（AIC）的逐步回归方法，构建了一个六种基因的预后风险评分模型，其中有最优化的候选基因（图3A,B）：IMP_Risk评分=MYC的表达水平*0.4998+COL13A1的表达水平*（0.2715）+UHRF2的表达水平*（0.3338）+MT1A的表达水平*（0.2558）+ACTB的表达水平*（-0.4997）+GBP1的表达水平*（-0.2012）。如随时间变化的ROC曲线（图3C）所示，IMP-Risk评分对1、3、5年总生存的预测AUC分别为0.827、0.822、0.858。对于训练队列的总体中位生存时间（3862天，图S2C），AUC为0.912，将IMP-Risk评分的0.751定义为划分高和低IMP Risk组的分界值。如图3D和图S2D所示，Kaplan-Meier生存分析显示，高IMP-Risk患者的总生存期明显较短。图S2E和图3E显示了高风险组和低风险组之间的风险评分、生存状态和六种基因表达谱的分布。

图3 骨肉瘤中IMP相关特征评分模型的构建和验证

        如图3F所示，IMP-Risk评分、组织学反应和诊断时的转移被确定为总生存率的独立预后因素。通过对1年、3年、5年总生存期和Kaplan-Meier生存分析的时间依赖性ROC分析，验证了IMP相关签名评分模型对无复发生存的良好预测能力（图S2F）。

        共有53名来自GSE21257的OS患者被用于验证IMP相关的特征评分模型，使用相同的公式和归一化的基因表达数据计算每个患者的IMP-Risk评分。通过对总体中位生存时间（5670天，图S2G）的时间依赖性ROC分析，确定了IMP-Risk评分的5.95作为划分高、低IMP Risk组的分界点（图S2H）。如图3G,H和图S2I所示，同样进行了时间依赖的ROC和Kaplan-Meier生存分析，发现IMP相关的特征评分模型对总生存期和无转移生存期都保持了良好的预后预测能力。

        对于TARGET-OS数据集，作者进一步整合了IMP Risk评分、组织学反应和诊断时的转移，构建了一个nomogram，结果显示IMP Risk评分是一个主要的风险贡献因素（图3I）。如校准曲线所示，该nomogram为OS患者的总生存期提供了理想的预测精度（图S2J）。此外，进行了决策曲线分析（DCA），结果显示，与诊断时的组织学反应和转移相比，IMP_Risk评分的净收益明显更高（图S2K）。

图S2

03 - 与IMP相关特征得分相关的分子特征

        使用Pearson相关分析，分析了IMP相关特征评分模型中包括的6个基因之间的正或负共表达（图4A）。作者发现MYC、COL13A1、UHRF2和MT1A的表达水平与特定启动子CpG位点的甲基化水平呈负相关（图4B，图S3A）。在高风险组和低风险组之间共发现746个上调和876个下调的DEGs，其中474个上调和579个下调的基因与免疫学特征基因组有关（图S3B）。通过预排序的GSEA，发现几个生物过程和途径在低IMP Risk组中被明显激活，而缺氧、MYC靶点、Wnt/β-catenin信号、酪氨酸代谢和其他肿瘤发生相关途径在高IMP Risk组中更为富集（图4C，图S3C）。

图S3

        通过Pearson相关分析，确定218个免疫DEGs为IMP Risk分数高的相关基因。热图显示，218个免疫DEGs的表达与IMP Risk评分相关（图4D）。过度代表分析确定了这些基因富集的生物功能和途径，如IFN-γ信号、T细胞增殖和活化、细胞-细胞粘附、SLIT/ROBO调节和NOTCH1/3信号途径（图4E,F）。此外，通过STRING数据库建立了一个PPI网络，并在这些基因之间建立了前瞻性的蛋白质-蛋白质相互作用。如图4G所示，通过cytoscapemcode和hubba插件确定了PPI网络的前5个枢纽基因和前3个枢纽集群。

图4 骨肉瘤分子特征的综合分析和IMP相关的特征评分模型

        作者通过GISTIC 2.0评估了高风险和低风险OS患者之间体细胞拷贝数交替（SCNAs）的分歧，并确定了不同IMP Risk组之间有明显扩增或缺失的特定染色体区域。如图4H所示，1号和8号染色体的扩增伴随着13号染色体的缺失在高风险组中富集，而14号和17号染色体的扩增伴随着3号染色体的缺失在低风险组中富集。在高危患者中发现了局灶性扩增峰，包括已被研究的癌症驱动基因MYC（8q.24.21）和几个抗凋亡基因（MCL1、HORMAD1和1q21.2的ECM1），以及13q14.2的局灶性缺失峰。值得注意的是，在低风险组中发现了包括多个TCR相关基因（14q11.2和17p11.2）的病灶扩增峰。

04 - 参与药物基因组学相互作用的IMP相关特征分数

        作者应用Spearman相关分析，根据GDSC数据库，确定了34个药物敏感性与IMP Risk评分相关的配对，包括PI3K/MTOR抑制剂AZD805，CDK抑制剂PHA-793887，以及RTK抑制剂舒尼替尼（图5A）。还确定了四对耐药性与IMP Risk评分相关的药物，包括Src和Abl激酶抑制剂saracatinib和EGFR抑制剂cetuximab。此外，询问了对IMP Risk评分敏感的药物的假定靶点和信号通路，发现它们大多针对JNK/p38、ERK/MAPK、PI3K/mTOR和细胞周期信号通路（图5B）。

图5 与IMP相关的特征分数和药物敏感性之间的关系

05 - 与IMP相关的特征分数和TME模式在OS中的关系

        Bagaev等人定义了四种不同的泛癌症TME亚型，即免疫富集、纤维化（IE/F）；免疫富集、非纤维化（IE）；纤维化（F）；和耗竭（D）。在四种TME模式中，两种免疫富集模式的患者预后相对较好，尤其是TME_IE模式的患者在特定类型的癌症中预后更好。然而，TME_D模式的病人在泛癌症数据集中的预后一直很差。同样，作者利用无监督聚类方法，通过使用29个功能基因表达特征（FGES）集的稳健标准化的GSVA富集分数，将TARGET OS患者分成四组（图S4A和图6A）。如图6B所示，t-SNE分析显示GSVA结果在四个TME模式中的明显分布。如热图所示，每个OS样本的详细标准化GSVA富集分数显示了四种TME模式中不同的FGES特征（图6C）。总的来说，在TARGET中，三个IMP群组、高/低IMP相关的特征评分风险水平和四个TME模式在OS患者中显示出明显的一致关系（图6D-G）。

图6 与IMP相关的特征得分揭示了骨肉瘤中不同的肿瘤微环境模式

06 - 利用IMP相关特征评分模型预测免疫疗法的反应

        作者试图研究IMP相关的特征评分模型是否能预测OS患者对免疫治疗的反应，利用TIDE模块来评估免疫疗法对不同IMP_风险组患者的潜在临床疗效。TIDE算法评估了T细胞功能障碍和T细胞排斥的表达特征，以评估肿瘤免疫规避，并将其整合为TIDE总分。较高的TIDE预测得分代表了较高的免疫规避潜力，表明患者不太可能从ICI治疗中获益。此外，TIDE模块分析多种特征来估计肿瘤免疫规避，如与MDSC、TAM或CAF特征的相关性。

        使用TARGET OS数据集，没有发现高风险组和低风险组之间的总体TIDE评分有任何明显的差异（图S4B）。尽管如此，作者发现高危组的T细胞排除分数相对较高，但T细胞功能障碍分数较低（图7A）。对于TIDE产生的其他特征，发现IFNG特征的低表达与高危组的M2 TAM、MDSCs和CAFs特征相关（图7B，C），较高比例的CTL与低风险组有关（图S4C）。此外，评估了T细胞炎症特征（TIS）的得分，发现高危组的TIS得分相对较低（图7D）。这些结果表明，高危组的OS患者可能具有更高的免疫抑制TME状态的潜力，对ICI治疗的反应更差。

        最后，使用IMvigor210队列，PD-L1治疗的随访数据集，评估IMP_Risk评分在预测ICI治疗反应方面的能力，发现低IMP_Risk的患者在PD-L1治疗后获得了明显的临床获益，治疗反应更好，总生存期明显延长（图S4D和图7E-H）。还发现低IMP_Risk组的肿瘤突变负担（TMB）和肿瘤新抗原负担（TNB）明显较高（图7I,J），这至少可以部分解释低IMP_Risk组的优势和ICI治疗的更大获益。

图7 与IMP相关的特征评分模型对免疫疗法反应的预测价值

四、结论

        作者对免疫相关的DNA甲基组和转录组的综合分析显示，TME的影响及其与OS预后的关系存在广泛的调控关系。构建并验证了IMP相关的特征评分模型，记录了特征基因在转录和体细胞拷贝数改变中的串扰和调节作用，并确定了它们在靶向治疗和免疫治疗中的潜在作用。本研究工作还强调了DNA甲基化在肿瘤细胞和TME之间串联的重要临床意义，并应有助于为OS患者开发个性化的免疫治疗策略。