芮 斌，金杰軍，范亞東，文 漢

（安徽農業(yè)大學生命科學學院，合肥230036）

大腸桿菌E.coli是進行微生物遺傳研究的重要材料，如局限性轉導就是1954年在大腸桿菌K12菌株中發(fā)現的.萊德伯格（Lederberg）采用兩株大腸桿菌的營養(yǎng)缺陷型進行實驗，奠定了研究細菌接合方法學以及基因工程研究的基礎.大腸桿菌作為外源基因表達的宿主，遺傳背景清楚，技術操作簡單，培養(yǎng)條件簡單，大規(guī)模發(fā)酵經濟，所以倍受遺傳工程專家的重視.目前大腸桿菌是應用最廣泛，最成功的表達體系，常做高效表達的首選體系.

系統(tǒng)生物學著重以系統(tǒng)的思路和方法定量研究生物過程.由于生物系統(tǒng)的組成關系主要用網絡來表示，代謝網絡研究就自然成為熱門領域［1－2］.最初的網絡理論往往只涉及一個底物和一個產物，故而與實際網絡不符.FELL等學者致力于線性代數的方法來更便利地計算代謝途徑，其主要思想是先假設代謝網絡的途徑是代謝網絡計量系數矩陣的線性無關解，然后通過對這些途徑的線性組合或者枚舉來產生存在的（但未必是熱力學可行的）代謝途徑.最近，凸分析被運用到代謝途徑分析中，代謝網絡理論于是得到了更進一步的發(fā)展，漸漸形成了普遍認可的兩種方法：基元模式分析（elementary flux mode analysis）和極端途徑分析（extreme pathway analysis）.同時在另一個層面上，對代謝網絡的研究進一步逐漸演化出單獨的靜態(tài)分析和動態(tài)分析.靜態(tài)分析主要包括Schilling等學者提出的極端途徑（EPS），基元模式（EFMS），FBA流量平衡分析和jeong等學者提出的圖論法（graph theory）［3－6］.而動態(tài)分析方法主要是基于動力學信息的，它包括代謝控制分析（metabolic control analysis），控制論模擬方法（cybernetic modeling）.目前對基因組規(guī)模的代謝網絡來說，許多反應的動力學參數和反應物濃度參數都是不清楚的，故靜態(tài)的研究方法占據了目前代謝網絡研究的主流位置［7－8］.

本文在matlab cobra工具箱的基礎上對大腸桿菌代謝系統(tǒng)進行FBA和FVA分析［9－12］，將分析結果進行匯總后，為顯示在敲除和不敲除各種酶基因的結果變化，對數據進行繪圖，通過直觀的圖形對比，顯示出敲除各種酶基因對大腸桿菌代謝系統(tǒng)中每一個反應的影響，以及在敲除大腸桿菌某個酶基因后代謝系統(tǒng)的全部反應的FBA和FVA的變化之間是否存在共性，以此為依據，大致判斷各個酶在大腸桿菌代謝系統(tǒng)中所起的作用.大腸桿菌代謝系統(tǒng)中有很多的酶基因參與，本研究以25個酶為例，分別單獨將每一個酶相應的酶基因從代謝系統(tǒng)中敲除后，對代謝系統(tǒng)進行FBA和FVA分析［13－15］，將結果與未敲除任何酶基因的代謝系統(tǒng)的FVA和FBA結果進行對比，以確定所敲除的酶在代謝系統(tǒng)中所起的作用，這提供了一個全新的酶分析方法.這樣的結果雖說只是猜測，但可對以后的實驗起到指導性作用.

1 材料與方法

1.1 所需操作軟件以及文件

Matlab版本7.0以上.cobra toolbox版本2.0以上.libsbml版本4.0.1以上.Sbml toolbox版本3.1.1以上.gurobi 5.0以上.

大腸桿菌小型中心代謝模型‘ecoli＿core＿model.xml’，此模型中包含大腸桿菌中心代謝的全部信息［6］.

1.2 方法

1.2.1 代謝模型讀取 在matlab命令窗口中執(zhí)行指令：model＝readCb Model（‘文件名’）其中文件名為模型保存目錄［11，15］.

