DSS 实跑指南——差异位点、效应阈值与区域合并，怎样避免把结果读错

差异甲基化分析常常从一张位点表开始，最后变成一张区域表。这个过程看似只是把相邻显著位置连起来，实际上涉及差异方向、平滑范围、效应门槛和区域组成等多个判断。若只关注最后找到多少个区域，很容易把参数造成的边界变化误当成新的生物学发现。

本文安装并运行 DSS 2.58.0，配合 bsseq 1.46.0 和 R 4.5.2，先用六个模拟计数样本完成模型拟合，再用另一张人工结果表单独检查区域合并规则。两部分数据用途不同，全文明确区分。所有数字来自实际执行，不代表真实疾病、组织或处理条件下的甲基化变化，也不构成工具准确率比较。

一、输入应保留覆盖和甲基化计数

官方使用指南将甲基化分析输入组织为染色体、位置、总覆盖和甲基化计数。总覆盖通常记为 N，甲基化计数记为 X。必须满足甲基化计数不大于总覆盖，不能把百分比或者经过任意缩放的连续数直接当作原始计数。

只提供比例会丢失深度信息。例如，一次覆盖中的一次甲基化，与一百次覆盖中的一百次甲基化都等于百分之百，但测量的不确定性不同。另一方面，覆盖再高也不能取代生物学重复，因为读段数量没有告诉我们不同个体或独立样本之间如何变化。

本例设置两组各三个样本，一千二百个等间距 CpG，每个位点每个样本有三十到六十次覆盖。先生成样本间有波动的甲基化概率，再生成整数计数。大部分位置没有预设组间均值变化，两个连续区块分别设置较大的升高和降低，另一个区块只设置十个百分点的较小变化。

两个大幅变化区块共七十个位点，小幅变化区块有三十个位点。它们是便于核验输出的构造标签，不是从结果倒推出来的“真实发现”。随机数据中也会出现未预设的波动，因此检出的数量与预设数量不同，需要结合规则解释，不能直接换算成真实应用中的灵敏度。

二、先固定组别方向，再看正负号

我们将三个对照设为第一组，三个处理设为第二组，分别运行不平滑和平滑两种分析。平滑版本使用五百碱基的范围，并将其他主要设置保持一致。本文不会把某一组设置称为适合所有甲基化数据的默认最佳方案。

bs <- makeBSseqData(dat, ids)
plain <- DMLtest(
  bs, group1 = ids[1:3], group2 = ids[4:6],
  smoothing = FALSE, ncores = 1
)
smooth <- DMLtest(
  bs, group1 = ids[1:3], group2 = ids[4:6],
  smoothing = TRUE, smoothing.span = 500, ncores = 1
)

这里报告的差异是第一组减第二组。因此，处理组更高的区块出现负值，处理组更低的区块出现正值。为了验证方向，我们又交换两组运行一次，确认所有位点的差异恰好反号，而相应 P 值保持一致。这个检查比根据某个工具的习惯猜测正负号更直接。

实际报告中，应把组名与方向写在列名说明或者图注中。不同工具可能使用处理减对照，也可能使用第一组减第二组，不能因结果列都叫差异就默认方向一致。后续基因注释、区域着色或富集分类一旦沿用错误方向，整套解释会被倒转。

三、差异不为零与差异超过门槛是不同问题

我们对同一份不平滑结果分别运行两种位点筛选。第一种把效应门槛设为零，按 P 值小于零点零零一选择。第二种把门槛设为零点二，并使用零点零五的筛选阈值。零点二在这里表示二十个百分点的甲基化比例差异，不是相对增加百分之二十。

dml_zero <- callDML(plain, delta = 0, p.threshold = 0.001)
dml_effect <- callDML(plain, delta = 0.2, p.threshold = 0.05)

安装版函数的实际处理显示，当效应门槛大于零时，筛选依据变成估计的绝对差异超过该门槛的后验概率。此时参数名虽然仍然带有 P，但它用于比较的是一减该概率。阈值零点零五对应超过门槛的后验概率大于零点九五，不能再把它描述成普通零差异检验的原始 P 值阈值。

这个概率建立在当前拟合结果及其不确定性估计之上，不是样本发生某种生物学机制的概率。它回答的是在模型下差异是否有足够证据超过指定幅度。样本设计有混杂、输入计数有问题或者模型假设不合适时，数字本身不会自动消除这些问题。

还可以采用另一种筛选方式，即先按位点校正值限制错误发现率，再要求估计差异超过二十个百分点。本例这种组合得到七十五个位点，而上述后验门槛选择了七十个。两者不是同一个筛选规则，因此数量不同并不矛盾。报告必须说明具体使用哪一种，而不能只笼统写“差异大于二十且显著”。

四、平滑既可能增强连续信号，也可能改变边界

不平滑时，零差异检验筛出七十三个位点，其中七十个位于预设大幅变化区块，另外三个在这些区块之外。小幅变化区块没有位点通过该门槛。采用二十个百分点的后验门槛时，留下七十个位点，均来自两个大幅变化区块。

平滑后，零差异检验筛出一百三十八个位点，包括七十个大幅变化位点、三十个小幅变化位点以及三十八个其他位置。后验效应门槛则留下七十三个，其中七十个来自大幅变化区块，另外三个位于其他位置。平滑增加的位点不应全部解释为额外发现，部分变化与邻近位置共享信息有关。

这个实验揭示了一个容易忽略的区别。连续但幅度较小的变化，可以在平滑后获得较强的统计证据，却仍达不到预先指定的效应幅度。显著性、效应大小和空间连续性是三个维度，不能用其中一个代替另外两个。

