R/modules.R

In MetaNet: Network Analysis for Omics Data

# ========6.modules=========

#' Generate a n-modules network
#'
#' @param module_number number of modules
#' @param n_node_in_module number of nodes in each modules
#' @param intra_module_density intra_module_density, recommend bigger than 20*inter_module_density, default:0.3
#' @param inter_module_density inter_module_density, default:0.01
#'
#' @description this is just a random generation method, the module number of result is not exactly the module_number, you can change the inter_module_density and intra_module_density to get the proper result.
#' @return n-modules metanet
#' @export
#'
#' @family module
#' @examples
#' g1 <- module_net()
#' get_n(g1)
#' plot(g1, mark_module = TRUE)
#' plot(g1, coors = g_layout(g1, zoom2 = 20))
#' plot(g1, coors = g_layout_polyarc(g1, group = "module"))
#' plot(g1, coors = g_layout_polygon(g1, group = "module"))
module_net <- function(module_number = 3, n_node_in_module = 30,
                       intra_module_density = 0.3,
                       inter_module_density = 0.01) {
  n_node_in_module <- rep(n_node_in_module, length = module_number)
  mat <- matrix(0, nrow = sum(n_node_in_module), ncol = sum(n_node_in_module))
  start <- c(1, cumsum(n_node_in_module[-length(n_node_in_module)]) + 1)
  end <- cumsum(n_node_in_module)
  seqls <- lapply(1:module_number, \(i)(start[i]:end[i]))
  # generate intra_module_edges
  for (i in seq_len(module_number)) {
    seq <- seqls[[i]]
    mat[seq, seq] <- as.matrix(igraph::get.adjacency(igraph::erdos.renyi.game(n_node_in_module[i], intra_module_density)))
  }
  # generate inter_module_edges
  idls <- combn(1:module_number, 2) %>% split(., col(.))
  for (ids in idls) {
    seq1 <- seqls[[ids[1]]]
    seq2 <- seqls[[ids[2]]]
    mat[seq1, seq2] <- sample(c(1, 0), length(seq1) * length(seq2),
      replace = TRUE,
      prob = c(inter_module_density, 1 - inter_module_density)
    )
  }

  g <- igraph::graph.adjacency(mat, mode = "undirected")
  # plot(g)
  c_net_update(g) -> g1
  g1 <- module_detect(g1)
  g1 <- to_module_net(g1)
  g1
}

#' Detect the modules
#'
#' @param go an igraph object
#' @param method cluster_method: "cluster_walktrap", "cluster_edge_betweenness", "cluster_fast_greedy", "cluster_spinglass"
#' @param n_node_in_module transfer the modules less than n_node_in_module to "others"
#' @param delete logical, delete others modules? default:FALSE, the others module will be "others".
#'
#' @return an igraph object
#' @export
#' @aliases c_net_module
#' @family module
#' @examples
#' data("c_net")
#' module_detect(co_net) -> co_net_modu
module_detect <- function(go, method = "cluster_fast_greedy", n_node_in_module = 0, delete = FALSE) {
  stopifnot(is_igraph(go))
  if ("original_module" %in% vertex_attr_names(go)) message("'module' already exsited, start a new module detection!")
  ms <- c("cluster_walktrap", "cluster_edge_betweenness", "cluster_fast_greedy", "cluster_spinglass")
  method <- match.arg(method, ms)
  if ("weight" %in% graph_attr_names(go)) E(go)$weight <- abs(igraph::E(go)$weight)

  switch(method,
    "cluster_walktrap" = {
      wc <- igraph::cluster_walktrap(go, weights = NULL)
    },
    "cluster_edge_betweenness" = {
      wc <- igraph::cluster_edge_betweenness(go, weights = NULL)
    },
    "cluster_fast_greedy" = {
      wc <- igraph::cluster_fast_greedy(go, weights = NULL)
    },
    "cluster_spinglass" = {
      wc <- igraph::cluster_spinglass(go, weights = NULL)
    }
  )

  # add components
  V(go)$components <- igraph::components(go)$membership %>% as.character()
  V(go)$module <- igraph::membership(wc) %>% as.character()
  V(go)$original_module <- V(go)$module

  go <- filter_n_module(go, n_node_in_module = n_node_in_module, delete = delete)
  go <- c_net_update(go)

  graph_attr(go)$communities <- wc
  graph_attr(go)$modularity <- modularity(wc)
  rand.g <- erdos.renyi.game(length(V(go)), length(E(go)), type = "gnm")
  rand_m <- modularity(cluster_fast_greedy(rand.g))
  relative_modularity <- (modularity(wc) - rand_m) / rand_m # relative modularity
  graph_attr(go)$relative_modularity <- relative_modularity

  return(go)
}

