基于M iRf ilter系統的毛果楊 m iRNA預測

趙潔苑, 龔云路, 王翼飛

(上海大學 理學院,上海 200444)

基于M iRf ilter系統的毛果楊 m iRNA預測

趙潔苑, 龔云路, 王翼飛

(上海大學 理學院,上海 200444)

從參數訓練、參數范圍訓練、候選成熟體打分等方面改進 miRNA預測系統M iRfilter,使其適應擁有更長前體的植物miRNA的預測.預測毛果楊基因組上的 miRNA,并對系統進行精度檢驗.利用M iRfilter系統共預測出 3 860條候選 miRNA;在 110個正樣本中,正確識別 91條前體和 80條成熟體,前體預測精度為 82.73%,成熟體預測精度為 72.73%;在毛果楊第 4號染色體 (LG_Ⅳ)上得到的 1 968個負樣本中,有 12個數據可認為是 miRNA,假陽性率為0.61%.

植物 miRNA;M iRfilter;毛果楊基因組;一類分類法;K-最近鄰分類器 (KNN)

Abstract:Thispaper improves theMiRfilter system by parameter,range of parameter and score in mature miRNAs in order to predict miRNAs in plants.We predict miRNAs in Populus trichocarpa genome and use itspositive and negative samples to examine the accuracy of MiRfilter’sp rediction.M iRfilter predicts 3 860 Populus trichocarpa miRNA candidates in all.It correctly identifies91 p recursors and 80 matures in 110 positive samp les.Accuracies of the precursor and mature prediction reach 82.73%and 72.73%,respectively.We find 12 false positive miRNA s from 1 968 negative samp les in LG_Ⅳchromosome,and the false positive rate reaches0.61%.

Key words:plant miRNA;MiRfilter;Populus trichocarpa genome;one-class methods;K-nearest neighbor(KNN)

miRNA是一類長度約為 20～24 nt(少數小于20 nt)的內源性非編碼調控單鏈小分子 RNA,由一段具有發夾結構的單鏈 RNA前體 (pre-miRNA)剪切后生成.成熟的 miRNA 5′端為單磷酸基,3′端為羥基,通過與其靶 mRNA分子的 3′端非編碼區域 (3′-untranslated region,3′UTR)互補匹配來抑制該mRNA分子的翻譯.miRNA基因以單拷貝、多拷貝或基因簇等多種形式存在于基因組中,而且絕大部分定位于基因間隔區.miRNA在生物體中的基因表達具有進化保守性、時序性和組織特異性等特點,顯示了其在控制個體發育、決定細胞命運和分化中的特定功能[1].到 2009年 2月為止,miRBase(http:∥microrna.sanger.ac.uk/cgi-bin/sequences/)數據庫已發布了 8 619種 miRNA.

除一些基本特征外,植物 miRNA和動物miRNA有明顯區別[2-5].植物 miRNA前體比動物的更長、更復雜,通常為數百核苷酸;不同于動物 miRNA的加工來自于蛋白質編碼基因內含子,大部分植物miRNA前體產生于其自身的轉錄單元;植物 miRNA以單基因形式為多,基因簇內的 miRNA排列也相對較松散;植物 miRNA的靶序列還包括蛋白質編碼區[2].研究植物 miRNA的功能可以通過以下兩種方法:①導入抗miRNA的靶基因或上調 miRNA的表達,分析植物出現的表型變化;②用生化與分子生物學方法,如 cDNA末端快速擴增 (rapid amplification of cDNA ends,RACE)技術,測定 miRNA指導下靶mRNA剪切反應的精確位點[2].研究表明,植物miRNA主要通過以下 3種途徑來調節基因的表達:①通過堿基間互補配對直接結合于靶基因 mRNA上,從而導致靶基因的特異性剪切;②miRNA介導的翻譯阻遏;③miRNA介導的翻譯沉默[4].由此可見,植物miRNA的正常表達是植物正常生長發育所必需的,它所調節的靶基因控制著植物生長發育的各個方面,包括根、葉、花等形態發生、細胞分化、疏導組織形成等,也在調節植物對環境脅迫如干旱、鹽害和養分脅迫反應等方面起著重要的作用[2-3].

