简介：这份资源是面向生物信息学与统计分析学习者的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值很大，说明即使有未观测混杂，也很难完全推翻结论。
从那以后我每次跑因果推断，都强制走一遍：先检查倾向性得分重叠，再诊断协变量平衡，最后做敏感性分析。这三步缺一个，结果都不敢往外报。希望帮到你。
本文还有配套的精品资源，点击获取
