簡介SmokingGeneExpression 是一份面向生物信息學與統(tǒng)計分析學習者的 R 語言項目資料圍繞吸煙對健康肺組織基因表達的短期因果效應展開。內容將因果推斷思路落到基因表達場景涉及工具變量、傾向得分匹配、有向無環(huán)圖等混雜控制方法以及 RNA-seq 差異表達分析、PCA 可視化與生物學解釋適合希望理解觀察性數(shù)據(jù)因果分析流程并借助 R 復現(xiàn)實戰(zhàn)的讀者。壓縮包共 45 個文件約 3.1MB以 png 圖表、txt 結果說明、xlsx 與 csv 數(shù)據(jù)表為主另含 R 腳本和 pdf 文檔覆蓋分析代碼、數(shù)據(jù)樣例與結果呈現(xiàn)等用途。已有 148 人學習可作為課程作業(yè)、課題預研或方法對照的參考。讀者可從中獲得較完整的分析素材按目錄梳理研究邏輯、復現(xiàn)關鍵步驟并積累 R 語言實踐。1. 項目定位用 R 把吸煙對肺組織基因表達的因果效應算清楚一個做肺組織轉錄組的同行跟我吐槽過差異表達分析跑完吸煙組和從不吸煙組之間掛著三千多個顯著基因圖表也漂亮結果報告交上去第一個問題就是“你怎么知道這些差異不是吸煙者普遍年齡更大、BMI 更高造成的”。相關性分析能告訴你兩組人不一樣但不能告訴你是吸煙把某個基因的表達抬上去的。SmokingGeneExpression 這個 R 項目要解決的正是這一層把“吸煙”當作處理變量把“某個基因的表達量”當作結局用因果推斷的估計量把混雜因素壓下去給出每個基因在健康肺組織里對吸煙暴露的短期因果效應估計。適合誰來讀手上有隊列或公開表達資源、會寫基本 R 腳本、但還沒系統(tǒng)跑過傾向評分加權和雙重穩(wěn)健估計的人。項目的技術棧是純 R——數(shù)據(jù)層用 data.table 和基礎的 Bioconductor 表達處理估計層主要靠 WeightIt、grf、sandwich 這幾個包可視化走 ggplot2。它不訓練模型、不做深度學習整套流程跑完就是一張帶置信區(qū)間和 q 值的基因效應表外加幾張診斷圖。對做生信但被審稿意見追著問因果的人這個骨架的參考價值比多跑十個差異分析大。2. 工程結構R 腳本怎么分工、依賴怎么裝、第一次怎么跑通這個項目不是包package是一套按分析階段切開的腳本工程。這意味著它沒有R CMD check那套約束好處是改起來自由壞處是環(huán)境不一致時排錯特別費勁——這是后面避坑章節(jié)的重要伏筆。2.1 依賴清單哪些包是必需的哪些是可選的依賴分成三檔。第一檔是數(shù)據(jù)處理data.table、dplyr、tidyr、readr第二檔是因果估計核心WeightIt傾向評分與權重、grf因果森林與 AIPW 的機器學習 nuisance 估計、sandwich加lmtest聚類穩(wěn)健標準誤第三檔是多重檢驗與出圖qvalue、ggplot2、patchwork。并行那塊用futurefurrr規(guī)?；苋蚪M時會開關它。環(huán)境鎖定用renv項目根目錄有一份renv.lock。恢復環(huán)境的標準動作# 在項目根目錄執(zhí)行renv 會按 lock 文件把包版本拉回一致狀態(tài) Rscript -e if (!requireNamespace(renv, quietly TRUE)) install.packages(renv); renv::restore()這條命令的重點不是“安裝依賴”而是“安裝指定版本的依賴”。grf和WeightIt這兩個包在版本升級時函數(shù)簽名改過尤其是WeightIt的estimand參數(shù)和grf的honesty默認值跨版本跑出來的數(shù)會對不上。我第一次拿到這套腳本時偷懶沒跑 restore直接用系統(tǒng)里現(xiàn)裝的包結果 IPW 和 AIPW 的效應估計差了將近 8%排查了兩天才發(fā)現(xiàn)是權重歸一化方式變了。血淚經(jīng)驗這類因果估計工程包版本必須鎖死。如果機器上裝不了qvalueBioconductor 包退路是用p.adjust(method BH)兩者在大多數(shù)情況下結論一致只是 q 值的具體數(shù)值略有差異報告里注明用的是哪一種即可。2.2 目錄分層data、R、scripts、output 各管什么工程目錄大致是這樣切的路徑內容說明data/raw/原始表達矩陣、表型表、協(xié)變量表只讀任何腳本不得寫回data/derived/標準化后的表達、構建好的分析集中間產(chǎn)物可全部刪除重跑R/函數(shù)層io.R、preprocess.R、estimators.R、sensitivity.R、utils.R只定義函數(shù)不執(zhí)行分析scripts/執(zhí)行層01_prepare_data.R到05_report.R有編號按順序跑output/tables/基因效應表、診斷統(tǒng)計表交付物output/figures/森林圖、權重分布圖、QQ 圖交付物這個分層是有講究的。R/里的函數(shù)不碰全局變量參數(shù)從外面?zhèn)鬟M去方便單獨做單元測試scripts/里的腳本只負責“讀什么、調哪個函數(shù)、寫到哪里”。很多生信分析腳本寫著寫著就變成一個兩千行的巨型文件中間變量到處飛改一個參數(shù)要重跑全流程這種結構就是為了防這個。常見的做法是每個scripts/開頭統(tǒng)一 source 一遍函數(shù)層# scripts/01_prepare_data.R 的開頭固定寫這一塊 suppressPackageStartupMessages({ library(data.table) library(dplyr) }) # 統(tǒng)一從項目根目錄解析路徑避免依賴 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)) # 全局參數(shù)集中在這一處后面所有腳本讀同一個配置 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因為因果森林內部有隨機抽樣不固定的話兩次跑出來的效應估計在小數(shù)點后幾位會飄。路徑用here::here()而不是相對路徑是為了讓腳本從 RStudio、命令行、還是調度系統(tǒng)提交都能找到文件這個小習慣能省掉大量“在我機器上能跑”的扯皮。2.3 第一次跑通從 01 到 05 的執(zhí)行順序五個腳本是一條串行鏈但中間結果都落盤所以可以從任意一步重跑# 全流程按編號順序執(zhí)行 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做數(shù)據(jù)對齊和表達預處理輸出data/derived/analysis_set.rds02和03分別在同一個分析集上跑 IPW 和 AIPW各自產(chǎn)出一張基因級效應表04只在候選基因子集上跑因果森林用來觀察異質性05把三張表按基因 ID 合并、統(tǒng)一多重檢驗校正、出圖和寫報告。第一次跑建議先在01里把基因數(shù)砍到幾百個做冒煙測試確認整條鏈沒斷再放開全量?；驍?shù)上萬的表達矩陣走完全套單核大概要幾十分鐘到幾小時取決于樣本量和是否開并行。冒煙測試這一步別省04的因果森林在小樣本上最容易因為某個基因方差為零直接拋錯早發(fā)現(xiàn)早改。3. 數(shù)據(jù)層拆解暴露、結局與協(xié)變量怎么對齊到同一張表因果推斷里最容易翻車的從來不是估計量寫錯而是數(shù)據(jù)對齊錯了。這個項目的核心數(shù)據(jù)結構只有一張寬表每一行是一個樣本列包括樣本 ID、吸煙暴露、幾十上百個協(xié)變量、以及每個基因的表達值。所有估計函數(shù)都只認這張表保證口徑一致。3.1 吸煙暴露變量的編碼二值化之前的三個判斷原始表型里吸煙狀態(tài)通常是分好幾檔的從不、曾經(jīng)、當前、偶爾、已戒若干年甚至還有每天支數(shù)。直接二值化會丟信息但項目走的是“當前吸煙 vs 從不吸煙”這條二值路線理由有三短期因果效應下“曾經(jīng)吸煙”人群的暴露已經(jīng)停止把它塞進處理組會稀釋效應每天支數(shù)在多數(shù)隊列里缺失率太高做連續(xù)暴露需要另一套劑量-反應模型而二值暴露能直接套用傾向評分框架診斷工具也成熟。編碼這一步要把“曾經(jīng)吸煙”和缺失顯式剔出去而不是當成對照組# 把原始吸煙編碼映射成 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)), # 只保留兩種明確狀態(tài)其余一律置 NA W case_when( smoke_status %in% c(current, current smoker, daily) ~ 1L, smoke_status %in% c(never, never smoker) ~ 0L, TRUE ~ NA_integer_ ), # 單獨打個標記方便統(tǒng)計各剔除原因的數(shù)量 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%里放了多種寫法是為了兼容不同來源的表型表真實數(shù)據(jù)里 “Current Smoker” 和 “current” 混著出現(xiàn)是常態(tài)。drop_reason這一列看起來多余但它是后面寫方法學部分的證據(jù)——審稿人一定會問你為什么最終樣本是 N你得能答出“剔除了多少既往吸煙者、多少狀態(tài)缺失”。attr掛上去的計數(shù)會被01腳本打印到日志里養(yǎng)成這個習慣能省掉回頭重算的麻煩。還有一層容易忽略的暴露的時間窗。項目說的是“短期”那就得確認表型采集時間和組織取樣時間接近如果表型是問卷回顧性的“過去一年是否吸煙”而樣本是五年前存的這個短期假設就不成立。數(shù)據(jù)里如果有采樣日期字段務必在預處理階段做個時間差統(tǒng)計。3.2 表達矩陣的預處理過濾、標準化、以及要不要去批次表達矩陣進模型前要過三道關。第一道是低表達過濾用min_expr和min_samples兩個參數(shù)控制目的是減少檢驗數(shù)量、降低多重校正的懲罰。第二道是分布變換如果是原始 count 就走log2(CPM 1)或vst如果上游已經(jīng)給了 log2 尺度就跳過。第三道是批次這一步是分水嶺批次如果和吸煙狀態(tài)強相關去批次會把效應一起減掉如果批次獨立于暴露去不去影響都不大。保守做法是先檢查批次與暴露的關聯(lián)再做決定# 過濾 標準化并在返回前檢查批次是否與暴露混雜 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 變換已經(jīng)在 log 尺度上的矩陣跳過 if (max(mat, na.rm TRUE) 50) mat - log2(mat 1) # 3) 分位數(shù)標準化抹平樣本間整體分布差異 mat - limma::normalizeBetweenArrays(mat, method quantile) # 4) 批次與暴露的關聯(lián)檢查卡方檢驗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(批次-暴露關聯(lián)檢驗 p %.4g, p_batch)) if (!is.na(p_batch) p_batch 0.05) { warning(批次與吸煙狀態(tài)顯著關聯(lián)去批次可能吸收掉部分處理效應慎用 ComBat) } } list(expr mat, n_gene nrow(mat), n_sample ncol(mat)) }rowMeans(mat cfg$min_expr)這一行是過濾的核心表達式返回每個基因超過閾值的樣本比例再和min_samples比。normalizeBetweenArrays走 quantile 方法對偏態(tài)分布的表達值比中位數(shù)標準化更穩(wěn)。最值得看的是第 4 段——批次卡方檢驗。如果 p 值很小說明這批數(shù)據(jù)里“哪臺機器測的”和“這個人吸不吸煙”高度相關此時上 ComBat 或limma::removeBatchEffect等于把處理效應的一部分當成批次削掉了后面估計出來的效應會系統(tǒng)性偏小。這種情況我一般不去批次而是把批次當作協(xié)變量放進傾向評分模型讓加權去平衡它而不是從數(shù)據(jù)里刪掉它。矩陣的樣本順序和表型的樣本順序必須嚴格一致這是最容易出錯的地方。prep_expression返回前應該強制按meta$sample_id重排或者用stopifnot(identical(colnames(mat), meta$sample_id))做一次斷言別指望兩邊碰巧一致。3.3 協(xié)變量集與設計矩陣共線性檢查協(xié)變量集的選擇原則是只放暴露之前就已經(jīng)確定的變量年齡、性別、遺傳背景主成分、測序深度、批次不要放可能被暴露影響的中介變量比如某些炎癥指標、肺功能。放中介進去會把總效應部分遮掉估出來的是直接效應解釋起來是另一回事。確定變量集之后必須檢查設計矩陣的秩。協(xié)變量之間完全共線比如兩個高度相關的技術指標或者某個分類變量某一水平只有一個人會讓傾向評分模型直接崩掉glm可能不報錯但給出荒謬的系數(shù)。# 構建協(xié)變量設計矩陣并做秩與稀疏水平檢查 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(設計矩陣秩虧存在完全共線協(xié)變量需手動剔除) # 分類變量稀疏水平某個水平樣本過少會讓該水平的傾向評分極端 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小于列數(shù)就是秩虧最常見的成因是分類變量做了獨熱編碼后某一列全是零或者兩個連續(xù)變量做了線性組合。tb 20這個閾值不是硬標準20 只是個提醒線真正要看的是該水平在兩組里的分布——如果某個批次只有 3 個人且全在吸煙組傾向評分在這個批次上會趨近 1加權時必然出極端權重這就是第 5 章要展開的問題。提前在這里發(fā)現(xiàn)比跑到加權階段再看權重分布圖省事得多。4. 因果估計核心三種估計量在單個基因上的完整鏈路數(shù)據(jù)對齊之后估計層的輸入就是一個向量結局Y某個基因的表達值、處理W0/1、協(xié)變量矩陣X。項目給三種估計量不是讓你三選一而是讓它們互相印證IPW 看加權后的可比人群AIPW 用兩個模型互相兜底因果森林看效應是不是在亞群里不一致。三者結論方向一致才敢往報告里寫。4.1 傾向評分與 IPW把不可比的人群“造”成可比傾向評分e(X) P(W1|X)的直覺是給每個樣本算一個“被吸煙的概率”然后用這個概率的倒數(shù)做權重讓協(xié)變量分布在兩組之間拉平。吸煙概率高的人如果實際沒吸煙權重就大代表他“替”了很多吸煙者。實現(xiàn)上用WeightIt它比手寫glm多了一層診斷# 單個基因的 IPW 估計輸入是結局向量 處理向量 協(xié)變量數(shù)據(jù)框 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 拿穩(wěn)健標準誤 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的區(qū)別要拎清ATE 估的是“如果全人群都吸煙 vs 全人群都不吸煙”的平均差異ATT 只估“吸煙者如果沒吸煙會怎樣”。做暴露危害評估一般用 ATE。stabilize TRUE把權重除以處理組的邊際概率能顯著降低權重的方差。overlap里那兩個n_lt01/n_gt99計數(shù)是診斷的核心——如果有一批人傾向評分大于 0.99說明模型認為他們“幾乎肯定吸煙”而那批人恰好沒吸煙權重會飆到 100 以上一個樣本就能主導整個估計。trim 0.01是在 1% 和 99% 分位處截尾代價是引入輕微偏差但換來方差可控實踐中這筆交易通常劃算。ess有效樣本量是另一個必看指標sum(w)^2 / sum(w^2)給出加權后相當于多少“等效樣本”。如果原始 500 例加權后 ESS 只剩 60說明權重分布極其不均這個估計不該信。我的經(jīng)驗閾值是 ESS 低于原始的 30% 就要回頭查協(xié)變量集。4.2 AIPW 雙重穩(wěn)健兩個模型互相兜底IPW 只依賴傾向評分模型正確一旦模型設定錯了比如該放交互項沒放估計就有偏。AIPW 加了結局模型兩個模型只要有一個對估計就一致這叫雙重穩(wěn)健。代價是計算量大一些但在這個場景里值得。結局模型可以用簡單的線性回歸也可以用grf的回歸森林。項目里默認走線性因為基因表達和協(xié)變量的關系大多近似線性森林在小樣本上反而方差大# 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 影響函數(shù)逐樣本取值 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]是必須的現(xiàn)實數(shù)據(jù)里fitted(ps_fit)出現(xiàn) 1.0000 是常見的除法直接炸出 Inf。extrap那幾個極值是我加的診斷看mu1的最大值有沒有遠超實際觀測到的表達量范圍——如果超過了說明結局模型在做荒謬的外推這條估計不能采信先回頭檢查協(xié)變量。use_forest TRUE那條分支用grf::regression_forest把W當成一個普通特征塞進X預測時分別把W置 1 和 0。這種寫法比訓練兩個獨立森林更高效但對樣本量要求更高小于 200 例時慎用。4.3 因果森林平均效應之外的異質性IPW 和 AIPW 給的是全人群平均效應但吸煙對基因表達的影響很可能在亞群里不一致——比如年齡大的人反應更強或者只在不吸煙的遺傳背景子群里顯現(xiàn)。因果森林直接估每個樣本的條件平均處理效應tau(X)然后看tau有沒有系統(tǒng)性的分化。# 因果森林估條件平均處理效應看異質性 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, # 誠實劃分避免用同一批數(shù)據(jù)既分裂又估值 min.node.size 10, # 葉子最小樣本太小會過擬合 sample.fraction 0.5, # 子采樣比例 seed 20240101 ) # 全人群平均效應target.sample overlap 只在有重疊的人群上平均 ate - grf::average_treatment_effect(cf, target.sample overlap) # 變量重要性看哪些協(xié)變量在驅動異質性 vimp - grf::variable_importance(cf) list(cf cf, ate ate, vimp vimp) }honesty TRUE是grf的默認值作用是樣本分裂和效應估計用不同的子樣本防止過擬合把噪聲當成異質性。這個參數(shù)不要關掉關了之后tau的分布會變得極其“好看”但全是假的。min.node.size 10是葉子節(jié)點的最小樣本量基因表達數(shù)據(jù)的噪聲大設太小會切出一堆只有幾個樣本的葉子估出來的tau方差爆炸。target.sample overlap這一項值得單獨說。默認是all在全員上平均改成overlap只在傾向評分有重疊的區(qū)域上平均更穩(wěn)健但也更保守。做暴露效應的基因篩查時用overlap結論不容易被少數(shù)極端個體帶跑。variable_importance的輸出是一個矩陣第一列是原始重要性第二列是歸一化到 0-1 的版本?？辞拔迕蛪颉绻挲g或某個技術協(xié)變量排第一說明異質性主要來自技術因素而非生物學這時候要警惕別急著往機制上解釋。在基因層面因果森林只在候選基因子集上跑不是全基因組。原因很實際2000 棵樹乘以幾百個基因單機跑一晚上都不一定完而且大多數(shù)基因的異質性檢驗并不顯著跑全量是浪費算力。4.4 從 p 值到 q 值結果表的組裝三種估計跑完每個基因會拿到三套est和se。組裝結果表要解決兩個問題一是統(tǒng)一效應方向基因表達上調還是下調二是多重檢驗。基因數(shù)量上千直接拿 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 作為主估計它雙重穩(wěn)健 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這一列是經(jīng)驗產(chǎn)物。IPW 和 AIPW 的方向如果相反說明兩個模型的假設沖突得很厲害可能是傾向評分模型設定有問題也可能是結局模型的線性假設在這個基因上不成立。這類基因即使 q 值很小也不該寫進結論最多放進補充材料標注方法間不一致等著用敏感性分析去解釋。主估計選 AIPW 而不是 IPW理由是雙重穩(wěn)健只要傾向評分和結局模型有一個對估計就一致比單靠一個模型的 IPW 更扛得住設定誤差。5. 避坑清單因果推斷流水線里最容易翻車的五處這套流程跑通不難難的是跑出來的數(shù)能不能信。下面五條是我在這類項目上反復踩到的每一條都按“現(xiàn)象 → 原因 → 解決”整理配上排查代碼。5.1 權重爆炸與重疊性缺失現(xiàn)象IPW 估出來的效應比 AIPW 大好幾倍標準誤也大得離譜ESS只剩原始樣本的 10% 不到。原因傾向評分模型太“自信”給某些樣本賦了接近 0 或 1 的概率。常見于某個協(xié)變量幾乎能完全預測處理狀態(tài)——比如某一批次的所有樣本都來自吸煙組模型在這個批次上學到“批次 吸煙”權重自然爆掉。另一個成因是分類變量的稀疏水平某一類只有幾個人邏輯回歸在該維度上給出極端的系數(shù)。解決先看重疊性診斷再決定是修剪樣本還是截尾權重。# 權重診斷分布、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) }打印極端權重樣本的協(xié)變量是定位問題最快的路徑。如果這十個樣本全在同一個批次那就是批次導致的如果全是高齡那是年齡分布尾部的問題。ESS低于 30% 時我一般先嘗試trim 0.01還不行就改用重疊權重WeightIt里method glm, estimand ATO重疊權重對極端傾向評分不敏感代價是估計目標變成“重疊人群的效應”解釋時要寫清楚。5.2 批次、共線性與“靜默”的樣本錯位現(xiàn)象腳本一聲不吭跑完結果表也有但效應估計整體偏小或者某個已知的強效應基因死活不顯著。原因三種可能。批次與暴露強相關去批次時把處理效應一起減掉了協(xié)變量共線導致glm系數(shù)不穩(wěn)定傾向評分本身就不準表達矩陣和表型的樣本順序沒對齊merge或cbind之后列和行的對應關系已經(jīng)亂了但 R 不會報錯。解決三個檢查必須都在預處理階段做完——批次卡方檢驗3.2 節(jié)、設計矩陣秩檢查3.3 節(jié)、樣本 ID 斷言。# 樣本對齊的最后一道保險放在所有合并操作之后 assert_alignment - function(expr, meta) { stopifnot( 表達矩陣列名缺失 !is.null(colnames(expr)), 表型樣本 ID 缺失 sample_id %in% names(meta), 樣本數(shù)不一致 ncol(expr) nrow(meta), 樣本順序錯位 identical(colnames(expr), as.character(meta$sample_id)) ) invisible(TRUE) }identical這一條是硬斷言不做容錯順序不一致直接拋錯。我寧可讓腳本在第一步就炸也不想在寫了兩周報告之后發(fā)現(xiàn)樣本對錯了。5.3 因果森林在小樣本上的“看起來很準”現(xiàn)象因果森林跑出來的tau分布特別漂亮亞群間差異明顯但換一個隨機種子結果就變了。原因honesty FALSE或者min.node.size設得太小樹在噪聲上過擬合?；虮磉_數(shù)據(jù)本身噪聲大樣本量幾百例時2000 棵樹的平均還壓不住方差tau的分布就是隨機的。解決三個參數(shù)一起調——保持honesty TRUEmin.node.size不低于 10sample.fraction取 0.4 到 0.6 之間然后用兩個不同種子各跑一遍看tau的相關系數(shù)。# 穩(wěn)定性檢驗換種子重跑看條件效應是否可復現(xiàn) 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 相關系數(shù) %.3f\n, r)) if (r 0.7) warning(條件效應不穩(wěn)定異質性結論不可靠) invisible(r) }相關系數(shù)低于 0.7 就別寫異質性結論了老老實實報平均效應。這個檢驗成本不高值得每個候選基因都過一遍。5.4 多重檢驗的兩種誤用現(xiàn)象q 值算出來了顯著基因上千個但和已知的吸煙相關基因對不上。原因要么用了 Bonferroni在基因表達場景過于嚴苛真正的信號被壓掉要么在篩選和驗證之間反復調閾值最后報的是“調出來的”顯著性。還有一種隱蔽的IPW 和 AIPW 各算一遍 q 值挑顯著的那個報這是變相的多重比較。解決主分析只跑一套AIPWq 值只算一次閾值預先定死比如 q 0.05敏感性分析的結果單獨列不參與篩選。方法學部分要寫清楚“主估計用的是哪一個”。5.5 用同一份數(shù)據(jù)調參又評估現(xiàn)象所有診斷都完美重疊性好、ESS 高、q 值分布漂亮但換一個獨立隊列就復現(xiàn)不出來。原因傾向評分模型的協(xié)變量集、trim閾值、因果森林的參數(shù)如果在看到結果之后反復調整等于把結局信息泄漏進了模型選擇。這類偏差不會體現(xiàn)在任何一張診斷圖上。解決把樣本按 7:3 隨機切開前 70% 用來定協(xié)變量集和超參數(shù)后 30% 只跑一次不看結果不回頭。如果樣本量不夠切就在分析計劃里把參數(shù)全部預先寫死中途改了要記錄在日志里。# 預注冊式的參數(shù)凍結把最終參數(shù)寫進文件后續(xù)腳本只讀不猜 freeze_params - function(cfg, path output/tables/frozen_params.rds) { saveRDS(cfg, path) message(參數(shù)已凍結至 , path, 后續(xù)分析請勿修改) }這條看著像流程潔癖但它是唯一能防住“調參調到顯著”的辦法。玄學的地方在于你自己調的時候完全不覺得在作弊每一步都有“合理”的理由。6. 二次開發(fā)把單基因流程改成批量篩選與可復現(xiàn)報告拿到這套腳本的人十有八九第一件事是想跑全基因組。單基因函數(shù)寫好了批量跑不是簡單套個for循環(huán)——直接循環(huán)會撞上三個問題內存、算力和失敗恢復。我給的做法是把估計函數(shù)包一層加上并行和緩存。先看批量層怎么寫# 批量估計并行 單基因失敗不影響整體 結果落盤緩存 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 的隨機數(shù)流可復現(xiàn)不加這一句因果森林那部分每次跑出來都不一樣。tryCatch把單基因的失敗隔離掉glm在完全分離的數(shù)據(jù)上會拋 warning 甚至 error一個基因炸掉不該拖垮整批。零方差判斷放在最前面是最常見的崩潰點——某些基因在所有樣本里表達值都一樣lm直接返回 NA 系數(shù)。n_workers不建議超過物理核心數(shù)減一grf內部本身也會占內存開太多反而會因為內存交換變慢。第二個改造點是報告層。原始腳本輸出的是 CSV二次開發(fā)可以把它變成參數(shù)化的 Rmd 或 Quarto 文檔把 q 值分布、火山圖、前若干個基因的森林圖一次性生成換一批數(shù)據(jù)只改cfg里的路徑。這里的關鍵是把“結果表”和“圖”解耦圖只讀表不重算估計。第三個改造點是加一層敏感性分析。做法很簡單把trim從 0.01 換成 0.05 和 0.1各跑一遍 AIPW看顯著基因的重疊比例。敏感性維度參數(shù)觀察指標警戒線權重截尾trim∈ {0.01, 0.05, 0.1}顯著基因重疊率低于 70% 說明結論依賴閾值結局模型線性 vs 回歸森林效應估計相關系數(shù)低于 0.8 說明模型設定敏感協(xié)變量集加/減批次變量q 0.05 基因數(shù)變化數(shù)量差一倍以上要解釋這張表建議直接貼進補充材料比任何文字描述都有說服力。審稿人問“你的結論穩(wěn)健嗎”把重疊率甩過去就行。最后一個技巧是日志。把sessionInfo()、cfg、每個腳本的運行時間和輸入文件的行列數(shù)統(tǒng)一寫進output/logs/一個腳本一個文件。這個習慣我是被坑出來的有次報告里的數(shù)和三個月前跑的對不上翻遍代碼也沒改過最后發(fā)現(xiàn)是上游給的表達矩陣悄悄更新了一版樣本數(shù)從 512 變成了 508。從那以后我每次跑因果估計流水線都強制在第一步把sessionInfo()和輸入文件的 md5 一起寫進日志跑完再核對一遍才敢往報告里放數(shù)。這類工程的可靠性不在估計量寫得多漂亮而在每一步都能回溯。希望幫到你。本文還有配套的精品資源點擊獲取 R語言因果推斷：吸煙對肺組織基因表達的IPW與AIPW分析
簡介SmokingGeneExpression 是一份面向生物信息學與統(tǒng)計分析學習者的 R 語言項目資料圍繞吸煙對健康肺組織基因表達的短期因果效應展開。內容將因果推斷思路落到基因表達場景涉及工具變量、傾向得分匹配、有向無環(huán)圖等混雜控制方法以及 RNA-seq 差異表達分析、PCA 可視化與生物學解釋適合希望理解觀察性數(shù)據(jù)因果分析流程并借助 R 復現(xiàn)實戰(zhàn)的讀者。壓縮包共 45 個文件約 3.1MB以 png 圖表、txt 結果說明、xlsx 與 csv 數(shù)據(jù)表為主另含 R 腳本和 pdf 文檔覆蓋分析代碼、數(shù)據(jù)樣例與結果呈現(xiàn)等用途。已有 148 人學習可作為課程作業(yè)、課題預研或方法對照的參考。讀者可從中獲得較完整的分析素材按目錄梳理研究邏輯、復現(xiàn)關鍵步驟并積累 R 語言實踐。1. 項目定位用 R 把吸煙對肺組織基因表達的因果效應算清楚一個做肺組織轉錄組的同行跟我吐槽過差異表達分析跑完吸煙組和從不吸煙組之間掛著三千多個顯著基因圖表也漂亮結果報告交上去第一個問題就是“你怎么知道這些差異不是吸煙者普遍年齡更大、BMI 更高造成的”。相關性分析能告訴你兩組人不一樣但不能告訴你是吸煙把某個基因的表達抬上去的。SmokingGeneExpression 這個 R 項目要解決的正是這一層把“吸煙”當作處理變量把“某個基因的表達量”當作結局用因果推斷的估計量把混雜因素壓下去給出每個基因在健康肺組織里對吸煙暴露的短期因果效應估計。適合誰來讀手上有隊列或公開表達資源、會寫基本 R 腳本、但還沒系統(tǒng)跑過傾向評分加權和雙重穩(wěn)健估計的人。項目的技術棧是純 R——數(shù)據(jù)層用 data.table 和基礎的 Bioconductor 表達處理估計層主要靠 WeightIt、grf、sandwich 這幾個包可視化走 ggplot2。它不訓練模型、不做深度學習整套流程跑完就是一張帶置信區(qū)間和 q 值的基因效應表外加幾張診斷圖。對做生信但被審稿意見追著問因果的人這個骨架的參考價值比多跑十個差異分析大。2. 工程結構R 腳本怎么分工、依賴怎么裝、第一次怎么跑通這個項目不是包package是一套按分析階段切開的腳本工程。這意味著它沒有R CMD check那套約束好處是改起來自由壞處是環(huán)境不一致時排錯特別費勁——這是后面避坑章節(jié)的重要伏筆。2.1 依賴清單哪些包是必需的哪些是可選的依賴分成三檔。第一檔是數(shù)據(jù)處理data.table、dplyr、tidyr、readr第二檔是因果估計核心WeightIt傾向評分與權重、grf因果森林與 AIPW 的機器學習 nuisance 估計、sandwich加lmtest聚類穩(wěn)健標準誤第三檔是多重檢驗與出圖qvalue、ggplot2、patchwork。并行那塊用futurefurrr規(guī)?；苋蚪M時會開關它。環(huán)境鎖定用renv項目根目錄有一份renv.lock。恢復環(huán)境的標準動作# 在項目根目錄執(zhí)行renv 會按 lock 文件把包版本拉回一致狀態(tài) Rscript -e if (!requireNamespace(renv, quietly TRUE)) install.packages(renv); renv::restore()這條命令的重點不是“安裝依賴”而是“安裝指定版本的依賴”。grf和WeightIt這兩個包在版本升級時函數(shù)簽名改過尤其是WeightIt的estimand參數(shù)和grf的honesty默認值跨版本跑出來的數(shù)會對不上。我第一次拿到這套腳本時偷懶沒跑 restore直接用系統(tǒng)里現(xiàn)裝的包結果 IPW 和 AIPW 的效應估計差了將近 8%排查了兩天才發(fā)現(xiàn)是權重歸一化方式變了。血淚經(jīng)驗這類因果估計工程包版本必須鎖死。如果機器上裝不了qvalueBioconductor 包退路是用p.adjust(method BH)兩者在大多數(shù)情況下結論一致只是 q 值的具體數(shù)值略有差異報告里注明用的是哪一種即可。2.2 目錄分層data、R、scripts、output 各管什么工程目錄大致是這樣切的路徑內容說明data/raw/原始表達矩陣、表型表、協(xié)變量表只讀任何腳本不得寫回data/derived/標準化后的表達、構建好的分析集中間產(chǎn)物可全部刪除重跑R/函數(shù)層io.R、preprocess.R、estimators.R、sensitivity.R、utils.R只定義函數(shù)不執(zhí)行分析scripts/執(zhí)行層01_prepare_data.R到05_report.R有編號按順序跑output/tables/基因效應表、診斷統(tǒng)計表交付物output/figures/森林圖、權重分布圖、QQ 圖交付物這個分層是有講究的。R/里的函數(shù)不碰全局變量參數(shù)從外面?zhèn)鬟M去方便單獨做單元測試scripts/里的腳本只負責“讀什么、調哪個函數(shù)、寫到哪里”。很多生信分析腳本寫著寫著就變成一個兩千行的巨型文件中間變量到處飛改一個參數(shù)要重跑全流程這種結構就是為了防這個。常見的做法是每個scripts/開頭統(tǒng)一 source 一遍函數(shù)層# scripts/01_prepare_data.R 的開頭固定寫這一塊 suppressPackageStartupMessages({ library(data.table) library(dplyr) }) # 統(tǒng)一從項目根目錄解析路徑避免依賴 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)) # 全局參數(shù)集中在這一處后面所有腳本讀同一個配置 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因為因果森林內部有隨機抽樣不固定的話兩次跑出來的效應估計在小數(shù)點后幾位會飄。路徑用here::here()而不是相對路徑是為了讓腳本從 RStudio、命令行、還是調度系統(tǒng)提交都能找到文件這個小習慣能省掉大量“在我機器上能跑”的扯皮。2.3 第一次跑通從 01 到 05 的執(zhí)行順序五個腳本是一條串行鏈但中間結果都落盤所以可以從任意一步重跑# 全流程按編號順序執(zhí)行 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做數(shù)據(jù)對齊和表達預處理輸出data/derived/analysis_set.rds02和03分別在同一個分析集上跑 IPW 和 AIPW各自產(chǎn)出一張基因級效應表04只在候選基因子集上跑因果森林用來觀察異質性05把三張表按基因 ID 合并、統(tǒng)一多重檢驗校正、出圖和寫報告。第一次跑建議先在01里把基因數(shù)砍到幾百個做冒煙測試確認整條鏈沒斷再放開全量?；驍?shù)上萬的表達矩陣走完全套單核大概要幾十分鐘到幾小時取決于樣本量和是否開并行。冒煙測試這一步別省04的因果森林在小樣本上最容易因為某個基因方差為零直接拋錯早發(fā)現(xiàn)早改。3. 數(shù)據(jù)層拆解暴露、結局與協(xié)變量怎么對齊到同一張表因果推斷里最容易翻車的從來不是估計量寫錯而是數(shù)據(jù)對齊錯了。這個項目的核心數(shù)據(jù)結構只有一張寬表每一行是一個樣本列包括樣本 ID、吸煙暴露、幾十上百個協(xié)變量、以及每個基因的表達值。所有估計函數(shù)都只認這張表保證口徑一致。3.1 吸煙暴露變量的編碼二值化之前的三個判斷原始表型里吸煙狀態(tài)通常是分好幾檔的從不、曾經(jīng)、當前、偶爾、已戒若干年甚至還有每天支數(shù)。直接二值化會丟信息但項目走的是“當前吸煙 vs 從不吸煙”這條二值路線理由有三短期因果效應下“曾經(jīng)吸煙”人群的暴露已經(jīng)停止把它塞進處理組會稀釋效應每天支數(shù)在多數(shù)隊列里缺失率太高做連續(xù)暴露需要另一套劑量-反應模型而二值暴露能直接套用傾向評分框架診斷工具也成熟。編碼這一步要把“曾經(jīng)吸煙”和缺失顯式剔出去而不是當成對照組# 把原始吸煙編碼映射成 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)), # 只保留兩種明確狀態(tài)其余一律置 NA W case_when( smoke_status %in% c(current, current smoker, daily) ~ 1L, smoke_status %in% c(never, never smoker) ~ 0L, TRUE ~ NA_integer_ ), # 單獨打個標記方便統(tǒng)計各剔除原因的數(shù)量 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%里放了多種寫法是為了兼容不同來源的表型表真實數(shù)據(jù)里 “Current Smoker” 和 “current” 混著出現(xiàn)是常態(tài)。drop_reason這一列看起來多余但它是后面寫方法學部分的證據(jù)——審稿人一定會問你為什么最終樣本是 N你得能答出“剔除了多少既往吸煙者、多少狀態(tài)缺失”。attr掛上去的計數(shù)會被01腳本打印到日志里養(yǎng)成這個習慣能省掉回頭重算的麻煩。還有一層容易忽略的暴露的時間窗。項目說的是“短期”那就得確認表型采集時間和組織取樣時間接近如果表型是問卷回顧性的“過去一年是否吸煙”而樣本是五年前存的這個短期假設就不成立。數(shù)據(jù)里如果有采樣日期字段務必在預處理階段做個時間差統(tǒng)計。3.2 表達矩陣的預處理過濾、標準化、以及要不要去批次表達矩陣進模型前要過三道關。第一道是低表達過濾用min_expr和min_samples兩個參數(shù)控制目的是減少檢驗數(shù)量、降低多重校正的懲罰。第二道是分布變換如果是原始 count 就走log2(CPM 1)或vst如果上游已經(jīng)給了 log2 尺度就跳過。第三道是批次這一步是分水嶺批次如果和吸煙狀態(tài)強相關去批次會把效應一起減掉如果批次獨立于暴露去不去影響都不大。保守做法是先檢查批次與暴露的關聯(lián)再做決定# 過濾 標準化并在返回前檢查批次是否與暴露混雜 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 變換已經(jīng)在 log 尺度上的矩陣跳過 if (max(mat, na.rm TRUE) 50) mat - log2(mat 1) # 3) 分位數(shù)標準化抹平樣本間整體分布差異 mat - limma::normalizeBetweenArrays(mat, method quantile) # 4) 批次與暴露的關聯(lián)檢查卡方檢驗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(批次-暴露關聯(lián)檢驗 p %.4g, p_batch)) if (!is.na(p_batch) p_batch 0.05) { warning(批次與吸煙狀態(tài)顯著關聯(lián)去批次可能吸收掉部分處理效應慎用 ComBat) } } list(expr mat, n_gene nrow(mat), n_sample ncol(mat)) }rowMeans(mat cfg$min_expr)這一行是過濾的核心表達式返回每個基因超過閾值的樣本比例再和min_samples比。normalizeBetweenArrays走 quantile 方法對偏態(tài)分布的表達值比中位數(shù)標準化更穩(wěn)。最值得看的是第 4 段——批次卡方檢驗。如果 p 值很小說明這批數(shù)據(jù)里“哪臺機器測的”和“這個人吸不吸煙”高度相關此時上 ComBat 或limma::removeBatchEffect等于把處理效應的一部分當成批次削掉了后面估計出來的效應會系統(tǒng)性偏小。這種情況我一般不去批次而是把批次當作協(xié)變量放進傾向評分模型讓加權去平衡它而不是從數(shù)據(jù)里刪掉它。矩陣的樣本順序和表型的樣本順序必須嚴格一致這是最容易出錯的地方。prep_expression返回前應該強制按meta$sample_id重排或者用stopifnot(identical(colnames(mat), meta$sample_id))做一次斷言別指望兩邊碰巧一致。3.3 協(xié)變量集與設計矩陣共線性檢查協(xié)變量集的選擇原則是只放暴露之前就已經(jīng)確定的變量年齡、性別、遺傳背景主成分、測序深度、批次不要放可能被暴露影響的中介變量比如某些炎癥指標、肺功能。放中介進去會把總效應部分遮掉估出來的是直接效應解釋起來是另一回事。確定變量集之后必須檢查設計矩陣的秩。協(xié)變量之間完全共線比如兩個高度相關的技術指標或者某個分類變量某一水平只有一個人會讓傾向評分模型直接崩掉glm可能不報錯但給出荒謬的系數(shù)。# 構建協(xié)變量設計矩陣并做秩與稀疏水平檢查 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(設計矩陣秩虧存在完全共線協(xié)變量需手動剔除) # 分類變量稀疏水平某個水平樣本過少會讓該水平的傾向評分極端 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小于列數(shù)就是秩虧最常見的成因是分類變量做了獨熱編碼后某一列全是零或者兩個連續(xù)變量做了線性組合。tb 20這個閾值不是硬標準20 只是個提醒線真正要看的是該水平在兩組里的分布——如果某個批次只有 3 個人且全在吸煙組傾向評分在這個批次上會趨近 1加權時必然出極端權重這就是第 5 章要展開的問題。提前在這里發(fā)現(xiàn)比跑到加權階段再看權重分布圖省事得多。4. 因果估計核心三種估計量在單個基因上的完整鏈路數(shù)據(jù)對齊之后估計層的輸入就是一個向量結局Y某個基因的表達值、處理W0/1、協(xié)變量矩陣X。項目給三種估計量不是讓你三選一而是讓它們互相印證IPW 看加權后的可比人群AIPW 用兩個模型互相兜底因果森林看效應是不是在亞群里不一致。三者結論方向一致才敢往報告里寫。4.1 傾向評分與 IPW把不可比的人群“造”成可比傾向評分e(X) P(W1|X)的直覺是給每個樣本算一個“被吸煙的概率”然后用這個概率的倒數(shù)做權重讓協(xié)變量分布在兩組之間拉平。吸煙概率高的人如果實際沒吸煙權重就大代表他“替”了很多吸煙者。實現(xiàn)上用WeightIt它比手寫glm多了一層診斷# 單個基因的 IPW 估計輸入是結局向量 處理向量 協(xié)變量數(shù)據(jù)框 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 拿穩(wěn)健標準誤 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的區(qū)別要拎清ATE 估的是“如果全人群都吸煙 vs 全人群都不吸煙”的平均差異ATT 只估“吸煙者如果沒吸煙會怎樣”。做暴露危害評估一般用 ATE。stabilize TRUE把權重除以處理組的邊際概率能顯著降低權重的方差。overlap里那兩個n_lt01/n_gt99計數(shù)是診斷的核心——如果有一批人傾向評分大于 0.99說明模型認為他們“幾乎肯定吸煙”而那批人恰好沒吸煙權重會飆到 100 以上一個樣本就能主導整個估計。trim 0.01是在 1% 和 99% 分位處截尾代價是引入輕微偏差但換來方差可控實踐中這筆交易通常劃算。ess有效樣本量是另一個必看指標sum(w)^2 / sum(w^2)給出加權后相當于多少“等效樣本”。如果原始 500 例加權后 ESS 只剩 60說明權重分布極其不均這個估計不該信。我的經(jīng)驗閾值是 ESS 低于原始的 30% 就要回頭查協(xié)變量集。4.2 AIPW 雙重穩(wěn)健兩個模型互相兜底IPW 只依賴傾向評分模型正確一旦模型設定錯了比如該放交互項沒放估計就有偏。AIPW 加了結局模型兩個模型只要有一個對估計就一致這叫雙重穩(wěn)健。代價是計算量大一些但在這個場景里值得。結局模型可以用簡單的線性回歸也可以用grf的回歸森林。項目里默認走線性因為基因表達和協(xié)變量的關系大多近似線性森林在小樣本上反而方差大# 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 影響函數(shù)逐樣本取值 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]是必須的現(xiàn)實數(shù)據(jù)里fitted(ps_fit)出現(xiàn) 1.0000 是常見的除法直接炸出 Inf。extrap那幾個極值是我加的診斷看mu1的最大值有沒有遠超實際觀測到的表達量范圍——如果超過了說明結局模型在做荒謬的外推這條估計不能采信先回頭檢查協(xié)變量。use_forest TRUE那條分支用grf::regression_forest把W當成一個普通特征塞進X預測時分別把W置 1 和 0。這種寫法比訓練兩個獨立森林更高效但對樣本量要求更高小于 200 例時慎用。4.3 因果森林平均效應之外的異質性IPW 和 AIPW 給的是全人群平均效應但吸煙對基因表達的影響很可能在亞群里不一致——比如年齡大的人反應更強或者只在不吸煙的遺傳背景子群里顯現(xiàn)。因果森林直接估每個樣本的條件平均處理效應tau(X)然后看tau有沒有系統(tǒng)性的分化。# 因果森林估條件平均處理效應看異質性 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, # 誠實劃分避免用同一批數(shù)據(jù)既分裂又估值 min.node.size 10, # 葉子最小樣本太小會過擬合 sample.fraction 0.5, # 子采樣比例 seed 20240101 ) # 全人群平均效應target.sample overlap 只在有重疊的人群上平均 ate - grf::average_treatment_effect(cf, target.sample overlap) # 變量重要性看哪些協(xié)變量在驅動異質性 vimp - grf::variable_importance(cf) list(cf cf, ate ate, vimp vimp) }honesty TRUE是grf的默認值作用是樣本分裂和效應估計用不同的子樣本防止過擬合把噪聲當成異質性。這個參數(shù)不要關掉關了之后tau的分布會變得極其“好看”但全是假的。min.node.size 10是葉子節(jié)點的最小樣本量基因表達數(shù)據(jù)的噪聲大設太小會切出一堆只有幾個樣本的葉子估出來的tau方差爆炸。target.sample overlap這一項值得單獨說。默認是all在全員上平均改成overlap只在傾向評分有重疊的區(qū)域上平均更穩(wěn)健但也更保守。做暴露效應的基因篩查時用overlap結論不容易被少數(shù)極端個體帶跑。variable_importance的輸出是一個矩陣第一列是原始重要性第二列是歸一化到 0-1 的版本?？辞拔迕蛪颉绻挲g或某個技術協(xié)變量排第一說明異質性主要來自技術因素而非生物學這時候要警惕別急著往機制上解釋。在基因層面因果森林只在候選基因子集上跑不是全基因組。原因很實際2000 棵樹乘以幾百個基因單機跑一晚上都不一定完而且大多數(shù)基因的異質性檢驗并不顯著跑全量是浪費算力。4.4 從 p 值到 q 值結果表的組裝三種估計跑完每個基因會拿到三套est和se。組裝結果表要解決兩個問題一是統(tǒng)一效應方向基因表達上調還是下調二是多重檢驗。基因數(shù)量上千直接拿 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 作為主估計它雙重穩(wěn)健 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這一列是經(jīng)驗產(chǎn)物。IPW 和 AIPW 的方向如果相反說明兩個模型的假設沖突得很厲害可能是傾向評分模型設定有問題也可能是結局模型的線性假設在這個基因上不成立。這類基因即使 q 值很小也不該寫進結論最多放進補充材料標注方法間不一致等著用敏感性分析去解釋。主估計選 AIPW 而不是 IPW理由是雙重穩(wěn)健只要傾向評分和結局模型有一個對估計就一致比單靠一個模型的 IPW 更扛得住設定誤差。5. 避坑清單因果推斷流水線里最容易翻車的五處這套流程跑通不難難的是跑出來的數(shù)能不能信。下面五條是我在這類項目上反復踩到的每一條都按“現(xiàn)象 → 原因 → 解決”整理配上排查代碼。5.1 權重爆炸與重疊性缺失現(xiàn)象IPW 估出來的效應比 AIPW 大好幾倍標準誤也大得離譜ESS只剩原始樣本的 10% 不到。原因傾向評分模型太“自信”給某些樣本賦了接近 0 或 1 的概率。常見于某個協(xié)變量幾乎能完全預測處理狀態(tài)——比如某一批次的所有樣本都來自吸煙組模型在這個批次上學到“批次 吸煙”權重自然爆掉。另一個成因是分類變量的稀疏水平某一類只有幾個人邏輯回歸在該維度上給出極端的系數(shù)。解決先看重疊性診斷再決定是修剪樣本還是截尾權重。# 權重診斷分布、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) }打印極端權重樣本的協(xié)變量是定位問題最快的路徑。如果這十個樣本全在同一個批次那就是批次導致的如果全是高齡那是年齡分布尾部的問題。ESS低于 30% 時我一般先嘗試trim 0.01還不行就改用重疊權重WeightIt里method glm, estimand ATO重疊權重對極端傾向評分不敏感代價是估計目標變成“重疊人群的效應”解釋時要寫清楚。5.2 批次、共線性與“靜默”的樣本錯位現(xiàn)象腳本一聲不吭跑完結果表也有但效應估計整體偏小或者某個已知的強效應基因死活不顯著。原因三種可能。批次與暴露強相關去批次時把處理效應一起減掉了協(xié)變量共線導致glm系數(shù)不穩(wěn)定傾向評分本身就不準表達矩陣和表型的樣本順序沒對齊merge或cbind之后列和行的對應關系已經(jīng)亂了但 R 不會報錯。解決三個檢查必須都在預處理階段做完——批次卡方檢驗3.2 節(jié)、設計矩陣秩檢查3.3 節(jié)、樣本 ID 斷言。# 樣本對齊的最后一道保險放在所有合并操作之后 assert_alignment - function(expr, meta) { stopifnot( 表達矩陣列名缺失 !is.null(colnames(expr)), 表型樣本 ID 缺失 sample_id %in% names(meta), 樣本數(shù)不一致 ncol(expr) nrow(meta), 樣本順序錯位 identical(colnames(expr), as.character(meta$sample_id)) ) invisible(TRUE) }identical這一條是硬斷言不做容錯順序不一致直接拋錯。我寧可讓腳本在第一步就炸也不想在寫了兩周報告之后發(fā)現(xiàn)樣本對錯了。5.3 因果森林在小樣本上的“看起來很準”現(xiàn)象因果森林跑出來的tau分布特別漂亮亞群間差異明顯但換一個隨機種子結果就變了。原因honesty FALSE或者min.node.size設得太小樹在噪聲上過擬合?；虮磉_數(shù)據(jù)本身噪聲大樣本量幾百例時2000 棵樹的平均還壓不住方差tau的分布就是隨機的。解決三個參數(shù)一起調——保持honesty TRUEmin.node.size不低于 10sample.fraction取 0.4 到 0.6 之間然后用兩個不同種子各跑一遍看tau的相關系數(shù)。# 穩(wěn)定性檢驗換種子重跑看條件效應是否可復現(xiàn) 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 相關系數(shù) %.3f\n, r)) if (r 0.7) warning(條件效應不穩(wěn)定異質性結論不可靠) invisible(r) }相關系數(shù)低于 0.7 就別寫異質性結論了老老實實報平均效應。這個檢驗成本不高值得每個候選基因都過一遍。5.4 多重檢驗的兩種誤用現(xiàn)象q 值算出來了顯著基因上千個但和已知的吸煙相關基因對不上。原因要么用了 Bonferroni在基因表達場景過于嚴苛真正的信號被壓掉要么在篩選和驗證之間反復調閾值最后報的是“調出來的”顯著性。還有一種隱蔽的IPW 和 AIPW 各算一遍 q 值挑顯著的那個報這是變相的多重比較。解決主分析只跑一套AIPWq 值只算一次閾值預先定死比如 q 0.05敏感性分析的結果單獨列不參與篩選。方法學部分要寫清楚“主估計用的是哪一個”。5.5 用同一份數(shù)據(jù)調參又評估現(xiàn)象所有診斷都完美重疊性好、ESS 高、q 值分布漂亮但換一個獨立隊列就復現(xiàn)不出來。原因傾向評分模型的協(xié)變量集、trim閾值、因果森林的參數(shù)如果在看到結果之后反復調整等于把結局信息泄漏進了模型選擇。這類偏差不會體現(xiàn)在任何一張診斷圖上。解決把樣本按 7:3 隨機切開前 70% 用來定協(xié)變量集和超參數(shù)后 30% 只跑一次不看結果不回頭。如果樣本量不夠切就在分析計劃里把參數(shù)全部預先寫死中途改了要記錄在日志里。# 預注冊式的參數(shù)凍結把最終參數(shù)寫進文件后續(xù)腳本只讀不猜 freeze_params - function(cfg, path output/tables/frozen_params.rds) { saveRDS(cfg, path) message(參數(shù)已凍結至 , path, 后續(xù)分析請勿修改) }這條看著像流程潔癖但它是唯一能防住“調參調到顯著”的辦法。玄學的地方在于你自己調的時候完全不覺得在作弊每一步都有“合理”的理由。6. 二次開發(fā)把單基因流程改成批量篩選與可復現(xiàn)報告拿到這套腳本的人十有八九第一件事是想跑全基因組。單基因函數(shù)寫好了批量跑不是簡單套個for循環(huán)——直接循環(huán)會撞上三個問題內存、算力和失敗恢復。我給的做法是把估計函數(shù)包一層加上并行和緩存。先看批量層怎么寫# 批量估計并行 單基因失敗不影響整體 結果落盤緩存 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 的隨機數(shù)流可復現(xiàn)不加這一句因果森林那部分每次跑出來都不一樣。tryCatch把單基因的失敗隔離掉glm在完全分離的數(shù)據(jù)上會拋 warning 甚至 error一個基因炸掉不該拖垮整批。零方差判斷放在最前面是最常見的崩潰點——某些基因在所有樣本里表達值都一樣lm直接返回 NA 系數(shù)。n_workers不建議超過物理核心數(shù)減一grf內部本身也會占內存開太多反而會因為內存交換變慢。第二個改造點是報告層。原始腳本輸出的是 CSV二次開發(fā)可以把它變成參數(shù)化的 Rmd 或 Quarto 文檔把 q 值分布、火山圖、前若干個基因的森林圖一次性生成換一批數(shù)據(jù)只改cfg里的路徑。這里的關鍵是把“結果表”和“圖”解耦圖只讀表不重算估計。第三個改造點是加一層敏感性分析。做法很簡單把trim從 0.01 換成 0.05 和 0.1各跑一遍 AIPW看顯著基因的重疊比例。敏感性維度參數(shù)觀察指標警戒線權重截尾trim∈ {0.01, 0.05, 0.1}顯著基因重疊率低于 70% 說明結論依賴閾值結局模型線性 vs 回歸森林效應估計相關系數(shù)低于 0.8 說明模型設定敏感協(xié)變量集加/減批次變量q 0.05 基因數(shù)變化數(shù)量差一倍以上要解釋這張表建議直接貼進補充材料比任何文字描述都有說服力。審稿人問“你的結論穩(wěn)健嗎”把重疊率甩過去就行。最后一個技巧是日志。把sessionInfo()、cfg、每個腳本的運行時間和輸入文件的行列數(shù)統(tǒng)一寫進output/logs/一個腳本一個文件。這個習慣我是被坑出來的有次報告里的數(shù)和三個月前跑的對不上翻遍代碼也沒改過最后發(fā)現(xiàn)是上游給的表達矩陣悄悄更新了一版樣本數(shù)從 512 變成了 508。從那以后我每次跑因果估計流水線都強制在第一步把sessionInfo()和輸入文件的 md5 一起寫進日志跑完再核對一遍才敢往報告里放數(shù)。這類工程的可靠性不在估計量寫得多漂亮而在每一步都能回溯。希望幫到你。本文還有配套的精品資源點擊獲取