CN108319817B - 循环肿瘤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的情况也会做去重处理。

但在基因组相同位置上，往往有可能会存在多个原始分子，这些原始分子并不是通常意义上的PCR重复，有可能存在有意义的突变(例如ctDNA中就是肿瘤相关变异)，但在上述的去重方法中，对于这种情况的判断，samtools rmdup和Picard’s MarkDuplicates会错误的认为是同一个原始分子，仅保留1对reads，导致过度去重，浪费了部分有意义的数据量。并且由于PCR过程可能伴有随机错误的产生，这些错误最后很可能造成假阳性的突变检出，而samtools rmdup和Picard’s MarkDuplicates并没有方法对其进行校正。无论是Samtools rmdup还是Picard’s MarkDuplicates都只考虑了read的某种质量值，而并没有考虑read上的具体序列上的差异，导致过去重或错误选取保留的unique reads的发生。

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

发明内容

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

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

进一步地，对每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络包括：将每个序列集合中的至少一对双端序列进行比较，得到至少一个结点以及每个结点的成员数，其中，同一个结点对应的双端序列的碱基序列相同，成员数用于表征同一个结点对应的双端序列的数量；获取每个序列集合中任意两个结点之间的编辑距离；判断任意两个结点之间的编辑距离是否小于预设距离；如果任意两个结点之间的编辑距离小于预设距离，则确定任意两个结点属于同一个网络，得到至少一个第一网络。

进一步地，获取每个第一网络中数量最多的双端序列，得到每个第一网络对应的独立序列包括：步骤A，获取每个第一网络中成员数最多的第一结点，计算第一结点对应的每对双端序列包含的所有碱基的碱基质量值之和，得到第一结点对应的每对双端序列的碱基质量和，并将最大碱基质量和对应的双端序列作为每个第一网络对应的独立序列；步骤B，从每个第一网络中将第一结点和与第一结点相邻的结点进行删除，得到至少一个第二网络；步骤C，判断任意一个第二网络是否包含结点；步骤D，如果任意一个第二网络包含结点，则将任意一个第二网络作为第一网络，并循环执行步骤A和步骤B。

进一步地，在确定任意两个结点属于同一个网络，得到至少一个第一网络之前，上述方法还包括：判断任意两个结点的成员数是否满足预设条件；如果任意两个结点的成员数满足预设条件，则确定任意两个结点属于同一个网络，得到至少一个第一网络。

进一步地，判断任意两个结点的成员数是否满足预设条件包括：获取任意两个结点中第一结点的成员数的预设倍数，得到第一数值；获取第一数值与预设值的差值；判断任意两个结点中第二结点的成员数是否大于差值；如果第二结点的成员数大于差值，则确定任意两个结点的成员数满足预设条件。

进一步地，在获取每个第一网络中数量最多的双端序列，得到每个第一网络对应的独立序列之后，上述方法还包括：将第一比对结果中除每个第一网络对应的独立序列之外的其他双端序列的比对结果删除，得到第二比对结果。

进一步地，在将第一比对结果中除每个第一网络对应的独立序列之外的其他双端序列的比对结果删除，得到第二比对结果之后，上述方法还包括：按照每个独立序列的基因组位置，对第二比对结果进行排序，得到第三比对结果。

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

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

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

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

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

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

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

附图说明

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

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

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

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

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

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

图6是根据本发明实施例的一种循环肿瘤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所示，该方法包括如下步骤：

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

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

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

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

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

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

步骤S104，基于第一比对结果，得到第一序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同。

步骤S106，对每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络。

步骤S108，获取每个第一网络中数量最多的双端序列，得到每个第一网络对应的独立序列。

