CN109887548A

CN109887548A - 基于捕获测序的ctDNA占比的检测方法及检测装置

Info

Publication number: CN109887548A
Application number: CN201910049677.7A
Authority: CN
Inventors: 韩天澄; 于佳宁; 侯军艳; 林小静; 陈维之; 杜波; 何骥
Original assignee: Yue Yue Biotechnology Jiangsu Co Ltd
Current assignee: Yue Yue Biotechnology Jiangsu Co Ltd
Priority date: 2019-01-18
Filing date: 2019-01-18
Publication date: 2019-06-14
Anticipated expiration: 2039-01-18
Also published as: CN109887548B

Abstract

本发明公开了一种基于捕获测序的ctDNA占比的检测方法及检测装置。该检测方法包括以下步骤：S1，获取基线样本DNA和待测cfDNA的捕获测序的基因数据；S2，同时使用基线样本DNA中的纯合与杂合位点，挑选出基线样本DNA和待测cfDNA之间突变频率存在显著不同且满足预设过滤条件的位点作为候选SNP位点；S3，结合候选SNP位点所在区域拷贝数突变的情况判断正常细胞DNA与ctDNA的候选SNP位点的突变类型；以及S4，使用极大似然法建立概率模型，通过概率模型计算得到ctDNA占比。应用本发明的技术方案能够同时从多个方面提高血浆ctDNA的CNV检测灵敏性和准确性。

Description

基于捕获测序的ctDNA占比的检测方法及检测装置

技术领域

本发明涉及基因工程技术领域，具体而言，涉及一种基于捕获测序的ctDNA占比的检测方法及检测装置。

背景技术

循环肿瘤DNA(circulating tumor ctDNA)是肿瘤细胞凋亡、坏死或分泌产生的，是循环游离DNA(circulating cell-free DNA，cfDNA)中的一种。ctDNA在血液中的半衰期短，可用于肿瘤患者的实时监控。除ctDNA的单核苷酸多态性(Single NucleotidePolymorphisms，SNP)、插入缺失标记(insertion-deletion，InDel)与拷贝数变异(copynumber variation，CNV)外，ctDNA在cfDNA中的占比，同样可以作为肿瘤进展、预后的一项指标。根据ctDNA的占比，也可以对SNP、InDel与CNV的检测结果进行校正。

ctDNA无法像肿瘤组织一样，采用病理切片染色的方式进行肿瘤细胞占比的检测，主要的检测方法为二代基因测序(Next Generation Sequencing，NGS)等。考虑到ctDNA的特性，一般采取捕获测序的方式：通过打断DNA、PCR扩增、捕获DNA片段与荧光测序等一系列流程获得海量的DNA片段信息，再回贴获得待测样品在捕获区域的DNA序列并应用软件或算法进行进一步分析。目前，针对二代测序数据的肿瘤细胞占比检测软件的检测主要利用分为两步：第一步是获取待测样本与参考样本在各位点的深度与碱基信息，提取可用的SNP，并使用在各目标捕获区间的平均测序深度来量化各捕获区域的拷贝数；第二步是利用第一步获得的信息，通过算法对肿瘤细胞的占比进行估计。常用的算法包括：Sequenza、FACETS、PureCN等。

其中，Sequenza使用Python和R语言，针对全外显子组测序(Whole ExomeSequencing，WES)设计，不适用于捕获测序。要求提供的输入文件包括：待测样本与配对的参考样本由SAMtools产生的pileup文件与人类基因组参考文件。参考样本必须为待测样本的配对正常细胞样本。其根据每个区域内位点的平均测序深度进行GC含量校正，计算拷贝数，结合胚系杂合突变的BAF(B Allele Frequency)建立贝叶斯概率模型，对肿瘤细胞占比进行估计，并输出根据肿瘤细胞占比校正后的拷贝数变异检测结果。

FACETS的预处理模块使用perl与c++语言编写，同时调用了SAMtools软件；统计分析模块采用R语言进行编写，要求的输入文件为待测样本与配对的参考样本测序数据经过比对后生成的BAM文件。拷贝数的计算与Sequenza类似。但在使用SNP信息时，对所有dbSNP与1000genome数据库有记载的SNP计算了待测组织相对于配对正常组织的log-ratio，并对杂合SNP额外计算了log-odds-ratio，以校正BAF的偏差。结合BAF与前一步获得的拷贝数信息，建立概率模型。最终输出结果为：待测样本校正后的拷贝数变异结果，以及肿瘤细胞占比、染色体倍性、肿瘤异质性检测结果与图形化展示。

PureCN使用R语言编写，可选择使用配对的正常样本或一组非配对的正常样本作为参考样本。当使用非配对样本时，需要利用dbSNP与COSMIC公共数据库。要求的输入文件为待测样本与参考样本的比对后BAM文件、MuTect或其他SNP检测软件输出的成对/单个样本的突变列表VCF文件。拷贝数计算过程与Sequenza和FACETS类似。当获得拷贝数检测结果后，结合胚系杂合突变在待测样本中的突变情况，采用极大似然法建立模型。该软件的输出结果包括：待测样本校正后的拷贝数变异结果，以及肿瘤细胞占比、染色体倍性的检测结果与图形化展示。

上述方法均存在一定的技术缺陷，例如：病理切片染色的分辨率低，无法获得确切的肿瘤细胞占比，仅可给出大致的范围。该方法仅能够用于肿瘤组织中肿瘤细胞占比的检测，不适用于ctDNA。捕获测序分辨率较高，但分析过程较为复杂。目前针对捕获数据的分析软件，部分仅适用于全外显子测序。适用于非全外显子捕获数据的软件各自也存在一些缺陷，对肿瘤细胞占比较低样本，灵敏度较低，且都不是针对ctDNA设计的。

Sequenza：针对全外显子测序数据设计，不适用于捕获测序。

FACETS：无法检测低浓度存在于血浆或其他类型体液中的ctDNA。

PureCN：没有考虑到肿瘤异质性，假定仅存在一种克隆，对DNA片段来源于不同组织的cfDNA样本，估计不准确。染色体倍性的估计不准，故而影响了肿瘤细胞占比的准确性。此外，与FACETS一样，无法检测低浓度存在于血浆或其他类型体液中的ctDNA。

发明内容

本发明旨在提供一种基于捕获测序的ctDNA占比的检测方法及检测装置，以解决现有技术中ctDNA占比检测存在的准确度低的技术问题。

为了实现上述目的，根据本发明的一个方面，提供了一种基于捕获测序的ctDNA占比的检测方法。该检测方法包括以下步骤：S1，获取基线样本DNA和待测cfDNA的捕获测序的基因数据，进行基因数据的处理，获得基线样本DNA注释后突变列表、待测cfDNA注释后突变列表和待测cfDNA拷贝数变异列表；S2，同时使用基线样本DNA中的纯合与杂合位点，对待测cfDNA中的每个突变进行统计检验，挑选出基线样本DNA和待测cfDNA之间突变频率存在显著不同且满足预设过滤条件的位点作为候选SNP位点；S3，结合候选SNP位点所在区域拷贝数突变的情况判断正常细胞DNA与ctDNA的候选SNP位点的突变类型；以及S4，利用候选SNP位点不同碱基的支持reads数量与拷贝数，使用极大似然法建立概率模型，通过概率模型计算得到ctDNA占比。

