一種基于代謝網(wǎng)絡(luò)分析最小化基因組的方法及其在大腸桿菌中的應(yīng)用

湯彬彩，郝彤,4，袁倩倩,3，陳濤，馬紅武,3

1天津大學(xué)化工學(xué)院生物工程系，天津 300072

2天津大學(xué) 教育部系統(tǒng)生物工程重點(diǎn)實(shí)驗(yàn)室，天津 300072

3中國(guó)科學(xué)院天津工業(yè)生物技術(shù)研究所 中國(guó)科學(xué)院系統(tǒng)微生物技術(shù)重點(diǎn)實(shí)驗(yàn)室，天津 300308

4天津師范大學(xué)生命科學(xué)學(xué)院，天津 300387

隨著2010年5月美國(guó)生物學(xué)家文特爾等在Science雜志發(fā)表了關(guān)于由化學(xué)合成的基因組控制的細(xì)胞的文章[1]，人工合成生命的研究引起了廣泛的關(guān)注。構(gòu)建最小生命體一直是人工合成生命的重要內(nèi)容之一，文特爾的工作再次引發(fā)了研究者們對(duì)最小基因組的研究興趣。

最小基因組是指在最適宜的條件下，維持細(xì)胞生長(zhǎng)繁殖所必需的最小數(shù)目的基因，在最小基因組中，敲除任何基因后細(xì)胞都會(huì)死亡。在不同培養(yǎng)條件下，得到的最小基因組可能是不一樣的，這與生物本身的合成能力有關(guān)，所以最小基因組是相對(duì)而言的。在代謝工程研究領(lǐng)域中，人們關(guān)心更多的是如何得到一個(gè)最小化的基因組而同時(shí)保持細(xì)胞的生長(zhǎng)不受影響。

目前對(duì)最小基因組的研究在濕實(shí)驗(yàn)方面已經(jīng)取得了一些進(jìn)展，1995年起研究者開(kāi)始對(duì)生殖支原體進(jìn)行研究[2]，因?yàn)樯持гw是能夠自主復(fù)制并生存的最小細(xì)胞生物，其基因組只有580 kb，2007年Lartigue等將生殖支原體整個(gè)基因組轉(zhuǎn)移到掏空的山羊支原體里[3]并使之成活，2008年Gibson等人工合成整個(gè)生殖支原體[4]，這些工作體現(xiàn)著研究者們對(duì)獲得人工合成最小生命的不斷探索。除了“自下而上”地從頭合成基因組外，另一個(gè)研究思路是從已有物種出發(fā)，通過(guò)“自上而下”地逐步敲除基因，去除細(xì)胞中的冗余基因，得到最小基因組。例如，Posfai等通過(guò)大片段基因敲除成功將Escherichia coliK-12基因組減少了15%[5]。同時(shí)這一小基因組大腸桿菌還具有一些意想不到的特性：比如高效的電轉(zhuǎn)化率，準(zhǔn)確的復(fù)制重組基因和穩(wěn)定的質(zhì)粒復(fù)制等。這揭示了在大腸桿菌中并非所有的基因都是必需的，通過(guò)大量刪除基因，可以獲得具有新性能的微生物。根據(jù)這一思路，必需基因[6-7]的確定是獲得最小基因組的重要依據(jù)。但最小基因組中除了這些必需基因外還需要包括其他的基因，這些基因雖然單獨(dú)敲除不會(huì)對(duì)細(xì)胞產(chǎn)生影響，但兩個(gè)或多個(gè)同時(shí)敲除時(shí)則可能導(dǎo)致細(xì)胞死亡。如何找出這些非必需但又需要包含在最小基因組中的基因是目前研究中面臨的一項(xiàng)重要挑戰(zhàn)。而基因組規(guī)模代謝網(wǎng)絡(luò)則有可能為解決這一問(wèn)題提供可能性。

