基于蚁群算法的肿瘤驱动通路搜索方法研究

潘 蕾， 秦 攀， 顾 宏

( 大连理工大学 电子信息与电气工程学部, 辽宁 大连 116024 )

0 引 言

随着人们对肿瘤机制的深入研究，发现肿瘤的发生和发展与细胞信号通路有着密切的关系．细胞内癌基因突变的累积导致各种信号通路的紊乱，从而影响细胞的增殖、分化和凋亡，最终引起肿瘤的发生[1-2]．因此，信号通路的搜索不仅可以对肿瘤的形成机制有进一步的了解，还为肿瘤的治疗提供了新的分子靶点，具有重要的研究意义．

目前对于信号通路的研究方法主要有两类．一类是针对信号通路中基因集的覆盖性和排他性这两个特性．例如，Vandin等[3]首先提出了最大权重子矩阵模型，在提高覆盖性的同时抑制重叠，保证了排他性．该方法使用MCMC优化算法虽然提高了数据适应度，但随机搜索的迭代方式也带来了容易产生局部最优解的问题．Leiserson等[4]对Vandin的方法进行了改进，提出了同时识别多个满足覆盖性和排他性通路的模型，并使用线性规划的优化思想来提高算法精度和收敛速度．但该方法并没有解决模型的NP难题．而且，上述方法都只在本地突变数据的基础上研究信号通路的两个共性特征，没有考虑基因自身的突变异质性．另一类是基于信号通路的先验知识．Babur等[5]使用通路数据库中已知的基因集数据，提出了一个用来量化基因间互斥性的度量，以此来构建信号网络用于筛选驱动通路基因集．由于目前信号通路先验知识数据库还不够完善，信号网络尚不能明确地表明基因间的相互作用关系．针对以上方法的缺陷，本文在不需要信号通路先验知识的基础上，将基因复制时间作为影响突变频率的重要协变量，重新定义最大权重子矩阵函数．在保证高覆盖性和高排他性的同时，充分考虑基因本身协变量对突变频率的影响，从而在更大程度上搜寻到癌症中的信号通路．

1 模型建立及优化

1.1 基因协变量与基因突变频率的相关性

基因突变异质性是肿瘤的特征之一，主要体现在两个方面．第一，不同肿瘤类型之间基因突变类型的异质性．例如，肺癌患者体内基因突变多为C→A类型的突变，而胃肠道肿瘤的突变则多为C→G类型的突变[6]．第二，肿瘤内部基因组的区域异质性，即基因不同时，突变频率也有很大的区别[7]．而造成区域异质性的基因协变量主要有3个：基因表达水平、复制时间和染色体状态[8-9]．这3个协变量的影响，会导致每个基因发生突变的频率有所不同．

本文中，针对这3个基因协变量相互之间的关系，和它们对基因突变频率的影响程度进行了数值实验分析．结果如图1所示，其中基因协变量数据来自于肿瘤基因组图谱数据库(The Cancer Genome Atlas，TCGA)，数据信息详见本文2.1节．该协变量数据已经证明适用于其他癌症实验分析[10]．

图1 基因协变量与基因突变频率的相关性Fig.1 The relationship between gene covariates and gene mutation frequency

根据图1，发现基因表达水平和染色体状态与基因突变频率之间均呈负相关关系，相关系数分别为-0.216 226 1和-0.275 197 7．基因复制时间和基因突变频率呈正相关关系，相关系数为0.352 776 1．而且，3个协变量之间的互相关关系图显示，各个基因协变量之间也存在着很强的关联性．相关研究表明，基因组不同区域的复制时间与基因表达水平、染色体状态有着密切的关系．复制时间较早的基因，染色体高度螺旋，基因表达水平较高．而复制时间较晚的基因，染色体状态疏松，基因表达水平较低，甚至不表达[11]．因此，为了减少算法的复杂度，选择基因复制时间作为对基因突变频率影响最为重要的协变量而结合到本文方法中．

1.2 方法介绍

