通過區域篩選和簡潔de Bruijn圖比對高重復短序列

黎 瑤，鐘 誠

(廣西大學 計算機與電子信息學院 廣西高校并行分布式計算技術重點實驗室，南寧 530004)

E-mail:1176423710@qq.com

1 引 言

序列比對是生物大數據分析應用中的基礎性工作,有利于后續的基因組信息解讀.短序列比對(short-read alignment)是序列比對中的一種，它將測序平臺產生的大量的短序列與參考基因組進行對準定位[1].細菌、哺乳動物等物種中存在著許多重復子序列，這給序列比對提出了新的挑戰.Illumina機器產生的高通量測序序列相對較短(100～300bp)、錯誤率相對較低(1%～2%)、但重復率高的特點使得序列在很多位置都將被“匹配”[2].如何找到唯一的最佳匹配相當困難，這使得序列比對更具挑戰性.處理高通量測序數據使用的關鍵數據結構之一是de Bruijn圖[3].De Bruijn圖在基因組裝[4]、轉錄組裝配[5]、宏基因組學[6]、變體識別和結構變異檢測[7]、序列比對[8]等領域得到了應用.

已有的大多數序列比對方法都是基于種子-擴展(seed-and-extending)策略進行搜索比對[9-13].基于種子-擴展策略的序列比對主要有4個步驟：種子生成、種子匹配、種子擴展(擴展每個匹配的種子)以及最終對準(測序序列與參考基因組的最終映射定位).基于種子-擴展策略的序列比對面臨的主要瓶頸之一是如何處理基因組中的重復區域(例如：人類基因組含有超過50%的重復子序列[14]).當將測序序列(read)與包含許多“重復”區域的參考基因組進行比對時，將會出現許多種子在參考基因組中有多個命中(hit)的情況，從而導致所有種子的擴展計算成本非常高.近幾年來，人們提出了啟發式的過濾命中[12,15]、限制每個種子的命中數[9,11]、在擴展前合并命中[9-10,13]等策略，以期望有效地處理種子命中.然而，這些策略均是獨立地處理每個命中.評估、排序和處理大量命中的成本高.在評估、排序和處理步驟完成之后，由于子序列的重復可能仍有許多候選位置留下來需進行擴展，這導致將大量測序序列與參考基因組每個候選位置逐個對準的成本仍然較高.

為解決大量的種子命中和擴展候選位置代價高的問題，BGREAT[16]和deBGA[17]采用早期停止機制：當種子不匹配的數量超過某個閾值時，則停止探索de Bruijn圖的節點.采用此機制可以減少搜索空間，但可能無法返回最佳比對結果.Heydari等人提出一種將大量的短序列與參考基因組生成的de Bruijn圖進行比對的算法BrownieAligner[18].該算法計算de Bruijn圖的每個分支對準計分上限，直到無法改進目前找到的最佳解決方案時，則丟棄分支.在搜索比對de Bruijn圖過程中，BrownieAligner算法采用的深度優先搜索可能無法在找到最佳路徑之前訪問更多的結點，這將導致正確比對百分比不夠高.針對將含有重復子序列的測序短序列與參考基因組進行比對的問題，本文采用基于空位種子查找策略的區域選擇方法，依據種子預定義由給定的shape布局中生成的關鍵字(key)建立hash索引，通過查找hash索引篩選候選位置以減少種子的候選位置個數，并運用Hough變換將種子命中進行分組操作，將其聚集為粗對準的形式，然后采用簡潔de Bruijn圖結構壓縮存儲和索引k-mer，以使得在確保正確比對百分比的同時，減少比對時間和所需的存儲空間.

2 算 法

2.1 算法思想

