1.本发明涉及生物信息学领域，具体涉及一种检测样本污染率的方法及装置。
背景技术：
2.dna甲基化是最早发现的基因表观修饰方式之一，可能存在于所有高等生物中，能够在不改变基因序列的前提下，改变遗传现象。它是基因调控的手段之一，即通过对位于启动子及第一外显子区的cpg岛的甲基化而抑制基因的表达，对生命活动非常重要。相比于普通的测序技术，甲基化测序会导致dna上面的碱基信息改变。主流的甲基化测序方法是使用重亚硫酸盐处理，会导致非甲基化的c碱基变换成t。而新的甲基化测序方法，tet酶和吡啶硼烷结合处理的方法(taps)会导致甲基化的c碱基变换成t。
3.目前，甲基化测序在肿瘤基因组中的应用越来越多，因此对于甲基化测序污染率的探索极其重要。甲基化测序样本污染主要存在三种情况，即个体内、不同个体间以及跨物种间的污染。跨物种间的污染比较容易解决，因为可以通过评估样本比对到物种参考基因组的情况，推测污染率的占比。但是同一物种不同个体间的污染则难以发现，因为污染可能发生在许多意想不到的场合，比如：样品存储过程中，收集样本的容器被污染；样品运输过程中，容器密封不严导致样本外溢；以及实验室人为操作制备过程中，不同样本移液时忘记更换枪尖或未使用带滤芯枪尖。
4.在现有的甲基化测序过程中，检测和计算样本间污染率对于甲基化测序的下游分析非常重要，即使只是少量的污染，也会导致分析结果出现许多假阳性或者假阴性，特别是在肿瘤与正常甲基化测序样本的比对研究中。因此，甲基化测序数据的样本间污染需要严格控制，但是遗憾的是，目前并没有相关软件或者流程能够实现对甲基化测序数据的污染率评估。如何实现对甲基化测序数据的污染率评估是目前亟待解决的问题。
技术实现要素：
5.根据第一方面，在一实施例中，提供一种检测样本污染率的方法，包括：
6.位点maf提取步骤，包括提取待测样本的测序数据中的snp位点在数据库中的最小等位基因频率(maf)；
7.过滤步骤，包括过滤去除不符合条件的snp位点；
8.错误率计算步骤，包括计算不同碱基替换的错误率；
9.似然值计算步骤，包括计算待测样本在不同污染率下的似然值；
10.候选污染率计算步骤，包括根据每个snp位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
11.优化步骤，包括根据优化函数优化所述候选污染率，获得最终的样本污染率。
12.根据第二方面，在一实施例中，提供一种检测样本污染率的装置，包括：
13.位点maf提取模块，用于提取待测样本的测序数据中的位点在数据库中的最小等位基因频率(maf)；
14.过滤模块，用于过滤去除不符合条件的位点；
15.错误率计算模块，用于计算不同碱基替换的错误率；
16.似然值计算模块，用于计算待测样本在不同污染率下的似然值；
17.候选污染率计算模块，用于根据每个snp位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
18.优化模块，用于根据优化函数优化所述候选污染率，获得最终的样本污染率。
19.根据第三方面，在一实施例中，提供一种检测样本污染率的装置，包括：
20.存储器，用于存储程序；
21.处理器，用于通过执行所述存储器存储的程序以实现第一方面的方法。
22.根据第四方面，在一实施例中，提供一种计算机可读存储介质，所述介质上存储有程序，所述程序能够被处理器执行以实现第一方面的方法。
23.依据上述实施例的一种检测样本污染率的方法及装置，该方法的分析结果可信度高。
24.在一实施例中，该方法在实际应用过程中，可通过总体评估甲基化测序样本在不同污染率下的评估污染率集合，在总体情况下，确定污染率阈值，能进一步增加下游分析结果的可靠性。
具体实施方式
25.下面通过具体实施方式对本发明作进一步详细说明。其中不同实施方式中类似元件采用了相关联的类似的元件标号。在以下的实施方式中，很多细节描述是为了使得本技术能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。在某些情况下，本技术相关的一些操作并没有在说明书中显示或者描述，这是为了避免本技术的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。
26.另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。
27.本文中为部件所编序号本身，例如“第一”、“第二”等，仅用于区分所描述的对象，不具有任何顺序或技术含义。
28.如本文所用，dbsnp是指ncbi于1998年建立的主要存储单核苷酸多态性(snp)的免费公共数据库。该数据库包含多种模式生物。虽然其名称为dbsnp，但该数据库实际上包括多种分子变异，具体如下：
29.单核苷酸多态性snp；
30.短缺失和插入多态性short deletion and insertion polymorphisms(indels/dips)；
31.微卫星标记或短串联重复microsatellite markers or short tandem repeats(strs)；
32.多核苷酸多态性multinucleotide polymorphisms(mnps)；
33.杂合序列heterozygous sequences；
34.命名变体named variants。
35.如本文所用，“maf”是指最小等位基因频率，通常是指在给定人群中的不常见的等位基因发生频率，例如tt、tc、cc三个基因型，在人群中c的频率＝0.36，t的频率＝0.64，则等位基因c就为最小等位基因频率，maf＝0.36。
36.如本文所用，“杂合子”是指同一位点上的两个等位基因不相同的基因型。
37.如本文所用，“纯合子”是指同一位点上的两个等位基因相同的基因型。
38.根据第一方面，在一实施例中，提供一种检测样本污染率的方法，包括：
39.位点maf提取步骤，包括提取待测样本的测序数据中的snp位点在数据库中的最小等位基因频率(maf)；
40.过滤步骤，包括过滤去除不符合条件的snp位点；
41.错误率计算步骤，包括计算不同碱基替换的错误率；
42.似然值计算步骤，包括计算待测样本在不同污染率下的似然值；
43.候选污染率计算步骤，包括根据每个snp位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
44.优化步骤，包括根据优化函数优化所述候选污染率，获得最终的样本污染率。
45.在一实施例中，所述位点maf提取步骤中，将snp位点映射到数据库的文件中，如果该位点存在，则保留该位点，并提取数据库中该位点的等位基因信息；如果该位点不存在，则删除该位点。
46.在一实施例中，所述位点maf提取步骤中，如果该位点存在，则提取其在数据库中的等位基因信息以及人群频率。
47.在一实施例中，所述位点maf提取步骤中，如果该位点存在，则提取其在数据库中的最小等位基因频率(maf)。
48.在一实施例中，所述位点maf提取步骤中，所述数据库包括但不限于dbsnp数据库、hapmap数据库中的至少一种。
49.在一实施例中，所述过滤步骤中，包括过滤并确定snp位点深度、基因型、位点先验污染概率和背景噪音读段(read)数。
50.在一实施例中，所述过滤步骤中，如果snp位点的等位基因不是t(胸腺嘧啶)、c(胞嘧啶)、g(鸟嘌呤)、a(腺嘌呤)或者n(未知碱基)，则过滤去除该位点。
51.在一实施例中，所述过滤步骤中，计算该位点主要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为主要占比；计算该位点的次要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为次要占比；根据主/次要占比判断样本基因型。
52.在一实施例中，所述过滤步骤中，如果主要占比小于杂合子位点的限制阈值，则将相应基因型判定为“1/1”；如果次要占比小于杂合子位点的限制阈值，则将相应基因型判定为“0/0”；其他则为“0/1”。
53.在一实施例中，所述限制阈值可以为0.25。
54.在一实施例中，所述过滤步骤中，如果基因型为“0/0”，背景噪音read计数为次要
sequencing)，该方法简称taps，无需亚硫酸氢盐，利用tet(ten-eleven translocation，简称tet)酶将5mc(5-甲基胞嘧啶)和5hmc(5-羟甲基胞嘧啶)氧化为5cac(5-羧甲基胞嘧啶)，随后使用有机硼烷(包括但不限于吡啶硼烷、2-甲基吡啶硼烷等等)将5cac还原为二氢尿嘧啶(dihydrouracil，dhu)，而后的pcr再将dhu转化为胸腺嘧啶(t)，对目标序列直接进行dna甲基化测序。本发明中，在进行污染率检测时，先通过中国专利《一种基于甲基化测序数据进行变异检测的方法及装置》(申请号：202110960293.8，公开号：cn113674802a)中的方法，将转化后的碱基修正为转化之前的碱基，具体是将t修正为c。
74.在一实施例中，甲基化测序方法不受限制，现有的亚硫酸氢盐测序法等其他甲基化测序方法获得的测序数据也适用于本发明。
75.根据第二方面，在一实施例中，提供一种检测样本污染率的装置，包括：
76.位点maf提取模块，用于提取待测样本的测序数据中的snp位点在dbsnp中的最小等位基因频率(maf)；
77.过滤模块，用于过滤去除不符合条件的snp位点；
78.错误率计算模块，用于计算不同碱基替换的错误率；
79.似然值计算模块，用于计算待测样本在不同污染率下的似然值；
80.候选污染率计算模块，用于根据每个snp位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
81.优化模块，用于根据优化函数优化所述候选污染率，获得最终的样本污染率。
82.根据第三方面，在一实施例中，提供一种检测样本污染率的装置，包括：
83.存储器，用于存储程序；
84.处理器，用于通过执行所述存储器存储的程序以实现第一方面任意一项的方法。
85.根据第四方面，在一实施例中，提供一种计算机可读存储介质，所述介质上存储有程序，所述程序能够被处理器执行以实现第一方面任意一项的方法。
86.在一实施例中，提供一种基于甲基化测序数据进行样本污染率检测的方法，包括：快速提取甲基化测序样本中的未转化位点的碱基信息；确定未转化位点在dbsnp中的最小等位基因频率；过滤并确定snp位点深度、基因型、位点污染先验概率和背景噪音read数；计算不同碱基替换的错误率；计算样本不同污染率下的似然值；选取最大似然值；根据优化函数优化最大似然对应的污染率。该方法可以应用于全基因组甲基化测序数据。
87.在一实施例中，提供一种基于甲基化测序数据进行样本污染率检测的方法，包括如下步骤：
88.(1)快速提取甲基化测序样本中的未转化位点的碱基信息。使用软件处理甲基化测序的压缩比bam文件，该软件能够检测甲基化测序样本的snp位点，并且能够将snp位点的碱基区分为经过转化或者未经过转化的碱基(具体方法参见申请号为202110960293.8的中国专利《一种基于甲基化测序数据进行变异检测的方法及装置》说明书第105～109段)，对于经过甲基化转化的碱基，修正为转化之前的碱基，即为未转化的碱基，通过提取未转化的碱基信息作为污染率检测软件的输入文件，评估样本污染率。
89.(2)确定未转化位点在dbsnp中的最小等位基因频率(maf)。将未转化的位点映射到dbsnp数据库的vcf文件中，如果该位点存在，则保留该位点并提取dbsnp中该位点的maf；如果该位点不存在，则删除该位点。
90.(3)过滤并确定snp位点深度、基因型、位点先验污染概率和背景噪音read数。如果snp位点的等位基因不是t、c、g、a或者n，则过滤去除该位点。
91.(4)计算不同碱基替换的错误率。每种碱基均有3种碱基替换形式，比如a碱基，可能的碱基替换形式有：a》t，a》g，a》c。例如要计算碱基替换a》t的错误率，仅需对样本中纯合子基因型进行分析，计算这部分基因型中主等位基因为a、次等位基因不为t的reads占参考碱基为a的reads的比值，即为a》t的错误率。
92.(5)计算样本在不同污染率下的似然值。根据给定的初始污染率集合，计算不同初始污染率下每个snp位点的对数似然值的比值，即用该位点的污染的似然值比去非污染的似然值，再取对数。
93.(6)选取最大似然值。将每个snp位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率。
94.(7)根据优化函数优化最大似然值对应的污染率。取候选污染率前后各一位污染率为优化函数的取值参数，再求全局最优解作为最终的样本对应的污染率。
95.实施例1
96.下述实施例中涉及的肿瘤样本由北京吉因加医学检验实验室有限公司提供。样本进行taps测序，下机数据预处理和基因组比对(使用bwa软件)得到压缩比对文件bam。bam文件为本方法的输入文件。
97.检测方法如下：
98.(1)读取bam文件，使用变异检测软件(具体方法参见申请号为202110960293.8的中国专利《一种基于甲基化测序数据进行变异检测的方法及装置》说明书第105～109段)处理，该软件能够检测甲基化测序样本的snp位点，并且能够将snp位点的碱基区分为经过转化或者未经过转化的碱基，通过提取未转化的碱基信息作为污染率软件的输入文件。该文件每一行为一个snp位点上未经过转化的碱基信息，具体包括：染色体、位点位置、位点深度、参考碱基、a、t、g、c、n碱基个数、ins个数、del个数。
99.(2)将上述文件snp位点映射到dbsnp数据库中。dbsnp数据库是单核苷酸多态性数据库，收录了snp、短插入缺失多态性和短重复序列等数据。将映射到的snp位点保留，并提取其在dbsnp数据库中的等位基因信息以及人群频率；在dbsnp数据库中未收录的位点，则删除该位点。
100.(3)过滤snp位点，如果snp位点的等位基因都不是单碱基或者n，则将这种类型的snp位点过滤掉。计算该位点主要等位基因的读段数(即reads个数)占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为主要占比；计算该位点的次要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为次要占比。根据主/次要占比判断样本基因型。如果主要占比小于杂合子位点的限制阈值(0.25)，则相应基因型判定为“1/1”；如果次要占比小于杂合子位点的限制阈值(0.25)，则相应基因型判定为“0/0”；其他则为“0/1”。如果基因型为“0/0”，背景噪音read计数为次要占比乘以位点深度，该位点的先验污染率为1-[1-(maf)]2；如果基因型为“1/1”，背景噪音read计数为主要占比乘以位点深度，该位点的先验污染率为1-(maf)2；如果基因型为“0/1”，背景噪音read计数为主要占比、次要占比中的最小值乘以位点深度，该位点的先验污染率为1-[2
×
(maf)
×
(1-maf)]。例如，等位基因a/t，a的计数是3，t的计数是
2，那么a为主要等位基因，b为次要等位基因。主要占比＝a count/(a count+b count)＝3/(3+2)＝0.6，次要占比＝0.4。因此与0.25比较，两者都大于0.25，因此基因型为0/1。
[0101]
(4)计算不同碱基替换的错误率。因为要计算碱基替换的错误率，所以仅需考虑基因型是纯合子的位点。假设计算a》t的错误率，错误率＝(样本中主等位基因是a、次等位基因不是t的read数和)
÷
(样本中参考碱基为a的全部read数
×
2/3)。这里之所以分母要乘以2/3，是因为计算a》t的错误率会和a》c或者a》g重复。其他的碱基替换方法相同。
[0102]
(5)给定样本初始污染率集合，具体为0.001、0.003、0.005、0.007、0.01、0.03、0.05、0.07、0.1、0.2、0.3。根据污染率集合中的每个值分别计算样本的污染率的概率值。首先移除杂合子snp位点；然后使用以下公式(1)分别计算每一个snp位点的污染似然值：
[0103][0104]
其中μ为位点先验污染率，ε为该位点的碱基替换错误率，α为样本污染率，n为位点深度，k为位点背景噪音read计数。
[0105]
再使用公式(2)计算每一个snp位点非污染似然值。最后将snp位点的污染似然值比去非污染似然值，再取对数(以自然数e为底)，即为snp位点的对数似然比值。
[0106]
公式(2)如下：
[0107][0108]
公式(2)中各符号的含义同公式(1)。
[0109]
(6)将上述所有snp位点的对数似然比值与位点深度计算加权平均值，则相应的11个初始污染率均会计算出对应的污染率概率值，选取其中最大的污染率概率值对应的污染率作为候选污染率。
[0110]
(7)挑选上述候选污染率。因为该污染率还不够精确，真实的污染率可能在候选污染率的前后区间内，因此将其前后各一位的初始污染率作为优化函数的取值范围，代入上述公式(1)，并根据公式(2)继续计算对数似然比值，计算全局最优解，即作为样本最终的污染率。
[0111]
本实施例对6例样本进行干实验模拟污染率情况，具体是将其他甲基化测序样本按一定比例混到当前样本中，其中每一例样本都分别模拟了0.001、0.005、0.01、0.05、0.1、0.2、0.3污染率的情况。按上述所示方法，对模拟样本进行污染率检测，具体检测结果如表1所示。
[0112]
表1六例模拟污染率样本的污染率检测结果
[0113]
污染率mix1mix2mix3mix4mix5mix6误差0.0010.00140.00150.00160.00140.00150.00140.047％0.0050.00510.00520.00520.00500.00530.00480.017％0.0100.01000.00990.00990.00950.01000.00920.029％0.0500.05650.05870.05280.05320.05510.05240.477％0.1000.12910.12890.12370.12430.13020.11892.585％0.2000.23100.22730.23110.22750.23390.22352.904％0.3000.28020.27600.28280.27780.28450.27582.048％
[0114]
通过上述结果可知，当模拟污染率为0.001时，6个样本计算的污染率与模拟污染
率的平均误差仅为0.047％；
[0115]
当模拟污染率为0.005时，6个样本计算的污染率与模拟污染率的平均误差仅为0.017％；
[0116]
当模拟污染率为0.01时，6个样本计算的污染率与模拟污染率的平均误差仅为0.029％；
[0117]
当模拟污染率为0.05时，6个样本计算的污染率与模拟污染率的平均误差仅为0.477％；
[0118]
当模拟污染率为0.1时，6个样本计算的污染率与模拟污染率的平均误差仅为2.585％；
[0119]
当模拟污染率为0.2时，6个样本计算的污染率与模拟污染率的平均误差仅为2.904％；
[0120]
当模拟污染率为0.3时，6个样本计算的污染率与模拟污染率的平均误差仅为2.048％。
[0121]
7种模拟污染率的平均误差仅为1.158％。显然，本方法的分析结果可信度高，并且在实际应用过程中，可通过总体评估甲基化测序样本在不同污染率下的评估污染率集合，在总体情况下，确定污染率阈值，能进一步增加下游分析结果的可靠性。
[0122]
在一实施例中，本发明可应用于甲基化测序样本中。
[0123]
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。
[0124]
以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。