企业官网建设流程全解析

做转录组数据分析的朋友，对WGCNA这个名字应该都不陌生。这是一套基于R语言的加权基因共表达网络分析方法，全称Weighted Gene Co-expression Network Analysis。简单说，它把你手上成百上千个基因的表达量数据，按照表达模式的相关性组织成网络，再从中找到表达行为相似的基因模块，进一步和样本的性状、临床信息、处理条件等关联起来，最终锁定那些可能起关键作用的基因。这套方法在医学、农学、植物学等领域都应用得非常广泛，尤其在寻找生物标志物、候选基因、关键调控因子的场景里，基本属于必跑的分析流程。

这篇文章是系列的第一篇，我会从拿到表达矩阵之后的第一步开始，把数据预处理、样本聚类检查和软阈值选择这三个最基础也最容易卡住的环节，用R代码加逐段解读的方式完整过一遍。内容适合已经有R语言基础、能看懂数据框操作、但对WGCNA流程还不熟的朋友。如果你是完全零基础，建议先把R的基本数据结构、tidyverse常用的几个函数跑熟，再回来看这篇会更顺手。第二篇再讲模块识别、模块与性状关联、hub基因筛选这些下游分析。

1. WGCNA到底在做什么：核心思路与应用场景

1.1 从相关矩阵到共表达网络的底层逻辑

WGCNA的中心思想听起来不复杂：如果两个基因在不同样本里的表达量变化趋势高度一致，那它们很可能参与同一个生物学过程或被同一个调控机制控制。传统做法是算皮尔逊相关系数，设置一个阈值，相关系数超过阈值的基因对就算“有关系”，这是硬阈值。问题在于阈值怎么定都有点拍脑袋，而且信息丢失严重——相关系数0.79和0.81的差别可能本身没有生物学意义，却因为一条线被划成了完全不同的结果。

WGCNA改用了软阈值的策略，不把相关关系硬切掉，而是把相关系数做幂指数运算，权重按照相关系数的大小连续变化。公式是连接强度 (a_{ij} = |cor(i,j)|^\beta)，这个(\beta)就是软阈值。这样做的好处是保留了基因间关系的强度信息，网络结构更加稳定，也更符合生物系统连续渐变的特性。后面代码里pickSoftThreshold就是在帮你自动挑选合适的(\beta)值。

另外一个关键点是，WGCNA关注的不是单个基因对，而是模块（module）。模块是一组高度互联的基因集合，你可以把它理解为转录组里的“功能单元”。分析的核心产出之一，就是把几千个基因归并成几十个以内可解释的模块，再用模块特征基因（module eigengene，简称ME）来代表每个模块的表达模式，跟表型数据做关联分析。这一步能大幅降低分析的维度，把“几千个基因”变成“十来个模块”，生物学解释的可行性一下就上来了。

1.2 什么时候该用WGCNA：适用场景与前置条件

WGCNA不是万能工具，它有自己的适用边界。最理想的应用场景是：样本量在15到50个之间，每个样本都有对应的表型数据（比如疾病组和对照组、不同发育时期、不同处理浓度等），你有全转录组或全基因组的表达谱数据。样本量太小，相关性估计不稳定；样本量太大，计算时间会非常感人，不过WGCNA包对大数据集做了分块处理，几百个样本也能跑，只是后面我会提到内存和时间的代价。

还有一类场景特别适合WGCNA，那就是你不满足于只做差异表达分析，想知道哪些基因“协同变化”，想从系统层面理解转录调控。差异表达分析回答的是“哪些基因变了”，WGCNA回答的是“哪些基因一起变，这些基因和什么性状有关”，两者互补，不冲突。实操中我一般两个都做，先跑差异表达拿到候选基因列表，再用WGCNA看这些基因在共表达网络里的位置和归属。

前置条件里最容易被忽略的是数据质量。WGCNA对缺失值敏感，对极端离群样本敏感，对批次效应也很敏感。正式跑网络构建前，花点时间做数据清洗和样本聚类检查，后面会省下大量排查问题的时间。这就是接下来第二部分要详聊的内容。

2. 环境准备与数据预处理：从表达矩阵到干净数据

2.1 R包安装与版本选择

WGCNA包在CRAN上可以直接安装，但声明一下，从CRAN安装的版本足够日常使用。安装命令如下：

# 安装BiocManager（如果还没有） if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") # 安装WGCNA BiocManager::install("WGCNA") # 加载 library(WGCNA)

