vcftools是用來處理vcf和bcf文件的工具集，功能涵蓋了過濾，數(shù)據(jù)格式轉(zhuǎn)換，一些指標(biāo)的統(tǒng)計計算，vcf文件之間變異信息的差異比較等等。
1.基本用法
vcftools [ --vcf FILE | --gzvcf FILE | --bcf FILE] \
[ --out OUTPUT PREFIX ] [ FILTERING OPTIONS ] [ OUTPUT OPTIONS ]
2.使用
官網(wǎng)提供了幾個常用示例，下面用我的vcf文件測試一下，我的vcf文件共249個樣本，27012個SNP標(biāo)記
2.1 對來自chr1的每一個位點統(tǒng)計其基因頻率
vcftools --gzvcf combined200.vcf.gz --freq --chr chr1 --out chr1_analysis
因為沒有加任何過濾選項，所以沒有過濾任何記錄
$ less chr1_analysis.log | tail -n 4
After filtering, kept 249 out of 249 Individuals
Outputting Frequency Statistics...
After filtering, kept 3453 out of a possible 27012 Sites #chr1上原本就有3453個SNP位點
Run Time = 44.00 seconds
第三列表示在所有樣本中，該位點出現(xiàn)了幾種堿基，大多數(shù)情況下為2，也可能是3、4；第四列表示該位點出現(xiàn)的堿基總數(shù)，一個樣本貢獻(xiàn)2個。
$ head -n 5 chr1_analysis.frq
CHROM POS N_ALLELES N_CHR {ALLELE:FREQ}
chr1 1146710 2 494 G:0.983806 A:0.0161943
chr1 1146714 2 494 G:0.989879 A:0.0101215
chr1 1146763 2 494 G:0.995951 A:0.00404858
chr1 1146783 2 494 C:0.995951 T:0.00404858
2.2 過濾掉vcf文件中的InDel記錄，只保留SNP記錄
vcftools --gzvcf combined200.vcf.gz --remove-indels --recode --recode-INFO-all --out SNPs_only
--recode
表示過濾之后會生成一個新文件，以.recode.vcf為后綴
--recode-INFO-all
因為在過濾之后，原先存在的INFO列的注釋信息可能不對，
比如剔除了一些樣本，那么AN就需要重新計算。
所以過濾之后默認(rèn)情況下是不含有INFO列的。
該選項表示將原來vcf文件中所有INFO信息保留下來。
--recode-INFO <string>
比如--recode-INFO AC表示只將AC保留
--keep-only-indels和--remove-indels作用相反
根據(jù)vcf文件的FILTER列來過濾
--remove-filtered-all 過濾掉FILTER列不是“PASS”的記錄，這一步grep也能做。
vcftools --vcf combined3000_head3k.vcf --remove-filtered-all --recode --out SNPs_filter
--max-missing
--max-missing的取值是0-1，為1時表示某個位點上所有的樣本必須都有基因型，一個樣本的基因型都不能缺。所以這個選項可以理解為：能分型的樣本占總樣本的比例至少為多少。
vcftools --gzvcf combined200.vcf.gz --max-missing 0.99 --recode --out output_noMissing
有過濾操作都要加上--recode
2.3 比較兩個vcf文件的變異位點
vcftools --gzvcf combined200.vcf.gz --diff output_noMissing.recode.vcf --diff-site --out in1_v_in2
下面的表頭的1、2分別表示文件1、文件2，第4列1或2表示僅文件1或文件2有，B表示都有
$ head in1_v_in2.diff.sites_in_files
CHROM POS1 POS2 IN_FILE REF1 REF2 ALT1 ALT2
chr1 1146710 1146710 B G G A A
chr1 1146714 1146714 B G G A A
chr1 1146763 1146763 B G G A A
chr1 1146783 1146783 B C C T T
chr1 1146785 1146785 B G G C C
chr1 1146852 1146852 B C C T T
chr1 1146853 1146853 B G G A A
chr1 1147062 . 1 G . A .
chr1 1147243 . 1 A . C .
2.4 結(jié)合管道操作符
用--stdout代替--out即可結(jié)合管道進(jìn)行其他處理
2.5 計算Hardy-Weinberg p-value
Hardy-Weinberg遺傳平衡檢驗用到的統(tǒng)計方法是卡方測驗，詳細(xì)過程參考：http://www.dxy.cn/bbs/thread/223387#223387以及https://wenku.baidu.com/view/9820fc2258fb770bf78a5571.html
vcftools --gzvcf combined200.vcf.gz --hardy --out each_site_hardy
# Only using biallelic SNPs
第6列即為所求
$ lsx each_site_hardy.hwe | head -n 2
CHR POS OBS(HOM1/HET/HOM2) E(HOM1/HET/HOM2) ChiSq_HWE P_HWE P_HET_DEFICIT P_HET_EXCESS
chr1 1146710 239/8/0 239.06/7.87/0.06 6.692747e-02 1.000000e+00 1.000000e+00 9.440689e-01
3.其它參數(shù)說明
3.1 基本選項
--vcf | --gzvcf | --bcf 輸入格式，三選一即可
--out <輸出文件前綴>
--stdout 和 -c 輸出到標(biāo)準(zhǔn)輸出，能通過管道與其它命令結(jié)合
3.2 過濾SNP的選項
POSITION FILTERING
可以提供一些位置信息，或是能表示位置的文件，具體參數(shù)見官方文檔
ALLELE FILTERING
--maf 和 --max-maf用來限定最小等位基因頻率(MAF)的范圍
--mac 和 --max-mac和上面類似，用來限定最小等位基因的數(shù)量
--min-alleles 和 --max-alleles用來限定等位基因的數(shù)量，比如這條命令只會保留vcf文件ALT列只有一種堿基的情況。
vcftools --gzvcf combined200.vcf.gz --min-alleles 2 --max-alleles 2
GENOTYPE VALUE FILTERING
--min-meanDP 和 --max-meanDP用來限定所有樣本DP的平均值，DP表示某一個樣本某一位點所有allele的總深度
--hwe 2.5 計算Hardy-Weinberg p-value講到如何求p值，這個參數(shù)就是根據(jù)p值來過濾的，小于閾值則被過濾掉
--max-missing 前面已經(jīng)舉例；--max-missing-count 某個位點缺失樣本個數(shù)多于某個閾值則過濾掉
--phased 某個位點如果含有未定相的基因型則過濾
--minQ 根據(jù)vcf文件的QUAL列來過濾，比如
vcftools --gzvcf combined200.vcf.gz --minQ 100 --recode --out minQ_100
3.3 過濾樣本的選項
--indv 和 --remove-indv保留或去除哪一個樣本，后接樣本ID
--keep 和 --remove保留或去除哪些樣本，后接文件名，文件中一行一個樣本ID
3.4 輸出選項
等位基因統(tǒng)計
--freq 和 --counts用于統(tǒng)計等位基因的頻率和數(shù)量
深度統(tǒng)計
每個樣本所有位點的平均深度、單個位點所有樣本的深度加和（或平均）、每個樣本每個位點基因型的深度。感覺沒什么用
LD統(tǒng)計
Ti/Tv統(tǒng)計
NUCLEOTIDE DIVERGENCE STATISTICS
FST STATISTICS
OTHER STATISTICS
格式轉(zhuǎn)換
3.5 比較選項
reference
The variant call format and VCFtools: https://academic.oup.com/bioinformatics/article/27/15/2156/402296
VCFtools--A set of tools written in Perl and C++ for working with VCF files:
https://vcftools.github.io/man_latest.html#EXAMPLES
