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命令有兩個用途：

（１）注釋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 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參數通過一個表達式來指定輸出格式敦锌，其中變量的寫法如下

%CHROM 代表VCF文件中染色體那一列，其他的列佳簸，比如POS, ID, REF, ALT, QUAL, FILTER也是類似的寫法
[] 對于FORMAT字段的信息乙墙，必須要中括號括起來
%SAMPLE 代表樣本名稱
%GT 代表FORMAT字段中genotype的信息
\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指定規(guī)則铛纬。
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

# 創(chuàng)建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的開發(fā)團隊將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命令才是真正的執(zhí)行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
http://www.reibang.com/p/b3a0d1448a36