CN114093428A - 一种ctDNA超高测序深度下低丰度突变的检测系统和方法 - Google Patents

Abstract

本发明公开一种NGS双端分子标签测序技术和背景降噪算法的检测血浆ctDNA超高测序深度下低丰度突变的检测方法，包括：(1)对所有测序读段携带的双端分子标签序列收集并标注分子标签序列和组合方式；(2)标注后的读段与参考序列进行比对，将分子标签组合序列和比对位置相同的测序读段归类至单分子共识序列；(3)将存在的分子标签组互补且读段序列互补的单分子共识序列进一步归类至双链共识序列；(4)对包含单链共识序列和双链共识序列的比对结果进行突变检测，注释和过滤；(5)对突变检测结果，利用健康人检测结果使用零膨胀泊松分布算法进行背景降噪。在保证100％特异性的前提下提升了检测灵敏度的，ctDNA中丰度0.1％突变的检测敏感性达到95％以上。

Description

一种ctDNA超高测序深度下低丰度突变的检测系统和方法

技术领域

本发明涉及一种ctDNA超高深度测序深度下低丰度突变的检测系统和方法，具体地，本发明涉及对测序读段进行聚类的方法和利用健康人群检测结果进行背景降噪。

背景技术

循环肿瘤DNA(circulating tumor DNA,ctDNA)来源于凋亡、坏死的肿瘤细胞或肿瘤细胞分泌释放产生的小片段DNA，是人体循环游离DNA(circulating cell-free DNA,cfDNA)的一部分。

CN113373524A专利披露了一种ctDNA测序方法以及该方法中涉及的标签接头，从而提高检测灵敏度和特异性，实现提升检出率的同时提高对假阳性突变的分辨率，结合分子标签技术、样本标签多样化技术及超高深度测序对ctDNA进行检测，可以有效提高测序灵敏度。但是，由于在超高深度测序的过程中，分子标签仅修正了在测序过程中产生的随机误差，对于背景基线噪音，如基因组上部分特定转座子或区域受建库试剂，DNA损伤修复等混合因素影响发生的碱基偏好性突变(G>T)，仍然无法消除，所以对检测准确性仍然产生影响；并且在对于检测数据与基线样本进行比对时，由于在对照样本中绝大多数健康人在同一位点上的突变丰度为0的频率较高，数据的分布与标准泊松分布有较大偏差，导致了对比时存在一定的误差。

发明内容

本发明所要实际解决的技术问题是：在对cfDNA进行超高深度的测序时，由于基线样本数据存在的噪音、基线样本中较多的位点存在突变丰度为0导致的不符合标准泊松分布，使得数据比对后存在着检测数据存在误差。

本发明提出了一种基于双端分子标签ctDNA超高深度测序下，通过归类具有相同分子标签组合序列和健康人群结果背景降噪的方法，从而提高检测灵敏度和特异性，实现提升检出率的同时提高对假阳性突变的分辨率，最终提高检测的敏感性和准确性。对突变检测结果，利用健康人检测结果使用零膨胀泊松分布算法进行背景降噪。在保证100％特异性的前提下提升了检测灵敏度的，ctDNA中丰度0.1％突变的检测敏感性达到95％以上。

技术方案是：

一种ctDNA超高测序深度下低丰度突变的检测方法，包括如下步骤：

步骤1，对待测样本进行高通量测序并获得下机数据；

步骤2，对下机数据中的读段数据进行校验后，计算出各个突变的丰度；

步骤3，获得对照样本中各个位点上存在的突变，并在对对照样本数据集中每个位点的是否发生突变的情况进行零膨胀矫正后，计算出每个位点上的零膨胀矫正后预期平均突变丰度；

步骤4，将待测样本的突变丰度与零膨胀矫正后预期平均突变丰度进行比较，进行差异显著性判定。

所述的步骤3中，还包括：获得对照样本上的各个位点的突变人数发生百分比；并且，在步骤4中，若突变人数发生百分比小于设定阈值时，将待测样本的突变丰度与零膨胀矫正后预期平均突变丰度进行单侧检验对比，若大于设定阈值时，进行累计分布频数检验。

