基于RNA-seq的基因訓(xùn)練集構(gòu)建方法

段榮靜 劉金定

(1. 南京農(nóng)業(yè)大學(xué)領(lǐng)域知識(shí)關(guān)聯(lián)研究中心，南京，210095； 2. 南京農(nóng)業(yè)大學(xué)生物信息學(xué)中心，南京，210095)

引 言

從頭基因預(yù)測軟件(下簡稱預(yù)測軟件)常用于識(shí)別基因組上的編碼基因結(jié)構(gòu)，在基因組注釋中發(fā)揮著非常重要的作用。預(yù)測軟件使用前必須優(yōu)化其模型參數(shù)使之適應(yīng)被測序物種。除了幾個(gè)非監(jiān)督預(yù)測軟件外[1]，大多數(shù)預(yù)測軟件都屬于監(jiān)督訓(xùn)練軟件[2-5]，其訓(xùn)練過程必須建立在足夠數(shù)量的可靠基因結(jié)構(gòu)上。對于新測序的基因組而言，其研究基礎(chǔ)薄弱，往往缺乏足夠數(shù)量的編碼基因結(jié)構(gòu)用于訓(xùn)練，這大大限制了預(yù)測軟件的使用。因此，構(gòu)建高質(zhì)量的基因訓(xùn)練集是一個(gè)非常有價(jià)值的工作。

為新測序物種構(gòu)建基因訓(xùn)練集(下簡稱訓(xùn)練集)是一個(gè)非常耗時(shí)、費(fèi)力的過程，有時(shí)甚至比基因組測序本身還繁瑣。構(gòu)建訓(xùn)練集不僅需要克隆全長基因，而且還需要專家人工矯正基因結(jié)構(gòu)。另外，由于高表達(dá)和高保守的基因容易被發(fā)現(xiàn)和驗(yàn)證，所以人工克隆方法構(gòu)建的訓(xùn)練集通常會(huì)出現(xiàn)樣本不平衡的現(xiàn)象?？紤]到人工準(zhǔn)備基因訓(xùn)練集的繁瑣和缺陷，很多基因組注釋項(xiàng)目直接用其他近緣物種基因?qū)︻A(yù)測軟件進(jìn)行訓(xùn)練，或者直接使用其他物種模型參數(shù)跳過預(yù)測軟件訓(xùn)練環(huán)節(jié)。預(yù)測軟件性能依賴物種特異性的模型參數(shù)，而物種間的差異性又是必然存在的，因此這種做法必然會(huì)導(dǎo)致預(yù)測性能下降，尤其在親緣關(guān)系較遠(yuǎn)物種之間更是如此。為此，研究人員提出在基因組上識(shí)別保守基因結(jié)構(gòu)用于構(gòu)建訓(xùn)練集的方法[6-8]。這些方法定義了真核生物共同擁有的核心蛋白家族并提供了識(shí)別這些蛋白對應(yīng)基因結(jié)構(gòu)的計(jì)算方法。這些方法存在兩個(gè)缺陷：(1)真核生物核心基因過于保守,導(dǎo)致訓(xùn)練集不能特異優(yōu)化被測序物種的模型參數(shù)；(2)測序基因組不完整導(dǎo)致基因組上核心基因丟失，從而引發(fā)訓(xùn)練基因數(shù)量不足的現(xiàn)象。RNA-seq是快速獲得基因表達(dá)數(shù)據(jù)的有效方法，已成為基因組注釋項(xiàng)目中必不可少的步驟。理論上，RNA-seq可以獲得全部表達(dá)的基因序列，而且通過組裝讀段可以獲得全長序列[9-11]。將測序讀段比對到基因組上，可以獲得大量外顯子內(nèi)含子位置信息，進(jìn)一步組裝將得到外顯子-內(nèi)含子相互串聯(lián)的基因結(jié)構(gòu)。但是由于測序和組裝錯(cuò)誤的存在，導(dǎo)致獲得的基因結(jié)構(gòu)質(zhì)量不能滿足構(gòu)建訓(xùn)練集的要求[12]。另一方面，基于RNA-seq組裝的基因結(jié)構(gòu)無法直接給出編碼區(qū)的結(jié)構(gòu)以及起始密碼子和終止密碼子位置信息，因此進(jìn)一步降低了組裝結(jié)果的可用性。

