vcf文件

生信喵

一、背景
       VCF 是生物信息分析中非常重要的一种格式。主要用来描述基因组突变的信息，无论是检测出来的 SNP，indel，cnv，还是 SV，都可以存储格式都为 vcf 格式。从比对生成的 bam 文件中，将潜在变异信息筛选出来，就是 vcf 格式。vcf 是一种列表格式，里面包含很多的内容。需要掌握每一列的信息，并能使用相对应的软件对 vcf 进行处理是。处理 VCF 格式软件主要包括 bcftools，vcftools，gatk，python pyvcf，plink 等。

二、vcf 文件格式介绍
2.1 vcf 简介
       VCF 是 Variant Call Format 的简称，是一种定义的专门用于存储基因序列突变信息的文本格式。在生物信息分析中会大量用到 VCF 格式。例如基因组中的单碱基突变,SNP，插入/缺失INDEL, 拷贝数变异 CNV，和结构变异 SV 等，都是利用 VCF 格式来存储的。vcf 是一种文本格式，可以直接查看。将其存储为二进制格式就是 BCF，二进制格式节省更多存储，vcf 与bcf 的关系类似 sam 与 bam 的关系。需要特别之处的是，不同软件产生的 vcf 会有很大的不同，有时候同样的操作命令在不同的 vcf 中会出错。当前 vcf 的版本为 4.3，可以参考下面的帮助文档，格式说明：

1. https://samtools.github.io/hts-specs/

复制代码

2.2 vcf 文件格式
       vcf 是一种表格是格式，主要分为三部分，第一部分为双井号注释的部分，为文件头信息，主要介绍文件内容以及 INFO 部分的详细解释；
第二部分单井号注释，为表头信息，基本内容分为 8 列，对于多样品可以继续添加列。前 8列信息分别为：

1. 1.CHROM [chromosome]： 染色体名称，
   
2. 2.POS [position]： 参考基因组突变碱基位置，如果是 INDEL，位置是 INDEL 的第一个碱基位
   
3. 置。
   
4. 3.ID [identifier]： 突变的名称，
   
5. 4.REF [reference base(s)]：参考染色体的碱基
   
6. 5.ALT [alternate base(s)]： 与参考序列比较，发生突变的碱基，
   
7. 6.QUAL [quality]： Phred 标准下的质量值
   
8. 7.FILTER [filter status]：使用其它的方法进行过滤后得到的过滤结果
   
9. 8.INFO

复制代码

      vcf 中可以保存多个样品的信息，当文件中包含多个样品时，就会出现“FORMAT” 一列，用于提示后续不同样品中展示的信息。由于样品多，不同样品不能展示全部 INFO 的信息，通常 FORMAT 只展示少部分信息，例如“Genotype”一个信息。接下来是样品名。 每个样品在后面增加一列即可，展示全部 INFO 的信息，通常FORMAT 中及介绍的内容，这样就能构成一个很大的矩阵，可以用于统计检验。

2.3 INFO 信息
       vcf 中的 INFO 关键字非常多，而且每个软件生成的 vcf 文件都可以单独自定义关键字。都是以 “TAG=Value”,并使用”;”分隔的形式。其中很多的 TAG 含义在 VCF 文件的头部注释信息##INFO 中已给出。这些关键字信息包含了非常多的内容，描述了每一个突变详细的信息。例如突变的类型，SNP 还是 SV，如果是 SNP 杂合还是纯合，如果是 SV，具体哪种类型，发生变化的长度是多少，有多少条 reads 支持等信息。这些信息根据不同的需求可以从中提取。在 vcf 文件格式的详细描述中会有介绍，无需掌握全部内容，只需要了解一些常见的关键信息即可。
       下面介绍一些常见关键字内容，更多内容可以查看表头信息，或者查看 vcf 详细文档。
       GT：GT 为基因型（genotype）的简写，是 vcf 中非常重要的关键字。在做变异检测之后通常还要做一步 genotype 就是为了得到这个信息。两个数字中间用’/"分开，这两个数字表示双倍体的 sample 的基因型，
       0 表示样品的基因型与 ref 的 allele 相同；
       1 表示样品中基因型与 alt variant 的相同；
       2 表示有第二个 variant 的 allele（和 ALT 的第二种碱基相同）
       所以，我们就可以根据 GT 关键字判断出样品的基因型。假设这个突变中 REF 为 A 碱基， ALT 为 T 碱基。如果 GT 为：

1. 0/0 ：表示 sample 中该位点为纯合位点，和 REF 的碱基类型一致，那么就是 AA。
   
2. 0/1：表示 sample 中该位点为杂合突变，有 REF 和 ALT 两个基因型，那么就是 AT；如果是 GT 是
   
3. 1/0，则表示 TA。
   
