李培順，馬紅武，趙學明，陳濤

1 天津大學化工學院生物工程系，天津 3000722 天津大學教育部系統(tǒng)生物工程重點實驗室，天津 3000723 中國科學院天津工業(yè)生物技術(shù)研究所 中國科學院系統(tǒng)微生物技術(shù)重點實驗室，天津 300308

?

基于代謝網(wǎng)絡(luò)預測菌種基因改造靶點方法的研究進展

李培順1,2,3，馬紅武3，趙學明1,2，陳濤1,2

1 天津大學化工學院生物工程系，天津 300072
2 天津大學教育部系統(tǒng)生物工程重點實驗室，天津 300072
3 中國科學院天津工業(yè)生物技術(shù)研究所 中國科學院系統(tǒng)微生物技術(shù)重點實驗室，天津 300308

李培順, 馬紅武, 趙學明, 等. 基于代謝網(wǎng)絡(luò)預測菌種基因改造靶點方法的研究進展. 生物工程學報, 2016, 32(1): 1–13.

Li PS, Ma HW, Zhao XM, et al. Predicting genetic modification targets based on metabolic network analysis–a review. Chin J Biotech, 2016, 32(1): 1–13.

摘 要:高產(chǎn)特定產(chǎn)品的人工細胞工廠的構(gòu)建需要對野生菌株進行大量的基因工程改造，近年來隨著大量基因組尺度代謝網(wǎng)絡(luò)模型的構(gòu)建，人們提出了多種基于代謝網(wǎng)絡(luò)分析預測基因改造靶點以使某一目標化合物合成最優(yōu)的方法。這些方法利用基因組尺度代謝網(wǎng)絡(luò)模型中的反應計量關(guān)系約束和反應不可逆性約束等，通過約束優(yōu)化的方法預測可使產(chǎn)物合成最大化的改造靶點，避免了傳統(tǒng)的通過相關(guān)途徑的直觀分析確定靶點的方法的局限性和主觀性，為細胞工廠的理性設(shè)計提供了新的思路。以下結(jié)合作者的實際研究經(jīng)驗，對這些菌種優(yōu)化方法的原理、優(yōu)缺點及適用性等進行詳細介紹，并討論了目前存在的主要問題和未來的研究方向，為人們針對不同目標產(chǎn)品選擇合適的方法及預測結(jié)果的可靠性評估提供了指導。

關(guān)鍵詞:基因組尺度，代謝網(wǎng)絡(luò)，菌種優(yōu)化，系統(tǒng)生物學，代謝工程

Received: March 19, 2015; Accepted: May 20, 2015

Supported by: National Basic Research Program of China (973 Program) (Nos. 2012CB725203, 2011CBA00804), National High Technology Research and Development Program of China (863 Program) (No. 2012AA022103), Applied Basic Research Program of Tianjin (No. 12JCYBJC33000).

Tao Chen. Tel: +86-22-27406770; E-mail: taochen@tju.edu.cn

國家重點基礎(chǔ)研究計劃 (973計劃) (Nos. 2012CB725203, 2011CBA00804)，國家高技術(shù)研究發(fā)展計劃 (863 計劃) (No. 2012AA022103)，天津市應用基礎(chǔ)及前沿技術(shù)研究計劃項目 (No. 12JCYBJC33000) 資助。

生物技術(shù)正在被廣泛地應用于一些工業(yè)產(chǎn)品例如抗生素、生物能源、生物材料、大宗化學品等的生產(chǎn)。生物法相對于傳統(tǒng)的化學法具有反應條件溫和、節(jié)約能源、轉(zhuǎn)化率高和環(huán)境友好等優(yōu)勢。為了進行高效生產(chǎn)，通常需要對工業(yè)微生物進行一系列基因改造以提高目標產(chǎn)物的生產(chǎn)速率及得率。代謝工程主要任務(wù)之一就是識別工業(yè)微生物基因改造的靶點以達到高產(chǎn)目標化合物的目的。早期基因改造靶點的確定主要依靠對局部代謝途徑的分析和實驗的經(jīng)驗。隨著系統(tǒng)生物學和合成生物學的發(fā)展，在基因組測序和注釋海量數(shù)據(jù)的基礎(chǔ)上各種微生物的基因組尺度代謝網(wǎng)絡(luò)模型相繼被重構(gòu)出來[1-5]。它們已經(jīng)成為研究生物代謝系統(tǒng)不可缺少的工具，并廣泛應用于預測能提高目標化合物生產(chǎn)的微生物代謝工程改造靶點?；诩s束的表型模擬方法，比如通量平衡分析 (FBA)[6-7]，能計算出特定網(wǎng)絡(luò)中使目標產(chǎn)品得率最大的最優(yōu)途徑，從而確定基因過表達靶點。相比之下，通過代謝網(wǎng)絡(luò)分析確定基因敲除靶點要更困難一些。

最早用來預測基因敲除靶點的方法是最小化代謝調(diào)整 (MOMA) 算法[8]，它假設(shè)經(jīng)過基因敲除的突變菌在達到穩(wěn)態(tài)過程中會受到最小化代謝調(diào)整的影響而使新代謝通量分布與野生型菌株的通量分布之間的歐幾里得距離最小，由它求得的解為滿足約束條件的次優(yōu)解。MOMA算法預測到的基因敲除靶點成功應用于畢赤酵母Pichia pastoris的基因工程改造以高效生產(chǎn)胞質(zhì)人類超氧化物歧化酶 (hSOD)。靶點乙醇脫氫酶基因adh2的敲除使hSOD的產(chǎn)量提高了20%，并且沒有減少細胞生長。另外，該方法預測到的敲除靶點成功地指導了大腸桿菌Escherichia coli代謝工程改造以過量生產(chǎn)番茄紅素、纈氨酸及聚乳酸等產(chǎn)品。MOMA算法需要先選擇幾個基因或基因組合進行敲除然后評估其對生長和產(chǎn)物合成的影響。對包含上千個基因的代謝網(wǎng)絡(luò)多基因組合靶點的計算非常困難，因此MOMA更多用于在確定敲除靶點后評估其對代謝通量分布的影響，并不適用于從一復雜網(wǎng)絡(luò)中找出使產(chǎn)物生成最大化的敲除靶點組合。針對這一問題人們提出了一些優(yōu)化目標化合物生產(chǎn)的敲除靶點預測方法，例如OptKnock[9]、OptGene[10]等。大多數(shù)方法是通過細胞生長和產(chǎn)物合成的雙層優(yōu)化找到將生長和產(chǎn)物合成相偶聯(lián)的最優(yōu)基因敲除策略，已成功用于指導大腸桿菌和釀酒酵母Saccharomyces cerevisiae的代謝工程改造以優(yōu)化部分重要生化產(chǎn)品的生產(chǎn)。本文將結(jié)合作者的研究工作，重點針對這類方法詳細闡述其基本原理、特征及優(yōu)缺點，并介紹幾種在優(yōu)化中常用的軟件。

