methylKit 实跑指南——覆盖筛选、加权差异与合并重复，为什么会改变结果

甲基化差异分析中，一个位点为什么消失，往往比它最后的显著性更值得先问。它可能在读取时就被跳过，可能在某个样本中过滤掉，也可能因为要求所有样本共同覆盖而无法进入检验。若只比较最终差异位点的数量，这些性质完全不同的变化很容易混在一起。

本文实际安装 methylKit 1.36.0，在 R 4.5.2 中完成计数表读取、覆盖筛选、位点合并、覆盖标准化、差异检验和磁盘数据库操作。数据由脚本生成，包括两组三个重复、五百个候选 CpG，以及可追溯的缺失和差异设计。这里没有真实测序样本，也没有运行比对或甲基化提取，目的在于核验分析规则，不是比较工具准确率。

一、先让计数表保留实验信息

methylKit 可以接收多种甲基化结果格式，官方教程介绍了从计数导入到差异分析的主要步骤。本例选择 Bismark coverage 格式，每行保留染色体、位置、百分比、甲基化计数和非甲基化计数。这里的文件由我们编写，并不意味着本次运行了 Bismark。

六个样本依次为三个对照和三个处理，分组向量与文件顺序一一对应。普通位点的覆盖数按样本设置不同倍数，同时加入样本间比例波动，再按覆盖数生成整数计数。这样，同一组内不会被人为设成完全一致，后续模型也有实际需要处理的重复间变化。

我们另外设置了几类可核对的异常。十个位点在全部处理样本中缺失，二十个位点在第三个对照和第三个处理样本中缺失。另有十个位点在第二个对照中只有三次覆盖，以及一个位点在第一个对照中达到十万次覆盖。它们的位置互不重叠，因此每一步减少多少行都可以独立算清。

输入检查应当同时覆盖非负整数计数、位置唯一性和分组对应关系。样本列错位可能不会让模型报错，却会改变研究问题。百分比还应与两种计数之和匹配，不能一边使用某套过滤后的百分比，一边搭配另一套未过滤计数。保留原始提取文件与导入参数，通常比事后追查一张最终结果表更省力。

二、这次真正遇到的问题是首行被跳过

最初读取无表头文件时，我们没有显式指定是否存在表头。程序能够完成读取，但六个文件各少了一行，首个坐标在所有样本中消失。这个错误没有表现为无法解析文件，反而产生了看起来正常的对象，随后覆盖数核对才暴露出不一致。

显式设定没有表头后，六个文件分别读入五百、五百、四百八十、四百九十、四百九十和四百七十行。我们又将默认读取与明确无表头读取作为两次独立实验，确认差异正好是每个文件的第一行。这个细节说明，文件格式名称不能代替所有导入参数。

raw <- methRead(
  as.list(files),
  sample.id = as.list(ids),
  assembly = "toy",
  treatment = c(0, 0, 0, 1, 1, 1),
  context = "CpG",
  pipeline = "bismarkCoverage",
  header = FALSE,
  mincov = 1
)

上面的低读取门槛只是为了保留本例故意设置的低覆盖记录，便于观察下一步筛选。真实任务中，读取阶段已经剔除的位点不会在后续放宽参数时自动回来。因此，报告覆盖过滤规则时，既要记录导入门槛，也要记录之后的显式筛选，不能只写最后一个参数。

三、逐样本过滤与共同位点筛选要分开统计

本例使用最低覆盖十、最高覆盖五百进行筛选。第一个对照因极端高覆盖减少一行，第二个对照因低覆盖减少十行，其余样本不变。这一步是在各样本内部完成的，还没有要求某个位点在其他样本中存在。

filt <- filterByCoverage(raw, lo.count = 10, hi.count = 500)
strict <- unite(filt, destrand = FALSE)
relaxed <- unite(filt, destrand = FALSE, min.per.group = 2L)

默认要求所有六个样本共同覆盖时，最终留下四百五十九个位点。计算过程是从五百个候选中，去掉整组缺失的十个、两个样本缺失的二十个、一个样本低覆盖的十个，以及一个样本极高覆盖的一个。这不是五百减去所有样本损失行数的简单总和，因为同一位置可能在多个样本中缺失。

每组至少两个样本覆盖时，只有全部处理样本缺失的十个位点无法进入，最终留下四百九十个位置。其他位置中缺失的样本仍然是缺失值，并没有补成零甲基化。这个区别直接影响解释，因为没有读段证据与观察到全部未甲基化的读段，并不是同一种数据。

放宽共同覆盖要求可以保留更多位置，但不同位点使用的有效样本数也会不同。真实项目应检查缺失是否集中于某组、某批次或特定区域。若缺失本身与处理有关，增加保留数量并不能消除选择偏差。本文后续主要检验使用严格集合，避免把不同位点样本数变化与其他实验同时混合。

四、覆盖标准化不是让甲基化比例变得相同

我们对过滤后的六个样本独立运行中位数覆盖标准化。处理前的中位覆盖依次为三十、五十八、一百二十、三十、六十一和一百一十六，处理后均变为一百二十。这个操作改变的是计数的尺度，并不要求每个样本的甲基化比例分布一致。

第一个对照的示例位点原来是十次覆盖，其中一次甲基化。标准化后变成四十次覆盖、四次甲基化，比例仍为百分之十。第二个对照的同一位点原来是百分之九十，经过非整数倍缩放和取整后，覆盖为二百零七、甲基化计数为一百八十六，比例出现小幅变化。

本次最大比例变化不到一个百分点，来自计数取整。还有个别位置的总覆盖与两类计数分别取整后的和相差一。这提醒我们，标准化后的整数是用于分析的变换结果，不应重新包装成原始测序分子数，也不宜与未经变换的计数混合后再送入另一套流程。

