bcftools使用总结

使用bcftools进行SNP calling

bcftools也能够进行SNP calling。在以前的版本中，一般都是和samtools的mpileup命令结合使用, 命令以下

samtools mpileup -uf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf

因为samtools和bcftools更新得都很快，只要有一个版本不对，采用上面的pipeline就会报错。为了减小版本不合适带来的问题，bcftools的开发团队将mpileup这个功能添加到了bcftools中。

在最新版的bcftools 中，只须要使用bcftools这一个工具就能够实现SNP calling， 用法以下

bcftools mpileup mpileup.1.bam --fasta-ref mpileup.ref.fa | bcftools call -mv -o raw.vcf

--fasta-ref参数指定参考序列的fasta文件，mpileup.bam是输入文件，一般都是GATK 标准预处理流程获得的bam文件。

须要注意的是mpileup命令虽然也会输出VCF格式的文件，可是并不直接进行snp calling。下面的命令能够生成VCF格式的文件

bcftools mpileup mpileup.1.bam --fasta-ref mpileup.ref.fa >mpileup.vcf

直接生成的VCF文件内容以下

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100
17 1 . A <*> 0 . DP=5; PL
17 2 . A <*> 0 . DP=5; PL
17 3 . G <*> 0 . DP=5; PL
17 4 . C <*> 0 . DP=5; PL
17 5 . T <*> 0 . DP=5; PL

里面的每一条记录并非一个SNP位点，而是染色体上每一个碱基的比对状况的汇总。这种信息官方称之为genotype likelihoods。

call命令才是真正的执行SNP calling的程序，基本用法以下

bcftools call mpileup.vcf -c  -v -o variants.vcf

在进行SNP calling 时，必须选择一种算法，有两种calling算法可供选择，分别对应-c和-m参数。-c参数对应consensus-caller算法， -m参数对应multiallelic-caller算法，后者更适合多种allel和罕见变异的calling。

-v参数也是经常使用参数，做用是只输出变异位点的信息，若是一个位点不是snp/indel, 不会输出。

注：新版本bcftools中 call命令可替代view命令

--------------------------------------------------------------------------------------------------------------

其余命令

1. annotate

annotate命令有两个用途，第一个是用于注释VCF文件，用法以下

bcftools annotate -a db.vcf -c ID,QUAL,+TAG view.vcf -o annotate.vcf

-a参数指定注释用的数据库文件，格式能够是VCF, BED, 或者是\t分隔的自定义文件。在\t分隔的自定义文件中，必须包含CHROM, POS字段；-c参数指定将数据库的哪些信息添加到输出文件中。

第二个用途是编辑VCF文件，好比去除其中的某些注释信息，或者去除某些样本，用法以下

bcftools annotate -x ID,INFO/DP,FORMAT/DP  view.vcf -o remove.id.vcf

-x参数表示去除VCF文件中的注释信息，能够是其中的某一列，好比ID, 也能够是某些字段，好比INFO/DP，多个字段的信息用逗号分隔；去除以后，这些信息所在的列并不会去除，而是用.填充。

2. concat

concat命令能够将多个VCF文件合并为一个VCF文件，要求输入的VCF文件必须是排序以后的，若是包含多个样本的信息，样本的顺序也必须一致。经典的应用场景包括合并不一样染色体上的VCF文件，合并SNP和INDEL 两种类型的VCF文件，用法以下

bcftools concat merge.2.a.vcf.gz merge.3.a.vcf.gz -o -o merge.vcf

还须要注意一点，输入的VCF文件必须是通过bgzip压缩的文件。

3. merge

merge命令也是用于合并VCF文件，主要用于将单个样本的VCF文件合并成一个多个样本的VCF文件。用法以下

bcftools merge merge.a.vcf.gz merge.b.vcf.gz -o merge.vcf

该命令要求输入文件必须是通过bgzip压缩的文件， 并且还须要有.tbi的索引。

4. isec

isec用于在多个VCF文件之间取交集，差集，并集等操做，经典的应用场景是对多种软件的SNP calling 结果进行venn 分析。用法以下

bcftools isec A.vcf.gz B.vcf.gz -p dir

默认参数就是取交集，更多高级用法请参考帮助文档。

5. stats

stats命令用于统计VCF文件的基本信息，好比突变位点的总数，不一样类型突变位点的个数等。用法以下

bcftools stats view.vcf >  view.stats

输出文件中记录了不少类型的统计数据，重点介绍如下几种

基本信息：

SN 0 number of samples: 3
SN 0 number of records: 15
SN 0 number of no-ALTs: 1
SN 0 number of SNPs: 11
SN 0 number of MNPs: 0
SN 0 number of indels: 3
SN 0 number of others: 0
SN 0 number of multiallelic sites: 1
SN 0 number of multiallelic SNP sites: 0

颠换和转换的比例：

# TSTV, transitions/transversions:
# TSTV  [2]id  [3]ts  [4]tv  [5]ts/tv  [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV  0  8  3  2.67  8  3  2.67

Indel 长度分布:

# IDD, InDel distribution:
# IDD [2]id [3]length (deletions negative) [4]count
IDD 0 -2 1
IDD 0 1 2
IDD 0 3 1

碱基替换类型：

# ST, Substitution types:
# ST [2]id [3]type [4]count
ST 0 A>C 0
ST 0 A>G 0
ST 0 A>T 0
ST 0 C>A 1
ST 0 C>G 0
ST 0 C>T 4
ST 0 G>A 1
ST 0 G>C 1
ST 0 G>T 1
ST 0 T>A 0
ST 0 T>C 3
ST 0 T>G 0

输出文件能够用于plot-vcfstats命令，进行可视化, 这个脚本位于bcftools安装目录的misc目录下。用法以下

plot-vcfstats view.stats -p output

-p参数指定输出结果的目录，这个脚本依赖latex 生成pdf 文件，因此系统上的latext 必定要安装好。

输出目录下文件不少，详细列表以下

├── counts_by_af.indels.dat
├── counts_by_af.snps.dat
├── depth.0.dat
├── depth.0.pdf
├── depth.0.png
├── indels.0.dat
├── indels.0.pdf
├── indels.0.png
├── plot.py
├── plot-vcfstats.log
├── substitutions.0.pdf
├── substitutions.0.png
├── summary.aux
├── summary.log
├── summary.pdf
├── summary.tex
├── tstv_by_af.0.dat
└── tstv_by_qual.0.dat

主要看summary.pdf文件就能够了，该文件包含了不少信息

1.不一样类型的突变位点汇总

2.插入缺失长度分布图

3.测序深度分布

4.碱基转换类型分布

6. index

index命令用于对VCF文件创建索引，要求输入的VCF文件必须是使用bgzip压缩以后的文件，支持.csi和.tbi两种索引，默认状况下创建的索引是.csi格式， 用法以下

bgzip view.vcf
bcftools index view.vcf.gz

运行成功后，会生成索引文件view.vcf.gz.csi。若是须要创建.tbi格式的索引，用法以下

bcftools index -t view.vcf.gz

tbi索引文件为view.vcf.gz.tbi。

7.  view

view命令能够用于VCF和BCF格式的转换，用法以下

bcftools view view.vcf.gz -O u -o view.bcf

-O参数指定输出文件的类型，b表明压缩后的BCF文件，u表明未经压缩的BCF文件，z表明压缩后的VCF文件，v表明未经压缩的VCF文件；-o参数指定输出文件的名字。

还能够根据样本筛选VCF文件，用法以下

bcftools view view.vcf.gz -s NA00001,NA00002  -o subset.vcf

-s参数指定想要保留的样本信息，多个样本用逗号分隔。

若是样本名称添加了^前缀，表明去除这些样本，好比-s ^NA00001,NA00002，这个写法表示从VCF文件中去除NA00001,NA00002这两个样本的信息。

还能够过滤突变位点，过滤的条件很是多，能够根据突变位点的类型，基因型类型等等条件进行过滤，详细的参数能够参考软件的帮助文档，这里只作一个基本示例

bcftools view view.vcf.gz -k -o known.vcf

-k参数表示筛选已知的突变位点，即ID那一列值不是.的突变位点。

8. query

query命令也是用于格式转换，和view命令不一样，query经过表达式来指定输出格式，可定制化程度更高。用法以下

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' view.vcf.gz

-f参数经过一个表达式来指定输出格式，其中变量的写法以下

1. %CHROM
   表明VCF文件中染色体那一列，其余的列，好比POS, ID, REF, ALT, QUAL, FILTER也是相似的写法

2. []
   对于FORMAT字段的信息，必需要中括号括起来

3. %SAMPLE
   表明样本名称

4. %GT
   表明FORMAT字段中genotype的信息

5. \t
   表明制表符分隔

6. \n
   表明新的一行

    

输出文件以下

11 2343543 A . NA00001=0/0 NA00002=0/0 NA00003=0/0
11 5464562 C T NA00001=./. NA00002=./. NA00003=./.
20 76962 T C NA00001=0/1 NA00002=1/1 NA00003=1/1

更多变量的写法请参考官方文档。

9. sort

sort 命令用于对VCF文件排序， 按照染色体位置进行排序，用法以下

bcftools sort view.vcf.gz -o sort.view.vcf

10. reheader

reheader命令有两个用途，第一用途用于编辑VCF文件的头部，第二个用途用于替换VCF文件中的样本名。

替换样本的用法以下

bcftools reheader -s sample.file view.vcf -o new.sample.vcf

-s参数指定须要替换的样本名，内容以下

NA00001 NA1
NA00002 NA2
NA00003 NA3

第一列表明VCF文件中原始的样本名称，第二列表明替换后的样本名称，两类之间用空格分隔，须要注意的是，样本名不容许有空格。

编辑VCF文件的头部的用法以下

bcftools reheader -h header.file  view.vcf -o new.header.vcf

-h参数指定新的header文件，内容以下

##fileformat=VCFv4.3
##reference=file:///seq/references/1000Genomes-NCBI37.fasta
##contig=<ID=11,length=135006516>
##contig=<ID=20,length=63025520>
....