柴文平， 薛 衛(wèi)， 張 梁*， 石貴陽

（1.糧食發(fā)酵工藝與技術(shù)國家工程實(shí)驗(yàn)室 江南大學(xué)，江蘇 無錫214122；2.南京農(nóng)業(yè)大學(xué) 信息科學(xué)技術(shù)學(xué)院，江蘇 南京 210095）

隨著物種基因組測(cè)序的完成以及大量生物學(xué)數(shù)據(jù)的產(chǎn)生，系統(tǒng)生物學(xué)研究技術(shù)也日益成熟。系統(tǒng)生物學(xué)能夠模擬和推測(cè)復(fù)雜生物體行為，其中網(wǎng)絡(luò)模型模擬是最主要的模擬方法[1-2]。系統(tǒng)生物學(xué)的網(wǎng)絡(luò)模型種類包含代謝網(wǎng)絡(luò)[3]、轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)[4]、信號(hào)轉(zhuǎn)導(dǎo)網(wǎng)絡(luò)[5]和轉(zhuǎn)錄翻譯網(wǎng)絡(luò)[6]等。而其中的基因組規(guī)模代謝網(wǎng)絡(luò)模型（Genome Scale Metabolic Model，GSMM）已經(jīng)成為系統(tǒng)生物學(xué)不可或缺的研究工具。它通過整合基因組學(xué)、文獻(xiàn)組學(xué)、蛋白質(zhì)組學(xué)等組學(xué)數(shù)據(jù)，建立由基因-蛋白質(zhì)-生化反應(yīng)（G-P-R）關(guān)聯(lián)組成的特定生物代謝網(wǎng)絡(luò)，是從全局深刻理解其生理特性與定向調(diào)控工業(yè)微生物生理功能的重要平臺(tái)[7]。

截止到2013年11月，已經(jīng)公布了98個(gè)生物全基因組代謝網(wǎng)絡(luò)模型，其中細(xì)菌65個(gè)，古細(xì)菌6個(gè)，真核生物28個(gè)[8]。而GOLD數(shù)據(jù)庫公布的已完成測(cè)序物種包含2 698個(gè)，其中細(xì)菌2 384個(gè)，古細(xì)菌163個(gè)，真核生物151個(gè)[9]。前者數(shù)量遠(yuǎn)遠(yuǎn)小于后者，造成這種情況的原因除了對(duì)很多物種生理生化機(jī)制了解較少之外，更重要的是重構(gòu)代謝網(wǎng)絡(luò)過程需要大量的人工，非常耗時(shí)耗力[10]。雖然已經(jīng)出現(xiàn)了一些通過自動(dòng)獲取信息初步重構(gòu)網(wǎng)絡(luò)的軟件平臺(tái)，如SEED服務(wù)器與GEM System軟件等[11-13]，但其中仍然需要大量的人工操作與修正工作，如SEED注釋的基因序列需要人工進(jìn)行，同NCBI注釋的基因序列相一致，以便后續(xù)其他構(gòu)建方法的補(bǔ)充與修正。此外，大量生物信息學(xué)數(shù)據(jù)庫的出現(xiàn)也對(duì)進(jìn)一步實(shí)現(xiàn)網(wǎng)絡(luò)自動(dòng)化重構(gòu)提出了挑戰(zhàn)。因此，代謝網(wǎng)絡(luò)重構(gòu)的自動(dòng)化研究，已成為推動(dòng)代謝網(wǎng)絡(luò)發(fā)展的重大課題。

文中以構(gòu)建樹干畢赤酵母（Scheffersomyces stipitis或者 Pichia stipitis）[14-15]CBS 6054的基因組規(guī)模代謝網(wǎng)絡(luò)模型為例，以簡單、面向?qū)ο蟮腏ava語言為基礎(chǔ)，對(duì)代謝網(wǎng)絡(luò)自動(dòng)化重構(gòu)的方法進(jìn)行了研討，提出了一種基于KEGG在線數(shù)據(jù)庫來自動(dòng)化構(gòu)建初模型的方法，并對(duì)基于Uniprot-MetaCyc本地?cái)?shù)據(jù)庫，以及親緣物種同源比對(duì)構(gòu)建初模型的方法，和整合過程進(jìn)行了自動(dòng)化研究，達(dá)到代謝網(wǎng)絡(luò)模型構(gòu)建的最大程度自動(dòng)化。自動(dòng)重構(gòu)流程如圖1所示。

圖1 模型自動(dòng)重構(gòu)流程Fig.1 Process of the auto-reconstruction of model

1 代謝網(wǎng)絡(luò)自動(dòng)重構(gòu)

基因組規(guī)模代謝網(wǎng)絡(luò)模型構(gòu)建過程主要涉及[16]：初模型反應(yīng)列表的獲得；模型精細(xì)化；轉(zhuǎn)換數(shù)學(xué)模型；模擬與應(yīng)用。一般認(rèn)為構(gòu)建出G-P-R反應(yīng)列表就等于完成了代謝網(wǎng)絡(luò)的初模型[10]。雖然構(gòu)建G-P-R反應(yīng)列表過程比較容易理解，但由于基因組規(guī)模的大量生物信息學(xué)數(shù)據(jù)涌現(xiàn)，反應(yīng)列表的構(gòu)建反而成為構(gòu)建過程中最繁瑣耗時(shí)耗力的一部分[17]。所以借助計(jì)算機(jī)技術(shù)來達(dá)到高效構(gòu)建反應(yīng)列表、提高代謝網(wǎng)絡(luò)模型構(gòu)建效率的作用。

