ChromHMM 实跑复盘，状态标签全变了，染色质分段就不稳定吗

用 ChromHMM 比较两次分析时，很容易先看状态编号。一份结果的某段是第四状态，另一份变成第一状态，看起来像完全变了。但编号只是模型内部的名称，没有跨模型固定的生物学身份。本次教学实验恰好得到一个很直观的结果，两次四状态训练的逐格原始标签一致率为零，把标签对应关系对齐以后，一致率却为百分之百。真正需要比较的是标记组合和区间边界，而不是颜色或者编号是否相同。

本文实际安装并运行 ChromHMM 1.28，使用自行构造的三种匿名标记、一条虚构染色体，完成训练、固定模型分段、重命名及计数二值化。所有数据都不是实验样本，本文也没有把任何状态命名为启动子或增强子。小数据的作用是让每一次改变都能追溯原因，不用于评价真实数据的注释准确率。软件定位和当前版本可以在官方项目页核查。

一、先把输入的含义固定下来

本次构造了一千二百个连续格子，每格二百碱基，总长度二十四万碱基。三种标记的主要组合依次为仅甲出现、仅乙出现、仅丙出现、全部未出现。每种组合连续维持五十格，四种组合循环六轮，再对每个二值位置以固定随机种子独立加入百分之三的翻转。这样既有空间连续性，也有少量不符合主模式的观测，不是把答案直接写成完全整齐的分段文件。

构造时保留了原始模式表，但它没有参与模型训练。这个区分很重要，教学者知道生成规律，软件只接收到二值矩阵。检查结果时可以研究模型如何概括这些规律，却不能把这次人为容易区分的模式当作真实组蛋白标记的全部复杂性。真实数据还涉及抗体质量、测序覆盖、重复一致性和可比性，这些问题在本例里都没有被模拟。

二值文件前两行分别标识细胞与染色体、标记名称，之后每行对应一个连续格子。下面用 \t 显示制表符的位置，实际文件使用真正的制表符，由复现脚本生成，不应把这段显示文本直接存成输入。

CellA\tchrToy
MarkA\tMarkB\tMarkC
1\t0\t0
1\t0\t0

在这里，零是经过前面处理后判为未出现，一是判为出现，二才表示缺失。缺失意味着这个位置的观测没有提供相同的信息，不能因为表格需要填满就改成零。我们后面专门用同一个训练模型验证这个差别。文件组织和取值定义可核对官方手册，但理解定义之后，仍应检查自己生成的文件有没有把缺失误写成阴性。

二、两次训练的标签全不同，边界却完全相同

实验在固定摘要的容器中运行，Java 版本为 25.0.2-internal。每个容器限制两核、四 GiB 内存，Java 堆上限一 GiB。这些是本例的环境记录，不是全基因组任务的资源建议。下面展示其中一次训练命令，其他参数相同的第二次运行只把随机种子改成一百六十六。

ChromHMM.sh -Xmx1g -Djava.awt.headless=true LearnModel \
  -b 200 -l genome.txt -init random -s 165 -p 2 -r 100 \
  -noautoopen -nobrowser -noenrich -noimage \
  -printposterior -printstatebyline \
  binary learn4 4 toy

第一轮四状态训练得到二十四段，第二轮也是二十四段。程序按每个格子的状态检查，发现对应位置的数字标签一个都不相同。但把第二轮的一、二、三、四状态分别映射到第一轮的四、三、二、一以后，所有格子的分段状态都一致。此时如果直接统计数字相等率，就会错误地把一次标签置换描述为模型不稳定。

实际项目不能只找一种能提高一致率的配对，然后宣布模型等价。至少要同时检查发射参数是否代表相似的标记组合、各状态覆盖多少区域、区间边界是否一致。若不同状态数量的模型出现一对多关系，还应该保留这种拆分信息。两个精细状态同时对应一个粗状态，并不意味着必须强行合并成一个相同名称。

