基于約束建模法的結(jié)核菌H37Rv代謝網(wǎng)絡(luò)分析

杜志成，關(guān) 鵬，黃德生

(1.中國醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)教研室，沈陽110001;2.中國醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院數(shù)學(xué)教研室，沈陽110001)

世界衛(wèi)生組織的發(fā)布的《2012年全球結(jié)核病報(bào)告》顯示，在2011年里有870萬新發(fā)病例，估計(jì)有140萬人死于結(jié)核病，其中包括50萬名婦女，該病已經(jīng)成為全世界婦女的主要?dú)⑹种唬?］。多重藥物耐受結(jié)核病和廣泛藥物耐受結(jié)核病的出現(xiàn)［2］提示新的抗結(jié)核病藥物靶標(biāo)相關(guān)研究迫在眉睫。

隨著測序技術(shù)和高通量生物技術(shù)的迅速發(fā)展，大量的生物體已經(jīng)完成全基因組測序、基因功能注釋。研究人員結(jié)合生化、生理數(shù)據(jù)，系統(tǒng)、合理、有效地重建了多種生物體的基因組尺度代謝網(wǎng)絡(luò)(Genome-scale metabolic network，GSMN)(包括金黃色葡萄球菌iSB619、結(jié)核分枝桿菌iNJ661、巴克紅曲霉iAF692和綠膿桿菌iJN746等)，從而得以在計(jì)算機(jī)上進(jìn)行相關(guān)的模擬預(yù)測以供實(shí)驗(yàn)室和臨床研究參考［3-4］。本研究采用基于約束的建模方法對結(jié)核分枝桿菌的基因組尺度代謝網(wǎng)絡(luò)iNJ661［5］進(jìn)行重構(gòu)和仿真分析以進(jìn)一步探索該代謝網(wǎng)絡(luò)的屬性。

1 材料與方法

1.1 數(shù)據(jù)獲取

