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命令有兩個(gè)用途：
（１）注釋VCF文件，用法如下
$ bcftools annotate -a db.vcf -c ID,QUAL,+TAG view.vcf -o annotate.vcf
-a參數(shù)指定注釋用的數(shù)據(jù)庫(kù)文件，格式可以是VCF, BED, 或者是\t分隔的自定義文件。在\t分隔的自定義文件中，必須包含CHROM, POS字段；
-c參數(shù)指定將數(shù)據(jù)庫(kù)的哪些信息添加到輸出文件中。
（２）編輯VCF文件，比如去除其中的某些注釋信息，或者去除某些樣本，用法如下
$ bcftools annotate -x ID,INFO/DP,FORMAT/DP view.vcf -o removed.id.vcf
-x參數(shù)表示去除VCF文件中的注釋信息，可以是其中的某一列，比如ID, 也可以是某些字段，比如INFO/DP，多個(gè)字段的信息用逗號(hào)分隔；去除之后，這些信息所在的列并不會(huì)去除，而是用.填充。
2. concat
concat命令可以將多個(gè)VCF文件合并為一個(gè)VCF文件，要求輸入的VCF文件必須是排序之后的，如果包含多個(gè)樣本的信息，樣本的順序也必須一致。經(jīng)典的應(yīng)用場(chǎng)景包括合并不同染色體上的VCF文件，合并SNP和INDEL 兩種類型的VCF文件，用法如下
$ bcftools concat a.vcf.gz b.vcf.gz -o merge.vcf
注意：輸入的VCF文件必須是經(jīng)過(guò)bgzip壓縮的文件。
3. merge
merge命令也是用于合并VCF文件，主要用于將單個(gè)樣本的VCF文件合并成一個(gè)多個(gè)樣本的VCF文件。用法如下
$ bcftools merge a.vcf.gz b.vcf.gz -o merge.vcf
注意：輸入文件必須是經(jīng)過(guò)bgzip壓縮的文件， 而且還需要有.tbi的索引。
concat可以進(jìn)行vcf的“縱”向合并，比如不同染色體的vcf文件，或者snp和indel進(jìn)行的合并。但是樣品順序必須保持一致。
merge可以進(jìn)行vcf的“橫”向合并，比如單個(gè)樣本的vcf文件的合并。
concat和merge的共同點(diǎn)是輸入文件必須是bgzip壓縮，且有索引。
4. isec
isec用于在多個(gè)VCF文件之間取交集，差集，并集等操作，經(jīng)典的應(yīng)用場(chǎng)景是對(duì)多種軟件的SNP calling 結(jié)果進(jìn)行venn 分析。用法如下
$ bcftools isec a.vcf.gz b.vcf.gz -p dir
默認(rèn)參數(shù)就是取交集，更多高級(jí)用法請(qǐng)參考幫助文檔。
5. stats
stats命令用于統(tǒng)計(jì)VCF文件的基本信息。比如，突變個(gè)數(shù)、突變類型的個(gè)數(shù)、轉(zhuǎn)換顛換個(gè)數(shù)、測(cè)序深度、Indel長(zhǎng)度。還可以利用plot-vcfstats進(jìn)行可視化處理。用法如下
$ bcftools stats view.vcf > view.stats
輸出文件中記錄了很多類型的統(tǒng)計(jì)數(shù)據(jù)，重點(diǎn)介紹以下幾種
基本信息：
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
顛換和轉(zhuǎn)換的比例：
# 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 長(zhǎng)度分布:
# 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命令，進(jìn)行可視化, 這個(gè)腳本位于bcftools安裝目錄的misc目錄下。用法如下
$ plot-vcfstats view.stats -p output
-p參數(shù)指定輸出結(jié)果的目錄，這個(gè)腳本依賴latex 生成pdf 文件，所以系統(tǒng)上的latext 一定要安裝好。
輸出目錄下文件很多，詳細(xì)列表如下
├── 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.不同類型的突變位點(diǎn)匯總
2.插入缺失長(zhǎng)度分布圖
3.測(cè)序深度分布
4.堿基轉(zhuǎn)換類型分布
6. index
index命令用于對(duì)VCF文件建立索引，要求輸入的VCF文件必須是使用bgzip壓縮之后的文件，支持.csi和.tbi兩種索引，默認(rèn)情況下建立的索引是.csi格式， 用法如下
$ bgzip view.vcf
$ bcftools index view.vcf.gz
運(yùn)行成功后，會(huì)生成索引文件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)文件。前者為文本文件，后者為其二進(jìn)制文件。最主要的命令是view命令來(lái)進(jìn)行SNP和Indel calling。
$ bcftools view view.vcf.gz -O u -o view.bcf
-O參數(shù)指定輸出文件的類型；
-o參數(shù)指定輸出文件的名字；
u代表未經(jīng)壓縮的BCF文件；
b代表壓縮后的BCF文件；
v代表未經(jīng)壓縮的VCF文件；
z代表壓縮后的VCF文件；
還可以根據(jù)樣本篩選VCF文件，用法如下-s select
$ bcftools view view.vcf.gz -s NA00001,NA00002 -o subset.vcf
-s參數(shù)指定想要保留的樣本信息，多個(gè)樣本用逗號(hào)分隔。如果樣本名稱添加了^前綴，代表去除這些樣本，比如-s ^NA00001,NA00002，這個(gè)寫(xiě)法表示從VCF文件中去除NA00001,NA00002這兩個(gè)樣本的信息。
還可以過(guò)濾突變位點(diǎn)，過(guò)濾的條件非常多，可以根據(jù)突變位點(diǎn)的類型，基因型類型等等條件進(jìn)行過(guò)濾，詳細(xì)的參數(shù)可以參考軟件的幫助文檔，這里只做一個(gè)基本示例
$ bcftools view view.vcf.gz -k -o known.vcf
-k參數(shù)表示篩選已知的突變位點(diǎn)，即ID那一列值不是.的突變位點(diǎn)。
8. query
query命令也是用于格式轉(zhuǎn)換，和view命令不同，query通過(guò)表達(dá)式來(lái)指定輸出格式，可定制化程度更高。用法如下
$ bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' view.vcf.gz
-f參數(shù)通過(guò)一個(gè)表達(dá)式來(lái)指定輸出格式，其中變量的寫(xiě)法如下
- %CHROM 代表VCF文件中染色體那一列，其他的列，比如
POS,ID,REF,ALT,QUAL,FILTER也是類似的寫(xiě)法 - [] 對(duì)于FORMAT字段的信息，必須要中括號(hào)括起來(lái)
- %SAMPLE 代表樣本名稱
- %GT 代表FORMAT字段中g(shù)enotype的信息
- \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
更多變量的寫(xiě)法請(qǐng)參考官方文檔。
9. sort
sort 命令用于對(duì)VCF文件排序， 按照染色體位置進(jìn)行排序，用法如下
$ bcftools sort view.vcf.gz -o sort.view.vcf
10. reheader
reheader命令有兩個(gè)用途，第一用途用于編輯VCF文件的頭部，第二個(gè)用途用于替換VCF文件中的樣本名。
(1) 替換樣本的用法如下
$ bcftools reheader -s sample.file view.vcf -o new.sample.vcf
-s參數(shù)指定需要替換的樣本名，內(nèi)容如下
NA00001 NA1
NA00002 NA2
NA00003 NA3
第一列代表VCF文件中原始的樣本名稱，第二列代表替換后的樣本名稱，兩類之間用空格分隔，需要注意的是，樣本名不允許有空格。
(2) 編輯VCF文件的頭部的用法如下
$ bcftools reheader -h header.file view.vcf -o new.header.vcf
-h參數(shù)指定新的header文件，內(nèi)容如下
##fileformat=VCFv4.3
##reference=file:///seq/references/1000Genomes-NCBI37.fasta
##contig=<ID=11,length=135006516>
##contig=<ID=20,length=63025520>
....
11.call、cnv
變異檢測(cè)和cnv檢測(cè)
常用例子
(1) 左對(duì)齊標(biāo)準(zhǔn)化Indel，對(duì)于多等位基因位點(diǎn)進(jìn)行拆分（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"
(2) 根據(jù)條件進(jìn)行篩選，比如：篩選出FILTER列是PASS的，DP >20 , GQ >100, QUAL >100的位點(diǎn)：
$ bcftools filter -i ' FILTER=="PASS" && DP>20 && GQ>100 && QUAL>100' "$outDIR"_norm/"$out_basen
(3) 提取 突變位點(diǎn)的AD和DP，順便可以計(jì)算 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
# 計(jì)算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
# 過(guò)濾VAF(獲得VAF > 0.3的位點(diǎn))
$ 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) 合并多個(gè)樣本的VCF文件, 注意需要每個(gè)文件的樣本名唯一，如果不唯一使用--force-samples 將自動(dòng)重命名。
# 定義用于存儲(chǔ)待合并的vcf文件路徑：
combine_vars=""
for((i=0;i<${#filtered_vcf[*]};i++)){
#為每個(gè)文件建立索引，如果沒(méi)壓縮要先壓縮
tabix -p vcf "${filtered_vcf[$i]}"
echo "完成了第"$i"個(gè)"
combine_vars="${combine_vars}"" ""${filtered_vcf[$i]}"
}
# --missing-to-ref表示缺失的GT表示為0/0，-merge操作多等位基因位點(diǎn)，這里表示不產(chǎn)生多等位基因位點(diǎn)。
# 當(dāng)然對(duì)于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) 上面是合并多個(gè)樣本，如果是相同樣本的位點(diǎn)合并呢，比如一個(gè)樣本的SNP和INDEL進(jìn)行合并，首先必須對(duì)待合并文件進(jìn)行排序：
# 排序位點(diǎn)：
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) 提取指定染色體上的位點(diǎn)：
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來(lái)注釋 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分割的文件進(jìn)行注釋VCF(1-bae)：
# 需要 1-base的坐標(biāo)系并且建立索引：
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文件進(jìn)行注釋(0-base)：
bcftools annotate -a annots.bed.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf
(12) 提取指定樣本的vcf文件
準(zhǔn)備樣本ID文件，這里命名為samplelistname.txt，一個(gè)樣本一行，如下所示：
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
# 還是輸出上面的結(jié)果，但用\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
# 輸出樣本的突變位點(diǎn)信息和GT：
bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]' -i'GT="alt"
（14） 使用bcftools進(jìn)行SNP calling
#一開(kāi)始寫(xiě)好引用，方便以后
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比對(duì)，samtools排序并構(gòu)建索引，為了之后更快調(diào)用比對(duì)文件
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之前，需要先建立兩個(gè)參考基因組的索引文件.dict和.fai【具體參見(jiàn)https://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference】
#.dict中包含了基因組中contigs的名字，也就是一個(gè)字典；
#構(gòu)建.dict文件（原來(lái)要使用picard的CreateSequenceDictionary模塊，但是現(xiàn)在gatk整合了此模塊，可以直接使用）
gatk CreateSequenceDictionary -R $REF -O $ebola/ref/$ACC.dict
#.fai也就是fasta index file，索引文件，可以快速找出參考基因組的堿基
#構(gòu)建
samtools faidx $REF
#gatk開(kāi)始：
#必選 -I -O -R，代表輸入、輸出、參考
#接下來(lái)可以按照字母順序依次寫(xiě)出來(lái)，這樣比較清晰
#-bamout：將一整套經(jīng)過(guò)gatk程序重新組裝的單倍體基因型（haplotypes）輸出到文件
#-stand-call-conf :低于這個(gè)數(shù)字的變異位點(diǎn)被忽略，可以設(shè)成標(biāo)準(zhǔn)30（默認(rèn)是10）
gatk HaplotypeCaller -R $REF -I $BAM -O $ebola/align/HaplotypeCaller.vcf \
-bamout $ebola/align/$SRR.gatk.bam \
-stand-call-conf 30
# gatk用時(shí)3.95 minutes.
#<gatk補(bǔ)充>GATK進(jìn)行變異檢測(cè)的時(shí)候，是按照染色體排序順序進(jìn)行的，并非多條染色體并行檢測(cè)的。因此，如果樣本數(shù)據(jù)量比較大的話，一般多個(gè)染色體并行。
bcftools也可以進(jìn)行SNP calling。在之前的版本中，通常都是和samtools的mpileup命令結(jié)合使用。首先，對(duì)排序好的bam數(shù)據(jù)用samtools生成bcf文件。然后，由于生成的是二進(jìn)制格式的數(shù)據(jù)，需要進(jìn)行解析或者轉(zhuǎn)換成vcf：
samtools mpileup -uf ref.fa aln.bam | bcftools view var.raw.vcf
由于samtools和bcftools更新得都很快，只要有一個(gè)版本不對(duì)，采用上面的pipeline就會(huì)報(bào)錯(cuò)。為了減少版本不合適帶來(lái)的問(wèn)題，bcftools的開(kāi)發(fā)團(tuán)隊(duì)將mpileup這個(gè)功能添加到了bcftools中。
在最新版的bcftools 中，只需要使用bcftools這一個(gè)工具就可以實(shí)現(xiàn)SNP calling， 用法如下
bcftools mpileup mpileup.1.bam --fasta-ref mpileup.ref.fa | bcftools call -mv -o raw.vcf
--fasta-ref參數(shù)指定參考序列的fasta文件，mpileup.bam是輸入文件，通常都是GATK 標(biāo)準(zhǔn)預(yù)處理流程得到的bam文件。
需要注意的是mpileup命令雖然也會(huì)輸出VCF格式的文件，但是并不直接進(jìn)行snp calling。下面的命令可以生成VCF格式的文件
bcftools mpileup mpileup.1.bam --fasta-ref mpileup.ref.fa >mpileup.vcf
直接生成的VCF文件內(nèi)容如下
#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
里面的每一條記錄并不是一個(gè)SNP位點(diǎn)，而是染色體上每個(gè)堿基的比對(duì)情況的匯總。這種信息官方稱之為genotype likelihoods。
call命令才是真正的執(zhí)行SNP calling的程序，基本用法如下
bcftools call mpileup.vcf -c -v -o variants.vcf
在進(jìn)行SNP calling 時(shí)，必須選擇一種算法，有兩種calling算法可供選擇，分別對(duì)應(yīng)-c和-m參數(shù)。-c參數(shù)對(duì)應(yīng)consensus-caller算法， -m參數(shù)對(duì)應(yīng)multiallelic-caller算法，后者更適合多種allel和罕見(jiàn)變異的calling。
-v參數(shù)也是常用參數(shù)，作用是只輸出變異位點(diǎn)的信息，如果一個(gè)位點(diǎn)不是snp/indel, 不會(huì)輸出。
注：新版本bcftools中 call命令可替代view命令
REF:
https://www.cnblogs.com/emanlee/p/4316581.html
https://msd.misuland.com/pd/3255818135034402688
http://www.itdecent.cn/p/b3a0d1448a36
