謝雨絲,盧吉曉,黃毅然,鐘 誠(chéng)

1(廣西大學(xué) 計(jì)算機(jī)與電子信息學(xué)院,南寧 530004)

2(廣西醫(yī)科大學(xué)附屬腫瘤醫(yī)院 放療技術(shù)中心,南寧 530021)

1 引 言

代謝路徑功能上是細(xì)胞內(nèi)參與特定功能的一組反應(yīng)或酶,給定起始化合物和目標(biāo)化合物,代謝路徑預(yù)測(cè)是使用計(jì)算方法在代謝數(shù)據(jù)庫(kù)中預(yù)測(cè)和分析從起始化合物到目標(biāo)化合物的新的替代路徑[1,2].人們可以通過(guò)這些替代路徑發(fā)現(xiàn)具有應(yīng)用價(jià)值的新的合成代謝路徑[3,4].

代謝路徑預(yù)測(cè)方法尋找新替代路徑的一種主要方式是以化合物為節(jié)點(diǎn),以連接化合物的代謝反應(yīng)為邊構(gòu)建代謝網(wǎng)絡(luò)圖,并在代謝網(wǎng)絡(luò)圖中搜索連接起始化合物與目標(biāo)化合物的路徑[5,6].由于化合物和代謝反應(yīng)種類繁多,由它們構(gòu)成的代謝網(wǎng)絡(luò)規(guī)模巨大,代謝路徑預(yù)測(cè)方法在代謝網(wǎng)絡(luò)中兩節(jié)點(diǎn)間可能會(huì)找到數(shù)以千計(jì)的結(jié)果路徑[7,8],因此需要通過(guò)合適的代謝路徑排序方法對(duì)返回的結(jié)果路徑進(jìn)行排序,以使用戶可以從數(shù)量眾多的路徑搜索結(jié)果中快速找到感興趣且具有生化相關(guān)性的替代路徑[9].

現(xiàn)有的代謝路徑預(yù)測(cè)方法通過(guò)選取代謝路徑的相關(guān)局部特征作為權(quán)重來(lái)為路徑排序.一些代謝路徑預(yù)測(cè)方法[10]通過(guò)代謝反應(yīng)的吉布斯自由能判斷路徑中的化學(xué)反應(yīng)是否熱力學(xué)可行,并結(jié)合路徑中的反應(yīng)自由能對(duì)結(jié)果路徑進(jìn)行排序以獲取熱力學(xué)可行性好的替代路徑[11].例如,代謝路徑預(yù)測(cè)方法Phdseeker利用結(jié)果路徑中的代謝反應(yīng)的吉布斯自由能為每一個(gè)代謝反應(yīng)評(píng)分,將所有代謝反應(yīng)評(píng)分匯總值作為結(jié)果路徑排序的因素[12].除了路徑中代謝反應(yīng)的吉布斯自由能特性,化合物對(duì)宿主的毒性可能會(huì)影響代謝路徑的產(chǎn)物產(chǎn)量,因此也常作為路徑排序的一個(gè)因素.例如,Carbonell等人設(shè)計(jì)的代謝路徑預(yù)測(cè)方法XTMS將結(jié)果路徑中的每個(gè)中間底物對(duì)宿主的毒性匯總,并用毒性匯總值作為結(jié)果路徑排序的依據(jù)[13].代謝反應(yīng)數(shù)量少的代謝路徑在代謝工程實(shí)踐中通常比代謝反應(yīng)數(shù)量多的路徑更容易合成.為此,Ravikrishnan等人設(shè)計(jì)的代謝路徑預(yù)測(cè)方法Metquest采用代謝路徑的長(zhǎng)度作為路徑排序因素,以選取代謝路徑中的反應(yīng)數(shù)量相對(duì)較少的替代路徑[7].

此外,通過(guò)追蹤由起始化合物轉(zhuǎn)移到目標(biāo)化合物的原子或原子團(tuán)來(lái)尋找代謝路徑能夠有效避免簇化合物[14].一些代謝路徑預(yù)測(cè)方法利用代謝路徑中從起始化合物到目標(biāo)產(chǎn)物間的原子或原子團(tuán)的轉(zhuǎn)移數(shù)量對(duì)代謝路徑進(jìn)行排序,以獲得生化相關(guān)性較好的路徑.例如,代謝路徑預(yù)測(cè)方法BPAT-M[15]和Retrace[16]利用代謝路徑長(zhǎng)度、代謝路徑中化合物間的碳原子轉(zhuǎn)移數(shù)量來(lái)為代謝路徑排序.BPFinder[17]和AGPathFinder[18]方法則結(jié)合反應(yīng)底物與產(chǎn)物的分子結(jié)構(gòu)相似程度、路徑中化合物間的原子團(tuán)轉(zhuǎn)移數(shù)量來(lái)為代謝路徑排序.

