Cooler 实跑指南——稀疏计数、矩阵对称、平衡权重与多分辨率，怎样正确读取 Hi-C 数据

第一次读取 Hi-C 接触矩阵时，常见的困惑是文件里记录的总计数与展开成矩阵后的总和不一致。另一类问题更隐蔽，平衡函数已经运行结束，但所有权重都是缺失值。若只确认文件存在，就把它交给下游绘图或结构识别，可能会把存储规则和质量过滤造成的变化误认为生物学信号。

本文实际安装 Cooler 0.12.1，在十个基因组分箱上完成建库、重复像素处理、区间查询、矩阵平衡和多分辨率生成。数据由脚本编写，没有真实 Hi-C 读段，也没有运行比对、去重、结构域或环检测。我们关注的是一个接触矩阵应当如何被保存和解释，以及怎样用可以手算的小案例检查处理流程。

一、先分清分箱表和像素表

Cooler 用分箱表描述基因组区间，用像素表保存分箱之间的接触数。官方概念说明介绍了这种坐标与矩阵的对应方式。本例两条虚构染色体分别长六千和四千碱基，每一千碱基一个分箱，共十个分箱。最后一个分箱故意没有任何接触，用于检查缺失权重的表现。

前九个分箱之间的接触数按一组预先指定的倍数构造，同时在第三个分箱的对角线上放入七次接触。最终像素表有三十七行。像素表一行表示某一对分箱的汇总，不表示一个独立读段，也不必对应一条未经处理的分子记录。

创建前应核对分箱编号是否越界、分箱顺序是否与染色体顺序一致、同一像素是否重复，以及输入究竟保存完整矩阵还是只保存上三角。若这些信息没有约定清楚，文件可以包含正确的整数，却表达错误的接触关系。实际任务还应保存参考版本、分箱大小和上游过滤来源。

cooler.create_cooler(
    "toy.cool",
    bins,
    pixels,
    ordered=True,
    assembly="toy"
)
clr = cooler.Cooler("toy.cool")
raw = clr.matrix(balance=False)[:]

代码明确读取未平衡矩阵，避免把默认权重应用与原始计数混在一起。真实分析中，建议在变量名或输出说明中区分原始、平衡后和其他归一化结果，而不是将它们都简称为接触矩阵。

二、上三角计数与完整矩阵总和为什么不同

默认的对称上三角存储，只需要保存一对分箱的一个方向。读取完整矩阵时，非对角像素会被映射到两个对称位置，而对角像素只出现一次。因此，对展开矩阵的所有元素直接求和，会重复计算非对角接触。

本次像素计数总和为一百三十一，对角线总和为七，展开矩阵的总和是二百五十五。它恰好等于两倍的一百三十一减去七。这三个数都正确，只是统计对象不同。若把完整矩阵的总和直接报告为有效接触数量，就改变了上游计数口径。

这个关系也能用来发现输入错误。如果原本是完整对称矩阵，却被当作上三角再次镜像，非对角接触可能被重复处理。若本来只存一半，却将存储模式标为普通方阵，另一半可能看起来全是零。排查时应先确认存储模式和三角关系，再检查后续绘图函数是否又额外做了对称化。

本例还尝试在要求上三角的文件中写入一个下三角像素，程序明确拒绝。保留这种输入检查很有价值，因为它能在数据进入分析前暴露约定不一致。为了让文件写成功而关闭检查，应当有明确的格式理由，不能只是绕过报错。

三、重复像素合并不是测序读段去重

我们给分箱零与分箱一的像素追加一条计数为五的记录。原像素计数为二，在非有序输入模式下，默认合并后计数变成七。这里发生的是同一个矩阵单元的数值相加，而不是判断两个读段是否来自同一个原始分子。

cooler.create_cooler(
    "duplicates.cool",
    bins,
    duplicate_pixels,
    ordered=False
)

