是否想知道你的基因表达是否真的直接导致了表型的改变；——中介分析；宏基因组、转录组、蛋白质组、代谢组的均适用

比如这篇文章的结论：

图C解读（E. rectale 的 InDel 变异，可以理解成基因突变或表达情况）

各元素解读

左侧：自变量（X）—— E. rectale 的 InDel 变异

指直肠真杆菌中一个 4bp 的插入变异，该变异会导致维生素 B12 ABC 转运蛋白 BtuD 提前截短（从 285 个氨基酸缩短至 24 个氨基酸），直接破坏 B12 的转运功能（B12 是微生物和宿主神经递质合成的关键物质）。
图中 “a = -1.09; p=0.004”：a 路径（X→M）的标准化回归系数和 p 值，表明该 InDel 变异与咖啡因水平显著相关（负向效应）。

中间：中介变量（M）—— 咖啡因（Caffeine）

咖啡因是具有神经保护作用的代谢物，文章提到其在 ASD 患儿粪便中呈显著 depletion（下调）。
该 InDel 变异通过破坏 B12 转运，影响微生物对神经活性化合物的降解 / 转化，最终导致咖啡因水平改变。

右侧：因变量（Y）—— ASD 表型（自闭症）

图中 “beta = -1.35; p=0.03”：b 路径（M→Y）的回归系数和 p 值，表明咖啡因水平与 ASD 状态显著相关（负向效应，即咖啡因水平越低，ASD 风险越高）。
中心 “28.25%”：中介效应占总效应的比例，意味着该 InDel 变异对 ASD 的影响中，有 28.25% 是通过调控咖啡因代谢实现的。

中介分析（Mediation Analysis）是统计学与因果推断领域中一种重要的量化分析方法，核心目的是揭示变量间“原因→中介→结果”的间接作用机制，明确自变量（预测变量）对因变量（结果变量）的影响，究竟是直接产生的，还是通过一个或多个中间变量（中介变量）间接传递的。它打破了传统“自变量→因变量”的简单线性认知，能够更精细地拆解因果路径，是探索复杂变量关系、挖掘内在作用机理的核心工具之一。

一、核心逻辑与核心问题

中介分析的核心逻辑基于因果链条的拆解：当我们观察到自变量（X）与因变量（Y）存在关联时，并非所有关联都是X对Y的直接作用。可能存在这样的路径：X先影响某个中间变量（M），再由M进一步影响Y。此时，X对Y的总效应（Total Effect）可分解为两部分：一是X直接作用于Y的直接效应（Direct Effect）；二是X通过M间接作用于Y的中介效应（Mediated Effect，也常称为间接效应）。

其核心要回答的问题包括：① 自变量对因变量的影响是否存在中介效应？② 中介效应的大小与方向如何？③ 中介效应占总效应的比例是多少（即中介效应的解释力）？④ 该中介路径是否具有统计显著性？

二、关键要素

完整的中介分析需明确三个核心变量，三者构成清晰的因果路径链条：

自变量（X，也称为处理变量/暴露变量）：引发后续效应的“原因”变量，是分析的起点。例如，在医学研究中可能是“药物干预”“基因表达水平”，在社会科学中可能是“教育投入”“政策实施”。
中介变量（M）：连接自变量与因变量的“桥梁”变量，是自变量影响因变量的中间载体。它会受到X的影响，同时又会影响Y。例如，“炎症因子”可能是“药物干预”与“疾病康复”之间的中介变量，“学习动机”可能是“教育投入”与“学业成绩”之间的中介变量。
因变量（Y，也称为结局变量）：最终被影响的“结果”变量，是分析的终点。例如，“疾病严重程度”“治愈率”“收入水平”等。

此外，分析中通常还会纳入协变量（Covariates）（如年龄、性别、基线状态等），用于控制混杂因素的干扰，确保中介效应估计的准确性。

三、核心分析步骤（经典逐步法框架）

传统中介分析采用逐步回归的思路验证中介效应，虽现代方法更注重直接检验中介效应的显著性，但逐步法仍是理解中介逻辑的基础，核心步骤如下：

