注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Bioinformatics home

 
 
 

日志

 
 

基于WGS和CBC测序策略的DNA序列拼接算法研究(六)  

2010-05-24 11:27:08|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

 

4.3.1  程序的结构

Atlas-prep-read

1create index file

atlas-createindex</perl/bin >(perl程序)

2split fasta and quality files into 20,000 reads sub files if necessary then screen and trim

       1atlas-divide-fafile </perl/bin/>(perl程序)

       2atlas-screen-trim-file</perl/bin>(perl程序)

              Across_match</local/>(可执行程序)

              Batlas-screen-window< /bin >(可执行程序)

3.  32mer count

       1)run_mer_count</perl/lib/PrepRead.pm>(PrepRead.pm包中定义的函数)

       2)find_kmer_peak</perl/lib/PrepRead.pm>(PrepRead.pm包中定义的函数)

       3)eliminate_low_freq</perl/lib/PrepRead.pm>(PrepRead.pm包中定义的函数)

atlas-asm-wgs

1Run overlapper all vs all

1)run_overlapper</perl/lib/Atlas/AsmWgs.pm>( AsmWgs.pm包中定义的函数)

A) atlas-overlapper< /bin >(可执行程序)

2. run binner

       A) atlas-binner< /bin >(可执行程序)

3.separate bins and run Phrap assembly for each bin

       1)atlas-separate-bin-assemble</perl/bin>( perl程序)

A)read_in_fonfile</perl/lib/Atlas/AsmWgs.pm>(AsmWgs.pm包中定义的函数)

B)mk_bin_dirs</perl/lib/Atlas/AsmWgs.pm>( AsmWgs.pm包中定义的函数)

(I)extract_bin_file</perl/lib/Atlas/AsmWgs.pm>(AsmWgs.pm包中定义的函数)

                            (a)   atlas-extractbins</bin>(可执行程序)

(II)run_asm_bin</perl/lib/Atlas/AsmWgs.pm>(AsmWgs.pm包中定义的函数)

                            (a)   cross_match</local/>(可执行程序)

                            (b)   atlas-screen-window< /bin >(可执行程序)

                            (c)   Phrap</local>(可执行程序)

4. combine ace and contig files and run scaffold

       1) atlas-build-scaffold-file</perl/bin>( perl程序)

              A) create a scaffold object first

              B) build the internal graph structure for the project

              C) delete some edges

              D) add 50kb stuff in

              E) fit gap

              F) delete some edges

              G) bac end pairs

              以上各个阶段都调用了/perl/lib/Atlas/Scaffold.pm中定义的一系列函数。

 

4.3.2  程序的运行流程

原始数据通过trimmingk-mer过程以后,送给overlapper进行处理,同时对原始的数据进行索引。Binner程序把overlapper的结果编成bin文件,然后使用Extract把这些bin文件中的序列片段提取出来(这里就利用了前面的所构造的索引)。把这些read送到Phrap中完成局部的拼接,形成contigs。然后把这些contigs连接成scafflod

 

4.4     关键技术研究

       在上面的分析中,已经对trimming,k-mer,的技术进行了说明,下面还将对几个个关键的算法进行详细的分析。

 

4.4.1  overlap的检测

       两个片段可以拼接的必要条件是,这两个片段之间存在一个重叠部分Overlap,并且在重叠部分有一个最小长度为min_word的精确匹配区域Match。如果Match的长度超过max_word,则认为这两个片段是几乎相同的而不适合拼接。因此,算法首先需要找出所有满足条件的片段对(称为Read Pair)。另外,如果两个片段可以进行拼接,那么只有那些首、尾相接的Read Pair才有意义,否则由于Overlap长度太长,可以认为两个片段是同一条序列。因此可以设置一个参数max_overlap,首先把每个片段划分为首、尾两部分(对于那些长度较短的片段,首、尾可以重叠),然后只需拿一个片段的首或尾与另一个片段的尾或首进行比对即可。

       它对每个base进行打分,完全匹配+1,可以替代-1,完全不匹配的-2

       例如如果一对重叠区域片段如图4.7,它的得分为:

       1+1+1+1-2+1+1+1+1+1-2-1-1+1+1+1+1+1=8

 

根据打分,可以得出这段overlap的相似度,分数越高证明相似度越大。

