基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法

摘要

本发明提供了一种基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法，装置中包括：FASTQ文件处理模块，用于获取待测ctDNA样本捕获测序的FASTQ文件，并处理得到过滤后的FASTQ文件；待测样本比对模块，用于将得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，得到对应的Bam文件；reads水平过滤模块，用于根据预先设定的C‑T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；甲基化水平预测模块，用于根据目标区域Bed文件及预先设定的各reads中覆盖CpG位点的数量进一步对Bam文件进行过滤，并根据剩余reads对CpG位点的甲基化水平进行预测。

说明书

技术领域

本发明涉及生物医学技术领域，尤其涉及一种ctDNA甲基化水平预测装置及方法。

背景技术

循环肿瘤DNA(circulating tumor，ctDNA)是一类来源于肿瘤细胞凋亡、坏死的DNA小片段，由肿瘤细胞释放到外周血循环后形成内源性单链或者双链DNA，携带有与原发肿瘤组织相一致的分子突变信息。因此，ctDNA样本检测可作为临床上组织样本基因检测的替代样本。

研究表明，表观遗传学变化是肿瘤形成最常见的分子变化之一。DNA甲基化是一种被广泛研究的表观遗传修饰方式，在调控基因表达等方面发挥了重要作用。通常地，DNA甲基化是指5-甲基胞嘧啶(5mC)在DNA甲基转移酶(DNMT)的作用下将甲基基团添加到胞嘧啶的5’C上形成的结构。研究表明，DNA甲基化参与细胞分化、组织特异性基因表达等细胞活动，异常的DNA甲基化会导致发育异常和肿瘤等疾病的发生。因此，DNA甲基化对个体发育和肿瘤的发生发展机制都具有重要意义。

随着二代测序技术的不断发展，其在人类遗传病和癌症诊断领域的应用越来越普遍，ctDNA的甲基化测序已经成为研究肿瘤发生发展机制的重要手段。然而，人类参考基因组大小为3G，进行全基因组甲基化测序成本过高，数据量较大。因此，目标区域捕获测序已成为科学研究中较为理想的方法。

当前传统的DNA甲基化捕获数据的质量检测过程一般为：将FASTQ格式的数据与人类参考基因组进行比对，保留高质量的唯一比对reads，并去除重复的reads，之后评估保留下的reads的碱基含量比例、捕获效率和测序深度，得到待测ctDNA样本的Bam文件，最后利用第三方软件对Bam文件进行分析得到待测ctDNA样本于CpG位点(胞嘧啶-磷酸-鸟嘌呤位点，即DNA序列中胞嘧啶后紧连鸟嘌呤的位点)的甲基化水平数据，直接用于后续的科学研究分析中。

在上述目标区域DNA甲基化捕获测序过程中需要进行重亚硫酸盐处理，将所有未甲基化的胞嘧啶(C)转化为尿嘧啶(U)及将尿嘧啶经过PCR(聚合酶链式反应，一种用于放大扩增特定的DNA片段的技术)扩增转变成胸腺嘧啶(T)，但是发生甲基化的胞嘧啶在这个过程中不会发生改变。可知，这一过程中很可能出现未甲基化的胞嘧啶转化不完全的现象，进而导致待测ctDNA样本甲基化水平出现预测偏差。且由于ctDNA的含量很低，ctDNA样本的甲基化水平更容易受到C-T转化率的影响，进而影响检测结果的准确性。

发明内容

针对上述问题，本发明提供了一种基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法，有效解决现有ctDNA甲基化水平预测中存在的准确性低、数据质量偏差大等缺陷。

本发明提供的技术方案如下：

一方面，本发明提供了一种基于目标区域捕获测序的ctDNA甲基化水平预测装置，包括：

FASTQ文件处理模块，用于获取待测ctDNA样本捕获测序的FASTQ文件，并对其进行预处理操作得到过滤后的FASTQ文件；

待测样本比对模块，用于将所述FASTQ文件处理模块得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，得到对应的Bam文件；

reads水平过滤模块，用于根据预先设定的C-T转化率对所述待测样本比对模块生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；

甲基化水平预测模块，用于根据目标区域Bed文件及预先设定的各reads中覆盖CpG位点的数量进一步对所述reads水平过滤模块输出的Bam文件进行过滤，并根据剩余reads对CpG位点的甲基化水平进行预测。

在本技术方案中，FASTQ为一种常见的高通量测序文件类型。reads为测序读长，测序仪测到的基因组或转录组序列片段。根据甲基化的C碱基所处的上下文环境，分为三种类型CpG、CHG和CHH，H代表除了G碱基之外的其他碱基，即A、C、T中的任意一种；CpG为甲基化的C的下游是1个G碱基，CHG代表甲基化的C下游的2个碱基是H和G，CHH表示甲基化的C下游的两个碱基都是H，CHG和CHH可以合称为non CpG context。Bam文件用来存储测序序列回贴到参考基因组的结果。C-T转化率为原始序列非CpG位点的C碱基转化为T碱基的比例。

进一步优选地，在所述FASTQ文件处理模块中，对获取的FASTQ文件进行的预处理操作包括：去除接头和低质量reads；和/或，