细胞信号通路中的基因集有两个特性，即高覆盖性和高排他性[12]．其中高覆盖性是指一个通路中的基因应该尽可能多地覆盖样本，高排他性则是指每个通路中基因的突变对于每个病人来说，要尽可能呈现唯一性．根据通路中基因集的这两个特性，Vandin等[3]提出了Dendrix(De novo driver exclusivity)方法，该方法的缺陷是忽略了基因自身协变量对突变频率的影响．为了解决上述问题，本文基于基因协变量的影响，提出了一种改进方法(ACO covariant driver pathway，ACDP)．方法模型建立步骤如下：

(1)构造突变矩阵Am×n，m是独立病人样本编号，n是基因名，如图2所示．aij=1表明第i个病人的第j个基因发生了突变．

图2 突变矩阵Fig.2 Mutation matrix

(2)定义基于基因协变量影响的最大权重子矩阵函数：

W(M)=Γ(M)-ω(M)=

(1)

上述模型也可转化为一个二元线性规划问题进行求解：

(2)

xi,yj∈{0,1}；i=1,…,m，j=1,…,n

式中：N表示M矩阵中基因集的元素个数；xi∈{0,1}代表第i个病人样本落入M矩阵中的基因是否发生了突变，突变为1，否则记为0；yj∈{0,1}代表第j个基因是否落入M矩阵，落入则为1，否则记为0；x与y分别是由xi与yj构成的向量．

1.3 目标函数优化

由于上述组合优化是一个NP难题[13]，本文使用启发式的蚁群算法来求解组合最优化问题．蚁群算法是一种随机优化算法，受到自然界真实蚂蚁的行为启发而提出的模拟进化算法，其根据蚂蚁在路径上释放的信息素指导蚁群进行路径寻优[14]．蚁群算法不仅可以解决局部最优问题，而且其后期快速收敛的特性适合解决这种大规模数据的优化．本文提出的优化目标函数实际上是一个0-1背包问题，即当限制背包的承重时，寻找使得背包中总价值最大物品问题．针对细胞信号通路搜索问题，将每个基因的“质量”w设置为1，用承重值控制基因集中基因个数N；每个基因的“价值”c是在协变量影响下的基因突变次数；把基因是否落入限制大小的基因集中描述为某个物品是否装入限定承重的背包．

本文使用蚁群算法对目标函数进行优化，当某个基因上累积的信息素越来越多时，这个基因最终落入结果基因集的概率就越大．在一次迭代中，蚁群中的每只蚂蚁都是按基因选择概率的大小来决定要选择的基因．下式表示第k只蚂蚁对基因g的选择概率：

(3)

式中：T(k)是禁忌表，是第k只蚂蚁在一次迭代中选择基因的历史记录表，作用是避免重复选择已落入基因集的基因；τg(t)是基因g在t次迭代过程中的信息素强度；ηg(t)是启发函数，在解决背包问题时，常令ηg(t)=cg/wg，cg是基因g的“价值”，wg是基因g的“质量”，则ηg代表基因g的“单位价值”；α是信息启发因子，决定信息素的重要性；β是期望启发因子，决定启发函数的重要性．

每只蚂蚁在选择一个基因后，需要判断此时背包的质量是否超出承重值，也就是判断已选择基因的个数是否超出设置的基因集大小N．当蚁群中每只蚂蚁都完成一次迭代中的所有选择后，每个基因上累积的信息素要进行一次调整，调整公式为

(4)

(5)

式中：Q表示信息强度，是一个常数；gk表示第k只蚂蚁在本次迭代中选择的基因列表；ck是第k只蚂蚁所选基因的“总价值”．本文中实验参数设置分别为α=2，β=5，ρ=0.5，Q=100，蚁群规模为30．

当完成所有迭代过程后，计算每只蚂蚁的背包价值，价值最大的背包中对应的基因则为选择的信号通路基因集．

1.4 ACDP方法的伪代码实现

采用Matlab实现本文ACDP方法，表1是对程序中部分变量的解释．ACDP方法的具体伪代码如下：

表1 程序中部分变量解释Tab.1 The explanation of part of the variables in the program

输入：w，c，m，n，N，Sm，na

输出：Mv，Cm

