VCF格式的学习及对VCF文件的统计

部分摘自# VincentLuo91的博客

Part 1 VCF格式的学习

1.什么是vcf？
VCF是用于描述SNP，INDEL和SV结果的文本文件。在GATK软件中得到最好的支持，当然SAMtools得到的结果也是VCF格式，和GATK的VCF格式有点差别。
2. VCF的主体结构
VCF文件分为两部分内容：以“#”开头的注释部分；没有“#”开头的主体部分
3.VCF的10列的意义
1CHROM ： 参考序列名称
2POS:variant的位置；如果是INDEL的话，位置是INDEL的第一个碱基位置
3ID:variant的ID;比如在dbSNP中有该SNP的id，则会在此行给出；若没有，则用’.'表示其为一个novel variant
4REF:参考序列的碱基
5ALT:variant的碱基
6QUAL:Phred格式(Phred_scaled)的质量值，表 示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%,qual值与p成正比例
7FILTER:使用上一个QUAL值来进行过滤的话，是不够的。GATK能使用其它的方法来进行过滤，过滤结果中通过则该值为”PASS”;若variant不可靠，则该项不为”PASS”或”.”。
8INFO:这一行是variant的详细信息，内容很多，以下再具体详述。
9FORMAT:variants的格式，例如GT:AD:DP:GQ:PL
10_SAMPLES _:各个Sample的值，由BAM文件中的@RG下的SM标签所决定
4.vcf文件的基因型信息
GT:样品的基因型（genotype）。两个数字中间用’/'分 开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele； 1 表示样品中variant的allele； 2表示有第二个variant的allele。所以：
0/0表示sample中该位点为纯合位点，和REF的碱基类型一致
0/1表示sample中该位点为杂合突变，有REF和ALT两个基因型（部分碱基和REF碱基类型一致，部分碱基和ALT碱基类型一致）
1/1表示sample中该位点为纯合突变，总体突变类型和ALT碱基类型一致
1/2表示sample中该位点为杂合突变，有ALT1和ALT2两个基因型（部分和ALT1碱基类型一致，部分和ALT2碱基类型一致）
AD和DP:AD(Allele Depth)为sample中每一种allele的reads覆盖度,在diploid（二倍体，或可指代多倍型）中则是用逗号分隔的两个值，前者对应REF基因，后者对应ALT基因型
DP(Depth)为sample中该位点的覆盖度，是所支持的两个AD值（逗号前和逗号后）的加和
例如：
1/1:0,175:175—GT:AD(REF),AD(ALT):DP
0/1:79,96:175
1/2:0,20,56:76
这里的三种类型对应的DP值均是其对应的AD值的加和，1/1的175是0+175，0/1的175是79+96，1/2的76是0+20+56
GQ:基因型的质量值（Genotype Quality）。Phred格式（Phred_scaled）的质量值，表示在该位点该基因型存在的可能性；该值越高，则Genotype的可能性越大；计算方法：Phred值=-10log(1-P)，P为基因型存在的概率。（一般在final.snp.vcf文件中，该值为99，为99时，其可能性最大）
PL:指定的三种基因型的质量值（provieds the likelihoods of the given genotypes）；这三种指定的基因型为（0/0，0/1，1/1），这三种基因型的概率总和为1。该值越大，表明为该种基因型的可能性越小。Phred值=-10log(P)**，P为基因型存在的概率。最有可能的genotype的值为0
5. VCF第8列的信息
第8列的信息包括18种，都是以“TAG=Value”，并使用分号分隔的形式，其中很多的注释信息在VCF文件的头部注释中给出，下面对常用的TAG进行解释:
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
DP（reads覆盖度）
表示reads被过滤后的覆盖度
FS
FisherStrand的缩写，表示使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值，该值越小越好；如果该值较大，表示strand bias（正负链偏移）越严重，即所检测到的variants位点上，reads比对到正负义链上的比例不均衡。一般进行filter的时候，推荐保留FS<10~20的variants位点。GATK可设定FS参数。
ReadPosRandSum
Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias.当variants出现在reads尾部的时候，其结果可能不准确。该值用于衡量alternative allele（变异的等位基因）相比于reference allele（参考基因组等位基因），其variant位点是否匹配到reads更靠中部的位置。因此只有基因型是杂合且有一个allele和参考基因组一致的时候，才能计算该值。若该值为正值，表明和alternative allele相当于reference allele，落来reads更靠中部的位置；若该值是负值，则表示alternative allele相比于reference allele落在reads更靠尾部的位置。
进行filter的之后，推荐保留ReadPosRankSum>-1.65~-3.0的variant位点

