CN113249453A - 一种检测拷贝数变化的方法 - Google Patents

Abstract

本发明公开了一种检测拷贝数变化的方法，包括以下步骤：样本基于2X设计的多重PCR建库测序，获得带UMI标签的样本测序数据，进行质控后，再与参考基因组进行比对，获得比对数据；样本包括测试样本和对照样本；计算样本测序数据中的每个扩增子的测序序列数目，得到矩阵文件，矩阵文件中每列代表一个样本，每行代表一个扩增子的测序序列数目；对矩阵文件进行校正；校正后的矩阵文件，将矩阵文件每行取平均数或中位数，用测试样本的某个扩增子所在行的平均数或中位数，除以对照样本对应于测试样本的某个扩增子所在行的平均数或中位数，得到比值；根据比值结果判断目标区域CNV情况。该方法降低了数据噪音，提高了CNV检测的准确性。

Description

一种检测拷贝数变化的方法

技术领域

本发明涉及生物信息学与精准医学全基因组变异检测技术领域，具体涉及一种基于2X设计多重PCR高通量测序的检测拷贝数变化的方法。

背景技术

拷贝数变异(copy number variation，CNV)是一种与基因功能和人类疾病发生密切相关的遗传突变，因此CNV的研究目前在新生儿遗传性疾病检测、遗传疾病的诊断和某些常见病的辅助诊断方面应用较为广泛。对于某些疾病，只关注特定基因即可达到同样效果。靶向测序技术可以将感兴趣的基因组区域富集出来测序，单个样本测序数据产出少且分析速度较快，因此更能经济高效地发挥NGS技术的优势，广泛应用到临床检测、健康筛查等众多领域。另外，靶向测序可以对目标区域进行深度测序，增加了目标区域内遗传变异检测的灵敏度和准确性。

靶向测序的方法主要分为两种：杂交捕获测序和多重扩增子测序。多重PCR(multiplex PCR)，又称多重引物PCR或复合PCR，多重扩增子测序即针对感兴趣的目标区域，设计多重PCR引物进行扩增富集并进行测序的技术。通常适用于检测几十到几千个位点，或几十kb以下的区域。杂交捕获测序，目前应用的主要是液相杂交捕获测序，即基于碱基互补配对原理，设计合成核酸探针，对DNA文库进行基于液相环境的目标区域杂交富集，并进行测序。但液相杂交捕获操作难度较高、操作时间较长、且容易受到探针捕获效率的影响，因此扩增子测序相比来说，更适合非专业技术人员操作。多重PCR作为一种快速构建靶向测序文库的方法，由于其高效性、系统性和经济简便性在目前的临床基因检测及研究领域中发挥越来越大的作用。

多重PCR目前的技术局限性包括等位基因脱扣（Allele drop-out, ADO），即杂合子的2个等位基因中的1个优势扩增，而另1个完全扩增失败。为解决ADO问题，有文献提供2X设计的解决方法（如图1所示），即在目标区域设计两对引物扩增，即使其中一对引物受到ADO影响，另一对引物还可以将正常等位信息扩增出来。该方案虽然可以避免检测点突变时的漏检误检，但对于拷贝数检测还是会造成影响，当其中一对引物受到ADO影响，基于测序深度(read depth) 或序列数目(read count) 算法检测CNV时会有一定数据波动，引入噪音，造成检测假阳性。

发明内容

本发明的目的在于，提供一种检测结果稳定、准确的检测拷贝数变化的方法。

本发明检测拷贝数变化的方法，是基于2X设计的多重PCR扩增体系，利用带分子标签（Unique Molecular Identifiers, UMI）的高通量测序数据，对每对引物的扩增效率进行评估，并对测试样本扩增子CNV进行检测。

UMI（unique molecular identifiers）是一种分子条形码，可以在测序过程中错误校正，提高准确性。这些分子条形码均为短序列，可特异性地标记样本文库中的每个分子，UMI测序可以降低假阳性变异检出的概率，同时能提高变异检测的灵敏度。由于起始材料中的每个核酸都有唯一的分子条形码，因此，生物信息学软件可以高度精确地过滤出重复的read和PCR错误，报告唯一read，从而在最终数据分析之前消除已识别的错误。UMI是一段随机化或特定的核苷酸短序列，在建库的过程中通过连接的方式导入，从而如同分子条形码一般特异性地标识每个模板，用以在高通量测序的海量数据中区分同一来源的DNA片段。

