多组学数据整合和建模揭示了胰腺癌的新机制并改善了预后预测

Multi-omics data integration and modeling unravels new mechanisms for pancreatic cancer and improves prognostic prediction

多组学数据整合和建模揭示了胰腺癌的新机制并改善了预后预测

发表期刊：NPJ Precis Oncol

发表日期：2022 Aug 17

影响因子：10.092

DOI:  10.1038/s41698-022-00299-z

一、背景

        胰腺导管腺癌（PDAC）是最具侵略性的肿瘤之一，远端转移的患者预后最差。定义PDAC预后和治疗反应的标志是由肿瘤细胞及其微环境部分的进展和串联决定的，上皮肿瘤细胞已被广泛地在分子水平上进行分型。

        基于上皮细胞的特征，一组中心转录因子（TFs）定义了肿瘤细胞的表型。经典亚型的特征是导管细胞和生殖系TFs，如PDX1、HNF4A、HNF1A和GATA6，而基底样亚型是由SNAI2、SIX1/4和TP63调节。其他TFs的诱导取决于肿瘤内和肿瘤外的因素，如缺氧和上皮-间质转化（EMT）相关的蛋白，这有助于肿瘤的侵略性。

二、材料与方法

1.数据来源

1） 注册号为NCT01692873，从三个专家临床中心获得，使用了90个异种移植（PDX）

2） RNA-Seq表达数据集：TCGA-PAAD、ICGC-PACA-AU Seq和PDX

2.实验流程

1） 低等级独立成分分析（LRICA）和预后相关成分提取：应用了改良版的低秩独立成分分析（LRICA）；PCA被用来将综合矩阵分解成低级矩阵（L）和稀疏矩阵（S），分别捕捉原始数据集的构成动态和去噪；用ProDenICA R软件包中的ProDenICA算法将L矩阵去卷积成独立的因素；ICA去卷积的结果是W和S矩阵，分别代表样本的正态矩阵和每个成分的基因贡献矩阵；在S矩阵上对每个成分的生物相关性和方向性进行加权，分别使用e1071和fgsea R包计算峰度过剩和基因集富集分析（GSEA），选择具有高峰度和富集在一组与生物相关的途径中的成分；对W矩阵的KM分析，按照预后标准进行了进一步的成分区分；应用主成分分析（PCA）确定预后估计的成分权重，然后在PCA样本维度坐标上采用Cox比例危险模型

2） 从大量RNA推断基质细胞类型的丰度及其与预后的关系：应用MCP-count算法来估计样本的免疫浸润和转录组数据的基质细胞丰度；特定的淋巴细胞T和B标记物的子集被用来确认细胞亚型的特异性

3） 主调节器（MR）梯度的生成和验证：将每个预后相关成分的转录因子（TFs）的核心作为表型驱动因素，将其贡献加权到全局分数中；MR-梯度的估计和可靠性在Puleo队列和ICGC-PACA-AU阵列中得到验证

4） 与MR-梯度分层相关的病人来源的异种移植（PDX）的分析：使用Limma R软件包对PDX的上皮肿瘤细胞（人类）和基质矩阵进行差异表达分析，比较特定的正交基因以估计MR-梯度转录因子的富集情况；将PDX的基质表达矩阵（小鼠）与上皮肿瘤细胞矩阵（人类）相加，计算出MR-梯度模型；应用MineICA R软件包中的ICA JADE算法和spearman相关性来提取与MR-Gradient相关的PDX表型；将TCGA-PAAD/ICGC-PACA-AU Seq综合矩阵中的ICA成分与PDX矩阵相关联，以评估PDX模型的再现程度

5） DNA甲基化分析：CpGs基因集富集分析（GSEA）是使用missMethyl R软件包进行的

6） SUV39H1/2对PDPCC的抑制和RNA-seq分析、KAT2B siRNA 转染和 RNA-seq 分析

