Mercurial > repos > ecology > cb_dissim

changeset 1:ffaa9a71b2d2 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/champ_blocs commit 0d86db7d42b608c386a54500064f5f9c9d7019a4
author ecology
date Wed, 04 Jan 2023 13:21:53 +0000
parents 7e6cc3da1189
children 8dc082da41c1
files cb_dissimilarity.xml div_index.r macro.xml test-data/Clean_Rdata.rdata
diffstat 4 files changed, 197 insertions(+), 9 deletions(-) [+]
line wrap: on
line diff
--- a/cb_dissimilarity.xml	Fri Dec 02 16:13:18 2022 +0000
+++ b/cb_dissimilarity.xml	Wed Jan 04 13:21:53 2023 +0000
@@ -1,14 +1,14 @@
-<tool id="cb_dissim" name="Dissimilarity" version="@VERSION@" profile = "20.01">
+<tool id="cb_dissim" name="QECB-Dissimilarity" version="@VERSION@" profile = "20.01">
     <description>indicator</description>
     <macros>
         <import>macro.xml</import>
     </macros>
     <expand macro="cb_requirements">
-        <requirement type="package" version="2.16">r-rmarkdown</requirement>
+        <requirement type="package" version="2.19">r-rmarkdown</requirement>
         <requirement type="package" version="0.2.4">r-officedown</requirement>
-        <requirement type="package" version="2.6_2">r-vegan</requirement>
+        <requirement type="package" version="2.6_4">r-vegan</requirement>
         <requirement type="package" version="1.3.2">r-tidyverse</requirement>
-        <requirement type="package" version="3.3.6">r-ggplot2</requirement>
+        <requirement type="package" version="3.4">r-ggplot2</requirement>
     </expand>
     <required_files>
         <include type="literal" path="test-data/champbloc_qecb.csv"/>
@@ -40,9 +40,9 @@
         <param name="choice_date" type="integer" value="2021" label="Until when do you have data (write only the YEAR) ?" help="Example write : 2021 "/>
     </inputs>
     <outputs>
-        <data name="output_qecb" from_work_dir="qecbnato0.RDS" format="Rdata" label="Clean Rdata"/>
-        <data name="output_eco" from_work_dir="bret_egmp_basq_qecb.RDS" format="Rdata" label="Bret Rdata"/>
-        <data name="output_matri" from_work_dir="matri_full.RDS" format="Rdata" label="Matrice full Rdata"/>
+        <data name="output_qecb" from_work_dir="qecbnato0.RDS" format="rdata" label="Clean Rdata"/>
+        <data name="output_eco" from_work_dir="bret_egmp_basq_qecb.RDS" format="rdata" label="Bret Rdata"/>
+        <data name="output_matri" from_work_dir="matri_full.RDS" format="rdata" label="Matrice full Rdata"/>
         <collection type="list" name="plots_bf" label="Bloc fixe et bloc mobile face sup plot">
             <discover_datasets pattern="bm(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
         </collection>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/div_index.r	Wed Jan 04 13:21:53 2023 +0000
