基于調(diào)控-代謝網(wǎng)絡(luò)整合模型預(yù)測(cè)微生物生長(zhǎng)表型的方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于微生物技術(shù)領(lǐng)域，具體的說，涉及一種基因調(diào)控網(wǎng)絡(luò)與代謝網(wǎng)絡(luò)聯(lián)合 建模的方法，其可用于預(yù)測(cè)微生物的生長(zhǎng)表型。
【背景技術(shù)】
[0002] 目前，代謝網(wǎng)絡(luò)建模是預(yù)測(cè)基因敲除對(duì)代謝表型影響的主要方法。代謝網(wǎng)絡(luò)包含 所有酶、代謝物和生化反應(yīng)之間的相互作用，全基因組范圍代謝網(wǎng)絡(luò)重構(gòu)與分析能夠發(fā)現(xiàn) 基因敲除、插入、異常表達(dá)、環(huán)境變化等對(duì)生物系統(tǒng)表型的影響。
[0003] 動(dòng)力學(xué)模擬是分析代謝網(wǎng)絡(luò)最直接有效的手段，但由于很多動(dòng)力學(xué)參數(shù)未知使 得大規(guī)模網(wǎng)絡(luò)的動(dòng)力學(xué)分析受到限制。而基于約束的代謝流分析、基元模式和極端途經(jīng) 分析、最優(yōu)化方法等得到更加廣泛的應(yīng)用，其中流平衡分析最為常用。流平衡分析（Flux Balance Analysis)是一種基于約束的代謝網(wǎng)絡(luò)分析方法，通過加入化學(xué)計(jì)量信息、熱力學(xué) 信息、代謝能力約束等，將系統(tǒng)行為約束在一個(gè)封閉的解空間內(nèi)，然后設(shè)置目標(biāo)函數(shù)通過 線性規(guī)劃得到最優(yōu)解。通常目標(biāo)函數(shù)設(shè)為細(xì)胞生物質(zhì)合成最大、生長(zhǎng)速率最快或某個(gè)目標(biāo) 代謝物產(chǎn)量最高。
[0004] 以大腸桿菌和酵母為代表的微生物代謝模型研究中，流平衡分析能夠比較準(zhǔn)確地 模擬基因敲除對(duì)細(xì)胞生長(zhǎng)的影響。流平衡分析可以不依賴于反應(yīng)動(dòng)力學(xué)參數(shù)等信息，快速 識(shí)別對(duì)生長(zhǎng)或目標(biāo)代謝物合成具有重要影響的必需基因，如果某個(gè)基因敲除造成生長(zhǎng)速率 為零，則屬于致死基因，不能改動(dòng)；如果某個(gè)基因敲除或過表達(dá)能夠提高目標(biāo)代謝物產(chǎn)量， 則成為候選的改造基因。通過這種計(jì)算模擬的手段能夠?qū)ふ铱尚械母脑焱緩?，?duì)于代謝工 程高產(chǎn)菌株設(shè)計(jì)具有重要的指導(dǎo)意義。
[0005] 如果能夠?qū)⒒蛘{(diào)控網(wǎng)絡(luò)整合進(jìn)來，采用合理方式設(shè)定約束，將有助于提高代謝 流的預(yù)測(cè)，發(fā)現(xiàn)生物在不同環(huán)境下的響應(yīng)變化，使得現(xiàn)有的代謝模型能更加精確地模擬基 因敲除對(duì)生長(zhǎng)表型的影響?；蛘{(diào)控網(wǎng)絡(luò)描述調(diào)控因子與目標(biāo)基因之間的相互作用關(guān)系， 近年來有很多研究分別對(duì)兩種網(wǎng)絡(luò)進(jìn)行構(gòu)建和分析，但是將二者整合起來揭示轉(zhuǎn)錄調(diào)控對(duì) 代謝表型影響的研究較少。這是由于調(diào)控網(wǎng)絡(luò)是基于統(tǒng)計(jì)推斷的，而代謝網(wǎng)絡(luò)是基于生化 知識(shí)的，如何實(shí)現(xiàn)基因組范圍的調(diào)控網(wǎng)絡(luò)與代謝網(wǎng)絡(luò)整合是極富挑戰(zhàn)的問題。

【發(fā)明內(nèi)容】

[0006] 針對(duì)現(xiàn)有技術(shù)的不足，本發(fā)明的目的在于提供一種基于調(diào)控-代謝網(wǎng)絡(luò)預(yù)測(cè)微生 物在不同環(huán)境下生長(zhǎng)表型的方法。該方法簡(jiǎn)單、高效，預(yù)測(cè)結(jié)果準(zhǔn)確。
[0007] 本發(fā)明提供一種基于調(diào)控-代謝網(wǎng)絡(luò)整合模型預(yù)測(cè)微生物生長(zhǎng)表型的方法，具體 步驟如下：
[0008] (1)首先構(gòu)建基因調(diào)控網(wǎng)絡(luò)，基于大量基因表達(dá)譜數(shù)據(jù)，用線性回歸推斷每個(gè) 目標(biāo)基因隨轉(zhuǎn)錄因子表達(dá)變化的線性方程，然后隨機(jī)抽取部分表達(dá)譜數(shù)據(jù)進(jìn)行若干次 bootstrap線性回歸，根據(jù)轉(zhuǎn)錄因子和目標(biāo)基因的調(diào)控關(guān)系在這若干個(gè)線性方程中出現(xiàn)與 否的概率，計(jì)算錯(cuò)誤發(fā)現(xiàn)率FDR (False Discovery Rate)，F(xiàn)DR表明該轉(zhuǎn)錄因子沒有調(diào)控該 基因的比例。
[0009] (2)然后取FDR〈 = 0· 05的調(diào)控關(guān)系作為全局的調(diào)控網(wǎng)絡(luò)；
[0010] (3)接著找出調(diào)控網(wǎng)絡(luò)中調(diào)控代謝基因的轉(zhuǎn)錄因子，根據(jù)轉(zhuǎn)錄因子的類型，設(shè)定敲 除該轉(zhuǎn)錄因子時(shí)目標(biāo)代謝基因表達(dá)的概率P (Gene = ON I Factor = OFF)，如果是激活因子， 則P = FDR，如果是抑制因子，則P = I-FDR ;
[0011] (4)由代謝網(wǎng)絡(luò)的流變分析FVA(Flux Variability Analysis)可得每個(gè)反應(yīng)的最 大流值Vmax，對(duì)每個(gè)轉(zhuǎn)錄因子進(jìn)行模擬敲除時(shí)，對(duì)每個(gè)反應(yīng)施加相應(yīng)的約束為Vmax*P，接 著進(jìn)行目標(biāo)函數(shù)為生長(zhǎng)速率最大的流平衡分析，得出細(xì)胞生長(zhǎng)速率對(duì)應(yīng)的反應(yīng)流值F。
[0012] (5)對(duì)于沒有任何基因敲除的野生型代謝網(wǎng)絡(luò)進(jìn)行同樣的流平衡分析得出最大細(xì) 胞生長(zhǎng)速率Fmax，然后計(jì)算突變體相對(duì)于野生型的生長(zhǎng)速率之比F/Fmax，即代表生長(zhǎng)表型 的變化。
[0013] 本發(fā)明中，隨機(jī)抽取部分表達(dá)譜數(shù)據(jù)進(jìn)行150-300次bootstrap線性回歸。
[0014] 本發(fā)明的有益效果在于，對(duì)比不施加其他約束的流平衡分析，結(jié)合了基因調(diào)控網(wǎng) 絡(luò)的信息，使得分析的精確度提高，得以更好地預(yù)測(cè)微生物的生長(zhǎng)表型。
【附圖說明】
[0015] 圖1是本發(fā)明方法預(yù)測(cè)生長(zhǎng)表型的總體流程示意圖。
[0016] 圖2是以葡萄糖一銨為細(xì)胞生長(zhǎng)底物時(shí)的ROC曲線的結(jié)果圖示。
[0017] 圖3是以半乳糖一銨為細(xì)胞生長(zhǎng)底物時(shí)的ROC曲線的結(jié)果圖示。
[0018] 圖4是以葡萄糖一尿素為細(xì)胞生長(zhǎng)底物時(shí)的ROC曲線的結(jié)果圖示。
【具體實(shí)施方式】
[0019] 以下結(jié)合附圖和實(shí)施例對(duì)本發(fā)明的技術(shù)方案作進(jìn)一步描述。
[0020] 本發(fā)明的實(shí)施例中，采用的算法流程圖如圖1所示。本實(shí)施例以酵母為例，首先基 于已搜集的2929組酵母的基因表達(dá)譜數(shù)據(jù)，用線性回歸推斷每個(gè)目標(biāo)基因隨轉(zhuǎn)錄因子表 達(dá)變化的線性方程，若某轉(zhuǎn)錄因子在方程中的系數(shù)為正，則代表有激活作用，若系數(shù)為負(fù)， 則代表有抑制作用，若系數(shù)為零，則代表沒有調(diào)控作用。然后隨機(jī)抽取2929組表達(dá)譜數(shù)據(jù) 中的子集進(jìn)行200次bootstrap線性回歸，根據(jù)轉(zhuǎn)錄因子和目標(biāo)基因的調(diào)控關(guān)系在這200 個(gè)線性方程中出現(xiàn)與否的概率，計(jì)算錯(cuò)誤發(fā)現(xiàn)率FDR ;然后取FDR〈 = 0. 05的可靠調(diào)控關(guān) 系，形成包含2588對(duì)轉(zhuǎn)錄因子與目標(biāo)基因的相互作用的酵母全局調(diào)控網(wǎng)絡(luò)。
[0021] 接下來輸入酵母的代謝網(wǎng)絡(luò)，以Yeast6模型為例，該模型包含1888個(gè)反應(yīng)，1458 個(gè)代謝物，900個(gè)編碼酶的基因，用流平衡分析得出最大細(xì)胞生長(zhǎng)速率Fmax。找出其中有調(diào) 控關(guān)系的代謝基因，根據(jù)調(diào)控其的轉(zhuǎn)錄因子類型，設(shè)定轉(zhuǎn)錄因子敲除時(shí)，該代謝基因表達(dá)的 概率P (Gene = ON I Factor = OFF)，如果是激活因子，則P = FDR，如果是抑制因子，則P =I-FDR ;再對(duì)完整的Yeast6代謝網(wǎng)絡(luò)進(jìn)行流變分析，可得每個(gè)反應(yīng)的最大流值Vmax，則對(duì) 每個(gè)轉(zhuǎn)錄因子進(jìn)行模擬敲除時(shí)，將對(duì)每個(gè)反應(yīng)施加新的約束為Vmax*P ;接著進(jìn)行目標(biāo)函數(shù) 為生長(zhǎng)速率最大的流平衡分析，得出細(xì)胞生長(zhǎng)速率的反應(yīng)流值F ;最后計(jì)算突變體相對(duì)于 野生型的生長(zhǎng)速率之比F/Fmax。
[0022] 下面驗(yàn)證預(yù)測(cè)結(jié)果的有效性，實(shí)驗(yàn)數(shù)據(jù)來源于文獻(xiàn)i，共有三種環(huán)境下的實(shí)驗(yàn)數(shù) 據(jù)，分別以葡萄糖一銨，半乳糖一銨，葡萄糖一尿素為細(xì)胞生長(zhǎng)的底物，對(duì)此我們也預(yù)測(cè)了 這三種情況下敲除轉(zhuǎn)錄因子的生長(zhǎng)表型。將我們預(yù)測(cè)的突變體相對(duì)于野生型的生長(zhǎng)速率之 比與實(shí)驗(yàn)中對(duì)應(yīng)的比值進(jìn)行比較，與實(shí)驗(yàn)對(duì)應(yīng)上的有51個(gè)轉(zhuǎn)錄因子敲除的突變體，對(duì)每個(gè) 環(huán)境計(jì)算預(yù)測(cè)值與實(shí)驗(yàn)值這兩組數(shù)據(jù)的Pearson相關(guān)性，結(jié)果如表1所示。
[0023] 表1預(yù)測(cè)值與實(shí)驗(yàn)值這兩組數(shù)據(jù)的pearson相關(guān)性結(jié)果
[0025] 結(jié)果顯示p值均小于0.05,可以認(rèn)為我們的預(yù)測(cè)值與實(shí)驗(yàn)值顯著線性相關(guān)。另外， 對(duì)實(shí)驗(yàn)值設(shè)定閾值，將數(shù)據(jù)二元化，表示是否影響生長(zhǎng)；以此為標(biāo)準(zhǔn)，通過ROC曲線判斷預(yù) 測(cè)結(jié)果是否能正確反映轉(zhuǎn)錄因子敲除對(duì)突變體生長(zhǎng)的影響。
[0026] 圖2、圖3、圖4分別顯示了以葡萄糖一銨，半乳糖一銨，葡萄糖一尿素為細(xì)胞生長(zhǎng) 底物的三種環(huán)境下ROC曲線的結(jié)果，其中對(duì)實(shí)驗(yàn)值設(shè)定了三個(gè)閾值分別為0. 2、0. 5和0. 8 ; 總體來說，所有的ROC曲線都在對(duì)角線之上，曲線下面積值如各圖所示，因此說明預(yù)測(cè)的結(jié) 果可以有效判斷轉(zhuǎn)錄因子敲除對(duì)生長(zhǎng)的影響。綜上，我們通過DIRAM方法可以有效地預(yù)測(cè) 酵母在三種不同環(huán)境下的生長(zhǎng)表型。
[0027] 參考文獻(xiàn)：
[0028] l、Fendt S Mj Oliveira A P，Christen S，et al. Unraveling condition - dependent networks of transcription factors that control metabolic pathway activity in yeast[J]. Molecular systems biology,2010, 6(I):432.
【主權(quán)項(xiàng)】
1. 一種基于調(diào)控-代謝網(wǎng)絡(luò)整合模型預(yù)測(cè)微生物生長(zhǎng)表型的方法，其特征在于將調(diào)控 與代謝兩個(gè)層次有機(jī)結(jié)合，預(yù)測(cè)轉(zhuǎn)錄調(diào)控對(duì)生長(zhǎng)表型的影響，具體步驟如下： (1) 首先構(gòu)建基因調(diào)控網(wǎng)絡(luò)，基于大量基因表達(dá)譜數(shù)據(jù)，用線性回歸推斷每個(gè)目標(biāo)基因 隨轉(zhuǎn)錄因子表達(dá)變化的線性方程，然后隨機(jī)抽取部分表達(dá)譜數(shù)據(jù)進(jìn)行若干次bootstrap線 性回歸，根據(jù)轉(zhuǎn)錄因子和目標(biāo)基因的調(diào)控關(guān)系在這若干個(gè)線性方程中出現(xiàn)與否的概率，計(jì) 算錯(cuò)誤發(fā)現(xiàn)率FDR，F(xiàn)DR表明該轉(zhuǎn)錄因子沒有調(diào)控該基因的比例； (2) 然后取FDR〈 = 0. 05的調(diào)控關(guān)系作為全局的調(diào)控網(wǎng)絡(luò)； (3) 接著找出調(diào)控網(wǎng)絡(luò)中調(diào)控代謝基因的轉(zhuǎn)錄因子，根據(jù)轉(zhuǎn)錄因子的類型，設(shè)定敲除該 轉(zhuǎn)錄因子時(shí)目標(biāo)代謝基因表達(dá)的概率P(Gene= 0N|Factor=OFF)，如果是激活因子，則P =FDR，如果是抑制因子，則P=I-FDR; (4) 由代謝網(wǎng)絡(luò)的流變分析FVA可得每個(gè)反應(yīng)的最大流值Vmax，對(duì)每個(gè)轉(zhuǎn)錄因子進(jìn)行 模擬敲除時(shí)，對(duì)每個(gè)反應(yīng)施加相應(yīng)的約束為Vmax*P，接著進(jìn)行目標(biāo)函數(shù)為生長(zhǎng)速率最大的 流平衡分析，得出細(xì)胞生長(zhǎng)速率對(duì)應(yīng)的反應(yīng)流值F; (5) 對(duì)于沒有任何基因敲除的野生型代謝網(wǎng)絡(luò)進(jìn)行同樣的流平衡分析得出最大細(xì)胞生 長(zhǎng)速率Fmax，然后計(jì)算突變體相對(duì)于野生型的生長(zhǎng)速率之比F/Fmax，即代表生長(zhǎng)表型的變 化。2. 根據(jù)權(quán)利要求1所述的方法，其特征在于，步驟（1)中，隨機(jī)抽取部分表達(dá)譜數(shù)據(jù)進(jìn) 行150-300次bootstrap線性回歸。
【專利摘要】本發(fā)明公開了一種基于調(diào)控-代謝網(wǎng)絡(luò)整合模型預(yù)測(cè)微生物生長(zhǎng)表型的方法。該方法首先構(gòu)建基因調(diào)控網(wǎng)絡(luò)，基于大量基因表達(dá)譜數(shù)據(jù)，運(yùn)行多次線性回歸推斷每個(gè)目標(biāo)基因隨轉(zhuǎn)錄因子表達(dá)變化的線性方程，計(jì)算錯(cuò)誤發(fā)現(xiàn)率FDR；然后取FDR<＝0.05的調(diào)控關(guān)系作為全局的調(diào)控網(wǎng)絡(luò)，找出調(diào)控網(wǎng)絡(luò)中調(diào)控代謝基因的轉(zhuǎn)錄因子，根據(jù)轉(zhuǎn)錄因子的類型，計(jì)算轉(zhuǎn)錄因子敲除時(shí)生長(zhǎng)速率對(duì)應(yīng)的反應(yīng)流值F，并由原始代謝網(wǎng)絡(luò)進(jìn)行同樣的流平衡分析得到最大細(xì)胞生長(zhǎng)速率Fmax，通過計(jì)算F/Fmax預(yù)測(cè)微生物在轉(zhuǎn)錄因子敲除時(shí)的生長(zhǎng)表型的變化。本發(fā)明方法使得分析的精確度提高，得以更好地預(yù)測(cè)微生物的生長(zhǎng)表型。
【IPC分類】G06F19/00
【公開號(hào)】CN105184049
【申請(qǐng)?zhí)枴緾N201510484159
【發(fā)明人】王卓, 沈方舟
【申請(qǐng)人】上海交通大學(xué)
【公開日】2015年12月23日
【申請(qǐng)日】2015年8月10日