需要注意的一点是，WGCNA对R版本有要求，建议安装前先更新R到当前的最新稳定版，而不是用系统自带的老版本。装完之后建议跑一下:

# 开启多线程（后面计算会快很多） allowWGCNAThreads()

这条命令会调用机器上可用的CPU线程，对后续的网络构建加速非常明显。如果你在Windows上遇到线程相关的报错，可以先运行disableWGCNAThreads()确认问题，再决定是否启用。实际分析中，多线程带来的提速在大样本量时候特别明显，我曾经用24线程跑过一个200个样本的分析，耗时从单线程的40分钟缩短到了不到10分钟，所以这个细节建议不要忽略。

2.2 输入数据格式与标准化的正确姿势

WGCNA标准输入是一个表达矩阵，行是基因，列是样本，值是表达量。至于这个表达量是RPKM、FPKM还是TPM，理论上有争议，但实操中大部分人直接采用经过标准化的表达值。如果做的是芯片数据，用RMA或MASS归一化后的结果；如果是RNA-seq，建议用DESeq2的varianceStabilizingTransformation或EdgeR的CPM/logCPM转换，把数据拉到一个近似正态分布的尺度上。

数据导入之后的第一个正经步骤，是过滤掉低表达基因。低表达基因的技术噪音大，相关性不可靠，留着只会干扰网络构建。常见的过滤策略有两种：一是保留在所有样本中表达量大于某个阈值的基因，二是保留方差排前一定比例的基因。下面是我常用的代码：

# 读取表达矩阵，行为基因，列为样本 expr <- read.csv("expression_matrix.csv", row.names = 1, check.names = FALSE) # 过滤低表达基因：至少在80%样本中表达量大于1 keep <- rowSums(expr > 1) >= 0.8 * ncol(expr) expr_filtered <- expr[keep, ] # 过滤低变异基因：保留MAD（绝对中位差）前75%的基因 mads <- apply(expr_filtered, 1, mad) expr_filtered <- expr_filtered[order(mads, decreasing = TRUE)[1:ceiling(nrow(expr_filtered) * 0.75)], ] dim(expr_filtered)

最后一步用mad排序再截取前75%，是为了在保留生物学信号的同时尽量去掉噪音基因。你也可以保留更多，比如前5000个高变基因，这在芯片时代很常见。我个人更推荐从“所有表达达标的基因”出发，而不是预先砍到几千个，因为WGCNA本身就能通过模块化把基因分组，你把输入砍得太狠反而可能丢掉与表型相关的基因。

对于基因注释信息，建议在表达矩阵里保留一列基因Symbol或Entrez ID，后面做富集分析和数据库注释要用。分析过程中可能会遇到基因名重复的问题，记得先去重，否则相关性计算会出错：

# 假设第二列是基因名 merged <- aggregate(expr_filtered[, -1], by = list(gene = expr_filtered$gene), FUN = mean) rownames(merged) <- merged$gene merged <- merged[, -1]

这种按基因名求平均的去重做法，比简单保留第一个出现的基因更稳，因为它把相同基因名下的多条探针或转录本表达量合并了。

2.3 样本聚类与离群样本处理：别急着建网络

样本质量检查这一步，很多人会跳过，但我强烈建议不要跳过。样本聚类能直观地显示出是否有离群样本，尤其是技术重复、批次差异导致的异常样本。如果再配合表型信息，你还能看出聚类结果和分组是否一致。不一致才是正常的，转录组数据的样本聚类本身就能反映样本的生物学特征。但如果同一个样本明显偏离所有其他样本，通常意味着这个样本的数据质量有问题，或者它收集时混入了一些非目标组织类型的细胞。

操作方法很直接，对样本做层级聚类，画个树状图看看：

# 转置矩阵，让行变成样本，列变成基因 sampleTree <- hclust(dist(expr_filtered), method = "average") pdf("sample_clustering.pdf", width = 12, height = 6) plot(sampleTree, main = "Sample clustering to detect outliers", sub = "", xlab = "", cex.lab = 1.5, cex.axis = 1.5) abline(h = 阈值, col = "red") dev.off()

这里dist()默认计算欧氏距离，把每个样本的基因表达向量当成一个高维空间里的点，计算点与点之间的距离。两个样本距离越近，表示它们的基因表达模式越相似。method = "average"表示聚类时类间距离取平均距离，这种连接方式相对稳定，是WGCNA官方教程里推荐的选择。

