六倍體小麥基因組注釋流程構建與優化①

祝海棟， 李瑞琳， 何小雨， 趙 丹， 韓鑫胤， 牛北方

(中國科學院 計算機網絡信息中心，北京 100190)

(中國科學院大學 計算機與控制學院，北京 100190)

小麥生產對保證糧食安全和農業可持續發展具有重要的現實意義，促進小麥的增產和品質改良成為當前小麥育種研究的前沿熱點. 為了培育具有優良性狀的新品種，首先要定位控制目標性狀的基因，因此建立一套完整準確的大尺度基因組注釋流程成為培育新品種過程中的難點之一. 基因組注釋主要包括基因識別和基因功能標注兩個方面[1]，本文的主要研究方向是基因識別，主要目標是準確定位基因位置及發現物種特異性基因.

近些年，基因組測序技術突飛猛進，其發展過程包含三個階段：1975年由桑格和考爾森開創的鏈終止法標志著第一代DNA測序技術的誕生，但測序成本高、通量低等缺點嚴重影響了其大規模的應用; 第二代測序建立在聚合酶鏈式反應擴增的基礎上，主要特點是為邊合成邊測序，測序結果讀長短、測序速度快、吞吐量大[2]; 第三代測序技術的核心是以單分子為目標，旨在解決第二代測序在準確性和組裝困難方面的問題.測序技術的高速發展，大大滿足了測序深度、重測序等大規?；蚪M的研究需求，改變了生命科學諸多領域的研究面貌，也給小麥等大尺度基因組的注釋研究奠定了重要基礎.

1 小麥基因組注釋流程研究現狀

傳統的基因注釋方法主要為數據庫比對，通過把基因組片段與已有的親緣物種基因數據庫比對，得到目標基因. 這種方法較為簡便，但具有三個明顯的缺點：一是對比速度慢，原因是該方法中需要與較多的數據庫進行比對分析，因此耗時長，尤其是用于小麥等較大基因組時該缺點更為明顯; 二是難以發現新的基因，由于依賴數據庫比對得到的基因都是目前相近物種中廣泛存在的基因，物種特有的基因不會被識別，造成注釋的不完整. 轉錄組測序可以全面快速的獲取物種在某一時期和特定組織中所有表達的基因序列，常被用于研究物種基因結構和基因功能[3]. 但是轉錄組分析軟件繁多，缺乏統一的選擇標準，且分析過程中涉及多個軟件配合完成，分析流程中不可避免地會存在軟件間銜接困難、格式轉換和大量數據重復讀寫等問題. 另外，由于各種軟件在內存、CPU等資源利用方面存在較大差異，且多數情況下生物信息學中的分析過程依賴于腳本生成的流程，沒有并行優化，因此資源利用率和分析效率較低. 針對上述問題，本文提出了整合基因組和轉錄組數據進行基因注釋的分析流程，以提高注釋的完整性和準確性.

2 實驗數據與測試環境

本次研究中使用的測試數據包括：科農9204小麥基因組組裝數據，數據大小14.24 GB; 二代轉錄組測序樣本77個，單樣本大小約為17 GB; 三代全長轉錄組測序樣本2個，單樣本大小約為40 GB. 測試環境為超級計算系統“元”. 其包含270臺計算節點，每節點采用2個Intel Xeon E5-2680V3處理器(2.5 GHz、12核)，單節點CPU計算能力 0.96Tflops，配備256 GB內存.操作系統為Linux version 2.6.32-358.el6.x86_64，CentOS release 6.4 (Final). 系統中配置Python、Perl、C++等基本編譯和運行環境.

3 基因組注釋分析軟件流程

本節分為3個部分，建立小麥基因組注釋分析流程，并對部分環節實現優化.

3.1 數據庫比對注釋

數據庫比對注釋是最傳統和最常用的注釋方式.其主要方法是把待注釋的基因組逐一與各個近親物種已有基因比對，獲取注釋結果. TriAnnot[4]是為解讀小麥基因組而開發的一個流程，集合了Blast、Repeat Masker等開源軟件，比對了NCBI[5]、TAIR10[6]等開放數據庫，對轉座子、編碼基因、非編碼序列、分子標記進行了多步的處理分析，可以得到比較完整的注釋結果. 因此，本文對TriAnnot注釋軟件進行優化并對科農9204小麥基因組進行初步注釋.

