CN112397142A

CN112397142A - 面向多核处理器的基因变异检测方法及系统

Info

Publication number: CN112397142A
Application number: CN202011090874.2A
Authority: CN
Inventors: 刘卫国; 张�浩; 宋红磊; 殷泽坤
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2020-10-13
Filing date: 2020-10-13
Publication date: 2021-02-23
Anticipated expiration: 2040-10-13
Also published as: CN112397142B

Abstract

本发明公开了面向多核处理器的基因变异检测方法，包括：对输入数据进行预处理；从预处理后得到的测序序列read提取简要比对信息表达式CIGAR信息；输入数据，是指：将待查询序列与参考序列进行比对得到的文件；对读取的测序序列read的简要比对信息表达式CIGAR信息进行修改；对修改后的测序序列read的简要比对信息表达式CIGAR信息进行处理，处理过程中从内存池中进行候选变异数据的调取，得到候选变异集合；对候选变异集合中变异的基因进行局部重比对，以降低假阳性变异的检测；对局部重比对后的变异的基因进行格式化处理，将格式化处理后的变异基因输出到输出文件中，并且将内存池中数据进行重置以便进行反复的使用。

Description

面向多核处理器的基因变异检测方法及系统

技术领域

本申请涉及基因变异检测技术领域，特别是涉及面向多核处理器的基因变异检测方法及系统。

背景技术

本部分的陈述仅仅是提到了与本申请相关的背景技术，并不必然构成现有技术。

随着测序技术的发展，大量的测序数据爆炸式生长，对变异检测的计算能力提出了挑战。目前有非常多的变异检测工具被提出，像基于haplotype的GATKHaplotypeCaller、Mutect2、FreeBayes，基于启发式的方法VarDict、Lofreq和VarScan2；基于深度学习的DeepVariant等。

发明人发现受限于java的虚拟机机制和java虚拟机的额外开销，许多基于Java的变异检测工具，例如VarDict，在多和平台上的性能表现不佳，不能充分发挥现代多核处理器的性能优势，不能够快速精准的实现基因变异的检测。

发明内容

为了解决现有技术的不足，本申请提供了面向多核处理器的基因变异检测方法及系统；本申请能快速精准的实现基因变异的检测。

第一方面，本申请提供了面向多核处理器的基因变异检测方法；

面向多核处理器的基因变异检测方法，包括：

对输入数据进行预处理；从预处理后得到的测序序列read提取简要比对信息表达式CIGAR信息；所述输入数据，是指：将待查询序列与参考序列进行比对得到的文件；

对读取的测序序列read的简要比对信息表达式CIGAR信息进行修改；对修改后的测序序列read的简要比对信息表达式CIGAR信息进行处理，处理过程中从内存池中进行候选变异数据的调取，得到候选变异集合；

对候选变异集合中变异的基因进行局部重比对，以降低假阳性变异的检测；

对局部重比对后的变异的基因进行格式化处理，将格式化处理后的变异基因输出到输出文件中，并且将内存池中数据进行重置以便进行反复的使用。

第二方面，本申请提供了面向多核处理器的基因变异检测系统；

面向多核处理器的基因变异检测系统，包括：

预处理模块，其被配置为：对输入数据进行预处理；从预处理后得到的测序序列read提取简要比对信息表达式CIGAR信息；所述输入数据，是指：将待查询序列与参考序列进行比对得到的文件；

候选变异集合生成模块，其被配置为：对读取的测序序列read的简要比对信息表达式CIGAR信息进行修改；对修改后的测序序列read的简要比对信息表达式CIGAR信息进行处理，处理过程中从内存池中进行候选变异数据的调取，得到候选变异集合；

局部重比对模块，其被配置为：对候选变异集合中变异的基因进行局部重比对，以降低假阳性变异的检测；

输出模块，其被配置为：对局部重比对后的变异的基因进行格式化处理，将格式化处理后的变异基因输出到输出文件中，并且将内存池中数据进行重置以便进行反复的使用。