1　OptKnock及雙層優(yōu)化基本思路

2003年Burgard等[9]提出了第一個基于基因組尺度代謝網(wǎng)絡(luò)模型進行敲除靶點預測的雙層優(yōu)化方法OptKnock。該方法利用雙層混合整型線性規(guī)劃 (MILP) 來進行優(yōu)化計算，如圖1所示。該雙層結(jié)構(gòu)中內(nèi)層問題是生物量最大化時的約束線性規(guī)劃問題，包括化學計量學、反應通量上下限、熱力學、反應敲除等約束條件。外層問題則是以目標化合物合成最大化為目標的最優(yōu)化問題，確定哪些反應敲除可以在滿足內(nèi)層優(yōu)化的同時也可使產(chǎn)物生成最大化。

圖1　OptKnock方法的雙層結(jié)構(gòu)Fig. 1 The bilevel structure of OptKnock.

為了減少計算時間和提高預測敲除靶點的效率，OptKnock方法預測的敲除靶點數(shù)目一般限制在5個以內(nèi)，并且在進行計算之前需要對基因組尺度代謝網(wǎng)絡(luò)進行預處理[11]，主要包括以下步驟：1) 通過通量可變性分析 (Flux variability analysis，F(xiàn)VA)[12]，除去模型中通量始終為零的反應，并縮小反應的通量變化范圍。2) 為了減小候選敲除靶點的范圍，需除去以下4類反應：a) 干、濕實驗中對于生物量生長是必需的；b) 非基因關(guān)聯(lián)的反應、自發(fā)反應和擴散反應，實際上這些反應在濕實驗中是無法被敲除的；c) 特定子系統(tǒng)的反應，比如細胞膜生物合成、膜脂質(zhì)代謝、無機離子運輸?shù)?；d) 對于偶聯(lián)反應子集，只對其中一個反應進行分析，因為敲除該子集中任何一個反應都是等效的。

基于大腸桿菌代謝網(wǎng)絡(luò)模型，利用OptKnock方法以琥珀酸、乳酸及1,3-丙二醇等產(chǎn)品為目標進行優(yōu)化計算，能得到使生物量生長和目標產(chǎn)品合成相偶聯(lián)的靶點組合。以琥珀酸為例，OptKnock預測到的敲除靶點基因組合為乙醇脫氫酶、乳酸脫氫酶、磷酸葡萄糖異構(gòu)酶和丙酮酸甲酸裂解酶。敲除這些靶點后在無氧條件最大化細胞生長時必須生成琥珀酸，即通過基因敲除使琥珀酸成為細胞生長時的一種強制性產(chǎn)物，并且不會導致細胞無法生長。對此可以由圖2中琥珀酸合成速率與生長速率間的關(guān)系予以說明。圖中虛線為野生菌在設(shè)定生長速率為最大生長速率范圍內(nèi)任一值時通過FBA計算得到的琥珀酸最大和最小生成速率。其最小速率始終為0，而最大速率則隨著生長速率的增加而不斷減小，這表明此時琥珀酸生成和生長是一種競爭關(guān)系，因此琥珀酸生成對生長不利，細胞最優(yōu)化生長的結(jié)果將使得副產(chǎn)物琥珀酸生成為0。而在敲除OptKnock計算得到的4個靶點基因后再進行計算的結(jié)果如圖中實線所示。此時最優(yōu)生長速率降低，并且在生長速率為0.052 h–1到0.09 h–1范圍內(nèi)琥珀酸最小生成速率隨生長速率的增大而增大，在最優(yōu)生長速率時琥珀酸最大和最小生成速率曲線交匯到一點 (圓圈所示)。這意味著要最大化細胞生長就必須要生成琥珀酸，琥珀酸生產(chǎn)與細胞生長相偶聯(lián)。因此對突變株細胞自身調(diào)控優(yōu)化生長的同時也就優(yōu)化了琥珀酸的生成。這正是OptKnock等雙層優(yōu)化方法進行靶點預測的生物學基礎(chǔ)。對應用OptKnock方法針對一些目標產(chǎn)品，如琥珀酸、乳酸、1,3-丙二醇等預測得到的部分敲除靶點人們已經(jīng)進行了濕實驗驗證，能夠顯著提高目標化合物的產(chǎn)量[13]。

圖2　琥珀酸最大和最小生產(chǎn)速率隨生長速率的變化Fig. 2 Maximal and minimal succinate production rates at different gowth rates.

2　OptKnock衍生的菌種優(yōu)化方法

2.1預測敲除靶點的方法

2.1.1OptGene

