R/utils_summarize.R

Defines functions stats_plot_force stats_plot_time_off stats_plot_time_on summarize_trap_data

Documented in stats_plot_force stats_plot_time_off stats_plot_time_on summarize_trap_data

#' Summarize trap project
#'
#' @param f reactive folder list
#' @param hz sampling frequency
#' @param factor_order user supplied order of factors
#' @export
summarize_trap_data <- function(f, hz, factor_order, is_shiny = T){
  
  trap_selected_project <- f$project$path
  
  trap_data_paths <- 
    list_files(trap_selected_project,
               pattern = "options.csv",
               recursive = TRUE)
  
 if(is_shiny) setProgress(0.05, detail = 'Reading Data')
  
  trap_data <-
    lapply(trap_data_paths$path, data.table::fread, nrows = 1) %>% 
    data.table::rbindlist(fill = TRUE) %>% 
    dplyr::filter(report == "success" & review == TRUE & include == TRUE)
  
  
  event_file_paths <- 
    unique(trap_data[, .(project, conditions, date, obs)]) %>% 
    unite('path', c(project, conditions, date, obs), sep = "/") %>% 
    pull(path) %>% 
    paste0("~/lasertrapr/", ., '/measured-events.csv')
  
  event_files_filtered <- 
    lapply(event_file_paths, data.table::fread) %>%
    data.table::rbindlist(fill = TRUE) %>% 
    dplyr::filter(keep == TRUE & event_user_excluded == FALSE)
  
  if(!is.null(factor_order)){
  event_files_filtered$conditions <- 
    factor(event_files_filtered$conditions,
           levels = factor_order)
  }
  
  trap_data_filtered_paths <- str_replace(event_file_paths, "measured-events.csv", "trap-data.csv")
  
  get_time <- 
    lapply(trap_data_filtered_paths,
           function(x){ 
             if(is_shiny)   incProgress(0.0025)
              data.table::fread(x) } ) %>%
              data.table::rbindlist(fill = TRUE) %>% 
              .[, .N, by = conditions] %>% 
              dplyr::mutate(minutes = round((N/hz)/60, 2))
  
  summarize_trap <- event_files_filtered %>%
    dplyr::group_by(conditions) %>%
    dplyr::summarize(time_on_avg = mean(time_on_ms, na.rm = TRUE),
                     time_on_se = plotrix::std.error(time_on_ms, na.rm = TRUE),
                     time_on_sd = sd(time_on_ms, na.rm = TRUE),
                     time_on_median = median(time_on_ms, na.rm = TRUE),
                     time_off_avg = mean(time_off_ms, na.rm = TRUE),
                     time_off_se = plotrix::std.error(time_off_ms, na.rm = TRUE),
                     time_off_sd = sd(time_off_ms, na.rm = TRUE),
                     time_off_median = median(time_off_ms, na.rm = TRUE), 
                     displacement_avg = mean(displacement_nm, na.rm = TRUE),
                     displacement_se = plotrix::std.error(displacement_nm, na.rm = TRUE),
                     displacement_sd = sd(displacement_nm, na.rm = TRUE),
                     force_avg = mean(force, na.rm = TRUE),
                     force_se = plotrix::std.error(force, na.rm = TRUE),
                     force_sd = sd(force, na.rm = TRUE),
                     trap_stiffness = mean(trap_stiffness, na.rm = T),
                     myo_stiffness = mean(myo_stiffness, na.rm = T),
                     num_events = n()) %>%
    dplyr::right_join(get_time)# %>% 
  
  return(list(event_files_filtered = event_files_filtered,
               summary = summarize_trap,
               event_file_paths = event_file_paths))
}

