software

DiffBind 实跑复盘,共识峰门槛提高后,为什么单组特有的信号消失了

作者:发布于 2026-11-15约 13 分钟阅读
门槛提高,哪些信号被提前排除了。六个模拟样本 · 共识区域、计数与报告尺度。300 → 260,支持门槛由两个样本提高到四个,本例 40 个单组出现区域不再进入集合。77 / 96 / 105,未被调用为峰,仍然有比对计数,共同区域需要在所有样本重新计数。报告尺度,原始计数导出也会改变 Fold 含义,P 值仍来自原来的归一化分析。

差异结合分析中,一个区域没有进入最终差异表,可能是它确实缺少统计证据,也可能是它在计数和检验之前就被共识峰规则排除了。后者与“不显著”含义不同。如果只检查最后一张结果表,容易把分析范围的变化误解成生物学信号不存在。

本文实际安装并运行 DiffBind 3.20.0,构造六个样本的单端比对文件和峰列表,完成共识区域比较、重新计数、归一化、设计对比与差异分析。底层使用 DESeq2 1.50.2,运行环境为 R 4.5.2。所有数据均为教学模拟,没有真实样本,没有实际调用峰,也没有做实验层面的验证。

门槛提高,哪些信号被提前排除了。六个模拟样本 · 共识区域、计数与报告尺度。300 → 260,支持门槛由两个样本提高到四个,本例 40 个单组出现区域不再进入集合。77 / 96 / 105,未被调用为峰,仍然有比对计数,共同区域需要在所有样本重新计数。报告尺度,原始计数导出也会改变 Fold 含义,P 值仍来自原来的归一化分析。
六个模拟样本 · 共识区域、计数与报告尺度

一、先让实验设计和输入身份对得上

六个样本分为对照与处理两组,每组三个。脚本名称中,C1 到 C3 是对照,T1 到 T3 是处理。每个样本都有自己的峰列表、比对文件、索引和样本标识。样本表明确记录条件与重复编号,不靠文件夹顺序推断哪几个属于同一组。

模拟一共定义三百个候选区间,每个长五百碱基,相邻区间充分分开。前二百六十个区间在所有样本的峰列表中出现。随后二十个只写入处理组的三个峰列表,最后二十个只写入对照组的三个峰列表。因此,每个单独的峰文件都有二百八十行,但合起来的候选集合是三百个。

比对数量由固定随机种子生成。后四十个区域分别设置了处理组升高或降低的均值,其余区域不设置组间均值改变。随后根据每个区域的计数写入相应数量的人工比对,另外给每个样本加入一千条位于所有峰之外的背景记录。我们因此同时知道输入数量和应当被统计的区域内数量。

峰列表采用脚本声明的四列原始文本格式,依次是染色体、起点、终点和分数,起止按本例的一起点闭区间定义。读取后逐个核对范围宽度为五百。不要把这里的原始格式直接改名为 BED 后交给其他工具,因为文件后缀不会自动转换坐标体系。

二、至少两个样本有峰,不表示两个条件都必须有峰

我们先取全部六个峰集,要求至少两个峰集支持,实际得到三百个共识区域。将条件提高到至少四个峰集支持,结果变为二百六十个。被移除的四十个区域,恰好就是只在单组的三个峰列表中出现的那两批。

区域类别 支持峰集数 至少两个 至少四个
两组共同出现 6 保留 260 保留 260
仅处理组出现 3 保留 20 排除 20
仅对照组出现 3 保留 20 排除 20
合计 不同 300 260
cons2 <- dba.peakset(d, 1:6, minOverlap=2, bRetrieve=TRUE)
cons4 <- dba.peakset(d, 1:6, minOverlap=4, bRetrieve=TRUE)

这里的门槛施加在传入的全部峰集上,不会自动改写成每个条件内部的重复要求。在六个样本的整体集合中要求四份支持,本例中的任何单组特有区域都不可能通过,即使它在该组的三个重复中全部被调用为峰。

因此,提高门槛并不总是单纯提高可信度,也会改变分析对象。若研究问题关注某条件获得或丢失的信号,应先思考是否希望保留单组稳定出现的区域。若问题只关注两组共同存在的区域中信号强度如何变化,则较严格的交集可能符合目的,但报告应说明这个范围限制。

真实数据的峰边界还会有所偏移,重叠和合并会进一步影响共识区间。本例让边界完全一致,只用于隔离支持数量这个因素,没有评价不同峰调用参数、宽峰连接或边界模糊时哪种共识策略更优。实际项目可以导出两种候选集合,检查丢失的是孤立噪声,还是正好集中在某个条件的重复峰。

