本發(fā)明屬于植物分子生物技術(shù)和基因工程領(lǐng)域，具體涉及一種基于基因表達(dá)量與性狀動態(tài)相關(guān)性預(yù)測玉米未知基因功能的方法。

背景技術(shù)：

玉米是世界三大主要農(nóng)作物之一，20世紀(jì)90年代以來，世界玉米總產(chǎn)量首次超越水稻和小麥，成為第一位的糧食作物。玉米籽粒中累積了大量的儲存物質(zhì)包括淀粉、油份和蛋白。隨著人們生活水平的提高和膳食結(jié)構(gòu)的變化，以及淀粉和油脂加工業(yè)的發(fā)展，玉米品種由產(chǎn)量型逐漸向質(zhì)量型轉(zhuǎn)變，玉米品質(zhì)及其專用性變得越來越重要。

復(fù)雜的性狀受到多個基因位點的調(diào)控，而基因之間相互作用又形成復(fù)雜的基因調(diào)控網(wǎng)絡(luò)，來控制細(xì)胞中各種生物學(xué)反應(yīng)的進(jìn)行。高通量測序技術(shù)的發(fā)展使我們獲得了大規(guī)模海量的組學(xué)數(shù)據(jù)，如基因型數(shù)據(jù)、基因表達(dá)量數(shù)據(jù)、蛋白相互作用數(shù)據(jù)等。研究表明，功能相近基因，其表達(dá)模式相關(guān)。因此構(gòu)建共表達(dá)網(wǎng)絡(luò)為預(yù)測未知基因功能提供了思路。然而，在構(gòu)建共表達(dá)網(wǎng)絡(luò)過程中，我們發(fā)現(xiàn)很多功能相近基因，其表達(dá)模式并不相關(guān)。因此，利用共表達(dá)分析預(yù)測未知基因功能具有一定的局限性。研究表明單個基因/蛋白本身對復(fù)雜數(shù)量性狀的影響有限，往往需要通過高階的細(xì)胞組織形式來行使功能，很多功能相關(guān)基因的表達(dá)量并不相關(guān)，挖掘控制表型性狀的遺傳位點，而這些遺傳位點之間相對獨立，它們之間的調(diào)控關(guān)系未知，并且傳統(tǒng)分析方法需要多年多點的表型鑒定，費時費力。

技術(shù)實現(xiàn)要素：

針對現(xiàn)有技術(shù)存在的問題，本發(fā)明提供了一種基于基因表達(dá)量與性狀動態(tài)相關(guān)性預(yù)測玉米未知基因功能的方法，該方法通過動態(tài)關(guān)聯(lián)分析，鑒定玉米基因組基因共表達(dá)模式的動態(tài)關(guān)聯(lián)，根據(jù)la結(jié)果，預(yù)測未知基因功能。

本發(fā)明是通過以下技術(shù)方案來實現(xiàn)的：

本發(fā)明提供了一種基于基因表達(dá)量與性狀動態(tài)相關(guān)性預(yù)測玉米未知基因功能的方法，包括以下步驟：

（1）收集玉米自交系授粉后15天的籽粒轉(zhuǎn)錄本測序獲得基因表達(dá)量數(shù)據(jù)；

（2）動態(tài)關(guān)聯(lián)分析la模型的建立；

（3）la顯著性評估；

（4）挖掘玉米全基因組基因共表達(dá)模式的動態(tài)關(guān)聯(lián)；

（5）對顯著la結(jié)果的基因進(jìn)行功能注釋，預(yù)測未知基因的功能。

進(jìn)一步的，所述玉米自交系分成了2組：熱帶和亞熱帶、溫帶，小組內(nèi)采用完全隨機(jī)區(qū)組法，設(shè)2個重復(fù)，每個自交系每個重復(fù)播種1行，所有材料均進(jìn)行自交，收獲未成熟的授粉后15天的籽粒，每個自交系的兩個重復(fù)各取3-4穗，每穗取1-2粒籽粒，混合提取籽?？俽na，隨機(jī)選擇自交系數(shù)量的樣品用于rna-seq。

