做群体遗传学或者重测序数据分析的人大概率都有过这样一个时刻比对、变异检测一路跑下来终于拿到一个几十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时很容易报“Cant locate Vcf.pm in INC”的错。我第一次编译装到自定义目录时就栽在这个上面后来找到原因之后每次装都记得顺手把PERL5LIB配上。2.4 macOS和Windows环境怎么处理macOS用户建议优先用Homebrewbrew install vcftools也可以走conda两者都试过没遇到什么坑。Windows上就比较尴尬了VCFtools官方并没有原生Windows版我一般不建议把时间花在折腾MSYS2或Cygwin上最省力的路线是装WSL2在WSL的Ubuntu环境里走一遍Linux安装流程后面所有命令都按Linux思路跑体验和服务器上完全一致。如果只是简单处理小文件Windows里也可以用Anaconda Prompt兼容层试试装conda版但文件大了性能不太行还是WSL2最靠谱。3. 第一次上手读文件、看内容、跑通基础过滤3.1 准备一个小型测试VCF别一上来就压全基因组初学VCFtools最忌拿着全基因组级别的文件直接试等程序跑个几分钟半天不知道对不对心态很容易崩。建议先从真实数据里截一个区域或者干脆手写一个微型VCF来练手。手写是最快的几行就够##fileformatVCFv4.2 ##FORMATIDGT,Number1,TypeString,DescriptionGenotype ##FORMATIDDP,Number1,TypeInteger,DescriptionRead 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过滤的是每个样本基因型质量GQGenotype 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_checkhet_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作为这个流程的“入口关卡”一直很可靠。
VCFtools实战指南：从VCF过滤、统计到群体遗传学分析
做群体遗传学或者重测序数据分析的人大概率都有过这样一个时刻比对、变异检测一路跑下来终于拿到一个几十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时很容易报“Cant locate Vcf.pm in INC”的错。我第一次编译装到自定义目录时就栽在这个上面后来找到原因之后每次装都记得顺手把PERL5LIB配上。2.4 macOS和Windows环境怎么处理macOS用户建议优先用Homebrewbrew install vcftools也可以走conda两者都试过没遇到什么坑。Windows上就比较尴尬了VCFtools官方并没有原生Windows版我一般不建议把时间花在折腾MSYS2或Cygwin上最省力的路线是装WSL2在WSL的Ubuntu环境里走一遍Linux安装流程后面所有命令都按Linux思路跑体验和服务器上完全一致。如果只是简单处理小文件Windows里也可以用Anaconda Prompt兼容层试试装conda版但文件大了性能不太行还是WSL2最靠谱。3. 第一次上手读文件、看内容、跑通基础过滤3.1 准备一个小型测试VCF别一上来就压全基因组初学VCFtools最忌拿着全基因组级别的文件直接试等程序跑个几分钟半天不知道对不对心态很容易崩。建议先从真实数据里截一个区域或者干脆手写一个微型VCF来练手。手写是最快的几行就够##fileformatVCFv4.2 ##FORMATIDGT,Number1,TypeString,DescriptionGenotype ##FORMATIDDP,Number1,TypeInteger,DescriptionRead 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过滤的是每个样本基因型质量GQGenotype 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_checkhet_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作为这个流程的“入口关卡”一直很可靠。