在所述待测样本比对模块中，将所述FASTQ文件处理模块得到的FASTQ文件中的基因序列分别与人类参考基因组和内参lambda DNA参考基因组进行比对并去重，生成人类参考基因组的Bam文件、去重前的比对报告和去重后的比对报告，及内参lambda DNA参考基因组Bam文件、去重前的比对报告和去重后的比对报告。

进一步优选地，在所述reads水平过滤模块中，包括：

甲基化数量统计单元，用于逐行读取所述待测样本比对模块生成的Bam文件中的reads，并对其non-CpG context模式下甲基化和非甲基化的碱基数量进行统计；

C-T转化率计算单元，用于根据发生甲基化的non-CpG context的碱基数及non-CpGcontext碱基数总和对每条reads的C-T转化率进行计算；

第一过滤单元，用于将Bam文件中C-T转化率小于预先设定的C-T转化率的reads滤除，得到滤后的Bam文件。

进一步优选地，在所述甲基化水平预测模块中，包括：

第二过滤单元，用于根据目标区域Bed文件对dbSNP数据库中已知的SNP位点及因特定变异原因产生的SNP位点进行过滤得到待测ctDNA样本的CpG位点；及用于根据过滤得到的CpG位点及预先设定的各reads中覆盖CpG位点的数量进一步对所述reads水平过滤模块输出的Bam文件进行过滤；

甲基化水平计算单元，用于根据第二过滤单元过滤后Bam文件剩余的reads计算CpG位点的甲基化水平。

另一方面，本发明提供了一种基于目标区域捕获测序的ctDNA甲基化水平预测方法，包括：

获取待测ctDNA样本捕获测序的FASTQ文件，并对其进行预处理操作得到过滤后的FASTQ文件；

将得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，得到对应的Bam文件；

根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；

根据目标区域Bed文件及预先设定的各reads中覆盖CpG位点的数量对过滤后的Bam文件进一步过滤，并根据剩余reads对CpG位点的甲基化水平进行预测。

进一步优选地，在所述获取待测ctDNA样本捕获测序的FASTQ文件，并对其进行预处理操作得到过滤后的FASTQ文件中，包括：对获取的FASTQ文件进行去除接头和低质量reads操作；和/或，

在将得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，得到对应的Bam文件中，包括：将所述FASTQ文件处理模块得到的FASTQ文件中的基因序列分别与人类参考基因组和内参lambda DNA参考基因组进行比对并去重，生成人类参考基因组的Bam文件、去重前的比对报告和去重后的比对报告，及内参lambda DNA参考基因组Bam文件、去重前的比对报告和去重后的比对报告。

进一步优选地，在所述根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件中，包括：

逐行读取Bam文件中的reads，并对其non-CpG context模式下甲基化和非甲基化的碱基数量进行统计；

根据发生甲基化的non-CpG context的碱基数及non-CpG context碱基数总和对每条reads的C-T转化率进行计算；

将Bam文件中C-T转化率小于预先设定的C-T转化率的reads滤除，得到滤后的Bam文件。

进一步优选地，所述根据目标区域Bed文件及预先设定的各reads中覆盖CpG位点的数量对过滤后的Bam文件进一步过滤，并根据剩余reads对CpG位点的甲基化水平进行预测中，包括：

根据目标区域Bed文件对dbSNP数据库中已知的SNP位点及因特定变异原因产生的SNP位点进行过滤得到待测ctDNA样本的CpG位点；

根据过滤得到的CpG位点及预先设定的各reads中覆盖CpG位点的数量进一步对Bam文件进行过滤；

根据滤后Bam文件剩余的reads计算CpG位点的甲基化水平。

另一方面，本发明提供了一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器运行所述计算机程序时实现上述基于目标区域捕获测序的ctDNA甲基化水平预测方法的步骤。

另一方面，本发明提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现上述任一项所述基于目标区域捕获测序的ctDNA甲基化水平预测方法的步骤。

本发明提供的基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法中，至少能够带来以下有益效果：

1.在传统甲基化数据质量检测流程的基础上，考虑了C-T转化率对后续预测甲基化水平产生的影响，使用严格的筛选标准，确保过滤后的甲基化数据具有较高的可靠性。具体，考虑到待测ctDNA样本-的特殊性，对每条reads的C-T转化率进行统计，并过滤因C-T转化率低产生的reads噪音，大大提高了甲基化数据的可靠性，为后续甲基化水平预测奠定了基础。

2.在常见的CpG位点甲基化水平预测方法的基础上，采用更为严格的甲基化水平预测标准，使甲基化水平预测更为准确。具体，考虑覆盖CpG位点的reads的甲基化状态，过滤掉可信度不高的reads，使甲基化水平预测更为准确，为科学研究提供可靠的数据依据。

附图说明

下面将以明确易懂的方式，结合附图说明优选实施方式，对上述特性、技术特征、优点及其实现方式予以进一步说明。

图1为本发明中基于目标区域捕获测序的ctDNA甲基化水平预测装置结构示意图；

图2为本发明中基于目标区域捕获测序的ctDNA甲基化水平预测方法流程示意图；

图3为本发明中终端设备结构示意图。

附图标记：

100-ctDNA甲基化水平预测装置，110-FASTQ文件处理模块，120-待测样本比对模块，130-reads水平过滤模块，140-甲基化水平预测模块。

