简介这份资源是面向生物信息学与统计分析学习者的R语言实战项目聚焦用因果推断方法评估吸烟对健康肺组织中基因表达的短期因果效应适合具备一定R基础、希望理解观察性数据因果分析流程的中高级读者。压缩包共45个文件、约3.1MB包含19个png结果图、13个txt文本、8个xlsx表格、2个csv数据、1个pdf文档及R脚本等覆盖数据表、分析结果与可视化输出便于对照复现。项目涉及工具变量、倾向得分匹配、DAG等因果推断思路并结合RNA-seq差异表达、质量控制与通路富集等环节展示了从原始数据到因果结论的完整分析链路。目前已有148人学习下载读者可借此掌握R包在因果推断与基因表达分析中的实际用法理解混杂控制与结果验证的关键细节为吸烟相关疾病机制研究提供可借鉴的分析框架。1. 因果推断做基因表达为什么相关性研究总在“翻车”你手头有一堆吸烟者和非吸烟者的肺组织转录组数据跑完差异表达分析得到几百个“显著”基因。然后呢你大概率会写一句“吸烟与这些基因的表达变化相关”。但审稿人一句“相关不等于因果”就能把你堵回去。更麻烦的是吸烟者往往年龄更大、饮酒更多、BMI更高这些混杂因素全都会影响基因表达。你看到的差异可能根本不是烟草烟雾造成的。这就是为什么需要因果推断。它要回答的不是“吸烟的人哪些基因变了”而是“如果让同一个人从不吸烟变成吸烟他的基因表达会怎么变”。前者是观测关联后者是反事实干预。在肺组织这种无法做随机对照试验的场景里因果推断提供了一套数学框架把观测数据“伪装”成接近随机实验的结构。这个项目做的事情就是把这套框架落到吸烟与健康肺组织基因表达的短期因果效应评估上。搞懂它你就知道为什么单纯跑个limma或DESeq2远远不够以及因果推断到底在算什么。2. 因果推断的底层逻辑从“看到什么”到“如果怎样”2.1 潜在结果框架每个样本都有两个平行宇宙因果推断的核心是潜在结果模型。对每个肺组织样本定义两个潜在结果Y(1) 表示这个人在吸烟状态下的基因表达Y(0) 表示同一个人在不吸烟状态下的基因表达。我们真正想估计的是个体因果效应 Y(1) - Y(0)。但现实中每个人要么吸烟要么不吸烟你只能观测到一个潜在结果另一个永远缺失。这就是因果推断的根本难题。常见的实现思路是把问题转化为估计平均因果效应ATE或处理组的平均因果效应ATT。在基因表达场景下通常关注的是对于实际吸烟的那群人吸烟对他们的基因表达造成了多大改变。这个量叫ATT它比ATE更贴近“吸烟者如果戒烟会怎样”的临床问题。2.2 混杂控制为什么回归调整不够用你可能会想把年龄、性别、BMI作为协变量放进线性模型不就控制了混杂吗问题在于回归调整假设了正确的函数形式。如果年龄和基因表达是非线性关系或者吸烟概率本身依赖于年龄倾向性得分不恒定回归系数就会有偏。更隐蔽的是如果存在未观测混杂——比如遗传背景同时影响吸烟行为和基因表达——任何基于观测协变量的调整都无法消除偏倚。因果推断的做法不同。它不直接建模Y对T和X的回归而是先建模处理分配机制。具体来说估计每个样本的倾向性得分在给定协变量X的条件下该样本接受处理吸烟的概率。然后通过匹配、加权或分层构造一个“伪总体”使得在这个伪总体里处理组和对照组的协变量分布是平衡的。这样一来处理分配就近似随机了。2.3 识别假设三个不能少的条件要让因果效应从观测数据中可识别必须满足三个假设。第一条件可交换性给定协变量X处理分配与潜在结果独立。换句话说所有同时影响吸烟和基因表达的因素都被X捕获了。第二正值性对任意X接受处理的概率在0到1之间不能有某些协变量组合下所有人都吸烟或都不吸烟。第三一致性观测到的结果等于对应处理下的潜在结果。这三个假设里条件可交换性最要命因为它本质上不可检验。你只能通过敏感性分析来评估违反它的后果。注意在基因表达因果推断中条件可交换性尤其脆弱。因为基因表达受遗传变异调控而遗传变异又可能影响吸烟行为基因-环境相关性。如果协变量里没有包含这些遗传因素估计就是有偏的。3. 从倾向性得分到因果效应源码级拆解3.1 倾向性得分估计逻辑回归还是机器学习项目里估计倾向性得分的常见做法是逻辑回归把吸烟状态作为因变量年龄、性别、测序批次、细胞类型比例等作为自变量。逻辑回归的好处是可解释、稳定但在高维协变量场景下容易欠拟合。更灵活的做法是用梯度提升机或随机森林来估计倾向性得分但要注意过拟合会导致极端权重。# 倾向性得分估计逻辑回归与GBM对比 library(gbm) # 逻辑回归 ps_logit - glm(smoking ~ age sex bmi batch cell_type_pct, data covariates, family binomial(link logit)) ps_logit_pred - predict(ps_logit, type response) # 梯度提升机 ps_gbm - gbm(smoking ~ age sex bmi batch cell_type_pct, data covariates, distribution bernoulli, n.trees 5000, interaction.depth 3, shrinkage 0.001, cv.folds 5) ps_gbm_pred - predict(ps_gbm, type response, n.trees 3000)这段代码里smoking是二值处理变量age、sex、bmi、batch、cell_type_pct是协变量。逻辑回归用glm拟合type response输出概率值。GBM用gbm函数distribution bernoulli指定二分类cv.folds 5做交叉验证选树数。关键参数是interaction.depth控制树的分裂深度深度越大越容易过拟合倾向性得分。我一般会先用逻辑回归跑一版再用GBM对比如果两者估计的倾向性得分分布差异很大说明模型设定有问题。3.2 逆概率加权给每个样本重新算“话语权”得到倾向性得分后逆概率加权IPW是最直接的因果效应估计方法。它的逻辑是如果一个吸烟者被观测到的概率很低倾向性得分小那他就代表了很多和他协变量相似但没被观测到的吸烟者所以给他更大的权重。权重公式是吸烟者权重为1/ps非吸烟者权重为1/(1-ps)。# 逆概率加权计算 weights - ifelse(covariates$smoking 1, 1 / ps_logit_pred, 1 / (1 - ps_logit_pred)) # 检查权重分布 summary(weights) # 截断极端权重避免方差爆炸 weights_trunc - pmin(weights, quantile(weights, 0.99)) # 加权后的基因表达差异 library(survey) design - svydesign(ids ~1, weights ~weights_trunc, data expr_data) ate_estimate - svyglm(gene_expression ~ smoking, design design)ifelse那行根据吸烟状态分配权重。summary(weights)用来检查权重是否有极端值如果最大权重超过10通常需要截断。pmin把权重限制在99%分位数以内这是血泪经验——不截断的话一个倾向性得分0.01的样本权重会到100直接主导整个估计。svydesign和svyglm来自survey包专门处理加权数据能给出正确的标准误。3.3 协变量平衡诊断加权后到底平不平衡加权之后必须检查协变量平衡。如果加权后处理组和对照组的协变量均值仍有显著差异说明倾向性得分模型没拟合好。常用指标是标准化均值差SMD一般要求加权后SMD小于0.1。# 标准化均值差计算 smd_before - sapply(covariates[, c(age, sex, bmi, batch)], function(x) { mean(x[covariates$smoking 1]) - mean(x[covariates$smoking 0]) }) / sapply(covariates[, c(age, sex, bmi, batch)], sd) smd_after - sapply(covariates[, c(age, sex, bmi, batch)], function(x) { weighted.mean(x[covariates$smoking 1], weights_trunc[covariates$smoking 1]) - weighted.mean(x[covariates$smoking 0], weights_trunc[covariates$smoking 0]) }) / sapply(covariates[, c(age, sex, bmi, batch)], sd)smd_before计算加权前的标准化均值差smd_after计算加权后的。分母用原始数据的标准差这是标准做法。如果某个协变量的SMD在加权后反而变大了说明该协变量的倾向性得分模型设定有误可能需要加入交互项或非线性项。这一步是因果推断里最容易被跳过、也最容易翻车的地方。4. 避坑与排查因果推断在基因表达里的五个常见问题4.1 倾向性得分重叠不足现象加权后权重分布极端最大权重超过50估计结果标准误巨大。原因某些协变量组合下吸烟者或非吸烟者几乎不存在导致倾向性得分接近0或1。比如年轻非吸烟者样本极少他们的倾向性得分接近0加权后权重爆炸。解决检查倾向性得分分布如果重叠区域太窄考虑限制研究人群到重叠区域或者用稳定权重stabilized weights替代原始IPW权重。4.2 协变量平衡诊断被忽略现象加权后基因表达差异显著但审稿人要求提供平衡诊断发现年龄的SMD仍有0.3。原因倾向性得分模型只包含了线性项但年龄与吸烟概率是非线性关系比如U型导致模型欠拟合。解决在倾向性得分模型中加入年龄的样条项或分段项重新估计并检查SMD。我一般会强制要求所有协变量加权后SMD小于0.1才继续下游分析。4.3 基因表达数据未做批次校正现象因果效应估计出大量基因显著但火山图显示差异基因与批次高度相关。原因测序批次同时影响基因表达和吸烟状态如果吸烟样本和对照样本在不同批次测序批次就成了混杂因素。如果没放进倾向性得分模型估计有偏。解决把批次作为协变量放入倾向性得分模型或者先用ComBat等方法校正批次效应再用校正后的表达数据做因果推断。注意校正后的数据不能直接用于差异表达分析但可以用于因果效应估计。4.4 用错误的标准误做推断现象加权后回归的p值极小但置换检验显示结果不显著。原因IPW估计的标准误需要考虑权重估计的不确定性。如果直接用加权回归的标准误会低估真实标准误导致假阳性。解决用bootstrap重抽样估计标准误或者用survey包的svyglm它内部考虑了权重。更严谨的做法是双重稳健估计同时建模倾向性得分和结果回归。4.5 把ATT当成ATE解释现象估计出的因果效应被解释为“所有人吸烟的效应”但实际估计的是“吸烟者如果戒烟的反事实效应”。原因IPW默认估计的是ATE但加权到处理组协变量分布上估计的是ATT。两者在存在异质性处理效应时不同。解决明确研究问题。如果关注的是吸烟者的戒烟效应用ATT如果关注的是人群干预效应用ATE。在基因表达场景下ATT通常更有临床意义因为吸烟者才是干预对象。5. 设计取舍为什么不用孟德尔随机化以及敏感性分析的边界因果推断在基因表达里的应用绕不开一个根本矛盾条件可交换性假设不可检验。你永远不知道是否遗漏了同时影响吸烟和基因表达的混杂因素。孟德尔随机化用遗传变异作为工具变量天然规避了这个问题但它要求遗传变异只通过吸烟影响基因表达这个排他性假设在肺组织里同样难以满足——很多遗传变异直接调控基因表达。这个项目的设计取舍在于它选择了倾向性得分框架而不是工具变量法。原因很实际——肺组织转录组数据里吸烟状态是明确可观测的处理变量而遗传数据往往缺失或维度太高。倾向性得分框架能利用现有的临床协变量快速给出因果效应的点估计和区间估计。代价是它依赖于无未观测混杂的假设而这个假设在生物学上几乎不可能完全成立。我一般会这么理解倾向性得分方法给出的是一个“在观测协变量范围内”的因果效应估计。它的价值不在于证明因果而在于把分析从“相关”推进到“在控制可观测混杂后的因果效应”。如果你要做更严格的因果推断需要结合敏感性分析——比如计算E值评估需要多强的未观测混杂才能把估计效应解释掉。如果E值很小说明结果脆弱如果E值很大说明即使有未观测混杂也很难完全推翻结论。从那以后我每次跑因果推断都强制走一遍先检查倾向性得分重叠再诊断协变量平衡最后做敏感性分析。这三步缺一个结果都不敢往外报。希望帮到你。本文还有配套的精品资源点击获取