一芝此、 基本介紹

VCF格式（Variant Call Format）是存儲(chǔ)變異位點(diǎn)的標(biāo)準(zhǔn)格式，用于記錄variants（SNP / InDel）因痛。BCF是VCF的二進(jìn)制文件婚苹。bcftools可以實(shí)現(xiàn)SNP calling。

二鸵膏、 背景知識(shí)

(1) 變異和突變

變異膊升，指的是實(shí)際測(cè)序數(shù)據(jù)與國(guó)際規(guī)定的參考基因組之間的區(qū)別。很多變異其實(shí)只是造成人類多樣性的原因谭企。突變廓译，指的是那些與疾病相關(guān)的變異。舉個(gè)例子：ENSEMBL等規(guī)定的人類參考基因組文件某位置是AAAAA债查，然后一個(gè)人實(shí)際測(cè)序得到的序列為AGCAA非区，那么相比于參考基因組，這個(gè)人就有2個(gè)變異位點(diǎn)攀操。對(duì)于第2個(gè)位置院仿，如果查看所有已知的測(cè)序，絕大部分人都是G，說(shuō)明是參考基因組出現(xiàn)了問(wèn)題歹垫，這個(gè)變異就不能稱作突變剥汤。對(duì)于第3個(gè)位置，如果查看所有已知的測(cè)序排惨，絕大部分人都是A吭敢，而恰好有一個(gè)人不是A，但他是個(gè)患者暮芭，那么這個(gè)變異就是突變了鹿驼。

(2) 變異類型

• SNP（single nucleotide polymorphism）：?jiǎn)魏塑账岫鄳B(tài)性。個(gè)體間基因組DNA序列同一位置單個(gè)核苷酸變異（替換辕宏、插入或缺失）所引起的多態(tài)性畜晰。在人類基因組中SNP分布普遍并且密度較大，總數(shù)超過(guò)10^7瑞筐， 平均每300bp（也有說(shuō)1kbp）就有一個(gè)SNP凄鼻。或稱單核苷酸位點(diǎn)變異SNV聚假。
• INDEL（insertion-deletion）：插入和缺失块蚌。基因組上小片段（>50bp）的插入或缺失膘格。
• CNV（copy number variation）：基因組拷貝數(shù)變異峭范。基因組中大片段的DNA形成非正常的拷貝數(shù)量瘪贱。比如一個(gè)基因在染色體的一條染色單體上的數(shù)目為1纱控，但是在染色體復(fù)制過(guò)程中，復(fù)制結(jié)束后該基因在染色單體數(shù)目由1變成了2或者n菜秦。它發(fā)生的頻率遠(yuǎn)遠(yuǎn)高于染色體結(jié)構(gòu)變異其徙，并且整個(gè)基因組中覆蓋的核苷酸總數(shù)大大超過(guò)SNP的總數(shù)。
• SV（structure variation）：結(jié)構(gòu)變異喷户。染色體大片段的插入與缺失唾那，染色體內(nèi)部的某區(qū)域發(fā)生翻轉(zhuǎn)顛換，兩條染色體之間發(fā)生重組褪尝。

(3) SNP

一般情況下只分析SNP闹获，其它類型的變異分析有難度或不準(zhǔn)確。來(lái)自兩個(gè)不同個(gè)體的DNA片段AAGCCTA和AAGCTTA為等位基因河哑。幾乎所有常見(jiàn)的SNP位點(diǎn)只有兩個(gè)等位基因避诽。在人體中，SNP的發(fā)生機(jī)率大約是0.1%璃谨，也就是每1000個(gè)堿基對(duì)就可能有一個(gè)SNP（密度高）沙庐。對(duì)疾病發(fā)生和藥物治療有重大影響的SNP鲤妥，估計(jì)只占數(shù)以百萬(wàn)計(jì)SNP的很小一部分。SNP位點(diǎn)的分布是不均勻的拱雏，在非轉(zhuǎn)錄序列比在轉(zhuǎn)錄序列更常見(jiàn)棉安。編碼區(qū)的單核苷酸多態(tài)性——編碼 SNP（coding SNP，cSNP）也有同義和非同義兩種類型铸抑，非同義SNP會(huì)改變蛋白質(zhì)的氨基酸序列贡耽。基因非編碼區(qū)鹊汛、基因間隔區(qū)的SNP仍然可能影響轉(zhuǎn)錄因子結(jié)合蒲赂、剪接等過(guò)程。從演化的觀點(diǎn)來(lái)看刁憋，SNP具有相當(dāng)程度的穩(wěn)定性滥嘴，即使經(jīng)過(guò)代代相傳，SNP所引起的改變卻不大至耻，因此可用以研究族群演化氏涩。