所述的设定阈值是1-10％。

所述的步骤1中，对待测样本进行高通量测序时，通过带有样本标签和分子标签的接头进行建库；接头包括第一接头以及第二接头；所述的第一接头以及第二接头都为Y型接头，每个Y型接头包括两条链，其中包括互补段和非互补段；在所述的互补段中包括分子标签，在非互补段中包括样本标签。

所述的步骤2中，对下机数据进行校验分析的过程包括如下步骤：S1，根据样本标签对数据进行合并；S2，对于含有相同的样本标签的数据，根据分子标签进行归类；对于含有相同的分子标签的读段，再根据读段的序列进行聚类，在一类中的读段中，只有当全部的读段的碱基完全一致时，则将该类的读段数据标记为真实数据；S3，对于S2中得到的真实数据结果，选取含有相互互补的分子标签的两类读段数据，若读段数据也在两类中相互互补，则标记为正确读段数据。

所述的第一接头以及第二接头中的样本标签相互构成一个特异性组合。

所述的样本标签的长度为6-20bp。

分子标签的长度为3-10bp。

步骤2中，还包括对数据进行注释的步骤，注释软件可以是Annover，VEP，GATKFuncotator，SnpEff软件的一种或多种联合使用。

步骤2中，各个突变是通过GATKMutact2，VarScan2，Strelka，Vardict，SomaticSniper，Samtools软件中的一种或多种联合使用而得到。

步骤2中，还包括通过比对于数据库过滤常见突变的步骤。

所述的数据库是公共数据库或者自建数据库；所述的公共数据库是如dbSNP，gnomAD，ExAC等。

一种ctDNA超高测序深度下低丰度突变的检测系统，其特征在于，包括：

测序模块，用于对待测样本进行高通量测序并获得下机数据；

数据校验模块，用于对下机数据中的读段数据进行校验；

丰度计算模块，用于对数据校验模块中得到的合格数据进行各个突变的丰度的计算；

零膨胀矫正处理模块，用于获得对照样本中各个位点上存在的突变，并在对对照样本数据集中每个位点的是否发生突变的情况进行零膨胀矫正后，计算出每个位点上的零膨胀矫正后预期平均突变丰度；

判定模块，用于将待测样本的突变丰度与零膨胀矫正后预期平均突变丰度进行比较，进行差异显著性判定。

一种计算机可读取介质，其特征在于，其记载有可以运行权利要求1-8任一项所述的检测方法的计算机程序。

有益效果

本发明通过归类具有相同分子标签的读段和健康人群背景降噪，在30000倍的超高深度测序下，提升ctDNA检测的敏感性并提高对假阳性突变的分辨率，从而整体提高ctDNA检测的灵敏度和特异性。在使用HD参考品进行性能对比的是严重，应用本发明技术可有效减小背景突变的频率；在使用Horizon商品化cfDNA参考品进行性能验证的实验中，应用本发明技术可以在0.1％丰度下，实现95％的灵敏度和100％的特异性。

附图说明

图1示出了根据本发明实施方式的分析方法原理示意图。

图2示出了根据本发明实施方式的标签接头分子结构示意图。

图3示出了根据本发明实施方式的归类单分子共识序列和双链共识序列的原理示意图。

图4示出了根据本发明实施方式的健康人群背景降噪前后对比图。

具体实施方式

本发明的数据的处理过程，可以沿用专利CN113373524A中的建库和数据处理，现简述如下：

首先，需要对数据质量控制。对原始下机数据进行质量分析，如果读段的质量指标低于阈值，例长度过短或测序碱基质量过低，则删除该读段；如果读段中包含测序接头序列，则切除测序接头序列，保留剩余部分。

接下来，需要进行测序接头的设计，测序接头包括第一接头以及第二接头；如图2所示，所述的第一接头以及第二接头都为Y型接头，每个Y型接头包括两条链，其中包括互补段和非互补段；在所述的互补段中包括分子标签，在非互补段中包括样本标签。

