一种基于新一代测序数据的Indel检测方法与流程

本发明属于基因工程技术领域，尤其涉及一种基于新一代测序数据的indel检测方法。

背景技术：

新一代测序是一种测序dna序列的技术。在测序过程中，将完整的样本dna序列打碎，从中筛选出满足特定长度(通常为数百bp)的片段，在每个片段的一端或两端各读取一段长度为数十至数百bp的序列。读取出的序列长度通常远远小于被测样本dna序列的长度，但是新一代测序技术可以同时读取大量这样的短序列，使得全部短序列的总长度达到样本dna长度的数倍至数十倍，从而使获得样本dna序列成为可能。indel(insertionanddeletion)变异是基因组中的一种重要的变异现象。主要表现为插入和删除两种状态，并且与人类的疾病发生相关。目前主要有4种检测基因组上indel变异的策略，分别为：(1)readpair(也称为pair-endmapping，简称pem，双端映射)；(2)splitread(简称sr)：分离读段。splitread是一类特殊的read，其出现通常是由基因组中的结构变异造成的。这类read在映射中不再保持连续序列的形式，而是包含了一定长度的空位，因此具有较高的映射难度；(3)readdepth(简称rd，读段覆盖深度)和(4)基于denovo组装的方法。(pem)将pair-endreads比对到参考序列上，若某一对reads插入长度小于映射长度，则这一对reads可以确定一个删除(deletion)；反之，若某一对reads插入长度大于映射长度，则可以确定一个插入(insertion)；对于序列删除的检测，其所能检测的片段长度受插入片段长度的标准差(sd)所影响(这里的插入片段长度指的是测序之前在构建dna测序文库阶段，所选取的经由超声波打断的dna片段长度，这些片段也称之为测序片段，这是实验过程中的操作，并不是指基因组的变异)，并且越大的序列删除越容易被检测到；对于序列插入的检测，长度只能在插入片段长度的范围内，并且最大长度也受限于测序的插入片段长度的标准差；这种检测方法的缺点是检测到的变异位置不够精确，不能达到bp级。sr首先提取具有以下特点的pair-endreads，一条正常比对到参考序列上，另外一条不能比对，然后利用正常比对的read位置和插入长度确定一个查找范围，在这个范围内寻找未比对上的read与参考序列的最佳匹配，通过最佳匹配点把未匹配的read分割成两段或者三段，从而确定deletion和insertion的位置；pindel是一个使用sr方法进行变异检测的软件。它在千人基因组计划和生物信息分析人员中被广泛使用。pindel理论上能够检测所有长度范围内的deletion，和小片段的insertion。pindel方法的一个优势在于它们能够精确到单碱基，但是在变异区域内若存在重复序列，pindel有可能会遗漏掉这些变异。rd通过samtools可以测得各位点的覆盖度，将测序reads比对到参考序列上，若某一段的覆盖深度低于平均覆盖深度很多，则可以确定这一段是一个deletion；缺点在于只能检测deletion，而不能检测insertion，并且检测位置也不够精确。denovoassembly的方法能够提供对于longinsertion的最好检测方法，但是组装仍然是一件棘手的事情，基因组上所存在的重复性序列会严重影响组装的质量，也在很大程度上阻碍了利用组装的方法在基因组变异检测方面的应用。

综上所述，现有技术存在的问题是：现有检测方法存在检测位置不精确；变异区域遇到重复序列容易遗漏；只依靠一条比对上的read和插入长度确定一个检测范围也容易造成变异检测的遗漏。

技术实现要素：

针对现有技术存在的问题，本发明提供了一种基于新一代测序数据的indel检测方法。

本发明是这样实现的，一种基于新一代测序数据的indel检测方法，所述基于新一代测序数据的indel检测方法包括以下步骤：

步骤一，利用bwa比对软件对原始的fastq数据做比对，生成sam文件；

步骤二，对discordant.sam文件中的每一对reads提取出比对位置当作一个二维点坐标，对二维点根据设定的阈值(此处阈值设定为插入长度)进行层次聚类，详细的聚类过程如下：