检验总效应：构建回归模型Y = f(X + 协变量)，验证自变量X对因变量Y的总效应是否显著。这是中介分析的前提——若总效应不显著，后续中介效应分析的意义通常有限（特殊情况除外，如存在抑制性中介）。
检验X对M的效应（a路径）：构建回归模型M = f(X + 协变量)，验证自变量X对中介变量M的影响是否显著。这一步对应“原因→中介”的路径，是中介效应存在的核心前提之一。
检验M对Y的效应（b路径）：构建回归模型Y = f(X + M + 协变量)，在控制X的前提下，验证中介变量M对因变量Y的影响是否显著。这一步对应“中介→结果”的路径。
判断中介效应类型：观察第三步模型中X对Y的系数（直接效应，c’路径）：① 若c’不显著，且a、b路径均显著，说明存在“完全中介效应”（X对Y的影响完全通过M传递）；② 若c’显著，且a、b路径均显著，说明存在“部分中介效应”（X对Y的影响部分通过M传递，部分直接产生）。
中介效应的量化与显著性检验：中介效应的大小通常为a路径系数与b路径系数的乘积（a×b）。现代分析中，更常用Bootstrap法、Sobel法等直接检验a×b的显著性——若置信区间不包含0，则中介效应显著。

代码：

日常好的代码已放入免💰共享服务器中（人人皆可用）：https://vip.r-py.com/

今日的代码日期是code260115