第三方面，本申请还提供了一种电子设备，包括：一个或多个处理器、一个或多个存储器、以及一个或多个计算机程序；其中，处理器与存储器连接，上述一个或多个计算机程序被存储在存储器中，当电子设备运行时，该处理器执行该存储器存储的一个或多个计算机程序，以使电子设备执行上述第一方面所述的方法。

第四方面，本申请还提供了一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成第一方面所述的方法。

第五方面，本申请还提供了一种计算机程序(产品)，包括计算机程序，所述计算机程序当在一个或多个处理器上运行的时候用于实现前述第一方面任意一项的方法。

与现有技术相比，本申请的有益效果是：

本申请通过使用多线程处理，线程拓展性好；使用内存池策略，可以减少频繁的内存分配带来的额外时间消耗，能够快速进行基因变异检测；

本申请通过采用变异检测工具Vardict的算法流程，并且与vardict的检测结果保持很高的一致性，能够精准的实现变异检测。

附图说明

构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1为第一个实施例的方法流程图；

图2(a)和图2(b)为第一个实施例的OpenMP动态任务调度；

图3为第一个实施例的germline mode下本申请和VarDict的运行时间图；

图4为第一个实施例的germline mode下本申请和VarDict的加速比示意图；

图5为第一个实施例的somatic mode下本申请和VarDict的运行时间图；

图6为第一个实施例的somatic mode下本申请和VarDict的加速比示意图。

具体实施方式

应该指出，以下详细说明都是示例性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

术语解释：

1.基因组：每个人都有一个基因组，这里的“基因组”并不只是“基因”的集合，基因是控制性状的遗传单元(性状也可以狭义的理解为个体的各种外在和内在特征，比如头发和眼睛颜色，高矮胖瘦，抵抗力强等)，但是基因组所指的其实是我们的所有遗传信息，而不单单只是一些外在和内在特征，也包含很多目前而言不明其功能性(或者被认为无功能)的DNA序列。其实说白了就是整一个的DNA序列！因而，基因也只是基因组的一个子集。此外，需要特别指出的是，我们虽都为“人”，但人与人之间的基因组是不一样的(即是多态的)，彼此之间都存在着一些差异，即使是和父母或是兄弟姐妹之间去比较。这些差异也是造成我们彼此之间为何如此这般不同的一个重要原因。而这些差异也是基因组多态性的来源。

2.Reads：这里的reads是一个基因组测序中的名词，指的就是一段特定长度的DNA片段，这个长度取决于测序仪的读长。

3.Germline和Somatic：根据发生变异的位置，变异可以分为Germline(种系变异)和Somatic(体细胞变异)变异。种系变异是发生在生殖细胞中的变异，体细胞变异顾名思义是发生在体细胞中的变异，通常癌细胞检验属于somatic变异，通常需要两个测序文件：正常细胞(normal)和癌细胞(tumor)的测序文件；而检测germline变异通常只需要一个测序文件即可。

3.变异是一个相对的概念，只有在彼此的比较中才有存在的意义。目前关于人类基因组变异的讨论，都是以“人类基因组计划”中所组装出来的人类基因组作为参照物。以下谈到的涉及比对过程所用的基因组指的就是这个人类参考基因组。

4.CIGAR：(Compact Idiosyncratic Gapped Alignment Report),简要比对信息表达式，其以参考序列为基础，使用数字加字母表示比对结果，比对结果信息,匹配碱基数，可变剪接等等。

5.SoftClip是指虽然比对不到基因组，但是还是存在于SEQ(segment SEQuence)中的序列，此时CIGAR列对应的S(Soft)的符号。直白点说，就是虽然比对不上参考基因组，但是在BAM/SAM文件中的reads上还是存在的序列(并没有被截断扔掉的序列)。只要一条序列两端有比对不上的序列部分，就是Soft Clip，这个一条序列上有比对不上的部分的现象是必然存在的(比如结构变异的断点的部分)，这种两端比对不上的read的特殊的表示方法，就是Soft Clip。Soft Clip是可以独立存在的。

