吸烟肺组织转录组的因果推断：DAG、倾向评分与FDR敏感性分析

简介：SmokingGeneExpression 面向生物信息学、统计分析与因果推断学习者，围绕吸烟对健康肺组织基因表达的短期因果效应，提供一套可复现的 R 语言项目材料。压缩包共 45 个文件，约 3.1MB，包含 19 个 png 结果图、13 个 txt 数据/说明、8 个 xlsx 表格、2 个 csv、1 个 pdf 及 R 脚本等，覆盖样本 PCA、基因表达表、差异基因计数与辅助数据；目录按数据、结果和脚本组织，便于按模块检索。项目涉及工具变量、倾向得分匹配、DAG 等因果识别思路，并结合 RNA-seq 预处理、质量控制、映射计数和差异表达分析，帮助读者理解从数据整理到结果解释的完整流程。已有 148 人学习，适合作为因果推断与基因表达分析交叉练习的参考案例，也可用于复现实验、拆解分析步骤和积累 R 语言生物信息学实践，并提升分析复现与结果解释的可靠性。

1. 相关性的天花板：吸烟与肺组织基因表达之间隔着三重识别难题

把吸烟者和非吸烟者的肺组织转录组并排跑一次差异表达，几乎总能得到一份长长的基因清单。这份资源要把问题往回退一步：不问「哪些基因和吸烟同时出现」，而问「同一个人的吸烟状态被改变，表达会往哪边走」。因果推断在这里只是换个名字说相关系数还不够，它得先证明自己配得上「因果」二字。R 在这条链路里也不只是画图——从倾向评分到多重检验校正，整条流程都落在同一套脚本生态里。值得拆开看，是因为这类分析最容易在三个地方悄悄失效：混杂没调干净、权重把样本量吃光、把一个横断面的短期效应讲成了终身结论。

2. 把混杂画成图：DAG、后门准则与「哪些变量碰不得」

差异表达分析把「吸烟」当成一个分组标签，因果推断把它当成一次干预。这个视角切换的代价是，你必须先写下自己的因果假设，而不能让模型自己去猜。DAG 就是这份假设的书面形式：哪些箭头存在、哪些箭头不存在，一旦画出来，调什么、不调什么就不再是主观选择，而是图上的可计算问题。

2.1 潜在结局与三个识别假设：可交换性到底在赌什么

潜在结局框架的写法很朴素：每个样本有两个潜在表达值，Y¹ 表示「如果这个人吸烟」，Y⁰ 表示「如果这个人不吸烟」，我们想看的是 E[Y¹] − E[Y⁰]。问题在于每个人只能落在其中一个世界里，观测到的永远只有 Y = Y^X。所以估计量能成立，靠的是三个假设撑住。

第一个是一致性，意思是「观测到的暴露水平」和「潜在结局上标记的暴露水平」必须对得上。这条在转录组数据里比在临床试验里麻烦得多：吸烟不是一个良定义的二元处理，剂量、年限、最后一支烟到采样的间隔都会改变暴露的实际含义。常见的做法是把暴露二值化成「采样前 12 个月内是否规律吸烟」，并在方法里把这个窗口写死。窗口越模糊，一致性假设越容易被质疑，后面所有数字的可信度也跟着打折。

第二个是条件可交换性，也就是在给定一组协变量 Z 之后，吸烟组和非吸烟组在潜在结局上没有系统差异。这是整套方法里唯一无法用数据验证的假设，只能用领域知识和敏感性分析去支撑。我一般会这么理解：它不是「证明没有混杂」，而是「声明我把哪些混杂处理掉了，剩下的留给偏倚分析去量化」。

第三个是正性，即每个协变量组合下，既有人吸烟也有人不吸烟。听起来是技术细节，实际是转录组因果分析里翻车最多的地方：如果某一批样本几乎全是吸烟者，那么这一批里根本不存在可比较的反事实，倾向评分会逼近 0 或 1，权重随之发散。正性一旦违背，再精巧的估计量也救不回来，只能缩小目标人群或做限制。

还有一个常被忽略的前提是无干扰性，即一个人的暴露状态不会影响另一个人的表达。在个体层面的组织样本里这一条基本成立，但在同一批处理、同一张芯片上的样本之间，技术层面的「干扰」是真实存在的，这也是批次必须进图的原因。

