CN104794371A - 检测逆转座子插入多态性的方法和装置 - Google Patents

Abstract

本发明公开一种检测逆转座子插入多态性的方法，包括：获取目标个体基因组测序结果；将测序结果与参考序列比对，获得异常匹配集，异常匹配集包括第一类读段对，第一类读段对中的每对读段中的两个读段中的一个至少能够与基因组参考序列匹配，另一个至少能够与TE参考序列匹配；按照匹配位置将异常匹配集中的第一类读段中的能够匹配到TE参考序列的读段聚类成簇；对聚类得到的簇进行处理，其中包括，过滤掉包含的读段的数目不大于1的簇；基于获得的处理后的簇，检测逆转座子插入多态性。本发明还提供一种检测逆转座子插入多态性的装置。本发明的方法和/或装置，能够快速、简便和准确的鉴定TE插入或者发现新的TE插入。

Description

检测逆转座子插入多态性的方法和装置

技术领域

本发明涉及生物信息领域，具体的，涉及检测逆转座子多态性的方法和装置。

背景技术

转座子(TE，Transposable elements)又被称为“跳跃基因”，它们通过复制-粘贴的增殖机制引起基因序列的插入，删除和重排。随着人类基因组测序的完成，人们发现有接近一半的基因组是由TE组成的。而超过90％的TE是逆转录转座子(retrotransposons)，逆转录转座子又分为包含长末端重复序列的LTR(Long Terminal Repeats)和不包含长末端重复的non-LTR，绝大多数的人类的TE主要来自于non-LTR逆转录转座子的活性，non-LTR包括LINE-1(L1；long interspersed element 1)，Alu和SVA，它们共同构成了人类基因组的三分之一。正是因为TE的转座特性，它可以引起个体基因组的结构变异，而并不是所有TE都具有活性，从能否编码蛋白质来看，只有non-LTR逆转录转座子具有活性，比如只有一小部分L1有高活性，它们通过编码特定的酶来进行转座，但是有些TE在生殖细胞和早期胚胎的生长过程中，这种抑制机制的短暂释放间隙会逃脱这种抑制并产生新的多态性的插入，从而出现了逆转录转座子的插入多态性(RIPs，retrotransposon insertion polymorphisms)。

有很多关于RIPs的研究，已经有一些方法可以检测TE插入，目前，用于检测可移动元件插入的方法主要分为两类，一是目标法，比如Transposon-Seq，ME-Scan，RC-Seq等，需要测序前对与TE相关的DNA片段做PCR实验丰富它的序列信息；二是后测序生物信息学方法，比如Variation Hunter，RetroSeq等，用全基因组测序数据来鉴定TE插入的多态性。

随着高通量测序技术的发展，能够利用高通量测序数据快速、简便、准确的鉴定TE插入或者发现新的TE插入的方法仍亟需开发。

发明内容

本发明旨在至少解决上述问题之一或者提出一种商业选择手段。

依据本发明的一方面，本发明提供一种检测逆转座子插入多态性的方法，包括以下步骤：获取目标个体基因组测序结果，所述测序结果包括多对读段对，每对读段对由两个读段组成，分别来自一条染色体片段的两端，每对读段对分别来自所述染色体片段的正链和负链，或者，每对读段对同时来自所述染色体片段的正链或所述染色体片段的负链；将所述测序结果与参考序列进行比对，获得异常匹配集，所述异常匹配集包括第一类读段对，所述第一类读段对中的每对读段中的两个读段中的一个至少能够与基因组参考序列匹配，另一个至少能够与TE参考序列匹配，所述TE参考序列包括以下四种类型中的至少一种：Alu参考序列、L1参考序列、SVA参考序列和LTR参考序列；按照匹配位置将所述异常匹配集中的第一类读段中的能够匹配到所述TE参考序列的读段聚类成簇(block)，所述簇的大小不大于一个TE插入的大小；对聚类得到的簇进行处理，其中包括，过滤掉包含的读段的数目不大于1的簇，以及任选的，将相邻的、最大距离不大于D_b且方向相反的簇合并为一个簇，其中，D_b为两个TE的大小；基于获得的处理后的簇，检测所述逆转座子插入多态性。

所说的染色体片段通常是将来自目标个体的基因组核酸打断获得的，根据所选用的测序方法进行相应的文库(library)制备，可选用的测序方法根据来自的测序平台包括但不限于CG(Complete Genomics)、Illumina/Solexa、Life Technologies/Ion Torrent和Roche 454，依据所选测序平台进行单端或双端测序文库的制备。在本发明一个实施例中进行双末端(Pair-end read)测序，获得多对读段对，每对读段对中的两个读段(reads)，可表示为reads1和reads2，reads1和reads2可能都来自相应染色体片段的正链或负链，也可能分别来自染色体片段的正链和负链。当然，若使用的单端(single-read)测序方法能够完整获得整个染色体片段的序列，从完整获得的序列的两端分别截取适当长度的序列来构成一对reads、或者将获得的序列截成两部分序列来构成一对reads也是可行的。本实施例对所选用的具体测序方法不作限定。所说的正链和负链是相对的，称一条双链序列的一条单链为正链，就可称另一条单链为负链，在本发明的一个实施例中，将染色体片段的两条链中的与基因组参考序列相同的那条链称为正链。

比对可以利用SOAP(Short Oligonucleotide Analysis Package)，BWA，Samtools等软件进行，本实施对此不作限制。所说的参考序列是预先确定的序列，可以是预先获得的目标个体所属生物类别中的任意的参考模板，例如，同一生物类别的已公开的基因组组装序列，若目标个体是人类，其基因组参考序列(也称为参考基因组)可选择NCBI数据库提供的HG19。进一步地，也可以预先配置包含更多参考序列的资源库，在进行序列比对前，先依据目标个体的性别、人种、地域等因素选择或是测定组装出更接近的序列来作为参考序列，有助于获得更准确的检测结果。所说的TE参考序列为包含已知的TE特异序列的序列，TE特异性序列包括Alu、L1、SVA和LTR至少之一。在比对过程中，根据比对参数的设置，例如设置测序结果中的每条或每对读段最多允许有n个碱基错配(mismatch)，n优选为1或2，若reads中有超过n个碱基发生错配或者比对质量值小于预设值，则视为该条/对reads无法比对上参考序列。一般利用比对软件进行比对后，都可获得诸如是否为唯一比对即是否为只比对到参考序列的一个位置、比对上参考序列多个位置的各个比对位置的比对质量值等评估比对情况的参数。在本发明的一个实施例中，所述比对利用BWA软件进行，设置每条reads允许的最多错配数为2且比对质量值不小于10，获得异常匹配集包括：将所述测序结果与所述基因组参考序列比对，获得初级异常匹配集，所述初级异常匹配集包括符合以下(i)-(iii)至少之一的读段：(i)匹配到所述基因组参考序列的多个位置，(ii)匹配到所述基因组参考序列的唯一位置，并且比对质量值小于10，(iii)匹配到所述基因组参考序列的唯一位置，比对质量值大于10，并且匹配到所述基因组参考序列的至少一个次优比对位置，所称的次优比对为不满足比对上的所有条件但至少满足其中之一的比对情况；将所述初级异常匹配集比对到所述TE参考序列，获得包含比对上TE参考序列的读段的比对结果，所述比对结果构成所述异常匹配集。所说的“匹配”或“匹配到”同“比对上”。这样，将可能支持TE插入的reads都筛选出来。