as_module <- \(x){
  x <- list(membership = x)
  vids <- names(x$membership)
  modus <- tapply(vids, x$membership, simplify = FALSE, function(x) x)
  return(modus)
}

#' Filter some modules as others
#'
#' @param go_m metanet with modules
#' @param n_node_in_module transfer the modules less than n_node_in_module to "others"
#' @param keep_id keep modules ids, will not be "others"
#' @param delete logical, delete others modules? default:FALSE, the others module will be "others".
#' @family module
#' @return metanet with modules
#' @export
#' @examples
#' data("c_net")
#' module_detect(co_net) -> co_net_modu
#' filter_n_module(co_net_modu, n_node_in_module = 30) -> co_net_modu
#' if (requireNamespace("ggtree") && requireNamespace("treeio")) plot_module_tree(co_net_modu)
#' combine_n_module(co_net_modu, 20) -> co_net_modu1
#' if (requireNamespace("ggtree") && requireNamespace("treeio")) plot_module_tree(co_net_modu1)
filter_n_module <- function(go_m, n_node_in_module = 0, keep_id = NULL, delete = FALSE) {
  module <- NULL
  if (!"original_module" %in% vertex_attr_names(go_m)) stop("'module' do not exsited, please do a `module_detect` first!")
  members <- V(go_m)$original_module
  table(members) %>% sort(decreasing = TRUE) -> s_members

  # filter modules whose nodes bigger than n_node_in_module
  keep_id1 <- names(s_members[s_members >= n_node_in_module])
  keep_id <- base::union(as.character(keep_id), as.character(keep_id1))
  new_module <- ifelse(members %in% keep_id, members, "others")

  V(go_m)$module <- new_module
  if (delete) go_m <- c_net_filter(go_m, module != "others")
  return(go_m)
}

#' Combine or cut modules to module_number
#'
#' @param module_number number of modules
#' @export
#' @rdname filter_n_module
combine_n_module <- function(go_m, module_number = 5) {
  get_community(go_m) -> comm
  igraph::cut_at(comm, module_number) -> new_modu
  V(go_m)$module <- as.character(new_modu)
  graph.attributes(go_m)$communities$membership <- as.numeric(new_modu)
  go_m
}

#' Transformation a network to a module network
#'
#' @param go metanet
#' @export
#' @family module
#' @return metanet with modules
to_module_net <- function(go) {
  if (!"module" %in% vertex_attr_names(go)) stop("no 'module', please `module_detect()` first or set the V(net)$module.")
  go <- anno_edge(go, get_v(go)[, c("name", "module")], verbose = FALSE)
  tmp_e <- igraph::edge.attributes(go)
  E(go)$e_type <- ifelse(tmp_e$module_from == tmp_e$module_to, "intra-module", "inter-module")
  # 刷新颜色
  # go=delete_edge_attr(go,"color")
  go <- c_net_set(go, vertex_class = "module")
  V(go)$color <- ifelse(V(go)$module == "others", "grey", V(go)$color)
  n_mod <- unique(V(go)$module)
  igraph::graph.attributes(go)$n_type <- "module"
  igraph::graph.attributes(go)$n_modules <- length(n_mod[n_mod != "others"])
  go
}

#' Get community
#' @param go_m module metanet
#'
#' @export
#' @family module
#' @return community
get_community <- function(go_m) {
  if (is.null(igraph::graph_attr(go_m)$communities)) stop("No community find, please do module_net() first.")
  igraph::graph_attr(go_m)$communities
}

#' Get module
#' @param go_m module metanet
#'
#' @export
#' @family module
#' @return module
get_module <- function(go_m) {
  if (!"module" %in% vertex_attr_names(go_m)) stop("no modules, please `module_detect()` first")
  setNames(V(go_m)$module, V(go_m)$name)
}