MQRankSum
该值用于衡量alternative allele上reads的mapping quality与reference allele上reads的mapping quality的差异。若该值是负数值，则表明alternative allele比reference allele的reads mapping quality差。进行filter的时候，推荐保留MQRankSum>-1.65~-3.0的variant位点。

Part 2 对VCF文件的统计

1.统计一下vcf文件variant的质量值的分布
统计的VCF文件是在GATK HaplotypeCaller call出来的VCF 文件
KPGP-00001.HC.vcf

把KPGP-00001.HC.vcf中以#开头的头文件部分除掉，保留主文件,并把第6列variant质量值提出来存到KPGP-00001.HC.QUAL.txt中。
``
grep -v '^#' KPGP-00001.HC.vcf |cut -f 6 >KPGP-00001.HC.QUAL.txt

将KPGP-00001.HC.QUAL.txt按第一列数值排序后取最后10行
``
sort -k1,1n KPGP-00001.HC.QUAL.txt|tail

结果：

177246.77
184866.77
185614.77
186061.77
187101.77
190499.77
192745.77
193674.77
195675.77
196247.77

所以,质量值的最大为196247.77
以10-1000区间，步长为10的区间进行统计，剩下的按大于1000算
python脚本如下

import sys
args=sys.argv
filename=args[1]
numDict={}
OT="over_1000"
numDict[OT]=0
for line in open (filename):
 lineL=float(line.strip())
 for i in range(10,1000,10):
    if i-10 <= lineL <= i:
      if i not in numDict:
         numDict[i]=1
      else:
         numDict[i]+=1
 if lineL > 1000:
        numDict[OT]+=1

for k,v in numDict.items():
   print(k,v)

运行python脚本并将结果按第一列数值大小排序

python3 distribution.py KPGP-00001.HC.QUAL.txt >KPGP-00001.HC.QUAL.distribution.txt

sort -k1,1n KPGP-00001.HC.QUAL.distribution.txt >KPGP-00001.HC.QUAL.distribution.sort.txt

结果如下：

over_1000 848602
20 41126
30 35794
40 34196
50 33893
60 32039
70 32968
80 29885
90 31482
100 31294
110 32549
120 31780
130 33910
140 33254
150 34696
160 37889
170 37403
180 37569
190 41049
200 43558
210 42034
220 44231
230 46817
240 48619
250 50107
260 51146
270 52414
280 54788
290 56195
300 55299
310 57655
320 59335
330 58915
340 59563
350 60876
360 60676
370 60336
380 61357
390 61234
400 59689
410 59058
420 60001
430 59384
440 56861
450 55951
460 56236
470 54859
480 51336
490 50367
500 50711
510 48943
520 47230
530 45798
540 43991
550 42699
560 41013
570 39158
580 37561
590 36790
600 36122
610 34610
620 32404
630 31108
640 32030
650 30382
660 27705
670 27306
680 28260
690 27063
700 25393
710 24746
720 25701
730 24971
740 23089
750 22921
760 23543
770 22872
780 22520
790 22975
800 22249
810 22231
820 22805
830 22918
840 21630
850 21104
860 22502
870 23507
880 21907
890 20746
900 22812
910 23912
920 21318
930 21153
940 22891
950 22529
960 21540
970 22029
980 22957
990 21794

可以看出质量值低于20的很少，一般情况下，软件call到的variation的质量值低于20，我们会选择舍弃掉，因为这样的variation可信度不高

再分别对VQSR之后的KPGP-00001.HC.snps.VQSR.vcf和KPGP-00001.HC.indels.VQSR.vcf进行统计

首先是SNP：KPGP-00001.HC.snps.VQSR.vcf

#首先取出VQSR之后PASS的变异

grep "PASS" KPGP-00001.HC.snps.VQSR.vcf > KPGP-00001.HC.snps.VQSR.PASS.vcf
 
#提取出第6列质量值
cut -f 6 KPGP-00001.HC.snps.VQSR.PASS.vcf >KPGP-00001.HC.snps.VQSR.PASS.QUAL.txt

