基于生物信息学分析肝细胞癌中lncRNA-miRNA-mRNA的调控网络

张 政，卢鸿健，李明慧，高文博，田 宇，王嘉祺，张荣花，章广玲

(1.华北理工大学临床医学院/河北省医工融合精准医疗重点实验室，河北 唐山 063000；2.华北理工大学基础医学院/河北省慢性疾病基础医学重点实验室，河北 唐山 063210)

肝细胞癌（hepatocellular carcinoma，HCC）是世界上常见的恶性肿瘤之一，占原发性肝癌的75%～85%。GLOBOCAN 的最新统计数据显示[1]，2020 年全球HCC 新发病例达91 万，死亡病例83 万，分别居肿瘤新发病率第6 位、死亡率第3 位；我国人口占世界的18%，但新发病例占全球的45%。HCC 的发生发展过程复杂，是一种异质性较高的癌症，大多数发现时已是中晚期，错过了手术的最佳时间且预后较差，其5 年总生存率只有12%[2]。因此，寻找HCC 发展过程中有效的靶点和信号通路并进行相应干预，对延缓HCC 的进展或者逆转具有重要意义。研究表明[3]，非编码RNA（non-coding RNA，ncRNA）包括长链非编码RNA（long non-coding RNA，lncRNA）和微小RNA（microRNA，miRNA）等，其可能在HCC 的发生发展中发挥了重要的作用。miRNAs 分子由18～25个bp 组成的非编码RNA，能够通过碱基互补配对与靶基因的3' 端非编码区（3'untranslated region，3'UTR）特异性结合，沉默或降解靶mRNAs 从而影响其翻译功能。lncRNAs 分子是一种长度超过200 bp 的非编码核苷酸，含有miRNAs 的结合位点，可以作为miRNAs 的“海绵”吸附miRNAs，从而竞争性抑制miRNAs 对靶mRNAs 的作用[4]。本研究拟通过生物信息学方法探究癌症基因组图谱（TCGA）数据库中HCC 组织中差异表达的基因，并构建lncRNAmiRNA-mRNA 调控网络，为靶向诊断和治疗HCC提供科学依据。

1 资料与方法

1.1 HCC 组织中差异表达基因的筛选 利用R 软件DESeq2 包分析TCGA 数据库中374 例HCC 组织与50 例正常肝组织中lncRNAs、miRNAs、mRNAs 的表达谱。使用R 软件ggplot2 包对表达谱数据绘制火山图，以|Log2FC|＞1（fold change 差异倍数大于2 倍）且P＜0.05 作为筛选标准。

1.2 miRNAs 转录因子的预测 利用FunRich 软件预测差异表达miRNAs 的转录因子，根据其所调控的基因数和值进行排序，选取富集程度最显著的前10 个转录因子进行展示。

1.3 基因功能（GO）和京都基因与基因组百科全书（KEGG）富集分析 利用R 软件clusterProfiler 包分别分析miRNAs、mRNAs 各基因之间潜在的影响，并以调整后的P＜0.05 为标准，进行GO 和KEGG分析，分别预测其主要生物功能和参与的主要信号通路。

1.4 miRNA-mRNA 调控网络 使用FunRich 软件预测miRNA 的候选靶基因，结合HCC 组织中差异表达的miRNAs 和mRNAs 的FC 值，筛选具有相互作用的miRNA-mRNA 关系。

1.5 PPI 网络构建 将筛选出的基因导入STRING 数据库，分析蛋白之间的相互作用（protein-protein interaction，PPI），并将分析结果导入Cytoscape 中，构建可视化PPI 蛋白互作网络图。

1.6 lncRNA-miRNA-mRNA 调控网络构建 利用Starbase 数据库，预测miRNAs 的上游lncRNAs，结合HCC 中差异表达lncRNAs 的FC 值，筛选具有相互作用的lncRNA-miRNA 关系对。基于共享的miRNAs 通过Cytoscape 软件对lncRNA-miRNAmRNA 调控网络进行可视化。

2 结果

2.1 差异表达lncRNAs、miRNAs、mRNAs 的筛选 利用R 软件DESeq2 包分析374 例HCC 组织和50 例正常肝 组织中lncRNAs、miRNAs、mRNAs 的表达谱，根据P＜0.05 和|Log2FC|＞1 的标准筛选筛选出HCC 组织中差异表达的lncRNAs 2850 个，其中上调2475 个，下调375 个；差异表达的miRNAs 38个，其中上调31 个，下调7 个；差异表达的mRNAs 4455 个，其中上调3199 个，下调1256 个，见图1。

图1 差异表达基因的火山图

2.2 转录因子-miRNA 调控网络 利用FunRich 软件对38 个差异表达miRNAs 进行转录因子预测，共得到201 个转录因子，根据其所调控的基因数和P值，选取富集程度显著的前10 个转录因子进行展示，见图2，其基本信息见表1，前10 个转录因子分别为EGR1、SP1、NKX6-1、POU2F1、LHX3、SP4、MEF2A、ARID3A、MSX2、HOXA9。

