生物基因資料檔案——vcf格式詳解
1.什麼是vcf檔案
VCF是用於描述SNP（單個鹼基上的變異），INDEL（插入缺失標記）和SV（結構變異位點）結果的文字檔案。在GATK軟體中得到最好的支援，當然SAMtools得到的結果也是VCF格式，和GATK的CVF格式有點差別。
2.VCF的主體結構
##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> 從範例上看，VCF檔案分為兩部分內容：以“#”開頭的註釋部分；沒有“#”開頭的主體部分。
值得注意的是，註釋部分有很多對VCF的介紹資訊。實際上不需要本文章，只是看看這個註釋部分就完全明白了VCF各行各列代表的意義。
主體部分中每一行代表一個Variant的資訊。
3.怎麼解釋Variation
CHROM：表示變異位點是在哪個contig 裡call出來的，如果是人類全基因組的話那就是chr1…chr22，chrX,Y,M。
POS： 變異位點相對於參考基因組所在的位置，如果是indel，就是第一個鹼基所在的位置。
ID： variant的ID。 如果call出來的SNP存在於dbSNP資料庫裡，就會顯示相應的dbSNP裡的rs編號；若沒有，則用’.’表示其為一個novel variant。
REF和ALT： 在這個變異位點處，參考基因組中所對應的鹼基和研究物件基因組（Variant）中所對應的鹼基。
QUAL： Phred格式(Phred_scaled)的質量值，可以理解為所call出來的變異位點的質量值。表 示在該位點存在variant的可能性；該值越高，則variant的可能性越大；
計算方法：① Q=-10*lgP，Q表示質量值；P表示這個位點發生錯誤的概率。
②Phred值Q = -10 * lg (1-p) ，p為variant存在的概率;
通過計算公式可以看出值為10的表示錯誤概率為0.1，該位點為variant的概率為90%。
同理，當Q=20時，錯誤率就控制在了0.01。
FILTER： 使用上一個QUAL值來進行過濾的話，是不夠的。理想情況下，QUAL這個值應該是用所有的錯誤模型算出來的，這個值就可以代表正確的變異位點了，但是事實是做不到的。因此，還需要對原始變異位點做進一步的過濾。無論你用什麼方法對變異位點進行過濾，過濾完了之後，在FILTER一欄都會留下過濾記錄，如果是通過了過濾標準，那麼這些通過標準的好的變異位點的FILTER一欄就會註釋一個PASS，如果沒有通過過濾，就會在FILTER這一欄提示除了PASS的其他資訊。如果這一欄是一個“.”的話，就說明沒有進行過任何過濾。
INFO： 這一行是variant的詳細資訊，內容很多，以下再具體詳述。
例子：
##fileformat=VCFv4.0
##FILTER= ##FORMAT= ##FORMAT= ##FORMAT= ##FORMAT= ##FORMAT= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##UnifiedGenotyperV2="analysis_type=UnifiedGenotyperV2 input_file=[TEXT CLIPPED FOR CLARITY]"
#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”，在後續分析中可以被過濾掉。
FORMAT 和 NA12878：這兩行合起來提供了’NA12878′這個sample的基因型的資訊。’NA12878′代表這該名稱的樣品，是由BAM檔案中的@RG下的 SM 標籤決定的。
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）。這三種指定的基因型(0/0,0/1,1/1)的概率總和為1。如果轉換成支援該基因型概率（P）的話，由於L=-10lgP，那麼P=10^（-L/10），因此，當L值為0時，P=10^0=1。因此，這個值越小，支援概率就越大，也就是說是這個基因型的可能性越大。
GQ： 表示最可能的基因型的質量值。表示的意義同QUAL。Phred格式(Phred_scaled)的質量值，表示在該位點該基因型存在的可能性；該值越高，則Genotype的可能性越 大；計算方法：Phred值 = -10 * log (1-p) p為基因型存在的概率。
舉個例子說明一下：
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第8列的資訊
該列資訊最多了，都是以 “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被過濾掉後的覆蓋度。
Dels： Fraction of Reads Containing Spanning Deletions。進行SNP和INDEL calling的結果中，有該TAG並且值為0表示該位點為SNP，沒有則為INDEL。
FS：使用Fisher’s精確檢驗來檢測strand bias而得到的Fhred格式的p值。該值越小越好。一般進行filter的時候，可以設定 FS < 10～20。
HaplotypeScore： Consistency of the site with at most two segregating haplotypes
InbreedingCoeff： Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hard-Weinberg expectation
MLEAC： Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed
MLEAF： Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT alle in the same order as listed
MQ： RMS Mapping Quality
MQ0： Total Mapping Quality Zero Reads
MQRankSum： Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities
QD： Variant Confidence/Quality by Depth
RPA： Number of times tandem repeat unit is repeated, for each allele (including reference)
RU： Tandem repeat unit (bases)
ReadPosRankSum： Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias
STR： Variant is a short tandem repeat
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
This is under continued development, please check the hts-specs page for the most recent specification
REF: