基于双重测序数据检测低频突变的方法、装置及存储介质

1.本发明涉及生物信息学领域，特别涉及一种基于双重测序数据检测低频突变的方法和装置。


背景技术：

2.在肿瘤活检和循环无细胞核酸等dna样本中，突变可能以极低的频率(小于0.01％)存在于所测体细胞dna分子中。检测这些极低频率的体细胞突变在肿瘤早期诊断、监测和预后、法医鉴定、产前诊断等方面具有广阔的应用前景。新一代测序(next-generation sequencing,ngs)技术的发展改变了生物和医学科学领域的研究规模和深度。由于ngs具有大规模、高通量、低成本等特点，它不仅能实现大型基因组的分析，而且能有效地识别体细胞变异。然而，ngs的高错误率(错误率约为10-3-10-2
)掩盖了频率低于测序错误率的真实突变，使低频突变的检测仍然是具有如下挑战。第一，低频体细胞突变的检测需要深度测序。但是，增加测序深度进而增加了测序成本。第二，在测序模板量充足和测序深度足够的情况下，由于ngs工作流程中积累的伪影，使极低频率的突变仍然难以检测。这些伪影可能来源于样品制备过程中的dna碱基损伤、富集和文库扩增过程中dna聚合酶的错误碱基掺入以及最终测序读数的错误。为了提高低频变异的识别能力，科学家提出了一系列ngs纠错方法。基于分子标签(unique molecule identifier,umi)的双重测序能有效抑制高通量测序错误，是一种能够检测和量化极低频突变的方法。当文库制备时，该方法在原始dna模板的两端加上一段特有的标签序列，文库经pcr扩增和ngs测序后，进行测序数据分析。测序数据分析时，由正义链和负义链中相同的标签序列识别同一dna模板扩展出的多个读段(reads)分别组合聚集在一起，形成正义链和反义链的单链一致性序列(single-strand consensus sequence,sscs)；将生成的正义链sscs与互补的反义链sscs进行比较，进一步生成双链一致性序列(duplex consensus sequence，dcs)，将dcs与参考基因组再次进行比较，进行突变或者测序错误的识别。由于基于umi的双重测序方法利用正义链和反义链的配对原则进行进一步纠错，大幅提高了ngs测序错误的抑制效果。然而，由于仅保留含有读段大于等于3的read families用于生成sscs，因此造成sscs生成dcs的效率低，导致测序数据利用率低，而且与传统ngs相比，该方法需要更高的测序深度。


技术实现要素：

