鄧世果, 吳干華,2, 楊會(huì)杰
(1上海理工大學(xué) 管理學(xué)院,上海 200090;2華南師范大學(xué) 南海學(xué)院,佛山 528225)
系統(tǒng)生物學(xué)的一個(gè)重要任務(wù)是認(rèn)識(shí)細(xì)胞內(nèi)調(diào)控和代謝關(guān)系.代謝網(wǎng)絡(luò)依據(jù)基因組數(shù)據(jù)將代謝系統(tǒng)表示為一個(gè)圖,為代謝組學(xué)數(shù)據(jù)的分析提供可視化的研究平臺(tái).隨著基因組測序技術(shù)的迅速發(fā)展,到目前為止,包括細(xì)菌、古生菌、真核生物在內(nèi)的數(shù)百種生物全基因組序列先后測序完成.根據(jù)基因組注釋信息可以按功能的不同將基因分成不同的組,其中較大的組是可指導(dǎo)合成酶的基因組.酶可催化細(xì)胞內(nèi)代謝反應(yīng),依據(jù)代謝反應(yīng)可構(gòu)造代謝網(wǎng)絡(luò).人們構(gòu)建了專門的代謝反應(yīng)數(shù)據(jù)庫,如KEGG、BioCyc、PUMA2等,這些數(shù)據(jù)可以用于構(gòu)建代謝網(wǎng)絡(luò),并且可以對(duì)代謝網(wǎng)絡(luò)作進(jìn)一步分析.作為復(fù)雜非線性動(dòng)態(tài)調(diào)控系統(tǒng),代謝網(wǎng)絡(luò)的研究越來越受到重視[1].
基因組范圍的代謝網(wǎng)絡(luò)包含的代謝物很多,傳統(tǒng)的方法(如代謝控制)不太實(shí)用.Ma等[2-3]首次提出具有可操作性的圖論方法用以構(gòu)造和分析代謝網(wǎng)絡(luò),得到普遍認(rèn)可,成為當(dāng)前對(duì)代謝系統(tǒng)的整體結(jié)構(gòu)進(jìn)行分析的主要手段.構(gòu)造的代謝網(wǎng)絡(luò)有3類:用節(jié)點(diǎn)表示代謝物,用節(jié)點(diǎn)間的連線表示代謝反應(yīng),構(gòu)成代謝物網(wǎng)絡(luò);以代謝反應(yīng)為節(jié)點(diǎn),以節(jié)點(diǎn)之間連線表示代謝物,就是反應(yīng)網(wǎng)絡(luò)[4-5];以反應(yīng)所需的酶為節(jié)點(diǎn),以節(jié)點(diǎn)間連線表示代謝物,得到酶網(wǎng)絡(luò)[6-8].以上構(gòu)建代謝網(wǎng)絡(luò)方法的基本步驟是:a.從KEGG數(shù)據(jù)庫取得某一個(gè)特定物種細(xì)胞內(nèi)全基因組注釋信息,得到所有可能的酶;b.根據(jù)酶催化的代謝反應(yīng)數(shù)據(jù)庫,得到選取的物種細(xì)胞內(nèi)所有可能的化學(xué)反應(yīng),確定出這些化學(xué)反應(yīng)的底物和產(chǎn)物;c.確定出每個(gè)化學(xué)反應(yīng)的反應(yīng)方向;d.將每個(gè)化學(xué)反應(yīng)的反應(yīng)物和底物用有向邊連接起來,得到代謝物網(wǎng)絡(luò).其它類代謝網(wǎng)絡(luò)可以采取相似的辦法得到.
上述步驟中的關(guān)鍵問題和需要深入考慮的問題包括:a.網(wǎng)絡(luò)節(jié)點(diǎn)的確定.因?yàn)楹芏嗌磻?yīng)涉及到水、氧氣、三磷酸腺苷(ATP)、二磷酸腺苷(ADP)等物質(zhì),如果在代謝網(wǎng)絡(luò)中保留這些物質(zhì)的話,大部分的節(jié)點(diǎn)之間通過這類節(jié)點(diǎn)就彼此相連了,整個(gè)網(wǎng)絡(luò)將是以這些節(jié)點(diǎn)為核心的致密結(jié)構(gòu),核心致密的網(wǎng)絡(luò)結(jié)構(gòu)掩蓋了人們感興趣的網(wǎng)絡(luò)結(jié)構(gòu).同時(shí),ATP和ADP等只是在反應(yīng)中提供能量,并不參與物質(zhì)的合成和分解.因此,這類物質(zhì)應(yīng)該去掉,以凸顯代謝網(wǎng)絡(luò)的非平凡結(jié)構(gòu).有些文獻(xiàn)依據(jù)經(jīng)驗(yàn),簡單地將節(jié)點(diǎn)度(也就是與其它節(jié)點(diǎn)聯(lián)結(jié)的數(shù)目)大于10的節(jié)點(diǎn)刪除掉;有些文獻(xiàn)則將節(jié)點(diǎn)按照度排序,將排在前20位的節(jié)點(diǎn)刪除掉;有些文獻(xiàn)比較保守,保留了這些度大的節(jié)點(diǎn),同時(shí)也列出了這些特殊節(jié)點(diǎn)作為分析代謝網(wǎng)絡(luò)的參考[9].b.生化反應(yīng)方向的確定.Ma等[2]提出的方案中,系統(tǒng)地列舉了多種情況,如消耗氧氣、產(chǎn)生二氧化碳等,認(rèn)為這些化學(xué)反應(yīng)是不可逆的.并且對(duì)整個(gè)化學(xué)反應(yīng)庫進(jìn)行了系統(tǒng)的鑒別.但是,這一方法本身需要經(jīng)驗(yàn)的累計(jì),不可避免地有人為因素起作用.特別是他們的工作是在本世紀(jì)初完成的,而近10年來代謝數(shù)據(jù)庫不斷地進(jìn)行更新,需要有更具可操作性的代謝網(wǎng)絡(luò)的構(gòu)造方法.本文在對(duì)前人文獻(xiàn)總結(jié)的基礎(chǔ)上,以構(gòu)建魚腥藻代謝網(wǎng)絡(luò)為例,試圖對(duì)這兩個(gè)關(guān)鍵問題提出一些可行的解決方案.
通過基因組注釋信息可以識(shí)別出能指導(dǎo)合成催化代謝反應(yīng)的酶的基因.到目前為止出現(xiàn)了多種用于預(yù)測物種特異的酶基因、酶、以及酶催化反應(yīng)的方法,由此產(chǎn)生了許多優(yōu)秀的代謝數(shù)據(jù)庫,如表1所示.這些數(shù)據(jù)庫是代謝網(wǎng)絡(luò)重建的必要資源[10-12].
表1 常用生物數(shù)據(jù)庫Tab.1 Common biological databases
在圖論中,圖表示元素與元素之間的二元關(guān)系,其中,元素表示為圖的頂點(diǎn),元素之間的關(guān)系表示為頂點(diǎn)之間的連線.1個(gè)無向圖G= (V,E),由頂點(diǎn)集合V和邊集合E構(gòu)成,每條邊代表1個(gè)頂點(diǎn)對(duì)(u,v)間的無方向連線.1個(gè)有向圖D= (V,A),由頂點(diǎn)集合V和弧集合A構(gòu)成,其中,每條弧代表1個(gè)頂點(diǎn)對(duì)(u,v)間從u到v的有向邊.如果忽略其中所有弧的方向,則一個(gè)有向圖就成為無向圖[13].復(fù)雜系統(tǒng)諸多元素作為節(jié)點(diǎn),元素之間的關(guān)系作為邊,構(gòu)造出來的圖就是復(fù)雜系統(tǒng)的復(fù)雜網(wǎng)絡(luò)描述,如描寫代謝系統(tǒng)結(jié)構(gòu)的代謝網(wǎng)絡(luò).
3.1.1 反應(yīng)數(shù)據(jù)
KEGG數(shù)據(jù)庫中文件reaction.list包含迄今發(fā)現(xiàn)的所有生化反應(yīng).每個(gè)反應(yīng)都有各自的編號(hào),以R開頭后跟5位數(shù)字;每個(gè)化合物也都有各自的編號(hào),以C開頭后跟5位數(shù)字.如化學(xué)反應(yīng)R00480:
L-Aspartate+ATP=4-Phospho-L-aspartate+ADP
在reaction.list文件中表示為R00480:
C00049+C00002→C03082+C00008
生物體內(nèi)的很多代謝反應(yīng)都是有方向的,不可逆性是代謝反應(yīng)的一個(gè)本質(zhì)特點(diǎn).而KEGG并沒有給出這些信息,Ma和Zeng整理并提出了11種不可逆反應(yīng)[2]:
a.有氧生成的反應(yīng);
b.大多數(shù)有二氧化碳生成的反應(yīng);
c.大多數(shù)有氨氣生成的反應(yīng);
d.大多數(shù)有磷酸鹽生成的反應(yīng);
e.由S-Adenosyl-L-methionine生成 S-Adenosyl-L-h(huán)omocystine,提供一個(gè)甲基的反應(yīng);
f.有四氫葉酸生成,轉(zhuǎn)移一個(gè)碳原的反應(yīng);
g.大多數(shù)消耗ATP,且沒有其它高能代謝物(GTP、CAP等)參與的反應(yīng);
h.消耗UDP-sugar,轉(zhuǎn)移糖的反應(yīng);
i.消耗 CDP-diacylglycerol,轉(zhuǎn)移磷脂酰基的反應(yīng);
j.類似于PAPS+A=PAP+B的反應(yīng);
k.幾種水解反應(yīng).
然后Ma和Zeng又對(duì)所有的生化反應(yīng)進(jìn)行了判別,給出化學(xué)反應(yīng)的方向.
實(shí)際上,文獻(xiàn)[2]中對(duì)代謝途徑有著更加全面深入的研究,對(duì)代謝途徑中發(fā)生的化學(xué)反應(yīng)方向等有著明確的結(jié)論,這些信息保存在KEGG數(shù)據(jù)庫的reaction_mapformula.lst文件中.這一文件收集了文獻(xiàn)中可以找到的生化反應(yīng)的方向、化學(xué)反應(yīng)發(fā)生所在的代謝途徑等信息.因此,本文采取一個(gè)新的解決方案,即可以采取一個(gè)組合的策略確定生化反應(yīng)的方向:a.在reaction_mapformula.lst里邊能查到的化學(xué)反應(yīng),其反應(yīng)方向采用這里邊已經(jīng)明確的反應(yīng)方向;b.剩余的化學(xué)反應(yīng),能采用Ma等提出的11條規(guī)則判斷反應(yīng)方向;c.前兩個(gè)策略不能判斷的其它化學(xué)反應(yīng),設(shè)為雙向的[2,14-15].
3.1.2 反應(yīng)與酶的關(guān)系
從KEGG數(shù)據(jù)庫提取酶與反應(yīng)的有關(guān)信息,見“reaction”文件,其中,給出了每個(gè)反應(yīng)以及催化該反應(yīng)所需的酶的信息,酶都有各自的EC編號(hào),以EC開頭,后跟4個(gè)整數(shù),整數(shù)之間用點(diǎn)隔開,4個(gè)整數(shù)依次對(duì)酶的功能類別進(jìn)行詮釋,越后面的整數(shù)對(duì)酶的功能描述得越詳盡.如R00480:2.7.2.4,表示EC2.7.2.4催化反應(yīng)R00480.
3.1.3 酶與基因的對(duì)應(yīng)關(guān)系
KEGG數(shù)據(jù)庫中含有enzyme文件,enzyme文件對(duì)酶的相關(guān)信息描述得非常詳盡,包括名稱、參與生化反應(yīng)過程和方式、所在的代謝路徑、不同物種中對(duì)應(yīng)的基因名稱、數(shù)據(jù)來源等.可由該數(shù)據(jù)得到某一特定物種里邊存在的所有的酶.
連接初步代謝反應(yīng)網(wǎng)絡(luò)圖的具體做法為:對(duì)每個(gè)反應(yīng),以代謝物和代謝底物為節(jié)點(diǎn),以反應(yīng)為有向邊,從代謝物指向代謝底物.如化學(xué)反應(yīng)R00480:
C00049+C00002→C03082+C00008
如果不可逆,從該反應(yīng)可獲得4條有向邊:
如果可逆,則可將反應(yīng)拆成方向相反的兩個(gè)反應(yīng),方法同上.把所有反應(yīng)都轉(zhuǎn)化一遍后就構(gòu)建完成了一個(gè)初步代謝圖,以糖酵解過程為例,如圖1所示(見下頁).
在圖1(b)中可以發(fā)現(xiàn)藍(lán)色框中的代謝物(CO2、H2O2等)的度比較大,稱為通用代謝物(currency metabolite)[4].然而這些通用代謝物一般是電子轉(zhuǎn)移或某些功能基團(tuán)(磷酸基、氨基、一碳單位、甲基等)轉(zhuǎn)移的攜帶者,只是協(xié)助代謝底物生成代謝產(chǎn)物,并未參與代謝產(chǎn)物的合成.包含有通用代謝物的代謝網(wǎng)絡(luò)表現(xiàn)出與真實(shí)生物意義不符的結(jié)論,使得代謝網(wǎng)絡(luò)的平均最短路徑縮短[5,8].如從草酸到氧氣只需要兩步反應(yīng)(見下頁圖1),顯然這在細(xì)胞內(nèi)是無法完成的,與細(xì)胞內(nèi)的真實(shí)生化反應(yīng)不符.
因此,為了使細(xì)胞中主要化合物的轉(zhuǎn)化表現(xiàn)得更明顯,在代謝網(wǎng)絡(luò)中確切地顯示生化反應(yīng)的步驟,人們通常將這些通用代謝物及一些小分子化合物(如水、氨氣、二氧化碳、氧等)從代謝網(wǎng)絡(luò)中移出.現(xiàn)介紹一種對(duì)此問題的處理方法.
文獻(xiàn)[2]在KEGG代謝反應(yīng)數(shù)據(jù)庫的基礎(chǔ)上進(jìn)行手工修正補(bǔ)充后得到一個(gè)新的數(shù)據(jù)庫,去掉了每個(gè)反應(yīng)中的通用代謝物及小分子化合物,并明確地給出了每個(gè)反應(yīng)的可逆性信息.在這個(gè)數(shù)據(jù)庫中,他們并不是統(tǒng)一地將通用代謝物都去掉,而是根據(jù)每個(gè)反應(yīng)來確定其中的通用代謝物.例如,在許多反應(yīng)中,谷胺酸(glutamate,GLU)和α-酮戊二酸(2-oxoglutarate,AKG)都是用于轉(zhuǎn)移氨基的通用代謝物,然而在反應(yīng)中AKG參與了合成GLU,因此
AKG+NH8+NADPH→GLU+NADP++H2O+AKG在此反應(yīng)中它們都是主要代謝物,它們之間的連接應(yīng)保留[9].由于該數(shù)據(jù)庫幾乎全部用手工構(gòu)建,因此,質(zhì)量有所保證,已被許多研究者采用.但是,他們的處理依據(jù)經(jīng)驗(yàn),帶有人為性質(zhì),并且也已經(jīng)有近10年的歷史,面臨數(shù)據(jù)更新問題.
圖1 糖酵解過程Fig.1 Glycolysis process
因此,本文提出另外一個(gè)可行的辦法.在KEGG數(shù)據(jù)庫的LIGAND部分包含的文件reaction_mapformula.lst中給出了代謝路徑中的化學(xué)反應(yīng),即包含了每個(gè)反應(yīng)的方向以及主要代謝物,而省略了如ATP、NADH這一類的輔助因子(cofactor).確定這些主要代謝物采用的標(biāo)準(zhǔn)和文獻(xiàn)[2]使用的標(biāo)準(zhǔn)是相同的.上述采用的標(biāo)準(zhǔn)只是針對(duì)代謝路徑中存在的化學(xué)反應(yīng)進(jìn)行了處理,而不是所有的化學(xué)反應(yīng)都處理過.因此,本文中采用的策略:a.在reaction_mapformula.lst中出現(xiàn)的化學(xué)反應(yīng),采用該文件中的信息;b.對(duì)比reaction_mapformula.lst與reaction中相同的化學(xué)反應(yīng),確定出所有可能的通用代謝物.將這些通用代謝物從a不能判斷的化學(xué)反應(yīng)中剔除掉.通過以上方法修正后,之前的R00480:
就從4條邊減少為只有1條邊:C00049→C03082這是因?yàn)镃0002和C0008均為通用代謝物,這2個(gè)節(jié)點(diǎn)和與之相連的邊都要去掉.這樣就可以得到一個(gè)完整的代謝物網(wǎng)絡(luò),可在此基礎(chǔ)上開展進(jìn)一步的研究,如網(wǎng)絡(luò)結(jié)構(gòu)、動(dòng)力學(xué)、代謝功能之間的關(guān)系的討論等.
圖2為將代謝物網(wǎng)絡(luò)圖轉(zhuǎn)化為反應(yīng)圖的示意圖.反應(yīng)圖就是以反應(yīng)為節(jié)點(diǎn),如果兩個(gè)反應(yīng)(反應(yīng)A和反應(yīng)B)共用一個(gè)或多個(gè)代謝物或代謝底物,且反應(yīng)A得到的產(chǎn)物是反應(yīng)B的底物,就將這兩個(gè)反應(yīng)連起來,由A指向B,這樣就構(gòu)成了一個(gè)有向的反應(yīng)圖,如圖2(c)所示.相對(duì)代謝圖而言,反應(yīng)圖離酶圖關(guān)系更近了些.具體細(xì)節(jié)和需要注意的問題與代謝圖類似.
圖2 代謝反應(yīng)網(wǎng)絡(luò)圖Fig.2 Metabolic reaction network
酶圖與反應(yīng)圖也并非一一對(duì)應(yīng),因?yàn)椋环N酶可以催化不同的反應(yīng),一種反應(yīng)也可能需要多種酶共同參與.酶圖以酶為節(jié)點(diǎn),以酶A和酶B為例,如果酶A參與的反應(yīng)產(chǎn)生的產(chǎn)物正好是酶B參與的反應(yīng)的底物,則就連接酶A與酶B,且由酶A指向酶B.
通過以上代謝網(wǎng)絡(luò)的構(gòu)造方法,構(gòu)造出魚腥藻的代謝網(wǎng)絡(luò),并對(duì)該網(wǎng)絡(luò)進(jìn)行了社團(tuán)劃分,如圖3所示,當(dāng)分為23個(gè)社團(tuán)時(shí),聚類系數(shù)Q值最大,此時(shí)Q=0.799.
為了更清楚地看清網(wǎng)絡(luò)結(jié)構(gòu),圖3沒有顯示度是1的節(jié)點(diǎn),也刪除掉了2個(gè)孤立的節(jié)點(diǎn).不同顏色表示不同的社團(tuán)結(jié)構(gòu).
本文計(jì)算了網(wǎng)絡(luò)度分布(圖4)和平均聚類系數(shù).由圖4可知,該網(wǎng)絡(luò)的度分布De服從冪律分布,冪律指數(shù)α=-3.0,其聚集系數(shù)C(k)與節(jié)點(diǎn)度k的關(guān)系如圖5所示,因此,不滿足,C(k)~k∝k-α,即該網(wǎng)絡(luò)不存在層次結(jié)構(gòu).經(jīng)計(jì)算,該代謝網(wǎng)絡(luò)的網(wǎng)絡(luò)直徑為D=21,平均最短路經(jīng)約為8.從圖6(見下頁)可知,最可幾分布Me發(fā)生在最短路徑d=7的位置,最可幾概率為0.12.以上結(jié)論與一般生物代謝網(wǎng)絡(luò)的性質(zhì)大體上一致,因此,該重建代謝網(wǎng)絡(luò)的構(gòu)造方法是可行的.
圖3 魚腥藻代謝圖Fig.3 Topological structure for anabaena metabolic network
圖4 魚腥藻代謝網(wǎng)絡(luò)的度分布Fig.4 Degree distribution of anabaena metabolic network
圖5 平均集聚系數(shù)與度的關(guān)系圖Fig.5 Relationship of average clustering coefficient versus degree
圖6 最短路經(jīng)分布Fig.6 Distribution function for shortest paths
代謝網(wǎng)絡(luò)重建是系統(tǒng)生物學(xué)的基本任務(wù),是采用復(fù)雜網(wǎng)絡(luò)觀點(diǎn)分析代謝數(shù)據(jù)的基礎(chǔ).隨著數(shù)據(jù)的日益增多和不斷更新,設(shè)計(jì)合理可行的代謝網(wǎng)絡(luò)重建方法成為迫切需要解決的課題.
本文在總結(jié)前人工作基礎(chǔ)上,就構(gòu)建代謝網(wǎng)絡(luò)步驟中的兩個(gè)核心問題,提出了可行的解決方法.依據(jù)KEGG中文件reaction_mapformula.lst提供的化學(xué)反應(yīng)的方向信息,提出確定代謝網(wǎng)絡(luò)中邊的方向的規(guī)則;并且與reaction文件中的化學(xué)反應(yīng)比較,提出確定通用代謝物的規(guī)則,用于確定代謝物網(wǎng)絡(luò)里的節(jié)點(diǎn).以魚腥藻為例,對(duì)構(gòu)造出來的代謝網(wǎng)絡(luò)進(jìn)行了討論.
代謝網(wǎng)絡(luò)的重建是一個(gè)復(fù)雜的過程,也因研究目的的不同而不同.本文討論的方法適合于復(fù)雜網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)及其與代謝功能相關(guān)性分析.當(dāng)討論代謝網(wǎng)絡(luò)動(dòng)力學(xué),特別是代謝網(wǎng)絡(luò)上的流問題的時(shí)候,則需要更加準(zhǔn)確的網(wǎng)絡(luò)結(jié)構(gòu),必須考慮到化學(xué)反應(yīng)在細(xì)胞內(nèi)發(fā)生的位置信息、特定環(huán)境和細(xì)胞生命周期中基因表達(dá)量等.因此,高精度的代謝網(wǎng)絡(luò)涉及脫氧核糖核酸DNA測序和注釋、化學(xué)反應(yīng)位置確定、基因表達(dá)、化學(xué)反應(yīng)方向、通用代謝物確定等多個(gè)復(fù)雜環(huán)節(jié),是一個(gè)復(fù)雜的系統(tǒng)工程.本文的構(gòu)造方法也可用于構(gòu)造人際關(guān)系等網(wǎng)絡(luò).
[1]Kanehisa M.Post-genome informatics[M]:Oxford:Oxford University Press,2000.
[2]Ma H W,Zeng A P.Reconstruction of metabolic networks from genome data and analysis of their global structure for various organisms[J].Bioinformatics,2003,19(2):270-277.
[3]Ma H W,Zeng A P.The connectivity structure,giant strong component and centrality of metabolic networks[J].Bioinformatics,2003,19(11):1423-1430.
[4]Wagner A,F(xiàn)ell D A.The small world inside large metabolic networks[J].Proc R Soc Lond B,2001,268(1478):1803-1810.
[5]Light S,Kraulis P,Elofsson A.Preferential attachment in the evolution of metabolic networks[J].BMC Bioinformatics,2005,6(1471):159-169.
[6]Horne A B,Hodgman T C,Spence H D,et al.Constructing an enzyme-centric view of metabolism[J].Bioinformatics,2004,20(13):2050-2055.
[7]Heymans M,Singh A K.Deriving phylogenetic trees from the similarity analysis of metabolic pathways[J].Bioinformatics,2003,19(suppl_1):138-146.
[8]Light S,Kraulis P.Network analysis of metabolic enzyme evolution in Escherichia coli[J].BMC Bioinformatics,2004,5(1):15-20.
[9]Huss M, Holme P.Currency and commodity metabolites:their identification and relation to the modularity of metabolic networks[J].IET Systems Biology,2007,1(5):280-285.
[10]Goto S,Nishioka T,Kanehisa M.LIGAND:chemical database for enzyme reactions[J].Bioinformatics,1998,14(7):591-599.
[11]Maltsev N,Glass E,Sulakhe D,et al.PUMA2——grid-based high-throughput analysis of genomes and metabolic pathways[J].Nucl Acids Res,2006,34(suppl_1):369-372.
[12]Karp P D,Ouzounis C A,Moore-Kochlacs C,et al.Expansion of he BioCyc collection of pathway/genome databases to 160genomes[J].Nucl Acids Res,2005,33(19):6083-6089.
[13]Bondy J A,Murty U S R.Graph theory with applications[M].London:Macmillan,1976.
[14]Guimera R,Amaral L A N.Functional cartography of complex metabolic networks[J].Nature,2005,433(7028):895-900.
[15]Zhao J,Tao L,Yu H,et al.Bow-tie topological features of metabolic networks and the functional significance[J].Chinese Science Bulletin,2007,52(8):1036-1045.