4. 1/1：表示 sample 中该位点为纯合突变，总体突变类型和 ALT 碱基类型一致，为 TT。
   
5. 1/2：表示 sample 中该位点为杂合突变，有 ALT1 和 ALT2 两个基因型，可能是 TC 或者 TG。

复制代码

      AD：Allele Depth：为 sample 中每一种 allele（等位碱基）的 reads 覆盖度，在 diploid（二倍体，或可指代多倍型）中则是用逗号分隔的两个值，前者对应 REF 基因，后者对应 ALT基因型；
       DP：Depth：为 sample 中该位点的覆盖度，是所支持的两个 AD 值（逗号前和逗号后）的加和，支持数越高，结果越可信，通常可以用于 DP 进行突变结果过滤，例如将小于 5 条 reads支持的过滤掉。
       GQ：基因型的质量值(Genotype Quality)。Phred 格式(Phred_scaled)的质量值，表示在该位点该基因型存在的可能性；该值越高，则 Genotype 的可能性越 大；计算方法：Phred 值 = -10 * log (1-p) p 为基因型存在的概率。和 DP 类似，也可以用于突变结果过滤。
       PL：指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为 1。和之前不一致，该值越大，表明为该种基因型的可能性越小。 Phred 值 = -10 * log (p) p 为基因型存在的概率。
       AC：（Allele Count）：Alt Allele 的数目（即当前位点等位基因数量）。
       AN：（Allele Number）：表示 alt Allele 的总数目。
       AF：（Allele Frequency）：表示 alt Allele 的频率，AF 值=AC 值/AN 值。
       FS：FS 是一个通过 Fisher 检验的 p-value 转换而来的值，它要描述的是测序或者比对时对于只含有变异的 read 以及只含有参考序列碱基的 read 是否存在着明显的正负链特异性（Strand bias，或者说是差异性）。这个差异反应了测序过程不够随机，或者是比对算法在基因组的某些区域存在一定的选择偏向。如果测序过程是随机的，比对是没问题的，那么不管 read 是否含有变异，以及是否来自基因组的正链或者负链，只要是真实的它们就都应该是比较均匀的，也就是说，不会出现链特异的比对结果，FS 应该接近于零。
       MLEAC：AC（Allele counts）的极大似然期望值。
       MLEAF：AF（Allele Frequency）的极大似然期望值。
       MQ：（Mapping quality）：比对质量值。

2.4 vcf 文件中如何描述 SV
       在 vcf 文件中，SV 通常可以通过 SVTYPE 关键字进行描述，然后用 SVLEN 关键字描述具体发生 SV 的长度。由于 SV 包含多种类型，SVTYPE 关键字也可以分成以下几种：

1. DEL ：缺失，Deletion relative to the reference
   
2. INS： 插入，Insertion of novel sequence relative to the reference
   
3. DUP ：倍增，Region of elevated copy number relative to the reference
   
4. INV ：翻转，Inversion of reference sequence
   
5. CNV ：拷贝数变异，Copy number variable region (may be both deletion and duplication)
   
6. BND ：断开，Breakend
   

复制代码

      CNV 拷贝数变化属于插入和缺失范畴，有些软件没有被单独列为一个类目，被包含在倍增DUP，缺失 DEL 以及插入 INS 中。

2.5 注意事项
1、不同版本的 vcf 格式文件有所差别，在使用过程中需要注意 vcf 的版本；
2、不同软件生成的 vcf 有很大的差别，尤其是对于 SV 的描述方法，例如 gatk，freebayes，sniffles 等，有很大的不同。
3、不同软件生成的 vcf 文件，INFO 部分会有很大的不同，在使用过程中要根据具体的内容修改代码。

三、利用 bcftools 处理 vcf 文件
       处理 VCF 格式软件：bcftools，vcftools，gatk，python pyvcf，plink 等。最常用的为 bcftools。这里我们分别介绍 bcftools 与 vcftools。bcfools 可以直接处理 vcf 文件，也可以处理二进制bcf 文件。
      
       bcftools 文件处理
3.1 软件安装
       软件说明文档：

1. http://www.htslib.org/doc/bcftools.html

复制代码

      bcftools 是专门用来处理 vcf 以及 bcf 文件格式的工具，使用起来类似与 samtools。软件的安装和使用非常简单。

1. conda install -y bcftools

复制代码

生信喵

3.2 软件介绍
      主要分为三大功能类。
      Indexing 建立索引；
      VCF/BCF manipulation ：vcf 和 bcf 文件操作；
      VCF/BCF analysis ：vcf 和 bcf 文件分析；
1、Indexing

1. index，就是为文件建立索引；

复制代码

2、manipulation 大项主要是对文件的一些操作，包括查看，注释，合并，转换，排序等。

1. annotate 注释，用于增加或者移除注释信息；
   
2. concat 连接同一个样品或者同一个数据集的 vcf 或者 bcf 文件；
   
