bcftools常用命令详解

http://samtools.github.io/bcftools/bcftools.html

https://www.biostars.org/p/95013/

下载

wget https://github.com/samtools/bcftools/releases/download/1.10/bcftools-1.10.tar.bz2
tar xjvf bcftools-1.10.tar.bz2
cd bcftools-1.10/
./configure --prefix=安装地址
make
make install

1. annotate

annotate命令有两个用途:

(1)注释VCF文件,用法如下

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

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

(2)编辑VCF文件,比如去除其中的某些注释信息,或者去除某些样本,用法如下

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

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

2. concat

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

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

注意:输入的VCF文件必须是经过bgzip压缩的文件。

3. merge

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

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

注意:输入文件必须是经过bgzip压缩的文件, 而且还需要有.tbi的索引。

concat可以进行vcf的“纵”向合并,比如不同染色体的vcf文件,或者snp和indel进行的合并。但是样品顺序必须保持一致。
merge可以进行vcf的“横”向合并,比如单个样本的vcf文件的合并。
concat和merge的共同点是输入文件必须是bgzip压缩,且有索引。

4. isec

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

$ bcftools isec a.vcf.gz b.vcf.gz -p dir

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

5. stats

stats命令用于统计VCF文件的基本信息。比如,突变个数、突变类型的个数、转换颠换个数、测序深度、Indel长度。还可以利用plot-vcfstats进行可视化处理。用法如下

$ 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(variant call format)文件和bcf(binary call format)文件。前者为文本文件,后者为其二进制文件。最主要的命令是view命令来进行SNP和Indel calling。

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

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

还可以根据样本筛选VCF文件,用法如下-s select

$ 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 制表符分隔,\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文件中的样本名。

(1) 替换样本的用法如下

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

-s参数指定需要替换的样本名,内容如下

NA00001 NA1
NA00002 NA2
NA00003 NA3

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

(2) 编辑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>
....

11.call、cnv

变异检测和cnv检测

常用例子

(1) 左对齐标准化Indel,对于多等位基因位点进行拆分(annovar注释必须的)。

#首先需要压缩VCF并建立索引。
$ bgzip -f "$outDIR"_tmp/"$out_basename".vcf -@ 10
$ tabix -p vcf "$outDIR"_tmp/"$out_basename".vcf.gz
$ bcftools norm --fasta-ref "$fasta_file" --multiallelics -both --output "$outDIR"_norm/"$out_basename" --output-type z "$inputVCF"
image.png

(2) 根据条件进行筛选,比如:筛选出FILTER列是PASS的,DP >20 , GQ >100, QUAL >100的位点:

$ bcftools filter  -i ' FILTER=="PASS" &&  DP>20 &&  GQ>100 && QUAL>100' "$outDIR"_norm/"$out_basen

(3) 提取 突变位点的AD和DP,顺便可以计算 VAF:

$ bcftools query -f '[%AD]\n' "$outDIR"/"$out_basename" |awk 'BEGIN{FS=","}{ print $2 }'  >  > "$outDIR"_tmp/"$out_basename"_AD_tmp.txt
$ bcftools query -f '[%DP]\n' "$outDIR"/"$out_basename"  > "$outDIR"_tmp/"$out_basename"_DP_tmp.txt

# 计算VAF:
$ awk 'BEGIN{OFS="\t"}{ if(NR==FNR)AD[NR]=$1; if(NR>FNR)DP[FNR]=$1; }END{ for(i=1;i<=length(AD);i++){ print AD[i],DP[i],AD[i]/DP[i] } }' $"$outDIR"_tmp/"$out_basename"_AD_tmp.txt "$outDIR"_tmp/"$out_basename"_DP_tmp.txt > "$outDIR"_tmp/"$out_basename"_AD_DP_VAF_tmp.txt

