Vcftools及vcf

一. VCF文件

Pasted image 20221026181450.png

1.1 简介

以“#”开头的注释部分：VCF的介绍信息，通常以##作为起始，其后一般接以FILTER，INFO，FORMAT等字样。

以##FILTER开头的行，表示注释VCF文件当中第7列中缩写词的说明，比如q10为Quality below 10；##INFO开头的行注释VCF第8列中的缩写字母说明，比如AF代表Allele Frequency也就是等位基因频率；##FILTER开头的行注释VCF第9列中的缩写字母说明；另外还有其他的一些信息，文件版本"fileformat=VCFv4.0"等等。

没有“#”开头的主体部分：每一行代表一个variant，各列之间用tab空白隔开，前面9列为固定列，第10列开始为样品信息列，可以无限多个

# 主体部分10列的范例
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  MP16
13      7       13:6    C       T       .       PASS    NS=9;AF=0.500   GT:DP:AD:GQ:GL  0/1:122:81,41:40:-70.65,0.00,-174.69

1.2 10列分别代表的意义

1.CHROM：染色体编号

2.POS：variant所在的left-most位置（1-base position，发生变异的位置的第一个碱基所在的位置）

3.ID：variant的ID，对应着dbSNP数据库中的ID。（SNP/INDEL的dbSNP编号通常以rs开头，一般只有人类基因组才有dbSNP编号，若没有，则默认使用‘.’）

4.REF：参考序列的Allele。（等位碱基，即参考序列该位置的碱基类型及碱基数量，必须是A,C,G,T,N且都大写）

5.ALT：variant的Allele，若有多个，则使用逗号分隔。（变异的碱基类型及碱基数量）这里的碱基类型和碱基数量，对于SNP来说是单个碱基类型的编号，而对于Indel来说是指碱基个数的添加或缺失，以及碱基类型的变化

6.QUAL：variants的质量。PPhred格式(Phred_scaled)的质量值，代表着此位点是纯合的概率，此值越大，则纯合概率越低，表示突变可能性越大。（表示变异碱基的可能性）
计算方法：Phred值 = -10 * log (1-p)， p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%。

7.FILTER：使用上一个QUAL值来进行过滤的话，是不够的。GATK能使用其它的方法来进行过滤，过滤结果中通过则该值为”PASS”;若variant不可靠，则该项不为”PASS”或为”.”。

8.INFO：variant的相关信息

9.FORMAT：variants的格式，例如GT:AD:DP:GQ:PL

10.SAMPLES：各个Sample的值，由BAM文件中的@RG下的SM标签所决定，这些值对应着第9列的各个格式，不同格式的值用冒号分开，每一个sample对应着1列；多个samples则对应着多列，这种情况下列的数多余10列。

1.3 FORMAT 及对应 SAMPLES值

FORMAT和最后一列（最后一列一般为样品名），两者和一起则为基因型信息，前者为格式，后者为对应的数据，如：

GT:AD:DP:GQ:PL    0/1:6,5:11:99:138,0,153
  1. GT（GenoType）：表示样品的基因型，通常用”/” or “|”分隔两个数字，“|”phase过也就是杂合的两个等位基因知道哪个等位基因来自哪条染色体

对于二倍体生物，GT值表示的是样本在这个位点所携带的两个等位基因。0代表参考基因组的碱基类型；1代表ALT碱基类型的第一个碱基（多个碱基用","分隔），2代表ALT第二个碱基，以此类推。0表示跟REF一样，1表示跟ALT一样，2表示有第二个ALT。

