diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e031132 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +tables/ +*.png +*.pdf +*.jpeg +*.jpg +*.svg diff --git a/R/.DS_Store b/R/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/R/.DS_Store differ diff --git a/R/network_multiRegion.R b/R/network_multiRegion.R new file mode 100644 index 0000000..1466721 --- /dev/null +++ b/R/network_multiRegion.R @@ -0,0 +1,532 @@ +# module UI function +networkMultiUI <- function(id, label = "Network Multiple Regions"){ + + ns <- NS(id) + + tagList(sidebarLayout( + + sidebarPanel( + pickerInput( + inputId = ns("network_multireg"), + label = "Select brain regions", + choices = c("Amygdala" = "AMY", + "Cerebellum" = "CER", + "Dorsal DG of Hippocampus" = "dDG", + "Dorsal CA1 of Hippocampus" = "dCA1", + "Prefrontal Cortex" = "PFC", + "PVN of Hypothalamus" = "PVN", + "Ventral DG of Hippocampus" = "vDG", + "Ventral CA1 of Hippocampus" = "vCA1"), + options = list( + `actions-box` = TRUE, + size = 8, + `selected-text-format` = "count > 3" + ), + multiple = TRUE + ), + searchInput( + inputId = ns("network_gene"), + label = "Enter comma-separated list of genes: ", + placeholder = "e.g. Nr3c1,Fkbp5", + btnSearch = icon("search", class = "fas fa-search", lib = "font-awesome"), + btnReset = icon("remove", class = "fas fa-remove", lib = "font-awesome"), + width = "100%" + ), + fileInput( + inputId = ns("file1"), + label = "Upload file with list of genes", + accept = "text/plain", + buttonLabel = "Upload", + placeholder = "Upload file with list of genes"), + radioButtons( + ns("network_type"), + label = "Select type of network to display:", + choiceNames = list("Vehicle", "Differential", "Treatment"), + choiceValues = list("baseline", "diff", "treatment"), + selected = "diff" + ), + radioButtons( + ns("vis_neighbours"), + label = "Include gene neighbourhood", + choiceNames = c("Yes", "No"), + choiceValues = c(TRUE, FALSE), + selected = FALSE + ), + radioButtons( + ns("id_type"), + label = "Select type of gene ID for plot and table:", + choices = list("Ensembl", "Gene Symbol"), + selected = "Gene Symbol" + ), + checkboxGroupInput( + ns("tableContent"), + label = "Table content (diff. network):", + choices = list( + "z-scores" = "z", + "p-values" = "p", + "beta values" = "beta" + ), + selected = "z" + ), + downloadButton(ns("download2"),"Download (filtered) table as csv"), + width = 3 + ), + mainPanel( + # visNetwork + fluidPage( + br(), + tags$style(".fa-project-diagram {color:#2980B9}"), + h3(p( + em("Network analysis of gene expression "), + icon("project-diagram", lib = "font-awesome"), + style = "color:black;text-align:center" + )), + hr(), + h4(span(textOutput(ns("network_type")), style = "color:black;text-align:center")), + strong("Visualize"),span(" a network in "), strong("multiple brain regions"), + span(" after vehicle treatment, Dexamethasone treatment or on a differential level. + Please select the brain regions you are interested in on the left panel. + You can enter your "), strong("genes of interest"), span(" either as a comma-separated list or upload a file + with the respective gene IDs in gene symbol or ENSEMBL format. Furthermore, you have the + option to "), strong("display neighbouring genes"), span(" and the respective connections in addition to + your genes of interest. Gene IDs can be displayed in the network as gene symbol or + ENSEMBL id. You have the option to "), strong("filter and download"), span(" the + network data table."), + br(),br(), + span(em("Hint: If no network shows up without including the neighbourhood, there are + probably no direct connections between the genes you entered. Please try to + include the gene neighbourhood."),br(), + em("However, if you include the gene neighbourhood, networks become quickly very large + and computationally expensive . We recommend to include the gene neighbourhood + only for up to 4 input genes.")), + br(),br(), + visNetworkOutput(ns("network_diff_multi"), height = "600px") + %>% withSpinner(color="#0dc5c1"), + br(),br(), + span(em("Hint: By hovering over the column names you get a more detailed description of the columns.")), + br(),br(), + DT::dataTableOutput(ns("network_table")) + ) + ) + + + )) +} + + + +# module server function +networkMultiServer <- function(id) { + moduleServer( + id, + + function(input, output, session) { + + ### MULTI REGION NETWORK ------------------------------- + + # Display correct type of network in title + output$network_type <- renderText({ + v <- c("Vehicle Network", "Differential Network", "Treatment Network") + k <- c("baseline", "diff", "treatment") + l <- setNames(as.list(v), k) + + n <- l[[input$network_type]] + }) + + # Load data from multiple brain region + network_data <- reactive({ + + req(input$network_multireg) + + # Read data tables from all required regions + list_network <- list() + + for (reg in input$network_multireg){ + file_path <- paste0( + "tables/network/", + "04_singleRegion_", + reg, + "_filtered_", + input$network_type, + "Network.csv" + ) + # Read data + if (input$network_type == "diff") { + data <- + fread(file = file_path) %>% + dplyr::select(target, predictor, beta_mean.base, beta_mean.dex, z, p_adj) + #dplyr::select(target, predictor, z) %>% + #dplyr::rename_with(.fn = ~ reg, .cols = z) + + cols <- colnames(data)[3:6] + data <- data %>% + dplyr::rename_with(.fn = ~paste0(., ".", reg), .cols = all_of(cols) ) + } else { + data <- + fread(file = file_path) %>% + dplyr::select(target, predictor, beta_mean) + # cols <- colnames(data)[3] + data <- data %>% + dplyr::rename_with(.fn = ~ reg, .cols = beta_mean) + # dplyr::rename_with(.fn = ~paste0(., ".", reg), .cols = all_of(cols) ) + } + list_network[[reg]] <- data + } + + data_joined <- list_network %>% + purrr::reduce(full_join, by = c("target","predictor")) + + return(data_joined) + }) + + + + # DE genes of required brain regions + de_genes <- reactive({ + + list_de <- list() + # Read and subset DE genes for coloring + for (reg in input$network_multireg){ + df_de <- fread(paste0("tables/de/02_", + reg,"_deseq2_Dex_1_vs_0_lfcShrink.txt")) + na_indices <- which(is.na(df_de$padj)) + df_de$padj[na_indices] <- 1 + df_de <- df_de[df_de$padj <= 0.1,] + + df_de <- df_de %>% + dplyr::select(Ensembl_ID, padj) %>% + dplyr::rename_with(.fn = ~ reg, .cols = padj) + + list_de[[reg]] <- df_de + } + + de_joined <- list_de %>% + purrr::reduce(full_join, by = c("Ensembl_ID")) + + return(de_joined) + }) + + + # hub genes of required brain regions + hub_genes <- reactive({ + + list_hub <- list() + # Read and subset hub genes for coloring + for (reg in input$network_multireg) { + df_hub <- + fread( + paste0( + "tables/network/04_", + reg, + "_funcoup_differential", + "_nodebetweennessNorm_betacutoff0.01.csv" + ) + ) %>% + filter(nodebetweenness_norm >= 1) %>% + dplyr::select(ensembl_id, nodebetweenness_norm) %>% + dplyr::rename_with(.fn = ~ reg, .cols = nodebetweenness_norm) + + list_hub[[reg]] <- df_hub + } + + hub_joined <- list_hub %>% + purrr::reduce(full_join, by = c("ensembl_id")) + + return(hub_joined) + }) + + + # input genes + input_genes <- reactive({ + + if (isTruthy(input$network_gene)){ + genes <- input$network_gene + genes <- stringr::str_split(genes, ",\\s?")[[1]] + genes <- reformat_genes(genes) + } else if (isTruthy(input$file1)){ + genes <- read.csv(input$file1$datapath, header = FALSE)$V1 + genes <- reformat_genes(genes) + } + }) + + + # bring input genes to correct format + reformat_genes <- function(list_genes){ + if (!startsWith(list_genes[1], "ENSMUSG")){ + format_gene <- sapply(list_genes, stringr::str_to_title) + format_gene <- mapIds(org.Mm.eg.db, keys = format_gene, column = "ENSEMBL", keytype = "SYMBOL") + ids_na <- names(format_gene)[which(is.na(format_gene))] + #print(ids_na) + if (length(ids_na) > 0) { + showNotification(paste("No Ensembl ID found for following genes: ", + ids_na), type = "message", duration = 5) + } + format_gene <- format_gene[which(!is.na(format_gene))] + } else { + format_gene <- list_genes + } + return(format_gene) + } + + + # Network visualization + output$network_diff_multi <- renderVisNetwork({ + + req(input_genes()) + req(network_data()) + req(de_genes()) + req(hub_genes()) + + # Get Nodes and Edges + network <- network(network_data(), + de_genes(), + hub_genes(), + input_genes(), + input$network_type, + input$id_type, + input$vis_neighbours) + + visNetwork(network$nodes, network$edges) %>% + visEdges(arrows = "to") %>% + visOptions(highlightNearest = list(enabled = T, degree = 1, hover = T)) %>% + visLegend(addEdges = network$ledges, addNodes = network$lnodes, + useGroups = FALSE) %>% + visIgraphLayout() + }) + + + # Data table with network results + output$network_table <- DT::renderDataTable({ + + n_reg <- length(input$network_multireg) + + hover_info <- c("Gene ID of the dependent variable", + "Gene ID of the independent variable") + + if(input$network_type == "diff"){ + add_hover <- c() + + if ("beta" %in% input$tableContent){ + add_hover <- c(add_hover, "Mean beta weight for the predictor across multiple refits of the model in the specified brain region after vehicle treatment", + "Mean beta weight for the predictor across multiple refits of the model in the specified brain region after Dexamethasone treatment") + } + if ("z" %in% input$tableContent){ + add_hover <- c(add_hover, "Z-value calculated from mean beta values and their standard errors in vehicle and treatment condition in the specified brain region") + } + if("p" %in% input$tableContent){ + add_hover <- c(add_hover, "FDR corrected p-value for the z-value in the specified region") + } + + hover_info <- c(hover_info, + rep(add_hover, n_reg)) + + } else { + hover_info <- c(hover_info, + rep("Mean beta weight for the predictor across multiple refits of the model + in the specified brain region", n_reg)) + } + + network_table() %>% + datatable(filter = list(position = 'top'), + options = list( + headerCallback= JS(callHeaderCallback(hover_info)) + )) + }, server = TRUE) # FALSE to enable download of all pages with button + + + # Download of (filtered) network table + output$download2 <- downloadHandler( + filename = function() { + paste0("network_dexStimMouse_multiRegion.csv") + }, + content = function(file) { + indices <- input$network_table_rows_all + write.csv(network_table()[indices, ], file) + } + ) + + + # prepare network for visualization + network <- function(data, de_genes, hub_genes, input_genes, mode, id_type, neighbours){ + + # input genes + print(input_genes) + + # filter data + data <- data %>% + dplyr::rename(from = predictor, to = target) %>% + dplyr::filter(from != to) + data <- subset_network(data, input_genes, neighbours) + + # color palettes for edges and nodes + # edge palette --> assumes that data stores only target, predictor and z scores + # edge_colors <- brewer.pal(ncol(data)-1, "Blues")[2:(ncol(data)-1)] + + # Edges + if(mode == "baseline" | mode == "treatment"){ + # color palettes for edges and nodes + edge_colors <- brewer.pal(ncol(data)-1, "Blues")[2:(ncol(data)-1)] + + # count regions per diff edge that are not na + data$edge_reg <- apply(data[,3:ncol(data)], 1, function(x) names(which(!is.na(x))) ) + data$count_reg <- sapply(data$edge_reg, length) + data$title <- paste0("