KEGG （KyotoEncyclopediaofGenesand Genomes）與 MetaCyc（http：//metacyc.org/）是模型構(gòu)建中最常用的網(wǎng)絡(luò)數(shù)據(jù)庫，包含了物種的基因與基因組、酶促反應(yīng)及其代謝途徑和化合物等相關(guān)信息[18-19]，是網(wǎng)絡(luò)構(gòu)建過程中網(wǎng)絡(luò)數(shù)據(jù)信息的主要來源。而同源比對(duì)這一比較基因組學(xué)策略可以快速找到親緣菌株之間的遺傳關(guān)系及對(duì)應(yīng)的生化信息，這一策略是構(gòu)建全新微生物代謝網(wǎng)絡(luò)模型的可靠信息來源。

1.1 基于KEGG在線數(shù)據(jù)庫的模型構(gòu)建

KEGG作為代謝網(wǎng)絡(luò)構(gòu)建常用數(shù)據(jù)庫，其包含有多個(gè)在線子數(shù)據(jù)庫，其中REACTION數(shù)據(jù)庫包含迄今為止發(fā)現(xiàn)的所有生化反應(yīng)[18]。各個(gè)子數(shù)據(jù)庫的網(wǎng)頁數(shù)據(jù)格式比較統(tǒng)一明確，方便人們進(jìn)行遠(yuǎn)程服務(wù)器訪問。但是，KEGG數(shù)據(jù)庫更新頻繁，各個(gè)子數(shù)據(jù)庫不能夠免費(fèi)下載，需要付費(fèi)使用。而在重構(gòu)基因組代謝網(wǎng)絡(luò)過程中，因?yàn)閿?shù)據(jù)信息量浩大，頻繁訪問遠(yuǎn)程服務(wù)器比較耗時(shí)耗力。因此，實(shí)現(xiàn)一種批量在線獲取并存取數(shù)據(jù)的方法意義重大。

1.1.1 方法概述 超文本轉(zhuǎn)移協(xié)議（Hypertext transfer protocol，HTTP）是一種詳細(xì)規(guī)定了瀏覽器和萬維網(wǎng)服務(wù)器之間互相通信的規(guī)則，通過因特網(wǎng)傳送萬維網(wǎng)文檔的數(shù)據(jù)傳送協(xié)議。KEGG提供物種特異性基因組信息以及所有反應(yīng)式信息查詢網(wǎng)頁，通過一定的 URL（Uniform Resource Locator，統(tǒng)一資源定位符）格式地址發(fā)送HTTP請(qǐng)求，返回網(wǎng)頁html腳本含有基因組信息或者反應(yīng)式相關(guān)信息。html腳本由標(biāo)題、js代碼、正文、相關(guān)鏈接、聲明等區(qū)域組成，而有用信息只出現(xiàn)在正文欄的＜table＞標(biāo)記內(nèi)。如果對(duì)于每個(gè)查詢網(wǎng)頁均全面分析，將大大降低效率。因此作者提出基于正文先驗(yàn)位置的網(wǎng)頁分析方法，獲取html的正文信息中的＜form＞標(biāo)志后第一個(gè)＜table＞在html腳本字符串中的起始位置begin_pos，＜/form＞標(biāo)記前最近一個(gè)＜/table＞結(jié)束位置end_pos，begin_pos至end_pos即正文先驗(yàn)位置，只處理begin_pos至end_ps腳本串信息即可獲得反應(yīng)等信息。具體可采用JAVA控件NekoHTML分析html腳本中每個(gè)節(jié)點(diǎn)數(shù)據(jù)，并用正則表達(dá)式提取相關(guān)信息。

1.1.2 算法實(shí)現(xiàn)

1）確定KO以及R號(hào)：提交物種基因組蛋白質(zhì)序列至KAAS自動(dòng)注釋服務(wù)器，下載KEGG BRITE Database中KO-R列表至本地為K-R.xls，依次讀取返回的KO號(hào)，自K-R.xls確定相應(yīng)的R號(hào)，寫入KR.xls中。

2）依據(jù)KO獲得基因蛋白質(zhì)信息：

向反應(yīng)式信息查詢網(wǎng)頁URL地址發(fā)送HTTP請(qǐng)求。

服務(wù)器響應(yīng)代碼串為Gene_string，分析提取html正文信息中的＜form＞標(biāo)志后的＜DIV＞節(jié)點(diǎn)，節(jié)點(diǎn)數(shù)為DIV_number。

設(shè)n=1。對(duì)第n個(gè)節(jié)點(diǎn)，依據(jù)JAVA正則表達(dá)式"K+\d{5}"提取對(duì)應(yīng)的KO號(hào)。

