CN108595918A - 循环肿瘤dna重复序列的处理方法及装置 - Google Patents

Abstract

本发明公开了一种循环肿瘤DNA重复序列的处理方法及装置。其中，该方法包括：获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，测序数据包括：多对双端序列；将测序数据和参考基因组序列进行比对，得到第一比对结果，其中，第一比对结果至少包括：多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列；基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。本发明解决了现有技术中测序数据的处理方法对样本测序进行重复序列删除或标记，准确度低的技术问题。

Description

循环肿瘤DNA重复序列的处理方法及装置

技术领域

本发明涉及遗传工程领域，具体而言，涉及一种循环肿瘤DNA重复序列的处理方法及装置。

背景技术

肿瘤细胞在进行分裂增殖过程当中，会凋亡、死亡、坏死，也会主动向体液中释放携带有肿瘤突变的DNA碎片，也即循环肿瘤DNA(circulating tumor DNA，简称为ctDNA)，多存在于血液、滑膜液和脑脊液等体液中，尤其是血浆游离DNA(cell-free DNA，cfDNA)中。通过对ctDNA的测序，检测肿瘤细胞DNA分子上发生碱基序列改变(突变)的基因组区域，能够有效反应病人对治疗的响应；在检测到药物响应之后，肿瘤有可能对药物治疗产生耐药，ctDNA检测也可以追踪耐药突变的产生，定性定量；检测手术后是否存在残余组织，判断预后效果以及早期肿瘤的筛查。不同于常规的基因组DNA，ctDNA片段较短，通常只有100～400bp，而且在血液中含量较少，所以实际中能提取到的ctDNA量含量很低。

由于ctDNA片段较短且含量较低的特点，因此在提取量较少时，需要在建库阶段进行多轮聚合酶链式反应(Polymerase Chain Reaction，简称为PCR)，扩大原始提取DNA的含量，以产生足够的DNA分子数目做高通量测序(High-throughput sequencing，简称为NGS测序)和后续生物信息学分析。由于PCR扩增导致对一个分子进行多次镜像复制，产生重复序列(Duplicated reads)，这些无效的重复数据对于检测变异极容易引入人工误差。对于最理想的NGS数据分析流程中，都需要尽可能把所有通过PCR获得的测序数据全部去除，还原到没有PCR的状态。

现有技术中，提供了两种重复序列的去重方法，samtools rmdup和Picard’sMarkDuplicates。其中，samtools rmdup的工作原理为：NGS测序得到的序列(read)通过与人类参考基因组比对(mapping)，得到这条read的比对位置，如果不同的reads比对到相同的基因组位置，则认为这部分的reads是通过PCR产生的多个重复序列，只保留mapping质量最高的read，删除其余的重复序列。对于PE reads，如果两端的read比多到基因组的不同染色体上或者两者之前的距离过长(即不是Proper Paired)，则不作去重考虑。Picard’sMarkDuplicates的基本思路与samtools rmdup相同，通过比较reads中5'端的mapping位置，对于具有相同5'位置的序列，选取测序质量最高的reads作为去重后保留的唯一reads，且对于PE reads不是Proper Paired的情况也会做去重处理。

由于血液中真正含有肿瘤信息的ctDNA含量很少，浓度很低，所以在血液中最终提取的DNA中大部分是并不包含肿瘤信息的基因组DNA(也即，gDNA)。这就导致最终检测到的肿瘤相关变异的突变频率很低，其中，突变频率＝覆盖到该位点的含有该变异的reads个数/覆盖到该位点上全部reads个数。由于PCR过程可能伴有随机错误的产生，这些错误最后很可能造成突变检出准确率下降。由于Picard和Samtools去重仅根据质量值保留最终unique reads，而这种策略必将导致下述三种情况中的任意一种发生：正确保留了含有真实突变的read作为最终的uniqe read；错误保留了含有假阳性突变的read；错误保留了未含有真实突变的假阴性的read，而samtools rmdup和Picard’s MarkDuplicates并不能有效控制后两种情况的发生，也即没有方法对PCR过程中产生的错误进行校正。通常为了减少假阳性位点的检出，而采用的策略是提高检出线(即最低可准确检出的突变频率)，但是，这样会导致很多真阳性的低频突变(即低于检出线的突变)位点不能有效检出。

针对现有技术中测序数据的处理方法对样本测序进行重复序列删除或标记，准确度低的问题，目前尚未提出有效的解决方案。

发明内容

本发明实施例提供了一种循环肿瘤DNA重复序列的处理方法及装置，以至少解决现有技术中测序数据的处理方法对样本测序进行重复序列删除或标记，准确度低的技术问题。

根据本发明实施例的一个方面，提供了一种循环肿瘤DNA重复序列的处理方法，包括：获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，测序数据包括：多对双端序列；将测序数据和参考基因组序列进行比对，得到第一比对结果，其中，第一比对结果至少包括：所述多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列；基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

进一步地，基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列包括：将多对双端序列划分为至少一个序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同；将每个序列集合中的每对双端序列逐碱基进行比较，得到每个序列集合中每个碱基位置上每个碱基类型的数量；获取每个序列集合中每个碱基位置上数量超过预设数量的第一碱基类型；将每个序列集合中每个碱基位置上的第一碱基类型进行组合，得到每个序列集合对应的一致性序列；根据每个一致性序列的每个碱基位置上第一碱基类型的碱基质量值及除第一碱基类型之外的其他碱基类型的碱基质量值，得到每个一致性序列对应的碱基质量值序列。

