微生物组时间序列分析¶

# QIIME2 纵向数据预处理
# 导入所有时间点的数据
qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path manifest_all_timepoints.tsv \
  --output-path all_reads.qza \
  --input-format PairedEndFastqManifestPhred33V2

# DADA2去噪（统一参数）
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs all_reads.qza \
  --p-trunc-len-f 250 \
  --p-trunc-len-r 200 \
  --p-n-threads 16 \
  --o-table table.qza \
  --o-representative-sequences rep-seqs.qza \
  --o-denoising-stats denoise-stats.qza

library(phyloseq)
library(tidyverse)

# 读取phyloseq对象
ps <- readRDS("phyloseq_object.rds")

# 检查纵向结构
meta <- sample_data(ps) %>% as.data.frame()
meta$subject_id <- as.factor(meta$subject_id)
meta$timepoint <- as.numeric(meta$day)

# 数据完整性汇总
completeness <- meta %>%
  group_by(subject_id) %>%
  summarise(
    n_timepoints = n(),
    first_time = min(timepoint),
    last_time = max(timepoint),
    group = first(treatment_group)
  )

cat("受试者数:", n_distinct(meta$subject_id), "\n")
cat("每人平均时间点数:", mean(completeness$n_timepoints), "\n")
cat("完整随访率:", mean(completeness$n_timepoints == max(completeness$n_timepoints)), "\n")

# 转换为长格式丰度表
otu_long <- psmelt(ps) %>%
  select(Sample, OTU, Abundance, subject_id, timepoint, treatment_group) %>%
  group_by(Sample) %>%
  mutate(rel_abundance = Abundance / sum(Abundance))

library(ggplot2)
library(gganimate)

# 1. α多样性轨迹图
alpha_div <- estimate_richness(ps, measures = c("Shannon", "Observed"))
alpha_div$subject_id <- meta$subject_id
alpha_div$timepoint <- meta$timepoint
alpha_div$group <- meta$treatment_group

ggplot(alpha_div, aes(x=timepoint, y=Shannon, group=subject_id, color=group)) +
  geom_line(alpha=0.3) +
  geom_smooth(aes(group=group), method="loess", se=TRUE, linewidth=1.5) +
  labs(x="时间（天）", y="Shannon多样性", title="微生物组α多样性时间轨迹") +
  scale_color_brewer(palette="Set1") +
  theme_minimal(base_size=14)

# 2. 关键分类群丰度变化
top_genera <- otu_long %>%
  group_by(OTU) %>%
  summarise(mean_abund = mean(rel_abundance)) %>%
  top_n(10, mean_abund) %>%
  pull(OTU)

otu_long %>%
  filter(OTU %in% top_genera) %>%
  ggplot(aes(x=timepoint, y=rel_abundance, color=treatment_group)) +
  geom_smooth(method="loess", se=TRUE) +
  facet_wrap(~OTU, scales="free_y", ncol=5) +
  labs(x="时间（天）", y="相对丰度") +
  theme_minimal()

# 3. 群落组成堆叠面积图（按时间）
otu_long %>%
  filter(OTU %in% top_genera) %>%
  group_by(timepoint, treatment_group, OTU) %>%
  summarise(mean_abund = mean(rel_abundance), .groups="drop") %>%
  ggplot(aes(x=timepoint, y=mean_abund, fill=OTU)) +
  geom_area() +
  facet_wrap(~treatment_group) +
  labs(x="时间（天）", y="平均相对丰度") +
  theme_minimal()

library(lme4)
library(lmerTest)  # 提供p值
library(emmeans)

# 拟合线性混合模型
model_shannon <- lmer(
  Shannon ~ timepoint * treatment_group + (1 + timepoint | subject_id),
  data = alpha_div
)

# 模型摘要
summary(model_shannon)

# ANOVA表（Type III检验）
anova(model_shannon, type = 3)

# 如果交互作用显著，进行简单斜率分析
emtrends(model_shannon, pairwise ~ treatment_group, var = "timepoint")

# 可视化模型拟合
library(effects)
plot(allEffects(model_shannon))

# QIIME2 longitudinal 分析
qiime longitudinal linear-mixed-effects \
  --m-metadata-file metadata.tsv \
  --m-metadata-file alpha_diversity.qza \
  --p-metric shannon_entropy \
  --p-state-column day \
  --p-individual-id-column subject_id \
  --p-group-columns treatment_group \
  --o-visualization lme_shannon.qzv

