发布时间 : 星期三 文章组装组参考文档version1更新完毕开始阅读6320b48f9b89680202d82527
的选取，当数据量充足时（>=40X），植物基因组一般采用大kmer会有比较好的效果，而对于动物基因组，k值一般多取27和29则足够。kmer removed表示的 –d 参数所去除的低频的kmer。
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构建很少有出问题的。
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为准。
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个数充分利用计算资源，如果仅申请一个节点的部分内存则根据实际情况考虑。对于大kmer(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文件，程序退出。
8 参考文献
参考http://soap.genomics.org.cn/soapdenovo.html相关内容。
附录五 项目数据备份和目录结构 1 项目需要备份数据
每个项目需要备份的数据，我们觉得主要就是脚本，配置文件和最优的组装结果需要保存 如目录：/nas/GAG_01D/chenyan/backup
00.data里面就是对数据处理的文件，lane.lst（主要是用杨林峰的流程处理）,lib.lst, kmer.lst（纠错时建kmer频数表的数据）和corr.lst（用于纠错的数据有哪些，主要是根据以前的经验某个短片段数据加入后组装结果反倒更差，所以没有必要再对它纠错）。 02.assembly 里面主要是最好的组装版本的组装脚本，组装的配置文件，补洞的脚本和配置文件。还有最好的组装版本的所有文件（项目一般很难得彻底结项，看项目情况吧，项目不多，基因组不大，就最好保留这些文件），统计结果，中期的报告等等，最好在里面写一个readme文件，就当是写给自己看的。
这些要备份的东西并不多，所占用空间也不大，主要是保存这些文件后，即使盘阵坏掉，直接重新跑一下，数据和组装版本就都可以重现。
还有就是每个项目需要备份的以上数据最好在当前run项目的06.backup里面备份一份，然后在自己的个人目录下（/nas/GAG_01D/）copy一份。这样应该比较安全了。
2 项目目录结构
对于项目的目录结构，还是按照以前王志文提议的那样，每个项目下分6个子目录。如：[chenyan@login-0-2 ]$ ls /nas/GAG_01A/assembly/cucumber/
data 01.analysis 02.assembly 03.soap 04.evaule 05.super-scaffold 06.backup
data主要放置原始数据链接，处理后数据，Sanger测序数据
01.analysis主要是针对初步产出的20X的分析，包括Kmer分析和细菌污染比对分析 assembly 主要是保存组装的版本
soap主要是需要分析文库的插入片段，大片段数据去除小峰，reads数据去除污染，还有GC—Depth的分布等等
evaluate 主要是对组装结果的评价，包括EST，BAC评价，全基因组线性化的比对（有reference序列） 和两个组装版本间的比较
super-scaffold 主要是针对有fosmid-end和bac-end数据的项目
backup 主要是备份中期传给客户的报告，序列，还有上述需要备份的脚本和配置文件。
3 项目存储
经过讨论，一个在线项目需要最大存储空间2T 。项目结束后需要存储的数据的存储空间100G。