第一次把 BAM 丢给 bcftools 想换成 VCF 的人几乎都会在头十分钟里被同一类报错拦住屏幕上刷出一行[E::fai_build3_core] Failed to open the file或者[E::hts_open_format] Failed to open file然后进程安静地退出只留下一个 0 字节的输出文件。更让人抓狂的是这些错误和你写的 mpileup 命令本身没什么关系问题往往出在那个你以为早就准备好了的参考基因组索引上。这篇内容就围绕bcftools这个工具把安装、bam to vcf的完整链路、参数背后的原理以及我自己踩过的坑一次说清楚从没装过的新手到想把流程跑稳的老手都能直接拿去用。1. 先弄清 bcftools 在变异检测链路里负责哪一段很多人对 bcftools 的定位是模糊的觉得它和 samtools、GATK 差不多都是处理测序数据的工具。这个理解会让你在参数选择上一直踩坑所以先花两分钟把边界划清楚。1.1 从 FASTQ 到 VCFbcftools 只吃中间那一段一条典型的重测序分析链路是这样走的FASTQ 经过比对bwa、bowtie2、minimap2 之类产出 SAMSAM 排序去重之后变成 BAM到这里都属于读段层面的工作。再往后从 BAM 里把每个位点上所有读段的碱基堆叠起来、判断这个位点有没有变异、变异是什么类型、质量如何这一层叫变异检测层bcftools 就活在这一层。它做两件核心的事前半段是bcftools mpileup把 BAM 里每个位点的碱基堆叠信息pileup算出来输出一份中间格式后半段是bcftools call拿着这份堆叠信息做统计模型判断输出 VCF。这两个子命令用管道连起来就是你到处都能搜到的那条bam to vcf经典命令。理解这个分工非常关键因为后面所有参数都归到这两类里mpileup的参数基本都在控制怎么堆叠、堆叠哪些读段、堆叠时怎么修正call的参数基本都在控制用什么模型、先验设多少、只输出变异位点还是全输出。你把参数贴在错误的子命令上bcftools 会直接报unrecognized option这也是新手最常见的第一类挫败。1.2 为什么不用 GATK几个真实的取舍点这里不讨论谁更好只讲我实际做项目时的取舍逻辑。第一个是部署成本。bcftools 是单文件可执行程序conda 装完就是一个二进制没有 Java 运行时、没有 GATK 那一堆 bundle 依赖也不存在 JVM 内存调参的问题。你在一个只有 16G 内存的服务器上跑 WGSGATK 的 HaplotypeCaller 需要小心翼翼地调-Xmx而 bcftools 基本是给多少用多少不会因为堆内存炸掉。对于需要批量部署几十个节点的场景这个差异是决定性的。第二个是流程的可解释性。bcftools mpileup | bcftools call这条链路的每一步都是透明的文本或 BCF 流你可以随时在中间插一个bcftools view看看数据长什么样。GATK 的局部重组装local reassembly在 indel 附近确实更强但它的中间状态对使用者是黑盒。做方法学验证或者需要精细控制每一步的时候透明比更强更重要。第三个是适用物种和场景。bcftools 的--ploidy参数允许你直接指定倍性二倍体、单倍体、甚至通过--ploidy-file自定义非整倍体对微生物、植物、非模式生物非常友好。而 GATK 的很多默认假设是围绕人类二倍体设计的挪到多倍体植物上要动的地方不少。一句话总结人类 WGS 追求最高灵敏度GATK 有优势细菌、植物、群体样本量大、需要自动化批处理、服务器资源受限的场景bcftools 往往是更务实的选择。我个人做微生物重测序和植物群体项目时默认就是 bcftools。1.3 版本这件事比你想的重要bcftools 和 samtools、htslib 是同一套生态三者版本必须匹配。尤其是 samtools 生成的 BAM 和 bcftools 读取 BAM 时用的 htslib 版本如果差太远会撞上一些莫名其妙的读取错误。我的做法是永远让 samtools 和 bcftools 用同一个渠道、同一批安装conda 里一条命令同时装两个让 solver 去保证版本一致。另外 bcftools 1.10 之前和之后的参数行为有变化1.12 之后 FORMAT 字段的默认输出集合也动过。你在网上抄到的老教程命令里可能省略了-a注解参数换到新版本上跑出来的 VCF 就少了AD、DP这些字段下游过滤脚本直接崩。所以我后面给的所有命令都会显式写全注解参数不依赖默认值。2. 装 bcftools 的三条路以及装完必须做的自检安装本身不复杂复杂的是装完之后环境里混着好几个版本的 htslib。2.1 conda / mamba 安装绝大多数人的最优解直接上命令mamba create -n bcf -c conda-forge -c bioconda bcftools samtools mamba activate bcf bcftools --version这里有三个细节值得单独强调因为它们决定了你会不会遇到装上了但跑不动的玄学问题。第一channel 顺序必须是 conda-forge 在前、bioconda 在后。bioconda 的包依赖 conda-forge 里的基础库如果 order 反了solver 有时候会从 bioconda 里挑一个老版本的 zlib 或 libcurl装完能跑但是压缩解压会出问题。可以在~/.condarc里锁死channels: - conda-forge - bioconda - defaults channel_priority: strictchannel_priority: strict这一行很关键它让 conda 严格按顺序取包避免混装。第二不要pip install任何和 htslib 相关的东西。Python 生态里有些包会顺手装一个自己的 htslib装完之后bcftools命令行调用的是 conda 那份但某些插件又去加载 pip 那份症状是undefined symbol之类的链接错误。你如果同时用 pysam 做下游分析最稳的做法是把 pysam 也装在同一个 conda 环境里让它们共享同一份 htslib。第三环境名别用base。不是洁癖是因为你迟早会有第二个项目需要不同版本的 bcftools全塞在 base 里最后一定互相打架。2.2 源码编译什么时候必须走这条路三种情况我会选择源码编译服务器完全没有联网权限、需要启用 conda 包没开的功能、或者需要给某个特定架构比如 ARM 服务器做适配。wget https://github.com/samtools/bcftools/releases/download/1.19/bcftools-1.19.tar.bz2 tar -xjf bcftools-1.19.tar.bz2 cd bcftools-1.19 ./configure --prefix/opt/bcf --with-libdeflate make -j 8 make install几个点解释一下。bcftools 的源码包自带 htslib不需要你先单独编译 htslib这一点和 samtools 一样很多人不知道白白折腾半天。--prefix指定安装路径装到/opt下方便多用户共享装完记得把/opt/bcf/bin加进 PATH。--with-libdeflate是启用 libdeflate 加速压缩解压如果你的流程里有大量 VCF.gz 读写这个开关能带来肉眼可见的提速大概 20% 到 40% 的 IO 时间节省。前提是系统里得先有 libdeflate 的开发包CentOS 系是libdeflate-develDebian 系是libdeflate-dev。如果 configure 报找不到就老老实实去掉这个参数别硬上。还有--with-libcurl如果你需要直接从 HTTP/FTP 读取参考基因组或者公开数据集的 VCF这个要开。但大多数内网环境用不上开了反而可能因为 libcurl 版本问题引入麻烦。编译过程中如果报zlib.h not found装zlib-devel报bzlib.h not found装bzip2-devel报lzma.h not found装xz-devel。这三个是必备依赖缺一个都编译不过。2.3 装完必做的三步自检装完别急着跑数据花两分钟做这三件事能省掉后面两小时的排查。第一步确认版本和编译信息bcftools --version # 会输出 bcftools 1.19, Using htslib 1.19注意看 htslib 那一行它告诉你当前链接的是哪个 htslib。如果这个版本号和你samtools --version里显示的 htslib 版本不一致说明环境里有多份 htslib赶紧处理。第二步确认子命令都在bcftools plugin -lplugin是 bcftools 的插件系统很多高级功能比如fill-tags、setGT都在这里。如果这个命令报错说明插件目录没装对。第三步跑一条最小命令验证管道能通bcftools view -h输出 VCF 头模板就说明基础功能正常。这一步看似废话但它能验证动态库链接是否完整——见过太多次bcftools --version能跑但view直接symbol lookup error的情况原因就是 htslib 版本错配。2.4 别忘了 samtools 那边的配对bcftools 只负责变异检测上游的 BAM 准备全靠 samtools。这两个工具在同一台机器上的版本一致性我前面说过一次这里再强调一遍具体影响samtools 1.9 生成的 BAM 用 bcftools 1.19 读绝大概率没问题因为 BAM 格式本身很稳定但如果你用了比较新的 CRAM 格式或者 BAM 里带了特殊的辅助标签跨大版本读取就可能报[E::bam_read1] the file is corrupted。统一版本是一劳永逸的解法。3. bam to vcf 全链路mpileup 和 call 到底在做什么这条管道是本文的核心我会把它拆成前置条件—堆叠阶段—判定阶段—输出阶段四块讲每一块的参数都给出选择的理由。3.1 两个前置索引缺一个都跑不起来bam to vcf需要两份索引而且这两份索引的缺失报错信息都很不直观。参考基因组需要.fai索引。用samtools faidx生成samtools faidx ref.fa # 生成 ref.fa.faimpileup要随机访问参考序列来获取每个位点的参考碱基没有.fai它就没法定位。这就是文章开头那个[E::fai_build3_core] Failed to open the file的来源。注意 bcftools 在较新版本里确实会尝试自动构建.fai但如果目录没有写权限自动构建就会失败报错信息还是那句话非常具有误导性。所以养成习惯永远手动先跑samtools faidx。BAM 需要.bai索引。用samtools index生成samtools index sample.bam这一步的前提是 BAM 已经按坐标排序。如果你是samtools sort出来的排序没问题如果是某些比对工具直接输出的未排序 BAM先排序再建索引samtools sort - 8 -o sample.sorted.bam sample.bam samtools index sample.sorted.bam关于索引格式还想补一句.bai和.csi的区别在于能覆盖的染色体长度上限。.bai只能索引长度小于 512Mb 的染色体人类、小鼠、拟南芥都没问题但如果你做的是某些超大基因组的物种需要用.csi。bcftools 从 1.x 某个版本开始对超长染色体的支持做了改进但为了保险超大基因组建议全流程用 CSI。一个快速自检的方法跑任何 mpileup 之前先确认两件事对得上——BAM 里的染色体命名和参考基因组的命名。samtools view -H sample.bam | grep ^SQ | cut -f2,3 | head cut -f1,2 ref.fa.fai | head第一列是染色体名SN第二列是长度LN。两边的名字和长度必须完全一致。我见过的真实翻车案例参考基因组来自 Ensembl染色体名是1、2、3BAM 比对时用的是 UCSC 版本名字是chr1、chr2、chr3。mpileup 跑完不报任何错输出的 VCF 里一个变异都没有因为它在 BAM 里找不到名为1的染色体读到 0 条读段自然什么也检测不出来。这种错误最难查因为流程成功了。3.2 mpileup 的参数分组过滤、修正、输出mpileup的参数可以分成三组来理解这样你就不需要死记硬背。第一组是读段和碱基的过滤阈值。最常用的是这两个-q, --min-MQ INT最小比对质量默认 0。意思是比对质量低于这个值的读段直接忽略。人类 WGS 常用 20微生物数据质量普遍好一些可以用 20 到 30。-Q, --min-BQ INT最小碱基质量默认 13。意思是单个碱基质量低于这个值的在堆叠时被当作噪声丢掉。默认 13 对应错误率约 5%是比较保守的设定一般不用动。这两个参数的区别新手最容易搞混比对质量是这条读段比在这个位置靠不靠谱碱基质量是这个位置的这一个碱基读得准不准。前者是整条读段的属性后者是单个碱基的属性。类比一下比对质量像是这本书是不是放对了书架碱基质量像是这一页上的字有没有印糊。第二组是读段堆叠时的修正选项。这里面最重要的是 BAQBase Alignment Quality碱基比对质量修正-B, --no-BAQ禁用 BAQ。默认是启用 BAQ 的。BAQ 是做什么的简单说它能识别出那些看起来比对上了其实是因为附近有 indel 导致的错位比对的碱基并把这些碱基的质量分数调低。它的作用是在 indel 附近压制假阳性 SNP。代价是计算量实测下来启用 BAQ 会让 mpileup 慢 30% 到 40%。所以在快速出结果预览和正式产出之间我的做法是预览时加-B提速正式跑一定开着。还有一个-C, --adjust-MQ INT默认 0不调整。这是老版本比对工具时代的遗留修正现代比对工具bwa-mem、minimap2一般不需要保持默认就行。如果你是拿很老的 BWA aln 结果来跑可以试试设 50。第三组是深度上限和输出格式。-d, --max-depth INT每个输入文件在每个位点最多用多少条读段默认 250最大 1000000。-a, --annotate LIST输出哪些注解字段。-Ou / -Oz / -Ov输出格式。-d这个参数值得单独说。默认 250意味着如果你测序深度有 500x微生物、扩增子、靶向捕获经常到这个深度超过 250 的那部分读段会被直接丢掉不参与统计。这会导致高覆盖区域反而检测不出变异——因为它只看了一半的数据。所以深度超过 250x 的数据一定要把-d提上去提到你实际深度的 1.5 倍左右比较保险。但也不要无脑设成 1000000因为内存占用和-d成正比设太高会把机器吃满。-a注解列表我建议永远显式写出来。常用的字段注解字段含义为什么需要FORMAT/AD每个样本的等位基因深度下游过滤和杂合判断的核心依据FORMAT/DP每个样本的总深度判断覆盖度是否足够INFO/AD所有样本汇总的等位基因深度群体层面快速筛查FORMAT/SP链偏好性统计排查假阳性过滤常见指标不写-a新版 bcftools 输出的 VCF 里就没有AD你后面想按支持变异的读段数至少 3 条来过滤就没法做。这个坑我踩过跑到过滤阶段才发现字段缺失只能从头重跑一遍 mpileup非常浪费时间。3.3 call 的参数模型选择与倍性设定call有两种调用模型-m, --multiallelic-caller多位点调用模型适合已知样本的变异检测就是我们常规做重测序时用的。-c, --consensus-caller一致性调用模型适合已知位点的基因分型比如你有一份已知变异位点列表想在新的样本里验证这些位点。绝大多数场景用-m。-v, --variants-only表示只输出有变异的位点不加这个参数会把所有位点都输出包括纯合参考位点VCF 体积能大十倍以上除非你明确需要全位点信息否则一定加-v。--ploidy是倍性设定。默认是 2二倍体。做细菌这类单倍体生物要设--ploidy 1做多倍体植物要么用--ploidy 3之类要么准备一个--ploidy-file描述非整倍体情况。这个参数设错了所有基因型判断都是错的但流程不会报错只会安静地给你一堆看起来正常的错误结果。-P, --prior FLOAT是突变率先验默认约 1.1e-3也就是假设每个位点有千分之一左右的概率是变异位点。这个值对多数物种都合适不要随便动。只有在做特殊场景比如肿瘤样本变异频率差异很大时才考虑调整。一个容易被忽略的点call也有-v之外的输出控制参数-A保留所有可能的等位基因但常规流程用不到知道有这么个东西就行。3.4 一条能跑通的完整命令以及耗时预估把上面所有东西拼起来我实际生产环境用的基础命令是这样bcftools mpileup \ -f ref.fa \ -q 20 -Q 20 -d 500 \ -a FORMAT/AD,FORMAT/DP,INFO/AD \ -Ou \ sample.bam \ | bcftools call \ -m -v \ --ploidy 2 \ -Oz -o sample.raw.vcf.gz跑完之后立刻建索引bcftools index -t sample.raw.vcf.gz关于耗时给你一个参考30x 的人类 WGS 全基因组8 线程mpileup加call全程大约 6 到 12 小时具体取决于硬盘 IO 和 BAQ 是否启用。如果只是测一个几千碱基的扩增子或者细菌基因组5Mb 左右单线程几分钟就出结果。-Ou这个输出格式的选择是有讲究的。它表示未压缩的 BCF 二进制流是管道传输最快的格式。对比一下-Ov是纯文本 VCF体积大、传输慢-Oz是 bgzip 压缩的 VCF适合落盘但不适合管道中间环节因为压缩解压本身有开销而-Ou不压缩但保持二进制体积比文本 VCF 小很多是管道中转的最优解。管道中间一律用-Ou最后落盘再用-Oz这是性能调优里最容易拿到的一个免费提升。4. 实测踩坑报错背后的真实原因这一节是我自己排错的完整记录报错信息我都保留原样方便你直接搜到。4.1 参考基因组相关的报错九成是索引问题报错长这样[E::fai_build3_core] Failed to open the file ref.fa [E::hts_open_format] Failed to open file ref.fa : No such file or directory看到这个按顺序排查三件事路径是不是相对路径mpileup 的工作目录和你写命令时所在的目录可能不一样尤其是写在脚本里被调度系统调用的时候。统一用绝对路径能避免这一类问题。目录有没有写权限bcftools 检测到没有.fai时会尝试自己建失败就报上面这个错。手动samtools faidx ref.fa跑一遍问题立刻暴露。.fai是不是过期了如果你的参考基因组文件被替换过内容但.fai还是旧的会报一些很难懂的坐标越界错误。ls -la ref.fa ref.fa.fai看时间戳.fai必须比.fa新。4.2 输出空 VCF染色体命名不一致是头号嫌疑这个坑我在 3.1 提过这里把它作为独立案例展开因为它太常见了。现象是命令跑完没有任何报错退出码 0但输出的 VCF 里只有头信息一条记录都没有或者只有零星几条。排查步骤是这样的。第一步看 BAM 里的染色体名samtools view -H sample.bam | grep ^SQ第二步看参考基因组的染色体名cut -f1,2 ref.fa.fai第三步对比。如果一边是chr1另一边是1就是这个原因。修复方式有两种要么重新比对最干净但费时要么在 mpileup 时用--rename-chrs做名称映射快但是要保证映射文件正确。--rename-chrs需要你准备一个两列的映射文件第一列是原始名第二列是目标名。顺便说还有一种情况是染色体名一样但长度不一样比如同一个物种的不同组装版本。这时候 mpileup 不会报错但坐标对不上结果全错。所以长度也必须核对。4.3 样本名变成文件名或者 unknown如果你跑完发现 VCF 的样本列名是sample.bam这种文件名而不是你期望的实际样本 ID原因在于 BAM 里缺少RG头信息或者RG里没有SM标签。检查一下samtools view -H sample.bam | grep ^RG正常的输出应该包含SM:样本ID。如果没有这行或者有RG但里面没有SMbcftools 就只能退而用文件名当样本名。这个问题的严重性在于单样本分析时无所谓但多样本联合调用joint calling时样本名错了会导致后面合并 VCF 的时候样本对不上或者出现重复样本。修复方法是在比对阶段就正确设置-R参数给 bwa 写入 RG 信息已经比对完的可以用samtools addreplacerg补但要注意补完之后需要重新排序和建索引。4.4-多线程没生效其实是管道在骗你很多人会这样写bcftools mpileup - 8 -f ref.fa sample.bam | bcftools call -m -v - 8 -Oz -o out.vcf.gz然后发现 CPU 占用很平稳根本没用满 8 核。原因有几个。第一-在 bcftools 里主要作用于压缩解压线程不是核心计算线程。mpileup的堆叠计算本质上是按染色体分块的单条染色体内部很难并行。所以你的数据如果染色体数很少比如细菌只有一条染色体多线程几乎帮不上忙。人类有 24 条主染色体多线程效果会好一些但也不是线性的。第二管道里加-Ou才是关键。如果你两边都用-Oz压缩解压的开销会吃掉多线程带来的收益而且压缩是串行的瓶颈。中间用-Ou把压缩推迟到最后一步才能真正让多线程发挥作用。第三要真正并行得按染色体拆分。做法是用bcftools mpileup -r chr1分别跑每条染色体然后bcftools concat合并。这样每个进程处理一条染色体线性加速。代价是需要临时文件管理和磁盘空间还要保证所有分片完成后正确合并。对于单次分析这个复杂度不太值得对于需要批量跑几百个样本的群体项目这个优化能省下大量时间。4.5 BAQ 带来的变异消失这个坑很隐蔽。现象是加-B禁用 BAQ跑一遍变异数量正常不加-B跑一遍某些区域的变异明显变少尤其是 indel 附近的 SNP。这不是 bug是 BAQ 在正常工作——它认为那些碱基是错位比对产生的假象所以压制了它们。但如果你的真实变异恰好就在 indel 附近BAQ 有可能误伤。我的处理方式是正式流程保留 BAQ但对于 BAQ 压制后消失的可疑位点用-B单独跑一遍对照人工确认。如果某批数据的 indel 附近变异特别重要比如做耐药位点检测很多耐药突变紧邻 indel我会考虑关闭 BAQ然后靠更严格的后续过滤来控制假阳性。这是一个典型的没有最优解只有权衡的场景。5. 从原始 VCF 到能用的 VCF归一化和过滤call出来的raw.vcf.gz只是半成品直接拿去做下游分析会有一堆问题。5.1 先 norm 再 filter顺序不能反bcftools norm做两件事左对齐left-align和拆分多等位基因位点。bcftools norm \ -f ref.fa \ -m -both \ -Oz -o sample.norm.vcf.gz \ sample.raw.vcf.gz bcftools index -t sample.norm.vcf.gz为什么要左对齐因为同一个 indel 在比对结果里可能被表示成多个位置。比如参考序列是AAAAA实际有一个 A 缺失那么缺失第一个 A、缺失第二个 A、缺失第三个 A在生物学上是同一个变异但在 VCF 里坐标不同。如果不做左对齐你和其他人的结果比对时就会对不上和已知数据库dbSNP 之类的注释也会漏掉。-f ref.fa提供参考序列norm 通过向左平移把表示统一化。为什么要拆分多等位基因一个位点如果有两个不同的 ALT 等位基因VCF 里会写成A - G,T这样的多等位记录。很多下游工具不认这种格式或者只能处理第一个 ALT。-m -both把它拆成两条独立的二态记录。-m的三个取值要记清楚参数值作用-m -拆分多等位位点SNP 和 indel 都拆-m 合并二态位点为多等位反向操作-m -both对 SNP 和 indel 都执行拆分还有一个容易被忽略的点norm之后如果做了拆分坐标会变化必须重新建索引。忘了这一步后面所有按区间查询的操作都会报索引过期或者给出错误结果。5.2 硬过滤阈值怎么定SNP 和 INDEL 分开谈硬过滤没有万能阈值但有一组经过大量实践验证的起点。对 SNP我会用这样的条件bcftools filter \ -s LowQual \ -e QUAL30 || INFO/DP10 || INFO/DP200 || MQBZ-3 \ -g 5 \ -Oz -o sample.snp.filtered.vcf.gz \ sample.norm.vcf.gz逐项解释。QUAL30是最基本的质量门槛对应大约千分之一的错误概率。INFO/DP10过滤覆盖度不足的位点10x 以下的位点在二倍体上很难可靠判断杂合。INFO/DP200过滤异常高覆盖这类位点往往是重复序列或者比对错误聚集区。MQBZ-3是比对质量偏倚的 Z 值负值大意味着变异读段的比对质量系统性低于参考读段是假阳性的典型特征。对 INDEL测试集要额外收紧bcftools filter \ -s LowQual \ -e QUAL40 || INFO/DP15 || IDV3 || IMF0.2 \ -Oz -o sample.indel.filtered.vcf.gz \ sample.norm.vcf.gzIDV是支持 indel 的读段数IMF是支持 indel 的读段比例。INDEL 的假阳性通常来自比对错误所以要求更多读段支持和更高的比例会更稳。-g 5这个参数是间隙作用是当一个变异被过滤时它周围 5bp 内的其他变异也一起标记。这是为了避免 indel 附近一串互相依赖的假阳性只被过滤掉一部分。这个参数的来历是VCF 里的连续位点往往不是独立的单独过滤某一个意义不大。关于软过滤和硬过滤的区别-s LowQual是软过滤它在 FILTER 列写一个标签但不删除记录。这样你后面还可以用bcftools view -f PASS来取回硬过滤结果保留灵活性。我强烈建议永远用软过滤因为过滤阈值是很主观的保留原始记录让你有机会回头调整。5.3 stats、isec、query三个高频配套命令bcftools stats是最快的质量体检bcftools stats sample.norm.vcf.gz sample.stats重点看几个指标SN 行里的 ts/tv转换/颠换比人类全基因组一般在 2.0 到 2.2 之间低于 1.8 说明假阳性偏多number of SNPs 和 number of indels 的比例人类一般在 4:1 到 5:1还有 singletons 数量如果单例变异比例过高可能是污染或者错误累积。配套的plot-vcfstats能把这些数字画成图plot-vcfstats -p stats_plot/ sample.statsbcftools isec用来做交集和差集做流程验证时非常有用。比如你想知道自己这条 bcftools 管道和参考结果差多少bcftools isec -p isec_dir sample.norm.vcf.gz reference.vcf.gz它会输出四个文件只在 A 里的、只在 B 里的、两者共有的、以及合并的。isec的前提是两个 VCF 都建了索引。bcftools query是提取字段的万金油bcftools query -f %CHROM\t%POS\t%REF\t%ALT\t[%GT\t%AD\t%DP]\n sample.norm.vcf.gz variants.tsv方括号里的部分是按样本重复输出的所以多样本 VCF 会为每个样本输出一组%GT %AD %DP。这个语法新手最容易漏漏了方括号就只能拿到第一个样本的值。query还支持-i条件筛选比如只取杂合位点-i GThet比先view再query效率高很多。5.4 生成样本一致序列做微生物或者小基因组的时候我经常需要把 VCF 转回 FASTA 做下游分析bcftools consensus -f ref.fa -s sample1 sample.norm.vcf.gz sample1.consensus.fa-s指定用哪个样本的基因型来构建一致序列多样本 VCF 必须写。还支持-i条件表达比如只应用高质量变异-i QUAL30这样你能生成高置信度版本和完整版本两份一致序列做对比。这个功能在病毒基因组组装、细菌基因组完成图这类场景里非常好用省掉了自己写脚本处理 VCF 的功夫。6. 让流程跑稳的几个个人习惯写脚本的时候我会在每一步之后加set -e和显式检查因为 bcftools 在某些错误场景下返回码是 0静默失败最可怕。具体做法是在管道之后检查输出文件是否存在且非空用bcftools view -H file.vcf.gz | head -1确认至少有一条记录。临时文件的命名我会带上染色体或者样本 ID避免多进程跑的时候互相覆盖。跑之前用bcftools view -h检查一下 VCF 头的contig行是否完整有些工具生成的 VCF 头里缺少 contig 定义norm的时候会报警告。还有一个几乎没人提但很实用的点保留中间文件。raw.vcf.gz、norm.vcf.gz这些不要跑完就删磁盘现在很便宜但重跑一次 mpileup 可能要好几个小时。我一般会保留 raw 和 norm 两个版本以及它们的 stats 文件出问题的时候能快速定位是哪一步的事情。最后分享一个我自己用的小脚本骨架把整个流程串起来#!/bin/bash set -euo pipefail REF$1 BAM$2 SAMPLE$3 THREADS${4:-8} samtools faidx $REF samtools index $BAM bcftools mpileup -f $REF -q 20 -Q 20 -d 500 \ -a FORMAT/AD,FORMAT/DP,INFO/AD -Ou $BAM \ | bcftools call -m -v --ploidy 2 -Oz -o ${SAMPLE}.raw.vcf.gz bcftools index -t ${SAMPLE}.raw.vcf.gz bcftools norm -f $REF -m -both -Oz -o ${SAMPLE}.norm.vcf.gz ${SAMPLE}.raw.vcf.gz bcftools index -t ${SAMPLE}.norm.vcf.gz bcftools filter -s LowQual \ -e QUAL30 || INFO/DP10 || INFO/DP200 -g 5 \ -Oz -o ${SAMPLE}.filtered.vcf.gz ${SAMPLE}.norm.vcf.gz bcftools index -t ${SAMPLE}.filtered.vcf.gz bcftools stats ${SAMPLE}.norm.vcf.gz ${SAMPLE}.statsset -euo pipefail里的pipefail是关键它让管道中任何一个环节失败都能被捕获到。默认情况下 bash 只看最后一个命令的返回码mpileup失败了但call正常退出你的脚本会以为成功了。这个参数值 10 秒钟输入能省掉无数次为什么输出是空的的困惑。跑熟之后你会发现bam to vcf这条链路真正难的不是命令本身——命令就那么长——难的是前置条件的检查和对参数行为的理解。索引齐不齐、染色体名对不对、深度上限够不够、注解字段全不全这四个问题解决了99% 的翻车都能避免。