代謝路徑的生化相關(guān)性主要體現(xiàn)在路徑中的反應(yīng)段表現(xiàn)出的特定功能,如糖酵解路徑通過(guò)將葡萄糖轉(zhuǎn)化為丙酮酸釋放和吸收能量[19],這種功能性路徑由特定的反應(yīng)段構(gòu)成,在代謝網(wǎng)絡(luò)中經(jīng)常出現(xiàn)[20].雖然現(xiàn)有代謝路徑預(yù)測(cè)方法通過(guò)利用代謝路徑的局部生化特征,可以準(zhǔn)確選取滿足特定需求的代謝路徑,但是通過(guò)路徑的局部特征對(duì)結(jié)果路徑排序沒(méi)有考慮多個(gè)化合物之間及多個(gè)代謝反應(yīng)段之間的相互作用和聯(lián)系,不易于發(fā)現(xiàn)具有結(jié)構(gòu)功能特征的代謝路徑.針對(duì)這些問(wèn)題,本文提出了一種新的代謝路徑排序方法KPRank以提高代謝路徑預(yù)測(cè)結(jié)果質(zhì)量,主要貢獻(xiàn)如下:

1)提出代謝路徑關(guān)鍵反應(yīng)段的概念,以捕捉和度量代謝路徑的結(jié)構(gòu)功能特征.

2)提出一種融合代謝路徑關(guān)鍵反應(yīng)段的信息熵與互信息對(duì)代謝路徑進(jìn)行評(píng)分和排序的模型,以充分利用代謝路徑的結(jié)構(gòu)功能特征信息.

在KEGG代謝數(shù)據(jù)庫(kù)[21-23]上的實(shí)驗(yàn)結(jié)果表明,通過(guò)KPRank方法對(duì)代謝路徑預(yù)測(cè)算法Retrace,BPFinder和MetQuest的路徑預(yù)測(cè)結(jié)果進(jìn)行排序可以發(fā)現(xiàn)生化相關(guān)性更好的代謝路徑,有效提高了代謝路徑預(yù)測(cè)結(jié)果的質(zhì)量.

2 方 法

2.1 問(wèn)題定義

代謝路徑排序,即對(duì)代謝路徑預(yù)測(cè)算法的路徑搜索結(jié)果中的代謝路徑按路徑生化相關(guān)性特征或用戶感興趣的特征對(duì)路徑進(jìn)行量化評(píng)分和排序,為用戶快速準(zhǔn)確提供生化相關(guān)性好的替代路徑.代謝路徑排序方法之間的差異主要是為路徑排序所選取的路徑特征不同,以及路徑特征的量化評(píng)分方式的不同.在代謝路徑排序中選擇合適的路徑排序特征,設(shè)計(jì)合理的路徑特征的量化評(píng)分方法,是準(zhǔn)確地發(fā)現(xiàn)生化相關(guān)的替代路徑的關(guān)鍵.

2.2 已知代謝路徑集合構(gòu)建

本文通過(guò)KEGG的PATHWAY數(shù)據(jù)庫(kù)[24]的67個(gè)已知代謝網(wǎng)絡(luò)來(lái)建立已知代謝路徑集合P,以便從已知代謝路徑集合P中選取合適的路徑排序特征對(duì)代謝路徑進(jìn)行排序.

對(duì)于一個(gè)已知代謝網(wǎng)絡(luò),首先枚舉該代謝網(wǎng)絡(luò)的所有化合物的成對(duì)組合,以每個(gè)化合物對(duì)組合中的兩個(gè)化合物分別作為起始化合物與目標(biāo)化合物,在這個(gè)代謝網(wǎng)絡(luò)里采用Yen等人設(shè)計(jì)的最短路徑搜索算法[25]搜索連接這對(duì)化合物的前5條線性路徑,然后將所有找到的線性路徑保存到已知代謝路徑集合P里.接下來(lái)通過(guò)同樣的方式從67個(gè)已知代謝網(wǎng)絡(luò)中的每一個(gè)代謝網(wǎng)絡(luò)中提取線性路徑并保存到已知代謝路徑集合P里,并去除已知代謝路徑集合P中重復(fù)的路徑.

圖1展示了已知代謝路徑集合構(gòu)建過(guò)程的一個(gè)例子,其中cn1、cn2、cn3和cn4是代謝網(wǎng)絡(luò)rn1中的化合物節(jié)點(diǎn),R20～R25是連接化合物節(jié)點(diǎn)的代謝反應(yīng).

圖1 構(gòu)建已知代謝路徑集合示例

由圖1可以看到,在構(gòu)建已知代謝路徑集合時(shí),首先選出代謝網(wǎng)絡(luò)rn1的所有化合物cn1、cn2、cn3和cn4,枚舉它們的成對(duì)組合,并以每對(duì)化合物組合中的兩個(gè)化合物分別作為起始化合物和目標(biāo)化合物,然后在代謝網(wǎng)絡(luò)rn1中通過(guò)Yen算法[25]搜索連接這對(duì)化合物的前5條線性路徑,最后將找到的線性路徑存入已知代謝路徑集合,同時(shí)去掉重復(fù)的代謝路徑.

