2019, Vol. 40
Seed-and-extend启发式方法是近些年被广泛研究的简化大规模基因数据库搜索工作的主要方法, 高效的种子选择算法可以大大加速read mapping的处理过程.现有基于seed-and-extend的映射器集成了大量有效的子选择机制, 如Hobbel Mapper[1], 带有自适应种子过滤器的GEM Mapper[2], 支持K-mer的FashHASH[3], 最优种子选择方法OSS[4]等.这些映射器都采用了鸽笼原理的种子过滤技术, 通常把read划分成不重叠的多个短的片段，称为种子(seeds).种子可以用于索引基因序列来减少搜索空间大小和提高映射处理速度.将种子在基因序列中的匹配位置预先存储到种子数据库中, 从而可以通过数据库的查找来快速检索匹配位置.常用的存储方法包括Hash表和BWT FM-Index方法[5-6].构建高效映射器的关键问题就是如何从基因序列中选择具有最小频率的种子集合来保证快速地检索操作.然而, 从read中选择具有最小频率的种子集合是一项困难的工作, 原因是一个任意长度的种子可以从read的任意位置选择, 导致相应的搜索空间庞大并且随着种子数量的增长呈指数增长.
本文基于seed-and-extend的read mapping映射器提出一种基于频率的合并种子选择FMSS(frequency based merge seed selecting)算法, 该算法通过区分位置频率来确定最小频率种子, 只需要扫描read的每一个位置一遍, 然后通过合并具有高频率种子为低频率种子并调整种子分隔位置来减少种子集合的总体频率.最后与目前最好的最优种子选择算法(OSS)[4, 7]进行比较，验证了FMSS的准确性和有效性.
1 问题定义定义1 最优种子选择问题:给定长度为L的read R, 参考基因序列ref, 编辑距离门限τ.设freq(r)表示任意字符串r在参考基因组序列ref中出现的频率.设π∈Π, π={seedi|i=1, 2, …, τ+1}表示read R拆分的一组种子集合, 其中Π表示种子选择的空间.最优种子选择问题就是找到一组种子满足:
|
(1) |
式中, 一组种子π是readR的一个划分, 即顺序连接seedi就会构成readR，也就是将readR拆成τ+1段.将readR拆成τ+1段共有CL－1τ种方法，因此, readR的种子选择共有CL－1τ种方法.最优选择方法就是在这些方法中选择频率和最小的作为最优的种子.
目前, 常见的方法有Average种子选择方法[8].例如已知长度为15的read R, 编辑距离τ为3, 也就是将readR划分成4段, 那么每段seed的长度就是4, 4, 4, 3.然而, 当数据集和种子数量较大时, 这种方法的代价是十分高昂的, 原因是需要考虑所有子串中的全部候选种子的组合.
文献[4]提出了OSS种子选择方法, 它是一种动态规划算法, 是目前所有基于鸽笼原理的seed-and-extend种子选择策略中最优种子选择方法.最优种子选择意味着超矩阵F=∑(τ+1)freq(segR(i))上的所有K-mer的总匹配位置的数量是最小的.现在要将read R划分为k=τ+1个段，使用H(i, j)表示通过把R[0, j]划分成i段的最小频数，那么H[k, L]表示将R分成k段的最小频数.为了计算H[k, L], 通过以下动态规划公式计算矩阵Hk×L确认方法.
|
(2) |
|
(3) |
但是该方法的计算复杂度高达O(kL2), 运行的时间代价很高.
为了提高种子选择方法的性能, 提出一种新的算法FMSS, 该方法能够高效地选择接近最小频率的种子集合.FMSS算法分为两个步骤:首先, 算法通过计算R中每一个位置频率识别出频繁出现种子的位置, 通过合并相邻频繁位置的高频率种子来产生具有较低频率的新种子, 直到产生x个种子为止; 然后, 通过调整选择的种子之间的分隔位置来进一步减少选择种子的总频率.
1.1 种子位置频率通常从一个给定read中选择具有最小频率的种子集合的困难在于随着种子数量的增加需要在一个指数空间内搜索候选种子集合[8], 由于一个种子可以从read中的任意位置选取, 观察到从read中的不同位置选取的种子的频率会有很大差异.为了描述这种差异性, 定义read中任意位置的位置频率见定义2.
定义2 位置频率:给定一个长度为L的read R, 让pos表示R上的位置, 且满足0≤pos≤L-1.设(pos, m)表示通过从位置pos向左或向右扩展m个位置选取的有效子串(即pos－m≥0和pos+m+1≤|R|), k表示扩展子串的个数, r代表一个整数.求位置pos的频率就是基于位置pos进行扩展子串, 求出所有子串频率之和的平均数:
|
定义3 种子分段频率:给定一个长度为L的readR, 让pos表示R上的位置, 且满足0≤pos≤L-1.SEG(i)代表readR的分段, R={SEG(i)|i=1, 2, …, τ+1}, 求位置pos的segment频率意味着求位置pos所占的分段i在参考序列出现的频数.所以定义位置pos的分段频率为
|
例1 给定一个read R=CCAGTGCATA TACGACTT，选定参考基因序列为:CCAGTGCTAG CTGCCCTG CCAGTGCAATCCAGGCT CCAG TGCGCATGCTAGCT CCAGTGCGCTAGCTACA CCAGTGCAAATACGACTT CCAGTGCAAGCCGG GGT CCAGTGCTACGTACGAGT CCAGTGCTA TTCGATAT CCAGTGCAAGCTAGCATG.
要想计算位置5的频率就是在位置5周围进行扩展子串, 如图 1所示, 可以得到n1, …, n6子串, 那么位置5的频数计算如下:
|
图 1 扩展位置得到的子串 Fig.1 Substrings obtained by position extension |
|
若将read R先进行平均分段, 也就是分成"CCAGT", "GCATA", "TACG", "ACTT"4段.这个分段策略的匹配ref的频数为:(9, 1, 3, 1), 总频数为14.其中"CCAGT"明显是一个频数高的分段.如果利用最优分段技术, 就会分成"CCAGTGCAT", "ATAC", "GAC", "TT".这个分段策略的匹配频数为:(1, 1, 1, 2)，总频数为5, 大大降低了种子segment的频数.图 2a表示read R上18个点的位置频数, 平均分段策略和最优分段策略的频数.
|
图 2 种子位置频数和三种合并策略 Fig.2 Seed location frequency and three merge strategies (a)—平均和最优分段策略的频率；(b)—低频率与低频率；(c)—低频率与高频率；(d)—高频率与高频率 |
根据定义3, F(pos, AVE, r)给出了查询字符串上每一个点处的一个频率值, 这个值越高说明在该点周围的segment在ref上出现频数会越大, 反之该点周围的segment出现频数会越小[9].如图 2所示, 在一个read R上, 各个位置的PAF值差异很大.利用图 2, 可以将R分为N段, 总频数用式(4)计算.
|
(4) |
也就是F(pos, AVE, r)曲线与x轴围成的区域的面积.而分段后的分段总频数可以通过F(pos, SEG, r)曲线的积分面积来计算[10].
一个很好的性质是, 两段相邻的segments合并后的segment的频数会小于等于合并前的任意segment的频数.那么有3种合并策略:①低频segment和低频segement合并, 找到当前频数和最小的两个相邻segment(如位置17和18)并将它们合并, 接着在合并后的分段中继续找当前频数和最小的两个segment, 直到产生x个种子为止.②低频segment和高频segment合并, 找到当前频数差距最大的两个相邻segment合并(如位置7和8), 接着在合并后的分段中继续执行上面的操作, 直到产生x个种子为止.③高频segment和高频segment合并, 通过每次找到当前频数和最大的相邻的两段segments进行合并(如位置6和7), 逐步将所有的高频的segments都合并为尽量低频数的segments, 直到产生x个种子为止.种子初始化分段会将种子划分为每段只有一个字符的segment, 每个segment的位置频率如图 2a所示, 就初始化分段对segment进行3种合并策略进行讨论.
图 2b~图 2d分别表示对3种策略执行7次, 11次和最终合并的结果.如图 2b所示, 可以观察到策略①是不断地将低频逐步合并, 一直利用的是低频的分段, 合并最终原来高频的segment频数依旧很高, 并且高频率和低频率位置的PAF值差异很大, 浪费了低频资源, 显然这种方法是不可行的.如图 2c所示, 策略②这样的合并方法反而会使高频和低频的差值变大, 每一次都要将合并后的低频与相邻高频进行合并, 浪费低频资源, 时间代价也很大.而图 2d所示的策略③可以很快降低高频率的segment, 也没有浪费低频资源.相比图 2b和图 2c, 图 2d各个位置的PAF值差异也很小, 曲线与x轴围成的面积也是最小的.综上, 高频segment与高频segment合并是最好的策略.
MergeSeed算法实现了合并种子策略(算法1).为从read R中选择x个种子, 算法1分为3个步骤.首先, 计算readR中全部的位置频率; 然后, 算法根据位置频率找出具有最高频率的相邻种子; 最后, 通过合并具有最高频率的相邻种子为一个新的低频率种子直到获得read的x个种子终止.给定长度为L的read, 算法1至多迭代L－x次终止计算.
定理1 算法MergeSeed复杂性.从长度为R的read中确定x个种子, 算法1最坏情况下需要Ο((L/k－x)×L/k)步操作.算法首先计算R的全部位置频率, R中至多有L/k个位置, 接下来算法迭代Ο(L/k－x)次来确定x个种子, 在每一次迭代中, 算法检查Ο(L－2k)个种子对并进行合并种子操作.综上, 整个算法的复杂性为Ο((L/k－x)×L/k).
| 算法1 种子合并 |
| 输入:read R和种子数目x 输出:read R的x个种子集合S，频率和接近于最小频率 1:初始化种子位置S={|R|/k种子的长度相等} 2:For R中每一个种子位置p 3: 计算R的位置频率 4:end for 5:For i=1 to |R|/k 6： 从p位置选取频率最高的两个相邻种子si和sj 7:If freq(si)+freq(sj)>freq(sk) then 8: 合并si andsj称为种子sk并且把sk加入S集合 9:else 10:考虑下一组种子对 11:End If 12:End For 13:返回种子集合S |
1.3 调整种子分隔算法
合并相邻高频率种子之后, 可以获得一个相对较低频率的种子集合, 还可以通过调整已选择种子之间的分隔位置来进一步降低种子集合的总体频率.对种子之间的频率差值，通过调整位置p左移1个位置获得新的种子对为(s′i, s′j), 则
|
如果Δ(f)大于0, 可以通过调整位置p左移1个位置来把种子对(si, sj)优化成(si′, sj′).算法通过不断调整具有最大频率差值的相邻种子对, 可以进一步降低read中已选择的种子集合的总体频率.
ShiftDivider算法实现了种子移动位置策略(算法2).为进一步降低R中已选择的种子集合S的总体频率, 算法2首先计算S中相邻种子对的频率差值, 然后调整种子之间的分隔位置来减少种子对的频率和, 通过至多L次迭代, 将获得具有更低总体频率的种子集合.
定理2 算法ShiftDivider复杂性.为减少长度为L的read的x个种子集合的总体频率, 算法2最坏情况下需要Ο((x/2)×x)步操作.算法迭代检查R中具有最大频率差值的种子对并调整其分隔位置, 至多有x/2个种子对需要检查和x个位置需要调整.当不存在两个相邻种子频率差值大于0时, 算法终止并返回具有更低总体频率的种子集合.
| 算法2 ShiftDivider |
| 输入:x个种子的集合S 输出:x个种子的新集合S′，总频率小于S 1:初始化S′=S 2:计算S中所有相邻种子的频率差值 3:for最大的Δ(f)两个相邻的种子si和sj do 4:if Δ(f)>0 then 5: 将si和sj的分隔位置pij向左移动1-bp的位置 6: 移除种子si和sj并且将新种子si′和sj′加入到S′中 7: 考虑其他最大频率差值的种子 8:else 9: 返回种子集合S′ 10:end if 11:end for 12:返回种子的集合S′ |
2 实验结果
本文比较了提出的算法FMSS、采用平均长度种子选择策略AVE最优种子选择算法OSS的准确性和有效性.准确性采用给定read的选取种子集合的总体频率的平均值度量, 有效性采用选取种子集合的平均时间代价度量.全部算法采用C++语言实现, 运行在一台CPU为Intel (R) Core (TM) i7-2600 3.40 GHz, 主存大小为32 GB的PC机上, Ubuntu 16的操作系统.
实验数据来源于真实的公开DNA数据库和DBLP数据库, 随机选取4 000条read作为查询串, read的长度为108 bp.在DBLP上选取200条片段作为查询, 每条长度为78 bp.将每条read在平均长度种子选择策略Average, 最优种子选择算法OSS和本文提出的FMSS 3种策略下进行分段.图 3表示两个数据库下每条read在3种策略下分成的segment在4 000条read中出现的频数和.其中, “-a”表示AVE方法, “-b”表示OSS方法, “-c”表示本文提出的FMSS方法.可以发现, 对于不同种子数量, 在DNA数据库下的FMSS选择的平均种子频率总和接近于算法OSS, 但远远小于算法AVE.在DBLP数据库的FMSS选择的平均种子频率总和远远小于OSS和AVE.
|
图 3 实验结果 Fig.3 Experimental results |
同时, 3种策略的运行时间如表 1所示, 可以发现在DBLP数据库下FMSS消耗的运行时间基本与算法AVE接近, 远远低于OSS算法.在DNA数据库下FMSS消耗的运行时间高于算法AVE, 远远低于OSS算法.表明FMSS能够用更少的时间, 来选择接近最优频率的种子集合, 因此, 可用于处理较大规模的read mapping问题.
|
表 1 运行时间 Table 1 Running time |
1) 本文提出了一种有效的基于频率的种子选择算法FMSS, 解决了read mapping过程中低频数种子的选择问题.
2) 实验结果表明FMSS生成种子的频数与最优种子频数相差很小.同时FMSS算法的时间花费远远低于最优种子选择算法OSS.表明FMSS算法更适用于处理较大规模的read mapping问题.
| [1] |
Paolo F, Giovanni M, Veli M, et al.
Compressed representations of sequences and full-text indexes[J]. ACM Transactions on Algorithms, 2007, 3(2): 1–25.
|
| [2] |
Xin H Y, Nahar S, Zhu R, et al.
Optimal seed solver:optimizing seed selection in read mapping[J]. Bioinformatics, 2016, 32: 1632–1642.
DOI:10.1093/bioinformatics/btv670 |
| [3] |
Flicek P, Birney E.
Sense from sequence reads:methods for alignment and assembly[J]. Nature Methods, 2009, 6(1): 6–12.
|
| [4] |
Bin M, John T, Ming L.
Patternhunter:faster and more ensitive homology search[J]. Bioinformatics, 2002, 18: 440–445.
DOI:10.1093/bioinformatics/18.3.440 |
| [5] |
Xin H Y, Lee D, Hormozdiari F, et al.
Accelerating read mapping with FastHASH[J]. BMC Genomics, 2013, 14(1): 1–13.
DOI:10.1186/1471-2164-14-1 |
| [6] |
Santiago M, Michael S, Roderic G, et al.
The GEM mapper:fast, accurate and versatile alignment by filtration[J]. Nature Methods, 2012, 9(12): 1185–1187.
DOI:10.1038/nmeth.2221 |
| [7] |
Suzuki S, Kakuta M, Ishida T, et al.
Faster sequence homology searches by clustering subsequences[J]. Bioinformatics, 2015, 31: 1183–1190.
DOI:10.1093/bioinformatics/btu780 |
| [8] |
Hao W, Jeffrey X, Can L.
String similarity search:a Hash-based approach[J]. Transactions on Knowledge & Data Engineering, 2018, 99: 170–184.
|
| [9] |
Heng L, Richard D.
Fast and accurate read alignment with burrows-wheeler transform[J]. Bioinformatics, 2009, 25: 1754–1760.
DOI:10.1093/bioinformatics/btp324 |
| [10] |
Ahmadi A, Behm A, Honnalli N, et al.
Hobbes:optimized gram-based methods for efficient read alignment[J]. Nucleic Acids Research, 2012, 40(6): 41.
DOI:10.1093/nar/gkr1246 |