bcftools+vcftools VCF过滤：从原始变异到干净SNP位点集-北京尧图网络

手上拿到一份 HaplotypeCaller 或者 DeepVariant 吐出来的原始 VCF打开一看几十万到上千万个位点密密麻麻铺满整条染色体直接丢给下游做 GWAS、群体结构分析或者亲缘关系推断结果大概率是先跑出一堆假阳性再花两三天回头排查是哪个位点把 PCA 图拽歪了。所以基于 bcftools vcftools 的 vcf 文件快速过滤流程本质上是每一个做变异分析的人都绕不过去的第一道闸门——它决定你后面所有分析可信度的上限。这篇内容我想按实际干活的顺序把从原始 VCF 到干净位点集的完整链路拆开讲清楚bcftools 和 vcftools 这对组合各自擅长什么、硬过滤阈值怎么定出来、哪些参数是拍脑袋抄的、哪些参数必须自己算。适合已经跑过变异检测、手里有 VCF 但不确定该砍掉哪些位点的朋友也适合刚入行想一次性把过滤逻辑搞明白的初学者。1. 变异过滤到底在解决什么问题从原始 VCF 到可用位点集1.1 原始 VCF 里到底有哪些不能直接用的东西先说清楚 VCF 这个格式。它是一个纯文本的制表符分隔文件头部以 ## 开头记录元信息包括参考基因组版本、contig 列表、INFO 字段定义、FORMAT 字段定义真正的数据行每行一个变异位点前面是 CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO 八列固定字段后面从 FORMAT 开始每个样本占一列。看起来规整但能解析和能用完全是两码事。一份未经处理的原始 VCF 里通常混着这几类垃圾位点一是低质量位点QUAL 值只有个位数本质上是测序噪声被错误识别成了变异二是深度严重不足的位点比如全基因组 30x 测序某个位点只有 2x 覆盖却被判成杂合这几乎肯定是错的三是深度异常高的位点常见于重复序列区域或者比对上发生了多位置匹配的片段深度可能是平均值的十几倍四是大片段缺失或者复杂结构变异断点附近的位点比对本身就不靠谱五是多等位位点一个位置有三个以上的等位基因很多下游工具尤其是老的 PLINK 流程处理不了。还有一类更隐蔽的问题REF 碱基和参考基因组对不上。这种情况多半是参考基因组版本搞错了比如变异检测用的是 hg19 但过滤时加载了 hg38 的 fasta一旦不一致后面所有基于 REF 的比对全部错位但程序不一定报错只会默默给你一个看起来正常的 VCF。1.2 为什么过滤是绕不开的一步有人会问能不能不做过硬过滤直接交给下游工具的贝叶斯模型去处理噪声答案是交给下游处理的前提是下游有能力处理。GWAS 的关联检验对每个位点独立做回归它没有任何机制去识别这个位点整体深度只有 1x它只会把这个位点的基因型当成真实观测值代入模型。真实观测值里混进系统性误差关联结果的膨胀系数就会上去QQ 图尾巴翘得老高你再去用基因组控制Genomic Control硬压只是把信号一起压没了。从数据量角度看过滤也是降本增效的刚需。一份 1000 样本的全基因组 VCF未过滤前可能有 500 万到 1000 万个位点压缩后动辄上百 GB。经过合理的硬过滤通常能保留 300 万到 800 万个高质量 SNP文件体积能砍掉一半下游每次读取、每次 LD 计算的时间成倍减少。我在实际项目里做过对比同一个 800 样本的数据集过滤前跑一次全基因组 LD 剪枝要 4 个多小时过滤后压到 1 小时出头省下来的时间成本非常实在。所以过滤这件事做得对是提纯做得糙就是自残。阈值太松假阳性横行阈值太狠真信号被误杀尤其是低频变异和稀有变异本来频率就低再被 MAF 阈值一刀切掉后续做负荷检验burden test就没数据了。整套流程的价值就在于用一套可解释、可复现的参数把噪声和信号分开。1.3 bcftools 与 vcftools 的分工逻辑这两个工具经常被放在一起用但它们的设计哲学差别挺大。bcftools 属于 htslib 生态底层用 C 实现支持流式处理、管道串联、多核压缩输出速度非常快适合做粗筛和大批量位点操作比如规范化、拆多等位、排序、去重、基于 INFO 字段的快速过滤。它的表达式引擎filter expressions功能很强可以写QUAL30 INFO/DP10这类条件还内置了 F_MISSING、MAF、AC、AN 这些计算标签很多原来要靠脚本统计的量它直接给你算好了。vcftools 是老的 Perl/C 混合工具作者已经多年不再更新速度慢、单线程、内存占用高但它的优势在于基因型级别的过滤极其直观--minDP、--maxDP、--minQ这些参数直接作用在 FORMAT 字段上逐个样本逐个基因型地判断语义清晰写起来短。它还有一大堆输出统计的功能--freq出等位基因频率、--missing出缺失率、--hardy出哈迪温伯格平衡检验结果、--het出个体杂合度都是查问题时的常用工具。我个人的分工习惯是bcftools 打头阵做规范化排序去重位点级粗过滤把数据量迅速压下来vcftools 收尾做基因型级细过滤同时输出统计用于校验。这个顺序的好处是让慢的 vcftools 处理的已经是压缩过的数据整体耗时能省很多。反过来的话你让 vcftools 先去啃几十 GB 的原始文件等着看进度条慢慢爬吧。2. 工具选型与运行环境准备2.1 为什么是 bcftools vcftools 这个组合市面上的 VCF 处理工具不少GATK 的 VariantFiltration 也能做硬过滤PLINK 的--maf --geno --hwe也能过滤为什么偏偏选这两个我说说实际考量。GATK VariantFiltration 的过滤逻辑是打标签而不是删位点它把不满足条件的位点在 FILTER 列写上标记真正的删除要靠 SelectVariants 再走一步流程偏长而且 GATK 对 Java 环境、内存配置有要求跑大文件比较吃资源。PLINK 确实快但它的输入是 bed/bim/fam 二进制格式你需要先转换一次转换过程本身就可能丢信息比如多等位位点会被强制处理INDEL 处理也比较粗糙。bcftools vcftools 的组合优势在于都是直接吃 VCF不需要格式转换中间产物可以用管道直接对接每一步的输入输出都是标准 VCF出了问题容易定位在哪一步。而且这两个工具对 VCF 规范的支持都很完整INFO、FORMAT 字段能原样保留vcftools 记得加--recode-INFO-all不会像某些工具那样把注释信息吃掉。需要提醒的是vcftools 从 2018 年之后基本停止维护官方仓库已经归档。这意味着它对超大样本比如几万样本的队列数据支持不好容易爆内存。如果你的数据规模到了这个量级建议直接用 bcftools 的表达式把 vcftools 那部分功能替代掉或者换用更新的 cyvcf2、pysam 写脚本。中小规模数据几百到几千样本继续用 vcftools 完全没问题熟练度高、语义清楚效率反而更高。2.2 安装与版本核对安装最省事的方式是走 conda把依赖一起装好避免 htslib 版本冲突带来的那些莫名其妙的报错。conda create -n vcf-filter -c conda-forge -c bioconda \ bcftools1.19 vcftools0.1.16 htslib1.19 tabix1.19 -y conda activate vcf-filter装完第一件事是核对版本因为不同版本的参数行为会变。bcftools 1.10 之前和之后在-i表达式的标签支持上有差别vcftools 0.1.14 和 0.1.16 在--max-missing的边界处理上也略有不同。bcftools --version vcftools --version tabix --version提示生产流程里一定要把版本号写进脚本注释或者日志里。我踩过一次坑同一套过滤脚本在本地跑出来的位点数和服务器上差了十几万查了半天发现是服务器上装的是旧版 bcftoolsF_MISSING标签还没支持条件表达式被静默忽略了一部分。如果不用 conda从源码编译也可以但要注意先编译 htslib再编译 bcftools并且用--with-htslib指定路径否则容易出现运行时报缺少某个符号这种错误。2.3 输入文件的预处理为什么要先规范化在动过滤参数之前有一个几乎所有人都会忽略但影响极大的步骤规范化normalization。规范化包含三件事每一件都有实际必要性。第一件向左对齐left-align。INDEL 在 VCF 里的表示方式不唯一比如一个缺失可以有多种等价写法ATA和AAT描述的是不同位置的变异。如果不统一到最左侧的等价表示同一个变异在不同样本或不同批次里会被当成不同位点做合并或者比较时就会出错。bcftools norm 默认就做左对齐。第二件拆多等位。一个位点如果 ALT 列写的是A,G说明这里有两个非参考等位基因。拆成两个独立的记录后下游做频率统计和关联分析会简单很多。用-m-any参数拆分拆分后 INFO 里的 AC、AN 这些计数会被自动重算。第三件REF 校验。用--check-ref w参数让 bcftools 拿参考基因组去核对每个位点的 REF 是否匹配。不匹配的位点会被警告甚至修正。这一步的价值在于如果参考基因组版本搞错了你在这里就能发现而不是等到下游分析结果诡异再去回溯。bcftools norm -f genome.fa -m-any --check-ref w \ -Oz -o 01.norm.vcf.gz raw.vcf.gz bcftools index -t 01.norm.vcf.gz-Oz表示输出 bgzip 压缩的 VCF-t表示建立 tabix 索引。索引这个事必须记牢任何按区域查询、任何需要随机访问的操作都依赖索引没有索引时 bcftools 会直接报错。国内不少同行习惯性地输出未压缩的.vcf然后抱怨bcftools view -r chr1用不了其实就是缺索引。3. 核心过滤参数的原理与取值依据3.1 位点层面过滤QUAL、DP、缺失率怎么定位点层面的过滤是最粗的一刀但也是最有效的。三个核心指标QUAL、深度、缺失率。QUAL 是 Phred 化的质量值定义是-10 * log10(P)P 是这个位点被错误识别的概率。所以 QUAL30 意味着错误概率约 1/1000QUAL20 意味着 1/100QUAL50 意味着 1/100000。常规全基因组测序数据我一般设 QUAL≥30。如果是目标区域捕获或者低深度数据可以适当放宽到 20但要接受更高的假阳性率。别小看这个换算很多人的过滤器里 QUAL 阈值写 20 却以为是 20% 的准确率实际上已经是 99% 了。深度 DP 需要分两层理解。INFO 里的 DP 是位点所有样本的覆盖总和FORMAT 里的 DP 是单个样本在该位点的覆盖深度。bcftools 的表达式里INFO/DP和FORMAT/DP是两个不同的东西这个区别后面讲 vcftools 的时候还要再强调一次。深度阈值怎么定假设你的数据平均测序深度是 30x二倍体。理论上每个等位基因平均有 15x 覆盖。用泊松分布算一下λ15 时观测到覆盖数小于等于 5 的概率大概是 0.0008 左右非常低。也就是说一个真实的杂合位点深度掉到 5 以下的概率极小。所以 minDP10 是一个比较稳妥的下限它既过滤掉了绝大多数假阳性又不会误杀真实变异。当然如果你的数据平均深度只有 10x比如低覆盖全基因组那 minDP 得降到 3 到 5同时用缺失率来兜底。深度上限同样重要。一个位点深度达到平均深度的 2.5 到 3 倍以上就要怀疑是不是比对到了重复区域或者存在拷贝数变异这类位点的基因型判读非常不可靠。实际做法是用bcftools query -f %INFO/DP\n把深度分布导出来取 99 分位数作为上限参考而不是凭空写个 100。我在一个人类全基因组项目里量过深度 99 分位在 78 左右所以 maxDP 设 80 比较合理如果草率设成 200那些 150x 的重复区域位点就全留下了。缺失率指的是有多少比例的样本在这个位点没有有效基因型基因型是./.。F_MISSING是 bcftools 内置标签F_MISSING0.2表示保留至少 80% 样本有基因型的位点。缺失率高通常意味着这个区域比对困难、或者测序覆盖有系统性空洞保留下来会严重影响下游分析的样本量。3.2 基因型层面过滤GQ、个体深度、杂合度位点级过滤是整列砍基因型级过滤是单元格砍。后者更精细但也更慢所以放在流程后半段。GQ 是基因型质量值同样是 Phred 标度。GQ20 意味着这个样本的基因型判读有 99% 的可信度。注意 GQ 是 FORMAT 字段每个样本一个值如果某个样本的 GQ 低于阈值处理方式有两种一是把该基因型置为缺失保留位点二是直接丢弃整个位点。vcftools 的--minQ采用的是前一种思路把低质量基因型改成缺失位点本身保留。这个处理更温和也更合理——毕竟一个位点只要还有足够多的样本基因型合格它依然有价值。个体深度--minDP/--maxDP作用在 FORMAT/DP 上逐样本判断。这里有个容易被忽略的细节当一个样本的基因型深度太低被置为缺失后位点的缺失率会上升可能触发后面的--max-missing再次过滤。所以参数的先后顺序是有讲究的vcftools 内部的处理顺序会影响最终结果建议先用小样本测试确认两个阈值叠加后的效果符合预期。bt 杂合度过滤是个进阶操作。有些样本存在样本污染或者近交导致的全基因组杂合度异常可以用vcftools --het计算每个个体的杂合度然后剔除偏离均值 3 个标准差以外的个体。这一步不属于标准流程但在做群体遗传分析时很有必要因为一两个异常样本足以把群体结构分析的主成分轴搅乱。3.3 等位基因频率与哈迪温伯格平衡怎么用MAF次等位基因频率阈值是最有争议的一个参数。设 0.05 是 GWAS 的常见做法理由是低频变异的统计功效不足强行纳入只会增加多重检验负担。但这个阈值对稀有变异研究是致命的做负荷检验、做罕见病关联MAF 阈值得放到 0.01 甚至更低甚至根本不用 MAF 过滤改用其他方式控制质量。这里有个坑必须说清楚bcftools 的 MAF 和 vcftools 的 MAF 计算基础不完全一样。vcftools 的--maf默认基于实际观测到的基因型计数计算缺失基因型不参与可以用--maf配合--max-missing控制bcftools 的MAF标签基于等位基因数 AN 计算而 AN 本身已经排除了缺失基因型。两者在缺失率高的位点上会给出不一致的结果。如果你在 bcftools 里已经用 MAF 过滤过一遍再在 vcftools 里过滤一次理论上不会误杀太多但如果两个阈值设得很接近边界位点会出现此消彼长的现象。我的建议是 MAF 只在一个工具里过滤不要两边都设。哈迪温伯格平衡HWE检验的用法更需要谨慎。--hwe 1e-6表示剔除 HWE 检验 p 值小于 1e-6 的位点。这个阈值的逻辑是严重偏离 HWE 的位点通常来自测序错误或者比对错误而不是真实的生物学现象。在 GWAS 的质控中这是标准操作因为病例和对照的比例偏差会造成 HWE 偏离而测序错误造成的偏离比真实生物学偏离更常见。但如果你是做群体遗传学、做自然群体重测序HWE 过滤要极其小心甚至完全不用。原因很简单自然群体普遍存在群体亚结构、近交、选择作用这些都会造成真实的 HWE 偏离。你用一个防错的过滤器去砍真实信号等于主动放弃了研究价值。我在一个野生群体项目里就犯过这个错用了 1e-6 的 HWE 阈值结果过滤掉的位点里有一批恰好落在受选择区域后来只能重跑。所以这个参数的默认建议是GWAS 用自然群体慎用且阈值不要低于 1e-6。3.4 阈值推导一个可以照着算的实例光讲原理不够我给一个完整的阈值推导过程。假设项目背景人类全基因组测序样本量 500测序策略 PE150平均深度 32x使用 GATK 做变异检测。目标构建用于 GWAS 的高质量 SNP 集。第一步算深度上下限。用 bcftools 统计全样本所有位点的深度分布bcftools query -f %INFO/DP\n raw.norm.vcf.gz \ | sort -n | awk {a[NR]$1} END {print mediana[int(NR*0.5)], p1a[int(NR*0.01)], p99a[int(NR*0.99)]}假设输出 median620500 样本 × 32x ≈ 16000 的理论总和但实际这里单位是位点总深度按需换算p1 和 p99 用来定上下限。实际更常用的做法是按平均深度的倍数来定下限取平均深度的 1/3约 10x上限取 2.5 倍约 80x。第二步算 QUAL 阈值。看 QUAL 的分布bcftools query -f %QUAL\n raw.norm.vcf.gz \ | awk {if($130) a; else b} END {print qual30:, a, qual30:, b, ratio:, b/(ab)}如果 QUAL30 的比例超过 20%说明数据质量本身有问题需要回头查比对和变异检测环节而不是靠过滤硬压。第三步定缺失率和 MAF。GWAS 常规F_MISSING0.05保留 95% 样本有位点的位点MAF0.05。样本量 500 的情况下MAF0.05 意味着次等位基因至少要出现约 50 次500 样本 × 2 拷贝 × 0.05统计功效基本够用。第四步定 HWE。GWAS 用p1e-6。把这些数字串起来就得到一套完整可执行的阈值。关键点在于每一个数字背后都有推导依据不是抄来的。抄来的阈值最怕的就是数据特征不匹配——你用别人 30x 数据的阈值去处理自己 10x 的数据过滤结果必然是灾难。4. 完整实操流程从原始 VCF 到干净位点集4.1 第一步和第二步规范化、排序、去重先把完整脚本的前半段给出来每一步的作用都有注释。#!/usr/bin/env bash set -euo pipefail REFref/GRCh38.fa RAWraw/cohort.raw.vcf.gz OUTwork THREADS8 mkdir -p $OUT # 步骤1标定参考基因组拆多等位向左对齐REF校验 bcftools norm -f $REF -m-any --check-ref w \ -Oz -o $OUT/01.norm.vcf.gz $RAW bcftools index -t --threads $THREADS $OUT/01.norm.vcf.gz # 步骤2排序规范化可能打乱原有顺序 bcftools sort -Oz -o $OUT/02.sorted.vcf.gz $OUT/01.norm.vcf.gz bcftools index -t --threads $THREADS $OUT/02.sorted.vcf.gz # 步骤3去重完全相同的变异记录只留一条 bcftools norm -d exact -Oz -o $OUT/03.dedup.vcf.gz $OUT/02.sorted.vcf.gz bcftools index -t --threads $THREADS $OUT/03.dedup.vcf.gz关于排序这一步有个细节值得说。bcftools norm 在拆多等位和左对齐之后位点位置会发生变化——左对齐可能把变异往左移动若干碱基拆多等位会插入新的记录行。这两件事都会打破原来的位置顺序所以排序不是可选项。如果不排序后面建索引会失败或者索引建出来是错的按区域查询会返回不完整的结果。bcftools index在遇到位置乱序的文件时一般会报错但某些版本会静默跳过非常危险。关于去重-d exact表示只删除完全相同的记录CHROM、POS、REF、ALT 全部相同。还有-d snps、-d indels等模式分别针对不同变异类型。全基因组数据里重复记录主要来自合并多个批次的 VCF或者变异检测工具的某些 bug。去重这一步不做下游做频率统计时会发现某个位点的等位基因数异常高。4.2 第三步bcftools 位点级快速过滤这一步是整个流程里性价比最高的几行命令就能把数据量砍掉一大半。# 只保留双等位 SNP应用位点级硬过滤 bcftools view \ -m2 -M2 \ -v snps \ -i QUAL30 INFO/DP600 INFO/DP4800 F_MISSING0.10 \ -Oz -o $OUT/04.snp.hard.vcf.gz $OUT/03.dedup.vcf.gz bcftools index -t --threads $THREADS $OUT/04.snp.hard.vcf.gz参数逐条解释-m2 -M2最小等位基因数 2最大等位基因数 2也就是只保留双等位位点。多等位位点在 PLINK、GEMMA 等下游工具里处理方式不一致直接排除最省心。如果确实需要研究多等位单独跑一条分支流程。-v snps只保留 SNP排除 INDEL 和其他类型。INDEL 的基因型判读误差更大常规 GWAS 一般不纳入。-i ...include 表达式满足条件的保留。注意这里 INFO/DP 的阈值要根据样本数换算。500 样本 × 平均每样本 12x考虑覆盖不均≈ 6000那下限设 600相当于每样本平均 1.2x 的下限太松了实际应该更严格。我这里举例的 600 到 4800 是按 500 样本、平均每样本 6x 到 48x 的区间来理解的你在实际使用时必须按自己数据的真实分布重新计算千万别直接抄。F_MISSING0.10缺失率低于 10%。这一条在粗过滤阶段设置得比最终阈值松一些因为后面 vcftools 的基因型过滤会进一步增加缺失率留点余量。注意-i表达式里的字段名大小写敏感INFO/DP和FORMAT/DP完全不同。还有一点bcftools view的-i和-e是互斥的前者是保留满足条件的后者是排除满足条件的写反了会把好位点全删掉。这种错误在数据量大的时候不容易发现因为过滤完还有位点只是数量不对。养成习惯过滤前后都跑一次bcftools view -H file.vcf.gz | wc -l记录位点数。4.3 第四步vcftools 基因型级精细过滤粗过滤完成后数据量已经降下来了这时候用 vcftools 做精细处理。# 基因型级过滤同时保留 INFO 字段 vcftools --gzvcf $OUT/04.snp.hard.vcf.gz \ --minDP 8 \ --maxDP 80 \ --minQ 20 \ --max-missing 0.90 \ --maf 0.05 \ --hwe 1e-6 \ --recode --recode-INFO-all \ --stdout $OUT/05.final.vcf # 压缩并建索引 bgzip - $THREADS -c $OUT/05.final.vcf $OUT/05.final.vcf.gz tabix -p vcf $OUT/05.final.vcf.gz几个关键细节第一--minDP和--maxDP这里作用的是 FORMAT/DP逐样本判断低于下限或高于上限的基因型会被置为缺失不是删除位点。这和 bcftools 那一步的 INFO/DP 完全是两回事两步配合起来既控制了位点整体深度也控制了单个样本的基因型质量。第二--minQ 20对应 FORMAT/GQ。同样是把不合格的基因型置为缺失。GQ 这个字段不是所有变异检测工具都会输出DeepVariant 默认就不写 GQ这种情况下要跳过这个参数或者改用其他质量指标。第三--max-missing 0.90放在基因型过滤之后正好把因为低深度、低 GQ 被置为缺失后缺失率超标的位点再筛一遍。这个顺序是我试过比较合理的反过来先做缺失率过滤再置缺失会漏掉一部分过滤后缺失率超标的位点。第四--recode-INFO-all这个参数必须加。vcftools 的--recode默认只保留少量核心 INFO 字段AC、AN、AF 这些常用注释会被丢掉下游分析会缺数据。我见过不止一个人因为这个参数没加跑到下游发现 INFO 全是空的重头来一遍。第五--stdout配合重定向避免 vcftools 生成一堆中间文件。vcftools 默认在--recode时会写一个.recode.vcf文件多批次跑的时候容易覆盖或者混乱用--stdout更干净。4.4 第五步结果校验与统计输出过滤完不是就结束了必须做校验确认过滤行为符合预期。我通常跑这几项统计# 位点数和样本数 bcftools stats $OUT/05.final.vcf.gz $OUT/final.stats.txt # 每位点的等位基因频率分布 vcftools --gzvcf $OUT/05.final.vcf.gz --freq --out $OUT/final.freq # 缺失率分布 vcftools --gzvcf $OUT/05.final.vcf.gz --missing --out $OUT/final.missing # 哈迪温伯格平衡 vcftools --gzvcf $OUT/05.final.vcf.gz --hardy --out $OUT/final.hardybcftools stats的输出里有几个关键指标要重点看number of SNPs最终保留位点数、number of samples、ts/tv转换/颠换比。ts/tv 是衡量 SNP 质量最直观的指标人类全基因组数据经过合理过滤后ts/tv 应该落在 2.0 到 2.1 之间。如果低于 1.8说明假阳性残留还比较多需要收紧过滤条件如果高得离谱比如 2.5 以上可能是过滤过度或者数据本身有问题。这个指标在 vcftools 的--TsTv输出里也能看到而且能按不同 MAF 区间分开统计做精细评估更有用。再检查一下缺失率分布awk {if(NR1) print $5} $OUT/final.missing | sort -n | tail -1--missing输出里第 5 列是每位点的缺失比例看看最大值是不是还在可接受范围应该刚好卡在 0.10 附近因为--max-missing 0.90已经约束过了。如果最大值明显超过 0.10说明过滤参数没生效去检查参数拼写。最后一步把关键数字记录到一个日志文件里方便日后复现和对比。echo date$(date -Iseconds) $OUT/filter.log echo input_sites$(bcftools view -H $RAW | wc -l) $OUT/filter.log echo output_sites$(bcftools view -H $OUT/05.final.vcf.gz | wc -l) $OUT/filter.log4.5 一条完整的可复现脚本把上面的步骤串成一条完整脚本可以直接改路径复用。#!/usr/bin/env bash set -euo pipefail REFref/GRCh38.fa RAWraw/cohort.raw.vcf.gz OUTwork THREADS8 mkdir -p $OUT echo [1/6] normalize bcftools norm -f $REF -m-any --check-ref w \ -Oz -o $OUT/01.norm.vcf.gz $RAW bcftools index -t --threads $THREADS $OUT/01.norm.vcf.gz echo [2/6] sort bcftools sort -Oz -o $OUT/02.sorted.vcf.gz $OUT/01.norm.vcf.gz bcftools index -t --threads $THREADS $OUT/02.sorted.vcf.gz echo [3/6] dedup bcftools norm -d exact -Oz -o $OUT/03.dedup.vcf.gz $OUT/02.sorted.vcf.gz bcftools index -t --threads $THREADS $OUT/03.dedup.vcf.gz echo [4/6] hard filter (site-level) bcftools view -m2 -M2 -v snps \ -i QUAL30 INFO/DP600 INFO/DP4800 F_MISSING0.10 \ -Oz -o $OUT/04.hard.vcf.gz $OUT/03.dedup.vcf.gz bcftools index -t --threads $THREADS $OUT/04.hard.vcf.gz echo [5/6] genotype-level filter (vcftools) vcftools --gzvcf $OUT/04.hard.vcf.gz \ --minDP 8 --maxDP 80 --minQ 20 \ --max-missing 0.90 --maf 0.05 --hwe 1e-6 \ --recode --recode-INFO-all --stdout $OUT/05.final.vcf echo [6/6] compress index stats bgzip - $THREADS -c $OUT/05.final.vcf $OUT/05.final.vcf.gz tabix -p vcf $OUT/05.final.vcf.gz bcftools stats $OUT/05.final.vcf.gz $OUT/final.stats.txt echo done. sites: $(bcftools view -H $OUT/05.final.vcf.gz | wc -l)这套脚本我在好几个项目里跑过500 到 2000 样本规模的数据从原始 VCF 到最终干净位点集耗时大概在 1 到 3 小时主要瓶颈在 vcftools 那一步。如果数据量更大建议把 vcftools 那段用 bcftools 表达式替代。5. 常见问题与排查技巧实录5.1 常见报错与快速排查对照表实际跑流程时遇到的大部分问题都集中在几类我整理成表方便对号入座。报错或现象常见原因排查与解决index file is older than vcf fileVCF 被修改过但索引没重建重新执行tabix -p vcf或bcftools index -t -fcontig chr1 is not defined in the header染色体命名不统一或 header 缺 contig 行用bcftools annotate --rename-chrs统一命名或用bcftools reheader补 contig过滤后位点数几乎没变-i表达式里的字段名写错条件被忽略用bcftools view -H抽样检查确认用INFO/DP还是FORMAT/DP过滤后位点数几乎归零阈值过严或深度单位搞错先不加条件跑一次统计分布再定阈值vcftools 报VCF file has no genotype data输入是 sites-only VCF检查原始文件是否包含 FORMAT 列和样本列输出 VCF 的 INFO 字段大量缺失忘了加--recode-INFO-all重新跑参数补上内存爆掉OOM killed样本量太大vcftools 单线程加载全量数据按染色体拆分并行处理或改用 bcftoolsts/tv 低于 1.8假阳性残留多提高 QUAL 阈值收紧 DP 上限加大--minQREF allele mismatch警告大量出现参考基因组版本不匹配核对变异检测使用的参考版本重新规范化多等位拆分后位点数翻倍正常现象拆分后每条记录是独立位点位点数增加是预期行为5.2 怎么判断过滤过度还是过滤不足这是最考验经验的一步因为没有绝对标准只能靠几个交叉验证的信号。判断过滤不足的信号ts/tv 偏低QQ 图在尾部明显膨胀PCA 图上前几个主成分解释的方差比例异常高可能是错误位点主导个体杂合度分布有明显的离群点。这时候的处理方向是收紧 QUAL 和 DP 上限加严 HWE 阈值。判断过滤过度的信号ts/tv 异常偏高最终保留位点数远低于预期比如 500 样本全基因组只剩 50 万位点明显太少MAF 分布整体右移低频变异几乎被清空LD 衰减曲线异常陡峭说明位点密度不足无法有效捕捉连锁信息。这时候要放宽 MAF 阈值和 HWE 阈值。一个很实用的交叉验证方法把过滤前后的 VCF 各跑一次--freq对比等位基因频率分布的直方图。合理的过滤应该只是削掉分布两端的部分中间的峰形保持完整。如果峰形被明显改变说明过滤参数动了不该动的东西。还有一个技巧做一次过滤敏感性分析用三套阈值宽松、标准、严格分别过滤然后各自跑一遍下游的初步分析比如 PCA 或者 IBD 计算看结果的一致性。如果三套结果高度一致说明你的数据质量好阈值选择不敏感如果差异很大说明数据本身有质量问题需要回到上游排查而不是靠调阈值。5.3 性能与内存优化的一些实操心得第一批经验bcftools 的多核加速主要体现在压缩和解压上-Oz配合--threads能明显提速但过滤本身的逻辑计算是单线程的。真正吃时间的是 IO 和压缩所以把中间文件放在 SSD 上收益立竿见影。我在机械盘上跑同样的流程比在 NVMe 上慢了将近一倍。第二批经验按染色体拆分并行处理。人类的 22 条常染色体加性染色体可以拆成 24 个任务并行跑每个任务独立过滤最后用bcftools concat合并。这个方式对 vcftools 这种单线程工具尤其有效24 核跑满整体耗时能压到原来的十分之一。拆分命令for chr in $(seq 1 22) X Y MT; do bcftools view -r chr${chr} $OUT/04.hard.vcf.gz \ -Oz -o $OUT/split/chr${chr}.vcf.gz done wait合并的时候要注意顺序bcftools concat需要按染色体顺序提供文件列表否则合并出来的 VCF 位置乱序索引会失败。第三批经验vcftools 的内存占用和样本数的平方大致成比例。1000 样本时可能占 10 到 20 GB5000 样本时轻松超过 100 GB。如果你的机器内存有限要么按染色体拆要么按样本拆做成两个子集分别过滤再用bcftools merge合并要么干脆换用 bcftools 表达式。第四批经验善用中间文件的检查点。整个流程不要写成一个大脚本一气呵成中间加几个bcftools view -H | wc -l记录位点数。一旦最终结果不对你能立刻定位是哪一步出了问题比从头重跑快得多。我一般的习惯是每一步都输出位点数和文件大小到日志里。6. 进阶扩展让过滤流程更贴合实际项目6.1 不同研究目标对应不同的过滤策略同一份原始数据做 GWAS 和做群体遗传分析过滤策略应该完全不一样。GWAS 场景下质控的核心目标是去掉会造成假关联的位点。策略是严格的 MAF 阈值0.05 或更高、严格的缺失率小于 0.05、启用 HWE 过滤1e-6、严格控制个体杂合度离群样本。样本层面的质控同样重要剔除个体缺失率过高的样本、剔除亲缘关系过近的样本用--relatedness或者 KING 计算。群体遗传场景下目标变成尽可能保留真实变异同时控制错误率。策略是MAF 阈值放宽到 0.01 或不做、缺失率放宽到 0.2、完全不用 HWE 过滤或者只保留 p 值极低的结果作为排查线索而不删除、更多依赖深度和质量的硬指标。如果做的是低频变异或者结构变异相关的研究过滤策略还要再放宽。做基因分型一致性验证时思路又不一样这时候你要保留尽可能多的位点用重复样本之间的基因型一致性作为质量指标过滤重点放在重复样本间基因型冲突的位点上而不是常规的 QUAL、MAF。所以快速过滤流程这四个字里快速是手段真正的核心是匹配研究目标。一套脚本走天下的做法在真正严谨的项目里是行不通的。6.2 把流程固化成可复用的模块跑过几次之后你会发现很多参数是固定的变的只是输入文件和少数阈值。这时候值得把流程封装成带参数的 Shell 脚本或者用 Snakemake、Nextflow 写成工作流。用 Snakemake 的好处是自动处理依赖关系和并行调度每条规则只描述输入、输出和命令Snakemake 自己判断哪一步该跑、哪些可以并行。对多批次数据尤其方便加个expand就能生成所有批次的过滤任务。如果不想引入额外框架最简做法是写一个接受参数的 bash 脚本把阈值作为变量放在开头#!/usr/bin/env bash set -euo pipefail RAW$1 OUTDIR$2 REF${3:-ref/GRCh38.fa} MIN_QUAL${MIN_QUAL:-30} MIN_DP${MIN_DP:-8} MAX_DP${MAX_DP:-80} MIN_GQ${MIN_GQ:-20} MAX_MISSING${MAX_MISSING:-0.90} MAF${MAF:-0.05} HWE${HWE:-1e-6}这样每个项目只需要改环境变量不用改脚本主体。配合set -euo pipefail任何一步失败都会立即终止不会带着错误结果往下跑。关于版本控制我的建议是把过滤脚本和参数日志一起纳入 git 管理。生物信息分析的可复现性很大程度上依赖于三个月后你还能记起当时用了什么参数。日志文件里记录工具版本、参数值、输入文件的 md5这三点齐全基本可以做到完全复现。6.3 与下游工具的衔接过滤完的 VCF 是终点也是起点。往下游走最常见的几个方向转 PLINK 格式做 GWASplink --vcf final.vcf.gz --make-bed --out gwas。注意 PLINK 默认会丢弃 MAF 低于 0.01 的位点如果你的数据是稀有变异研究要加--maf 0.001显式放宽。转换后建议立刻跑一次plink --bfile gwas --freq和--missing复核确认转换过程没有引入意外变化。转 VCF 到二进制做快速查询bcftools view -Ob输出 BCF 格式二进制格式比压缩 VCF 更紧凑读取更快适合作为长期存储格式。BCF 同样需要索引用bcftools index --csi。做注释过滤后的位点集交给 VEP、ANNOVAR 或者 SnpEff 做功能注释。注释这一步最好放在过滤之后做因为注释很耗时对已经被过滤掉的位点做注释纯属浪费算力。做位点子集提取如果只关心某些基因或者某些区域用 bed 文件提取bcftools view -R target.bed或者vcftools --bed target.bed。这一步可以在过滤之前做也可以在之后做取决于你的目标是全基因组质控后提取子集还是对子集做精细质控前者更常见也更省事。我个人在实际操作中的体会是整个流程最容易出错的不是命令本身而是参数与数据特征的匹配。同一套参数A 项目跑得好好的B 项目直接翻车原因往往是两个项目的测序深度或者样本量差了一倍。所以每次开新项目我都会先花二十分钟做一轮分布统计看看 QUAL、DP、缺失率的实际分布长什么样再定阈值。这二十分钟的投入能省下后面反复调试的几天时间。另外一个小技巧把每次设定的阈值和最终的 ts/tv、位点数一起记在一个表格里跑的项目多了你会发现自己数据集的经验阈值区间下一次定参数会快很多。

bcftools+vcftools VCF过滤：从原始变异到干净SNP位点集

相关推荐

L型3麦音源定位：TDOA物理原理与工程最优解

游戏引擎入门：Cocos 引擎如何让你 5 分钟跑通第一个跨平台项目

Inngest Core API 深度解析：GraphQL 远程管理服务的架构与开发工作流

更多资讯

DeepSeek-Coder-6.7B本地部署全指南：硬件适配、GGUF格式与llama.cpp实战

NoETL明细语义层：让AI Agent真正读懂业务数据

Cesium自定义指南针：从坐标系原理到Canvas高性能实现

以太网温湿度传感器通信中CRC16与CRC32选型实战指南

COMSOL多物理场耦合在水力压裂模拟中的应用与优化

Spring Boot 实战：流浪宠物管理系统开发与部署全流程

今日更新

AWS无服务器应用开发指南：从Lambda到SAM的架构与实践

SpringBoot+Vue人事系统实战：解决HR数据断点与事务一致性

MATLAB手写紧束缚模型计算石墨烯能带结构

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？