Microbiome-GWAS关联¶

# 研究设计概览
# 样本量要求：一般需要 > 1000人（最好 > 5000）
# 数据需求：
# 1. 宿主基因组数据（GWAS芯片或WGS，imputation后）
# 2. 微生物组数据（16S rRNA或宏基因组shotgun）
# 3. 协变量信息（年龄、性别、BMI、饮食、用药等）

# 已知的强关联位点示例：
# - LCT (chr2): 影响Bifidobacterium丰度
# - FUT2 (chr19): 影响Ruminococcus等菌属
# - ABO (chr9): 影响多个菌属
# - VDR (chr12): 影响微生物多样性

# === 微生物表型准备 ===
library(phyloseq)
library(compositions)

# 读取微生物组数据
ps <- readRDS("phyloseq_object.rds")

# 1. 相对丰度（属水平）
ps_genus <- tax_glom(ps, taxrank = "Genus")
genus_ra <- transform_sample_counts(ps_genus, function(x) x/sum(x))
genus_matrix <- as.data.frame(otu_table(genus_ra))

# 2. CLR变换（处理组成性）
genus_clr <- as.data.frame(clr(t(genus_matrix) + 0.5))  # 加pseudocount避免log(0)

# 3. 过滤低流行度物种（保留在>10%样本中存在的属）
prevalence <- colMeans(genus_matrix > 0)
genus_filtered <- genus_clr[, prevalence > 0.10]

# 4. 逆正态变换（Inverse Normal Transformation）
# 将每个微生物丰度转换为正态分布
int_transform <- function(x) {
  qnorm((rank(x, na.last = "keep") - 0.5) / sum(!is.na(x)))
}
genus_int <- as.data.frame(apply(genus_filtered, 2, int_transform))

# 5. Alpha多样性作为表型
alpha_div <- estimate_richness(ps, measures = c("Shannon", "Simpson", "Observed"))
alpha_int <- as.data.frame(apply(alpha_div, 2, int_transform))

# 6. Beta多样性PCoA坐标
bray_dist <- distance(ps_genus, method = "bray")
pcoa_res <- ordinate(ps_genus, method = "PCoA", distance = bray_dist)
pcoa_axes <- as.data.frame(pcoa_res$vectors[, 1:10])  # 前10个PCoA轴

# 7. 合并所有表型
phenotypes <- cbind(genus_int, alpha_int, pcoa_axes)
# 保存为PLINK格式
write.table(phenotypes, "microbiome_phenotypes.txt",
            sep = "\t", quote = FALSE, row.names = TRUE)

# === GWAS数据质控（PLINK2）===

# 基本QC
plink2 --bfile raw_genotypes \
    --mind 0.05 \
    --geno 0.02 \
    --maf 0.01 \
    --hwe 1e-6 \
    --make-bed --out qc_step1

# 去除亲缘关系过近的个体
plink2 --bfile qc_step1 \
    --king-cutoff 0.0884 \
    --make-bed --out qc_step2

# PCA计算（用于人群分层校正）
plink2 --bfile qc_step2 \
    --pca 20 \
    --out pca_results

# Imputation前准备（phasing + imputation）
# 使用Michigan Imputation Server或TOPMed Server
# 或本地使用 SHAPEIT5 + IMPUTE5/Minimac4

# Imputation后QC
plink2 --pfile imputed_data \
    --maf 0.01 \
    --mac 20 \
    --geno 0.05 \
    --hwe 1e-6 \
    --make-pgen --out imputed_qc

# === PLINK2 线性关联分析 ===

# 准备协变量文件（PCA + 环境因素）
# covar.txt: FID IID PC1 PC2 ... PC10 age sex BMI

# 对单个微生物表型运行GWAS
plink2 --pfile imputed_qc \
    --pheno microbiome_phenotypes.txt \
    --pheno-name Bifidobacterium \
    --covar covariates.txt \
    --covar-name PC1-PC10,age,sex,BMI \
    --glm hide-covar \
    --threads 16 \
    --out gwas_Bifidobacterium

# 批量运行所有微生物表型
for PHENO in $(head -1 microbiome_phenotypes.txt | tr '\t' '\n' | tail -n+2); do
    plink2 --pfile imputed_qc \
        --pheno microbiome_phenotypes.txt \
        --pheno-name ${PHENO} \
        --covar covariates.txt \
        --covar-name PC1-PC10,age,sex,BMI \
        --glm hide-covar \
        --threads 16 \
        --out gwas_${PHENO}
done

# === 使用BOLT-LMM（混合线性模型）===
# BOLT-LMM比普通线性回归更好地控制人群分层

# 命令行
# bolt --bfile=qc_data --phenoFile=pheno.txt --phenoCol=Bifidobacterium \
#      --covarFile=covar.txt --qCovarCol=PC{1:10} --qCovarCol=age --qCovarCol=BMI \
#      --covarCol=sex --lmmForceNonInf --LDscoresFile=LDSCORE.1000G \
#      --statsFile=bolt_results.txt --numThreads=16