基因組尺度的代謝網(wǎng)絡(luò)模型指的是以注釋信息及基因序列為基礎(chǔ)，利用基因-蛋白質(zhì)-反應(yīng)的關(guān)系為主線，成功模擬代謝物在生物機(jī)體內(nèi)的過(guò)程的一種包含大量信息的模型，到2012年6月為止已經(jīng)構(gòu)建完成135個(gè)基因組尺度代謝網(wǎng)絡(luò)模型，涉及64種原核生物、19種真核生物和4種古生菌 (GSMNDB：http://synbio.tju.edu.cn/GSMNDB/gsmndb.htm)。國(guó)內(nèi)外已經(jīng)有數(shù)篇關(guān)于基因組尺度代謝網(wǎng)絡(luò)模型重構(gòu)以及模擬方法的綜述[8-12]，模型在模擬基因敲除、靶點(diǎn)藥物發(fā)現(xiàn)、擾動(dòng)實(shí)驗(yàn)預(yù)測(cè)以及物種進(jìn)化研究方面都有著重要作用[13]。

雖然最小基因組的研究在濕實(shí)驗(yàn)方面已經(jīng)有了一些進(jìn)展，但是到目前為止還沒(méi)有一個(gè)完整的關(guān)于最小基因組的代謝網(wǎng)絡(luò)模型和計(jì)算方法，Yizhak等[14]注重研究基因敲除與合成生物量之間的關(guān)系，Pa等[15]構(gòu)建的模型只是對(duì)于模型的簡(jiǎn)化，并沒(méi)有得到一個(gè)不能刪除任何基因的最小基因組模型。在本文中我們按照“自上而下”的研究思路，通過(guò)對(duì)基因組尺度代謝網(wǎng)絡(luò)模型的逐步精簡(jiǎn)第一次獲得最小基因組代謝網(wǎng)絡(luò)模型。

1 方法

1.1 通量平衡分析方法

在基因組規(guī)模代謝網(wǎng)絡(luò)模型分析中較為廣泛應(yīng)用的是基于約束的最優(yōu)化方法，這一方法主要由決策變量、約束條件和目標(biāo)函數(shù)三部分組成，本文采用其中的通量平衡分析 (Flux balance analysis，F(xiàn)BA)[16-17]方法進(jìn)行模擬計(jì)算，在FBA中，決策變量就是模型中所有反應(yīng)的通量值，不同的通量分配方案決定著不同的代謝狀態(tài)，限制條件是實(shí)現(xiàn)最優(yōu)化模擬必須滿足的基本條件，包括等式限制和不等式限制，等式限制條件基于擬穩(wěn)態(tài)條件下所有中間代謝物滿足物質(zhì)守恒的假設(shè)，即代謝物的生成與消耗相等，而不等式限制基于模型中每個(gè)反應(yīng)速度的極限值。約束條件表示為數(shù)學(xué)形式如下：

其中，vmin和vmax分別表示反應(yīng)速度的最小和最大值，在本文中，在沒(méi)有特殊限制條件下，對(duì)于可逆反應(yīng)，其最小值和最大值分別為?1 000和1 000 mmol g/(DW·h)，而對(duì)于不可逆反應(yīng)，其最小值和最大值則分別為0和1 000 mmol g/(DW·h)。目標(biāo)函數(shù)在數(shù)學(xué)上可以表示為決策變量的解析函數(shù) (即一系列決策變量的線性或非線性組合)，根據(jù)不同的模擬目的，目標(biāo)函數(shù)既可以只有一個(gè)，也可以同時(shí)有多個(gè)，本文的目標(biāo)是在保持細(xì)胞生長(zhǎng)不受影響的情況下將基因組最小化，因此以生物量的最大速率為目標(biāo)函數(shù)。在本文中當(dāng)一個(gè)基因被敲除后，生物合成量的模擬值小于最初模型的生物合成量，這樣的基因被定義為必需基因。我們利用的計(jì)算工具是Matlab平臺(tái)下運(yùn)行的COBRA[18]工具包。

1.2 通量可變性分析

通量可變性分析 (Flux variability analysis，F(xiàn)VA)[19]是在FBA基礎(chǔ)上發(fā)展的通量分析方法。FBA計(jì)算出來(lái)的最優(yōu)解不是唯一的，有多種通量分布可以使目標(biāo)反應(yīng)達(dá)到最優(yōu)值。FVA就是在保證目標(biāo)反應(yīng)在最優(yōu)值的情況下得到每個(gè)反應(yīng)的通量的可能最大值和最小值，從而得到每個(gè)反應(yīng)的通量變化范圍。我們定義那些通量變化范圍總是大于0或總是小于0的反應(yīng)為達(dá)到最優(yōu)解的必需反應(yīng)，因?yàn)檫@些反應(yīng)對(duì)每個(gè)最優(yōu)解都是必須要存在的。這里要說(shuō)明的是必需反應(yīng)對(duì)應(yīng)的基因并非都是必需基因，因?yàn)橐粋€(gè)反應(yīng)可能有多個(gè)同工酶基因?qū)?yīng)。除了必需反應(yīng)外的其他反應(yīng)定義為非必需反應(yīng)，這些非必需反應(yīng)可以單個(gè)敲除但并不能同時(shí)敲除。

1.3 組合刪除法