# 过滤VAF(获得VAF > 0.3的位点)
$ num_anno=`zcat "$outDIR"/"$out_basename" |awk '{if($0~"^#")print 1 }' |wc -l`
$ zcat "$outDIR"/"$out_basename" |awk -v num_anno="$num_anno" 'BEGIN{FS="\t"}{ if(NR==FNR && $0~"^#"){print $0}; if(NR==FNR && $0!~"^#")vcf[NR-num_anno]=$0; if(NR>FNR && $3>0.3)print vcf[FNR] }'  - "$outDIR"_tmp/"$out_basename"_AD_DP_VAF_tmp.txt > "$outDIR"_tmp/"$out_basename"_VAF_filtered.vcf

(4) 合并多个样本的VCF文件, 注意需要每个文件的样本名唯一,如果不唯一使用--force-samples 将自动重命名。

# 定义用于存储待合并的vcf文件路径:
combine_vars=""
for((i=0;i<${#filtered_vcf[*]};i++)){
  #为每个文件建立索引,如果没压缩要先压缩
  tabix -p vcf "${filtered_vcf[$i]}"
  echo "完成了第"$i"个"
  combine_vars="${combine_vars}"" ""${filtered_vcf[$i]}"
}
#  --missing-to-ref表示缺失的GT表示为0/0,-merge操作多等位基因位点,这里表示不产生多等位基因位点。
# 当然对于header使用--use-header指定,info列的合并使用--info-rules指定规则。
bcftools merge --missing-to-ref --merge none --output "$workdir"/test/combined.vcf.gz --output-type z --threads 30  

(5) 上面是合并多个样本,如果是相同样本的位点合并呢,比如一个样本的SNP和INDEL进行合并,首先必须对待合并文件进行排序:

# 排序位点:
bcftools sort SNP_filtered.vcf -O z -o SNP_filtered_sorted.vcf.gz
bcftools sort INDEL_filtered.vcf -O z -o INDEL_filtered_sorted.vcf.gz
# 合并:
bcftools concat SNP_filtered_sorted.vcf.gzINDEL_filtered_sorted.vcf.gz  -a -O z -o ALL_filtered_sorted.vcf.gz

(6) 提取指定染色体上的位点:

bcftools filter -t chr1,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr2,chr20,chr21,chr22,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chrM,chrX,chrY  "$workdir"/test/combined_split.vcf.gz  --output "$workdir"/test/combined_split_chr.vcf.gz --output-type z 

(7) 移除INFO和FORMAT中除了GT和PL的列:

bcftools annotate -x INFO,^FORMAT/GT,FORMAT/PL file.vcf

(8) 使用 src.bcf来注释 dst.bcf,只注释ID,QUAL和TAG,如果TAG存在则不覆盖。

bcftools annotate -a src.bcf -c ID,QUAL,+TAG dst.bcf

(9) 除了FORMAT的GT列外,注释所有的INFO和FORMAT.

bcftools annotate -a src.bcf -c INFO,^FORMAT/GT dst.bcf

(10) 使用TAB分割的文件进行注释VCF(1-bae):

# 需要 1-base的坐标系并且建立索引:
tabix -s1 -b2 -e2 annots.tab.gz
bcftools annotate -a annots.tab.gz -h annots.hdr -c CHROM,POS,REF,ALT,-,TAG file.vcf
bcftools annotate -a annots.tab.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf

(11) 使用bed文件进行注释(0-base):

bcftools annotate -a annots.bed.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf

(12) 提取指定样本的vcf文件
准备样本ID文件,这里命名为samplelistname.txt,一个样本一行,如下所示:

sample1
sample2
sample3
bcftools view -S samplelistname.txt  /genomes/ALL..genotypes.vcf.gz -Ov > samplelist_1000Genomes.vcf`

(13) 常用查询命令:

# 输出染色体、位置、REF、ALT:
bcftools query -f '%CHROM  %POS  %REF  %ALT{0}\n' file.vcf.gz

# 还是输出上面的结果,但用\t代替空格并输出样本名和基因型
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' file.vcf.gz

# 输出GQ和GT:
bcftools query -f 'GQ:[ %GQ] \t GT:[ %GT]\n' file.vcf

# 创建bed文件: chr, pos (0-based), end pos (1-based), id
bcftools query -f'%CHROM\t%POS0\t%END\t%ID\n' file.bcf

# 输出样本的突变位点信息和GT:
bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]' -i'GT="alt"

(14) 使用bcftools进行SNP calling

#一开始写好引用,方便以后
ACC=AF086833
ebola=/vol2/agis/xiaoyutao_group/liuyunze/project/ebola
REF=$ebola/ref/$ACC.fa
SRR=SRR1553500
BAM=$ebola/align/$SRR.bam
R1=$ebola/raw/${SRR}_1.fastq
R2=$ebola/raw/${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
VARI=$ebola/variant

##bwa比对,samtools排序并构建索引,为了之后更快调用比对文件
mkdir -p $ebola/align && cd $ebola/align
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM

mkdir -p $VARI
samtools faidx $REF

##第一种方法:bcftools召唤变异
samtools mpileup -uvf $REF $BAM | bcftools call -vm -Oz > bcftools.vcf.gz

##第二种方法:freebayes
freebayes -f $REF $BAM > $ebola/align/freebayes.vcf

##第三种方法:GATK(版本:4.0.7.0)
#注意:在使用GATK之前,需要先建立两个参考基因组的索引文件.dict和.fai【具体参见https://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference】
#.dict中包含了基因组中contigs的名字,也就是一个字典;
#构建.dict文件(原来要使用picard的CreateSequenceDictionary模块,但是现在gatk整合了此模块,可以直接使用)
gatk CreateSequenceDictionary -R $REF -O $ebola/ref/$ACC.dict
#.fai也就是fasta index file,索引文件,可以快速找出参考基因组的碱基
#构建
samtools faidx $REF
#gatk开始:
#必选 -I -O -R,代表输入、输出、参考
#接下来可以按照字母顺序依次写出来,这样比较清晰
#-bamout:将一整套经过gatk程序重新组装的单倍体基因型(haplotypes)输出到文件
#-stand-call-conf :低于这个数字的变异位点被忽略,可以设成标准30(默认是10)
gatk HaplotypeCaller -R $REF -I $BAM -O $ebola/align/HaplotypeCaller.vcf \
-bamout $ebola/align/$SRR.gatk.bam \
-stand-call-conf 30 
# gatk用时3.95 minutes.
#<gatk补充>GATK进行变异检测的时候,是按照染色体排序顺序进行的,并非多条染色体并行检测的。因此,如果样本数据量比较大的话,一般多个染色体并行。

bcftools也可以进行SNP calling。在之前的版本中,通常都是和samtools的mpileup命令结合使用。首先,对排序好的bam数据用samtools生成bcf文件。然后,由于生成的是二进制格式的数据,需要进行解析或者转换成vcf:

samtools mpileup -uf ref.fa aln.bam | bcftools view var.raw.vcf

由于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命令

REF:
https://www.cnblogs.com/emanlee/p/4316581.html
https://msd.misuland.com/pd/3255818135034402688
https://www.jianshu.com/p/b3a0d1448a36

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 143,540评论 1 302
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 61,539评论 1 258
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 95,005评论 0 213
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 41,159评论 0 180
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 48,946评论 1 258
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 38,864评论 1 178
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 30,468评论 2 273
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,218评论 0 167
  • 想象着我的养父在大火中拼命挣扎,窒息,最后皮肤化为焦炭。我心中就已经是抑制不住地欢快,这就叫做以其人之道,还治其人...
    爱写小说的胖达阅读 29,070评论 6 234
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 32,574评论 0 213
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,353评论 2 215
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 30,683评论 1 232
  • 白月光回国,霸总把我这个替身辞退。还一脸阴沉的警告我。[不要出现在思思面前, 不然我有一百种方法让你生不如死。]我...
    爱写小说的胖达阅读 24,270评论 0 32
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,168评论 2 214
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 31,591评论 3 210
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,638评论 0 9
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,035评论 0 166
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 33,590评论 2 232
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 33,685评论 2 233