☰
VCFtools实战指南:从VCF过滤、统计到群体遗传学分析
2026/9/29 8:07:45 网站建设 项目流程

做群体遗传学或者重测序数据分析的人,大概率都有过这样一个时刻:比对、变异检测一路跑下来,终于拿到一个几十GB的VCF文件,想筛一筛高质量的位点,算一算样本缺失率,结果对着文件里的几十列信息发懵,写awk脚本处理几百万行又慢又容易翻车。我最早遇到VCFtools是在一个水稻群体重测序项目里,当时需要从SNP calling结果里按深度、质量、缺失率、次等位基因频率做一套过滤,再转成下游软件需要的格式,VCFtools帮了大忙,而且一套参数拿走就能复用到别的项目,省下的时间相当可观。

这篇内容不打算写成官方文档的搬运工,而是从实际使用角度把VCFtools的安装、基础过滤、常用统计和典型报错讲清楚。哪怕你对VCF格式还不太熟悉,只要会敲Linux命令,跟着操作就能跑起来。适合刚入门的生信学生,也适合想系统梳理VCFtools常用功能的研究人员。

1. VCFtools到底解决了什么问题,装之前先想明白

1.1 我是在什么场景下第一次用到它的

我记得很清楚,当时手头是一个包含120个样本的全基因组重测序数据,经过GATK HaplotypeCaller之后产出了原始的VCF文件,里面大概有1800多万个变异位点,文件体积接近40GB。这个规模直接用文本工具去处理是不现实的,而且我需要的不是简单的“按列筛一筛”,而是一套组合条件:每个位点的质量值、每个样本的测序深度、缺失率、等位基因频率都要同时满足要求,最好还能顺手把基因型转成0/1/2的矩阵,方便后续做PCA。VCFtools是我当时找了一圈后觉得最合适的选择。

后来陆陆续续接触了bcftools、PLINK、vcflib这些工具,但VCFtools仍然是我做“快速质检+常规过滤”时的首选。它的核心定位一句话就能讲清楚:它是一个专门针对VCF/VCF.gz文件做过滤、统计和格式转换的命令行工具集,目标用户就是做群体遗传学、进化生物学和分子育种相关分析的人。

1.2 工具定位:过滤、统计、格式转换三件事

VCFtools从诞生到现在已经十几年了,底层用C++实现,支持流式读取,能处理比较大的VCF文件。它做的事情可以概括成三大类:

  • 过滤:按位点质量(--minQ)、测序深度(--minDP/--maxDP)、基因型质量(--minGQ)、缺失率(--max-missing)、等位基因频率(--maf/--min-allele-count)、双等位/多等位(--min-alleles/--max-alleles)等条件筛掉不合格的变异。
  • 统计:计算等位基因频率、样本和位点缺失率、杂合度、亲缘关系指数、连锁不平衡等,输出结果通常是文本表格,方便导入Excel或R做进一步分析。
  • 格式转换:在VCF、PLINK、012矩阵、BEAST、Structure、Ped等格式之间来回切换,这个功能在实际项目中极其常用。

我自己的体会是,VCFtools最能帮上忙的阶段是“拿到VCF之后、做下游分析之前”。这个阶段看起来不起眼,但如果处理不当,后面所有的群体结构、选择压力分析都会受影响。

1.3 和bcftools、PLINK这些工具的分工

经常会有人问:已经有bcftools了,为什么还要用VCFtools?两者确实有功能重叠,但侧重点不一样。bcftools在速度和大数据量处理上有优势,它对VCF规范的支持也更现代,比如能直接处理多等位位点、能灵活操作INFO字段。而VCFtools的统计功能更“开箱即用”,很多遗传学分析需要的统计量,比如--relatedness、--het、--site-quality,在VCFtools里一条命令就能出结果,不用自己写复杂表达式。PLINK则更偏向于关联分析和群体分层场景,它的强项在格式转换和LD计算,但在VCF的精细过滤上不如VCFtools灵活。

所以我自己现在的习惯是:批量大、流程化操作优先用bcftools,快速质检和统计出报告用VCFtools,要做GWAS或IBD分析再转PLINK。它们不是替代关系,是配合关系。这篇内容先聚焦VCFtools,等后面有机会再单独聊聊bcftools的使用技巧。

2. 安装的三种姿势与亲测踩坑记录