7） 功能分析：为了确定与所选ICA成分和差异表达分析有关的途径，使用fgsea R软件包进行了基因集富集分析（GSEA），该软件包在预先排名的基因列表和MsigDB信号数据库上进行GSEA

8） PDX的CNV分析、蛋白质提取和Western Blot、免疫荧光和核定位的定量分析、免疫荧光和信号量化

9） 脂质组学分析：共检测到28个亚家族，每个亚家族的每个样品的中位数被计算出来，选择与MR-梯度模型有显著统计相关性的脂质亚家族；用主成分分析法（PCA）对所选亚家族中与MR-梯度相关度较高的代谢物进行分析，通过PCA维度坐标应用于Cox比例风险模型，对代谢物对PDAC预后的贡献进行加权

10）葡萄糖和谷氨酰胺的代谢：3名高危和3名低危患者身上提取的1mm2的PDX外植体

三、实验结果

01 - 通过转化细胞和微环境成分的解卷转录组分层提高了PDAC的预后预测能力

        来自三个RNA表达数据集的数据，产生了一个可以捕捉PDAC异质性的发现队列，本研究实验数据来自90个病人衍生的异种移植（PDXs）和两个来自公共数据集（TCGA-PAAD和ICGC-PACA-AU Seq）。作者将低秩ICA（LRICA）应用于发现队列的分析，还使用稳健的主成分分析（RPCA）将表达矩阵分解为低秩（L）矩阵和稀疏（S）矩阵（图1a）。这种方法能够分离肿瘤细胞群的底层生物学（L）和噪声（S），而不改变表达矩阵的整体结构和数据集之间的接近性（图1b）。随后应用了ProDenICA，目的是增加可以从L矩阵中获得的信息。在这种方法中，每个成分的选择和方向性是由超额峰度（评估正常值分布）和GSEA测量的生物相关性决定的，这导致了五个成分（图1c）。支持这种方法的验证，作者发现了一个被称为PAMG的成分，它捕捉到了肿瘤上皮细胞的表型，显示出对原代细胞和鳞状表型的强烈极化（图1c）。此外，该方法区分了两个微环境特异性成分，包括一个主要由成纤维细胞相关特征定义的成分，如MYCAF和ICAF（补充图1a）。第二个是免疫学成分，它捕获了造血系和炎症过程（图1c）的信息。最后，还确定了一个神经分泌物和一个细胞周期成分。因此，这种方法产生了更合适的加权信息，即在数据的数学结构中存在的生物信息类型，并有助于胰腺癌相关过程的基于生物信息学的建模。

        接下来，作者确定了上述组件与患者总生存期（OS）的关联，使用TCGA-PAAD和ICGC-PACA-AU Seq作为下游分析的发现队列。排除了人类异种移植数据集以避免与缺乏微环境区间有关的任何偏差。KM生存分析显示，患者的OS与PAMG、基质、免疫学和细胞周期成分高度相关，不良预后与细胞周期和EMT途径的激活相关（补充图1b）。另一方面，预后良好的患者与脂质代谢和免疫学途径的富集高度相关（补充图1b）。免疫学相关途径决定了基质成分内的预后情况，强调了将微环境成分作为预后标记的重要性。

        为了权衡KM分析的重要组成部分对预后的贡献，在PCA坐标上应用了PCA和Cox比例-hazards模型。该分析表明，维度3（Dim 3）和维度4（Dim 4）都与OS明显相关，其中Dim 3的贡献最强（图1c，补充图1c）。这一观察结果被多变量cox回归分析所证实，Dim 3仍与OS明显相关。维度3主要由免疫学成分代表，相关系数为0.68（图1d）。此外，对基质区间的进一步反卷积表明，有利的预后与T细胞和B细胞呈现正相关，而与成纤维细胞的丰度呈现负相关（图1d; 补充图1d）。具体来说，细胞毒性细胞标志物CD8A和两个浆细胞标志物CD27和CD38与Dim 3成正相关。综上所述， PDAC的预后不是由肿瘤细胞或微环境单独决定的，而是通过它们的综合贡献决定的。因此，通过强调这种相互依存关系，应有助于更好地概念化地寻找和开发标记物和靶向药物。最后，发现根据PAMG将PDAC分层为一个分子梯度，当免疫学成分的贡献加强时，可以明显区分患者的临床结果（图1e）。

