Mercurial > repos > ecology > pampa_communitymetrics

comparison FunctExeCalcCommIndexesGalaxy.r @ 5:cc0b32aa574f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/ColineRoyaux/PAMPA-Galaxy commit 04381ca7162ec3ec68419e308194b91d11cacb04"
author ecology
date Mon, 16 Nov 2020 11:01:21 +0000
parents fd0c1dd5db33
children a6729da1c623
comparison
equal deleted inserted replaced
4:b79530a322a9 5:cc0b32aa574f
1 #Rscript 1 #Rscript
2 2
3 ##################################################################################################################### 3 #####################################################################################################################
4 ##################################################################################################################### 4 #####################################################################################################################
5 ################################# Calculate community indexes from observation data ################################# 5 ################################# Calculate community indexes from observation data #################################
6 ##################################################################################################################### 6 #####################################################################################################################
7 ##################################################################################################################### 7 #####################################################################################################################
8 8
9 ###################### Packages R 9 ###################### Packages R
10 10
11 suppressMessages(library(tidyr)) 11 suppressMessages(library(tidyr))
12 12
13 ###################### Load arguments and declaring variables 13 ###################### Load arguments and declaring variables
14 14
15 args = commandArgs(trailingOnly=TRUE) 15 args <- commandArgs(trailingOnly = TRUE)
16 #options(encoding = "UTF-8")
17 16
18 if (length(args) < 4) { 17 if (length(args) < 4) {
19 stop("At least one argument must be supplied, an input dataset file (.tabular).", call.=FALSE) #si pas d'arguments -> affiche erreur et quitte / if no args -> error and exit1 18 stop("At least one argument must be supplied, an input dataset file (.tabular).", call. = FALSE) # if no args -> error and exit1
20 19
21 } else { 20 } else {
22 Importdata<-args[1] ###### Nom du fichier importé avec son extension / file name imported with the file type ".filetype" 21 import_data <- args[1] ###### Nom du fichier importé avec son extension / file name imported with the file type ".filetype"
23 index <- args[2] ###### List of selected metrics to calculate 22 index <- args[2] ###### List of selected metrics to calculate
24 source(args[3]) ###### Import functions 23 source(args[3]) ###### Import functions
25
26 } 24 }
27 #### Data must be a dataframe with at least 3 variables : unitobs representing location and year ("observation.unit"), species code ("species.code") and abundance ("number") 25 #### d_ata must be a dataframe with at least 3 variables : unitobs representing location and year ("observation.unit"), species code ("species.code") and abundance ("number")
28 26
29 27
30 #Import des données / Import data 28 #Import des données / Import data
31 obs<- read.table(Importdata,sep="\t",dec=".",header=TRUE,encoding="UTF-8") # 29 obs <- read.table(import_data, sep = "\t", dec = ".", header = TRUE, encoding = "UTF-8") #
32 obs[obs == -999] <- NA 30 obs[obs == -999] <- NA
33 factors <- fact.det.f(Obs=obs) 31 factors <- fact_det_f(obs = obs)
34 ObsType <- def.typeobs.f(Obs=obs) 32 obs_type <- def_typeobs_f(obs = obs)
35 obs <- create.unitobs(data=obs) 33 obs <- create_unitobs(data = obs)
36 34
37 vars_data<-c("observation.unit","species.code","number") 35 vars_data <- c("observation.unit", "species.code", "number")
38 err_msg_data<-"The input dataset doesn't have the right format. It need to have at least the following 3 variables :\n- observation.unit (or point and year)\n- species.code\n- number\n" 36 err_msg_data <- "The input dataset doesn't have the right format. It need to have at least the following 3 variables :\n- observation.unit (or location and year)\n- species.code\n- number\n"
39 check_file(obs,err_msg_data,vars_data,3) 37 check_file(obs, err_msg_data, vars_data, 3)
40 38
41 39
42 40
43 #################################################################################################### 41 ####################################################################################################
44 ################## create community metrics table ## Function : calcBiodiv.f ####################### 42 ################# create community metrics table ## Function : calc_biodiv_f #######################
45 #################################################################################################### 43 ####################################################################################################
46 44
47 ######################################################################################################################## 45 ########################################################################################################################
48 calcBiodiv.f <- function(Data, MPA, unitobs="observation.unit", code.especes="species.code", nombres="number", 46 calc_biodiv_f <- function(d_ata, unitobs = "observation.unit", code_species = "species.code", nombres = "number",
49 indices=index) 47 indices = index) {
50 { 48 ## Purpose: compute biodiversity indexes
51 ## Purpose: calcul des indices de biodiversité
52 ## ---------------------------------------------------------------------- 49 ## ----------------------------------------------------------------------
53 ## Arguments: Data : les données à partir desquelles calculer les 50 ## Arguments: d_ata : input observation file
54 ## indices. Doivent comporter au minimum (colones) : 51 ## unitobs : name of column observation unit
55 ## * unités d'observations/sites 52 ## code_species : name of species column
56 ## * espèces présentes 53 ## nombres : name of abundance column
57 ## * nombre d'individus /espèce/unitobs. 54 ## indices : list of indexes to compute
58 ## refesp : le référentiel espèces.
59 ## MPA : l'AMP (chaîne de charactères).
60 ## unitobs : nom de la colone d'unités d'observation.
61 ## code.especes : nom de la colone d'espèces.
62 ## nombres : nom de la colone de nombres.
63 ## indices : liste des indices à calculer
64 ## (vecteur de caractères)
65 ## ---------------------------------------------------------------------- 55 ## ----------------------------------------------------------------------
66 ## Author: Yves Reecht, Date: 29 oct. 2010, 08:58 56 ## Author: Yves Reecht, Date: 29 oct. 2010, 08:58 modified by Coline ROYAUX in june 2020
67 57
68 ## Supression de tout ce qui n'a pas d'espèce précisee (peut être du non biotique ou identification >= genre) : 58 ## Supress lines that doesn't represent a species :
69 59
70 notspline <- grep("(sp\\.)$|([1-9])$|^(Absencemacrofaune)$|^(NoID)$|^(Acrobranc)$|^(Acrodigit)$|^(Acroencr)$|^(Acrosubm)$|^(Acrotabu)$|^(Adredure)$|^(Adremoll)$|^(Algaturf)$|^(Balimona)$|^(Corablan)$|^(CoradurV)$|^(Coraenal)$|^(Coramor1)$|^(Coramor2)$|^(Coramou)$|^( Dallcora)$|^(Debrcora)$|^(Debris)$|^(Hare)$|^(HexaChar)$|^(MuraCong)$|^(Nacrbran)$|^(Nacrcham)$|^(Nacrencr)$|^(Nacrfoli)$|^(Nacrmass)$|^(Nacrsubm)$|^(Recrcora)$|^(Roche)$|^(Sable)$|^(Vase)$",Data[, code.especes], value=FALSE) 60 notspline <- grep("(sp\\.)$|([1-9])$|^(Absencemacrofaune)$|^(NoID)$|^(Acrobranc)$|^(Acrodigit)$|^(Acroencr)$|^(Acrosubm)$|^(Acrotabu)$|^(Adredure)$|^(Adremoll)$|^(Algaturf)$|^(Balimona)$|^(Corablan)$|^(CoradurV)$|^(Coraenal)$|^(Coramor1)$|^(Coramor2)$|^(Coramou)$|^( Dallcora)$|^(Debrcora)$|^(Debris)$|^(Hare)$|^(HexaChar)$|^(MuraCong)$|^(Nacrbran)$|^(Nacrcham)$|^(Nacrencr)$|^(Nacrfoli)$|^(Nacrmass)$|^(Nacrsubm)$|^(Recrcora)$|^(Roche)$|^(Sable)$|^(Vase)$", d_ata[, code_species], value = FALSE)
71 if (length(notspline) != 0) 61 if (length(notspline) != 0) {
72 { 62 d_ata <- d_ata[-notspline, ]
73 Data <- Data[-notspline, ] 63 }
74 }else{}
75 64
76 ## Suppression des niveaux de facteur inutilisés : 65 ## Suppress unused factor levels :
77 Data <- dropLevels.f(df=Data) 66 d_ata <- .GlobalEnv$drop_levels_f(df = d_ata)
78 67
79 68
80 ## Si les données ne sont pas encore agrégées /espèce/unitobs on le fait ici : 69 ## aggregation of data if not already done :
81 if (nrow(Data) > nrow(expand.grid(unique(Data[ , unitobs]), unique(Data[ , code.especes])))) 70 if (nrow(d_ata) > nrow(expand.grid(unique(d_ata[, unitobs]), unique(d_ata[, code_species])))) {
82 { 71 d_ata <- agregations_generic_f(d_ata = d_ata, metrics = nombres,
83 Data <- agregations.generic.f(Data=Data, metrics=nombres, 72 factors = c(unitobs, code_species),
84 factors=c(unitobs, code.especes), 73 list_fact = NULL)
85 listFact=NULL) 74 }
86 }else{}
87 75
88 df.biodiv <- as.data.frame(as.table(tapply(Data[ , nombres], 76 df_biodiv <- as.data.frame(as.table(tapply(d_ata[, nombres],
89 Data[ , unitobs], 77 d_ata[, unitobs],
90 sum, na.rm=TRUE))) 78 sum, na.rm = TRUE)))
91 79
92 colnames(df.biodiv) <- c(unitobs, nombres) 80 colnames(df_biodiv) <- c(unitobs, nombres)
93 81
94 ## ################################################## 82 ## ##################################################
95 ## Richesse spécifique : 83 ## species richness :
96 Data$pres.abs <- presAbs.f(nombres=Data[ , nombres], logical = FALSE) 84 d_ata$presence_absence <- .GlobalEnv$pres_abs_f(nombres = d_ata[, nombres], logical = FALSE)
97 85
98 df.biodiv$species.richness <- as.vector(tapply(Data$pres.abs, 86 df_biodiv$species_richness <- as.vector(tapply(d_ata$presence_absence,
99 Data[ , unitobs], sum, na.rm=TRUE), 87 d_ata[, unitobs], sum, na.rm = TRUE),
100 "integer") 88 "integer")
101 ## ... as.vector to avoid the class "array". 89 ## ... as.vector to avoid the class "array".
102 90
103 ## ################################################## 91 ## ##################################################
104 ## Indices de Simpson et Shannon et dérivés : 92 ## Simpson, Shannon indexes and derivatives :
105 93
106 matNombres <- tapply(Data[ , nombres], # Matrice de nombres d'individus /espèce/unitobs. 94 mat_nb <- tapply(d_ata[, nombres], # Matrix of individual count /species/unitobs.
107 list(Data[ , unitobs], Data[ , code.especes]), 95 list(d_ata[, unitobs], d_ata[, code_species]),
108 sum, na.rm=TRUE) 96 sum, na.rm = TRUE)
109 97
110 matNombres[is.na(matNombres)] <- 0 # Vrais zéros 98 mat_nb[is.na(mat_nb)] <- 0 # Vrais zéros
111 99
112 ## Proportion d'individus de chaque espèce dans l'unitobs : 100 ## each species individual proportion in the dataset :
113 propIndiv <- sweep(matNombres, 1, 101 prop_indiv <- sweep(mat_nb, 1,
114 apply(matNombres, 1, sum, na.rm = TRUE), # Nombre d'individus / unitobs ; équiv df.biodiv$nombre. 102 apply(mat_nb, 1, sum, na.rm = TRUE), # individual count / unitobs ; equiv df_biodiv$nombre.
115 FUN="/") 103 FUN = "/")
116 104
117 ## Indices de Simpson : 105 ## Simpson indexes :
118 df.biodiv$simpson <- apply(propIndiv^2, 1, sum, na.rm=TRUE) 106 df_biodiv$simpson <- apply(prop_indiv^2, 1, sum, na.rm = TRUE)
119 107
120 if (any(is.element(c("all", "simpson.l"), indices))) 108 if (any(is.element(c("all", "simpson.l"), indices))) {
121 { 109 df_biodiv$simpson_l <- 1 - df_biodiv$simpson
122 df.biodiv$simpson.l <- 1 - df.biodiv$simpson
123 } 110 }
124 111
125 ## calcul de l'indice de Shannon : 112 ## Shannon index :
126 df.biodiv$shannon <- -1 * apply(propIndiv * log(propIndiv), 1, sum, na.rm=TRUE) 113 df_biodiv$shannon <- -1 * apply(prop_indiv * log(prop_indiv), 1, sum, na.rm = TRUE)
127 114
128 ## calcul de l'indice de Pielou : 115 ## Pielou index :
129 if (any(is.element(c("all", "pielou"), indices))) 116 if (any(is.element(c("all", "pielou"), indices))) {
130 { 117 df_biodiv$pielou <- df_biodiv$shannon / log(df_biodiv$species_richness)
131 df.biodiv$pielou <- df.biodiv$shannon / log(df.biodiv$species.richness)
132 } 118 }
133 119
134 ## calcul de l'indice de Hill : 120 ## Hill index :
135 if (any(is.element(c("all", "hill"), indices))) 121 if (any(is.element(c("all", "hill"), indices))) {
136 { 122 df_biodiv$hill <- (1 - df_biodiv$simpson) / exp(df_biodiv$shannon) # equiv df_biodiv$l.simpson / exp(df_biodiv$shannon)
137 df.biodiv$hill <- (1 - df.biodiv$simpson) / exp(df.biodiv$shannon)
138 # équiv df.biodiv$l.simpson / exp(df.biodiv$shannon)
139 } 123 }
140 124
141 125
142 return(df.biodiv) 126 return(df_biodiv)
143 } 127 }
144 128
145 ################# Analysis 129 ################# Analysis
146 130
147 res <- calc.numbers.f(obs, ObsType=ObsType , factors=factors, nbName="number") 131 res <- calc_numbers_f(obs, obs_type = obs_type, factors = factors, nb_name = "number")
148 132
149 tableCommunityIndexes <- calcBiodiv.f(res, MPA, unitobs="observation.unit", code.especes="species.code", nombres="number", 133 table_comm_indexes <- calc_biodiv_f(res, unitobs = "observation.unit", code_species = "species.code", nombres = "number",
150 indices=index) 134 indices = index)
151 tableCommunityIndexes <- create.year.point(tableCommunityIndexes) 135 table_comm_indexes <- create_year_location(table_comm_indexes)
152 #Save dataframe in a tabular format 136 #Save dataframe in a tabular format
153 137
154 filenameComm <- "TabCommunityIndexes.tabular" 138 filename_comm <- "TabCommunityIndexes.tabular"
155 write.table(tableCommunityIndexes, filenameComm, row.names=FALSE, sep="\t", dec=".",fileEncoding="UTF-8") 139 write.table(table_comm_indexes, filename_comm, row.names = FALSE, sep = "\t", dec = ".", fileEncoding = "UTF-8")
156 cat(paste("\nWrite table with Community indexes. \n--> \"",filenameComm,"\"\n",sep="")) 140 cat(paste("\nWrite table with Community indexes. \n--> \"", filename_comm, "\"\n", sep = ""))
157