MethylDackel 实跑指南——重叠读段、质量阈值与 CpG 合并，怎样改变甲基化计数

拿到甲基化提取结果后，最容易直接查看的是百分比。然而，同样一个百分之百，可以来自一个通过筛选的片段，也可以来自大量独立片段。它们在数值上相同，证据强度却相差很大。若连分母中的“一个”代表什么都没有确定，后续差异检验就可能建立在错误的计数单位上。

MethylDackel 位于比对与统计分析之间，负责根据参考序列、比对记录和过滤条件汇总甲基化证据。本文安装并运行了 0.6.2 版本，用十六条人工构造的比对记录分别检查配对重叠、碱基质量、比对质量、重复标记和链合并。这里没有运行读段比对器，没有真实样本，也不评价甲基化检测准确率。小数据的价值，是让每一条输出都能追溯到明确的输入条件。

一、先把输入边界检查清楚

官方项目说明要求提供按坐标排序并建立索引的比对文件，以及对应的参考序列索引。实际工作中，文件能够打开只是第一步，还要核对染色体名称、长度、参考版本和比对器使用的转换链信息。参考版本错误未必总会导致程序退出，却可能让位点上下文发生变化。

本例构造八条各三百碱基的参考序列，每条只有一个用于检查的 CpG。正链胞嘧啶的零起点坐标是第一百位，相邻鸟嘌呤位于第一百零一位。七个场景各放入一对读段，第八个场景放入两个来自相反转换链的独立单端记录。因此，输入共有七对加两个单端记录，而不是八对双端数据。

比对文件中的位置、配对标记、质量值和转换链标签均由脚本明确写入，再实际使用 Samtools 排序和建立索引。这个设计刻意绕开比对不确定性，只观察提取程序如何处理已给定的记录。将这种结果解释成某种比对流程的准确率，超出了实验能支持的范围。

samtools faidx reference.fa
samtools sort -o input.bam input.sam
samtools index input.bam
MethylDackel extract -@ 2 -o default reference.fa input.bam

二、先看实际列数，再写读取脚本

本次默认 CpG 文件实际有六列，依次是染色体、零起点区间起点、区间终点、甲基化百分比、甲基化计数、非甲基化计数。对于一个未合并的胞嘧啶，区间是一碱基长。下面以空格展示字段关系，真实文件使用制表符分隔。

agree 100 101 100 1 0
strand_pair 100 101 100 1 0
strand_pair 101 102 0 0 1

不能只凭扩展名把第四列读成覆盖度。本次百分比模式下的第四列是一百或者零，后两列才是计数。我们还发现，项目页面中部分输出说明与这个安装版本的实际文件列数并不完全一致。因此，导入规则应同时依据安装版本、帮助信息和少量实测行进行确认，不能把某段网页描述长期写死在数据处理脚本中。

最实用的检查是选取一个有明确计数的位点，手算甲基化计数除以总计数，再与输出百分比核对。如果文件换成比例模式，第四列又会进入零到一的尺度。字段意义和单位是两个独立问题，列数正确并不意味着数值尺度也正确。

三、一对重叠读段只提供一次位点计数

第一对读段在目标位置都报告胞嘧啶，碱基质量均为三十。它们来自同一片段，重叠部分包含目标位点。默认输出是一个甲基化计数、零个非甲基化计数，而不是两个甲基化计数。这个结果说明，直接对重叠读段逐条累加会把同一片段的证据重复使用。

第二对读段在同一位置互相冲突，一端是胞嘧啶，一端是胸腺嘧啶，而且质量同为三十。该位置没有留下合格调用。不能把它写成一个甲基化和一个非甲基化计数，再报告百分之五十，因为这样会把同一片段内部的测量冲突解释成两个独立分子的混合。

第三对也互相冲突，但胞嘧啶质量为四十，胸腺嘧啶质量为三十五。按本次处理结果，高质量一端留下的有效质量为五。默认碱基门槛下保留一个甲基化计数，把门槛提高到十后则不再输出该位置。质量值原本都不低，却仍可能在重叠冲突处理后无法通过筛选。

这里的经验不是一律把门槛调低，而是遇到覆盖数减少时，先检查减少发生在哪个环节。原始读段覆盖、比对筛选后覆盖、重叠处理后有效调用可以不同。把这几个量分别记录，才能解释为什么一个看起来深度充足的位置，在甲基化表里只剩很少证据。

四、比对质量、碱基质量和重复标记各管一件事

我们另设三个场景，分别降低比对质量、降低目标碱基质量，或者给读段设置重复标记。它们在默认结果中均没有留下目标位点。随后每次只修改一种条件，观察哪一个场景重新出现，从而避免同时放宽多个参数后无法判断原因。

比对质量为五的记录，在设置最低比对质量为零后进入结果。目标碱基质量为三的记录，在最低碱基质量改为一后进入结果。带重复标记的记录，则在启用保留重复的选项后进入结果。三个变化对应三个不同的筛选层次，不能用一个“数据质量较低”概括所有原因。

MethylDackel extract -q 0 -o mapq0 reference.fa input.bam
MethylDackel extract -p 1 -o base1 reference.fa input.bam
MethylDackel extract --keepDupes -o keepdupes reference.fa input.bam

保留重复标记的选项，并不负责判断哪些记录应该被标记为重复。本文也没有运行重复识别程序，只是在人工记录上设置了这个标志。真实项目需要先说明重复标记从哪里来、是否符合文库类型，再决定提取时是否排除，不能把“保留重复”写成某种增加真实覆盖的通用技巧。

