0. 该软件原理

它以kerm为节点单位，利用de Bruijn图的方法实现全基因组的组装。何为de Bruijn...............

contig 的构建过程：

 （1）选取初始Kmer， 满足两个条件。

　　①给定阈值，min_read_num， 出现在几条reads上面；

　　②Kmer出现在reads的第一位，就可以开始参与拼接；如图2-4。 表中，Kmer为4， reads长12，出现在5条reads上

    

（2）选取后继Kmer

满足一下条件：

　　①后继Kmer的前k-1个碱基与当前Kmer的后k-1碱基相同；

　　②后继Kmer尽可能出现在正在参与拼接的read中，且出现位置为read的当前pos+1， 这时称Kmer出现在该read的合适位置

　　③后继Kmer要使尽可能多的reads参与拼接。

给定当前Kmer后，后继Kmer有四种可能，如图2-5。

 

 根据正在参与拼接的reads，选择ACCG作为后继Kmer， 如图2-6，此时，reads1到reads4 当前pos为2， reads5处于拼接中断状态。此外，也要让该后继Kmer处于pos1的其他reads参与拼接，即reads6，reads7，reads8。

接下来，一直重复步骤二，直到没有合适的Kmer可以选择。当出现如下情况时，则表明没有合适的Kmer可供选择：

　　①后继Kmer不出现在当前正在参与拼接的reads中；

　　②后继Kmer，对于任意条read库中的reads，该Kmer不出现在该read上pos为1的位置，contig拼接结束，若该contig大于100bp，则保留，否则舍弃。

1.使用程序及参数

    可以一步跑完，也可以分4步

   一步跑完命令： SOAPdenove all -s example.config -K 29 -D 1 -o test >>ass.log

   四步跑完脚本：

2 参数说明：

　　-s str  配置文件

　　-o str 输出文件名前缀

　　-g str 输入文件名前缀

　　-K int 输入Kmer大小，默认23

　　-p int 线程数，默认8

　　-R  利用reads先别短的重复序列，默认不进行

　　-d int 去除频数小于该值的kmer， 默认值0

　　-D 去除频数小于该值 由kmer连接诶的边，默认值为1， 即该边上每个点的频数都小于等于1时才去除

　　-M int Contig 合并相似序列的等级， 默认值为1，max为3

　　-F 利用reads对scaffold进行填补，默认不执行

　　-u 构建scaffold前不屏蔽高覆盖度的contig， 高频书覆盖度指平均contig覆盖深度的2倍。默认屏蔽

　　-G int 估计Gap大小 实际补Gap的大小差异，默认为50bp

　　-L 用于构建scaffold的contig的最短长度，默认为：Kmer参数值×2

 

3配置文件构成

其中，最重要的是配置文件的设置。其参数如下；

#maximal read length  最大长度， 应比实际reads 稍短一点

max_rd_len=100

[LIB]   每个文库信息以此为分割，每设置一个sample， 用改标记来分割。

#average insert size  平均插入片段

avg_ins=200

#if sequence needs to be reversed  反转为1， 否则为0， 一般片段>=2K的，采用环化，需要反转

reverse_seq=0

#in which part(s) the reads are used  reads用处？ 1：用于构建contig；2：构建scaffold；3：scaffold 和contig； 4;用于补洞    插入片段>=2K设为3， <2K设为3

asm_flags=3

#use only first 100 bps of each read 

rd_len_cutoff=100

#in which order the reads are used while scaffolding    reads构建scaffold的次序，值越低越优先。一般设为1， 2K设为2， 5K设为3， 10k设为4.

rank=1

# cutoff of pair number for a reliable connection (at least 3 for short insert size) 可选参数， 规定两个contig or scaffold可信连接阈值，大于该值，连接才有效。默认为3

pair_num_cutoff=3   

#minimum aligned length to contigs for a reliable read location (at least 32 for short insert size) map过程中，reads和contig 的比对长度达到该值， 该比对有效，2K设为32

map_len=32

#a pair of fastq file, read 1 file should always be followed by read 2 file     文件的路径。可以识别fasta，fastq以及bam文件。一下采用fastq格式文件输入

q1=/path/**LIBNAMEA**/fastq1_read_1.fq

q2=/path/**LIBNAMEA**/fastq1_read_2.fq

4 输出文件及说明

SOAPdenovo 分四部分别对应的输出文件：

  ● pregraph  生成7个文件 *.kmerFreq  *.edge  *.preArc  *.markOnEdge  *.path  *.vertex  *.preGraphBasic

  ● contig       生成4个文件 *.contig  *.ContigIndex  *.updated.edge  *.Arc

  ● map          生成3个文件 *.readOnContig  *.peGrads  *.readInGap

  ● scaff          生成6个文件 *.newContigIndex  *.links  *.scaf  *.scaf_gap  *.scafSeq  *.gapSeq

*.contig：contig序列文件，fasta格式；

