CN114550820B - 一种基于WFA算法的第三代测序RNA-seq比对方法 - Google Patents

Abstract

本发明公开了一种基于WFA算法的第三代测序RNA‑seq比对方法，该方法通过获取包含目标序列和查询序列的数据集；然后对参考基因组进行索引；以及区域选择和图映射；查找k长子字符串的最长公共子序列LCSk，再进行锚点过滤和锚点比对；引入注释文件获得比对结果，最后对比对结果进行评估。实验证明，本发明的方法有效地提高了序列比对的准确性，尤其是剪接位点比对的准确率，同时一定程度减少比对时间。

Description

一种基于WFA算法的第三代测序RNA-seq比对方法

技术领域

本发明涉及生物信息学中序列比对技术领域，具体是一种基于WFA算法的第三代测序RNA-seq比对方法。

背景技术

第三代RNA测序技术能确定几万到十几万碱基长度的RNA序列，因而被广泛地应用到转录体的检测、基因表达量估计以及剪接异构体的识别等转录组的相关研究中。然而第三代测序技术测序长度太长，且RNA-seq存在剪接位点，这都为算法设计带来了挑战。

现有的第三代RNA-seq比对算法主要有STAR、BBMap、GMAP、Minimap2、GraphMap2。这些方法都使用种子-扩展策略，且都有处理剪接连接的能力。但是，大多数算法准确率较低，不能很好地处理短外显子，且剪接识别准确率较低。因此，我们提出一种新的针对第三代测序技术的RNA-seq的比对算法，使用WFA算法进行精确比对，该算法利用序列相似性，沿对角线进行比对，从而提高比对的效率；同时，引入注释文件准确识别剪接位点，进而提高剪接位点识别准确率。

发明内容

本发明的目的在于解决三代测序RNA-seq比对算法准确率问题，而提供一种基于WFA算法的第三代测序RNA-seq比对方法，该方法通过WFA算法降低算法时间成本，通过引入注释文件提高比对算法准确率，尤其是在剪接位点方面。

实现本发明目的的技术方案是：

一种基于WFA算法的第三代测序RNA-seq比对方法，包括如下步骤：

1)获取数据集，数据集包括4个模拟数据集和3个真实数据集，每个数据集包括目标序列和查询序列，其中目标序列为参考基因组序列，记为ref；查询序列为第三代测序的RNA-seq读段，记为reads；

2)参考基因组ref为DNA序列，采用minimizer索引方法对参考基因组ref进行索引，将参考基因组序列ref的minimizer存储到哈希表中，其中minimizer为一段读段内最小哈希值的种子，然后将第三代测序RNA-seq序列与参考基因组的minimizer进行比对找到其所有的minimizer，通过哈希表找到其在参考基因组的位置，记为锚点；minimizer种子的最小哈希值计算公式如下：

其中序列s＝a₁a₂…a_n，长度|s|＝n，序列s的反补 函数π：∑^*×{0，1}→∑^*，使得π(s，0)＝s，/>

假设且对于k-mer s＝a₁…a_k的哈希函数/>定义为：

3)根据步骤2)产生的索引位置，将读段reads与ref的索引进行比对，初步确定读段reads在参考基因组ref上的比对区域，并将比对的结果聚类成候选区域；

4)通过读段reads的种子匹配构建对比链或锚点细化步骤3)得到的候选区域；

5)查找k长子字符串的最长公共子序列LCSk：每个锚点表示目标ref与查询序列reads之间的共享段，两个序列的起始和结束位置都是已知的，由于重复锚点的存在，步骤4)获得的锚点集不一定是单调递增的，因此选择满足单调性条件的锚点子集表示为在长度为k的子字符串中寻找最长公共子序列；

6)对步骤5)产生的锚点子集中的每个锚点进行过滤，使用背包算法的变体寻找最优锚点集合，具体是：

将参考基因组锚的起始和结束位置分别表示为x_s和x_e，读段起始和结束位置分别表示为y_s和y_e，把每一个锚点看作一条二维的直线，直线的起点在这条直线上表示为T_s(x_s，y_s)，终点表示为T_e(x_e，y_e)，且每个锚点都有其适度f，即锚覆盖的碱基d的数量，将寻找最优锚点集的问题形式化如下：