在代謝網(wǎng)絡(luò)中從代謝物A到B可能有多個(gè)不同的代謝途徑，如圖1所示，從A到B可以通過(guò)A-X-B、A-Y-Z-B和A-M-N-O-B三條路線實(shí)現(xiàn)，顯然A-X-B途徑最短。因而需要考慮的問(wèn)題是如何保證在精簡(jiǎn)的代謝網(wǎng)絡(luò)中從代謝物A到B需要通過(guò)的途徑是基因X而不是Y-Z或者M(jìn)-N-O。我們采用的方法是將參與從A到B所有途徑的基因{X,Y,Z,M,N,O}進(jìn)行分步組合模擬刪除，首先進(jìn)行單基因敲除實(shí)驗(yàn)，發(fā)現(xiàn)可被敲除的基因數(shù)目即非必需基因數(shù)目為6個(gè)，第一步從這6個(gè)基因里取出6個(gè)基因組合模擬刪除，看能否得到代謝物B，只有一種取法C(6,6)=1，很明顯不能，那我們就從6個(gè)基因里選出5個(gè)基因進(jìn)行組合，共有C(5,6)=6種選法，每種選法都進(jìn)行一次模擬刪除，看能否得到代謝物B，在本例中當(dāng)選取的組合為{Y,Z,M,N,O}的時(shí)候，代謝物A能夠生成B，最后保留的基因數(shù)目就為6(所有非必需基因數(shù)目)?5(選取組合包含基因的數(shù)目)=1(非必需基因保留的數(shù)目)，我們將這種思想運(yùn)用到整個(gè)模型上，以葡萄糖為底物，生物量合成為目標(biāo)函數(shù)，利用程序?qū)Ψ潜匦杌蜻M(jìn)行逐步遞減組合模擬刪除，在保證生物量不變的情況下，保留刪除基因最多的反應(yīng)途徑在模型中，最后得到基因組最小化的代謝網(wǎng)絡(luò)模型。

圖1 從代謝物A到B的代謝途徑示例Fig.1 Metabolic pathways from metabolite A to B.

1.4 基因組簡(jiǎn)化的步驟

整個(gè)基因組簡(jiǎn)化步驟如圖2所示，主要包括以下步驟：

1)對(duì)初始代謝網(wǎng)絡(luò)模型進(jìn)行FVA分析，其中有很多反應(yīng)的最大通量和最小通量都為0，即這些反應(yīng)與細(xì)胞最優(yōu)生長(zhǎng)無(wú)關(guān)，因而首先將這些反應(yīng)從網(wǎng)絡(luò)中刪除。

圖2 基因組最小化的步驟 (圖中數(shù)字為各步處理后的基因數(shù))Fig.2 Procedure of genome minimization.The numbers of genes after every step are shown in the diagram.

2)處理同功酶基因。在模型中有很多同功酶基因催化同一反應(yīng)的情況，這些基因只需存在一個(gè) (組)，反應(yīng)就能夠進(jìn)行，因而多余的基因可以去除。需要注意的是在這些同功能基因中，有些同時(shí)催化多個(gè)反應(yīng)，其中有些反應(yīng)是菌體生長(zhǎng)所必需的，因而簡(jiǎn)單地去除每個(gè)同功能基因可能導(dǎo)致菌體的死亡，因而我們僅將只催化單個(gè) (一類)反應(yīng)的同功能基因去除，而同時(shí)催化多個(gè)反應(yīng)的同功能基因則保留下來(lái)進(jìn)行進(jìn)一步計(jì)算。

3)經(jīng)過(guò)上述兩步我們得到了部分縮減的代謝網(wǎng)絡(luò)模型，通過(guò)FVA分析模型中的反應(yīng)分為必需反應(yīng)和非必需反應(yīng)兩部分，在此基礎(chǔ)上我們將該模型中的基因分為兩個(gè)子集：保留基因集和候選基因集。所有必需反應(yīng)中的基因構(gòu)成保留基因集，要提出的是有的必需反應(yīng)是有幾個(gè)同功能基因催化的，只要保證該必需反應(yīng)能夠進(jìn)行，多余的基因是可以被刪除的，所以保留基因集并不等于所有基因都保留，除了保留基因以外的基因則為候選基因集。

4)保留基因集的基因組最小化，利用單基因敲除函數(shù)找出必需反應(yīng)部分可以敲除的基因，通過(guò)組合刪除方法進(jìn)行篩選，找到必需反應(yīng)部分對(duì)應(yīng)的最小基因組。

