CN104298892B - 基因融合的检测装置和方法 - Google Patents

Abstract

本发明公开了一种基因融合的检测装置和方法。该检测方法在序列比对和局部组装验证之间，还包括序列局部聚类的步骤，该步骤包括将与第一位置部分匹配且与参考基因组第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；将一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；将组装序列再次比对到参考基因组上，查看组装序列的比对结果是否与第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明第一位置和第二位置即为基因融合位置。该检测方法真阳性率高，结果更可靠。

Description

基因融合的检测装置和方法

技术领域

本发明涉及生物信息领域，具体而言，涉及一种基因融合的检测装置和方法。

背景技术

在农业经济物种中，特别是具有重要经济价值的农业动植物，基因组结构的差异性导致了不同个体性状的差异性，其中，基因融合现象扮演着重要的角色。基因融合现象是指染色体上两个易位的基因嵌合在一起，形成一个嵌合基因的现象。这种现象一般是由于染色体发生易位、缺失或者倒置造成的。基因融合是生物体内重要的变异类型，这种变异可能会引起基因的过度表达等从而导致生物体性状发生极大变化，并且这些变化通常是不利的。

因此，准确的检测个体基因组中的基因融合现象能够为农业动植物遗传种质资源的进一步发掘利用以及分子育种的实现提供新的平台，尤其是对于多基因关联的复杂性状(如数量性状)的研究具有重要的资源性意义和巨大的应用潜力。

传统基因融合研究方法存在通量低、操作复杂、不便于大规模样品筛查的缺点，而借助于具有通量高、成本低、检测精度高和检测范围广的DNA高通量测序技术，能够低成本、快速检测高深度测序序列支持的基因重排导致的基因融合。因此，目前出现了两种主要的基于高通量测序数据的基因融合的检测方法。

第一种，基于双末端(Pair End，PE)关系的检测方法：由于高通量测序文库构建时插入大小是确定的，那么如果根据PE测序所得到的序列(reads)的比对位置所判定的插入大小，显著偏离了测序文库构建时的插入大小的平均值(例如，一对reads分别比对到不同的染色体上)，则有可能是发生了基因融合。此类方法主要利用这样的双末端关系来判断基因重排导致的异常双末端比对序列(reads)，根据这些序列(reads)的比对位置、插入大小等信息来检测融合。

第二种，基于截断比对(split-mapping)的检测方法：主要利用非完全比对序列(soft-clipped reads—软截断序列)的序列信息进行融合断点识别，然后对断点上下游比对的reads做聚类分析以及拼接组装，最后重新对序列做定位分析，进而检测基因融合现象。

上述两种检测方法中，基于PE关系的检测方法，是根据异常双末端比对序列(reads)的信息来进行融合检测(例如，Break Dancer)，只能大致给出融合位置，而不能确定准确的断点信息，并且此类方法的假阳性较高。如图1所示，每对方向相对的黑线及中间的虚线表示测序文库中插入片段的大小，箭头方向表示文库测序时的方向，方向相对的一对黑线表示双端测序所得的一对序列(reads)，图1中显示了测序所得的多对序列(reads)中每一对序列(reads)的分别与参考基因组序列进行比对时，能够比对到参考基因组不同的染色体位置上，但无法确定准确的断点位置。

而基于截断比对的检测方法，是利用非完全比对的序列信息确定断点位置进而判定融合的方法，或者并未进行局部组装从而没有对融合断点进行验证(例如，Pindel)，或者没有充分利用序列比对信息造成组装耗时(例如，CREST)。因此，仍需对上述基因融合的检测方法进行改进，以提高检测的精确度。

发明内容

本发明旨在提供一种基因融合的检测的装置和方法，以改善现有技术中检测假阳性率较高的缺陷。

为了实现上述目的，根据本发明的一个方面，提供了一种基因融合的检测方法，该检测方法包括序列比对步骤、局部组装验证步骤，检测方法在序列比对步骤和局部组装验证步骤之间，还包括序列局部聚类的步骤，序列局部聚类的步骤包括：根据序列比对步骤得到的与参考基因组上第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列，将第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；局部组装验证的步骤包括：将一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；将组装序列再次比对到参考基因组上，查看组装序列的比对结果是否与第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明第一位置和第二位置即为基因融合位置。

进一步地，在序列局部聚类步骤中，将第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇的步骤为：将第二A类序列和第二B类序列中与第一位置偏差不超过3bp的序列进行聚类，得到基于第一位置的一系列第一位置聚类簇；将第二A类序列中与第二位置偏差不超过3bp的序列进行聚类，得到基于第二位置的一系列第二位置聚类簇，同时将第二B类序列中不能与第二位置进行匹配的部分序列分别与一系列第二位置聚类簇中的序列进行比对，若能与一系列第二位置聚类簇中的某一个第二位置聚类簇的序列相匹配，则将第二B类序列归于某一个第二位置聚类簇中。

进一步地，该检测方法在序列局部聚类之后，以及在进行局部组装步骤之前，还包括对一系列第二位置聚类簇进行过滤的步骤，过滤步骤包括：将一系列第二位置聚类簇中所包含的第二A类序列和第二B类序列的数目之和小于设定值的第二位置聚类簇过滤掉，得到一次过滤后的一系列第二位置聚类簇；根据第二A类序列中的第一位置和第二位置的比对结果所判定的第一位置和第二位置之间的连接方式，将一次过滤后的一系列第二位置聚类簇中与连接方式不一致的第二位置聚类簇过滤掉，得到二次过滤后的一系列第二位置聚类簇。

进一步地，在对一系列第二位置聚类序列进行过滤步骤后，以及进行局部组装的步骤之前，还包括PE验证步骤，PE验证步骤包括：对于二次过滤后的一系列第二位置聚类簇中的任一第二位置聚类簇来说，提取任一第二位置聚类簇中与其所对应的第一位置和第二位置存在PE关系的序列，若根据存在PE关系的序列与参考基因组的比对结果所确定的第一位置和第二位置之间的连接方式，与任一第二位置聚类簇中的序列所确定的第一位置和第二位置之间的连接方式一致，则将存在PE关系的序列记为结构异常PE支持簇；当结构异常PE支持簇中序列的数目达到设定阈值时，进一步验证所二次过滤后的一系列第二位置聚类簇中的序列确实支持第一位置和第二位置为基因融合位置。

进一步地，在局部组装验证步骤后还包括对基因融合的变异类型进行解析的步骤，对基因融合的变异类型进行解析的步骤包括：将二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，第一位置和第二位置分别记为异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端的位置相距不超过10bp，则将这两个异常连接记为相关异常连接；将多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇；在多组异常连接簇中寻找异常种子，异常种子是指两个位于不同异常连接的距离不超过10bp的两个端，当异常连接簇中没有异常种子，则基因融合的变异的类型为缺失；当异常连接簇中存在一对异常种子且一对异常种子分别来自两条异常连接的两端，则基因融合的变异类型为转向；当异常连接簇中存在一个异常种子且形成异常种子的一对异常连接的未形成种子的一端位于与异常种子相同的染色体上，则基因融合的变异的类型为插入、易位或拷贝数目变异。

进一步地，将多组不同的相关异常连接分配到不同的异常连接簇的步骤中，包括过滤属于同一种异常连接的冗余的异常连接，从而得到多组具有变异结构信息的异常连接簇。

进一步地，检测方法在序列比对的步骤之前，还包括序列连接步骤，序列连接步骤包括对双端测序得到的序列中具有重叠序列的片段进行连接，得到较长的连接片段。

进一步地，序列比对的步骤包括：将连接片段和双端测序得到的序列中无重叠序列的非连接片段均与参考基因组序列进行比对，得到连接片段比对结果和非连接片段比对结果；将连接片段比对结果和非连接片段比对结果进行合并，得到总比对结果；将总比对结果进行排序，得到与参考基因组序列第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列。

进一步地，连接片段基于单末端序列与参考基因组序列进行比对，非连接片段基于双末端序列与参考基因组序列进行比对。

