一种样本交叉污染的检测方法和装置

技术领域

本发明涉及生物技术领域，具体涉及一种样本交叉污染的检测方法，针对DNA甲基化二代测序样本的单核苷酸多态性等位基因频率进行监控，提供样本是否有交叉污染的判断。

背景技术

重亚硫酸盐甲基化测序(BS-seq,bisulfite sequencing)作为甲基化测序的金标准，以其单碱基分辨率、高通量的特点，其作用在癌症筛查、诊断，以及监控的作用越来越被认识。在高通量二代测序(NGS,Next generation sequencing)检测中，由于多个样本是并行处理，所以在样本储存、制备等过程中导致的相邻样本之间异源DNA交叉污染的风险是难以排除的。而这个风险在早期肿瘤的诊断筛查中，一旦发生后果更加严重，因为早期肿瘤的血液样本中肿瘤组分通常占比很低(<0.001)，血液样本的痕量污染即可造成筛查或者诊断结果错误，而目前的NGS污染检测方法，往往无法达到对<0.001的污染比例的检测敏感度。而且目前常用的NGS的样本污染判定中，通常会在各个批次的样本中设计阳性参考品和阴性参考品，然而在真实的临床实践中，由于对于成本的控制和考量而忽略参考品的设置，这也提高了样本发生交叉污染而不能准确识别的风险。

发明内容

本发明为克服现有技术中的不足，本发明提供了一种样本交叉污染的检测方法和装置，该检测方法和相应的装置可以低成本、高精度的对血液游离DNA样本中是否存在来自其他样本的污染进行判断。

在一方面，本发明提供了一种用于检测甲基化测序中样本污染的单核苷酸多态性(SNP)位点的筛选方法，包括以下步骤：

S1：选取预设人群中频率在0.3～0.7之间的SNP位点；

S2：选取突变方向为腺嘌呤(A)突变到胸腺嘧啶(T)或胸腺嘧啶(T)突变到腺嘌呤(A)的SNP位点；

S3：选取重复区域以外的SNP位点；

S4：选取相互之间物理距离大于预设长度的SNP位点；

任选地，S2和S3的顺序互换。

另一方面，本发明提供了一种用于检测甲基化测序中样本污染的方法，包括以下步骤：

(1)获取对待测样本进行甲基化测序后得到的测序信息；

(2)根据如上述方法筛选到的用于检测甲基化测序中样本污染的SNP位点确定样本污染状态。

另一方面，本发明提供了一种用于检测甲基化测序中样本污染的装置，包括：

测序信息获取模块，被配置成获取对待测样本进行甲基化测序后得到的测序信息；

样本状态确定模块，被配置成根据如上述方法筛选到的用于检测甲基化测序中样本污染的SNP位点确定样本污染状态。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本说明书的实施例，并与说明书一起用于解释本说明书的原理。

图1示出了不同SNPs数量场景下，模拟污染样本在取不同范围的SNP位点时所对应的不同阈值条件的样本污染指数(SCS)，针对单一SCS阈值，模拟污染的掺比梯度分别为：0005pct(0.005％)、001pct(0.01％)、005pct(0.05％)、001pct(0.01％)、05pct(0.5％)和1pct(1％)。其中，图1a为50SNPs的情况；图1b为100SNPs的情况；图1c为200SNPs的情况；图1d为300SNPs的情况；图1e为400SNPs的情况；图1f为500SNPs的情况；图1g为600SNPs的情况。

图2示出了模拟污染样本在取相同范围(AR值波动前20％)的SNP位点时，不同SNPs数量场景下多种掺比的模拟污染样本的样本污染指数。

图3示出了不同cfDNA掺比情况下，污染样本的AR频率分布图，其中横坐标为AR值，纵坐标为SNP位点数量。其中，图3a为无掺比的本底样本；图3b为掺比万分之一的污染样本；图3c为掺比万分之五的污染样本；图3d为掺比千分之一的污染样本；图3e为掺比千分之五的污染样本。