python3 distribution.py KPGP-00001.HC.snps.VQSR.PASS.QUAL.txt >KPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.txt
sort -k1,1n KPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.txt >KPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.sort.txt
cat KPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.sort.txt
20 7586
30 6978
40 7114
50 8739
60 8036
70 11621
80 8613
90 10986
100 10779
110 12786
120 12336
130 14106
140 14009
150 15026
160 19027
170 18099
180 18743
190 21315
200 24440
210 23806
220 25568
230 27696
240 29416
250 31908
260 32480
270 34080
280 35897
290 37987
300 37837
310 40023
320 41743
330 41722
340 42906
350 44285
360 44345
370 43904
380 45876
390 45802
400 44924
410 44462
420 45175
430 45410
440 43349
450 42603
460 43145
470 42375
480 39487
490 39193
500 39057
510 37388
520 36891
530 35623
540 34056
550 32612
560 31811
570 30527
580 29294
590 28211
600 27463
610 27180
620 25359
630 23971
640 24313
650 24081
660 21905
670 21021
680 21742
690 21399
700 20357
710 19156
720 20183
730 19629
740 18785
750 18204
760 19113
770 18403
780 17932
790 18952
800 18514
810 17973
820 18393
830 19646
840 18406
850 17691
860 18685
870 19866
880 19042
890 17716
900 19044
910 20381
920 19043
930 18190
940 19804
950 19731
960 18517
970 19511
980 20326
990 19140

可见，VQSR之后低质量值的变异就更少了。
INDEL同理，直接给结果

20 25080
30 21471
40 20634
50 21342
60 20962
70 22412
80 19481
90 21583
100 21830
110 23103
120 22728
130 24317
140 24168
150 26206
160 28923
170 28230
180 28881
190 32201
200 34403
210 33287
220 35318
230 37679
240 39749
250 41159
260 42035
270 43483
280 45919
290 47002
300 46716
310 49116
320 50613
330 50536
340 51441
350 52712
360 52555
370 52505
380 53665
390 53737
400 52396
410 51932
420 53208
430 52683
440 50411
450 49896
460 50297
470 49165
480 45873
490 45110
500 45542
510 44209
520 42649
530 41291
540 39920
550 38712
560 37208
570 35687
580 34286
590 33613
600 33056
610 31768
620 29761
630 28634
640 29663
650 28096
660 25537
670 25378
680 26364
690 25347
700 23682
710 23067
720 24142
730 23504
740 21770
750 21745
760 22306
770 21678
780 21417
790 21886
800 21158
810 21262
820 21900
830 22006
840 20846
850 20290
860 21667
870 22772
880 21170
890 20036
900 22105
910 23200
920 20650
930 20509
940 22271
950 21931
960 20946
970 21412
980 22444
990 21245

2.统计一下vcf文件杂合变异位点跟纯合变异位点的分布
首先是SNP

#取出第10行
cut -f 10 KPGP-00001.HC.snps.VQSR.PASS.vcf >KPGP-00001.HC.snps.VQSR.PASS.samples.txt

写一个python脚本

import sys
args=sys.argv
filename=args[1]
varDict={}
for line in open(filename):
  lineL=line.strip().split(":")
  variant=lineL[0]
  if variant not in varDict:
      varDict[variant]=1
  else:
      varDict[variant]+=1
for k,v in varDict.items():
   print(k,v,sep=":")

运行

python3 sample_count.py KPGP-00001.HC.snps.VQSR.PASS.samples.txt
1/1:1559238
1/2:1173
0/1:1709664

接下来是INDEL

cut -f 10 KPGP-00001.HC.indels.VQSR.PASS.vcf >KPGP-00001.HC.indels.VQSR.PASS.samples.txt

python3 sample_count.py KPGP-00001.HC.indels.VQSR.PASS.samples.txt

1/2:26765
1/1:1847390
0/1:2129613

可以看出，不论是SNP还是INDEL，杂合和纯合突变的比例大致为1：1，根据Jimmy老师的说法，本次call variation的步骤还算合理。

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,230评论 5赞 476
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,261评论 2赞 380
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,089评论 0赞 336
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,542评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,542评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,544评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,922评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,578评论 0赞 257
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,816评论 1赞 296
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,576评论 2赞 320
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,658评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,359评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,937评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,920评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,859评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,381评论 2赞 342

推荐阅读更多精彩内容