VCF文件示例(VCFv4.2)
##fileformat=VCFv4.2 ##fileDate=20090805 ##source=myImputationProgramV3.1 ##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x> ##phasing=partial ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data"> ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth"> ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency"> ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele"> ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129"> ##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership"> ##FILTER=<ID=q10,Description="Quality below 10"> ##FILTER=<ID=s50,Description="Less than 50% of samples have data"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality"> ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth"> ##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003 20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,. 20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3 20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4 20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2 20 1234567 microsat1 GTC G,GTCT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3
CHROM: 表示變異位點是在哪個contig 里call出來的,如果是人類全基因組的話那就是chr1…chr22,chrX,Y,M。
POS: 變異位點相對於參考基因組所在的位置,如果是indel,就是第一個鹼基所在的位置。
ID: 如果call出來的SNP存在於dbSNP數據庫里,就會顯示相應的dbSNP里的rs編號。
REF和REF: 在這個變異位點處,參考基因組中所對應的鹼基和研究對象基因組中所對應的鹼基。
QUAL: 可以理解為所call出來的變異位點的質量值。Q=-10lgP,Q表示質量值;P表示這個位點發生錯誤的概率。因此,如果想把錯誤率從控制在90%以上,P的閾值就是1/10,那lg(1/10)=-1,Q=(-10)*(-1)=10。同理,當Q=20時,錯誤率就控制在了0.01。
FILTER: 理想情況下,QUAL這個值應該是用所有的錯誤模型算出來的,這個值就可以代表正確的變異位點了,但是事實是做不到的。因此,還需要對原始變異位點做進一步的過濾。無論你用什么方法對變異位點進行過濾,過濾完了之后,在FILTER一欄都會留下過濾記錄,如果是通過了過濾標准,那么這些通過標准的好的變異位點的FILTER一欄就會注釋一個PASS,如果沒有通過過濾,就會在FILTER這一欄提示除了PASS的其他信息。如果這一欄是一個“.”的話,就說明沒有進行過任何過濾。
例子:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878 chr1 873762 . T G 5231.78 PASS AC=1;AF=0.50;AN=2;DP=315;Dels=0.00;HRun=2;HaplotypeScore=15.11;MQ=91.05;MQ0=15;QD=16.61;SB=-1533.02;VQSLOD=-1.5473 GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255 chr1 877664 rs3828047 A G 3931.66 PASS AC=2;AF=1.00;AN=2;DB;DP=105;Dels=0.00;HRun=1;HaplotypeScore=1.59;MQ=92.52;MQ0=4;QD=37.44;SB=-1152.13;VQSLOD= 0.1185 GT:AD:DP:GQ:PL 1/1:0,105:94:99:255,255,0 chr1 899282 rs28548431 C T 71.77 PASS AC=1;AF=0.50;AN=2;DB;DP=4;Dels=0.00;HRun=0;HaplotypeScore=0.00;MQ=99.00;MQ0=0;QD=17.94;SB=-46.55;VQSLOD=-1.9148 GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26 chr1 974165 rs9442391 T C 29.84 LowQual AC=1;AF=0.50;AN=2;DB;DP=18;Dels=0.00;HRun=1;HaplotypeScore=0.16;MQ=95.26;MQ0=0;QD=1.66;SB=-0.98 GT:AD:DP:GQ:PL 0/1:14,4:14:60.91:61,0,255
到現在,我們就可以解釋上面的例子:
chr1:873762 是一個新發現的T/G變異,並且有很高的可信度(qual=5231.78)。
chr1:877664 是一個已知的變異為A/G 的SNP位點,名字rs3828047,並且具有很高的可信度(qual=3931.66)。
chr1:899282 是一個已知的變異為C/T的SNP位點,名字rs28548431,但可信度較低(qual=71.77)。
chr1:974165 是一個已知的變異為T/C的SNP位點,名字rs9442391,但是這個位點的質量值很低,被標 成了“LowQual”,在后續分析中可以被過濾掉。
Vcf文件看起來很復雜,挺嚇人的樣子,但是里面大部分都是一些tags,而這些tags基本上都是在VASR中過濾用的,能夠理解每個tags的意思最好,如果實在不理解也就不用管了。其實最關鍵的信息也就是那么幾列:
chr1 873762 . T G [CLIPPED] GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255
chr1 877664 rs3828047 A G [CLIPPED] GT:AD:DP:GQ:PL 1/1:0,105:94:99:255,255,0
chr1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26
其中最后面兩列是相對應的,每一個tag對應一個或者一組值,如:
chr1:873762,GT對應0/1;AD對應173,141;DP對應282;GQ對應99;PL對應255,0,255。
GT: 表示這個樣本的基因型,對於一個二倍體生物,GT值表示的是這個樣本在這個位點所攜帶的兩個等位基因。0表示跟REF一樣;1表示表示跟ALT一樣;2表示第二個ALT。當只有一個ALT 等位基因的時候,0/0表示純和且跟REF一致;0/1表示雜合,兩個allele一個是ALT一個是REF;1/1表示純和且都為ALT; The most common format subfield is GT (genotype) data. If the GT subfield is present, it must be the first subfield. In the sample data, genotype alleles are numeric: the REF allele is 0, the first ALT allele is 1, and so on. The allele separator is '/' for unphased genotypes and '|' for phased genotypes.
0 - reference call
1 - alternative call 1
2 - alternative call 2
AD: 對應兩個以逗號隔開的值,這兩個值分別表示覆蓋到REF和ALT鹼基的reads數,相當於支持REF和支持ALT的測序深度。
DP: 覆蓋到這個位點的總的reads數量,相當於這個位點的深度(並不是多有的reads數量,而是大概一定質量值要求的reads數)。
PL: 對應3個以逗號隔開的值,這三個值分別表示該位點基因型是0/0,0/1,1/1的沒經過先驗的標准化Phred-scaled似然值(L)。如果轉換成支持該基因型概率(P)的話,由於L=-10lgP,那么P=10^(-L/10),因此,當L值為0時,P=10^0=1。因此,這個值越小,支持概率就越大,也就是說是這個基因型的可能性越大。
GQ: 表示最可能的基因型的質量值。表示的意義同QUAL。
舉個例子說明一下:
chr1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26
在這個位點,GT=0/1,也就是說這個位點的基因型是C/T;GQ=25.92,質量值並不算太高,可能是因為cover到這個位點的reads數太少,DP=4,也就是說只有4條reads支持這個地方的變異;AD=1,3,也就是說支持REF的read有一條,支持ALT的有3條;在PL里,這個位點基因型的不確定性就表現的更突出了,0/1的PL值為0,雖然支持0/1的概率很高;但是1/1的PL值只有26,也就是說還有10^(-2.6)=0.25%的可能性是1/1;但幾乎不可能是0/0,因為支持0/0的概率只有10^(-10.3)=5*10-11。
VCF (Variant Call Format) version 4.1
The VCF specification is no longer maintained by the 1000 Genomes Project. The group leading the management and expansion of the format is the Global Alliance for Genomics and Health Data Working group file format team, http://ga4gh.org/#/fileformats-team
The main version of the specification can be found on https://github.com/samtools/hts-specs
This is under continued development, please check the hts-specs page for the most recent specification
A PDF of the v4.1 spec is http://samtools.github.io/hts-specs/VCFv4.1.pdf
A PDF of the v4.2 spec is http://samtools.github.io/hts-specs/VCFv4.2.pdf
VCFTools host a discussion list about the specification called vcf-spec http://sourceforge.net/p/vcftools/mailman/
REF:
http://blog.sina.com.cn/s/blog_12d5e3d3c0101qv1u.html