平滑也会将区块内部与边缘的估计联系起来。本例第一个大幅变化区块位于五千零一到五千五百八十一，未平滑的区域结果覆盖这个范围。平滑后按零差异规则识别的相应区域扩展到四千八百二十一至五千七百四十一。边界越宽并不自动意味着真实调控区更长，它也反映所选平滑与区域调用规则的共同作用。

五、区域结果需要与位点结果分别审阅

在未平滑结果上，我们使用最小长度五十、位点数门槛三、合并距离一百和显著比例门槛零点五进行区域调用。零差异规则与二十个百分点规则都得到两个区域，分别包含三十和四十个 CpG。这与本次构造的两个大幅变化区块对应，但不能据此推断真实样本也会如此整齐。

dmr <- callDMR(
  plain, delta = 0, p.threshold = 0.001,
  minlen = 50, minCG = 3,
  dis.merge = 100, pct.sig = 0.5
)

区域表包含起止位置、长度、位点数量、两组均值、差异以及区域统计量。本次输出的长度按终点减起点再加一计算。例如五千零一到五千五百八十一对应五百八十一碱基。若要转换为零起点、右侧开放的 BED 区间，应对起点做相应转换，而不是只将两列复制过去。

区域内均值与位点平均有关，不能自动理解为把全部样本读段汇总后得到的覆盖加权比例。区域统计量也不是显著性概率。本次输出中没有区域层面的 FDR 列，不能因为输入位点表有校正值，就把结果文件标题写成“区域 FDR 小于某值”。

六、用人工结果表看清合并规则

为了将区域规则与前面的模型估计分开，我们另写一张十七个位点的结果表，位置从一百到二百六十，每十个碱基一个。前两段各四个位点为正向差异，最后一段四个位点为负向差异，段间五个位点只有很弱的统计证据。这里的均值、差异和统计量是人工指定的，不是前面六个样本拟合得到的输出。

当最小长度为三十、合并距离为二十时，三个区段分别保留，范围是一百到一百三十、一百六十到一百九十，以及二百三十到二百六十。每个区域包含四个位点，长度三十一。把合并距离改为三十后，前两段连接成一百到一百九十，连同中间两个弱位点共十个位点。

本次还发现一个不能只看参数名猜测的细节。保持最小长度三十，将合并距离改成一百，结果仍与合并距离三十完全相同。检查安装版函数实现可见，它会在合并距离超过最小长度时，将前者限制到后者。因此，名义上把距离调大，并不总意味着实际允许更远的区域连接。

这里最后一行同时改变了两个参数，不能把全部变化单独归因于合并距离。它的作用是展示一组较宽松连接规则的实际后果。查看函数实现还发现，该版本最终长度和位点数的判断使用严格大于门槛。因此，临界位置是否保留，应核对安装版实现并用小例子验证，不能只按参数名字推断。

七、一个区域也可能含有相反方向的位点

在最小长度和合并距离均为一百时，人工表中的十七个位点合成一个长一百六十一碱基的区域。区域内既有八个正向显著位点，也有四个负向显著位点，还包含五个弱位点。区域平均差异只有约六点一八个百分点，明显小于各显著区段分别设定的二十五个百分点。

区域统计量也受到方向抵消影响。八个正向位点各贡献五，四个负向位点各贡献负五，再加上五个弱位点各零点二，合计为二十一。若只看区域统计量排序，很难直接知道内部有这样明显的方向分裂。因此，对候选区域应同时查看位点方向、各样本比例和覆盖分布。

将显著位点比例门槛提高到零点八后，这个合并区域不再输出。因为十七个位点中只有十二个达到设定的位点证据门槛，比例约为百分之七十点六。这个比例限制是对区域内部组成的要求，并不是区域层面的错误发现率控制，两者不能混用名称。

这也是为什么区域注释不应立即等同于机制结论。一个区域与启动子重叠，只说明空间关系。区域是否方向一致、是否由少数离群样本驱动、不同 CpG 的覆盖是否可比，都需要继续检查。把区域先贴到基因名上，再忽略内部位点结构，容易让一个复杂结果看起来过于整齐。

八、DSS 的优势与适用边界

DSS 的优势是将计数深度、重复间变化和差异估计结合起来，并提供从位点到区域的操作接口。平滑对连续变化的分析也很实用。其使用代价是用户必须分清模型输出、位点筛选和区域构造三个阶段，尤其不能把平滑后的更多位点自动视为更好的结果。

复杂实验设计还需要考虑批次、配对关系和其他协变量。本文只运行了简单两组比较，没有验证多因素模型，也没有验证无生物学重复的分析。如果研究问题需要这些功能，应重新核对相应模型的输入和支持范围，不能直接把两组比较的参数解释复制过去。

对于正式研究，建议预先写明组别方向、是否平滑、效应门槛、位点筛选规则及区域参数，再检查合理参数范围内结论是否稳定。稳定性检查的目标是辨认哪些结论依赖特定设定，而不是不断修改阈值直到得到满意的区域数量。

九、怎样复现实跑并检查自己的结果

下载本组完整脚本与实跑记录，在本篇目录执行下面命令，可重新生成六个计数文件、两种模型结果及人工区域练习。脚本固定随机种子和容器摘要，保留安装版本、分组交换检查、区域表、函数实现与日志。容器限制为两核、四 GiB 内存，没有进行大规模耗时或内存基准测试。

python3 run.py output

运行后可以先确认交换组别是否只改变差异方向，再对照零差异与效应门槛的入选数量，最后检查人工区域内部的十七个位点。把这三步应用到自己的项目，能够分别发现方向错误、筛选含义混淆和区域内部不一致的问题。它们比单独输出一个差异区域总数，更能帮助读者判断结果是否值得继续研究。

当分析扩展到更多样本和全基因组区域，需要保留模型、参数与位点审阅记录时，可通过云生信生物分析服务器服务开展统一的批量分析。