3.1.1 優化方法

為了提高注釋效率，本研究的主要貢獻是實現TriAnnot注釋軟件的優化，重點分為3個方面：單任務多實例并行優化、多核計算并行優化，多數據庫查找并行優化，下面給出具體的方法與實現.

首先，對TriAnnot注釋軟件的單任務多實例并行優化. 六倍體小麥基因組較大，每條染色體的平均長度接近700 MB，給序列比對帶來許多困難. 為了便于比對分析，在注釋過程中，必須把染色體切分成小的片段，本研究中選擇的切分大小為1 MB，切分時的保留的重復長度為50 KB. 切分后的每個片段即為每個實例，實例之間相對獨立. 為了提高注釋速度，本研究采用了多實例并行，即在每個時刻都有多個實例同時執行. 因為每個步驟的CPU和內存使用率各不相同，該優化策略可以實現資源的充分利用.

其次，實現TriAnnot注釋軟件的單任務多實例并行優化. 在實驗中，針對每個軟件的特點，本研究采用相應的調度方式優化. RepeatMasker通過相似性比對來識別重復序列，可以屏蔽序列中轉座子重復序列和低復雜度序列[7]. 本研究在流程中加入了RepeatMasker的多核心并行，可以根據機器硬件情況指定4至24核心實現并行運行，并且可以通過使用-qq指令加快比對效率.

最后，實現對TriAnnot注釋軟件的多數據庫查找并行優化. SIMSearch軟件通過使用多個同源數據庫進行序列比對找到親緣關系較近的基因序列. 為了加快同源基因的查找速度，研究采用了多數據庫并行的方案，把多個同源數據庫同時讀取到內存中，將每個基因片段在多個核心上與不同的數據庫進行對比.

3.1.2 軟件與數據庫

TriAnnot依賴的軟件及其下載地址如表1所示，數據庫及其下載地址如表2所示.

表1 TriAnnot主要依賴軟件

表2 TriAnnot主要數據庫

3.1.3 分析流程

該步驟的輸入為基因組裝得到的KN9204小麥基因序列文件. 六倍體小麥有21條染色體，此外還有少量未有效定位到染色體上的基因片段，在其中加入100個未知堿基標識“N”，構成未分組染色體，共22條fasta序列，每條序列單獨輸入.

TriAnnot軟件運行前需要下載完整的基因數據庫.主要參數包括：-W指定工作目錄，-s指定輸入的fasta文件，-t指定注釋流程xml文件，--type設置輸入為核酸，--maxlength設置最大序列長度，--splitseq設置超過最大長度的序列自動切分，--overlap設置切分時冗余長度.

軟件的輸出為gff文件，包含了詳細的內含子、外顯子、編碼區、轉座子等注釋.

3.2 轉錄組高通量測序

為了準確注釋物種特異性基因，本研究結合了轉錄組高通量測序數據，選取了苗期、孕穗期、7天、14天等不同時期的根、葉、穗等不同組織的樣本，測序深度約為30 X. 常用的轉錄組分析工具有HISAT、SATR、StringTie、Cufflinks等. 使用不同的分析工具和方法對分析結果的準確度和耗時影響較大，需要根據特定的數據集及特定的研究目標選擇合適的分析工具和方法. HISAT解決了轉錄組中僅有不連續的外顯子難以比對的問題，對比上代主流轉錄組比對工具Tophat效率高50倍，且內存需求更少[8]. StringTie繼承于Cufflinks，在準確性方面有了較大提升，且可以通過輸入數據庫比對注釋結果提高在已知基因區域的準確性，在組裝的過程中會計算每個基因及可變剪切的表達水平. 綜合以上優點，對于復雜的小麥基因組，本文使用HISAT[9]和StringTie[10]工具進行轉錄組組裝.主要分為以下四個步驟

(1) 建立HISAT2基因組索引. 轉錄組數據分析過程遇到的第一個問題就是，小麥上億條reads如何在保證錯誤率在可接受的范圍內，高效率地比對到基因組上. 針對上述問題，需要根據基因組序列使用hisat2-build命令建立索引.

