WGCNA 加权基因共表达网络分析
WGCNA（Weighted Gene Co-expression Network Analysis）根据基因在样本间的表达相关性构建加权网络，把共同升降的基因聚成模块，用模块特征基因（模块表达的第一主成分）与表型做关联，找出和性状最相关的模块及其中的 hub 基因。它关注的是「哪一群基因一起变化、整体和什么性状有关」，适合样本较多、表型连续或多分组的数据。典型输入是 bulk RNA-seq 或芯片的标准化表达矩阵加样本性状表，标准实现是 R 包 WGCNA。
什么时候用
- 样本数较多的转录组，希望找到与表型变量相关的基因群
- 表型是连续变量，如体重、某项生化指标、药物浓度或时间点
- 差异基因很多，需要压缩成少数几个功能单元来解释
- 从与性状最相关的模块中挑选 hub 基因作为实验候选
- 比较两种条件或两个数据集之间模块是否保守
不适合的情况
- 样本少于 15 个，官方 FAQ 不建议使用，相关系数噪声过大
- 只用差异基因建网，会破坏网络结构并夸大模块与分组的相关
- 直接用单细胞稀疏计数矩阵，零值主导相关，应改用基于 metacell 的方法
需要准备的数据
- 表达矩阵
- 样本 × 基因，已标准化的 log 尺度数据，如 DESeq2 的 vst 结果或芯片 RMA 值，不用原始 count
- 性状表
- 行名与表达矩阵样本一致；分类变量拆成 0/1 哑变量，连续变量保留原值
- 基因筛选规则
- 常按方差或中位绝对偏差保留前 5000–10000 个基因，并去除低表达基因
- 批次信息
- 存在批次时需在建网前校正，否则模块会按批次划分
分析步骤
- 1
过滤与样本质控
去除低表达基因，按变异度保留数千个基因；用 goodSamplesGenes() 检查缺失值和零方差基因。对样本做层次聚类，剔除明显离群的样本，并在方法中写明剔除依据。
- 2
选择软阈值
用 pickSoftThreshold() 计算不同 power 下的无标度拟合 R² 和平均连接度，选 R² 首次达到约 0.8–0.85 的最小 power。若始终达不到，先排查批次效应或离群样本。
- 3
构建网络识别模块
常用 blockwiseModules() 一步完成邻接矩阵、TOM、动态剪切树和模块合并。检查树状图下的颜色条是否与分支对应；灰色是未分配基因，占比过高说明参数或数据有问题。
- 4
模块-性状关联
计算每个模块特征基因与各性状的相关系数及 P 值并画热图。挑选相关绝对值高且 P 值小的模块，同时排除那些与批次、RNA 质量等技术变量强相关的模块。
- 5
计算 GS 与 MM
基因显著性（GS）是基因与性状的相关，模块成员度（MM，即 kME）是基因与模块特征基因的相关。在目标模块画 MM-GS 散点图，二者正相关说明模块核心基因也与性状相关。
- 6
筛选 hub 基因
常用经验阈值是 |MM|>0.8 且 |GS|>0.2，也可按模块内连接度取前若干个。hub 列表再与差异分析、PPI 网络或外部数据取交集，缩小到少数优先验证的候选。
- 7
注释与导出网络
对目标模块基因做 GO/KEGG 富集，确认模块有一致的功能主题；用 exportNetworkToCytoscape() 按 TOM 阈值导出边表和节点表，用于网络可视化。
关键参数与决策
软阈值 power
R² 达不到时，WGCNA FAQ 给出按样本数取值的经验表，如无符号网络在 40 个样本以上约取 6，样本越少取值越高，signed 网络约为其两倍。
networkType
unsigned 把正负相关都视为连接，signed 只把正相关视为强连接。若希望模块内基因同向变化、便于解释，常选 signed 或 signed hybrid，并相应提高 power。
minModuleSize
最小模块基因数，常设 20–30。设得过小会产生大量碎片化小模块，设得过大则会把有意义的小模块并入灰色。
mergeCutHeight
特征基因相异度低于该值的模块会被合并，blockwiseModules 默认 0.15，即相关约 0.85 以上合并。模块过多且彼此相似时可调到 0.25 左右。
maxBlockSize
分块计算的基因上限，默认 5000。内存允许时设为大于基因总数，避免分块造成同一共表达群被拆进不同块。
deepSplit
动态剪切树的灵敏度，取 0–4，默认 2。值越大模块越多越小，一般在 2 附近调整，并观察模块数量和灰色比例的变化。
常见错误
- 用未经 log 转换的 count 或 TPM 建网：高表达基因主导相关性，模块反映的是表达量高低，而不是共调控关系。
- 先筛差异基因再做 WGCNA：网络被分组信号主导，几乎必然得到与分组强相关的模块，属于循环论证，官方 FAQ 明确不建议。
- cor 函数被其他包覆盖：加载多个包后 cor 可能冲突导致报错，常见做法是在脚本开头写 cor <- WGCNA::cor，结束后恢复。
- 分类性状用 1、2、3 编码直接做相关：数值大小被当成顺序，结果没有意义，应拆成多个 0/1 变量分别计算。
- 只看模块-性状热图上 P 值最小的一格：样本少时显著性很脆弱，应同时看相关系数大小以及 MM-GS 关系是否一致。
- 把 hub 基因直接写成关键调控因子：hub 只是共表达网络中连接度高，不代表因果或调控关系，需要实验证据支持。
结果怎么看
软阈值选择图
左图是不同 power 下的无标度拟合 R²，右图是平均连接度。选 R² 达标的最小 power；若 R² 一直很低，或高 power 下平均连接度仍极高，提示存在离群样本或强批次效应。
基因聚类树与模块颜色
树状图下方的颜色条代表模块。理想情况下同一分支对应同一颜色；灰色占比过大，如超过一半，说明多数基因没有形成稳定模块，需要检查过滤规则和 power。
模块-性状相关热图
每格给出模块特征基因与性状的相关系数和 P 值。重点关注相关绝对值高且 P 值小的模块，同时看它与技术变量的关联；与技术变量也强相关的模块，结论要保守。
MM-GS 散点图
横轴为模块成员度，纵轴为基因显著性。点云明显正相关，说明越核心的基因与性状越相关，模块-性状关联更可信；若无相关，热图上的显著可能只是偶然。
模块特征基因关系图
用 plotEigengeneNetworks() 画模块特征基因的聚类树和相关热图，可把性状一起放入。彼此高度相关的模块往往代表同一过程，可考虑合并；与性状聚在同一分支的模块是重点。
和相近方法怎么选
WGCNA 与差异表达分析
差异分析逐个基因比较两组，得到单基因的变化；WGCNA 利用全部样本的相关结构找协同变化的基因群，适合连续性状和多分组。两者常联用：模块基因与差异基因取交集，缩小候选范围。
WGCNA 与 PPI 网络
WGCNA 的边来自本数据的表达相关，反映共表达；PPI 的边来自数据库收录的蛋白互作，与本数据无关。前者适合发现数据驱动的模块，后者适合为候选基因补充已知互作证据。
WGCNA 与 hdWGCNA
单细胞数据过于稀疏，直接计算相关会被零值主导。hdWGCNA 先把相似细胞聚合为 metacell，再沿用 WGCNA 框架建网，适合在特定细胞类型内寻找共表达模块。
常见问题
WGCNA 需要多少样本？
WGCNA 官方 FAQ 认为少于 15 个样本时不建议使用，20 个以上结果更稳定。样本越少，相关系数越不可靠，模块越容易被噪声或个别离群样本驱动。
软阈值 R² 达不到 0.8 怎么办？
先检查数据：是否有离群样本、批次效应，或强烈的分组效应让大量基因同时升降。排除这些原因后仍达不到，可参照官方 FAQ 按样本数给出的经验 power，并在方法中说明。
WGCNA 的 hub 基因怎么筛选？
常见做法是在目标模块中按 |MM|>0.8 且 |GS|>0.2 筛选，或按模块内连接度取前 10 个。这些是经验值，可按模块大小调整；筛出后最好与 PPI、差异分析结果取交集。
灰色模块是什么意思？
灰色模块收纳的是没有被分配到任何模块的基因，并不是真正的共表达模块，不应拿去做性状关联或富集分析。灰色基因比例很高时，需要回头调整过滤规则和网络参数。
做 WGCNA 之前需要去批次吗？
需要。批次带来的系统差异会让大量基因同步升降，被识别成一个或多个模块，并与批次分布相同的性状产生虚假相关。可在建网前用 ComBat 或 limma 的 removeBatchEffect 处理。
在 Codar 里怎么做
在 Codar 中，上传标准化表达矩阵和样本性状表，Codar 会先确认基因过滤规则、网络类型和要关联的性状，再完成离群样本检查、软阈值选择、模块识别、模块-性状关联、GS 与 MM 计算和 hub 基因筛选。交付软阈值图、模块树、模块-性状热图、MM-GS 散点图、hub 基因表、报告和脚本。
我有 48 个小鼠肝脏样本的 vst 表达矩阵，以及体重、空腹血糖、肝脏甘油三酯等表型，想用 WGCNA 找与肝脂最相关的模块和 hub 基因。