1.2.2 模型的修飾 執(zhí)行指令：

model＝changeRxnBounds（model，＇EX＿glc（e）＇，－16.7，＇l＇）；其中＇EX＿glc（e）＇代表了底物攝取速率，‘l’是‘lower’的縮寫，表示模型的下限［2］.1.2.3酶基因的敲除 執(zhí)行指令：

model＝removeRxns（model，＇rxnRemove List＇）；＇rxn Remove List＇代表的是含有所敲除反應的反應信息的單元矩陣，酶基因的敲除中，rxnRemove List可直接用所需敲除酶的英文縮寫代替.如敲除烯醇化酶即ENO酶，直接在命令窗口中輸入指令model＝removeRxns（model，＇ENO＇）即可，執(zhí)行命令后表示該酶所直接參與的反應的電子信息從模型中去除，該反應在模型中不再起作用.

1.2.4 FBA最優(yōu)化值的獲取 執(zhí)行指令：change－CobraSolver（＇gurobi5＇，＇LP＇）；Fluxdistribution＝optimizeCb Model（model）；flux＝geometricFBA（model），

其中第一步指令是修改cobra solver；第二步是執(zhí)行標準FBA命令.

1.2.5 FVA（Flux Variability Analysis）值的獲取 執(zhí)行命令：

［minFlux，max Flux］＝flux Variability（model，90）；其中‘90’代表一個百分比，為模型所能達到的目標函數相對于最優(yōu)值的比例.FVA命令的結果有min Flux和max Flux，分別代表模型中95個反應的通量分布情況.

1.3 作圖說明

作圖從兩個方向進行，分別為橫向和縱向，橫向作圖在于顯示模型中每個單獨反應在不同酶基因被敲除后受到的影響，作圖時橫坐標代表被敲除的酶基因，為更好的顯示結果，分別用數字代表，對應情況為：1.沒有敲除任何酶基因.2.EDD.3.GLUUPT.4.GDH.5.GND.6.TK2.7.PFK.8.HX I.9.TAL.10.TK1.11.SER.12.SHMT.13.GEV.14.ENO.15.PYK.16.PDH.17.ACEX.18.GLTA.19.ICD.20.AKD.21.FUM.22.MDH.23.GS1.24.GS2.25.PPC.26.MEZ［4］.橫向作圖時FBA的結果使用折線圖顯示，FVA的結果選用柱形圖來顯示區(qū)間的變化.縱向作圖在于顯示敲除某個酶基因對整個代謝系統(tǒng)的影響，作圖時橫坐標的數字分別代表大腸桿菌代謝模型中的95個反應.縱向作圖時由于反應有95個，在進行FVA分析時，用柱形圖難以顯示結果的變化，所以將結果進行拆分，分別對minflux和maxflux作圖.

2 結果與分析

2.1 橫向分析

文中只將反應ATP synthase，Biomass，glutamate transport的流量變化圖附上（圖1～圖14），其他反應可以同理進行分析.

2.2 縱向分析

橫向分析后，會發(fā)現在酶基因14（ENO），15（PYK），16（PDH），21（FUM），25（PPC）等被敲除后，反應受到的影響最大，本文只將酶ENO即烯醇化酶，酶PDH即丙酮酸脫氫酶，酶FUM即延胡索酸酶，酶PPC即磷酸烯醇式丙酮酸羧化酶的縱向分析圖附上，其它酶可以用同種方法進行分析.

圖1 反應ATP synthase在敲除不同的酶基因后最優(yōu)流量值的變化Fig.1 Changes in reaction ATP synthase after the knockout of different enzyme gene

圖2 反應ATP synthase在敲除不同的酶基因后流量范圍的變化Fig.2 Changes in flux range of reaction ATP synthase after the knock out of different enzyme gene

圖3 反應Biomass在敲除不同的酶基因后最優(yōu)流量值的變化Fig.3 Changes in reaction Biomass after the knockout of different enzyme gene

圖4 反應Biomass在敲除不同的酶基因后流量范圍的變化Fig.4 Changes in the flux range of reaction Biomass after the knockout of different enzyme gene

