部分摘自# 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的步骤还算合理。
网友评论