5)候選基因集的基因組最小化，在4的基礎(chǔ)上，同樣利用單基因敲除函數(shù)找出可以被敲除的基因，因?yàn)榉潜匦杌驍?shù)目較多，但由于4步驟中對(duì)應(yīng)必需反應(yīng)的非必需基因已經(jīng)在前一步進(jìn)行了單獨(dú)處理，大部分必需基因得到保留，因而在做組合刪除的時(shí)候運(yùn)算量大大減少。最后得到的新模型就是該物種模型的基因組最小化代謝網(wǎng)絡(luò)模型。

2 結(jié)果與分析

2.1 由代謝網(wǎng)絡(luò)分析進(jìn)行大腸桿菌基因組最小化

大腸桿菌是目前了解得最清楚、分析得最透徹的微生物之一，是遺傳學(xué)研究、生化研究和代謝模擬研究的重要平臺(tái)，大腸桿菌作為遺傳重組領(lǐng)域廣泛應(yīng)用的宿主之一，在醫(yī)藥和發(fā)酵工業(yè)中可用于生產(chǎn)重組蛋白、氨基酸和其他化學(xué)品。因?yàn)榇竽c桿菌在動(dòng)物腸道和環(huán)境中進(jìn)化其部分基因組對(duì)某些用途是不需要的，甚至可能是有害的，因此，通過(guò)刪除盡可能多的基因片段，構(gòu)建遺傳穩(wěn)定的菌株，可使大腸桿菌具有健壯的代謝行為，同時(shí)也可以添加具有實(shí)際用途的基因。

目前大腸桿菌的基因組尺度代謝網(wǎng)絡(luò)模型共有 4 個(gè)，包括iJE660[20]、iJR904[21]、iAF1260[22]和iJO1366[23]，iAF1260在必需基因預(yù)測(cè)及與實(shí)驗(yàn)結(jié)果的擬合度上準(zhǔn)確率較好，所以本文采用iAF1260模型作為出發(fā)模型，培養(yǎng)條件設(shè)定為最小培養(yǎng)基。葡萄糖為唯一碳源，其通量設(shè)定為?8 mmol g/(DW·h)。最優(yōu)生物合成速率為0.7367 mmol g/(DW·h)。

通過(guò)前文方法部分介紹的步驟，我們對(duì)基因組尺度代謝網(wǎng)絡(luò)模型分步驟進(jìn)行簡(jiǎn)化，各步驟結(jié)果 (圖2)如下：

1)我們首先刪除模型中的經(jīng)過(guò)FVA分析后的0通量反應(yīng)，在2 382個(gè)反應(yīng)中共有1 535個(gè)反應(yīng)在FVA結(jié)果中的代謝通量為0，最后將0通量反應(yīng)刪除后，得到847個(gè)反應(yīng)，對(duì)應(yīng)731個(gè)基因。我們將由這些反應(yīng)構(gòu)成的新模型稱為model1，將其轉(zhuǎn)化為新的SBML文件 (數(shù)據(jù)未顯示)，以便于進(jìn)行下一步的計(jì)算。

2)優(yōu)化網(wǎng)絡(luò)，去除僅催化一個(gè)或一類反應(yīng)的同功酶基因，這樣的基因共有100個(gè)，因而得到的模型 (model2)包含847個(gè)反應(yīng)，對(duì)應(yīng)631個(gè)基因，將其重新轉(zhuǎn)化為新的SBML文件 (數(shù)據(jù)未顯示)，以便于進(jìn)行下一步計(jì)算。

3)通過(guò)對(duì)FVA結(jié)果的分析，model2模型中必需反應(yīng)部分包含366個(gè)反應(yīng)，對(duì)應(yīng)有305個(gè)基因，利用單基因敲除函數(shù)，找出非必需基因40個(gè) (具體數(shù)據(jù)作為附件，可在網(wǎng)絡(luò)版中下載。下文簡(jiǎn)稱“附件”)，對(duì)這40個(gè)基因進(jìn)行逐步遞減組合模擬刪除，最后確定其中26個(gè)基因 (附件)敲除后對(duì)菌體生長(zhǎng)沒(méi)有影響，這部分基因敲除后整個(gè)模型包含841個(gè)反應(yīng)，605個(gè)基因。此時(shí)必需反應(yīng)部分中無(wú)任何基因可以被敲除，規(guī)模達(dá)到最小。我們將由這些反應(yīng)構(gòu)成的新模型稱為(model3)，將其轉(zhuǎn)化為新的SBML文件 (數(shù)據(jù)未顯示)，以便于進(jìn)行下一步的計(jì)算。