2.3 關(guān)鍵反應(yīng)段

特定功能路徑是表現(xiàn)出特定生化功能的連續(xù)固定的代謝反應(yīng)組合,它們通常大量出現(xiàn)在不同物種的代謝網(wǎng)絡(luò)中[20].例如,三羧酸循環(huán)路徑是各種需氧生物體內(nèi)普遍存在的連接糖類、脂類及氨基酸類化合物的代謝路徑[19].為更充分的挖掘路徑反應(yīng)的結(jié)構(gòu)功能特征,本文提出代謝路徑關(guān)鍵反應(yīng)段KRS(Key Reaction Set)概念,以捕捉和度量代謝路徑的反應(yīng)結(jié)構(gòu)功能.

對(duì)于已知代謝路徑集合P中的每一條代謝路徑,每個(gè)反應(yīng)段可由路徑中的連續(xù)2個(gè)、3個(gè)或者4個(gè)反應(yīng)構(gòu)成.可以統(tǒng)計(jì)每個(gè)反應(yīng)段在已知代謝路徑集合P中出現(xiàn)的次數(shù).關(guān)鍵反應(yīng)段是在已知代謝路徑集合P中出現(xiàn)次數(shù)大于100次的反應(yīng)段.關(guān)鍵反應(yīng)段通常是已知代謝網(wǎng)絡(luò)中高頻出現(xiàn)的連續(xù)反應(yīng)組合,這些反應(yīng)組合前后的代謝反應(yīng)搭配豐富,反應(yīng)組合間的關(guān)聯(lián)程度較高.

本文引入兩個(gè)評(píng)估關(guān)鍵反應(yīng)段的指標(biāo)以評(píng)估關(guān)鍵反應(yīng)段中反應(yīng)間的關(guān)聯(lián)程度:關(guān)鍵反應(yīng)段的信息熵和互信息值.信息熵是某條消息所含的信息量,它反映的是獲知某個(gè)消息后,關(guān)于該事件的不確定性的減少量[26].對(duì)已知代謝路徑集合P,關(guān)鍵反應(yīng)段的左反應(yīng)是關(guān)鍵反應(yīng)段最左邊的反應(yīng)在關(guān)鍵反應(yīng)段出現(xiàn)的路徑p∈P里的前一個(gè)反應(yīng),關(guān)鍵反應(yīng)段的右反應(yīng)是關(guān)鍵反應(yīng)段最右邊的反應(yīng)在關(guān)鍵反應(yīng)段出現(xiàn)的路徑p∈P里的后一個(gè)反應(yīng).設(shè)Se是由m個(gè)連續(xù)反應(yīng)組成的關(guān)鍵反應(yīng)段,L是關(guān)鍵反應(yīng)段Se的左反應(yīng)集合,Y是關(guān)鍵反應(yīng)段Se的右反應(yīng)集合,關(guān)鍵反應(yīng)段Se的左反應(yīng)信息熵Hl(Se)和右反應(yīng)信息熵Hy(Se)分別可以通過(guò)式(1)和式(2)計(jì)算[26]:

Hl(Se)=-∑x∈Lf(x)log2f(x)

(1)

Hy(Se)=-∑y∈Yf(y)log2f(y)

(2)

式(1)和式(2)中的f(x)和f(y)分別是左反應(yīng)x∈L和右反應(yīng)y∈Y在已知代謝路徑集合P里出現(xiàn)的頻率.關(guān)鍵反應(yīng)段的左、右信息熵越大,說(shuō)明該關(guān)鍵反應(yīng)段前后反應(yīng)搭配越豐富,通過(guò)關(guān)鍵反應(yīng)段信息熵可以評(píng)估關(guān)鍵反應(yīng)段的外部反應(yīng)搭配的豐富程度.互信息是一種度量離散型隨機(jī)變量相關(guān)程度的方法.為了進(jìn)一步評(píng)估關(guān)鍵反應(yīng)段內(nèi)部多個(gè)反應(yīng)之間的相關(guān)程度,本文采用改進(jìn)的互信息計(jì)算方法[26]來(lái)計(jì)算關(guān)鍵反應(yīng)段的互信息值以評(píng)估關(guān)鍵反應(yīng)段內(nèi)部反應(yīng)間的關(guān)聯(lián)程度:

(3)

(4)

(5)