#' Title
#'
#' @param event_files_filtered 
#' @param plot_colors 
stats_plot_step <- function(event_files_filtered, plot_colors){
  # event_files_filtered = rv$data$event_files_filtered
  # plot_colors = plot_colors
  step_histo <- ggplot(data = event_files_filtered,
                       aes(x = displacement_nm,
                           fill = conditions))+
    geom_histogram(aes(y = stat(density)),
                   binwidth = 3,
                   color = "black",
                   size = 0.9)+
    facet_wrap(~conditions, ncol = 2)+
    xlab("Step Size (nm)")+
    scale_y_continuous(expand = c(0,0))+
    scale_x_continuous(breaks = seq(-100, 100, by = 10))+
    scale_fill_manual(values = plot_colors)+
    #scale_fill_brewer(palette = "Dark2")+
    theme_linedraw(base_size = 11)+
    theme(panel.grid = element_blank(),
          legend.position = 'none')
  
  step_aov <-  event_files_filtered %>%
    rstatix::anova_test(displacement_nm ~ conditions)
  
  step_tukey <-  event_files_filtered %>%
    stats::aov(displacement_nm ~ conditions, data = .) %>%
    rstatix::tukey_hsd()
  
  step_is_sig <- which(step_tukey$p.adj < 0.05)
  
  step_table <- step_tukey %>%
    dplyr::select(group1, group2, p.adj) %>%
    dplyr::mutate_if(is.numeric, round, digits = 3) %>%
    ggpubr::ggtexttable (theme = ggpubr::ttheme(base_style = 'light', base_size = 9), rows = NULL)
  
  if(any(step_is_sig))  step_table %<>% ggpubr::table_cell_bg(row = step_is_sig + 1, column = 3, fill = '#ffa9a6')
  
  # step_tukey <- rstatix::add_xy_position(step_tukey,
  #                                        data = event_files_filtered,
  #                                        formula = displacement_nm ~ conditions)
  # 
  
  step_violin <- ggpubr::ggviolin(event_files_filtered,
                                  x = "conditions",
                                  y = "displacement_nm",
                                  fill = "conditions",
                                  palette = plot_colors,
                                  add = "boxplot",
                                  xlab = '',
                                  ylab = 'Displacement (nm)',
                                  ggtheme = ggpubr::theme_pubr(base_size=11))+
    ggpubr::stat_pvalue_manual(step_tukey, bracket.nudge.y = seq(10, by = 4, length.out = nrow(step_tukey)))+
    scale_y_continuous(breaks = seq(-50, 1000, by = 15))+
    rotate_x_text(angle = 45)+
    # rotate()+
    #yscale('log10', .format = T)+
    labs(
      subtitle = rstatix::get_test_label(step_aov, detailed = TRUE),
      caption = rstatix::get_pwc_label(step_tukey)
    ) +
    theme(legend.position = 'none')
  
  
  
  #  step_top <- cowplot::plot_grid(step_violin, step_histo, align = "hv", axis = "bt", rel_widths = c(1,1))
  step_side <- cowplot::plot_grid(step_histo, step_table, ncol = 1, nrow = 2,  rel_heights = c(0.75,0.25))
  step_main <-  cowplot::plot_grid(step_violin, step_side, ncol = 2, nrow = 1,  rel_heights = c(1, 1))
  # now add the title
  step_title <- cowplot::ggdraw() +
    cowplot::draw_label(
      "Displacement Distributions, ANOVA, & Tukey Table",
      fontface = 'bold',
      size = 18,
      x = 0.5
    )
  
  cowplot::plot_grid(step_title, 
                     step_main, #step_table,
                     ncol = 1, 
                     nrow =2,
                     rel_heights = c(0.05, 1))
  
}

#' Title
#'
#' @param event_files_filtered 
#' @param plot_colors 
#' @param p_adj 

stats_plot_time_on <- function(event_files_filtered, plot_colors, p_adj){
  
  ton_histo <- ggplot(data = event_files_filtered,
                      aes(x = time_on_ms,
                          fill = conditions))+
    geom_histogram(aes(y = stat(density)),
                   binwidth = 20,
                   color = "black",
                   size = 0.9)+
    facet_wrap(~conditions, ncol = 2)+
    xlab("Time on (ms)")+
    scale_y_continuous(expand = c(0,0))+
    #scale_x_continuous(breaks = seq(0, 6000, by = 100))+
    coord_cartesian(c(0, 1000))+
    scale_fill_manual(values = plot_colors)+
    #scale_fill_brewer(palette = "Dark2")+
    theme_linedraw(base_size = 11)+
    theme(panel.grid = element_blank(),
          legend.position = "none")
  
  ton_kruskal <-  event_files_filtered %>%
    rstatix::kruskal_test(time_on_ms ~ conditions)
  
  ton_wilcox <-   event_files_filtered %>%
    rstatix::wilcox_test(time_on_ms ~ conditions,  p.adjust.method = p_adj)
  
  ton_is_sig <- which(ton_wilcox$p.adj < 0.05)
  
  ton_table <- ton_wilcox %>%
    dplyr::select(group1, group2, p.adj) %>%
    dplyr::mutate_if(is.numeric, round, digits = 3) %>%
    ggpubr::ggtexttable (theme = ggpubr::ttheme(base_style = 'light', base_size = 9), rows = NULL)
  
  if(any(ton_is_sig)) ton_table %<>% ggpubr::table_cell_bg(row = ton_is_sig + 1, column = 3, fill = '#ffa9a6')
  
  ton_wilcox %<>% rstatix::add_xy_position(x = "conditions", y.trans = log10)
  
  ton_violin <- ggpubr::ggviolin(event_files_filtered,
                                 x = "conditions",
                                 y = "time_on_ms",
                                 fill = "conditions",
                                 palette = plot_colors,
                                 add = "boxplot",
                                 xlab = '',
                                 ylab = 'Time On (ms)',
                                 ggtheme = ggpubr::theme_pubr(base_size=11))+
    ggpubr::stat_pvalue_manual(ton_wilcox, bracket.nudge.y = seq(0.5, by = 0.3, length.out = nrow(ton_wilcox)))+
    yscale('log10', .format = T)+
    rotate_x_text(angle=45)+
    labs(
      subtitle = rstatix::get_test_label(ton_kruskal, detailed = TRUE),
      caption = rstatix::get_pwc_label(ton_wilcox)
    ) +
    theme(legend.position = 'none')
  
  ton_side <- cowplot::plot_grid(ton_histo, ton_table, ncol = 1, nrow = 2,  rel_heights = c(0.75,0.25))
  ton_main <-  cowplot::plot_grid(ton_violin, ton_side, ncol = 2, nrow = 1)
  
  # now add the title
  ton_title <- cowplot::ggdraw() +
    cowplot::draw_label(
      "Time On Distributions, Kruskal-Wallis, & Pairwise Wilcox Table",
      fontface = 'bold',
      size = 18,
      x = 0.5
    )
  
  cowplot::plot_grid(ton_title,
                     ton_main, #ton_table,
                     ncol = 1,
                     nrow =2,
                     rel_heights = c(0.05, 1))
}

