Merge pull request #98 from hadexversum/dev_1.5

Dev 1.5

NAMESPACE

@@ -9,18 +9,19 @@ export(calculate_diff_uptake)
 export(calculate_exp_masses)
 export(calculate_exp_masses_per_replicate)
 export(calculate_kinetics)
+export(calculate_p_value)
 export(calculate_peptide_kinetics)
 export(calculate_state_uptake)
 export(create_control_dataset)
 export(create_diff_uptake_dataset)
 export(create_kinetic_dataset)
 export(create_overlap_distribution_dataset)
+export(create_p_diff_uptake_dataset)
 export(create_quality_control_dataset)
 export(create_replicate_dataset)
 export(create_state_comparison_dataset)
 export(create_state_uptake_dataset)
 export(create_uptake_dataset)
-export(create_volcano_dataset)
 export(get_protein_coverage)
 export(get_protein_redundancy)
 export(plot_amino_distribution)
@@ -30,7 +31,9 @@ export(plot_coverage)
 export(plot_differential)
 export(plot_differential_butterfly)
 export(plot_differential_chiclet)
+export(plot_differential_uptake_curve)
 export(plot_kinetics)
+export(plot_manhattan)
 export(plot_overlap)
 export(plot_overlap_distribution)
 export(plot_peptide_charge_measurement)

R/calculations.R