細菌、哺乳動物等物種中存在大量重復的DNA子序列，尤其是在人類基因中超過50%的位置出現重復子序列.為能夠使用盡可能少的存儲空間，快速準確地將存在較多重復子序列的測序序列和參考基因組進行比對，本文給出一種改進的短序列精確比對算法.算法的關鍵是要快速地確定出參考基因組上可能對準序列的區域.區域選擇依賴于在測序序列和參考基因組之間找到種子，并將其聚類為粗對準形式.算法的第一階段，采用基于空位種子(gapped seeds)[19]搜索策略的區域選擇方法來篩選比對的候選位置，以盡可能地減少候選位置[20]，從而有效地減少搜索空間；然后，運用Hough變換將種子命中聚集為粗對準形式.第二階段，生成參考基因組對應的簡潔de Bruijn圖，將測序序列與生成的de Bruijn圖進行映射定位.第三階段，驗證粗對準中的候選位置子序列，以生成最終的序列對準結果.

2.2 算法描述

為便于算法描述，首先說明算法使用的參數的含義.其中全局參數bin_id表示所在區域下標，bin_value為所在區域命中個數，bin_value_size表示區域命中個數的規模，|r|為序列r的長度，|k_index|表示索引的規模，hit_counts表示hit計數，total_num_hits為總的種子命中數量，G表示參考基因組.局部參數x表示測序序列read坐標，y表示參考基因組坐標，l_local表示參考基因組與序列read的相對位置坐標.

為方便理解比對算法，下面介紹一些相關概念和知識.空位種子搜索策略由Burkhardt等人提出，其中q-gapped seeds表示非連續shape(模式)的長度為q的字符子串[19].例如，字符串ACAGCT中的3-gapped seeds的shape“1101”是ACG，CAC和AGT，其中“0”被定義為“do not care(DC)”字符的位置.文獻[19]的研究表明，gapped seeds允許快速地對不精確匹配進行敏感查找，在種子的預定義“do not care(DC)”堿基位置中允許變化(包含插入和缺失).空位種子(gapped seeds)查找策略是基于Levenshtein距離的種子查找方法，在序列比對中它比核酸和氨基酸序列分析中的連續種子(連續q-grams，連續k-mers)更敏感[19].

通過構造參考基因組的hash索引來實現基于Levenshtein距離的空位種子查找策略，其中q-gapped 種子的子序列是由shape布局(layout)中生成的關鍵字進行hash索引的，該索引僅包含用于構造關鍵字的堿基，而不包含DC堿基(見圖1中的“Do not care”堿基).

圖1 用于索引構建和查找的插入空位的種子結構示例Fig.1 Example of gapped seeds structure for index creating and search

在索引的查找階段，為每個shape構造多個關鍵字并用于檢索.對于每個DC堿基，構造3個查找關鍵字：

1)DC堿基匹配或誤配所有種子位置的關鍵字.見圖1中的匹配或誤配種子.

2)包含DC堿基的關鍵字.該關鍵字在指定位置最多允許一個刪除(因為indel的長度一般為1bp).見圖1中的刪除種子.

3)不包含DC堿基及其后堿基的關鍵字.該關鍵字最多允許一個插入堿基和一次匹配或者不匹配堿基.見圖1中的插入種子.

于是，每個shape建立d3個關鍵字，d表示shape中DC堿基的個數.基于Levenshtein距離的最佳shape的計算困難性[21,22]和對一系列組合的經驗評估，本文選擇以下的互補shape組合進行區域選擇：1111110111111(6-1-6 shape)和11110111101111(4-1-4-1-4 shape)，其中1表示包含的堿基，0表示DC的位置.本文為每個shape分別建立自己的hash索引.在參考基因組中每個種子出現的位置，搜索找到相關的兩個互補shape，并將所有的命中執行下一步的分組操作，從而將候選位置上的種子聚類得到粗對準形式.

本文運用Hough變換方法獲得分組種子命中對角線，以進一步完成區域選擇(見圖2)，其中Hough變換用于檢測諸如線、圓和橢圓的形狀的處理，它定義從圖像點到累積器(accumulator，基地址寄存器)的空間稱為Hough空間的映射.對于線檢測的情形，若笛卡爾空間中給定的一組點(一組種子)是共線的，則它們的關系可以用具有公共斜率m和截距c的線性方程表示：

y=mx+c

(1)

其中x和y是直角坐標系中的坐標(x,y)的元素.使用Hough變換確定描述給定點集(種子集)所在直線的參數m和c.本文可使用線性回歸技術來確定給定點集的共線問題，則式(1)可轉換為雙參數空間：

c=-mx+y

