由于BWA能够支持sam输出格式,且产生较少的错误匹配,因此在进行变异检测时,BWA是常用的比对工具,前两期也介绍了基于bwa比对结果,利用samtools进行变异检测。而SOAPaligner是速度上较有优势的短序列比对工具,它很好的解决了在小内存情况下将短序列比对到参考基因组上的问题,因此SOAPaligner的作用也不容小觑。以宏基因组研究为例,在进行前期的质控时,需要考虑宿主基因组污染的影响,因而SOAPaligner可作为测序read与宿主参考基因组的比对工具,达到去除宿主污染的目的。
cd soap2.21release
1.输入文件
:以人肠道宏基因组数据为例,获取测序下机数据fastq格式双端测序数据(read1.fastq,read2.fastq),和人的参考基因组序列(以hg19为例,hg19.fa)。前两期已经对fastq格式及参考基因组获取介绍过,这里就不做介绍。
2.主要流程
A. 对参考基因组序列构建索引
./2bwt-builder ~/hg19.fa
产生13个索引文件,文件前
缀为hg19.fa.index,后缀为 .amb, .ann, .bwt, .fmv,.hot, .lkt, .pac, .rev.bwt, .rev,.fmv, .rev.lkt, .rev.pac, .sa, 和.sai。
B. 如果是单端测序,则执行该条命令,即单端比对
./soap -a read.fastq -D hg19.fa.index -o read_output
C. 如果是双末端测序,则执行该条命令,即双端比对
./soap -a read1.fastq -b read2.fastq -D hg19.fa.index -o PE_output -2 SE_output
-uunmapping.fa
3.参数说明
-a,-b 表示输入的测序文件。单端测序,仅用-a;双 端测序,使用-a,-b分别表示双端测序的两个文件。
-D 表示参考基因组索引前缀名
-o 表示输出文件
-2 仅用于双端比对。表示输出比对上但不是成对reads比对上的序列的文件
-u 输出没有被比对上的reads
其他常用参数:
-m 表示最小允许插入片段长度[默认值400]
-x 表示最大允许插入片段长度[默认值600]
-r 比对到多个位置的序列的输出,1-不输出,2-任意输出一次,3-全部输出[默认值1]
-n 低质量reads阈值[默认值5]。即当read中碱基为N的数量小于阈值时,该序列在比对前被筛掉
Note: 双端比对时-2并不是可选参数,必须加上。
4.输出文件
输出文件包括:PE_output(双端比对时,2个read都比对上的信息)
SE_output(双端比对时,单个read比对上的信息)
unmapping.fa(没有比对上的read信息)
PE_output与SE_output文件结构相同,因此用SE_output中一条比对信息进行详细说明:
结果文件格式说明(从左到右):