m5C 相关 lncRNA 风险模型鉴定新的PDAC预后相关基因

Prognostic Risk Model and Tumor Immune Environment Modulation of m5C-Related LncRNAs in Pancreatic Ductal Adenocarcinoma

胰腺导管腺癌中与m5C相关lncRNAs的预后风险模型和肿瘤免疫环境调控

发表期刊：Front Immunol

发表日期：2021 Dec 8

DOI:  10.3389/fimmu.2021.800268

期刊相关信息

一、背景

        胰腺导管腺癌（PDAC）是消化系统常见的恶性肿瘤，其特征是恶性程度高，易复发和转移。近年来，PDAC的发病率逐年上升。由于缺乏早期诊断方法，大多数患者在晚期被诊断出来，失去了根治性切除的机会。针对特定肿瘤相关基因的个性化分子疗法和免疫检查点抑制剂的应用有效地改善了晚期癌症患者的预后效果。然而胰腺癌由于其特定的肿瘤微环境而具有高耐药性和转移率，目前的治疗尚未有效地促进PDAC患者的预后，仍然需要更多特异性生物标志物来开发更有效的诊断和治疗策略。因此，必须找到新的分子生物标志物来早期诊断PDAC。

        RNA m5C甲基化是由甲基转移酶复合物催化的，该复合物由三种甲基化相关的酶组成，包括甲基转移酶（"Writer"）、脱甲基酶（"Erase"）和m5C结合蛋白（"Reader"）。与m6A甲基化的功能相似，m5C甲基化主要通过调节RNA稳定性、影响转录效率和介导RNA定位来发挥生物效应。

        长链非编码RNA由相应的基因转录，具有与mRNA相似的结构，包括polyA尾巴和启动子结构。由于不同的剪接方法，在分化过程中形成了多个lncRNAs。一些研究发现，lncRNAs可以在表观遗传、转录和转录后水平上调控下游基因，包括基因沉默、组蛋白加工、转录调控、转录干扰和核转运，这与人类各种疾病的发生密切相关。

二、材料与方法

1.数据来源

从TCGA数据网站上共下载了182个患者数据，包括临床和转录组表达的原始数据；与Ensemble Genes ID相比，在TCGA数据集中发现了14086个lncRNAs

2.实验流程

1）数据收集和处理：在线生物信息学工具GEPIA被用来检测m5C相关基因在胰腺正常和肿瘤组织中的表达

2）LASSO分析：LASSO回归方法构建m5C调控基因相关的lncRNA预测风险模型

3）GSEA：根据m5C相关lncRNAs的风险得分中值将所有样本分为两组

4）实验：细胞培养、RNA分离和定量实时PCR

5）鉴定肿瘤浸润性免疫细胞：CIBERSORT、ESTIMATE R软件包

流程图

三、实验结果

01 - 鉴定PDAC患者中与m5C调节相关的LncRNAs

        根据已发表的文章筛选了m5C甲基化修饰相关的基因，共筛选出13个m5C调节因子，即YBX1、ALYREF、DNMT1、NSUN4、TRDMT1、TET2、NSUN7、NSUN6、NSUN5、NSUN3、NSUN2、DNMT3a和DNMT3b。作者总结了m5C调节基因，然后利用在线数据库GEPIA比较了179个正常胰腺组织和171个胰腺癌组织中m5C调节基因的表达。发现在m5C调节因子中，YBX1、ALYREF、DNMT1和NSUN4在PDAC患者组织中的表达明显较高，而其它基因没有明显的差异。此外，根据177个PDAC样本中m5C调节因子和lncRNAs的表达情况，通过相关性分析来鉴定m5C相关的lncRNAs，筛选出242个m5C相关的lncRNAs。

        作者还构建了m5C相关基因与其共同表达的lncRNAs之间的网络，以显示m5C相关lncRNA的共同表达关系（图1A）。根据筛选出的m5C相关lncRNAs，采用单变量Cox分析，结合患者的生存数据，进一步筛选出PDAC患者的预后性m5C相关lncRNA。结果表明，共筛选出17个PDAC预后lncRNAs，大部分m5C相关lncRNAs是保护性因素（HR<1），只有CASC8是PDAC的危险因素（HR>1）。此外，保护性的m5C相关lncRNAs的表达是正相关的，而危险因素CASC8与其他因素是负相关的（图1B, C）。作者试图根据TCGA数据库来阐明所筛选的预后性m5C相关lncRNAs在PDAC患者中的表达。热图显示，所有的lncRNAs在正常和肿瘤的胰腺组织中都有统计学差异，这表明m5C相关的lncRNAs可能在PDAC的进展中发挥关键作用（图1D，E）。

 图1    在PDAC患者中识别具有预后价值的m5C调节相关lncRNAs

