跳转至

微生物组纵向研究设计

一句话概述

微生物组纵向研究需要使用零膨胀负二项回归(ZINB)、线性混合效应模型(LME)、负二项广义线性混合模型(NB-GLMM)和时间序列分析等专门统计方法，正确处理重复测量、时间相关性和零膨胀等数据特征，是研究微生物群落动态变化的规范分析框架。

核心知识点表格

知识点关键内容常用工具
纵向设计重复测量、时间序列、面板数据-
零膨胀模型ZINB/ZIP处理过多零值ZIBR, pscl, glmmTMB
线性混合模型随机截距/斜率处理个体差异lme4, nlme
负二项GLMM计数数据+随机效应glmmTMB, NBZIMM
时间序列方法DTW/ARIMA/状态空间模型tseries, dtw
差异丰度纵向差异分析MaAsLin2, LinDA, NBZIMM
多变量方法社区水平纵向变化PERMANOVA(adonis2)
因果推断Granger因果/中介分析mediation, vars

各步骤详解

第一步：理解纵向研究设计

白话解释： 纵向研究就是对同一批个体在不同时间点反复采样测量。相比横截面研究（一次性采样），纵向研究能看到微生物组随时间的变化轨迹——是谁在变、怎么变、什么时候变。但统计分析更复杂，因为同一个人的多次测量不是独立的（相关性需要纳入模型）。

技术细节： - 重复测量相关性：同一个体的样本间相关（ICC通常0.3-0.7） - 零膨胀：结构性零（物种确实不存在）+ 抽样零（存在但没检出） - 不平衡设计：不同个体的采样时间点可能不同 - 时间自相关：相邻时间点比远间隔时间点更相关 - 丢失随访(dropout)：部分个体中途退出

第二步：ZIBR零膨胀贝叶斯回归

白话解释： ZIBR (Zero-Inflated Beta Regression with Random effects) 是专门为微生物组纵向数据设计的方法。它用两部分模型处理零值问题：第一部分（logistic）建模物种"存在/不存在"的概率；第二部分（beta回归）建模在存在的样本中的相对丰度。两部分都包含随机效应处理个体差异。

代码示例：

# 安装ZIBR（已在CRAN上，当前版本1.0.2）
# install.packages("ZIBR")  # 或 devtools::install_github("PennChopMicrobiomeProgram/ZIBR")
library(ZIBR)

# 准备数据
# 需要: 长格式数据框
# subject_id, time, abundance(相对丰度0-1), covariates
data <- data.frame(
    subject = rep(1:50, each=5),           # 50个个体
    time = rep(c(0, 1, 3, 6, 12), 50),    # 5个时间点(月)
    abundance = rbeta(250, 0.5, 5),        # 相对丰度
    treatment = rep(c("A","B"), each=125), # 处理组
    age = rnorm(250, 50, 10)
)
# 添加结构性零
data$abundance[sample(1:250, 80)] <- 0

# 对每个目标OTU运行ZIBR
zibr_result <- zibr(
    logistic_cov = data[, c("time", "treatment")],  # 零值模型的协变量
    beta_cov = data[, c("time", "treatment")],       # 丰度模型的协变量
    Y = data$abundance,                               # 响应变量(相对丰度)
    subject_ind = data$subject,                        # 个体ID
    time_ind = data$time                               # 时间变量
)

# 查看结果
# logistic部分: 处理对存在/缺失的影响
cat("Zero-inflation logistic model:\n")
print(zibr_result$logistic_est_table)

# beta部分: 处理对非零丰度的影响
cat("Beta regression model:\n")
print(zibr_result$beta_est_table)

# 联合检验
cat("Joint p-value (treatment effect):", 
    zibr_result$joint_p, "\n")

第三步：线性混合效应模型（LME）

白话解释： 线性混合效应模型是纵向分析的基础。"固定效应"表示所有个体共有的趋势（如治疗效果），"随机效应"表示个体间的差异（如每个人的基线水平和变化速率不同）。对微生物组数据，通常先做CLR转化再使用LME。

代码示例：

library(lme4)
library(lmerTest)
library(ggplot2)

# 数据准备（CLR转化后的丰度）
library(compositions)
# 假设 otu_table: samples × OTUs
otu_clr <- as.data.frame(clr(otu_table + 0.5))

