1.本发明属于生物信息
技术领域：
：，尤其涉及一种生物测序序列快速修剪方法及系统。
背景技术：
：：2.本部分的陈述仅仅是提供了与本发明相关的
背景技术：
：信息，不必然构成在先技术。3.在新一代测序中，要进行测序的核酸序列与接头序列（adapter）连接，以便被测序仪识别，然而当核酸序列的长度短于测序平台运行的读取长度时，测序得到的基因序列片段（称之为read）将同时包含需要进行测序的核酸序列以及全部或者部分接头序列。除此之外，在ngs（nextgenerationsequencing：下一代测序）测序中，测序结果的可信度在末尾循环（tailcycle）会变得较低，得到的是一些低质量的测序序列。被测序接头或者低质量的测序过程污染的测序序列经常导致不能满意的下游分析（如基因比对工作等）结果，因此修剪（trim）测序中的接头以及低质量的数据成为下游分析任务之前不可缺少的环节。4.随着现代测序仪的进步，不断增长的吞吐量和序列长度，为修剪工作提出了新的挑战，当前的一些处理工具的性能（比如速度和准确率）已经难以满足如此规模的测序数据，使得预处理步骤成为了数据分析的瓶颈，一种用于ngs数据预处理的超快的、准确的测序数据接头质量修剪工具仍是迫切需求。5.目前有许多的工具被用于新一代测序数据中的测序接头和低质量数据的修剪工作，比如：trimmomatic、fastp、ktrim等，这些工具采用了不同的trim算法，实现了ngs测序数据中的adapter-trim和quality-trim。6.发明人发现，ktrim相比其他几种工具，在ngs短读测序数据上有更好的表现，但是，随着固态硬盘技术的发展和普及，以及磁盘阵列技术的进步，ktrim当前的处理速率距离硬盘的读写峰值差距较大，无法满足现在对生物数据预处理的性能要求，而且经过测试，ktrim的线程扩展性并不好，使用四个以上的线程数量程序的处理速度难以继续提高。技术实现要素：7.本发明为了解决上述问题，提供了一种生物测序序列快速修剪方法及系统，所述方案通过采用轻量级的i/o框架，大幅度提升了i/o的速率；同时，在此基础上采用了向量化的方式优化了adapter的寻找过程，有效提高了生物测序序列快速修剪的处理效率。8.根据本发明实施例的第一个方面，提供了一种生物测序序列快速修剪方法，包括：获取待修剪的生物测序序列；对所述生物测序序列进行读操作、修剪操作以及写操作；其中，基于生产者—消费者模型对所述读操作、修剪操作以及写操作进行解耦，实现异步执行；且所述生物测序序列的格式化过程从读操作中转移到修剪操作中。9.进一步的，所述读操作、修剪操作以及写操作分别采用独立的线程进行实现，其中，读线程和写线程均设置有一个，所述修剪线程设置有一个或多个。10.进一步的，所述读操作用于通过读线程对所述生物测序序列按照块方式进行读取，并将读取的块对象存储入第一数据队列中。11.进一步的，所述修剪操作用于通过修剪线程从所述第一数据队列中获取数据，对所述生物测序序列进行格式化，去除生物测序序列中低质量碱基序列和接头序列；同时将处理后的序列存储入第二数据队列中。12.进一步的，在所述修剪线程中获取接头序列包括：将所述生物测序序列中的每个碱基作为一个字符；基于向量寄存器，采用若干次位运算获得预设长度序列数据中接头序列的位置。13.进一步的，所述写操作用于通过写线程从所述第二数据队列中获取处理后的生物测序序列，并进行存储。14.根据本发明实施例的第二个方面，提供了一种生物测序序列快速修剪系统，包括：数据获取单元，其用于获取待修剪的生物测序序列；数据处理单元，其用于对所述生物测序序列进行读操作、修剪操作以及写操作；其中，基于生产者—消费者模型对所述读操作、修剪操作以及写操作进行解耦，实现异步执行；且所述生物测序序列的格式化过程从读操作中转移到修剪操作中。15.与现有技术相比，本发明的有益效果是：（1）本发明提供了一种生物测序序列快速修剪方法及系统，所述方案通过采用轻量级的i/o框架，大幅度提升了i/o的速率；同时，在此基础上采用了向量化的方式优化了adapter的寻找过程，有效提高了生物测序序列快速修剪的处理效率。16.（2）本发明解决了ktrim的io效率低的问题，将文件数据的输入输出速率改进达到或接近磁盘读写的性能峰值。所述方案通过采用生产者-消费者模型实现读线程、工作线程以及写线程之间解耦、异步和速度平衡；同时，将数据格式化的任务从读线程转移到工作线程，将写回数据准备的任务从写线程转移到工作线程来最大可能地保证读线程和写线程的速度。17.（3）本发明所述方案使用数据池datapool重复利用chunk对象，减少对象的创建与销毁，减少创建对象的开销；同时，在读线程、工作线程以及写线程之间传递数据时，尽可能使用数据的指针，而不进行数据复制，有效减少了内存拷贝的开销；（4）本发明所述方案在整个工作流程中的读线程、工作线程（即修剪线程）以及写线程均只创建一次，直到完成其所有的任务后才销毁，有效减少了线程创建的开销。18.本发明附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。附图说明19.构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。20.图1为本发明实施例中所述的读线程和工作线程（即trim线程）之间chunk对象流转示意图；图2为本发明实施例中所述的工作线程（即trim线程）和写线程之间的关系示意图；图3为本发明实施例中所述的轻量级io框架结构示意图；图4为本发明实施例中所述的整体框架示意图。具体实施方式21.下面结合附图与实施例对本发明做进一步说明。22.应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属
技术领域：
：的普通技术人员通常理解的相同含义。23.需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。24.在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。25.术语解释：adapter：接头序列，是指在生物序列测序过程中因测序技术或者测序平台原因测得的额外序列片段。26.trim：测序数据质量控制的重要环节，主要是去除低质量碱基序列和adapter序列。27.实施例一：本实施例的目的是提供一种生物测序序列快速修剪方法。28.首先，需要说明的是，本发明是以ktrim工具为基础进行的改进，ktrim相比现有的其他几种工具（如trimmomatic、fastp），在ngs短读测序数据上有更好的表现，通过在配备英特尔至强cpu和标准64位centos系统的计算服务器设备上进行测试，使用单端测试数据（一千万条fastq格式的测序数据，大小约2.1gb），在四个线程上该工具可以达到其最佳性能，处理时间是7.26秒，处理速率大约是300mb/秒；使用双端测试数据（两个单端测序数据文件，共两千万条fastq格式的数据，大小约4.2gb），在四个线程上达到最佳性能，处理时间是11.29秒，最佳性能处理速率大约是370mb/秒。ktrim工具的作者在其论文（sunk.ktrim:anextra-fastandaccurateadapter-andquality-trimmerforsequencingdata）中提到ktrim在保持较高的准确率的同时，处理速度是其它几种工具（trimgalore、trimmomatic、seqpurgea）的~2-18倍。29.尽管，ktrim相比其它工具有着更好的性能，但是随着固态硬盘技术的发展和普及，以及磁盘阵列技术的进步，ktrim当前的处理速率距离硬盘的读写峰值差距较大，无法满足当前对生物数据预处理的性能要求，而且经过测试，ktrim的线程扩展性并不好，使用四个以上的线程数量程序的处理速度难以继续提高，因此，ktrim在计算性能上仍有较大的进步空间。30.通过分析和测试，我们发现ktrim的性能瓶颈主要在fastq文件的解析上，低效率的io方式导致程序的处理性能收到了限制，是典型的的io密集型应用，同时，在ktrim中寻找adapter的过程也相对耗时，针对以上情况，我们提出了一种生物测序序列快速修剪方法（以下称之为rabbittrim工具），该工具在ktrim的基础上使用了轻量级的fastq数据io框架（如图3所示，展示了该轻量级io框架结构图），大幅度提升io的速率，在此基础上又使用向量化的方式优化了adapter（接头序列）的寻找过程，使得程序的处理效率得到了明显的提升，可接近我们测试平台的硬盘读写峰值（约2gb/秒）。31.具体的，ktrim工具的性能瓶颈主要在数据文件的io上，包括：（1）读线程压力大。在ktrim中使用一个读线程读取输入数据文件，读线程不仅需要读取字符串数据，还需要按照fastq格式进行数据解析，数据格式化的工作大大加重了读线程的负担，导致生产fastq格式的数据效率很低。32.（2）数据预取数量少。ktrim支持多线程模式，但在多线程模式下也仅是会多预取一块数据（一个batch数据），当计算部分执行较快时输入数据供给依然存在问题。33.（3）工作流程上各部分间存在依赖关系。读操作和trim操作（即修剪操作）之间实现了部分并行，但当两个操作的执行速率不匹配的时候，会造成一方在等待另一方的问题，而trim操作和写操作之间是串行的方式，处理完一块数据并输出完成后才会处理下一块，读操作、trim操作、写操作之间存在较大的依赖。34.（4）寻找接头的过程较慢，ktrim工具使用《string.h》库中的char*strstr(constchar*haystack,constchar*needle)函数来寻找接头的“种子”，这种方式相对较慢。35.为了解决上述问题，本发明提供了一种生物测序序列快速修剪方法，包括：获取待修剪的生物测序序列；对所述生物测序序列进行读操作、修剪操作以及写操作；其中，基于生产者—消费者模型对所述读操作、修剪操作以及写操作进行解耦，实现异步执行；且所述生物测序序列的格式化过程从读操作中转移到修剪操作中。36.进一步的，所述读操作、修剪操作以及写操作分别采用独立的线程进行实现，其中，读线程和写线程均设置有一个，所述修剪线程设置有一个或多个。37.进一步的，所述读操作用于通过读线程对所述生物测序序列按照块方式进行读取，并将读取的块对象存储入第一数据队列中。38.进一步的，所述块对象的创建引入数据池思想，仅创建预设数量的块对象进行重复使用。39.进一步的，所述修剪操作用于通过修剪线程从所述第一数据队列中获取数据，对所述生物测序序列进行格式化，去除生物测序序列中低质量碱基序列和接头序列；同时将处理后的序列存储入第二数据队列中。40.进一步的，在所述修剪线程中获取接头序列包括：将所述生物测序序列中的每个碱基作为一个字符；基于向量寄存器，采用若干次位运算获得预设长度序列数据中接头序列的位置。41.进一步的，所述写操作用于通过写线程从所述第二数据队列中获取处理后的生物测序序列，并进行存储。42.进一步的，所述读操作、修剪操作以及写操作所对应的线程仅创建一次，直到处理任务完成后进行销毁。43.进一步的，所述格式化具体包括按照fastq格式进行数据解析。44.具体的，为了便于理解，以下结合附图对本发明所述方案进行详细说明：本发明所述方案主要针对ktrim工具中存在的问题进行改进，包括以下两个方面：输入输出过程的改进和adapter查询过程的改进，具体为：（一）输入输出过程的改进针对输入输出过程，本发明提出的rabbittrim工具使用了轻量级的io框架来完成数据的输入、解析和输出，其中，所述框架是基于生产者-消费者模型，将读操作、trim操作、写操作三者进行解耦，实现异步执行；具体的，其处理过程如下：在读操作和trim操作之间建立生产者-消费者模型，读线程作为生产者而trim线程作为消费者，为了保证读取数据的有序性和正确性，生产者设置为一个，而消费者可以设置为一个或多个。在上文提到，ktrim中读线程不仅负责读取数据也负责格式化数据，这给读线程带去了较大的压力。因为格式化数据时构造数据对象相比从文件中读取指定大小的字符串的速度是比较慢的，所以本发明将格式化的工作转移到了消费者（即trim线程上），有效提高生产者（即读线程）的效率；同时，本发明可以通过增加trim线程的数量来平衡消费者和消费者的速度。45.本发明所述方案按照块（chunk）的方式读取数据，读线程将数据读取到chunk对象中，然后将chunk对象存储进队列供trim线程使用，为了减少创建反复chunk对象和为其申请内存带来的开销，同时采用了数据池的思想，只创建固定数据的chunk对象反复使用。在读线程、trim线程之间chunk对象流转过程如图1所示。46.其次，在trim线程和写线程之间，同样建立了生产者-消费者模型，trim线程作为生产者，写线程是消费者，为了保证写回数据的正确性，程序中只使用了一个写线程，因为写磁盘的速度较慢，程序中将fastq对象数据转化回字符串数据的过程放在了trim线程中，其关系如图2所示。47.综上，针对输入输出过程，本发明所述方案相对于现有的ktrim工具，具有如下改进：（1）有效解决了ktrim的io效率低的问题，将文件数据的输入输出速率改进达到或接近磁盘读写的性能峰值：（2）使用生产者-消费者模型实现读线程、工作线程以及写线程之间解耦、异步和速度平衡；（3）将数据格式化的任务从读线程转移到工作线程，将写回数据准备的任务从写线程转移到工作线程来最大可能地保证读线程和写线程的速度；（4）使用数据池datapool重复利用chunk对象，减少了对象的创建与销毁，减少了创建对象的开销；（5）在读线程、工作线程（即trim线程）以及写线程之间传递数据时，尽可能使用数据的指针，而不进行数据复制，有效减少内存拷贝的开销；（6）整个工作流程中读线程、每个工作线程以及写线程只会创建一次，直到完成其所有的任务后才销毁，减少了线程创建的开销。48.（二）adapter查询过程的改进在ktrim工具中的核心部分为寻找可能的adapter的位置，即adapter的“seed”，它是一个长度为3的核苷酸序列，寻找“seed”的过程可以抽象为在一个长度为几十到几百的字符串中寻找一个长度为3的子串；现有的ktrim使用了《string.h》库中的char*strstr(constchar*haystack,constchar*needle)函数，在优化了io模块后，这种简单的方式效率低下，成为了程序运行时的热点。49.本发明所述方案基于现代计算机中的向量化部件，使用向量化的方式来寻找所有“seed”的位置，具体的：将测序序列中的每个碱基可以看成一个字符，我们使用avx512向量寄存器可以一次性处理64字节的数据，即64个碱基，通过avx512向量寄存器和指令集，使用多次位运算（三次异或运算、两次左移运算以及两次或运算）我们就可以得到长度为64的序列数据中“seed”所在的位置，我们将长度为l的测序的序列数据分为n组，每组包含64个碱基，进行上述位运算就可以得到整条序列中所有“seed”的位置。其中，符号表示取上整。其具体过程如下述伪代码所示：算法1vectorized-find-seed输入:seq测序序列，index1种子1，index2种子2，seed1种子1出现的所有位置，seed2种子2出现的所有位置，seedtable位置查询表输出:seed1andseed21：functionvectorized-find-seed(seq,index1,index2,seed1,seed2,seedtable)2:len←length(seq)▷序列长度3:num←ceil(len−2/62)▷需要迭代的次数4:i←05:v11←init(index1[0])▷将avx512寄存器以“种子”的某个碱基初始化6:v12←init(index1[1])7:v13←init(index1[2])8:v21←init(index2[0])9:v22←init(index2[1])10:v23←init(index2[2])11:whilei《numdo12:v1←load(seq+62∗i)▷加载第i次迭代需要的序列数据到avx512寄存器v1中13:res11←xor(v1,v11)▷异或运算14:res21←xor(v1,v21)15:v1←leftshift(v1,8)▷v1寄存器数据左移8bit16:res12←xor(v1,v12)17:res22←xor(v1,v22)18:res12←or(res11,res12)▷将两次异或的结果进行或运算19:res22←or(res21,res22)20:v1←leftshift(v1,8)21:res13←xor(v1,v13)22:res23←xor(v1,v23)23:res13←or(res12,res13)24:res23←or(res22,res23)25:res1←mask(res13)▷使用mask将或运算后的结果转化为mask6426:res2←mask(res23)27:j←028:whilej《8do29:t1←res1and255▷mask64的后8位表示的值30:t2←res2and25531:res1←rightshift(res1,8)▷mask64右移8bit32:res2←rightshift(res2,8)33:putseedtable[i][j][t1]intoseed1▷使用查表法去数组seed-table中查询“种子”位置34:putseedtable[i][j][t2]intoseed235:j++36:endwhile37:i++38:endwhile39:endfunction其中，按照62而不是64划分是因为我们要寻找的是长度为3的碱基序列，当一个64个碱基的序列最后两个（第63、64个）碱基与“seed”的前两个碱基完全匹配，但是由于下一个碱基数据没法读进来，最终得到的结果则认为不是“seed”的开始的位置，但是实际上如果下一个碱基恰好与“seed”的第三个碱基匹配，则该位置就是“seed”的位置。所以我们在进行划分时，需要考虑到上一个64个碱基序列的最后两位，下一个64个碱基序列的前两个应该是它们，那两个碱基进行了冗余读取，相当于第一次真正读取了64个碱基，从第二次开始都是只读取了62个新的碱基，所以是按照62个碱基的大小进行划分，l-2是第一次读取64个碱基，取上整是因为l-2并非恰好是62的倍数，如果最后剩余的碱基不足62个，也需要新的一轮来完成。[0050]进一步的，上文中给出的采用avx512向量寄存器是为了更好的进行解释说明，可以理解的，在具体实施过程中还可以采用其他存储容量的向量寄存器。[0051]进一步的，为了证明本发明所述方案（如图4所示为其整体框架结构示意图）的有效性，以下进行了实验测试：其中，图4中的low-quality-trim表示低质量碱基裁剪工作。根据fastq格式数据中的质量分数，使用一个固定大小的滑动窗口，从测序数据的3’端到5’端进行扫描，当滑动窗口内的碱基平均质量分数小于用户指定的阈值，将滑动窗口左移一个碱基并继续判断，直到滑动窗口内的碱基平均质量分数满足阈值，滑动窗口停止滑动，然后将滑动窗口所在位置及其后面（靠近3’端）的碱基序列裁剪掉。[0052]图4中的adapter-trim表示裁剪测序序列中的adapter即测序中使用的接头，adapter本质也是一段核苷酸序列。先根据用户指定的adapter序列的前三个碱基，将其作为“种子”，在测序序列中找到“种子”出现的所有位置，将这些位置作为候选位置。然后从左到右遍历序列中的候选位置，将候选位置处的序列与用户指定的adapter序列进行比较，计算海明编辑距离，距离满足用户指定阈值的第一个候选位置作为裁剪位置，将该候选位置及其后面（靠近3’端）的碱基序列裁剪掉。这个过程使用了向量化的方式，加快了寻找候选位置的速度，提高工作线程的处理效率，向量化过程采用上面的算法1。[0053]为了评估rabbittrim程序的正确性和性能，我们在不同的模拟数据和真实数据上进行了一系列的实验和分析，具体实验如下：（1）正确性测试本实施例中使用模拟数据进行正确性验证。我们使用模拟数据生成软件生成了在不同的测序错误率下的单端数据和双端数据，经过实验测试，我们的rabbittrim在正确性上和原始的ktrim保持一致。[0054]（2）性能测试本实施例中使用了多个模拟数据和真实数据进行了性能测试的实验，记录了rabbittrim、ktirm以及trimmomatic的运行时间和线程扩展性。下面以某个真实数据为例，数据来自ncbi(nationalcenterforbiotechnologyinformation)数据库的geo（gse81178），使用gsm2144218~gsm2144224并且将文件名后缀为read1.fastq的合并为pooled.read1.fq，将文件后缀名为read2.fastq合并为pooled.read2.fq，每个文件的大小为19gb。[0055]表1单端数据运行时间（秒）：线程数工具名称12468rabbittrim42.6821.411.149.529.58ktrim51.8638.7626.5826.8626.69trimmomatic167.2386.5287.5184.2586.64表2单端数据加速比：线程数工具名称12468rabbittrim11.993.834.494.46ktrim11.341.951.931.94trimmomatic11.931.911.981.93表3双端数据运行时间（秒）：线程数工具名称12468rabbittrim62.8832.7722.6922.7322.51ktrim96.0468.5669.0566.8169.32trimmomatic364.38149.72124.76124.6123.19表4双端数据加速比：线程数工具名称12468rabbittrim11.922.772.772.79ktrim11.41.391.441.39trimmomatic12.432.922.922.96通过表1至表4的实验结果我们可以看出，rabbittrim在性能上相比ktrim有了较大的提升，单线程执行速度提升，线程扩展性更好，最佳处理速度是ktrim的3倍左右，约为2gb/秒（在相同平台下ktrim约为730mb/秒），达到了测试平台的磁盘读写性能峰值（iobound）。[0056]本发明所述方案通过对输入输出模块以及向量化模块的优化，解决了io效率低的问题，提高了寻找adapter的速率。经过实验测试，rabbittrim在正确性上和ktrim保持一致，但是性能得到了较大的提升，数据处理的速度可接近磁盘读写的性能峰值（~2gb/秒）。在相同的实验平台上，rabbittrim处理数据的最佳性能峰值是原始ktrim的3倍左右，较大幅度地提升了对新一代测序数据预处理的速度，对于加快下游分析任务有重要意义，为生物信息领域的工作者提供了一种更加快速的quality-adapter-trim工具。[0057]实施例二：本实施例的目的是提供一种生物测序序列快速修剪系统。[0058]一种生物测序序列快速修剪系统，包括：数据获取单元，其用于获取待修剪的生物测序序列；数据处理单元，其用于对所述生物测序序列进行读操作、修剪操作以及写操作；其中，基于生产者—消费者模型对所述读操作、修剪操作以及写操作进行解耦，实现异步执行；且所述生物测序序列的格式化过程从读操作中转移到修剪操作中。[0059]在更多实施例中，还提供：一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成实施例一中所述的方法。为了简洁，在此不再赘述。[0060]应理解，本实施例中，处理器可以是中央处理单元cpu，处理器还可以是其他通用处理器、数字信号处理器dsp、专用集成电路asic，现成可编程门阵列fpga或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。[0061]存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据、存储器的一部分还可以包括非易失性随机存储器。例如，存储器还可以存储设备类型的信息。[0062]一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成实施例一中所述的方法。[0063]实施例一中的方法可以直接体现为硬件处理器执行完成，或者用处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器、闪存、只读存储器、可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成上述方法的步骤。为避免重复，这里不再详细描述。[0064]本领域普通技术人员可以意识到，结合本实施例描述的各示例的单元即算法步骤，能够以电子硬件或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。[0065]上述实施例提供的一种生物测序序列快速修剪方法及系统可以实现，具有广阔的应用前景。[0066]以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。当前第1页12当前第1页12