基于转录组测序构建食管腺癌预后相关的内源竞争RNA网络

陈超 李冬 童国俊 胡四平

食管腺癌已成为发病率最高的食管恶性肿瘤，也是生长速度最快、最具侵袭性的恶性肿瘤之一[1-2]。由于目前食管腺癌的预后仍极不理想，探索食管腺癌的发生、发展及预后影响因素是食管腺癌治疗的基础。微小非编码 RNA（microRNA，miRNA）可通过识别 mRNA上的miRNA应答元件抑制mRNA的翻译或导致mRNA降解，同时，长链非编码RNA（long non-coding RNA，lncRNA）上也存在miRNA应答元件并竞争性结合miRNA[3]。因此，lncRNA、miRNA、mRNA之间形成的内源竞争RNA（competing endogenous RNA，ceRNA）网络是食管癌预后机制研究的重要突破点，可能为食管腺癌患者的治疗提供新的方法。癌症和肿瘤基因图谱计划（The Cancer Genome Atlas，TCGA）数据库中包含了大量的肿瘤标本的数据，包括数百例食管恶性肿瘤患者的转录组测序数据及相对应的临床基本信息。然而，关于食管腺癌中的ceRNA网络及其影响预后的机制目前仍知之甚少。本研究提取食管腺癌组织的转录组数据，利用生物信息分析手段构建食管腺癌预后相关的ceRNA网络，以进一步探索食管腺癌预后影响机制，现将结果报道如下。

1 对象和方法

1.1 对象 选择截至2019年7月从TCGA获取78例食管腺癌组织和13例正常食管组织的lncRNA、miRNA、mRNA表达量数据，以及食管腺癌组织患者性别、年龄、病理分期等临床数据。食管腺癌患者男67例，女11例，年龄 68.4（58.0～77.1岁），生存期为 2.15（1.01～4.88）年；按肿瘤病理分期，处于Ⅰ期、Ⅳ期各10例（12.8%），Ⅱ期25例（32.1%），Ⅲ期 33例（42.3%），13例正常食管组织TCGA未提供临床信息。

1.2 方法

1.2.1 差异基因筛选 采用R软件分析食管腺癌组织和正常食管组织中存在差异表达的lncRNA、miRNA、mRNA。差异基因的筛选条件：（1）lncRNA、miRNA的差异倍数（fold change，FC）以2为底的对数的绝对值（|log2FC|）均为≥1；（2）mRNA 的|log2FC|≥0.6；（3）P＜0.05。以lncRNA、miRNA、mRNA的表达量的中位值将食管腺癌样本分为各基因所对应的高表达组和低表达组，并统计每个差异表达的lncRNA、miRNA、mRNA的基因数。

1.2.2 生存分析及权重基因共表达网络分析 用Kaplan-Meier法分析食管腺癌组织中各lncRNA、miRNA、mRNA的表达量的多少与生存时间的关系，若P＜0.05即说明该基因为食管腺癌预后相关基因。用R软件对差异表达且与预后相关的lncRNA和miRNA、miRNA和mRNA构建无尺度网络，以无尺度网络指数＞0.85为符合无尺度网络条件并确定软阈值。计算并获得lncRNA与miRNA、miRNA与mRNA的权重基因共表达网络。

1.2.3 预后相关的ceRNA网络构建及可视化 以miRNA为中心，在权重基因共表达网络中筛选与miRNA存在共表达关系的mRNA，同时逆向筛选可吸附miRNA并与预后相关的lncRNA，构建lncRNA/miRNA/mRNA内源竞争网络，并在Cytoscape软件中可视化；再次采用生存分析评估ceRNA网络中各基因的表达水平与预后的关系，分析ceRNA网络中各种RNA在食管腺癌组织、正常食管组织中的表达差异，以及miRNA与lncRNA、mRNA表达的相关性。

1.2.4 基因集富集分析（gene set enrichment analysis,GSEA） 将食管腺癌组织分别分为细胞周期蛋白激酶1（cyclin dependent kinase 1，CDK1）高表达组、CDK1 低表达组、微管成核因子（microtubule nucleation factor,TPX2）高表达组、TPX2低表达组，随后进行GSEA研究京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）结果。对富集分析结果根据标准化富集分数（normalized enrichment score，NES）进行排序，对各组中NES排名前5的KEGG通路结果进行可视化。

1.3 统计学处理 数据统计分析及可视化采用R软件（3.6.1版）、Cytoscape软件（3.7.2版）。计量资料如呈正态分布以表示，组间比较采用t检验，否则以四分位数表示，组间比较采用非参数检验。分类变量采用例数（百分率）表示，生存资料采用log-rank检验，相关性分析采用Pearson相关，P＜0.05为差异有统计学意义。

2 结果

2.1 差异基因筛选结果 与正常食管组织比较，食管腺癌组织中符合筛选条件被分到高表达组的lncRNA、miRNA、mRNA分别有678、131、1 862个，被分到低表达组的 lncRNA、miRNA、mRNA分别有 774、59、2 825个。

2.2 生存分析与权重基因共表达网络分析结果 （1）在所有表达有差异的lncRNA、miRNA、mRNA中，生存分析结果提示与预后相关的RNA数目分别有81、20、374个。（2）由这些与预后相关的 lncRNA、miRNA、mRNA组成的权重基因共表达网络中，lncRNA共有7个，分别为 AC131009.3、AC024337.2、AC245407.1、BDNF-AS、AL136172.1、AC087388.1、AC010261.2；miRNA 共有 3个，分别为 miR-125b-5p、miR-143-3p、miR-29c-3p，而与这3个miRNA又存在共表达关系的预后相关的mRNA有7个，分别为CDK1、TPX2、细丝蛋白A（filamin A，FLNA），肌肉相关受体酪氨酸激酶（muscle associated receptor tyrosine kinase，MUSK）、蛋白激酶camp依赖的Ⅱ型调节亚基β（protein kinase cAMP-dependent type Ⅱ regulatory subunit beta，PRKAR2B）、非受体5型蛋白酪氨酸磷酸酶（protein tyrosine phosphatase，non-receptor type 5，PTPN5）、血管活性肠肽受体 2（vasoactive intestinal peptide receptor 2，VIPR2）。

2.3 食管腺癌预后相关的ceRNA网络构建

2.3.1 预后相关的ceRNA网络 见图1。

由图1可见，在该ceRNA网络中AC131009.3、AC087388.1、miR-29c-3p、CDK1、TPX2 互相之间存在内源竞争网络关系，即构成了ACO87388.1/miR-29c-3p/CDK1、ACO87388.1/miR-29c-3p/TPX2、AC131009.3/miR-29c-3p/CDK1、AC131009.3/miR-29c-3p/TPX2 的ceRNA调控网络。

图1 食管腺癌的预后ceRNA网络图，菱形为lncRNA，三角形为miRNA，椭圆形为mRNA；深灰色为与正常食管组织比较，在食管腺癌组织中高表达，浅灰色为食管腺癌组织中低表达（TPX2为微管成核因子，CDK1为细胞周期蛋白激酶1，FLNA为细丝蛋白A，MUSK为肌肉相关受体酪氨酸激酶，PRKAR2B为蛋白激酶camp依赖的Ⅱ型调节亚基β，PTPN5为非受体5型蛋白酪氨酸磷酸酶，VIPR2为血管活性肠肽受体2）

2.3.2 ceRNA网络中各基因表达水平与生存时间的关系 见图2。

由图 2 可见，AC131009.3、AC087388.1、CDK1、TPX2在食管腺癌组织中的表达水平越高，患者的生存时间越短，生存率越低；miR-29c-3p在食管腺癌组织中的表达水平越高，患者的生存时间越长，生存率越高。即miR-29c-3p在食管腺癌中可能是肿瘤抑制基因，AC131009.3、AC087388.1、CDK1、TPX2 在食管腺癌中可能是促肿瘤生长基因。

图2 ceRNA网络中各基因表达水平与生存时间的关系

2.3.3 ceRNA网络中各基因在食管腺癌组织、正常食管组织中的表达差异 见图3。

由图3可见，与正常食管组织比较，在食管腺癌组织中 ACO87388.1、AC131009.3、CDK1、TPX2 的表达水平更高，而miR-29c-3p在食管腺癌组织中的表达水平更低，差异有统计学意义。ceRNA网络中各基因在食管腺癌组织、正常食管组织中的表达差异趋势与食管腺癌患者生存分析结果相符。

图3 ceRNA网络中各基因在食管腺癌组织、正常食管组织中的表达差异

2.3.4 miR-29c-3p 表达量与 ACO87388.1、AC131009.3、CDK1、TPX2表达量的相关性分析 见图4。

由图4可见，在共表达分析中，miR-29c-3p的表达量与ceRNA网络中的lncRNA（AC087388.1、AC131009.3）、mRNA（CDK1、TPX2）的表达量存在明显的Pearson相关，相关系数分别为-0.473、-0.564、-0.655、-0.647。

图4 miR-29c-3p表达量与ACO87388.1、AC131009.3、CDK1、TPX2表达量的相关性分析

2.4 GSEA中的可能与肿瘤相关的KEGG富集结果CDK1、TPX2基因参与的可能与肿瘤预后相关的KEGG见表1。

由表1可见，KEGG显示的是CDK1高表达组、CDK1低表达组、TPX2高表达组、TPX2低表达组中NSE排名前5的KEGG通路；CDK1基因可能参与肿瘤的免疫、细胞周期等生物学过程有关，TPX2基因主要与肿瘤的代谢等。对表1进行可视化，CDK1和TPX2基因参与的KEGG结果见图5（插页）。