进一步地，概率模型的评估指标包括：对每一个候选SNP位点进行ctDNA的占比估计、对ctDNA来源的异同进行判断以及对待测cfDNA整体的ctDNA占比情况进行检测。

进一步地，在概率模型计算得到ctDNA占比后通过可视化模块绘制散点图和概率密度分布图对结果进行展示。

进一步地，S3中所用的判断标准如下表1所示：

表1

其中，突变类型中的A代表参考碱基类型，B代表突变碱基类型，待测cfDNA拷贝数缺失的判断标准为：

上式中CN_i是样本第i个候选SNP位点的拷贝数，CCF_i是样本第i个候选SNP位点对应的ctDNA占比，即根据最终预测的ctDNA占比的变化，待测cfDNA拷贝数缺失的标准也随之变化。

进一步地，针对候选SNP位点的突变类型的各种可能性，待测cfDNA在给定CCF_i时的理论突变频率按照以下表2中的公式进行计算：

表2

候选SNP位点的总测序深度、突变碱基支持reads数与突变频率满足如下公式所代表的分布：

AD_i～B(DP_i，VAF_i)

其中，AD_i是样本第i个候选SNP位点的突变碱基支持reads数；DP_i是样本第i个候选SNP位点的总测序深度；VAF_i是样本第i个候选SNP位点的突变频率；对于某一给定CCF_i，代入分布与理论突变频率的公式，计算各位点在该CCF_i下获得相应AD_i与DP_i的似然值；对每一个候选SNP位点使用极大似然法，使位点i似然值最大化的CCF_i占比即为该位点的CCF_i预测值，具体的log极大似然公式为：

进一步地，检测方法还包括：对基线样本DNA的候选SNP位点为纯合且待测cfDNA拷贝数无缺失的位点的突变类型和CCF_i进行预测，具体包括：对基线样本DNA的候选SNP位点为纯合且待测cfDNA拷贝数无缺失的位点外的其他位点的CCF进行原假设为单峰分布的Hartigans单峰检验：若检验接受原假设，认为不存在多峰分布，则假定这部分位点均来自于同一肿瘤组织来源，对所有N个候选SNP位点的CCF进行极大似然估计，作为CCF_temp：

若检验拒绝原假设，认为存在多峰分布，则对这些位点的CCF_i进行聚类，假定聚到同一簇的候选SNP位点都来自同一肿瘤组织来源，且具有同一CCF_j；聚类的簇数J由gap法确认；完成聚类后，对每一簇的N_j个候选SNP位点的CCF_j分别进行极大似然估计，作为CCF_{j_temp}：

完成其他位点的肿瘤组织来源聚类与CCF_{j_temp}的计算后，对未能确定突变类型的候选SNP位点，分别计算其作为纯合突变、杂合突变且同源于已知聚类簇的似然值，并以最大化似然值的标准对候选SNP位点进行分类：

完成上述过程后，所有候选SNP位点的突变类型与来源均得到确定；加入新确定的突变类型与来源的候选SNP位点，对每一来源簇的候选SNP位点重复上述的极大似然估计过程，获得ctDNA来源数量J、分来源的占比预测结果CCF_cluster_j。

进一步地，检测方法还包括对占比预测结果进行检验的步骤，具体包括：如待测cfDNA包含多个肿瘤组织来源，占比最高的肿瘤组织来源，包含的候选SNP位点数量低于2个或低于总候选SNP位点数量的10％，则认为该组织来源的估计结果不可靠，剔除该估计结果；如待测cfDNA最终的最高占比估计结果低于15％，或候选SNP位点数量过少，或基线样本DNA纯合位点与杂合位点最终的估计结果存在明显差异，则认为待测cfDNA存在较高的低ctDNA占比可能，剔除所有基线样本DNA的杂合位点，并放宽杂合位点的过滤条件，重新提取候选SNP位点的列表并进行估计；对于部分待测cfDNA，重新提取候选SNP位点列表后仍未获得足够可用于预测的候选SNP位点时，认为待测cfDNA的ctDNA占比极低，极低为低于2％。

进一步地，S2中，对待测cfDNA中的每个SNP位点的突变进行统计检验包括：

SNP位点的总测序深度、突变碱基支持reads数与突变频率满足如下公式所代表的分布：

AD_i～B(DP_i，VAF_i)

其中，AD_i是样本第i个SNP位点的突变碱基支持reads数；DP_i是样本第i个SNP位点的总测序深度；VAF_i是样本第i个SNP位点的突变频率；上式表示AD_i服从n＝DP_i，p＝VAF_i的二项分布；

根据上述概率分布，分别对每个SNP位点使用二项检验进行如下的验证：

1)在基线样本DNA中是否为纯合，基线样本突变频率＝0％，

假定背景误差为0.01％，检验基线样本DNA；

原假设：VAF_i＜0.0001；

备择假设：VAF_i≥0.0001；

如上述检验p-value大于0.05，接受原假设，基线样本DNA在该位点为纯合；

2)在基线样本DNA中是否为杂合，基线样本突变频率＝50％；

杂合突变频率为50％；

原假设：VAF_i＝0.5；

备择假设：VAF_i≠0.5；

如上述检验p-value大于0.05，接受原假设，基线样本DNA在该位点为杂合；

3)在待测cfDNA中是否为突变，频率是否高于背景误差与测序错误，假定背景误差为0.01％，待测cfDNA的样本测序错误为该样本所有突变频率低于1％的SNP位点突变频率的中位数；

原假设：VAF_i＜max{0.0001，median(all SNP V AF)}；

备择假设：V AF_i≥max{0.0001，median(all SNP V AF)}；

如上述检验p-value小于0.05，拒绝原假设，待测cfDNA在该位点为突变；

4)在待测cfDNA中是否为杂合，基线样本DNA突变频率＝50％，

杂合突变频率为50％；

原假设：VAF_i＝0.5；

备择假设：VAF_i≠0.5；

如上述检验p-value大于0.05，接受原假设，待测cfDNA在在该位点为杂合；

另外，使用Fisher精确检验，检验待测cfDNA中的突变频率是否与基线样本DNA存在显著差异，待测cfDNA突变频率≠基线样本DNA突变频率：

原假设：VAF_cfDNA_i＝VAF_baseline_i；

备择假设：VAF_cfDNA_i≠VAF_baseline_i。

如上述检验p-value小于0.05，拒绝原假设，待测cfDNA的突变频率与基线DNA样本存在显著差异；

优选的，S2中预设过滤条件如下：

读取算法中内置的多态性基因列表，如待过滤SNP总表中的位点所在基因存在多态性，剔除该SNP；