(4) vcf格式

vcf格式（Variant Call Format）是存儲(chǔ)變異位點(diǎn)的標(biāo)準(zhǔn)格式，用于記錄variants（SNP / InDel）有梆。BCF是VCF的二進(jìn)制文件。

• 以#開頭的注釋部分：

##fileformat：VCF格式版本號(hào)意系。
##reference & contig：使用的參考基因組信息及參考基因組contig信息泥耀。
##INFO行：是堿基位點(diǎn)的注釋。每一行必須的四個(gè)標(biāo)簽是：ID蛔添、Number痰催、Type、Description迎瞧。

• 沒(méi)有#開頭的主體部分：
  包含10列數(shù)據(jù)夸溶，每一行代表一個(gè)variant的信息。
  主體部分10列的范例：CHROM凶硅、POS缝裁、ID、REF足绅、ALT捷绑、QUAL、FILTER氢妈、INFO粹污、FORAMT、SAMPLE（前8列必須要有）首量。

例如：

chrM（CHROM染色體）
150（POS變異的第一個(gè)位置壮吩，1-based coordinate system）
.（ID變異位點(diǎn)名稱进苍，在dbSNP數(shù)據(jù)庫(kù)中的ID以rs開頭 ，一般只有人類基因組才有dbSNP編號(hào)鸭叙，如果沒(méi)有則為點(diǎn)）
T（REF參考序列該位置堿基類型和數(shù)目）
C（ALT該位置變異的堿基類型和數(shù)目觉啊，點(diǎn)代表缺失，多個(gè)用逗號(hào)分隔）
7766.77（QUAL變異的質(zhì)量值递雀，Q=-10lgP柄延，值越大是變異的可能性越大）
PASS（FILTER是否要被過(guò)濾掉，為PASS表示可能是變異缀程，點(diǎn)代表沒(méi)有進(jìn)行任何過(guò)濾）
AC=2; AF=1.00; AN=2; DP=199; ExcessHet=3.0103; FS=0.000; MLEAC=2; MLEAF=1.00; MQ=49.78; QD=32.91; SOR=0.904（INFO：variant的相關(guān)信息搜吧。）
GT:AD:DP:GQ:PL（FORAMT：variants的格式）
1/1:0,175:175:99:7795,531,0（一個(gè)SAMPLE為1列，總列數(shù)可以多于10杨凑，每列分別對(duì)應(yīng)第9列的各個(gè)格式滤奈，由bam文件中@RG的SM標(biāo)簽決定）

第8列INFO：
variant的相關(guān)信息。

AC（Allele Count）變異的等位基因數(shù)目
AF（Allel Frequency）等位基因頻率
AN（Allel Number）等位基因總數(shù)目
DP撩满，reads覆蓋度
FS蜒程，F(xiàn)ishers精確檢驗(yàn)的p值

