--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/graph_lcbd.r	Fri Aug 13 18:16:26 2021 +0000
@@ -0,0 +1,225 @@
+#Rscript
+
+#########################
+##    Beta diversity   ##
+#########################
+
+#####Packages : ggplot2
+#               vegan
+#               adespatial
+#               dplyr
+#               tibble
+#               tdyr
+
+#####Load arguments
+
+args <- commandArgs(trailingOnly = TRUE)
+
+if (length(args) < 2) {
+    stop("This tool needs at least 2 arguments")
+}else{
+    table <- args[1]
+    hr <- args[2]
+    abund <- as.numeric(args[3])
+    loc <- as.numeric(args[4])
+    spe <- as.numeric(args[5])
+    date <- as.numeric(args[6])
+    map <- as.logical(args[7])
+    sepa <- as.logical(args[8])
+    not <- as.logical(args[9])
+    lat <- as.numeric(args[10])
+    long <- as.numeric(args[11])
+    var <- as.numeric(args[12])
+    source(args[13])
+}
+
+if (hr == "false") {
+  hr <- FALSE
+}else{
+  hr <- TRUE
+}
+
+#####Import data
+data <- read.table(table, sep = "\t", dec = ".", header = hr, fill = TRUE, encoding = "UTF-8")
+colabund <- colnames(data)[abund]
+colloc <- colnames(data)[loc]
+if (map) {
+  collat <- colnames(data)[lat]
+  collong <- colnames(data)[long]
+}
+colspe <- colnames(data)[spe]
+coldate <- colnames(data)[date]
+data[, coldate] <- as.factor(data[, coldate])
+
+data <- data[grep("^$", data[, spe], invert = TRUE), ]
+
+if (sepa) {
+colvar <- colnames(data)[var]
+}
+
+# Data for species
+data_num <- make_table_analyse(data, colabund, colspe, colloc, coldate)
+nb_spe <- length(unique(data[, spe]))
+nb_col <- ncol(data_num) - nb_spe + 1
+
+#Data with coordinates and environmental
+if (map) {
+  data_xy <- data_num[, c(collat, collong)]
+  colnames(data_xy) <- c("latitude", "longitude")
+  # Data for environment
+  data_env <- data_num[, c(colloc, collat, collong)]
+  colnames(data_env) <- c("site", "latitude", "longitude")
+}
+
+# Data with only species and their abundance
+data_spe <- data_num[, nb_col:ncol(data_num)]
+rownames(data_spe) <- paste0(data_num[, colloc], " - ", data_num[, coldate])
+
+#####Your analysis
+
+# Computation beta.div {adespatial}
+# Beta.div on Hellinger-transformed species data
+data_beta <- adespatial::beta.div(data_spe, method = "hellinger", nperm = 9999)
+
+save(data_beta, file = "beta_diversity.Rdata")
+cat("##############################################################################",
+    "\n########################### Beta Diversity Summary ###########################",
+    "\n##############################################################################",
+    "\n\n### All data ###",
+    "\nBeta diversity: ", data_beta$beta[[2]],
+    "\nSum of Squares: ", data_beta$beta[[1]],
+    "\n\n### Vector of Local Contributions to Beta Diversity (LCBD) for the sites each date ###",
+    "\n", capture.output(data_beta$LCBD),
+    "\n\n### Vector of P-values associated with the LCBD indices ###",
+    "\n", capture.output(data_beta$p.LCBD),
+    "\n\n### Vector of Corrected P-values for the LCBD indices, Holm correction ###",
+    "\n", capture.output(data_beta$p.adj),
+    "\n\n### Vector of Species contributions to beta diversity (SCBD) ###",
+    "\n", capture.output(data_beta$SCBD), file = "LCBD.txt", fill = 1, append = TRUE)
+
+# Which species have a SCBD larger than the mean SCBD?
+scbd <- capture.output(data_beta$SCBD[data_beta$SCBD >= mean(data_beta$SCBD)])
+write(scbd, "SCBD.txt")
+
+##1st fonction
+beta_div_ext <- function(data_beta, data_xy, data_env) {
+   data_beta_ext <- data.frame(data_xy, data_env, LCBD = data_beta$LCBD * 100, p.LCBD = data_beta$p.LCBD, signif = data_beta$p.LCBD < 0.05)
+
+  graph_beta_ext <- ggplot2::ggplot(data = data_beta_ext, ggplot2::aes(x = latitude, y = longitude, size = LCBD, col = signif)) +
+  ggplot2::geom_point() +
+  ggplot2::scale_colour_manual(values = c("#57bce0", "#ce0b0b"), labels = c("Non significant", "Significant"), name = "Significance at 0.05") +
+  ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude")
+
+  ggplot2::ggsave("Beta_diversity_through_space.png", graph_beta_ext)
+}
+
+## Boyé et al. 2017 JSR Fig R
+####################################################
+
+####LCBD####
+lcbd_site <- adespatial::beta.div(data_spe, "hellinger", nperm = 999)
+
+compute_lcbd <- function(data_beta, data_spe, data_num) {
+
+#############
+  mat_lcbd_site <- data.frame(data_spe, LCBD = data_beta$LCBD * 100, p.LCBD = data_beta$p.LCBD, signif = data_beta$p.LCBD < 0.05, site = data_num[, colloc], date = data_num[, coldate])
+
+## Map spatio-temp
+##################
+  p1 <- ggplot2::qplot(date, site, size = LCBD, col = signif, data = mat_lcbd_site)
+  p1 <- p1 + ggplot2::scale_colour_manual(values = c("#57bce0", "#ce0b0b"), labels = c("Non significant", "Significant"), name = "Significance at 0.05")
+  p1 <- p1 + ggplot2::theme_bw() + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) + ggplot2::xlab("Date") + ggplot2::ylab("Site")
+
+  ggplot2::ggsave("LCBD_sites_time.png", p1)
+
+
+## Par années
+#############
+  mean_time <- tapply(mat_lcbd_site$LCBD, mat_lcbd_site$date, mean)
+  sd_time <- tapply(mat_lcbd_site$LCBD, mat_lcbd_site$date, sd)
+  date <- unique(mat_lcbd_site$date)
+
+  data <- data.frame(date, mean_time, sd_time)
+
+  time <- ggplot2::ggplot() + ggplot2::geom_pointrange(ggplot2::aes(x = date, y = mean_time, ymin = mean_time - sd_time, ymax = mean_time + sd_time), data = data)
+  time <- time + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90), axis.line.y = ggplot2::element_line(size = 0.5)) + ggplot2::ylab("mean LCBD")
+
+  ggplot2::ggsave("Mean_LCBD_through_time.png", time)
+}
+
+## Choose another graph
+#######################
+compute_lcbd2 <- function(data_beta, data_spe, data_num) {
+
+#############
+  mat_lcbd_site <- data.frame(data_spe, LCBD = data_beta$LCBD * 100, p.LCBD = data_beta$p.LCBD, signif = data_beta$p.LCBD < 0.05, site = data_num[, colloc], date = data_num[, coldate], variable = data_num[, colvar])
+
+  p1 <- ggplot2::qplot(date, variable, size = LCBD, col = signif, data = mat_lcbd_site)
+  p1 <- p1 + ggplot2::scale_colour_manual(values = c("#57bce0", "#ce0b0b"), labels = c("Non significant", "Significant"), name = "Significance at 0.05")
+  p1 <- p1 + ggplot2::theme_bw() + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) + ggplot2::xlab("Date") + ggplot2::ylab(colvar)
+
+  ggplot2::ggsave(paste0("LCBD_per_", colvar, "_through_time.png"), p1)
+}
+
+####SCBD###
+# Function to compute SCBD
+library(dplyr)
+make_scbd_uvc <- function(data_spe, z, data_beta) {
+  # Computation using beta.div {adespatial} on
+  # Hellinger-transformed species data
+
+  # Which species have a SCBD larger than the mean SCBD?
+  spe_scbd <-  data_beta$SCBD[data_beta$SCBD >= mean(data_beta$SCBD)] %>%
+    as.data.frame() %>%
+    tibble::rownames_to_column(var = "Taxon") %>%
+    dplyr::mutate("Methode" = z)
+
+  return(spe_scbd)
+}
+
+# Function to make a radar plot
+
+coord_radar <- function(theta = "x", start = 0, direction = 1) {
+  theta <- match.arg(theta, c("x", "y"))
+  r <- if (theta == "x") "y" else "x"
+  ggplot2::ggproto("CordRadar", ggplot2::coord_polar(theta = theta, start = start,
+          direction = sign(direction)),
+          is_linear = function(coord) TRUE)
+}
+
+# Make the radar plot
+radar_plot <- function(scbd_uvc_tc) {
+  uvc_rd_plot_data <- scbd_uvc_tc %>%
+    rename(scbd_score = ".")
+
+  rad_uvc <- ggplot2::ggplot(uvc_rd_plot_data, ggplot2::aes(x = Taxon, y = scbd_score, group = Methode)) +
+    ggplot2::geom_line() +
+    ggplot2::geom_point(size = 3) +
+    coord_radar() +
+    ggplot2::theme_bw() +
+    ggplot2::theme(axis.text.x = ggplot2::element_text(size = 10),
+        legend.position = "bottom")
+
+  ggplot2::ggsave("SCBD_Species_Radar_plot.png", rad_uvc)
+}
+
+## LCBD
+
+if (map) {
+  #Beta diversity
+  beta_div_ext(data_beta, data_xy, data_env)
+}
+
+#Lcbd per places and time
+compute_lcbd(data_beta, data_spe, data_num)
+
+#Lcbd of your choice
+if (sepa) {
+  compute_lcbd2(data_beta, data_spe, data_num)
+}
+
+##SCBD
+
+scbd_uvc_tc <- make_scbd_uvc(data_spe, z = "TC", data_beta)
+
+radar_plot(scbd_uvc_tc)