由于基因組尺度代謝網(wǎng)絡(luò)的復雜性、冗余性和雙層結(jié)構(gòu)的混合整型線性規(guī)劃問題求解面臨著巨大的挑戰(zhàn)，導致OptKnock方法搜索到敲除靶點的計算時間過長，甚至搜索不到結(jié)果?；诖?，Patil等[10]利用遺傳算法 (GA) 提出了一種新的菌種優(yōu)化方法OptGene。遺傳算法是一種通過模擬達爾文生物進化論的原理搜索最優(yōu)解的方法，能夠在相對較短的時間內(nèi)為超大問題提供近似最優(yōu)解，但是很多時候無法得到全局最優(yōu)解，而是得到在預設(shè)的時間或者突變次數(shù)達到之后的局部最優(yōu)解，并且OptGene較OptKnock沒有速度上的提升。不過，OptGene相對OptKnock的優(yōu)勢是它既能夠進行反應水平的敲除，也能夠進行基因水平的敲除。因此，它不需要像OptKnock一樣排除自發(fā)反應，只需要確定哪些基因是可以敲除的。另外，OptGene還能優(yōu)化包含非線性目標函數(shù)的問題。

基于釀酒酵母代謝網(wǎng)絡(luò)模型，利用OptGene方法以琥珀酸、甘油及香草醛為目標進行優(yōu)化計算，預測到的靶點能將生物量生長和這些目標產(chǎn)品的生產(chǎn)進行偶聯(lián)。

2.1.2RobustKnock

隨著基因組尺度代謝網(wǎng)絡(luò)維度的增大，代謝網(wǎng)絡(luò)中會存在與目標化合物最優(yōu)合成相競爭的途徑，它們的存在能夠使代謝通量偏離目標化合物，從而導致較低的生產(chǎn)速率，甚至得不到目標產(chǎn)物。例如，OptKnock方法基于大腸桿菌基因組尺度代謝網(wǎng)絡(luò)模型iJO1366[14]以乳酸為目標進行優(yōu)化計算，由于競爭性途徑的存在，預測到的靶點會導致乳酸的生產(chǎn)速率依然為0。2009年，Tepper等提出了能夠識別代謝網(wǎng)絡(luò)中競爭性途徑的菌種優(yōu)化方法RobustKnock[15]。與OptKnock不同，RobustKnock外層問題是最大化目標化合物最小生成速率的優(yōu)化問題，從而識別到的敲除靶點能夠使目標化合物的生產(chǎn)成為生物量生長時的一種強制性產(chǎn)物。RobustKnock與OptKnock優(yōu)化方法預測敲除靶點的差異可以由圖3中的簡單例子說明，該圖中簡化的網(wǎng)絡(luò)模型由5個代謝物 (M1–M5) 組成，Vuptake、Vbiomass、Vproduct和Vby-product分別表示底物、生物量、目標產(chǎn)物和副產(chǎn)物的速率，由M1到M3的反應計量系數(shù)為2∶1，即2分子M1生成1分子M3，其他反應的計量系數(shù)關(guān)系都是1∶1。當以生物量生長速率Vbiomass為目標對該網(wǎng)絡(luò)模型進行通量平衡分析時，通量分布結(jié)果是生物量最優(yōu)生長速率Vbiomass等于Vuptake，而目標產(chǎn)物生成速率Vproduct為0，即此時代謝通量只經(jīng)過反應M1→M2。OptKnock基于該網(wǎng)絡(luò)模型以Vproduct為目標進行優(yōu)化計算，預測到的敲除靶點 (反應M1→M2) 能夠使目標化合物的最大生產(chǎn)速率達到0.5×Vuptake，但是當敲除反應M1→M2并以生物量生長速率為目標進行通量平衡分析時，目標化合物的生產(chǎn)速率可能為0，因為此時代謝通量會流向競爭性產(chǎn)物M5。而應用RobustKnock基于該網(wǎng)絡(luò)模型進行優(yōu)化計算預測到的敲除靶點則包括M1→M2和M4→M5兩個反應，能夠使目標產(chǎn)物的最小生產(chǎn)速率達到0.5×Vuptake，從而保證了目標產(chǎn)物的生成與生物量生長相偶聯(lián)。

圖3　RobustKnock與OptKnock預測敲除策略的差異Fig. 3 Difference of knockout strategies for RobustKnock and OptKnock.

基于大腸桿菌代謝網(wǎng)絡(luò)模型，應用RobustKnock方法以乙酸、甲酸、延胡索酸、乙醇、琥珀酸和乳酸這些產(chǎn)品為目標進行優(yōu)化計算，預測到的靶點都能將生物量生長和這些目標產(chǎn)品合成進行偶聯(lián)，這表明與OptKnock方法相比，RobustKnock方法預測的敲除靶點更具有魯棒性。

2.1.3OptSwap

微生物代謝過程中氧化還原酶 (如脫氫酶,還原酶) 通常與輔因子NAD (H) 和NADP (H)這兩個主要流通代謝物其中的一種有結(jié)合特異性。正是氧化還原酶的這種結(jié)合特異性才導致代謝功能的分工：將NAD+還原為NADH的酶能驅(qū)動氧化磷酸化，而將NADP+還原為NADPH的酶能驅(qū)動合成代謝反應，這樣生物系統(tǒng)才能調(diào)控資源流向能量生成或合成代謝。因此，優(yōu)化輔因子生成的策略在菌種優(yōu)化方法中也應該考慮。King 等開發(fā)的OptSwap[16]方法能通過預測氧化還原酶輔因子交換位點和敲除靶點來優(yōu)化目標產(chǎn)物的生產(chǎn)。該方法是基于RobustKnock的雙層MILP問題，并且添加氧化還原酶輔因子特異性交換作為約束條件，結(jié)合輔因子交換和反應敲除來獲得優(yōu)化目標產(chǎn)物的菌種改造靶點。

基于大腸桿菌基因組尺度代謝網(wǎng)絡(luò)模型iJO1366，OptSwap以丙氨酸、琥珀酸、乙酸及乳酸這些產(chǎn)品為目標進行優(yōu)化計算，預測到的改造靶點能將生物量生長和這些產(chǎn)品的生成相偶聯(lián)。針對丙氨酸和乳酸進行優(yōu)化計算時，同時敲除3個反應和交換1個氧化還原酶輔因子結(jié)合位點才能將這兩種產(chǎn)物的生成與生物量生長相偶聯(lián)，而只敲除4個或以下反應時都不能將兩者偶聯(lián)起來。針對琥珀酸和乙酸的優(yōu)化計算表明，同時考慮輔因子交換位點和敲除靶點的菌種優(yōu)化策略比僅僅考慮敲除靶點的效果要好。OptSwap方法拓展了能與生物量生長相偶聯(lián)的目標產(chǎn)品的范圍，但是比RobustKnock計算要求更高[16]。

