简介SmokingGeneExpression 是一份面向生物信息学与统计分析学习者的 R 语言项目资料围绕吸烟对健康肺组织基因表达的短期因果效应展开。内容将因果推断思路落到基因表达场景涉及工具变量、倾向得分匹配、有向无环图等混杂控制方法以及 RNA-seq 差异表达分析、PCA 可视化与生物学解释适合希望理解观察性数据因果分析流程并借助 R 复现实战的读者。压缩包共 45 个文件约 3.1MB以 png 图表、txt 结果说明、xlsx 与 csv 数据表为主另含 R 脚本和 pdf 文档覆盖分析代码、数据样例与结果呈现等用途。已有 148 人学习可作为课程作业、课题预研或方法对照的参考。读者可从中获得较完整的分析素材按目录梳理研究逻辑、复现关键步骤并积累 R 语言实践。1. 项目定位用 R 把吸烟对肺组织基因表达的因果效应算清楚一个做肺组织转录组的同行跟我吐槽过差异表达分析跑完吸烟组和从不吸烟组之间挂着三千多个显著基因图表也漂亮结果报告交上去第一个问题就是“你怎么知道这些差异不是吸烟者普遍年龄更大、BMI 更高造成的”。相关性分析能告诉你两组人不一样但不能告诉你是吸烟把某个基因的表达抬上去的。SmokingGeneExpression 这个 R 项目要解决的正是这一层把“吸烟”当作处理变量把“某个基因的表达量”当作结局用因果推断的估计量把混杂因素压下去给出每个基因在健康肺组织里对吸烟暴露的短期因果效应估计。适合谁来读手上有队列或公开表达资源、会写基本 R 脚本、但还没系统跑过倾向评分加权和双重稳健估计的人。项目的技术栈是纯 R——数据层用 data.table 和基础的 Bioconductor 表达处理估计层主要靠 WeightIt、grf、sandwich 这几个包可视化走 ggplot2。它不训练模型、不做深度学习整套流程跑完就是一张带置信区间和 q 值的基因效应表外加几张诊断图。对做生信但被审稿意见追着问因果的人这个骨架的参考价值比多跑十个差异分析大。2. 工程结构R 脚本怎么分工、依赖怎么装、第一次怎么跑通这个项目不是包package是一套按分析阶段切开的脚本工程。这意味着它没有R CMD check那套约束好处是改起来自由坏处是环境不一致时排错特别费劲——这是后面避坑章节的重要伏笔。2.1 依赖清单哪些包是必需的哪些是可选的依赖分成三档。第一档是数据处理data.table、dplyr、tidyr、readr第二档是因果估计核心WeightIt倾向评分与权重、grf因果森林与 AIPW 的机器学习 nuisance 估计、sandwich加lmtest聚类稳健标准误第三档是多重检验与出图qvalue、ggplot2、patchwork。并行那块用futurefurrr规模化跑全基因组时会开关它。环境锁定用renv项目根目录有一份renv.lock。恢复环境的标准动作# 在项目根目录执行renv 会按 lock 文件把包版本拉回一致状态 Rscript -e if (!requireNamespace(renv, quietly TRUE)) install.packages(renv); renv::restore()这条命令的重点不是“安装依赖”而是“安装指定版本的依赖”。grf和WeightIt这两个包在版本升级时函数签名改过尤其是WeightIt的estimand参数和grf的honesty默认值跨版本跑出来的数会对不上。我第一次拿到这套脚本时偷懒没跑 restore直接用系统里现装的包结果 IPW 和 AIPW 的效应估计差了将近 8%排查了两天才发现是权重归一化方式变了。血泪经验这类因果估计工程包版本必须锁死。如果机器上装不了qvalueBioconductor 包退路是用p.adjust(method BH)两者在大多数情况下结论一致只是 q 值的具体数值略有差异报告里注明用的是哪一种即可。2.2 目录分层data、R、scripts、output 各管什么工程目录大致是这样切的路径内容说明data/raw/原始表达矩阵、表型表、协变量表只读任何脚本不得写回data/derived/标准化后的表达、构建好的分析集中间产物可全部删除重跑R/函数层io.R、preprocess.R、estimators.R、sensitivity.R、utils.R只定义函数不执行分析scripts/执行层01_prepare_data.R到05_report.R有编号按顺序跑output/tables/基因效应表、诊断统计表交付物output/figures/森林图、权重分布图、QQ 图交付物这个分层是有讲究的。R/里的函数不碰全局变量参数从外面传进去方便单独做单元测试scripts/里的脚本只负责“读什么、调哪个函数、写到哪里”。很多生信分析脚本写着写着就变成一个两千行的巨型文件中间变量到处飞改一个参数要重跑全流程这种结构就是为了防这个。常见的做法是每个scripts/开头统一 source 一遍函数层# scripts/01_prepare_data.R 的开头固定写这一块 suppressPackageStartupMessages({ library(data.table) library(dplyr) }) # 统一从项目根目录解析路径避免依赖 setwd proj_root - here::here() source(file.path(proj_root, R, io.R)) source(file.path(proj_root, R, preprocess.R)) source(file.path(proj_root, R, estimators.R)) # 全局参数集中在这一处后面所有脚本读同一个配置 cfg - list( expr_file file.path(proj_root, data/raw/expression_matrix.tsv), pheno_file file.path(proj_root, data/raw/phenotype.tsv), covar_file file.path(proj_root, data/raw/covariates.tsv), min_expr 1, # 表达过滤阈值log2 尺度下取 1 大致对应低表达噪声层 min_samples 0.3, # 一个基因至少在 30% 样本中超过阈值才保留 seed 20240101 )cfg这个列表是整条流水线的黑匣子开关。min_expr和min_samples直接决定后面有多少基因进入多重检验改动了结论就得重出seed因为因果森林内部有随机抽样不固定的话两次跑出来的效应估计在小数点后几位会飘。路径用here::here()而不是相对路径是为了让脚本从 RStudio、命令行、还是调度系统提交都能找到文件这个小习惯能省掉大量“在我机器上能跑”的扯皮。2.3 第一次跑通从 01 到 05 的执行顺序五个脚本是一条串行链但中间结果都落盘所以可以从任意一步重跑# 全流程按编号顺序执行 Rscript scripts/01_prepare_data.R Rscript scripts/02_fit_ipw.R Rscript scripts/03_fit_aipw.R Rscript scripts/04_causal_forest.R Rscript scripts/05_report.R01做数据对齐和表达预处理输出data/derived/analysis_set.rds02和03分别在同一个分析集上跑 IPW 和 AIPW各自产出一张基因级效应表04只在候选基因子集上跑因果森林用来观察异质性05把三张表按基因 ID 合并、统一多重检验校正、出图和写报告。第一次跑建议先在01里把基因数砍到几百个做冒烟测试确认整条链没断再放开全量。基因数上万的表达矩阵走完全套单核大概要几十分钟到几小时取决于样本量和是否开并行。冒烟测试这一步别省04的因果森林在小样本上最容易因为某个基因方差为零直接抛错早发现早改。3. 数据层拆解暴露、结局与协变量怎么对齐到同一张表因果推断里最容易翻车的从来不是估计量写错而是数据对齐错了。这个项目的核心数据结构只有一张宽表每一行是一个样本列包括样本 ID、吸烟暴露、几十上百个协变量、以及每个基因的表达值。所有估计函数都只认这张表保证口径一致。3.1 吸烟暴露变量的编码二值化之前的三个判断原始表型里吸烟状态通常是分好几档的从不、曾经、当前、偶尔、已戒若干年甚至还有每天支数。直接二值化会丢信息但项目走的是“当前吸烟 vs 从不吸烟”这条二值路线理由有三短期因果效应下“曾经吸烟”人群的暴露已经停止把它塞进处理组会稀释效应每天支数在多数队列里缺失率太高做连续暴露需要另一套剂量-反应模型而二值暴露能直接套用倾向评分框架诊断工具也成熟。编码这一步要把“曾经吸烟”和缺失显式剔出去而不是当成对照组# 把原始吸烟编码映射成 0/1 的处理变量并显式记录剔除原因 build_exposure - function(pheno) { stopifnot(all(c(sample_id, smoke_status) %in% names(pheno))) out - pheno %% mutate( smoke_status tolower(trimws(smoke_status)), # 只保留两种明确状态其余一律置 NA W case_when( smoke_status %in% c(current, current smoker, daily) ~ 1L, smoke_status %in% c(never, never smoker) ~ 0L, TRUE ~ NA_integer_ ), # 单独打个标记方便统计各剔除原因的数量 drop_reason case_when( is.na(W) smoke_status %in% c(former, past) ~ former_smoker, is.na(W) is.na(smoke_status) ~ missing_status, is.na(W) smoke_status ~ blank_status, is.na(W) ~ other_status, TRUE ~ NA_character_ ) ) attr(out, exposure_counts) - table(out$W, useNA always) attr(out, drop_counts) - table(out$drop_reason, useNA no) out }case_when的匹配顺序很重要%in%里放了多种写法是为了兼容不同来源的表型表真实数据里 “Current Smoker” 和 “current” 混着出现是常态。drop_reason这一列看起来多余但它是后面写方法学部分的证据——审稿人一定会问你为什么最终样本是 N你得能答出“剔除了多少既往吸烟者、多少状态缺失”。attr挂上去的计数会被01脚本打印到日志里养成这个习惯能省掉回头重算的麻烦。还有一层容易忽略的暴露的时间窗。项目说的是“短期”那就得确认表型采集时间和组织取样时间接近如果表型是问卷回顾性的“过去一年是否吸烟”而样本是五年前存的这个短期假设就不成立。数据里如果有采样日期字段务必在预处理阶段做个时间差统计。3.2 表达矩阵的预处理过滤、标准化、以及要不要去批次表达矩阵进模型前要过三道关。第一道是低表达过滤用min_expr和min_samples两个参数控制目的是减少检验数量、降低多重校正的惩罚。第二道是分布变换如果是原始 count 就走log2(CPM 1)或vst如果上游已经给了 log2 尺度就跳过。第三道是批次这一步是分水岭批次如果和吸烟状态强相关去批次会把效应一起减掉如果批次独立于暴露去不去影响都不大。保守做法是先检查批次与暴露的关联再做决定# 过滤 标准化并在返回前检查批次是否与暴露混杂 prep_expression - function(expr, meta, cfg) { mat - as.matrix(expr[, -1, drop FALSE]) rownames(mat) - expr[[1]] # 1) 低表达过滤两个条件同时满足才保留 keep - rowMeans(mat cfg$min_expr) cfg$min_samples mat - mat[keep, , drop FALSE] # 2) log2 变换已经在 log 尺度上的矩阵跳过 if (max(mat, na.rm TRUE) 50) mat - log2(mat 1) # 3) 分位数标准化抹平样本间整体分布差异 mat - limma::normalizeBetweenArrays(mat, method quantile) # 4) 批次与暴露的关联检查卡方检验p 小说明批次和吸烟纠缠在一起 if (!is.null(meta$batch) !all(is.na(meta$batch))) { tab - table(meta$batch, meta$W) p_batch - suppressWarnings(chisq.test(tab)$p.value) message(sprintf(批次-暴露关联检验 p %.4g, p_batch)) if (!is.na(p_batch) p_batch 0.05) { warning(批次与吸烟状态显著关联去批次可能吸收掉部分处理效应慎用 ComBat) } } list(expr mat, n_gene nrow(mat), n_sample ncol(mat)) }rowMeans(mat cfg$min_expr)这一行是过滤的核心表达式返回每个基因超过阈值的样本比例再和min_samples比。normalizeBetweenArrays走 quantile 方法对偏态分布的表达值比中位数标准化更稳。最值得看的是第 4 段——批次卡方检验。如果 p 值很小说明这批数据里“哪台机器测的”和“这个人吸不吸烟”高度相关此时上 ComBat 或limma::removeBatchEffect等于把处理效应的一部分当成批次削掉了后面估计出来的效应会系统性偏小。这种情况我一般不去批次而是把批次当作协变量放进倾向评分模型让加权去平衡它而不是从数据里删掉它。矩阵的样本顺序和表型的样本顺序必须严格一致这是最容易出错的地方。prep_expression返回前应该强制按meta$sample_id重排或者用stopifnot(identical(colnames(mat), meta$sample_id))做一次断言别指望两边碰巧一致。3.3 协变量集与设计矩阵共线性检查协变量集的选择原则是只放暴露之前就已经确定的变量年龄、性别、遗传背景主成分、测序深度、批次不要放可能被暴露影响的中介变量比如某些炎症指标、肺功能。放中介进去会把总效应部分遮掉估出来的是直接效应解释起来是另一回事。确定变量集之后必须检查设计矩阵的秩。协变量之间完全共线比如两个高度相关的技术指标或者某个分类变量某一水平只有一个人会让倾向评分模型直接崩掉glm可能不报错但给出荒谬的系数。# 构建协变量设计矩阵并做秩与稀疏水平检查 check_design - function(meta, covars) { f - reformulate(covars) mm - model.matrix(f, data meta) # 秩亏说明存在完全共线glm 的结果不可信 r - qr(mm)$rank message(sprintf(设计矩阵维度 %d x %d秩 %d, nrow(mm), ncol(mm), r)) if (r ncol(mm)) warning(设计矩阵秩亏存在完全共线协变量需手动剔除) # 分类变量稀疏水平某个水平样本过少会让该水平的倾向评分极端 for (v in covars) { if (is.character(meta[[v]]) || is.factor(meta[[v]])) { tb - table(meta[[v]]) if (any(tb 20)) warning(sprintf(变量 %s 存在样本量 20 的水平%s, v, paste(names(tb)[tb 20], collapse , ))) } } invisible(mm) }qr(mm)$rank小于列数就是秩亏最常见的成因是分类变量做了独热编码后某一列全是零或者两个连续变量做了线性组合。tb 20这个阈值不是硬标准20 只是个提醒线真正要看的是该水平在两组里的分布——如果某个批次只有 3 个人且全在吸烟组倾向评分在这个批次上会趋近 1加权时必然出极端权重这就是第 5 章要展开的问题。提前在这里发现比跑到加权阶段再看权重分布图省事得多。4. 因果估计核心三种估计量在单个基因上的完整链路数据对齐之后估计层的输入就是一个向量结局Y某个基因的表达值、处理W0/1、协变量矩阵X。项目给三种估计量不是让你三选一而是让它们互相印证IPW 看加权后的可比人群AIPW 用两个模型互相兜底因果森林看效应是不是在亚群里不一致。三者结论方向一致才敢往报告里写。4.1 倾向评分与 IPW把不可比的人群“造”成可比倾向评分e(X) P(W1|X)的直觉是给每个样本算一个“被吸烟的概率”然后用这个概率的倒数做权重让协变量分布在两组之间拉平。吸烟概率高的人如果实际没吸烟权重就大代表他“替”了很多吸烟者。实现上用WeightIt它比手写glm多了一层诊断# 单个基因的 IPW 估计输入是结局向量 处理向量 协变量数据框 estimate_ipw - function(y, w, x, trim 0.01) { df - data.frame(y y, W w, x) # 用 logistic 回归估倾向评分estimand ATE 表示估的是全人群平均效应 wt - WeightIt::weightit(W ~ ., data df, method glm, estimand ATE, stabilize TRUE) ps - wt$ps # 重叠性诊断倾向评分贴边说明两组人根本不可比 overlap - c(min min(ps), max max(ps), n_lt01 sum(ps 0.01), n_gt99 sum(ps 0.99)) # 权重截尾压掉极端权重的杠杆 wts - wt$weights if (!is.null(trim) trim 0) { q - quantile(wts, probs c(trim, 1 - trim)) wts - pmin(pmax(wts, q[1]), q[2]) } # 加权最小二乘给 ATErobust TRUE 拿稳健标准误 fit - lm(y ~ W, data df, weights wts) se - sqrt(sandwich::vcovHC(fit, type HC1)[W, W]) list(est coef(fit)[[W]], se se, ps ps, overlap overlap, ess sum(wts)^2 / sum(wts^2)) }estimand ATE和ATT的区别要拎清ATE 估的是“如果全人群都吸烟 vs 全人群都不吸烟”的平均差异ATT 只估“吸烟者如果没吸烟会怎样”。做暴露危害评估一般用 ATE。stabilize TRUE把权重除以处理组的边际概率能显著降低权重的方差。overlap里那两个n_lt01/n_gt99计数是诊断的核心——如果有一批人倾向评分大于 0.99说明模型认为他们“几乎肯定吸烟”而那批人恰好没吸烟权重会飙到 100 以上一个样本就能主导整个估计。trim 0.01是在 1% 和 99% 分位处截尾代价是引入轻微偏差但换来方差可控实践中这笔交易通常划算。ess有效样本量是另一个必看指标sum(w)^2 / sum(w^2)给出加权后相当于多少“等效样本”。如果原始 500 例加权后 ESS 只剩 60说明权重分布极其不均这个估计不该信。我的经验阈值是 ESS 低于原始的 30% 就要回头查协变量集。4.2 AIPW 双重稳健两个模型互相兜底IPW 只依赖倾向评分模型正确一旦模型设定错了比如该放交互项没放估计就有偏。AIPW 加了结局模型两个模型只要有一个对估计就一致这叫双重稳健。代价是计算量大一些但在这个场景里值得。结局模型可以用简单的线性回归也可以用grf的回归森林。项目里默认走线性因为基因表达和协变量的关系大多近似线性森林在小样本上反而方差大# AIPW 估计需要倾向评分 e(X) 和两个结局预测 mu1(X), mu0(X) estimate_aipw - function(y, w, x, use_forest FALSE) { df - data.frame(y y, W w, x) # 倾向评分模型 ps_fit - glm(W ~ ., data df, family binomial()) e - fitted(ps_fit) # 夹一下避免 e 正好等于 0 或 1 导致除零 e - pmin(pmax(e, 1e-6), 1 - 1e-6) # 结局模型分两组各拟合一个预测反事实结局 if (use_forest) { cf - grf::regression_forest(cbind(df$x, W df$W), df$y) mu1 - predict(cf, cbind(df$x, W 1))$predictions mu0 - predict(cf, cbind(df$x, W 0))$predictions } else { m1 - lm(y ~ ., data subset(df, W 1)) m0 - lm(y ~ ., data subset(df, W 0)) mu1 - predict(m1, newdata df) mu0 - predict(m0, newdata df) } # AIPW 影响函数逐样本取值 psi - (mu1 - mu0) w * (y - mu1) / e - (1 - w) * (y - mu0) / (1 - e) est - mean(psi) se - sd(psi) / sqrt(length(psi)) list(est est, se se, mu1 mean(mu1), mu0 mean(mu0), # 结局模型在两组上的外推幅度过说明模型在硬撑 extrap c(max_mu1 max(mu1), min_mu1 min(mu1), max_mu0 max(mu0), min_mu0 min(mu0))) }psi那四行是整套方法的核心。前两项mu1 - mu0是结局模型的预测效应后面两个减项是“模型残差乘以逆倾向评分”的修正项。直觉上如果结局模型预测得准修正项接近零估计主要由mu1 - mu0决定如果倾向评分模型更准修正项会把结局模型的偏差拉回来。e夹到[1e-6, 1-1e-6]是必须的现实数据里fitted(ps_fit)出现 1.0000 是常见的除法直接炸出 Inf。extrap那几个极值是我加的诊断看mu1的最大值有没有远超实际观测到的表达量范围——如果超过了说明结局模型在做荒谬的外推这条估计不能采信先回头检查协变量。use_forest TRUE那条分支用grf::regression_forest把W当成一个普通特征塞进X预测时分别把W置 1 和 0。这种写法比训练两个独立森林更高效但对样本量要求更高小于 200 例时慎用。4.3 因果森林平均效应之外的异质性IPW 和 AIPW 给的是全人群平均效应但吸烟对基因表达的影响很可能在亚群里不一致——比如年龄大的人反应更强或者只在不吸烟的遗传背景子群里显现。因果森林直接估每个样本的条件平均处理效应tau(X)然后看tau有没有系统性的分化。# 因果森林估条件平均处理效应看异质性 fit_causal_forest - function(y, w, x, num_trees 2000) { X - cbind(as.data.frame(x)) cf - grf::causal_forest( X X, Y y, W w, num.trees num_trees, honesty TRUE, # 诚实划分避免用同一批数据既分裂又估值 min.node.size 10, # 叶子最小样本太小会过拟合 sample.fraction 0.5, # 子采样比例 seed 20240101 ) # 全人群平均效应target.sample overlap 只在有重叠的人群上平均 ate - grf::average_treatment_effect(cf, target.sample overlap) # 变量重要性看哪些协变量在驱动异质性 vimp - grf::variable_importance(cf) list(cf cf, ate ate, vimp vimp) }honesty TRUE是grf的默认值作用是样本分裂和效应估计用不同的子样本防止过拟合把噪声当成异质性。这个参数不要关掉关了之后tau的分布会变得极其“好看”但全是假的。min.node.size 10是叶子节点的最小样本量基因表达数据的噪声大设太小会切出一堆只有几个样本的叶子估出来的tau方差爆炸。target.sample overlap这一项值得单独说。默认是all在全员上平均改成overlap只在倾向评分有重叠的区域上平均更稳健但也更保守。做暴露效应的基因筛查时用overlap结论不容易被少数极端个体带跑。variable_importance的输出是一个矩阵第一列是原始重要性第二列是归一化到 0-1 的版本。看前五名就够——如果年龄或某个技术协变量排第一说明异质性主要来自技术因素而非生物学这时候要警惕别急着往机制上解释。在基因层面因果森林只在候选基因子集上跑不是全基因组。原因很实际2000 棵树乘以几百个基因单机跑一晚上都不一定完而且大多数基因的异质性检验并不显著跑全量是浪费算力。4.4 从 p 值到 q 值结果表的组装三种估计跑完每个基因会拿到三套est和se。组装结果表要解决两个问题一是统一效应方向基因表达上调还是下调二是多重检验。基因数量上千直接拿 p 0.05 筛出来的基本都是假阳性。# 合并三种估计量的结果做 BH 校正输出可读的结果表 assemble_results - function(ipw_list, aipw_list, cf_list NULL) { build - function(lst, tag) { data.frame( gene_id names(lst), est vapply(lst, [[, numeric(1), est), se vapply(lst, [[, numeric(1), se), stringsAsFactors FALSE ) %% setNames(c(gene_id, paste0(tag, c(_est, _se)))) } res - Reduce(function(a, b) merge(a, b, by gene_id), list(build(ipw_list, ipw), build(aipw_list, aipw))) # 用 AIPW 作为主估计它双重稳健 res$z - res$aipw_est / res$aipw_se res$p - 2 * pnorm(-abs(res$z)) res$q - p.adjust(res$p, method BH) # 方向标记方便画火山图 res$dir - ifelse(res$q 0.05 res$aipw_est 0, up, ifelse(res$q 0.05 res$aipw_est 0, down, ns)) # 一致性旗标两种方法方向相反时标出来这类基因不写进结论 res$flag_inconsistent - sign(res$ipw_est) ! sign(res$aipw_est) res[order(res$q), ] }p.adjust(method BH)是 Benjamini-Hochberg控制 FDR比 Bonferroni 温和得多在基因表达场景是标准做法。dir那一列把显著基因分成上调和下调画图直接映射颜色。flag_inconsistent这一列是经验产物。IPW 和 AIPW 的方向如果相反说明两个模型的假设冲突得很厉害可能是倾向评分模型设定有问题也可能是结局模型的线性假设在这个基因上不成立。这类基因即使 q 值很小也不该写进结论最多放进补充材料标注方法间不一致等着用敏感性分析去解释。主估计选 AIPW 而不是 IPW理由是双重稳健只要倾向评分和结局模型有一个对估计就一致比单靠一个模型的 IPW 更扛得住设定误差。5. 避坑清单因果推断流水线里最容易翻车的五处这套流程跑通不难难的是跑出来的数能不能信。下面五条是我在这类项目上反复踩到的每一条都按“现象 → 原因 → 解决”整理配上排查代码。5.1 权重爆炸与重叠性缺失现象IPW 估出来的效应比 AIPW 大好几倍标准误也大得离谱ESS只剩原始样本的 10% 不到。原因倾向评分模型太“自信”给某些样本赋了接近 0 或 1 的概率。常见于某个协变量几乎能完全预测处理状态——比如某一批次的所有样本都来自吸烟组模型在这个批次上学到“批次 吸烟”权重自然爆掉。另一个成因是分类变量的稀疏水平某一类只有几个人逻辑回归在该维度上给出极端的系数。解决先看重叠性诊断再决定是修剪样本还是截尾权重。# 权重诊断分布、ESS、以及极端权重来自哪些样本 diagnose_weights - function(wts, meta) { ess - sum(wts)^2 / sum(wts^2) cat(sprintf(有效样本量 ESS %.1f / %d (%.1f%%)\n, ess, length(wts), 100 * ess / length(wts))) # 找出权重最大的 10 个样本看看它们有什么共同点 idx - order(wts, decreasing TRUE)[1:min(10, length(wts))] print(cbind(meta[idx, , drop FALSE], weight wts[idx])) invisible(ess) }打印极端权重样本的协变量是定位问题最快的路径。如果这十个样本全在同一个批次那就是批次导致的如果全是高龄那是年龄分布尾部的问题。ESS低于 30% 时我一般先尝试trim 0.01还不行就改用重叠权重WeightIt里method glm, estimand ATO重叠权重对极端倾向评分不敏感代价是估计目标变成“重叠人群的效应”解释时要写清楚。5.2 批次、共线性与“静默”的样本错位现象脚本一声不吭跑完结果表也有但效应估计整体偏小或者某个已知的强效应基因死活不显著。原因三种可能。批次与暴露强相关去批次时把处理效应一起减掉了协变量共线导致glm系数不稳定倾向评分本身就不准表达矩阵和表型的样本顺序没对齐merge或cbind之后列和行的对应关系已经乱了但 R 不会报错。解决三个检查必须都在预处理阶段做完——批次卡方检验3.2 节、设计矩阵秩检查3.3 节、样本 ID 断言。# 样本对齐的最后一道保险放在所有合并操作之后 assert_alignment - function(expr, meta) { stopifnot( 表达矩阵列名缺失 !is.null(colnames(expr)), 表型样本 ID 缺失 sample_id %in% names(meta), 样本数不一致 ncol(expr) nrow(meta), 样本顺序错位 identical(colnames(expr), as.character(meta$sample_id)) ) invisible(TRUE) }identical这一条是硬断言不做容错顺序不一致直接抛错。我宁可让脚本在第一步就炸也不想在写了两周报告之后发现样本对错了。5.3 因果森林在小样本上的“看起来很准”现象因果森林跑出来的tau分布特别漂亮亚群间差异明显但换一个随机种子结果就变了。原因honesty FALSE或者min.node.size设得太小树在噪声上过拟合。基因表达数据本身噪声大样本量几百例时2000 棵树的平均还压不住方差tau的分布就是随机的。解决三个参数一起调——保持honesty TRUEmin.node.size不低于 10sample.fraction取 0.4 到 0.6 之间然后用两个不同种子各跑一遍看tau的相关系数。# 稳定性检验换种子重跑看条件效应是否可复现 check_cf_stability - function(y, w, x, seeds c(1, 2024)) { taus - lapply(seeds, function(s) { cf - grf::causal_forest(X as.data.frame(x), Y y, W w, num.trees 2000, honesty TRUE, min.node.size 10, seed s) predict(cf)$predictions }) r - cor(taus[[1]], taus[[2]]) cat(sprintf(跨种子 tau 相关系数 %.3f\n, r)) if (r 0.7) warning(条件效应不稳定异质性结论不可靠) invisible(r) }相关系数低于 0.7 就别写异质性结论了老老实实报平均效应。这个检验成本不高值得每个候选基因都过一遍。5.4 多重检验的两种误用现象q 值算出来了显著基因上千个但和已知的吸烟相关基因对不上。原因要么用了 Bonferroni在基因表达场景过于严苛真正的信号被压掉要么在筛选和验证之间反复调阈值最后报的是“调出来的”显著性。还有一种隐蔽的IPW 和 AIPW 各算一遍 q 值挑显著的那个报这是变相的多重比较。解决主分析只跑一套AIPWq 值只算一次阈值预先定死比如 q 0.05敏感性分析的结果单独列不参与筛选。方法学部分要写清楚“主估计用的是哪一个”。5.5 用同一份数据调参又评估现象所有诊断都完美重叠性好、ESS 高、q 值分布漂亮但换一个独立队列就复现不出来。原因倾向评分模型的协变量集、trim阈值、因果森林的参数如果在看到结果之后反复调整等于把结局信息泄漏进了模型选择。这类偏差不会体现在任何一张诊断图上。解决把样本按 7:3 随机切开前 70% 用来定协变量集和超参数后 30% 只跑一次不看结果不回头。如果样本量不够切就在分析计划里把参数全部预先写死中途改了要记录在日志里。# 预注册式的参数冻结把最终参数写进文件后续脚本只读不猜 freeze_params - function(cfg, path output/tables/frozen_params.rds) { saveRDS(cfg, path) message(参数已冻结至 , path, 后续分析请勿修改) }这条看着像流程洁癖但它是唯一能防住“调参调到显著”的办法。玄学的地方在于你自己调的时候完全不觉得在作弊每一步都有“合理”的理由。6. 二次开发把单基因流程改成批量筛选与可复现报告拿到这套脚本的人十有八九第一件事是想跑全基因组。单基因函数写好了批量跑不是简单套个for循环——直接循环会撞上三个问题内存、算力和失败恢复。我给的做法是把估计函数包一层加上并行和缓存。先看批量层怎么写# 批量估计并行 单基因失败不影响整体 结果落盘缓存 run_batch - function(expr_mat, w, x, method c(ipw, aipw), n_workers 4) { method - match.arg(method) future::plan(future::multisession, workers n_workers) on.exit(future::plan(future::sequential), add TRUE) gene_ids - rownames(expr_mat) est_fun - if (method ipw) estimate_ipw else estimate_aipw res - furrr::future_map(gene_ids, function(g) { y - as.numeric(expr_mat[g, ]) # 零方差基因直接跳过否则 glm 会抛错并中断整个批次 if (stats::var(y, na.rm TRUE) 1e-8) return(NULL) tryCatch( est_fun(y y, w w, x x), error function(e) { structure(list(err conditionMessage(e)), class fail) } ) }, .options furrr::furrr_options(seed TRUE)) names(res) - gene_ids fails - vapply(res, inherits, logical(1), fail) if (any(fails)) message(sprintf(%d 个基因估计失败已跳过, sum(fails))) res[!fails] }furrr_options(seed TRUE)保证每个 worker 的随机数流可复现不加这一句因果森林那部分每次跑出来都不一样。tryCatch把单基因的失败隔离掉glm在完全分离的数据上会抛 warning 甚至 error一个基因炸掉不该拖垮整批。零方差判断放在最前面是最常见的崩溃点——某些基因在所有样本里表达值都一样lm直接返回 NA 系数。n_workers不建议超过物理核心数减一grf内部本身也会占内存开太多反而会因为内存交换变慢。第二个改造点是报告层。原始脚本输出的是 CSV二次开发可以把它变成参数化的 Rmd 或 Quarto 文档把 q 值分布、火山图、前若干个基因的森林图一次性生成换一批数据只改cfg里的路径。这里的关键是把“结果表”和“图”解耦图只读表不重算估计。第三个改造点是加一层敏感性分析。做法很简单把trim从 0.01 换成 0.05 和 0.1各跑一遍 AIPW看显著基因的重叠比例。敏感性维度参数观察指标警戒线权重截尾trim∈ {0.01, 0.05, 0.1}显著基因重叠率低于 70% 说明结论依赖阈值结局模型线性 vs 回归森林效应估计相关系数低于 0.8 说明模型设定敏感协变量集加/减批次变量q 0.05 基因数变化数量差一倍以上要解释这张表建议直接贴进补充材料比任何文字描述都有说服力。审稿人问“你的结论稳健吗”把重叠率甩过去就行。最后一个技巧是日志。把sessionInfo()、cfg、每个脚本的运行时间和输入文件的行列数统一写进output/logs/一个脚本一个文件。这个习惯我是被坑出来的有次报告里的数和三个月前跑的对不上翻遍代码也没改过最后发现是上游给的表达矩阵悄悄更新了一版样本数从 512 变成了 508。从那以后我每次跑因果估计流水线都强制在第一步把sessionInfo()和输入文件的 md5 一起写进日志跑完再核对一遍才敢往报告里放数。这类工程的可靠性不在估计量写得多漂亮而在每一步都能回溯。希望帮到你。本文还有配套的精品资源点击获取
简介SmokingGeneExpression 是一份面向生物信息学与统计分析学习者的 R 语言项目资料围绕吸烟对健康肺组织基因表达的短期因果效应展开。内容将因果推断思路落到基因表达场景涉及工具变量、倾向得分匹配、有向无环图等混杂控制方法以及 RNA-seq 差异表达分析、PCA 可视化与生物学解释适合希望理解观察性数据因果分析流程并借助 R 复现实战的读者。压缩包共 45 个文件约 3.1MB以 png 图表、txt 结果说明、xlsx 与 csv 数据表为主另含 R 脚本和 pdf 文档覆盖分析代码、数据样例与结果呈现等用途。已有 148 人学习可作为课程作业、课题预研或方法对照的参考。读者可从中获得较完整的分析素材按目录梳理研究逻辑、复现关键步骤并积累 R 语言实践。1. 项目定位用 R 把吸烟对肺组织基因表达的因果效应算清楚一个做肺组织转录组的同行跟我吐槽过差异表达分析跑完吸烟组和从不吸烟组之间挂着三千多个显著基因图表也漂亮结果报告交上去第一个问题就是“你怎么知道这些差异不是吸烟者普遍年龄更大、BMI 更高造成的”。相关性分析能告诉你两组人不一样但不能告诉你是吸烟把某个基因的表达抬上去的。SmokingGeneExpression 这个 R 项目要解决的正是这一层把“吸烟”当作处理变量把“某个基因的表达量”当作结局用因果推断的估计量把混杂因素压下去给出每个基因在健康肺组织里对吸烟暴露的短期因果效应估计。适合谁来读手上有队列或公开表达资源、会写基本 R 脚本、但还没系统跑过倾向评分加权和双重稳健估计的人。项目的技术栈是纯 R——数据层用 data.table 和基础的 Bioconductor 表达处理估计层主要靠 WeightIt、grf、sandwich 这几个包可视化走 ggplot2。它不训练模型、不做深度学习整套流程跑完就是一张带置信区间和 q 值的基因效应表外加几张诊断图。对做生信但被审稿意见追着问因果的人这个骨架的参考价值比多跑十个差异分析大。2. 工程结构R 脚本怎么分工、依赖怎么装、第一次怎么跑通这个项目不是包package是一套按分析阶段切开的脚本工程。这意味着它没有R CMD check那套约束好处是改起来自由坏处是环境不一致时排错特别费劲——这是后面避坑章节的重要伏笔。2.1 依赖清单哪些包是必需的哪些是可选的依赖分成三档。第一档是数据处理data.table、dplyr、tidyr、readr第二档是因果估计核心WeightIt倾向评分与权重、grf因果森林与 AIPW 的机器学习 nuisance 估计、sandwich加lmtest聚类稳健标准误第三档是多重检验与出图qvalue、ggplot2、patchwork。并行那块用futurefurrr规模化跑全基因组时会开关它。环境锁定用renv项目根目录有一份renv.lock。恢复环境的标准动作# 在项目根目录执行renv 会按 lock 文件把包版本拉回一致状态 Rscript -e if (!requireNamespace(renv, quietly TRUE)) install.packages(renv); renv::restore()这条命令的重点不是“安装依赖”而是“安装指定版本的依赖”。grf和WeightIt这两个包在版本升级时函数签名改过尤其是WeightIt的estimand参数和grf的honesty默认值跨版本跑出来的数会对不上。我第一次拿到这套脚本时偷懒没跑 restore直接用系统里现装的包结果 IPW 和 AIPW 的效应估计差了将近 8%排查了两天才发现是权重归一化方式变了。血泪经验这类因果估计工程包版本必须锁死。如果机器上装不了qvalueBioconductor 包退路是用p.adjust(method BH)两者在大多数情况下结论一致只是 q 值的具体数值略有差异报告里注明用的是哪一种即可。2.2 目录分层data、R、scripts、output 各管什么工程目录大致是这样切的路径内容说明data/raw/原始表达矩阵、表型表、协变量表只读任何脚本不得写回data/derived/标准化后的表达、构建好的分析集中间产物可全部删除重跑R/函数层io.R、preprocess.R、estimators.R、sensitivity.R、utils.R只定义函数不执行分析scripts/执行层01_prepare_data.R到05_report.R有编号按顺序跑output/tables/基因效应表、诊断统计表交付物output/figures/森林图、权重分布图、QQ 图交付物这个分层是有讲究的。R/里的函数不碰全局变量参数从外面传进去方便单独做单元测试scripts/里的脚本只负责“读什么、调哪个函数、写到哪里”。很多生信分析脚本写着写着就变成一个两千行的巨型文件中间变量到处飞改一个参数要重跑全流程这种结构就是为了防这个。常见的做法是每个scripts/开头统一 source 一遍函数层# scripts/01_prepare_data.R 的开头固定写这一块 suppressPackageStartupMessages({ library(data.table) library(dplyr) }) # 统一从项目根目录解析路径避免依赖 setwd proj_root - here::here() source(file.path(proj_root, R, io.R)) source(file.path(proj_root, R, preprocess.R)) source(file.path(proj_root, R, estimators.R)) # 全局参数集中在这一处后面所有脚本读同一个配置 cfg - list( expr_file file.path(proj_root, data/raw/expression_matrix.tsv), pheno_file file.path(proj_root, data/raw/phenotype.tsv), covar_file file.path(proj_root, data/raw/covariates.tsv), min_expr 1, # 表达过滤阈值log2 尺度下取 1 大致对应低表达噪声层 min_samples 0.3, # 一个基因至少在 30% 样本中超过阈值才保留 seed 20240101 )cfg这个列表是整条流水线的黑匣子开关。min_expr和min_samples直接决定后面有多少基因进入多重检验改动了结论就得重出seed因为因果森林内部有随机抽样不固定的话两次跑出来的效应估计在小数点后几位会飘。路径用here::here()而不是相对路径是为了让脚本从 RStudio、命令行、还是调度系统提交都能找到文件这个小习惯能省掉大量“在我机器上能跑”的扯皮。2.3 第一次跑通从 01 到 05 的执行顺序五个脚本是一条串行链但中间结果都落盘所以可以从任意一步重跑# 全流程按编号顺序执行 Rscript scripts/01_prepare_data.R Rscript scripts/02_fit_ipw.R Rscript scripts/03_fit_aipw.R Rscript scripts/04_causal_forest.R Rscript scripts/05_report.R01做数据对齐和表达预处理输出data/derived/analysis_set.rds02和03分别在同一个分析集上跑 IPW 和 AIPW各自产出一张基因级效应表04只在候选基因子集上跑因果森林用来观察异质性05把三张表按基因 ID 合并、统一多重检验校正、出图和写报告。第一次跑建议先在01里把基因数砍到几百个做冒烟测试确认整条链没断再放开全量。基因数上万的表达矩阵走完全套单核大概要几十分钟到几小时取决于样本量和是否开并行。冒烟测试这一步别省04的因果森林在小样本上最容易因为某个基因方差为零直接抛错早发现早改。3. 数据层拆解暴露、结局与协变量怎么对齐到同一张表因果推断里最容易翻车的从来不是估计量写错而是数据对齐错了。这个项目的核心数据结构只有一张宽表每一行是一个样本列包括样本 ID、吸烟暴露、几十上百个协变量、以及每个基因的表达值。所有估计函数都只认这张表保证口径一致。3.1 吸烟暴露变量的编码二值化之前的三个判断原始表型里吸烟状态通常是分好几档的从不、曾经、当前、偶尔、已戒若干年甚至还有每天支数。直接二值化会丢信息但项目走的是“当前吸烟 vs 从不吸烟”这条二值路线理由有三短期因果效应下“曾经吸烟”人群的暴露已经停止把它塞进处理组会稀释效应每天支数在多数队列里缺失率太高做连续暴露需要另一套剂量-反应模型而二值暴露能直接套用倾向评分框架诊断工具也成熟。编码这一步要把“曾经吸烟”和缺失显式剔出去而不是当成对照组# 把原始吸烟编码映射成 0/1 的处理变量并显式记录剔除原因 build_exposure - function(pheno) { stopifnot(all(c(sample_id, smoke_status) %in% names(pheno))) out - pheno %% mutate( smoke_status tolower(trimws(smoke_status)), # 只保留两种明确状态其余一律置 NA W case_when( smoke_status %in% c(current, current smoker, daily) ~ 1L, smoke_status %in% c(never, never smoker) ~ 0L, TRUE ~ NA_integer_ ), # 单独打个标记方便统计各剔除原因的数量 drop_reason case_when( is.na(W) smoke_status %in% c(former, past) ~ former_smoker, is.na(W) is.na(smoke_status) ~ missing_status, is.na(W) smoke_status ~ blank_status, is.na(W) ~ other_status, TRUE ~ NA_character_ ) ) attr(out, exposure_counts) - table(out$W, useNA always) attr(out, drop_counts) - table(out$drop_reason, useNA no) out }case_when的匹配顺序很重要%in%里放了多种写法是为了兼容不同来源的表型表真实数据里 “Current Smoker” 和 “current” 混着出现是常态。drop_reason这一列看起来多余但它是后面写方法学部分的证据——审稿人一定会问你为什么最终样本是 N你得能答出“剔除了多少既往吸烟者、多少状态缺失”。attr挂上去的计数会被01脚本打印到日志里养成这个习惯能省掉回头重算的麻烦。还有一层容易忽略的暴露的时间窗。项目说的是“短期”那就得确认表型采集时间和组织取样时间接近如果表型是问卷回顾性的“过去一年是否吸烟”而样本是五年前存的这个短期假设就不成立。数据里如果有采样日期字段务必在预处理阶段做个时间差统计。3.2 表达矩阵的预处理过滤、标准化、以及要不要去批次表达矩阵进模型前要过三道关。第一道是低表达过滤用min_expr和min_samples两个参数控制目的是减少检验数量、降低多重校正的惩罚。第二道是分布变换如果是原始 count 就走log2(CPM 1)或vst如果上游已经给了 log2 尺度就跳过。第三道是批次这一步是分水岭批次如果和吸烟状态强相关去批次会把效应一起减掉如果批次独立于暴露去不去影响都不大。保守做法是先检查批次与暴露的关联再做决定# 过滤 标准化并在返回前检查批次是否与暴露混杂 prep_expression - function(expr, meta, cfg) { mat - as.matrix(expr[, -1, drop FALSE]) rownames(mat) - expr[[1]] # 1) 低表达过滤两个条件同时满足才保留 keep - rowMeans(mat cfg$min_expr) cfg$min_samples mat - mat[keep, , drop FALSE] # 2) log2 变换已经在 log 尺度上的矩阵跳过 if (max(mat, na.rm TRUE) 50) mat - log2(mat 1) # 3) 分位数标准化抹平样本间整体分布差异 mat - limma::normalizeBetweenArrays(mat, method quantile) # 4) 批次与暴露的关联检查卡方检验p 小说明批次和吸烟纠缠在一起 if (!is.null(meta$batch) !all(is.na(meta$batch))) { tab - table(meta$batch, meta$W) p_batch - suppressWarnings(chisq.test(tab)$p.value) message(sprintf(批次-暴露关联检验 p %.4g, p_batch)) if (!is.na(p_batch) p_batch 0.05) { warning(批次与吸烟状态显著关联去批次可能吸收掉部分处理效应慎用 ComBat) } } list(expr mat, n_gene nrow(mat), n_sample ncol(mat)) }rowMeans(mat cfg$min_expr)这一行是过滤的核心表达式返回每个基因超过阈值的样本比例再和min_samples比。normalizeBetweenArrays走 quantile 方法对偏态分布的表达值比中位数标准化更稳。最值得看的是第 4 段——批次卡方检验。如果 p 值很小说明这批数据里“哪台机器测的”和“这个人吸不吸烟”高度相关此时上 ComBat 或limma::removeBatchEffect等于把处理效应的一部分当成批次削掉了后面估计出来的效应会系统性偏小。这种情况我一般不去批次而是把批次当作协变量放进倾向评分模型让加权去平衡它而不是从数据里删掉它。矩阵的样本顺序和表型的样本顺序必须严格一致这是最容易出错的地方。prep_expression返回前应该强制按meta$sample_id重排或者用stopifnot(identical(colnames(mat), meta$sample_id))做一次断言别指望两边碰巧一致。3.3 协变量集与设计矩阵共线性检查协变量集的选择原则是只放暴露之前就已经确定的变量年龄、性别、遗传背景主成分、测序深度、批次不要放可能被暴露影响的中介变量比如某些炎症指标、肺功能。放中介进去会把总效应部分遮掉估出来的是直接效应解释起来是另一回事。确定变量集之后必须检查设计矩阵的秩。协变量之间完全共线比如两个高度相关的技术指标或者某个分类变量某一水平只有一个人会让倾向评分模型直接崩掉glm可能不报错但给出荒谬的系数。# 构建协变量设计矩阵并做秩与稀疏水平检查 check_design - function(meta, covars) { f - reformulate(covars) mm - model.matrix(f, data meta) # 秩亏说明存在完全共线glm 的结果不可信 r - qr(mm)$rank message(sprintf(设计矩阵维度 %d x %d秩 %d, nrow(mm), ncol(mm), r)) if (r ncol(mm)) warning(设计矩阵秩亏存在完全共线协变量需手动剔除) # 分类变量稀疏水平某个水平样本过少会让该水平的倾向评分极端 for (v in covars) { if (is.character(meta[[v]]) || is.factor(meta[[v]])) { tb - table(meta[[v]]) if (any(tb 20)) warning(sprintf(变量 %s 存在样本量 20 的水平%s, v, paste(names(tb)[tb 20], collapse , ))) } } invisible(mm) }qr(mm)$rank小于列数就是秩亏最常见的成因是分类变量做了独热编码后某一列全是零或者两个连续变量做了线性组合。tb 20这个阈值不是硬标准20 只是个提醒线真正要看的是该水平在两组里的分布——如果某个批次只有 3 个人且全在吸烟组倾向评分在这个批次上会趋近 1加权时必然出极端权重这就是第 5 章要展开的问题。提前在这里发现比跑到加权阶段再看权重分布图省事得多。4. 因果估计核心三种估计量在单个基因上的完整链路数据对齐之后估计层的输入就是一个向量结局Y某个基因的表达值、处理W0/1、协变量矩阵X。项目给三种估计量不是让你三选一而是让它们互相印证IPW 看加权后的可比人群AIPW 用两个模型互相兜底因果森林看效应是不是在亚群里不一致。三者结论方向一致才敢往报告里写。4.1 倾向评分与 IPW把不可比的人群“造”成可比倾向评分e(X) P(W1|X)的直觉是给每个样本算一个“被吸烟的概率”然后用这个概率的倒数做权重让协变量分布在两组之间拉平。吸烟概率高的人如果实际没吸烟权重就大代表他“替”了很多吸烟者。实现上用WeightIt它比手写glm多了一层诊断# 单个基因的 IPW 估计输入是结局向量 处理向量 协变量数据框 estimate_ipw - function(y, w, x, trim 0.01) { df - data.frame(y y, W w, x) # 用 logistic 回归估倾向评分estimand ATE 表示估的是全人群平均效应 wt - WeightIt::weightit(W ~ ., data df, method glm, estimand ATE, stabilize TRUE) ps - wt$ps # 重叠性诊断倾向评分贴边说明两组人根本不可比 overlap - c(min min(ps), max max(ps), n_lt01 sum(ps 0.01), n_gt99 sum(ps 0.99)) # 权重截尾压掉极端权重的杠杆 wts - wt$weights if (!is.null(trim) trim 0) { q - quantile(wts, probs c(trim, 1 - trim)) wts - pmin(pmax(wts, q[1]), q[2]) } # 加权最小二乘给 ATErobust TRUE 拿稳健标准误 fit - lm(y ~ W, data df, weights wts) se - sqrt(sandwich::vcovHC(fit, type HC1)[W, W]) list(est coef(fit)[[W]], se se, ps ps, overlap overlap, ess sum(wts)^2 / sum(wts^2)) }estimand ATE和ATT的区别要拎清ATE 估的是“如果全人群都吸烟 vs 全人群都不吸烟”的平均差异ATT 只估“吸烟者如果没吸烟会怎样”。做暴露危害评估一般用 ATE。stabilize TRUE把权重除以处理组的边际概率能显著降低权重的方差。overlap里那两个n_lt01/n_gt99计数是诊断的核心——如果有一批人倾向评分大于 0.99说明模型认为他们“几乎肯定吸烟”而那批人恰好没吸烟权重会飙到 100 以上一个样本就能主导整个估计。trim 0.01是在 1% 和 99% 分位处截尾代价是引入轻微偏差但换来方差可控实践中这笔交易通常划算。ess有效样本量是另一个必看指标sum(w)^2 / sum(w^2)给出加权后相当于多少“等效样本”。如果原始 500 例加权后 ESS 只剩 60说明权重分布极其不均这个估计不该信。我的经验阈值是 ESS 低于原始的 30% 就要回头查协变量集。4.2 AIPW 双重稳健两个模型互相兜底IPW 只依赖倾向评分模型正确一旦模型设定错了比如该放交互项没放估计就有偏。AIPW 加了结局模型两个模型只要有一个对估计就一致这叫双重稳健。代价是计算量大一些但在这个场景里值得。结局模型可以用简单的线性回归也可以用grf的回归森林。项目里默认走线性因为基因表达和协变量的关系大多近似线性森林在小样本上反而方差大# AIPW 估计需要倾向评分 e(X) 和两个结局预测 mu1(X), mu0(X) estimate_aipw - function(y, w, x, use_forest FALSE) { df - data.frame(y y, W w, x) # 倾向评分模型 ps_fit - glm(W ~ ., data df, family binomial()) e - fitted(ps_fit) # 夹一下避免 e 正好等于 0 或 1 导致除零 e - pmin(pmax(e, 1e-6), 1 - 1e-6) # 结局模型分两组各拟合一个预测反事实结局 if (use_forest) { cf - grf::regression_forest(cbind(df$x, W df$W), df$y) mu1 - predict(cf, cbind(df$x, W 1))$predictions mu0 - predict(cf, cbind(df$x, W 0))$predictions } else { m1 - lm(y ~ ., data subset(df, W 1)) m0 - lm(y ~ ., data subset(df, W 0)) mu1 - predict(m1, newdata df) mu0 - predict(m0, newdata df) } # AIPW 影响函数逐样本取值 psi - (mu1 - mu0) w * (y - mu1) / e - (1 - w) * (y - mu0) / (1 - e) est - mean(psi) se - sd(psi) / sqrt(length(psi)) list(est est, se se, mu1 mean(mu1), mu0 mean(mu0), # 结局模型在两组上的外推幅度过说明模型在硬撑 extrap c(max_mu1 max(mu1), min_mu1 min(mu1), max_mu0 max(mu0), min_mu0 min(mu0))) }psi那四行是整套方法的核心。前两项mu1 - mu0是结局模型的预测效应后面两个减项是“模型残差乘以逆倾向评分”的修正项。直觉上如果结局模型预测得准修正项接近零估计主要由mu1 - mu0决定如果倾向评分模型更准修正项会把结局模型的偏差拉回来。e夹到[1e-6, 1-1e-6]是必须的现实数据里fitted(ps_fit)出现 1.0000 是常见的除法直接炸出 Inf。extrap那几个极值是我加的诊断看mu1的最大值有没有远超实际观测到的表达量范围——如果超过了说明结局模型在做荒谬的外推这条估计不能采信先回头检查协变量。use_forest TRUE那条分支用grf::regression_forest把W当成一个普通特征塞进X预测时分别把W置 1 和 0。这种写法比训练两个独立森林更高效但对样本量要求更高小于 200 例时慎用。4.3 因果森林平均效应之外的异质性IPW 和 AIPW 给的是全人群平均效应但吸烟对基因表达的影响很可能在亚群里不一致——比如年龄大的人反应更强或者只在不吸烟的遗传背景子群里显现。因果森林直接估每个样本的条件平均处理效应tau(X)然后看tau有没有系统性的分化。# 因果森林估条件平均处理效应看异质性 fit_causal_forest - function(y, w, x, num_trees 2000) { X - cbind(as.data.frame(x)) cf - grf::causal_forest( X X, Y y, W w, num.trees num_trees, honesty TRUE, # 诚实划分避免用同一批数据既分裂又估值 min.node.size 10, # 叶子最小样本太小会过拟合 sample.fraction 0.5, # 子采样比例 seed 20240101 ) # 全人群平均效应target.sample overlap 只在有重叠的人群上平均 ate - grf::average_treatment_effect(cf, target.sample overlap) # 变量重要性看哪些协变量在驱动异质性 vimp - grf::variable_importance(cf) list(cf cf, ate ate, vimp vimp) }honesty TRUE是grf的默认值作用是样本分裂和效应估计用不同的子样本防止过拟合把噪声当成异质性。这个参数不要关掉关了之后tau的分布会变得极其“好看”但全是假的。min.node.size 10是叶子节点的最小样本量基因表达数据的噪声大设太小会切出一堆只有几个样本的叶子估出来的tau方差爆炸。target.sample overlap这一项值得单独说。默认是all在全员上平均改成overlap只在倾向评分有重叠的区域上平均更稳健但也更保守。做暴露效应的基因筛查时用overlap结论不容易被少数极端个体带跑。variable_importance的输出是一个矩阵第一列是原始重要性第二列是归一化到 0-1 的版本。看前五名就够——如果年龄或某个技术协变量排第一说明异质性主要来自技术因素而非生物学这时候要警惕别急着往机制上解释。在基因层面因果森林只在候选基因子集上跑不是全基因组。原因很实际2000 棵树乘以几百个基因单机跑一晚上都不一定完而且大多数基因的异质性检验并不显著跑全量是浪费算力。4.4 从 p 值到 q 值结果表的组装三种估计跑完每个基因会拿到三套est和se。组装结果表要解决两个问题一是统一效应方向基因表达上调还是下调二是多重检验。基因数量上千直接拿 p 0.05 筛出来的基本都是假阳性。# 合并三种估计量的结果做 BH 校正输出可读的结果表 assemble_results - function(ipw_list, aipw_list, cf_list NULL) { build - function(lst, tag) { data.frame( gene_id names(lst), est vapply(lst, [[, numeric(1), est), se vapply(lst, [[, numeric(1), se), stringsAsFactors FALSE ) %% setNames(c(gene_id, paste0(tag, c(_est, _se)))) } res - Reduce(function(a, b) merge(a, b, by gene_id), list(build(ipw_list, ipw), build(aipw_list, aipw))) # 用 AIPW 作为主估计它双重稳健 res$z - res$aipw_est / res$aipw_se res$p - 2 * pnorm(-abs(res$z)) res$q - p.adjust(res$p, method BH) # 方向标记方便画火山图 res$dir - ifelse(res$q 0.05 res$aipw_est 0, up, ifelse(res$q 0.05 res$aipw_est 0, down, ns)) # 一致性旗标两种方法方向相反时标出来这类基因不写进结论 res$flag_inconsistent - sign(res$ipw_est) ! sign(res$aipw_est) res[order(res$q), ] }p.adjust(method BH)是 Benjamini-Hochberg控制 FDR比 Bonferroni 温和得多在基因表达场景是标准做法。dir那一列把显著基因分成上调和下调画图直接映射颜色。flag_inconsistent这一列是经验产物。IPW 和 AIPW 的方向如果相反说明两个模型的假设冲突得很厉害可能是倾向评分模型设定有问题也可能是结局模型的线性假设在这个基因上不成立。这类基因即使 q 值很小也不该写进结论最多放进补充材料标注方法间不一致等着用敏感性分析去解释。主估计选 AIPW 而不是 IPW理由是双重稳健只要倾向评分和结局模型有一个对估计就一致比单靠一个模型的 IPW 更扛得住设定误差。5. 避坑清单因果推断流水线里最容易翻车的五处这套流程跑通不难难的是跑出来的数能不能信。下面五条是我在这类项目上反复踩到的每一条都按“现象 → 原因 → 解决”整理配上排查代码。5.1 权重爆炸与重叠性缺失现象IPW 估出来的效应比 AIPW 大好几倍标准误也大得离谱ESS只剩原始样本的 10% 不到。原因倾向评分模型太“自信”给某些样本赋了接近 0 或 1 的概率。常见于某个协变量几乎能完全预测处理状态——比如某一批次的所有样本都来自吸烟组模型在这个批次上学到“批次 吸烟”权重自然爆掉。另一个成因是分类变量的稀疏水平某一类只有几个人逻辑回归在该维度上给出极端的系数。解决先看重叠性诊断再决定是修剪样本还是截尾权重。# 权重诊断分布、ESS、以及极端权重来自哪些样本 diagnose_weights - function(wts, meta) { ess - sum(wts)^2 / sum(wts^2) cat(sprintf(有效样本量 ESS %.1f / %d (%.1f%%)\n, ess, length(wts), 100 * ess / length(wts))) # 找出权重最大的 10 个样本看看它们有什么共同点 idx - order(wts, decreasing TRUE)[1:min(10, length(wts))] print(cbind(meta[idx, , drop FALSE], weight wts[idx])) invisible(ess) }打印极端权重样本的协变量是定位问题最快的路径。如果这十个样本全在同一个批次那就是批次导致的如果全是高龄那是年龄分布尾部的问题。ESS低于 30% 时我一般先尝试trim 0.01还不行就改用重叠权重WeightIt里method glm, estimand ATO重叠权重对极端倾向评分不敏感代价是估计目标变成“重叠人群的效应”解释时要写清楚。5.2 批次、共线性与“静默”的样本错位现象脚本一声不吭跑完结果表也有但效应估计整体偏小或者某个已知的强效应基因死活不显著。原因三种可能。批次与暴露强相关去批次时把处理效应一起减掉了协变量共线导致glm系数不稳定倾向评分本身就不准表达矩阵和表型的样本顺序没对齐merge或cbind之后列和行的对应关系已经乱了但 R 不会报错。解决三个检查必须都在预处理阶段做完——批次卡方检验3.2 节、设计矩阵秩检查3.3 节、样本 ID 断言。# 样本对齐的最后一道保险放在所有合并操作之后 assert_alignment - function(expr, meta) { stopifnot( 表达矩阵列名缺失 !is.null(colnames(expr)), 表型样本 ID 缺失 sample_id %in% names(meta), 样本数不一致 ncol(expr) nrow(meta), 样本顺序错位 identical(colnames(expr), as.character(meta$sample_id)) ) invisible(TRUE) }identical这一条是硬断言不做容错顺序不一致直接抛错。我宁可让脚本在第一步就炸也不想在写了两周报告之后发现样本对错了。5.3 因果森林在小样本上的“看起来很准”现象因果森林跑出来的tau分布特别漂亮亚群间差异明显但换一个随机种子结果就变了。原因honesty FALSE或者min.node.size设得太小树在噪声上过拟合。基因表达数据本身噪声大样本量几百例时2000 棵树的平均还压不住方差tau的分布就是随机的。解决三个参数一起调——保持honesty TRUEmin.node.size不低于 10sample.fraction取 0.4 到 0.6 之间然后用两个不同种子各跑一遍看tau的相关系数。# 稳定性检验换种子重跑看条件效应是否可复现 check_cf_stability - function(y, w, x, seeds c(1, 2024)) { taus - lapply(seeds, function(s) { cf - grf::causal_forest(X as.data.frame(x), Y y, W w, num.trees 2000, honesty TRUE, min.node.size 10, seed s) predict(cf)$predictions }) r - cor(taus[[1]], taus[[2]]) cat(sprintf(跨种子 tau 相关系数 %.3f\n, r)) if (r 0.7) warning(条件效应不稳定异质性结论不可靠) invisible(r) }相关系数低于 0.7 就别写异质性结论了老老实实报平均效应。这个检验成本不高值得每个候选基因都过一遍。5.4 多重检验的两种误用现象q 值算出来了显著基因上千个但和已知的吸烟相关基因对不上。原因要么用了 Bonferroni在基因表达场景过于严苛真正的信号被压掉要么在筛选和验证之间反复调阈值最后报的是“调出来的”显著性。还有一种隐蔽的IPW 和 AIPW 各算一遍 q 值挑显著的那个报这是变相的多重比较。解决主分析只跑一套AIPWq 值只算一次阈值预先定死比如 q 0.05敏感性分析的结果单独列不参与筛选。方法学部分要写清楚“主估计用的是哪一个”。5.5 用同一份数据调参又评估现象所有诊断都完美重叠性好、ESS 高、q 值分布漂亮但换一个独立队列就复现不出来。原因倾向评分模型的协变量集、trim阈值、因果森林的参数如果在看到结果之后反复调整等于把结局信息泄漏进了模型选择。这类偏差不会体现在任何一张诊断图上。解决把样本按 7:3 随机切开前 70% 用来定协变量集和超参数后 30% 只跑一次不看结果不回头。如果样本量不够切就在分析计划里把参数全部预先写死中途改了要记录在日志里。# 预注册式的参数冻结把最终参数写进文件后续脚本只读不猜 freeze_params - function(cfg, path output/tables/frozen_params.rds) { saveRDS(cfg, path) message(参数已冻结至 , path, 后续分析请勿修改) }这条看着像流程洁癖但它是唯一能防住“调参调到显著”的办法。玄学的地方在于你自己调的时候完全不觉得在作弊每一步都有“合理”的理由。6. 二次开发把单基因流程改成批量筛选与可复现报告拿到这套脚本的人十有八九第一件事是想跑全基因组。单基因函数写好了批量跑不是简单套个for循环——直接循环会撞上三个问题内存、算力和失败恢复。我给的做法是把估计函数包一层加上并行和缓存。先看批量层怎么写# 批量估计并行 单基因失败不影响整体 结果落盘缓存 run_batch - function(expr_mat, w, x, method c(ipw, aipw), n_workers 4) { method - match.arg(method) future::plan(future::multisession, workers n_workers) on.exit(future::plan(future::sequential), add TRUE) gene_ids - rownames(expr_mat) est_fun - if (method ipw) estimate_ipw else estimate_aipw res - furrr::future_map(gene_ids, function(g) { y - as.numeric(expr_mat[g, ]) # 零方差基因直接跳过否则 glm 会抛错并中断整个批次 if (stats::var(y, na.rm TRUE) 1e-8) return(NULL) tryCatch( est_fun(y y, w w, x x), error function(e) { structure(list(err conditionMessage(e)), class fail) } ) }, .options furrr::furrr_options(seed TRUE)) names(res) - gene_ids fails - vapply(res, inherits, logical(1), fail) if (any(fails)) message(sprintf(%d 个基因估计失败已跳过, sum(fails))) res[!fails] }furrr_options(seed TRUE)保证每个 worker 的随机数流可复现不加这一句因果森林那部分每次跑出来都不一样。tryCatch把单基因的失败隔离掉glm在完全分离的数据上会抛 warning 甚至 error一个基因炸掉不该拖垮整批。零方差判断放在最前面是最常见的崩溃点——某些基因在所有样本里表达值都一样lm直接返回 NA 系数。n_workers不建议超过物理核心数减一grf内部本身也会占内存开太多反而会因为内存交换变慢。第二个改造点是报告层。原始脚本输出的是 CSV二次开发可以把它变成参数化的 Rmd 或 Quarto 文档把 q 值分布、火山图、前若干个基因的森林图一次性生成换一批数据只改cfg里的路径。这里的关键是把“结果表”和“图”解耦图只读表不重算估计。第三个改造点是加一层敏感性分析。做法很简单把trim从 0.01 换成 0.05 和 0.1各跑一遍 AIPW看显著基因的重叠比例。敏感性维度参数观察指标警戒线权重截尾trim∈ {0.01, 0.05, 0.1}显著基因重叠率低于 70% 说明结论依赖阈值结局模型线性 vs 回归森林效应估计相关系数低于 0.8 说明模型设定敏感协变量集加/减批次变量q 0.05 基因数变化数量差一倍以上要解释这张表建议直接贴进补充材料比任何文字描述都有说服力。审稿人问“你的结论稳健吗”把重叠率甩过去就行。最后一个技巧是日志。把sessionInfo()、cfg、每个脚本的运行时间和输入文件的行列数统一写进output/logs/一个脚本一个文件。这个习惯我是被坑出来的有次报告里的数和三个月前跑的对不上翻遍代码也没改过最后发现是上游给的表达矩阵悄悄更新了一版样本数从 512 变成了 508。从那以后我每次跑因果估计流水线都强制在第一步把sessionInfo()和输入文件的 md5 一起写进日志跑完再核对一遍才敢往报告里放数。这类工程的可靠性不在估计量写得多漂亮而在每一步都能回溯。希望帮到你。本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