根据本发明的另一方面，提供了一种基因融合的检测装置，该检测装置包括比对模块和局部组装模块，检测装置还包括局部聚类模块，其中，局部聚类模块：与比对模块相连，用于将比对模块输出的第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；局部组装模块：与局部聚类模块相连，用于将一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；并将组装序列再次比对到参考基因组上，查看组装序列的比对结果是否与一系列第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明第一位置和第二位置即为基因融合位置。

进一步地，局部聚类模块包括：第一位置聚类簇单元：与比对模块相连，通过将比对模块输出的第二A类序列和第二B类序列中与第一位置偏差不超过3bp的序列进行聚类，得到基于第一位置的一系列第一位置聚类簇；第二位置聚类簇单元：与比对模块相连，通过将比对模块输出的第二A类序列中与第二位置偏差不超过3bp的序列进行聚类，得到基于第二位置的一系列第二位置聚类簇，同时将第二B类序列中不能与第二位置进行匹配的部分序列分别与一系列第二位置聚类簇中的序列进行比对，若第二B类序列能与一系列第二位置聚类簇中的某一个第二位置聚类簇的序列相匹配，则将第二B类序列归于某一个第二位置聚类簇中。

进一步地，检测装置还包括过滤模块，过滤模块包括：一次过滤单元：与局部聚类模块相连，通过将局部聚类模块输出的一系列第二位置聚类簇中所包含的第二A类序列和第二B类序列的数目之和小于设定值的第二位置聚类簇过滤掉，得到一次过滤后的一系列第二位置聚类簇；二次过滤单元：与一次过滤单元相连，并根据第二A类序列中的第一位置和第二位置的比对结果所判定的第一位置和第二位置之间的连接方式，将一次过滤后的一系列第二位置聚类簇中与连接方式不一致的第二位置聚类簇过滤掉，得到二次过滤后的一系列第二位置聚类簇，并将二次过滤后的一系列第二位置聚类簇输入局部组装模块。

进一步地，检测装置还包括PE验证模块，与过滤模块相连，通过读取二次过滤单元输出的二次过滤后的一系列第二位置聚类簇中在第一位置和第二位置存在PE关系的序列，若根据存在PE关系的序列与参考基因组的比对结果所确定的第一位置和第二位置之间的连接方式，与存在PE关系的序列所在的第二位置聚类簇中的序列所确定的第一位置和第二位置之间的连接方式一致，则将存在PE关系的序列记为结构异常PE支持簇；当结构异常PE支持簇中序列的数目达到设定阈值时，进一步验证二次过滤后的一系列第二位置聚类簇中的序列确实支持第一位置和第二位置为基因融合位置；并将二次过滤后的一系列第二位位置聚类簇的序列输入局部组装模块。

进一步地，检测装置还包括基因融合的变异类型解析模块，基因融合的变异类型解析模块包括：异常连接分组单元：与局部组装模块相连，通过将局部组装模块输出的二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，第一位置和第二位置分别记为异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端位置相距不超过10bp，则将这两个异常连接记为相关异常连接；将多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇；变异类型解析单元：与异常连接分组单元相连，通过从异常连接分组单元输出的多组具有变异结构信息的异常连接簇中寻找异常种子，异常种子是指两个位于不同异常连接的距离不超过10bp的两个端，当异常连接簇中没有异常种子，则基因融合的变异的类型为缺失；当异常连接簇中存在一对异常种子且一对异常种子分别来自两条异常连接的两端，则基因融合的变异类型为转向；当异常连接簇中存在一个异常种子且形成异常种子的一对异常连接的各自未形成种子的一端位于相同的染色体上，则基因融合的变异的类型为插入、易位或拷贝数目变异。

进一步地，异常连接分组单元中还包括：相关异常连接子单元：与局部组装模块相连，通过将局部组装模块输出的二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，第一位置和第二位置分别记为异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端位置相距不超过10bp，则将这两个异常连接记为相关异常连接；冗余异常连接过滤子单元，与相关异常连接子单元相连，通过读取相关异常连接子单元中属于同一种异常连接的冗余异常连接，从而得到多组不同的相关异常连接；异常连接分组子单元：与冗余异常连接过滤子单元相连，通过将冗余异常连接过滤子单元输出的多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇，并将多组具有变异结构信息的异常连接簇输入变异类型解析单元。

进一步地，检测装置还包括连接模块：与序列比对模块相连，用于将具有重叠序列的双端测序所得的序列进行连接，得到较长的连接片段，并将连接片段输入序列比对模块中。

进一步地，比对模块包括：比对单元：通过将双端测序得到的序列中无重叠序列的非连接片段和连接模块得到的连接片段分别与参考基因组序列进行比对，得到非连接片段比对结果和连接片段比对结果；合并单元：与比对单元相连，用于将非连接片段比对结果和连接片段比对结果进行合并，得到总比对结果；排序单元：与合并单元相连，用于将总比对结果进行排序，得到与参考基因组序列第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列。

进一步地，比对模块中的比对单元包括：单末端序列比对子单元：与连接模块相连，用于将连接片段与参考基因组序列进行比对，得到连接片段比对结果，并将连接片段比对结果输入合并单元；双末端序列比对子单元：用于将双端测序得到的序列中无重叠序列的非连接片段与参考基因组序列进行比对，得到非连接片段比对结果，并将非连接片段比对结果输入合并单元。

应用本发明的技术方案，通过对现有的基因融合的检测方法进行改进，通过在序列比对步骤后增加序列局部聚类的步骤，充分利用了所有测序所得序列并将这些序列根据比对得到的匹配位置的不同进行聚类，然后将聚类的序列进行组装得到长的组装序列，利用聚类序列组装得到的长的组装序列在与参考基因组再次比对验证基因融合位置是否与序列比对步骤中得到的基因融合位置更准确，能够有效降低假阳性率，为后续分析基因融合的具体变异类型提供了可靠的依据。

附图说明

构成本申请的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1示出了现有技术一种基因融合的检测方法；

图2示出了本发明实施例1中基因融合的检测方法的流程图；

图3示出了本发明实施例2中基因融合的检测方法的流程图；

图4示出了本发明一种优选的实施例中基因融合的检测装置；

图5示出了本发明另一种优选的实施例中基因融合的检测装置；

图6a和图6b示出了本发明的实施例1中比对结果示意图；以及

图7示出了本发明的实施例2中的第二类序列聚类簇示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。

本发明中“参考基因组”是指所待测样品相应物种中已经公开发表的全基因组序列信息；“双末端比对异常的序列”是指根据双末端的序列所得到的片段大小与实际检测到的双末端位置之间的序列大小不一致的序列；“完全匹配”是指双端测序所得的序列(reads)能够全部比对到参考基因组的某一位置上；“部分匹配”是指双端测序所得的序列(reads)比对到参考基因组两个不同的位置上。

正如背景技术部分所提到的，现有技术中检测基因融合的方法存在假阳性率较高、检测耗时长的缺陷，为了改善这一缺陷，在本发明一种典型的实施方式中，提供了一种基因融合的检测方法，如图2所示，该检测方法包括序列比对步骤、局部组装验证步骤，在序列比对步骤和局部组装验证步骤之间，还包括序列局部聚类的步骤，序列局部聚类的步骤包括：根据序列比对步骤得到的与参考基因组上第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列，将第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；局部组装验证的步骤包括：将一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；将组装序列再次比对到参考基因组上，查看组装序列的比对结果是否与第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明第一位置和第二位置即为基因融合位置。

本发明的上述检测方法通过对现有的基因融合的检测方法进行改进，通过在序列比对步骤后增加序列局部聚类的步骤，充分利用了所有测序所得序列并将这些序列根据比对得到的匹配位置的不同进行聚类，然后将聚类的序列进行组装得到长的组装序列，利用聚类序列组装得到的长的组装序列在与参考基因组再次比对验证基因融合位置是否与序列比对步骤中得到的基因融合位置更准确，为后续分析基因融合的具体变异类型的分析提供了可靠的依据。