2.1 conda安装:实测最省心的一种

如果你是生信环境重度用户,我强烈建议直接用conda。它会连依赖一起装好,基本不会出现“装完了运行报错找不到某个库”的问题。要装最新版或者指定版本都很方便:

conda install -c bioconda -c conda-forge vcftools

如果想快一点,可以先把conda换成mamba,再用mamba安装:

mamba install -c bioconda -c conda-forge vcftools

装完之后验证一下:

vcftools --version

正常会输出版本号,比如VCFtools (0.1.16)。我在三台服务器上用conda装过VCFtools,没有一次失败的,唯一需要留意的是base环境和项目环境的隔离。建议单独为生信分析建一个环境,比如conda create -n bioinfo,再在这个环境里装VCFtools,避免几十个包之间互相抢依赖。

2.2 apt安装:快但版本可能偏老

Ubuntu和Debian系的用户可以直接用apt安装:

sudo apt-get update sudo apt-get install vcftools

这个方式胜在简单,装完就能用,系统会把vcftools、vcf-validator、vcf-stats这些辅助脚本一起装好。但我要提醒一点:apt仓库里的VCFtools版本往往比官方GitHub上的滞后不少。倒不是说老版本不能用,基础的过滤和统计都支持,但当你想用一些后来新增的功能时,可能就会因为版本问题报“unrecognized option”。

另外,apt install装的是系统全局环境,如果你在服务器上没有root权限,这条命令就行不通了,还是老老实实用conda或者源码编译。

2.3 源码编译:最折腾,但能保证功能完整

源码编译适合两种情况:一是需要最新开发版的功能,二是服务器网络环境差、conda和apt都用不了。VCFtools的源码托管在GitHub上,编译过程不算复杂,但依赖需要提前装全:

git clone https://github.com/vcftools/vcftools.git cd vcftools ./autogen.sh ./configure make sudo make install

从个人经验看,编译过程最容易出问题的是缺autoconf、automake、libtool和zlib的开发包。Debian系可以先补齐这些依赖:

sudo apt-get install autoconf automake make g++ zlib1g-dev libpcre3-dev

如果要安装到用户目录而不是系统目录,在./configure时指定prefix:

./configure --prefix=$HOME/software/vcftools make && make install

然后在~/.bashrc里追加环境变量:

export PATH=$HOME/software/vcftools/bin:$PATH export PERL5LIB=$HOME/software/vcftools/lib/perl5:$PERL5LIB

这里要提醒一下PERL5LIB,因为VCFtools附带了很多Perl写的辅助脚本,这些脚本在运行时会调用它的Perl模块,如果不把lib路径加进去,后面跑vcf-validator时很容易报“Can't locate Vcf.pm in @INC”的错。我第一次编译装到自定义目录时就栽在这个上面,后来找到原因之后,每次装都记得顺手把PERL5LIB配上。

2.4 macOS和Windows环境怎么处理

macOS用户建议优先用Homebrew:

brew install vcftools

也可以走conda,两者都试过,没遇到什么坑。Windows上就比较尴尬了,VCFtools官方并没有原生Windows版,我一般不建议把时间花在折腾MSYS2或Cygwin上,最省力的路线是装WSL2,在WSL的Ubuntu环境里走一遍Linux安装流程,后面所有命令都按Linux思路跑,体验和服务器上完全一致。如果只是简单处理小文件,Windows里也可以用Anaconda Prompt兼容层试试装conda版,但文件大了性能不太行,还是WSL2最靠谱。

3. 第一次上手:读文件、看内容、跑通基础过滤

3.1 准备一个小型测试VCF,别一上来就压全基因组

初学VCFtools最忌拿着全基因组级别的文件直接试,等程序跑个几分钟半天不知道对不对,心态很容易崩。建议先从真实数据里截一个区域,或者干脆手写一个微型VCF来练手。

手写是最快的,几行就够:

##fileformat=VCFv4.2 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S1 S2 S3 chr1 100 . A G 50 PASS . GT:DP 0/1:12 1/1:25 0/0:8 chr1 200 . C T 20 q10 . GT:DP 1/1:30 0/0:10 ./.:. chr1 300 . G A 80 PASS . GT:DP 0/1:18 0/1:22 1/1:31

保存成test.vcf备用。如果你有真实数据,也可以直接从大VCF里截取一段:

bcftools view big.vcf.gz chr1:1000000-1010000 > test.vcf

或者用VCFtools自己的--min-alleles配合--chr拿一个染色体的部分位点,本质上都是让我的测试文件足够小,方便观察每条命令的输出变化。

3.2 先确认文件能被正常读入

装好VCFtools后的第一个动作,不是急着过滤,而是确认它能正确读取文件。最简单的方式是跑一个统计命令:

vcftools --vcf test.vcf --out test_basic --missing

跑完会看到终端输出一段日志,里面有“Parameters as interpreted”和样本数、位点数的汇总信息。更重要的是,目录会生成test_basic.imiss和test_basic.lmiss两个文件,一行一行看过去,能直观看到每个样本的缺失率、每个位点的缺失率。如果文件读取失败,日志里会在最显眼的位置写清楚错误原因。

VCFtools还会在每个输出文件名后面自动加后缀,比如--missing会生成.imiss和.lmiss,--freq生成.frq,--recode生成.recode.vcf。搞清楚这个命名规则,能少踩很多“我输出的文件去哪了”的坑。

3.3 第一套过滤命令:从QC到MAF

数据读取正常后,就可以按项目需求做过滤了。我这里给出一套非常典型的过滤命令,它在很多群体遗传学分析里可以作为通用起点:

vcftools --vcf test.vcf \ --minQ 30 \ --minDP 5 \ --max-missing 0.8 \ --maf 0.05 \ --min-alleles 2 \ --max-alleles 2 \ --recode \ --recode-INFO-all \ --out test_filtered

逐条解释一下:

  • --minQ 30:保留QUAL值不低于30的位点,对应Phred质量值,相当于错误率低于千分之一。
  • --minDP 5:去掉测序深度低于5x的基因型。注意这个参数作用于单个样本的基因型水平,低深度的基因型会被标记为缺失,它不会直接删除整个位点。
  • --max-missing 0.8:位点缺失率不超过0.8,也就是至少有80%的样本在该位点有非缺失基因型。这个参数在过滤群体数据时特别有用,能剔除只在极少数样本里出现的位点。
  • --maf 0.05:只保留最小等位基因频率不低于5%的位点,这个条件能过滤掉大量低频稀有变异,减少下游分析的噪音。
  • --min-alleles 2 --max-alleles 2:只保留双等位基因位点。很多下游分析软件并不支持多等位位点,这个条件在实战中几乎必加。
  • --recode:生成过滤后的VCF文件,这是“保留结果”的开关。
  • --recode-INFO-all:把原始INFO列里的注释信息也一起保留,如果省略,INFO列大部分内容会被丢弃。

跑完之后,test_filtered.recode.vcf就是过滤后的新文件。我建议每次跑完先看日志里剩余位点数,判断过滤条件是否过严或过松。如果原来100万个位点过滤完只剩1万个,那大概率是--max-missing或者--maf定得太苛刻了。

3.4 深度解析几个高频过滤参数

VCFtools的过滤参数很多,刚接触容易看得眼花缭乱。我根据自己的使用频率,把最有用的几个挑出来单独说说。

--minQ 和 --minGQ 的区别,这是个容易混淆的点。--minQ过滤的是VCF中QUAL列的整体位点质量值,作用于所有样本;--minGQ过滤的是每个样本基因型质量(GQ,Genotype Quality),是逐样本判断的。实际中我一般两个都会用,位点层面用--minQ,样本基因型层面用--minGQ,该设多少要根据你的测序深度和变异检测软件来定,没有绝对的通用值。

--minDP 和 --maxDP,控制的是深度。深度过低可能是假阳性,深度过高往往意味着重复区域或拷贝数变异区域,这些区域的变异可靠性也不高。我在人类全外显子组数据分析里常用--minDP 8 --maxDP 100,这只是经验值,具体情况要结合测序深度分布调整。有一个好习惯是先用--depth统计一下样本深度分布,再定阈值。

--remove-filtered-all,这个参数作用是移除所有FILTER列不是PASS的位点。大家在跑GATK时,通常会在joint calling之后对$FILTER列做一些标注,比如关于链特异性、HaplotypeScore等注释可能被放进FILTER列。如果在过滤时想彻底避开这部分位点,就加上--remove-filtered-all。但要注意,如果你的原始VCF里FILTER列是空的或不是PASS,也会被它丢掉,处理前最好先看看FILTER列的分布。

