软件12 —— bcftools

一、 基本介绍

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 文件）