上述rna-seq包括以下步驟：首先，用ploy（t）寡聚核苷酸從總rna中抽取全部帶ploy（a）尾的rna，主要為mrna，然后將截獲的mrna隨機(jī)打斷成片段，用六堿基隨機(jī)引物合成cdna第一鏈，并加入逆轉(zhuǎn)錄酶合成cdna第二鏈，經(jīng)過試劑盒純化并對cdna片段進(jìn)行末端修飾，連接測序接頭，再經(jīng)瓊脂糖凝膠電泳回收目的大小片段，進(jìn)行pcr擴(kuò)增，用illuminagaⅱ基因分析系統(tǒng)進(jìn)行序列測定及分析，獲得基因的表達(dá)量數(shù)據(jù)。

進(jìn)一步的，所述動態(tài)關(guān)聯(lián)分析la模型具體采用以下方法建立：la的數(shù)學(xué)定義如下：

la(x,y|z)=eg'(z)公式1

其中，所述x、y和z為玉米籽粒中基因表達(dá)量數(shù)據(jù)；假設(shè)x，y，z是均值為0，方差為1的連續(xù)隨機(jī)變量，則x，y的相關(guān)性表示為e(xy)；當(dāng)z=z時，g(z)=e(xy|z=z)，g(z)檢測的是當(dāng)z=z時，xy基因?qū)Φ墓脖磉_(dá)模式。g(z)的導(dǎo)數(shù)表示為g'(z)，該值可用于共表達(dá)模式變化的期望測定，

當(dāng)z符合標(biāo)準(zhǔn)正態(tài)分布時，la值可簡單的表示為la(x,y|z)=e(xyz)。

x，y，z代表具有正態(tài)分布表達(dá)譜的三個基因，則la(x,y|z)表示為：e(xyz)=(x1y1z1+x2y2z2+...+xmymzm)/m公式2

la用來反應(yīng)基因?qū)脖磉_(dá)模式的動態(tài)變化，即當(dāng)z基因表達(dá)量較高時，xy基因?qū)Φ谋磉_(dá)量呈正相關(guān)（co-regulated），e(xy|z=1)為正數(shù)；當(dāng)z基因表達(dá)量較低時，xy基因?qū)Φ谋磉_(dá)量呈負(fù)相關(guān)（contra-regulated，），e(xy|z=0)為負(fù)數(shù)，因此基因?qū)Φ谋磉_(dá)調(diào)控模式由正相關(guān)（co-regulated）轉(zhuǎn)變?yōu)樨?fù)相關(guān)（contra-regulated），la值記為正；相反，基因?qū)Φ谋磉_(dá)調(diào)控模式由負(fù)相關(guān)（contra-regulated）轉(zhuǎn)變?yōu)檎嚓P(guān)（co-regulated），la值記為負(fù)。

本發(fā)明所建立的動態(tài)關(guān)聯(lián)分析模型的評估步驟如下：混合所有基因的表達(dá)量值；在每次模擬中，用放回隨機(jī)抽樣法隨機(jī)抽取一對基因（x，y）的表達(dá)量值,z基因取全基因組所有基因，計算xy基因?qū)υ谌蚪M的la值，可分別得到la的正極大值和負(fù)極小值；重復(fù)模擬一百萬次，分別得到la的正值參考分布和負(fù)值參考分布，用la正負(fù)參考分布的99%分位數(shù)作為la正負(fù)顯著性閾值。

進(jìn)一步的，所述全基因組動態(tài)關(guān)聯(lián)分析的結(jié)果按照la值的大小進(jìn)行過濾，對顯著la的基因進(jìn)行功能注釋，預(yù)測未知基因功能。

研究表明功能相近基因表達(dá)模式不相關(guān)的原因主要包含以下兩個假設(shè)，一是這些功能相近基因的表達(dá)調(diào)控不在mrna水平上，二是功能相近基因的表達(dá)模式只在特定的細(xì)胞環(huán)境中才相關(guān)，即共表達(dá)模式的動態(tài)關(guān)聯(lián)，動態(tài)關(guān)聯(lián)分析（liquidassociation,la）為驗證第二種假設(shè)提供了理論支持。本發(fā)明基于功能相近基因，表達(dá)模式相關(guān)的科學(xué)假設(shè)，采用la方法鑒定玉米全基因組基因共表達(dá)模式的動態(tài)關(guān)聯(lián)，根據(jù)顯著la結(jié)果中基因的功能注釋，預(yù)測未知基因功能，并根據(jù)未知基因在擬南芥中的同源基因功能，驗證la預(yù)測結(jié)果，思路創(chuàng)新，在植物學(xué)領(lǐng)域該項研究尚無報道。