#' Title
#'
#' @param event_files_filtered 
#' @param plot_colors 

stats_plot_time_off <- function(event_files_filtered, plot_colors){
  toff_histo <- 
    ggplot(data = event_files_filtered,
           aes(x = time_off_ms,
               fill = conditions))+
    geom_histogram(aes(y = stat(density)),
                   binwidth = 100,
                   color = "black",
                   size = 0.9)+
    facet_wrap(~conditions, ncol=2)+
    xlab("Time Off (ms)")+
    scale_y_continuous(expand = c(0,0))+
    #scale_x_continuous(breaks = seq(0, 20000, by = 500))+
    coord_cartesian(c(0, 2000))+
    scale_fill_manual(values = plot_colors)+
    #scale_fill_brewer(palette = "Dark2")+
    theme_linedraw(base_size = 11)+
    theme(panel.grid = element_blank(),
          legend.position = "none")
  
  
  toff_kruskal <-  
    event_files_filtered %>%
    rstatix::kruskal_test(time_off_ms ~ conditions)
  
  toff_wilcox <-   
    event_files_filtered %>%
    wilcox_test(time_off_ms ~ conditions,  p.adjust.method = "holm")
  
  toff_is_sig <- which(toff_wilcox$p.adj < 0.05)
  
  toff_table <- 
    toff_wilcox %>%
    dplyr::select(group1, group2, p.adj) %>%
    dplyr::mutate_if(is.numeric, round, digits = 3) %>%
    ggpubr::ggtexttable (theme = ttheme(base_style = 'light', base_size = 9), rows = NULL)
  
  if(any(toff_is_sig)) toff_table %<>%  ggpubr::table_cell_bg(row = toff_is_sig + 1, column = 3, fill = '#ffa9a6')
  
  toff_wilcox %<>% add_xy_position(x = "conditions", y.trans = log10)
  
  toff_violin <- 
    ggviolin(event_files_filtered,
             x = "conditions",
             y = "time_off_ms",
             fill = "conditions",
             palette = plot_colors,
             add = "boxplot",
             xlab = '',
             ylab = 'Time Off (ms)',
             ggtheme = theme_pubr(base_size=11))+
    stat_pvalue_manual(toff_wilcox, bracket.nudge.y = seq(0.5, by = 0.3, length.out = nrow(toff_wilcox)))+
    yscale('log10', .format = T)+
    rotate_x_text(angle=45)+
    labs(
      subtitle = get_test_label(toff_kruskal, detailed = TRUE),
      caption = get_pwc_label(toff_wilcox)
    ) +
    theme(legend.position = 'none')
  
  
  toff_side <- cowplot::plot_grid(toff_histo, toff_table, ncol = 1, nrow = 2,  rel_heights = c(0.75,0.25))
  toff_main <-  cowplot::plot_grid(toff_violin, toff_side, ncol = 2, nrow = 1)
  
  # now add the title
  toff_title <- ggdraw() +
    draw_label(
      "Time Off Distributions, Kruskal-Wallis, & Pairwise Wilcox Table",
      fontface = 'bold',
      size = 18,
      x = 0.5
    )
  
  plot_grid(toff_title, 
             toff_main, #toff_table,
             ncol = 1,
             nrow =3,
             rel_heights = c(0.05, 1))
}

#' Title
#'
#' @param event_files_filtered 
#' @param plot_colors 