从N锚点的集合C中，C_i＝(T_si，T_ei，d_i)∈C，目标是找到k锚点的最优集合C′＝{(T_s1，T_e1，d₁)，(T_s2，T_e2，d₂)，…，(T_sk，T_ek，d_k)}，其中T_e1≤T_s2，T_e2≤T_s3，…，T_e(k-1)≤T_sk，且最大化/>

该问题等价于0-1背包问题，用相同的算法求解，0-1背包问题描述为：

对于给定的N元素集合E，其中每个元素e都有其权重w和适应度f，且有最大权重T的限制，需要找到其权重和不大于T，适应度和存在最大值，公式：

对于当前的锚点集，元素的权重是在锚点|y_e-y_s|中的查询长度，锚的适应度是覆盖碱基d的数量，总权重T是处理过的读段的总长度，通过优化问题和寻找到构建读段比对的锚的结果集x_i.N∈{0，1}；

7)使用WFA算法对步骤6)得到的最优锚点集中的锚点进行比对，将动态规划对角线得分递增的思想扩展到间隙-放射模型中，利用序列之间的同源性优化比对过程；具体是：

最远到达f.r.点表示在对角线为k和分数为s时的动态规划单元，距离对角线的起点更远，对于三个SWG(Smith Waterman Gotoh，SWG)矩阵I、D、M，定义/>为在对角线上到f.r.点/>的偏移；因此，对于给定的分数，定义s-波前WF_s为在分数为s时，所有f.r.点的集合，目标是计算最小s，使得WF_s的任何f.r.点到达(m，n)，对于每个s，分数只能来自s-o，s-e，s-x或沿对角线的前一个匹配分数；/>的表达式为：

只考虑插入、删除和不匹配，定义到f.r.点的偏移公式如(4)所示，且初始化

8)在步骤7)进行序列比对时，引入注释文件，进一步识别剪接位点，得到比对结果。

9)利用步骤1)中的4个模拟数据集和3个真实数据集对步骤8)得到的比对结果进行评估，模拟数据集主要使用了比对率(Aligned)、匹配率(Match rate)、正确率(Correct)、外显子全部覆盖率(HitAll)、至少覆盖一个外显子率(HitOne)、剪接读段(Split reads)、正确的剪接读段(Correct，split)、外显子全部覆盖的剪接读段(Splithit all)、至少覆盖一个外显子的剪接读段(Split hit one)评估指标进行评估；真实数据集主要使用了比对率(Aligned)、匹配率(Match rate)、表达基因数量(No.expressedgenes)、外显子覆盖率(Exon hit)、连续的比对率(Contiguous alignment)评估指标进行评估。

本发明提供的一种基于WFA算法的第三代测序RNA-seq比对方法，该方法通过WFA精确的空位仿射算法，利用序列之间的同源区域来加速比对过程，此算法在使用较少内存的同时，性能优于其他最先进的方法。

(2)引入注释文件可以更为精确地识别剪接位点，从而提高比对的准确率，尤其是剪接位点识别的准确率。

附图说明

图1为一种基于WFA算法的第三代测序RNA-seq比对方法的流程图；

图2为WFA算法流程图；

图3为识别剪接位点流程图；

图4为酿酒酵母比对结果图。

具体实施方式

下面结合附图和实施例对本发明内容做进一步阐述，但不是对本发明的限定。

实施例：

一种基于WFA算法的第三代测序RNA-seq比对方法，包括如下步骤：

(1)获取评估所需数据集，包括4个模拟数据集和3个真实数据集，4个模拟数据集使用PBSIM工具模拟生成，分别模拟了PacBio ROI酿酒酵母、黑腹果蝇、人类第19号染色体和ONT R2 2D黑腹果蝇数据；3个真实数据集包括RacBio ROI黑腹果蝇、纠错后的PacBioROI黑腹果蝇和PacBio子读段的黑腹果蝇。这里，使用不同物种、不同技术下的数据，更为全面地评估算法的性能。表1表示的是第三代测序RNA-seq读段reads的信息，针对不同物种，每一个物种对应一个参考基因组ref，使用ref和reads进行比对，找到reads在ref中的位置；

