微生物组纵向数据分析¶

一句话概述¶

微生物组纵向(Longitudinal)数据分析研究微生物群落随时间的变化规律——白话说就是：不是只拍一张"快照"看微生物组成，而是连续多个时间点追踪同一批人的肠道菌群变化，看看菌群是怎么随治疗/饮食/疾病发展而改变的。

核心知识点表格¶

知识点	说明
纵向研究	同一批对象在多个时间点重复采样
SplinectomeR	R包，用平滑样条比较两组的时间趋势
MicrobiomeStat	纵向微生物组数据分析的综合R包
线性混合效应模型(LMM)	考虑个体间差异的统计模型
ZIBR	零膨胀Beta回归，处理组成型+零膨胀数据
ZIGMM	零膨胀广义混合模型
q2-longitudinal	QIIME2的纵向分析插件
组成型数据	微生物丰度是相对比例，总和=1的约束
挑战	不规则采样、缺失数据、dropout、零膨胀
volatility	QIIME2中追踪个体微生物组波动的方法

各步骤详解¶

第一步：理解纵向数据的特殊性¶

白话解释： 纵向数据和横截面数据（只采一次样）的最大区别是：同一个人的多次测量是"相关"的——你今天的肠道菌群和昨天的很像。普通统计方法（如t检验）会忽略这种相关性，导致假阳性或假阴性。

纵向数据的三大挑战：
1. 个体内相关性 → 需要混合效应模型
2. 不规则采样 → 每个人的采样时间点可能不同
3. 零膨胀 → 很多菌在很多时间点检测不到（值为0）

常用研究设计：
- 前瞻性队列：健康人随访，观察谁生病
- 干预研究：用药/饮食干预前后对比
- 发育研究：婴儿出生后菌群建立过程

第二步：SplinectomeR分析¶

白话解释： SplinectomeR用平滑样条（spline）来拟合每组的时间趋势线，然后通过置换检验比较两条线是否有显著差异。

# 安装SplinectomeR
# devtools::install_github("RRShieldsCutler/splinectomeR")
library(splinectomeR)                    # 加载包

# 准备数据格式
# 数据框需包含：sample_id, subject_id, group, time, value
# value可以是alpha多样性指标或某个菌的丰度

# 示例数据
long_data <- data.frame(
  subject = rep(paste0("S", 1:20), each = 5),  # 20个受试者
  group = rep(c("Treatment", "Control"), each = 50),  # 两组
  time = rep(c(0, 7, 14, 28, 56), 20),   # 5个时间点（天）
  shannon = rnorm(100, mean = 3, sd = 0.5)  # Shannon多样性（示例）
)

# ===== permuspliner：比较两组的总体趋势 =====
# 核心问题：两组的时间趋势是否显著不同？
result_perm <- permuspliner(
  data = long_data,                      # 输入数据
  x = "time",                            # x轴：时间
  y = "shannon",                         # y轴：响应变量
  cases = "subject",                     # 受试者ID
  category = "group",                    # 分组变量
  groups = c("Treatment", "Control"),    # 比较的两组
  perms = 999,                           # 置换次数
  cut_low = 0,                           # 时间范围下限
  cut_high = 56                          # 时间范围上限
)
# p值 < 0.05 → 两组趋势显著不同

# ===== sliding_spliner：在哪些时间段有差异？ =====
# 核心问题：两组在哪些具体时间点开始出现差异？
result_slide <- sliding_spliner(
  data = long_data,                      # 输入数据
  x = "time",                            # 时间变量
  y = "shannon",                         # 响应变量
  cases = "subject",                     # 受试者ID
  category = "group",                    # 分组变量
  groups = c("Treatment", "Control"),    # 两组
  perms = 999,                           # 置换次数
  ints = 100                             # 检测的时间间隔数
)
# 输出：每个时间点的p值 → 找到差异出现的具体时间段

# 可视化
# SplinectomeR内置了可视化函数
# 展示两组的平滑趋势线及置信区间

第三步：MicrobiomeStat综合分析¶

白话解释： MicrobiomeStat是一个"全家桶"式的R包，专门针对纵向微生物组数据，提供从数据导入到可视化的完整流程。

# 安装MicrobiomeStat
# install.packages("MicrobiomeStat")
library(MicrobiomeStat)                  # 加载包

