diff visualize_pore_diameter_aqp.R @ 2:c574ada16e76 draft

"planemo upload for repository https://github.com/mesocentre-clermont-auvergne/aubi_piaf commit 48a10de1b21f94ab8019d9d0e4a43e0bd9d0c31e-dirty"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/visualize_pore_diameter_aqp.R	Wed May 25 09:42:17 2022 +0000
@@ -0,0 +1,113 @@
+#!/usr/bin/env Rscript
+
+# install and/or load necessary packages
+useful_packages <- c("conflicted", "getopt", "tidyverse", "ggplot2", "ggpubr")
+uninstalled_packages <- setdiff(useful_packages, rownames(installed.packages()))
+invisible(lapply(useful_packages, require, character.only = TRUE, warn.conflicts = TRUE, quietly = TRUE))
+
+spec <- matrix(c(
+  "input_file", "i", 1, "character",
+  "aqp_distribution", "a", 2, "logical",
+  "protomer_distribution", "p", 2, "logical",
+  "all_distribution", "d", 2, "logical",
+  "pdf", "f", 2, "logical"
+), byrow = TRUE, ncol = 4)
+opt <- getopt(spec)
+
+if (is.null(opt$input_file)) {
+  print("A file containing an array must be given as input with the -f argument")
+  quit(status = 1)
+}
+
+if (is.null(opt$aqp_distribution)) {
+  opt$aqp_distribution <- TRUE
+}
+
+if (is.null(opt$protomer_distribution)) {
+  opt$protomer_distribution <- TRUE
+}
+
+if (is.null(opt$all_distribution)) {
+  opt$all_distribution <- TRUE
+}
+
+if (is.null(opt$pdf)) {
+  opt$pdf <- TRUE
+}
+
+tibble_sort <- read.delim(opt$input_file)
+tibble_sort <- as_tibble(tibble_sort)
+colnames(tibble_sort) <- colnames(tibble_sort) %>%
+  as_tibble() %>%
+  mutate(value = str_replace_all(value, "\\.", "_"), value = gsub("_+Angstroms_", "", value), value = str_replace(value, "Time__ps_", "time"), value = str_replace(value, "_", "."), value = str_replace(value, "_", ".")) %>%
+  unlist()
+tibble_sort$time <- as.numeric(gsub("-[0-9]+.[0-9]+", "", tibble_sort$time))
+
+tibble_sort_mean <- tibble_sort[, 1:9]
+tibble_sort_std <-  tibble_sort[, c(1, 10:17)]
+colnames(tibble_sort_mean) <- gsub("_mean", "", colnames(tibble_sort_mean))
+colnames(tibble_sort_std) <- gsub("_std", "", colnames(tibble_sort_std))
+
+tibble_sort_mean_long <- tibble_sort_mean %>%
+  pivot_longer(cols = contains("AQP"), values_to = "distance") %>%
+  separate(name, into = c("aqp", "protomer", "Couple"), sep = "[.]") %>%
+  mutate(protomer = str_replace(protomer, "P", "Protomer "))
+
+tibble_sort_std_long <- tibble_sort_std %>%
+  pivot_longer(cols = contains("AQP"), values_to = "distance") %>%
+  separate(name, into = c("aqp", "protomer", "Couple"), sep = "[.]") %>%
+  mutate(protomer = str_replace(protomer, "P", "Protomer "))
+
+create_ggplot <- function(tibble_mean, tibble_std, group, color, wrap, title) {
+  tibble_mean_distance <- tibble_mean %>%
+    group_by_at(group) %>%
+    summarise(mean_distance = (mean(distance)))
+  tibble_std_distance <- tibble_std %>%
+    group_by_at(group) %>%
+    summarise(std_distance = (mean(distance)))
+  tibble_mean_std_distance <- inner_join(tibble_mean_distance, tibble_std_distance, by = group)
+  g_distribution <- ggplot(tibble_mean_std_distance, aes(x = time / 1000, y = mean_distance)) +
+    geom_line(aes(color = .data[[color]])) +
+    geom_ribbon(aes(ymin = mean_distance - std_distance, ymax = mean_distance
+                    + std_distance, fill = .data[[color]]), alpha = 0.3) +
+    facet_wrap(wrap) + theme(legend.position = "top") + guides(color = "none", fill = "none") +
+    ggtitle(label = title, subtitle = "The envelope represents the standard deviation") +
+    ylab(label = "ArR-ArR distance (Ångströms)") + xlab(label = "Time (ns)")
+  return(g_distribution)
+}
+
+list_ggplot <- list()
+
+if (opt$aqp_distribution == TRUE) {
+  group_aqp <- c("time", "aqp")
+  color_aqp <- "aqp"
+  wrap_aqp <- c("aqp")
+  title_aqp <- "Average distance ArR-ArR by AQP (4 protomers)"
+  g_aqp_distribution <- create_ggplot(tibble_sort_mean_long, tibble_sort_std_long, group_aqp, color_aqp, wrap_aqp, title_aqp)
+  ggsave("Distance_distribution_by_aquaporin.png", g_aqp_distribution, width = 20, height = 18, units = "cm")
+  list_ggplot[[1]] <- g_aqp_distribution
+}
+
+if (opt$protomer_distribution == TRUE) {
+  group_protomer <- c("time", "protomer")
+  color_protomer <- "protomer"
+  wrap_protomer <- c("protomer")
+  title_protomer <- "Average distance ArR-ArR by by protomer (2 aquaporins)"
+  g_protomer_distribution <- create_ggplot(tibble_sort_mean_long, tibble_sort_std_long, group_protomer, color_protomer, wrap_protomer, title_protomer)
+  ggsave("Distance_distribution_by_protomer.png", g_protomer_distribution, width = 20, height = 18, units = "cm")
+  list_ggplot[[2]] <- g_protomer_distribution
+}
+
+if (opt$all_distribution == TRUE) {
+  group_all <- c("time", "aqp", "protomer")
+  color_all <- "protomer"
+  wrap_all <- c("aqp", "protomer")
+  title_all <- "Average distance ArR-ArR by protomer anq AQP"
+  g_all_distribution <- create_ggplot(tibble_sort_mean_long, tibble_sort_std_long, group_all, color_all, wrap_all, title_all)
+  ggsave("Distance_distribution_on_all_protomers.png", g_all_distribution, width = 20, height = 18, units = "cm")
+  list_ggplot[[3]] <- g_all_distribution
+}
+
+if (opt$pdf == TRUE) {
+  ggexport(list_ggplot, filename = "all_graphics_distribution.pdf")
+}