Mercurial > repos > ecology > ecology_stat_presence_abs

comparison graph_pres_abs_abund.r @ 0:e5552099d0e3 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/Data_explo_tools commit 2f883743403105d9cac6d267496d985100da3958"
author ecology
date Tue, 27 Jul 2021 16:57:02 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e5552099d0e3
1 #Rscript
2
3 #########################################################
4 ## Presence abscence and abundance in environment ##
5 #########################################################
6
7 #####Packages : ggplot2
8 # vegan
9
10 #####Load arguments
11
12 args <- commandArgs(trailingOnly = TRUE)
13
14 if (length(args) < 5) {
15 stop("This tool needs at least 5 arguments")
16 }else{
17 table <- args[1]
18 hr <- args[2]
19 abundance <- as.logical(args[3])
20 presabs <- as.logical(args[4])
21 rarefaction <- as.logical(args[5])
22 lat <- as.numeric(args[6])
23 long <- as.numeric(args[7])
24 ind <- as.character(args[8])
25 loc <- as.numeric(args[9])
26 num <- as.character(args[10])
27 spe <- as.numeric(args[11])
28 abond <- as.numeric(args[12])
29 }
30
31 if (hr == "false") {
32 hr <- FALSE
33 }else{
34 hr <- TRUE
35 }
36
37 #####Import data
38 data <- read.table(table, sep = "\t", dec = ".", header = hr, fill = TRUE, encoding = "UTF-8")
39
40 if (abundance) {
41 collat <- colnames(data)[lat]
42 collong <- colnames(data)[long]
43 }
44
45 if (presabs) {
46 colloc <- colnames(data)[loc]
47 }
48
49 if (presabs | rarefaction | abundance) {
50 colabond <- colnames(data)[abond]
51 colspe <- colnames(data)[spe]
52 data <- data[grep("^$", data[, colspe], invert = TRUE), ]
53 }
54
55 #####Your analysis
56
57 ####The abundance in the environment####
58
59 ##Representation of the environment##
60
61 ## Mapping
62 graph_map <- function(data, collong, collat, colabond, ind, colspe) {
63 cat("\nCoordinates range\n\nLatitude from ", min(data[, collat], na.rm = TRUE), " to ", max(data[, collat], na.rm = TRUE), "\nLongitude from ", min(data[, collong], na.rm = TRUE), " to ", max(data[, collong], na.rm = TRUE), file = "Data_abund.txt", fill = 1, append = TRUE)
64 if (mult0) {
65 mappy <- ggplot2::ggplot(data, ggplot2::aes_string(x = collong, y = collat, cex = colabond, color = colspe)) +
66 ggplot2::geom_point() + ggplot2::ggtitle(paste("Abundance of", ind, "in the environment")) + ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), legend.text = ggplot2::element_text(size = 8)) + ggplot2::guides(cex = ggplot2::guide_legend(reverse = TRUE))
67
68 }else{
69 mappy <- ggplot2::ggplot(data, ggplot2::aes_string(x = collong, y = collat, cex = colabond, color = colabond)) +
70 ggplot2::geom_point() + ggplot2::ggtitle(paste("Abundance of", ind, "in the environment")) + ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), legend.text = ggplot2::element_text(size = 8)) + ggplot2::guides(cex = ggplot2::guide_legend(reverse = TRUE))
71 }
72 ggplot2::ggsave("mappy.png", mappy, width = 20, height = 9, units = "cm")
73 }
74
75 ####Presence absence abundance####
76
77 ## Histogram
78 graph_hist <- function(data, col1, col2, col3) {
79 cat("\nLocations\n", unique(data[, col1]), file = "Locations.txt", fill = 1, append = TRUE)
80 if (mult1) {
81 for (loc in unique(data[, col1])) {
82 data_cut <- data[data[, col1] == loc, ]
83 data_cut <- data_cut[data_cut[, col3] > 0, ]
84 if (length(unique(data_cut[, col2])) <= 40) {
85 top <- nrow(data_cut)
86 var <- nchar(as.character(round(top * 0.1, digits = 0)))
87 step <- round(top * 0.1, digits = ifelse(var == 1, 1, -var + 1))
88 graph <- ggplot2::ggplot(data_cut) +
89 ggplot2::geom_bar(ggplot2::aes_string(x = col1, fill = col2)) +
90 ggplot2::scale_y_continuous(breaks = seq(from = 0, to = top, by = step)) +
91 ggplot2::theme(plot.title = ggplot2::element_text(color = "black", size = 12, face = "bold")) +
92 ggplot2::ggtitle("Number of presence")
93
94 ggplot2::ggsave(paste("number_in_", loc, ".png"), graph)
95 }else{
96 cat(paste0("\n", loc, " had more than 40 species and plot isn't readable please select a higher taxon level or cut your data"))
97 }
98 }
99 }else{
100 top <- nrow(data)
101 var <- nchar(as.character(round(top * 0.1, digits = 0)))
102 step <- round(top * 0.1, digits = ifelse(var == 1, 1, -var + 1))
103 graph <- ggplot2::ggplot(data) +
104 ggplot2::geom_bar(ggplot2::aes_string(x = col1, fill = col2)) +
105 ggplot2::scale_y_continuous(breaks = seq(from = 0, to = top, by = step)) +
106 ggplot2::theme(plot.title = ggplot2::element_text(color = "black", size = 12, face = "bold")) +
107 ggplot2::ggtitle("Number of individuals")
108
109 ggplot2::ggsave("number.png", graph)
110 }
111 }
112
113 rare <- function(data, spe, abond, rare, num) {
114 # Put the data in form
115 new_data <- as.data.frame(data[, spe])
116 colnames(new_data) <- c("Species")
117 new_data$total <- data[, abond]
118
119 new_data$rarefaction <- as.integer(rare)
120
121 end <- length(unique(new_data$Species))
122 out <- vegan::rarecurve(new_data[, 2:3], step = 10, sample = rarefy_sample, label = FALSE)
123 names(out) <- paste(unique(new_data$Species), sep = "")
124
125 # Coerce data into "long" form.
126 protox <- mapply(FUN = function(x, y) {
127 mydf <- as.data.frame(x)
128 colnames(mydf) <- "value"
129 mydf$species <- y
130 mydf$subsample <- attr(x, "Subsample")
131 mydf <- na.omit(mydf)
132 mydf
133 }, x = out, y = as.list(names(out)), SIMPLIFY = FALSE)
134
135 xy <- do.call(rbind, protox)
136 rownames(xy) <- NULL # pretty
137
138 # Plot.
139
140 if (mult2) {
141 for (spe in unique(data[, spe])) {
142 xy_cut <- xy[xy$species == spe, ]
143 top <- max(xy_cut$subsample)
144 var <- nchar(as.character(round(top * 0.1, digits = 0)))
145 step <- round(top * 0.1, digits = ifelse(var == 1, 1, -var + 1))
146 courbe <- ggplot2::ggplot(xy_cut, ggplot2::aes(x = subsample, y = value)) +
147 ggplot2::theme_gray() +
148 ggplot2::geom_line(size = 1) +
149 ggplot2::scale_x_continuous(breaks = seq(from = 0, to = top, by = step)) +
150 ggplot2::xlab("Abundance") + ggplot2::ylab("Value") +
151 ggplot2::ggtitle("rarefaction curve")
152
153 ggplot2::ggsave(paste("rarefaction_of_", spe, ".png"), courbe)
154 }
155 }else{
156 top <- max(xy$subsample)
157 var <- nchar(as.character(round(top * 0.1, digits = 0)))
158 step <- round(top * 0.1, digits = ifelse(var == 1, 1, -var + 1))
159 courbe <- ggplot2::ggplot(xy, ggplot2::aes(x = subsample, y = value, color = species)) +
160 ggplot2::theme_gray() +
161 ggplot2::geom_line(size = 1) +
162 ggplot2::scale_x_continuous(breaks = seq(from = 0, to = top, by = step)) +
163 ggplot2::xlab("Subsample") + ggplot2::ylab("Value") +
164 ggplot2::ggtitle("rarefaction curves")
165
166 ggplot2::ggsave("rarefaction.png", courbe)
167 }
168 }
169
170 if (abundance) {
171 #Maps
172 mult0 <- ifelse(length(unique(data[, colspe])) > 10, FALSE, TRUE)
173 graph_map(data, collong = collong, collat = collat, colabond = colabond, ind = ind, colspe = colspe)
174 }
175
176 if (presabs) {
177 #Presence absence count
178 mult1 <- ifelse(length(unique(data[, colloc])) <= 10, FALSE, TRUE)
179 graph_hist(data, col1 = colloc, col2 = colspe, col3 = colabond)
180 }
181
182 if (rarefaction) {
183 #rarefaction
184
185 #### rarefaction indice ####
186 rarefy_sample <- as.numeric(num)
187
188 ## Calcul de la rarefaction
189 rarefaction <- vegan::rarefy(data[, abond], rarefy_sample)
190
191 write.table(rarefaction, "rare.tabular")
192
193 mult2 <- ifelse(length(unique(data[, colspe])) <= 30, FALSE, TRUE)
194 rare(data, spe = colspe, abond = colabond, rare = rarefaction, num = rarefy_sample)
195 }