单细胞技术评估胰腺导管腺癌转移风险

小编今天给大家带来的是2022年11月发表在International Journal of Molecular Sciences（IF=6.2） 杂志上的一篇利用单细胞技术开发评估导管腺癌转移风险的方法。(Estimating Metastatic Risk of Pancreatic Ductal Adenocarcinoma at Single-Cell Resolution)。

摘要：

胰腺导管腺癌(PDAC)具有瘤内异质性，患者往往在发生转移后才被确诊。因此，研究如何有效地评估胰腺导管腺癌的转移风险具有重要意义。本研究中，我们提出了scMetR方法来评估胰腺导管肿瘤细胞的转移风险，该方法主要基于单细胞RNA测序进行展开。首先，我们确定了多种细胞类型，包括肿瘤细胞和其他细胞类型。然后，我们根据scMetR评分将肿瘤细胞分为三个亚群，包括转移性肿瘤细胞(MFTC)，移行性转移肿瘤细胞(TransMTC)和常规肿瘤细胞(ConvTC)，通过比较转移性肿瘤细胞(MFTC)和常规肿瘤细胞(ConvTC)，我们确定了转移特征基因 (MSGs)。功能富集分析显示，上调的转移特征基因 (MSGs)富集于多个转移相关通路。我们还发现，高表达上调的转移特征基因 (MSGs)的患者预后更差。转移性肿瘤细胞(MFTC)空间定位显示，它们偏向定位于癌和导管上皮区域，同时这些区域富集了导管细胞相关的炎症。我们进一步分析了细胞间的相互作用，观察到转移性肿瘤细胞(MFTC)中与转移相关的ADGRE5信号通路之间的细胞作用比其他肿瘤亚群中增加。最后，我们预测了12种具有逆转转移特征基因 (MSGs)表达潜力的候选药物。综上所述，我们借助单细胞技术开发了scMetR方法，该方法可能有助于分析肿瘤的异质性。

介绍

胰腺癌(PC)是最危险的癌症之一，其5年生存率低于5%。大约90%的胰腺癌患者是胰腺导管腺癌(PDAC)。胰腺导管腺癌具有较高转移扩散能力，这往往会导致不良预后。手术切除仍然是唯一可能治愈的治疗方法。当胰腺导管癌发生远端迁移的时候。

肿瘤细胞远端转移是一个多步骤的过程，包括局部侵袭、体内转移、血液循环中存活、外渗、适应新环境生存，最终在远处器官定植和生长。上皮间充质转化(EMT)在肿瘤转移中发挥重要作用。在皮间充质转化过程中，肿瘤细胞的细胞黏附分子的表达逐渐减少，并获得间充质表型以及迁移和侵袭能力，能够使癌细胞脱离原发肿瘤并在其他地方形成转移瘤。

在远端器官中具有明显的间充质-上皮转化(MET)过程。间充质-上皮转化(MET)可使肿瘤细胞恢复上皮形态，有利于胰腺导管腺癌(PDAC)细胞侵袭周围组织。因此，原发组织和远处器官中可能存在具有转移特征的细胞。以往研究报道肿瘤细胞分化状态的差异也与转移相关。不成熟的肿瘤比分化程度更高的肿瘤更具侵袭性。低分化肿瘤的扩散速度比高分化肿瘤快。胰腺导管腺癌(PDAC)具有不同的瘤内异质性，因此评估其转移风险具有一定的挑战性。

单细胞RNA测序(scRNA-seq)是可以区分肿瘤中不同的细胞类型。一些研究基于临床、组织学和基因组特征来评估肿瘤的转移风险。然而，这些研究不能应用于单细胞转录组学。

我们对之前使用原发性和转移性胰腺导管腺癌(PDAC)的scRNA-seq数据进行分析，我们确定了胰腺导管腺癌(PDAC)中的主要细胞类型，并剖析了胰腺导管腺癌(PDAC)的异质性，建立了用于评估肿瘤细胞的转移风险的scMetR方法。研究过程中我们发现肿瘤细胞中异质性的亚群具有不同转移风险。我们进一步鉴定了转移特征基因 (MSGs)，并对其从功能富集、患者预后、转移相关基因的空间表达模式、肿瘤亚群的空间位置和细胞间通讯等多个角度进行验证，探讨肿瘤亚群的转移风险。上述发现在一个独立的scRNA-seq数据集中得到重现。最后，通过基于转移特征基因 (MSGs)的药物预测，我们发现了一些转移的候选抑制剂。

