CN106599616A - 基于duplex‑seq的超低频突变位点检测分析方法 - Google Patents

Abstract

本发明公开的一种基于duplex‑seq的超低频突变位点检测分析方法，包括如下步骤：1)对原始测序数据质量进行评估，降低数据噪声，为后续分析提供有效数据；2)把随机barcode提取到序列文件的每一条序列的标题行，方便后续对barcode进行快速检索并创建一致性序列；3)根据family barcode和duplex barcode创建一致性序列，排除由于建库过程或者PCR过程中引入的突变；4)根据duplex‑tag构建双链一致性序列，进一步排除序列中的非对称突变位点；5)对比对后的数据进行局部质量矫正，并进行低频变异位点检测；将变异位点进行基因结构、功能、及临床表型三个层次的注释；6)统计SSCS、DCS序列数目、比对结果、变异位点信息，并输出可视化图表。

Description

基于duplex-seq的超低频突变位点检测分析方法

技术领域

本发明属于生物信息数据处理方法，特别涉及一种基于duplex-seq的超低频突变位点检测分析方法，其主要应用与第二代高通量测序领域，基于duplex-seq全外显子组测序，对ctDNA的超低频变异位点进行检测分析。

背景技术

下一代测序技术的发展如火如荼，其正在以势如破竹的力量深刻变革着传统遗传学的研究，并因此催生了精准医学的萌芽。相比与传统的实验技术，该技术可以一次性检测成千上万的遗传突变。然而美中不足的是下一代测序技术仍然存在相对较高的错误率(0.1～1％)。对于高频的遗传突变检测，这个错误是可以接受的，但是对于一些低频突变的检测，该测序方法存在很大的局限性。

低频突变的检测现在应用越来越广泛，特别是在肿瘤细胞的检测中的应用倍受关注。肿瘤组织中出现突变的细胞比例小于1％，对于这样的低频突变，一般的方法很难检测到。因此有学者提出了双重标记的测序方法(Duplex Sequencing，Duplex-seq)(KennedySR,Schmitt MW,Fox EJ,Kohrn BF,Salk JJ,Ahn EH,et al.Detecting ultralow-frequency mutations by Duplex Sequencing.Nature protocols.2014；9:2586-606.)和(Newman AM,Bratman SV,To J,Wynne JF,Eclov NC,Modlin LA,et al.Anultrasensitive method for quantitating circulating tumor DNA with broadpatient coverage.Nature medicine.2014；20:548-54.)。该方法能应用于任何双链DNA样本，通过对双链加入随机的双重标签序列，使得每条链都具有唯一的分子标记序列(UniqueMolecular index，UMI)，可来排除由于建库过程或者测序过程中引入的误差。对于此类测序数据的分析方法已有相关报道(Newman AM,Lovejoy AF,Klass DM,Kurtz DM,ChabonJJ,Scherer F,et al.Integrated digital error suppression for improveddetection of circulating tumor DNA.Nature biotechnology.2016；34:547-55)，但是都存在一定的不足，因此，本发明拟通过开发新的方法来对Duplex-seq测序数据进行更加系统专业的分析，以用来指导相关疾病患者在临床中选择更好的药物和剂量，达到精准、高效治疗的目的，以及对健康人群的患病风险的评估提供更为准确的参考。

Duplex-seq技术在高通量测序已得到广泛应用，基于Duplex-seq对低频突变位点进行检测的分析方法也越来越多，但是目前存在以下几点问题：

1.数据质量控制：现有Duplex-seq的数据分析方法流程中，未见对前期数据质量的系统分析，如数据重复率、UMI的种类、数量、比例、R1R2的平衡性等。

2.UMI的差异分析：现有Duplex-seq数据分析方法中对于单链特异性UMI和双链互补的UMI的综合和独立分析还未见报道。

3.变异位点注释流程：基于Duplex-seq数据拷出的突变位点属于低频突变位点，因此对于变异位点的检测相关参数，现有方法还可以进行优化，后续注释流程和相关统计还不够系统。

4.结果可读性:现有Duplex-seq数据分析方法中，结果中只有一些简单的图表文件和文本报告，还有很多数据信息没有有效的呈现。且现有的图表在可视化、形象化角度来说还可以进一步提升。

