微生物组季节性分析¶

# === 数据准备 ===
library(phyloseq)
library(dplyr)

# 读取纵向微生物组数据
ps <- readRDS("longitudinal_phyloseq.rds")

# 提取元数据
meta <- sample_data(ps) %>% as.data.frame()
meta$date <- as.Date(meta$collection_date)
meta$day_of_year <- as.numeric(format(meta$date, "%j"))  # 1-365
meta$month <- as.numeric(format(meta$date, "%m"))
meta$season <- ifelse(meta$month %in% c(12, 1, 2), "Winter",
               ifelse(meta$month %in% c(3, 4, 5), "Spring",
               ifelse(meta$month %in% c(6, 7, 8), "Summer", "Autumn")))

# 检查采样分布
table(meta$month)  # 各月样本数
range(meta$date)   # 时间跨度
cat(sprintf("Duration: %d days\n", as.numeric(diff(range(meta$date)))))

# === Lomb-Scargle 周期图 ===
library(lomb)

# 对单个菌属的丰度做Lomb-Scargle分析
genus_ra <- transform_sample_counts(ps, function(x) x/sum(x))
genus_table <- as.data.frame(otu_table(genus_ra))

# 选择一个菌属
target_genus <- "Prevotella"
abundance <- as.numeric(genus_table[target_genus, ])
time_points <- as.numeric(meta$date - min(meta$date))  # 天数

# Lomb-Scargle周期图
ls_result <- lsp(abundance, times = time_points,
                  from = 7, to = 400,  # 搜索7-400天的周期
                  type = "period",
                  ofac = 10,  # 过采样因子
                  plot = TRUE)

# 查看最显著周期
cat(sprintf("Peak period: %.1f days\n", ls_result$peak.at[1]))
cat(sprintf("P-value: %.4f\n", ls_result$p.value))

# 如果peak在365天附近且p<0.05，说明有年周期性

# 批量分析所有菌属
all_genera <- rownames(genus_table)
ls_results <- data.frame()

for (genus in all_genera) {
  abd <- as.numeric(genus_table[genus, ])
  if (sum(abd > 0) < 10) next  # 跳过低流行度菌

  tryCatch({
    ls <- lsp(abd, times = time_points, from = 7, to = 400,
              type = "period", ofac = 10, plot = FALSE)
    ls_results <- rbind(ls_results, data.frame(
      genus = genus,
      peak_period = ls$peak.at[1],
      power = ls$peak,
      p_value = ls$p.value
    ))
  }, error = function(e) NULL)
}

# 筛选有年周期性的菌属（周期在300-400天且显著）
annual_periodic <- ls_results %>%
  filter(peak_period > 300 & peak_period < 400 & p_value < 0.05) %>%
  arrange(p_value)

cat(sprintf("Genera with annual periodicity: %d\n", nrow(annual_periodic)))

# === 谐波回归 ===

# 年周期的谐波回归
fit_harmonic <- function(abundance, day_of_year, subject_id = NULL) {
  # 定义正弦余弦项
  sin_term <- sin(2 * pi * day_of_year / 365.25)
  cos_term <- cos(2 * pi * day_of_year / 365.25)

  # 如果有重复测量，加入个体随机效应
  if (!is.null(subject_id)) {
    library(lme4)
    df <- data.frame(y = abundance, sin_t = sin_term, cos_t = cos_term,
                     subject = subject_id)
    model <- lmer(y ~ sin_t + cos_t + (1|subject), data = df)
  } else {
    df <- data.frame(y = abundance, sin_t = sin_term, cos_t = cos_term)
    model <- lm(y ~ sin_t + cos_t, data = df)
  }

  # 提取振幅和相位
  # 注意：R没有Python式的三元表达式，必须用if/else语句
  if (!is.null(subject_id)) {
    coefs <- fixef(model)
  } else {
    coefs <- coef(model)
  }
  A <- coefs["sin_t"]
  B <- coefs["cos_t"]
  amplitude <- sqrt(A^2 + B^2)
  phase <- atan2(B, A) * 365.25 / (2 * pi)  # 转换为天数
  if (phase < 0) phase <- phase + 365.25

  # 模型显著性与R-squared
  if (!is.null(subject_id)) {
    null_model <- lmer(y ~ (1|subject), data = df)
    p_val <- anova(null_model, model)$`Pr(>Chisq)`[2]
    # lmer模型没有r.squared，使用MuMIn::r.squaredGLMM或手动计算
    r_sq <- NA  # 混合模型的R^2需要额外计算
  } else {
    p_val <- summary(model)$fstatistic
    p_val <- pf(p_val[1], p_val[2], p_val[3], lower.tail = FALSE)
    r_sq <- summary(model)$r.squared
  }

  return(list(amplitude = amplitude, phase = phase, p_value = p_val,
              model = model, r_squared = r_sq))
}