图1. 胰腺导管腺癌(PDAC)数据分析流程。

结果

胰腺导管腺癌原发和转移患者的主要细胞类型

我们使用10个原发胰腺导管腺癌患者和转移胰腺导管腺癌患者的单细胞数据进行分析。质量控制之后，我们获得了单细胞转录组，其中原发癌症样本一共有8201个细胞，转移样本一共有7332个细胞。我们对所有细胞的高变异基因(HVGs)进行了主成分分析(PCA)。采用T-分布随机邻近嵌入(t-SNE)技术对细胞进行可视化。我们发现细胞主要根据样本分组，这提示样本之间的批次效应存在于细胞中(图2A)，并且会混淆PDAC的关键生物学变异。因此，我们整合所有的样本并进行批次矫正。然后，我们检查了整合数据的样本分布，发现大多数细胞的分布是无偏的(图2A)。使用Seurat (V3.2.2)进行的无监督聚类显示了7个簇，这些簇被注释为肿瘤细胞、癌症相关成纤维细胞(CAF)、T细胞和巨噬细胞(图2B)。为了确定每个细胞簇的身份，我们通过进行差异基因表达分析产生了细胞簇特异性的标记基因。我们使用了一些知名的标记基因来鉴定细胞类型，如EPCAM和KRT19代表肿瘤细胞，COL1A1和ACTA2代表癌症相关成纤维细胞(CAF), CD3D和CD3E代表T细胞，CD68和CD14代表巨噬细胞(图2C,D)。当这些标记基因在对应的细胞中高表达时，说明注释的标记基因是正确的。

为了进一步验证注释基因的准确度，我们使用inferCNV工具对注释基因进行基因拷贝数分析。使用Genome Sequence Archive的正常上皮细胞表达谱作为参考基因，分析结果显示肿瘤细胞比正常上皮细胞具有较多的拷贝数变异(copy number variation, CNV)(图2E)。这与胰腺导管腺癌(PDAC)具有大量拷贝数变异(copy number variations, CNA)的结果一致。

图2。胰腺导管腺癌(PDAC)的细胞组成(A)不同患者的t分布随机邻域嵌入(t-SNE)图。细胞根据患者进行着色,整合前(左)和整合后(右)。(B)原发性和转移性胰腺导管腺癌(PDAC)中4种主要细胞类型的t-SNE图。(C)各细胞类型特异性标记的热图。(D)小提琴表达图。y轴表示表达水平。(E)推断的正常上皮细胞和肿瘤细胞的拷贝数变异(CNV)热图。

2.2. 肿瘤细胞的异质性

目前认为，胰腺导管腺癌(PDAC)可以由多种肿瘤细胞亚群组成，这些肿瘤细胞在许多特性(如转移风险)方面存在差异。为了探索肿瘤细胞的异质性，我们利用基于平均剪影宽度(silwidth)的分辨率将肿瘤细胞重新聚类成亚簇。silwidth代表一个簇的相似性和不同簇的差异性。silwidth值越大，说明同一簇的细胞具有较高的相似性。silwidth值越大，说明同一簇的细胞具有较高的相似性。我们通过0.1的步骤计算了0.1到1不同分辨率下的平均silwidth宽度，观察到silwidth宽度的平均值在0.1时最大(图3A)。最后，在0.1的分辨率下将肿瘤细胞划分为6个亚簇(图3B)。为了估计每个亚簇的转移风险，我们使用S_EMT和S_CyTo方法进行评估，其中S_EMT是评估间充质-上皮转化(MET)相关基因表达的评分，S_CyTo是评估肿瘤细胞差异潜能的评分。之前的一项研究发现，分化良好的肿瘤细胞与正常细胞更相似，与分化较差或未分化的肿瘤细胞相比，倾向于生长和扩散更慢。此外，间充质-上皮转化(MET)也与肿瘤转移相关。为了准确评估肿瘤细胞亚群的转移风险，我们通过整合S_EMT和S_CyTo，提出了转移风险评分scMetR(图3C)。结果显示，亚群0、1和2的scMetR评分有显著差异(Wilcoxon检验:p < 0.01;图3D)。因此，我们确定了三个具有不同转移风险的亚群。我们将亚群0和2定义为具有转移特征的肿瘤细胞(MFTC)，亚群3、4和5定义为转移性肿瘤细胞(TransMTC)，亚群4和5定义为转移性肿瘤细胞(TransMTC)。亚簇1作为传统肿瘤细胞(ConvTC)。我们比较了不同肿瘤亚群的scMetR评分，发现具有转移特征的肿瘤细胞(MFTC)的scMetR评分显著高于移行转移性肿瘤细胞(TransMTC), 移行转移性肿瘤细胞(TransMTC)的scMetR评分显著高于传统肿瘤细胞(ConvTC)(图3E)。然后，我们可视化了原发和转移样本中肿瘤亚群的数量和比例。虽然具有转移特征的肿瘤细胞(MFTC)在转移性样本中较多，但其比例在原发样本中较高(图3F和G)。