stats_plot_force <- function(event_files_filtered, plot_colors){
  
  force_histo <- ggplot(data = event_files_filtered,
                        aes(x = force,
                            fill = conditions))+
    geom_histogram(aes(y = stat(density)),
                   binwidth = 0.25,
                   color = "black",
                   size = 0.9)+
    facet_wrap(~conditions, ncol=2)+
    xlab("Force (pN)")+
    scale_y_continuous(expand = c(0,0))+
    #scale_x_continuous(breaks = seq(0, 20000, by = 500))+
    scale_fill_manual(values = plot_colors)+
    #scale_fill_brewer(palette = "Dark2")+
    theme_linedraw(base_size = 11)+
    theme(panel.grid = element_blank(),
          legend.position = "none")
  
  
  force_aov <-  event_files_filtered %>%
    anova_test(force ~ conditions)
  
  force_tukey <-  event_files_filtered %>%
    aov(force~conditions, data = .) %>%
    rstatix::tukey_hsd()
  
  force_is_sig <- which(force_tukey$p.adj < 0.05)
  
  force_table <- force_tukey %>%
    dplyr::select(group1, group2, p.adj) %>%
    dplyr::mutate_if(is.numeric, round, digits = 3) %>%
    ggpubr::ggtexttable (theme = ggpubr::ttheme(base_style = 'light', base_size = 9), rows = NULL)
  
  if(any(force_is_sig)) force_table %<>% ggpubr::table_cell_bg(row = force_is_sig + 1, column = 3, fill = '#ffa9a6')
  
  force_tukey <- add_xy_position(force_tukey, data = event_files_filtered, formula = force ~ conditions)
  
  force_violin <- ggviolin(event_files_filtered,
                           x = "conditions",
                           y = "force",
                           fill = "conditions",
                           palette = plot_colors,
                           add = "boxplot",
                           xlab = '',
                           ylab = 'Time Off (ms)',
                           ggtheme = theme_pubr(base_size=11))+
    stat_pvalue_manual(force_tukey, bracket.nudge.y = seq(0.5, by = 0.3, length.out = nrow(force_tukey)))+
    scale_y_continuous(breaks = seq(-3, 5, by = 1))+
    rotate_x_text(angle=45)+
    #yscale('log10', .format = T)+
    labs(
      subtitle = get_test_label(force_aov, detailed = TRUE),
      caption = get_pwc_label(force_tukey)
    ) +
    theme(legend.position = 'none')
  
  force_side <- cowplot::plot_grid(force_histo, force_table, ncol = 1, nrow = 2,  rel_heights = c(0.75,0.25))
  force_main <-  cowplot::plot_grid(force_violin, force_side, ncol = 2, nrow = 1)
  
  force_title <- cowplot::ggdraw() +
    cowplot::draw_label(
      "Force Distributions, ANOVA, & Tukey Table",
      fontface = 'bold',
      size = 18,
      x = 0.5
    )
  
  cowplot::plot_grid(force_title,
                     force_main, #force_table,
                     ncol = 1,
                     nrow =2,
                     rel_heights = c(0.05, 1))
}



#' Time On ECDF
#'
#' @param event_files_filtered event data
#' @param plot_colors user colors
#'
#' @return a ggplot
time_on_ecdf <- function(event_files_filtered, plot_colors){
  cum_freq_curves <-
    event_files_filtered %>% 
    split(.$conditions) %>% 
    imap_dfr(~ .x %>% dplyr::summarize(conditions = .y,
                                       time_on = unique(time_on_ms), 
                                       ecdf = ecdf(time_on_ms)(unique(time_on_ms)))) 
 
  
    ggplot(cum_freq_curves, aes(time_on, ecdf)) + 
      geom_step(aes(color = conditions), size = 0.75)+
      #coord_cartesian(c(0, 1000))+
      #facet_grid(~myo)+
      ggtitle("Time On Cumulative Event Distributions")+
      ylab("Cumulative Distribution")+
      xlab("Time (ms)")+
      scale_color_manual(values = plot_colors)+
      theme_cowplot()+
      theme(panel.grid = element_blank(),
            legend.position = "none")
  

}
#' Time Off ECDF
#'
#' @param event_files_filtered event data
#' @param plot_colors user colors
#'
#' @return a ggplot
time_off_ecdf <- function(event_files_filtered, plot_colors){
  cum_freq_curves <-
    event_files_filtered %>% 
    split(.$conditions) %>% 
    imap_dfr(~ .x %>% dplyr::summarize(conditions = .y,
                                       time_off = unique(time_off_ms), 
                                       ecdf = ecdf(time_off_ms)(unique(time_off_ms)))) 
  
  
  ggplot(cum_freq_curves, aes(time_off, ecdf)) + 
    geom_step(aes(color = conditions), size = 0.75)+
    #coord_cartesian(c(0, 1000))+
    #facet_grid(~myo)+
    ggtitle("Time Off Cumulative Event Distributions")+
    ylab("Cumulative Distribution")+
    xlab("Time (ms)")+
    scale_color_manual(values = plot_colors)+
    theme_cowplot()+
    theme(panel.grid = element_blank())
  
  
}
#' Time On Survival Analysis
#'
#' @param event_files_filtered 
#' @param plot_colors 
#' @return
#' @export
#'
#' @noRd
#' @import survminer survival
survival_analysis <- function(event_files_filtered, plot_colors){
  
  fit <- survival::survfit(survival::Surv(time_on_ms) ~ conditions, data = event_files_filtered)
  
  ggsurvival <-  
    survminer::ggsurvplot(fit,
                          data = event_files_filtered,
                          ggtheme = theme_linedraw(base_size = 12), # Change ggplot2 theme 
                          palette = plot_colors,
                          xlim = c(0, 1000),
                          break.time.by = 150,
                          fun = 'event'
                           # pval = T
  )
  
  all_surv <- ggsurvival$plot +
    xlab('Time (ms)')+
    xlab("")+
    theme(panel.grid = element_blank(),
          legend.title = element_blank(),
    )
  
  facet_surv <- ggsurvival$plot +
    facet_wrap(~conditions)+
    xlab('Time (ms)')+
    theme(panel.grid = element_blank(),
          legend.title = element_blank(),
          legend.position = 'none')
  
  ton_survival_plots <- all_surv/facet_surv
  ton_log_rank_test <- survival::survdiff(survival::Surv(time_on_ms) ~ conditions, data = event_files_filtered)
  ton_cox <- survival::coxph(survival::Surv(time_on_ms) ~ conditions, data = event_files_filtered)
  
  list(ton_survival_plots = ton_survival_plots,
       ton_log_rank_test = ton_log_rank_test,
       ton_cox = ton_cox)
       
}

