## ----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()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.