进一步地，在基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列之后，上述方法还包括：将至少一个一致性序列和参考基因组序列进行比对，得到第二比对结果；按照每个一致性序列的基因组位置，对第二比对结果进行排序，得到第三比对结果。

进一步地，在按照每个一致性序列的基因组位置，对第二比对结果进行排序，得到第三比对结果之后，上述方法还包括：根据第三比对结果，显示每个基因组位置对应的双端序列的比对信息和碱基质量值；对比对质量满足预设条件的双端序列进行过滤

进一步地，在按照每个一致性序列的基因组位置，对第二比对结果进行排序，得到第三比对结果之后，上述方法还包括：获取捕获测序区间；根据捕获测序区间，对每个一致性序列进行单核苷酸变异检测和插入缺失检测，得到检测结果。

进一步地，在基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列之前，上述方法还包括：按照多对双端序列的基因组位置，对第一比对结果进行排序，得到第四比对结果，并为第四比对结果建立索引；对第五比对结果进行过滤，得到第五比对结果；基于第五比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

进一步地，将测序数据和参考基因组序列进行比对，得到第一比对结果包括：获取多对双端序列中每条序列和参考基因组序列中的每段序列的匹配度；获取最高匹配度对应的至少一段序列，得到每条序列的匹配序列；根据每条序列的匹配序列，确定每条序列的基因组位置。

根据本发明实施例的另一方面，还提供了一种循环肿瘤DNA重复序列的处理装置，包括：获取模块，用于获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，测序数据包括：多对双端序列；比对模块，用于对测序数据和参考基因组序列进行比对，得到第一比对结果，其中，第一比对结果至少包括：多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列；处理模块，用于基于第一比对结果和每对双端序列中每个碱基位置上的碱基类型，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

根据本发明实施例的另一方面，还提供了一种存储介质，存储介质包括存储的程序，其中，在程序运行时控制存储介质所在设备执行上述的循环肿瘤DNA重复序列的处理方法。

根据本发明实施例的另一方面，还提供了一种处理器，处理器用于运行程序，其中，程序运行时执行上述的循环肿瘤DNA重复序列的处理方法。

在本发明实施例中，获取待检测循环肿瘤DNA的测序数据和参考基因组序列，将测序数据和参考基因组序列进行比对，得到第一比对结果，进一步基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列，从而实现重复序列的去重处理。容易注意到的是，由于在将测序数据和参考基因组序列进行比对之后，需要结合每对双端序列中每个碱基位置上的碱基类型，得到一致性序列，从而实现在考虑到序列质量值的同时，考虑到具体序列上的差异，达到减少人工错误，提高处理准确度的技术效果，进而解决了现有技术中测序数据的处理方法对样本测序进行重复序列删除或标记，准确度低的技术问题。

附图说明

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

图1是根据本发明实施例的一种循环肿瘤DNA重复序列的处理方法的流程图；

图2是根据本发明实施例的一种可选的循环肿瘤DNA重复序列的处理方法的示意图；

图3是根据本发明实施例的一种可选的循环肿瘤DNA重复序列的处理方法的流程图；以及

图4是根据本发明实施例的一种循环肿瘤DNA重复序列的处理装置的示意图。

具体实施方式

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

下面，首先对本发明实施例中出现的技术名词进行解释如下：

NGS测序：高通量测序(High-throughput sequencing)，又称“下一代”测序("Next-generation"sequencing，简称为NGS)，是相对于传统的桑格测序(SangerSequencing)而言的，以能一次并行对几十万到几百万条DNA分子进行序列测定和一般读长较短等特点为标志。NGS技术目前已经在诊断遗传病、测量基因表达水平、构建进化树、区分形态学上近似物种、对非模式生物进行重头测序头(de novo sequencing)，构建其参考基因组序列等多个领域得到广泛应用。NGS测序有单端(Single-end，简称为SE)和双端(Paird-end，PE)测序之分，目前大部分使用的都是PE测序技术，PE即指同一个DNA分子的两头，测序时先测一端，然后反过来再测另一端。测的这两条序列就叫做Paired ends Reads(PE Reads)，很多时候，它又被叫做"Mate Paires"。

捕获测序：将感兴趣的基因组区域定制成特异性探针与基因组DNA在序列捕获芯片(或溶液)进行杂交，将目标基因组区域的DNA片段进行富集后再利用NGS测序技术进行测序的研究策略。由于其捕获区域短，大大降低测序成本，被各领域广泛应用，肿瘤检测领域也多使用捕获测序的策略。

ctDNA：循环肿瘤DNA，肿瘤细胞在进行分裂增值过程当中，主动向体液中分泌的已经经历过基因突变的DNA片段。

PCR：聚合酶链式反应，一种用于放大扩增特定的DNA片段的技术。

重复序列：由于PCR扩增导致对一个分子进行多次镜像复制的后果。

reads：测序读长，测序仪测到的基因组或转录组序列片段。

fasta：一种基于文本用于表示核苷酸序列或氨基酸序列的格式。

fastq：一种常见的高通量测序文件类型，通常原始测序数据都是以该文件类型储存的。

bwa：一种比对方法软件，用于查找测序序列在人类基因参考序列中的位置,可输出bam格式结果文件。

sam：一种序列比对格式，用来存储测序序列回贴到参考基因组的结果

bam：sam文件的二进制压缩格式，用来存储测序序列回贴到参考基因组的结果。

samtools：一种处理bam/sam文件的工具。

picard：一种处理高通量测序数据的工具，可用于处理sam/bam等比对结果文件。