圖5 反應glutamate transport在敲除不同的酶基因后最優(yōu)流量值的變化Fig.5 Changes in reaction glutamate transport after the knockout of different enzyme gene

圖6 反應glutamate transport在敲除不同的酶基因后流量范圍的變化Fig.6 Changes in flux range of reaction glutamate transport after the knockout of different enzyme gene

圖7 敲除酶基因ENO的代謝系統(tǒng)中全部反應的FBA結果變化Fig.7 The FBA results change of all reaction in ENO knockout metabolic system

圖8 敲除酶ENO的代謝系統(tǒng)中全部反應FVA結果的變化Fig.8 The FVA results change of all reaction in ENO knockout metabolic system

圖9 敲除酶PDH的代謝系統(tǒng)中全部反應的FBA結果變化Fig.9 The FBA results change of all reaction in PDH knockout metabolic system

圖10 敲除酶PDH的代謝系統(tǒng)中全部反應的FVA結果變化Fig.10 The FVA results change of all reaction in ENO knockout metabolic system

圖11 敲除酶FUM的代謝系統(tǒng)中全部反應的FBA結果變化Fig.11 The FBA results change of all reaction in FUM knockout metabolic system

圖12 敲除酶FUM的代謝系統(tǒng)中全部反應的FVA結果變化Fig.12 The FVA results change of all reaction in FUM knockout metabolic system

圖13 敲除酶PPC的代謝系統(tǒng)中全部反應的FBA結果變化Fig.13 The FBA results change of all reaction in PPC knockout metabolic system

3 結論與討論

綜合橫向分析和縱向分析，可以確定酶ENO即烯醇化酶，酶FUM即延胡索酸酶，酶PDH即丙酮酸脫氫酶等對整個代謝系統(tǒng)的影響較大，所以以酶ENO、FUM和PDH為例進行討論.

3.1 關于ENO

圖14 敲除酶PPC的代謝系統(tǒng)中全部反應的FVA結果變化Fig.14 The FVA results change of all reaction in PPC knockout metabolic system

通過縱向分析，發(fā)現ENO敲除后，在大腸桿菌代謝系統(tǒng)中有相當多的反應都受到了影響，而酶的催化作用具有專一性，即只催化一種或一類底物，因此我們可以推測出酶ENO所參與的反應應該是位于大腸桿菌代謝系統(tǒng)比較靠前的位置，可以為代謝系統(tǒng)之后的反應提供底物等.再返回匯總數據，可以發(fā)現，酶ENO被敲除后反應代謝系統(tǒng)中反應就會少一個，說明在這些酶所直接參與的反應中，它們的作用是不可替代的，即大腸桿菌代謝系統(tǒng)中ENO酶無同工酶，由此可以基本判定酶ENO所催化的反應.再觀察橫向分析圖與匯總的數據，會發(fā)現酶ENO對反應的影響可以分為兩類：

（1）ENO酶敲除后，反應被抑制，可分為兩類：A.ENO酶敲出后該反應被抑制，但是并未停止，如反應ATP synthease；B.ENO酶被敲出后該反應完全停止，如反應Biomass.

這說明代謝系統(tǒng)中反應ATP synthease和反應Biomass在酶ENO直接參與的反應后面進行.而且反應Biomass離ENO酶直接參與的反應可能更近一點，這樣有利于我們對大腸桿菌代謝網絡結構的理解.

（2）ENO酶被敲除后，反應被促進，如反應glutamate transport.考慮到大腸桿菌是兼性厭氧性微生物，因此猜測：酶ENO可抑制反應glutamate transport的進行或者是抑制反應glutamate transport底物的合成，也可能是酶ENO的產物對反應有抑制作用，即反應glutamate transport可能是大腸桿菌在另一種不同的環(huán)境才能進行的.

3.2 關于FUM

（1）FUM酶敲除后，會有小部分反應被抑制.可分為兩類：A.FUM酶敲除后該反應被抑制，但是并未停止，如enolase催化的反應；B.FUM酶被敲出后該反應幾乎停止，如由aconitase催化citrate生成isocitrate的反應.這說明代謝系統(tǒng)中FUM反應生成物對enolase和aconitase催化的反應產生了約束，當FUM被敲除后，這二者都受到影響.