# 设置工作目录，定位到中介分析的数据文件夹setwd('/home/data/CWN/ASD/mediation')# 获取命令行传入的参数（用于批量处理，指定循环的起始、结束位置和输出文件编号）args <- commandArgs(trailingOnly = TRUE)first <- args[1]    # 暴露变量循环起始索引last <- args[2]     # 暴露变量循环结束索引file_num <- args[3] # 输出文件的编号后缀# 读取合并后的差异分析表格数据# 参数说明：header=T保留表头，sep=","按逗号分隔，comment.char="#"忽略注释行，row.names=NULL不使用行名，check.names=F不检查列名合法性final <- read.table('WGS/2-merge/3.mediation/1.diff_table_combine.csv',                    header = T,                     sep = ",",                     comment.char = "#",                    row.names = NULL,                    check.names = F)# 定义变量分组：mediator <- colnames(final)[10:108]    # 第10-108列为中介变量（微生物特征）treators <- colnames(final)[109:2181]  # 第109-2181列为暴露变量（表达量数据）covars <- c("Age","Gender")            # 定义协变量：年龄和性别# 提取协变量数据并处理性别变量（将字符型转换为数值型：F=1，M=0）covar_df <- final %>% dplyr::select(row.names, all_of(covars))covar_df$Gender <- as.integer(sub( "F", 1, sub("M",0,covar_df$Gender)))# ===================== 核心：中介分析 =====================# 初始化存储中介分析结果的空对象Mediation.forward <- NULL# 提取结局变量（Group：ASD/TD）并转换为数值型：ASD=1，TD=0meta.lf = final %>% dplyr::select(row.names, all_of("Group"))meta.lf$Group<- as.integer(sub( "ASD", 1, sub("TD",0,meta.lf$Group)))lf='Group'  # 定义结局变量名称# 提取所有中介变量数据mf.data <- final %>% dplyr::select(row.names, all_of(mediator))# 安全校验：如果结束索引超过暴露变量总数（2037），则重置为2037if(as.integer(last)>2037){  last=2037}# 外层循环：遍历指定范围的暴露变量（从first到last）for(exp in treators[as.integer(first):as.integer(last)]){  # 提取当前暴露变量的数据  meta.expo <- final %>% dplyr::select(row.names, all_of(exp))
  # 内层循环：遍历每个中介变量（从第2列开始，第1列是row.names）  for(i in 2:ncol(mf.data)){    # 提取当前中介变量的临时数据    mf.dat.tmp <- mf.data[,c(1,i), drop=F]    mf = colnames(mf.data)[i]  # 获取当前中介变量名称
    # 数据合并：将暴露变量、结局变量、协变量、中介变量按行名合并    dat <- left_join(left_join(left_join(meta.expo, meta.lf, by = "row.names"),                       covar_df, by="row.names"),                 mf.dat.tmp, by="row.names")
    # 移除含有缺失值的行，保证后续建模数据完整    dat <- dat[complete.cases(dat),]
    # -------------------- 向前中介分析 --------------------    # 生成唯一标识：暴露变量|中介变量|结局变量    id = paste(exp, mf, lf, sep = "|")
    # 1. 模型m：中介变量 ~ 暴露变量 + 协变量（线性回归）    # 构建公式：mf ~ exp + 协变量（排除row.names、lf、mf、exp列）    model.m = lm(as.formula( paste(mf, " ~ ", exp, " + ",                                    paste(colnames(dat)[!(colnames(dat) %in% c("row.names",lf, mf, exp))], collapse = " + "),                                   sep = "") ), data = dat)    # 计算标准化beta系数（a路径：暴露对中介的效应）和p值    m_beta <- summary(model.m)$coefficients[2,1] * (sd(dat[,exp])/sd(dat[,mf]))    m_p <- summary(model.m)$coefficients[2,4]
    # 2. 模型y：结局变量 ~ 暴露变量 + 中介变量 + 协变量（二项逻辑回归，因结局是分类变量）    model.y = glm(as.formula( paste(lf, " ~ ", exp," + ", mf," + ",                                    paste(colnames(dat)[!(colnames(dat) %in% c("row.names",lf, mf, exp))], collapse = " + "),                                   sep = "") ), data = dat, family = binomial)    # 计算各路径标准化beta系数和p值：    y_beta <- summary(model.y)$coefficients[2,1] * (sd(dat[,exp])/sd(dat[,lf]))  # c'路径：暴露对结局的直接效应    y_p <- summary(model.y)$coefficients[2,4]    my_beta <- summary(model.y)$coefficients[3,1]* (sd(dat[,mf])/sd(dat[,lf]))   # b路径：中介对结局的效应    my_p <- summary(model.y)$coefficients[3,4]
    # 3. 使用mediate包进行中介效应检验（非参数bootstrap，1000次模拟）    summary = summary(mediate(model.m,model.y,treat=exp,mediator=mf,boot=F,sims=1000))    # 捕获mediate输出结果（用于后续提取ACME/ADE等指标）    res <- capture.output(summary,append=FALSE)
    # 4. 提取ACME（平均因果中介效应）及其95%CI、p值    tmp <-  base::strsplit(res[grep("ACME",res)],"\\s")[[3]]  # 按空格分割ACME行    tmp <- tmp[tmp != "" & tmp!="."]                         # 移除空值和点    tmp <- tmp[!grepl("*",tmp,fixed = T)]                    # 移除显著性星号    ACME.p <- tmp[length(tmp)]                               # ACME的p值    ACME <- tmp[3]                                           # ACME值    ACME_low <- tmp[4]                                       # ACME 95%CI下限    ACME_high <- tmp[5]                                      # ACME 95%CI上限
    # 5. 提取ADE（平均直接效应）及其95%CI、p值    tmp <- base::strsplit(res[grep("ADE",res)],"\\s")[[3]]    tmp <- tmp[tmp != "" & tmp!="."]    tmp <- tmp[!grepl("*",tmp,fixed = T) ]    ADE.p <- tmp[length(tmp)]    ADE <- tmp[3]    ADE_low <- tmp[4]    ADE_high <- tmp[5]
    # 6. 提取中介效应比例（Prop. Mediated）及其p值    tmp <- base::strsplit(res[grep("Prop. Mediated", res)],"\\s")[[3]]    tmp <- tmp[tmp != "" ]    i_str = which(grepl("Mediated", tmp))    prop.mediated <- tmp[(i_str + 2)]    # 中介效应比例    prop.mediated.p <- tmp[(i_str + 5)]  # 中介效应比例的p值
    # 7. 整合当前组合（暴露-中介-结局）的所有结果    forw_vec = c(id, m_beta, m_p, my_beta, my_p,y_beta, y_p,                  ACME,ACME_low, ACME_high, ACME.p,                  ADE, ADE_low, ADE_high, ADE.p,                  prop.mediated,prop.mediated.p)    # 为结果向量命名（方便后续合并和识别）    names(forw_vec) <- c("Treat_Mediator_Y","a-Treater_Mediator_beta","a-Treater_Mediator_p",                         "b-Mediator_y_beta","b-Mediator_y_p","c'-Treater_y_beta","c'-Treater_y_p",                         "ACME","ACME(95%CI low)","ACME(95%CI high)", "ACME.p",                          "ADE", "ADE(95%CI low)", "ADE(95%CI high)","ADE.p",                          "prop.mediated","prop.mediated.p")
    # 将当前结果添加到总结果集中（按行合并）    Mediation.forward <- bind_rows(Mediation.forward, forw_vec)  }# 内层循环结束：单个微生物特征遍历完成}# 外层循环结束：指定范围的暴露变量遍历完成# 将最终中介分析结果写入CSV文件（按指定编号命名）# 参数说明：row.names=F不输出行名，col.names=TRUE输出列名，sep=","逗号分隔，quote=FALSE不添加引号，na=''缺失值为空write.table(Mediation.forward,            file = paste0("3.mediation_",file_num,'.csv'),            row.names = F,            col.names = TRUE,             sep = ',',             quote = FALSE,             na = '')