比对质量：用于量化比对到错误位置的可能性，值越高表示可能性越低。

CIGAR：简要比对信息表达式，其以参考序列为基础，使用数据加字母表示比对结果。

实施例1

根据本发明实施例，提供了一种循环肿瘤DNA重复序列的处理方法的实施例，需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

图1是根据本发明实施例的一种循环肿瘤DNA重复序列的处理方法的流程图，如图1所示，该方法包括如下步骤：

步骤S102，获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，测序数据包括：多对双端序列。

具体地，上述的待检测循环肿瘤DNA可以从病人的血液、淋巴液、组织间隙液、脑髓液等体液中提取得到的基因序列，在本发明实施例中以血液中提取到的ctDNA为例进行说明；上述的测序数据可以是对待检测ctDNA进行NGS测序得到的ctDNA样本捕获测序fastq数据；上述的参考基因组序列可以是从公开数据库NCBI等网站下载的人类参考基因组fasta数据。

步骤S104，将测序数据和参考基因组序列进行比对，得到第一比对结果，其中，第一比对结果至少包括：多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列。

具体地，上述的基因组位置可以是每对PE reads比对到参考基因组序列中的位置，不同的PE reads可以比对到相同的位置；上述的碱基质量值可以是通过NGS测序得到的测序质量，用于衡量每个碱基位置上碱基类型测量的准确度，碱基质量值越大，说明碱基类型测量的准确度越高；上述的碱基序列可以是每对双端序列中每个碱基位置上的碱基类型，DNA序列中包含四种类型的碱基，分别为G、C、T、A，在NGS测序过程中，可以确定每个碱基位置上的碱基类型，并得到该碱基类型的碱基质量值。

在一种可选的方案中，可以获取人类参考基因组fasta数据和ctDNA样本捕获测序fastq数据，利用基因组比对工具bwa mem进行序列比对，得到比对结果文件(.bam)，也即，得到上述的第一比对结果，比对结果文件为bam格式，包含每对PE reads的名称、位置信息、SAM标记、比对质量信息、CIGAR字串、mate pair信息、片段序列、测序质量等)。

需要说明的是，第一比对结果中的多对双端序列的碱基序列和碱基质量值是从NGS测序的测序数据中直接继承过来的数据，第一比对结果包含基因组比对位置及比对情况的信息的同时，还存储了多对双端序列的碱基序列和碱基质量值，方便后续的其他分析，不再使用fastq文件。

步骤S106，基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

在一种可选的方案中，由于PCR和测序错误是随机发生的小概率事件，所以其产生的碱基类型相比正确的碱基占比较少，因此，本发明实施例提供了一种MostP算法，基于上述基本事实，比对到基因组相同位置就可以认为来自同一个原始的fragment，且同一组PEreads内，占大多数的碱基类型构成新的一对reads(也即上述的一致性序列)，并将一致性序列作为最终代表的unique reads。

需要说明的是，无论是通过本发明提供的去重方法进行去重之前，还是通过本发明提供的去重方法进行去重之后，所有的bam文件均包含下列信息：每对PE reads的名称、位置信息、SAM标记、比对质量信息、CIGAR字串、mate pair信息、片段序列、测序质量等。

根据本发明上述实施例，获取待检测循环肿瘤DNA的测序数据和参考基因组序列，将测序数据和参考基因组序列进行比对，得到第一比对结果，进一步基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列，从而实现重复序列的去重处理。容易注意到的是，由于在将测序数据和参考基因组序列进行比对之后，需要结合每对双端序列中每个碱基位置上的碱基类型，得到一致性序列，从而实现在考虑到序列质量值的同时，考虑到具体序列上的差异，达到减少人工错误，提高处理准确度的技术效果，进而解决了现有技术中测序数据的处理方法对样本测序进行重复序列删除或标记，准确度低的技术问题。

可选地，在本发明上述实施例中，步骤S106，基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列包括：

步骤S1062，将多对双端序列划分为至少一个序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同。

具体地，上述的每个序列集合用于表征同一个fragment，每个序列集合中包含的PE reads为同一个fragment经过PCR产生的多个重复序列。

步骤S1064，将每个序列集合中的每对双端序列逐碱基进行比较，得到每个序列集合中每个碱基位置上每个碱基类型的数量。

步骤S1066，获取每个序列集合中每个碱基位置上数量超过预设数量的第一碱基类型。

具体地，上述的预设数量可以根据实际实验确定的数量，也可以是根据实际需要进行设定的数量。

步骤S1068，将每个序列集合中每个碱基位置上的第一碱基类型进行组合，得到每个序列集合对应的一致性序列。

在一种可选的方案中，可以将比对到相同基因组位置的全部PE reads中逐碱基的比较，挑选占绝大多数的碱基类型组成一对新的PE reads，也即得到每个序列集合对应的一致性序列。

例如，如图2所示，对于包含如2所示的8对PE reads的测序数据，假设8对PE reads的基因组位置相同，则可以将8对PE reads逐碱基进行比较，挑选每个碱基位置上数量最多的碱基类型，从而得到一致性序列，将该一致性序列作为独立序列。

步骤S1060，根据每个一致性序列的每个碱基位置上第一碱基类型的碱基质量值及除第一碱基类型之外的其他碱基类型的碱基质量值，得到每个一致性序列对应的碱基质量值序列。

