608 微生物组因果推断¶

========== 因果推断方法层次 ==========

金标准：随机对照实验（RCT）
  └── 随机给人吃不同的菌，看谁得T2D → 不可能！

银标准：孟德尔随机化（MR）
  └── 用基因变异代替随机分配 → 目前最可行的方法

铜标准：纵向队列 + 调解分析
  └── 先测菌群，后看疾病发生 → 有用但仍无法排除所有混杂

观察性研究（横断面）
  └── 同时测菌群和疾病 → 只能说相关，不能说因果

传统观察性研究（有混杂风险）：
菌群丰度 ----?---→ T2D
    ↑                ↑
    └── 混杂因素 ──┘（如饮食）

孟德尔随机化（排除混杂）：
基因变异(SNP) ──→ 菌群丰度 ──→ T2D
    ↑                              ↑
    X 基因不受混杂因素影响         X 基因不直接影响T2D

========== MR三大假设（必须全部满足） ==========

假设1：相关性（Relevance）
  SNP必须与暴露（菌群丰度）显著相关
  → 用F统计量检验，F > 10才合格

假设2：独立性（Independence）
  SNP与混杂因素无关
  → 基因在受精时随机分配，理论上满足

假设3：排他性（Exclusion restriction）
  SNP只通过暴露影响结局，不能直接影响结局
  → 最难验证！需要用多种MR方法交叉验证（敏感性分析）

========== 微生物组MR的GWAS数据来源 ==========

暴露数据（菌群GWAS）：
├── MiBioGen联盟（最大的肠道菌群GWAS，N=18,340）
│   └── 分类水平：门/纲/目/科/属（211个分类单元）
│   └── 下载：https://mibiogen.gcc.rug.nl/
├── FinnGen + UK Biobank（部分有菌群GWAS数据）
└── DutchMicrobiomeProject（荷兰人群）

结局数据（疾病GWAS）：
├── UK Biobank（最大的生物样本库）
├── FinnGen（芬兰人群）
├── DIAGRAM联盟（T2D专用GWAS）
└── 各种疾病专用GWAS汇总统计数据

获取GWAS汇总数据的数据库：
├── IEU OpenGWAS：https://gwas.mrcieu.ac.uk/
├── GWAS Catalog：https://www.ebi.ac.uk/gwas/
└── PhenoScanner：http://www.phenoscanner.medschl.cam.ac.uk/

# ========== 步骤1：安装和加载R包 ==========
# install.packages("remotes")
# remotes::install_github("MRCIEU/TwoSampleMR")  # 安装TwoSampleMR

library(TwoSampleMR)  # 两样本MR核心包
library(dplyr)         # 数据操作
library(ggplot2)       # 可视化

# ========== 步骤2：获取暴露数据（菌群SNP） ==========
# 方法1：从IEU OpenGWAS数据库获取
# 搜索MiBioGen的菌群GWAS数据
available_outcomes()  # 查看可用的GWAS数据集

# 提取特定菌属的GWAS显著SNP（如Bifidobacterium）
exposure_dat <- extract_instruments(
  outcomes = "ebi-a-GCST90027446",  # Bifidobacterium属的GWAS ID（示例）
  p1 = 1e-5,                         # P值阈值（微生物GWAS用宽松阈值，因为位点少）
  clump = TRUE,                      # LD剪枝（去除连锁不平衡的SNP）
  r2 = 0.001,                        # LD r2阈值（<0.001表示几乎无连锁）
  kb = 10000                          # 剪枝窗口10Mb
)

# 方法2：从本地文件读取（如果GWAS数据已下载）
exposure_dat <- read_exposure_data(
  filename = "bifidobacterium_gwas.txt",  # 本地GWAS文件
  sep = "\t",                              # 分隔符
  snp_col = "SNP",                         # SNP列名
  beta_col = "BETA",                       # 效应值列名
  se_col = "SE",                           # 标准误列名
  effect_allele_col = "A1",                # 效应等位基因列名
  other_allele_col = "A2",                 # 参考等位基因列名
  pval_col = "P",                          # P值列名
  eaf_col = "EAF"                          # 效应等位基因频率列名
)

# 检查工具变量强度
cat("SNP数量:", nrow(exposure_dat), "\n")  # 打印SNP数

# 计算F统计量（衡量工具变量强度）
exposure_dat$F_stat <- (exposure_dat$beta.exposure)^2 / (exposure_dat$se.exposure)^2
cat("平均F统计量:", mean(exposure_dat$F_stat), "\n")  # F>10才是强工具变量
cat("弱工具变量(F<10)数:", sum(exposure_dat$F_stat < 10), "\n")

# 去除弱工具变量
exposure_dat <- exposure_dat[exposure_dat$F_stat >= 10, ]  # 只保留F>=10的SNP