3. convert 格式转换 VCF 与 BCF 之间转换，或者转换为其他格式；
   
4. isec 取交集，多个 vcf 文件文件可以取共有的突变；
   
5. merge 和 samtools 的 merge 类似，合并文件，
   
6. norm normalize 标准化 indels
   
7. plugin 用户自定义功能
   
8. query 将 vcf 或 bcf 转换为用户自定义格式；
   
9. reheader 修改 header 文件，修改样品名等；
   
10. sort 对文件排序；
    
11. view 查看文件
    

复制代码

3、VCF/BCF analysis

1. call Call SNP 和 Indel；
   
2. consensus 将所有变异位点合并成一个一致性序列
   
3. cnv Copy Number Variation caller，检测拷贝数变异
   
4. csq 检测变异的一致性序列；
   
5. filter 过滤 vcf
   
6. gtcheck 检查样本一致性，检测样品交换和污染
   
7. mpileup 和 samtools 的 mpileup 类似；
   
8. roh 识别 autozygosity
   
9. stats 进行统计

复制代码

4、常用选项参数

1. output options:
   
2. --no-version
   
3. -o, --output <file> ：输出文件名，最好与文件类型保持一致，例如 vcf，vcf.gz，bcf，
   
4. bcf.gz
   
5. -O, --output-type <b|u|z|v>：输出文件类型，b,u,z,v 分别对应 bcf 与 vcf，以及是否压缩
   
6. --threads <int> ：线程数
   
7. -e, --exclude <expr> ：排除
   
8. -i, --include <expr> ：包括
   
9. -r, --regions <region> ：区域，给定表达式，染色体：起始位置-终止位置
   
10. -R, --regions-file <file> ：区域，bed 文件
    
11. -s, --samples <list> ：样品名，多个样品之间用逗号分隔
    
12. -S, --samples-file <file> ：样品名文件
    
13. -t, --targets <region> ：目标，染色体名字，例如 chr1，多个染色体之间用逗号分隔
    
14. -T, --targets-file <file> ：目标，写到文件中
    
15. -l, --compression-level [0-9] ：压缩等级，数字越大，压缩率越高，压缩时间越长

复制代码

5、表达式

1. http://www.htslib.org/doc/bcftools.html#expressions

复制代码

3.3 bcftools 使用案例
1、格式转换
      格式转换是将 vcf 与 bcf 之间进行格式转换，并同时进行压缩，bcf 为二进制格式，无法使用 less 等命令直接查看，但更加节约存储。如果是 bcf 格式文件，可以使用 view 功能进行查看。vcf 或者 bcf 必须使用 bgzip 压缩。

1. #格式转换
   
2. bcftools view file.vcf -O b -o file.bcf.gz
   

复制代码

2、建立索引
      对于行数较多的数据，需要建立索引。这样便于快速检索。索引需要基于排序后的结果。后面很多应用都需要使用到索引。排序可以使用 sort 功能，索引使用 index 功能。index 命令用于对 VCF 文件建立索引，要求输入的 VCF 文件必须是使用 bgzip 压缩之后的文件，支持.csi和.tbi 两种索引，默认情况下建立的索引是.csi 格式。

1. #排序
   
2. bcftools sort file.bcf.gz
   
3. #建立索引
   
4. bcftools index file.bcf.gz

复制代码

3、统计绘图
      如果想知道 vcf 中包含多少中突变，以及每种突变对应的数据，就需要对文件进行统计，stats选项可以用于 vcf 文件的统计。可以直接统计突变个数、突变类型的个数、转换颠换个数、测序深度、Indel 长度。统计完成之后可以使用 plot-vcfstats 进行可视化绘图。

1. #数据统计与绘图
   
2. bcftools stats file.bcf.gz >view.stats
   
3. plot-vcfstats view.stats -p output
   

复制代码

4、查看固定区域
      建立索引之后，就可以快速检索目标区域，例如得到结果之后，想快速查看某一区域，或者某个基因的情况，可以直接填写目标区域，格式为“染色体：起始位点-终止位点”。如果一次要查看多个区域，则可以添加一个 bed 文件。

1. #查看固定区域
   
2. bcftools view file.bcf.gz
   
3. bcftools view all.bcf.gz 12:200000-300000
   
4. #选出位于 bed 文件中的所有区域的突变位点
   
5. bcftools view -R Target.bed file.bcf.gz >filt.onTarget.vcf

复制代码

5、拆分不同类型

1. #提取 SNP
   
2. bcftools view -v snps chr22.vcf >chr22.snp.vcf
   
3. #提取 InDel
   
4. bcftools view -v indels chr22.vcf >chr22.indel.vcf
   
5. #提取 SV
   
6. bcftools view -v other chr22.vcf >chr22.sv.vcf

复制代码