#' Title
#'
#' @param event_file_paths 
#' @param factor_order 
#' @param plot_colors 
stats_plot_event_frequency <- function(event_file_paths, factor_order, plot_colors){
  ef_paths <- str_replace(event_file_paths, "measured-events", "event-frequency")
  
  ef_read <- function(path){
    data.table::fread(path) %>%
      dplyr::mutate(path = path)
  }
  ef <- map_df(ef_paths, ef_read) %>%
    separate(path, c('root', 'id_info'), sep = '/project_') %>% 
    separate(id_info, c("project", 'conditions', 'date', 'obs', 'file'), sep = '/')
  
  ef$conditions <- factor(ef$conditions, levels = factor_order)
  
  ef_kruskal <- ef %>% rstatix::kruskal_test(freq_start ~ conditions)
  
  ef_wilcox <-  ef %>%
    wilcox_test(freq_start ~ conditions, p.adjust.method = "holm")
  
  ef_is_sig <- which(ef_wilcox$p.adj < 0.05)
  
  ef_table <- ef_wilcox %>%
    dplyr::select(group1, group2, p.adj) %>%
    dplyr::mutate_if(is.numeric, round, digits = 3) %>%
    ggtexttable (theme = ttheme(base_style = 'light', base_size = 10), rows = NULL) 
  
  if(any(ef_is_sig)) ef_table %<>%table_cell_bg(row = ef_is_sig + 1, column = 3, fill = '#ffa9a6')
  
  ef_wilcox %<>% add_xy_position(x = "conditions")
  
  ef_plot <- ggplot(data = ef, aes(conditions, freq_start))+
    geom_jitter(aes(color = conditions), size = 2, alpha = 0.8)+
    stat_pvalue_manual(ef_wilcox) +
    stat_summary(
      fun.data = mean_sdl,  fun.args = list(mult = 1),
      geom = "pointrange", color = "black"
    )+
    scale_color_manual(values = plot_colors)+
    scale_y_continuous(breaks = 0:15)+
    xlab('')+
    ylab('Events/sec')+
    labs(
      title = 'Event Frequency',
      subtitle = get_test_label(ef_kruskal, detailed = TRUE),
      caption = get_pwc_label(ef_wilcox)
    )+
    theme_cowplot(font_size = 11)+
    theme(legend.position = 'none')#+

  ef_plot <- plot_grid(ef_plot, ef_table, ncol = 2, nrow = 1, rel_widths = c(2, 1))
  
  #no zero
  ef_no0 <- dplyr::filter(ef, freq_start != 0 )
  ef_kruskal_no_zero <- ef_no0 %>%
    kruskal_test(freq_start ~ conditions)

  ef_wilcox_no_zero <-  ef_no0 %>%
    wilcox_test(freq_start ~ conditions, p.adjust.method = "holm")

  ef_is_sig_no_zero <- which(ef_wilcox_no_zero$p.adj < 0.05)

  ef_table_no_zero <- ef_wilcox_no_zero %>%
    dplyr::select(group1, group2, p.adj) %>%
    dplyr::mutate_if(is.numeric, round, digits = 3) %>%
    ggpubr::ggtexttable(theme = ggpubr::ttheme(base_style = 'light', base_size = 10), rows = NULL) 
  
  if(any(ef_is_sig_no_zero)) ef_table_no_zero %<>% table_cell_bg(row = ef_is_sig_no_zero + 1, column = 3, fill = '#ffa9a6')

  ef_wilcox_no_zero %<>% add_xy_position(x = "conditions")

  ef_plot_no_zero <- ggplot(data = ef_no0, aes(conditions, freq_start))+
    geom_jitter(aes(color = conditions), size = 2, alpha = 0.8)+
    stat_pvalue_manual(ef_wilcox_no_zero) +
    stat_summary(
      fun.data = mean_sdl,  fun.args = list(mult = 1),
      geom = "pointrange", color = "black"
    )+
    scale_color_manual(values = plot_colors)+
    scale_y_continuous(breaks = 0:15)+
    xlab('')+
    ylab('Events/sec')+
    labs(
      title = 'Event Frequency (zero counts eliminated)',
      subtitle = get_test_label(ef_kruskal_no_zero, detailed = TRUE),
      caption = get_pwc_label(ef_wilcox_no_zero)
    )+
    theme_cowplot(font_size = 11)+
    theme(legend.position = 'none')#+

  ef_plot_no_zero <- plot_grid(ef_plot_no_zero, ef_table_no_zero, ncol = 2, nrow = 1, rel_widths = c(2, 1))
  
  list(ef_plot = ef_plot,
       ef_plot_no_zero = ef_plot_no_zero)
}