library(vegan)
library(phyloseq)

# 计算Bray-Curtis距离
bc_dist <- distance(ps, method="bray")

# 1. PERMANOVA with strata（考虑重复测量）
adonis_result <- adonis2(
  bc_dist ~ timepoint * treatment_group,
  data = meta,
  strata = meta$subject_id,
  permutations = 999
)
print(adonis_result)

# 2. 计算个体内时间变化量（turnover）
calc_turnover <- function(ps, meta) {
  subjects <- unique(meta$subject_id)
  turnover_list <- list()

  for (subj in subjects) {
    subj_samples <- meta %>% filter(subject_id == subj) %>% arrange(timepoint)
    if (nrow(subj_samples) < 2) next

    for (i in 2:nrow(subj_samples)) {
      s1 <- subj_samples$sample_id[i-1]
      s2 <- subj_samples$sample_id[i]
      bc <- as.matrix(bc_dist)[s1, s2]
      turnover_list[[length(turnover_list)+1]] <- data.frame(
        subject_id = subj,
        from_time = subj_samples$timepoint[i-1],
        to_time = subj_samples$timepoint[i],
        turnover = bc,
        group = subj_samples$treatment_group[1]
      )
    }
  }
  bind_rows(turnover_list)
}

turnover_df <- calc_turnover(ps, meta)

# 比较组间turnover差异
wilcox.test(turnover ~ group, data = turnover_df)

# 可视化
ggplot(turnover_df, aes(x=to_time, y=turnover, color=group)) +
  geom_smooth(method="loess", se=TRUE) +
  labs(x="时间（天）", y="Bray-Curtis变化量", title="群落组成变化速率") +
  theme_minimal()

library(Maaslin2)

# 准备输入数据
features <- as.data.frame(otu_table(ps))  # 特征表（行=样本, 列=特征）
metadata <- as.data.frame(sample_data(ps))

# MaAsLin2纵向分析
fit_data <- Maaslin2(
  input_data = features,
  input_metadata = metadata,
  output = "maaslin2_longitudinal",
  fixed_effects = c("timepoint", "treatment_group", "timepoint:treatment_group"),
  random_effects = c("subject_id"),
  normalization = "TSS",
  transform = "LOG",
  min_abundance = 0.001,
  min_prevalence = 0.1,
  max_significance = 0.05,
  cores = 8
)

# 查看显著的时间×组别交互效应
sig_interaction <- fit_data$results %>%
  filter(metadata == "timepoint:treatment_group" & qval < 0.05) %>%
  arrange(qval)

print(sig_interaction)

# QIIME2 波动性分析
qiime longitudinal volatility \
  --m-metadata-file metadata.tsv \
  --m-metadata-file shannon.qza \
  --p-state-column day \
  --p-individual-id-column subject_id \
  --p-default-group-column treatment_group \
  --p-default-metric shannon_entropy \
  --o-visualization volatility.qzv

# R中手动计算波动性指标
volatility_metrics <- meta %>%
  left_join(alpha_div, by=c("sample_id")) %>%
  group_by(subject_id) %>%
  arrange(timepoint) %>%
  summarise(
    group = first(treatment_group),
    mean_shannon = mean(Shannon),
    sd_shannon = sd(Shannon),
    cv_shannon = sd(Shannon) / mean(Shannon),  # 变异系数
    range_shannon = max(Shannon) - min(Shannon),
    # 相邻时间点变化量的均值
    mean_delta = mean(abs(diff(Shannon))),
    n_timepoints = n()
  )

# 比较组间波动性
wilcox.test(cv_shannon ~ group, data = volatility_metrics)

# 可视化
ggplot(volatility_metrics, aes(x=group, y=cv_shannon, fill=group)) +
  geom_boxplot() +
  geom_jitter(width=0.2, alpha=0.5) +
  labs(x="组别", y="Shannon变异系数", title="微生物组稳定性比较") +
  theme_minimal()

library(DirichletMultinomial)
library(markovchain)

# 1. DMM聚类确定群落状态
count_matrix <- as.matrix(otu_table(ps))  # 样本×OTU
# 尝试不同的k值
fit_list <- lapply(1:6, function(k) dmn(count_matrix, k))

# 选择最优k（最小Laplace近似）
laplace_vals <- sapply(fit_list, laplace)
best_k <- which.min(laplace_vals)
cat("最优状态数:", best_k, "\n")

