diff --git a/06_Shiny/.DS_Store b/06_Shiny/.DS_Store deleted file mode 100644 index b455ec8..0000000 Binary files a/06_Shiny/.DS_Store and /dev/null differ diff --git a/06_Shiny/.Rproj.user/.DS_Store b/06_Shiny/.Rproj.user/.DS_Store deleted file mode 100644 index 7651473..0000000 Binary files a/06_Shiny/.Rproj.user/.DS_Store and /dev/null differ diff --git a/06_Shiny/.Rproj.user/94131CCE/pcs/files-pane.pper b/06_Shiny/.Rproj.user/94131CCE/pcs/files-pane.pper deleted file mode 100644 index 7935ff1..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/pcs/files-pane.pper +++ /dev/null @@ -1,9 +0,0 @@ -{ - "sortOrder": [ - { - "columnIndex": 2, - "ascending": true - } - ], - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/pcs/source-pane.pper b/06_Shiny/.Rproj.user/94131CCE/pcs/source-pane.pper deleted file mode 100644 index bc7dc99..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/pcs/source-pane.pper +++ /dev/null @@ -1,3 +0,0 @@ -{ - "activeTab": 4 -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/pcs/windowlayoutstate.pper b/06_Shiny/.Rproj.user/94131CCE/pcs/windowlayoutstate.pper deleted file mode 100644 index 104979b..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/pcs/windowlayoutstate.pper +++ /dev/null @@ -1,14 +0,0 @@ -{ - "left": { - "splitterpos": 312, - "topwindowstate": "NORMAL", - "panelheight": 743, - "windowheight": 781 - }, - "right": { - "splitterpos": 468, - "topwindowstate": "NORMAL", - "panelheight": 743, - "windowheight": 781 - } -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/pcs/workbench-pane.pper b/06_Shiny/.Rproj.user/94131CCE/pcs/workbench-pane.pper deleted file mode 100644 index 07157f3..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/pcs/workbench-pane.pper +++ /dev/null @@ -1,5 +0,0 @@ -{ - "TabSet1": 0, - "TabSet2": 4, - "TabZoom": {} -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/rmd-outputs b/06_Shiny/.Rproj.user/94131CCE/rmd-outputs deleted file mode 100644 index 3f2ff2d..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/rmd-outputs +++ /dev/null @@ -1,5 +0,0 @@ - - - - - diff --git a/06_Shiny/.Rproj.user/94131CCE/saved_source_markers b/06_Shiny/.Rproj.user/94131CCE/saved_source_markers deleted file mode 100644 index 2b1bef1..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/saved_source_markers +++ /dev/null @@ -1 +0,0 @@ -{"active_set":"","sets":[]} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/1532F1BD b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/1532F1BD deleted file mode 100644 index 0a412e4..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/1532F1BD +++ /dev/null @@ -1,24 +0,0 @@ -{ - "id": "1532F1BD", - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/R/network_multiRegion.R", - "project_path": "R/network_multiRegion.R", - "type": "r_source", - "hash": "4109658605", - "contents": "", - "dirty": false, - "created": 1636637444100.0, - "source_on_save": false, - "relative_order": 5, - "properties": { - "cursorPosition": "219,31", - "scrollLine": "211" - }, - "folds": "", - "lastKnownWriteTime": 1636637543, - "encoding": "UTF-8", - "collab_server": "", - "source_window": "", - "last_content_update": 1636637543028, - "read_only": false, - "read_only_alternatives": [] -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/1532F1BD-contents b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/1532F1BD-contents deleted file mode 100644 index 402cda3..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/1532F1BD-contents +++ /dev/null @@ -1,494 +0,0 @@ -# module UI function -networkMultiUI <- function(id, label = "Network Multiple Regions"){ - - ns <- NS(id) - - tagList( - - 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("Baseline", "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:", - choices = list("Ensembl", "Gene Symbol"), - selected = "Gene Symbol" - ), - checkboxGroupInput( - ns("tableContent"), - label = "Table content:", - 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(p("Differential network", style = "color:black;text-align:center")), - visNetworkOutput(ns("network_diff_multi"), height = "600px") - %>% withSpinner(color="#0dc5c1"), - DT::dataTableOutput(ns("network_table")) - ) - ) - - - ) -} - - - -# module server function -networkMultiServer <- function(id) { - moduleServer( - id, - - function(input, output, session) { - - ### MULTI REGION NETWORK ------------------------------- - - # 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) - }) - - - # #### TEST - # - # list_network <- list() - # - # for (reg in c("AMY", "CER", "PFC")){ - # file_path <- paste0( - # basepath, - # "tables/coExpression_kimono/03_AnalysisFuncoup/", - # "04_singleRegion_", - # reg, - # "_filtered_", - # "baseline", - # "Network.csv" - # ) - # # Read data - # data <- - # fread(file = file_path) %>% - # #dplyr::select(target, predictor, beta_mean.base, beta_mean.dex, z, p_adj) - # dplyr::select(target, predictor, beta_mean, beta_mean.dex, z, p_adj) - # - # cols <- colnames(data)[3:6] - # #cols <- colnames(data)[3] - # data <- data %>% - # 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")) - - - # 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({ - - network_table() %>% - datatable(filter = list(position = 'top')) - }, 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 = c, - title = data$title) - # df for edge legend - ledges <- data.frame(color = edge_colors, - label = 1:(ncol(data)-6), - font.align = "top") - } 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 = 1:(ncol(data)-6), - font.align = "top") - } - 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"), - shape = c( "ellipse"), color = c("orange", "darkred", "purple", "darkblue"), - title = "Informations", id = 1:4) - - - # return(list("edges" = relations, "nodes" = nodes, "ledges" = ledges, "lnodes" = lnodes)) - return(list("edges" = relations, "nodes" = nodes, "ledges" = ledges)) - - } - - - - # prepare network table for visualization - network_table <- reactive({ - - req(network_data()) - req(input_genes()) - - print(input$tableContent) - - - # 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) - - # 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")) - } - - } - - data - - }) - - - } - - ) -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/3F51A6BD b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/3F51A6BD deleted file mode 100644 index 91cbb3f..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/3F51A6BD +++ /dev/null @@ -1,24 +0,0 @@ -{ - "id": "3F51A6BD", - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/R/network_singleRegion.R", - "project_path": "R/network_singleRegion.R", - "type": "r_source", - "hash": "1168392640", - "contents": "", - "dirty": false, - "created": 1636637354945.0, - "source_on_save": false, - "relative_order": 4, - "properties": { - "cursorPosition": "129,46", - "scrollLine": "115" - }, - "folds": "", - "lastKnownWriteTime": 1636637425, - "encoding": "UTF-8", - "collab_server": "", - "source_window": "", - "last_content_update": 1636637425157, - "read_only": false, - "read_only_alternatives": [] -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/3F51A6BD-contents b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/3F51A6BD-contents deleted file mode 100644 index 9ef0406..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/3F51A6BD-contents +++ /dev/null @@ -1,240 +0,0 @@ -# module UI function -networkSingleUI <- function(id, label = "Network Single Region"){ - - ns <- NS(id) - - tagList( - - 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("Baseline", "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:", - 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(p("Differential network", style = "color:black;text-align:center")), - visNetworkOutput(ns("network"), height = "600px") %>% - withSpinner(color="#0dc5c1"), - DT::dataTableOutput(ns("network_table")) - ) - ) - - ) - -} - - - -# module server function -networkSingleServer <- function(id) { - moduleServer( - id, - - function(input, output, session) { - - ### SINGLE REGION NETWORK ------------------------------- - - # 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({ - - network_data() %>% - datatable(filter = list(position = 'top')) - }, 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/06_Shiny/.Rproj.user/94131CCE/sources/per/t/65326343 b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/65326343 deleted file mode 100644 index e5619a0..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/65326343 +++ /dev/null @@ -1,21 +0,0 @@ -{ - "id": "65326343", - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/app/ui.R", - "project_path": null, - "type": "r_source", - "hash": "2781797785", - "contents": "", - "dirty": false, - "created": 1636637158435.0, - "source_on_save": false, - "relative_order": 2, - "properties": {}, - "folds": "", - "lastKnownWriteTime": 1635167585, - "encoding": "UTF-8", - "collab_server": "", - "source_window": "", - "last_content_update": 1635167585, - "read_only": false, - "read_only_alternatives": [] -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/65326343-contents b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/65326343-contents deleted file mode 100644 index 43ba7e8..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/65326343-contents +++ /dev/null @@ -1,183 +0,0 @@ - -ui <- fluidPage(theme = shinytheme("flatly"), -titlePanel("Glucocorticoid receptor regulated gene expression in the mouse brain"), -navbarPage("Explore!", - - tabPanel(icon("home"), - fluidRow(column(tags$img(src="mousejavi_reversebrain.png", - width="450px",height="320px"), - width=3), - column( - p(strong("Description of project (e.g. paper abstract)."), "Lorem ipsum dolor sit amet, - consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna - aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut - aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate - velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non - proident, sunt in culpa qui officia deserunt mollit anim id est laborum.", - style="text-align:justify;color:black;background-color:#2980B9;padding:15px;border-radius:10px"), - br(), - p(strong("Data and code availability. Reference to paper."), "Sed ut perspiciatis unde omnis iste - natus error sit voluptatem accusantium doloremque laudantium, totam rem aperiam, eaque ipsa quae - ab illo inventore veritatis et quasi architecto beatae vitae dicta sunt explicabo. Nemo enim - ipsam voluptatem quia voluptas sit aspernatur aut odit aut fugit, sed quia consequuntur magni - dolores eos qui ratione voluptatem sequi nesciunt. Neque porro quisquam est, qui dolorem ipsum - quia dolor sit amet, consectetur, adipisci velit, sed quia non numquam eius modi tempora incidunt - ut labore et dolore magnam aliquam quaerat voluptatem. Ut enim ad minima veniam, quis nostrum - exercitationem ullam corporis suscipit laboriosam, nisi ut aliquid ex ea commodi consequatur? - Quis autem vel eum iure reprehenderit qui in ea voluptate velit esse quam nihil molestiae - consequatur, vel illum qui dolorem eum fugiat quo voluptas nulla pariatur?", - style="text-align:justify;color:black;background-color:#AED6F1;padding:15px;border-radius:10px"), - width = 7 - ), - column( - br(), - tags$img(src="mpilogo.png", width="250px", height="60px"), - br(), - br(), - p("For more information please visit the website of", em("the MPI of Psychiatry"), - br(), - a(href="https://www.psych.mpg.de/", "Here", target="_blank"), - style="text-align:center;color:black"), - width=2 - ))), - - navbarMenu("Diff. Gene Expression", - 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(), - splitLayout(cellWidths = c("55%", "45%"), - plotlyOutput("volcano_plot"), - plotlyOutput("exp_plot")), - DT::dataTableOutput("de_table") - ) - ) - ) - ), - tabPanel("Comparison Brain Regions", - # UI defined in upsetPlot module - upsetPlotUI("upsetDE") - ) - ), - - navbarMenu( - "Network Analysis", - tabPanel( - "Overview", - 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("Baseline" = "baseline", - "Differential" = "differential"), - 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(), - splitLayout(cellWidths = c("50%", "50%"), - h4( - p("Baseline network", style = "color:black;text-align:center") - ), - h4( - p("Differential network", style = "color:black;text-align:center") - )), - # plotlyOutput("histogram_network") - splitLayout( - cellWidths = c("50%", "50%"), - plotlyOutput("histogram_network"), - plotlyOutput("barplot_network") - ) - # DT::dataTableOutput("network_table") - ) - ) - ), - - 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("About", - # 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") - # ) - # ) - # ) - ) - ) -) -) diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/951FD869 b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/951FD869 deleted file mode 100644 index a6c3cdb..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/951FD869 +++ /dev/null @@ -1,24 +0,0 @@ -{ - "id": "951FD869", - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/app/global.R", - "project_path": null, - "type": "r_source", - "hash": "2758396765", - "contents": "", - "dirty": false, - "created": 1636637160030.0, - "source_on_save": false, - "relative_order": 3, - "properties": { - "cursorPosition": "20,29", - "scrollLine": "0" - }, - "folds": "", - "lastKnownWriteTime": 1636637200, - "encoding": "UTF-8", - "collab_server": "", - "source_window": "", - "last_content_update": 1636637200978, - "read_only": false, - "read_only_alternatives": [] -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/951FD869-contents b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/951FD869-contents deleted file mode 100644 index 3ea3837..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/951FD869-contents +++ /dev/null @@ -1,21 +0,0 @@ -library(shiny) -library(shinythemes) -library(markdown) -library(plotly) -library(DT) -library(data.table) -library(stringr) -library(dplyr) -library(ggplot2) -library(igraph) -library(visNetwork) -library(shinyWidgets) -library(UpSetR) -library(org.Mm.eg.db) -library(scales) -library(shinycssloaders) -library(RColorBrewer) - - -#basepath <- "~/Documents/ownCloud/DexStim_RNAseq_Mouse/" -source("utilities_network.R") \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/95FC14D0 b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/95FC14D0 deleted file mode 100644 index 95179d0..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/95FC14D0 +++ /dev/null @@ -1,24 +0,0 @@ -{ - "id": "95FC14D0", - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/app/server.R", - "project_path": null, - "type": "r_source", - "hash": "898784622", - "contents": "", - "dirty": false, - "created": 1636637147649.0, - "source_on_save": false, - "relative_order": 4, - "properties": { - "cursorPosition": "130,56", - "scrollLine": "120" - }, - "folds": "", - "lastKnownWriteTime": 1636637341, - "encoding": "UTF-8", - "collab_server": "", - "source_window": "", - "last_content_update": 1636637341969, - "read_only": false, - "read_only_alternatives": [] -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/95FC14D0-contents b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/95FC14D0-contents deleted file mode 100644 index 192f63d..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/95FC14D0-contents +++ /dev/null @@ -1,189 +0,0 @@ -#source(paste0(basepath, "scripts/06_Shiny/utilities_network.R")) -#source(paste0(basepath, "scripts/06_Shiny/utilities_de.R")) - -# Define server logic required to generate and plot -server <- shinyServer(function(input, output) { - - ### 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, "p_adj <= 0.1", "p_adj > 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_brewer(palette = "Set1") + - 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)) + - geom_boxplot() + - geom_jitter(color = "black", - size = 0.4, - alpha = 0.9, - aes(text = sprintf('mouse_ID: %s', mouse_id))) + - # scale_fill_manual(values = c("#1F6D84", "#E6A54D")) + - 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({ - de_data() %>% - # dplyr::rename("Ensembl_ID" = "V1") %>% - datatable(filter = list(position = 'top')) - }, 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/", - "03_AnalysisFuncoup"), - pattern = paste0("*_", input$overview_region, - "_funcoup_", input$overview_network, "_", - input$overview_metric,"Norm_betacutoff0.01.csv"), - full.names = TRUE) - data <- fread(file = filename) - - }) - - - # Histogram network overview - output$histogram_network <- renderPlotly({ - - req(overview_data()) - - # Histogram of nodedegree/nodebetweenness - if (input$overview_metric == "nodedegrees"){ - histogram_network <- ggplot(overview_data(), aes(x=nodedegree)) + - geom_histogram() - } else { - histogram_network <- ggplot(overview_data(), aes(x=nodebetweenness)) + - geom_histogram() - } - ggplotly(histogram_network) - }) - - # Barplot network overview - output$barplot_network <- renderPlotly({ - - req(overview_data()) - - # Histogram of nodedegree/nodebetweenness - if (input$overview_metric == "nodedegrees"){ - barplot_network <- ggplot(overview_data(), aes(x=nodedegree)) - } else { - barplot_network <- ggplot(overview_data(), aes(x=nodebetweenness)) + - geom_histogram() - } - - barplot_network <- barplot_network + - geom_bar() - ggplotly(barplot_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/06_Shiny/.Rproj.user/94131CCE/sources/per/t/EB039850 b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/EB039850 deleted file mode 100644 index df31e4c..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/EB039850 +++ /dev/null @@ -1,24 +0,0 @@ -{ - "id": "EB039850", - "path": "~/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/R/upset_plot.R", - "project_path": "R/upset_plot.R", - "type": "r_source", - "hash": "204956232", - "contents": "", - "dirty": false, - "created": 1636637642141.0, - "source_on_save": false, - "relative_order": 6, - "properties": { - "cursorPosition": "104,46", - "scrollLine": "92" - }, - "folds": "", - "lastKnownWriteTime": 1636637694, - "encoding": "UTF-8", - "collab_server": "", - "source_window": "", - "last_content_update": 1636637694314, - "read_only": false, - "read_only_alternatives": [] -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/EB039850-contents b/06_Shiny/.Rproj.user/94131CCE/sources/per/t/EB039850-contents deleted file mode 100644 index b11b130..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/per/t/EB039850-contents +++ /dev/null @@ -1,421 +0,0 @@ -# module UI function -upsetPlotUI <- function(id, label = "Upset Plot"){ - - ns <- NS(id) - - #tabPanel("Comparison Brain Regions", - tagList( - 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 - ), - 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; }" - ), - # plotlyOutput("upset_de") - #plotOutput("upset_de"), - plotlyOutput(ns("plotly_upset"), height = "600px"), - DT::dataTableOutput(ns("de_comp_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", "Gene_Symbol", "region")] %>% - group_by(Ensembl_ID, Gene_Symbol) %>% - dplyr::summarise(region = list(region)) %>% - dplyr::select(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) - } - - - # Data table below upset plot - output$de_comp_table <- DT::renderDataTable({ - - req(input$myPicker) - data <- upset_datatable(basepath, input$myPicker, input$padj_upset, input$FC_upset) - }) - - - # 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 = "black", width = 3), - marker = list(color = "black", - 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 = "black"), - 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 = "black", - 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/06_Shiny/.Rproj.user/94131CCE/sources/prop/0C9B8A0A b/06_Shiny/.Rproj.user/94131CCE/sources/prop/0C9B8A0A deleted file mode 100644 index 9e26dfe..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/0C9B8A0A +++ /dev/null @@ -1 +0,0 @@ -{} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/0DDD84F5 b/06_Shiny/.Rproj.user/94131CCE/sources/prop/0DDD84F5 deleted file mode 100644 index 2b1fe59..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/0DDD84F5 +++ /dev/null @@ -1,4 +0,0 @@ -{ - "cursorPosition": "104,46", - "scrollLine": "92" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/0E49E7B8 b/06_Shiny/.Rproj.user/94131CCE/sources/prop/0E49E7B8 deleted file mode 100644 index 9e26dfe..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/0E49E7B8 +++ /dev/null @@ -1 +0,0 @@ -{} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/1532276C b/06_Shiny/.Rproj.user/94131CCE/sources/prop/1532276C deleted file mode 100644 index 0af093c..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/1532276C +++ /dev/null @@ -1,4 +0,0 @@ -{ - "cursorPosition": "129,46", - "scrollLine": "115" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/15E20D6C b/06_Shiny/.Rproj.user/94131CCE/sources/prop/15E20D6C deleted file mode 100644 index 9e26dfe..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/15E20D6C +++ /dev/null @@ -1 +0,0 @@ -{} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/8BA97ED7 b/06_Shiny/.Rproj.user/94131CCE/sources/prop/8BA97ED7 deleted file mode 100644 index 8096394..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/8BA97ED7 +++ /dev/null @@ -1,4 +0,0 @@ -{ - "cursorPosition": "20,29", - "scrollLine": "0" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/8D825508 b/06_Shiny/.Rproj.user/94131CCE/sources/prop/8D825508 deleted file mode 100644 index 5aa7743..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/8D825508 +++ /dev/null @@ -1,4 +0,0 @@ -{ - "cursorPosition": "124,2", - "scrollLine": "122" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/9889CFFF b/06_Shiny/.Rproj.user/94131CCE/sources/prop/9889CFFF deleted file mode 100644 index 6426554..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/9889CFFF +++ /dev/null @@ -1,4 +0,0 @@ -{ - "cursorPosition": "219,31", - "scrollLine": "211" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/CD86E44A b/06_Shiny/.Rproj.user/94131CCE/sources/prop/CD86E44A deleted file mode 100644 index 0bcef40..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/CD86E44A +++ /dev/null @@ -1,4 +0,0 @@ -{ - "cursorPosition": "130,56", - "scrollLine": "120" -} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/CFBC27D1 b/06_Shiny/.Rproj.user/94131CCE/sources/prop/CFBC27D1 deleted file mode 100644 index 9e26dfe..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/CFBC27D1 +++ /dev/null @@ -1 +0,0 @@ -{} \ No newline at end of file diff --git a/06_Shiny/.Rproj.user/94131CCE/sources/prop/INDEX b/06_Shiny/.Rproj.user/94131CCE/sources/prop/INDEX deleted file mode 100644 index c0e95cf..0000000 --- a/06_Shiny/.Rproj.user/94131CCE/sources/prop/INDEX +++ /dev/null @@ -1,10 +0,0 @@ -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fapp%2Fglobal.R="8BA97ED7" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fapp%2Fserver.R="CD86E44A" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fapp%2Fui.R="0E49E7B8" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2FR%2Fnetwork_multiRegion.R="9889CFFF" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2FR%2Fnetwork_singleRegion.R="1532276C" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2FR%2Fupset_plot.R="0DDD84F5" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2Fapp.R="0C9B8A0A" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2Fserver.R="8D825508" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2Fui.R="15E20D6C" -~%2FDocuments%2FownCloud%2FDexStim_RNAseq_Mouse%2Fscripts%2F06_Shiny%2Fupset_plot.R="CFBC27D1" diff --git a/06_Shiny/.Rproj.user/shared/notebooks/patch-chunk-names b/06_Shiny/.Rproj.user/shared/notebooks/patch-chunk-names deleted file mode 100644 index e69de29..0000000 diff --git a/06_Shiny/.Rproj.user/shared/notebooks/paths b/06_Shiny/.Rproj.user/shared/notebooks/paths deleted file mode 100644 index d9e2a9e..0000000 --- a/06_Shiny/.Rproj.user/shared/notebooks/paths +++ /dev/null @@ -1,5 +0,0 @@ -/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/R/network_singleRegion.R="2E73C7DD" -/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/app.R="1F978958" -/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/server.R="9E289B0A" -/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/ui.R="D511866B" -/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/06_Shiny/upset_plot.R="91BB5B50" diff --git a/06_Shiny/R/.DS_Store b/06_Shiny/R/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/06_Shiny/R/.DS_Store and /dev/null differ diff --git a/06_Shiny/R/network_multiRegion.R b/06_Shiny/R/network_multiRegion.R deleted file mode 100644 index 0a2e053..0000000 --- a/06_Shiny/R/network_multiRegion.R +++ /dev/null @@ -1,515 +0,0 @@ -# module UI function -networkMultiUI <- function(id, label = "Network Multiple Regions"){ - - ns <- NS(id) - - tagList( - - 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("Baseline", "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:", - choices = list("Ensembl", "Gene Symbol"), - selected = "Gene Symbol" - ), - checkboxGroupInput( - ns("tableContent"), - label = "Table content:", - 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(" at baseline, after 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(), - visNetworkOutput(ns("network_diff_multi"), height = "600px") - %>% withSpinner(color="#0dc5c1"), - 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("Baseline 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( - basepath, - "tables/coExpression_kimono/03_AnalysisFuncoup/", - "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) - }) - - - # #### TEST - # - # list_network <- list() - # - # for (reg in c("AMY", "CER", "PFC")){ - # file_path <- paste0( - # basepath, - # "tables/coExpression_kimono/03_AnalysisFuncoup/", - # "04_singleRegion_", - # reg, - # "_filtered_", - # "baseline", - # "Network.csv" - # ) - # # Read data - # data <- - # fread(file = file_path) %>% - # #dplyr::select(target, predictor, beta_mean.base, beta_mean.dex, z, p_adj) - # dplyr::select(target, predictor, beta_mean, beta_mean.dex, z, p_adj) - # - # cols <- colnames(data)[3:6] - # #cols <- colnames(data)[3] - # data <- data %>% - # 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")) - - - # 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(basepath, "tables/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( - basepath, - "tables/coExpression_kimono/03_AnalysisFuncoup/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({ - - network_table() %>% - datatable(filter = list(position = 'top')) - }, 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 = c, - title = data$title) - # df for edge legend - ledges <- data.frame(color = edge_colors, - label = 1:(ncol(data)-6), - font.align = "top") - } 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 = 1:(ncol(data)-6), - font.align = "top") - } - 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"), - shape = c( "ellipse"), color = c("orange", "darkred", "purple", "darkblue"), - title = "Informations", id = 1:4) - - - # return(list("edges" = relations, "nodes" = nodes, "ledges" = ledges, "lnodes" = lnodes)) - return(list("edges" = relations, "nodes" = nodes, "ledges" = ledges)) - - } - - - - # prepare network table for visualization - network_table <- reactive({ - - req(network_data()) - req(input_genes()) - - print(input$tableContent) - - - # 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) - - # 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")) - } - - } - - data - - }) - - - } - - ) -} \ No newline at end of file diff --git a/06_Shiny/R/network_singleRegion.R b/06_Shiny/R/network_singleRegion.R deleted file mode 100644 index af3b095..0000000 --- a/06_Shiny/R/network_singleRegion.R +++ /dev/null @@ -1,260 +0,0 @@ -# module UI function -networkSingleUI <- function(id, label = "Network Single Region"){ - - ns <- NS(id) - - tagList( - - 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("Baseline", "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:", - 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(" at baseline, after 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(), - visNetworkOutput(ns("network"), height = "600px") %>% - withSpinner(color="#0dc5c1"), - 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("Baseline 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( - basepath, - "tables/coExpression_kimono/03_AnalysisFuncoup/", - "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(basepath, "tables/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(basepath, "tables/coExpression_kimono/03_AnalysisFuncoup/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({ - - network_data() %>% - datatable(filter = list(position = 'top')) - }, 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/06_Shiny/R/upsetPlot_hub.R b/06_Shiny/R/upsetPlot_hub.R deleted file mode 100644 index a4ac022..0000000 --- a/06_Shiny/R/upsetPlot_hub.R +++ /dev/null @@ -1,448 +0,0 @@ -# module UI function -comparisonHubGenesUI <- function(id, label = "Upset Plot Hub Genes"){ - - ns <- NS(id) - - #tabPanel("Comparison Brain Regions", - tagList( - 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(basepath, "tables/coExpression_kimono/", - "03_AnalysisFuncoup/"), - 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(basepath, "tables/coExpression_kimono/", - "03_AnalysisFuncoup/"), - 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 = "black", width = 3), - marker = list(color = "black", - 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 = "black"), - 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 = "black", - 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/06_Shiny/R/upset_plot.R b/06_Shiny/R/upset_plot.R deleted file mode 100644 index a190690..0000000 --- a/06_Shiny/R/upset_plot.R +++ /dev/null @@ -1,454 +0,0 @@ -# module UI function -upsetPlotUI <- function(id, label = "Upset Plot"){ - - ns <- NS(id) - - #tabPanel("Comparison Brain Regions", - tagList( - 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"), - 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(basepath, "tables/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(basepath, "tables/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 - - }) - - - # 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 = "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 = "black", width = 3), - marker = list(color = "black", - 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 = "black"), - 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 = "black", - 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/06_Shiny/app.R b/06_Shiny/app.R deleted file mode 100644 index 84b346c..0000000 --- a/06_Shiny/app.R +++ /dev/null @@ -1,2 +0,0 @@ - -shinyApp(ui,server) \ No newline at end of file diff --git a/06_Shiny/global.R b/06_Shiny/global.R deleted file mode 100644 index 935eb3d..0000000 --- a/06_Shiny/global.R +++ /dev/null @@ -1,21 +0,0 @@ -library(shiny) -library(shinythemes) -library(markdown) -library(plotly) -library(DT) -library(data.table) -library(stringr) -library(dplyr) -library(ggplot2) -#library(igraph) -library(visNetwork) -library(shinyWidgets) -#library(UpSetR) -library(org.Mm.eg.db) -library(scales) -library(shinycssloaders) -library(RColorBrewer) - - -basepath <- "~/Documents/ownCloud/DexStim_RNAseq_Mouse/" -source(paste0(basepath, "scripts/06_Shiny/utilities_network.R")) \ No newline at end of file diff --git a/06_Shiny/server.R b/06_Shiny/server.R deleted file mode 100644 index 9988378..0000000 --- a/06_Shiny/server.R +++ /dev/null @@ -1,200 +0,0 @@ -source(paste0(basepath, "scripts/06_Shiny/utilities_network.R")) -#source(paste0(basepath, "scripts/06_Shiny/utilities_de.R")) - -# Define server logic required to generate and plot -server <- shinyServer(function(input, output) { - - ### DIFFERENTIAL EXPRESSION ---------------------------------- - - # DE data set - de_data <- reactive({ - # read DE table to create volcano plot - table <- fread(paste0(basepath, "tables/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(basepath, "tables/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("Baseline", "Treatment"), - values = c("#B0BFBB", "#46866E")) + - scale_x_discrete(breaks = c("CNTRL", "DEX"), - labels = c("Baseline", "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({ - de_data() %>% - # dplyr::rename("Ensembl_ID" = "V1") %>% - datatable(filter = list(position = 'top')) - }, 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(basepath, "tables/coExpression_kimono/", - "03_AnalysisFuncoup"), - 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"){ - histogram_network <- ggplot(overview_data(), aes(x=nodedegree)) + - geom_histogram() - } else { - histogram_network <- ggplot(overview_data(), aes(x=nodebetweenness)) + - geom_histogram() - } - ggplotly(histogram_network) - }) - - # Barplot network overview - output$barplot_network <- renderPlotly({ - - req(overview_data()) - - # Histogram of nodedegree/nodebetweenness - if (input$overview_metric == "nodedegrees"){ - barplot_network <- ggplot(overview_data(), aes(x=nodedegree)) - } else { - barplot_network <- ggplot(overview_data(), aes(x=nodebetweenness)) + - geom_histogram() - } - - barplot_network <- barplot_network + - geom_bar() - ggplotly(barplot_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/06_Shiny/ui.R b/06_Shiny/ui.R deleted file mode 100644 index f46823e..0000000 --- a/06_Shiny/ui.R +++ /dev/null @@ -1,239 +0,0 @@ - -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"), - fluidRow(column(tags$img(src="mousejavi_reversebrain.png", - width="100%", height= "100%"), - #width="450px",height="320px"), - width=floor(0.4*12)), - column( - p(strong("Abstract."), "Network analysis can identify the molecular connectivity - that is underpinning function at baseline, 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 baseline and treatment levels using KiMONo (Ogris et al.) and at - the differential level using DiffGRN (Kim et al.). DiffGRN uses the regression coefficients - of the baseline and treatment networks to calculate differential connections, thus providing - us with the actual functional couplings of the genes after the stimulus that differ from the - baseline ones.",br(), - "As a stimulus we used dexamethasone. Dexamethasone is 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. - Glucocorticoid receptors chromatin binding is one of the main results of stress-axis activation, - which in turn is an important component of the biology of stress-related psychiatric disorders. - Thus, dexamethasone mimics some aspects of the altered transcriptomic landscape that is associated - with stress-related psychiatric disorders. 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 baseline, dexamethasone-treatment and differential levels. 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 baseline and - differential levels.", - style="text-align:justify;color:black;background-color:#2980B9;padding:15px;border-radius:10px"), - br(), - p(strong("Data and code availability. Reference to paper."), "Sed ut perspiciatis unde omnis iste - natus error sit voluptatem accusantium doloremque laudantium, totam rem aperiam, eaque ipsa quae - ab illo inventore veritatis et quasi architecto beatae vitae dicta sunt explicabo. Nemo enim - ipsam voluptatem quia voluptas sit aspernatur aut odit aut fugit, sed quia consequuntur magni - dolores eos qui ratione voluptatem sequi nesciunt. Neque porro quisquam est, qui dolorem ipsum - quia dolor sit amet, consectetur, adipisci velit, sed quia non numquam eius modi tempora incidunt - ut labore et dolore magnam aliquam quaerat voluptatem. Ut enim ad minima veniam, quis nostrum - exercitationem ullam corporis suscipit laboriosam, nisi ut aliquid ex ea commodi consequatur? - Quis autem vel eum iure reprehenderit qui in ea voluptate velit esse quam nihil molestiae - consequatur, vel illum qui dolorem eum fugiat quo voluptas nulla pariatur?", - 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.png", width="250px", height="60px"), - br(), - br(), - p("For more information please visit the website of", em("the MPI of Psychiatry"), - br(), - a(href="https://www.psych.mpg.de/", "Here", target="_blank"), - style="text-align:center;color:black"), - width=ceiling(0.15*12) - ))), - - navbarMenu("Diff. Gene Expression", - 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 before and after Dexamethasone treatment. 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")), - DT::dataTableOutput("de_table") - ) - ) - ) - ), - tabPanel("Comparison Brain Regions", - # UI defined in upsetPlot module - upsetPlotUI("upsetDE") - ) - ), - - navbarMenu( - "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", - 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("Baseline" = "baseline", - "Differential" = "differential"), - 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(), - splitLayout(cellWidths = c("50%", "50%"), - h4( - p("Baseline network", style = "color:black;text-align:center") - ), - h4( - p("Differential network", style = "color:black;text-align:center") - )), - # plotlyOutput("histogram_network") - splitLayout( - cellWidths = c("50%", "50%"), - plotlyOutput("histogram_network"), - plotlyOutput("barplot_network") - ) - # DT::dataTableOutput("network_table") - ) - ) - ), - - 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("About", - # 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") - # ) - # ) - # ) - ) - ) -) -) diff --git a/06_Shiny/utilities_network.R b/06_Shiny/utilities_network.R deleted file mode 100644 index 0d449fc..0000000 --- a/06_Shiny/utilities_network.R +++ /dev/null @@ -1,178 +0,0 @@ -### 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/06_Shiny/v2_dashboard/server.R b/06_Shiny/v2_dashboard/server.R deleted file mode 100644 index e69de29..0000000 diff --git a/06_Shiny/v2_dashboard/ui.R b/06_Shiny/v2_dashboard/ui.R deleted file mode 100644 index ec90b2f..0000000 --- a/06_Shiny/v2_dashboard/ui.R +++ /dev/null @@ -1,91 +0,0 @@ -library(semantic.dashboard) -library(shiny) -library(markdown) -library(plotly) - -ui <- dashboardPage(dashboardHeader("Explore!"), - ## Sidebar content - dashboardSidebar(sidebarMenu( - menuItem(tabName = "home", text = "Home", icon = icon("home")), - menuItem(tabName = "de", text = "Diff. Gene Expression", icon = "heart")) - ), - dashboardBody( - tabItems( - # First tab content - tabItem(tabName = "de", - fluidRow( - selectInput("region", "Choose a brain region:", - choices = c("AMY", "CER", "dDG", - "dCA1", "PFC", - "vDG", "vCA1")), - plotlyOutput("volcano_output"), - plotlyOutput("exp_plot") - )) - ) - )) - - -library(data.table) -library(stringr) -library(ggplot2) - - -# Define server logic required to generate and plot -server <- function(input, output) { - - output$volcano_plot <- renderPlotly({ - - table <- fread(paste0("/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - input$region,"_deseq2_Dex_1_vs_0_lfcShrink.txt")) - # table <- fread(paste0("/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - # "AMY","_deseq2_Dex_1_vs_0_lfcShrink.txt")) - table$sig <- as.factor(ifelse(table$padj <= 0.1, "p_adj <= 0.1", "p_adj > 0.1")) - - # plot volcano plot - volcano_plot <- plot_ly(data = table, x = ~log2FoldChange, y = ~-log10(padj), - color = ~sig, colors = "Set1", - text = ~paste0("Gene: ", V1)) %>% - layout(title = paste0("Volcano Plot ", input$region)) - }) - - output$exp_plot <- renderPlotly({ - d <- event_data("plotly_click") - print(d) - if (is.null(d)) return(NULL) - - table <- fread(paste0("/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - input$region,"_deseq2_Dex_1_vs_0_lfcShrink.txt")) - ensembl_id <- table$V1[d$pointNumber + 1] - - exp_table <- fread(paste0("/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - input$region,"_deseq2_expression_vsd.txt")) - # exp_table <- fread(paste0("/Users/nathalie_gerstner/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - # "AMY","_deseq2_expression_vsd.txt")) - - exp_gene <- transpose(exp_table[V1 == ensembl_id,2:ncol(exp_table)], keep.names = "col") - exp_gene$col <- str_replace(exp_gene$col, ".*\\_","") - exp_gene$mouse_id <- str_extract(exp_gene$col, pattern = "[0-9]+") - exp_gene$col <- as.factor(str_replace(exp_gene$col, "[0-9]+","")) - - exp_plot <- ggplot(exp_gene, aes(x = col, y = V1)) + - geom_boxplot() + - geom_jitter(color="black", size=0.4, alpha=0.9) + - theme_light() + - theme( - legend.position="none", - plot.title = element_text(size=11) - ) + - ggtitle("A boxplot with jitter") + - xlab("") - ggplotly(exp_plot) - - # gSel <- tg %>% filter(dose %in% d$y) %>% group_by(supp) %>% mutate(newLen=floor(len)) %>% - # ggplot(aes(x=supp, fill=as.factor(newLen))) + geom_bar() - # ggplotly(gSel) - - - }) -} - -shinyApp(ui, server) - \ No newline at end of file diff --git a/06_Shiny/www/DiffBrainNet_logo.png b/06_Shiny/www/DiffBrainNet_logo.png deleted file mode 100644 index 607caed..0000000 Binary files a/06_Shiny/www/DiffBrainNet_logo.png and /dev/null differ diff --git a/06_Shiny/www/DiffNetworks.png b/06_Shiny/www/DiffNetworks.png deleted file mode 100644 index 5da4451..0000000 Binary files a/06_Shiny/www/DiffNetworks.png and /dev/null differ diff --git a/06_Shiny/www/DiffNetworks.pptx b/06_Shiny/www/DiffNetworks.pptx deleted file mode 100644 index 9c154e3..0000000 Binary files a/06_Shiny/www/DiffNetworks.pptx and /dev/null differ diff --git a/06_Shiny/www/TermsOfUse_v1.docx b/06_Shiny/www/TermsOfUse_v1.docx deleted file mode 100644 index 957f540..0000000 Binary files a/06_Shiny/www/TermsOfUse_v1.docx and /dev/null differ diff --git a/06_Shiny/www/WhatsApp Image 2021-02-09 at 22.17.04.jpeg b/06_Shiny/www/WhatsApp Image 2021-02-09 at 22.17.04.jpeg deleted file mode 100644 index e487ff9..0000000 Binary files a/06_Shiny/www/WhatsApp Image 2021-02-09 at 22.17.04.jpeg and /dev/null differ diff --git a/06_Shiny/www/mouse_design.jpg b/06_Shiny/www/mouse_design.jpg deleted file mode 100644 index 03e0719..0000000 Binary files a/06_Shiny/www/mouse_design.jpg and /dev/null differ diff --git a/06_Shiny/www/mouse_design.png b/06_Shiny/www/mouse_design.png deleted file mode 100644 index 79bee00..0000000 Binary files a/06_Shiny/www/mouse_design.png and /dev/null differ diff --git a/06_Shiny/www/mousejavi b/06_Shiny/www/mousejavi deleted file mode 100644 index 78c90f1..0000000 Binary files a/06_Shiny/www/mousejavi and /dev/null differ diff --git a/06_Shiny/www/mousejavi.png b/06_Shiny/www/mousejavi.png deleted file mode 100644 index 78c90f1..0000000 Binary files a/06_Shiny/www/mousejavi.png and /dev/null differ diff --git a/06_Shiny/www/mousejavi_reversebrain b/06_Shiny/www/mousejavi_reversebrain deleted file mode 100644 index 166e1ec..0000000 Binary files a/06_Shiny/www/mousejavi_reversebrain and /dev/null differ diff --git a/06_Shiny/www/mousejavi_reversebrain.png b/06_Shiny/www/mousejavi_reversebrain.png deleted file mode 100644 index f0384e2..0000000 Binary files a/06_Shiny/www/mousejavi_reversebrain.png and /dev/null differ diff --git a/06_Shiny/www/mpilogo.png b/06_Shiny/www/mpilogo.png deleted file mode 100644 index 74c563d..0000000 Binary files a/06_Shiny/www/mpilogo.png and /dev/null differ diff --git a/06_Shiny/www/~$rmsOfUse_v1.docx b/06_Shiny/www/~$rmsOfUse_v1.docx deleted file mode 100644 index c75703e..0000000 Binary files a/06_Shiny/www/~$rmsOfUse_v1.docx and /dev/null differ diff --git a/07_PlotsManuscript/.DS_Store b/07_PlotsManuscript/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/07_PlotsManuscript/.DS_Store and /dev/null differ diff --git a/07_PlotsManuscript/01_FigureII_C.R b/07_PlotsManuscript/01_FigureII_C.R deleted file mode 100644 index 69f79c8..0000000 --- a/07_PlotsManuscript/01_FigureII_C.R +++ /dev/null @@ -1,89 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 21.11.2021 -## Author: Nathalie -################################################## -# GO Enrichment plot for manuscript -# Unique DE and hub genes in PFC - -# set working directory to source file location -setwd("~/Documents/ownCloud/DexStim_RNAseq_Mouse/scripts/07_PlotsManuscript") - -library(readxl) -library(dplyr) -library(stringr) -library(ggplot2) - - -# 1. GO enrichment unique DE and hub genes in PFC - -data_hub <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure II_PFC/FUMA_diff unique hubgenes_gene2func67664/PFC diff unique hubgenes_FUMA.xlsx", - sheet = "GO bp") - -data_DE <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure II_PFC/DE PFC/FUMA_gene2func66314/PFC unique DE genes _FUMA.xlsx", - sheet = "GO bp") - - -# data for left panel with top 14 terms for DE -top_DE <- data_DE %>% - dplyr::filter(N_overlap_A > 24) -top_DE <- top_DE[order(top_DE$`my adj p BH`, decreasing = FALSE),][1:14,] - -top_DE_hub <- data_hub %>% - filter(GeneSet %in% top_DE$GeneSet) - -top_DE_comb <- bind_rows("de" = top_DE, "hub" = top_DE_hub, .id = "group") - - -# data for right panel with top 14 terms for hubs -top_hub <- data_hub %>% - dplyr::filter(N_overlap_A > 2) -top_hub <- top_hub[order(top_hub$`my adj p BH`, decreasing = FALSE),][1:14,] - -top_hub_DE <- data_DE %>% - filter(GeneSet %in% top_hub$GeneSet) - -missing_term <- setdiff(top_hub$GeneSet, top_hub_DE$GeneSet) -add_entry <- data.frame("GeneSet" = missing_term) -top_hub_DE <- bind_rows(top_hub_DE, add_entry) - -top_hub_comb <- bind_rows("hub" = top_hub, "de" = top_hub_DE, .id = "group") - - -# combine everything into one df for plotting -top_data <- bind_rows("de" = top_DE_comb, "hub" = top_hub_comb, .id = "panel") -# remove the "GO_" from GO terms -top_data$GeneSet <- str_replace_all(top_data$GeneSet, '^GO_', ' ') -# remove the "_" from GO terms -top_data$GeneSet <- str_replace_all(top_data$GeneSet, "_", " ") -# only capitalize first letter per word -top_data$GeneSet <- str_to_title(top_data$GeneSet) -# of and to should start with lower case -top_data$GeneSet <- str_replace_all(top_data$GeneSet, "Of", "of") -top_data$GeneSet <- str_replace_all(top_data$GeneSet, "To", "to") - -# labeller function for facet titles -facet_names <- c( - 'de'="Top 14 GO terms for DE genes", - 'hub'="Top 14 GO terms for hub genes" -) - -# barplot -g <- ggplot(top_data, aes(x = GeneSet, y = `[-log10 fdr`, fill = group)) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity") + - - scale_colour_manual(values = c("grey", "black"), guide = FALSE) + - scale_fill_manual(name = "Geneset", - values = c("orange", "darkred"), - labels = c("DE genes", "hub genes")) + - coord_flip() + - scale_x_discrete(limits = rev) + - xlab("GOterm") + - ylab("-log10(FDR)") + - ggtitle("GO terms enriched for unique DE and hub genes in PFC") + - facet_wrap(~panel, scales="free", labeller = as_labeller(facet_names)) + - theme_light() + - theme(text = element_text(size= 14)) -print(g) -ggsave("01_FigureII_C.pdf", g, width = 14, height = 8) - diff --git a/07_PlotsManuscript/02_FigureII_F.R b/07_PlotsManuscript/02_FigureII_F.R deleted file mode 100644 index 6b2283e..0000000 --- a/07_PlotsManuscript/02_FigureII_F.R +++ /dev/null @@ -1,83 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 22.11.2021 -## Author: Nathalie -################################################## -# Pathway enrichment plot for manuscript -# Abcd1 and neighbourhood in diff network of PFC - -library(readxl) - -# 1. Pathway enrichment for Abcd1 and neighbours - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure II_PFC/Abcd1 gene neighbourhood_FUMA_gene2func67210/Abcd1_diffnetwork_FUMA.xlsx", - sheet = "KEGG AND REACTOME") - -data <- arrange(data, desc(`[-log10FDR]`)) -data <- data[1:20,] - - -# remove the "_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, "_", " ") -# only capitalize first letter per word -data$GeneSet <- str_to_title(data$GeneSet) -# of and to should start with lower case -data$GeneSet <- str_replace_all(data$GeneSet, "Of", "of") -data$GeneSet <- str_replace_all(data$GeneSet, "To", "to") -data$GeneSet <- str_replace_all(data$GeneSet, " In ", " in ") -data$GeneSet <- str_replace_all(data$GeneSet, " Into ", " into ") -# Ii from "Polymerase II" back to upper case -data$GeneSet <- str_replace_all(data$GeneSet, "Reactome", "Reactome:") -data$GeneSet <- str_replace_all(data$GeneSet, "Kegg", "KEGG:") -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - - - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-log10FDR]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.9,1,1)) - -ggexport(comb, filename = "02_FigureII_F.pdf", width = 14, height = 8) diff --git a/07_PlotsManuscript/03_FigureIII_B.R b/07_PlotsManuscript/03_FigureIII_B.R deleted file mode 100644 index 977faf3..0000000 --- a/07_PlotsManuscript/03_FigureIII_B.R +++ /dev/null @@ -1,83 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 29.11.2021 -## Author: Nathalie -################################################## -# GO enrichment plot for manuscript -# vCA1 unique DE genes enrichments - -library(readxl) -library(dplyr) -library(stringr) -library(ggpubr) - -# 1. GO enrichment for vCA unique DE genes - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure III_vCA1/DE vCA1/FUMA_gene2func67310 (2)/vCA1 Unique DE genes_FUMA enrichments.xlsx", - sheet = "GO bp") - -data <- arrange(data, desc(`[-log10 of FDR`)) -data <- data[1:20,] - - -# remove the "GO_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, '^GO_', ' ') -# remove the "_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, "_", " ") -# only capitalize first letter per word -data$GeneSet <- str_to_title(data$GeneSet) -# of and to should start with lower case -data$GeneSet <- str_replace_all(data$GeneSet, "Of", "of") -data$GeneSet <- str_replace_all(data$GeneSet, "To", "to") -data$GeneSet <- str_replace_all(data$GeneSet, " In ", " in ") -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-log10 of FDR`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(2.5,1,1)) - -ggexport(comb, filename = "03_FigureIII_B.pdf", width = 14, height = 8) diff --git a/07_PlotsManuscript/04_FigureIII_C.R b/07_PlotsManuscript/04_FigureIII_C.R deleted file mode 100644 index ba9e55f..0000000 --- a/07_PlotsManuscript/04_FigureIII_C.R +++ /dev/null @@ -1,84 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 29.11.2021 -## Author: Nathalie -################################################## -# GO enrichment plot for manuscript -# vCA1 unique DE genes with neighbours enrichments - - -library(readxl) -library(dplyr) -library(stringr) -library(ggpubr) - -# 1. GO enrichment for vCA unique DE genes and their diff neighbours - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure III_vCA1/network vCA1/Unique DEs and their diff neighbours_FUMA_gene2func67313/vCA1 Unique DE genes+diff neighbours_FUMA enrichments.xlsx", - sheet = "GO bp") - -data <- arrange(data, desc(`[-LOG10 of FDR`)) -data <- data[1:20,] - -# remove the "GO_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, '^GO_', ' ') -# remove the "_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, "_", " ") -# only capitalize first letter per word -data$GeneSet <- str_to_title(data$GeneSet) -# of and to should start with lower case -data$GeneSet <- str_replace_all(data$GeneSet, "Of", "of") -data$GeneSet <- str_replace_all(data$GeneSet, "To", "to") -data$GeneSet <- str_replace_all(data$GeneSet, " In ", " in ") -data$GeneSet <- str_replace_all(data$GeneSet, " Into ", " into ") -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-LOG10 of FDR`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.8,1,1)) - -ggexport(comb, filename = "03_FigureIII_C.pdf", width = 12, height = 8) diff --git a/07_PlotsManuscript/05_FigureIV_Bbase.R b/07_PlotsManuscript/05_FigureIV_Bbase.R deleted file mode 100644 index f996410..0000000 --- a/07_PlotsManuscript/05_FigureIV_Bbase.R +++ /dev/null @@ -1,90 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 29.11.2021 -## Author: Nathalie -################################################## -# GO enrichment plot for manuscript -# Tcf4 baseline network in PFC - -library(readxl) -library(dplyr) -library(stringr) -library(ggpubr) - -# 1. GO enrichment for vCA unique DE genes and their diff neighbours - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure IV_GWAS/TCF4/FUMA_baselin Tcf4PFCnetwork_gene2func68958/FUMA Tcf4BaselinePFCnetwork.xlsx", - sheet = "GO bp") - -data <- arrange(data, desc(`[-log10]`)) -data <- data[1:25,] - -# remove the "GO_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, '^GO_', ' ') -# remove the "_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, "_", " ") -# only capitalize first letter per word -data$GeneSet <- str_to_title(data$GeneSet) -# of and to should start with lower case -data$GeneSet <- str_replace_all(data$GeneSet, "Of", "of") -data$GeneSet <- str_replace_all(data$GeneSet, "To", "to") -data$GeneSet <- str_replace_all(data$GeneSet, " In ", " in ") -data$GeneSet <- str_replace_all(data$GeneSet, " Into ", " into ") -# make some terms shorter with abbreviations -data$GeneSet <- str_replace_all(data$GeneSet, "Positive Regulation", "Pos. Reg.") -data$GeneSet <- str_replace_all(data$GeneSet, "Negative Regulation", "Neg. Reg.") -# Ii from "Polymerase II" back to upper case -data$GeneSet <- str_replace_all(data$GeneSet, "Ii", "II") -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-log10]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.9,1,1)) - -ggexport(comb, filename = "05_FigureIV_Bbase.pdf", width = 14, height = 10) - - diff --git a/07_PlotsManuscript/05_FigureIV_Bdiff.R b/07_PlotsManuscript/05_FigureIV_Bdiff.R deleted file mode 100644 index ddd3f64..0000000 --- a/07_PlotsManuscript/05_FigureIV_Bdiff.R +++ /dev/null @@ -1,88 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 29.11.2021 -## Author: Nathalie -################################################## -# GO enrichment plot for manuscript -# Tcf4 differential network in PFC - -library(readxl) -library(dplyr) -library(stringr) -library(ggpubr) - -# 1. GO enrichment for vCA unique DE genes and their diff neighbours - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure IV_GWAS/TCF4/FUMA_diffTcf4PFCnetwork_gene2func68788/FUMA_diffPFCTcf4.xlsx", - sheet = "GO bp") - -data <- arrange(data, desc(`[-log10]`)) -data <- data[1:25,] - -# remove the "GO_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, '^GO_', ' ') -# remove the "_" from GO terms -data$GeneSet <- str_replace_all(data$GeneSet, "_", " ") -# only capitalize first letter per word -data$GeneSet <- str_to_title(data$GeneSet) -# of and to should start with lower case -data$GeneSet <- str_replace_all(data$GeneSet, "Of", "of") -data$GeneSet <- str_replace_all(data$GeneSet, "To", "to") -data$GeneSet <- str_replace_all(data$GeneSet, " In ", " in ") -data$GeneSet <- str_replace_all(data$GeneSet, " Into ", " into ") -# make some terms shorter with abbreviations -data$GeneSet <- str_replace_all(data$GeneSet, "Positive Regulation", "Pos. Reg.") -data$GeneSet <- str_replace_all(data$GeneSet, "Negative Regulation", "Neg. Reg.") -# Ii from "Polymerase II" back to upper case -data$GeneSet <- str_replace_all(data$GeneSet, "Ii", "II") -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-log10]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.9,1,1)) - -ggexport(comb, filename = "05_FigureIV_Bdiff.pdf", width = 14, height = 10) \ No newline at end of file diff --git a/07_PlotsManuscript/06_FigureSIII_A.R b/07_PlotsManuscript/06_FigureSIII_A.R deleted file mode 100644 index 01122cd..0000000 --- a/07_PlotsManuscript/06_FigureSIII_A.R +++ /dev/null @@ -1,72 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 29.11.2021 -## Author: Nathalie -################################################## -# GWAS enrichment plot for manuscript -# vCA1 - -library(readxl) -library(dplyr) -library(ggplot2) -library(ggpubr) - -# 1. Pathway enrichment for Abcd1 and neighbours - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure III_vCA1/network vCA1/Unique DEs and their diff neighbours_FUMA_gene2func67313/vCA1 Unique DE genes+diff neighbours_FUMA enrichments.xlsx", - sheet = "GWAS") - -data <- arrange(data, desc(`[-log10 of FDR`)) -data <- data[data$`my adj p BH` <= 0.05,] -data <- data[!is.na(data$`my adj p BH`),] - -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-log10 of FDR`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.8,1,1)) - -ggexport(comb, filename = "06_FigureSIII_A.pdf", width = 14, height = 8) diff --git a/07_PlotsManuscript/07_FigureSIV.R b/07_PlotsManuscript/07_FigureSIV.R deleted file mode 100644 index 9afb5a5..0000000 --- a/07_PlotsManuscript/07_FigureSIV.R +++ /dev/null @@ -1,120 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 29.11.2021 -## Author: Nathalie -################################################## -# Pathway enrichment plot for manuscript -# Tcf4 Differential Expression - -library(stringr) -library(ggplot2) - -regions <- c("AMY", "CER", "dCA1", "dDG", "PFC", "PVN", "vCA1", "vDG") - -ensembl_tcf4 <- "ENSMUSG00000053477" - -# Panel A: Expression values -df <- data.frame("region" = character(), - "treatment" = character(), - "sample" = character(), - "expression" = numeric()) - -for (reg in regions){ - - # read vsd normalized data --> better raw data? - data_exp <- read.table(paste0("~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - reg, "_deseq2_expression_vsd.txt"), - header = TRUE, sep = "\t") - # subset Tcf4 row - tcf4 <- data_exp[ensembl_tcf4,] - - # get column indices for cntrl samples - indices_cntrl <- str_detect(colnames(data_exp), "CNTRL") - - # df for cntrl samples - exp_cntrl <- data.frame("expression" = t(tcf4[,indices_cntrl]), - "sample" = rownames(t(tcf4[,indices_cntrl]))) %>% - rename(expression = ENSMUSG00000053477) %>% - mutate("treatment" = "CNTRL", "region" = reg) - df <- rbind(df, exp_cntrl) - - # df for dex samples - exp_dex <- data.frame("expression" = t(tcf4[,!indices_cntrl]), - "sample" = rownames(t(tcf4[,!indices_cntrl]))) %>% - rename(expression = ENSMUSG00000053477) %>% - mutate("treatment" = "DEX", "region" = reg) - df <- rbind(df, exp_dex) - -} - -df$treatment <- factor(df$treatment) -df$region <- factor(df$region) - -ggplot(df, aes(x = region, y = expression, fill = treatment)) + - geom_boxplot() + - scale_fill_manual("", - breaks = c("CNTRL", "DEX"), - labels = c("Baseline", "Treatment"), - values = c("#B0BFBB", "#46866E")) + - theme_light() + - theme(text = element_text(size= 14)) + - xlab("Brain Region") + - ylab("Norm. Expression Level") - -ggsave(filename = "07_FigureSIV_A.pdf", width = 12, height = 8) - - -# Panel B: Fold Change in AMY, vDG and dDG - -list_tcf4 <- list() - -for (reg in c("AMY", "dDG", "vDG")){ -# for (reg in regions){ - - # read DE results - data_exp <- read.table(paste0("~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - reg, "_deseq2_Dex_1_vs_0_lfcShrink.txt"), - header = TRUE, sep = "\t") - tcf4 <- as.data.frame(data_exp[data_exp$Ensembl_ID == ensembl_tcf4,]) - print(tcf4) - - list_tcf4[[reg]] <- tcf4 - -} - -df_tcf4 <- bind_rows(list_tcf4, .id = "region") - - -# bar plot Fold Change -fc <- ggplot(df_tcf4, aes(x = region, y = log2FoldChange) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "darkgrey") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("log2(FoldChange)") + - xlab("Brain Region") - -# barplot fdr p-value -fdr <- ggplot(df_tcf4, aes(x = region, y = -log10(padj)) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.1),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - xlab("") - -# combined barplot -comb <- ggarrange(fc, - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.3,1)) - -# ggexport(comb, filename = "07_FigureSIV_B_allRegions.pdf", width = 12, height = 8) -ggexport(comb, filename = "07_FigureSIV_B.pdf", width = 12, height = 8) - diff --git a/07_PlotsManuscript/08_FigureIII_A.R b/07_PlotsManuscript/08_FigureIII_A.R deleted file mode 100644 index f71384a..0000000 --- a/07_PlotsManuscript/08_FigureIII_A.R +++ /dev/null @@ -1,46 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# Numbers/Percentages of DE and hub genes in vCA1 -# including also number of unique ones - -library(ggplot2) - -# correct numbers? --> DE percentage different than before -DE_genes <- 465 -hub_genes <- 260 - -unique_DE <- 25 -unique_hub <- 24 - -# data frame for plotting -df <- data.frame("type" = c("DE", "DE", "hub", "hub"), - "unique" = c(TRUE, FALSE, TRUE, FALSE), - "number" = c(unique_DE, DE_genes - unique_DE, - unique_hub, hub_genes - unique_hub), - "text" = c(paste0(round(unique_DE/DE_genes*100,1),"%"), "", - paste0(round(unique_hub/hub_genes*100,1), "%"), "")) - -# barplot -ggplot(df, aes(x = type, y = number, fill = type)) + - geom_bar(position = "stack", stat="identity", aes(alpha = unique)) + - geom_text(aes(y = number, label = text), vjust = 1.5) + - scale_alpha_manual("", - breaks = c(FALSE, TRUE), - values = c(1.0,0.5), - labels = c("Shared genes", "Unique genes")) + - xlab("") + - ylab("# DE and hub genes in vCA1") + - scale_fill_manual("", - breaks = c("DE", "hub"), - labels = c("DE genes", "Hub genes"), - values = c("orange", "darkred")) + - scale_x_discrete(breaks = c("DE", "hub"), - labels = c("DE genes", "Hub genes")) + - theme_light() + - theme(text = element_text(size= 14)) + - guides(fill = "none") - -ggsave(filename = "08_FigureIII_A.pdf", width = 8, height = 6) diff --git a/07_PlotsManuscript/09_FigureII_DandE.R b/07_PlotsManuscript/09_FigureII_DandE.R deleted file mode 100644 index 2d195b9..0000000 --- a/07_PlotsManuscript/09_FigureII_DandE.R +++ /dev/null @@ -1,104 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# Expression level plots Syn3 and Abcd1 in PFC - -library(stringr) -library(ggplot2) - -reg <- "PFC" - -ensembl_syn3 <- "ENSMUSG00000059602" -ensembl_abcd1 <- "ENSMUSG00000031378" - -# read vsd normalized data -data_exp <- read.table(paste0("~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/02_", - reg, "_deseq2_expression_vsd.txt"), - header = TRUE, sep = "\t") - - -### Syn3 ---------------- - -# subset Syn3 row -syn3 <- data_exp[ensembl_syn3,] - -# get column indices for cntrl samples -indices_cntrl <- str_detect(colnames(data_exp), "CNTRL") - -# df for cntrl samples -exp_cntrl <- data.frame("expression" = t(syn3[,indices_cntrl]), - "sample" = rownames(t(syn3[,indices_cntrl]))) %>% - rename(expression = all_of(ensembl_syn3)) %>% - mutate("treatment" = "CNTRL") - -# df for dex samples -exp_dex <- data.frame("expression" = t(syn3[,!indices_cntrl]), - "sample" = rownames(t(syn3[,!indices_cntrl]))) %>% - rename(expression = all_of(ensembl_syn3)) %>% - mutate("treatment" = "DEX") - -# combine cntrl and dex df -df <- rbind(exp_cntrl, exp_dex) -df$treatment <- factor(df$treatment) - -ggplot(df, aes(x = treatment, y = expression, fill = treatment)) + - geom_boxplot() + - geom_jitter() + - scale_fill_manual("", - breaks = c("CNTRL", "DEX"), - labels = c("Baseline", "Treatment"), - values = c("#B0BFBB", "#46866E")) + - scale_x_discrete(breaks = c("CNTRL", "DEX"), - labels = c("Baseline", "Treatment")) + - theme_light() + - theme(text = element_text(size= 16)) + - xlab("") + - ylab("Norm. Expression Level") + - guides(fill = "none") - -ggsave(filename = "09_FigureII_D.pdf", width = 7, height = 6) - - - -### Abcd1 ---------------- - -# subset Abcd1 row -abcd1 <- data_exp[ensembl_abcd1,] - -# get column indices for cntrl samples -indices_cntrl <- str_detect(colnames(data_exp), "CNTRL") - -# df for cntrl samples -exp_cntrl <- data.frame("expression" = t(abcd1[,indices_cntrl]), - "sample" = rownames(t(abcd1[,indices_cntrl]))) %>% - rename(expression = all_of(ensembl_abcd1)) %>% - mutate("treatment" = "CNTRL") - -# df for dex samples -exp_dex <- data.frame("expression" = t(abcd1[,!indices_cntrl]), - "sample" = rownames(t(abcd1[,!indices_cntrl]))) %>% - rename(expression = all_of(ensembl_abcd1)) %>% - mutate("treatment" = "DEX") - -# combine cntrl and dex df -df <- rbind(exp_cntrl, exp_dex) -df$treatment <- factor(df$treatment) - -ggplot(df, aes(x = treatment, y = expression, fill = treatment)) + - geom_boxplot() + - geom_jitter() + - scale_fill_manual("", - breaks = c("CNTRL", "DEX"), - labels = c("Baseline", "Treatment"), - values = c("#B0BFBB", "#46866E")) + - scale_x_discrete(breaks = c("CNTRL", "DEX"), - labels = c("Baseline", "Treatment")) + - theme_light() + - theme(text = element_text(size= 16)) + - xlab("") + - ylab("Norm. Expression Level") + - guides(fill = "none") - -ggsave(filename = "09_FigureII_E.pdf", width = 7, height = 6) diff --git a/07_PlotsManuscript/10_FigureSII_A.R b/07_PlotsManuscript/10_FigureSII_A.R deleted file mode 100644 index 9430e6f..0000000 --- a/07_PlotsManuscript/10_FigureSII_A.R +++ /dev/null @@ -1,90 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# Number of DE genes and unique percentage - -library(data.table) -library(ggplot2) -library(dplyr) -library(tidyverse) - -regions <- - c("AMY", "PFC", "PVN", "CER", "vDG", "dDG", "vCA1", "dCA1") - - -# 1. Read DE tables from all regions ---------- - -list_reg_sig <- list() -list_genes_sig <- list() - -for (reg in regions) { - res <- - fread( - file = paste0( - "~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables", - "/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 <= 0.1, ] - # res_sig <- res[res$log2FoldChange >= 1] - list_reg_sig[[reg]] <- res_sig - list_genes_sig[[reg]] <- rownames(res_sig) -} - -# 2. Concatenate DE tables ----------------- - -data <- bind_rows(list_reg_sig, .id = "region") %>% - group_by(Ensembl_ID) %>% - summarise(region = list(region)) - -data_unique <- data %>% - mutate(nr_regions = lengths(region)) %>% - mutate(unique = (nr_regions == 1)) %>% - unnest(cols = c(region)) %>% - mutate("combined_id" = paste0(region, "-", Ensembl_ID)) - -data_barplot <- data_unique %>% - group_by(region, unique) %>% - count() %>% - group_by(region) %>% - mutate(sum = sum(n)) - - -# 3. Stacked barplot ------------------------- - -ggplot(data_barplot, aes(x = region, y = n, alpha = unique)) + - geom_bar(position = "stack", stat = "identity", fill = "orange") + - scale_alpha_manual( - name = "", - labels = c("DE in multiple regions", "DE unique"), - values = c(1.0, 0.5) - ) + - xlab("Brain region") + - ylab("# DE genes") + - theme_light() + - theme( - axis.title.x = element_text(size = 15), - axis.title.y = element_text(size = 15), - axis.text.x = element_text(size = 12), - axis.text.y = element_text(size = 12), - legend.text = element_text(size = 12), - # legend.title = element_blank(), - legend.position = "top" - ) + - geom_text(aes(label = paste0(round((n / sum) * 100, digits = 1 - ), "%")), - position = position_stack(vjust = 0.5), - size = 4, - show.legend = FALSE) -ggsave( - "10_FigureSII_A.pdf", - width = 8, - height = 6 -) diff --git a/07_PlotsManuscript/11_FigureSII_B.R b/07_PlotsManuscript/11_FigureSII_B.R deleted file mode 100644 index 3858189..0000000 --- a/07_PlotsManuscript/11_FigureSII_B.R +++ /dev/null @@ -1,92 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# Number of hub genes and unique percentage - -library(data.table) -library(ggplot2) -library(dplyr) -library(tidyverse) - -regions <- - c("AMY", "PFC", "PVN", "CER", "vDG", "dDG", "vCA1", "dCA1") - - -# 1. Read DE tables from all regions ---------- - -list_reg_sig <- list() -list_genes_sig <- list() - -for (reg in regions) { - res <- - fread( - file = paste0( - "~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/coExpression_kimono/03_AnalysisFuncoup/", - "/04_", - reg, - "_funcoup_differential_nodebetweennessNorm_betacutoff0.01.csv" - ), - sep = "," - ) - na_indices <- which(is.na(res$nodebetweenness_norm)) - res$padj[na_indices] <- 0 - res_sig <- res[res$nodebetweenness_norm >= 1.0, ] - list_reg_sig[[reg]] <- res_sig - list_genes_sig[[reg]] <- rownames(res_sig) -} - - -# 2. Concatenate hub tables ----------------- - -data <- bind_rows(list_reg_sig, .id = "region") %>% - group_by(ensembl_id) %>% - summarise(region = list(region)) - -data_unique <- data %>% - mutate(nr_regions = lengths(region)) %>% - mutate(unique = (nr_regions == 1)) %>% - unnest(cols = c(region)) %>% - mutate("combined_id" = paste0(region, "-", ensembl_id)) - -data_barplot <- data_unique %>% - group_by(region, unique) %>% - count() %>% - group_by(region) %>% - mutate(sum = sum(n)) - - -# 3. Stacked barplot ------------------------- - -ggplot(data_barplot, aes(x = region, y = n, alpha = unique)) + - geom_bar(position = "stack", stat = "identity", fill = "darkred") + - scale_alpha_manual( - name = "", - labels = c("Hub in multiple regions", "Hub unique"), - values = c(1.0, 0.5) - ) + - xlab("Brain region") + - ylab("# Hub genes") + - theme_light() + - theme( - axis.title.x = element_text(size = 15), - axis.title.y = element_text(size = 15), - axis.text.x = element_text(size = 12), - axis.text.y = element_text(size = 12), - legend.text = element_text(size = 12), - # legend.title = element_blank(), - legend.position = "top" - ) + - geom_text(aes(label = paste0(round((n / sum) * 100, digits = 1 - ), "%")), - position = position_stack(vjust = 0.5), - size = 4, - color = "white", - show.legend = FALSE) - -ggsave( - "11_FigureSII_B.pdf", - width = 8, - height = 6 -) diff --git a/07_PlotsManuscript/12_FigureSII_C.R b/07_PlotsManuscript/12_FigureSII_C.R deleted file mode 100644 index 9083131..0000000 --- a/07_PlotsManuscript/12_FigureSII_C.R +++ /dev/null @@ -1,248 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# GO enrichment for DE genes across all regions - -library(clusterProfiler) -library(org.Mm.eg.db) -library(ggplot2) -library(dplyr) -library(stringr) -library(writexl) - - -basepath <- "~/Documents/ownCloud/DexStim_RNAseq_Mouse/" - - -# 0. Read genes DE in all regions and background ----------------- -genes <- read.table(file = paste0(basepath, "tables/06_overlap_AMY-CER-PFC-PVN-dDG-vDG-dCA1-vCA1_entrezID.txt"), - header = FALSE)[,1] - -# background are all genes in out dataset -background <- read.table(file = paste0(basepath, "tables/06_background_entrezID.txt"), - header = FALSE)[,1] - - - -# 1.1 GO enrichment for genes DE in all regions --------------------- - -# GO enrichment -min_genes <- round(length(genes)/10) -ego <- enrichGO(gene = as.character(genes), - universe = as.character(background), - OrgDb = org.Mm.eg.db, - ont = "BP", - pAdjustMethod = "BH", - minGSSize = min_genes, # min number of genes associated with GO term - maxGSSize = 10000, # max number of genes associated with GO term - readable = TRUE)@result -head(ego, n = 20) - - -# 1.2 GO enrichment for upregulated genes DE in all regions --------------------- -genes_up <- read.table(file = paste0(basepath, "tables/06_overlap_AMY-CER-PFC-PVN-dDG-vDG-dCA1-vCA1_up_entrezID.txt"), - header = FALSE)[,1] - -# GO enrichment -min_genes <- round(length(genes_up)/10) -ego_up <- enrichGO(gene = as.character(genes_up), - universe = as.character(background), - OrgDb = org.Mm.eg.db, - ont = "BP", - pAdjustMethod = "BH", - minGSSize = min_genes, # min number of genes associated with GO term - maxGSSize = 10000, # max number of genes associated with GO term - readable = TRUE)@result -head(ego_up, n = 20) - - - - -# 1.3 GO enrichment for downregulated genes DE in all regions --------------------- -genes_down <- read.table( - file = paste0(basepath, - "tables/06_overlap_AMY-CER-PFC-PVN-dDG-vDG-dCA1-vCA1_down_entrezID.txt"), - header = FALSE)[,1] - -# GO enrichment -min_genes <- length(genes_down)/10 -ego_down <- enrichGO(gene = as.character(genes_down), - universe = as.character(background), - OrgDb = org.Mm.eg.db, - ont = "BP", - pAdjustMethod = "BH", - minGSSize = min_genes, # min number of genes associated with GO term - maxGSSize = 10000, # max number of genes associated with GO term - readable = TRUE)@result -head(ego_down, n = 20) - - -# 1.4 Merge dataframes from all, up and downregulated genes -------------------- - -data_go <- left_join(ego, ego_down, by = c("ID", "Description"), - suffix = c(".all", ".down")) -data_go <- left_join(data_go, ego_up, by = c("ID", "Description"), - suffix = c("", ".up")) -# apparently not all entrez ids are valid --> not all taken into account in enrichment analysis -ngenes_rel <- 165 -nback_rel <- 12089 -# add columns relevant for gene ratio and odds ratio -data_go <- data_go %>% - mutate(n_genes = as.numeric(str_extract(data_go$BgRatio.all, "^\\d+"))) %>% - mutate(gr = Count.all/n_genes*100) %>% - mutate(b = ngenes_rel - Count.all, - c = n_genes - Count.all, - d = nback_rel - Count.all, - or = log2(Count.all*d)-log2(b*c)) - -data_heat <- data_go[1:30,c("Description", "p.adjust.down", "p.adjust")] %>% - tidyr::pivot_longer(cols = p.adjust.down:p.adjust) - - -# 1.5 Barplot ------------------------------- -# gene ratio -bp.1 <- - ggplot(data = data_go[1:25,], aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = gr - )) + - geom_bar(stat = "identity", position = "stack", - fill = "#0F5057") + - ylab("% Gene Ratio") + - xlab("") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - -# odds ratio -bp.2 <- - ggplot(data = data_go[1:25,], aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = or - )) + - geom_bar(stat = "identity", position = "stack", - fill = "#FAA916") + - ylab("Odds Ratio") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - -# p-value -bp.3 <- - ggplot(data = data_go[1:25,], aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = -log10(p.adjust.all) - )) + - geom_bar(stat = "identity", position = "stack", - fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - ylab("-log10(FDR)") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - -hm.1 <- - ggplot(data = data_heat, aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = name, - fill = value <= 0.01 - )) + - geom_tile() + - scale_y_discrete(name ="sig. GO term", - limits=c("p.adjust","p.adjust.down"), - labels=c("upreg.", "downreg.")) + - scale_fill_manual( - name = "GO term sig.", - values = c("darkgrey", "black") - ) + - ylab("sig. GO term") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - #legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - -# combined plot -comb <- ggarrange(bp.1, - bp.2 + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - bp.3 + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - hm.1 + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(2.5,1,1,1)) - -# Save plot -ggexport( - comb, - filename = "12_FigureSII_C.pdf", - height = 10, - width = 15 -) - - -# bring table to a format similar as other tables -tf <- data_go[1:25,] - -# add Category column -tf$Category <- "GO_bp" -# add GeneSet column -tf$GeneSet <- tf$Description %>% - stringr::str_replace_all(" ", "_") %>% - stringr::str_to_upper() -tf$GeneSet <- paste0("GO_", tf$GeneSet) -# add -log10 transformed FDR pvalue -tf$`[-log10 FDR]` <- -log10(tf$p.adjust.all) - -# subset to columns that should be printed -tf <- tf %>% - dplyr::select(Category, GeneSet, n_genes, Count, b, c, d, gr, or, pvalue.all, - p.adjust.all, `[-log10 FDR]`, geneID.all) %>% - dplyr::rename(N_genes = n_genes, - N_overlap_A = Count, - `Input that does not overlap_B` = b, - `GO term genes - overlap_C` = c, - `Not GO term genes and not in my geneset_D` = d, - `Genes ratio` = gr, - `OR[log2(a*d)-log2(b*c)]` = or, - p = pvalue.all, - adjP = p.adjust.all, - genes = geneID.all) - -# save to a excel file -write_xlsx(tf, "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure II_PFC/GOterms_DEgenesAllRegions.xlsx") diff --git a/07_PlotsManuscript/13_FigureSII_D.R b/07_PlotsManuscript/13_FigureSII_D.R deleted file mode 100644 index 8bc2fca..0000000 --- a/07_PlotsManuscript/13_FigureSII_D.R +++ /dev/null @@ -1,199 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# GO enrichment for hub genes across all regions - - -library(data.table) -library(ggplot2) -library(dplyr) -library(tidyverse) -library(org.Mm.eg.db) -library(clusterProfiler) -library(ggpubr) -library(writexl) - -regions <- - c("AMY", "PFC", "PVN", "CER", "vDG", "dDG", "vCA1", "dCA1") - - -# 1. Read DE tables from all regions ---------- - -list_reg_sig <- list() -list_genes_sig <- list() - -for (reg in regions) { - res <- - fread( - file = paste0( - "~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/coExpression_kimono/03_AnalysisFuncoup/", - "/04_", - reg, - "_funcoup_differential_nodebetweennessNorm_betacutoff0.01.csv" - ), - sep = "," - ) - na_indices <- which(is.na(res$nodebetweenness_norm)) - res$padj[na_indices] <- 0 - res_sig <- res[res$nodebetweenness_norm >= 1.0, ] - list_reg_sig[[reg]] <- res_sig - list_genes_sig[[reg]] <- res_sig$ensembl_id -} - -hub_shared <- Reduce(intersect, list_genes_sig) -genes <- mapIds(org.Mm.eg.db, keys = hub_shared, keytype = "ENSEMBL", column="ENTREZID") - -# background are all genes in out dataset -background <- read.table(file = paste0("~/Documents/ownCloud/DexStim_RNAseq_Mouse/tables/06_background_entrezID.txt"), - header = FALSE)[,1] - - -# 1.1 GO enrichment for genes DE in all regions --------------------- - -# GO enrichment -min_genes <- round(length(genes)/10) -ego <- enrichGO(gene = as.character(genes), - universe = as.character(background), - OrgDb = org.Mm.eg.db, - ont = "BP", - pAdjustMethod = "BH", - minGSSize = min_genes, # min number of genes associated with GO term - maxGSSize = 10000, # max number of genes associated with GO term - readable = TRUE)@result -head(ego, n = 20) - - -# 1.4 Merge dataframes from all, up and downregulated genes -------------------- - -# apparently not all entrez ids are valid --> not all taken into account in enrichment analysis -ngenes_rel <- 7 -nback_rel <- 12089 -# add columns relevant for gene ratio and odds ratio -data_go <- ego %>% - mutate(n_genes = as.numeric(str_extract(ego$BgRatio, "^\\d+"))) %>% - mutate(gr = Count/n_genes*100) %>% - mutate(b = ngenes_rel - Count, - c = n_genes - Count, - d = nback_rel - Count, - or = log2(Count*d)-log2(b*c)) - - -# 1.5 Barplot ------------------------------- -# gene ratio -bp.1 <- - ggplot(data = data_go[1:25,], aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = gr - )) + - geom_bar(stat = "identity", position = "stack", - fill = "#0F5057") + - ylab("% Gene Ratio") + - xlab("") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - -# odds ratio -bp.2 <- - ggplot(data = data_go[1:25,], aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = or - )) + - geom_bar(stat = "identity", position = "stack", - fill = "#FAA916") + - ylab("Odds Ratio") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - -# p-value -bp.3 <- - ggplot(data = data_go[1:25,], aes( - x = factor(Description, levels = rev(data_go$Description[1:25])), - y = -log10(p.adjust) - )) + - geom_bar(stat = "identity", position = "stack", - fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - ylab("-log10(FDR)") + - theme_light() + - theme( - axis.title.y = element_text(size = 14), - axis.text.x = element_text(size = 12), - axis.title.x = element_text(size =14), - axis.text.y = element_text(size = 12), - legend.position = "top", - legend.text = element_text(size = 10) - ) + - coord_flip() - - -# combined plot -comb <- ggarrange(bp.1, - bp.2 + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - bp.3 + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(2.5,1,1)) - -# Save plot -ggexport( - comb, - filename = "13_FigureSII_D.pdf", - height = 10, - width = 13 -) - - -# bring table to a format similar as other tables -tf <- data_go[1:25,] - -# add Category column -tf$Category <- "GO_bp" -# add GeneSet column -tf$GeneSet <- tf$Description %>% - stringr::str_replace_all(" ", "_") %>% - stringr::str_to_upper() -tf$GeneSet <- paste0("GO_", tf$GeneSet) -# add -log10 transformed FDR pvalue -tf$`[-log10 FDR]` <- -log10(tf$p.adjust) - -# subset to columns that should be printed -tf <- tf %>% - dplyr::select(Category, GeneSet, n_genes, Count, b, c, d, gr, or, pvalue, - p.adjust, `[-log10 FDR]`, geneID) %>% - dplyr::rename(N_genes = n_genes, - N_overlap_A = Count, - `Input that does not overlap_B` = b, - `GO term genes - overlap_C` = c, - `Not GO term genes and not in my geneset_D` = d, - `Genes ratio` = gr, - `OR[log2(a*d)-log2(b*c)]` = or, - p = pvalue, - adjP = p.adjust, - genes = geneID) - -# save to a excel file -write_xlsx(tf, "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure II_PFC/GOterms_hubsAllRegions.xlsx") - diff --git a/07_PlotsManuscript/14_FigureSIV_gwas.R b/07_PlotsManuscript/14_FigureSIV_gwas.R deleted file mode 100644 index f7fa6a2..0000000 --- a/07_PlotsManuscript/14_FigureSIV_gwas.R +++ /dev/null @@ -1,75 +0,0 @@ -################################################## -## Project: DexStim Mouse Brain -## Date: 30.11.2021 -## Author: Nathalie -################################################## -# GWAS enrichment plot for manuscript -# Tcf4 neighbours? - -library(readxl) -library(dplyr) -library(ggplot2) -library(ggpubr) - -# 1. Pathway enrichment for Abcd1 and neighbours - -data <- read_xlsx(path = "~/Documents/ownCloud/DexStim_RNAseq_Mouse/manuscript/Figures/Figure IV_GWAS/TCF4/FUMA_diffTcf4PFCnetwork_gene2func68788/FUMA_diffPFCTcf4.xlsx", - sheet = "GWAS") - -data <- arrange(data, desc(`[-log10]`)) -data$`my adj p` <- as.numeric(data$`my adj p`) -data <- data[data$`my adj p` <= 0.05,] -data <- data[!is.na(data$`my adj p`),] - -# GeneSet as factor -data$GeneSet <- factor(data$GeneSet, levels = data$GeneSet) - - -# bar plot gene ratio -gr <- ggplot(data, aes(x = GeneSet, y = `Genes ratio`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#0F5057") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("% Gene Ratio") + - labs(x = NULL) - -# bar plot odds ratio -or <- ggplot(data, aes(x = GeneSet, y = `OR[log2(a*d)-log2(b*c)]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#FAA916") + - scale_fill_manual() + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("Odds Ratio") + - labs(x = NULL) - -# barplot fdr p-value -fdr <- ggplot(data, aes(x = GeneSet, y = `[-log10]`) ) + - geom_bar(position = position_dodge2(reverse=TRUE), stat="identity", fill = "#D45E60") + - geom_hline(yintercept = -log10(0.05),linetype="dashed", color = "red") + - theme_light() + - coord_flip() + - scale_x_discrete(limits = rev) + - theme(text = element_text(size= 14)) + - ylab("-log10(FDR)") + - labs(x = NULL) - -# combined barplot -comb <- ggarrange(gr, - or + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - fdr + - theme(axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() ), - nrow = 1, - widths = c(1.8,1,1)) - -ggexport(comb, filename = "14_FigureSIV_gwas.pdf", width = 14, height = 8) -