1. 初始化t=1，T=∅，Cm=0

2. whilet≤Smdo

3. fork=1:nado

4.T(k)={gi}

5. forj=2:Ndo

8. end

9. end

11.Cm=max (Cm,Cmt)

12.Mv={Gs|Cm}

13.τg(t+1)=(1-ρ)τg(t)+Δτg(t)

14.T=∅

15.t=t+1

16. end

2 实验结果与分析

2.1 实验数据集

2.2 实验结果对比

本文实验是在肺腺癌突变数据上分别运行ACDP、Dendrix[3]、Multi-Dendrix[4]和Mutex[5]，并对4种方法的通路搜索结果在样本中的覆盖性、在肺腺癌相关通路中存在性，以及互斥性进行了对比分析．其中，基因对的互斥性由费希尔精确检验得到的P值来度量，P值越小，互斥性越显著．对于不同的基因集大小N(4～10)，实验结果精度均在78%以上，以N取9为例，实验结果如表2所示．其中括号中表示显著性很高的基因对，加粗部分表示该基因在肺腺癌相关信号通路中．

表2 4种方法实验结果Tab.2 The experimental results of four methods

表2结果显示，本文方法比另外3种方法找到了更多在肺腺癌通路中的癌基因，其中包括NF1、STK11、APC等基因．很多相关的医学研究表明，这些基因在肺腺癌中有着很高的突变频率，会影响机体对细胞增殖、分化和凋亡过程的正常调控，从而导致肿瘤的发生[16-20]．更重要的是，本文方法可以直接搜索到多个互斥基因对，而Dendrix方法则需要在数据中删除找到的基因对，才能搜索其他的基因对，这样做忽略了已搜索的互斥基因对和其他基因之间的互斥性．如表2所示，在搜索到的基因集结果中，有两个互斥性较为显著的基因对(EGFR KRAS)和(TP53 ATM)．医学研究人员通过实验在两个基因对中发现了很明显的负相关性，充分说明了基因对中任意一个基因发生突变，足以使得相关通路所控制的生物功能失控．另一方面，也指出了这两个基因对分别在MAPK信号通路和细胞循环通路中有着直接影响的作用关系[21]．如ATM基因是一个重要的抑癌基因，负责细胞循环检测点酶蛋白的编码，参与细胞周期调控．有实验证据显示DNA损失程度信号是通过ATM来传递给下游的TP53基因，再由TP53来修复受损的DNA，促进癌细胞的凋亡，发挥其抑癌基因的作用[22]．本文还对两个互斥性显著基因对的覆盖性和在细胞信号通路中作用关系做了可视化表达，更形象地体现本文方法结果的可靠性及生物意义，结果分别如图3、4所示．

(a) 基因对(EGFR KRAS)的覆盖性55%(90/163)

(b) 基因对(TP53 ATM)的覆盖性47%(76/163)

图3 基因对的覆盖性统计

Fig.3 The coverage statistic of gene pairs

图4 肺腺癌信号通路Fig.4 Signal pathways in lung adenocarcinoma

3 结 语

本文提出了一种新的用于细胞信号通路中基因集搜索的方法．在考虑基因集覆盖性和排他性的基础上，用基因复制时间作为影响基因突变频率的权重协变量，结合到方法当中．由于本文提出的目标函数属于NP难题，使用蚁群算法进行优化．实验结果表明，方法不仅在癌症相关通路中搜索到更多的癌基因，更重要的是比现有方法能找到更多的已经证实在细胞信号通路中存在直接作用的互斥基因对．

[1] WEINSTEIN I B. Addiction to oncogenes:The Achilles heal of cancer [J].Science, 2002,297(5578):63-64.

[2]Cancer Genome Atlas Research Network. Comprehensive genomic characterization defines human glioblastoma genes and core pathways [J].Nature, 2008,455(7216):1061-1068.

[3]VANDIN F, UPFAL E, RAPHAEL B J. De novo discovery of mutated driver pathways in cancer [J].GenomeResearch, 2012,22(2):375-385.

[4]LEISERSON M D M, BLOKH D, SHARAN R,etal. Simultaneous identification of multiple driver pathways in cancer [J].PLoSComputationalBiology, 2013,9(5):e1003054.