# 从phyloseq对象转换
# mstat_obj <- mStat_convert_phyloseq_to_data_obj(ps_obj)

# 或直接构建数据对象
mstat_obj <- list(
  feature.tab = otu_table,               # OTU/ASV丰度矩阵
  meta.dat = sample_metadata,            # 样本元数据（含时间、分组、受试者ID）
  feature.ann = taxonomy_table           # 分类注释
)

# ===== 纵向Alpha多样性变化 =====
generate_alpha_change_test_long(
  data.obj = mstat_obj,                  # 数据对象
  alpha.name = c("shannon", "observed_species"),  # 多样性指标
  time.var = "Timepoint",                # 时间变量名
  subject.var = "SubjectID",             # 受试者ID
  group.var = "Group",                   # 分组变量
  adj.vars = c("Age", "Sex"),            # 协变量校正
  change.base = "Baseline"              # 基线时间点
)

# ===== 纵向Beta多样性变化 =====
generate_beta_change_test_long(
  data.obj = mstat_obj,                  # 数据对象
  dist.name = "bray",                    # 距离度量
  time.var = "Timepoint",                # 时间变量
  subject.var = "SubjectID",             # 受试者ID
  group.var = "Group",                   # 分组变量
  change.base = "Baseline"              # 基线
)

# ===== 纵向差异丰度分析 =====
generate_taxa_change_test_long(
  data.obj = mstat_obj,                  # 数据对象
  time.var = "Timepoint",                # 时间变量
  subject.var = "SubjectID",             # 受试者ID
  group.var = "Group",                   # 分组变量
  feature.level = "Genus",               # 分析级别：属
  change.base = "Baseline",             # 基线
  feature.dat.type = "proportion"        # 数据类型：比例
)

第四步：线性混合效应模型¶

白话解释： LMM是纵向数据分析的"标准武器"。它把每个人当作一个"随机效应"，允许每个人有不同的基线水平和变化速度。

library(lme4)                            # 混合效应模型
library(lmerTest)                        # 提供p值

# 对单个菌属的纵向变化建模
# 因变量：菌属丰度（CLR转换后）
# 固定效应：时间、分组、时间×分组交互
# 随机效应：受试者（截距和斜率）

model <- lmer(
  clr_abundance ~                        # CLR转换后的丰度
    Time * Group +                       # 时间×分组交互效应（核心！）
    Age + Sex +                          # 协变量
    (1 + Time | SubjectID),              # 随机截距+随机斜率
  data = long_data                       # 数据
)

summary(model)                           # 查看结果

# 关键看Time:Group交互效应的p值
# p < 0.05 → 两组的时间变化趋势显著不同

# 对所有菌属批量运行（多重检验校正）
all_genera <- unique(long_data$Genus)    # 所有菌属
results <- data.frame()                  # 存结果

for (genus in all_genera) {              # 遍历每个属
  subset_data <- long_data[long_data$Genus == genus, ]
  tryCatch({
    m <- lmer(clr_abundance ~ Time * Group + (1 | SubjectID),
              data = subset_data)
    coefs <- summary(m)$coefficients     # 提取系数
    interaction_p <- coefs["Time:GroupTreatment", "Pr(>|t|)"]  # 交互p值
    results <- rbind(results, data.frame(
      Genus = genus,
      p_value = interaction_p
    ))
  }, error = function(e) NULL)           # 模型不收敛则跳过
}

results$fdr <- p.adjust(results$p_value, method = "BH")  # FDR校正
significant <- results[results$fdr < 0.05, ]  # 显著结果

第五步：QIIME2纵向分析¶

白话解释： QIIME2的q2-longitudinal插件提供了简洁的命令行接口做纵向分析。

# QIIME2纵向分析（q2-longitudinal）

# 挥发性分析（volatility）：可视化个体水平的时间变化
qiime longitudinal volatility \          # 波动性分析
  --m-metadata-file metadata.tsv \       # 元数据（含时间和分组）
  --p-state-column Timepoint \           # 时间列
  --p-individual-id-column SubjectID \   # 个体ID列
  --p-default-group-column Group \       # 分组列
  --p-default-metric shannon \           # 度量指标
  --o-visualization volatility.qzv       # 输出可视化