2.1.4GDBB

為了解決雙層優(yōu)化方法計算時間過長的問題，2012年Egen等利用分支定界 (Truncated branch and bound) 算法提出了一種新的菌種優(yōu)化方法GDBB[17]。與OptKnock相比，該方法能夠在更短時間內(nèi)計算得到更多敲除數(shù)的基因敲除策略，但是它只能找到優(yōu)化問題的近似最優(yōu)解，很多時候無法得到全局最優(yōu)解。

基于大腸桿菌基因組尺度代謝網(wǎng)絡(luò)模型iAF1260，利用GDBB以乙酸和琥珀酸為目標進行優(yōu)化計算，預測到的靶點能將生物量生長和目標產(chǎn)物的生產(chǎn)進行偶聯(lián)，并且近似最優(yōu)的敲除策略分別在81 s和345 s時被找到。

2.1.5其他預測敲除靶點的方法

與OptGene方法的求解方式不同，GDLS[18]方法采用局部搜索的啟發(fā)式算法，在求解空間中進行有效、低復雜度的多路徑搜索，這樣可以降低計算時間。有研究表明經(jīng)過基因改造的突變菌事實上在達到穩(wěn)態(tài)過程中會受到最小化代謝調(diào)整的影響，從而使調(diào)整之后的代謝通量分布與野生型菌株的通量分布之間的改變最小，而并非以生物量生長為目標大幅調(diào)整通量分布。2013年Ren等基于此提出了一種新的雙層優(yōu)化方法MOMAKnock[19]，它的內(nèi)層問題是以最小化代謝調(diào)整為目標時的二次規(guī)劃問題。結(jié)果表明該方法能在最小化代謝調(diào)整通量分布情況下找出具有魯棒性的敲除靶點。另外，Xu等開發(fā)的ReacKnock[20]方法采用Karush-Kuhn-Tucker (KKT) 算法求解雙層混合整型線性規(guī)劃問題，這為雙層規(guī)劃問題的求解提供了一種新思路。最近，Zhang等結(jié)合LTM (Logic transformation of model) 方法與預測反應敲除的OptKnock雙層優(yōu)化方法開發(fā)了OptGeneKnock[21]，并且利用分支定界算法進行求解，不僅能直接進行基因水平的預測，還能較快地計算得到使目標化合物過量生產(chǎn)的近似最優(yōu)的敲除策略。

2.2預測敲除、上調(diào)和下調(diào)靶點的方法

2.2.1OptReg

由于微生物菌種的基因改造手段并不局限于基因的添加和完全敲除，而且研究表明基因表達的上調(diào)或下調(diào)會對胞內(nèi)代謝產(chǎn)生重要影響[22]，因此有必要開發(fā)能夠識別基因上調(diào)或下調(diào)靶點的菌種優(yōu)化方法。2006年，Pharkya等開發(fā)了一種新的菌種優(yōu)化方法OptReg[23]，該方法在OptKnock基礎(chǔ)上進一步整合了反應通量的調(diào)節(jié)機制，能夠同時識別令目標產(chǎn)物過量生產(chǎn)的上調(diào)、下調(diào)和敲除靶點，但是由于拓展了識別基因改造靶點的范圍，導致雙層規(guī)劃問題含有更多變量和非線性關(guān)系，進一步增加了轉(zhuǎn)化為可求解的單層優(yōu)化問題的難度，計算上面臨著巨大挑戰(zhàn)。假設(shè)某一反應Vj在模型中通量范圍為[0, 1 000]，而依據(jù)通量可變性分析 (FVA) 確定的反應上下限為[Vj,L, Vj,U]。在OptReg方法中，當限制Vj在 (Vj,U, 1 000] 區(qū)間時，表示該反應被上調(diào)；當限制Vj在 (0, Vj,L) 區(qū)間時，表示該反應被下調(diào)；當限制Vj為0時，則表示該反應被敲除，如圖4所示。

圖4　OptReg方法中反應上調(diào)、下調(diào)及敲除的定義[23]Fig. 4 Definitions of up/down regulation and knockout in OptReg[23].

基于大腸桿菌代謝網(wǎng)絡(luò)模型，OptReg方法以乙醇為目標進行優(yōu)化計算，預測到的靶點能將生物量生長和乙醇進行偶聯(lián)，研究表明菌種改造過程中反應敲除和調(diào)節(jié)之間存在著協(xié)同作用。

2.2.2OptORF

由于多功能酶、多亞基酶復合物及同工酶的存在，基因與反應之間的關(guān)系并非總是一一映射，導致基于反應敲除的突變菌株有時很難在實驗過程中被構(gòu)建。此外，微生物細胞代謝還會受到轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的影響。而大多數(shù)菌種優(yōu)化方法是基于反應敲除，而且沒有考慮轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)。為了克服以上缺點，2010年Kim等開發(fā)了OptORF[24]，該方法能利用基因組尺度代謝與轉(zhuǎn)錄調(diào)控整合網(wǎng)絡(luò)模型去識別代謝基因敲除、過表達和轉(zhuǎn)錄因子敲除的菌種改造靶點。OptORF方法應用于一個簡化的代謝與轉(zhuǎn)錄調(diào)控整合網(wǎng)絡(luò)的實例如圖5所示。網(wǎng)絡(luò)中底物S通過中間代謝物M1或M2能被轉(zhuǎn)化為生物量B。假設(shè)反應R2能將M1轉(zhuǎn)化為目標產(chǎn)物P1和0.08 B，而反應R5能將M2轉(zhuǎn)化為副產(chǎn)物P2和0.12 B。當?shù)孜颯存在，轉(zhuǎn)錄因子TF1呈激活狀態(tài)，由此激活基因G3和G5的表達，并抑制G1A的表達。最大化生物量生長時，代謝通量會流過反應R3和R5，因此會生成副產(chǎn)物P2，而不會生成目標產(chǎn)物P1。以P1為目標產(chǎn)品，利用OptORF方法對該網(wǎng)絡(luò)進行優(yōu)化計算，預測到的基因改造靶點包括敲除基因G3和G4 (關(guān)聯(lián)反應R3、R4) 或G5和G6 (關(guān)聯(lián)反應R5)。由于基因G1A的表達被轉(zhuǎn)錄因子TF1所抑制，導致反應R1無法進行。因此，OptORF在預測到上述敲除靶點的同時還會識別到過表達基因G1A。除此之外，OptORF也能預測到轉(zhuǎn)錄因子敲除靶點TF1，從而解除它對基因G1A的抑制作用，同時使基因G3和基因G5失活。