三、峰列表里没有,不等于比对计数为零

区域 P261 只存在于处理组三个峰列表中,但在对照组的比对文件里仍然有信号。实际重新计数得到,对照三个样本分别为七十七、九十六、一百零五,处理三个样本分别为二百四十二、二百四十一、二百四十五。

如果在合并峰表时,将未被调用为峰的样本填成零,就会把连续的覆盖强度变成错误的零值,夸大组间差别。峰调用结果表示一个区域是否通过了该次调用的规则,重新计数则询问同一套区域中各样本有多少记录。两者承担不同任务,不能互相替代。

我们在三百个共识区域上对六份比对重新计数,并逐个检查得到的整数与生成时保存的数量矩阵完全一致。这一步涵盖一千八百个区域与样本组合。核验对象是每一格原始计数,而不是只检查总和相等,因为行顺序错位时,总和仍可能看起来正常。

d <- dba.count(
  d, minOverlap=2, score=DBA_SCORE_READS,
  summits=FALSE, filter=0,
  bUseSummarizeOverlaps=FALSE, fragmentSize=50,
  bRemoveDuplicates=FALSE, bSubControl=FALSE,
  bParallel=FALSE
)
raw <- as.data.frame(dba.peakset(d, bRetrieve=TRUE))

本例使用内部计数路径,读段和指定片段长度均为五十,且都放在目标区间内部,不进行峰顶重新定位。这些设置是为了让输入计数可核对,不是对真实双端数据、重复过滤或宽峰分析的统一建议。尤其不能将本例的单端读取行为直接当成真实双端片段计数的验证。

四、归一化的分母要能追溯到输入

计数后,我们明确选择按完整文库大小进行缩放,记录每个样本用于归一化的总读段数。本次前两个对照样本分别为 32440 和 32504,第一份处理样本为 36754。这些数字与生成时区域内读段总数加上一千条背景记录相等。

处理组文库大小较大,不是因为它的峰列表行数更多。两组每个峰列表都是二百八十行,而比对数量受到各区间模拟计数的影响。把峰数量当作测序深度分母,会混淆二值的峰调用结果与连续的文库信息。

d <- dba.contrast(
  d, design="~Condition",
  contrast=c("Condition", "treated", "control")
)
d <- dba.normalize(
  d, method=DBA_DESEQ2,
  normalize=DBA_NORM_LIB,
  library=DBA_LIBSIZE_FULL,
  background=FALSE
)
norm <- dba.normalize(d, method=DBA_DESEQ2, bRetrieve=TRUE)
d <- dba.analyze(d, method=DBA_DESEQ2, bParallel=FALSE)

归一化记录还保存了实际缩放因子,六个值约在零点九四到一点零八之间。这样,发现某个样本异常时,可以从分母和因子入手排查,而不是仅凭最终显著区域的数量选择方法。本次没有运行外源参照或背景分箱归一化,也没有比较不同归一化方法的优劣。

如果真实实验预期存在整体信号升降,单纯依靠总量缩放是否适合,就需要结合实验设计与参照信息讨论。本文的小数据无法回答这一问题。它能给出的实用经验是,选择方法后应把方法名、分母来源及实际因子一起保存,使结论依赖的尺度假设能够被检查。

五、明确比较方向,再解释正负号

本次设计显式写成条件因素,并明确比较“处理相对于对照”。导出的对比摘要显示,两边各有三个样本。这一步值得单独检查,因为样本表中的组名存在,并不意味着结果表的第一组必定是你想象的处理组。

在归一化报告中,P261 的效应约为正一点二,P281 约为负二点三七。前者在处理组更高,后者在对照组更高。我们又建立反向对比,使用相同计数范围和归一化,重新运行分析。结果中所有区域的效应符号反转,而原始概率值保持一致,脚本对这一关系作了逐行核对。

这里的数值是以二为底的对数效应,不是普通倍数。一点二对应的倍数需要通过二的相应次方还原,不能写成“提高一点二倍”。判断方向时,也不要只看列名中有正负号,而应把对比定义、组别浓度和一个已知方向的区域放在一起检查。

导出包含样本计数的报告时,样本列按比较两组组织。我们的原始矩阵先放对照、再放处理,而正向报告先列处理、再列对照。如果后续用固定列号挑选样本,很容易在没有报错的情况下交换分组。按明确的样本名称取列,通常比依靠列位置更稳妥。

六、报告原始计数时,Fold 列也会改变