#' Get module_eigen
#' @param go_m module metanet
#'
#' @export
#' @family module
#' @return module_eigen
get_module_eigen <- function(go_m) {
  graph_attr(go_m, "module_eigen")
}

#' Summary module index
#' @param go_m module metanet
#' @param var variable name
#' @param module which column name is module. default: "module"
#' @param ... add
#' @family module
#' @export
#' @return ggplot
#' @examples
#' data("c_net")
#' module_detect(co_net, n_node_in_module = 30) -> co_net_modu
#' summary_module(co_net_modu, var = "v_class", module = "module")
#' summary_module(co_net_modu, var = "Abundance", module = "module")
summary_module <- function(go_m, var = "v_class", module = "module", ...) {
  tmp_v <- get_v(go_m)
  if ((length(module) > 1) | (length(var) > 1)) stop("var or module should be one column!")
  a <- tmp_v %>% dplyr::select(!!module, !!var)
  colnames(a)[1] <- "module"

  i <- var
  if (is.numeric(a[, i])) {
    pcutils::group_box(a[i], group = "module", metadata = a, ...)
  } else {
    table(a[, i], a$module) %>%
      as.data.frame() %>%
      reshape2::acast(Var1 ~ Var2, value.var = "Freq") %>%
      as.data.frame() -> tab
    pcutils::stackplot(tab, legend_title = var, ...) + labs(x = "Module")
  }
}

#' Plot module tree
#' @param go_m module metanet
#' @param module which column name is module. default: "module"
#' @param community community object, default: NULL, use the community of go_m
#' @param label.size label.size
#' @return ggplot
#'
#' @export
#' @rdname filter_n_module
plot_module_tree <- function(go_m, module = "module", community = NULL, label.size = 2) {
  tmp_v <- get_v(go_m)
  mdata <- tmp_v[, c("name", module)]

  lib_ps("ggtree", "treeio", library = FALSE)
  # modules tree
  if (is.null(community)) {
    get_community(go_m) %>% treeio::as.phylo() -> mcl
  } else {
    community %>% treeio::as.phylo() -> mcl
  }
  mcl <- dplyr::left_join(mcl, mdata, by = c("label" = "name"))
  p <- ggtree::ggtree(mcl, size = 0.3) +
    ggtree::geom_tiplab(aes(color = module), show.legend = FALSE, size = label.size) +
    scale_color_manual(values = pcutils::get_cols(length(unique(mdata$module)), "col3"))

  ggtree::gheatmap(p, mdata %>% tibble::column_to_rownames("name")) +
    scale_fill_manual(values = pcutils::get_cols(length(unique(mdata$module)), "col3"), name = NULL)
}

#' Calculate the eigenvalue of each module and correlation of nodes and eigenvalue (node_eigen_cor).
#'
#' @param go_m module metanet
#' @param totu original abundance table
#' @param cor_method "pearson", "kendall", "spearman"
#'
#' @export
#' @return module metanet with module_eigen
#' @rdname module_expression
module_eigen <- function(go_m, totu, cor_method = "spearman") {
  modules <- get_module(go_m)
  totu <- totu[, names(modules)]

  res <- lapply(levels(factor(modules)), \(i){
    if (i == "others") {
      return(NULL)
    }
    totu1 <- totu[, modules == i]
    # PCA
    # pc <- prcomp(totu1)
    # pc$x[, 1]
    # cor(pc$x[, 1],rowMeans(totu1))

    # SVD
    totu_scale <- scale(totu1)
    svd <- svd(totu_scale)
    if (cor(svd$u[, 1], rowMeans(totu_scale)) < 0) {
      return(-svd$u[, 1])
    } else {
      return(svd$u[, 1])
    }
  })
  names(res) <- levels(factor(modules))
  # simplify method
  eigen_res <- do.call(cbind, res) %>% data.frame(., check.names = FALSE)
  # colnames(eigen_res)=paste0("ME_",colnames(eigen_res))
  rownames(eigen_res) <- rownames(totu)

  lapply(colnames(eigen_res), \(i)cor(totu[, modules == i], eigen_res[, i], method = cor_method)) %>%
    do.call(rbind, .) %>%
    as.data.frame() -> node_eigen_cor
  colnames(node_eigen_cor) <- "node_eigen_cor"
  go_m <- anno_vertex(go_m, node_eigen_cor, verbose = FALSE)
  graph.attributes(go_m)$module_eigen <- eigen_res
  go_m
}