# 对所有菌属运行谐波回归
harmonic_results <- data.frame()
for (genus in rownames(genus_table)) {
  abd <- as.numeric(genus_table[genus, ])
  if (mean(abd > 0) < 0.1) next

  result <- fit_harmonic(abd, meta$day_of_year, meta$subject_id)
  harmonic_results <- rbind(harmonic_results, data.frame(
    genus = genus,
    amplitude = result$amplitude,
    peak_day = result$phase,
    p_value = result$p_value,
    r_squared = result$r_squared
  ))
}

harmonic_results$padj <- p.adjust(harmonic_results$p_value, method = "BH")
seasonal_genera <- harmonic_results %>% filter(padj < 0.05) %>% arrange(padj)

# === Alpha多样性季节性 ===
library(vegan)

# 计算alpha多样性
alpha_div <- data.frame(
  Shannon = diversity(t(otu_table(ps)), index = "shannon"),
  Richness = specnumber(t(otu_table(ps))),
  Evenness = diversity(t(otu_table(ps))) / log(specnumber(t(otu_table(ps))))
)
alpha_div$month <- meta$month
alpha_div$day_of_year <- meta$day_of_year
alpha_div$subject <- meta$subject_id

# 谐波回归
diversity_seasonal <- fit_harmonic(alpha_div$Shannon, alpha_div$day_of_year,
                                    alpha_div$subject)
cat(sprintf("Shannon diversity seasonality p = %.4f\n", diversity_seasonal$p_value))
cat(sprintf("Amplitude = %.3f, Peak day = %.0f\n",
            diversity_seasonal$amplitude, diversity_seasonal$phase))

# 可视化月度变化（箱线图）
ggplot(alpha_div, aes(x = factor(month), y = Shannon)) +
  geom_boxplot(fill = "lightblue") +
  labs(x = "Month", y = "Shannon Diversity",
       title = "Monthly Shannon Diversity") +
  theme_minimal()

# === Beta多样性季节性 ===
# 使用PERMANOVA检验季节是否解释群落组成变异
bray_dist <- phyloseq::distance(ps, method = "bray")
adonis_season <- adonis2(bray_dist ~ season + subject_id,
                          data = meta, permutations = 999)
print(adonis_season)

# PCoA按季节着色
ordinate_result <- ordinate(ps, method = "PCoA", distance = "bray")
plot_ordination(ps, ordinate_result, color = "season") +
  stat_ellipse() +
  theme_minimal()

# === 环境驱动因素分析 ===
library(Maaslin2)

# 准备环境变量
env_data <- data.frame(
  sample_id = meta$sample_id,
  temperature = meta$avg_temperature,  # 月均温
  daylight_hours = meta$daylight,       # 日照时长
  fiber_intake = meta$dietary_fiber,    # 纤维摄入
  day_of_year = meta$day_of_year,
  subject_id = meta$subject_id
)