RGSP中,它的设置是:最大种子频率是12,带宽为3,最大的错误匹配为3%。虽然所有的WGS-WGSBAC-WGS重叠区会被计算,但是Atlas拼接需要依赖最初的BAC-WGS的重叠区。每一个overlap被保存为一个directed edge,并且指明跨度,得分,左延伸,右延伸,strand,全局的k-mer频率。如图4.8

Span是简单的两个overlap read base的数目,得分是banded 序列在跨度范围内的得分,左右延伸是两个read之间没有重叠的的区域大小,正的值说明origin read比较长,负的值说明sink read比较长。Strand”f”说明是同一个strand,为”r”说明是相对的strand重叠区。

 

4.4.2  选择WGS reads 进入BAC 文库里

在程序中,有一个binner过程,目的就是把那些拥有overlap的片段添加到一个BAC文库里面,用于后面的拼接过程。Overlap的结果将被用做选择最好的6个重复区间的WGS read6是在程序开始就被选择的,因为它在RGSC项目中是WGS覆盖的平均值的1.5倍。怎样才算得上最好的呢?这里将有一个公式来计算权值:

 

RGSP中。权值小于35的将被去掉,权值最高的6read将被保存下来。

4.4.3  Enriched BAC(eBAC)拼接的调整

我们需要调整一系列的参数,如overlapper,binner ,Phrap等等,来保证在这最好的6BACs中产生最好的序列。于是,我们将要做以下测试:

Phrap 选项设定:一般缺省值是“-shatter_greedy, 或“-trim_qual 16 -penalty -5

124-mer32-mer的选择。

2K-mer种子频率的限定。

3BAC最大overlap的限定。

使用一组参数8101214对每一个BAC read ,超过这个限定的overlap read要被忽视。

还需要选择最好的N overlapN的选择为4,6,8.

经过反复的测试和使用相关工具配合,得到最好的参数设置:

使用32-mer24-mer都是可以的。

最大的种子频率为12

最小的overlap权值为35

12为每一个BACoverlap的极限。.

 

4.4.4  rolling-phrap

       在前面的程序的分析中,已经对它的运算流程进行了分析。在这里,将从数据这个层次对它进行分析。

把每一个bactig中最左的BAC作为一个起始的source-BAC,来自source-BACread在这个bactig中和来自overlap BACread完成拼接。这里的拼接的规模很小,小于10000read,我们称这个小数目的overlapping BACsquery-BAC。完成这种规模的拼接最好使用Phrap,因为可以使用品质信息,从而充分发挥Phrap的优势。

当拼接程序开始运行后,结果contigs将做以下检查:第一,看看结果中是不是只包含source-BAC read;第二,是不是混合了source-BAC readquery-BAC

没有overlapquery-BAC我们称为Source-only contigs。这些contigs将被证明是拼接的左端,并且在rolling-phrap中没有右边的重叠区。这些contigs将不会写入”ace”文件中,并且与这些contigs关联的read将回从这组read中删除。接下来那些留下来的来自query-BACcontigs和混合的source/query-BACcontigs将变成新的source-BAC reads。通过这个bactigtill path可以决定那些BACs与这些新的source-BAC有直接的覆盖关系。这些有重叠区域的BACs将变成新的query-BACs,它们的reads将会与新的source-BAC继续拼接。

这样产生的结果又是一组contigs,一些只包括source-BAC reads,一些包括混合source/query-BACreads,还有一些只有query-BAC构成。那些只包含source-BAC readscontig将被写到一个增长中的ace文件中,那些和它相关联的read在后面的过程中不被考虑,query-BACs变成新的source-BAC。这个过程重复执行,只到这一拼接波覆盖了全部的bactig

 

4.4.5  Atlas-scaffold and split-scaffold

split-scaffold同时应用在BAC 拼接(Phrap程序之后)和bactig拼接(rolling-phrap之后)之中。使用一个与早期公布的一个贪婪算法类似的方法,来完成所有的contigs拼接和为每个bactig拼接产生scaffolds(Myers et al. 2000; Batzoglou et al. 2002).首先被识别的是连接2contigs的所有的read paris,而缝隙大小的估算则是由插入的和连接中的read区域的大小决定的。

为了排斥contigs之间零星的连接,对于缝隙大小来说的一个适当的窗口,将被用于寻找的在contigs之中最好的一组链。提供在窗口大小中缝隙大小的Paired ends被绑定,而且还决定了决定了在二个contigs 中占优势地位的连接的绑定。通过计算平均通过read-pair链,可以得到contigs 之间的缝隙大小和它的标准偏离。插入大小小于10kb,50kb250kb的文库将继续使用这个过程处理,它们相应的窗口大小为3kb10kb50kb。估计缝隙的大小小于-4标准的那些contig pairs被丢弃。

