本發(fā)明涉及分子生物學(xué)技術(shù)領(lǐng)域，具體涉及一種植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋的方法。

背景技術(shù)：

高通量rna測(cè)序技術(shù)是現(xiàn)代基因組學(xué)研究最重要的實(shí)驗(yàn)技術(shù)，在整個(gè)生物學(xué)領(lǐng)域具有廣泛的應(yīng)用。隨著大量測(cè)序數(shù)據(jù)的建立，基因組學(xué)研究得到快速發(fā)展。在利用高通量rna測(cè)序技術(shù)運(yùn)用于差異表達(dá)基因即編碼rna的研究基礎(chǔ)上，越來(lái)越多的研究關(guān)注到基因組的非編碼rna上。這些非編碼rna包括mirnas、sirnas、lncrnas、circularrna等，其中mirnas由于其具有獨(dú)特序列模塊以及與靶基因的作用方式已經(jīng)被廣泛地證明在植物生長(zhǎng)發(fā)育和響應(yīng)逆境脅迫的轉(zhuǎn)錄調(diào)控方面具有重要的作用。而lncrnas是一類長(zhǎng)度超過(guò)200nt,不具有蛋白編碼能力的非編碼rna，可以通過(guò)多種方式影響靶基因的轉(zhuǎn)錄調(diào)控。因此，高通量的開展lncrnas的功能注釋工作對(duì)于今后非編碼rna的功能研究具有重要意義。

目前，已有的研究主要通過(guò)計(jì)算非編碼rna與編碼rna共表達(dá)關(guān)系的方法將lncrnas與共表達(dá)的mrna分為一組，利用mrna的注釋信息來(lái)完成lncrnas的功能注釋。該方法一是需要大量的表達(dá)數(shù)據(jù)用于共表達(dá)分析，同時(shí)還將缺失脅迫響應(yīng)特異表達(dá)lncrnas的功能注釋信息。

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

有鑒于此，本發(fā)明的第一目的在于提供了一種植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法，包括以下步驟：

步驟s1，篩選植物響應(yīng)逆境脅迫的lncrnas及其靶基因；

步驟s2，再次篩選植物響應(yīng)逆境脅迫的lncrnas及其靶基因表達(dá)模式分析；

步驟s3，按照靶基因表達(dá)模式進(jìn)行植物響應(yīng)逆境脅迫的lncrnas分組；

步驟s4，響應(yīng)逆境脅迫的lncrnas特異序列模塊富集分析；

步驟s5，響應(yīng)逆境脅迫的lncrnas特異序列模塊功能注釋。

優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述步驟s1中所述逆境脅迫為低溫處理或高鹽處理；優(yōu)選地，所述低溫處理為4℃處理6小時(shí)；所述高鹽處理為150mmnacl處理6小時(shí)。

優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述步驟s1中l(wèi)ncrnas篩選標(biāo)準(zhǔn)為①長(zhǎng)度大于200nt；②最小讀長(zhǎng)覆蓋率為3；③開放閱讀框小于300nt；④cpcscore<0,cnciscore<0。

優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述步驟s2中所述再次篩選響應(yīng)逆境脅迫差異表達(dá)的lncrnas，篩選最小閾值為：差異倍數(shù)>2或<0.5，p-值<0.05，q-值<0.05。

優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述步驟s3中所述分組方法為根據(jù)響應(yīng)逆境脅迫的lncrnas表達(dá)模式或靶基因功能富集結(jié)果，將響應(yīng)逆境脅迫的lncrnas進(jìn)行分組。