# 长格式数据
library(tidyverse)
long_data <- otu_clr %>%
    mutate(sample_id = rownames(.)) %>%
    pivot_longer(-sample_id, names_to="taxon", values_to="clr_abundance") %>%
    left_join(metadata, by="sample_id")

# 对单个分类单元拟合LME
# 随机截距+随机斜率模型
model <- lmer(clr_abundance ~ time * treatment + age + sex +
              (1 + time | subject_id),   # 随机截距+随机斜率
              data = long_data %>% filter(taxon == "Bacteroides"),
              REML = TRUE)

summary(model)
# 关注: time:treatment 交互项 → 处理组的变化速率是否不同

# 模型诊断
plot(model)  # 残差图
qqnorm(resid(model)); qqline(resid(model))

# 可视化拟合曲线
library(effects)
eff <- effect("time:treatment", model)
plot(eff, multiline=TRUE, ci.style="bands")

# 批量分析所有分类单元
taxa_list <- colnames(otu_clr)
lme_results <- lapply(taxa_list, function(tax) {
    sub_data <- long_data %>% filter(taxon == tax)
    tryCatch({
        m <- lmer(clr_abundance ~ time * treatment + age + sex +
                  (1 | subject_id),
                  data = sub_data)
        coefs <- summary(m)$coefficients
        data.frame(
            taxon = tax,
            effect = rownames(coefs),
            estimate = coefs[, "Estimate"],
            pvalue = coefs[, "Pr(>|t|)"]
        )
    }, error = function(e) NULL)
})

lme_df <- do.call(rbind, lme_results)

# 筛选时间×治疗交互显著的分类单元
sig_interaction <- lme_df %>%
    filter(effect == "time:treatmentB") %>%
    mutate(padj = p.adjust(pvalue, method="BH")) %>%
    filter(padj < 0.05)

第四步：负二项GLMM（直接用计数）

白话解释： 如果不想做CLR转化，可以直接对计数数据使用负二项GLMM。负二项分布天然适合处理微生物计数数据的过度离散(overdispersion)问题，GLMM的随机效应处理重复测量。

代码示例：

library(glmmTMB)
library(NBZIMM)

# 方法1: glmmTMB (灵活，推荐)
# 零膨胀负二项混合模型
model_zinb <- glmmTMB(
    count ~ time * treatment + age + offset(log(total_reads)) +
            (1 + time | subject_id),
    ziformula = ~ time + treatment,  # 零膨胀部分的协变量
    family = nbinom2(),               # 负二项分布
    data = count_data)

summary(model_zinb)

# 方法2: NBZIMM (专门为微生物组设计)
# devtools::install_github("nyiuab/NBZIMM")
library(NBZIMM)

# 使用NBZIMM进行纵向差异丰度分析
# mms可同时分析多个分类单元（y为矩阵，列为不同taxa）
nbzimm_result <- mms(
    y = otu_count_matrix,              # taxa计数矩阵（样本×分类单元）
    fixed = ~ time * treatment + age + offset(log(total_reads)),
    random = ~ 1 | subject_id,
    zi_fixed = ~ 1,
    data = count_data,
    method = "zinb",                    # 可选: "nb", "lme", "zinb", "zig"
    min.p = 0.2)                        # 非零值最小比例

# 方法3: MaAsLin2 (Huttenhower lab，推荐用于微生物组)
library(Maaslin2)

# MaAsLin2直接处理纵向设计
# analysis_method可选: "LM"(线性模型), "CPLM"(复合泊松), "NEGBIN"(负二项), "ZINB"(零膨胀负二项)
# 注意: NEGBIN/ZINB只支持normalization="CSS"/"NONE"/"TMM"
maaslin_result <- Maaslin2(
    input_data = otu_table,           # OTU表 (samples × taxa)
    input_metadata = metadata,         # 元数据
    output = "maaslin2_output/",
    fixed_effects = c("time", "treatment", "age"),
    random_effects = c("subject_id"),  # 随机效应处理重复测量
    normalization = "TSS",
    transform = "AST",                 # arcsin sqrt 转化
    analysis_method = "LM",            # 线性模型（配合TSS+AST最常用）
    min_abundance = 0.001,
    min_prevalence = 0.1,
    cores = 8)