針對基因組新測序物種基因匱乏，難以滿足預(yù)測軟件訓(xùn)練需要的現(xiàn)狀，本文提出了一種以新測序物種自身RNA-seq組裝為基礎(chǔ)，結(jié)合蛋白同源證據(jù)的可靠基因訓(xùn)練集構(gòu)建方法(Building reliable training gene set，BRTGS)。該方法利用蛋白同源證據(jù)排除不完整和錯(cuò)誤的RNA-seq組裝結(jié)構(gòu)，利用RNA-seq組裝結(jié)構(gòu)彌補(bǔ)蛋白同源證據(jù)難以獲得同源區(qū)外基因結(jié)構(gòu)以及起始密碼子和終止密碼子位置的缺陷。將該方法應(yīng)用于果蠅、人、擬南芥和水稻等幾種模式生物基因組上，實(shí)驗(yàn)結(jié)果表明該方法不僅可以為各種組裝水平的基因組構(gòu)建足夠數(shù)量的高質(zhì)量訓(xùn)練基因，而且訓(xùn)練后的預(yù)測軟件可以獲得很好的預(yù)測性能。

1 從頭基因預(yù)測軟件訓(xùn)練以及相關(guān)數(shù)據(jù)基礎(chǔ)

1.1 從頭基因預(yù)測軟件訓(xùn)練

從頭預(yù)測軟件訓(xùn)練是指利用物種已有的基因結(jié)構(gòu)信息優(yōu)化其模型參數(shù)，使之適應(yīng)該物種。雖然預(yù)測軟件訓(xùn)練效果會(huì)隨著訓(xùn)練基因的數(shù)量增加而逐步提高，但實(shí)踐表明當(dāng)訓(xùn)練基因數(shù)量超過1 000 h時(shí)，訓(xùn)練效果的提升不再明顯。因此，構(gòu)建基因訓(xùn)練集時(shí)更注重基因質(zhì)量而不是數(shù)量。高質(zhì)量的訓(xùn)練基因不僅要具有完整的外顯子-內(nèi)含子串聯(lián)結(jié)構(gòu)，而且要有正確的起始密碼子和終止密碼子位置信息(圖1)。RNA-seq組裝可以快速獲得大量基因外顯子-內(nèi)含子串聯(lián)結(jié)構(gòu)信息，但是這些結(jié)構(gòu)信息不僅存在組裝錯(cuò)誤，而且缺少起始密碼子和終止密碼子位置信息，因此不能直接用于構(gòu)建基因訓(xùn)練集。

圖1 真核基因結(jié)構(gòu)Fig．1 Structure of eukaryotic genes

1.2 RNA-seq組裝中的編碼基因結(jié)構(gòu)

目前公共數(shù)據(jù)庫中存在大量低N50基因組，甚至存在N50為10K的基因組。從頭預(yù)測軟件訓(xùn)練一般需200～1 000個(gè)可靠的基因結(jié)構(gòu)。為了模擬各種組裝水平的基因組，果蠅、人、擬南芥和水稻等幾種模式生物基因組被打碎隨機(jī)抽取序列構(gòu)成N50分別為2K，5K，10K，20K，40K，80K，160K，320K，640K和1M的片段化基因組(這些基因組在下文被分別引用為G2K，G5K，G10K，G20K，G40K，G80K，G160K，G320K，G640K和G1M)。在低組裝水平的基因組中大量基因組序列沒有被組裝到一起，導(dǎo)致跨度較大的基因被分散在多個(gè)基因組序列中而被重復(fù)計(jì)數(shù)，因此隨著N50下降，基因總數(shù)不斷提高(圖2(a))，但完整基因數(shù)量不斷下降(圖2(b))。在4個(gè)物種中，人類基因平均跨度最大，所以基因被碎片化的程度最高。在G2K基因組中，人類基因總數(shù)上升比例最高達(dá)到了460%，而完整基因數(shù)量下降最多達(dá)到42%。從數(shù)字上看，在G2K基因組中，4個(gè)物種分別保留了8 228，10 003，16 349和22 431個(gè)完整基因，這表明組裝水平較低的基因組中仍然保留著大量完整基因資源可用于構(gòu)建基因訓(xùn)練集。