02 - 建立8-m5C-LncRNA风险模型

        根据m5C相关lncRNAs的初筛结果，进一步随机进行LASSO回归构建PDAC的预后风险模型，结果显示8个m5C相关lncRNAs适合构建预后风险模型（图2A，B）。风险得分计算如下：风险得分 =（-0.780839063578865 * AC022098.1.表达）+（-0.220638925265728 * AL031775.1表达）+（-0.0579614996945241 * AC005332.6表达）+（-0。367271578600146 * AC096733.3表达）+（-0.0490022123517058 * AC025165.1表达）+（0.0600252490282948 * CASC8表达）+（-0.528438814531151 * AC009974.1表达）+（-0.113641932851614 * PAN3-AS1表达）。Sankey图显示了6个m5C调节mRNAs和8个被筛选的lncRNAs之间的关系，其中CASC8属于危险因素，其他属于保护因素（图2C）。此外，根据上述公式计算的风险中位数，将所有PDAC患者分为低风险或高风险组。主成分分析（PCA）和三维PCA显示，不同风险的患者被很好地分成了两个群组（图2D，E）。高危组患者的生存时间明显短于低危组的患者。通过Kaplan-Meier生存分析，发现低风险组和高风险组之间的OS时间有明显差异（图2F）。

图2    在PDAC队列中构建m5C相关的lncRNA风险模型

03 - m5C相关LncRNA与临床病理参数的关系

        作者进一步检测了所选8个lncRNAs对PDAC患者总生存时间的影响。m5C相关lncRNAs的OS曲线显示，风险性lncRNA高表达（CASC8）的患者生存时间较短，而保护性lncRNA高表达的患者生存时间更长（图3A）。热图显示，胰腺肿瘤的大小在高风险组和低风险组之间有明显差异。尽管如此，其他临床因素包括患者年龄、性别、肿瘤分期和等级都没有统计学差异（图3B）。进一步将这些临床指标逐一细分，并分析了各亚组的风险分数。KM生存曲线显示，在男性患者组、年轻患者（年龄小于65岁）、I-II期组、I-II级或III-IV级组、T I-II或TIII-IV组、N0或N I-III组以及无任何转移的患者中，高风险评分的患者的OS较短（图3C）。

图3    m5C相关lncRNAs与临床病理特征的关系

04 - m5C相关的LncRNA风险评分是PDAC的独立预后因素

        作者根据风险得分对患者进行排序。热图显示，8个m5C相关lncRNAs的表达在不同风险组的PDAC患者中存在很大差异。此外，散点图也揭示了PDAC患者的死亡率随着风险分数的增加而增加，风险分数低的患者表现出更长的生存时间（图4A）。接下来，进行单变量和多变量的Cox回归分析，以验证由m5C相关lncRNA风险模型计算的风险分数是否可以作为PDAC患者的独立预后因素。如图4B所示，单变量Cox回归分析表明，只有m5C相关lncRNA风险评分与OS明显正相关。多变量分析还显示，m5C相关lncRNAs预后风险评分与PDAC患者的OS有明显关系，可以作为一个独立的预后因素（图4C）。结果表明，m5C相关lncRNA风险模型是PDAC的独立预后因素，并利用了传统的临床病理指标，包括患者年龄、性别和肿瘤分期。此外，1年ROC曲线证明，m5C相关lncRNA风险模型的AUC值为0.716，优于传统的临床因素如年龄，性别，stage，AJCC阶段，T阶段，N阶段和M阶段（图4D）。此外，3年ROC曲线分析也表明，风险评分AUC值是PDAC患者的一个优越的预后因素（图4E）。

图4    m5C相关lncRNA风险模型作为PDAC独立预后因素的评估

05 - 模拟m5C相关LncRNA风险模型的预测值

        为了准确预测PDAC患者的总生存期，作者根据m5C相关lncRNAs的表达风险评分和临床病理特征（包括年龄、性别、等级、T期、M期和N期）构建了一个nomogram，以揭示1、3和5年的生存率（图5A）。使用校准曲线来比较实际和预测的1年、1.5年、2年和3年患者生存率的一致性。发现，实际和预测线在3年内几乎一致（图5B）。上述结果证明，通过m5C相关lncRNA预后风险评分生成的nomogram是可靠的。此外，将TCGA的所有PDAC患者按1：1的比例随机分为两个亚组（A组和B组），对m5C相关lncRNA风险模型进行了内部验证。对每个亚组的KM生存曲线和5年ROC曲线进行了研究。结果显示，A组中m5C相关lncRNA风险评分较高的患者OS较短，5年ROC曲线的AUC值为0.814（图5C，D）。B组的患者有类似的OS趋势，AUC值为0.903（图5E，F）。此外，散点图也显示了所筛选的lncRNAs在高风险或低风险组中的差异表达显著，并且在所有亚组中，患者的生存时间与风险评分呈正相关（补充图3）。上述结果表明，m5C相关lncRNA风险模型是PDAC患者的可靠预测因素。

