172_微生物组网络稳定性¶

library(SpiecEasi)  # 推荐的微生物网络推断方法
library(igraph)

# === 使用SpiecEasi构建网络(推荐) ===
# SpiecEasi使用稀疏逆协方差(GLASSO)或MB方法
# 比简单相关更适合组成性数据

# 输入: OTU/ASV counts矩阵 (样本×物种)
otu_matrix <- as.matrix(otu_counts)

# 方法1: GLASSO (graphical lasso)
se_result <- spiec.easi(
  otu_matrix,
  method = "glasso",
  lambda.min.ratio = 0.01,
  nlambda = 50,
  sel.criterion = "stars",  # StARS模型选择
  pulsar.params = list(rep.num = 50, ncores = 8)
)

# 方法2: MB (Meinshausen-Bühlmann)
se_mb <- spiec.easi(
  otu_matrix,
  method = "mb",
  lambda.min.ratio = 0.01,
  nlambda = 50,
  sel.criterion = "stars",
  pulsar.params = list(rep.num = 50, ncores = 8)
)

# 提取邻接矩阵和igraph对象
adj_matrix <- getRefit(se_result)
ig <- adj2igraph(adj_matrix, vertex.attr = list(name = colnames(otu_matrix)))

# === 替代方法: SparCC (组成性数据相关) ===
# sparcc结果转为网络
# library(SpiecEasi)  # 包含sparcc
# sparcc_result <- sparcc(otu_matrix)
# 使用bootstrap P值过滤显著边

# === 基本网络统计 ===
cat(sprintf("节点数: %d\n", vcount(ig)))
cat(sprintf("边数: %d\n", ecount(ig)))
cat(sprintf("平均度: %.2f\n", mean(degree(ig))))
cat(sprintf("聚类系数: %.3f\n", transitivity(ig, type = "global")))
cat(sprintf("模块化度: %.3f\n", modularity(cluster_louvain(ig))))
cat(sprintf("图密度: %.4f\n", graph.density(ig)))

# === 网络鲁棒性: 节点移除模拟 ===

# 自然连通度计算
natural_connectivity <- function(graph) {
  adj <- as_adjacency_matrix(graph, sparse = FALSE)
  eigenvalues <- eigen(adj, only.values = TRUE)$values
  nc <- log(mean(exp(eigenvalues)))
  return(nc)
}

# === 随机攻击(Random failure) ===
robustness_random <- function(graph, n_reps = 100) {
  n_nodes <- vcount(graph)
  nc_original <- natural_connectivity(graph)

  results <- matrix(0, nrow = n_nodes, ncol = n_reps)

  for (rep in 1:n_reps) {
    g_temp <- graph
    removal_order <- sample(V(g_temp)$name)

    for (i in 1:(n_nodes - 1)) {
      g_temp <- delete_vertices(g_temp, removal_order[i])
      if (vcount(g_temp) > 0) {
        results[i, rep] <- natural_connectivity(g_temp) / nc_original
      }
    }
  }

  # 返回平均鲁棒性曲线
  return(rowMeans(results))
}

# === 定向攻击(Targeted attack: 按度从高到低移除) ===
robustness_targeted <- function(graph) {
  n_nodes <- vcount(graph)
  nc_original <- natural_connectivity(graph)
  nc_curve <- numeric(n_nodes)

  g_temp <- graph
  for (i in 1:(n_nodes - 1)) {
    # 移除当前度最高的节点
    degrees <- degree(g_temp)
    hub <- names(which.max(degrees))
    g_temp <- delete_vertices(g_temp, hub)

    if (vcount(g_temp) > 0) {
      nc_curve[i] <- natural_connectivity(g_temp) / nc_original
    }
  }
  return(nc_curve)
}

# 运行模拟
rob_random <- robustness_random(ig, n_reps = 50)
rob_targeted <- robustness_targeted(ig)

# 可视化
fraction_removed <- seq(0, 1, length.out = length(rob_random))
plot(fraction_removed, rob_random, type = "l", col = "blue", lwd = 2,
     xlab = "Fraction of nodes removed",
     ylab = "Relative natural connectivity",
     main = "Network Robustness")
lines(fraction_removed, rob_targeted, col = "red", lwd = 2)
legend("topright", c("Random failure", "Targeted attack"),
       col = c("blue", "red"), lwd = 2)

# 鲁棒性指数: 曲线下面积(AUC)
robustness_index_random <- sum(rob_random) / length(rob_random)
robustness_index_targeted <- sum(rob_targeted) / length(rob_targeted)
cat(sprintf("随机攻击鲁棒性: %.3f\n", robustness_index_random))
cat(sprintf("定向攻击鲁棒性: %.3f\n", robustness_index_targeted))

