Snippy 实跑教程｜少检出一个 SNP，可能只是覆盖不足，怎样避免共识序列制造假一致

同一条细菌序列，改变测序覆盖后，变异表里的 SNP 数量会不会变化。答案通常容易想到，但真正容易出错的是下一步。分析者把变异表应用到参考上生成共识序列，再把没有变化的位置都当作真实参考等位基因。这样一来，数据不足就可能被包装成序列一致。

本篇实际安装 Snippy 4.6.0，用同一个已知序列真值生成三组读段，分别保留正常覆盖、挖去一个覆盖区间，以及在该区间补回少量读段。我们沿着变异表、比对序列、共识序列和多样本核心比对逐层检查，看看一个位点为什么会变成替代碱基、缺口或未知字符。

先看输入条件，再决定能否使用这类流程

Snippy 面向单倍体变异分析，常用于细菌分离株读段相对参考的比较。它整合读段比对、变异调用及结果整理，让多个样本可以进入一致的分析流程。官方仓库提供输入格式、结果文件和核心比对说明。它的便利之处是减少手工拼接步骤，但不会替分析者决定参考是否合适、样本是否混合，以及缺失证据应该怎样进入下游统计。

参考选择直接决定哪些位置有机会被可靠比较。若参考与样本差异较大，一些区域可能无法正常比对，输出中缺失位置增加。这时只看检出的 SNP 总数，可能得到相反的印象，因为许多本应讨论的位置根本没有通过分析条件。真实项目中，应先阅读覆盖、可比对长度和质量，再解释变异数量。

本次刻意使用同一参考片段及其确定的替换，排除了复杂参考差异，也没有引入测序错误、混合菌株和结构变异。这种简化适合验证文件语义，不代表实际样本也会像示例一样干净。文章后面的观察只对本次条件成立，不用于宣称某一调用阈值在所有项目中最优。

构造三组只有覆盖条件不同的数据

参考来自公开的大肠杆菌 NC_000913.3，取前面三万个碱基。我们在从一开始计数的第五千、第一万五千和第二万五千位分别设置明确替换。三组读段都来自这条改变后的序列，所以三处替换是共同真值，不是三组样本之间的生物学差异。

第一组生成八千对长度为一百五十碱基的双端读段，片段长度在三百五十至五百五十之间，使用固定随机种子，并排除重复的起点与长度组合。第二组移除两端读段与指定区间重叠的片段，使中间目标位置没有覆盖。第三组在第二组基础上补回三个不同片段，让目标位置只有三条读段支持。

这里删除的是测序观察，不是把目标变异从样本基因组中改回参考。这个区别非常关键。只有保持真值不变，我们才能把输出差异归因于覆盖和过滤条件，而不是同时改变输入基因型造成无法解释的混合对照。

深度是运行后从 BAM 中按质量条件读取的结果，不是只根据生成读段数量推算。三组在第五千和第二万五千位的深度均分别为七十五和八十二，因而这两个位置还能用来确认实验没有顺便改变所有区域的观察条件。

安装与调用时把关键阈值写明

本次固定镜像中的 Snippy 版本为 4.6.0。每组使用两个线程，显式指定最低覆盖十和变异支持比例零点九。镜像实际帮助中，支持比例参数的默认值为零，表示自动模式。因此本文的零点九是实验设置，不能写成软件默认值。

IMAGE=quay.io/biocontainers/snippy:4.6.0--hdfd78af_6
docker pull "$IMAGE"
docker run --rm "$IMAGE" snippy --version

docker run --rm --cpus 2 --memory 6g \
  -u "$(id -u):$(id -g)" -w /work -v "$(pwd):/work" "$IMAGE" \
  snippy --ref /work/reference.fna \
  --R1 /work/covered_R1.fq --R2 /work/covered_R2.fq \
  --outdir /work/covered --cpus 2 --ram 4 \
  --mincov 10 --minfrac 0.9

比对质量、碱基质量、变异质量和支持比例约束的是不同问题。覆盖足够但比对不唯一，不能通过增加读段数量自动解决。支持比例很高但只有三条读段，也不等于满足最低覆盖条件。阅读命令时，要把这些条件看作共同作用的一组过滤规则，而不是一个可随意替换的质量分数。

复现实验包会自动构造三组读段、运行三次调用并收集结果。新项目建议选择空输出目录，并记录参考摘要、软件版本和参数。尤其在反复调参时，不应只保存最终那次变异表，否则无法解释是哪一步过滤改变了结论。

第一个发现，变异表缺少记录不等于参考一致

正常覆盖组在三个目标位置都检出了预期替换。无覆盖组和低覆盖组都只保留另外两个位置，中间目标没有进入最终变异表。由于我们事先知道三组真值相同，这个差别说明调用结果不仅反映样本序列，也反映当前读段是否提供足够证据。

对于低覆盖组，三条读段支持的都是替代碱基，但它仍未达到本文设置的最低覆盖十。不能因为支持比例看起来是百分之百，就把这个位点视为通过全部质量条件。比例的分母很小，和几十条独立观察得到同一比例的含义不同。

真实分析遇到一个样本比另一个少几个 SNP，应先检查这些位置的深度、比对质量和过滤原因。把未报告位点统一补成参考，会把不可判断状态转换成一个确定答案。更合理的结果表应允许可靠参考、可靠替代和无法判断三类状态同时存在。