UMI碱基，即由碱基短序列组成的分子条形码。

本发明技术方案详述如下：

一种检测拷贝数变化的方法，包括以下步骤：

（1）将样本DNA加上UMI碱基，进行基于2X设计的多重PCR建库测序，获得带UMI标签的样本测序数据，对样本测序数据进行质控后，再与参考基因组进行比对，获得比对数据；所述样本包括测试样本和对照样本；

（2）根据样本测序数据及引物坐标信息计算样本测序数据中的每个扩增子的测序序列数目，得到矩阵文件，矩阵文件中每列代表一个样本，每行代表一个扩增子的测序序列数目；对矩阵文件进行校正；

校正后的矩阵文件，将矩阵文件每行取平均数或中位数，用测试样本的某个扩增子所在行的平均数或中位数，除以对照样本对应于测试样本的某个扩增子所在行的平均数或中位数，得到比值；

（3）根据比值结果进行以下判断：

根据测试样本的扩增子的分析结果判断该扩增子目标区域CNV情况：

其中，第1X和第2X表示同一扩增区域中的两对引物，不分先后顺序。

可选或优选的，上述检测拷贝数变化的方法中，步骤（2）中对矩阵文件进行校正的方法如下：

先按行进行几何平均：

得到如下矩阵：

再计算列的中位数Corr：

校正

其中，M为校正前矩阵文件，M*为校正后矩阵文件，a代表每个样本在每个引物扩增区间计算得到的测序深度，g代表每行即每个引物扩增区间所有样本取几何平均数，Test为待测样本，Ctrl为对照样本，b代表每个样本在每个引物扩增区间的测序深度校正后的数值，med为每列取中位数即每个样本所有引物扩增区间的深度取中位数，mean为取平均数，即对所有对照样本校正后的深度按行（每个引物扩增区间）计算平均数。

可选或优选的，上述检测拷贝数变化的方法中，步骤（1）中对样本测序数据进行质控，包括数据量合格、平均测序深度达到要求、数据质量Q20>90%，Q30> 85%。

数据量合格，是指测序数据量达到1G以上。

平均测序深度达到要求，是指平均测序深度达到3000X以上。

Q20和Q30：测序数据中每个碱基都有对应的质量值，质量值是Q20，则错误识别的概率是1%，即错误率1%，或者正确率是99%；质量值是Q30，则错误识别的概率是0.1%，即错误率0.1%，或者正确率是99.9%。根据前面的解释，Q20与Q30则表示质量值≥20或30的碱基所占百分比。Q20>90%：质量值≥20的碱基所占百分比需大于90%；Q30> 85%：质量值≥30的碱基所占百分比需大于85%。

可选或优选的，上述检测拷贝数变化的方法中，步骤（1）中样本测序数据是使用sentieon软件的UMI的流程得到consensus 的fastq文件，所述比对数据是fastq文件再经过bwa 软件和参考基因组进行比对，经过sentieon 的UMI 处理后得到consensus的bam；

步骤（2）中，根据样本测序数据及引物坐标信息计算样本测序数据中的每个扩增子的测序序列数目的具体方法如下：

使用python中的pysam模块对经过sentieon 的UMI 处理后得到consensus的bam进行处理，得到每条序列比对到参考基因组的起始终止坐标，将序列坐标与引物的坐标一致的进行计数，得到该引物扩增的序列数目，其他引物以此计算方式进行扩增区域序列数目的计算。

当面对多个目标区域拷贝数变化检测时，每个目标区域都会设计2X引物，每个区域也按照类似的处理方式进行检测。

名词释义：

测序序列数目：read count，利用二代测序技术对基因或者转录本进行测序，测序测到的每条序列为read，通过统计某一区域测到的reads 数目即read count。

与现有技术相比，本发明具有如下有益效果：

以往基于2X设计的多重扩增高通量测序，虽然能避免ADO对于点突变检测的影响，但对于CNV检测来说引入了很大的噪音。本发明的方法根据2X设计得到的高通量数据，引物扩增坐标计算每个扩增子的测序序列 count，并以历史阴性样本作为对照，评估了每个扩增子的扩增效果，对待测样本CNV进行检测，一定程度上降低了数据的噪音，提高了对于2X设计下的CNV检测的准确性。