聚类可采用各种聚类算法，本发明对此不作限定。例如，一种简单的做法是，按照设置的簇的大小不大于1个插入的TE的大小，一般一个插入TE的大小约为1000bp，可设置一个簇的大小为不大于1000bp，依据匹配位置对reads进行排序，第一条reads的第一个碱基与某一条reads的最后一个碱基之间的距离不大于1000bp，就可将所说的第一条reads和该某一条reads以及其间的所有reads聚类成一个簇，聚类时可不考虑reads的方向，这里的“第一”、“最后”是相对而言的，指在参考序列上相距最远的。任选的，若相邻的两个簇其最大距离不大于D_b，且方向是相反的，说明这两个簇包含同一个TE插入，可将这两个相邻的簇合并为一个簇。这里，所说的两个簇之间的最大距离为不考虑簇的方向，为这两个簇中的在参考序列上距离最远的两个碱基之间的距离。结合簇的大小，可知相邻两个簇之间的距离是非负值，即最小为0，即相连。簇的方向为其所包含的方向一致的多数派读段的方向，例如一个簇中包含的大于一半的reads与参考序列同向，这该簇为正向。

在本发明的一个实施例中，所述TE参考序列包括以下四种类型中的至少两种：Alu参考序列、L1参考序列、SVA参考序列和LTR参考序列。所说的聚类得的簇不大于一个TE插入大小，基于获得的处理后的簇，检测所述逆转座子插入多态性包括，基于处理后的簇中的最多读段匹配的TE参考序列类型，确定所述逆转座子插入的类型，以及任选的，基于处理后的簇中支持该类型逆转座子插入的读段中的方向一致的多数派读段的方向是否与所述TE参考序列一致，判断所述逆转座子插入的方向。例如，处理后的簇中匹配到Alu参考序列的读段最多，则判定所发生的TE插入为Alu插入。进一步可选的，当匹配到Alu参考序列的读段中的方向一致的多数派读段与所述Alu参考序列的方向相反，则判定该Alu插入突变是逆向插入的。读段的方向是相对的，若一读段与参考序列比对上，即与参考序列的方向相同，则定义为正向；若一读段的反向互补链与参考序列比对上，即与参考序列的方向相反，则定义为反向(逆向)。

在本发明的一个实施例中，检测所述逆转座子插入多态性还包括确定所述逆转座子插入的位置，包括获取候选断点集和筛选所述候选断点集，所述候选断点集包含多个候选断点，获取所述候选断点集包括：将所述处理后的簇中的读段与所述基因组参考序列比对，将获得的比对结果中的包含至少一个长度不小于10bp的割裂部分的一级读段中的割裂位置加入到所述候选断点集，定义所述一级读段的大于10bp的割裂部分为二级读段，将所述二级读段比对到所述基因组参考序列，将获得的比对结果中的包含至少一个长度不小于10bp的割裂部分的二级读段中的割裂位置加入到所述候选断点集，定义所述二级读段的大于10bp的割裂部分为三级读段，如此，直至N级读段的任一割裂部分小于10bp，其中，N为自然数，N≥2。所说的候选断点为可能的TE插入的位置，所称的断点(breakpoint)指TE插入的边界点。其中，所说的一级读段为处理后的簇中的割裂读段(split-reads或者soft-clipped reads)，割裂读段指一条reads被切成至少两段，匹配到基因组参考序列的不同区域的读段。这里，将割裂读段的每段都称为割裂部分，所定义的二级读段实为前述割裂读段的一部分，所定义的N级读段实为(N-1)级读段的一部分。如此，能够获得割裂读段支持的所有可能的插入位置边界点，以确定所述逆转座子插入的确切位置。

在本发明的一个实施例中，获取所述候选断点集还包括，将i级读段与所述基因组参考序列的比对结果中的满足任一以下(a)-(c)的错配位置加入到所述候选断点集，(a)所述错配位置包含多聚序列(poly)，所述多聚序列包含的碱基满足：6个碱基中至少有5个是一样的碱基，且所述一样的碱基为A或T，(b)在所述比对结果中支持所述错配的读段的比例小于1/2，(c)增加和/或删除所述错配位置中的预定碱基数后，将所述i级读段与所述基因组参考序列重新比对，获得的比对结果显示所述i级读段为割裂的，所述预定碱基数为所述比对结果中的插入和/或缺失的碱基的数目；其中，所述i级读段为所述一级读段、所述二级读段至所述N级读段中的至少之一。依据(a)、(b)和/或(c)，可以将在所述比对结果中显示为错配但实际可能是断点的给包含进候选断点集中。其中，(a)是根据TE特异序列polyA/T尾来设定的；(b)是利用匹配上的读段的支持比例来判定该错配是不是真实的来设定的，所说的支持所述错配的读段为该读段的错配位置上的碱基与参考序列的对应位置的碱基不同，如参考序列上的碱基为A，比对上该位点的读段的该位置上的碱基为非A，则该读段为支持所述错配的读段，所说的支持所述错配的读段的比例为支持所述错配的读段的数目与比对上该位置的读段的总数的比值；(c)是将比对结果中显示错配的位置包含的插入/缺失消除后，再次检测判断该错配位置为一个可能的断点的给包含进所述候选断点集。

在本发明的一个实施例中，筛选所述候选断点集包括：过滤掉处理后的簇中匹配读段的数目不大于2的候选断点，从所述候选断点的位置开始所述匹配读段的长度不小于10bp，所述匹配读段为匹配到所述候选断点的读段，当任一所述匹配读段都不包含所述多聚序列时，过滤掉处理后的簇中支持读段的比例不大于0.5的候选断点，所述支持读段的比例＝支持所述候选断点的读段数目/匹配到所述断点的读段数目。经此筛选得的候选断点集中的候选断点为真实断点。

在本发明的一个实施例中，所说方法还包括对所述TE插入所带的靶位点重复序列的类型进行判定，将满足以下的逆转座子插入判定为带双端靶位点重复序列，两相邻j级读段在所述基因组参考序列上的距离小于30bp，所述j级读段为所述所述二级读段至所述N级读段中的至少之一。

在本发明的一个实施例中，所述方法还包括，通过以下公式计算所述逆转座子插入的杂合概率 $P\left(\mathrm{Heter}\right|N_C,N_S)=\frac{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)}{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)+B(N_S,N_C+N_S,P_{\mathrm{HO}})F\left(\mathrm{HO}\right)},$ 当MSL大于1则判定所述逆转座子插入为杂合，当MSL小于1则判定所述逆转座子插入为纯合，以判断所述逆转座子插入的类型，其中，MSL＝-log₂(1-P(Heter|N_C，N_S))，NC为比对上所述逆转座子插入位置中的割裂读段的数目，NS比对上所述逆转座子插入位置中的跨越读段的数目，P_HE＝0.6和P_HO＝0.1分别是期望杂合概率和期望纯合概率，F(HE)＝F(HO)＝0.5，F(HE)和F(HO)分别杂合和纯合的权重值，B(N_S，N_C+N_S，P_HE)是二项分布的概率，即为概率质量函数，表示指定事件发生的概率是P_HE、在N_C+N_S次独立重复试验中该指定事件发生P_HE次的概率，其计算公式可参考概率论等教科书，在概率论中，概率质量函数(probability massfunction，简写为pmf)是离散随机变量在各特定取值上的概率。所称的割裂读段，也称为截断读段，指一条reads能够分成匹配到(比对上)不同基因组位置的至少两部分。这里，所称的跨越读段(跨越reads)指一对reads中的至少一个能够比对上基因组序列上，另一个能够比对上TE参考序列。

本发明的方法能够对RIPs进行鉴定与挖掘，能够基于核酸序列测序数据中的跨越读段快速、简便、准确的鉴定TE插入或者发现新的TE插入，而且本领域普通技术人员能够知道，开发程序或者组合程序子模块、借助机器运行程序能够快速实现本发明或者各个实施例中的方法。例如，通过程序指令相关硬件来完成上述实施例中各种方法的全部或部分步骤，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

