黃毅然，萬志遠，鐘 誠

(1.廣西大學(xué) 計算機與電子信息學(xué)院，廣西 南寧 530004；2.廣西大學(xué) 廣西高校并行與分布式計算技術(shù)重點實驗室，廣西 南寧 530004；3.廣西大學(xué) 廣西多媒體通信與網(wǎng)絡(luò)技術(shù)重點實驗室，廣西 南寧 530004)

0 引 言

代謝路徑是生物體通過一組代謝反應(yīng)將起始代謝物轉(zhuǎn)化為目標(biāo)代謝物的路徑。代謝路徑預(yù)測是生物合成路徑的設(shè)計-構(gòu)建-測試循環(huán)中的重要步驟之一，主要目標(biāo)是通過計算方法在代謝數(shù)據(jù)庫中尋找到新的有價值的合成路徑[1-3]。

在代謝數(shù)據(jù)庫中找到的合成代謝路徑通常包含連接度很高的簇代謝物[4，5]。簇代謝物的出現(xiàn)可能會產(chǎn)生生物學(xué)上沒有意義的代謝路徑。原子交換是代謝路徑中的代謝反應(yīng)中的底物和產(chǎn)物之間存在原子轉(zhuǎn)移的現(xiàn)象，研究結(jié)果表明通過利用原子交換信息搜索路徑可以有效排除合成路徑中的簇代謝物[6-8]。

一些基于約束的代謝路徑預(yù)測方法利用代謝網(wǎng)絡(luò)中代謝物和代謝反應(yīng)具有的生化特性，構(gòu)建化學(xué)計量MILP模型以尋找代謝路徑，并通過在MILP模型中利用原子交換信息來尋找沒有簇代謝物的替代路徑以提高路徑搜索質(zhì)量[9]。NetFlow方法通過檢測到的碳通量搜索具有生化相關(guān)性的代謝路徑[10]。NFP方法[11，12]結(jié)合代謝網(wǎng)絡(luò)中的拓撲結(jié)構(gòu)和化學(xué)計量信息以構(gòu)建特定的MILP模型，通過MILP模型實現(xiàn)從任意起始代謝物到目標(biāo)代謝物的代謝路徑搜索。雖然NFP方法可以搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑，但它沒有將代謝反應(yīng)中的原子轉(zhuǎn)移信息融入MILP模型中以尋找生化可行的代謝路徑，使得該方法獲得的代謝路徑的生化相關(guān)性較弱。

本文研究提出一種基于約束的代謝路徑預(yù)測方法PVA，以搜索出從任意起始代謝物到目標(biāo)代謝物的反應(yīng)中包含特定原子交換的代謝路徑。不同于搜索給定的起始代謝物和目標(biāo)代謝物之間代謝路徑的路徑預(yù)測方法[6-8]，本文提出的方法可以搜索從任意起始代謝物到給定目標(biāo)代謝物的代謝路徑。此外，與NFP方法不同，本文提出的方法將代謝反應(yīng)中特定原子交換信息與反應(yīng)化學(xué)計量相結(jié)合，建立基于約束的MILP模型以找到含有特定原子交換的代謝路徑，以進一步提高從任意起始代謝物開始搜索路徑所找到的路徑的生化相關(guān)性，更好地發(fā)現(xiàn)與已知代謝路徑相近的替代路徑。

1 方 法

本文方法的主要思想是：首先，從KEGG數(shù)據(jù)庫中提取在輸入代謝物和輸出代謝物之間存在特定原子交換的代謝反應(yīng)。然后，將反應(yīng)的化學(xué)計量信息與包含特定原子交換信息的代謝反應(yīng)結(jié)合起來，構(gòu)建基于約束的MILP模型。最后，通過求解MILP模型來找到從任意起始代謝物到目標(biāo)代謝物的反應(yīng)中包含特定原子交換的代謝路徑。

1.1 代謝路徑預(yù)測模型