图3。不同转移风险的肿瘤亚群。(A)每个分辨率下的平均剪影宽度(silwidth)图。(B)分辨率为0.1的肿瘤亚簇的t-分布随机邻域嵌入(t-SNE)图。(C)计算scMetR评分的示意图。(D)每个肿瘤亚簇的scMetR评分的小提琴图。(E)每个肿瘤亚群的scMetR评分的小提琴图。(F)原发灶和转移灶肿瘤亚群的t-SNE图样本，基于scMetR进行识别。右上:主要样本。右下:转移性样本。(G)肿瘤原发和转移患者比例条形图亚类。使用Wilcoxon检验*** p < 0.001

2.3. 与转移相关的特征基因

为了揭示转移在癌症中的作用，我们进行了进一步的分析。我们发现了具有转移特征的肿瘤细胞(MFTC)和传统肿瘤细胞(ConvTC)之间的差异表达基因(differentially expressed genes, DEGs)，并将其定义为转移特征基因（MSGs），这些基因一共有431个。为了了解潜在的相关功能，我们对上调的转移特征基因（MSGs）进行了功能富集分析，并确定了几种与转移相关的功能，包括细胞趋化、缺氧反应、粒细胞趋化和髓样白细胞迁移(图4A)。这表示上调的转移特征基因（MSGs）可能促进肿瘤转移。研究发现，转移性胰腺导管腺癌(PDAC)患者的临终结果比原发胰腺导管腺癌(PDAC)患者差，表明高转移风险患者的生存比低转移风险患者的预后要差。因此，我们使用了癌症基因组(the Cancer Genome)中患者的生存数据，通过TCGA和ICGC数据库分析患者预后，验证具有转移特征的肿瘤细胞(MFTC)是否具有较高的转移风险。研究发现转移特征基因（MSGs）高表达的患者可能具有较高的转移风险。因此，我们进行了单样本基因集富集分析(ssGSEA)来评估上调的转移特征基因（MSGs）的表达水平，并根据ssGSEA评分将患者分为高组和低组。生存分析结果表明，高表达上调MSGs的患者生存期更短(图4B,C)。进一步，我们对已发表的2例原发性具有转移特征的肿瘤细胞(MFTC)患者的数据进行了空间转录组学(ST)分析，以探索具有转移特征的肿瘤细胞(MFTC)的空间特征。利用PCA对空间转录组学(ST)数据进行高变异基因(HVGs)分区。进行无监督聚类后，我们发现识别出的聚类与定义明确的组织学注释一致(图4D)。我们观察到一些上调的转移特征基因（MSGs）,这些基因在之前研究中被证实可以促进胰腺导管腺癌(PDAC)细胞的迁移，如TFF1, TFF2和HMGA1。它们在导管上皮区域和癌区高表达(图4E和4F)。通过绘制胰腺导管腺癌(PDAC)切片中的细胞类型，我们发现转移性肿瘤细胞(MFTC)优先定位于导管上皮和癌区域(图4G)。Moncada等人研究发现与炎症相关的导管细胞在导管上皮区域富集。结果表明，转移相关基因在导管上皮区域高表达，移性肿瘤细胞(MFTC)主要分布在导管上皮区域，这与之前关于炎症促进胰腺导管腺癌(PDAC)转移的报道一致。