3.本发明旨在解决高通量测序错误率高导致低频突变检测困难的问题，提供一种基于双重测序数据检测低频突变的方法、装置及存储介质。该方法通过与参考基因组比对分别获取正义链和负义链的read families，不仅保留正义链和负义链数目大于等于2的read families，而且保留只含有1条读段的read family，且该读段能与上述读段数目大于等于2的read families生成的sscs序列互补的read family，以及保留均只含有一条读段，但所含读段之间互补的两组read families，以提高数据利用率。同时，通过采用贝叶斯定理计算每个位置上每种碱基是真实碱基的概率，选取概率最大的碱基作为一致性碱基并生成单
链一致性序列sscs，再根据此概率重新计算该碱基对应的质量，以进一步提高测序数据分析的准确性。该方法能有效移除提高测序原始数据的利用率，减少数据浪费以及提高测序准确率、降低测序深度，为低频突变的检测提供一种可靠的生物信息分析流程和工具，有望为医学检测和临床诊断提供一种有效的检测手段，以推动个体化医疗的发展。
4.鉴于此，本发明采用的技术方案是，基于双重测序数据检测低频突变的方法，包括以下步骤：
5.(1)对测序原始数据进行清洗处理，包括移除低质量读段和碱基、重复的读段和污染的接头序列，得到清洗后的数据。
6.(2)将清洗后的数据比对到参考基因组上，并根据比对结果分别建立正义链和反义链的read family，降低read family所包含的读段数目的阈值为1，并利用读段互补的特性筛选read family，构建有效read families。
7.(3)对构建的有效read families采用贝叶斯算法计算每个位置上每个碱基是真实碱基的概率，选取概率最大的碱基作为一致性碱基并生成单链一致性序列sscs，再根据此概率重新计算该碱基对应的质量。
8.(4)将步骤(3)生成的sscs进一步形成dcs。
9.(5)将所述dcs与参考基因组再次比对，根据比对结果，识别读段上的低频突变以及测序错误。
10.具体地，步骤(2)所述有效read families包括以下三种：1)读段数目大于等于2的read family；2)只含有1条读段的read family，且该读段能与所述读段数目大于等于2的read families生成的sscs序列互补；3)均只含有一条读段，但所含读段之间互补的两组read families。
11.具体地，步骤(3)所述贝叶斯算法首先精确计算每个位置上每种碱基是真实碱基的概率，选取概率最大的碱基作为一致性碱基并生成单链一致性序列sscs，再根据此概率重新计算该碱基对应的质量。
12.更优选，所述贝叶斯算法具体包括以下步骤：确定先验概率的方法为，如果比对到的碱基与可能的真实碱基一致，则先验概率为1-10-q/10
；如果比对到的碱基与可能的真实碱基不一致，则先验概率为10-q/10
/3，其中，q为碱基质量值，此分布用p(b,bi,qi)描述；对于4种可能的碱基a、g、c或者t，根据下述公式(1)计算后验概率：
[0013][0014]
对于sscs中每个位置的碱基，计算一致性碱基i为b时的概率(b∈{a,g,c,t}),概率值最大的碱基类型即为真实碱基；同时，根据公式(2)，使用上述计算所得的后验概率值重新计算碱基质量，得到纠错后的一致性读段；
[0015]
qc＝-10log
10
(1-p[i＝bc|{(bi,qi)}])
ꢀꢀ
公式(2)
[0016]
具体地，所述生成的单链一致性序列sscs包括：(ⅰ)如果read family的测序读段数目大于等于2，则保留该组read family，并采用贝叶斯算法精确计算每个位置上每种碱基是真实碱基的概率，选取概率最大的碱基作为一致性碱基并生成单链一致性序列sscs，再根据此概率重新计算该碱基对应的质量。
[0017]
(ⅱ)如果read family的测序读段数目等于1，且上述读段能与其他读段数目大于
等于2的read families生成的sscs序列互补，则保留该组read family，并用贝叶斯算法确定该条读段上每个位置的一致性碱基并生成单链一致性序列sscs，再根据概率重新计算碱基质量。(ⅲ)如果两组read family都只含有一条读段，但读段序列之间互补，则保留该组read family，并用贝叶斯算法确定该条读段上每个位置的一致性碱基并生成单链一致性序列sscs，再根据概率重新计算碱基质量。
[0018]
进一步，所述生成dcs包括：
[0019]
对于(ⅰ)中生成的正义链sscs与反义链sscs序列进行比较分析，如果相同位置上的碱基相匹配，则该位置上的碱基不变，并计算正义链和反义链上该位置的平均碱基质量作为该碱基的质量；如果相同位置上的碱基不匹配，则将该位置上的碱基用
‘
n’表示，并分配一个低的碱基质量(如10)，最终形成双链一致性序列dcs，但如果dcs中
‘
n’出现的比例大于50％，则去除该条dcs；如果生成sscs没有互补的链，但是由于其read family中含有两条或两条以上的读段，则该条sscs也被继续保留。对于(ⅱ)中生成dcs的形成方法如下：如果正义链sscs与反义链sscs序列相同位置上的碱基相匹配，则该位置上的碱基不变，并计算正义链和反义链上该位置的平均碱基质量作为该碱基的质量；如果相同位置上的碱基不匹配，则将该位置上的碱基用
‘
n’表示，并分配一个低的碱基质量(如10)，最终形成双链一致性序列dcs，但如果dcs中
‘
n’出现的比例大于50％，则去除该条dcs；对于(ⅲ)中将两组read family的sscs序列进行比较分析，如果正义链sscs与反义链sscs序列相同位置上的碱基相匹配，则该位置上的碱基不变，并计算正义链和反义链上该位置的平均碱基质量作为该碱基的质量；如果相同位置上的碱基不匹配，则将该位置上的碱基用
‘
n’表示，并分配一个低的碱基质量(如10)，最终形成双链一致性序列dcs，但如果dcs中
‘
n’出现的比例大于50％，则去除该条dcs。
[0020]
本发明还提供一种计算机可读存储介质，其存储有能实现上述基于双重测序数据检测低频突变的计算机程序。
[0021]
本发明的主要优点包括：
[0022]
1.本发明不仅保留读段数目大于等于2的read family，而且保留含有1条读段，且上述读段能与上述读段数目大于等于2的read families生成的sscs序列互补的read family，以及保留均只含有一条读段，但所含读段之间互补的两组read families。通过将read families包含的读段数目的阈值降低为1以及充分利用读段互补的特性来筛选read families，大幅地提高了原始测序数据的利用率。
[0023]
2.本发明贝叶斯算法重新计算每个位置上每种碱基是真实碱基的概率，选取概率最大的碱基作为该位置的真实碱基，再根据此概率重新计算该碱基对应的质量分数，该方法能非常有效地降低测序错误和文库制备过程很重引入的错误。
[0024]
3.本发明为微量模板以及低丰度基因突变提供一种灵敏度更高、准确性更好的高通量测序数据分析方法及工具。
附图说明
[0025]
图1为本发明基于双重测序数据检测低频突变的方法的原理图；(a)将原始测序数据比对到参考基因组；(b)根据比对结果生成read families，并对生成的read families进行过滤，保留读段数目大于等于2的有效read families，保留测序读段数目等于1，且上述
读段能与上述读段数目大于等于2的read families生成的sscs序列互补的read family，以及保留均只含有一条读段，但所含读段之间互补的两组read families；(c)根据贝叶斯算法使read families生成sscs；(d)sscs进一步生成dcs；(e)突变识别，获得双重测序数据中的突变和测序错误信息。
[0026]
图2为采用本发明的方法以及传统双重测序数据分析方法的错误率对比图。
具体实施方式
[0027]
为了能清晰地描述本发明的技术内容，下面结合实施例子来进一步描述。
[0028]
从美国国立生物技术信息中心ncbi(https://www.ncbi.nlm.nih.gov/)网站下载双重测序数据(accession number:srp140497)以及人类基因组参考序列hg19。本实施例中，对数据分别进行(a)本发明方法的基于双重测序数据检测低频突变的方法，也即保留3种有效read families，并采用贝叶斯算法确定每个位置上的一致性碱基及其碱基质量以生成sscs以及dcs和(b)传统双重测序数据的分析方法，也即仅保留读段read数目大于等于3的read families，并采用q30作为阈值进行质量控制。最后，对两种方法处理后的结果进行比较。具体方法为：
[0029]
a、对数据采用本发明方法进行分析
[0030]
(1)数据清洗
[0031]
对测序原始数据首先采用soapnuke或者fastp去除低质量碱基、污染的接头序列和低质量读段。此步骤包括，移除原始数据集中每个序列的前3bp(该数据集在文库制备的过程中添加了3bp，其中包含2bp的标签序列和单个固定碱基t)，编写python程序将标签序列添加到read id中，以便后续分类。
[0032]
(2)生成read families
[0033]
采用软件bwa(v 0.7.15)将测序数据与参考基因组序列hg19进行比对，并使用samtools(v 1.3)对序列进行排序，生成一个含有序列比对数据的有序bam文件。根据bwa的比对结果，去除未比对，以及比对上多个位置的读段。根据比对位置、分子标签序列、cigar标签和比对方向生成read families。而且，通过降低read families包含的读段数目的阈值为1以及有效利用读段互补的特性来生成有效read families。本发明中生成的有效read families包含以下3种：(1)如果read family的测序读段数目大于等于2；(2)如果read family的测序读段数目等于1，且该读段能与上述测序读段数目大于等于2的read families生成的sscs序列互补；(3)如果两组read families都只含有一条读段，但读段序列之间互补。
[0034]
(3)生成sscs
[0035]
在每个read family生成sscs的过程中，每个位置上碱基的确定，一般方法采用大多数规则，即计算此位置上每种碱基的比例，保留比例大于70％的碱基，认为此位置的真实碱基即为该碱基，同时使用其中较高的碱基质量值作为最终碱基质量值。此方法简单易行，但是对于测序仪器中存在的随机测序错误和pcr扩增引入的错误，仍然需要进一步提高碱基质量值的准确率。本发明为了进一步对测序数据进行校正，采用贝叶斯定理进一步确定一致性碱基，以便纠正测序的错误率。实施方式为：根据贝叶斯定理计算每个位置上每种碱基是真实碱基的概率，选取概率最大的作为一致性碱基，根据此概率计算一致性读段上每
个碱基的质量，使一致性读段的每个碱基更加准确可靠。由于此案例保留了3种read families生成sscs，因此分三种情况进行说明。
[0036]
第一种情况，如果read family的测序读段数目大于等于2，则保留该组read family，并采用贝叶斯算法精确计算每个位置是真实碱基的概率，且最大概率的碱基即指定为该位置的真实碱基，得到单链一致性序列sscs，然后再重新计算该碱基质量。第二种情况，如果read family的测序读段数目等于1，且该读段能与上述读段数目大于等于2的read families生成的sscs序列互补，则保留该组read family，并用贝叶斯算法确定该条读段上每个位置的一致性碱基并生成单链一致性序列sscs，再根据概率重新计算碱基质量。第三种情况，如果两组read families均只含有一条读段，但读段序列之间互补，则保留该组read family，并用贝叶斯算法确定该条读段上每个位置的一致性碱基并生成单链一致性序列sscs，再根据概率重新计算碱基质量。将上述三种情况生成的sscs写到一个bam1文件中。
[0037]
(4)生成dcs
[0038]
对于上述第(3)步测序读段数目大于等于2的read families，将上述所生成的正义链sscs与反义链sscs序列进行比较分析。如果相同位置上的碱基相匹配，则该位置上的碱基不变，并计算正义链和反义链上该位置的平均碱基质量作为该碱基的质量；如果相同位置上的碱基不匹配，则将该位置上的碱基用
‘
n’表示，并分配一个低的碱基质量(如10)，最终形成双链一致性序列dcs，但如果dcs中
‘
n’出现的比例大于50％，则去除该条dcs；如果生成sscs没有互补的链，但是由于其read family中含有两条或两条以上的读段，则该条sscs也被继续保留。对于上述第(3)步含有1条测序读段，且该读段能与上述读段数目大于等于2的read families生成的sscs序列互补的read family，寻找与上述所生成的sscs互补的sscs，如果序列相同位置上的碱基相匹配，则该位置上的碱基不变，并计算正义链和反义链上该位置的平均碱基质量作为该碱基的质量；如果相同位置上的碱基不匹配，则将该位置上的碱基用
‘
n’表示，并分配一个低的碱基质量(如10)，最终形成双链一致性序列dcs，但如果dcs中
‘
n’出现的比例大于50％，则去除该条dcs。对于上述第(3)步均只含有一条测序读段，但所含读段之间互补的两组read families，将上述两组read family的sscs序列进行比较分析。如果正义链sscs与反义链sscs序列相同位置上的碱基相匹配，则该位置上的碱基不变，并计算正义链和反义链上该位置的平均碱基质量作为该碱基的质量；如果相同位置上的碱基不匹配，则将该位置上的碱基用
‘
n’表示，并分配一个低的碱基质量(如10)，最终形成双链一致性序列dcs，但如果dcs中
‘
n’出现的比例大于50％，则去除该条dcs。将上述生成的dcs写到上述的bam1文件中。
[0039]
(5)突变识别
[0040]
采用软件bwa(v 0.7.15)将生成的dcs集合进一步与参考基因组进行比对，以获得双重测序数据中的突变和测序错误信息。
[0041]
b、对数据采用传统多重测序方法进行分析
[0042]
(1)数据清洗
[0043]
对测序原始数据首先采用soapnuke或者fastp去除低质量碱基、污染的接头序列和低质量读段。此步骤包括，移除原始数据集中每个序列的前3bp(该数据集在文库制备的过程中添加了3bp，其中包含2bp的标签序列和单个固定碱基t)，编写python程序将标签序
列添加到read id中，以便后续分类。
[0044]
(2)生成read families
[0045]
采用软件bwa(v0.7.15)将测序数据与参考基因组序列hg19进行比对，并使用samtools(v1.3)对序列进行排序，生成一个含有序列比对数据的有序bam文件。根据bwa的比对结果，去除未比对，以及比对多个位置的读段。根据比对位置、分子标签序列、cigar标签和比对方向生成read families。由于本方法中read families包含的读段数目的阈值为3。因此，本方法中生成的有效read families为测序读段数目大于等于3的read families。
[0046]
(3)生成sscs
[0047]
在每个read family生成sscs时，每个位置上碱基的确定，采用大多数规则，即计算此位置上每种碱基的比例，保留比例大于70％的碱基，认为此位置的真实碱基即为该碱基，同时使用其中较高的碱基质量值作为最终碱基质量值。对任意一个read family中的所有读段进行一致性处理，生成一致性读段。方法如下：每个read family中至少要有3条及3条以上的读段，且对每个read family中的每条读段上的碱基进行一致性比较，大于等于70％的一致位点为真，否则用n代替碱基信息。去除n的数目大于等于30％的读段。最终分别形成正义链和反义链的sscs。
[0048]
(4)生成dcs
[0049]
将正义链和反义链生成的sscs序列进一步进行比对分析，形成dcs。
[0050]
(5)突变识别
[0051]
采用软件bwa(v 0.7.15)将生成的dcs集合进一步与参考基因组进行比对，以获得双重测序数据中突变和测序错误的信息。
[0052]
对上述两种分析方法的错误率进行了比较，结果如图2所示，采用本发明方法分析上述测序数据后获得测序错误的百分比为0.000497183％；同时，采用传统多重测序方法分析上述测序数据后获得测序错误的百分比为0.000512317％。这表明本发明能够实现高质量的测序错误抑制。