162_微生物组与宿主遗传互作¶

library(OpenMx)

# === ACE模型: 估计微生物丰度的遗传度 ===
# 数据: 同卵(MZ)和异卵(DZ)双胞胎的微生物丰度

# 准备数据
twin_data <- data.frame(
  MZ_twin1 = rnorm(100),   # MZ双胞胎1的某菌属CLR丰度
  MZ_twin2 = rnorm(100),   # MZ双胞胎2的某菌属CLR丰度
  DZ_twin1 = rnorm(100),   # DZ双胞胎1
  DZ_twin2 = rnorm(100)    # DZ双胞胎2
)

# 计算组内相关系数(ICC)
mz_cor <- cor(twin_data$MZ_twin1, twin_data$MZ_twin2)
dz_cor <- cor(twin_data$DZ_twin1, twin_data$DZ_twin2)

# Falconer公式估计遗传度
h2_falconer <- 2 * (mz_cor - dz_cor)
c2 <- 2 * dz_cor - mz_cor  # 共享环境
e2 <- 1 - h2_falconer - c2  # 非共享环境

cat(sprintf("遗传度 h2: %.3f\n共享环境 c2: %.3f\n非共享环境 e2: %.3f\n",
            h2_falconer, c2, e2))

# === 批量计算所有菌属的遗传度 ===
# 对每个菌属(genus)运行ACE模型
genera_list <- colnames(microbiome_clr)
heritability_results <- data.frame(
  Genus = genera_list,
  h2 = numeric(length(genera_list)),
  h2_CI_lower = numeric(length(genera_list)),
  h2_CI_upper = numeric(length(genera_list)),
  p_value = numeric(length(genera_list))
)

# 使用更精确的SEM模型 (OpenMx)
for (i in seq_along(genera_list)) {
  # ... OpenMx ACE模型拟合 ...
  # 省略详细代码, 参考OpenMx twin modeling vignette
}

# 已知高遗传度菌属:
# Christensenellaceae: h2 ~ 0.40
# Bifidobacterium: h2 ~ 0.30
# Methanobrevibacter: h2 ~ 0.35

# === 数据准备 ===
# 1. 宿主基因型数据 (PLINK格式: .bed/.bim/.fam)
# 标准GWAS QC
plink --bfile genotypes \
  --mind 0.05 \
  --geno 0.05 \
  --maf 0.05 \
  --hwe 1e-6 \
  --make-bed \
  --out genotypes_qc

# 2. 微生物组表型数据 (CLR转换后的丰度)
# 每个菌属作为一个"表型"

# === 使用PLINK或R进行mbQTL ===
library(MatrixEQTL)
# MatrixEQTL适合大规模QTL分析

# 加载数据
snps <- SlicedData$new()
snps$LoadFile("genotypes_matrix.txt")  # SNP矩阵 (SNP × 样本)

microbiome <- SlicedData$new()
microbiome$LoadFile("microbiome_clr.txt")  # 微生物CLR丰度 (taxa × 样本)

covariates <- SlicedData$new()
covariates$LoadFile("covariates.txt")  # 协变量: 年龄、性别、BMI、PCA前10个主成分

# 运行mbQTL
me <- Matrix_eQTL_main(
  snps = snps,
  gene = microbiome,  # 这里"gene"就是微生物丰度
  cvrt = covariates,
  output_file_name = "mbqtl_results.txt",
  pvOutputThreshold = 1e-5,  # 存储的P值阈值
  useModel = modelLINEAR,
  errorCovariance = numeric(),
  verbose = TRUE,
  pvalue.hist = "qqplot"
)

# 结果
cat(sprintf("检测到 %d 个显著mbQTL (FDR<0.05)\n", 
            sum(me$all$eqtls$FDR < 0.05)))

# === 多重检验校正 ===
# SNP数 × taxa数 的多重检验负担很重
# 建议: 先对每个taxon做GWAS, 各自Bonferroni校正(5e-8)
# 或使用FDR < 0.05跨所有测试

# === 已知的mbQTL举例 ===
# LCT (chr2): 影响Bifidobacterium丰度 (乳糖酶持久性)
# FUT2 (chr19): 影响Ruminococcus torques (分泌型状态)
# ABO (chr9): 影响肠道菌群组成 (血型)
# HLA region (chr6): 广泛影响免疫相关菌属