实际核对中,一个容易忽略的细节是,关闭报告归一化并不只改变末尾的样本计数列。本次同时导出归一化报告和原始尺度报告,P261 的 Fold 从约一点二变为一点三八九,P281 从约负二点三七四变为负二点二五八。两份报告对应区域的原始概率值和校正值却逐一相同。

区域 归一化报告效应 原始尺度报告效应
P001 0.075 0.379
P261 1.200 1.389
P281 −2.374 −2.258
normalized <- dba.report(
  d, th=1, bCounts=TRUE, bNormalized=TRUE,
  bCalled=TRUE, precision=0
)
raw_scale <- dba.report(
  d, th=1, bCounts=TRUE, bNormalized=FALSE,
  bCalled=TRUE, precision=0
)

原因在于原始尺度报告根据原始计数计算组别浓度及其对数差,而归一化报告使用底层分析提供的效应。因此,导出原始尺度之后,不能把其中的效应继续描述成与另一张表完全相同的模型效应。统计显著性仍取自既有的归一化分析,并没有因为显示原始计数而重新进行未归一化检验。

以 P261 为例,处理组原始计数均值约二百四十二点七,对照组约九十二点七,两者之比取以二为底的对数,得到约一点三八九,与原始尺度报告对应。这个可以手算的核对,比仅凭表头含有同一个 Fold 名称就认为含义相同更可靠。

在交付中,可以同时保留一张用于核对原始计数的表和一张用于解释模型效应的表,并在文件名中标明尺度。比较两个分析版本时,也应确认报告参数一致。否则,一次单纯改变导出选项的操作,就可能被误写成重新归一化之后效应发生变化。

七、通过统计门槛,仍然不是模拟真值的证明

在三百个候选区域上,本次有四十四个区域通过百分之五的校正阈值。其中四十个属于生成时人为设置变化的区域,另外四个来自未设置均值变化的区域。这个结果应如实保留,不能为了让演示看起来整齐,把额外通过的区域删除。

这只是一套固定随机数据和一个明确分析设置的结果,不是多次模拟的误差率评估。不能用一次实验中四十四这个数字宣称工具准确率,也不能把显著性阈值理解成每次结果名单中必定只有固定比例的错误。若要评价统计性质,需要另行设计重复模拟、参数范围与评价标准。

我们还保留了全部三百个区域的报告,而不只保存通过阈值的四十四行。这样可以区分效应小、变异大、校正后证据不足,以及压根没有进入候选集合等情况。本次四支持门槛的二百六十个区域只用于共识集合对照,没有在那个集合上重新做一套差异检验,因此不能声称其最终显著区域数是多少。

八、把分析拆成能够检查的几层

这次实跑形成了一个可迁移的检查顺序。先检查样本表,确认哪个重复属于哪个条件。再检查共识范围,确认研究关心的单组区域是否仍在。随后检查原始计数,确认没有把未调用峰错误填零。最后检查归一化、比较方向及报告尺度,再阅读显著区域名单。

DiffBind 的优势是把峰集、比对计数、实验设计和差异结果组织在同一流程中,能够保留原始峰是否出现的信息。它的局限也在于前面的选择会持续传到后面的结果。一个区域被选入,表示它有资格接受当前分析,不表示它必定差异显著。一个区域被排除,也不表示已经证明两组相同。

本例没有匹配的输入对照比对文件,也没有适用于虚构染色体的黑名单或灰名单,因此关闭了相应自动处理。原始计数与生成数量完全一致,是本例计数核验的证据。结果记录中的某个配置字段不能独立证明做过某项实验处理,例如本次归一化元数据中的控制扣除标记为真,也不能据此宣称分析实际使用了不存在的输入对照。

下载本组三篇复现脚本与实测记录,在有 Docker 的 Linux 环境运行下面的命令,即可生成六份比对与峰列表,重跑全部步骤。脚本固定镜像摘要,并保存样本表、预期计数、共识集合、原始计数、两种报告尺度、反向对比和环境信息。固定环境用于复现本次结果,不代表最新软件组合。

python3 run.py output

进一步可查阅 对应 Bioconductor 版本的软件页、DiffBind 分析教程 和 参数手册。后两份发布版文档可能继续更新,本文的实测数字以随附版本和记录为准。

当候选区域、计数规则和比较设计确定后,若要开展更多样本的重复计数与差异分析,可携带样本表、峰文件格式和任务规模,通过云生信生物分析服务器服务评估运行环境与批量处理安排。

把知识用于当前任务

下一步可以直接开始分析

按研究目标查找在线工具、课程、计算资源或完整分析服务。

按目标找方案仍有疑问,联系客服 ↗