具体实施方式

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。

在现有技术中，针对目标区域甲基化捕获数据的质量检测，主要集中在与参考基因组的比对率、碱基分布比例、碱基含量比例、捕获效率和测序深度，对待测ctDNA样本的C-T转化率未做考虑，由于目标区域甲基化测序的建库方式和ctDNA样本的特殊性，C-T转化率可能会导致未甲基化的胞嘧啶转化不完全，从而对ctDNA样本甲基化水平产生较大的预测偏差。另外，当前各种预测CpG位点甲基化水平的软件参差不齐，针对CpG位点甲基化水平预测的算法大多数集中于发生甲基化的reads数除以未发生甲基化和发生甲基化的reads数之和，而未考虑每条reads上含有的CpG位点的数量及状态，同样会导致ctDNA样本甲基化水平的预测偏差，不能保证数据结果的准确性和可靠性，从而对数据解读产生影响。基于此，本发明提供了一种全新的基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法，提高甲基化水平预测的准确性，为科学研究提供可靠的数据依据。

本发明的第一实施例，如图1所示，一种基于目标区域捕获测序的ctDNA甲基化水平预测装置100，包括：FASTQ文件处理模块110，用于获取待测ctDNA样本捕获测序的FASTQ文件，并对其进行预处理操作得到过滤后的FASTQ文件；待测样本比对模块120，用于将FASTQ文件处理模块110得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，得到对应的Bam文件；reads水平过滤模块130，用于根据预先设定的C-T转化率对待测样本比对模块120生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；甲基化水平预测模块140，用于根据目标区域Bed文件及预先设定的各reads中覆盖CpG位点的数量进一步对reads水平过滤模块130输出的Bam文件进行过滤，并根据剩余reads对CpG位点的甲基化水平进行预测。

在该ctDNA甲基化水平预测装置100中，首先，FASTQ文件处理模块110对获取的FASTQ文件进行去除接头和低质量reads操作，得到不含有接头和低质量碱基的FASTQ格式数据。之后，待测样本比对模块120将FASTQ文件处理模块110得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，保留高质量没有重复的唯一比对的reads，得到待测ctDNA样本的Bam文件。接着，reads水平过滤模块130对得到的Bam文件中reads的non-CpGcontext进行C-T转化率的评估，过滤掉不合格的reads，得到可用于后续分析的Bam文件；最后通过甲基化水平预测模块140对Bam文件的过滤和分析，得到待测ctDNA样本的CpG位点较为准确的甲基化水平数据。

FASTQ文件处理模块110获取待测ctDNA样本捕获测序的FASTQ文件之后，利用去接头软件Trimmomatic去除接头和低质量reads，得到过滤后的FASTQ文件，并利用FASTQC（一种高通量测序数据的质控软件，用于数据质量的评估）软件对待测ctDNA样本的数据量、碱基质量分布、碱基含量比例进行统计分析。具体，在切除接头序列后，切去剩余部分开头和结尾处碱基质量低于20的碱基，从reads的5’端开始，以大小为5的窗口进行划窗计算平均质量，如果窗口内平均碱基质量低于20，则切除该窗口，并要求切除后剩余碱基数量超过75。

待测样本比对模块120利用基因组比对工具Bismark（一种比对方法软件，用于查找测序序列在基因参考序列中的位置，可输出Bam格式结果文件）将FASTQ文件处理模块110得到的FASTQ文件中的基因序列分别与人类参考基因组和内参lambda DNA参考基因组进行比对并去重，生成人类参考基因组的Bam文件、去重前的比对报告和去重后的比对报告，及内参lambda DNA参考基因组Bam文件、去重前的比对报告和去重后的比对报告；并利用SAMtools和Picard工具对比对后的Bam文件进行排序和标记重复处理。这一过程中，输入为待测ctDNA样本原始数据路径及待测ctDNA样本名称。

reads水平过滤模块130的输入为待测ctDNA样本比对到参考基因组且去除重复的Bam文件路径及对non-CpG context C-T转化率的最低要求。在过滤过程中，首先甲基化数量统计单元逐行读取待测样本比对模块120生成的Bam文件中的reads，根据Bam文件中每条read中原始序列为C碱基的位点的实际碱基情况，对其non-CpG context模式下甲基化和非甲基化的碱基数量进行统计；之后，C-T转化率计算单元根据发生甲基化的non-CpGcontext的碱基数及non-CpG context碱基数总和（甲基化和非甲基化的碱基数量总和）对每条reads的C-T转化率进行计算；最后，第一过滤单元将Bam文件中C-T转化率小于预先设定的C-T转化率的reads滤除，以此过滤掉不符合non-CpG context C-T转化率最低要求的reads，输出过滤后的Bam格式文件、过滤后待测ctDNA样本C-T转化率及过滤后待测ctDNA样本的reads数据量。