#' Title
#'
#' @param event_files_filtered 
#' @param plot_colors 
correlations <- function(event_files_filtered, plot_colors){
  step_vs_on <- ggscatter(event_files_filtered,
                          y = "time_on_ms", 
                          x = "displacement_nm",
                          color = 'conditions',
                          palette = plot_colors, 
                          alpha = 0.25,
                          add = "reg.line",  # Add regressin line
                          add.params = list(color = "conditions", fill = "grey", size = 1.25), # Customize reg. line
                          conf.int = TRUE,
  )+
    facet_wrap(~conditions, ncol = 1)+
    scale_y_log10()+
    scale_x_continuous(breaks = seq(-100, 100, by = 10))+
    ylab('Time On (ms)')+
    xlab('Displacement (nm)')+
    ggtitle('Displacement vs Time On')+
    stat_cor(method = "pearson", label.x.npc = "left", label.y.npc = "top")+
    theme_linedraw(base_size = 12)+
    theme(panel.grid  = element_blank(),
          legend.position = 'none')
  
  
  
  on_vs_off <- ggscatter(event_files_filtered,
                         y = "time_on_ms", 
                         x = "time_off_ms",
                         color = 'conditions',
                         palette = plot_colors, 
                         alpha = 0.25,
                         #  facet.by = "conditions",
                         add = "reg.line",  # Add regressin line
                         add.params = list(color = "conditions", fill = "grey", size = 1.25), # Customize reg. line
                         conf.int = TRUE)+
                         facet_wrap(~conditions, ncol = 1)+
                         scale_y_log10()+
                         scale_x_log10()+
                         ylab('Time On (ms)')+
                         xlab('Time Off (ms)')+
                         ggtitle('Time On vs Time Off')+
                         stat_cor(method = "pearson", label.x.npc = "left", label.y.npc = "top")+
                         theme_linedraw(base_size = 12)+
                         theme(panel.grid  = element_blank(),
                                legend.position = 'none')
  
  
  step_vs_off <- ggscatter(event_files_filtered,
                           x = "displacement_nm", 
                           y = "time_off_ms",
                           color = 'conditions',
                           alpha = 0.25,
                           palette = plot_colors, 
                           #facet.by = "conditions",
                           add = "reg.line",  # Add regressin line
                           add.params = list(color = "conditions", fill = "grey", size = 1.25), # Customize reg. line
                           conf.int = TRUE,
  )+
    facet_wrap(~conditions, ncol = 1)+
    scale_y_log10()+
    scale_x_continuous(breaks = seq(-100, 100, by = 10))+
    xlab('Displacement (nm)')+
    ylab('Time Off (ms)')+
    ggtitle('Displacement vs Time Off')+
    stat_cor(method = "pearson", label.x.npc = "left", label.y.npc = "top")+
    theme_linedraw(base_size = 12)+
    theme(panel.grid  = element_blank(),
          legend.position = 'none')
  
  
  ktrap_vs_kmyo <- ggscatter(event_files_filtered,
                             x = "trap_stiffness", 
                             y = "myo_stiffness",
                             color = 'conditions',
                             alpha = 0.25,
                             palette = plot_colors,
                             # facet.by = "conditions",
                             add = "reg.line",  # Add regressin line
                             add.params = list(color = "conditions", fill = "grey", size = 1.25), # Customize reg. line
                             conf.int = TRUE,
  )+
    facet_wrap(~conditions, ncol = 1)+
    xlab(expression(k[trap]))+
    ylab(expression(k[myo]))+
    ggtitle('ktrap vs kmyo')+
    scale_x_log10()+
    stat_cor(method = "pearson", label.x.npc = "left", label.y.npc = "top")+
    theme_linedraw(base_size = 12)+
    theme(panel.grid  = element_blank(),
          legend.position = 'none')
  
  step_vs_kmyo <- ggscatter(event_files_filtered,
                            y = "displacement_nm", 
                            x = "myo_stiffness",
                            color = 'conditions',
                            alpha = 0.25,
                            palette = plot_colors, 
                            # facet.by = "conditions",
                            add = "reg.line",  # Add regressin line
                            add.params = list(color = "conditions", fill = "grey", size = 1.25), # Customize reg. line
                            conf.int = TRUE,
  )+
    facet_wrap(~conditions, ncol = 1)+
    ylab('Displacement (nm)')+
    xlab(expression(k[myo]))+
    scale_y_continuous(breaks = seq(-100, 100, by = 10))+
    ggtitle('Displacement vs kmyo')+
    #scale_x_log10()+
    stat_cor(method = "pearson", label.x.npc = "left", label.y.npc = "top")+
    theme_linedraw(base_size = 12)+
    theme(panel.grid  = element_blank(),
          legend.position = 'none')
  
  
  plot_grid(step_vs_on, on_vs_off, step_vs_kmyo, ktrap_vs_kmyo, step_vs_off, nrow = 1) 
}

