diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..63d1874 Binary files /dev/null and b/.DS_Store differ diff --git a/Scripts/delta_list.R b/Scripts/delta_list.R new file mode 100644 index 0000000..9f299aa --- /dev/null +++ b/Scripts/delta_list.R @@ -0,0 +1,88 @@ +to_matrix <- function(df) { + + nm <- df[[1]] + + df %>% select(-1) %>% as.matrix() %>% set_rownames(nm) + +} + +my_quick_mass <- function(x) { + check_chemform(isotopes = isotopes, x)$monoisotopic_mass +} + +K = my_quick_mass("K") +Na = my_quick_mass("Na") +H = my_quick_mass("H") +C13 = 1.003355 + +NH3 = my_quick_mass("NH3") +H2O = my_quick_mass("H2O") + + +delta_list = + c(`Na-H` = Na - H, + `K-H` = K-H, + NH3 = NH3, + `Na-NH4` = Na - NH3 - H, + `K-NH4` = K - NH3 - H, + H2O = H2O, + C13 = C13, + C1H2O2 = my_quick_mass("C1H2O2"), + C2H7N1 = my_quick_mass("C2H7N1"), + C8H6O3 = my_quick_mass("C8H6O3"), + CO2 = my_quick_mass("C1O2"), + Bz = my_quick_mass("C7H5O1"), + Acetonitrile = my_quick_mass("C2H3N1"), + DMSO = my_quick_mass("C2H6O1S1") + ) + + +d1 <- tibble( + adduct = c("H", "Na", "K", "NH4"), + delta1 = c(0, Na-H, K-H, NH3) +) + + + +d2 <- tibble( + gain_loss = c("None", "+DMSO", "+Acetonitrile", "+C2H7N1", "+H2O", "-H2O", "-CO2", "-C1H2O2", "-C8H6O3", "-Bz", "-NH3"), + delta2 = c(0, delta_list[c("DMSO", "Acetonitrile", "C2H7N1", "H2O")], -delta_list[c("H2O", "CO2", "C1H2O2", "C8H6O3", "Bz", "NH3")]) +) + +d3 <- tibble( + isotope = c("M", "M+1", "M+2", "M+3"), + delta3 = C13 * c(0:3) +) + + +format_species <- function(a, gl, iso) { + + case_when( + gl == "None" & iso == "M" ~ sprintf("[%s + %s]", iso, a), + gl == "None" & iso != "M" ~ sprintf("[(%s) + %s]", iso, a), + gl != "None" & iso == "M" ~ sprintf("[%s + %s] %s", iso, a, gl), + TRUE ~ sprintf("[(%s) + %s] %s", iso, a, gl) + ) + +} + +delta_list2 <- + expand_grid( + adduct = c("H", "Na", "K", "NH4"), + gain_loss = c("None", "+DMSO", "+Acetonitrile", "+C2H7N1", "+H2O", "-H2O", "-CO2", "-C1H2O2", "-C8H6O3", "-Bz", "-NH3"), + isotope = c("M", "M+1", "M+2", "M+3") + ) %>% + left_join( + d1, by = "adduct" + ) %>% + left_join( + d2, by = "gain_loss" + ) %>% + left_join( + d3, by = "isotope" + ) %>% + filter(!(gain_loss == "-NH3" & adduct == "NH4")) %>% + mutate( + name2 = format_species(adduct, gain_loss, isotope), + delta = delta1 + delta2 + delta3 + ) diff --git a/Scripts/get_network_layout.R b/Scripts/get_network_layout.R new file mode 100644 index 0000000..58ec930 --- /dev/null +++ b/Scripts/get_network_layout.R @@ -0,0 +1,45 @@ +get_network_layout <- function(obj, clust) { + + # node_cat, edge_cat, node_size, edge_size + + node <- obj$node %>% + filter(clust_id == clust) + + edge <- obj$edge %>% + filter(FEATURE_ID_1 %in% node$FEATURE_ID & FEATURE_ID_2 %in% node$FEATURE_ID) + + coord <- get_layout(edge) + + node %<>% + left_join(coord, by = "FEATURE_ID") + + edge %<>% + left_join(coord, by = c("FEATURE_ID_1" = "FEATURE_ID")) %>% + left_join(coord, by = c("FEATURE_ID_2" = "FEATURE_ID"), suffix = c("", "end")) + + deg <- c(edge$FEATURE_ID_1, edge$FEATURE_ID_2) %>% table() %>% {tibble(FEATURE_ID = names(.), Degree = unname(.))} + + node %<>% + left_join(deg, by = "FEATURE_ID") + + list( + edge = edge, + node = node + ) + +} + + + + + + + + + + + + + + + diff --git a/Scripts/network_plotter.R b/Scripts/network_plotter.R new file mode 100644 index 0000000..faf69f5 --- /dev/null +++ b/Scripts/network_plotter.R @@ -0,0 +1,71 @@ + +gg_network_plotter <- function( + obj, + node_size1 = 10, + node_colour_var = "adduct", + show_text = F, + text_colour_var = "name2", + text_size = 8, + text_colour = "cyan", + edge_colour_var = "Signature", + edge_size = 1 +) + +{ + + node <- obj$node + edge <- obj$edge + + save(list = ls(), file = "bip.RData") + + gg <- ggplot() + + geom_segment( + mapping = aes(x=x, y=y, xend=xend, yend=yend, colour = !!as.name(edge_colour_var)), + data = edge, + size = edge_size + ) + + + geom_point( + data = node, + mapping=aes(x=x, y=y, fill = !!as.name(node_colour_var)), + shape = 21, + size = node_size1, + ) + + + theme_void() + + labs( + title = "COMPOUND_ID %s MZ %.4f RT %.2f" %>% sprintf(node$COMPOUND_ID[1], node$MZ[1], node$RT[1]) + ) + + if(show_text) { + + gg <- gg + + geom_text( + data = node, + mapping=aes(x=x, y=y, label = !!as.name(text_colour_var)), + size = text_size, + colour = text_colour, + ) + } + + + return(gg) + +} + + +plot_cluster <- function(obj, clust) { + + obj %>% + get_network_layout(clust) %>% + gg_network_plotter(show_text = T) + +} + + + + + + + + diff --git a/Scripts/propagate_label.R b/Scripts/propagate_label.R new file mode 100644 index 0000000..c6968e8 --- /dev/null +++ b/Scripts/propagate_label.R @@ -0,0 +1,128 @@ +delta_converter_list <- read_csv("delta_converter_list.csv", show_col_types = FALSE) %>% to_matrix() + +my_label_propagation2 <- function(obj) { + + node <- obj$node %>% + arrange(-mean) + + edge <- obj$edge + + # anchor_node: character[1] (FEATURE_ID) + anchor_node <- node$FEATURE_ID %>% first() + + # initiliase tibbles + nodelist <- anchor_node + done <- character(0) + + + # matrix + node_assign <- + bind_cols( + node %>% select(FEATURE_ID), + map_df(as_tibble(delta_converter_list), ~ rep(0L, nrow(node))) + ) %>% + to_matrix() + + # assume principal node is M+H + node_assign[anchor_node, "H"] <- 1 + + # neighborhood tibble + N <- get_adjacency_tibble(edge) + + + while (length(nodelist) > 0) { + + new_nodelist <- character(0) + + # i: character[1] + for (i in nodelist) { + + # print(paste("i", i)) + + # j character[1] + for(j in N %>% filter(n1==i) %$% n2 %>% setdiff(done)) { + + # print(paste("j", j)) + + node_assign[j,] <- node_assign[i,] + convert_delta_code(N %>% filter(n1==i, n2==j) %$% code) + + new_nodelist <- c(new_nodelist, j) + } + + done <- c(done, i) + } + + done %<>% unique() + + nodelist <- new_nodelist %>% unique() + + + + } + + node_assign %<>% + as_tibble(rownames = "FEATURE_ID") %>% + get_species() + + list( + edge = edge, + node = node %>% left_join(node_assign %>% select(FEATURE_ID, M, adduct, gain_loss, name2), by = "FEATURE_ID"), + decomposition_table = nodelist, + anchor_node = anchor_node + + ) + +} + +get_species <- function(x) { + + + x %>% + mutate( + adduct = sprintf("_Na_%d__K_%d__H_%d_", Na, K, H) %>% + gsub("_[A-Za-z0-9]+_0_", "", .) %>% + gsub("(_1)|(^_)|(_$)", "", .) %>% + # sub("H__NH3", "NH4", .) %>% + sub("__H_-([1-9])", "-\\1H", .) %>% + sub("-1H", "-H", .) %>% + gsub("__", "+", .), + + M = if_else(C13 == 0, "M", sprintf("(M+%d)", C13)), + + gain_loss = sprintf("_NH3_%d__H2O_%d__FA_%d__C2H7N1_%d__C8H6O3_%d__CO2_%d__Bz_%d__ACN_%d__DMSO_%d_", NH3, H2O, C1H2O2, C2H7N1, C8H6O3, CO2, Bz, Acetonitrile, DMSO) %>% + + gsub("_[A-Za-z0-9]+_0_", "", .) %>% + gsub("_([A-Za-z0-9]+)_-([1-9])", "-\\2.\\1", .) %>% + gsub("-1\\.", "-", .) %>% + gsub("(_1)|(^_)|(_$)", "", .) %>% + gsub("__", "+", .) %>% + gsub("_-", "-", .), + + name2 = sprintf("[%s + %s] +%s", M, adduct, gain_loss) %>% + gsub("\\+-", "-", .) %>% + gsub("\\+$", "", .) %>% + trimws() + + + ) + +} + + + +get_adjacency_tibble <- function(edge) { + + adj <- get_adjacency_matrix_delta(edge) %>% t() + + N <- lapply(seq(nrow(adj)), function(i) adj[i,] %>% keep(.!=0)) + + tibble( + n1 = rownames(adj), + n2 = map(N, ~names(.x)), + code = map(N, ~unname(.x)) + ) %>% + unnest(c(n2, code)) + +} + + diff --git a/cor_cluster_mix.RDS b/cor_cluster_mix.RDS new file mode 100644 index 0000000..1089330 Binary files /dev/null and b/cor_cluster_mix.RDS differ diff --git a/data20.RDS b/data20.RDS new file mode 100644 index 0000000..509d4e5 Binary files /dev/null and b/data20.RDS differ diff --git a/delta_converter_list.csv b/delta_converter_list.csv new file mode 100644 index 0000000..f3c7b0d --- /dev/null +++ b/delta_converter_list.csv @@ -0,0 +1,15 @@ +Name,Na,K,H,NH3,H2O,C13,C1H2O2,C2H7N1,C8H6O3,CO2,Bz,Acetonitrile,DMSO +Na-H,1,0,-1,0,0,0,0,0,0,0,0,0,0 +K-H,0,1,-1,0,0,0,0,0,0,0,0,0,0 +NH3,0,0,0,1,0,0,0,0,0,0,0,0,0 +Na-NH4,1,0,-1,-1,0,0,0,0,0,0,0,0,0 +K-NH4,0,1,-1,-1,0,0,0,0,0,0,0,0,0 +H2O,0,0,0,0,1,0,0,0,0,0,0,0,0 +C13,0,0,0,0,0,1,0,0,0,0,0,0,0 +C1H2O2,0,0,0,0,0,0,1,0,0,0,0,0,0 +C2H7N1,0,0,0,0,0,0,0,1,0,0,0,0,0 +C8H6O3,0,0,0,0,0,0,0,0,1,0,0,0,0 +CO2,0,0,0,0,0,0,0,0,0,1,0,0,0 +Bz,0,0,0,0,0,0,0,0,0,0,1,0,0 +Acetonitrile,0,0,0,0,0,0,0,0,0,0,0,1,0 +DMSO,0,0,0,0,0,0,0,0,0,0,0,0,1 diff --git a/make_data20.R b/make_data20.R new file mode 100644 index 0000000..7e97a8b --- /dev/null +++ b/make_data20.R @@ -0,0 +1,7 @@ + +cor_cluster_mix <- readRDS("cor_cluster_mix.RDS") + +cl <- cor_cluster_mix$node %>% arrange(-clust_size) %$% clust_id %>% unique() %>% head(100) + +purrr::map(cl %>% head(20), ~get_network_layout(cor_cluster_mix, .x) %>% my_label_propagation2()) %>% + saveRDS("data20.RDS") diff --git a/shiny_uNetwork.R b/shiny_uNetwork.R new file mode 100644 index 0000000..ea14ad1 --- /dev/null +++ b/shiny_uNetwork.R @@ -0,0 +1,386 @@ +rm(list=ls(envir = .GlobalEnv), envir = .GlobalEnv) + +library(magrittr) +library(tidyverse) +library(igraph) +library(shiny) +library(shinythemes) +library(Matrix) +library(readxl) +library(enviPat) + +data("isotopes") + + +dir1 <- "Scripts" + +for (f in list.files(dir1, pattern = "\\.R$")) { + source(file.path(dir1, f)) + +} + +# source("Scripts/make_network.R") +# source("Scripts/delta_list.R") +# source("Scripts/network_plotter.R") + + +# source("Scripts/calc_layout.R") +# source("Scripts/network_bipartite.R") +# source("Scripts/convert_to_adjacency.R") +# +assign("colours1", colours() %>% keep(!grepl("^gr[ea]y", .) & !grepl("[0-9]$", .)) %>% sort(), envir = .GlobalEnv) + + + +ui <- + + fluidPage( + title = "Feature matrix plotter", + sidebarLayout( + sidebarPanel( + h4("Inputs"), + br(), + tabsetPanel( + tabPanel( + title = "Cluster", + br(), + + uiOutput("select_cluster") + ), + + tabPanel( + title = "Nodes", + br(), + + sliderInput( + inputId = "node_size1", + label = "Node size", + value = 10, + min = 1, + step = 1, + max = 20 + ), + + + selectizeInput( + inputId = "node_colour_var", + label = "Choose node colour variable", + choices = c("adduct", "gain_loss"), + selected = c("adduct", "gain_loss")[1], + multiple = F, + width = "200px" + ) + + + ), # end node + + tabPanel( + title = "Edges", + br(), + + sliderInput( + inputId = "edge_size", + label = "Edge size (thickness)", + value = 1, + min = 0.5, + step = 0.5, + max = 5 + ), + + selectizeInput( + inputId = "edge_colour_var", + label = "Choose node colour variable", + choices = c("Signature", "COR"), + selected = c("Signature", "COR")[1], + multiple = F, + width = "200px" + ) + + + ), # end edge + + tabPanel( + title = "Label", + br(), + + radioButtons( + inputId = "show_text", + label = "Show labels", + choices = c("yes", "no"), + selected = "yes", + inline = T + ), + + sliderInput( + inputId = "text_size", + label = "Label size", + value = 5, + min = 1, + step = 0.5, + max = 10 + ), + + selectizeInput( + inputId = "text_var", + label = "Choose node colour variable", + choices = c("FEATURE_ID", "name2"), + selected = "name2", + multiple = F, + width = "200px" + ), + + selectizeInput( + inputId = "text_colour", + label = "Label colour", + choices = colours1, + selected = "black", + multiple = F, + width = "200px" + ) + + + ), # end text + + tabPanel( + title = "Canvas", + br(), + + sliderInput( + inputId = "plot_width", + label = "Plot width (pixels)", + value = 600, + min = 500, + step = 100, + max = 1500 + ), + + sliderInput( + inputId = "plot_height", + label = "Plot height (pixels)", + value = 600, + min = 500, + step = 100, + max = 1500 + ) + + ), # end canvas + + + + tabPanel( + title = "Pdf", + br(), + + sliderInput( + inputId = "plot_width_cm", + label = "Pdf plot widht (cm)", + value = 12, + min = 10, + step = 1, + max = 20 + ), + + sliderInput( + inputId = "plot_height_cm", + label = "Pdf plot height (cm)", + value = 12, + min = 10, + step = 1, + max = 20 + ), + + br(), + + downloadButton("download_plot", "Save plot as pdf") + + ), # end pdf + + tabPanel( + title = "Layout", + br(), + + # selectizeInput( + # inputId = "layout_method", + # label = "Igraph layout algorithm", + # choices = igraph_layout_algos, + # selected = "layout_nicely", + # multiple = F, + # width = "200px" + # ), + + + br(), + + # fileInput( + # inputId = "matrix_infile", + # label = "Upload feature matrix" + # ), + + br(), + + actionButton("recalc_layout", "Recalculate layout") + + ) # end layout + + ) + + + + ), + mainPanel( + br(), + + uiOutput("main_header"), + + br(), + uiOutput("bipartite_plot_wrapper"), + + br(), + + tableOutput("plot_click") + + ) + ) + + + ) + + +server <- function(input, output, session) { + + my_plot <- reactiveValues() + + my_data <- readRDS("data20.RDS") + + + output$select_cluster <- renderUI({ + + req(my_data) + + ch <- names(my_data) + + nm <- sprintf("%s (%s)", ch, map_int(my_data, ~nrow(.x$node))) + map_int(my_data, ~nrow(.x$node)) + + selectizeInput( + inputId = "cluster", + label = "Choose cluster", + choices = ch %>% set_names(nm), + selected = ch[1], + multiple = F, + width = "200px" + ) + + }) + + observeEvent( + input$cluster, + { + my_plot$data <- my_data[[input$cluster]] + } + ) + + output$bipartite_plot_wrapper <- renderUI({ + + plotOutput( + "bipartite_plot", + click = "plot_click", + dblclick ="plot_dbl_click", + hover = "plot_hover", + brush = "plot_brush", + height = ifelse(is.null(input$plot_height), 600, input$plot_height), + width = ifelse(is.null(input$plot_width), 600, input$plot_width) + ) + }) + + + + + output$bipartite_plot <- renderPlot({ + + req(my_plot$data) + + my_plot$gg <- + gg_network_plotter( + my_plot$data, + node_size1 = input$node_size1, + node_colour_var = input$node_colour_var, + show_text = input$show_text == "yes", + text_colour_var = input$text_var, + text_size = input$text_size, + text_colour = input$text_colour, + edge_colour_var = input$edge_colour_var, + edge_size = input$edge_size + ) + + my_plot$gg + + }) + + output$plot_click <- renderTable({ + + req(input$plot_click) + + nearPoints( + my_plot$data$node, + input$plot_click, + maxpoints = 1, + threshold = 300 + ) + }) + + # observeEvent( + # input$plot_dbl_click, + # { + # my_plot$highlight_nodes <- NA + # }) + + output$plot_hover <- renderTable({ + + req(input$plot_hover) + + nearPoints( + my_plot$data$node, + input$plot_hover, + maxpoints = 1, + threshold = 100 + ) + + }) + + output$download_plot <- downloadHandler( + + filename = "fbf_bipartite_plot.pdf", + content = function(file) { + ggsave(file, height = input$plot_height_cm, width = input$plot_width_cm, units = "cm") + } + ) + + + onStop(function() { + + env = parent.env(environment()) + + ls(envir = env) %>% + ## .. excluding e.g. session" and "output" + setdiff(c("session", "output", "my_pref_update")) %>% + ## object in a list + mget(envir = env) %>% + ## keep only reactive values + subset(., subset = sapply(., is.reactivevalues)) %>% + ## convert reactive values + sapply(function(z) isolate(reactiveValuesToList(z)), simplify=F) %>% + ## assign to global environment + list2env(envir = .GlobalEnv) + + + cat("Session stopped\n") + cat("Type `save(list=ls(), file='tmp.RData')` to save state variables in your current working directory\n") + cat("Find out what current working directory by typing `getwd()`\n") + }) + + +} + + +shinyApp(ui, server) \ No newline at end of file