黃毅然，萬志遠，鐘 誠

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

0 引 言

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

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

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

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

1 方 法

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

1.1 代謝路徑預測模型

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

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

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

(1)

(2)

(3)

(4)

(5)

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

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

(6)

(7)

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

(8)

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

(9)

(10)

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

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

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

(11)

1.2 算 法

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

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

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

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

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

算法1：PVA

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

輸出：任意起始代謝物到目標代謝物T的反應中含有原子AM交換的代謝路徑

Begin

(1) for allrin reactions do

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

end for

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

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

(4) for allrinRSdo

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

end for

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

(6) 輸出預測路徑結果path；

End.

2 實 驗

2.1 實驗環(huán)境

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

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

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

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

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

2.2 性能指標

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

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

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

2.3 預測路徑與已知路徑中代謝反應和代謝物的比較

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

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

表1 預測路徑的平均代謝反應準確率R_ACC值

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

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

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

表2 預測路徑的平均代謝物準確率C_ACC值

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

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

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

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

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

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

圖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到目標代謝物UDP-glucose的預測路徑b4和d4，這兩條代謝路徑都是大腸桿菌和酵母菌中常用于合成目標代謝物UDP-glucose的代謝路徑[14]。對比圖1(b)和圖1(c)還可看出，PVA-C和追蹤氧原子交換的方法獲得的前5條預測路徑是相同，只是返回的路徑排序不同。這些結果表明通過追蹤碳原子和氧原子在代謝路徑反應中的交換能較好地尋找到與合成代謝物UDP-glucose的已知路徑相似的替代路徑，在一定程度上可以幫助代謝工程研究人員更方便地設計合成目標代謝物UDP-glucose的替代路徑。

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

3 結束語

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