一种对多序列bwt索引构建进行并行加速的方法
一种对多序列bwt索引构建进行并行加速的方法
【专利说明】
[0001]
技术领域：本发明涉及生物信息领域全基因组的组装方法，尤其是在全基因组组 装过程中大规模短序列集合（1亿条以上序列）的Burrows-Wheeler变换（以下简称BWT) 索引构建的并行加速方法。
【背景技术】：
[0002] 全基因组组装是生物信息学领域的核心问题，是基因组学其他相关研宄的基础 与前提。一般生物的基因组包含数百万乃至数十亿个碱基，而当前的基因测序技术一次 只能测得包含数百个碱基的序列片段，根据测序所得到的短序列之间的重叠关系来将短 序列还原成原基因组的过程称为基因组组装。对于N条序列片段，直接计算它们两两之间 的重叠关系需要〇(N2)的时间复杂度，而对真核生物测序得到的序列片段的数量可高达数 亿条，无法在有效时间内完成序列片段重叠关系的计算。研宄发现，在已知序列片段集合的 BWT索引的如提下，序列片段间的重萱关系计算可在数小时内完成。序列集合的BWT索引定 义如下文所示。
[0003] 令2 = {Cl，c2，...，c。}为一个有限的字母表，满足Cl<c2<…<c。，其中'〈' 表示字典序，〇表示字母表2中的字符个数，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条字符串，长度为='$'。为了区分不同的字符串，定义Sjk] <Sj[k]， 对于1 <i<j<m。对序列集合R中序列的所有后缀按字典序排序，然后依次取各后缀的 BWT字符所组成的字符串就称作序列集合R的BWT。
[0004] 由BWT的定义可以看出，构建BWT的主要步骤是对序列集合中所有序列的后缀按 字典序进行排序。然而，对于大规模短序列集合，直接对其包含的序列的所有后缀进行，排 序所需的主存大小高达TB数量级。以人类为例，人类基因组包含约30亿个碱基，每个碱基 可用A、C、G、T中的一个字符表示。典型的深度为30X的测序将产生约10亿条长约100个 碱基的序列，仅列举这些序列的所有后缀就需要1. 25TB的空间，远远超过了现有计算设备 的内存大小。
[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索引的方法解决了直接对所有序列的 后缀进行排序内存需求过大的问题，然而，由于这种方法是串行执行的，时间效率差，无法 满足大规模短序列集合BWT构建的时效性要求。
[0006] 为了加速BWT索引的构建过程，有研宄者以上述分块-合并构建BWT索引的方法 为基础，提出了一种并行加速方法。仍以字符串集合R= {3。$，...，$}的BWT索引构建 为例。首先配备包含4个节点PpP2,P3,P4的机群系统，每个节点的主存大小可满足两条序 列的所有后缀的排序。然后将R分为4块1?1，1?2，1?3，1? 4，其中氏={52卜1，52丄1彡1彡4。 同时在节点Pi上对Ri中的所有后缀进行排序得到Sorti，1彡i彡4 ;然后在节点Pi上对 SortpSort2合并排序得到Sort12，同时在节点P3上对Sort3,Sort4合并排序得到Sort34，之 后在节点PJtSort12和Sort34合并排序得到Sort1234，即得到了R中所有后缀的字典序，最 后按照从小到大的顺序依次取3〇竹1234各后缀的BWT字符，所得到的字符串即为R的BWT索 弓丨。可以看出，在初始阶段，这种策略的并行度为分块的总数目4,随着合并的进行并行度每 次折半，到合并的最后一步并行度降为1，也就是完全串行执行，整体的平均并行度低下。 经实验验证加速比只有三分之一左右，仍旧无法满足大规模序列BWT索引构建的时效性要 求。
[0007] 上述两种方法都是直接对大规模序列集合进行分块然后对各分块进行排序，这种 分块方法能解决直接排序内存需求过大的问题。然而由于各分块的后缀之间没有特定的大 小关系，在递归地合并分块的过程中仍要对来自不同分块的后缀的大小进行比较，极大地 降低了整体效率。
【发明内容】
：
[0008] 本发明要解决的技术问题是现有大规模序列集合BWT索引构建中采取对序列集 合进行分块排序，然后递归地两两合并排序的方式，造成大规模序列集BWT索引构建速度 较慢，效率低下的问题。
[0009] 解决本发明技术问题所采用的技术方案是：首先遍历序列集合R中各序列的所有 后缀，检查各后缀的前1个字符，把具有相同前1个字符的后缀划分到同一个内存分块；之 后并行地在各个分块内部独立地对该分块所包含的后缀进行字典序排序；然后将各排好序 的分块按照分块相应前1个字符字典序从小到大的顺序拼接起来，得到序列集合R中所有 后缀的字典序；最后按字典序从小到大的顺序依次取各后缀的BWT字符连接起来得到序列 集合R的BWT索引。
[0010] 具体技术方案如下：
[0011] 步骤1:根据序列规模和处理器内存大小确定对后缀分块时所用分隔字符串的长 度1。对于序列集合R={Si，…，Si，…Sm}，其中长度为k，l彡i彡m，字母表2包 含的字符个数为〇，处理器内存为M(byte)的情形，
[0012] 步骤2:构建含有o1个处理器（CPU)的机群系统，依次编号为
[0013] 步骤3:在机群系统主存中开辟〇 1个动态内存分块（以下简称桶），初始大小为 mk2/(4〇 4字节，标号依次为1到〇 \
[0014] 步骤4:对序列集合R中包含的mXk条后缀进行分区。
[0015] 步骤 4.1:置i=l;
[0016] 步骤4. 2 :对第i条序列的后缀进行分区。
[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)字节。
[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 t;Suffix2<~<SuffixmXk。
[0025] 步骤7:依次取SuffiXl，Suffix2,. ..，SuffixmXk的BWT字符，连接起来得到R的 BWT索引。
[0026] 步骤8:将BWT索引输出。
[0027] 本发明的有益效果是，设计了一种BWT索引构建的并行加速方法，从而有效提高 多序列BWT索引的构建过程，减少全基因组组装所需时间约90%。而且该方法还可利用到 其他涉及大规模排序的应用中，易于移植和推广。
【附图说明】：
[0028] 图1是本发明总体流程图。
【具体实施方式】：
[0029] 下面结合附图1，以在每个节点64GB内存的机群上构建10亿条长度为100的DNA 序列的BWT索引为例（以下简称本例），对本发明进行进一步的详细说明。DNA序列的字母 表 2 = {A，C，G，T}，大小为 4,即 〇 = 4。
[0030] 如图1所示，本发明提出的新型BWT索引并行加速算法主要包括8个步骤。
[0031] 步骤1 :根据序列规模和处理器内存大小确定对后缀分块时所用分隔字符串的长 度 1，取
对于本例，m= 109,k= 100，〇 = 4，M= 64X23°，计算 得到 / =「3.1439^ = 4。
[0032] 步骤2 :计算得到〇 ^ 4 4= 256,我们来配备含有256个处理器（CPU)的机群系 统，分别编号为PuP2, . . .，P256。
[0033] 步骤3 :在机群系统内存中开辟〇 256个桶，标号依次为1到256,初始大小为 109X10〇V(4X256)字节，约为10GB，分别用来存放以AAAA到TTTT开头的后缀。
[0034] 步骤4 :扫描10亿条长度为100的DNA序列的所有后缀（一共有109X100 = 1011 条后缀），检查各后缀的前4个字符。
[0035] 步骤 4. 1 :置i= 1 ;
[0036] 步骤4. 2 :对第i条序列的后缀进行分区。
[0037] 步骤 4. 2. 1 :置j= 1 ;
[0038] 步骤4. 2. 2 :检查第i条序列的后缀Si[j，100]的前4个字符Si[j，j+3]进行分区。 对于DNA，字母表 2 = {A，C，G，T}，所以Cl=，A'，c2=，C'，c3=，G'，c4= 'T'。字符'A' 的顺序号为1，字符'C'的顺序号为2,字符'G'的顺序号为3,字符'T'的顺序号为4。对 于长度不足4的后缀，在其末尾添加字符'A'直至长度达到4。以AAAA开头的后缀放到h =(1-1)X43+(l-l)X42+(l-l)Xf+d-l)X4°+l= 1 号桶，以AAAC开头的后缀放到h= (1-1)X43+(l-l)X42+(l-l)X4i+(2-l)X4°+l= 2 号桶，依次类推，以TTTT开头的后缀 放到h= (4-1)X43+(4-l)X42+(4-l)X4i+(4-l)X4°+l= 256 号桶。若桶中内存空间不足 以存放新的后缀，则将其内存空间扩展l〇9X10〇V(16X256)字节，约2. 5GB。
[0039] 步骤 4. 2. 3 :置j=j+1 ;
[0040] 步骤4. 2. 4 :若j彡100,则转步骤4. 2. 2,否则转步骤4. 3.
[0041] 步骤 4. 3 :置i=i+1 ;
[0042] 步骤4. 4 :若i彡109,则转步骤4. 2,否则转步骤5.
[0043] 步骤5 :在256个处理器并行地对256个桶中的后缀进行排序，具体方式为第t号 处理器pt对第t号桶中的后缀进行排序，1 <t< 256。
[0044] 步骤6 :按从1号桶到256号桶的顺序将排序的结果拼接起来，得到所有后缀的次 序。对于1〇9条长100的序列集合，m= 109，k= 100,其后缀一共有109X100 = 1011条，设 其字典序为SuffiXl<Suffix2<~<Suffix1(111。
[0045] 步骤 7 :依次取SuffiXl，Suffix2,. . .，Suffix1(ll^BWT字符，连接起来组成 10 亿 条
[0046] DNA序列的BWT索引。
[0047] 步骤8 :将BWT索引输出。
【主权项】
1. 一种对多序列BWT索引构建进行并行加速方法，其特征在于包括以下步骤： 步骤1:根据序列规模和处理器内存大小确定对后缀分块时所用分隔字符串的长度1 ; 令2 = {〇1，〇2，...，〇。}为一个有限的字母表，满足〇1〈〇 2〈~〈〇。，其中'〈'表示字典序，〇表 示字母表Σ中的字符个数；令S = S1S2-Sn1Sk为一个有限的字符串，其中S k= 为字符串结束标志，SiE Σ，1彡i〈k ;S[i，j] = s iSi+1. . . Sj表示S的从第i个字符到第j个 字符构成的子串，其中1彡i彡j彡k ;形如S[j，k]的子串称作S的后缀，其中1彡j彡k ; 称s[j-l]为后缀S[j，k]的BWT字符，其中l〈j彡k ;令后缀S[l，k]的BWT字符为' $' ；令 R = (S1, S2,. . .，SJ表示字母表Σ上的m条字符串，长度为k且S i [k] = ' $' ；对于处 理器内存大小为M字节的情形，I = [(log。（mkV(2M)))]; 步骤2:构建含有〇 1个处理器的机群系统，依次编号为A，A，···，Pfji ; 步骤3:在机群系统主存中开辟σ1个动态桶，初始大小为mkVG。1)字节，标号依次 为1到σ S桶为内存分块； 步骤4:对序列集合R中包含的m*k条后缀进行分区， 步骤4. 1 :置i = 1 ; 步骤4. 2 :对第i条序列的后缀进行分区， 步骤 4. 2. 1 :置 j = 1 ; 步骤4. 2. 2 :检查第i条序列的后缀Si [ j，k]的前1个字符Si [ j，j+1-l]，对于长度不足 1的后缀，在其末尾添加字符(^直至达到1长度，c i e Σ ;若S i [j，j+1-1] = cnci2··· Cil，则 将后缀SiU, k]放入编号为h的桶中，其中h = (i「l)* σ H+h-l)* σ K+ -" + (ifl)* σ 14+1，11，12，一，11分别是31[」，」+1-1]中包含的1个字符在字母表2 = {〇1，(：2，...'。}中 的顺序号，若桶中内存空间不足以存放新的后缀，则将其内存空间扩展mkVae。1)字节； 步骤 4. 2. 3 :置 j = j+Ι ; 步骤4. 2. 4 :若j彡k，则转步骤4. 2. 2,否则转步骤4. 3 ; 步骤 4. 3 :置 i = i+Ι ; 步骤4. 4 :若i彡m，则转步骤4. 2,否则转步骤5 ; 步骤5:在ο1个处理器上并行地对〇1个桶内的后缀分别进行字典序排序，处理器pt对编号为t的桶内的后缀进行排序，K t < σ S 步骤6:按照编号从1到σ1的顺序将各个桶内已排好序的后缀拼接起来，得到R中所 有后缀的次序（字典序），序列集合R中包含m条序列，每条序列有k个后缀，序列集合R中 一共包含m*k条后缀，设这些后缀的字典序为SufTiXl〈SufTix2〈…〈Suffix- 步骤7:依次取Suff iXl，Suff ix2，…，Suff丨1_的BWT字符，连接起来得到R的BWT索 引； 步骤8:将BWT索引输出。2. 根据权利要求1所述的一种对多序列BWT索引构建进行并行加速方法，其特征在于 步骤2中所述的ο1个处理器指〇1ACPl
【专利摘要】本发明公开了一种多序列BWT索引构建的并行加速方法，目的是解决现有大规模序列集合BWT索引构建中对序列集合进行分块排序，然后通过两两合并再排序的方式不断递归合并排序的过程，造成BWT索引构建速度较慢，效率低下的问题。采用的技术方案是首先遍历序列集合R中各序列的所有后缀，检查各后缀的前l个字符,把具有相同前l个字符的后缀划分到同一个内存分块；之后并行地在各个分块内部独立地对该分块所包含的后缀进行排序；然后将各排好序的分块拼接起来，得到序列集合R中所有后缀的次序；最后按字典序从小到大的顺序依次取各后缀的BWT字符连接起来得到序列集合R的BWT索引。有益效果是有效提高多序列BWT索引的构建过程，减少全基因组组装所需时间约90％。
【IPC分类】G06F19/22
【公开号】CN104899476
【申请号】CN201510328718
【发明人】彭绍亮, 朱小谦, 王恒, 卢宇彤, 杨灿群, 吴诚堃, 崔英博, 刘欣, 王海强, 程乾, 夏徐伟
【申请人】中国人民解放军国防科学技术大学
【公开日】2015年9月9日
【申请日】2015年6月15日
【专利说明】
[0001]
技术领域：本发明涉及生物信息领域全基因组的组装方法，尤其是在全基因组组 装过程中大规模短序列集合（1亿条以上序列）的Burrows-Wheeler变换（以下简称BWT) 索引构建的并行加速方法。
【背景技术】：
[0002] 全基因组组装是生物信息学领域的核心问题，是基因组学其他相关研宄的基础 与前提。一般生物的基因组包含数百万乃至数十亿个碱基，而当前的基因测序技术一次 只能测得包含数百个碱基的序列片段，根据测序所得到的短序列之间的重叠关系来将短 序列还原成原基因组的过程称为基因组组装。对于N条序列片段，直接计算它们两两之间 的重叠关系需要〇(N2)的时间复杂度，而对真核生物测序得到的序列片段的数量可高达数 亿条，无法在有效时间内完成序列片段重叠关系的计算。研宄发现，在已知序列片段集合的 BWT索引的如提下，序列片段间的重萱关系计算可在数小时内完成。序列集合的BWT索引定 义如下文所示。
[0003] 令2 = {Cl，c2，...，c。}为一个有限的字母表，满足Cl<c2<…<c。，其中'〈' 表示字典序，〇表示字母表2中的字符个数，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条字符串，长度为='$'。为了区分不同的字符串，定义Sjk] <Sj[k]， 对于1 <i<j<m。对序列集合R中序列的所有后缀按字典序排序，然后依次取各后缀的 BWT字符所组成的字符串就称作序列集合R的BWT。
[0004] 由BWT的定义可以看出，构建BWT的主要步骤是对序列集合中所有序列的后缀按 字典序进行排序。然而，对于大规模短序列集合，直接对其包含的序列的所有后缀进行，排 序所需的主存大小高达TB数量级。以人类为例，人类基因组包含约30亿个碱基，每个碱基 可用A、C、G、T中的一个字符表示。典型的深度为30X的测序将产生约10亿条长约100个 碱基的序列，仅列举这些序列的所有后缀就需要1. 25TB的空间，远远超过了现有计算设备 的内存大小。
[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索引的方法解决了直接对所有序列的 后缀进行排序内存需求过大的问题，然而，由于这种方法是串行执行的，时间效率差，无法 满足大规模短序列集合BWT构建的时效性要求。
[0006] 为了加速BWT索引的构建过程，有研宄者以上述分块-合并构建BWT索引的方法 为基础，提出了一种并行加速方法。仍以字符串集合R= {3。$，...，$}的BWT索引构建 为例。首先配备包含4个节点PpP2,P3,P4的机群系统，每个节点的主存大小可满足两条序 列的所有后缀的排序。然后将R分为4块1?1，1?2，1?3，1? 4，其中氏={52卜1，52丄1彡1彡4。 同时在节点Pi上对Ri中的所有后缀进行排序得到Sorti，1彡i彡4 ;然后在节点Pi上对 SortpSort2合并排序得到Sort12，同时在节点P3上对Sort3,Sort4合并排序得到Sort34，之 后在节点PJtSort12和Sort34合并排序得到Sort1234，即得到了R中所有后缀的字典序，最 后按照从小到大的顺序依次取3〇竹1234各后缀的BWT字符，所得到的字符串即为R的BWT索 弓丨。可以看出，在初始阶段，这种策略的并行度为分块的总数目4,随着合并的进行并行度每 次折半，到合并的最后一步并行度降为1，也就是完全串行执行，整体的平均并行度低下。 经实验验证加速比只有三分之一左右，仍旧无法满足大规模序列BWT索引构建的时效性要 求。
[0007] 上述两种方法都是直接对大规模序列集合进行分块然后对各分块进行排序，这种 分块方法能解决直接排序内存需求过大的问题。然而由于各分块的后缀之间没有特定的大 小关系，在递归地合并分块的过程中仍要对来自不同分块的后缀的大小进行比较，极大地 降低了整体效率。
【发明内容】
：
[0008] 本发明要解决的技术问题是现有大规模序列集合BWT索引构建中采取对序列集 合进行分块排序，然后递归地两两合并排序的方式，造成大规模序列集BWT索引构建速度 较慢，效率低下的问题。
[0009] 解决本发明技术问题所采用的技术方案是：首先遍历序列集合R中各序列的所有 后缀，检查各后缀的前1个字符，把具有相同前1个字符的后缀划分到同一个内存分块；之 后并行地在各个分块内部独立地对该分块所包含的后缀进行字典序排序；然后将各排好序 的分块按照分块相应前1个字符字典序从小到大的顺序拼接起来，得到序列集合R中所有 后缀的字典序；最后按字典序从小到大的顺序依次取各后缀的BWT字符连接起来得到序列 集合R的BWT索引。
[0010] 具体技术方案如下：
[0011] 步骤1:根据序列规模和处理器内存大小确定对后缀分块时所用分隔字符串的长 度1。对于序列集合R={Si，…，Si，…Sm}，其中长度为k，l彡i彡m，字母表2包 含的字符个数为〇，处理器内存为M(byte)的情形，
[0012] 步骤2:构建含有o1个处理器（CPU)的机群系统，依次编号为
[0013] 步骤3:在机群系统主存中开辟〇 1个动态内存分块（以下简称桶），初始大小为 mk2/(4〇 4字节，标号依次为1到〇 \
[0014] 步骤4:对序列集合R中包含的mXk条后缀进行分区。
[0015] 步骤 4.1:置i=l;
[0016] 步骤4. 2 :对第i条序列的后缀进行分区。
[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)字节。
[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 t;Suffix2<~<SuffixmXk。
[0025] 步骤7:依次取SuffiXl，Suffix2,. ..，SuffixmXk的BWT字符，连接起来得到R的 BWT索引。
[0026] 步骤8:将BWT索引输出。
[0027] 本发明的有益效果是，设计了一种BWT索引构建的并行加速方法，从而有效提高 多序列BWT索引的构建过程，减少全基因组组装所需时间约90%。而且该方法还可利用到 其他涉及大规模排序的应用中，易于移植和推广。
【附图说明】：
[0028] 图1是本发明总体流程图。
【具体实施方式】：
[0029] 下面结合附图1，以在每个节点64GB内存的机群上构建10亿条长度为100的DNA 序列的BWT索引为例（以下简称本例），对本发明进行进一步的详细说明。DNA序列的字母 表 2 = {A，C，G，T}，大小为 4,即 〇 = 4。
[0030] 如图1所示，本发明提出的新型BWT索引并行加速算法主要包括8个步骤。
[0031] 步骤1 :根据序列规模和处理器内存大小确定对后缀分块时所用分隔字符串的长 度 1，取
对于本例，m= 109,k= 100，〇 = 4，M= 64X23°，计算 得到 / =「3.1439^ = 4。
[0032] 步骤2 :计算得到〇 ^ 4 4= 256,我们来配备含有256个处理器（CPU)的机群系 统，分别编号为PuP2, . . .，P256。
[0033] 步骤3 :在机群系统内存中开辟〇 256个桶，标号依次为1到256,初始大小为 109X10〇V(4X256)字节，约为10GB，分别用来存放以AAAA到TTTT开头的后缀。
[0034] 步骤4 :扫描10亿条长度为100的DNA序列的所有后缀（一共有109X100 = 1011 条后缀），检查各后缀的前4个字符。
[0035] 步骤 4. 1 :置i= 1 ;
[0036] 步骤4. 2 :对第i条序列的后缀进行分区。
[0037] 步骤 4. 2. 1 :置j= 1 ;
[0038] 步骤4. 2. 2 :检查第i条序列的后缀Si[j，100]的前4个字符Si[j，j+3]进行分区。 对于DNA，字母表 2 = {A，C，G，T}，所以Cl=，A'，c2=，C'，c3=，G'，c4= 'T'。字符'A' 的顺序号为1，字符'C'的顺序号为2,字符'G'的顺序号为3,字符'T'的顺序号为4。对 于长度不足4的后缀，在其末尾添加字符'A'直至长度达到4。以AAAA开头的后缀放到h =(1-1)X43+(l-l)X42+(l-l)Xf+d-l)X4°+l= 1 号桶，以AAAC开头的后缀放到h= (1-1)X43+(l-l)X42+(l-l)X4i+(2-l)X4°+l= 2 号桶，依次类推，以TTTT开头的后缀 放到h= (4-1)X43+(4-l)X42+(4-l)X4i+(4-l)X4°+l= 256 号桶。若桶中内存空间不足 以存放新的后缀，则将其内存空间扩展l〇9X10〇V(16X256)字节，约2. 5GB。
[0039] 步骤 4. 2. 3 :置j=j+1 ;
[0040] 步骤4. 2. 4 :若j彡100,则转步骤4. 2. 2,否则转步骤4. 3.
[0041] 步骤 4. 3 :置i=i+1 ;
[0042] 步骤4. 4 :若i彡109,则转步骤4. 2,否则转步骤5.
[0043] 步骤5 :在256个处理器并行地对256个桶中的后缀进行排序，具体方式为第t号 处理器pt对第t号桶中的后缀进行排序，1 <t< 256。
[0044] 步骤6 :按从1号桶到256号桶的顺序将排序的结果拼接起来，得到所有后缀的次 序。对于1〇9条长100的序列集合，m= 109，k= 100,其后缀一共有109X100 = 1011条，设 其字典序为SuffiXl<Suffix2<~<Suffix1(111。
[0045] 步骤 7 :依次取SuffiXl，Suffix2,. . .，Suffix1(ll^BWT字符，连接起来组成 10 亿 条
[0046] DNA序列的BWT索引。
[0047] 步骤8 :将BWT索引输出。
【主权项】
1. 一种对多序列BWT索引构建进行并行加速方法，其特征在于包括以下步骤： 步骤1:根据序列规模和处理器内存大小确定对后缀分块时所用分隔字符串的长度1 ; 令2 = {〇1，〇2，...，〇。}为一个有限的字母表，满足〇1〈〇 2〈~〈〇。，其中'〈'表示字典序，〇表 示字母表Σ中的字符个数；令S = S1S2-Sn1Sk为一个有限的字符串，其中S k= 为字符串结束标志，SiE Σ，1彡i〈k ;S[i，j] = s iSi+1. . . Sj表示S的从第i个字符到第j个 字符构成的子串，其中1彡i彡j彡k ;形如S[j，k]的子串称作S的后缀，其中1彡j彡k ; 称s[j-l]为后缀S[j，k]的BWT字符，其中l〈j彡k ;令后缀S[l，k]的BWT字符为' $' ；令 R = (S1, S2,. . .，SJ表示字母表Σ上的m条字符串，长度为k且S i [k] = ' $' ；对于处 理器内存大小为M字节的情形，I = [(log。（mkV(2M)))]; 步骤2:构建含有〇 1个处理器的机群系统，依次编号为A，A，···，Pfji ; 步骤3:在机群系统主存中开辟σ1个动态桶，初始大小为mkVG。1)字节，标号依次 为1到σ S桶为内存分块； 步骤4:对序列集合R中包含的m*k条后缀进行分区， 步骤4. 1 :置i = 1 ; 步骤4. 2 :对第i条序列的后缀进行分区， 步骤 4. 2. 1 :置 j = 1 ; 步骤4. 2. 2 :检查第i条序列的后缀Si [ j，k]的前1个字符Si [ j，j+1-l]，对于长度不足 1的后缀，在其末尾添加字符(^直至达到1长度，c i e Σ ;若S i [j，j+1-1] = cnci2··· Cil，则 将后缀SiU, k]放入编号为h的桶中，其中h = (i「l)* σ H+h-l)* σ K+ -" + (ifl)* σ 14+1，11，12，一，11分别是31[」，」+1-1]中包含的1个字符在字母表2 = {〇1，(：2，...'。}中 的顺序号，若桶中内存空间不足以存放新的后缀，则将其内存空间扩展mkVae。1)字节； 步骤 4. 2. 3 :置 j = j+Ι ; 步骤4. 2. 4 :若j彡k，则转步骤4. 2. 2,否则转步骤4. 3 ; 步骤 4. 3 :置 i = i+Ι ; 步骤4. 4 :若i彡m，则转步骤4. 2,否则转步骤5 ; 步骤5:在ο1个处理器上并行地对〇1个桶内的后缀分别进行字典序排序，处理器pt对编号为t的桶内的后缀进行排序，K t < σ S 步骤6:按照编号从1到σ1的顺序将各个桶内已排好序的后缀拼接起来，得到R中所 有后缀的次序（字典序），序列集合R中包含m条序列，每条序列有k个后缀，序列集合R中 一共包含m*k条后缀，设这些后缀的字典序为SufTiXl〈SufTix2〈…〈Suffix- 步骤7:依次取Suff iXl，Suff ix2，…，Suff丨1_的BWT字符，连接起来得到R的BWT索 引； 步骤8:将BWT索引输出。2. 根据权利要求1所述的一种对多序列BWT索引构建进行并行加速方法，其特征在于 步骤2中所述的ο1个处理器指〇1ACPl
【专利摘要】本发明公开了一种多序列BWT索引构建的并行加速方法，目的是解决现有大规模序列集合BWT索引构建中对序列集合进行分块排序，然后通过两两合并再排序的方式不断递归合并排序的过程，造成BWT索引构建速度较慢，效率低下的问题。采用的技术方案是首先遍历序列集合R中各序列的所有后缀，检查各后缀的前l个字符,把具有相同前l个字符的后缀划分到同一个内存分块；之后并行地在各个分块内部独立地对该分块所包含的后缀进行排序；然后将各排好序的分块拼接起来，得到序列集合R中所有后缀的次序；最后按字典序从小到大的顺序依次取各后缀的BWT字符连接起来得到序列集合R的BWT索引。有益效果是有效提高多序列BWT索引的构建过程，减少全基因组组装所需时间约90％。
【IPC分类】G06F19/22
【公开号】CN104899476
【申请号】CN201510328718
【发明人】彭绍亮, 朱小谦, 王恒, 卢宇彤, 杨灿群, 吴诚堃, 崔英博, 刘欣, 王海强, 程乾, 夏徐伟
【申请人】中国人民解放军国防科学技术大学
【公开日】2015年9月9日
【申请日】2015年6月15日
最新回复(0)