甲基化水平预测模块140的输入为reads水平过滤模块130过滤后得到的Bam文件路径、目标区域Bed文件及每条reads上覆盖CpG位点的最低要求。在预测过程中，首先，第二过滤单元利用BisSNP软件（一种分析甲基化数据的软件，可用于鉴定甲基化位点和预测甲基化水平）根据目标区域Bed文件对dbSNP数据库中已知的SNP位点及因特定变异原因（如，结构变异、染色体拷贝数变异等）产生的SNP位点进行过滤得到待测ctDNA样本的CpG位点；之后，根据过滤得到的CpG位点及预先设定的各reads中覆盖CpG位点的数量（即上述每条reads上覆盖CpG位点的最低要求）进一步对reads水平过滤模块130输出的Bam文件进行过滤，将不满足覆盖CpG位点数量最低要求的reads滤除；最后，甲基化水平计算单元根据第二过滤单元过滤后Bam文件剩余的reads计算CpG位点的甲基化水平，每个CpG位点甲基化水平计算公式为覆盖到该CpG位点上发生甲基化的CpG数量符合最低要求的reads数除以覆盖到该位点上的CpG数量符合最低要求的所有reads数。于此同时，利用Bedtools软件（一种处理基因组算法的工具），结合Bed文件，对reads水平过滤模块130过滤后的Bam文件进行处理，得到待测ctDNA样本的捕获效率；利用SAMtools软件（一种处理bam/sam文件的工具）对过滤后的Bam文件进行处理，得到待测ctDNA样本在目标区域每个位点的测序深度，并统计待测ctDNA样本的平均测序深度等数据。

在实际应用中，上述FASTQ文件处理模块110、待测样本比对模块120、reads水平过滤模块130及甲基化水平预测模块140的处理过程可以分别进行，即独立模块化进行，也可以整合在一起自动化完成所有过程。在自动化甲基化数据质量检测和甲基化水平预测过程中：一次性输入：待测ctDNA样本进行甲基化目标捕获测序的FASTQ文件和目标区域Bed文件（包含染色体、起始点、终止点三列信息）。输出文件包括：待测ctDNA样本数据统计表（包含原始碱基数据量、原始reads数据量、过滤后碱基数据量、过滤后reads数据量、比对到参考基因组reads数据量和比例、去重比例、去重后reads数据量、Total C碱基含量、甲基化C碱基含量，非甲基化C碱基含量、CpG context和non-CpG context中甲基化的C碱基含量、CpGcontext和non-CpG context中非甲基化的C碱基含量、reads水平过滤前样本C-T转化率、内参样本lambda C-T转化率、reads水平过滤后样本C-T转化率、reads水平过滤后数据量、目标区域碱基数量、目标区域捕获的数据量和比例、不同测序深度下目标区域捕获的碱基数量及比例及平均测序深度）、及待测ctDNA样本目标区域CpG位点的甲基化水平（包含染色体、起始点、终止点、甲基化水平及测序深度五列信息）。

相对应地，本发明还提供了一种基于目标区域捕获测序的ctDNA甲基化水平预测方法，应用于上述ctDNA甲基化水平预测装置，如图2所示，该ctDNA甲基化水平预测方法包括：S10 获取待测ctDNA样本捕获测序的FASTQ文件，并对其进行预处理操作得到过滤后的FASTQ文件；S20 将得到的FASTQ文件中的基因序列与参考基因组进行比对并去重，得到对应的Bam文件；S30 根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；S40 根据目标区域Bed文件及预先设定的各reads中覆盖CpG位点的数量对过滤后的Bam文件进一步过滤，并根据剩余reads对CpG位点的甲基化水平进行预测。

具体，在步骤S20中，包括：将FASTQ文件处理模块得到的FASTQ文件中的基因序列分别与人类参考基因组和内参lambda DNA参考基因组进行比对并去重，生成人类参考基因组的Bam文件、去重前的比对报告和去重后的比对报告，及内参lambda DNA参考基因组Bam文件、去重前的比对报告和去重后的比对报告。在步骤S30中包括：逐行读取Bam文件中的reads，并对其non-CpG context模式下甲基化和非甲基化的碱基数量进行统计；根据发生甲基化的non-CpG context的碱基数及non-CpG context碱基数总和对每条reads的C-T转化率进行计算；将Bam文件中C-T转化率小于预先设定的C-T转化率的reads滤除，得到滤后的Bam文件。在步骤S40中包括：根据目标区域Bed文件对dbSNP数据库中已知的SNP位点及因特定变异原因产生的SNP位点进行过滤得到待测ctDNA样本的CpG位点；根据过滤得到的CpG位点及预先设定的各reads中覆盖CpG位点的数量进一步对Bam文件进行过滤；根据滤后Bam文件剩余的reads计算CpG位点的甲基化水平。

以下通过一实例对基于目标区域捕获测序的ctDNA甲基化水平预测方法及其有益效果进行说明：

1.样本准备

选取6名肿瘤患者的ctDNA样本进行文库构建、目标区域捕获和测序，每名患者重复2次，分别进行如下操作：

1.1 处理血浆

1.1.1样本融化后，每1mL样本中加入15μL蛋白酶K（Proteinase K）(20mg/mL)和50 μL十二烷基硫酸钠（SDS）溶液(20%)。若血浆量不足4mL，用磷酸缓冲盐溶液（PBS）补足。翻转混匀，60℃孵育20min，然后冰浴5min。

1.1.2 向深孔板中加入试剂，各深孔板中添加的试剂及对应的量如表1所示：

表1：深孔板中加入的试剂列表

