1. 先想清楚过滤目标再谈工具组合说句实在话从测序公司、公共数据库或者自己服务器上拿到的 VCF 文件几乎没有一份是干净的。GATK、DeepVariant、Sentieon 各家流程吐出来的原始 VCF 里低深度位点、测序错误堆出来的假阳性、只在个别样本里出现的基因型全都混在一起。直接拿去跑群体结构、GWAS 或者亲缘关系分析结果大概率没法看。我这几年处理过的 VCF 少说也有几百份最后稳定下来的一套组合就是bcftools vcftools前者负责高速的流式切片、归一化和重格式化后者负责统计量充足、参数直观的群体级位点过滤。两者串起来一个中等规模的 VCF 十几分钟内就能从毛坯变成可用件。这套流程适合谁如果你手上有几十到几千个样本的 SNP/INDEL 数据集准备做群体遗传、关联分析、进化树或者 PCA那基本都能用。如果你是刚接触生信的新人只要会敲命令行、能看懂 VCF 的表头也能照着往下走。我不打算把参数堆成一份手册而是把每个参数为什么定成这个值讲清楚——因为真正让人踩坑的从来不是命令写错而是参数拍脑袋。1.1 VCF 里的脏数据到底长什么样先把敌人看清楚。VCF 文件里的噪声大致分三层处理手段也完全不同。第一层是位点级问题整个位点所有样本都不靠谱。典型表现是 QUAL 很低比如测序深度只有 2 条 read 支持、INFO 里的 DP 极低或极高极高通常是重复序列或者比对堆叠区域、大量样本在这个位点上是./.缺失率过高。这类问题的处理方式是整条记录剔除。第二层是基因型级问题位点本身没问题但某个样本在这个位点上的基因型不可信。比如某个样本 DP1就算叫出来0/1也基本是噪声再比如 GQ基因型质量只有 3说明 caller 自己都不确定。这类问题的处理方式是把该样本该位点的基因型置为缺失而不是把整个位点删掉——这一点非常关键很多人一上来就把整条记录删了最后样本间共有的位点数掉得惨不忍睹。第三层是群体统计层面的问题位点本身数据完整但从群体角度看有偏。比如次等位基因频率MAF过低在几百个样本里只出现一两次这种位点对群体分析几乎没有信息量反而会引入噪声再比如严重偏离哈迪-温伯格平衡HWE的位点往往暗示着比对错误或者分型错误。这类问题只有在样本集合确定之后才能计算顺序绝对不能搞反。搞清楚这三层后面的参数选择就有了依据。我见过太多人把三层问题混在一起用一条命令解决最后既不知道删掉了什么也不知道为什么删。1.2 bcftools 与 vcftools 的分工为什么这样切有人会问bcftools 和 vcftools 功能大量重叠为什么还要两个一起用只用其中一个不行吗答案是能用但各有各的别扭。bcftools 的核心优势是流式和内存效率。它基于 htslib读写 BCF/VCF 走的是管道式处理几百 GB 的文件也能用很低的常驻内存跑完-i/-e表达式系统表达能力极强可以同时写 QUAL、DP、FMT/DP、MAF、F_MISSING 各种条件还支持-Ou直接输出未压缩 BCF 接到下一个命令省掉中间文件的磁盘 IO。做归一化bcftools norm、多等位拆分-m -any、类型筛选-v snps这些活它是目前最好用的工具没有之一。vcftools 的优势是群体统计参数的直觉性和输出的完整性。它内置了--max-missing、--maf、--hwe、--min-meanDP、--minDP、--maxDP这一整套现成的开关参数名一看就懂还顺带输出.log文件记录每一步过滤前后的位点数。做群体遗传分析的人最需要的过滤前后还剩多少位点、平均深度是多少、缺失率分布如何这些信息vcftools 给你的东西比 bcftools 更顺手。我的实际做法是这样切的环节工具理由格式转换、压缩、索引bcftools / bgzip / tabix快支持-Ou管道参考基因组左对齐、多等位拆分bcftools normvcftools 没有等价功能SNP/INDEL 类型筛选bcftools view -v一条参数搞定语义清晰位点级硬过滤QUAL/DPbcftools view -i表达式灵活速度快基因型级过滤DP/GQ 置缺失vcftools --minDP/--maxDP语义就是置缺失不删位点缺失率、MAF、HWE 过滤vcftools参数直观日志完整需要复杂组合条件时bcftools fill-tags view -i一次到位避免多轮 IO这套分工不是教条。如果你只做小规模数据几百个位点、几十个样本两个工具随便挑一个都行但当数据量上到百万位点、上千样本上面这个分工能帮你省掉大量时间和内存。注意vcftools 在处理超大 VCF 时会把不少信息读进内存我遇到过 200 万位点、2000 样本的数据直接把 64G 内存吃满的情况。这种规模建议先用 bcftools 把位点砍到百万以内再用 vcftools 做群体级过滤。1.3 过滤目标和参数先对齐动手之前先回答三个问题这三个问题的答案直接决定参数怎么定。第一这个数据集最终要干什么做 PCA 和群体结构MAF 阈值可以放到 0.05位点数砍到几万到几十万就够用做全基因组关联分析MAF 可以放宽到 0.01因为低频变异也可能有表型效应做系统发育树通常不需要 MAF 过滤但缺失率要求要更严。第二样本集合是什么这个最容易被忽略。等位基因频率、缺失率、HWE 全都是基于当前样本集合计算出来的。如果你打算只保留某一个亚群做分析那必须在过滤之前先把样本定下来否则算出来的 MAF 和 HWE 都是错的删掉的位点在子集里可能是完全合格的多态位点。第三能容忍丢掉多少位点我一般会在跑正式流程前先跑一遍保守参数看看剩多少再决定要不要收紧。不要一上来就用 QUAL≥50、MAF≥0.1、max-missing 0.99 这种组合最后可能只剩几千个位点做啥都不够。把这三点想清楚后面所有参数就都有了落脚点。接下来我按实际操作顺序把整个流程拆开讲。2. 动手前必须做的三件事环境、摸底、对齐很多人拿到 VCF 第一反应就是敲过滤命令这是最容易出岔子的地方。我自己的习惯是三步走装好工具确认版本、用 stats 给数据做体检、核对参考基因组与染色体命名。这三步加起来不到十分钟但能省掉后面几小时的排查。2.1 工具安装与版本确认最省事的方式还是 conda/mamba。我的环境一般是这么建的mamba create -n vcf clean -c conda-forge -c bioconda \ bcftools vcftools htslib tabix datamash plink2 mamba activate vcf为什么把 plink2 也装上因为后面做 LD 剪枝、PCA 的时候迟早要用装在一个环境里省得来回切。datamash是个小工具做中位数、均值统计特别方便配合bcftools query用起来很爽。装完第一件事是确认版本这非常重要bcftools --version vcftools --version为什么强调版本两个原因。一是bcftools view的-q/--min-afMAF 过滤参数、-i表达式里的FMT/前缀支持、fill-tags插件的可用性都是较新版本才有的老版本1.9 以前会直接报错。二是bcftools norm的行为在不同版本间有过调整尤其是-m的默认值和--check-ref的处理逻辑。我一般要求 bcftools ≥ 1.17、vcftools ≥ 0.1.16这两个版本用下来最稳。fill-tags是插件有些系统里默认没编译进去装完可以先测一下bcftools fill-tags 21 | head -5如果提示 Could not load plugin换 bioconda 的 bcftools 包基本就好了。2.2 用 bcftools stats 给数据做一次体检这一步是整个流程里我最看重的。不做体检就开始过滤等于闭着眼睛开车。bcftools stats -s - input.vcf.gz logs/raw.stats.txt bt-view raw.stats.txt # 如果有 plot-vcfstats 也可以直接出图输出的stats.txt分很多段每段以两三个字母开头我要重点看这几段SN 段Summary Numbers记录总记录数、SNP 数、INDEL 数、多等位位点数、无变异位点数、样本数。这里最重要的是多等位位点占比——如果超过 5%基本可以确定需要先做bcftools norm -m -any拆分。TSTV 段Transition/Transversion ratio全基因组水平的 Ts/Tv 对 SNP 来说正常在 2.0 到 2.1 之间人类。如果只有 1.5 甚至更低说明假阳性很多位点级过滤要收紧。QUAL 段QUAL 值的分布。如果大量位点 QUAL 在 20 以下那大概率是低深度区域的噪声。DP 段INFO/DP 的分布。这个直接决定你--maxDP该定多少。PSC 段Per-Sample Counts每个样本的位点数、SNP 数、杂合比例。如果某个样本的杂合比例明显偏离其他人比如别人 0.6它 0.3这个样本很可能有污染或者近交问题得单独处理。有一个特别实用的技巧是直接算深度中位数用来定maxDPbcftools query -f %INFO/DP\n input.vcf.gz \ | datamash median 1 mean 1假设输出中位数是 14、均值是 18那么合理的maxDP大概是中位数的 3 倍左右也就是 40~45。为什么是 3 倍因为超过 3 倍中位深度的位点绝大多数落在重复序列、旁系同源区域或者拷贝数变异区域这些地方 call 出来的基因型极不可靠。这个倍数不是硬规定如果数据是外显子捕获的深度分布更均匀2 倍就够了。还想看单样本深度分布的话bcftools query -f [%DP\n] input.vcf.gz | datamash median 1 mean 1 max 1注意[%DP]加了方括号表示取 FORMAT/DP 并且每个样本一行不加方括号的%DP取的是 INFO/DP这个区别我第一次用的时候栽过跟头取出来数不对还纳闷了半天。2.3 参考基因组与染色体命名一致性核对这一步是纯粹的防雷但雷一旦踩上就很烦。# 看 VCF 里声明的 contig bcftools view -h input.vcf.gz | grep ^##contig | head -30 # 看参考基因组的染色体列表 cut -f1,2 ref.fa.fai要核对两件事命名风格和长度。命名风格是指chr1还是1chrM还是MT。有些流程尤其从公共数据库下载的 VCF用的是1而你的 BED 文件或者参考基因组用的是chr1一旦不一致bcftools norm -f会报 reference allele mismatchvcftools --bed会静默地一个位点都不保留——这个静默失败是最坑的命令跑完了没报错结果文件是空的。长度也要对一下。同一条染色体在不同参考版本里长度差异可能很大比如 GRCh37 和 GRCh38 的 chr1 相差几十 Mb。如果 VCF 是 GRCh37 的而参考是 GRCh38左对齐和 REF 校验都会出问题。如果确实需要改命名最快的办法是bcftools annotate --rename-chrs chr_name_map.txt \ -Oz -o renamed.vcf.gz input.vcf.gzchr_name_map.txt是一个两列文件第一列是旧名第二列是新名比如1\tchr1。这个命令是纯字符串替换速度快也不改动坐标比重新比对靠谱得多。3. 位点级与基因型级过滤的参数怎么定工具和数据都准备好了接下来是核心环节。我把过滤按位点级 → 基因型级 → 群体统计级的顺序讲每层都给出参数的推导过程和常见取值。3.1 QUAL 与深度最容易拍脑袋拍错的两个参数QUAL是位点质量Phred 标度。QUAL30 意味着这个位点分型错误的概率大约是千分之一10^-3QUAL20 是百分之一QUAL50 是十万分之一。做群体研究我一般用 30做临床相关的分析会提到 50。但这里有个大坑有些 VCF 的 QUAL 字段是.。某些流程尤其是只输出基因型的流程或者经过某些转换工具之后会把 QUAL 写成缺失值。这时候如果你用bcftools view -i QUAL30表达式对缺失值求值为假结果就是所有位点全被删掉。我遇到过不止一次有人跑完发现输出是空文件查了半天以为是命令写错其实是 QUAL 全是点。所以过滤前一定先确认bcftools query -f %QUAL\n input.vcf.gz | head -20如果全是点那就别用 QUAL 过滤改用 DP 加基因型质量。或者用bcftools fill-tags之类的方式补但那通常没有意义因为 QUAL 本身就没算。DP深度分两个层面INFO/DP 是位点级别的总深度或者平均深度不同 caller 定义不一样FORMAT/DP 是样本级别的深度。位点级过滤用INFO/DP样本级过滤用FMT/DP。下限怎么定经验值是≥ 8~10。为什么不是 5因为二倍体基因型分型至少需要两条 read 支持等位基因 A、两条支持等位基因 B再加上测序错误和比对不确定低于 8 条的位点假阴性率明显上升。我做过的几组对比里DP 从 10 降到 5杂合位点会多出 15% 左右但这些多出来的位点里很大一部分在独立数据集上没法重复。上限怎么定前面已经说了中位深度的 2~3 倍。用datamash算出来的中位数乘以 3往上取整到整数。这个值我一般会稍微放宽一点因为过滤太紧会损失掉真实的 CNV 区域信号。比如中位数 14我可能定maxDP 50而不是 42。用 bcftools 一次性搞定位点级的 QUAL 和 DPbcftools view \ -i QUAL30 INFO/DP10 INFO/DP50 \ -Oz -o 02_site/site_filtered.vcf.gz \ 01_norm/normalized.vcf.gz注意如果你的 VCF 里没有 INFO/DP 这个键有些 caller 只在 FORMAT 里写 DP上面这条命令会因为表达式里的键不存在而报错或者更糟——静默地一个位点都不保留。跑之前用bcftools view -h检查一下 INFO 行里有没有IDDP。3.2 缺失率与等位基因频率缺失率用--max-missing控制取值是 0~1 的小数表示允许的最大缺失比例。注意它是缺失不是保留--max-missing 0.9意思是缺失率不能超过 10%也就是至少 90% 的样本要有基因型。这个方向经常有人搞反写成--max-missing 0.1结果是只保留缺失率低于 10% 的位点——虽然严格说这个写法也是合法的只是你的意图可能完全不是这样。我做群体分析常用0.9做小样本比如几十个样本时会提到0.95 甚至 1.0因为小样本里缺失率高的位点对统计量的扰动特别大。做全基因组关联时一般用0.95宁可位点少一点也要保证每个位点的样本量一致。MAF次等位基因频率是群体遗传里最有争议的参数之一。取值习惯从 0.05、0.01 到 0.001 都有完全取决于研究目的研究目的建议 MAF说明PCA / 群体结构0.05常见的多态位点足够刻画结构常规 GWAS0.01兼顾低频变异的效应罕见变异关联不设或 0.001需要保留低频位点系统发育树不设只关心是否有变异不关心频率亲缘关系推断0.05低频位点的基因型误差会干扰 IBD 估计vcftools 里就是--maf 0.05。bcftools 也可以直接做bcftools view -q 0.05:minor input.vcf.gz -Oz -o maf_filtered.vcf.gz-q是--min-af的简写:minor后缀表示按次等位基因频率计算另外还有:nref按参考等位基因和:alt1按第一个替代等位基因两种模式。这个参数在 bcftools 1.11 之后才有用之前确认版本。这里必须强调一个顺序问题MAF 和缺失率必须在样本集合确定之后计算。假如你有一千个样本来自三个不同地区你打算只分析其中一个地区的三百个样本那一定要先用bcftools view -S samples.txt或vcftools --keep把样本裁出来再做 MAF 过滤。否则你会删掉那些在全球样本里频率低、但在你的目标群体里频率正常的位点这种损失是不可逆的。3.3 基因型质量 GQ 与哈迪-温伯格平衡GQGenotype Quality是 caller 对单个样本基因型判断的置信度也是 Phred 标度。GQ20 表示这个基因型分错的概率大约 1%。我一般用GQ≥20作为阈值比较严格的场景用 30。这里要区分清楚操作方式bcftools 的-e FMT/GQ20会删掉整个位点只要有一个样本不满足就删——这通常不是我们想要的行为。真正想要的是把 GQ 低的那个样本的基因型置成./.保留位点。这个操作用 vcftools 更顺手vcftools --gzvcf in.vcf.gz \ --minGQ 20 \ --minDP 8 --maxDP 50 \ --recode --recode-INFO-all --stdout \ | bgzip -c 03_geno/gq_dp_filtered.vcf.gz--minGQ会把低于阈值的基因型置为缺失而不是删除位点。--minDP和--maxDP同理作用在 FORMAT/DP 上。这几个参数是 vcftools 相对 bcftools 最有价值的地方语义清晰、行为符合直觉。需要注意--minDP/--maxDP生效之后某些位点可能会因为大量基因型被置缺失而变得缺失率很高。所以基因型级过滤要放在缺失率过滤之前这样缺失率过滤才能看到真实的缺失情况。如果你反过来先做缺失率过滤、再做基因型过滤最后的结果里会残留一堆过滤后变成高缺失的位点。这个顺序问题我在早期项目里吃过亏写在这里提醒一下。哈迪-温伯格平衡HWE是判断位点是否靠谱的经典手段参数是--hwe值是一个 p 值阈值。常规做法是--hwe 1e-6严格的用1e-3。但 HWE 过滤有一个非常重要的前提它只适合在对照组或者非病例样本上做。为什么因为如果一个位点真的和疾病相关那么在病例群体里它就必然偏离 HWE你要是拿全部样本病例对照算 HWE 然后过滤等于在把真实信号当噪声删掉。这是群体遗传里一个经典的翻车点。所以我的做法是如果数据集有明确的分组HWE 过滤只在对照组做或者干脆不做如果没有分组、纯粹是自然群体样本可以用1e-6这个比较宽松的阈值做一次轻度清理。vcftools 的--hwe默认只在有基因型的样本上计算还有一个--hwe的过滤选项可以通过--keep指定子集样本具体做法是先把对照样本 ID 写成文件然后--keep controls.txt --hwe 1e-6。3.4 多等位位点、SNP/INDEL 类型与归一化顺序多等位位点同一位置有 3 个或更多等位基因在原始 VCF 里很常见尤其是串联重复区域。很多下游工具处理不了多等位或者处理方式不一致有的只取前两个等位基因所以通常在流程开头就拆开bcftools norm -m -any -f ref.fa -Oz -o normalized.vcf.gz raw.vcf.gz-m -any表示把多等位拆成多条双等位记录-f ref.fa同时做参考基因组左对齐left-align和 REF 校验。左对齐为什么重要因为 INDEL 在不同样本里的表示方式可能不同比如ATA和AAT描述的是同一个变异如果不统一同一个变异会被当成两个不同的位点影响后续合并和频率计算。归一化必须在所有过滤之前做这是硬规则。原因有三一是拆分多等位会改变坐标和等位基因编号如果先过滤后拆分你过滤时看到的 MAF 和拆分后的 MAF 不是一回事二是左对齐会改变 INDEL 的位置如果先按 BED 区域过滤再左对齐边界附近的 INDEL 可能会被移出目标区域三是拆分之后每个记录才是真正的双等位MAF、HWE 这些统计量才有明确定义。类型筛选用 bcftools 最省事# 只留 SNP bcftools view -v snps -Oz -o snps_only.vcf.gz normalized.vcf.gz # 只留 INDEL bcftools view -v indels -Oz -o indels_only.vcf.gz normalized.vcf.gz # 排除 INDEL bcftools view -V indels -Oz -o no_indels.vcf.gz normalized.vcf.gz注意-vinclude variant type和-Vexclude variant type是方向相反的很容易搞混。-v snps是只要 SNP-V snps是不要 SNP。我在脚本里一般会加注释写清楚不然过两个月自己都得愣一下。另外提一句如果只做 SNP 分析INDEL 直接排除掉就行因为 INDEL 的分型准确性普遍低于 SNP在群体分析里引入的噪声比贡献的信息多。但如果是做结构变异或者功能注释INDEL 必须保留那就得单独设一套更严的参数。4. 把命令串成一条可复现的流水线单个命令学会了接下来是把它们串起来。这一步的关键不是命令本身而是目录结构、中间文件策略和可复现性。我见过太多项目半年后想重跑发现当时的命令散落在几个终端历史里参数也不知道是哪个版本。4.1 目录结构与中间文件策略我的目录习惯是这样的project/ ├── raw/ # 原始 VCF只读绝不修改 ├── 01_norm/ # 归一化、拆分多等位后的文件 ├── 02_samples/ # 样本裁剪后的文件 ├── 03_site/ # 位点级过滤后 ├── 04_geno/ # 基因型级过滤后 ├── 05_final/ # 最终分析就绪文件 ├── logs/ # stats、日志 └── scripts/ # 脚本为什么按数字编号因为这一步的顺序是有严格依赖的编号能强制你和别人按顺序理解。为什么raw/设为只读因为一旦你误改了原始文件整个流程就再也无法复现了。我一般直接chmod -R a-w raw/。中间文件格式的选择也很关键。我的经验是中间步骤一律用压缩 BCF-Ob或者管道传输的未压缩 BCF-Ou不要用文本 VCF。原因很直接——BCF 是二进制格式体积小、读写快而且不会在文本解析上浪费时间。同样一份 100 万位点的数据VCF.gz 和 BCF 的处理时间能差出 2 到 3 倍。管道式串联的典型写法bcftools view -Ou 01_norm/normalized.vcf.gz \ | bcftools view -i QUAL30 INFO/DP10 -Ou \ | bcftools norm -d all -Ou \ | bcftools view -Oz -o 03_site/site.vcf.gz这里的bcftools norm -d all是用来去除完全重复的位点记录。有些 VCF 里会出现同一坐标、同一等位基因的多条记录多来自合并流程不去重的话后面统计位点数会虚高。-Ou传给下一个命令全程不落盘速度快很多。最后一步才用-Oz输出压缩 VCF方便给其他工具读。4.2 完整脚本从原始 VCF 到分析就绪文件下面是我实际用的一套脚本参数都带了注释直接改改就能用#!/usr/bin/env bash set -euo pipefail RAW${1:?usage: pipeline.sh raw.vcf.gz ref.fa prefix} REF$2 PREFIX${3:-clean} mkdir -p 01_norm 02_samples 03_site 04_geno 05_final logs # ---------- Step 1: 归一化 多等位拆分 去重 ---------- bcftools norm -f $REF -m -any -Ou $RAW \ | bcftools norm -d all -Oz -o 01_norm/${PREFIX}.norm.vcf.gz bcftools index -t 01_norm/${PREFIX}.norm.vcf.gz # ---------- Step 2: 样本裁剪可选如果没有 samples.txt 就跳过---------- if [ -s samples.txt ]; then bcftools view -S samples.txt -Oz -o 02_samples/${PREFIX}.samples.vcf.gz \ 01_norm/${PREFIX}.norm.vcf.gz bcftools index -t 02_samples/${PREFIX}.samples.vcf.gz WORK02_samples/${PREFIX}.samples.vcf.gz else WORK01_norm/${PREFIX}.norm.vcf.gz fi # ---------- Step 3: 只留 SNP ---------- bcftools view -v snps -Oz -o 03_site/${PREFIX}.snp.vcf.gz $WORK bcftools index -t 03_site/${PREFIX}.snp.vcf.gz # ---------- Step 4: 位点级硬过滤 ---------- bcftools view \ -i QUAL30 INFO/DP10 INFO/DP50 \ -m2 -M2 \ -Oz -o 03_site/${PREFIX}.site.vcf.gz \ 03_site/${PREFIX}.snp.vcf.gz bcftools index -t 03_site/${PREFIX}.site.vcf.gz # ---------- Step 5: 基因型级过滤 群体统计过滤 ---------- vcftools --gzvcf 03_site/${PREFIX}.site.vcf.gz \ --minDP 8 --maxDP 50 --minGQ 20 \ --max-missing 0.9 --maf 0.05 \ --recode --recode-INFO-all --stdout \ 04_geno/${PREFIX}.geno.vcf bgzip -f 04_geno/${PREFIX}.geno.vcf tabix -p vcf 04_geno/${PREFIX}.geno.vcf.gz # ---------- Step 6: 输出最终文件与统计 ---------- cp 04_geno/${PREFIX}.geno.vcf.gz 05_final/${PREFIX}.final.vcf.gz tabix -p vcf 05_final/${PREFIX}.final.vcf.gz bcftools stats 05_final/${PREFIX}.final.vcf.gz logs/${PREFIX}.final.stats.txt echo done. sites: bcftools view -H 05_final/${PREFIX}.final.vcf.gz | wc -l几个细节值得展开说。set -euo pipefail是必须的。-e让任何命令失败就退出-u让未定义变量报错-o pipefail让管道中任何一环失败都能被捕获。没有pipefail的时候管道中间的命令失败了你根本不知道最后拿到一个残缺的结果文件还以为一切正常。-m2 -M2限制至少 2 个、最多 2 个等位基因。经过 Step 1 的拆分之后所有记录本来都是双等位这个参数更多是防御性的——万一有记录没拆干净这里能挡住。vcftools 的--stdout配合 file.vcf是输出未压缩 VCF然后再bgzip。有人会问为什么不直接--recode --out prefix那样 vcftools 会自己生成prefix.recode.vcf。原因是--out只能输出到固定文件名不好接管道而且生成的是未压缩文件大文件时会额外占磁盘。--stdout更灵活。但要注意vcftools 输出的是普通 gzip 还是 bgzip 得看清楚——--stdout出来的是纯文本你用bgzip压就没问题如果直接gzip后面tabix会报错因为 tabix 只认 bgzip 的块压缩格式。这个坑我踩过报错信息还挺隐晦。--recode-INFO-all一定要加。vcftools 默认在 recode 的时候不输出 INFO 字段老版本行为不加这个参数你的 INFO 信息全没了后面想按 INFO 过滤或者做注释就抓瞎。4.3 用 fill-tags 把 vcftools 的活搬到 bcftools前面说过 vcftools 在大数据集上吃内存。如果你的数据规模到了千万位点级别有个替代方案是把群体统计量提前算好写进 INFO然后用 bcftools 的表达式过滤。fill-tags插件就是干这个的bcftools fill-tags 03_site/site.vcf.gz -Oz -o 03_site/tagged.vcf.gz \ -- -t AN,AC,AF,MAF,NS,ExcHet,HWE,F_MISSING bcftools index -t 03_site/tagged.vcf.gz这一步会往 INFO 里写入AN有效等位基因数、AC等位基因计数、AF等位基因频率、MAF次等位基因频率、NS有基因型的样本数、ExcHet杂合子过剩检验、HWE哈迪-温伯格检验 p 值、F_MISSING缺失率。有了这些字段一条表达式就能把 vcftools 那一堆参数全做完bcftools view \ -i MAF0.05 F_MISSING0.1 HWE1e-6 AN18 \ -Oz -o 04_geno/expr_filtered.vcf.gz \ 03_site/tagged.vcf.gz这条命令的意思和前面 vcftools 那套参数基本等价MAF≥0.05、缺失率≤10%、HWE p 值大于 1e-6、有效等位基因数至少 18对应 9 个样本可以理解为一个硬性的样本量下限。速度上fill-tags加表达式过滤比 vcftools 快大概 3 到 5 倍内存占用只有十分之一左右。代价是 HWE 的计算方式可能和 vcftools 略有差异一个是精确检验一个是对数似然近似在极端阈值下结果会有些出入。注意fill-tags是插件--后面的-t是插件自己的参数不是 bcftools 主程序的参数。忘了写--会报参数错误。另外不同版本支持的可填字段不同跑之前用bcftools fill-tags -- -l看一下支持列表。4.4 运行耗时与内存实测参考给几组我自己跑出来的数据方便你估资源。数据集规模bcftools 位点过滤vcftools 群体过滤fill-tags 方案峰值内存1 万位点 × 100 样本 5 秒约 20 秒 5 秒 200 MB100 万位点 × 500 样本约 40 秒约 6 分钟约 90 秒3~6 GB500 万位点 × 1000 样本约 4 分钟约 45 分钟约 8 分钟25~40 GB2000 万位点 × 2000 样本约 18 分钟建议放弃约 35 分钟80 GB几组数据横向看下来结论很清楚位点越密、样本越多vcftools 的劣势越明显。一百万位点、五百样本这个量级vcftools 还能接受但一上到五百万位点时间差距就拉到了十倍以上。所以我的建议是数据规模超过三百万位点就切到 fill-tags 方案。另外一个小技巧加线程。bcftools view --threads 4对压缩输出有加速效果但瓶颈通常在解压和表达式求值上加线程的收益不如想象中大。真正提速的手段还是减少中间 IO也就是用-Ou管道。5. 报错排查与踩坑实录前面讲的是应该怎么做这一节讲做错了会怎样。我把这些年遇到过的典型问题整理成一张速查表再补充几个不容易发现的隐性错误。5.1 常见报错速查表报错/现象根本原因解决办法reference allele mismatch at chr1:12345REF 与参考基因组不一致或染色体命名不匹配核对命名风格用bcftools annotate --rename-chrs统一后重跑Could not read from ...vcf.gz文件不是 bgzip 压缩的或索引缺失用bgzip重新压缩再tabix -p vcf建索引tabix: could not find index file.tbi或.csi不存在tabix -p vcf file.vcf.gz或bcftools index -t file.vcf.gzError: the tag INFO/DP is not definedVCF 里没有 INFO/DP 字段先bcftools view -h检查改用FMT/DP或去掉该条件过滤后位点数为 0表达式里的字段全是缺失值常见于 QUAL 为.先bcftools query打样本看取值vcftools 报Maximum number of alleles多等位位点未拆分先跑bcftools norm -m -anyvcftools 内存溢出被 kill数据集过大vcftools 全量载入先裁位点或改用 fill-tags 方案--max-missing 90报参数错误该参数需要 0~1 的小数改成0.9输出 VCF 丢失 INFO 字段没加--recode-INFO-all加上该参数重跑bcftools norm输出行数比输入多多等位被正常拆分这是预期行为无需处理但要知道位点计数定义变了5.2 过滤完数据看起来不对的几种表现有些问题不会报错只能靠观察。下面几种是我遇到过最多的。第一样本数莫名其妙变少了。用bcftools view -H file | head -1 | cut -f10- | wc -w数一下样本列数和原始文件对比。样本数变化通常是因为某个中间步骤用了--keep或者-S但样本 ID 文件里有多余的空行或制表符。我遇到过一次是因为样本 ID 文件是 Windows 换行符\r被当成了 ID 的一部分结果只匹配上了几个样本。第二杂合率异常偏低。跑完过滤后重新算一下每个样本的杂合比例bcftools stats -s - filtered.vcf.gz | grep ^PSC | awk {print $3, $5/($4$5)} | head如果大量样本的杂合比例掉到 0.1 以下说明--minDP设得太高把真实的杂合位点都置成缺失了。杂合位点理论上需要两条 read 分别支持两个等位基因深度要求本来就比纯合位点高--minDP一旦超过 15 就会明显伤到杂合位点。这是个隐性的偏倚很容易被忽略。第三MAF 分布形态不对。一个健康的自然群体数据集MAF 分布应该是低频位点多、高频位点多、中间少的 U 型或者 L 型。如果你过滤后发现 MAF 集中在 0.2~0.5 之间说明过滤太狠把大量真实的多态位点删掉了。命令行看一眼bcftools fill-tags filtered.vcf.gz -Ou -- -t MAF \ | bcftools query -f %INFO/MAF\n \ | awk {bint($1*10); c[b]} END {for (i in c) print i/10, c[i]} | sort -n这个分布画出来心里就有数了。第四位点坐标和原来对不上。如果流程里做了左对齐INDEL 的坐标会移动这是正常的。但如果你发现 SNP 的坐标也变了那说明中间某一步用了错误的参考基因组。SNP 的左对齐不会改变坐标因为 SNP 长度是 1坐标变了就说明有问题。5.3 几条我自己反复用得上的经验经验一永远先跑一版保守参数看结果再收紧。我的习惯是先跑QUAL20、max-missing 0.8、MAF 0.01这一版最宽松的看看剩多少位点。如果剩几十万那就可以收紧到QUAL30、0.9、0.05如果只剩几千那说明原始数据质量本身就有问题得回头查上游流程而不是靠过滤硬救。经验二每一步都留日志记录过滤前后的位点数。vcftools 会自动生成.log里面清清楚楚写着After filtering, kept X out of a possible Y Sites。bcftools 没有这个需要自己加before$(bcftools view -H in.vcf.gz | wc -l) # ... 过滤 ... after$(bcftools view -H out.vcf.gz | wc -l) echo sites: $before - $after logs/filter_counts.txt这份计数记录在写论文方法部分的时候就是现成的材料不用回头翻终端历史。经验三不要在流程中间临时改参数。我见过有人跑到一半觉得 MAF 太严手动改了一下继续跑结果不同批次的数据用了不同参数最后合并的时候完全没法解释。参数要么写死在脚本里要么做成配置文件的变量所有批次统一读同一份配置。经验四最终文件一定要建索引。这个听起来是废话但我确实见过交付出去的数据集没有.tbi文件下游同事用 IGV 打不开来回沟通了半天。养成习惯流程最后一步统一跑tabix -p vcf。经验五不要迷信任何一套标准参数。网上流传的QUAL30, DP10, MAF0.05, max-missing 0.9这套组合对很多数据集确实好用但它是从人类重测序数据总结出来的。如果你是植物、微生物或者非模式生物深度分布、杂合率、连锁不平衡衰减速度都不一样参数必须重新调。我自己做过的几套植物数据里maxDP定到中位深度的 5 倍才合理因为多倍体和重复序列多深度分布的尾部特别长。最后再分享一个排查思路当你怀疑某个过滤步骤有问题时不要一次改一堆参数而是把数据集切小——取一条染色体、前 10000 个位点跑一遍完整流程用bcftools query把每个位点在过滤前后的 QUAL、DP、MAF 都 dump 出来然后在表格里对比。到底是哪一步把哪些位点删掉了一目了然。这比对着几百万位点猜要高效得多我用这招定位过好几次位点莫名消失的问题其中一次就是因为 VCF 里 INFO/DP 的定义是所有样本 DP 之和而不是平均深度导致maxDP阈值把几乎全部位点都挡掉了——这种问题只看命令是看不出来的必须看数据本身。