图5    检测m5C相关lncRNA风险模式的预测值 补充图3 热图显示8个m5C相关lncRNAs在亚组中的差异表达

06  - 检测LnRNA的体外表达和功能富集分析

        来自TCGA数据库的m5C相关lncRNA表达谱显示，AC022098.1、AL031775.1、AC005332.6、AC096733.3、AC025165.1、AC009974.1和PAN3-AS1在PDAC中下调，而CASC8过度表达。为了检查lncRNAs在体外的表达水平，作者使用了三个胰腺癌细胞系和一个正常的胰腺导管细胞系来进行qRT-PCR实验。体外结果与TCGA的数据不完全一致（图6A），发现CASC8在胰腺癌细胞（Mia-PaCa-2、CFPAC-1和Panc-1细胞）中上调，而AC096733.3与HPNE细胞相比表达量较低，这与TCGA数据表达情况相同。AC096733.3、AC025165.1、PAN3-AS1和AC009974.1在至少两个癌细胞系中下调，可能在PDAC中发挥肿瘤抑制基因的作用。上述结果与TCGA的患者OS分析数据部分一致。然而，AC022098.1和AC005332.6在Mia-PaCa-2中被抑制，而在CFPAC-1中过表达，它们在不同细胞中的表达并不一致。此外，AL031775.1只在Panc-1细胞中被下调。因此，m5C相关的lncRNAs调控PDAC患者OS时间的内在机制仍需进一步探讨。

图6    体外验证m5C相关lncRNAs的表达水平和功能富集分析

        为了寻找m5C相关lncRNAs在低风险或高风险组可能涉及的信号通路，作者进行了GSEA发现MTORC1信号通路在高危组被激活；然而，MYOGENESIS和KRAS信号通路在低危组被激活（图6B，C）。GO富集度最高的五个生物过程是T细胞激活、钙离子平衡、激素运输、化学突触传递的调节和跨突触信号通路的调节（图6D-F）。此外，KEGG分析显示，一些免疫相关的途径被富集，如原发性免疫缺陷、T细胞受体信号通路、白细胞跨内皮迁移和Th1/Th2细胞分化（图6G，H）。

07 - 与m5C相关的LncRNAs和肿瘤浸润淋巴细胞之间的关系

        作者分析了高风险或低风险组的肿瘤微环境的差异。通过使用CIBERSORT工具，共筛选出22种肿瘤浸润性免疫细胞（图7A）。结果显示， naïve B细胞、CD8+T细胞、调节性T（Treg）细胞和静息NK细胞在低风险组表现出更高的表达量，而M0和M2巨噬细胞在高风险组有更高的表达量。M2巨噬细胞已被证明是胰腺癌的一个致癌因素。因此，m5C相关的lncRNAs可能促进M2巨噬细胞在PDAC中的极化或浸润。此外，进一步检测了m5C相关lncRNAs与免疫相关基因（LMTK3、LAG3、CD27、CD28、CD86和BTLA）之间的关系，结果表明，除了AC025165.1，大多数lncRNAs与免疫相关基因有统计学相关性（图7B）。基因表达分析也显示，上述免疫相关基因的表达在低风险或高风险组中是不同的（图7C）。此外，通过ESTIMATE R软件包进行肿瘤微环境评分，分析肿瘤环境中基质细胞和免疫细胞的比例。明显发现，m5C相关lncRNA风险评分低的PDAC组织具有较高的免疫评分、基质评分和ESTIMATE总体评分，这表明低风险组的PDAC具有较低的胰腺肿瘤比例（图7D）。最后，作者还研究了m5C相关lncRNA风险评分与肿瘤淋巴细胞之间的关联。Spearman相关分析显示，风险评分与3种肿瘤浸润淋巴细胞（NK细胞、M0和M2巨噬细胞）呈正相关，但与调节性T细胞、CD8+T细胞、活化记忆CD4+T细胞、naïve B细胞和浆细胞呈负相关（图7E）。这些结果显示，m5C相关的风险评分可以区分PDAC中肿瘤浸润淋巴细胞的不同特征。

图7    m5C相关lncRNAs与肿瘤浸润淋巴细胞的关系

四、结论

        作者基于TCGA的转录组表达和临床数据，构建了PDAC患者的8-m5C相关lncRNA预后风险模型。该m5C相关lncRNA风险模型被证明具有独立的预后价值，并为PDAC患者提供了准确的生存预测。此外，本研究还为提供了对PDAC中肿瘤浸润淋巴细胞调控的更好理解。简而言之，m5C相关的lncRNA风险模型可以为我们提供潜在的PDAC生物标志物或治疗目标。