讀取KR.xls中的每個(gè)KO，與上步中KO比較，若相等，設(shè)置JAVA正則表達(dá)式"sign+\_+\d+"提取 GENE，設(shè)置 JAVA 正則表達(dá)式"\s+（\w+）+\；|\s+（\w++\.+\w+）+\；"提取 PROTEIN。一并寫入表中KO對(duì)應(yīng)行。

n增加1，如果n≤DIV_number，重新執(zhí)行；否則，結(jié)束。

3）依據(jù)R號(hào)獲得反應(yīng)式及相關(guān)信息：

讀取上步生成的KR.xls。

設(shè)i=1。

讀取KR.xls的第i行R值，設(shè)置KEGG服務(wù)器訪 問 地 址 為 “http：//www.genome.jp/dbget-bin/www_bget?+R值”，發(fā)送http請(qǐng)求。

獲取html格式腳本，服務(wù)器響應(yīng)代碼串為reaction_string。如果i=1，首先計(jì)算正文先驗(yàn)位置，在reaction_string中分析獲取html的正文信息中的＜form＞標(biāo)志后第一個(gè)＜table＞在html腳本字符串中的起始位置begin_pos并保存。在reaction_string串中查詢獲得 end_pos值，得到＜table＞…＜/table＞內(nèi)容字符串content_string。

用NekoHTML讀取content_string中Name、Definition、Equation、Enzyme、Pathway、Orthology 字段值并寫入該表對(duì)應(yīng)行中。

i增 1，如果 i≤“P-K.xls表行數(shù)”，重新執(zhí)行；否則，結(jié)束。

4）實(shí)例結(jié)果：以樹干畢赤酵母為例，實(shí)現(xiàn)流程如圖2所示。系統(tǒng)運(yùn)行大約30 min，得出的初模型包含786條反應(yīng)及關(guān)聯(lián)基因（549個(gè)）、酶及代謝途徑等信息，不僅節(jié)省了大量的人力和時(shí)間，而且保證了數(shù)據(jù)的最新性。

圖2 KEGG網(wǎng)頁挖掘數(shù)據(jù)實(shí)現(xiàn)流程Fig.2 Procedure of the excavation of date based on the web page of KEGG

1.2 基于Uniprot-MetaCyc本地?cái)?shù)據(jù)庫的模型構(gòu)建

MetaCyc數(shù)據(jù)庫具有可信度高、信息全面、易訪問、免費(fèi)用于學(xué)術(shù)研究等優(yōu)點(diǎn)[20]。相較于KEGG來說，也存在格式繁復(fù)、網(wǎng)頁打開速度慢等缺點(diǎn)，所以不容易網(wǎng)頁查找與數(shù)據(jù)抓取。因此可構(gòu)建MetaCyc本地?cái)?shù)據(jù)庫進(jìn)行數(shù)據(jù)提取。MetaCyc數(shù)據(jù)同樣存在無物種特異性問題，只能依據(jù)物種特異性信息從數(shù)據(jù)庫中查找提取特定的反應(yīng)。UniProt[21]是信息最豐富、資源最廣的蛋白質(zhì)數(shù)據(jù)庫。它包含高質(zhì)量的、手工注釋的、非冗余的數(shù)據(jù)集?？梢栽赨niProtKB中找到物種特異性的基因組注釋的蛋白質(zhì)相關(guān)詳細(xì)信息。

1.2.1 方法概述 通過UniProt中的酶號(hào)（EC number）搭建一條連接UniProt與MetaCyc數(shù)據(jù)庫的橋梁，從而獲取物種特異性的基因組反應(yīng)信息。UniProt數(shù)據(jù)庫（http：//www.uniprot.org/）中下載物種特異性基因組蛋白質(zhì)注釋信息至本地，建立Excel數(shù)據(jù)庫DBU。下載MetaCyc數(shù)據(jù)庫至本地，建立本地Excel數(shù)據(jù)庫，其中包含有EC和酶學(xué)反應(yīng)號(hào)（ERN）的DBMR子數(shù)據(jù)庫和有ERN和對(duì)應(yīng)酶學(xué)反應(yīng)式（ER）的DBME子數(shù)據(jù)庫。

設(shè)DBU中有N條有EC的蛋白質(zhì)信息，N＞1。DBMReaction中有M條有EC和對(duì)應(yīng)ERN的信息，M＞1。DBMEnzymes有 K條 ERN對(duì)應(yīng)的 ER的信息，K＞1。 Ei=（ECi，Gi，Pi）表示第 i條 EC 對(duì)應(yīng)的基因（G）和蛋白質(zhì)（P）信息，1≤i≤N；Fi=（ECi，ERNi）表示第 i條 EC 對(duì)應(yīng)的 ERN 信息；Hi=（ERNi，ERi）表示第i 條 ERN 對(duì)應(yīng)的 ER 信息；EFHi=（ECi，Gi，Pi，ERNi，ERi）表示第i條初模型反應(yīng)信息（酶號(hào)、基因、蛋白質(zhì)、酶學(xué)反應(yīng)編號(hào)、酶學(xué)反應(yīng)式）。求取某個(gè)EC對(duì)應(yīng)的反應(yīng)信息表示為 EFHi=Ei∪Fi∪Hi，1≤i≤N。 設(shè)EFH0=（EC0，G0，P0，ERN0，ER0） 表示整個(gè)初模型，則將問題轉(zhuǎn)化為求解各反應(yīng)信息的合集，即求解EFH0=∑EFHi。

