SOAPdenovo

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)是打断后直接建库，测两端，没有反转；（>=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,这种就是比较严重的错误了。



SOAPdenovo虽然有文档，不过文档太简单，很多参数解释的比较模糊，经过一段时间摸索，对一些参数略有了解，要感谢张MM的仔细调试跟研究。
SOAPdenovo拼接参数相关参数分为两类，一类参数是直接可以在命令行上指定的，一类参数则需要在config文件里面进行修改。有些参数我至今还是不明白它的用处。


1.先看看命令行参数：
-o STR output graph file prefix 
-g STR input graph file prefix -K INT K-mer size [default 23] 
-p INT multithreads, n threads [default 8] 
-R use reads to solve tiny repeats [default no]
-d remove low-frequency K-mers with single occurrence [default no]
 -D remove edges comprised by entirely single frequency K-mers [default no]
 -F intra-scaffold gap closure [default no]
 -L minimum contigs length used for scaffolding 
-a initMemoryAssumption: Initiate the memory assumption to avoid further reallocation
-M INT strength of merging similar sequences during contiging [default 1, max 3]
-G gapLenDiff(default 50): allowed length difference between estimated and filled gap
 -u (optional): un-mask contigs with high coverage before scaffolding (default mask)


下面分别解释各参数：
-s STR configuration file
就是你的config文件了，config文件包含了所有的测序数据信息与及部分拼接设置
-o STR output graph file prefix 
所有输出文件的前缀名
-g STR input graph file prefix -K INT K-mer size [default 23] 
Kmer值，默认是23
-p INT multithreads, n threads [default 8] 
拼接使用的进程数，默认是8
-R use reads to solve tiny repeats [default no] 
利用read对很短的repeat进行分解，有些repeat非常短，例如十几bp得tandem repeat,很容易导致大量的kmer集中在一个节点（node)，这时候如果有read可以支持这个repeat两端的信息，就可以把这个kmer集中点分解成多个节点。
-d remove low-frequency K-mers with single occurrence [default no] 
删除只出现一次的kmer
-D remove edges comprised by entirely single frequency K-mers [default no] 
删除只由单个kmer建立的连接
这个-d 和-D我有点晕，如果使用了-d,按照字面上的解释，那么-D这个参数似乎完全没有存在的必要了，因为单个的kmer已经被删除掉了。
 -F intra-scaffold gap closure [default no] 
对scaffold内部的gap进行填充，这个参数现在似乎没什么用，因为SOAPdenovo附带了一个Gapcloser工具，就是用于scaffold内部填充的。
 -L minimum contigs length used for scaffolding 
用于建立scaffold的contig的最短长度,默认是kmer值+2，如果-g=25,则27bp以上长度的contig都会被用于scaffolding
-a initMemoryAssumption: Initiate the memory assumption to avoid further reallocation 
内存控制参数，-a 加个数值，这个参数可以限制你所需的内存。
例如你的机器有512G内存，但是有很多用户在用，这时你可以使用-a 500来占用所有的内存，防止其他用户的干扰,详细可见SOAP的文档：
"New parameter added in "pregraph" module. This parameter initiate the memory assumption to avoid further reallocation. Unit of the parameter is GB. Without further reallocation, SOAPdenovo runs faster and provide the potential to eat up all the memory of the machine. For example, if the workstation provides 50g free memory, use -a 50 in pregraph step, then a static amount of 50g memory would be allocated before processing reads. This can also avoid being interrupted by other users sharing the same machine."
-G gapLenDiff(default 50): allowed length difference between estimated and filled gap 
这个参数也有点含糊，可能是在填充scaffold的过程中，允许之前估算出来的gap跟填充进序列之后，该gap的总长度（包括填充序列）之间有50bp的差异？
 -u (optional): un-mask contigs with high coverage before scaffolding (default mask) 
在scaffolding过程中是否要将高覆盖度的区域（很可能是repeat)屏蔽掉，默认情况下是mask，也就是不要用这个参数。


2. config文件的参数
下面是SOAP提供的一个配置文件样本，我想英文普通人都看得懂，就不一一解释了
APPENDIX: example.config
#maximal read length 
max_rd_len=50 [LIB] 
#average insert size 
avg_ins=200 
#if sequence needs to be 
reversed reverse_seq=0 
#in which part(s) the reads are used 
asm_flags=3 
#in which order the reads are used while scaffolding 
rank=1 
#fastq file for read 1 
q1=/pathfastq_read_1.fq 
#fastq file for read 2 always follows fastq file for read 1 
q2=/pathfastq_read_2.fq 
#fasta file for read 1 
f1=/pathfasta_read_1.fa 
#fastq file for read 2 always follows fastq file for read 1 
f2=/pathfasta_read_2.fa 
#fastq file for single reads 
q=/pathfastq_read_single.fq 
#fasta file for single reads 
f=/pathfasta_read_single.fa 
#a single fasta file for paired reads 
p=/pathpairs_in_one_file.fa 
[LIB] avg_ins=2000
 reverse_seq=1 
asm_flags=2
 rank=2 
q1=/pathfastq_read_1.fq 
q2=/pathfastq_read_2.fq 
q=/pathfastq_read_single.fq 
f=/pathfasta_read_single.fa


