## ---- echo=FALSE-------------------------------------------------------------- knitr::opts_chunk$set(fig.align = "center") ## ----eval=FALSE--------------------------------------------------------------- # # install the package via BioConductor # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install("netOmics") ## ----eval=FALSE--------------------------------------------------------------- # # install the package via github # library(devtools) # install_github("abodein/netOmics") ## ---- eval=TRUE, message=FALSE------------------------------------------------ # load the package library(netOmics) ## ---- eval=TRUE, message=FALSE------------------------------------------------ # usefull packages to build this vignette library(timeOmics) library(tidyverse) library(igraph) ## ----load_data---------------------------------------------------------------- # load data data("hmp_T2D") ## ----timeOmics_1, eval=FALSE-------------------------------------------------- # # not evaluated in this vignette # # #1 filter fold-change # remove.low.cv <- function(X, cutoff = 0.5){ # # var.coef # cv <- unlist(lapply(as.data.frame(X), # function(x) abs(sd(x, na.rm = TRUE)/mean(x, na.rm= TRUE)))) # return(X[,cv > cutoff]) # } # fc.threshold <- list("RNA"= 1.5, "CLINICAL"=0.2, "GUT"=1.5, "METAB"=1.5, # "PROT" = 1.5, "CYTO" = 1) # # # --> hmp_T2D$raw # data.filter <- imap(raw, ~{remove.low.cv(.x, cutoff = fc.threshold[[.y]])}) # # #2 scale # data <- lapply(data.filter, function(x) log(x+1)) # # --> hmp_T2D$data # # # #3 modelling # lmms.func <- function(X){ # time <- rownames(X) %>% str_split("_") %>% # map_chr(~.x[[2]]) %>% as.numeric() # lmms.output <- lmms::lmmSpline(data = X, time = time, # sampleID = rownames(X), deri = FALSE, # basis = "p-spline", numCores = 4, # keepModels = TRUE) # return(lmms.output) # } # data.modelled <- lapply(data, function(x) lmms.func(x)) # # # 4 clustering # block.res <- block.pls(data.modelled, indY = 1, ncomp = 1) # getCluster.res <- getCluster(block.res) # # --> hmp_T2D$getCluster.res # # # # 5 signature # list.keepX <- list("CLINICAL" = 4, "CYTO" = 3, "GUT" = 10, "METAB" = 3, # "PROT" = 2,"RNA" = 34) # sparse.block.res <- block.spls(data.modelled, indY = 1, ncomp = 1, scale =TRUE, # keepX =list.keepX) # getCluster.sparse.res <- getCluster(sparse.block.res) # # --> hmp_T2D$getCluster.sparse.res ## ----timeOmics_2-------------------------------------------------------------- # clustering results cluster.info <- hmp_T2D$getCluster.res ## ----graph.rna, warning=FALSE------------------------------------------------- cluster.info.RNA <- timeOmics::getCluster(cluster.info, user.block = "RNA") graph.rna <- get_grn(X = hmp_T2D$data$RNA, cluster = cluster.info.RNA) # to get info about the network get_graph_stats(graph.rna) ## ----PROT_graph, warning=FALSE------------------------------------------------ # Utility function to get the molecules by cluster get_list_mol_cluster <- function(cluster.info, user.block){ require(timeOmics) tmp <- timeOmics::getCluster(cluster.info, user.block) res <- tmp %>% split(.$cluster) %>% lapply(function(x) x$molecule) res[["All"]] <- tmp$molecule return(res) } cluster.info.prot <- get_list_mol_cluster(cluster.info, user.block = 'PROT') graph.prot <- get_interaction_from_database(X = cluster.info.prot, db = hmp_T2D$interaction.biogrid, type = "PROT", user.ego = TRUE) # get_graph_stats(graph.prot) ## ----GUT_graph, eval = FALSE-------------------------------------------------- # # not evaluated in this vignette # library(SpiecEasi) # # get_sparcc_graph <- function(X, threshold = 0.3){ # res.sparcc <- sparcc(data = X) # sparcc.graph <- abs(res.sparcc$Cor) >= threshold # colnames(sparcc.graph) <- colnames(X) # rownames(sparcc.graph) <- colnames(X) # res.graph <- graph_from_adjacency_matrix(sparcc.graph, # mode = "undirected") %>% simplify # return(res.graph) # } # # gut_list <- get_list_mol_cluster(cluster.info, user.block = 'GUT') # # graph.gut <- list() # graph.gut[["All"]] <- get_sparcc_graph(hmp_T2D$raw$GUT, threshold = 0.3) # graph.gut[["1"]] <- get_sparcc_graph(hmp_T2D$raw$GUT %>% # dplyr::select(gut_list[["1"]]), # threshold = 0.3) # graph.gut[["-1"]] <- get_sparcc_graph(hmp_T2D$raw$GUT %>% # dplyr::select(gut_list[["-1"]]), # threshold = 0.3) # class(graph.gut) <- "list.igraph" ## ----GUT---------------------------------------------------------------------- graph.gut <- hmp_T2D$graph.gut # get_graph_stats(graph.gut) ## ----CYTO_graph, warning=FALSE------------------------------------------------ # CYTO -> from database (biogrid) cyto_list = get_list_mol_cluster(cluster.info = cluster.info, user.block = "CYTO") graph.cyto <- get_interaction_from_database(X = cyto_list, db = hmp_T2D$interaction.biogrid, type = "CYTO", user.ego = TRUE) # get_graph_stats(graph.cyto) # METAB -> inference cluster.info.metab <- timeOmics::getCluster(X = cluster.info, user.block = "METAB") graph.metab <- get_grn(X = hmp_T2D$data$METAB, cluster = cluster.info.metab) # get_graph_stats(graph.metab) # CLINICAL -> inference cluster.info.clinical <- timeOmics::getCluster(X = cluster.info, user.block = 'CLINICAL') graph.clinical <- get_grn(X = hmp_T2D$data$CLINICAL, cluster = cluster.info.clinical) # get_graph_stats(graph.clinical) ## ---- merged_0---------------------------------------------------------------- full.graph <- combine_layers(graph1 = graph.rna, graph2 = graph.prot) full.graph <- combine_layers(graph1 = full.graph, graph2 = graph.cyto) full.graph <- combine_layers(graph1 = full.graph, graph2 = hmp_T2D$interaction.TF) # get_graph_stats(full.graph) ## ----merged_1_gut, warning=FALSE---------------------------------------------- all_data <- reduce(hmp_T2D$data, cbind) # omic = gut gut_list <- get_list_mol_cluster(cluster.info, user.block = "GUT") omic_data <- lapply(gut_list, function(x)dplyr::select(hmp_T2D$data$GUT, x)) # other data = "RNA", "PROT", "CYTO" other_data_list <- get_list_mol_cluster(cluster.info, user.block = c("RNA", "PROT", "CYTO")) other_data <- lapply(other_data_list, function(x)dplyr::select(all_data, x)) # get interaction between gut data and other data interaction_df_gut <- get_interaction_from_correlation(X = omic_data, Y = other_data, threshold = 0.99) # and merge with full graph full.graph <- combine_layers(graph1 = full.graph, graph2 = graph.gut, interaction.df = interaction_df_gut$All) ## ---- merged_2_clinical, warning=FALSE---------------------------------------- # omic = Clinical clinical_list <- get_list_mol_cluster(cluster.info, user.block = "CLINICAL") omic_data <- lapply(clinical_list, function(x)dplyr::select(hmp_T2D$data$CLINICAL, x)) # other data = "RNA", "PROT", "CYTO", "GUT" other_data_list <- get_list_mol_cluster(cluster.info, user.block = c("RNA", "PROT", "CYTO", "GUT")) other_data <- lapply(other_data_list, function(x)dplyr::select(all_data, x)) # get interaction between gut data and other data interaction_df_clinical <- get_interaction_from_correlation(X = omic_data , Y = other_data, threshold = 0.99) # and merge with full graph full.graph <- combine_layers(graph1 = full.graph, graph2 = graph.clinical, interaction.df = interaction_df_clinical$All) ## ---- merged_3_metab, warning=FALSE------------------------------------------- # omic = Metab metab_list <- get_list_mol_cluster(cluster.info, user.block = "METAB") omic_data <- lapply(metab_list, function(x)dplyr::select(hmp_T2D$data$METAB, x)) # other data = "RNA", "PROT", "CYTO", "GUT", "CLINICAL" other_data_list <- get_list_mol_cluster(cluster.info, user.block = c("RNA", "PROT", "CYTO", "GUT", "CLINICAL")) other_data <- lapply(other_data_list, function(x)dplyr::select(all_data, x)) # get interaction between gut data and other data interaction_df_metab <- get_interaction_from_correlation(X = omic_data, Y = other_data, threshold = 0.99) # and merge with full graph full.graph <- combine_layers(graph1 = full.graph, graph2 = graph.metab, interaction.df = interaction_df_metab$All) ## ----------------------------------------------------------------------------- # ORA by cluster/All mol_ora <- get_list_mol_cluster(cluster.info, user.block = c("RNA", "PROT", "CYTO")) # get ORA interaction graph by cluster graph.go <- get_interaction_from_ORA(query = mol_ora, sources = "GO", organism = "hsapiens", signif.value = TRUE) # merge full.graph <- combine_layers(graph1 = full.graph, graph2 = graph.go) ## ----------------------------------------------------------------------------- # medlineRanker -> database medlineranker.res.df <- hmp_T2D$medlineranker.res.df %>% dplyr::select(Disease, symbol) %>% set_names(c("from", "to")) mol_list <- get_list_mol_cluster(cluster.info = cluster.info, user.block = c("RNA", "PROT", "CYTO")) graph.medlineranker <- get_interaction_from_database(X = mol_list, db = medlineranker.res.df, type = "Disease", user.ego = TRUE) # get_graph_stats(graph.medlineranker) # merging full.graph <- combine_layers(graph1 = full.graph, graph2 = graph.medlineranker) ## ----------------------------------------------------------------------------- # graph cleaning graph_cleaning <- function(X, cluster.info){ # no reusability X <- igraph::simplify(X) va <- vertex_attr(X) viewed_mol <- c() for(omic in unique(cluster.info$block)){ mol <- intersect(cluster.info %>% dplyr::filter(.$block == omic) %>% pull(molecule), V(X)$name) viewed_mol <- c(viewed_mol, mol) X <- set_vertex_attr(graph = X, name = "type", index = mol, value = omic) X <- set_vertex_attr(graph = X, name = "mode", index = mol, value = "core") } # add medline ranker and go mol <- intersect(map(graph.go, ~ as_data_frame(.x)$to) %>% unlist %>% unique(), V(X)$name) # only GO terms viewed_mol <- c(viewed_mol, mol) X <- set_vertex_attr(graph = X, name = "type", index = mol, value = "GO") X <- set_vertex_attr(graph = X, name = "mode", index = mol, value = "extended") mol <- intersect(as.character(medlineranker.res.df$from), V(X)$name) viewed_mol <- c(viewed_mol, mol) X <- set_vertex_attr(graph = X, name = "type", index = mol, value = "Disease") X <- set_vertex_attr(graph = X, name = "mode", index = mol, value = "extended") other_mol <- setdiff(V(X), viewed_mol) if(!is_empty(other_mol)){ X <- set_vertex_attr(graph = X, name = "mode", index = other_mol, value = "extended") } X <- set_vertex_attr(graph = X, name = "mode", index = intersect(cluster.info$molecule, V(X)$name), value = "core") # signature mol <- intersect(V(X)$name, hmp_T2D$getCluster.sparse.res$molecule) X <- set_vertex_attr(graph = X, name = "sparse", index = mol, value = TRUE) mol <- setdiff(V(X)$name, hmp_T2D$getCluster.sparse.res$molecule) X <- set_vertex_attr(graph = X, name = "sparse", index = mol, value = FALSE) return(X) } ## ----------------------------------------------------------------------------- FULL <- lapply(full.graph, function(x) graph_cleaning(x, cluster.info)) get_graph_stats(FULL) ## ---- eval = FALSE------------------------------------------------------------ # # degree analysis # d <- degree(FULL$All) # hist(d) # d[max(d)] # # # modularity # Warnings: can take several minutes # res.mod <- walktrap.community(FULL$All) # # ... # # # modularity # sp <- shortest.paths(FULL$All) ## ----------------------------------------------------------------------------- # seeds = all vertices -> takes 5 minutes to run on regular computer # seeds <- V(FULL$All)$name # rwr_res <- random_walk_restart(FULL, seeds) # seed = some GO terms seeds <- head(V(FULL$All)$name[V(FULL$All)$type == "GO"]) rwr_res <- random_walk_restart(FULL, seeds) ## ----------------------------------------------------------------------------- rwr_type_k15 <- rwr_find_seeds_between_attributes(X = rwr_res, attribute = "type", k = 15) # a summary plot function summary_plot_rwr_attributes(rwr_type_k15) summary_plot_rwr_attributes(rwr_type_k15$All) ## ----------------------------------------------------------------------------- rwr_type_k15 <- rwr_find_seeds_between_attributes(X = rwr_res$All, attribute = "cluster", k = 15) summary_plot_rwr_attributes(rwr_type_k15) ## ----------------------------------------------------------------------------- sub_res <- rwr_type_k15$`GO:0005737` sub <- plot_rwr_subnetwork(sub_res, legend = TRUE, plot = TRUE) ## ----------------------------------------------------------------------------- rwr_res <- random_walk_restart(FULL$All, seed = "ZNF263") # closest GO term rwr_find_closest_type(rwr_res, seed = "ZNF263", attribute = "type", value = "GO", top = 5) # closest Disease rwr_find_closest_type(rwr_res, seed = "ZNF263", attribute = "type", value = "Disease", top = 5) # closest nodes with an attribute "cluster" and the value "-1" rwr_find_closest_type(rwr_res, seed = "ZNF263", attribute = "cluster", value = "-1", top = 5) ## ---- eval = FALSE------------------------------------------------------------ # seeds <- V(FULL$All)$name[V(FULL$All)$type %in% c("GO", "Disease")] ## ----------------------------------------------------------------------------- sessionInfo()