引言

结直肠癌(Colorectal cancer，CRC)包括结肠癌和直肠癌，是目前全球范围内较为常见的消化系统恶性肿瘤之一^[1]。近几年来，随着经济的发展、人们生活水平的提高、饮食习惯的改变以及遗传因素的影响等，我国结直肠癌的发病率和死亡率逐年上升^[2]，其病死率居世界范围内恶性肿瘤的第4位，发病率居第3位^[3]，已严重威胁人类的健康。大多数研究者认为结直肠癌是一种复杂疾病，其发生与发展是在遗传因素、环境因素以及肿瘤发展过程的个体反应等背景下，多因素参与相互作用的结果。因此，预测并鉴定CRC的致病基因对深入理解其致病机制有着重要的意义。

由多个基因相互作用所产生的疾病通常被认为是复杂疾病，而挖掘与其相关的致病基因并阐释它们在疾病发展过程中所产生的作用已逐渐成为人们研究复杂疾病的主要目标之一。随着生物信息学分析的发展，基于网络的方法已经成为研究疾病机制的有力工具^[4-6]，其主要包括共表达网络^[7]、蛋白质互作网络^[8]、蛋白质磷酸化网络^[9]和DNA甲基化网络^[10]等生物分子网络。而随着蛋白质相互作用(Protein-protein interaction, PPI)数据的丰富，基于蛋白质互作网络^[11]的方式逐渐成为挖掘复杂疾病候选基因的主要方法。其中，模块分析或聚类分析作为一种工具和算法已经广泛应用到蛋白质互作网络分析中，例如从局部扩展来探测社区的重叠社区发现算法(Local and wave-like extension algorithm of detecting overlapping community, LWS-OCD)算法^[12]能够有效地发现无向图中的重叠区域。而拓扑分析作为另外一种分析复杂网络的工具也被广泛应用，如其中的中心性分析中图中心的措施如度、介数和接近中心性对识别在网络中具有关键作用的节点十分有用。在蛋白质互作网络中，通常称度数高的节点为中心节点，介数高的节点为瓶颈节点，这些节点在网络中都发挥着举足轻重的作用。

疾病基因筛选面临的主要问题是数据高维且样本量少，因此其数据处理主要侧重于数据挖掘方法的多重筛选与验证。本研究是基于结直肠癌的公开表达谱数据，进行结直肠癌差异表达基因分析以获取疾病基因的初选；通过这些初选的疾病基因与结直肠癌已知致病基因之间的重要关系来组建互作网络，进行结直肠癌致病的候选基因筛选；最后，利用功能富集分析和子网络的拓扑分析进行结直肠癌疾病候选基因的筛选与验证。

1 实验材料与方法 1.1 材料来源