依据本发明的另一方面，本发明提供一种检测逆转座子插入多态性的装置，该装置能够用以实施上述本发明任一实施例中的检测逆转座子插入多态性的方法，该装置包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与所述数据输入单元、数据输出单元及存储单元数据连接，用于执行所述可执行的程序，所述程序的执行包括完成上述任一实施例中的检测逆转座子插入多态性的方法的全部或部分步骤。前述对本发明一方面的RIPs的检测方法的技术特征和优点的描述，同样适用本发明这一方面的检测装置，在此不再赘述。所称的处理器一般是负责运算和处理，存储单元中的至少一部分为内存，主要为交换数据。当程序或者操作者对处理器发出指令，这些指令和输入数据暂存在内存里，在处理器空闲时传送给处理器，处理器处理后把结果输出到诸如显示器、打印机等输出设备上。在输出完成之前，这些结果也保存在内存里，如果内存不足，结果数据读取的速度会变慢很多。执行该程序能快速实现本发明一方面的方法的全部或部分步骤，对内存的需求比目前的方法低。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是本发明的一个实施例中的检测TE插入多态性的方法的示意图。

图2是本发明的一个实施例中的检测TE插入多态性的方法的示意图。

图3是本发明的一个实施例中的测序深度对TE检测结果的影响的示意图。

具体实施方式

图1显示本发明的一个实施方式中的检测TE插入多态性的方法的步骤流程，该方法包括以下步骤：S10获取目标个体基因组测序结果，测序结果包括多对读段对，每对读段对由两个读段组成，分别来自一条染色体片段的两端，每对读段对分别来自一条染色体片段的正链和负链，或者，每对读段对同时来自一条染色体片段的正链或该染色体片段的负链；S20将测序结果与参考序列进行比对，获得异常匹配集，异常匹配集包括第一类读段对，第一类读段对中的每对读段对中的两个读段中的一个至少能够与基因组参考序列匹配，另一个至少能够与TE参考序列匹配，TE参考序列包括以下四种类型中的至少一种：Alu参考序列、L1参考序列、SVA参考序列和LTR参考序列；S30按照匹配位置将异常匹配集中的第一类读段中的能够匹配到TE参考序列的读段聚类成簇(block)，所述簇的大小不大于一个TE插入的大小；S40对聚类得到的簇进行处理，其中包括，过滤掉包含的读段的数目不大于1的簇，以及任选的，将相邻的、最大距离不大于D_b且方向相反的簇合并为一个簇，其中，D_b为两个TE的大小；S50基于获得的处理后的簇，检测所说的逆转座子插入多态性。所说的染色体片段通常是将来自目标个体的基因组核酸打断获得的，根据所选用的测序方法进行相应的文库(library)制备，可选用的测序方法根据来自的测序平台包括但不限于CG(Complete Genomics)、Illumina/Solexa、Life Technologies/Ion Torrent和Roche454，依据所选测序平台进行单端或双端测序文库的制备。在本发明一个实施例中进行双末端(Pair-end read)测序，获得多对读段对，每对读段对中的两个读段(reads)，可分别表示为reads1和reads2，reads1和reads2可能都来自相应染色体片段的正链或负链，也可能分别来自染色体片段的正链和负链。当然，若使用的单端(single-read)测序方法能够完整获得整个染色体片段的序列，从完整获得的序列的两端分别截取适当长度的序列来构成一对reads、或者将获得的序列截成两部分序列来构成一对reads也是可行的。本实施例对所选用的具体测序方法不作限定。所说的正链和负链是相对的，称一条双链序列的一条单链为正链，就可称另一条单链为负链，在本发明的一个实施例中，将染色体片段的两条链中的与基因组参考序列相同的那条链称为正链。所说的比对可以利用SOAP(ShortOligonucleotide Analysis Package)，BWA，Samtools等软件进行，本实施对此不作限制。所说的参考序列是预先确定的序列，可以是预先获得的目标个体所属生物类别中的任意的参考模板，例如，同一生物类别的已公开的基因组组装序列，若目标个体是人类，其基因组参考序列(也称为参考基因组)可选择NCBI数据库提供的HG19。进一步地，也可以预先配置包含更多参考序列的资源库，在进行序列比对前，先依据目标个体的性别、人种、地域等因素选择或是测定组装出更接近的序列来作为参考序列，有助于获得更准确的检测结果。所说的TE参考序列包括已知的TE特异序列的序列，TE特异性序列包括Alu、L1、SVA和LTR中的至少之一。在比对过程中，根据比对参数的设置，例如设置测序结果中的每条或每对读段最多允许有n个碱基错配(mismatch)，n优选为1或2，若reads中有超过n个碱基发生错配或者比对质量值小于预设值，则视为该条/对reads无法比对上参考序列。一般利用比对软件进行比对后，都可获得诸如是否为唯一比对即是否为只比对到参考序列的一个位置、比对上参考序列多个位置的各个比对位置的比对质量值等评估比对情况的参数。在本发明的一个实施例中，所述比对利用BWA软件进行，设置每条reads允许的最多错配数为2且比对质量值不小于10，获得异常匹配集包括：将所述测序结果与所述基因组参考序列比对，获得初级异常匹配集，所述初级异常匹配集包括符合以下(i)-(iii)至少之一的读段：(i)匹配到基因组参考序列的多个位置，(ii)匹配到基因组参考序列的唯一位置，并且比对质量值小于10，(iii)匹配到基因组参考序列的唯一位置，比对质量值大于10，并且匹配到基因组参考序列的至少一个次优比对位置，所称的次优比对为不满足比对上的所有条件但至少满足其中之一的比对情况；将初级异常匹配集比对到所述TE参考序列，获得包含比对上TE参考序列的读段的比对结果，所称比对结果构成所说的异常匹配集。所说的“匹配”或“匹配到”同“比对上”。这样，将可能支持TE插入的reads都筛选出来。

聚类可采用各种聚类算法，本发明对此不作限定。例如，一种简单的做法是，按照设置的簇的大小不大于1个插入的TE的大小，一般一个插入TE的大小约为1000bp，可设置一个簇的大小为不大于1000bp，依据匹配位置对reads进行排序，第一条reads的第一个碱基与某一条reads的最后一个碱基之间的距离不大于1000bp，就可将所说的第一条reads和该某一条reads以及其间的所有reads聚类成一个簇，聚类时可不考虑reads的方向，这里的“第一”、“最后”是相对而言的，指在参考序列上相距最远的。任选的，若相邻的两个簇其最大距离不大于D_b，且方向是相反的，说明这两个簇包含同一个TE插入，可将这两个相邻的簇合并为一个簇。这里，所说的两个簇之间的最大距离为不考虑簇的方向，为这两个簇中的在参考序列上距离最远的两个碱基之间的距离。结合簇的大小，可知相邻两个簇之间的距离是非负值，即最小为0，即相连。簇的方向为其所包含的方向一致的多数派读段的方向，例如一个簇中包含的大于一半的reads与参考序列同向，这该簇为正向。

在本发明的一些实施例中，TE参考序列包括以下四种类型中的至少两种：Alu参考序列、L1参考序列、SVA参考序列和LTR参考序列。所说的聚类得到的簇不大于一个TE插入大小，基于获得的处理后的簇，检测RIPs包括：基于处理后的簇中的最多读段匹配的TE参考序列类型，确定逆转座子插入的类型。和/或，基于处理后的簇中支持该类型逆转座子插入的读段中的方向一致的多数派读段的方向是否与TE参考序列一致，判断该逆转座子插入的方向。例如，处理后的簇中匹配到Alu参考序列的读段最多，则判定所发生的TE插入为Alu插入。进一步可选的，当匹配到Alu参考序列的读段中的方向一致的多数派读段与所述Alu参考序列的方向相反，则判定该Alu插入突变是逆向插入的。读段的方向是相对的，若一读段与参考序列比对上，即与参考序列的方向相同，则定义为正向；若一读段的反向互补链与参考序列比对上，即与参考序列的方向相反，则定义为反向(逆向)。