我们还实际执行了重排序命令，将第一轮的编号反向对应，并逐行核对输出。二十四个区间的染色体、起点和终点全部不变，变化的是标签名称，输出前缀也从默认的字母变为重新排序后的形式。这个检查比只看重画后的热图更可靠，因为热图的颜色和排列都可以变，坐标是否改变却是另一个问题。

这一结果给比较工作一个很实用的次序，先确认输入矩阵和格子宽度一致，再比较发射模式，最后才比较重命名后的分段。如果直接用状态编号连接两份表，很可能在没有报错的情况下得到一份逻辑错误的变化矩阵。格式正确并不能保证比较对象正确。

三、看发射参数时，先读列名再读颜色

第一轮四状态的发射表有一个容易忽略的细节。输入列顺序是甲、乙、丙，输出展示顺序却变成甲、丙、乙。如果照着输入顺序给热图列贴标签，原本代表乙标记的状态就会被解释成丙标记。脚本因此按表头读取数值，而不依赖第三列在视觉上处于哪个位置。

本次第一状态的乙标记发射概率约为零点九七七，第二状态的丙标记约为零点九六四，第四状态的甲标记约为零点九五四。第三状态三种标记均低，约为零点零二至零点零五。这些数值支持对本例状态进行“甲高”“乙高”“丙高”“全部低”这样的描述，但并没有提供启动子、增强子或抑制区的生物学证据。

发射概率也不是一条区域属于某个状态的概率。前者描述在给定状态下观察到某种标记的倾向，后者是结合观测与模型后对位置状态的推断。把两者混在一起，会出现一种常见误读，某状态的甲标记概率很高，于是所有被分到该状态的格子都被认为必定观测到了甲。我们故意加入的翻转已经说明，单个格子可以偏离状态的主要组合。

对于真实分析，状态名称应当是后续解释的结果。先依据标记组合提出有限的描述，再借助外部注释检查与转录起始点、基因体或其他区域的关系，并记录支持和不支持的证据。名称越具体，需要的证据越多。只凭一幅色彩清晰的发射热图，就给每个状态分配确定功能，会把统计概括误写成直接测量。

四、全缺失和全零，为什么都能得到标签

接下来固定第一轮四状态模型，不再训练，分别对两份新的矩阵分段。一份所有位置都写二，另一份全部写零。它们的大小和坐标完全相同，唯一差别就是缺失与未出现。

全零输入得到一个覆盖整条虚构染色体的第三状态区间。全缺失输入仍然得到四段，而且也有大量位置被分为第三状态。如果只看最终标签，可能会觉得两者没有多大区别。但检查后验分布，差异就非常明显。全零输入所有格子的最大后验概率平均约为零点九九九九七，全缺失约为零点三零五。

以从零开始编号的第六百格，也就是自然计数的第六百零一格为例，全缺失的四个后验值约为零点二三九、零点二三六、零点二八八和零点二三七，分布相当分散。全零输入在同一位置的输出则几乎全部落在第三状态。输出保留了有限小数，显示为一不能解释成脱离数值精度的绝对确定。

这个对照说明，软件能输出一个标签，并不表示输入提供了充分证据。在观测全部缺失时，初始状态和状态转移结构仍然会参与推断，所以模型并不会因为没有标记观测就自动输出空白。此时标签来自模型已有约束，与看到三种标记均未出现后得到相同标签，证据来源不同。

本次两种输入的硬标签有一百四十二格不同，其他格子即使标签相同，也不能说信息等价。实际审阅时应把缺失掩码与后验分布一起保存，尤其不能用填零操作把低覆盖、未测量和阴性结果混为一谈。若一个区域主要依靠邻近状态延续得到标签，报告中应当能够识别这种情况，而不是让一张完整染色图掩盖它。

五、后验更接近一，不代表状态数量选得更好

本次另外训练了一个两状态模型，使用同一份输入和相同的初始化种子。四状态模型可以分别概括三种单标记模式与全部低的模式，两状态模型则只能压缩这些区别。输出结果里，一种状态主要保留丙高的特征，另一种把甲高、乙高和全部低等情况混在了一起。

