潘 蕾， 秦 攀， 顧 宏

( 大連理工大學(xué) 電子信息與電氣工程學(xué)部, 遼寧 大連 116024 )

0 引 言

隨著人們對(duì)腫瘤機(jī)制的深入研究，發(fā)現(xiàn)腫瘤的發(fā)生和發(fā)展與細(xì)胞信號(hào)通路有著密切的關(guān)系．細(xì)胞內(nèi)癌基因突變的累積導(dǎo)致各種信號(hào)通路的紊亂，從而影響細(xì)胞的增殖、分化和凋亡，最終引起腫瘤的發(fā)生[1-2]．因此，信號(hào)通路的搜索不僅可以對(duì)腫瘤的形成機(jī)制有進(jìn)一步的了解，還為腫瘤的治療提供了新的分子靶點(diǎn)，具有重要的研究意義．

目前對(duì)于信號(hào)通路的研究方法主要有兩類．一類是針對(duì)信號(hào)通路中基因集的覆蓋性和排他性這兩個(gè)特性．例如，Vandin等[3]首先提出了最大權(quán)重子矩陣模型，在提高覆蓋性的同時(shí)抑制重疊，保證了排他性．該方法使用MCMC優(yōu)化算法雖然提高了數(shù)據(jù)適應(yīng)度，但隨機(jī)搜索的迭代方式也帶來了容易產(chǎn)生局部最優(yōu)解的問題．Leiserson等[4]對(duì)Vandin的方法進(jìn)行了改進(jìn)，提出了同時(shí)識(shí)別多個(gè)滿足覆蓋性和排他性通路的模型，并使用線性規(guī)劃的優(yōu)化思想來提高算法精度和收斂速度．但該方法并沒有解決模型的NP難題．而且，上述方法都只在本地突變數(shù)據(jù)的基礎(chǔ)上研究信號(hào)通路的兩個(gè)共性特征，沒有考慮基因自身的突變異質(zhì)性．另一類是基于信號(hào)通路的先驗(yàn)知識(shí)．Babur等[5]使用通路數(shù)據(jù)庫中已知的基因集數(shù)據(jù)，提出了一個(gè)用來量化基因間互斥性的度量，以此來構(gòu)建信號(hào)網(wǎng)絡(luò)用于篩選驅(qū)動(dòng)通路基因集．由于目前信號(hào)通路先驗(yàn)知識(shí)數(shù)據(jù)庫還不夠完善，信號(hào)網(wǎng)絡(luò)尚不能明確地表明基因間的相互作用關(guān)系．針對(duì)以上方法的缺陷，本文在不需要信號(hào)通路先驗(yàn)知識(shí)的基礎(chǔ)上，將基因復(fù)制時(shí)間作為影響突變頻率的重要協(xié)變量，重新定義最大權(quán)重子矩陣函數(shù)．在保證高覆蓋性和高排他性的同時(shí)，充分考慮基因本身協(xié)變量對(duì)突變頻率的影響，從而在更大程度上搜尋到癌癥中的信號(hào)通路．

1 模型建立及優(yōu)化

1.1 基因協(xié)變量與基因突變頻率的相關(guān)性

基因突變異質(zhì)性是腫瘤的特征之一，主要體現(xiàn)在兩個(gè)方面．第一，不同腫瘤類型之間基因突變類型的異質(zhì)性．例如，肺癌患者體內(nèi)基因突變多為C→A類型的突變，而胃腸道腫瘤的突變則多為C→G類型的突變[6]．第二，腫瘤內(nèi)部基因組的區(qū)域異質(zhì)性，即基因不同時(shí)，突變頻率也有很大的區(qū)別[7]．而造成區(qū)域異質(zhì)性的基因協(xié)變量主要有3個(gè)：基因表達(dá)水平、復(fù)制時(shí)間和染色體狀態(tài)[8-9]．這3個(gè)協(xié)變量的影響，會(huì)導(dǎo)致每個(gè)基因發(fā)生突變的頻率有所不同．