聚类图里如果出现一个样本单独挂在一根很长的分支下面，距离其他样本都很远，就需要考虑是否剔除。判断标准一般看距离阈值，比如把聚类高度大于某个值的样本视为离群。我在实际项目中遇到过不止一次因为某个样本的RNA质量差导致全基因组表达量异常偏低的情况，这种样本不做剔除，后续模块检测结果会被明显带偏。判断的时候要结合样本的QC指标综合判断，不要仅凭聚类图就做决定。如果剔除样本后分组样本量不够，可以考虑补充样本或者谨慎保留，但一定要在文章里如实说明样本筛选过程。

剔除离群样本后的矩阵，建议重新保存一份用于后续分析：

# 假设sampleTree识别出的离群样本是 "SAMPLE_017" keep_samples <- !colnames(expr_filtered) %in% "SAMPLE_017" expr_clean <- expr_filtered[, keep_samples]

3. 软阈值选择：网络构建的关键一步

3.1 无标度网络与软阈值的来龙去脉

进入网络构建之前，先解决一个绕不开的概念问题：为什么WGCNA要用软阈值。

WGCNA假设基因调控网络具有无标度网络的特征。无标度网络的特点是少数节点拥有大量连接（hub），大多数节点只有少数连接，生物学网络比如蛋白质互作网络、代谢网络都近似符合这种特性。如果用硬阈值来构建网络，相关性刚好超过阈值的基因对全部保留、刚好低于阈值的全部丢弃，网络结构容易变得碎片化，不稳定。更麻烦的是硬阈值选择的“拍脑袋”属性太强，阈值从0.8改成0.85，网络结构可能面目全非。

软阈值的作用是给每对基因的连接强度做一个幂函数变换，在保留网络连续性的同时，让网络结构尽量逼近无标度特性。你的任务是找一个合适的(\beta)值，使得加权后的网络连接度分布尽可能拟合无标度分布。WGCNA包里的pickSoftThreshold函数专门做这件事。它会计算一系列候选(\beta)值对应的拟合指数（scale-free fit index，就是R方），同时统计网络的平均连接度，然后帮你画一张图，你根据这张图来选定最终的(\beta)。

通常标准是选择第一个R方达到0.85以上的(\beta)值。但注意，这个0.85不是绝对的。有些数据集比如来自复杂组织的转录组数据，天然就很难达到高R方，这时候综合考虑平均连接度的下降趋势，选一个折中的(\beta)值就能继续分析。另外(\beta)也不需要设得过高，(\beta)太高意味着对强相关的依赖度更大，网络会变得稀疏且过于集中，下游模块检测的稳定性反而受影响。

3.2 用pickSoftThreshold挑选参数并解读输出

下面这段代码是软阈值选择的模板：

# 允许多线程加速 allowWGCNAThreads() # 候选软阈值向量 powers <- c(1:10, seq(from = 12, to = 30, by = 2)) # 计算软阈值 sft <- pickSoftThreshold( datExpr = t(expr_clean), # WGCNA要求样本为行、基因为列 powerVector = powers, verbose = 5 ) # 查看软阈值评估结果 print(sft$fitIndices)

这里有一个非常容易踩坑的地方：datExpr参数需要的是行是样本、列是基因的矩阵，和正常的表达矩阵（行是基因、列是样本）正好相反。所以传参时一定要转置。很多人第一步就在这里报错，要么提示“Error: datExpr must be a data frame or matrix”，要么后续结果乱套，原因就是忘了转置。

pickSoftThreshold的输出里有两列最关键：SFT.R.sq和mean.k.。前者是候选阈值下网络拟合无标度分布的R方，后者是平均连接度。画图可以用WGCNA自带的绘图函数：

# 设置画布分两栏 pdf("soft_threshold.pdf", width = 12, height = 6) par(mfrow = c(1, 2)) # 左图：无标度拟合指数 plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit (Signed R^2)", type = "n", main = "Scale independence") text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels = powers, col = "red") abline(h = 0.85, col = "red") # 右图：平均连接度 plot(sft$fitIndices[, 1], sft$fitIndices[, 5], xlab = "Soft Threshold (power)", ylab = "Mean Connectivity", type = "n", main = "Mean connectivity") text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, col = "red") dev.off()