結(jié)核分枝桿菌國際標(biāo)準(zhǔn)強(qiáng)毒株H37Rv菌株的GSMN模型iNJ661于2013年4月18日從BIGG(http://bigg.ucsd.edu/)下載，格式為 Systems Biology Markup Language(SBML)［6］。它包含了 661個(gè)基因、540個(gè)基因編碼蛋白、826個(gè)代謝物和1 025個(gè)代謝反應(yīng)(其中939個(gè)發(fā)生在細(xì)胞內(nèi)的)。同時(shí)從 KEGG(http://www.genome.jp/kegg)、PubMed(www.ncbi.nlm.nih.gov/pubmed)、TubercuList(http://genolist.pasteur.fr/TubercuList)獲得詳細(xì)的關(guān)于結(jié)核桿菌H37Rv菌株的反應(yīng)、蛋白、基因組等多源化信息。

1.2 COBRA 軟件

約束的重建和分析軟件 COBRA［7］(Constraintbased reconstruction and analysis)是在美國加利福尼亞大學(xué)圣地亞哥分校Palsson和Herrgard領(lǐng)導(dǎo)下開發(fā)的，它是一個(gè)運(yùn)行在 Matlab(Mathworks Inc.，Natick，MA，USA)環(huán)境下專門用于代謝網(wǎng)絡(luò)基于約束的模擬分析工具，主要功能包括:通量平衡分析(Flux balance analysis，F(xiàn)BA)、通量平衡幾何分析、必要基因分析法、最小代謝調(diào)整分析、碳十三追蹤、差異填充(Gap Filling)、代謝工程和計(jì)算機(jī)模擬可視化。該軟件的最新版本可以從Palsson小組的主頁免費(fèi)下載(http://systemsbiology.ucsd.edu/)。

1.3 通量平衡分析

FBA［8-9］是最基本的一種約束建模分析方法，可以計(jì)算細(xì)胞的生長情況。它可以使用化學(xué)計(jì)量矩陣(S)進(jìn)行穩(wěn)態(tài)分析(Sv=0)處理一些系統(tǒng)性問題，其中S是代謝網(wǎng)絡(luò)的化學(xué)計(jì)量矩陣，矩陣中的每一行表示一個(gè)代謝物，每一欄表示一個(gè)代謝反應(yīng)，每個(gè)元素表示該代謝物在該代謝反應(yīng)中的化學(xué)計(jì)量系數(shù)。當(dāng)未施加約束條件時(shí)，代謝網(wǎng)絡(luò)的通量可分布于解空間內(nèi)的任意點(diǎn)，若設(shè)定反應(yīng)邊界值作為約束條件后，即可得到一個(gè)有限的解空間。最后給定某目標(biāo)函數(shù)，用線性規(guī)劃的方法尋找其最優(yōu)解。本研究進(jìn)行的最優(yōu)生長率、模擬動(dòng)態(tài)生長情況、基因刪除、功能都通過FBA過程實(shí)現(xiàn)。FBA過程只返回最優(yōu)解，從而使有關(guān)整個(gè)解決方案的空間未能體現(xiàn)。為解決該問題，COBRA軟件使用抽樣對解空間進(jìn)行描述。

1.4 OptKnock 過程［10］

類似“生物工程”，COBRA軟件通過對整個(gè)代謝網(wǎng)絡(luò)進(jìn)行分析，設(shè)定某個(gè)代謝物為目標(biāo)產(chǎn)物和允許刪除反應(yīng)的數(shù)量為約束條件，獲得為使目標(biāo)產(chǎn)物產(chǎn)量最大，何種或哪些特定反應(yīng)可刪除的解決方案。

2 結(jié)果

2.1 結(jié)核桿菌H37Rv菌株的GSMN模型iNJ661描述

模型中每個(gè)代謝反應(yīng)均有一個(gè)范圍從0到4的得分來反映可用的信息和證據(jù)，其中4分為最高，代表有實(shí)驗(yàn)的生化數(shù)據(jù)來支持。一個(gè)模型反應(yīng)集合的平均得分水平代表該模型的真實(shí)性，iNJ661的平均得分為2.31分。該模型基因 -蛋白質(zhì) -反應(yīng)(Gene-protein-reaction，GPR)相互聯(lián)系的比例占76.76%，表明該模型完整性較好。應(yīng)用 Gap Filling［11］過程對該模型進(jìn)行檢測，僅發(fā)現(xiàn)37個(gè)不能被生成的代謝物和48個(gè)缺乏所需反應(yīng)物的代謝物，這將有待于研究人員進(jìn)行完善。

2.2 iNJ661的生長表型

iNJ661模型提供了87個(gè)可供調(diào)整的培養(yǎng)基成分，其預(yù)設(shè)的培養(yǎng)基為Middlebrook 7H9。結(jié)核桿菌生長緩慢且難于培養(yǎng)，選擇或配制一種促進(jìn)結(jié)核桿菌生長的培養(yǎng)基顯得尤其重要。模擬菌株在Middlebrook 7H9［12］、Youmans［13］、CAMR［14］三種培養(yǎng)基生長，培養(yǎng)基成分和對應(yīng)的生長率見表1。

表1 Middlebrook 7H9、Youmans和 CAMR培養(yǎng)基成分和生長率Table 1 Composition and biomass accumulation of Middlebrook 7H9，Youmans and CAMR

分別向模型中加入不同氨基酸和碳源來觀察其對菌株生長率的影響，結(jié)果顯示在加入絲氨酸(Ser)、天冬氨酸(Asn)、谷氨酸(Glu)足夠讓菌株快速生長;以麥芽糖(Maltose)和海藻糖(Trehalose)作為碳源加入到培養(yǎng)基中，菌株生長最快(圖1和圖2)。

圖1 單獨(dú)加入不同氨基酸時(shí)iNJ661的生長率Fig.1 Biomass accumulation of iNJ661 model when different amino acids added singly

同理，通過改變預(yù)設(shè)培養(yǎng)基中每個(gè)成分來觀察相應(yīng)的生長率，結(jié)果顯示去掉銨鹽(Ammonium)、三價(jià)鐵鹽(Ferric iron)、磷酸鹽(Phosphate)、硫酸鹽(Sulfate)、甘油(Glycerol)，iNJ661生長將會(huì)受到阻礙甚至無法生長(表2)，這提示以上5種代謝物的轉(zhuǎn)運(yùn)體和激酶都可能成為藥物的靶標(biāo)，達(dá)到抑菌甚至殺菌的效果。

表2 單獨(dú)改變Middlebrook 7H9每一成分對應(yīng)的iNJ661生長率Table 2 Corresponding biomass accumulation of iNJ661 model with single change in Middlebrook 7H9

2.3 解空間抽樣

在模擬H37Rv菌株生長表型的步驟中，已經(jīng)發(fā)現(xiàn)從培養(yǎng)基中去掉磷酸鹽會(huì)使模型的生長受到明顯的阻礙，現(xiàn)通過抽樣直方圖來觀察一系列酶的活性(圖3)。結(jié)果顯示十種酶都在不同程度上受到了抑制，其中丙糖磷酸異構(gòu)酶(Triosephosphate isomerase，TPI)、3-磷酸甘油醛脫氫酶(Glyceraldehyde-3-phosphate dehydrogenase，GAPD)、磷酸甘油酸變位酶 (Phosphoglycerate mutase， PGM)、烯 醇 酶(Enolase，ENO)受限明顯。

圖3 培養(yǎng)基中去掉磷酸鹽前后不同反應(yīng)通量下蛋白酶濃度的抽樣直方圖Fig.3 Sample histogram of proteinase concentration under differentreaction flux in medium with phosphate vs.without phosphate

2.4 基因刪除

隨著多重耐藥菌株的不斷出現(xiàn)，識別更多靶標(biāo)迫在眉睫。通過FBA過程完成iNJ661的基因刪除學(xué)習(xí)，結(jié)果顯示模型包含188個(gè)必要基因，其中ilvD(Singh V，et al.2011)、rmlA(Qu H，et al.2007)，inhA(Banerjee A，et al.1994)、kasA(Lee W，et al.2011)、embC/embA/embB(Goude R，et al.2008/Amin A，et al.2008/Li W，et al.2010)等已被研究證實(shí)。另外在非必要基因的雙致死性基因?qū)W習(xí)中獲得16個(gè)基因?qū)?表3)。