本文中，針對(duì)這3個(gè)基因協(xié)變量相互之間的關(guān)系，和它們對(duì)基因突變頻率的影響程度進(jìn)行了數(shù)值實(shí)驗(yàn)分析．結(jié)果如圖1所示，其中基因協(xié)變量數(shù)據(jù)來自于腫瘤基因組圖譜數(shù)據(jù)庫(The Cancer Genome Atlas，TCGA)，數(shù)據(jù)信息詳見本文2.1節(jié)．該協(xié)變量數(shù)據(jù)已經(jīng)證明適用于其他癌癥實(shí)驗(yàn)分析[10]．

圖1 基因協(xié)變量與基因突變頻率的相關(guān)性Fig.1 The relationship between gene covariates and gene mutation frequency

根據(jù)圖1，發(fā)現(xiàn)基因表達(dá)水平和染色體狀態(tài)與基因突變頻率之間均呈負(fù)相關(guān)關(guān)系，相關(guān)系數(shù)分別為-0.216 226 1和-0.275 197 7．基因復(fù)制時(shí)間和基因突變頻率呈正相關(guān)關(guān)系，相關(guān)系數(shù)為0.352 776 1．而且，3個(gè)協(xié)變量之間的互相關(guān)關(guān)系圖顯示，各個(gè)基因協(xié)變量之間也存在著很強(qiáng)的關(guān)聯(lián)性．相關(guān)研究表明，基因組不同區(qū)域的復(fù)制時(shí)間與基因表達(dá)水平、染色體狀態(tài)有著密切的關(guān)系．復(fù)制時(shí)間較早的基因，染色體高度螺旋，基因表達(dá)水平較高．而復(fù)制時(shí)間較晚的基因，染色體狀態(tài)疏松，基因表達(dá)水平較低，甚至不表達(dá)[11]．因此，為了減少算法的復(fù)雜度，選擇基因復(fù)制時(shí)間作為對(duì)基因突變頻率影響最為重要的協(xié)變量而結(jié)合到本文方法中．

1.2 方法介紹

細(xì)胞信號(hào)通路中的基因集有兩個(gè)特性，即高覆蓋性和高排他性[12]．其中高覆蓋性是指一個(gè)通路中的基因應(yīng)該盡可能多地覆蓋樣本，高排他性則是指每個(gè)通路中基因的突變對(duì)于每個(gè)病人來說，要盡可能呈現(xiàn)唯一性．根據(jù)通路中基因集的這兩個(gè)特性，Vandin等[3]提出了Dendrix(De novo driver exclusivity)方法，該方法的缺陷是忽略了基因自身協(xié)變量對(duì)突變頻率的影響．為了解決上述問題，本文基于基因協(xié)變量的影響，提出了一種改進(jìn)方法(ACO covariant driver pathway，ACDP)．方法模型建立步驟如下：

(1)構(gòu)造突變矩陣Am×n，m是獨(dú)立病人樣本編號(hào)，n是基因名，如圖2所示．a(chǎn)ij=1表明第i個(gè)病人的第j個(gè)基因發(fā)生了突變．

圖2 突變矩陣Fig.2 Mutation matrix

(2)定義基于基因協(xié)變量影響的最大權(quán)重子矩陣函數(shù)：

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

(1)

上述模型也可轉(zhuǎn)化為一個(gè)二元線性規(guī)劃問題進(jìn)行求解：

(2)

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

式中：N表示M矩陣中基因集的元素個(gè)數(shù)；xi∈{0,1}代表第i個(gè)病人樣本落入M矩陣中的基因是否發(fā)生了突變，突變?yōu)?，否則記為0；yj∈{0,1}代表第j個(gè)基因是否落入M矩陣，落入則為1，否則記為0；x與y分別是由xi與yj構(gòu)成的向量．

1.3 目標(biāo)函數(shù)優(yōu)化