1.2.2 方法實(shí)現(xiàn) Uniprot與MetaCyc數(shù)據(jù)庫數(shù)據(jù)均有文獻(xiàn)支持、手工注釋錄入的特點(diǎn)，可信度較高。其中，MetaCyc數(shù)據(jù)庫數(shù)據(jù)為無物種特異性的酶及酶學(xué)反應(yīng)信息。為了獲得物種特異性的G-P-R關(guān)系列表，則需由Uniprot獲得物種特異性基因組蛋白質(zhì)注釋信息，依據(jù)信息中的酶號(hào)信息自MetaCyc數(shù)據(jù)庫獲得相應(yīng)的酶學(xué)反應(yīng)信息。初模型G-P-R列表獲得流程如圖3所示。

圖3 基于UniProt-MetaCyc的G-P-R獲取流程Fig.3 Process of getting the G-P-R list based on UniProt-MetaCyc

1.3 基于同源物種的模型構(gòu)建

比較基因組學(xué)的策略可以快速找到兩物種間的遺傳關(guān)系與對(duì)應(yīng)的生化反應(yīng)信息?；诰植勘葘?duì)算法的搜索工具 （Basic Local Alignment Search Tool，BLAST）[22]對(duì)目的生物與其親緣物種的基因組序列信息進(jìn)行雙向比對(duì)，可推測(cè)大部分同源序列的基因功能[23]。在比對(duì)過程中蛋白質(zhì)序列比基因序列有更高的保守性，所以基于蛋白質(zhì)序列的同源比對(duì)有效性也更高。在蛋白質(zhì)序列雙向比對(duì)（Blastp）過程中，通過設(shè)置一定的期望值（E）、相似度（Identity）與匹配序列長度來確定兩者是否為同源序列。當(dāng)確定兩個(gè)蛋白質(zhì)序列為同源序列后，可推測(cè)兩者具有相似的酶學(xué)功能，進(jìn)而推測(cè)其可以作用于同一反應(yīng)。

1.3.1 方法概述 通過親緣關(guān)系樹的建立，查找目的生物的親緣物種，并下載親緣物種的高質(zhì)量代謝網(wǎng)絡(luò)模型，建立本地?cái)?shù)據(jù)庫。依據(jù)目的生物與親緣物種的基因組蛋白質(zhì)序列雙向比對(duì)結(jié)果，和每條序列的登錄號(hào)（Accession Number），自動(dòng)獲得蛋白質(zhì)對(duì)應(yīng)的基因，進(jìn)而根據(jù)基因-酶-反應(yīng)的關(guān)聯(lián)與親緣物種的全基因組規(guī)模代謝網(wǎng)絡(luò)模型，獲得目的生物的規(guī)模代謝網(wǎng)絡(luò)模型。利用Java分別調(diào)用Poi_3.7.jar和tm_extractors_0.4.jar工具包函數(shù)，對(duì)excel和word進(jìn)行操作。

1.3.2 方法實(shí)現(xiàn) 以樹干畢赤酵母CBS 6054為例，選擇親緣菌株為巴斯德畢赤酵母（Pichia pastoris）和釀酒酵母（Saccharomyces cerevisiae），并分別以釀酒酵母的iMM904[24]與iLL672[25]、巴斯德畢赤酵母的iPP668[26]與iLC915[27]4個(gè)高質(zhì)量模型數(shù)據(jù)（excel）為同源比對(duì)模型。利用BLAST軟件對(duì)本地?cái)?shù)據(jù)庫蛋白質(zhì)序列文件進(jìn)行雙向Blastp，設(shè)置篩選E值為1e-40。利用JAVA語言編寫程序?qū)π畔⒘烤薮蟮膚ord結(jié)果文件進(jìn)行篩選處理。同時(shí)，基于登錄號(hào)自對(duì)應(yīng)的基因蛋白質(zhì)列表中獲得該序列對(duì)應(yīng)的基因與蛋白質(zhì)信息。依據(jù)比對(duì)上的基因，進(jìn)一步替換iMM904、iLL672、iPP668與 iLC915 4個(gè)模型中的基因，獲得4個(gè)樹干畢赤酵母基因組規(guī)模代謝網(wǎng)絡(luò)初模型，分別包含有809、613、1 105、1 413條反應(yīng)。

1.4 初模型的比較和整合方法簡述

上述3種自動(dòng)化重構(gòu)方法構(gòu)建的樹干畢赤酵母各個(gè)初模型，除了包含G-P-R列表之外，還包含一些構(gòu)建精細(xì)模型的附加信息。不同初模型內(nèi)容比較如表1所示。

表1 自動(dòng)化重構(gòu)的不同模型比較Table 1 Comparison of the different auto-reconstruction models