在一种可选的方案中，由于在没有PCR和测序错误时，同一个DNA分子片段(fragment)经过PCR产生多个完全一样的fragment，这组fragment均可以比对到参考基因组同一个位置，且序列相同；PCR和测序错误是随机发生的小概率事件，所以其产生的碱基序列相比正确的碱基序列占比较少；而来自不同fragment的DNA分子，虽然可能比对到基因组的相同位置，但由于其可能分别属于不同的DNA分子(例如，血液中提取到的DNA包含两种类型，一种是包含肿瘤信息的ctDNA分子；另一种是在血液中游离的自身DNA，多是从身体的细胞或者白血球破裂释放出来的，一般认为是无害的，不用多久会被人体自身清理掉，两种DNA携带的信息不同；又例如，不同的ctDNA分子)，所以其序列可能并不相同。基于上述基本事实，可以将比对到基因组相同位置的所有PE reads的碱基序列做网络聚类，两两之间的编辑距离越小的越容易聚到同一个类别中，而编辑距离大的更可能是来自另两个不同的fragment的数据，从而被分到两个不同的类中，从而得到至少一个第一网络，进一步可以选择每个第一网络中选择最多的一对序列作为最终代表的unique reads，其余的认为是重复序列。通过上述方案在保留更多的原始分子的时候也降低了最终选择带有PCR或测序错误序列的可能性，提高了结果的准确性。

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

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

可选地，在本发明上述实施例中，步骤S106，对每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络包括：

步骤S1061，将每个序列集合中的至少一对双端序列的碱基序列进行比较，得到至少一个结点以及每个结点的成员数，其中，同一个结点对应的双端序列的碱基序列相同，成员数用于表征同一个结点对应的双端序列的数量。

具体地，可以以PE reads为结点，将序列相同的PE reads确定为一个结点，并根据序列相同的PE reads的数量得到该结点的成员数，从而进一步可以对每个序列集合中的结点进行网络聚类。

步骤S1062，获取每个序列集合中任意两个结点之间的编辑距离。

具体地，任意两个结点之间的编辑距离可以表征任意两个结点的相似性，编辑距离越小，相似性越高，任意两个结点越有可能来自同一个fragment；编辑距离越大，相似度越低，越有可能是来自两个不同的fragment。

需要说明的是，可以采用现有的编辑距离的处理方法获取任意两个结点之间的编辑距离，本发明对此不做具体限定。

步骤S1063，判断任意两个结点之间的编辑距离是否小于预设距离。

具体地，上述的预设距离可以根据实际需要设定的阈值，例如，预设距离可以是1。

步骤S1064，如果任意两个结点之间的编辑距离小于预设距离，则确定任意两个结点列属于同一个网络，得到至少一个第一网络。

在一种可选的方案中，本发明实施例提供了第一种网络聚类方法，也即Cluster算法，可以将比对到基因组相同位置且序列间编辑距离小于阈值(默认为1)的所有序列连接成一个网络，以序列为结点，每个结点的成员数即代表这种序列在这个位置上出现的次数，整个网络就可以认为来自同一个原始的fragment，且同一组内，选择成员数最多的一种序列，并确定这种序列中碱基质量和最高的一个序列作为最终的代表unique reads。

例如，如图2所示，对于如图2所示的6种PE reads，分别为AGTC、CGTC、TGTC、TGAC、TAAC、TGAG，假设6种PE reads的基因组位置相同，且编辑距离小于1，该6种PE reads连接成一个网络，进一步根据每个结点的成员数，也即AGTC的成员数为2，CGTC的成员数为2，TGTC的成员数为365，TGAC的成员数为62，TAAC的成员数为80，TGAG的成员数为1，可以确定最终的代表unique reads为碱基质量和最高的TGTC。

可选地，在本发明上述实施例中，步骤S108，获取每个第一网络中数量最多的双端序列，得到每个第一网络对应的独立序列包括：

步骤A，获取每个第一网络中成员数最多的第一结点，计算第一结点对应的每对双端序列包含的所有碱基的碱基质量值之和，得到第一结点对应的每对双端序列的碱基质量和，并将最大碱基质量和对应的双端序列作为每个第一网络对应的独立序列。

步骤B，从每个第一网络中将第一结点和与第一结点相邻的结点进行删除，得到至少一个第二网络。

步骤C，判断任意一个第二网络是否包含结点。

具体地，当第二网络不包含结点时，可以确定第一网络中所有的结点均遍历完全。

步骤D，如果任意一个第二网络包含结点，则将任意一个第二网络作为第一网络，并循环执行步骤A和步骤B。