library(TwoSampleMR)

# === 微生物组 → 疾病 的因果推断 ===
# 工具变量: mbQTL中发现的显著SNP

# 1. 暴露数据(微生物): 使用mbQTL结果作为工具变量
exposure_dat <- data.frame(
  SNP = c("rs4988235", "rs601338", "rs505922"),  # mbQTL SNPs
  beta = c(0.15, 0.12, 0.08),  # 对微生物丰度的效应
  se = c(0.03, 0.02, 0.02),
  pval = c(1e-7, 1e-8, 1e-5),
  effect_allele = c("T", "A", "G"),
  other_allele = c("C", "G", "A")
)
exposure_dat$exposure <- "Bifidobacterium"
exposure_dat$id.exposure <- "bifidobacterium"

# 格式化
exposure_formatted <- format_data(exposure_dat, type = "exposure")

# 2. 结局数据(疾病): 从GWAS summary statistics获取
outcome_dat <- extract_outcome_data(
  snps = exposure_formatted$SNP,
  outcomes = "ieu-a-1100"  # 如T2D的GWAS ID (MR-Base)
)

# 3. 协调数据
dat <- harmonise_data(exposure_formatted, outcome_dat)

# 4. MR分析
mr_results <- mr(dat, method_list = c(
  "mr_ivw",             # 逆方差加权
  "mr_egger_regression", # MR-Egger (检测pleiotropy)
  "mr_weighted_median"   # 加权中位数
))

print(mr_results)

# 5. 敏感性分析
mr_pleiotropy_test(dat)  # 多效性检验
mr_heterogeneity(dat)    # 异质性检验
mr_scatter_plot(mr_results, dat)  # 散点图

# === 反向MR: 疾病 → 微生物组 ===
# 使用疾病GWAS的SNP作为工具变量
# 检验是否是疾病导致微生物改变(而非反方向)

# === 宿主基因表达与微生物组的关联 ===
library(Maaslin2)

# Maaslin2: 多变量关联分析
# 在关联分析中同时纳入宿主基因型、表达量和微生物丰度

# 分析宿主基因表达与微生物丰度的关联
fit_expression <- Maaslin2(
  input_data = microbiome_clr_df,     # 微生物CLR丰度
  input_metadata = host_metadata,      # 含基因表达和临床变量
  output = "maaslin2_output/",
  fixed_effects = c("host_gene_expr", "age", "sex", "BMI"),
  random_effects = c("batch"),
  normalization = "NONE",  # 已做CLR
  transform = "NONE",
  analysis_method = "LM"
)

# === 中介效应分析 ===
# 测试: 宿主SNP → 宿主基因表达 → 微生物丰度
library(mediation)

# 模型1: SNP → 基因表达(中介变量)
mediator_model <- lm(gene_expression ~ SNP + age + sex, data = combined_data)

# 模型2: SNP + 基因表达 → 微生物丰度(结局)
outcome_model <- lm(microbe_abundance ~ SNP + gene_expression + age + sex, 
                     data = combined_data)

# 中介效应
med_result <- mediate(mediator_model, outcome_model,
                       treat = "SNP", mediator = "gene_expression",
                       boot = TRUE, sims = 1000)
summary(med_result)
# ACME: 平均因果中介效应
# ADE: 平均直接效应
# Total Effect: 总效应

# === mbQTL分析完整流程 ===
# 1. 基因型QC
plink --bfile raw_genotypes --mind 0.05 --geno 0.05 --maf 0.05 \
  --hwe 1e-6 --make-bed --out qc_genotypes
plink --bfile qc_genotypes --indep-pairwise 50 5 0.2 --out pruned
plink --bfile qc_genotypes --pca 10 --out pca

# 2. 微生物组数据准备 (QIIME2 → CLR)
qiime diversity core-metrics --i-table table.qza --p-sampling-depth 10000 --output-dir diversity/

# 3. 关联分析
Rscript run_mbqtl.R --geno qc_genotypes --micro microbiome_clr.csv \
  --covariates covariates.csv --output mbqtl_results/

# 4. 孟德尔随机化
Rscript run_mr.R --exposure mbqtl_instruments.csv \
  --outcome disease_gwas.csv --output mr_results/