[5]BABUR Ö, GÖNEN M, AKSOY B A,etal. Systematic identification of cancer driving signaling pathways based on mutual exclusivity of genomic alterations [J].GenomeBiology, 2015,16(1):45.

[6]PLEASANCE E D, STEPHENS P J, O′MEARA S,etal. A small-cell lung cancer genome with complex signatures of tobacco exposure [J].Nature, 2010,463(7278):184-190.

[7]Cancer Genome Atlas Network. Comprehensive molecular characterization of human colon and rectal cancer [J].Nature, 2015,487(7407):330-337.

[8]PLEASANCE E D, CHEETHAM R K, STEPHENS P J,etal. A comprehensive catalogue of somatic mutations from a human cancer genome [J].Nature, 2010,463(7278):191-196.

[9]STAMATOYANNOPOULOS J A, ADZHUBEI I, THURMAN R E,etal. Human mutation rate associated with DNA replication timing [J].NatureGenetics, 2009,41(4):393-395.

[10]LAWRENCE M S, STOJANOV P, POLAK P,etal. Mutational heterogeneity in cancer and the search for new cancer-associated genes [J].Nature, 2013,499(7457):214-218.

[11]NEPH S, VIERSTRA J, STERGACHIS A B,etal. An expansive human regulatory lexicon encoded in transcription factor footprints [J].Nature, 2012,489(7414):83-90.

[12]VOGELSTEIN B, KINZLER K W. Cancer genes and the pathways they control [J].NatureMedicine, 2004,10(8):789-799.

[13]高海昌,冯博琴,朱 利. 智能优化算法求解TSP问题[J]. 控制与决策, 2006,21(3):241-247, 252.

GAO Haichang, FENG Boqin, ZHU Li. Reviews of the meta-heuristic algorithms for TSP [J].ControlandDecision, 2006,21(3):241-247, 252. (in Chinese)

[14]丁建立,陈增强,袁著祉. 基于自适应蚂蚁算法的动态最优路由选择[J]. 控制与决策, 2003,18(6):751-753, 757.

DING Jianli, CHEN Zengqiang, YUAN Zhuzhi. Dynamic optimization routing method based on ant adaptive algorithm [J].ControlandDecision, 2003,18(6):751-753,757. (in Chinese)

[15]TAKAHASHI T, NAU M M, CHIBA I,etal. p53:A frequent target for genetic abnormalities in lung cancer [J].Science, 1989,246(4929):491-494.

[16]AHRENDT S A, HU Y C, BUTA M,etal. p53 mutations and survival in stage Ⅰ non-small-cell lung cancer:results of a prospective study [J].JournaloftheNationalCancerInstitute, 2003,95(13):961-970.

[17]VELCHETI V, GOVINDAN R. Hedgehog signaling pathway and lung cancer [J].JournalofThoracicOncology, 2007,2(1):7-10.

[18]STEWART D J. Wnt signaling pathway in non-small cell lung cancer [J].JournaloftheNationalCancerInstitute, 2014,106(1):djt356.

[19]TOMASINI P, WALIA P, LABBE C,etal. Targeting the KRAS pathway in non-small cell lung cancer [J].Oncologist, 2016,21(12):1450-1460.

[20]PACKENHAM J P, TAYLOR J A, WHITE C M,etal. Homozygous deletions at chromosome 9p21 and mutation analysis of p16 and p15 in microdissected primary non-small cell lung cancers [J].ClinicalCancerResearch, 1995,1(7):687-690.

[21]DING Li, GETZ G, WHEELER D A,etal. Somatic mutations affect key pathways in lung adenocarcinoma [J].Nature, 2008,455(7216):1069-1075.

[22]唐光明. p53及其上游基因ATM、下游基因PUMA在大肠癌中的表达及意义[D]. 南充:川北医学院, 2015.

TANG Guangming. The expression and significance of p53 and upstream gene ATM, downstream gene PUMA in colorectal cancer [D]. Nanchong:North Sichuan Medical College, 2015. (in Chinese)