Connection in: ", sapply(data$edge_reg, paste, collapse = ", "), + "

") + # assign color to edges according to number of regions + data$c <- sapply(data$count_reg, function(x) edge_colors[x]) + + #print(head(data)) + # df for edges + relations <- data.table(from=data$from, + to=data$to, + #beta = data$beta_mean, + color = data$c, + title = data$title) + # df for edge legend + ledges <- data.frame(color = edge_colors, + label = as.character(1:(ncol(data)-6)), + font.align = "top") + #print(ledges) + } else { + # remove columns + data <- data %>% dplyr::select(-contains("beta"), -contains("p_adj")) + + # color palettes for edges and nodes + # edge palette --> assumes that data stores only target, predictor and z scores + edge_colors <- brewer.pal(ncol(data)-1, "Blues")[2:(ncol(data)-1)] + + # count regions per diff edge that are not na + data$edge_reg <- apply(data[,3:ncol(data)], 1, function(x) names(which(!is.na(x))) ) + data$count_reg <- sapply(data$edge_reg, length) + data$title <- paste0("

Diff. co-expressed in: ", sapply(data$edge_reg, paste, collapse = ", "), + "

") + # assign color to edges according to number of regions + data$c <- sapply(data$count_reg, function(x) edge_colors[x]) + + #print(head(data)) + # df for edges + relations <- data.table(from=data$from, + to=data$to, + #z = data$z, + #p_adj = data$p_adj, + color = data$c, + title = data$title) + # df for edge legend + ledges <- data.frame(color = edge_colors, + label = as.character(1:(ncol(data)-6)), + font.align = "top") + #print(ledges) + } + #print(nrow(relations)) + + + # Nodes + # get unique nodes with correct id + nodes <- data.frame("id" = unique(union( + c(relations$from, relations$to), + input_genes + )), stringsAsFactors = FALSE) + if (id_type == "Gene Symbol"){ + #print(nodes$id) + nodes$label <- mapIds(org.Mm.eg.db, keys = nodes$id, + column = "SYMBOL", keytype = "ENSEMBL") + } else { + nodes$label <- nodes$id + } + + # count regions where gene is DE or/and hub + nodes_de <- left_join(x = nodes, y = de_genes, + by = c("id" = "Ensembl_ID")) + nodes_hub <- left_join(x = nodes, y = hub_genes, + by = c("id" = "ensembl_id")) + + nodes$de_reg <- apply(nodes_de[,3:ncol(nodes_de)], 1, + function(x) names(which(!is.na(x))) ) + nodes$de_count <- sapply(nodes$de_reg, length) + + nodes$hub_reg <- apply(nodes_hub[,3:ncol(nodes_hub)], 1, + function(x) names(which(!is.na(x))) ) + nodes$hub_count <- sapply(nodes$hub_reg, length) + + # set color according to DE/hub status + nodes$color <- rep("darkblue", nrow(nodes_de)) + nodes$color[nodes$de_count > 0] <- "orange" + nodes$color[nodes$hub_count > 0] <- "darkred" + nodes$color[nodes$de_count > 0 & nodes$hub_count > 0] <- "purple" + + #nodes$opacity <- nodes$de_count + nodes$hub_count + #nodes$opacity <- nodes$opacity/max(nodes$opacity) + + nodes$title <- paste0("

",nodes$label,"", + "
DE in regions: ", sapply(nodes$de_reg, paste, collapse = ", "), + "
Hub in regions: ", sapply(nodes$hub_reg, paste, collapse = ", "),"