在一种可选的方案中，对于每个一致性序列，可以计算每个碱基位置上第一碱基类型的所有碱基质量值的和，得到第一和值，同时计算每个碱基位置上除第一碱基类型之外的其他碱基类型的所有碱基质量值的和，得到第二和值，通过计算第一和值和第二和值的差值，得到每个碱基位置上的碱基质量值，从而得到每个一致性序列对应的碱基质量值序列。在此过程中，如果某个碱基位置上的碱基质量值>90，则确定该一致性序列在该碱基位置上的碱基质量值为90；，如果某个碱基位置上的碱基质量值<0，则确定该一致性序列在该碱基位置上的碱基质量值为0。

可选地，在本发明上述实施例中，在步骤S106，基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列之后，该方法还包括：

步骤S108，将至少一个一致性序列和参考基因组序列进行比对，得到第二比对结果。

在一种可选的方案中，由于一致性序列是重新生成的，不存在包含一致性序列的比对结果文件(.bam)，因此，在得到一致性序列之后，可以利用基因组比对工具bwa mem将一致性序列和人类参考基因组序列重新进行序列比对，得到重新比对结果文件(.bam)，也即，得到上述的第二比对文件。

步骤S110，按照每个一致性序列的基因组位置，对第二比对结果进行排序，得到第三比对结果。

在一种可选的方案中，为了使得相同位置的PE reads相邻，方便后续对PE reads进行去重处理，可以调用Picard’s SortSam模块将重新比对结果文件(.bam)(也即上述的第二比对结果)按比对位置进行排序。

可选地，在本发明上述实施例中，在步骤S110，按照每个一致性序列的基因组位置，对第二比对结果进行排序，得到第三比对结果之后，该方法还包括：

步骤S112，根据第三比对结果，显示每个基因组位置对应的双端序列的比对信息和碱基质量值。

具体地，上述的比对信息可以是是否与参考基因组的每个位置上的碱基相同，是否发生插入缺失，比对到参考基因组的正链还是负链等，本发明对此不作具体限定。

步骤S114，对比对质量满足预设条件的双端序列进行过滤。

具体地，上述的预设条件可以是比对质量为0。

在一种可选的方案中，可以调用samtools-1.3mpileup根据去重后的bam文件，按基因组位置展示比对到该位置的所有reads的比对信息和质量值，每一行记录一个基因组位置的比对情况，也即，每一行记录比对到每个基因组位置的所有reads的比对信息和质量值。可以认为是一种比对结果纵向的堆叠式的统计和展示。

可选地，在本发明上述实施例中，在步骤S110，按照每个一致性序列的基因组位置，对第二比对结果进行排序，得到第三比对结果之后，该方法还包括：

步骤S116，获取捕获测序区间。

步骤S118，根据捕获测序区间，对每个一致性序列进行单核苷酸变异检测和插入缺失检测，得到检测结果。

在一种可选的方案中，可以获取捕获测序区间Bed文件，并调用varscan2mpileup2snp模块检测单核苷酸变异(SNV)，mpileup2indel模块检测插入缺失(INDEL)，其中，单核苷酸变异是指参考基因组的某个位置上发生碱基类型的改变，插入缺失是指在参考基因组的某段序列上插入了一小段新的序列或缺失了某段序列。

可选地，在本发明上述实施例中，在步骤S106，基于第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列之前，该方法还包括：

步骤S120，按照多对双端序列的基因组位置，对第一比对结果进行排序，得到第四比对结果，并为第四比对结果建立索引。

在一种可选的方案中，可以调用Picard’s SortSam模块将比对结果文件(.bam)(也即上述的第一比对结果)按比对位置排序，同时建立bam文件的索引文件(.bai)。通过比对结果文件按比对位置排序，从而使得相同位置的PE reads相邻，方便后续对PE reads进行去重处理。

步骤S122，对第四比对结果进行过滤，得到第五比对结果。

在一种可选的方案中，由于同一个PE reads可能会对比到多个基因组位置，在进行去重处理之前，首先需要对比对结果文件(.bam)进行过滤，具体可以调用samtools view模块对排序后的bam文件进行筛选，得到第五比对结果。

步骤S124，基于第五比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

在一种可选的方案中，可以将第五比对结果中，比对到基因组相同位置就可以认为来自同一个原始的fragment，且同一组PE reads内，占大多数的碱基类型构成新的一对reads(也即上述的一致性序列)，并将一致性序列作为最终代表的unique reads。

可选地，在本发明上述实施例中，步骤S104，将测序数据和参考基因组序列进行比对，得到第一比对结果包括：

步骤S1042，获取多对双端序列中每条序列和参考基因组序列中的每段序列的匹配度。

步骤S1044，获取最高匹配度对应的至少一段序列，得到每条序列的匹配序列。

具体地，上述的预设相似度可以根据实际检测需求进行设定，本发明对此不做具体限定。

步骤S1046，根据每条序列的匹配序列，确定每条序列的基因组位置。

在一种可选的方案中，可以计算每一对PE reads中每条reads与人类参考基因组序列的匹配度，通过匹配度判断每一条reads是否来自人类参考基因组序列中某一段序列，匹配度越高，每一条reads是自人类参考基因组序列中该序列的可能性越大，可以将每条reads比对到最高匹配度的序列，从而根据该序列的位置，可以得到该条reads的基因组位置。

需要说明的是，在本发明实施例中，可以采用现有技术中提供的比对算法进行比对，本发明对此不做具体限定。