由于上述組合優(yōu)化是一個(gè)NP難題[13]，本文使用啟發(fā)式的蟻群算法來求解組合最優(yōu)化問題．蟻群算法是一種隨機(jī)優(yōu)化算法，受到自然界真實(shí)螞蟻的行為啟發(fā)而提出的模擬進(jìn)化算法，其根據(jù)螞蟻在路徑上釋放的信息素指導(dǎo)蟻群進(jìn)行路徑尋優(yōu)[14]．蟻群算法不僅可以解決局部最優(yōu)問題，而且其后期快速收斂的特性適合解決這種大規(guī)模數(shù)據(jù)的優(yōu)化．本文提出的優(yōu)化目標(biāo)函數(shù)實(shí)際上是一個(gè)0-1背包問題，即當(dāng)限制背包的承重時(shí)，尋找使得背包中總價(jià)值最大物品問題．針對(duì)細(xì)胞信號(hào)通路搜索問題，將每個(gè)基因的“質(zhì)量”w設(shè)置為1，用承重值控制基因集中基因個(gè)數(shù)N；每個(gè)基因的“價(jià)值”c是在協(xié)變量影響下的基因突變次數(shù)；把基因是否落入限制大小的基因集中描述為某個(gè)物品是否裝入限定承重的背包．

本文使用蟻群算法對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，當(dāng)某個(gè)基因上累積的信息素越來越多時(shí)，這個(gè)基因最終落入結(jié)果基因集的概率就越大．在一次迭代中，蟻群中的每只螞蟻都是按基因選擇概率的大小來決定要選擇的基因．下式表示第k只螞蟻對(duì)基因g的選擇概率：

(3)

式中：T(k)是禁忌表，是第k只螞蟻在一次迭代中選擇基因的歷史記錄表，作用是避免重復(fù)選擇已落入基因集的基因；τg(t)是基因g在t次迭代過程中的信息素強(qiáng)度；ηg(t)是啟發(fā)函數(shù)，在解決背包問題時(shí)，常令ηg(t)=cg/wg，cg是基因g的“價(jià)值”，wg是基因g的“質(zhì)量”，則ηg代表基因g的“單位價(jià)值”；α是信息啟發(fā)因子，決定信息素的重要性；β是期望啟發(fā)因子，決定啟發(fā)函數(shù)的重要性．

每只螞蟻在選擇一個(gè)基因后，需要判斷此時(shí)背包的質(zhì)量是否超出承重值，也就是判斷已選擇基因的個(gè)數(shù)是否超出設(shè)置的基因集大小N．當(dāng)蟻群中每只螞蟻都完成一次迭代中的所有選擇后，每個(gè)基因上累積的信息素要進(jìn)行一次調(diào)整，調(diào)整公式為

(4)

(5)

式中：Q表示信息強(qiáng)度，是一個(gè)常數(shù)；gk表示第k只螞蟻在本次迭代中選擇的基因列表；ck是第k只螞蟻所選基因的“總價(jià)值”．本文中實(shí)驗(yàn)參數(shù)設(shè)置分別為α=2，β=5，ρ=0.5，Q=100，蟻群規(guī)模為30．

當(dāng)完成所有迭代過程后，計(jì)算每只螞蟻的背包價(jià)值，價(jià)值最大的背包中對(duì)應(yīng)的基因則為選擇的信號(hào)通路基因集．

1.4 ACDP方法的偽代碼實(shí)現(xiàn)

采用Matlab實(shí)現(xiàn)本文ACDP方法，表1是對(duì)程序中部分變量的解釋．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 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)集

2.2 實(shí)驗(yàn)結(jié)果對(duì)比

本文實(shí)驗(yàn)是在肺腺癌突變數(shù)據(jù)上分別運(yùn)行ACDP、Dendrix[3]、Multi-Dendrix[4]和Mutex[5]，并對(duì)4種方法的通路搜索結(jié)果在樣本中的覆蓋性、在肺腺癌相關(guān)通路中存在性，以及互斥性進(jìn)行了對(duì)比分析．其中，基因?qū)Φ幕コ庑杂少M(fèi)希爾精確檢驗(yàn)得到的P值來度量，P值越小，互斥性越顯著．對(duì)于不同的基因集大小N(4～10)，實(shí)驗(yàn)結(jié)果精度均在78%以上，以N取9為例，實(shí)驗(yàn)結(jié)果如表2所示．其中括號(hào)中表示顯著性很高的基因?qū)Γ哟植糠直硎驹摶蛟诜蜗侔┫嚓P(guān)信號(hào)通路中．