# 分配状态标签
best_fit <- fit_list[[best_k]]
state_assignments <- apply(mixture(best_fit), 1, which.max)
meta$community_state <- paste0("CST", state_assignments)

# 2. 构建转移矩阵
transitions <- meta %>%
  arrange(subject_id, timepoint) %>%
  group_by(subject_id) %>%
  mutate(next_state = lead(community_state)) %>%
  filter(!is.na(next_state))

# 估计转移概率
trans_matrix <- table(transitions$community_state, transitions$next_state)
trans_prob <- prop.table(trans_matrix, margin=1)
print(round(trans_prob, 3))

# 3. 拟合马尔可夫链
mc <- new("markovchain",
  states = colnames(trans_prob),
  transitionMatrix = as.matrix(trans_prob)
)

# 稳态分布
steady <- steadyStates(mc)
cat("稳态分布:\n")
print(round(steady, 3))

# 差异中的差异（DID）分析
library(lme4)

# 标记干预前后
meta$post_intervention <- ifelse(meta$timepoint > intervention_day, 1, 0)

# DID模型
did_model <- lmer(
  Shannon ~ post_intervention * treatment_group + (1 | subject_id),
  data = alpha_div_meta
)
summary(did_model)

# 恢复时间分析
library(survival)

# 定义恢复事件
recovery_data <- meta %>%
  filter(treatment_group == "antibiotic", timepoint > intervention_day) %>%
  group_by(subject_id) %>%
  mutate(
    baseline_dist = bc_to_baseline,  # 与基线的Bray-Curtis距离
    recovered = baseline_dist < recovery_threshold,
    time_since_intervention = timepoint - intervention_day
  ) %>%
  summarise(
    recovery_time = min(time_since_intervention[recovered], na.rm=TRUE),
    censored = all(!recovered)
  )

# Kaplan-Meier恢复曲线
surv_obj <- Surv(recovery_data$recovery_time, !recovery_data$censored)
km_fit <- survfit(surv_obj ~ 1)
plot(km_fit, xlab="干预后天数", ylab="未恢复比例", main="微生物组恢复曲线")

library(vars)
library(SparseTSCGM)

# 选择几个关键属进行VAR分析
key_genera <- c("Bacteroides", "Faecalibacterium", "Prevotella", "Ruminococcus")

# 对单个受试者拟合VAR模型（需要足够多的时间点）
subj_data <- otu_long %>%
  filter(subject_id == "S001", OTU %in% key_genera) %>%
  pivot_wider(names_from=OTU, values_from=rel_abundance) %>%
  arrange(timepoint) %>%
  select(all_of(key_genera))

# CLR变换（处理组合数据）
subj_clr <- compositions::clr(subj_data + 1e-6)

# 选择最优滞后阶数
var_select <- VARselect(subj_clr, lag.max=5, type="const")
optimal_lag <- var_select$selection["AIC(n)"]

# 拟合VAR模型
var_model <- VAR(subj_clr, p=optimal_lag, type="const")

# Granger因果检验
causality(var_model, cause="Bacteroides")

# 脉冲响应函数（一个属的扰动如何影响其他属）
irf_result <- irf(var_model, impulse="Bacteroides", response="Faecalibacterium", n.ahead=10)
plot(irf_result)

#!/usr/bin/env Rscript
# =====================================================
# 微生物组时间序列分析完整流程
# 场景：抗生素干预研究，追踪肠道菌群变化与恢复
# =====================================================

# ---- 加载包 ----
library(phyloseq)
library(tidyverse)
library(lme4)
library(lmerTest)
library(vegan)
library(Maaslin2)
library(ggplot2)

# ---- 1. 数据加载 ----
ps <- readRDS("data/phyloseq_longitudinal.rds")
meta <- as.data.frame(sample_data(ps))

# 确保关键变量类型正确
meta$subject_id <- as.factor(meta$subject_id)
meta$day <- as.numeric(meta$day)
meta$group <- factor(meta$group, levels=c("control","antibiotic"))
meta$phase <- factor(
  ifelse(meta$day < 0, "baseline",
    ifelse(meta$day <= 7, "treatment", "recovery")),
  levels=c("baseline","treatment","recovery")
)
sample_data(ps) <- sample_data(meta)

