WGCNA 分块后，原来的共表达模块还在吗？用基因名单核对两次建网

第一次运行 WGCNA 得到一个和表型有关的模块。后来更改了 maxBlockSize，第二次仍然看到一个大小相近的模块。能沿用上次的解释吗？先问一个更容易回答的问题。原来属于这个模块的基因，这次分别去了哪里？只看颜色、编号或模块总数，都可能错过少量关键基因的移动，也可能把模块的重新编号误读成结果变化。

下面做一个可以完整复现的对照。模拟 48 个样本、220 个基因，固定输入、软阈值和聚类参数，只改变分块上限。然后用交叉表和成员重叠度找到模块的对应关系。模拟数据让我们知道哪些基因共享表达模式，适合练习读表。得到的模块和数值只用于说明方法。

先固定输入，再让两次建网只差一个参数

WGCNA 用样本之间的表达变化计算基因相关性，交给 blockwiseModules() 的矩阵应当一行一个样本、一列一个基因。真实 RNA-seq 项目要先归一化并处理明显的低表达或近乎无变化的基因。不要先按病例与对照的差异检验挑一批显著基因再建网络。这样会提前保留随分组一起变化的基因，使后面的模块与分组相关性更难解释。WGCNA 作者关于基因筛选的说明

下面生成三组带共同信号的变量，另加 45 个独立噪声变量。第二组与第三组的潜在信号有一定相关性，这样可以观察边界上的基因。把随机种子、样本名和基因名固定下来，后面即使两次输出的模块顺序不同，也能按名字对齐。

library(WGCNA)
set.seed(17)
n <- 48
z <- matrix(rnorm(n * 3), nrow = n)
z[, 3] <- 0.55 * z[, 2] + sqrt(1 - 0.55^2) * z[, 3]
signal <- function(x, k, noise) {
  matrix(rep(x, k), nrow = n) + matrix(rnorm(n * k, sd = noise), nrow = n)
}
datExpr <- cbind(signal(z[, 1], 70, 0.7), signal(z[, 2], 60, 0.8),
                 signal(z[, 3], 45, 0.8), matrix(rnorm(n * 45), nrow = n))
colnames(datExpr) <- sprintf("g%03d", seq_len(ncol(datExpr)))
rownames(datExpr) <- sprintf("s%02d", seq_len(nrow(datExpr)))
stopifnot(goodSamplesGenes(datExpr, verbose = 0)$allOK)

自己换表达矩阵时，用两个检查代替凭印象转置。先将行名与样本表中已知的样本编号匹配，再选一个熟悉的基因，看它是否恰好占一列。两次比较要保留相同的样本、基因和顺序。如果一轮删了异常样本，另一轮又按不同门槛过滤表达量，最后就无法区分变化来自分块，还是来自输入。含缺失值或零方差基因时，先看 goodSamplesGenes() 返回的标记，记录实际删除的行列，再保存共同输入。

建一次单块网络，再建一次多块网络

fit_net <- function(block_size) {
  blockwiseModules(datExpr, power = 6, networkType = "signed", TOMType = "signed",
                   maxBlockSize = block_size, minModuleSize = 15,
                   mergeCutHeight = 0.15, numericLabels = TRUE,
                   randomSeed = 17, verbose = 0)
}
a <- fit_net(250)
b <- fit_net(80)
lengths(a$blockGenes)
lengths(b$blockGenes)

本次运行中，第一次的 220 个基因都在一块，第二次形成了 78、73、69 个基因的三个块。给定上限和实际块大小是两个概念，blockGenes 能直接告诉你这次到底怎样切了块。块内先依据邻接与拓扑重叠识别模块，随后可以依据模块特征基因之间的相似性合并模块。两个基因分在不同初始块时，它们不会经历同样的块内聚类过程，因此有必要检查分块前后的对应关系。blockwiseModules 函数手册

