Mercurial > repos > ecology > ecology_presence_abs_abund

comparison graph_lcbd.r @ 1:1cde7ebc234f draft default tip

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 60627aba07951226c8fd6bb3115be4bd118edd4e"
author ecology
date Fri, 13 Aug 2021 18:12:46 +0000
parents
children
comparison
equal deleted inserted replaced
0:9b5a034eadd0 1:1cde7ebc234f
1 #Rscript
2
3 #########################
4 ## Beta diversity ##
5 #########################
6
7 #####Packages : ggplot2
8 # vegan
9 # adespatial
10 # dplyr
11 # tibble
12 # tdyr
13
14 #####Load arguments
15
16 args <- commandArgs(trailingOnly = TRUE)
17
18 if (length(args) < 2) {
19 stop("This tool needs at least 2 arguments")
20 }else{
21 table <- args[1]
22 hr <- args[2]
23 abund <- as.numeric(args[3])
24 loc <- as.numeric(args[4])
25 spe <- as.numeric(args[5])
26 date <- as.numeric(args[6])
27 map <- as.logical(args[7])
28 sepa <- as.logical(args[8])
29 not <- as.logical(args[9])
30 lat <- as.numeric(args[10])
31 long <- as.numeric(args[11])
32 var <- as.numeric(args[12])
33 source(args[13])
34 }
35
36 if (hr == "false") {
37 hr <- FALSE
38 }else{
39 hr <- TRUE
40 }
41
42 #####Import data
43 data <- read.table(table, sep = "\t", dec = ".", header = hr, fill = TRUE, encoding = "UTF-8")
44 colabund <- colnames(data)[abund]
45 colloc <- colnames(data)[loc]
46 if (map) {
47 collat <- colnames(data)[lat]
48 collong <- colnames(data)[long]
49 }
50 colspe <- colnames(data)[spe]
51 coldate <- colnames(data)[date]
52 data[, coldate] <- as.factor(data[, coldate])
53
54 data <- data[grep("^$", data[, spe], invert = TRUE), ]
55
56 if (sepa) {
57 colvar <- colnames(data)[var]
58 }
59
60 # Data for species
61 data_num <- make_table_analyse(data, colabund, colspe, colloc, coldate)
62 nb_spe <- length(unique(data[, spe]))
63 nb_col <- ncol(data_num) - nb_spe + 1
64
65 #Data with coordinates and environmental
66 if (map) {
67 data_xy <- data_num[, c(collat, collong)]
68 colnames(data_xy) <- c("latitude", "longitude")
69 # Data for environment
70 data_env <- data_num[, c(colloc, collat, collong)]
71 colnames(data_env) <- c("site", "latitude", "longitude")
72 }
73
74 # Data with only species and their abundance
75 data_spe <- data_num[, nb_col:ncol(data_num)]
76 rownames(data_spe) <- paste0(data_num[, colloc], " - ", data_num[, coldate])
77
78 #####Your analysis
79
80 # Computation beta.div {adespatial}
81 # Beta.div on Hellinger-transformed species data
82 data_beta <- adespatial::beta.div(data_spe, method = "hellinger", nperm = 9999)
83
84 save(data_beta, file = "beta_diversity.Rdata")
85 cat("##############################################################################",
86 "\n########################### Beta Diversity Summary ###########################",
87 "\n##############################################################################",
88 "\n\n### All data ###",
89 "\nBeta diversity: ", data_beta$beta[[2]],
90 "\nSum of Squares: ", data_beta$beta[[1]],
91 "\n\n### Vector of Local Contributions to Beta Diversity (LCBD) for the sites each date ###",
92 "\n", capture.output(data_beta$LCBD),
93 "\n\n### Vector of P-values associated with the LCBD indices ###",
94 "\n", capture.output(data_beta$p.LCBD),
95 "\n\n### Vector of Corrected P-values for the LCBD indices, Holm correction ###",
96 "\n", capture.output(data_beta$p.adj),
97 "\n\n### Vector of Species contributions to beta diversity (SCBD) ###",
98 "\n", capture.output(data_beta$SCBD), file = "LCBD.txt", fill = 1, append = TRUE)
99
100 # Which species have a SCBD larger than the mean SCBD?
101 scbd <- capture.output(data_beta$SCBD[data_beta$SCBD >= mean(data_beta$SCBD)])
102 write(scbd, "SCBD.txt")
103
104 ##1st fonction
105 beta_div_ext <- function(data_beta, data_xy, data_env) {
106 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)
107
108 graph_beta_ext <- ggplot2::ggplot(data = data_beta_ext, ggplot2::aes(x = latitude, y = longitude, size = LCBD, col = signif)) +
109 ggplot2::geom_point() +
110 ggplot2::scale_colour_manual(values = c("#57bce0", "#ce0b0b"), labels = c("Non significant", "Significant"), name = "Significance at 0.05") +
111 ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude")
112
113 ggplot2::ggsave("Beta_diversity_through_space.png", graph_beta_ext)
114 }
115
116 ## Boyé et al. 2017 JSR Fig R
117 ####################################################
118
119 ####LCBD####
120 lcbd_site <- adespatial::beta.div(data_spe, "hellinger", nperm = 999)
121
122 compute_lcbd <- function(data_beta, data_spe, data_num) {
123
124 #############
125 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])
126
127 ## Map spatio-temp
128 ##################
129 p1 <- ggplot2::qplot(date, site, size = LCBD, col = signif, data = mat_lcbd_site)
130 p1 <- p1 + ggplot2::scale_colour_manual(values = c("#57bce0", "#ce0b0b"), labels = c("Non significant", "Significant"), name = "Significance at 0.05")
131 p1 <- p1 + ggplot2::theme_bw() + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) + ggplot2::xlab("Date") + ggplot2::ylab("Site")
132
133 ggplot2::ggsave("LCBD_sites_time.png", p1)
134
135
136 ## Par années
137 #############
138 mean_time <- tapply(mat_lcbd_site$LCBD, mat_lcbd_site$date, mean)
139 sd_time <- tapply(mat_lcbd_site$LCBD, mat_lcbd_site$date, sd)
140 date <- unique(mat_lcbd_site$date)
141
142 data <- data.frame(date, mean_time, sd_time)
143
144 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)
145 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")
146
147 ggplot2::ggsave("Mean_LCBD_through_time.png", time)
148 }
149
150 ## Choose another graph
151 #######################
152 compute_lcbd2 <- function(data_beta, data_spe, data_num) {
153
154 #############
155 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])
156
157 p1 <- ggplot2::qplot(date, variable, size = LCBD, col = signif, data = mat_lcbd_site)
158 p1 <- p1 + ggplot2::scale_colour_manual(values = c("#57bce0", "#ce0b0b"), labels = c("Non significant", "Significant"), name = "Significance at 0.05")
159 p1 <- p1 + ggplot2::theme_bw() + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) + ggplot2::xlab("Date") + ggplot2::ylab(colvar)
160
161 ggplot2::ggsave(paste0("LCBD_per_", colvar, "_through_time.png"), p1)
162 }
163
164 ####SCBD###
165 # Function to compute SCBD
166 library(dplyr)
167 make_scbd_uvc <- function(data_spe, z, data_beta) {
168 # Computation using beta.div {adespatial} on
169 # Hellinger-transformed species data
170
171 # Which species have a SCBD larger than the mean SCBD?
172 spe_scbd <- data_beta$SCBD[data_beta$SCBD >= mean(data_beta$SCBD)] %>%
173 as.data.frame() %>%
174 tibble::rownames_to_column(var = "Taxon") %>%
175 dplyr::mutate("Methode" = z)
176
177 return(spe_scbd)
178 }
179
180 # Function to make a radar plot
181
182 coord_radar <- function(theta = "x", start = 0, direction = 1) {
183 theta <- match.arg(theta, c("x", "y"))
184 r <- if (theta == "x") "y" else "x"
185 ggplot2::ggproto("CordRadar", ggplot2::coord_polar(theta = theta, start = start,
186 direction = sign(direction)),
187 is_linear = function(coord) TRUE)
188 }
189
190 # Make the radar plot
191 radar_plot <- function(scbd_uvc_tc) {
192 uvc_rd_plot_data <- scbd_uvc_tc %>%
193 rename(scbd_score = ".")
194
195 rad_uvc <- ggplot2::ggplot(uvc_rd_plot_data, ggplot2::aes(x = Taxon, y = scbd_score, group = Methode)) +
196 ggplot2::geom_line() +
197 ggplot2::geom_point(size = 3) +
198 coord_radar() +
199 ggplot2::theme_bw() +
200 ggplot2::theme(axis.text.x = ggplot2::element_text(size = 10),
201 legend.position = "bottom")
202
203 ggplot2::ggsave("SCBD_Species_Radar_plot.png", rad_uvc)
204 }
205
206 ## LCBD
207
208 if (map) {
209 #Beta diversity
210 beta_div_ext(data_beta, data_xy, data_env)
211 }
212
213 #Lcbd per places and time
214 compute_lcbd(data_beta, data_spe, data_num)
215
216 #Lcbd of your choice
217 if (sepa) {
218 compute_lcbd2(data_beta, data_spe, data_num)
219 }
220
221 ##SCBD
222
223 scbd_uvc_tc <- make_scbd_uvc(data_spe, z = "TC", data_beta)
224
225 radar_plot(scbd_uvc_tc)