4)將model3的非必需反應(yīng)部分進(jìn)行基因組最小化，同樣利用單基因敲除函數(shù)找出非必需基因，其數(shù)目為304個(gè) (附件)，對(duì)這些非必需基因，我們采用組合法進(jìn)行刪除，在保證合成生物量不變的情況下，依次從里面取出304個(gè)基因、303個(gè)基因、302個(gè)基因……，以此類推，進(jìn)行模擬刪除，當(dāng)取出293個(gè)基因 (附件)刪除的時(shí)候得到一種組合方式，滿足條件，所以非必需基因保留的最小值就是304?293=11(個(gè))，最終得到含有403個(gè)反應(yīng)，對(duì)應(yīng)312個(gè)基因的模型model4，將該模型轉(zhuǎn)化為SBML文件 (數(shù)據(jù)未顯示)。

5)將新模型model4再次進(jìn)行單基因敲除實(shí)驗(yàn)，發(fā)現(xiàn)可以敲除的基因數(shù)目為0，這驗(yàn)證了利用這一方法得到的確實(shí)是最小基因組。利用該基因組最小化的代謝網(wǎng)絡(luò)模型，在葡萄糖輸入設(shè)定為?8 mmol g/(DW·h)時(shí)求得最優(yōu)生物質(zhì)合成速率為0.7367 mmol g/(DW·h)，與原始模型相同，表明基因組精簡(jiǎn)后，菌體生長(zhǎng)不受影響。整個(gè)模型簡(jiǎn)化過(guò)程中各個(gè)步驟的結(jié)果如表1、圖2、圖3及附件所示。

從表1、圖3可以看出在簡(jiǎn)化基因組的過(guò)程中，F(xiàn)VA分析和非必需反應(yīng)部分刪除的基因和反應(yīng)占所有刪除的基因和反應(yīng)數(shù)目的大部分，這符合生物學(xué)的特征，最后整個(gè)基因組縮小為原來(lái)的24.7%，這也為濕實(shí)驗(yàn)尋求敲除目標(biāo)基因部分提供參考。

2.2 最小基因組代謝網(wǎng)絡(luò)模型途徑分析

為了進(jìn)一步了解得到的基因組最小化代謝網(wǎng)絡(luò)模型的特征，我們對(duì)模型中反應(yīng)的代謝途徑分布進(jìn)行了分析。模型中的403個(gè)反應(yīng)分布在32個(gè)代謝途徑中，具體分布情況如表2所示。

表1 基于代謝網(wǎng)絡(luò)模型的大腸桿菌基因組最小化結(jié)果Table 1 Minimization of E.coli genome based on metabolic network analysis

圖3 每個(gè)步驟后的基因數(shù)目和反應(yīng)數(shù)與初始模型的比較結(jié)果Fig.3 Comparison of the gene number and the reactions number with the initial model after every step.

從表2可以看出在基因組最小模型中，氨基酸代謝和脂類代謝包含的反應(yīng)數(shù)最多，分別為112個(gè)和80個(gè)，這表明大腸桿菌的生長(zhǎng)過(guò)程中，氨基酸、脂類代謝具有最為重要的作用，這與氨基酸和脂類物質(zhì)在大腸桿菌細(xì)胞組成中占有較高比例相符。另外，輔因子代謝包含77個(gè)反應(yīng)，在整個(gè)網(wǎng)絡(luò)中也具有很大的比重，這表明輔因子是合成生物量不可缺少的重要成分。另外要注意的一點(diǎn)是交換反應(yīng)在模型同樣具有重要的作用，任何基因組最小化的模型都是相對(duì)的，當(dāng)培養(yǎng)基條件發(fā)生改變時(shí)，模型包含的交換也會(huì)發(fā)生改變，交換反應(yīng)越多時(shí)，表明外界提供的營(yíng)養(yǎng)物質(zhì)越多，則菌體自身合成的物質(zhì)越少，從而模型所包含的反應(yīng)數(shù)和基因數(shù)都會(huì)相應(yīng)的減少。

另外，從碳源的利用上來(lái)看，在基因組最小化模型模擬中乳酸的合成途徑完全被刪除，其他副產(chǎn)物乙酸、琥珀酸合成途徑通量也為零，表明底物充分用于生物量的合成，去除了分流到其他副產(chǎn)物而產(chǎn)生的碳源消耗，該結(jié)果顯示了基因組最小化代謝網(wǎng)絡(luò)模型對(duì)生物代謝途徑的優(yōu)化。

2.3 與濕實(shí)驗(yàn)結(jié)果的比較

