黃祖成 沈夢圓 侯至丞* TOKUYASU Taku Andrew 蒙海林*

1(廣州中國科學(xué)院先進(jìn)技術(shù)研究所 機(jī)器人與智能裝備研究中心 廣州 511458)

2(廣州中國科學(xué)院先進(jìn)技術(shù)研究所 生物工程研究中心 廣州 511458)

3(中國科學(xué)院深圳先進(jìn)技術(shù)研究院 深圳合成生物學(xué)創(chuàng)新研究院 中國科學(xué)院定量工程生物學(xué)重點(diǎn)實(shí)驗(yàn)室廣東省合成基因組學(xué)重點(diǎn)實(shí)驗(yàn)室 深圳 518055)

1 引 言

代謝通路(Metabolic Pathways)是指生物體內(nèi)將源代謝物轉(zhuǎn)化為目標(biāo)代謝物的一系列連續(xù)的酶催化反應(yīng)[1]，如糖酵解通路、三羧酸循環(huán)、梭狀芽胞桿菌固定二氧化碳的 Wood-Ljungdahl 通路[2]、藍(lán)藻吸收氨的谷氨酰胺合成酶循環(huán)(GSGOGAT)[3]等。在合成生物學(xué)領(lǐng)域，常常需要對代謝通路進(jìn)行設(shè)計(jì)與分析。代謝通路設(shè)計(jì)的目標(biāo)是給定一個(gè)或者若干個(gè)起始和目標(biāo)化合物，從代謝數(shù)據(jù)中找到并返回生物相關(guān)的、由酶催化反應(yīng)構(gòu)成的、有特定功能和作用的、從起始化合物到目標(biāo)化合物的可行路徑[4]。在過去的十多年里，隨著實(shí)驗(yàn)技術(shù)的進(jìn)步，高通量組學(xué)數(shù)據(jù)日益增長，越來越多的代謝通路得以解析。以 KEGG[5]和 MetaCyc[6]為代表的代謝數(shù)據(jù)庫發(fā)展迅速，為識別和發(fā)現(xiàn)新的合成替代路徑帶來了新的機(jī)遇和挑戰(zhàn)。研究人員根據(jù)不同物種的代謝數(shù)據(jù)特點(diǎn)和生物合成的實(shí)際應(yīng)用需求，開發(fā)了多個(gè)代謝通路設(shè)計(jì)算法和工具，如 MAPPS[7]、路徑搜索系統(tǒng)[8]、novoPathFinder[9]、MRSD[4]、FogLight[10]和 Metabolic Tinker[11]等。這些軟件從算法選擇上可以歸為基于圖論[12]、基于化學(xué)計(jì)量學(xué)[13]和基于逆合成[14]等不同類別。利用這些工具，研究人員可方便地從代謝數(shù)據(jù)中識別、設(shè)計(jì)出具有潛在價(jià)值的代謝通路[15]。

搜索算法是代謝通路設(shè)計(jì)工具的核心步驟。在代謝網(wǎng)絡(luò)較大時(shí)，現(xiàn)有的代謝設(shè)計(jì)工具存在搜索時(shí)間相對較長的問題。為此，現(xiàn)有的代謝設(shè)計(jì)工具嘗試從數(shù)據(jù)結(jié)構(gòu)上對算法進(jìn)行優(yōu)化[16-17]，或在路徑選取策略上進(jìn)行優(yōu)化[18]。此外，在多路徑搜索算法方面，傳統(tǒng)的前K條最短路徑算法(K-Shortest Path，KSP)是以 Yen’s 算法為代表。該算法是基于 Dijkstra 算法進(jìn)行演化的版本，在搜索功能上實(shí)現(xiàn)了較為有效的前K條最短路徑搜索，但其算法過程從第 3 條路徑開始則需要提升。另外，綜合來看 KSP 算法的核心思想都是邊刪除選擇法，而過去眾多算法圍繞數(shù)據(jù)結(jié)構(gòu)及路徑選取策略上進(jìn)行過優(yōu)化，也取得了一定的性能提升效果。近年來，由于計(jì)算機(jī)硬件性能的提升，傳統(tǒng) KSP 算法優(yōu)勢變得不再明顯。因此，本文基于 Yen’s 算法進(jìn)行優(yōu)化改進(jìn)的 KSP 算法，在對路徑選取策略進(jìn)行優(yōu)化的同時(shí)，結(jié)合多核心CPU 計(jì)算機(jī)硬件特點(diǎn)，使用并行計(jì)算方式來進(jìn)一步提升代謝通路搜索算法的運(yùn)算性能。