当我们声明输入已经有序并直接写入同样的重复像素时，程序拒绝写入。改用非有序模式但明确要求重复时报错，也同样拒绝。这些行为都已在本次安装版本中验证。参数声明代表对输入性质的承诺，不只是一个提高速度的开关。

如果上游存在重复测序分子，将它们先聚合到像素表后，再让 Cooler 合并重复像素，重复分子的贡献仍然保留。真实的分子去重通常需要读段位置、方向或相关标识等信息，而这些信息在汇总矩阵中可能已经丢失。因此，报告中应将读段层面的去重与像素层面的归并分别描述。

新版本的输入处理细节可能变化，官方接口说明列出了有序输入、重复处理和三角检查之间的关系。复用旧脚本时，最好保留一个含已知重复像素的小例子，核对当前版本到底是相加、拒绝还是按其他规则处理。

四、请求一个区间，不代表只返回区间中的碱基

我们请求虚构染色体上两千五百到三千五百的区间。这个范围横跨两个一千碱基分箱，因此返回的分箱实际覆盖两千到四千，而不是将分箱内部裁切成精确的一千碱基区域。返回矩阵是二乘二，因为查询以现有分箱为基本单位。

region = "chrA:2500-3500"
selected_bins = clr.bins().fetch(region)
selected_matrix = clr.matrix(balance=False).fetch(region)

对于高分辨率图形，这种向分箱边界扩展可能不明显，却会影响与基因、增强子或变异位点的重叠解释。报告候选区间时，应同时给出原始查询范围和实际选中的分箱边界。细化查询字符串本身不会增加矩阵分辨率，也不会创造新的接触证据。

同样需要区分像素表查询与矩阵查询。对一个分箱查询像素表，本例返回它作为第一端存储的七行，其中包括与更后面分箱的接触。要求两端都在该分箱内的矩阵像素查询，则只返回一个对角像素。读取完整矩阵的这一行，又可以看到它与前面分箱的接触，因为对称关系已被展开。

这三种结果都来自同一个文件，却回答不同问题。想查看一个区域与全基因组的互作，不能只取该区域作为第一端存储的上三角记录。想查看区域内部接触，也不能把第一端位于区域、第二端远在区域外的像素全部算进去。

五、平衡运行结束，不代表权重可以使用

本次首先使用默认设置执行矩阵平衡。函数正常返回，但十个分箱中没有一个得到有限权重，收敛状态也为否。原因与这个极小案例的规模有关，默认的非零接触数量门槛为十，而这里只有九个参与接触的分箱，无法满足相应要求。

随后我们为教学案例显式设置较低的非零数量门槛，关闭相应的低覆盖异常筛选，并忽略主对角线参与拟合。这次九个有接触的分箱得到有限权重，最后一个无接触分箱仍然是缺失。这里放宽设置是为了观察算法行为，不是建议真实低覆盖数据一律关闭质量过滤。

weights, stats = cooler.balance_cooler(
    clr,
    ignore_diags=1,
    min_nnz=1,
    mad_max=0,
    min_count=0,
    tol=1e-10,
    max_iters=1000,
    store=True,
    store_name="weight"
)

平衡后，在有效分箱之间去掉主对角线计算行和，九行都接近一，偏差小于本例的核验容差。这个结果与输入故意设置的分箱倍数结构一致，说明计算按预期完成。但真实 Hi-C 数据还受到染色体组成、可比对性和接触结构等因素影响，不能仅凭一个行和检查就宣称所有偏差都被校正。

应同时查看收敛状态、有效权重数量、过滤参数和被排除分箱的位置分布。某一条染色体大量失去权重，与零星低覆盖区域被排除的含义不同。函数返回值、文件写入成功和统计结果可用，是三个需要分别核验的条件。

六、忽略对角线拟合，不会自动删除查询中的对角线

我们在平衡时明确忽略主对角线，但随后读取平衡矩阵，第三个分箱的对角元素仍然存在，约为零点八七五。它由原计数七乘以对应的两个权重得到。忽略对角线约束的是权重拟合使用哪些数据，而不是要求结果文件永久移除这些像素。