对于样本标签，其符合i5/i7规则，用于对每个样本来源的读段进行标记和校验，在下机数据分析过程中，通过对i5-i7双端标签进行分析，由于一个样本中的配对是特定的，因此可以将样本标记并且可以去除标签跳跃引起的错误。

对于分子标签，即分子标签barcode，该分子标签不仅能够识别单分子错误突变，也能够识别双链DNA错误突变，含有4bp的单分子和双链标签序列。通过将下机读段中的数据按照分子标签进行归类，再进一步在同一类下按照读段序列进行聚类后，获得小类，每一个小类中完全一致的读段数据才被认为是单链共识序列。进一步地，由于分子标签是在互补段上，因此将上述步骤中获得的读段进行双链还原后，只有在互补的两条链上的碱基完全一致时，才被认为是正确的读段数据，就获得了双链共识序列。

在获得了上述的校验后的读段数据后，与参考基因组进行对比时，仍然会存在着由基线噪音导致的数据判定准确性方面的问题，如建库试剂引起的热点区域碱基偏好性突变。本专利的主要改进方面在于对背景基线噪音的处理，本方法中通过健康人群的下机数据，对每个位点上存在的突变情况进行检测，考虑到绝大多数健康人在同一位点上的突变丰度为0的频率较高，数据的分布与标准泊松分布有较大偏差，无法使用普通的泊松分布解释，因此在对检测结果进行矫正时，采用零膨胀柏松回归模型消除高频率健康人突变丰度为0的影响，获得位点零膨胀矫正后突变丰度的平均值和方差；通过将测试样本的突变丰度与这个基线值对比后，可以判定出测试样本的突变情况，该方法解决了基线噪音对测试结果的干扰作用，避免了在超高深度测序条件下的背景基线噪音对结果的影响。

实施例1

(1)数据质量控制模块

对原始下机数据进行质量分析，如果读段的质量指标低于阈值，例长度过短或测序碱基质量过低，则删除该读段；如果读段中包含测序接头序列，则切除测序接头序列，保留剩余部分。

(2)分子标签识别标记模块

对每条读段由5‘端到3’端开始识别两端分子，将5‘端和3’端分子标签序列分别标记为UMI1和UMI2，将分子标签组合形式UMI1-UMI2记录在FASTQ文件中读段识别码中，在碱基序列和测序质量分数中切除分子标签对应序列和质量结果。

(3)读段比对和单双链共识序列归类模块

a.对标记后的所有读段，使用BWA MEM比对到人类参考基因组(hs37d5)上，并使用picard按照序列名对比，对数据进行排序。

b.对a得到的排序比对数据进行共识序列合并，对于某一个样本中得到的读段数据，根据其读段上的分子标签组合和参考基因组上的比对位置进行进一步分类，在每一个分子标签组合下，再进一步地根据读段碱基情况聚类，获得一个小类，在一个小类中，要求所有序列至少大于等于2条；如果在一个小类中读段之间的数据是完全一致的，则认为该碱基突变为真实的单链DNA突变，则进入步骤a中继续分析处理；当只要有一个读段与其它的读段数据不完全一致时，则将这个小类的数据全部标记，并不进入后续步骤；该步骤示意图如图3所示。

c.对于b得到的分类数据，当存在反向互补的两条正义和反义链时，对两条链序列构建双链共识序列，碱基一致的序列仍保持原碱基，认为该碱基突变为真实的原双链DNA突变；双链碱基不一致时，则进行标记后，不进入后续的步骤中使用。

(4)突变检测和注释过滤模块

分别对构建好的单双链共识序列数据进行合并并排序，对数据进行突变分析，然后使用注释软件进行注释。对比公共数据库和自建数据库，过滤掉人群常见胚系突变和流程引入常见假突变。公共数据库，如dbSNP，gnomAD，ExAC等，收录了大量胚系突变信息。自建数据库，是使用了100例健康人血液样本，用既定相同流程进行突变检测，收集3次以上出现的突变结果；这些突变中除去胚系突变外，包含由于软件版本引进的习惯性假突变。