其中MMI(Se)是關(guān)鍵反應(yīng)段Se的互信息值,Count(Se)是關(guān)鍵反應(yīng)段Se在已知代謝路徑集合P中出現(xiàn)的次數(shù),fr(Se)是關(guān)鍵反應(yīng)段Se出現(xiàn)在已知代謝路徑集合P中的概率,N表示已知代謝路徑集合P的路徑數(shù)量,Avg(Se)是關(guān)鍵反應(yīng)段Se中不同的反應(yīng)組合在已知代謝路徑集合P中的平均出現(xiàn)概率.通過(guò)關(guān)鍵反應(yīng)段的互信息值MMI(Se)可以反映出關(guān)鍵反應(yīng)段Se內(nèi)部多個(gè)反應(yīng)之間的相關(guān)程度.關(guān)鍵反應(yīng)段Se出現(xiàn)在代謝路徑集合P中的概率越高,同時(shí)關(guān)鍵反應(yīng)段Se中不同的反應(yīng)組合在已知代謝路徑集合P中出現(xiàn)越少,MMI(Se)也越高,說(shuō)明關(guān)鍵反應(yīng)段Se內(nèi)部的反應(yīng)組合的關(guān)聯(lián)程度也較高.

在文本信息處理中,互信息可以用于句子分詞,將句子中的字符視為隨機(jī)變量,利用字符在文本句子中出現(xiàn)的概率計(jì)算該字符與前后字符之間的互信息值來(lái)度量這些字符之間的關(guān)聯(lián)程度[27].類似地,在本文的代謝路徑的反應(yīng)關(guān)聯(lián)分析中,將代謝路徑中的代謝反應(yīng)視為隨機(jī)變量,通過(guò)利用反應(yīng)在代謝路徑集合中的代謝路徑里出現(xiàn)的概率計(jì)算該反應(yīng)與前后反應(yīng)之間的互信息值來(lái)度量這些反應(yīng)之間的關(guān)聯(lián)程度.文本信息處理和代謝反應(yīng)關(guān)聯(lián)分析中的互信息計(jì)算,都是在特定數(shù)據(jù)集內(nèi)根據(jù)隨機(jī)變量的出現(xiàn)概率來(lái)度量隨機(jī)變量之間的關(guān)聯(lián)程度,與隨機(jī)變量的具體形式是字符或代謝反應(yīng)實(shí)際上并沒(méi)有關(guān)系.

2.4 基于關(guān)鍵反應(yīng)段的代謝路徑評(píng)分方法

綜合關(guān)鍵反應(yīng)段Se的信息熵和互信息值,可以得到代謝路徑的關(guān)鍵反應(yīng)段Se的評(píng)分方法:

Scorekr(Se)=H(Se)+I(Se)

(6)

(7)

(8)

其中Score_kr(Se)是關(guān)鍵反應(yīng)段Se的評(píng)分.H(Se)是關(guān)鍵反應(yīng)段Se的左信息熵和右信息熵的綜合值,其中Hl(Se)是通過(guò)式(1)計(jì)算得到的關(guān)鍵反應(yīng)段Se的左信息熵,Hy(Se)是通過(guò)式(2)計(jì)算得到的關(guān)鍵反應(yīng)段Se的右信息熵.minHl表示已知代謝路徑集合P的所有關(guān)鍵反應(yīng)段的左信息熵的最小值,maxHl表示已知代謝路徑集合P的所有關(guān)鍵反應(yīng)段的左信息熵的最大值.minHy表示已知代謝路徑集合P的所有關(guān)鍵反應(yīng)段的右信息熵的最小值,maxHy表示已知代謝路徑集合P的所有關(guān)鍵反應(yīng)段的右信息熵的最大值.I(Se)表示關(guān)鍵反應(yīng)段Se的歸一化互信息值,其中minMMI是已知代謝路徑集合P的所有關(guān)鍵反應(yīng)段的互信息值的最小值,maxMMI是已知代謝路徑集合P的所有關(guān)鍵反應(yīng)段的互信息值的最大值.

式(6)結(jié)合代謝路徑關(guān)鍵反應(yīng)段的信息熵和互信息來(lái)計(jì)算關(guān)鍵反應(yīng)段評(píng)分以評(píng)估關(guān)鍵反應(yīng)段的反應(yīng)搭配的豐富程度和反應(yīng)間的關(guān)聯(lián)程度,式(6)的評(píng)分越高說(shuō)明關(guān)鍵反應(yīng)段的反應(yīng)搭配越豐富,反應(yīng)間的關(guān)聯(lián)程度越高,生化相關(guān)性越好.代謝路徑的關(guān)鍵反應(yīng)段的信息熵與互信息值之間存在一定差別.為均衡評(píng)估這兩個(gè)路徑特征對(duì)關(guān)鍵反應(yīng)段評(píng)分的作用,本文分別采用式(7)和式(8)對(duì)代謝路徑的關(guān)鍵反應(yīng)段的信息熵和互信息值做歸一化處理.

基于線性路徑pa的所有關(guān)鍵反應(yīng)段ks,線性路徑pa的評(píng)分可以通過(guò)式(9)計(jì)算:

(9)

其中ad為關(guān)鍵反應(yīng)段權(quán)重,其取值范圍為[0,1];npa為線性路徑pa的長(zhǎng)度,ap為線性路徑pa的長(zhǎng)度權(quán)重.在式(9)中,線性路徑pa的所有關(guān)鍵反應(yīng)段評(píng)分之和越高、路徑長(zhǎng)度越短,說(shuō)明線性路徑pa的生化相關(guān)性越好.

