基于蟻群算法的腫瘤驅動通路搜索方法研究

潘 蕾， 秦 攀， 顧 宏

( 大連理工大學 電子信息與電氣工程學部, 遼寧 大連 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)