AC（Allele Count）變異的等位基因數(shù)目。AF（Allele Frequency）等位基因頻率伺帘，AN（Allele Number）等位基因總數(shù)目昭躺。（0/1：AC=1，AF=0.5伪嫁，AN=2）
DP领炫，是一些reads被過(guò)濾掉后的覆蓋度。
DP4张咳，高質(zhì)量測(cè)序堿基帝洪，在ref或alt前后。
Dels脚猾，有這個(gè)tag且為0時(shí)表示該位點(diǎn)是SNV葱峡，沒(méi)有就是InDel。
FS龙助，使用Fisher精確檢驗(yàn)來(lái)檢測(cè)strand bias而得到的Fhred格式的p值砰奕。該值越小越好。一般進(jìn)行filter的時(shí)候提鸟，可以設(shè)置 FS < 10～20脆淹。

第9列FORMAT：
variants的格式。

GT（genotype）沽一，0代表樣本中ref的allel盖溺，1代表樣本variant的allel，2表示有第二個(gè)variant的allel铣缠；0/0表示樣品中該位點(diǎn)為純合位點(diǎn)烘嘱，和REF的堿基類型一致昆禽；0/1表示sample中該位點(diǎn)為雜合突變，AC=1蝇庭，AF=0.5哮内，AN=2北发；1/1表示為變異純合子，AC=2瞭恰，AF=1惊畏，AN=2颜启。
AD（Allele Depth）為sample中每一種allele（等位堿基）的reads覆蓋度。
DP（Depth）為sample中該位點(diǎn)的覆蓋度。
GQ（Genotype Quality）基因型的質(zhì)量值合呐，基因型存在的概率淌实。
PL（likelihood genotypes）指定的三種基因型的質(zhì)量值（0/0，0/1放坏，1/1）淤年，對(duì)應(yīng)的值越大溉苛，表示這種基因型的可能性越小愚战。

GT（genotype），0代表樣本中ref的allele敢茁，1代表樣本variant的allele彰檬，2表示有第二個(gè)variant的allele（幾乎所有常見(jiàn)的SNP位點(diǎn)只有兩個(gè)等位基因）。0/0表示樣品中該位點(diǎn)為純合位點(diǎn)较雕，和REF的堿基類型一致；0/1表示sample中該位點(diǎn)為雜合突變慎玖，AC=1，AF=0.5润努，AN=2；1/1表示為變異純合子随抠，AC=2二驰，AF=1，AN=2矗积。
AD（Allele Depth）對(duì)應(yīng)兩個(gè)以逗號(hào)隔開的值，這兩個(gè)值分別表示覆蓋到REF和ALT堿基的reads數(shù)，相當(dāng)于支持REF和支持ALT的測(cè)序深度茵烈。
DP（Depth）覆蓋到這個(gè)位點(diǎn)的總的reads數(shù)量，相當(dāng)于這個(gè)位點(diǎn)的深度存璃。
PL（likelihood genotypes）對(duì)應(yīng)3個(gè)以逗號(hào)隔開的值仑荐，指定的三種基因型（0/0，0/1，1/1）沒(méi)經(jīng)過(guò)先驗(yàn)的標(biāo)準(zhǔn)化Phred-scaled似然值篮迎。對(duì)應(yīng)的值越大示姿，表示這種基因型的可能性越小。
GQ（Genotype Quality）最可能的基因型的質(zhì)量值乃戈。

例如：
chr1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26

GT=0/1谍憔，也就是說(shuō)這個(gè)位點(diǎn)的基因型是C/T习贫。AD=1,3，也就是說(shuō)支持REF的read有一條千元，支持ALT的有3條祟身。DP=4，也就是說(shuō)只有4條reads支持這個(gè)地方的變異涕烧，cover到這個(gè)位點(diǎn)的reads數(shù)太少月而。GQ=25.92，質(zhì)量值并不算太高议纯。在PL里父款，這個(gè)位點(diǎn)基因型的不確定性就表現(xiàn)的更突出了，0/1的PL值為0瞻凤，雖然支持0/1的概率很高憨攒；但是1/1的PL值只有26，也就是說(shuō)還有10^(-2.6)=0.25%的可能性是1/1阀参；但幾乎不可能是0/0肝集，因?yàn)橹С?/0的概率只有10^(-10.3)=5*10^-11。

