基于多路徑廣度優(yōu)先搜索算法的代謝通路設計與實現(xiàn)

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

1(廣州中國科學院先進技術研究所 機器人與智能裝備研究中心 廣州 511458)

2(廣州中國科學院先進技術研究所 生物工程研究中心 廣州 511458)

3(中國科學院深圳先進技術研究院 深圳合成生物學創(chuàng)新研究院 中國科學院定量工程生物學重點實驗室廣東省合成基因組學重點實驗室 深圳 518055)

1 引 言

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

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

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

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

3 算法設計

3.1 傳統(tǒng) KSP 算法

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

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

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

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

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

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

3.2 加速選擇法

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

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

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

3.3 刪除重復候選路徑

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

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

3.4 記錄已選邊

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

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

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

3.5 記錄關鍵邊

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

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

圖4 記錄關鍵邊Fig.4 Recording the key edges

4 平臺實現(xiàn)

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

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

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

圖6 改進的 KSP 算法流程Fig.6 The flow chart of improved KSP algorithm

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

5 算法性能測試

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

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

表1 改進的 KSP 算法性能對比Table 1 The performance of improved KSP algorithm

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

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

6 討論與分析

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

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

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

7 結(jié)論

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