如图2所示，在本发明的一个实施例中，检测逆转座子插入多态性还包括：S60确定逆转座子插入的位置，即确定TE插入的发生位置，发生位置也称为断点。S60包括获取候选断点集和筛选所述候选断点集，候选断点集包含多个候选断点，所称候选断点集的获取包括：将处理后的簇中的读段与基因组参考序列比对，将获得的比对结果中的包含至少一个长度不小于10bp的割裂部分的一级读段中的割裂位置加入到所述候选断点集，定义所说的一级读段的大于10bp的割裂部分为二级读段，将二级读段比对到基因组参考序列，将获得的比对结果中的包含至少一个长度不小于10bp的割裂部分的二级读段中的割裂位置加入到候选断点集，定义二级读段的大于10bp的割裂部分为三级读段，如此，直至N级读段的任一割裂部分小于10bp，其中，N为自然数，N≥2。所说的候选断点为可能的TE插入的位置，所称的断点(breakpoint)指TE插入的边界点。其中，所说的一级读段为处理后的簇中的割裂读段(split-reads或者soft-clipped reads)，割裂读段指比对时一条reads被切成至少两段、且各段能够匹配到基因组参考序列的不同区域的读段。这里，将割裂读段的每段都称为割裂部分，所定义的二级读段实为前述割裂读段的一部分(一段)，所定义的N级读段实为(N-1)级读段的一部分。如此，能够获得割裂读段支持的所有可能的插入位置边界点，用以确定逆转座子插入的确切位置。

在本发明的一个实施例中，获取候选断点集还包括，将i级读段与基因组参考序列的比对结果中的满足任一以下(a)-(c)的错配位置加入到候选断点集中，(a)错配位置包含多聚序列(poly)，多聚序列包含的碱基满足：6个碱基中至少有5个是一样的碱基，且一样的碱基为A或T，(b)在比对结果中支持所述错配的读段的比例小于1/2，(c)增加和/或删除所述错配位置中的预定碱基数后，将该i级读段与基因组参考序列重新比对，获得的比对结果显示该i级读段为割裂的，所说的预定碱基数为比对结果中的插入和/或缺失的碱基的数目；其中，所述i级读段为所述一级读段、所述二级读段至所述N级读段中的至少之一。依据(a)、(b)和/或(c)，可以将在所述比对结果中显示为错配但实际可能是断点的给包含进候选断点集中。TE的存在也可以通过明显的TE的结构特点来鉴定，比如TSD(target site duplications)和poly(A)延伸。其中，(a)是根据TE特异序列polyA/T尾来设定的；(b)是利用匹配上的读段的支持比例来判定该错配是不是真实的来设定的，所说的支持所述错配的读段为该读段的错配位置上的碱基与参考序列的对应位置的碱基不同，如参考序列上的碱基为A，比对上该位点的读段的该位置上的碱基为非A，则该读段为支持所述错配的读段，所说的支持所述错配的读段的比例为支持所述错配的读段的数目与比对上该位置的读段的总数的比值；(c)是将比对结果中显示错配的位置包含的插入/缺失消除后，再次检测判断该错配位置为一个可能的断点的给包含进所述候选断点集。

在本发明的一个实施例中，筛选所述候选断点集包括：过滤掉处理后的簇中匹配读段的数目不大于2的候选断点，从所述候选断点的位置开始所述匹配读段的长度不小于10bp，所述匹配读段为匹配到所述候选断点的读段，当任一所述匹配读段都不包含所述多聚序列时，过滤掉处理后的簇中支持读段的比例不大于0.5的候选断点，所述支持读段的比例＝支持所述候选断点的读段数目/匹配到所述断点的读段数目。经此筛选得的候选断点集中的候选断点可认定为真实断点。

在本发明的一个实施例中，所说方法还包括对所述TE插入所带的靶位点重复序列的类型进行判定，将满足以下的逆转座子插入判定为带双端靶位点重复序列，两相邻j级读段在基因组参考序列上的距离小于30bp，所述j级读段为所述所述二级读段至所述N级读段中的至少之一。

高深度测序数据促进了TE的检测，因为高深度数据可以检测到更多TSD周围的截断reads，就能够提高检测的敏感度，但是因为测序和映射误差使得紧邻这些位置的测序深度明显低于平均测序深度，这就不可避免的增加了检测难度和假阳性率。在传统检测程序中，截断reads是唯一可以确定断点位置的信息，但实际上截断信息并不是一直可靠的，比如Indel(即缺失)，如果在比对时将Indel区域错判则会导致假的断点出现，出现假阳性结果，而断点附近的Indel，则容易引起假阴性结果，比如比对时错误的将原本截断的reads判断为存在indel的情况。截断点附近的错配也会引起假阴性的判断结果，因为组装软件通常都把这些碱基当作是错配而不是截断，所以一些包含错配碱基的reads末端也可能是截断部分。在实际比对时，截断点位置会和实际情况有偏差，导致断点错位等问题。另外截断位置reads支持数过低或过高也会在实际检测中会出现影响准确性的结果，为了检测出某些断点支持度很低的TE插入，对于断点支持数的标准较低，由此会带来一些假阳性；而支持度过高，则可能是重复区域，存在映射错误的风险。在实际检测中还会出现一种影响判断结果的情况，当存在poly T或poly A序列的时候，截断部分的一致性有时会很差，附近的错配碱基有时也会变得很多，给检测带来难度，如ACE基因区域。以上各个实施例中的检测方法中的设置或者处理，能够避免或者消除部分或全部以上一些影响，提高检测准确度。

在本发明的一个实施例中，所述方法还包括，通过以下公式计算所述逆转座子插入的杂合概率 $P\left(\mathrm{Heter}\right|N_C,N_S)=\frac{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)}{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)+B(N_S,N_C+N_S,P_{\mathrm{HO}})F\left(\mathrm{HO}\right)},$ 当MSL大于1则判定所述逆转座子插入为杂合，当MSL小于1则判定所述逆转座子插入为纯合，以判断所述逆转座子插入的类型，其中，MSL＝-log₂(1-P(Heter|N_C，N_S))，NC为比对上所述逆转座子插入位置中的割裂读段的数目，NS比对上所述逆转座子插入位置中的跨越读段的数目，P_HE＝0.6和P_HO＝0.1分别是期望杂合概率和期望纯合概率，F(HE)＝F(HO)＝0.5，F(HE)和F(HO)分别杂合和纯合的权重值，B(N_S，N_C+N_S，P_HE)是二项分布的概率，即为概率质量函数，表示指定事件发生的概率是P_HE、在N_C+N_S次独立重复试验中该指定事件发生P_HE次的概率，这里的指定事件为TE插入杂合的发生概率，B(N_S，N_C+N_S，P_HE)的计算公式可参考概率论等教科书，在概率论中，概率质量函数(probability mass function，简写为pmf)是离散随机变量在各特定取值上的概率。所称的割裂读段，也称为截断读段，指一条reads能够分成匹配到(比对上)不同参考序列位置的至少两部分。这里，所称的跨越读段(跨越reads)指一对reads中的至少一个能够比对上基因组序列上，另一个能够比对上TE参考序列。

本发明的方法能够对RIPs进行鉴定与挖掘，能够基于核酸序列测序数据中的跨越读段快速、简便、准确的鉴定TE插入或者发现新的TE插入，而且本领域普通技术人员能够知道，开发程序或者组合程序子模块、借助机器运行程序能够快速实现本发明或者各个实施例中的方法。例如，通过程序指令相关硬件来完成上述实施例中各种方法的全部或部分步骤，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