正常处理CNV流程是根据设计的目标区域计算深度（例如基因的每个exon），2X设计下，使用设计时的目标区域坐标来计算时会受到深度干扰，所以从每对引物扩增的位置去计算，先对引物扩增效果进行评估，再进行CNV的检测。

多重体系下，正常扩增引物应该是和引物设计时的起始终止坐标一致，即测序得到的成对的序列比对到参考基因组上的起始终止坐标和引物设计时的坐标一致，不会出现捕获测序中每条序列比对坐标不一致。参见图7。

附图说明

图1为公开文献中提供的解决ADO问题的2X设计方法示意图；

图2为实施例1第一部分构建样本的基因组比对与校正后文件流程示意图；

图3为实施例1第二部分构建扩增区域文件，根据扩增区域计算每个扩增子深度流程示意图；

图4为实施例1第三部分基于扩增子检出结果对CNV进行检测流程示意图；

图5为正常样本目标区域2X扩增结果；

图6为发生了ADO的样本目标区域2X扩增结果；

图7为多重测序序列比对结果和捕获测序序列比对结果。

具体实施方式

下面结合具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。

实施例1 基于2X设计的多重PCR扩增测序数据进行引物扩增效率的评估和后续CNV的检测

第一部分：构建样本的基因组比对与校正后文件（umi_consensus.bam）

基于UMI技术原理，将样本DNA加上UMI碱基，进行基于2X设计的多重PCR建库测序，经过高通量测序获得带UMI标签的样本测序数据。样本测序数据是若干带标签的扩增子序列的集合数据库。

首先对样本测序数据进行基本质控，数据质控包括数据量合格、平均测序深度达到要求、数据质量Q20>90%,Q30> 85%。

Q20和Q30：测序数据中每个碱基都有对应的质量值，质量值是Q20，则错误识别的概率是1%，即错误率1%，或者正确率是99%；质量值是Q30，则错误识别的概率是0.1%，即错误率0.1%，或者正确率是99.9%；

然后使用sentieon软件（NGS基因数据分析加速软件），根据该软件的UMI流程（Sentieon UMI consensus 模块）提取UMI信息，处理原始的fastq文件，联合sentieon软件的bwa模块得到consensus 的fastq文件。然后利用sentieon软件的bwa模块和参考基因组进行比对，获得比对数据，再使用sentieon软件中的UMI处理模块对比对后的bam文件进行后续处理，得到最终的bam文件，即比对数据。包括测试样本bam文件和对照样本bam文件。

本部分方法流程如图2所示。

第二部分：构建扩增区域文件，根据扩增区域计算每个扩增子深度，详细步骤如下：

1、根据引物设计文件，统计每对扩增子的起始坐标与终止坐标。

引物设计文件，是指在多重PCR建库测序时，使用的每个目标扩增子的扩增引物。因基于2X，所以是每个目标扩增子设计两对扩增引物。

引物设计文件中包括正向引物的起始坐标和终止坐标，反向引物的起始坐标和终止坐标。以正向引物的起始坐标和反向引物的终止坐标构建得到引物扩增区域文件。

2、获得测试样本和对照样本扩增子位置信息和测序序列数目

根据第一部分得到的测试样本bam文件中，每对引物扩增的序列可拼接得到一条序列，扩增片段较大的，则在bam中显示为paired 测序序列。

拼接的序列比对到参考基因组的起始位置和终止位置即为引物扩增的正向 引物的5’端坐标和反向引物的5’端坐标。

Paired 测序序列比对到参考基因组的起始位置和终止位置也为引物扩增的正向引物的5’端坐标和反向引物的5’端坐标。

通过步骤1中引物设计文件中的正反向引物起始终止坐标信息（已知的引物设计的坐标信息），可根据正向引物的5’端坐标和反向引物的5’端坐标去计算每对引物扩增出的扩增子测序序列数目，将测序序列比对的坐标（该测序序列read所对应的扩增子在参考基因组中的位置信息）与引物坐标（引物设计文件中针对已知目标扩增子设计的引物的位置信息）进行比较，计算比对到该对引物对应扩增子的测序序列数目，即可计算该扩增子对应的测序序列数目，得到样本的深度文件。

每对扩增子的起始终止坐标信息为输入文件，根据第一部分得到的bam文件，每个测试样本可得到一个每对扩增子测序序列数目的文件。