图3是根据本发明实施例的一种可选的循环肿瘤DNA重复序列的处理方法的流程图，下面结合图3对本发明一种优选的实施例进行详细说明。如图3所示，该方法可以包括如下步骤：输入cfDNA样本捕获测序fastq文件和人类参考基因组fasta文件，利用bwa mem软件进行基因组比对；调用Picard软件进行reads排序；调用samtools软件进行reads过滤；利用本发明上述实施例提供的MostP算法进行去重，得到一致性序列；利用bwa mem软件进行基因组比对；调用Picard软件进行reads排序，得到cfDNA样本去除重复后的bam文件；输入捕获测序区间Bed文件，调用samtools mpileup对标记重复后的bam文件按位置展示所有reads的比对情况和质量值；调用varscan2mpileup2snp模块鉴定SNV，mpileup2indel模块鉴定INDEL。

需要说明的是，上述的cfDNA样本也可以是其他含有ctDNA的体液样本。

本发明输入文件包括：待测样本经过比对、排序、过滤等步骤后生成的测序数据文件(bam格式，包含每条测序片段的名称、SAM标记、位置信息、比对质量信息、CIGAR字串、mate pair信息、片段序列、测序质量等)、人类参考基因组序列(fasta格式)；

本发明的输出文件包括：待测样本标记重复后的比对结果文件(bam格式)以及检测到的SNV和INDEL的vcf格式文件。

通过上述方案，本发明通过对相同位置上所有PE reads逐个碱基的比较，并挑选每个位置上比例最多的碱基类型，构建一个新的一致性序列，最大限度的矫正了大部分的随机错误，保证变异检测的准确性。对于DNA分子碎片化严重、覆盖基因组范围小、经过多轮PCR的样本或测序方案，尤其是血浆ctDNA样本的捕获测序数据可以有效利用碱基序列，提高了原始数据的利用率，和最终变异检测的准确性。

下面通过单碱基变异(SNV)梯度稀释细胞系测试实验验证对上述实施例进行验证。

1、细胞系培养

细胞系HCT116、KYSE450、NCI-H1573、NCI-H1975、NCI-H441、PC-9、SK-HEP-1、SW48、THP-1、BEAS-2B购买自南京科佰生物科技有限公司，按照提供的说明书进行细胞培养，即RPMI-1640培养基中加入10％胎牛血清，在37度条件下进行培养。

2、细胞DNA提取

收集细胞悬液后，常温300g离心5分钟后弃上清，用200uLPBS重悬细胞，然后用QIAamp DNA Mini Kit(货号为51304；Qiagen,Germany)进行基因组DNA提取。经过裂解后过柱纯化，最后用low-TE缓冲液洗脱DNA。

3、用ddPCR的方法确定以上细胞系中突变位点的理论VAF

用细胞提取的基因组DNA作为模板，进行ddPCR的实验，以上细胞系中突变位点的理论VAF如表1所示。ddPCR用伯乐的仪器、商品化探针和反应体系。反应体系组成为：10ulddPCR supermix for probes(no dUTP)，1ul突变探针，1ul野生型探针，以及20ng待测DNA。配制好反应体系后，按照仪器使用方法进行乳糜生成，吸取乳糜至96孔PCR板，用Pierceable Foil Heat Seal进行热封。PCR反应的条件为：酶激活95度，8min；94度30s解链，55度1min退火延伸，共39个循环；酶失活98度10min；4度保温。PCR扩增之后，伯乐的微滴读取仪读取每个反应孔中的带有荧光的微滴数目。每批次反应用超纯水作为阴性对照。每个待测DNA做三个复孔作为技术重复。

表1

4、含有11个突变位点的样本制备

按照下表2中的质量百分比混合上表中的10种细胞系，制备成1个样本，并计算预期的VAF值。

表2

5、样本的ddPCR结果

用ddPCR实验的方法检测样本中以上列表中各个位点的VAF值，如表3所示，每个反应体系中加入20ng样本DNA，每个样本做三个复孔作为技术重复。

表3

基因	突变	DDPCR VAF
			KRAS	G13D	0.53
PIK3CA	H1047R	1.06
			EGFR	G719S	0.88
NRAS	Q61K	1.80
			EGFR	L858R	1.26
EGFR	T790M	1.52
			KRAS	G12V	1.43
EGFR	E746_A750del	4.76
			BRAF	V600E	0.92
EGFR	S768I	2.42
			NRAS	G12D	4.48

6、样本的文库构建、捕获和测序

将混合的细胞系样本DNA首先用covaris超声打断成200bp左右的DNA片段，qubit荧光定量后，如表4所示，用不同的起始量DNA，不足50ul用无酶水补平，采用KAPA hyperpreparation kit(罗氏公司)进行文库构建，经过末端修复、3’端加polyA、连接测序接头、进行无偏向扩增，之后进行纯化获得文库。

表4

样本	起始量DNA(ng)	PCR循环数
			样本1	20	6
样本2	5	8
			样本3	5	8

详述如下：

1)末端平齐并在3’末端加A：反应体系如下表5所示：

表5

试剂	体积
		Fragmented，double-stranded DNA	50μL
End Repair&A-Tailing Buffer	7μL
		End Repair&A-Tailing Enzyme Mix	3μL
总体积	60μL

Buffer和酶应预先在EP管中混匀，与DNA涡旋混匀后按以下反应进行。反应步骤如下表6所示：

表6

该步操作将PCR管盖温度设为85℃，而非105℃。若该操作结束后立即进行下步实验，应将终止温度设为20℃，而非4℃。

2)连接接头：根据建库说明书的指导，20ng DNA应该采用7.5uM接头。按照下表7所示配制反应体系：

表7

试剂	体积
		反应产物	60μL
接头体积	5μL
		超纯水	5μL
连接Buffer	30μL
		DNA连接酶	10μL
总体积	110μL