根据本发明的另一实施方式，本发明提供一种检测逆转座子插入多态性的装置，该装置能够用以实施上述本发明任一实施例中的检测逆转座子插入多态性的方法，该装置包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与所述数据输入单元、数据输出单元及存储单元数据连接，用于执行所述可执行的程序，所述程序的执行包括完成上述任一实施例中的检测逆转座子插入多态性的方法的全部或部分步骤。前述对本发明一方面的RIPs的检测方法的技术特征和优点的描述，同样适用本发明这一方面的检测装置，在此不再赘述。所称的处理器一般是负责运算和处理，存储单元中的至少一部分为内存，主要为交换数据。当程序或者操作者对处理器发出指令，这些指令和输入数据暂存在内存里，在处理器空闲时传送给处理器，处理器处理后把结果输出到诸如显示器、打印机等输出设备上。在输出完成之前，这些结果也保存在内存里，如果内存不足，结果数据读取的速度会变慢很多。实行本发明的方法时可以将整个分析聚类过程采用异步扫描法，即在第一步检测聚类block后，对于每个block，立即检测候选断点/断点，如果一个block经断点检测后不支持TE插入，则删除这个block及所有这个block区域内的reads存储，以避免重复扫描及减小信息存储量。这样，执行该方法程序能快速实现本发明一方面的方法的全部或部分步骤，对内存的需求比目前的方法低。

下面详细描述本发明的实施例，所述实施例中的示例图，其自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。需要说明的是在本文中所使用的术语“第一”、“第二”等仅用于方便描述目的，而不能理解为指示或暗示相对重要性，也不能理解为之间有先后顺序关系。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。

除另有交待，以下实施例中涉及的试剂、仪器或者软件，都是常规市售产品或者开源的，比如从Illumina公司购买测序文库制备试剂盒，依照试剂盒说明书进行建库等。

一般方法包括的内容或模块：

1、筛选可能支持插入信息的reads

1)筛选不和谐的reads

用BWA软件将深度基因组测序结果比对到人类参考基因组上，得到bam文件，从中提取出其中不和谐的reads作为疑似支持TE插入的reads，具体提取reads方法如下：根据BWA比对结果中的标签标记(tag)进行筛选：筛选出tag为XT:A:R的reads；如果标签为XT:A:U，筛选MAPQ小于10的；如果MAPQ大于10，筛选X1：i标签值大于0的。各个标签及其数值含义可参见BWA软件使用说明，例如其中的XT:A这个标签用于标记reads是唯一比对还是多处比对，它有4个值U、R、M、N，U表示uniq比对(唯一比对),R表示repeat比对(重复比对)，M表示这个比对是根据Mate-SW算法得出的，N表示XN:i里记录的N的个数大于10。MAPQ(MAPping Quality)为比对质量。X1标签为BWA找到的次优比对位置的个数。

2)将疑似reads通过与TE库中的参考序列进行blast比对，blast设置期望值为-e 0.2，-F参数设置为F，即不屏蔽简单重复和低复杂度序列，提取可以比对到TE库中的reads作为潜在存在TE的reads。

2、聚类分析

1)将潜在reads进行聚类，聚类单位称之为block，聚类过程用异步扫描法来实现，即对所有潜在reads进行扫描，将reads聚类为单边block，并以比对到基因组上的第一个碱基坐标作为block的标签位置进行排序，将支持reads数为1的block删除，即弃去只包含1条reads的block。

2)确定block的类型，始末端位置及TE类型的注释

根据上一步得到的单边block的方向性和相邻block的位置，以及相邻block之间的距离，判断block是否双端block。如果相邻block方向相反，且相邻block距离小于阈值2000，则判断相邻block为一个双端block，进一步确定这个双端block的5’端和3’端。如果不是双端block，单端block的reads支持数大于等于阈值5，则认为是一个单端block，再根据该block的方向，判断其另一端的位置，给定block长度的阈值为1000bp。

在与TE库进行blast比对结果中，通过四个亚类Alu、L1、SVA和LTR各自的一致性序列即特异序列比对，注释出每个潜在reads所支持的TE类型。统计block内所有支持的reads被注释的TE类型，如果TE插入的类型能匹配到两个或两个以上的亚类，那就定类型为匹配到最多的亚类。

通过reads的方向判断TE插入的方向，如果某个TE插入与reads的方向不同，就判定TE插入的方向是逆向的，插入的方向是由更多与它方向相同的reads来确定的。

3、检测可能的断点位置，即获得候选断点集合

1)进行断点位置修正，对于某一个block，从所有比对上的reads(mapping reads)中找到在该block中的每条reads，从bam文件的CIGAR字符中提取截断长度，根据截断长度初步确定断点，设定TSD长度的阈值为10bp，与参考基因组序列进行比较，根据reads的方向在断点10bp范围内搜索是否存在错配，挑出错配位置后重新确定断点及截断长度。CIAGR字符为读段碱基序列中的符号^标记，可用于提取该读段中的的子序列(subsequence)。

2)考虑TSD附近的indel，从bam文件中确定该条reads是否存在indel，从错配数当中删掉Indel的长度。判断该断点是否是一个poly或者截断reads，如果是poly或者错配数低于阈值5，并且是一个截断reads，认定该位置是一个断点，统计支持该断点的reads数。

3)对TSD附近的错配进行筛选，对于block区域内的reads，如果mapping结果中它被认定为非截断，而是存在错配，检测它是否为是一个截断reads，重新确定断点，统计支持该断点的reads数。

我们通过截断reads来确定断点，这种reads的比对情况可表示为mSnM或者nMmS，m和n分别表示截断长度和比对长度，大写字母表示序列，S表示TSD开始的位置，从比对结果来看，m开始的位置就是TSD开始的位置也就是所说的断点位置，但是为了更精确TSD位置，我们进一步将reads中n部分截取出来与相应位置上的基因组去逐个碱基比较，确定实际的m开始的位置，那么这里所说的“TSD附近”就指的是一条reads中的m部分，也就是说它最大不会超过一条reads的长度。

4、筛选候选断点

对上个步骤所有检测到的候选断点进行筛选，设定断点reads支持数阈值2，即一个候选断点至少得有2条reads比对上才可能得以保留，TSD的长度阈值10，过滤这些不能支持TE插入的block。检测支持的reads的一致性，为不同reads支持数水平设置不同的一致性判断标准，最大程度的降低由低reads支持度带来的假阳性结果，过滤掉支持reads非一致性的断点，同时检测是否有在参考序列中不存在的特异性ploy A/T序列，对于这些特异性ploy A/T序列，跳过一致性检验。删除两头都是截断的reads。poly的检测标准如下：如果支持reads在尾端的六个碱基中有五个或以上的连续碱基A，就注释插入具有poly-A尾，相反，如果在头端的六个碱基中有五个或五个以上的连续碱基T，就注释插入具有poly-T。在这些情况下，TE插入既不是poly-A也不是poly-T：(1)插入既有支持reads是poly-A，也有支持reads是poly-T，(2)插入的方向与poly-A或者poly-T不一致，比如，插入的方向是正的，却具有poly-T。

5、TSD检测

