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

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

1.什么是vcf垃僚？
VCF是用于描述SNP，INDEL和SV結(jié)果的文本文件认烁。在GATK軟件中得到最好的支持，當(dāng)然SAMtools得到的結(jié)果也是VCF格式，和GATK的VCF格式有點(diǎn)差別蒸甜。
2. VCF的主體結(jié)構(gòu)
VCF文件分為兩部分內(nèi)容：以“#”開頭的注釋部分；沒有“#”開頭的主體部分
3.VCF的10列的意義
1CHROM ：參考序列名稱
2POS:variant的位置余佛；如果是INDEL的話柠新，位置是INDEL的第一個(gè)堿基位置
3ID:variant的ID;比如在dbSNP中有該SNP的id，則會(huì)在此行給出辉巡；若沒有恨憎，則用’.'表示其為一個(gè)novel variant
4REF:參考序列的堿基
5ALT:variant的堿基
6QUAL:Phred格式(Phred_scaled)的質(zhì)量值，表示在該位點(diǎn)存在variant的可能性郊楣；該值越高憔恳，則variant的可能性越大；計(jì)算方法：Phred值 = -10 * log (1-p) p為variant存在的概率; 通過計(jì)算公式可以看出值為10的表示錯(cuò)誤概率為0.1净蚤，該位點(diǎn)為variant的概率為90%,qual值與p成正比例
7FILTER:使用上一個(gè)QUAL值來進(jìn)行過濾的話钥组，是不夠的。GATK能使用其它的方法來進(jìn)行過濾今瀑，過濾結(jié)果中通過則該值為”PASS”;若variant不可靠程梦，則該項(xiàng)不為”PASS”或”.”。
8INFO:這一行是variant的詳細(xì)信息橘荠，內(nèi)容很多屿附，以下再具體詳述。
9FORMAT:variants的格式哥童，例如GT:AD:DP:GQ:PL
10_SAMPLES _:各個(gè)Sample的值拿撩，由BAM文件中的@RG下的SM標(biāo)簽所決定
4.vcf文件的基因型信息
GT:樣品的基因型（genotype）。兩個(gè)數(shù)字中間用’/'分開如蚜，這兩個(gè)數(shù)字表示雙倍體的sample的基因型压恒。0 表示樣品中有ref的allele影暴； 1 表示樣品中variant的allele； 2表示有第二個(gè)variant的allele探赫。所以：
0/0表示sample中該位點(diǎn)為純合位點(diǎn)型宙，和REF的堿基類型一致
0/1表示sample中該位點(diǎn)為雜合突變，有REF和ALT兩個(gè)基因型（部分堿基和REF堿基類型一致伦吠，部分堿基和ALT堿基類型一致）
1/1表示sample中該位點(diǎn)為純合突變妆兑，總體突變類型和ALT堿基類型一致
1/2表示sample中該位點(diǎn)為雜合突變，有ALT1和ALT2兩個(gè)基因型（部分和ALT1堿基類型一致毛仪，部分和ALT2堿基類型一致）
AD和DP:AD(Allele Depth)為sample中每一種allele的reads覆蓋度,在diploid（二倍體搁嗓，或可指代多倍型）中則是用逗號(hào)分隔的兩個(gè)值，前者對(duì)應(yīng)REF基因箱靴，后者對(duì)應(yīng)ALT基因型
DP(Depth)為sample中該位點(diǎn)的覆蓋度腺逛，是所支持的兩個(gè)AD值（逗號(hào)前和逗號(hào)后）的加和
例如：
1/1:0,175:175—GT:AD(REF),AD(ALT):DP
0/1:79,96:175
1/2:0,20,56:76
這里的三種類型對(duì)應(yīng)的DP值均是其對(duì)應(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ì)量值抛杨，表示在該位點(diǎn)該基因型存在的可能性够委；該值越高，則Genotype的可能性越大怖现；計(jì)算方法：Phred值=-10log(1-P)茁帽，P為基因型存在的概率。（一般在final.snp.vcf文件中屈嗤，該值為99脐雪，為99時(shí)，其可能性最大）
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”埋泵，并使用分號(hào)分隔的形式，其中很多的注釋信息在VCF文件的頭部注釋中給出，下面對(duì)常用的TAG進(jìn)行解釋:
AC丽声，AF和AN
AC(Allele Count) 表示該Allele的數(shù)目礁蔗；AF(Allele Frequency) 表示Allele的頻率； AN(Allele Number) 表示Allele的總數(shù)目雁社。對(duì)于1個(gè)diploid sample而言：則基因型 0/1 表示sample為雜合子浴井，Allele數(shù)為1(雙倍體的sample在該位點(diǎn)只有1個(gè)等位基因發(fā)生了突變)，Allele的頻率為0.5(雙倍體的 sample在該位點(diǎn)只有50%的等位基因發(fā)生了突變)霉撵，總的Allele為2磺浙；基因型 1/1 則表示sample為純合的，Allele數(shù)為2徒坡，Allele的頻率為1撕氧，總的Allele為2
DP（reads覆蓋度）
表示reads被過濾后的覆蓋度
FS
FisherStrand的縮寫，表示使用Fisher’s精確檢驗(yàn)來檢測(cè)strand bias而得到的Fhred格式的p值喇完，該值越小越好伦泥；如果該值較大，表示strand bias（正負(fù)鏈偏移）越嚴(yán)重何暮，即所檢測(cè)到的variants位點(diǎn)上，reads比對(duì)到正負(fù)義鏈上的比例不均衡铐殃。一般進(jìn)行filter的時(shí)候海洼，推薦保留FS<10~20的variants位點(diǎn)。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尾部的時(shí)候坏逢，其結(jié)果可能不準(zhǔn)確。該值用于衡量alternative allele（變異的等位基因）相比于reference allele（參考基因組等位基因）赘被，其variant位點(diǎn)是否匹配到reads更靠中部的位置是整。因此只有基因型是雜合且有一個(gè)allele和參考基因組一致的時(shí)候，才能計(jì)算該值民假。若該值為正值浮入，表明和alternative allele相當(dāng)于reference allele，落來reads更靠中部的位置羊异；若該值是負(fù)值事秀，則表示alternative allele相比于reference allele落在reads更靠尾部的位置。
進(jìn)行filter的之后野舶，推薦保留ReadPosRankSum>-1.65~-3.0的variant位點(diǎn)

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

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

