1、 1 / 8 S-W 算法在 NGS 应用下的测试与分析 摘要： 生物序列比对是生物信息学中 最常见 的问题之一，它 承载着 正确还原和判断生物序列的重要工作 。 NGS 在 吞吐 率 提高的同时， 比对的错误率也有明显上升， 如何 提高 比对命中率成了 工程师们不得不面对的一个问题。本文通过实验样品的特异性差异， 针对 性的 对不同样本 状况 调整 碱基匹配的 优先级策略， 整体上 提高了比对的命中 率 。 关 键词： S-W 算法、 NGS 应用 、测试 、 分析 生物信息学中，生物 信息 需要借助计算机编码 来表达 遗传信息 ， 基于 已有的技术手段 ， NGS(New Generati
2、on Sequencing-新一代 基因测序 )已经足够成熟 ， 由 NGS 测序仪 上 能 输出 大量 长度较短的碱基序列， 它们 以 固定的 字符编码。 为了 还原 原始 遗传 信息 ， 大量的短碱基序列需要和 权威性 较高的参考基因组序列进行比对。 由于 参考基因序列长度 远远 长于测序仪中产生的碱基序列， 比对 过程中，短 序列需要不断 地 在参考序列的不同位置 “游荡 ”， 尽量 找到属于自己 正确的 位置 。 但是由于人体各异， 基因组 状况各不相同， 不同算法策略下的比对结果也不尽相同。此时， 如何评估 序列比对结果的科学性与准确 性 成 为当务之急。本文 通过 不同类型 样本
3、间不同表达倾向的差异， 测试了 S-W 算法在 不同 策略条件下对序列比对结果的影响，从而给出了序列比对流程 中 的一些优化性 建议和 参考 。 1 序列 比对 算法 1.1 序列 相似 性 比较 在 生物 信息 学中， 对 各种生物大分子 序列进行分析 是 一件基本工作。 从序列 的判断 、测定、 拼接 、基因的表达分析， 到 RNA(Ribonucleic Acid)和蛋白质的结构 功能预测 、 物种亲缘 树的构建都需要进行生物分子序列的相似性比较。 在遗传物质 长期演化过程中， 原本 相同的 DNA(Deoxyribonucleic acid)序列 由于 其中 一条序列缺失或 增加了 几
4、个片段， 或者 某段子序列发生了位置的变化等， 从而 导致他们 发生变化 ，不 能 保证 其 精确匹配， 但它 们 有一定的 相似度 。 我们应该如何判定序列之间的这种相似性呢？ 对于这种情况， 生物学家 提出了一种 评定 序列相似 性 的方法 记分 函数法。 为了后文 阐述方便， 先看如下几个定义： 定义 1:如果 S 是一个序列， 那么 |S|表示 S 中的字符 长度 ， Si表示 序列的第 i 个字符。 如果 序列 S 和序列 T 相同， 必须 满足如下条件： (1).|S|=|T| (2).Si=Ti,(00) if(Mi,j=Mi-1,j-1+(,) i-,j-; else if(M
5、i,j=Mi-1,j+(,) i-;insert(-,T,j); else if(Mi,j=Mi,j-1+(,) j-;insert(-,T,i); 其中 insert(a,b,c)表示 在 序列 b 的 第 c 个 位置 插入 a。 我们假设 在 S-W 算法 中 |S|=m,|T|=n,初始化 得分数组时需要计算二维数组中每一个元素的 得分 ，其 时间复杂 度 为 O(m*n)。 由 回溯算法得知其时间复杂度为 O(max(x,y)。由此 得知 ， S-W 算法 的 时间 复杂度为 O(m*n) ， 算法 的主要工 作 量 都 集 中 在 了计 算 回 溯 数 组 的 得分上。当Mi-1,
6、j-1,Mi-1,j,Mi,j-1都 已知时，我们才 认为 Mi,j是 可计算的。 2 影响 S-W 算法执行效果的因素探究 2.1 NGS 业务 背景 NGS(New Generation Sequencing)是相对于第一代基因测序（ Sigar 测序）来说有较 大 优势的测序4 / 8 方法， 优点 有：成本低， 高通量 ， 周期短 等特点 ； 其缺点 就是后期需要经过大量的数据处理增 大 了数据还原不一致的风险 ， 因为 基因序列的 碱基字符 集 中 字符数量少， 只有 A,G,C,T 四种 ，这时就需要 更加 优秀的序列比对算法实现更低错误 率 的序列判断。 虽然 NGS 有高通量、
7、周期短等特点，但是 它 与 第一代 测序技术相比也有一些不可回避的问题。其中包括 很重要的一点就是 NGS 可能 得出 与供体 DNA 不 完全相同的序列 ， 但是 好的策略还是可以得到 高还原 度的比对结果的。 2.2 NGS 测序 原理 生物体 组织中， 具有 遗传性质的物质为脱氧核糖核酸（也就是 DNA）， 它 存在于生物体的各个组织器官中， 甚至 是人的唾液里。如果想得到这些遗传 信息 的 字符 序列表示， 需要 经过一下几个 大致 步骤： 首先 ， 要将 DNA 样品 打断 称 300-800bp 的 短 片断； 然后 ， 对短 片段加接头 并 结合 到 固相 表面使之相互独立的扩增
8、 ， 每次 只复制一个碱基并检测信号； 最后 ， 高分辨率 成像， 将 化学信号或者光信号转化为信号编码交由计算机存储。 通过 以上步骤我们就得到了大量短序列的 字符串集合， 但是这些集合都是 DNA 序列中的其中一小段，相对于 整条 DNA 片段来说犹如九牛一毛， 不能 得出任何 信息 ， 那么 如何才能得出可以 反映 整条序列信息的整体片段呢？ 这时 我们就需要根据现有的公认参考序列和 测序仪 中产生的大量经过扩增的序列来还 原 原始序列，那么还原的过程就是考验算法策略 优劣 的时候 了 。 在还原 原始 序列的过程中， 大量的短序列需要找到他们 在 原始序列中 对应的正确位置。 整个 寻
9、找过程 是一个 不断试探的过程， 每当 片段与参考序列开始新一轮的匹配的时候，当前位置就是我们怀疑的位置对象。 比对 完成后会得出 一个 可以量化的结果，之后我们对这个结果做一个 评估 ， 判断它到底是不是当前位置的片段， 如果 量化后的结果高于我们设置的阀值， 我们 就认为 他 可信， 是 当前位置的片段； 如果低于 阀值， 我们 会认为不可信， 此轮 比对到此为止， 准备 下一轮比对。 2.3 S-W 算法中可控因子对 NGS 业务的 影响 判断 一个算法的好坏有几个很重要的维度， 首先 ， 就是 看算法是否满足 算法 本身的特性： 有穷性， 确定性 ， 可行性 ， 有 输入， 有 输出；
10、其次， 就是 看算法在应用环境下对业务需求的契合程度。 NGS 的应用场景相对于一般的序列比对有其特殊性， 这种 特殊性表现为 字符集合的大小 为 4，5 / 8 元素 个数极为有限， 这种情况 下， 序列 中元素匹配的可能性很高， 这样 就 很容易 造成序列中字符的大面积匹配。 如果 此时设置的阀值不高，就很 容易 造成短序列在参考序列中几轮比对周期内就可以找到策略中承认的 片段 位置， 造成 错配； 如果阀值 过高， 也就是 我们要求供体基因序列需要和参考序列有很高的相似性 ， 对 不匹配 字符 的宽容度很低， 就 很有可能出现 供体片段比对完整个参考序列也没能 找到 对应位置， 这是 因
11、为供体序列不可能和参考序列完全相似， 要 容许潜在的差异性， 者 中差异性 就 包括 本文 比较关注的 SNP 和 Indel。 在 NGS 中， Indel 是仅次于 SNP 的变异类型， SNP（即： Single Nucleotide Polymorphisms， 单个碱基 突变） 广泛 存在 的 生物 体 基因变异， 目前的研究并没有证明 SNP 有 致病性 和致 缺陷 证明，它一般只能被认为是个体 间差异 的 基因 表达。 与 SNP 不同， Indel（ Insertion & Deletion， 即插入 与删除）是出现频率低于 SNP 的变异类型， 它 的出现 往往 导致生物体会
12、发生比较大的变动 ，有好的方面， 比如 ： 身高 比上一代人高出很多、牙齿防蛀能力很强； 但是 不幸的 是 大多数情况下都是坏的方面， 大量 Indel 很 容易 造成一定的功能缺陷，比如： 白化病 ， 癌症 等。 那么针对不同生物 体 基因序列，如果我们 对 对 它们 进行分类划分， 基本上分为 强变异型和弱变异型两类， 分别 对其 序列 进行比对测试， 从中得出 匹配 效果差异。 2.3.1 不同 阀值下的比对结果测试 除了 试验样本 分为 两个大类分别统计， 其他 因素尽量控制为基本相同， 以期待 在 单个 变量下的序列比对效果， 测试 环境和样本控制如下： A:相同 的计算机环境和 S
13、-W 序列比对程序； B:相同的 参考序列； C:两组 不同的序列样本， 分别 为正常样本和变异样本，每个样本中的片段 个数 控制为 10000 个 ，尽量 均匀的分布在参考基因组的各个位置，每个片段的长度控制在 500bp 左右， 这些 样本我们 事先已经知道了他们所在序列的正确位置。 针对 两种不同样本， 分别 测试 其 在 不 同阀值下 的比对结果。 本文中的阀值指的是在两条序列比对 成功 的情况下， 必须 满足的对应位置碱基的匹配百分比，就是 匹配的字符占片段长度字符 数 的比例 。 因为 基因序列的特殊性以及数据量大的特性， 设置 的阀值较低可导致 匹配 失误率较高， 不能 很好的保
14、证比对结果的可参考性； 设置 较高的阀值又 会 造成大量片段脱靶（没有命中任何位置就返回）。经过 大量 实验 ，本文选用 95%、 96%、 97%、 97.5% 、 98%、 98.5%这 六个阀值， 得出 的比对正确率6 / 8 如 图 3 所示： 图 3.不同 阀值下的命中率对比 从 上图 的 数据中我们可以得到如下结论： 第一 ： 总体 命中 率 上来看， 变异样本的命中率 整体上低于正常样本 ， 这个 不难理解， 变异 样本相对于正常样本 Indel 比例更高， 脱靶 可能性高于正常样本 ； 第二 ： 从 命中趋势上看， 两组 样本都表现为 类似于 驼峰 的比例分布， 中间 的命中比
15、例 高于 两端的命中 率 。出现这种情况的原因可能在于 阀值 较低是会 造成一定数量 的片段提前匹配 从而形成错配 ，导致 匹配生命期的 过早 结束；阀值 较高 时 ， 会造成脱靶 的可能性 变大 ， 即 寻找整个序列都没有找 可能 的位置 ； 第三 ： 从 命中率的变化幅度上看， 两组 数据都在阀值为 98%-98.5%出现 断崖式下跌， 这 说明了两组 样本 相对 于参考序列来说都存在一定比例的 SNP， 阀值 98.5%是一个较高的阀值， 容忍度 相对较低， 此时的 命中 率 下降说明了样本片段中普遍存在的 SNP 比例超过了 算法 的容忍度， 从而 造成了大量脱靶现象， 使得 命中率下
16、降。 实验过程 中 ，将阀值 设置 为最高时， 程序 的执行时间相对变长，侧面 也佐证了脱靶的猜测。 2.3.2 不同罚分 策略下的比对结果测试 由于 两种样本有一定差异， 变异样本 中存在一定 Indel 变异 ， 那么 在进行算法设计 时 ， 设置相同的流水线和优先匹配策略 可能会 不合时宜 。 为了 验证上面的猜测， 本节中将 测试不同罚分 策略 下分别在两组样本中设置不同优先级的情况下 的 比对情况对比 ， 在 正确比对为前提的情况下， 探索 不 同罚分策略下对比对判断的倾向性的影响 。 实验 前， 我们 已知的 样本 Indel 所在位点 ， 选取了 变异样本中特定的样本 集合 进行
17、试验， 样本 所包含 的 位点位于 hg19-chr13 的 序列上 ，与之相对应的是包含可疑 位点 的正常样本 短序列 的集合， 正常样本 和 变异样本 中 都存在包含此位点 的 多个短序列。 7 / 8 为了 保证实验的有更多的可对比性， 我们 对可疑位点的 Indel 进行了 认为 的干预， 分别 设置为10+、 20+、 30+、 40+、 50+的 五类 Insertion（ Deletion 的情况大体相同， 不做更多 介绍）， 分别 观察实验结果 实验过程 中， 保证 实验环境如下： a:包含 此位点的 reads 控制为 200 个 ， 正常样本 和变异样本各一半； b:相同
18、的算法策略、匹配 得分 （ 2 分）、 错配 罚分（ -1）， 不同 的 Indel 罚分， 分别为 ： -1,-2,-3,-4； c:位点 Insetion 控制为 10+、 20+、 30+、 40+、 50+五类； d:其他 环境基本相同； 实验中 ， 根据 Indel 大小不同划分了从 10+到 50 不同 大小的 Indel 分类 。 实验流程中 ， 每一次对一条的序列进行比对时， 罚分 从 -1 开始 递减， 每次 递减 0.1， 直到 设置的最大罚分（ -3） 为止 ， 一旦 达到正确比对结果就不再递减， 记录 当前罚分并更改当前罚分通过百分比， 当 当前 Indel 大小的通过
19、比例 100%时 ， 确定 全部通过情况下的最低罚分为 Indel 罚分阀值。 以上 的实验针对变异样本， 正常样本由于 Indel 长度 一般 不超过 10，不对其进行测试， 设置 Indel 罚分为 -1 时绝大多数 reads 都能通过测试，所以罚分都设置为 1， 只做对比使用， 不具有 实际意义。 实验 结果如 图 4 所示。 图 4.不同 Indel 长度下成功判断的平均阀值对比 从上图 中能够得出以下几点结论： 第一 ： 当 Indel 长度 越大时 ， 罚分 要求降低， 也就是 说， 如果 Indel 的目标足够大， 设置比较低的罚分就可以 对 其的筛查； 反之 ， 如果 Ind
20、el 长度不足够长， 我们 建议 要 相应的提高罚分 值 ， 这样 可以增大 Indel 检测 可能性； 第二 ： 从 阀值变化趋势上看， 当 Indel 长度变长时， 平均 阀值加快降低。 Indel 长度在 10+ 30+时 ， 阀值还 徘徊在 2.5 左右 的水平， 当 Indel 长度在 40+ 50+时 ， 阀值 快速回落到 1.5 左右 的水平。也就是 说， Indel 长度越大， 检测出 Indel 需要的罚分成本 加速 减少。 结束语： 本文 以算法为背景， 在 考虑了 分子 生物 学 应用的基础上， 探究 更改算法 策略 对 模式匹配结果 的8 / 8 影响， 进而 给出 该
21、影响下对 NGS 中 还原 样本 原始序列 和 Indel 判断的一些影响。 由于 实验数据有限，得出的结果 和 实际 业务数据量庞大 的 背景下得出的数据并不能保证一致， 但 大体趋势基本不变。 本文 通过还原序列命中率和 Indel判断命中率的两个实验，希望 能 给 NGS流程优化 提供 有价值 的参考。 参考文献 ： 1 陈科 , 朱清新 , 杨曦 . 最优 搜索机制下寻找最优插入 - 删除 种子 J. 电子科技大学学报 ,2011,40(2):292-295 2 汪东 ,唐敏 .Smith-Waterman 算法 在脉动阵列 上 的实现及分析 J.计算科学报 ,2004,27(1):1
22、2-20 3 张洪礼 ,燕翠霞 ,王常武 .多 序列比对问题的概率统计粒子群算法求解 J.小型微型计算机系统 ,2010,31(9):1833-1837 4 尚 婧 .下一代测序短序列比对软件算法比较及评价 ，苏州： 苏州大学 ， 2013 5 Mills R E,walter K,Stewart C,et al.Mapping copy number variation by population-scale genome sequencingJ.Nature,2011,470(7332):59-65. 6 Mills R E,Walter K,Stewart C,et al.Mapping
23、 copy number wariation by population-scale genome sequencingJ.Nature,2011,470(7332):59-65. 7 Li H,Ruan J,Durbin R.Mapping short DNA sequencing reads and calling variants using mapping quality scores J.Genome Res,2008,18(11):1851-8. 8 Chen K Y,Chao K M.A fully compressed algorithm for computing the e
24、dit distance of run-length encoded strings.Algorithmica,2013,65(2):354-370 9 张法 ， 乔香珍 ， 刘志勇 。 基于 Smith-Waterman 算法的并行分而治之生物序列比对算法。 中国 科学 E 辑技术科学， 2004， 34（ 2）： 190-199 10 Khajeh-Saeed A,Poole S,Blair Perot J.Acceleration of the Smith-Waterman algorithm using single and nultiple graphics processors.Journal of Compputational Physics,2010,229(11):4247-4258