濕實(shí)驗(yàn)方面，Posfai等[5]通過(guò)大片段基因敲除成功將E.coliK-12 MG1655基因組減少了15%，構(gòu)建了 MDS40、41、42和 43菌株，Kolisnychenko等將大腸桿菌[24]基因組減少了 8.1%，構(gòu)建了MDS12菌株，Hashimoto等[25]將大腸桿菌基因組減少了29.7%構(gòu)建了△16菌株，但其類核組織表現(xiàn)異常，生長(zhǎng)速率較低，在M9培養(yǎng)基上無(wú)法生長(zhǎng)。同時(shí)Mizoguchi等[26]通過(guò)無(wú)痕敲除方法刪除了103個(gè)候選區(qū)域，共刪除1 080個(gè)基因，有84個(gè)刪除突變株的生長(zhǎng)行為與野生株類似，將這些刪除累積起來(lái)構(gòu)建了一個(gè)基因組減小菌株——MGF-01。Hirokawa等[27]在MGF-01的基礎(chǔ)上繼續(xù)刪減得到結(jié)構(gòu)穩(wěn)定，沒(méi)有任何營(yíng)養(yǎng)缺陷性的，可以更高效生長(zhǎng)的菌株——DGF-327和DGF-298。我們綜合考慮選用MDS43和DGF-298菌株與模型相比較，因?yàn)檫@兩種菌株的生產(chǎn)速率與野生型類似，且穩(wěn)定性強(qiáng)，電轉(zhuǎn)化效率高。

表2 大腸桿菌基因組最小化代謝網(wǎng)絡(luò)模型的代謝途徑分布Table 2 Metabolism pathways of the minimized metabolic network in E.coli

Posfai等[5]通過(guò)構(gòu)建的基因組縮減菌株MDS43共敲除了743個(gè)基因，iAF1260模型中包含其中的82個(gè)基因，這82個(gè)基因在基因組最小化模型并沒(méi)有被全部去除 (附件)，其中5個(gè)基因仍然保留在模型中，這5個(gè)基因是b0339(cynt)、b2976(glcb)、b2979(glcd)、b4467(glcf) 和 b4468(glce)。b0339、b2976兩個(gè)基因在原有模型當(dāng)中都具有同功酶，即處于“or”的關(guān)系，由于組合刪除的隨機(jī)性，導(dǎo)致刪除了與其同功能的基因，而保留了這兩個(gè)基因，這一點(diǎn)對(duì)于菌體生長(zhǎng)并沒(méi)有影響。而剩下的3個(gè)基因 (b2979、b4467和b4468)構(gòu)成一酶復(fù)合物催化反應(yīng)GLYCTO2(Glycolate+Q8=Glyoxylate+Q8H2)。該反應(yīng)可以將代謝網(wǎng)絡(luò)中某些反應(yīng)過(guò)程生成的乙醇酸轉(zhuǎn)化為乙醛酸，進(jìn)而通過(guò)乙醛酸循環(huán)轉(zhuǎn)化為生物質(zhì)合成需要的前體。我們注意到大腸桿菌代謝網(wǎng)絡(luò)模型中還有另外兩個(gè)基因 (b1033、b3553)催化乙醇酸和乙醛酸間的轉(zhuǎn)化但所用輔酶為NADH或NADPH(R_GLYCLTD：Glycolate+NAD(P)<==>Glyoxylate+NAD(P)H2)，但該反應(yīng)被定義為不可逆且反應(yīng)方向?yàn)閺囊胰┧岬揭掖妓帷．?dāng)我們將該反應(yīng)改為可逆后進(jìn)行基因組最小化計(jì)算得到的最小基因組中就不包括b2979、b4467和b4468而是只包括b1033或b3553，基因個(gè)數(shù)也從312個(gè)減少到了310個(gè)。這表明大腸桿菌中很可能反應(yīng)R_GLYCLTD也是可逆的，因此敲除b2979、b4467和b4468這3個(gè)基因?qū)ιL(zhǎng)沒(méi)有影響。