最后，附上文档里对文库rank设置的说法。
2. How to set library rank?
SOAPdenovo will use the pair-end libraries with insert size from smaller to larger to construct scaffolds. Libraries with the same rank would be used at the same time. For example, in a dataset of a human genome, we set five ranks for five libraries with insert size 200-bp, 500-bp, 2-Kb, 5-Kb and 10-Kb, separately. It is desired that the pairs in each rank provide adequate physical coverage of the genome.



1 使用程序及参数：
SOAPdenovo可以一步跑完，也可以分成四步单独跑
一步跑完的脚本:
./ SOAPdenovo all -s lib.cfg -K 29 -D 1 -o ant >>ass.log
四步单独跑的脚本:
./ SOAPdenovo pregraph-s lib.cfg -d 1  -K 29 -o ant>pregraph.log
./ SOAPdenovo contig -g ant -D 1 -M 3>contig.log
./ SOAPdenovo map -s lib23.cfg -g ant>map.log
./ SOAPdenovo scaff -g ant -F>scaff.log
2 参数说明
    -s  STR     配置文件
    -o  STR     输出文件的文件名前缀
    -g  STR     输入文件的文件名前缀
    -K  INT     输入的K-mer值大小，默认值23，取值范围 13-63
    -p  INT     程序运行时设定的线程数，默认值8
    -R          利用read鉴别短的重复序列，默认值不进行此操作
    -d  INT     去除频数不大于该值的k-mer，默认值为0
    -D  INT     去除频数不大于该值的由k-mer连接的边，默认值为1，即该边上每个点的频数都小于等于1时才去除
    -M  INT     连接contig时合并相似序列的等级，默认值为1，最大值3。
    -F          利用read对scaffold中的gap进行填补，默认不执行
    -u          构建scaffold前不屏蔽高覆盖度的contig，这里高频率覆盖度指平均contig覆盖深度的2倍。默认屏蔽
    -G   INT    估计gap的大小和实际补gap的大小的差异，默认值为50bp。
    -L          用于构建scaffold的contig的最短长度，默认为：Kmer参数值×2
3 使用方法及示例
1）示例
SOAPdenovo_Release1.0/SOAPdenovo all -sData/HCB.lib -K 25 -d  -o test
2) 输入文件
configFile (配置文件内容如下，非程序生成，需要软件使用者自己配置)
#maximal read length （read的最大长度）
以“#”开头的行是注释内容
max_rd_len=50   #该值一般设置的比实际read读长稍微短一些，截去测序最后的部分，具体长度看测序质量
[LIB]           #文库信息以此开头
avg_ins=200     #文库平均插入长度，一般取插入片段分布图中给出的文库大小
reverse_seq=0   #序列是否需要被反转，目前的测序技术，插入片段大于等于2k的采用了环化，所以对于插入长度大于等    于2k文库，序列需要反转，reverse_seq＝1，小片段设为0
asm_flags=3     #该文库中的read序列在组装的哪些过程（contig/scaff/fill）中用到
                设为1：只用于构建contig；
                设为2：只用于构建scaffold；
                设为3：同时用于构建contig和scaffold；
                设为4：只用于补洞
短插入片段(<2K)的设为3，同时用于构建contig和scaffold，长插入片段(>=2k)设为2，不用于构建contig，只用于构建scaffold，454single 长reads只用于补洞。
rank=1           #rank该 值取整数，决定了reads用于构建scaffold的次序，值越低，数据越优先用于构建scaffold。设置了同样rank的文库数据会同时用于组装 scaffold。一般将短插入片段设为1；2k设为2；5k设为3；10k设为4；当某个档的数据量较大时，也可以将其分为多个档，同样，当某档数据量 不足够时，可以将多个档的数据合在一起构建scaffold。这里说的数据量够与不够是从该档的测序覆盖度和物理覆盖度两个方面来考虑的。
pair_num_cutoff=3 #可选参数，pair_num_cutoff该参数规定了连接两个contig或者是pre-scaffold 的可信连接的阈值，即，当连接数大于该值，连接才算有效。短插入片段(<2k)默认值为3，长插入长度序列默认值为5
map_len=32        #map_len该参数规定了在map过程中 reads和contig的比对长度必须达到该值（比对不容mismacth和gap），该比对才能作为一个可信的比对。可选参数，短插入片段(<2k)一般设置为32，长插入片段设置为35，默认值是K＋2。


q1=/path/**LIBNAMEA**/fastq_read_1.fq
#read 1的fastq格式的序列文件，“/path/**LIBNAMEA**/fastq_read_1.fq”为read的存储路径
q2=/path/**LIBNAMEA**/fastq_read_2.fq
#read 2的fastq格式的序列文件，与read1对应的read2文件紧接在read1之后）
f1=/path/**LIBNAMEA**/fasta_read_1.fa
#read 1的fasta格式的序列文件
f2=/path/**LIBNAMEA**/fasta_read_2.fa
#read 2的fasta格式的序列文件
q=/path/**LIBNAMEA**/fastq_read_single.fq
#单向测序得到的fastq格式的序列文件
f=/path/**LIBNAMEA**/fasta_read_single.fa
#单向测序得到的fasta格式的序列文件
p=/path/**LIBNAMEA**/pairs_in_one_file.fa
#双向测序得到的一个fasta格式的序列文件