184_微生物群落组装分析¶

library(picante)
library(vegan)
library(ape)

# 读取数据
otu_table <- read.csv("otu_table.csv", row.names = 1)
tree <- read.tree("phylogenetic_tree.nwk")

# 确保OTU表和树匹配
common_taxa <- intersect(colnames(otu_table), tree$tip.label)
otu_table <- otu_table[, common_taxa]
tree <- keep.tip(tree, common_taxa)

# 计算βMNTD（Beta Mean Nearest Taxon Distance）
# 首先计算系统发育距离矩阵
phylo_dist <- cophenetic(tree)

# 观测的βMNTD
beta_mntd_obs <- comdistnt(otu_table, phylo_dist, abundance.weighted = TRUE)

# 零模型：999次随机化
n_rand <- 999
beta_mntd_rand <- array(NA, dim = c(nrow(otu_table), nrow(otu_table), n_rand))

set.seed(42)
for (i in 1:n_rand) {
  # 随机打乱物种在树上的位置
  rand_tree <- tree
  rand_tree$tip.label <- sample(rand_tree$tip.label)
  rand_phylo_dist <- cophenetic(rand_tree)
  beta_mntd_rand[, , i] <- as.matrix(comdistnt(otu_table, rand_phylo_dist,
                                                  abundance.weighted = TRUE))
  if (i %% 100 == 0) cat("Randomization:", i, "/", n_rand, "\n")
}

# 计算βNTI
beta_nti <- matrix(NA, nrow(otu_table), nrow(otu_table))
for (j in 1:(nrow(otu_table) - 1)) {
  for (k in (j + 1):nrow(otu_table)) {
    rand_vals <- beta_mntd_rand[j, k, ]
    beta_nti[j, k] <- (beta_mntd_obs[j, k] - mean(rand_vals)) / sd(rand_vals)
    beta_nti[k, j] <- beta_nti[j, k]
  }
}

rownames(beta_nti) <- colnames(beta_nti) <- rownames(otu_table)

# βNTI解读
cat("同质选择 (βNTI < -2):", sum(beta_nti[upper.tri(beta_nti)] < -2, na.rm = TRUE), "\n")
cat("异质选择 (βNTI > 2):", sum(beta_nti[upper.tri(beta_nti)] > 2, na.rm = TRUE), "\n")
cat("随机过程 (|βNTI| <= 2):", sum(abs(beta_nti[upper.tri(beta_nti)]) <= 2, na.rm = TRUE), "\n")

# 计算RCbray
# Raup-Crick with Bray-Curtis modification
calc_rcbray <- function(otu_table, n_rand = 999) {
  n_samples <- nrow(otu_table)
  rc_matrix <- matrix(NA, n_samples, n_samples)
  obs_bray <- as.matrix(vegdist(otu_table, method = "bray"))

  # 零模型
  gamma_richness <- ncol(otu_table)
  species_pool <- colSums(otu_table > 0) / nrow(otu_table)  # 物种出现频率

  rand_bray <- array(NA, dim = c(n_samples, n_samples, n_rand))

  for (i in 1:n_rand) {
    rand_comm <- matrix(0, n_samples, ncol(otu_table))
    for (s in 1:n_samples) {
      n_species <- sum(otu_table[s, ] > 0)
      selected <- sample(1:gamma_richness, n_species, prob = species_pool)
      rand_comm[s, selected] <- 1
      # 重新分配丰度
      rand_comm[s, selected] <- rmultinom(1, sum(otu_table[s, ]), rep(1, n_species))
    }
    rand_bray[, , i] <- as.matrix(vegdist(rand_comm, method = "bray"))
  }

  # 计算RCbray
  for (j in 1:(n_samples - 1)) {
    for (k in (j + 1):n_samples) {
      p <- sum(rand_bray[j, k, ] >= obs_bray[j, k]) / n_rand
      rc_matrix[j, k] <- (p - 0.5) * 2
      rc_matrix[k, j] <- rc_matrix[j, k]
    }
  }

  rownames(rc_matrix) <- colnames(rc_matrix) <- rownames(otu_table)
  return(rc_matrix)
}

rcbray <- calc_rcbray(otu_table, n_rand = 999)

library(ggplot2)
library(reshape2)

# 提取上三角
nti_vals <- beta_nti[upper.tri(beta_nti)]
rc_vals <- rcbray[upper.tri(rcbray)]

# 分类各过程
process <- rep(NA, length(nti_vals))
process[nti_vals > 2] <- "Variable Selection"
process[nti_vals < -2] <- "Homogeneous Selection"
process[abs(nti_vals) <= 2 & rc_vals > 0.95] <- "Dispersal Limitation"
process[abs(nti_vals) <= 2 & rc_vals < -0.95] <- "Homogenizing Dispersal"
process[abs(nti_vals) <= 2 & abs(rc_vals) <= 0.95] <- "Drift"

# 统计各过程比例
process_table <- table(process)
process_pct <- round(prop.table(process_table) * 100, 1)
print(process_pct)

# 饼图
df_process <- data.frame(
  Process = names(process_pct),
  Percentage = as.numeric(process_pct)
)

ggplot(df_process, aes(x = "", y = Percentage, fill = Process)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +
  theme_void() +
  geom_text(aes(label = paste0(Percentage, "%")),
            position = position_stack(vjust = 0.5)) +
  ggtitle("Community Assembly Processes") +
  scale_fill_brewer(palette = "Set2")

# βNTI分布直方图
ggplot(data.frame(bNTI = nti_vals), aes(x = bNTI)) +
  geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
  geom_vline(xintercept = c(-2, 2), linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(x = "βNTI", y = "Count", title = "βNTI Distribution")

library(iCAMP)

# 运行iCAMP
icamp_result <- icamp.big(
  comm = otu_table,
  tree = tree,
  pd.wd = "icamp_pd/",      # 工作目录
  rand = 999,                 # 随机化次数
  nworker = 8,                # 并行核心数
  memory.G = 16,              # 内存限制(GB)
  bin.size.limit = 24,        # 每个bin的最小物种数
  sig.index = c("Confidence", "SES.RC")
)

# 提取各过程比例
process_importance <- icamp_result$CbMPDiCBraya$processes

# 可视化
barplot(t(as.matrix(process_importance)),
        beside = FALSE,
        col = c("blue", "red", "green", "orange", "grey"),
        legend = TRUE,
        main = "Community Assembly Processes (iCAMP)")

# βNTI快速计算（使用MicEco包）
library(MicEco)
bnti <- bNTI(otu_table, tree, runs = 999, cores = 8)

# iCAMP一行运行
result <- iCAMP::icamp.big(comm, tree, pd.wd = "pd/", rand = 999)

# NST（归一化随机比率）
library(NST)
nst_result <- tNST(comm = otu_table, group = metadata$group,
                    dist.method = "bray", abundance.weighted = TRUE,
                    rand = 999, nworker = 8)