1、首先初始化每个二维点为一类，计算每两个类之间的距离，也就是每两个类之间的相似度，此处的距离为欧式距离，设有两点a[a1,a2]，b[b1,b2]，则a与b的欧式距离为

2、设定阈值为插入长度，找两个距离最近的类，若两个距离最近的类之间的距离小于此阈值，则将这两个类归为一类，转到3；若两个距离最近的类之间的距离大于等于此阈值，则终止迭代，转4；

3、重新计算新生成的这个类与各个旧类之间的相似度，转2；

4、结束；

步骤三，对hang.sam文件中的每一对reads，取出未正常比对的read，若其比对上的部分在某个聚类单元所代表的范围内，则将read插入聚类单元，经从hang.sam文件中提取出含有变异信息的read；

步骤四，每个聚类单元确定一个变异范围(对二维点数据进行层次聚类之后，每个聚类单元里有一些二维点，对这些二维点求均值，最终可以确定一个点[a,b]，从a到b即是变异发生的范围)，提取出此变异范围内含有的携带变异信息的reads，根据每一条read比对上的位置和变异的范围截取参考序列上的一段序列，将read和截取下来的参考序列做比对即可确定变异类型，变异位置，以及变异大小；

步骤五，将变异类型deletion记为1，insertion记为2，某一个确定的变异可以表示为“1_变异位置_变异大小”；然后利用现有技术哈希结构(哈希结构是计算机数据结构中的一种)来存储变异；对于某个变异，根据测序的覆盖度设置阈值，若支持某变异的reads数大于阈值，则可确定一个变异。

进一步，所述步骤一的sam文件中包括正常比对上pair-end数据、未正常比对的pair-end数据。

进一步，所述步骤二通过层次聚类并且设定阈值可以自动完成聚类，对每个聚类单元中的所有点求其平均值，聚类单元中含有的点如下：a[a1，b1]，b[a2，b2]，得到每个聚类单元含有一个范围[a,b]，a＝(a1+a2)/2，b＝(b1+b2)/2，即范围[a,b]内含有变异。

进一步，所述步骤四具体包括：read序列为序列a，截取的参考序列为序列b；从a和b的左端开始比较，遇到第一个碱基不相同的位置即为变异位置，记为q，然后从不相同的位置截取a序列，以a序列剩下的部分作为窗口开始滑动，起始位置为变异位置，每次向右滑动一个距离，窗口的得分函数为窗口内比对上的碱基的个数，若某个位置窗口的得分大于窗口内碱基的总数乘以0.95，则停止滑动，确定变异类型为deletion，记录此时的位置，记为w，w-q即为变异的大小，同理，在参考序列b上同样滑动，若存在某个位置匹配，即可确定变异类型为insertion。

本发明的优点及积极效果为：通过聚类确定一个变异的范围，提取splitread与变异范围内的参考序列进行比对，使得比对的过程变得简单，同时比对的范围更加精确，防止遗漏的发生；使用层次聚类，突破提前设置聚类个数的限制，并且使用层次聚类也会更符合此问题的需求；操作简单。

本发明只需要比对之后的sam文件和参考序列即可完成检测，同时，通过仿真数据的测试(仿真结果如表1所示)，本发明的测试结果相对于其它现有方法来说更为准确。本发明能够解决pem技术检测indel变异位置不精确的问题；解决sr方法检测indel变异造成遗漏的问题；解决现有技术遇到重复序列可能会错失变异点的问题。

表1：仿真数据实验结果

其中deletion为加入的5个变异，insertion为加入的1个变异，position为变异位置，size为变异大小，c_position为本方法检测到的变异位置，c_size为本方法检测到的变异大小，p_position为pindel方法检测到的变异位置，p_size为pindel方法检测到的变异大小。

附图说明

图1是本发明实施例提供的基于新一代测序数据的indel检测方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例提供的基于新一代测序数据的indel检测方法包括以下步骤：

