GTDB-Tk 的比对矩阵怎样检查？实跑六种掩码，找出正常结束却少了一半列的原因

拿到一份微生物基因组分类结果以后，很多分析会直接提取物种名称，再合并到丰度表。真正需要追问的还有参考版本、分类依据和参与比较的有效位点。尤其在整理标记蛋白比对、追加样本或构建自定义系统树时，一处掩码格式变化，就可能改变每个样本保留下来的序列信息。

本文安装 GTDB-Tk 2.8.0，实际调用该版本的原版比对裁列方法，比较六种掩码输入。实验采用自建的四条短蛋白比对，运行范围覆盖矩阵裁列与输出核对。我们同时结合官方分类文档，说明这些检查应当放在实际分析的什么位置，以及分类表里哪些字段值得一起保存。

先确定自己在检查哪一段分析

GTDB-Tk 将基因组与参考分类体系联系起来。当前分类流程会先进行核苷酸相似性筛查，满足相应条件的基因组可以得到物种归属，其余基因组继续进入标记识别、蛋白比对和参考树放置等环节。流程会根据证据路径写入不同的分类方法，因此读取结果时需要保留方法列和相关参考信息。官方分类流程

如果一个基因组已经通过相似性筛查获得分类，后面没有对应的树放置记录，可以先检查分类方法及是否启用了额外放置选项。反过来，一个进入标记比对的基因组，即使序列文件生成成功，也需要检查有效氨基酸比例。把“文件存在”“矩阵有列”“样本有足够可比较残基”分成三项检查，能够更快判断问题发生在数据准备还是分类证据阶段。

本文专门检查第二条路径中的矩阵处理。GTDB-Tk 提供掩码裁列功能，可以使用参考包的掩码，也可以指定自己的掩码文件。官方示例中的每个零或一，对应原始比对的一列，一代表保留。掩码索引针对对齐后的坐标，缺口同样占据一个位置。官方裁列命令

用四条短序列让变化看得见

我们构造一个四行四十列的蛋白比对。第一行全部是明确的氨基酸，第二行只有左半段存在残基，第三行只有右半段存在残基，第四行左半段全部为未知氨基酸符号。它们的行长度相同，信息分布明显不同，适合观察掩码位置和缺失值的影响。

为方便手工核对，每二十列采用二十种常见氨基酸按固定次序排列。这个矩阵只用于验证程序对文件内容的处理规则。读者可以逐字数出留下的位点，也可以把实验脚本中的矩阵替换成自己的标记比对，继续执行相同的格式检查。

>complete
ACDEFGHIKLMNPQRSTVWYACDEFGHIKLMNPQRSTVWY
>left_only
ACDEFGHIKLMNPQRSTVWY--------------------
>right_only
--------------------ACDEFGHIKLMNPQRSTVWY
>unknown_left
XXXXXXXXXXXXXXXXXXXXACDEFGHIKLMNPQRSTVWY

本次环境使用 Python 3.12、GTDB-Tk 2.8.0、NumPy 1.26.4 和 DendroPy 4.6.4。脚本调用安装包中的原版方法，完整保存输入、掩码、输出以及异常记录。这个方法也是裁列功能实际使用的实现，便于在固定的微型输入上观察行为，完整基因组分类另外需要匹配的参考数据及相关运行依赖。本次版本源码

from gtdbtk.misc import Misc

Misc().trim_msa(
    "aligned.faa",
    "file",
    "balanced.mask",
    "balanced.faa",
)

相同的列数，也可能保留完全不同的证据

第一份掩码采用一零交替，共保留二十列。完整序列留下二十个明确残基，左侧有序列和右侧有序列的两个样本各留下十个明确残基，未知符号样本也只有十个明确残基。

第二份掩码只保留左半段，也得到二十列。左侧样本的明确残基占比升到百分之百，右侧样本降到零，未知符号样本也降到零。两份输出的文件结构都正常，每行长度都相同，完整样本也都保留二十个残基，但它们对其他样本提供的比较条件已经发生变化。

这里的比例由实验脚本按二十种明确氨基酸重新计算，分母是裁列后的列数。它适合核对本实验的有效信息量。实际分类表中的字段有自己的程序定义，引用时应核对相应版本的计算方式，避免把两种统计口径混写。

这一组对照提示我们，在调整掩码时，单看矩阵总长度很容易遗漏样本层面的变化。可以同时输出每行明确残基数、缺口数、未知符号数及保留比例。对于某一批样本集中丢失特定标记的情况，整体平均占比还可能掩盖组间差异，因此需要按研究分组检查分布。

为什么把掩码折成两行会改变输出

第三次运行只改变掩码的文件排版。原本四十个字符写在一行，我们将它们折成两行，每行二十个字符。运行正常结束，输出却只剩十列，内容与只提供前二十位掩码的第四次运行完全相同。

检查本次安装的实现可以看到，它只读取掩码文件的第一行，再取出其中标记为一的位置。因此，后半行没有参与裁列。对于这个输入，程序的完成状态与研究者预想的完整掩码长度发生了偏差。这个观察来自上述版本的原版方法，换版本后应重新运行微型对照确认。

这类问题常在复制粘贴、人工编辑或通过文字处理软件保存时出现。序列文件可以折行，并不代表掩码也能沿用同样的文件组织方式。实际协作时，最好保存掩码原文件的校验值，同时把它对应的未裁列比对文件一起归档。收到修改后的掩码，先比对字符数与保留位置，再重新计算下游矩阵。

