01 背景
VCFtools
VCFtools 是一套用于处理 VCF 文件(如 1000 基因组计划生成的文件)的工具,主要用 Perl 和 C++ 编写。
02 参考
https://github.com/vcftools/vcftools #官网
03 安装
构建 VCFtools
关于构建过程中的配置步骤的常规帮助,可以通过以下命令获取:
./configure --help
从发布的 Tarball 构建
./configure
make
make install
你可能需要使用 sudo 权限来运行 make install
。
从 GitHub 构建
wget -c https://github.com/vcftools/vcftools/releases/download/v0.1.16/vcftools-0.1.16.tar.gz
wget -c https://github.com/vcftools/vcftools/archive/refs/tags/v0.1.16.tar.gz
git clone https://github.com/vcftools/vcftools.git
cd vcftools
./autogen.sh
./configure
make
make install
你可能需要使用 sudo 权限来运行 make install
。
04 使用
文档和使用示例可以在这里找到:
VCFtools
手册页面也可以访问。如果将前缀设置为 /usr
或者 MANPATH 指向 $prefix/share/man
,你可以通过以下命令访问手册页面:
man vcftools
获取帮助
关于 VCFtools 的最佳帮助方式是通过电子邮件向邮件列表发送问题:
vcftools-help@lists.sourceforge.net
05 常用命令行
优化变异筛选流程:
在变异筛选中,保留了以下条件的变异:至少50%的个体成功进行了 SNP calling,最低质量分数为30,次要等位基因计数(MAC)为3的位点。
第一步:初步过滤
vcftools --gzvcf raw.vcf.gz --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out raw.g5mac3
此命令的作用如下:
--gzvcf raw.vcf.gz
:指定输入的压缩 VCF 文件。--max-missing 0.5
:过滤掉基因型缺失率超过50%的位点。--mac 3
:过滤掉次要等位基因计数小于3的SNP。--minQ 30
:仅保留质量分数大于30的位点。--recode
:输出过滤后的新VCF文件。--recode-INFO-all
:保留原VCF文件中的所有INFO字段。--out raw.g5mac3
:指定输出文件的前缀。
运行此命令后,生成的输出会显示:
- 所有40个个体均被保留下来。
- 经过过滤,147540个变异位点中,78434个被保留下来(约一半的位点被过滤掉)。
- 运行时间为17秒。
第二步:根据深度进一步过滤
继续根据最小基因型深度进行过滤:
vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3
该命令执行后,输出结果如下:
- 所有40个个体保留。
- 78434个位点全部被保留,说明这些位点都满足最小深度(大于3)的条件。
- 运行时间为14秒。
关于最小深度的讨论:
- 设置最小深度为3是为了确保变异的可靠性。虽然理论上深度越高越好,但过高的深度可能会引入误差,因此这个标准在大多数情况下是合理的。
第三步:评估低深度带来的潜在错误
为了进一步检查低深度可能导致的基因分型错误,使用以下脚本进行分析:
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/ErrorCount.sh
chmod +x ErrorCount.sh
./ErrorCount.sh raw.g5mac3dp3.recode.vcf
运行结果显示:
- 对于只有1到3个读取的变异,潜在的基因分型错误从15986到53714之间不等。
- 总体的错误率在0.010至0.043之间。
根据这些结果,最低深度为3的位点的错误率小于5%,表明该深度是合适的。
第四步:筛除数据丢失过多的个体
接下来,我们基于个体间数据丢失水平来进行筛选:
vcftools --vcf raw.g5mac3dp3.recode.vcf --missing-indv
此命令生成的 out.imiss
文件显示每个个体的缺失数据百分比。对于缺失率高于50%的个体,我们可以通过以下步骤进行筛选:
mawk '$5 > 0.5' out.imiss | cut -f1 > lowDP.indv
vcftools --vcf raw.g5mac3dp3.recode.vcf --remove lowDP.indv --recode --recode-INFO-all --out raw.g5mac3dplm
此命令将移除那些缺失数据过多的个体。
第五步:进一步过滤基因座的变异
在此步骤中,进一步根据变异的缺失率和次要等位基因频率(MAF)进行筛选:
vcftools --vcf raw.g5mac3dplm.recode.vcf --max-missing 0.95 --maf 0.05 --recode --recode-INFO-all --out DP3g95maf05 --min-meanDP 20
--max-missing 0.95
:要求每个基因座至少在95%的个体中有数据。--maf 0.05
:只保留次要等位基因频率大于5%的变异。--min-meanDP 20
:要求基因型的平均深度大于20。
第六步:按群体筛选变异
若数据来自多个采样地点,可以根据不同群体的缺失数据进行筛选。首先创建群体信息文件(popmap
):
cat popmap
接下来,创建每个群体的个体列表并分别计算缺失数据:
mawk '$2 == "BR"' popmap > 1.keep
mawk '$2 == "WL"' popmap > 2.keep
vcftools --vcf DP3g95maf05.recode.vcf --keep 1.keep --missing-site --out 1
vcftools --vcf DP3g95maf05.recode.vcf --keep 2.keep --missing-site --out 2
合并缺失数据文件并按10%缺失数据筛选变异:
cat 1.lmiss 2.lmiss | mawk '!/CHR/' | mawk '$6 > 0.1' | cut -f1,2 >> badloci
vcftools --vcf DP3g95maf05.recode.vcf --exclude-positions badloci --recode --recode-INFO-all --out DP3g95p5maf05
以上流程确保了通过多步骤过滤后的高质量SNP集,适合后续的分析。
06 引用
The Variant Call Format and VCFtools, Petr Danecek, Adam Auton, Goncalo Abecasis, Cornelis A. Albers, Eric Banks, Mark A. DePristo, Robert Handsaker, Gerton Lunter, Gabor Marth, Stephen T. Sherry, Gilean McVean, Richard Durbin and 1000 Genomes Project Analysis Group, Bioinformatics, 2011 http://dx.doi.org/10.1093/bioinformatics/btr330