常见的基因组变异一般可以归为如下几类：

SNP，单核苷酸多态性,一个碱基的变化

INDEL，插入或缺失,一个碱基的增加或移除

SNV，单核苷酸变异，一个碱基的改变，可以是SNP，也可以是INDEL

MNP，多核苷酸多态性，一个区块中有多个SNP

MNV，多核苷酸变异，一个区块中有多个SNP或INDEL

SV,结构变异，通常是上千个碱基，甚至是染色体级别上的变异

实施例一

本实施例提供了面向多核处理器的基因变异检测方法；

面向多核处理器的基因变异检测方法，包括：

S101：对输入数据进行预处理；从预处理后得到的测序序列read提取简要比对信息表达式CIGAR信息；所述输入数据，是指：将待查询序列与参考序列进行比对得到的文件；

S102：对读取的测序序列read的简要比对信息表达式CIGAR信息进行修改；对修改后的测序序列read的简要比对信息表达式CIGAR信息进行处理，处理过程中从内存池中进行候选变异数据的调取，得到候选变异集合；

S103：对候选变异集合中变异的基因进行局部重比对，以降低假阳性变异的检测；

S104：对局部重比对后的变异的基因进行格式化处理，将格式化处理后的变异基因输出到输出文件中，并且将内存池中数据进行重置以便进行反复的使用。

作为一个或多个实施例，所述S101：对输入数据进行预处理，具体步骤包括：

S1011：利用多线程获取用户指定的待进行变异检测的位点region信息并多线程并行地进行后续处理；

S1012：根据用户指定的待进行变异检测的位点region信息，对参考基因序列进行存储；

S1013：读取测序序列的read，过滤掉低质量的测序序列read；得到处理后的测序序列read；其中，测序序列read表示设定长度的DNA片段。

应理解的，Region信息，就是指待进行变异检测的位点，比如说“ch1:10000-20000”就是指要检测染色体一号的第10000个位点到第20000位点的所有变异。

示例性的，所述S1012根据用户指定的region信息，对参考基因序列进行处理；具体步骤包括：

比如用户指定的region信息是“chr：10000-20000”我们需要这部分位点的参考序列的信息，比如是“ACCTC...”，这些序列信息会被存到一个map数据结构里，比如<10000,‘A’>,<10001,‘C’>。

进一步地，S1013：read的质量的判断步骤包括：

S10131：通过测序过程，产生测序质量；

S10132：通过比对过程，产生比对质量和未比对上的碱基的个数；

S10133：根据测序质量、比对质量和未比对上的碱基的个数，计算出测序序列read的质量值，将测序序列read的质量值与设定阈值进行比较，高于设定阈值，则表示当前测序序列read为高质量read；否则，表示当前read为低质量read。

示例性的，所述测序过程，是指：对一个生物样本进行DNA序列的测序，是由illumia测序仪进行测序的。illumia测序仪输出测出的序列以及序列中每个位点的质量(可以认为质量越高可信度越大)。

示例性的，所述比对过程，是指：测序数据产生的是一段段的小序列(比如第二代测序数据大概是100bp～500bp左右)。比对过程就是通过BWA-MEM或者bowtie2，找到每一段小序列对应的参考序列的具体位置。

进一步地，测序序列read的质量值，获取步骤包括：

1.测序序列read的质量值mapping quality，从bam文件中获得，如果低于用户设置的值(未设置的话默认为0)则不是高质量read；

2.未比对上的碱基数也可以由bam文件中获得(samtools提供一个c语言的api，里面可以读取bam文件中record的各种信息，就是利用这个api获取信息进行过滤操作和其他处理的)，获得未比对上的碱基数如果大于用户设置的最大未比对上碱基(mismatch)值，则不是高质量的read；

