简介这份资源是面向生物信息学与统计分析学习者的R语言实战项目聚焦用因果推断方法评估吸烟对健康肺组织基因表达的短期因果效应适合具备一定R基础、想了解观察性数据因果分析流程的研究者与高年级学生。压缩包共45个文件约3.1MB包含19个png结果图、13个txt文本、8个xlsx表格、2个csv数据、1个pdf文档及R脚本等覆盖数据表、分析结果与可视化输出便于对照复现。项目涉及工具变量、倾向得分匹配、DAG等因果推断思路并结合RNA-seq数据的质控、比对、计数与差异表达分析以及通路富集和结果验证环节能帮助读者理解从原始表达数据到因果效应量化的完整链路。目前已有148人学习可作为因果推断与基因表达分析结合的小型案例参考。1. 吸烟与基因表达因果推断能回答相关分析答不了的问题做转录组分析的人多半遇到过这种尴尬吸烟样本和健康肺组织样本的差异表达基因列表拉出来几百个富集分析也漂亮但审稿人一句“这是相关还是因果”就把人问住了。常规差异表达分析DEG只能告诉你“吸烟者肺组织里基因A表达量低”但没法排除反向因果——比如某个基因低表达的人本身更容易吸烟或者存在一个共同混杂因素年龄、炎症基线、组织细胞比例同时影响吸烟行为和基因表达。SmokingGeneExpression 这个方向要做的就是用因果推断框架去估计吸烟对健康肺组织中基因表达的短期因果效应而不是停留在相关性层面。它适合三类人手里有吸烟/非吸烟人群肺组织转录组数据但被审稿人追问因果性的分析者想入门孟德尔随机化MR或工具变量法但缺少真实场景练手的生信人以及需要为公共卫生政策提供更硬证据的研究者。核心思路是把遗传变异作为工具变量或者用倾向性评分、断点回归等准实验设计把“吸烟”这个暴露变量从混杂网络中剥离出来再去看它对每个基因表达量的因果效应量。2. 因果推断选型为什么工具变量法在肺组织转录组里更稳2.1 从相关到因果三个必须跨过的门槛在肺组织转录组场景里做因果推断要跨过三个门槛。第一是混杂控制吸烟者往往年龄偏大、炎症水平更高、气道细胞比例不同这些都会直接影响基因表达如果不控制估计出的效应就是偏的。第二是反向因果某些基因的表达状态可能通过神经递质通路影响尼古丁依赖倾向导致“基因表达低→更容易吸烟”而不是反过来。第三是测量误差吸烟量通常靠问卷自报存在回忆偏倚而基因表达又有技术噪声两者叠加会让朴素回归的系数向零衰减。常规做法是用多元回归把年龄、性别、BMI、细胞类型比例作为协变量塞进去但这只能控制已测混杂未测混杂比如遗传背景、早年环境暴露依然会污染估计。更稳的做法是引入工具变量找一个与吸烟行为强相关、但与肺组织基因表达只能通过吸烟间接相关的遗传位点通常是 rs 开头的 SNP。这就是孟德尔随机化MR的核心逻辑——遗传变异在受精时随机分配相当于自然界的随机对照试验能绕过大部分混杂和反向因果。2.2 工具变量法的三个核心假设与检验方法工具变量法不是随便找个 SNP 就能用它依赖三个假设每个假设都有对应的检验手段。假设含义检验方法不满足的后果相关性工具变量与吸烟行为强相关第一阶段 F 统计量 10弱工具偏倚估计量向朴素回归靠拢独立性工具变量与混杂因素无关检查与年龄、性别、细胞比例的关联混杂重新引入因果解释失效排他性工具变量只通过吸烟影响基因表达多效性检验MR-Egger 截距、MR-PRESSO高估或低估因果效应第一阶段 F 统计量是最先要看的。实操中我一般用TwoSampleMR包跑F 值低于 10 的 SNP 直接剔除不要心疼。排他性假设最难验证通常用 MR-Egger 回归的截距项做粗略判断截距显著偏离零说明存在方向性多效性这时候要么换工具变量子集要么改用加权中位数法。2.3 短期因果效应的定义与时间窗口选择标题里“短期因果关系”是个容易被忽略的限定词。肺组织基因表达是动态的吸烟后几小时到几天的早期响应基因比如 CYP1A1、AHRR 通路和吸烟几十年的累积损伤基因比如与 COPD 相关的 MMP 家族因果结构完全不同。短期因果效应通常指暴露后 24 小时到 4 周内的表达变化对应的是急性炎症、氧化应激和代谢激活通路。在数据层面如果你用的是横断面人群样本其实很难严格区分短期和长期。常见做法是优先选择有明确吸烟状态记录当前吸烟/已戒烟/从不吸烟且戒烟时间在 4 周以内的样本作为“短期暴露组”或者用干预实验比如体外气液界面暴露模型的数据做验证。如果只有横断面数据就在论文里明确写“本研究估计的是当前吸烟状态对基因表达的因果效应时间尺度受限于横断面设计”不要硬说成短期。3. 从原始数据到因果估计一套可复现的 MR 分析流程3.1 数据准备暴露、结局与协变量三张表假设你手里有一份肺组织转录组表达矩阵行是基因列是样本和一份样本元数据吸烟状态、年龄、性别、细胞比例。第一步是把数据整理成 MR 分析需要的三张表暴露表吸烟状态和工具变量 SNP、结局表每个基因的表达量、协变量表年龄、性别、测序批次、细胞比例。import pandas as pd import numpy as np # 读取表达矩阵和元数据 expr pd.read_csv(lung_expr_matrix.csv, index_col0) # 行基因列样本 meta pd.read_csv(sample_metadata.csv, index_col0) # 行样本 # 对齐样本顺序这一步不做后面全错 common_samples expr.columns.intersection(meta.index) expr expr[common_samples] meta meta.loc[common_samples] # 构建暴露表吸烟状态编码为 0/1 exposure meta[[smoking_status]].copy() exposure[smoking_status] exposure[smoking_status].map( {never: 0, former: 0, current: 1} ) # 协变量表年龄、性别、批次、细胞比例 covariates meta[[age, sex, batch, cell_prop_epithelial, cell_prop_immune]] # 结局表选一个目标基因比如 CYP1A1 target_gene CYP1A1 outcome expr.loc[target_gene].to_frame(expression) print(f样本数: {len(common_samples)}, 暴露分布:\n{exposure[smoking_status].value_counts()})这段代码的关键在样本对齐。表达矩阵的列名和元数据的行名如果不一致后续所有合并都会错位而且 pandas 不会报错只会静默产生 NaN。我一般会在合并后打印样本数做二次确认。暴露编码里把“已戒烟”归到 0是因为短期因果效应关注当前暴露已戒烟者的基因表达可能已经部分恢复混在一起会稀释效应。3.2 工具变量筛选从 GWAS 汇总数据到独立 SNP 集工具变量通常来自公开的吸烟行为 GWAS 汇总数据。筛选标准一般设三条与吸烟行为关联的 p 值小于 5e-8、SNP 之间连锁不平衡 r² 小于 0.001、距离至少 10000 bp。下面是一个筛选脚本的骨架。# 假设 gwas_df 包含列: snp, beta, se, pval, chr, pos gwas_df pd.read_csv(smoking_gwas_summary.csv) # 第一步全基因组显著性筛选 sig_snps gwas_df[gwas_df[pval] 5e-8].copy() # 第二步按染色体和位置排序做 LD 剪枝简化版实际用 plink sig_snps sig_snps.sort_values([chr, pos]) pruned [] last_pos {} for _, row in sig_snps.iterrows(): key row[chr] if key not in last_pos or (row[pos] - last_pos[key]) 10000: pruned.append(row) last_pos[key] row[pos] iv_df pd.DataFrame(pruned) print(f独立工具变量数量: {len(iv_df)}) # 第三步计算第一阶段 F 统计量 iv_df[F_stat] (iv_df[beta] / iv_df[se]) ** 2 iv_df iv_df[iv_df[F_stat] 10] print(f通过 F10 的工具变量: {len(iv_df)})LD 剪枝这里用的是简化版距离法实际分析中应该用 PLINK 的--indep-pairwise 1000 10 0.001做标准剪枝。F 统计量公式是 (beta/se)²等价于 t 统计量的平方。F 大于 10 是经验阈值低于这个值说明工具变量解释的暴露变异太少估计结果不可信。注意这里 beta 和 se 都来自 GWAS 汇总数据不是你自己样本算的所以要求 GWAS 和你的转录组样本来自同一祖先群体否则等位基因频率不匹配会导致工具变量失效。3.3 因果效应估计Wald 比值法与 MR-Egger 回归拿到工具变量后最基础的因果效应估计是 Wald 比值法用工具变量对结局的效应除以工具变量对暴露的效应。单个 SNP 的 Wald 估计量方差可以用 delta 方法近似。多个 SNP 时用逆方差加权IVW合并。# 假设 iv_df 有 beta_exposure, se_exposure # outcome_df 有每个 SNP 对基因表达的 beta_outcome, se_outcome merged iv_df.merge(outcome_df, onsnp) # Wald 比值法 merged[wald_ratio] merged[beta_outcome] / merged[beta_exposure] merged[wald_se] np.sqrt( (merged[se_outcome] ** 2) / (merged[beta_exposure] ** 2) (merged[beta_outcome] ** 2) * (merged[se_exposure] ** 2) / (merged[beta_exposure] ** 4) ) # IVW 合并 weights 1 / merged[wald_se] ** 2 ivw_estimate np.sum(weights * merged[wald_ratio]) / np.sum(weights) ivw_se np.sqrt(1 / np.sum(weights)) print(fIVW 因果效应估计: {ivw_estimate:.4f} (SE: {ivw_se:.4f})) print(f95% CI: [{ivw_estimate - 1.96*ivw_se:.4f}, {ivw_estimate 1.96*ivw_se:.4f}])IVW 的假设是所有工具变量都是有效的如果存在多效性IVW 会被带偏。这时候要跑 MR-Egger 回归做敏感性分析Egger 回归的截距项如果显著不为零说明存在方向性多效性IVW 结果不可信。另一个常用方法是加权中位数法它在部分工具变量无效时仍然稳健但精度略低。实操中我一般三个方法都跑如果 IVW 和加权中位数方向一致且 Egger 截距不显著才认为因果证据比较硬。3.4 多基因并行把单基因流程封装成批量函数肺组织转录组有两万多个基因不可能一个一个手动跑。把上面的流程封装成函数对每个基因循环执行同时记录效应量、标准误、p 值和多效性检验结果。def mr_analysis_for_gene(gene_name, expr_df, iv_df, covariates_df): 对单个基因跑 MR 分析返回因果效应估计和诊断指标 outcome expr_df.loc[gene_name].to_frame(expression) merged iv_df.merge(outcome, left_onsnp, right_indexTrue, howinner) if len(merged) 3: return None # 工具变量太少跳过 # IVW merged[wald_ratio] merged[beta_outcome] / merged[beta_exposure] merged[wald_se] np.sqrt( (merged[se_outcome] ** 2) / (merged[beta_exposure] ** 2) (merged[beta_outcome] ** 2) * (merged[se_exposure] ** 2) / (merged[beta_exposure] ** 4) ) weights 1 / merged[wald_se] ** 2 ivw_est np.sum(weights * merged[wald_ratio]) / np.sum(weights) ivw_se np.sqrt(1 / np.sum(weights)) # MR-Egger 截距 from scipy import stats slope, intercept, r, p, se stats.linregress( merged[beta_exposure], merged[beta_outcome] ) return { gene: gene_name, ivw_beta: ivw_est, ivw_se: ivw_se, ivw_pval: 2 * stats.norm.sf(abs(ivw_est / ivw_se)), egger_intercept: intercept, egger_pval: p, n_iv: len(merged) } # 批量跑前 100 个高表达基因做测试 high_expr_genes expr.mean(axis1).sort_values(ascendingFalse).head(100).index results [mr_analysis_for_gene(g, expr, iv_df, covariates) for g in high_expr_genes] results_df pd.DataFrame([r for r in results if r is not None]) print(results_df.sort_values(ivw_pval).head(10))这个批量函数里有两个细节值得注意。一是howinner合并确保只有同时出现在工具变量表和表达矩阵里的 SNP 才参与分析避免 NaN 污染。二是len(merged) 3时直接跳过工具变量少于 3 个时 MR-Egger 回归没有意义IVW 也不稳定。跑完批量分析后p 值要做多重检验校正Benjamini-Hochberg FDR否则两万个基因里按 p0.05 筛会得到上千个假阳性。4. 避坑与排查因果推断在转录组里最容易翻车的五个地方4.1 弱工具变量导致估计量向朴素回归靠拢现象IVW 估计的因果效应和普通线性回归的系数几乎一样F 统计量在 10 附近徘徊。原因工具变量对吸烟行为的解释力太弱估计量被弱工具偏倚拉向朴素回归结果因果解释不成立。解决把 F 阈值提到 20 甚至 30 做敏感性分析如果合格 SNP 太少考虑用多基因风险评分PRS作为工具变量把多个弱效应 SNP 合并成一个强工具。4.2 群体分层导致工具变量与结局的虚假关联现象MR-Egger 截距显著多效性检验通不过但换一个祖先群体数据后结果完全反转。原因GWAS 汇总数据和你的转录组样本来自不同祖先群体等位基因频率和 LD 结构不匹配工具变量与结局的关联被群体分层污染。解决用与转录组样本祖先一致的 GWAS 数据如果做不到至少在分析中纳入前几个遗传主成分作为协变量并在论文里明确写清群体来源。4.3 细胞类型比例作为碰撞变量被错误调整现象调整细胞比例后因果效应消失不调整时效应显著。原因细胞比例可能是吸烟的下游后果吸烟改变气道细胞组成把它作为协变量调整会阻断部分因果路径导致效应被低估。这在因果图里叫“过度调整偏倚”。解决先画因果有向无环图DAG判断细胞比例是混杂还是中介。如果是中介不要调整如果既有混杂又有中介成分考虑用中介分析拆解直接和间接效应。4.4 表达矩阵批次效应未校正导致假阳性现象因果效应显著的基因富集在某个测序批次换批次后消失。原因批次与吸烟状态偶然相关比如吸烟样本恰好集中在某个批次批次效应被误读为吸烟的因果效应。解决用 ComBat 或 limma 的removeBatchEffect做批次校正校正后再跑 MR同时在协变量里保留批次作为固定效应双重保险。4.5 工具变量多效性通过非吸烟通路影响基因表达现象某个 SNP 已知与肺功能或炎症直接相关同时又被选为吸烟的工具变量。原因排他性假设不成立SNP 通过独立于吸烟的通路影响基因表达因果估计被高估。解决查 GWAS Catalog 或 PhenoScanner把与肺功能、炎症、自身免疫病显著相关的 SNP 剔除用 MR-PRESSO 做全局多效性检验把离群 SNP 剔除后重跑。5. 进阶技巧用中介分析和组织特异性验证把因果链补完整跑完 MR 得到一批因果效应显著的基因后下一步要回答的是“吸烟通过什么通路影响这些基因”。中介分析可以把因果链拆成“吸烟→中介变量→基因表达”的间接效应和“吸烟→基因表达”的直接效应。常用中介变量包括 DNA 甲基化、炎症因子、氧化应激标志物。实操中我一般用mediation包做先跑暴露对中介的回归再跑中介对结局的回归最后用 bootstrap 估计间接效应的置信区间。组织特异性验证是另一个容易被忽略的环节。MR 估计的是吸烟对肺组织基因表达的因果效应但工具变量来自全血或全基因组的 GWAS如果某个基因只在肺上皮细胞表达全血数据里的工具变量可能根本不调控这个基因。验证方法是查 GTEx 或单细胞转录组数据确认工具变量附近的 SNP 在肺组织中是 eQTL表达数量性状位点。如果工具变量在肺组织里没有 eQTL 证据因果估计的可信度要打折扣。# 用 pyGTEx 或本地 GTEx 数据查 eQTL 证据示意 # 假设 gtex_eqtl 有列: snp, gene, tissue, pval lung_eqtl gtex_eqtl[ (gtex_eqtl[tissue] Lung) (gtex_eqtl[gene].isin(results_df[gene])) (gtex_eqtl[pval] 0.05) ] # 统计有多少因果显著基因在肺组织中有 eQTL 支持 supported results_df[results_df[gene].isin(lung_eqtl[gene])] print(f因果显著基因数: {len(results_df)}, 肺组织 eQTL 支持: {len(supported)})最后说一个我自己的习惯每次跑完 MR先不看 p 值最小的基因而是先看 F 统计量分布和 MR-Egger 截距的 QQ 图。如果 F 分布左尾太厚说明弱工具问题严重后面所有结果都要打问号如果 Egger 截距的 QQ 图偏离对角线说明多效性普遍存在IVW 结果不能直接信。这个检查顺序帮我省过好几次返工。希望帮到你。本文还有配套的精品资源点击获取
简介这份资源是面向生物信息学与统计分析学习者的R语言实战项目聚焦用因果推断方法评估吸烟对健康肺组织基因表达的短期因果效应适合具备一定R基础、想了解观察性数据因果分析流程的研究者与高年级学生。压缩包共45个文件约3.1MB包含19个png结果图、13个txt文本、8个xlsx表格、2个csv数据、1个pdf文档及R脚本等覆盖数据表、分析结果与可视化输出便于对照复现。项目涉及工具变量、倾向得分匹配、DAG等因果推断思路并结合RNA-seq数据的质控、比对、计数与差异表达分析以及通路富集和结果验证环节能帮助读者理解从原始表达数据到因果效应量化的完整链路。目前已有148人学习可作为因果推断与基因表达分析结合的小型案例参考。1. 吸烟与基因表达因果推断能回答相关分析答不了的问题做转录组分析的人多半遇到过这种尴尬吸烟样本和健康肺组织样本的差异表达基因列表拉出来几百个富集分析也漂亮但审稿人一句“这是相关还是因果”就把人问住了。常规差异表达分析DEG只能告诉你“吸烟者肺组织里基因A表达量低”但没法排除反向因果——比如某个基因低表达的人本身更容易吸烟或者存在一个共同混杂因素年龄、炎症基线、组织细胞比例同时影响吸烟行为和基因表达。SmokingGeneExpression 这个方向要做的就是用因果推断框架去估计吸烟对健康肺组织中基因表达的短期因果效应而不是停留在相关性层面。它适合三类人手里有吸烟/非吸烟人群肺组织转录组数据但被审稿人追问因果性的分析者想入门孟德尔随机化MR或工具变量法但缺少真实场景练手的生信人以及需要为公共卫生政策提供更硬证据的研究者。核心思路是把遗传变异作为工具变量或者用倾向性评分、断点回归等准实验设计把“吸烟”这个暴露变量从混杂网络中剥离出来再去看它对每个基因表达量的因果效应量。2. 因果推断选型为什么工具变量法在肺组织转录组里更稳2.1 从相关到因果三个必须跨过的门槛在肺组织转录组场景里做因果推断要跨过三个门槛。第一是混杂控制吸烟者往往年龄偏大、炎症水平更高、气道细胞比例不同这些都会直接影响基因表达如果不控制估计出的效应就是偏的。第二是反向因果某些基因的表达状态可能通过神经递质通路影响尼古丁依赖倾向导致“基因表达低→更容易吸烟”而不是反过来。第三是测量误差吸烟量通常靠问卷自报存在回忆偏倚而基因表达又有技术噪声两者叠加会让朴素回归的系数向零衰减。常规做法是用多元回归把年龄、性别、BMI、细胞类型比例作为协变量塞进去但这只能控制已测混杂未测混杂比如遗传背景、早年环境暴露依然会污染估计。更稳的做法是引入工具变量找一个与吸烟行为强相关、但与肺组织基因表达只能通过吸烟间接相关的遗传位点通常是 rs 开头的 SNP。这就是孟德尔随机化MR的核心逻辑——遗传变异在受精时随机分配相当于自然界的随机对照试验能绕过大部分混杂和反向因果。2.2 工具变量法的三个核心假设与检验方法工具变量法不是随便找个 SNP 就能用它依赖三个假设每个假设都有对应的检验手段。假设含义检验方法不满足的后果相关性工具变量与吸烟行为强相关第一阶段 F 统计量 10弱工具偏倚估计量向朴素回归靠拢独立性工具变量与混杂因素无关检查与年龄、性别、细胞比例的关联混杂重新引入因果解释失效排他性工具变量只通过吸烟影响基因表达多效性检验MR-Egger 截距、MR-PRESSO高估或低估因果效应第一阶段 F 统计量是最先要看的。实操中我一般用TwoSampleMR包跑F 值低于 10 的 SNP 直接剔除不要心疼。排他性假设最难验证通常用 MR-Egger 回归的截距项做粗略判断截距显著偏离零说明存在方向性多效性这时候要么换工具变量子集要么改用加权中位数法。2.3 短期因果效应的定义与时间窗口选择标题里“短期因果关系”是个容易被忽略的限定词。肺组织基因表达是动态的吸烟后几小时到几天的早期响应基因比如 CYP1A1、AHRR 通路和吸烟几十年的累积损伤基因比如与 COPD 相关的 MMP 家族因果结构完全不同。短期因果效应通常指暴露后 24 小时到 4 周内的表达变化对应的是急性炎症、氧化应激和代谢激活通路。在数据层面如果你用的是横断面人群样本其实很难严格区分短期和长期。常见做法是优先选择有明确吸烟状态记录当前吸烟/已戒烟/从不吸烟且戒烟时间在 4 周以内的样本作为“短期暴露组”或者用干预实验比如体外气液界面暴露模型的数据做验证。如果只有横断面数据就在论文里明确写“本研究估计的是当前吸烟状态对基因表达的因果效应时间尺度受限于横断面设计”不要硬说成短期。3. 从原始数据到因果估计一套可复现的 MR 分析流程3.1 数据准备暴露、结局与协变量三张表假设你手里有一份肺组织转录组表达矩阵行是基因列是样本和一份样本元数据吸烟状态、年龄、性别、细胞比例。第一步是把数据整理成 MR 分析需要的三张表暴露表吸烟状态和工具变量 SNP、结局表每个基因的表达量、协变量表年龄、性别、测序批次、细胞比例。import pandas as pd import numpy as np # 读取表达矩阵和元数据 expr pd.read_csv(lung_expr_matrix.csv, index_col0) # 行基因列样本 meta pd.read_csv(sample_metadata.csv, index_col0) # 行样本 # 对齐样本顺序这一步不做后面全错 common_samples expr.columns.intersection(meta.index) expr expr[common_samples] meta meta.loc[common_samples] # 构建暴露表吸烟状态编码为 0/1 exposure meta[[smoking_status]].copy() exposure[smoking_status] exposure[smoking_status].map( {never: 0, former: 0, current: 1} ) # 协变量表年龄、性别、批次、细胞比例 covariates meta[[age, sex, batch, cell_prop_epithelial, cell_prop_immune]] # 结局表选一个目标基因比如 CYP1A1 target_gene CYP1A1 outcome expr.loc[target_gene].to_frame(expression) print(f样本数: {len(common_samples)}, 暴露分布:\n{exposure[smoking_status].value_counts()})这段代码的关键在样本对齐。表达矩阵的列名和元数据的行名如果不一致后续所有合并都会错位而且 pandas 不会报错只会静默产生 NaN。我一般会在合并后打印样本数做二次确认。暴露编码里把“已戒烟”归到 0是因为短期因果效应关注当前暴露已戒烟者的基因表达可能已经部分恢复混在一起会稀释效应。3.2 工具变量筛选从 GWAS 汇总数据到独立 SNP 集工具变量通常来自公开的吸烟行为 GWAS 汇总数据。筛选标准一般设三条与吸烟行为关联的 p 值小于 5e-8、SNP 之间连锁不平衡 r² 小于 0.001、距离至少 10000 bp。下面是一个筛选脚本的骨架。# 假设 gwas_df 包含列: snp, beta, se, pval, chr, pos gwas_df pd.read_csv(smoking_gwas_summary.csv) # 第一步全基因组显著性筛选 sig_snps gwas_df[gwas_df[pval] 5e-8].copy() # 第二步按染色体和位置排序做 LD 剪枝简化版实际用 plink sig_snps sig_snps.sort_values([chr, pos]) pruned [] last_pos {} for _, row in sig_snps.iterrows(): key row[chr] if key not in last_pos or (row[pos] - last_pos[key]) 10000: pruned.append(row) last_pos[key] row[pos] iv_df pd.DataFrame(pruned) print(f独立工具变量数量: {len(iv_df)}) # 第三步计算第一阶段 F 统计量 iv_df[F_stat] (iv_df[beta] / iv_df[se]) ** 2 iv_df iv_df[iv_df[F_stat] 10] print(f通过 F10 的工具变量: {len(iv_df)})LD 剪枝这里用的是简化版距离法实际分析中应该用 PLINK 的--indep-pairwise 1000 10 0.001做标准剪枝。F 统计量公式是 (beta/se)²等价于 t 统计量的平方。F 大于 10 是经验阈值低于这个值说明工具变量解释的暴露变异太少估计结果不可信。注意这里 beta 和 se 都来自 GWAS 汇总数据不是你自己样本算的所以要求 GWAS 和你的转录组样本来自同一祖先群体否则等位基因频率不匹配会导致工具变量失效。3.3 因果效应估计Wald 比值法与 MR-Egger 回归拿到工具变量后最基础的因果效应估计是 Wald 比值法用工具变量对结局的效应除以工具变量对暴露的效应。单个 SNP 的 Wald 估计量方差可以用 delta 方法近似。多个 SNP 时用逆方差加权IVW合并。# 假设 iv_df 有 beta_exposure, se_exposure # outcome_df 有每个 SNP 对基因表达的 beta_outcome, se_outcome merged iv_df.merge(outcome_df, onsnp) # Wald 比值法 merged[wald_ratio] merged[beta_outcome] / merged[beta_exposure] merged[wald_se] np.sqrt( (merged[se_outcome] ** 2) / (merged[beta_exposure] ** 2) (merged[beta_outcome] ** 2) * (merged[se_exposure] ** 2) / (merged[beta_exposure] ** 4) ) # IVW 合并 weights 1 / merged[wald_se] ** 2 ivw_estimate np.sum(weights * merged[wald_ratio]) / np.sum(weights) ivw_se np.sqrt(1 / np.sum(weights)) print(fIVW 因果效应估计: {ivw_estimate:.4f} (SE: {ivw_se:.4f})) print(f95% CI: [{ivw_estimate - 1.96*ivw_se:.4f}, {ivw_estimate 1.96*ivw_se:.4f}])IVW 的假设是所有工具变量都是有效的如果存在多效性IVW 会被带偏。这时候要跑 MR-Egger 回归做敏感性分析Egger 回归的截距项如果显著不为零说明存在方向性多效性IVW 结果不可信。另一个常用方法是加权中位数法它在部分工具变量无效时仍然稳健但精度略低。实操中我一般三个方法都跑如果 IVW 和加权中位数方向一致且 Egger 截距不显著才认为因果证据比较硬。3.4 多基因并行把单基因流程封装成批量函数肺组织转录组有两万多个基因不可能一个一个手动跑。把上面的流程封装成函数对每个基因循环执行同时记录效应量、标准误、p 值和多效性检验结果。def mr_analysis_for_gene(gene_name, expr_df, iv_df, covariates_df): 对单个基因跑 MR 分析返回因果效应估计和诊断指标 outcome expr_df.loc[gene_name].to_frame(expression) merged iv_df.merge(outcome, left_onsnp, right_indexTrue, howinner) if len(merged) 3: return None # 工具变量太少跳过 # IVW merged[wald_ratio] merged[beta_outcome] / merged[beta_exposure] merged[wald_se] np.sqrt( (merged[se_outcome] ** 2) / (merged[beta_exposure] ** 2) (merged[beta_outcome] ** 2) * (merged[se_exposure] ** 2) / (merged[beta_exposure] ** 4) ) weights 1 / merged[wald_se] ** 2 ivw_est np.sum(weights * merged[wald_ratio]) / np.sum(weights) ivw_se np.sqrt(1 / np.sum(weights)) # MR-Egger 截距 from scipy import stats slope, intercept, r, p, se stats.linregress( merged[beta_exposure], merged[beta_outcome] ) return { gene: gene_name, ivw_beta: ivw_est, ivw_se: ivw_se, ivw_pval: 2 * stats.norm.sf(abs(ivw_est / ivw_se)), egger_intercept: intercept, egger_pval: p, n_iv: len(merged) } # 批量跑前 100 个高表达基因做测试 high_expr_genes expr.mean(axis1).sort_values(ascendingFalse).head(100).index results [mr_analysis_for_gene(g, expr, iv_df, covariates) for g in high_expr_genes] results_df pd.DataFrame([r for r in results if r is not None]) print(results_df.sort_values(ivw_pval).head(10))这个批量函数里有两个细节值得注意。一是howinner合并确保只有同时出现在工具变量表和表达矩阵里的 SNP 才参与分析避免 NaN 污染。二是len(merged) 3时直接跳过工具变量少于 3 个时 MR-Egger 回归没有意义IVW 也不稳定。跑完批量分析后p 值要做多重检验校正Benjamini-Hochberg FDR否则两万个基因里按 p0.05 筛会得到上千个假阳性。4. 避坑与排查因果推断在转录组里最容易翻车的五个地方4.1 弱工具变量导致估计量向朴素回归靠拢现象IVW 估计的因果效应和普通线性回归的系数几乎一样F 统计量在 10 附近徘徊。原因工具变量对吸烟行为的解释力太弱估计量被弱工具偏倚拉向朴素回归结果因果解释不成立。解决把 F 阈值提到 20 甚至 30 做敏感性分析如果合格 SNP 太少考虑用多基因风险评分PRS作为工具变量把多个弱效应 SNP 合并成一个强工具。4.2 群体分层导致工具变量与结局的虚假关联现象MR-Egger 截距显著多效性检验通不过但换一个祖先群体数据后结果完全反转。原因GWAS 汇总数据和你的转录组样本来自不同祖先群体等位基因频率和 LD 结构不匹配工具变量与结局的关联被群体分层污染。解决用与转录组样本祖先一致的 GWAS 数据如果做不到至少在分析中纳入前几个遗传主成分作为协变量并在论文里明确写清群体来源。4.3 细胞类型比例作为碰撞变量被错误调整现象调整细胞比例后因果效应消失不调整时效应显著。原因细胞比例可能是吸烟的下游后果吸烟改变气道细胞组成把它作为协变量调整会阻断部分因果路径导致效应被低估。这在因果图里叫“过度调整偏倚”。解决先画因果有向无环图DAG判断细胞比例是混杂还是中介。如果是中介不要调整如果既有混杂又有中介成分考虑用中介分析拆解直接和间接效应。4.4 表达矩阵批次效应未校正导致假阳性现象因果效应显著的基因富集在某个测序批次换批次后消失。原因批次与吸烟状态偶然相关比如吸烟样本恰好集中在某个批次批次效应被误读为吸烟的因果效应。解决用 ComBat 或 limma 的removeBatchEffect做批次校正校正后再跑 MR同时在协变量里保留批次作为固定效应双重保险。4.5 工具变量多效性通过非吸烟通路影响基因表达现象某个 SNP 已知与肺功能或炎症直接相关同时又被选为吸烟的工具变量。原因排他性假设不成立SNP 通过独立于吸烟的通路影响基因表达因果估计被高估。解决查 GWAS Catalog 或 PhenoScanner把与肺功能、炎症、自身免疫病显著相关的 SNP 剔除用 MR-PRESSO 做全局多效性检验把离群 SNP 剔除后重跑。5. 进阶技巧用中介分析和组织特异性验证把因果链补完整跑完 MR 得到一批因果效应显著的基因后下一步要回答的是“吸烟通过什么通路影响这些基因”。中介分析可以把因果链拆成“吸烟→中介变量→基因表达”的间接效应和“吸烟→基因表达”的直接效应。常用中介变量包括 DNA 甲基化、炎症因子、氧化应激标志物。实操中我一般用mediation包做先跑暴露对中介的回归再跑中介对结局的回归最后用 bootstrap 估计间接效应的置信区间。组织特异性验证是另一个容易被忽略的环节。MR 估计的是吸烟对肺组织基因表达的因果效应但工具变量来自全血或全基因组的 GWAS如果某个基因只在肺上皮细胞表达全血数据里的工具变量可能根本不调控这个基因。验证方法是查 GTEx 或单细胞转录组数据确认工具变量附近的 SNP 在肺组织中是 eQTL表达数量性状位点。如果工具变量在肺组织里没有 eQTL 证据因果估计的可信度要打折扣。# 用 pyGTEx 或本地 GTEx 数据查 eQTL 证据示意 # 假设 gtex_eqtl 有列: snp, gene, tissue, pval lung_eqtl gtex_eqtl[ (gtex_eqtl[tissue] Lung) (gtex_eqtl[gene].isin(results_df[gene])) (gtex_eqtl[pval] 0.05) ] # 统计有多少因果显著基因在肺组织中有 eQTL 支持 supported results_df[results_df[gene].isin(lung_eqtl[gene])] print(f因果显著基因数: {len(results_df)}, 肺组织 eQTL 支持: {len(supported)})最后说一个我自己的习惯每次跑完 MR先不看 p 值最小的基因而是先看 F 统计量分布和 MR-Egger 截距的 QQ 图。如果 F 分布左尾太厚说明弱工具问题严重后面所有结果都要打问号如果 Egger 截距的 QQ 图偏离对角线说明多效性普遍存在IVW 结果不能直接信。这个检查顺序帮我省过好几次返工。希望帮到你。本文还有配套的精品资源点击获取