#' Plot the expression of each modules
#'
#' @param go_m module metanet
#' @param totu original abundance table used for module_eigen().
#' @param group group variable for totu
#' @param r_threshold the threshold for node_eigen_cor, default: 0.6.
#' @param x_order order the x axis.
#' @param facet_param parameters parse to \code{\link[ggplot2]{facet_wrap}}, e.g. nrow=2.
#' @param plot_eigen plot the eigen value line.
#' @family module
#' @export
#'
#' @examples
#' data("otutab", package = "pcutils")
#' t(otutab) -> totu
#' data("c_net")
#' module_detect(co_net, n_node_in_module = 30) -> co_net_modu
#' module_eigen(co_net_modu, totu) -> co_net_modu
#' module_expression(co_net_modu, totu)
module_expression <- function(go_m, totu, group = NULL, r_threshold = 0.6,
                              x_order = NULL, facet_param = NULL, plot_eigen = FALSE) {
  node_eigen_cor <- variable <- value <- name <- module <- rowname <- NULL
  if (is.null(graph_attr(go_m, "module_eigen"))) stop("Please do module_eigen() first")

  graph_attr(go_m, "module_eigen") %>%
    rownames_to_column() %>%
    reshape2::melt(id.vars = "rowname", variable.name = "module") -> module_eigen

  get_v(go_m) -> tmp_v

  if (!is.null(group)) totu <- pcutils::hebing(totu, group, 1)

  totu_scale <- scale(totu)
  pdat <- cbind(tmp_v[, c("name", "module", "node_eigen_cor")], t(totu_scale[, tmp_v$name]))
  pdat <- dplyr::filter(pdat, node_eigen_cor > r_threshold)
  pdat_m <- reshape2::melt(pdat, id.vars = c("name", "module", "node_eigen_cor"))
  if (!is.null(x_order)) pdat_m$variable <- pcutils::change_fac_lev(pdat_m$variable, level = x_order)
  pdat_m$module <- factor(pdat_m$module)

  p1 <- ggplot() +
    geom_line(data = pdat_m, aes(
      x = variable, y = value, group = name, color = module,
      size = node_eigen_cor, alpha = node_eigen_cor
    ), size = 0.8) +
    MetaNet_theme +
    theme(plot.margin = unit(c(1, 2, 1, 1), "lines")) +
    do.call(facet_wrap, append(list(facets = ~module), pcutils::update_param(list(scales = "free_y", ncol = 2), facet_param))) +
    labs(x = NULL, y = NULL) +
    # facet_wrap(facets = ~module,nrow = nrow,scales = "free_y")+
    scale_alpha_continuous(range = c(0, 0.6)) +
    scale_x_discrete(expand = c(0, 0)) +
    scale_color_manual(values = get_cols(nlevels(pdat_m$module)))

  if (plot_eigen) p1 <- p1 + geom_line(data = module_eigen, aes(x = rowname, y = value, col = module, group = module), size = 2, alpha = 1)

  return(p1)
}

