文章來(lái)源 http://www.reibang.com/p/38f734ae47f5
?
?0.9?2018.08.03 14:52*?字?jǐn)?shù) 2060?閱讀 2137評(píng)論 0喜歡 15
部分摘自#?VincentLuo91的博客
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)容：以“#”開(kāi)頭的注釋部分；沒(méi)有“#”開(kāi)頭的主體部分
3.VCF的10列的意義
1CHROM?： 參考序列名稱(chēng)
2POS:variant的位置懦铺；如果是INDEL的話，位置是INDEL的第一個(gè)堿基位置
3ID:variant的ID;比如在dbSNP中有該SNP的id支鸡，則會(huì)在此行給出冬念；若沒(méi)有趁窃，則用’.'表示其為一個(gè)novel variant
4REF:參考序列的堿基
5ALT:variant的堿基
6QUAL:Phred格式(Phred_scaled)的質(zhì)量值，表 示在該位點(diǎn)存在variant的可能性急前；該值越高醒陆，則variant的可能性越大；計(jì)算方法：Phred值 = -10 * log (1-p) p為variant存在的概率; 通過(guò)計(jì)算公式可以看出值為10的表示錯(cuò)誤概率為0.1裆针，該位點(diǎn)為variant的概率為90%,qual值與p成正比例
7FILTER:使用上一個(gè)QUAL值來(lái)進(jìn)行過(guò)濾的話刨摩，是不夠的。GATK能使用其它的方法來(lái)進(jìn)行過(guò)濾据块，過(guò)濾結(jié)果中通過(guò)則該值為”P(pán)ASS”;若variant不可靠码邻，則該項(xiàng)不為”P(pán)ASS”或”.”。
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ù)字中間用’/'分 開(kāi)戈轿，這兩個(gè)數(shù)字表示雙倍體的sample的基因型。0 表示樣品中有ref的allele阵子； 1 表示樣品中variant的allele思杯； 2表示有第二個(gè)variant的allele。所以：
0/0表示sample中該位點(diǎn)為純合位點(diǎn)挠进，和REF的堿基類(lèi)型一致
0/1表示sample中該位點(diǎn)為雜合突變色乾，有REF和ALT兩個(gè)基因型（部分堿基和REF堿基類(lèi)型一致，部分堿基和ALT堿基類(lèi)型一致）
1/1表示sample中該位點(diǎn)為純合突變领突，總體突變類(lèi)型和ALT堿基類(lèi)型一致
1/2表示sample中該位點(diǎn)為雜合突變暖璧，有ALT1和ALT2兩個(gè)基因型（部分和ALT1堿基類(lèi)型一致，部分和ALT2堿基類(lèi)型一致）
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
這里的三種類(lèi)型對(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被過(guò)濾后的覆蓋度
FS
FisherStrand的縮寫(xiě)屎债，表示使用Fisher’s精確檢驗(yàn)來(lái)檢測(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，落來(lái)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出來(lái)的VCF 文件
KPGP-00001.HC.vcf
把KPGP-00001.HC.vcf中以#開(kāi)頭的頭文件部分除掉秋柄，保留主文件,并把第6列variant質(zhì)量值提出來(lái)存到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.77184866.77185614.77186061.77187101.77190499.77192745.77193674.77195675.77196247.77
所以,質(zhì)量值的最大為196247.77
以10-1000區(qū)間蠢正，步長(zhǎng)為10的區(qū)間進(jìn)行統(tǒng)計(jì)骇笔，剩下的按大于1000算
python腳本如下
importsysargs=sys.argvfilename=args[1]numDict={}OT="over_1000"numDict[OT]=0forlineinopen (filename): lineL=float(line.strip())foriinrange(10,1000,10):ifi-10<= lineL <= i:ifinotinnumDict:? ? ? ? numDict[i]=1else:? ? ? ? numDict[i]+=1iflineL >1000:? ? ? ? numDict[OT]+=1fork,vinnumDict.items():? print(k,v)
運(yùn)行python腳本并將結(jié)果按第一列數(shù)值大小排序
python3distribution.pyKPGP-00001.HC.QUAL.txt>KPGP-00001.HC.QUAL.distribution.txtsort-k1,1nKPGP-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-f6KPGP-00001.HC.snps.VQSR.PASS.vcf>KPGP-00001.HC.snps.VQSR.PASS.QUAL.txtpython3distribution.pyKPGP-00001.HC.snps.VQSR.PASS.QUAL.txt>KPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.txtsort-k1,1nKPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.txt>KPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.sort.txtcatKPGP-00001.HC.snps.VQSR.PASS.QUAL.distr.sort.txt20 758630 697840 711450 873960 803670 1162180 861390 10986100 10779110 12786120 12336130 14106140 14009150 15026160 19027170 18099180 18743190 21315200 24440210 23806220 25568230 27696240 29416250 31908260 32480270 34080280 35897290 37987300 37837310 40023320 41743330 41722340 42906350 44285360 44345370 43904380 45876390 45802400 44924410 44462420 45175430 45410440 43349450 42603460 43145470 42375480 39487490 39193500 39057510 37388520 36891530 35623540 34056550 32612560 31811570 30527580 29294590 28211600 27463610 27180620 25359630 23971640 24313650 24081660 21905670 21021680 21742690 21399700 20357710 19156720 20183730 19629740 18785750 18204760 19113770 18403780 17932790 18952800 18514810 17973820 18393830 19646840 18406850 17691860 18685870 19866880 19042890 17716900 19044910 20381920 19043930 18190940 19804950 19731960 18517970 19511980 20326990 19140
可見(jiàn)雹舀，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)計(jì)一下vcf文件雜合變異位點(diǎn)跟純合變異位點(diǎn)的分布
首先是SNP
#取出第10行cut-f10KPGP-00001.HC.snps.VQSR.PASS.vcf>KPGP-00001.HC.snps.VQSR.PASS.samples.txt
寫(xiě)一個(gè)python腳本
importsysargs=sys.argvfilename=args[1]varDict={}forlineinopen(filename):? lineL=line.strip().split(":")? variant=lineL[0]ifvariantnotinvarDict:? ? ? varDict[variant]=1else:? ? ? varDict[variant]+=1fork,vinvarDict.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
接下來(lái)是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老師的說(shuō)法签财，本次call variation的步驟還算合理串慰。