首先利用KEGG RPAIR數(shù)據(jù)庫中輸入代謝物和輸出代謝物之間存在特定原子交換的代謝反應(yīng)來重構(gòu)代謝網(wǎng)絡(luò)。在重構(gòu)的代謝網(wǎng)絡(luò)中，頂點代表代謝物，弧表示代謝物之間的代謝反應(yīng)，弧中的反應(yīng)包含輸入代謝物和輸出代謝物之間的原子交換特征信息。

重構(gòu)的代謝網(wǎng)絡(luò)由一組反應(yīng)R和一組代謝物M組成。本文利用化學(xué)計量矩陣S、二元系數(shù)G、一組內(nèi)部代謝物I、一組外部代謝物E、一組起始代謝物SM和一組可逆反應(yīng)RB來刻畫重構(gòu)的代謝網(wǎng)絡(luò)。在化學(xué)計量矩陣S中，行代表代謝物，列代表代謝反應(yīng)。Smr表示反應(yīng)r中代謝物m的化學(xué)計量值。Smr<0表示m作為反應(yīng)r的輸入代謝物，Smr>0表示m作為反應(yīng)r的輸出代謝物。內(nèi)部代謝物集合I給出了代謝網(wǎng)絡(luò)中的內(nèi)部代謝物列表。外部代謝物集合E給出了需要通過預(yù)測路徑中的代謝反應(yīng)產(chǎn)生的外部代謝物列表。起始代謝物集合SM給出了可作為預(yù)測路徑中的第一個代謝物的代謝物列表?？赡娣磻?yīng)集RB給出了重構(gòu)的代謝網(wǎng)絡(luò)中的可逆反應(yīng)。

依據(jù)上述模型信息，本文使用如下代謝路徑預(yù)測的約束方程[11]

(1)

(2)

(3)

(4)

(5)

在上述約束方程中，二元變量aij表示輸入代謝物i和輸出代謝物j之間的弧，其中|A|表示代謝網(wǎng)絡(luò)中有|A|個不同的代謝物。如果輸入代謝物i到輸出代謝物j的弧aij存在于預(yù)測路徑中，則aij=1，否則為0，i，j=1，2，…，|A|。 約束方程(1)選擇目標(biāo)代謝物p為預(yù)測路徑的終點[11]。約束方程(2)確保目標(biāo)代謝物p不會是預(yù)測路徑的起點[11]。約束方程(3)確保起始代謝物集合SM中的代謝物l或者是中間代謝物，或者是起始代謝物[11]。約束方程(4)確保起始代謝物集合以外的其它代謝物都可以作為代謝路徑的中間代謝物[11]。約束方程(5)確保預(yù)測路徑中的代謝物都是唯一的[11]。約束方程(1)～方程(5)確保代謝路徑預(yù)測方法可以從任意起始代謝物開始搜索，其中起始代謝物可以通過約束方程(3)從起始代謝物集合SM中隨機選擇。

通過約束方程(1)～方程(5)只是獲得代謝路徑中簡單的代謝物序列。下面將通過約束方程(6)～方程(10)定義代謝路徑中的反應(yīng)通量

(6)

(7)

cr≤vr≤Max*crr=1，…，|N|

(8)

ca+cb≤1 ?(a，b)∈RB

(9)

(10)

在約束方程(6)～方程(10)中，vr是反應(yīng)r的通量，它為每個反應(yīng)r分配一個非負通量，|N|表示代謝網(wǎng)絡(luò)中有|N|個的不同代謝反應(yīng)，r=1，2，…，|N|。 二元變量cr用作控制反應(yīng)通量vr的開閉，cr=0使得相關(guān)通量vr為0，而cr=1使vr為非零值，r=1，2，…，|N|。Max大于或等于1。

約束方程(6)～方程(10)定義了預(yù)測路徑的有效通量分布以確保找到的預(yù)測路徑在生物學(xué)上是可行的，即該預(yù)測路徑可以經(jīng)過多個代謝反應(yīng)由起始代謝物合成目標(biāo)代謝物。約束方程(6)確保外部代謝物只能產(chǎn)生，而不能消耗[11]。約束方程(7)確保目標(biāo)代謝物p作為產(chǎn)物存在于預(yù)測路徑中[11]。約束方程(8)通過二元變量cr控制反應(yīng)通量vr的切換[11]。約束方程(9)約束代謝網(wǎng)絡(luò)中的可逆反應(yīng)不能同時出現(xiàn)在預(yù)測路徑中[11]。