如位点在待测cfDNA与基线样本DNA中深度低于50x，剔除该SNP；

如位点所在基因的拷贝数高于2.2或低于1，剔除该SNP；

如位点位于InDel列表中任意一个InDel所在位置的上游或下游50bp内，剔除该SNP；

如位点在基线样本DNA中为纯合，要求待测cfDNA突变频率高于基线样本DNA，并且在检验中与背景噪音或基线样本DNA存在明显差异；

如位点在基线样本DNA中为杂合，要求待测cfDNA非杂合，且突变频率与待测cfDNA存在明显差异。

根据本发明的另一个方面，提供一种基于捕获测序的ctDNA占比的检测装置。该装置用于存储或者运行模块，或者模块为装置的组成部分；其中，模块为软件模块，软件模块为一个或多个，软件模块用于执行上述任一种检测方法。

进一步地，软件模块包括：待测cfDNA SNP过滤模块、待测cfDNA的ctDNA占比估计模块和ctDNA占比检测结果可视化模块；或软件模块包括整合了待测cfDNASNP过滤模块、待测cfDNA的ctDNA占比估计模块和ctDNA占比检测结果可视化模块的自动化模块。

应用本发明的技术方案，同时使用待测样本在基线样本中纯合与杂合位点上的突变，增加了可使用的突变位点信息；通过统计检验方法对基线样本与待测样本突变的纯合、杂合与拷贝数变异情况进行检验与过滤，保证最终用于估计占比的位点的可靠性；进一步的，本发明考虑ctDNA来自于多个肿瘤组织来源的可能性，并对此进行验证，同时从多个方面提高血浆ctDNA的CNV检测灵敏性和准确性。

更进一步的，本发明对于低肿瘤细胞比例的特殊样本，尤其是血浆ctDNA样本的捕获测序数据做了有针对性SNP位点筛选，可以更灵敏更准确地用于低纯度组织样本或者血浆样本的肿瘤细胞或ctDNA占比。通过ctDNA梯度稀释细胞系测试实验验证，本发明可以准确的检测混样样本中的ctDNA占比，与稀释梯度具有极高的线性相关性，且能够区分不同的ctDNA来源，给出各来源的占比。

附图说明

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

图1示出了本发明一实施方式的基于捕获测序的ctDNA占比的检测方法流程示意图；以及

图2示出了实施例1中使用plot模块，对sample F-1与sample F-2的鉴定结果展示图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。

本发明中涉及的缩写或术语解释如下：

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

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

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

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

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

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

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

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

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

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

VarScan：一种处理高通量测序数据的工具，可用于鉴定SNP与InDel。

Annovar：一种针对SNP或InDel的注释工具。

cnvkit：一种处理高通量测序数据的工具，可用于鉴定拷贝数变异。

cfDNA梯度稀释样本：由混合血细胞DNA样本与包含ctDNA的血浆cfDNA样本所获得，通过调整样本的混合比例，以获得具有不同cfDNA占比与ctDNA占比的样本组。

基线样本：仅包含正常DNA的样本，通常为血细胞DNA样本或正常组织DNA样本。本发明实例中为血细胞DNA样本。

针对背景技术中提到的技术问题，本发明的目的在于提供以下适用于ctDNA检测的方法或装置，克服现有软件存在的灵敏度低、准确度低、稳定性差等缺陷，为ctDNA占比的检测提供一种可靠的检测方法，并提供易于解读的结果展示方案。

本发明的主要技术原理如下：基于捕获测序能够在发生体细胞突变的位置在cfDNA与基线DNA中捕获到不同比例的突变型与野生型支持reads。首先对突变位点进行过滤，提取出配对样本之间突变频率存在显著不同的位点作为候选SNP位点，结合位点所在区域拷贝数突变的情况，确定正常细胞DNA与ctDNA的突变类型。然后利用上述候选SNP位点不同碱基的支持reads数量与拷贝数，使用极大似然法建立概率模型。建立模型的同时，考虑ctDNA来自多个肿瘤组织来源的可能性，最终输出ctDNA的占比检测结果、ctDNA样本中以突变为单位的ctDNA占比与结果的图形化展示。

根据本发明一种典型的实施方式，提供一种基于捕获测序的ctDNA占比的检测方法。该检测方法包括以下步骤：S1，获取基线样本DNA和待测cfDNA的捕获测序的基因数据，进行基因数据的处理，获得基线样本DNA注释后突变列表、待测cfDNA注释后突变列表和待测cfDNA拷贝数变异列表；S2，同时使用基线样本DNA中的纯合与杂合位点，对待测cfDNA中的每个突变进行统计检验，挑选出基线样本DNA和待测cfDNA之间突变频率存在显著不同且满足预设过滤条件的位点作为候选SNP位点；S3，结合候选SNP位点所在区域拷贝数突变的情况判断正常细胞DNA与ctDNA的候选SNP位点的突变类型；以及S4，利用候选SNP位点不同碱基的支持reads数量与拷贝数，使用极大似然法建立概率模型，通过概率模型计算得到ctDNA占比。其中“预设过滤条件”可以根据实际需求进行设定。

优选的，概率模型的评估指标包括：对每一个候选SNP位点进行ctDNA的占比估计、对ctDNA来源的异同进行判断以及对待测cfDNA整体的ctDNA占比情况进行检测。

根据本发明一种典型的实施方式，在概率模型计算得到ctDNA占比后通过可视化模块绘制散点图和概率密度分布图对结果进行展示。在本发明一实施方式中，基于捕获测序的ctDNA占比的检测方法的主流程如图1所示，主要包括以下几个步骤：1.获取基线样本DNA(也可称基线DNA)和待测cfDNA(可以是血浆cfDNA)的捕获测序fastq文件，利用基因组比对工具bwa mem进行序列比对，利用SAMtools和Picard工具对比对后的bam文件进行排序和标记重复处理；2.利用SAMtools获取处理后bam文献按位置整理的read信息，利用VarScan与Annovar对整理出的mpileup文件进行SNP与InDel的检测与注释，获得基线样本DNA和待测cfDNA样本的注释后突变列表；3.利用cnvkit检测待测cfDNA的拷贝数，获得拷贝数变异列表；4.利用本发明的SNP过滤模块，整理基线样本DNA和待测cfDNA注释后的突变列表与待测cfDNA的拷贝数变异列表，并挑选出可用的突变位点，输出整理后的候选SNP位点列表list文件；5.利用本发明的ctDNA占比估计模块，分析候选SNP位点，输出ctDNA的占比检测结果，与ctDNA中各候选SNP位点的来源与占比情况；6.利用本发明的检测结果可视化模块，根据检测结果绘制散点图与概率密度分布图对结果进行展示。

