★ FEATURED ARTICLE 简介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改成directdagitty 会把中介加进调整集——但那时候你得到的就不再是「吸烟对表达的总影响」而是「扣掉炎症与细胞组成之后剩下的直接影响」两者是完全不同的问题报告时必须写明白。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 是具体的基因每个基因一个模型必须逐个拟合。这个不对称性直接决定了流水线的结构环节依赖结局拟合次数可否复用倾向评分模型否1全矩阵复用逆概率加权估计否0纯加权计算全矩阵复用g-computation 结局模型是2 × 基因数不可复用AIPW 校正项是2 × 基因数不可复用理解这张表之后很多工程上的取舍就自然了如果样本量很大、基因很多纯 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 值带着走等发现问题时已经写进结论了。这份资源的价值也在这里它把这条容易走偏的链路固化成了一套可复查的步骤而不是一组可以直接抄的数字。想动手的人从倾向评分那一段开始对照自己的数据改会比从头搭一遍省很多时间。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