基于上述3種自動(dòng)化方法構(gòu)建的初模型中包含有大量重復(fù)反應(yīng)信息，需對(duì)初模型進(jìn)行整合獲得格式統(tǒng)一、信息量充足、無重復(fù)反應(yīng)的反應(yīng)列表[28]。整合過程中遇到的最大問題就是反應(yīng)式中化合物格式不統(tǒng)一，一般表示為4種情況：親緣物種替換的模型中同一種化合物以不同的簡寫形式表達(dá)；配位化合物由于配位單元順序不同，造成了名稱的不同；某些化合物為同一種物質(zhì)，可能有不同的名稱；化合物中出現(xiàn)一些符號(hào)，使得不能夠識(shí)別為同一種化合物。由此，作者提出了基于化合物數(shù)據(jù)庫和反應(yīng)式字符頻度直方圖特征比對(duì)兩種方法來達(dá)到模型反應(yīng)自動(dòng)化整合。

1.4.1 基于化合物數(shù)據(jù)庫的模型整合 構(gòu)建了整合KEGG和MetaCyc的化合物數(shù)據(jù)庫。KEGG Compound列表中包括每個(gè)化合物的KEGG ID，即C號(hào)（字母C+5個(gè)數(shù)字）以及對(duì)應(yīng)的化合物不同的名稱。MetaCyc數(shù)據(jù)庫中也包含了每個(gè)化合物不同的表達(dá)形式，以及與KEGG數(shù)據(jù)庫相對(duì)應(yīng)的C號(hào)。兩者中部分化合物有數(shù)據(jù)庫特異性，所以整合兩個(gè)數(shù)據(jù)庫建立本地化合物數(shù)據(jù)庫。

編寫程序?qū)φ夏Ｐ椭械拿總€(gè)化合物，建立對(duì)化合物數(shù)據(jù)庫的映射關(guān)系。即查找替換對(duì)應(yīng)的C號(hào)，然后比對(duì)C號(hào)形式的反應(yīng)式中反應(yīng)物與生成物的異同來確定是否為同一反應(yīng)（反應(yīng)式中化合物均為簡寫表達(dá)形式的模型數(shù)據(jù)，需先替換為代謝物列表中對(duì)應(yīng)的全稱形式）。

基于化合物數(shù)據(jù)庫的模型整合并不能映射全部相同反應(yīng)，原因可分為以下3種：相同反應(yīng)因H或H2O的缺失或冗余而不能判定為同一反應(yīng)；某些配位化合物由于配位單元順序不同，或者化合物中一些符號(hào)而造成了不能成功替換C號(hào)；反應(yīng)方向不同造成反應(yīng)物與產(chǎn)物識(shí)別錯(cuò)誤。

1.4.2 基于字符頻度直方圖特征的模型整合 針對(duì)上述情況，作者提出了一種判斷兩個(gè)反應(yīng)是否為相同反應(yīng)的新方法。該方法通過提取反應(yīng)式字符頻度直方圖特征，進(jìn)一步計(jì)算直方圖間的馬氏距離[29]來實(shí)現(xiàn)?；衔锏幕瘜W(xué)式核心由英文字母與阿拉伯?dāng)?shù)字組成，建立每個(gè)反應(yīng)式的直方圖，26個(gè)英文字母與阿拉伯?dāng)?shù)字（0—9）為橫坐標(biāo)，36個(gè)元素出現(xiàn)頻次為縱坐標(biāo)，計(jì)算直方圖之間的距離（馬氏距離），設(shè)置一定的閾值，當(dāng)距離小于該閾值時(shí)表示為相同反應(yīng)。

在進(jìn)行判斷之前，程序先對(duì)化合物區(qū)室標(biāo)志如[c]或[m]進(jìn)行比較，因?yàn)橥环磻?yīng)可能存在不同的細(xì)胞區(qū)室，若區(qū)室相同依據(jù)反應(yīng)式直方圖特征進(jìn)行馬氏距離計(jì)算。在統(tǒng)計(jì)過程中，自動(dòng)將大寫字母轉(zhuǎn)化為小寫字母進(jìn)行統(tǒng)計(jì)，且自動(dòng)忽略標(biāo)點(diǎn)符號(hào)或特殊字符（如，、-等）等。在計(jì)算之前還應(yīng)移去無意義的干擾詞。干擾詞分為3類：不存在生物意義的詞匯（如 alpha，beta等）；出現(xiàn)頻度非常大的詞匯（如 a、the等）；單詞長度小于2的詞匯。

馬氏距離表示數(shù)據(jù)的協(xié)方差距離，它是一種有效的計(jì)算兩個(gè)未知樣本集的相似度的方法。設(shè)X{X1，X2…Xn}和 Y{Y1，Y2…Yn}為總體中抽取的樣本，則X、Y兩組樣本之間的馬氏距離為[30]

