一种对多序列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