还在为微生物重测序变异检测发愁？samtools帮助你！

Samtools作为一款操作序列比对结果文件（SAM/BAM）的工具，能够灵活转换sam/bam，并且能基于参考序列和Sam/bam文件进行变异位点检测，并通过bcftools进行变异结果统计。本期将基于上期bwa软件的比对结果，利用samtools进行变异检测。

一、软件安装

下载地址http://www.htslib.org/download/，利用wget下载得到samtools-1.5.tar.bz2，bcftools-1.5.tar.bz2；
解压缩：tar -jxvf bcftools-1.5.tar.bz2

tar -jxvf samtools-1.5.tar.bz2

安装（详细的安装信息可查看bcftools-1.5(samtools-1.5)文件夹下的INSTALL）：

cd bcftools-1.5 ( samtools-1.5)

./configure –prefix=/where/to/install#设置安装路径

Make

Make install

安装完成后，我们后面所用到的可执行程序都在上面设置的安装路径（/where/to/install）的bin文件夹下；
配置环境变量

a.若要临时修改环境变量，可直接在终端输入下面一行命令：

Export PATH=/where/to/install/bin:$PATH

b.要永久修改环境变量可将下面第一行添加到~/.bash_profile(针对当前用户)或者/etc/profile(针对所有用户)文件的末尾，再执行第二行命令即可：

Export PATH=$PATH: /where/to/install/bin；

source ~.bash_profile 或者source /etc/profile

现在samtools 和bcftools已经安装好了，下面我们就进行calling SNPs和Indels了。

二、使用流程

1. 输入文件

比对结果文件sample.sort..bam和参考基因组序列ref.fa（上一期已对其建立索引，此处不再建立索引）;

2. Variant Calling主要流程

Samtools index sample.sort.bam#对bam文件建立索引,默认生成文件为bam文件加.bai,此处生成sample.sort.bam.bai;

Samtools rmdup sample.sort.bam sample.dedup.bam#去除PCR等实验过程中产生的多余duplications;

Samtools mpileup –q 20 –d 8000 –ugo sample.bcf –f ref.fa sample.dedup.bam#variant calls

参数说明：

-q, --min-MQ mapQ最小值，低于这个值则跳过；

-d, --max-depth 最大覆盖深度；

-g, --BCF 生成bcf格式文件；

-o, --output FILE 输出文件；

-f , --fasta-ref FILE参考序列(fasta format)的索引文件名；

其他常用参数：

-A --count-orphans 不丢弃异常配对reads；

-l --positions FILE 包含区域位点的位置列表文件或者bed区域文件；

-r --region REG 在指定区域产生pipeup，和-l一起使用；

-I --skip-indel 不检测INDEL;

-m --min-ireads 候选INDEL的最小间隔的reads数；

-v --VCF 生成VCF格式文件；

-u --uncompressed 产生未压缩的vcf/bcf文件；

Bcftools call –vmO z –o sample.vcf.gz sample.bcf #将bcf转化为vcf文件，bcf为vcf的二进制格式

参数说明：

-v --variants-only 仅输出变异位点，不加此参数，则输入所有位点信息；

-c --consensus-caller the original calling method;

-m -- multiallelic-caller alternative model for multiallelic and rare-variant calling，calling的方法，与-c不能同时使用；

-O --output-type <b|u|z|v> 指定输出文件类型， 'b'表示压缩后BCF文件; 'u'表示未压缩的BCF文件; 'z'表示压缩的VCF文件; 'v' 表示未压缩的VCF文件；

三、vcf文件的过滤统计

bcftools filter -O z -o sample.filtered.vcf.gz -i 'QUAL>20 && DP>5' sample.vcf.gz#过滤质量值低于20并且深度小于5的变异位点

参数说明：

-O --output-type <b|u|z|v> 指定输出文件类型， 'b'表示压缩后BCF文件; 'u'表示未压缩的BCF文件; 'z'表示压缩的VCF文件; 'v' 表示未压缩的VCF文件；

-o, --output FILE 输出文件；

-i --Include <exp> 筛选出满足exp条件的变异位点；

-e --exclude<exp> 剔除满足exp条件的变异位点，与-i相反；

--threads 线程数；

Note: 上面的过滤条件中可以加入DP4参数，过滤出高质量的变异位点

Vcf 文件格式：由#开头的注释部分（图1）和变异主体部分（图2，图3），注释部分主要包括软件版本，命令行，参考序列信息以及主体部分参数的注释，主体部分由以tab隔开的8列组成，依次为参考序列名，variant的位置，variant的ID，参考序列的碱基，Phred格式(Phred_scaled)的质量值，过滤结果和第8列variant 的详细信息（图3），采用tag=value形式，tag间用“;”隔开；

图1 VCF文件注释部分

图2 VCF文件变异主体部分前7列

图3 VCF文件变异主体部分详细信息和基因型列

bcftools stats sample.filtered.vcf.gz > sample.stat.file# 统计SNPs和Indels变异数目等，包括转换，颠换数目；

sample.stat.file部分内容如下图所示，number of SNPs为SNPs总数目，number of indels为发生indels的总数目，ts表示转换数目，tv表示颠换数目：

bcftools view –v snps –g hom sample.filtered.vcf.gz# 仅输出纯合SNPs，通过改变参数可以输出杂合SNPs，纯合/杂合Indels；

参数说明：

-h 仅输出头部注释文件；

-H 不输出头部注释文件；

-O --output-type <b|u|z|v> 指定输出文件类型， 'b'表示压缩后BCF文件; 'u'表示未压缩的BCF文件; 'z'表示压缩的VCF文件; 'v' 表示未压缩的VCF文件；

-T, --threads 设置线程数目

-i --Include <exp> 筛选出满足exp条件的变异位点；

-e --exclude<exp> 剔除满足exp条件的变异位点，与-i相反；

-v --types [snps|indels|mps|other]，仅输出指定类型的变异位点信息；

-V --exclude-type [snps|indels|mps|other]，仅输出不包含指定类型的变异位点信息；

-g --genotype [hom|het|miss]，仅输出指定基因型的位点。

四、参考

Samtools官方说明文档：http://www.htslib.org/doc/#manual-pages

Bcftools官方说明文档：http://samtools.github.io/bcftools/bcftools.html

参考文献：

Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup (2009)The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics, 25, 2078-9. [PMID: 19505943]
Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011 Nov 1;27(21):2987-93. Epub 2011 Sep 8. [PMID: 21903627]doi:10.1093/bioinformatics/btr509;
Danecek P., Schiffels S., Durbin R. Multiallelic calling model in bcftools (-m) [link]
Li H. Improving SNP discovery by base alignment quality. Bioinformatics. 2011 Apr 15;27(8):1157-8. doi: 10.1093/bioinformatics/btr076. Epub 2011 Feb 13. [PMID: 21320865]
Durbin R. Segregation based metric for variant call QC [link]
Li H, Mathematical Notes on SAMtools Algorithms [link]
The Genomes Project Consortium. An integrated map of genetic variation from 1,092 human genomes[J]. Nature, 2012, 491(7422):56-65. doi: 10.1038/nature11632
Wu Y, Jing R, Dong Y, et al. Functional annotation of sixty-five type-2 diabetes risk SNPs and its application in risk prediction[J]. Scientific Reports, 2017, 7.