在一种可选的方案中，本发明实施例提供了第二种网络聚类方法，也即Adjacency算法，可以将比对到基因组相同位置且序列间编辑距离小于阈值(默认为1)的所有序列连接成一个网络，以序列为结点，每个结点的成员数即代表这种序列在这个位置上出现的次数，整个网络中的序列可能来自多个原始的fragment，进一步地，uniqe reads的挑选规则包含下面几步：首选从网络中挑选成员数最多的结点，从该节点中选择碱基质量和最高的双端序列作为第一个uniqe reads；然后从网络中删掉所有与这个结点直接相连的结点；进一步挑选剩余结点中成员数最多的结点，从该节点中选择碱基质量和最高的双端序列作为第二个代表reads，再删掉所有与它直接相连的结点；重复上述方法直到遍历完全部的结点，从而得到所有的独立序列。

例如，如图3所示，对于如图3所示的6种PE reads，分别为AGTC、CGTC、TGTC、TGAC、TAAC、TGAG，假设6种PE reads的基因组位置相同，且编辑距离小于1，该6种PE reads连接成一个网络，进一步根据每个结点的成员数，也即AGTC的成员数为2，CGTC的成员数为2，TGTC的成员数为365，TGAC的成员数为62，TAAC的成员数为80，TGAG的成员数为1，可以获取最大成员数且碱基质量和最高的序列，也即，获取碱基质量和最高的TGTC作为独立序列，进一步删除与TGTC直接相连的序列，也即删除AGTC、CGTC和TGAC，剩余的序列为TAAC和TGAG，获取最大成员数且碱基质量和最高的序列，也即，获取碱基质量和最高的TAAC作为独立序列，由于TAAC和TGAG并未直接相连，因此，可以将碱基质量和最高的TAAC和碱基质量和最高的TGAG均作为独立序列，因此，可以得到三个独立序列，分别为碱基质量和最高的TGTC、碱基质量和最高的TAAC和碱基质量和最高的TGAG。

可选地，在本发明上述实施例中，在步骤S1064，确定任意两个结点列属于同一个网络，得到至少一个第一网络之前，该方法还包括：

步骤S1065，判断任意两个结点的成员数是否满足预设条件。

具体地，上述的预设条件用于表征任意两个结点可以连接成一个网络的条件。

步骤S1066，如果任意两个结点的成员数满足预设条件，则确定任意两个结点属于同一个网络，得到至少一个第一网络。

在一种可选的方案中，本发明实施例提供了第三种网络聚类方法，也即Direct算法，可以将比对到基因组相同位置且序列间编辑距离小于阈值(默认为1)且相邻结点a和结点b满足一定的条件的所有序列连接成一个网络，以序列为结点，每个结点的成员数即代表这种序列在这个位置上出现的次数，这样在这个位置处可以形成一个或多个网络，每个网络中的序列来自通一个原始的fragment，选择成员数最多的一种序列，并确定这种序列中碱基质量和最高的一个序列作为最终的代表。

例如，如图4所示，对于如图4所示的6种PE reads，分别为AGTC、CGTC、TGTC、TGAC、TAAC、TGAG，其中，AGTC和CGTC相邻，CGTC和TGTC相邻，TGTC和AGTC相邻，TGTC和TGAC相邻，TGAC和TAAC相邻，TGAC和TGAG相邻，假设6种PE reads的基因组位置相同，且编辑距离小于1，而且AGTC、CGTC、TGTC、TGAC和TGAG中相邻两结点满足一定的条件，则可以将AGTC、CGTC、TGTC、TGAC和TGAG连接成一个网络，TAAC为单独的网络。进一步根据每个结点的成员数，也即AGTC的成员数为2，CGTC的成员数为2，TGTC的成员数为365，TGAC的成员数为62，TAAC的成员数为80，TGAG的成员数为1，可以获取每个网络中最大成员数且碱基质量和最高的序列，也即，获取碱基质量和最高的TGTC和碱基质量和最高的TAAC作为独立序列。

可选地，在本发明上述实施例中，步骤S1065，判断任意两个结点的成员数是否满足预设条件包括：

步骤S10652，获取任意两个结点中第一结点的成员数的预设倍数，得到第一数值。

步骤S10654，获取第一数值与预设值的差值。

具体地，上述的预设值可以根据实际需要进行设定，例如，可以是1。