1.1.3 运行KingFisher FLEX磁珠提取仪

程序运行前将干净磁头套放入检测程序指定位置，运行程序检测磁头套是否会掉落。深孔板加好后，点击自动提取仪上的SATRT键，按照显示屏要求依次放入磁头套和对应的深孔板。再次点击SATRT键，自动提取仪开始运行。

1.1.4 吸出DNA样品：

自动提取仪运行结束后，先取出7号深孔板，然后点击STOP键。用移液器将DNA样本吸出至对应的贴标标签的离心管中。

1.2 cfDNA文库构建

1.2.1 内参准备

取Lamdba DNA加入50uL打断管中，使用M220打断仪打断，将打断的内参DNA稀释，建库时加入样本中。Lamdba为参考品，用于确定样本的转化情况。

1.2.2 DNA样本的准备

将提取的6例肿瘤患者血浆按照10ng总量分成2份，加入打断的参考品，准备建库，这里cfDNA样本不需要打断。样本操作信息如表2所示。

表2：样本操作信息列表

1.3文库制备步骤：

1.3.1使用EZ DNA Methylation-Lightning

1.3.1.1 样品起始体积为20μL。当不足20μL时，用水补足。

1.3.1.2 取130μL试剂盒中的Lightning Conversion Reagent加入DNA样本中，震荡混匀，短暂离心，置于PCR仪上，按表3的条件进行PCR反应。

表3：PCR反应的条件

1.3.1.3 向试剂盒中的Zymo-Spin

1.3.1.4 将收集管中的液体重新加回吸附柱中，静置2min，12000rpm离心1min，弃废液。

1.3.1.5 加入100μL试剂盒中的M-Wash Buffer，12000rpm离心1min，弃废液。

1.3.1.6 加入200μL试剂盒中的L-Desulphonation Buffer室温(20-30°C)孵育15-20min，孵育完成后，12000rpm离心1min，弃废液。

1.3.1.7 加入200μL试剂盒中的M-Wash Buffer，12000rpm离心1min，弃废液，重复两次。

1.3.1.8 将吸附柱放回收集管中，12000rpm离心2min，倒掉废液。将吸附柱开盖置于室温放置2-5min，以彻底晾干吸附材料中残余的漂洗液。

1.3.1.9 将吸附柱转入一个干净的离心管中，向吸附膜的中间部位悬空滴加20μL洗脱缓冲液TE洗脱，室温放置2-5min，12000 rpm离心1min。

1.3.1.10 将收集管中的液体重新加回吸附柱中，室温放置2-5min，12000rpm离心1min，将收集有转化后DNA的离心管-20℃保存（转化后DNA尽快使用）。

1.3.2 DNA预处理

1.3.2.1 PCR仪提前95℃预热，热盖温度105℃。

1.3.2.2 取转化后的片段化DNA放入0.2ml的PCR管中，加入低浓度乙二胺四乙酸TE缓冲液（Low EDTA TE）稀释总体积到15μL。

1.3.2.3 将PCR管放入PCR仪中，进行95℃孵育2min后，立即放置到冰上，静置2min。

1.3.3 加T7接头

1.3.3.1 PCR仪提前37℃预热，热盖温度105℃。

1.3.3.2 按照表4配置反应体系，表格中的试剂均来自ACCEL-NGS® METHYL-SEQDNA LIBRARY KIT试剂盒（Swift Biosciences公司生产）。

表4：反应试剂列表

1.3.3.3 加25μL试剂到冰上放置的预处理DNA样本PCR管中，使用移液器进行吹打混匀，瞬时离心。

1.3.3.4 将PCR管置于PCR仪中，进行反应，条件如表5所示。

表5：反应条件

1.3.4二链合成反应（Second strand synthesis reaction）

1.3.4.1 PCR仪提前98℃预热，热盖温度105℃。

1.3.4.2 按照表6配置反应试剂，表格中的试剂均来自剂来自ACCEL-NGS®METHYL-SEQ DNA LIBRARY KIT试剂盒（Swift Biosciences公司生产）。

表6：反应试剂列表

1.3.4.3 加44μL表6试剂到上一步反应体系中，使用移液器进行吹打混匀，瞬时离心。

1.3.4.4 将PCR管置于PCR仪中，进行二链合成反应，反应条件如表7所示。

表7：二链合成反应条件

1.3.4.5 提前将纯化磁珠从4℃取出，室温平衡半小时。

1.3.4.6 待上一步反应结束后，在产物中加入101μL磁珠，吹打混匀。

1.3.4.7 室温静置5min，置于磁力架上至液体澄清，弃去上清。

1.3.4.8 加入200μL 80％乙醇孵育30sec后弃去。这里80%乙醇现用现配。重复一次200μL 80％乙醇清洗步骤。

1.3.4.9 用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。

1.3.4.10 从磁力架取下离心管，加入16μL超纯水，振荡混匀。室温孵育2min。

1.3.4.11 短暂离心，置于磁力架上至液体澄清，将15μL样本转入新的离心管中。

1.3.5 加T5接头

1.3.5.1 按照表8配置反应试剂，表格中的试剂均来自ACCEL-NGS® METHYL-SEQ DNALIBRARY KIT试剂盒（Swift Biosciences公司生产）。加15μL反应体系到上一步的样本中，使用移液器进行吹打混匀，瞬时离心。