# 配对差异检验：基线 vs 终点
qiime longitudinal pairwise-differences \  # 配对检验
  --m-metadata-file metadata.tsv \       # 元数据
  --p-metric shannon \                   # 检测指标
  --p-state-column Timepoint \           # 时间列
  --p-state-1 Baseline \                 # 时间点1
  --p-state-2 Week12 \                   # 时间点2
  --p-individual-id-column SubjectID \   # 个体ID
  --p-group-column Group \               # 分组
  --o-visualization pairwise.qzv         # 输出

# 线性混合效应模型
qiime longitudinal linear-mixed-effects \  # LMM
  --m-metadata-file metadata.tsv \       # 元数据
  --p-metric shannon \                   # 度量
  --p-state-column Timepoint \           # 时间
  --p-individual-id-column SubjectID \   # 个体ID
  --p-group-columns Group \              # 分组
  --o-visualization lme.qzv              # 输出

常见报错与解决¶

报错	原因	解决方案
`lmer: model failed to converge`	数据太稀疏或模型太复杂	简化随机效应结构（去掉随机斜率）
`permuspliner: too few observations`	每组受试者太少	至少每组10个受试者
`singular fit`	随机效应方差估计为0	数据变异太小，简化模型
`QIIME2: metadata type error`	时间列被当作分类变量	在metadata中确保时间列是数值型
`compositional data warning`	直接用比例数据建模	使用CLR转换或专用组成型方法（ZIBR）

速查表¶

# ========== 纵向微生物组分析速查 ==========

# SplinectomeR（两组趋势比较）
permuspliner(data, x="time", y="value", cases="subject",
  category="group", groups=c("A","B"), perms=999)

# MicrobiomeStat（综合纵向分析）
generate_alpha_change_test_long(data.obj, alpha.name="shannon",
  time.var="Time", subject.var="ID", group.var="Group")

# 线性混合效应模型
lmer(value ~ Time * Group + (1 + Time | SubjectID), data)

# QIIME2纵向分析
qiime longitudinal volatility --p-state-column Time --p-individual-id-column ID
qiime longitudinal linear-mixed-effects --p-metric shannon

# 方法选择指南
简单两组比较    → SplinectomeR
多因素复杂设计  → LMM / NBMM
零膨胀数据      → ZIBR / ZIGMM
探索性分析      → q2-longitudinal volatility
综合分析        → MicrobiomeStat

面试高频问题¶

Q1: 纵向微生物组分析和横截面分析有什么区别？

横截面分析只在一个时间点采样，只能看到"快照"；纵向分析在多个时间点追踪同一批人，能看到动态变化。纵向数据的优势：(1) 可以检测时间趋势；(2) 每个人做自己的对照（减少个体差异干扰）；(3) 可以发现因果关系线索。挑战：需要处理个体内相关性、缺失数据、dropout。

Q2: 为什么不能直接用t检验分析纵向数据？

因为同一个人的多次测量是相关的（不独立），而t检验假设样本独立。忽略这种相关性会低估标准误，导致假阳性率膨胀。必须用混合效应模型（random effect处理个体差异）或重复测量方差分析。

Q3: SplinectomeR的permuspliner和sliding_spliner有什么区别？

permuspliner回答"两组的整体趋势是否不同？"——给出一个总的p值。sliding_spliner回答"在哪些时间段两组出现差异？"——给出每个时间点的p值。前者看全局，后者看局部。两者都用置换检验，对异常值和不规则采样有一定鲁棒性。

Q4: 微生物组的组成型数据问题在纵向分析中怎么处理？

微生物丰度是相对比例（总和=1），一个菌增多必然导致其他菌"看起来"减少。处理方法：(1) CLR转换（中心对数比）——最推荐；(2) ZIBR专门处理零膨胀的Beta分布数据；(3) ANCOM-BC/LinDA在组成型框架下做差异分析。不建议直接用原始比例或raw counts建模。

Q5: 纵向研究中受试者dropout怎么处理？

dropout（中途退出）导致数据不完整。处理策略：(1) 混合效应模型天然可以处理不等长数据（MAR假设下无偏）；(2) 多重插补（MI）估计缺失值；(3) SplinectomeR用样条拟合对不规则采样有天然鲁棒性；(4) 敏感性分析：比较完整数据和不完整数据的结果是否一致。