在本发明的上述检测方法中，上述序列局部聚类步骤中，将第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇的步骤，为降低检测结果的假阳性，上述局部聚类所依据的距离可根据测序深度的不同，设置不同的距离进行聚类。测序深度越深，该距离可设置得相对较大；测序深度越浅，该距离设置相对较小。在本发明一种优选的实施例中，上述序列局部聚类步骤是将第二A类序列和第二B类序列中与第一位置相距不超过3bp的序列进行聚类，得到基于第一位置的一系列第一位置聚类簇；将第二A类序列中与第二位置偏差不超过3bp的序列进行聚类，得到基于第二位置的一系列第二位置聚类簇，同时将第二B类序列中不能与第二位置进行匹配的部分序列分别与一系列第二位置聚类簇中的序列进行比对，若能与一系列第二位置聚类簇中的某一个第二位置聚类簇的序列相匹配，则将第二B类序列归于某一个第二位置聚类簇中。

上述优选实施例中，通过将与第一位置和第二位置偏差不超过3bp的序列进行聚类，使得聚类序列所聚类的位置基本与第一位置和第二位置没有偏差，从而使所支持第一位置和第二位置为基因融合的位置可信度更高。而且本发明的上述序列局部聚类步骤中，不仅对能够比对到第一位置和第二位置的第二A类序列在第一位置和第二位置进行了聚类，而且还对不能比对到第二位置的部分第二B类序列在第二位置进行了聚类，从而使得形成第一位置聚类簇和第二位置聚类簇的序列均是比对的全长序列而非不仅仅是比对匹配上的那部分序列，这样使得参与聚类的序列比较长，从而将聚类序列组装后得到的组装序列较现有技术中的参与组装的序列要更长，进而使比对的结果更准确，假阳性率更低。

在本发明的检测方法中，在序列局部聚类之后，以及在进行局部组装步骤之前，还包括对一系列第二位置聚类簇进行过滤的步骤，设置该步骤的目的是为了使所得到的一系列第二位置聚类簇中的序列的可靠性相对较好，从而将一些可信度较低的序列排除在第二位置聚类簇之外，因此，本领域技术人员可以在现有技术的基础上，设置相应的过滤步骤以提高所得到的第二位置聚类簇的可靠性。

在本发明中，上述过滤步骤包括：将一系列第二位置聚类簇中所包含的第二A类序列和第二B类序列的数目之和小于设定值的第二位置聚类簇过滤掉，得到一次过滤后的一系列第二位置聚类簇；根据第二A类序列中的第一位置和第二位置的比对结果所判定的第一位置和第二位置之间的连接方式，将一次过滤后的一系列第二位置聚类簇中与连接方式不一致的第二位置聚类簇过滤掉，得到二次过滤后的一系列第二位置聚类簇。

上述步骤中的设定值不是一个固定值，而是随不同测序样品的测序深度的不同而有所不同。因此，该过滤步骤的过滤条件在设置时也有所不同。根据测序深度和变异细胞的纯度可以估计出变异序列的覆盖深度，然后根据估计出的覆盖深度选择设定值。例如，在测序深度达到100X的时候，当检测种系突变时可以认为变异细胞的纯度为1，估计出的纯合突变的覆盖深度为100x，杂合突变的覆盖深度为50x，上述设定值可为20到40之间的一个值；同样在测序深度达到100x时，当检测体细胞突变时假设估计出变异细胞的纯度为50％，估计出纯合突变的覆盖深度约为50，杂合突变的覆盖深度约为25，上述设定值可以为10到20之间的一个值。同时，由于真正发生基因融合的两个位置在连接方式是一定的，那么基于第二A类序列比对结果所判定的第一位置和第二位置的连接方式应当与上述过滤后的一系列第二位置聚类簇中的序列所判定的第一位置和第二位置的连接方式一致，而不一致的有可能是测序错误或其他原因所导致的不一致，因而，也需要将这些序列过滤掉，从而得到更接近真实情况的一系列第二位置聚类簇。

相比现有技术，本发明的上述检测方法所检测到的基因融合的位置已经更准确，为了进一步提高本发明的检测方法的可靠性和准确度，在本发明又一种优选的实施例中，上述检测方法在对一系列第二位置聚类序列进行过滤步骤后，以及进行局部组装的步骤之前，还包括PE验证步骤，PE验证步骤包括：对于二次过滤后的一系列第二位置聚类簇中的任一第二位置聚类簇来说，提取任一第二位置聚类簇中与其所对应的第一位置和第二位置存在PE关系的序列，若根据存在PE关系的序列与参考基因组的比对结果所确定的第一位置和第二位置之间的连接方式，与任一第二位置聚类簇中的序列所确定的第一位置和第二位置之间的连接方式一致，则将存在PE关系的序列记为结构异常PE支持簇；当结构异常PE支持簇中序列的数目达到设定阈值时，进一步验证所二次过滤后的一系列第二位置聚类簇中的序列确实支持第一位置和第二位置为基因融合位置。

本发明的上述PE验证步骤，通过从二次过滤后的一系列第二位置聚类簇中抽提所对应的第一位置和第二位置存在PE关系的序列，并将这类存在PE关系的序列所比对到参考基因组上确定的第一位置和第二位置的连接方式与其对应的第二位置聚类簇中的序列所确定的第一位置和第二位置之间连接方式进行相比，若连接方式一致则统计为结构异常PE支持簇，通过统计结构异常PE支持簇中的序列的数目来验证所确定的第一位置和第二位置为基因融合位置的准确性。上述统计得到的结构异常PE支持簇中的序列数目不同，该数目针对不同的测序数据有不同的设定阈值。该设定阈值随测序数据的测序深度的不同而不同，也随测序样品中存在基因融合现象的细胞所占的比例的不同而不同。在测序深度达到100X的时候，当检测种系突变时可以认为变异细胞的纯度为1，估计出的纯合突变的覆盖深度为100x，杂合突变的覆盖深度为50x，上述设定值可为20到40之间的一个值；同样在测序深度达到100x时，当检测体细胞突变时假设估计出变异细胞的纯度为50％，估计出纯合突变的覆盖深度约为50，杂合突变的覆盖深度约为25，上述设定值可以为10到20之间的一个值。

在本发明的检测方法在上述局部组装验证步骤后还包括对基因融合的变异类型进行解析的步骤，该步骤是根据上述步骤所确定的基因融合位置进一步对发生基因融合的具体的结构变异类型进行解析，任何能够根据本发明的上述基因融合位置解析出基因融合的具体结构变异类型的步骤或操作均适用于本发明。在本发明一种优选的实施例中，上述对基因融合的变异类型进行解析的步骤包括：将二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，第一位置和第二位置分别记为异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端的位置相距不超过10bp，则将这两个异常连接记为相关异常连接；将多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇；在多组异常连接簇中寻找异常种子，异常种子是指两个位于不同异常连接的距离不超过10bp的两个端，当异常连接簇中没有异常种子，则基因融合的变异的类型为缺失；当异常连接簇中存在一对异常种子且一对异常种子分别来自两条异常连接的两端，则基因融合的变异类型为转向；当异常连接簇中存在一个异常种子且形成异常种子的一对异常连接的未形成种子的一端位于与异常种子相同的染色体上，则基因融合的变异的类型为插入、易位或拷贝数目变异。

本发明的上述对基因融合的变异类型进行解析的步骤不仅能够提供变异位点的连接情况，而且能够提供变异结构的各种精确信息，包括变异的种类、发生的位点以及长度，当上述基因融合的变异类型是拷贝数目变异、插入或者易位的话，还能提供供体的精确信息。

在上述将多组不同的相关异常连接分配到不同的异常连接簇的步骤中，包括过滤属于同一种异常连接的冗余的异常连接，从而得到多组具有变异结构信息的异常连接簇。该步骤是过滤部分重复的异常连接，以减少数据量，从而提高检测速度。