以二倍体生物为例，基因型由两条染色体上的allel构成。当我们知道每一个allel来自于具体哪条染色体时，这种genotype叫做Phased genotype, 用|连接，1|0和0|1代表两种不同的基因型；不清楚allel对应的染色体的时， genotype叫做unphased genotype, 用/连接，0/1和1/0这两种写法是等价的。目前高通量分析鉴定到的基因型，大多数都是unphased genotype。

  1. AD(Allele Depth)：sample中每一种allele（等位碱基）的reads覆盖度，在diploid（二倍体，或可指代多倍型）中则是用逗号分隔的两个值，前者对应REF基因，后者对应ALT基因型
  2. DP(Depth)：表示覆盖在这个位点的总reads数，也就是这个位点的测序深度（并不是指具体有多少个reads数量，而是大概满足一定质量值要求的reads数），是所支持的两个AD值（逗号前和逗号后）的加和
  3. PL（likelihood genotypes）：指定的三种基因型的质量值（provieds the likelihoods of the given genotypes），分别对应该位点的三个基因型0/0，0/1，1/1的没经过先验的标准化Phred-scaled似然值（L）
  1. GQ(Genotype Quality)：基因型的质量值。Phred格式（Phred_scaled）的质量值，表示在该位点该基因型存在的可能性

1.4 INFO 信息列

AC=1;AF=0.500;AN=2;BaseQRankSum=0.748;ClippingRankSum=0.000;DB;DP=34;ExcessHet=3.0103;FS=3.424;MLEAC=1;MLEAF=0.500;MQ=31.07;MQRankSum=-0.087;QD=11.87;ReadPosRankSum=-1.349;SOR=2.636
AC=2;AF=1.00;AN=2;DB;DP=14;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=31.60;QD=29.36;SOR=5.421

以 “TAG=Value”,并使用”;”分隔的形式。其中很多的注释信息在VCF文件的头部注释中给出。以下是这些TAG的解释：

  1. AC(Allele Count)：表示基因型为与variant一致的Allele的数目，Allele数目为1表示双倍体的样本在该位点只有1个等位基因发生了突变
  2. AF(Allele Frequency)：表示Allele的频率，AF值=AC值/AN值，Allele频率为0.5表示双倍体的样本在该位点只有50%的等位基因发生了突变
  3. AN(Allele Number)：表示Allele的总数目
  1. DP：样本在这个位置的reads覆盖度，是一些reads被过滤掉后的覆盖度（跟上面提到的DP类似）
  2. Dels：Fraction of Reads Containing Spanning Deletions。进行SNP和INDEL calling的结果中，有该TAG并且值为0表示该位点为SNP，没有则为INDEL。
  3. FS(FisherStrand)：使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值，值越小越好。
  1. HaplotypeScore：Consistency of the site with at most two segregating haplotypes.最多有2个分离的单倍型的一致性。
  2. InbreedingCoeff：Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hard-Weinberg expectation.与哈代温伯格的期望相比，近亲繁殖估计每个样品基因型的可能性。
  3. MLEAC：Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed.对于每个ALT等位基因，等位基因计数（不一定与AC相同）的最大似然期望（MLE），顺序与列出的顺序相同。
  4. MLEAF：Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT alle in the same order as listed.对于每个ALT等位基因，等位基因频率（不一定与AF相同）的最大似然期望（MLE），顺序与列出的顺序相同。
  5. MQ：RMS Mapping Quality 表示覆盖序列质量的均方值
  6. MQ0：Total Mapping Quality Zero Reads.总的Mapping 质量 零Reads 。
  7. MQRankSum：Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities.对Alt vs 参考片段映射质量的Wilcoxon秩和检验的z 分数。
  1. BaseQRankSum：Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities.来自Wilcoxon的Z分数 Alt与Ref基本质量的秩和测试
  1. ClippingRankSum：Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases.Z 得分来自 Wilcoxon 的 Alt 与 Ref 硬剪切基数的秩和检验
  2. ExcessHet：过量Het.Phred-scaled p-value for exact test of excess heterozygosity.用于精确检验过量杂合度的Phred标度p值。
  1. QD：Variant Confidence/Quality by Depth.变异置信度/深度质量。
  1. RPA：Number of times tandem repeat unit is repeated, for each allele (including reference).对于每个等位基因（包括参考），大量的串联重复序列单位被重复。
  2. RU：Tandem repeat unit (bases).串联重复序列单元（基础）。
  3. ReadPosRankSum：Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias.对Alt vs 参考片段位置偏差的Wilcoxon秩和检验的z 分数。
  1. SOR：Symmetric Odds Ratio of 2x2 contingency table to detect strand bias. 2x2 列联表的对称比值比，用于检测链偏置。
  1. STR：Variant is a short tandem repeat.Variant是一个短的串联重复。

