看懂变异记录结果文件（VCF）
VCF文件介绍:
做过DNA重测序，群体遗传进化，BSA，GWAS等项目的人都会遇到VCF文件，这个文件记录了所有样品全基因组中所有位置变异（主要包括SNP和InDel）信息。后续几乎所有的分析内容都是基于此文件，比如进化树分析、群体结构分析、PCA分析、GWAS关联分析等等，如果不懂VCF文件就无法进行后续分析。

因此了解VCF文件格式及其记录结果的意义非常重要。VCF文件其实是文本文件，可以用Windows当中文本编辑器软件打开，比如editplus等。由于VCF文件往往很大（通常超过1G），在Windows系统下直接打开会消耗大量内存进而造成卡死的现象。如果想顺利打开的话，可以使用cmd命令行more命令（more E:\vcf\merge.final.vcf）打开查看或使用pilotedit（http://www.pilotedit.com/），这里建议使用pilotedit。

官方说明：http://www.internationalgenome.org/wiki/Analysis/vcf4.0/

下面是一个典型VCF文件的示例（部分）：

整体说明信息（Meta-information lines）

VCF文件的开头是整体注释信息，通常以##作为起始，其后一般接以FILTER，INFO，FORMAT等字样。

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

VCF各列意义说明

各列之间用tab空白隔开；前面9列为固定列，第10列开始为样品信息列，可以无限多个；

具体说明如下

1. CHROM 记录染色体编号

2. POS 记录染色体位置信息

3. ID SNP/INDEL的dbSNP编号通常以rs开头，一般只有人类基因组才有dbSNP编号

4. REF 参考基因组碱基类型，必须是A,C,G,T,N且都大写。

5. ALT 变异碱基类型，必须是A,C,G,T,N，. 且都大写，多个用逗号分割。"."表示这个地方没有reads覆盖为缺失。* 号参考：https://www.omicsclass.com/question/2230

6. QUAL 变异信息的检测质量值，越高越可靠。

7. FILTER 标记过滤结果的列，通常我们把VCF文件中的变异信息进行质控，过滤掉低质量的变异位点，如果该位点通过过滤标准那么我们可以在该列标记为"PASS",说明该列质量值高。标记完之后我们就可以用其他工具，把标记为"PASS"的列给筛选出来，这样方便后续分析。如果没有应用缺失值"."代替。

8. INFO 为附加信息列，一般以
   =；形式添加额外的注释信息列，常见的如DP=18 表示该位点测序深度为18X；AF=0.1表示等位基因频率为0.1；

BaseQRankSum：比较支持变异的碱基和支持参考基因组的碱基的质量，负值表示支持变异的碱基质量值不及支持参考基因组的。

DP: read depth。样本在这个位置的reads覆盖度。

ExcessHet：检测样本的相关性，与InbreedingCoeff相似，值越大越可能是错误。

FS：使用F检验来检验测序是否存在链偏好性（？）。链偏好性可能会导致变异等位基因检测出现错误。输出值Phred-scaled p-value，值越大越可能出现链偏好性。

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

MLEAF：Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed

MQ：表示覆盖序列质量的均方值RMS Mapping Quality

MQRankSum：比较支持变异的序列和支持参考基因组的序列的质量，负值表示支持变异的碱基质量值不及支持参考基因组的，只针对杂合。正值则相反，支持变异的质量值好于参考基因组的。0表示两者无明显差异。实际应用中一般过滤掉较小的负值。

QD：通过深度来评估一个变异的可信度。

ReadPosRankSum：检测变异位点是否有位置偏好性（是否存在于序列末端，此时往往容易出错）。最佳值为0，表示变异与其在序列上的位置无关。负值表示变异位点更容易在末端出现，正值表示参考基因组中的等位基因更容易在末端出现。

SOR：也是一个用来评估是否存在链偏向性的参数，相当于FS的升级版。The StrandOddsRatio annotation is one of several methods that aims to evaluate whether there is strand bias in the data. It is an updated form of the Fisher Strand Test that is better at taking into account large amounts of data in high coverage situations. It is used to determine if there is strand bias between forward and reverse strands for the reference or alternate allele. The reported value is ln-scaled.

AC，AF & AN：AC（Allele Count）表示Allele的数目；AF（Allele Frequency）表示Allele的频率；AN（Allele Number） 表示Allele的总数目；对于1个diploid sample而言：则基因型 0/1 表示sample为杂合子，Allele数为1(双倍体的sample在该位点只有1个等位基因发生了突变)，Allele的频率为0.5(双倍体的sample在该位点只有50%的等位基因发生了突变)，总的Allele为2；基因型 1/1 则表示sample为纯合的，Allele数为2，Allele的频率为1，总的Allele为2。以我的结果为例:

9. FORMAT 为后面10列信息的说明列，通常以":"隔开各个缩写词。不同的变异检测软件可能会有差异，以下用GATK的检测结果为例：

10. 10列（包含）以后为样品基因型列，各信息以":"分隔与FORMAT列一一对应；

GT 表示genotype，通常用”/” or “|”分隔两个数字，“|”phase过也就是杂合的两个等位基因知道哪个等位基因来自哪条染色体；0代表参考基因组的碱基类型；1代表ALT碱基类型的第一个碱基（多个碱基用","分隔），2代表ALT第二个碱基，以此类推；比如 REF列为：A， ALT列为G,T；那么0/1基因型为AG 杂合，1/1基因型为GG纯合SNP；1/2代表GT基因型；./.表示缺失；

AD 两种碱基各自支持的碱基数量，用","分开两个数据，分别代表两个等位基因的深度；

DP 该样品该变异位点的测序深度总和，也就是AD两个数字的和；

PL 归一化后各基因型的可能性，通常有三个数字用’,'隔开，顺序对应AA,AB,BB基因型，A代表REF，B代表ALT(也就是0/0, 0/1, and 1/1)，由于是归一化之后，数值越小代表基因型越可靠；那么最小的数字对应的基因型判读为该样品的最可能的基因型；(provieds the likelihoods of the given genotypes)

GQ 针对PL的判读得到的基因型的质量值，此值越大基因型质量值越好。由于PL归一化之后通常最小的数字为0；那么基因型的质量值取PL中第二小的数字，如果第二小的数字大于99，我们只取99，因为在GATK中再大的值是没有意义的，第二小的数大于99的话一般说明基因型的判读是很可靠的，只有当第二小的数小于99的时候，才有必要怀疑基因型的可靠性；

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

使用方法：

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

执行以下命令：

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

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

执行以下命令：

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

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

筛选指定位置变异位点
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：没有此参数则不会输出

过滤指定缺失率的变异位点
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代表允许全部缺失。

计算snp缺失率
vcftools中有两个参数可以计算vcf文件中snp的缺失率。

分别是：

–missing-indv：生成一个文件，报告每个样品的缺失情况，该文件的后缀为“.imiss”。

–missing-site：生成一个文件，报告每个snp位点的缺失情况，该文件的后缀为“.lmiss”。

使用方法：

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列为缺失的总数，最后一列为缺失率。

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列为缺失的总数，最后一列为缺失率。

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

使用方法：

随机抽取5个样品，执行以下代码：

vcftools --vcf snp.vcf --max-indv 5 --remove-indels --recode --out outfilename

彩蛋：vcftools 超详细帮助文档：http://vcftools.sourceforge.net/man_latest.html