3. 选取历史阴性样本（对照样本），使用上述计算方法计算对照样本深度文件。

4. 将测试样本的测序序列数目文件与对照样本的文件按照相同扩增引物位置进行合并，可得到一个矩阵文件。矩阵文件中每列代表一个样本，每行代表一个扩增子的测序序列数目文件。对最终得到的包含测试样本和对照样本的深度统计文件（即矩阵文件）进行数据校正，具体校正方式如下：

先按行(基因)几何平均：

得到如下矩阵：

然后计算列的中位数：

校正

a代表每个样本在每个引物扩增区间计算得到的测序深度，g代表每行即每个引物扩增区间所有样本取几何平均数，Test 为待测样本，Ctrl为对照样本，

b代表每个样本在每个引物扩增区间的测序深度校正后的数值，

med为每列取中位数即每个样本所有引物扩增区间的深度取中位数，mean为取平均数，即对所有对照样本校正后的深度按行（每个引物扩增区间）计算平均数。

5、根据校正后的矩阵文件 ，将对照样本的矩阵每行取平均数（或中位数），将测试样本的数值除以对照样本的中位数或平均数（如果对照样本取平均数，则测试样本的数值除以对照样本的平均数，如果对照样本取中位数，则测试样本的数值除以对照样本的中位数），可得到对照样本的比值情况。

6、本部分方法流程如图3所示。

第三部分：基于扩增子检出结果对CNV进行检测，详细步骤如下：

1、根据第二部分得到的比值结果，可以判断测试样本中每个扩增子的扩增情况：

2、基于2X 设计的多重PCR，每个目标区域都有2X的引物覆盖。根据每个目标区域的扩增子的分析结果判断该区域的CNV情况。具体结果详见下表，表中1X 和2X引物指同一扩增区域中两对引物，不分先后顺序。即其中一对为正常，另一对引物也正常，即拷贝数没有变化。

针对上述表格中的特殊情况进行说明：

当目标区域中一对引物显示无扩增，另一对引物正常扩增或减半扩增或加倍扩增时，若无扩增的引物的测序序列深度较低或接近于0，可认为该引物失效，以另一对引物检出的结果作为检测CNV的依据。

3、单个目标区域的拷贝数结果可根据上述表格得到，后续分析为segmentation，即将相邻区域拷贝数一致的区域进行合并。根据每个目标区域的引物情况计算每个目标区域的ratio，并将ratio值作log2后，将坐标信息以及log2值作为输入传入DNAcopy R包中使用循环二元切割算法（circular binary segmentation，CBS)做segmentation，并将segmentation的结果进行过滤（去除拷贝数正常的区域）。最后将拷贝数变化的区域输出。

每个目标区域的拷贝数结果可通过上述步骤得到。当目标区域为多个时候，需要将相邻区域拷贝数一致的区域进行合并，即Segmentation。当1X引物和2X引物的扩增情况一致时，目标区域的比值情况为1X和2X引物比值的平均值；当1X引物扩增情况为减半扩增/无扩增，2X扩增情况为正常或加倍扩增时，目标区域的比值情况为2X引物的比值；当1X引物扩增情况为减半扩增/无扩增，2X扩增情况为正常或加倍扩增时，目标区域的比值情况为2X引物的比值。通过上述方式计算目标区域的比值，例如：

将待测样本的比值作为输入，使用循环二元切割算法（circular binarysegmentation，CBS)对拷贝数一致的区域进行合并（DNAcopy R包）得到片段化的结果。正常人的拷贝数为2，当发生杂合缺失，即拷贝数为1，比值则为0.5，发生重复，即拷贝数为3，比值则为1.5，拷贝数正常区域比值则为1。根据片段化后的结果，根据拷贝数变化的阈值将拷贝数正常区域过滤后，即可得到待测样本的CNV检测结果。

例如上表中三个区域，segmentation的结果即chr1: 76190463-76190512，CNV结果为杂合缺失。chr1:76194064-76194183 和chr1: 76199203-76199323 两个区域的比值一致，两个区域可以合并，segmentation的结果为chr1: 76194064-76199323，根据ratio值判断CNV结果为正常。

本部分方法流程如图4所示。

