CN113555061B

CN113555061B - 一种无参考基因组的变异检测的数据工作流处理方法

Info

Publication number: CN113555061B
Application number: CN202110835302.0A
Authority: CN
Inventors: 刘珍; 刘志岩; 王海宁
Original assignee: Harbin Yinji Technology Co ltd
Current assignee: Harbin Yinji Technology Co ltd
Priority date: 2021-07-23
Filing date: 2021-07-23
Publication date: 2023-03-14
Anticipated expiration: 2041-07-23
Also published as: CN113555061A

Abstract

本发明提出了一种无参考基因组的变异检测的数据工作流处理方法，包括如下步骤：提取原始测序数据FASTQ文件，建立存储库，将FASTQ文件切割为多个数据工作流，并将所述多个数据工作流存储在存储库中；采用基于最近投影算法将多个数据工作流进行投影变换，得到多个虚拟投影序列；利用变换空间索引算法，对所述虚拟投影序列进行字符串的搜索，使用解码操作，查找出所述字符串在数据工作流上的测序序列上的位置；当字符串的搜索步骤中没有找到合适的区间的时候，修改查找的字符串以命中不同的区间，然后在新的区间上继续进行搜索操作。

Description

一种无参考基因组的变异检测的数据工作流处理方法

技术领域

本发明涉及生物信息领域，具体涉及一种无参考基因组的变异检测的数据工作流处理方法。

背景技术

新一代测序仪现在以可承受的价格产生大量测序数据。最近的系统在36h的单次运行中产生超过60亿个150个核苷酸长的序列，足以对20个完整的人类基因组进行测序。这为遗传疾病的诊断和个性化医疗的发展开辟了许多新的视角，旨在基于人基因组特异性调整治疗。然而，这也带来了新的挑战，特别是与处理大量数据工作流相关的成本。

现有技术中，对于变异检测的数据工作流处理方法分为基于参考基因组的方法或基于无参考基因组的方法。虽然基于参考基因组的方法具有良好的处理精度但处理速度仍比较慢，而基于无参考基因组的处理方法更快但具有更低的处理精度。

发明内容

为了解决上述现有技术中对于变异检测的数据工作流处理方法中存在的弊端，本发明提出了一种无参考基因组的变异检测的数据工作流处理方法。包括如下步骤：

S1、提取原始测序数据FASTQ文件，建立存储库，将FASTQ文件切割为多个数据工作流，并将所述多个数据工作流存储在存储库中；

S2、采用基于最近投影算法将多个数据工作流进行投影变换，得到多个虚拟投影序列；

S3、利用变换空间索引算法，对所述虚拟投影序列进行字符串的搜索，使用解码操作，查找出所述字符串在数据工作流上的测序序列上的位置；

S4、当字符串的搜索步骤中没有找到合适的区间的时候，修改查找的字符串以命中不同的区间，然后在新的区间上继续进行搜索操作。

进一步地，步骤S1中，采用多层级多节点分配和排序的方法，将存储库划分到多个区域，并将所述多个数据工作流发送到相应的多个区域，接着在多个区域上完成排序。

进一步地，步骤S4中，还包括差异查找步骤，获取搜索匹配的最小次数，利用搜索匹配的最小次数对差异读段区间，实现快速地搜索操作与解码操作。

进一步地，所述差异读段包括未投影和不完全投影的读段。

进一步地，步骤S2具体包括：

S21、假设每个数据工作流上的测序序列为T，长度为n，在测序序列T的末尾添加一个投影标识符，得到带有投影标识符的测序序列T'；

S22、对带有投影标识符的测序序列T’进行循环右移，每次右移一位，循环n次后得到n个长度为n的字符串，对n个长度为n的字符串按照字典顺序排序，得矩阵M；

进一步地，还原虚拟投影序列得到数据工作流上的测序序列的步骤如下:

S23、求矩阵M的第一列F列，按照字典顺序重新排列虚拟投影序列；

S24、从矩阵M的第一行开始，找到其F列对应的字符，并输出该字符。

S25、找到其前面的同一行L列的字符。

S26、根据L列的字符，找到与其在同一行的F列字符；

S27、重复步骤S24-S26，直到L列的字符为投影标识符，输出的结果即为数据工作流上的测序序列的逆序，再一次逆序输出结果，并去掉标识符，得数据工作流上的测序序列。

进一步地，对所述虚拟投影序列进行字符串的搜索，具体包括：

S31、字符串设为P，长度为m，P[i]表示P上的第i个字符，1≤i≤_m，设在虚拟投影序列上查找到区间为[start,end]；