这段绘图代码里有个细节，-sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2]是在处理无标度拟合的判定系数。WGCNA包输出里第三列是斜率，斜率为负代表网络满足无标度特性（即高连接度节点数量少），所以取负号把R方修正为正值。如果某一行的斜率是正的，说明网络结构发散，不满足无标度假设，这个候选阈值就不适合。

实际判断的时候，我会先看左图，找R方首次超过0.85的(\beta)，然后看右图确认平均连接度没有跌到特别低。比如之前分析一个植物转录组数据，R方在(\beta=6)时达到了0.91，但(\beta=6)时平均连接度还有70多，网络不会太稀疏，我就选了6。另一个动物组织数据，R方在(\beta=12)时是0.87，但平均连接度只有6，模块后续非常碎，这时候我会下调到(\beta=10)来平衡，保证模块的规模能用生物学语言解释。

4. 网络构建与模块识别实操

4.1 一步法构建网络与参数含义

选定软阈值之后，就可以构建网络并识别模块了。有两种路径：一种是完全手动分布执行，先算邻接矩阵、再算TOM矩阵、再做层级聚类、最后动态剪枝；另一种是直接用blockwiseModules函数一步完成。日常分析推荐直接使用后者，它把全套流程封装好了，参数也足够灵活。

# 网络构建与模块识别 net <- blockwiseModules( datExpr = t(expr_clean), power = 6, # 上一步选定的软阈值 TOMType = "unsigned", # TOM的计算类型 minModuleSize = 30, # 最小模块基因数 reassignThreshold = 0, # 重分配阈值 mergeCutHeight = 0.25, # 模块合并的聚类高度阈值 numericLabels = TRUE, # 模块用数字标识 pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "TOM_data", # TOM矩阵保存的文件名 verbose = 3 )

这段代码里最容易让人犹豫的是TOMType参数。它有三个取值：unsigned、signed和signed hybrid。unsigned只考虑相关性的绝对值，基因正相关和负相关都会算作连接；signed只把正相关算作连接，负相关不计入。转录组分析中，共调控的基因通常是正相关关系，signed更符合生物学预期，但负相关的基因也有可能属于同一个调控通路上的抑制关系。我个人做转录组默认用signed，如果做的是芯片数据且表达谱分布比较宽，可以考虑unsigned。关键是前后一致，不要一会儿signed一会儿unsigned导致结果对不上。

mergeCutHeight是模块合并的阈值。因为直接聚类出来的模块可能非常碎，把聚类树上距离接近的模块合并掉，可以减少下游分析中模块太多、过度碎片化的问题。值越小合并越少，常见范围0.15到0.3。我一般先用默认0.25跑一遍，看模块数量和合并结果，如果模块数量太多、每个模块基因太少，就调大一点；如果模块数太少导致很多表达模式差异明显的模块被合并吞掉，就调小一点。这一点需要反复尝试，结合你后续模块与性状关联的结果来最终确定。

4.2 识别模块结果与模块特征基因

运行结束后，net对象里包含多个结果组件。最常用的是net$colors，它给每个基因打了一个模块编号的标签。数字0表示没有被分配到任何模块的基因，这些基因会被后面的分析忽略。可以用代码快速看下模块的基因数分布：

# 查看模块颜色与基因数 moduleColors <- labels2colors(net$colors) table(moduleColors)

labels2colors函数会把数字标签转成颜色名，方便画图时使用。模块基因数如果出现很多个低于30的小模块，说明minModuleSize设小了或者mergeCutHeight设小了。对于下游分析，模块大小在50到2000之间比较理想。太小的模块生物学重复性差，太大的模块内部异质性高，难以解读。

模块识别完之后，WGCNA会为每个模块计算模块特征基因（module eigengene），用来代表这个模块的整体表达模式：

# 计算模块特征基因 MEs <- moduleEigengenes(t(expr_clean), moduleColors)$eigengenes # 查看模块特征基因的样本×模块矩阵 head(MEs)

特征基因本质上是对模块内所有基因的表达矩阵做PCA，取第一主成分。它用一组数值概括了每个样本在该模块上的整体表达水平。有了这个值，后面跟表型做相关性分析就变得非常直接。一个模块的特征基因如果与某个细胞类型、临床指标高度相关，那这个模块很可能反映了该表型相关的生物学过程。

查看模块特征基因之间的相关性，以及样本间的模块表达模式，通常用moduleEigengenes之后接着做heatmap或者plotEigengeneNetworks。这一步能帮你快速发现哪些模块之间表达模式相似，可能共同参与某个通路。我在实践中经常发现两个模块的特征基因相关系数超过0.8，这就是你后续做模块合并或者解读功能时要特别注意的信号。