表8：反应试剂列表

1.3.5.2 将PCR管置于PCR仪中，按表9的条件进行PCR反应。

表9：PCR反应的条件

1.3.5.3 提前将纯化磁珠从4℃取出，室温平衡半小时。

1.3.5.4 连接反应结束后，加入36μL磁珠，吹打混匀。

1.3.5.5 室温静置5min，置于磁力架上至液体澄清，弃去上清。

1.3.5.6 加入200μL 80％乙醇孵育30sec后弃去。这里80%乙醇现用现配。重复一次200μL 80％乙醇清洗步骤。

1.3.5.7 用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。

1.3.5.8 从磁力架取下离心管，加入20μL超纯水，振荡混匀。室温孵育2min。

1.3.5.9 短暂离心，置于磁力架上至液体澄清，将20μL样本转入新的离心管中。

1.3.6 扩增

1.3.6.1 按照表10配置反应试剂，加30μL反应体系到上一步的样本中，使用移液器进行吹打混匀，瞬时离心，表格中的试剂来自ACCEL-NGS® METHYL-SEQ DNA LIBRARY KIT试剂盒（Swift Biosciences公司生产）。

表10：反应试剂列表

1.3.6.2 将PCR管置于PCR仪中，按表11的条件进行PCR反应。

表11：PCR反应的条件

1.3.6.3 提前将纯化磁珠从4℃取出，室温平衡半小时。

1.3.6.4 连接反应结束后，加入60μL磁珠，吹打混匀。

1.3.6.5 室温静置5min，置于磁力架上至液体澄清，弃去上清。

1.3.6.6 加入200μL80％乙醇孵育30sec后弃去。这里80%乙醇现用现配。重复一次200μL 80％乙醇清洗步骤。

1.3.6.7 用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。

1.3.6.8 从磁力架取下离心管，加入50μL超纯水，振荡混匀。室温孵育2min。

1.3.6.9 短暂离心，置于磁力架上至液体澄清，将50μL样本转入新的离心管中。

1.4 文库捕获

1.4.1 混合文库：

按每个捕获总量1ug捕获。向上述体系中加入杂交试剂，震荡混匀，短暂离心。

1.4.2 用封口膜封住EP管，放入真空离心浓缩仪中蒸干（60℃，约20min-1hr）。注意随时查看是否已蒸干。

1.4.3 DNA变性：

样本完全蒸干后，每个捕获中加入7.5μL 2×Hybridization Buffer (vial5)和3μLHybridization Component A (vial 6)，震荡混匀，短暂离心，置于95℃变性10min。该步骤中的两种试剂都来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。

1.4.4 文库与探针杂交：

1.4.4.1 取出探针，短暂离心。

1.4.4.2 短暂离心，将变性的DNA（始终保持在95℃）快速转移至含有探针的PCR管中，震荡混匀，短暂离心。

1.4.4.3 置于PCR仪中，47℃杂交。

1.4.5 配制纯化试剂

1.4.5.1 一个捕获所需纯化试剂的配制方法如表12所示，根据捕获的个数按下表配制缓冲液。表格中试剂均来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。

表12：捕获所需纯化试剂的配制试剂列表

1.4.5.2 孵育捕获磁珠（Capture Beads）和清洗缓冲液（Wash Buffer）工作液：

l Capture Beads使用前须室温平衡30min。

l Wash Buffer使用前须47℃孵育2hr。

1.4.6 杂交后纯化

1.4.6.1 每个捕获分装100μL捕获磁珠，将100μL捕获磁珠置于磁力架上至液体澄清，弃去上清。

1.4.6.2 加入200μL 1×Bead Wash Buffer (vial 7)，震荡混匀，置于磁力架上至液体澄清，弃去上清，重复两次。再次加入100μL 1×Bead Wash Buffer (vial 7)，震荡混匀，置于磁力架上至液体澄清，彻底弃去上清。磁珠预处理完成，立即进行下一步试验。

1.4.6.3 将捕获过夜的杂交液体转入清洗好的磁珠中，移液器吹打十次。置于PCR仪中47℃孵育45min（PCR热盖温度设为57℃），每隔15min震荡一次保证磁珠悬浮。1×BeadWash Buffer (vial 7)来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。

1.4.7使用SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）进行清洗

1.4.7.1 孵育完成后，每管加入100μL 47℃预热的1×Wash Buffer I(vial 1)，震荡混匀。置于磁力架上至液体澄清，弃去上清。该步骤至1.4.7.6中所使用的试剂均来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。

1.4.7.2 加入200μL 47℃预热的1×Stringent Wash Buffer (vial 4)，移液器吹打十次混匀。47℃孵育5min，置于磁力架上至液体澄清，弃去上清。

1.4.7.3 加入200μL 47℃预热的1×Stringent Wash Buffer (vial 4)，移液器吹打十次混匀。47℃孵育5min，置于磁力架上至液体澄清，弃去上清。

1.4.7.4 加入200μL室温放置的1×Wash Buffer I (vial 1)，振荡2min，短暂离心，置于磁力架上至液体澄清，弃去上清。