2.2 后门路径的枚举：从 DAG 到最小充分调整集

把假设画成图之后，「该调什么」就变成一个图论问题。后门准则说的是：调整集 Z 必须阻断从暴露到结局的所有后门路径，同时不能包含暴露的后代。前者防止漏调，后者防止过度控制。这两个方向的错误在转录组分析里都极其常见。

下面是一张最小可用的肺组织转录组 DAG，用 R 里的 dagitty 写出来：

library(dagitty) g <- dagitty('dag { smoke [exposure] expr [outcome] age -> smoke ; age -> expr sex -> smoke ; sex -> expr bmi -> smoke ; bmi -> expr batch -> smoke ; batch -> expr purity -> smoke ; purity -> expr smoke -> expr smoke -> inflam -> expr smoke -> cellcomp -> expr }') # 枚举所有满足后门准则的最小调整集 adjustmentSets(g, exposure = "smoke", outcome = "expr", effect = "total")

这段代码的价值不在输出本身，而在它逼你把「为什么调这个变量」说清楚。effect = "total"表示要估计总效应，此时inflam（炎症通路活性）和cellcomp（细胞组成）作为中介变量不进调整集，因为它们处在smoke -> ... -> expr的路径上。如果你把effect改成"direct"，dagitty 会把中介加进调整集——但那时候你得到的就不再是「吸烟对表达的总影响」，而是「扣掉炎症与细胞组成之后剩下的直接影响」，两者是完全不同的问题，报告时必须写明白。

purity（组织纯度）和batch是两个容易争议的节点。组织纯度既是样本处理的产物，也会影响可测量的表达谱，同时在某些队列里与吸烟状态相关，进调整集是合理的。批次同理，但如果批次与暴露高度共线，调整会撞上正性违背，这时候不是「调不调」的问题，而是这批数据不该被用来回答这个问题。

我通常会把adjustmentSets()的结果和领域知识对一遍，而不是照单全收。图是我画的，图错了，后门准则只会帮我把错误执行得更彻底。

2.3 中介与对撞：两类「看起来该调、其实不能调」的变量

中介变量是第一个陷阱。吸烟会诱导气道炎症，炎症通路活化又会改变一大批基因的表达。如果你把炎症相关的通路活性放进调整集，问的就不再是「吸烟会不会改变表达」，而是「在不改变炎症的前提下吸烟会不会改变表达」——这个反事实在现实里几乎不存在。过度控制偏倚的机制就是这样：你从总效应里切掉了一段真实存在的通路，剩下的估计量既不是总效应也不是干净的直接效应，而是一个没有明确解释的混合体。

这里有一个判断标准：这个变量是不是在暴露发生之后才被改变的？如果是，它大概率是中介或者中介的下游，默认不该进总效应的调整集。要做中介分析，就得换一套更强的识别假设，而不是简单地在回归里多加一列。

对撞变量是第二个陷阱，也是更隐蔽的一个。对撞的意思是，某个变量同时被暴露和结局（或结局的原因）影响。按对撞分层，会在暴露与结局之间凭空造出一条关联。举个具体的场景：假设是否进入某项后续检查同时受吸烟状态和某些基因表达水平影响，而你为了「控制数据质量」只保留了做过这项检查的样本。这一步筛选就把对撞条件化了，剩下的样本里，吸烟与表达之间会出现一个在真实人群中不存在的关联，方向还可能完全反掉。这类偏倚在生信流程里极难察觉，因为筛选步骤通常被写在「数据清洗」里，没人会把它和因果结构联系起来。

还有一类变量介于两者之间：细胞类型组成。吸烟会改变肺组织的细胞构成（比如免疫细胞浸润比例上升），而不同细胞类型的表达谱差异巨大。所以细胞组成既可能是混杂（如果它影响采样与暴露的关联），也可能是中介（如果它确实是吸烟导致的下游改变）。这个岔路口值得单独拿出来讲，后面第 5 章会展开。

3. 估计量的取舍：倾向评分与 g-computation 在表达矩阵上怎么落地