图1 PDAC生物相关成分测定 补充图1 选定成分的特征

02 - 转录调控网络分析提供了病理生物学信息，并为病人分层产生了有用的分子标志物

        作者对上述两个公共数据集上进行了转录调控网络分析，以揭示支配每个有助于预后预测的成分的关键上游调节器。将分析的重点放在对LRICA成分有高贡献的相关转录因子（TFs）上，即PAMG以及免疫和基质成分，这些转录因子显示出与OS的显著关联。首先，使用ARACNe算法为每个成分构建了一个调控转录网络（RTN），用于通过GO注释识别的TFs。使用这种方法，共检测到113个TFs，代表121个调节子。随后，通过测试与特定成分相关的每个调节子的富集程度，确定了TF集内的主调节器（MR），检测到54个绝对富集分数>1的MR。PAMG显示了一个紧凑的相互作用网络（图2a），主要由祖先相关的调节子驱动，其中HNF4A、NR1I2和GATA6的贡献最大。鳞状MR网络与SNAI2、MYBL1和HMGA2有关。此外，观察到免疫学MR网络极化为调节性和促炎性节点（图2a），其特点是与Treg细胞（FOXP3和STAT5）和B/T细胞激活（IKZF1和NFATC2）相关的TFs。最后，基质成分的特点是代表多种微环境细胞类型的多态性TFs。然而，观察到免疫学相关调节因子的富集，如MAFB、BCL6B、IKZF3和SP1（图2a）。

        一旦建立了每个预后相关成分的转录调控网络，作者假设MR可以准确推断病人的预后，以无偏见的方式捕捉细胞的整体表型。为了测试这个想法的有效性，对每个MR应用了Cox单变量比例风险模型来评估它们的预测能力。发现HMGA2、SNAI2、GATA6和ZFPM1与预后的关联度最高，与评估所用的队列无关（补充图2）。为了生成捕获上皮和微环境特征的一致分层，用从转录网络分析中提取的MR建立了一个统一的梯度，每个梯度的计算，用每个MR（i）和病人（j）的富集分数（ES）来加权基因表达（GE），然后按比例求和。该分析显示，来自PAMG和免疫学转录因子的综合贡献优于其他成分，即使合并起来，也能估计发现队列（补充图2a），ICGC-Array（补充图2b）和Puleo（补充图2c）中患者的预后。作者将这种新的连续分层方法称为MR-梯度。值得注意的是，这种MR-梯度法简化了预后估计，它使用了一组40个MR，同时捕捉了上皮和微环境特征。

补充图2MR梯度在PDAC人类队列中的验证

        然后，对发现的队列实施了ICA JADE算法，以揭示与MR-梯度相关的全局表型。ICA2显示出与MR-梯度的正相关，显示出CHOLESTEROL_HOMEOSTASIS和FATTY_ACID_METABOLISM途径的预后良好，而不利的表型与HYPOXIA、EMT和CELL_CYCLE途径的上调相关联(图2c)。最后，评估了这个PDAC预后价值，以捕捉先前建立的亚型。发现ICA2包含了在肿瘤细胞和微环境水平上决定患者预后的关键特征，分别代表了原基/鳞状谱和激活的基质（图2d）。这些结果共同表明，基于MR的精细PDAC梯度是一个强大的临床可操作工具，可用于患者分层。

图2 转录主调节器（MR）的梯度和肿瘤表型的特征