基于代謝路徑PB的每條線性路徑pi的評(píng)分Score′(pi),代謝路徑PB的評(píng)分可由式(10)計(jì)算得出:

(10)

其中l(wèi)b是構(gòu)成代謝路徑PB的線性分支路徑個(gè)數(shù),pi表示構(gòu)成代謝路徑PB的一條線性路徑.最后KPRank方法通過(guò)式(10)計(jì)算代謝路徑PB的所有線性分支路徑評(píng)分的平均值得到代謝路徑PB的最終評(píng)分.

2.5 算法描述

對(duì)于待排序代謝路徑集PR,使用本文方法KPRank對(duì)待排序代謝路徑集合PR中的路徑進(jìn)行排序的基本步驟如下:

Step 1.對(duì)已知代謝路徑集合P中每一條線性路徑,枚舉其中由2至4個(gè)連續(xù)反應(yīng)構(gòu)成的反應(yīng)段,計(jì)算每個(gè)反應(yīng)段rs在已知代謝路徑集合P中的出現(xiàn)次數(shù),將出現(xiàn)次數(shù)大于100次的關(guān)鍵反應(yīng)段存入關(guān)鍵反應(yīng)段集合HRS.

Step 2.對(duì)于待排序路徑集PR中的每一條代謝路徑pl中的每條線性分支路徑pa,使用逆向最長(zhǎng)匹配算法[28]從路徑pa中提取反應(yīng)段,并保留路徑pa中屬于HRS集合的關(guān)鍵反應(yīng)段,通過(guò)式(6)計(jì)算路徑pa中的屬于HRS集合的每個(gè)關(guān)鍵反應(yīng)段評(píng)分,并通過(guò)式(9)得到線性路徑pa的評(píng)分.

Step 3.利用Step 2得到的代謝路徑pl包含的所有線性路徑的評(píng)分,采用式(10)計(jì)算出代謝路徑pl的評(píng)分.

Step 4.重復(fù)Step 2和Step 3計(jì)算待排序代謝路徑集合PR中的每條代謝路徑的評(píng)分.最后根據(jù)PR中的每條路徑的評(píng)分以從高到低的順序輸出路徑排序結(jié)果.

算法1形式描述了本文提出的方法KPRank.

算法1.KPRank

輸入:已知代謝路徑集合P,待排序路徑集合PR

輸出:待排序路徑集合PR的路徑排序結(jié)果

Begin

for allpathinPdo

RS←枚舉路徑path中由2至4個(gè)連續(xù)反應(yīng)構(gòu)成的反應(yīng)段;

end for

for allrsinRSdo

統(tǒng)計(jì)反應(yīng)段rs在P中出現(xiàn)次數(shù);

HRS←在已知代謝路徑集合P中出現(xiàn)次數(shù)>100的反應(yīng)段;

end for

for allplinPRdo /*待排序路徑集合PR中的每條代謝路徑pl*/

for allpainpldo /*代謝路徑pl中的每條線性分支路徑pa*/

采用逆向最長(zhǎng)匹配算法從路徑pa中提取屬于HRS集合的關(guān)鍵反應(yīng)段;

利用式(6)計(jì)算路徑pa的每個(gè)關(guān)鍵反應(yīng)段評(píng)分;

基于路徑pa的每個(gè)關(guān)鍵反應(yīng)段評(píng)分,通過(guò)式(9)計(jì)算得到線性路徑pa的評(píng)分;

end for

基于代謝路徑pl包含的每條線性路徑評(píng)分,采用式(10)計(jì)算出代謝路徑pl的評(píng)分;

end for

以待排序路徑集合PR中的路徑評(píng)分降序輸出路徑排序結(jié)果;

End.

3 實(shí) 驗(yàn)

3.1 實(shí)驗(yàn)環(huán)境

為測(cè)試KPRank的路徑排序效果,本文使用KEGG PATHWAY數(shù)據(jù)庫(kù)中的67個(gè)已知代謝網(wǎng)絡(luò)來(lái)構(gòu)建已知代謝路徑集合P.