3.碱基的测序质量一般不会直接用在过滤read，而是在进行变异检测过程中，如果检测到的变异的碱基质量值过低就不认为它是一个可靠的变异(比如检测到又一个位点的碱基由A变成了C，但是这个C碱基的质量小于15，那么就不认为它是一个变异，而是当作可能测错了)。

作为一个或多个实施例，所述S102：对读取的测序序列read的简要比对信息表达式CIGAR信息进行修改；具体步骤包括：

S1021：将查询序列与参考序列相比的改变值与改变值近邻的SoftClip一起替换为一个SoftClip；所述改变值，包括：插入值或删除值；

S1022：如果查询序列与参考序列相比的改变值与SoftClip相隔的匹配碱基数不超过10个，则查询序列与参考序列相比的改变值、碱基数与SoftClip也会被处理为一个SoftClip；

S1023：将长度小于设定阈值的匹配序列与开始或结尾处的插入处理成SoftClip；

S1024：将两个相邻近的删除合并为一个删除；将两个相近邻的插入合并为一个插入；

S1025：将三个相邻近的删除合三为一；

S1026：重新比对下SoftClip部分，把和参考序列匹配的部分更改为匹配的序列。

示例性的，所述S1021中，例如：12S2I78M会被处理成14S78M，12S2D78M被处理成12S78M。

示例性的，所述S1022中，例如12S5M2I78M会被处理成19S78M，12S5M2D78M会被处理成17S78M。

示例性的，所述S1023中，2M8I78M会被处理成10S78M；8I78M会被处理成8S78M。

示例性的，所述S1024中，相邻近的意思是中间隔着<＝2的匹配的碱基数，2这个值可以用户自己设置例如50M8I2M2I20M合并为50M12I50M。

示例性的，所述S1026中，12S78M代表的这个read，如果我们发现12S中靠近78M的部分由2个碱基是和参考基因相匹配的，那么这个cigar信息就会更改为10S80M。

应理解的，CIGAR＝56M1I30M；它的含义就是：这条序列与参考基因组相比；前56bp能够match上；中间有1bp的insertion(相比于参考基因组有1bp的插入)；最后是30bp的match；CIGAR＝145M，含义就是：这条序列与参考基因组比对的结果是145bp完全match上。

作为一个或多个实施例，所述S102：对修改后的测序序列read的简要比对信息表达式CIGAR信息进行处理，处理过程中从内存池中进行候选变异数据的调取，得到候选变异集合具体步骤包括：

S102a1：对修改后的测序序列read的简要比对信息表达式CIGAR信息，找到不匹配的位点，从内存池中进行数据的调取用以存储候选变异，并加入候选变异集合，如果是有多个连续的不匹配，那这个变异可能是多核苷酸变异MNV，把这多个连续的不匹配序列加入到候选变异集合里面；

S102a2：对于SoftClip的处理，使用SoftClip类来存储与SoftClip有关的变异信息，遍历这块序列，将不匹配的位点和此位点的相关信息存到SoftClip实例里；

S102a3：对Insertion插入和deletion删除的处理，将其加入候选变异集合里。

示例性的，12S67M2I43M，对每个CIGAR部分(12S、67M、2I、43M)分别进行处理，对于处理M(M包括match和mismatch，也就是匹配和不匹配)遍历这块序列，

作为一个或多个实施例，所述S103：对候选变异集合中变异的基因进行局部重比对，以降低假阳性变异的检测；具体步骤包括：

S1031：对read末段的Indels进行重新检验；如果发现有因测序过程被切断的变异，则将被切断的变异合并为同一个变异；

Indels指的是在基因组的某个位置上所发生的小片段序列的插入或者删除，其长度通常在50bp以下；

S1032：对碱基的插入、删除和多核苷酸变异MNV进行重新检验；如果发现有因测序过程被切断的变异，则将被切断的变异合并为同一个变异。

示例性的，所述S1031：对read末段的Indels进行重新检验；具体步骤包括：

