modkit 读表指南——修饰概率、有效覆盖与单倍型汇总，怎样避免误读甲基化比例

纳米孔修饰碱基分析的结果里，一个百分比往往经历了多层选择。先由模型给出逐读段概率，再将概率转换成类别，排除不够可信的调用，最后按位点、修饰类型或单倍型汇总。直接把最终百分比理解为“所有读段中有多少被修饰”，有时会用错分母。

本文采用官方文档解读与原创算术练习，没有安装或运行 modkit，也没有处理真实 BAM 或原始电信号。我们实际运行自编 Python 脚本，核验概率门槛、bedMethyl 计数、位点聚合和单倍型分组中的数字关系。练习用于理解结果含义，不是 modkit 的完整实现，也不能用于评价修饰识别准确率。

一、先确认自己拿到的是哪一层结果

modkit 处理带修饰信息的比对结果，可以把逐读段修饰信息整理为位点汇总或其他分析结果。它不能从一张没有修饰标签的普通比对表中凭空恢复修饰证据。审阅项目时，应先说明上游识别模型、模型适用的修饰类型、参考版本和数据处理方式，再讨论汇总参数。

修饰标签检查文档介绍了检查无效标签与相关错误的用途。实际工作中，应关注被忽略记录的数量和原因，而不只是输出文件是否生成。若一批样本丢失的记录比例明显高于另一批，最后的位点覆盖和百分比就可能受到不对称影响。

SAM 的修饰标签有自己的编码规则，格式规范区分了修饰位置、概率以及生成这些信息时的序列长度。不要把修饰位置编码直接当作参考基因组坐标，也不要把标签缺失简单理解为该位置未修饰。序列裁切、比对方向和未明确列出位置的解释，都需要与标签约定一致。

在交接数据时，最有用的说明并不是“已经生成甲基化 BAM”，而是具体保存了哪些修饰类型、采用哪套模型、是否做过影响序列或标签的处理，以及有多少记录通过检查。这样的记录能让后续汇总建立在明确的输入边界上。

二、概率高低与最终类别不是一回事

为便于理解，我们构造八条假想调用，每条包含普通胞嘧啶、甲基胞嘧啶和羟甲基胞嘧啶三种概率，三者之和为一。所有概率都是人工指定，不来自任何真实模型，也没有进行标签概率解码。练习只使用三类中的最大概率作为候选类别，再按统一门槛判断是否保留。

其中一条调用的三类概率为零点零三、零点九二和零点零五，因此最可能的类别是甲基胞嘧啶。另一条是零点二、零点六八和零点一二，虽然甲基胞嘧啶仍是最大值，却无法通过零点七的门槛。被排除的调用应标记为未通过，而不是转成普通胞嘧啶。

这个例子说明，最大概率分类与置信门槛分别回答“哪种类别最可能”和“证据是否足够强”。概率最大的类别不一定值得保留，未保留也不等于已经证明另一类别正确。将未通过门槛的调用补成未修饰，会系统性改变分母与比例。

官方门槛示例也将类别选择和通过条件分开说明。本文的脚本没有处理类别并列、特定修饰门槛、编码精度或其他实际实现细节，因此只能用来验证这里明确写出的简单规则，不能替代软件对真实数据的判断。

三、调高门槛，修饰比例不一定单调变化

对同样八条人工调用，门槛设为零点七时保留七条，其中三条为甲基胞嘧啶、两条为羟甲基胞嘧啶、两条为普通胞嘧啶。门槛提高到零点八后保留五条，再提高到零点九后只剩三条。通过数量持续减少，但最终比例并不只有一种变化方向。

在这个练习中，甲基胞嘧啶比例逐步降低，两类修饰合计比例却先降后升。原因是不同类别被移除的数量不同。门槛变化既改变分子，也改变分母，不能仅凭“更严格”推断结果一定更低或者更高。

八条调用的甲基胞嘧啶概率平均值是零点四三一二五，也不等于某个门槛下按类别计数得到的比例。概率平均保留了不确定性数值，硬分类比例统计的是通过门槛后的类别数量。两者各有用途，但不能在比较样本时一边使用前者、一边使用后者，再把差值解释为生物学变化。

真实项目应报告实际生效的门槛及其确定方式。过滤说明指出门槛可以由数据估计，也可以明确指定。因此，两个样本命令行看起来相同，并不足以证明它们最终使用相同的数值门槛。需要结合日志与运行设置核对，尤其是在跨批次比较时。

四、bedMethyl 的有效覆盖要按定义计算

下面进入独立的计数表练习。我们设一个位置有四次甲基胞嘧啶调用、两次羟甲基胞嘧啶调用和三次普通胞嘧啶调用，另外有三次未通过门槛、一次缺失、两次其他碱基和四次没有修饰调用的记录。它们只是用于读表的人工分类计数，不是上面八条概率的汇总。

按照bedMethyl 字段说明，这一位置的有效覆盖为九。甲基胞嘧啶所在行的目标修饰计数为四，其他修饰计数为二，普通碱基计数为三，因此比例是四除以九，约为百分之四十四点四四。

如果把所有类别加起来得到十九，并使用四除以十九，就会变成约百分之二十一点零五。这个分母混入了未通过、缺失、其他碱基及无调用记录，回答的已经不是该字段所定义的有效调用比例。文件中存在这些辅助计数，正是为了让读者检查证据如何损失，而不是让它们自动进入所有百分比的分母。

valid = modified + canonical + other_modified
percent_modified = 100 * modified / valid if valid else None

这个公式只处理我们选定的普通 bedMethyl 语义。不同输出规范或版本可能调整字段解释，不能仅凭某个文件也有十八列就直接套用同一解析器。读表前应固定输出类型，再验证少量记录中的计数恒等式，最后才批量计算。