图4示出了参考示例，上方显示的一段预期进行甲基化检测的双链DNA片段，按箭头方向排序，其包含原始上链(CCGGCATGTTTAAACGCT)和原始下链(AGCGTTTAAACATGCCGG)，其中部分假定所有CpG中的胞嘧啶(C)都发生了甲基化，以-mC标识出。上述双链DNA片段经过变性解旋为单链形式后，经过重亚硫酸盐转化处理，原始上链和原始下链中未被甲基化(-mC)修饰的C被转化为尿嘧啶(U)，而甲基化修饰的C则依然保持为C。在随后的PCR扩增过程中，由于尿嘧啶(U)与腺嘌呤(A)互补配对，而DNA的PCR扩增中引入的与腺嘌呤(A)配对碱基为胸腺嘧啶(T)。在PCR扩增中，首先形成了与重亚硫酸盐处理后的带有尿嘧啶(U)的原始上链互补的目标上链互补链(CTOT)，以及与重亚硫酸盐处理后的带有尿嘧啶(U)的原始下链互补的目标下链互补链(CTOB)。在之后的PCR扩增过程中，形成了由原始上链转化的与CTOT互补的目标上链(OT)，以及由原始下链转化的与CTOB互补的目标下链(OB)。对比可知，原始上链和原始下链中的未被甲基化修饰的C在目标上链和目标下链中被T取代，而甲基化修饰的C(以下划线标识出)则保持不变。根据这一特点，可以通过测定重亚硫酸盐转化处理后的C来识别经甲基化修饰的C的数量与位置，从而实现DNA甲基化检测的目的。

具体实施方式

I.定义

在本发明中，除非另有说明，否则本文中使用的科学和技术名词具有本领域技术人员所通常理解的含义。并且，本文中所用的蛋白质和核酸化学、分子生物学、细胞和组织培养、微生物学、免疫学相关术语和实验室操作步骤均为相应领域内广泛使用的术语和常规步骤。同时，为了更好地理解本发明，下面提供相关术语的定义和解释。

术语“SNP”(Single Nucleotide Polymorphism，单核苷酸多态性)主要是指基因组水平上由单个核苷酸的变异所引起的DNA(Deoxyribo Nucleic Acid，脱氧核糖核酸)序列多态性。SNP位点所表现的多态性只涉及到单个碱基的变异，这种变异可由单个碱基的转换或颠换所引起，也可以是碱基的插入或缺失所致。

术语“纯合SNP位点”是指一种SNP位点，在该位点上，所有与参考基因组进行比对的序列上的该位点都显示相同的碱基，且该碱基与参考基因组序列上该位点的碱基不同。例如，若参考基因组序列上某SNP位点上的碱基为G，而所有与参考基因组进行比对的序列上的该SNP位点上的碱基都为A，则该SNP位点则被称为纯合SNP位点。

术语“等位基因”(allele)，是指位于一对同源染色体相同位置上控制同一性状不同形态的基因。当一个生物体带有一对完全相同的等位基因时，则该生物体就该基因而言是纯合的(homozygous)；反之，如果一对等位基因不相同，则该生物体就该基因而言是杂合的(heterozygous)。等位基因各自编码蛋白质产物，决定某一性状，并可因突变而失去功能。

术语“等位基因比率”(AR，allelic ratio)，是指突变型的等位基因与野生型的等位基因的比值，在NGS中，可等价于突变型序列数与野生型序列数的比值。

术语“野生型”(wide type)，是指基因或生物体在自然界中常见的或非突变型的形式。也就是指在野生群体中观察到的最高频率的表型，或具有这种表型的系统、生物或基因。

术语“突变”(mutation)，是指基因的结构发生改变而导致细胞、病毒或微生物的基因型发生稳定的、可遗传的变化过程。

术语“突变型”(mutant)，是指为DNA分子中发生碱基对的增添、缺失或改变，而引起的基因结构的改变。