# === 使用MaAsLin2/3 进行微生物组关联（非GWAS，用于验证）===
# MaAsLin3（2025年发布）相比MaAsLin2的改进：
# - 同时检测丰度(abundance)和流行度(prevalence)关联
# - 可校正组成性(compositionality)——通过定量PCR或spike-in或计算方法
# - 在IBD等数据集中，77%的关联来自流行度而非丰度
library(Maaslin2)  # 或 library(MaAsLin3) 使用新版

fit_data <- Maaslin2(
  input_data = genus_matrix,
  input_metadata = metadata,  # 包含遗传PRS或候选SNP基因型
  output = "maaslin2_output",
  fixed_effects = c("PRS_score", "age", "sex", "BMI"),
  normalization = "CLR",
  transform = "NONE",
  analysis_method = "LM"
)

# === 结果整理与校正 ===

# 读取所有表型的GWAS结果
library(data.table)
library(qqman)

# 全基因组显著性阈值
# 方法1: Bonferroni for multiple phenotypes
n_phenotypes <- 200  # 假设200个微生物表型
bonf_threshold <- 5e-8 / n_phenotypes  # ~2.5e-10

# 方法2: 有效独立表型数（考虑表型间相关性）
# 使用eigenvalue分解估计有效检验数
pheno_cor <- cor(phenotypes, use = "pairwise.complete")
eigen_vals <- eigen(pheno_cor)$values
n_effective <- sum(eigen_vals > 1)  # 或其他阈值
adjusted_threshold <- 5e-8 / n_effective

# 读取和合并结果
results_list <- list()
for (pheno in phenotype_names) {
  res <- fread(paste0("gwas_", pheno, ".Bifidobacterium.glm.linear"))
  res$phenotype <- pheno
  results_list[[pheno]] <- res[res$P < 1e-5, ]  # 只保留suggestive
}
all_results <- rbindlist(results_list)

# Manhattan plot
manhattan(all_results[all_results$phenotype == "Bifidobacterium", ],
          chr = "CHROM", bp = "POS", snp = "ID", p = "P",
          suggestiveline = -log10(1e-5),
          genomewideline = -log10(5e-8))

# QQ plot
qq(all_results[all_results$phenotype == "Bifidobacterium", ]$P)

# === METAL 荟萃分析 ===
# 准备METAL脚本
cat > metal_script.txt << 'EOF'
SCHEME STDERR
AVERAGEFREQ ON
MINMAXFREQ ON

# Cohort 1
MARKER ID
ALLELE A1 A2
FREQ A1_FREQ
EFFECT BETA
STDERR SE
PVALUE P
PROCESS cohort1_gwas.txt

# Cohort 2
MARKER ID
ALLELE A1 A2
FREQ A1_FREQ
EFFECT BETA
STDERR SE
PVALUE P
PROCESS cohort2_gwas.txt

# Cohort 3
PROCESS cohort3_gwas.txt

OUTFILE meta_Bifidobacterium_ .txt
ANALYZE HETEROGENEITY
QUIT
EOF

metal metal_script.txt

# === R中meta分析 ===
library(meta)

# 对显著位点做森林图
snp_data <- data.frame(
  cohort = c("Dutch", "German", "UK", "Finnish"),
  beta = c(0.15, 0.12, 0.18, 0.10),
  se = c(0.03, 0.04, 0.03, 0.05)
)

meta_result <- metagen(TE = snp_data$beta, seTE = snp_data$se,
                        studlab = snp_data$cohort,
                        sm = "MD", random = TRUE)
forest(meta_result, main = "rs4988235 - Bifidobacterium association")

# === 双向孟德尔随机化 ===
library(TwoSampleMR)

# 方向1：微生物→疾病
# exposure: 微生物丰度的遗传工具变量
exposure_dat <- read_exposure_data(
  filename = "microbiome_gwas_instruments.txt",
  sep = "\t",
  snp_col = "SNP", beta_col = "BETA", se_col = "SE",
  eaf_col = "EAF", effect_allele_col = "A1",
  other_allele_col = "A2", pval_col = "P"
)

# outcome: 疾病GWAS结果
outcome_dat <- extract_outcome_data(
  snps = exposure_dat$SNP,
  outcomes = "ieu-b-2"  # IBD GWAS
)

# harmonize
dat <- harmonise_data(exposure_dat, outcome_dat)

# MR分析
mr_results <- mr(dat, method_list = c("mr_ivw", "mr_egger_regression",
                                       "mr_weighted_median"))
mr_scatter_plot(mr_results, dat)

# 方向2：疾病→微生物
# 使用疾病GWAS位点作为IV，微生物丰度GWAS作为outcome
# 比较两方向结果判断因果

# 标准mbQTL流程
# 1. 基因分型QC + imputation
plink2 --bfile raw --mind 0.05 --geno 0.02 --maf 0.01 --hwe 1e-6 --make-bed --out qc
# 2. PCA
plink2 --bfile qc --pca 20 --out pca
# 3. 微生物表型准备（R中CLR+INT变换）
# 4. GWAS
plink2 --pfile imputed --pheno pheno.txt --covar covar.txt --glm --out results
# 5. Meta-analysis
metal metal_script.txt