- 确定ActiveRegion的范围
- 进行局部组装,进而得到Haplotype信息
- 等位基因的注释以及变异结果的过滤
1.确定ActiveRegion的范围
2.进行局部组装,进而得到Haplotype信息
HC的第二部就是将第一个部分得到的ActiveRegion区间里面的信息进行计算,通过第一部分我们得到了存在”疑似”变异位点的区间,但是这个判断是基于比对软件所得到的位置信息,我们的目的是减少比对软件对后面结果的干扰,所以在第二部使用局部组装将read重新进行分析,这样会提高INDEL的准确性.
首先我们将ActiveRegion中对应的reference序列作为组装的De Bruijn图的模板,方法是将reference序列按照每个位置开始向后取k个长度的序列作为一个小片段,这样我们得到了一些长度为k的序列,并且每两个相邻的序列会有k-1个长度的完全一样的片段,我们也称这样的k长度的短序列为k-mer.我们将这些k-mer作为De Bruijn图的点,如果两个k-mer之间存在k-1长度的重叠那么对应的图中他们之间存在边,也就是说如果两个位点是相邻的,那么以他们为起点得到的k-mer对应的点也是相邻的.这是我们将得到的所有的边权重赋予为0.
同时我们建立了一个哈希表用来存储每一个不同的k-mer,哈希表的key值选取为k-mer的序列,value对应为在图中表示的顶点,正如上面提到第一阶段是用reference来构建De Bruijn图,所以表中只有由参考序列组生成的k-mer.默认情况下,程序会分别用长度为10和25构建两个不同的图,并且也可以使用参数–kmerSize去指定k-mer的长度,最后的haplotype也是从所有不同k-mer长度的图中得到的.下面一步将之前得到的由reference构建的图变成read-threading graph.所谓的RT-graph就是我们将每一个read匹配到之前的模板图中所得到的的图.构建的方法如下:
我们将read序列按照与模板图相同的k-mer长度进行处理得到n条k-mer,然后将每一个k-mer在之前的hash表中进行查询,如果存在我们就记录此k-mer在图中的位置,如果不存在我们将在图中创建一个新的点代表此k-mer并且将其插入到hash表中.然后我们遍历这一个read中的相邻位点,如果两个位点他们已经在图中是相邻的了(有一条边连接这两个点)那么将这条边的权重加一,如果他们之间是不相邻的,那么我们在图中增加一条边,然后将此边的权重赋值为1,这样我们依次对每一条read的每两个相邻位点进行操作,最终得到一个带权重的read-threading graph,我们之后按照权重信息来获得可能性最高的路径.在生成read-threading graph图的过程中我们要注意对重复的k-mer的处理,首先我们得到的图应该是足够复杂的(也就是说唯一的k-mer数量最够多,多余非唯一k-mer数量的四倍)并且图中也没有环的存在,所以在基因组的某些高重复区域,会照成图中存在环结构,并且唯一的k-mer数量也会减少.
如果之前所选择的k-mer长度无法得到满足以上条件的图,程序会自动的延长k-mer的长度以减少重复性,默认的每一次将k-mer的长度增加10,直到我们获得复合条件图或者达到迭代的上限次数(HC中默认6次),如果达到上限还不能得到满足条件的图,那么整个的assembly region会被直接丢弃.
经过以上的计算我们得到了不同的k-mer对应的图,如果直接遍历图得到haplotype的话,过多的分支和测序的随机性错误会导致时间复杂度很高,所以我们要对图进行剪枝,所谓的剪枝就是去掉图中支持度比较低的分支减少计算量,过度的剪枝会损失灵敏性,所以剪枝的边界条件就成为了一个很重要的指标,我们这里以每一个边的权值作为剪枝的根据,如果一个边的权值(也就是说支持这个边的read数量)小于2的话,我们就去掉这个边.我们可以使用–minPruning来控制这个值,提高阈值的会加快计算速度和特异性,但是会降低敏感性.降低阈值则会起到相反的效果.
下面将会对所得到的的图进行微调,