Gubbins 实跑教程｜为什么必须保留全长比对坐标，SNP 矩阵跑通了也可能没有意义

同样包含三百八十个多态位点，两份比对交给 Gubbins，为什么一份得到一个候选区域，另一份一个也没有。因为这类分析不仅关心有多少变异，还关心变异出现在哪里、集中在哪条分支，以及局部密度与背景有什么差别。把不变位点全部删除，可能同时删掉方法赖以判断的空间信息。

本篇实际安装 Gubbins 3.4.3，构造八条序列的可控比对，比较变异集中、变异分散和仅保留 SNP 列三种输入。三组都真实运行，保存区域文件、过滤后比对和停止日志。实验只验证输入语义与结果解释，不把人为密集替换称为已经发生的生物学重组。

Gubbins 判断的是局部替换密度异常

Gubbins 通过迭代的系统发育与祖先状态分析，寻找特定分支上局部替换密度升高的区域，再利用过滤后的信息更新分析。这个机制要求输入保留变异之间的距离关系。官方说明将适用对象限定在具有近期共同祖先、差异相对有限的样本集合，并明确指出不能直接用核心基因拼接比对替代全基因组坐标。项目主页

理解这一点，才能看出它与普通 SNP 计数的区别。两个样本即使总差异数相同，如果一个把差异集中在短片段，另一个均匀散布，局部密度证据也不同。只有位点列表而没有原始间隔，就无法还原同一个统计问题。

另一方面，密集替换并不只有重组这一种来源。局部比对错误、组装问题或其他异常也可能产生类似模式。因此输出更适合被称为推断的候选区域，并结合读段、序列质量和研究背景复核。不能把算法识别出的矩形区间直接当作观察到的重组实验事件。

构造保持总变异数一致的两组对照

本次从固定随机种子生成五万列的核酸序列，并建立八个设计样本。先沿着预设的分组关系安排两百个背景变异位点，再让 S1 和 S2 共同增加一百八十个替换。第一组把新增替换集中到一基坐标第一万零一至第一万二千的区间内，第二组将相同数量的新增替换分散到整个比对。

两组都具有八条等长序列和三百八十个多态列，但局部空间分布不同。新增变异关联的样本也都相同，因此比较重点是位置变化，而不是同时改变样本数量、总变异数和支持样本。固定种子与完整脚本保存在复现包中，读者可以检查每一个人为设置的坐标。

这里使用随机核酸背景，而非真实细菌基因组。它没有基因结构、重复区域和真实进化过程，不能用于证明算法在生物样本中的召回率或误报率。它的价值是把一个方法前提单独拿出来演示，使输入改变和结果改变之间的关系更容易观察。

第三组从第一组提取全部三百八十个多态列，删除其余位置，作为故意设置的错误输入。它在格式上仍然是多序列 FASTA，各条序列也等长，但原本五万列中的间隔已被压缩。我们运行它，是为了检验“程序会不会自动阻止这种误用”，并不把它作为可接受的分析方案。

安装后先检查依赖和可写缓存

本次使用固定版本容器，构树和模型拟合均选择 FastTree，使用两个线程与固定随机种子，最多运行五轮。教程不是不同构树软件的性能比较，不能从这一设置推导其他构树方式的速度或稳定性。

IMAGE=quay.io/biocontainers/gubbins:3.4.3--py310hfc0ef84_1
docker pull "$IMAGE"
mkdir -p numba-cache

docker run --rm --cpus 2 --memory 6g \
  -u "$(id -u):$(id -g)" -w /work \
  -e NUMBA_CACHE_DIR=/work/numba-cache \
  -e NUMBA_NUM_THREADS=2 -e OMP_NUM_THREADS=2 \
  -e OPENBLAS_NUM_THREADS=2 \
  -v "$(pwd):/work" "$IMAGE" \
  run_gubbins.py --prefix result --threads 2 \
  --tree-builder fasttree --model-fitter fasttree \
  --iterations 5 --seed 133 input.aln

第一次以普通用户运行时，程序在导入祖先重建模块的过程中失败，错误指向 Numba 无法为相关函数建立缓存。这个阶段甚至还没有开始读取用于实验的比对，因此反复修改输入文件不会解决问题。为缓存指定一个可写目录后，运行才继续完成。

这类安装经验值得保留，因为容器启动成功不代表所有依赖都能在当前权限下工作。检查版本时若以默认身份运行，而正式计算改为普通用户，二者环境条件也可能不同。复现脚本把缓存目录放到输出目录下，让同一个权限设置覆盖检查与正式计算。

正确全长输入的结果怎样读

集中组得到一个候选区域，坐标为一基闭区间 10001–11987，关联样本为 S1 和 S2，区域内记录一百八十个替换。分散组没有候选区域。两组输入都为五万列且含三百八十个多态位点，所以这一差别不能用总 SNP 数不同解释。

集中组过滤后的多态比对剩两百列，与本次设计的背景位点数量一致。这里说的是过滤后的多态位点矩阵长度，并不是声称整个基因组只剩两百个碱基。若在报告中只写“序列由五万缩短到两百”，就会混淆全长比对与多态位点输出两种文件。

候选区间也不必与设计时允许放置替换的整个两千碱基窗口端点完全一致。设计窗口规定的是新增变异可以出现的位置，而程序输出边界来自实际变异分布与判断规则。本次区间结束于第一万一千九百八十七位，因此应原样报告，不应为了让结果整齐而手工改成第一万二千。