針對(duì)上述基于化合物數(shù)據(jù)庫不能匹配的相同反應(yīng)出現(xiàn)的原因，文中在計(jì)算馬氏距離的時(shí)候采用模糊匹配值。當(dāng)距離值相差2.828之內(nèi)被認(rèn)為是相同的，否則是不同的。當(dāng)判斷兩個(gè)反應(yīng)為相同反應(yīng)后，程序接著判斷該反應(yīng)對(duì)應(yīng)的基因是否相同，若不同，則將不同的基因保存于同一反應(yīng)后，方便下一步進(jìn)行人工判斷取舍。

1.4.3 實(shí)施實(shí)例 初模型整合過程中，對(duì)相同反應(yīng)的取舍依據(jù)各個(gè)初模型信息全面的高低：同源比對(duì)模型＞KEGG模型＞MetaCyc模型，但整合后的化合物格式統(tǒng)一以KEGG為標(biāo)準(zhǔn)，以同一C號(hào)對(duì)應(yīng)不同的化合物表達(dá)形式的第一個(gè)為準(zhǔn)。以樹干畢赤酵母為例，6個(gè)初模型基于化合物數(shù)據(jù)庫整合后的模型包含1 878條反應(yīng)以及對(duì)應(yīng)的956個(gè)基因，經(jīng)字符頻度直方圖特征的模型整合方法補(bǔ)充，最終獲得1 531條反應(yīng)。經(jīng)人工逐一檢查，發(fā)現(xiàn)無重復(fù)反應(yīng)，結(jié)果理想。結(jié)果模型中大部分反應(yīng)來自基于同源比對(duì)構(gòu)建的模型，所以其中囊括了精細(xì)模型包含的基因、反應(yīng)、酶、酶號(hào)、代謝途徑、反應(yīng)方向、亞細(xì)胞定位，以及代謝物的附加信息等內(nèi)容。整合流程如圖4所示。

圖4 初模型自動(dòng)整合流程Fig.4 Auto-integration process of the draft model

2 核心反應(yīng)的識(shí)別

構(gòu)建基因組規(guī)模代謝網(wǎng)絡(luò)模型的核心代謝途徑就是糖代謝，不管是脂肪酸還是蛋白質(zhì)，或者是多糖，最終都要轉(zhuǎn)化為糖類進(jìn)行能量代謝與產(chǎn)物合成。EMP糖酵解（Glycolysis）是糖代謝過程的第一步，存在于所有生物體內(nèi)；TCA三羧酸循環(huán)（Tricarboxylic acid cycle）是需氧生物體內(nèi)普遍存在的代謝途徑；PPP磷酸戊糖途徑 （Pentose phosphate pathway）不是機(jī)體產(chǎn)能的方式，但生成具有重要生理功能的NADPH和5-磷酸核糖。文中基于上述模型整合的兩種方法，提出了一種基于基本核心代謝途徑[31]（EMP、TCA、PPP）完善初模型的方法，該方法通過已有核心反應(yīng)列表對(duì)模型數(shù)據(jù)的識(shí)別，可以快速確定模型是否包含完整的基本核心代謝途徑。該方法通用于絕大多數(shù)生物體。核心反應(yīng)識(shí)別流程如圖5所示，基本核心代謝反應(yīng)列表見表2。

圖5 核心反應(yīng)識(shí)別流程Fig.5 Process of Identification of core reaction

運(yùn)用本方法對(duì)上述構(gòu)建出的樹干畢赤酵母代謝網(wǎng)絡(luò)模型進(jìn)行核心反應(yīng)識(shí)別，結(jié)果無核心斷點(diǎn)，每條核心反應(yīng)都可以在模型反應(yīng)列表中找到相同反應(yīng)。進(jìn)一步人工查看發(fā)現(xiàn)無誤，說明了該方法的有效性，同時(shí)也說明上述3種構(gòu)建初模型并整合為相對(duì)精細(xì)模型方法的數(shù)據(jù)完善性。結(jié)果中，基于化合物數(shù)據(jù)庫方法識(shí)別的核心反應(yīng)數(shù)為20條，馬氏距離計(jì)算反應(yīng)式字符頻度直方圖補(bǔ)充識(shí)別了剩余5條反應(yīng)，而這5條反應(yīng)未被第一步識(shí)別的主要原因?yàn)榉磻?yīng)式中出現(xiàn)了H。

3 討論

基于KEGG網(wǎng)頁抓取數(shù)據(jù)在提高了模型構(gòu)建效率的同時(shí)，保證了數(shù)據(jù)的最新性和全面性（包含反應(yīng)式、酶、基因、代謝途徑等），但提取的反應(yīng)無方向性，全部表示為可逆?；贘ava實(shí)現(xiàn)的Uniprot-MetaCyc模型構(gòu)建的方法保證了程序的跨平臺(tái)通用性；將兩個(gè)數(shù)據(jù)庫本地化也降低了遠(yuǎn)程web訪問的時(shí)間。但構(gòu)建出的模型出現(xiàn)一個(gè)基因?qū)?yīng)多個(gè)反應(yīng)的情況，需要進(jìn)一步進(jìn)行核對(duì)篩選，同時(shí)確定同工酶或聚合酶?；谕淳昴Ｐ蜆?gòu)建的代謝網(wǎng)絡(luò)模型的真實(shí)性更高，里面包含信息比較全面，包括代謝反應(yīng)區(qū)室的劃分，代謝途徑的確定，化合物完整信息列表，方向的可逆性等。所以整合上述3種方法構(gòu)建出的初模型，以同源比對(duì)結(jié)果模型為目標(biāo)進(jìn)行整合，當(dāng)前兩種反應(yīng)列表中有相同反應(yīng)時(shí)，以同源菌株模型構(gòu)建的代謝網(wǎng)絡(luò)模型為準(zhǔn)，這樣同時(shí)也確定了大部分反應(yīng)與化合物的附加信息，減少了模型修正的大量工作。