cat("样本数:", nsamples(ps), "\n")
cat("受试者数:", n_distinct(meta$subject_id), "\n")
cat("时间点范围:", range(meta$day), "\n")

# ---- 2. α多样性时序分析 ----
alpha <- estimate_richness(ps, measures=c("Shannon","Observed","Simpson"))
alpha$sample_id <- rownames(alpha)
alpha_meta <- cbind(alpha, meta)

# 线性混合效应模型
lme_shannon <- lmer(
  Shannon ~ day * group + phase * group + (1 + day | subject_id),
  data = alpha_meta
)
cat("\n=== Shannon多样性混合模型 ===\n")
print(anova(lme_shannon, type=3))

# 可视化
p_alpha <- ggplot(alpha_meta, aes(x=day, y=Shannon, color=group)) +
  geom_line(aes(group=subject_id), alpha=0.2) +
  geom_smooth(method="loess", linewidth=1.5, se=TRUE) +
  geom_vline(xintercept=c(0, 7), linetype="dashed", color="grey50") +
  annotate("rect", xmin=0, xmax=7, ymin=-Inf, ymax=Inf, alpha=0.1, fill="red") +
  labs(x="天数（Day 0 = 干预开始）", y="Shannon多样性",
       title="肠道微生物组α多样性纵向变化") +
  scale_color_brewer(palette="Set1", name="组别") +
  theme_minimal(base_size=13)
ggsave("figures/alpha_longitudinal.pdf", p_alpha, width=10, height=6)

# ---- 3. β多样性时序分析 ----
bc_dist <- distance(ps, method="bray")

# 按时间切片的PERMANOVA
timepoints <- sort(unique(meta$day))
permanova_results <- lapply(timepoints, function(t) {
  samps <- meta$sample_id[meta$day == t]
  if (length(samps) < 10) return(NULL)
  sub_dist <- usedist::dist_subset(bc_dist, samps)
  sub_meta <- meta[meta$day == t, ]
  result <- adonis2(sub_dist ~ group, data=sub_meta, permutations=999)
  data.frame(day=t, R2=result$R2[1], pvalue=result$`Pr(>F)`[1])
}) %>% bind_rows()

# ---- 4. 差异丰度分析（MaAsLin2） ----
features_df <- as.data.frame(t(otu_table(ps)))
metadata_df <- as.data.frame(sample_data(ps))

maaslin_result <- Maaslin2(
  input_data = features_df,
  input_metadata = metadata_df,
  output = "results/maaslin2_longitudinal",
  fixed_effects = c("day", "group", "day:group", "phase"),
  random_effects = c("subject_id"),
  normalization = "TSS",
  transform = "LOG",
  min_abundance = 0.001,
  min_prevalence = 0.1,
  cores = 8
)

# ---- 5. 波动性计算 ----
volatility <- alpha_meta %>%
  group_by(subject_id, group) %>%
  arrange(day) %>%
  summarise(
    cv_shannon = sd(Shannon) / mean(Shannon),
    mean_delta = mean(abs(diff(Shannon))),
    .groups = "drop"
  )

vol_test <- wilcox.test(cv_shannon ~ group, data=volatility)
cat(sprintf("\n波动性比较 (Wilcoxon): p = %.4f\n", vol_test$p.value))

# ---- 6. 恢复时间分析 ----
# 计算每个样本与该受试者基线的距离
baseline_samples <- meta %>% filter(phase == "baseline") %>%
  group_by(subject_id) %>% slice_tail(n=1) %>% pull(sample_id)

recovery_dist <- meta %>%
  filter(group == "antibiotic", day > 7) %>%
  rowwise() %>%
  mutate(
    baseline_sample = baseline_samples[match(subject_id, 
      meta$subject_id[meta$sample_id %in% baseline_samples])],
    dist_to_baseline = as.matrix(bc_dist)[sample_id, baseline_sample]
  )

p_recovery <- ggplot(recovery_dist, aes(x=day, y=dist_to_baseline)) +
  geom_line(aes(group=subject_id), alpha=0.3) +
  geom_smooth(method="loess", color="red") +
  geom_hline(yintercept=0.3, linetype="dashed", color="blue") +
  labs(x="天数", y="与基线的Bray-Curtis距离", title="抗生素组微生物组恢复轨迹") +
  theme_minimal()
ggsave("figures/recovery_trajectory.pdf", p_recovery, width=8, height=5)

cat("\n=== 分析完成 ===\n")