03 - 患者来源的异种移植（PDXs）再现了PDAC预后相关的关键特征

        为了在多组学水平上扩展PDAC特征，作者使用90个PDXs，这些PDX已经证明了它们作为推导具有重要医学相关性的分子特征的工具的实用性。由于PDX队列经历了数据驱动的解构，以明确代表上皮细胞（人类）和微环境（小鼠），这可能是由宿主贡献的，作者进行了差异表达分析，以验证表征这两个隔间的MR。在人类肿瘤细胞和小鼠微环境中，代表PAMG和免疫学区间的MR分别显示出较高的表达水平（图3a）。GSEA显示代谢和上皮细胞分化途径在PDXs的人类部分富集，以及造血和ECM途径在小鼠基质区块的上调（图3a）。

        随后，重组了PDX表达矩阵，加入了人类和小鼠区间的表达矩阵，以估计MR-梯度，从而权衡其对患者结果的预测能力。单变量Cox回归分析显示，PDX衍生的MR-Gradient与患者OS之间存在明显的正相关关系（图3b）。此外，测量了PDX和发现队列之间的转录组学的相似性，以验证其作为定义预后特征的可靠代表的用途。使用ICA JADE算法将PDX混合矩阵解读为潜在的生物空间，确定ICA2与MR-梯度高度相关（图3b）。该成分显示了在发现队列中观察到的相同的表型极化，捕捉到原基和鳞状的特征，加上微环境衍生的ICAF和活化基质（图3c）的预后不良特征。此外，来自发现队列的ICA2和来自PDX组的ICA2在转化细胞和基质水平（图3d），保持了PDAC关键细胞特征的基因贡献的方向性，例如细胞骨架蛋白和代谢转运体等等（图3e）。

图3 患者衍生的异种移植再现了PDAC预后的关键决定因素

04 - CpG甲基化谱有助于定义PDAC的转录组表型

        作者首先分析了MR的DNA甲基化水平。大体上，决定祖先表型的关键MR显示出强烈的高甲基化，并与不良预后有关，特别是ZFPM1、GATA6和HNF4A（图4a）。此外，还进行了ICA，以捕捉与患者预后相关的甲基组图谱。共有12162个重要的（SD≥3）CpG被选入组件，其甲基化程度被分析为每个病人的β值的中位数。与祖先相关的TR一样，观察到DNA甲基化水平的增加与MR-Gradient的减少有关（图4b）。这一成分富集了与脂质代谢途径有关的CpGs，包括GLYCEROPHOSPHOLIPID METABOLISM和FATTY ACID TRIACYLGLYCEROLISM。这些结果在TCGA-PAAD队列上得到证实，其中23,448个CpGs在所选成分中显示出高贡献度（图4d-f）。因此，MR的DNA甲基化状态作为一种潜在的表观遗传机制，有助于PDAC预后相关的表型。

图4 CpGs甲基化概况与MR-梯度的关系

05 - SUV39H1 / 2和KAT2B是两种基于组蛋白的拮抗途径，有助于建立PDAC转录组学谱

        作者探索了多个组蛋白修饰因子和读取物（作为额外的表观遗传调节因子）与我们的表型类别的相关性，以调节转录结果的机制，发现149个蛋白质与MR-梯度有明显的相关性。发现SUV39H1、SUV39H2和KAT2B（图5a），它们是对H3K9残基有明显拮抗作用的writer，即甲基化的抑制（SUV39H1/2）与乙酰化的激活（KAT2B）。值得注意的是，与SUV39H1/2相反，发现KAT2B的基因组缺失发生在25%的队列中（图5b），这一数据被TCGA证实（图5c）。这种染色体缺失显示出与KAT2B启动子甲基化的平衡，调节其表达，从而调节预后表型。此外，KAT2B下调是鳞状表型的一个重要特征（图5c）。因此，作者在一组代表MR-梯度极值的PDX样本中量化了特定的表观遗传标记，即SUV39H1和SUV39H2的H3K9me3和KAT2B的H3K9ac。此外，用两个著名的激活标记H3K4me3和H3K27ac作为参考来补充组蛋白标记分析，这两个标记分别显示了与PDAC表型相关的一般和极化表达模式。H3K9me3和H3K9ac显示出相反的模式，其中三甲基化标记在高危患者的表观遗传学景观中占主导地位，大约有60%的阳性核（图5d）。相反，K9乙酰化在结果良好的群体中普遍存在，同时还有高水平的H3K27ac。用H3K4me3染色作为对照标记，因为它显示了独立于表型的同质表达水平（图5d）。