*.scafSeq：fasta格式的scaffold序列文件，contig之间的gap用N填充；

对于得到的*.scafSeq文件还需要用GapCloser去合并其中的gap，最后的contig文件则是对补洞之后的scaffold文件通过打断N区的方法得到。

以上两个文件是组装结果中最主要的输出。

*.scaf：包括scaffold中contig的详细信息；

在scaffold行中包括scaffold名字、contig长度和该scaffold长度。

在contig行包括contig名字、contig在scaffold上的起始位置、正反链、长度和contig间的链接信息

*.links：contig间的pair-end连接信息

*.readOnContig：reads在contig上的位置

*.peGrads： 主要可以通过调整本文件中的参数来显示构建scaffold所用到的插入片段库的个数，总共要到的read数，最长的read的长度，每个库对应的哪些reads，rank设置，pair_num_cutoff设置。例如：

grads&num: 10   522083934       70

323     104577616       1       3

334     180770522       1       3

345     226070520       1       3

486     361955834       2       3

2200    392088076       3       5

2290    422272580       3       5

2400    445522690       3       5

4870    475666064       4       5

9000    511030930       5       8

9110    522083934       5       5

该文件中共分成4列。组装的配置文件中有n个文库，该文件则有n+1行，且按照文库大小顺序排列。

第1行中，第二三四列分别是 所用文库，reads总数和组装中用到的最长的reads长度。

第2行中，四列分别是文库大小，文库中的reads数目，该文库reads用到的rank等级和该文库中reads用到的pair_num_cutoff。

第3～n+1行，四列分别是文库大小，文库中的reads数目加上前面的文库中的reads总数，该文库reads用到的rank等级和该文库中reads用到的pair_num_cutoff。

如果配置文件中没有设置pair_num_cutoff，即使用默认参数，则最后一列显示为0。

对于SOAPdenovo的每个步骤都有日志文件输出，要保存好日志文件，日志文件中包含有很多有用的信息。

SOAPdenovo日志输出说明：

1）pregraph.log:  其中有很多的统计信息，包括构建debruijn-graph时用到多少reads数，构图中生成了多少uniq的kmer以及设置-d参数后去除了多少kmer。

在pregraph中，可选参数有 –R –K –d  结果如：

5467781332 nodes allocated, 70662750348 kmer in reads, 70662750348 kmer processed

3283081670 kmer removed

其中Kmer 数是取决于所设k值大小以及数据量，nodes数即特异性的kmer数目，当nodes数目过高（一般和基因组大小差不多大小），可能是数据的错误率比较高，也可能是存在杂合。若nodes数目偏小，并且kmer数目很多，则基因组本身可能存在一定的重复度。对于k值的选取，当数据量充足时（>=40X），植物基因组一般采用大kmer会有比较好的效果，而对于动物基因组，k值一般多取27和29则足够。kmer removed表示的 –d 参数所去除的低频的kmer。

2）contig.log: contig 中，可选参数 –R –D –M，注，-R 参数的选定，必须pregraph和contig中同时选择才有效。结果例子：

16430183 pairs found, 2334584 pairs of paths compared, 1674493 pairs merged

从merged的数量可以作为估计杂合以及测序错误的程度。

sum up 1932549703bp, with average length 1170

the longest is 36165bp, contig N50 is 2871 bp,contig N90 is 553 bp

相关的统计信息，一般情况下，植物基因组的contig N90都在200bp左右，若N90高于200bp，则该基因组的scaffold构建都不会有太大的问题。动物基因组，scaffold构建很少有出问题的。

3）map.log: 

Output 415219610 out of 1956217742 (21.2)% reads in gaps

1661094582 out of 1956217742 (84.9)% reads mapped to contigs

一般情况下，reads in gap的比例和map to contig 的比例总和大于1。可能是因为reads map到多个地方都被算在其中的原因。当map to contig的比例很高（80%左右时），但是组装效果并不很好，可能是重复序列比较多。reads in gap比例较高（大于40%），是因为基因组组装的较碎，gap区域较多。

map_len 默认值=K+5，当默认值大于设置的map_len时，以默认值为准，当默认值小于map_len值时，设置的map_len为准。

4）scaff.log:

average contig coverage is 23, 5832270 contig masked

构建scaffold是对高频覆盖的contig进行屏蔽（即频率高于average contig coverage的两倍的contig不用于构建scaffold），从这里可以看出组装的基因组一定的重复情况。

estimated PE size 162, by 40034765 pairs

on contigs longer than 173, 38257479 pairs found,SD=8, insert_size estimated: 163