圖2 N50對保留在片段化基因組中基因數(shù)量的影響Fig．2 Influence of genome N50 on the numbers of genes retained in fragmented genomes

雖然低水平組裝的基因組中仍然存在大量完整基因結(jié)構(gòu)，但需要確定RNA-seq組裝結(jié)果是否能夠提供足夠數(shù)量的正確外顯子-內(nèi)含子串聯(lián)結(jié)構(gòu)(簡稱基因結(jié)構(gòu))。鑒于生物學(xué)實(shí)驗(yàn)通常3個(gè)重復(fù)的要求，本文隨機(jī)地為每個(gè)物種準(zhǔn)備了3個(gè)RNA-seq數(shù)據(jù)。首先將每個(gè)RNA-seq數(shù)據(jù)比對到基因組上，然后組裝比對結(jié)果獲得基因結(jié)構(gòu)。一般而言，將多個(gè)RNA-seq組裝結(jié)果合并將得到更多基因結(jié)構(gòu)(基因結(jié)構(gòu)數(shù)量并不總是隨著RNA-seq數(shù)量增加而增長，當(dāng)RNA-seq測序總量足夠飽和時(shí)，基因結(jié)構(gòu)數(shù)將不再增長)，因此，本文進(jìn)一步利用CuffMerge將3個(gè)RNA-seq組裝結(jié)果合并在一起產(chǎn)生1個(gè)非冗余組裝結(jié)果，這樣每個(gè)基因組上得到4個(gè)RNA-seq組裝結(jié)果。由于絕大多數(shù)預(yù)測軟件只預(yù)測編碼區(qū)部分，所以本文只統(tǒng)計(jì)編碼區(qū)部分的結(jié)構(gòu)被完整且正確組裝的數(shù)量結(jié)果如表1所示。第1個(gè)樣品RNA-seq數(shù)據(jù)在G2K基因組上組裝獲得正確編碼基因結(jié)構(gòu)數(shù)量最少，組合3個(gè)RNA-seq數(shù)據(jù)在G1M上組裝獲得正確編碼基因結(jié)構(gòu)數(shù)量最多。由于這兩種情況分別代表了最差組裝結(jié)果和最優(yōu)組裝結(jié)果，因此本文后續(xù)研究只關(guān)注這兩種RNA-seq組裝結(jié)果。在G2K基因組上，第1個(gè)樣品RNA-seq組裝結(jié)果分別為4個(gè)物種提供了2 482，6 627，5 649 和7 928個(gè)正確編碼基因結(jié)構(gòu)。將3個(gè)樣品RNA-seq組裝結(jié)果合并起來，為G2K基因組分別提供4 933, 10 941, 8 764和11 846個(gè)正確編碼基因結(jié)構(gòu)。對于G1M基因組而言，無論是第1個(gè)樣品RNA-seq還是組合RNA-seq都組裝出非常多的正確編碼基因結(jié)構(gòu)。果蠅的第1個(gè)樣品RNA-seq數(shù)據(jù)在G1M上組裝的正確編碼基因結(jié)構(gòu)數(shù)量最少為3 866個(gè)，人類的組合RNA-seq數(shù)據(jù)在G1M上組裝的正確編碼基因結(jié)構(gòu)數(shù)量最多達(dá)到19 535個(gè)。預(yù)測軟件的訓(xùn)練只需幾百個(gè)正確編碼基因結(jié)構(gòu)，從RNA-sq組裝結(jié)果上看，即使在低組裝水平的基因組(G2K)上采用單個(gè)RNA-seq數(shù)據(jù)，也能組裝出近2 500個(gè)正確編碼基因結(jié)構(gòu)，這表明基于RNA-seq的組裝結(jié)果可為預(yù)測軟件訓(xùn)練提供充足的編碼基因結(jié)構(gòu)。

