基因組規(guī)模代謝網(wǎng)絡(luò)模型自動化修正
2017-11-01 08:12吳曉紅石貴陽
食品與生物技術(shù)學(xué)報 2017年9期
吳曉紅， 薛 衛(wèi)， 張 梁*， 石貴陽
（1.糧食發(fā)酵工藝與技術(shù)國家工程實驗室，江南大學(xué)，江蘇 無錫 214122；2.南京農(nóng)業(yè)大學(xué) 信息科學(xué)技術(shù)學(xué)院 江蘇 南京 210095）
基因組規(guī)模代謝網(wǎng)絡(luò)模型自動化修正
吳曉紅1， 薛 衛(wèi)2， 張 梁*1， 石貴陽1
（1.糧食發(fā)酵工藝與技術(shù)國家工程實驗室，江南大學(xué)，江蘇 無錫 214122；2.南京農(nóng)業(yè)大學(xué) 信息科學(xué)技術(shù)學(xué)院 江蘇 南京 210095）
基于KEGG在線數(shù)據(jù)庫以及6個蛋白質(zhì)區(qū)間預(yù)測數(shù)據(jù)庫，對基因組規(guī)模代謝網(wǎng)絡(luò)模型進(jìn)行了自動化修正。作者提出了蛋白質(zhì)區(qū)間預(yù)測結(jié)果的權(quán)重打分機(jī)制，同時利用圖像處理算法確定可信度高的特異性反應(yīng)。上述修正的研究均在Spathaspora passalidarum NRRL Y-27907基因組規(guī)模代謝網(wǎng)絡(luò)精煉過程中得到運用實施，對于提高模型構(gòu)建效率意義重大。
基因組規(guī)模；代謝網(wǎng)絡(luò)；斷點補(bǔ)齊；圖像處理；區(qū)間預(yù)測
隨著基因組高通量測序數(shù)據(jù)的涌現(xiàn)以及大量的生物學(xué)數(shù)據(jù)的產(chǎn)生，代謝網(wǎng)絡(luò)模型構(gòu)建成為研究生物信息學(xué)的熱點之一。代謝網(wǎng)絡(luò)構(gòu)建是一個耗時費力的過程，因此許多自動化構(gòu)建的工具隨之應(yīng)運而生。通常這些自動化工具側(cè)重關(guān)注代謝網(wǎng)絡(luò)粗模型的構(gòu)建如 metaSHARK[1]和 AUTOGRAPH[2]，其次關(guān)注代謝網(wǎng)絡(luò)模型的模擬過程，如CellNetAnalyzer[3]、OptFlux[4]和 COBRA Toolbox[5]，只有少量的自動化工具是針對代謝網(wǎng)絡(luò)模型的精煉過程。目前能夠提供代謝網(wǎng)絡(luò)模型自動化精煉過程的工具有Model SEED、Pathway Tools、RAVEN 和 SuBliMinaL。
代謝網(wǎng)絡(luò)的模型構(gòu)建包括粗模型的構(gòu)建、模型的精煉、數(shù)學(xué)模型的轉(zhuǎn)換、模型的預(yù)測驗證四個過程。一個高質(zhì)量的代謝網(wǎng)絡(luò)模型，應(yīng)達(dá)到模型模擬結(jié)果和生物實際生長表型一致，否則要不斷的重復(fù)精煉修正過程，直到模擬與表型一致。模型的精煉修正無疑是代謝網(wǎng)絡(luò)模型構(gòu)建過程中最耗時耗力的過程，現(xiàn)有模型精煉工具并不能真正實現(xiàn)真菌代謝網(wǎng)絡(luò)模型精煉過程的自動化。模型的精煉過程必須包括漏洞代謝的填補(bǔ)、反應(yīng)區(qū)間定位等。Model SEED[6]和Pathway Tools[7]只能提供原核生物的代謝網(wǎng)絡(luò)模型的精煉自動化過程，不能提供反應(yīng)區(qū)間的定位。RAVEN[8]和SuBliMinaL[9]是基于Wolf PSORT蛋白質(zhì)區(qū)間預(yù)測數(shù)據(jù)庫實現(xiàn)自動化定位區(qū)間的程序。但是Wolf PSORT[10]只是基于氨基酸組成特征的在線預(yù)測數(shù)據(jù)庫。研究表明，基于氨基酸組成、二肽和物理化學(xué)三種綜合特征的蛋白質(zhì)區(qū)間定位預(yù)測結(jié)果更為準(zhǔn)確[11]。
利用作者所在實驗室自動化構(gòu)建全基因組代謝網(wǎng)絡(luò)模型的程序，自動構(gòu)建了Spathasporapassalidarum NRRL Y-27907全基因組規(guī)模代謝的粗模型。以S.passalidarum NRRL Y-27907的基因組規(guī)模代謝網(wǎng)絡(luò)模型的精煉過程為例，以簡單、面向?qū)ο蟮腏ava語言為基礎(chǔ)，對精煉過程中人工冗雜的斷點補(bǔ)齊的方法進(jìn)行了研究，提出了一種基于KEGG[12]在線數(shù)據(jù)庫自動化填補(bǔ)漏洞反應(yīng)的方法，并利用權(quán)重打分機(jī)制分析，6個真菌蛋白質(zhì)定位數(shù)據(jù)庫預(yù)測S.passalidarum NRRL Y-27907的結(jié)果，在保證模型中反應(yīng)的物種特異性的同時，實現(xiàn)了真菌代謝網(wǎng)絡(luò)模型精煉的自動化。自動化修正的流程見圖1。圖中進(jìn)程g、進(jìn)程n、進(jìn)程o為一個小的流程循環(huán)。進(jìn)程g中判斷反應(yīng)包含斷點，則進(jìn)入進(jìn)程h，查找該反應(yīng)在注釋圖譜中對應(yīng)的坐標(biāo)，并在進(jìn)程i中讀取此坐標(biāo)，在進(jìn)程j中判斷此坐標(biāo)是否為特異性反應(yīng)，如果是，則在進(jìn)程p中記錄該反應(yīng)。如果不是，則在進(jìn)程l中判斷此坐標(biāo)是否為最后一個坐標(biāo)，如果是最后一個坐標(biāo)，則進(jìn)入進(jìn)程n，即進(jìn)入進(jìn)程g、進(jìn)程n、進(jìn)程o該流程循環(huán)。如果不是最后一個坐標(biāo)，則進(jìn)入進(jìn)程m，讀取下一個坐標(biāo)，判斷此坐標(biāo)是否為特異性反應(yīng)，重復(fù)此循環(huán)直至將所有的特異性反應(yīng)都被找出，進(jìn)入進(jìn)程q，進(jìn)行模型修正。在進(jìn)程r中判斷模型中是否已經(jīng)包含此反應(yīng)，若已經(jīng)包含，則回到進(jìn)程n，即進(jìn)入進(jìn)程g、進(jìn)程n、進(jìn)程o該流程循環(huán)，檢查下一條反應(yīng)。若不包含此反應(yīng)，則進(jìn)入進(jìn)程s，將此反應(yīng)加入到模型中。
圖1 自動補(bǔ)齊斷點流程Fig.1 Process of the auto-refinement of gap
1 自動填補(bǔ)網(wǎng)絡(luò)漏洞
采用柴文平[13]等人的方法構(gòu)建了S.passalidarum NRRL Y-27907代謝網(wǎng)絡(luò)粗模型。構(gòu)建的代謝網(wǎng)絡(luò)粗模型需要進(jìn)一步精細(xì)化與修正，最終完成一個高質(zhì)量的基因組規(guī)模代謝網(wǎng)絡(luò)模型。
1.1 代謝網(wǎng)絡(luò)漏洞查找
模型導(dǎo)入到裝有COBRA工具包和GLPK線性規(guī)劃器的Matlab中，將模型轉(zhuǎn)化為計算機(jī)可讀的格式 （SBML）才能進(jìn)行代謝網(wǎng)絡(luò)漏洞查找。通過xls2model程序?qū)⒛Ｐ虴xcel表讀取為計量學(xué)S矩陣。S矩陣（828×984）表示該模型由828個代謝物和984個反應(yīng)組成。同時通過GapFind程序完成代謝漏洞的查找，其中上游漏洞代謝物有為44個，下游漏洞代謝物有128個。
1.2 基于KEGG網(wǎng)絡(luò)爬蟲反應(yīng)
KEGG是代謝網(wǎng)絡(luò)構(gòu)建常用數(shù)據(jù)庫，含有多個在線子數(shù)據(jù)庫，其中REACTION數(shù)據(jù)庫包含迄今為止發(fā)現(xiàn)的所有生化反應(yīng)。各個子數(shù)據(jù)庫的網(wǎng)頁數(shù)據(jù)格式比較統(tǒng)一明確，方便人們進(jìn)行遠(yuǎn)程服務(wù)器訪問。但是，KEGG數(shù)據(jù)庫更新頻繁，各個子數(shù)據(jù)庫不能夠免費下載，需要付費使用。而在基因組代謝網(wǎng)絡(luò)斷點補(bǔ)齊過程中，因為數(shù)據(jù)信息量浩大，頻繁訪問遠(yuǎn)程服務(wù)器比較耗時耗力。因此，實現(xiàn)一種批量在線獲取并存取數(shù)據(jù)的方法意義重大。
1.2.1 方法概述 利用超文本轉(zhuǎn)移協(xié)議和Java控件HttpClient相結(jié)合，實現(xiàn)對網(wǎng)頁中特定信息的抓取KEGG提供物種特異性基因組信息以及所有反應(yīng)式信息查詢網(wǎng)頁，通過一定的URL（Uniform Resource Locator，統(tǒng)一資源定位符）格式地址發(fā)送HTTP請求并獲取網(wǎng)頁中的基因信息。在漏洞填補(bǔ)的過程中需要訪問大量不同的網(wǎng)絡(luò)資源，獲取相關(guān)的基因信息，由于數(shù)據(jù)量較大且人工操作比較繁瑣，這里利用Java控件HttpClient實現(xiàn)爬蟲技術(shù)，抓去符合特定條件的網(wǎng)絡(luò)資源。HttpClient是Apache Jakarta Common下的子項目，可以用來提供高效的、最新的、功能豐富的支持HTTP協(xié)議的客戶端編程工具包，并且它支持HTTP協(xié)議最新的版本和建議。利用HttpClient訪問具體的URL地址，獲取服務(wù)器端返回的獲取html內(nèi)容，html內(nèi)容由標(biāo)題、js代碼、正文、相關(guān)鏈接、聲明等區(qū)域組成，而有用信息只出現(xiàn)在正文中的各種html標(biāo)簽標(biāo)記內(nèi)，分析html標(biāo)簽并獲取特定的網(wǎng)頁信息。
1.2.2 漏洞填補(bǔ)算法實現(xiàn)
1）獲取注釋圖譜：提交物種基因組蛋白質(zhì)序列至KAAS自動注釋服務(wù)器，獲取注釋信息，下載html和text格式。
2）查找包含斷點的注釋圖譜：根據(jù)Matlab軟件中GapFind程序返回的漏洞代謝物列表，在代謝網(wǎng)絡(luò)模型Excel格式中確定代謝物的反應(yīng)途徑，依據(jù)KASS注釋返回的途徑圖譜找到包含漏洞代謝物的所有反應(yīng)。
注釋返回的KEGG代謝途徑為包含糖代謝等在內(nèi)的110個途徑。查找包含斷點的代謝圖譜的流程見圖2。具體思路和偽代碼步驟如下：
A：獲取斷點化合物所對應(yīng)的Subsystem信息，記為sub。
B：向注釋查詢網(wǎng)頁URL地址發(fā)送HTTP請求。
C：如果服務(wù)器端響應(yīng)代碼為HTTPStatus.SC_OK則正常響應(yīng)，否則繼續(xù)請求，獲取html正文內(nèi)容。
D：分析html內(nèi)容，設(shè)i為行號，由第一行開始遍歷標(biāo)簽對中的每一行，
For i from 1 to n
｛
if（該行中第二個標(biāo)簽中的內(nèi)容與sub相等）
｛
提取對應(yīng)的第一個標(biāo)簽中的內(nèi)容，記為KO；
｝
else
忽略該行，遍歷下一行；
｝
E：根據(jù)D中的KO號得到滿足條件圖譜的URL地址，向URL地址發(fā)送HTTP請求得到服務(wù)器端響應(yīng)的網(wǎng)頁圖片記為T1，T1即為整個網(wǎng)絡(luò)結(jié)構(gòu)圖，其中綠色酶號表示包含斷點的特異性反應(yīng)。
F：點擊T1左上角途徑方框，進(jìn)去包含所有反應(yīng)頁面page1，網(wǎng)頁中每一個EC號對應(yīng)圖譜中的一個具體反應(yīng)，它的URL地址指向具體的反應(yīng)方程式。
G：獲取page1中所有EC號對應(yīng)的反應(yīng)，設(shè)ec_num為每一個EC號，從第一個開始
EC_K_Break.txt保存包含斷點化合物的EC，K號的信息。3）查找EC_K_Break.txt中每個K對應(yīng)的坐標(biāo)根據(jù)K號獲取其在T1中對應(yīng)的坐標(biāo)，判斷特異性反應(yīng)。
圖2 斷點代謝途徑定位Fig.2 Orientation of gap metabolic pathway
1.3 判斷特異性反應(yīng)
KEGG所有的反應(yīng)都包含在通路數(shù)據(jù)庫（PATHWAY database）中，PATHWAY 圖譜上有顏色標(biāo)記的酶號是指這個物種特定的基因或酶，只有有顏色標(biāo)記的酶號表示的反應(yīng)才是具有該物種特異性的反應(yīng)，也才能添加到代謝網(wǎng)絡(luò)模型中。在代謝網(wǎng)絡(luò)模型中添加非特異性的反應(yīng)會改變整個代謝途徑和代謝物流量，進(jìn)而使模型模擬的結(jié)果偏離實驗數(shù)據(jù)，影響模型的準(zhǔn)確性和可信度。
構(gòu)建代謝網(wǎng)絡(luò)模型需要提取代謝途徑中的特異性反應(yīng)，圖中特異性反應(yīng)對應(yīng)的酶號所在的方形框有顏色標(biāo)記。因此通過網(wǎng)絡(luò)爬蟲技術(shù)獲得方形框的位置列表，定位到某酶號所在的方形框后需要選取框內(nèi)的像素點，讀取其顏色值，如果顏色分量RGB均為0或255，則沒有顏色標(biāo)記，反之則有。代謝網(wǎng)絡(luò)特異性反應(yīng)獲取流程見圖3。
基本思路為：
根據(jù)得到的position坐標(biāo)讀取T1對應(yīng)點的RGB色彩值。
Picture（Key：酶號；Value：代謝網(wǎng)絡(luò)圖中所有方形框的坐標(biāo)向量集｛V1，V2，……，Vn｝）
For i from 1 to n
{
If（某酶號所在的方形框）
｛
沿方形框的長邊內(nèi)側(cè)逐一選取像素點，讀取其顏色值；
If顏色分量RGB均為0或255 then沒有顏色標(biāo)記
else有顏色標(biāo)記；
If有顏色標(biāo)記then該酶號對應(yīng)的是特異性反應(yīng)
do將反應(yīng)加入菌的代謝網(wǎng)絡(luò)模型中；
else舍棄該酶號對應(yīng)的反應(yīng)。
｝｝
反應(yīng)式漏洞填補(bǔ)
遍歷new_rec.TXT中每一個反應(yīng)，查看模型中是否存在，存在則不處理，否則添加。
A：讀取new_rec.TXT中每行反應(yīng)記為new_rec，i為行號
For i from 1 to n
｛
if（模型中不包含 new_rec）
｛
將new_rec添加到模型中；
｝
else
忽略該反應(yīng)，查找下一條反應(yīng)；｝
圖3 特異性反應(yīng)獲取流程Fig.3 Process of getting the pecificreaction
2 獲取反應(yīng)區(qū)間定位
細(xì)胞是生命活動的基本單位，它由執(zhí)行不同機(jī)體功能的稱為亞細(xì)胞的各部分組成，如細(xì)胞膜、細(xì)胞核、線粒體、高爾基體、內(nèi)質(zhì)網(wǎng)等。亞細(xì)胞功能是由位于其中的蛋白質(zhì)執(zhí)行的，蛋白質(zhì)所在的亞細(xì)胞稱為蛋白質(zhì)的亞細(xì)胞位置[14]。蛋白質(zhì)必須轉(zhuǎn)運到其應(yīng)在的亞細(xì)胞位置上才能正確行使其功能，否則就會出現(xiàn)機(jī)體功能紊亂，正確合理的蛋白區(qū)間定位是高質(zhì)量模型構(gòu)建的基礎(chǔ)，見表1。
表1 真菌蛋白質(zhì)亞細(xì)胞預(yù)測數(shù)據(jù)庫Table 1 Database for subcellular localization of fungal proteins
確定一條蛋白質(zhì)的亞細(xì)胞位置稱為蛋白質(zhì)亞細(xì)胞定位[15]。蛋白質(zhì)亞細(xì)胞定位的傳統(tǒng)方法是通過生物化學(xué)實驗，如射線晶體衍射電子顯微鏡核磁共振等方法進(jìn)行測定[16]。實驗方法精確度高，但費時耗力代價昂貴，而且對難于結(jié)晶的蛋白質(zhì)來說，實驗方法不再有效。借助于先進(jìn)高效的計算機(jī)自動化數(shù)據(jù)處理技術(shù)，出現(xiàn)了一些蛋白質(zhì)定位預(yù)測網(wǎng)站。結(jié)合Spathasporapassalidarum NRRL Y-27907的生理生化性質(zhì)和蛋白質(zhì)特征提取方法、算法和準(zhǔn)確性等，選取了6個真菌生物蛋白質(zhì)區(qū)間預(yù)測網(wǎng)站，自動化提取分析網(wǎng)站的預(yù)測結(jié)果，在權(quán)重打分機(jī)制的基礎(chǔ)上得到最佳的蛋白質(zhì)定位區(qū)間。這6個網(wǎng)站是基于蛋白質(zhì)的氨基酸組成、偽氨基酸組成、二肽、生物化學(xué)特征或是四種特征的綜合。
2.1 區(qū)間定位算法實現(xiàn)
A：對每條反應(yīng)獲取對應(yīng)的KO號。
B：將A中的KO號在KASS注釋結(jié)果中查找基因號，并在本地下載Spathasporapassalidarum NRRL Y-27907蛋白質(zhì)序列庫提取其對應(yīng)的蛋白質(zhì)序列。
C：將蛋白質(zhì)序列提交到對應(yīng)網(wǎng)站的表單中，獲取返回的定位信息。
D：獲取定位區(qū)間的信息并填入反應(yīng)式中。
在獲取具體反應(yīng)的區(qū)間信息過程中，需要將反應(yīng)所對應(yīng)的蛋白質(zhì)序列提交到網(wǎng)頁的表單中，提交后返回具體的區(qū)間定位信息，此處會遇到兩個問題：1）表單提交過程中不支持大量蛋白質(zhì)序列自動提交。由于模型中蛋白質(zhì)序列數(shù)量較大，在有的網(wǎng)站中獲取定位信息時不支持大量序列的一次性提交而只能分別提交單個序列獲得定位信息，在提交過程中任務(wù)量大且人工耗費時間長。2）大量蛋白質(zhì)序列提交耗費時間長，在網(wǎng)站中提交多個序列后等待服務(wù)器端反饋的定位信息耗費時間太長，甚至?xí)l(fā)生無響應(yīng)等問題，見圖4。
圖4 反應(yīng)亞細(xì)胞定位流程Fig.4 Process of subcellular localization
HttpClient支持訪問特定的URL地址，獲取服務(wù)器端返回的html信息，并且能夠分析html中form表單中的信息，實現(xiàn)內(nèi)容的自動提交。由于涉及到的定位頁面所有的表單提交方式都是POST提交，利用HttpClient中的PostMethod方法實現(xiàn)post提交。表單中的元素賦值過程：獲取表單中需要賦值的元素標(biāo)簽，以蛋白質(zhì)序列元素賦值標(biāo)簽為例，標(biāo)簽為