(2)

給定直角坐標系中的點(x,y)，其參數空間表達被定位為一條直線.若有多個坐標位置，則每個位置被轉換為參數空間中不同的直線.它們的交叉點在方程(1)的坐標系中指定了候選線的位置.在Hough空間的映射定義的累積器空間中，參數m和c轉化為坐標，從而獲取各自取值的有限區間，通過查找每個(m,c)坐標并增加其計數，最后高于閾值的坐標都可被選為候選線.

根據確定分組種子命中對角線方法來獲得種子粗對準的形式(見圖2)，結合參考基因組和測序序列組成的平面直角坐標系，單個種子的命中在此坐標系中用k-point(p,t)表示，其中p是種子在測序序列read上的位置，t是種子在參考基因組上的位置.當序列沒有包含錯誤且是由序列模擬器生成時，k-points集合在直角坐標系中是完全共線的.m是已知的，只需要確定截距參數c，即可找到精確的映射位置.同時，c對應于參考基因組上(已經離散)的坐標，因此可確定分組種子對角線位置(見圖2).對于每個k-point(p,t)，它的c參數值可由(3)式確定：

c=t-p

(3)

圖2 對候選種子進行區域選擇Fig.2 Region selection by clustering of candidate seeds on the reference

算法1描述了基于區域選擇和簡潔de Bruijn圖的短序列精確比對算法(Short-read alignment using region filtering and compact de Bruijn graph，簡記為SRA-RFCdBG)

算法1.SRA-RFCdBG

輸入：參考基因組G，測序序列read集合R；

輸出：包含R在G上的正確比對百分比的結果results；

Begin

1.解析參考基因組G建立k-mer文件；

2.依據k-mer文件生成overlap表格；

3.index←對參考基因組G建立k-mer索引；

4.fori=0 to |R|-1 do

5.regionsi←region_selection(Ri)；//對Ri使用算法2進行

區域選擇

6.end for

7.dBG←對regionsi生成簡潔de Bruijn圖；

8.results←測序序列read集合R與dBG的對準結果；

End

算法2.region_selection

輸入：參考基因組G，read序列r，k-mer長度k，基于G建立的索引k_index,設定的區域大小bin_size，最小閾值min_percent；

輸出：輸出篩選后的種子區域向量bins；

Begin

1.bins←利用k_index和G創建初始區域向量；

2.total_num_hits←0;

3.fori=1 to |r|-14 + 1 do //選取14-mers作為種子

4. fork_index_id= 1 to |k_index| do

5.index←k_index[k_index_id]；

6. ifindex≠NULL then

7. 根據(6-1-6)模式建立shape，獲得SubSeed1和

SubSeed2

8. 將SubSeed1和SubSeed2根據DC分別構建三種key

對應誤配，插入種子，刪除種子，進行位置查找，獲得

命中集Hitseti

9. 根據(4-1-4-1-4)模式建立shape，獲得SubSeed3、SubSeed4和SubSeed5

10. 將SubSeed3、SubSeed4和SubSeed5根據DC分別構建三種key對應誤配，插入種子，刪除種子，進行位置查找，獲得互補命中集reHitseti

11. 統計Hitseti和reHitseti中相同的候選位置作為候選區域，進行區域選擇，并統計命中數，記作hit_counts

12.total_num_hits←total_num_hits+hit_counts；

13. forj=1 tohit_countsdo //Hough變換

14.x←i；//序列r的坐標x

15.y←seed_position(hit[j])；//G的坐標y

16.l_local←y-x；//絕對坐標轉化為hit本地坐標

17.bin_pos←?l_local/bin_size」;

//bin_pos為種子區域bin所屬的位置

18.bin_value[bin_pos]←bin_value[bin_pos]+1；

//第bin_pos個區域的命中數

19. end for

20. end if

21. end for

22.end for

23.max_value←max_element(bin_value)//bin_value最大值

24.fori=0 tobin_value_size-1 do

25. ifbin_value[i] >min_percent*max_valuethen //保

留大于閾值的區域bin

26.bin.bin_id←i;

27.bin.bin_value←bin_value[i];//區域bini命中數

