企业官网建设流程全解析

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 \ bcftools=1.19 samtools=1.19 htslib=1.19 conda activate varcall bcftools --version

把版本号写死这件事，很多人觉得没必要，但我强烈建议写。原因很实际：bcftools 1.15 之前的bcftools norm -m用的是-both这类写法，1.17 之后改成了-any，跨版本迁移脚本时这点差异足够让你浪费半天。项目开始前把版本钉住，写进流程文档里，后面复现的时候不会互相甩锅。

如果解析依赖慢得受不了，用 mamba 替代 conda 的 solver：

mamba create -n varcall -c conda-forge -c bioconda bcftools=1.19 samtools=1.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-htslib=system只在系统装了 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 --version

Singularity/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，它会按SM（sample）标签把多个文件合并成同一个样本列，SM缺失或写错，结果里就会出现奇怪的样本名。

注意：不要拿已经被 hard-clip 处理过的 BAM 去做 indel 检测。很多比对软件在末端会做 soft-clip，这是正常的；但如果你的上游流程做过激进的 hard-clip，indel 附近的证据已经丢了，bcftools 再厉害也补不回来。

3.2 参考基因组索引：一个字符都不能差

bcftools mpileup的-f参数要求参考序列有.fai索引，bcftools norm的左对齐也依赖它：

samtools faidx ref.fa head -3 ref.fa.fai

这里最经典的坑是contig 命名不一致。UCSC 下载的参考是chr1、chr2，Ensembl 下载的是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

加上LB（library）标签是个好习惯，后面做 duplicate marking 或者排查文库特异性偏倚时用得上。PL（platform）不写，bcftools 会默认按 ILLUMINA 处理，用 PacBio 或 ONT 数据的时候最好显式写上，因为不同平台的错误模型不一样。

4. mpileup + call：一条命令拆成五段讲

4.1 第一段：先看帮助，再写命令

这一步我建议每次都做，尤其是换了一台机器、换了一个版本之后：

bcftools mpileup 2>&1 | head -60 bcftools call 2>&1 | head -40

原因很直接：-C（BAQ 相关的 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 'QUAL>=20 && INFO/DP>=10' \ -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设得很大不会让程序崩溃，只会让内存占用上去，因为每个位点要保留的中间数据结构变大了。所以宁可设大一点，也不要卡在默认值上。

5.3 -C 和 -B：BAQ 到底开不开

BAQ（Base Alignment Quality）是一个"把 indel 附近碱基的比对质量人为调低"的机制。它的初衷是好的：indel 附近的比对位置往往有歧义，一个真的 SNP 可能只是因为附近的 indel 导致比对偏了一位，看起来像变异，其实是假的。BAQ 通过调低这些碱基的质量，让它们对变异判定的贡献降低。

-C就是设置这个调整的强度，常见推荐值是 50。设成 0 就相当于关闭调整。-B则是直接禁用 BAQ 计算。

什么时候开、什么时候关：

• SNP 为主、且样本里 indel 不多：-C 50是好选择，能明显降低 indel 附近的假阳性 SNP。
• 本来就是要找 indel：BAQ 会把 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/AD和FORMAT/AD，前面说过两者不能同时要，我一般选FORMAT/AD，因为样本级的深度信息信息量更大。

5.5 call 侧的 -m、-v、-P 与 --ploidy

call的参数不多，但每一个都影响结果形态。

-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 索引、路径与格式类报错

这一类最好排查，因为程序会直接把原因说出来。

最后一条特别值得说。普通gzip压出来的.vcf.gz看起来和bgzip压的一模一样，文件名也一样，但它是块不可随机访问的，tabix和bcftools index都会拒绝。判断方法很简单：

file yourfile.vcf.gz # bgzip 的结果里会提到 BGZF

6.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设成 0，indel 附近的假 SNP 大量涌入。
• 重复序列区域没有屏蔽。BAM 里如果还带着 PCR 重复而没有标记，这些重复 read 会人为堆高支持变异的 read 数。用samtools markdup处理完再 call，效果立竿见影。
• -d截断导致频率估计失真，等位基因深度比例被扭曲，某些本该被判为纯合的位点被叫成杂合。
• read group 缺失或样本串了。多个样本的 BAM 用了同一个SM，mpileup 把它们合并成一个样本，导致深度异常高、等位深度异常大。

排查假阳性最有效的手段是拿几个已知的真变异位点去 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 'QUAL>=20 && INFO/DP>=10 && MQ>=30' -Oz -o pass.vcf.gz raw.vcf.gz # 软过滤：保留记录但打标签 bcftools filter \ -s LowQual -e 'QUAL<20 || INFO/DP<10' \ -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 annotate：

bcftools 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 # 取交集，-n=2 表示两套结果都有的位点 bcftools isec -n=2 -w1 setA.vcf.gz setB.vcf.gz > common.vcf

merge要求所有输入在同一套参考坐标系下，并且样本名不冲突。isec的逻辑比较绕，-n=2表示"至少在两套结果里出现"的位点，配合-w1指定输出格式属于哪一套结果的文件。建议先用小数据集跑一遍看清楚输出结构，再上大批量。

9. 我个人在实际操作中的几点体会

跑这套流程跑得多了，慢慢会形成一些不太写在文档里的习惯。

一个是"先小后大"。任何一批新数据，我都会先挑一条染色体、甚至一段一 Mb 的区域跑通整条链路，看一眼 Ts/Tv、SNP/indel 比例、平均深度这些宏观指标合不合理，然后再放全量。全量跑一次可能几个小时，中途发现参数不对重来一次代价太大。

另一个是"参数写进脚本，不写在命令行历史里"。所有用到的参数都固化成一个可执行的 shell 脚本或者配置文件，包括版本号。命令行敲的那一次只是调试，正式跑的一定是脚本。

再一个就是不要迷信过滤阈值。网上流传的QUAL>20、DP>10这类阈值都是从别人的数据里总结出来的，你的建库方式、测序平台、比对参数都不一样，直接套用大概率不是最优。花点时间拿真阳性对照样本（如果有）算一下灵敏度，比什么都强。

最后一个小技巧：bcftools mpileup和bcftools call之间的管道如果中断，你只会看到一个不完整的输出文件，程序不一定报错。所以批量跑的时候，务必在脚本里加上set -o pipefail，让管道里任何一个环节失败都能被捕获到，不然你会收获一批"看起来跑完了"但其实是半成品的 VCF。