表1 差异表达前10 名转录因子的基本信息

图2 差异表达miRNAs 转录因子的预测

2.3 miRNAs 的GO 和KEGG 富集分析 对38 个差异表达的miRNAs 进行GO 和KEGG 功能富集分析，其生物过程主要涉及血管内皮细胞迁移的调节、血管内皮细胞增殖调控与新生血管生成、内皮细胞增殖的调节、血管内皮细胞增殖与新生血管生成；分子功能主要涉及mRNA 结合参与转录后基因沉默、RNA 聚合酶Ⅱ复合物结合、RNA 聚合酶核心酶结合、基础转录机制结合；KEGG 分析主要富集于miRNAs 在癌症中的作用，见图3。

图3 差异表达miRNAs 的GO 和KEGG 功能富集分析

2.4 miRNA-mRNA 调控网络 通过FunRich 软件预测差异表达的38 个miRNAs 的下游靶基因，共得到1530 个靶基因，然后与差异表达的4455 个mRNAs取交集获得286 个候选靶基因，见图4。通过Starbase 和MiRDB 数据库比对分析，获得14 302 对miRNA-mRNA，与差异表达的miRNA-mRNA 匹配后，最后筛选出交叉表达负调控的12 个miRNAs 和108 个mRNAs 作为后续研究对象。

图4 差异表达miRNAs 的下游靶基因与差异表达mRNAs 的韦恩图

2.5 108 个mRNAs 的蛋白PPI 调控网络及富集分析利用STRING 数据库和Cytoscape 软件构建筛选出的108 个基因的PPI 网络，见图5。图中节点代表mRNAs，有相互作用关系的mRNAs 之间有连线，共有44 个节点和134 条边；节点越大、颜色越深则介数中心度越大，边的粗细代表连接评分，边越粗，mRNAs 间的互作关系越大。

图5 108 个mRNAs 的蛋白互作网络图

2.6 108 个mRNAs 的GO 与KEGG 富集分析 利用GO 功能分析筛选出的108 个mRNAs 的生物功能，生物过程主要涉及肌肉适应、MAP 激酶活性的调节、胚胎器官发育、蛋白丝氨酸/苏氨酸激酶活性的正调控；分子功能主要涉及MAP 激酶酪氨酸/丝氨酸/苏氨酸磷酸酶活性、MAP 激酶磷酸酶活性、激酶调节活性、DNA 结合转录激活物活性和RNA 聚合酶Ⅱ特异性。KEGG 富集分析108 个mRNAs 参与的信号通路，主要涉及调节干细胞多能性的信号通路、催乳素信号通路、MAPK 信号通路、EGFR 酪氨酸激酶抑制剂耐药性，见图6。

图6 108 个mRNAs 的GO 和KEGG 富集分析图

2.7 构建lncRNA-miRNA-mRNA 调控网络 利用Starbase 数据库，分析筛选出的12 个miRNAs 的上游lncRNAs 共221 个，结合差异表达lncRNAs 的FC 值，筛选出5 对lncRNA-miRNA，包括2 个miRNAs（miR-23A 和miR-146B）与5 个lncRNAs（GAS5、MCM3AP-AS1、LINC00173、KCNQ1OT1、SNHG7）；结合筛选出的2 个miRNAs 的与组织中差异表达mRNAs 的FC 值，共匹配42 对miRNAmRNA（包括40 个mRNAs）。通过Cytoscape 软件对lncRNA-miRNA-mRNA 调控网络进行可视化处理，见图7。

图7 HCC 组织中lncRNA-miRNA-mRNA 的调控网络

3 讨论

HCC 是由肝细胞恶性病变引起的癌症，患者发现时大多数已是中晚期，具有发病隐匿、进展快、易耐药、预后差以及易复发等特点，且现有的治疗效果并不理想[5]。因此寻找开发有效的HCC 诊断和预后的生物标记物具有重要意义。越来越多的研究表明，lncRNAs 和miRNAs 等非编码RNAs 与HCC 的发病密切相关。LncRNA n339260 通过降低miR-30e-5p的表达使TP53INP1 的表达增加，从而促进肝细胞癌的进展[6]。LncRNA ZEB1-AS1 通过调节miR-1224-5p/MAP4K4 轴，调节持续沙眼衣原体感染中线粒体介导的HeLa 细胞凋亡[7]。抑制lncRNA GAS6-AS2 通过调节miR-136-5p/OXSR1 轴能够在体外和体内水平减轻败血症相关的急性肾损伤[8]。然而在HCC 中尚没有报道完整的lncRNA-miRNA-mRNA 调控网络，因此深入研究lncRNAs 和miRNAs 在HCC 发生发展中的作用及调控机制具有重要意义。