# MaAsLin2 纵向模型
fit_maaslin <- Maaslin2(
  input_data = as.data.frame(t(genus_table)),
  input_metadata = env_data,
  output = "maaslin2_seasonal",
  fixed_effects = c("temperature", "daylight_hours", "fiber_intake"),
  random_effects = c("subject_id"),
  normalization = "CLR",
  transform = "NONE",
  analysis_method = "LM"
)

# 分析哪些菌与温度显著关联
temp_assoc <- fit_maaslin$results %>%
  filter(metadata == "temperature" & qval < 0.1) %>%
  arrange(qval)
print(temp_assoc)

# 可视化温度-菌群关系
ggplot(data.frame(temp = env_data$temperature,
                   prevotella = as.numeric(genus_table["Prevotella", ])),
       aes(x = temp, y = prevotella)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess") +
  labs(x = "Temperature (°C)", y = "Prevotella relative abundance") +
  theme_minimal()

# === 时间序列分解 ===

# 对于等间隔数据，使用STL分解
# 首先需要创建ts对象（假设月度汇总数据）
monthly_mean <- meta %>%
  mutate(year_month = format(date, "%Y-%m")) %>%
  group_by(year_month) %>%
  summarise(
    Shannon = mean(Shannon),
    Prevotella = mean(genus_table["Prevotella", cur_group_rows()])
  ) %>%
  arrange(year_month)

# 创建时间序列对象
ts_shannon <- ts(monthly_mean$Shannon, frequency = 12, start = c(2020, 1))

# STL分解
stl_result <- stl(ts_shannon, s.window = "periodic")
plot(stl_result, main = "STL Decomposition of Shannon Diversity")

# 提取各成分
seasonal_component <- stl_result$time.series[, "seasonal"]
trend_component <- stl_result$time.series[, "trend"]
remainder <- stl_result$time.series[, "remainder"]

# 季节性强度指标（Wang, Smith & Hyndman, 2006）
# 正确公式: F_S = max(0, 1 - Var(R) / Var(S + R))
# 分母是S+R的方差，不是Var(S)+Var(R)
var_sr <- var(seasonal_component + remainder)
var_remainder <- var(remainder)
seasonality_strength <- max(0, 1 - var_remainder / var_sr)
cat(sprintf("Seasonality strength: %.3f\n", seasonality_strength))
# 接近1表示季节性主导, 接近0表示无季节性, > 0.7认为有较强季节性

# === 个体水平季节性分析 ===

# 对每个个体分别做谐波回归
subjects <- unique(meta$subject_id)
individual_seasonality <- data.frame()

for (subj in subjects) {
  subj_idx <- meta$subject_id == subj
  subj_data <- genus_table[, subj_idx]
  subj_days <- meta$day_of_year[subj_idx]

  if (sum(subj_idx) < 8) next  # 至少8个时间点

  # 对Shannon多样性
  shannon_vals <- diversity(t(subj_data), index = "shannon")

  result <- tryCatch(
    fit_harmonic(shannon_vals, subj_days, subject_id = NULL),
    error = function(e) NULL
  )

  if (!is.null(result)) {
    individual_seasonality <- rbind(individual_seasonality, data.frame(
      subject = subj,
      amplitude = result$amplitude,
      peak_day = result$phase,
      p_value = result$p_value
    ))
  }
}

# 个体季节性振幅分布
hist(individual_seasonality$amplitude, breaks = 20,
     main = "Individual Seasonality Amplitude",
     xlab = "Amplitude of annual cycle")

# 哪些因素与个体季节性敏感度相关？
# 合并个体元数据（如饮食多样性、BMI、运动量）
merged <- merge(individual_seasonality, subject_metadata, by = "subject")
cor.test(merged$amplitude, merged$dietary_diversity)
cor.test(merged$amplitude, merged$outdoor_hours)

# 快速季节性分析
library(lomb)
# Lomb-Scargle
lsp(abundance, times=time_days, from=7, to=400, type="period")
# 谐波回归
lm(abundance ~ sin(2*pi*day/365) + cos(2*pi*day/365))
# PERMANOVA
adonis2(dist ~ season + subject, data=meta)