发明内容

本发明的目的是针对上述存在的不足，提出一种基于duplex-seq的超低频突变位点检测分析方法。其对ctDNA的duplex-seq测序数据进行高效、快捷、准确、深入的信息挖掘分析，并给出可视化的数据结果。

为了实现上述目的，本发明所采用的技术方案如下：

一种基于duplex-seq的超低频突变位点检测分析方法，包括如下步骤：

1)对原始测序数据质量进行评估，降低数据噪声，为后续分析提供有效数据；

2)把随机barcode提取到序列文件的每一条序列的标题行，方便后续对barcode进行快速检索并创建一致性序列；

3)根据family barcode和duplex barcode创建一致性序列，排除由于建库过程或者PCR过程中引入的突变；

4)根据duplex-tag构建双链一致性序列，进一步排除序列中的非对称突变位点；

5)对比对后的数据进行局部质量矫正，并进行低频变异位点检测；将变异位点进行基因结构、功能、及临床表型三个层次的注释；

6)统计SSCS、DCS序列数目、比对结果、变异位点信息，并输出可视化图表。

在本发明的一个优选实施例中，所述步骤1)包括如下步骤：

1.1)查看测序质量结果，包括碱基得分值、GC分布、碱基平衡性、重复率等；统计所有序列的总Total Reads，Total bases、Q20、Q30、GC含量、N字符等数量及相关比例；

1.2)去除序列中含有的接头序列；

1.3)对去除接头后的序列计算其重复度；

1.4)提取每条序列的family barcode和duplex barcode对其种类、数量、比例进行相应统计。

在本发明的一个优选实施例中，所述步骤2)包括如下步骤：

提取配对序列的Duplex barcode、family barcode；

对上述提取的barcode进行过滤，如有非碱基字符则丢弃。

把通过步骤2.2)过滤后的barcode构建duplex-tag。

2.4)根据duplex-tag分别输出配对序列R1，R2；

2.5)对下一条序列重复上述2.1-2.4步骤。

在本发明的一个优选实施例中，所述步骤3)包括如下步骤：

3.1)把步骤(2)中获得的序列，比对到相应的基因组获得比对文件；

3.2)获取比对文件中基因组的第一个位点；

3.3)对所有比对到步骤3.2)提取的位点上的reads进行flag字段的过滤，对于未比对上的reads()则另存为NM文件；其中flag字段为77和141；

3.4)对所有通过步骤3.3)过滤后的reads根据duplex-tag进行排序，并进行分组；

3.5)提取步骤3.4)分组中的第一组Duplex-tag及其相关序列；

3.6)对步骤3.5)中提取的一组序列根据CIGAR string进一步分组归类，对于含有相同CIGAR string的序列则进行下一步分析，对于不含有共同CIGAR string的序列则另存为LCC文件；

3.7)对步骤3.6)中分组的序列，计算其family size，如果family size小于3则丢弃该组序列，通过则进行下一步的分析；

3.8)对步骤3.7)中通过的一组序列创建单链一致性序列；对于碱基一致性较高的位点则该位点归一为含量较高的碱基，对于一致性不够好的序列则该位点以N代替；一致性的值根据用户自行定义，设置为70％；

3.9)通过3.8构建的单链一致性序列，过滤掉含有30％以上N的序列，并输出最终合格的序列到单链一致性(SSCS)文件；

3.10)创建完SSCS后，如果含有更多的duplex-tag，则重复上述步骤3.6)-3.9)，如果没有则进入3.11)步骤；

3.11)若果还有更多的位点，则重复上述步骤3.3)-3.10)，否则结束该模块分析。

在本发明的一个优选实施例中，所述步骤4)包括如下步骤：

4.1)将上述步骤3)中获得的单链一致性序列比对文件转化为sam格式文件，并去除family barcode序列，方便下面创建双链一致性序列；

4.2)提取比对文件中的第一个基因组位点信息；

4.3)提取步骤4.2)中基因组位点对应的第一个duplex tag；