表1 RNA-seq組裝的完整編碼基因結(jié)構(gòu)的數(shù)量

注：每個(gè)物種的樣品1、樣品2和樣品3對應(yīng)的具體RNA-seq見3.1節(jié)實(shí)驗(yàn)數(shù)據(jù)。

2 基于RNA-seq組裝的基因訓(xùn)練集構(gòu)建方法

2.1 包含完整編碼區(qū)的RNA-seq轉(zhuǎn)錄本識(shí)別

首先用Tophat將RNA-seq短讀段序列比對到基因組上，然后用Cufflinks對比對結(jié)果進(jìn)行組裝獲得基因結(jié)構(gòu)，最后利用CuffMerge對每個(gè)基因組的3個(gè)RNA-seq組裝結(jié)果進(jìn)行合并獲得非冗余基因結(jié)構(gòu)。根據(jù)RNA-seq組裝的基因結(jié)構(gòu)獲得對應(yīng)的轉(zhuǎn)錄本序列后，用局部比對搜索工具X(Basic local alignment search tool X, BLASTX)將轉(zhuǎn)錄本序列比對到參考蛋白數(shù)據(jù)庫中(E-value設(shè)為1e-10)，只保留前50個(gè)洞(Gap)百分比小于3%且對齊區(qū)域覆蓋參考蛋白序列75%以上的比對結(jié)果用于后續(xù)分析。如果一個(gè)基因結(jié)構(gòu)符合上述比對結(jié)果，那么該基因結(jié)構(gòu)將被認(rèn)為可能包含完整編碼區(qū)。

2.2 起始密碼子和終止密碼子位置范圍計(jì)算

正確識(shí)別起始密碼子和終止密碼子位置對提升訓(xùn)練集質(zhì)量具有重要的意義，為此BRTGS進(jìn)一步利用同源參考蛋白序列初步確定起始密碼子和終止密碼子位置范圍，以提高后續(xù)起始密碼子和終止密碼子的預(yù)測可靠性,結(jié)果如圖3所示。轉(zhuǎn)錄本上的每個(gè)比對結(jié)果都可以推導(dǎo)一對起始密碼子和終止密碼子位置。在轉(zhuǎn)錄本對齊起始位置向5端延伸一定長度可以推導(dǎo)出起始密碼子位置，這個(gè)長度為參考蛋白序列上N端非對齊氨基酸長度的3倍。同樣，在轉(zhuǎn)錄本對齊終止位置向3端延伸一定長度可以推導(dǎo)出終止密碼子位置。一個(gè)轉(zhuǎn)錄本所有比對結(jié)果推導(dǎo)出的起始密碼子和終止密碼子位置收集在一起就構(gòu)成起始密碼子和終止密碼子位置范圍。

圖3 計(jì)算起始密碼子和終止密碼子范圍Fig．3 Calculation of site range of start codon and stop codon

2.3 轉(zhuǎn)錄本起始密碼子和終止密碼子位置確定

根據(jù)起始密碼子和終止密碼子確定方法，可將包含完整編碼結(jié)構(gòu)的轉(zhuǎn)錄本(Complete coding transcripts, CCTs)分為以下4類：

(1)CCT I :如果最佳比對結(jié)果的對齊區(qū)域完全覆蓋參考蛋白序列，對齊區(qū)域5端的3個(gè)堿基為起始密碼子，3端向后3個(gè)堿基為終止密碼子，那么這個(gè)轉(zhuǎn)錄本屬于CCT I，轉(zhuǎn)錄本上的對齊區(qū)域?yàn)榫幋a區(qū)。

(2)CCT II:如果最佳比對結(jié)果的對齊區(qū)域覆蓋參考蛋白序列的N端，對齊區(qū)域5端的3個(gè)堿基為起始密碼子，3端向后出現(xiàn)的第1個(gè)終止密碼子位于終止密碼子范圍內(nèi)，那么這個(gè)轉(zhuǎn)錄本屬于CCT II。轉(zhuǎn)錄本的編碼區(qū)從對齊區(qū)域開始到推導(dǎo)的終止密碼子上游位置結(jié)束。

