WGCNA+机器学习+孟德尔随机化：从共表达网络到因果验证的生信分析组合

如果你做的生信项目还停留在“跑个差异基因、画个火山图”的阶段那这套组合拳——整合WGCNA、机器学习和孟德尔随机化——我强烈建议你花一个周末学一下。这不是普通的发文章套路而是把“数据挖掘”和“因果推断”接成一条完整的证据链先用WGCNA从全局网络里抓出和疾病相关的共表达模块再用机器学习从模块里筛出真正有区分度的核心基因最后用孟德尔随机化给这些基因补上因果层面的验证。适合做肿瘤、代谢病、免疫性疾病等方向的研究生和科研人员尤其适合手上有公共数据、缺一个新分析亮点的朋友。1. 为什么把WGCNA、机器学习和孟德尔随机化放在一起用1.1 三种工具各解决什么问题先说WGCNA。它的全称是Weighted Gene Co-expression Network Analysis中文叫加权基因共表达网络分析。它做的是无监督的“分群”把表达模式相似的基因聚成一堆每个堆叫一个模块module。为什么要这么做因为基因不是孤立工作的功能相关的基因往往共表达。你拿到一个转录组矩阵动辄上万个基因直接做差异分析只能看到单基因的上下调看不到“整组基因协同变化”的规律。WGCNA把上万个基因压缩成几十个模块再拿模块的特征基因module eigengene去跟临床性状做关联一下子就能锁定跟疾病最相关的“基因集团”。机器学习在这里的作用完全不一样。它是监督学习是在你有明确分组标签比如病例/对照的情况下从一堆候选基因里找到那组组合起来预测效果最好的特征。为什么要放在WGCNA后面因为直接拿全部基因扔给LASSO或随机森林噪声太大容易出现“用一堆不相关基因碰巧拟合”的过拟合情况。WGCNA先帮你把范围缩到那几个跟性状相关的模块里机器学习的任务就从“大海捞针”变成了“从几百个备选里挑精兵”。孟德尔随机化Mendelian Randomization, MR是压轴的因果验证。前面两种方法说到底都是相关性分析——基因表达和疾病相关但不能证明因果关系。MR利用遗传变异通常用eQTL作为工具变量模拟随机分组的效果。因为等位基因在受精时就随机分配不受环境、生活方式等混杂因素影响所以如果某个基因的表达水平受某个SNP的调控而这个SNP又和疾病显著关联那就等于给“基因表达→疾病”这条关系线上了一层因果可信度。1.2 整合逻辑从相关性到因果性的证据链这套流程的本质是在一条链上不断过滤和升级证据。第一步表达数据筛选。你用GEO、TCGA这些数据跑WGCNA找到与疾病相关的模块得到一批候选基因。这批基因是有生物学背景的因为它们处在同一个共表达网络里可能是同一通路或同一细胞类型的标记。第二步机器学习定位核心基因。同一个模块里可能有几百个基因不是每个都值得进入后续验证。机器学习做的事是看哪些基因的组合能最准确地把病例和对照分开。这里我推荐至少用两种算法交叉验证比如LASSO和随机森林取交集或者SVM-RFE和逻辑回归取交集。这样做的好处是避免单一算法对数据分布的偏好。第三步MR做因果兜底。机器学习筛出来的基因放到MR框架里。你去找这些基因在组织中最常用血液或肝脏的cis-eQTL数据然后拿这些SNP作为工具变量从GWAS summary数据里提取对应信息看它们对目标疾病的因果效应。为什么这个顺序不能乱因为如果你先跑MR候选基因通常来自差异分析或文献数量大且缺少表达网络背景如果你先跑机器学习再倒回去做WGCNA也说得通但逻辑上不如先利用无监督网络获取全局结构、再利用监督学习精挑、最后用遗传学证据兜底来得顺畅。审稿人很喜欢这种“层层递进、逻辑自洽”的分析链。2. WGCNA实操细节拿到表达矩阵后先做这些2.1 数据清洗与样本过滤WGCNA对数据质量要求很高这一步偷懒后面全白做。我习惯的步骤是去低表达基因。在大多数表达芯片里我会设定一个标准至少有80%的样本中表达量大于0或log2CPM大于1否则剔除。RNA-seq数据建议用CPM或TPM转换一般会过滤掉那些在所有样本里都接近0的基因不然它们在后续相关性计算里只会制造噪声。样本聚类看离群值。这一步用hclust对所有样本做层次聚类。如果一个样本离群太远往往是芯片质量差、批次效应严重或样本标注错误。我一般会用1.5倍四分位距IQR作为初筛标准见到明显离群样本就剔除。如果你连自己的样本离群都没发现后续模块分析可能直接被一个坏样本带偏。批次效应校正。如果样本来自不同批次或不同平台必须用limma包的removeBatchEffect或sva包的ComBat函数做校正。注意校正时不要混入结局变量信息否则会引入数据泄露。这里很多人容易犯错拿着包括结局在内的完整矩阵去调批次结果下游差异分析和机器学习全部虚高。标准化。使用z-score或标准化到均值为0、标准差为1。WGCNA本身在算相关系数时相当于自动标准化了但后续模块特征基因的数值尺度统一对模块-性状关联会友好一些。2.2 软阈值选择与模块识别WGCNA的核心思想是构建“软阈值”网络——不简单地把基因间相关系数是否大于某个阈值变成0/1连接而是对相关系数做幂指数变换让网络更接近无标度拓扑scale-free topology。这个无标度性质是什么意思简单说真实生物网络里少数节点连接很多基因多数节点连接很少基因度分布近似幂律。我们要选一个power值让网络尽可能接近这个状态。实操上我一般用pickSoftThreshold函数计算一系列power通常1到20下的无标度拟合指数R²和平均连接度。我选power的原则R²尽量大于0.8越高越好同时平均连接度不要太高否则网络太密模块分不开如果多个power都满足R²0.8我倾向于选最小的那个因为幂越小网络保留的信息越多。选完power后运行blockwiseModules或更手工的步骤。blockwiseModules一步到位很省事但我更推荐拆开做方便调试。模块识别用动态剪枝法dynamic tree cut。关键参数是deepSplit控制分支切割的敏感度。范围从0到4数值越大模块越多越细。我通常从deepSplit2开始如果模块数量太多比如多于50个就调小一点如果模块太粗糙就调大到3试。还有一个minModuleSize最小模块基因数默认30。小样本时建议设成40或50避免出现太多小碎块。模块识别完之后立刻算每个模块的特征基因ME。特征基因相当于这个模块里所有基因表达模式的“主成分第一轴”。然后计算每个模块ME和临床性状比如患病与否、肿瘤分期的相关系数和P值。重点关注显著模块的基因集。2.3 我踩过的坑第一个坑是样本量不足。如果样本少于15个WGCNA的模块稳定性很差你今天跑出来的模块换一批样本可能完全不一样。一般建议至少20个样本最好30个以上。如果你手里只有十几例能不做WGCNA就别做强行做出来的结果自己都不信。第二个坑是软阈值选不出来。有时候扫完1到20最大R²也到不了0.8。这时候先检查是不是有离群样本没干掉再看数据是否过于同质。如果还不行考虑用signed network代替unsigned。signed network只保留正相关基因的连接更符合“共表达”的生物学含义通常更容易达到无标度拟合。第三个坑是模块过多或过少。模块太多后面机器学习筛选时容易过拟合模块太少又会丢掉很多候选基因。我每次拿到blockwiseModules结果都会先看模块数量和每个模块的基因数分布。如果看到一堆只有四五十个基因的碎模块我会调大deepSplit和minModuleSize再做一次聚类合并。WGCNA自带mergeCloseModules按模块特征基因相关性阈值合并模块我用0.25作为默认阈值调低到0.2能并得更少。3. 机器学习二次筛选怎么从模块里选出真正有区分度的基因3.1 算法选型LASSO、随机森林、SVM-RFEWGCNA给你的是一个候选基因集合——通常是某个或某几个疾病相关模块里的全部基因数量可能还有几百个。这时候直接拿去做湿实验验证成本和风险都很高。机器学习在这里的价值是压缩候选基因数量同时保证分类能力。我经常用的三种算法各有特点LASSOL1正则化逻辑回归用L1惩罚让不重要的特征系数变成0适合在高维低样本情境下筛选特征。它的输出是一组非零系数的基因可解释性强系数大小也代表贡献程度。随机森林Random Forest基于决策树的集成学习能捕捉非线性关系自带特征重要性排序。它对多重共线性的耐受度比较好而WGCNA模块里基因之间往往高度共表达这恰恰是随机森林的优点。SVM-RFE支持向量机递归特征消除用SVM反复训练模型每轮去掉最不重要的特征。效果通常不错但计算量大在几百个基因时还能跑上千个就有点慢。我的经验是不要只用一种方法。选择两种或三种取交集或排名靠前的基因。一般来说LASSO和随机森林的结果重合度不会太高因为二者逻辑不同——一个偏线性稀疏一个偏非线性交互。取交集能筛出那些“无论从哪个角度看都重要”的基因可靠性更高。3.2 实操流程训练集/测试集切分机器学习的部分我会用R语言的caret、glmnet和randomForest包也可以用Python的scikit-learn原理都一样。关键是流程要规范否则结果虚高。标准化特征。先把候选基因的表达量做z-score标准化。这一步不能漏因为LASSO和SVM对特征尺度敏感。划分训练集和测试集。一般70%/30%。注意分层抽样让病例和对照的比例在训练集测试集中保持一致。在训练集上用交叉验证调参。LASSO我直接用cv.glmnet默认是10折交叉验证lambda选择标准用lambda.min或lambda.1se。我偏好lambda.1se它给出的模型更稀疏可重复性更高。随机森林用ranger包fast或randomForestntree设1000mtry调优交叉验证可以在caret里做。性能评估。用测试集计算AUC、准确率、敏感度和特异度。AUC低于0.7的模型说实话在生信文章里说服力有限。我之前遇到过一个数据集测试集AUC有0.95看起来很漂亮但后来发现是训练集测试集划分时没有按样本来源分层导致相同病人的不同样本同时跑到了两边。这个教训很痛——一定要保证样本独立性。设定随机种子。所有随机过程都要seed不然结果无法复现。我通常用set.seed(2024)这样的统一种子跑完记录一下。LASSO代码大致是library(glmnet) # expr.matrix: 候选基因表达矩阵已标准化 # group: 0/1分组向量 set.seed(2024) cv.fit - cv.glmnet(as.matrix(expr.matrix), group, family binomial, alpha 1, type.measure auc) plot(cv.fit) coef.min - coef(cv.fit, s lambda.1se) # 提取非零系数基因 selected_genes - rownames(coef.min)[which(coef.min[,1] ! 0)][-1]随机森林代码大致是library(randomForest) set.seed(2024) rf.fit - randomForest(as.factor(group) ~ ., data data.frame(expr.matrix), ntree 1000, importance TRUE) importance(rf.fit) # 按MeanDecreaseGini排序后取Top503.3 小心过拟合和数据泄露这是机器学习应用里最致命的问题。WGCNA模块基因之间高度相关LASSO虽然能解决部分共线性但随机森林和SVM-RFE在高度相关的特征下有可能把“一组基因的重要性”分配到其中某一个导致你选出一种基因完全取决于随机种子。所以我的习惯是至少做三次不同种子的重复筛选只保留那些在多次筛选里稳定出现的基因。还有一个常见坑是特征选择的嵌套问题。有人直接在全体数据上做特征选择再做交叉验证评估AUC这实际上已经泄露了测试集信息。正确做法是把特征选择放进每一折交叉验证里做或者把特征选择和模型调参都限制在训练集内。如果你用的是caret的rfe函数它会自动处理嵌套但如果你手工选基因再扔进模型注意别用全样本信息做选择。4. 孟德尔随机化验证给候选基因补上因果一环4.1 MR基础概念与假设孟德尔随机化本质上是用遗传变异作为工具变量去检验“暴露因素这里是基因表达量”和“结局疾病”之间的因果关系。核心基于三大假设关联性假设工具变量遗传变异与暴露因素强相关。实际操作中我们使用基因的cis-eQTL——即距离基因座一定范围内与基因表达水平显著相关的SNP。F统计量一般要求大于10说明工具变量不是弱工具。独立性假设工具变量与混杂因素独立。遗传变异是出生时随机分配的所以理论上不受后天环境、生活习惯等常见混杂因素影响。这也是MR相比传统观察性研究的最大优势。排他性假设工具变量只通过暴露因素影响结局没有其他路径。这是最难验证的假设。如果某个SNP不仅影响基因表达还通过其他通路直接致病就违反了排他性假设。这也是为什么我们要做敏感性分析。4.2 基于GTEx eQTL和GWAS数据的两样本MR流程当你通过机器学习筛出若干个候选基因后接下来面对的任务是为每个基因找到合适的工具变量并用它来做MR分析。具体操作我会用TwoSampleMR这个R包配合openGWAS数据库或者本地的eQTL和GWAS目录。获取基因的cis-eQTL。我通常用GTEx V8数据选择与你的研究组织最接近的组织类型。比如研究脂质代谢疾病用肝脏的eQTL研究神经系统疾病用大脑的eQTL。表达量化用“normalized expression”提取基因附近的显著eQTL一般以转录起始位点TSS为中心上下各500kb或1Mb并保留P值小于5e-8的SNP。对eQTL进行LD clumping。同一区域内的SNP之间往往存在连锁不平衡不能都作为独立工具变量。使用PLINK的clump方法r²阈值设为0.1窗口设为500kb保留P值最显著的SNP。这一步的目的是尽可能去除工具变量之间的相关性避免违反工具变量独立性假设。提取GWAS数据。从GWAS Catalog或者专门的疾病GWAS summary数据中提取与你关注的疾病比如冠心病、2型糖尿病相关的对应SNP的效应值、效应等位基因等。注意数据格式对齐效应等位基因effect allele和参考等位基因other allele必须和eQTL数据中的方向一致否则要回补。运行MR分析。TwoSampleMR中先把工具变量数据合并好再调用mr()函数使用IVW逆方差加权作为主分析。同时我还会跑MR-Egger、加权中位数法、简单众数法、加权众数法作为敏感性分析。IVW是假设所有工具变量都有效时使用的方法相当于每个工具的Wald比率做倒方差加权平均。如果结果在IVW中显著但MR-Egger方向不一致或者MR-Egger截距显著说明可能存在多效性。4.3 怎么解读MR结果当拿到MR结果时首先看IVW的P值。如果P值小于0.05说明基因表达水平与疾病风险存在因果效应。接着看OR值或beta值判断方向——OR大于1表示基因表达升高会增加疾病风险小于1表示有保护作用。再看敏感性分析是否一致。如果只有IVW显著其他方法不显著结果的可信度就一般。如果IVW、加权中位数、MR-Egger方向一致且至少有两种方法P值显著这个因果关系的证据就相对扎实。此外要做一个MR-PRESSO检验检测和修正多效性异常值。MR-PRESSO可以识别出那些对结果影响很大的SNP多效性SNP剔除后重新分析看结果是否稳健。我会把这些结果整理成一张表格方法OR95%CIP值IVW1.121.03-1.220.008MR-Egger1.180.99-1.400.067加权中位数1.101.00-1.210.045加权众数1.090.96-1.240.187如果IVW和加权中位数都显著MR-Egger方向一致但没到显著也可以认为因果证据相对支持。但如果MR-Egger方向和IVW完全相反那就要小心了——多效性可能使IVW结果失真。5. 三个环节无缝衔接一个完整的示例流程5.1 示例数据背景为了让你更直观地理解整个流程我以“某种常见慢性炎症性疾病”为例设计一个标准分析路径。假设你有一组疾病和对照的血液转录组表达数据60个样本疾病30对照30公开的血液eQTL数据比如GTEx全血该疾病的公开GWAS summary数据。目标是通过WGCNA机器学习锁定疾病相关的核心基因再用MR验证这些基因与该疾病的因果关联。5.2 从表达矩阵到MR结果的代码级走读第一步数据导入与清洗。假设你已经获得了表达矩阵expr行是基因列是样本和一个临床信息表meta里面有sample、group1表示疾病0表示对照。用limma包做标准化和去批次然后过滤低表达基因library(limma) library(WGCNA) allowWGCNAThreads() # 过滤低表达基因至少在50%样本中log2CPM1 expr - expr[rowSums(expr 1) 0.5 * ncol(expr), ] # 对样本做聚类看离群 sampleTree - hclust(dist(t(expr)), method average) plot(sampleTree)如果发现离群样本根据聚类结果移除然后重新聚类。第二步软阈值选择与模块识别。datExpr - t(expr) # WGCNA要求样本在行基因在列 powers - c(1:20) sft - pickSoftThreshold(datExpr, powerVector powers, verbose 5) # 看R²是否达到0.8 par(mfrow c(1,2)) plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3]) * sft$fitIndices[,2], xlabSoft Threshold (power), ylabScale Free Topology Model Fit (R²)) abline(h0.8, colred)假设在power为8时R²达到0.85我就接着跑net - blockwiseModules(datExpr, power 8, networkType signed, TOMType signed, minModuleSize 40, deepSplit 2, mergeCutHeight 0.25, numericLabels TRUE, verbose 3)拿到的net$colors就是每个基因的模块编号0代表未分配到任何模块。我再计算模块特征基因与分组关联MEs - moduleEigengenes(datExpr, colors net$colors)$eigengenes modTraitCor - cor(MEs, meta$group, use p) modTraitP - corPvalueStudent(modTraitCor, nrow(datExpr))选出P值小于0.05的模块提取对应基因作为机器学习的输入。第三步机器学习筛选。将提取的模块基因表达矩阵整理成datML配合group标签用LASSO和随机森林分别筛选library(glmnet) library(randomForest) set.seed(2024) # 划分训练集和测试集 idx - createDataPartition(group, p 0.7, list FALSE) trainX - as.matrix(datML[idx, ]) testX - as.matrix(datML[-idx, ]) trainY - group[idx] testY - group[-idx] # LASSO cv.fit - cv.glmnet(trainX, trainY, family binomial, alpha 1, type.measure auc) plot(cv.fit) lasso_genes - rownames(coef(cv.fit, s lambda.1se))[which(coef(cv.fit, s lambda.1se)[,1] ! 0)][-1] # 随机森林 rf.fit - randomForest(as.factor(trainY) ~ ., data data.frame(trainX), ntree 1000, importance TRUE) rf_imp - importance(rf.fit) rf_genes - rownames(rf_imp)[order(rf_imp[, MeanDecreaseGini], decreasing TRUE)[1:30]]最终取交集或前后排名都靠前的基因作为候选基因列表。第四步MR验证。对每个候选基因从GTEx全血eQTL提取cis-eQTLLD clumping后用TwoSampleMR做MR。library(TwoSampleMR) # 读入eQTL和GWAS数据假设格式已经处理 # eQTL_data 包含 SNP, effect_allele, other_allele, beta, se, pval, gene # gwas_data 包含 SNP, effect_allele, other_allele, beta, se, pval, disease # 选某个基因的eQTL执行clump exposure_dat - format_data(eQTL_data, type exposure) # 使用本地或API进行LD clumping可以参考ieugwasr::ld_clump exposure_dat - clump_data(exposure_dat, clump_kb 500, clump_r2 0.1) outcome_dat - format_data(gwas_data, type outcome) dat - harmonise_data(exposure_dat, outcome_dat) res - mr(dat, method_list c(mr_ivw, mr_egger_regression, mr_weighted_median, mr_weighted_mode)) res如果工具变量数量足够再加一个MR-PRESSO检测。5.3 结果整合与报告到了这一步你手里有三大块结果WGCNA找到的疾病相关模块基因机器学习筛选出的核心候选基因可能有10到30个MR验证后具有显著因果效应的基因可能只有2到5个。最终结论里我会以MR显著的基因作为“重点基因”来展示用维恩图展示三者的交集。实际汇报中可以展示一个表格包含每个基因在WGCNA模块中的身份、机器学习排序、MR的OR值、P值和敏感性分析结果。审稿人看这种表会非常舒服因为你把证据链写得清清楚楚。6. 常见问题与排查技巧实录6.1 常见问题速查表问题可能原因解决思路软阈值R²始终低于0.8样本离群、数据噪声大、芯片类型不适合去除离群样本、检查是否要校正批次、改用signed network模块数量太多或太碎deepSplit过大、minModuleSize太小调低deepSplit到1或0调高minModuleSize到50LASSO筛选结果不稳定特征高度共线性、样本量小设多个种子重复跑只保留重复出现的基因改用弹性网络alpha0.5随机森林训练集AUC很高测试集不高过拟合缩小特征数量、加正则化、换更简单的模型、检查数据泄露MR结果IVW显著但MR-Egger方向相反存在多效性用MR-PRESSO剔除异常SNP报告多效性检验结果解释要谨慎eQTL工具变量F值不够高选择的SNP距离基因太远或效应弱使用更严格的P值阈值、扩大cis区域到1Mb、检查组织特异性多个基因MR结果都不显著候选基因范围太窄或表型异质性可以尝试扩大机器学习候选基因范围换GWAS人群或调整疾病定义6.2 我的一些独家心得第一WGCNA的模块颜色只是个标签别纠结于颜色本身。关键是模块与表型的关联以及模块内基因的通路富集。如果某个模块与疾病强相关但富集出来的通路全是核糖体或线粒体那要警惕这是不是平台批次效应导致的伪关联。我通常会在富集结果里加一个对照用无关的基因集跑一遍同样的富集看是不是也有类似结果。第二机器学习筛选到的基因一定要做表达方向的一致性检查。比如某基因在疾病组上调那么MR结果里的基因表达增加方向也应与疾病风险增加方向一致这才符合生物学逻辑。如果表达谱说上调MR却显示表达升高是保护性的那就要想想是不是用了不同组织的eQTL导致的——基因在不同组织的表达调控模式可能不同甚至方向相反。第三MR验证部分如果某个基因的可用工具变量很少比如只有1个SNP那么证据强度会大打折扣。我会建议在这种情况下列出方向一致但没有统计学意义的结果作为“趋势性证据”不要强行下因果结论。第四整套流程的时间成本比单跑一个分析高得多。我建议你在项目开始前就把这部分工作写进研究方案的“数据挖掘与验证策略”里一边跑一边记录参数最后写论文时能省很多事。最后再分享一个小技巧分析过程中所有中间文件都要保留好版本。我习惯用一个文件夹按步骤存放比如WGCNA_output、ML_output、MR_output代码里用相对路径设置工作目录为项目根目录。这样万一哪个环节有改动不需要从头跑改一步保存一步即可。这套流程跑通一遍之后换一个疾病、换一套数据你就是批量生产高质量生信结论的节奏了。