我們還將得到的最小基因組與Hirokawa等[27]構(gòu)建的 DGM-298菌株進(jìn)行了比較，DGF-298一共刪除了1 429個(gè)基因，其中iAF1260模型中包含被實(shí)驗(yàn)證明可被刪除的基因數(shù)為298個(gè) (附件)，同樣這298個(gè)基因在基因組最小化模型中并沒(méi)有被全部去除，其中有12個(gè)基因仍保留在最后的模型當(dāng)中，分別是b0149(mrcb)、b0339(cynt)、b0381(ddla)、b0521(ybcf)、b1415(alda)、b2097(fbab)、b2708(gutq)、b2797(sdab)、b2979(glcd)、b3744(asna)、b4094(phnn)、b4467(glcf) 和 b4468(glce)(附件)。其中 b0149、b0339、b0381、b4383、b0521、b2097、b2708 和 b2797七個(gè)基因是因?yàn)橥δ芑?，即“or”關(guān)系的選擇而被保留在模型中。b2979、b4467和b4468三個(gè)基因在上一部分已經(jīng)說(shuō)明是因?yàn)槟Ｐ椭幸恍┓磻?yīng)的可逆性設(shè)置的問(wèn)題。b1415和b3744兩個(gè)基因所催化的反應(yīng)分別是gcald+h2o+nad<==>glyclt+nadh和asp-L+atp+nh4<==>amp+asn-L+ppi，這兩個(gè)反應(yīng)在模型中以葡萄糖為唯一碳源時(shí)為必需反應(yīng)。而Hirokawa等采用的是富培養(yǎng)基，可以從培養(yǎng)基中直接得到天門冬酰胺等而不需要從葡萄糖生成，因此在其濕實(shí)驗(yàn)中可以敲除而在我們的模擬計(jì)算中保留在最小基因組中。

從上面的結(jié)果可以看出，當(dāng)需要敲除某些基因的時(shí)候，可以利用得到的最小基因組模型反饋于濕實(shí)驗(yàn)，確定哪些基因可以嘗試敲除，哪些途徑可以得到保留，模擬的結(jié)果為濕實(shí)驗(yàn)的設(shè)計(jì)提供了有價(jià)值的參考。

3 結(jié)論與展望

現(xiàn)代生物技術(shù)的高速發(fā)展，特別是組學(xué)技術(shù)的興起，高通量數(shù)據(jù)的迅速積累，促進(jìn)了生物代謝網(wǎng)絡(luò)模型的快速發(fā)展，使網(wǎng)絡(luò)構(gòu)建的理論逐漸完善，質(zhì)量大幅度提升，應(yīng)用范圍同樣也顯著擴(kuò)大，在此基礎(chǔ)上，本文第一次提出了基于基因組尺度代謝網(wǎng)絡(luò)模型，通過(guò)模型簡(jiǎn)化構(gòu)建基因組最小化代謝網(wǎng)絡(luò)模型的方法，該方法通過(guò)零通量反應(yīng)刪除和對(duì)非必需基因的組合刪除模擬計(jì)算，確定能夠計(jì)算的基因組最小化代謝網(wǎng)絡(luò)模型。該方法依托于超級(jí)計(jì)算機(jī)的計(jì)算能力，運(yùn)用Matlab中的COBRA工具包，以大腸桿菌的iAF1260模型為例，將模型從最初的1 260個(gè)基因縮減到312個(gè)基因，模型簡(jiǎn)化后生物量合成速率保持不變。通過(guò)該方法得到的基因組最小化代謝網(wǎng)絡(luò)模型中氨基酸代謝及脂類代謝仍為生物體生長(zhǎng)的主要需求，符合生物代謝規(guī)律，證明了基于非必需基因組合的刪除方法的合理性。

通過(guò)對(duì)基因組尺度代謝網(wǎng)絡(luò)模型的精簡(jiǎn)，我們可以確定對(duì)于生物生長(zhǎng)最為重要的反應(yīng)和基因，并且得到相應(yīng)代謝途徑、代謝系統(tǒng)的生物信息，通過(guò)計(jì)算機(jī)模擬得到的結(jié)果可以為濕實(shí)驗(yàn)的方案設(shè)計(jì)提供重要的參考。由于該方法是在基因組尺度代謝網(wǎng)絡(luò)模型的基礎(chǔ)上實(shí)施的，因而高質(zhì)量的模型對(duì)于最終得到的基因組最小化代謝網(wǎng)絡(luò)模型非常重要，基因組最小化代謝網(wǎng)絡(luò)模型對(duì)于今后專一生產(chǎn)某種代謝物的途徑設(shè)計(jì)也將提供更好思路。另外，隨著網(wǎng)絡(luò)模型的發(fā)展，基因組最小化代謝網(wǎng)絡(luò)模型中也可能會(huì)結(jié)合調(diào)控、信號(hào)轉(zhuǎn)導(dǎo)等信息，使基因組最小化的生物網(wǎng)絡(luò)向更加趨近于最小生命體的方向發(fā)展。

[1]Gibson DG,Glass JI,Lartigue C,et al.Creation of a bacterialcell controlled by a chemically synthesized genome.Science,2010,329(5987):52?56.

[2]Fraser CM,Gocayne JD,White O,et al.The minimal gene complement ofMycoplasma genitalium.Science,1995,270(5235):397?403.

[3]Lartigue C,Glass JI,Alperovich N,et al.Genome transplantation in bacteria:changing one species to another.Science,2007,317(5838):632?638.