已有研究表明，转录因子调控miRNAs 的表达，且这些转录因子在HCC 的发病过程中具有重要的调控作用。本研究预测了调控HCC 中差异表达miRNAs 的上游转录因子共201 个，根据值筛选出差异显著的前10 个转录因子分别为EGR1、SP1、NKX6-1、POU2F1、LHX3、SP4、MEF2A、ARID3A、MSX2、HOXA9，其中EGR1 差异最明显。EGR1 是一种转录因子，可能参与肿瘤细胞的增殖、侵袭和转移以及肿瘤血管生成，与癌症的发生和发展密切相关[9]。有报道表明[10]，EGR1 可能以转录因子方式调控hsa-miR-146b 的表达。EGR1 在miR-195 的上游启动子区有结合位点，且调控miR-195 的表达，EGR1 可抑制miR-195/AKT3 轴并抑制胃癌细胞凋亡[11]。EGR1 还能够通过负调控miR-491-5p 促进细胞增殖、生长和迁移，并抑制细胞凋亡，从而抑制胃癌的进程[12]。因此，研究转录因子对miRNAs 的调控作用可在系统层次上对疾病的发生提供理论依据。

近年来，伴随着高通量测序技术的发展产生了大量的基因数据，这为进一步研究HCC 的发生发展及预后监测提供了有力的基础。本研究首先分析TCGA 数据库中差异表达基因的表达谱，得到的差异表达的lncRNAs 2850 个、miRNAs 38 个、mRNAs 4455 个。对38 个差异表达的miRNAs 进行GO 和KEGG 功能富集分析，显示主要涉及血管内皮细胞迁移和增殖、新生血管生成、内皮细胞增殖的调节等生物过程，mRNA 结合参与转录后基因沉默等分子功能，且KEGG 分析显示这些miRNAs 主要在癌症中发挥作用。通过FunRich 软件预测得到的miRNAs 靶基因与差异表达的mRNAs 交叉匹配，对筛选出的108 个mRNAs 进行GO 分析发现，主要在胚胎器官发育、MAP 激酶活性的调节等生物学过程富集，DNA 结合转录激活物活性和RNA 聚合酶Ⅱ特异性、激酶调节活性等分子功能富集。KEGG 分析发现，靶基因参与的信号通路主要涉及调节干细胞多能性的信号通路、MAPK 信号通路等。因此可以得出在HCC 中差异表达的mRNAs 其生物功能主要与细胞增殖、形态分化和细胞迁移有关。同样地，在PPI 网络的Top 10 基因中KLF4、NRG1、FOSB、MAP2K1 也参与了上述生物学功能以及相关信号通路。研究发现[13]，miR-146B 通过靶向KLF4 调节肝星状细胞活化。NRG1 是鳞状细胞癌分化的关键调节因子，抑制肿瘤细胞角化和终末鳞状分化并促进肿瘤细胞增殖[14]。TGF-β1通过诱导FOSB 促进前列腺癌细胞迁移和侵袭[15]。CircRNA ZFR 通过上调MAP2K1 促进HCC 的增殖[16]。

本研究寻找与miRNAs 相匹配的上游lncRNAs和下游mRNAs 构建HCC 中lncRNA-miRNA-mRNA 的调控网络，其中包含2 个关键的miRNAs（miR-23A、miR-146B）以及与其匹配的5 个lncRNAs（GAS5、MCM3AP-AS1、LINC00173、KCNQ1OT1、SNHG7）和40 个mRNAs（ROBO1、LIN28B 等）。研究表明[17]，miR-23A 是HCC 中下调最显著的miRNA，可能是改善肝癌患者索拉非尼反应性的潜在靶点。miR-146B 通过靶定PTP1B 抑制胃癌细胞的增殖并促进其凋亡[18]。miR-146B 通过NF-κB（核因子κB）途径靶向TLR4（Toll 样受体4），抑制胆囊癌肿瘤的发生和转移[19]。以上研究表明，miR-23A、miR-146B具有抑制肿瘤的作用。此外，lncRNA GAS5 靶向miR-23A 抑制CCl4 诱导的肝纤维化[20]。LncRNA SNHG7 通过调节miR-146B-5p/ROBO1 轴促进胰腺癌进展[21]。LncRNA SNHG17 通过抑制miR-23A-3p 来调节CXCL12 介导的血管生成，从而促进结直肠腺癌细胞的增殖和迁移[22]。LncRNA KCNQ1OT1通过miR-146A-5p/ACER3 轴抑制HCC 的放射敏感性并促进肿瘤发生[23]。

综上所述，lncRNA-miRNA-mRNA 调控网络的构建，有望阐明HCC 的发生发展机制，为HCC 的临床诊断和预后监测提供有价值的标记物。但是lncRNA-miRNA-mRNA 调控网络的功能以及各差异基因的功能，还有待通过体内外实验和临床实验进一步证实。