5. 常见报错与排查技巧：实操中踩过的坑

5.1 数据格式与数据处理类报错

从一开始就把格式处理好，能避免很多不必要的折腾。但我还是把实际中遇到过的典型报错整理出来，方便你排查对照。

报错一：Error: 'datExpr' must be a data frame or matrix
几乎都是因为没转置矩阵。WGCNA整个流程里，datExpr要求的是行是样本、列是基因。如果你读入的表格是行基因列表，记得在传给pickSoftThreshold或blockwiseModules前运行t()转置。还有个容易忽略的点，t()之后原本的data.frame会变成matrix，有些后续操作需要再转回data.frame，否则某些函数会报类型错误。

报错二：Error: 'y' must be numeric或聚类时报相似性矩阵找不到
这通常是因为表达矩阵里有非数值列，比如基因注释列没有去掉。dist()函数只能对数值矩阵计算距离，如果你的矩阵里混入了字符型数据，就会报错。解决方法是在读取数据后，把基因注释列拆出来单独存，表达量部分的矩阵确保每个列都是numeric类型。

报错三：运行blockwiseModules时内存不足或卡死
当基因数量超过2万、样本数量超过100时，TOM矩阵的计算量迅速增长，内存消耗会非常惊人。一种方案是设置maxBlockSize参数，让WGCNA自动分块计算。另一种是在服务器上跑，增加内存。如果只能在本机跑，我建议把输入基因先做个过滤，比如只保留MAD前10000的基因，能大幅降低计算量。实际分析中影响不大，但内存问题能很快缓解。

5.2 软阈值选择与模块识别的思路调整

软阈值选定后，模块识别结果不理想的情况也经常出现。比较常见的情形有两种：一是没有达到0.85的R方，二是模块太多太碎、难以解读。这两种情况处理思路不太一样。

先看第一种，没有达到0.85。这时需要回到软阈值图，找一个拐点，也就是R方上升变缓的位置。比如候选阈值从6到8，R方从0.82涨到0.86，但平均连接度从80跌到30，我一般会选6而不是8。因为更高的(\beta)让网络过于稀疏，模块会变得碎片化，后续结果稳定性差。0.85是参考值而不是绝对标准，判断时看到的是网络拓扑特性整体是否合理。

再看第二种，模块太多太碎。最常见的原因是mergeCutHeight太小，模块合并不够。我的习惯是先输出不同mergeCutHeight值下的模块数量，画个对照表，再从中挑一个模块数量适中、模块大小分布合理的值。另一个原因是minModuleSize设太小，如果设置15，很容易产生大量几十个基因的小模块。建议至少设到30，如果模块数量还是多，可以试试50。当然，如果你的研究目标就是找特定的小模块，那可以适当调低阈值，但要做好结果稳定性检查。

操作中还有一个小细节值得注意：WGCNA的运行结果和随机种子有关。blockwiseModules内部用到的一些随机过程不是完全确定的，如果你的分析文章需要可重复性，建议在代码开头设置随机种子：

set.seed(2024)

这样别人按同一份代码跑，结果能完全复现。不同版本R或不同平台下分块方式不同会导致结果细微差异，但设置种子能保证绝大多数情况下的一致性，写文章和做报告时这个细节会显得你很专业。

实操总结与两点补充心得

到这里，WGCNA分析流程的第一阶段就跑通了。梳理一下第一篇的内容：数据清洗与低表达基因过滤、样本聚类与离群样本剔除、软阈值选择、一步法网络构建与模块识别、模块特征基因计算。跑完这套流程，你已经能看到基因模块的概貌了。第二篇会继续往前走，做模块与性状的关联分析、模块内基因的功能富集解读、hub基因筛选，以及最后那些漂亮的网络可视化图是怎么画出来的。

最后分享两个实际操作中积累的小经验。第一个是关于数据分析记录的，WGCNA分析涉及大量参数选择和反复尝试，建议每次调参都记录对应的结果，比如模块数量、模块大小分布等，方便事后追溯和复现，也方便写文章时说明参数选择的依据。第二个是关于结果的解释，WGCNA跑出来只是第一步，模块里富集到哪些通路、与哪些性状显著相关，才决定模块能否转化为一个有故事可讲的生物学发现。这两点在第二篇里都会有更详细的实际操作演示，到时候结合代码一起看会更有感觉。