# ========== 步骤3：获取结局数据（疾病SNP） ==========
# 从IEU数据库提取T2D的GWAS数据
outcome_dat <- extract_outcome_data(
  snps = exposure_dat$SNP,            # 用暴露的SNP列表
  outcomes = "ebi-a-GCST006867",      # T2D GWAS ID（DIAGRAM联盟）
  proxies = TRUE,                      # 如果SNP不存在，寻找代理SNP
  rsq = 0.8,                          # 代理SNP的LD r2阈值
  maf_threshold = 0.01                 # 最小等位基因频率>1%
)

# ========== 步骤4：数据协调（Harmonisation） ==========
# 确保暴露和结局的效应方向一致
harmonised_dat <- harmonise_data(
  exposure_dat = exposure_dat,  # 暴露数据
  outcome_dat = outcome_dat,    # 结局数据
  action = 2                    # action=2：尝试推断正向链（推荐）
)

# 检查协调后的数据
cat("协调后可用SNP数:", nrow(harmonised_dat), "\n")

# 去除回文SNP（palindromic SNPs，如A/T或C/G）
harmonised_dat <- harmonised_dat[harmonised_dat$palindromic == FALSE |
                                 harmonised_dat$ambiguous == FALSE, ]

# ========== 步骤5：MR分析（核心步骤） ==========
# 运行多种MR方法
mr_results <- mr(
  harmonised_dat,              # 协调后的数据
  method_list = c(
    "mr_ivw",                  # 逆方差加权法（主要方法）
    "mr_egger_regression",     # MR-Egger回归（检测水平多效性）
    "mr_weighted_median",      # 加权中位数法（允许50%无效IV）
    "mr_weighted_mode",        # 加权众数法
    "mr_simple_mode"           # 简单众数法
  )
)

# 查看结果
print(mr_results[, c("method", "nsnp", "b", "se", "pval")])
# b: 效应值（正=风险因素，负=保护因素）
# pval: P<0.05表示因果关系显著

# 结果解读：
# 如果IVW和其他方法方向一致且IVW P<0.05 → 支持因果关系
# 如果各方法结果不一致 → 可能存在多效性，结果不可靠

# OR值（优势比）
mr_results$OR <- exp(mr_results$b)       # 计算OR值
mr_results$OR_lci <- exp(mr_results$b - 1.96 * mr_results$se)  # OR 95%CI下限
mr_results$OR_uci <- exp(mr_results$b + 1.96 * mr_results$se)  # OR 95%CI上限
print(mr_results[, c("method", "OR", "OR_lci", "OR_uci", "pval")])

# ========== 步骤6：敏感性分析（必须做！） ==========

# 1. 异质性检验（Cochran's Q）
heterogeneity <- mr_heterogeneity(harmonised_dat)  # 检验SNP间异质性
print(heterogeneity)
# Q_pval < 0.05 → 存在异质性，需要用随机效应模型

# 2. 水平多效性检验（MR-Egger截距）
pleiotropy <- mr_pleiotropy_test(harmonised_dat)  # 检验水平多效性
print(pleiotropy)
# 截距P值 < 0.05 → 存在水平多效性，IVW结果可能有偏

# 3. 留一法（Leave-one-out）
loo <- mr_leaveoneout(harmonised_dat)  # 逐一去掉每个SNP重新分析
# 如果去掉某个SNP后结果反转 → 该SNP可能是异常值

# 4. 漏斗图（Funnel plot）
# 用于检测发表偏倚和异常SNP
mr_funnel_plot(
  mr_singlesnp(harmonised_dat)  # 单SNP分析结果
)

# ========== 步骤7：MR结果可视化 ==========

# 1. 散点图（最重要的MR可视化）
p_scatter <- mr_scatter_plot(
  mr_results,       # MR结果
  harmonised_dat     # 协调后的数据
)
print(p_scatter[[1]])  # 显示散点图
# X轴：SNP对菌群的效应
# Y轴：SNP对T2D的效应
# 如果因果关系存在，点应该沿着回归线分布

# 2. 森林图
p_forest <- mr_forest_plot(
  mr_singlesnp(harmonised_dat)  # 单SNP分析
)
print(p_forest[[1]])  # 显示森林图

# 3. 漏斗图
p_funnel <- mr_funnel_plot(
  mr_singlesnp(harmonised_dat)  # 单SNP分析
)
print(p_funnel[[1]])

# 4. 留一法图
p_loo <- mr_leaveoneout_plot(
  mr_leaveoneout(harmonised_dat)  # 留一法分析
)
print(p_loo[[1]])

# 保存所有图形
ggsave("mr_scatter.pdf", p_scatter[[1]], width = 8, height = 6)    # 散点图
ggsave("mr_forest.pdf", p_forest[[1]], width = 8, height = 10)     # 森林图
ggsave("mr_funnel.pdf", p_funnel[[1]], width = 8, height = 6)      # 漏斗图
ggsave("mr_leaveoneout.pdf", p_loo[[1]], width = 8, height = 10)   # 留一法图

# ========== 步骤8：反向MR（检验反向因果） ==========
# 正向：菌群 → T2D
# 反向：T2D → 菌群

