簡(jiǎn)介SmokingGeneExpression 面向生物信息學(xué)、統(tǒng)計(jì)分析與因果推斷學(xué)習(xí)者圍繞吸煙對(duì)健康肺組織基因表達(dá)的短期因果效應(yīng)提供一套可復(fù)現(xiàn)的 R 語(yǔ)言項(xiàng)目材料。壓縮包共 45 個(gè)文件約 3.1MB包含 19 個(gè) png 結(jié)果圖、13 個(gè) txt 數(shù)據(jù)/說明、8 個(gè) xlsx 表格、2 個(gè) csv、1 個(gè) pdf 及 R 腳本等覆蓋樣本 PCA、基因表達(dá)表、差異基因計(jì)數(shù)與輔助數(shù)據(jù)目錄按數(shù)據(jù)、結(jié)果和腳本組織便于按模塊檢索。項(xiàng)目涉及工具變量、傾向得分匹配、DAG 等因果識(shí)別思路并結(jié)合 RNA-seq 預(yù)處理、質(zhì)量控制、映射計(jì)數(shù)和差異表達(dá)分析幫助讀者理解從數(shù)據(jù)整理到結(jié)果解釋的完整流程。已有 148 人學(xué)習(xí)適合作為因果推斷與基因表達(dá)分析交叉練習(xí)的參考案例也可用于復(fù)現(xiàn)實(shí)驗(yàn)、拆解分析步驟和積累 R 語(yǔ)言生物信息學(xué)實(shí)踐并提升分析復(fù)現(xiàn)與結(jié)果解釋的可靠性。1. 相關(guān)性的天花板吸煙與肺組織基因表達(dá)之間隔著三重識(shí)別難題把吸煙者和非吸煙者的肺組織轉(zhuǎn)錄組并排跑一次差異表達(dá)幾乎總能得到一份長(zhǎng)長(zhǎng)的基因清單。這份資源要把問題往回退一步不問「哪些基因和吸煙同時(shí)出現(xiàn)」而問「同一個(gè)人的吸煙狀態(tài)被改變表達(dá)會(huì)往哪邊走」。因果推斷在這里只是換個(gè)名字說相關(guān)系數(shù)還不夠它得先證明自己配得上「因果」二字。R 在這條鏈路里也不只是畫圖——從傾向評(píng)分到多重檢驗(yàn)校正整條流程都落在同一套腳本生態(tài)里。值得拆開看是因?yàn)檫@類分析最容易在三個(gè)地方悄悄失效混雜沒調(diào)干凈、權(quán)重把樣本量吃光、把一個(gè)橫斷面的短期效應(yīng)講成了終身結(jié)論。2. 把混雜畫成圖DAG、后門準(zhǔn)則與「哪些變量碰不得」差異表達(dá)分析把「吸煙」當(dāng)成一個(gè)分組標(biāo)簽因果推斷把它當(dāng)成一次干預(yù)。這個(gè)視角切換的代價(jià)是你必須先寫下自己的因果假設(shè)而不能讓模型自己去猜。DAG 就是這份假設(shè)的書面形式哪些箭頭存在、哪些箭頭不存在一旦畫出來調(diào)什么、不調(diào)什么就不再是主觀選擇而是圖上的可計(jì)算問題。2.1 潛在結(jié)局與三個(gè)識(shí)別假設(shè)可交換性到底在賭什么潛在結(jié)局框架的寫法很樸素每個(gè)樣本有兩個(gè)潛在表達(dá)值Y1 表示「如果這個(gè)人吸煙」Y? 表示「如果這個(gè)人不吸煙」我們想看的是 E[Y1] ? E[Y?]。問題在于每個(gè)人只能落在其中一個(gè)世界里觀測(cè)到的永遠(yuǎn)只有 Y Y^X。所以估計(jì)量能成立靠的是三個(gè)假設(shè)撐住。第一個(gè)是一致性意思是「觀測(cè)到的暴露水平」和「潛在結(jié)局上標(biāo)記的暴露水平」必須對(duì)得上。這條在轉(zhuǎn)錄組數(shù)據(jù)里比在臨床試驗(yàn)里麻煩得多吸煙不是一個(gè)良定義的二元處理劑量、年限、最后一支煙到采樣的間隔都會(huì)改變暴露的實(shí)際含義。常見的做法是把暴露二值化成「采樣前 12 個(gè)月內(nèi)是否規(guī)律吸煙」并在方法里把這個(gè)窗口寫死。窗口越模糊一致性假設(shè)越容易被質(zhì)疑后面所有數(shù)字的可信度也跟著打折。第二個(gè)是條件可交換性也就是在給定一組協(xié)變量 Z 之后吸煙組和非吸煙組在潛在結(jié)局上沒有系統(tǒng)差異。這是整套方法里唯一無法用數(shù)據(jù)驗(yàn)證的假設(shè)只能用領(lǐng)域知識(shí)和敏感性分析去支撐。我一般會(huì)這么理解它不是「證明沒有混雜」而是「聲明我把哪些混雜處理掉了剩下的留給偏倚分析去量化」。第三個(gè)是正性即每個(gè)協(xié)變量組合下既有人吸煙也有人不吸煙。聽起來是技術(shù)細(xì)節(jié)實(shí)際是轉(zhuǎn)錄組因果分析里翻車最多的地方如果某一批樣本幾乎全是吸煙者那么這一批里根本不存在可比較的反事實(shí)傾向評(píng)分會(huì)逼近 0 或 1權(quán)重隨之發(fā)散。正性一旦違背再精巧的估計(jì)量也救不回來只能縮小目標(biāo)人群或做限制。還有一個(gè)常被忽略的前提是無干擾性即一個(gè)人的暴露狀態(tài)不會(huì)影響另一個(gè)人的表達(dá)。在個(gè)體層面的組織樣本里這一條基本成立但在同一批處理、同一張芯片上的樣本之間技術(shù)層面的「干擾」是真實(shí)存在的這也是批次必須進(jìn)圖的原因。2.2 后門路徑的枚舉從 DAG 到最小充分調(diào)整集把假設(shè)畫成圖之后「該調(diào)什么」就變成一個(gè)圖論問題。后門準(zhǔn)則說的是調(diào)整集 Z 必須阻斷從暴露到結(jié)局的所有后門路徑同時(shí)不能包含暴露的后代。前者防止漏調(diào)后者防止過度控制。這兩個(gè)方向的錯(cuò)誤在轉(zhuǎn)錄組分析里都極其常見。下面是一張最小可用的肺組織轉(zhuǎn)錄組 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 }) # 枚舉所有滿足后門準(zhǔn)則的最小調(diào)整集 adjustmentSets(g, exposure smoke, outcome expr, effect total)這段代碼的價(jià)值不在輸出本身而在它逼你把「為什么調(diào)這個(gè)變量」說清楚。effect total表示要估計(jì)總效應(yīng)此時(shí)inflam炎癥通路活性和cellcomp細(xì)胞組成作為中介變量不進(jìn)調(diào)整集因?yàn)樗鼈兲幵趕moke - ... - expr的路徑上。如果你把effect改成directdagitty 會(huì)把中介加進(jìn)調(diào)整集——但那時(shí)候你得到的就不再是「吸煙對(duì)表達(dá)的總影響」而是「扣掉炎癥與細(xì)胞組成之后剩下的直接影響」兩者是完全不同的問題報(bào)告時(shí)必須寫明白。purity組織純度和batch是兩個(gè)容易爭(zhēng)議的節(jié)點(diǎn)。組織純度既是樣本處理的產(chǎn)物也會(huì)影響可測(cè)量的表達(dá)譜同時(shí)在某些隊(duì)列里與吸煙狀態(tài)相關(guān)進(jìn)調(diào)整集是合理的。批次同理但如果批次與暴露高度共線調(diào)整會(huì)撞上正性違背這時(shí)候不是「調(diào)不調(diào)」的問題而是這批數(shù)據(jù)不該被用來回答這個(gè)問題。我通常會(huì)把a(bǔ)djustmentSets()的結(jié)果和領(lǐng)域知識(shí)對(duì)一遍而不是照單全收。圖是我畫的圖錯(cuò)了后門準(zhǔn)則只會(huì)幫我把錯(cuò)誤執(zhí)行得更徹底。2.3 中介與對(duì)撞兩類「看起來該調(diào)、其實(shí)不能調(diào)」的變量中介變量是第一個(gè)陷阱。吸煙會(huì)誘導(dǎo)氣道炎癥炎癥通路活化又會(huì)改變一大批基因的表達(dá)。如果你把炎癥相關(guān)的通路活性放進(jìn)調(diào)整集問的就不再是「吸煙會(huì)不會(huì)改變表達(dá)」而是「在不改變炎癥的前提下吸煙會(huì)不會(huì)改變表達(dá)」——這個(gè)反事實(shí)在現(xiàn)實(shí)里幾乎不存在。過度控制偏倚的機(jī)制就是這樣你從總效應(yīng)里切掉了一段真實(shí)存在的通路剩下的估計(jì)量既不是總效應(yīng)也不是干凈的直接效應(yīng)而是一個(gè)沒有明確解釋的混合體。這里有一個(gè)判斷標(biāo)準(zhǔn)這個(gè)變量是不是在暴露發(fā)生之后才被改變的如果是它大概率是中介或者中介的下游默認(rèn)不該進(jìn)總效應(yīng)的調(diào)整集。要做中介分析就得換一套更強(qiáng)的識(shí)別假設(shè)而不是簡(jiǎn)單地在回歸里多加一列。對(duì)撞變量是第二個(gè)陷阱也是更隱蔽的一個(gè)。對(duì)撞的意思是某個(gè)變量同時(shí)被暴露和結(jié)局或結(jié)局的原因影響。按對(duì)撞分層會(huì)在暴露與結(jié)局之間憑空造出一條關(guān)聯(lián)。舉個(gè)具體的場(chǎng)景假設(shè)是否進(jìn)入某項(xiàng)后續(xù)檢查同時(shí)受吸煙狀態(tài)和某些基因表達(dá)水平影響而你為了「控制數(shù)據(jù)質(zhì)量」只保留了做過這項(xiàng)檢查的樣本。這一步篩選就把對(duì)撞條件化了剩下的樣本里吸煙與表達(dá)之間會(huì)出現(xiàn)一個(gè)在真實(shí)人群中不存在的關(guān)聯(lián)方向還可能完全反掉。這類偏倚在生信流程里極難察覺因?yàn)楹Y選步驟通常被寫在「數(shù)據(jù)清洗」里沒人會(huì)把它和因果結(jié)構(gòu)聯(lián)系起來。還有一類變量介于兩者之間細(xì)胞類型組成。吸煙會(huì)改變肺組織的細(xì)胞構(gòu)成比如免疫細(xì)胞浸潤(rùn)比例上升而不同細(xì)胞類型的表達(dá)譜差異巨大。所以細(xì)胞組成既可能是混雜如果它影響采樣與暴露的關(guān)聯(lián)也可能是中介如果它確實(shí)是吸煙導(dǎo)致的下游改變。這個(gè)岔路口值得單獨(dú)拿出來講后面第 5 章會(huì)展開。3. 估計(jì)量的取舍傾向評(píng)分與 g-computation 在表達(dá)矩陣上怎么落地有了調(diào)整集剩下的問題是用什么估計(jì)量把它落實(shí)。轉(zhuǎn)錄組數(shù)據(jù)有幾個(gè)很實(shí)際的特征樣本量通常在幾百這個(gè)量級(jí)協(xié)變量維度不高但彼此相關(guān)結(jié)局維度是兩萬(wàn)左右。這三條決定了估計(jì)量不能照搬流行病學(xué)教科書里的默認(rèn)選項(xiàng)。3.1 傾向評(píng)分的直覺每個(gè)樣本代表多少人傾向評(píng)分 e(Z) P(X 1 | Z) 的定義很直白在給定協(xié)變量的情況下這個(gè)人吸煙的概率。它的作用是把一堆協(xié)變量壓縮成一個(gè)標(biāo)量讓「調(diào)整」變成「加權(quán)」。直觀理解是逆概率加權(quán)給每個(gè)樣本發(fā)一張票票的權(quán)重表示「這個(gè)人在目標(biāo)人群里代表多少人」。一個(gè)吸煙概率只有 0.1 卻實(shí)際吸煙的人代表了 10 個(gè)像他這樣的人所以權(quán)重要放大到 10 倍。這種加權(quán)在協(xié)變量平衡上很有效但代價(jià)是方差。權(quán)重分布的尾部決定了估計(jì)量的穩(wěn)定性所以任何一份嚴(yán)肅的 IPW 分析都應(yīng)該報(bào)告有效樣本量ESS (Σw)2 / Σw2ESS 掉到原始樣本量的一半以下時(shí)估計(jì)量的方差會(huì)明顯膨脹。我在實(shí)際操作里會(huì)把 ESS 當(dāng)成一個(gè)硬性診斷指標(biāo)而不是等結(jié)果跑完再去解釋為什么置信區(qū)間那么寬。傾向評(píng)分模型本身也有選擇logistic 回歸簡(jiǎn)單、可解釋、系數(shù)能看但對(duì)非線性關(guān)系不敏感梯度提升或超級(jí)學(xué)習(xí)器擬合能力更強(qiáng)但容易過擬合在小樣本上反而不穩(wěn)。常見的折中是先用 logistic 回歸加必要的樣條項(xiàng)再用平衡診斷標(biāo)準(zhǔn)化均值差驗(yàn)證加權(quán)后的協(xié)變量是否真的對(duì)齊而不是看模型 AUC 有多高。AUC 高不代表平衡好這一點(diǎn)很多人會(huì)搞混。3.2 為什么兩萬(wàn)個(gè)基因只需要一套權(quán)重這是這套分析里最值得講清楚的一個(gè)設(shè)計(jì)點(diǎn)也是計(jì)算成本能壓下來的關(guān)鍵。傾向評(píng)分模型是 P(X | Z)它只涉及暴露和協(xié)變量完全不碰任何基因的表達(dá)值。換句話說無論你分析的是第一個(gè)基因還是第一萬(wàn)個(gè)基因暴露模型都是同一個(gè)。所以最優(yōu)做法是擬合一次傾向評(píng)分得到一組權(quán)重然后在整個(gè)表達(dá)矩陣上復(fù)用。結(jié)局模型則相反E[Y | X, Z] 里 Y 是具體的基因每個(gè)基因一個(gè)模型必須逐個(gè)擬合。這個(gè)不對(duì)稱性直接決定了流水線的結(jié)構(gòu)環(huán)節(jié)依賴結(jié)局?jǐn)M合次數(shù)可否復(fù)用傾向評(píng)分模型否1全矩陣復(fù)用逆概率加權(quán)估計(jì)否0純加權(quán)計(jì)算全矩陣復(fù)用g-computation 結(jié)局模型是2 × 基因數(shù)不可復(fù)用AIPW 校正項(xiàng)是2 × 基因數(shù)不可復(fù)用理解這張表之后很多工程上的取舍就自然了如果樣本量很大、基因很多純 IPW 的加權(quán)均值差計(jì)算可以完全向量化成本幾乎可以忽略而 AIPW 因?yàn)橐o每個(gè)基因擬合兩個(gè)結(jié)局模型計(jì)算量會(huì)隨基因數(shù)線性增長(zhǎng)。這也是為什么在探索階段大家更愿意先跑 IPW 拿到一個(gè)粗篩結(jié)果再對(duì)顯著的那部分基因補(bǔ) AIPW 校正。順帶說一句這個(gè)不對(duì)稱性還有一個(gè)統(tǒng)計(jì)上的好處暴露模型只擬合一次意味著權(quán)重只估計(jì)一次不會(huì)因?yàn)榛虿煌腩~外的模型選擇變異性。用一個(gè)固定的權(quán)重去掃全矩陣至少保證了所有基因是在同一個(gè)「?jìng)稳巳骸估锉容^的。3.3 代碼骨架從 PS 到 AIPW 的最小可運(yùn)行鏈路下面這段代碼把上面的思路落成一個(gè)可以跑通的鏈路。協(xié)變量表cov_df每行一個(gè)樣本表達(dá)矩陣expr_mat每行一個(gè)基因列順序與協(xié)變量表一致。# ---- 1. 暴露模型一次擬合全基因復(fù)用 ---- ps_fit - glm(smoke ~ age sex bmi purity batch, family binomial(link logit), data cov_df) cov_df$ps - predict(ps_fit, type response) # 截尾把極端傾向評(píng)分壓到 [0.01, 0.99]避免單個(gè)樣本權(quán)重爆炸 cov_df$ps_t - pmin(pmax(cov_df$ps, 0.01), 0.99) # ---- 2. 穩(wěn)定化權(quán)重 ---- # 分子用邊際暴露概率防止權(quán)重整體漂移導(dǎo)致方差虛高 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)) # 有效樣本量權(quán)重越不均勻ESS 越小方差越大 ess - sum(cov_df$w)^2 / sum(cov_df$w^2) cat(sprintf(n %d, ESS %.1f\n, nrow(cov_df), ess)) # ---- 3. 加權(quán)均值差估計(jì)單個(gè)基因的 ATE ---- ate_ipw - function(y, x, w) { # y 應(yīng)先做 log2(CPM 1) 轉(zhuǎn)換并做分位數(shù)標(biāo)準(zhǔn)化 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這一步需要解釋清楚截尾不是「修正」正性違背而是把極端權(quán)重壓下來?yè)Q取方差穩(wěn)定代價(jià)是引入一點(diǎn)點(diǎn)偏倚。閾值取 0.01/0.99 是常見起點(diǎn)實(shí)際應(yīng)該看傾向評(píng)分的分布形狀決定如果兩端本來就很稀疏壓到 0.05/0.95 也合理。ate_ipw是加權(quán)均值差本質(zhì)上是 Horvitz-Thompson 型估計(jì)量的穩(wěn)定化版本。它的假設(shè)是「加權(quán)后兩組可比」但沒有用結(jié)局模型做任何校正所以在協(xié)變量與結(jié)局關(guān)系復(fù)雜時(shí)會(huì)有殘余偏倚。補(bǔ)一個(gè) AIPW 版本可以拿到雙重穩(wěn)健性ate_aipw - function(y, d, ps) { # 結(jié)局模型分別在兩組內(nèi)擬合再對(duì)全體樣本做預(yù)測(cè) 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) # 雙重穩(wěn)健結(jié)局模型或傾向評(píng)分模型有一個(gè)對(duì)估計(jì)量就一致 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 的預(yù)測(cè)差第二段和第三段是殘差校正項(xiàng)用傾向評(píng)分的倒數(shù)給殘差加權(quán)。雙重穩(wěn)健的意思是f1/f0設(shè)定正確時(shí)即使傾向評(píng)分模型錯(cuò)了估計(jì)也一致反過來傾向評(píng)分對(duì)、結(jié)局模型錯(cuò)估計(jì)仍然一致。兩個(gè)都錯(cuò)才會(huì)失效這比單靠一種模型實(shí)在得多。代價(jià)也很明確每個(gè)基因要擬合兩個(gè)線性模型兩萬(wàn)個(gè)基因就是四萬(wàn)次擬合。在幾百個(gè)樣本的規(guī)模上R 里的lm會(huì)慢但能忍受如果基因數(shù)再往上走就得換limma的矩陣化擬合或者speedglm把循環(huán)壓掉。還有一點(diǎn)容易被忽略ps必須和w使用同一個(gè)截尾后的版本否則 AIPW 的校正項(xiàng)和加權(quán)項(xiàng)不在同一個(gè)尺度上估計(jì)量會(huì)偏。這類細(xì)節(jié)不會(huì)報(bào)錯(cuò)只會(huì)讓結(jié)果默默錯(cuò)掉。4. 兩萬(wàn)個(gè)基因同時(shí)檢驗(yàn)FDR、陰性對(duì)照與未測(cè)混雜的敏感性分析單個(gè)基因的因果效應(yīng)估計(jì)出來之后問題并沒有結(jié)束。多重檢驗(yàn)、殘余混雜、模型設(shè)定錯(cuò)誤這三件事會(huì)在規(guī)模效應(yīng)下同時(shí)放大。一份可信的分析必須在這三處都給出交代而不是只報(bào)一個(gè) q 值。4.1 BH 校正解決的是規(guī)模不是偏倚兩萬(wàn)個(gè)基因同時(shí)檢驗(yàn)即使全部為零效應(yīng)按 α 0.05 也會(huì)有一千個(gè)「顯著」。Benjamini-Hochberg 過程把假發(fā)現(xiàn)率的期望控制在設(shè)定水平這是標(biāo)準(zhǔn)做法res$p_value - 2 * pnorm(-abs(res$ate / res$se)) res$q_value - p.adjust(res$p_value, method BH) # 報(bào)告時(shí)同時(shí)給出效應(yīng)量與校正后顯著性 sig - subset(res, q_value 0.05 abs(ate) 0.2)這里有兩個(gè)容易被忽略的點(diǎn)。第一p 值來自漸近正態(tài)近似在 ESS 掉到三位數(shù)時(shí)這個(gè)近似很不靠譜置換檢驗(yàn)會(huì)更穩(wěn)set.seed(2024) perm_null - replicate(200, { # 打亂暴露標(biāo)簽保留協(xié)變量與權(quán)重的聯(lián)合結(jié)構(gòu) idx - sample(nrow(cov_df)) ate_ipw(expr_mat[1, ], cov_df$smoke[idx], cov_df$w) })第二也是更根本的BH 控制的是「零假設(shè)為真時(shí)的誤報(bào)比例」它假設(shè)你的 p 值本身是有效的。如果因果估計(jì)存在系統(tǒng)性偏倚校正后的 q 值只會(huì)讓偏倚變得更「顯著」。FDR 處理的是規(guī)模問題不是識(shí)別問題這兩件事必須分開看。4.2 陰性對(duì)照用「不該有信號(hào)」的地方校驗(yàn)流程陰性對(duì)照是這套流程里性價(jià)比最高的驗(yàn)證手段因?yàn)樗恍枰~外的數(shù)據(jù)只需要你事先想清楚「哪里本來就不該有信號(hào)」。陰性對(duì)照結(jié)局的做法是挑一組已知與吸煙沒有因果關(guān)系的基因比如某些管家基因或者與暴露完全無關(guān)的組織特異性基因跑完整條流水線。如果它們?cè)?q 0.05 上大量顯著說明有殘余混雜或技術(shù)偏倚在推動(dòng)整個(gè)分布這時(shí)候主結(jié)果里的信號(hào)也要打折扣。這個(gè)診斷的說服力很強(qiáng)因?yàn)樗鼨z驗(yàn)的不是「我的方法有沒有檢出真信號(hào)」而是「我的方法會(huì)不會(huì)制造假信號(hào)」。陰性對(duì)照暴露的思路類似找一個(gè)理論上不可能影響基因表達(dá)、但在數(shù)據(jù)里可能與表達(dá)相關(guān)的變量比如樣本采集的時(shí)間編號(hào)或者某個(gè)隨機(jī)分配的技術(shù)參數(shù)。如果它跑出了「顯著效應(yīng)」幾乎可以肯定流程里混進(jìn)了虛假關(guān)聯(lián)。我一般會(huì)把這組對(duì)照放在主分析之前跑。先確認(rèn)系統(tǒng)是干凈的再看主結(jié)果比事后解釋一堆異常要省力得多。4.3 E-value未測(cè)混雜需要多強(qiáng)才能掀翻結(jié)論條件可交換性無法用數(shù)據(jù)驗(yàn)證但可以量化它有多脆弱。E-value 回答的問題是需要多強(qiáng)的未測(cè)混雜與暴露和結(jié)局的關(guān)聯(lián)強(qiáng)度都達(dá)到這個(gè)水平才能把觀測(cè)到的效應(yīng)完全解釋掉。E-value 越大結(jié)論越穩(wěn)。對(duì)連續(xù)結(jié)局需要先把標(biāo)準(zhǔn)化均值差 d 轉(zhuǎn)成近似的風(fēng)險(xiǎn)比evalue_from_d - function(d) { rr - exp(0.91 * abs(d)) # 標(biāo)準(zhǔn)化均值差到風(fēng)險(xiǎn)比的常用近似 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 只是一個(gè)「脆弱性指標(biāo)」不是「混雜存在的概率」。一個(gè) E-value 為 1.5 的結(jié)論意味著只需要一個(gè)中等強(qiáng)度的未測(cè)混雜就能翻盤這種結(jié)果在報(bào)告里應(yīng)該被明確標(biāo)成「不穩(wěn)健」而不是和 E-value 為 4 的結(jié)果混在一起講。0.91這個(gè)系數(shù)來自標(biāo)準(zhǔn)化均值差與風(fēng)險(xiǎn)比的近似轉(zhuǎn)換量綱上沒問題但前提是結(jié)局近似正態(tài)表達(dá)值做過 log 轉(zhuǎn)換之后這個(gè)前提大體成立。把 FDR、陰性對(duì)照、E-value 三者放在一起看你會(huì)得到一張比「顯著基因列表」信息量大得多的表哪些基因既有統(tǒng)計(jì)顯著性又有足夠大的效應(yīng)量哪些基因的結(jié)論經(jīng)得起中等強(qiáng)度混雜的沖擊。這份表才是這份資源真正值得深挖的部分。5. 排查清單五個(gè)把因果結(jié)論做廢的典型現(xiàn)場(chǎng)下面這五條都是我在實(shí)際流程里反復(fù)遇到的。每一條都按「現(xiàn)象 → 原因 → 解決」寫方便直接對(duì)照自己跑出來的結(jié)果。5.1 權(quán)重極端與 ESS 崩塌現(xiàn)象權(quán)重最大值到了幾十某個(gè)基因的 ATE 在去掉一兩個(gè)樣本后符號(hào)翻轉(zhuǎn)置信區(qū)間寬得沒有意義。原因傾向評(píng)分的尾部重疊很差。常見誘因是把「吸煙包年數(shù)」這類直接定義暴露的變量塞進(jìn)暴露模型或者協(xié)變量里有一個(gè)與暴露幾乎一一對(duì)應(yīng)的變量比如某個(gè)只在吸煙者中采集的樣本批次。這兩種情況都會(huì)把 ps 推到 0 或 1 附近。解決先畫兩組傾向評(píng)分的疊加密度圖看重疊區(qū)域。然后做兩件事——把權(quán)重截尾到 [0.01, 0.99]報(bào)告截尾前后的 ESS 變化如果 ESS 仍然低于樣本量的一半改用重疊權(quán)重overlap weights它把目標(biāo)人群收縮到「兩組都可能出現(xiàn)」的區(qū)域代價(jià)是改變了估計(jì)目標(biāo)需要在方法里寫明。5.2 細(xì)胞組成的岔路口現(xiàn)象調(diào)整細(xì)胞類型比例之后一批免疫相關(guān)基因的效應(yīng)幾乎消失而調(diào)整之前它們是全矩陣?yán)镒铒@著的。原因bulk 組織的表達(dá)值是「細(xì)胞組成」和「細(xì)胞內(nèi)表達(dá)」的混合。吸煙會(huì)改變肺組織的細(xì)胞構(gòu)成所以組成變化既是真實(shí)的效應(yīng)通路也是觀測(cè)表達(dá)的驅(qū)動(dòng)因素。把它當(dāng)混雜調(diào)掉等于把一部分真實(shí)效應(yīng)切走不調(diào)又會(huì)把組成差異誤讀成細(xì)胞內(nèi)的調(diào)控變化。解決先做去卷積拿到每個(gè)樣本的細(xì)胞類型比例然后分兩條線報(bào)告一條是總效應(yīng)不調(diào)組成一條是調(diào)整組成后的效應(yīng)近似「細(xì)胞內(nèi)表達(dá)變化」。兩條線的差異本身就是信息它告訴你這個(gè)基因的效應(yīng)有多少來自組織構(gòu)成重塑。只報(bào)其中一條都是在隱藏關(guān)鍵結(jié)論。5.3 批次與暴露共線現(xiàn)象傾向評(píng)分模型里批次變量的系數(shù)很大加入批次后 ESS 從 400 掉到 150一批基因的效應(yīng)方向整體反轉(zhuǎn)。原因吸煙樣本集中在某幾個(gè)批次里。這可能是樣本收集年份不同、也可能是樣本來源不同。此時(shí)批次和暴露在數(shù)據(jù)里幾乎不可分正性假設(shè)在批次層面被違背。解決先做 batch × smoke 的列聯(lián)表看每個(gè)批次里兩組的樣本數(shù)。如果某些批次只有一組調(diào)整解決不了問題只能限制到兩組都有的批次或者把批次作為分層變量做分層估計(jì)再合并。硬調(diào)一個(gè)共線的協(xié)變量得到的不是「控制」是把暴露的變異也一起控掉了。5.4 選擇性推斷先篩后估現(xiàn)象先在全部基因上做關(guān)聯(lián)篩選再對(duì)篩出來的基因做因果估計(jì)q 值一片漂亮。原因同一份數(shù)據(jù)被用了兩次——一次用來選基因一次用來做推斷。篩選步驟本身消耗了數(shù)據(jù)里的信息后續(xù)的 p 值不再服從零假設(shè)下的均勻分布FDR 控制失效。解決樣本分裂是最直接的辦法一半選基因、一半估效應(yīng)。樣本量不夠時(shí)至少做一次獨(dú)立性檢查把數(shù)據(jù)隨機(jī)分成兩半各自跑完整的篩選加估計(jì)看兩次的顯著基因重合度。重合度低就說明結(jié)果高度依賴具體的篩選路徑。5.5 把「FDR 顯著」當(dāng)成機(jī)制證據(jù)現(xiàn)象報(bào)告里寫「吸煙通過某炎癥通路導(dǎo)致這批基因表達(dá)改變」。原因主效應(yīng)顯著只說明暴露和結(jié)局之間有關(guān)聯(lián)完全不涉及中介。要聲稱「通過某通路」需要正式的中介分析而中介分析要求跨世界反事實(shí)假設(shè)比如「吸煙但通路活性被固定在不吸煙水平」這種情況是否存在橫斷面數(shù)據(jù)幾乎無法滿足。解決把結(jié)論降級(jí)。能說的是「這批基因的表達(dá)變化與炎癥通路活性變化在方向上一致」不能說「由炎癥通路介導(dǎo)」。真要做中介就得換數(shù)據(jù)設(shè)計(jì)或者至少用敏感性分析把中介假設(shè)的脆弱程度量化出來。6. 短期因果的邊界這份資源能回答什么、不該被拿來回答什么寫到這里有一個(gè)判斷需要擺到明面上這套流程估計(jì)的是「短期」因果而這個(gè)「短期」不是修辭是數(shù)據(jù)結(jié)構(gòu)決定的。暴露是采樣時(shí)點(diǎn)的近期吸煙狀態(tài)結(jié)局是同一時(shí)點(diǎn)的轉(zhuǎn)錄狀態(tài)兩者之間沒有時(shí)間維度。所以它能回答的是「當(dāng)前吸煙狀態(tài)與當(dāng)前表達(dá)之間的差異在多大程度上不能被已測(cè)混雜解釋掉」而不是「吸煙若干年后表達(dá)會(huì)怎么變」。把能答和不能答列清楚比任何模型細(xì)節(jié)都重要問題這套結(jié)構(gòu)能否回答卡在哪當(dāng)前吸煙與當(dāng)前表達(dá)的差異是否被已測(cè)混雜解釋掉能在三個(gè)識(shí)別假設(shè)成立時(shí)假設(shè)本身無法用數(shù)據(jù)驗(yàn)證只能靠敏感性分析戒煙若干年后表達(dá)能否回到非吸煙水平不能橫斷面沒有隨訪也沒有戒煙時(shí)長(zhǎng)表達(dá)變化有多少來自細(xì)胞組成改變部分能需要去卷積且組成變量的中介/混雜身份需逐基因判斷某基因的效應(yīng)是否由某通路介導(dǎo)不能直接答中介分析需要更強(qiáng)的跨世界假設(shè)效應(yīng)在某個(gè)疾病亞組里是否不同通常不能分層后每格樣本量太小交互檢驗(yàn)功效不足這張表也解釋了為什么第 4 章的敏感性分析值得花那么多篇幅。橫斷面數(shù)據(jù)里條件可交換性永遠(yuǎn)是一個(gè)聲明而不是一個(gè)事實(shí)你能做的只是把「這個(gè)聲明有多脆弱」量化出來。E-value 和陰性對(duì)照不能證明因果但它們能在結(jié)論站不住的時(shí)候提前報(bào)警。還有一個(gè)更微妙的取舍這套流程把「總效應(yīng)」作為默認(rèn)目標(biāo)也就是不調(diào)整中介變量。這個(gè)選擇的代價(jià)是估計(jì)量里混進(jìn)了細(xì)胞組成變化、炎癥通路活化等下游過程效應(yīng)的生物學(xué)解釋會(huì)變得模糊好處是它對(duì)應(yīng)一個(gè)真實(shí)存在、可被干預(yù)的問題——「讓人不吸煙表達(dá)譜會(huì)怎樣變化」。如果改成估計(jì)直接效應(yīng)數(shù)字會(huì)更「干凈」但那個(gè)反事實(shí)在現(xiàn)實(shí)里幾乎無法實(shí)現(xiàn)。我傾向前者因?yàn)樗辽賹?duì)應(yīng)一個(gè)能落地的干預(yù)。從那以后我每次跑這類腳本都會(huì)先強(qiáng)制走一遍「傾向評(píng)分重疊診斷 ESS 陰性對(duì)照」三件套再去看主結(jié)果。順序反過來人很容易被漂亮的 q 值帶著走等發(fā)現(xiàn)問題時(shí)已經(jīng)寫進(jìn)結(jié)論了。這份資源的價(jià)值也在這里它把這條容易走偏的鏈路固化成了一套可復(fù)查的步驟而不是一組可以直接抄的數(shù)字。想動(dòng)手的人從傾向評(píng)分那一段開始對(duì)照自己的數(shù)據(jù)改會(huì)比從頭搭一遍省很多時(shí)間。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取