一、 基本介绍
VCF格式（Variant Call Format）是存储变异位点的标准格式，用于记录variants（SNP / InDel）。BCF是VCF的二进制文件。bcftools可以实现SNP calling。
二、 背景知识
(1) 变异和突变
变异，指的是实际测序数据与国际规定的参考基因组之间的区别。很多变异其实只是造成人类多样性的原因。突变，指的是那些与疾病相关的变异。举个例子：ENSEMBL等规定的人类参考基因组文件某位置是AAAAA，然后一个人实际测序得到的序列为AGCAA，那么相比于参考基因组，这个人就有2个变异位点。对于第2个位置，如果查看所有已知的测序，绝大部分人都是G，说明是参考基因组出现了问题，这个变异就不能称作突变。对于第3个位置，如果查看所有已知的测序，绝大部分人都是A，而恰好有一个人不是A，但他是个患者，那么这个变异就是突变了。
(2) 变异类型
- SNP（single nucleotide polymorphism）：单核苷酸多态性。个体间基因组DNA序列同一位置单个核苷酸变异（替换、插入或缺失）所引起的多态性。在人类基因组中SNP分布普遍并且密度较大，总数超过10^7， 平均每300bp（也有说1kbp）就有一个SNP。或称单核苷酸位点变异SNV。
- INDEL（insertion-deletion）：插入和缺失。基因组上小片段（>50bp）的插入或缺失。
- CNV（copy number variation）：基因组拷贝数变异。基因组中大片段的DNA形成非正常的拷贝数量。比如一个基因在染色体的一条染色单体上的数目为1，但是在染色体复制过程中，复制结束后该基因在染色单体数目由1变成了2或者n。它发生的频率远远高于染色体结构变异，并且整个基因组中覆盖的核苷酸总数大大超过SNP的总数。
- SV（structure variation）：结构变异。染色体大片段的插入与缺失，染色体内部的某区域发生翻转颠换，两条染色体之间发生重组。
(3) SNP
一般情况下只分析SNP，其它类型的变异分析有难度或不准确。来自两个不同个体的DNA片段AAGCCTA和AAGCTTA为等位基因。几乎所有常见的SNP位点只有两个等位基因。在人体中，SNP的发生机率大约是0.1%，也就是每1000个碱基对就可能有一个SNP（密度高）。对疾病发生和药物治疗有重大影响的SNP，估计只占数以百万计SNP的很小一部分。SNP位点的分布是不均匀的，在非转录序列比在转录序列更常见。编码区的单核苷酸多态性——编码 SNP（coding SNP，cSNP）也有同义和非同义两种类型，非同义SNP会改变蛋白质的氨基酸序列。基因非编码区、基因间隔区的SNP仍然可能影响转录因子结合、剪接等过程。从演化的观点来看，SNP具有相当程度的稳定性，即使经过代代相传，SNP所引起的改变却不大，因此可用以研究族群演化。
(4) vcf格式
vcf格式（Variant Call Format）是存储变异位点的标准格式，用于记录variants（SNP / InDel）。BCF是VCF的二进制文件。
- 以#开头的注释部分：
##fileformat：VCF格式版本号。
##reference & contig：使用的参考基因组信息及参考基因组contig信息。
##INFO行：是碱基位点的注释。每一行必须的四个标签是：ID、Number、Type、Description。
- 没有#开头的主体部分：
包含10列数据，每一行代表一个variant的信息。
主体部分10列的范例：CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO、FORAMT、SAMPLE（前8列必须要有）。
例如：
chrM（CHROM染色体）
150（POS变异的第一个位置，1-based coordinate system）
.（ID变异位点名称，在dbSNP数据库中的ID以rs开头 ，一般只有人类基因组才有dbSNP编号，如果没有则为点）
T（REF参考序列该位置碱基类型和数目）
C（ALT该位置变异的碱基类型和数目，点代表缺失，多个用逗号分隔）
7766.77（QUAL变异的质量值，Q=-10lgP，值越大是变异的可能性越大）
PASS（FILTER是否要被过滤掉，为PASS表示可能是变异，点代表没有进行任何过滤）
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的相关信息。）
GT:AD:DP:GQ:PL（FORAMT：variants的格式）
1/1:0,175:175:99:7795,531,0（一个SAMPLE为1列，总列数可以多于10，每列分别对应第9列的各个格式，由bam文件中@RG的SM标签决定）
第8列INFO：
variant的相关信息。
AC（Allele Count）变异的等位基因数目
AF（Allel Frequency）等位基因频率
AN（Allel Number）等位基因总数目
DP，reads覆盖度
FS，Fishers精确检验的p值
AC（Allele Count）变异的等位基因数目。AF（Allele Frequency）等位基因频率，AN（Allele Number）等位基因总数目。（0/1：AC=1，AF=0.5，AN=2）
DP，是一些reads被过滤掉后的覆盖度。
DP4，高质量测序碱基，在ref或alt前后。
Dels，有这个tag且为0时表示该位点是SNV，没有就是InDel。
FS，使用Fisher精确检验来检测strand bias而得到的Fhred格式的p值。该值越小越好。一般进行filter的时候，可以设置 FS < 10～20。
第9列FORMAT：
variants的格式。
GT（genotype），0代表样本中ref的allel，1代表样本variant的allel，2表示有第二个variant的allel；0/0表示样品中该位点为纯合位点，和REF的碱基类型一致；0/1表示sample中该位点为杂合突变，AC=1，AF=0.5，AN=2；1/1表示为变异纯合子，AC=2，AF=1，AN=2。
AD（Allele Depth）为sample中每一种allele（等位碱基）的reads覆盖度。
DP（Depth）为sample中该位点的覆盖度。
GQ（Genotype Quality）基因型的质量值，基因型存在的概率。
PL（likelihood genotypes）指定的三种基因型的质量值（0/0，0/1，1/1），对应的值越大，表示这种基因型的可能性越小。
GT（genotype），0代表样本中ref的allele，1代表样本variant的allele，2表示有第二个variant的allele（几乎所有常见的SNP位点只有两个等位基因）。0/0表示样品中该位点为纯合位点，和REF的碱基类型一致；0/1表示sample中该位点为杂合突变，AC=1，AF=0.5，AN=2；1/1表示为变异纯合子，AC=2，AF=1，AN=2。
AD（Allele Depth）对应两个以逗号隔开的值，这两个值分别表示覆盖到REF和ALT碱基的reads数，相当于支持REF和支持ALT的测序深度。
DP（Depth）覆盖到这个位点的总的reads数量，相当于这个位点的深度。
PL（likelihood genotypes）对应3个以逗号隔开的值，指定的三种基因型（0/0，0/1，1/1）没经过先验的标准化Phred-scaled似然值。对应的值越大，表示这种基因型的可能性越小。
GQ（Genotype Quality）最可能的基因型的质量值。
例如：
chr1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26
GT=0/1，也就是说这个位点的基因型是C/T。AD=1,3，也就是说支持REF的read有一条，支持ALT的有3条。DP=4，也就是说只有4条reads支持这个地方的变异，cover到这个位点的reads数太少。GQ=25.92，质量值并不算太高。在PL里，这个位点基因型的不确定性就表现的更突出了，0/1的PL值为0，虽然支持0/1的概率很高；但是1/1的PL值只有26，也就是说还有10^(-2.6)=0.25%的可能性是1/1；但几乎不可能是0/0，因为支持0/0的概率只有10^(-10.3)=5*10^-11。
- 1-based coordinate system：序列的第一个碱基设为数字1，如SAM, VCF, GFF, wiggle格式
- 0-based coordinate system：序列的第一个碱基设为数字0，如BAM, BCFv2, BED, PSL格式
三、 用法和参数
(1) SNP calling
mpileup命令：得到染色体上每个碱基的比对情况的汇总（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命令：执行SNP calling
bcftools call -vmO z -o sample.vcf.gz sample.bcf
-v：只输出变异位点的信息，如果一个位点不是snp/indel则不会输出
两种calling算法：-c参数对应consensus-caller算法，-m参数对应multiallelic-caller算法，后者更适合多种allel和罕见变异的calling
(2) tabix
tabix -p vcf sample.vcf.gz
输入为压缩文件vcf.gz，生成的索引文件为sample.vcf.gz.tbi
(3) index对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通过表达式来指定输出格式
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：代表新的一行
bcftools query -r '19:400300-400800' -f '%CHROM\t%POS\t%REF\t%ALT\n' hg19.vcf.gz | head -3
-r：从特定区域提取varients信息，格式 chr|chr:pos|chr:from-to|chr:from-[,…]
bcftools query -t ^'19:400300-400800' -f '%CHROM\t%POS\t%REF\t%ALT\n'
排除特定区域
(5) sort按照染色体位置进行排序
bcftools sort view.vcf.gz -o sort.view.vcf
(6) filter过滤不可靠位点
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：将过滤掉的位点用字符串注释
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这列信息，过滤掉的信息为FLTER，保留的为PASS
-i：筛选条件，筛选出QUAL大于20、DP大于4、MQ大于30或者GT等于0/0的位点
-Ov：输出文件为未压缩vcf格式
(7) view命令用于VCF和BCF格式的转换
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：想要保留的样本信息，多个样本用逗号分隔；如果样本名称添加了^前缀，代表去除这些样本，比如-s ^NA00001,NA00002
-k：表示筛选已知的突变位点，即ID那一列值不是.的突变位点
bcftools view -i 'SAO=1' b37.vcf > b37.germline.vcf #选出INFO中SAO=1的所有位点
bcftools view -i "AC>=2" vep.vcf.gz > vep.vcf #选出INFO中AC>2的所有位点
bcftools view -i 'INHERITANCE[*] = "AR" || INHERITANCE[*] = "XR"' ar.vcf > ar.AR.vcf #选出遗传方式是AR或XR的位点（需INFO字段中已有INHERITANCE注释）
bcftools view b37.vcf.bgz X:31136335-33358725 > DMD.vcf #选出位于区域X:31136335-33358725的所有位点
bcftools view -e 'CLINSIG~"Benign"' Fun.vcf > Fun.exBenign.vcf #选出除了INFO字段CLINSIG匹配"Benign"以外的所有位点
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的所有位点
拆分snp和indel数据：
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命令用于统计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 统计的样本列表，在输出结果中显示所有的样本名称
(9) plot-vcfstats命令进行可视化
plot-vcfstats sample.vcf.gz.stats -p plots/sample.vcf.gz.stats
-p：指定输出结果的目录
（这个脚本位于bcftools安装目录的misc目录下，依赖latex 生成pdf 文件）
网友评论