图4。分析上调的转移特征基因(msg)。(A)最重要的20个基因本体论(GO)术语的条形图，按- log10p值递减排序。x轴显示对应GO项富集的基因数量。红色的项与转移相关。(B,C)癌症基因组图谱(TCGA)和国际癌症基因组联盟(ICGC)数据集中患者的Kaplan-Meier生存曲线。(D)胰腺导管腺癌(PDAC)空间转录组学(ST)的无监督聚类分析。颜色表示不同的区域。(E,F) TFF1和TFF2表达模式的空间解析热图。(G)转移特征肿瘤细胞比例的空间分辨热图(MFTC)。比例由颜色表示。

2.4. 胰腺导管腺癌(PDAC)细胞之间的相互作用

除了细胞内在的信息外，细胞间的信息交流也可能促进肿瘤的转移。我们使用CellChat (V1.5.0)研究所有主要细胞类型和肿瘤细胞亚群之间的信号相互作用。我们比较了在胰腺导管腺癌(PDAC)中通过CellChat推断肿瘤亚群总数和相互作用的强度。我们发现，相互作用的总数和强度在转移性肿瘤细胞(MFTC)中最高，在移行性转移肿瘤细胞(TransMTC)中次之(图5A和5B)。为了确定促进转移的信号通路，我们预测了41个信号通路中的206个显著的配体-受体相互作用，包括ADGRE5, COLLAGEN, SPP1和ICAM。先前研究表明，ADGRE5信号通路的激活与淋巴结介入、转移和血管侵犯相关。我们推断ADGRE5信号通路细胞-细胞网络发现，转移性肿瘤细胞(MFTC)肿瘤亚群之间的相互作用强度最高，移行性转移肿瘤细胞(TransMTC)次之(图5C,D)。这些结果表明，与其他肿瘤亚群相比，MFTC中ADGRE5信号通路的相互作用的总数和强度以及相互作用的强度均增强，提示转移性肿瘤细胞(MFTC)在肿瘤亚群中具有最高的转移风险，这与肿瘤亚群的scMetR评分一致。

图5。主要细胞类型和肿瘤亚群的细胞间通讯。(A)来自肿瘤亚群的交互作用总数。(B)不同细胞类型之间相互作用的数量。连接单元格的线的粗细表示交互的数量。圆圈的大小表示单元格的数量。(C) ADGRE5信号通路中任意两个细胞群之间相互作用的热图。右侧和上方的条形图分别代表总输出和输入信号强度。(D)推断的ADGRE5信令网络。左、右侧图分别显示了对肿瘤亚群和其他细胞类型的自分泌和旁分泌信号。实圆和开圆分别代表源和目标。圆圈的大小表示单元格的数量。连接细胞的线条的粗细表示相互作用的强度。

2.5. 潜在与转移相关的抑制剂

手术切除胰腺导管腺癌(PDAC)是唯一可能治愈的治疗方法，但胰腺导管腺癌(PDAC)术后远处常常发生复发。因此，找到与转移相关的抑制剂具有重要作用。为了寻找潜在的治疗药物，我们分析了来自代表药物治疗前后基因表达变化的综合网络细胞标签库(LINCS)的数据。使用基因集富集分析(GSEA)，我们检测了药物使用后上调或下调的基因中是否存在转移特征基因（MSGs）。如果在药物使用后下调的基因中存在上调的转移特征基因（MSGs），则表明该药物可能逆转上调的转移特征基因（MSGs）的表达，成为潜在的转移抑制因子(图6)。

图6。基于转移特征基因（MSGs）的药物预测(A)药物预测工作流程。(B,C)各药物调控的基因网络。显示上调(B)或下调(C) msg中调节前10%基因的药物。矩形表示药物，圆圈表示基因。每一个线和节点颜色表示一种药物。如果一个基因可以被一种以上的药物调节，这个节点就会被涂成红色。(D)描述dinaciclib在肿瘤生长和转移中对细胞周期蛋白依赖性激酶5 (CDK5)的调节。

材料与方法

3.1. 单细胞RNA测序数据处理

