扩增子或混池测序中,常见的提问是“两条替代读段能不能算一个变异”。先看分母。两条来自四条总深度,观察比例是二分之一,两条来自四十条总深度,就只占百分之五。我们在同一份 BAM 里放入这两种情况,又添加一个四十条中有八条替代的位点,用 FreeBayes 的连续混池模式逐次调整候选门槛。四组输出不仅展示哪些位置被留下,还留下一个很好的反例,有一行明明写出了替代碱基,基因型字段却是参考纯合。
候选发现要回答的两个问题
FreeBayes 从读段比对中寻找替代等位的证据。最低替代观测数约束的是“至少要看到几条”,最低替代比例约束的是“相对于该样本参与判断的观测,它占多少”。这两个参数可以同时起作用。只限制条数时,深度很高的位点即使只有很小一部分读段出现替代,也有可能进入候选。只限制比例时,低深度位置的少数几条读段就可能使比例和预期相当。建模前先把两项各自的分析目的写清楚,才能解释为什么改变一个参数会专门影响某类位点。FreeBayes 作者对候选与混池分析的说明
本例主动选用连续混池模式,以读段支持和观察比例为主要教学对象。在这样的模式下,输出某行不等于已经证明存在一个高质量二倍体杂合。样本里有几个个体、每个个体贡献是否相同、读段是否来自独立分子,都会影响观察比例与生物学比例的关系。如果要把结果解释为个体基因型,应先选择与实际倍性和实验单位相匹配的模型。下面的三个位置都来自已知构造,便于把候选是否出现与最终结论分开检验。FreeBayes 官方版本与可执行程序
实跑经验|用三个位点把两道门槛拆开看
我们安装 FreeBayes 1.3.10,人工生成一条一千九百碱基的参考。每条读段长八十个碱基,统一定位质量和碱基质量,在每个目标附近错开读段起点,正反链均有观测。只在这三个互不相邻的位置插入指定数量的单碱基差异。输入先经排序和建索引,然后固定参考、样本与单倍型长度设置,在四次运行间只改变最低条数或最低比例。下表的“入选”是指在实际输出 VCF 里出现这一行。
| 最低替代条数与比例 | 301 位,2/4 | 701 位,2/40 | 1101 位,8/40 |
|---|---|---|---|
| 2 条,20% | 入选 | 无记录 | 入选 |
| 2 条,5% | 入选 | 入选 | 入选 |
| 3 条,5% | 无记录 | 无记录 | 入选 |
| 2 条,21% | 入选 | 无记录 | 无记录 |

当最低条数保持两条,最低比例从二成降到百分之五,新增的是第七百零一位。它的四十条读段中有三十八条参考、两条替代,刚好达到较低比例。第一位原本只有四条读段,其中两条替代,早已超过较高比例。反过来,固定百分之五的比例要求,把最低条数从两条提到三条,前两处都消失,即使第一处的观察比例高达一半。最后把比例从百分之二十提到百分之二十一,四十条中八条替代的位置退出,展示了“刚到边界”和“稍高于边界”的实际区别。
第七百零一位尤其值得打开整行看。放宽到两条、百分之五以后,替代碱基 A 被写出,样本列仍给出零斜线零,位点质量约为二乘十的负十四次方。两条替代观测并没有因为出现在 VCF 中就自动得到高置信分数。第一个位置在相同参数下输出杂合样本字段,位点质量约二十七点六。第三处八条替代的位点质量约八十八点八。它们让读者在一份文件里同时看到候选门槛、观测数量和模型评分回答的是不同问题。把所有有输出行的位置都当作已确认阳性,会抹去这几项差异。
freebayes -f ref.fa --pooled-continuous --haplotype-length 0 \
--min-alternate-count 2 --min-alternate-fraction 0.05 input.bam > loose.vcf
bcftools query -f '%POS\t%REF\t%ALT\t%QUAL[\t%GT:%RO:%AO:%DP]\n' loose.vcf
拿到一个低频候选,下一步查哪些证据
先回到原始比对,核对替代观测是否来自同一扩增家族、同一起点、同一方向,是否都落在读段末端或插入缺失附近。两条支持若来自同一模板的重叠双端,独立证据量与两份独立模板不同。若实验带有分子标签,可在与建库方法相符的纠错或共识步骤后再核对。若同一替代在多个阴性对照同一位置反复出现,还需要查看序列背景、引物与批次。单纯降低频率阈值无法区分一个稳定的实验偏差和偶然观察到的真实低频等位。
复核时不要先把所有样本的读段混到一行。一个样本中有两条替代和另一个样本中有两条替代,是否满足检测流程的条件,要按相应选项和样本划分来解释。若文件头把多个样本误写为相同样本,候选的分子、分母和后续模型输入都会发生变化。另一个实用核算是单独列出深度不足的位点。四条里的两条虽然比例高,采样不确定性也大,报告总深度与替代条数能让复核者一眼区分四分之二和四十分之二十。
要选适合本实验的门槛,可以用已知梯度的混合材料或独立技术重复,对每个稀释点统计能否被找回,同时查看空白和阴性对照的背景。把候选发现、质量复核、最终报告分成可追踪的几列,记录哪一步丢失了感兴趣的位点。这个小例子没有模拟测序平台的全部错配来源,也不能据此给真实实验指定最低检测限,三个已知比例提供的是一套可以自己验证阈值行为的方法。
当重复样本和稀释梯度增加,需要用一致参数批量比较候选列表时,可使用云生信的生物信息分析服务器保留各轮变异结果与证据统计。
