基于WGS和CBC测序策略的DNA序列拼接算法研究（五）

第四章  Atlas程序的分析与改进

 

4.1    数据的准备

 

一个完整的拼接程序应该分为两个部分，一个是数据的准备，另一个是数据的拼接。

一般数据的准备分为以下几个过程

4.1.1  片段的调整（trimming reads）

 

片段之所以要调整，是因为拼接的数据源非常容易被污染。无论是全序列shotgun(WGS)还是BAC,都无法避免的引入污染。所以必须对DNA片段进行处理。

首先是去掉重复的和错误的片段，使用程序cross_match和univec.fa数据进行比对，把重复的去掉，把可能是错误和污染的片段去掉。

Cross_match是核酸和蛋白质序列快速比对和数据库搜索的工具。是由Phil Green编写，基于高效的Smith-Waterman-Gotoh算法。除了用于找到重复片段以外，cross_match还用来比对一组reads和另一组带菌序列，产生一个标记带菌版本的reads；用于比对由Phrap产生的contig序列等等。

 

UniVec是一个用来快速定义可能是带菌序列的数据库。

 

在atlas-screen-trim-file中，$CROSSMATCH_EXE $fafile $screen_file -minmatch 12 -minscore 20 -penalty -2 -screen >$fafile.screen.stdout")就完成了对read的扫描比对过程。

其次，要做的是把很短的片段去掉，因为这样可以减少扫描片段的数目，提高寻找overlapper的速度。

在程序中使用的程序是atlas-screen-window，图4.1 是一个简单的模型，有两条DNA序列，一条是原始的，另一条是经过cross_match处理的后的片段，扫描窗口的大小为4，最小品质的分界限为2。窗口首先找到第一个满足条件的一个位置，然后开始往后搜索到最后一个类似的窗口。一直搜索下来，得到的就是一个调整后的序列（trimming reads）。

在实际应用中，trimming windows大小的选择和最小品质的下限会根据基因个体的不同而有不同的选择。 RGSP（Rat Genome Sequencing Project）中选择的就是trimming windows大小为50, 品质的下限为20。调整以后，小于100 base的WGS read通过交替分析而被淘汰。

 

4.1.2    k-mer的计算

 

通过分析这些基因的频率能够提供对几因大小和覆盖序列的真实长度的估计，在拼接最初阶段可以抑制重复拼接现象的出现。

分别计算出每一个read中，连续长度k的字符串的数目和频率。把这些数据生成一张表，它最主要的目的是为在接下来的overlapper过程中使用。k串的计算还有一个目的是为了粗略的估算整个DNA片段的长度。

对于k值的大小的确定要根据具体片段的大小确定，在这个程序中，我们的k=32。也就是说我们把每个read中连续32个字符串的频率计算出来。

产生的这张k-mer表和它们的频率接下来是用于指导拼接。为了限制重叠，采用基于重复序列和避免测试所有的基因对，复杂度是 , overlapper 将只比较一个共享“rare k-mer”。并且定义一个WGS固定最小频率为R。在内存里面需要很小心维护一张频率的表，大约每一个k-mer需要10个字节。通过记录在WGS中最小的拷贝的k-mer，我们需要保持这张表很小。在RGSP，使用R=10来产生这个表来描述59,000,000。图4.2是一个k-mer分析的部分结果（mer.o.results）：

在实际应用中，数据准备阶段还需要做其他的工作。在这个程序中，有一个子程序atlas-divide-fafile,他的作用是把FASTA