步骤S10656，判断任意两个结点中第二结点的成员数是否大于差值。

步骤S10658，如果第二结点的成员数大于差值，则确定任意两个结点的成员数满足预设条件。

在一种可选的方案中，可以将将比对到基因组相同位置且序列间编辑距离小于某一阈值(默认为1)的且相邻结点a和结点b的成员数满足：结点a的成员数≥2*结点b的成员数-1的所有结点连接成一个网络。

可选地，在本发明上述实施例中，在步骤S108，获取每个第一网络中数量最多的双端序列，得到每个第一网络对应的独立序列之后，该方法还包括：

步骤S110，将第一比对结果中除每个第一网络对应的独立序列之外的其他双端序列的比对结果删除，得到第二比对结果。

在一种可选的方案中，可以将所有PE reads中除独立序列之外的重复序列全部删除，仅保留独立序列，避免重复序列引起的人工误差。

可选地，在本发明上述实施例中，在步骤S110，将第一比对结果中除每个第一网络对应的独立序列之外的其他双端序列的比对结果删除，得到第二比对结果之后，该方法还包括：

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

在一种可选的方案中，由于根据编辑距离进行网络聚类，比对结果文件(.bam)中的序列位置发生变化，为了使得相同位置的PE reads相邻，方便后续对PE reads进行去重处理，可以调用Picard’s SortSam模块将比对结果文件(.bam)(也即上述的第二比对结果)按比对位置排序，得到第三比对结果。

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

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

步骤S116，根据捕获测序区间，对每个独立序列进行单核苷酸变异检测和***缺失检测，得到检测结果。

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

可选地，在本发明上述实施例中，在步骤S104，基于第一比对结果，得到第一序列集合之前，该方法还包括：

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

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

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

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

步骤S122，基于第五比对结果，得到至少一个序列集合。

在一种可选的方案中，可以将第五比对结果中，比对到基因组相同位置划分为一个序列集合，从而方便后续对每个序列集合中的所有PE reads的碱基序列做网络聚类。

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

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

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

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

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

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

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

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

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

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

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

通过上述方案，通过网络聚类方法不仅可以保留更多的原始分子，而且可以杜绝大部分随机错误的影响，从而在一定程序上减少随机错误在最终变异检测时候的影响，使得变异检测的更准确、可靠。对于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根据本发明上述实施例提供的方法，分别调用"Cluster"、"Adjacency"和"Direct"方法对过滤后的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方法保留的数据量更多，提高了数据的有效利用率，并且，满足Adjacency>Direct>Cluster>Picard。

表10

样本	Picard	Cluster	Adjacency	Direct
					样本1	24872747	51683842	51930648	51698546
样本2	13687626	46207524	46730986	46247784
					样本3	14290322	44097044	44541914	44129370

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	Cluster	Adjacency	Direct
					样本1	1625.370	4111.270	4130.270	4112.240
样本2	533.496	3251.800	3302.600	3258.130
					样本3	627.380	3084.870	3121.860	3087.800

本发明提供的方法均比Picard的方法在positive.bed区间平均深度更高，且满足Adjacency>Direct>Cluster>Picard。

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	Cluster	Adjacency	Direct
						NRAS	p.Q61K	0	1	1.05	1
PIK3CA	p.H1047R	0.96	1.24	1.23	1.24
						BRAF	p.V600E	0.83	0.81	0.8	0.81
NRAS	p.G12D	3.87	4.22	4.24	4.21
						EGFR	p.G719S	0.88	0.77	0.77	0.77
EGFR	p.L858R	1.64	1.98	1.97	1.98
						EGFR	p.S768I	2.15	2.36	2.34	2.36
KRAS	p.G13D	0.6	0.52	0.51	0.52
						EGFR	p.745_750del	3.05	2.79	2.78	2.79
KRAS	p.G12V	1.02	1.16	1.15	1.15
						EGFR	p.T790M	1.39	1.1	1.09	1.1

表14

基因	Aachange	Picard	Cluster	Adjacency	Direct
						NRAS	p.Q61K	4.22	4.3	4.28	4.3
PIK3CA	p.H1047R	0	1.68	1.67	1.68
						BRAF	p.V600E	0	1.03	1.01	1.02