这里固定软阈值 6 是为了让对照只改变分块。实际研究应结合数据质量、网络类型、软阈值诊断及平均连接度选定参数。若软阈值、网络类型和分块都想评估，先给每项变化各保留一组可比结果，否则一次重跑改变三个条件，很难解释究竟是哪一个影响了目标模块。

用一张交叉表找出分走的基因

A <- setNames(a$colors, colnames(datExpr))
B <- setNames(b$colors, colnames(datExpr))
stopifnot(!anyDuplicated(names(A)), setequal(names(A), names(B)))
tab <- table(A, B[names(A)])
print(tab)

按行读这张表，A2 原有 61 个基因，其中 60 个进入 B2，另一个进入 B3。按列读，B3 有 48 个基因，其中 47 个来自 A3，一个来自 A2。这两种读法分别回答“原名单保留多少”和“新名单多了谁”。模块 0 表示未分配，表里应当保留这一行和这一列，这样才能看出哪些基因退出或进入了正式模块。只有在寻找正式模块的最佳匹配时才把 0 排除。

颜色和编号由每次运行分配，碰巧相同的编号不能作为代码中的配对依据。可以对一个原模块试遍新网络中的所有正式模块，再按交集除以并集计算 Jaccard 系数。它同时考虑了原模块丢掉的基因和候选模块多出来的基因。如果一个很小的模块被完全并入一个大模块，只看原成员保留率会得到百分之百，此时 Jaccard 和列方向的交叉表更能看出合并。

x <- names(A)[A == 2]
targets <- sort(unique(B[B != 0]))
j <- vapply(targets, function(k) {
  y <- names(B)[B == k]
  length(intersect(x, y)) / length(union(x, y))
}, numeric(1))
best <- targets[which.max(j)]
y <- names(B)[B == best]
c(best_match = best, shared = length(intersect(x, y)), J = max(j))
setdiff(x, y)  # 原模块中这次没有进入最佳匹配模块的基因
setdiff(y, x)  # 新模块中多出来的基因

本例 A2 与 B2 的系数是 60/61，约 0.984。A3 与 B3 则是 47/48，约 0.979。运行第一条 setdiff() 会得到 g176。回看模拟设计，前三组共同信号只覆盖前 175 列，g176 恰好是第一列独立噪声。有限样本中，一个噪声变量也可能偶然跟某个模块共享表达趋势。这个已知答案解释了为什么即使模块整体很稳定，具体候选基因也值得单独核查。

真实数据里没有这张答案表。如果准备验证的恰好是移动的基因，可以画出它和两个候选模块特征基因在各样本中的表达关系，再按批次着色，查看趋势是否由一个异常样本拉动。若交叉表某一行的成员大量分散到两列，就把拆分的两个子模块都列出来，保留一对多的对应。本例在 R 4.5.2、WGCNA 1.73 下实算，另用 R 4.3.3、WGCNA 1.71 复算，得到了同一张交叉表。

最后把证据分层记录。这个对照支持“在同一矩阵上，所选基因模块对这一次分块调整较稳定”。它仍然没有检验在另一批受试者中能否复现。涉及独立队列时，应重新核对表达尺度、共同基因和样本结构，再做模块保留性分析。将这两项结果分别写入分析记录，读者就能分清参数敏感性和跨队列重复性。

哪一种重叠能支持下一步实验？

假设研究者要从一个旧模块挑候选基因，而新分析把这个模块的二十个成员全部收入一个有两百个基因的大模块。从旧名单看，保留率是百分之百，从新名单看，旧成员只占十分之一，交集除以并集也是十分之一。这种情形很适合在报告中加一列“新模块中原成员比例”。若研究的问题是原来的候选能不能找回来，第一列有意义。若想沿用旧模块的功能命名来解释新模块中的全部两百个基因，另外一百八十个成员就必须解释。单个高百分数无法同时回答这两个问题。