S32、令初始i＝m，根据start＝C[P[i]]+1,end＝C[P[i]+1]计算最初的start与end，然后i自减1；

S33、从后往前搜索P，按照公式(1)、(2)进行计算获得start与end，

start＝C[P[i]]+Occ(P[i]，start-1)+1 (1)

end＝C[P[i]]+Occ(P[i]，end) (2)

S34、完成计算后i自减1，重新计算start与end，直到遍历P，即i为0；

S35、遍历结束后，如果end<start，则字符串P没有在虚拟投影序列中出现；否则，输出end-start+l作为出现次数；

S36、进行虚拟投影序列到测序序列的解码操作；从位置k开始执行解码操作，直到得到的字符为标识符，解码过程中得到的字符个数就是i。

附图说明

图1是本发明的无参考基因组的变异检测的数据工作流处理方法流程图；

图2是本发明步骤2中将多个数据工作流进行投影变换和还原变化的具体步骤流程图；

图3是本发明步骤3中利用变换空间索引算法，对虚拟投影序列进行字符串的搜索的具体步骤流程图。

具体实施方式

结合附图并参考以下详细说明，本发明各实施方式的特征、优点及其他方面将变得更加明显。虽然附图中显示了本方案的相对合适的实施方案，但部分细节并未清晰描述，接下来将结合附图对这个方案进行更加详尽的描述，以便将本方案更完整的传达给本领域的相关技术人员。

提取原始测序数据FASTQ文件，建立存储库，将FASTQ文件切割为多个数据工作流，并将这些数据工作流存储在存储库中。

具体地，采用多层级多节点分配和排序的方法，并调整数据切割和分配的顺序。首先，将存储库划分到多个区域，并将多个数据工作流发送到相应的多个区域，接着在多个区域上完成排序。

由于多个数据工作流的数据量庞大，给传输带来严峻的挑战，因此需要对多个数据工作流进行高效变换。

本发明采用基于最近投影算法进行数据空间变换：算法的构思为转换一个文本为另一虚拟投影文本，使得转换后相同的字符尽可能投影对应，转换后的结果比使用其他压缩算法实现更有效地压缩。

最近投影算法的具体步骤为：

建立多个虚拟投影序列，每个虚拟投影序列均为所对应的数据工作流按照投影规则进行投影所得；多个虚拟投影序列的长度和数量与多个数据工作流上的测序序列长度和数量一致；保存原始测序数据与其对应的虚拟投影序列之间的投影标识符。

具体地，假设输入的多个数据工作流的每个数据工作流上的测序序列为T，长度为n，在测序序列T的末尾添加一个投影标识符，得到带有投影标识符的测序序列T'；

对带有投影标识符的测序序列T’进行循环右移，每次右移一位，循环n次后得到n个长度为n的字符串；

对n个长度为n的字符串按照字典顺序排序，得矩阵M。

矩阵M的第一列为F列，最后一列为L列，L列即为变换后的虚拟投影序列B(T)。矩阵M有两个特性。一是在同一行中，F列的字符是L列的字符在原数据工作流上的测序序列T上的后一个字符；二是F列与L列中同一字符出现的相对顺序不变。

可逆地，还原虚拟投影序列B(T)得到数据工作流上的测序序列T的步骤如下:

1、求矩阵M的第一列F列，按照字典顺序重新排列虚拟投影序列B(T)即可。

2、从矩阵M的第一行开始，找到其F列对应的字符，并输出该字符。

3、应用特性一，找到其前面的一个字符，即同一行L列的字符。

4、应用特性二，根据L列的字符，找到与其在同一行的F列字符。

5、重复步骤2-4，直到L列的字符为投影标识符，输出的结果即为原数据工作流上的测序序列T的逆序，再一次逆序输出结果，并去掉标识符，即可得数据工作流上的测序序列T。

其次，利用变换空间索引算法，对变换后的虚拟投影序列B(T)进行字符串的搜索，为高效实现搜索操作，变换空间索引算法引入两个辅助查询表:C(c)以及Occ(c,k)，c表示数据工作流上的测序序列T中出现的字符；k表示虚拟投影序列B(T)中第k位；C(c)表记录字典顺序小于字符c的所有字符个数和；Occ(c,k)表记录字符从在L列的前缀L[1...k]出现的次数。

搜索操作以模式串为输入，查找模式串在数据工作流上的测序序列T上出现的次数，其搜索步骤如下:

(1)、模式串设为P，长度为m。P[i]表示P上的第i个字符，1≤i≤_m。设在虚拟投影序列B(T)上查找到区间为[start,end]。

(2)、初始化start与end。令初始i＝m，根据start＝C[P[i]]+1,end＝C[P[i]+1]计算最初的start与end，然后i自减1。

(3)、从后往前搜索P，根据P[i]计算start与end的值。即从i＝m-1开始，对于0≤i<m,按照公式(1)、(2)进行计算获得start与end。完成计算后i自减1，重新计算start与end，直到遍历P，即i为0。

start＝C[P[i]]+Occ(P[i]，start-1)+1 (1)

end＝C[P[i]]+Occ(P[i]，end) (2)

(4)、遍历结束后，如果end<start，则模式串P没有在虚拟投影序列B(T)中出现；否则，输出end-start+l作为出现次数。

(5)、进行虚拟投影序列B(T)到测序序列T的解码操作；接受L列的位置k作为输入，求出字符L[k]在原数据工作流上的测序序列T上的位置i。具体地，从位置k开始执行解码操作，直到得到的字符为标识符，解码过程中得到的字符个数就是i。

通过以上的搜索操作以及解码操作的结合，可以实现模式串在原数据工作流上的测序序列T上的准确查找。首先通过搜索操作找出模式串P在虚拟投影序列B(T)的区间，然后对于区间中每个字符，使用解码操作，查找出该字符在原数据工作流上的测序序列T上的位置。

由于在建立虚拟投影序列的过程中，会存在对应于“未投影”和“不完全投影”的读段情况，“不完全投影”指读段包含至少一个错配，而读段的至少一部分与虚拟投影序列的一部分匹配。

因此，在优选实施例中还包括对数据工作流测序序列的虚拟投影序列进行差异查找步骤，实际上是对搜索操作进行修改，添加查找失败时执行的回溯功能。

当搜索操作结束而没有找到合适的区间的时候，搜索算法回溯到上一步，并修改查找的字符以命中不同的区间，然后在新的区间上继续进行搜索操作。

进一步地，在新的区间上继续进行搜索操作中，需要获取搜索匹配的最小次数d，进一步设置差异查找步骤中，设搜索过程中进行搜索匹配的最小次数d，本实施例采用如下流程进行获取：

(1)、设有搜索字符串为A，长度为a；投影匹配字符串为B，长度为b；

(2)、构造行数是a+1，列数是b+1的矩阵M’。初始化矩阵M’，第零行元素为0到b，第零列元素为0到a。

(3)、M_h,g’表示矩阵上第h行第g列的元素，其中1≤h≤a,1≤g≤b。从第一行第一列开始，根据公式(3)求矩阵剩余元素的值。公式(3)表示M_h,g’的值是由该元素前一列元素的值加上1、该元素前一行元素的值加上1、以及该元素的对角线上的元素加上temp三者确定，三者中之中最小的为M_h,g’。其中，temp取决于字符串A上第h个字符与字符串B上第g个字符是否相等，相等则为0，不相等则为1。

M_h,g’＝min(M_h-1,g’+1,M_h,g-1’+1,M_h-1,g-1’+temp) (3)

(4)、矩阵填充完成后，输出元素M_a+1,b+1’，作为字符串A与B的搜索匹配次数。

(5)、从M_a+1,b+1开始回溯求A变换为B的字符串S_AB，以及B变换A的字符串S_BA。设一开始M_a+1,b+1为当前元素M_h,g’，寻找M_h,g’的前一列元素M_h,g-1’，前一行元素M_h-1,g’，对角线元素M_h-1,g-1’中最小的值。如果三者中最小值有相同，则按照对角、行、列的优先顺序选取。根据选择元素与原本元素的位置关系，在S_AB与S_BA从后往前添加元素。

如果选择M_h,g-1’，则S_AB加上字符A_i,S_BA加上字符“-”；如果选择M_h-1,g’，则S_AB加上字符“-”，S_BA加上字符B_j；如果选择M_h-1,g-1’，则S_AB加上字符A_i,S_BA加上字符B_j；

选择的元素作为下一个当前元素，重复执行该步，直到回溯到M_a,b，根据字符串A,S_AB,B,S_BA即可以知道搜索匹配操作的次数。利用该投影匹配操作的次数即可以对“未投影”和“不完全投影”的读段情况，获取差异查找的捷径，辅助快速地实现搜索操作与解码操作。