本發(fā)明的有益效果為：

（1）本發(fā)明以玉米籽粒中基因?qū)脖磉_(dá)模式動態(tài)關(guān)聯(lián)這一現(xiàn)象為突破口，預(yù)測未知基因功能。相比較于傳統(tǒng)的共表達(dá)網(wǎng)絡(luò)構(gòu)建，動態(tài)關(guān)聯(lián)分析可以快速找到調(diào)控共表達(dá)模式的調(diào)控基因；

（2）本發(fā)明通過對顯著la結(jié)果的基因進(jìn)行功能注釋，推測未知基因功能，并通過同源基因的功能驗證預(yù)測結(jié)果，是預(yù)測未知基因功能的有效方法。

附圖說明

圖1為隨機(jī)模擬生成la值評估la分析的顯著性。

具體實施方式

下面結(jié)合附圖和具體實施例對本發(fā)明作進(jìn)一步說明，下述說明僅是實例性的，不限定本發(fā)明的保護(hù)范圍。

實施例1

一種本發(fā)明所述基于動態(tài)關(guān)聯(lián)分析預(yù)測玉米未知基因功能的方法，主要包括三步，基因表達(dá)量數(shù)據(jù)的收集、la顯著性評估和全基因組la分析。

（1）基因表達(dá)量數(shù)據(jù)的收集：

368份自交系（本發(fā)明所使用的玉米品種可為任意品種，包括中國農(nóng)業(yè)大學(xué)宋同明教授培育的35份高油玉米自交系（yang等，2010b））于2010年在湖北荊州種植，根據(jù)系譜信息分成了2組（熱帶和亞熱帶、溫帶），小組內(nèi)采用完全隨機(jī)區(qū)組法，設(shè)2個重復(fù)，每個自交系每個重復(fù)播種1行。所有材料均進(jìn)行自交，收獲未成熟的授粉后15天（15dap）的籽粒，每個自交系的兩個重復(fù)各取3-4穗，每穗取1-2粒籽粒，混合提取籽?？俽na，隨機(jī)選擇368個樣品用于rna-seq。樣品的rna-seq工作是由深圳華大基因研究院（beijinggenomicsinstitute，bgi）完成，測序方法簡要描述如下：首先，用ploy（t）寡聚核苷酸從總rna中抽取全部帶ploy（a）尾的rna，主要為mrna，然后將截獲的mrna隨機(jī)打斷成片段，用六堿基隨機(jī)引物（randomhexamers）合成cdna第一鏈，并加入逆轉(zhuǎn)錄酶等合成cdna第二鏈，經(jīng)過試劑盒（ampurexpbeads）純化并對cdna片段進(jìn)行末端修飾，連接測序接頭，再經(jīng)瓊脂糖凝膠電泳回收目的大小片段，進(jìn)行pcr擴(kuò)增，從而完成整個文庫構(gòu)建工作，構(gòu)建好的文庫用illuminagaⅱ基因分析系統(tǒng)進(jìn)行序列測定及分析。轉(zhuǎn)錄本測序獲得的368個玉米自交系中28769個基因的表達(dá)量數(shù)據(jù)，對基因表達(dá)量數(shù)據(jù)集進(jìn)行的缺失值預(yù)處理如下：基因表達(dá)數(shù)據(jù)因為實驗中的噪聲、檢測技術(shù)等原因而存在缺失。對于數(shù)據(jù)集中的每個基因，如果其表達(dá)值在高于30%的樣本中缺失，則在后續(xù)的分析中舍棄該基因。共獲得24,907個基因表達(dá)量數(shù)據(jù)（該部分?jǐn)?shù)據(jù)根據(jù)需要可以直接從數(shù)據(jù)庫下載）用于后續(xù)的全基因組la分析；

（2）動態(tài)關(guān)聯(lián)分析la模型的建立：

所述動態(tài)關(guān)聯(lián)分析la模型具體采用以下方法建立：la的數(shù)學(xué)定義如下：

la(x,y|z)=eg'(z)公式1

