SOAP2是SOAP的升级版本,提高了短序列比对的运行速度和精度,同时SOAP2的一个重要改进是支持不同长度的读长。
使用步骤:
1.用2bwt-builder对fa文件建立索引
使用方法:2bwt-builder <sequence file>
2.将reads与序列进行比对
SE:/Soap/soap2.21release/soap –a <reads_a> -D -o <output>
PE:/Soap/soap2.21release/soap–a
<reads_a> -b <reads_b> -D <index.files> -o <PE_output> -2 <SE_output> -m <min_insert_size> -x <max_insert_size>
参数说明:
-D STR 基因组建索引时候的前缀名称(注意,不是基因组文件的名称)
-a STR 如果是single-end就输入single-end文件路径,paired-end则输入其中一端的reads文件路径
-b STR single-end就不用输入了,paired-end的时候输入另外一端reads文件路径
-o STR 输出比对结果文件的路径(包括文件名称)
-2 STR 输出比对上但是不是成对reads比对上的序列的文件路径(包含文件名称)。Single-end不必输入,这个文件包含的序列都是一端被比对上的序列。
-u STR 在输出文件中输出没有被比对上的reads
-m INT 最小插入片段,single-end不需要 [400]
-x INT 最大插入片段 ,single-end不需要[600]
-n INT 丢掉包含指定数量N的序列,默认5个。在质量控制的时候一般都做过N去除,这里就可以不管了
-t 用序列输出的id代替序列名称,(慎重使用)
-r INT 比对到多个位置的序列如何输出,1-不输出,2-任意输出一次,3-全部输出。全部输出的时候,同一条reads会有多个记录,之前看到有朋友说比对结果的reads统计数量变多,就是这个原因
[0,1,2] how to report repeat hits, 0=none; 1=random one; 2=all, [1]
-R 插入片段大于2K设置,低于不设置
-l INT 长reads的时候使用。从5‘端取指定长度碱基作为种子先比对,如果能比上再比对全长的reads
-v INT 一条read允许的错配数量[5]
-s INT 最小比对的长度[255]
-g INT 在一条read上准许一个连续gap大小 [5]
-M INT read或者read种植匹配的方式(不超过2个错配)
0: exact match only 只取完全匹配的
1: 1 mismatch match only 取1个错配的
2: 2 mismatch match only 取2个错配的
4: find the best hits 取最佳的
-p INT 线程, [1]
3.格式说明
I333_2_FC30JNFAAXX:7:1:3:1635/2 CTCAGCTCACCTCTCACATCTCAAGAACAGCCCATTTGATGCTG hhhH;h^hBg\PhhcQhMhhZdQAO\UZQSQL_LUQX\MW`OaP 1 b 44 +scaffold5000 16591 1 A->11T16 44M 11A32
从左到右
1. 编号: read 的编号
2. read的序列.如果read比对上参考序列的负链,会被反向互补为正链
3. 质量值:序列的质量值,和序列顺序一致,如果read反向互补,质量值也会随着改变
4. 比对上的次数: 最优比对的次数。没有比对上的read将被忽略
5. a/b:pair-end比对的标记, 表示read属于来自哪个文件
6. 长度: read长度,如果是容缺失的比对,长度将是加上缺失片断的长度
7. +/-: 比对上参考序列的正链或负链
8. 染色体名称:参考序列的染色体名称
9. 位点:第一个碱基在染色体上的位置,从1开始
10. 错配的个数
11. 错配的详细信息("A->11T16" 意思是一个错配,在参考序列的位置是第11个(从0开始),在参考序列上是A,rea上是T,质量值是16)
12. 比对上的数目 ("44M" 意思是44个碱基比对上了)
13. 对比的细节 ("11A32"意思是前1个比对上了,第11个是错配(从0开始),后面32个还是比对上了)
参考博客:http://blog.csdn.net/u014182497/article/details/51604867