miRNA可以通過 cDNA克隆測序的方法加以識別.一些實驗室已經建立了不同組織、不同發育時期或不同生長條件下的 miRNA基因文庫.然而,真核生物組織中有些 miRNA的豐度較低,其表達又具有時序性和組織特異性,使得克隆法分離 miRNA十分困難.通過計算預測miRNA成為miRNA發現的一個有效方法,該方法以基因組序列和計算機程序鑒定為基礎進行科學預測和鑒定,彌補了 cDNA克隆測序方法中的不足.根據計算預測方法的本質,可分為以下 5種[6]:同源片段搜索方法、基于比較基因組學的預測方法[7-8]、基于序列和結構特征打分的預測方法[9]、結合作用靶標的預測方法、基于機器學習的預測方法[10-11].

MiRfilter系統是本實驗室自主開發的一個用于預測miRNA的自動化軟件,不依賴物種的同源性,僅根據物種已知 miRNA的固有信息進行預測,屬于基于序列和結構特征打分的預測方法,在病毒miRNA識別中具有良好的預測精度[12].然而,由于植物 miRNA的序列結構特征較病毒來說更為多樣和復雜,直接使用M iRfilter系統進行預測的效果并不理想.為了使MiRfilter系統能夠適應植物 miRNA的預測,本研究從參數訓練、參數范圍訓練、候選成熟體打分等方面對系統進行改進,并應用于毛果楊(Populus trichocarpa)的 miRNA預測.

1 材料與方法

1.1 M iRf ilter系統簡介

M iRfilter系統的預測步驟層次分明,整個過程分為 4個階段[12]:①在預測之前,對待測物種的基因組序列及已知 miRNA序列進行預處理;②根據預處理后的訓練集界定訓練參數和參數范圍;③對預測區域作二級結構模擬,從中提取合格的發夾結構;④根據訓練得到的參數范圍,從合格的發夾結構中篩選候選miRNA成熟體序列和前體序列.系統具體流程如圖1所示.

圖1 M iRf ilter系統流程圖Fig.1 Flow char t of the M iRf ilter system

1.2 M iRf ilter系統改進

1.2.1 參數訓練

在前體參數中,由于植物前體長度跨度較大,預測出的前體與miRBase數據庫給出的前體在序列兩端可能存在一定的差異,該差異會同時影響最小自由能的數值,使依賴于前體長度和自由能這兩個參數的MFEL[13]參數變化過大.用MFEL篩選前體,有時會將已知的 miRNA前體排除出去.因此,本研究去掉MFEL參數,另外增加 3個新參數,即定位前體(不包括發夾結構尾部的前體,如圖2所示)的長度、莖區配對個數以及它們的比值,使篩選不受前體序列兩端差異的影響,提高預測精度;還統計了兩條莖上配對堿基的個數,以此替代原來用于判斷合格發夾結構的標準——18 nt.

圖2 m iRNA前體二級結構示意圖(ptc-M IR156a發夾結構)Fig.2 Pre-m iRNA stem-loop(ptc-M IR 156a ha irp in structure)

在成熟體參數中,miRNA序列的首字母特征明顯,傾向于以 U開始[2],因此,本研究添加成熟體序列首字母參數,作為參與成熟體打分的一個變量;補充成熟體所在臂參數、發夾環與成熟體之間的配對堿基個數、成熟體與前體尾端之間的距離,這些參數與原有的發夾環與成熟體之間的距離參數一起界定候選成熟體在前體上的位置.

另外,本研究增加了一組成熟體互補序列參數,包括互補序列的長度、最大內環大小、平均內環大小、內環個數.這些參數能更好地反映成熟體互補序列本身的序列結構特征,以及成熟體序列與其在結構上的對稱性質.

最終,本研究確定以下三類描述 miRNA前體及成熟體特征的參數.

(1)前體參數:H1為前體序列長度 (Prelen);H2為發夾環大小 (Hplen);H3為發夾結構最小自由能(Energy);H4為前體莖區配對堿基個數 (paNum);H5為定位前體序列長度 (MarPrelen);H6為定位前體莖區配對堿基個數 (MarPaNum);H7為定位前體莖區配對堿基個數與其長度之比 (MarPAPL=MarPaNum/MarPrelen).

