# ========5.topological=======
#' Extract each sample network from the whole network
#'
#' @param whole_net the whole network
#' @param otutab otutab, columns are samples, these columns will be extract
#' @param threads threads, default: 1
#' @param save_net should save these sub_nets? FALSE or a filename
#' @param fast less indexes for faster calculate ?
#' @param verbose verbose
#' @param remove_negative remove negative edge or not? default: FALSE
#'
#' @return a dataframe contains all sub_net parameters
#' @export
#' @family topological
#' @examples
#' data(otutab, package = "pcutils")
#' extract_sample_net(co_net, otutab) -> sub_net_pars
extract_sample_net <- function(whole_net, otutab, threads = 1, save_net = FALSE, fast = TRUE, remove_negative = FALSE, verbose = TRUE) {
i <- NULL
V(whole_net)$name -> v_name
reps <- ncol(otutab)
if (verbose) message("extracting")
sub_nets <- lapply(1:reps, \(i){
rownames(otutab)[otutab[, i] > 0] -> exist_sp
subgraph(whole_net, which(v_name %in% exist_sp)) -> spe_sub
class(spe_sub) <- c("metanet", "igraph")
return(spe_sub)
})
names(sub_nets) <- colnames(otutab)
if (verbose) message("calculating topological indexes")
# parallel
# main function
loop <- function(i) {
spe_sub <- sub_nets[[i]]
indexs <- net_par(spe_sub, mode = "n", fast = fast, remove_negative = remove_negative)[["n_index"]]
wc <- igraph::cluster_fast_greedy(spe_sub, weights = abs(igraph::E(spe_sub)$weight))
indexs$modularity <- igraph::modularity(wc)
indexs
}
{
if (threads > 1) {
pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)
if (verbose) {
pb <- utils::txtProgressBar(max = reps, style = 3)
opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))
} else {
opts <- NULL
}
cl <- snow::makeCluster(threads)
doSNOW::registerDoSNOW(cl)
res <- foreach::`%dopar%`(
foreach::foreach(i = 1:reps, .options.snow = opts),
loop(i)
)
snow::stopCluster(cl)
gc()
} else {
res <- lapply(1:reps, loop)
}
}
# simplify method
sub_net_pars <- do.call(rbind, res)
rownames(sub_net_pars) <- colnames(otutab)
if (is.logical(save_net)) {
if (save_net) save_net <- paste0("sub_net_", date())
}
if (is.character(save_net)) {
saveRDS(sub_nets, file = paste0(save_net, ".RDS"))
}
sub_net_pars
}
#' Calculate natural_connectivity
#'
#' @param p an igraph or metanet object
#' @return natural_connectivity (numeric)
#' @export
#' @references \code{`nc` in `ggClusterNet`}
#' @family topological
#' @examples
#' igraph::make_ring(10) %>% nc()
nc <- function(p) {
adj_matrix <- as.matrix(igraph::as_adj(p, sparse = FALSE))
adj_matrix[abs(adj_matrix) != 0] <- 1
lambda <- eigen(adj_matrix, only.values = TRUE)$values
lambda <- sort(lambda, decreasing = TRUE)
lambda_sum <- 0
N <- length(lambda)
for (i in 1:N) lambda_sum <- lambda_sum + exp(lambda[i])
lambda_average <- log(lambda_sum / N, base = exp(1))
lambda_average
}
#' Calculate all topological indexes of a network
#'
#' @param go an igraph or metanet object
#' @param fast less indexes for faster calculate ?
#' @param mode calculate what? c("v", "e", "n", "all")
#' @param remove_negative remove negative edge or not? default: FALSE
#'
#' @return a 3-elements list
#' \item{n_index}{indexs of the whole network}
#' \item{v_index}{indexs of each vertex}
#' \item{e_index}{indexs of each edge}
#' @export
#' @family topological
#' @examples
#' igraph::make_graph("Walther") %>% net_par()
#' c_net_index(co_net) -> co_net_with_par
net_par <- function(go, mode = c("v", "e", "n", "all"), fast = TRUE, remove_negative = FALSE) {
from <- to <- NULL
stopifnot(is_igraph(go))
if ("all" %in% mode) mode <- c("v", "e", "n")
n_index <- NULL
v_index <- NULL
e_index <- NULL
Negative_percentage <- ifelse(!is.null(E(go)$cor), sum(igraph::E(go)$cor < 0) / length(igraph::E(go)), NA)
# remove negative weight
if (remove_negative) {
if (!is.null(E(go)$cor)) {
# message("Remove negative correlation edges")
c_net_filter(go, cor > 0, mode = "e") -> go
}
}
# non-weighted network
up <- go
if (!is.null(igraph::edge_attr(up)[["weight"]])) up <- igraph::delete_edge_attr(up, "weight")
if ("n" %in% mode) {
# Calculate Network Parameters
n_index <- data.frame(
check.names = F,
`Node_number` = length(igraph::V(go)), # number of nodes
`Edge_number` = length(igraph::E(go)), # number of edges
`Edge_density` = igraph::edge_density(go), # density of network, connectance
`Negative_percentage` = Negative_percentage, # negative edges percentage
`Average_path_length` = igraph::average.path.length(up), # Average path length
`Global_efficiency` = igraph::global_efficiency(up),
`Average_degree` = mean(igraph::degree(go)), # Average degree
`Average_weighted_degree` = ifelse(is.null(igraph::E(go)$weight), mean(igraph::degree(go)), sum(igraph::E(go)$weight) / length(igraph::V(go))), # weighted degree
Diameter = igraph::diameter(up), # network diameter
`Clustering_coefficient` = igraph::transitivity(go), # Clustering coefficient
`Centralized_betweenness` = igraph::centralization.betweenness(go)$centralization, # Betweenness centralization
`Natural_connectivity` = nc(go) # natural
)
if (!fast) {
# mean_dist=mean_distance(go)#
# w_mean_dist=ifelse(is.null(E(go)$weight),mean_dist,mean_distance(go))
# v_conn= vertex.connectivity(go) #
# e_conn= edge.connectivity(go) #
# components= count_components(go) #
modularity <- igraph::modularity(igraph::cluster_fast_greedy(go)) #
rand.g <- igraph::erdos.renyi.game(length(V(go)), length(E(go)), type = "gnm")
rand_m <- igraph::modularity(igraph::cluster_fast_greedy(rand.g))
relative_modularity <- (modularity - rand_m) / rand_m #
n_index <- data.frame(
check.names = F,
n_index,
Modularity = modularity,
`Relative_modularity` = relative_modularity,
`Centralized_closeness` = igraph::centralization.closeness(go)$centralization, # Closeness centralization
`Centralized_degree` = igraph::centralization.degree(go)$centralization, # Degree centralization
`Centralized_eigenvector` = igraph::centralization.evcent(go)$centralization # eigenvector centralization
)
}
n_index <- apply(n_index, 1, FUN = \(x)replace(x, is.nan(x), 0)) %>%
t() %>%
as.data.frame()
n_index <- cbind_new(get_n(go, simple = TRUE), n_index)
}
if ("v" %in% mode) {
# Calculate Vertices Parameters
v_index <- data.frame(
check.names = F,
Degree = igraph::degree(go),
`Clustering_coefficient` = igraph::transitivity(go, type = "local"), # local clustering coefficient
Betweenness = igraph::betweenness(go), # betweenness
Eccentricity = igraph::eccentricity(go),
Closeness = igraph::closeness(go),
`Hub_score` = igraph::hub_score(go)[["vector"]]
# page_rank = page.rank(go)$vector
# igraph::evcent(go)[["vector"]]
# igraph::local_efficiency(go)
)
# weighted degree
if (!is.null(E(go)$cor)) {
get_e(go) -> edge_list
edge_list %>%
dplyr::select(from, cor) %>%
rbind(., dplyr::select(edge_list, to, cor) %>% dplyr::rename(from = to)) %>%
dplyr::group_by(from) %>%
dplyr::summarise(w_degree = sum(cor)) -> w_degree
v_index$`Average_weighted_degree` <- w_degree[match(rownames(v_index), w_degree$from), "w_degree"] %>% unlist()
}
v_index <- apply(v_index, 1, FUN = \(x)replace(x, is.nan(x), 0)) %>%
t() %>%
as.data.frame()
v_index <- cbind_new(get_v(go), v_index)
}
if ("e" %in% mode) {
# Calculate Edges Parameters
e_index <- get_e(go)
# if(!(edge_attr(go)%>%unlist()%>%is.null()))e_index=data.frame(edge_attr(go),e_index)
}
return(list(n_index = n_index, v_index = v_index, e_index = e_index))
}
#' Add topological indexes for a network
#' @param go igraph or metanet
#' @param force replace existed net_par
#'
#' @export
#' @rdname net_par
c_net_index <- function(go, force = FALSE) {
if (!force) {
if (!is.null(graph_attr(go)[["net_par"]])) stop("Already calculated net_pars, set `force = TRUE to replace existed net_par")
}
net_par(go, fast = FALSE) -> res
graph_attr(go) <- as.list(res$n_index)
graph_attr(go)[["net_par"]] <- TRUE
vertex_attr(go) <- as.list(res$v_index)
edge_attr(go) <- as.list(res$e_index)
go
}
#' Fit power-law distribution for an igraph
#'
#' @param go igraph
#' @param p.value calculate p.value
#'
#' @return ggplot
#' @export
#' @family topological
#' @examples
#' fit_power(co_net)
fit_power <- function(go, p.value = FALSE) {
x <- y <- formula <- NULL
# igraph::degree distribution
degree_dist <- table(igraph::degree(go))
dat <- data.frame(degree = as.numeric(names(degree_dist)), count = as.numeric(degree_dist))
# fit, set the original a & b
mod <- stats::nls(count ~ a * degree^b, data = dat, start = list(a = 2, b = 1.5))
summary(mod)
# extract the coefficient
a <- round(coef(mod)[1], 3)
b <- round(coef(mod)[2], 3)
fit <- fitted(mod)
SSre <- sum((dat$count - fit)^2)
SStot <- sum((dat$count - mean(dat$count))^2)
R2 <- round(1 - SSre / SStot, 3)
# bootstrap t get p.value
if (p.value) {
dat_rand <- dat
p_num <- lapply(seq_len(999), \(i){
dat_rand$count <- sample(dat_rand$count)
SSre_rand <- sum((dat_rand$count - fit)^2)
SStot_rand <- sum((dat_rand$count - mean(dat_rand$count))^2)
R2_rand <- 1 - SSre_rand / SStot_rand
R2_rand > R2
})
p_value <- (sum(unlist(p_num)) + 1) / (999 + 1)
}
p <- ggplot(dat, aes(x = degree, y = count)) +
geom_point() +
theme_bw() +
stat_smooth(method = "nls", formula = y ~ a * x^b, method.args = list(start = list(a = 2, b = 1.5)), se = FALSE) +
labs(x = "Degree", y = "Count")
if (p.value) {
label <- data.frame(
x = 0.8 * max(dat$degree),
y = c(0.9, 0.8, 0.7) * max(dat$count),
formula = c(
sprintf("italic(Y) == %.3f*italic(X)^%.3f", a, b),
sprintf("italic(R^2) == %.3f", R2),
sprintf("italic(P) < %.3f", p_value)
)
)
} else {
label <- data.frame(
x = 0.8 * max(dat$degree),
y = c(0.9, 0.8) * max(dat$count),
formula = c(
sprintf("italic(Y) == %.3f*italic(X)^%.3f", a, b),
sprintf("italic(R^2) == %.3f", R2)
)
)
}
p + geom_text(aes(x = x, y = y, label = formula), data = label, parse = TRUE)
}
#' Degree distribution comparison with random network
#'
#' @param go igraph object
#' @param plot plot or not
#'
#' @return ggplot
#' @export
#' @family topological
#' @examples
#' rand_net(co_net)
rand_net <- function(go = go, plot = TRUE) {
freq <- net <- NULL
# generate a random network
rand.g <- igraph::erdos.renyi.game(length(V(go)), length(E(go)), type = "gnm")
if (!plot) {
return(rand.g)
}
data1 <- data.frame(
freq = igraph::degree_distribution(go), net = "Network",
degree = 0:(length(degree_distribution(go)) - 1)
)
data2 <- data.frame(
freq = igraph::degree_distribution(rand.g), net = "Random E-R",
degree = 0:(length(degree_distribution(rand.g)) - 1)
)
# if data1[1,1]=0, it'is delete single vertex
if (data1[1, 1] == 0) data1 <- data1[-1, ]
data <- rbind(data1, data2)
p1 <- ggplot(data) +
geom_point(aes(x = degree, y = freq, group = net, fill = net), pch = 21, size = 2) +
geom_smooth(aes(x = degree, y = freq, group = net, color = net), se = FALSE, method = "loess", formula = "y ~ x") +
labs(x = "Degree", y = "Proportion") +
scale_color_manual(values = c("#F58B8B", "#7AADF0")) +
scale_fill_manual(values = c("#F58B8B", "#7AADF0")) +
MetaNet_theme +
theme(legend.position = c(0.8, 0.9), legend.title = element_blank())
print(p1)
return(rand.g)
}
#' Net_pars of many random network
#'
#' @param go igraph
#' @param reps simulation time
#' @param threads threads
#' @param verbose verbose
#'
#' @export
#' @rdname compare_rand
rand_net_par <- function(go, reps = 99, threads = 1, verbose = TRUE) {
i <- NULL
# parallel
# main function
loop <- function(i) {
# generate a random network
rand.g <- igraph::erdos.renyi.game(length(igraph::V(go)),
length(igraph::E(go)),
type = "gnm"
)
indexs <- net_par(rand.g, mode = "n")[["n_index"]]
wc <- igraph::cluster_fast_greedy(rand.g)
indexs$modularity <- igraph::modularity(wc)
indexs
}
{
if (threads > 1) {
pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)
if (verbose) {
pb <- utils::txtProgressBar(max = reps, style = 3)
opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))
} else {
opts <- NULL
}
cl <- snow::makeCluster(threads)
doSNOW::registerDoSNOW(cl)
res <- foreach::`%dopar%`(
foreach::foreach(i = 1:reps, .options.snow = opts),
loop(i)
)
snow::stopCluster(cl)
gc()
} else {
res <- lapply(1:reps, loop)
}
}
# simplify method
rand_net_pars <- do.call(rbind, res)
rand_net_pars
}
#' Compare some indexes between your net with random networks
#'
#' @param pars your net pars resulted by net_pars()
#' @param randp random networks pars resulted by rand_net_par()
#' @param index compared indexes: "Average_path_length","Clustering_coefficient" or else
#'
#' @return ggplot
#' @export
#' @family topological
#' @examples
#' data("c_net")
#' rand_net_par(co_net_rmt, reps = 30) -> randp
#' net_par(co_net_rmt, fast = FALSE) -> pars
#' compare_rand(pars, randp)
compare_rand <- function(pars, randp, index = c("Average_path_length", "Clustering_coefficient")) {
V1 <- NULL
labss <- t(pars$n_index[, index, drop = FALSE]) %>% as.data.frame()
rownames(labss) -> labss$indexes
p <- pcutils::group_box(randp[, index, drop = FALSE])
p <- p +
geom_hline(data = labss, aes(yintercept = V1), linetype = 2, color = "blue3") +
geom_text(
data = labss, aes(x = 1, y = V1 * 1.05, label = paste0("Network: ", round(V1, 3))),
color = "blue3"
) +
MetaNet_theme +
theme(legend.position = "none", axis.text.x = element_blank())
p
}
#' Calculate small-world coefficient
#'
#' @param go igraph or metanet
#' @param reps simulation time
#' @param threads threads
#' @param verbose verbose
#'
#' @return number
#' @export
#' @family topological
#' @examples
#' \donttest{
#' # set reps at least 99 when you run.
#' smallworldness(co_net, reps = 9)
#' }
smallworldness <- function(go, reps = 99, threads = 1, verbose = TRUE) {
rand_net_par(go, reps = reps, threads = threads, verbose = verbose) -> rands
small_world_coefficient <- (igraph::transitivity(go) / mean(rands$Clustering_coefficient)) /
(igraph::average.path.length(go) / mean(rands$`Average_path_length`))
small_world_coefficient
}
Datasets
Models