28.bins.push(bin);//滿足條件的區域bin增加到種子區域向量bins

29. end if

30.end for

31.依據bin_value將種子區域bins按降序排序;

End

2.3 算法分析

本文給出的序列比對方法同時考慮了“替換”和“indel(插刪)”錯誤類型，其中“替換”類型錯誤會導致區域選擇時滿足條件的c值的最大計數減少，而“indel”類型錯誤會影響種子命中的對角線位置并使得正確位置c值的最大計數更準確.此外，對于算法2中種子區域bin的每個種子命中，它增加了使用式(3)確定的參數c值所對應的bin_value，然后按照命中數降序對篩選后的種子區域向量bins排序.為將搜索限制為最可能的區域，僅選擇大于閾值的區域進行比對，以保持或提高正確比對百分比，且減少比對時間.

設read序列集R共有|R|=n條長度為l的序列.算法2在序列r中搜索索引index對應的參考基因組與種子的命中hit，所需時間為O((l-k+1)×|index|×|hit|),其中|index|表示index的大小，|hit|表示種子命中的規模.算法1中步4迭代n次；步5所需時間為O((l-k+1)×|index|×|hit|)；步8中，序列r在比對中搜索遇到的簡潔de Bruijn圖的分支節點個數最多為b，每個節點有4個輸出弧，最壞情形下要搜索4b個節點，搜索時間為|seeds|×O(4b+l-k-1+b)，其中|seeds|為區域選擇中的種子個數.算法1所需時間如下：O(n×(l-k+1)×|index|×|hit|)+O(n×|seeds|×(4b+l-k-1+b))=O(n((l-k+1)×|index|×|hit|+|seeds|×(4b+l-k-1+b))).