本文結(jié)合約束方程(1)～方程(10)來求解MILP模型中的目標(biāo)函數(shù)式(11)，以獲得從任意起始代謝物到目標(biāo)代謝物p的反應(yīng)中含有特定原子交換的預(yù)測路徑

(11)

1.2 算 法

對于給定目標(biāo)代謝物T以及被追蹤的原子AM，使用融合原子交換特征信息的代謝路徑預(yù)測算法PVA，搜索從任意起始代謝物到目標(biāo)代謝物T的反應(yīng)中含有原子AM交換的代謝路徑的步驟如下：

步驟1 從KEGG RPAIR數(shù)據(jù)庫中提取輸入代謝物和輸出代謝物之間存在原子AM交換的反應(yīng)，利用存在原子AM交換的代謝反應(yīng)以建立MILP模型。

步驟2 對于目標(biāo)代謝物T，通過MILP模型中的約束方程(1)～方程(5)找到從任意起始代謝物到目標(biāo)代謝物T的代謝路徑的代謝物序列。通過MILP模型中的約束方程(6)～方程(9)獲取從任意起始代謝物到目標(biāo)代謝物T的代謝路徑的代謝反應(yīng)通量。通過MILP模型中的約束方程(10)從代謝路徑中篩選出存在原子AM交換的代謝反應(yīng)。

步驟3 結(jié)合MILP模型約束方程(1)～方程(10)求解方程(11)的MILP模型中的目標(biāo)函數(shù)，以獲得從任意起始代謝物到目標(biāo)代謝物T的反應(yīng)中含原子AM交換的代謝路徑。

算法1形式描述了本文提出的算法PVA。

算法1：PVA

輸入：目標(biāo)代謝物T，追蹤的目標(biāo)原子AM，KEGG RPAIR數(shù)據(jù)庫的代謝反應(yīng)集合reactions

輸出：任意起始代謝物到目標(biāo)代謝物T的反應(yīng)中含有原子AM交換的代謝路徑

Begin

(1) for allrin reactions do

RAM← 遍歷reactions中每個反應(yīng)，獲取代謝反應(yīng)r中具有原子AM交換特征信息的反應(yīng)；

end for

(2)POM←利用MILP模型中的約束方程(1)～(5)找到任意起始代謝物到目標(biāo)代謝物T之間的代謝物序列；

(3)RS← 利用MILP模型中的約束方程(6)～(9)找到與代謝物序列POM相匹配的代謝反應(yīng)通量；

(4) for allrinRSdo

POR←通過MILP模型中的約束方程(10)從RS中篩選出代謝反應(yīng)r中存在原子AM交換的反應(yīng)；

end for

(5)path←結(jié)合任意起始代謝物到目標(biāo)代謝物T之間的代謝物序列POM和存在原子AM交換的反應(yīng)集合POR，求解MILP模型中目標(biāo)函數(shù)(方程(11))以獲得從任意起始代謝物到目標(biāo)代謝物T的反應(yīng)中含有原子AM交換的代謝路徑；

(6) 輸出預(yù)測路徑結(jié)果path；

End.

2 實 驗

2.1 實驗環(huán)境

通過比較預(yù)測路徑與來自KEGG的大腸桿菌數(shù)據(jù)庫中的20條已知路徑來評估代謝路徑預(yù)測算法的有效性。

在實驗中，本文方法PVA采用4種代謝路徑預(yù)測模式來尋找從任意起始代謝物到目標(biāo)代謝物的反應(yīng)中含有特定原子交換的代謝路徑。PVA-C表示通過追蹤反應(yīng)中的碳原子交換來搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑。PVA-O表示通過追蹤反應(yīng)中的氧原子交換來搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑。PVA-N表示通過追蹤反應(yīng)中的氮原子交換來搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑。PVA-P表示通過追蹤反應(yīng)中的磷原子交換來搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑。