五、两种修饰各一行，不代表两批独立读段

同一个位置的羟甲基胞嘧啶行，目标计数为二、其他修饰计数为四、普通碱基计数仍为三，有效覆盖仍是九。两行共享同一批有效调用，只是把不同修饰设为当前目标。将两行有效覆盖直接相加成十八，会重复计算这批证据。

如果研究问题关注两类修饰合计，人工表的合计比例是六除以九，约为百分之六十六点六七，而不是六除以十八。是否适合合并两类修饰取决于研究问题，合并后将无法从单个比例中恢复各自贡献。因此，原始分类计数值得保留，避免下游需要区分时已经失去信息。

同样的位置若还有不同链、不同序列模式或不同分组，也不能只按染色体和起点直接去重。首先要确定每一行的统计单位，哪些行共享读段、哪些行表示不同子集。否则，删除所谓重复行可能丢失修饰类型，直接相加又可能重复计数。

缺乏有效覆盖的位置也不能填成零比例。在练习脚本中，零修饰除以零有效覆盖返回空值。零修饰但有十次普通碱基调用则可以得到零比例。一个表示没有足够数据定义比例，另一个表示有效观测中未见目标修饰，应该在表格和图形中保留这种区别。

六、区域平均会因统计单位不同而改变方向

我们再构造两个样本、两个位点的例子。样本甲在第一个位点是九次修饰除以十次有效覆盖，在第二个位点是一除以一百。样本乙分别是四除以十和四十除以一百。两个样本的覆盖结构相同，只有修饰计数不同。

若先计算每个位点比例，再对位点等权平均，样本甲是百分之四十五点五，样本乙是百分之四十，甲高于乙。若把区域内的计数先合并，样本甲是十除以一百一十，约百分之九点零九，样本乙仍是百分之四十，甲反而低于乙。

这不是算术错误。前者让每个位点权重相同，后者让覆盖更高的位点影响更大。研究问题究竟关注区域中一个典型位点的平均状态，还是所有有效读段调用的总体构成，需要事先说明。不能在看到两种结果后，挑选更符合预期的那一个作为唯一报告。

真实区域中的覆盖差异还可能与序列特征、可比对性和模型适用范围有关。除了报告一个区域比例，保留位点数、覆盖分布、位点间差异和缺失比例，会让读者更容易判断汇总是否被少数位置主导。遇到方向翻转时，先拆回位点层面查看，往往比继续调整统计阈值更有效。

七、单倍型分组后，合计不一定等于两个分组相加

单倍型相关分析需要区分已分组读段与未分组读段。我们构造第一单倍型四次有效调用，其中三次修饰，第二单倍型也有四次，其中一次修饰，另外还有两次未分组的普通碱基调用。两个单倍型的修饰比例分别为百分之七十五和百分之二十五。

将所有记录合在一起，有效覆盖是十、修饰计数是四，总体比例为百分之四十。若只把两个单倍型相加，则覆盖只有八，比例为百分之五十。这个差别来自未分组记录是否纳入，而不是比例计算本身出了问题。

审阅真实输出时，要检查合并结果是否包含未分组记录，并核对单倍型标签来自什么分相流程。一个单倍型覆盖明显低于另一个时，估计的不确定性也不同。看到两个百分比不同，不能立即宣称存在可靠的等位基因特异修饰，更不能仅凭标签编号指定父源或母源。

链合并与单倍型合并也不是同一操作。链描述参考或修饰调用的方向关系，单倍型描述读段所属的遗传分组。把两个概念混用，会让样本注释看起来完整，却无法解释究竟合并了哪些记录。整理分析表时，建议分别保留链、单倍型、修饰类型和位点坐标。

八、怎样把读表检查放进实际流程

一个可审阅的流程应当先检查输入标签与模型信息，再确认目标修饰类型和筛选门槛，随后验证计数表中的有效覆盖关系。之后才进行位点筛选、区域汇总和组间比较。这样能把输入错误、质量选择和统计问题分开，而不是把所有异常都归到“样本差异较大”。

质量检查不必一开始就覆盖所有图形。可以先选取覆盖充分、覆盖不足和比例极端的几个位置，逐个核对分母与类别。再检查有无整批样本共享的异常，例如未通过比例偏高、未分组比例偏高或某种修饰类型完全缺失。少量能解释的核验，比大量没有明确判断依据的图更有帮助。

modkit 的实用价值在于连接逐读段修饰信息与多种汇总视角，适合在保留细节的同时形成可读结果。使用代价则是必须理解标签、分类与汇总之间的关系。上游识别模型的适用性、样本重复和实验混杂，仍然需要研究者独立判断，位点表本身不会替代这些证据。

九、复现本文练习与使用范围

下载本组原创练习脚本与核验记录，进入本篇目录运行下面命令即可复算八条概率、两行人工 bedMethyl、区域平均和单倍型例子。脚本只依赖 Python 标准库，不调用 modkit，不读取真实修饰标签，也不执行碱基识别。

python3 run.py output

本篇选择读表练习，是因为项目公开许可对使用范围有明确限制，我们没有把本次内容制作写成获授权的软件实测。读者可以独立复现这里的算术关系，但这些结果不能证明某个安装版本在所有输入条件下的行为，更不能用来替代真实模型的验证。

复算后最值得保留的经验，是每次看到百分比都同时询问统计单位、有效分母和筛选条件。把这三项说明清楚，很多看似矛盾的结果就能被定位到具体环节，也能让跨样本和跨批次比较更有依据。

当需要统一整理多批修饰结果、运行原创审阅脚本并保存复现记录时，可通过云生信生物分析服务器服务组织分析环境。