典型的，本发明所要求的输入文件包括：待测样本(待测cfDNA)与基线样本DNA经过SNP与InDel检测、注释后的突变列表文件(txt格式，要求必须包含每个突变所在的染色体、位置、参考碱基类型、突变碱基类型、位点测序深度、突变支持reads数量、基因名称，其余可选)、待测cfDNA经过cnvkit分析后获得的拷贝数变异列表(cns格式，包含染色体、拷贝数变异起始点、终止点、区间内基因名称、拷贝数变异log2ratio、区间平均深度、区间内探针设计数量、权重八列信息)。

典型的，本发明的输出文件包括：候选突变列表(包含候选突变在基线样本与待测样本中的深度、突变支持reads数量、CNV结果)、ctDNA检测结果文件(包含待测样本ctDNA占比检测结果)、候选SNP位点检测结果文件(包含待测样本(待测cfDNA)部分体细胞突的克隆与占比检测结果)。

在本发明一典型的实施方式中，基于捕获测序的ctDNA占比的检测装置主要包括以下几个模块：

待测样本SNP过滤模块(filter模块)：该模块要求输入输出文件前缀、输出路径、对SNP最低突变频率的要求、待测样本与基线样本注释后的突变列表txt文件、待测样本经过cnvkit分析后获得的拷贝数变异列表cns路径。

本模块首先根据输入的待测样本突变列表，整理出待测样本(待测cfDNA)所有的SNP突变所在的染色体、位置、参考碱基类型、突变碱基类型、突变所在基因名称。再从输入的基线样本突变列表中提取每个突变在基线样本中相应位点的总测序深度、突变碱基支持reads数、突变频率，并从输入的拷贝数变异列表中提取每个突变在待测样本中所属基因的拷贝数，生成包含上述所有信息的待过滤SNP总表。

除上述待过滤SNP总表外，另外提取出输入的待测样本与基线样本DNA突变列表中所有的InDel，将其对应的染色体、位置进行整理，生成InDel列表。

完成上述步骤后，针对SNP总表中的位点进行统计检验。

AD_i～B(DP_i，VAF_i)

其中，AD_i是样本第i个SNP位点的突变碱基支持reads数；DP_i是样本第i个SNP位点的总测序深度；VAF_i是样本第i个SNP位点的突变频率；上式表示AD_i服从n＝DP_i，p＝VAF_i的二项分布。

(1)在基线样本DNA中是否为纯合(基线样本突变频率＝0％)

假定背景误差为0.01％，检验基线样本DNA。

原假设：VAF_i＜0.0001；

备择假设：VAF_i≥0.0001。

如上述检验p-value大于0.05，接受原假设，基线样本在该位点为纯合。

(2)在基线样本中是否为杂合(基线样本突变频率＝50％)

杂合突变频率为50％。

原假设：VAF_i＝0.5；

备择假设：VAF_i≠0.5。

如上述检验p-value大于0.05，接受原假设，基线样本在该位点为杂合。

(3)在待测样本中是否为突变而非背景噪音或样本测序错误(待测样本突变频率高于背景噪音/样本测序错误)

假定背景误差为0.01％，待测样本的样本测序错误为该样本所有突变频率低于1％的SNP位点突变频率的中位数。

原假设：VAF_i＜max{0.0001，median(all SNP V AF)}；

备择假设：VAF_i≥max{0.0001，median(all SNP V AF)}。

如上述检验p-value小于0.05，拒绝原假设，待测样本在该位点为突变。

(4)在待测样本中是否为杂合(基线样本突变频率＝50％)

杂合突变频率为50％。

原假设：VAF_i＝0.5；

备择假设：VAF_i≠0.5。

如上述检验p-value大于0.05，接受原假设，待测样本在该位点为杂合。

另外使用Fisher精确检验，检验待测样本中的突变频率是否与基线样本存在显著差异(待测样本突变频率≠基线样本突变频率)：

原假设：VAF_cfDNA_i＝VAF_baselinei；

备择假设：VAF_cfDNA_i≠VAF_baseline_i。

如上述检验p-value小于0.05，拒绝原假设，待测样本中的突变频率与基线样本存在显著差异。

过滤条件共计6组：

(1)读取算法中内置的多态性基因列表，如待过滤SNP总表中的位点所在基因存在多态性，剔除该SNP；

(2)如位点在待测样本与基线样本中深度低于50x，剔除该SNP；

(3)如位点所在基因的拷贝数高于2.2或低于1，剔除该SNP；

(4)如位点位于InDel列表中任意一个InDel所在位置的上游或下游50bp内，剔除该SNP；

(5)如位点在基线样本中为纯合，要求待测样本突变频率高于基线样本，并且在检验中与背景噪音或基线样本存在明显差异；

(6)如位点在基线样本中为杂合，要求待测样本非杂合，且突变频率与待测样本存在明显差异。

完成上述过滤步骤后，输出过滤后的候选SNP位点总表。输出路径为给定的输出路径，输出文件前缀由输入参数制定；

待测样本ctDNA占比估计模块(estimate模块)：该模块要求输入输出文件前缀、输出路径、待测样本的候选突变列表list文件路径。

该模块根据过滤模块输出的候选SNP位点总表，首先根据基线样本的纯合/杂合情况、基线样本与待测样本的突变频率、待测样本的拷贝数，对位点在基线样本DNA与ctDNA中的突变类型进行判断。具体的判断条件如下表1：

表1

其中，突变类型中的A代表参考碱基类型，B代表突变碱基类型。待测样本拷贝数缺失的判断标准为：

上式中CN_i是样本第i个候选SNP位点的拷贝数，CCF_i是样本第i个候选SNP位点对应的ctDNA占比。即根据最终预测的占比的变化，拷贝数缺失的标准也随之变化。

针对候选SNP位点突变类型的各种可能性，待测样本在给定CCF_i时的理论突变频率VAF_CCFi可以按照以下表2中的公式进行计算：

表2

考虑到各候选SNP位点满足AD_i～B(DP_i，VAF_i)所代表的分布，对于某一给定CCF_i，可代入上述分布与理论突变频率的公式，计算各位点在该CCF_i下获得相应AD_i与DP_i的似然值。对每一候选SNP位点使用极大似然法，使位点i似然值最大化的CCF_i占比即为该位点的CCF_i预测值。具体的log极大似然公式为：

由于对于基线样本为纯合且待测样本拷贝数无缺失的位点，无法确定待测样本中所含的ctDNA包含杂合突变还是纯合突变，仅根据该位点信息，无法准确计算得该位点的CCF_i。针对该问题，本发明综合该类位点外的其他位点，对这类位点的突变类型和CCF_i进行预测。

首先对其他位点的CCF进行原假设为单峰分布的Hartigans单峰检验：

若检验接受原假设，认为不存在多峰分布。则假定这部分位点均来自于同一肿瘤组织来源，对所有N个SNP的CCF进行极大似然估计，作为CCF_temp：

若检验拒绝原假设，认为存在多峰分布，则对这些位点的CCF_i进行聚类。假定聚到同一簇的候选SNP位点都来自同一肿瘤组织来源，且具有同一CCF_j。聚类的簇数J由gap法确认。完成聚类后，对每一簇的N_j个候选SNP位点的CCF_j分别进行极大似然估计，作为CCF_{j_temp}：