(3)CCT III:如果最佳比對結(jié)果的對齊區(qū)域覆蓋參考蛋白序列的C端，3端向后3個(gè)堿基正好是終止密碼子，用ATGpr[13]計(jì)算推導(dǎo)出的起始密碼子位于起始密碼子位置范圍內(nèi)，那么這個(gè)轉(zhuǎn)錄本屬于CCT III。轉(zhuǎn)錄本的編碼區(qū)為從ATGpr計(jì)算的起始密碼子位置開始到終止密碼子上游位置結(jié)束。

(4)CCT IV：如果最佳比對結(jié)果的對齊區(qū)域既不能覆蓋參考蛋白序列的N端，也不能覆蓋C端，但是用ATGpr計(jì)算推導(dǎo)出的起始密碼子位于起始密碼子位置范圍內(nèi)，下游第1個(gè)終止密碼子位于終止密碼子范圍內(nèi)，那么這個(gè)轉(zhuǎn)錄本屬于CCT IV。轉(zhuǎn)錄本的編碼區(qū)為從ATGpr計(jì)算的密碼子位置開始到推導(dǎo)的終止密碼子上游位置結(jié)束。

在確定起始密碼子和終止密碼位置后，再次計(jì)算轉(zhuǎn)錄本編碼區(qū)的長度，如果編碼區(qū)長度大于最佳比對的參考蛋白序列長度130%，那么這個(gè)轉(zhuǎn)錄本將被淘汰。

2.4 同源聚類后排序生成基因訓(xùn)練集

在構(gòu)建的4類CCTs中，CCT I的起始密碼子和終止密碼同時(shí)具有表達(dá)、同源和預(yù)測證據(jù)支持，其可靠性最高；CCT II的起始密碼子同時(shí)有轉(zhuǎn)錄、同源和預(yù)測證據(jù)支持，終止密碼子只有表達(dá)和預(yù)測證據(jù)支持，其可靠性次之；CCT III的起始密碼子只有表達(dá)和預(yù)測證據(jù)支持，終止密碼子同時(shí)有表達(dá)、同源和預(yù)測證據(jù)支持，其可靠性再次；CCT IV的起始密碼子和終止密碼子都只有表達(dá)和預(yù)測證據(jù)支持，其可靠性在4類CCTs中最低。構(gòu)建基因訓(xùn)練集不僅要考慮基因結(jié)構(gòu)可靠性，而且還要考慮基因的同源性，應(yīng)盡量避免使用過多同源基因構(gòu)造訓(xùn)練集，以防止產(chǎn)生訓(xùn)練偏好性。圖4給出了根據(jù)可靠性和同源性排序生成基因訓(xùn)練集的示意圖，由圖4可知，BRTGS對4類CCTs進(jìn)行排序處理的過程為：

圖4 根據(jù)可靠性和同源性排序生成基因訓(xùn)練集示意圖Fig．4 Demonstration of building training gene set according to reliability and homology

(1)用OrthoMcl[14]對所有的CCT s進(jìn)行同源聚類，將同源的CCTs聚類到同1個(gè)同源分組中，然后對每個(gè)同源分組內(nèi)的CCTs按照“CCT I > CCT II > CCT III > CCT IV”進(jìn)行可靠性排序，以確保在每個(gè)同源分組中高可靠的CCTs被優(yōu)先選擇用于構(gòu)建訓(xùn)練集。

(2)對同源分組按照包含CCTs類別進(jìn)行排序。排序后包含CCT I的同源分組將排在最前面，如果同源分組中沒有CCT I，那么包含CCT II的同源分組排在最前面，以此類推。如圖4所示,在第1輪排序中Group 4不包含CCT I被排到最后；從4個(gè)Group中取走第1個(gè)CCT后進(jìn)入第2輪排序，Group 1具有CCT I被排在最前，Group 4和Group 3都沒有CCT I和II，所以被排在最后；再次從4個(gè)Group中取走第1個(gè)CCT后進(jìn)入第3輪排序，Group 2不包含CCT III被排到最后。

(3)根據(jù)步驟(2)對同源分組排序的結(jié)果，依次將每個(gè)同源分組中的第1個(gè)CCT移到訓(xùn)練集中。