在本发明的上述检测方法中，上述序列比对步骤中直接利用双端测序得到的序列与参考基因组进行比对也能实现上述比对结果。在本发明中，连接片段并非是来自两对不同的双末端序列，而是同一对双末端序列的两条序列基于其之间的重叠部分进行连接，主要考虑到如果测序文库的插入片段较小，PE读长较长，则存在测通的情况(例如文库150bp，PE100测序)，如此便可基于中间的重叠序列将同一对测序序列(reads)连接成更长的序列，从而之后进行比对和组装。因此，在本发明的上述序列比对的步骤之前，还包括序列连接步骤，该步骤包括对双端测序得到的序列中具有重叠序列的片段进行连接，得到较长的连接片段。基于较长的连接片段比对得到结果更可靠。

在本发明另一种优选的实施例中，上述序列比对的步骤包括：将连接片段和双端测序得到的序列中无重叠序列的非连接片段均与参考基因组序列进行比对，得到连接片段比对结果和非连接片段比对结果；将连接片段比对结果和非连接片段比对结果进行合并，得到总比对结果；将总比对结果进行排序，得到与参考基因组序列第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列。

本发明的上述对比步骤通过利用部分连接片段进行比对，使得用于比对的序列更长，使比对结果更准确，便于融合位置的检测；再对比对结果进行合并和排序，并可更直观、更清晰地体现与第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列。

在本发明的上述检测方法中，由于连接片段是基于同一对双末端序列的两条序列之间的重叠部分进行连接。此处主要是考虑到如果文库插入片段较小，PE读长较长，则存在测通的情况(例如文库150bp，PE100测序)，如此便可基于中间的重叠序列将同一对reads连接成更长的序列，从而之后进行比对、组装等。因而对连接片段基于该一对双末端序列中的其中一条序列进行比对，而无重叠序列仍基于其测序所得的双末端序列进行比对。

在本发明另一种典型的实施方式中，提供了一种基因融合的检测装置，如图4所示，该检测装置包括比对模块和局部组装模块，还包括局部聚类模块，其中，局部聚类模块：与比对模块相连，用于将比对模块输出的第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；局部组装模块：与局部聚类模块相连，用于将一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；并将组装序列再次比对到参考基因组上，查看组装序列的比对结果是否与一系列第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明第一位置和第二位置即为基因融合位置。

本发明的上述检测装置，通过对现有的基因融合的检测装置进行改进，通过在序列比对模块和局部组装模块之间增加局部聚类模块，充分利用序列比对模块得到的比对信息，并根据比对得到的基因融合的第一位置和第二位置，进行有针对性地进行局部聚类，然后将聚类的序列进行组装得到长的组装序列，并利用该长的组装序列与参考基因组再次比对验证基因融合位置是否与序列比对步骤中得到的基因融合位置更准确，为后续分析基因融合的具体变异类型的分析提供了可靠的依据。

在本发明的上述局部聚类模块中，为降低检测结果的假阳性，上述局部聚类装置中所依据的距离可根据测序深度的不同，设置不同的距离进行聚类。测序深度越深，该距离可设置得相对较大；测序深度越浅，该距离设置相对较小。在本发明一种优选的实施例中，上述局部聚类模块包括：第一位置聚类簇单元：与比对模块相连，通过将比对模块输出的第二A类序列和第二B类序列中与第一位置偏差不超过3bp的序列进行聚类，得到基于第一位置的一系列第一位置聚类簇；第二位置聚类簇单元：与比对模块相连，通过将比对模块输出的第二A类序列中与第二位置偏差不超过3bp的序列进行聚类，得到基于第二位置的一系列第二位置聚类簇，同时将第二B类序列中不能与第二位置进行匹配的部分序列分别与一系列第二位置聚类簇中的序列进行比对，若第二B类序列能与一系列第二位置聚类簇中的某一个第二位置聚类簇的序列相匹配，则将第二B类序列归于某一个第二位置聚类簇中。

上述优选实施例中，通过将与第一位置和第二位置偏差不超过3bp的序列进行聚类，使得聚类序列所聚类的位置基本与第一位置和第二位置没有偏差，从而使所支持第一位置和第二位置为基因融合的位置可信度更高。而且本发明的上述序列局部聚类模块中，不仅对能够比对到第一位置和第二位置的第二A类序列在第一位置和第二位置进行了聚类，而且还对不能比对到第二位置的部分第二B类序列在第二位置进行了聚类，从而使得形成第一位置聚类簇和第二位置聚类簇的序列均是比对的全长序列而非不仅仅是比对匹配上的那部分序列，这样使得参与聚类的序列比较长，从而使后续的局部组装模块得到的组装序列较现有技术中的组装序列更长，进而使比对的结果更准确，假阳性率更低。

在本发明的检测装置中，在局部聚类模块之后及局部组装模块之前，还可以设置对一系列第二位置聚类簇进行过滤的模块，设置该模块的目的是为了使所得到的一系列第二位置聚类簇中的序列的可靠性相对较好，从而将一些可信度较低的序列排除在第二位置聚类簇之外，因此，本领域技术人员可以在现有技术的基础上，设置相应的过滤模块以提高所得到的第二位置聚类簇的可靠性。

在本发明中，上述过滤模块包括：一次过滤单元和二次过滤单元，一次过滤单元与局部聚类模块相连，通过将局部聚类模块输出的一系列第二位置聚类簇中所包含的第二A类序列和第二B类序列的数目之和小于设定值的第二位置聚类簇过滤掉，得到一次过滤后的一系列第二位置聚类簇；二次过滤单元与一次过滤单元相连，并根据第二A类序列中的第一位置和第二位置的比对结果所判定的第一位置和第二位置之间的连接方式，将一次过滤后的一系列第二位置聚类簇中与连接方式不一致的第二位置聚类簇过滤掉，得到二次过滤后的一系列第二位置聚类簇，并将二次过滤后的一系列第二位置聚类簇输入局部组装模块。

上述过滤模块中的设定值不是一个固定值，而是随不同测序样品的测序深度的不同而有所不同。因此，该过滤模块的过滤条件在设置时也有所不同。例如，在测序深度达到100X的时候，当检测种系突变时可以认为变异细胞的纯度为1，估计出的纯合突变的覆盖深度为100x，杂合突变的覆盖深度为50x，上述设定值可为20到40之间的一个值；同样在测序深度达到100x时，当检测体细胞突变时假设估计出变异细胞的纯度为50％，估计出纯合突变的覆盖深度约为50，杂合突变的覆盖深度约为25，上述设定值可以为10到20之间的一个值。

同时，由于真正发生基因融合的两个位置在连接方式是一定的，那么基于第二A类序列比对结果所判定的第一位置和第二位置的连接方式应当与上述过滤后的一系列第二位置聚类簇中的序列所判定的第一位置和第二位置的连接方式一致，而不一致的有可能是测序错误或其他原因所导致的不一致，因而，也需要将这些序列过滤掉，从而得到更接近真实情况的一系列第二位置聚类簇。

相比现有技术，本发明的上述检测装置所检测到的基因融合的位置已经更准确，为了进一步提高本发明的检测装置的可靠性和准确度，在本发明又一种优选的实施例中，如图5所示，上述检测装置还包括PE验证模块，该PE验证模块与过滤模块相连，通过读取二次过滤单元输出的二次过滤后的一系列第二位置聚类簇中在第一位置和第二位置存在PE关系的序列，若根据存在PE关系的序列与参考基因组的比对结果所确定的第一位置和第二位置之间的连接方式，与存在PE关系的序列所在的第二位置聚类簇中的序列所确定的第一位置和第二位置之间的连接方式一致，则将存在PE关系的序列记为结构异常PE支持簇；当结构异常PE支持簇中序列的数目达到设定阈值时，进一步验证所二次过滤后的一系列第二位置聚类簇中的序列确实支持第一位置和第二位置为基因融合位置；并将二次过滤后的一系列第二位位置聚类簇的序列输入局部组装模块。