") + + #print(head(nodes)) + + # df for node data frame + lnodes <- data.frame(label = c("DE", "hub", "DE & hub", "no DE & no hub"), + font.color = c("black", "white", "white", "white"), + shape = c( "ellipse"), + color = c("orange", "darkred", "purple", "darkblue"), + title = "Informations", id = 1:4) + + + return(list("edges" = relations, "nodes" = nodes, "ledges" = ledges, "lnodes" = lnodes)) + + } + + + + # prepare network table for visualization + network_table <- reactive({ + + req(network_data()) + req(input_genes()) + + # filter data + # remove self connections + data <- network_data() %>% + dplyr::rename(from = predictor, to = target) %>% + dplyr::filter(from != to) + data <- subset_network(data, input_genes(), input$vis_neighbours) %>% + dplyr::rename(predictor = from, target = to) + + # Network table + if(input$network_type == "baseline" | input$network_type == "treatment"){ + + # add beta to column name + cols <- colnames(data)[3:ncol(data)] + data <- data %>% + # dplyr::rename_with(.fn = ~ reg, .cols = beta_mean) + dplyr::rename_with(.fn = ~paste0("beta.", .), .cols = all_of(cols) ) + + } else { + + if (!("beta" %in% input$tableContent)){ + data <- data %>% dplyr::select(-contains("beta")) + } + if (!("z" %in% input$tableContent)){ + data <- data %>% dplyr::select(-contains("z")) + } + if(!("p" %in% input$tableContent)){ + data <- data %>% dplyr::select(-contains("p_adj")) + } + + } + + if (input$id_type == "Gene Symbol"){ + data$predictor <- mapIds(org.Mm.eg.db, keys = data$predictor, + column = "SYMBOL", keytype = "ENSEMBL") + data$target <- mapIds(org.Mm.eg.db, keys = data$target, + column = "SYMBOL", keytype = "ENSEMBL") + } + + data + + }) + + + } + + ) +} \ No newline at end of file diff --git a/R/network_singleRegion.R b/R/network_singleRegion.R new file mode 100644 index 0000000..73ae020 --- /dev/null +++ b/R/network_singleRegion.R @@ -0,0 +1,301 @@ +# module UI function +networkSingleUI <- function(id, label = "Network Single Region") { + ns <- NS(id) + + tagList(sidebarLayout( + sidebarPanel( + selectInput( + ns("network_region"), + "Choose a brain region:", + choices = c( + "Amygdala" = "AMY", + "Cerebellum" = "CER", + "Dorsal DG of Hippocampus" = "dDG", + "Dorsal CA1 of Hippocampus" = "dCA1", + "Prefrontal Cortex" = "PFC", + "PVN of Hypothalamus" = "PVN", + "Ventral DG of Hippocampus" = "vDG", + "Ventral CA1 of Hippocampus" = "vCA1" + ) + ), + searchInput( + inputId = ns("network_gene"), + label = "Enter comma-separated list of genes: ", + placeholder = "e.g. Nr3c1,Fkbp5", + btnSearch = icon("search", class = "fas fa-search", lib = "font-awesome"), + btnReset = icon("remove", class = "fas fa-remove", lib = "font-awesome"), + width = "100%" + ), + fileInput( + inputId = ns("file1"), + label = "Upload file with list of genes", + accept = "text/plain", + buttonLabel = "Upload", + placeholder = "Upload file with list of genes" + ), + radioButtons( + ns("network_type"), + label = "Select type of network to display:", + choiceNames = list("Vehicle", "Differential", "Treatment"), + choiceValues = list("baseline", "diff", "treatment"), + selected = "diff" + ), + radioButtons( + ns("vis_neighbours"), + label = "Include gene neighbourhood", + choiceNames = c("Yes", "No"), + choiceValues = c(TRUE, FALSE), + selected = FALSE + ), + radioButtons( + ns("id_type"), + label = "Select type of gene ID for plot and table:", + choices = list("Ensembl", "Gene Symbol"), + selected = "Gene Symbol" + ), + downloadButton(ns("download2"), "Download (filtered) table as csv"), + width = 3 + ), + mainPanel( + # visNetwork + fluidPage( + br(), + tags$style(".fa-project-diagram {color:#2980B9}"), + h3(p( + em("Network analysis of gene expression "), + icon("project-diagram", lib = "font-awesome"), + style = "color:black;text-align:center" + )), + hr(), + h4(span(textOutput( + ns("network_type") + ), style = "color:black;text-align:center")), + strong("Visualize"), + span(" a network in a "), + strong("single brain region"), + span( + " after vehicle treatment, Dexamethasone treatment or on a differential level. + Please choose the brain region you are interested in on the left panel. + You can enter your " + ), + strong("genes of interest"), + span( + " either as a comma-separated list or upload a file + with the respective gene IDs in gene symbol or ENSEMBL format. Furthermore, you have the + option to " + ), + strong("display neighbouring genes"), + span( + " and the respective connections in addition to + your genes of interest. Gene IDs can be displayed in the network as gene symbol or + ENSEMBL id. You have the option to " + ), + strong("filter and download"), + span(" the + network data table."), + br(),br(), + span(em("Hint: If no network shows up without including the neighbourhood, there are + probably no direct connections between the genes you entered. Please try to + include the gene neighbourhood."), br(), + em("However, if you include the gene neighbourhood, networks become quickly very large + and computationally expensive . We recommend to include the gene neighbourhood + only for up to 4 input genes.")), + br(), + br(), + visNetworkOutput(ns("network"), height = "600px") %>% + withSpinner(color = "#0dc5c1"), + br(),br(), + span(em("Hint: By hovering over the column names you get a more detailed description of the columns.")), + br(),br(), + DT::dataTableOutput(ns("network_table")) + ) + ) + )) + +} + + + +# module server function +networkSingleServer <- function(id) { + moduleServer( + id, + + function(input, output, session) { + + ### SINGLE REGION NETWORK ------------------------------- + + # Display correct type of network in title + output$network_type <- renderText({ + v <- c("Vehicle Network", "Differential Network", "Treatment Network") + k <- c("baseline", "diff", "treatment") + l <- setNames(as.list(v), k) + + n <- l[[input$network_type]] + }) + + # Load data from one brain region + network_df <- reactive({ + file_path <- paste0( + "tables/network/", + "04_singleRegion_", + input$network_region, + "_filtered_", + input$network_type, + "Network.csv" + ) + # Read data + if (input$network_type == "diff") { + data <- + fread(file = file_path) %>% + dplyr::select(target, predictor, beta_mean.base, beta_mean.dex, z, p_adj) + } else { + data <- + fread(file = file_path) %>% + dplyr::select(target, predictor, beta_mean) + } + + }) + + # DE genes of region + de_genes <- reactive({ + # Read and subset DE genes for coloring + df_de <- fread(paste0( "tables/de/02_", + input$network_region,"_deseq2_Dex_1_vs_0_lfcShrink.txt")) + na_indices <- which(is.na(df_de$padj)) + df_de$padj[na_indices] <- 1 + df_de <- df_de[df_de$padj <= 0.1,] + df_de$Ensembl_ID + }) + + # hub genes of region + hub_genes <- reactive({ + # Read and subset hub genes for coloring + df_hub <- fread(paste0("tables/network/04_", + input$network_region,"_funcoup_differential", + "_nodebetweennessNorm_betacutoff0.01.csv")) %>% + filter(nodebetweenness_norm >= 1) + df_hub$ensembl_id + }) + + + # input genes + input_genes <- reactive({ + + if (isTruthy(input$network_gene)){ + genes <- input$network_gene + genes <- stringr::str_split(genes, ",\\s?")[[1]] + genes <- reformat_genes(genes) + } else if (isTruthy(input$file1)){ + genes <- read.csv(input$file1$datapath, header = FALSE)$V1 + genes <- reformat_genes(genes) + } + #genes + }) + + + # bring input genes to correct format + reformat_genes <- function(list_genes){ + if (!startsWith(list_genes[1], "ENSMUSG")){ + format_gene <- sapply(list_genes, stringr::str_to_title) + format_gene <- mapIds(org.Mm.eg.db, keys = format_gene, column = "ENSEMBL", keytype = "SYMBOL") + } else { + format_gene <- list_genes + } + return(format_gene) + } + + + # Network visualization + output$network <- renderVisNetwork({ + + req(input_genes()) + req(network_df()) + req(de_genes()) + req(hub_genes()) + + # Get Nodes and Edges + network <- read_network(network_df(), + de_genes(), + hub_genes(), + input_genes(), + input$network_type, + input$id_type, + input$vis_neighbours) + + visNetwork(network$nodes, network$edges) %>% + visEdges(arrows = "to") %>% + visOptions(highlightNearest = TRUE) %>% + visLegend(addEdges = network$ledges, addNodes = network$lnodes, + useGroups = FALSE) + #visIgraphLayout() + }) + + # network table + network_data <- reactive({ + + req(input_genes()) + req(network_df()) + + # input genes + gene <- input_genes() + + # Get data and subset to neighbours of gene of interest + data <- network_df() %>% + dplyr::rename("from" = "target", "to" = "predictor") + data <- simplify_network_df(data) + data <- subset_network(data, gene, input$vis_neighbours) + + # ID type to gene symbol + if(input$id_type == "Gene Symbol"){ + data$from <- mapIds(org.Mm.eg.db, keys = data$from, + column = "SYMBOL", keytype = "ENSEMBL") + data$to <- mapIds(org.Mm.eg.db, keys = data$to, + column = "SYMBOL", keytype = "ENSEMBL") + } + + data %>% + dplyr::rename("target" = "from", "predictor" = "to") + }) + + # Data table with network results + output$network_table <- DT::renderDataTable({ + + if(input$network_type == "diff"){ + hover_info <- c("Gene ID of the dependent variable", + "Gene ID of the independent variable", + "Mean beta weight for the predictor across multiple refits of the model after vehicle treatment", + "Mean beta weight for the predictor across multiple refits of the model after Dexamethasone treatment", + "Z-value calculated from mean beta values and their standard errors in vehicle and treatment condition", + "FDR corrected p-value for the z-value") + + }else{ + hover_info <- c("Gene ID of the dependent variable", + "Gene ID of the independent variable", + "Mean beta weight for the predictor across multiple refits of the model") + } + + network_data() %>% + datatable(filter = list(position = 'top'), + options = list( + headerCallback= JS(callHeaderCallback(hover_info)) + )) + }, server = TRUE) # FALSE to enable download of all pages with button + + + # Download of (filtered) network table + output$download2 <- downloadHandler( + filename = function() { + paste0("network_dexStimMouse_", input$region, ".csv") + }, + content = function(file) { + indices <- input$network_table_rows_all + write.csv(network_data()[indices, ], file) + } + ) + + } + + ) + +} \ No newline at end of file diff --git a/R/upsetPlot_hub.R b/R/upsetPlot_hub.R new file mode 100644 index 0000000..1416680 --- /dev/null +++ b/R/upsetPlot_hub.R @@ -0,0 +1,446 @@ +# module UI function +comparisonHubGenesUI <- function(id, label = "Upset Plot Hub Genes"){ + + ns <- NS(id) + + #tabPanel("Comparison Brain Regions", + tagList(sidebarLayout( + sidebarPanel( + pickerInput( + inputId = ns("myPicker"), + label = "Select brain regions", + choices = c("Amygdala" = "AMY", + "Cerebellum" = "CER", + "Dorsal DG of Hippocampus" = "dDG", + "Dorsal CA1 of Hippocampus" = "dCA1", + "Prefrontal Cortex" = "PFC", + "PVN of Hypothalamus" = "PVN", + "Ventral DG of Hippocampus" = "vDG", + "Ventral CA1 of Hippocampus" = "vCA1"), + options = list( + `actions-box` = TRUE, + size = 8, + `selected-text-format` = "count > 3" + ), + multiple = TRUE + ), + radioButtons( + ns("network_type"), + label = "Select type of network to display:", + choiceNames = list("Baseline", "Differential", "Treatment"), + choiceValues = list("baseline", "differential", "treatment"), + selected = "differential" + ), + sliderInput( + inputId = ns("normBetween"), + label = "Choose a threshold for the normalized nodebetweenness:", + min = 0.5, + max = 1.5, + step = 0.05, + value = 1.0 + ), + sliderInput( + ns("nintersects"), + label = "Number of intersections", + min = 2, + max = 40, + step = 1, + value = 20 + ), + downloadButton(ns("download"),"Download (filtered) table as csv"), + width = 3 + ), + mainPanel( + br(), + tags$style(".fa-brain {color:#2980B9}"), + h3(p(em("Comparsion of hub genes across brain regions "),icon("brain",lib = "font-awesome"),style="color:black;text-align:center")), + hr(), + + tags$style(type="text/css", + ".shiny-output-error { visibility: hidden; }", + ".shiny-output-error:before { visibility: hidden; }" + ), + strong("Compare"),span(" hub genes either at baseline, after Dexamethasone treatment or at the differential level"), + strong(" between different brain regions"), span(". Please choose the brain regions you are interested + in on the left panel. You can also select a different threshold for the normalized + nodebetweenness. The table on the bottom displays all genes in the dataset with the brain regions + they are hub genes in."), + br(),br(), + plotlyOutput(ns("plotly_upset"), height = "600px"), + DT::dataTableOutput(ns("upset_table")) + #h3("Intersection plots") + )) + ) +} + +# module server function +comparisonHubGenesServer <- function(id) { + moduleServer( + id, + + function(input, output, session) { + + # Read data required for upset plot + read_upset_data <- function(basepath, regions, between_cutoff = 1.0){ + + # Read DE tables from all required regions + list_genes_sig <- list() + + for (reg in regions){ + f <- Sys.glob(file.path(paste0("tables/network/"), + paste0("*_", reg, "_funcoup_", input$network_type, + "_nodebetweennessNorm_betacutoff0.01.csv"))) + res <- fread(file=f, sep=",") + #na_indices <- which(is.na(res$padj)) + #res$padj[na_indices] <- 1 + res_sig <- res[res$nodebetweenness_norm >= between_cutoff,] + list_genes_sig[[reg]] <- res_sig$ensembl_id + } + + return(list_genes_sig) + } + + # Read data required for data table below upset plot + upset_datatable <- function(basepath, regions, between_cutoff = 1.0){ + + # READ DATA + + # Read DE tables from all required regions + list_reg_sig <- list() + + for (reg in regions){ + f <- Sys.glob(file.path(paste0("tables/network/"), + paste0("*_", reg, "_funcoup_", input$network_type, + "_nodebetweennessNorm_betacutoff0.01.csv"))) + res <- fread(file=f, sep=",") + #na_indices <- which(is.na(res$padj)) + #res$padj[na_indices] <- 1 + res_sig <- res[res$nodebetweenness_norm >= between_cutoff,] + list_reg_sig[[reg]] <- res_sig + } + + df <- bind_rows(list_reg_sig, .id="region") + df <- df[,c("ensembl_id", "region")] %>% + group_by(ensembl_id) %>% + dplyr::summarise(region = list(region)) %>% + dplyr::select(ensembl_id, region) + df$gene_symbol <- mapIds(org.Mm.eg.db, keys = df$ensembl_id, + column = "SYMBOL", keytype = "ENSEMBL") + df <- dplyr::select(df, gene_symbol, ensembl_id, region) + + return(df) + } + + + # reactive table + upset_data <- reactive({ + + req(input$myPicker) + data <- upset_datatable(basepath, input$myPicker, input$normBetween) %>% + mutate(region = sapply(region, toString)) + data + + }) + + + # Data table below upset plot + output$upset_table <- DT::renderDataTable({ + + upset_data() %>% + datatable(filter = list(position = 'top')) + }, server = TRUE) + + # Download of (filtered) network table + output$download <- downloadHandler( + filename = "hubGenes.csv", + content = function(file) { + indices <- input$upset_table_rows_all + print(upset_data()) + write.csv(upset_data()[indices, ], file) + } + ) + + + # UPSET PLOT + + # TODO: mark here with a comment (SPDX) that following code snippet is + # under AGPL license https://github.com/pinin4fjords/shinyngs/blob/master/R/upset.R + # --> whole repo needs to be AGPL but everything that is my own can be MIT + + # Accessor for user-selected sets + + getSelectedSetNames <- reactive({ + req(input$myPicker) + input$myPicker + }) + + + # Get the sets we're going to use based on nsets + + getSets <- reactive({ + selected_set_names <- getSelectedSetNames() + req(length(selected_set_names) > 0) + + selected_sets <- read_upset_data(basepath, selected_set_names, input$normBetween) + }) + + + # Accessor for the nintersections parameter + + getNintersections <- reactive({ + validate(need(!is.null(input$nintersects), "Waiting for nintersects")) + input$nintersects + }) + + # Calculate intersections between sets + + calculateIntersections <- reactive({ + selected_sets <- getSets() + + withProgress(message = "Calculating set intersections", value = 0, { + sets <- getSets() + nsets <- length(sets) + + # Get all possible combinations of sets + + combinations <- function(items, pick) { + x <- combn(items, pick) + lapply(seq_len(ncol(x)), function(i) + x[, i]) + } + + combos <- lapply(1:nsets, function(x) { + combinations(1:length(selected_sets), x) + }) + combos <- lapply(1:nsets, function(x) { + combinations(1:nsets, x) + }) + + # Calculate the intersections of all these combinations + + withProgress(message = "Running intersect()", value = 0, { + intersects <- lapply(combos, function(combonos) { + lapply(combonos, function(combo) { + Reduce(intersect, selected_sets[combo]) + }) + }) + }) + + # For UpSet-ness, membership of higher-order intersections takes priority + # Otherwise just return the number of entries in each intersection + + intersects <- lapply(1:length(intersects), function(i) { + intersectno <- intersects[[i]] + if (i != length(intersects)){ + members_in_higher_levels <- + unlist(intersects[(i + 1):length(intersects)]) + } else { + members_in_higher_levels <- NULL + } + lapply(intersectno, function(intersect) { + length(setdiff(intersect, members_in_higher_levels)) + }) + }) + + combos <- unlist(combos, recursive = FALSE) + intersects <- unlist(intersects) + + + # Sort by intersect size + + combos <- combos[order(intersects, decreasing = TRUE)] + intersects <- + intersects[order(intersects, decreasing = TRUE)] + list(combinations = combos, intersections = intersects) + + }) + + }) + + output$plotly_upset <- renderPlotly({ + grid <- upsetGrid() + set_size_chart <- upsetSetSizeBarChart() + intersect_size_chart <- upsetIntersectSizeBarChart() + + # add axis labels + intersect_size_chart <- + intersect_size_chart %>% layout(yaxis = list(title = "Intersections size")) + + set_size_chart <- + set_size_chart %>% layout(xaxis = list(title = "Set size"), + yaxis = list(title = "Brain region")) + + # subplots + s1 <- + subplot( + plotly_empty(type = "scatter", mode = "markers"), + set_size_chart, + nrows = 2, + # widths = c(0.6, 0.4), + titleX = TRUE, + titleY = TRUE + ) %>% layout(showlegend = FALSE) + s2 <- + subplot(intersect_size_chart, + grid, + nrows = 2, + shareX = TRUE, titleY = TRUE) %>% layout(showlegend = FALSE) + + subplot(s1, s2, widths = c(0.25, 0.75), + shareX=F, + shareY=F, + titleX=T, + titleY=T) + + + }) + + # Add some line returns to contrast names + + getSetNames <- reactive({ + selected_sets <- getSets() + names(selected_sets) + }) + + # Make the grid of points indicating set membership in intersections + + upsetGrid <- reactive({ + selected_sets <- getSets() + ints <- calculateIntersections() + + intersects <- ints$intersections + combos <- ints$combinations + + # Reduce the maximum number of intersections if we don't have that many + nintersections <- getNintersections() + nintersections <- min(nintersections, length(combos)) + + # Fetch the number of sets + nsets <- length(getSelectedSetNames()) + setnames <- getSelectedSetNames() + + lines <- + data.table::rbindlist(lapply(1:nintersections, function(combono) { + data.frame( + combo = combono, + x = rep(combono, max(2, length(combos[[combono]]))), + y = (nsets - combos[[combono]]) + 1, + name = setnames[combos[[combono]]] + ) + })) + + plot_ly( + type = "scatter", + mode = "markers", + marker = list(color = "lightgrey", size = 8), + customdata = "grid", + source = "upset_hub" + ) %>% add_trace( + type = "scatter", + x = rep(1:nintersections, + length(selected_sets)), + y = unlist(lapply(1:length(selected_sets), function(x) + rep(x - 0.5, nintersections))), + hoverinfo = "none" + ) %>% add_trace( + type = "scatter", + data = group_by(lines, combo), + mode = "lines+markers", + x = lines$x, + y = lines$y - 0.5, + line = list(color = "#2980B9", width = 3), + marker = list(color = "#2980B9", + size = 10), + hoverinfo = "text", + text = ~ name + ) %>% layout( + xaxis = list( + showticklabels = FALSE, + showgrid = FALSE, + zeroline = FALSE + ), + yaxis = list( + showticklabels = FALSE, + showgrid = TRUE, + range = c(0, nsets), + zeroline = FALSE, + range = 1:nsets + ), + margin = list(t = 0, b = 40) + ) + + }) + + # Make the bar chart illustrating set sizes + + upsetSetSizeBarChart <- reactive({ + setnames <- getSelectedSetNames() + selected_sets <- getSets() + + plot_ly( + x = unlist(lapply(selected_sets, length)), + y = setnames, + type = "bar", + orientation = "h", + marker = list(color = "#2c3e50"), + customdata = "setsize", + source = "upset_hub" + ) %>% layout( + bargap = 0.4, + yaxis = list( + categoryarray = rev(setnames), + categoryorder = "array" + ) + ) + }) + + # Make the bar chart illustrating intersect size + + upsetIntersectSizeBarChart <- reactive({ + ints <- calculateIntersections() + intersects <- ints$intersections + combos <- ints$combinations + nintersections <- getNintersections() + + p <- + plot_ly(showlegend = FALSE, + customdata = "intersectsize", + source = "upset_hub") %>% add_trace( + x = 1:nintersections, + y = unlist(intersects[1:nintersections]), + type = "bar", + marker = list(color = "#2c3e50", + hoverinfo = "none") + ) + + bar_numbers <- TRUE + + if (bar_numbers) { + p <- + p %>% add_trace( + type = "scatter", + mode = "text", + x = 1:nintersections, + y = unlist(intersects[1:nintersections]) + (max(intersects) * 0.05), + text = unlist(intersects[1:nintersections]), + textfont = list(color = "black") + ) + } + + p + }) + + observe({ + click <- event_data("plotly_click", source = "upset_hub") + req(click) + print(click) + }) + + return(NULL) + } + + + + + ) +} \ No newline at end of file diff --git a/R/upset_plot.R b/R/upset_plot.R new file mode 100644 index 0000000..9e4cf78 --- /dev/null +++ b/R/upset_plot.R @@ -0,0 +1,469 @@ +# module UI function +upsetPlotUI <- function(id, label = "Upset Plot"){ + + ns <- NS(id) + + #tabPanel("Comparison Brain Regions", + tagList(sidebarLayout( + sidebarPanel( + pickerInput( + inputId = ns("myPicker"), + label = "Select brain regions", + choices = c("Amygdala" = "AMY", + "Cerebellum" = "CER", + "Dorsal DG of Hippocampus" = "dDG", + "Dorsal CA1 of Hippocampus" = "dCA1", + "Prefrontal Cortex" = "PFC", + "PVN of Hypothalamus" = "PVN", + "Ventral DG of Hippocampus" = "vDG", + "Ventral CA1 of Hippocampus" = "vCA1"), + options = list( + `actions-box` = TRUE, + size = 8, + `selected-text-format` = "count > 3" + ), + multiple = TRUE + ), + sliderTextInput( + inputId = ns("padj_upset"), + label = "Choose a threshold for FDR p-value:", + choices = c(0.01, 0.05, 0.1), + selected = 0.1, + grid = TRUE + ), + sliderTextInput( + inputId = ns("FC_upset"), + label = "Choose a threshold for logFC:", + choices = c(0.0, 0.5, 1.0, 2.0), + selected = 0.0, + grid = TRUE + ), + sliderInput( + ns("nintersects"), + label = "Number of intersections", + min = 2, + max = 40, + step = 1, + value = 20 + ), + downloadButton(ns("download"),"Download (filtered) table as csv"), + width = 3 + ), + mainPanel( + br(), + tags$style(".fa-brain {color:#2980B9}"), + h3(p(em("Comparsion of DE genes across brain regions "),icon("brain",lib = "font-awesome"),style="color:black;text-align:center")), + hr(), + tags$style(type="text/css", + ".shiny-output-error { visibility: hidden; }", + ".shiny-output-error:before { visibility: hidden; }" + ), + strong("Compare"),span(" the differentially expressed (DE) genes after Dexamethasone treatment"), + strong(" between different brain regions"), span(". Please choose the brain regions you are interested + in on the left panel. You can also select a different FDR p-value threshold and a different fold + change cutoff. The table on the bottom displays all genes in the dataset with the brain regions + they are differentially expressed in."), + br(),br(), + # plotlyOutput("upset_de") + #plotOutput("upset_de"), + plotlyOutput(ns("plotly_upset"), height = "600px"), + br(),br(), + span(em("Hint: By hovering over the column names you get a more detailed description of the columns.")), + br(),br(), + DT::dataTableOutput(ns("upset_table")) + #h3("Intersection plots") + )) + ) +} + +# module server function +upsetPlotServer <- function(id) { + moduleServer( + id, + + function(input, output, session) { + + # Read data required for upset plot + read_upset_data <- function(basepath, regions, padj_cutoff = 0.1, fc_cutoff = 0){ + + # Read DE tables from all required regions + list_genes_sig <- list() + + for (reg in regions){ + res <- fread(file=paste0("tables/de/02_", reg, + "_deseq2_Dex_1_vs_0_lfcShrink.txt"), + sep="\t") + na_indices <- which(is.na(res$padj)) + res$padj[na_indices] <- 1 + res_sig <- res[res$padj <= padj_cutoff,] + res_sig <- res_sig[abs(res_sig$log2FoldChange) >= fc_cutoff,] + list_genes_sig[[reg]] <- res_sig$Ensembl_ID + } + + return(list_genes_sig) + } + + # Read data required for data table below upset plot + upset_datatable <- function(basepath, regions, padj_cutoff = 0.1, fc_cutoff = 0){ + + # READ DATA + + # Read DE tables from all required regions + list_reg_sig <- list() + + for (reg in regions){ + res <- fread(file=paste0("tables/de/02_", reg, + "_deseq2_Dex_1_vs_0_lfcShrink.txt"), + sep="\t") + na_indices <- which(is.na(res$padj)) + res$padj[na_indices] <- 1 + res_sig <- res[res$padj <= padj_cutoff,] + res_sig <- res_sig[abs(res_sig$log2FoldChange) >= fc_cutoff,] + list_reg_sig[[reg]] <- res_sig + } + + df <- bind_rows(list_reg_sig, .id="region") + df <- df[,c("Ensembl_ID", "region")] %>% + group_by(Ensembl_ID) %>% + dplyr::summarise(region = list(region)) %>% + dplyr::select(Ensembl_ID, region) %>% + dplyr::distinct(Ensembl_ID, .keep_all = TRUE) + df$Gene_Symbol <- mapIds(org.Mm.eg.db, keys = df$Ensembl_ID, + column = "SYMBOL", keytype = "ENSEMBL") + df <- dplyr::select(df, Gene_Symbol, Ensembl_ID, region) + + # df <- df[,c("Ensembl_ID", "Gene_Symbol", "region", "padj", "log2FoldChange")] %>% + # group_by(Ensembl_ID, Gene_Symbol) %>% + # dplyr::summarise(region = list(region), + # p_adjusted = list(padj), + # log2FC = list(log2FoldChange)) %>% + # dplyr::select(Gene_Symbol, Ensembl_ID, region, p_adjusted, log2FC) + + return(df) + } + + # reactive table + upset_data <- reactive({ + + req(input$myPicker) + data <- upset_datatable(basepath, input$myPicker, input$padj_upset, input$FC_upset) %>% + mutate(region = sapply(region, toString)) + data + + }) + + + # display information when hovering column names + hover_info <- c("Gene ID in gene symbol format", + "Gene ID in Ensembl format", + "List of brain regions this gene is differentially expressed in") + + # Data table below upset plot + output$upset_table <- DT::renderDataTable({ + + hover_info <- c("Gene ID in gene symbol format", + "Gene ID in Ensembl format", + "List of brain regions this gene is differentially expressed in") + + upset_data() %>% + datatable(filter = list(position = 'top'), + options = list( + headerCallback= JS(callHeaderCallback(hover_info)) + )) + }, server = TRUE) + + # Download of (filtered) network table + output$download <- downloadHandler( + filename = "DEGenes.csv", + content = function(file) { + indices <- input$upset_table_rows_all + print(upset_data()) + write.csv(upset_data()[indices, ], file) + } + ) + + + + # UPSET PLOT + + # TODO: mark here with a comment (SPDX) that following code snippet is + # under AGPL license https://github.com/pinin4fjords/shinyngs/blob/master/R/upset.R + # --> whole repo needs to be AGPL but everything that is my own can be MIT + + # Accessor for user-selected sets + + getSelectedSetNames <- reactive({ + req(input$myPicker) + input$myPicker + }) + + + # Get the sets we're going to use based on nsets + + getSets <- reactive({ + selected_set_names <- getSelectedSetNames() + req(length(selected_set_names) > 0) + + selected_sets <- read_upset_data(basepath, selected_set_names, input$padj_upset, input$FC_upset) + }) + + + # Accessor for the nintersections parameter + + getNintersections <- reactive({ + validate(need(!is.null(input$nintersects), "Waiting for nintersects")) + input$nintersects + }) + + # Calculate intersections between sets + + calculateIntersections <- reactive({ + selected_sets <- getSets() + + withProgress(message = "Calculating set intersections", value = 0, { + sets <- getSets() + nsets <- length(sets) + + # Get all possible combinations of sets + + combinations <- function(items, pick) { + x <- combn(items, pick) + lapply(seq_len(ncol(x)), function(i) + x[, i]) + } + + combos <- lapply(1:nsets, function(x) { + combinations(1:length(selected_sets), x) + }) + combos <- lapply(1:nsets, function(x) { + combinations(1:nsets, x) + }) + + # Calculate the intersections of all these combinations + + withProgress(message = "Running intersect()", value = 0, { + intersects <- lapply(combos, function(combonos) { + lapply(combonos, function(combo) { + Reduce(intersect, selected_sets[combo]) + }) + }) + }) + + # For UpSet-ness, membership of higher-order intersections takes priority + # Otherwise just return the number of entries in each intersection + + intersects <- lapply(1:length(intersects), function(i) { + intersectno <- intersects[[i]] + if (i != length(intersects)){ + members_in_higher_levels <- + unlist(intersects[(i + 1):length(intersects)]) + } else { + members_in_higher_levels <- NULL + } + lapply(intersectno, function(intersect) { + length(setdiff(intersect, members_in_higher_levels)) + }) + }) + + combos <- unlist(combos, recursive = FALSE) + intersects <- unlist(intersects) + + + # Sort by intersect size + + combos <- combos[order(intersects, decreasing = TRUE)] + intersects <- + intersects[order(intersects, decreasing = TRUE)] + list(combinations = combos, intersections = intersects) + + }) + + }) + + output$plotly_upset <- renderPlotly({ + grid <- upsetGrid() + set_size_chart <- upsetSetSizeBarChart() + intersect_size_chart <- upsetIntersectSizeBarChart() + + # add axis labels + intersect_size_chart <- + intersect_size_chart %>% layout(yaxis = list(title = "Intersections size")) + + set_size_chart <- + set_size_chart %>% layout(xaxis = list(title = "Set size"), + yaxis = list(title = "Brain region")) + + # subplots + s1 <- + subplot( + plotly_empty(type = "scatter", mode = "markers"), + set_size_chart, + nrows = 2, + # widths = c(0.6, 0.4), + titleX = TRUE, + titleY = TRUE + ) %>% layout(showlegend = FALSE) + s2 <- + subplot(intersect_size_chart, + grid, + nrows = 2, + shareX = TRUE, titleY = TRUE) %>% layout(showlegend = FALSE) + + subplot(s1, s2, widths = c(0.25, 0.75), + shareX=F, + shareY=F, + titleX=T, + titleY=T) + + + }) + + # Add some line returns to contrast names + + getSetNames <- reactive({ + selected_sets <- getSets() + names(selected_sets) + }) + + # Make the grid of points indicating set membership in intersections + + upsetGrid <- reactive({ + selected_sets <- getSets() + ints <- calculateIntersections() + + intersects <- ints$intersections + combos <- ints$combinations + + # Reduce the maximum number of intersections if we don't have that many + nintersections <- getNintersections() + nintersections <- min(nintersections, length(combos)) + + # Fetch the number of sets + nsets <- length(getSelectedSetNames()) + setnames <- getSelectedSetNames() + + lines <- + data.table::rbindlist(lapply(1:nintersections, function(combono) { + data.frame( + combo = combono, + x = rep(combono, max(2, length(combos[[combono]]))), + y = (nsets - combos[[combono]]) + 1, + name = setnames[combos[[combono]]] + ) + })) + + plot_ly( + type = "scatter", + mode = "markers", + marker = list(color = "lightgrey", size = 8), + customdata = "grid", + source = "upset_de" + ) %>% add_trace( + type = "scatter", + x = rep(1:nintersections, + length(selected_sets)), + y = unlist(lapply(1:length(selected_sets), function(x) + rep(x - 0.5, nintersections))), + hoverinfo = "none" + ) %>% add_trace( + type = "scatter", + data = group_by(lines, combo), + mode = "lines+markers", + x = lines$x, + y = lines$y - 0.5, + line = list(color = "#2980B9", width = 3), + marker = list(color = "#2980B9", + size = 10), + hoverinfo = "text", + text = ~ name + ) %>% layout( + xaxis = list( + showticklabels = FALSE, + showgrid = FALSE, + zeroline = FALSE + ), + yaxis = list( + showticklabels = FALSE, + showgrid = TRUE, + range = c(0, nsets), + zeroline = FALSE, + range = 1:nsets + ), + margin = list(t = 0, b = 40) + ) + + }) + + # Make the bar chart illustrating set sizes + + upsetSetSizeBarChart <- reactive({ + setnames <- getSelectedSetNames() + selected_sets <- getSets() + + plot_ly( + x = unlist(lapply(selected_sets, length)), + y = setnames, + type = "bar", + orientation = "h", + marker = list(color = "#2c3e50"), + customdata = "setsize", + source = "upset_de" + ) %>% layout( + bargap = 0.4, + yaxis = list( + categoryarray = rev(setnames), + categoryorder = "array" + ) + ) + }) + + # Make the bar chart illustrating intersect size + + upsetIntersectSizeBarChart <- reactive({ + ints <- calculateIntersections() + intersects <- ints$intersections + combos <- ints$combinations + nintersections <- getNintersections() + + p <- + plot_ly(showlegend = FALSE, + customdata = "intersectsize", + source = "upset_de") %>% add_trace( + x = 1:nintersections, + y = unlist(intersects[1:nintersections]), + type = "bar", + marker = list(color = "#2c3e50", + hoverinfo = "none") + ) + + bar_numbers <- TRUE + + if (bar_numbers) { + p <- + p %>% add_trace( + type = "scatter", + mode = "text", + x = 1:nintersections, + y = unlist(intersects[1:nintersections]) + (max(intersects) * 0.05), + text = unlist(intersects[1:nintersections]), + textfont = list(color = "black") + ) + } + + p + }) + + observe({ + click <- event_data("plotly_click", source = "upset_de") + req(click) + print(click) + }) + + return(NULL) + } + + + + + ) +} \ No newline at end of file diff --git a/app.R b/app.R new file mode 100644 index 0000000..84b346c --- /dev/null +++ b/app.R @@ -0,0 +1,2 @@ + +shinyApp(ui,server) \ No newline at end of file diff --git a/docker/Dockerfile b/docker/Dockerfile new file mode 100644 index 0000000..33bf732 --- /dev/null +++ b/docker/Dockerfile @@ -0,0 +1,55 @@ +FROM openanalytics/r-base + +MAINTAINER Nathalie Gerstner "nathalie_gerstner@psych.mpg.de" + +# system libraries of general use +RUN apt-get update && apt-get install -y \ + sudo \ + pandoc \ + pandoc-citeproc \ + libcurl4-gnutls-dev \ + libcairo2-dev \ + libxt-dev \ + libssl-dev \ + libssh2-1-dev \ + libssl1.1 \ + texlive-latex-recommended \ + texlive-latex-extra \ + language-pack-de \ + && rm -rf /var/lib/apt/lists/* + +# system library dependency for the app +# RUN apt-get update && apt-get install -y \ +# libmpfr-dev \ +# && rm -rf /var/lib/apt/lists/* + +# Docker inheritance +FROM bioconductor/bioconductor_docker:devel + +RUN apt-get update + +RUN R -e 'BiocManager::install(ask = F)' && R -e 'BiocManager::install(c("org.Mm.eg.db", \ + ask = F))' + +# basic shiny functionality +RUN R -e "install.packages(c('shiny', 'markdown', 'shinythemes', 'shinyWidgets', 'shinycssloaders'), \ + repos = 'https://cloud.r-project.org/', \ + Ncpus = parallel::detectCores())" + +# install dependencies of the diffbrainnet app +RUN R -e "install.packages(c('tidyverse', 'plotly', 'DT', 'data.table', 'scales', 'visNetwork', 'igraph'), \ + repos='https://cloud.r-project.org/', \ + Ncpus = parallel::detectCores())" + + +# copy the app to the image +RUN mkdir /root/diffbrainnet +COPY diffbrainnet /root/diffbrainnet + +#COPY Rprofile.site /usr/lib/R/etc/ + +EXPOSE 3838 + +CMD ["R", "-e", "shiny::runApp('/root/diffbrainnet', host = '0.0.0.0', port = 3838)"] + +#CMD ["R", "-e", "shiny::runApp('/root/diffbrainnet')"] diff --git a/docker/documentation.md b/docker/documentation.md new file mode 100644 index 0000000..950539c --- /dev/null +++ b/docker/documentation.md @@ -0,0 +1,9 @@ +## Documentation: Dockerizing the app + +- start Docker engine on the computer +- check existing images with `docker images` -> the image, your image is based on, needs to be there +- from project folder (one folder above diffbrainnet folder) where the *Dockerfile* needs to be stored: `docker build -t ngerst/diffbrainnet:v6 .` (or the version one wants to build) +- run container to try it out with `docker run -d --rm -p 3838:3838 ngerst/diffbrainnet:v6`` +- in your browser open *localhost:3838* to see the app +- run `docker tag ngerst/diffbrainnet:v6 ngerst/diffbrainnet:latest` to set the latest tag to the current image +- run `docker push ngerst/diffbrainnet:v6` and 'docker push ngerst/diffbrainnet:latest` to publish image on *Docker Hub* \ No newline at end of file diff --git a/global.R b/global.R new file mode 100644 index 0000000..9f38464 --- /dev/null +++ b/global.R @@ -0,0 +1,34 @@ +library(shiny) +library(shinythemes) +library(shinyWidgets) # custom input widgets for shiny +library(shinycssloaders) # add loading animations +library(markdown) +library(plotly) +library(DT) # visualization of data tables +library(data.table) +library(stringr) +library(dplyr) +library(ggplot2) +library(visNetwork) # visualization of networks +library(org.Mm.eg.db) +library(scales) # graphical scales (map data to colors etc.) +library(RColorBrewer) + +source("utilities_network.R") + + +# function necessary to display column name information when hovering +callHeaderCallback <- function(hover_info){ + + headerCallback <- c( + "function(thead, data, start, end, display){", + sprintf(" var tooltips = [%s];", toString(paste0("'", hover_info, "'"))), + " for(var i = 1; i <= tooltips.length; i++){", + " $('th:eq('+i+')',thead).attr('title', tooltips[i-1]);", + " }", + "}" + ) + + return(headerCallback) + +} \ No newline at end of file diff --git a/server.R b/server.R new file mode 100644 index 0000000..7a5d278 --- /dev/null +++ b/server.R @@ -0,0 +1,379 @@ +# Define server logic required to generate and plot +server <- shinyServer(function(input, output) { + + ### TERMS OF USE & IMPRINT ---------------------------------- + + observeEvent(input$termsOfUse, { + showModal(modalDialog( + h4("Terms of Use, v1.0 (03/2022)"), + tags$ol( + tags$p("By using the software application DiffBrainNet (hereinafter “DiffBrainNet”), + you acknowledge that you have read these terms and conditions, understand + them, and agree to be bound by them. If you do not agree with these + terms and conditions, you must not use DiffBrainNet. Please read the + following Terms of Use and any accompanying documentation carefully + before you use the application DiffBrainNet."), + tags$li( + tags$h5("Contents and availability"), + tags$p( + "DiffBrainNet consists of the entire web application, i.e., the + texts, figures, and data shown on this page, as well as the programming + source code with which the visualizations and results are generated. " + ), + tags$p( + "The use of DiffBrainNet is licensed to You by the Max-Planck-Gesellschaft + zur Förderung der Wissenschaften e.V., here represented by the managing + director Elisabeth Binder, acting for the Max Planck Institute of Psychiatry, + with offices at Kraepelinstr. 2-10, 80804 Munich, Germany (hereinafter: + Licensor) subject to these Terms of Use." + ), + tags$p( + "All figures, texts, and data shown in DiffBrainNet are licensed under + the Creative Commons Attribution license (", + a(href="https://creativecommons.org/licenses/by/4.0/", + "https://creativecommons.org/licenses/by/4.0/", target="_blank"),"). + The use of the DiffBrainNet software is licensed under the GNU Affero + General Public License v3.0 (", + a(href="http://www.gnu.org./licenses/", "http://www.gnu.org./licenses/", target="_blank"), + "). " + ), + tags$p( + "The Licensor does not bear any responsibility for i) third-party contents, + ii) for the use of DiffBrainNet different from the Purpose set out below, + iii) for the accessibility via external links, iv) nor for its availability. + Furthermore, the Licensor waives any responsibility for contents that + may be illegal or violate common decency." + ) + ), + tags$li( + tags$h5("Purpose"), + tags$p( + "The Licensor has built DiffBrainNet in order to allow the scientific + community to explore the transcriptional landscape of 8 mouse brain + regions before and after treatment with dexamethasone. " + ), + tags$p( + "The use and output of DiffBrainNet can only be complementary to and + neither replace informed judgement, nor can it be used as a stand-alone + recommendation without further validation. " + ) + ), + tags$li( + tags$h5("Rules of Use"), + tags$p( + "You accept to use DiffBrainNet non-commercially and for the aforementioned + Purpose. You are not authorized to make the contents of the use of DiffBrainNet + available to third parties. In duly substantiated exceptional cases, + the Licensor may restrict Your right to use DiffBrainNet." + ) + ), + tags$li( + tags$h5("Proprietary rights to the contents of DiffBrainNet"), + tags$p( + "You agree to respect copyrights, rights to names, trademarks and other + intellectual property rights, if any, of the Licensor and of third parties + when making use of DiffBrainNet. By enabling access to the use of DiffBrainNet, + the Licensor does not grant any licence nor any other right of use." + ) + ), + tags$li( + tags$h5("Improper use of DiffBrainNet"), + tags$p( + "You agree to refrain from any improper use of the use of DiffBrainNet; + in particular, no security precautions must be circumvented. DiffBrainNet + may not be used to create fake, libellous, misleading, or defamatory + content of any kind." + ), + tags$p( + "Furthermore, no facilities may be used nor may any applications be + run that could lead to a damage or a performance failure of any of the + facilities from which the use of DiffBrainNet is provided, in particular + through changes in the physical or logical structure of the servers + or its network or of any other network. No commercial, systematic use + of DiffBrainNet is permitted without the Licensor's consent." + ) + ), + tags$li( + tags$h5("Limitation of Liability (for Civil Law Countries)"), + tags$p( + "Without prejudice to Licensor´s responsibility for tort or for violation + of mandatory laws of the Federal Republic of Germany, the Licensor shall + not be liable for any damages You incur when using DiffBrainNet." + ) + ), + tags$li( + tags$h5( + "Representation, Warranty, Limitation of Liability, Indemnification + (for Common Law Countries)" + ), + tags$ol( + tags$li( + tags$h5("Representation and Warranty"), + tags$p( + "LICENSOR REPRESENTS THAT LICENSOR HAS ALL RIGHTS REQUIRED TO MAKE + AVAILABLE AND DISTRIBUTE THE MATERIALS. EXCEPT FOR SUCH REPRESENTATION, + THE MATERIALS ARE PROVIDED “AS IS” AND “AS AVAILABLE” AND WITHOUT + WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED + TO, NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR + PURPOSE, AND ANY WARRANTIES IMPLIED BY ANY COURSE OF PERFORMANCE OR + USAGE OF TRADE, ALL OF WHICH ARE EXPRESSLY DISCLAIMED." + ), + tags$p( + "WITHOUT LIMITING THE FOREGOING, LICENSOR DOES NOT WARRANT THAT: (A) + THE MATERIALS ARE ACCURATE, COMPLETE, RELIABLE OR CORRECT; (B) THE + MATERIALS FILES WILL BE SECURE; (C) THE MATERIALS WILL BE AVAILABLE + AT ANY PARTICULAR TIME OR LOCATION; (D) ANY DEFECTS OR ERRORS WILL + BE CORRECTED; (E) THE MATERIALS AND ACCOMPANYING FILES ARE FREE OF + VIRUSES OR OTHER HARMFUL COMPONENTS; OR (F) THE RESULTS OF USING THE + MATERIALS WILL MEET DOWNLOADER’S REQUIREMENTS. DOWNLOADER’S USE OF + THE MATERIALS IS SOLELY AT DOWNLOADER’S OWN RISK." + ) + ), + tags$li( + tags$h5("Limitation of Liability"), + tags$p( + "IN NO EVENT SHALL LICENSOR BE LIABLE UNDER CONTRACT, TORT, STRICT + LIABILITY, NEGLIGENCE OR ANY OTHER LEGAL THEORY WITH RESPECT TO THE + MATERIALS (I) FOR ANY DIRECT DAMAGES, OR (II) FOR ANY LOST PROFITS + OR SPECIAL, INDIRECT, INCIDENTAL, PUNITIVE, OR CONSEQUENTIAL DAMAGES + OF ANY KIND WHATSOEVER." + ) + ), + tags$li( + tags$h5("Indemnification"), + tags$p( + "You will indemnify and hold Licensor harmless from and against any + and all loss, cost, expense, liability, or damage, including, without + limitation, all reasonable attorneys’ fees and court costs, arising + from i) Your misuse of the use of the Application DiffBrainNet; (ii) + Your violation of these Terms of Use; or (iii) infringement by You + or any third party of any intellectual property or other right of + any person or entity contained in the use of the Application DiffBrainNet. + Such losses, costs, expenses, damages, or liabilities shall include, + without limitation, all actual, general, special, and consequential + damages." + ) + ) + ) + ) + ), easyClose = TRUE, size = "l", footer = NULL) + ) + }) + + observeEvent(input$show, { + showModal(modalDialog( + #HTML('Graphical Abstract'), + # tags$iframe(style="height:600px; width:100%", src="graphicalabstract_new.pdf"), + tags$img(src="graphicalabstract_new.svg", + width="100%", height= "100%"), + footer=tags$div( + tags$b("Graphical Abstract:"), + tags$p("Schematic representation of experimental and analytical steps. + DiffBrainNet is a resource of differential expression and differential + networks in 8 mouse brain regions. (Experiment) C57Bl/6 mice were treated + intraperitoneally with 10mg/kg Dexamethasone or 0.9% saline as vehicle + for 4hours. Eight different brain regions were isolated: amygdala – AMY, + cerebellar cortex – CER, prefrontal cortex – PFC, paraventricular nucleus + of the hypothalamus – PVN, dorsal Cornu Ammonis 1 – dCA1, ventral Cornu + Ammonis 1 – vCA1, dorsal dentate gyrus – dDG, ventral dentate gyrus – vDG. + (Analysis) We performed RNA sequencing in the 8 brain regions, followed + by differential expression analysis (DE) and differential prior-knowledge-based + genome-wide network analysis (DN). (Results) DiffBrainNet includes + differential expression results and network results for all brain regions.")), + size="l", + easyClose = TRUE + )) + }) + + + + ### DIFFERENTIAL EXPRESSION ---------------------------------- + + # DE data set + de_data <- reactive({ + # read DE table to create volcano plot + table <- fread(paste0("tables/de/02_", + input$region,"_deseq2_Dex_1_vs_0_lfcShrink.txt")) + # table$Gene_Symbol <- mapIds(org.Mm.eg.db, keys = table$V1, + # keytype = "ENSEMBL", column="SYMBOL") + table <- table %>% + dplyr::select(Gene_Symbol, Ensembl_ID:padj) %>% + dplyr::mutate_at(vars(baseMean, log2FoldChange, lfcSE), ~round(.,4)) %>% + dplyr::mutate_at(vars(pvalue, padj), ~signif(.,6)) + }) + + + # Volcano Plot displaying DE results + output$volcano_plot <- renderPlotly({ + + # req(de_data()) + + # add column that indicates if gene is significant + indices <- input$de_table_rows_all + table <- de_data()[indices,] %>% + dplyr::mutate(sig = as.factor(ifelse(padj <= 0.1, "FDR <= 0.1", "FDR > 0.1"))) + + # volcano plot + # volcano_plot <- plot_ly(data = table, x = ~log2FoldChange, y = ~-log10(padj), + # color = ~sig, colors = "Set1", + # text = ~paste0("Gene: ", V1), + # source = "volcano_plot") %>% + # layout(title = paste0("Volcano Plot ", input$region)) + volcano_plot <- + ggplot(data = table, aes( + x = log2FoldChange, + y = -log10(padj), + col = sig, + text = paste0("Ensembl_ID: ", Ensembl_ID, "\nGene_Symbol: ", Gene_Symbol) + )) + + geom_point() + + scale_colour_manual(breaks = c("FDR <= 0.1", "FDR > 0.1"), + values = c("orange", "darkgrey")) + + ylab("-log10(FDR)")+ + theme_light() + + theme(plot.title = element_text(size = 11), + legend.title = element_blank()) + + ggtitle(label = paste0("Volcano Plot ", input$region)) + ggplotly(volcano_plot, + source = "volcano_plot") + }) + + # Boxplot displaying norm expression values + output$exp_plot <- renderPlotly({ + + # req(de_data()) + + # access data from click event + d <- event_data("plotly_click", source = "volcano_plot", + priority = "event") + # print(d) + # don't show anythiny if no point was clicked + if (is.null(d)) return(NULL) + + # identify ensembl ID using the DE data + ensembl_id <- de_data()[input$de_table_rows_all,]$Ensembl_ID[d$pointNumber + 1] # use pointNumber/rowNumber of clicked point + gene_symbol <- de_data()[input$de_table_rows_all,]$Gene_Symbol[d$pointNumber + 1] + + # read table with normalized expression values + exp_table <- fread(paste0("tables/de/02_", + input$region,"_deseq2_expression_vsd.txt")) + + # subset to row of clicked gene and reformat + exp_gene <- data.table::transpose(exp_table[V1 == ensembl_id,2:ncol(exp_table)], keep.names = "condition") %>% + dplyr::rename("expression" = "V1") + exp_gene$condition <- str_replace(exp_gene$condition, ".*\\_","") + exp_gene$mouse_id <- str_extract(exp_gene$condition, pattern = "[0-9]+") + exp_gene$condition <- as.factor(str_replace(exp_gene$condition, "[0-9]+","")) + + # boxplot with jitter points of expression values in CNTRL and DEX + exp_plot <- + ggplot(exp_gene, aes(x = condition, y = expression, fill = condition)) + + geom_boxplot() + + geom_jitter(color = "black", + size = 0.4, + alpha = 0.9, + aes(text = sprintf('mouse_ID: %s', mouse_id))) + + scale_fill_manual("", + breaks = c("CNTRL", "DEX"), + labels = c("Vehicle", "Treatment"), + values = c("#B0BFBB", "#46866E")) + + scale_x_discrete(breaks = c("CNTRL", "DEX"), + labels = c("Vehicle", "Treatment")) + + theme_light() + + theme(legend.position = "none", + plot.title = element_text(size = 11)) + + ggtitle(paste0("Gene expression of ", gene_symbol, " in ", input$region)) + + xlab("") + + ylab("Normalized expression") + # ggplot to plotly + ggplotly(exp_plot) + + }) + + + # Data table with DE results + output$de_table <- DT::renderDataTable({ + + # display information when hovering column names + hover_info <- c("Gene ID in gene symbol format", + "Gene ID in Ensembl format", + "Mean expression in vehicle condition", + "Logarithmized fold change of gene expression between control and treatment condition", + "Standard Error of the fold change", + "Respective p-value of the fold change", + "FDR corrected p-value") + + de_data() %>% + # dplyr::rename("Ensembl_ID" = "V1") %>% + datatable(filter = list(position = 'top'), + options = list( + headerCallback= JS(callHeaderCallback(hover_info)) + )) + }, server = TRUE) # FALSE to enable download of all pages with button + + # Download of (filtered) DE results + output$download1 <- downloadHandler( + filename = function() { + paste0("DE_dexStimMouse_", input$region, ".csv") + }, + content = function(file) { + indices <- input$de_table_rows_all + write.csv(de_data()[indices, ], file) + } + ) + + # Upset plot and table from upsetPlot module + upset_de <- upsetPlotServer("upsetDE") + + + ############# NETWORK ANALYSIS ############################## + + + # Network data (nodedegree/nodebetweenness) + overview_data <- reactive({ + + # Read nodedegrees/nodebetweenness + filename <- list.files(path = paste0("tables/network"), + pattern = paste0("*_", input$overview_region, + "_funcoup_", input$overview_network, "_", + input$overview_metric,"Norm_betacutoff0.01.csv"), + full.names = TRUE) + data <- fread(file = filename) + + }) + + # Upset plot and table from upsetPlot module + upset_hub <- comparisonHubGenesServer("compHubGenes") + + + # Histogram network overview + output$histogram_network <- renderPlotly({ + + req(overview_data()) + + # Histogram of nodedegree/nodebetweenness + if (input$overview_metric == "nodedegrees"){ + hist_network <- ggplot(overview_data(), aes(x=nodedegree)) + + geom_histogram(bins=80) + } else { + hist_network <- ggplot(overview_data(), aes(x=nodebetweenness)) + + geom_histogram(bins=80) + } + + ggplotly(hist_network) + }) + + ### SINGLE REGION NETWORK ------------------------------- + + # Single region network and table from networkSingle module + net_single <- networkSingleServer("singleVisualization") + + + ### MULTI REGION NETWORK --------------------------------------- + + # Multi region network and table from networkMulti module + net_multi <- networkMultiServer("multiVisualization") + + +}) + diff --git a/ui.R b/ui.R new file mode 100644 index 0000000..29b228c --- /dev/null +++ b/ui.R @@ -0,0 +1,233 @@ + +ui <- fluidPage(theme = shinytheme("flatly"), +titlePanel(title=div(img(src="DiffBrainNet_logo.png", height = 50, width = 200), "Glucocorticoid receptor regulated gene expression in the mouse brain")), +navbarPage("Explore!", + + tabPanel(icon("home"), + tags$head(tags$style(".modal-footer {text-align: left;}")), + fluidRow(column(tags$img(src="mouse_javi_mpi.png", + width="100%", height= "100%"), + width=floor(0.4*12)), + column( + p(strong("Abstract."), "Network analysis can identify the molecular connectivity + that is underpinning function at the control level, after a stimulus or a disease state. We inferred + regression-based prior-knowledge-guided gene networks in 8 brain regions of the mouse brain: + the prefrontal cortex, the amygdala, the paraventricular nucleus of the hypothalamus, the dorsal + and ventral Cornu ammonis 1, the dorsal and ventral dentate gyrus and the cerebellar cortex. + We constructed networks at the control and treatment level using KiMONo (Ogris et al.) and at + the differential level using DiffGRN (Kim et al.).",br(), + "As a stimulus we used dexamethasone, a synthetic glucocorticoid that is used + to activate the glucocorticoid receptors. Glucocorticoid receptors, when coupled with glucocorticoids + like dexamethasone, act as transcription factors modulating the transcriptional landscape. + We provide differential networks and differential + expression analysis (DESeq2, Love et al.) that can be used to analyse the effects of dexamethasone + both at the molecular connectivity and at the gene level in each brain region.", br(), + "DiffBrainNet is an analysis framework and a resource for studying the transcriptional landscape + of 8 mouse brain regions at the control, dexamethasone-treatment and differential level. It can be + used to pinpoint molecular pathways important for the basic function and response to glucocorticoids + in a brain-region specific manner. DiffBrainNet can also support the identification and analysis + of biological processes regulated by brain and psychiatric diseases risk genes at the control and + differential levels.", br(), + actionButton("show", "Show graphical abstract", + style="color: #fff; background-color: #2980B9; border-color: #2e6da4"), + style="text-align:justify;color:#fff;background-color:#2980B9;padding:15px;border-radius:10px"), + br(), + p(strong("Data and code availability."), "Raw and normalized gene expression data generated in this + study are provided at GEO under", a(href="https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190712", "GSE190712", target="_blank"), + ". Differential expression and differential network data can be downloaded from this resource.", + "Data analysis scripts and source code of this shiny app is available on our", + a(href="https://github.molgen.mpg.de/mpip/DiffBrainNet", "github page", target="_blank"),".", + style="text-align:justify;color:black;background-color:#AED6F1;padding:15px;border-radius:10px"), + width = floor(0.45*12) + ), + column( + br(), + tags$img(src="mpilogo.jpeg", width="80%", height= "80%"), + br(), + br(), + p("Department of Translational Research in Psychiatry"), + p("Medical Genomics Group"), + p("For more information please visit the website of the", #em("the MPI of Psychiatry"), + a(href="https://www.psych.mpg.de/", em("MPI of Psychiatry"), target="_blank")), + width=ceiling(0.15*12)+1 + ))), + + navbarMenu("Diff. Expression Analysis", + tabPanel("Single Brain Region", + sidebarLayout( + sidebarPanel( + selectInput("region", "Choose a brain region:", + choices = c("Amygdala" = "AMY", + "Cerebellum" = "CER", + "Dorsal DG of Hippocampus" = "dDG", + "Dorsal CA1 of Hippocampus" = "dCA1", + "Prefrontal Cortex" = "PFC", + "PVN of Hypothalamus" = "PVN", + "Ventral DG of Hippocampus" = "vDG", + "Ventral CA1 of Hippocampus" = "vCA1")), + downloadButton("download1","Download (filtered) table as csv"), + width = 3 + ), + mainPanel( + fluidPage( + br(), + tags$style(".fa-chart-bar {color:#2980B9}"), + h3(p(em("Differential Gene Expression "),icon("chart-bar",lib = "font-awesome"),style="color:black;text-align:center")), + hr(), + span("Explore the transcriptomic response to Dexamethasone treatment of a brain region of your + interest on a "), strong("differential expression level"), span(". You can select one of 8 + different brain regions on the left panel. The volcano plot shows the log2-transformed fold change + and the -log10-transformed FDR corrected p-value of the genes detected in our dataset. Please "), + strong("click on a point/gene"), span(" in the volcano plot to see its + normalized expression levels at the vehicle and treatment condition. You can filter and + download the table of the differentially expressed genes."), + br(),br(), + # style="text-align:left;color:black"), + splitLayout(cellWidths = c("55%", "45%"), + plotlyOutput("volcano_plot"), + plotlyOutput("exp_plot")), + br(),br(), + span(em("Hint: By hovering over the column names you get a more detailed description of the columns.")), + br(),br(), + DT::dataTableOutput("de_table") + ) + ) + ) + ), + tabPanel("Comparison Brain Regions", + # UI defined in upsetPlot module + upsetPlotUI("upsetDE") + ) + ), + + navbarMenu( + "Diff. Network Analysis", + tabPanel( + "Introduction diff. networks", + fluidPage( + br(), + tags$style(".fa-project-diagram {color:#2980B9}"), + h3(p( + em("Introduction: Differential network analysis"), + icon("project-diagram", lib = "font-awesome"), + style = "color:black;text-align:center" + )), + hr(), + + tags$img(src="DiffNetworks.png", + width="60%", height= "60%", + style="display: block; margin-left: auto; margin-right: auto;") + ) + + + + ), + + tabPanel( + "Network overview", + sidebarLayout( + sidebarPanel( + selectInput( + "overview_region", + "Choose a brain region:", + choices = c( + "Amygdala" = "AMY", + "Cerebellum" = "CER", + "Dorsal DG of Hippocampus" = "dDG", + "Dorsal CA1 of Hippocampus" = "dCA1", + "Prefrontal Cortex" = "PFC", + "PVN of Hypothalamus" = "PVN", + "Ventral DG of Hippocampus" = "vDG", + "Ventral CA1 of Hippocampus" = "vCA1" + ) + ), + radioButtons( + "overview_metric", + label = "Select network metric:", + choices = list("Nodedegree" = "nodedegrees", + "Nodebetweenness" = "nodebetweenness"), + selected = "nodebetweenness" + ), + radioButtons( + "overview_network", + label = "Select network:", + choices = list("Vehicle" = "baseline", + "Differential" = "differential", + "Treatment" = "treatment"), + selected = "differential" + ), + width = 3 + ), + mainPanel( + fluidPage( + br(), + tags$style(".fa-project-diagram {color:#2980B9}"), + h3(p( + em("Network analysis of gene expression: Overview "), + icon("project-diagram", lib = "font-awesome"), + style = "color:black;text-align:center" + )), + hr(), + h4( + p("Distribution nodedegrees/nodebetweenness", + style = "color:black;text-align:center") + ), + plotlyOutput("histogram_network") + ) + ) + )), + + tabPanel( + "Comparison hub genes", + comparisonHubGenesUI("compHubGenes") + ), + + tabPanel( + "Network visualization", + # UI from networkSingle module + networkSingleUI("singleVisualization") + ), + + tabPanel( + "Multi-region network visualization", + # UI from networkMulti module + networkMultiUI("multiVisualization") + ) + ), + + + # navbarMenu("More", + # tabPanel("Data download", + # # DT::dataTableOutput("table") + # ), + # tabPanel("Terms of Use", + # + # # fluidRow( + # # column(6, + # # includeMarkdown("about.md") + # # ), + # # column(3, + # # img(class="img-polaroid", + # # src=paste0("http://upload.wikimedia.org/", + # # "wikipedia/commons/9/92/", + # # "1919_Ford_Model_T_Highboy_Coupe.jpg")), + # # tags$small( + # # "Source: Photographed at the Bay State Antique ", + # # "Automobile Club's July 10, 2005 show at the ", + # # "Endicott Estate in Dedham, MA by ", + # # a(href="http://commons.wikimedia.org/wiki/User:Sfoskett", + # # "User:Sfoskett") + # # ) + # # ) + # # ) + # ) + # ), + footer = tags$footer( + p(tags$span(actionLink(inputId = "termsOfUse", label = "Terms of Use")), + " ", + a(href="https://www.psych.mpg.de/2354/impress", "Imprint", target="_blank"), + " ", + " © 2022 MPI of Psychiatry"), + align="center") +) +) diff --git a/utilities_network.R b/utilities_network.R new file mode 100644 index 0000000..0d449fc --- /dev/null +++ b/utilities_network.R @@ -0,0 +1,178 @@ +### FUNCTIONS ------------------------------------- + + +# remove duplicated edges from network dataframe +simplify_network_df <- function(network_dt){ + + network_dt <- network_dt %>% + dplyr::filter(from != to) #%>% + # dplyr::mutate(from_sort = pmin(from, to), to_sort = pmax(from, to)) %>% + # dplyr::distinct(from_sort, to_sort, .keep_all = TRUE) %>% + # dplyr::select(!c(from_sort, to_sort)) +} + + +# subset network to gene of interest and if required also +# neighbours of genes +subset_network <- function(network_dt, subset_gene, neighbours){ + + if (neighbours) { + e_interest <- network_dt %>% + filter(from %in% subset_gene | to %in% subset_gene) + e_interest <- unique(c(e_interest$from, e_interest$to)) + network_subset <- network_dt %>% + filter(from %in% e_interest & to %in% e_interest) + } else { + network_subset <- network_dt %>% + filter(from %in% subset_gene & to %in% subset_gene) + } + + return(network_subset) +} + + +# function to generate baseline/diff network, single regions +read_network <- function(data, de_genes, hub_genes, gene, mode, id_type, neighbours){ + + # input genes + print(gene) + + # Edges + if(mode == "baseline" | mode == "treatment"){ + c <- rep("grey", nrow(data)) + relations <- data.table(from=data$predictor, + to=data$target, + beta = data$beta_mean, + color = c) + # width = rescale(abs(data$beta_mean), to = c(0,1))) + ledges <- data.frame(color = c("grey"), + label = c("beta"), + #arrows = c("right", "right"), + font.align = "top") + } else { + c <- rep("#db9512", nrow(data)) + c[data$z > 0] <- "#128bdb" + relations <- data.table(from=data$predictor, + to=data$target, + z = data$z, + p_adj = data$p_adj, + color = c) + ledges <- data.frame(color = c("#db9512", "#128bdb"), + label = c("z < 0", "z > 0"), + arrows = c("right", "right"), + font.align = "top") + } + relations <- simplify_network_df(relations) + relations <- subset_network(relations, gene, neighbours) + print(nrow(relations)) + + # Nodes + # Generate dataframe for node layout + node_vec <- unique(c(relations$from, relations$to)) + if (id_type == "Gene Symbol"){ + labels <- mapIds(org.Mm.eg.db, keys = node_vec, + column = "SYMBOL", keytype = "ENSEMBL") + } else { + labels <- node_vec + } + + # set colours according to DE status of gene + col <- rep("darkblue", length(node_vec)) + col[node_vec %in% de_genes] <- "orange" + col[node_vec %in% hub_genes] <- "darkred" + col[node_vec %in% de_genes & node_vec %in% hub_genes] <- "purple" + nodes <- data.table(id = node_vec, + label = labels, + color = col) + + # nodes data.frame for legend + lnodes <- data.frame(label = c("DE", "hub", "DE & hub", "no DE & no hub"), + font.color = c("black", "white", "white", "white"), + shape = c( "ellipse"), + color = c("orange", "darkred", "purple", "darkblue"), + title = "Informations", id = 1:4) + + return(list("edges" = relations, "nodes" = nodes, "ledges" = ledges, "lnodes" = lnodes)) + +} + + +# function to generate baseline/diff network, multiple regions +read_network_multi <- function(regions, gene_unsplit, mode, id_type){ + + list_reg <- list() + + for (reg in regions){ + + # Read data + data <- + fread( + file = paste0( + basepath, + "tables/coExpression_kimono/03_AnalysisFuncoup/", + "04_singleRegion_", + reg, + "_filtered_diffNetwork.csv" + ) + ) %>% + dplyr::select(target, predictor, beta_mean.base, beta_mean.dex, z, p_adj) + + cols <- colnames(data)[3:6] + data <- data %>% + dplyr::rename_with(.fn = ~paste0(., ".", reg), .cols = all_of(cols) ) + + list_reg[[reg]] <- data + + } + + data_joined <- list_reg %>% + purrr::reduce(full_join, by = c("target","predictor")) + + # # Read data + # data <- fread(file = paste0(basepath, "tables/coExpression_kimono/03_AnalysisFuncoup/", + # "04_singleRegion_",region,"_filtered_", mode, "Network.csv")) + gene <- stringr::str_split(gene_unsplit, ",\\s?")[[1]] + print(gene) + if (!startsWith(gene_unsplit, "ENSMUSG")){ + gene <- mapIds(org.Mm.eg.db, keys = gene, column = "ENSEMBL", keytype = "SYMBOL") + } + + # Edges + c <- rep("#db9512", nrow(data)) + c[data$z > 0] <- "#128bdb" + relations <- data.table(from=data$target, + to=data$predictor, + z = data$z, + p_adj = data$p_adj, + color = c) + # print(relations) + relations <- simplify_network_df(relations) + relations <- subset_network(relations, gene) + print(nrow(relations)) + + # Nodes + node_vec <- unique(c(relations$from, relations$to)) + if (id_type == "Gene Symbol"){ + labels <- mapIds(org.Mm.eg.db, keys = node_vec, + column = "SYMBOL", keytype = "ENSEMBL") + } else { + labels <- node_vec + } + nodes <- data.table(id = node_vec, + label = labels) + # title = node_vec) + return(list("edges" = relations, "nodes" = nodes)) + +} + + +# function to generate network stats +network_stats <- function(df){ + + # Initialize data frame for stats + stats <- data.frame(#network = character(), + type = character(), + value = numeric()) + + stats <- rbind(stats, list("type" = "edges", "value" = nr_edges_base)) +} diff --git a/www/.DS_Store b/www/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/www/.DS_Store differ diff --git a/www/DiffNetworks.pptx b/www/DiffNetworks.pptx new file mode 100644 index 0000000..9c154e3 Binary files /dev/null and b/www/DiffNetworks.pptx differ diff --git a/www/geo_main.gif b/www/geo_main.gif new file mode 100644 index 0000000..2ea819a Binary files /dev/null and b/www/geo_main.gif differ diff --git a/www/mousejavi b/www/mousejavi new file mode 100644 index 0000000..78c90f1 Binary files /dev/null and b/www/mousejavi differ diff --git a/www/mousejavi_reversebrain b/www/mousejavi_reversebrain new file mode 100644 index 0000000..166e1ec Binary files /dev/null and b/www/mousejavi_reversebrain differ