本文从美国生物技术中心GEO(http://www.ncbi.nlm.nih.gov/geo/)数据库下载结直肠癌基因表达谱数据，登录号为GSE9348，平台为Affymetrix Plus 2.0。下载其中的82例包含癌旁样本12例和癌样本70例作为样本数据，其中数据以DNA微阵列的表达谱数据形式给出，且维数是5万多。这一数据作为本研究的实验数据，记为Ω；从在线孟德尔人类遗传学数据库(Online Mendelian inheritance in man, OMIM)中下载135个结直肠癌的疾病基因，这些疾病基因是已被证实的结直肠癌致病基因，该数据在本文中当作验证数据，记为Ω₀。

1.2 方法 1.2.1 数据处理

在数据Ω的基础上，首先通过软件R对原始数据进行重复值合并、数据转换、标准化和缺失数据删除等处理。然后，利用T检验法进行基因筛选实现降维，挑选出样本间有显著性差异的基因，即筛选出结直肠癌的差异表达基因。这一过程是通过R中的LIMMA包(http://www.bioconductor.org/packages/release/bioc/html/limma.html)来实现，其中涉及参数P和Fold change。参数P为基因间差异的显著性水平，参数Fold change为实验组与对照组的差异表达的比值，且P < α和Fold change >β。一般地，参数P越小且Fold change越大，则筛选的基因就越多。在具体实验中，依据研究需要选取阈值α和β，一般分别取0.05和2。通过数据处理获得的差异表达基因记为Ω₁。

1.2.2 蛋白质互作网络的构建

将基因集Ω₀和基因集Ω₁上传到STRING(Search tool for the retrieval of interacting genes)数据库中，通过SRTING在线分析工具获得疾病基因与差异基因之间的相互作用^[13]，然后选取可靠性指数(Confidence scores)大于0.9的基因对，以及它们之间的作用数据，得到数据集Ω₂。将Ω₂导入Cytoscape软件中，对结直肠癌的已知疾病与差异表达基因之间的相互作用进行可视化^[14]，得到相互作用网络图。

1.2.3 子网络的识别

重叠领域扩展聚类(Clustering with overlap neighborhood expansion，ClusterONE)是一种从加权或非加权的蛋白质网络中挖掘重叠的密集连接区域的复合物识别算法^[15]，已有很多研究者利用其识别能力在各种生物网络中检测有意义的局部结构^[16-17]。因此，本文可通过这种图聚类算法来挖掘蛋白质互作网络中节点高度连接的重叠区域，其算法的聚类凝聚力定义为

$ f\left( V \right) = \frac{{{W^{{\rm{in}}}}\left( V \right)}}{{{W^{{\rm{in}}}}\left( V \right) + {W^{{\rm{boun}}\;{\rm{d}}}}\left( V \right) + P\left| V \right|}} $

(1)

式中：Wⁱⁿ(V)表示一组顶点V内的边的总重量；W^bound(V)表示连接此组顶点到图的其余部分边的总重量；P|V|为惩罚项。

本文通过Cytoscape的ClusterONE插件来进行子网络的挖掘，其中，涉及到参数minimum density、degree和P的设定。通常情况下，参数P越大，一般被挖掘的聚类(模块)数量越少，minimum density和degree越大，反之亦然。因此，P值的选择应在一个合理范围内，一般认为不能太小。为便于分析，被挖掘的子网络模块的数量应在适宜的范围，因此参数minimum density和degree也有阈值，一般分别不小于0.5和6。

1.2.4 子网络的拓扑分析与富集分析

在无标度网络中^[18]的中心节点通常由大量互作连线的节点来代表，而中心节点对应的蛋白质为核心蛋白质(基因)，并且这些基因在生理调节的过程中扮演着重要的角色。因此，本文通过Cytoscape的NetworkAnalyzer插件来筛选子网络中的中心节点，NetworkAnalyzer插件可对有向网络和无向网络进行拓扑分析，得到拓扑属性值如度、介数和接近中心性等，定义如下：

(1) 度(Degree)

$ \deg \left( v \right) = \left| {N\left( v \right)} \right| $

(2)

式中：v为节点，N(v)表示节点v邻节点的集合；deg(v)表示v与网络中其他节点的关联性，当deg(v)越大时，说明v在网络中越重要。

(2) 接近中心性(Closeness centrality)

$ {C_{\rm{c}}}\left( v \right) = \sum\limits_{w \in V} {\frac{1}{{{\rm{dist}}\left( {v, w} \right)}}} $

(3)

式中：dist(v, w)表示节点v, w间的最短路径长度；C_c(v)表示v接近中心位置的程度，当C_c(v)越小时，越说明v是网络的核心点，在网络中作用越重要。

(3) 介数(Betweenness centrality)

$ BC\left( v \right) = \sum\limits_{s \ne t \ne v \in C\left( v \right)} {\frac{{{\sigma _{st}}\left( v \right)}}{{{\sigma _{st}}}}} $

(4)

式中：C(v)为包含节点v的组成部分；σ_st表示从节点s到节点t的最短路径数；BC(v)表示v对网络中其他节点之间通信连接的影响程度。当BC(v)越大时，表明v的重要度越高。

CytoHubba软件(http://hub.iis.sinica.edu.tw/cytoHubba/index.html)通过网络特征对网络中的节点进行排名，它实现了11个节点的排名方法来评估生物网络中节点的重要性，包括度、边缘渗出成分、最大社区成分、最大邻域成分和最大团中心等，其中最大团中心(Maximal clique centrality，MCC)方式能够更准确地预测蛋白质互作网络中的重要基因^[19]，其定义为

$ {\rm{MCC}}\left( v \right) = \sum\limits_{C \in S\left( v \right)} {\left( {\left| C \right|-1} \right)!} $

(5)

式中：S(v)是包含节点v的最大团簇，(|C|－1)!是小于|C|的所有正整数乘积。如果节点v的邻节点间没有边界，则节点v的最大团中心就是它的度。

子网络的拓扑分析就是利用网络节点度、介数、接近中心性等指标进行差异表达基因分析比较，筛选出子网络的中心节点，即结直肠癌的致病基因。子网络富集分析的目的是用来注释这些差异表达基因参与的生物学过程^[20]，且通过选取基因本体论(Gene ontology, GO)条目来实现。在GO条目中涉及错误发现率(False discovery rate, FDR)，且FDR < γ。阈值γ越小表示GO条目样本变化判断越准确, 试验中取γ=0.05。

综上，本文所采用方法的流程为

(1) 在Ω₀基础上，利用R软件和LIMMA包进行降维和筛选，获得差异表达基因，即基因集Ω₁；

(2) 将基因集Ω₀和基因集Ω₁上传到STRING在线工具，获得可靠性大于0.9的基因对，并从中删除差异表达基因之间的基因对，得到数据集Ω₂，利用Cytoscape得到Ω₂蛋白质互作网络；

(3) 在步骤(2)的基础上，通过ClusterONE聚类算法对蛋白质互作网络进行分析，获取其子网络，记子网络的基因集为Ω₃；

(4) 在步骤(3)的基础上，进行子网络拓扑分析，列出节点度、介数、接近中心性排名前10的差异表达基因，并对这些差异表达基因进行分析比较，筛选出子网络中的中心节点，即结直肠癌的致病基因。

2 实验结果与分析 2.1 实验结果

在原始数据Ω的基础上，通过一系列的预处理和基因的筛选，共获得了339个差异基因Ω₁，其中α=0.05，β=2。

采用STRING在线分析工具可获得268个可靠性指数大于0.9的基因对，即|Ω₂|=268。利用Cytoscape对Ω₂可视化的结果如图 1所示。在图 1中，蛋白质相互作用网路由132个基因对组建。在此基础上，通过ClusterONE插件，取minimum density>0.5，degree > 6，P < 0.001，可挖掘出一个含有53个节点的子网络模块，如图 2所示, 其中图 1，2中的粉色为识别差异表达基因，青色为结直肠癌的已知疾病基因。

为了识别子网络中潜在的枢纽基因，NetworkAnalyzer插件统计网络中节点的度、介数和接近中心性，且节点度、介数和接近中心性最高的点被认为是中心节点，统计结果如表 1所示，表中所列的是子网络中节点度、介数和接近中心性分别排名前10的差异表达基因。从表 1可知，节点度、介数和接近中心性都排在前10的差异基因有FOS，CCND1，CEBPB，EGR1和NOS3。因此，FOS，CCND1，CEBPB，EGR1和NOS3在被识别的网络中具有重要作用，即为网络的中心节点。

此外，通过Cytoscape的CytoHubba插件获得了排名前15的差异表达基因，包括TP53，FOS，EP300，CCND1，PPARG，AKT1，SMAD3，SERPINE1，CEBPB，MYC，SMAD4，CDC25A，EGR1和NOS3，其中FOS，CCND1，SERPINE1，CEBPB，CDC25A，EGR1和NOS3为差异表达基因。结合筛选的中心节点，可获得结直肠癌候选疾病基因FOS，CCND1，CEBPB，EGR1，NOS3。

进一步，对含有53个节点的子网络模块进行富集分析。依据FDR < 0.05筛选出的子网络差异表达基因的9条GO条目测试，结果见表 2，其中Count记录的是子网络模块中参与相应生物学过程的基因数，FDR是检验的错误发生率。由表 2可知，子网络模块中差异基因参与大分子代谢负调控过程、基因表达调控和细胞代谢负调控过程等9个生物学过程。

2.2 结果分析

对于实验所获结果，细胞周期蛋白D1(Cyclin D1，CCND1)为细胞周期蛋白的编码基因，参与细胞周期G1-S期转换的调控，目前已被公认为是一种原癌基因^[21]。它的过度表达能够使细胞G1期缩短，促进G1期向S期过渡^[22]，从而导致DNA合成及增殖，引起细胞增殖失控和癌变。此外，已有研究表明，Cyclin D1异常表达与包括结直肠癌在内的多种肿瘤发生密切相关^[23]，例如Tan等^[24]发现CCND1 870A能够增加结直肠癌的发病风险；Toncheva等^[25]发现Cycling D1在结直肠癌中的蛋白表达阳性率明显高于正常结直肠组织。

c-fos为早期反应因子，能够对各种刺激做出迅速的反应。作为原癌基因的c-fos参与了细胞增殖、侵袭、凋亡的调控和血管生成、分化和凋亡等过程，并在肿瘤细胞的运动、细胞外基质的降解、异常粘附及转移灶新生血管生长等多个环节中有着重要作用。此外，有研究表明，c-fos异常表达对肿瘤细胞的生长产生重要的影响^[26]，例如调节肿瘤细胞的侵袭性生长^[27]。

EGR1为早期生长反应因子1，是一种对细胞生长起调节作用的转录因子。在细胞的增殖、分化、凋亡、信号传递和细胞生长调控等过程中产生了重要影响，并且能够促进细胞增生和组织修复。近年来，有研究表明外源性Egr-1基因有抑制肿瘤细胞的生长和恶性转化的作用，其异常表达可引起细胞增殖与凋亡的失衡，对乳腺癌、食管癌、胃癌、前列腺癌等癌细胞有着十分重要的影响。

NOS3为人血管内皮型一氧化氮合酶(endothelial Nitric oxide synthase，eNOS)，其释放的内皮型一氧化氮在血管生成过程中能够促进血管舒张，调控血管的生成。事实上，实体肿瘤的固有特征之一就是异常的血管增生，而肿瘤细胞或激活的免疫细胞都能产生促血管因子。例如，巨噬细胞、中性粒细胞、肥大细胞以及骨髓造血祖细胞等参与固有免疫的细胞在血管生成转换的过程中都能产生重要的作用。促进肿瘤细胞的增殖和迁移的NO/NOS3信号转导同时可以增强血管通透性、诱导细胞外基质降解以及肿瘤血管新生^[28]。目前，已有多个研究发现NOS 3 894G>T增加了膀胱癌、结直肠癌和乳腺癌的发病风险^[29]。

转录因子CEBPB(CCAAT enhancer-binding protein，CCAAT增强结合蛋白)，参与许多生物过程，包括细胞分化、代谢平衡、增殖，肿瘤发生、凋亡以及机体的免疫、应激反应，能量代谢和血液生成^[30]。有研究表明C/EBP在肿瘤发生发展中通过与其他基因相互作用形成调控网络而影响肿瘤发生^[31]，也有研究显示C/EBP转录因子的LIP亚型可以诱导人乳腺癌细胞凋亡并诱发其自身的吞噬作用，分析可能是一种诱发肿瘤自我吞噬的抑癌因袭。此外，PLAC1基因在人类大部分恶性肿瘤中都是异常表达。例如，在乳腺癌组织中，PLAC1能够被C/EBP激活，并在乳腺癌的恶性转化过程中产生重要的作用^[32]。另外，通过检测得到C/EBP在肝癌、卵巢癌、淋巴细胞白血病、淋巴瘤、皮肤癌等多种肿瘤中均表达上调，但其在结直肠癌中具体的表达机制尚不清楚，需要进一步的深入研究。

综上，本文获得的结直肠癌致病基因FOS，CCND1，CEBPB，EGR1和NOS3在结直肠癌的发展过程中起着重要的作用。

3 结束语

本研究通过生物学分析对结直肠癌致病基因进行了识别。首先，基于GEO中GSE9348基因表达数据集，利用R语言的LIMMA包筛选出P < 0.05，Fold change>2的结直肠癌差异基因339个，并从OMIM数据库中下载结直肠癌已知的致病基因135个；其次，将这339个差异基因和135个已知的致病基因上传到STRING数据库中，获得了228个可靠性指数大于0.9的基因对，并通过这些基因对构建了差异表达基因与致病基因的蛋白质互作网络；进一步，利用Cytoscape软件的ClusterONE插件进行了蛋白质互作网络模块分析，获得了一个含有53个基因的子网络；最后，通过对子网络的拓扑分析，获得了FOS，CCND1，CEBPB，EGR1和NOS3等5个新结直肠癌致病基因。同时，通过功能富集分析和文献挖掘对新发现的致病基因进行了验证，数据试验结果显示本文的研究方法是行之有效的。本文提供的方法具有通用性，它对癌症发病机制的阐述以及分子靶向寻找等研究具有重要意义。