(2)成熟體參數:M1為成熟體序列長度(Marlen);M2為成熟體序列首字母 (Marst1);M3為成熟體序列中 C+G含量 (cgCon);M4為成熟體序列在發夾結構中的位置,即左臂或右臂 (Arm);M5為發夾環與成熟體之間的距離 (Dist1);M6為發夾環與成熟體之間的配對堿基個數(Dist1P);M7為成熟體與前體尾端之間的距離 (Dist2);M8為發夾結構莖區內成熟序列中不配對堿基的個數 (upNum);M9為發夾結構莖區內成熟序列的兩端處不配對堿基的個數(TerUpNum);M10為發夾結構莖區內成熟序列中最大內環的大小(InlpMax);M11為發夾結構莖區內成熟序列中內環的平均大小 (InlpAvg);M12為發夾結構莖區內成熟序列中內環的個數(InlpNum).

(3)成熟體互補序列參數:P1為成熟體互補序列長度(Parlen);P2為發夾結構莖區內成熟體互補序列中最大內環的大小 (Par InlpMax);P3為發夾結構莖區內成熟體互補序列中內環的平均大小(Par InlpAvg);P4為發夾結構莖區內成熟體互補序列中內環的個數(Par InlpNum).

1.2.2 參數范圍訓練

在對 1.2.1節參數進行范圍界定時,依然采用按總數據量的 3%刪除最大或最小參數值的方法.根據每個參數的實際意義選擇適合的刪除原則,使得到的范圍更具有針對性.具體可分為以下 4種情況:

(1)同時刪除最大值和最小值,得到一個范圍,適合一般參數,包括 H1,H2,H4,H5,H6,H7,M3,M7,P1.

(2)只刪除最大值,得到一個范圍,適合一般認為參數值越小越好的參數,包括 H3,M8,M9,M10,M11,M12,P2,P3,P4.

(3)取中位數,適合只需用一個均值描述整體的參數,包括M1(使用中位數是為了避免再次取整).

(4)將左臂和右臂作為兩類分別統計,適合在這兩類中范圍相差較大的參數,避免其中一類范圍擴大,包括 M5,M6.

表 1列出了毛果楊各參數的范圍.

1.2.3 M iRNA預測

本研究根據修改后的新參數調整了預測前體和成熟體的篩選標準.在最后預測成熟體的過程中,一個候選前體上可能會有多個符合標準的成熟體被保留下來.為此,引入一個打分機制,為每一個候選前體上預測出的成熟體打分,從中挑選出得分最佳的成熟體作為該前體的候選成熟體.具體打分方法采用最近鄰一類分類法.

一般而言,兩類分類法 (two-classmethods)需要考慮正樣本和負樣本兩組數據,通過一定的算法學習這兩類樣本,從而構建一個能夠區分它們的分類器.使用兩類分類法識別miRNA,是將已知的miRNA作為正樣本的同時,還需要人為地構造一組非miRNA的負樣本數據.但是負樣本的選擇具有一定的難度,如果選出的負樣本并不適合,就會顯著影響分類器的表現或者產生巨大誤差.另一方面,一類分類法 (one-class methods)只需要考慮目標類 (正樣本)的信息,就可以構建一個能夠識別目標類樣本并丟棄其他非目標類樣本的分類器,避免了人為構造負樣本數據[14].因此,在無法確定負樣本的情況下,采用一類分類法識別新的miRNA.

最近鄰一類分類法(one-class K-nearest neighbor classifier,OC-KNN)是一種修正了已知的最近鄰兩類分類法,使其只學習正樣本數據的分類方法[14].該算法存儲所有的訓練樣本(正樣本)y,將其作為鄰居集;對于一個給定的測試樣本 z,計算 z到鄰居集中所有鄰居 y的距離 d(z,y);將 k個最近鄰居距離的平均值作為 z的得分,當得分滿足一定條件時,將z歸為目標類.

在實際應用中,將已知的 miRNA作為訓練樣本y,將預測出的 miRNA作為測試樣本 z,每個樣本包含成熟體及其互補序列的序列參數和結構參數,即H7,M2,M3,M8,M9,M10,M11,M12,P1,P2,P3,P4共12個變量;取 k=1,保留滿足以下打分公式的測試樣本:

