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?FF scaffold_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$,...,..... !0000000000 scaffold_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