[4]Gibson DG,Benders GA,Andrews-Pfannkoch C,et al.Complete chemical synthesis,assembly,and cloning ofaMycoplasmagenitaliumgenome.Science,2008,319(5867):1215?1220.

[5]Posfai G,Plunkett G,3rd,Feher T,et al.Emergent properties of reduced-genomeEscherichia coli.Science,2006,312(5776):1044?1046.

[6]Juhas M,Eberl L,Glass JI.Essence of life:essential genes of minimal genomes.Trends Cell Biol,2011,21(10):562?568.

[7]Qiu DR.Essential genes,minimal genome and synthetic cell of bacteria:a review.Chin J Biotech,2012,28(5):540?549(in Chinese).

邱東茹.細(xì)菌必需基因、最小基因組和合成細(xì)胞.生物工程學(xué)報(bào),2012,28(5):540?549.

[8]Durot M, Bourguignon PY, Schachter V.Genome-scale modelsofbacterialmetabolism:reconstruction and applications.FEMS Microbiol Rev,2009,33(1):164?190.

[9]Thiele I,Palsson BO.A protocol for generating a high-quality genome-scale metabolic reconstruction.Nat Protoc,2010,5(1):93?121.

[10]Feist AM, Herrgard MJ, Thiele I, et al.Reconstruction of biochemical networks in microorganisms.Nat Rev Microbiol,2009,7(2):129?143.

[11]Park JM,Kim TY,Lee SY.Constraints-based genome-scale metabolic simulation for systems metabolic engineering.BiotechnolAdv,2009,27(6):979?988.

[12]Hao T,Ma HW,Zhao XM.Progress in the automatic reconstruction of genome-scale metabolic network.Chin J Biotech,2012,28(6):661?670(in Chinese).

郝彤,馬紅武,趙學(xué)明.基因組尺度代謝網(wǎng)絡(luò)自動(dòng)重構(gòu)及分析工具研究進(jìn)展.生物工程學(xué)報(bào),2012,28(6):661?670.

[13]Wang H,Ma HW,Zhao XM.Progressin genome-scale metabolic network:a review.Chin J Biotech,2010,26(10):1340?1348(in Chinese).

王暉,馬紅武,趙學(xué)明.基因組尺度代謝網(wǎng)絡(luò)研究進(jìn)展.生物工程學(xué)報(bào),2010,26(10):1340?1348.

[14]Yizhak K,Tuller T,Papp B,et al.Metabolic modeling of endosymbiont genome reduction on a temporal scale.Mol Syst Biol,2011,7:479.

[15]Pal C,Papp B,Lercher MJ,et al.Chance and necessity in the evolution of minimal metabolic networks.Nature,2006,440(7084):667?670.

[16]Orth JD,Thiele I,Palsson BO.What is flux balance analysis?Nat Biotechnol,2010,28(3):245?248.

[17]Edwards JS,Palsson BO.Metabolic flux balance analysis and the in silico analysis of Escherichia coli K-12 gene deletions.BMC Bioinformatics,2000,1:1.

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

[19]Muller AC,Bockmayr A.Fast thermodynamically constrained flux variability analysis.Bioinformatics,2013,29(7):903?909.

[20]Edwards JS,Palsson BO.TheEscherichia coliMG1655 in silico metabolic genotype: its definition,characteristics,and capabilities.Proc Natl Acad Sci USA,2000,97(10):5528?5533.

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

[22]FeistAM,Henry CS,Reed JL,etal.A genome-scale metabolic reconstruction forEscherichia coliK-12 MG1655 that accounts for 1260 ORFs and thermodynamic information.Mol Systems Biol,2007,3(121):1?18.

[23]Orth JD,Conrad TM,Na J,et al.A comprehensive genome-scale reconstruction ofEscherichia colimetabolism—2011. Mol Systems Biol, 2011,7(535):1?9.

[24]Kolisnychenko V,Plunkett G 3rd,Herring CD,et al. Engineering a reducedEscherichiacoligenome.Genome Res,2002,12(4):640-647.

[25]Hashimoto M,Ichimura T,Mizoguchi H,et al.Cell size and nucleoid organization ofengineeredEscherichia colicells with a reduced genome.Mol Microbiol,2005,55(1):137?149.

[26]Mizoguchi H,Mori H,Fujio T.Escherichia coliminimum genome factory. Biotechnol Appl Biochem,2007,46(Pt 3):157?167.

[27]Hirokawa Y,Kawano H,Tanaka-Masuda K,et al.Genetic manipulations restored the growth fitness ofreduced-genomeEscherichiacoli.JBiosci Bioeng,2013,1(10):1?7.