因此，算法SRA-RFCdBG的時間復雜度為O(n((l-k+1)×|index|×|hit|+|seeds|×(4b+l-k-1+b)))，其中l>k,b≤l-k+1.而算法BrownieAligner[18]比對時間與正確比對百分比均與測序序列的長度l和分支節點數b′有關，其時間復雜度為O(n×|seeds′|×(4b′+l-k-1+b′),其中|seeds′|為種子個數，b′為比對搜索中遇到的分支節點個數.算法SRA-RFCdBG采用區域選擇方法篩選候選位置來減少候選位置個數，使得種子數量和節點數有所減少，即|seeds|小于|seeds′|、b小于b′、4b遠遠小于4b′.因此，算法SRA-RFCdBG所需時間少于算法BrownieAligner.

本文采用簡潔de Bruijn圖[23]結構來完成比對.該圖結構是一種壓縮和索引標記樹，類似于XBW-transform[24]結構，是壓縮和索引字符串的BWT[25]的擴展.簡潔de Bruijn圖是de Bruijn圖的一種緊湊表示，其需要m(2+logσ)個二進制位的存儲空間，(2+logσ)的值不受k值的影響，且σ是字母表的大小(在蛋白質序列比對中，σ=4).而算法BrownieAligner的空間復雜度O(mk)，其使用的de Bruijn圖中的k取值為31.這表明本文算法SRA-RFCdBG所需的內存空間小于算法BrownieAligner.

3 實 驗

3.1 實驗環境與數據

實驗使用的計算機為8核Intel(R)Core i7-6700k CPU@4.00GHz處理器，內存容量為32GB，操作系統為64位Ubuntu18.04LTS,采用C++語言編程實現算法.采用模擬數據集與真實數據集進行實驗測試.對于3個參考基因組：大腸桿菌(E.coli str.K-12 substr.DH10B，E.coli)、人類第21條染色體(Humanchr-21，Hc21)和果蠅(Drosophila melanogaster，DM)，使用ART[26]模擬器在兩種不同的Illumina平臺HiSeq 2000(100bp)和HiSeq 2500(150bp)上，生成長度分別為100bp和150bp的模擬數據序列(reads)，見下一頁的表1.真實數據采用Illumina平臺上的數據集，參考基因組和測序序列從NCBI1網站獲得，參考基因組規模從2Mbp到116Mbp，見下一頁的表2.依據實驗測試結果將所允許的種子區域向量bins最小閾值min_percent設為75%.

對于真實序列與參考基因組比對，首先使用BWA[27]將序列reads與參考基因組比對生成bwa.sam文件，其次使用samtools[28]從bwa.sam文件中提取出通過雙端(paired-end)序列(將序列存儲為正反兩條序列形成的序列)比對方式正確映射到參考基因組,并存儲到mappedPairs.sam文件中.接著使用基于CIGAR字符串和MD標記的sam2pairwise2重建每個序列的雙端序列.最后使用Python腳本從雙端序列和初始的真實數據中提取出無錯誤的參考基因組，以評估真實數據集序列與參考基因組比對的正確比對百分比.實驗將本文算法SRA-RFCdBG與有其他代表性的同類代表性算法BGREAT[16]、deBGA[17]和BrownieAligner[18]進行性能比較.首先從參考基因組構建簡潔de Bruijn圖，然后將測序序列(reads)與生成的簡潔de Bruijn圖進行比對.我們使用GNU時間命令測量算法運行的(wall clock)時間.

3.2 模擬數據集上的實驗

表3給出了在模擬數據集S1～S6上，4種算法將序列與參考基因組進行映射比對獲得的正確對準百分比.

從表3可見，對于6種模擬數據集，與其他3種算法相比，本文算法SRA-RFCdBG均獲得最高的正確對準百分比.這是由于SRA-RFCdBG算法采用的區域選擇方法進行篩選候選位置時考慮了空格符號對比對的影響，這樣可以減少序列與參考基因組比對的候選位置個數.此外，Illumina平臺產生的序列存在的主要錯誤類型為“替換”，但也存在“插入和刪除(indel)”錯誤.SRA-RFCdBG算法通過同時對這兩種類型錯誤的處理，提高了種子命中率，從而獲得更高的正確對準百分比.從表4也看到：每種算法在人類基因組數據集(S3,S4)上的正確對準百分比都低于其他4個數據集(S1，S2，S5,S6)上的正確對準百分比.這是由于人類基因組(S3,S4)中堿基序列重復率較高，導致算法將S3、S4數據集中序列與參考基因組比對時，獲得正確對準的百分比較低.

圖3給出了在6種模擬數據集S1～S6上，4種算法比對一百萬條序列所用的運行時間.

從圖3可以看出，在S1～S5數據集上，與其他3種算法相比，本文算法SRA-RFCdBG所需的時間最少.之所以獲得這樣的結果，是因為算法SRA-RFCdBG采用區域選擇方法篩選搜索區域，減少了搜索空間，從而減少了比對的候選位置個數，減少比對計算量.對于S6數據集，算法SRA-RFCdBG所需的比對時間稍多于BGREAT所需的時間，但少于其他兩個算法所需的時間.這是由于BGREAT采用早期停止機制，對數據集S6上的序列比對時，當序列堿基不匹配的數量快達到設置的閾值時，算法停止探索de Bruijn圖的節點；而SRA-RFCdBG算法為了獲得更高的正確對準百分比，采取在篩選后的區域中計算每個分支的對準計分上限，直到無法改進當前找到的最佳解決方案才丟棄分支，其所需的計算量稍高于BGREAT算法的計算量，但低于deBGA和BrownieAligner算法的計算量.也就是說，對于重復率率高達50%的大數據集S6的序列比對，本文算法SRA-RFCdBG通過增強一定的計算量來獲得提升正確對準百分比.

圖3 算法在模擬數據集上比對一百萬條reads所用的時間Fig.3 Average runtime of tools to align 1M reads for the simulated datasets

表1 實驗使用的模擬數據集Table 1 Artificial datasets used for the evaluation of graph aligner tools

表2 實驗使用的真實數據集Table 2 Real datasets used for the evaluation of graph aligner tools

表3 模擬數據集上正確對準百分比Table 3 Percentage of correct aligned reads on simulateddatasets

3.3 真實數據集上的實驗

為評測算法在真實數據集上進行序列比對的效果，表4給出4種算法在8個真實數據集R1～R8上正確對準測序序列與參考基因組的百分比.表4的結果表明，本文算法SRA-RFCdBG在細菌基因組(R1-R5)和真核基因組(R7，R8)上正確對準序列與參考基因組的百分比均高于其他3種算法.這是因為算法SRA-RFCdBG采用基于空位種子的區域選擇方法對種子進行篩選，根據DC堿基構造的不同類型的查找關鍵字構建hash索引，索引查找兩個滿足條件的互補shape，從而完成近似匹配的敏感查找，這樣即使在空位中存在“indel”的情形，也可以找到對應匹配的種子；而其他3種算法未考慮該情形對比對結果的影響，未能找出空位中存在“indel”錯誤的匹配的種子，這導致了較低的正確比對百分比.對于人類基因組的真實數據集R6，算法SRA-RFCdBG的正確對準百分比略微比deBGA的正確對準百分比低0.10%，但高于其他2種算法BrownieAligner和BGREAT的正確對準百分比.整體而言，相較于其他3種算法，算法SRA-RFCdBG在8個數據集上獲得了更高的正確對準百分比.

表4 真實數據集上正確對準的百分比Table 4 Percentage of correct aligned reads on real datasets

為評估算法在真實數據集上的運行效率，圖4給出了4種算法在真實數據集上將一百萬條序列reads與參考基因組進行比對所需的運行時間.

圖4 在真實數據集上比對一百萬條reads所需的時間Fig.4 Average run time of aligning 1,000,000 readsfor the real datasets

從圖4中可以看出，對于R1、R4和R5這3個數據集，4種算法中，算法SRA-RFCdBG所用的時間最少.對于3個較大的數據集R6、R7和R8，算法SRA-RFCdBG的運行時間明顯比BrownieAligner低.這是因為算法SRA-RFCdBG采用基于空位種子策略的區域選擇方法進行篩選，對種子進行“預定義”，利用從shape布局中生成的關鍵字hash索引，快速地進行近似匹配的敏感查找，從而減少比對的候選位置的個數，尤其是對于重復率較高的基因組，可減少比對的候選位置的個數更多,有效減少了需搜索的空間，同時提高了正確比對的百分比.對于R6、R7和R8數據集，算法SRA-RFCdBG所需的時間比deBGA稍多一點.這是因為deBGA利用RdBG-index(Reference de Bruijn Graph-index)實現基于de Bruijn圖的種子擴展方法，可以更快地處理參考基因組中的序列重復[17].對于數據集R2、R3和R8，算法SRA-RFCdBG所需的時間比BGREAT稍多一些，這是因為BGREAT算法的實現采用最小完美hash索引，它用來索引以(k-1)-mer開始或結束的序列片段的每個重疊部分，減少了查詢時間.需要指出的是，雖然deBGA和BGREAT在某些數據集上完成比對所需的時間少一些，但它們獲得的正確比對百分比低于本文的算法SRA-RFCdBG.

綜合模擬數據集和真實數據集的序列比對實驗結果表明，除了個別情形之外，本文算法SRA-RFCdBG獲得的序列正確比對百分比均高于其他3種算法；在比對效率方面，SRA-RFCdBG優于算法BrownieAligner，與算法BGREAT和deBGA相比各有優勢.

4 結束語

對于海量的測序序列reads與存在許多重復子序列的參考基因組的比對問題，本文采用基于空位種子搜索策略的區域選擇方法來篩選候選位置，可以減少比對的候選位置，有效減少了搜索空間，在提高或保持正確比對百分比的同時減少了比對計算時間，尤其是對高重復率的序列進行比對，本文算法獲得了更高的正確對準百分比；此外，通過采用簡潔de Bruijn圖僅存儲由參考基因組生成的de Bruijn圖中m條邊的字符串和長度為σ的數組，減少了比對過程所需的內存容量，使得算法所需的存儲空間下降到m(2+logσ)個二進制位，比對所需的內存容量不受k-mer長度k的影響.

近年來，隨著更多的基因組可用，高通量測序技術的新應用需要將海量的測序序列與多個參考基因組對準[29].如何在既確保正確比對百分比，又能快速完成海量測序序列與多個參考基因組的比對工作，是一個值得繼續研究的問題[17].下一步工作將本文算法研究擴展到海量測序序列與多個參考基因組同時比對的情形.