表2 4種方法實(shí)驗(yàn)結(jié)果Tab.2 The experimental results of four methods

表2結(jié)果顯示，本文方法比另外3種方法找到了更多在肺腺癌通路中的癌基因，其中包括NF1、STK11、APC等基因．很多相關(guān)的醫(yī)學(xué)研究表明，這些基因在肺腺癌中有著很高的突變頻率，會(huì)影響機(jī)體對(duì)細(xì)胞增殖、分化和凋亡過程的正常調(diào)控，從而導(dǎo)致腫瘤的發(fā)生[16-20]．更重要的是，本文方法可以直接搜索到多個(gè)互斥基因?qū)?，而Dendrix方法則需要在數(shù)據(jù)中刪除找到的基因?qū)Γ拍芩阉髌渌幕驅(qū)?，這樣做忽略了已搜索的互斥基因?qū)推渌蛑g的互斥性．如表2所示，在搜索到的基因集結(jié)果中，有兩個(gè)互斥性較為顯著的基因?qū)?EGFR KRAS)和(TP53 ATM)．醫(yī)學(xué)研究人員通過實(shí)驗(yàn)在兩個(gè)基因?qū)χ邪l(fā)現(xiàn)了很明顯的負(fù)相關(guān)性，充分說明了基因?qū)χ腥我庖粋€(gè)基因發(fā)生突變，足以使得相關(guān)通路所控制的生物功能失控．另一方面，也指出了這兩個(gè)基因?qū)Ψ謩e在MAPK信號(hào)通路和細(xì)胞循環(huán)通路中有著直接影響的作用關(guān)系[21]．如ATM基因是一個(gè)重要的抑癌基因，負(fù)責(zé)細(xì)胞循環(huán)檢測(cè)點(diǎn)酶蛋白的編碼，參與細(xì)胞周期調(diào)控．有實(shí)驗(yàn)證據(jù)顯示DNA損失程度信號(hào)是通過ATM來傳遞給下游的TP53基因，再由TP53來修復(fù)受損的DNA，促進(jìn)癌細(xì)胞的凋亡，發(fā)揮其抑癌基因的作用[22]．本文還對(duì)兩個(gè)互斥性顯著基因?qū)Φ母采w性和在細(xì)胞信號(hào)通路中作用關(guān)系做了可視化表達(dá)，更形象地體現(xiàn)本文方法結(jié)果的可靠性及生物意義，結(jié)果分別如圖3、4所示．

(a) 基因?qū)?EGFR KRAS)的覆蓋性55%(90/163)

(b) 基因?qū)?TP53 ATM)的覆蓋性47%(76/163)

圖3 基因?qū)Φ母采w性統(tǒng)計(jì)

Fig.3 The coverage statistic of gene pairs

圖4 肺腺癌信號(hào)通路Fig.4 Signal pathways in lung adenocarcinoma

3 結(jié) 語

本文提出了一種新的用于細(xì)胞信號(hào)通路中基因集搜索的方法．在考慮基因集覆蓋性和排他性的基礎(chǔ)上，用基因復(fù)制時(shí)間作為影響基因突變頻率的權(quán)重協(xié)變量，結(jié)合到方法當(dāng)中．由于本文提出的目標(biāo)函數(shù)屬于NP難題，使用蟻群算法進(jìn)行優(yōu)化．實(shí)驗(yàn)結(jié)果表明，方法不僅在癌癥相關(guān)通路中搜索到更多的癌基因，更重要的是比現(xiàn)有方法能找到更多的已經(jīng)證實(shí)在細(xì)胞信號(hào)通路中存在直接作用的互斥基因?qū)Γ?/p>

[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]高海昌,馮博琴,朱 利. 智能優(yōu)化算法求解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]丁建立,陳增強(qiáng),袁著祉. 基于自適應(yīng)螞蟻算法的動(dòng)態(tài)最優(yōu)路由選擇[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á)及意義[D]. 南充:川北醫(yī)學(xué)院, 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)