#' Zi-Pi calculate
#'
#' @param go_m igraph object after `module_detect()`
#' @param mode use 7-group (mode=1) or 4-group (mode=2), default: mode=2
#' @param use_origin use original_module, default:TRUE, if FALSE, use module
#'
#' @return igraph
#' @export
#' @family module
#' @references 1. Guimerà, R. & Amaral, L. Functional cartography of complex metabolic networks. (2005) doi:10.1038/nature03288.
#' @examples
#' data("c_net")
#' module_detect(co_net) -> co_net_modu
#' zp_analyse(co_net_modu) -> co_net_modu
#' if (requireNamespace("ggrepel")) {
#'   zp_plot(co_net_modu)
#'   zp_plot(co_net_modu, mode = 3)
#' }
zp_analyse <- function(go_m, mode = 2, use_origin = TRUE) {
  go_m -> go1
  v_index <- get_v(go_m)
  if (!"module" %in% names(v_index)) stop("no modules, please `module_detect()` first")
  if ("roles" %in% names(v_index)) message("areadly has roles, overwrite!")

  # use original_module to do zp_analyse
  if (use_origin) {
    {
      go1 <- anno_vertex(go1, data.frame(row.names = V(go1)$name, module = V(go1)$original_module %>% as.numeric()), verbose = FALSE)
    } %>% suppressMessages()
  } else {
    if ("others" %in% v_index$module) message("Consider others as one module!")
    {
      go1 <- anno_vertex(go1, data.frame(
        row.names = V(go1)$name,
        module = tidai(v_index$module, seq_along(unique(v_index$module)))
      ), verbose = FALSE)
    } %>% suppressMessages()
  }

  within <- within_module_deg_z_score(go1)
  v_index$Ki <- within$Ki
  v_index$Zi <- within$Zi
  pc <- part_coeff(go1)
  v_index$Pi <- pc$Pi

  if (mode == 1) {
    lab <- c("Ultra-peripherals", "Peripherals", "Non-hub connectors", "Non-hub kinless nodes", "Provincial hubs", "Connector hubs", "Kinless hubs")
    backs <- data.frame(
      x1 = c(0, 0.05, 0.62, 0.8, 0, 0.3, 0.75),
      x2 = c(0.05, 0.62, 0.8, 1, 0.3, 0.75, 1),
      y1 = c(-Inf, -Inf, -Inf, -Inf, 2.5, 2.5, 2.5),
      y2 = c(2.5, 2.5, 2.5, 2.5, Inf, Inf, Inf),
      lab = factor(lab, levels = lab)
    )
  } else if (mode == 2) {
    lab <- c("Peripherals", "Network hubs", "Module hubs", "Connectors")
    backs <- data.frame(
      x1 = c(0, 0.62, 0, 0.62),
      x2 = c(0.62, 1, 0.62, 1),
      y1 = c(-Inf, 2.5, 2.5, -Inf),
      y2 = c(2.5, Inf, Inf, 2.5),
      lab = factor(lab, levels = lab)
    )
  }
  deter_role <- \(x, y, backs = backs){
    for (i in 1:nrow(backs)) {
      if (dplyr::between(as.numeric(x), backs$x1[i], backs$x2[i]) & dplyr::between(as.numeric(y), backs$y1[i], backs$y2[i])) {
        # if((backs$x1[i]<=x)&(backs$x2[i]>=x)&(backs$y1[i]<=y)&(backs$y2[i]>=y)){
        role <- backs$lab[i]
        break
      }
    }
    return(role)
  }
  v_index$roles <- apply(v_index, 1, \(x)deter_role(x["Pi"], x["Zi"], backs))

  vertex.attributes(go_m) <- as.list(v_index)
  return(go_m)
}

# calculate Zi
within_module_deg_z_score <- function(g, A = NULL, weighted = FALSE) {
  stopifnot(is_igraph(g))
  if (is.null(A)) {
    if (isTRUE(weighted)) {
      A <- as_adj(g, sparse = FALSE, names = TRUE, attr = "weight")
    } else {
      A <- as_adj(g, sparse = FALSE, names = TRUE)
    }
  }
  memb <- vertex_attr(g, "module")
  N <- max(memb)
  nS <- tabulate(memb)
  z <- Ki <- rep.int(0, dim(A)[1L])
  Ksi <- sigKsi <- rep.int(0, N)
  names(z) <- names(Ki) <- rownames(A)
  for (S in seq_len(N)) {
    x <- rowSums(as.matrix(A[memb == S, memb == S]))
    Ki[memb == S] <- x
    Ksi[S] <- sum(x) / nS[S]
    sigKsi[S] <- sqrt(sum((x - Ksi[S])^2) / (nS[S] - 1))
  }
  z <- (Ki - Ksi[memb]) / sigKsi[memb]
  z[is.infinite(z)] <- 0
  z[is.nan(z)] <- 0
  Zi <- z
  df <- data.frame(Ki, Zi, row.names = names(Ki))
  return(df)
}