表1数据集

(2)使用minimizer的索引方法对参考基因组进行索引。基于minimizer的索引方法是将参考基因组序列的minimizer存储到哈希表中，其中，minimizer是一段读段内最小哈希值的种子。minimizer计算公式如(1)所示：

其中，对于序列s＝a₁a₂…a_n，长度|s|＝n，序列s的反补 函数π：∑^*×{0，1}→∑^*，使得π(s，0)＝s，/>

假设且对于k-mer s＝a₁…a_k的哈希函数/>定义为：

(3)区域选择。区域选择依赖于在查询序列reads和参考序列ref之间找到种子，然后将它们聚类为候选区域。种子查找采用了一种间隔种子的形式，类似于Levenshtein距离的空位q-gram滤波器。空位q-gram是一种播种策略，它允许对不精确匹配进行快速和敏感的查找，并允许在预定义的种子“不关心”(DC)位置上进行变化。有空位的q-gram允许形状中的DC位置也包含插入和删除。实现Levenshtein空位q-gram是基于构建参考序列的散列索引，其中q-gram位置通过形状布局所构造的关键字进行散列，仅仅包含构建关键的碱基，而简单地跳过DC碱基。

为了推导出种子命中的一般方法，借鉴了Hough变换(HT)的概念，这是图像处理中常用的一种方法，用于检测线、圆和椭圆等形状。HT定义了从图像点到累加器空间的映射，形成霍夫空间。在直线检测的情况下，如果笛卡尔空间中给定的一组点是共线的，那么它们之间的关系可以用具有公共斜率m和截距c的线性方程来表示：

y＝mx+c (5)

其中(x，y)是点在2D空间中的坐标。HT尝试确定描述给定点集合直线的参数m和c。然而，HT不是只检测一条最佳的线，公式(6)可以在参数空间中转化为它的对偶：

c＝-mx+y (6)

HT定义了一个累加器空间，在该空间中，其中m和c被栅格化，以便只取有限范围的值。然后，HT通过追踪笛卡尔空间中每个点的所有对偶线，并增加每个(m，c)坐标的投票计数，简单地计算累加空间中的所有潜在解。所有在定义阈值以上的HT空间坐标都可以被认为是原始笛卡尔空间中的候选线。

由于c对应于参考序列上的坐标，一个简单的整数数组的参考长度可用于计票计数。对于每个k点(单粒种子命中)，它的c参数值是用一个简单的表达式确定的：

c＝t-q (7)

(4)图映射。这一步骤，通过从段种子匹配中构建比对链或锚点来细化上一步的候选区域。我们使用了“kmer图映射”的概念。给定一序列(目标和查询序列)，首先从目标序列构造一个kmer映射图。对于图的构建，使用从目标动态构建的索引，为了灵敏度使用较小的连续种子。

(5)查找k长子字符串的最长公共子序列。每个锚表示目标和查询序列之间的一个共享段，两个序列的起始和结束位置都是已知的。由于重复锚点的存在，所获得的锚的集合在目标坐标和查询坐标中并不一定是单调递增的。因此，需要选择满足单调性条件的锚的子集。识别一个这样子集的问题可以表示为k长子字符串的最长公共子序列。在算法中，允许可变长度的子字符串，而不是使用固定大小k的子字符串，每个子串的大小等于两个序列中对应锚的长度。

(6)锚点过滤。完成步骤5)LCSk后，对生成的锚点进行处理，以获得用于构建比对的最优锚点集合。参考基因组锚的起始和结束位置分别表示为x_s和x_e，读段起始和结束位置分别表示为y_s和y_e，我们把每一个锚点看作一条二维的直线，直线的起点在这条直线上表示为T_s(x_s，y_s)，终点表示为T_e(x_e，y_e)。且每个锚点都有其适度f，最简单的形式就说锚覆盖的碱基d的数量。将寻找最优锚点集的问题形式化如下：从N锚点的集合C中，C_i＝(T_si，T_ei，d_i)∈C，目标是找到k锚点的最优集合