本发明通过上述数据工作流的处理过程，我们只需要对原始基因数据切割为多个数据工作流，并将这些数据工作流存储在存储库，对多个数据工作流进行高效变换，提高存储空间及传输效率，变换后的虚拟投影序列的实际内容相比原始基因数据工作流的序列会小很多；需要时，利用变换空间索引算法，对变换后的虚拟投影序列进行字符串的搜索和解码以及差异查找步骤，实现高效的字符串的定位。

本发明的数据工作流处理方法可以利用软件进行实现，当以软件实现时，计算机程序代码可以存储在计算机介质上并由包括一个或多个处理器的硬件处理单元执行。应当理解，本文使用的术语“处理器”旨在包括一个或多个处理设备，包括信号处理器、微处理器、微控制器、专用集成电路(ASIC)、现场可编程门阵列(FPGA)或其他类型的处理电路，以及这些电路元件的部分或组合。此外，如本文所使用的术语“存储器”旨在包括与处理器相关联的电子存储器，诸如随机存取存储器(RAM)、只读存储器(ROM)或其他类型的存储器，以任何组合。

因此，用于执行本文描述的方法和协议的软件指令或代码可以存储在一个或多个相关联的存储器设备(例如，ROM、固定或可移动存储器)中，并且当准备好被利用时，加载到RAM中并由处理器执行。

本公开的技术可以在各种各样的设备或装置中实现，包括例如移动电话、计算机、服务器、平板电脑和类似设备。

尽管本文已经参考附图描述了本发明的说明性实施例，但是应当理解，本发明不限于那些精确的实施例，并且在不脱离本发明的范围或精神的情况下，本领域技术人员可以进行各种其他改变和修改。

Claims

1.一种无参考基因组的变异检测的数据工作流处理方法，其特征在于，包括如下步骤：

S2、采用基于最近投影算法将多个数据工作流进行投影变换，得到多个虚拟投影序列；建立多个虚拟投影序列，每个虚拟投影序列均为所对应的数据工作流按照投影规则进行投影所得；多个虚拟投影序列的长度和数量与多个数据工作流上的测序序列长度和数量一致；保存原始测序数据与其对应的虚拟投影序列之间的投影标识符；

对所述虚拟投影序列进行字符串的搜索，具体包括：

S31、字符串设为P，长度为m，P[i]表示P上的第i个字符，1≤i≤_m，设在虚拟投影序列上查找到区间为[start, end]；

S32、令初始i=m，根据start= C[P[i]]+1, end= C[P[i]+1]计算最初的start与end，然后i自减1；

S33、从后往前搜索P，按照公式（1）、（2）进行计算获得start与end，

start= C[P[i]]+ Occ（P[i]，start-1）+1 （1）；

end= C[P[i]]+ Occ（P[i]，end）（2）；

其中，C[ ]和Occ（）为辅助查询表；

S35、遍历结束后，如果end < start，则字符串P没有在虚拟投影序列中出现；否则，输出end-start+l作为出现次数；

S36、进行虚拟投影序列到测序序列的解码操作；从位置k开始执行解码操作，直到得到的字符为标识符，解码过程中得到的字符个数就是i；

S4、当字符串的搜索步骤中没有找到合适的区间的时候，修改查找的字符串以命中不同的区间，然后在新的区间上继续进行搜索操作；还包括差异查找步骤，获取搜索匹配的最小次数，利用搜索匹配的最小次数对差异读段区间，差异读段包括未投影和不完全投影的读段。

2.根据权利要求1所述的数据工作流处理方法，其特征在于，步骤S1中，采用多层级多节点分配和排序的方法，将存储库划分到多个区域，并将所述多个数据工作流发送到相应的多个区域，接着在多个区域上完成排序。

3.根据权利要求1所述的数据工作流处理方法，其特征在于，步骤S2具体包括：

S22、对带有投影标识符的测序序列T’进行循环右移，每次右移一位，循环n次后得到n个长度为n的字符串，对n个长度为n的字符串按照字典顺序排序，得矩阵M。

4.根据权利要求3所述的数据工作流处理方法，其特征在于，还原虚拟投影序列得到数据工作流上的测序序列的步骤如下:

S24、从矩阵M的第一行开始，找到其F列对应的字符，并输出该字符；

S25、找到其前面的同一行L列的字符；

S26、根据L列的字符，找到与其在同一行的F列字符；

S27、重复步骤S24- S26，直到L列的字符为投影标识符，输出的结果为数据工作流上的测序序列的逆序，再一次逆序输出结果，并去掉标识符，得数据工作流上的测序序列。