1.4.7.5 加入200μL室温放置的1×Wash Buffer II (vial 2)，震荡1min，短暂离心，放置磁力架上至液体澄清，弃去上清。

1.4.7.6 加入200μL室温放置的1×Wash Buffer III (vial 3)，震荡30sec，短暂离心，放置磁力架上至液体澄清，弃去上清。

1.4.7.7 向离心管中加入36μL超纯水洗脱，震荡混匀，进行下一步扩增试验。

1.4.8 PCR反应

1.4.8.1 根据捕获个数，按照表13配制混合液，震荡混匀。表中试剂均来自SeqCap®Hyb and Wash Kit试剂盒（Roche公司生产）。

表13：混合液的配置试剂列表

1.4.8.2 短暂离心，将混合液分装至PCR管中，每管30μL。每个捕获样本分为两管进行PCR扩增，每管样本20μL。

1.4.8.3 上述样本转入PCR反应中，震荡混匀，短暂离心。

1.4.8.4 置于PCR仪上，按表14的条件进行PCR反应。

表14：PCR反应的条件

1.4.9 扩增后纯化

1.4.9.1 取出纯化磁珠，室温平衡30min备用。

1.4.9.2 取180μL纯化磁珠于1.5mL离心管中，加入100μL扩增后的捕获DNA文库，振荡混匀，室温孵育15min。

1.4.9.3 置于磁力架上至液体澄清，弃去上清。

1.4.9.4 加入200μL 80％乙醇孵育30sec后弃去。这里80%乙醇现用现配。

重复一次200μL 80％乙醇清洗步骤。

1.4.9.5 用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。

1.4.9.6 从磁力架取下离心管，加入120μL超纯水，振荡混匀。室温孵育2min。

1.4.9.7 短暂离心，置于磁力架上至液体澄清，将捕获样本转入新的离心管中。

1.5 文库混库和测序

将上述每个捕获按照数据量比例计算混库质量，按照数据量比例将不同捕获混合成一个样本。加入Phix文库混合成上机样本，进行测序。Phix为一种噬菌体，能够改善碱基不平衡，其作为参考品也可以对测序质量进行评估。

将下机FASTQ文件处理为各模块和软件可使用的输入文件

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

2.1 去接头

调用Trimmomatic-0.36将每一对FASTQ文件都作为配对序列（paired reads）进行去除接头和低质量碱基处理，生成去接头后的FASTQ文件。具体，在切除接头序列后，切去剩余部分开头和结尾处碱基质量低于20的碱基，从reads的5’端开始，以大小为5的窗口进行划窗计算平均质量，如果窗口内平均碱基质量低于20，则切除该窗口，并要求切除后剩余碱基数量超过75。

2.2 比对

调用Bismark-v0.19.0将去接头后的FASTQ文件作为paired reads比对到hg19人类参考基因组序列和lambda DNA参考基因组序列，生成初始Bam文件和比对报告。

2.3 去重

调用Bismark-v0.19.0的deduplicate模块，对初始Bam文件进行去重复处理，并生成去重后的Bam文件和去重后的结果报告。

2.4 排序标记

调用SAMtools-1.3的sort模块，对去重后的Bam文件进行排序，生成排序后的Bam文件；并调用Picard-2.1.0（一种处理高通量测序数据的工具，可用于处理sam/bam等比对结果文件）的AddOrReplaceReadGroups模块，对排序后的Bam文件进行标记分组。

2.5 筛选

调用BamUtil-1.0.14的clipOverlap模块对标记分组后的Bam文件进行筛选，对Bam文件与配对序列的碱基存在重叠且配对序列比对到参考序列负链的reads进行cigar值转换处理，生成Bam文件；并调用SAMtools-1.3 view对去除重叠序列的Bam文件的比对质量（用于量化比对到错误位置的可能性，值越高表示可能性越低）进行过滤，要求比对质量超过20，生成最终Bam文件。Cigar值报告了Bam文件中每条read的相关比对信息。

2.6 建立索引

调用SAMtools-1.3的index模块对最终生成的Bam文件建立索引，生成与最终Bam文件配对的bai文件。

2.7 数据统计

调用FASTQC-0.11.3统计去接头前后的FASTQ文件的碱基数据量、reads数据量以及碱基分布等；统计比对过程中生成的人类参考基因组比对报告中Total C碱基含量、甲基化C碱基含量，非甲基化C碱基含量、CpG context和non-CpG context中甲基化的C碱基含量、CpG context和non-CpG context中非甲基化的C碱基含量；调用Bedtools- v2.26.0的intersect模块统计最终生成的Bam文件中目标区域碱基数量、目标区域捕获的数据量和比例；调用SAMtools-1.3统计最终生成的Bam文件的测序深度、平均测序深度以及不同测序深度下目标区域捕获的碱基数量和比例。

用传统方法直接对ctDNA样本的CpG甲基化水平进行鉴定

利用BisSNP软件对最终生成的Bam文件进行处理：首先调用BisSNP-0.82.2的BisulfiteCountCovariates和BisulfiteTableRecalibration模块进行碱基质量校正，生成校正后的csv文件和Bam文件；之后，利用BisulfiteGenotyper模块和目标区域Bed文件，鉴定待测样本的SNP位点和CpG位点，生成SNP和CpG的原始VCF文件；最后，根据生成的VCF文件，调用VCFpostprocess模块对CpG位点进行过滤，得到最终的CpG位点及其甲基化水平。

