同一个特征在十二个样本中的原始计数都等于一百,加入 TMM 因子后没有明显差异,只按文库总量处理时却出现显著下降。这并非改写了原始计数,而是同样的一百,被放在了不同的有效文库背景中解释。本次用明确构造的组成偏移,实际运行 edgeR,检查标准化到底改变了什么。
文章还安排了一个低表达过滤对照。相同输入只改变过滤时使用的分组或设计信息,就会出现不同的保留结果。读者可以据此检查自己的分析流程,避免把过滤、标准化和模型检验混成一条不可解释的黑箱命令。

这次实际运行了哪些部分
我们安装 edgeR 4.8.2 和 limma 3.66.0,在 R 4.5.2 环境中运行计数过滤、TMM 因子估计、准似然模型和效应门槛检验。输入矩阵、分组信息及结果表均由脚本生成或导出,没有真实测序样本,也没有拿模拟结果评价方法在真实项目中的准确率。
主实验有十二个样本,分成两个批次,每批各有三个对照和三个处理。这样可以在每个批次内部进行条件比较,同时在模型里包含批次项。矩阵有两千四百二十一行,其中两千行按负二项分布生成,包含预先设定的部分处理效应和批次效应。
另外加入二十条组成偏移特征,它们在对照中计数较低,在处理中占据大量计数。再加入一条所有样本都为一百的对照行,以及四百条极稀疏行。这种设计有意把文库组成改变得明显,便于理解总量与相对背景的区别。它没有模拟所有真实实验因素,不适合据此推算测序深度或检验功效。
edgeR 的版本、依赖和计数模型背景可在 Bioconductor 3.22 软件页面 查阅。本文使用的是固定版本,后续版本的默认参数变化应单独核对,不能只凭函数名相同就认为结果必然一致。
一个独立的小表,检查低表达过滤
在主模型之前,我们先构造六条计数模式。为单独观察过滤逻辑,小表把每个样本的文库大小人为固定为一百万,使默认计数门槛对应的每百万计数值容易核对。这是阈值演示,不是由这六行计数计算出的真实文库总量,也没有把这份人为分母用于主实验拟合。
按处理条件分组时,每组有六个样本。若传入同时包含批次与条件的设计矩阵,设计杠杆值对应的有效最小样本数在本例中变成四。两种输入的信息不同,程序使用的保留要求也不同。不能看到一个函数返回布尔值,就认为它执行的是永远不变的低表达阈值。
| 计数模式 | 仅按条件分组 | 传入完整设计 | 本例中看到的现象 |
|---|---|---|---|
| 十二个样本各有 12 | 保留 | 保留 | 普遍有表达 |
| 四个样本各有 12 | 不保留 | 保留 | 保留要求不同 |
| 三个样本各有 12 | 不保留 | 不保留 | 支持样本不足 |
| 六个样本各有 12 | 保留 | 保留 | 达到所需数量 |
| 一个样本有 1000 | 不保留 | 不保留 | 单点高计数不够 |
| 十二个样本各有 1 | 不保留 | 不保留 | 每个样本都太低 |
表中的四样本模式跨越对照和处理,并非某一组独有。六样本模式也包含两种条件。这说明过滤使用分组来确定所需的样本数量,不等于要求某一指定组内所有样本都达到门槛。解释过滤结果时,要分清用于确定门槛的信息和逐行计数检查本身。
单个样本计数很高的那行没有保留,也值得注意。总计数高并不一定说明这个特征适合进行组间比较。支持分散在哪些样本,与总数同样重要。遇到一行主要由单个样本贡献的情况,应查看数据分布,不宜因为它的总计数超过其他行就跳过过滤审阅。
过滤时不要机械地塞入所有协变量
完整模型需要包含哪些因素,与低表达过滤应依据哪些研究分组,不一定完全相同。官方指南特别提醒,作为阻断或校正使用、但并非研究目标的变量,不应机械地决定过滤分组,并用配对设计举例说明。参见 edgeR 用户指南的过滤部分。
因此,主实验根据处理条件执行过滤,拟合时仍包含批次。小表中传入完整设计的结果只是用于展示差异,没有被宣传成更优选择。分析人员应先说清楚希望在哪些条件之间获得可解释的比较,再决定过滤规则,而不是以保留更多行作为选规则的理由。
主矩阵过滤后保留两千零二十一行,四百条稀疏行被去除。随后重算保留矩阵的文库大小,再估计标准化因子。结果记录保存了过滤向量对应的数量与样本表,使别人可以追查某一行是没有进入模型,还是进入模型后没有达到显著门槛。
design <- model.matrix(~ batch + condition, samples)
y <- DGEList(cts)
keep <- filterByExpr(y, group = samples$condition)
y <- y[keep, , keep.lib.sizes = FALSE]
y <- calcNormFactors(y, method = "TMM")
TMM 没有改写这张原始计数矩阵
本次在调用因子估计前后,用严格相等检查比较计数矩阵,结果为真。所有原始计数保持不变。变化保存在样本信息中的标准化因子,随后与文库大小共同参与模型对背景规模的解释。把这一点验证出来,比只说已经标准化更具体。
可以将本实验中的有效文库大小理解为保留矩阵的文库大小乘以对应因子。十二个因子的乘积在本次输出中为一。这是尺度约束下的结果,不表示每个样本都必须有同样的因子,更不意味着估计出了每个样本的细胞数量或总核糖核酸分子数。
为了形成对照,我们复制同一份过滤后的对象,一份保留默认的单位因子,另一份加入 TMM 因子。前者仍然使用文库总量,并非完全忽略测序深度。文章把它称为只按文库总量处理,避免让读者误以为未加入 TMM 就等于模型直接比较未经任何尺度解释的数字。
两种对象都使用相同的设计和准似然检验。这样最终差异来自有效文库背景的处理,而不是更换了样本、过滤了不同的基因,或者修改了比较方向。控制住这些条件后,结果对照才具有清楚的解释。
同样的一百,在两种背景下得到不同结论
人工对照行在十二个样本中全部等于一百。加入 TMM 后,处理相对对照的对数倍数变化约为负零点零一一,原始 P 值约零点九五五,校正值约零点九九四。这个结果与我们设计的计数模式相符,未出现明显的相对变化证据。
只按文库总量处理时,同一行的对数倍数变化约为负一点二一零,原始 P 值约三点六一乘以十的负九次方,校正值约七点九八乘以十的负九次方。原始计数一项未改,统计结论却明显不同。原因是处理组中少数高计数特征抬高了文库总量,相同的一百在该背景下变得相对更小。
| 核对项目 | 加入 TMM 因子 | 单位因子对照 |
|---|---|---|
| 对照行原始计数 | 每个样本 100 | 每个样本 100 |
| 对数倍数变化 | −0.011 | −1.210 |
| 原始 P 值 | 0.955 | 3.6e-9 |
| 校正值 | 0.994 | 8.0e-9 |
这里不能进一步声称 TMM 能修复所有组成问题。我们有意构造了大多数背景特征没有共同大幅变化的情境,方法在这组数据里提供了更合适的相对参照。若研究问题涉及大范围整体表达改变、样本组成差异或绝对分子量变化,还需要额外设计和证据,不能从这一张表直接推到绝对表达结论。
另一条没有预设处理效应的随机特征也显示类似现象。加入因子时估计接近零,单位因子对照则偏向下降。这个例子说明问题会影响一批背景行,不只是专门加入的恒定行。但由于这里存在随机抽样误差,不能把每条未设效应的特征都要求为精确零值。
准似然检验与效应门槛检验回答不同问题
在当前固定版本中,我们使用准似然拟合,再检验处理条件的系数。默认问题是该效应是否偏离零。在校正值小于零点零五时,得到一百六十四条。这个数量包含特意设置的组成偏移特征,不能只和最初那部分处理效应行数比较后,就立即计算准确率。
接着对同一拟合对象运行效应门槛检验,把门槛设为以二为底的对数变化一,对应两倍变化的界限。本次有一百五十四条达到校正后的显著门槛。两次数量不同,是因为所问的问题变了,并不代表后一次程序漏掉了结果。
fit <- glmQLFit(y, design, robust = TRUE)
qlf <- glmQLFTest(fit, coef = "conditiontreated")
treat <- glmTreat(fit, coef = "conditiontreated", lfc = 1)
第一条特征的报告效应约为一点八四六,在普通检验中原始 P 值极小,而在效应门槛检验中约为零点零零零零一一六一。估计值看起来差不多,并不意味着两个检验的证据强度必须一样。不能把普通显著性列复制到一个标为超过两倍变化的报告中,就当作已经进行了相应检验。
结果中还出现了报告效应与未收缩效应并列的字段。它们应保留原名和来源,不宜为统一表格而强行覆盖。本实验没有执行跨方法的效应收缩比较,因此不会把不同列值之间的微小差异解释成某种方法优劣。处理软件输出时,首先应维护字段含义,再考虑展示格式。
遇到零值 P 值和很强的信号怎么办
人为设置的组成偏移特征非常极端,在本次输出里出现数值上的零 P 值。这不是事件为真的概率等于百分之百,也不是已经获得无限强的实验保证。计算机表示极小概率有数值范围限制,解读时应保留它与其他证据之间的区别。
画负对数图时,直接对零取对数会产生无穷值。展示程序可以设置明确的显示上限,但应另存原始结果,并在图注说明显示处理。不能为了让图看起来整齐,直接把原表中的零改成一个随意的小数,再让后续分析误以为这个数字由模型计算得出。
同样,极大的倍数变化应回到各样本计数分布核查。本次知道这些行来自人工设定,因此可以解释它们为何特别突出。真实项目没有这样的预先答案,必须检查样本结构、计数来源和异常模式,不能因为数值极端就省略复查。
复现时应保留哪些文件
完整资料见 本组实跑脚本与结果。驱动脚本会安装校验过的固定版本软件,在独立库中运行数据生成、过滤与两种标准化对照,并保存所有结果。进入 edgeR 对应目录后,可在有 Docker 的 Linux 环境中执行下面的命令。
python3 run.py run159
复现材料包括计数矩阵、样本表、过滤练习表、文库大小及因子,以及普通检验、效应门槛检验和单位因子对照结果。不要只留下排序后的前二十行,因为那样很难检查一个候选是在哪个阶段消失的,也难以确认结果排序是否混用了不同检验的字段。
本次依赖环境中出现了包由较新 R 补丁版本构建的提示,加载与实际拟合均完成,提示仍保留在日志里。这提醒读者,记录软件名称还不够,需要保存版本和会话信息。若更换运行环境,先核对关键计数、保留行数和代表性结果,再判断是否可以与旧分析直接比较。
优点、局限和适合的使用方式
edgeR 的优势是能把过滤、样本尺度、设计和对比明确地表达出来,并提供可核查的中间对象。对有合理重复、计数来源清楚、需要表达复杂比较的任务,它允许分析人员把统计问题写得具体。本篇中,计数不变但结论改变的过程,就能通过这些中间对象逐步追踪。
它的局限来自数据与问题本身。缺少有效重复时,模型不能自动补足实验信息。分组混杂、错标样本和不合理的参照,也不会因为函数正常结束而消失。稳健选项有明确作用范围,不代表所有异常样本都被识别,更不能代替实验质量判断。
读者最值得检查的三件事,是过滤依据是否与研究比较一致,标准化后保存的究竟是什么,以及显著性检验是否回答了真正关心的效应问题。把这些信息写进分析记录,能让后续结果讨论围绕事实展开,也能减少仅凭差异条数评价分析好坏的误区。
当这些计数分析流程需要在多项目中稳定运行时,可按实际数据规模了解 云生信生物信息学分析服务器服务。
