基于生物信息学的结肠癌核心基因筛选及验证

曾成， 祁国萍， 沈颖，2， 陆文斌，2， 邓建忠，2， 刘迁，2， 金建华，2

(1. 江苏大学附属武进医院肿瘤科，江苏 常州 213017；2. 徐州医科大学武进临床学院肿瘤科，江苏 常州 213017)

结直肠癌是常见的消化道恶性肿瘤，其全球发病率仅次于乳腺癌、肺癌，居第3位，死亡率居癌症相关死亡的第2位[1]。目前认为结肠癌与遗传、环境、饮食、炎症性肠病、性别、种族等相关[2]，但其发病确切分子机制尚不清楚，临床亦缺乏早期诊断和预后判断的分子标志物。基于高通量测序及基因芯片技术的快速发展，有研究者通过生物信息学方法筛选结肠癌的潜在核心基因[3-4]，但由于数据集和筛选方法不同，筛选结果存在差异。本研究通过分析GEO数据库中结肠癌数据集，筛选结肠癌核心基因，利用TCGA数据库验证核心基因表达水平并做生存分析，最后通过体外实验探讨预后相关基因对结肠癌细胞增殖的影响，旨在获得更多与结肠癌发生发展相关的潜在核心基因。

1 材料和方法

1.1 数据集来源

从GEO数据库(www.ncbi.nlm.nih.gov/geo/)中下载数据集，筛选条件：人类结肠癌全基因组表达谱数据，样本中包含结肠癌组织及癌旁组织，样本量>50。本研究选择GSE23878、GSE37182和GSE74602数据集进行后续分析。GSE23878数据集基于GPL570平台，包含结肠癌组织35例，癌旁组织24例；GSE37182数据集基于GPL6947平台，包含结肠癌组织84例，癌旁组织88例；GSE74602数据集基于GPL6104平台，包含结肠癌组织30例，癌旁组织30例。

1.2 细胞系及试剂

人结肠癌HCT116、SW620细胞购自中国科学院上海生物化学与细胞研究所。DMEM、RPMI 1640培养基、胎牛血清均为美国Gibco公司产品；pECMV、pECMV-AURKA、pEGFP、pEGFP-TIMP1质粒由淼灵质粒平台构建；转染试剂TurboFect Transfection、兔抗人Aurora A蛋白激酶(AURKA)一抗、兔抗人β-肌动蛋白一抗、山羊抗兔二抗、MTT试剂盒、5×上样缓冲液、细胞裂解液均购自上海生工生物工程有限公司；兔抗人基质金属蛋白酶抑制物1(TIMP1)一抗购自武汉三鹰生物技术有限公司；ECL化学发光液购自合肥兰杰柯科技有限公司；BCA蛋白定量试剂盒购自上海碧云天生物技术有限公司。

1.3 筛选差异表达基因

在R软件(版本：4.0.2)中利用ggplot2、limma、pheatmap等软件包处理3个数据集。差异表达基因(differentially expressed genes，DEGs)筛选标准如下：错误发现率(false discovery rate，FDR)校正后P值<0.05，且|log2FC|>1。倍数变化(fold change，FC)表示DEGs的差值倍数。对3个数据集筛选出的DEGs取交集，利用VennDiagram软件包作韦恩图，得到共有DEGs。

1.4 共有DEGs的基因本体论富集分析和京都基因与基因组百科全书信号通路分析

在R软件(版本：4.0.2)中利用clusterProfiler、org.Hs.eg.db、enrichplot、ggplot2等软件包对共有DEGs进行基因本体论(gene ontology，GO)富集分析和京都基因与基因组百科全书(kyoto encyclopedia of genes and genomes，KEGG)信号通路分析，GO富集分析包括分子功能、生物学过程及细胞组成。

1.5 构建蛋白质相互作用网络及筛选核心基因