173 是配置文件中该文库的insertsize，163 是根据reads  map到contig上的距离的估计值，8是这个分布的标准偏差。一般考虑 比对上去的pair数目和SD值。若pair对数很多且SD值很小（小片段文库数据不超过三位数，大片段文库数据部超过500），那我们一般可以将配置文件中的文库插入片段的值改对短插入片段文库（<1k）的大小估计值，一般是比较准确的，下次组装以及补洞时应根据这个值对原来配置文件中的insertsize信息做修正。对于大片段文库（>=2K），因为是把reads map到contig上，若最长contig较短时，可能找不到成pair比对上去的reads，这时，无法估计文库大小，需要自己将大片段一级一级的map到前一级的组装结果上，然后再分析大片段文库的插入片段大小。注，需要调整insertsize信息时，只需要修改* .peGrads文件中的第一列，然后删除*.links文件，重新跑scaff这一步即可。即构建scaffold时，主要是根据*.links文件的信息进行连接。

Cutoff for number of pairs to make a reliable connection: 3

1124104 weak connects removed (there were 4773564 active cnnects))

Cutoff for number是在配置文件中设的pair_num_cutoff值，weak connects是低于这个值被认定为无效的连接数，active connects是满足cutoff的连接数，根据这个数值可对pair_num_cutoff做调整

Picked  25241 subgraphs,4 have conflicting connections

conflicting connections 是表示构建scaffold时的矛盾数，矛盾数比较高（>100）时，可根据前面的有效连接数，适当提高pair_num_cutoff值，即提高scaffold连接要求的最少关系数

182483 scaffolds&singleton sum up 1990259817bp, with average length 10906

the longest is 6561520bp,scaffold N50 is 836795 bp, scaffold N90 is 157667 bp

scaffold 统计信息，将是根据rank分梯度的统计

Done with 13301 scaffolds, 2161915 gaps finished, 2527441 gaps overall

-F 参数补洞的统计信息。

5 参数调整

一般组装时需要调整的参数，主要分两种：

一种是针对脚本中的参数改动：如调整  -K  -R  -d  -D  -M

-K 值一般与基因组的特性和数据量相关，目前用到的SOAPdenovo软件主要有两个版本，grape1123和grape63mer，其中grape1123是最新版的组装软件，K值范围13-31，grape63mer是可以使用大kmer的组装版本，K值范围13-63。 

【经验】：植物基因组的组装采用大kmer效果会比较好（要求短片段reads长度75bp），动物基因组很少有用到大kmer后有明显改进效果的，且动物基因组的组装K值一般设置为27和29较多。

-R参数，对于动物基因组，R参数一般不设置，植物基因组由于较多的repeat区，则设置R参数后，效果更好。注意，设置-R时，一般使用-M 的默认值。（熊猫基因组组装时得出的结论）

-M 参数，0-3,默认值1。一般杂合率为千分之几就设为几。熊猫基因组组装时-M 2 。

-d 参数，对于没有纠错，没有处理的质量又较差的原始数据，kmer的频数为1的很多的数据的组装，一般设置为-d 1 则足够。对于处理过，或者是测序质量较好的数据，可以不用设置。数据量很多时，也可以以-d 参数去除部分质量稍差的数据。

-D 参数，默认为1，一般不用另行设置。

第二种，从map这一过程去调节参数。可以调整配置文件的map_len的值和调整文件*.peGrads。

当文库插入片段分布图中文库大小与实验给出的文库大小差异很大时，调整*.peGrads文件中的插入片段大小。

根据每一档数据的数据量去调整文库的rank等级。当该文库的数据量很多或者是在构建scaffold的过程中的冲突数很多时，可是适当的调大第四列 的pair_num_cutoff，把条件设置的更严一些。

6 内存估计

SOAPdenovo的四个步骤消耗的内存是不一样的，其中第一步消耗的内存最多，使用没有纠错的的reads，(K<=31)第一步消耗的内存在基因组大小的80－100倍左右，纠过错则在40－50倍左右，第二步相对消耗的内存会少很多，第三步消耗的内存是仅次于第一步的，在第一步的一半左右，第四步消耗的内存也会比较少。对于CPU的使用，默认是8个，如果申请内存时申请一个计算节点的所有内存测将CPU就设置为该计算节点的CPU个数充分利用计算资源，如果仅申请一个节点的部分内存则根据实际情况考虑。对于大kemr(K>31)其内存使用是(k<=31)的1.5倍左右，有时甚至更多，要充分估计内存的使用，在第一次运行的时候考虑不能太保守。

    

7 常见错误

1）配置文件中read存储路径错误

只输出日志文件。

pregraph.log中的错误信息：“Cannot open /path/**LIBNAMEA**/fastq_read_1.fq. Now exit to system...”

2）-g 后所跟参数与pregraph（第一步） -o  后所跟参数名不一致

contig  map  scaff 这三个步骤都只是输出日志文件。

contig.log中的错误信息：“Cannot open *.preGraphBasic. Now exit to system...”

map.log中的错误信息：“Cannot open *.contig. Now exit to system...”

scaff.log中的错误信息：“Cannot open *.preGraphBasic. Now exit to system...”

3）从map开始重新跑时，需要删除*.links文件，否则会生成core文件，程序退出。