基于化合物數(shù)據(jù)庫與字符頻度直方圖特征兩種方法的整合，彌補(bǔ)了彼此整合過程中的缺憾，提高了反應(yīng)式的匹配率。比如，前者不能判別電荷不平衡引起的相同反應(yīng)，以及因分子式中特殊符號(hào)或分子機(jī)構(gòu)順序不同造成的同種化合物不能識(shí)別；后者通過計(jì)算每個(gè)反應(yīng)式的字符頻度直方圖特征，設(shè)定一定的閾值，快速找到前者不能匹配的相同反應(yīng)。反過來，前者能夠識(shí)別某些只相差幾個(gè)字母或數(shù)字的不同反應(yīng)。比如D-果糖-6-磷酸轉(zhuǎn)化為D-果糖-1，6-二磷酸，即第三條核心反應(yīng) （ATP+DFructose-6-phosphate-＞ ADP+D-Fructose-1，6-bisphosphate+H+），這條反應(yīng)的缺失直接導(dǎo)致糖酵解途徑中斷，進(jìn)而使得模型不能夠生長。但模型中也包含D-果糖-6-磷酸轉(zhuǎn)化為D-果糖-2，6-二磷酸的反應(yīng) （ATP+D-Fructose-6-phosphate-＞ADP+D-Fructose-2，6-bisphosphate+H+）， 由于兩個(gè)反應(yīng)只有一個(gè)數(shù)字的不同，基于反應(yīng)式字符頻度直方圖特征的整合，兩者被識(shí)別為相同反應(yīng)。

表2 基本核心代謝反應(yīng)列表Tabel 2 List of the basic core metabolic reactions

整合后的樹干畢赤酵母模型需要進(jìn)一步精細(xì)化與修正，通過文獻(xiàn)數(shù)據(jù)以及其他網(wǎng)絡(luò)數(shù)據(jù)信息進(jìn)行完善。基于Matlab平臺(tái)的COBRA工具包函數(shù)，可以對(duì)其進(jìn)行反應(yīng)式電荷平衡檢測(cè)、代謝漏洞查找與填補(bǔ)，以及后續(xù)的完整模型分析等，最終完成一個(gè)高質(zhì)量的基因組規(guī)模代謝網(wǎng)絡(luò)模型。

4 結(jié)語

國內(nèi)外雖然對(duì)代謝網(wǎng)絡(luò)自動(dòng)化重構(gòu)也做了多方面的研究，但仍然不能夠?qū)崿F(xiàn)完全自動(dòng)化構(gòu)建代謝網(wǎng)絡(luò)，代謝網(wǎng)絡(luò)構(gòu)建過程中仍然存在一些無法避免的人工操作與修正工作。作者提出的3種自動(dòng)化重構(gòu)代謝網(wǎng)絡(luò)并整合為一個(gè)相對(duì)精細(xì)模型的方法，以及依據(jù)化合物數(shù)據(jù)庫映射和字符頻度直方圖特征整合初模型并識(shí)別核心反應(yīng)的應(yīng)用，能夠在構(gòu)建一個(gè)相對(duì)精細(xì)模型的過程中實(shí)現(xiàn)計(jì)算機(jī)技術(shù)與代謝網(wǎng)絡(luò)構(gòu)建的最大化結(jié)合，減免了構(gòu)建過程中大量的人力與時(shí)間，提高了代謝網(wǎng)絡(luò)構(gòu)建的效率及其精確性。

[1]Veliz-Cuba A，Jarrah A S，Laubenbacher R.Polynomial algebra of discrete models in systems biology[J].Bioinformatics，2010，26（13）：1637-1643.

[2]Van Norman J M，Benfey P N.Arabidopsis thaliana as a model organism in systems biology[J].Wiley Interdisciplinary Reviews：Systems Biology and Medicine，2009，1（3）：372-379.

[3]Reed J L，Vo T D，Schilling C H，et al.An expanded genome-scale model of Escherichia coli K-12 （iJR904 GSM/GPR）[J].Genome Biol，2003，4（9）：R54.

[4]Shlomi T，Eisenberg Y，Sharan R，et al.A genome-scale computational study of the interplay between transcriptional regulation and metabolism[J].Molecular Systems Biology，2007，3（1）：1.

[5]Hyduke D R，Palsson B ?.Towards genome-scale signalling-network reconstructions[J].Nature Reviews Genetics，2010，11（4）：297-307.