表3 INJ661模型非必要基因中16個(gè)致死基因?qū)able 3 16 pairs of non-essential genes in iNJ661 model whose simultaneous deletion is lethal

2.5 OptKnock 過程

研究發(fā)現(xiàn)結(jié)核桿菌的致病性主要與占細(xì)胞壁干重 60%以上的脂質(zhì)成分(如 Phthiocerol dimycocerosate，PDIM) 密 切相 關(guān)［15］。本 文 應(yīng) 用OptKnock過程，將原有的培養(yǎng)基葡萄糖供給調(diào)整為10 mmol·gDW-1·h-1，分別將目標(biāo)設(shè)定為 PPDIM(Phenol phthiocerol dimycocerosate)和PDIM以求其最大產(chǎn)量(表4)，結(jié)果顯示抑制ACGS一個(gè)反應(yīng)可以讓PPDIM產(chǎn)量最大化，抑制ACDS、GLCS1、ICL三個(gè)反應(yīng)可讓PDIM產(chǎn)量最大。通過以上模擬實(shí)驗(yàn)，提示可能短時(shí)間獲得更多的PPDIM和PDIM，以進(jìn)行其致病和毒力損傷機(jī)制研究，不僅加深對病原體-宿主互作網(wǎng)絡(luò)中的理解，還有望為結(jié)核病的治療帶來新的希望與突破。

表4 OptKnock過程Table 4 OptKnock process model whose simultaneous deletion is lethal

3 討論

本研究通過仿真分析發(fā)現(xiàn)iNJ661模型的動(dòng)態(tài)生長情況與H37Rv菌株的實(shí)驗(yàn)數(shù)據(jù)相一致。以絲氨酸、天冬氨酸、谷氨酸作為氨基酸來源，以麥芽糖、海藻糖作為碳源，菌株增長明顯。去除培養(yǎng)基中的銨鹽、三價(jià)鐵鹽、磷酸鹽、硫酸鹽和甘油，菌株的生長將受限。解空間抽樣分析在蛋白酶水平解釋關(guān)鍵的營養(yǎng)目標(biāo)(如磷酸鹽)對菌株生長的影響，后續(xù)研究若設(shè)計(jì)相關(guān)的酶抑制劑很有可能達(dá)到抑菌甚至殺菌的效果。