另外,在这个过程中有多重overlapcontigs连接的也同样被丢弃。每个连接的权值是依据绑定序列中read pairs的数量来计算的。如果最高权值的连接能满足它的位置约束,则被添加到一个sacffold的末尾。新形成的scaffold和其他contigs的大小将通过read pair重新计算。这个过程反复执行,只到没有链留下两个以上的read pair

错误的连接是被一个多步骤程序所改正。这个过程的本质是在scaffold中使用相同的read-pairing引擎去发现带有不协调read pairscontigs区域。有些矛盾是应归于小部分不包括在任何contig中的WGS reads。在这些情况下,整体上而言,contig不需要被修改,只需要移除没有关系的reads。所以,我们不能简单假设那些简单的分裂contigs就是一个处理矛盾的正确做法。我们用这些矛盾的read pairs 去鉴别那些可疑的区域,那些区域是用来被用来检查错误的连接。一个在contig内部的一端而pair另一端在另一个contig中的这些read paired被认为是可疑的。如果每个内部contig 端点的左边,右边或是两边都有 read-pair可以连接到另一个contig,我们将它们定义为一个检查点。基于经验的测试,我们识别到3种规则用来决定何时分裂contigs

1) 如果read-pair覆盖的深度小于0,而且靠近一个检查点,在这个区域的分离contig

2) 如果在检查点附近BAC覆盖大于0的区域大于2kb,分离这个contig

3) 如果BAC覆盖大于0的区域小于2kb,并且左、右两边都是检查点。在检查点周围分离这个区域。

通过这一轮精练的拼接,尽可能多的contigs将基于最后两个read-pair连接生成scaffold。在每一个bactig中有两种版本的scaffold:一种是分离以后的contig组成的;另一种是splittingtrimming后的contig。最终最长的scaffold将作为最后的版本。

 

 

 

 

4.5    Atlas并行化探讨

为了提高程序的运行速度和效率,使整个程序并行运算是一种切实可行的办法。

对于程序的并行化改造,可以从两个层次上进行。一个是模块间的层次;另一个是模块内部并行化改造。

首先,介绍模块间的并行性的改造。在整个程序的运行过程中,有一个结合点是值得我们关注的,就是trimmingoverlapper的过程。这个过程是把原始的数据经过trimming以后,进行k-mer计算,得到k-mer结果后,送到overlapper中进行计算。程序中虽然也把序列文件分块了,但仍然是一个串行的过程。如果把这个过程能够并行起来,即一开始就把原始的数据进行分块,把这几块数据同时进行trimmingk-mer过程,如图4.9,这样应该可以提高整个程序的运算速度。

 

 

       其次,对于某些模块的内部,仍然可以并行化改造。本程序的一个核心运算Phrap,就是一个可以并行优化的部分。现在对于phap的并行研究已经有很多成果,例如中科院计算所生物信息课题组开发的p_Phrap,就已经能够高效的完成拼接任务。

对于Phrap并行化的关键是Layout的过程。在Layout阶段,算法的任务是通过对Overlap集合 的处理来确定所有片段间的组合关系。Phrap算法是基于Hamilton路径方法的拼接算法,算法是首先对 中所有元素按照某种评估分值进行排序,然后串行计算每个元素所涉及到的两个片段间的组合关系,最后得出所有片段间的组合关系。以Phrap为例,假设有如下的OverlapO(f1, f2)O(f1, f3) O(f4, f5),其中O(f1, f2)的分值最大,O(f4, f5)的分值最小。则Phrap将首先处理O(f1, f2),然后处理O(f1, f3),最后处理O(f4, f5)。我们可以看到,O(f1, f2)O(f1, f3)通过f1相关连,因此它们的处理不能够并行进行。但是O(f4, f5)O(f1, f2)O(f1, f3)是无关的,因而可以并行处理。我们的思想是首先对Overlap集合中全体元素进行划分,然后对划分后的元素所形成的多个独立的集合分别进行处理。当然,对于overlapper过程也可以并行处理,只是在本程序中overlapper过程已经单独分离出来成为可执行文件,所以没有办法进行修改。

  评论这张
 
阅读(840)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017