基于大腸桿菌代謝與轉(zhuǎn)錄調(diào)控整合網(wǎng)絡(luò)模型iMC1010v2[25]，OptORF方法以乙醇為目標產(chǎn)品進行優(yōu)化計算，預測到的改造靶點 (比如敲除基因pgi和過表達基因edd) 能將生物量生長和乙醇的生成進行偶聯(lián)。與OptKnock方法比較結(jié)果表明，基于反應敲除的策略通常要求敲除更多相關(guān)聯(lián)的基因以除去被識別的反應，并且當考慮轉(zhuǎn)錄調(diào)控時還可能會導致致死性生長表型。

圖5　OptORF應用于簡化的代謝與轉(zhuǎn)錄調(diào)控整合網(wǎng)絡(luò)Fig. 5 Application of OptORF to an example metabolic and regulatory network.

2.2.3OptForce

上述的雙層菌種優(yōu)化算法在設(shè)計時均假設(shè)基因改造后的代謝網(wǎng)絡(luò)模型會根據(jù)最大化細胞目標的思想去分配代謝通量。而OptForce[26]采用逆向的計算方法，假設(shè)當目標產(chǎn)物生產(chǎn)速率達到某一預設(shè)值時，突變型代謝網(wǎng)絡(luò)中反應的通量范圍與野生型進行比較，進而確定那些通量范圍偏差 (上調(diào)、下調(diào)或變?yōu)?) 較大的反應作為基因改造的靶點。與OptReg方法一樣，OptForce確定反應通量范圍也是基于FVA算法，依次以代謝網(wǎng)絡(luò)中每個反應為優(yōu)化目標，計算滿足目標產(chǎn)物生產(chǎn)速率和其他約束條件下每個反應通量的最大值和最小值。

2.2.4其他預測改造靶點的方法

2004年P(guān)harkya等提出了一種整合的菌種優(yōu)化方法OptStrain[27]，該方法首先通過構(gòu)建反應數(shù)據(jù)庫為微生物細胞添加最優(yōu)外源途徑從而賦予該微生物生產(chǎn)某一非天然化合物的能力，然后在加入外源合成途徑的代謝網(wǎng)絡(luò)中通過OptKnock算法預測到使該目標化合物過量生產(chǎn)的敲除靶點。

除了上述由OptKnock衍生過來的菌種優(yōu)化方法，還有一些基于FBA或FVA的代謝網(wǎng)絡(luò)分析預測改造靶點的方法，比如通量響應分析(Flux response analysis)[28]、FSEOF[29]、FVSEOF[30]和RobOKoD[31]等。其中，通量響應分析、FSEOF和FVSEOF方法僅用來預測過表達靶點。通量響應分析方法通過系統(tǒng)地考察目標產(chǎn)物反應通量隨代謝網(wǎng)絡(luò)模型中其他反應通量的改變的變化趨勢來確定反應之間的相互關(guān)系，從而識別出能夠提高目標產(chǎn)物產(chǎn)量的過表達靶點。類似地，F(xiàn)SEOF通過系統(tǒng)地考察代謝網(wǎng)絡(luò)模型中那些隨目標產(chǎn)物反應通量的增加而增加的反應作為過表達靶點。另外，F(xiàn)VSEOF方法利用FVA算法考慮了基于約束問題的多種最優(yōu)解和代謝網(wǎng)絡(luò)中各反應速率的范圍，克服了FSEOF方法在這方面的缺陷。同時該方法通過添加反應簇群 (GR) 約束來模擬細胞的生理狀態(tài)，有效提高了預測準確性。而RobOKoD最近由Stanford等[31]提出，該方法利用FVA針對生物量生長和目標化合物最大化時的反應進行分級評價，找出能夠使目標化合物過量生產(chǎn)的敲除、過表達和抑制靶點。

3　雙層優(yōu)化方法存在的問題

3.1雙層結(jié)構(gòu)及求解問題

微生物代謝網(wǎng)絡(luò)的行為會受到內(nèi)部細胞目標的調(diào)控，而一般情況下內(nèi)部細胞目標會與目標產(chǎn)物的過量生產(chǎn)相競爭。不過，基于基因組尺度代謝網(wǎng)絡(luò)模型，雙層菌種優(yōu)化方法能夠通過靶點的改造，比如反應或基因的敲除、上調(diào)及下調(diào)，將內(nèi)部細胞目標與部分目標產(chǎn)物的生產(chǎn)偶聯(lián)起來，從而實現(xiàn)目標產(chǎn)物的過量生產(chǎn)。這類方法主要通過雙層混合整型線性或非線性規(guī)劃來實現(xiàn)。正是由于此類雙層結(jié)構(gòu)的特點和基因組尺度代謝網(wǎng)絡(luò)的復雜性及冗余性，導致雙層菌種優(yōu)化方法直接求解是不可行的。因此，若要求解它們都是先利用線性規(guī)劃的對偶理論將內(nèi)層的原始問題轉(zhuǎn)化為對偶問題，然后再將雙層問題轉(zhuǎn)化為標準的單層混合整型線性規(guī)劃問題。而這類單層問題能利用各種數(shù)學優(yōu)化軟件中MILP求解器進行求解。由于求解空間的維度極大 (例如從150個候選目標反應集合中搜索10個敲除靶點，其組合數(shù)超過1015，實際上遠遠不止如此)，導致求解時間隨敲除數(shù)的增加而呈指數(shù)型增長。盡管GDBB方法能夠在較短時間內(nèi)給出計算結(jié)果，但是不能得到全局最優(yōu)解。大多數(shù)雙層菌種優(yōu)化方法能預測到的改造靶點數(shù)目一般不會超過5個。