完成其他位点的肿瘤组织来源聚类与CCF_{j_temp}的计算后，对未能确定突变类型的候选SNP位点，分别计算其作为纯合突变、杂合突变且同源于已知聚类簇的似然值，并以最大化似然值的标准对SNP进行分类：

完成上述过程后，所有SNP的突变类型与来源均得到确定。

加入新确定突变类型与来源的候选SNP位点位点，对每一来源簇的候选SNP位点P重复上述的极大似然估计过程，获得ctDNA来源数量J、分来源的占比预测结果CCF_cluster_j。

对占比预测结果进行检验，判断估计结果是否可靠、样本是否存在低ctDNA占比的可能性。具体的判断方法与解决方案为：

(1)如样本包含多个肿瘤组织来源，占比最高的肿瘤组织来源，包含的SNP位点数量低于2个或低于总SNP数量的10％，则认为该组织来源的估计结果不可靠。剔除该估计结果；

(2)如样本最终的最高占比估计结果低于15％，或SNP数量过少(其中，候选SNP位点“数量过少”的判断标准可根据测序数据的不同进行调整，例如，在本发明一实施例中，候选SNP位点数量过少是指当主要肿瘤组织来源包含的候选SNP位点，在基线样本DNA中情况为杂合的数量低于2个，或在基线样本DNA中为纯合的位点数量有任意一个低于2个)，或基线样本纯合位点与杂合位点最终的估计结果存在明显差异，则认为待测样本存在较高的低ctDNA占比可能。剔除所有基线样本的杂合位点，并适当放宽杂合位点的过滤条件，重新提取SNP列表并进行估计。

(3)对于部分样本，重新提取SNP列表后仍未获得足够可用于预测的SNP时，认为该样本的ctDNA占比极低(低于2％)。

完成上述全部分析步骤后，输出最终的ctDNA来源数量J、分来源的占比预测结果CCF_cluster_j为预测结果文件一，输出按照染色体与位置排序的SNP列表为预测结果文件二，其中包含每个SNP对应所属的ctDNA来源与CCF_i。输出路径为给定的输出路径，输出文件前缀由输入参数指定；

待测样本检测结果可视化模块(plot模块)：该模块要求输入输出文件前缀、候选突变位点对应的ctDNA占比估计结果列表的路径与包含待测样本ctDNA来源数量、ctDNA占比的结果文件路径。

该模块首先根据ctDNA占比估计模块输出文件二中提供的各SNP的所属来源与其CCF_i分来源绘制CCF_i概率密度分布图，其中，不同的来源使用不同的颜色进行展示。

再根据该结果，以SNP突变位置为横轴(按染色体与位置排序)、位点对应的ctDNA占比为纵轴，绘制散点图展示各SNP的CCF_i，并根据ctDNA占比估计模块输出文件一提供的每个ctDNA来源的CCF_cluster绘制出水平直线，反应各来源占比高低。同样，该图中使用不同性状或色彩展示不同的来源。

最终，两张图将以左右并列展示的方式输出为pdf格式的图片，对检测结果进行直观展示，方便对检测结果进行解释。输出路径与输入文件的路径一致，输出文件前缀由输入参数指定；自动化ctDNA占比分析模块(autocall)：其整合了上述的待测样本SNP过滤模块(filter模块)、待测样本ctDNA占比估计模块(estimate模块)与待测样本检测结果可视化模块(plot模块)，通过一次性输入所有上述三个模块所要求的信息，可以一次性完成所有的检测步骤，输出全部输出文件。

在本发明一实施方式中，提供两种检测运行模式：第一种即分别运行待测样本SNP过滤模块(filter模块)、待测样本ctDNA占比估计模块(estimate模块)与ctDNA占比检测结果可视化模块(plot模块)；第二种为运行整合了上述三个模块的自动化ctDNA占比检测模块(autocall模块)。

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

实施例1

1.样本的准备

选取6名肿瘤患者(分别记作样本A、样本B、样本C、样本D、样本E和样本F)的血细胞DNA与cfDNA样本，分别进行如下操作：

1)cfDNA的纯化与定量

将6对待测样本中的cfDNA进行2100质检，看是否有大片段的存在。对含有大片段的cfDNA样本进行纯化，由于磁珠具有先吸附大片段的特性，先用0.5倍磁珠纯化，此时大片段吸附在磁珠上，吸取上清液到2.5倍磁珠中，回收所有除大片段后的产物。对所有样本进行qubit定量。

2)血细胞DNA的纯化与定量

取足够量的6对待测样本中的血细胞DNA，用covaris超声打断成200bp左右的DNA片段，用2100质检片段大小是否合适，对片段化的血细胞DNA进行qubit定量。

2.cfDNA梯度样本的制备

使用上述的6对样本，按照下表3中所述的梯度进行稀释，共计32个样本，其中除作为参考样本的血细胞DNA样本外，其余原始cfDNA样本与梯度稀释样本均进行2次重复：

表3

由于建库所需起始量为20ng，为了保证所配样本能够进行重复实验，则每个浓度梯度至少配制了60ng，并平行稀释两份。每次稀释的样本均需充分混匀后，再进行后续样本的稀释。

稀释过程如下所述：

1)1/3原始cfDNA突变频率样本的配制

在30ng cfDNA中加入60ng打断后的血细胞DNA，即将原液稀释了3倍，得到90ng 1/3原始突变频率样本，由于后续稀释需要用到30ng此频率样本，故最后剩余60ng 1/3原始突变频率样本。

2)1/9原始cfDNA突变频率样本的配制

在30ng 1/3原始突变频率样本中加入60ng打断后的血细胞DNA，即稀释了3倍，得到90ng 1/9原始突变频率样本，由于后续稀释需要用到30ng此频率样本，故最后剩余60ng1/9原始突变频率样本。

3)1/27原始cfDNA突变频率样本的配制

在30ng 1/9原始突变频率样本中加入60ng打断后的血细胞DNA，即稀释了3倍，得到90ng 1/27原始突变频率样本，由于后续稀释需要用到20ng此频率样本，故最后剩余70ng1/27原始突变频率样本。

4)1/81原始cfDNA突变频率样本的配制

在20ng 1/27原始突变频率样本中加入40ng打断后的血细胞DNA，即稀释了3倍，得到60ng 1/81原始突变频率样本。

3.cfDNA梯度样本的文库构建、捕获和测序

取用20ng片段化的DNA作为起始量采用KAPA hyper preparation kit(罗氏公司)进行文库构建，经过末端修复、3’端加polyA、连接测序接头、进行无偏向扩增，之后进行纯化获得文库。详述如下：

1)末端平齐并在3’末端加A

反应体系如表4所示：

表4

缓冲液和酶应预先在EP管中混匀，与DNA涡旋混匀后按表5所示反应进行：

表5

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

2)连接接头

对20ng DNA采用7.5μM接头。按照下表6配制反应体系：

表6

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

