VCF格式的學(xué)習(xí)及對VCF文件的統(tǒng)計

部分摘自# VincentLuo91的博客

Part 1 VCF格式的學(xué)習(xí)

1.什么是vcf？
VCF是用于描述SNP，INDEL和SV結(jié)果的文本文件。在GATK軟件中得到最好的支持，當(dāng)然SAMtools得到的結(jié)果也是VCF格式，和GATK的VCF格式有點差別。
2. VCF的主體結(jié)構(gòu)
VCF文件分為兩部分內(nèi)容：以“#”開頭的注釋部分；沒有“#”開頭的主體部分
3.VCF的10列的意義
1CHROM ： 參考序列名稱
2POS:variant的位置；如果是INDEL的話，位置是INDEL的第一個堿基位置
3ID:variant的ID;比如在dbSNP中有該SNP的id，則會在此行給出；若沒有，則用’.'表示其為一個novel variant
4REF:參考序列的堿基
5ALT:variant的堿基
6QUAL:Phred格式(Phred_scaled)的質(zhì)量值，表 示在該位點存在variant的可能性；該值越高，則variant的可能性越大；計算方法：Phred值 = -10 * log (1-p) p為variant存在的概率; 通過計算公式可以看出值為10的表示錯誤概率為0.1，該位點為variant的概率為90%,qual值與p成正比例
7FILTER:使用上一個QUAL值來進行過濾的話，是不夠的。GATK能使用其它的方法來進行過濾，過濾結(jié)果中通過則該值為”PASS”;若variant不可靠，則該項不為”PASS”或”.”。
8INFO:這一行是variant的詳細(xì)信息，內(nèi)容很多，以下再具體詳述。
9FORMAT:variants的格式，例如GT:AD:DP:GQ:PL
10_SAMPLES _:各個Sample的值，由BAM文件中的@RG下的SM標(biāo)簽所決定
4.vcf文件的基因型信息
GT:樣品的基因型（genotype）。兩個數(shù)字中間用’/'分 開，這兩個數(shù)字表示雙倍體的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（二倍體，或可指代多倍型）中則是用逗號分隔的兩個值，前者對應(yīng)REF基因，后者對應(yīng)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
這里的三種類型對應(yīng)的DP值均是其對應(yīng)的AD值的加和，1/1的175是0+175，0/1的175是79+96，1/2的76是0+20+56
GQ:基因型的質(zhì)量值（Genotype Quality）。Phred格式（Phred_scaled）的質(zhì)量值，表示在該位點該基因型存在的可能性；該值越高，則Genotype的可能性越大；計算方法：Phred值=-10log(1-P)，P為基因型存在的概率。（一般在final.snp.vcf文件中，該值為99，為99時，其可能性最大）
PL:指定的三種基因型的質(zhì)量值（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的數(shù)目；AF(Allele Frequency) 表示Allele的頻率； AN(Allele Number) 表示Allele的總數(shù)目。對于1個diploid sample而言：則基因型 0/1 表示sample為雜合子，Allele數(shù)為1(雙倍體的sample在該位點只有1個等位基因發(fā)生了突變)，Allele的頻率為0.5(雙倍體的 sample在該位點只有50%的等位基因發(fā)生了突變)，總的Allele為2； 基因型 1/1 則表示sample為純合的，Allele數(shù)為2，Allele的頻率為1，總的Allele為2
DP（reads覆蓋度）
表示reads被過濾后的覆蓋度
FS
FisherStrand的縮寫，表示使用Fisher’s精確檢驗來檢測strand bias而得到的Fhred格式的p值，該值越小越好；如果該值較大，表示strand bias（正負(fù)鏈偏移）越嚴(yán)重，即所檢測到的variants位點上，reads比對到正負(fù)義鏈上的比例不均衡。一般進行filter的時候，推薦保留FS<10~20的variants位點。GATK可設(shè)定FS參數(shù)。
ReadPosRandSum
Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias.當(dāng)variants出現(xiàn)在reads尾部的時候，其結(jié)果可能不準(zhǔn)確。該值用于衡量alternative allele（變異的等位基因）相比于reference allele（參考基因組等位基因），其variant位點是否匹配到reads更靠中部的位置。因此只有基因型是雜合且有一個allele和參考基因組一致的時候，才能計算該值。若該值為正值，表明和alternative allele相當(dāng)于reference allele，落來reads更靠中部的位置；若該值是負(fù)值，則表示alternative allele相比于reference allele落在reads更靠尾部的位置。
進行filter的之后，推薦保留ReadPosRankSum>-1.65~-3.0的variant位點