本发明的上述PE验证模块，通过从二次过滤后的一系列第二位置聚类簇中抽提所对应的第一位置和第二位置存在PE关系的序列，并将这类存在PE关系的序列所比对到参考基因组上确定的第一位置和第二位置的连接方式与其对应的第二位置聚类簇中的序列所确定的第一位置和第二位置之间连接方式进行相比，若连接方式一致则统计为结构异常PE支持簇，通过统计结构异常PE支持簇中的序列的数目来验证所确定的第一位置和第二位置为基因融合位置的准确性。上述统计得到的结构异常PE支持簇中的序列数目不同，该数目针对不同的测序数据有不同的设定阈值。该设定阈值随测序数据的测序深度的不同而不同，也与测序样品中存在基因融合现象的细胞所占的比例的不同而不同。在测序深度达到100X的时候，当检测种系突变时可以认为变异细胞的纯度为1，估计出的纯合突变的覆盖深度为100x，杂合突变的覆盖深度为50x，上述设定值可为20到40之间的一个值；同样在测序深度达到100x时，当检测体细胞突变时假设估计出变异细胞的纯度为50％，估计出纯合突变的覆盖深度约为50，杂合突变的覆盖深度约为25，上述设定值可以为10到20之间的一个值。

本发明的检测装置还包括基因融合的变异类型解析模块，该模块是根据上述模块中所确定的基因融合位置进一步对发生基因融合的具体的结构变异类型进行解析，任何能够根据本发明的上述基因融合位置解析出基因融合的具体结构变异类型的模块均适用于本发明。在本发明一种优选的实施例中，上述基因融合的变异类型解析模块包括：异常连接分组单元：与局部组装模块相连，通过将局部组装模块输出的二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，第一位置和第二位置分别记为异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端位置相距不超过10bp，则将这两个异常连接记为相关异常连接；将多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇；变异类型解析单元：与异常连接分组单元相连，通过从异常连接分组单元输出的多组具有变异结构信息的异常连接簇中寻找异常种子，异常种子是指两个位于不同异常连接的距离不超过10bp的两个端，当异常连接簇中没有异常种子，则基因融合的变异的类型为缺失；当异常连接簇中存在一对异常种子且一对异常种子分别来自两条异常连接的两端，则基因融合的变异类型为转向；当异常连接簇中存在一个异常种子且形成异常种子的一对异常连接的各自未形成种子的一端位于相同的染色体上，则基因融合的变异的类型为插入、易位或拷贝数目变异。

本发明的上述基因融合的变异类型解析模块不仅能够提供变异位点的连接情况，而且能够提供变异结构的各种精确信息，包括变异的种类、发生的位点以及长度，当上述基因融合的变异类型是拷贝数目变异、插入或者易位的话，还能提供供体的精确信息。

在本发明另一种优选的实施例中，上述异常连接分组单元中还包括：相关异常连接子单元、冗余异常连接过滤子单元和异常连接分组子单元，相关异常连接子单元与局部组装模块相连，通过将局部组装模块输出的二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，第一位置和第二位置分别记为异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端位置相距不超过10bp，则将这两个异常连接记为相关异常连接；冗余异常连接过滤子单元与相关异常连接子单元相连，通过读取相关异常连接子单元中属于同一种异常连接的冗余异常连接，从而得到多组不同的相关异常连接；异常连接分组子单元与冗余异常连接过滤子单元相连，通过将冗余异常连接过滤子单元输出的多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇，并将多组具有变异结构信息的异常连接簇输入变异类型解析单元。

在本发明上述优选的实施例中，上述异常连接分组单元通过设置上述三个子单元，利用冗余异常连接过滤子单元将相关异常连接子单元中所找到的属于同一种异常连接的重复的异常连接进行过滤，以减少数据量，提高检测装置的检测速度。

在本发明的上述检测装置中还包括连接模块，该连接模块与序列比对模块相连，用于将具有重叠序列的双端测序所得的序列进行连接，得到较长的连接片段，并将连接片段输入序列比对模块中。通过利用该连接模块，可以将测序数据中的某些具有重叠序列的测序序列进行连接，以形成更长的连接片段，基于更长的连接片段使得与其相连的序列比对模块中所得到比对结果更可靠。

在本发明另一种优选的实施例中，上述比对模块包括：比对单元、合并单元和排序单元，比对单元通过将双端测序得到的序列中无重叠序列的非连接片段和连接模块得到的连接片段分别与参考基因组序列进行比对，得到非连接片段比对结果和连接片段比对结果；合并单元与比对单元相连，用于将非连接片段比对结果和连接片段比对结果进行合并，得到总比对结果；排序单元与合并单元相连，用于将总比对结果进行排序，得到与参考基因组序列第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列。

本发明的上述对比子单元通过利用部分较长的连接片段进行比对，使得比对结果更准确，便于融合位置的检测；再对比对结果进行合并和排序，并可更直观、更清晰地体现与第一位置完全匹配的第一类序列、与第一位置部分匹配且与参考基因组上第二位置部分匹配的第二A类序列和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列。

由于连接片段是基于同一对双末端序列的两条序列之间的重叠部分进行连接。本发明中，考虑到如果文库插入片段较小，PE读长较长，则存在测通的情况(例如文库150bp，PE100测序)，这样便可基于中间的重叠序列将同一对reads连接成更长的序列，从而利用这种较长的连接片段进行比对、组装等。因而在本发明又一种优选的实施例中，上述比对单元又包括单末端序列比对子单元和双末端序列比对子单元，单末端序列比对子单元与连接模块相连，用于将连接片段与参考基因组序列进行比对，得到连接片段比对结果，并将连接片段比对结果输入合并单元；双末端序列比对子单元用于将双端测序得到的序列中无重叠序列的非连接片段与参考基因组序列进行比对，得到非连接片段比对结果，并将非连接片段比对结果输入合并单元。上述两个子单元分别对较长的连接片段基于某一末端序列进行比对和对无重叠序列的测序所得序列(reads)仍基于其测序所得的双末端序列进行比对，比对速度更快。

下面将结合具体的实施例来进一步说明本发明的有益效果。

下列实施例对模拟的水稻的全基因组测序数据进行基因融合的检测。

实施例1

实施例1按照图2所示的流程，利用BWA序列比对软件对模拟的水稻的全基因组测序数据基于双端序列与水稻参考基因组进行比对，如图6a所示，得到BAM格式的比对结果，示例见图6b。

其中，图6a中的coor表示标尺；ref表示参考序列；r001至r004表示待比对的测序所得序列；图6b表示r001至r004比对到参考序列上之后的结果，第一行和第二行为BAM文件的标题(title)。余下的各列分别表示：1：序列的编号；2：序列比对的FLAG值，用二进制的方式表示，体现序列比对到参考序列上的各种信息，比如可以体现出序列比对的是否成功，比对的正反链信息，配对序列的比对情况等；3：参考序列的信息；4：比对的位置，以r0001序列为例，7表示该序列比对到参考序列的第7个碱基位置处；5：比对质量值，体现与参考序列的匹配情况；6：CIGAR值，“M”表示比对成功，“S”和“H”都表示没有比对在该位置上，“I”表示插入；“D”表示缺失；字母前面的数字表示序列的长度。例如roo1待对比序列在参考序列上的比对状况为：“8M2I4M1D3M”表示序列的前8个碱基比对成功，后2个碱基相比参考序列为插入碱基，再后面的4个碱基与参考序列匹配，然后1个碱基相比参考序列缺失，最后3个碱基也与参考序列匹配；7：配对序列的染色体位置，“＝”表示与参考序列在同一条染色体上，“*”表示配对序列(PE序列)没有比对在参考序列上；8：配对的序列(PE序列如r001/2序列)的比对位置；9：插入片段长度；10：表示测序所得的序列。11.测序质量值，*表示此处没有测序质量值。余下各列随比对结果不同而不同，其中第12列中SA：Z表示软截断(soft-clip)的比对情况。

从上图6a和6b的比对结果可以得到与水稻参考基因组序列第一位置完全匹配的第一类序列r001/2、与第一位置部分匹配且与水稻参考基因组上第二位置部分匹配的第二A类序列r003和与第一位置部分匹配而不与第二位置部分匹配的第二B类序列r002；

根据上述序列比对结果中的第二A类序列所比对到水稻参考基因组的第一位置和第二位置，将上述第二A类序列和第二B类序列基于与第一位置和第二位置的距离不超过3bp进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；