更关键的是，标准化不能增加实验观测，也不能制造生物学重复。本文将它作为独立演示，后面的差异检验仍使用未做此缩放的严格共同集合。若自己的流程采用标准化，应该明确记录它发生在什么阶段，并在合理范围内检查结果对该选择是否敏感。

五、加权均值与样本均值回答的问题不同

我们专门设置一个异质位点。三个对照的甲基化计数与总覆盖分别为一除以十、九十除以一百、九十除以一百，三个处理则分别为五除以十、五十除以一百、五十除以一百。处理组在每个样本中都是百分之五十，对照组则是百分之十、九十和九十。

若按覆盖加权，对照组是总甲基化计数一百八十一除以总覆盖二百一十，约为百分之八十六点一九。处理减对照得到负三十六点一九个百分点。如果先算每个样本比例，再对三个样本取普通平均，对照组约为百分之六十三点三三，差异就成为负十三点三三个百分点。

weighted <- calculateDiffMeth(
  strict, overdispersion = "MN",
  test = "F", adjust = "BH", effect = "wmean", mc.cores = 1
)
average <- calculateDiffMeth(
  strict, overdispersion = "MN",
  test = "F", adjust = "BH", effect = "mean", mc.cores = 1
)

两次运行的 P 值和校正值逐项相同，改变的是报告的差异幅度。对于这个位点，校正值约为零点三八四，两种效应定义都没有通过本例的显著性门槛。不能因为加权差异超过二十五个百分点，就把它单独认定为可靠差异，更不能将效应定义变化误写成统计模型重新证明了差异。

百分比之差的单位是百分点。由百分之二十变为百分之三十，差十个百分点，与相对增加百分之五十是两种表达。设置差异幅度阈值时必须确认软件使用哪一种。实际报告也应说明处理减对照的方向，避免读者把负号解释成相反的生物学变化。

六、先合并重复，会丢失组内波动信息

接着我们将每组的三个重复合并计数，再运行适合每组一个合并样本的 Fisher 检验。异质位点的加权差异仍然是负三十六点一九个百分点，但校正值从约零点三八四变成约八乘十的负十五次方。差异幅度没有改变，统计证据却看起来强了许多。

原因可以直接从输入理解。合并后只剩对照组一百八十一比二十九、处理组一百零五比一百零五这两组计数，原来对照样本之间百分之十与百分之九十的差别已经看不见。模型获得了更多汇总读段，却失去了哪些读段来自哪个生物学重复的信息。

在本例预设的六十个大幅变化位点中，保留重复的流程按校正值小于零点零五且绝对差异超过二十五个百分点筛出五十八个。合并流程选出六十个预设位置，还额外选中上述异质位点。这些数量只描述本次构造数据，不能据此认定某种方法普遍更敏感或更准确。

合并技术重复与合并独立生物学重复也需要区别。前者是否合理，应结合实验设计与文库关系判断，后者会直接压缩组内变异信息。分析时不应为了增加显著位点数量而把独立重复提前合成一个样本。更有价值的做法是保留重复身份，并检查离群样本、批次及模型假设。

七、磁盘数据库模式改变存储方式，不替代结果核验

我们还把同一批输入写为 tabix 支持的磁盘对象，再执行相同过滤与共同位点合并。结果得到四百五十九行，逐列与内存模式一致。这个检查证明本次两条处理路径在给定输入和参数下产生相同数据，而不是证明数据库模式在任何规模下都更快。

磁盘方式适合在数据规模较大时按块处理，但压缩文件、索引和对象引用需要一起管理。如果只复制对象而没有保留它引用的数据文件，换一个目录或环境后可能无法继续使用。归档时应把这些依赖作为同一个分析产物保存，并用一次实际读取检查归档完整性。

每次演示使用新目录也有意义。沿用旧数据库文件名可能让人难以分辨读取的是本次生成的数据还是已有文件。我们保存了输出文件清单，便于核对实际产生的压缩文件与索引。这个例子规模很小，没有测量内存峰值或者吞吐量，因此不把存储模式的选择写成未经测量的性能结论。

八、工具的优势与使用代价

methylKit 的优势是从样本计数、覆盖过滤到差异检验提供连贯的数据结构，也允许保留重复、加入协变量并使用磁盘对象。对需要逐步审阅甲基化流程的读者，这种明确的中间对象很有帮助。它的使用代价，是很多看似简单的参数其实会改变分析集合或者效应含义。

真实项目中，覆盖门槛、共同位点规则、差异幅度定义和多重检验范围应在分析记录中一起说明。本文显式使用 BH 校正，校正所对应的是进入本次检验的四百五十九个位点。若只导出显著位置再重新校正，或者不同方法各自在不同位点集合上校正后直接比较数字，得到的结果不再回答同一个问题。

九、复现后最应该保留哪些经验

下载本组完整脚本与实跑记录，在本篇目录运行下面命令，可重新生成六个计数文件，并完成内存和磁盘两条路径。脚本固定镜像摘要与随机种子，保存 R 代码、版本、读取数量、统计结果和原始日志。执行环境限制为两核、四 GiB 内存，未据此推断真实项目的资源需求。

python3 run.py output

复现时先核对首行坐标与输入行数，再核对逐样本过滤数量、共同集合大小和示例位点的手算差异。最后再看显著位点列表。这个顺序能把文件读取、样本缺失、效应定义和统计推断分开检查，也能让分析报告解释结果为什么改变，而不是只展示哪次运行得到了更多位点。

当样本规模扩大，需要持续保存多批计数、模型与复现记录时，可通过云生信生物分析服务器服务组织统一的分析环境。