Scientify输入数据
WGCNA 的结果主要由输入矩阵决定。官方 FAQ 不建议在少于 15 个样本的数据上做 WGCNA，尽量有 20 个以上；输入是行为样本、列为基因、已经标准化并取对数或做过方差稳定变换的矩阵。
# R >= 4.3；WGCNA 1.74 在 CRAN，依赖的 impute、preprocessCore、GO.db 在 Bioconductor
install.packages("BiocManager")
BiocManager::install(c("impute", "preprocessCore", "GO.db", "AnnotationDbi",
"DESeq2", "org.Hs.eg.db"))
install.packages("WGCNA")# 先加载 DESeq2 等 Bioconductor 包，最后加载 WGCNA，避免 cor 被 S4Vectors 覆盖
suppressMessages({library(DESeq2); library(org.Hs.eg.db); library(WGCNA)})
options(stringsAsFactors = FALSE, timeout = 600)
acc <- "GSE130970" # 78 例 NAFLD 肝活检 RNA-seq，带纤维化分期等临床评分
dir.create("geo", showWarnings = FALSE)
f_cnt <- file.path("geo", paste0(acc, "_raw_counts.tsv.gz"))
f_mat <- file.path("geo", paste0(acc, "_series_matrix.txt.gz"))
if (!file.exists(f_cnt)) download.file(paste0(
"https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=", acc,
"&format=file&file=", acc, "_raw_counts_GRCh38.p13_NCBI.tsv.gz"), f_cnt, mode = "wb")
if (!file.exists(f_mat)) download.file(paste0(
"https://ftp.ncbi.nlm.nih.gov/geo/series/GSE130nnn/", acc, "/matrix/", acc,
"_series_matrix.txt.gz"), f_mat, mode = "wb")
cnt <- as.matrix(read.delim(f_cnt, row.names = 1, check.names = FALSE)) # 39376 x 78，行名是 Entrez ID
# 1. 去掉低表达基因：在 90% 以上样本中 count < 10 的基因（官方 FAQ 的例子）
keep <- rowSums(cnt >= 10) >= 0.1 * ncol(cnt)
sym <- mapIds(org.Hs.eg.db, rownames(cnt), "SYMBOL", "ENTREZID")
cnt <- cnt[keep & !is.na(sym), ]
rownames(cnt) <- make.unique(sym[rownames(cnt)]) # 23337 个基因
# 2. 方差稳定变换；不能把 counts 或未取对数的 FPKM 直接送进 WGCNA
dds <- DESeqDataSetFromMatrix(cnt, data.frame(row.names = colnames(cnt), x = rep(1, ncol(cnt))), ~1)
expr <- assay(vst(dds, blind = TRUE))
# 3. 按 MAD 取前 5000 个基因，转置成 行 = 样本、列 = 基因
datExpr <- t(expr[order(apply(expr, 1, mad), decreasing = TRUE)[1:5000], ])
dim(datExpr) # 78 5000
# 临床性状：从 Series Matrix 的 characteristics 行解析
sm <- readLines(gzfile(f_mat))
row <- function(p) lapply(sm[startsWith(sm, p)], function(x) gsub('"', "", strsplit(x, "\t")[[1]][-1]))
ch <- row("!Sample_characteristics_ch1")
tr <- as.data.frame(lapply(ch, function(v) sub("^[^:]+: ", "", v)))
names(tr) <- sapply(ch, function(v) sub(":.*", "", v[1]))
rownames(tr) <- row("!Sample_geo_accession")[[1]]
tr <- tr[rownames(datExpr), ]
traits <- data.frame(fibrosis = as.numeric(tr$`fibrosis stage`),
NAS = as.numeric(tr$`nafld activity score`),
steatosis = as.numeric(tr$`steatosis grade`),
ballooning = as.numeric(tr$`cytological ballooning grade`),
inflammation = as.numeric(tr$`lobular inflammation grade`),
age = as.numeric(tr$`age at biopsy`),
male = as.numeric(tr$Sex == "M"), row.names = rownames(tr))直接用 counts 的后果
本页把 GSE130970 的 raw counts 不做变换直接按 MAD 取前 5000 个基因：与 vst 后选出的 5000 个只重叠 1296 个，被选中的基因 count 中位数为 1456，vst 选出的为 91.5，筛选被高表达基因主导。signed 网络的 R² 到 β = 30 仍只有 0.81。整数矩阵直接给 blockwiseModules 会报错 REAL() can only be applied to a 'numeric', not a 'integer'。
按性状相关筛基因的后果
按与纤维化分期的相关性取前 2000 个基因建网：signed 网络的 R² 到 β = 22 才超过 0.85；grey 只剩 10 个基因，几乎所有基因都被分进模块；最强的模块-纤维化相关从 0.58 升到 0.73。这个 0.73 来自筛选本身，不能作为结果。先取差异基因再做 WGCNA 的流程有同样问题。
| 环节 | 做法 | 依据或实测 |
|---|---|---|
| 样本数 | 少于 15 个不做；20 个以上；按性别、组织等类别分别建网（consensus 分析）时每类约 30 个以上 | 官方 FAQ（2020-06-10 更新） |
| 低表达过滤（RNA-seq） | 去掉在 90% 以上样本中 count < 10 的基因 | FAQ 给出的例子；GSE130970 从 39376 个降到 23337 个 |
| 变换 | counts 用 DESeq2 的 vst 或 varianceStabilizingTransformation；FPKM、TPM 或标准化 counts 用 log2(x + 1) | FAQ；本页 vst 与 log2 CPM 的逐样本相关中位数 0.9985 |
| 基因数 | 按 MAD 或方差取前 5000 个左右（常见 3000–10000）；筛掉的是低变异基因 | FAQ：均值和方差筛选结果相近；官方肝脏教程用 3600 个 |
| 不要按差异表达筛选 | WGCNA 是无监督方法，用差异基因建网会只剩一个或几个高度相关的模块，无标度拟合失效 | FAQ；本页实测见下文 |
| 批次 | 先检查批次；类别型批次用 sva::ComBat，连续型技术变量用线性回归去除 | FAQ |
离群样本
官方教程用样本聚类树手动定切割高度（雌鼠肝脏数据切在 15，去掉 1 个样本）。样本较多、树上没有明显孤立分支时，用标准化连接度 Z.k 判断更客观，常用阈值是 Z.k < -2.5。
GSE130970 的聚类树最高合并高度为 84.4，没有单个孤立样本；Z.k 标出 GSM3758009（-6.14）和 GSM3758047 两个样本。去掉后剩 76 例。本数据中保留这两个样本时结果相近（β 同为 7，14 个模块，grey 1334 个基因），离群样本的影响在小样本数据中更大。
聚类树上出现两大分支时，先把临床和技术变量（批次、性别、测序日期）画在树下方对照，判断分支来自生物分组还是技术因素，再决定是去除样本、校正还是分组建网。
忘记转置是不报错的错误：基因在行、样本在列的矩阵送进 goodSamplesGenes 同样返回 TRUE，之后会把样本当成基因建网。进入分析前检查 dim(datExpr)，行数应等于样本数。
gsg <- goodSamplesGenes(datExpr, verbose = 0); gsg$allOK # 有缺失或零方差基因时为 FALSE
sampleTree <- hclust(dist(datExpr), method = "average")
plot(sampleTree, cex = 0.6, main = "Sample clustering") # 先看有没有单独挂在外面的样本
# 标准化连接度 Z.k：样本网络中连接度过低的样本视为离群（常用阈值 -2.5）
A <- adjacency(t(datExpr), type = "distance")
Zk <- as.numeric(scale(colSums(A) - 1))
rownames(datExpr)[Zk < -2.5] # GSE130970: GSM3758009 GSM3758047
datExpr <- datExpr[Zk >= -2.5, ]; traits <- traits[rownames(datExpr), ]
nSamples <- nrow(datExpr) # 76软阈值
pickSoftThreshold 返回第一个使无标度拟合 R²（signed R²）超过 RsquaredCut 的 β，默认 RsquaredCut = 0.85、默认 networkType = "unsigned"。选 β 时的 networkType 和相关系数要与后面 blockwiseModules 完全一致。
达不到阈值时，官方 FAQ（2017 年 12 月改为更保守的数值）给出的经验 β：unsigned 与 signed hybrid 网络，样本数小于 20 取 9，20–30 取 8，30–40 取 7，40 以上取 6；signed 网络分别为 18、16、14、12。使用前提是已经排除批次、离群样本和按差异表达筛基因这些原因。本页在 76 例上用 FAQ 的 β = 12 建 signed 网络：13 个模块，grey 1661 个基因，与纤维化最强的相关为 0.53（β = 7 时为 0.58）。
FAQ 认为“合理”的 β 范围：unsigned 和 signed hybrid 小于 15，signed 小于 30。在这个范围内 R² 都到不了 0.8、平均连接度仍在几百以上时，数据中通常有一个把部分样本整体拉开的因素。本页模拟这种情况（一半样本的 40% 基因整体加 1.5）：R² 到 β = 30 只有 0.68，β = 30 时平均连接度仍为 26，β = 12 时出现一个 2088 个基因的大模块。
signed 还是 unsigned：FAQ 推荐 signed 或 signed hybrid。unsigned 网络把正相关和负相关的基因放进同一模块，模块特征基因的方向就不再对应模块中每个基因的方向。相关系数推荐 bicor（双权中值相关），必须加 maxPOutliers = 0.05 或 0.10，否则当基因表达受二分类变量（疾病状态、基因型）强烈影响时，bicor 会把一组样本当作离群值。
allowWGCNAThreads(4) # RStudio 中并行出错时改用 disableWGCNAThreads()
powers <- c(1:10, seq(12, 20, 2))
sft <- pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed",
corFnc = bicor, corOptions = list(maxPOutliers = 0.05),
RsquaredCut = 0.85, verbose = 0)
fi <- sft$fitIndices
data.frame(power = fi$Power, R2 = round(-sign(fi$slope) * fi$SFT.R.sq, 3), mean.k = round(fi$mean.k., 1))
beta <- sft$powerEstimate # GSE130970（76 例）: 7
# 达不到阈值时用官方 FAQ 的经验值（2017-12 更新版）
if (is.na(beta)) {
n <- nSamples; signed <- TRUE
beta <- if (n < 20) 9 else if (n < 30) 8 else if (n < 40) 7 else 6
if (signed) beta <- beta * 2
}
par(mfrow = c(1, 2))
plot(fi$Power, -sign(fi$slope) * fi$SFT.R.sq, type = "n", xlab = "power", ylab = "signed R^2")
text(fi$Power, -sign(fi$slope) * fi$SFT.R.sq, fi$Power, col = "red"); abline(h = 0.85, col = "red")
plot(fi$Power, fi$mean.k., type = "n", xlab = "power", ylab = "mean connectivity")
text(fi$Power, fi$mean.k., fi$Power, col = "red")R² 取 0.8、0.85 还是 0.9
0.85 是 pickSoftThreshold 的默认值；官方教程图上画的是 0.90 线，在雌鼠肝脏数据上选 β = 6（R² 0.902）；FAQ 判断数据是否有问题时用 0.8。选择原则：取 R² 曲线进入平台的第一个 β，不取更高的 β。本页数据在 β = 7 后 R² 基本平稳，β 越高平均连接度越低，到 β = 20 时只有 2.0，网络过于稀疏。
R² 后段的波动
平均连接度降到 1 以下后，拟合只靠少数基因，R² 会大幅跳动。本页 78 例 unsigned-pearson 的 R² 在 β = 22 为 0.918、β = 24 骤降到 0.340、β = 26 又回到 0.937。只看曲线前段第一个平台。
斜率不必接近 -1
有的帖子要求 slope 接近 -1。signed 网络的斜率通常更陡：本页 signed + bicor 在 β = 7 时 slope 为 -2.67。官方函数和教程都只用 R² 选 β。
四种组合的 β 不同
同一数据（78 例）的 powerEstimate：unsigned + pearson 6，unsigned + bicor 7，signed + pearson 9，signed + bicor 8。signed 网络的邻接是 (0.5 + 0.5 × cor)^β，同样的 β 下比 unsigned 稠密，需要更大的 β。
| β | signed R² | 平均连接度 |
|---|---|---|
| 5 | 0.754 | 238.7 |
| 6 | 0.819 | 143.6 |
| 7 | 0.865 | 89.2 |
| 8 | 0.892 | 57.2 |
| 10 | 0.878 | 25.7 |
| 14 | 0.916 | 7.2 |
| 20 | 0.874 | 2.0 |
模块识别
一步法 blockwiseModules 的默认值与官方教程写法并不一致。下面列出 WGCNA 1.74 的默认值、教程用值和本页实测的影响，基线为 76 例、5000 个基因、signed + bicor、β = 7（14 个模块，grey 1369 个基因）。
本页复现官方教程 I：雌鼠肝脏 3600 个探针、134 个样本，β = 6，TOMType = "unsigned"、minModuleSize = 30、mergeCutHeight = 0.25，得到 18 个模块和 99 个 grey 基因，各模块大小与教程表格逐一相同，单线程用时 7.9 秒。
β 也影响模块数。同一数据 β = 2、3、4、7 时分别得到 8、8、12、14 个模块。模块数没有“正确值”，论文中报告所用参数并做敏感性分析，比反复调参追求某个模块与性状高度相关更站得住。
net <- blockwiseModules(datExpr, power = beta,
networkType = "signed", TOMType = "signed", # networkType 要与 pickSoftThreshold 一致
corType = "bicor", maxPOutliers = 0.05,
maxBlockSize = 6000, # 大于基因数，保证只有一个块
minModuleSize = 30, deepSplit = 2, mergeCutHeight = 0.25,
pamRespectsDendro = FALSE, numericLabels = TRUE, verbose = 3)
moduleColors <- labels2colors(net$colors)
table(moduleColors) # GSE130970: 14 个模块，grey 1369 个基因
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], "Module",
dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)| 参数 | 1.74 默认值 | 常用值 | 含义与实测影响 |
|---|---|---|---|
| networkType | "unsigned" | "signed" | 必须与 pickSoftThreshold 一致。把 signed 选出的 β（14、20）用于 unsigned 网络：grey 从 1903 个增至 4121、4659 个（共 5000） |
| TOMType | "signed" | 教程写 "unsigned" | 对 unsigned 网络，signed TOM 保留相关符号再取绝对值。实测与 "unsigned" 几乎相同：雌鼠肝脏数据 18 个模块完全一致；GSE130970 的 grey 1903 与 1903、最大模块 623 与 622 |
| corType / maxPOutliers | "pearson" / 1 | "bicor" / 0.05 | bicor 不设 maxPOutliers 时，二分类驱动的表达模式会被当作离群值 |
| minModuleSize | min(20, 基因数/2) | 30 | 20/30/50/100 → 16/14/11/9 个模块 |
| deepSplit | 2 | 2 | 0–4 → 11/12/14/19/20 个模块；越大切得越细，grey 越少 |
| mergeCutHeight | 0.15 | 0.25 | 特征基因相关 > 1 − mergeCutHeight 的模块合并。0.15/0.25/0.35 → 14/14/12 个模块 |
| maxBlockSize | 5000 | 大于基因数 | 基因数超过它就分块，见下一节 |
| pamRespectsDendro | TRUE | 教程写 FALSE | 本数据中 grey 1360 与 1369，差别很小 |
| randomSeed | 54321 | 保持默认 | 分块的预聚类使用随机数；R 3.6.0 改变了随机数生成方式，复现旧结果用 RNGkind("Mersenne-Twister", "Inversion", "Rounding") |
内存
TOM 是基因数 × 基因数的矩阵，内存随基因数平方增长。基因数超过 maxBlockSize（默认 5000）时，blockwiseModules 先用 k-means 式预聚类把基因分块，每块单独建网，最后合并相近模块。分块在跨块基因的模块归属上会出错，官方教程建议尽量用一个块。
分块对结果的影响：15000 个基因分 4 块与 1 块相比，模块划分的调整兰德指数（ARI）为 0.42；单块的 27 个模块中，有 15 个在分块结果里留在同一模块的基因不到 70%；主分析中与纤维化相关的 yellow 模块基因，在 15000 个基因的单块结果中主要落在 2 个模块（190 和 128 个），在分块结果中分散到 3 个模块（123、92 和 70 个）。官方教程在 3600 个探针上比较分块与单块，结论是“非常相似”，这个结论在上万个基因时不成立。
官方教程给出的内存经验：4 GB 内存约 8000–10000 个探针，16 GB 约 20000 个，32 GB 约 30000 个。本页 16 GB 机器上 20000 个基因单块完成、峰值 6.5 GB，约 22500 个基因单块崩溃（同时有其他程序占用内存）。maxBlockSize 的硬上限是 sqrt(2^31) ≈ 46340。
内存不够时优先减少基因数（MAD 前 5000–10000 个通常已覆盖主要的共表达结构），让所有基因在一个块中；确实需要全部基因时，换大内存机器。分块结果的模块边界要在论文中说明。
| 基因数 | maxBlockSize | 块数 | blockwiseModules 用时 | 进程峰值内存 | 结果 |
|---|---|---|---|---|---|
| 5000 | 6000 | 1 | 5.9 秒 | 1.6 GB | 14 个模块 |
| 10000 | 10000 | 1 | 18.3 秒 | 2.5 GB | 21 个模块 |
| 15000 | 15000 | 1 | 72 秒 | 3.9 GB | 27 个模块，grey 3136 |
| 15000 | 5000 | 4（4996/4819/2739/2446） | 29.8 秒 | 1.6 GB | 29 个模块，grey 3864 |
| 20000 | 20000 | 1 | 196 秒 | 6.5 GB | 31 个模块 |
| 22490 | 23337 | 1 | — | 崩溃时 5.5–7.1 GB | TOM 计算完成后崩溃（segfault / bad binding access） |
模块-性状
每个模块用模块特征基因（ME，模块表达矩阵的第一主成分）代表，再与临床性状做相关。热图每格是相关系数和 P 值，红色为正相关，蓝色为负相关。
MEs <- orderMEs(moduleEigengenes(datExpr, moduleColors)$eigengenes)
moduleTraitCor <- cor(MEs, traits, use = "p")
moduleTraitP <- corPvalueStudent(moduleTraitCor, nSamples)
textMatrix <- paste0(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitP, 1), ")")
par(mar = c(6, 9, 3, 3))
labeledHeatmap(moduleTraitCor, xLabels = names(traits), yLabels = names(MEs), ySymbols = names(MEs),
colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix,
setStdMargins = FALSE, cex.text = 0.6, zlim = c(-1, 1), main = "Module-trait relationships")
# 二分类或等级性状也可用 bicor，但要关掉对性状的稳健处理：
# bicor(MEs, traits, maxPOutliers = 0.05, robustY = FALSE)负相关模块同样有效
turquoise 与 NAS 的相关为 -0.67，含义是该模块整体随病变加重而下调。按相关系数的绝对值挑模块。
把混杂变量放进热图
pink 与 tan 两个模块分别与性别相关 -0.76 和 0.77（P = 2.8e-16），是性别模块。性状表中加入性别、年龄、批次，能识别出这类模块，避免把它们解释成疾病相关。
多重比较
14 个模块 × 7 个性状共 98 个检验，本页 P < 0.05 的有 33 个，Bonferroni 校正（0.05/98）后剩 17 个。挑模块时同时看相关强度和校正后的显著性，并预先确定关注的主要性状。
grey 行不解释
grey 是未分入任何模块的基因集合，它的 ME 与性状的相关没有生物学含义，作图时可以去掉。
二分类与等级性状
性别、分组等 0/1 变量和分期这类等级变量可以直接用 Pearson 相关。改用 bicor 时加 robustY = FALSE；大多数样本取同一个值时，bicor 会提示 zero MAD in variable 'y' 并对该列退回 Pearson。本页 Pearson 与 bicor(robustY = FALSE) 的结果最大相差 0.15。
| 模块 | 纤维化分期 | NAS 评分 | 气球样变 | 男性 |
|---|---|---|---|---|
| yellow（405 个基因） | 0.58（P = 3.1e-8） | 0.69（P = 8.4e-12） | 0.57 | -0.12 |
| turquoise（798 个基因） | -0.52 | -0.67 | -0.54 | -0.07 |
| pink（164 个基因） | -0.02 | -0.16 | -0.11 | -0.76 |
| tan（98 个基因） | -0.26 | -0.03 | -0.06 | 0.77 |
hub 基因
GS（基因显著性）是基因表达与性状的相关，MM（模块隶属度，即 kME）是基因表达与模块特征基因的相关，kWithin 是模块内连接度。在与性状显著相关的模块里，GS 与 MM 也呈显著正相关时，模块内的核心基因更可能与性状有关。
|GS| > 0.2 的显著性取决于样本量。双侧 P < 0.05 所需的 |GS|：20 个样本 0.444，30 个 0.361，50 个 0.279，76 个 0.226，100 个 0.197；样本数达到 97 时 0.2 才显著。样本量小于 97 时，用 |GS| > 0.2 筛出的基因可能与性状没有显著相关，阈值宜按样本量换算，或直接要求 GS 的 P 值显著。
WGCNA 作者 Peter Langfelder 在 Bioconductor 支持站的回答：kME 与模块内连接度的排序通常非常接近，作者团队常用 kME，因为它易于计算并有对应的 P 值；排序只作为后续验证的优先级参考，多个基因 kME 都很高且接近时都可视为 hub，第 3 名与第 7 名没有实质区别。本页 yellow 模块按 kWithin 排序的前 10 个为 TNFRSF12A、YWHAH、KPNA2、ARPC5、THBS2、ANXA2、TIGAR、ANXA2P2、LGALS3、CCL20，与 |MM| > 0.8 且 |GS| > 0.2 的 12 个基因基本重合。
yellow 模块内 |MM| 与 |GS| 的相关为 0.715（P = 1.5e-64）。这一相关高，说明模块的核心基因同时也是与纤维化相关最强的基因。相关接近 0 时，即使 ME 与性状相关显著，模块内的 hub 也不一定与性状有关。
Langfelder 等 2013 年在 PLoS One 比较了 hub 基因选择与标准 meta 分析：共识模块中的 hub 基因在生物学解释上更有用，验证成功率上标准 meta 分析不差于网络方法。目标是找生物标志物时，把 GS 的显著性放在 MM 前面。
module <- "yellow"; trait <- "fibrosis"
MM <- cor(datExpr, MEs[, paste0("ME", module)], use = "p")[, 1] # 模块隶属度（kME）
GS <- cor(datExpr, traits[[trait]], use = "p")[, 1] # 基因显著性
inMod <- moduleColors == module
verboseScatterplot(abs(MM[inMod]), abs(GS[inMod]), xlab = paste("MM in", module),
ylab = paste("GS for", trait), col = module, abline = TRUE)
abline(v = 0.8, h = 0.2, lty = 2)
# 样本量 n 下 |GS| 的显著性门槛（双侧 P < 0.05）
tq <- qt(0.975, nSamples - 2); sqrt(tq^2 / (tq^2 + nSamples - 2)) # n = 76 时为 0.226
# 模块内连接度 kWithin
adj <- adjacency(datExpr, power = beta, type = "signed",
corFnc = "bicor", corOptions = "maxPOutliers = 0.05")
kIM <- intramodularConnectivity(adj, moduleColors)
hub <- data.frame(gene = colnames(datExpr), MM = MM, GS = GS, kWithin = kIM$kWithin)[inMod, ]
hub <- hub[order(-hub$kWithin), ]
head(hub, 10)
subset(hub, abs(MM) > 0.8 & abs(GS) > 0.2) # GSE130970: 12 个
cor(hub$kWithin, hub$MM, method = "spearman") # 0.95：两种排序几乎一致| 筛选标准 | yellow 模块（405 个基因）中入选数 | 说明 |
|---|---|---|
| |MM| > 0.8 且 |GS| > 0.2 | 12 | 国内论文最常用；WGCNA 官方教程和 FAQ 中没有这个阈值 |
| |MM| > 0.9 且 |GS| > 0.2 | 0 | 阈值稍严就可能一个都没有 |
| |MM| > 0.8 且 |GS| > 0.4 | 8 | |
| |MM| > 0.8 且 |GS| > 0.5 | 5 | |
| kWithin 前 30 | 30，包含上面全部 12 个 | kWithin 与 MM 的 Spearman 相关 0.95 |
报错与坑
下表的报错原文多数在 WGCNA 1.74 上复现，其余来自官方 FAQ。
# 报错：unused arguments (weights.x = NULL, weights.y = NULL, cosine = FALSE)
find("cor") # 若 IRanges / S4Vectors 排在 WGCNA 前面，就会出错
cor <- WGCNA::cor # 临时指定，跑完 blockwiseModules 后再恢复
net <- blockwiseModules(datExpr, power = beta, networkType = "signed", numericLabels = TRUE)
cor <- stats::cor
# 或者：重开 R，先加载 DESeq2 等包，最后 library(WGCNA)
# 报错：REAL() can only be applied to a 'numeric', not a 'integer'
storage.mode(datExpr) <- "double" # 整数矩阵（如直接用 counts）会触发；但 counts 本身也不该直接用| 现象或报错原文 | 原因 | 处理 |
|---|---|---|
| unused arguments (weights.x = NULL, weights.y = NULL, cosine = FALSE) | WGCNA 自带的 cor 被其他包的 cor 覆盖。本页实测：先 library(WGCNA) 再 library(DESeq2) 时，IRanges / S4Vectors 的 cor 泛型排在前面，报错；顺序反过来不报错 | 最后加载 WGCNA；或在 blockwiseModules 前 cor <- WGCNA::cor，之后改回 stats::cor |
| REAL() can only be applied to a 'numeric', not a 'integer' | 输入是整数矩阵（通常是直接用了 counts） | 先做 vst 或 log 变换；只想转类型用 storage.mode(x) <- "double" |
| could not find function 或 GOenrichmentAnalysis 报 deprecated | WGCNA 没有成功加载；或调用了 1.74 已移除的 GOenrichmentAnalysis（现在只返回提示信息） | library(WGCNA) 看报错；富集分析改用 clusterProfiler 或作者的 anRichment |
| pickSoftThreshold 在 RStudio 中并行出错 | 第三方 GUI 中的多线程问题 | disableWGCNAThreads() 后重跑 |
| thread 0 could not be started successfully. Error code: 11 | 集群每个任务只分配 1 个核 | disableWGCNAThreads() |
| malloc: *** mmap(size=...) failed ... can't allocate region（Mac） | FAQ 认为是无害信息 | 可忽略；若随后 R 崩溃，按内存不足处理 |
| R 崩溃（segfault、bad binding access） | 单块基因太多。本页 16 GB 机器上约 22500 个基因单块时在 TOM 计算后崩溃 | 减少基因数或换大内存机器；分块要说明代价 |
| 探针名前多了 X | 矩阵转成 data.frame 后，数字开头的列名被加上 X | 全程用 matrix，或读入时 check.names = FALSE |
| grey 占一半以上 | β 与 networkType 不匹配；β 过高网络太稀疏；minModuleSize 过大；样本少、噪声大 | 检查 pickSoftThreshold 与 blockwiseModules 的 networkType 一致；本页正确设置时 grey 约 27% |
| 只有 1–2 个模块或一个巨大模块 | 按差异基因或性状筛选了基因；批次或组织差异主导表达；β 太低 | 改用 MAD/方差筛选；画聚类树查批次并校正；按 FAQ 判断是否改用经验 β |
审稿
WGCNA 论文被质疑最多的是样本量和 hub 基因没有验证。验证可以完全用公开数据完成：用 modulePreservation 检验模块在独立数据集中是否保留，再看 hub 基因与同一性状的相关。
# 外部数据（同样处理成 vst 矩阵，行 = 样本，列 = 基因）：GSE135251，216 例
common <- intersect(colnames(datExpr), colnames(datExpr2))
mp <- modulePreservation(
list(ref = list(data = datExpr[, common]), test = list(data = datExpr2[, common])),
list(ref = setNames(moduleColors, colnames(datExpr))[common]),
referenceNetworks = 1, nPermutations = 50, networkType = "signed",
corFnc = "bicor", randomSeed = 1, maxGoldModuleSize = 300, maxModuleSize = 1000, verbose = 0)
z <- mp$preservation$Z$ref.ref$inColumnsAlsoPresentIn.test
z[, c("moduleSize", "Zsummary.pres")] # Zsummary > 10 强保留，2-10 中等，< 2 不保留| 质疑 | 应对 |
|---|---|
| 样本量小 | 写明样本数与 FAQ 的建议（≥ 15，最好 ≥ 20）；β 达不到阈值时说明采用 FAQ 经验值；把 GS 阈值按样本量换算成显著性门槛 |
| hub 基因没有验证 | 在独立队列中做模块保留分析并检验 hub 基因与性状的相关。本页用 GSE135251（216 例）验证：yellow 模块 Zsummary 20.4（> 10 为强保留）；11 个候选 hub 中 10 个与纤维化分期显著相关（Spearman rho 0.21–0.60），ARPC5 不显著（P = 0.096）；用时 157 秒 |
| 先取差异基因再做 WGCNA | 改为按 MAD 或方差筛选全部表达基因；差异基因与模块基因取交集放在 WGCNA 之后 |
| β 的选择依据不清 | 给出 R² 与平均连接度两张图，写明 RsquaredCut、networkType 和相关系数 |
| 参数随意 | 报告 minModuleSize、deepSplit、mergeCutHeight、maxBlockSize 与是否分块；补充 deepSplit 或 β 的敏感性分析，说明关键模块是否稳定 |
| 模块与性状的相关可能来自混杂 | 热图中加入性别、年龄、批次；对关键模块做调整混杂后的回归 |
| 只有生信分析 | hub 基因用 qPCR、免疫组化或外部数据（TCGA、HPA、其他 GEO 队列）验证 |
国内经验
以下来自 CSDN、简书、腾讯云和知乎专栏的 WGCNA 教程，每条都有报错原文或源码，并在本页复现或核对。
cor 冲突的报错与临时替换
CSDN 帖子（2019）贴出 blockwiseModules 报错 Error in (new("standardGeneric", .Data = function (x, y = NULL, ...: unused arguments (weights.x = NULL, weights.y = NULL, cosine = FALSE)，用 cor <- WGCNA::cor 后运行成功，并提醒用完改回 stats::cor。本页复现到同一报错，原因是 DESeq2 依赖的 IRanges / S4Vectors 在 WGCNA 之后加载。
MAD 筛选与经验 β
生信宝典在简书的 WGCNA 教程（2018）去掉 MAD 最低的 25% 基因（且 MAD 至少 0.01），无合适 β 时按 FAQ 表格取经验值，并推荐 signed 网络和 bicor。注意同文模块-性状代码写成 if (corType == "pearsoon")，拼写错误使 Pearson 分支永远不执行，照抄时即使选 pearson 也会走 bicor。
β 与网络类型不一致
一篇 2026 年的 CSDN 结果解读帖用 networkType = "signed" 选 β，blockwiseModules 却没写 networkType（默认 unsigned）且 TOMType = "unsigned"。本页实测这种错配：signed 的 β 用在 unsigned 网络上，grey 基因占 82%–93%。
FPKM 要先取对数
腾讯云「RNA-seq 入门实战」对 FPKM 做 log2(x + 1) 后按 MAD 取前 5000 个基因，maxBlockSize 设为基因数，做法与 FAQ 一致；同文正文写 power = 16、代码是 15，照抄时以代码为准并自己重跑 pickSoftThreshold。知乎一篇高排名教程对 FPKM 只做均值过滤、不取对数，属于反例。
交给 Agent
数据下载、变换、选 β、建网、作图和外部验证的步骤固定，适合交给科学智能体执行；选数据集、定义性状和解释模块仍由你负责。
一句话指令示例：「用 GSE130970 做 WGCNA：NCBI raw counts 过滤低表达后做 vst，取 MAD 前 5000 个基因，用 Z.k < -2.5 去离群样本；signed 网络加 bicor（maxPOutliers = 0.05）选 β，单块 blockwiseModules（minModuleSize 30、mergeCutHeight 0.25）；与纤维化、NAS、性别、年龄做模块-性状热图；对最相关模块输出 GS-MM 散点图和按 kWithin 排序的 hub 表；再用 GSE135251 做 modulePreservation 和 hub 基因与纤维化的相关验证；附 deepSplit 0–4 的敏感性分析。」
- 01
数据与预处理
下载 counts 和 Series Matrix，解析性状，做低表达过滤、vst 和 MAD 筛选，输出聚类树和离群样本列表。
- 02
选 β 与建网
运行 pickSoftThreshold 并作图，确认 networkType 一致；按内存决定单块，记录用时和内存。
- 03
模块-性状与 hub
输出热图（含校正后 P 值）、GS-MM 散点图、hub 基因表和模块基因列表。
- 04
验证与敏感性分析
在独立数据集中做模块保留分析和 hub 基因相关；比较不同 deepSplit 和 β 下关键模块是否稳定。
- 05
对抗审阅
检查是否按差异基因筛选、β 与网络类型是否一致、是否发生分块、|GS| 阈值是否达到显著、热图中是否有性别等混杂模块。
- 工作区保留：原始数据、R 脚本、sessionInfo、β 选择图、模块分配表、热图和验证结果。
- 你仍需核对：性状编码是否正确、离群样本的去除是否合理、所选模块的生物学解释。
- qPCR、免疫组化等实验验证由你完成。
参考资料
- CRAN: WGCNA 1.74 — 当前版本与发布日期（2026-01-30）
- WGCNA ChangeLog — 1.74 移除 GOenrichmentAnalysis；maxBlockSize 上限；预聚类默认值变化
- WGCNA package FAQ（Langfelder & Horvath） — 样本数、基因筛选、RNA-seq 变换、signed 与 bicor、经验 β 表、多线程与常见报错
- WGCNA 官方教程（雌鼠肝脏数据） — 教程 I 的数据、一步法与分块法、maxBlockSize 内存经验、模块-性状与 GS/MM
- Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 2008 — 方法原文
- Langfelder P, Luo R, Oldham MC, Horvath S. Is my network module preserved and reproducible? PLoS Comput Biol, 2011 — modulePreservation 与 Zsummary 阈值
- Langfelder P, Mischel PS, Horvath S. When is hub gene selection better than standard meta-analysis? PLoS One, 2013 — hub 基因选择与边际分析的比较
- Bioconductor support: WGCNA hub gene selection method — Langfelder 关于 kME 与 kIM、hub 排序的回答
- Hoang SA, et al. Gene expression predicts histological severity and reveals distinct molecular profiles of nonalcoholic fatty liver disease. Sci Rep, 2019（GSE130970） — 本页实测数据集
- GSE135251（216 例 NAFLD 肝活检） — 本页外部验证数据集
- CSDN：转录组 WGCNA 包使用报错 — 经验帖：cor 冲突的报错原文与 WGCNA::cor 替换
- 简书（生信宝典）：WGCNA 分析，简单全面的最新教程 — 经验帖：MAD 筛选、经验 β、signed 与 bicor；pearsoon 拼写错误的反例
- CSDN：WGCNA 代码及结果解读 — 经验帖：β 与网络类型不一致的反例
- 腾讯云：RNA-seq 入门实战（十一）WGCNA 加权基因共表达网络分析 — 经验帖：FPKM 取 log2、MAD 前 5000、单块设置
- 知乎专栏：加权基因共表达网络分析（WGCNA） — 经验帖：FPKM 未取对数的反例