C′＝{(T_s1，T_e1，d₁)，(T_s2，T_e2，d₂)，…，(T_sk，T_ek，d_k)}，

这里，T_e1≤T_s2，T_e2≤T_s3，…，T_e(k-1)≤T_sk， (8)

且尝试最大化

这样定义的问题等价于0-1背包问题，可以用相同的算法求解。0-1背包问题可以描述为：对于给定的N元素集合E，其中每个元素e都有其权重w和适应度f，且有最大权重T的限制，需要找到其权重和不大于T，适应度和存在最大值，公式：

对于当前的锚点集，元素的权重是在锚点|y_e-y_s|中的查询长度，锚的适应度是覆盖剪辑d的数量。总权重T是处理过的读段的总长度。通过解决这个优化问题和寻找，我们找到了构建读段比对的锚的结果集x_i.N∈{0，1}。

(7)锚点比对。在找到读段在参考基因组上的大概位置后，需要进一步对其进行更精细的比对，从而得到更为精确的结果。这里主要使用WFA(wavefront alignment)算法进行更为精细的比对，但是由于RNA-seq存在剪接位点的问题，因此，引入注释文件识别剪接位点，从而提高比对的准确性。

WFA算法是将动态规划对角线得分递增等思想扩展到间隙-仿射模型中，利用序列之间的同源性来达到更好的比对结果。

最远到达(furthest-reaching，表示为f.r.)点表示在对角线为k和分数为s时的动态规划单元，距离对角线的起点更远。对于三个SWG(Smith-Waterman-Gotoh)矩阵I，D，M，定义/>为在对角线上到f.r.点/>的偏移。因此，对于给定的分数，定义s-波前(WF_s)为在分数为s时，所有f.r.点的集合。我们的目标是计算最小s，使得WF_s的任何f.r.点到达(m，n)。对于每个s，分数只能来自s-o，s-e，s-x或沿对角线的前一个匹配分数。

这里，只考虑插入、删除和不匹配，定义到f.r.的偏移公式如(4)所示，且初始化

WFA算法逐步计算增长分数的波前，直到到达(n，m)。算法流程图如图2所示，首先，对于每个分数s，使用波前扩展算法，沿着对角线跟踪匹配字符来扩展点然后，检查波前WF_s的每个f.r.点是否到达(n，m)。如果没有，继续计算下一个波前WF_s+1，然后再一次迭代。最后，通过波前的偏移来进行回溯，从而得到比对结果。

(8)引入注释文件由于剪接位点的存在，增加了读段与参考基因组进行比对的难度。因此，引入注释文件来识别剪接位点，从而提高比对准确率。识别剪接位点流程如图3所示。注释文件中包含参考基因组外显子起始和终止的位置信息，由此，可以找到剪接位点的位置与内含子的起始和终止位置。在两个锚点进行比对时，根据注释文件提供的信息，识别是否含有内含子，如果存在，就跳过内含子，使用WFA算法进行比对；如果没有，就直接使用WFA算法进行比对。这样，可以更准确地进行比对。

(9)评估结果

在模拟数据集上将我们的算法与四个现有的表现较好的算法(BBMap、GMAP、Minimap2和GraphMap2)进行比较，并生成评估结果。评估主要使用了比对率(Aligned)、正确率(Correct)、外显子全部覆盖率(HitAll)和至少覆盖一个外显子率(HitOne)四个指标进行评估。这里使用酿酒酵母数据作为例子，可以从图4中看出我们的算法表现出了较好的结果，四个指标分别达到98.08％、85.41％、95.56％和96.83％。说明算法在数据集上可以表现出比较好的比对结果。

将四个模拟数据集按照表1顺序分别编号为S1、D2、H3和D4，表2我们的算法在这四个数据机上与GraphMap2、GraphMap2_origin和Minimap2的运行时间比较结果，其中GraphMap2_origin是GraphMap2算法中的ksw_extd2_sse函数未使用SSE加速的情况。从结果可以看出，我们的方法在S1数据集上比GraphMap2运行时间短，在H3数据集上比GraphMap2_origin运行时间短，但是在D2和D4数据集上运行时间较长。同时，可以清楚地看出Minimap2在运行时间上仍然具有很大的优势。