3.2解的特點及問題

圖6　PHB最大和最小生產(chǎn)速率隨生長速率的變化Fig. 6 Maximal and minimal PHB production rates at different gowth rates.

雙層優(yōu)化方法試圖通過識別令目標化合物過量生產(chǎn)的改造靶點來實現(xiàn)代謝網(wǎng)絡(luò)的重新布局，將生物量生長與目標產(chǎn)物相偶聯(lián)，使目標化合物的生產(chǎn)成為生物量生長時的一種強制性副產(chǎn)物，即該產(chǎn)物的生成是平衡細胞生長的還原力和能量需求的唯一途徑，通過這類方法預測到的改造靶點不會使細胞無法生長而只是平衡生長與產(chǎn)物生成的需求。但是并非所有的代謝物其生成都可以和細胞生長偶聯(lián)，在某些情況下，雖然用雙層優(yōu)化方法可以計算得到敲除靶點，但對改造后的代謝網(wǎng)絡(luò)進行模擬的結(jié)果卻表明該敲除策略并不能使產(chǎn)物和生長相偶聯(lián)。例如，我們應用OptKnock針對聚羥基脂肪酸脂 (PHB) 進行優(yōu)化計算，得到的一組敲除靶點基因組合為乙酸激酶、琥珀酰輔酶A合成酶和磷酸丙糖異構(gòu)酶。然而敲除這些靶點后進行計算并不能得到類似圖2針對琥珀酸計算的結(jié)果，而是如圖6所示。在敲除改造的菌株中生長速率最大時PHB生成速率為0，并沒有與生長偶聯(lián)。我們還針對蘇氨酸、核黃素等產(chǎn)品進行了計算，發(fā)現(xiàn)也無法得到真正將產(chǎn)物生成和細胞生長相偶聯(lián)的敲除結(jié)果，因此這類雙層優(yōu)化方法僅適用于一些可以與生長相偶聯(lián)的代謝物的敲除靶點預測，以同時實現(xiàn)產(chǎn)物合成和細胞生長的最大化，對于無法與生長相偶聯(lián)的代謝物則需要從其它角度尋找新的敲除靶點預測方法。

4　菌種優(yōu)化軟件及工具包

早期的菌種優(yōu)化方法，比如OptKnock、OptStrain、OptReg以及后來的OptORF，都是基于通用代數(shù)建模系統(tǒng) (The general algebraic modeling system，GAMS) 這一計算平臺來實現(xiàn)運算，但都沒有提供源代碼。隨后在Matlab平臺的基礎(chǔ)上，一些菌種優(yōu)化方法相繼實現(xiàn)運算。其中，基于約束的重構(gòu)及分析工具，COBRA (Constraints based reconstruction and analysis)[32-33]工具箱是目前應用最為廣泛的分析軟件，2007年由Palsson課題組設(shè)計，不僅能實現(xiàn)基因組尺度代謝網(wǎng)絡(luò)的構(gòu)建，也能進一步利用代謝網(wǎng)絡(luò)對微生物的表型進行模擬、分析及預測。COBRA工具箱在菌種優(yōu)化設(shè)計方面集成的算法有MOMA、OptKnock、OptGene、GDLS等。另一個基于Matlab用于代謝網(wǎng)絡(luò)半自動化構(gòu)建、分析和可視化的工具包是RAVEN[34]，可以讀寫SBML和Excel格式的模型文件，集成的算法有FBA、MOMA等。此外，GDLS、GDBB、RobustKnock和OptSwap等方法也是基于Matlab平臺實現(xiàn)運算，并提供了工具包。2013年，Palsson研究組基于Python開發(fā)了COBRApy[35]工具包。該工具包支持基本的COBRA方法，具有面向?qū)ο蟮男问剑苓m合處理像代謝和基因表達這類復雜的生物過程，并且在數(shù)據(jù)讀寫以及FBA計算速度上更有優(yōu)勢。另一個基于Python用于代謝網(wǎng)絡(luò)構(gòu)建和分析的工具包是Metano[36]，集成的算法有FBA、FVA、MOMA等，不僅可以讀取SBML格式的文件，也可以讀寫txt文件。2013年，Gelius-Dietrich等基于R開發(fā)了Sybil[37]工具包，集成的算法有FBA、FVA、MOMA、ROOM等，能夠快速地進行基因組尺度計算分析，并且具有面向?qū)ο蟮男问?，很容易被用戶拓展。但是，COBRApy、Sybil、Metano和RAVEN等工具包只集成有FBA和幾種常用的代謝網(wǎng)絡(luò)構(gòu)建及分析方法，而預測靶點的方法只有MOMA，并沒有類似OptKnock這種雙層優(yōu)化方法。另外，OptFlux[38]是一種開放型資源和模塊化的軟件，并且是第一個整合了菌種優(yōu)化任務(wù)的工具，其中集成有OptKnock 和OptGene方法，具有用戶友好性，可以讓用戶在可視化界面中執(zhí)行某種操作，不需要在計算平臺上輸入指令，但是它并不能正確地預測到敲除靶點，因此有很大的局限性。常用的菌種優(yōu)化軟件及工具包如表1所示。雖然目前已有以上幾種菌種優(yōu)化的工具包及軟件被開發(fā)，并為代謝工程改造微生物工廠提供指導，但是無論是計算時間及可行性，還是預測的準確性都存在很大問題，亟待進一步改進和開發(fā)新型的工具。