这解释了为什么有人在完成平衡后重新计算完整行和，发现某些行不相等。若拟合时排除了主对角线，而核验时又把它加入，就改变了计算范围。应按照拟合所使用的范围检查行和，同时保留原始像素，以便明确哪些数据参与估计、哪些数据只是在查询时应用权重。

无接触分箱也展示了零与缺失的区别。原始矩阵中，该分箱对应的整行都是零，表示存储中没有观测接触。平衡矩阵中，它的整行都是缺失值，因为没有有效权重。把这些缺失直接填零，会让“无法进行这项归一化”看起来像一个经过归一化后确定的零信号。

实际绘图可以对缺失使用独立颜色，也可以在分析中保留有效分箱掩码。无论选择哪种方式，都应让后续步骤能够区分没有观测、被过滤和有效零值，尤其不要在不同样本的有效区域不同的情况下直接比较所有像素。

七、降低分辨率时，应合并原始计数

我们使用多分辨率生成功能，将一千碱基分箱聚合为两千碱基分箱。输出文件中有对应的两个分辨率组，粗分辨率的像素计数总和仍然是一百三十一。脚本按照原上三角像素逐条映射到新分箱，手工汇总的结果与 Cooler 输出完全一致。

cooler.zoomify_cooler(
    "toy.cool", "toy.mcool",
    resolutions=[1000, 2000],
    chunksize=1000, nproc=1
)
coarse = cooler.Cooler("toy.mcool::resolutions/2000")

一个很容易写出的错误做法，是先将矩阵展开成完整对称形式，再把每个二乘二区块直接求和。非对角细分箱之间的接触，一旦落入同一个粗分箱，就会在新的对角单元中被计算两次。原来用于展示的对称镜像，此时被误当成独立接触参与了聚合。

本次正确粗矩阵的对角线是二、十一、二、三、零，而对完整矩阵直接做区块求和得到四、十五、四、六、零。两者在非对角细接触进入粗对角单元时出现差别。这个对照说明，多分辨率构建应理解存储对称性，不能机械套用普通图像缩小的方法。

本次粗分辨率也没有自动复制细分辨率的平衡权重。分箱集合和接触分布改变之后，直接平均旧权重通常不能替代在新分辨率上重新检查和平衡。每个分辨率应有自己的过滤与收敛记录，不能只因为它们放在同一个文件中，就假定分析性质完全相同。

八、这个工具的优势与需要承担的检查

Cooler 的优势是将稀疏接触数据、基因组坐标和分辨率组织在一致的结构中，便于按区域读取，也便于分块处理大矩阵。它的局限是存储本身不等于质量保证，数据的上游筛选、分子重复、参考一致性和结构解释仍需要额外核验。

规模扩大时，优先避免无必要地展开全基因组稠密矩阵。这个建议来自矩阵存储方式本身，不是本次小实验的性能测量结果。本例没有记录吞吐量或者峰值内存，因此不提供未经实测的加速倍数。真正选择分辨率时，也应考虑有效覆盖和研究问题，而不只是文件能否生成。

九、复现时按计数、查询、权重的顺序检查

下载本组完整脚本与实跑记录，在本篇目录执行下面命令，可以重新构造所有像素与输出文件。脚本固定容器摘要，限制为两核、四 GiB 内存，保存版本、预期报错、权重状态、区间查询和粗分辨率核验结果。

python3 run.py output

先确认上三角总和与完整矩阵总和的关系，再比较三种区域查询，最后检查平衡权重与分辨率聚合。这个顺序能把数据结构问题与统计处理问题分开，也能帮助判断一次结果变化到底来自新的接触证据，还是来自读取、筛选和聚合规则的改变。

当项目需要管理多样本、多分辨率 Hi-C 矩阵及其核验记录时，可通过云生信生物分析服务器服务建立可复现的批量处理环境。