vcf文件与vcftools（二）
vcftools是用来处理vcf和bcf文件的工具集，功能涵盖了过滤，数据格式转换，一些指标的统计计算，vcf文件之间变异信息的差异比较等等。
1.基本用法
vcftools [ --vcf FILE | --gzvcf FILE | --bcf FILE] \
[ --out OUTPUT PREFIX ] [ FILTERING OPTIONS ] [ OUTPUT OPTIONS ]
2.使用
官网提供了几个常用示例，下面用我的vcf文件测试一下，我的vcf文件共249个样本，27012个SNP标记
2.1 对来自chr1的每一个位点统计其基因频率
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
第三列表示在所有样本中，该位点出现了几种碱基，大多数情况下为2，也可能是3、4；第四列表示该位点出现的碱基总数，一个样本贡献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就需要重新计算。
所以过滤之后默认情况下是不含有INFO列的。
该选项表示将原来vcf文件中所有INFO信息保留下来。
--recode-INFO <string>
比如--recode-INFO AC表示只将AC保留
--keep-only-indels和--remove-indels作用相反
根据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 结合管道操作符
用--stdout代替--out即可结合管道进行其他处理
2.5 计算Hardy-Weinberg p-value
Hardy-Weinberg遗传平衡检验用到的统计方法是卡方测验，详细过程参考：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.其它参数说明
3.1 基本选项
--vcf | --gzvcf | --bcf 输入格式，三选一即可
--out <输出文件前缀>
--stdout 和 -c 输出到标准输出，能通过管道与其它命令结合
3.2 过滤SNP的选项
POSITION FILTERING
可以提供一些位置信息，或是能表示位置的文件，具体参数见官方文档
ALLELE FILTERING
--maf 和 --max-maf用来限定最小等位基因频率(MAF)的范围
--mac 和 --max-mac和上面类似，用来限定最小等位基因的数量
--min-alleles 和 --max-alleles用来限定等位基因的数量，比如这条命令只会保留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值，这个参数就是根据p值来过滤的，小于阈值则被过滤掉
--max-missing 前面已经举例；--max-missing-count 某个位点缺失样本个数多于某个阈值则过滤掉
--phased 某个位点如果含有未定相的基因型则过滤
--minQ 根据vcf文件的QUAL列来过滤，比如
vcftools --gzvcf combined200.vcf.gz --minQ 100 --recode --out minQ_100
3.3 过滤样本的选项
--indv 和 --remove-indv保留或去除哪一个样本，后接样本ID
--keep 和 --remove保留或去除哪些样本，后接文件名，文件中一行一个样本ID
3.4 输出选项
等位基因统计
--freq 和 --counts用于统计等位基因的频率和数量
深度统计
每个样本所有位点的平均深度、单个位点所有样本的深度加和（或平均）、每个样本每个位点基因型的深度。感觉没什么用
LD统计
Ti/Tv统计
NUCLEOTIDE DIVERGENCE STATISTICS
FST STATISTICS
OTHER STATISTICS
格式转换
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