表1 CDK1、TPX2基因参与的可能与肿瘤预后相关的KEGG

图5 细胞周期蛋白激酶1（CDK1）和微管成核因子（TPX2）基因参与的KEGG结果（a：CDK1；b：TPX2）

由图5可见，GSEA结果提示 CDK1、TPX2基因在食管腺癌中可参与多种高级生物学功能，以及各生物学功能在CDK1高表达组、CDK1低表达组、TPX2高表达组、TPX2低表达组中的分布情况。

3 讨论

食管腺癌是常见的食管恶性肿瘤，肿瘤细胞对常规抗癌治疗的耐药率很高，晚期食管癌的5年生存率仅为0.9%[4]。食管腺癌表观遗传现象异常引起的生物学功能改变是导致肿瘤患者预后不佳的重要因素之一，而基因表达水平的改变是表现遗传现象建立的核心。基因表达水平的异常是导致肿瘤预后差异的起源，分子间的信号传导是预后改变的重要机制。Zhang等[5]就曾利用TCGA中的食管腺癌数据并采用权重基因共表达网络分析方法探索食管腺癌中的预后基因。近年来已有大量文献表明，具有差异表达的lncRNA、miRNA可能在食管癌患者的预后中扮演者重要的作用[6-10]。非编码RNA的生物学作用主要是通过复杂的互作网络最终作用于具有蛋白编码功能的mRNA而发挥作用[11]。因此，研究食管腺癌中lncRNA、miRNA、mRNA之间的竞争网络来探索肿瘤的生物学价值可能为肿瘤的治疗提供新的治疗手段[12]。

生物信息分析技术是一种高效率、精准确率、高可信度的数据挖掘技术，已被《Nature》等知名期刊认可[13-14]。本研究中笔者通过生物信息分析筛选食管腺癌中存在差异表达且与预后相关的lncRNA、miRNA、mRNA，再挖掘到以miRNA为中心的ceRNA网络。ceRNA网络构建方式是多样化的，例如可以通过靶基因预测工具或权重基因共表达网络分析等方法来实现[15-17]，本研究中的ceRNA是通过分析相应样本中lncRNA、miRNA、mRNA的共表达关系的权重系数获得的，研究结果可能更贴近实际，可靠性更强。该ceRNA网络的核心是miR-29c-3p，多项研究表明miR-29c-3p在ceRNA中可以抑制卵巢癌生长[18]，在室管膜瘤中可以改变肿瘤细胞对化疗药物的敏感性等[19]。在食管腺癌中，Craig等[20]学者报道miR-29c-3p出现了异常表达，与本研究结果一致。本研究共表达分析结果提示miR-29c-3p与ACO87388.1、AC131009.3、CDK1、TPX2 呈负相关关系，与CDK1、TPX2的相关性最为密切（Cor＞-0.6）。然而，食管腺癌中关于miR-29c-3p在ceRNA网络中的作用及其机制报道甚少。在现有研究中，ceRNA网络中的ACO87388.1、AC131009.3与恶性肿瘤的相关性尚不明确，但是已有证据表明CDK1、TPX2是重要的肿瘤预后相关分子[21-23]。因此，进一步地研究ceRNA网络在食管腺癌中的生物学作用显得极为迫切。

lncRNA、miRNA最终是通过影响mRNA的翻译或使mRNA降解而发挥生物学作用，因此CDK1、TPX2是本研究的核心分子。基于CDK1、TPX2基因RNA表达量进行的GSEA结果认为，本ceRNA网络可能参与了食管腺癌的多种生物学功能，例如影响食管腺癌的离子通到、细胞代谢、肿瘤微环境等途径从而调控肿瘤细胞的增殖、侵袭等。在恶性肿瘤中，离子通道的转运因改变了肿瘤的代谢而影响了患者的预后，是恶性肿瘤靶向治疗的靶点[24-25]。此外，肿瘤微环境可通过表达异常的分子建立强大的肿瘤免疫逃逸系统，在食管鳞癌中肿瘤的微环境可影响肿瘤的生长[26]，在食管腺癌也可能存在类似的机制。

综上所述，食管腺癌是一种高度恶性的肿瘤，本研究发现 ACO87388.1、AC131009.3、miR-29c-3p、CDK1、TPX2是影响食管腺癌患者的编码、非编码基因，以miR-29c-3p为中心的ceRNA可能参与了食管腺癌肿瘤细胞的发生、发展过程，在未来食管腺癌的预后分析中存在重要的研究价值。本研究仅局限于生物信息的分析研究，研究结论尚需进一步的基础实验、临床研究去证实，但研究结果为未来的研究探索提供了思路。