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

王朝other·作者佚名 2006-01-31

宽屏版字体: 小|中|大|超大

第四章 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文件和品质文件进行分割，把大于20000个碱基以上的片段进行分割，这样才能进行trimmed过程。还需要对FASTA文件和品质创建索引文件，其目的是为了在后面的extract过程中使用。

4.2 数据的拼接

数据的拼接是程序的重点，主要有以下几个步骤

4.2.1 Overlapper

Overlapper是定义哪些WGS reads 是属于同一个基因区域的第一步处理过程。它产生了一个重叠图，记录高质量队列在共享rare序列的read。虽然overlapper对于read 排列的主要求是all-against-all read 的WGS-only拼接，但是它们的不同点是值得关注的。最值得关注的是，因为每一个BAC局部的拼接将重新测试每一个read是否有效，少量的每一个BAC的重叠需要被储存。

Overlapper的关键问题是：第一，在有限内存的运行时间问题；第二，鉴定正确的overlapper基因（算法的灵敏性）；第三，丢弃因为重复序列引起的排列（算法的独特性）；第四，详细资料的记录问题，为了在接下来帮助完成边界的化定。

这个算法成功的关键是对k-mer的保存。对WGS base 的k-mer计算能够提供对于重复序列的基因大小的视图，这些序列是能够被载入到每一个overlapper作业中。Overlapper还需要两个参数，一个是R(最小的出现次数，已经保存在k-mer频率表内部)；另一个是Y(是有可能成为overlap最大k-mer数)。每一对read的共享一个k-mer频率小于或等于的Y的将被排一列，那些小于R的都看作同一频率，绑在一起的将被强行分开。这样做可以避免对所有的read对进行比对，可以简化一个复杂度为

的工作。

这个方法虽然直接，但是在接下来的阶段中并不能很好的保证它的十分的灵敏和特殊性。在RGSP中，选取Y为8的时候，几乎有5%的基因片段有特别的k-mer数，这样影响了灵敏性；如果为了要有高的灵敏性而是R到6的话，几乎有2%的重复k-mer出现，导致大比例错误的overlaps。一个额外的banded alignment处理过程能够帮助我们采用一个大的值来定义候选overlap，但是要排除错误overlap还是需要品质序列。Overlapper的品质控制是在每一个候选overlap中计算一个端对端的banded alignment。通过限制在一个重复read和其他的片段的连续插入和删除的数目，使得banded alignment能更加直接和更加合适进行高严密的比对。

4.2.2 生成eBAC

基于每一个经过筛选的最好的overlaps，接下来的选择中把WGS read加入到每一个BAC中，包括没有通过trimming的。这些高重复的片段是可以用来拼接的，但是不能用与overlapper的分析。只有是低品质的重复片段才会被舍弃。这些BAC文件接下来将被Phrap拼接成为contigs。使用read-pair信息可以改进contig并且去掉错误，可是Phrap没有使用read-pair信息，所以需要一个工具来定义contig的错误连接并且分离它们，这就是split-scaffold。经过splitting以后，那些与scaffold有冲突的和那些缺少read pair mate的read仍然要从contig中删除。最后，split- scaffold工具创造了scaffold：因为那些contigs中每一个中都含有read-pair中的一个read。这个线形化的scaffold保存为FASTA格式，它们之间的间隙用Ns填充。

4.2.3 生成bactigs

接下来的处理中，将采用一个两步走的策略来完成eBAC和scaffold的局部拼接过程。这里称每一组overlaping克隆为一个bactig。这些bactig是使用rolling-phrap把eBAC拼接而成。为了减少错误，bactig和它的相关品质的检查必需有计划进行。它的正确性是有FPC和放射性混合STS标记图决定。潜在的错误能够在以前被确定和修正。

4.2.4 Bactig的拼接