# 反向MR：以T2D为暴露，菌群为结局
exposure_reverse <- extract_instruments(
  outcomes = "ebi-a-GCST006867",  # T2D GWAS作为暴露
  p1 = 5e-8,                       # 标准GWAS阈值（T2D的位点充足）
  clump = TRUE,
  r2 = 0.001,
  kb = 10000
)

outcome_reverse <- extract_outcome_data(
  snps = exposure_reverse$SNP,
  outcomes = "ebi-a-GCST90027446"  # Bifidobacterium GWAS作为结局
)

harmonised_reverse <- harmonise_data(
  exposure_dat = exposure_reverse,
  outcome_dat = outcome_reverse,
  action = 2
)

# 运行反向MR
mr_reverse <- mr(harmonised_reverse)
print(mr_reverse[, c("method", "b", "se", "pval")])

# 结果解读：
# 正向MR显著 + 反向MR不显著 → 菌群 → T2D（单向因果）
# 正向MR显著 + 反向MR显著   → 双向因果关系
# 正向MR不显著 + 反向MR显著 → T2D → 菌群（反向因果！）
# 两者都不显著              → 无因果关系（可能只是共有混杂）

# ========== 纵向数据的调解分析（Mediation Analysis） ==========
# 适用场景：有时间序列的队列数据

library(mediation)  # 调解分析包

# 假设有纵向数据：
# T0: 基线菌群丰度、协变量
# T1: 代谢指标（如HOMA-IR）
# T2: T2D发病情况

# 模型1：暴露 → 调解变量
mediator_model <- lm(
  HOMA_IR_T1 ~ Bacteroides_T0 + age + sex + BMI,  # 菌群预测HOMA-IR
  data = longitudinal_data
)

# 模型2：暴露 + 调解变量 → 结局
outcome_model <- glm(
  T2D_T2 ~ Bacteroides_T0 + HOMA_IR_T1 + age + sex + BMI,  # 菌群+HOMA-IR预测T2D
  data = longitudinal_data,
  family = binomial  # 二分类结局用logistic回归
)

# 调解分析
med_result <- mediate(
  mediator_model,         # 调解变量模型
  outcome_model,          # 结局模型
  treat = "Bacteroides_T0",  # 处理变量（暴露）
  mediator = "HOMA_IR_T1",   # 调解变量
  boot = TRUE,               # 用bootstrap估计CI
  sims = 1000                # bootstrap次数
)

summary(med_result)
# ACME: 平均因果调解效应（通过HOMA-IR介导的间接效应）
# ADE:  平均直接效应（不通过HOMA-IR的直接效应）
# Total Effect: 总效应 = ACME + ADE
# Prop. Mediated: 调解比例（ACME/Total Effect）

# ========== 格兰杰因果检验（适用于纵向微生物组数据） ==========
library(lmtest)  # 格兰杰因果检验

# 格兰杰因果：如果菌群的过去值能预测未来的代谢指标，
# 且代谢指标的过去值无法解释这种预测力，则菌群"格兰杰因果"代谢指标

# 准备时间序列数据
ts_data <- data.frame(
  time = 1:20,                           # 20个时间点
  bacteroides = rnorm(20, mean = 5),     # Bacteroides丰度（示例）
  glucose = rnorm(20, mean = 100)        # 血糖值（示例）
)

# 格兰杰因果检验
granger_test <- grangertest(
  glucose ~ bacteroides,  # 检验bacteroides是否格兰杰因果glucose
  order = 2,               # 滞后阶数
  data = ts_data
)
print(granger_test)
# P < 0.05 → bacteroides的过去值能显著预测glucose

========== 两样本MR标准流程 ==========
选择暴露/结局GWAS → 提取工具变量(SNP)
    → F统计量检验(F>10) → 提取结局数据
    → 数据协调(harmonise) → MR分析(IVW+敏感性)
    → 异质性检验 → 多效性检验 → 留一法
    → 反向MR → 结果可视化

========== MR方法对比 ==========
方法              假设           推荐度   解读
IVW              无多效性       主要方法  最强统计效力
MR-Egger         允许多效性     敏感性    截距检验多效性
加权中位数        ≤50%无效IV    敏感性    稳健
加权众数         多数IV有效     补充      最稳健

========== 四张必出图 ==========
1. 散点图：SNP的暴露效应 vs 结局效应（核心图）
2. 森林图：单个SNP和汇总的效应值
3. 漏斗图：检测异常SNP和发表偏倚
4. 留一法图：去掉每个SNP后的汇总效应

========== 核心R包 ==========
TwoSampleMR   ：两样本MR全流程（MRCIEU团队）
MendelianRandomization ：备选MR包
MR-PRESSO     ：检测异常SNP（水平多效性outlier）
MVMR          ：多变量MR
mediation     ：因果调解分析

========== GWAS数据获取 ==========
IEU OpenGWAS ：https://gwas.mrcieu.ac.uk/（最全面）
MiBioGen     ：https://mibiogen.gcc.rug.nl/（菌群GWAS）
FinnGen      ：https://finngen.gitbook.io/（芬兰人群）
GWAS Catalog ：https://www.ebi.ac.uk/gwas/（EBI维护）