1.統(tǒng)計(jì)一下vcf文件variant的質(zhì)量值的分布
統(tǒng)計(jì)的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é)果：

所以,質(zhì)量值的最大為196247.77
以10-1000區(qū)間袋哼，步長(zhǎng)為10的區(qū)間進(jìn)行統(tǒng)計(jì)，剩下的按大于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)

運(yùn)行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，我們會(huì)選擇舍棄掉楚堤，因?yàn)檫@樣的variation可信度不高

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

首先是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é)果

2.統(tǒng)計(jì)一下vcf文件雜合變異位點(diǎn)跟純合變異位點(diǎn)的分布
首先是SNP

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

寫一個(gè)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=":")

運(yùn)行

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的步驟還算合理。

最后編輯于：2018.08.04 11:47:50

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末嘿歌，一起剝皮案震驚了整個(gè)濱河市掸掏，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌宙帝，老刑警劉巖丧凤，帶你破解...
沈念sama閱讀 219,427評(píng)論 6贊 508
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場(chǎng)離奇詭異步脓，居然都是意外死亡愿待，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,551評(píng)論 3贊 395
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門靴患，熙熙樓的掌柜王于貴愁眉苦臉地迎上來仍侥，“玉大人，你說我怎么就攤上這事鸳君∨┰ǎ” “怎么了？”我有些...
開封第一講書人閱讀 165,747評(píng)論 0贊 356
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵或颊，是天一觀的道長(zhǎng)腿时。經(jīng)常有香客問我，道長(zhǎng)饭宾，這世上最難降的妖魔是什么批糟？我笑而不...
開封第一講書人閱讀 58,939評(píng)論 1贊 295
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮看铆，結(jié)果婚禮上徽鼎，老公的妹妹穿的比我還像新娘。我一直安慰自己，他們只是感情好否淤，可當(dāng)我...
茶點(diǎn)故事閱讀 67,955評(píng)論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布悄但。她就那樣靜靜地躺著，像睡著了一般石抡。火紅的嫁衣襯著肌膚如雪檐嚣。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 51,737評(píng)論 1贊 305
城市分裂傳說
那天啰扛，我揣著相機(jī)與錄音嚎京，去河邊找鬼。笑死隐解，一個(gè)胖子當(dāng)著我的面吹牛鞍帝，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播煞茫，決...
沈念sama閱讀 40,448評(píng)論 3贊 420
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼帕涌，長(zhǎng)吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼！你這毒婦竟也來了续徽？” 一聲冷哼從身側(cè)響起蚓曼，我...
開封第一講書人閱讀 39,352評(píng)論 0贊 276
萬榮殺人案實(shí)錄
序言：老撾萬榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎钦扭，沒想到半個(gè)月后纫版，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 45,834評(píng)論 1贊 317
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡土全，尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 37,992評(píng)論 3贊 338
?白月光啟示錄
正文我和宋清朗相戀三年捎琐，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了会涎。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片裹匙。...
茶點(diǎn)故事閱讀 40,133評(píng)論 1贊 351
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖末秃，靈堂內(nèi)的尸體忽然破棺而出概页，到底是詐尸還是另有隱情，我是刑警寧澤练慕，帶...
沈念sama閱讀 35,815評(píng)論 5贊 346
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布惰匙，位于F島的核電站，受9級(jí)特大地震影響铃将，放射性物質(zhì)發(fā)生泄漏项鬼。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 41,477評(píng)論 3贊 331
男人毒藥：我在死后第九天來索命
文/蒙蒙一劲阎、第九天我趴在偏房一處隱蔽的房頂上張望绘盟。院中可真熱鬧，春花似錦、人聲如沸龄毡。這莊子的主人今日做“春日...
開封第一講書人閱讀 32,022評(píng)論 0贊 22
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽(yáng)沦零。三九已至祭隔，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間路操，已是汗流浹背疾渴。一陣腳步聲響...
開封第一講書人閱讀 33,147評(píng)論 1贊 272
情欲美人皮
我被黑心中介騙來泰國(guó)打工，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留寻拂，地道東北人程奠。一個(gè)月前我還...
沈念sama閱讀 48,398評(píng)論 3贊 373
代替公主和親
正文我出身青樓，卻偏偏與公主長(zhǎng)得像祭钉，于是被迫代替她去往敵國(guó)和親瞄沙。傳聞我的和親對(duì)象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 45,077評(píng)論 2贊 355

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

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

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

推薦閱讀更多精彩內(nèi)容