Buffer和酶应预先在EP管中混匀，涡旋震荡后离心，20℃孵育15分钟。

3)连接后纯化：在上一步反应体系(110ul)中加入Agencourt AMPure XP纯化磁珠88ul。

充分涡旋振荡，轻微离心。室温吸附5-15分钟，使DNA与磁珠充分结合EP管放至磁力架吸附至液体澄清缓慢吸取EP管中上清并丢弃。EP管中加入200μL 80％乙醇孵育30秒缓慢吸取EP管中乙醇并丢弃。重复一次乙醇洗磁珠。EP管室温干燥3-5分钟至乙醇完全挥发。从磁力架取下EP管，加入22μL超纯水，涡旋振荡，轻微离心室温孵育2分钟洗脱DNA，EP管放至磁力架吸附至液体澄清,上清转移至新的EP管，取1μL上清测DNA浓度，剩余的进行扩增。

4)PCR扩增：按照下表8所示配制PCR体系。

表8

试剂	体积
		KAPA HiFi HotStart ReadyMix(2X)	25μL
KAPA Library Amplification Primer Mix(10X)*	5μL
		接头连接文库	20μL
总体积	50μL

充分震荡后快速离心，按照下表9所示条件进行PCR反应。

表9

5)扩增后纯化：加入与PCR反应体系同等体积的Agencourt AMPure XP纯化磁珠(50μl)。

充分涡旋振荡，轻微离心，室温吸附5-15分钟，使DNA与磁珠充分结合。EP管放至磁力架吸附至液体澄清，缓慢吸取EP管中上清并丢弃。EP管中加入200μL 80％乙醇孵育30秒，缓慢吸取EP管中乙醇并丢弃。重复一次乙醇洗磁珠。EP管室温干燥3-5分钟至乙醇完全挥发。从磁力架取下EP管，加入52μL超纯水，涡旋振荡，轻微离心。室温孵育2分钟洗脱DNA，EP管放至磁力架吸附至液体澄清，上清转移至新的EP管，取1μL上清测DNA浓度。

6)在测序前采用探针捕获的方法，用Roche NimbleGen探针将包含11个突变位点的目的区域进行富集和进一步扩增，获得目的区域的文库。经过q-PCR定量后进行上机测序。

7、处理下机fastq数据为各软件可使用的输入文件。

数据下机后，首先将下机数据从fastq文件处理成bam文件，具体使用的软件和步骤如下：

7.1比对

调用bwa-0.7.12mem将每一对fastq文件都作为PE reads比对到hg19人类参考基因组序列，除-M参数与指定Reads Group的ID外，不使用其余参数选项，生成初始bam文件。

7.2排序

调用picard-2.1.0的SortSam模块，对初始bam文件按照染色体位置进行排序，参数设置为“SORT_ORDER＝coordinate”。

7.3筛选

调用samtools-1.3view对排序后的bam文件进行筛选，采用“-F 0x900”作为参数。

7.4建立索引

调用samtools-1.3的index模块对最终生成的bam文件建立索引，生成与过滤后的bam文件配对的bai文件。

8、标记重复

8.1使用Picard’s MarkDuplicates模块标记重复，后续的变异检测时，会自动过滤这部分重复序列，再进行分析。

8.2根据本发明上述实施例提供的方法(MostP)对过滤后的bam文件去除重复序列，生成去除重复的bam文件。

8.3统计比对情况：

调用samtools-1.3的flagstat模块对最终生成的bam文件进行统计，生成去除重复后的bam文件的比对情况文件，包括总reads的数量、重复reads的数量、比对到参考基因组上的reads数量、成对的reads数据数量、read1的数量、read2的数量、完美匹配到参考序列的reads数量(properly paired)、一对reads都比对到了参考序列上的数量、一对reads中只有一条与参考序列相匹配的数量、一对reads比对到不同染色体的数量、一对reads比对到不同染色体的且比对质量值大于5的数量等。

8.4结果比较：

本发明上述实施例提供的算法与Picard方法的数据量统计结果如下表10所示，从下表10中可以看出，本发明提供的算法比Picard方法保留的数据量更多，提高了数据的有效利用率。

表10

样本	Picard	MostP
			样本1	24872747	24821496
样本2	13687626	14068170
			样本3	14290322	14477716

9、变异检测

9.1堆叠

调用samtools-1.3mpileup对去除重复后的bam文件按位置展示所有reads的比对情况和质量值，参数设置为“q＝1”，mpileup的结果文件(mpileup文件)包含染色体、基因组位置、参考基因组碱基类型、该位点测序深度、全部覆盖该位点reads的比对情况和质量值。

由于ddPCR验证阳性位点有限，仅对下列区间做mpileup处理，使用参数“-lpositive.bed”，positive.bed文件如表11所示。

表11

染色体	起始位置	结束位置	基因
				chr1	115256527	115256530	NRAS
chr1	115258745	115258748	NRAS
				chr3	178952083	178952086	PIK3CA
chr12	25398279	25398282	KRAS
				chr12	25398282	25398285	KRAS
chr7	140453134	140453137	BRAF
				chr7	55241706	55241709	EGFR
chr7	55242414	55242513	EGFR
				chr7	55249003	55249006	EGFR
chr7	55249069	55249072	EGFR
				chr7	55259513	55259516	EGFR

9.2统计positive.bed区间的平均测序深度

使用简单的脚本或bash命令根据mpileup文件统计不同去除重复序列方法在positive.bed区间的测序深度的平均值，结果见表12。

表12

