专利名称：一种高杂合二倍体基因组支架序列组装策略的制作方法
技术领域：
本发明适用于生物信息学领域，具体而言，涉及一种高杂合二倍体基因组支架序列组装策略。
背景技术：
第二代测序技术高速度、低成本的推动了动植物基因组的研究进程。现今基于第二代测序技术组装的生物信息学算法是针对序列情况简单的基因组开发的，组装结果在某些生物基因组的序列复杂程度较低的情况下会取得比较好的结果，然而当基因组较为复杂(如杂合度高于O. 5%)时组装结果会很不理想，达不到后续分析所需的最低要求。组装是将测序中得到的读长短序列(Read)，通过一些信息的处理，尽可能的得到线性的基因组序列的过程，一般得到的最终结果是中间有空白未知区域(Gap)的支架序列(Scaffold)。在主流的第二代测序技术短序列组装软件中，大多采用将读长短序列(Read)分解成定长为指定值k的短序列(Kmer)，使用这些短序列的的重叠关系来构建德布鲁因图(de Brujin graph)的方法来组装重叠群序列(Contig),再加上使用插入片段长度库(PairEnd, PE) Read信息将Contig连接成中间有空白区域的支架序列(Scaffold),再根据deBrujin graph和PE信息将Gap尽量填补或缩小。对于情况比较复杂的基因组来说，由于较多杂合和重复信息导致德布鲁因图过于复杂而难以线性化，得到的Contig序列长度较短(N50低于1K)，而在比较简单的基因组中则由于杂合度和重复度较低而不会有过多难以线性化的情况出现，Contig序列长度足够长(N50 高于 3K)。对于不同复杂程度的基因组得到不同程度Contig序列的情况，现有的组装软件并没有组装Scaffold不同针对性的策略，这样得到的Scaffold序列在基因组杂合度较高(大于O. 5%)的情况下会长度短(N50低于20K)且准确性低。
发明内容
本发明的主要目的在于提供一种高杂合二倍体基因组支架序列组装策略，以解决现有技术中基因组本身的杂合度较高而使得支架组装的效果差的问题。为了实现上述目的，本发明提供了一种高杂合二倍体基因组支架序列组装策略，所述的方法包括下述步骤
将Read比对到Contig上得到所需要映射信息Arc和Link，对应于附图5中比对映射单元；
根据Contig的长度和覆盖深度(Coverage Depth)设定阈值，将短的和高覆盖深度的Contig过滤,对应于附图5中过滤单元I ;
先由Contig两两间的Arc关系,构建Contig之间的有向连接图，使用寻找泡状结构过滤单路径算法单元对图进行处理，对应于附图5中过滤单元2 ； 使用Contig之间的Link关系，构建Contig和Temp Scaffold之间的有向连接图,对图进行线性化处理，对应于附图5中线性化单元1，线性化单元2和过滤单元3 ；
在所有插入片段库都被遍历使用后，得到最终的Temp Scaffold即为最终的Scaffold ；
根据在过滤单元2过程中储存的信息,对应Scaffold中Contig路径补回被过滤的杂合单路径，对应于附图5中补回单元。本发明所述的一种高杂合二倍体基因组支架序列组装策略，根据Reads信息寻找由于杂合导致Contig和Temp Scaffold图中出现的复杂结构并对图进行化简，与现有技术相比优点在于能够在高杂合二倍体基因组的Scaffold组装中大幅度提高组装效果(直接体现在Scaffold N50等组装指标上)，最终得到符合后续分析要求的结果。
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中
图I是根据本发明所需的两类映射信息的示意 图2是本发明寻找泡状结构过滤单路径单元的示意 图3是Scaffold处理模块中大插入片段长度单元的示意 图4是本发明过滤非关键节点单元的示意 图5是本发明采用的Scaffold策略流程图。
具体实施例方式为了使本发明的目的、技术方案及优点更加清楚明白，下面结合附图和实施例对本发明进一步详细说明，应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。支架序列(Scaffold)组装所需要的包括通过线性化德布鲁因图得到的连续的包含覆盖度信息的重叠群序列(Contig)文件和各插入长度片段库的Read文件，经过一系列的单元处理过程最终得到的结构为线性化连接的Contig，Contig之间由一定长度的空白区域填充，表示为Contig之间的距离,该线性化连接的Contig结构即为Scaffold。由于插入长度片段库的插入长度会相差较大，处理的过程也会由插入长度的等级(Rank)分开，有些处理单元只要调用一次，有些处理单元需要根据插入长度库分的等级数多次调用。参考附图5的处理流程图，具体如下
第一步，将Read比对到Contig上得到所需要映射信息(Arc和Link)。下面具体说明将一个单端读长序列(SE Read)比对到Contig上，同时比对到两条Contig的终止和起始位置的Read构建了两条Contig的无缝连接关系(Arc);将插入片段成对读长序列(PERead)比对到Contig上,两端Read之间的距离由插入片段库长度和Read长度决定，实际值与给定值有一定的标准偏差，如果各有一端Read比对到Contig上的PE Read构建了两条Contig的在一定长度范围的空白区域的连接关系(Link)。比对可采用生物信息学常见的比对软件如Blat，soap等，由于真实的Read数据上会有一定的几率出现错误，所以根据实际情况可以允许一定长度的错误匹配，具体的如IOObp长度的Read可以允许5bp长度的错误匹配。比对映射的目的是使用Read信息建立起Contig之间的有向连接关系。此步为对应比对映射单元,参考附图I。第二步,根据Contig的长度和覆盖深度(Coverage Depth)分布设定阈值,将短的和过高覆盖深度的Contig过滤。具体的如将Contig长度阈值设置为lOObp，覆盖深度阈值设置为2倍的平均覆盖深度阈值，则我们遍历所有的Contig，只要满足其长度低于IOObp或覆盖深度高于2倍的平均覆盖深度这两条件之一的，就将该Contig过滤，即该Contig与其上的Arc和Link信息都不使用到Scaffold的构建中，此步对应过滤单元I。第三步,先由Contig两两间的Arc关系，使用寻找泡状结构过滤单路径算法单元进行处理。寻找类似附图2中的Bubble结构，将Bubble两边的路径分为快路径(FastPath)和慢路径(Slow Path),将Slow Path合并到Fast Path上；过程如下我们对每条路径的每个节点算出一个时间(Time)值，用所有节点的平均Time值表示路径的Time值，Time值较小的定义为Fast Path。初始的每个Contig的Time值都设置为-I,遍历过的Contig的Time值都会经过计算而改变，Time值的计算如下
ContigB的Time值与其相连接的前一个ContigA有关，同时与ContigB的长度(Length[ContigB])和覆盖深度(Coverage[ContigB])有关；
被遍历过的Contig的Time值都会由改变，因此可以在一条路径延伸到一个Time值不为-I的Contig时,表示此Contig已经被遍历过了，由此Contig回溯上游,会有两条路径回到之前同一个Contig,这样可以寻找到一个Bubble结构。进一步的，根据这两条路径中的节点平均Time值确定一条Fast Path,Fast Path中的节点平均Time值较低，并将Slow Path所有的节点标记过滤；同时将Slow Path上的节点与Fast Path上的节点对应标记储存，此步为后续补回杂合二倍体两条单倍体的其中被过滤的单路径。由于基因组上的重复序列也会在Contig图中产生Bubble的结果，但Bubble的两条路径上的节点平均覆盖深度会位于纯合深度的范围，而由杂合形成的Bubble旁边的两条路径上的节点的平均覆盖深度则会位于杂合深度的范围，因此先决条件是如果这两条路径上的平均覆盖深度在一个合适阈值内，具体如小于O. 8倍的纯合深度，我们才将两条路径合并。合并的时候选择将Slow Path合并到Fast Path上,根据Time值计算公式，即是选择其中节点平均长度(Average Length)较长和平均覆盖深度(Average Coverage Depth)较高的一条，目的是将杂合部分组装的比较完整的一条单倍体路径保留下来。对每个寻找到Bubble结构进行处理,直到遍历过所有的Contig有向连接图中的Bubble结构,此步对应过滤单元2。参考附图2,在由Arc关系构建的Contig子图中寻找到了一个从ContigA出发，分别经过 ContigB, ContigC 和 ContigB，, ContigC’ 两条路径到达 ContigD 的 Bubble 结构。这时候每条路径上的节点的Time值都已经在寻找Bubble的过程中计算得到，确定了ContigB’，ContigC’这一条路径为 Fast Path, ContigB, ContigC 这一条路径为 Slow Path。我们再计算两条路径的平均覆盖深度，即ContigB, ContigC, ContigB’和ContigC’的平均覆盖深度，得到的值低于O. 8倍的纯合深度，这时我们才将Slow Path合并到Fast Path上。将ContigB, ContigC标记过滤，同时将这个Bubble结构信息存储到ContigA和ContigD中。第四步,使用Contig之间的Link关系,构建Contig和Temp Scaffold之间的有向连接图，对图进行线性化处理。由于PE Read是分多个长度插入片段库，由此得到的Contig之间的Link关系也是有不同插入长度之分，一般来说，插入长度越长，Contig之间的距离也会越大。小插入长度片段的插入长度低于lkbp, Contig之间的Link—般比较紧凑，即连接的Contig之间的空白区域较短；大插入长度片段的插入长度一般为2kbp至40kbp,相对的Contig间的空白区域较大。对于一个插入片段库，其插入长度服从在一个给定值为峰值的泊松分布，每一对插入长度库中的PE Read的距离都为一个峰值附近的不精确估算。插入片段长度越大，PE Read可能偏离峰值的差值会越大。因此从PE Read映射得到的Link连接的Contig之间的空白区域也会有未知的偏差。为了对Contig之间的空白区域的长度有个较为精确的估计，对于小插入长度片段和大插入长度片段的处理有一定的差异。对于小插入长度片段而言，我们将Contig作为节点，Contig之间的Link作为关系构建有向图，进行线性化处理；对于大插入长度片段而言，我们将使用之前由Contig或上一个等级Temp Scaffold线性化得到的Temp Scaffold作为节点,将两个Temp Scaffold中Contig间的Link关系转化为Temp Scaffold之间的Link关系,以Temp Scaffold之间的Link关系构建有向图，进行线性化处理。这样的处理会弱化插入片段长度对组装的影响，参考附图3。进一步的，我们将来自同一或相近长度插入片段库的Link关系放在一个等级(Rank)进行处理。对于一个小插入片段库的Rank处理，首先设定阈值过滤权重小于阈值的弱Link关系，遍历由Rank内的Link关系连接的Contig构建的有向连接图，寻找到一个范围内的泡状结构的子图，内部节点之间有严重的位置重叠，先决条件也是先对内部节点的平均覆盖深度在一个合适的阈值内(具体如小于O. 8倍的纯合深度)，然后以此子图为单元，使用过滤非关键节点算法单元对子图内部的节点进行处理直到节点位置间没有严重重叠
先在一个子图中确定其中的关键节点(Key Node), Key Node是一个子图中去掉后会影响子图连通性的节点，Key Node使用Tarjan的强连通分量算法(Tarjan’ s stronglyconnected components algorithm)进行寻找；
寻找到Key Node之后，将这些Node标记排除，将子图中的其它Node做以下处理将Node长度进行排序，将最短的Node过滤后检查子图的内部Node位置的重叠程度，如果已经没有严重的重叠程度则结束，如果仍然有严重的重叠程度则重复此过程直到子图没有严重的重叠程度。参考附图4，在图中有3个关键节点，其余4个为非关键节点。我们对这4个非关键节点按照长度从小到大排序。首先将最短的非关键节点标记过滤，再检查图中的重叠程度，发现依然有严重的重叠程度；再将现在图中的最短的非关键节点标记过滤，既是刚才倒数第二长的非关键节点，发现图中已经没有严重的重叠程度，节点之间可以线性化的相连了，我们就此结束此步处理。对每个寻找到的子图进行处理,直到遍历过此Rank内的所有Link连接的Contig。由使用Arc和Link信息过滤之后的Contig中，杂合二倍体中对应的杂合Contig已经大部分被去掉一边的杂合单倍体(由Arc信息找到的做了标记)，再按照现有技术的Scaffold策略去逐步对Contig图进行处理
对于位置关系没有冲突的有上下游的多条Contig的Link关系的Contig，采用去除跨度较大Link建立过渡Link的策略。对上下游只有一个入Link和一个出Link的子图进行线性化处理，进一步过滤由Link和位置关系判断可能是重复序列构成的Contig节点和一些造成子图混乱的问题Contig。对子图进行过滤处理后，相关信息已经改变，再重新寻找子图直到所有的可能是重复序列或是问题Contig都被过滤,这时所有Contig都保持线性化的连接。将所有小插入片段库中Rank遍历之后,将所有无间断的线性化连接的Contig转换为Temp Scaffold,此步对应线性化单元I。在引入了大插入片段库的Link之后，首先将Temp Scaffold中Contig之间的Link转换为Temp Scaffold之间的Link,再同样设定设定阈值过滤权重小于阈值的弱Link关系，遍历由Rank内的Link关系连接的Temp Scaffold构建的有向连接图，寻找到一个范围内的泡状结构的子图，内部节点之间有严重的位置重叠，先决条件也是先对内部节点的平均覆盖深度在一个合适的阈值内(具体如小于O. 8倍的纯合深度)，然后以此子图为单元，同样使用过滤非关键节点算法单元对子图内部的节点进行处理直到节点位置间没有严重重叠。这里与小插入片段库的区别在于小插入片段库子图中的节点为Contig，以Contig之间的Link构建有向图，大插入片段库子图中的节点为Temp Scaffold,以Temp Scaffold之间的Link构建有向图。进一步的,我们再按照现有技术的Scaffold策略逐步对Temp Scaffold图进行处理
对上下游只有一个入Link和一个出Link的子图进行线性化处理,进一步过滤由Link和位置关系判断可能是重复序列构成的Temp Scaffold节点和一些造成子图混乱的问题Temp Scaffold。对子图进行过滤处理后,相关信息已经改变,再重新寻找子图直到所有的可能是重复序列或是问题Temp Scaffold都被过滤,这时所有Temp Scaffold都保持线性化的连接。遍历完此Rank中的所有Temp Scaffold之后,将所有无间断的线性化连接的Temp Scaffold转换为新的更大范围里的Scaffold,此步对应线性化单元2。由于对子图内部有严重重叠无法直接线性化的Contig子图和Temp Scaffold子图的数量会明显降低，子图中从矛盾的Link关系中间断开的部分会减少，相对来说Scaffold会连接的更长。直接表现在结果中就是ScaffoldN50和N90增加。第五步，在所有插入片段库都被遍历使用后，得到最终的Temp Scaffold即为最终的Scaffold ;我们再根据保存的信息将之前使用Arc信息过滤掉的杂合路径对应补回，并最终展示在结果中，此步对应补回单元。具体的,我们再最终的同一个Scaffold中搜索到一个Bubble结构对应的Bubble两端的Contig，并且两端Contig之间存在其中保留的杂合路径，这时我们将储存在两端Contig中的被过滤的杂合路径中的Contig还原，与保留的杂合路径中Contig —起并列的输出在最后的序列结果中。以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
权利要求
1.一种高杂合二倍体支架序列组装策略，其特征在于，包括如下步骤 将Reads比对到Contig上得到所需要映射信息Arc和Link ； 根据Contig的长度和覆盖深度(Coverage Depth)设定阈值，将短的和过高覆盖深度的Contig过滤； 先由Contig之间的Arc关系,构建Contig之间的有向连接图，使用寻找泡状结构过滤单路径算法单元对图进行处理； 使用Contig之间的Link关系,构建Contig和Temp Scaffold之间的有向连接图,对图进行线性化处理； 在所有插入片段库都被遍历使用后，得到最终的Temp Scaffold即为最终的Scaffold ； 根据在寻找泡状结构过滤单路径算法单元的过程中储存的信息，对应最终Scaffold中Contig路径补回被过滤的杂合单路径，并将其输出在结果序列中。
2.根据权利要求1所述的组装策略，其特征在于，将Reads比对到Contig上得到所需要映射信息Arc和Link，包括 将一个单端读长序列(SE Read)比对到Contig上，同时比对到两条Contig的终止和起始位置的Read构建了两条Contig的无缝连接关系(Arc);将插入片段成对读长序列(PERead)比对到Contig上,两端Read之间的距离由插入片段库长度和Read长度决定,实际值与给定值有一定的标准偏差，如果各有一端Read比对到Contig上的PE Read构建了两条Contig的在一定长度范围的空白区域的连接关系(Link)。
3.根据权利要求1所述的组装策略，其特征在于，根据Contig的长度和覆盖深度设定阈值，将短的和高覆盖深度的Contig过滤,其中，被过滤的Contig与其上的Arc和Link信息都不使用到Scaffold的构建中。
4.根据权利要求1所述的组装策略，其特征在于，先由Contig两两间的Arc关系，构建Contig之间的有向连接图，使用寻找泡状结构过滤单路径算法单元对图进行处理，包括 寻找Bubble结构，将Bubble两边的路径分为快路径(Fast Path)和慢路径(SlowPath),将Slow Path合并到Fast Path上的方法。
5.根据权利要求4所述的组装策略,其特征在于,寻找Bubble结构,将SlowPath合并到Fast Path上的方法包括如下步骤对每条路径的每个节点算出一个Time值,用所有节点的平均Time值表示路径的Time值，Time值较小的定义为Fast Path ;初始的每个Contig的Time值都设置为_1，遍历过的Contig的Time值都会经过计算而改变，Time值的计算如下
6.根据权利要求4或5所述的组装策略,其特征在于，所述将SlowPath合并到FastPath上的方法还包括Fast Path和Slow Path两条路径上的节点平均覆盖深度在一个合适阈值内，才将两条路径合并；合并的时候选择路径中节点平均长度较长的一条，即将杂合部分组装的比较完整的一条单倍体路径保留下来；对每个寻找到Bubble结构进行处理，直到遍历过所有的Contig有向连接图中的Bubble结构。
7.根据权利要求I所述的组装策略，其特征在于，使用Contig之间的Link关系,构建Contig和Temp Scaffold之间的有向连接图，对图进行线性化处理包括 对于小插入长度片段而言，将Contig作为节点，Contig之间的Link作为关系构建有向图，进行线性化处理；对于大插入长度片段而言，将使用由Contig或上一个等级TempScaffold线性化得到的Temp Scaffold作为节点，将两个Temp Scaffold中Contig间的Link关系转化为Temp Scaffold之间的Link关系，以Temp Scaffold之间的Link关系构建有向图，进行线性化处理；将来自同一或相近长度插入片段库的Link关系放在一个等级(Rank)进行处理。
8.根据权利要求7的所述的组装策略，其特征在于，将来自同一或相近长度插入片段库的Link关系放在一个等级(Rank)进行处理，该处理包括如下步骤 设定阈值过滤权重小于阈值的弱Link关系,遍历由Rank内的Link关系连接的Contig和Temp Scaffold构建的有向连接图，寻找到一个范围内的泡状结构的子图，内部节点之间有严重的位置重叠，然后以此子图为单元，使用过滤非关键节点算法单元对子图内部的节点进行处理直到节点位置间没有严重重叠。
9.根据权利要求8所述的组装策略,其特征在于，由Rank内的Link关系连接的Contig和Temp Scaffold构建有向连接图,其包括如下步骤 先在一个子图中确定其中的Key Node,Key Node是一个子图中去掉后会影响子图连通性的节点，Key Node使用Tarjan的强连通分量算法进行寻找；寻找到Key Node之后,将这些Node标记排除,将子图中的其它Node做以下处理J^Node长度进行排序,将最短的Node过滤后检查子图的内部Node位置的重叠程度，如果已经没有严重的重叠程度则结束，如果仍然有严重的重叠程度则重复此过程直到子图没有严重的重叠程度。
全文摘要
本发明适用于生物信息领域，提供了一种高杂合二倍体基因组支架序列组装策略。具体包括将Reads比对到Contig上得到所需要映射信息Arc和Link；根据Contig的长度和覆盖深度(CoverageDepth)设定阈值，将短的和高覆盖深度的Contig过滤；先由Contig两两间的Arc关系，构建Contig之间的有向连接图，使用寻找泡状结构过滤单路径算法单元对图进行处理；使用Contig之间的Link关系，构建Contig和TempScaffold之间的有向连接图，对图进行线性化处理；在所有插入片段库都被遍历使用后，得到最终的TempScaffold即为最终的Scaffold；根据保存的信息将之前使用Arc信息过滤掉的杂合单路径对应补回，并最终展示在结果中。通过本发明，能够在高杂合二倍体基因组的Scaffold组装中起到至关重要的作用，最终得到符合后续分析要求的结果。
文档编号G06F19/20GK102982252SQ20121051543
公开日2013年3月20日 申请日期2012年12月5日 优先权日2012年12月5日
发明者阮航, 王海龙, 朱红梅, 李瑞强 申请人:北京诺禾致源生物信息科技有限公司