对read末端的indels，我们在之前的步骤已经转化为softclip，对于每个位点上的softclip，我们用一个数据结构来存，里面包括了从这个位点开始前面的(5’端)或者后面的(3‘端)几个位点(取决于softclip的长度)的碱基情况：有几个A、G、C、T，然后通过类似“投票”的机制，如果某个位点上有3个A，1个T，那么这个位点碱基就定为A。

通过上面这种方式得到该位点上的一个序列，然后这个序列与其他已知的该位点上的其他变异比较，如果有的话就将这两个变异合并。

示例性的，所述S1032：有些MNV在测序的过程中被切断，比如说有个MNV是“ATAC”，如果这个变异被切断的话，会生成其他的变异，比如“A”、“AT”，“ATA”，“TAC”等等，对于这个样一个MNV会便利上述情况并且合并。因为本质上这些都是同一个变异。

作为一个或多个实施例，所述S104：对局部重比对后的变异的基因进行格式化处理，将格式化处理后的变异基因输出到输出文件中；具体步骤包括：

S1041：对变异的基因通过变异频率进行过滤；

S1042：对变异的基因通过位置信息进行过滤；

S1043：对变异的基因通过比对质量进行过滤。

示例性的，S1041：对变异的基因通过变异频率进行过滤；具体步骤包括：

变异的频率(frequency)是由varsCount/positionCoverage得到的，二者都是在在处理变异信息的时候进行计数。

变异频率的过滤条件是frequency<config.freq，config.freq可以由用户指定，不指定的话默认值是0.01.

示例性的，S1042：对变异的基因通过位置信息进行过滤；具体步骤包括：

根据位置信息过滤的过滤条件是：

meanPosition<config.readPosFilter。

同config.freq，config.readPosFilter也是可以由用户指定，未指定的话默认为5。

config.readPosFilter：变异距离测序序列(read)两端最小距离的平均值，小于该值的过滤掉，留下大于等于该值的。

meanPositon记录的是该变异距离测序序列(read)两端最小距离的平均值。

计算方法是在处理变异信息的过程记录该变异距离包含该变异的record的两端最小值的和，最后再除以包含这个变异的record的个数(varsCount)。

示例性的，S1043：对变异的基因通过比对质量进行过滤；具体步骤包括：

根据位置信息过滤的过滤条件是：

meanMappingQuality<config.mapq，

meanMappingQuality：变异的平均比对质量，小于config.mapq的要被过滤掉，留下大于等于该值的。

默认值是0(即如果用户不设置的话不会使用这个条件来过滤)。

近年来，随着芯片技术和第二代高通量测序技术的发展，人类基因组上的结构性变异图谱才被真正全面而又集中地进行了研究。生物信息研究人员已针对这两种不同的技术开发了许多相对应的软件用于检测基因组的变异信息。相比较而言，虽然成本较高，但是基于测序的方法要明显优于芯片的检测，其中最重要的一个方面是，高通量测序技术能够在单碱基精度之下对全基因组范围内所有类型的变异进行检测，而芯片技术实际上只对大片段的序列删除比较敏感。

为了实现在多核平台的性能优势，本申请采用了和经典的变异检测工具VarDict相同的流程，主要分为如下四个模块：

预处理模块，目的是根据用户指定的region信息对参考序列的序列信息进行处理，同时过滤掉输入比对文件中低质量的read。read的质量通常由测序质量(由测序过程产生)、比对质量(由比对过程中产生)和比对过程中未比对上上的碱基的个数等共同计算得来。除此之外，预处理模块还会进行比对文件的读取以及修改，用于下一个模块的分析。预处理模块还支持使用用户指定的过滤规则来过滤掉一些read。

定位候选变异。对读取的read的CIGAR信息进行修改，以便进行后续的变异检测处理；然后，结合修改后的CIGAR信息进行处理，得到候选变异集合。在此步骤中，直接对record的CIGAR数组进行处理而不是现转化成字符串然后对字符串进行处理。