@@ -0,0 +1,187 @@
+# author: "Jonathan Richir"
+# date: "01 October 2022"
+
+
+#Rscript
+
+###############################
+##  ##
+###############################
+
+#####Packages : dplyr
+#               tidyr
+#               readr
+#               writexl
+#               stringr
+#               readxl
+#               tibble
+#               lubridate
+#               cowplot
+#               magrittr
+#               rmarkdown
+library(magrittr)
+library(dplyr)
+#####Load arguments
+
+args <- commandArgs(trailingOnly = TRUE)
+
+#####Import data
+
+if (length(args) < 1) {
+    stop("This tool needs at least 1 argument")
+}else {
+    qecnato0 <- args[1]
+
+}
+
+qecnato0 <- readRDS(qecnato0)
+
+
+# first, create vector (4) for qecb and fishing by region (same as above)
+
+bret_egmp_basq_qecb <- c(
+  "X..algues.brunes",
+  "X..algues.rouges",
+  "X..algues.vertes",
+  "X..Cladophora",
+  "X..Lithophyllum",
+  "Nb.Littorina.obtusata",
+  "Nb.Gibbula.cineraria",
+  "Nb.Gibbula.pennanti",
+  "Nb.Gibbula.umbilicalis",
+  "Nb.Phallusia.mamillata",
+  "Nb.Tethya.aurantium",
+  "Nb.Spirobranchus.lamarckii.total",
+  "Nb.spirorbis.total",
+  "X..Eponges",
+  "X..Ascidies.Coloniales",
+  "X..Ascidies.Solitaires",
+  "X..Bryozoaires.Dresses",
+  "X..Balanes.Vivantes"
+  #, "X..Recouvrement.Sediment"
+  #, "X..Roche.Nue"
+  #, "X..Surface.Accolement"
+  )
+
+egmp_basq_qecb <- c("Nb.Crassostrea.gigas", "Nb.Ostrea.edulis", "X..Mytilus.sp.", "X..Hermelles", "X..Hydraires")
+
+bret_egmp_basq_fishing <- c("Nb.Cancer.pagurus..Tourteau.",
+                            "Nb.Necora.puber..Etrille.",
+                            "Nb.Carcinus.maenas..Crabe.vert.",
+                            "Nb.Nucella.lapilus..Pourpre.",
+                            "Nb.Galathea..Galathées.",
+                            "Nb.Lophozozymus.incisus..ancien.Xantho.incisus.",
+                            "Nb.Palaemon.sp..Crevette.bouquet.ou.crevette.rose.",
+                            "Nb.Haliotis.tuberculata..Ormeau.",
+                            "Nb.Littorina.littorea..Bigorneau.",
+                            "Nb.Xantho.pilipes..Xanthe.poilu.",
+                            "Nb.Mimachlamys.varia..Pétoncle.noir.")
+
+egmp_basq_fishing <- c("Nb.Eriphia.verrucosa..Crabe.verruqueux.", "Nb.Octopus.vulgaris..Poulpe.", "Nb.Paracentrotus.lividus..Oursin.", "Nb.Stramonita.haemastoma..Pourpre.bouche.de.sang.")
+
+# here I can choose to either replace spirorbis and/or spirobranchus by their log10 transformation in bret_egmp_basq_qecb vector
+bret_egmp_basq_qecb <- replace(bret_egmp_basq_qecb, bret_egmp_basq_qecb == "Nb.spirorbis.total", "log10.Nb.spirorbis.total")
+
+## Diversity index
+
+
+# adiv contains two main functions for species diversity indices: speciesdiv, which includes widely used indices such as species richness and the Shannon index, and divparam, which includes indices that have a parameter to control the importance given to rare versus abundant species in diversity measurements (Pavoine (2020) - adiv: An r package to analyse biodiversity in ecology).
+
+# NB: just like for dissimilarity distance matrices, no sense to use the "fishing" variable lists, because either they are present for the bloc mobile and not for the bloc fixe (therefore false higher diversity for bloc mobile), either they are repeated between face supérieure and face inférieure of bloc mobile.
+
+# function in a loop
+
+row.names(qecnato0) <- c(paste0(qecnato0$region.site_year_month_day, "_", qecnato0$Quadrat.bis, "_", qecnato0$Type.Bloc, "_", qecnato0$Numéro.Bloc.échantillon, "_", qecnato0$Face))
+
+# later on I can copy-paste above code to recreate variable names vector
+#bret_egmp_basq_qecb
+#egmp_basq_qecb
+#Bret_EGMP.BASQ_fishing
+#EGMP.BASQ_fishing
+
+# remove boulder variables
+bret_egmp_basq_qecb <- bret_egmp_basq_qecb[! bret_egmp_basq_qecb %in% c("X..Recouvrement.Sediment", "X..Roche.Nue", "X..Surface.Accolement")]
+
+qecnato0$period <- as.character(qecnato0$period)
+qecnato0$Face <- as.character(qecnato0$Face)
+
+div_list <- vector("list", length(unique(qecnato0$site_year_month_day)))
+
+for (i in c(1:nrow(qecnato0))) {
+  div_i <- dplyr::filter(qecnato0, site_year_month_day == qecnato0$site_year_month_day[i])
+
+  ifelse(unique(div_i$region) == "Bretagne", var. <- c(bret_egmp_basq_qecb), var. <- c(bret_egmp_basq_qecb, egmp_basq_qecb)) # Qu. : Why can't R's ifelse statements return vectors? => you can circumvent the problem if you assign the result inside the ifelse.
+
+  #8 remove empty row cfr: In speciesdiv(div_i[, var.]) & divparam(div_i[, var.]) : empty communities should be discarded
+  div_i <- dplyr::filter(div_i, rowSums(div_i[, var.]) > 0)
+
+  div_i_speciesdiv <- adiv::speciesdiv(div_i[, var.])
+  adiv_i_df <- data.frame(div_i_speciesdiv)
+  
+  div_i_divparam <- adiv::divparam(div_i[, var.], q = c(0, 0.25, 0.5, 1, 2, 4, 8)) # When q increases, abundant species are overweighted compared to rare species, we thus expect that the evenness in species weights decreases.
+
+
+  par(mfrow = c (1, 1))
+  plot(adiv::divparam(div_i[, var.], q = 0), main = unique(div_i$site_year_month_day))
+  plot(adiv::divparam(div_i[, var.], q = 0:10), legend = FALSE, main = unique(div_i$site_year_month_day))
+  
+  adiv_i_df$x <- div_i_divparam$div$`1`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[1], " (equi. richness)")
+  adiv_i_df$x <- div_i_divparam$div$`2`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[2])
+  adiv_i_df$x <- div_i_divparam$div$`3`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[3])
+  adiv_i_df$x <- div_i_divparam$div$`4`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[4])
+  adiv_i_df$x <- div_i_divparam$div$`5`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[5], " (equi. Simpson)")
+  adiv_i_df$x <- div_i_divparam$div$`6`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[6])
+  adiv_i_df$x <- div_i_divparam$div$`7`
+  colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[7])
+  
+  # plot
+  par(mfrow = c(3, 2))
+  sapply(names(adiv_i_df[, c(1, 8, 2, 3, 12, 4:7, 9:11, 13:ncol(adiv_i_df))]), 
+         function(cname) {
+           png(paste0(cname, "_histo.png"))
+           hist(adiv_i_df[, c(1, 8, 2, 3, 12, 4:7, 9:11, 13:ncol(adiv_i_df))][[cname]], main = "", xlab = cname, breaks = length(unique(adiv_i_df[, c(1, 8, 2, 3, 12, 4:7, 9:11, 13:ncol(adiv_i_df))][[cname]])))
+dev.off()
+         }
+
+)
+  par(mfrow = c(1,1))
+  
+  div_list[[i]] <- adiv_i_df
+  
+  rm(div_i, adiv_i_df, div_i_speciesdiv, div_i_divparam)
+  
+}
+
+# for the error message due to richness NA data => 35 observations in 21 surveys; no reason to remove these data "...remo", was checked in the complete script.
+
+
+div_df <- do.call("rbind", div_list)
+
+
+# There is an issue with region.terri that are merged with no "_" ...
+
+div_df <- tibble::add_column(div_df, rownames. = rownames(div_df), .before = "richness")
+div_df <- tidyr::separate(div_df, rownames., into = c("region.terri.", "site_year_month_day", "Quadrat.bis", "Type.Bloc", "Numéro.Bloc.échantillon", "Face"), sep = "_")
+
+# I therefore add these lines to solve that issue
+
+div_df <- tibble::add_column(div_df, terri. = substring(div_df$region.terri., nchar(div_df$region.terri.)-3), .after = "region.terri.")
+
+div_df$region.terri. <- substring(div_df$region.terri., 1, nchar(div_df$region.terri)-4)
+div_df <- dplyr::rename(div_df, region = region.terri.)
+
+div_df$site_year_month_day <- paste0(div_df$terri., "_", div_df$site_year_month_day)
+div_df <- subset(div_df, select = -c(terri.))
+
+div_df$Type.Bloc <- as.factor(div_df$Type.Bloc)
+div_df$Face <- as.factor(div_df$Face)
+div_df$Numéro.Bloc.échantillon <- as.integer(div_df$Numéro.Bloc.échantillon)
+
+saveRDS(div_df, "div_df.RDS")
+write.table(div_df, "Valeurs_stat.tabular", row.names = FALSE, quote = FALSE, sep = "\t", dec = ".", fileEncoding = "UTF-8")
--- a/macro.xml	Fri Dec 02 16:13:18 2022 +0000
+++ b/macro.xml	Wed Jan 04 13:21:53 2023 +0000
@@ -2,11 +2,12 @@
     <token name="@VERSION@">0.0.0</token>
     <xml name="cb_requirements">
         <requirements>
+            <requirement type="package" version="4.2.2">r-base</requirement>
             <requirement type="package" version="1.0.10">r-dplyr</requirement>
             <requirement type="package" version="1.2.1">r-tidyr</requirement>
-            <requirement type="package" version="1.4.1">r-stringr</requirement>
+            <requirement type="package" version="1.5.0">r-stringr</requirement>
             <requirement type="package" version="3.1.8">r-tibble</requirement>
-            <requirement type="package" version="1.40">r-knitr</requirement>
+            <requirement type="package" version="1.41">r-knitr</requirement>
             <requirement type="package" version="2.0.3">r-magrittr</requirement>
            <yield/>
         </requirements>
Binary file test-data/Clean_Rdata.rdata has changed