優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述步驟s4中響應(yīng)逆境脅迫的lncrnas特異序列模塊富集分析方法為利用meme(http://meme-suite.org/tools/meme)對(duì)分組后的響應(yīng)逆境脅迫的lncrnas特異序列模塊進(jìn)行富集分析。

更優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述分析中篩選的參數(shù)設(shè)置為：①每個(gè)lncrnas上至少預(yù)測(cè)到一個(gè)序列模塊；②預(yù)測(cè)序列模塊數(shù)量3-5個(gè)；③模塊長(zhǎng)度為6bp-15bp；④分布模式為正反義兩條鏈；⑤不允許軟件對(duì)序列進(jìn)行重排。

優(yōu)選地，本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述步驟s5中所述響應(yīng)逆境脅迫的lncrnas特異序列模塊功能注釋為利用gomo(http://meme-suite.org//tools/gomo)，對(duì)富集的lncrnas特異序列模塊功能進(jìn)行功能預(yù)測(cè)。

更優(yōu)選地，期望值(e-value)表示由于隨機(jī)性造成獲得與數(shù)據(jù)庫(kù)比對(duì)結(jié)果的可能次數(shù)。期望值越小，發(fā)生這一事件的概率越低，比對(duì)結(jié)果越顯著。假設(shè)機(jī)率(p-value)代表給定原假設(shè)為真時(shí)樣本結(jié)果出現(xiàn)的概率。p值越小代表原假設(shè)不成立的概率越大。q值(q-value)代表經(jīng)過(guò)假陽(yáng)性率校正之后的p值，該數(shù)值越低則代表假陽(yáng)性率越小。特異性(specificity)代表該模塊的獨(dú)特性，數(shù)值越高則獨(dú)特性越強(qiáng)。本發(fā)明所述的植物響應(yīng)逆境脅迫的lncrnas序列模塊功能注釋方法中，所述注釋方法中篩選參數(shù)設(shè)置如下：①顯著性參數(shù)設(shè)置為p-value<0.01并且q-value<0.01；②計(jì)算次數(shù)5000次；③專一性(specificity)參數(shù)大于80％；④期望值e-value<0.001。序列模塊1(e-value＝2.7e-007；p-value＝2.6e-07；q-value＝1.1e-04；specificity＝83％)獲得goterm注釋轉(zhuǎn)錄因子活性(見圖2)。

因此，針對(duì)現(xiàn)有l(wèi)ncrnas功能預(yù)測(cè)方法依賴于需要大量的表達(dá)數(shù)據(jù)用于共表達(dá)分析，以及還沒有響應(yīng)逆境脅迫的特異表達(dá)lncrnas的功能注釋信息。本方案提供一種響應(yīng)逆境脅迫的植物lncrnas序列模塊功能注釋方法，結(jié)合生物信息學(xué)與差異表達(dá)分析對(duì)植物響應(yīng)逆境脅迫的lncrnas序列模塊進(jìn)行注釋，將為系統(tǒng)解析lncrns生物學(xué)功能以及構(gòu)建其轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)提供了技術(shù)支持。

即本發(fā)明提供一種植物lncrnas逆境脅迫響應(yīng)序列模塊功能注釋的方法，結(jié)合生物信息學(xué)與差異表達(dá)分析對(duì)植物響應(yīng)逆境脅迫的lncrnas序列模塊進(jìn)行注釋，不但極大地提高了實(shí)驗(yàn)的效率、精準(zhǔn)性以及靈活性，并顯著地降低了實(shí)驗(yàn)成本。

附圖說(shuō)明

圖1為響應(yīng)逆境脅迫的植物lncrnas序列模塊功能注釋的方法流程示意圖；

圖2毛白楊響應(yīng)低溫lncrnas序列模塊功能注釋結(jié)果圖；

圖3小葉楊響應(yīng)高鹽脅迫lncrnas序列模塊功能注釋結(jié)果圖。

具體實(shí)施方式

根據(jù)本發(fā)明一個(gè)典型的實(shí)施方式，待測(cè)樣本為毛白楊1年生植株。對(duì)其進(jìn)行低溫(4℃,6小時(shí))處理，立即收集葉片用于提取總rnas。利用ribo-zerorrna試劑盒對(duì)核糖體rna進(jìn)行去除。利用smart試劑盒進(jìn)行鏈特異性cdna文庫(kù)構(gòu)建。利用illuminahiseq^tm2500測(cè)序平臺(tái)完成cdna文庫(kù)測(cè)序，測(cè)序深度為10×。去除接頭以及冗余序列，通過(guò)cufflinks軟件拼接轉(zhuǎn)錄本，篩選長(zhǎng)度大于200nt、最小讀長(zhǎng)覆蓋率為5、cpcscore<0、cnciscore<0以及與對(duì)照差異表達(dá)倍數(shù)大于2(p<0.05)的lncrnas。預(yù)測(cè)差異表達(dá)的lncrnas順式及反式作用靶基因，并對(duì)靶基因的表達(dá)模式進(jìn)行解析，篩選差異表達(dá)的靶基因作為候選基因(最小閾值為：差異倍數(shù)>2或<0.5,p-值<0.05,q-值<0.05)。對(duì)于候選基因，利用ncbi核酸數(shù)據(jù)庫(kù)(https://blast.ncbi.nlm.nih.gov/blast.cgi)進(jìn)行功能注釋。利用agrigo(http://bioinfo.cau.edu.cn/agrigo/)對(duì)候選基因進(jìn)行g(shù)oterm富集分析。利用meme(http://meme-suite.org/tools/meme)對(duì)分組后的響應(yīng)逆境脅迫的lncrnas特異序列模塊進(jìn)行富集分析。利用gomo(http://meme-suite.org//tools/gomo)，對(duì)富集的lncrnas特異序列模塊功能進(jìn)行功能預(yù)測(cè)。

以下通過(guò)具體實(shí)施例進(jìn)一步對(duì)本發(fā)明的技術(shù)方案進(jìn)行說(shuō)明，應(yīng)理解以下僅為本發(fā)明的示例性說(shuō)明，并不用于限制本發(fā)明權(quán)利要求的保護(hù)范圍。

實(shí)施例1

對(duì)毛白楊1年生植株進(jìn)行低溫(4℃,6小時(shí))處理，提取其總rna用于對(duì)響應(yīng)低溫脅迫的lncrnas序列模塊進(jìn)行功能注釋。

s1，利用ribo-zerorrna試劑盒對(duì)核糖體rna進(jìn)行去除。利用smart試劑盒進(jìn)行鏈特異性cdna文庫(kù)構(gòu)建。利用illuminahiseq^tm2500測(cè)序平臺(tái)完成cdna文庫(kù)測(cè)序，測(cè)序深度為10×。去除接頭以及冗余序列，通過(guò)cufflinks軟件拼接轉(zhuǎn)錄本，篩選長(zhǎng)度大于200nt、最小讀長(zhǎng)覆蓋率為5、cpcscore<0、cnciscore<0lncrnas，共獲得4218個(gè)。篩選差異倍數(shù)小于0.35且大于0.2(p-值<0.05,q-值<0.05)的lncrnas，共17個(gè)(見表1)。

表1毛白楊響應(yīng)低溫逆境脅迫的lncrnas

在14個(gè)響應(yīng)高溫脅迫的lncrna上下游10kb范圍內(nèi)篩選順式作用靶基因，共獲得26個(gè)(見表2)。利用blast進(jìn)行序列互補(bǔ)計(jì)算，參數(shù)設(shè)置為e-value＝1e-10，identity＝90％以及利用rnaplex進(jìn)行熱力學(xué)上的互補(bǔ)計(jì)算，參數(shù)設(shè)置為e＝-70。共獲得反式作用靶基因25個(gè)(見表3)。根據(jù)差異表達(dá)靶基因篩選標(biāo)準(zhǔn)(p-值<0.05,q-值<0.05)，篩選獲得候選靶基因42個(gè)(見表4)

表2毛白楊響應(yīng)低溫逆境脅迫的lncrnas順式作用靶基因

表3毛白楊響應(yīng)低溫逆境脅迫的lncrnas反式作用靶基因

表4毛白楊響應(yīng)低溫逆境脅迫的lncrnas靶基因功能注釋

s2，對(duì)于候選基因，利用ncbi核酸數(shù)據(jù)庫(kù)(https://blast.ncbi.nlm.nih.gov/blast.cgi)進(jìn)行功能注釋(見表4)。利用agrigo(http://bioinfo.cau.edu.cn/agrigo/)對(duì)候選基因進(jìn)行g(shù)oterm富集分析(見表5)。

表5毛白楊響應(yīng)低溫逆境脅迫的lncrnas靶基因功能富集分析

s3，植物響應(yīng)逆境脅迫的lncrnas分組。根據(jù)逆境脅迫響應(yīng)lncrnas表達(dá)模式或靶基因功能富集結(jié)果，將響應(yīng)逆境脅迫的lncrnas進(jìn)行分組(見表6)。

表6毛白楊響應(yīng)逆境脅迫的lncrnas進(jìn)行分組

s4,利用meme(http://meme-suite.org/tools/meme)對(duì)分組后的響應(yīng)逆境脅迫的lncrnas特異序列模塊進(jìn)行富集分析。篩選參數(shù)設(shè)置如下：①每個(gè)lncrnas上至少預(yù)測(cè)到一個(gè)序列模塊；②預(yù)測(cè)序列模塊數(shù)量3-5個(gè)；③模塊長(zhǎng)度為6bp-15bp；④分布模式為正反義兩條鏈；⑤不允許軟件對(duì)序列進(jìn)行重排。第3分組共獲得富集的lncrnas序列模塊3個(gè)(見圖2)

s5，利用gomo(http://meme-suite.org//tools/gomo)，對(duì)富集的lncrnas特異序列模塊功能進(jìn)行功能預(yù)測(cè)。篩選參數(shù)設(shè)置如下：①顯著性參數(shù)設(shè)置為q-value<0.01；②計(jì)算次數(shù)5000次；③專一性(specificity)參數(shù)大于80％。序列模塊1獲得goterm注釋(見圖2)。

如圖2所示，序列模塊1期望值e-value＝2.7e-007；注釋條目：go：0003700，功能預(yù)測(cè)結(jié)果：轉(zhuǎn)錄因子活性(p-value＝2.6e-07；q-value＝1.1e-04；specificity＝83％)；跨膜受體蛋白酪氨酸激酶信號(hào)通路(p-value＝2.6e-07；q-value＝1.1e-04；specificity＝82％)；序列模塊2、序列模塊3無(wú)顯著富集注釋條目。期望值(e-value)表示由于隨機(jī)性造成獲得與數(shù)據(jù)庫(kù)比對(duì)結(jié)果的可能次數(shù)。期望值越小，發(fā)生這一事件的概率越低，比對(duì)結(jié)果越顯著。假設(shè)機(jī)率(p-value)代表給定原假設(shè)為真時(shí)樣本結(jié)果出現(xiàn)的概率。p值越小代表原假設(shè)不成立的概率越大。q值(q-value)代表經(jīng)過(guò)假陽(yáng)性率校正之后的p值，該數(shù)值越低則代表假陽(yáng)性率越小。特異性(specificity)代表該模塊的獨(dú)特性，數(shù)值越高則獨(dú)特性越強(qiáng)。

實(shí)施例2

對(duì)小葉楊1年生植株進(jìn)行高鹽處理(150mmnacl,6小時(shí))，提取其總rna用于對(duì)響應(yīng)滲透脅迫的lncrnas序列模塊進(jìn)行功能注釋。

s1，利用ribo-zerorrna試劑盒對(duì)核糖體rna進(jìn)行去除。利用smart試劑盒進(jìn)行鏈特異性cdna文庫(kù)構(gòu)建。利用illuminahiseqtm2500測(cè)序平臺(tái)完成cdna文庫(kù)測(cè)序，測(cè)序深度為10×。去除接頭以及冗余序列，通過(guò)cufflinks軟件拼接轉(zhuǎn)錄本，篩選長(zhǎng)度大于200nt、最小讀長(zhǎng)覆蓋率為5、cpcscore<0、cnciscore<0lncrnas，共獲得4241個(gè)。篩選差異倍數(shù)大于3且小于13(p-值<0.05,q-值<0.05)共13個(gè)lncrnas(見表7)。

表7小葉楊響應(yīng)高鹽逆境脅迫lncrnas

在13個(gè)高鹽脅迫響應(yīng)的lncrna上下游10kb范圍內(nèi)篩選順式作用靶基因，共獲得44個(gè)(見表8)。利用blast進(jìn)行序列互補(bǔ)計(jì)算，參數(shù)設(shè)置為e-value＝1e-10，identity＝90％以及利用rnaplex進(jìn)行熱力學(xué)上的互補(bǔ)計(jì)算，參數(shù)設(shè)置為e＝-70。共獲得反式作用靶基因59個(gè)(見表9)。根據(jù)差異表達(dá)靶基因篩選標(biāo)準(zhǔn)(p-值<0.05,q-值<0.05)篩選獲得候選靶基因73個(gè)(見表10)

表8小葉楊響應(yīng)高鹽逆境脅迫的lncrnas順式作用靶基因

表9小葉楊響應(yīng)高鹽逆境脅迫的lncrnas反式作用靶基因

表10小葉楊響應(yīng)高鹽逆境脅迫的lncrnas靶基因功能注釋

s2,對(duì)于候選基因，利用ncbi核酸數(shù)據(jù)庫(kù)(https://blast.ncbi.nlm.nih.gov/blast.cgi)進(jìn)行功能注釋(見表10)。利用agrigo(http://bioinfo.cau.edu.cn/agrigo/)對(duì)候選基因進(jìn)行g(shù)oterm富集分析(見表11)。

表11小葉楊響應(yīng)高鹽逆境脅迫的lncrnas靶基因功能富集分析

s3，植物響應(yīng)逆境脅迫的lncrnas分組。根據(jù)響應(yīng)逆境脅迫的lncrnas表達(dá)模式或靶基因功能富集結(jié)果，將響應(yīng)逆境脅迫的lncrnas進(jìn)行分組(見表12)。

表12小葉楊響應(yīng)逆境脅迫的lncrnas分組

s4,利用meme(http://meme-suite.org/tools/meme)對(duì)分組后的響應(yīng)逆境脅迫的lncrnas特異序列模塊進(jìn)行富集分析。篩選參數(shù)設(shè)置如下：①每個(gè)lncrnas上至少預(yù)測(cè)到一個(gè)序列模塊；②預(yù)測(cè)序列模塊數(shù)量3-5個(gè)；③模塊長(zhǎng)度為6bp-15bp；④分布模式為正反義兩條鏈；⑤不允許meme軟件對(duì)序列進(jìn)行重排。第3組共獲得富集的lncrnas序列模塊3個(gè)(見圖3)。

s5，利用gomo(http://meme-suite.org//tools/gomo)，對(duì)富集的lncrnas特異序列模塊進(jìn)行功能預(yù)測(cè)。篩選參數(shù)設(shè)置如下：①顯著性參數(shù)設(shè)置為q-value<0.01；②計(jì)算次數(shù)5000次；③專一性(specificity)參數(shù)大于80％。三個(gè)序列模塊分別獲得goterm注釋(見圖3)。

如圖3所示，序列模塊1期望值e-value＝1.5e-004；注釋條目：go：0009507，功能預(yù)測(cè)結(jié)果：葉綠體(p-value＝3.447e-07；q-value＝1.300e-02；specificity＝85％)；序列模塊2期望值e-value＝8.1e-003；注釋條目：go：0010287，功能預(yù)測(cè)結(jié)果：質(zhì)體球(p-value＝1.856e-07；q-value＝3.500e-03；specificity＝100％)；注釋條目：go：0009535，功能預(yù)測(cè)結(jié)果：葉綠體類囊體膜(p-value＝1.856e-07；q-value＝3.500e-03；specificity＝91％)。序列模塊3期望值e-value＝1.6e-004；注釋條目：go：0005739，功能預(yù)測(cè)結(jié)果：線粒體(p-value＝5.569e-06；q-value＝1.050e-02；specificity＝82％)；期望值(e-value)表示由于隨機(jī)性造成獲得與數(shù)據(jù)庫(kù)比對(duì)結(jié)果的可能次數(shù)。期望值越小，發(fā)生這一事件的概率越低，比對(duì)結(jié)果越顯著。假設(shè)機(jī)率(p-value)代表給定原假設(shè)為真時(shí)樣本結(jié)果出現(xiàn)的概率。p值越小代表原假設(shè)不成立的概率越大。q值(q-value)代表經(jīng)過(guò)假陽(yáng)性率校正之后的p值，該數(shù)值越低則代表假陽(yáng)性率越小。特異性(specificity)代表該模塊的獨(dú)特性，數(shù)值越高則獨(dú)特性越強(qiáng)。

以上所述僅是本發(fā)明的優(yōu)選實(shí)施方式，應(yīng)當(dāng)指出，對(duì)于本技術(shù)領(lǐng)域的普通技術(shù)人員來(lái)說(shuō)，在不脫離本發(fā)明原理的前提下，還可以做出若干改進(jìn)和潤(rùn)飾，這些改進(jìn)和潤(rùn)飾也應(yīng)視為本發(fā)明的保護(hù)范圍。