3)连接后纯化

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

4)PCR扩增

按照下表7配制PCR体系：

表7

充分震荡后快速离心，按照下表8条件进行PCR反应：

表8

5)扩增后纯化

在上一步反应体系中加入50μL Agencourt AMPure XP纯化磁珠，充分涡旋振荡，轻微离心。室温吸附5～15分钟，使DNA与磁珠充分结合。EP管放到磁力架上至液体澄清，缓慢弃上清。加入200μL 80％乙醇，孵育30秒，缓慢丢弃EP管中乙醇。重复用乙醇清洗一次。EP管室温干燥3～5分钟至乙醇完全挥发。从磁力架取下EP管，加入52μL超纯水，涡旋振荡，室温孵育2分钟洗脱DNA。将EP管放至磁力架上吸附至液体澄清，上清转移至新的EP管，取1μL测DNA浓度，剩余的即为所得文库。

4.cfDNA梯度样本的捕获和测序

1)文库捕获

按下表9要求依次加入试剂于新的1.5ml离心管中：

表9

根据文库个数计算样本量，若1个捕获样本加入6个文库，则每个文库需加入167ng。

用移液器吹打混匀，封口膜封住EP管，在膜上插若干小孔，用真空离心浓缩仪在60℃、1350r/min下进行干燥，直至液体完全蒸干。待液体蒸干后，加入如下表10所示的组分：

表10

涡旋震荡混匀，短暂离心以去除管壁残留。于恒温金属浴仪95℃孵育10分钟使DNA变性，短暂离心以去除管壁残留。用移液器将杂交混合液转移至新的PCR管中，加入4.5μl探针，涡旋震荡混匀，短暂离心以去除管壁残留。于PCR仪47℃孵育16～20小时，同时PCR仪加热盖温度设置为57℃以上。

2)捕获产物漂洗

按下表11稀释洗脱缓冲液：

表11

吸取100μl 1×洗脱缓冲液I和400μl 1×洗脱缓冲液IV在47℃预热至少2小时。捕获磁珠室温放置30分钟后使用。取100μl捕获磁珠于新的1.5ml离心管中，将EP管放至磁力架上吸附至液体澄清，用移液器吸去上清。从磁力架上取下离心管，加入200μl 1×磁珠洗脱缓冲液，涡旋震荡混匀。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。

重复一次上述步骤。

向离心管加入100μl 1×磁珠洗脱缓冲液，涡旋震荡混匀。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。

重复一次上述步骤。

向离心管加入100μl 1×磁珠洗脱缓冲液，涡旋震荡混匀。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。将15μl捕获产物加入到磁珠离心管中，用移液器吹打混匀，于47℃孵育45分钟。每间隔15分钟涡旋震荡3秒，使磁珠处于悬浮状态。离心管中加入100μl 47℃预热的1×洗脱缓冲液I，涡旋震荡混匀。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。从磁力架上取下离心管，加入200μl 47℃预热的1×洗脱缓冲液IV，用移液器吹打混匀。于恒温金属浴仪47℃孵育5分钟。

重复一次上述步骤。

将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。取下离心管，每个离心管中分别依次加入200μl未加热的1×洗脱缓冲液I，涡旋震荡2分钟。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。上取下离心管，每个离心管中分别依次加入200μl 1×洗脱缓冲液II，涡旋震荡1分钟。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。取下离心管，每个离心管中分别依次加入200μl 1×洗脱缓冲液III，涡旋震荡30秒。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。取下离心管，加入40μl水，用移液器吹打混匀。

3)捕获产物扩增

在上述40μL混合液中加入如下表12中所示组分：

表12

涡旋震荡混匀，按50μl/管分装量分装到两个新的PCR管中，按以下表13中所示的反应程序扩增：

表13

4)捕获产物纯化

将100μl扩增产物转移至新的1.5ml离心管中，加入180μl纯化磁珠，涡旋震荡混匀。室温静置15分钟。将EP管放至磁力架吸附至液体澄清，用移液器吸去上清。向离心管中加入200μl 80％乙醇，室温静置30秒，用移液器吸去上清。

重复一次上述步骤，室温静置3～5分钟至乙醇完全挥发。

从磁力架取下EP管，加入52μL超纯水，涡旋振荡，室温孵育2分钟洗脱DNA。将EP管放至磁力架吸附至液体澄清，上清转移至新的EP管，即为所得捕获后产物。

取1μL测DNA浓度。

5)测序

1.处理下机fastq数据为各软件可使用的输入文件

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

1)比对

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

2)排序

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

3)筛选

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

4)标记重复

调用Picard-2.1.0的MarkDuplicates模块，对筛选后bam文件中的重复序列进行标记，后续的分析时，会过滤这部分重复序列，采用去重后的数据进行分析；

5)建立索引

调用SAMtools-1.3的index模块对最终生成的bam文件建立索引，生成与标记重复后的bam文件配对的bai文件。

6.用其他软件对cfDNA梯度样本的ctDNA占比进行检测

1)FACETS

在检测时，使用样本A-0、样本B-0、样本C-0、样本D-0、样本E-0、样本F-0(即6名患者的血细胞DNA样本)分别作为参考样本，采用待测样本与该样本对应参考样本的bam文件、预先根据记录捕获区间的bed文件、人类参考基因组序列的fasta文件作为输入文件。首先使用SAMtools的mpileup模块，根据各样本的bam文件、bed文件、人类参考基因组序列的fasta文件生成mpileup文件；再利用VarScan的mpileup2cns模块，根据mpileup文件，生成各样本的突变列表vcf文件；使用待测样本与对应参考样本的vcf文件，生成可供FACETS使用的snpmat文件；最终利用FACETS对样本的肿瘤纯度(在本实例中代表ctDNA占待测样本的比例)与染色体倍性进行预测。设置参数时均使用该软件的默认参数。各样本具体的肿瘤纯度与染色体倍性结果如下表14：

表14

其中理论ctDNA占比由cfDNA占比为1的两个样本肿瘤纯度的均值乘cfDNA占比计算获得；未能检测出肿瘤纯度的样本，用“-”进行标记，并在计算R-square时，将其按照0％进行处理，与cfDNA占比进行相关性计算。

此外，同一样本的两次技术重复，经FACETS估计所得的ctDNA占比结果，相关系数为0.8951。

FACETS对染色体倍性的检测结果相对较为准确，但对肿瘤纯度的检测准确度极低。在检测高ctDNA占比的样本时，FACETS的检测结果偏低，同时，其在检测低ctDNA占比的样本时，多数样本无法给出检测结果，少数样本的检测结果偏高。

2)PureCN

