From e106c56a2150b86304de5ce8e8b4e6b3dffbd499 Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Wed, 30 Jun 2021 17:21:24 +0200 Subject: [PATCH 1/7] Started work on ttest and mann-whitney u test. --- DESCRIPTION | 1 + NAMESPACE | 4 ++ R/shinyAppServer.R | 93 +++++++++++++++++++++++++++++++++++++++++++--- R/shinyAppUI.R | 23 +++++++++++- 4 files changed, 114 insertions(+), 7 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2288d6d..b7a81b2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,6 +10,7 @@ License: GPL (>= 3) Encoding: UTF-8 LazyData: true Imports: + broom, dplyr, DT, ggCPM, diff --git a/NAMESPACE b/NAMESPACE index 5f3f721..d23f4e8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ import(shiny) import(tidyselect) importFrom(DT,DTOutput) importFrom(DT,renderDT) +importFrom(broom,tidy) importFrom(dplyr,across) importFrom(dplyr,arrange) importFrom(dplyr,bind_rows) @@ -91,7 +92,9 @@ importFrom(recipes,step_pca) importFrom(recipes,step_scale) importFrom(recipes,tidy) importFrom(recipes,update_role) +importFrom(rlang,"!!") importFrom(rlang,.data) +importFrom(rlang,sym) importFrom(scales,scientific) importFrom(sessioninfo,session_info) importFrom(shiny,NS) @@ -110,6 +113,7 @@ importFrom(stringr,str_split) importFrom(tibble,column_to_rownames) importFrom(tibble,tibble) importFrom(tidyr,everything) +importFrom(tidyr,nest) importFrom(tidyr,pivot_longer) importFrom(tidyr,pivot_wider) importFrom(tidyr,separate) diff --git a/R/shinyAppServer.R b/R/shinyAppServer.R index 8b0e239..7392e1e 100644 --- a/R/shinyAppServer.R +++ b/R/shinyAppServer.R @@ -11,10 +11,10 @@ #' @importFrom sessioninfo session_info #' @importFrom tibble tibble #' @importFrom dplyr filter mutate select pull distinct case_when relocate bind_rows across -#' @importFrom rlang .data +#' @importFrom rlang .data sym !! #' @importFrom purrr map #' @importFrom magrittr %>% -#' @importFrom tidyr unnest pivot_wider +#' @importFrom tidyr unnest pivot_wider nest #' @importFrom tidyselect last_col everything matches #' @importFrom utils head #' @importFrom tools file_ext @@ -23,6 +23,7 @@ #' @importFrom plotly renderPlotly plotlyOutput plot_ly add_markers event_data #' @importFrom shinycssloaders withSpinner #' @importFrom openxlsx write.xlsx +#' @importFrom broom tidy #' #' @author Rico Derks @@ -1331,6 +1332,7 @@ shinyAppServer <- function(input, output, session) { } }) + #### heatmap # update color selection observeEvent(input$select_group_column, { req(input$select_group_column) @@ -1344,7 +1346,6 @@ shinyAppServer <- function(input, output, session) { } }) - #### compare samples output$compare_samples <- renderPlotly({ req(all_data$lipid_data_filter, input$select_z_heatmap) @@ -1361,9 +1362,91 @@ shinyAppServer <- function(input, output, session) { z = input$select_z_heatmap, clust = input$heatmap_use_clust, sample_group = input$select_heatmap_group) - } + } }) - #### + #### end heatmap + + #### compare samples + # create some ui output + output$test_group_selection <- renderUI({ + req(all_data$analysis_data) + + if(all_data$merged_data == TRUE & !is.null(input$select_group_column)) { + tagList( + selectInput(inputId = "test_select_group", + label = "Select a group:", + choices = c("none", input$select_group_column), + selected = "none"), + uiOutput(outputId = "test_vs_groups") + ) + } + }) + + output$test_vs_groups <- renderUI({ + tagList( + selectInput(inputId = "test_group1", + label = "Group 1:", + choices = "none"), + HTML("
vs
"), + selectInput(inputId = "test_group2", + label = "Group 2:", + choices = "none") + ) + }) + + observeEvent(input$test_select_group, { + req(all_data$meta_data) + + if(input$test_select_group != "none") { + # get the groups + group_options <- all_data$meta_data() %>% + mutate(across(everything(), as.character)) %>% + select(matches(paste0("^", input$test_select_group, "$"))) %>% + pull() %>% + unique() + + # remove NA + group_options <- group_options[!is.na(group_options)] + + updateSelectInput(inputId = "test_group1", + label = "Group 1:", + choices = c("none", group_options)) + + updateSelectInput(inputId = "test_group2", + label = "Group 2:", + choices = c("none", group_options)) + } + }) + + observeEvent({ + input$test_group1 + input$test_group2 + }, { + # check if something is selected and not the same thing + if(input$test_group1 != "none" & + input$test_group2 != "none" & + input$test_group1 != input$test_group2) { + # get the column name + my_column <- input$test_select_group + + tmp <- all_data$analysis_data %>% + rename(my_group_info = !!sym(my_column)) %>% + filter(.data$my_group_info == input$test_group1 | + .data$my_group_info == input$test_group2) %>% + select(.data$my_id, .data$ShortLipidName, .data$LipidClass, .data$sample_name, .data$my_group_info, .data$area) %>% + nest(test_data = c(.data$sample_name, .data$my_group_info, .data$area)) %>% + mutate(model_ttest = map(.x = .data$test_data, + .f = ~ broom::tidy(t.test(area ~ my_group_info, + data = .x))), + model_mwtest = map(.x = .data$test_data, + .f = ~ broom::tidy(wilcox.test(area ~ my_group_info, + data = .x)))) + + print(tmp$model_ttest[[1]]) + } + }) + + #### end compare samples #### PCA # do the PCA analysis diff --git a/R/shinyAppUI.R b/R/shinyAppUI.R index 541470a..dced606 100644 --- a/R/shinyAppUI.R +++ b/R/shinyAppUI.R @@ -350,8 +350,8 @@ shinyAppUI <- fluidPage( ), # end tabpanel issues # start navbarMenu analysis navbarMenu(title = "Analysis", - # start tabpanel compare samples - tabPanel(title = "Compare samples", + # start tabpanel heatmap + tabPanel(title = "Heatmap", fluidPage( sidebarPanel(width = 2, radioButtons(inputId = "select_z_heatmap", @@ -375,6 +375,25 @@ shinyAppUI <- fluidPage( height = "900px"), type = 5)) ) + ), # end tabpanel heatmap + # start compare samples + tabPanel(title = "Compare samples", + fluidPage( + sidebarPanel(width = 3, + radioButtons(inputId = "select_test", + label = "Test:", + choices = c("t-test" = "ttest", + "Mann-Whitney U test" = "mwtest"), + selected = "ttest"), + uiOutput(outputId = "test_group_selection") + ), + mainPanel(width = 9, + p("After merge do some group comparissons") + # shinycssloaders::withSpinner(plotlyOutput(outputId = "volcano_plot", + # height = "900px"), + # type = 5) + ) + ) ), # end tabpanel compare samples # start tabPanel pca analysis tabPanel(title = "PCA", From 2040265f35f0e664f53ae0eafe9b408fa475dc58 Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Thu, 1 Jul 2021 09:30:31 +0200 Subject: [PATCH 2/7] T-test volcano plot is there. --- NAMESPACE | 2 ++ R/do_ttest.R | 46 ++++++++++++++++++++++++++++++++++++++ R/shinyAppServer.R | 35 +++++++++++++++++++---------- R/shinyAppUI.R | 7 +++--- R/volcano_plot.R | 54 +++++++++++++++++++++++++++++++++++++++++++++ man/do_ttest.Rd | 28 +++++++++++++++++++++++ man/volcano_plot.Rd | 20 +++++++++++++++++ 7 files changed, 176 insertions(+), 16 deletions(-) create mode 100644 R/do_ttest.R create mode 100644 R/volcano_plot.R create mode 100644 man/do_ttest.Rd create mode 100644 man/volcano_plot.Rd diff --git a/NAMESPACE b/NAMESPACE index d23f4e8..34b1c19 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -75,6 +75,7 @@ importFrom(plotly,plot_ly) importFrom(plotly,plotlyOutput) importFrom(plotly,renderPlotly) importFrom(purrr,map) +importFrom(purrr,map_dbl) importFrom(readr,col_character) importFrom(readr,col_double) importFrom(readr,col_integer) @@ -106,6 +107,7 @@ importFrom(shinyjs,hidden) importFrom(shinyjs,toggle) importFrom(shinyjs,useShinyjs) importFrom(stats,cor) +importFrom(stats,p.adjust) importFrom(stats,sd) importFrom(stringr,str_extract) importFrom(stringr,str_replace) diff --git a/R/do_ttest.R b/R/do_ttest.R new file mode 100644 index 0000000..86af020 --- /dev/null +++ b/R/do_ttest.R @@ -0,0 +1,46 @@ +#' @title Do a t-test on all lipids +#' +#' @description Do a t-test on all lipids. +#' +#' @param lipid_data tibble in tidy format, already nested +#' +#' @details A t-test will be done for each lipid. Also the p-value will be corrected +#' for multiple testing. +#' +#' @return a tibble with the results and the following columns +#' model_test contains the model information for each lipid (nested) +#' pvalue contains the uncorrected pvalue +#' pvalue_cor contains the corrected pvalue +#' fc the fold change estimate1/estimate2 +#' +#' @import tidyselect +#' @importFrom dplyr select filter mutate left_join case_when distinct across +#' @importFrom tidyr pivot_wider everything +#' @importFrom magrittr %>% +#' @importFrom rlang .data +#' @importFrom purrr map map_dbl +#' @importFrom broom tidy +#' @importFrom stats p.adjust +#' +#' @author Rico Derks +#' +do_ttest <- function(lipid_data) { + results <- lipid_data %>% + mutate(model_test = map(.x = .data$test_data, + .f = ~ broom::tidy(t.test(formula = area ~ my_group_info, + data = .x)) %>% + # calculate the fold change + mutate(fc = estimate1 / estimate2)), + pvalue = map_dbl(.x = .data$model_test, + .f = ~ .x$p.value), + pvalue_adj = p.adjust(.data$pvalue, + method = "BH"), + fc = map_dbl(.x = .data$model_test, + .f = ~ .x$fc), + p_log10 = -log10(.data$pvalue), + p_log10_adj = -log10(.data$pvalue_adj), + fc_log2 = log2(.data$fc)) + + return(results) +} + diff --git a/R/shinyAppServer.R b/R/shinyAppServer.R index 7392e1e..b1f54fc 100644 --- a/R/shinyAppServer.R +++ b/R/shinyAppServer.R @@ -23,7 +23,6 @@ #' @importFrom plotly renderPlotly plotlyOutput plot_ly add_markers event_data #' @importFrom shinycssloaders withSpinner #' @importFrom openxlsx write.xlsx -#' @importFrom broom tidy #' #' @author Rico Derks @@ -1418,7 +1417,7 @@ shinyAppServer <- function(input, output, session) { } }) - observeEvent({ + test_result <- eventReactive({ input$test_group1 input$test_group2 }, { @@ -1429,20 +1428,32 @@ shinyAppServer <- function(input, output, session) { # get the column name my_column <- input$test_select_group - tmp <- all_data$analysis_data %>% + # prepare the data for the testing + prep_test_data <- all_data$analysis_data %>% rename(my_group_info = !!sym(my_column)) %>% filter(.data$my_group_info == input$test_group1 | .data$my_group_info == input$test_group2) %>% select(.data$my_id, .data$ShortLipidName, .data$LipidClass, .data$sample_name, .data$my_group_info, .data$area) %>% - nest(test_data = c(.data$sample_name, .data$my_group_info, .data$area)) %>% - mutate(model_ttest = map(.x = .data$test_data, - .f = ~ broom::tidy(t.test(area ~ my_group_info, - data = .x))), - model_mwtest = map(.x = .data$test_data, - .f = ~ broom::tidy(wilcox.test(area ~ my_group_info, - data = .x)))) - - print(tmp$model_ttest[[1]]) + nest(test_data = c(.data$sample_name, .data$my_group_info, .data$area)) + + if(input$select_test == "ttest") { + # do the t-test + results_test <- do_ttest(lipid_data = prep_test_data) + } + # model_mwtest = map(.x = .data$test_data, + # .f = ~ broom::tidy(wilcox.test(area ~ my_group_info, + # data = .x)))) + return(results_test) + } else { + return(NULL) + } + }) + + output$volcano_plot <- renderPlotly({ + req(test_result) + + if(!is.null(test_result())) { + volcano_plot(lipid_data = test_result()) } }) diff --git a/R/shinyAppUI.R b/R/shinyAppUI.R index dced606..5ee9b3a 100644 --- a/R/shinyAppUI.R +++ b/R/shinyAppUI.R @@ -388,10 +388,9 @@ shinyAppUI <- fluidPage( uiOutput(outputId = "test_group_selection") ), mainPanel(width = 9, - p("After merge do some group comparissons") - # shinycssloaders::withSpinner(plotlyOutput(outputId = "volcano_plot", - # height = "900px"), - # type = 5) + shinycssloaders::withSpinner(plotlyOutput(outputId = "volcano_plot", + height = "500"), + type = 5) ) ) ), # end tabpanel compare samples diff --git a/R/volcano_plot.R b/R/volcano_plot.R new file mode 100644 index 0000000..1fbdc8d --- /dev/null +++ b/R/volcano_plot.R @@ -0,0 +1,54 @@ +#' @title Create volcano plot +#' +#' @description Create volcano plot. +#' +#' @param lipid_data tibble with all the lipid data and test data +#' +#' @return plotly object +#' +#' @importFrom magrittr %>% +#' @importFrom plotly plot_ly add_markers layout +#' +#' @author Rico Derks +#' +volcano_plot <- function(lipid_data) { + p <- lipid_data %>% + plot_ly(x = ~fc_log2, + y = ~p_log10) %>% + add_markers() %>% + layout(xaxis = list(zeroline = FALSE), + #yaxis = list(zeroline = FALSE), + shapes = list(vline(-1), + vline(1), + hline(-log10(0.05)))) + + return(p) +} + +vline <- function(x = 0, color = "blue") { + list( + type = "line", + y0 = 0, + y1 = 1, + yref = "paper", + x0 = x, + x1 = x, + line = list(color = color, + width = 1, + dash = "dash") + ) +} + +hline <- function(y = 0, color = "blue") { + list( + type = "line", + x0 = 0, + x1 = 1, + xref = "paper", + y0 = y, + y1 = y, + line = list(color = color, + width = 1, + dash = "dash") + ) +} diff --git a/man/do_ttest.Rd b/man/do_ttest.Rd new file mode 100644 index 0000000..214667a --- /dev/null +++ b/man/do_ttest.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/do_ttest.R +\name{do_ttest} +\alias{do_ttest} +\title{Do a t-test on all lipids} +\usage{ +do_ttest(lipid_data) +} +\arguments{ +\item{lipid_data}{tibble in tidy format, already nested} +} +\value{ +a tibble with the results and the following columns + model_test contains the model information for each lipid (nested) + pvalue contains the uncorrected pvalue + pvalue_cor contains the corrected pvalue + fc the fold change estimate1/estimate2 +} +\description{ +Do a t-test on all lipids. +} +\details{ +A t-test will be done for each lipid. Also the p-value will be corrected + for multiple testing. +} +\author{ +Rico Derks +} diff --git a/man/volcano_plot.Rd b/man/volcano_plot.Rd new file mode 100644 index 0000000..b81f217 --- /dev/null +++ b/man/volcano_plot.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/volcano_plot.R +\name{volcano_plot} +\alias{volcano_plot} +\title{Create volcano plot} +\usage{ +volcano_plot(lipid_data) +} +\arguments{ +\item{lipid_data}{tibble with all the lipid data and test data} +} +\value{ +plotly object +} +\description{ +Create volcano plot. +} +\author{ +Rico Derks +} From a1fdcc78e8bc3359c5e87e47474e6c80ac4c5cda Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Thu, 1 Jul 2021 11:06:06 +0200 Subject: [PATCH 3/7] Added mann-whitney u test. --- R/do_mwtest.R | 41 +++++++++++++++++++++++++++++++++++++++++ R/do_ttest.R | 8 ++------ R/shinyAppServer.R | 27 +++++++++++++++++---------- R/shinyAppUI.R | 3 +++ R/volcano_plot.R | 24 +++++++++++++++++++----- man/do_mwtest.Rd | 28 ++++++++++++++++++++++++++++ man/volcano_plot.Rd | 4 +++- 7 files changed, 113 insertions(+), 22 deletions(-) create mode 100644 R/do_mwtest.R create mode 100644 man/do_mwtest.Rd diff --git a/R/do_mwtest.R b/R/do_mwtest.R new file mode 100644 index 0000000..5825f52 --- /dev/null +++ b/R/do_mwtest.R @@ -0,0 +1,41 @@ +#' @title Do a Mann-Whitney U test on all lipids +#' +#' @description Do a Mann-Whitney U test on all lipids +#' +#' @param lipid_data tibble in tidy format, already nested +#' +#' @details A Mann-Whitney U test will be done for each lipid. Also the p-value will be corrected +#' for multiple testing. +#' +#' @return a tibble with the results and the following columns +#' model_test contains the model information for each lipid (nested) +#' pvalue contains the uncorrected pvalue +#' pvalue_cor contains the corrected pvalue +#' fc the fold change estimate1/estimate2 +#' +#' @import tidyselect +#' @importFrom dplyr mutate group_by summarise pull +#' @importFrom tidyr pivot_wider +#' @importFrom magrittr %>% +#' @importFrom rlang .data +#' @importFrom purrr map map_dbl +#' @importFrom broom tidy +#' @importFrom stats p.adjust +#' +#' @author Rico Derks +#' +do_mwtest <- function(lipid_data) { + results <- lipid_data %>% + mutate(model_test = map(.x = .data$test_data, + .f = ~ broom::tidy(wilcox.test(formula = area ~ my_group_info, + data = .x))), + pvalue = map_dbl(.x = .data$model_test, + .f = ~ .x$p.value), + pvalue_adj = p.adjust(.data$pvalue, + method = "BH"), + p_log10 = -log10(.data$pvalue), + p_log10_adj = -log10(.data$pvalue_adj)) + + return(results) +} + diff --git a/R/do_ttest.R b/R/do_ttest.R index 86af020..a9aeb44 100644 --- a/R/do_ttest.R +++ b/R/do_ttest.R @@ -14,8 +14,7 @@ #' fc the fold change estimate1/estimate2 #' #' @import tidyselect -#' @importFrom dplyr select filter mutate left_join case_when distinct across -#' @importFrom tidyr pivot_wider everything +#' @importFrom dplyr mutate #' @importFrom magrittr %>% #' @importFrom rlang .data #' @importFrom purrr map map_dbl @@ -35,11 +34,8 @@ do_ttest <- function(lipid_data) { .f = ~ .x$p.value), pvalue_adj = p.adjust(.data$pvalue, method = "BH"), - fc = map_dbl(.x = .data$model_test, - .f = ~ .x$fc), p_log10 = -log10(.data$pvalue), - p_log10_adj = -log10(.data$pvalue_adj), - fc_log2 = log2(.data$fc)) + p_log10_adj = -log10(.data$pvalue_adj)) return(results) } diff --git a/R/shinyAppServer.R b/R/shinyAppServer.R index b1f54fc..43eaa74 100644 --- a/R/shinyAppServer.R +++ b/R/shinyAppServer.R @@ -12,7 +12,7 @@ #' @importFrom tibble tibble #' @importFrom dplyr filter mutate select pull distinct case_when relocate bind_rows across #' @importFrom rlang .data sym !! -#' @importFrom purrr map +#' @importFrom purrr map map_dbl #' @importFrom magrittr %>% #' @importFrom tidyr unnest pivot_wider nest #' @importFrom tidyselect last_col everything matches @@ -1420,7 +1420,12 @@ shinyAppServer <- function(input, output, session) { test_result <- eventReactive({ input$test_group1 input$test_group2 + input$select_test }, { + req(input$test_group1, + input$test_group2, + input$select_test) + # check if something is selected and not the same thing if(input$test_group1 != "none" & input$test_group2 != "none" & @@ -1434,15 +1439,16 @@ shinyAppServer <- function(input, output, session) { filter(.data$my_group_info == input$test_group1 | .data$my_group_info == input$test_group2) %>% select(.data$my_id, .data$ShortLipidName, .data$LipidClass, .data$sample_name, .data$my_group_info, .data$area) %>% - nest(test_data = c(.data$sample_name, .data$my_group_info, .data$area)) + nest(test_data = c(.data$sample_name, .data$my_group_info, .data$area)) %>% + mutate(fc = map_dbl(.x = .data$test_data, + .f = ~ mean(.x$area[.x$my_group_info == input$test_group1]) / mean(.x$area[.x$my_group_info == input$test_group2])), + fc_log2 = log2(.data$fc)) + + # what test to do + results_test <- switch(input$select_test, + "ttest" = do_ttest(lipid_data = prep_test_data), + "mwtest" = do_mwtest(lipid_data = prep_test_data)) - if(input$select_test == "ttest") { - # do the t-test - results_test <- do_ttest(lipid_data = prep_test_data) - } - # model_mwtest = map(.x = .data$test_data, - # .f = ~ broom::tidy(wilcox.test(area ~ my_group_info, - # data = .x)))) return(results_test) } else { return(NULL) @@ -1453,7 +1459,8 @@ shinyAppServer <- function(input, output, session) { req(test_result) if(!is.null(test_result())) { - volcano_plot(lipid_data = test_result()) + volcano_plot(lipid_data = test_result(), + pvalue_adjust = input$test_cor_pvalue) } }) diff --git a/R/shinyAppUI.R b/R/shinyAppUI.R index 5ee9b3a..c95feba 100644 --- a/R/shinyAppUI.R +++ b/R/shinyAppUI.R @@ -385,6 +385,9 @@ shinyAppUI <- fluidPage( choices = c("t-test" = "ttest", "Mann-Whitney U test" = "mwtest"), selected = "ttest"), + checkboxInput(inputId = "test_cor_pvalue", + label = "Show corrected p-value", + value = FALSE), uiOutput(outputId = "test_group_selection") ), mainPanel(width = 9, diff --git a/R/volcano_plot.R b/R/volcano_plot.R index 1fbdc8d..b3cb8d6 100644 --- a/R/volcano_plot.R +++ b/R/volcano_plot.R @@ -3,21 +3,35 @@ #' @description Create volcano plot. #' #' @param lipid_data tibble with all the lipid data and test data +#' @param pvalue_adjust show the corrected p value, default is FALSE #' #' @return plotly object #' #' @importFrom magrittr %>% +#' @importFrom dplyr mutate case_when +#' @importFrom rlang .data #' @importFrom plotly plot_ly add_markers layout #' #' @author Rico Derks #' -volcano_plot <- function(lipid_data) { +volcano_plot <- function(lipid_data, pvalue_adjust = FALSE) { + # create y-axis title + y_title <- ifelse(pvalue_adjust == FALSE, + "-log10(p value)", + "-log10(cor. p value)") + + # create the plot p <- lipid_data %>% + mutate(show_p = case_when( + pvalue_adjust == FALSE ~ .data$p_log10, + pvalue_adjust == TRUE ~ .data$p_log10_adj + )) %>% plot_ly(x = ~fc_log2, - y = ~p_log10) %>% - add_markers() %>% - layout(xaxis = list(zeroline = FALSE), - #yaxis = list(zeroline = FALSE), + y = ~show_p) %>% + add_markers(color = ~LipidClass) %>% + layout(xaxis = list(zeroline = FALSE, + title = "log2(fold change"), + yaxis = list(title = y_title), shapes = list(vline(-1), vline(1), hline(-log10(0.05)))) diff --git a/man/do_mwtest.Rd b/man/do_mwtest.Rd new file mode 100644 index 0000000..b4dbe99 --- /dev/null +++ b/man/do_mwtest.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/do_mwtest.R +\name{do_mwtest} +\alias{do_mwtest} +\title{Do a Mann-Whitney U test on all lipids} +\usage{ +do_mwtest(lipid_data) +} +\arguments{ +\item{lipid_data}{tibble in tidy format, already nested} +} +\value{ +a tibble with the results and the following columns + model_test contains the model information for each lipid (nested) + pvalue contains the uncorrected pvalue + pvalue_cor contains the corrected pvalue + fc the fold change estimate1/estimate2 +} +\description{ +Do a Mann-Whitney U test on all lipids +} +\details{ +A Mann-Whitney U test will be done for each lipid. Also the p-value will be corrected + for multiple testing. +} +\author{ +Rico Derks +} diff --git a/man/volcano_plot.Rd b/man/volcano_plot.Rd index b81f217..189f798 100644 --- a/man/volcano_plot.Rd +++ b/man/volcano_plot.Rd @@ -4,10 +4,12 @@ \alias{volcano_plot} \title{Create volcano plot} \usage{ -volcano_plot(lipid_data) +volcano_plot(lipid_data, pvalue_adjust = FALSE) } \arguments{ \item{lipid_data}{tibble with all the lipid data and test data} + +\item{pvalue_adjust}{show the corrected p value, default is FALSE} } \value{ plotly object From d262ebfb2609375701a93ebdec0a53b8ceb0d0ba Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Thu, 1 Jul 2021 11:49:55 +0200 Subject: [PATCH 4/7] Some work improving several plots. --- R/pca_observation_plot.R | 3 ++- R/pca_scores_plot.R | 3 ++- R/pca_variable_plot.R | 3 ++- R/volcano_plot.R | 16 ++++++++++++---- 4 files changed, 18 insertions(+), 7 deletions(-) diff --git a/R/pca_observation_plot.R b/R/pca_observation_plot.R index d245c73..70fe9f2 100644 --- a/R/pca_observation_plot.R +++ b/R/pca_observation_plot.R @@ -29,7 +29,8 @@ pca_observation_plot <- function(obs_data, var_name, color_by = "none") { color = ~color_group) %>% add_bars() %>% layout(title = list(text = paste("Observation plot for", var_name), - x = 0)) %>% + x = 0), + xaxis = list(title = "Sample name")) %>% hide_legend() %>% config(displayModeBar = FALSE) diff --git a/R/pca_scores_plot.R b/R/pca_scores_plot.R index 702609c..2a57a0c 100644 --- a/R/pca_scores_plot.R +++ b/R/pca_scores_plot.R @@ -37,7 +37,8 @@ pca_scores_plot <- function(scores_data, xaxis = "PC1", yaxis = "PC2", color_by layout(title = list(text = "Scores plot", x = 0), xaxis = list(title = xaxis), - yaxis = list(title = yaxis)) %>% + yaxis = list(title = yaxis), + legend = list(orientation = "h")) %>% # config(displayModeBar = FALSE) %>% event_register(event = "plotly_click") diff --git a/R/pca_variable_plot.R b/R/pca_variable_plot.R index 5d5eab3..f599f23 100644 --- a/R/pca_variable_plot.R +++ b/R/pca_variable_plot.R @@ -32,7 +32,8 @@ pca_variable_plot <- function(var_data, sample_name) { categoryorder = "array", categoryarray = ~order_x)) %>% layout(title = list(text = paste("Variable plot", sample_name), - x = 0)) %>% + x = 0), + xaxis = list(title = "Lipid")) %>% hide_legend() return(p) diff --git a/R/volcano_plot.R b/R/volcano_plot.R index b3cb8d6..d538e72 100644 --- a/R/volcano_plot.R +++ b/R/volcano_plot.R @@ -10,7 +10,8 @@ #' @importFrom magrittr %>% #' @importFrom dplyr mutate case_when #' @importFrom rlang .data -#' @importFrom plotly plot_ly add_markers layout +#' @importFrom plotly plot_ly add_markers layout event_register +#' @importFrom grDevices rainbow #' #' @author Rico Derks #' @@ -27,14 +28,21 @@ volcano_plot <- function(lipid_data, pvalue_adjust = FALSE) { pvalue_adjust == TRUE ~ .data$p_log10_adj )) %>% plot_ly(x = ~fc_log2, - y = ~show_p) %>% - add_markers(color = ~LipidClass) %>% + y = ~show_p, + colors = rainbow(n = 100), + customdata = lipid_data$ShortLipidName, + source = "volcano_plot") %>% + add_markers(color = ~LipidClass, + size = 3) %>% layout(xaxis = list(zeroline = FALSE, title = "log2(fold change"), yaxis = list(title = y_title), shapes = list(vline(-1), vline(1), - hline(-log10(0.05)))) + hline(-log10(0.05))), + legend = list(orientation = "h", + size = 1)) %>% + event_register(event = "plotly_click") return(p) } From 2d61a36912a49a3af476f3cd13c1294cd22c58f6 Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Thu, 1 Jul 2021 13:45:17 +0200 Subject: [PATCH 5/7] Added boxplot --- R/box_plot.R | 36 ++++++++++++++++++++++++++++++++++++ R/shinyAppServer.R | 25 +++++++++++++++++++++++++ R/shinyAppUI.R | 18 +++++++++++------- R/volcano_plot.R | 6 +++--- man/box_plot.Rd | 22 ++++++++++++++++++++++ 5 files changed, 97 insertions(+), 10 deletions(-) create mode 100644 R/box_plot.R create mode 100644 man/box_plot.Rd diff --git a/R/box_plot.R b/R/box_plot.R new file mode 100644 index 0000000..df62d76 --- /dev/null +++ b/R/box_plot.R @@ -0,0 +1,36 @@ +#' @title Create box plot +#' +#' @description Create box plot. +#' +#' @param lipid_data tibble with all the lipid data and test data +#' @param title title of the plot +#' +#' @return plotly object +#' +#' @importFrom magrittr %>% +#' @importFrom dplyr filter select +#' @importFrom tidyr unnest +#' @importFrom rlang .data +#' @importFrom plotly plot_ly layout hide_legend +#' +#' @author Rico Derks +#' +box_plot <- function(lipid_data, title = "") { + # create the plot + p <- lipid_data %>% + plot_ly(x = ~my_group_info, + y = ~area, + text = ~sample_name, + color = ~my_group_info, + type = "box", + boxpoints = "all", + jitter = 0.4, + pointpos = 0) %>% + layout(title = list(text = title, + x = 0), + yaxis = list(title = "Value"), + xaxis = list(title = "Group")) %>% + hide_legend() + + return(p) +} diff --git a/R/shinyAppServer.R b/R/shinyAppServer.R index 43eaa74..6ab0378 100644 --- a/R/shinyAppServer.R +++ b/R/shinyAppServer.R @@ -1464,6 +1464,31 @@ shinyAppServer <- function(input, output, session) { } }) + output$test_boxplot <- renderPlotly({ + req(test_result) + + if(!is.null(test_result())) { + # capture the click event + # this contains a column with the shortlipidname (column name: customdata) + my_data <- event_data(event = "plotly_click", + source = "volcano_plot_click") + + if(!is.null(my_data)) { + # restructure data + plot_data <- test_result() %>% + filter(.data$ShortLipidName %in% my_data$customdata) %>% + select(.data$test_data) %>% + unnest(.data$test_data) + + # show the boxplot + box_plot(lipid_data = plot_data, + title = paste0("Lipid: ", my_data$customdata)) + } + } + + + }) + #### end compare samples #### PCA diff --git a/R/shinyAppUI.R b/R/shinyAppUI.R index c95feba..9004c76 100644 --- a/R/shinyAppUI.R +++ b/R/shinyAppUI.R @@ -379,21 +379,25 @@ shinyAppUI <- fluidPage( # start compare samples tabPanel(title = "Compare samples", fluidPage( - sidebarPanel(width = 3, + sidebarPanel(width = 2, radioButtons(inputId = "select_test", label = "Test:", choices = c("t-test" = "ttest", "Mann-Whitney U test" = "mwtest"), selected = "ttest"), checkboxInput(inputId = "test_cor_pvalue", - label = "Show corrected p-value", - value = FALSE), + label = "Show corrected p-value", + value = FALSE), uiOutput(outputId = "test_group_selection") ), - mainPanel(width = 9, - shinycssloaders::withSpinner(plotlyOutput(outputId = "volcano_plot", - height = "500"), - type = 5) + mainPanel(width = 10, + column(width = 6, + shinycssloaders::withSpinner(plotlyOutput(outputId = "volcano_plot", + height = "900px"), + type = 5)), + column(width = 6, + shinycssloaders::withSpinner(plotlyOutput(outputId = "test_boxplot"), + type = 5)) ) ) ), # end tabpanel compare samples diff --git a/R/volcano_plot.R b/R/volcano_plot.R index d538e72..0a6f806 100644 --- a/R/volcano_plot.R +++ b/R/volcano_plot.R @@ -29,9 +29,10 @@ volcano_plot <- function(lipid_data, pvalue_adjust = FALSE) { )) %>% plot_ly(x = ~fc_log2, y = ~show_p, + text = ~ShortLipidName, colors = rainbow(n = 100), customdata = lipid_data$ShortLipidName, - source = "volcano_plot") %>% + source = "volcano_plot_click") %>% add_markers(color = ~LipidClass, size = 3) %>% layout(xaxis = list(zeroline = FALSE, @@ -40,8 +41,7 @@ volcano_plot <- function(lipid_data, pvalue_adjust = FALSE) { shapes = list(vline(-1), vline(1), hline(-log10(0.05))), - legend = list(orientation = "h", - size = 1)) %>% + legend = list(orientation = "h")) %>% event_register(event = "plotly_click") return(p) diff --git a/man/box_plot.Rd b/man/box_plot.Rd new file mode 100644 index 0000000..08fc9b9 --- /dev/null +++ b/man/box_plot.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/box_plot.R +\name{box_plot} +\alias{box_plot} +\title{Create box plot} +\usage{ +box_plot(lipid_data, title = "") +} +\arguments{ +\item{lipid_data}{tibble with all the lipid data and test data} + +\item{title}{title of the plot} +} +\value{ +plotly object +} +\description{ +Create box plot. +} +\author{ +Rico Derks +} From a18f1084274f9e1c13031dd2ebe5fd408bfaf260 Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Thu, 1 Jul 2021 14:41:29 +0200 Subject: [PATCH 6/7] Added normalization and transformation for statistical tests. --- R/box_plot.R | 5 +++-- R/do_mwtest.R | 2 +- R/do_ttest.R | 2 +- R/shinyAppServer.R | 33 +++++++++++++++++++++++++++++---- R/shinyAppUI.R | 10 ++++++++++ R/volcano_plot.R | 9 ++++++--- man/volcano_plot.Rd | 4 +++- 7 files changed, 53 insertions(+), 12 deletions(-) diff --git a/R/box_plot.R b/R/box_plot.R index df62d76..12a136a 100644 --- a/R/box_plot.R +++ b/R/box_plot.R @@ -19,7 +19,7 @@ box_plot <- function(lipid_data, title = "") { # create the plot p <- lipid_data %>% plot_ly(x = ~my_group_info, - y = ~area, + y = ~value, text = ~sample_name, color = ~my_group_info, type = "box", @@ -28,7 +28,8 @@ box_plot <- function(lipid_data, title = "") { pointpos = 0) %>% layout(title = list(text = title, x = 0), - yaxis = list(title = "Value"), + yaxis = list(title = "Value", + exponentformat = "E"), xaxis = list(title = "Group")) %>% hide_legend() diff --git a/R/do_mwtest.R b/R/do_mwtest.R index 5825f52..e829054 100644 --- a/R/do_mwtest.R +++ b/R/do_mwtest.R @@ -27,7 +27,7 @@ do_mwtest <- function(lipid_data) { results <- lipid_data %>% mutate(model_test = map(.x = .data$test_data, - .f = ~ broom::tidy(wilcox.test(formula = area ~ my_group_info, + .f = ~ broom::tidy(wilcox.test(formula = value ~ my_group_info, data = .x))), pvalue = map_dbl(.x = .data$model_test, .f = ~ .x$p.value), diff --git a/R/do_ttest.R b/R/do_ttest.R index a9aeb44..b670e89 100644 --- a/R/do_ttest.R +++ b/R/do_ttest.R @@ -26,7 +26,7 @@ do_ttest <- function(lipid_data) { results <- lipid_data %>% mutate(model_test = map(.x = .data$test_data, - .f = ~ broom::tidy(t.test(formula = area ~ my_group_info, + .f = ~ broom::tidy(t.test(formula = value ~ my_group_info, data = .x)) %>% # calculate the fold change mutate(fc = estimate1 / estimate2)), diff --git a/R/shinyAppServer.R b/R/shinyAppServer.R index 6ab0378..2bafe16 100644 --- a/R/shinyAppServer.R +++ b/R/shinyAppServer.R @@ -1421,10 +1421,14 @@ shinyAppServer <- function(input, output, session) { input$test_group1 input$test_group2 input$select_test + input$select_test_normalization + input$select_test_transformation }, { req(input$test_group1, input$test_group2, - input$select_test) + input$select_test, + input$select_test_normalization, + input$select_test_transformation) # check if something is selected and not the same thing if(input$test_group1 != "none" & @@ -1432,6 +1436,10 @@ shinyAppServer <- function(input, output, session) { input$test_group1 != input$test_group2) { # get the column name my_column <- input$test_select_group + # get the normalization + normalization <- input$select_test_normalization + # get the transformation + transformation <- input$select_test_transformation # prepare the data for the testing prep_test_data <- all_data$analysis_data %>% @@ -1439,9 +1447,25 @@ shinyAppServer <- function(input, output, session) { filter(.data$my_group_info == input$test_group1 | .data$my_group_info == input$test_group2) %>% select(.data$my_id, .data$ShortLipidName, .data$LipidClass, .data$sample_name, .data$my_group_info, .data$area) %>% - nest(test_data = c(.data$sample_name, .data$my_group_info, .data$area)) %>% + # total area normalisation + group_by(.data$sample_name) %>% + mutate(norm_area = .data$area / sum(.data$area)) %>% + ungroup() %>% + # select which normalization to use for PCA + mutate(value = case_when( + normalization == "raw" ~ .data$area, + normalization == "tot_area" ~ .data$norm_area + )) %>% + # do transformations and select which transformation to keep + mutate(value = case_when( + transformation == "none" ~ .data$value, + transformation == "log10" ~ log10(.data$value + 1) # the +1 is correct for any zero's + )) %>% + # remove the 2 area columns + select(-.data$area, -.data$norm_area) %>% + nest(test_data = c(.data$sample_name, .data$my_group_info, .data$value)) %>% mutate(fc = map_dbl(.x = .data$test_data, - .f = ~ mean(.x$area[.x$my_group_info == input$test_group1]) / mean(.x$area[.x$my_group_info == input$test_group2])), + .f = ~ mean(.x$value[.x$my_group_info == input$test_group1]) / mean(.x$value[.x$my_group_info == input$test_group2])), fc_log2 = log2(.data$fc)) # what test to do @@ -1460,7 +1484,8 @@ shinyAppServer <- function(input, output, session) { if(!is.null(test_result())) { volcano_plot(lipid_data = test_result(), - pvalue_adjust = input$test_cor_pvalue) + pvalue_adjust = input$test_cor_pvalue, + title = paste0(input$test_group1, " vs ", input$test_group2)) } }) diff --git a/R/shinyAppUI.R b/R/shinyAppUI.R index 9004c76..a495063 100644 --- a/R/shinyAppUI.R +++ b/R/shinyAppUI.R @@ -385,6 +385,16 @@ shinyAppUI <- fluidPage( choices = c("t-test" = "ttest", "Mann-Whitney U test" = "mwtest"), selected = "ttest"), + radioButtons(inputId = "select_test_normalization", + label = "Normalization:", + choices = c("Raw data" = "raw", + "Total area normalization" = "tot_area"), + selected = "tot_area"), + radioButtons(inputId = "select_test_transformation", + label = "Transformation:", + choices = c("No transformation" = "none", + "Log10 transformation" = "log10"), + selected = "none"), checkboxInput(inputId = "test_cor_pvalue", label = "Show corrected p-value", value = FALSE), diff --git a/R/volcano_plot.R b/R/volcano_plot.R index 0a6f806..4793640 100644 --- a/R/volcano_plot.R +++ b/R/volcano_plot.R @@ -4,6 +4,7 @@ #' #' @param lipid_data tibble with all the lipid data and test data #' @param pvalue_adjust show the corrected p value, default is FALSE +#' @param title title of the plot #' #' @return plotly object #' @@ -15,7 +16,7 @@ #' #' @author Rico Derks #' -volcano_plot <- function(lipid_data, pvalue_adjust = FALSE) { +volcano_plot <- function(lipid_data, pvalue_adjust = FALSE, title = "") { # create y-axis title y_title <- ifelse(pvalue_adjust == FALSE, "-log10(p value)", @@ -36,12 +37,14 @@ volcano_plot <- function(lipid_data, pvalue_adjust = FALSE) { add_markers(color = ~LipidClass, size = 3) %>% layout(xaxis = list(zeroline = FALSE, - title = "log2(fold change"), + title = "log2(fold change)"), yaxis = list(title = y_title), shapes = list(vline(-1), vline(1), hline(-log10(0.05))), - legend = list(orientation = "h")) %>% + legend = list(orientation = "h"), + title = list(text = title, + x = 0)) %>% event_register(event = "plotly_click") return(p) diff --git a/man/volcano_plot.Rd b/man/volcano_plot.Rd index 189f798..bd8c0e6 100644 --- a/man/volcano_plot.Rd +++ b/man/volcano_plot.Rd @@ -4,12 +4,14 @@ \alias{volcano_plot} \title{Create volcano plot} \usage{ -volcano_plot(lipid_data, pvalue_adjust = FALSE) +volcano_plot(lipid_data, pvalue_adjust = FALSE, title = "") } \arguments{ \item{lipid_data}{tibble with all the lipid data and test data} \item{pvalue_adjust}{show the corrected p value, default is FALSE} + +\item{title}{title of the plot} } \value{ plotly object From 5225c3d174ab5af6c16d7bead2758c0eb39a3644 Mon Sep 17 00:00:00 2001 From: Rico Derks Date: Thu, 1 Jul 2021 14:49:52 +0200 Subject: [PATCH 7/7] bumped version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index b7a81b2..52d7a9a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: lipidomics Type: Package Title: Lipidomics workflow -Version: 0.5 +Version: 0.6.0 Author: Rico Derks Maintainer: Rico Derks Description: Lipidomics workflow to proces the result files (export alignment)