胰腺癌关键基因的筛选和分析：多数据联合生物信息学分析

洪坤巧，邱 虎，殷海森，杨益波，王 琰，余保平

武汉大学人民医院消化内科，湖北 武汉 430060

胰腺癌是一种常见的恶性肿瘤，每年造成超过331 000人死亡，在癌症死亡原因里排名第七，胰腺癌发病率和死亡率的趋势在世界范围内差异很大，总体5年生存率约6%(2%～9%)[1]。现在认为，胰腺癌通常伴随着KRAS癌基因的激活突变，这是肿瘤发生和发展的有效驱动力，抑癌基因(如CDKN2A/p16、TP53和SMAD4)的失活突变与KRAS突变协同作用，导致胰腺癌的侵袭性生长[2]，但胰腺癌的具体发生机制仍不十分明确。胰腺癌中90%的病理类型是腺癌，在早期无有效的无创诊断方法，影像学检查难以发现早期微小病变，并且疾病进展速度较快，迅速扩散到周围脏器，是最致命的癌症之一[3]。早期患者往往无明显临床症状，多数患者确诊时已发生局部器官浸润和全身转移，失去手术根治性切除的机会。因此，探讨胰腺癌发病的分子机制，特别是寻找在胰腺癌组织中高度特异性表达的基因，对于胰腺癌的早期诊断和治疗具有重要意义。

在过去的几十年里，微阵列技术和生物信息学分析被广泛应用于筛选基因转录组水平上的改变，这帮助我们识别与胰腺癌的癌变和进展相关的差异表达基因和功能通路。然而，独立芯片分析的假阳性率使其难以获得可靠的结果。因此，本研究从Gene Expression Omnibus(GEO)下载3个mRNA芯片数据集进行分析，获得胰腺癌组织与癌旁组织的差异表达基因。随后，我们进行了GO、KEGG途径富集分析和蛋白质相互作用(protein-protein interactio，PPI)网络分析，以帮助我们了解潜在的致癌和促癌的分子机制。

1 材料与方法