首先检测TSD的类型，设定阈值为截断reads的碱基数大于10，即截断长度大于10bp，若TSD距离小于30，判定为双端TSD。然后检测跨TSD区域的比对质量，质控条件包括平均错配数小于3，reads支持数大于或等于3，上下游平均qual(quality)分别大于25，和大于65以及TSD附近是否有indel。保留质控合格的断点输出其TSD区域的信息。对于不满足双端TSD条件的断点，再检测其是否为合格的单端TSD，设定单侧reads支持数阈值5，碱基数阈值5。对于单端TSD，根据错配数和reads支持数判断截断方向，与另一侧25个碱基范围内检测跨TSD区域的比对质量，同样保留质控合格的断点输出其TSD区域的信息。在bam文件中qual列表示序列质量(query QUALity)，可以认为是reads本身的quality，那么这里的“上下游平均qual分别大于25”指的是控制支持断点的reads本身的quality，这个条件值可以根据测序方法等原因去自行设定，此为一建议值。

进一步的，还可以计算插入位点的杂合度，根据插入位点截断reads数N_C和跨越reads数N_S构建贝叶斯模型预测插入位点的杂合度。P_HE＝0.6和P_HO＝0.1分别是期望杂合概率和期望纯和概率，F(HE)＝F(HO)＝0.5是权重值，B(k,n,P)是二项分布的概率函数，即概率质量函数，其计算公式可参考《概率论》。P(Heter|N_C，N_S)表示TE插入的杂合概率， $P\left(\mathrm{Heter}\right|N_C,N_S)=\frac{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)}{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)+B(N_S,N_C+N_S,P_{\mathrm{HO}})F\left(\mathrm{HO}\right)},$ 进而，通过MSL值判断TE插入的杂合性，MSL大于1则为杂合，MSL小于1则为纯合，MSL＝-log₂(1-P(Heter|N_C，N_S))。

6、单边再搜索

过滤两边都是poly和TSD为重复A/T序列的插入位点，过滤至少有一端深度过高，例如，大于150X的序列，保留单边TSD进行再搜索，过滤深度低于60X，杂合度小于等于10，单边clip reads计算得的深度/单边深度>1/4且无ploy的TSD。

优选的，上述整个分析聚类过程采用异步扫描法进行检测，即在第一步检测聚类block后，对于每个block，立即检测断点，如果一个block经断点检测后不支持TE插入，则删除这个block及所有这个block区域内的reads存储，以避免重复扫描及减小信息存储量。

通过上述过程检测RIPs，检测准确性高，可靠程度高，检测速度快以及适用于高深度测序数据。

上述过程的通过一系列对断点附近的mismatch、indel及reads支持数的特异性分析，能够极大的降低了复杂的RIPs检测过程中的假阳性和假阴性率。以下实施例通过对不同数据集中的RIPs进行检测，并且与其他软件检测结果的比较，来示意本发明的方法和/装置的实施，并且证明本发明方法的检测准确性和敏感度都高。

实施例一

对在Broad Institute基于Illumina HiSeq平台的高深度(>75X)测序数据CEU trio data进行TE插入检测。这个CEU trio数据集包含一个父亲NA12891，一个母亲NA12892和他们的女儿NA12878，他们都是千人基因组计划中的样本，测序数据可从http://hapmap.ncbi.nlm.nih.gov/citinghapmap.html获得。

1、筛选可能支持插入信息的reads

1)筛选不和谐的reads

用BWA将三个个体的测序结果比对到人类参考基因组hg19上，得到bam文件，从中提取出其中不和谐的reads作为疑似支持TE插入的reads，具体提取reads方法如下：根据BWA比对结果产生的标签进行筛选，首先筛选XT:A:R的reads，如果标签为XT:A:U，筛选MAPQ小于10的，如果MAPQ大于10，筛选X1：i标签值大于0的。

2)将疑似reads通过与TE库中的参考序列进行blast比对，期望值设置为-e 0.2，-F参数设置为F，即不屏蔽简单重复和低复杂度序列，提取可以比对到TE库中的reads作为潜在存在TE的reads。

2、聚类分析

1)将潜在的可能支持TE插入的reads进行聚类，聚类单位称之为block，聚类过程用异步扫描法来实现，对所有reads进行扫描，将reads聚类为单边block，并以比对到基因组上的第一个碱基坐标作为block的标签位置进行排序，支持reads数为1的block被删除。

2)确定block的类型，始末端位置、TE类型的注释和/或TE插入的方向

根据上一步得到的单边block的方向性和相邻block的位置，以及相邻block之间的距离，判断block是否为双端block，如果相邻block方向相反，相邻block距离小于阈值2000，则判断为双端block，确定block的5’端和3’端。如果不是双端block，单端block的reads支持数大于等于阈值5，则认为是一个真的单端block，再根据该block的方向，判断其另一端的位置，给定block长度的阈值为1000bp。

在与TE库进行blast比对结果中，已经通过四个亚类Alu、L1、SVA和LTR各自的一致性序列比对，注释出每个潜在reads所支持的TE类型。统计block内所有支持reads被注释的TE类型，如果TE插入的类型能匹配到两个或两个以上的亚类，那就是被匹配到最多的亚类。

通过reads的方向判断TE插入的方向，如果某个TE插入与reads的方向不同，就判定TE插入的方向是逆向的，插入的方向是由更多与它方向相同的reads判定的

3、检测可能的断点的位置

首先进行断点位置修正，对于某一个block，从所有mapping reads中找到在block中的每条reads，从bam文件的CIGAR字符中提取截断长度，根据截断长度初步确定断点，给定TSD长度的阈值为10bp，与参考基因组序列进行比较，根据read的方向在断点10bp范围内搜索是否存在错配，挑出错配位置后重新确定断点及截断长度。

然后考虑TSD附近的indel，从bam文件中确定该条read是否存在indel，从错配数当中删掉Indel的长度。判断该断点是否是一个poly或者截断reads，如果是poly或者错配数低于阈值5，并且是一个截断reads，认定该位置是一个断点，统计支持该断点的reads数。

再对TSD附近的错配进行筛选，对于block区域内的reads，如果mapping结果中它被认定为非截断，而是存在错配，检测它是否为一个截断reads，若是截断reads则依据上述过程重新确定断点，统计支持该断点的reads数。

4、筛选出断点

对所有检测到的断点进行筛选，设定断点reads支持数阈值2，及TSD的长度阈值10，过滤这些不能支持断点插入的block。检测支持reads的一致性，为不同reads支持数水平设置不同的一致性判断标准，最大程度的降低由低reads支持度带来的假阳性结果，过滤掉支持reads非一致性的断点，同时检测是否有在参考序列中不存在的特异性ploy A/T序列，对于这些特异性ploy A/T序列，跳过一致性检验。删除两头都是截断的reads。poly的检测标准如下：如果支持reads在尾端的六个碱基中有五个或五个以上的连续碱基A，就注释插入具有poly-A尾，相反，如果支持reads在头端的六个碱基中有五个或五个以上的连续碱基T，就注释插入具有poly-T。在这些情况下，TE插入既不是poly-A也不是poly-T：1)插入既有支持reads是poly-A，也有支持reads是poly-T；2)插入的方向与poly-A或者poly-T不一致，比如，插入的方向是正的，却具有poly-T。

5.TSD检测

首先检测TSD的类型，设定阈值为截断reads的碱基数大于10，TSD距离小于30，判定为双端TSD。然后检测跨TSD区域的map质量，质控条件包括平均错配数小于3，reads支持数大于等于3，上下游平均qual分别大于25，和大于65以及TSD附近是否有indel。保留质控合格的断点输出其TSD区域的信息。对于不满足双端TSD条件的断点，再检测其是否为合格的单端TSD，设定单侧reads支持数阈值5，碱基数阈值5。对于单端TSD，根据错配数和reads支持数判断截断方向，与另一侧25个碱基范围内检测跨TSD区域的map质量，同样保留质控合格的断点输出其TSD区域的信息。

计算插入位点的杂合度，根据插入位点截断reads数N_C和跨越reads数N_S构建贝叶斯模型预测插入位点的杂合度。P_HE＝0.6和P_HO＝0.1分别是期望杂合概率和期望纯和概率，F(HE)＝F(HO)＝0.5是权重值，B(k,n,P)是二项分布的概率。