4.4)寻找与步骤4.3)中的duplex-tag进行互补配对的duplex-tag，如果没有对应的duplex-tag与其进行匹配则该序列丢弃，如果有匹配的tag则进入步骤4.5)；

4.5)对含有相同duplex-tag的序列构建双链一致性序列；

4.6)对步骤4.5)创建的双链一致性序列进行过滤，如果序列中含有大量的N则该序列丢弃，否则进入步骤4.7)；

4.7)输出为双链一致性文件(DCS),并输出对应的配对序列(DCS R1 R2)；

4.8在步骤4.5)创建完双链一致性序列后，如果还有更多的duplex-tag则重复步骤4.4)-4.7)，如果没有duplex-tag则进入步骤4.9)；

4.9)如果还有更多的位点需要分析，则重复步骤4.3)-4.8)；否则结束分析。

在本发明的一个优选实施例中，所述步骤5)包括如下步骤：

将步骤4)中生存的双链一致性序列比对到基因组；5.2)对比对结果进行优化，如去除未比对上的序列，提取唯一比对上参考基因组的序列。

本发明在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。通过在indel位点附近进行局部重排比对可以大大减少由于indel而导致的很多SNP的假阳性；

本发明基于duplex-seq，并且经过了上述步骤的严谨处理，可以大大减少由于测序原因引起的误差，因此本方法最后对所有碱基的得分值赋值为ASCII字符“J”；

本发明通过降低突变位点的检测阈值，以此来检测更多的低频突变位点；对于测序深度小于2的变异位点进行过滤；结合已有的变异位点数据库dbSNP和千人基因组标准INDELs，以及体细胞突变数据库COSMIC等数据库对变异位点进行注释。

本发明方法主要开发应用于基于Duplex-seq测序数据对低频突变位点进行检测分析。本方法主要有以下几大优势：

1.数据质控系统化：本软件对于Duplex-seq的测序数据在正式分析之前，对数据质量会进行一个系统分析，如数据Q30、GC、Insert size、Duplication rate、UMI的种类、数量、比例、R1R2的平衡性等。通过各个指标来确保数据是否能够满足后续分析的要求。

2.UMI分析差异化：本软件对于单链特异性UMI和双链互补的UMI进行了差异化使用，对于构建单链一致性序列时，本软件同时利用了单链特异性UMI和双链互补的UMI，以增加UMI的多样性，可有效提高单链一致性序列的数量。另外在构建双链一致性序列时，本软件去除了单链特异性UMI，以提高双链UMI匹配的比例。

3.变异位点分析流程化：对于Duplex-seq的测序数据变异位点的分析，为了提高分析效率，本软件把变异位点锁定在目标区域内进行分析，这个区域可以根据不同的需要进行相应的调整。除此之外，对于变异位点的最大测序深度、比对得分值本软件也进行了优化。另外，对于变异位点本软件也进行了个性化的统计和注释。

4.数据结果可视化:本分析方法结果中除了提供有效的表格文件外，还生成了形象的信息图，使得更多的数据信息得以展示，也使得结果一目了然。

附图说明

图1为本发明的主要原理示意图。

图2为本发明的主要分析流程示意图。

图3为本发明数据质量控制流程图。

图4为本发明构建Duplex-tag序列文件流程图。

图5为本发明构建单链一致性序列流程图。

图6为本发明构建双链一致性序列流程图。

图7为本发明变异位点检测分析流程示意图。

图8为本发明family barcode分布情况示意图。

图9为本发明R1R2duplex barcode的均一性情况示意图。

图10为本发明碱基分布情况示意图。

图11为本发明Duplex barcode比例分布图。

图12为本发明目标区域覆盖度示例图。

图13为发明构建一致性序列统计列表示例图。

具体实施方式

为了实现本发明目的，如图1和图2所示，本发明基于duplex-seq的超低频突变位点检测分析方法，包括如下步骤：

1)对原始测序数据质量进行评估，降低数据噪声，为后续分析提供有效数据；

2)把随机barcode提取到序列文件的每一条序列的标题行，方便后续对barcode进行快速检索并创建一致性序列；

3)根据family barcode和duplex barcode创建一致性序列，排除由于建库过程或者PCR过程中引入的突变；