图5 鉴定与PDAC预后表型相关的关键表观遗传修饰因子

        功能验证是通过分别抑制或耗尽SUV39H1/2或KAT2B进行的。用10nM的Chaetocin，对6个PDX衍生的原始细胞培养（PDPCC），降低了H3K9me3水平，并伴随着祖细胞相关基因的上调（图6a）。通过使用一组特定的KAT2B siRNA，H3K9ac的水平被降低（图6b）。KAT2B的下调导致鳞状表型，这在具有中性CNV的PDPCC中更为明显。

图6 SUV39H1/H2和KAT2B的调控决定了PDX衍生的原始细胞培养（PDPCC）的表型

06 - 将MR-Gradient与支撑预后相关表型的代谢组功能联系起来

        作者分析了代谢网络在我们的MR-梯度的背景下是如何表现的，以及它与预后的关系。利用PDX的表达和中心代谢途径中关键酶的编码基因的甲基组图谱建立了一个代谢图（补充图4）。发现诸如ACSS1、ACACB和HMGCR等分别参与醋酸代谢、FA合成和胆固醇代谢的酶与MR-梯度呈正相关。这种亲脂代谢反映在与甘油磷脂和鞘脂途径的复杂脂质生物合成有关的酶的高表达上。相反，不利的表型增强了以Warburg效应和谷氨酰胺分解为中心的OXPHOS独立代谢，其中氨基酸的合成和甘油三酯（TG）在脂滴中的积累起到了核心作用。

补充图4 MR-梯度后的PDAC代谢图

        这些观察结果通过72个PDXs的脂质组学分析得到证实，其中检测到28个亚家族；然而，只有6个显示出与MR-Gradient的显著关联（图7a）。具体来说，磷酸甘油酯与MR-梯度显示正相关，而TG和神经酰胺代谢物，如单己基甘油酰胺和鞘氨醇，与预后不良有关。对所选亚家族高度相关的代谢物进行PCA后的Cox回归分析（图7b）显示TG和神经酰胺代谢物与预后有很大关系，分别占维度1贡献的53.84%和24.6%。此外，在PDXs中评估了SPHK1和PLIN2的表达，以验证这些蛋白分别作为与鞘脂代谢和TG积累相关的不良预后的标志物。值得注意的是，SPHK1和PLIN2在高风险（76%-85%）中显示出更多的阳性细胞，而它们在低风险（2%-12%）患者中的代表性也与EMT标志物vimentin正相关（图7c）。最后，通过测量PDX外植体上清液中的相应代谢物，证明高危患者表现出对无氧糖酵解和谷氨酰胺酵解的依赖。高危患者与低危组相比，葡萄糖的消耗和乳酸的产生分别高出1.4和2.6倍（图7d）。一致的是，高风险样本中谷氨酰胺的消耗和谷氨酸的产生分别增加了0.4和0.7倍（图7d）。这些结果强调ATP来源和脂质代谢是预后相关表型的决定因素。

图7 MR-Gradient捕获与PDAC预后相关的关键代谢特征

四、结论

        目前的研究为通过转录网络、DNA甲基化、表观基因组调节器和代谢组学对患者预后进行有价值的预测提供了强有力的整合，这些机制具有预后和机理价值，并揭示了对抗这种疾病的潜在治疗目标。