1.1 芯片数据下载和整理GEO(https://www.ncbi.nlm.nih.gov/geo/)[4]是公共的功能基因组数据库，从GEO下载GSE15471、GSE16515、GSE28735三个数据集和相对应的平台文件GPL7691、GPL570、GPL1261。GSE15471包含39个正常胰腺样本和39个胰腺癌样本，GSE16515包含15个正常胰腺样本，37个胰腺癌样本，GSE28735包含45个正常胰腺样本和45个胰腺癌样本。将数据进行归一化，标准化，ID转换处理，最终实现差异表达基因在多个联合数据的综合分析。

1.2 筛选差异表达基因筛选胰腺肿瘤组织和非肿瘤组织的差异表达基因，利用GEO在线分析工具GEO2R(https://www.ncbi.nlm.nih.gov/geo/geo2r)进行分析，这是一个基于R的Web应用程序，可以帮助用户分析GEO数据，根据实验条件对样本进行分组，用于比较两组或两组以上不同样本之间差异有统计学意义的基因。筛选出的TOP250基因包括P值，|log2FC|以及相关的基因信息，筛选出的差异表达基因导入EXCEL表格，根据|log2FC|>1，P-value<0.01条件对数据进行过滤。利用WEEN在线工具(http://bioinformatics.psb.ugent.be/beg/)[5]，找出三个数据集的共同差异表达基因。

1.3 PPI网络构建与模块分析采用STRING(https://string-db.org/)[6]数据库预测PPI网络，STRING数据库是一个搜寻已知蛋白质和预测蛋白质之间相互作用的系统，建立一个全面、客观的蛋白网络，包括直接(物理)和间接(功能)交互。分析蛋白之间的相互作用和关系，为疾病的发生发展机制研究，甚至诊断及预后提供思路。差异表达基因之间的网络图通过STRING数据库进行构建，导出网络图数据在Cytoscape进一步优化和分析。Cytoscape是一款图形化显示网络并进行分析和编辑的软件，安装Cytoscape的插件MCODE聚类网络中的基因富集模块。使用Cytoscape绘制PPI网络，MCODE识别PPI网络中最重要的模块，选择标准是：MCODE scores >5、degree cut-off=2、node score cut-off=0.2、Max depth=100、k-score=2。

1.4 差异表达基因的KEGG、GO通路富集分析功能富集和通路分析在高通量组学数据的生物学解释中发挥了关键作用，利用在线分析工具WebGestalt(http://www.webgestalt.org/)[7]对差异表达基因进行KEGG富集分析和GO富集分析，WebGestalt是一个专注于富集分析的在线网站，支持多种富集分析算法，而且涵盖的功能注释数据库较为全面。KEGG是一个从分子水平信息，特别是基因组测序和其他高通量实验技术产生的大规模分子数据集，了解细胞、有机体、生态系统等生物系统的高级功能和效用的数据库资源[8]。GO知识库是世界上最大的基因功能信息来源，这些知识是人类可读和机器可读的，是生物医学研究中大规模分子生物学和遗传学实验计算分析的基础[9]。利用WebGestalt在线数据库进行GO和KEGG生物分析，实现基因与生物学过程可视化的连接，P<0.05为差异有统计学意义。

1.5 关键基因筛选和分析在network中MCODE的每个模块的seed基因为关键基因。利用UCSC Xena构建关键基因的层次聚类分析，UCSC Xena(http://genome-cancer.ucsc.edu)[10]是一个在线公共工具数据库，支持几乎所有功能基因组学数据，包括SNV、INDEL、大型结构变异、CNV、表达、DNA甲基化、ATAC-seq信号和表型注释，该数据库包含多基因组和临床/表型数据，通过探索功能基因组数据集，以寻找基因组和表型变量之间的相关性。关键基因的生物学过程利用iDEP进行分析，cBioPotal分析8个关键基因的突变和缺失，探索可视化和分析多维度癌症基因组数据，cBioPortal(https://www.cbioportal.org/)[11]是一个开放平台，用于在临床数据和生物学途径的背景下交互式探索多维癌症基因组学数据集，基因组数据类型包括体细胞突变、DNA拷贝数改变、mRNA和miRNA表达、DNA甲基化、蛋白质富集、磷酸化蛋白富集。GEPIA(http://gepia.cancer-pku.cn/)[12]是一个癌症数据挖掘网站，主要基于TGCA和GTExprojects，GEPIA分析8个关键基因在胰腺癌和正常组织中的差异表达和生存预后的关系，DDX60、ERP27在多种不同肿瘤组织中的表达情况。

2 结果

2.1 筛选差异表达基因从GEO2R在线下载的3个数据集按照|log2FC|>1，P-value<0.01条件进行筛选后GSE15471、GSE16515、GSE28735分别包含1 761、1 694、407个基因，其中取交集后相同基因247个，其中上调基因190个，下调基因57个(见图1A)。

2.2 PPI网络构建与模块分析利用STRING数据库预测PPI网络，在Cytoscape进一步加工，190个上调基因被标记为红色，57个下调基因被标记为绿色，节点的直径较大者为8个hub基因(见图1B)，PPI网络图中最突出的富集团块含有13个节点，46条边(见图1C)。

图1 韦恩图(A)、差异表达基因最显著模块(B)、PPI网络图(C)Fig 1 Venn diagram (A), the most module of differentially expressed genes (B), PPI network (C)

2.3 差异表达基因的KEGG、GO富集分析通过WebGestalt对差异表达基因进行生物学分类，功能和通路的富集分析。GO富集分析结果显示，差异表达基因生物学过程(biological processes，BP)明显富集在生物调节、代谢过程、刺激反应、细胞信号传导、细胞增殖、生长。细胞成分(cellular component，CC)主要富集在细胞膜、内膜系统、核、胞质、内质网。差异表达基因分子功能(molecular function，MF)的改变主要富集在蛋白结合、离子结合、水解酶活性、结构分子活性、转移酶活性、核苷酸结合。KEGG通路分析显示，上调基因主要富集在ECM-receptor相互作用、阿米巴病、局部粘附、蛋白消化吸收、PI3K-Akt信号通路，而下调基因主要富集在胰液分泌、肾素-血管紧张素系统、补体和凝血级联(见图2)。

注：每个生物过程、细胞成分和分子功能类别分别用红、蓝、绿条表示。每一个柱子代表一个种类，柱子的高度表示基因ID的数量，FDR<0.05。

2.4 关键基因表达与生存预后分析在这些基因中，COMP、ANLN、VCAN、DDX60、MET呈高表达；ERP27、NR5A2、C5呈低表达(见图3)。Oncomin分析每个关键基因的突变：COMP：4%，ANLN：10%，VCAN：6%、DDX60：7%、ERP27：5%、MET：8%、NR5A2：4%、C5：6%(见图4A)。通过GEPIA数据库我们发现，ANLN、MET基因组学的改变使总体生存率和无病生存率显著降低(ANLN：P=0.015，P=0.00046；MET：P=0.00034，P=0.00067)，DDX60、ERP27的改变降低总体生存率，但对无病生存率无明显影响(DDX60：P=0.0019，P=0.26；ERP27：P=0.042，P=0.38)(见图5)。DDX60、ERP27在不同组织中的表达分析显示：DDX60在弥漫性大B细胞淋巴瘤、食管癌、多形成性胶质细胞瘤、头颈鳞状细胞癌、急性髓细胞样白血病、脑低级别胶质瘤、胰腺癌、胸腺癌中呈高表达；ERP27在乳腺浸润癌、结肠癌、肾嫌色细胞癌、急性髓细胞样白血病、卵巢浆液性囊腺癌、甲状腺癌、子宫内膜癌等多种肿瘤中呈高表达，但在胰腺癌中明显呈低表达(见图4B～4C)。利用CCLE(https://portals.broadinstitute.org/ccle/)分析ERP27和DDX60在不同组织细胞中的表达以及在胰腺癌细胞不同株系中的表达(见图4D～4E)，DDX60在绝大多数胰腺癌细胞株系中呈高表达，而ERP27则呈低表达。

注：*P<0.01。

图4 关键基因在不同组织中的表达和突变分析 A：cBioportal分析8个关键基因的突变分析，OncoPrint概述了影响特定单个样品(列)的特定基因(行)的基因组改变；B～C：GEPIA分析DDX60和ERP27在不同组织中的表达，横轴代表不同的组织，肿瘤组织标记为红色，非肿瘤组织标记为黑色，纵轴代表基因的表达量。D～E：来源于Dccle对DDX60和ERP27在胰腺癌不同细胞系中的表达

图5 GEPIA数据库分析8个关键基因对生存预后的影响：利用GEPIA在线平台对关键基因进行总体生存分析和无病生存分析Fig 5 GEPIA database analysis of the effects of eight key genes on survival and prognosis: the overall survival and disease free survival of key genes were analyzed using the GEPIA online platform

3 讨论

胰腺癌在世界癌症死亡原因中排第七，并且死亡率与不良预后息息相关[13]，胰腺癌每年的患病率呈持续性增长，由于胰腺癌发病与年龄相关，因此人口老龄化将导致胰腺癌的发病率增加[14]。胰腺癌恶性程度较高，进展迅速，但起病隐匿，早期症状不典型，大多数患者被诊断为中-晚期，失去手术机会，因此癌前病变或肿瘤早期的识别至关重要，然而合适的筛选方法目前还未找到[15]。尽管随着腹腔镜手术技术、新辅助放化疗、靶向及免疫治疗的不断进展，但仍未取得较好的治疗结果。识别新的生物标志物是迈向精准医疗时代的迫切需求，避免对患者生存质量无效的医疗浪费[16]。鉴于胰腺癌的不良预后和对其的认知不足，进一步的研究对于更好地了解这种疾病，早期诊断和提高生存率至关重要。近年分子生物学和高通量测序技术的发展，在全基因组、转录组、蛋白质组水平筛选胰腺癌分子标志物成为了可能，有望提高该病的早期诊断率。

胰腺癌的发病机制目前仍知之甚少，约97%的胰腺癌具有基因异常，如点突变、扩增、缺失、易位和倒位[17]，2/3的肿瘤具有与有丝分裂错误相关的复杂重排模式[18]。本研究通过Oncomin分析发现,每个关键基因均有不同程度的突变，成为胰腺癌发生相关的潜在基因(COMP：4%、ANLN：10%、VCAN：6%、DDX60：7%、ERP27：5%、MET：8%、NR5A2：4%、C5：6%)。胰腺癌常见的驱动基因：KRAS、CDKN2A、TP53和SMAD4，致癌基因KRAS，通过激活MAPK依赖的信号通路MYC上调和无氧氧化戊糖磷酸途径的转录，从而促进核酸生物合成，重新连接细胞代谢，促进肿瘤细胞的增殖[19]，TP53功能改变导致G1细胞周期阻滞，促进肿瘤淋巴管数量的增加和转移，TP53在细胞周期和凋亡中起着重要的作用[20-21]。SMAD4作为转录生长因子(TGF-β)信号通路的关键因子调控细胞增殖、粘附和运动[22]。CDKN2A在胰腺癌中高表达，甲基化与胰腺癌风险增加相关[23]。然而，胰腺癌具体发生机制仍不清楚，胰腺癌的高度恶性和不良预后成为一个亟需解决的难题。

本研究通过分析3个mRNA微阵列数据集，获得胰腺癌组织与非癌组织间的差异基因，3个数据集中共检测到247个差异表达基因，其中上调的有190个，下调的有57个。进行GO、KEGG富集分析，以探索差异表达基因之间的相互作用，上调基因主要富集在细胞外基质相互作用、局部粘附、蛋白消化吸收、PI3K-Akt信号通路，而下调基因主要富集在胰液分泌、肾素-血管紧张素系统、补体和凝血级联。新的研究发现，细胞外基质在肿瘤微环境的形成中起着核心的、动态的作用，癌细胞对细胞外基质重构的失调促进了不可逆的蛋白水解和交联，进而影响细胞信号、微环境、血管生成和组织生长[24]。PI3K-Akt信号通路参与多种癌症的发生，如乳腺癌[25]、膀胱癌[26]、肺癌[27]等。

筛选出8个关键基因(深度≥10)，分别是COMP、ANLN、VCAN、DDX60、ERP27、MET、NR5A2、C5，通过GEPIA数据库我们发现ANLN、MET、DDX60、ERP27对总体生存率有影响，提示这些基因可能在胰腺癌的发生、发展、侵袭或复发等过程中发挥重要作用，其中DDX60、ERP27在胰腺癌中的作用鲜有报道。DDX60是一种DEXD/H box RNA解旋酶，作为抗病毒因子具有与病毒RNA和DNA结合的活性，促进RIG-I样受体介导的信号传导，对Ⅰ型干扰素(IFN)和其他炎性细胞因子的产生具有重要作用[28]，也可能通过TBK1信号通路激活IFN[29]，在大鼠的黑色素肿瘤中，IFNβ信使RNA和蛋白质水平显著提高，DNA传感器DDX60有上调的趋势，这些影响伴随着抑制肿瘤的生长和增加肿瘤的坏死[30]。在针对胆囊收缩素-2受体(cholecystokinin-2 receptor，CCK2R)介导胰腺癌化学预防的实验中，CCK2R的拮抗剂YF476能有效地抑制胰腺癌，并且与未处理组相比，DDX60呈明显的低表达，这提示DDX60在胰腺癌中与非癌组织中呈现明显的差异表达[31]。最新研究发现，DDX60与结直肠癌的发生和发展相关[32]，然而DDX60在胰腺癌中的作用机制到目前为止不确切。本次研究通过对DDX60在组织的中的差异表达分析，可以看出，与其他肿瘤组织相比DDX60在胰腺癌中呈高表达，排名第二，并且在胰腺癌的不同细胞系中，几乎均呈高表达，DDX60可作为胰腺癌诊断和预后的指标。ERP27是定位于内质网的蛋白质二硫异构酶的同源物，且它是PDI家族中一个不具有氧化还原活性的成员，ERP27能够区分折叠和未折叠的底物，只是与后者相互作用[33]。人ERP27在许多不同的人体组织中均有表达，如骨髓、肺、肾脏和脾脏，但最显著的表达部位是胰腺[34]。ERP27域的功能尚不清楚，C端域已被证明与Δ-生长抑素肽和ERP57相互作用，ERP27在内质网应激过程中表达上调，从而可能使其与积累的错误折叠的底物结合，并将其呈现给ERP57进行催化[35]。ERP57促进多种癌症的发生，如促进乳腺癌细胞的生长[36]，增强喉癌细胞的耐辐射，并作为喉癌的预后标志物[37]。内质网参与蛋白的合成与加工，与稳态相关，蛋白质在内质网中的错误折叠会导致内质网应激和未折叠蛋白反应的激活，在大多数主要类型的人类癌症，包括乳腺癌、胰腺癌、肺癌、皮肤癌、前列腺癌、脑癌，甚至是液体恶性肿瘤中，均发现了强内子网应激反应，此外，原位内质网应激与晚期疾病和耐药密切相关[38]。本研究发现，ERP27在乳腺浸润癌、结肠癌、子宫内膜癌等多种肿瘤中呈高表达，但在胰腺癌中明显呈低表达，可能作为机体的一个保护因子，抑制肿瘤的发生发展，ERP27的预后生存分析发现，总体生存率差异有统计学意义，可以作为胰腺癌预后的生物标志物。

软骨寡聚基质蛋白(cartilage oligomeric matrix protein，COMP)是一种524 KDa的可溶性五聚糖蛋白，是一种很有价值的软骨细胞再生物质[39]，在早期结肠癌组织mRNA和蛋白水平明显升高，在肿瘤发生中促进细胞增殖[40]，肝星状细胞来源的COMP有助于肝细胞癌的进展[41]。ANLN是一种肌动蛋白结合蛋白，在许多恶性肿瘤中起重要作用，如膀胱癌[42]、乳腺癌[43]，ANLN通过调控EZH2/miR-218-5p/LASP1信号轴促进胰腺癌的进展[44]。VCAN基因在多个乳腺癌细胞中编码了硫酸软骨素蛋白聚糖，促进肿瘤的迁袭和转移[45]，VCAN表达可作为Ⅱ～Ⅲ期结肠癌复发的潜在预后生物标志物[46]。跨膜酪氨酸激酶受体(MET)由肝细胞生长因子激活，大量研究表明，MET在胰腺癌中高表达，参与肿瘤的发生发展，并且可作为胰腺癌的诊断预后指标和治疗靶点[47-50]。孤核受体(orphan nuclear receptor，NR5A2)属于核受体的NR5A亚家族，在内胚层起源的发育和成体组织中均有表达[51]，一项荟萃分析显示，NR5A2可能对胰腺癌有保护作用[51]，在胰腺癌中可能通过抑制胰腺癌干细胞和上皮间质转换发挥作用[51]。C5基因编码了补体系统的一个组成部分，这是先天免疫系统的一部分，在炎症、宿主稳态和宿主防御病原体方面起着重要作用[52]，补体的激活系统已与肿瘤的发生相关[53]。

总之，ANLN、MET、DDX60、ERP27等基因的改变会影响胰腺癌患者的整体生存质量，可作为生存预后判断的生物标志物，其中DDX60在胰腺癌组织和胰腺癌各细胞系中几乎均呈高表达，可作为胰腺癌诊断的特异性标志物，ERP27在胰腺组织中呈特异性的高表达，在胰腺癌中显著降低，提示可能成为胰腺的保护因子、胰腺癌治疗的新靶标和预后生物分子。