用本发明的方法对ctDNA样本CpG甲基化水平进行鉴定

采用最终生成的Bam文件作为输入文件，以non-CpG context的C-T转化率95%为最低要求，调用本发明的reads水平过滤模块，对Bam文件进行逐行读取，判断每条reads的non-CpGcontext是否满足C-T转化率的最低要求，筛选出符合要求的reads，并生成过滤后的Bam文件；再采用过滤后的Bam文件和BisSNP-0.82.2软件鉴定出的CpG位点为输入文件，要求每条reads上至少包含3个CpG位点，调用甲基化水平预测模块，过滤掉Bam文件中不符合要求的reads，然后计算每个CpG位点的甲基化水平。

比较传统方法和本发明方法的甲基化水平预测结果

分别利用不同的甲基化水平预测方法，对6对重复样本的甲基化水平进行样本间相关性比较，结果如下：

5.1 不同方法，各重复样本间CpG位点甲基化水平同时小于1的位点甲基化水平预测结果一致性如表15所示，其中，表15的Sample列表示用于计算相关性的配对重复样本，non-C-T-BisSNP列表示未经过C-T转化率过滤，利用BisSNP-0.82.2软件计算甲基化水平的方法，即传统方法的相关系数；C-T-BisSNP列表示经过C-T转化率过滤，利用BisSNP-0.82.2软件计算甲基化水平的方法的相关系数；C-T-estimate列表示本发明方法的相关系数。

表15：不同方法下各重复样本甲基化水平预测结果相关系数列表（全部位点）

5.2 不同方法，各重复样本间CpG位点甲基化水平同时小于0.02的位点甲基化水平预测结果一致性如表16所示，其中，表16的Sample列表示用于计算相关性的配对重复样本，non-C-T-BisSNP列表示未经过C-T转化率过滤，利用BisSNP-0.82.2软件计算甲基化水平的方法，即传统方法的相关系数；C-T-BisSNP列表示经过C-T转化率过滤，利用BisSNP-0.82.2软件计算甲基化水平的方法的相关系数；C-T-estimate列表示本发明方法的相关系数。

表16：不同方法下各重复样本甲基化水平预测结果相关系数列表（低甲基化水平位点）

从表中可以看出，相较于non-C-T-BisSNP和C-T-BisSNP方法，本发明中增加的reads水平过滤模块和甲基化水平预测模块均使重复样本间低甲基化水平的相关性得到了提高，更适用于ctDNA的甲基化水平预测。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各程序模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的程序模块完成，即将装置的内部结构划分成不同的程序单元或模块，以完成以上描述的全部或者部分功能。实施例中的各程序模块可以集成在一个处理单元中，也可是各个单元单独物理存在，也可以两个或两个以上单元集成在一个处理单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件程序单元的形式实现。另外，各程序模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。

图3是本发明一个实施例中提供的终端设备的结构示意图，如所示，该终端设备200包括：处理器220、存储器210以及存储在存储器210中并可在处理器220上运行的计算机程序211，例如：基于目标区域捕获测序的ctDNA甲基化水平预测关联程序。处理器220执行计算机程序211时实现上述各个基于目标区域捕获测序的ctDNA甲基化水平预测方法实施例中的步骤，或者，处理器220执行计算机程序211时实现上述基于目标区域捕获测序的ctDNA甲基化水平预测装置实施例中各模块的功能。

终端设备200可以为笔记本、掌上电脑、平板型计算机、手机等设备。终端设备200可包括，但不仅限于处理器220、存储器210。本领域技术人员可以理解，图3仅仅是终端设备200的示例，并不构成对终端设备200的限定，可以包括比图示更多或更少的部件，或者组合某些部件，或者不同的部件，例如：终端设备200还可以包括输入输出设备、显示设备、网络接入设备、总线等。

处理器220可以是中央处理单元(Central Processing Unit，CPU)，还可以是其他通用处理器、数字信号处理器 (Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器220可以是微处理器或者该处理器也可以是任何常规的处理器等。

存储器210可以是终端设备200的内部存储单元，例如：终端设备200的硬盘或内存。存储器210也可以是终端设备200的外部存储设备，例如：终端设备200上配备的插接式硬盘，智能TF存储卡(Smart Media Card，SMC)，安全数字(Secure Digital，SD)卡，闪存卡(Flash Card)等。进一步地，存储器210还可以既包括终端设备200的内部存储单元也包括外部存储设备。存储器210用于存储计算机程序211以及终端设备200所需要的其他程序和数据。存储器210还可以用于暂时地存储已经输出或者将要输出的数据。

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

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。

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

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

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

集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序211发送指令给相关的硬件完成，的计算机程序211可存储于一计算机可读存储介质中，该计算机程序211在被处理器220执行时，可实现上述各个方法实施例的步骤。其中，计算机程序211包括：计算机程序代码，计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。计算机可读存储介质可以包括：能够携带计算机程序211代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器 (ROM，Read-Only Memory)、随机存取存储器(RAM，RandomAccess Memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，计算机可读存储介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如：在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。

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