$P\left(\mathrm{Heter}\right|N_C,N_S)=\frac{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)}{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)+B(N_S,N_C+N_S,P_{\mathrm{HO}})F\left(\mathrm{HO}\right)}$

最后，通过MSL值判断TE插入的杂合性，MSL大于1则为杂合，MSL小于1则为纯合。

MSL＝-log₂(1-P(Heter|N_C，N_S))

6、单边再搜索

过滤两边都是poly和TSD为重复A/T序列的插入位点，过滤至少有一端深度过高(>150)的序列，保留单边TSD进行再搜索，过滤深度低于60、杂合度小于等于10、利用单边clipreads计算出的深度/单边深度>1/4且无ploy的TSD。

整个分析聚类过程采用异步扫描法进行检测，即在第一步检测聚类block后，对于每个block，立即检测断点，如果一个block经断点检测后不支持TE插入，则删除这个block及所有这个block区域内的reads存储，以避免重复扫描及减小信息存储量。

7、检测结果

统计结果如表1所示。表1显示CEU trio数据集利用本发明方法检测RIPs的结果及敏感度。将以上步骤方法开发成程序，命名程序为Sid。

表1

注：敏感度是根据对经过PCR验证的数据集GP1000(Stewart,C.,et al.,A comprehensive map of mobileelement insertion polymorphisms in humans.PLoS Genet,2011.7(8):p.e1002236.)的覆盖度计算。

表1显示了用Sid软件检测到的RIPs可以很高的覆盖到已经经过PCR验证的TE插入，尤其是Alu类型，覆盖度高达97％。

实施例二

同上个实施例，NA18571，NA18572和NA18537也为千人计划中个体，各测序数据仍旧可以通过千人基因组计划数据库获得，经过重测序后平均测序深度为69X，用BWA比对人类参考基因组hg19上，又用GATK标准流程进行合并，本地重排，标记重复和比对校正后得到bam格式的测序数据。

1、筛选可能支持插入信息的reads

1)筛选不和谐的reads

提取出其中不和谐的reads作为疑似支持TE插入的reads，具体提取reads方法如下：根据BWA比对结果产生的标签进行筛选，首先筛选XT:A:R的reads，如果标签为XT:A:U，筛选MAPQ小于10的，如果MAPQ大于10，筛选X1：i标签值大于0的。

2)将疑似reads通过与TE库中的参考序列进行blast比对，期望值设置为-e 0.2，-F参数设置为F，即不屏蔽简单重复和低复杂度序列，提取可以比对到TE库中的reads作为潜在存在TE的reads。

2、聚类分析

1)将潜在reads进行聚类，聚类单位称之为block，聚类过程用异步扫描法来实现，对所有reads进行扫描，将reads聚类为单边block，并以比对到基因组上的第一个碱基坐标作为block的标签位置进行排序，支持reads数为1的block被删除。

2)确定block的类型，始末端位置及TE类型的注释

根据上一步得到的单边block的方向性和相邻标签位置(标签为指示区分block)，以及相邻标签之间的距离，判断block是否双端block，如果相邻block方向相反，相邻标签距离小于阈值2000，则判断为双端block，确定block的5’端和3’端。如果不是双端block，单端block的reads支持数大于或等于阈值5，则认为是一个单端block，再根据该block的方向，判断其另一端的位置，给定block长度的阈值为1000bp。

在与TE库进行blast比对结果中，已经通过四个亚类Alu、L1、SVA和LTR各自的一致性序列比对，注释出每个潜在reads所支持的TE类型。统计block内所有支持reads被注释的TE类型，如果TE插入的类型能匹配到两个或两个以上的亚类，那就是被匹配到最多的亚类。

通过reads的方向判断TE插入的方向，如果某个TE插入与reads的方向不同，就判定TE插入的方向是逆向的，插入的方向是由更多与它方向相同的reads判定的。

3、检测断点位置

首先进行断点位置修正，对于某一个block，从所有mapping reads中找到在block中的每条read，从bam文件的CIGAR字符中提取截断长度，根据截断长度初步确定断点，给定TSD长度的阈值为10bp，与参考基因组序列进行比较，根据read的方向在断点10bp范围内搜索是否存在错配，挑出错配位置后重新确定断点及截断长度。

然后考虑TSD附近的indel，从bam文件中确定该条read是否存在indel，从错配数当中删掉Indel的长度。判断该断点是否是一个poly或者截断reads，如果是poly或者错配数低于阈值5，并且是一个截断reads，认定该位置是一个端点，统计支持该断点的reads数。

再对TSD附近的错配进行筛选，对于block区域内的reads，如果mapping结果中它被认定为非截断，而是存在错配，检测它是否为是一个截断read，重新确定断点，统计支持该断点的reads数。

4、筛选断点

对所有检测到的断点进行筛选，设定断点reads支持数阈值2，及TSD的长度阈值10，过滤这些不能支持断点插入的block。检测支持reads的一致性，为不同reads支持数水平设置不同的一致性判断标准，最大程度的降低由低reads支持度带来的假阳性结果，过滤掉支持reads非一致性的断点，同时检测是否有在参考序列中不存在的特异性ploy A/T序列，对于这些特异性ploy A/T序列，跳过一致性检验。删除两头都是截断的reads。poly的检测标准如下：如果支持reads在尾端的六个碱基中有五个或以上的连续碱基A，就注释插入具有poly-A尾，相反，如果在头端的六个碱基中有五个或五个以上的连续碱基T，就注释插入具有poly-T。在这些情况下，TE插入既不是poly-A也不是poly-T：1)插入既有支持reads是poly-A，也有支持reads是poly-T；2)插入的方向与poly-A或者poly-T不一致，比如，插入的方向是正的，却具有poly-T。

5、TSD检测

首先检测TSD的类型，设定阈值为截断reads的碱基数大于10，TSD距离小于30，判定为双端TSD。然后检测跨TSD区域的map质量，质控条件包括平均错配数小于3，reads支持数大于等于3，上下游平均qual分别大于25，和大于65以及TSD附近是否有indel。保留质控合格的断点输出其TSD区域的信息。对于不满足双端TSD条件的断点，再检测其是否为合格的单端TSD，设定单侧reads支持数阈值5，碱基数阈值5。对于单端TSD，根据错配数和reads支持数判断截断方向，与另一侧25个碱基范围内检测跨TSD区域的map质量，同样保留质控合格的断点输出其TSD区域的信息。

计算插入位点的杂合度，根据插入位点截断reads数N_C和跨越reads数N_S构建贝叶斯模型预测插入位点的杂合度。P_HE＝0.6和P_HO＝0.1分别是期望杂合概率和期望纯和概率，F(HE)＝F(HO)＝0.5是权重值，B(k,n,P)是二项分布的密度函数。

$P\left(\mathrm{Heter}\right|N_C,N_S)=\frac{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)}{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)+B(N_S,N_C+N_S,P_{\mathrm{HO}})F\left(\mathrm{HO}\right)}$

最后，通过MSL值判断TE插入的杂合性，MSL大于1则为杂合，MSL小于1则为纯和。

MSL＝-log₂(1-P(Heter|N_C，N_S))

6、单边再搜索

过滤两边都是poly和TSD为重复A/T序列的插入位点，过滤至少有一端深度过高(>150)的序列，保留单边TSD进行再搜索，过滤深度低于60，杂合度小于等于10，单边clip深度/单边深度>1/4且无ploy的TSD。

整个分析聚类过程采用异步扫描法进行检测，即在第一步检测聚类block后，对于每个block，立即检测断点，如果一个block经断点检测后不支持TE插入，则删除这个block及所有这个block区域内的reads存储，以避免重复扫描及减小信息存储量。