5　菌種優(yōu)化方法總結(jié)及展望

近幾年，菌種優(yōu)化方法的研究發(fā)展非常迅速，不僅能基于基因組尺度代謝網(wǎng)絡(luò)預測到令目標化合物過量生產(chǎn)的基因改造靶點，而且相關(guān)技術(shù)、算法、各種軟件工具都有了很大的進步。然而目前仍然存在一些關(guān)鍵的問題暫時無法解決。

表1　常用的菌種優(yōu)化軟件和工具包Table 1 Software and toolboxes used for strain optimization

首先，雙層優(yōu)化方法只能預測到能夠與生物量生長相偶聯(lián)的目標產(chǎn)物的改造靶點，以實現(xiàn)生物量生長和產(chǎn)物合成這兩個目標的同時優(yōu)化，而大部分代謝物實際上無法與生物量生長相偶聯(lián)，所以針對這些代謝物該類方法無法預測到有效的改造靶點。

其次，現(xiàn)有的菌種優(yōu)化方法很少利用基因組尺度代謝網(wǎng)絡(luò)與其他組學網(wǎng)絡(luò)整合在一起的生物大網(wǎng)絡(luò)。在細胞內(nèi)存在著很多其他的重要機制，如轉(zhuǎn)錄調(diào)控、信號轉(zhuǎn)導等，如果整合這些組學網(wǎng)絡(luò)形成全細胞網(wǎng)絡(luò)，將大大提高菌種優(yōu)化方法對生物表型的預測能力，并成為代謝工程決策的有力武器。

再次，大多數(shù)雙層菌種優(yōu)化方法求解時間過長。由于基因組尺度代謝網(wǎng)絡(luò)的復雜性及冗余性，導致菌種優(yōu)化方法在計算上面臨著巨大的挑戰(zhàn)。盡管GDBB能夠在較短的時間內(nèi)給出計算結(jié)果，卻不能得到全局最優(yōu)解。然而，隨著并行計算技術(shù)的發(fā)展，相信未來能夠解決更大維度的問題，搜索到更多敲除數(shù)的改造策略。

最后，目前大多數(shù)菌種優(yōu)化方法是基于反應敲除來實現(xiàn)代謝網(wǎng)絡(luò)的重新布局。然而，細胞內(nèi)存在復雜的基因-蛋白質(zhì)-反應 (GPR) 映射關(guān)系，這樣從基因改造層面有時很難實現(xiàn)代謝網(wǎng)絡(luò)重新布局這一目標，甚至會導致產(chǎn)量降低或致死性表型。不過，OptGene、OptORF和OptGeneKnock等方法能夠直接做基因水平的預測，解決了這個問題。

隨著細胞整合網(wǎng)絡(luò)、高性能計算機和并行計算技術(shù)的不斷發(fā)展，菌種優(yōu)化方法預測的結(jié)果將更準確，更符合生物學意義，計算速度更快，這將會使基因改造的范圍變廣，更趨向于理性，大大拓展能夠與生物量生長相偶聯(lián)的化合物范圍。

REFERENCES

[1] Edwards JS, Palsson BO. The Escherichia coli MG1655 in silico metabolic genotype: its definition, characteristics, and capabilities. Proc Natl Acad Sci USA, 2000, 97(10): 5528–5533.

[2] Reed JL, Vo TD, Schilling CH, et al. An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR). Genome Biol, 2003, 4(9): R54.

[3] Feist AM, Henry CS, Reed JL, et al. A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Mol Syst Biol, 2007, 3(1): 121.

[4] Zhou MD, Zou W, Liu LM, et al. Reconstruction and analysis of Bacillus megaterium metabolic model based on literature study. Acta Microbiol Sin, 2012, 52(4): 457–465 (in Chinese).周冒達, 鄒偉, 劉立明, 等. 基于文獻挖掘的巨大芽胞桿菌代謝網(wǎng)絡(luò)模型的構(gòu)建與分析. 微生物學報, 2012, 52(4): 457–465.

[5] Chai WP, Xue W, Zhang L, et al. Research on the auto-reconstruction of genome-scale metabolic network model. J Food Sci Biotechnol, 2014, 33(9): 957–965 (in Chinese).柴文平, 薛衛(wèi), 張梁, 等. 基因組規(guī)模代謝網(wǎng)絡(luò)模型的自動化重構(gòu). 食品與生物技術(shù)學報, 2014, 33(9): 957–965.

[6] Kauffman KJ, Prakash P, Edwards JS. Advances in flux balance analysis. Curr Opin Biotechnol, 2003, 14(5): 491–496.

[7] Price ND, Papin JA, Schilling CH, et al. Genome-scale microbial in silico models: the constraints-based approach. Trends Biotechnol, 2003, 21(4): 162–169.

[8] Segrè D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA, 2002, 99(23): 15112–15117.

[9] Burgard AP, Pharkya P, Maranas CD. Optknock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng, 2003, 84(6): 647–657.

[10] Patil KR, Rocha I, F?rster J, et al. Evolutionary programming as a platform for in silico metabolic engineering. BMC Bioinformatics, 2005, 6: 308. [11] Feist AM, Zielinski DC, Orth JD, et al. Model-driven evaluation of the production potential for growth-coupled products of Escherichia coli. Metab Eng, 2010, 12(3): 173–186.

[12] Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng, 2003, 5(4): 264–276.

[13] Fong SS, Burgard AP, Herring CD, et al. In silico design and adaptive evolution of Escherichia coli for production of lactic acid. Biotechnol Bioeng, 2005, 91(5): 643–648.

[14] Orth JD, Conrad TM, Na J, et al. A comprehensive genome-scale reconstruction of Escherichia coli metabolism—2011. Mol Syst Biol, 2011, 7(1): 535.

[15] Tepper N, Shlomi T. Predicting metabolic engineering knockout strategies for chemical production: accounting for competing pathways. Bioinformatics, 2010, 26(4): 536–543.