2 數(shù)據(jù)建模

本文主要是針對合成生物學(xué)的代謝網(wǎng)絡(luò)路徑搜索進(jìn)行優(yōu)化，代謝網(wǎng)絡(luò)數(shù)據(jù)來源于 KEGG 通路數(shù)據(jù)庫(https://www.kegg.jp/kegg/pathway.html)。通過 KEGG API 接口批量下載不同物種及參考的代謝反應(yīng)圖的 KGML 文件，并進(jìn)行解析來提取化合物和反應(yīng)方向。將代謝通路表示為圖，圖中的節(jié)點(diǎn)表示為代謝通路中的化合物(代謝物)，邊表示化合物之間的轉(zhuǎn)化(化學(xué)反應(yīng)或代謝反應(yīng))。本文設(shè)計(jì)的代謝網(wǎng)絡(luò)路徑搜索模型針對 KEGG數(shù)據(jù)格式的特點(diǎn)進(jìn)行了專門優(yōu)化。在不考慮能耗和確保路徑正確的前提下，針對可能存在的代謝通路，盡可能搜索多條路徑，把代謝網(wǎng)絡(luò)數(shù)據(jù)簡化以提高搜索的效率。具體處理方法為：把具有通路的化合物之間的 Cost 值設(shè)為 1，不連通的化合物之間的 Cost 值設(shè)為 max(實(shí)驗(yàn)時(shí)，考慮到KEGG 代謝網(wǎng)絡(luò)參考圖中化合物數(shù)量約為 3 000個(gè)，設(shè) max＝9999)。所生成的代謝網(wǎng)絡(luò)圖的鄰接表如下：

3 算法設(shè)計(jì)

3.1 傳統(tǒng) KSP 算法

傳統(tǒng) KSP 算法主要以 Yen’s 算法為代表。其中，Yen’s 算法是 Yen 在 1971 年提出并以其名字命名的算法。Yen’s 算法采用了遞推法中的偏離路徑算法思想，適用于非負(fù)權(quán)邊的有向無環(huán)圖結(jié)構(gòu)。其主算法步驟如下：

(1)在非負(fù)權(quán)邊的有向無環(huán)圖G中，使用Dijkstra 算法找出第一條從起點(diǎn)S到終點(diǎn)D的最短路徑P1，加入集合 A(已選路徑集合)；

(2)從集合 A 中選擇最后加入的一條路徑Pn，在G中依次刪除路徑Pn中的邊(把邊的 Cost值設(shè)為最大)，并使用 Dijkstra 算法計(jì)算出對應(yīng)的最短路徑PNi加入集合 B(候選路徑集合)；

(3)從集合 B 中選出一條 Cost 最小的路徑Pm，加入集合 A 并從集合 B 中刪除Pm；

(4)重復(fù)步驟(2)～(3)直到集合 A 的數(shù)量達(dá)到K，則集合 A 為所求的前K條最短路徑。

本文將步驟(2)～(3)定義為次短路徑搜索(Second Shortest Path Search，SSPS)過程。在KSP 算法中，最耗時(shí)的運(yùn)算是 Dijkstra 運(yùn)算，而該算法中的 SSPS 過程存在較多重復(fù)的 Dijkstra計(jì)算。因此，在節(jié)點(diǎn)數(shù)量較多的網(wǎng)絡(luò)中進(jìn)行 KSP運(yùn)算將會耗費(fèi)較長時(shí)間。本文主要研究如何對SSPS 過程進(jìn)行優(yōu)化以減少不必要的 Dijkstra 運(yùn)算，從而提高 KSP 算法性能。

3.2 加速選擇法

在 Dijkstra 路徑搜索算法中，一次只能找出一條最短路徑，但 Cost 相同的最短路徑可能并不唯一。在 SSPS 過程前，集合 B 中可能已經(jīng)存在所需要的最短路徑。在合成生物代謝網(wǎng)絡(luò)中，所有連通的邊的 Cost 值設(shè)為 1，此時(shí) Cost 值相同的路徑都認(rèn)可為是最短路徑，但并不設(shè)先后順序。因此，若在集合 B 中存在與集合 A 最后一條路徑 Cost 值相同的路徑，則可認(rèn)為是在集合 A之后的最短路徑。此時(shí)，可直接從集合 B 中選擇出最短路徑Pn，則省去了 SSPS 的計(jì)算過程。

如圖 1 所示，當(dāng)前需要搜索第 4 條最短路徑時(shí)集合 A 中已有 3 條路徑，集合 B 中 Cost 值最小的路徑與集合 A 的最后一條路徑 Cost 值相等，此時(shí)直接從集合 B 中選取 Cost 值為 5 的路徑加入集合 A 中作為第 4 條最短路徑，從而節(jié)省了一輪 SSPS 運(yùn)算(此例子中可節(jié)省 4 次 Dijkstra運(yùn)算)。

圖1 加速選擇法Fig.1 Accelerated selection method

3.3 刪除重復(fù)候選路徑

在一個(gè) SSPS 過程中，得出的候選路徑可能與上一輪 SSPS 得到的候選路徑存在相同的，此時(shí)需要判斷，當(dāng)路徑已存在集合 B 中則不加入到集合 B。如圖 2 所示，刪除邊BC或CD所產(chǎn)生的候選路徑是相同的，此時(shí)只選擇一條加入集合B 即可。

圖2 相同的候選路徑Fig.2 The same candidate pathways

3.4 記錄已選邊

在已找到的路徑集合 A 中，除第一條路徑外，其他路徑都有對應(yīng)刪除的邊，這些邊的集合記為 EA。在進(jìn)行 SSPS 過程前，先在圖G中刪除 EA 里的邊，可減少重復(fù)的 Dijkstra 運(yùn)算。

如圖 3 所示，集合 A 中第二條最短路徑所對應(yīng)的邊為BD，在進(jìn)行下一輪 SSPS 運(yùn)算前把邊BD刪除以減少 SSPS 的運(yùn)算時(shí)間(此例子中減少了 1 次 Dijkstra 運(yùn)算)。

圖3 記錄已選邊Fig.3 Recording the selected edges

3.5 記錄關(guān)鍵邊

關(guān)鍵邊(Critical Edge)是從起點(diǎn)到終點(diǎn)的必經(jīng)之路，在一個(gè) SSPS 過程中跳過關(guān)鍵邊的 Dijkstra運(yùn)算可節(jié)省運(yùn)算時(shí)間。具體做法為：在 SSPS 過程中，若 Dijkstra 運(yùn)算的結(jié)果 Cost 是最大值，則把對應(yīng)的邊加入到關(guān)鍵邊的集合 C 中；隨后在下一輪 SSPS 過程中跳過關(guān)鍵邊的 Dijkstra 運(yùn)算，可節(jié)省運(yùn)算時(shí)間。

如圖 4 所示，當(dāng)經(jīng)過一輪 SSPS 運(yùn)算后，可以得到邊AB和DH為關(guān)鍵邊，并在下一輪 SSPS運(yùn)算時(shí)把關(guān)鍵邊排除，以減少運(yùn)算量(此例子中節(jié)省 2 次 Dijkstra 運(yùn)算)。

圖4 記錄關(guān)鍵邊Fig.4 Recording the key edges

4 平臺實(shí)現(xiàn)

平臺采用 Flask＋Vue 框架，實(shí)現(xiàn)前后端分離。核心算法代碼采用 C++實(shí)現(xiàn)以提升運(yùn)算速率。代謝網(wǎng)絡(luò)數(shù)據(jù)主要來源于 KEGG 數(shù)據(jù)庫，在路徑搜索前把 KEGG 數(shù)據(jù)進(jìn)行初始化(KEGG 數(shù)據(jù)庫轉(zhuǎn)化為鄰接表)以便進(jìn)行 KSP 算法運(yùn)算。系統(tǒng)主要程序流程如圖 5 所示。KEGG 數(shù)據(jù)的初始化在系統(tǒng)啟動(dòng)前期執(zhí)行，此部分程序運(yùn)行時(shí)間約為 3～5 s，這對用戶來說并不產(chǎn)生影響。對系統(tǒng)用戶產(chǎn)生影響的過程為 KSP 運(yùn)算部分，這也是本文主要的研究工作。

圖5 程序主流程Fig.5 The flow chart of the program

圖6 所示為改進(jìn)的 KSP 算法流程，其中SSPS 過程采用了并行方式同時(shí)進(jìn)行多個(gè) Dijkstra運(yùn)算，算法時(shí)間復(fù)雜度從O(n×m)減少到O(n)，算法結(jié)束后會產(chǎn)生K條最短路徑。為提高用戶體驗(yàn)，在實(shí)際應(yīng)用系統(tǒng)中每產(chǎn)生一條最短路徑則返回到用戶界面上顯示，這樣可以減少用戶主觀的等待時(shí)間。圖 7 所示用戶系統(tǒng)界面，搜索得到的路徑為按短到長依次顯示。

圖6 改進(jìn)的 KSP 算法流程Fig.6 The flow chart of improved KSP algorithm

圖7 系統(tǒng)界面Fig.7 The system interface

5 算法性能測試

在實(shí)驗(yàn)數(shù)據(jù)中，隨機(jī)選取起點(diǎn)和終點(diǎn)，分別設(shè)定最短路徑數(shù)量為 1～7 條，觀察算法所調(diào)用 Dijkstra 的次數(shù)來分析算法的性能。具體數(shù)據(jù)如表 1 所示。

從表 1 可知，當(dāng)路徑數(shù)量K大于 2 時(shí)，改進(jìn)的 KSP 算法在性能上有 2～3 倍的提升，并行的KSP 算法在性能上有 5～9 倍的提升，并隨著路徑數(shù)量的增加而不斷提升。本次改進(jìn)的 KSP 算法性能提升較為顯著。

表1 改進(jìn)的 KSP 算法性能對比Table 1 The performance of improved KSP algorithm

表2 數(shù)據(jù)為隨機(jī)選取 4 條目標(biāo)路徑進(jìn)行前 7條最短路徑的搜索，每條路徑的平均搜索時(shí)間為0.27～0.39 s。其中，所選取的 KEGG 參考代謝網(wǎng)絡(luò)圖節(jié)點(diǎn)總數(shù)為 2 828，邊總數(shù)為 4 051。即使在 KSP 算法運(yùn)算時(shí)間上增加 0.5 s 的 http 請求及獲取化合物及反應(yīng)信息所需要的時(shí)間，在用戶端平均每條路徑的搜索總時(shí)間仍然小于 1 s，提升了系統(tǒng)的用戶體驗(yàn)。

表2 改進(jìn)的 KSP 算法性能實(shí)測數(shù)據(jù)Table 2 The performance data of improved KSP algorithm

6 討論與分析

在合成生物學(xué)設(shè)計(jì)尤其是細(xì)胞工廠/底盤細(xì)胞設(shè)計(jì)中，代謝通路設(shè)計(jì)工具可有效輔助研究人員快速找到出發(fā)底物及目標(biāo)產(chǎn)物之間的連接路徑，提高設(shè)計(jì)效率。例如，Yim 等[19]對 1,4-丁二醇(Butane-1,4-diol，BDO)的合成代謝通路進(jìn)行 1萬多次的計(jì)算調(diào)查后，設(shè)計(jì)獲得最佳途徑，并在大腸桿菌中獲得高達(dá) 18 g/L 的 BDO。隨著下游酶的改善，BDO 滴度提高到 110 g/L。該案例成功展示了代謝通路設(shè)計(jì)算法在生物合成中的巨大應(yīng)用潛力。

本文主要針對基于圖論的路徑搜索算法進(jìn)行性能優(yōu)化，用于代謝通路的搜索和發(fā)掘。該類算法通過代謝網(wǎng)絡(luò)中的化合物和化學(xué)反應(yīng)之間的連接關(guān)系來搜索可行的代謝通路。目前，該類算法的代表軟件主要有 MRSD[4]、FogLight[10]和Metabolic Tinker[11]等。其主要優(yōu)點(diǎn)是可以在單一物種或者跨物種中尋找可行的代謝通路，不受物種和反應(yīng)流平衡等約束；缺點(diǎn)是在找到的代謝通路中容易出現(xiàn)一些連通度高的簇代謝物，而這些簇代謝物的存在會影響整體路徑的生物化學(xué)意義[11]。

在圖論中，路徑搜索算法主要有 Dijkstra 算法、A*算法、Bellman-Ford 算法、Floyd-Warshall算法和 Johnson 算法等。其中，A*算法在有估價(jià)函數(shù)的條件下可以快速搜索目標(biāo)路徑；Bellman-Ford 算法可以處理含負(fù)邊值的路徑；Floyd-Warshall 算法實(shí)現(xiàn)代碼簡單；Johnson 算法在 Bellman-Ford 算法的基礎(chǔ)上優(yōu)化并提高了稀疏圖的運(yùn)算效率；結(jié)合生物學(xué)代謝網(wǎng)絡(luò)中為搜索單源非負(fù)邊的路徑，Dijkstra 算法在時(shí)間復(fù)雜度上更有優(yōu)勢。自 Dijkstra 算法提出以來，許多學(xué)者都對它進(jìn)行過不同程度的優(yōu)化以提高其性能。在后來的多路徑搜索(KSP)算法上，大部分研究圍繞刪除路徑核心思想進(jìn)行了許多的改進(jìn)[16-18]。這些改進(jìn)的算法在過去計(jì)算機(jī)性能有限的條件下取得了較好的效果，算法時(shí)間復(fù)雜度從O(n×n)→O(n×m)→O(n×logm)不斷提升。但隨著計(jì)算機(jī)性能的提升，特別是多核心 CPU 以及多CPU 架構(gòu)的計(jì)算機(jī)系統(tǒng)的廣泛應(yīng)用，傳統(tǒng)以單線程進(jìn)行算法迭代運(yùn)算的方式并不能發(fā)揮很好的效果。本文正是利用了多核心 CPU 的硬件特點(diǎn)，在算法優(yōu)化的同時(shí)采用并行編程方式，把算法的時(shí)間復(fù)雜度提升到O(n)水平，大幅提升了運(yùn)算性能。

7 結(jié)論

本文針對合成生物學(xué)代謝網(wǎng)絡(luò)中代謝通路非唯一以及傳統(tǒng) KSP 算法效率偏低的問題，對 KSP算法進(jìn)行改進(jìn)優(yōu)化以提升運(yùn)算性能。在對 KSP 算法策略上優(yōu)化的同時(shí)，采用并行計(jì)算方式進(jìn)一步提升算法的性能。使用 Python 實(shí)現(xiàn)代謝通路設(shè)計(jì)Web 平臺，并采用 C++編寫核心算法的代碼。通過引入 KEGG 代謝網(wǎng)絡(luò)圖，驗(yàn)證改進(jìn)的 KSP 算法比傳統(tǒng) KSP 算法有較大的性能提升，在合成生物學(xué)代謝通路設(shè)計(jì)上具有一定的應(yīng)用價(jià)值。然而，由于代謝反應(yīng)在不同物種的代謝網(wǎng)絡(luò)中會有差別，本文基于 KEGG 參考圖上進(jìn)行了搜索算法的研究，并未對物種的特性進(jìn)行區(qū)分。因此，后續(xù)工作可在 KEGG 參考圖基礎(chǔ)上針對不同物種的約束條件進(jìn)行路徑算法研究，以適應(yīng)合成生物學(xué)對不同物種代謝網(wǎng)絡(luò)通路設(shè)計(jì)的特定需求。