样本	Picard	MostP
			样本1	1625.370	1716.500
样本2	533.496	652.767
			样本3	627.380	717.372

本发明提供的方法比Picard的方法在positive.bed区间平均深度略高。

9.3变异检测

调用varscan2mpileup2snp模块检测单核苷酸变异(SNV)，mpileup2indel模块检测插入缺失标记(INDEL)，参数设置“--min-coverage 100--min-reads2 2--min-var-freq0.001--p-value 0.05--min-avg-qual 20”。

对上述3个样本的ddPCR验证为阳性的位点用不同去重方法之后统计的变异结果如下表13至15所示(表格中数值为突变频率)，其中，表13示出样本1的变异结果，表14示出样本2的变异结果，表15示出样本3的变异结果。

表13

基因	Aachange	Picard	MostP
				NRAS	p.Q61K	0	1.29
PIK3CA	p.H1047R	0.96	1.08
				BRAF	p.V600E	0.83	0.82
NRAS	p.G12D	3.87	3.78
				EGFR	p.G719S	0.88	0.83
EGFR	p.L858R	1.64	1.6
				EGFR	p.S768I	2.15	2.31
KRAS	p.G13D	0.6	0.52
				EGFR	p.745_750del	3.05	3.1
KRAS	p.G12V	1.02	1.24
				EGFR	p.T790M	1.39	1.22

表14

基因	Aachange	Picard	MostP
				NRAS	p.Q61K	4.22	3.24
PIK3CA	p.H1047R	0	0.94
				BRAF	p.V600E	0	0
NRAS	p.G12D	0	0
				EGFR	p.G719S	0.93	1.08
EGFR	p.L858R	2.3	2.33
				EGFR	p.S768I	1.04	1.1
KRAS	p.G13D	1.07	0.98
				EGFR	p.745_750del	2.92	2.67
KRAS	p.G12V	1.34	1.3
				EGFR	p.T790M	0.96	1.01

表15

基因	Aachange	Picard	MostP
				NRAS	p.Q61K	0	1.41
PIK3CA	p.H1047R	0	0
				BRAF	p.V600E	0.99	1.21
NRAS	p.G12D	5.45	6.2
				EGFR	p.G719S	0	1.07
EGFR	p.L858R	0.76	0.64
				EGFR	p.S768I	1.66	1.93
KRAS	p.G13D	0	0
				EGFR	p.745_750del	2.56	2.93
KRAS	p.G12V	1.54	2.11
				EGFR	p.T790M	0	0.71

Picard在多处阳性位点检测的突变频率为0(频率>0为阳性，频率＝0为阴性)，MostP在部分位点检测为阴性，发生漏检，但个数少于Picard漏检个数。综上可以看出使用本发明相比Picard去重可以检测更多的阳性位点。

实施例2

根据本发明实施例，提供了一种循环肿瘤DNA重复序列的处理装置的实施例。

图4是根据本发明实施例的一种循环肿瘤DNA重复序列的处理装置的示意图，如图4所示，该装置包括：

获取模块42，用于获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，测序数据包括：多对双端序列。

具体地，上述的待检测循环肿瘤DNA可以从病人的血液、淋巴液、组织间隙液、脑髓液等体液中提取得到的基因序列，在本发明实施例中以血液中提取到的ctDNA为例进行说明；上述的测序数据可以是对待检测ctDNA进行NGS测序得到的ctDNA样本捕获测序fastq数据；上述的参考基因组序列可以是从公开数据库NCBI等网站下载的人类参考基因组fasta数据。

比对模块44，用于将测序数据和参考基因组序列进行比对，得到第一比对结果，其中，第一比对结果至少包括：多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列。

具体地，上述的基因组位置可以是每对PE reads比对到参考基因组序列中的位置，不同的PE reads可以比对到相同的位置；上述的碱基质量值可以是通过NGS测序得到的测序质量，用于衡量每个碱基位置上碱基类型测量的准确度，碱基质量值越大，说明碱基类型测量的准确度越高；上述的碱基序列可以是每对双端序列中每个碱基位置上的碱基类型，DNA序列中包含四种类型的碱基，分别为G、C、T、A，在NGS测序过程中，可以确定每个碱基位置上的碱基类型，并得到该碱基类型的碱基质量值。

在一种可选的方案中，可以获取人类参考基因组fasta数据和ctDNA样本捕获测序fastq数据，利用基因组比对工具bwa mem进行序列比对，得到比对结果文件(.bam)，也即，得到上述的第一比对结果，比对结果文件为bam格式，包含每对PE reads的名称、位置信息、SAM标记、比对质量信息、CIGAR字串、mate pair信息、片段序列、测序质量等)。

需要说明的是，第一比对结果中的多对双端序列的碱基序列和碱基质量值是从NGS测序的测序数据中直接继承过来的数据，第一比对结果包含基因组比对位置及比对情况的信息的同时，还存储了多对双端序列的碱基序列和碱基质量值，方便后续的其他分析，不再使用fastq文件。

处理模块46，用于基于第一比对结果和每对双端序列中每个碱基位置上的碱基类型，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

在一种可选的方案中，由于PCR和测序错误是随机发生的小概率事件，所以其产生的碱基类型相比正确的碱基占比较少，因此，本发明实施例提供了一种MostP算法，基于上述基本事实，比对到基因组相同位置就可以认为来自同一个原始的fragment，且同一组PEreads内，占大多数的碱基类型构成新的一对reads(也即上述的一致性序列)，并将一致性序列作为最终代表的unique reads。