4)根据duplex-tag构建双链一致性序列，进一步排除序列中的非对称突变位点；

5)对比对后的数据进行局部质量矫正，并进行低频变异位点检测；将变异位点进行基因结构、功能、及临床表型三个层次的注释；

6)统计SSCS、DCS序列数目、比对结果、变异位点信息，并输出可视化图表。

具体步骤和方法如下：

1).质量控制：

1.1)利用FastQC查看测序质量结果，包括碱基得分值、GC分布、碱基平衡性、重复率等；本发明利用更加快速、高效的perl程序，统计所有序列的总Total Reads，Totalbases、Q20、Q30、GC含量、N字符等数量及相关比例。

1.2)利用AdapterRemoval^[4]去除序列3`端含有的接头序列，对于paie-end测序数据并保留配对的序列，此步骤具有较高的特异性和敏感性。

1.3)对步骤1.2)中去除接头后的序列，利用perl程序计算其重复度。

1.4)本步骤是本分析方法最为关键的步骤之一，这一步骤的相关指标决定了数据是否合格，能否进行下一步的分析。本方法利用perl程序结合shell程序提取每条序列的family barcode和duplex barcode，并利用R语言对其种类、数量、比例进行相应统计。如本实施案例中family barcode为4位碱基，则有256种组合，其比例分布比较均匀可以参与后续的分析(附图8)。对于duplex barcode本实施案例中的为双N，Spacer为2个碱基GT，可以据此提取有效的duplex barcode，其组合有16种，其分布比例也较均匀(附图9)，且其总比例达到98.04％，数据质量较好，可以进行后续分析。

2).标记序列提取：

2.1本实施案例中Duplex barcode为2个碱基，提取配对序列的duplex barcode、和4个碱基的family barcode。

2.2)对上述提取的barcode进行过滤，如有非碱基字符则丢弃。

2.3)把通过步骤2.2)过滤后的barcode构建duplex-tag，如R1 R2的duplexbarcode分别为AC、GA，family barcode为ACTA，则最后构建的duplex-tag为ACTAACGA。本步骤是本方法的核心创新之一。传统方法仅仅考虑了duplex barcode序列，使得单链的tag多样性较少。且本方法采用family barcode预先提取到单独文件的方法，使得程序读入数据仅为原始数据的四分之一，可以有效降低服务器的内存占比，提高分析效率。

2.4)根据duplex-tag分别输出配对序列(R1，R2)

2.5)对下一条序列重复上述2.1-2.4步骤。

3)构建单链一致性序列：

3.1)把步骤2)中获得的序列，利用BWA^[5]软件比对到相应的基因组获得比对文件。

3.2)获取比对文件中基因组的第一个位点。

3.3)对所有比对到步骤3.2)提取的位点上的reads进行flag字段的过滤，对于未比对上的reads(flag值为77和141)则另存为NM文件。

3.4)对所有通过步骤3.3)过滤后的reads根据duplex-tag进行排序，并进行分组。

3.5)提取步骤3.4)分组中的第一组Duplex-tag及其相关序列。

3.6)对步骤3.5)中提取的一组序列根据CIGAR string进一步分组归类，对于含有相同CIGAR string的序列则进行下一步分析，对于不含有共同CIGAR string的序列则另存为LCC文件。

3.7)对步骤3.6)中分组的序列，计算其family size，如果family size小于3则丢弃该组序列则，否则进行下一步的分析。

3.8)对步骤3.7)中通过的一组序列创建单链一致性序列。对于碱基一致性较高的位点则该位点归一为含量较高的碱基，对于一致性不够好的序列则该位点以N代替。一致性的值根据用户自行定义，本实施案例设置为70％。

3.9)通过步骤3.8)构建的单链一致性序列，过滤掉含有大量N的序列(比例大于30％)，并输出最终合格的序列到单链一致性(SSCS)文件。

3.10)创建完SSCS后，如果含有更多的duplex-tag，则重复上述步骤3.6)-3.9)，如果没有则进入步骤3.11)。

3.11)若果还有更多的位点，则重复上述步骤3.3)-3.10)，否则结束该模块分析。

4)构建双链一致性序列：

4.1)本步骤也是本分析方法的主要核心创新之一，在family barcode行使了其功能后应该将其剔除，使得后续分析更加准确。因此本实施案例中将上述模块(3)中获得的单链一致性序列比对文件由二进制的bam格式转化为可阅读的sam格式文件，并去除每条记录中4个碱基的family barcode序列，方便下一步骤更加准确的创建双链一致性序列。去除完family barcode序列后再将sam文件转换回bam文件。

4.2)提取比对文件中的第一个基因组位点信息。

4.3)提取步骤4.2)中基因组位点对应的第一个duplex tag。

4.4)寻找与步骤4.3)中的duplex-tag进行互补配对的duplex-tag，如果没有对应的duplex-tag与其进行匹配则该序列丢弃，如果有匹配的tag则进入步骤4.5)。

4.5)对含有相同duplex-tag的序列构建双链一致性序列。

4.6)对步骤4.5)创建的双链一致性序列进行过滤，如果序列中含有大量的N则该序列丢弃，本实施案例中如果N的含量超过30％，则舍弃改DCS。如果序列合格则进入步骤4.7)。

4.7)输出为双链一致性文件(DCS),并输出对应的配对序列(DCSR1R2)。

4.8)在步骤4.5)创建完双链一致性序列后，如果还有更多的duplex-tag则重复步骤4.4)-4.7)，如果没有duplex-tag则进入步骤4.9)。

4.9)如果还有更多的位点需要分析，则重复步骤4.3-4.8。否则结束分析。

5.变异分析：

5.1)将步骤4)中生存的双链一致性序列比对到基因组。

5.2)对比对结果进行优化，去除未比对上的序列，提取唯一比对上参考基因组的序列。

6)统计报告

本分析模块对数据质量、碱基平衡性、barcode的分布、构建的SSCS\DC、以及目标区域的覆盖度，变异位点频率进行了相关统计并输出可视化图表，如附图10-13。

本具体实施方式在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。通过在indel位点附近进行局部重排比对可以大大减少由于indel而导致的很多SNP的假阳性。因而本发明利用GATK的Realigner-TargetCreator来确定在INDEL附近需要进行重比对的区域，利用IndelRealigner在确定的区域内进行重新比对。

本具体实施方式基于duplex-seq，并且经过了上述步骤的严谨处理，可以大大减少由于测序原因引起的误差，因此本方法最后对所有碱基的得分值赋值为ASCII字符“J”。在传统分析过程中一般保留reads碱基质量值在Q20以上的碱基，这些碱基的错误率在～1％。因此本方法可以有效的排除测序误差，使得分析更加准确。

本具体实施方式通过降低突变位点的检测阈值，以此来检测更多的低频突变位点。

本具体实施方式对于测序深度小于2的变异位点进行过滤。传统分析方法中一般要求测序深度要达到8X以上，结果才会更加可靠，本分析方法降低了测序深度要求，更加容易检测出低频突变位点。

本具体实施方式结合已有的变异位点数据库dbSNP、千人基因组标准INDELs、体细胞突变数据库COSMIC、以及dbNSFP等数据库对变异位点进行注释。

Claims (6)

1.一种基于duplex-seq的超低频突变位点检测分析方法，其特征在于，包括如下步骤：

1)对原始测序数据质量进行评估，降低数据噪声，为后续分析提供有效数据；

2)把随机barcode提取到序列文件的每一条序列的标题行，方便后续对barcode进行快速检索并创建一致性序列；

3)根据family barcode和duplex barcode创建一致性序列，排除由于建库过程或者PCR过程中引入的突变；

4)根据duplex-tag构建双链一致性序列，进一步排除序列中的非对称突变位点；

5)对比对后的数据进行局部质量矫正，并进行低频变异位点检测；将变异位点进行基因结构、功能、及临床表型三个层次的注释；

6)统计SSCS、DCS序列数目、比对结果、变异位点信息，并输出可视化图表。

2.如权利要求1所述的一种基于duplex-seq的超低频突变位点检测分析方法，其特征在于，所述步骤1)包括如下步骤：

1.1)查看测序质量结果，包括碱基得分值、GC分布、碱基平衡性、重复率等；统计所有序列的总Total Reads，Total bases、Q20、Q30、GC含量、N字符等数量及相关比例；

1.2)去除序列中含有的接头序列；

1.3)对去除接头后的序列计算其重复度；

1.4)提取每条序列的family barcode和duplex barcode对其种类、数量、比例进行相应统计。

3.如权利要求1所述的一种基于duplex-seq的超低频突变位点检测分析方法，其特征在于，所述步骤2)包括如下步骤：

提取配对序列的Duplex barcode、family barcode；

对上述提取的barcode进行过滤，如有非碱基字符则丢弃。

把通过步骤2.2)过滤后的barcode构建duplex-tag。

2.4)根据duplex-tag分别输出配对序列R1，R2；

2.5)对下一条序列重复上述2.1-2.4步骤。

4.如权利要求1所述的一种基于duplex-seq的超低频突变位点检测分析方法，其特征在于，所述步骤3)包括如下步骤：

3.1)把步骤(2)中获得的序列，比对到相应的基因组获得比对文件；

3.2)获取比对文件中基因组的第一个位点；

3.3)对所有比对到步骤3.2)提取的位点上的reads进行flag字段的过滤，对于未比对上的reads()则另存为NM文件；其中flag字段为77和141；

3.4)对所有通过步骤3.3)过滤后的reads根据duplex-tag进行排序，并进行分组；

3.5)提取步骤3.4)分组中的第一组Duplex-tag及其相关序列；

3.6)对步骤3.5)中提取的一组序列根据CIGAR string进一步分组归类，对于含有相同CIGAR string的序列则进行下一步分析，对于不含有共同CIGARstring的序列则另存为LCC文件；

3.7)对步骤3.6)中分组的序列，计算其family size，如果family size小于3则丢弃该组序列，通过则进行下一步的分析；

3.8)对步骤3.7)中通过的一组序列创建单链一致性序列；对于碱基一致性较高的位点则该位点归一为含量较高的碱基，对于一致性不够好的序列则该位点以N代替；一致性的值根据用户自行定义，设置为70％；

3.9)通过3.8构建的单链一致性序列，过滤掉含有30％以上N的序列，并输出最终合格的序列到单链一致性(SSCS)文件；

3.10)创建完SSCS后，如果含有更多的duplex-tag，则重复上述步骤3.6)-3.9)，如果没有则进入3.11)步骤；

3.11)若果还有更多的位点，则重复上述步骤3.3)-3.10)，否则结束该模块分析。

5.如权利要求1所述的一种基于duplex-seq的超低频突变位点检测分析方法，其特征在于，所述步骤4)包括如下步骤：

4.1)将上述步骤3)中获得的单链一致性序列比对文件转化为sam格式文件，并去除family barcode序列，方便下面创建双链一致性序列；

4.2)提取比对文件中的第一个基因组位点信息；

4.3)提取步骤4.2)中基因组位点对应的第一个duplex tag；

4.4)寻找与步骤4.3)中的duplex-tag进行互补配对的duplex-tag，如果没有对应的duplex-tag与其进行匹配则该序列丢弃，如果有匹配的tag则进入步骤4.5)；

4.5)对含有相同duplex-tag的序列构建双链一致性序列；

4.6)对步骤4.5)创建的双链一致性序列进行过滤，如果序列中含有大量的N则该序列丢弃，否则进入步骤4.7)；

4.7)输出为双链一致性文件(DCS),并输出对应的配对序列(DCS R1 R2)；

4.8在步骤4.5)创建完双链一致性序列后，如果还有更多的duplex-tag则重复步骤4.4)-4.7)，如果没有duplex-tag则进入步骤4.9)；

4.9)如果还有更多的位点需要分析，则重复步骤4.3)-4.8)；否则结束分析。

6.如权利要求1所述的一种基于duplex-seq的超低频突变位点检测分析方法，其特征在于，所述步骤5)包括如下步骤：

将步骤4)中生存的双链一致性序列比对到基因组；5.2)对比对结果进行优化，如去除未比对上的序列，提取唯一比对上参考基因组的序列。