BPFinder是基于原子團(tuán)追蹤的代謝路徑預(yù)測(cè)算法[17],ReTrace是基于原子追蹤的代謝路徑預(yù)測(cè)算法[16],MetQuest是基于動(dòng)態(tài)規(guī)劃的代謝路徑預(yù)測(cè)算法[7].本文分別使用3種預(yù)測(cè)代謝路徑的方法Retrace,BPFinder和MetQuest在KEGG PATHWAY數(shù)據(jù)庫(kù)中搜索26條已知代謝路徑.分別從https://www.cs.helsinki.fi/group/sysfys/software/retrace/下載了Retrace,從https://github.com/RamanLab/metquest下載了MetQuest.在實(shí)驗(yàn)中,每個(gè)代謝路徑預(yù)測(cè)方法尋找每條已知路徑的起始化合物和目標(biāo)化合物之間的代謝路徑,并返回路徑搜索結(jié)果的前5條路徑,然后利用KPRank方法分別為Retrace,BPFinder和MetQuest找到的代謝路徑排序,并分別對(duì)經(jīng)過(guò)KPRank排序后和排序前的代謝路徑與已知代謝路徑進(jìn)行比較來(lái)評(píng)估路徑排序結(jié)果的質(zhì)量.采用Java語(yǔ)言編程實(shí)現(xiàn)KPRank算法.KPRank的實(shí)驗(yàn)參數(shù)配置為:ad與ap都是0.5.方法Retrace、BPFinder和MetQuest均采用默認(rèn)運(yùn)行參數(shù).實(shí)驗(yàn)的計(jì)算機(jī)配置為Intel(R)Xeon(R)CPU 6130 @ 2.10GHz和40GB RAM.

3.2 評(píng)估指標(biāo)

通過(guò)計(jì)算工具得到的替代路徑的生化相關(guān)性會(huì)影響替代路徑在實(shí)際生物合成應(yīng)用的可用性.人們一般通過(guò)比較已知路徑與計(jì)算路徑來(lái)評(píng)估計(jì)算路徑的生化相關(guān)性.本文通過(guò)比較計(jì)算路徑與已知路徑得到的7個(gè)指標(biāo)來(lái)評(píng)估計(jì)算路徑的生化相關(guān)性[29]:

1)最大連通子圖邊數(shù)占比(RELCCS)是已知路徑與計(jì)算路徑同構(gòu)的子圖里最大的聯(lián)通子圖的邊數(shù)與已知路徑的邊數(shù)之比[29].RELCCS值越大,表明計(jì)算路徑與已知路徑越相似,也說(shuō)明計(jì)算路徑具有更好的生化相關(guān)性.RELCCS的平均值的計(jì)算如下[29]:

(11)

N為計(jì)算路徑的個(gè)數(shù),RELCCSi指第i條計(jì)算路徑與對(duì)應(yīng)這條計(jì)算路徑的已知路徑同構(gòu)的子圖里最大的連通子圖邊數(shù)與已知路徑邊數(shù)之比[29].

2)化合物靈敏度Sn=tp/(tp+fn),真陽(yáng)性tp是在計(jì)算路徑與已知路徑中同時(shí)出現(xiàn)的化合物,且這些化合物在已知路徑中的出現(xiàn)順序與它們?cè)谟?jì)算路徑中的出現(xiàn)順序一樣[12,18].假陰性fn是在計(jì)算路徑中不存在,但在已知的路徑中出現(xiàn)的化合物[12,18].

3)化合物陽(yáng)性預(yù)測(cè)值PPV=tp/(tp+fp),假陽(yáng)性fp是在計(jì)算路徑中存在但在已知路徑中不存在的化合物[12,18].

4)化合物準(zhǔn)確率AC=(SN+PPV)/2,AC值越大,說(shuō)明算法具有更好的恢復(fù)已知路徑的化合物的能力[12,18].

5)反應(yīng)靈敏度R_Sn=r_tp/(r_tp+r_fn),真陽(yáng)性r_tp是在已知路徑和計(jì)算路徑中同時(shí)出現(xiàn)的反應(yīng),并且這些反應(yīng)在已知路徑中的出現(xiàn)順序與在計(jì)算路徑中的出現(xiàn)順序一樣.假陰性r_fn是在計(jì)算路徑中不存在,但在已知路徑中存在的反應(yīng)[12,18].

6)反應(yīng)陽(yáng)性預(yù)測(cè)值R_PPV=r_tp/(r_tp+r_fp),假陽(yáng)性r_fp是在計(jì)算路徑中存在,但在已知路徑中不存在的反應(yīng)[12,18].

7)反應(yīng)準(zhǔn)確率R_AC=(R_SN+R_PPV)/2,R_AC值越大,說(shuō)明算法具有更好的恢復(fù)已知路徑反應(yīng)的能力[12,18].

3.3 計(jì)算路徑和已知路徑結(jié)構(gòu)的比較

表1顯示了各種方法尋找26條已知路徑所找到的代謝路徑的平均RELCCS值.其中KPRank(BPFinder)、KPRank(Retrace)和KPRank(MetQuest)分別表示用KPRank方法對(duì)BPFinder、Retrace、MetQuest方法所找到的代謝路徑進(jìn)行排序以后得到的結(jié)果.

表1 26條計(jì)算路徑的平均RELCCS值和前5條結(jié)果路徑中分支路徑的數(shù)量