第二个发现，共识文件在缺证据的位置保留了参考碱基

本次检查的关键位置是第一万五千位，参考为 G，三个样本的设计真值都是 T。正常覆盖组在比对序列和共识序列中都为 T。无覆盖组的比对序列为缺口，低覆盖组的比对序列为 N，但二者共识序列都为参考碱基 G。

这张小表展示了为什么“共识”这个名字不能代替文件用途说明。对本次输出而言，共识序列应用了通过过滤的变异，但没有在这个缺证据位置保留相同的缺失标记。若直接比较三条共识序列，就会把已知相同的样本真值误读出差异，或者在其他比较里制造参考一致的假象。

也不要把本次比对序列解释为仅将参考打上覆盖标记。实际文件已经在通过条件的位置带入了替代碱基。第五千和第二万五千位在三组中都是 A，正常覆盖组的第一万五千位是 T。应以所用版本的真实输出及对应实现核对文件语义，而不是凭文件名或旧印象做选择。

如果研究问题是样本之间的可比较差异，应选择能保留位置对应和缺失状态的结果，并明确统计时排除哪些字符。若研究问题需要完整的替代后序列，则还必须说明哪些位置没有可靠观察。没有一种文件可以不加解释地同时满足所有后续目的。

第三个发现，三万个位置最后只剩两列核心 SNP

将三组结果交给核心比对程序后，每条全长比对序列仍有三万列，而核心 SNP 比对每条只有两列。这两列对应三组都可以可靠比较的第五千与第二万五千位。正常覆盖组独有证据支持的第一万五千位，没有进入共同核心 SNP 集合。

snippy-core --ref reference.fna --prefix core covered gap low
snippy-clean_full_aln core.full.aln > clean.full.aln

这一现象解释了一个常见项目问题。增加低质量样本后，核心 SNP 数量可能减少，并不表示前面样本的真实差异突然消失，而是共同可比较集合发生了改变。样本筛选与核心定义共同决定了分析分母，不能把不同样本集合生成的矩阵当成同一尺度的直接比较。

因此报告 SNP 距离时，最好同时说明可比较位点数量或使用的共同核心范围。两个样本在很少的可比较位置上没有差异，和在大部分基因组上都没有差异，证据强度显然不同。单独一个零差异数字，会隐藏这种差别。

本次三组其实拥有完全相同的设计基因型，只是测序观察不同。核心 SNP 比对中还包含参考序列，所以保留了样本共同相对参考的两个替换。不要看到矩阵存在两列，就误写成三个设计样本彼此有两个差异。解释距离时还需要看清比较对象是否包含参考。

清理全长比对时保留坐标

核心全长比对中的替代碱基以小写出现，本次正常覆盖组目标位置为小写 t。清理程序先转成大写，再替换不允许的字符，因此该位置变为 T。原有 N 和缺口符号保持其相应形式，清理后的每条序列仍为三万列。

这里的清理不是把所有缺失位置删除，也不是将所有缺口都变为 N。该版本清理脚本保留 A、C、G、T、N 和缺口字符。运行后同时检查长度和字符集合，是一种成本很低但很有效的核查方式。

若下游方法需要变异在基因组中的空间分布，例如随后分析局部变异聚集，就应保留这种全长坐标。只有两列的核心 SNP 矩阵已经丢失原来的间隔，不能用改文件名的方法恢复。清理字符和删除不变位点是两种完全不同的操作，需要分别记录。

从这次对照形成实际项目的检查顺序

先查看每个样本能可靠覆盖多少参考位置，再看低覆盖、无覆盖及异常状态分布。若某个样本明显偏离其他样本，应检查原始数据质量、参考适配性和比对情况。此时直接扩大样本集合生成核心矩阵，可能让所有样本共同可用的范围都被拉低。

接着针对研究结论依赖的具体位点核对读段证据。尤其是某条分支上的关键替换，或用于区分两个近缘样本的位置，不宜只引用汇总表。将位点深度、替代支持、过滤状态和参考坐标一起保存，后续复核才能回答为什么某个差异被保留或排除。

最后再决定下游输入。做核心 SNP 比较、全长坐标分析、功能影响注释和导出替代后序列，需要的文件各有用途。给每个输出添加一句用途说明，比把所有文件打包给读者更有帮助，也能减少下一位分析者误用共识序列的机会。

Snippy 的优势是读段到核心比对的衔接清楚，常用结果容易追溯。限制在于它仍然受参考和调用条件约束，对混合样本、多拷贝区域及复杂变异不能套用本次理想对照的结论。本文没有通过降低阈值把缺失位置硬补为变异，因为实验目的正是观察证据不足如何被表达。

下载脚本并核对三个关键结果

下载实跑脚本与本次结果。驱动包含参考校验、读段生成、三组调用、深度读取和核心比对检查，结果文件保留了原始变异表、字符统计与核心汇总。读者可以先核对目标位点的六十九、零、三覆盖，再检查 T、缺口、N 的差别，最后确认全长三万列与核心两列的区别。

本次经验可以直接用于审查已有项目。只要某个结论依赖“没有检出”，就继续追问对应位置有没有足够证据。这个问题能帮助区分真实参考一致、数据不足和过滤排除，避免一张简洁的变异表把三种不同情况压成同一个答案。

在明确参考、覆盖检查和核心比对规则后，如需组织更多分离株的批量分析，可了解云生信生物分析与计算服务。