式中,d(z,y)采用歐拉距離,且變量在計算之前先標準化;N(z)為測試樣本 z所在的候選前體上所有預測出的miRNA的個數;閾值δ可根據已知 miRNA的得分進行選取.在每個候選前體上選取得分最低的,即與已知 miRNA相似度最高的成熟體,將其作為該候選前體的候選成熟體.表 2列出了毛果楊各條染色體選取的δ值.

表 1 毛果楊參數范圍Table 1 Ranges of Populus trichocarpa’s param eter s

表 2 毛果楊各染色體δ值Table 2 δof Populus trichocarpa’s each chrom osom e

1.3 數據集

本研究使用改進的MiRfilter系統在毛果楊基因組序列中預測miRNA,并根據預測結果對該系統的預測精度進行檢驗.毛果楊基因組 19對染色體 4.8億個堿基的測序工作已于 2004年 9月 21日完成,這是林木上第一個、植物上繼擬南芥和水稻之后第三個進行基因組測序的物種[15],其基因組序列 (版本 1.1)及相關注釋文件可從楊樹基因網站 JGI(http://www.jgi.doe.gov/poplar/)上獲得.毛果楊已知的 miRNA數據取自 miRBase數據庫 (2008年 8月).本研究保留 19對染色體中前體二級結構只含有一個發夾環、成熟體長度為 21 nt的 miRNA序列,共 110條前體上的 110條成熟體,將其作為正樣本數據;負樣本數據選用在毛果楊第 4號染色體 (LG_Ⅳ)的外顯子部分中預測出的 1 968條可能的成熟體.

1.4 評價標準

對于每一個測試樣本,只可能屬于以下 4種類型之一:正確識別的正樣本 TP、正確識別的負樣本TN、本來是負樣本卻被識別為正樣本 (假陽性樣本)FP、本來是正樣本卻被識別為負樣本 (假陰性樣本)FN.用 N表示樣本總數,Q表示總預測精度,QP表示正樣本的預測精度,QN表示負樣本的預測精度,FPR表示假陽性預測率,FNR表示假陰性預測率,MCC表示Matthew相關系數,分別定義如下[12]:

2 結果與討論

本研究在毛果楊非外顯子序列中預測出 3 860條成熟體,對應 3 860條前體;在 110個正樣本中正確識別出 91條已知前體和 80條已知成熟體,前體預測精度達 82.73%,成熟體預測精度達 72.73%,表 3為各條染色體的預測情況.在未被識別出的 30個miRNA中,有 16個 miRNA的前體已被預測出來,但由于存在得分更低的成熟體序列而被排除;有14個miRNA因沒有預測出其前體而被排除,表 4列出了未被成功識別出的毛果楊 miRNA.根據毛果楊第 4號染色體的閾值δ,對 1 968個負樣本數據進行篩選,最終有 12個數據被認為是 miRNA,假陽性率為 0.61%,具體假陽性數據見表 5.表 6為改進后的M iRfilter系統具體的預測精度.

表 3 毛果楊各染色體的預測結果和預測精度Table 3 Pred iction resultsand accuracy of Populus trichocarpa’s each chromosome

序列分析發現,病毒 miRNA之間的序列相似性很低,很少存在同源序列.對于很多病毒而言,它們只存在進化距離很遠的直系同源成員.類似的問題也發生在高等真核生物中.迄今為止,具有完整基因組序列且與擬南芥進化距離相對最近的物種是水稻,而水稻與擬南芥基因組早在 2億年前就已經分化.具有完整基因組序列且與人類進化距離相對最近的物種是黑猩猩,而黑猩猩與人類的基因組也早在 4百萬年前就已經分化[6].因此,不依賴序列保守性的 M iRfilter系統適用于各種生物的miRNA預測,它是發現非同源、物種特異 miRNA的有效途徑.

表 4 未成功識別的 m iRNATable 4 Un identif ied m iRNAs

雖然小分子的 miRNA可能以幾乎任意序列存在,其前體發卡環形的二級結構和它在前體的位置卻呈現出十分固定的特點,可以說,相對于序列間的相似性,miRNA更具備結構上的相似性[2].改進M iRfilter的時候,更注重加強 miRNA結構信息的描述.在 23個訓練參數中,有 4個序列參數和 19個結構參數;在 12個用于打分的變量中,也只有 2個與序列信息有關.對 miRNA結構特征的關注使得M iRfilter可以在不搜索同源片段的情況下仍然表具有良數好據的預測精度.

Table 5 False positivem iRNAs

表 6 M iRf ilter系統預測精度Table 6 Pred iction accuracy of M iRf ilter

雖然制定了嚴格的篩選標準來降低假陽性率,仍然會得到大量的預測結果,為了確認其中的 miRNA,需要進行實驗驗證.目前,miRNA的實驗檢測方法主要有:RNA印跡 (Northem blot)、實時熒光定量 PCR(real-time PCR)、芯片技術 (microarray)等[16].每種方法都有其優缺點,可以相互結合進行檢測.

miRNA的確定需要計算預測和實驗檢測共同完成.像M iRfilter這樣的計算預測工具可以從海量的數據中篩選出合理的潛在對象,彌補實驗方法效率低、成本高的缺點,已然成為實驗檢測不可或缺的前提條件.

[1] 華友佳,肖華勝.microRNA研究進展 [J].生命科學,2005,17(3):278-281.

[2] 金由辛.核糖核酸與核糖核酸組學[M].北京:科學出版社,2005:106-133.

[3] 王磊,范云六.植物微小 RNA(microRNA)研究進展[J].中國農業科技導報,2007,9(3):18-23.

[4] 李培旺,盧向陽,李昌珠,等.植物 microRNA s研究進展 [J].遺傳,2007,29(3):283-288.

[5] DUGAS D V,BARTEL B.M icroRNA regulation of gene expression in plants[J]. Current Opinion in Plant Biology,2004,7(5):512-520.

[6] 侯妍妍,應曉敏,李伍舉.M icroRNA計算發現方法的研究進展[J].遺傳,2008,30(6):687-696.

[7] GLAZOV E A,COTTEE P A,BARRISW C,et al.A microRNA catalog of the developing chicken embryo identified by a deep sequencing app roach[J].Genome Research,2008,18(6):957-964.

[8] RITCHIEW,THEODULE FX,GAUTHERETD.M ireval:a web tool for simple microRNA p rediction in genome sequences[J].Bioinformatics,2008,24(11):1394-1396.

[9] WANG X J,REYES J L,CHUA N H,et al.Prediction and identification of Arabidopsis thaliana microRNAs and theirmRNA targets[J].Genome Biology,2004,5(9):R65.

[10] XUE C H,L I F,HE T,et al.Classification of real and pseudo microRNA p recursors using local structuresequence features and support vector machine[J].BMC Bioinformatics,2005,6:310.

[11] HERTEL J,STADLER P F.Hairp ins in a Haystack:recognizing microRNA precursors in comparative genomics data[J].Bioinformatics,2006,22(14):197-202.

[12] 張玉濱.基于生物信息學方法預測 microRNA的研究[D].上海:上海大學,2007:36-49.

[13] 陳薇,譚軍,何晨.植物 miRNAs前體的生物信息分析[J].重慶郵電學院學報:自然科學版,2006,18(6):803-808.

[14] YOUSEFM,JUNG S,SHOWE L C,et al.Learning from positive examples when the negative class is undetermined-microRNA gene identification [J].Algorithms for Molecular Biology,2008,3:2.

[15] 甘四明,蘇曉華.林木基因組學研究進展 [J].植物生理與分子生物學學報,2006,32(2):133-142.

[16] 王旭丹.M icroRNA檢測方法的發展現狀[J].國際內科學雜志,2007,34(11):679-682.

(編輯:劉志強)

Pred iction of Populus trichocarpa m iRNAs w ith Im proved M iRf ilter System

ZHAO Jie-yuan, GONG Yun-lu, WANG Yi-fei
(College of Sciences,Shanghai University,Shanghai200444,China)

O 224

A

1007-2861(2010)04-0397-07

10.3969/j.issn.1007-2861.2010.04.014

2009-02-20

國家自然科學基金資助項目 (30871341);上海市重點學科建設資助項目 (S30104);上海市教委重點學科建設資助項目(J50101)

王翼飛 (1948～),男,教授,博士生導師,研究方向為計算分子生物學.E-mail:yifei_wang@staff.shu.edu.cn