[6]Thiele I，Jamshidi N，F(xiàn)leming R M T，et al.Genome-scale reconstruction of Escherichia coli's transcriptional and translational machinery：a knowledge base，its mathe-matical formulation，and its functional characterization[J].PLoS Computational Biology，2009，5（3）：e1000312.

[7]Oberhardt M A，Palsson B O，Papin J A.Applications of genome-scale metabolic reconstructions[J].Mol Syst Biol，2009（5）：320.

[8]Systemsbiology[EB/OL].[2013-11-21].http：//systems biology.ucsd.edu/InSilicoOrganisms/OtherOrganisms

[9]Genomes OnLine Database[EB/OL].[2013-11-21].http：//www.genomesonline.org/cgi-bin/GOLD/bin/gold.cgi.

[10]Thiele I，Palsson B O.A protocol for generating a high-quality genome-scale metabolic reconstruction[J].Nat Protoc，2010，5（1）：93-121.

[11]Henry C S，DeJongh M，Best A A，et al.High-throughput generation，optimization and analysis of genome-scale metabolic models[J].Nature Biotechnology，2010，28（9）：977-982.

[12]DeJongh M，F(xiàn)ormsma K，Boillot P，et al.Toward the automated generation of genome-scale metabolic networks in the SEED[J].BMC Bioinformatics，2007，8（1）：139.

[13]Arakawa K，Yamada Y，Shinoda K，et al.GEM System：Automatic prototyping of cell-wide metabolic pathway models from genomes[J].BMC Bioinformatics，2006，7（1）：168.

[14]Agbogbo F K，Coward-Kelly G.Cellulosic ethanol production using the naturally occurring xylose-fermenting yeast，Pichia stipitis[J].Biotechnology Letters，2008，30（9）：1515-1524.

[15]Jeffries T W，Van Vleet J R H.Pichia stipitis genomics，transcriptomics，and gene clusters[J].FEMS Yeast Research，2009，9（6）：793-807.

[16]Zou W，Zhou M，Liu L，et al.Reconstruction and analysis of the industrial strain Bacillus megaterium WSH002 genome-scale in silico metabolic model[J].Journal of Biotechnology，2013：1.

[17]Dreyfuss J M，Zucker J D，Hood H M，et al.Reconstruction and validation of a genome-scale metabolic model for the filamentous fungus Neurospora crassa using FARM[J].PLoS Computational Biology，2013，9（7）：e1003126.

[18]Kanehisa M，Goto S.KEGG：Kyoto encyclopedia of genes and genomes[J].Nucleic Acids Research，2000，28（1）：27-30.

[19]Caspi R，F(xiàn)oerster H，F(xiàn)ulcher C A，et al.The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases[J].Nucleic Acids Research，2008，36（suppl 1）：D623-D631.

[20]Caspi R，Altman T，Dale J M，et al.The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases[J].Nucleic Acids Research，2010，38（suppl 1）：D473-D479.

[21]Bairoch A，Apweiler R，Wu C H，et al.The universal protein resource（UniProt）[J].Nucleic Acids Research，2005，33（suppl 1）：D154-D159.

[22]Pertsemlidis A，F(xiàn)ondon J W，John W.Having a BLAST with bioinformatics and avoiding BLASTphemy[J].Genome Biol，2001，2（10）：1.

[23]Becker S A，Palsson B ?.Genome-scale reconstruction of the metabolic network in Staphylococcus aureus N315：An initial draft to the two-dimensional annotation[J].BMC Microbiology，2005，5（1）：8.

[24]Mo M L，Palsson B ?，Herrgard M J.Connecting extra-cellular metabolomic measurements to intracellular flux states in yeast[J].BMC Systems Biology，2009，3（1）：37.

[25]Kuepfer L，Sauer U，Blank L M.Metabolic functions of duplicate genes in Saccharomyces cerevisiae[J].Genome Research，2005，15（10）：1421-1430.

[26]Chung B K S，Selvarasu S，Camattari A，et al.Research Genome-scale metabolic reconstruction and in silico analysis of methylotrophic yeast Pichia pastoris for strain improvement[J].Microbial Cell Factories，2010（9）：50

[27]Caspeta L，Shoaie S，Agren R，et al.Genome-scale metabolic reconstructions of Pichia stipitis and Pichia pastoris and in silico evaluation of their potentials[J].BMC Systems Biology，2012，6（1）：24.

[28]Dreyfuss J M，Zucker J D，Hood H M，et al.Reconstru-ction and validation of a genome-scale metabolic model for the filamentous fungus Neurospora crassa using FARM[J].PLoS Computational Biology，2013，9（7）：e1003126.

[29]Bedrick E J.Graphical modelling and the Mahalanobis distance[J].Journal of Applied Statistics，2005，32（9）：959-967.

[30]Leuven E，Sianesi B.PSMATCH2：Stata module to perform full mahalanobis and propensity score matching，common support graphing，and covariate imbalance testing[J].Statistical Software Components，2012：1.

[31]Riemer S A，Rex R，Schomburg D.A metabolite-centric view on flux distributions in genome-scale metabolic models[J].BMC Systems Biology，2013，7（1）：33.