(2) 將所有二代測序reads比對到基因組. 使用HISAT2利用基因組索引將高通量測序reads比對到基因組上. 參數-p指定并行核心數，-x指定索引位置，--dta為組裝提供錨點. 使用samtools將比對結果按染色體和起始位點排序.

(3) 使用StringTie對排序完成的reads進行組裝.不同組織中表達數據差異相對較大，比對到基因組的reads也各有不同，這些因素都會影響組裝的效率.

(4) 將所有轉錄本的組裝結果使用StringTie的merge模塊合并. 由于不同組織和不同時期表達的基因各不相同，為了獲取更加完整的注釋，需要對多個測序樣本合并. merge步驟可以跨多個測序樣本生成統一的轉錄本. 首先要創建一個文本文件，該文件包含所有轉錄本組裝結果路徑，文本的每行是單個樣本組裝結果文件路徑. 參數設置為：--merge指定使用合并模塊，-p指定并行核心數，輸入上述文本文件，即可得到最終的二代轉錄組組裝結果.

3.3 全長轉錄組單分子測序數據處理

二代測序可以準確地進行基因定量分析研究，但是受讀長限制，不能得到全轉錄本的信息. 全長轉錄組采用單分子實時測序技術，通過構建啞鈴型文庫，以環形方式循環測序[11]. 因此，通過全長轉錄組單分子測序可以不經過組裝，準確、直接地獲取整個轉錄本. 三代測序存在單堿基錯誤率較高的問題[12]，本研究使用PacBio公司發布的SMRTLINK Pipeline[13]，對三代測序得到的數據進行過濾與質量控制. 由于全長轉錄組測序成本相對較高，本次研究采取了常用的組織混合測序方式. 選取了葉、穗、幼葉、幼根四種組織混合，設置兩個生物學重復，共得到兩組測序數據. 數據處理過程主要分為以下3個步驟：

(1) 使用SMRTLINK進行三代測序數據的清洗.主要分為三個步驟，首先召回環形一致性序列，包括單堿基糾錯和序列過濾; 然后對序列分類，包括去除接頭、polyA尾部和串聯子; 最后進行迭代的聚類糾錯，主要是合并相似的序列，形成全長轉錄本. 該軟件提供了用戶可視化接口，安裝后使用瀏覽器訪問服務器地址的對應端口即可進入管理界面. 在管理界面中，使用“數據管理”選項導入原始測序結果文件，然后使用“SMRT分析”選項，選擇分析流程為“Iso-Seq”，設置相關參數，選取對應的樣本即可開始全長轉錄組的糾錯.在本次研究中，我們設置的參數主要有以下幾個：By Strand CCS：OFF; Maximum Dropped Fraction：0.8;Maximum Subread Length：15000; Minimum Predicted Accuracy：0.75; Minimum SNR：3.75; Polish CCS：ON;

其余參數均為默認值.

(2) 使用GMAP[14]比對全長轉錄本到基因組.GMAP具有一次對多條reads同時進行比對的優點，比對結果較為可靠，因此，本文采用GMAP將全長轉錄本比對到基因組上. 為了提高運算速度，GMAP比對階段對全長轉錄本序列進行數據分割，將分割后的多個數據進行并行處理. 首先使用gmap-build建立索引，由于小麥基因組較大，會自動使用長索引. 使用-D參數指定索引存儲位置，-d參數指定索引前綴，輸入基因組fasta文件即可開始建立索引，然后使用gmapl命令開始比對. 指定的索引存儲位置和前綴需與上述過程中對應參數相同，-B指定批處理個數，-t指定并行核心數，-f指定輸出格式，-O指定順序輸出. 使用samtools將bam文件按染色體和起始位點排序.

(3) 合并多個樣本的全長轉錄組結果. 合并時使用TAMA軟件，共分為兩個步驟. 首先根據比對到基因組上的位置情況合并可變剪切，然后合并多個測序樣本的轉錄本.

3.4 合并注釋結果

