# library(r4projects)
# setwd(get_project_wd())
# source("R/enrich_pathways.R")
# source("R/7_utils.R")
#
# setwd("demo_data/2_smartd_project")
# load("peak_marker")
# load("urine_metabolomics_data.rda")
# # load("hmdb_compound_ms1.rda")
# load("kegg_compound_ms1.rda")
# # load("metabolic_network.rda")
# # edge_data <- tidygraph::as_tibble(metabolic_network, active = "edges")
# # node_data <- tidygraph::as_tibble(metabolic_network, active = "nodes")
# #
# # node_data <-
# # node_data %>%
# # dplyr::filter(!is.na(KEGG_ID))
# #
# # edge_data <-
# # edge_data %>%
# # dplyr::filter(!is.na(from_compound_KEGG_ID) &
# # !is.na(to_compound_KEGG_ID)) %>%
# # dplyr::filter(
# # from_compound_KEGG_ID %in% node_data$KEGG_ID &
# # to_compound_KEGG_ID %in% node_data$KEGG_ID
# # )
# #
# # node_data <-
# # node_data %>%
# # dplyr::filter(KEGG_ID %in% unique(
# # c(
# # edge_data$from_compound_KEGG_ID,
# # edge_data$to_compound_KEGG_ID
# # )
# # ))
# #
# # node_data <-
# # node_data %>%
# # dplyr::mutate(name = KEGG_ID) %>%
# # dplyr::distinct(name, .keep_all = TRUE)
# #
# # edge_data <-
# # edge_data %>%
# # dplyr::mutate(from = from_compound_KEGG_ID, to = to_compound_KEGG_ID) %>%
# # dplyr::filter(from %in% node_data$name &
# # to %in% node_data$name) %>%
# # dplyr::distinct(from, to, .keep_all = TRUE)
# #
# # node_data <-
# # node_data %>%
# # dplyr::filter(name %in% unique(c(edge_data$from, edge_data$to)))
# #
# # metabolic_network <-
# # tidygraph::tbl_graph(nodes = node_data,
# # edges = edge_data,
# # directed = FALSE)
# #
# # save(metabolic_network, file = "metabolic_network.rda")
# # load("metabolic_network.rda")
# #
# ###remove duplicated edges
# # node_data <- tidygraph::as_tibble(metabolic_network, active = "nodes")
# # edge_data <- tidygraph::as_tibble(metabolic_network, active = "edges")
# # edge_data <-
# # edge_data %>%
# # dplyr::mutate(from = from_compound_KEGG_ID, to = to_compound_KEGG_ID) %>%
# # dplyr::filter(from != to) %>%
# # dplyr::distinct(from, to, .keep_all = TRUE)
# #
# # edge_data <-
# # edge_data %>%
# # dplyr::mutate(rpair_id = apply(cbind(from, to), 1, function(x)
# # paste(sort(as.character(
# # x
# # )), collapse = "_"))) %>%
# # dplyr::distinct(rpair_id, .keep_all = TRUE)
# #
# #
# # metabolic_network <-
# # tidygraph::tbl_graph(nodes = node_data,
# # edges = edge_data,
# # directed = FALSE)
# #
# # save(metabolic_network, file = "metabolic_network.rda")
# load("metabolic_network.rda")
# # load("hmdb_pathway.rda")
# load("kegg_hsa_pathway.rda")
#
# # ###for HMDB, only remain the primary_pathway
# # hmdb_pathway <-
# # hmdb_pathway %>%
# # dplyr::filter(stringr::str_detect(pathway_class, "primary_pathway"))
#
# library(tidymass)
# library(tidyverse)
# library(igraph)
# library(ggraph)
#
# up_marker <-
# peak_marker %>%
# dplyr::filter(score > 0 & correlation1 > 0.3 & p1_adj < 0.05) %>%
# dplyr::select("Gene ID", score, p1_adj)
#
# down_marker <-
# peak_marker %>%
# dplyr::filter(score < 0 &
# correlation1 < -0.3 & p1_adj < 0.05) %>%
# dplyr::select("Gene ID", score, p1_adj)
#
# up_marker <-
# up_marker %>%
# dplyr::left_join(urine_metabolomics_data@variable_info,
# by = c("Gene ID" = "variable_id"))
#
# down_marker <-
# down_marker %>%
# dplyr::left_join(urine_metabolomics_data@variable_info,
# by = c("Gene ID" = "variable_id"))
#
# feature_table_marker_up <-
# up_marker %>%
# dplyr::select("Gene ID", mz, rt, score, p1_adj) %>%
# dplyr::rename(variable_id = "Gene ID",
# degree = score,
# p_value = p1_adj) %>%
# dplyr::mutate(polarity = case_when(
# stringr::str_detect(variable_id, "POS") ~ "positive",
# stringr::str_detect(variable_id, "NEG") ~ "negative"
# ))
#
# feature_table_marker_down <-
# down_marker %>%
# dplyr::select("Gene ID", mz, rt, score, p1_adj) %>%
# dplyr::rename(variable_id = "Gene ID",
# degree = score,
# p_value = p1_adj) %>%
# dplyr::mutate(polarity = case_when(
# stringr::str_detect(variable_id, "POS") ~ "positive",
# stringr::str_detect(variable_id, "NEG") ~ "negative"
# ))
#
# feature_table_marker <-
# rbind(feature_table_marker_up, feature_table_marker_down)
#
# feature_table_all <-
# urine_metabolomics_data@variable_info %>%
# dplyr::mutate(polarity = case_when(
# stringr::str_detect(variable_id, "POS") ~ "positive",
# stringr::str_detect(variable_id, "NEG") ~ "negative"
# )) %>%
# dplyr::select(variable_id, mz, rt, polarity)
#
# metabolite_database <-
# kegg_compound_ms1
#
# pathway_database <-
# kegg_hsa_pathway
#
# column <- "rp"
# adduct.table <- NULL
# adduct.table = NULL
# ms1.match.ppm = 15
# rt.match.tol = 5
# mz.ppm.thr = 400
# threads = 3
# include_hidden_metabolites = FALSE
#
# fpa_result <-
# perform_fpa(
# feature_table_marker = feature_table_marker,
# feature_table_all = feature_table_all,
# metabolite_database = metabolite_database,
# column = column,
# adduct.table = adduct.table,
# ms1.match.ppm = ms1.match.ppm,
# rt.match.tol = rt.match.tol,
# mz.ppm.thr = mz.ppm.thr,
# threads = threads,
# include_hidden_metabolites = include_hidden_metabolites,
# metabolic_network = metabolic_network,
# pathway_database = pathway_database
# )
#
#
# #####network visualization
# result$enriched_pathways
#
# plot_metabolic_network_fpa(
# fpa_result,
# feature_table_marker = feature_table_marker,
# include_feature = FALSE,
# include_hidden_metabolites = FALSE,
# add_compound_name = TRUE
# )
#
# plot_metabolic_module_fpa(
# fpa_result,
# feature_table_marker = feature_table_marker,
# include_feature = FALSE,
# include_hidden_metabolites = FALSE,
# add_compound_name = TRUE,
# metabolic_module_index = 6
# )
#
#' Perform Feature-based Pathway Analysis (FPA)
#'
#' This function performs feature-based pathway analysis. It contains
#' metabolite annotation based on ms1, network-based metabolite filtering,
#' metabolic module identification.
#'
#' @param feature_table_marker A data frame of marker features.
#' @param feature_table_all A data frame containing all features.
#' @param metabolite_database A metabolite database object containing metabolite spectra.
#' @param column Character vector specifying the chromatography column type. Default is `c("rp", "hilic")`.
#' @param adduct.table Optional data frame containing adduct information. Default is `NULL`.
#' @param ms1.match.ppm Numeric; mass tolerance for MS1 matching in parts per million (ppm). Default is `10`.
#' @param rt.match.tol Numeric; retention time tolerance threshold in seconds. Default is `5`.
#' @param mz.ppm.thr Numeric; M/Z tolerance threshold for filtering in ppm. Default is `400`.
#' @param threads Integer; number of threads to use for parallel processing. Default is `3`.
#' @param include_hidden_metabolites Logical; whether to include hidden metabolites. Default is `FALSE`.
#' @param metabolic_network A metabolic network object representing metabolic reactions.
#' @param pathway_database A pathway database object containing metabolic pathways.
#'
#' @return A list containing:
#' - `metabolic_module_result`: A data frame with metabolic module information.
#' - `final_metabolic_network`: A `tbl_graph` object representing the final metabolic network.
#' - `annotation_table_final`: A data frame with final annotated features.
#'
#' @export
#'
<-
(feature_table_marker,
feature_table_all,
metabolite_database,
column = ("rp", "hilic"),
adduct.table = ,
ms1.match.ppm = 25,
rt.match.tol = 5,
##unit is second
mz.ppm.thr = 400,
threads = 3,
include_hidden_metabolites = ,
metabolic_network,
pathway_database) {
###check the feature_table_marker and feature_table_all
check_feature_table_marker(feature_table_marker)
check_feature_table_all(feature_table_all)
###metabolite annotation
("Annotating metabolites for feature table...\n")
# Extract metabolic network data: nodes (metabolites) and edges (connections)
edge_data <-
tidygraph::as_tibble(metabolic_network, active = "edges")
node_data <-
tidygraph::as_tibble(metabolic_network, active = "nodes")
# Remove entries with missing KEGG IDs from the metabolite database
metabolite_database@spectra.info <-
metabolite_database@spectra.info %>%
dplyr::(!(KEGG.ID)) %>%
dplyr::(KEGG.ID %in% node_data$KEGG_ID)
# Perform metabolite annotation based on ms1
annotation_table_all <-
(
feature_table = feature_table_all,
metabolite_database = metabolite_database,
column = column,
adduct.table = adduct.table,
ms1.match.ppm = ms1.match.ppm,
rt.match.tol = rt.match.tol,
mz.ppm.thr = mz.ppm.thr,
threads = threads
)
# Remove pathways that contain no compound
remain_idx <-
pathway_database@ %>%
purrr::map((x) {
(x)
}) %>%
() %>%
`>`(0) %>%
()
pathway_database <-
(object = pathway_database, remain_idx = remain_idx)
# Remove pathway databases that have no pathways left after filtering
if ((pathway_database) == 0) {
("All pathway database length is zero.\n")
}
# Combine peaks into metabolite classes based on retention time similarity
####
("Identifying metablite classes...\n")
annotation_table_all <-
(annotation_table_all$Lab.ID) %>%
furrr::future_map(
.f = (temp_id) {
x = annotation_table_all %>%
dplyr::(Lab.ID == temp_id)
x <- x %>%
dplyr::mutate( = ()) %>%
dplyr::()
rt_class <-
( = x$, rt.tol = rt.match.tol) %>%
dplyr::()
rt_class <- (x$Lab.ID[1], rt_class$, sep = "_")
x <-
(x,
compound_class = rt_class,
stringsAsFactors = )
x <-
(x$compound_class) %>%
purrr::map((y) {
z <-
x[x$compound_class == y, , = ]
score <- (z)
z <- (z, score, stringsAsFactors = )
z
}) %>%
(, .) %>%
()
x
},
.progress =
) %>%
(, .) %>%
()
# ####if our method can be consistent with standard method
# annotation_standard <-
# urine_metabolomics_data@variable_info %>%
# dplyr::filter(!is.na(SS)) %>%
# dplyr::filter(!is.na(KEGG.ID)) %>%
# dplyr::filter(KEGG.ID %in% node_data$KEGG_ID) %>%
# dplyr::filter(KEGG.ID %in% metabolite_database@spectra.info$KEGG.ID)
#
# compared_result <-
# annotation_standard %>%
# dplyr::select(variable_id, Compound.name, KEGG.ID, SS) %>%
# dplyr::left_join(annotation_table_all[, c("variable_id",
# "Compound.name",
# "KEGG.ID",
# "compound_class",
# "score")], by = "variable_id")
#
# write.csv(compared_result, "compared_result.csv")
#
# temp <-
# unique(compared_result$variable_id) %>%
# purrr::map(function(x) {
# temp <-
# compared_result %>%
# dplyr::filter(variable_id == x)
#
# accurate_score <-
# temp %>%
# dplyr::filter(KEGG.ID.x == KEGG.ID.y) %>%
# pull(score)
# max_score <- max(temp$score)
# c(accurate_score = accurate_score,
# max_score = max_score,
# diff_score = accurate_score - max_score)
# }) %>%
# do.call(rbind, .) %>%
# as.data.frame() %>%
# dplyr::filter(diff_score != 0)
#
# temp %>%
# dplyr::filter(max_score >= 80) %>%
# pull(accurate_score)
#
#
# idx = temp %>% lapply(length) %>% unlist() %>% `==`(0) %>% which()
#
# temp2 = urine_metabolomics_data@variable_info %>%
# dplyr::filter(variable_id %in% unique(compared_result$variable_id)[idx]) %>%
# dplyr::left_join(metabolite_database@spectra.info[,c("KEGG.ID", "Formula")],
# by = "KEGG.ID")
# write.csv(temp2, "temp2.csv")
#
# x = masstools::sum_formula(formula = "C9H18NO4", adduct = "M+H")
# Rdisop::getMass(molecule = Rdisop::getMolecule("C9H18NO4"))
# (156.040183 - 156.0514)*10^6/400
# Remove redundant metabolite annotations
# calculate_redundance(annotation_table = annotation_table_all)
annotation_table_all <-
(annotation_table = annotation_table_all)
# calculate_redundance(annotation_table = annotation_table_all2)
# Filter the marker feature table to include only annotated markers
feature_table_marker <-
feature_table_marker %>%
dplyr::(variable_id %in% annotation_table_all$variable_id)
annotation_table_marker <-
annotation_table_all %>%
dplyr::(variable_id %in% feature_table_marker$variable_id)
("Detecting metabolic modules from metabolic reaction network.\n")
###detected_metabolites are IDs of metabolites that matched with the metabolic features
detected_metabolites <-
(annotation_table_marker$KEGG.ID)
###hidden metabolites are metabolites that are not matched with the metabolic features,
###and are connected to the detected_metabolites in the metabolic network within three steps
("Detecting hidden metabolites...\n")
if (include_hidden_metabolites) {
hidden_metabolites <-
(
metabolic_network = metabolic_network,
detected_metabolites = detected_metabolites,
threads = threads,
max.reaction = 3
)
hidden_metabolites <-
hidden_metabolites %>%
() %>%
()
} {
hidden_metabolites <-
}
total_metabolites <-
(detected_metabolites, hidden_metabolites)
####extract subnetwork
sub_metabolic_network <-
metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% total_metabolites) %>%
dplyr::mutate(node_class = ( %in% detected_metabolites, "Detected", "Hidden")) %>%
tidygraph::activate(what = "nodes") %>%
dplyr::mutate(degree2 = tidygraph::centrality_degree())
# sub_metabolic_network2 <-
# igraph::subgraph(graph = metabolic_network,
# vids = total_metabolites) %>%
# tidygraph::as_tbl_graph() %>%
# dplyr::mutate(node_class = ifelse(name %in% detected_metabolites, "Detected", "Hidden")) %>%
# tidygraph::activate(what = "nodes") %>%
# dplyr::mutate(degree2 = tidygraph::centrality_degree())
###metabolic_module detection
#####detecting metabolic modules
("Detecting metabolic modules...\n")
metabolic_modules <-
(
sub_metabolic_network = sub_metabolic_network,
detected_metabolites = detected_metabolites,
hidden_metabolites = hidden_metabolites
)
###calculate activity scores and impacts of metabolic modules
("Calculating quality scores of metabolic modules.\n")
module_quality_score <-
(
metabolic_modules = metabolic_modules,
detected_metabolites = detected_metabolites,
hidden_metabolites = hidden_metabolites,
sub_metabolic_network = sub_metabolic_network,
threads = threads
)
# module_quality_score <-
# metabolic_modules %>%
# purrr::map(function(x) {
# calculate_module_quality(graph = sub_metabolic_network, nodes = x)$conductance
# }) %>%
# unlist()
##Module impact
("Calculating impacts of metabolic modules...\n")
module_impact <-
(metabolic_modules, (temp_group) {
temp_graph <-
igraph::subgraph(graph = sub_metabolic_network, v = temp_group)
centrality <- (graph = temp_graph, = "d")
temp_idx <- (temp_group %in% detected_metabolites)
module_impact <- (centrality[temp_idx]) / (centrality)
module_impact
})
module_impact <- (module_impact)
####Hub Dominance Score (HDS) or Dominant Edge Ratio
##DER = \frac{\sum_{\text{top } K} e_i}{|E|}
# where:
# K is a small fraction (e.g., top 10-20%) of nodes sorted by degree.
# e_i is the number of edges connected to those nodes.
# |E| is the total number of edges.
#
# Rationale: If the top nodes account for almost all edges, DER → 1.
("Calculating Dominant Edge Ratio (DER) of metabolic modules...\n")
module_dominant_edge_rate <-
(metabolic_modules, (temp_group) {
temp_graph <-
igraph::subgraph(graph = sub_metabolic_network, v = temp_group)
all_degree <- igraph::degree(temp_graph)
# Sort nodes by degree (descending)
sorted_degrees <- (all_degree, decreasing = )
# Define top K nodes (e.g., top 20% most connected nodes)
# K <- ceiling(0.2 * length(V(temp_graph))) # 20% of nodes
# top_nodes <- names(sorted_degrees)[1:K]
top_nodes <-
(all_degree[all_degree > (sorted_degrees, probs = 0.9)])
if ((top_nodes) == 0) {
top_nodes <- (sorted_degrees)[1]
}
all_edges <-
igraph::as_edgelist(temp_graph)
# Count edges where at least one node is in the top nodes
dominant_edges <-
(all_edges, 1] %in% top_nodes |
all_edges, 2] %in% top_nodes)
# Compute Dominant Edge Ratio (DER)
total_edges <- igraph::gsize(temp_graph)
dominant_edges / total_edges
}) %>%
()
##get NULL distribution
("Identify null distributation of metabolic module quality scores...\n")
null_quality_score <-
(
annotation_table_all = annotation_table_all,
feature_table_marker = feature_table_marker,
metabolite_database = metabolite_database,
metabolic_network = metabolic_network,
permutation_times = 5,
threads = threads,
adduct.table = adduct.table,
column = column,
ms1.match.ppm = ms1.match.ppm,
mz.ppm.thr = mz.ppm.thr,
include_hidden_metabolites = include_hidden_metabolites
)
# null_quality_score <-
# generate_null_conductance_distribution(
# annotation_table_all = annotation_table_all,
# feature_table_marker = feature_table_marker,
# metabolite_database = metabolite_database,
# metabolic_network = metabolic_network,
# permutation_times = 5,
# threads = threads,
# adduct.table = adduct.table,
# column = column,
# ms1.match.ppm = ms1.match.ppm,
# mz.ppm.thr = mz.ppm.thr,
# include_hidden_metabolites = include_hidden_metabolites
# )
null_quality_score <-
(null_quality_score)
# mean(module_quality_score)
# mean(null_quality_score)
# median(mean(module_quality_score))
# median(mean(null_quality_score))
####calculate the p values for each metabolic module
module_quality_score <- module_quality_score + 0.00000001
###
para <- MASS::fitdistr(x = null_quality_score, densfun = "gamma")[1]]
# standard.distribution <- rgamma(n = length(null_quality_score), shape = para[1], rate = para[2])
standard.distribution <- (n = 100000,
shape = para[1],
rate = para[2])
# ks.test(x = null_quality_score, y = standard.distribution)
cdf <- (x = standard.distribution)
metabolic_module_p <- 1 - cdf(module_quality_score)
#order module according to p value
# metabolic_modules <- metabolic_modules[order(metabolic_module_p)]
# module_impact <- module_impact[order(metabolic_module_p)]
# metabolic_module_p <- metabolic_module_p[order(metabolic_module_p)]
# module_dominant_edge_rate <- module_dominant_edge_rate[order(metabolic_module_p)]
## get information of each metabolic module
metabolic_module_result <-
(metabolic_modules, (x) {
temp_graph <-
igraph::subgraph(graph = metabolic_network, v = x)
detected_metabolite_number <- (x %in% detected_metabolites)
hidden_metabolite_number <- (x %in% hidden_metabolites)
total_metabolite_number <- (x)
i.id <- x[which(x %in% detected_metabolites)]
h.id <- x[which(x %in% hidden_metabolites)]
detected_metabolite_name <-
(KEGG.ID = i.id) %>%
dplyr::left_join(
annotation_table_all %>%
dplyr::select(KEGG.ID, Compound.name) %>%
dplyr::distinct(KEGG.ID, .keep_all = ),
= "KEGG.ID"
) %>%
pull(Compound.name) %>%
(collapse = "{}")
hidden_metabolite_name <-
(KEGG.ID = h.id) %>%
dplyr::left_join(
annotation_table_all %>%
dplyr::select(KEGG.ID, Compound.name) %>%
dplyr::distinct(KEGG.ID, .keep_all = ),
= "KEGG.ID"
) %>%
pull(Compound.name) %>%
(collapse = "{}")
total_metabolite_name <-
(KEGG.ID = x) %>%
dplyr::left_join(
annotation_table_all %>%
dplyr::select(KEGG.ID, Compound.name) %>%
dplyr::distinct(KEGG.ID, .keep_all = ),
= "KEGG.ID"
) %>%
pull(Compound.name) %>%
(collapse = "{}")
i.id <- (i.id, collapse = "{}")
h.id <- (h.id, collapse = "{}")
t.id <- (x, collapse = "{}")
temp_result <-
(
total_metabolite_number,
detected_metabolite_number,
hidden_metabolite_number,
t.id,
i.id,
h.id,
total_metabolite_name,
detected_metabolite_name,
hidden_metabolite_name
)
(temp_result)
})
metabolic_module_result <-
(, metabolic_module_result)
metabolic_module_result <-
(
module_impact,
module_dominant_edge_rate,
metabolic_module_p,
metabolic_module_result,
stringsAsFactors =
) %>%
dplyr::(metabolic_module_p) %>%
dplyr::mutate(metabolic_module_name = ("metabolic_module", 1:(.), sep = "_")) %>%
dplyr::select(metabolic_module_name, dplyr::everything())
(metabolic_module_result) <- (
"Name",
"Impact",
"Dominant_edge_rate",
"p_value",
"Total_metabolite_number",
"Detected_metabolite_number",
"Hidden_metabolite_number",
"Total_metabolite_id",
"Detected_metabolite_id",
"hidden_metabolite_id",
"Total_metabolite_name",
"Detected_metabolite_name",
"hidden_metabolite_name"
)
metabolic_module_result$Total_metabolite_number <-
(metabolic_module_result$Total_metabolite_number)
metabolic_module_result$Detected_metabolite_number <-
(metabolic_module_result$Detected_metabolite_number)
metabolic_module_result$Hidden_metabolite_number <-
(metabolic_module_result$Hidden_metabolite_number)
#####extract the final dysregulated networks
final_metabolite_id <-
metabolic_module_result %>%
dplyr::(p_value < 0.05) %>%
dplyr::pull(Total_metabolite_id) %>%
stringr::str_split( = "\\{\\}") %>%
() %>%
()
annotation_table_final <-
annotation_table_all %>%
dplyr::(
variable_id %in% feature_table_marker$variable_id &
KEGG.ID %in% final_metabolite_id
)
# ####not in the workflow
# calculate_redundance(annotation_table_final)
#
# annotation_standard <-
# urine_metabolomics_data@variable_info %>%
# dplyr::filter(!is.na(SS)) %>%
# dplyr::filter(!is.na(KEGG.ID)) %>%
# dplyr::filter(KEGG.ID %in% node_data$KEGG_ID) %>%
# dplyr::filter(KEGG.ID %in% metabolite_database@spectra.info$KEGG.ID)
#
# dim(annotation_standard)
#
# compared_result2 <-
# annotation_standard %>%
# dplyr::select(variable_id, Compound.name, KEGG.ID, SS) %>%
# dplyr::left_join(annotation_table_final[, c("variable_id",
# "Compound.name",
# "KEGG.ID",
# "compound_class",
# "score")], by = "variable_id")
#
# write.csv(compared_result2, file = "compared_result2.csv")
#
# unique(compared_result2$variable_id) %>%
# purrr::map(function(x) {
# temp <-
# compared_result2 %>%
# dplyr::filter(variable_id == x)
#
# if (all(is.na(temp$KEGG.ID.y))) {
# return("No annotation")
# }
#
# if (unique(temp$KEGG.ID.x) %in% temp$KEGG.ID.y) {
# return("Correct annotation")
# } else{
# return("Incorrect annotation")
# }
#
# }) %>%
# unlist() %>%
# table()
#
# compared_result3 <-
# compared_result2 %>%
# dplyr::filter(!is.na(KEGG.ID.y))
#
# write.csv(compared_result3, file = "compared_result3.csv")
final_metabolic_network <-
metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% final_metabolite_id) %>%
tidygraph::activate(what = "nodes") %>%
dplyr::mutate(node_class = ( %in% detected_metabolites, "Detected", "Hidden")) %>%
tidygraph::activate(what = "nodes") %>%
dplyr::mutate(degree2 = tidygraph::centrality_degree())
final_node_data <-
tidygraph::as_tibble(final_metabolic_network, active = "nodes") %>%
dplyr::select(, KEGG_ID, HMDB_ID, Compound_name, node_class)
final_edge_data <-
tidygraph::as_tibble(final_metabolic_network, active = "edges") %>%
dplyr::mutate(from = from_compound_KEGG_ID, to = to_compound_KEGG_ID) %>%
dplyr::select(from, to) %>%
dplyr::mutate(edge_class = "metabolic_reaction")
final_edge_data2 <-
annotation_table_final %>%
dplyr::(isotope == "[M]") %>%
dplyr::select(variable_id, KEGG.ID) %>%
dplyr::rename(from = variable_id, to = KEGG.ID) %>%
dplyr::mutate(edge_class = "feature_metabolite")
final_node_data2 <-
( = (final_edge_data2$from)) %>%
dplyr::mutate(
KEGG_ID = ,
HMDB_ID = ,
Compound_name = ,
node_class = "Feature"
)
###add feature information to the final metabolic network
final_node_data <-
(final_node_data, final_node_data2)
final_edge_data <-
(final_edge_data, final_edge_data2)
final_node_data <-
final_node_data %>%
dplyr::( %in% ((
final_edge_data$from, final_edge_data$to
))) %>%
dplyr::distinct(, .keep_all = )
final_edge_data <-
final_edge_data %>%
dplyr::(from %in% final_node_data$ &
to %in% final_node_data$)
final_metabolic_network <-
tidygraph::tbl_graph(nodes = final_node_data,
edges = final_edge_data,
directed = )
####Pathway enrichment for all the metabolites
kegg_id <-
annotation_table_final$KEGG.ID
("Enriching pathways for dysregulated metabolic network...\n")
enriched_pathways <-
(
query_id = kegg_id,
query_type = "compound",
id_type = "KEGG",
pathway_database = pathway_database,
p_cutoff = 0.05,
p_adjust_method = "fdr",
threads = threads
)
("Enriching pathways for each metabolic module...\n")
enriched_pathways_list <-
("list", = (metabolic_module_result))
(enriched_pathways_list) <- metabolic_module_result$Name
for (i 1:(metabolic_module_result$p_value < 0.05)) {
(i, " ")
x <-
metabolic_module_result$Total_metabolite_id[i] %>%
stringr::str_split( = "\\{\\}") %>%
() %>%
()
enriched_pathways_list[i]] <-
(
query_id = x,
query_type = "compound",
id_type = "KEGG",
pathway_database = pathway_database,
p_cutoff = 0.05,
p_adjust_method = "fdr",
threads = threads
)
}
(
(
dysregulated_metabolic_module = metabolic_module_result,
dysregulated_metabolic_network = final_metabolic_network,
annotation_table = annotation_table_final,
enriched_pathways = enriched_pathways,
enriched_pathways_list = enriched_pathways_list
)
)
}
#' Identify Metabolic Modules
#'
#' This function identifies metabolic modules from a metabolic reaction network.
#' It applies a clustering algorithm to detect modules based on network structure.
#' Large modules are recursively split, and weakly connected hidden metabolites are removed.
#'
#' @param sub_metabolic_network A metabolic network represented as a `tbl_graph` object.
#' @param detected_metabolites A character vector of detected metabolite IDs.
#' @param hidden_metabolites A character vector of hidden metabolite IDs.
#'
#' @return A list of metabolic modules, where each module is a character vector
#' containing the metabolite IDs within the module.
#'
#' @export
<-
(sub_metabolic_network,
detected_metabolites,
hidden_metabolites) {
# Apply Walktrap clustering to detect metabolic modules
metabolic_modules <- igraph::cluster_walktrap(graph = sub_metabolic_network)
# Extract module membership information
node_name <- metabolic_modules$
membership <- metabolic_modules$membership
membership_counts <- (membership)
# Group nodes by module membership
<- ((membership_counts), (x) {
node_name[which(membership == (x))]
})
# Remove modules with fewer than 3 nodes
<- [unlist((, )) > 3]
# Recursively split modules containing more than 100 nodes
large_module_idx <- (((, )) > 100)
while ((large_module_idx) > 0) {
large_group <- [large_module_idx]
<- -large_module_idx]
new_group <- (large_group, (x) {
# Extract subnetwork for the large module
temp_graph <- sub_metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% x)
# Apply clustering again within the subnetwork
temp_metabolic_modules <- igraph::cluster_walktrap(graph = temp_graph)
temp_node_name <- temp_metabolic_modules$
temp_membership <- temp_metabolic_modules$membership
temp_membership_counts <- (temp_membership)
# Split into smaller groups
temp_group <- ((temp_membership_counts), (y) {
temp_node_name[which(temp_membership == (y))]
})
# Remove small groups with fewer than 3 nodes
temp_group <- temp_group[unlist((temp_group, )) > 3]
temp_group
})
new_group <- (new_group, recursive = )
<- (, new_group)
# Identify any remaining large modules
large_module_idx <- (((new_group, )) > 100)
}
# Remove weakly connected groups where the number of edges is less than the number of nodes
remove_idx <- (((, (x) {
temp_graph <- sub_metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% x)
node_count <- igraph::gorder(temp_graph) # Number of nodes
edge_count <- igraph::gsize(temp_graph) # Number of edges
degree_values <- igraph::degree(graph = temp_graph, v = x)
# Remove groups where edges are fewer than nodes and all nodes have degree < 3
if (node_count - edge_count > 0 & (degree_values < 3)) {
()
} {
()
}
})))
if ((remove_idx) > 0) {
<- -remove_idx]
}
# Remove hidden metabolites with degree < 2
<- (, (x) {
temp_graph <- sub_metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% x)
degree_values <- igraph::degree(graph = temp_graph, v = x)
metabolite_class <- ((degree_values), (y) {
(y %in% detected_metabolites, "Detected", "Hidden")
})
remove_idx <- (degree_values == 1 &
metabolite_class == "Hidden")
while ((remove_idx) > 0) {
x <- x-remove_idx]
temp_graph <- sub_metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% x)
degree_values <- igraph::degree(graph = temp_graph, v = x)
metabolite_class <- ((degree_values), (y) {
(y %in% detected_metabolites, "Detected", "Hidden")
})
remove_idx <- (degree_values == 1 &
metabolite_class == "Hidden")
}
(x)
})
# Remove groups that still have fewer than 3 nodes
<- [unlist((, )) > 3]
()
}
#' Calculate Node Centrality in a Graph
#'
#' This function computes the centrality of nodes in a graph using either degree centrality or betweenness centrality.
#'
#' @param graph An `igraph` object representing the network.
#' @param type A character string specifying the type of centrality to calculate. Options are "degree" (default) or "betweenness".
#'
#' @return A named vector of centrality values for each node in the graph.
#'
#' @export
#'
<-
(graph, = ("degree", "betweenness")) {
# Ensure that the specified centrality type is valid
<- ()
# Extract node names from the graph
node <- igraph::V(graph)$
# Compute degree centrality if selected
if ( == "degree") {
node_degree <- igraph::degree(graph = graph, v = node)
(node_degree)
}
# Compute betweenness centrality
node_betweenness <- (node, (x) {
i <- (x, node)
if (i == (node))
()
from <- x
to <- node(i + 1):(node)]
# Compute shortest paths between the node and others
temp <- igraph::shortest_paths(graph = graph,
from = from,
to = to)[1]]
# Extract nodes that lie on shortest paths
temp <- (temp, (path) {
if ((path) <= 2)
()
(path)-(1, (path))]
})
(temp)
})
node_betweenness <- (node_betweenness)
node_betweenness <- (node_betweenness)
# Ensure that all nodes are included in the result
missing_nodes <- (node, (node_betweenness))
if ((missing_nodes) > 0) {
add <- (0, (missing_nodes))
(add) <- missing_nodes
node_betweenness <- (node_betweenness, add)
}
# Ensure correct ordering of results
node_betweenness <- node_betweenness[match(node, (node_betweenness))]
(node_betweenness)
}
#' Identify Hidden Metabolites in a Metabolic Network
#'
#' This function identifies hidden metabolites that are connected to detected metabolites
#' within a specified number of reaction steps in a metabolic network.
#'
#' @param metabolic_network An `igraph` object representing the metabolic reaction network.
#' @param detected_metabolites A character vector of detected metabolite IDs.
#' @param threads Integer; number of parallel threads to use. Default is `3`.
#' @param max.reaction Integer; maximum number of reaction steps allowed to define a hidden metabolite. Default is `3`.
#'
#' @return A list of hidden metabolites connected to detected metabolites within the specified reaction steps.
#'
#' @export
#'
<-
(metabolic_network,
detected_metabolites,
threads = 8,
max.reaction = 3) {
# Helper function to find hidden metabolites for each detected metabolite
temp_fun_hidden_metabolites <- (,
metabolic_network,
detected_metabolites,
max.reaction) {
(warn = -1) # Suppress warnings
# Identify remaining detected metabolites for distance calculation
idx <- (, detected_metabolites)
to <- detected_metabolites-(1:idx)]
# Compute shortest path distances from the current metabolite
distance <- igraph::distances(
graph = metabolic_network,
v = ,
to = to,
= "all"
)[1, ]
# Set infinite distances to exceed max.reaction threshold
distance[is.infinite(distance)] <- max.reaction + 10
to <- to[which(distance <= max.reaction)]
if ((to) == 0) {
() # No valid hidden metabolites found
}
# Extract shortest paths between the metabolite and its connected nodes
result <- igraph::shortest_paths(
graph = metabolic_network,
from = ,
to = to,
= "all"
)[1]]
# Extract node names from shortest paths and remove already detected metabolites
result <- (((result, )))
result <- (result, detected_metabolites)
(result)
}
# Select parallel processing backend based on OS
if (masstools::get_os() == "windows") {
bpparam <- BiocParallel::SnowParam(workers = threads, progressbar = )
} {
bpparam <- BiocParallel::MulticoreParam(workers = threads, progressbar = )
}
# Identify hidden metabolites using parallel processing
hidden_metabolites <- BiocParallel::bplapply(
detected_metabolites-(detected_metabolites)],
FUN = temp_fun_hidden_metabolites,
BPPARAM = bpparam,
metabolic_network = metabolic_network,
detected_metabolites = detected_metabolites,
max.reaction = max.reaction
)
(hidden_metabolites)
}
#' Calculate Activity Scores for Metabolic Modules
#'
#' This function calculates activity scores for metabolic modules based on their network properties,
#' including modularity, detected metabolite contributions, and overall network connectivity.
#'
#' @param metabolic_modules A list of metabolic modules, each containing metabolite IDs.
#' @param detected_metabolites A character vector of detected metabolite IDs.
#' @param hidden_metabolites A character vector of hidden metabolite IDs.
#' @param sub_metabolic_network An `igraph` object representing the metabolic subnetwork.
#' @param threads Integer; number of parallel threads to use. Default is `3`.
#'
#' @return A numeric vector of activity scores for each metabolic module.
#'
#' @export
#'
<-
(metabolic_modules,
detected_metabolites,
hidden_metabolites,
sub_metabolic_network,
threads = 8) {
# Helper function to compute activity scores for a given metabolic module
temp_fun_activity_scores <- (x,
detected_metabolites,
hidden_metabolites,
total_metabolite_number,
sub_metabolic_network) {
(warn = -1) # Suppress warnings
# Extract the subgraph for the given metabolic module
temp_graph <- igraph::subgraph(graph = sub_metabolic_network, v = x)
detected_number <- (x %in% detected_metabolites)
hidden_number <- (x %in% hidden_metabolites)
total_edges <- (igraph::E(sub_metabolic_network)) # Total edges in network
module_edges <- (igraph::E(temp_graph)) # Edges within the module
node_names <- (igraph::V(temp_graph))
# Compute modularity-based contribution
value <- (node_names, (node) {
temp_value <- (node_names, (x) {
degree_x <- igraph::degree(graph = sub_metabolic_network, v = x)
degree_node <- igraph::degree(graph = sub_metabolic_network, v = node)
(degree_x * degree_node) / (4 * total_edges^2)
})
(temp_value)
})
modularity <- module_edges / total_edges - (value)
# Compute activity score
<- modularity * (total_metabolite_number / (x))
activity <- detected_number * / (x)
(activity)
}
# Choose appropriate parallelization method based on OS
if (masstools::get_os() == "windows") {
bpparam <- BiocParallel::SnowParam(workers = threads, progressbar = )
} {
bpparam <- BiocParallel::MulticoreParam(workers = threads, progressbar = )
}
# Compute activity scores for all metabolic modules in parallel
activity_scores <- BiocParallel::bplapply(
metabolic_modules,
temp_fun_activity_scores,
BPPARAM = bpparam,
detected_metabolites = detected_metabolites,
hidden_metabolites = hidden_metabolites,
total_metabolite_number = (detected_metabolites),
sub_metabolic_network = sub_metabolic_network
)
((activity_scores))
}
#' Annotate Metabolites for Feature Pathway Analysis (FPA)
#'
#' This function performs metabolite annotation for a given feature table using an MS1-based approach.
#' It processes both positive and negative ionization modes separately and extracts annotation results.
#'
#' @param feature_table A data frame containing feature data, including `variable_id`, `mz`, `rt`, and `polarity` columns.
#' @param metabolite_database A metabolite database object used for annotation.
#' @param column A character string specifying the chromatography column type. Options are `"rp"` (reverse phase) or `"hilic"` (hydrophilic interaction chromatography). Default is `"rp"`.
#' @param adduct.table Optional data frame specifying possible adducts. Default is `NULL`.
#' @param ms1.match.ppm Numeric; mass tolerance for MS1 matching in parts per million (ppm). Default is `10`.
#' @param rt.match.tol Numeric; retention time tolerance threshold for matching in seconds. Default is `5`.
#' @param mz.ppm.thr Numeric; M/Z tolerance threshold for filtering in ppm. Default is `400`.
#' @param threads Integer; number of parallel threads to use for annotation. Default is `3`.
#'
#' @return A data frame containing the annotated metabolites with polarity information.
#'
#' @export
#'
<-
(feature_table,
metabolite_database,
column = ("rp", "hilic"),
adduct.table = ,
ms1.match.ppm = 10,
rt.match.tol = 5,
mz.ppm.thr = 400,
threads = 3) {
column <- (column) # Ensure valid column type
##### Separate features by polarity #####
feature_table_pos <- feature_table %>%
dplyr::(polarity == "positive")
feature_table_neg <-
feature_table %>% dplyr::(polarity == "negative")
##### Create mass dataset objects for positive and negative modes #####
# Positive mode processing
if ((feature_table_pos) == 0) {
object_pos <-
} {
expression_data_pos <- (sample1 = 1:(feature_table_pos))
(expression_data_pos) <- feature_table_pos$variable_id
sample_info_pos <- (sample_id = "sample1", = "Subject")
object_pos <-
massdataset::create_mass_dataset(
expression_data = expression_data_pos,
sample_info = sample_info_pos,
variable_info = feature_table_pos
)
}
# Negative mode processing
if ((feature_table_neg) == 0) {
object_neg <-
} {
expression_data_neg <- (sample1 = 1:(feature_table_neg))
(expression_data_neg) <- feature_table_neg$variable_id
sample_info_neg <- (sample_id = "sample1", = "Subject")
object_neg <- massdataset::create_mass_dataset(
expression_data = expression_data_neg,
sample_info = sample_info_neg,
variable_info = feature_table_neg
)
}
##### Perform metabolite annotation #####
("Annotating metabolites in positive mode...\n")
if (!(object_pos)) {
object_pos <- metid::annotate_metabolites(
object = object_pos,
database = metabolite_database,
based_on = "ms1",
polarity = "positive",
column = column,
adduct.table = adduct.table,
ms1.match.ppm = ms1.match.ppm,
candidate.num = 100,
threads = threads,
mz.ppm.thr = mz.ppm.thr
)
annotation_table_pos <-
massdataset::extract_annotation_table(object_pos) %>%
dplyr::mutate(polarity = "positive")
####add feature mz, rt and annotation formula to it
annotation_table_pos <-
annotation_table_pos %>%
dplyr::left_join(feature_table_pos, ("variable_id", "mz", "rt")], = ("variable_id" = "variable_id")) %>%
dplyr::left_join(metabolite_database@spectra.info, ("Lab.ID", "Formula")],
= ("Lab.ID" = "Lab.ID")) %>%
dplyr::select(variable_id, mz, , Formula, Adduct, dplyr::everything()) %>%
dplyr::mutate(isotope = "[M]")
####isotope annotation for annotations
("Annotating isotopes in positive mode...\n")
future::plan(future::multisession, workers = threads)
(
isotope_pos <-
furrr::future_map(
.x = ((annotation_table_pos)),
.f = (x) {
adduct <-
stringr::str_extract(x[5], "\\(.+\\)") %>%
stringr::str_replace("\\(", "") %>%
stringr::str_replace("\\)", "")
mz = (stringr::str_trim(x[2], side = "both"))
= (stringr::str_trim(x[3], side = "both"))
temp_iso <- ((
= stringr::str_trim(x[4], side = "both"),
adduct = adduct,
mz = mz,
= ,
peak.mz = feature_table_pos$mz,
peak.rt = feature_table_pos$,
rt.tol = rt.match.tol,
mz.tol = ms1.match.ppm,
max.isotope = 3
),
silent = )
if ((temp_iso, "try-error")) {
()
}
if ((temp_iso)) {
()
}
# return(temp_iso)
temp_iso <-
(feature_table_pos[temp_iso$peakIndex, ], temp_iso) %>%
dplyr::select(-(peakIndex))
(temp_iso) <- ("variable_id",
"mz",
"rt",
"polarity",
"mz.error",
"isotope",
"RT.error")
x <- (x, = 1) %>% ()
(x) <- (annotation_table_pos)
temp_iso$Compound.name <- x$Compound.name
temp_iso$Lab.ID <- x$Lab.ID
temp_iso$Adduct <- x$Adduct
temp_iso$Formula <- x$Formula
# temp_iso$rt <- x$rt
temp_iso$CAS.ID <- x$CAS.ID
temp_iso$HMDB.ID <- x$HMDB.ID
temp_iso$KEGG.ID <- x$KEGG.ID
temp_iso$Database <- x$Database
temp_iso$ms2_files_id <-
temp_iso$ms2_spectrum_id <-
temp_iso$CE <-
temp_iso$SS <-
temp_iso$Total.score <-
temp_iso$Level <-
temp_iso$mz.match.score <-
(ms1.match.ppm - temp_iso$mz.error) / ms1.match.ppm
temp_iso$RT.match.score <-
(rt.match.tol - temp_iso$RT.error) / rt.match.tol
temp_iso %>%
dplyr::select((x))
},
.progress =
)
)
isotope_pos <- isotope_pos %>%
(, .) %>%
()
annotation_table_pos <-
(annotation_table_pos, isotope_pos) %>%
dplyr::(Compound.name, isotope, )
} {
annotation_table_pos <-
}
("Annotating metabolites in negative mode...\n")
if (!(object_neg)) {
object_neg <- metid::annotate_metabolites(
object = object_neg,
database = metabolite_database,
based_on = "ms1",
polarity = "negative",
column = column,
adduct.table = adduct.table,
ms1.match.ppm = ms1.match.ppm,
candidate.num = 100,
threads = threads,
mz.ppm.thr = mz.ppm.thr
)
annotation_table_neg <-
massdataset::extract_annotation_table(object_neg) %>%
dplyr::mutate(polarity = "negative")
####add feature mz, rt and annotation formula to it
annotation_table_neg <-
annotation_table_neg %>%
dplyr::left_join(feature_table_neg, ("variable_id", "mz", "rt")], = ("variable_id" = "variable_id")) %>%
dplyr::left_join(metabolite_database@spectra.info, ("Lab.ID", "Formula")],
= ("Lab.ID" = "Lab.ID")) %>%
dplyr::select(variable_id, mz, , Formula, Adduct, dplyr::everything()) %>%
dplyr::mutate(isotope = "[M]")
####isotope annotation for annotations
("Annotating isotopes in negative mode...\n")
future::plan(future::multisession, workers = threads)
(
isotope_neg <-
furrr::future_map(
.x = ((annotation_table_neg)),
.f = (x) {
adduct <-
stringr::str_extract(x[5], "\\(.+\\)") %>%
stringr::str_replace("\\(", "") %>%
stringr::str_replace("\\)", "")
mz = (stringr::str_trim(x[2], side = "both"))
= (stringr::str_trim(x[3], side = "both"))
temp_iso <- ((
= stringr::str_trim(x[4], side = "both"),
adduct = adduct,
mz = mz,
= ,
peak.mz = feature_table_neg$mz,
peak.rt = feature_table_neg$,
rt.tol = rt.match.tol,
mz.tol = ms1.match.ppm,
max.isotope = 3
),
silent = )
if ((temp_iso, "try-error")) {
()
}
if ((temp_iso)) {
()
}
# return(temp_iso)
temp_iso <-
(feature_table_neg[temp_iso$peakIndex, ], temp_iso) %>%
dplyr::select(-(peakIndex))
(temp_iso) <- ("variable_id",
"mz",
"rt",
"polarity",
"mz.error",
"isotope",
"RT.error")
x <- (x, = 1) %>% ()
(x) <- (annotation_table_neg)
temp_iso$Compound.name <- x$Compound.name
temp_iso$Lab.ID <- x$Lab.ID
temp_iso$Adduct <- x$Adduct
temp_iso$Formula <- x$Formula
# temp_iso$rt <- x$rt
temp_iso$CAS.ID <- x$CAS.ID
temp_iso$HMDB.ID <- x$HMDB.ID
temp_iso$KEGG.ID <- x$KEGG.ID
temp_iso$Database <- x$Database
temp_iso$ms2_files_id <-
temp_iso$ms2_spectrum_id <-
temp_iso$CE <-
temp_iso$SS <-
temp_iso$Total.score <-
temp_iso$Level <-
temp_iso$mz.match.score <-
(ms1.match.ppm - temp_iso$mz.error) / ms1.match.ppm
temp_iso$RT.match.score <-
(rt.match.tol - temp_iso$RT.error) / rt.match.tol
temp_iso %>%
dplyr::select((x))
},
.progress =
)
)
isotope_neg <- isotope_neg %>%
(, .) %>%
()
annotation_table_neg <-
(annotation_table_neg, isotope_neg) %>%
dplyr::(Compound.name, isotope, )
} {
annotation_table_neg <-
}
##### Combine annotation tables #####
annotation_table <-
(annotation_table_pos, annotation_table_neg)
(annotation_table)
}
#' Generate Null Distribution of Activity Scores
#'
#' This function generates a null distribution of activity scores by performing random sampling
#' of marker metabolites and calculating activity scores through network analysis.
#'
#' @param annotation_table_all A data frame containing all annotated metabolites.
#' @param feature_table_marker A data frame containing marker features.
#' @param metabolite_database A metabolite database object used for annotation.
#' @param metabolic_network An `igraph` object representing the metabolic network.
#' @param permutation_times Integer; number of permutations to perform. Default is `20`.
#' @param threads Integer; number of parallel threads to use. Default is `8`.
#' @param adduct.table Optional data frame specifying possible adducts. Default is `NULL`.
#' @param column Character; chromatography column type. Options are `"rp"` or `"hilic"`. Default is `"rp"`.
#' @param ms1.match.ppm Numeric; MS1 matching tolerance in ppm. Default is `25`.
#' @param mz.ppm.thr Numeric; M/Z tolerance threshold in ppm. Default is `400`.
#' @param include_hidden_metabolites Logical; whether to include hidden metabolites in the network analysis. Default is `FALSE`.
#'
#' @return A list of numeric vectors representing null activity scores for each permutation.
#'
#' @export
#'
<-
(annotation_table_all,
feature_table_marker,
metabolite_database,
metabolic_network,
permutation_times = 20,
threads = 8,
adduct.table = ,
column = ("rp", "hilic"),
ms1.match.ppm = 25,
mz.ppm.thr = 400,
include_hidden_metabolites = ) {
column <- (column) # Ensure valid column type
# Initialize list to store activity scores for each permutation
null_activity_score <- ( = "list", = permutation_times)
marker_number <- (feature_table_marker) # Number of markers
for (i (permutation_times)) {
(i, " ")
# Randomly sample markers for permutation
if ((feature_table_marker) / ((annotation_table_all$variable_id)) < 0.2) {
temp_marker_variable_id <- (
x = (
(annotation_table_all$variable_id),
feature_table_marker$variable_id
),
size = marker_number,
=
)
} {
temp_marker_variable_id <- (
x = (annotation_table_all$variable_id),
size = marker_number,
=
)
}
# Extract annotation information for sampled markers
temp_annotation_table_marker <- annotation_table_all %>%
dplyr::(variable_id %in% temp_marker_variable_id)
# Extract detected metabolites from annotation table
temp_detected_metabolites <- (temp_annotation_table_marker$KEGG.ID)
# Identify hidden metabolites connected to detected metabolites
if (include_hidden_metabolites) {
temp_hidden_metabolites <- (
metabolic_network = metabolic_network,
detected_metabolites = temp_detected_metabolites,
threads = threads,
max.reaction = 3
) %>% () %>% ()
} {
temp_hidden_metabolites <-
}
temp_total_metabolites <- (temp_detected_metabolites, temp_hidden_metabolites)
# Extract subnetwork containing detected and hidden metabolites
temp_sub_metabolic_network <- metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% temp_total_metabolites) %>%
dplyr::mutate(node_class = ( %in% temp_detected_metabolites, "Detected", "Hidden")) %>%
tidygraph::mutate(degree2 = tidygraph::centrality_degree())
# Identify metabolic modules in the subnetwork
temp_metabolic_modules <- (
sub_metabolic_network = temp_sub_metabolic_network,
detected_metabolites = temp_detected_metabolites,
hidden_metabolites = temp_hidden_metabolites
)
# Calculate activity scores for the identified metabolic modules
temp_activity_score <- (
metabolic_modules = temp_metabolic_modules,
detected_metabolites = temp_detected_metabolites,
hidden_metabolites = temp_hidden_metabolites,
sub_metabolic_network = temp_sub_metabolic_network
)
# Store the computed activity score in the list
null_activity_score[i]] <- temp_activity_score
}
(null_activity_score)
}
#####g
#' Group Retention Time (RT) Peaks
#'
#' This function groups retention time (RT) values based on a specified tolerance,
#' clustering peaks that fall within the same RT range.
#'
#' @param rt A numeric vector of retention times.
#' @param rt.tol Numeric; retention time tolerance for grouping. Default is `10`.
#'
#' @return A data frame with two columns:
#' - `rt`: The retention times of grouped peaks.
#' - `class`: The assigned group ID.
#'
#' @export
#'
<-
(, rt.tol = 10) {
# Identify clusters of retention times within the specified tolerance
rt_class <- (, (x) {
( >= x & < x + rt.tol)
})
# Ensure each peak is uniquely assigned to a single class
rt_class <- ((rt_class)-1], (i) {
(rt_class[i]], ((rt_class[1:(i - 1)])))
}) %>% (rt_class[1], .)
# Remove empty groups
rt_class <- rt_class[which((rt_class, ) != 0)]
# Assign group names
(rt_class) <- (rt_class)
# Convert grouping information into a data frame
rt_class <- purrr::map2(
.x = rt_class,
.y = (rt_class),
.f = (x, y) {
( = [x],
= y,
stringsAsFactors = )
}
) %>% (, .)
(rt_class) <-
(rt_class)
}
#' Score a Peak Group Based on Adduct and Polarity Information
#'
#' This function assigns a score to a peak group based on the detected adduct types,
#' polarity, and isotopic information. Higher scores indicate a higher likelihood
#' of correct identification.
#'
#' @param peak_group A data frame containing peak group information, including `Adduct`, `polarity`, and `isotope` columns.
#'
#' @return A numeric score representing the confidence of the peak group assignment.
#'
#' @details
#' The scoring system is based on the following criteria:
#' - `+50` if the positive adduct is `(M+H)+`.
#' - `+20` if `(M+H)+` is present and is not the `[M]` isotope.
#' - `+20` if the polarity is positive but the adduct is not `(M+H)+`.
#' - `+10` if the polarity is positive, the adduct is not `(M+H)+`, and the isotope is not `[M]`.
#' - `+50` if the negative adduct is `(M-H)-`.
#' - `+20` if `(M-H)-` is present and is not the `[M]` isotope.
#' - `+20` if the polarity is negative but the adduct is not `(M-H)-`.
#' - `+10` if the polarity is negative, the adduct is not `(M-H)-`, and the isotope is not `[M]`.
#'
#' @export
#'
<-
(peak_group) {
score <- 0
## This should be optimized in the future
# Positive mode scoring
if ((peak_group$Adduct == "(M+H)+")) {
score <- score + 50
}
if ((peak_group$Adduct == "(M+H)+" &
peak_group$isotope != "[M]")) {
score <- score + 20
}
if ((peak_group$polarity == "positive" &
peak_group$Adduct != "(M+H)+")) {
score <- score + 20
}
if ((
peak_group$polarity == "positive" &
peak_group$Adduct != "(M+H)+" &
peak_group$isotope != "[M]"
)) {
score <- score + 10
}
# Negative mode scoring
if ((peak_group$Adduct == "(M-H)-")) {
score <- score + 50
}
if ((peak_group$Adduct == "(M-H)-" &
peak_group$isotope != "[M]")) {
score <- score + 20
}
if ((peak_group$polarity == "negative" &
peak_group$Adduct != "(M-H)-")) {
score <- score + 20
}
if ((
peak_group$polarity == "negative" &
peak_group$Adduct != "(M-H)-" &
peak_group$isotope != "[M]"
)) {
score <- score + 10
}
(score)
}
#' Calculate Redundancy in Metabolite Annotation
#'
#' This function calculates redundancy metrics in a metabolite annotation table.
#' It evaluates the number of distinct compound classes per compound and the number
#' of compounds assigned to each peak.
#'
#' @param annotation_table A data frame containing metabolite annotations, including `Lab.ID`, `compound_class`, and `variable_id`.
#'
#' @return A numeric vector of length two:
#' - The first value represents the average number of compound classes per compound (`redundancy1`).
#' - The second value represents the average number of compounds assigned to each peak (`redundancy2`).
#'
#' @export
#'
<-
(annotation_table) {
# Ensure the annotation table is in a consistent format
annotation_table <- dplyr::bind_rows(annotation_table)
# Redundancy 1: Average number of compound classes per compound
redundancy1 <- annotation_table %>%
dplyr::group_by(Lab.ID) %>%
dplyr::distinct(compound_class) %>%
dplyr::summarise(n = dplyr::n()) %>%
dplyr::pull(n) %>%
()
# Redundancy 2: Average number of compounds assigned to each peak
redundancy2 <- annotation_table %>%
dplyr::group_by(variable_id) %>%
dplyr::summarise(n = dplyr::n()) %>%
dplyr::pull(n) %>%
()
((redundancy1, redundancy2))
}
#' Remove Redundant Annotations in Metabolite Data
#'
#' This function removes redundant metabolite annotations by filtering out low-confidence matches.
#' It iteratively refines the annotation table based on score thresholds.
#'
#' @param annotation_table A data frame containing metabolite annotations,
#' including `Lab.ID`, `variable_id`, and `score` columns.
#'
#' @return A filtered version of the annotation table with reduced redundancy.
#'
#' @details
#' The function follows these rules:
#' - If a compound has an annotation with a score > 80, it removes annotations with scores ≤ 20.
#' - If a peak has an annotation with a score > 80, it removes other annotations with scores ≤ 20.
#' - Iterates until redundancy stabilizes, ensuring that lower-quality annotations are removed.
#'
#' @export
#'
<-
(annotation_table) {
## Track redundancy change across iterations
redundancy_diff <- (-1, -1)
while ((redundancy_diff < 0)) {
# Calculate redundancy before filtering
before_redundancy <- (annotation_table = annotation_table)
# Filter compound annotations: Remove low-scoring annotations if high-scoring ones exist
annotation_table <- annotation_table %>%
dplyr::group_by(Lab.ID) %>%
dplyr::(if ((score > 80)) {
score > 20
} {
score > 0
}) %>%
dplyr::ungroup()
# Filter peak annotations: Remove lower-scoring annotations when a high-scoring annotation exists
annotation_table <- annotation_table %>%
dplyr::group_by(variable_id) %>%
dplyr::(if ((score > 80)) {
score > 20
} {
score > 0
}) %>%
dplyr::ungroup()
# Calculate redundancy after filtering
after_redundancy <- (annotation_table = annotation_table)
# Compute change in redundancy
redundancy_diff <- after_redundancy - before_redundancy
}
(annotation_table)
}
#' Annotate Isotopes for a Given Molecular Formula
#'
#' This function identifies isotopic peaks for a given molecular formula and adduct,
#' matching them with observed peaks based on mass-to-charge ratio (m/z) and retention time (RT).
#'
#' @param formula Character. The molecular formula of the compound (e.g., `"C9H14N2O12P2"`).
#' @param adduct Character. The adduct ion type (e.g., `"M-H"`).
#' @param charge Numeric. The charge state of the ion. Default is `1`.
#' @param mz Numeric. The accurate mass-to-charge ratio (m/z) of the monoisotopic peak.
#' @param rt Numeric. The retention time (RT) of the compound in seconds.
#' @param peak.mz Numeric vector. Observed m/z values of peaks in the dataset.
#' @param peak.rt Numeric vector. Observed retention times (RT) of peaks in the dataset.
#' @param rt.tol Numeric. The retention time tolerance (in seconds) for matching peaks. Default is `5`.
#' @param mz.tol Numeric. The m/z tolerance (in ppm) for matching peaks. Default is `15`.
#' @param max.isotope Numeric. The maximum number of isotopes to consider. Default is `4`.
#'
#' @return A data frame containing:
#' \item{peakIndex}{Index of the matched peak in the dataset.}
#' \item{mzError.ppm}{Mass error (in ppm) between the observed and theoretical isotope peak.}
#' \item{isotopes}{Label of the identified isotope (e.g., `[M+1]`, `[M+2]`).}
#' \item{rtError.s}{Retention time error (in seconds).}
#'
#' If no isotopic peaks are matched within the given tolerances, the function returns `NULL`.
#'
#' @details
#' The function calculates the theoretical isotope distribution using `Rdisop::getMolecule()`
#' and then matches observed peaks based on m/z and RT tolerances. The function ensures that
#' at least the `[M+1]` isotope is present for a valid match.
#'
#' @importFrom masstools sum_formula
#' @importFrom Rdisop getMolecule getIsotope
#' @importFrom purrr map
#' @importFrom dplyr mutate
#' @export
<-
( = "C9H14N2O12P2",
adduct = "M-H",
charge = 1,
mz = 402.9998,
= 823.1462,
## peak information
peak.mz,
peak.rt,
## other parameters
rt.tol = 5,
mz.tol = 15,
max.isotope = 4) {
formula1 <- masstools::sum_formula( = , adduct = adduct)
###should be fix latter
if ((formula1)) {
formula1 <-
}
molecule <- Rdisop::getMolecule( = formula1, # z = charge,
maxisotopes = max.isotope + 1)
isotopes <- (Rdisop::getIsotope(molecule = molecule))
(isotopes) <-
("[M]", ("[M", "+", (1:((
isotopes
) - 1)), "]", sep = ""))
isotopes <- (isotopes, (isotopes), stringsAsFactors = )
(isotopes) <- ("mz", "intensity", "isotope")
accurate.mz <- mz
isotopes <- isotopes-1, , = ]
###rt filtering
rt.error <- ( - peak.rt)
index1 <- (rt.error <= rt.tol)
if ((index1) == 0)
()
iso.info <-
purrr::map(
((isotopes)),
.f = (x) {
temp.mz <- (x[1])
## calculate error
peak.mz.error <-
(temp.mz - peak.mz) * 10^6 / (temp.mz >= 400, temp.mz, 400)
##has peak matched the mz tolerance
idx <- (peak.mz.error <= mz.tol & rt.error < rt.tol)
###idx=0, no peaks matched
if ((idx) == 0) {
# return(c(NA, NA, NA))
()
}
## more than one matched, see the intensity ratio error
if ((idx) > 0) {
idx <- idx[which.min(peak.mz.error[idx])]
((idx, peak.mz.error[idx], x[3]))
}
}
) %>%
(, .) %>%
()
if ((iso.info) == 0) {
()
}
if (!"[M+1]" %in% iso.info$V3) {
()
}
(iso.info) <- ("peakIndex", "mzError.ppm", "isotopes")
if (!"[M+2]" %in% iso.info$isotopes) {
iso.info <- iso.info[1, , = ]
}
iso.info$`rtError.s` <-
rt.error[as.numeric(iso.info$peakIndex)]
iso.info <-
iso.info %>%
dplyr::mutate(
peakIndex = (peakIndex),
mzError.ppm = (mzError.ppm),
rtError.s = (rtError.s)
)
if (!"[M+2]" %in% iso.info$isotopes) {
iso.info <- iso.info[1, ]
}
# rm(list = c("peak.mz", "peak.rt", "peak.int", "cor"))
# gc()
(iso.info)
}
#' Plot a Metabolic Module from FPA Results
#'
#' This function visualizes a selected metabolic module from the results of a
#' Functional Pathway Analysis (FPA), displaying interactions between metabolites
#' and features in the metabolic network.
#'
#' @param fpa_result List. The output of the Functional Pathway Analysis (FPA),
#' containing `dysregulated_metabolic_module`, `dysregulated_metabolic_network`,
#' and `annotation_table`.
#' @param feature_table_marker Data frame. The feature table containing metabolic
#' markers. Default is `feature_table_marker`.
#' @param include_feature Logical. Whether to include detected metabolic features
#' in the plot. Default is `FALSE`.
#' @param include_hidden_metabolites Logical. Whether to include hidden metabolites
#' in the plot. Default is `FALSE`.
#' @param add_compound_name Logical. Whether to add compound names as labels in
#' the visualization. Default is `TRUE`.
#' @param metabolic_module_index Numeric. The index of the metabolic module to
#' visualize. Default is `1`.
#' @param layout The layout of the network, such as `kk` or `fr`.
#' @param add_pathways Add pathways beside of the network or not. Default is `FALSE`.
#'
#' @return A `ggplot2` object representing the metabolic module network.
#'
#' @details
#' The function extracts and filters a metabolic subnetwork based on the selected
#' module index. It then uses `ggraph` to visualize the metabolic interactions
#' while allowing customization, such as including metabolic features or hidden
#' metabolites.
#'
#' The network layout follows the Fruchterman-Reingold force-directed algorithm
#' (`layout = "fr"`). Nodes represent metabolites and features, while edges
#' represent different types of metabolic interactions.
#'
#' @importFrom stringr str_split
#' @importFrom dplyr filter pull mutate
#' @importFrom tidygraph activate centrality_degree
#' @importFrom ggraph ggraph geom_edge_arc geom_node_point geom_node_text theme_graph
#' @importFrom ggraph scale_edge_color_manual
#' @export
<-
(fpa_result,
feature_table_marker = feature_table_marker,
include_feature = ,
include_hidden_metabolites = ,
add_compound_name = ,
metabolic_module_index = 1,
= "fr",
add_pathways = ) {
dysregulated_metabolic_module <-
fpa_result$dysregulated_metabolic_module
dysregulated_metabolic_network <-
fpa_result$dysregulated_metabolic_network
annotation_table <-
fpa_result$annotation_table
remain_compound <-
dysregulated_metabolic_module$Total_metabolite_id[metabolic_module_index] %>%
stringr::str_split( = "\\{\\}") %>%
() %>%
()
remain_feature <-
annotation_table %>%
dplyr::(KEGG.ID %in% remain_compound) %>%
dplyr::pull(variable_id) %>%
()
remain_feature <-
remain_feature[remain_feature %in% igraph::V(dysregulated_metabolic_network)$]
temp_graph <-
dysregulated_metabolic_network %>%
dplyr::( %in% (remain_compound, remain_feature)) %>%
tidygraph::activate("nodes") %>%
dplyr::mutate(degree = tidygraph::centrality_degree()) %>%
tidygraph::activate("nodes") %>%
dplyr::(degree > 0)
if (!include_feature) {
temp_graph <-
temp_graph %>%
tidygraph::activate("nodes") %>%
dplyr::(node_class != "Feature")
}
if (!include_hidden_metabolites) {
temp_graph <-
temp_graph %>%
tidygraph::activate("nodes") %>%
dplyr::(node_class != "Hidden")
}
<-
temp_graph %>%
ggraph::ggraph( = ) +
ggraph::geom_edge_arc(strength = 0.1, aes(color = edge_class)) +
ggraph::scale_edge_color_manual(values = (
"feature_metabolite" = "#eaeaea",
"metabolic_reaction" = "#6e8da9"
)) +
ggraph::geom_node_point(aes(
color = node_class,
shape = node_class,
size = degree
)) +
scale_color_manual(values = (
"Feature" = "#ffca3a",
"Detected" = "#ff595e",
"Hidden" = "#6e8da9"
)) +
ggraph::theme_graph()
if (add_compound_name) {
<- +
ggraph::geom_node_text(aes(label = Compound_name), repel = )
}
Total_metabolite_number <-
fpa_result$dysregulated_metabolic_module$Total_metabolite_number[metabolic_module_index]
Impact <-
fpa_result$dysregulated_metabolic_module$Impact[metabolic_module_index]
Dominant_edge_rate <-
fpa_result$dysregulated_metabolic_module$Dominant_edge_rate[metabolic_module_index]
p_value <-
fpa_result$dysregulated_metabolic_module$p_value[metabolic_module_index]
Detected_metabolite_number <-
fpa_result$dysregulated_metabolic_module$Detected_metabolite_number[metabolic_module_index]
<-
+
labs(
= fpa_result$dysregulated_metabolic_module$Name[metabolic_module_index],
subtitle = (
"Total Metabolite Number: ",
Total_metabolite_number,
"\n",
"Detected Metabolite Number: ",
Detected_metabolite_number,
"\n",
"Impact: ",
Impact,
"\n",
"Dominant Edge Rate: ",
(Dominant_edge_rate, 3),
"\n",
"P-value: ",
p_value
)
)
if (add_pathways) {
pathway <-
fpa_result$enriched_pathways_list[metabolic_module_index]]
if ((pathway)) {
plot_pathway <-
ggplot() +
theme_bw() +
labs( = "No enriched pathways")
} {
plot_pathway <-
(object = pathway, label = )
}
<-
+ patchwork::wrap_plots(plot_pathway, = 1)
}
}
#' Plot the Dysregulated Metabolic Network from FPA Results
#'
#' This function visualizes the dysregulated metabolic network obtained from Functional Pathway Analysis (FPA),
#' highlighting interactions between metabolites and features in the network.
#'
#' @param fpa_result List. The output of the Functional Pathway Analysis (FPA),
#' containing `dysregulated_metabolic_module`, `dysregulated_metabolic_network`,
#' and `annotation_table`.
#' @param feature_table_marker Data frame. The feature table containing metabolic markers.
#' @param include_feature Logical. Whether to include detected metabolic features in the plot.
#' Default is `FALSE`.
#' @param node_color_by_module Logical. When doesn't include features (include_feature = FALSE),
#' whether set the colors of the nodes by metabolics modules.
#' @param include_hidden_metabolites Logical. Whether to include hidden metabolites in the plot.
#' Default is `FALSE`.
#' @param add_compound_name Logical. Whether to add compound names as labels in
#' the visualization. Default is `TRUE`.
#' @param layout The layout of the network, such as `kk` or `fr`.
#'
#' @return A `ggplot2` object representing the metabolic network.
#'
#' @details
#' The function processes the metabolic network by filtering and selecting relevant nodes and edges.
#' It allows the inclusion of metabolic features and hidden metabolites based on user preference.
#' The network is visualized using the Fruchterman-Reingold force-directed algorithm (`layout = "fr"`)
#' to clearly illustrate metabolic interactions.
#'
#' @importFrom tidygraph activate centrality_degree
#' @importFrom dplyr mutate filter
#' @importFrom ggraph ggraph geom_edge_arc geom_node_point geom_node_text theme_graph
#' @importFrom ggplot2 scale_color_manual
#' @importFrom grDevices colorRampPalette
#' @importFrom stats quantile
#' @export
#'
#' @examples
#' \dontrun{
#' # Assume `fpa_result` is a precomputed list from an FPA analysis
#' plot_metabolic_network_fpa(
#' fpa_result = fpa_result,
#' feature_table_marker = feature_table_marker,
#' include_feature = TRUE,
#' include_hidden_metabolites = FALSE,
#' add_compound_name = TRUE
#' )
#' }
<-
(fpa_result,
feature_table_marker,
include_feature = ,
node_color_by_module = ,
include_hidden_metabolites = ,
add_compound_name = ,
= "fr") {
dysregulated_metabolic_module <-
fpa_result$dysregulated_metabolic_module
dysregulated_metabolic_network <-
fpa_result$dysregulated_metabolic_network
annotation_table <-
fpa_result$annotation_table
temp_graph <-
dysregulated_metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::mutate(degree = tidygraph::centrality_degree())
if (!include_feature) {
temp_graph <-
temp_graph %>%
tidygraph::activate("nodes") %>%
dplyr::(node_class != "Feature")
if (node_color_by_module) {
temp_data <-
((dysregulated_metabolic_module)) %>%
purrr::map((i) {
(
module = dysregulated_metabolic_module$Name[i],
KEGG.ID = stringr::str_split(
dysregulated_metabolic_module$Total_metabolite_id[i],
"\\{\\}"
)[1]]
) %>%
dplyr::distinct(KEGG.ID, .keep_all = )
}) %>%
(, .) %>%
() %>%
dplyr::mutate(module = (module, = stringr::str_sort((module), = )))
temp_graph <-
temp_graph %>%
tidygraph::activate("nodes") %>%
dplyr::left_join(temp_data, = ("name" = "KEGG.ID"))
}
}
if (!include_hidden_metabolites) {
temp_graph <-
temp_graph %>%
tidygraph::activate("nodes") %>%
dplyr::(node_class != "Hidden")
}
if (!include_feature & node_color_by_module) {
<-
temp_graph %>%
ggraph::ggraph( = ) +
ggraph::geom_edge_arc(strength = 0.1, aes(color = edge_class)) +
ggraph::scale_edge_color_manual(values = (
"feature_metabolite" = "#eaeaea",
"metabolic_reaction" = "#6e8da9"
)) +
ggraph::geom_node_point(aes(
color = module,
shape = node_class,
size = degree
)) +
scale_color_manual(values =
(ggsci::pal_aaas("default")(10))(((
(
tidygraph::as_tibble(temp_graph, active = "nodes") %>% pull(module)
)
)))) +
ggraph::theme_graph()
} {
<-
temp_graph %>%
ggraph::ggraph( = ) +
ggraph::geom_edge_arc(strength = 0.1, aes(color = edge_class)) +
ggraph::scale_edge_color_manual(values = (
"feature_metabolite" = "#eaeaea",
"metabolic_reaction" = "#6e8da9"
)) +
ggraph::geom_node_point(aes(
color = node_class,
shape = node_class,
size = degree
)) +
scale_color_manual(values = (
"Feature" = "#ffca3a",
"Detected" = "#ff595e",
"Hidden" = "#6e8da9"
)) +
ggraph::theme_graph()
}
if (add_compound_name) {
<- +
ggraph::geom_node_text(aes(label = Compound_name), repel = )
}
}
#' Calculate Node Centrality in a Graph
#'
#' This function calculates the centrality of nodes in a given network graph based on either **degree centrality** or **betweenness centrality**.
#'
#' @param graph An `igraph` object representing the network.
#' @param type A character string specifying the type of centrality to compute. Options are `"degree"` (default) or `"betweenness"`.
#'
#' @return A named numeric vector where names represent node IDs and values represent centrality scores.
#'
#' @details
#' - **Degree centrality** measures the number of direct connections a node has.
#' - **Betweenness centrality** measures the number of shortest paths passing through a node, indicating its role as a bridge in the network.
#' - The function uses `igraph::degree()` for degree centrality and a custom approach for betweenness centrality using `igraph::shortest_paths()`.
#'
#' @importFrom igraph V degree shortest_paths
#'
#' @examples
#' \dontrun{
#' library(igraph)
#'
#' # Create a sample graph
#' g <- make_ring(10)
#' V(g)$name <- paste0("Node", 1:10)
#'
#' # Calculate degree centrality
#' degree_centrality <- calculate_centrality(g, type = "degree")
#' print(degree_centrality)
#'
#' # Calculate betweenness centrality
#' betweenness_centrality <- calculate_centrality(g, type = "betweenness")
#' print(betweenness_centrality)
#' }
#'
#' @export
<-
(graph, = ("degree", "betweenness")) {
= ()
node <- igraph::V(graph)$
if ( == "degree") {
node.degree <- igraph::degree(graph = graph, v = node)
(node.degree)
}
node.betweenness <- (node, (x) {
i <- (x, node)
if (i == (node))
()
from <- x
to <- node(i + 1):(node)]
temp <- igraph::shortest_paths(graph = graph,
from = from,
to = to)[1]]
temp <- (temp, (x) {
if ((x) <= 2)
()
(x)-(1, (x))]
})
(temp)
})
node.betweenness <- (node.betweenness)
node.betweenness <- (node.betweenness)
node0 <- (node, (node.betweenness))
if ((node0) == 0)
(node.betweenness)
add <- (0, (node0))
(add) <- node0
node.betweenness <- (node.betweenness, add)
node.betweenness <- node.betweenness[match(node, (node.betweenness))]
(node.betweenness)
}
#' Calculate Activity Score for Metabolic Modules
#'
#' This function calculates an activity score for each metabolic module based on its detected and hidden metabolites,
#' as well as its structure within a sub-metabolic network. The calculation incorporates modularity and degree-based
#' weighting to assess the significance of the module.
#'
#' @param metabolic_modules A list of metabolic modules, where each module is represented as a vector of metabolite IDs.
#' @param detected_metabolites A character vector of detected metabolite IDs.
#' @param hidden_metabolites A character vector of hidden metabolite IDs.
#' @param sub_metabolic_network An `igraph` object representing the sub-metabolic network.
#' @param threads An integer specifying the number of threads to use for parallel computation. Default is 3.
#'
#' @return A numeric vector where each value represents the activity score of a corresponding metabolic module.
#'
#' @details
#' The activity score is computed by assessing the modularity of each metabolic module within the sub-metabolic network.
#' The function extracts the subgraph corresponding to each module and computes a modularity-based score adjusted for the
#' total number of metabolites in the dataset.
#'
#' The computation is performed in parallel using `BiocParallel`, selecting either `MulticoreParam` (for Unix-based systems)
#' or `SnowParam` (for Windows) to handle multi-threaded execution efficiently.
#'
#' @import igraph
#' @importFrom BiocParallel bplapply MulticoreParam SnowParam
#' @importFrom masstools get_os
#'
#' @examples
#' \dontrun{
#' library(igraph)
#'
#' # Example metabolic network
#' g <- make_ring(10)
#' V(g)$name <- paste0("M", 1:10)
#'
#' # Example input data
#' metabolic_modules <- list(
#' c("M1", "M2", "M3"),
#' c("M4", "M5", "M6", "M7")
#' )
#' detected_metabolites <- c("M1", "M3", "M5", "M7")
#' hidden_metabolites <- c("M2", "M6")
#'
#' # Compute activity scores
#' scores <- calculate_activity_score(
#' metabolic_modules = metabolic_modules,
#' detected_metabolites = detected_metabolites,
#' hidden_metabolites = hidden_metabolites,
#' sub_metabolic_network = g,
#' threads = 2
#' )
#'
#' print(scores)
#' }
#'
#' @export
<-
(metabolic_modules,
detected_metabolites,
hidden_metabolites,
sub_metabolic_network,
threads = 3) {
temp_fun_activity_scores <-
(x,
detected_metabolites,
hidden_metabolites,
total_metabolite_number,
sub_metabolic_network = sub_metabolic_network) {
("igraph", quietly = )
# options(warn = -1)
temp_graph <- igraph::subgraph(graph = sub_metabolic_network, v = x)
detected_number <- (x %in% detected_metabolites)
hidden_number <- (x %in% hidden_metabolites)
m <- (igraph::E(sub_metabolic_network))
e <- (igraph::E(temp_graph))
node_name <- (igraph::V(temp_graph))
value <- (node_name, (node) {
temp_value <- (node_name, (x) {
temp1 <- igraph::degree(graph = sub_metabolic_network, v = x)
temp2 <- igraph::degree(graph = sub_metabolic_network, v = node)
temp1 * temp2 / (4 * m^2)
})
(temp_value)
})
modularity <- e / m - (value)
# modularity
<- modularity * (total_metabolite_number / (x))
a <- detected_number * / (x)
a
}
if (masstools::get_os() == "windows") {
bpparam = BiocParallel::SnowParam(workers = threads, progressbar = )
} {
bpparam = BiocParallel::MulticoreParam(workers = threads, progressbar = )
}
activity_score <-
BiocParallel::bplapply(
metabolic_modules,
temp_fun_activity_scores,
BPPARAM = bpparam,
detected_metabolites = detected_metabolites,
hidden_metabolites = hidden_metabolites,
total_metabolite_number = (detected_metabolites),
sub_metabolic_network = sub_metabolic_network
)
activity_score <- (activity_score)
(activity_score)
}
# Function to compute module quality metrics
calculate_module_quality <-
(graph, nodes) {
temp_graph <-
igraph::induced_subgraph(graph, vids = nodes)
internal_edges <- igraph::gsize(temp_graph)
external_edges <- (igraph::degree(graph, v = nodes)) - (2 * internal_edges)
# Conductance
conductance <- ((internal_edges + external_edges) == 0,
0,
external_edges / (internal_edges + external_edges)
)
# Internal Density
num_nodes <- (nodes)
internal_density <-
(num_nodes < 2, 0, (2 * internal_edges) / (num_nodes * (num_nodes - 1)))
# Edge Density Ratio (EDR)
edge_density_ratio <-
(external_edges == 0, 1, internal_edges / external_edges)
# Return metrics
(
(
module_size = num_nodes,
internal_edges = internal_edges,
external_edges = external_edges,
conductance = conductance,
internal_density = internal_density,
edge_density_ratio = edge_density_ratio
)
)
}
generate_null_conductance_distribution <-
(annotation_table_all,
feature_table_marker,
metabolite_database,
metabolic_network,
permutation_times = 5,
threads = 8,
adduct.table = ,
column = ("rp", "hilic"),
ms1.match.ppm = 25,
mz.ppm.thr = 400,
include_hidden_metabolites = ) {
column <- (column) # Ensure valid column type
# Initialize list to store activity scores for each permutation
null_conductance <- ( = "list", = permutation_times)
marker_number <- (feature_table_marker) # Number of markers
for (i (permutation_times)) {
(i, " ")
# Randomly sample markers for permutation
if ((feature_table_marker) / ((annotation_table_all$variable_id)) < 0.2) {
temp_marker_variable_id <- (
x = (
(annotation_table_all$variable_id),
feature_table_marker$variable_id
),
size = marker_number,
=
)
} {
temp_marker_variable_id <- (
x = (annotation_table_all$variable_id),
size = marker_number,
=
)
}
# Extract annotation information for sampled markers
temp_annotation_table_marker <- annotation_table_all %>%
dplyr::(variable_id %in% temp_marker_variable_id)
# Extract detected metabolites from annotation table
temp_detected_metabolites <- (temp_annotation_table_marker$KEGG.ID)
# Identify hidden metabolites connected to detected metabolites
if (include_hidden_metabolites) {
temp_hidden_metabolites <- (
metabolic_network = metabolic_network,
detected_metabolites = temp_detected_metabolites,
threads = threads,
max.reaction = 3
) %>% () %>% ()
} {
temp_hidden_metabolites <-
}
temp_total_metabolites <-
(temp_detected_metabolites, temp_hidden_metabolites)
# Extract subnetwork containing detected and hidden metabolites
temp_sub_metabolic_network <- metabolic_network %>%
tidygraph::activate("nodes") %>%
dplyr::( %in% temp_total_metabolites) %>%
dplyr::mutate(node_class = ( %in% temp_detected_metabolites, "Detected", "Hidden")) %>%
tidygraph::mutate(degree2 = tidygraph::centrality_degree())
# Identify metabolic modules in the subnetwork
temp_metabolic_modules <- (
sub_metabolic_network = temp_sub_metabolic_network,
detected_metabolites = temp_detected_metabolites,
hidden_metabolites = temp_hidden_metabolites
)
# Calculate activity scores for the identified metabolic modules
temp_conductance <-
temp_metabolic_modules %>%
purrr::map((x) {
calculate_module_quality(graph = temp_sub_metabolic_network, nodes = x)$conductance
}) %>%
()
# Store the computed activity score in the list
null_conductance[i]] <- temp_conductance
}
(null_conductance)
}
check_feature_table_all <-
(feature_table_all) {
####feature_table_all
# 1. Check if it's a data frame
if (!(feature_table_all)) {
("Error: feature_table_all must be a data frame.")
}
# 2. Check column names
expected_colnames <- ("variable_id", "mz", "rt", "polarity")
if (!((feature_table_all), expected_colnames)) {
("Error: Column names must be exactly: 'variable_id', 'mz', 'rt', 'polarity'.")
}
# 3. Check if variable_id is unique
if ((feature_table_all$variable_id) > 0) {
("Error: 'variable_id' column must have unique values.")
}
# 4. Check if mz is numeric and has no NA
if (!(feature_table_all$mz) ||
((feature_table_all$mz))) {
("Error: 'mz' column must be numeric and contain no missing values.")
}
# 5. Check if rt is numeric and has no NA
if (!(feature_table_all$) ||
((feature_table_all$))) {
("Error: 'rt' column must be numeric and contain no missing values.")
}
# 6. Check if rt is in seconds (not in minutes)
max_rt <- (feature_table_all$, na.rm = ) # Get the max retention time
if (max_rt < 100) {
# If max RT is too small, it might be in minutes
("Error: 'rt' values seem to be in minutes, but they should be in seconds.")
}
# 7. Check polarity column
valid_polarity <- ("positive", "negative")
if (((feature_table_all$polarity))) {
("Error: 'polarity' column must not contain NA values.")
}
if ((("\\s", feature_table_all$polarity))) {
("Error: 'polarity' column must not contain spaces.")
}
if (!(feature_table_all$polarity %in% valid_polarity)) {
("Error: 'polarity' must only contain 'positive' or 'negative'.")
}
# If all checks pass
("feature_table_all is correct.\n")
}
check_feature_table_marker <-
(feature_table_marker) {
# 1. Check if it's a data frame
if (!(feature_table_marker)) {
("Error: feature_table_marker must be a data frame.")
}
# 2. Check column names
expected_colnames <- ("variable_id", "mz", "rt", "degree", "p_value", "polarity")
if (!((feature_table_marker), expected_colnames)) {
(
"Error: Column names must be exactly: 'variable_id', 'mz', 'rt', 'degree', 'p_value', 'polarity'."
)
}
# 3. Check if variable_id is unique
if ((feature_table_marker$variable_id) > 0) {
("Error: 'variable_id' column must have unique values.")
}
# 4. Check if mz is numeric and has no NA
if (!(feature_table_marker$mz) ||
((feature_table_marker$mz))) {
("Error: 'mz' column must be numeric and contain no missing values.")
}
# 5. Check if rt is numeric, has no NA, and represents time in seconds
if (!(feature_table_marker$) ||
((feature_table_marker$))) {
("Error: 'rt' column must be numeric and contain no missing values.")
}
# 6. Ensure rt is in seconds (not minutes)
max_rt <- (feature_table_marker$, na.rm = ) # Get the max retention time
if (max_rt < 100) {
# If max RT is too small, it might be in minutes
("Error: 'rt' values seem to be in minutes, but they should be in seconds.")
}
# 7. Check if degree is numeric and has no NA
if (!(feature_table_marker$degree) ||
((feature_table_marker$degree))) {
("Error: 'degree' column must be numeric and contain no missing values.")
}
# 8. Check if p_value is numeric and has no NA
if (!(feature_table_marker$p_value) ||
((feature_table_marker$p_value))) {
("Error: 'p_value' column must be numeric and contain no missing values.")
}
# 9. Check polarity column
valid_polarity <- ("positive", "negative")
if (((feature_table_marker$polarity))) {
("Error: 'polarity' column must not contain NA values.")
}
if ((("\\s", feature_table_marker$polarity))) {
("Error: 'polarity' column must not contain spaces.")
}
if (!(feature_table_marker$polarity %in% valid_polarity)) {
("Error: 'polarity' must only contain 'positive' or 'negative'.")
}
# If all checks pass
("feature_table_marker is correct.")
}
# Example usage:
# check_feature_table_marker(feature_table_marker)