我们从基因表达综合(GEO)数据库(GSE154778)下载了已发表的胰腺导管腺癌(PDAC)数据，包括10例原发样本和6例转移样本。每例患者获得的细胞数量范围为:原发样本为585 ~ 1570个，转移样本为272 ~ 2905个。数据集包含9621个原代样本细胞和7465个转移样本细胞。在少于3个细胞中表达的基因被去除。我们筛选了具有40%以上线粒体基因计数、7500多个检测到的基因和100,000多个唯一分子标识符(UMI)的细胞。使用Seurat (V3.2.2) R软件包进行预处理。

3.2. 确定胰腺导管腺癌(PDAC)的主要细胞类型

使用带有默认参数的NormalizeData函数对所有通过质量控制的单元格进行合并和规范化。使用FindVariableGenes函数筛选出前2000个HVGs。然后，使用带有默认参数的ScaleData函数对这些数据进行缩放。使用Harmony软件进行批量校正。FindCluster函数用于查找集群。该函数返回了几个细胞簇，这些细胞簇被t-SNE修饰。利用FindMarkers功能识别集群特异性的标志基因和不同细胞类型的标志基因。调整后的p < 0.05和|logFC| > 0.25(FC, fold change)设置为阈值。根据经典标志物的表达鉴定细胞簇的特性:肿瘤细胞的EPCAM、KRT19和KRT7; 癌症相关成纤维细胞(CAF)的COL1A1和ACTA2;T细胞CD3D和CD3E;巨噬细胞的CD68、CD14和AIF1。我们使用inferCNV (inferCNV of the Trinity CTAT Project: https://github.com/broadinstitute/inferCNV，于2022年4月10日访问)推断拷贝数变异分析。以GSA来源的正常上皮细胞(CRA001160)作为参照细胞。对于推断cnv分析，使用以下参数:cut off = 0.1，和HMM_type = ' i6 '。因此，生成热图来说明每个染色体的相对表达强度。

3.3重聚集肿瘤细胞用于探索肿瘤细胞的异质性

我们对肿瘤细胞进行了重新聚类，以细化具有不同转移风险的亚群。在少于3个肿瘤细胞中表达的基因被切除。进一步分析包括对数据进行标准化和缩放、识别hvg、批次校正和重新聚类。使用的功能与细胞类型注释相同。通过0.1的步骤，我们以0.1到1不等的分辨率重新聚集了肿瘤细胞。在这里，我们通过silwidth来评估不同分辨率下的聚类结果，silwidth是一个衡量簇内细胞和不同簇间细胞相似性的指标。silwidth越大，表明同一簇的细胞具有较高的相似性。使用聚类计算每个分辨率下的平均横径。

3.4肿瘤细胞亚群的转移风险评估

为了估计肿瘤细胞亚群的转移风险，我们为每个细胞提出了一个转移风险评分，命名为scMetR。使用以下公式对每个细胞进行评分：

其中S_EMT是ssGSEA计算的分数，以估计来自EMT-associated gene resource (dbEMT 2.0)的基因表达;S_CyTo为CyToTRACE计算的细胞分化评分;r为S_EMT与S_CyTo的Pearson相关系数。S_EMT和S_CyTo被缩放到[0,1]。采用Wilcoxon检验计算每个肿瘤亚群与其余肿瘤亚群之间的显著性差异。以P < 0.01为阈值。根据scMetR评分，定义了3个不同转移风险的肿瘤细胞亚群，包括MFTC、TransMTC和ConvTC。

3.5功能和生存分析揭示MFTC的转移风险

使用FindMarkers函数鉴定MFTC和ConvTC之间差异表达基因的(DEGs)，并将该类基因被定义为转移特征基因（MSGs）。调整后的p < 0.05和|logFC| > 0.25被设置为阈值。利clusterProfiler (V4.2.1) R包对基因本体(GO)的生物过程(BP)类进行功能富集分析。选择最小计数≥1且调整后p < 0.05的项。我们从TCGA中获取了PDAC患者的批量RNA-seq数据(http://portal.gdc.cancer.gov/，于2022年4月8日访问)以及ICGC的AU和CA队列

(http://dcc.icgc.org/， 2022年4月8日访问)。对无生存数据的患者进行筛选。最终TCGA、ICGC AU和CA队列分别保留146、90和215个样本用于生存分析。我们使用ssGSEA评估患者上调的差异表达基因的(DEGs)表达水平。根据ssGSEA评分将患者分为高危组和低危组。使用survminer (V0.4.9) R包的surv_cutpoint函数确定cut-off值，采用log-rank检验计算p值。

3.6 分析空间转录组学数据以发现空间特征

我们从GEO数据库(GSE111672)中访问了PDAC的公共ST数据，该数据包含两个主要样本。PDAC-A样本包含428个空间点，PDAC-B样本包含224个空间点。在少于5个点表达的基因被去除。我们使用NormalizeData函数对数据进行归一化，RunPCA函数执行PCA, FindNeighbors和FindClusters对ST点进行聚类。根据组织学特征对每个聚类进行注释。为了揭示MFTC的分布，我们使用条件自回归反卷积(conditional autoregressive-based deconvolution, CARD)结合scRNA-seq数据。CARD需要包含不同细胞类型基因表达信息的scRNA-seq以及包含定位信息的ST数据。CARD通过非负因子分解框架进行反卷积，并输出跨空间位置的估计细胞类型组成，其中包括PDAC中所有细胞类型的单细胞转录谱和两个PDAC样本的ST数据。利用CARD (V1.0) R包的CARD_deconvolution函数计算每个空间位置的细胞类型比例。

3.7 揭示促进转移的细胞-细胞通讯模式

为了揭示主要细胞类型之间的通讯模式，我们使用CellChat (V1.5.0)进行细胞间通讯分析。使用标准化数据创建了一个CellChat对象，用于推断不同细胞类型之间的细胞通信。CellChat通过评估和整合CellChat管理的人类配体受体数据库的基因表达水平，计算了有生物学意义的通信模式的概率。我们遵循了官方的工作流程并应用了预处理步骤，包括函数identifyOverExpressedGenes、identifyOverExpressedInteractions和projectData。然后，我们基于computecommunprobb, computeCommunProbPathway和aggregateNet计算了细胞类型之间潜在配体-受体相互作用。所有函数都使用默认参数运行。

3.8 逆转转移特征基因（MSGs）表达的候选药物的预测

我们从LINCS (http://www.Linc-sproject.org/，2022年6月2日访问)下载了药物治疗后的转录应答数据集。对每种药物，转录反应使用FC表示治疗和未治疗的情况。药物的作用条件不同，如不同的胰腺癌细胞株、药物浓度和治疗次数。计算每个基因在治疗组和未治疗组之间的FC。然后，使用基于分层多数投票方案的原型排序表(Prototype ranking List, PRL)合并不同情况下使用相同药物的FC排序表。各排序表中上调或下调的基因分别位于合并排序表的顶部或底部。接下来，我们使用GSEA检测被药物上调或下调的基因中是否存在转移特征基因（MSGs），以评估其作为转移抑制剂的潜力。潜在转移抑制因子上调的基因中包含下调的转移特征基因（MSGs）。PRL通过GeneExpressionSignature (V1.38.0)包的rank merge功能实现。利用fgsea (V1.18.0) R包的fgsea功能进行GSEA分析。通过设置p < 0.05的阈值来选择候选药物。

4 讨论

胰腺导管腺癌(PDAC)是一种异质性疾病，在肿瘤细胞中具有不同的转移风险。先前研究表明转移会导致患者预后较差。因此，探索肿瘤的异质性和潜在的转移风险是非常必要的，这是改善胰腺导管腺癌(PDAC)预后的关键。在本研究中，我们基于scRNA-seq数据鉴定了胰腺导管腺癌(PDAC)中的不同细胞类型，包括肿瘤细胞，CAF, T细胞和巨噬细胞（图2B）。

总之该研究确定了原发性和转移性胰腺导管腺癌(PDAC)的主要细胞类型，我们提出了scMetR来揭示关于转移风险的肿瘤异质性。进一步分析发现，基于scMetR定义的MFTC具有较高的转移风险。同时预测潜在的肿瘤转移抑制剂。因此，我们的分析可能提供一种新的评估转移风险的方法，这将有助于转移的诊断和治疗。