NRAS	p.G12D	0	1.55	1.52	1.55
						EGFR	p.G719S	0.93	1.14	1.15	1.14
EGFR	p.L858R	2.3	2.23	2.23	2.26
						EGFR	p.S768I	1.04	0.93	0.91	0.93
KRAS	p.G13D	1.07	0.88	0.87	0.88
						EGFR	p.745_750del	2.92	2.05	2.05	2.05
KRAS	p.G12V	1.34	1.85	1.88	1.85
						EGFR	p.T790M	0.96	0.95	0.93	0.95

表15

基因	Aachange	Picard	Cluster	Adjacency	Direct
						NRAS	p.Q61K	0	1.1	1.09	1.1
PIK3CA	p.H1047R	0	0.46	0.46	0.46
						BRAF	p.V600E	0.99	0.98	0.97	0.98
NRAS	p.G12D	5.45	5.86	5.8	5.85
						EGFR	p.G719S	0	1.26	1.27	1.25
EGFR	p.L858R	0.76	0.94	0.97	0.94
						EGFR	p.S768I	1.66	1.55	1.56	1.55
KRAS	p.G13D	0	0.3	0.31	0.3
						EGFR	p.745_750del	2.56	2.52	2.49	2.52
KRAS	p.G12V	1.54	2.06	2.04	2.06
						EGFR	p.T790M	0	0.95	0.93	0.95

Picard在多处阳性位点检测的突变频率为0(频率>0为阳性，频率＝0为阴性)，而Adjacency、Direct和Cluster四种方法在全部11个位点都检测为阳性。综上可以看出使用本发明相比Picard去重可以检测更多的阳性位点，且在保证阳性检出的同时，Direct和Cluster方法去掉更多的重复序列。

实施例2

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

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

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

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

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

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

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

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

处理模块64，用于基于第一比对结果，得到至少一个序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同。

聚类模块66，用于对每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络。

第二获取模块68，用于获取每个第一网络中数量最多的双端序列，得到每个第一网络对应的独立序列。

在一种可选的方案中，由于在没有PCR和测序错误时，同一个DNA分子片段(fragment)经过PCR产生多个完全一样的fragment，这组fragment均可以比对到参考基因组同一个位置，且序列相同；PCR和测序错误是随机发生的小概率事件，所以其产生的碱基序列相比正确的碱基序列占比较少；而来自不同fragment的DNA分子，虽然可能比对到基因组的相同位置，但由于其可能分别属于不同的DNA分子(例如，血液中提取到的DNA包含两种类型，一种是包含肿瘤信息的ctDNA分子；另一种是在血液中游离的自身DNA，多是从身体的细胞或者白血球破裂释放出来的，一般认为是无害的，不用多久会被人体自身清理掉，两种DNA携带的信息不同；又例如，不同的ctDNA分子)，所以其序列可能并不相同。基于上述基本事实，可以将比对到基因组相同位置的所有PE reads的碱基序列做网络聚类，两两之间的编辑距离越小的越容易聚到同一个类别中，而编辑距离大的更可能是来自另两个不同的fragment的数据，从而被分到两个不同的类中，从而得到至少一个第一网络，进一步可以选择每个第一网络中选择最多的一对序列作为最终代表的unique reads，其余的认为是重复序列。通过上述方案在保留更多的原始分子的时候也降低了最终选择带有PCR或测序错误序列的可能性，提高了结果的准确性。

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

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

实施例3

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

实施例4

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

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

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

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

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

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

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

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

Claims (11)

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

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

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

基于所述第一比对结果，得到至少一个序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同；

对所述每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络；

获取每个第一网络中数量最多的双端序列，得到所述每个第一网络对应的独立序列；

其中，对所述每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络包括：将所述每个序列集合中的至少一对双端序列的碱基序列进行比较，得到至少一个结点以及每个结点的成员数，其中，同一个结点对应的双端序列的碱基序列相同，所述成员数用于表征同一个结点对应的双端序列的数量；获取所述每个序列集合中任意两个结点之间的编辑距离；判断所述任意两个结点之间的编辑距离是否小于预设距离；如果所述任意两个结点之间的编辑距离小于所述预设距离，则确定所述任意两个结点属于同一个网络，得到所述至少一个第一网络；