# 查看显著结果
sig_results <- maaslin_result$results %>%
    filter(qval < 0.05)

# 注意：MaAsLin3已发布，支持同时检测丰度和流行率关联
# 可通过 BiocManager::install("biobakery/MaAsLin3") 安装

第五步：时间序列特异方法

白话解释： 经典纵向统计方法假设时间只是一个协变量。但微生物组时间序列还有自身特点：自相关性（今天的组成与昨天相似）、周期性（日周期/月周期）、突变事件（如抗生素使用后的急剧变化和恢复）。专门的时间序列方法能捕获这些模式。

代码示例：

# 1. 动态时间规整（DTW）- 比较时间序列相似性
library(dtw)

# 比较两个个体的微生物组变化轨迹
ts_subject1 <- diversity_over_time[subject == "S001", "shannon"]
ts_subject2 <- diversity_over_time[subject == "S002", "shannon"]

# DTW距离
dtw_result <- dtw(ts_subject1, ts_subject2, 
    step.pattern = symmetric2)
plot(dtw_result, type="twoway")

# 2. 变化点检测
library(changepoint)

# 检测Alpha多样性的突变时间点
ts_data <- ts(shannon_diversity, frequency=1)
cpt_result <- cpt.mean(ts_data, method="PELT")
plot(cpt_result, main="Change points in Shannon diversity")
cat("变化点位置:", cpts(cpt_result), "\n")

# 3. Granger因果检验
library(vars)

# 检验微生物A的变化是否"导致"微生物B的变化
# 需要等间隔时间序列
granger_data <- data.frame(
    species_A = ts_speciesA,
    species_B = ts_speciesB)

var_model <- VAR(granger_data, p=2, type="const")
granger <- causality(var_model, cause="species_A")
cat("Granger causality p-value:", 
    granger$Granger$p.value, "\n")

# 4. 状态空间模型（检测不同"微生物组状态"及转换）
library(depmixS4)

# 隐马尔可夫模型检测微生物组状态转换
hmm_model <- depmix(
    list(PC1 ~ 1, PC2 ~ 1),  # 观测方程
    data = pca_timeseries,
    nstates = 3,              # 假设3种状态
    family = list(gaussian(), gaussian()),
    ntimes = rep(n_timepoints, n_subjects))

hmm_fit <- fit(hmm_model)
states <- posterior(hmm_fit)$state

# 可视化状态转换
plot(states, type="l", col="red", ylab="State", xlab="Time")

第六步：研究设计与样本量计算

白话解释： 纵向研究设计需要考虑：采样多少个体？每个个体采多少时间点？时间间隔多长？这些决策取决于研究问题和预期效果量。

代码示例：

# 功效分析/样本量计算
library(lmerTest)
library(simr)

# 基于LME的功效模拟
# 设定效果量和参数
pilot_model <- lmer(abundance ~ time * treatment + (1 | subject),
    data = pilot_data)

# 修改目标效果量
fixef(pilot_model)["time:treatmentB"] <- 0.5  # 设定期望效果

# 模拟功效
power_result <- powerSim(pilot_model, 
    test = fixed("time:treatmentB"),
    nsim = 200)
print(power_result)

# 增加样本量看功效变化
power_curve <- powerCurve(pilot_model,
    test = fixed("time:treatmentB"),
    along = "subject",         # 沿着个体数增加
    breaks = seq(20, 100, 10),
    nsim = 200)
plot(power_curve)

# 经验法则:
# 每组至少20-30个个体
# 每个个体至少3-5个时间点
# 关键时间点（如干预前后）必须覆盖
# 总样本量(个体×时间点) > 100

实战命令

#!/usr/bin/env Rscript
# === 微生物组纵向分析完整流程 ===

library(Maaslin2); library(vegan); library(lme4)

# 1. 数据加载
otu <- read.csv("otu_table.csv", row.names=1)
meta <- read.csv("metadata.csv")  # 含subject_id, timepoint, treatment

# 2. 纵向Alpha多样性分析
meta$shannon <- diversity(otu[meta$sample_id,], index="shannon")
lme_alpha <- lmer(shannon ~ timepoint * treatment + age + (1+timepoint|subject_id), data=meta)
summary(lme_alpha)