• 1-based coordinate system：序列的第一個(gè)堿基設(shè)為數(shù)字1蛛壳，如SAM, VCF, GFF, wiggle格式
• 0-based coordinate system：序列的第一個(gè)堿基設(shè)為數(shù)字0杏瞻，如BAM, BCFv2, BED, PSL格式

三、 用法和參數(shù)

(1) SNP calling

mpileup命令：得到染色體上每個(gè)堿基的比對(duì)情況的匯總（genotype likelihoods）

bcftools mpileup  -Ob  -o  sample.bcf  -f  dmel.genome.fa  sample.sorted.bam

輸入BAM文件sorted.bam
-f / --fasta-ref：指定參考序列的fasta文件
-O：指定輸出文件的類型衙荐，壓縮的BCF(b)捞挥，未壓縮的BCF(u)，壓縮的VCF(z)忧吟，未壓縮的VCF(v)
-o：指定輸出文件的名字sample.bcf

call命令：執(zhí)行SNP calling

bcftools call  -vmO  z  -o  sample.vcf.gz  sample.bcf

-v：只輸出變異位點(diǎn)的信息砌函，如果一個(gè)位點(diǎn)不是snp/indel則不會(huì)輸出
兩種calling算法：-c參數(shù)對(duì)應(yīng)consensus-caller算法，-m參數(shù)對(duì)應(yīng)multiallelic-caller算法，后者更適合多種allel和罕見(jiàn)變異的calling

(2) tabix

tabix  -p  vcf  sample.vcf.gz

輸入為壓縮文件vcf.gz讹俊，生成的索引文件為sample.vcf.gz.tbi

(3) index對(duì)vcf文件建立索引

bgzip  view.vcf  #輸入的VCF文件必須是bgzip壓縮后的文件
gunzip  view.vcf.gz  #解壓縮
bcftools index  view.vcf.gz  #生成索引文件view.vcf.gz.csi
bcftools index  -t  view.vcf.gz  #生成索引文件view.vcf.gz.tbi

(4) query通過(guò)表達(dá)式來(lái)指定輸出格式

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

-f：通過(guò)一個(gè)表達(dá)式來(lái)指定輸出格式
%CHROM：代表VCF文件中染色體那一列垦沉，其他的列POS, ID, REF, ALT, QUAL, FILTER也是類似的寫法
[]：對(duì)于FORMAT字段的信息用中括號(hào)括起來(lái)
%SAMPLE：代表樣本名稱
%GT：代表FORMAT字段中g(shù)enotype的信息
\t：代表制表符分隔
\n：代表新的一行

bcftools query  -r '19:400300-400800'  -f '%CHROM\t%POS\t%REF\t%ALT\n'  hg19.vcf.gz | head -3

-r：從特定區(qū)域提取varients信息，格式 chr|chr:pos|chr:from-to|chr:from-[,…]

bcftools query  -t ^'19:400300-400800'  -f '%CHROM\t%POS\t%REF\t%ALT\n'

排除特定區(qū)域

(5) sort按照染色體位置進(jìn)行排序

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

(6) filter過(guò)濾不可靠位點(diǎn)

bcftools filter  -O  z  -o  sample_filtered.vcf.gz  -s  LOWQUAL –i '%QUAL>10'  sample.vcf.gz

-O / --output-type：輸出的格式仍劈，z和v都行厕倍，壓縮的VCF(z)，未壓縮的VCF(v)
-o / --output：輸出文件的名稱
-s / --soft-filter：將過(guò)濾掉的位點(diǎn)用字符串注釋

bcftools filter  --no-version  -s FLTER  -i '(%QUAL>20 && format/DP>4 && MQ>30)||(GT="0/0")'  -Ov -o BL48384.Raw.flt.vcf  BL48384.Raw.vcf.gz
bcftools filter  -i 'FILTER=="PASS"'  --no-version  -Ov -o BL48384.Raw.flter.vcf  BL48384.Raw.flt.vcf