二. Vcftools用法

2.1 过滤变异类型

vcf文件中可能会同时包含snp以及indel两种变异类型，vcftools可以很快的将两者进行分离。

2.1.1 过滤掉indel，只保留snp，用到的命令选项：–remove-indels。

vcftools --remove-indels --recode --recode-INFO-all --vcf raw.vcf --stdout >raw.snp.vcf

2.1.2 过滤掉snp，只保留indel，用到的命令选项：–keep-only-indels。

vcftools --keep-only-indels --recode --recode-INFO-all --vcf raw.vcf --stdout >raw.indel.vcf

这样，就可以分别得到只包含snp和indel的vcf文件。

2.2 筛选指定位置变异位点

vcftools还可以挑选出基因组上某些区域的变异信息。

vcftools --vcf Variants.snp.unknown_multianno.vcf --chr A03 --from-bp 577700 --to-bp 607700 --out out_prefix --recode --recode-INFO-all

这里解释一下各个参数：
–vcf：后面跟的是vcf文件
–chr：后面跟筛选区域所在的染色体
–form-bp：后跟筛选区域的起始位置
–to-bp：后跟筛选区域的终止位置
–out：输出文件的前缀
–recode：没有此参数则不会输出

2.3 过滤指定缺失率的变异位点

vcf 文件中很多snp在某些样品中是缺失的，也就是基因型为 “./.” 。如果缺失率较高，这种snp位点在很多分析中是不能用的，需要去掉。这里用到的选项是 --max-missing。

vcftools --vcf snp.vcf --recode --recode-INFO-all --stdout --max-missing 1 > snp.new.vcf
--max-missing 后跟的值为 0-1 ，1代表不允许缺失，0代表允许全部缺失。

2.4 计算snp缺失率

vcftools中有两个参数可以计算vcf文件中snp的缺失率。
分别是：
–missing-indv：生成一个文件，报告每个样品的缺失情况，该文件的后缀为“.imiss”。
–missing-site：生成一个文件，报告每个snp位点的缺失情况，该文件的后缀为“.lmiss”。

2.4.1 –missing-site

vcftools --vcf snp.vcf. --missing-site

运行以上命令后会在当前目录生成一个 out.lmiss 文件，其格式如下：

CHR POS N_DATA N_GENOTYPE_FILTERED N_MISS F_MISS  
chr01 194921 988 0 368 0.37247  
chr01 384714 988 0 204 0.206478  
chr01 384719 988 0 202 0.204453  
chr01 518438 988 0 488 0.493927  
chr01 518473 988 0 452 0.45749  
chr01 518579 988 0 418 0.423077  
chr01 518635 988 0 428 0.433198  
chr01 680786 988 0 346 0.350202  
chr01 680834 988 0 412 0.417004  
# 前两列为snp所在位置，第三列为等位基因总数，第5列为缺失的总数，最后一列为缺失率。

2.4.2 –missing-indv

vcftools --vcf snp.vcf. --missing-indv

运行以上命令后会在当前目录生成一个 out.imiss 文件，其格式如下：

INDV N_DATA N_GENOTYPES_FILTERED N_MISS F_MISS  
1 8747 0 3632 0.415228  
10 8747 0 1264 0.144507  
102 8747 0 2016 0.230479  
105 8747 0 6322 0.722762  
106 8747 0 2365 0.270378  
107 8747 0 4376 0.500286  
108 8747 0 5682 0.649594  
109 8747 0 1877 0.214588  
11 8747 0 1039 0.118784  
# 第一列为样品名称，第二列为总的snp数，第4列为缺失的总数，最后一列为缺失率。

2.5 随机抽取指定个样品

vcftools可以随机抽取指定个样品的vcf文件，用到的选项为 --max-indv ，指定要从vcf文件中随机抽取指定个样品。

# 随机抽取5个样品，执行以下代码：
vcftools --vcf snp.vcf --max-indv 5 --remove-indels --recode --out outfilename