这里的直接处理是相对于将cigar信息转化为string再处理，fastvc对一些操作直接处理cigar数组(record的cigar信息是以uint32_t*的形式存储的，32位int中，后4位表示cigar类型，前28位表示cigar的长度，比如64M用uint32表示就是0x400)。主要操作有：

1.如果cigar是由”D”(deletion)开始或结尾的话，删除这部分即3D92M2D->92M；

2.如果cigar是由“I”开始或结尾，那么把这个部分由“I”转为“S”即2I92M3I->2S92M3S。

局部重比对。目的在于降低假阳性变异的检测。由于比对过程中的gap open的惩罚要比错误匹配的惩罚大，导致了对read末段的Indels的检验比较困难。同时一些Indels和MNV可能会在测序过程中被切断，所以，局部的重新比对可以有效的减少假阳性的出现。

后处理模块。是对检测出来的变异通过变异频率、位置信息、比对质量等进行过滤，然后格式化输出到输出文件中。

针对现代多核处理器做了很多优化。主要包括：

1.使用c++实现算法以便很进行更细粒度的优化。

2.手动设计使用内存池，优化内存使用，降低内存频繁分配的开销。

3.基于region的多线程机制。

4.针对程序中热点操作进行优化。

采用了和经典的变异检测工具VarDict相同的流程，VarDict由Java实现(最初版本有perl实现，java版本比perl版本号称快10倍)，囿于java的虚拟机机制，VarDict的线程拓展性并不是很好，导致VarDict并不能充分发挥现代多核处理的性能优势。我们使用c++实现，提升程序单线程性能的同时也方便我们进行更细粒度的优化。

使用手动设计的内存池来降低频繁的内存分配开销，整体设计流程如图1所示：内存池多线程开启时被初始化，每个线程维护自己的内存池，在处理过程中从用户指定的区域文件或者命令行中获取region信息，针对该region，进行CIGAR信息处理获取候选变异、局部重对齐、格式化输出等操作，其中，在生成候选变异过程中，从内存池中获取数据，对从内存池中获取的数据进行操作(统计候选变异信息，例如变异出现的次数，质量统计等等。)然后处理完当前region之后对内存池进行reset，然后处理下一个region，最后当所有region都处理完成之后再对内存池进行销毁。

变异检测属于易并行问题，因为不同的region之间通常是没有相互关联的。

使用基于region的多线程机制，使用OpenMP来实现多线程任务。以全基因组检测为例，变异检测的不同region写在区域文件(.bed)中，这些region通常是由用户指定的，用户可以对基因组的region进行任意划分。

使用多线程并行的处理这些region，同时用户也可以直接在命令行参数中指定一个region(使用-R参数)，此时只会使用单线程来处理该region。

多线程处理region的过程中，如果region的大小不一或者分配不均匀可能会导致多线程负载不均衡的问题，所以我们使用动态调度策略来解决负载不均衡的问题，如图2(a)和图2(b)所示，动态的任务调度可以有效减少程序的总运行时间。

动态调度就是指在多线程过程中，多个任务划分给多个线程的策略：比如有0-99号共100个任务，有0-9号共10个线程，如果是静态调度的话0-9号任务给0号线程，10-19号任务给1号线程....但是如果是动态调度的话，每个线程处理的任务是不一定的，而是一个线程执行完一个任务后系统再分配下一个任务。上述动态分配的策略就不会出现某个线程要处理的任务比其他的重很多的情况，每个线程处理的任务都是相当的。我们使用OpenMP来实现多线程，OpenMP的制导语句可以直接指定任务划分方式，比如#omp parallel xxschedule(dynamic)即可。

针对一些耗时的操作我们还使用更高效的第三方库以提升程序性能：

1).程序运行过程中需要频繁操作hash map，我们使用robin_hood_hashing，一种更高效、内存利用率更高的map实现，robin_hood_hashing还针对string类型做了针对性的优化。

2).htslib，一种处理测序数据(bam、sam、cram、vcf等文件格式)的c语言库。

耗时的操作具体有：