--no-version：不添加bcftools版本和命令到vcf頭文件
-s：注釋FILTER這列信息贩疙，過(guò)濾掉的信息為FLTER绑青，保留的為PASS
-i：篩選條件，篩選出QUAL大于20屋群、DP大于4闸婴、MQ大于30或者GT等于0/0的位點(diǎn)
-Ov：輸出文件為未壓縮vcf格式

(7) view命令用于VCF和BCF格式的轉(zhuǎn)換

bcftools view  view.vcf.gz  -O  u  -o  view.bcf
bcftools view  view.vcf.gz  -s  NA00001,NA00002  -o  subset.vcf
bcftools view  view.vcf.gz  -k  -o  known.vcf

-O：指定輸出文件的類型，壓縮的BCF(b)芍躏，未壓縮的BCF(u)邪乍，壓縮的VCF(z)，未壓縮的VCF(v)
-o：指定輸出文件的名字
-s：想要保留的樣本信息对竣，多個(gè)樣本用逗號(hào)分隔庇楞；如果樣本名稱添加了^前綴，代表去除這些樣本否纬，比如-s ^NA00001,NA00002
-k：表示篩選已知的突變位點(diǎn)吕晌，即ID那一列值不是.的突變位點(diǎn)

bcftools view  -i 'SAO=1'  b37.vcf  >  b37.germline.vcf  #選出INFO中SAO=1的所有位點(diǎn)
bcftools view  -i "AC>=2"  vep.vcf.gz  >  vep.vcf  #選出INFO中AC>2的所有位點(diǎn)
bcftools view  -i 'INHERITANCE[*] = "AR" || INHERITANCE[*] = "XR"'  ar.vcf  >  ar.AR.vcf  #選出遺傳方式是AR或XR的位點(diǎn)（需INFO字段中已有INHERITANCE注釋）
bcftools view  b37.vcf.bgz  X:31136335-33358725  >  DMD.vcf  #選出位于區(qū)域X:31136335-33358725的所有位點(diǎn)
bcftools view  -e 'CLINSIG~"Benign"'  Fun.vcf  >  Fun.exBenign.vcf  #選出除了INFO字段CLINSIG匹配"Benign"以外的所有位點(diǎn)
bcftools view  -i " MIN(FMT/DP)>500 && FORMAT/AD[1:1]/FORMAT/DP[1]>0.05 "  exTSG.vcf | sed '/^#/d' | less -S  #選出最小的depth>500而且，腫瘤樣品的VAF>0.05的所有位點(diǎn)

拆分snp和indel數(shù)據(jù)：

bcftools view  --no-version  -i '%TYPE=="snp"'  -Oz -o BL48384.snp.vcf.gz BL48384.Raw.fliter.vcf.gz
bcftools view  --no-version  -i '%TYPE=="indel"'  -Oz -o BL48384.indel.vcf.gz BL48384.Raw.fliter.vcf.gz
bcftools view  -v indels  hg19.vcf
bcftools view  -i 'TYPE="indel"'  hg19.vcf

-v/V, --types/--exclude-types <list>：select/exclude comma-separated list of variant types: snps,indels,mnps,ref,bnd,other

(8) stats命令用于統(tǒng)計(jì)VCF文件的基本信息

bcftools stats  view.vcf  >  view.stats
bcftools stats  -F  dmel.genome.fa  -s  -  sample.vcf.gz  >  sample.vcf.gz.stats

-F / --fasta-ref：faidx indexed reference sequence（參考基因組） file to determine INDEL context
-s：list of samples for sample stats, "-" to include all samples 統(tǒng)計(jì)的樣本列表临燃，在輸出結(jié)果中顯示所有的樣本名稱

(9) plot-vcfstats命令進(jìn)行可視化

plot-vcfstats  sample.vcf.gz.stats  -p  plots/sample.vcf.gz.stats

-p：指定輸出結(jié)果的目錄
（這個(gè)腳本位于bcftools安裝目錄的misc目錄下睛驳，依賴latex 生成pdf 文件）