--keep 和 --remove,这两个参数直接通过样本列表文件来筛选样本。文件里每行一个样本名。比如:

vcftools --vcf test.vcf --keep samples_keep.txt --recode --out test_subset

在处理多个分组比较的场景里,我先用--keep把每个组的样本拆出来,再分别做频率统计,比在原始文件里反复操作省事很多。

4. 简单统计:用VCFtools让样本和位点质量现出原形

4.1 等位基因频率统计,一张表看清位点分布

过滤完数据后,我最先跑的统计命令基本都是频率。

vcftools --vcf test_filtered.recode.vcf --freq --out final_freq

这条命令会生成final_freq.frq文件,里面包含CHROM、POS、N_ALLELES、N_CHR以及每个等位基因的频率。特别说明一下N_CHR,它表示该位点上实际覆盖到的染色体数量,等于样本数乘以2再减去缺失基因型的数目。如果某位点缺失率高,N_CHR会明显低于其它位点,这个值在后续算观测杂合度时非常关键。

如果VCFtools在处理时遇到多等位基因位点,频率输出表里会多出几列。要是下游分析只需要双等位位点,就在这一步加过滤条件,不要等统计完再人工去筛,非常麻烦。

4.2 样本缺失率与位点缺失率,两兄弟要一起看

缺失率统计是我每一次拿到新数据都会跑的命令:

vcftools --vcf test_filtered.recode.vcf --missing --out qc_summary

生成的两个文件,一个叫qc_summary.imiss,是每个样本的缺失基因型比例;另一个叫qc_summary.lmiss,是每个位点在不同样本间的缺失比例。样本缺失率可以用来判断测序质量,比如某个样本的缺失率高达30%,其他样本都只有5%,那这个样本测序过程大概率出了问题,要么DNA质量差,要么测序深度波动剧烈。位点缺失率则反映了这个位点在多少样本中能可靠检测出来,对于那些只在少数样本中存在的位点,后续做群体结构分析时很容易产生分组假象。

实际分析中这两个文件可以直接导进R或者Excel做分布图。我自己习惯看缺失率直方图,如果样本缺失率呈现明显的长尾分布,尾部那几个样本要特别警惕。曾经有个项目,一个样本缺失率0.28,其他样本都在0.1以下,查了测序记录才发现这个样本的文库测序数据量只有平均水平的一半。

4.3 杂合度、亲缘关系:初步的样本质量控制

除了频率和缺失率,VCFtools还能快速算杂合度和样本间亲缘关系,这些在检查样本污染或者意外重复时很有用。

vcftools --vcf test_filtered.recode.vcf --het --out het_check

het_check.het文件里,O(HOM)是观测纯合数,E(HOM)是期望纯合数,N_SITES是参与计算的位点总数。如果一个样本的观测纯合数明显高于期望,说明这个样本可能存在近亲繁殖或者群体分层;反过来,如果明显低于期望,可能是这个样本混入了其他个体或者存在污染。单独看一两个数字可能没感觉,把所有样本画成分布图,异常点就一目了然了。

再看亲缘关系:

vcftools --vcf test_filtered.recode.vcf --relatedness --out relatedness

输出的.relatedness矩阵可以告诉你样本间的亲缘系数。全同胞或亲子关系的系数会接近0.5,半同胞接近0.25,无关个体接近0。如果一对样本明明标记为不同个体,亲缘系数却接近0.5,多半是样本管理环节出了乌龙。这个检查在群体遗传学项目里属于标准动作。

5. 实战翻车现场:这些年我在VCFtools上踩过的坑

5.1 压缩格式没认出来,输出一堆乱码报错

有一次我把一个xxx.vcf.gz文件传到了新服务器,忘了这是压缩格式,直接写了--vcf xxx.vcf.gz,结果程序读了一小段就在屏幕上喷出大段乱码然后报错退出。后来才意识到VCFtools读取压缩文件必须显式指定--gzvcf参数:

vcftools --gzvcf xxx.vcf.gz --freq --out result

这个参数看起来只是多了个gz,但忘了加就会让程序去按纯文本方式解析二进制GZIP数据,百分百翻车。所以说,拿到一个VCF文件,第一件事应该是检查后缀,不能想当然。

5.2 多等位基因位点导致统计结果异常