#' Title
#'
#' @param event_files_filtered 
#' @param plot_colors 
stiffness <- function(event_files_filtered, plot_colors){
  
  stiff_data <- event_files_filtered %>% 
    dplyr::select(myo_stiffness, trap_stiffness, conditions) %>% 
    tidyr::pivot_longer(cols= c(myo_stiffness, trap_stiffness), names_to = 'stiffness')
  
  stiff_data$stiffness <- factor(stiff_data$stiffness, levels = c('trap_stiffness', 'myo_stiffness'))
  mean_trap <- stiff_data %>% 
    group_by(conditions, stiffness) %>% 
    dplyr::summarize(mean_val = mean(value, na.rm = T),
                     sd = sd(value, na.rm = T)) %>% 
    mutate(s2n = ifelse(stiffness == 'trap_stiffness', '1', round(max(mean_val)/min(mean_val), 2)),
           label = paste0("\u03bc = ", round(mean_val, 3)))
  
  
  
  g <- ggplot()+
    geom_density(data = filter(stiff_data, stiffness == "trap_stiffness"), 
                 aes(value), fill = "black", alpha = 0.4,
                 #color = 'blac k'
    )+
    geom_density(data = filter(stiff_data, stiffness == "myo_stiffness"), 
                 aes(value, fill = conditions), alpha = 0.4,
                 #color = 'blac k'
    )+
    # geom_vline(data = mean_trap, aes(xintercept=mean_val))+
    scale_y_continuous(breaks = NULL, expand = expansion(c(0, 0)))+
    # scale_x_continuous(breaks = seq(0, 5, by = 0.05))+
    scale_fill_manual(values = plot_colors, 
                      name = '',
                     # labels = c(bquote('k'[trap]), bquote('k'[myo]))
                     )+
    xlab('Stiffness (pN/nm)')+
    facet_wrap(~conditions, nrow=1)+
    coord_cartesian(c(0, 0.5))+
    theme_linedraw()+
    theme(legend.position = 'none',
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
         # strip.background = element_rect(fill = "white"),
          strip.text = element_text(face = "bold"))
  
  #g
   gb <- ggplot_build(g)
  # 
  # #get_y <- gb$data[[1]] %>% filter(abs(x - mean_val) == min(abs(x - mean_val)) %>% pull(y)
  # # mean_trap %<>%
  # #   arrange(conditions) %>%
  # #   mutate(gg_build = split(gb$data[[1]], gb$data[[1]]$PANEL))
  # 
  # 
  # by_conditions <- split(gb$data[[2]], gb$data[[2]]$PANEL)
  # by_fill <- map(by_conditions, ~split(., .$fill))
  # names(by_fill) <- levels(event_files_filtered$conditions)
  # gg_df <- unlist(by_fill, recursive = F)
  # myo_mean <- filter(mean_trap, stiffness == "myo_stiffness")
  # y_position_df <- map2(gg_df, myo_mean$mean_val,  ~dplyr::filter(.x, abs(x - .y) == min(abs(x - .y))))
  # get_y <- map_dbl(y_position_df, `[[`, 'y')
  # myo_mean$y  <- get_y
  # 
  # trap_by_conditions <- split(gb$data[[1]], gb$data[[1]]$PANEL)
  # trap_by_fill <- map(by_conditions, ~split(., .$fill))
  # names(trap_by_fill) <- levels(rv$data$event_files_filtered$conditions)
  # trap_gg_df <- unlist(trap_by_fill, recursive = F)
  # trap_mean <- filter(mean_trap, stiffness == "trap_stiffness")
  # trap_y_position_df <- map2(trap_gg_df, trap_mean$mean_val,  ~dplyr::filter(.x, abs(x - .y) == min(abs(x - .y))))
  # trap_get_y <- map_dbl(trap_y_position_df, `[[`, 'y')
  # trap_mean$y  <- get_y
  # 
  # 
  # mean_trap <- rbind(myo_mean, trap_mean)
  #
  #
  #

  # g <- g +
  #   geom_point(data = mean_trap, aes(x = mean_val, y = y), size = 2)+
  #   geom_segment(data = mean_trap, aes(x = mean_val, y = y, xend = mean_val+0.05, yend = y+2))+
  #   geom_label(data = mean_trap, aes(x = mean_val+0.05, y = y+2, label = label))
    # geom_col(data = mean_trap, aes(conditions, s2n, fill = stiffness))+
    #coord_cartesian(c(0, 0.3))

  # geom_point(data = blue_data, aes(x = x, y = density), size = 2)+
  # geom_segment(data = blue_data, aes(x = x, y = density, xend = x+0.05, yend = density+2))+
  # geom_label(data = blue_data, aes(x = x+0.05, y = density+2, label = paste0("\u03bc = ", blue_mean) ))
  annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data)
  {
    layer(data = data, stat = StatIdentity, position = PositionIdentity,
          geom = ggplot2:::GeomCustomAnn,
          inherit.aes = TRUE, params = list(grob = grob,
                                            xmin = xmin, xmax = xmax,
                                            ymin = ymin, ymax = ymax))
  }

  get_inset <- function(df, colorz){
    ggplot(df, aes(conditions, s2n, fill = stiffness))+
      geom_col(alpha = 0.5)+
      geom_text(aes(label = s2n), position = position_stack(vjust = 0.5), size = 4)+
      scale_fill_manual(values = c('grey30', colorz))+
      # scale_y_continuous(position = 'right')+
      coord_flip()+
      ggtitle('Signal-to-noise')+
      theme_nothing()+
      theme(legend.position = 'none',
            # axis.title.x  = element_text(size = 14, face = 'bold'),
            #axis.text.y = element_text(size = 8),
            plot.title = element_text(size = 10, face = 'bold'))

  }

  max_y <- max(gb[["data"]][[1]][["density"]])
  insets <- mean_trap %>%
    split(f = .$conditions) %>%
    purrr::map2(., plot_colors, ~annotation_custom2(
      grob = ggplotGrob(get_inset(.x, .y)),
      data = data.frame(conditions=unique(.$conditions)),
      ymin = max_y-25, ymax=max_y-5, xmin=0.3, xmax=0.5
      )
    )

  g <- g+insets
  
  
  trap_mean <-
    mean_trap %>% 
    filter(stiffness == "trap_stiffness") %>% 
    mutate(color = "black")
  
  myo_mean <-
    mean_trap %>% 
    filter(stiffness == "myo_stiffness") 
  
  myo_mean$color <- plot_colors
  
  mean_trap <- 
    rbind(trap_mean, myo_mean) %>% 
    arrange(conditions)
    
                    
  table <- 
    mean_trap %>% 
    select("Conditions" = conditions, 
           "Type" = stiffness, 
           "Mean" = mean_val,
           "SD" = sd) %>% 
    mutate_if(is.numeric, round, digits = 4) %>% 
    ggtexttable(theme = ttheme(
      colnames.style = colnames_style(fill = "black", color = "white"),
      tbody.style = tbody_style(fill = alpha(mean_trap$color, 0.5))))
  
  plot_grid(g, table, nrow = 2, rel_heights = c(0.6, 0.4))

}