利用CAP3序列组装软件，对上述一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；将组装序列再次比对到水稻参考基因组上，查看比对结果是否与第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明所述第一位置和第二位置即为基因融合位置。

实施例2

实施例2按照图3所示的流程，首先，将模拟的序列中具有多个不同的重叠序列的片段进行连接，得到多个较长的连接片段；将上述连接片段和不具有重叠片段的模拟测序序列均与水稻的参考基因组序列进行比对，得到连接片段比对结果和非连接片段比对结果；将连接片段比对结果和非连接片段比对结果进行合并，得到总比对结果；将总比对结果进行排序，得到与水稻的参考基因组序列第一位置完全匹配的第一类片段和与第一位置不完全匹配的第二类片段；第二类片段包括能够比对到第二位置的第二A类片段和未比对到第二位置的第二B类片段；

将所有的第二类片段的比对结果中按照与第一位置的距离偏差小于3bp进行聚类，得到一系列基于第一位置的第一类聚类簇(1st cluster)。在任意第一类聚类簇中，所有片段的第一位置位置均相近，但是第二A类片段的第二位置有可能相互之间的距离较远。因此需要根据第二A类序列的第二位置进行第二次聚类，生成第二类聚类簇(2ndcluster)。在任意第二类聚类簇中，所有序列的第一位置和第二位置距离均较近。对于任意第一类聚类簇中的第二B类序列，提取出未能成功比对的部分，分别于各个第二类聚类簇中的序列进行比对。假如第二B类的序列能与某个第二类聚类簇比对成功，则将该第二B类序列添加到该第二类聚类簇中。

在进行下一步操作之前，需要先过滤掉不符合条件的第二类聚类簇。首先，统计出每一个第二类聚类簇中所包含的第二类序列的数目，过滤掉序列数目小于3的第二类聚类簇。然后根据第二A类序列第一位置和第二位置的比对情况识别出两个位置之间的连接方式，并要求第二类聚类簇中的所有连接方式一致。经过过滤后的第二类聚类簇就被认为是由于结构变异所导致的异常连接，然后分别使用PE支持和局部组装的方法对该异常连接进行验证。

对任意第二类聚类簇，提取出第一位置和第二位置周围的存在PE关系的序列。假如一对PE序列的比对方式与该第二类聚类簇的连接方式一致，则将该PE序列归类到支持簇中，当支持簇中的PE序列的数目达到10时，则认为上述PE序列支持上述结构变异所导致的异常连接。

将第二类聚类簇中的所有第二类序列的全长提取出来，并使用CAP3软件进行局部组装，并将组装的结果序列比对到水稻参考基因组上。查看组装序列的比对结果是否与第二类聚类簇中第二A类序列的比对情况一致，并过滤掉比对位置不一致的第二类聚类簇。第二类聚类簇实例如下图7所示。

在图7中，第一行的信息分别表示：第一位置(其中序列比对的情况参照图6b中所列内容)、聚类簇中序列的数目、第一位置附近的测序深度、第二位置聚类簇的连接方式以及支持这种连接方式的PE序列的个数。其中，SSSSSSSS…SSSSSSSS表示比对到参考基因组上的序列，由于序列较长，用Sequence的首字母S来表示；QQQQQQQQ…QQQQQQQQ表示序列的测序质量，用Quality的首字母Q来表示。其余三行分别表示：第二A类序列、第二B类序列和PE支持的序列。聚类簇的连接方式用VCF(variat call format)格式的方式表示，其中r1和r2分别表示第一位置和第二位置。第一位置和第二位置之间一共有四种可能的连接方式，详细的含义见下表1。

表1：