还有一次是在跑--freq的时候,发现某个染色体上频率表有很多行格式乱掉,仔细一看,原来是那个文件里含有大量三等位和四等位位点。VCFtools虽然能读入这些位点,但某些统计在输出上并不会按照你的直觉组织。比如--freq表会列出多个ALLELE列,但后续用R处理的人往往默认只有两个等位基因,解析时就会出错。

我现在的习惯是:在做任何统计之前,先把多等位位点过滤掉:

vcftools --vcf input.vcf --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out biallelic

跑统计用biallelic.recode.vcf,这一步之后基本不会再遇到等位基因数量相关的解析问题。虽然有一些分析场景需要保留多等位位点,但建议单独建一个文件处理,不要在原始文件上频繁操作。

5.3 输出文件太多,对应关系搞不清楚

VCFtools每个参数会生成不同后缀的文件,刚上手时很容易搞晕。这里列一个我常用的对应表,方便查阅:

参数输出后缀内容说明
--missing.imiss/.lmiss样本缺失率 / 位点缺失率
--freq.frq等位基因频率
--het.het样本观测与期望纯合数
--relatedness.relatedness样本间亲缘关系矩阵
--depth.idepth/.ldepth样本平均深度 / 位点平均深度
--recode.recode.vcf过滤后的VCF文件
--012.012/.012.indv/.012.pos基因型0/1/2矩阵

拿--012举个实际例子,很多下游工具需要基因型数字矩阵,不需要额外写脚本解析VCF,直接用这个参数就能生成三件套:.012是真正的矩阵,.012.indv是样本名列表,.012.pos是位点坐标列表。这个功能在画PCA或计算个体间遗传距离时很常用,也是我最初接触VCFtools时觉得最惊喜的功能之一。

5.4 大文件处理太慢的应对思路

VCFtools是单线程程序,遇到几千万位点的全基因组VCF,一次--recode可能跑上很久。如果只是做质量检查,可以先用--chr限定一条染色体,比如:

vcftools --gzvcf whole_genome.vcf.gz --chr chr1 --freq --out chr1_freq

如果确实需要全基因组范围过滤,我一般会在VCFtools之前先用bcftools按染色体并行预处理,合并结果后再回到VCFtools做统一统计,时间能快不少。另一个思路是分样本处理:先用--keep拆出部分样本跑通流程,确认参数无误后再全量跑。这个步骤虽然看起来多此一举,但在大项目里能省掉大量反复调试的时间。

5.5 路径、权限与临时文件目录的小问题

还有一个特别容易忽略但实际很常见的坑:--out指定的路径目录不存在。VCFtools遇到这种情况不会自动创建目录,而是直接报错退出。所以每次写输出路径前,我先mkdir -p把目录建好。另外,把输出写到/tmp或者网络挂载盘有时也会遇到权限或IO慢的问题,建议在本地磁盘工作目录下操作,跑完再把结果同步走。

如果同时跑多个VCFtools任务,/tmp下会因为临时文件冲突导致任务失败,最稳妥的做法是给每个任务指定独立的--out前缀,而且尽量不要让两个任务共用同一个输出目录。

6. 写在最后:我对VCFtools的一些使用心得

用VCFtools这么多年,我越来越觉得它是一个典型的“老而弥坚”型工具:界面不炫,速度不算顶尖,但统计功能覆盖全面、参数设计直觉、结果稳定可靠。日常项目中我把它当成第一道质检关卡,样本一到手先跑缺失率、深度、杂合度,看数据顺不顺眼,能过滤的先过滤,然后才轮到PLINK、Structure或者R出场。

想给新接触的朋友一个建议:不要试图一次记住所有参数。先记住--vcf、--recode、--out、--maf、--max-missing这几个核心用法,能解决工作中80%的场景。其余参数用到再查,查完跑个小样本验证,多来几次自然就熟练了。另外就是处理正式项目数据前,一定先用小文件试跑一遍,哪怕宁可多花几分钟做测试,也比直接拿全量数据跑完才发现过滤条件写错了要省时间得多。

最后分享一个我常用的组合用法:先用VCFtools做完整QC并生成imiss和ldepth,再用R画分布图,把异常样本挑出来后用--keep生成新VCF,继续往下游分析。这套流程不依赖大型流程管理工具,一条条命令敲下来,思路非常清楚。VCFtools作为这个流程的“入口关卡”,一直很可靠。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询