1. 先把 bam to vcf 这条链路想明白干过几年测序分析的人大概都有个共识拿到一个比对好的 BAM 文件要把它变成能看、能筛、能注释的 VCF这一步是整个变异检测流程里最容易被低估的环节。看起来只是一条命令实际上它同时牵扯到参考基因组的版本、BAM 的质量、深度上限、并行策略、以及后面的过滤口径。而在这件事上bcftools几乎是绕不开的工具——它和 samtools、htslib 出自同一个团队命令行风格统一单文件部署管道友好跑起来不需要一堆配置文件几乎是我见过最干活的变异检测工具。我接触的很多项目里BAM 到 VCF 这段要么用 GATK HaplotypeCaller要么用 bcftools mpileup call。GATK 的局部重组装确实在 indel 上更讲究但代价是内存、时间和一堆 Java 参数bcftools 的路线是堆叠式的碱基级证据统计速度快、依赖少、脚本化极其方便对于 SNP 为主、样本量不大、或者只是想快速出一版变异列表的场景性价比高得多。尤其是它支持直接从标准输入读、往标准输出写用-Ou输出未压缩 BCF 串起管道中间不落盘这个设计在批量跑样本的时候省下的 I/O 时间非常可观。这篇内容主要面向三类人刚入门生信、第一次要把 BAM 转成 VCF 的同学手上有一批 WGS 或 panel 数据、想自己搭一套轻量 call 变异流程的分析人员以及用惯了 GATK、想换个方案做交叉验证的老手。我会从安装开始讲起把bcftools mpileup和bcftools call的每个关键参数掰开来说清楚它为什么在那儿、取值怎么定、代价是什么最后把实操中踩过的坑整理成一张排查表。2. bcftools 安装三条路线怎么选2.1 conda/mamba 一把梭最省事也最容易踩坑大多数人第一反应是 conda 装这没错但要注意两点channel 顺序和环境隔离。bioconda 依赖 conda-forge 的很多底层库如果顺序配反了依赖解析会绕得你怀疑人生。# 建一个独立环境别往 base 里塞会和系统 samtools 打架 conda create -n varcall -c conda-forge -c bioconda -c defaults \ bcftools1.19 samtools1.19 htslib1.19 conda activate varcall bcftools --version把版本号写死这件事很多人觉得没必要但我强烈建议写。原因很实际bcftools 1.15 之前的bcftools norm -m用的是-both这类写法1.17 之后改成了-any跨版本迁移脚本时这点差异足够让你浪费半天。项目开始前把版本钉住写进流程文档里后面复现的时候不会互相甩锅。如果解析依赖慢得受不了用 mamba 替代 conda 的 solvermamba create -n varcall -c conda-forge -c bioconda bcftools1.19 samtools1.19另外提一句别在同一个环境里同时装 GATK 和 bcftools 再共享同一个 htsjdk/htslib 生态历史上出现过因为 zlib、libcurl 版本冲突导致bcftools启动即崩的情况。隔离环境几秒钟的事。2.2 源码编译需要控版本、控依赖时选它有些集群不给联网、或者管理员要求所有软件统一编译到/opt下这时候源码编译是唯一选择。bcftools 的编译依赖不多但有一步容易卡住htslib。# 先把 htslib 编好 wget https://github.com/samtools/htslib/releases/download/1.19/htslib-1.19.tar.bz2 tar -xjf htslib-1.19.tar.bz2 cd htslib-1.19 ./configure --prefix/opt/htslib-1.19 make -j 8 make install # 再编 bcftools把头文件目录指过去 cd .. 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/bcftools-1.19 --with-htslib/opt/htslib-1.19 make -j 8 make install export PATH/opt/bcftools-1.19/bin:$PATH export LD_LIBRARY_PATH/opt/htslib-1.19/lib:$LD_LIBRARY_PATH这里有几个经验值。--with-htslibsystem只在系统装了 htslib 开发包libhts-dev之类时才成立很多服务器的包管理器给的 htslib 版本偏老编出来会出现编译通过但跑起来提示插件不匹配的怪问题所以宁可自己编一份。另外如果目标机器是别的架构比如 ARM 的服务器make之前记得检查./configure输出里有没有正确识别到-march相关的优化跨架构交叉编译是需要显式指定--host的。提示编译完之后别急着删源码目录。后面排查问题时经常要回去看一眼bcftools --version里带的 htslib 版本号以及bcftools plugin -l能不能列出插件。2.3 系统包和容器镜像适合够用就行的场景apt或yum直接装是最快的sudo apt-get update sudo apt-get install -y bcftools tabix bcftools --version | head -1但要有心理预期发行版仓库里的版本往往落后两三年Ubuntu 20.04 给的是 1.10 系列很多新参数比如 gVCF 输出相关的选项根本没有。如果只是做一个教学演示或者处理很小规模的数据用系统包完全可以一旦进入正式项目还是建议 conda 或源码。容器方案适合在集群上跑流程、需要保证环境完全一致的场景。用 biocontainers 提供的镜像注意 tag 里带有构建哈希直接抄别人的 tag 大概率拉不到去镜像列表里挑当前最新的那个docker run --rm -v $(pwd):/data -w /data \ quay.io/biocontainers/bcftools:1.19--h8b25389_0 \ bcftools --versionSingularity/Apptainer 在 HPC 上更常见用法就是把docker://换成oras://或者下载好的.sif文件其余一样。2.4 装完必须做的三件事第一件确认主程序、htslib、插件三个层面都没问题bcftools --version bcftools plugin -l samtools --version | head -2第二件确认bgzip和tabix也在 PATH 里。bcftools 的索引功能依赖它们虽然 bcftools 自己也有bcftools index但很多下游工具IGV、部分 R 包只认.tbi所以两个都要有。第三件在你的流程脚本里加一条版本断言。不然半年后回头复现数据发现结果对不上第一件事就是怀疑版本而那时候你已经记不清当时装的是哪个了。3. 上机之前BAM 与参考基因组的准备工作3.1 BAM 必须满足的四个硬条件bcftools mpileup对输入 BAM 的要求写在文档里但实际踩坑的时候往往不是报错而是跑通了但结果不对所以这四条值得单独拿出来说。第一坐标排序。必须是samtools sort之后的坐标排序 BAM不是 name-sorted。判断方法samtools quickcheck -v sample.bam # 不输出任何东西才算通过 samtools view -H sample.bam | head -5第二必须建索引。.bai或.csi都行-r区域参数完全依赖它samtools index sample.bam ls -lh sample.bam.bai第三SQ头信息必须和参考基因组完全一致。这里的一致性不只是名字还包括长度。名字对了长度错了bcftools 一样会报错或者悄悄跑出错误结果。第四read group 要规范。单样本场景下很多人忽略这一点但如果后面要用bcftools mpileup处理多个 BAM它会按SMsample标签把多个文件合并成同一个样本列SM缺失或写错结果里就会出现奇怪的样本名。注意不要拿已经被 hard-clip 处理过的 BAM 去做 indel 检测。很多比对软件在末端会做 soft-clip这是正常的但如果你的上游流程做过激进的 hard-clipindel 附近的证据已经丢了bcftools 再厉害也补不回来。3.2 参考基因组索引一个字符都不能差bcftools mpileup的-f参数要求参考序列有.fai索引bcftools norm的左对齐也依赖它samtools faidx ref.fa head -3 ref.fa.fai这里最经典的坑是contig 命名不一致。UCSC 下载的参考是chr1、chr2Ensembl 下载的是1、2而你的 BAM 用的是哪一套取决于比对时用的参考。两者不匹配时bcftools mpileup要么直接报错找不到染色体要么在-t流式模式下默默跑完但一条记录都不输出——后者更折磨人因为你以为跑成功了。先做一次对照检查# BAM 头里的 contig 名 samtools view -H sample.bam | awk /^SQ/{print $2} | sed s/SN:// | head # 参考里的 contig 名 cut -f1 ref.fa.fai | head如果确实不一致用 bcftools 自带的改名功能# rename.txt 两列旧名 新名 cat rename.txt EOF 1 chr1 2 chr2 EOF bcftools annotate --rename-chrs rename.txt -Oz -o renamed.vcf.gz in.vcf.gz如果要在 BAM 层面改名可以用samtools reheader配合samtools view -H重写头信息但改完必须重新建索引这一点经常被忘掉。3.3 read group 与多样本场景单样本直接跑就行多样本要注意两件事。一是每个 BAM 的SM必须唯一且正确二是bcftools mpileup的多文件输入会把它们当同一个样本的多个文库合并处理而不是当成多个独立样本。真正要叫多个样本各自的基因型正确姿势是先各自 call 出 VCF再用bcftools merge合并或者一开始就用-s指定样本子集。补齐 read group 可以用samtools addreplacerg \ -r ID:lib1 -r SM:sampleA -r LB:lib1 -r PL:ILLUMINA \ -o rg_sampleA.bam sampleA.bam samtools index rg_sampleA.bam加上LBlibrary标签是个好习惯后面做 duplicate marking 或者排查文库特异性偏倚时用得上。PLplatform不写bcftools 会默认按 ILLUMINA 处理用 PacBio 或 ONT 数据的时候最好显式写上因为不同平台的错误模型不一样。4. mpileup call一条命令拆成五段讲4.1 第一段先看帮助再写命令这一步我建议每次都做尤其是换了一台机器、换了一个版本之后bcftools mpileup 21 | head -60 bcftools call 21 | head -40原因很直接-CBAQ 相关的 mapping quality 调整这个参数在不同版本里的默认值处理不完全一样文档里写的和--help里[方括号]中显示的默认值才是真正生效的。花三十秒确认一下能省掉后面半天的困惑。4.2 第二段mpileup 生成 BCF 中间结果完整的核心命令长这样bcftools mpileup \ -f ref.fa \ -q 20 \ -Q 20 \ -d 1000 \ -C 50 \ -a FORMAT/AD,FORMAT/DP,FORMAT/SP \ -Ou \ -r chr1 \ sample.bam \ | bcftools call -m -v -P 1e-3 -Oz -o chr1.raw.vcf.gz逐段说。-f ref.fa指定参考必须有.fai。-q 20是比对质量下限低于 20 的 reads 直接不参与堆叠。-Q 20是碱基质量下限低于 20 的碱基按 N 处理。-d 1000把单碱基位点的最大深度上限抬到 1000默认是 250这个改动在深度较高的数据上非常关键。-a指定要输出的注释字段。-Ou输出未压缩的 BCF 直接进管道。-r chr1只处理 1 号染色体。这里有个值得展开的细节-a后面能写的字段很多FORMAT/AD等位基因深度、FORMAT/DP样本深度、FORMAT/SP链偏好是最常用的三个。文档里明确提到FORMAT/AD和INFO/AD不能同时请求因为两者语义重复、实现路径不同。如果你只要位点总深度用INFO/DP就够了能省下不少文件体积。具体支持哪些字段还是以本地bcftools mpileup --help的列表为准。4.3 第三段call 把堆叠证据变成基因型bcftools call是模型层它把 mpileup 输出的位点碱基计数用贝叶斯方法转成这个位点是不是变异、基因型是什么。bcftools call -m -v -P 1e-3 -Oz -o out.vcf.gz in.bcf-m表示使用多等位基因模型multiallelic caller它会同时考虑 SNP 和 indel比老的-c一致调用模型更适合一般的 germline 场景。-v表示只输出变异位点如果要做 gVCF 用于后续 joint calling就不能加-v同时还要用 gVCF 相关选项新版 bcftools 已经支持-g输出 gVCF具体行为建议对着版本帮助确认。-P 1e-3是位点变异先验概率默认值是约 1.1e-3把它调大比如 1e-2会在低深度数据上提高灵敏度但同时假阳性也会上升属于典型的灵敏度/特异度权衡旋钮。4.4 第四段norm 做左对齐和多等位拆分这一步很多人省掉了结果下游比对时发现同一个 indel 在不同样本里位置不同合并、取交集全乱套。bcftools norm \ -f ref.fa \ -m -any \ -Oz -o chr1.norm.vcf.gz \ chr1.raw.vcf.gz-f ref.fa是必须的norm用它来把 indel 左对齐并规范化表示。-m -any表示拆分所有多等位位点让每个 ALT 独立成行。老版本 bcftools 里这个参数写作-m -both如果你的版本比较老报错说参数不识别换成-both试试就对了。如果还想去掉完全重复的记录加上-d exact。提示norm的执行顺序建议是先 norm 再过滤。因为拆分多等位之后每个 ALT 的深度和频率才是有意义的独立数值用未拆分的记录去写过滤表达式AD数组的下标对应关系很容易搞错。4.5 第五段过滤、索引、统计# 硬过滤 bcftools view \ -i QUAL20 INFO/DP10 \ -Oz -o chr1.pass.vcf.gz \ chr1.norm.vcf.gz # 建索引-t 生成 .tbi不加 -t 默认生成 .csi bcftools index -t chr1.pass.vcf.gz # 出统计报告 bcftools stats chr1.pass.vcf.gz chr1.stats.txt grep ^SN chr1.stats.txt看一眼变异总数、SNP/indel 比例、转换/颠换比Ts/Tv。人类全基因组 germline 数据的 Ts/Tv 通常落在 2.0 到 2.1 之间如果明显偏低比如 1.5往往意味着假阳性偏多或者过滤太松如果高到 2.5 以上可能是过滤太狠把真 indel 都杀掉了。这个数值是一个非常廉价的结果是否合理的体检指标建议每次都看。5. 关键参数逐个掰开含义、取值与代价5.1 -q 和 -Q质量阈值的边界在哪里这两个参数看起来最简单实际最容易设错。-q是比对质量MAPQ阈值-Q是碱基质量BQ阈值。对比对质量来说20 意味着至少有三成把握这个 read 没比错位置对于唯一比对的 read 来说这个门槛很低大多数 read 都能过。真正需要提高-q的场景是重复区域多、或者用了允许大量多重比对的比对策略这时候-q 30甚至-q 40能把大量模糊比对过滤掉。但代价也很明显重复区、着丝粒附近的覆盖率会明显下降那里的变异基本就检不出来了。所以如果你的目标区域里有低复杂度序列-q千万别一刀切设太高考虑用 BED 文件分区处理。碱基质量阈值-Q的情况类似。默认值通常在 13 附近设到 20 是个常见起点。但要注意如果上游比对时已经做过 BQSR 之类的碱基质量重校准再叠加一个 20 的硬阈值有可能过度过滤。实际做法是先用samtools stats看一眼碱基质量分布再决定阈值落在哪个分位点上。5.2 -d覆盖度上限这个坑比想象中深这是我认为最容易造成结果静默出错的参数。bcftools mpileup对每个位点的深度有上限默认值在很多版本里是 250。如果你的数据是 30x 全基因组这个上限基本不会触发但如果做的是 panel、amplicon 或者靶向捕获某些位点的原始深度可能到几千甚至上万这时候超过 250 的部分会被直接丢弃。后果不是报错而是深度越高的位点被截断得越厉害等位基因频率的估计就越偏。对于一个真实频率 30% 的变异如果只统计了前 250 条 read而这个前 250 条的顺序又和 read 名也就是位置相关估计值可能偏出去好几个百分点甚至因为等位基因 read 数量低于阈值而直接漏检。我的经验取值场景典型深度建议 -d 取值全基因组 30x30500 到 1000全外显子60 到 1501000 到 2000靶向 panel500 到 20005000 到 10000amplicon 超深度上万50000 以上或改用其他策略把-d设得很大不会让程序崩溃只会让内存占用上去因为每个位点要保留的中间数据结构变大了。所以宁可设大一点也不要卡在默认值上。5.3 -C 和 -BBAQ 到底开不开BAQBase Alignment Quality是一个把 indel 附近碱基的比对质量人为调低的机制。它的初衷是好的indel 附近的比对位置往往有歧义一个真的 SNP 可能只是因为附近的 indel 导致比对偏了一位看起来像变异其实是假的。BAQ 通过调低这些碱基的质量让它们对变异判定的贡献降低。-C就是设置这个调整的强度常见推荐值是 50。设成 0 就相当于关闭调整。-B则是直接禁用 BAQ 计算。什么时候开、什么时候关SNP 为主、且样本里 indel 不多-C 50是好选择能明显降低 indel 附近的假阳性 SNP。本来就是要找 indelBAQ 会把 indel 附近的证据削弱可能造成真 indel 漏检。这时候可以试试-C 0或者--redo-BAQ对比一下结果差异。跑得慢、要提速-B直接跳过 BAQ 计算速度提升很明显但要接受假阳性上升。我个人的做法是先用小型测试数据比如一条染色体跑三组对比-C 50、-C 0、-B看差异位点的数量级和分布特征再决定正式跑用哪一组。这个对比实验花半小时能避免整批数据返工。5.4 -a注释字段怎么选直接影响下游-a决定 VCF 里带哪些信息。带多了文件巨大、读写变慢带少了后面想筛都没得筛。几个常用字段的用途字段含义什么时候需要FORMAT/AD每个样本各等位基因的深度做频率过滤、杂合/纯合判断FORMAT/DP每个样本的总深度按样本深度过滤FORMAT/SP链偏好排查链特异性假阳性INFO/DP位点总深度位点级快速过滤INFO/AD位点各等位深度与FORMAT/AD二选一我的习惯是至少带上FORMAT/AD和FORMAT/DP因为后面几乎所有过滤表达式都绕不开它们。FORMAT/SP在排查特定的假阳性时非常有用属于平时占地方、关键时刻救命的字段如果存储不是瓶颈建议带上。至于INFO/AD和FORMAT/AD前面说过两者不能同时要我一般选FORMAT/AD因为样本级的深度信息信息量更大。5.5 call 侧的 -m、-v、-P 与 --ploidycall的参数不多但每一个都影响结果形态。-m建议默认开着尤其是有 indel 需求的时候。-v决定是只输出变异位点还是输出所有位点做 gVCF 或者需要无损记录时要关掉。-P是变异先验概率。它的作用可以这么理解它告诉模型在没有任何数据的情况下我认为一个位点是变异的概率有多大。默认值约 1.1e-3 是基于人类基因组的经验统计。如果你的物种或数据来源差异很大比如细菌、植物重测序这个先验不一定合适可以适当调整。--ploidy是二倍体以外的场景必须设的。默认是二倍体但性染色体、线粒体、以及一些多倍体物种都不是二倍体。bcftools 内置了一些预设比如针对特定参考基因组的 ploidy 设置能处理 X/Y 染色体的性别差异也可以自己写 ploidy 文件# ploidy.txt 格式示意 X 1 2 Y 1 1 MT 1 1第一列染色体名第二列是样本序号或性别标记第三列是倍性。这个文件的准确格式在不同版本间有细微差异写完之后用一条小命令试跑一下确认没有报错再用到全量数据上。6. 常见报错与结果异常排查手册6.1 索引、路径与格式类报错这一类最好排查因为程序会直接把原因说出来。报错片段根本原因处理方式Could not retrieve index fileBAM 没索引samtools index x.bamFailed to open file ref.fa路径错或文件不在当前工作目录用绝对路径检查权限[E::fai_build3_core] Failed to open参考没有.faisamtools faidx ref.faunknown file type扩展名和实际格式不符用bcftools view -h看一眼头not compressed with bgzip用了普通 gzip 压缩 VCF用bgzip重新压或直接-Oz生成最后一条特别值得说。普通gzip压出来的.vcf.gz看起来和bgzip压的一模一样文件名也一样但它是块不可随机访问的tabix和bcftools index都会拒绝。判断方法很简单file yourfile.vcf.gz # bgzip 的结果里会提到 BGZF6.2 参考一致性类报错这一类最隐蔽。程序可能不报错只是结果里一条记录都没有或者报一个语焉不详的chromosome not found。排查顺序是先比 contig 名字再比 contig 长度最后比参考序列本身。# 名字和长度一起比 samtools view -H sample.bam | awk /^SQ/{print $2\t$3} awk {print SN:$1\tLN:$2} ref.fa.fai如果长度对不上那说明 BAM 和参考压根不是同一套别想着用改名糊过去老老实实回退到比对步骤重做。如果长度一样但叫法不同chr1vs1改名就能解决。6.3 跑完了但 VCF 几乎是空的这是最让人抓狂的一类。按可能性从高到低排第一区域参数写错。-r chr1写成-r chr01或者-r 1而参考里叫chr1结果就是零条记录。注意-r在小版本间的行为差异有些版本对不存在的区域会直接报错有些会静默返回空。用-r之前先grep一下参考的.fai。第二质量阈值太狠。-q 30 -Q 30在低质量数据上能把几乎所有证据滤掉。处理办法是把阈值降下来做一次对照跑看记录数是不是上来了。第三-d或者深度相关参数设置和实际数据不匹配。少见但确实遇到过尤其是深度极低比如 1x 到 3x的数据堆叠证据不足以触发变异判定。第四参考基因组选错了版本。BAM 是 hg19 的参考给了 hg38坐标能对上序号但对不上位置结果就是一堆噪声位点或者一片空白。6.4 假阳性偏多的几种典型原因Ts/Tv 明显偏低、indel 数量远超预期、或者同一位置在多个样本里都出现变异通常指向这几个原因BAQ 没开或者-C设成 0indel 附近的假 SNP 大量涌入。重复序列区域没有屏蔽。BAM 里如果还带着 PCR 重复而没有标记这些重复 read 会人为堆高支持变异的 read 数。用samtools markdup处理完再 call效果立竿见影。-d截断导致频率估计失真等位基因深度比例被扭曲某些本该被判为纯合的位点被叫成杂合。read group 缺失或样本串了。多个样本的 BAM 用了同一个SMmpileup 把它们合并成一个样本导致深度异常高、等位深度异常大。排查假阳性最有效的手段是拿几个已知的真变异位点去 VCF 里查一下看它们的 QUAL、AD、DP 落在什么区间然后反过来定过滤阈值。这比套用别人的过滤参数靠谱得多。7. 效率优化大样本怎么跑才不熬人7.1 分区并行然后用 concat 合回去bcftools mpileup对单条命令的多线程支持有限bcftools view的--threads主要用于压缩解压不是计算并行所以大批量数据最实际的提速方式是按染色体或窗口切分多进程并行最后合并。# 按染色体列表并行 cut -f1 ref.fa.fai chroms.txt cat chroms.txt | xargs -P 8 -I {} bash -c bcftools mpileup -f ref.fa -q 20 -Q 20 -d 1000 -C 50 \ -a FORMAT/AD,FORMAT/DP -Ou -r {} sample.bam \ | bcftools call -m -v -Oz -o {}.raw.vcf.gz bcftools index -t {}.raw.vcf.gz # 合并-a 允许区域之间有重叠 bcftools concat -a -Oz -o all.vcf.gz $(cut -f1 ref.fa.fai | sed s/$/.raw.vcf.gz/) bcftools index -t all.vcf.gz两个注意点。一是concat要求每个输入文件都是 bgzip 压缩并且已经建好索引这一点经常被漏掉报错信息也不会直接告诉你缺索引。二是concat会做重复区域的检查如果切分窗口之间有重叠不加-a会报错退出。7.2 管道与中间格式的选择-Ou未压缩 BCF和-Oz压缩 VCF之间的选择对性能影响很大。原则是管道内部一律用-Ou只在最后落盘的时候用-Oz。未压缩 BCF 是二进制格式读写都不需要做压缩解压速度快得多而且体积也不大。如果中间结果一定要落盘优先考虑-Ob压缩 BCF它比压缩 VCF 更紧凑且读的时候更快。7.3 资源估算的一点经验内存占用大致和同时处理的最长区域 × 平均深度成正相关。一条 250 Mb 的染色体、30x 深度单次 mpileup 的内存峰值通常在几 GB如果深度是 200x同样的染色体可能要吃几十 GB。所以在 HPC 上用-r按染色体切分之前先看一眼最长的那条染色体有多长、数据深度多少别一上来就把整个基因组丢进去。如果遇到超长染色体比如某些植物基因组单条染色体几百 Mb把窗口切到 5 到 10 Mb 一段更稳。切分的时候注意窗口之间留 0 重叠即可因为变异位点不会跨窗口norm的左对齐是在单个文件内做的跨窗口的 indel 可能会因为对齐位置不同而在合并时出现重复concat -a会做去重问题不大。8. 下游衔接从原始 VCF 到能用的结果8.1 过滤表达式怎么写bcftools view -i和bcftools filter -e是两套逻辑。view -i是保留满足条件的filter -e是剔除满足条件的后者还可以配合-s给被剔除的记录打一个软标签而不删除。# 硬过滤直接删除不合格记录 bcftools view -i QUAL20 INFO/DP10 MQ30 -Oz -o pass.vcf.gz raw.vcf.gz # 软过滤保留记录但打标签 bcftools filter \ -s LowQual -e QUAL20 || INFO/DP10 \ -g 5 \ -Oz -o soft.vcf.gz raw.vcf.gz-g 5这个参数值得单独提它会过滤掉距离 indel 5 bp 以内的 SNP因为这些位置的 SNP 有很大概率是 indel 附近的比对歧义造成的假阳性。这是降假阳性非常有效的一招代价是会丢掉一些真实的紧密连锁变异。8.2 注释与统计原始 VCF 只是一堆坐标和基因型要变成能解读的结果还需要注释和统计。# 统计报告 bcftools stats all.vcf.gz all.stats.txt plot-vcfstats -p stats_plots all.stats.txt # 常用查询只看关键列 bcftools query -f %CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%INFO/DP\n all.vcf.gz | head # 按样本拆出基因型矩阵 bcftools query -f %CHROM\t%POS[\t%GT]\n all.vcf.gz genotype_matrix.txt注释一般交给专门的工具snpEff 或 VEP 做功能注释加上 dbSNP 之类的已知变异库做 ID 回填。回填用bcftools annotatebcftools annotate \ -a dbSNP.vcf.gz \ -c ID \ -Oz -o annotated.vcf.gz \ all.pass.vcf.gz版本匹配这里又是一个坑dbSNP.vcf.gz的染色体命名必须和你的 VCF 一致坐标版本也要一致。不匹配的结果不是报错而是几乎所有的 ID 都填不上ID列全是.。跑完记得bcftools query -f %ID\n | grep -v ^\.$ | wc -l数一下有多少条被填上了比例太低就是版本对不上。8.3 多套结果的合并与比较不同样本、不同 call 流程出来的 VCF 经常需要放一起看。# 合并多个单样本 VCF-0 把缺失填充为 0/0 bcftools merge -0 -Oz -o merged.vcf.gz sampleA.vcf.gz sampleB.vcf.gz # 取交集-n2 表示两套结果都有的位点 bcftools isec -n2 -w1 setA.vcf.gz setB.vcf.gz common.vcfmerge要求所有输入在同一套参考坐标系下并且样本名不冲突。isec的逻辑比较绕-n2表示至少在两套结果里出现的位点配合-w1指定输出格式属于哪一套结果的文件。建议先用小数据集跑一遍看清楚输出结构再上大批量。9. 我个人在实际操作中的几点体会跑这套流程跑得多了慢慢会形成一些不太写在文档里的习惯。一个是先小后大。任何一批新数据我都会先挑一条染色体、甚至一段一 Mb 的区域跑通整条链路看一眼 Ts/Tv、SNP/indel 比例、平均深度这些宏观指标合不合理然后再放全量。全量跑一次可能几个小时中途发现参数不对重来一次代价太大。另一个是参数写进脚本不写在命令行历史里。所有用到的参数都固化成一个可执行的 shell 脚本或者配置文件包括版本号。命令行敲的那一次只是调试正式跑的一定是脚本。再一个就是不要迷信过滤阈值。网上流传的QUAL20、DP10这类阈值都是从别人的数据里总结出来的你的建库方式、测序平台、比对参数都不一样直接套用大概率不是最优。花点时间拿真阳性对照样本如果有算一下灵敏度比什么都强。最后一个小技巧bcftools mpileup和bcftools call之间的管道如果中断你只会看到一个不完整的输出文件程序不一定报错。所以批量跑的时候务必在脚本里加上set -o pipefail让管道里任何一个环节失败都能被捕获到不然你会收获一批看起来跑完了但其实是半成品的 VCF。