连接方式	含义
		r1[r2[	正向r1后面连着正向的r2
r1]r2]	正向r1后面连着反向的r2

]r2]r1	正向r2后面连着正向的r1
		[r2[r1	反向r2后面连着正向的r1

至此得到了所有高度可信的第二类聚类簇，这些聚类簇的实质是由于基因组变异而产生的不正常连接方式。可以产生不正常连接方式的基因组变异有缺失，插入，染色体间或染色体内易位，转向和拷贝数目变异。在该实施例中，无法直接通过某一个第二类聚类簇推断出产生该不正常连接的结构变异，但是可以根据相关的第二类聚类簇做出正确的推断。

第二类聚类簇的第一位置和第二位置可以分别被认为是一个不正常连接的两端，假如一个异常连接的任意一端与另一个异常连接的任意一端距离较近，则认为他们是相关异常连接。将多个不同的相关异常连接分配到不同的异常连接簇中，并且去除表示同一种连接的冗余异常连接。此时，异常连接簇中包含了解析出变异结构的所有信息，由于不同的结构变异具有不同的特征，结构解析时也要根据异常连接簇的特征进行。其中，得到部分异常连接簇如下表2。

表2：

从表2中可以看出，异常连接1对比的第一位置为1号染色体的3051951bp处，第二位置在1号染色体的86982034bp处，第一位置与第二位置的连接方式为r1[r2[，其中，r1代表比对到第一位置的序列，r2代表比对到第二位置的序列；r1[r2[代表比对到第二位置的序列正向连接在比对到第一位置的序列的右边，即：比对在86982034bp处的序列正向连接在比对到3051951bp处的序列之后。而异常连接4对比的第一位置为1号染色体的86982034bp处，第二位置在1号染色体的3051951bp处，第一位置与第二位置的连接方式为]r2]r1，说明比对到3051951bp处的序列正向连接在比对到86982034bp处的序列的前面。同理，异常连接2所判断的异常连接的第一位置和第二位置与异常连接3所判断的异常连接的第二位置和第二位置刚好相反，是从同一序列的两端分别比对的角度来判断得到的结果，其体现的仍是同一种异常连接。可见，异常连接簇中的所有异常连接都存在直接或者间接的关联，由于异常连接1和4，异常连接2和3互为冗余连接，在后续分析中只需保留一组即可。

在结构解析时，首先在异常连接簇中寻找异常种子，异常种子是指两个位于不同异常连接上的距离不超过10bp的两个端。异常种子的实质是结构变异发生时断开的正常连接。不同类型变异的异常种子出现的情况也会不同。假如某个异常连接簇中没有异常种子出现，那么变异类型就推断缺失。假如异常连接簇中存在一对异常种子，他们分别来自两条异常连接的两端，那变异类型就推断为转向。假如包含一个异常种子的一对异常连接的未形成种子的一端位于同一条染色体上，那就推断变异类型为插入、易位或拷贝数目变异。以上三种推断出的结构都要与异常连接簇中的异常连接相一致。对于未能推断出变异结构的变异连接簇标记为“未知(unknown)”。三种变异类型的解析结果分别以下列实例进行说明：

1.缺失，缺失时并不存在异常种子，示例异常连接簇见表3。

表3：

表3表示这个异常连接簇中有两条异常连接，它们都表示1号染色体的1706687号碱基后面正向连接着1号染色体的1843279号碱基。这条连接跨过了1706688bp位置到1843279bp位置之间的碱基，因此该变异类型为缺失。具体说来，在1号染色体上1706687bp位置处缺失了从1706688bp到1843279bp之间的136592bp长度的序列，此预测结果如下表4所示。

表4：

2.转向(又称倒置)，转向时在一对变异连接上出现一对变异种子，具体见表5。

表5：

从表5中可以看出，第一个异常连接表示1号染色体的3478027bp后连着反向的3650098bp，且第二个异常连接表示1号染色体的3478031bp处后连着反向的3650102bp，可以得出：‘1:3478027’和‘1:3478031’为一对变异种子，‘1:3650098’和‘1:3650102’为一对变异种子。根据变异种子的实质可以推断出，‘1:3478027’和‘1:3478031’之间的正常连接断开后分别和‘1:3650098’和‘1:3650102’相连。结合异常连接的方向可以推断出此异常连接簇的实质为1号染色体3478027bp位置到3650102位置之间的碱基发生了倒置的情况。具体的预测结果使用如下表6所示的方式表示：

表6：

3.易位(cut-paste)或插入，易位发生时会出现一个变异种子，具体见表7。

表7：

从表7中可以看出，1号染色体的3404996bp前面连接着正向的1号染色体上的226721129处，1号染色体的226721130bp前面连接着正向的1号染色体3447873。‘1:226721129’和‘1:226721130’组成一个变异种子，说明这两个位点之间的正常连接被断开，之后分别于‘1:3404996’和‘1:3447873’相连。这些信息说明1号染色体上的3404996bp处至3447873bp处的碱基插入到‘1:226721129’和‘1:226721130’之间。具体的预测结果使用如下表8所示方式表示。

表8：

经过结构解析之后，该实施例所提供的不仅是变异位点的连接情况，而且是变异结构的各种精确信息，包括变异的种类、发生的位点、长度，如果是拷贝数变异、插入或者易位的话，还能提供供体的精确信息。

由于本发明的实施例是经PE验证、局部组装验证、以及结构解析时的一致性验证三重验证步骤，因此，最终结果的准确率是非常高的。在使用水稻的模拟数据进行测试时，在预测转向和基因融合是的准确率，即真阳性率高达到87％以上，最高可达98.5％，检出率，即敏感性约为82％，具体数据见下表9，其中，准确率＝真阳性位点数/鉴定出的位点数×100％；检出率＝真阳性位点数/模拟位点数×100％。

表9：

变异类型	模拟位点数	鉴定出的位点数	真阳性位点数	检出率	准确率
						缺失	164	174	152	92.68％	87.36％
倒置	91	81	78	85.71％	96.30％
						插入	155	129	127	81.94％	98.45％

对比例1

CREST软件是本领域常用的寻找基因组结构变异的方法，它通过比对组装-比对-查找-再组装-再比对的方式来寻找结构变异。首先从BAM文件中提取出具有软截断的序列，假定这些序列的比对位置为结构变异的第一断点，并将软截断序列组装成第一重叠群(contig)。然后将第一重叠群比对到基因组上，第一重叠群比对的位置就认为可能是第二断点。提取出比对在第二断点上的具有软截断的序列，使用同样的方法将软截断的序列组装成第二重叠群，并将第二重叠群比对到基因组上。假如第二重叠群的比对位置与第一断点的位置一致，就认为此二位置上存在一个结构变异。

为了比较CREST和本文中介绍的方法的效果差别，我们使用相同的水稻模拟数据测试CREST软件得到的结果见下表10。

表10：

	模拟位点数	鉴定出位点数	真阳性位点数	检出率	准确率
						缺失	164	221	135	82.31％	61.09％
倒置	91	54	47	51.65％	87.03％
						插入	155	231	58	37.42％	25.11％

从以上的描述中，可以看出，与对比例1相比，本发明的实施例1和2通过将待测样本通过测序得到的序列(reads)比对到相应物种已发表的基因组上，将比对结果中不完全比对得上的序列(reads)根据其不能比对上的断点位置确定基因的候选融合位置，提高了检测的准确性(真阳性率)；实施例2进一步通过局部聚类(断点位置相近的reads为同一个类别)和局部组装的结果对候选融合位置进行校正或筛选，从而更准确地、更快速地确定基因的融合位置。本发明的检测方法相比传统的检测方法具有准确率高、敏感性强、精度高、速度快、资源消耗少等显著优势。

显然，本领域的技术人员应该明白，本发明的上述的一些模块或一些步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (18)

1.一种基因融合的检测方法，所述检测方法包括序列比对步骤、局部组装验证步骤，其特征在于，所述检测方法在所述序列比对步骤和所述局部组装验证步骤之间，还包括序列局部聚类的步骤，

所述序列局部聚类的步骤包括：根据所述序列比对步骤得到的与参考基因组上第一位置完全匹配的第一类序列、与所述第一位置部分匹配且与所述参考基因组上第二位置部分匹配的第二A类序列和与所述第一位置部分匹配而不与所述第二位置部分匹配的第二B类序列，将所述第二A类序列和第二B类序列基于与所述第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；

所述局部组装验证的步骤包括：将所述一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；将所述组装序列再次比对到所述参考基因组上，查看所述组装序列的比对结果是否与所述第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明所述第一位置和第二位置即为基因融合位置。

2.根据权利要求1所述的检测方法，其特征在于，在所述序列局部聚类步骤中，将所述第二A类序列和第二B类序列基于与所述第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇的步骤为：

将所述第二A类序列和所述第二B类序列中与所述第一位置偏差不超过3bp的序列进行聚类，得到基于所述第一位置的一系列第一位置聚类簇；将所述第二A类序列中与所述第二位置偏差不超过3bp的序列进行聚类，得到基于所述第二位置的一系列第二位置聚类簇，同时将第二B类序列中不能与所述第二位置进行匹配的部分序列分别与所述一系列第二位置聚类簇中的序列进行比对，若能与所述一系列第二位置聚类簇中的某一个第二位置聚类簇的序列相匹配，则将所述第二B类序列归于所述某一个第二位置聚类簇中。

3.根据权利要求2所述的检测方法，其特征在于，所述检测方法在所述序列局部聚类之后，以及在进行所述局部组装验证步骤之前，还包括对所述一系列第二位置聚类簇进行过滤的步骤，所述过滤步骤包括：

将所述一系列第二位置聚类簇中所包含的第二A类序列和第二B类序列的数目之和小于设定值的第二位置聚类簇过滤掉，得到一次过滤后的一系列第二位置聚类簇；

根据所述第二A类序列中的第一位置和第二位置的比对结果所判定的所述第一位置和所述第二位置之间的连接方式，将所述一次过滤后的一系列第二位置聚类簇中与所述连接方式不一致的第二位置聚类簇过滤掉，得到二次过滤后的一系列第二位置聚类簇。

4.根据权利要求3所述的检测方法，其特征在于，所述检测方法在对所述一系列第二位置聚类序列进行过滤步骤后，以及进行所述局部组装验证的步骤之前，还包括PE验证步骤，所述PE验证步骤包括：

对于所述二次过滤后的一系列第二位置聚类簇中的任一第二位置聚类簇来说，提取所述任一第二位置聚类簇中与其所对应的第一位置和第二位置存在PE关系的序列，若根据所述存在PE关系的序列与所述参考基因组的比对结果所确定的所述第一位置和第二位置之间的连接方式，与所述任一第二位置聚类簇中的序列所确定的所述第一位置和第二位置之间的连接方式一致，则将所述存在PE关系的序列记为结构异常PE支持簇；

当所述结构异常PE支持簇中序列的数目达到设定阈值时，进一步验证所述二次过滤后的一系列第二位置聚类簇中的序列确实支持所述第一位置和第二位置为基因融合位置。

5.根据权利要求4所述的检测方法，其特征在于，所述检测方法在局部组装验证步骤后还包括对基因融合的变异类型进行解析的步骤，所述对基因融合的变异类型进行解析的步骤包括：

将所述二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，所述第一位置和第二位置分别记为所述异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端的位置相距不超过10bp，则将这两个异常连接记为相关异常连接；

将多组不同的所述相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇；

在多组所述异常连接簇中寻找异常种子，所述异常种子是指两个位于不同异常连接的距离不超过10bp的两个端，

当所述异常连接簇中没有异常种子，则所述基因融合的变异的类型为缺失；

当所述异常连接簇中存在一对异常种子且所述一对异常种子分别来自两条异常连接的两端，则所述基因融合的变异类型为转向；

当所述异常连接簇中存在一个异常种子且形成所述异常种子的一对异常连接的未形成种子的一端位于与所述异常种子相同的染色体上，则所述基因融合的变异的类型为插入、易位或拷贝数目变异。

6.根据权利要求5所述的检测方法，其特征在于，将多组不同的所述相关异常连接分配到不同的异常连接簇的步骤中，包括过滤属于同一种异常连接的冗余的异常连接，从而得到多组具有变异结构信息的异常连接簇。

7.根据权利要求1所述的检测方法，其特征在于，所述检测方法在序列比对的步骤之前，还包括序列连接步骤，所述序列连接步骤包括对双端测序得到的序列中具有重叠序列的片段进行连接，得到较长的连接片段。

8.根据权利要求7所述的检测方法，其特征在于，所述序列比对的步骤包括：

将所述连接片段和双端测序得到的序列中无重叠序列的非连接片段均与参考基因组序列进行比对，得到连接片段比对结果和非连接片段比对结果；

将所述连接片段比对结果和所述非连接片段比对结果进行合并，得到总比对结果；

将所述总比对结果进行排序，得到与所述参考基因组序列第一位置完全匹配的第一类序列、与所述第一位置部分匹配且与所述参考基因组上第二位置部分匹配的第二A类序列和与所述第一位置部分匹配而不与所述第二位置部分匹配的第二B类序列。

9.根据权利要求8所述的检测方法，其特征在于，所述连接片段基于单末端序列与所述参考基因组序列进行比对，所述非连接片段基于双末端序列与所述参考基因组序列进行比对。

10.一种基因融合的检测装置，所述检测装置包括比对模块和局部组装模块，其特征在于，所述检测装置还包括局部聚类模块，其中，

所述局部聚类模块：与所述比对模块相连，用于将所述比对模块输出的第二A类序列和第二B类序列基于与第一位置和第二位置的距离进行局部聚类，得到一系列第一位置聚类簇和一系列第二位置聚类簇；

所述局部组装模块：与所述局部聚类模块相连，用于将所述一系列第二位置聚类簇中的第二A类序列和第二B类序列进行局部组装，得到组装序列；并将所述组装序列再次比对到参考基因组上，查看所述组装序列的比对结果是否与所述一系列第二位置聚类簇中第二A类序列的比对结果一致，若一致，则证明所述第一位置和第二位置即为基因融合位置。

11.根据权利要求10所述的检测装置，其特征在于，所述局部聚类模块包括：

第一位置聚类簇单元：与所述比对模块相连，通过将所述比对模块输出的所述第二A类序列和所述第二B类序列中与所述第一位置偏差不超过3bp的序列进行聚类，得到基于所述第一位置的所述一系列第一位置聚类簇；

第二位置聚类簇单元：与所述比对模块相连，通过将所述比对模块输出的所述第二A类序列中与所述第二位置偏差不超过3bp的序列进行聚类，得到基于所述第二位置的一系列第二位置聚类簇，同时将所述第二B类序列中不能与所述第二位置进行匹配的部分序列分别与所述一系列第二位置聚类簇中的序列进行比对，若所述第二B类序列能与所述一系列第二位置聚类簇中的某一个第二位置聚类簇的序列相匹配，则将所述第二B类序列归于所述某一个第二位置聚类簇中。

12.根据权利要求11所述的检测装置，其特征在于，所述检测装置还包括过滤模块，所述过滤模块包括：

一次过滤单元：与所述局部聚类模块相连，通过将所述局部聚类模块输出的所述一系列第二位置聚类簇中所包含的第二A类序列和第二B类序列的数目之和小于设定值的第二位置聚类簇过滤掉，得到一次过滤后的一系列第二位置聚类簇；

二次过滤单元：与所述一次过滤单元相连，并根据所述第二A类序列中的第一位置和第二位置的比对结果所判定的所述第一位置和所述第二位置之间的连接方式，将所述一次过滤后的一系列第二位置聚类簇中与所述连接方式不一致的第二位置聚类簇过滤掉，得到二次过滤后的一系列第二位置聚类簇，并将所述二次过滤后的一系列第二位置聚类簇输入所述局部组装模块。

13.根据权利要求12所述的检测装置，其特征在于，所述检测装置还包括PE验证模块，与所述过滤模块相连，通过读取所述二次过滤单元输出的二次过滤后的一系列第二位置聚类簇中在所述第一位置和第二位置存在PE关系的序列，若根据所述存在PE关系的序列与所述参考基因组的比对结果所确定的第一位置和第二位置之间的连接方式，与所述存在PE关系的序列所在的第二位置聚类簇中的序列所确定的第一位置和第二位置之间的连接方式一致，则将所述存在PE关系的序列记为结构异常PE支持簇；

当所述结构异常PE支持簇中序列的数目达到设定阈值时，进一步验证所述二次过滤后的一系列第二位置聚类簇中的序列确实支持所述第一位置和第二位置为基因融合位置；并将所述二次过滤后的一系列第二位置聚类簇的序列输入所述局部组装模块。

14.根据权利要求13所述的检测装置，其特征在于，所述检测装置还包括基因融合的变异类型解析模块，所述基因融合的变异类型解析模块包括：

异常连接分组单元：与所述局部组装模块相连，通过将所述局部组装模块输出的所述二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，所述第一位置和第二位置分别记为所述异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端位置相距不超过10bp，则将这两个异常连接记为相关异常连接；将多组不同的所述相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇；

变异类型解析单元：与所述异常连接分组单元相连，通过从所述异常连接分组单元输出的多组具有变异结构信息的所述异常连接簇中寻找异常种子，所述异常种子是指两个位于不同异常连接的距离不超过10bp的两个端，

当所述异常连接簇中没有异常种子，则所述基因融合的变异的类型为缺失；

当所述异常连接簇中存在一对异常种子且所述一对异常种子分别来自两条异常连接的两端，则所述基因融合的变异类型为转向；

当所述异常连接簇中存在一个异常种子且形成所述异常种子的一对异常连接的各自未形成种子的一端位于相同的染色体上，则所述基因融合的变异的类型为插入、易位或拷贝数目变异。

15.根据权利要求14所述的检测装置，其特征在于，所述异常连接分组单元中还包括：

相关异常连接子单元：与所述局部组装模块相连，通过将所述局部组装模块输出的所述二次过滤后的一系列第二位置聚类簇中具有第一位置和第二位置的序列记为一个异常连接，所述第一位置和第二位置分别记为所述异常连接的两端；若一个异常连接的任意一端与另外一个异常连接的任意一端位置相距不超过10bp，则将这两个异常连接记为相关异常连接；

冗余异常连接过滤子单元，与所述相关异常连接子单元相连，通过读取所述相关异常连接子单元中属于同一种异常连接的冗余异常连接，从而得到多组不同的相关异常连接；

异常连接分组子单元：与所述冗余异常连接过滤子单元相连，通过将所述冗余异常连接过滤子单元输出的多组不同的相关异常连接分配到不同的异常连接簇中，得到多组具有变异结构信息的异常连接簇，并将所述多组具有变异结构信息的异常连接簇输入所述变异类型解析单元。

16.根据权利要求10所述的检测装置，其特征在于，所述检测装置还包括连接模块：与所述序列比对模块相连，用于将具有重叠序列的双端测序所得的序列进行连接，得到较长的连接片段，并将所述连接片段输入所述序列比对模块中。

17.根据权利要求16所述的检测装置，其特征在于，所述比对模块包括：

比对单元：通过将双端测序得到的序列中无重叠序列的非连接片段和所述连接模块得到的连接片段分别与所述参考基因组序列进行比对，得到非连接片段比对结果和连接片段比对结果；

合并单元：与所述比对单元相连，用于将所述非连接片段比对结果和所述连接片段比对结果进行合并，得到总比对结果；

排序单元：与所述合并单元相连，用于将所述总比对结果进行排序，得到与所述参考基因组序列第一位置完全匹配的第一类序列、与所述第一位置部分匹配且与所述参考基因组上第二位置部分匹配的第二A类序列和与所述第一位置部分匹配而不与所述第二位置部分匹配的第二B类序列。

18.根据权利要求17所述的检测装置，其特征在于，所述比对模块中的比对单元包括：

单末端序列比对子单元：与所述连接模块相连，用于将所述连接片段与所述参考基因组序列进行比对，得到连接片段比对结果，并将所述连接片段比对结果输入合并单元；

双末端序列比对子单元：用于将所述双端测序得到的序列中无重叠序列的非连接片段与所述参考基因组序列进行比对，得到非连接片段比对结果，并将所述非连接片段比对结果输入合并单元。