由表1可知,對(duì)于方法BPFinder、Retrace、MetQuest所找到的前5條路徑里性能最好的路徑,使用方法KPRank為結(jié)果路徑排序后得到的路徑的RELCCS值皆優(yōu)于未使用方法KPRank排序得到的路徑的RELCCS值.這些結(jié)果表明經(jīng)過(guò)方法KPRank對(duì)路徑結(jié)果進(jìn)一步排序后,所得到的路徑結(jié)構(gòu)更接近已知路徑,路徑搜索質(zhì)量得到進(jìn)一步提高.另一方面,從表1可以看到,經(jīng)過(guò)執(zhí)行方法KPRank排序后,Retrace的首條路徑和前5條路徑里性能最好的路徑的RELCCS值均比原來(lái)的方法更好,而B(niǎo)PFinder和MetQuest的首條路徑的RELCCS值有一定程度的減少.結(jié)合找到的分支路徑數(shù)量來(lái)看,相比其他方法,Retrace找到的前5條路徑的分支路徑數(shù)量更多.這說(shuō)明方法KPRank更適用于協(xié)助分析和尋找分支多、反應(yīng)組合復(fù)雜的分支替代路徑.與此同時(shí),從表1還可以看到,3種方法的首條路徑的RELCCS值都不是前5條路徑里最好的,這表明要尋找拓?fù)浣Y(jié)構(gòu)上更接近已知路徑的替代路徑,首條路徑不一定是最好的選擇.

3.4 計(jì)算路徑和已知路徑中的化合物和反應(yīng)的比較

本文采用AC,SN和PPV來(lái)度量計(jì)算路徑的化合物精度、敏感度和陽(yáng)性預(yù)測(cè)值,采用R_AC,R_SN和R_PPV來(lái)度量計(jì)算路徑的反應(yīng)精度、敏感度和陽(yáng)性預(yù)測(cè)值.表2給出了各種方法的26條計(jì)算路徑的反應(yīng)的平均精度AC、平均敏感度SN和平均陽(yáng)性預(yù)測(cè)值PPV.

表2 26條計(jì)算路徑的反應(yīng)平均R_SN值、平均R_PPV值和平均R_AC值

由表2可以看到,對(duì)于BPFinder和Retrace的前5條路徑里性能最好的路徑,使用方法KPRank為結(jié)果路徑排序后得到的路徑的R_Sn值,R_PPV值和R_AC值都比未使用方法KPRank排序得到的路徑的R_Sn值,R_PPV值和R_AC值更好.對(duì)于BPFinder的首條路徑,經(jīng)過(guò)方法KPRank為結(jié)果路徑排序后得到的路徑的R_PPV值和R_AC值也比沒(méi)有使用方法KPRank排序得到的路徑的R_PPV值和R_AC值更好.對(duì)于MetQuest的首條路徑和前5條路徑里性能最好的路徑,沒(méi)有采用方法KPRank排序得到的路徑的R_PPV值和R_AC值相對(duì)更好一些.這些結(jié)果說(shuō)明KPRank更適用于BPFinder和Retrace,KPRank對(duì)BPFinder和Retrace的路徑結(jié)果進(jìn)一步排序可以搜索出路徑的反應(yīng)組成更接近已知路徑的替代路徑.

表3給出了各種方法的26條計(jì)算路徑的化合物的平均精度AC、平均敏感度SN和平均陽(yáng)性預(yù)測(cè)值PPV.

表3 26條計(jì)算路徑的化合物平均SN值、平均PPV值和平均AC值

從表3可以看到,經(jīng)過(guò)方法KPRank排序后,雖然BPFinder和Retrace的首條路徑的Sn值和AC值相對(duì)有所下降,MetQuest的首條路徑的Sn值和AC值比未使用方法KPRank排序得到的路徑的Sn值和AC值好一些.相對(duì)于首條路徑,前5條路徑里性能最好的路徑更有代表性,對(duì)于方法BPFinder、Retrace、MetQuest所找到的前5條路徑里性能最好的路徑,采用方法KPRank為結(jié)果路徑排序后得到的路徑的Sn值、PPV值和AC值都好于未使用方法KPRank排序得到的路徑的Sn值,PPV值和AC值.這些結(jié)果表明經(jīng)過(guò)方法KPRank對(duì)這些路徑預(yù)測(cè)方法的路徑搜索結(jié)果進(jìn)一步排序后,所得到的路徑中的化合物更接近已知路徑中的化合物,搜索得到的路徑具有更好的還原已知代謝路徑的化合物的能力.

表1～表3的結(jié)果說(shuō)明KPRank排序方法能夠提升代謝路徑預(yù)測(cè)方法搜索路徑的生化相關(guān)性,可以有效提高代謝路徑預(yù)測(cè)結(jié)果的質(zhì)量.

3.5 實(shí)例分析:磷酸戊糖合成路徑的搜索

上述實(shí)驗(yàn)結(jié)果表明基于路徑關(guān)鍵反應(yīng)段的路徑排序方法KPRank有效的提高了代謝路徑預(yù)測(cè)方法的路徑搜索結(jié)果的生化相關(guān)性.本文通過(guò)討論分析一個(gè)利用方法KPRank對(duì)BPFinder的路徑搜索結(jié)果經(jīng)過(guò)排序以后得到的路徑實(shí)例來(lái)進(jìn)一步了解方法KPRank的特點(diǎn).

