624_微生物组Meta分析¶

策略1：基于摘要统计量的Meta分析（传统）
- 从每个研究提取效应量（如log fold change、标准化均差）
- 用Meta分析统计方法合并效应量
- 优点：不需要原始数据，文献即可做
- 缺点：受限于各研究报告的信息

策略2：基于原始数据的跨研究分析（推荐）
- 获取各研究的原始OTU/ASV表
- 统一处理后合并分析
- 优点：可以探索更多问题，控制更多混杂
- 缺点：需要原始数据，批次效应处理困难

# === 准备上传到MicrobiomeAnalyst的数据 ===
# 网址：https://www.microbiomeanalyst.ca/

# 需要三个文件：
# 1. OTU/ASV丰度表（行=特征，列=样本）
# 2. 分类注释文件（特征ID到分类的映射）
# 3. 样本元数据（分组、临床信息等）

# === 在R中准备格式 ===
library(phyloseq)  # 加载phyloseq

# 假设已有phyloseq对象
# 导出OTU表
otu_table_export <- as.data.frame(otu_table(ps))  # 提取OTU表
write.csv(otu_table_export, "otu_table.csv")  # 保存

# 导出分类表
tax_table_export <- as.data.frame(tax_table(ps))  # 提取分类表
write.csv(tax_table_export, "taxonomy.csv")  # 保存

# 导出元数据
sample_data_export <- as.data.frame(sample_data(ps))  # 提取元数据
write.csv(sample_data_export, "metadata.csv")  # 保存

# MicrobiomeAnalyst 2.0 主要模块：
# 
# 1. MDP (Marker Data Profiling)
#    - 16S/ITS标记基因分析
#    - Alpha/Beta多样性
#    - 差异分析（DESeq2, LEfSe, ANCOM等）
#    - 功能预测（PICRUSt2）
#
# 2. SDP (Shotgun Data Profiling)
#    - 宏基因组数据分析
#    - 功能通路分析
#    - 代谢网络可视化
#
# 3. PPD (Patterns & Predictions)
#    - 多研究Meta分析
#    - 机器学习分类
#    - 生物标志物发现
#
# 4. MDEA (Meta-Data Enrichment Analysis)
#    - 微生物组富集分析
#    - 跨数据集比较

# === 安装必要R包 ===
install.packages("meta")           # Meta分析核心包
install.packages("metafor")        # 更灵活的Meta分析
install.packages("dmetar")         # Meta分析辅助工具
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("curatedMetagenomicData")  # 公共宏基因组数据

library(meta)      # 加载meta
library(metafor)   # 加载metafor
library(ggplot2)   # 加载ggplot2

# === 准备各研究的效应量数据 ===
# 假设从多个研究中提取了某细菌属在疾病vs对照中的差异
study_data <- data.frame(
    study = c("Study1_2022", "Study2_2023", "Study3_2024",
              "Study4_2024", "Study5_2025"),
    n_case = c(50, 80, 120, 65, 200),     # 病例组样本量
    n_control = c(50, 75, 110, 70, 195),   # 对照组样本量
    mean_case = c(0.15, 0.18, 0.12, 0.20, 0.14),    # 病例组均值（相对丰度）
    sd_case = c(0.08, 0.10, 0.06, 0.12, 0.07),      # 病例组标准差
    mean_control = c(0.05, 0.08, 0.04, 0.09, 0.06), # 对照组均值
    sd_control = c(0.03, 0.05, 0.03, 0.06, 0.04)    # 对照组标准差
)

# === 计算标准化均差（SMD）===
m <- metacont(
    n.e = n_case,           # 实验组（病例）样本量
    mean.e = mean_case,     # 实验组均值
    sd.e = sd_case,         # 实验组标准差
    n.c = n_control,        # 对照组样本量
    mean.c = mean_control,  # 对照组均值
    sd.c = sd_control,      # 对照组标准差
    studlab = study,        # 研究标签
    data = study_data,      # 数据框
    sm = "SMD",             # 效应量类型：标准化均差
    method.smd = "Hedges",  # 用Hedges' g（小样本校正）
    random = TRUE,          # 随机效应模型
    fixed = FALSE           # 不用固定效应模型
)

# === 输出结果 ===
print(summary(m))  # 打印摘要

# 关键结果解读：
# SMD > 0：病例组丰度高于对照组
# p < 0.05：差异显著
# I² > 50%：研究间异质性较高

# === 画森林图 ===
forest(m,
    sortvar = TE,                    # 按效应量排序
    prediction = TRUE,               # 显示预测区间
    print.tau2 = TRUE,               # 显示τ²（研究间方差）
    leftcols = c("studlab", "n.e", "n.c"),  # 左侧显示列
    leftlabs = c("Study", "Cases", "Controls"),  # 左侧标签
    rightcols = c("effect", "ci"),   # 右侧显示列
    rightlabs = c("SMD", "95% CI"),  # 右侧标签
    col.diamond = "steelblue",       # 合并效应菱形颜色
    col.square = "darkblue"          # 各研究方块颜色
)

# === 异质性检验 ===
cat("Cochran Q检验 p值:", m$pval.Q, "\n")  # Q检验
cat("I² 统计量:", m$I2, "%\n")              # I²
cat("τ² (研究间方差):", m$tau2, "\n")        # τ²