术语“样本污染指数”(sample contamination score，SCS)，是指用于表征待测样本受到其它样本的交叉污染情况的数值。

II.具体实施方案详述

在一方面，本发明提供了一种用于检测甲基化测序中样本污染的单核苷酸多态性(SNP)位点的筛选方法，包括以下步骤：

S1：选取预设人群中频率在0.3～0.7之间的SNP位点；

S2：选取突变方向为腺嘌呤(A)突变到胸腺嘧啶(T)或胸腺嘧啶(T)突变到腺嘌呤(A)的SNP位点；

S3：选取重复区域以外的SNP位点；

S4：选取相互之间物理距离大于预设长度的SNP位点。

在一些任选的实施方案中，S2和S3的顺序互换。

在一些实施方案中，S1中的SNP位点选取自预设数据库，预设数据库选自gnomAD数据库、1000Genome Project数据库、HapMap数据库、dbSNP数据库中的一个或多个。

在一些优选的实施方案中，预设数据库为gnomAD数据库；

其中，参考基因组为GRCh37/hg19人类参考基因组。

在一些实施方案中，预设人群选自东亚人群、非洲/非裔美国人群、拉丁美洲人群、非芬兰欧洲人群、芬兰欧洲人群、艾希肯纳兹犹太人群或西亚人群。

在一些优选的实施方案中，预设人群为东亚人群。

在一些实施方案中，S4中所述预设长度选自0.4～1Mb。

在一些优选的实施方案中，预设长度为1Mb。

另一方面，本发明提供了一种用于检测甲基化测序中样本污染的方法，包括以下步骤：

(1)获取对待测样本进行甲基化测序后得到的测序信息；

(2)根据如上述方法筛选到的用于检测甲基化测序中样本污染的SNP位点确定样本污染状态。

在一些实施方案中，步骤(2)包括：

确定待测样本中对应于用于检测甲基化测序中样本污染的SNP位点中的纯合SNP位点。

在一些优选的实施方案中，确定方法为通过计算用于检测甲基化测序中样本污染的SNP位点的突变等位基因比率(AR)。

在一些更优选的实施方案中，纯合SNP位点的AR值小于0.25或大于0.75。

在一些实施方案中，步骤(2)包括：

计算所述纯合SNP位点的样本污染指数，其中，所述样本污染指数为所述纯合SNP位点的标准AR值的中位数；

其中，所述标准AR值通过对所述纯合SNP位点的AR值进行标准化得到，包括：

当AR值小于等于0.5时，该SNP位点的标准AR值等于AR值；

当AR值大于0.5时，该SNP位点的标准AR值为1与所述AR值的差值。

在一些实施方案中，对纯合SNP位点的标准AR值进行排序，选取部分纯合SNP位点的标准AR值计算中位数。

在一些优选的实施方案中，排序为由大到小的排序。

在一些优选的实施方案中，选取由大到小排序后位于前5％、10％、15％、20％或25％的纯合SNP位点的标准AR值以计算样本污染指数。

在一些更优选的实施方案中，选取由大到小排序后位于前20％的纯合SNP位点的标准AR值以计算样本污染指数。

在一些实施方案中，当样本污染指数大于预设阈值时，判断待测样本受到其他样本的交叉污染。

在一些优选的实施方案中，预设阈值选自0.001～0.01。

在一些更优选的实施方案中，预设阈值为0.001。

在一些实施方案中，判断待测样本是否收到其他样本的交叉污染进一步包括以下步骤：

在对样本污染指数和预设阈值进行比较之前，先对样本污染指数与待测样本甲基化测序后的背景噪音进行比较，当样本污染指数大于背景噪音时，再比较样本污染指数和预设阈值。

另一方面，本发明提供了一种用于检测甲基化测序中样本污染的装置，包括：

测序信息获取模块，被配置成获取对待测样本进行甲基化测序后得到的测序信息；