一旦BAC布局成bactig被确定后，the reads在重叠BAC文件中用了一个叫做rolling-phrap的过程，把每个bactig聚集成一条单一的顺序。为了限制每个Phrap的范围和减小被重复的顺序引起的错误的拼接，我们把每个bactig用rolling-phrap组织和拼接成一系列的重叠的BACs，最基本的思想就是用Phrap去重新聚集一系列的重叠的BACs，为了保持合理的拼接大小，和在进程中输出的contig和 reads能沿着bactig移动。结果是个沿着bactig移动的拼接波，在一个大的拼接区域伴随着很少的错误连接。随意的选择bactig的一个终端最为左边的终端，在设计中第一波包括了所有的BACs，重叠了最左边的BACs。下一个波增加了另一个新的BACs，它覆盖了以前的波。当每个波都完成后，contig将会被删除不会与下一波重叠的BACs reads。使用这个技术，我们聚集了bactig相当于115个BACs那么大，包括了287000个reads，最后拼接的大小是12.1Mb。

如图4.3 所示，第一波包含了所以的重叠eBAC文件的最左的一个eBAC（箭头处），第二波包含所有eBAC还要增加与最左eBAC有重叠关系的所有的eBAC。这样一波一波的前进。Contig文件就保存在一个.ace文件里面（atlas.ace ），它相对应的read已经输入到Phrap而删除，直到没有或者几乎没有来自eBAC的read覆盖下一波的拼接。

4.2.5 生成superbactig

更高一级的结构称为superbactigs，它是由基于read-pair数据的bactig的连接而产生的。取代了contigs，在进程中，bactig的基本单位是scaffolds。和上面介绍的scaffolds的结构相似，两个或更多的read pairs将会连接两个scaffolds，我们能够把6247bactig scaffold基于Scaffold 融合成4584个新的scaffold ,bactig能够被连接而形成824superbactig和93 singleton，其中30个是人工合成。在superbactig中为了进一步降低scaffolds的数目，Scaffold的路径能够被用来决定scaffold的顺序和方向。Scaffold相对于clone tiling 路径的位置和方向能够通过检查BAC skin reads在scaffold中的分布情况来计算。如果他们都含有从相同BAC 克隆中继承来的BAC skim read ，Scaffold被认为是彼此临近的，.那些临近的scaffold能被连接，基于潜在的BAC细胞空隙的大小能估计出。结果，一个主要的scaffold能被superbactig获得，superbactig最后被放置在chromosome上。不能融合到主要的scaffold中的小的scaffold将不会被放在任意的chromosome文件中。在RGSP的137901个contig中，只有4289个contigs没有放置scaffolds,占有1.3%的比例。最后，每个superbactig的线性化序列都会生成。

在atlas具体实现的时候，数据的拼接部分是由5个部分构成：

1． Overlapper

使用一个过滤的算法比较所有的序列和其他的，定义出序列之间的overlapper.。少数样本的情况下使用的是一个快速和精确的比对算法。在这个程序中使用的是Hash表的方法。重复的将被选出来，结果产生一个图表。图表的每一行或者是空白的界限做为一个overlapper片段。产生的结果是直接用于指导Phrap拼接程序的运行。

2． Run binner

在进行拼接之前，为了能够是Phrap程序更好的运行，还需要做一些工作，atlas-binner就是必不可少的一步。他的作用是把每一个overlapper片段做为一个文件，产生一个文件名。核心的reads在一个二进制文件中只能出现一次，其他外围的可以重复出现。

最后的输出是一个.fon文件。如图4.4

从图上可以看出来，根据overlapper的结果，把拥有overlap的片段都组成一个bin文件，接下来这些文件被Extract分离出来。

1．把bin文件分离出来

在拼接的时候，在read片段池（reads pool）中很多文件，前面生成的一些文件例如trimming read等等，我们需要把能够用于拼接的片段选择出来。根据在数据准备阶段产生的索引文件和binner产生的bin文件清单，把相应的reads和品质文件分离出来，以便后面的拼接使用。值得注意的是在索引文件的列表上面的是真正的reads文件名，而不是索引的名字。

2．完成拼接。

这是整个程序的核心部分。在这个部分中，总共分为3个小的阶段。

首先是cross_match处理过程。

