基因组项目组装参考文档.doc

基因组项目组装参考文档基因组项目组装参考文档1一 原始数据路径2二 数据处理21 数据过滤22 数据纠错2三 根据少量数据对基因组进行分析21 基因组大小分析22 基因组杂合率模拟33 基因组是否有细菌或其他基因组污染3四 组装（SOAPdenovo）3五 基因组的soap比对分析31 文库插入片段大小的校正32 如果有其他污染，soap比对去除污染33 对于大片段文库数据，会有比较严重的小峰，需要去除小峰34 做基因组的GC-Depth分布图，查看是否有GC偏向性35 做基因组的单碱基分布图，看是否存在数据随机性不好的问题3六 基因组的组装版本的EST和BAC评价41 EST评价42 BAC评价4七 基因组的fosmid-end或bac-end数据构建super-scaffold4八 基因组数据的备份和目录结构4九 附录4附录一 filter_data_v1.241 程序简介42 基本用法：53 运行过程：54 注意事项以及常见错误：6附录二 纠错流程71 纠错策略72 纠错流程83 内存及时间损耗10附录三 Kmer分析总论101 引言 Kmerfreq程序说明102 基于kmer的基因组大小估计113 估计准确性影响因素134 结论18附录四 基因组de novo组装191 SOAPdenovo软件192 参数说明203 使用方法及示例204 输出文件及说明225 参数调整246 内存估计257 常见错误258 参考文献25附录五 项目数据备份和目录结构261 项目需要备份数据262 项目目录结构263 项目存储27一 原始数据路径/share/fqdata02solexa: 07年09年3月的数据（原raid9备份）/share/fqdata03solexa: 09年3月09年6月的数据/share/fqdata04solexa: 09年6月09年9月的数据/share/fqdata05solexa: 09年9月09年10月的数据/share/fqdata06solexa: 09年10月09年11月的数据/share/fqdata08solexa: 09年11月10年1月的数据/share/fqdata10solexa: 10年1月10年2月的数据/share/fqdata12solexa: 10年2月二 数据处理1 数据过滤filter_data流程，lane.lst 和lib.lst需要备份，生成的*.xls需要保存。文档说明见附录一。2 数据纠错石仲斌的kmerfreq，unicorn merge等程序，其中kmer.lst，kmer.sh ，corr.lst和corr.sh等文件需要备份。文档说明见附录二。三 根据少量数据对基因组进行分析1 基因组大小分析文档说明见附录三。2 基因组杂合率模拟文档说明可参考附件中的ppt3 基因组是否有细菌或其他基因组污染参考/ifs1/GAG/assemble/chenyan/Algal/01.analysis/work.sh四 组装（SOAPdenovo）文档说明见附录四五 基因组的soap比对分析1 文库插入片段大小的校正参考/nas/GAG_01C/chenyan/Penguin/03.soap/Insert/work.sh2 如果有其他污染，soap比对去除污染首先将过滤处理之后的数据与细菌基因组比对，然后去除与细菌基因组比对上的reads，再将过滤后的reads进行组装，做出GC-Depth分布图，找出其中的离群特异点，然后根据这些离群的scaffold序列去除一部分reads。3 对于大片段文库数据，会有比较严重的小峰，需要去除小峰参考/nas/GAG_01C/chenyan/Penguin/03.soap/Remove_small/work.sh4 做基因组的GC-Depth分布图，查看是否有GC偏向性/nas/GAG_01A/assembly/cucumber/03.soap/Cov2Depth/work2.sh5 做基因组的单碱基分布图，看是否存在数据随机性不好的问题参考/nas/GAG_01A/assembly/cucumber/03.soap/Cov2Depth/work.sh六 基因组的组装版本的EST和BAC评价1 EST评价参考/nas/GAG_01A/assembly/hm_Ant/04.evaule/EST/work.sh2 BAC评价参考/nas/GAG_01A/assembly/daoshuFF/04.evaule/BAC/work.sh七 基因组的fosmid-end或bac-end数据构建super-scaffold参考/nas/GAG_01A/assembly/daoshuFF/05.super-scaffold/work.sh八 基因组数据的备份和目录结构参考附录五九 附录附录一 filter_data_v1.21 程序简介数据过滤流程，主要包括以下5个步骤：filter reads with X percent base is N, set a cutoff, default 10 or polyAFilter reads with many low quality, set a cutoff , default 40Filter adapter contamination. (match length=10bp, mismatch=10bp, mismatch =10%).Filter PCR duplicate. ( read1 和read2完全一样才算是duplicate)程序实现：/nas/GAG_01A/assembly/Database/Assembly/Package/Filter_data/run_filter.pl2 基本用法：Usage:run_filter.pl maxjob, default 10（例子如下）nohup perl /nas/GAG_01A/assembly/Database/Assembly/yanglinfeng/filter_1.2_2/run_filter.pl lane.lst lib.lst 5 &保存需要处理的fq文件路径，按read1到read2的顺序。 第一行分别是fq路径，read开头需要截掉的碱基数，read末尾需要截掉的碱基数，以及对每一条reads中低质量(BCDEFG)碱基个数的cutoff(默认为40，=这个值时舍去一对reads)。第二行分别是fq路径，read开头需要截掉的碱基数，read末尾需要截掉的碱基数，以及对应的N含量的cutoff（默认为10，即一条read N含量=10%则舍去这一对）（格式如下）090709_I352_FC42CLNAAXX_L4_PHOlcpDAADEAAPE_1.fq 5 10 40090709_I352_FC42CLNAAXX_L4_PHOlcpDAADEAAPE_2.fq 5 10 20保存文库插入长度信息。分两列：文库ID以及插入长度。（格式如下）ORYcogDAIDBAAPE 170(这里当read1read2的长度+30插入长度时，Filter reads with small insert size这部将直接跳过，Small Insert %=0。这主要是针对目前测通的lane,+30是考虑到插入长度SD的问题)maxjob, default 10每次投的任务数，缺省是10。程序通过调用qsub-sge.pl 来投任务。由于程序需要读入原始fq文件以及输出处理好的结果，IO比较大，使用时注意控制任务数！3 运行过程：依次运行filter_data, duplication.pl 和stat.pl。 前两步是调用qsub-sge.pl往计算节点投任务，第三步在登陆节点进行（只是统计一下数据）。对应的分别生成*.filter.sh *.dup.sh *.stat.shfilter_data每次是读入一对reads，通过4个函数判断其是否输出，主要消耗IO，内存占用很低，所以申请的是200M，duplication.pl占用的内存比较大，短片段申请9G，大片段申请6G，stat.pl是在节点直接运行的，只是统计而以。对于101PE 6G的数据量所有的5个步骤都运行的话的运行时间在90分钟左右。输出后注意检查下结果。输出结果：在当前目录为每个库建一个子目录，最后处理好的fq文件放在相应文库的子目录下， 文件名为 “*.dup.clean”。在当前目录会生成一个 “*.stat.xls” 文件. 每一列依次是：1Library ID2Lane ID3Insert size (experiment) 4Read length5GC%6Q20%7Ns_num %8Low quality %9Adapter contamination % 10Small Insert % 11PCR duplicate % 12Raw Reads (M) 13Raw bases (Mb) 14Read length(trimmed)15Usable reads (M) 16Usable bases (Mb) 关于filter_data的代码：主函数在filter_data.cpp，调用一下三个函数实现过滤功能：1trim(string &seq1, string &seq2, string &qual1, string &qual2, int start_trim1, int end_trim1, int start_trim2, int end_trim2) 实现截取功能，具体代码见：filter_low_quality.cpp2filter_Ns(string &seq1, string &seq2, float N_rate) 过滤含N和ployA序列，具体代码见：filter_low_quality.cpp3. filter_low_qual(string &qual1, string &qual2, int Qual_rate) 过滤低质量reads，Qual_rate就是对应的cutoff，具体代码见：filter_low_quality.cpp4filter_adapter( string & seq1, string &seq2) 过滤adapter序列，具体代码见：filter_adapter.cpp5filter_small_size( string &seq1, string &seq2) 过滤small insert size序列，具体代码见：filter_small_size.cpp4 注意事项以及常见错误：lane.lst中fq文件路径是要求匹配 /Ld+_(_+)_1.fq/ ，100105_I638_FC61AWVAAXX_L2_ORYcogDAIDBAAPE_1.fq L2_和_1.fq之间的部分必须和lib.lst中的库名一致。在文件名上程序应对的是默认的命名方式，如果对输入文件有做另外的处理的话注意命名方式，建议加前缀以示区别。因为整个程序是分为三个部分的，每一步都是在前一部完整输出的前提下才能正确执行的，目前主要问题是将任务qsub到计算节点后运行可能会出现问题，这个时候qsub-sge.pl会重新投上去，第一步filter_data执行完后生成*.reads.stat以及对应的fq文件*.clean，第二步duplication.pl执行完后生成* .clean.dup.stat以及对应的fq文件* .clean.dup.clean并将第一步生成的*.clean删除，最后stat.pl生成的*.stat.xls是在fq文件的上级目录下。对于属于同一个库的不同批的数据，如果和第一批数据生成的*.stat.xls同名则在原来的*.stat.xls后加上后缀，为一个较大的数字，是对应的time生成的数值。附录二 纠错流程1 纠错策略建立高频kmer表.。遍历reads, 在一个read上找到一个连续高频kmer最多的一个区域，做为下一步动态规划的起始点S。从S向左向右遍历，如下图所示：Sreadkmer,X=A, X=Akmer,X=Akmer,X=T, X=Akmer,X=A, X=Tkmer,X=Tkmer,X=T, X=T以向左遍历为例：在S左端取n-1长度的序列N，其中n为kmer的长度。当X依次为A,C,G,T四种碱基时，考察kmer= X+N是否为高频kmer,， 若kmer是高频则存入table，初始化其上游指针为NULL，若当前X的值与对应位点read上的碱基一致，则其score=0，否则score=1。在kmer左端取n-1长度的序列N，考察kmer= X+N是否为高频 ,， 若kmer是高频则存入table，初始化其上游指针指向kmer，若当前X的值与对应位点read上的碱基一致，则其score= score，否则score= score+1。继续上述迭代，迭代结束后，遍历table找到一条分数最低的路径作为最优解。2 纠错流程程序路径：/share/data2/shizhb/workspace/MySoftware/assembly/correction/流程：kmerfreq - unicorn - merge_pair_lst.pl1. kmerfreq使用kmerfreq程序读入所有需纠错的文件（尽量挑取高质量无偏向性数据），统计kmer频数，生成频数表。Usage: kmerfreq options -i 需要纠错的reads文件列表 -o 输出的文件前缀 -q quality cutoff (默认为5,可以不变) 碱基质量值低于该值的kmer 去除 -s seed length (默认为17,需16G内存，可以不变) -n output kmer index?: 0: no, 1: yes. (默认为0,输出一个16G大小的kmer频数表，unicorn要用，如果是小数据集可以设为1，注意unicorn也有这个参数，必须一致。) -f file format: 1: fq, 2: fa. -l is there have file list?: 0: no, 1: yes. (默认为1,把需要纠错的文件全路径保存在这个file list文件中，如果只有一个文件需要纠错，可以设为0) -h -? Help脚本 /share/data2/shizhb/workspace/MySoftware/assembly/correction/kmerfreq -i kmer.lst -o ant -f 1 -l 1 -n 0 kmer.log 输出文件 ant.stat 和ant.freq2. unicorn使用unicorn程序读入需纠错的文件和频数表进行纠错。Usage: unicorn options -i 需要纠错的reads 文件list -r 指定kmerfreq 输出的kmer频数表*freq -n input kmer index?: 0: no, 1: yes. (和kmerfreq保持一致) -k start of kmer frequence cutoff (默认为5，可以不变) 纠错前，去除kmer频数低于5的kmer -e end of kmer frequence cutoff (默认为5，可以不变) 纠错后去除kmer频数低于5的kmer -d set deltas value (默认为2，可以不变) reads中有超出2个以上碱基错误的不进行纠错 -s seed length (默认为17,16G内存，和kmerfreq保持一致，可以不变) -t thread number （默认为4个） -f file format: 1: fq, 2: fa. -l is there have file list?: 0: no, 1: yes. (这个file list可以和前面一样，也可以将前面的list拆分成几份，并行的纠错。) -h -? Help脚本：/share/data2/shizhb/workspace/MySoftware/assembly/correction/unicorn -i corr.lst -r ant.freq -n 0 -f 1 -l 1 corr.log输出文件：corr.lst中的文件纠错后的*.corr文件。3. 合并使用merge_pair.pl或merge_pair_lst.pl程序将pair reads合并成一个文件。用merge_pair_lst.pl进行合并时 merge_pair_lst.pl和merge_pair.pl须放在同一目录下。脚本：perl /share/data2/shizhb/workspace/MySoftware/assembly/correction/merge_pair_lst.pl corr.lstcorr.lst文件中，纠错后的reads1和reads2文件放在一起，read1在前，read2在后。3 内存及时间损耗kmerfreq程序kmer等于17mer的时候占用内存16G。unicorn程序kmer等于17mer的时候占用内存16G。另外每个线程在处理一个文件的时候需要将该文件的所有reads读入内存。现在reads一般长度75bp， reads name 75个字节，每个文件至少25M个reads，那么一个文件要占4G左右内存。每个线程还有单独的动态规划表占1G内存。一个线程5G，unicorn程序默认开设4个线程要占36G。merge_pair_lst.pl 程序所耗内存很少kmerfreq程序统计kmer频数，输出频数表的耗时跟文件的多少和io状况有关。处理一个文件约需100s，african总共606个文件耗时15h。unicorn程序可将所有需纠错文件拆分处理，处理一个文件约需1000s，100个文件4个线程耗时1000s*100/4 = 25000s = 7h。african 606个文件拆分成6份耗时7h。African纠错总共耗时22h。参考文献：/about.html附录三 Kmer分析总论1 引言 Kmerfreq程序说明程序功能：统计测序read中kmer频数。参数介绍：-k 设定kmer长度，建议设定为奇数，默认17，占用内存16G。-a 将read从头截取长度。-d 将read从尾截取长度。这两个参数根据实际测序质量进行调整。-r 设定固定read长度，每个read仅仅截取该固定长度，用于提取kmer。建议-r和-d不要同时使用。-n 设定输出最小kmer深度。只输出深度大于该值的深度频数列表。-t 设定总碱基数上限，读取碱基数达到该值时不再读取文件。基因组大小估计：获得深度分布曲线和深度乘积曲线各自峰值深度，综合考虑估计准确峰值位置，然后根据公式：G=kmer_num/kmer_depth，估计基因组大小。通过深度1处频数估计错误率，特异kmer数参考估计基因组大小。输出图表：1numkmer_numpkdepthgenome_sizeused_baseXnode_num Xxx 306358 999579000 8 124947375 1189975000 9.5101166945其中1num为深度为1处的kmer频数，其余均由程序自身生成（*.log文件）。2 基于kmer的基因组大小估计基于短片段估计整个序列长度问题可以抽象为如下问题：假设存在完整连续序列G，随机选取片段长度为k，该片段称为kmer。当达到一定覆盖度时，根据kmer数量和深度估计序列长度G。假设：kmer深度频数分布服从泊松分布(Havlak and Chen et al., 2004)。对于泊松分布，随机变量的概率为：均值为，众数等于均值取整（floor（），因此可将峰值对应深度作为kmer期望深度。假设：选取的kmer能够遍历整个基因组时，则根据Lander_waterman算法，基因组大小（G）满足如下公式：其中，为kmer个数，为kmer期望深度，为碱基个数，为碱基期望深度，为测序生成的read个数，为测序read平均长度。因此可以获得如下公式：从上述公式可知，若获得kmer期望深度，即可计算碱基期望深度以及基因组大小。kmer深度频数分布服从泊松分布，因此可将kmer深度曲线主峰处深度作为kmer期望深度，从而估计基因组大小。3 估计准确性影响因素用拟南芥基因组，生成10X 100bp read。设定如下图内参数，其中深度有10X，20X；杂合率有0和1%；错误率有0和0.003。完成如下图：1numkmer_numpkdepthgenome_sizeused_baseXnode_num10Xh0f0306358999579000 812494737511899750009.510116694510Xh10f02616526999579000 714279700011899750008.311753275110Xh0f0341062471999579000 714279700011899750008.314260035610Xh10f0343230395999579000 616659650011899750007.1158454797结论1：杂合与错误均引起1处峰高增加，但两者的影响有量级上差异。错误的影响更大。两者也均会增加特异kmer数。结论2：10X时杂合与错误均引起主峰位置前移，两者混合可加快主峰前移。而20X时，该杂合率下并不引起主峰前移。表明小数据量时仅仅依靠错误率估计主峰前移是不够的，杂合率也会影响主峰位置的准确估计，从而导致基因组大小估计偏大。结论3： 10X时，深度乘积曲线峰值在杂合率和错误率都有的情况下也前移，但慢于深度分布曲线峰值的前移速度。由以上分析得出错误率会动过移动主峰位置严重影响基因组大小估计准确性，杂合率本身对主峰位置影响不大。但错误率、杂合率和深度三者综合作用会在不同程度上影响基因组大小估计准确性。下面分别对三者的影响进行分析。（一）测序深度的影响定义基因组大小为G，观测值为，kmer深度曲线实际峰值位置：则基因组大小偏差，由此可知，深度期望值越高，在相同情况下，基因组大小估计值偏差越小。在不考虑其他因素的情况下，取整造成的与期望深度值无关，所以高深度并不意味着取整误差一定就小。用拟南芥基因组数据（nogap），生成100bp长的read，分别生成10X，15X，30X，40X。用17kmer分析，分别获得kmer深度分布曲线和深度乘积曲线。depthKmer numKmer depthGenome sizeUsed baseXnode_numexpect_Xerror rate 109995432168124942902118993240091011676948.3996527014.9956591514993054161212494211817848874001410120891912.599399994.9952524988573682012494286829748302002310121017720.99912614.9956313029986224242511994489635697886002910121021425.198897390.795594039981656403312115653447597210003910121024133.59855011.813788结论：随着深度增加，整体估计准确性增加。但并不意味着高深度估计的一定比低深度估计的准。乘积曲线峰值位置比深度分布曲线峰值位置大1，因为深度往往估计偏小，因此推荐采用乘积曲线峰值位置也估计基因组大小，从而反映基因组大小范围。（二）错误率的影响测序错误对kmer曲线的影响体现在两个方面，一方面unique kmer（node number）增加，甚至大于基因组；另一方面，深度为1的kmer频率增加。从深度为1的kmer频率出发，对于长度为k的kmer来讲，假设1个错误碱基平均造成个特异kmer，则有如下等式：其中为深度为1的频率观测值，为不受错误影响的情况下的实际值，f为错误率，为碱基总数，为kmer总数， 为kmer深度为1的kmer个数。因为kmer深度频率服从泊松分布，则，其中即为kmer的期望峰值深度。因为往往小于100，因此随着增加，深度为1处的概率减小，当时，值为0.045%。而错误率往往导致1出频率达到40%，相差1000倍。因此可以不用考虑这个影响。从而获得等式：错误率：其中与该序列的重复特征，边际等有关，与深度也有一定关系但影响不大，拟南芥的估计值大约是0.5。用拟南芥基因组数据（nogap），生成100bp read，错误率为0，0.01，0.03，0.06和0.08的20X read。用17kmer分析，分别获得kmer深度分布曲线和深度乘积曲线。false ratio1numkmer_numpkdepthgenome_sizeused_baseused_readXnode_numFALSE0.03428281320199915800011181741636237995000023799500135881992440.0211710.02333409169199915800013153781384237995000023799500154739180870.0164810.015274594648199915800013153781384237995000023799500154039413120.0135740.01205106324199915800014142797000237995000023799500163229716110.0101390.005119303181199915800015133277200237995000023799500172265555920.0058970.00378728716199915800015133277200237995000023799500171826685690.0038920321199915800016124947375237995000023799500191012115000False ratio 是模拟数据时加入的错误率，FALSE为估计出来的错误率。结论：错误率会导致峰位置前移，使得基因组大小估计偏大，且峰值移动是离散的。这个导致的主峰位置估计不准所带来的基因组大小估计误差远大于取整造成的误差。从深度乘积曲线来看，乘积曲线受到错误率的影响较小。当错误率达到0.003时，乘积曲线保持稳定，而深度分布曲线主峰前移一位；达到0.005时，乘积曲线主峰前移一位。错误率会导致unique kmer数量增加和深度为1的kmer频数增加。可以根据1处频数，估计错误率，然后用仿真数据分析该错误率下主峰前移程度，从而修正基因组大小估计值。不同深度和错误率的综合影响：为了弄清楚深度与的关系。设定错误率为0.003，分别看10X，20X，25X，30X时获得如下统计表。1numkmer_numpkdepthgenome_sizeused_baseXnode_numFALSE10Xf034103650199957900071.43E+081.19E+098.331.43E+080.00312120Xf03787287161999158000151.33E+082.38E+0917.861.83E+080.00299425Xf03932229602398989600191.26E+082.86E+0922.621.98E+080.00295430Xf031.15E+082998737000231.3E+083.57E+0927.382.21E+080.002903其中FALSE为估计出的错误率，其中设定为0.65。当错误率1%时，设定为0.5合适。应用：根据kmer分析获得的数据和统计表格，估计错误率，在该错误率指导下仿真生成数据，分析该错误率下主峰位置的移动距离，从而调整基因组估计大小。（三）杂合率的影响用拟南芥基因组，分别加入0.001，0.003，0.005和0.01的杂合率，错误率统一设定为0.003，20X。随着杂合率的增加，在1/2主峰位置处逐渐出现杂合峰。hybridkmer_numpkdepthgenome_sizeused_baseused_readXnode_num01999158000151332772002.38E+0923799500171826685690.0011999158000151332772002.38E+092379950017.91843386650.0031999158000151332772002.38E+092379950017.91878314080.0051999158000151332772002.38E+092379950017.91911161480.011999158000151332772002.38E+092379950017.9198896818可以发现，随着杂合率的增加，特异kmer数量也会增加，和错误率相比，增加缓慢。不同深度和杂合率的综合影响：随着深度增加，杂合峰趋向于明显。应用： 估计基因组杂合率，首先需要达到一定覆盖深度，能够将杂合峰明显呈现出来。（四）不同基因组的差别用人的X染色体和酵母基因组做同上相似分析。人X染色体序列真实长度：151,058,754。Yeast序列真实长度12,070,899。仿真生成20X，长度为100bp read，设定错误率(f)为0, 0.003,0.005三个水平，杂合率(h）为0, 0.005两个水平。获得如下表格。HomoX1numkmer_numpkdepthgenome_sizeused_baseused_readXnode_numFALSEf05h01.31E+082.54E+09151691856323.02E+0930211720172.5E+080.005106f03h5866518982.54E+09151691856153.02E+0930211717172.11E+080.003374f03h0866856662.54E+09151691856323.02E+0930211720172.01E+080.003376f0h5581072.54E+09161586115303.02E+0930211720191.21E+08f0h03432.54E+09161586115303.02E+0930211720191.11E+08Yeast1numkmer_numpkdepthgenome_sizeused_baseXnode_numFALSEf05h013929396202789440151351929624141600017255557620.004849f03h09000235202789440151351929624141600017204711300.003133f03h58997648202789440151351929624141600017214315160.003132f0h5615920278944016126743402414160001912317353f0h066202789440161267434024141600019113413061num为深度为1处的kmer频数，node_num为特异kmer数，FALSE为估计出的错误率。两个序列的区别在于序列重复特征和序列长度（homoX约为yeast的10倍）。相同点：错误率均导致深度为1处kmer数增加，主峰前移，特异kmer数增加等。主峰前移速度一致。不同点： 错误率估计，homoX与拟南芥均为0.5左右，yeast为0.7。在这个参数下能够比较准确的估计出错误率。从而校正基因组大小估计值。4 结论对于真实数据，因为错误和杂合等因素会导致主峰前移，推荐使用乘积曲线峰值作为深度期望值，用深度分布曲线估计出基因组大小偏大。为了估计杂合率或者调整基因组大小估计值必须先估计错误率。估计出错误率后，调整重复参数，即可进一步估计出杂合率。附录四 基因组de novo组装基因测序仍然是基因研究的重点。其应用包括鉴定一个新物种的基因组序列，鉴定一个种群中个体的基因组序列，测定一个特定的样品的RNA分子序列，还有把DNA序列作为一个解析出来的资料应用在分子生物学技术里。尽管人类、老鼠和很多其他的基因组已经被测定出来了，鉴定一个物种的全部基因组序列依然是测序的一个重要应用。它们只是生物界中百万种生物的一个很小的部分。随着新一代高通量测序仪器的使用，测序的费用得到相当大的降低。但是，所测得到的数据读取的长度很短，这对 de novo 组装提出了非常严峻的挑战，下面将主要介绍SOAPdenovo组装软件的使用。1 SOAPdenovo软件背景SOAPdenovo是利用一种新的组装短read的方法，它能为诸如人类基因组大小(3G)的基因组de novo组装。这个程序是为组装Illumina GA short reads特别设计的。它为构建全基因组参考序列序列和以低测序成本对未知基因组实施精确分析创造了可能。程序的下载及安装： 下载地址：/soapdenovo.html安装：a 下载SOAPdenovo的压缩包 b 解压缩 c 将得到可执行文件SOAPdenovo和一个配置文件的模板example.contig使用程序及参数： 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.log2 参数说明-sSTR配置文件-oSTR输出文件的文件名前缀-gSTR输入文件的文件名前缀-KINT输入的K-mer值大小，默认值23，取值范围 13-63-pINT程序运行时设定的线程数，默认值8-R利用read鉴别短的重复序列，默认值不进行此操作-dINT去除频数不大于该值的k-mer，默认值为0-DINT 去除频数不大于该值的由k-mer连接的边，默认值为1，即该边上每个点的频数都小于等于1时才去除-MINT连接contig时合并相似序列的等级，默认值为1，最大值3。-F利用read对scaffold中的gap进行填补，默认不执行-u构建scaffold前不屏蔽高覆盖度的contig，这里高频率覆盖度指平均contig覆盖深度的2倍。默认屏蔽-G INT 估计gap的大小和实际补gap的大小的差异，默认值为50bp。-L用于构建scaffold的contig的最短长度，默认为：Kmer参数值 23 使用方法及示例1）示例SOAPdenovo_Release1.0/SOAPdenovo all -s Data/HCB.lib -K 25 -d -o test2) 输入文件configFile (配置文件内容如下，非程序生成，需要软件使用者自己配置)#maximal read length （read的最大长度）以“#”开头的行是注释内容max_rd_len=50 #该值一般设置的比实际read读长稍微短一些，截去测序最后的部分，具体长度看测序质量LIB#文库信息以此开头avg_ins=200#文库平均插入长度，一般取插入片段分布图中给出的文库大小reverse_seq=0#序列是否需要被反转，目前的测序技术，插入片段大于等于2k的采用了环化，所以对于插入长度大于等于2k文库，序列需要反转，reverse_seq1，小片段设为0asm_flags=3#该文库中的read序列在组装的哪些过程（contig/scaff/fill）中用到设为1：只用于构建contig；设为2：只用于构建scaffold；设为3：同时用于构建contig和scaffold；设为4：只用于补洞短插入片段(=2k)设为2，不用于构建contig，只用于构建scaffold，454single 长reads只用于补洞。rank=1#rank该值取整数，决定了reads用于构建scaffold的次序，值越低，数据越优先用于构建scaffold。设置了同样rank的文库数据会同时用于组装scaffold。一般将短插入片段设为1；2k设为2；5k设为3；10k设为4；当某个档的数据量较大时，也可以将其分为多个档，同样，当某档数据量不足够时，可以将多个档的数据合在一起构建scaffold。这里说的数据量够与不够是从该档的测序覆盖度和物