# 3. 纵向Beta多样性分析
dist_bc <- vegdist(otu, "bray")
# 注意：adonis2的strata参数只提供简单分层，推荐用how()设置置换方案
library(permute)
perm_design <- with(meta, how(nperm=999, blocks=subject_id))
adonis2(dist_bc ~ timepoint * treatment + age, data=meta, permutations=perm_design)

# 4. MaAsLin2差异丰度
Maaslin2(otu, meta, "maaslin2_results/",
    fixed_effects=c("timepoint","treatment","age"),
    random_effects=c("subject_id"),
    normalization="TSS", transform="AST")

面试常问点

Q1: 纵向研究为什么不能用普通t-test或Wilcoxon比较时间点？ A: 因为同一个体的重复测量不独立，违反了这些检验的独立性假设。忽略个体内相关性会导致：1）标准误被低估→假阳性增加；2）损失检测效力（未利用配对信息）。必须用混合效应模型或GEE等方法处理。

Q2: 为什么微生物组数据需要零膨胀模型？ A: 微生物组数据的零值通常远超标准分布（如负二项）的预期。零有两种来源：1）结构性零（该物种在该样本中确实不存在）；2）抽样零（存在但测序深度不够未检出）。零膨胀模型用两个子模型分别处理这两类零，更准确反映数据生成过程。

Q3: 随机截距和随机斜率分别表示什么？ A: 随机截距：每个个体有不同的基线水平（有人基线高，有人低）。随机斜率：每个个体的变化速率不同（有人变化快，有人慢）。通常先试随机截距模型，如果有理论依据或模型比较支持再加随机斜率。

Q4: 纵向微生物组研究的最少样本量是多少？ A: 取决于效果量和设计。经验法则：每组至少20-30个个体，每人3-5个时间点。如果期望效果小（如ΔPSI<0.1），需要更多。功效分析（如simr包）可以给出精确估计。总体建议：比横截面研究多50-100%的个体数。

易错点

  1. 忽略个体内相关性：使用普通回归或t-test分析纵向数据会严重膨胀I类错误率。

  2. 时间点编码错误：时间应作为连续变量（月/天）而非分类变量处理（除非时间效应明显非线性）。

  3. 不平衡数据处理不当：不同个体的采样时间点不同时，不能简单按时间点合并。混合模型能自然处理不平衡设计。

  4. 忽略dropout模式：如果丢失随访不是完全随机的（如病情严重的人更容易退出），会导致偏倚。需要检查MCAR/MAR/MNAR假设。

  5. 过度参数化：时间点少而随机效应结构复杂会导致模型不收敛。根据数据量选择合理的随机效应结构。

  6. 多重比较：同时检验数百个分类单元，必须做FDR校正。MaAsLin2默认BH校正。

补充知识

纵向分析方法选择决策树

数据类型是什么？
├── 相对丰度(0-1) → Beta回归/ZIBR
├── CLR转化后 → LME (lmer)
├── 原始计数 → NB-GLMM (glmmTMB)
└── 二分类(存在/缺失) → 逻辑回归GLMM

是否有大量零值？
├── 是(>50%零) → 零膨胀模型(ZINB/ZIBR)
└── 否 → 标准模型即可

样本量大小？
├── 小(n<30) → 非参数/permutation方法
├── 中(30-200) → LME/GLMM
└── 大(>200) → GEE也可以

分析目标？
├── 单分类单元 → LME/ZIBR
├── 多分类单元 → MaAsLin2/NBZIMM
├── 群落水平 → PERMANOVA(分层)
└── 时间模式 → DTW/HMM/变化点检测

常用纵向分析R包

R包方法适用
lme4/lmerTestLMECLR转化后连续数据
glmmTMBGLMM (NB/ZIP/ZINB)计数数据，零膨胀
ZIBR零膨胀Beta回归相对丰度，纵向
NBZIMMNB-GLMM/ZINB-GLMM微生物计数，纵向
MaAsLin2线性模型+随机效应一站式微生物组
simr功效分析混合模型样本量

📚 相关文章推荐

微生物组纵向数据分析
38% 相关
微生物组时间序列分析
35% 相关
混合效应模型纵向数据
35% 相关
676 微生物组零膨胀模型
33% 相关
零膨胀模型微生物组
29% 相关