# === 模块检测 ===
# Louvain算法(快速)
modules_louvain <- cluster_louvain(ig)
# Greedy优化
modules_greedy <- cluster_fast_greedy(ig)

# 模块化度
Q_louvain <- modularity(modules_louvain)
cat(sprintf("Modularity (Louvain): %.3f\n", Q_louvain))
# Q > 0.4 通常认为有显著模块结构

# 模块成员
membership <- membership(modules_louvain)
V(ig)$module <- membership

# 模块内/模块间连接比
# Zi-Pi分析(Within-module connectivity vs Among-module connectivity)
zi_pi <- function(graph, modules) {
  n <- vcount(graph)
  adj <- as_adjacency_matrix(graph, sparse = FALSE)
  mem <- membership(modules)

  results <- data.frame(node = V(graph)$name, Zi = 0, Pi = 0)

  for (i in 1:n) {
    ki <- degree(graph, V(graph)[i])
    if (ki == 0) next

    # 模块内连接
    same_module <- which(mem == mem[i])
    ki_s <- sum(adj[i, same_module])

    # Zi: within-module degree z-score
    all_in_module <- which(mem == mem[i])
    k_module <- sapply(all_in_module, function(j) sum(adj[j, same_module]))
    results$Zi[i] <- (ki_s - mean(k_module)) / (sd(k_module) + 1e-10)

    # Pi: participation coefficient
    unique_modules <- unique(mem)
    ki_t <- sapply(unique_modules, function(m) sum(adj[i, which(mem == m)]))
    results$Pi[i] <- 1 - sum((ki_t / ki)^2)
  }
  return(results)
}

zi_pi_results <- zi_pi(ig, modules_louvain)

# 节点角色分类 (Guimerà & Amaral 2005)
# Hub连接器: Zi>2.5, Pi>0.62 (Module hubs that connect modules)
# 模块Hub: Zi>2.5, Pi<0.62 (Hubs within their module)
# 连接器: Zi<2.5, Pi>0.62 (Connectors between modules)
# 外围: Zi<2.5, Pi<0.62 (Peripheral nodes)

zi_pi_results$Role <- case_when(
  zi_pi_results$Zi > 2.5 & zi_pi_results$Pi > 0.62 ~ "Network Hub",
  zi_pi_results$Zi > 2.5 & zi_pi_results$Pi <= 0.62 ~ "Module Hub",
  zi_pi_results$Zi <= 2.5 & zi_pi_results$Pi > 0.62 ~ "Connector",
  TRUE ~ "Peripheral"
)

table(zi_pi_results$Role)

# === 综合指标鉴定Keystone种 ===
keystone_score <- data.frame(
  Node = V(ig)$name,
  Degree = degree(ig),
  Betweenness = betweenness(ig, normalized = TRUE),
  Closeness = closeness(ig, normalized = TRUE),
  Eigenvector = eigen_centrality(ig)$vector,
  Zi = zi_pi_results$Zi,
  Pi = zi_pi_results$Pi,
  Role = zi_pi_results$Role
)

# 综合得分(多个中心性指标的标准化加权和)
keystone_score$composite <- scale(keystone_score$Degree) * 0.25 +
  scale(keystone_score$Betweenness) * 0.30 +
  scale(keystone_score$Closeness) * 0.15 +
  scale(keystone_score$Eigenvector) * 0.30

# Top keystone候选
keystones <- keystone_score[order(-keystone_score$composite), ][1:10, ]
print(keystones[, c("Node", "Degree", "Betweenness", "Role", "composite")])

# === 组间网络稳定性比较 ===
# 比较健康vs疾病的网络稳定性
# 分别构建两组的网络并比较拓扑指标
compare_networks <- data.frame(
  Metric = c("Nodes", "Edges", "Avg_Degree", "Clustering", "Modularity",
             "Robustness_random", "Robustness_targeted", "Neg_edge_ratio"),
  Healthy = c(vcount(ig_healthy), ecount(ig_healthy), 
              mean(degree(ig_healthy)), transitivity(ig_healthy),
              modularity(cluster_louvain(ig_healthy)),
              rob_healthy_random, rob_healthy_targeted,
              neg_ratio_healthy),
  Disease = c(vcount(ig_disease), ecount(ig_disease),
              mean(degree(ig_disease)), transitivity(ig_disease),
              modularity(cluster_louvain(ig_disease)),
              rob_disease_random, rob_disease_targeted,
              neg_ratio_disease)
)
print(compare_networks)

# === 微生物组网络分析流程 ===
Rscript network_analysis.R --otu_table otu_table.csv \
  --metadata metadata.csv \
  --method spieceasi \
  --groups "Healthy,Disease" \
  --output network_results/