为进一步修正cfDNA测序中比对引起的碱基偏好背景影响，使用相同流程分析100例健康人样本，获得100个健康BAM文件。对每一个BAM文件，使用samtools mpileup统计捕获探针覆盖范围内每一个位点的信息，统计内容如下表格：

使用R包pscl的zeroinfl功能的默认参数(https://github.com/atahk/pscl)对捕获探针覆盖范围内的每一个位点的100个健康人突变结果进行零膨胀泊松分布分析，获得每一个位点的零膨胀矫正后预测突变丰度的平均值与标准差的健康集合背景参照，结果汇总格式如下表格：

列名	描述
		CHR	参考基因组染色体号
POS	参考基因组上的位置
		REF	参考基因组上该位置的碱基
ALT	该位置的突变碱基
		TotalSamples	健康组人数
ZIPFracSamples	零膨胀矫正后健康组突变人数发生百分比
		ZIPmeanAF	零膨胀矫正后预期平均突变丰度
ZIPsdAF	零膨胀矫正后预期突变丰度标准差

表中的突变人数发生百分比是指：对于某一个位点，在全部的对照样本中发生了突变的样本数量的占比。

对ctDNA的检测结果中的所有位点的突变丰度AF，与健康集合背景参照进行检验：(1)ZIPFracSamples<＝5％时，说明样本量分布不足，使用one-sided z test验证；(2)ZIPFracSamples>5％时，使用累计分布函数(cumulative Distribution Function，CDF)检验；(3)保留验证结果p-value<0.05，保留该突变结果。

下面将结合实例进一步说明本发明的有益效果。

本发明的实现步骤主要包括：

(1)单双链共识序列配合背景降的检测方法对背景突变的影响

选取四例HD标准品(HD776,HD777,HD778,HD779),突变丰度分别为0％，5％，1％，0.1％，使用不同的检测方法进行突变检测，结果如图4可知，未构建共识序列的样本有较高的背景突变，构建共识序列数据和背景降噪均有效减少假阳突变发生频率，同时使用效果最佳，证明本发明的检测方法确实有效。

(2)Horizon商品化cfDNA标准品验证

选取了不同丰度的6个基因突变位点标准品，同时设立空白对照组，6个基因突变位点分别为NRAS Q61K、NRAS A59T、PIK3CA E545K、EGFR T790M、EGFR L858R和KRAS G12D，对应的丰度依次为0.5％，0.2％，0.1％和0％空白对照。每个位点每个稀释浓度分别设置4个重复，在不出现假阳性的前提下设定阳性判断值，分别评估本发明的检测方法和常规的检测方法在不同丰度对应突变位点的检出情况，具体检测结果见下表：

根据本实施例中的检测结果，在不出现假阳性的前提下：对照组实验在0.5％、0.2％和0.1％的标准品检测中灵敏度分别为91.7％、79.2％和41.7％，而使用本发明分析方法后，在0.5％、0.2％和0.1％的标准品检测中灵敏度分别高达100％、100％和95.8％。

(3)不同细胞系混合进一步验证阳性预测值

上步验证使用的是商品化样本，位点数量较为有限，无法评估大量位点的检测阳性预测值，进一步采用临床上常用的两个NA细胞系进行梯度稀释，验证本发明标签接头+背景降噪技术的检测性能：分别选择NA18535和NA19240作为背景细胞系和检测细胞系，按照0.5％、0.2％和0.1％三个梯度进行稀释，即NA18535与NA19240的DNA量配比分别为99.9:0.1、99.8:0.2和99.5:0.5。分别采用了三种突变检测技术进行对比：

第一组：未采用的未进行专利CN113373524A中双链碱基互补分析的检测；

第二组：采用了专利CN113373524A中的标签接头以及双链碱基互补分析的检测；

第三组：在第二组的基础上还增加了对基线零膨胀矫正的处理；

对上述混合后的样本进行高通量测序，通过生信分析获得两个NA样本差异的纯合和杂合SNP位点，考察三种技术对于杂合和纯合SNP位点的检测阳性预测值(PPV＝TP/(TP+FP))，具体检测结果见下表：

根据本实施例中的检测结果，随着检测丰度的逐渐降低，阳性预测值大幅度下降，至0.1％丰度时，阳性预测值只有64％左右；标签接头技术修正了测序流程中的随机误差，即使是0.1％丰度时，对于SNP的检测阳性预测值可提高至88％；而使用本发明标签接头+背景降噪技术，进一步修正背景碱基偏好误差，在0.1％丰度时，进一步提高阳性预测值至95％。

Claims (10)

1.一种ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，包括如下步骤：

步骤1，对待测样本进行高通量测序并获得下机数据；

步骤2，对下机数据中的读段数据进行校验后，计算出各个突变的丰度；

步骤3，获得对照样本中各个位点上存在的突变，并在对对照样本数据集中每个位点的是否发生突变的情况进行零膨胀矫正后，计算出每个位点上的零膨胀矫正后预期平均突变丰度；

步骤4，将待测样本的突变丰度与零膨胀矫正后预期平均突变丰度进行比较，进行差异显著性判定。

2.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，所述的步骤3中，还包括：获得对照样本上的各个位点的突变人数发生百分比；并且，在步骤4中，若突变人数发生百分比小于设定阈值时，将待测样本的突变丰度与零膨胀矫正后预期平均突变丰度进行单侧检验对比，若大于设定阈值时，进行累计分布频数检验。

3.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，所述的步骤1中，对待测样本进行高通量测序时，通过带有样本标签和分子标签的接头进行建库；接头包括第一接头以及第二接头；所述的第一接头以及第二接头都为Y型接头，每个Y型接头包括两条链，其中包括互补段和非互补段；在所述的互补段中包括分子标签，在非互补段中包括样本标签。

4.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，所述的步骤2中，对下机数据进行校验分析的过程包括如下步骤：S1，根据样本标签对数据进行合并；S2，对于含有相同的样本标签的数据，根据分子标签进行归类；对于含有相同的分子标签的读段，再根据读段的序列进行聚类，在一类中的读段中，只有当全部的读段的碱基完全一致时，则将该类的读段数据标记为真实数据；S3，对于S2中得到的真实数据结果，选取含有相互互补的分子标签的两类读段数据，若读段数据也在两类中相互互补，则标记为正确读段数据。

5.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，所述的第一接头以及第二接头中的样本标签相互构成一个特异性组合。

6.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，所述的样本标签的长度为6-20bp；分子标签的长度为3-10bp。

7.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，步骤2中，还包括对数据进行注释的步骤，注释软件可以是Annover，VEP，GATKFuncotator，SnpEff软件的一种或多种联合使用。

8.根据权利要求1所述的ctDNA超高测序深度下低丰度突变的检测方法，其特征在于，在一个实施方式中，步骤2中，各个突变是通过GATKMutact2，VarScan2，Strelka，Vardict，SomaticSniper，Samtools软件中的一种或多种联合使用而得到；在一个实施方式中，步骤2中，还包括通过比对于数据库过滤常见突变的步骤；在一个实施方式中，所述的数据库是公共数据库或者自建数据库；所述的公共数据库是如dbSNP，gnomAD，ExAC等。

9.一种ctDNA超高测序深度下低丰度突变的检测系统，其特征在于，包括：

测序模块，用于对待测样本进行高通量测序并获得下机数据；

数据校验模块，用于对下机数据中的读段数据进行校验；

丰度计算模块，用于对数据校验模块中得到的合格数据进行各个突变的丰度的计算；

零膨胀矫正处理模块，用于获得对照样本中各个位点上存在的突变，并在对对照样本数据集中每个位点的是否发生突变的情况进行零膨胀矫正后，计算出每个位点上的零膨胀矫正后预期平均突变丰度；

判定模块，用于将待测样本的突变丰度与零膨胀矫正后预期平均突变丰度进行比较，进行差异显著性判定。

10.一种计算机可读取介质，其特征在于，其记载有可以运行权利要求1-8任一项所述的检测方法的计算机程序。

Images