一種對多序列bwt索引構建進行并行加速的方法
【專利說明】
[0001]
技術領域:本發(fā)明涉及生物信息領域全基因組的組裝方法,尤其是在全基因組組 裝過程中大規(guī)模短序列集合(1億條以上序列)的Burrows-Wheeler變換(以下簡稱BWT) 索引構建的并行加速方法。
【背景技術】:
[0002] 全基因組組裝是生物信息學領域的核心問題,是基因組學其他相關研宄的基礎 與前提。一般生物的基因組包含數(shù)百萬乃至數(shù)十億個堿基,而當前的基因測序技術一次 只能測得包含數(shù)百個堿基的序列片段,根據(jù)測序所得到的短序列之間的重疊關系來將短 序列還原成原基因組的過程稱為基因組組裝。對于N條序列片段,直接計算它們兩兩之間 的重疊關系需要〇(N2)的時間復雜度,而對真核生物測序得到的序列片段的數(shù)量可高達數(shù) 億條,無法在有效時間內完成序列片段重疊關系的計算。研宄發(fā)現(xiàn),在已知序列片段集合的 BWT索引的如提下,序列片段間的重萱關系計算可在數(shù)小時內完成。序列集合的BWT索引定 義如下文所示。
[0003] 令2 = {Cl,c2,...,c。}為一個有限的字母表,滿足Cl<c2<…<c。,其中'〈' 表示字典序,〇表示字母表2中的字符個數(shù),1,2…〇為字母表中字母的順序號。令S= Sls2. ..Si. ..Sh為一個有限的字符串,其中SiG2 ;另外,定義字符串的末尾字符用'$' 表示,'$'在字典序上小于字母表2中的任一個字符。于是,S可以寫為長度為k的字符 串,S=SiS2. ? ?Sj^Sk,其中sk= $。我們用S[i,j] =s而+1. ? ?Sj表不S的從弟i個子付 到第j個字符構成的子串,其中1彡i彡j彡k。形如S[l,i]的子串稱作S的前綴,形如 S[j,k]的子串稱作S的后綴,其中1彡i,j彡k。稱s[j-l]為后綴S[j,k]的BWT字符,其 中1 <j彡k;令后綴S[l,k]的BWT字符為' $'。令R=以,S2,. . .,SJ表示字母表2上 的m條字符串,長度為='$'。為了區(qū)分不同的字符串,定義Sjk] <Sj[k], 對于1 <i<j<m。對序列集合R中序列的所有后綴按字典序排序,然后依次取各后綴的 BWT字符所組成的字符串就稱作序列集合R的BWT。
[0004] 由BWT的定義可以看出,構建BWT的主要步驟是對序列集合中所有序列的后綴按 字典序進行排序。然而,對于大規(guī)模短序列集合,直接對其包含的序列的所有后綴進行,排 序所需的主存大小高達TB數(shù)量級。以人類為例,人類基因組包含約30億個堿基,每個堿基 可用A、C、G、T中的一個字符表示。典型的深度為30X的測序將產(chǎn)生約10億條長約100個 堿基的序列,僅列舉這些序列的所有后綴就需要1. 25TB的空間,遠遠超過了現(xiàn)有計算設備 的內存大小。
[0005] 為此,研宄者提出了一種分塊排序而后再遞歸地兩兩合并的方法。以序列集合R =的BWT索引構建為例,假設計算設備的主存大小可滿足兩條序列的所有 后綴的排序,則將1?分為4塊1? 1,1?2,1?3,1?4,其中氏={521_ 1,心},1彡1彡4。首先依次對 RpR2,R3,R4中的所有后綴進行排序分別得到Sort^Sort2,Sort3,Sort4然后采用合并排序 算法對SortpSort2合并排序得到Sort12,再對Sort3,Sort4合并排序得到Sort34,之后對 SortiJPSort34合并排序得到Sort1234,即得到了R中所有后綴的字典序,最后按照從小到 大的字典順序依次取3〇竹1234各后綴的BWT字符,所組成的字符串即為序列集合R的BWT索 弓丨。這種分塊排序而后遞歸合并排序結果來構建BWT索引的方法解決了直接對所有序列的 后綴進行排序內存需求過大的問題,然而,由于這種方法是串行執(zhí)行的,時間效率差,無法 滿足大規(guī)模短序列集合BWT構建的時效性要求。
[0006] 為了加速BWT索引的構建過程,有研宄者以上述分塊-合并構建BWT索引的方法 為基礎,提出了一種并行加速方法。仍以字符串集合R= {3。$,...,$}的BWT索引構建 為例。首先配備包含4個節(jié)點PpP2,P3,P4的機群系統(tǒng),每個節(jié)點的主存大小可滿足兩條序 列的所有后綴的排序。然后將R分為4塊1?1,1?2,1?3,1? 4,其中氏={52卜1,52丄1彡1彡4。 同時在節(jié)點Pi上對Ri中的所有后綴進行排序得到Sorti,1彡i彡4 ;然后在節(jié)點Pi上對 SortpSort2合并排序得到Sort12,同時在節(jié)點P3上對Sort3,Sort4合并排序得到Sort34,之 后在節(jié)點PJtSort12和Sort34合并排序得到Sort1234,即得到了R中所有后綴的字典序,最 后按照從小到大的順序依次取3〇竹1234各后綴的BWT字符,所得到的字符串即為R的BWT索 弓丨??梢钥闯?,在初始階段,這種策略的并行度為分塊的總數(shù)目4,隨著合并的進行并行度每 次折半,到合并的最后一步并行度降為1,也就是完全串行執(zhí)行,整體的平均并行度低下。 經(jīng)實驗驗證加速比只有三分之一左右,仍舊無法滿足大規(guī)模序列BWT索引構建的時效性要 求。
[0007] 上述兩種方法都是直接對大規(guī)模序列集合進行分塊然后對各分塊進行排序,這種 分塊方法能解決直接排序內存需求過大的問題。然而由于各分塊的后綴之間沒有特定的大 小關系,在遞歸地合并分塊的過程中仍要對來自不同分塊的后綴的大小進行比較,極大地 降低了整體效率。
【發(fā)明內容】
:
[0008] 本發(fā)明要解決的技術問題是現(xiàn)有大規(guī)模序列集合BWT索引構建中采取對序列集 合進行分塊排序,然后遞歸地兩兩合并排序的方式,造成大規(guī)模序列集BWT索引構建速度 較慢,效率低下的問題。
[0009] 解決本發(fā)明技術問題所采用的技術方案是:首先遍歷序列集合R中各序列的所有 后綴,檢查各后綴的前1個字符,把具有相同前1個字符的后綴劃分到同一個內存分塊;之 后并行地在各個分塊內部獨立地對該分塊所包含的后綴進行字典序排序;然后將各排好序 的分塊按照分塊相應前1個字符字典序從小到大的順序拼接起來,得到序列集合R中所有 后綴的字典序;最后按字典序從小到大的順序依次取各后綴的BWT字符連接起來得到序列 集合R的BWT索引。
[0010] 具體技術方案如下:
[0011] 步驟1:根據(jù)序列規(guī)模和處理器內存大小確定對后綴分塊時所用分隔字符串的長 度1。對于序列集合R={Si,…,Si,…Sm},其中長度為k,l彡i彡m,字母表2包 含的字符個數(shù)為〇,處理器內存為M(byte)的情形,
[0012] 步驟2:構建含有o1個處理器(CPU)的機群系統(tǒng),依次編號為
[0013] 步驟3:在機群系統(tǒng)主存中開辟〇 1個動態(tài)內存分塊(以下簡稱桶),初始大小為 mk2/(4〇 4字節(jié),標號依次為1到〇 \
[0014] 步驟4:對序列集合R中包含的mXk條后綴進行分區(qū)。
[0015] 步驟 4.1:置i=l;
[0016] 步驟4. 2 :對第i條序列的后綴進行分區(qū)。
[0017] 步驟 4.2. 1 :置j= 1 ;
[0018] 步驟4. 2. 2 :檢查第i條序列的后綴Si[j,k]的前1個字符Si[j,j+1-l],對于長度 不足1的后綴,在其末尾添加字符Cl(如【背景技術】中所述,Cle= {(^,(^,...,(^:^序 列集合1?中的序列只包含〇1(:2...(3。這〇個字符,且按字典序(3 1<(32<~<(3。,1,2...〇 為字母表中字母的順序號)直至達到1長度。若Si[j,j+1-l] =Cilci2~Cil,其中 il,i2, . . .,il分別是Sjj,j+1-1]中包含的1個字符在字母表2 = {Cl,c2, . . .,c。}中的 順序號,則將后綴Si[j,k]放入編號為h的桶中,其中h= (il-l)X〇H+(i2_l)Xo1、-" + (il-l)Xoh+1的桶中;若桶中內存空間不足以存放新的后綴,則將其內存空間擴展mk2/ (16o1)字節(jié)。
[0019] 步驟 4.2.3 :置j=j+1 ;
[0020] 步驟4. 2. 4 :若j彡k,則轉步驟4. 2. 2,否則轉步驟4. 3.
[0021] 步驟 4. 3 :置i=i+1 ;
[0022] 步驟4. 4 :若i彡m,則轉步驟4. 2,否則轉步驟5.
[0023] 步驟5:在〇1個處理器上并行地對〇 1個桶內的后綴分別進行字典序排序,處理 器Pt對編號為t的桶內的后綴進行排序,1 <t< 〇 \
[0024] 步驟6:按照編號從1到〇1的順序將各個桶內已排好序的后綴拼接起來,得到R 中所有后綴的次序。序列集合R中包含m條序列,每條序列有k個后綴,序列集合R中一共 包含mXk條后綴,設這些后綴的字典序為SuffiXl&l