(4)重復(fù)步驟(2，3)，直到每個(gè)分組中的CCT都轉(zhuǎn)移到訓(xùn)練集中。

經(jīng)過上述排序處理后，訓(xùn)練集中的同源基因?qū)⒈淮蛏?。由于?xùn)練時(shí)，BRTGS只取訓(xùn)練集中前1 000個(gè)基因用于訓(xùn)練，因此不僅可以避免同一同源分組中的CCTs被過度采用的現(xiàn)象，而且確保每個(gè)同源分組中高可靠性CCTs被優(yōu)先采用。

3 實(shí)驗(yàn)與結(jié)果分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

果蠅、擬南芥、人和水稻4個(gè)模式物種基因組數(shù)據(jù)來自于Ensemble(http://asia.ensembl.org)，版本號分別為BDGP6.31，TAIR10.31，GRCh38.p2和IGSP-1.0.31。為了模擬不同組裝水平的基因組，按照N50定義將每個(gè)物種基因組隨機(jī)打碎。在產(chǎn)生低于N50值的序列時(shí)，如果序列長度小于100 bp時(shí)，那么這個(gè)序列將會(huì)被丟棄。在產(chǎn)生大于N50值的序列時(shí)，最長序列不超過1 Mb(G640K和G1M的最長序列分別為2 Mb和8 Mb)。4個(gè)物種RNA-seq數(shù)據(jù)都來源于NCBI的SRA(http://www.ncbi.nlm.nih.gov/sra)數(shù)據(jù)庫，訪問編號分別為“SRR3091999，SRR3091976，SRR3138705”，“ERR754089，SRR1792926，ERR1109342”，“SRR3184282，SRR3184286，SRR3184279”和“SRR2156305, SRR2047477, SRR2048540”。每個(gè)物種的3個(gè)RNA-seq數(shù)據(jù)都分別被依次標(biāo)為樣品1、樣品2和樣品3。RNA-seq比對和組裝分別采用Tophat和Cufflinks實(shí)現(xiàn)，參數(shù)使用默認(rèn)值[15-17]。參考蛋白序列來自NCBI-RefSeq數(shù)據(jù)庫，并去除掉來自這4個(gè)物種的蛋白序列。

3.2 實(shí)驗(yàn)方案

為了評估BRTGS的性能，G2K上最差RNA-seq組裝結(jié)果和G1M上最佳RNA-seq組裝結(jié)果(合并了3個(gè)RNA-Seq的組裝結(jié)果)分別用于構(gòu)建基因訓(xùn)練集。作為比較，核心真核基因定位方法(Core eukaryotic genes mapping approach, CEGMA)(v2.5)、通用單拷貝直系同源基因檢測(Benchmarking universal single-copy orthologs, BUSCO)(v1.22)和轉(zhuǎn)錄本解碼器(TransDecoder)(3.0)也用于構(gòu)建基因訓(xùn)練集。CEGMA和BUSCO利用保守的核心基因構(gòu)建訓(xùn)練基因集，因此直接用在G2K和G1M基因組上。TransDecoder會(huì)為每個(gè)轉(zhuǎn)錄本序列報(bào)道多個(gè)編碼區(qū)，這里只保留標(biāo)志為“Complete”且最長的結(jié)果作為基因編碼結(jié)構(gòu)。本文將對4種方法產(chǎn)生的訓(xùn)練基因數(shù)量和質(zhì)量進(jìn)行比較。

預(yù)測軟件的訓(xùn)練效果也是構(gòu)建基因訓(xùn)練集優(yōu)劣的重要評價(jià)指標(biāo)。本文用訓(xùn)練集中的前1 000個(gè)基因(如果不足，則訓(xùn)練集中的全部基因都用于訓(xùn)練)訓(xùn)練Augustus(v3.2.1)，并用其訓(xùn)練后的預(yù)測性能評估訓(xùn)練集的訓(xùn)練效果。測試集來自基因組上和訓(xùn)練集不重疊的參考基因。本文將用訓(xùn)練后Augustus在核酸水平、外顯子水平以及基因水平上的預(yù)測性能評估訓(xùn)練效果，訓(xùn)練和評估過程完全按照Augustus操作手冊進(jìn)行。

3.3 實(shí)驗(yàn)結(jié)果分析

首先對4種方法構(gòu)建訓(xùn)練基因數(shù)量情況進(jìn)行比較，結(jié)果如表2所示?？傮w上CEGMA和BUSCO構(gòu)建的基因數(shù)量明顯偏少，主要原因是CEGMA和BUSCO構(gòu)建基因訓(xùn)練集被限制在其定義的真核(或者某個(gè)譜系)核心基因集上。這些核心基因數(shù)量本來就偏少，在低水平組裝的基因組上構(gòu)建的訓(xùn)練基因則進(jìn)一步減少。BRTGS和TransDecoder構(gòu)建訓(xùn)練基因數(shù)量偏多，主要原因是BRTGS和TransDecoder構(gòu)建基因訓(xùn)練集來自于本物種RNA-seq測序數(shù)據(jù)。RNA-seq測序深度高，其組裝結(jié)果基本可以覆蓋基因組上全部基因(除了少數(shù)極低豐度表達(dá)或不表達(dá)的基因外)，因此，構(gòu)建的訓(xùn)練基因數(shù)量明顯比CEGMA和BUSCO多。和TransDecoder相比，BRTGS進(jìn)一步采用蛋白同源證據(jù)檢驗(yàn)，因此構(gòu)建的基因數(shù)量比較少。盡管如此，BRTGS在每個(gè)基因組上仍然構(gòu)建了超過1 600個(gè)基因的訓(xùn)練集。BRTGS構(gòu)建的訓(xùn)練基因并沒有完全來自CCT I，這說明其構(gòu)建的基因訓(xùn)練集并不完全來自跨物種保守基因。在擬南芥和人的基因組上大約64%和61%來自CCT I，在果蠅和水稻上只有26%和39%來自CCT I。此外，在BRTG構(gòu)建的基因集中，只有不足5%的基因被同源聚類到一起，因此訓(xùn)練集不具有同源偏好性。

表2 4種方法在不同基因組上構(gòu)建的訓(xùn)練基因數(shù)量

其次，對4種方法構(gòu)建訓(xùn)練基因質(zhì)量進(jìn)行比較,結(jié)果如表3所示。以4個(gè)物種基因組注釋的基因結(jié)構(gòu)為參考，統(tǒng)計(jì)4種方法構(gòu)建的基因訓(xùn)練集的準(zhǔn)確性。在外顯子水平上，當(dāng)編碼外顯子兩側(cè)邊界和參考基因外顯子完全重疊時(shí)，該外顯子結(jié)構(gòu)被認(rèn)為正確，否則為錯(cuò)誤。在基因水平上，構(gòu)成基因編碼結(jié)構(gòu)的外顯子和內(nèi)含子與參考基因完全一致時(shí)，該基因結(jié)構(gòu)被認(rèn)為正確，否則為錯(cuò)誤。準(zhǔn)確率為構(gòu)建訓(xùn)練集中正確外顯子(或基因)數(shù)與外顯子(或基因)總數(shù)的比值。在所有測試的基因組上，BRTGS構(gòu)建的基因結(jié)構(gòu)可靠性最高。BRTGS平均外顯子準(zhǔn)確性達(dá)到93%，分別比CEGMA，BUSCO和TransDecoder高16%，24%和20%；平均基因結(jié)構(gòu)準(zhǔn)確性為77%，分別比CEGMA，BUSCO和TransDecoder高40%，40% 和 31%。BRTGS構(gòu)建的基因結(jié)構(gòu)具有最高的可靠性，主要原因是其利用了表達(dá)、同源和預(yù)測3個(gè)方面證據(jù)確?；蚪Y(jié)構(gòu)的完整性和正確性。而CEGMA和BUSCO只利用跨物種保守基因的蛋白同源證據(jù)構(gòu)建基因訓(xùn)練集，難以準(zhǔn)確捕捉同源區(qū)域外的外顯子結(jié)構(gòu)以及起始密碼子和終止密碼子位置。雖然TransDecoder也利用了RNA-seq組裝數(shù)據(jù)(即RNA-seq表達(dá)組裝證據(jù))，但沒有對組裝結(jié)構(gòu)正確性以及完整性進(jìn)行檢查，因此基因結(jié)構(gòu)的可靠性也比BRTGS低。

最后，用訓(xùn)練效果(即預(yù)測軟件訓(xùn)練后的預(yù)測性能)對4種方法構(gòu)建的訓(xùn)練集進(jìn)行比較，結(jié)果如表4所示。在每個(gè)測試的基因組上，基于BRTGS基因集訓(xùn)練的Augustus，在核酸、外顯子和基因水平上獲得的靈敏度和特異度都比其他3種方法高(詳細(xì)數(shù)據(jù)未列出)?；?qū)用娴撵`敏度和特異度是檢查訓(xùn)練效果的最重要指標(biāo)。用BRTGS基因集訓(xùn)練的Augustus，平均基因靈敏度為33%，分別比CEGMA，BUSCO和TransDecoder高9%，12%和9%；平均基因特異度為22%，分別比CEGMA，BUSCO和TransDecoder高6%，8%和4%?？傮w上，BRTGS和TransDecoder構(gòu)建基因集的訓(xùn)練效果比CEGMA和BUSCO好。BRTGS構(gòu)建的基因訓(xùn)練集具有最好的訓(xùn)練效果主要得益于其充足的訓(xùn)練基因數(shù)量和可靠的基因結(jié)構(gòu)。CEGMA和BUSCO構(gòu)建基因集的訓(xùn)練效果差的原因來自兩個(gè)方面：(1)訓(xùn)練基因數(shù)量偏少、可靠性偏低；(2)在跨物種高度保守的核心基因上訓(xùn)練，難以獲得基因組新測序物種自身特性。

表3 4種方法構(gòu)建的基因訓(xùn)練集質(zhì)量

表4 基于4種方法基因訓(xùn)練集的預(yù)測軟件平均預(yù)測性能

4 結(jié)束語

長期以來，為從頭預(yù)測軟件準(zhǔn)備訓(xùn)練基因集需要專家人工矯正基因結(jié)構(gòu)，這是個(gè)耗時(shí)、費(fèi)力的工作。高通量測序技術(shù)能夠快速獲得大量基因表達(dá)序列，尤其是三代測序技術(shù)的應(yīng)用能夠產(chǎn)生更長甚至全長mRNA序列，這為構(gòu)造基因訓(xùn)練集進(jìn)一步提供了十分有價(jià)值的序列數(shù)據(jù)。然而從頭預(yù)測軟件訓(xùn)練不僅需要基因序列數(shù)據(jù)，而且還要基因元件(如起始密碼子、終止密碼子、內(nèi)含子以及外顯子等)在基因組上的位置信息。本文提出了一種利用物種自身RNA-seq組裝，結(jié)合同源證據(jù)構(gòu)建基因訓(xùn)練集的方法，該方法克服了單獨(dú)利用蛋白同源證據(jù)或RNA-seq表達(dá)證據(jù)構(gòu)建基因訓(xùn)練集的缺陷，提高了構(gòu)建基因訓(xùn)練集的質(zhì)量。通過比較分析，該方法構(gòu)建的基因訓(xùn)練集不僅具有更高質(zhì)量的基因結(jié)構(gòu)，而且具有很好的訓(xùn)練效果。值得注意的是，由于基因編碼區(qū)決定了基因功能，因此本文提出的方法只關(guān)注訓(xùn)練基因編碼區(qū)結(jié)構(gòu)的準(zhǔn)確性。事實(shí)上，大量研究表明，基因的非編碼區(qū)(Untranslated region，UTR)也非常重要，比如microRNA經(jīng)常結(jié)合在3UTR上實(shí)現(xiàn)對基因的調(diào)控。下一步研究重點(diǎn)是進(jìn)一步優(yōu)化本文提出的方法，利用RNA-seq測序數(shù)據(jù)尤其是三代測序長序列數(shù)據(jù)為基因訓(xùn)練集補(bǔ)充結(jié)構(gòu)完整的UTR。