仍然是使用cross_match程序，对选择出来的，有overlapper的片段进行比对。除去重复的read.

其次是screen-windows处理过程。

把read片段小于50大小和品质比较差的删掉。

这里的数据调整的过程（trimming），是对overlapper以后的数据进行优化，为了更好的使用Phrap进行拼接。

最后是Phrap处理过程。

Phrap是用于shotgun序列拼接的程序。同样是由Phil Green编写的，是Phrap/cross_match/swat软件包的一部分。也是现在最流行拼接程序。但是Phrap在使用的时候占用内存很大，程序的容错性也不强，很容易被错误的信息所误倒，造成冗余和错误的计算。所以前面还需要对overlapper以后的数据进行比对和淘汰处理。（cross_match和screen-windows）

3． Scaffold过程。

当序列从Phrap拼接出来后，形成了contig .在实际情况中，contig和contig之间并不能直接连接起来。很多情况下是它们之间只有通过它们内部的一些小的read之间的某些距离信息或者mate信息进行连接。这个就是scaffold所做的事情。它借助其他reads之间的关系信息，把contig直接的缝隙进行填充。

在实际程序中，总共有以下几个步骤：

首先创建一个scaffold项目。

my $scaffold=Atlas::Scaffold->create();

$scaffold->set_verbose(1) if ($opt_verbose);

$scaffold->set_error_rate($opt_error_rate) if($opt_error_rate);

$scaffold->set_min_contig_length(1000);

$scaffold->read_file($file);

@contigs=$scaffold->list_contigs();

foreach my $contig (@contigs) {

$contig->find_pair_samples();

}

$scaffold->cal_contig_pair();

my $lib=$scaffold->get_library_info();

第二，取得所有contig对的尾部基因。

$scaffold->find_inter_pair(0, 15000);

接着为这个项目构造一个内部的关联图

$scaffold->build_graph();

$scaffold->set_bundle_size(3000);

$scaffold->print_edges() if($opt_verbose);

$scaffold->incremental_find_path();

然后去除一些边角余料

print "first round of addding delets back\n";

my $delete;

$delete=$scaffold->{'delete_sample'};

if($#{$delete}>=0) {

$scaffold->add_sample_to_graph($delete, 1);

$scaffold->Clean_Graph();

$scaffold->Cal_Bundle();

$scaffold->incremental_find_path(2);

}

再然后把50kb左右的原料加入进来。

$scaffold->set_bundle_size(15000);

$scaffold->Reset_Inter_Contig_Sample();

$scaffold->find_inter_pair(15000, 70000);

$scaffold->build_graph(1);

$scaffold->incremental_find_path();

$scaffold->Unify_Delete_Sample();

$delete=$scaffold->{'delete_sample'};

$scaffold->add_sample_to_graph($delete, 1);

$scaffold->Cal_Bundle();

接着找到合适的间隙，删除边角余料

print "FIT GAP........\n";

$scaffold->print_edges() if($opt_verbose);

$scaffold->gap_fill();

print "FINISH FITTING.......\n";

$scaffold->print_edges() if($opt_verbose);

$scaffold->incremental_find_path(3,1);

$scaffold->incremental_find_path(2,1);

$scaffold->Unify_Delete_Sample();

$delete=$scaffold->{'delete_sample'};

$scaffold->add_sample_to_graph($delete, 1);

$scaffold->Clean_Graph();

$scaffold->Cal_Bundle();

$scaffold->incremental_find_path(2);

再然后形成BAC文件

unless($opt_nobac) {

$scaffold->set_bundle_size(60000);

$scaffold->Reset_Inter_Contig_Sample();

$scaffold->find_inter_pair(70000);

$scaffold->build_graph(1);

if($opt_verbose);

$scaffold->incremental_find_path(2);

}

$scaffold->incremental_find_path(1) if($opt_aggressive);

$scaffold->print_rep_node();

$scaffold->print_scaffold();

close OUT;

select STDOUT;

最后，把在ace文件和关联图文件中的read产生参考序列。