Diff of /R/2-1.build.R [000000] .. [13df9a]

Switch to side-by-side view

--- a
+++ b/R/2-1.build.R
@@ -0,0 +1,460 @@
+# =========2.build======
+
+#' Construct a metanet from a corr object
+#'
+#' @param corr corr object from `c_net_calculate()` or `read_corr()`.
+#' @param r_threshold r_threshold (default: >0.6).
+#' @param p_threshold p_threshold (default: <0.05).
+#' @param use_p_adj use the p.adjust instead of p.value (default: TRUE), if p.adjust not in the corr object, use p.value.
+#' @param delete_single should delete single vertexes?
+#'
+#' @return an metanet object
+#' @export
+#'
+#' @family build
+#' @examples
+#' data("otutab", package = "pcutils")
+#' t(otutab) -> totu
+#' metadata[, 3:10] -> env
+#' c_net_calculate(totu) -> corr
+#' c_net_build(corr, r_threshold = 0.65) -> co_net
+#'
+#' c_net_calculate(totu, env) -> corr2
+#' c_net_build(corr2) -> co_net2
+c_net_build <- function(corr, r_threshold = 0.6, p_threshold = 0.05, use_p_adj = TRUE, delete_single = TRUE) {
+  stopifnot(inherits(corr, "corr"))
+  # set thresholds to construct
+  if ("r" %in% names(corr)) {
+    occor.r <- corr$r
+  } else {
+    stop("No r in the input corr object.")
+  }
+
+  if (("p.adjust" %in% names(corr)) && use_p_adj) {
+    occor.p <- corr$p.adjust
+  } else {
+    if ("p.value" %in% names(corr)) {
+      occor.p <- corr$p.value
+      message("Have not do p-value adjustment! use the p.value to build network.")
+    } else {
+      occor.p <- occor.r
+      occor.p[occor.p != 0] <- 0
+      message("No p.value or p.adjust in the input corr object, use the r to build network.")
+    }
+  }
+
+  if (any(is.na(occor.r))) {
+    occor.r[is.na(occor.r)] <- 0
+    warning("NA values in the input corr object, set to 0.")
+  }
+  occor.r[occor.p > p_threshold | abs(occor.r) < r_threshold] <- 0
+
+  # make igraph
+  if (t_flag(occor.r)) {
+    go <- igraph::graph_from_adjacency_matrix(as.matrix(occor.r),
+      mode = "undirected", weighted = TRUE, diag = FALSE
+    )
+    igraph::graph.attributes(go)$n_type <- "single"
+  } else {
+    inset <- intersect(rownames(occor.r), colnames(occor.r))
+    if (length(inset) > 0) {
+      insets <- paste0(inset, collapse = ", ")
+      warning("some nodes (row and column) has same name, please check: ", insets)
+    }
+    go <- igraph::graph_from_incidence_matrix(as.matrix(occor.r),
+      directed = FALSE, weighted = TRUE
+    )
+    igraph::graph.attributes(go)$n_type <- "bipartite"
+  }
+
+  # delete single vertexes?
+  if (delete_single) go <- igraph::delete.vertices(go, V(go)[igraph::degree(go) == 0])
+
+  # set vertex attributes
+  V(go)$v_group <- ifelse(V(go)$name %in% rownames(occor.r), "v_group1", "v_group2")
+  V(go)$v_class <- ifelse(V(go)$name %in% rownames(occor.r), "v_class1", "v_class2")
+  V(go)$size <- ceiling(60 / sqrt(length(V(go)))) + 1
+  V(go)$label <- V(go)$name
+
+  # abs edges weight
+  go.weight <- E(go)$weight
+  E(go)$cor <- go.weight
+  E(go)$weight <- abs(go.weight)
+
+  # time-consuming
+  if (TRUE) {
+    tmp_e <- get_e(go)
+    if ("p.value" %in% names(corr)) {
+      # E(go)$p.value <-get_e(go)%>%dplyr::select(from,to)%>%apply(., 1, \(x)occor.r$p.value[x[1],x[2]])
+      tmp <- reshape2::melt(corr$p.value, varnames = c("from", "to"), value.name = "p.value", as.is = TRUE)
+      tmp_e <- dplyr::left_join(tmp_e, tmp, by = c("from", "to"))
+    }
+    if ("p.adjust" %in% names(corr)) {
+      # E(go)$p.adjust <-get_e(go)%>%dplyr::select(from,to)%>%apply(., 1, \(x)occor.r$p.adjust[x[1],x[2]])
+      tmp <- reshape2::melt(corr$p.adjust, varnames = c("from", "to"), value.name = "p.adjust", as.is = TRUE)
+      tmp_e <- dplyr::left_join(tmp_e, tmp, by = c("from", "to"))
+    }
+    # 应该直接expand,再left_join快很多
+    igraph::edge.attributes(go) <- as.list(tmp_e)
+  }
+
+  # set edges type
+  E(go)$e_type <- ifelse(go.weight > 0, "positive", "negative")
+  # set edges width
+  E(go)$width <- E(go)$weight
+
+  if (!t_flag(occor.r)) {
+    # set edges from_to
+    anno_edge(go, get_v(go)[, c("name", "v_group")], verbose = FALSE) -> go
+    # set edge intra-inter
+    tmp_e <- igraph::edge.attributes(go)
+    E(go)$e_class <- ifelse(tmp_e$v_group_from == tmp_e$v_group_to, "intra", "inter")
+  }
+
+  c_net_update(go, initialize = TRUE) -> go
+
+  return(go)
+}
+
+
+#' Multi-omics network build
+#'
+#' @param ... some omics abundance tables
+#' @param mode "full"
+#' @param method "spearman" or "pearson"
+#' @param filename the prefix of saved .corr file or FALSE
+#' @param p.adjust.method see \code{\link[stats]{p.adjust}}
+#' @param r_threshold r_threshold (default: >0.6)
+#' @param p_threshold p_threshold (default: <0.05)
+#' @param use_p_adj use the p.adjust instead of p-value (default: TRUE)
+#' @param delete_single should delete single vertexes?
+#'
+#' @return metanet
+#' @export
+#' @family build
+#' @examples
+#' data("multi_test")
+#' multi1 <- multi_net_build(list(Microbiome = micro, Metabolome = metab, Transcriptome = transc))
+#' multi1 <- c_net_set(multi1, micro_g, metab_g, transc_g,
+#'   vertex_class = c("Phylum", "kingdom", "type")
+#' )
+#' multi1 <- c_net_set(multi1, data.frame("Abundance1" = colSums(micro)),
+#'   data.frame("Abundance2" = colSums(metab)), data.frame("Abundance3" = colSums(transc)),
+#'   vertex_size = paste0("Abundance", 1:3)
+#' )
+#' c_net_plot(multi1)
+multi_net_build <- function(..., mode = "full", method = "spearman",
+                            filename = FALSE, p.adjust.method = NULL,
+                            r_threshold = 0.6, p_threshold = 0.05, use_p_adj = TRUE, delete_single = TRUE) {
+  tables <- list(...)
+  if (all(class(tables[[1]]) == "list")) tables <- tables[[1]]
+  tabs_name <- names(tables)
+
+  tables <- check_tabs(tables)
+  if (mode == "full") {
+    all_totu <- do.call(cbind, tables)
+    message("Calculating ", nrow(all_totu), " samples and ", ncol(all_totu), " features of ", length(tables), " groups.\n")
+    all_corr <- c_net_calculate(all_totu, method = method, filename = filename, p.adjust.method = p.adjust.method)
+
+    c_net_build(all_corr,
+      r_threshold = r_threshold, p_threshold = p_threshold,
+      use_p_adj = use_p_adj, delete_single = delete_single
+    ) -> multi_net
+
+    igraph::graph.attributes(multi_net)$n_type <- "multi_full"
+
+    get_v(multi_net) -> tmp_v
+    if (is.null(tabs_name)) tabs_name <- paste0("omic", seq_len(length(tables)))
+    position <- rep(tabs_name, vapply(tables, ncol, numeric(1)))
+    names(position) <- lapply(tables, colnames) %>% do.call(c, .)
+
+    tmp_v$v_class <- tmp_v$v_group <- vapply(tmp_v$name, \(x)position[x], character(1))
+    tmp_v %>% as.list() -> igraph::vertex_attr(multi_net)
+
+    # set edges from_to
+    suppressMessages(anno_edge(multi_net, tmp_v[, c("name", "v_group")], verbose = FALSE) -> multi_net)
+    # set edge intra-inter
+    tmp_e <- igraph::edge.attributes(multi_net)
+    E(multi_net)$e_class <- ifelse(tmp_e$v_group_from == tmp_e$v_group_to, "intra", "inter")
+    c_net_update(multi_net, initialize = TRUE, verbose = FALSE) -> multi_net
+  }
+  multi_net
+}
+
+
+default_v_color <- c(
+  "#a6bce3", "#fb9a99", "#fdbf6f", "#1f78b4", "#b2df8a",
+  "#cab2d6", "#33a02c", "#e31a1c", "#ff7f00", "#6a3d9a",
+  "#F8CC00", "#b15928"
+)
+default_v_color_num <- c(
+  "#FEE0D2", "#FCBBA1", "#FC9272", "#FB6A4A", "#EF3B2C", "#CB181D",
+  "#A50F15", "#67000D"
+)
+default_e_color <- c(
+  "#a6cee3", "#78c679", "#c2a5cf", "#ff7f00", "#1f78b4",
+  "#810f7c", "#F8CC00", "#006d2c", "#4d4d4d", "#8c510a",
+  "#d73027", "#7f0000", "#41b6c4", "#e7298a", "#54278f"
+)
+default_e_color_num <- c(
+  "#E5F5E0", "#C7E9C0", "#A1D99B", "#74C476", "#41AB5D", "#238B45",
+  "#006D2C", "#00441B"
+)
+
+default_v_shape <- c("circle" = 21, "square" = 22, "diamond" = 23, "triangle1" = 24, "triangle2" = 25) # "triangle", "diamond", "hexagon", "star"
+
+color_generate <- function(v_class, n_break = 5, mode = "v",
+                           v_color = default_v_color, v_color_num = default_v_color_num,
+                           e_color = default_e_color, e_color_num = default_e_color_num) {
+  # !!!现在的legend数字不好看,考虑pretty改进
+  if (is.numeric(v_class)) {
+    # n_color=n_break
+    # n_label=seq(min(v_class, na.rm = TRUE), max(v_class, na.rm = TRUE), length.out = n_break)
+
+    n_break <- pretty(v_class, n_break)
+    n_color <- length(n_break) - 1
+    n_label <- n_break[seq_len(n_color)] %>% format(., trim = TRUE)
+
+    if (mode == "v") {
+      return(as.character(cut(v_class,
+        breaks = n_break,
+        labels = pcutils::get_cols(n_color, v_color_num)
+      )))
+    } else if (mode == "e") {
+      return(as.character(cut(v_class,
+        breaks = n_break,
+        labels = pcutils::get_cols(n_color, e_color_num)
+      )))
+    } else if (mode == "label") {
+      return(as.character(cut(v_class,
+        breaks = n_break,
+        labels = n_label
+      )))
+    }
+  } else {
+    if (mode == "v") {
+      return(pcutils::tidai(v_class, pcutils::get_cols(nlevels(factor(v_class)), v_color), fac = TRUE))
+    } else if (mode == "e") {
+      v_class <- droplevels(as.factor(v_class))
+      if (all(levels(v_class) %in% c("negative", "positive"))) {
+        ncols <- c(negative = "#E85D5D", positive = "#48A4F0")[levels(v_class)]
+      } else if (all(levels(v_class) %in% c("inter-module", "intra-module"))) {
+        ncols <- c("inter-module" = "#FA789A", "intra-module" = "#A6CEE3")[levels(v_class)]
+      } else {
+        ncols <- e_color
+      }
+      return(pcutils::tidai(v_class, pcutils::get_cols(nlevels(factor(v_class)), ncols), fac = TRUE))
+    } else if (mode == "label") {
+      return(v_class)
+    }
+  }
+}
+
+#' Update a metanet object or transform igraph object to metanet object
+#'
+#' @param go a metanet object or igraph object
+#' @param node_break node_break if v_class is numeric, default: 5
+#' @param initialize initialize?
+#' @param edge_break edge_break if e_type is numeric, default: 5
+#' @param verbose verbose?
+#' @param uniq_v_class if TRUE, add prefix to v_class if multiple v_class belong to same v_group.
+#'
+#' @aliases as.metanet
+#' @export
+#' @family build
+#' @return metanet
+c_net_update <- function(go, node_break = 5, edge_break = 5, initialize = FALSE, verbose = TRUE, uniq_v_class = FALSE) {
+  stopifnot(inherits(go, "igraph"))
+  if (length(go) == 0) {
+    message("Empty graph, return empty metanet object.")
+    return(go)
+  }
+  # name
+  if (!"name" %in% igraph::vertex_attr_names(go)) V(go)$name <- as.character(seq_len(length(go)))
+  if (!"label" %in% igraph::vertex_attr_names(go)) V(go)$label <- V(go)$name
+  get_v(go) -> tmp_v
+  # v_size
+  if (!"size" %in% colnames(tmp_v)) tmp_v$size <- 1
+
+  # !!!从这里开始仔细修改,考虑对应关系
+  # v_shape
+  if (!"v_group" %in% colnames(tmp_v)) {
+    tmp_v$v_group <- "v_group1"
+  } else {
+    tmp_v$v_group <- as.character(tmp_v$v_group)
+    if (length(unique(tmp_v$v_group)) > 10) warning("Too many groups, please check the 'v_group'.")
+  }
+  flag <- TRUE
+  if ("shape" %in% colnames(tmp_v)) {
+    if (initialize) {
+      if (e_match(tmp_v[, c("v_group", "shape")])) {
+        flag <- FALSE
+      } else if (verbose) message("'v_group' and 'shape' not match one by one, update 'shape'.")
+    } else {
+      if (e_match(tmp_v[, c("v_group", "shape")], test = 1)) {
+        flag <- FALSE
+      } else if (verbose) message("a 'v_group' should not match multiple shapes, update 'shape'.")
+    }
+  }
+  if (flag) tmp_v$shape <- pcutils::tidai(tmp_v$v_group, names(default_v_shape), fac = TRUE)
+
+  # v_color
+  if (!"v_class" %in% colnames(tmp_v)) {
+    tmp_v$v_class <- "v_class1"
+  }
+
+  if (uniq_v_class) {
+    # 如果有同样的v_class对应不同v_group,则加上前缀
+    dplyr::distinct_all(tmp_v[, c("v_class", "v_group")]) %>% dplyr::count(v_class) -> group_class
+    group_class <- setNames(group_class$n, group_class$v_class)
+    if (any(group_class > 1)) {
+      la <- which(group_class > 1)
+      if (verbose) message("some 'v_class' belong to multiple 'v_group': \n", paste0(names(la), collapse = ", "))
+      tmp_v$v_class %in% names(la) -> r_index
+      tmp_v$v_class[r_index] <- paste0(tmp_v$v_group[r_index], "-", tmp_v$v_class[r_index])
+    }
+  }
+
+  flag <- TRUE
+  if ("color" %in% colnames(tmp_v)) {
+    if (initialize) {
+      if (e_match(tmp_v[, c("v_class", "color")])) {
+        flag <- FALSE
+      } else if (verbose) message("'v_class' and 'color' not match one by one, update 'color'.")
+    } else {
+      if (e_match(tmp_v[, c("v_class", "color")], test = 1)) {
+        flag <- FALSE
+      } else if (verbose) message("a 'v_class' should not match multiple colors, update 'color'.")
+    }
+  }
+
+  if (flag) {
+    tmp_v$color <- color_generate(tmp_v$v_class, node_break, mode = "v")
+    tmp_v$v_class <- color_generate(tmp_v$v_class, node_break, mode = "label")
+  }
+
+  as.list(tmp_v) -> igraph::vertex_attr(go)
+
+  # e_color
+  if (!"e_type" %in% igraph::edge_attr_names(go)) {
+    E(go)$e_type <- "e_type1"
+  }
+  flag <- TRUE
+  if ("color" %in% igraph::edge_attr_names(go)) {
+    if (initialize) {
+      if (e_match(data.frame(a = E(go)$e_type, b = E(go)$color))) {
+        flag <- FALSE
+      } else if (verbose) message("'e_type' and 'color' not match one by one, update 'color'.")
+    } else {
+      if (e_match(data.frame(a = E(go)$e_type, b = E(go)$color), test = 1)) {
+        flag <- FALSE
+      } else if (verbose) message("a 'e_type' should not match multiple colors, update 'color'.")
+    }
+  }
+  if (flag) {
+    E(go)$color <- color_generate(E(go)$e_type, edge_break, mode = "e")
+    E(go)$e_type <- color_generate(E(go)$e_type, edge_break, mode = "label")
+  }
+
+  # e_lty
+  if (!"e_class" %in% igraph::edge_attr_names(go)) {
+    E(go)$e_class <- "e_class1"
+  }
+  flag <- TRUE
+  if ("lty" %in% igraph::edge_attr_names(go)) {
+    if (initialize) {
+      if (e_match(data.frame(a = E(go)$e_class, b = E(go)$lty))) {
+        flag <- FALSE
+      } else if (verbose) message("'e_class' and 'lty' not match one by one, update 'lty'.")
+    } else {
+      if (e_match(data.frame(a = E(go)$e_class, b = E(go)$lty), test = 1)) {
+        flag <- FALSE
+      } else if (verbose) message("a 'e_class' should not match multiple linetypes, update 'lty'.")
+    }
+  }
+  if (flag) E(go)$lty <- pcutils::tidai(E(go)$e_class, 1:4, fac = TRUE)
+
+  # e_width
+  if (!"width" %in% igraph::edge_attr_names(go)) {
+    E(go)$width <- 1
+  }
+
+  class(go) <- c("metanet", "igraph")
+  return(go)
+}
+
+
+#' Clean a igraph object
+#'
+#' @param go igraph, metanet objects
+#' @param direct direct?
+#'
+#' @return a igraph object
+#' @export
+clean_igraph <- function(go, direct = TRUE) {
+  stopifnot(inherits(go, "igraph"))
+  go <- igraph::graph_from_data_frame(
+    d = get_e(go)[, c("from", "to")],
+    directed = direct,
+    vertices = get_v(go)["name"]
+  )
+  go
+}
+
+
+#' Construct a network from edge_list dataframe
+#'
+#' @param edgelist first is source, second is target, others are annotation
+#' @param vertex_df vertex metadata data.frame
+#' @param direct logical
+#' @param e_type set e_type
+#' @param e_class set e_class
+#' @return metanet
+#' @export
+#' @family build
+#' @examples
+#' data(edgelist)
+#' edge_net <- c_net_from_edgelist(arc_count, vertex_df = arc_taxonomy)
+#' edge_net <- c_net_set(edge_net, vertex_class = "Phylum", edge_width = "n")
+#' c_net_plot(edge_net)
+c_net_from_edgelist <- function(edgelist, vertex_df = NULL, direct = FALSE, e_type = NULL, e_class = NULL) {
+  if (!is.null(vertex_df)) {
+    if (!"name" %in% colnames(vertex_df)) {
+      vertex_df$name <- rownames(vertex_df)
+      message("No 'name' in the colnames(vertex_df), use rownames(vertex_df) as the 'name'.")
+    }
+    vertex_df <- data.frame(vertex_df[, "name", drop = FALSE], vertex_df[, setdiff(colnames(vertex_df), "name"), drop = FALSE])
+  }
+
+  if (!all(c("from", "to") %in% colnames(edgelist))) {
+    message("No 'from' and 'to' in the colnames(edgelist), use the first two columns as the 'from' and 'to'.")
+    colnames(edgelist)[1:2] <- c("from", "to")
+  }
+  edgelist <- data.frame(edgelist[, c("from", "to"), drop = FALSE], edgelist[, setdiff(colnames(edgelist), c("from", "to")), drop = FALSE])
+  go <- igraph::graph_from_data_frame(edgelist, directed = direct, vertices = vertex_df)
+  if (!is.null(e_type)) E(go)$e_type <- edgelist[, e_type]
+  if (!is.null(e_class)) E(go)$e_class <- edgelist[, e_class]
+  go <- c_net_update(go, initialize = TRUE)
+  go
+}
+
+update_c_net_rda <- function() {
+  data("otutab", package = "pcutils", envir = environment())
+  t(otutab) -> totu
+  metadata[, 3:10] -> env
+  c_net_calculate(totu) -> corr
+  c_net_build(corr, r_threshold = 0.65) -> co_net
+  c_net_build(corr, r_threshold = 0.69) -> co_net_rmt
+
+  c_net_calculate(totu, env) -> corr2
+  c_net_build(corr2) -> co_net2
+
+  co_net <- c_net_set(co_net, taxonomy, data.frame("Abundance" = colSums(totu)),
+    vertex_class = "Phylum", vertex_size = "Abundance"
+  )
+  co_net2 <- c_net_set(co_net2, taxonomy, data.frame(name = colnames(env), env = colnames(env)),
+    vertex_class = c("Phylum", "env")
+  )
+  co_net2 <- c_net_set(co_net2, data.frame("Abundance" = colSums(totu)), vertex_size = "Abundance")
+  save(co_net, co_net2, co_net_rmt, file = "data/c_net.rda")
+}