#' Contruct Cumulative Frequency Curves
#'
#' @param data df
#' @param type time_on_ms or time_off_ms
#'
#' @return a list with a df and plot
#' @export
#'
#' @examples plot_cum_freq_curves(event_files_filtered, time_on_ms)
plot_cum_freq_curves <- function(data, type){
  
  cum_freq_curves <- 
    data %>% 
    split(.$conditions) %>% 
    imap_dfr(~ .x %>% dplyr::summarize(conditions = .y,
                                       time = unique({{type}}), 
                                       ecdf = ecdf({{type}})(unique({{type}})))) 
  
  # You can reorder your conditiond factors like this
  # cum_freq_curves$conditions <- factor(cum_freq_curves$conditions, 
  #                                      levels =  c("myoV-WT_pH-7.0_0mM-Pi",
  #                                                  "myoV-WT_pH-7.0_30mM-Pi",
  #                                                  "myoV-S217A_pH-7.0_0mM-Pi",
  #                                                  "myoV-S217A_pH-7.0_30mM-Pi"))
  
  p <- ggplot(cum_freq_curves, aes(time, ecdf)) + 
    geom_step(aes(color = conditions))+
    #coord_cartesian(c(0, 1000))+
    #facet_grid(~myo)+
    #ggtitle("Time On")+
    ylab("ECDF")+
    xlab("Time (ms)")+
    #scale_color_manual(values = plot_colors)+
    #coord_cartesian(c(0, 1000))+
    theme_cowplot(font_size = 9)+
    theme(panel.grid = element_blank(),
          legend.position = 'none')
  
  return(list(data = cum_freq_curves,
              plot = p))
}
brentscott93/lasertrapr documentation built on Jan. 15, 2022, 8:21 p.m.