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

背景技術：

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

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

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

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

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

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

步驟s3，按照靶基因表達模式進行植物響應逆境脅迫的lncrnas分組；

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

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

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

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

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

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

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

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

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

更優(yōu)選地，期望值(e-value)表示由于隨機性造成獲得與數(shù)據(jù)庫比對結果的可能次數(shù)。期望值越小，發(fā)生這一事件的概率越低，比對結果越顯著。假設機率(p-value)代表給定原假設為真時樣本結果出現(xiàn)的概率。p值越小代表原假設不成立的概率越大。q值(q-value)代表經(jīng)過假陽性率校正之后的p值，該數(shù)值越低則代表假陽性率越小。特異性(specificity)代表該模塊的獨特性，數(shù)值越高則獨特性越強。本發(fā)明所述的植物響應逆境脅迫的lncrnas序列模塊功能注釋方法中，所述注釋方法中篩選參數(shù)設置如下：①顯著性參數(shù)設置為p-value<0.01并且q-value<0.01；②計算次數(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注釋轉錄因子活性(見圖2)。

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

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

附圖說明

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

圖2毛白楊響應低溫lncrnas序列模塊功能注釋結果圖；

圖3小葉楊響應高鹽脅迫lncrnas序列模塊功能注釋結果圖。

具體實施方式

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

以下通過具體實施例進一步對本發(fā)明的技術方案進行說明，應理解以下僅為本發(fā)明的示例性說明，并不用于限制本發(fā)明權利要求的保護范圍。

實施例1

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

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

表1毛白楊響應低溫逆境脅迫的lncrnas

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

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

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

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

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

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

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

表6毛白楊響應逆境脅迫的lncrnas進行分組

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

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

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

實施例2

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

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

表7小葉楊響應高鹽逆境脅迫lncrnas

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

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

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

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

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

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

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

表12小葉楊響應逆境脅迫的lncrnas分組

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

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

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

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