两状态模型得到十三段，四状态模型得到二十四段。但两状态模型每格最大后验概率的平均值约为零点九九七五五，反而略高于四状态的零点九九七二五。这个差别不说明两状态更接近生物学真相，它只说明模型在自己有限的类别之间作出了很集中的分配。选项更粗，也可能更容易作出确定选择。

因此，选择状态数不能只追求最清晰的热图、最高的后验值或者最少的区间数。应先说明研究需要区分到什么程度，再检查增加状态是否形成可解释、可重复的细分。一个状态被拆成两个但标记组合几乎相同，和原本混在一起的两种组合被分开，是不同的结果。外部注释可以帮助评价，但也要避免一边反复调整模型，一边把同一份注释当成独立验证。

本例提供的是概念上的反例，没有进行全范围状态数搜索，也没有建立真实数据上的最优模型。将模拟中的四状态原样推荐给其他研究，恰好违背了这个实验想说明的事情。模型复杂度需要跟输入信息和分析目标对应，不能跟一篇教程的演示参数绑定。

六、二值化是上游决策，不能当作无损转换

为独立检查二值化，我们另建了一百格、两种标记的计数输入。第一种只在其中二十格写入二十个计数，其余为零。第二种在全部格子写负一，表示该计数输入中的缺失。注意这里计数文件的缺失写法与训练用二值文件不同，不能互换文件格式。

默认运行后，第一种标记有二十格被判为出现，第二种的一百格都转为二值缺失。再增加一个三十的绝对计数门槛后，第一种没有任何格子被判为出现，而第二种仍然全部保持缺失。这两次运行使用的原始计数完全相同，变化来自判定条件，而不是样本生物学发生变化。

ChromHMM.sh BinarizeSignal signal binarize_default
ChromHMM.sh BinarizeSignal -g 30 signal binarize_high

这也提醒我们，不应把任意归一化后的连续轨道直接当成原始计数输入。已经进行对数转换、比例缩放或者其他校正的数值，与一个格子里观察到多少条片段不是同一对象。即使软件接受了数字文件，后续阈值的统计含义也可能已经改变。分析记录应该保留输入单位、二值化方法、控制数据和阈值，而不只是最后的零一矩阵。

本次安装与运行还留下两个复现经验。第一，设置了随机种子后，模型文件名含种子标识，不能在后续命令中盲写一个不存在的固定文件名，应读取实际产物并确认唯一性。第二，虽然传入两线程上限，本例只有一个输入序列，训练日志实际报告使用一个线程。线程参数不是已经获得相应加速的证据，这个小实验也没有测量并行性能。

七、拿到结果以后怎样逐层检查

建议先检查三件可以直接计数的事情，二值矩阵行数是否对应染色体长度、分段是否落在预期格子边界、每个状态累计长度是否与总覆盖长度相符。接着检查标记表头与发射参数，再检查缺失位置及后验分布。只有这些基础关系一致，才适合继续做功能命名和条件之间的比较。

ChromHMM 的优点是能够把多种标记的组合与空间连续性放在同一模型中，并保留可以审阅的参数与分段。限制也很明确，二值化会压缩强度信息，状态编号不能跨模型直接比较，功能标签需要额外证据。它适合用来提出和组织染色质状态的假设，不会替代原始实验质量控制，也不会把缺失数据变成真实观测。

下载本组脚本与本次结果记录，解压后在本篇目录执行下面命令即可重建全部模拟输入。脚本固定镜像摘要，记录命令、退出码、模型、后验和比较结果，既能核对本文数字，也能修改翻转比例或缺失位置继续练习。建议每次使用新的输出目录，避免把旧模型误认为当前结果。

python3 run.py output

当流程扩展到多细胞类型、多套标记和多轮模型比较时，可通过云生信生物分析服务器服务组织环境、数据与可追溯的分析任务。