vcf格式文件基础知识与编辑操作详解

vcf格式文件基础知识与编辑操作详解

1. VCF文件基础知识介绍：

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

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

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

下面是一个典型VCF文件的部分示例（可左右拖动）：


##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS     ID        REF ALT    QUAL FILTER INFO                              FORMAT      NA00001        NA00002        NA00003
20     14370   rs6054257 G      A       29   PASS   NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
20     17330   .         T      A       3    q10    NS=3;DP=11;AF=0.017               GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3   0/0:41:3
20     1110696 rs6040355 A      G,T     67   PASS   NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2   2/2:35:4
20     1230237 .         T      .       47   PASS   NS=3;DP=13;AA=T                   GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
20     1234567 microsat1 GTCT   G,GTACT 50   PASS   NS=3;DP=9;AA=G                    GT:GQ:DP    0/1:35:4       0/2:17:2       1/1:40:3

整体说明信息（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列开始为样品信息列，可以无限多个；

  • #CHROM

  • POS

  • ID

  • REF

  • ALT

  • QUAL

  • FILTER

  • INFO

  • FORMAT

  • 后面的列都为样品基因型信息列

具体说明如下

  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覆盖为缺失。

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

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

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

  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)，由于是归一化之后，数值越小代表基因型越可靠；那么最小的数字对应的基因型判读为该样品的最可能的基因型；

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

介绍到这里，大家应该对VCF文件有了一个初步理解，如有其他问题请下方留言。

2. VCF文件处理工具介绍：

筛选特定条件下的SNP/Indel信息是比较常见的个性化分析需求，不过VCF文件通常比较大，不建议在Windows当中处理。那么我们怎么筛选处理VCF文件呢？这里介绍一个Linux中非常好用的vcf处理工具，vcftools http://vcftools.sourceforge.net/man_latest.html。这是他的帮助链接; 功能非常丰富；vcftools是一款处理vcf文件的工具，体积小、运行速度快，可以按位点信息、碱基位置等信息筛选特定条件下的变异信息，也可以比较两个vcf文件之间变异信息的差异，还可以对VCF文件拆分，格式转换，质量过滤筛选等等功能；强烈推荐使用；这里介绍一下我在工作中用过的一些命令选项。

过滤变异类型

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使用手册：vcftools使用手册

vcftools就讲到这了，有兴趣的同学赶快试一下吧！


此外，我们在网易云课堂上有各种教学视频，有兴趣可以了解一下：

1. 文章越来越难发？是你没发现新思路，基因家族分析发2-4分文章简单快速，学习链接：基因家族分析实操课程

2. 转录组数据理解不深入？图表看不懂？点击链接学习深入解读数据结果文件，学习链接：转录组（有参）结果解读；转录组（无参）结果解读

3. 转录组数据深入挖掘技能-WGCNA，提升你的文章档次，学习链接：WGCNA-加权基因共表达网络分析

4. 转录组数据怎么挖掘？学习链接：转录组标准分析后的数据挖掘

5. 微生物16S/ITS/18S分析原理及结果解读

6. 更多学习内容：linux、perl、R语言画图，更多免费课程请点击以下链接：

https://study.omicsclass.com/

  • 发表于 2019-05-24 15:13
  • 阅读 ( 16327 )
  • 分类：软件工具

你可能感兴趣的文章

相关问题

0 条评论

请先 登录 后评论
安生水
安生水

348 篇文章

作家榜 »

  1. omicsgene 700 文章
  2. 安生水 348 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. xun 82 文章
  6. 红橙子 78 文章
  7. rzx 75 文章
  8. CORNERSTONE 72 文章

安生水

如果觉得我的文章对您有用，请随意打赏。你的支持将鼓励我继续创作！