SOAPdenovo拼接参数详解

评论5,605

1. kmer值

kmer值哪个比较好很难说,对不同的数据,用不同的kmer值会有很不同的结果。最好的办法就是测试不同的kmer,然后看结果的N50,找到N50最高的kmer。不过SOAPdenovo最新的版本已经支持最长127bp的kmer了,所以要从20多测试到127,显然不太可能。下面是文档中对kmer的说明。

How to set K-mer size?
The program accepts odd numbers between 13 and 127. Larger K-mers would have higher rate of uniqueness in the genome and would make the graph simpler, but it requires deep sequencing depth and longer read length to guarantee the overlap at any genomic location.

根据我的实际使用经验,如果你的read足够长,覆盖度足够高,kmer设的越高越好。

但是实际情况是,测序的覆盖度经常不够,或者用早期的GA平台测出来read长度只有35bp,或者为了节省成本,在mate-pair library(长片段insert的文库,一般>2kb)测序时双端只有70bp,甚至40bp之类的,情况比较复杂。

一般来说,我尽量使用更高的kmer,如果我有100bp的pair-end,50bp的mate-pair,而且覆盖度挺高,我就用到kmer=45左右,如果mate-pair只有40bp,kmer=35左右。如果mate-pair更短,只有35bp,kmer值就再降一点。

但是覆盖度不够时,我一般还是使用kmer=25来拼接。

SOAP推出了一个新的工具,prepare模块,似乎就是为了解决混合长度read的问题,你可以先用很高的kmer进行contig的拼接,只使用来自180bp,300bp,500bp双端100bp的pair-end文库的reads。之后使用这些contig进行scaffolding,你可以为这些contig重新构建一个短的kmer graph,然后整合来自mate-pair文库的短read进行scaffolding.只是遗憾的是,这个模块我尝试过,不行。这个模块的运行SOAP应该要更详细说明才行。

7) Data Preparation Module generates necessary data for SOAPdenovo to run "map" and "scaff" steps from Contigs generated by SOAPdenovo or other assemblers with various length of kmer.

options:

-g [string] Prefix of output.
-K [int] Kmer length.
-c [string] Input Contigs FASTA. (Filename cannot be prefix.contig)

2. config文件中的一个重要参数reverse_seq

这个参数在很多使用soapdenovo进行拼接的人当中都会设置错误,因为默认是0,下面是从其他人博客中看到的:
“reverse_seq
This option takes value 0 or 1. It tells the assembler if the read sequences need to be complementarily reversed.
Illumima GA produces two types of paired-end libraries: a) forward-reverse, generated from fragmented DNA ends with typical insert size less than 500 bp; b) forward-forward, generated from circularizing libraries with typical insert size greater than 2 Kb. The parameter “reverse_seq” should be set to indicate this: 0, forward-reverse; 1, forward-forward.

RF: first read of fragment pair is sequenced as anti-sense (reverse), and second read is in the sense strand (forward);

FR: first read is in the sense strand (forward);second read of fragment pair is sequenced as anti-sense (reverse)”

“read/1,read/2哪个是正向?哪个是反向的?

,这个是不能确定的,在华大这边建库小插入片段(=2000bp)的文库是打断后环化,环化后再打断测,这时称为reverse_seq,在soapdenovo里面reverse_seq=1

转录组的不用管这个了,都是小于2k的”

根据我们的实际经验,如果你有mate-pair的文库,那你得咨询建库的人,中间是否经过了环化这个步骤,如果有,则必须把revser_seq设置为1。之前有个孩子不看参数,全部默认0,拼出来的N50只有30k,惨不忍睹,把这个参数纠正后,N50就有150k了。

其实这个参数就算设置正确,还是存在一个很严重的问题,因为在illumina mate-pair library建库的过程中,由于实验方法上的技术缺陷,很多dna片断并没有被成功的环化,这些没有被环化的片断测序后实际上还是pair end,也就是说两个read的方向是FR,而非RF,这时你用了reverse_seq=1,这些read的方向就是错的。所以,如果有可能,特别是你已经有reference的时候,尽量先对mate-pair文库的read进行筛选,把那些insert-size远小于理论值的,方向不正确的read删掉,否则你的拼接将会引入大量的错误。

3.其他几个对N50有重要影响的参数

-M mergeLevel(default 1,min 0, max 3): the strength of merging similar sequences during contiging
这个参数似乎在上一篇博文中没有解释,并且SOAP的提示也很简单。默认情况下M=1,M的值可以设置,从0到3。这个参数的作用是在拼contig的过程中,对buble合并和分解的一个重要参数。详细可见下面一段文字:
We used Dijkstra’s algorithm to detect bubbles, which is similar to
the ‘‘Tour-bus’’ method in Velvet.We merged the detected bubbles
into a single path if the sequences of the parallel paths were very
similar; that is, only had a single base pair difference or had fewer
than four base pairs difference with >90% identity.

简单说,就是一些差别很小的kmer,例如只有1个mismatch,在M值高的时候,将会被合并到一个node里面。
这个对于杂合度高的基因组意义比较大,因为两个allel在拼接时可能会被独立的拼出来,如果你把M值调高,这些allel就可以被合并在一起。
但是如果一个基因组repeat很多的话,M值设高了就会把很多差别很小的repeat序列合并在一起,那么相当多的序列将无法被用于构建edges,这时你的N50就会比较差,降低这个值可以显著的提高N50,但是吧,由于测序错误的存在,或者repeat特别多时候,降低M值又会导致误拼的概率大大提升,这个问题很难说得清,有些人在拼接时候尽量提高M值,把尽可能多的repeat与及可能的测序错误都合并到一个区域,这样可以保证基因组其他区域拼接的准确率,但是代价是N50的显著降低。
一般情况下,我还是使用了默认的M=1。要是碰上repeat多,杂合度又高的基因组,我只能建议远离illumina,远离SOAPdenovo,呵呵。

-R 也是跟repeat分解有关的参数,用他可以把一些由短片段重复而被合并在一起node分解开,一般情况下也可以明显的提高N50,但是一样的,repeat很多的情况下,这个参数要慎用,误拼的概率会大大提升。

config文件里的asm_flag参数
1 (reads only used for contig assembly), 2 (only used for scaffold assembly) and 3 (used for both contig and scaffold assembly).

如果你的pair-end数据不够,那就让mate-pair文库的序列也用于contig拼接,asm_flag=3。也可以比较明显的提高N50。

4.提醒

以上这些都是我的个人经验,我相信大部分情况下还是适用的,但是每个基因组的情况的差别很大,所以一定还是要多测试,才能得到比较好的拼接结果。短read拼接最大的问题就是处理repeat序列,N50提高了,某种程度上拼接准确率也要下降。因为repeat被放到哪个scaffold,似乎在SOAP中没有特别优化,只是随机扔个位置,这样就会出现不该连的scaffold被连在一起,甚至在contig内部都会出现错误,于是人为造成了很多假的structural variation。不过一般来说,insertion,deletetion之类的错误还可以忍受,但是还是有不少inversion,translocation,这种就是比较严重的错误了。

本文转自:http://blog.sina.com.cn/s/blog_5542c24a0100uaby.html

发表评论

匿名网友