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