samtools及BCFtools常用参数和使用文档

1、概述

samtools的说明文档：http://samtools.sourceforge.net/samtools.shtml
samtools是一个用于操作sam和bam文件的工具合集。在前期测序数据的处理中比较常用。其包含许多命令，以下是常用命令的介绍。

2. view

view命令的主要功能：将sam文件转换成bam文件；然后对bam文件进行各种操作，比如数据的排序(其他命令功能)和提取(这些操作是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)；最后将排序或提取得到的数据输出为bam或sam格式。

bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。

View重要参数：

Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]#默认情况下不加 region，则是输出所有的 region.Options: -b     output BAM                #默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式         -h     print header for the SAM output                #默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息         -H     print header only (no alignments)         -S     input is SAM                #默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。         -u     uncompressed BAM output (force -b)                #该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。         -F INT   filtering flag, 0 for unset [0]                   Skip alignments with bits present in INT [0]                  #数字4代表该序列没有比对到参考序列上                  #数字8代表该序列的mate序列没有比对到参考序列上         -q INT   #最低比对质量 [0]         -l STR   only output reads in library STR [null]

常用命令示例

#将sam文件转换成bam文件$ samtools view -bS abc.sam > abc.bam$ samtools view -b -S abc.sam -o abc.bam#提取比对到参考序列上的比对结果$ samtools view -bF 4 abc.bam > abc.F.bam#提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可$ samtools view -bF 12 abc.bam > abc.F12.bam#提取没有比对到参考序列上的比对结果$ samtools view -bf 4 abc.bam > abc.f.bam#根据fasta文件，将 header 加入到 sam 或 bam 文件中$ samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam

3.sort

sort命令主要功能：对bam文件进行排序

Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix>  -m 参数默认下是 500,000,000 即500M（不支持K，M，G等缩写）以下数据进行处理。对于处理大数据时，如果内存够用，则设置较大值，以节约时间。-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。

4.index

注意：必须对bam文件进行默认情况下的排序后，才能进行index。否则会报错。

建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。

Usage: samtools index <in.bam> [out.index]例子：#以下两种命令结果一样$ samtools index abc.sort.bam$ samtools index abc.sort.bam abc.sort.bam.bai

5.mpileup(重要，此工具用于call snp)

mpileup用法：该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。

其用法和最简单的例子如下：

Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]$ samtools mpileup -f genome.fasta abc.bam > abc.txt$ samtools mpileup -gSDf genome.fasta abc.bam > abc.bcf$ samtools mpileup -guSDf genome.fasta abc.bam | \           bcftools view -cvNg - > abc.vcf

mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：

scaffold_1      2841    A       11      ,,,...,....     BHIGDGIJ?FFscaffold_1      2842    C       12      ,$,,...,....^I. CFGEGEGGCFF+scaffold_1      2843    G       11      ,,...,.....     FDDDDCD?DD+scaffold_1      2844    G       11      ,,...,.....     FA?AAAA<AA+scaffold_1      2845    G       11      ,,...,.....     F656666166*scaffold_1      2846    A       11      ,,...,.....     (1.1111)11*scaffold_1      2847    A       11      ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG       %.+....-..)scaffold_1      2848    N       11      agGGGgGGGGG     !!$!!!!!!!!scaffold_1      2849    A       11      c$,...,.....    !0000000000scaffold_1      2850    A       10      ,...,.....      353333333

mpileup生成的结果包含6行：参考序列名；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量。其中第5列比较复杂,做如下解释：

1 ‘.’代表与参考序列正链匹配。2 ‘,’代表与参考序列负链匹配。3 ‘ATCGN’代表在正链上的不匹配。4 ‘atcgn’代表在负链上的不匹配。5 ‘*’代表模糊碱基6 ‘^’代表匹配的碱基是一个read的开始；’^’后面紧跟的ascii码减去33代表比对质量；这两个符号修饰的是后面的碱基，其后紧跟的碱基(.,ATCGatcgNn)代表该read的第一个碱基。7 ‘$’代表一个read的结束，该符号修饰的是其前面的碱基。8 正则式’+[0-9]+[ACGTNacgtn]+’代表在该位点后插入的碱基；比如上例中在scaffold_1的2847后插入了9个长度的碱基acggtgaag。表明此处极可能是indel。9 正则式’-[0-9]+[ACGTNacgtn]+’代表在该位点后缺失的碱基；

6.bcftools

bcftools用法：用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件，后者为其二进制文件。最主要的命令是view命令来进行SNP和Indel calling。该命令的使用例子为：

$ bcftools view -cvNg abc.bcf > snp_indel.vcf

使用bcftools得到variant calling结果后。需要对结果再次进行过滤。主要依据比对结果中第8列信息。其中的 DP4 一行较为重要，它提供了4个数据：
1 比对结果和正链一致的reads数、
2 比对结果和负链一致的reads数、
3 比对结果在正链的variant上的reads数、
4 比对结果在负链的variant上的reads数。
可以设定 （value3 + value4）大于某一阈值，才算是variant。比如：

$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ && ($3+$4)>=10 && ($3+$4)/($1+$2+$3+$4)>=0.8' snp_indel.vcf > snp_indel.final.vcf