170_群体宏基因组分析¶

# === 使用统一流程处理所有队列的原始数据 ===
# 推荐: MetaPhlAn4 + HUMAnN3

for cohort in Study_A Study_B Study_C Study_D; do
  for sample in $(ls ${cohort}/fastq/*.fq.gz | sed 's/_[12].fq.gz//' | sort -u); do
    # MetaPhlAn4物种分类
    metaphlan ${sample}_1.fq.gz,${sample}_2.fq.gz \
      --bowtie2out ${sample}.bt2.bz2 \
      --input_type fastq \
      --nproc 8 \
      -o ${cohort}/profiles/$(basename ${sample})_profile.txt

    # HUMAnN3功能分析
    humann --input ${sample}_1.fq.gz \
      --output ${cohort}/humann/ \
      --threads 8 \
      --metaphlan-options "--bowtie2db /db/metaphlan"
  done

  # 合并该队列的所有profile
  merge_metaphlan_tables.py ${cohort}/profiles/*_profile.txt \
    > ${cohort}/merged_abundance.tsv
done

# === 使用curatedMetagenomicData获取已处理数据 ===
library(curatedMetagenomicData)

# 查询可用的T2D数据集
available_studies <- sampleMetadata %>%
  filter(disease == "T2D") %>%
  select(study_name, body_site, number_reads) %>%
  distinct(study_name)

# 下载多个队列
datasets <- curatedMetagenomicData(
  c("QinJ_2012.relative_abundance",
    "KarlssonFH_2013.relative_abundance",
    "LiJ_2017.relative_abundance"),
  dryrun = FALSE
)

# 合并为统一格式
library(TreeSummarizedExperiment)
merged_data <- mergeData(datasets)

library(vegan)
library(ggplot2)

# === 诊断批次效应 ===
# PCoA/PERMANOVA检测研究是否是主要变异来源
abundance_matrix <- assay(merged_data)  # taxa × samples
metadata <- colData(merged_data)

# Bray-Curtis距离
bc_dist <- vegdist(t(abundance_matrix), method = "bray")

# PERMANOVA: 疾病 vs 研究
adonis_result <- adonis2(bc_dist ~ disease + study_name + age + BMI,
                          data = as.data.frame(metadata),
                          permutations = 999)
print(adonis_result)
# 如果study_name解释的方差远大于disease → 批次效应严重

# PCoA可视化
pcoa_res <- cmdscale(bc_dist, k = 2)
plot_df <- data.frame(PC1 = pcoa_res[,1], PC2 = pcoa_res[,2],
                       Study = metadata$study_name,
                       Disease = metadata$disease)
ggplot(plot_df, aes(x = PC1, y = PC2, color = Study, shape = Disease)) +
  geom_point(alpha = 0.7) +
  theme_minimal() +
  labs(title = "PCoA - Before Batch Correction")

# === 批次校正方法 ===
# 方法1: ComBat (适合连续数据如CLR转换后)
library(sva)
abundance_clr <- compositions::clr(t(abundance_matrix + 0.5))
batch <- metadata$study_name
condition <- metadata$disease

combat_corrected <- ComBat(
  dat = t(abundance_clr),
  batch = batch,
  mod = model.matrix(~ condition)
)

# 方法2: MMUPHin (专为微生物组meta-analysis)
library(MMUPHin)
fit_adjust <- adjust_batch(
  feature_abd = abundance_matrix,
  batch = "study_name",
  covariates = "disease",
  data = as.data.frame(metadata)
)
adjusted_abundance <- fit_adjust$feature_abd_adj

# 方法3: 不做批次校正, 在每个研究内分析后合并效应量
# (最保守, 避免校正引入artifact)

library(meta)
library(metafor)

# === 按研究分别做差异丰度分析 ===
# 对每个研究独立运行DA分析
study_results <- list()

for (study in unique(metadata$study_name)) {
  study_data <- abundance_matrix[, metadata$study_name == study]
  study_meta <- metadata[metadata$study_name == study, ]

  # 使用Wilcoxon/ALDEx2/等做DA分析
  # 提取每个taxa的效应量(log2FC)和标准误
  results <- run_da_analysis(study_data, study_meta$disease)
  results$study <- study
  study_results[[study]] <- results
}

all_results <- do.call(rbind, study_results)

# === Random-effects meta-analysis ===
taxa_list <- unique(all_results$taxon)
meta_results <- data.frame()

for (taxon in taxa_list) {
  taxon_data <- all_results[all_results$taxon == taxon, ]

  if (nrow(taxon_data) >= 3) {  # 至少3个研究
    # Random-effects model
    meta_fit <- rma(yi = taxon_data$log2FC,
                     sei = taxon_data$SE,
                     method = "REML")

    meta_results <- rbind(meta_results, data.frame(
      taxon = taxon,
      overall_effect = meta_fit$beta[1],
      se = meta_fit$se,
      pvalue = meta_fit$pval,
      I2 = meta_fit$I2,      # 异质性
      n_studies = nrow(taxon_data)
    ))
  }
}

# FDR校正
meta_results$qvalue <- p.adjust(meta_results$pvalue, method = "BH")
sig_taxa <- meta_results[meta_results$qvalue < 0.05, ]
print(paste("显著差异taxa:", nrow(sig_taxa)))

# === 森林图可视化 ===
# 对top taxa绘制森林图
top_taxon <- sig_taxa$taxon[1]
taxon_data <- all_results[all_results$taxon == top_taxon, ]

forest_meta <- metagen(
  TE = taxon_data$log2FC,
  seTE = taxon_data$SE,
  studlab = taxon_data$study,
  sm = "MD"
)
forest(forest_meta, main = top_taxon)

# === Leave-One-Study-Out验证 ===
# 每次留出一个研究作为验证集，其余做训练
studies <- unique(metadata$study_name)
loso_results <- list()

for (held_out in studies) {
  # 训练集
  train_idx <- metadata$study_name != held_out
  train_data <- abundance_matrix[, train_idx]
  train_labels <- metadata$disease[train_idx]

  # 测试集
  test_idx <- metadata$study_name == held_out
  test_data <- abundance_matrix[, test_idx]
  test_labels <- metadata$disease[test_idx]

  # 训练分类器
  library(randomForest)
  rf_model <- randomForest(x = t(train_data), y = factor(train_labels), ntree = 500)

  # 预测
  predictions <- predict(rf_model, t(test_data), type = "prob")[, "T2D"]
  auc <- pROC::auc(test_labels, predictions)

  loso_results[[held_out]] <- data.frame(
    study = held_out,
    AUC = as.numeric(auc),
    n_test = sum(test_idx)
  )
}

loso_df <- do.call(rbind, loso_results)
cat(sprintf("LOSO平均AUC: %.3f ± %.3f\n", mean(loso_df$AUC), sd(loso_df$AUC)))

# === 群体宏基因组meta-analysis流程 ===
# 1. 统一处理所有队列
for cohort in A B C D; do
  metaphlan --input_type fastq -o ${cohort}_profile.txt ${cohort}_reads.fq.gz
done

# 2. R分析
Rscript meta_analysis.R --studies "A,B,C,D" --output results/