所述获取每个第一网络中数量最多的双端序列，得到所述每个第一网络对应的独立序列包括：步骤A，获取所述每个第一网络中成员数最多的第一结点，计算所述第一结点对应的每对双端序列包含的所有碱基的碱基质量值之和，得到所述第一结点对应的所述每对双端序列的碱基质量和，并将最大碱基质量和对应的双端序列作为所述每个第一网络对应的独立序列；步骤B，从所述每个第一网络中将所述第一结点和与所述第一结点相邻的结点进行删除，得到至少一个第二网络；步骤C，判断任意一个第二网络是否包含结点；步骤D，如果所述任意一个第二网络包含结点，则将所述任意一个第二网络作为所述第一网络，并循环执行所述步骤A和所述步骤B。

2.根据权利要求1所述的方法，其特征在于，在确定所述任意两个结点属于同一个网络，得到所述至少一个第一网络之前，所述方法还包括：

判断所述任意两个结点的成员数是否满足预设条件；

如果所述任意两个结点的成员数满足所述预设条件，则确定所述任意两个结点属于同一个网络，得到所述至少一个第一网络。

3.根据权利要求2所述的方法，其特征在于，判断所述任意两个结点的成员数是否满足预设条件包括：

获取所述任意两个结点中第一结点的成员数的预设倍数，得到第一数值；

获取所述第一数值与预设值的差值；

判断所述任意两个结点中第二结点的成员数是否大于所述差值；

如果所述第二结点的成员数大于所述差值，则确定所述任意两个结点的成员数满足所述预设条件。

4.根据权利要求1所述的方法，其特征在于，在获取每个第一网络中数量最多的双端序列，得到所述每个第一网络对应的独立序列之后，所述方法还包括：

将所述第一比对结果中除所述每个第一网络对应的独立序列之外的其他双端序列的比对结果删除，得到第二比对结果。

5.根据权利要求4所述的方法，其特征在于，在将所述第一比对结果中除所述每个第一网络对应的独立序列之外的其他双端序列的比对结果删除，得到第二比对结果之后，所述方法还包括：

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

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

获取捕获测序区间；

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

7.根据权利要求1所述的方法，其特征在于，在基于所述第一比对结果，得到至少一个序列集合之前，所述方法还包括：

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

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

基于所述第五比对结果，得到所述至少一个序列集合。

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

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

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

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

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

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

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

处理模块，用于基于所述第一比对结果，得到至少一个序列集合，其中，每个序列集合包括：至少一对双端序列，同一个序列集合中的双端序列的基因组位置相同；

聚类模块，用于对所述每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络；

第二获取模块，用于获取每个第一网络中数量最多的双端序列，得到所述每个第一网络对应的独立序列；

其中，对所述每个序列集合中的至少一对双端序列进行网络聚类，得到至少一个第一网络包括：将所述每个序列集合中的至少一对双端序列的碱基序列进行比较，得到至少一个结点以及每个结点的成员数，其中，同一个结点对应的双端序列的碱基序列相同，所述成员数用于表征同一个结点对应的双端序列的数量；获取所述每个序列集合中任意两个结点之间的编辑距离；判断所述任意两个结点之间的编辑距离是否小于预设距离；如果所述任意两个结点之间的编辑距离小于所述预设距离，则确定所述任意两个结点属于同一个网络，得到所述至少一个第一网络；

所述获取每个第一网络中数量最多的双端序列，得到所述每个第一网络对应的独立序列包括：步骤A，获取所述每个第一网络中成员数最多的第一结点，计算所述第一结点对应的每对双端序列包含的所有碱基的碱基质量值之和，得到所述第一结点对应的所述每对双端序列的碱基质量和，并将最大碱基质量和对应的双端序列作为所述每个第一网络对应的独立序列；步骤B，从所述每个第一网络中将所述第一结点和与所述第一结点相邻的结点进行删除，得到至少一个第二网络；步骤C，判断任意一个第二网络是否包含结点；步骤D，如果所述任意一个第二网络包含结点，则将所述任意一个第二网络作为所述第一网络，并循环执行所述步骤A和所述步骤B。

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

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

Images