需要说明的是，无论是通过本发明提供的去重方法进行去重之前，还是通过本发明提供的去重方法进行去重之后，所有的bam文件均包含下列信息：每对PE reads的名称、位置信息、SAM标记、比对质量信息、CIGAR字串、mate pair信息、片段序列、测序质量等。

根据本发明上述实施例，获取待检测循环肿瘤DNA的测序数据和参考基因组序列，将测序数据和参考基因组序列进行比对，得到第一比对结果，进一步基于第一比对结果和每对双端序列中每个碱基位置上的碱基类型，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列，从而实现重复序列的去重处理。容易注意到的是，由于在将测序数据和参考基因组序列进行比对之后，需要结合每对双端序列中每个碱基位置上的碱基类型，得到一致性序列，从而实现在考虑到序列质量值的同时，考虑到具体序列上的差异，达到减少人工错误，提高处理准确度的技术效果，进而解决了现有技术中测序数据的处理方法对样本测序进行重复序列删除或标记，准确度低的技术问题。

实施例3

根据本发明实施例，提供了一种存储介质的实施例，存储介质包括存储的程序，其中，在程序运行时控制存储介质所在设备执行上述的循环肿瘤DNA重复序列的处理方法。

实施例4

根据本发明实施例，提供了一种处理器的实施例，处理器用于运行程序，其中，程序运行时执行上述的循环肿瘤DNA重复序列的处理方法。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

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

在本申请所提供的几个实施例中，应该理解到，所揭露的技术内容，可通过其它的方式实现。其中，以上所描述的装置实施例仅仅是示意性的，例如所述单元的划分，可以为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，单元或模块的间接耦合或通信连接，可以是电性或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可为个人计算机、服务器或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (10)

1.一种循环肿瘤DNA重复序列的处理方法，其特征在于，包括：

获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，所述测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，所述测序数据包括：多对双端序列；

将所述测序数据和所述参考基因组序列进行比对，得到第一比对结果，其中，所述第一比对结果至少包括：所述多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列；

基于所述第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

2.根据权利要求1所述的方法，其特征在于，基于所述第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列包括：

将所述多对双端序列划分为至少一个序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同；

将每个序列集合中的每对双端序列逐碱基进行比较，得到所述每个序列集合中每个碱基位置上每个碱基类型的数量；

获取所述每个序列集合中每个碱基位置上数量超过预设数量的第一碱基类型；

将所述每个序列集合中所述每个碱基位置上的第一碱基类型进行组合，得到所述每个序列集合对应的一致性序列；

根据每个一致性序列的每个碱基位置上第一碱基类型的碱基质量值及除所述第一碱基类型之外的其他碱基类型的碱基质量值，得到所述每个一致性序列对应的碱基质量值序列。

3.根据权利要求1所述的方法，其特征在于，在基于所述第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列之后，所述方法还包括：

将所述至少一个一致性序列和所述参考基因组序列进行比对，得到第二比对结果；

按照每个一致性序列的基因组位置，对所述第二比对结果进行排序，得到第三比对结果。

4.根据权利要求3所述的方法，其特征在于，在按照每个一致性序列的基因组位置，对所述第二比对结果进行排序，得到第三比对结果之后，所述方法还包括：

根据所述第三比对结果，显示每个基因组位置对应的双端序列的比对信息和碱基质量值；

对比对质量满足预设条件的双端序列进行过滤。

5.根据权利要求3所述的方法，其特征在于，在按照每个一致性序列的基因组位置，对所述第二比对结果进行排序，得到第三比对结果之后，所述方法还包括：

获取捕获测序区间；

根据所述捕获测序区间，对所述每个一致性序列进行单核苷酸变异检测和插入缺失检测，得到检测结果。

6.根据权利要求1所述的方法，其特征在于，在基于所述第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列之前，所述方法还包括：

按照所述多对双端序列的基因组位置，对所述第一比对结果进行排序，得到第四比对结果，并为所述第四比对结果建立索引；

对所述第四比对结果进行过滤，得到第五比对结果；

基于所述第五比对结果，得到所述至少一个一致性序列和所述每个一致性序列对应的碱基质量值序列。

7.根据权利要求1所述的方法，其特征在于，将所述测序数据和所述参考基因组序列进行比对，得到第一比对结果包括：

获取所述多对双端序列中每条序列和所述参考基因组序列中的每段序列的匹配度；

获取最高匹配度对应的至少一段序列，得到所述每条序列的匹配序列；

根据所述每条序列的匹配序列，确定所述条序列的基因组位置。

8.一种循环肿瘤DNA重复序列的处理装置，其特征在于，包括：

获取模块，用于获取待检测循环肿瘤DNA的测序数据和参考基因组序列，其中，所述测序数据为对待检测循环肿瘤DNA进行高通量测序得到的数据，所述测序数据包括：多对双端序列；

比对模块，用于对所述测序数据和所述参考基因组序列进行比对，得到第一比对结果，其中，所述第一比对结果至少包括：所述多对双端序列的基因组位置、碱基序列和对应的碱基质量值序列；

处理模块，用于基于所述第一比对结果，得到至少一个一致性序列和每个一致性序列对应的碱基质量值序列。

9.一种存储介质，其特征在于，所述存储介质包括存储的程序，其中，在所述程序运行时控制所述存储介质所在设备执行权利要求1至7中任意一项所述的循环肿瘤DNA重复序列的处理方法。

10.一种处理器，其特征在于，所述处理器用于运行程序，其中，所述程序运行时执行权利要求1至7中任意一项所述的循环肿瘤DNA重复序列的处理方法。