本次还专门尝试把最低碱基质量设为零。程序没有以非零状态退出，而是输出提示，将该值重设为一，并继续完成运行。结果与显式设为一相同。这说明检查任务成功与否不能仅看退出码，还应保存标准错误并核对实际生效参数。自动化任务若只在退出失败时通知，可能漏掉这种参数被修正的情况。

五、CpG 链合并会改变深度的计算对象

第八个场景包含两个独立单端记录。一个支持正链胞嘧啶甲基化，另一个支持相反链对应胞嘧啶未甲基化。这里不是一对重叠读段，而是两个独立记录，必须与前面的配对冲突场景区分开。

默认输出将两条链对应的胞嘧啶分别列出，每条的有效深度都是一。启用 CpG 上下文合并后，输出区间覆盖相邻两碱基，甲基化与非甲基化计数各为一，百分比变成五十。变化来自计数对象由单个胞嘧啶转成合并后的 CpG，而不是原始测序数据增加了。

MethylDackel extract --mergeContext -o merged reference.fa input.bam
MethylDackel extract -d 2 -o depth2 reference.fa input.bam
MethylDackel extract --mergeContext -d 2 -o merged_depth2 reference.fa input.bam

未合并时，最低深度设为二，整个结果没有任何数据行。合并后再设同样门槛，第八个场景保留一行，因为两条链合起来达到两个有效调用。本次结果与安装版帮助中所述的合并后应用深度门槛一致。比较两个项目的“至少两次覆盖”时，如果链合并设置不同，实际上比较的可能不是同一个过滤标准。

还要注意，合并后两条链各一次支持，与同一条链有两个独立片段并不是完全相同的实验覆盖结构。是否适合合并，应结合所关注的甲基化上下文和后续模型要求决定。若分析需要保留链间不对称信息，过早合并会把这种差异压缩成一个平均值。

六、缺失行、零比例和零覆盖要分别解释

默认输出中没有某个位点，可能是没有读段覆盖，也可能是读段被比对质量、碱基质量或标记筛掉，还可能是配对冲突导致没有合格调用。单凭输出表缺一行，无法判定这个位点完全未甲基化，更不能在跨样本合并时直接补成百分之零。

本次胞嘧啶报告模式输出十六行，对应八个 CpG 的两条链。其中有些行的甲基化计数与非甲基化计数都是零。这样的零除以零没有可解释的甲基化比例，应当作为缺乏有效观测处理。与之不同，一个甲基化计数为零、非甲基化计数为一的位置，才具有一次支持未甲基化的调用。

比例模式又带来另一种区别。启用比例输出后，本次文件变为四列，第四列使用一与零，不再保留两种计数。这样的文件适合某些展示用途，但若后续统计需要计数及覆盖信息，仅保留比例会损失必要输入。百分比、一位小数的比例和整数计数之间，不能靠简单改列名互相替代。

七、M-bias 图应该服务于具体的修剪决定

本次实际运行了 M-bias 子命令并生成相应图文件，但十六条人工记录不足以估计真实文库的末端偏倚，所以本文不从这些图推导修剪长度。实跑一个命令，与得到足够支持某项实验判断的数据，是两件需要分别说明的事。

真实数据中的末端异常，应同时检查转换链、读段端别、覆盖数量和偏离是否在多个样本中一致。只有极少数读段贡献的末端波动，与大量覆盖下持续出现的系统偏移，处理依据不同。把某篇教程中的固定修剪长度直接套到所有文库，可能移除正常信号，也可能留下真正受影响的位置。

提取前后的位点数量可以帮助判断修剪影响，但不能只追求保留更多。更好的做法是预先记录修剪依据，使用同一套明确规则处理可比样本，然后检查变化是否集中于预期位置。若不同批次需要不同处理，也应保存理由，避免在看到差异结果后才反向挑选参数。

八、这个工具适合解决什么问题

MethylDackel 的实用之处，是把重叠、质量和上下文合并等关键规则集中到明确的提取步骤中，便于建立可检查的位点计数表。它也提供了多种输出方式，可以适配不同后续任务。代价是用户需要真正理解每种输出的单位，不能以为所有以甲基化为名的表都能直接互换。

它不会替代文库设计、参考核验、重复识别或者差异统计。一个格式完整的结果，也不会自动说明转换充分、样本可比或生物学重复足够。选工具时，先确认上游比对格式受支持、链信息正确、下游需要计数还是比例，再考虑具体参数，比先列出一长串优化开关更有效。

九、怎样复现并迁移到自己的数据

下载本组完整脚本与实跑记录，在本篇目录执行下面命令即可重新构造参考和比对记录。脚本固定容器镜像摘要，保存每条实际命令、标准输出、标准错误和提取结果。容器限制为两核、四 GiB 内存，这是本例的执行边界，不构成全基因组数据的资源建议。

python3 run.py output

迁移到真实项目时，可以先抽取少量已知位点，人工核对读段重叠和过滤原因，再扩展到整个样本。重点保留版本、参考来源、重复处理、质量阈值、链合并状态和输出单位。这些记录能让后续差异分析明确知道输入中的一个计数到底意味着什么，也能在结果不一致时缩小排查范围。

当这套提取与核验流程需要覆盖多批全基因组样本时，可通过云生信生物分析服务器服务统一运行环境与批量任务记录。