样本状态确定模块，被配置成根据如上述方法筛选到的用于检测甲基化测序中样本污染的SNP位点确定样本污染状态。

为了达到清楚和简洁描述的目的，本文中作为相同的或分开的一些实施方案的一部分来描述特征，然而，将要理解的是，本发明的范围可包括具有所描述的所有或一些特征的组合的一些实施方案。

实施例

数据准备：SNP筛选

评估样本是否受到交叉污染的第一步是从突变数据库筛选合适的SNP位点。首先根据hg19 gnomAD数据库结果筛选东方人群中频率在0.3～0.7之间的多核苷酸多态性(SNP)位点，人群频率过低或者过高均不适合作为评估位点。由于本发明的应用场景是甲基化测序数据，在DNA被重亚硫酸盐处理后，原始上链和原始下链非甲基化的C会被转换成T，而对应的互补链G转变成A。基于此，这一步筛选的SNPs的突变方向需要满足A->T或者T->A。如果该SNP位于重复区域，这个位点在测序过程中，更容易被测错，所以同样过滤去除位于重复区域的SNP。最后挑选的SNP相互之间物理距离大于0.4M。

基于监控的SNP位点搭建算法

1.监控SNP位点等位基因比率(AR，allelic ratio)的计算

在对于实际用于污染的监控的样本，计算实施例1中所筛选出的所有SNPs的AR值。AR的计算方式是含有突变SNP的reads数目除以覆盖该SNP的所有读长(reads)数目。理论上，纯合野生型SNP位点的AR为0，纯合突变型SNP位点的AR为1，杂合SNP位点的AR为0.5。在真实实验中，由于存在测序错误等一系列误差，实际上纯合SNP的AR在0或者1附近会出现细微的波动现象，杂合SNP的AR会在0.5附近波动。根据样本的测序结果，利用软件samtoolsmpileup模块对SNPs的检测，以及对实施例1中所筛选出的所有SNPs进行AR统计。当一个样本受到50％以下的污染，SNP的AR值在<0.25或者>0.75的范围波动。如果某个SNP的AR位于0.25～0.75之间，则认为该SNP在这个样本中为杂合SNP，不纳入后续计算。同时，用于后续计算的SNP位点满足测序深度大于50X。保留所有被认为是纯合SNP的位点进入下一步计算。

2.根据SNP位点等位基因频率评估样本整体受到污染的可能性

针对上述步骤1中所保留的纯合SNPs的AR进行标准化，以a值为标准化后的AR值：

1)如果SNP为纯合野生型，即AR≤0.5，则a＝AR；

2)如果SNP为纯合突变型，即AR>0.5，则a＝1-AR。对所有a值进行从大到小的排序，挑选位于前20％的a值计算中位数，设为样本污染指数(SCS，sample contaminationscore)。也就是挑选AR值波动最大的20％的SNP去计算SCS。公式如下：

SCS＝median(r)-

r为纯合SNP中a值位于前20％的a值列表。

在理想情况下，一个样本如果没有受到其他样本的交叉污染，那么它的SCS为0。但是由于测序错误等背景噪音的存在，SCS会在0附近出现细微的波动。如果一个样本的SCS值大于背景噪音的阈值，那么就判断这个样本很有可能受到了其他样本的交叉污染；否则认为这个样本受到其他样本污染的概率很低。综合考虑测序错误，实验误差等影响，将背景噪音的阈值设定为0.001。

实施例1：模拟数据验证

为了验证本方法学的可行性，对真实样本A和真实样本B的SNP位点进行人为修改，从而产生多个样本A、B之间差异纯合位点。例如SNP_1位点在样本A中的基因型是ref/ref，相同位点在样本B中基因型被修改为alt/alt型，那么则认为SNP_1位点在样本A与B之间是差异纯合位点。修改后为模拟真实的应用场景，保留了万分之四的测序背景噪音。构建多个梯度的reads数量掺比，将样本A的reads片段掺入样本B中，梯度为十万分之五、万分之一、万分之五、千分之一、千分之五和百分之一，分别对应图1横坐标0005pct、001pct、005pct、001pct、05pct和1pct。每个梯度掺比重复模拟了5次，每次均为随机抽取reads。同时构建了差异纯合位点个数的七个梯度，分别对应修改了50、100、200、300、400、500和600个SNPs，这些SNPs都随机选取自符合本发明筛选条件的1000个SNP位点(参见表3)。