6、提取某一条染色体

1. #提取 21 号染色体
   
2. bcftools view -r 21
   
3. ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz >chr21.vc
   
4. f
   
5. #提取 22 号染色体
   
6. bcftools view -r 22
   
7. ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz >chr22.vc
   
8. f

复制代码

7、concat 合并
      concat 命令用于合并 vcf 文件，为“横向合并“，是同一样品内部的合并，例如将同一样品不同染色体的突变信息进行合并，或者将同一样品的 SNP 结果与 InDel 结果进行合并。合并之前每个 VCF 文件必须是排序之后的，如果包含多个样本的信息，样本的顺序也必须一致。

1. #合并 21 和 22 号染色体的结果
   
2. bcftools concat chr21.vcf chr22.vcf -o merge.bcf.gz -O b
   
3. #合并 SNP 与 InDel
   
4. bcftools concat chr22.snp.vcf chr22.indel.vcf -o chr22.merge.vcf

复制代码

8、提取固定样品

1. #筛选样品 HG00100
   
2. bcftools view -s HG00100
   
3. ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz -O
   
4. b >HG00100.vcf
   
5. #筛选样品 HG00141
   
6. bcftools view -s HG00141
   
7. ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz -O
   
8. b >HG00141.vcf

复制代码

9. merge 合并
      vcf 中不仅可以包含单个样品，也可以同时包含多个样品的信息，只需要将多个样品的 vcf合并即可。可以使用 merge 功能进行合并，与 concat“横向合并”不同，merge 合并为“纵向合并”，是合并不同样品到同一个 vcf 文件中。注意合并之前需要对每个样品创建索引。

1. #合并样品 HG00100 与 HG00141
   
2. bcftools merge a.vcf.gz b.vcf.gz -o merge.vcf
   

复制代码

     注意：输入文件必须是经过 bgzip 压缩的文件， 而且还需要有.tbi 的索引。
10、集合运算
      isec 用于在多个 VCF 文件之间取交集，差集，并集等操作，经典的应用场景是对多种软件的SNP calling 结果进行 venn 分析。用法如下

1. #isec 取交集
   
2. bcftools isec a.vcf.gz b.vcf.gz -p dir

复制代码

     默认参数就是取交集，更多高级用法请参考帮助文档。
11、query 功能
      vcf 里面包含的信息非常多，比较混乱，如果只想从中筛选出需要的内容，例如只需要Genotype 信息，可以使用 bcftools 的 query 功能实现。query 中最重要的就是表达式的写法。

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

复制代码

     注意格式的熟悉，每个关键字前面使用%，“\t”或者“\n”代表制表符与换行符。
12、过滤
      变异检测的策略一般是先找全，然后在找准。也就是软件首先输出尽可能多的结果，保存到vcf 文件中，然后在采取不同的标准对 vcf 进行过滤。过滤可以采取很多的标准，一般包括测序深度，打分制，碱基质量值，先验概率等。可以使用 bcftools 进行过滤。bcftools 的 filter功能其实与 query，view 都类似，可以进行多种模式的过滤，关键是要掌握其表达式EXPRESSIONS 的写法。

1. #过滤掉等位频率>0.3 同时覆盖深度小于 10 的突变位点
   
2. bcftools filter -e "INFO/AF[0] > 0.3 & FORMAT/DP < 30" A1.bcf.gz
   
3. bcftools view -i '%TYPE!="snp" & MAF[0]<0.01'
   
4. ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.g

复制代码

13、注释
      vcf 的注释主要是将突变位点定位到基因组上，确定突变发生在哪个基因，因为不同的突变发生位置，会对基因产生不同的影响，例如同义突变，错误突变或者无义突变等。另外一种注释就是与已知突变位点进行比较，定位到已知的 rs number 号上面。
      annotate 命令有两个用途：

1. #注释 VCF 文件，用法如下
   
2. bcftools annotate -a /ifs1/Database/gatk/hg38/dbsnp_146.hg38.vcf.gz -c ID,QUAL
   
3. ,+TAG file.vcf -o annotate.vcf

复制代码

     -a 参数指定注释用的数据库文件，格式可以是 VCF, BED, 或者是\t 分隔的自定义文件。在\t 分隔的自定义文件中，必须包含 CHROM, POS 字段；
      -c 参数指定将数据库的哪些信息添加到输出文件中。

1. #编辑 VCF 文件，比如去除其中的某些注释信息，或者去除某些样本
   
2. bcftools annotate -x ID,INFO/DP,FORMAT/DP view.vcf -o removed.id.vcf

复制代码

     -x 参数表示去除 VCF 文件中的注释信息，可以是其中的某一列，比如 ID, 也可以是某些字段，比如 INFO/DP，多个字段的信息用逗号分隔；去除之后，这些信息所在的列并不会去除，而是用.填充。