software

BCFtools 实跑|拆开一个多等位位点,为什么同一批读段能算出两种变异比例?

作者:发布于 2026-11-14约 6 分钟阅读
拆分位点,分母去哪了。样本原有 A、C、G 三类读段,分别为 3、7、11 条。默认拆分,C 行为 3,7,G 的 11 条不再出现在这行的 AD 中。保留总和,C 行为 14,7,14 包含 A 和 G,目标 C 仍是 7 条。其他等位,0 或缺失,不同编码改变下游如何解释一条拷贝。

一个样本在同一位置测到了 A、C、G 三种碱基,其中支持 A 的有三条,支持 C 的有七条,支持 G 的有十一条。要问 C 占所有已分配碱基的几成,手算是七除以二十一。把多等位 VCF 拆成每种替代等位各一行以后,C 那行却只剩三和七两个深度值。此时把七除以两数之和就得到百分之七十。两个结果都可以从表里正确算出,但分母已经变了。我们用 BCFtools 对同一条 VCF 做三次拆分,同时追踪深度和基因型怎样被重新编码。

拆分前先读懂四个符号

VCF 以参考等位为零号,后面的替代等位按 ALT 列中的顺序从一开始编号。本例参考为 A,替代依次是 C 和 G。第一个样本基因型为 1/2,两条拷贝分别携带 C 和 G,斜线表示没有在字段中指定相位。第二个样本为 0/2,携带 A 和 G。格式头还声明等位深度依照参考与所有替代等位的顺序排列,所以样本一的三个数对应 A、C、G,样本二的三个数则为十、零、八。位点的序列名和位置相同,样本层面的证据仍要分开读。VCF 4.3 对基因型与等位深度的定义

深度字段有两个常见的分母来源。一个是输出等位深度向量自身的总和,另一个是样本总深度字段。真实调用器可能因为读段质量、重叠或等位归属规则,使两者在原始文件中就不完全相等。我们特意构造了两者相等的输入,这样拆分后出现的差值只能来自本次转换。生产数据应回查调用器对字段的定义,并在原始记录上先做一次核算。只有清楚哪些观测能进入分母,事后才知道这个比例能回答什么。

实跑经验|看着 C 行里的三个数怎样移动

安装 BCFtools 1.19,使用它的规范化命令把多等位记录拆成两条。实验的输入只有一个位置、两个样本,包含基因型、逐等位深度和总深度。第一次使用默认拆分,第二次增加保持等位深度总和的选项,第三次要求在拆分基因型时把另一个替代等位记为缺失。三个输出都重新用查询命令按样本导出,直接读取程序写下的字段值。BCFtools norm 官方选项说明

样本一的输出方式 C 行的 GT/AD/DP G 行的 GT/AD/DP
默认拆分 1/0,3,7,21 0/1,3,11,21
保持 AD 总和 1/0,14,7,21 0/1,10,11,21
其他替代等位设缺失 1/.,3,7,21 ./1,3,11,21
拆分位点,分母去哪了。样本原有 A、C、G 三类读段,分别为 3、7、11 条。默认拆分,C 行为 3,7,G 的 11 条不再出现在这行的 AD 中。保留总和,C 行为 14,7,14 包含 A 和 G,目标 C 仍是 7 条。其他等位,0 或缺失,不同编码改变下游如何解释一条拷贝。
样本原有 A、C、G 三类读段,分别为 3、7、11 条。

默认拆分的 C 行将原始参考深度三和目标 C 深度七保留,另一种替代等位 G 的十一条不在该行的等位深度向量中。总深度仍为二十一。第二次运行要求保持总和,C 行的参考栏变为十四,恰好是原来 A 的三加上 G 的十一。G 行的参考栏则变成十,等于 A 的三加 C 的七。这一操作保住了本例每条拆分记录的深度总数,却改变了输出参考栏的生物学解释。在需要报告原始参考碱基支持量时,不应把十四写成支持 A 的读段数。

第二个样本提供了更直观的检查。它原来有十条 A、零条 C、八条 G。默认拆分后,C 行是参考纯合的编码,等位深度十和零,总深度十八。保持总和后,C 行变成十八和零,而输入中真正支持 A 的仍只有十条。若面对的是要筛选“携带 C 的样本”,其结论恰好相同。若写成“十八条都观测到参考碱基 A”,证据就不支持。样本二的 G 行保留十和八,因为被去掉的 C 在这个样本里计数刚好为零。

bcftools norm -m -any -Ov -o split.vcf in.vcf
bcftools norm -m -any --keep-sum AD -Ov -o sum.vcf in.vcf
bcftools norm -m -any --multi-overlaps . -Ov -o missing.vcf in.vcf
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT:%AD:%DP]\n' split.vcf

另一个容易漏掉的细节是两个输出行都带着二十一的总深度。把 C 行和 G 行相加得到四十二,是把同一批观测重复计算了两次。若要统计原位点总深度,回到拆分前的那一行读取。若要统计每种替代的支持数,则分别取七和十一,并记住这两份证据来自同一个样本、同一个位置。对多等位位点很多的队列,这个差异会在按基因汇总时累积。

缺失、参考和“其他等位”该怎样传给下游

默认拆分把另一个替代等位在目标行的基因型中编码为零,所以样本一在 C 行写作一斜线零,在 G 行写作零斜线一。原始基因型却没有任何一条拷贝是 A。若只看其中一行,就很容易把零解释成原始参考纯合拷贝。选择将其他替代等位设为缺失后,C 行保留一条已知携带拷贝并留下一个点,G 行也呈现部分缺失。这个点来自我们的转换选择,并非重新测序后失去了观测。进入群体统计前,最好先查看下游如何计算部分缺失基因型的有效等位数。

我们可以在该位置同时保留两列比例。第一列把携带 C 的观测除以所有三类已分配观测,总共七除以二十一,约百分之三十三点三。第二列只比较 A 和 C 两类,用七除以十,是百分之七十。若研究某个扩增子里目标替代碱基的总体负担,其他替代等位所占的十一条不能悄悄从分母消失。若研究设计本来只想在可明确分入 A 或 C 的观测间比较,就可以使用第二列,列名要写清纳入范围。跨样本应用频率阈值之前,先拿至少一个有两种替代的位点做这个手算。

拆分也会让原来一个位置对应多行。与注释表连接时,应使用染色体、坐标、参考和目标替代的组合键,统计位置总数时则要另设去重口径。保存原始多等位文件、版本和转换命令,必要时用原始位点标识把两行追溯到同一个记录。本文的例子是长度相同的单碱基替换,若换成插入缺失,还要在明确的参考版本上核对左对齐和最简表示,然后再比较两份结果是否代表同一个事件。

当多个样本需要按同一套规则拆分、注释与复核时,可用云生信的生物信息分析服务器保留原始队列文件和各步转换结果。

把知识用于当前任务

下一步可以直接开始分析

按研究目标查找在线工具、课程、计算资源或完整分析服务。

按目标找方案仍有疑问,联系客服 ↗