如实施例2中的步骤，对所得到的SNPs的a值进行计算。

为了确定用于计算SCS所需要取的SNP位点的最佳范围，对构建的SNPs数量分别为50、100、200、300、400、500和600个SNPs的七个场景，分别将计算得到的a值从大到小进行排序，并分别挑选了排序为前5％、10％、15％、20％、25％的SNP位点的a值，对应图1横坐标0.05、0.1、0.15、0.2、0.25。图中虚线对应背景基线，即掺比为0％时本底样本对应的SCS值。由图1可知，在取不同SNPs数量得应用场景中，取前20％的SNP位点a值时，模拟掺比结果的SCS值均能很明显地与背景拉开距离，故用于计算SCS应取的SNP位点的最佳范围是前20％，即挑选AR波动前20％的位点进行SCS计算。

当取AR波动前20％的位点时，根据不同的掺比梯度和使用的不同差异SNPs个数所计算的SCS值箱线图如图2所示。从掺比结果可以推测出，可用于计算的SNPs个数越少，SCS波动越大，结果越不稳定。当应用场景为600个SNPs时，不同掺比样本的SCS值均能明显与基线拉开距离，且当掺比梯度越高时，SCS值上升的越明显。因此在本实施例中选择SNPs数目为600来进行下一步模拟。

当取AR波动前20％的位点时，不同掺比梯度在应用场景SNP数目为600SNPs时计算的SCS值如下表1。本底样本(无掺比)的SCS值为0。设计的5组重复组为分别随机抽取样本A的reads片段掺入样本B中构建得到，从模拟结果可以得出，以能与本底样本拉开明显差距为目的确定阈值，当SCS值的阈值cutoff定在0.001时，可以检出掺比在万分之一(0.01％)的样本污染。因此，当待测样本的SCS值大于0.001时，可以判别待测样本存在掺比在万分之一及以上的污染。

表1：应用场景为600SNPs时的SCS值

实施例2：真实样本实验验证

在实验过程中，将两份来自于不同供体的cfDNA血样以不同梯度互掺，梯度掺比分布为万分之一，万分之五，千分之一，千分之五。其中样本一为掺入样本，样本二为待掺入样本，将样本一按照不同掺比梯度(0、0.01％、0.05％、0.1％、0.5％)掺入样本二中。本实施例中检测的可用SNP位点为1000个。从表2中结果可知，当梯度掺比为万分之五(0.05％)时，可以看出SCS数值与本底样本(即掺比梯度为0)相比，有较为显著的变化，依据模拟数据所定的SCS值阈值cutoff＝0.001，从实验结果验证了本方法学在检测1000个SNPs(参见表3)的条件下，可以判别万分之五的样本交叉污染。

根据实验cfDNA不同掺比情况，挑选本底(无掺比)样本，掺比为万分之一，万分之五，千分之一和千分之五的样本分布做AR频率分布图，如图3所示。本底样本绝大部分位点的AR值集中在两端，随着掺比的增加，AR的分布呈现向中间扩散的趋势。符合前述实施例2中在算法搭建过程中对待测样本遭到其它样本交叉污染后的AR值模型的评估。

表2：不同实验cfDNA掺比检测情况

表3：1000个SNP位点清单

其中：CHROM＝染色体编号；POS＝位置；REF＝参考序列碱基(即SNP位点未发生变异时野生型的碱基)；

ALT＝变异序列碱基(即SNP位点发生变异之后的变异型的碱基)

/>

/>

/>

/>

/>

/>

/>

/>

/>