基因刪除學(xué)習(xí)中獲得的部分必要基因與Sassetti等［16］在體的轉(zhuǎn)座子定點(diǎn)雜交實(shí)驗(yàn)結(jié)果相符，并與目前H37Rv部分的藥物靶標(biāo)相一致，若結(jié)合與人類基因的同源性分析，研究發(fā)現(xiàn)的必要基因和基因?qū)τ锌赡転樾碌目菇Y(jié)核病藥物篩選提供提示和借鑒。OptKnock過程給出量產(chǎn)致病物質(zhì)PDIM的代謝工程途徑，進(jìn)而模擬其他重要物質(zhì)量產(chǎn)機(jī)制，指導(dǎo)疫苗開發(fā)和病原微生物研究人員的工作。需要注意的是建模進(jìn)行的實(shí)驗(yàn)是基于現(xiàn)有代謝網(wǎng)絡(luò)模型的，其陽性結(jié)果很有可能屬存?zhèn)巍?/p>

本研究簡要提供了對結(jié)核桿菌H37Rv菌株進(jìn)行約束建模分析的基本分析框架，可為后續(xù)相關(guān)研究提供參考和借鑒。

References)

［1］ DIAS H M，F(xiàn)ALZON D，F(xiàn)ITZPATRICK C，et al.Global tuberculosis report 2012［M］.Geneva:WHO Press，2012，8-28.

［2］ HINGLEY-WILSON S M，CASEY R，CONNELL D，et al.Undetected multidrug-resistant tuberculosis amplified by first-line therapy in mixed infection［J］.Emerging Infectious diseases，2013，19(7):1138-1141.

［3］ DUROT M，BOURGUIGNON P Y，SCHACHTER V.Genome-scale models of bacterial metabolism:reconstruction and applications［J］.FEMS Microbiology Reviews，2009，33(1):164-190.

［4］ SCHELLENBERGER J，PARK J O，CONRAD T M，et al.BiGG: a BiochemicalGenetic and Genomic knowledgebase of large scale metabolic reconstructions［J］.BMC Bioinformatics，2010，11:213.

［5］ FANG X，WALLQVIST A，REIFMAN J.Modeling phenotypic metabolic adaptations of Mycobacterium tuberculosis H37Rv under hypoxia［J］.PLoS Computational Biology，2012，8(9):e1002688.

［6］ HUCKA M，F(xiàn)INNEY A，SAURO H M，et al.SBML Forum.The systems biology markup language(SBML):a medium for representation and exchange of biochemical network models［J］.Bioinformatics，2003，19(4):524-531.

［7］ SCHELLENBERGER J，QUE R，F(xiàn)LEMING R M，et al.Quantitative prediction of cellular metabolism with constraint-based models:the COBRA Toolbox v2.0 ［J］.Nature Protocols，2011，6(9):1290-1307.

［8］ METRIS A，GEORGE S，BARANYI J.Modelling osmotic stress by Flux Balance Analysis at the genomic scale［J］.International Journal of Food Microbiology，2012，152(3):123-128.

［9］ RAMAN K，CHANDRA N.FLUX balance analysis of biological systems:applications and challenges ［J］.Briefings in Bioinformatics，2009，10(4):435-449.

［10］ BURGARD A P，PHARKYA P，MARANAS C D.Optknock:a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization［J］.Biotechnology and Bioengineering，2003，84(6):647-657.

［11］ ORTH J D，PALSSON B ?.Gap-filling analysis of the iJO1366 Escherichia coli metabolic network reconstruction for discovery of metabolic functions［J］.BMC Systems Biology，2012，6:30.

［12］ PARK Y K，KOH W J，KIM S O，et al.Clarithromycin susceptibility testing of Mycobacterium avium complex using 2，3-diphenyl-5-thienyl-(2)-tetrazolium chloride microplate assay with Middlebrook 7H9 broth ［J］.Journal of Korean Medical Science，2009，24(3):511-512.

［13］GADRE D V，TALWAR V.In vitro susceptibility testing of Mycobacterium tuberculosis strains to trifluoperazine［J］.Journal of Chemotherapy，1999，11(3):203-206.

［14］JAMES B W，WILLIAMS A，MARSH P D.The physiology and pathogenicity of Mycobacterium tuberculosis grown under controlled conditions in a defined medium［J］.Journal of Applied Microbiology，2000，88(4):669-677.

［15］ DOMENECH P，REED M B.Rapid and spontaneous loss of phthiocerol dimycocerosate(PDIM)from Mycobacterium tuberculosis grown in vitro:implications for virulence studies［J］.Microbiology，2009，155(Pt 11):3532-3543.

［16］ SASSETTI C M，RUBIN E J.Genetic requirements for mycobacterial survival during infection［J］.Proceedings of the National Academy of Sciences of the United States of America，2003，100(22):12989-12994.