在检测时，使用样本A-0、样本B-0、样本C-0、样本D-0、样本E-0、样本F-0(即6名患者的血细胞样本)分别作为参考样本，采用待测样本与该样本对应参考样本的bam文件、预先根据记录捕获区间的bed文件、人类参考基因组序列的fasta文件与bigWig文件生成的捕获信息interval文件作为输入文件。首先利用各样本的bam文件与捕获信息interval文件，使用PureCN的coverage模块生成coverage文件；再使用mutect软件，对各待测样本与参考样本的bam生成的SNP列表vcf文件；最终使用待测样本与对应参考样本的coverage文件、待测样本的vcf文件，通过PureCN对各待测样本的肿瘤纯度与染色体倍性进行检测，设置参数时均使用该软件的默认参数。各样本具体的肿瘤纯度与染色体倍性结果如下表15：

表15

其中理论ctDNA占比由cfDNA占比为1的两个样本肿瘤纯度的均值乘cfDNA占比计算获得。此外，同一样本的两次技术重复，经PureCN估计所得的ctDNA占比结果，相关系数为0.5850。

PureCN对所有样本的染色体倍性预测明显偏高，对同一样本检测结果的一致性较差，对低浓度血浆稀释样本的肿瘤纯度预测偏高，不适合对低ctDNA占比的样本进行检测。

7.用本发明对cfDNA梯度样本的ctDNA占比进行检测

在检测时，使用样本A-0、样本B-0、样本C-0、样本D-0、样本E-0、样本F-0(即6名患者的血细胞样本)分别作为参考样本，采用待测样本与该样本对应参考样本的bam文件、人类参考基因组序列的fasta文件作为输入文件。首先使用SAMtools的mpileup模块，根据各样本的bam文件、bed文件、人类参考基因组序列的fasta文件生成mpileup文件；再利用VarScan的mpileup2cns模块，根据mpileup文件，生成各样本的突变列表vcf文件，并使用Annovar软件对vcf文件进行注释；另外使用cnvkit，利用输入文件获得待测样本的拷贝数变异结果cns文件；整理注释后突变列表文件与cns文件为本发明所需格式，分别作为本发明的输入文件，使用autocall模块检测，各浓度梯度样本检测结果如下表16所示：

表16

其中理论ctDNA占比由cfDNA占比为1的两个样本主要肿瘤来源占比的均值乘cfDNA占比计算获得。

使用本发明的技术方案(具体可参见具体实施方式部分)，无论对高ctDNA占比或低占比的样本，结果都很准确，且与稀释梯度的线性关系非常明显，全部6组样本R-square均超过0.9。即使1/81稀释梯度样本的ctDNA占比也能稳定检出。此外，同一样本的两次技术重复，经本发明估计所得的ctDNA占比结果，相关系数为0.9985，检测结果一致性高。

与其余肿瘤细胞占比检测软件相比，本发明的技术方案更加稳定，无论是对高、低占比样本的检测，都具有更高的灵敏度，与理论ctDNA占比的一致性更高，检测下限也更低，且能够输出不同ctDNA来源的占比，更适合针对ctDNA的检测。

使用plot模块，对样本F-1与样本F-2的检测结果分别进行图片展示，如图2所示。

其中右侧散点图的横坐标代表了对应SNP的位置(按染色体与位置排序)，纵坐标代表SNP对应的占比检测结果，水平直线为该ctDNA来源(图中显示为cluster)的占比，颜色或形状代表不同的ctDNA来源。左图为根据右图中SNP与对应占比检测结果的概率密度分布图。

由图2可以看出1/3稀释梯度样本各ctDNA来源(图中显示为cluster)的占比检测结果与未稀释样本线性相关，且各来源包含的SNP分布也较为一致。

从以上的描述中，可以看出，本发明上述的实施例实现了如下技术效果：

1)同时使用基线样本DNA中的纯合与杂合位点，对每个突变进行统计检验，使用严格的筛选标准，确保挑选出用于计算ctDNA占比的位点均为可靠位点。目前针对肿瘤占比的检测主要通过使用基线样本与待测样本的突变信息进行。为避免低频体细胞突变与背景噪音混杂，基线样本中的纯合位点通常会被检测软件剔除。该做法虽然能够排除背景噪音对结果的影响，但也会浪费大量的可用信息，并且在检测低浓度样本时，会严重影响准确性。本发明通过使用统计检验方法，对低频突变与背景噪音进行了有效的区分，在避免背景噪音对检测结果影响的同时，保留了更多的有效信息，提高了检测方法的灵敏度。此外，本发明还通过使用内置的多态性基因列表对突变进行了进一步过滤，剔除了由于多态性而造成的突变。

2)仅使用拷贝数稳定或拷贝数缺失区域进行估计。现有检测软件，通常在检测肿瘤占比的同时，对突变区域的拷贝数不进行限定。但拷贝数扩增区域的情况复杂，会引入较大的不确定性，继而影响后续的概率模型的估计效果。本发明剔除了拷贝数扩增区域的突变，降低了不确定性，在不降低模型准确度的同时，使后续的检测步骤得以简化。

3)使用极大似然法，既分别对每一个有效突变进行ctDNA的占比估计与ctDNA来源的异同进行判断，也对待测样本整体的ctDNA占比情况进行检测。由于ctDNA来源于多个肿瘤组织，每个突变可能来自于不同的肿瘤组织，对应的ctDNA占比也各不相同。

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

Claims

1.一种基于捕获测序的ctDNA占比的检测方法，其特征在于，包括以下步骤：

S1，获取基线样本DNA和待测cfDNA的捕获测序的基因数据，进行所述基因数据的处理，获得基线样本DNA注释后突变列表、待测cfDNA注释后突变列表和待测cfDNA拷贝数变异列表；

S2，同时使用所述基线样本DNA中的纯合与杂合位点，对所述待测cfDNA中的每个突变进行统计检验，挑选出所述基线样本DNA和所述待测cfDNA之间突变频率存在显著不同且满足预设过滤条件的位点作为候选SNP位点；

S3，结合所述候选SNP位点所在区域拷贝数突变的情况判断正常细胞DNA与ctDNA的所述候选SNP位点的突变类型；以及

S4，利用所述候选SNP位点不同碱基的支持reads数量与拷贝数，使用极大似然法建立概率模型，通过所述概率模型计算得到ctDNA占比。

2.根据权利要求1所述的检测方法，其特征在于，所述概率模型的评估指标包括：对每一个所述候选SNP位点进行ctDNA的占比估计、对ctDNA来源的异同进行判断以及对所述待测cfDNA整体的ctDNA占比情况进行检测。

3.根据权利要求1所述的检测方法，其特征在于，在所述概率模型计算得到ctDNA占比后通过可视化模块绘制散点图和概率密度分布图对结果进行展示。

4.根据权利要求1所述的检测方法，其特征在于，所述S3中所用的判断标准如下表1所示：

表1

其中，突变类型中的A代表参考碱基类型，B代表突变碱基类型，所述待测cfDNA拷贝数缺失的判断标准为：

上式中CN_i是样本第i个所述候选SNP位点的拷贝数，CCF_i是样本第i个所述候选SNP位点对应的ctDNA占比，即根据最终预测的ctDNA占比的变化，所述待测cfDNA拷贝数缺失的标准也随之变化。