通過MILP模型搜索代謝路徑的預(yù)測方法有兩種。第一種方法(表示為“Topology”)在KEGG的原代謝網(wǎng)絡(luò)中使用MILP模型搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑。第二種方法(表示為“Hubs”)在調(diào)整后的代謝網(wǎng)絡(luò)中使用MILP模型搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑，調(diào)整后的代謝網(wǎng)絡(luò)是指刪除了高度連接的代謝物(Hubs)的代謝網(wǎng)絡(luò)[13]。

本文分別利用方法PVA-C、PVA-O、PVA-N、PVA-P、NFP[11]、Topology和Hubs[13]在KEGG數(shù)據(jù)庫中搜索20條已知代謝路徑。每種方法搜索從任意起始代謝物到20條已知代謝路徑中的目標(biāo)代謝物的代謝路徑，并通過將計算得到的替代路徑與已知代謝路徑進行比較來評估代謝路徑預(yù)測方法的有效性。

本文采用Java和GUROBI編程實現(xiàn)PVA算法。基于約束的代謝路徑預(yù)測方法NFP和其它兩種代謝路徑預(yù)測方法(“Topology”和“Hubs”)的MILP模型都使用GUROBI在Java環(huán)境中求解。所有實驗均在 Intel Xeon 6130和192 GB RAM的計算機上運行。

2.2 性能指標(biāo)

代謝路徑預(yù)測方法得到的預(yù)測路徑的生化相關(guān)性會影響預(yù)測路徑在代謝工程中生物合成應(yīng)用的可行性。本文通過將預(yù)測得到的替代路徑與已知路徑進行比較以評估預(yù)測路徑的生化相關(guān)性。與同類研究方法相同，本文以代謝反應(yīng)準(zhǔn)確率[7]和代謝物準(zhǔn)確率[7]作為評估預(yù)測路徑的性能指標(biāo)，具體計算方法如下：

(1)代謝反應(yīng)準(zhǔn)確率R_ACC=(R_Sn+R_PPV)/2， 其中R_S是代謝反應(yīng)的靈敏度，R_Sn=R_TP/(R_TP+R_FP)，R_PPV是代謝反應(yīng)的陽性預(yù)測值，R_PPV=R_TP/(R_TP+R_FN)[7]。 其中真陽性R_TP表示預(yù)測路徑中存在的代謝反應(yīng)同已知路徑中的代謝反應(yīng)相同，且在預(yù)測路徑和已知路徑中的代謝反應(yīng)出現(xiàn)的相對順序是相同的；假陽性R_FP表示存在于預(yù)測路徑中，但不存在于已知路徑中的反應(yīng)；假陰性R_FN表示不存在于預(yù)測路徑中，而是存在于已知路徑中的反應(yīng)[7]。代謝反應(yīng)準(zhǔn)確率R_ACC值越高說明代謝路徑預(yù)測算法返回的預(yù)測路徑中的代謝反應(yīng)同已知路徑中的代謝反應(yīng)越趨于一致[7]。

(2)代謝物準(zhǔn)確率C_ACC=(C_Sn+C_PPV)/2， 其中C_Sn是代謝物靈敏度，C_Sn=TP/(TP+FP)，C_PPV是代謝物的陽性預(yù)測值，C_PPV=TP/(TP+FN)[7]。 真陽性TP表示預(yù)測路徑中存在的代謝物同已知路徑中的代謝物相同，且在預(yù)測路徑和已知路徑中的代謝物出現(xiàn)的相對順序是相同的；假陽性FP表示存在于預(yù)測路徑中，但不存在于已知路徑中的代謝物；假陰性FN表示代謝物不存在于預(yù)測路徑中，而是存在于已知路徑中[7]。代謝物準(zhǔn)確率C_ACC值越高說明代謝路徑預(yù)測方法能夠更準(zhǔn)確地發(fā)現(xiàn)已知路徑中存在的代謝物[7]。