[16] King ZA, Feist AM. Optimizing cofactor specificity of oxidoreductase enzymes for the generation of microbial production strains—OptSwap. Ind Biotechnol, 2013, 9(4): 236–246.

[17] Egen D, Lun DS. Truncated branch and bound achieves efficient constraint-based genetic design. Bioinformatics, 2012, 28(12): 1619–1623.

[18] Lun DS, Rockwell G, Guido NJ, et al. Large-scale identification of genetic design strategies using local search. Mol Syst Biol, 2009, 5(1): 296.

[19] Ren SG, Zeng B, Qian XN. Adaptive bi-level programming for optimal gene knockouts for targeted overproduction under phenotypic constraints. BMC Bioinformatics, 2013, 14(Suppl 2): S17.

[20] Xu ZX, Zheng P, Sun JB, et al. ReacKnock: identifying reaction deletion strategies for microbial strain optimization based on genome-scale metabolic network. PLoS ONE, 2013, 8(12): e72150.

[21] Zhang C, Ji BY, Mardinoglu A, et al. Logical transformation of genome-scale metabolic models for gene level applications and analysis. Bioinformatics, 2015, 31(14): 2324–2331.

[22] Koffas MAG, Jung GY, Stephanopoulos G. Engineering metabolism and product formation in Corynebacterium glutamicum by coordinated gene overexpression. Metab Eng, 2003, 5(1): 32–41.

[23] Pharkya P, Maranas CD. An optimization framework for identifying reaction activation/inhibition or elimination candidates for overproduction in microbial systems. Metab Eng, 2006, 8(1): 1–13.

[24] Kim J, Reed JL. OptORF: optimal metabolic and regulatory perturbations for metabolic engineering of microbial strains. BMC Syst Biol, 2010, 4: 53.

[25] Covert MW, Knight EM, Reed JL, et al. Integrating high-throughput and computational data elucidates bacterial networks. Nature, 2004, 429(6987): 92–96.

[26] Ranganathan S, Suthers PF, Maranas CD. OptForce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions. PLoS Comput Biol, 2010, 6(4): e1000744.

[27] Pharkya P, Burgard AP, Maranas CD. OptStrain: a computational framework for redesign of microbial production systems. Genome Res, 2004, 14(11): 2367–2376.

[28] Song CW, Kim DI, Choi S, et al. Metabolic engineering of Escherichia coli for the production of fumaric acid. Biotechnol Bioeng, 2013, 110(7): 2025–2034.

[29] Choi HS, Lee SY, Kim TY, et al. In silico identification of gene amplification targets for improvement of lycopene production. Appl Environ Microbiol, 2010, 76(10): 3097–3105.

[30] Park JM, Park HM, Kim WJ, et al. Flux variability scanning based on enforced objective flux for identifying gene amplification targets. BMC Syst Biol, 2012, 6(1): 106.

[31] Stanford NJ, Millard P, Swainston N. RobOKoD: microbial strain design for (over) production of target compounds. Front Cell Dev Biol, 2015, 3: 17.

[32] Becker SA, Feist AM, Mo ML, et al. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA toolbox. Nat Protoc, 2007, 2(3): 727–738.

[33] Schellenberger J, Que R, Fleming RMT, et al. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA toolbox v2.0. Nat Protoc, 2011, 6(9): 1290–1307.

[34] Agren R, Liu LM, Shoaie S, et al. The RAVEN toolbox and its use for generating a genome-scale metabolic model for Penicillium chrysogenum. PLoS Comput Biol, 2013, 9(3): e1002980.

[35] Ebrahim A, Lerman JA, Palsson BO, et al. COBRApy: COnstraints-based reconstruction and analysis for python. BMC Syst Biol, 2013, 7(1): 74.

[36] Ulas T, Riemer SA, Zaparty M, et al. Genome-scale reconstruction and analysis of the metabolic network in the hyperthermophilic archaeon Sulfolobus solfataricus. PLoS ONE, 2012, 7(8): e43401.

[37] Gelius-Dietrich G, Desouki AA, Fritzemeier CJ, et al. Sybil-efficient constraint-based modelling in R. BMC Syst Biol, 2013, 7(1): 125.

[38] Rocha I, Maia P, Evangelista P, et al. OptFlux: an open-source software platform for in silico metabolic engineering. BMC Syst Biol, 2010, 4: 45.

(本文責編 郝麗芳)

Predicting genetic modification targets based on metabolic network analysis–a review

Peishun Li1,2,3, Hongwu Ma3, Xueming Zhao1,2, and Tao Chen1,2
1 Department of Biochemical Engineering, School of Chemical Engineering & Technology, Tianjin University, Tianjin 300072, China
2 Key Laboratory of Systems Bioengineering, Ministry of Education, Tianjin University, Tianjin 300072, China
3 Key Laboratory of Systems Microbial Biotechnology, Tianjin Institute of Industrial Biotechnology, Chinese Academy of Sciences, Tianjin 300308, China

Abstract:Construction of artificial cell factory to produce specific compounds of interest needs wild strain to be genetically engineered. In recent years, with the reconstruction of many genome-scale metabolic networks, a number of methods have been proposed based on metabolic network analysis for predicting genetic modification targets that lead to overproduction of compounds of interest. These approaches use constraints of stoichiometry and reaction reversibility in genome-scale models of metabolism and adopt different mathematical algorithms to predict modification targets, and thus can discover new targets that are difficult to find through traditional intuitive methods. In this review, we introduce the principle, merit, demerit and application of various strain optimization methods in detail. The main problems in existing methods and perspectives on this emerging research field are also discussed, aiming to provide guidance to choose the appropriate methods according to different types of products and the reliability of the predicted results.

Keywords:genome-scale, metabolic network, strain optimization, systems biology, metabolic engineering

DOI:10.13345/j.cjb.150118

Corresponding authors: Hongwu Ma. Tel: +86-22-24828735; E-mail: ma_hw@tib.cas.cn