# I²解读：
# <25%: 低异质性
# 25-50%: 中等异质性
# 50-75%: 较高异质性
# >75%: 高异质性

# === 漏斗图检查发表偏倚 ===
funnel(m,
    xlab = "Standardized Mean Difference",  # x轴标签
    studlab = TRUE                           # 显示研究标签
)

# === Egger检验（发表偏倚统计检验）===
metabias(m, method.bias = "Egger")  # Egger线性回归检验
# p < 0.05 提示可能存在发表偏倚

# === 使用curatedMetagenomicData获取公共数据 ===
library(curatedMetagenomicData)  # 加载包

# 查看可用数据集
available_datasets <- sampleMetadata  # 所有可用样本元数据
table(available_datasets$study_condition)  # 查看疾病类型

# 获取特定疾病的所有研究数据
t2d_data <- curatedMetagenomicData(
    "*.relative_abundance",  # 相对丰度数据
    dryrun = FALSE,          # 实际下载
    counts = TRUE            # 返回计数
) |> mergeData()             # 合并数据

# 过滤T2D相关样本
t2d_samples <- sampleMetadata |>
    dplyr::filter(study_condition %in% c("T2D", "control"))  # T2D和对照

# === 使用MMUPHin进行批次效应校正 ===
# MMUPHin专门为微生物组Meta分析设计
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("MMUPHin")  # 安装MMUPHin

library(MMUPHin)  # 加载MMUPHin

# 批次效应校正（对丰度表）
adjusted_abundance <- adjust_batch(
    feature_abd = abundance_matrix,  # 特征丰度矩阵
    batch = "study",                 # 批次变量（研究来源）
    covariates = "disease_status",   # 要保留的生物学变量
    data = metadata                  # 样本元数据
)

# 校正后的差异分析
lm_results <- lm_meta(
    feature_abd = adjusted_abundance$feature_abd_adj,  # 校正后丰度
    batch = "study",           # 批次变量
    exposure = "disease_status",  # 暴露/分组变量
    data = metadata,           # 元数据
    control = list(verbose = FALSE)
)

# 查看显著关联的微生物
sig_features <- lm_results$meta_fits |>
    dplyr::filter(qval.fdr < 0.05) |>  # FDR < 0.05
    dplyr::arrange(qval.fdr)            # 按q值排序

print(sig_features)  # 打印显著特征

# === Python跨疾病Meta分析 ===
import pandas as pd      # 数据处理
import numpy as np       # 数值计算
from scipy import stats  # 统计检验
import matplotlib.pyplot as plt  # 画图
import seaborn as sns    # 高级可视化

# 读取多个研究的差异分析结果
studies = {
    'IBD_study1': pd.read_csv('ibd_study1_deseq2.csv'),
    'IBD_study2': pd.read_csv('ibd_study2_deseq2.csv'),
    'T2D_study1': pd.read_csv('t2d_study1_deseq2.csv'),
    'T2D_study2': pd.read_csv('t2d_study2_deseq2.csv'),
    'CRC_study1': pd.read_csv('crc_study1_deseq2.csv'),
}

# 提取共同特征的效应量
def extract_effect_sizes(studies, feature):
    """提取某个微生物特征在所有研究中的效应量"""
    results = []  # 存储结果
    for study_name, df in studies.items():
        row = df[df['feature'] == feature]  # 查找特征
        if len(row) > 0:  # 如果找到了
            results.append({
                'study': study_name,
                'log2fc': row['log2FoldChange'].values[0],  # 效应量
                'se': row['lfcSE'].values[0],                # 标准误
                'pvalue': row['pvalue'].values[0]             # p值
            })
    return pd.DataFrame(results)  # 返回数据框

# 跨研究合并效应量（随机效应模型）
def random_effects_meta(effects, ses):
    """DerSimonian-Laird随机效应Meta分析"""
    weights_fixed = 1 / (ses ** 2)  # 固定效应权重

    # 计算Q统计量
    theta_fixed = np.sum(weights_fixed * effects) / np.sum(weights_fixed)
    Q = np.sum(weights_fixed * (effects - theta_fixed) ** 2)

    # 计算τ²
    k = len(effects)  # 研究数
    c = np.sum(weights_fixed) - np.sum(weights_fixed ** 2) / np.sum(weights_fixed)
    tau2 = max(0, (Q - (k - 1)) / c)  # τ²不能为负

    # 随机效应权重
    weights_random = 1 / (ses ** 2 + tau2)
    theta_random = np.sum(weights_random * effects) / np.sum(weights_random)
    se_random = np.sqrt(1 / np.sum(weights_random))

    # I²
    I2 = max(0, (Q - (k - 1)) / Q * 100)

    # z检验
    z = theta_random / se_random
    p = 2 * (1 - stats.norm.cdf(abs(z)))

    return {
        'effect': theta_random,  # 合并效应量
        'se': se_random,         # 标准误
        'p': p,                  # p值
        'I2': I2,                # I²
        'tau2': tau2             # τ²
    }

# 对每个微生物特征做Meta分析
print("跨研究Meta分析完成")