1.对map数据结构内部元素的频繁访问，比如说我们在发现一个潜在变异的时候会索引一下map表该位点上是否已经有相关的变异候选数据了，如果有的话返回这个候选变异的指针，然后进行一些操作；如果没有的的话就要新建一个。因为变异检测的数据量大，这个map会经常被访问到，所以有一个高效的map会降低程序运行时间

2.测序数据和比对数据通常比较大(上百G)能够有一种高效读取这些数据的工具可以有效降低程序运行时间。

我们测试了本申请与其他的变异检测工具的性能对比，结果如下表1所示：

表1本申请与其他目前主流变异检测工具的性能对比

可以看到，本申请无论在Germline mode还是在Somatic mode上的变异检测时间都是最短的，特别的，在germline mode下，CHM1_CHM13数据全基因组变异检测本申请只需要17分钟，而目前最常用的GATK需要816分钟，于此同时，在准确性方面，本申请也取得了与其他工具相当甚至更高的F1-Score。由于本申请和VarDict采用相同的pipeline，我们针对本申请和VarDict的程序运行时间和线程拓展性做了对比，结果显示本申请在somatic和germline mode下都具有较短的运行时间和较好的线程加速比。

于此同时，我们比较了本申请和VarDict对变异的检测结果，如表2所示。本申请和VarDict的检测结果在准确性、召回率和F1-score方面都具有很高的一致性。

表2本申请和VarDict的检测结果对比

图3为第一个实施例的germline mode下本申请和VarDict的运行时间图；图4为第一个实施例的germline mode下本申请和VarDict的加速比示意图；图5为第一个实施例的somatic mode下本申请和VarDict的运行时间图；图6为第一个实施例的somatic mode下本申请和VarDict的加速比示意图。

实施例二

本实施例提供了面向多核处理器的基因变异检测系统；

面向多核处理器的基因变异检测系统，包括：

此处需要说明的是，上述预处理模块、候选变异集合生成模块、局部重比对模块和输出模块对应于实施例一中的步骤S101至S104，上述模块与对应的步骤所实现的示例和应用场景相同，但不限于上述实施例一所公开的内容。需要说明的是，上述模块作为系统的一部分可以在诸如一组计算机可执行指令的计算机系统中执行。

上述实施例中对各个实施例的描述各有侧重，某个实施例中没有详述的部分可以参见其他实施例的相关描述。

所提出的系统，可以通过其他的方式实现。例如，以上所描述的系统实施例仅仅是示意性的，例如上述模块的划分，仅仅为一种逻辑功能划分，实际实现时，可以有另外的划分方式，例如多个模块可以结合或者可以集成到另外一个系统，或一些特征可以忽略，或不执行。

实施例三

本实施例还提供了一种电子设备，包括：一个或多个处理器、一个或多个存储器、以及一个或多个计算机程序；其中，处理器与存储器连接，上述一个或多个计算机程序被存储在存储器中，当电子设备运行时，该处理器执行该存储器存储的一个或多个计算机程序，以使电子设备执行上述实施例一所述的方法。

应理解，本实施例中，处理器可以是中央处理单元CPU，处理器还可以是其他通用处理器、数字信号处理器DSP、专用集成电路ASIC，现成可编程门阵列FPGA或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。

存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据、存储器的一部分还可以包括非易失性随机存储器。例如，存储器还可以存储设备类型的信息。

在实现过程中，上述方法的各步骤可以通过处理器中的硬件的集成逻辑电路或者软件形式的指令完成。

实施例一中的方法可以直接体现为硬件处理器执行完成，或者用处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器、闪存、只读存储器、可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成上述方法的步骤。为避免重复，这里不再详细描述。

本领域普通技术人员可以意识到，结合本实施例描述的各示例的单元即算法步骤，能够以电子硬件或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。

实施例四

本实施例还提供了一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成实施例一所述的方法。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims

1.面向多核处理器的基因变异检测方法，其特征是，包括：