有了调整集，剩下的问题是用什么估计量把它落实。转录组数据有几个很实际的特征：样本量通常在几百这个量级，协变量维度不高但彼此相关，结局维度是两万左右。这三条决定了估计量不能照搬流行病学教科书里的默认选项。

3.1 倾向评分的直觉：每个样本代表多少人

倾向评分 e(Z) = P(X = 1 | Z) 的定义很直白：在给定协变量的情况下，这个人吸烟的概率。它的作用是把一堆协变量压缩成一个标量，让「调整」变成「加权」。直观理解是，逆概率加权给每个样本发一张票，票的权重表示「这个人在目标人群里代表多少人」。一个吸烟概率只有 0.1 却实际吸烟的人，代表了 10 个像他这样的人，所以权重要放大到 10 倍。

这种加权在协变量平衡上很有效，但代价是方差。权重分布的尾部决定了估计量的稳定性，所以任何一份严肃的 IPW 分析都应该报告有效样本量：

ESS = (Σw)² / Σw²

ESS 掉到原始样本量的一半以下时，估计量的方差会明显膨胀。我在实际操作里会把 ESS 当成一个硬性诊断指标，而不是等结果跑完再去解释为什么置信区间那么宽。

倾向评分模型本身也有选择：logistic 回归简单、可解释、系数能看，但对非线性关系不敏感；梯度提升或超级学习器拟合能力更强，但容易过拟合，在小样本上反而不稳。常见的折中是先用 logistic 回归加必要的样条项，再用平衡诊断（标准化均值差）验证加权后的协变量是否真的对齐，而不是看模型 AUC 有多高。AUC 高不代表平衡好，这一点很多人会搞混。

3.2 为什么两万个基因只需要一套权重

这是这套分析里最值得讲清楚的一个设计点，也是计算成本能压下来的关键。

倾向评分模型是 P(X | Z)，它只涉及暴露和协变量，完全不碰任何基因的表达值。换句话说，无论你分析的是第一个基因还是第一万个基因，暴露模型都是同一个。所以最优做法是：拟合一次倾向评分，得到一组权重，然后在整个表达矩阵上复用。

结局模型则相反，E[Y | X, Z] 里 Y 是具体的基因，每个基因一个模型，必须逐个拟合。这个不对称性直接决定了流水线的结构：

理解这张表之后，很多工程上的取舍就自然了：如果样本量很大、基因很多，纯 IPW 的加权均值差计算可以完全向量化，成本几乎可以忽略；而 AIPW 因为要给每个基因拟合两个结局模型，计算量会随基因数线性增长。这也是为什么在探索阶段大家更愿意先跑 IPW 拿到一个粗筛结果，再对显著的那部分基因补 AIPW 校正。

顺带说一句，这个不对称性还有一个统计上的好处：暴露模型只拟合一次，意味着权重只估计一次，不会因为基因不同而引入额外的模型选择变异性。用一个固定的权重去扫全矩阵，至少保证了所有基因是在同一个「伪人群」里比较的。

3.3 代码骨架：从 PS 到 AIPW 的最小可运行链路

下面这段代码把上面的思路落成一个可以跑通的链路。协变量表cov_df每行一个样本，表达矩阵expr_mat每行一个基因，列顺序与协变量表一致。

# ---- 1. 暴露模型：一次拟合，全基因复用 ---- ps_fit <- glm(smoke ~ age + sex + bmi + purity + batch, family = binomial(link = "logit"), data = cov_df) cov_df$ps <- predict(ps_fit, type = "response") # 截尾：把极端倾向评分压到 [0.01, 0.99]，避免单个样本权重爆炸 cov_df$ps_t <- pmin(pmax(cov_df$ps, 0.01), 0.99) # ---- 2. 稳定化权重 ---- # 分子用边际暴露概率，防止权重整体漂移导致方差虚高 p_marg <- mean(cov_df$smoke) cov_df$w <- ifelse(cov_df$smoke == 1, p_marg / cov_df$ps_t, (1 - p_marg) / (1 - cov_df$ps_t)) # 有效样本量：权重越不均匀，ESS 越小，方差越大 ess <- sum(cov_df$w)^2 / sum(cov_df$w^2) cat(sprintf("n = %d, ESS = %.1f\n", nrow(cov_df), ess)) # ---- 3. 加权均值差：估计单个基因的 ATE ---- ate_ipw <- function(y, x, w) { # y 应先做 log2(CPM + 1) 转换并做分位数标准化 m1 <- sum(w[x == 1] * y[x == 1]) / sum(w[x == 1]) m0 <- sum(w[x == 0] * y[x == 0]) / sum(w[x == 0]) m1 - m0 } ate_all <- apply(expr_mat, 1, ate_ipw, x = cov_df$smoke, w = cov_df$w)