其中，所述x、y和z為玉米籽粒中基因表達(dá)量數(shù)據(jù)；假設(shè)x，y，z是均值為0，方差為1的連續(xù)隨機(jī)變量，則x，y的相關(guān)性表示為e(xy)；當(dāng)z=z時，g(z)=e(xy|z=z)，g(z)檢測的是當(dāng)z=z時，xy基因?qū)Φ墓脖磉_(dá)模式。g(z)的導(dǎo)數(shù)表示為g'(z)，該值可用于共表達(dá)模式變化的期望測定，

當(dāng)z符合標(biāo)準(zhǔn)正態(tài)分布時，la值可簡單的表示為la(x,y|z)=e(xyz)。

x，y，z代表具有正態(tài)分布表達(dá)譜的三個基因，則la(x,y|z)表示為：e(xyz)=(x1y1z1+x2y2z2+...+xmymzm)/m公式2

la用來反應(yīng)基因?qū)脖磉_(dá)模式的動態(tài)變化，即當(dāng)z基因表達(dá)量較高時，xy基因?qū)Φ谋磉_(dá)量呈正相關(guān)（co-regulated），e(xy|z=1)為正數(shù)；當(dāng)z基因表達(dá)量較低時，xy基因?qū)Φ谋磉_(dá)量呈負(fù)相關(guān)（contra-regulated，），e(xy|z=0)為負(fù)數(shù)，因此基因?qū)Φ谋磉_(dá)調(diào)控模式由正相關(guān)（co-regulated）轉(zhuǎn)變?yōu)樨?fù)相關(guān)（contra-regulated），la值記為正；相反，基因?qū)Φ谋磉_(dá)調(diào)控模式由負(fù)相關(guān)（contra-regulated）轉(zhuǎn)變?yōu)檎嚓P(guān)（co-regulated），la值記為負(fù)。

（3）la顯著性評估

混合所有基因的表達(dá)量值；在每次模擬中，用放回隨機(jī)抽樣法隨機(jī)抽取一對基因（x，y）的表達(dá)量值,z基因取全基因組所有基因，計算xy基因?qū)υ谌蚪M的la值，可分別得到la的正極大值和負(fù)極小值；重復(fù)模擬一百萬次，分別得到la的正值參考分布和負(fù)值參考分布，用la正負(fù)參考分布的99%分位數(shù)作為la正負(fù)顯著性閾值，具體結(jié)果見圖1。

（4）全基因組la分析

以x=全基因組基因，y=全基因組基因，z=全基因組基因進(jìn)行l(wèi)a分析，重點關(guān)注la絕對值最大的前50個共表達(dá)基因?qū)Γ╨iquidassociationpari,lap）列表。對x、y和z進(jìn)行功能注釋，表1為grmzm5g858880調(diào)控的參與蛋白翻譯過程的基因列表?；?i>grmzm5g858880調(diào)控多對共表達(dá)基因?qū)Γ╨iquidassociationpair,lap），玉米基因組數(shù)據(jù)庫（maizegenomedatabase）對該基因的功能注釋為“編碼包含ww功能域的蛋白”。在grmzm5g858880調(diào)控的lap列表中，發(fā)現(xiàn)部分基因參與蛋白質(zhì)翻譯過程，包括核糖體蛋白合成、蛋白翻譯起始以及蛋白磷酸化，并出現(xiàn)多次，grmzm2g092663(編碼核糖體s5蛋白家族，4次)，grmzm2g099352（編碼核糖體s3蛋白家族），grmzm2g168149（編碼核糖體s5蛋白家族），grmzm2g129015（編碼核糖體s26e蛋白家族，2次），grmzm2g164352（編碼蛋白磷酸酶2a亞基a2，4次），grmzm2g122135（編碼蛋白磷酸酶2a亞基a2，2次），grmzm2g064133（編碼真核生物翻譯起始因子3g1），因此推測調(diào)控基因grmzm5g858880也參與蛋白翻譯過程。研究表明，grmzm5g858880在擬南芥中的同源基因（at3g13225）通過核糖體減速和降低再起始效率來調(diào)控蛋白翻譯過程（tran等，bmcgenomics，2008）。

表1grmzm5g858880調(diào)控的參與蛋白翻譯過程的基因列表

以上的結(jié)果證明了本發(fā)明的有效性，通過全基因組基因?qū)脖磉_(dá)模式的動態(tài)關(guān)聯(lián)分析，并結(jié)合功能注釋，預(yù)測未知基因功能，為玉米功能基因組學(xué)研究提供了新的思路和方法。