表2算法运行时间

Claims (1)

1.一种基于WFA算法的第三代测序RNA-seq比对方法，其特征在于，包括如下步骤：

1)获取数据集，数据集包括4个模拟数据集和3个真实数据集，每个数据集包括目标序列和查询序列，其中目标序列为参考基因组序列，记为ref；查询序列为第三代测序的RNA-seq读段，记为reads；

2)参考基因组ref为DNA序列，采用minimizer索引方法对参考基因组ref进行索引，将参考基因组序列ref的minimizer存储到哈希表中，其中minimizer为一段读段内最小哈希值的种子，然后将第三代测序RNA-seq序列与参考基因组的minimizer进行比对找到其所有的minimizer，通过哈希表找到其在参考基因组的位置，记为锚点；minimizer种子的最小哈希值计算公式如下：

其中序列s＝a₁a₂…a_n，长度|s|＝n，序列s的反补 函数π：∑^*×{0，1}→∑^*，使得π(s，0)＝s，/>

假设且对于k-mer s＝a₁…a_k的哈希函数定义为：

3)根据步骤2)产生的索引位置，将读段reads与ref的索引进行比对，初步确定读段reads在参考基因组ref上的比对区域，并将比对的结果聚类成候选区域；

4)通过读段reads的种子匹配构建对比链或锚点细化步骤3)得到的候选区域；

5)查找k长子字符串的最长公共子序列LCSk：每个锚点表示目标ref与查询序列reads之间的共享段，两个序列的起始和结束位置都是已知的，由于重复锚点的存在，步骤4)获得的锚点集不一定是单调递增的，因此选择满足单调性条件的锚点子集表示为在长度为k的子字符串中寻找最长公共子序列；

6)对步骤5)产生的锚点子集中的每个锚点进行过滤，使用背包算法的变体寻找最优锚点集合，具体是：

将参考基因组锚的起始和结束位置分别表示为x_s和x_e，读段起始和结束位置分别表示为y_s和y_e，把每一个锚点看作一条二维的直线，直线的起点在这条直线上表示为T_s(x_s，y_s)，终点表示为T_e(x_e，y_e)，且每个锚点都有其适度f，即锚覆盖的碱基d的数量，将寻找最优锚点集的问题形式化如下：

从N锚点的集合C中，C_i＝(T_si，T_ei，d_i)∈C，目标是找到k锚点的最优集合其中T_e1≤T_s2，T_e2≤T_s3，…，T_e(k-1)≤T_sk，且最大化/>

该问题等价于0-1背包问题，用相同的算法求解，0-1背包问题描述为：

对于给定的N元素集合E，其中每个元素e都有其权重w和适应度f，且有最大权重T的限制，需要找至其权重和不大于T，适应度和存在最大值，公式：

对于当前的锚点集，元素的权重是在锚点|y_e-y_s|中的查询长度，锚的适应度是覆盖碱基d的数量，总权重T是处理过的读段的总长度，通过优化问题和寻找到构建读段比对的锚的结果集x_i.N∈{0，1}；

7)使用WFA算法对步骤6)得到的最优锚点集中的锚点进行比对，将动态规划对角线得分递增的思想扩展到间隙-放射模型中，利用序列之间的同源性优化比对过程；具体是：

最远到达f.r.点表示在对角线为k和分数为s时的动态规划单元，距离对角线的起点更远，对于三个SWG(Smith Waterman Gotoh，SWG)矩阵I、D、M，定义/>为在对角线上到f.r.点/>的偏移；因此，对于给定的分数，定义s-波前WF_s为在分数为s时，所有f.r.点的集合，目标是计算最小s，使得WF_s的任何f.r.点到达(m，n)，对于每个s，分数只能来自s-o，s-e，s-x或沿对角线的前一个匹配分数；/>的表达式为：

只考虑插入、删除和不匹配，定义到f.r.点的偏移公式如(4)所示，且初始化

8)在步骤7)进行序列比对时，引入注释文件，进一步识别剪接位点，得到比对结果。