利用STRING(https：//string-db.org/)在线数据库对共有DEGs构建蛋白质相互作用(protein-protein interaction， PPI)网络，设置有效结合分数>0.4。随后将PPI网络导入Cytoscape软件(版本：3.7.2)，利用CytoHubba插件筛选核心基因。CytoHubba插件中有EPC、BottleNeck、EcCentricity、Closeness、Radiality、Betweenness、Stress、ClusteringCoefficient、MCC、DMNC、MNC和Degree 共12种分析方法[5]。每种分析方法选择前10个基因，所有基因按照在12种分析方法中出现次数进行排序，最后选择出现次数最多的10个基因作为核心基因，并构建核心基因的PPI网络。

1.6 基于TCGA数据库的核心基因验证和生存分析

利用基于TCGA数据库的GEPIA在线网站(http：//gepia.cancer-pku.cn/ index.html)验证核心基因在结肠癌及癌旁组织表达水平，并对核心基因进行生存分析。

1.7 细胞培养与转染

在37 ℃、5% CO2条件下，将HCT116、SW620细胞分别置于含10%胎牛血清DMEM、RPMI 1640培养基中培养。取对数生长期细胞接种于6孔板，待细胞融合度达80%时，按照转染说明书，将pECMV(空载体)、pECMV-AURKA、pEGFP(空载体)、pEGFP-TIMP1质粒各2 μg分别转染HCT116和SW620细胞。

1.8 蛋白质印迹法检测细胞AURKA、TIMP1表达

收集“1.7”转染48 h后各组细胞，加入含1% PMSF的RIPA裂解液，冰上裂解30 min；4 ℃、15 000 r/min离心20 min后取上清液，BCA法行蛋白定量；加入5×上样缓冲液，100 ℃沸水浴5 min；分别取40 μg样品行SDS-PAGE，80 V分离蛋白；200 mA 90 min冰上转PVDF膜；5%脱脂牛奶室温封闭2 h；PBST洗膜3次，5 min/次；加入相应的用5% 脱脂牛奶稀释的一抗β-肌动蛋白(1 ∶4 000，内参)、AURKA和TIMP1(均为1 ∶1 000)， 4 ℃孵育过夜；PBST洗膜3次，5 min/次；加入用PBS稀释的二抗(1 ∶5 000)室温孵育1 h；PBST洗膜3次，5 min/次；用ECL化学发光液暗室曝光；用Image J 软件处理分析蛋白条带。

1.9 MTT法检测细胞增殖

取“1.7”转染后各组细胞悬液接种于96孔板，每孔2×103个，每组5个复孔。5% CO2、37 ℃分别培养24、48、72 h，每孔加入20 μL MTT(5 mg/mL)，37 ℃孵育4 h；弃上清液，每孔加入150 μL二甲基亚砜，置摇床上低速振荡10 min，待结晶物充分溶解，于酶联检测仪490 nm波长处测定光密度(D)值，计算细胞增殖率。细胞增殖率(%)=(实验组D值-空白组D值)/(对照组D值-空白组D值)×100%。

1.10 统计学方法

2 结果

2.1 DEGs筛选结果

通过R软件分析，GSE23878数据集DEGs 1 411个，其中上调基因502个，下调基因909个；GSE37182数据集DEGs 627个，其中上调基因261个，下调基因366个；GSE74602数据集DEGs 1 485个，其中上调基因816个，下调基因669个。3个数据集DEGs通过VennDiagram软件包取交集后得到270个共有DEGs(图1)。

图1 DEGs韦恩图

2.2 共有DEGs的GO富集分析及KEGG信号通路分析

GO富集分析结果显示，共有DEGs生物学过程主要与生长负调控、细胞外基质组织、细胞外结构组织等相关；细胞组成主要富集于含胶原的细胞外基质、收缩纤维、肌节等；分子功能主要富集于受体配体活性、信号转导受体激活剂、糖胺聚糖结合等。KEGG信号通路富集结果表明，共有DEGs主要与矿物质吸收、Wnt、紧密连接、NF-κB、细胞周期、细胞黏附分子等信号通路相关。见图2。

图2 共有DEGs的GO富集分析及KEGG信号通路分析

2.3 PPI网络构建及核心基因的筛选

利用STRING数据库对共有DEGs进行PPI网络分析，结果显示，PPI网络具有明显的交互作用(P<1.0×10-16)。将结果导入Cytoscape软件(图3A)，利用插件CytoHubba中的12种算法筛选出重复次数最多的前10个核心基因，由于第11个基因重复出现的次数和第10个基因重复出现的次数相同，本研究也将其纳入核心基因，故最终获得11个核心基因。核心基因在12种算法中重复出现次数如下：MYC8次、基质金属蛋白酶抑制物1(TIMP1)6次、泛素偶联酶E2C(UBE2C)5次、小窝蛋白1(CAV1)5次、Y 染色体中性别决定区相关的高迁移率组框9(SOX9)5次、C-X-C型趋化因子配体12(CXCL12) 5次、Aurora A蛋白激酶(AURKA) 4次、Ⅰ型胶原蛋白α1链(COL1A1)4次、细胞周期分裂蛋白20(CDC20)4次、DNA拓扑异构酶Ⅱα(TOP2A)4次、着丝粒蛋白F(CENPF)4次。11个核心基因重新导入STRING数据库后构建核心基因的PPI网络(图3B)，结果显示有明显交互作用(P=7.66×10-7)。

A:共有DEGs的PPI网络；B:核心基因的PPI网络

2.4 核心基因表达验证及生存分析

GEPIA在线网站包含TCGA数据库中275例结肠癌组织和41例癌旁组织相关信息。11个核心基因通过GEPIA网站验证显示，其中MYC、TIMP1、UBE2C、SOX9、AURKA、COL1A1、CDC20、TOP2A和CENPF共9个基因在结肠癌中呈高表达，而CAV1和CXCL12在结肠癌中呈低表达，差异有统计学意义(P<0.05)，见图4。

图4 核心基因在结肠癌和癌旁组织表达水平的验证

此外，本研究也利用该网站对核心基因进行生存分析，发现TIMP1表达与结肠癌患者总体生存期(overall survival，OS)呈显著负相关(P<0.05)，即TIMP1在结肠癌样本中表达越高，结肠癌患者OS越短，AURKA表达与结肠癌患者OS呈显著正相关(P<0.05)，而COL1A1表达与结肠癌患者无病生存期呈显著负相关(P<0.05)，见图5。

图5 3个核心基因的生存分析

2.5 TIMP1或AURKA对结肠癌细胞增殖能力的影响

为进一步验证影响结肠癌患者OS相关基因(TIMP1和AURKA)在结肠癌细胞中的功能，本研究分别将pECMV，pECMV-AURKA，pEGFP，pEGFP-TIMP1质粒转染入HCT116和SW620细胞。蛋白质印迹结果显示，HCT116、SW620细胞pECMV-AURKA组AURKA蛋白相对表达量显著高于相应pECMV组(P均<0.01)；HCT116、SW620细胞 pEGFP-TIMP1组TIMP1蛋白表达相对量显著高于相应pEGFP 组(P均<0.01)，见图6。MTT结果显示，HCT116、SW620细胞pECMV-AURKA 组72 h细胞增殖能力明显高于pECMV组(t=4.039，5.731，P均<0.05)，24、48 h两组间差异无统计学意义。HCT116细胞 pEGFP-TIMP1组48、72 h细胞增殖能力明显高于pEGFP组(t=11.716，5.673，P均<0.01)，24 h两组间差异无统计学意义。SW620细胞pEGFP-TIMP1组72 h细胞增殖能力明显高于pEGFP组(t=5.920，P<0.01)，24、48 h两组间差异无统计学意义(图7)。

图6 蛋白质印迹检测HCT116和SW620细胞中AURKA和TIMP1蛋白表达

a：P<0.05，与同时间点pECMV组比较；b：P<0.01， 与同时间点pEGFP组比较

3 讨论

本研究通过下载GEO数据库中GSE23878、GSE37182和GSE74602数据集，经分析获得270个共有DEGs。GO富集和KEGG信号通路分析显示这些基因主要与生长负调控、受体配体活性、信号转导受体激活剂、Wnt信号通路、NF-κB信号通路、细胞周期信号通路等有关。在筛选核心基因时，为了减少单一算法的局限性，本研究充分利用CytoHubba软件中12种算法，将每一种算法的前10个核心基因按照重复出现的次数进行排序，选取重复次数最多的11个核心基因(MYC、TIMP1、UBE2C、CAV1、SOX9、CXCL12、AURKA、COL1A1、CDC20、TOP2A和CENPF)进行后续验证分析。为了验证GEO数据集筛选结果的准确性，将11个核心基因通过基于TCGA数据库的GEPIA在线网站进行验证，结果显示11个核心基因在TCGA数据库的表达水平与GEO数据库筛选结果一致，其中MYC、TIMP1、UBE2C、SOX9、AURKA、COL1A1、CDC20、TOP2A和CENPF在结肠癌中高表达，而CAV1和CXCL12在结肠癌中低表达。

MYC蛋白是一种转录因子，在众多肿瘤中发挥重要作用。UBE2C蛋白是泛素-蛋白酶体系统的重要组成部分，UBE2C高表达与多种肿瘤不良预后相关[6]，但其在结肠癌中的分子机制研究较少。CAV1蛋白是细胞膜主要支架蛋白，CAV1在结肠癌细胞中过表达能够引起启动子CpG位点低甲基化，进而促进细胞增殖[7]。SOX9蛋白是转录因子SOX家族成员之一，陈玉昌等[8]研究显示，SOX9 mRNA和蛋白表达在结肠癌组织中均上调，且SOX9蛋白与肿瘤的分化程度、TNM分期和淋巴转移有关。CXCL12蛋白是与CXCR4(G蛋白偶联受体)特异性结合的一种趋化因子，CXCL12/CXCR4信号与多种肿瘤的侵袭、迁移密切相关[9-10]。TOP2A蛋白是一种酶蛋白，与细胞增殖、凋亡及有丝分裂相关。研究显示，TOP2A敲减可通过影响凋亡蛋白(Bcl-2、Bax)、侵袭相关蛋白(MMP2、MMP9)表达抑制结肠癌细胞的增殖和侵袭[11]。CENPF蛋白是着丝粒蛋白家族成员之一，在有丝分裂及肿瘤中起重要作用，被认为是结肠腺癌的新型分子标志物[12]。CDC20蛋白是细胞周期相关蛋白，其异常表达使有丝分裂发生错误，从而导致一些癌基因的过表达及抑癌基因的失活，最终抑制肿瘤细胞凋亡和促进肿瘤细胞增殖、侵袭、转移[13]。有研究显示，CDC20高表达与结肠癌的临床分期、病理分化程度和TNM分期有关，且CDC20高表达结肠癌患者OS较CDC20低表达患者短[14]。本研究显示，CDC20高表达与结肠癌患者OS无统计学意义(P>0.05)，这可能与采用不同的生存分析方法有关。因此，CDC20对结肠癌患者预后的影响有待进一步研究。

COL1A1蛋白是胶原蛋白家族成员之一，是细胞外基质重要组成成分，COL1A1过表达可导致结肠癌细胞上皮-间质转化，进而促进结肠癌细胞肝转移[15]。本研究显示COL1A1高表达与结肠癌患者无病生存期呈负相关，与丁志祥等[4]研究结果一致。AURKA蛋白在细胞有丝分裂过程中起重要作用，是调节细胞周期关键分子。研究显示，使用AURKA蛋白抑制剂Alisertib可显著抑制MYC驱动的结肠癌细胞增殖[16]，也能促进KRAS驱动的结肠癌细胞死亡[17]。本研究显示AURKA过表达能够促进结肠癌细胞增殖，但AURKA高表达结肠癌患者OS却更长，这可能与AURKA高表达增加结肠癌细胞化疗敏感性有关[18]。AURKA在结肠癌中呈高表达，可作为结肠癌的预后标志物[19]。TIMP1蛋白是组织金属蛋白酶抑制剂-1，在抑制金属蛋白酶介导的细胞外基质转化方面发挥重要作用，进而参与肿瘤的侵袭和转移。研究表明，TIMP1敲除后能够通过FAK-PI3K/AKT 和 MAPK信号通路抑制结肠癌细胞的增殖和转移，且与TNM分期、无病生存期、血管侵犯以及淋巴结转移有关[20]。本研究通过生物信息学发现TIMP1表达量越高的结肠癌患者OS越短，MTT实验显示结肠癌细胞中过表达TIMP1能够显著促进细胞增殖，该结果与相关报道[20]一致，故认为TIMP1蛋白可能为潜在的结肠癌分子标志物。

综上所述，本研究通过生物信息学方法筛选出可能的结肠癌核心基因，包括MYC、TIMP1、UBE2C、CAV1、SOX9、CXCL12、AURKA、COL1A1、CDC20、TOP2A和CENPF，其中AURKA和TIMP1表达变化与结肠癌患者预后相关，但仍需更多实验研究进行验证。