错误的 SNP 输入仍然成功结束

仅保留 SNP 列的第三组退出码同样为零，并输出两个候选区域。其中一个位于压缩后坐标 55–212，关联 S1 和 S2，记录一百五十八个替换。另一个位于 272–278，关联 S7 和 S8，记录四个替换。过滤后保留二百一十八列，已经不同于正确全长输入的两百列。

这些坐标是三百八十列压缩矩阵中的位置，不能直接当成原五万碱基序列上的坐标。即便另外保存了列到原始位置的映射，也不能简单映射回去后就把结果当作正确推断，因为算法计算局部密度时使用的间隔已经改变。坐标解释错误和统计输入错误是同时发生的。

这次运行最有价值的现象不是“错格式会报错”，而是格式有效、前提无效的输入可能照常产生整套结果。程序能验证字符和长度，却不能从一个看起来整齐的矩阵中知道中间曾经删掉四万多列。分析者必须在进入流程前核实比对是怎样构建的。

因此，不应把输出树存在、区域文件不为空、程序没有异常视为方法适用性的充分证据。一次完整的分析审查，需要同时检查输入来源、坐标保留方式、样本范围和输出解释。自动化流程可以检查其中一部分，但仍要把这些前提写入分析记录。

从 Snippy 衔接时应该选哪个文件

如果上游使用 Snippy，关键是区分核心 SNP 比对与核心全长比对。前者只保留共同可比较的多态列，适合某些距离和构树任务，但已经压缩了空间信息。后者保留参考坐标范围及缺失字符，经过适当清理后，才保有本类分析所需的位置关系。

snippy-clean_full_aln core.full.aln > clean.full.aln
run_gubbins.py --prefix lineage --threads 2 clean.full.aln

清理字符不等于删除列。运行后应核对各序列长度一致、标题唯一、字符集合符合要求，并确认全长长度与参考坐标约定相符。若输入来自拼接的多个片段，还要保留片段边界与坐标映射，避免把人为相邻的区域解释为基因组上的真实邻接。

同样，泛基因组软件输出的核心基因串联序列虽然常被用于系统发育分析，也不能因为它叫 core alignment 就直接拿来运行。基因顺序、间隔以及不在核心集合中的区域可能已经发生改变。选择输入要根据方法需要的信息，而不是根据文件名中的“核心”二字。

退出成功以后还要看停止原因

本次三组最终都达到设定的五轮上限后停止，日志没有证明已经按收敛标准结束。因此本文报告的是该明确迭代设置下的输出，不能描述成已经获得稳定收敛的最终系统发育关系。此前三轮试跑用于检查流程，正式复现包保存五轮运行结果。

对真实项目而言，达到上限和满足收敛条件需要区别记录。如果研究结论依赖树拓扑或某些边界，应继续检查迭代变化，判断是否需要更合适的样本分组、模型或迭代设置。不能只把上限随意调大，然后假定所有不稳定性都会消失。

同样，固定随机种子主要帮助复现随机过程，并不保证换了软件构建、依赖或数据顺序后所有内部节点编号都保持一致。审核时更应该核对区域坐标、关联叶节点、位点数量和主要关系，而不是只比较树文件中的内部编号字符串。

如何把候选区域变成值得相信的结果

先检查候选区域的原始证据。若它位于覆盖异常或比对质量明显下降的位置，应优先排查数据问题。若区域只在少数低质量样本中出现，也应检查缺失模式是否影响推断。区域图可以帮助定位问题，但不能代替读段和比对本身。

再检查样本群体是否过于分散。将差异很大的多个谱系放进同一次分析，可能改变背景替换水平及局部异常的解释。官方手册建议在合适的近缘群体范围内使用，并讨论了潜在混杂因素。Gubbins 使用手册

最后检查坐标能否稳定回到参考与注释。区域文件描述的是比对上的区间，叠加基因功能时还需要明确坐标系统。若参考来自多条 contig 或经过重排，不应把区间数字直接对照另一份未经映射的注释文件。很多看起来像生物学解释的问题，其实来自坐标对象不一致。

正式报告可以把输入范围、全长列数、多态位点数、候选区域数量、过滤后位点数量、构树设置和停止原因放在同一张摘要表中。这样读者能同时看到数据基础和分析状态。对于关键区域，再补充关联样本和原始证据，比分别堆放一张树和一张彩色热图更容易审核。

优势、边界与本次实跑带来的经验

Gubbins 的优势是把局部替换密度与系统发育背景结合，并提供区域、分支统计和过滤后比对，使分析者有机会追踪某段序列如何影响后续关系。它的代价是对输入比对、近缘范围和迭代状态有要求，不能把任何 SNP 矩阵都视为等价输入。

本次对照给出的经验很具体。总位点数相同不保证空间证据相同，运行成功不保证方法前提成立，过滤后的多态比对长度也不等于剩余基因组长度。把这三个问题写清楚，就能避免把一个看似顺利的流程解释成过强的生物学结论。

下载本次三组实验脚本与实测结果。包内保存设计位点、原始区域记录、过滤后比对、树和日志末尾。建议先检查集中组与分散组都为五万列，再看错误输入只有三百八十列，最后对照区域与保留位点数。复现时可在此基础上换成自己的近缘样本，但必须重新核实比对来源和数据质量。

当全长比对、样本分组与区域复核规则确定后，如需开展更多样本的批量分析，可了解云生信生物分析与计算服务。