2.如权利要求1所述的方法，其特征是，对输入数据进行预处理，具体步骤包括：

利用多线程获取用户指定的待进行变异检测的位点region信息并多线程并行地进行后续处理；

根据用户指定的待进行变异检测的位点region信息，对参考基因序列进行存储；

读取测序序列的read，过滤掉低质量的测序序列read；得到处理后的测序序列read；其中，测序序列read表示设定长度的DNA片段。

3.如权利要求2所述的方法，其特征是，read的质量的判断步骤包括：

通过测序过程，产生测序质量；

通过比对过程，产生比对质量和未比对上的碱基的个数；

根据测序质量、比对质量和未比对上的碱基的个数，计算出测序序列read的质量值，将测序序列read的质量值与设定阈值进行比较，高于设定阈值，则表示当前测序序列read为高质量read；否则，表示当前read为低质量read。

4.如权利要求1所述的方法，其特征是，对读取的测序序列read的简要比对信息表达式CIGAR信息进行修改；具体步骤包括：

将查询序列与参考序列相比的改变值与改变值近邻的SoftClip一起替换为一个SoftClip；所述改变值，包括：插入值或删除值；

如果查询序列与参考序列相比的改变值与SoftClip相隔的匹配碱基数不超过10个，则查询序列与参考序列相比的改变值、碱基数与SoftClip也会被处理为一个SoftClip；

将长度小于设定阈值的匹配序列与开始或结尾处的插入处理成SoftClip；

将两个相邻近的删除合并为一个删除；将两个相近邻的插入合并为一个插入；

将三个相邻近的删除合三为一；

重新比对下SoftClip部分，把和参考序列匹配的部分更改为匹配的序列。

5.如权利要求1所述的方法，其特征是，对修改后的测序序列read的简要比对信息表达式CIGAR信息进行处理，处理过程中从内存池中进行候选变异数据的调取，得到候选变异集合具体步骤包括：

对修改后的测序序列read的简要比对信息表达式CIGAR信息，找到不匹配的位点，从内存池中进行数据的调取用以存储候选变异，并加入候选变异集合，如果是有多个连续的不匹配，那这个变异可能是多核苷酸变异MNV，把这多个连续的不匹配序列加入到候选变异集合里面；

对于SoftClip的处理，使用SoftClip类来存储与SoftClip有关的变异信息，遍历这块序列，将不匹配的位点和此位点的相关信息存到SoftClip实例里；

对Insertion插入和deletion删除的处理，将其加入候选变异集合里。

6.如权利要求1所述的方法，其特征是，对候选变异集合中变异的基因进行局部重比对，以降低假阳性变异的检测；具体步骤包括：

对read末段的Indels进行重新检验；如果发现有因测序过程被切断的变异，则将被切断的变异合并为同一个变异；Indels指的是在基因组的某个位置上所发生的小片段序列的插入或者删除，其长度通常在50bp以下；

对碱基的插入、删除和多核苷酸变异MNV进行重新检验；如果发现有因测序过程被切断的变异，则将被切断的变异合并为同一个变异。

7.如权利要求1所述的方法，其特征是，对局部重比对后的变异的基因进行格式化处理，将格式化处理后的变异基因输出到输出文件中；具体步骤包括：

对变异的基因通过变异频率进行过滤；对变异的基因通过位置信息进行过滤；对变异的基因通过比对质量进行过滤。

8.面向多核处理器的基因变异检测系统，其特征是，包括：

9.一种电子设备，其特征是，包括：一个或多个处理器、一个或多个存储器、以及一个或多个计算机程序；其中，处理器与存储器连接，上述一个或多个计算机程序被存储在存储器中，当电子设备运行时，该处理器执行该存储器存储的一个或多个计算机程序，以使电子设备执行上述权利要求1-7任一项所述的方法。

10.一种计算机可读存储介质，其特征是，用于存储计算机指令，所述计算机指令被处理器执行时，完成权利要求1-7任一项所述的方法。