第五次运行把掩码中的一个位置设为数字二，程序仍然完成，只选择字符一对应的位置，最终保留十九列。第六次运行提供四十一位全一掩码，其中最后一位超过序列长度，出现索引越界异常。可见，异常退出能够发现部分格式错误，短掩码和额外字符仍需要由输入检查来发现。

给裁列过程补上可执行的检查

开始裁列前，先要求输入比对中的所有序列等长，序列标识唯一，掩码只包含零和一，掩码长度与比对列数完全相等。下面的代码适用于明确要求单行掩码的这份实验，也会把误折行的文件直接拦下来。

from pathlib import Path
from gtdbtk.biolib_lite.seq_io import read_fasta

rows = read_fasta("aligned.faa", False)
lengths = {len(seq) for seq in rows.values()}
assert rows and len(lengths) == 1
columns = lengths.pop()

mask_lines = Path("balanced.mask").read_text().splitlines()
assert len(mask_lines) == 1
mask = mask_lines[0]
assert set(mask) <= {"0", "1"}
assert len(mask) == columns
assert mask.count("1") > 0

序列标识唯一性应在原始文件读取时检查，因为一些读取函数会把同名序列放入字典，后出现的记录可能覆盖前面的记录。最简单的做法是先扫描所有大于号开头的行，提取实际使用的标识，再比较记录数和去重后的标识数。文章附带实验使用四个明确且不同的标识，便于把注意力放在掩码行为上。

输出端还可以增加三项核对。第一，输出标识集合与输入一致。第二，所有输出序列的长度等于掩码中一的数量。第三，对每条序列按保留位置重新取字符，并与程序输出逐字比较。第三项尤其适合升级软件或更换矩阵处理工具后使用，它能直接发现偏移一列、错误忽略缺口等问题。

还要留意一种更隐蔽的情况。掩码长度虽然正确，但它来自另一份列顺序不同的比对，此时前面的长度检查全部能够通过。解决办法是在保存掩码时同时记录输入比对的校验值、标记拼接顺序和参考版本，使掩码与生成它的矩阵成对管理。

回到分类表，要把名称与证据一起读取

完整流程产生的分类汇总表包含最终分类、分类方法、参考基因组、相似性及比对比例等信息。对于物种级归属，相似性和比对覆盖范围需要结合相应代表基因组的物种边界解释。只提取最终分类字符串，会丢失判断归属来源所需的上下文。官方汇总表说明

实践中可以把样本编号作为主键，保留最终分类、分类方法、最接近参考及警告列，并将实际用于该分类路径的相似性字段一起导出。不同分类路径可能填写不同参考列，整理时要检查空值含义。来自筛查阶段的结果和来自参考树放置的结果，可以各自追溯到对应的原始文件。

参考版本变化也值得认真对待。本文检查资料时，官方更新说明已经记录了分类工作流在特定物种边界判断及结果元数据方面的调整。复现旧研究时，应一并固定软件版本和参考数据版本，再比较具体哪些基因组的证据字段发生变化。官方更新记录

如果不同批次使用了不同参考包，可以先建立一张批次与版本对应表，定位需要重新分类的样本。盲目把名称替换成统一格式，会让原本不同的分类证据看起来像同一套结果。尤其在比较物种丰度时，同名合并、名称拆分以及代表基因组调整，都需要回到原始分类依据核对。

整理输出时，还应给“输入已有缺口”和“裁列丢弃位点”分别计数。某个样本原来就只有少量标记，即使保留所有列，有效残基也可能很少。另一个样本原始标记充足，却在特定掩码下丢失大量信息，需要回到掩码位置分布继续检查。把未裁列明确残基数、裁列后明确残基数和最终矩阵列数放在同一张表里，才能判断损失来自哪一步。

如果把多个标记先拼接再裁列，还要保存每个标记的起止列。这样可以把丢失的信息重新归属到具体蛋白，检查是否集中在某几个标记。两个输出矩阵总长度相同，仍可能覆盖不同的蛋白区域。对这种差异，逐标记的保留数量通常比一个总百分比更容易解释。

这次实跑可以怎样用于自己的项目

最容易复用的检查顺序，是先确认分类路径，再检查中间矩阵，最后核对结果表。遇到标记比对有效信息偏低，可以逐个检查基因组完整情况、标记识别结果、未裁列矩阵和掩码。这样能够把“标记缺失”“对齐缺口多”“保留位点不适合这一批样本”分别定位到可查看的文件。

本次实跑还说明，格式异常未必伴随报错。短掩码、折行掩码和非二进制字符，分别产生了可以继续被其他程序读取的输出。对这种情况，运行后的内容核对与运行前的格式检查同样重要。把几个十几字节的小型对照保存在项目目录里，后续环境迁移时重新执行一次，能够迅速确认裁列语义。

在正式分析中，自定义掩码的依据应来自预先确定的比对质量要求。先规定允许缺失的程度、关注的标记集合和比较对象，再观察结果是否满足这些条件，可以减少围绕预期树形反复调整参数带来的偏差。需要发表或共享的矩阵，还应附上掩码、参考版本以及保留位点数，方便别人复查。

本文的运行脚本会生成六份掩码及相应输出，并保存异常记录和按样本计算的残基比例。读者可以从原始四行矩阵开始逐项核对，再把同样的方法应用到自己的标记比对中。下载实跑脚本与本次结果

如果接下来需要整理多批次基因组的运行环境、参考数据和可追溯结果，可了解云生信的生物信息学分析计算服务。