换一种情形，旧模块的四十个成员等分给两个新模块，各自都没有混入其他基因。此时为旧模块只保留一个“最佳匹配”会遗漏一半成员，即使两个候选的系数并列。可先按照旧模块成员落入各列的比例排序，并完整保留第一名与第二名，再查看拆成两群的成员在同一批样本上的表达趋势。若一个小群主要与某个组织成分或技术批次同步，就带着分组变量去复核；若两个子群有稳定但各自不同的趋势，可以分别解释。基因本体富集所用背景也要与这次实际进入建网的基因集合一致。

overlap <- function(old, new) {
  common <- length(intersect(old, new))
  c(old_kept = common / length(old), new_from_old = common / length(new),
    Jaccard = common / length(union(old, new)))
}
overlap(1:20, 1:200)
rbind(first = overlap(1:40, 1:20), second = overlap(1:40, 21:40))

对一个按计划要验证的基因，再多查一层

有了变化名单后，可以检查一个基因在模块中的归属强度。模块特征基因概括一组成员在各样本上的共同变化。把候选基因表达与两个候选模块的特征基因分别画成散点，并始终用同一套样本编号和颜色。如果绝大多数样本之间几乎没有趋势，只有一个高表达样本在图的右上角，相关系数会对这个点很敏感。先查该样本的质量指标、来源和测序批次，再计算保留和去掉它时的数值，记录差异。这样能指出相关性依赖哪份观测，也避免看到一个碍眼的点就把它从主分析中删除。

还可以问一个更严格的问题。特征基因是由模块成员共同估计出来的，拿成员与包含自身的模块总趋势计算相关会有一定自我贡献，小模块尤为值得留意。针对准备验证的候选，可以在同一份样本上临时去掉它，用余下成员重新求第一主成分，再比较候选与这个留一趋势的关系。主成分的正负号可以整体翻转，跨两次计算解释方向前先用共同成员的平均表达或同一参考特征对齐方向。若候选只在包含自己的小模块中显得很强，报告中应把这个依赖关系讲清楚。

对于这份模拟数据，我们知道第一个噪声基因的来源。换到真实基因时，把三种证据放在同一行会更容易作决定。第一项记录跨分块运行时成员是否稳定，第二项记录散点是否由少数样本主导，第三项记录独立样本能否验证相似的表达关系。即使候选在第一项发生移动，它也未必应当剔除，尤其当两个候选模块本来就相关时。把有生物学依据但归属模糊的基因单列，并保留原始图和编号，会比把它强行归到某种颜色更方便后续实验设计。

下载代码和示例数据，在 Linux 上复现

下载本篇 R 代码（ZIP，5.0 KiB） · 下载配套示例数据（ZIP，80.6 KiB）

将两个 ZIP 解压到同一位置，进入 041-wgcna 目录，运行 Rscript --vanilla run.R。打开 results/metrics.tsv 查看结果，代码包中的 expected/metrics.tsv 可供逐行核对。

本次从下载包复现的步骤是“在同一表达矩阵上运行两次分块建网”。Linux x86_64、R 4.3.3、WGCNA 1.71 的输出为建网基因 220，最佳匹配共同基因 60，Jaccard 系数 0.98360656。

下载版另存两次分块建网的交叉表。找到共有 60 个基因的单元格，沿行、列分别检查剩余基因去了哪里。并集为 61，所以 Jaccard 为 60/61。把整张交叉表留下来，还能区分一个基因迁移与整个模块分裂。

实跑后的经验总结

在固定的模拟矩阵中，大多数候选基因在两种分块设置下保持归属，真正改变模块的基因可以按名字追到。最有用的记录是保存输入基因顺序和完整成员表，再按两次分析的交集和并集匹配模块。这样既能看到整体稳定，也能找出少数会影响实验候选的边界成员；仅保存最后的模块颜色就做不到这一点。

当基因名单和比较规则都确定，后面还要在多个队列上重复建网时，可以把同一份 R 脚本放到云生信的生物信息学分析服务器上运行。查看云生信服务器租赁