2.3 預(yù)測路徑與已知路徑中代謝反應(yīng)和代謝物的比較

對于20條已知路徑中的每個目標(biāo)代謝物，分別采用預(yù)測方法PVA-C、PVA-O、PVA-N、PVA-P、NFP、Topo-logy和Hubs搜索從任意起始代謝物到目標(biāo)代謝物的代謝路徑。通過將每種方法找到的替代路徑與已知路徑進行比較，以評估代謝路徑預(yù)測方法的性能。

表1給出了各算法預(yù)測得到的替代路徑中的首條路徑、前5條路徑、前10條路徑以及前100路徑的平均代謝反應(yīng)準(zhǔn)確率R_ACC值。

表1 預(yù)測路徑的平均代謝反應(yīng)準(zhǔn)確率R_ACC值

從表1可以看出，對于首條路徑、前5條路徑、前10條路徑以及前100路徑，PVA-O獲得的預(yù)測路徑的R_ACC值都比NFP、Topology、Hubs、PVA-C、PVA-N和PVA-P 獲得的預(yù)測路徑的R_ACC值更高；此外PVA-C和NFP獲得的預(yù)測路徑的R_ACC值與PVA-O獲得的預(yù)測路徑的R_ACC值相近。這些結(jié)果表明，相比其它基于約束的代謝路徑預(yù)測方法，PVA-O能夠更準(zhǔn)確地恢復(fù)已知路徑中的代謝反應(yīng)。

從表1中還可以看到，對于首條路徑、前5條路徑、前10條路徑以及前100路徑，PVA-P獲得的預(yù)測路徑的R_ACC值低于PVA-C、PVA-O、PVA-N和NFP獲得的預(yù)測路徑的R_ACC值；PVA-P獲得的預(yù)測路徑的R_ACC值與Topology和Hubs獲得的預(yù)測路徑的R_ACC值相當(dāng)。這些結(jié)果表明，在基于約束的MILP模型中追蹤不同的原子交換對發(fā)現(xiàn)已知路徑中的代謝反應(yīng)有較大的影響。

表2給出了各算法預(yù)測得到的替代路徑中的首條路徑、前5條路徑、前10條路徑以和前100路徑的平均代謝物準(zhǔn)確率C_ACC值。

表2 預(yù)測路徑的平均代謝物準(zhǔn)確率C_ACC值

從表2中可以觀察到，對于首條路徑、前5條路徑、前10條路徑以及前100路徑，PVA-O獲得的預(yù)測路徑的C_ACC值也都高于NFP、Topology、Hubs、PVA-C、PVA-N和PVA-P獲得的預(yù)測路徑的C_ACC值。這些結(jié)果表明PVA-O能夠更準(zhǔn)確地找到已知路徑中的代謝物。另一方面，從表2還可以看到，對于首條路徑、前5條路徑、前10條路徑以及前100路徑，與PVA的其它追蹤模式相比，PVA-N和PVA-P獲得的預(yù)測路徑的C_ACC值比較接近，但是都低于PVA-C、PVA-O獲得的預(yù)測路徑的C_ACC值。這些結(jié)果表明，在基于約束的MILP模型中通過追蹤磷原子和氮原子可能較難獲得生化相關(guān)較好的代謝路徑。

表1和表2的結(jié)果表明，融合原子交換特征信息的代謝路徑預(yù)測算法PVA能夠有效地找到與已知路徑更接近的預(yù)測路徑，其中PVA-O的預(yù)測路徑的反應(yīng)準(zhǔn)確率和代謝物準(zhǔn)確率都較高，說明已知路徑反應(yīng)中的氧原子交換特征更顯著，人們可以通過追蹤氧原子在代謝路徑反應(yīng)中的交換更容易地發(fā)現(xiàn)生化相關(guān)性較好的替代路徑。

2.4 實例分析：UDP-葡萄糖代謝路徑的搜索