较于其他基于分子标签测序的多重测序，2X引物设计可以避免对点突变的漏检，误检。比如若引物上有杂合突变，导致扩增区域内的杂合突变只扩增出一侧等位基因（allele）的序列，导致检出的突变为野生或纯合，因此2X引物设计可以一定程度上避免这种突变的漏检和误检。但2X引物设计的数据检测CNV时，引物扩增受到引物浓度等其他条件影响，扩增同一区域的扩增子（amplicon）测到的测序序列数目并不是完全一致的，此时若其中一对引物扩增测序序列数目较多，且受到了ADO影响，数据波动的影响会导致CNV检测的不准确性增加。比如：

正常样本目标区域深度为18（目标区域扩增子测序序列reads=18），由于引物扩增受到引物浓度等其他条件影响，扩增同一区域的扩增子（amplicon）测到的测序序列数目并不是完全一致的，一对引物扩增出的测序序列条目数为12，另一对引物扩增出的测序序列条目数为6。如图5所示。

当测试样本发生了ADO，并没有拷贝数变化。由于ADO影响，其中一对引物只扩增出一半的测序序列数目为6，另一对引物正常扩增出了6条测序序列，如图6所示。

此时，按照CNV检测原理，根据测试样本和对照样本该区域的深度进行比较，计算测试样本的比例为（6+6）/18，约为0.667。理想情况下CNV杂合缺失的比值应为0.5附近波动。0.667 的比例无法判断该区域是真实CNV还是受到了数据波动的影响。因此在2X引物设计的多重体系下，需要对每对引物扩增的情况进行判断，才能准确对CNV进行检测。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (4)

1.一种检测拷贝数变化的方法，其特征在于，包括以下步骤：

（1）将样本DNA加上UMI碱基，进行基于2X设计的多重PCR建库测序，获得带UMI标签的样本测序数据，对样本测序数据进行质控后，再与参考基因组进行比对，获得比对数据，所述样本包括测试样本和对照样本；UMI碱基，即由碱基短序列组成的分子条形码；

（2）根据样本测序数据及引物坐标信息计算样本测序数据中的每个扩增子的测序序列数目，得到矩阵文件，矩阵文件中每列代表一个样本，每行代表一个扩增子的测序序列数目；对矩阵文件进行校正；

校正后的矩阵文件，将矩阵文件每行取平均数或中位数，用测试样本的某个扩增子所在行的平均数或中位数，除以对照样本对应于测试样本的某个扩增子所在行的平均数或中位数，得到比值；

（3）根据比值结果进行以下判断：

根据测试样本的扩增子的分析结果判断该扩增子目标区域CNV情况：

其中，第1X和第2X表示同一扩增区域中的两对引物，不分先后顺序。

2.根据权利要求1所述的检测拷贝数变化的方法，其特征在于，步骤（2）中对矩阵文件进行校正的方法如下：

先按行进行几何平均：

得到如下矩阵：

再计算列的中位数Corr：

校正

其中，M为校正前矩阵文件，M*为校正后矩阵文件，a代表每个样本在每个引物扩增区间计算得到的测序深度，g代表每行即每个引物扩增区间所有样本取几何平均数，Test 为待测样本，Ctrl为对照样本，b代表每个样本在每个引物扩增区间的测序深度校正后的数值，med为每列取中位数即每个样本所有引物扩增区间的深度取中位数，mean为取平均数，即对所有对照样本校正后的深度按行计算平均数。

3.根据权利要求1所述的检测拷贝数变化的方法，其特征在于，步骤（1）中对样本测序数据进行质控，包括数据量合格、平均测序深度达到要求、数据质量Q20>90%，Q30> 85%。

4.根据权利要求1所述的检测拷贝数变化的方法，其特征在于，步骤（1）中样本测序数据是使用sentieon软件的UMI的流程得到consensus 的fastq文件，所述比对数据是fastq文件再经过bwa 软件和参考基因组进行比对，经过sentieon 的UMI 处理后得到consensus的bam；

步骤（2）中，根据样本测序数据及引物坐标信息计算样本测序数据中的每个扩增子的测序序列数目的具体方法如下：

使用python中的pysam模块对经过sentieon 的UMI 处理后得到consensus的bam进行处理，得到每条序列比对到参考基因组的起始终止坐标，将序列坐标与引物的坐标一致的进行计数，得到该引物扩增的序列数目，其他引物以此计算方式进行扩增区域序列数目的计算。

Images