@@ -42,7 +42,7 @@
 calculate_state_uptake <- function(dat,
                                    protein = unique(dat[["Protein"]])[1], 
                                    state = unique(dat[["State"]])[1], 
-                                   time_0 = min(dat[["Exposure"]]),
+                                   time_0 = min(dat[dat[["Exposure"]]>0, ][["Exposure"]]),
                                    time_t = unique(dat[["Exposure"]])[3], 
                                    time_100 = max(dat[["Exposure"]]),
                                    deut_part = 0.9){

R/differential_plots.R

@@ -379,107 +379,3 @@ plot_differential_chiclet <- function(diff_uptake_dat,
 
 }
 
-#' Volcano plot
-#'
-#' @importFrom ggplot2 coord_cartesian
-#'
-#' @param vol_data data produced by the \code{\link{create_volcano_dataset}}
-#' function.
-#' @param state_1 biological state for chosen protein. It is used in the title.
-#' @param state_2 biological state for chosen protein. It is used in the title.
-#' @param adjust_axes logical, indicating if the X-axis is symmetrical in
-#' relation to 0.
-#' @param show_confidence_limits logical, indicates if the hybrid testing
-#' confidence intervals are shown.
-#' @param confidence_level confidence level for the test, from range [0, 1]. It
-#' should be the same as used in \code{\link{create_volcano_dataset}} function.
-#' @param color_times logical, indicating if different time points on the plot
-#' are distinguishable by color. 
-#'
-#' @details The data produced by \code{\link{create_volcano_dataset}} are plotted
-#' in the form of a volcano plot. The generation of the data is described in the documentation
-#' of \code{\link{create_volcano_dataset}} function. The confidence limit on
-#' P-value is calculated based on the confidence level. The confidence limit on deuterium
-#' uptake difference is calculated using the Houde test for the time point of measurement
-#' from the provided data. The confidence limits are indicated by the red dotted
-#' lines. The points above confidence limits (upper right and left corner) are
-#' statistically significant in hybrid testing.
-#' This plot is visible in GUI.
-#'
-#' @references Hageman, T. S. & Weis, D. D. Reliable Identification of Significant
-#' Differences in Differential Hydrogen Exchange-Mass Spectrometry Measurements
-#' Using a Hybrid Significance Testing Approach. Anal Chem 91, 8008–8016 (2019).
-#' @references Houde, D., Berkowitz, S.A., and Engen, J.R. (2011).
-#' The Utility of Hydrogen/Deuterium Exchange Mass Spectrometry in
-#' Biopharmaceutical Comparability Studies. J Pharm Sci 100, 2071–2086.
-#'
-#' @return a \code{\link{ggplot}} object.
-#'
-#' @seealso
-#' \code{\link{create_volcano_dataset}}
-#' \code{\link{show_volcano_data}}
-#'
-#' @examples
-#' dat <-  read_hdx(system.file(package = "HaDeX", "HaDeX/data/KD_180110_CD160_HVEM.csv"))
-#' vol_dat <- create_volcano_dataset(dat)
-#' plot_volcano(vol_dat, show_confidence_limits = T)
-#'
-#' @export plot_volcano
-
-plot_volcano <- function(vol_data,
-                         state_1 = "",
-                         state_2 = "",
-                         adjust_axes = TRUE,
-                         show_confidence_limits = FALSE,
-                         confidence_level = 0.98,
-                         color_times = TRUE) {
-
-
-
-  if (color_times){
-
-    volcano_plot <- ggplot(vol_data, aes(x = D_diff, y = log_p_value)) +
-      geom_point(aes(color = as.factor(Exposure))) +
-      geom_errorbar(aes(xmin = D_diff - Uncertainty, xmax = D_diff + Uncertainty), alpha = 0.2) +
-      labs(title = paste0("Volcano Plot ", state_1, " " , state_2),
-           x = "Mass difference [Da]",
-           y = "-log(P value)") +
-      labs(color = "Exposure") 
-
-  } else {
-
-    volcano_plot <- ggplot(vol_data, aes(x = D_diff, y = log_p_value)) +
-      geom_point() +
-      geom_errorbar(aes(xmin = D_diff - Uncertainty, xmax = D_diff + Uncertainty), alpha = 0.2) +
-      labs(title = paste0("Volcano Plot ", state_1, " " , state_2),
-           x = "Mass difference [Da]",
-           y = "-log(P value)")
-  }
-
-  if(adjust_axes){
-
-    x_max <- ceiling(max(abs(vol_data[["D_diff"]])))
-    y_max <- ceiling(max(vol_data[["log_p_value"]])) + 2
-
-    volcano_plot <- volcano_plot +
-      coord_cartesian(xlim = c(-x_max, x_max), ylim = c(0, y_max), expand = FALSE)
-
-  }
-
-  if(show_confidence_limits){
-
-    y_threshold <- -log(1 - confidence_level)
-
-    t_value <- qt(c((1 - confidence_level)/2, 1-(1 - confidence_level)/2), df = 2)[2]
-    x_threshold <- t_value * mean(vol_data[["Uncertainty"]], na.rm = TRUE)/sqrt(length(vol_data))
-
-    volcano_plot <- volcano_plot +
-      geom_segment(aes(x = -x_threshold, xend = -x_threshold, y = y_threshold, yend = Inf), linetype = "dashed", color = "red") +
-      geom_segment(aes(x = x_threshold, xend = x_threshold, y = y_threshold, yend = Inf), linetype = "dashed", color = "red") +
-      geom_segment(aes(y = y_threshold, yend = y_threshold, x = -Inf, xend = -x_threshold), linetype = "dashed", color = "red") +
-      geom_segment(aes(y = y_threshold, yend = y_threshold, x = x_threshold, xend = Inf), linetype = "dashed", color = "red")
-
-  }
-
-  return(volcano_plot)
-}

R/kinetics.R

@@ -355,8 +355,8 @@ show_kinetic_data <- function(kin_dat,
 #'                            state = "CD160",
 #'                            start = 1, 
 #'                            end = 15,
-#'                            time_in = 0.001, 
-#'                            time_out = 1440)
+#'                            time_0 = 0.001, 
+#'                            time_100 = 1440)
 #' plot_kinetics(kin_dat = kin1, 
 #'               theoretical = FALSE, 
 #'               fractional = TRUE)
@@ -433,13 +433,15 @@ plot_kinetics <- function(kin_dat,
 
   kin_plot <- plot_dat %>% 
     ggplot(aes(x = time_chosen, y = value, group = prop)) +
-    geom_point(aes(color = prop)) + 
+    geom_point(aes(color = prop), size = 2) + 
     theme(legend.position = "bottom",
           legend.title = element_blank()) +
     labs(x = "Time points [min]", 
          y = y_label,
          title = title)
 
+  if(log_x){ err_width = 0.1 } else { err_width = 5 }
+
   if(uncertainty_type == "ribbon"){
 
     kin_plot <- kin_plot +
@@ -449,12 +451,14 @@ plot_kinetics <- function(kin_dat,
   } else if (uncertainty_type == "bars") {
 
     kin_plot <- kin_plot +
-      geom_errorbar(aes(x = time_chosen, ymin = value - err_value, ymax = value + err_value, color = prop))
+      geom_errorbar(aes(x = time_chosen, ymin = value - err_value, ymax = value + err_value, color = prop),
+                    width = err_width)
 
   } else if (uncertainty_type == "bars + line"){
 
     kin_plot <- kin_plot +
-      geom_errorbar(aes(x = time_chosen, ymin = value - err_value, ymax = value + err_value, color = prop)) + 
+      geom_errorbar(aes(x = time_chosen, ymin = value - err_value, ymax = value + err_value, color = prop),
+                    width = err_width) + 
       geom_line(aes(color = prop))
 
   }
@@ -467,4 +471,131 @@ plot_kinetics <- function(kin_dat,
   }
 
   kin_plot
+}
+
+#' Plot differential uptake curve
+#' 
+#' @description 
+#' 
+#' @param diff_uptake_dat produced by \code{\link{create_diff_uptake_dataset}} function
+#' @param sequence 
+#' @param theoretical \code{logical}, determines if plot shows theoretical values.
+#' @param fractional \code{logical}, determines if plot shows fractional values.
+#' @param uncertainty_type type of presenting uncertainty, possible values:
+#' "ribbon", "bars" or "bars + line".
+#' @param log_x \code{logical}, determines if x axis shows logarithmic values.
+#' 
+#' @details Currently there is no possibility to plot multiple peptides on the plot.
+#' 
+#' @return a \code{\link{ggplot2}} object.
+#' 
+#' @seealso 
+#' \code{\link{read_hdx}}
+#' \code{\link{create_diff_uptake_dataset}}
+#' 
+#' @examples 
+#' dat <- read_hdx(system.file(package = "HaDeX", "HaDeX/data/KD_180110_CD160_HVEM.csv"))
+#' diff_uptake_dat <- create_diff_uptake_dataset(dat)
+#' plot_differential_uptake_curve(diff_uptake_dat, sequence = "LCKDRSGDCSPETSLKQL")
+#' 
+#' @export plot_differential_uptake_curve
+
+plot_differential_uptake_curve <- function(diff_uptake_dat,
+                                           sequence = NULL,
+                                           theoretical = FALSE,
+                                           fractional = FALSE,
+                                           uncertainty_type = "ribbon",
+                                           log_x = TRUE){
+
+  if(is.null(sequence)){ sequence <- diff_uptake_dat[["Sequence"]][1] }
+
+  diff_uptake_dat <- diff_uptake_dat %>%
+    filter(Sequence == sequence)
+
+  if (theoretical){
+
+    title <- "Theoretical differential uptake curve"
+
+    if (fractional){
+
+      value <- "diff_theo_frac_deut_uptake"
+      err_value <- "err_diff_theo_frac_deut_uptake"
+      y_label <- "Fractional differential uptake [%]"
+
+    } else {
+
+      value <- "diff_theo_deut_uptake"
+      err_value <- "err_diff_theo_deut_uptake"
+      y_label <- "Differential uptake [Da]"
+
+    }
+
+  } else {
+
+    title <- "Differential uptake curve"
+
+    if (fractional){
+
+      value <- "diff_frac_deut_uptake"
+      err_value <- "err_diff_frac_deut_uptake"
+      y_label <- "Fractional differential uptake [%]"
+
+    } else {
+
+      value <- "diff_deut_uptake"
+      err_value <- "err_diff_deut_uptake"
+      y_label <- "Differential uptake [Da]"
+
+    }
+
+  }
+
+  plot_dat <- data.frame(Sequence = diff_uptake_dat[["Sequence"]],
+                         Start = diff_uptake_dat[["Start"]],
+                         End = diff_uptake_dat[["End"]],
+                         time_chosen = diff_uptake_dat[["Exposure"]],
+                         value = diff_uptake_dat[[value]],
+                         err_value = diff_uptake_dat[[err_value]])
+
+  diff_kin_plot <- plot_dat %>% 
+    ggplot(aes(x = time_chosen, y = value, group = Sequence)) +
+    geom_point(aes(color = Sequence), size = 2) + 
+    theme(legend.position = "bottom",
+          legend.title = element_blank()) +
+    labs(x = "Time points [min]", 
+         y = y_label,
+         title = title)
+
+  if(log_x){ err_width = 0.1 } else { err_width = 5 }
+
+  if(uncertainty_type == "ribbon"){
+
+    diff_kin_plot <- diff_kin_plot +
+      geom_ribbon(aes(ymin = value - err_value, ymax = value + err_value, fill = Sequence), alpha = 0.15) +
+      geom_line(aes(color = Sequence)) 
+
+  } else if (uncertainty_type == "bars") {
+
+    diff_kin_plot <- diff_kin_plot +
+      geom_errorbar(aes(x = time_chosen, ymin = value - err_value, ymax = value + err_value, color = Sequence),
+                    width = err_width)
+
+  } else if (uncertainty_type == "bars + line"){
+
+    diff_kin_plot <- diff_kin_plot +
+      geom_errorbar(aes(x = time_chosen, ymin = value - err_value, ymax = value + err_value, color = Sequence),
+                    width = err_width) + 
+      geom_line(aes(color = Sequence))
+
+  }
+
+  if(log_x){
+
+    diff_kin_plot <- diff_kin_plot + 
+      scale_x_log10()
+
+  }
+
+  diff_kin_plot
+
 }