為了得到高質量的注釋結果，需對上述結果進行合并和過濾，在本次研究中我們開發了一個自動化合并注釋的軟件Annotator，該軟件包含的功能模塊有格式轉換，結果合并，去除重復序列，過濾可變剪切，根據證據支持評價可信度，編碼區預測，蛋白翻譯等多個步驟，最終生成gff注釋文件和轉錄本序列、編碼區序列、編碼蛋白序列. Annotator詳細流程如圖1所示. 合并過程分為以下5個步驟：

(1) 轉換結果文件為bed12格式. 本步驟調用了cufflinks[15]軟件的gffread模塊和bedops[16]軟件，將數據庫比對注釋得到gff文件和二代轉錄組組裝得到的gtf文件轉換為bed文件. bed格式使用單行定義單個基因，具有簡單易讀的特點.

圖1 注釋合并流程

(2) 合并數據庫比對注釋、二代轉錄組組裝、三代全長轉錄組結果. 本步驟使用了TAMA的merge模塊，生成含有全部基因的bed文件. 根據每個基因的支持證據的不同，分為高可信度基因和低可信度基因.

(3) 過濾重復的可變剪切. 由于測序誤差或reads組裝錯誤的不可避免，測序結果中可變剪切會出現許多冗余，因此需要對重復的可變剪切進行過濾. 過濾過程中，保留的優先級依次為全長轉錄組得到的可變剪切結果、數據庫比對注釋結果中的可變剪切，由于二代轉錄組組裝有更多的錯誤可能，其優先級最低.

(4) 預測所有基因的編碼區. 該步驟使用三種可能的翻譯方式分別將基因翻譯為氨基酸序列，取最長的序列，得到基因的編碼區.

(5) 翻譯編碼區序列. 根據注釋結果中的編碼區位置，將核酸序列翻譯為氨基酸序列，生成序列文件.

4 實驗結果與分析

經過優化，使用TriAnnot注釋科農9204基因組的重復序列時，速度提升達到60%，在1號染色體上的測試結果如圖2所示.

在轉錄組高通量測序過程中，建立索引過程輸入全基因組大小約為14 GB，耗時為8122秒，最大內存使用約為144 GB. 序列比對時，輸入共77個樣本，雙端測序的單個fastaq文件大小約17 GB，比對耗時約640秒. StringTie組裝輸入bam文件大小約為7 GB，耗時在10小時至24小時不等.三代全長轉錄組測序中單樣本bam文件為38 GB，運行時間約為149小時. 最終分別輸出高質量和低質量的全長轉錄組fasta序列. 最終得到的環形一致性序列質量分布如圖3所示，超過50%的序列質量均在0.99以上，可信度較高.

圖2 注釋耗時變化

圖3 全長轉錄組質量分布

本流程在六倍體小麥科農9204基因組上完成測試，共注釋出110 326個高可信度基因. 對比同源的中國春小麥基因組，其注釋包含107 891個高可信度基因[17]，其中有102 413個基因匹配，占中國春基因總數的94.9%，占科農9204基因總數的92.8%，具有高度一致性，這說明了本流程注釋結果具有較高的準確性.

5 結論與展望

本文提出了一種綜合運用數據庫比對、二代轉錄組高通量測序、三代全長轉錄組測序技術獲得準確注釋的分析流程，并獨立研發了注釋軟件Annotator. 隨后對流程中用到的部分軟件進行了優化，大大提高了注釋效率，為大尺度多倍體基因組提供了一個較為成熟的注釋軟件流程.

當前流程也存在一些問題：(1) 注釋速度仍然較慢.數據庫比對注釋過程是性能提升的主要瓶頸，仍需優化. (2) 成本較高. 注釋的準確性依賴于較高的測序深度，這會帶來成本的大幅提高，尤其是三代測序更為如此，這大大限制了該流程的廣泛應用.

因此，在未來的工作中將嘗試解決上述問題，以進一步優化整個流程. 針對注釋速度問題，可以對整個基因組進行更細粒度的并行處理，提升比對過程中的并行效率; 此外可以使整個比對過程均在內存中進行，避免中間結果寫入硬盤，減少不必要的時間開銷. 針對注釋中測序成本問題，可以在成本較高的三代全長轉錄組測序中采取多組織混樣測序的方案，選取最關注的組織和時期的樣本混合，通過單次測序降低成本.