（2）在FUM被敲除后，還有一部分反應是發(fā)生了大程度的提升.比如pyruvate reversible transport和Fumarate transport發(fā)生明顯的上調.這個結果是由于fumarate的消除途徑被敲除后，fumarate累積之后外流到胞外.同時，整個TCA cycle的循環(huán)受到阻遏，使得pyruvate也發(fā)生累積并外流到胞外.這代表細菌代謝系統(tǒng)發(fā)生了非常明顯的轉化.

（3）還有部分反應在FUM敲除后，只發(fā)生了非常微弱的上調.這些反應是在野生型菌株中就有著很高流量的反應.這說明FUM的敲除對于代謝網絡的物質分配影響不大.

3.3 關于PDH

PDH敲除后對整體代謝有一定的影響，但是其影響程度并沒有前兩種酶基因敲除后的影響大.這表明，在氧氣充分的條件下，PDH通路在大腸桿菌體內有其替代途徑，存在冗余性.

同理可以應用于大腸桿菌代謝系統(tǒng)中的其他酶和反應的分析，這將有利于我們對大腸桿菌代謝系統(tǒng)的認識.在此所用的大腸桿菌小型中心代謝系統(tǒng)模型，是比較簡單的的代謝系統(tǒng)，這樣的分析方案不僅適用于這樣的小型代謝系統(tǒng)，也可用于對比較復雜的代謝網絡分析.

［1］ 花 強，楊 琛.代謝流量比率分析及其在代謝工程中的應用［J］.生物工程學報，2009，25（9）：1303－1311.

［2］ 史晉輝.細胞代謝過程分析及模型優(yōu)化［J］.生命的化學，2001，21（1）：71－73.

［3］ Tie Shen，Bin Rui，Hong Zhou，et al.Metabolic flux ratio analysis and multi－objective optimization revealed a globally conserved and coordinated metabolic response of E.coli to paraquat－induced oxidative stress［J］.Molecular BioSystems，2012，9（1）：121－132.

［4］ 楊志杰.基于代謝流量分析的微生物細胞代謝特性研究［D］.上海：華東理工大學，2011.

［5］ 陳 琦，王 卓，魏冬青.代謝網絡流分析進展及應用［J］.中國科學，2010，55（14）：1302－1309

［6］ 丁德武，何小青，陸克中，等.用約束法研究大腸桿菌代謝網絡核心模型［J］.計算機與應用化學，2010，27（6）：835－838.

［7］ 李 毅，鄭浩然，鈕俊清，等.VFA：一種可視化代謝網絡建模工具［J］.北京生物醫(yī)學工程，2008，27（5）：490－494.

［8］ 譚青喬，陳 寧，王東洋.代謝工程發(fā)展及應用［J］.生物技術通訊，2001，12（2）：123－126.

［9］ 魏 春，陳 寧.代謝網絡定量分析研究進展［J］.生物技術通訊，2002，13（3）：234－237.

［10］ 馬紅武，趙學明，遲萬忠.應用Excel處理生化過程數據（Ⅰ）代謝通量分析及代謝網絡優(yōu)化［J］.計算機與應用化學，1998，15（6）：274－278.

［11］ 蔣 達.代謝網絡通量和結構分析［M］.大連：大連理工大學出版，2006.

［12］ 趙 靜.細胞代謝網絡的結構、功能與進化研究［D］.上海：上海交通大學，2008.

［13］ Tassoula E A.Growth possibilities of E.Coli in natural waters［J］.International Journal of Environmental Studies，1997，52（1－4）：67－73.

［14］ Kostromins Andrejs，Stalidzans Egils，Paint Net.COBRA Toolbox extension for visualization of stoichiometric models of metabolism［J］.BioSystems，2012，109（2）：233－241.

［15］ Jan Schellenberger，Richard Que，Ronan M T Fleming，et al.Quantitative prediction of cellular metabolism with constraint－based models：the COBRA Toolbox v2.0［J］.Nature，2011，21（6）：1290－1307.