7、检测结果

统计结果如表2所示。表2中显示Sid与TEA检测RIPs结果的比较及敏感度，Sid检测到的RIPs，30％与GP1000数据集覆盖度达到90％外其余70％为检测到的新的RIPs，与TEA软件检测结果相比，覆盖到83％以上。异步扫描法的应用使得本方法检测RIPs的过程非常高效，表3显示利用不同软件检测该测序数据中RIPs所需的运行时间，依据本发明方法所开发的sid软件的速度较另外两种软件提高了5到7倍。

表2

表3

实施例三

模拟数据的结果显示，当测序深度为30X时，可以检测到95％的TE插入，假阳性率低于0.6％，如图3，图3显示测序深度对TE检测结果的影响，从图中可看出随着测序深度在一定范围内增加，检测敏感度也随之上升，假阳性率逐渐降低，测序深度高于30X时，检测的高敏感度和低假阳性率表现稳健。

Claims (10)

1.一种检测逆转座子插入多态性的方法，其特征在于，包括以下步骤：

获取目标个体基因组测序结果，所述测序结果包括多对读段对，每对读段对由两个读段组成，分别来自一条染色体片段的两端，每对读段对分别来自所述染色体片段的正链和负链，或者，每对读段对同时来自所述染色体片段的正链或所述染色体片段的负链；

将所述测序结果与参考序列进行比对，获得异常匹配集，所述异常匹配集包括第一类读段对，所述第一类读段对中的每对读段中的两个读段中的一个至少能够与基因组参考序列匹配，另一个至少能够与TE参考序列匹配，所述TE参考序列包括以下四种类型中的至少一种：Alu参考序列、L1参考序列、SVA参考序列和LTR参考序列；

按照匹配位置将所述异常匹配集中的第一类读段中的能够匹配到所述TE参考序列的读段聚类成簇，所述簇的大小不大于一个TE插入的大小；

对聚类得到的簇进行处理，其中包括，

过滤掉包含的读段的数目不大于1的簇，以及任选的，

将相邻的、最大距离不大于D_b且方向相反的簇合并为一个簇，其中，

D_b为两个TE插入的大小；

基于获得的处理后的簇，检测所述逆转座子插入多态性。

2.权利要求1的方法，其特征在于，所述比对利用BWA软件进行；

任选的，

获得异常匹配集，包括，

将所述测序结果与所述基因组参考序列比对，获得初级异常匹配集，所述初级异常匹配集包括符合以下(i)-(iii)至少之一的读段：

(i)匹配到所述基因组参考序列的多个位置，

(ii)匹配到所述基因组参考序列的唯一位置，并且比对质量值小于10，

(iii)匹配到所述基因组参考序列的唯一位置，比对质量值大于10，并且能够匹配到所述基因组参考序列的至少一个次优比对位置，

将所述初级异常匹配集比对到所述TE参考序列，获得包含比对上TE参考序列的读段的比对结果，所述比对结果构成所述异常匹配集。

3.权利要求1的方法，其特征在于，所述TE参考序列包括以下四种类型中的至少两种：Alu参考序列、L1参考序列、SVA参考序列和LTR参考序列。

4.权利要求3的方法，其特征在于，基于获得的处理后的簇，检测所述逆转座子插入多态性，包括，

基于处理后的簇中的最多读段匹配的TE参考序列类型，确定所述逆转座子插入的类型，以及任选的，

基于处理后的簇中支持该类型逆转座子插入的读段中的方向一致的多数派读段的方向是否与所述TE参考序列一致，判断所述逆转座子插入的方向。

5.权利要求1-4任一方法，其特征在于，检测所述逆转座子插入多态性还包括确定所述逆转座子插入的位置，确定所述逆转座子插入的位置包括获取候选断点集和筛选所述候选断点集，所述候选断点集包含多个候选断点，

获取所述候选断点集包括，

将所述处理后的簇中的读段与所述基因组参考序列比对，将获得的比对结果中的包含至少一个长度不小于10bp的割裂部分的一级读段中的割裂位置加入到所述候选断点集，定义所述一级读段的大于10bp的割裂部分为二级读段，

将所述二级读段比对到所述基因组参考序列，将获得的比对结果中的包含至少一个长度不小于10bp的割裂部分的二级读段中的割裂位置加入到所述候选断点集，定义所述二级读段的大于10bp的割裂部分为三级读段，

如此，直至N级读段的任一割裂部分小于10bp，其中，

N为自然数，N≥2。

6.权利要求5的方法，其特征在于，获取所述候选断点集还包括，

将i级读段与所述基因组参考序列的比对结果中的满足任一以下(a)-(c)的错配位置加入到所述候选断点集，

(a)所述错配位置包含多聚序列，所述多聚序列包含的碱基满足：6个碱基中至少有5个是一样的碱基，且所述一样的碱基为A或T，

(b)在所述比对结果中支持所述错配的读段的比例小于1/2，

(c)增加和/或删除所述错配位置中的预定碱基数后，将所述i级读段与所述基因组参考序列重新比对，获得的比对结果显示所述i级读段为割裂读段，所述预定碱基数为所述比对结果中的插入和/或缺失的碱基的数目，

所述i级读段包括所述一级读段、所述二级读段至所述N级读段中的至少之一。

7.权利要求6的方法，其特征在于，筛选所述候选断点集包括，

过滤掉处理后的簇中匹配读段的数目不大于2的候选断点，从所述候选断点的位置开始所述匹配读段的长度不小于10bp，所述匹配读段为匹配到所述候选断点的读段，

当任一所述匹配读段都不包含所述多聚序列时，过滤掉处理后的簇中支持读段的比例不大于0.5的候选断点，所述支持读段的比例＝支持所述候选断点的读段数目/匹配到所述断点的读段数目。

8.权利要求5-7任一方法，其特征在于，还包括，

将满足以下的逆转座子插入判定为带双端靶位点重复序列，

两相邻j级读段在所述基因组参考序列上的距离小于30bp，所述j级读段包括所述所述二级读段至所述N级读段中的至少之一。

9.权利要求5-8任一方法，其特征在于，还包括，

通过以下公式计算所述逆转座子插入的杂合概率P(Heter|N_C，N_S)，

$P\left(\mathrm{Heter}\right|N_C,N_S)=\frac{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)}{B(N_S,N_C+N_S,P_{\mathrm{HE}})F\left(\mathrm{HE}\right)+B(N_S,N_C+N_S,P_{\mathrm{HO}})F\left(\mathrm{HO}\right)},$

当MSL大于1则判定所述逆转座子插入为杂合，当MSL小于1则判定所述逆转座子插入为纯合，其中，

MSL＝-log₂(1-P(Heter|N_C，N_S))，

N_C为比对上所述逆转座子插入位置中的截断读段的数目，

N_S比对上所述逆转座子插入位置中的跨越读段的数目，

P_HE＝0.6和P_HO＝0.1分别是期望杂合概率和期望纯合概率，

F(HE)和F(HO)分别是杂合和纯合的权重值，F(HE)＝F(HO)＝0.5，

B(N_S，N_C+N_S，P_HE)是二项分布的概率，表示指定事件发生的概率为P_HE，在N_C+N_S次独立重复试验中该事件发生Ns次的概率。

10.一种检测逆转座子插入多态性的装置，其特征在于，包括，

数据输入单元，用于输入数据；

数据输出单元，用于输出数据；

存储单元，用于存储数据，其中包括可执行的程序；

处理器，与所述数据输入单元、数据输出单元及存储单元数据连接，用于执行所述可执行的程序，所述程序的执行包括完成权利要求1-9任一方法。