-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path220726_UpdateGRL_BacteriaPlotsPatch.R
211 lines (165 loc) · 8.06 KB
/
220726_UpdateGRL_BacteriaPlotsPatch.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
########### Update EP grassland datasets 27706 and 27707############
#This script creates datasets xx version xx and xx version xx in Bexis
#Script by Caterina Penone
#This script is not a full update but a fix to two issues:
#1. add two more plots for bacteria 2014 (after a mistake in script 210112_UpdateGRLdataset.R LINE:282)
#2. add dung beetles 2014
#3. remove protists 2011_2012, OTU names in the 3 regions were not comparable, add one missing plot (SEG16 in 2011)
require(data.table)
setwd("N:/")
setwd("C:/Users/Caterina/Dropbox/")
source("R/SCRIPTS UTILES/BE_plots_zero.R")
#this script is also here: https://github.com/biodiversity-exploratories-synthesis/Synthesis_useful_functions
#Read last version of grassland dataset
grl <- fread("Exploratories/Data/GRASSLANDS/210826_EP_species_diversity_GRL_BEXIS27707.txt")
tr <- fread("Exploratories/Data/GRASSLANDS/210112_EP_species_info_GRL_BEXIS.txt")
grl2 <- merge(grl, tr, by="Species")
#1. Add two plots for bacteria 2014 ###################################################
# There is a mistake in script 210112_UpdateGRLdataset.R LINE:282
# Bacteria have 2 missing plots in 2011, not 2014, add these plots in 2014
# This script adds these plots and also changes type from "OTU_number" to "ASV_number"
# bacteria 2011
bac <- fread("Exploratories/Data/GRASSLANDS/TEXTfiles/190319_Update/24866.txt")
length(unique(bac$Plot_ID))*length(unique(bac$Sequence_variant)) #zeros are missing
bac$DataID <- 24866
bac$Dataversion <- "2"
bac$Year <- 2011
bac[, c("kingdom", "phylum", "class", "order", "family", "genus", "species") := tstrsplit(Taxonomy, ", ", fixed=TRUE)]
bac[,c(4,9:14):=NULL]
plotIDs <- data.frame(Plot_ID = unique(bac$Plot_ID))
plotIDs <- data.table(BEplotZeros(plotIDs,"Plot_ID",plotnam = "Plot"))
bac <- merge(bac,plotIDs,by="Plot_ID")
gc()
# bacteria 2014
bac4 <- fread("Exploratories/Data/GRASSLANDS/TEXTfiles/190319_Update/25066.txt")
length(unique(bac4$Plot_ID))*length(unique(bac4$Sequence_variant)) #zeros are missing
bac4$DataID <- 25066
bac4$Dataversion <- "3"
bac4$Year <- 2014
bac4[, c("kingdom", "phylum", "class", "order", "family", "genus", "species") := tstrsplit(Taxonomy, ", ", fixed=TRUE)]
bac4[,c(4,9:14):=NULL]
plotIDs <- data.frame(Plot_ID = unique(bac4$Plot_ID))
plotIDs <- data.table(BEplotZeros(plotIDs,"Plot_ID", plotnam = "Plot"))
bac4 <- merge(bac4, plotIDs, by="Plot_ID")
gc()
#merge taxonomy
setkey(bac4, Sequence_variant)
setkey(bac, Sequence_variant)
bac[bac4, kingdom := i.kingdom, by=.EACHI]
bac <- rbind(bac, bac4)
length(unique(bac$Plot)) #150
rm(bac4, plotIDs); gc()
#change column names
setnames(bac, c("Sequence_variant","Read_count","Plot_ID","kingdom"),
c("Species","value","Plot_bexis","Group_fine"))
#add columns
bac$type <- "ASV_number"
#add info on trophic level etc.
names(grl2)
bac$Group_broad <- bac$Trophic_level <- bac$Fun_group_broad <- bac$Fun_group_fine <- "bacteria.RNA"
###check if two OTUs have same "name" but different taxonomy
temp <- unique(bac, by=c("Species","Group_fine"))
length(unique(bac$Species)) == nrow(temp)
rm(temp)
#remove old bact and rbind with main dataset
grl2 <- grl2[!Group_broad %in% "bacteria.RNA"]
sort(unique(grl2$Group_broad))
setdiff(names(grl2), names(bac))
grl2 <- rbind(grl2, bac, use.names=T)
rm(bac); gc()
#######################################################################################
#2. add dung beetles species ##########################################################
# 21207: Dungwebs Species List 2014 & 2015 (Invertebrates, Scarabaeoidea, Dung Beetles) - downloaded 27/06/18
dung <- fread("Exploratories/Data/GRASSLANDS/TEXTfiles/190319_Update/21207.txt")
dung <- dung[!grepl("W", dung$EP)] #remove forest data
length(unique(dung$EP)) #150
summary(dung$count) #dataset has zeros
unique(dung$month) #use only august, july june 2014 (rest is VIP)
unique(dung$dungtype) #use only cow, deer, fox, horse, sheep, wildboar
#see how many are already in the dataset to decide whether to include the dataset
dungsp <- unique(dung$species)
setdiff(dungsp, unique(grl2$Species)) #none of the 34 sp is in the dataset
rm(dungsp)
dung <- dung[month %in% c("August_2014","July_2014","June_2014")]
length(unique(dung$EP)) #150
#drop species that are only in forests
dung[,temp:=sum(count), by="species"]
summary(dung)
dung <- dung[temp>0]
unique(dung$species) #17 species now
length(unique(dung$species)) * length(unique(dung$EP)) #need to aggregate months
dung$temp <- NULL
#remove non needed columns
dung$date <- dung$month <- dung$exotic <- dung$chem <- dung$dungtype<-NULL
#aggregate months
dung[,value:=sum(count), by=c("species", "EP")]
dung <- unique(dung[, c("species", "EP", "value"), with=F])
length(unique(dung$species)) * length(unique(dung$EP)) == nrow(dung)
#homogenise names
setnames(dung, c("species", "EP"), c("Species","Plot_bexis"))
#add missing columns
dung$Year <- 2014
dung$type <- "abundance"
dung$DataID <- 21207
dung$Dataversion<-"2"
dung <- data.table(BEplotZeros(dung,"Plot_bexis",plotnam="Plot"))
sort(unique(grl2$Group_broad))
dung$Group_broad<-"Coleoptera"
dung$Group_fine<-"Scarabaeidae"
sort(unique(grl2$Trophic_level))
dung$Trophic_level<-"decomposer.arthropod"
sort(unique(grl2$Fun_group_broad))
dung$Fun_group_broad<-"decomposer.dungbeetle"
sort(unique(grl2$Fun_group_fine))
dung$Fun_group_fine<-"decomposer.dungbeetle"
length(unique(dung$Plot)) #150
#add to grassland dataset
setdiff(unique(dung$Species), unique(grl2$Species))
grl2 <- rbind(grl2, dung, use.names=T)
rm(dung)
#######################################################################################
#3. fix issue with protists 2011_2012 and add NA for SEG16 in 2011, dataset 24426 #####
#The 2011_2012 dataset of ""mainlybacterivore.protist" has two issues:
# 1-regions are treated separately in the pipeline
# 2- Acanthoamoeba are anecdotical and Myxomycetes are a group hard to classify (can switch from amoeba to flagellates)
# For those reasons, in particular issue 1, with Anna-Maria Fiore-Donno we decided to remove them from the dataset
unique(grl2$Year)
unique(grl2[Year=="2011_2012"]$Trophic_level) #safe to remove based on year
#which datasets ID to remove from metadata?
unique(grl2[Year=="2011_2012"]$DataID) #18166 18206 18226 18187 18208 18207
grl2 <- grl2[!Year=="2011_2012"]
grl2[Trophic_level=="mainlybacterivore.protist"] #ok
#Add missing plot: the 2011 dataset has one plot missing (SEG16 in 2011 for dataset 24426 )
length(unique(grl2[Group_broad=="Protists" & Year=="2011"]$Plot)) #149
pro <- grl2[Group_broad=="Protists"]
length(unique(pro[Year=="2011"]$Species))
length(unique(pro[Year=="2017"]$Species))
proadd <- pro[Year=="2011" & Plot=="SEG17"] #use another plot as model, ok this as 2004 species
proadd$Plot_bexis <- proadd$Plot <- "SEG16"
proadd$value <- NA
#add to main dataset
grl2 <- rbindlist(list(grl2, proadd))
length(unique(grl2[Group_broad=="Protists" & Year=="2011"]$Plot)) #150
rm(pro, proadd); gc()
#######################################################################################
###########Save the diversity and characteristics tables separately ###################
grl <- grl2[,.(Plot_bexis, Plot, Species, value, type, DataID, Year, Dataversion)]
tr <- grl2[,.(Species, Group_broad, Group_fine, Trophic_level, Fun_group_broad, Fun_group_fine)]
length(unique(tr$Species))
dim(unique(tr)) #no duplicates
length(unique(grl2$Species))
tr <- unique(tr)
sum(is.na(grl))
apply(grl,2,function(x)sum(is.na(x)))
unique(grl2[is.na(Dataversion)]$Trophic_level) #NAs in dataversion are normal (birds and bats)
sum(is.na(tr))
apply(tr,2,function(x)sum(is.na(x)))
tr[is.na(Fun_group_broad)] #two myriapods and one plant
unique(tr[Trophic_level=="autotroph"]$Fun_group_broad) #ok to keep NA, it's unkown anyway
unique(tr[Trophic_level=="decomposer.myriapod"]$Fun_group_broad) #ok to keep NA
summary(factor(tr$Trophic_level))
#reorder column names in grl
setcolorder(grl,c("Plot_bexis","Plot","Species","value","type","Year","DataID","Dataversion"))
#all good: save
fwrite(grl,"Exploratories/Data/GRASSLANDS/221102_EP_species_diversity_GRL_Patch22.txt",row.names=F,quote=F,sep=";",na=NA)
fwrite(tr,"Exploratories/Data/GRASSLANDS/221102_EP_species_info_GRL_patch22.txt",row.names=F,quote=F,sep=";",na=NA)