MQRankSum
該值用于衡量alternative allele上reads的mapping quality與reference allele上reads的mapping quality的差異。若該值是負(fù)數(shù)值，則表明alternative allele比reference allele的reads mapping quality差。進行filter的時候，推薦保留MQRankSum>-1.65~-3.0的variant位點。

Part 2 對VCF文件的統(tǒng)計

1.統(tǒng)計一下vcf文件variant的質(zhì)量值的分布
統(tǒng)計的VCF文件是在GATK HaplotypeCaller call出來的VCF 文件
KPGP-00001.HC.vcf

把KPGP-00001.HC.vcf中以#開頭的頭文件部分除掉，保留主文件,并把第6列variant質(zhì)量值提出來存到KPGP-00001.HC.QUAL.txt中。
``
grep -v '^#' KPGP-00001.HC.vcf |cut -f 6 >KPGP-00001.HC.QUAL.txt

將KPGP-00001.HC.QUAL.txt按第一列數(shù)值排序后取最后10行
``
sort -k1,1n KPGP-00001.HC.QUAL.txt|tail

結(jié)果：

177246.77
184866.77
185614.77
186061.77
187101.77
190499.77
192745.77
193674.77
195675.77
196247.77

所以,質(zhì)量值的最大為196247.77
以10-1000區(qū)間，步長為10的區(qū)間進行統(tǒng)計，剩下的按大于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腳本并將結(jié)果按第一列數(shù)值大小排序

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

結(jié)果如下：

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

可以看出質(zhì)量值低于20的很少，一般情況下，軟件call到的variation的質(zhì)量值低于20，我們會選擇舍棄掉，因為這樣的variation可信度不高

再分別對VQSR之后的KPGP-00001.HC.snps.VQSR.vcf和KPGP-00001.HC.indels.VQSR.vcf進行統(tǒng)計

首先是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列質(zhì)量值
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之后低質(zhì)量值的變異就更少了。
INDEL同理，直接給結(jié)果

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.統(tǒng)計一下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，根據(jù)Jimmy老師的說法，本次call variation的步驟還算合理。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 劉小澤寫于18.7.17所有的數(shù)據(jù)，一旦要找變異位點信息，就離不開VCF。豆豆也是在寫一個重測序的操作流程，遇到了...
    劉小澤閱讀 34,120評論 2贊 73
  • 轉(zhuǎn)自：https://blog.csdn.net/sinat_38163598/article/details/7...
    簡單點lili閱讀 4,484評論 0贊 9
  • 錢鐘書曾經(jīng)說過，要想結(jié)為夫妻的話，那么就先去旅行一次。旅行，是最考驗情感的時候，有的人因為短暫的旅行，收獲愛情，有...
    是芋頭仔啊閱讀 850評論 0贊 3
  • 睡懶覺的人（wm197049） 不知道究竟是種進步還是落后，睡懶覺始終具有無可抵擋的魅力。慵懶不是叫人厭惡憎恨的惡...
    阿尚青子自由寫作人閱讀 796評論 2贊 1
  • 美式風(fēng)格-大氣隨意 美式風(fēng)格，顧名思義是來自于美國的裝修和裝飾風(fēng)格。是殖民地風(fēng)格中最著名的代表風(fēng)格，某種意義上已經(jīng)...
    關(guān)山北X閱讀 207評論 0贊 0

友情鏈接更多精彩內(nèi)容