5-磷酸核糖(D-Ribose 5-phosphate)是核苷酸和核酸合成所必需的底物,它可以通過(guò)磷酸戊糖合成路徑中非氧化階段由β-D-果糖6-磷酸(beta-D-Fructose 6-phosphate)為起始化合物通過(guò)代謝合成[30].磷酸戊糖路徑存在于細(xì)胞的胞漿中,它同時(shí)是葡萄糖氧化的替代途徑[31],并可用于生產(chǎn)NADPH和5-磷酸核糖[30].

圖2(a)是KEGG數(shù)據(jù)庫(kù)的已知磷酸戊糖合成路徑,圖2(b)是使用方法KPRank對(duì)BPFinder的磷酸戊糖合成路徑的路徑搜索結(jié)果排序后得到的首條路徑,圖2(c)是BPFinder搜索磷酸戊糖合成路徑得到的首條路徑.在圖2中,虛線橢圓中的化合物表示代謝路徑中的分支化合物,計(jì)算路徑中與已知路徑相同的反應(yīng)和化合物均以斜體加粗表示.

圖2 磷酸戊糖合成路徑

從圖2(a)可以看到反應(yīng)R01529和R01056構(gòu)成的反應(yīng)段是已知磷酸戊糖合成路徑的一部分.基于式(6)的計(jì)算,由反應(yīng)R01529和R01056構(gòu)成的反應(yīng)段的關(guān)鍵反應(yīng)段評(píng)分是0.99,這說(shuō)明這個(gè)反應(yīng)段在從67個(gè)已知代謝網(wǎng)絡(luò)中提取的已知代謝路徑集合里出現(xiàn)頻率較高,反應(yīng)結(jié)構(gòu)功能特征明顯,因此在圖2(b)中,KPRank將包含這個(gè)反應(yīng)段的路徑作為首條路徑輸出.有趣的是,圖2(b)的這條經(jīng)過(guò)KPRank排序后得到的首條磷酸戊糖合成路徑位于未經(jīng)KPRank排序的BPFinder發(fā)現(xiàn)的磷酸戊糖合成路徑的第13位,這表明基于關(guān)鍵反應(yīng)段的評(píng)分策略可以更快速有效的發(fā)現(xiàn)包含功能反應(yīng)段的路徑.

對(duì)比圖2(a)和圖2(c)可以看到,雖然BPFinder搜索磷酸戊糖合成路徑得到的首條路徑中有部分化合物和反應(yīng)出現(xiàn)在已知路徑里,但是BPFinder的這條從起始化合物beta-D-Fructose 6-phosphate到目標(biāo)化合物D-Ribose 5-phosphate的路徑分支眾多、跨度大,有很大部分與已知路徑不同,這在一定程度上可能會(huì)增加設(shè)計(jì)實(shí)現(xiàn)新的磷酸戊糖合成路徑的難度.對(duì)比圖2(a)和圖2(b)還可以看到,KPRank發(fā)現(xiàn)了已知路徑中沒(méi)有的從起始化合物beta-D-Fructose 6-phosphate開(kāi)始,經(jīng)過(guò)D-Glyceraldehyde 3-phosphate最終到達(dá)目標(biāo)化合物D-Ribose 5-phosphate的新替代路徑,同時(shí)這條新替代路徑中從D-Glyceraldehyde 3-phosphate開(kāi)始到目標(biāo)化合物D-Ribose 5-phosphate之間的大部分路徑與已知路徑相同,這將在一定程度上降低設(shè)計(jì)實(shí)現(xiàn)新的磷酸戊糖合成路徑的難度.

4 結(jié)束語(yǔ)

本文提出了的代謝路徑排序方法的特色是:首先從67個(gè)已知的代謝網(wǎng)絡(luò)提取高頻關(guān)鍵反應(yīng)段,并融合代謝路徑的關(guān)鍵反應(yīng)段的信息熵與互信息對(duì)代謝路徑預(yù)測(cè)算法找到的路徑進(jìn)行評(píng)分和排序,以獲取生化相關(guān)性好的替代路徑.實(shí)驗(yàn)結(jié)果表明,利用KPRank對(duì)代謝路徑預(yù)測(cè)算法的結(jié)果路徑進(jìn)行排序可以找出生化相關(guān)性更好的路徑,進(jìn)一步提高了預(yù)測(cè)的結(jié)果路徑的質(zhì)量.

下一步工作將考慮利用代謝數(shù)據(jù)庫(kù)的反應(yīng)規(guī)則對(duì)已有代謝網(wǎng)絡(luò)中的代謝反應(yīng)進(jìn)行功能性分類,并結(jié)合關(guān)鍵反應(yīng)段和反應(yīng)功能分類對(duì)路徑預(yù)測(cè)結(jié)果進(jìn)行排序,以為代謝路徑的分析和合成提供更豐富的替代路徑選擇.