s101：利用bwa比对软件对原始的fastq数据做比对，生成sam文件；

s102：对discordant.sam文件中的每一对reads提取出它们的比对位置当作一个二维点坐标，对二维点根据自己设定的阈值进行层次聚类；

s103：对hang.sam文件中的每一对reads，取出未正常比对的read，若其比对上的部分在某个聚类单元所代表的范围内，则将此read插入此聚类单元，经从hang.sam文件中提取出含有变异信息的read；

s104：每个聚类单元可以确定一个变异范围，提取出此变异范围内含有的携带变异信息的reads，根据每一条read比对上的位置和变异的范围截取参考序列上的一段序列，将read和截取下来的参考序列做比对即可确定变异类型，变异位置，以及变异大小；

s105：将变异类型deletion记为1，insertion记为2，某一个确定的变异可以表示为“1_变异位置_变异大小”；然后利用哈希结构来存储变异；对于某个变异，根据测序的覆盖度设置阈值，某变异reads支持数大于阈值则可确定这个变异。

本发明实施例提供的基于新一代测序数据的indel检测方法具体包括以下步骤：

(1)sam数据的预处理

利用bwa比对软件对原始的fastq数据做比对，生成sam文件，此sam文件中既含有正常比对上pair-end数据，又含有未正常比对的pair-end数据，提取出每一对正常比对的reads的比对位置，求差值即可得到映射长度，将映射长度不符合条件(即映射长度大于插入长度加方差或者映射长度小于插入长度减方差)的reads提取出来得到discordant.sam，并且将一条正常比对，另外一条没有正常比对上的pair-endreads也提取出来生成hang.sam。

(2)层次聚类

对(1)中discordant.sam文件中的每一对reads提取出它们的比对位置当作一个二维点坐标，对这些二维点根据自己设定的阈值(阈值默认设定为插入长度的平方)进行层次聚类，由于检测前并不知道样本里含有多少个变异，所以通过层次聚类并且自己设定阈值可以自动完成聚类，然后对每个聚类单元中的所有点求其平均值，设某聚类单元中含有的点如下：a[a1，b1]，b[a2，b2]，最终得到每个聚类单元含有一个范围[a,b]，a＝(a1+a2)/2，b＝(b1+b2)/2，即范围[a,b]内含有变异。

(3)将splitread插入到相应的聚类单元中

对(1)中hang.sam文件中的每一对reads，取出未正常比对的read，若其比对上的部分在某个聚类单元所代表的范围内，则将此read插入此聚类单元，经过这一步骤的处理，可以从hang.sam文件中提取出含有变异信息的reads。

(4)设置滑动窗口进行比对

(2)步骤每个聚类单元可以确定一个变异范围，(3)步骤可以提取出此变异范围内含有的携带变异信息的reads，根据每一条read比对上的位置和变异的范围截取参考序列上的一段序列，然后将此read和截取下来的参考序列做比对即可确定变异类型，变异位置，以及变异大小。具体方法为：设read序列为序列a，截取的参考序列为序列b，首先从a和b的左端开始比较，遇到第一个碱基不相同的位置即为变异位置，记为q，然后从不相同的位置截取a序列，以a序列剩下的部分作为窗口开始滑动，起始位置为变异位置，每次向右滑动一个距离，窗口的得分函数为窗口内比对上的碱基的个数，若某个位置窗口的得分大于窗口内碱基的总数乘以0.95，则停止滑动，可以确定变异类型为deletion，记录此时的位置，记为w，w-q即为变异的大小，同理，在参考序列b上同样滑动，若存在某个位置匹配，即可确定变异类型为insertion。

(5)筛选与确定

将变异类型deletion记为1，insertion记为2，这样某一个确定的变异可以表示为“1_变异位置_变异大小”，然后利用哈希结构来存储变异，哈希键为字符串，哈希值为整数类型(代表支持这个变异的read个数)，最终对于某个变异，根据测序的覆盖度设置阈值(默认设置为0.02)，即若支持此变异的read个数大于总reads数的百分之二，则将其输出。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。