5.根据权利要求4所述的检测方法，其特征在于，针对所述候选SNP位点的突变类型的各种可能性，所述待测cfDNA在给定CCF_i时的理论突变频率按照以下表2中的公式进行计算：

表2

所述候选SNP位点的总测序深度、突变碱基支持reads数与突变频率满足如下公式所代表的分布：

AD_i～B(DP_i，VAF_i)

其中，AD_i是样本第i个所述候选SNP位点的突变碱基支持reads数；DP_i是样本第i个所述候选SNP位点的总测序深度；VAF_i是样本第i个所述候选SNP位点的突变频率；上式表示AD_i服从n＝DP_i，p＝VAF_i的二项分布；对于某一给定CCF_i，代入所述分布与理论突变频率的公式，计算各位点在该CCF_i下获得相应AD_i与DP_i的似然值；对每一个所述候选SNP位点使用极大似然法，使位点i似然值最大化的CCF_i占比即为该位点的CCF_i预测值，具体的log极大似然公式为：

6.根据权利要求5所述的检测方法，其特征在于，所述检测方法还包括：对所述基线样本DNA的所述候选SNP位点为纯合且所述待测cfDNA拷贝数无缺失的位点的突变类型和CCF_i进行预测，具体包括：

对所述基线样本DNA的所述候选SNP位点为纯合且所述待测cfDNA拷贝数无缺失的位点外的其他位点的CCF进行原假设为单峰分布的Hartigans单峰检验：

若检验接受原假设，认为不存在多峰分布，则假定这部分位点均来自于同一肿瘤组织来源，对所有N个所述候选SNP位点的CCF进行极大似然估计，作为CCF_temp：

若检验拒绝原假设，认为存在多峰分布，则对这些位点的CCF_i进行聚类，假定聚到同一簇的所述候选SNP位点都来自同一肿瘤组织来源，且具有同一CCF_j；聚类的簇数J由gap法确认；完成聚类后，对每一簇的N_j个所述候选SNP位点的CCF_j分别进行极大似然估计，作为CCF_{j_temp}：

完成其他位点的肿瘤组织来源聚类与CCF_{j_temp}的计算后，对未能确定突变类型的所述候选SNP位点，分别计算其作为纯合突变、杂合突变且同源于已知聚类簇的似然值，并以最大化似然值的标准对所述候选SNP位点进行分类：

完成上述过程后，所有所述候选SNP位点的突变类型与来源均得到确定；

加入新确定的突变类型与来源的所述候选SNP位点，对每一来源簇的所述候选SNP位点重复上述的极大似然估计过程，获得ctDNA来源数量J、分来源的占比预测结果CCF_cluster_j。

7.根据权利要求6所述的检测方法，其特征在于，所述检测方法还包括对占比预测结果进行检验的步骤，具体包括：

如所述待测cfDNA包含多个肿瘤组织来源，占比最高的肿瘤组织来源，包含的所述候选SNP位点数量低于2个或低于总所述候选SNP位点数量的10％，则认为该组织来源的估计结果不可靠，剔除该估计结果；

如所述待测cfDNA最终的最高占比估计结果低于15％，或所述候选SNP位点数量过少，或所述基线样本DNA纯合位点与杂合位点最终的估计结果存在明显差异，则认为所述待测cfDNA存在较高的低ctDNA占比可能，剔除所有基线样本DNA的杂合位点，并放宽杂合位点的过滤条件，重新提取所述候选SNP位点的列表并进行估计；

对于部分所述待测cfDNA，重新提取所述候选SNP位点列表后仍未获得足够可用于预测的所述候选SNP位点时，认为所述待测cfDNA的ctDNA占比极低，所述极低为低于2％。

8.根据权利要求1所述的检测方法，其特征在于，所述S2中，对所述待测cfDNA中的每个SNP位点的突变进行统计检验包括：

AD_i～B(DP_i，VAF_i)

1)在所述基线样本DNA中是否为纯合，基线样本突变频率＝0％，

假定背景误差为0.01％，检验所述基线样本DNA；

原假设：VAF_i＜0.0001；

备择假设：VAF_i≥0.0001；

如上述检验p-value大于0.05，接受原假设，所述基线样本DNA在该位点为纯合；

2)在所述基线样本DNA中是否为杂合，基线样本突变频率＝50％；

杂合突变频率为50％；

原假设：VAF_i＝0.5；

备择假设：VAF_i≠0.5；

如上述检验p-value大于0.05，接受原假设，所述基线样本DNA在该位点为杂合；

3)在所述待测cfDNA中是否为突变，频率是否高于背景误差与测序错误，假定背景误差为0.01％，所述待测cfDNA的样本测序错误为该样本所有突变频率低于1％的SNP位点突变频率的中位数；

原假设：VAF_i＜max{0.0001，median(all SNP VAF)}；

备择假设：VAF_i≥max{0.0001，median(all SNP VAF)}；

如上述检验p-value小于0.05，拒绝原假设，所述待测cfDNA在该位点为突变；

4)在所述待测cfDNA中是否为杂合，所述基线样本DNA突变频率＝50％，

杂合突变频率为50％；

原假设：VAF_i＝0.5；

备择假设：VAF_i≠0.5；

如上述检验p-value大于0.05，接受原假设，所述待测cfDNA在在该位点为杂合；

另外，使用Fisher精确检验，检验所述待测cfDNA中的突变频率是否与所述基线样本DNA存在显著差异，所述待测cfDNA突变频率≠所述基线样本DNA突变频率：

原假设：VAF_cfDNA_i＝VAF_baseline_i；

备择假设：VAF_cfDNA_i≠VAF_baseline_i；

如上述检验p-value小于0.05，拒绝原假设，所述待测cfDNA的突变频率与所述基线DNA样本存在显著差异；

优选的，所述S2中预设过滤条件如下：

如位点在所述待测cfDNA与所述基线样本DNA中深度低于50x，剔除该SNP；

如位点所在基因的拷贝数高于2.2或低于1，剔除该SNP；

如位点在所述基线样本DNA中为纯合，要求所述待测cfDNA突变频率高于所述基线样本DNA，并且在检验中与背景噪音或所述基线样本DNA存在明显差异；

如位点在所述基线样本DNA中为杂合，要求所述待测cfDNA非杂合，且突变频率与所述待测cfDNA存在明显差异。

9.一种基于捕获测序的ctDNA占比的检测装置，其特征在于，所述装置用于存储或者运行模块，或者所述模块为所述装置的组成部分；其中，所述模块为软件模块，所述软件模块为一个或多个，所述软件模块用于执行如权利要求1至8中任一项所述的检测方法。

10.根据权利要求9所述的检测装置，其特征在于，所述软件模块包括：待测cfDNA SNP过滤模块、待测cfDNA的ctDNA占比估计模块和ctDNA占比检测结果可视化模块；或

所述软件模块包括整合了所述待测cfDNA SNP过滤模块、待测cfDNA的ctDNA占比估计模块和ctDNA占比检测结果可视化模块的自动化模块。