# calculate Pi
part_coeff <- function(g, A = NULL, weighted = FALSE) {
  stopifnot(is_igraph(g))
  if (is.null(A)) {
    if (isTRUE(weighted)) {
      A <- as_adj(g, sparse = FALSE, attr = "weight")
    } else {
      A <- as_adj(g, sparse = FALSE)
    }
  }
  memb <- vertex_attr(g, "module")
  Ki <- colSums(A)
  N <- max(memb)
  Kis <- t(rowsum(A, memb))
  Pi <- 1 - ((1 / Ki^2) * rowSums(Kis^2))
  names(Pi) <- rownames(A)
  Pi <- data.frame(Pi)
  return(Pi)
}

#' Zi-Pi plot of vertexes
#'
#' @param go igraph object after zp_analyse()
#' @param label show label or not
#' @param mode plot style, 1~3
#'
#' @return a ggplot object
#' @export
#' @rdname zp_analyse
zp_plot <- function(go, label = TRUE, mode = 1) {
  v_class <- value <- size <- roles <- x1 <- x2 <- y1 <- y2 <- Pi <- Zi <- name <- NULL
  lib_ps("ggrepel", library = FALSE)
  get_v(go) -> taxa.roles
  if (!"roles" %in% names(taxa.roles)) stop("no roles, please zp_analyse() first")

  if (mode == 3) {
    reshape2::melt(taxa.roles, measure.vars = c("Zi", "Pi")) -> taxa.roles1
    p <- ggplot(taxa.roles1, aes(x = v_class, y = value, col = v_class, size = size, shape = roles)) +
      geom_point() +
      facet_grid(variable ~ ., scales = "free_y") +
      theme_bw() +
      labs(x = NULL, y = NULL) +
      guides(col = "none") +
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
      scale_color_manual(values = setNames(unique(taxa.roles$color), unique(taxa.roles$v_class)))
    return(p)
  }

  mode <- ifelse(nlevels(V(go)$roles) == 7, 1, 2)
  if (mode == 1) {
    lab <- c("Ultra-peripherals", "Peripherals", "Non-hub connectors", "Non-hub kinless nodes", "Provincial hubs", "Connector hubs", "Kinless hubs")
    CPCOLS <- c("#FCF6EFFC", "#EEBCF5", "#EDEDA4", "#FAA371", "#FC5D6096", "#9BC799B9", "#94CCF2AC")
    names(CPCOLS) <- lab
    backs <- data.frame(
      x1 = c(0, 0.05, 0.62, 0.8, 0, 0.3, 0.75),
      x2 = c(0.05, 0.62, 0.8, 1, 0.3, 0.75, 1),
      y1 = c(-Inf, -Inf, -Inf, -Inf, 2.5, 2.5, 2.5),
      y2 = c(2.5, 2.5, 2.5, 2.5, Inf, Inf, Inf),
      lab = factor(lab, levels = lab)
    )
  } else if (mode == 2) {
    lab <- c("Peripherals", "Network hubs", "Module hubs", "Connectors")
    CPCOLS <- c("#FCF6EFFC", "#FC5D6096", "#9BC799B9", "#94CCF2AC")
    names(CPCOLS) <- lab
    backs <- data.frame(
      x1 = c(0, 0.62, 0, 0.62),
      x2 = c(0.62, 1, 0.62, 1),
      y1 = c(-Inf, 2.5, 2.5, -Inf),
      y2 = c(2.5, Inf, Inf, 2.5),
      lab = factor(lab, levels = lab)
    )
  }
  p <- ggplot() +
    geom_rect(data = backs, mapping = aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = lab), alpha = 0.7) +
    guides(fill = guide_legend(title = "Topological roles")) +
    scale_fill_manual(values = CPCOLS) +
    geom_point(data = taxa.roles, aes(x = Pi, y = Zi, color = factor(v_class))) +
    scale_color_manual(values = setNames(unique(taxa.roles$color), unique(taxa.roles$v_class))) +
    MetaNet_theme +
    guides(colour = "none") +
    theme(strip.background = element_rect(fill = "white")) +
    xlab("Participation Coefficient") +
    ylab("Within-module connectivity z-score")

  if (label) {
    p <- p + ggrepel::geom_text_repel(
      data = taxa.roles[!taxa.roles$roles %in% c("Peripherals", "Ultra-peripherals"), ],
      aes(x = Pi, y = Zi, label = name), size = 3
    )
  }
  return(p)
}