基因组组装模型论文
基于reads引导的contigs构建过程如下：
1.轮拼接，初始化决策表为空，在de Bruijn图中选择拼接的初始kmer作为初始拼接的contigs，将参与拼接的reads信息添加进决策表；
2.扩展候选的k-mers，根据决策表中的reads信息，对这些kmer评分，选择得分最高的k-mer；
3.候选k-mers为空，并且已经处于第二轮拼接，则拼接停止，转到第6步；如果仍处于第一轮拼接，则标记为第二轮拼接，初始化决策表为空，将contigs反向互补，在新的3’端重新选择k-mers进行拼接；
4.选择的k-mer更新决策表中的reads的拼接信息，删除拼接成功或者失败的reads；
5.在de Bruijn图中标记拼接成功的reads的所有k-mers为删除状态，将成功的reads添加进contigs；
6.拼接结束，则保存contig及成功的reads。
该算法的拼接示意图如图4-4[5]所示。经过打分导航机制确定了一个带扩展的新kmer。该kmer在rid为3的当前read中没有出现，因此它的mismatches增1，last pos 变为0.而该kmer在rid为6的kmer中出现了，则它的matches增1，last pos 增1 ，变为14.如果此时rid为6的last pos不是13而是0，则last由0变为14.rid为2和8的read是新引进决策表的read。
图4-4contigs拼接示意（read长为36bp，k值取12）
9
基因组组装模型论文
4.2.4 contigs组装
本阶段以reads拼接阶段生成的contigs和配对数据文件作为输入，通过配对数据确定contigs之间的相对位置，连接contigs，填充相邻contigs间的gaps区域等过程，最终生成长度更长的scaffolds ，如图4-5[6]所示。
图4-5 contigs组装示意图
4.3.1模型配对文库参数的校正：
在contigs组装阶段，配对数据常用语确定相邻两个contigs的相对方向和位置关系，校正contigs中的序列重排，以及填充相邻的两个contigs间gaps区域，最终生成质量更高的超长序列scaffolds。在该阶段，reads被映射到contigs上，以获得配对reads在contigs中的位置。如果配对的两条reads分别出现在不同的contigs中，则说明这两个contigs在目标基因组中是相邻的。Contigs组装的关键是利用配对数据确定contigs的相对方向和位置关系，通常采用基于图的方法解决该问题。在图中，contigs作为顶点，contigs之间的配对连接作为边，配对俩接的数量作为边上的权重，这样就建立了任意两条contigs之间的配对关系，该问题就转化为在图中或子图中寻找遍历每个节点的最优路径。通常与给定的contigs具有配对关系的contigs数量是很少的，因此，该图并不复杂，处理起来相对简单。然而，在实际中，由于目标基因组中的重复片段以及可能测序错误的影响，contigs之间的连接有一些是不真实的（或错误的），为了保证结果的准确性，拼接算法通常忽略配对数量低于一定阀值的contigs之间的连接。
通常，配对文库的片段大小服从正态分布，并且标准差大约在均值的10%左右。在实际应用中，往往只给出配对文库的均值，而标准差是未知的。然而，配对文 库的片段尺寸的分布通常与事先给定的分布并不完全相同，总会有一些差距。因此，我们对配对数据的分布参数进行了校正。
10
基因组组装模型论文
首先，随机选取若干条contigs作为参数校正的contigs样本并对数据文件的reads进行编号。一般而言，配对数据会存储在两个文件里，对数据文件1的reads从1开始用奇数编号，对数据文件2的reads从2开始用偶数编号，这样，1、2号是配对的，3、4号是配对的，以此类推。
然后，对选取的contigs按read的长度建立索引。把两个文件里的reads往contigs上映射。通常，如果reads数据文件太大的话，可以截取文件1和文件2相对应的各一部分，再往选取的contigs上映射。
最后，分析映射结果，计算配对数据的分布参数，提取reads数据的映射信息，并按reads编号由小到大排列，这样便于分析映射到每一条contigs上的配对数据之间的相关信息。图3-5-1[7]展现了对两对配对数据映射到contigs上的情形。
图4-5-1配对数据映射到contigs
4.3.2 contigs相对位置的确定
在contigs组装之前，需要先无额定contigs之间的相对位置和方向。将配对数据往reads拼接阶段生成的contigs上映射，通常会有若干配对的reads映射到不同的contigs上，如果被配对的reads映射上的两条contigs之间有交叠，则看配对的两条reads在它们上的位置之间的距离是否在合理范围内。如果距离在合理范围内且配对reads数目达到设定的阀值，则认为两条contigs是相邻的，应该连接到一起。如果这两个contigs没有交叠，则认为它们之间存在间隙gaps，需要后续的gaps填充操作。
在contigs的构建阶段，拼接错误通常发生在contigs的末端，这是因为contigs是因找不到符合田间的下一个kmer才结束拼接的，contigs靠近末端的若干碱基正确率较低。因此，我们摒弃了contigs组装通常采用的把配对reads数据往整条contigs上映射的做法，而只取contigs的两端各Lbp长的序列片段，这样既减少了内存消耗，也有利于提高计算速度。如果contigs的长度小于L，则取其整个序列。L值的确定于配对文库的平均插入距离有关，一般取值要大于配对文库的平均插入距离，比如如果所用配对文库的平均插入距离为200bp的话，那么L就取值300。
（1）建立contigs索引；
（2）映射配对reads和contigs之间的匹配信息； （3）确定contigs的相对位置。
4.3.3 contigs连接
contigs相对位置的确定指的是contigs排放位置的确定，并没真正的连接在一起。contigs连接时需要先计算contigs之间的交叠overlap，由于contigs
11
基因组组装模型论文
之间的相对位置已经确定，我们可以根据相邻两条contigs上配对的reads来计算contigs之间的距离distance，从而确定overlap的大小。此时的distance大小可以大于0，也可以小于0。若计算出的distance值小于0，则表明这两条contigs之间可能有overlap存在；若计算出的distance值大于0，则表明这两条contigs之间可能有真正的gaps存在。
4.3.4 gaps的填充
contigs之间的相对位置确定后，在contigs连接阶段，有交叠的contigs被连接到一起，合并成长度更长的contigs，没有被连接的相邻contigs之间有空隙gaps的存在，需要进行gaps填充。gaps填充就是通过局部序列拼接的方式将contigs之间缺少的那部分碱基序列构造出来。配对数据往拼接生成的contigs上映射时，有这样的一部分配对的reads：它们只有其中的一条找到了匹配位置，而与之配对的另一条没有找到，我们称配对reads中在contigs上找到匹配位置的那些reads为悬浮hanging reads。针对每两条相邻的contigs，我们将悬浮的那些reads保存在一个表中，该表记为悬浮reads表HRL。与此同时，收集那些与hang reads对应的在contigs上未找到匹配位置的reads，并用它们构建出一个De Bruijn图。有了悬浮reads表HRL和所构建出来的De Bruijn图，我们就可以很快速的对contigs之间的gaps进行填充。如图4-5-4.1[8]所示。
图4-5-4.1gaps填充示意图
如图4-5-4.2[9]所示，选取contigs A末端READ_LEN长的序列处开始拼接，每次扩展一个碱基。待拼接kmer有效的条件：悬浮reads表中存在与kmer所在read的配对read，并且该kmer所在的read的方向与contigs的扩展方向相反以及该配对reads的距离在配对文库的误差范围之内。最优kmer的确定方法与reads拼接阶段的kmer选取方法一样，就是通过打分机制，对待选的kmer进行打分，取得分最大者。在reads拼接阶段，如果下一个kmer为空时，则停止contigs的首轮扩展，并进行contigs的反向扩展，而在填充gaps时，当下一个kmer为空时，我们把N作为扩展字符，只进行单向扩展。
12