本发明属于基于二代测序数据的变异检测方法，涉及一种基于二代测序数据的反转相关复杂变异检测方法。
背景技术：
与人类参考基因组(reference)相比，每个人的基因组都会存在一定的不同，这些不同被称为变异。变异因大小不同可以分为snp(singlenucleotidepolymorphism)和结构变异(structuralvariants)。
不同人体中广泛存在着各种类型的结构变异(structuralvariants)，其中简单变异主要分为删除(deletion)、插入(insertion)、重复(duplication)和反转(inversion)等类型。其中与本发明相关的删除(deletion)指的是与参考基因组相比缺少了一段本该有的基因，重复(duplication)指的是与参考基因组相比将某一段基因多复制了一次或者几次，反转(inversion)指的是与参考基因组相比将某一段基因的方向颠倒了。这些结构变异(structuralvariants)中，有些就目前已知信息评估尚不明确与什么相关，有些是与疾病有直接或者间接关系，例如小儿常见的猫叫综合征是由于删除(deletion)造成，甲型血友病与反转(inversion)相关。
dna序列由a、g、c、t四种碱基组成，测序得到的read实际上就是通过仪器得到由这四种碱基组成的一条一定长度的字符串，以成对的、朝向相反并且有一定范围距离的形式存在可读的bam文件中。由于测序技术的发展，read长度从早年的十几bp发展到一百多甚至两百多bp，想把不同人的很多条read和数量级为三十亿的参考基因组(reference)进行对比找到相符合的位置，在算法层面来看难度其实是相当高的。
目前以二代测序数据为基础检测变异的方法很多，主要使用的信号有assembly、readpair、readdepth和splitread信号等。其中，assembly通过将bam文件中的read进行组装重新比对得到相应的信息；readpair信号指的是在bam文件中每对read之间的距离信息与方向信息，即insertsize和pairorientation信息，能够反应出一些较大的结构变异(structuralvariants)的相关信息；readdepth信号指每个参考基因组位点上覆盖的read的数量信息，能够较为明显地体现某个区域中read数量会有明显变化的变异，例如删除(deletion)和重复(duplication)等；splitread信号指在bam文件中每对read中有一个不能完全比对到参考基因组时，将其分成两部分或者以上与参考基因组进行比对的位点信息，能够较为精确地体现变异的断点信息。通常，变异检测工具使用其中一种或者结合多种信息来检测变异。
目前，业界已经发现了反转(inversion)的及其相关复杂变异检测结果中变异类型错误、断点不准确甚至完全无法检测的问题，并针对此进行了相关的科学研究和算法开发，但国内外检测反转(inversion)及其相关复杂变异的工具较少，目前只有以下几种：
1)gridss:此方法的亮点在于使用了assembly来进行重比对，因此可以较为精确地报出断点，但是它的缺点也较为明显，它只能报出断点而无法明确地报出变异类型，必须要通过三代数据或者其他相关辅助手段才能进行变异类型的确认。
2)tardis：对readpair信号进行聚类，对splitread信号进行一定条件筛选并对readdepth信号提出纯合/杂合变异的假设，结合三种信号来报出变异类型和变异断点，但它目前只能报出其他类型的复杂变异，不能报出与反转(inversion)相关的复杂变异。
3)svelter:建立无效模型(nullmodel)，对检测到的断点进行聚类和整合确定断点，针对n个断点形成的n-1个block随机分配变异类型并迭代打分，输出打分最高的变异类型，并输出变异类型，但是它的变异类型经常出现错误。
技术实现要素：
为了克服以上技术的缺点，本发明目的是提供一种基于二代测序数据的反转相关复杂变异检测方法。
为实现上述目的，本发明是通过以下技术方案来实现：
一种基于二代测序数据的反转相关复杂变异检测方法，包括以下步骤：
步骤1，在滑动窗口内，根据给定的bam文件与选定的参考基因组进行比对，得到pairorientation异常或者insertsize异常的readpair信号，并以readpair信号对不能完全匹配的read进行splitread信号分析，得到对应的断点匹配情况；
步骤2，针对想要寻找的简单变异和复杂变异，建立splitread信号理论模型；包括反转splitread信号的模型，反转-删除splitread信号的模型以及反转-重复splitread信号的模型；
步骤3，将步骤1中得到的断点匹配情况经过步骤2中建立的模型，如果符合某个模型时，记录下相应的变异类型和位置，再判断是否是可信的变异。
本发明进一步的改进在于，步骤1中，用聚类算法进行readpair信号分析，得到pairorientation异常或者insertsize异常的readpair信号；
用模式增长算法进行splitread信号分析，得到不能完全匹配的read的断点匹配情况。
本发明进一步的改进在于，步骤1的具体过程为：
首先，在给定的bam文件中划定了一个100万bp大小的窗口；
然后，在这个100万bp大小的窗口中，以readpair为单位进行第一次扫描：如果一个readpair的pairorientation和/或insertsize信息异常，记录为一个未定的readpair信号，并进行聚类；
最后，在这个100万bp大小的窗口中，以单个read进行第二次扫描：例如，某个read不能完全比对到reference，则称为reada，那么分成两段，以reada两端到中间的方向在64bp范围内与reference比对，如果不能找到reada的两段和reference比较的最小和最大公共子串，则扩大范围为上次查找范围的四倍范围，并反复进行比对，以找到reada和reference比较的最小和最大公共子串为止，并记录对应的位置信息；如果没有，则不记录。
本发明进一步的改进在于，进行聚类的具体过程为：在未定的readpair附近确定是否有五个及以上和readpair信息一致且位置接近的readpair，如果有，则将这个readpair信号记录，具体包括它的位置和方向。
本发明进一步的改进在于，如果reada在确定的readpair信号中有记录，那么开始splitread分析的位置是相对应的readpair信号的位置；如果没有记录，那么开始splitread分析的位置是reada不能完全比对的位置。
本发明进一步的改进在于，步骤2的具体过程为：
根据反转、反转-删除和反转-重复的理论建立相对应的splitread信号的模型；包括反转splitread信号的模型，反转-删除splitread信号的模型以及反转-重复splitread信号的模型；
splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c两段的方向相反；reada和readd为一对readpair，若readd完全比对到reference上，方向为向前和向后的readd都至少有一个，且断点信息一致；
splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada和readd为一对readpair，若readd完全比对到reference上，方向为向前和向后的readd都至少有一个，且断点信息不一致；
splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada和readd为一对readpair，若readd完全比对到reference上，方向为向前和向后的readd都至少有一个，且断点信息不一致。
本发明进一步的改进在于，步骤3的具体过程为：
首先，将步骤1中最后得到的断点匹配情况经过步骤2中模型的检验，如果符合，记录变异类型、断点信息；
然后，判断可能是同个变异的多个splitread信号中，如果断点位置在该read的左半段的read数量和断点位置在该read的左半段的read数量都大于等于1，那么认为该变异报告是可信的；
最后，输出可信的变异。
与现有技术相比，本发明的有益效果在于：
1)基于二代数据的测序成本较低、正确率较高、数据量较大等特点，本发明以二代测序数据为基础，对bam文件中一定大小的移动窗口中所有read的readpair和splitread信息进行分析和筛选，并检验其是否符合构建的复杂变异模型，从而得到特定复杂变异的位置、类型等信息，同时也可以用三代数据画的dotplot图来对其进行检验和确认。
2)本发明明确定义了简单变异和复杂变异的变异类型，严格、准确地根据理论信号建立了变异模型信号，因此可以很准确地提出变异类型；
3)本发明使用splitread信号，以模式增长算法(patterngrowthalgorithm)寻找字符串的最大最小唯一子串，所以能够很精确地指出变异的位置信息。
附图说明
图1为本发明系统整体结构图。
图2为本发明readpair信号分析聚类图。
图3为本发明splitread信号分析举例图。
图4为本发明简单变异删除(deletion)的建模。
图5为本发明简单变异反转(inversion)的建模。
图6为本发明简单变异重复(duplication)的建模。
图7为本发明复杂变异反转-删除(inversion-deletion)的建模。
图8为本发明复杂变异反转-重复(inversion-duplication)的建模。
具体实施方式
为使本发明实施的目的、技术方案和优点更加清楚，下面结合附图和实施例详细说明本发明的实施方式。
本发明主要探究的是结构变异(structuralvariants)中的一部分，尤其以检测反转(inversion)，反转-删除(inversion-deletion)和反转-重复(inversion-duplication)为重点。
本发明包括以下步骤：
步骤1，在滑动窗口内，根据给定的bam文件与选定的参考基因组进行比对，用聚类算法进行readpair信号分析，得到pairorientation异常或者insertsize异常的readpair信号，为splitread信号分析提供前提信息；用模式增长算法(patterngrowthalgorithm)进行splitread信号分析，得到不能完全匹配的read的断点匹配情况即splitread信号，从而精确地确定变异断点信息；
步骤1的具体过程为：
首先，在给定的bam文件中划定了一个100万bp大小的窗口，每次只研究一个窗口里的read的信息，从而解决reference太大的问题；
然后，在这个100万bp大小的窗口中，以readpair为单位进行第一次扫描：如果一个readpair的pairorientation、insertsize等任一信息有异常，先记录下来为一个未定的readpair信号，并进行聚类，即在这个readpair附近确定是否有五个及以上的和readpair信息一致且位置接近的readpair，如果有，则把这个readpair信号确定地记录下来，具体包括它的位置、方向等，为后续的splitread信号分析提供前提信息；
最后，在这个100万bp大小的窗口中，以单个read进行第二次扫描：例如，某个read不能完全比对到reference，则称为reada，那么分成两段，以reada两端到中间的方向去与reference比对，如果不能找到reada的两段和reference比较的最小和最大公共子串，，则反复扩大范围为上次查找范围的四倍范围内去找，以找到reada和reference比较的最小和最大公共子串为止，并记录下来对应的位置信息；如果没有，则不记录。此时只能确定splitread开始比对的方向，不能确定开始比对的位置。如果reada在确定的readpair信号中有记录，那么开始splitread分析的位置就是相对应的readpair信号的位置；如果没有记录，那么开始splitread分析的位置就是reada不能完全比对的位置。最终得到的是reada分成两段去比对后两段的位置和方向，可以在步骤3作为输入。
步骤2，针对想要寻找的简单变异和复杂变异，即反转(inversion)，反转-删除(inversion-deletion)和反转-重复(inversion-duplication)在splitread信号分析的不同特点，建立精确的splitread信号理论模型；
步骤2的具体过程为：
根据反转(inversion)、反转-删除(inversion-deletion)和反转-重复(inversion-duplication)的理论来建立相对应的readpair信号和splitread信号的模型，即反转splitread信号的模型，反转-删除splitread信号的模型以及反转-重复splitread信号的模型。其中反转(inversion)的信号特点是：readpair信号中pairorientation有异常；splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada的matereadd可以完全比对到reference上，而方向为forward和reverse的readd都至少有一个，且断点信息一致，为反转(inversion)的准确断点。反转-删除(inversion-deletion)的信号特点是：readpair信号中pairorientation有异常；splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada的matereadd可以完全比对到reference上，而方向为forward和reverse的readd都至少有一个，且断点信息不一致，即靠近删除(deletion)一侧的断点信息为删除(deletion)远离反转(inversion)一侧的断点和反转(inversion)远离删除(deletion)的断点，靠近反转(inversion)一侧的断点信息为反转(inversion)的断点，它们有一个公共断点即反转(inversion)远离删除(deletion)的断点。反转-重复(inversion-duplication)的信号特点是readpair信号中pairorientation有异常；splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada的matereadd可以完全比对到reference上，而方向为forward和reverse的readd都至少有一个，且断点信息不一致，即靠近反转(inversion)一侧的断点是反转(inversion)靠近重复的断点和重复(duplication)远离反转(inversion)的断点，靠近重复(duplication)一侧的断点是反转(inversion)的断点，它们有一个公共断点即反转(inversion)靠近重复的断点。
步骤3，将步骤1中得到的即splitread信号经过步骤2中建立的模型，如果符合某个模型，记录下相应的变异类型和位置，再通过判断一定的标准判断是否是可信的变异，最后进行输出。
步骤3的具体过程为：
首先，将步骤1中最后得到的不能完全匹配的read的断点匹配情况即splitread信号经过步骤2中模型的检验，如果符合就记录下变异类型、断点等信息；
然后，判断可能是同个变异的多个splitread信号中，如果断点位置在该read的左半段的read数量和断点位置在该read的左半段的read数量都大于等于1，那么认为该变异报告是可信的；
最后，输出可信的变异。
如图1所示，本发明通过读取bam文件中包含各个位点信息的read，并对其进行readpair信号和splitread信号的分析和筛选，观察是否有满足要求的信号满足根据理论情况建立的模型，并对其进行筛选后得到输出。这个详细的过程下面会进行说明。
实施例1
首先获取输入文件即bam文件的信息。人类基因组的长度很长，大约30亿bp左右，即很长的字符串信息，例如图3中的dataset就是一个简化的示例。在使用二代测序技术对样本进行测量时，随着技术的进步，一次测量的长度由几十bp到现在的几百bp，图3中的pattern就是一个简化的示例。因此，可以把在bam文件中保存的信息看作是互相之间有重叠的、长度大小为几百的字符串信息。
然后对每个字符串进行readpair分析和splitread分析。
在100万bp的滑动窗口内，在正常无变异的理想情况下，两个read即字符串为一对，都可以完全比对到参考基因组上，方向相反，距离在某个正常范围内。如图2所示，如果在一定筛选条件下可以认定为是同位置的字符串有不少于5对的虽然可以完全比对到参考基因组(reference)上，但是它们之间的距离超过了正常范围或者方向出现同向的异常信号时，就把这个异常的readpair信号记录下来，作为后续的splitread分析的辅助信息，这个一定的筛选条件是指在read开始位置相差不超过read长度的一定百分比(例如10％)的情况下。
而当两个read中只有一个字符串可以完全比对到参考基因组(reference)时，将另一个不能完全比对上的字符串从两端向内开始用模式增长算法(patterngrowthalgorithm)分别与参考基因组(reference)进行比对，即split成两段进行比对，一直到找到有唯一子字符串的位置停止比对。如果在上一步readpair分析中有相应的信息，那么使用readpair中的位置信息去进行splitread信息的比对，如果有，记录下比对的位置和方向，如果没有，不进行记录；如果在上一步readpair分析中没有相应的信息，那么直接以正常比对的位置为中心，以64bp为半径去进行splitread信息的比对，如果没有满足条件的信息，将寻找半径扩大为上次比对的4倍去寻找，反复3次后，如果有，记录下比对的位置和方向，如果没有，不进行记录。如图3所示，模式增长算法是splitread信息分析的基础，目的是查找给定模式的最小和最大唯一子字符串，例如有数据集s:atcaagtaatgcttagc，模式p:atgca，那么模式p的最小唯一子串为’atg’，最大唯一子串为’atgc’，因为‘at’有两个比对位置，不符合唯一子字符串要求。
接着首先对简单变异进行建模。删除(deletion)、反转(inversion)和重复(duplication)的基本splitread模型如图4，图5和图6所示。此时将上一步得到的splitread信息进行排序、整合和筛选，对符合某种变异模型的字符串按类别进行存储后，根据断点位置等信息对它们继续进行整合判断是否符合模型的其他要求，例如同一个变异中splitread的断点位置是否在read的左半侧和右半侧都有，完全比对上的字符串的方向正负都需要有。最后可以得到简单变异的输出。
不同于简单变异要求两断点一致且完全比对上的字符串的正负方向都需要有，如图7所示为反转-删除(inversion-deletion)的模型，此时反转-删除(inversion-deletion)要求的信号是：两侧的信号类型为反转(inversion)，但大小不同；删除(deletion)一侧的变异大小是反转(inversion)加上删除(deletion)信号大小，反转(inversion)一侧的变异大小是反转(inversion)信号大小，且它们在反转(inversion)一侧有共同的断点。
如图8所示的是反转-重复(inversion-duplication)的模型，此时反转-重复(inversion-duplication)要求的信号是：两侧的信号类型为反转(inversion)，但大小不同；反转(inversion)一侧的变异大小是反转(inversion)和正常区域间的间隔大小，重复(duplication)一侧的变异大小是重复(duplication)的大小加上反转(inversion)和正常区域间的间隔大小，且它们在反转(inversion)一侧有共同的断点。
如果有符合模型的变异输出，就能得到该变异的染色体编号、位置、大小、类型、在左右断点支持的read的数量，从而也可以粗略地得到变异的可信与否。
本发明重点关注的是对反转(inversion)及其相关的复杂变异的检测。从算法角度来说，由于反转(inversion)信号本身难以发掘，加上其多发于基因组上的复杂区域或者重复区域，并且常与其他类型的变异同时出现形成复杂变异，所以对它的检测本身是很困难的；从临床应用角度来说，对反转(inversion)和相关复杂变异的检测与研究可能在疾病的预判和治疗中发挥重要的作用，因此如果能针对其发现重要靶点进行针对性治疗，甚至能够为临床相关应用打下必要的基础，这也是本发明工作的未来意义所在。
1.一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，包括以下步骤：
步骤1，在滑动窗口内，根据给定的bam文件与选定的参考基因组进行比对，得到pairorientation异常或者insertsize异常的readpair信号，并以readpair信号对不能完全匹配的read进行splitread信号分析，得到对应的断点匹配情况；
步骤2，针对想要寻找的简单变异和复杂变异，建立splitread信号理论模型；包括反转splitread信号的模型，反转-删除splitread信号的模型以及反转-重复splitread信号的模型；
步骤3，将步骤1中得到的断点匹配情况经过步骤2中建立的模型，如果符合某个模型时，记录下相应的变异类型和位置，再判断是否是可信的变异。
2.根据权利要求1所述的一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，步骤1中，用聚类算法进行readpair信号分析，得到pairorientation异常或者insertsize异常的readpair信号；
用模式增长算法进行splitread信号分析，得到不能完全匹配的read的断点匹配情况。
3.根据权利要求1所述的一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，步骤1的具体过程为：
首先，在给定的bam文件中划定了一个100万bp大小的窗口；
然后，在这个100万bp大小的窗口中，以readpair为单位进行第一次扫描：如果一个readpair的pairorientation和/或insertsize信息异常，记录为一个未定的readpair信号，并进行聚类；
最后，在这个100万bp大小的窗口中，以单个read进行第二次扫描：例如，某个read不能完全比对到reference，则称为reada，那么分成两段，以reada两端到中间的方向在64bp范围内与reference比对，如果不能找到reada的两段和reference比较的最小和最大公共子串，则扩大范围为上次查找范围的四倍范围，并反复进行比对，以找到reada和reference比较的最小和最大公共子串为止，并记录对应的位置信息；如果没有，则不记录。
4.根据权利要求3所述的一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，进行聚类的具体过程为：在未定的readpair附近确定是否有五个及以上和readpair信息一致且位置接近的readpair，如果有，则将这个readpair信号记录，具体包括它的位置和方向。
5.根据权利要求3所述的一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，如果reada在确定的readpair信号中有记录，那么开始splitread分析的位置是相对应的readpair信号的位置；如果没有记录，那么开始splitread分析的位置是reada不能完全比对的位置。
6.根据权利要求1所述的一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，步骤2的具体过程为：
根据反转、反转-删除和反转-重复的理论建立相对应的splitread信号的模型；包括反转splitread信号的模型，反转-删除splitread信号的模型以及反转-重复splitread信号的模型；
splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c两段的方向相反；reada和readd为一对readpair，若readd完全比对到reference上，方向为向前和向后的readd都至少有一个，且断点信息一致；
splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada和readd为一对readpair，若readd完全比对到reference上，方向为向前和向后的readd都至少有一个，且断点信息不一致；
splitread信号中，如果reada不能完全比对到reference，reada分成b、c两段后比对到reference上时，b和c的方向相反；reada和readd为一对readpair，若readd完全比对到reference上，方向为向前和向后的readd都至少有一个，且断点信息不一致。
7.根据权利要求1所述的一种基于二代测序数据的反转相关复杂变异检测方法，其特征在于，步骤3的具体过程为：
首先，将步骤1中最后得到的断点匹配情况经过步骤2中模型的检验，如果符合，记录变异类型、断点信息；
然后，判断可能是同个变异的多个splitread信号中，如果断点位置在该read的左半段的read数量和断点位置在该read的左半段的read数量都大于等于1，那么认为该变异报告是可信的；
最后，输出可信的变异。
技术总结