第2.3節(jié)的實驗結(jié)果表明融合原子交換特征信息的代謝路徑預(yù)測算法PVA可以有效發(fā)現(xiàn)生化相關(guān)性較好的代謝路徑。下面將通過一個實例來研究PVA所找到的路徑的生化特點。

UDP-葡萄糖(UDP-glucose)是一種核苷酸糖，在許多糖基化反應(yīng)中充當(dāng)葡萄糖殘基的供體，它可以通過以α-D-葡萄糖(alpha-D-Glucose)為起始代謝物經(jīng)過多個中間代謝物由多個代謝反應(yīng)合成獲得[14，15]。同時，UDP-葡萄糖在一系列代謝物的糖基化中必不可少[15]。

圖1(a)是KEGG數(shù)據(jù)庫中從起始代謝物alpha-D-Glucose到目標(biāo)代謝產(chǎn)物UDP-glucose的已知代謝路徑，該已知代謝路徑存在于大腸桿菌和酵母菌中。圖1(b)是PVA-C搜索從任意起始代謝物到目標(biāo)代謝產(chǎn)物UDP-glucose所獲得的前5條預(yù)測路徑，圖1(c)是PVA-O搜索從任意起始代謝物到目標(biāo)代謝產(chǎn)物UDP-glucose所獲得的前5條預(yù)測路徑，圖1(d)是PVA-P搜索從任意起始代謝物到目標(biāo)代謝產(chǎn)物UDP-glucose所獲得的前5條預(yù)測路徑。在圖1中，預(yù)測路徑中與已知路徑相同的代謝反應(yīng)和代謝物分別用虛線矩形和虛線橢圓繪制。

圖1 UDP-glucose合成路徑

從圖1中可以發(fā)現(xiàn)PVA-C、PVA-O和PVA-P選取了6種不同的代謝物作為起始代謝物搜索代謝路徑，其中就包含起始代謝物alpha-D-Glucose。

從圖1(b)和圖1(c)可以看到，PVA-C和PVA-O分別發(fā)現(xiàn)了從起始代謝物alpha-D-Glucose到目標(biāo)代謝物UDP-glucose的預(yù)測路徑b4和d4，這兩條代謝路徑都是大腸桿菌和酵母菌中常用于合成目標(biāo)代謝物UDP-glucose的代謝路徑[14]。對比圖1(b)和圖1(c)還可看出，PVA-C和追蹤氧原子交換的方法獲得的前5條預(yù)測路徑是相同，只是返回的路徑排序不同。這些結(jié)果表明通過追蹤碳原子和氧原子在代謝路徑反應(yīng)中的交換能較好地尋找到與合成代謝物UDP-glucose的已知路徑相似的替代路徑，在一定程度上可以幫助代謝工程研究人員更方便地設(shè)計合成目標(biāo)代謝物UDP-glucose的替代路徑。

從圖1(b)、圖1(c)和圖1(d)還可以看到，雖然PVA-P沒有返回已知路徑，但是PVA-P找到的從起始代謝物alpha-D-Glucose6-phosphate到目標(biāo)代謝物UDP-glucose的路徑f1有很大一部分與已知路徑相同。這些結(jié)果表明，通過追蹤磷原子在代謝路徑反應(yīng)中的交換也可以為研究人員發(fā)現(xiàn)合成目標(biāo)代謝物UDP-glucose提供潛在有用的替代路徑。

3 結(jié)束語

本文提出的代謝路徑預(yù)測算法的特色是：從KEGG代謝數(shù)據(jù)庫中獲取含有原子交換信息的代謝反應(yīng)和代謝物，并結(jié)合反應(yīng)化學(xué)計量信息構(gòu)建基于約束的MILP模型，通過求解該MILP模型以獲取不限定起始代謝物的生化相關(guān)性較好的代謝路徑。實驗結(jié)果表明，通過在基于約束的代謝路徑預(yù)測模型中追蹤原子交換，能夠獲得從任意起始代謝物到目標(biāo)代謝物的反應(yīng)中含有原子交換的生化相關(guān)性更好的替代路徑，為合成路徑的設(shè)計提供啟發(fā)和參考。