ps_t这一步需要解释清楚：截尾不是「修正」正性违背，而是把极端权重压下来换取方差稳定，代价是引入一点点偏倚。阈值取 0.01/0.99 是常见起点，实际应该看倾向评分的分布形状决定，如果两端本来就很稀疏，压到 0.05/0.95 也合理。

ate_ipw是加权均值差，本质上是 Horvitz-Thompson 型估计量的稳定化版本。它的假设是「加权后两组可比」，但没有用结局模型做任何校正，所以在协变量与结局关系复杂时会有残余偏倚。补一个 AIPW 版本可以拿到双重稳健性：

ate_aipw <- function(y, d, ps) { # 结局模型：分别在两组内拟合，再对全体样本做预测 f1 <- lm(y ~ age + sex + bmi + purity + batch, data = d, subset = d$smoke == 1) f0 <- lm(y ~ age + sex + bmi + purity + batch, data = d, subset = d$smoke == 0) m1 <- predict(f1, newdata = d) m0 <- predict(f0, newdata = d) # 双重稳健：结局模型或倾向评分模型有一个对，估计量就一致 mean(m1 - m0 + d$smoke * (y - m1) / ps - (1 - d$smoke) * (y - m0) / (1 - ps)) } ate_all_aipw <- apply(expr_mat, 1, ate_aipw, d = cov_df, ps = cov_df$ps_t)

AIPW 的公式看着吓人，拆开就三段：第一段m1 - m0是 g-computation 的预测差，第二段和第三段是残差校正项，用倾向评分的倒数给残差加权。双重稳健的意思是，f1/f0设定正确时，即使倾向评分模型错了估计也一致；反过来倾向评分对、结局模型错，估计仍然一致。两个都错才会失效，这比单靠一种模型实在得多。

代价也很明确：每个基因要拟合两个线性模型，两万个基因就是四万次拟合。在几百个样本的规模上，R 里的lm会慢，但能忍受；如果基因数再往上走，就得换limma的矩阵化拟合或者speedglm，把循环压掉。

还有一点容易被忽略：ps必须和w使用同一个截尾后的版本，否则 AIPW 的校正项和加权项不在同一个尺度上，估计量会偏。这类细节不会报错，只会让结果默默错掉。

4. 两万个基因同时检验：FDR、阴性对照与未测混杂的敏感性分析

单个基因的因果效应估计出来之后，问题并没有结束。多重检验、残余混杂、模型设定错误，这三件事会在规模效应下同时放大。一份可信的分析必须在这三处都给出交代，而不是只报一个 q 值。

4.1 BH 校正解决的是规模，不是偏倚

两万个基因同时检验，即使全部为零效应，按 α = 0.05 也会有一千个「显著」。Benjamini-Hochberg 过程把假发现率的期望控制在设定水平，这是标准做法：

res$p_value <- 2 * pnorm(-abs(res$ate / res$se)) res$q_value <- p.adjust(res$p_value, method = "BH") # 报告时同时给出效应量与校正后显著性 sig <- subset(res, q_value < 0.05 & abs(ate) > 0.2)

这里有两个容易被忽略的点。第一，p 值来自渐近正态近似，在 ESS 掉到三位数时这个近似很不靠谱，置换检验会更稳：

set.seed(2024) perm_null <- replicate(200, { # 打乱暴露标签，保留协变量与权重的联合结构 idx <- sample(nrow(cov_df)) ate_ipw(expr_mat[1, ], cov_df$smoke[idx], cov_df$w) })

第二，也是更根本的：BH 控制的是「零假设为真时的误报比例」，它假设你的 p 值本身是有效的。如果因果估计存在系统性偏倚，校正后的 q 值只会让偏倚变得更「显著」。FDR 处理的是规模问题，不是识别问题，这两件事必须分开看。

4.2 阴性对照：用「不该有信号」的地方校验流程

阴性对照是这套流程里性价比最高的验证手段，因为它不需要额外的数据，只需要你事先想清楚「哪里本来就不该有信号」。

阴性对照结局的做法是：挑一组已知与吸烟没有因果关系的基因，比如某些管家基因或者与暴露完全无关的组织特异性基因，跑完整条流水线。如果它们在 q < 0.05 上大量显著，说明有残余混杂或技术偏倚在推动整个分布，这时候主结果里的信号也要打折扣。这个诊断的说服力很强，因为它检验的不是「我的方法有没有检出真信号」，而是「我的方法会不会制造假信号」。

阴性对照暴露的思路类似：找一个理论上不可能影响基因表达、但在数据里可能与表达相关的变量，比如样本采集的时间编号或者某个随机分配的技术参数。如果它跑出了「显著效应」，几乎可以肯定流程里混进了虚假关联。

我一般会把这组对照放在主分析之前跑。先确认系统是干净的，再看主结果，比事后解释一堆异常要省力得多。

4.3 E-value：未测混杂需要多强才能掀翻结论

条件可交换性无法用数据验证，但可以量化它有多脆弱。E-value 回答的问题是：需要多强的未测混杂（与暴露和结局的关联强度都达到这个水平）才能把观测到的效应完全解释掉。E-value 越大，结论越稳。

对连续结局，需要先把标准化均值差 d 转成近似的风险比：

evalue_from_d <- function(d) { rr <- exp(0.91 * abs(d)) # 标准化均值差到风险比的常用近似 if (rr < 1) rr <- 1 / rr rr + sqrt(rr * (rr - 1)) # E-value 闭式解 } res$evalue <- evalue_from_d(res$ate / res$sd_resid)

解读方式要小心：E-value 只是一个「脆弱性指标」，不是「混杂存在的概率」。一个 E-value 为 1.5 的结论，意味着只需要一个中等强度的未测混杂就能翻盘，这种结果在报告里应该被明确标成「不稳健」，而不是和 E-value 为 4 的结果混在一起讲。0.91这个系数来自标准化均值差与风险比的近似转换，量纲上没问题，但前提是结局近似正态，表达值做过 log 转换之后这个前提大体成立。

把 FDR、阴性对照、E-value 三者放在一起看，你会得到一张比「显著基因列表」信息量大得多的表：哪些基因既有统计显著性又有足够大的效应量，哪些基因的结论经得起中等强度混杂的冲击。这份表才是这份资源真正值得深挖的部分。

5. 排查清单：五个把因果结论做废的典型现场

下面这五条都是我在实际流程里反复遇到的。每一条都按「现象 → 原因 → 解决」写，方便直接对照自己跑出来的结果。

5.1 权重极端与 ESS 崩塌

现象：权重最大值到了几十，某个基因的 ATE 在去掉一两个样本后符号翻转，置信区间宽得没有意义。

原因：倾向评分的尾部重叠很差。常见诱因是把「吸烟包年数」这类直接定义暴露的变量塞进暴露模型，或者协变量里有一个与暴露几乎一一对应的变量（比如某个只在吸烟者中采集的样本批次）。这两种情况都会把 ps 推到 0 或 1 附近。

解决：先画两组倾向评分的叠加密度图，看重叠区域。然后做两件事——把权重截尾到 [0.01, 0.99]，报告截尾前后的 ESS 变化；如果 ESS 仍然低于样本量的一半，改用重叠权重（overlap weights），它把目标人群收缩到「两组都可能出现」的区域，代价是改变了估计目标，需要在方法里写明。

5.2 细胞组成的岔路口

现象：调整细胞类型比例之后，一批免疫相关基因的效应几乎消失，而调整之前它们是全矩阵里最显著的。

原因：bulk 组织的表达值是「细胞组成」和「细胞内表达」的混合。吸烟会改变肺组织的细胞构成，所以组成变化既是真实的效应通路，也是观测表达的驱动因素。把它当混杂调掉，等于把一部分真实效应切走；不调，又会把组成差异误读成细胞内的调控变化。

解决：先做去卷积，拿到每个样本的细胞类型比例，然后分两条线报告：一条是总效应（不调组成），一条是调整组成后的效应（近似「细胞内表达变化」）。两条线的差异本身就是信息，它告诉你这个基因的效应有多少来自组织构成重塑。只报其中一条，都是在隐藏关键结论。

5.3 批次与暴露共线

现象：倾向评分模型里批次变量的系数很大，加入批次后 ESS 从 400 掉到 150，一批基因的效应方向整体反转。

原因：吸烟样本集中在某几个批次里。这可能是样本收集年份不同、也可能是样本来源不同。此时批次和暴露在数据里几乎不可分，正性假设在批次层面被违背。

解决：先做 batch × smoke 的列联表，看每个批次里两组的样本数。如果某些批次只有一组，调整解决不了问题，只能限制到两组都有的批次，或者把批次作为分层变量做分层估计再合并。硬调一个共线的协变量，得到的不是「控制」，是把暴露的变异也一起控掉了。

5.4 选择性推断：先筛后估

现象：先在全部基因上做关联筛选，再对筛出来的基因做因果估计，q 值一片漂亮。

原因：同一份数据被用了两次——一次用来选基因，一次用来做推断。筛选步骤本身消耗了数据里的信息，后续的 p 值不再服从零假设下的均匀分布，FDR 控制失效。

解决：样本分裂是最直接的办法，一半选基因、一半估效应。样本量不够时，至少做一次独立性检查：把数据随机分成两半，各自跑完整的筛选加估计，看两次的显著基因重合度。重合度低就说明结果高度依赖具体的筛选路径。

5.5 把「FDR 显著」当成机制证据

现象：报告里写「吸烟通过某炎症通路导致这批基因表达改变」。

原因：主效应显著只说明暴露和结局之间有关联，完全不涉及中介。要声称「通过某通路」，需要正式的中介分析，而中介分析要求跨世界反事实假设（比如「吸烟但通路活性被固定在不吸烟水平」这种情况是否存在），横断面数据几乎无法满足。

解决：把结论降级。能说的是「这批基因的表达变化与炎症通路活性变化在方向上一致」，不能说「由炎症通路介导」。真要做中介，就得换数据设计，或者至少用敏感性分析把中介假设的脆弱程度量化出来。

6. 短期因果的边界：这份资源能回答什么、不该被拿来回答什么

写到这里，有一个判断需要摆到明面上：这套流程估计的是「短期」因果，而这个「短期」不是修辞，是数据结构决定的。暴露是采样时点的近期吸烟状态，结局是同一时点的转录状态，两者之间没有时间维度。所以它能回答的是「当前吸烟状态与当前表达之间的差异，在多大程度上不能被已测混杂解释掉」，而不是「吸烟若干年后表达会怎么变」。

把能答和不能答列清楚，比任何模型细节都重要：

这张表也解释了为什么第 4 章的敏感性分析值得花那么多篇幅。横断面数据里，条件可交换性永远是一个声明而不是一个事实，你能做的只是把「这个声明有多脆弱」量化出来。E-value 和阴性对照不能证明因果，但它们能在结论站不住的时候提前报警。

还有一个更微妙的取舍：这套流程把「总效应」作为默认目标，也就是不调整中介变量。这个选择的代价是估计量里混进了细胞组成变化、炎症通路活化等下游过程，效应的生物学解释会变得模糊；好处是它对应一个真实存在、可被干预的问题——「让人不吸烟，表达谱会怎样变化」。如果改成估计直接效应，数字会更「干净」，但那个反事实在现实里几乎无法实现。我倾向前者，因为它至少对应一个能落地的干预。

从那以后我每次跑这类脚本，都会先强制走一遍「倾向评分重叠诊断 + ESS + 阴性对照」三件套，再去看主结果。顺序反过来，人很容易被漂亮的 q 值带着走，等发现问题时已经写进结论了。这份资源的价值也在这里：它把这条容易走偏的链路固化成了一套可复查的步骤，而不是一组可以直接抄的数字。想动手的人，从倾向评分那一段开始对照自己的数据改，会比从头搭一遍省很多时间。希望帮到你。

本文还有配套的精品资源，点击获取