20191029_SciData.Rmd

---
title: "Global nematode dataset"
author: "J. van den Hoogen et al., Scientific Data, 2020"
output:
  github_document: default
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load packages, results='hide', message=FALSE, warning=FALSE}
suppressPackageStartupMessages({
  library(cowplot)
  library(scales)
  library(plyr)
  library(reshape2)
  library(RColorBrewer)
  library(data.table)
  library(tidyverse)
})
```

# Initial data formating and collating
```{r Data import & general calculations}
# Load full dataset, including WWF_Biome and Pixel_Lat, Pixel_Long
full_data_wBiome <- fread("data/nematode_full_dataset_wBiome.csv") %>% 
  mutate(WWF_Biome = round(WWF_Biome, 0))

# Set WWF_Biome to 25 as arbitrary biome for Antarctic sites
full_data_wBiome[full_data_wBiome$Latitude < -60,]$WWF_Biome <- 25

# Aggregate onto pixel level
aggregated <- full_data_wBiome %>%  
  select(Bacterivores, Fungivores, Herbivores, Omnivores, Predators, Unidentified, Total_Number, Pixel_Lat, Pixel_Long) %>% 
  group_by(Pixel_Lat, Pixel_Long) %>% 
  summarise_all(mean) %>% 
  filter(!is.na(Pixel_Lat))

# Write to csv
# Uncomment if need be
# write.csv(aggregated, "data/nematodes_aggregated.csv", row.names = F)

# This dataset ("nematodes_aggregated.csv") was sampled to pull covariate data
sampled_data <- read.csv("data/nematode_aggregated_wCovariateData.csv") %>%
  mutate(WWF_Biome = round(WWF_Biome,0)) 

Data_tot <- sampled_data %>% 
  select(Total_Number, Bacterivores, Fungivores, Herbivores, Omnivores, Predators) %>% 
  pivot_longer(everything(),
               names_to = "Group",
               values_to = "Count")
```

```{r}
full_data_summary <- full_data_wBiome %>% 
  select(Bacterivores, Fungivores, Herbivores, Omnivores, Predators, Total_Number) %>% 
  pivot_longer(everything(),names_to = "Group", values_to = "Count") %>% 
  na.omit() %>% 
  group_by(Group) %>% 
  summarise(mean = mean(Count), median = median(Count), n = n()) %>% 
  mutate(mean = round(mean),median = round(median))

full_data_summary

# Write to file
write.csv(full_data_summary,"output/table1_full_data_sum.csv") 
```

# Sampling locations
```{r Observations per biome}
# Plot sampling points on map
nematode_pointmap <- ggplot() + geom_polygon(data = map_data("world"), 
                                             aes(x = long, y = lat, group = group),
                                             fill = "#bababa",
                                             color = NA,
                                             size = 0.1) + 
  coord_fixed(1.1) +
  geom_point(data = full_data_wBiome,
             aes(x = Longitude, y = Latitude),
             fill = "red",
             color = "black",
             pch = 21
  ) +
  scale_fill_gradientn(colors = brewer.pal(8, "YlOrRd"),
                       limits = c(0, 4000),
                       oob = scales::squish,
                       name = "Nematodes per 100g dry soil") +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.box="horizontal",
        panel.grid = element_blank(),
        axis.title=element_blank(),
        axis.text=element_blank()) +
  guides(fill = guide_colorbar(title.position = "top"))

nematode_pointmap

# Save to file
ggsave("output/figXX_nematode_pointmap.pdf", plot = nematode_pointmap)
```

# Summary per biome
```{r}
Data_biome <- full_data_wBiome %>% 
  select('Total_Number','Bacterivores','Fungivores',"Herbivores",'Omnivores','Predators','WWF_Biome') %>% 
  mutate(WWF_Biome = round(WWF_Biome)) %>% 
  filter(WWF_Biome != 98) %>%
  filter(WWF_Biome != 0) %>% 
  # na.omit() %>% 
  reshape2::melt(id.vars = "WWF_Biome") %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 1, "Tropical Moist Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 2, "Tropical Dry Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 3, "Tropical Coniferous Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 4, "Temperate Broadleaf Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 5, "Temperate Conifer Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 6, "Boreal Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 7, "Tropical Grasslands")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 8, "Temperate Grasslands")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 9, "Flooded Grasslands")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 10, "Montane Grasslands")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 11, "Tundra")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 12, "Mediterranean Forests")) %>%
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 13, "Deserts")) %>% 
  mutate(WWF_Biome = replace(WWF_Biome, WWF_Biome == 25, "Antarctica")) %>% 
  setNames(., c("Biome","Group","Count")) 

Data_biome_sum <- Data_biome %>%
  filter(Group == "Total_Number") %>% 
  group_by(Biome) %>% 
  dplyr::summarize(mean = mean(Count, na.rm=TRUE), 
                   median = median(Count, na.rm=TRUE),
                   n = n()) %>% 
  arrange(desc(median))  %>% 
  mutate(mean = round(mean), median = round(median))

Data_biome_sum

# Write to file
write.csv(Data_biome_sum, "output/table2_biome_sum.csv")
```

# Frequency plot per method
```{r}
method_plot_data <- full_data_wBiome %>% 
  group_by(sampling_ref) %>% 
  summarise(n = n()) %>% 
  na_if("") %>% na.omit() %>% 
  mutate(sampling_ref = replace(sampling_ref, sampling_ref == "Baermann", "Baermann funnel")) %>% 
  mutate(sampling_ref = replace(sampling_ref, sampling_ref == "Cobb", "Decanting and sieving")) %>% 
  mutate(sampling_ref = replace(sampling_ref, sampling_ref == "Jenkins/Freckman", "Sugar floatation")) %>% 
  mutate(sampling_ref = replace(sampling_ref, sampling_ref == "Oostenbrink", "Oostenbrink elutriator ")) %>% 
  mutate(sampling_ref = replace(sampling_ref, sampling_ref == "Seinhorst", "Seinhorst elutriation")) %>% 
  mutate(sampling_ref = replace(sampling_ref, sampling_ref == "Whitehead", "Whitehead tray"))

methods_barplot <- ggplot((method_plot_data %>% mutate(sampling_ref = fct_reorder(sampling_ref, desc(n)))), aes(x = sampling_ref, y = n)) +
  geom_bar(stat = "identity", fill = "#328da8", alpha = 0.6, width = 0.6) +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(),
    axis.ticks = element_line(),
    plot.title = element_text(hjust = 0.5, size = 10),
    axis.text.x = element_text(hjust = 1, angle = 60),
    legend.position = "none") +
  ylab("Number of samples") +
  xlab("")

methods_barplot

# Save to file
save_plot(methods_barplot, filename = "output/methods_barplot.pdf", device = "pdf")
```

```{r}
plot_data <- sampled_data %>% 
  select(Aridity_Index,WorldClim2_Annual_Precipitation, WorldClim2_Annual_Mean_Temperature, SG_Soil_pH_H2O_015cm, 
         SG_SOC_Content_015cm, Human_Development_Percentage, SG_CEC_015cm, SG_Sand_Content_015cm, Npp, EVI) %>% 
  mutate(SG_Soil_pH_H2O_015cm = SG_Soil_pH_H2O_015cm/10) %>% 
  rename("Soil pH" = SG_Soil_pH_H2O_015cm,
         "SOC content (g/kg)" = SG_SOC_Content_015cm,
         "Mean annual temperature (°C)" = WorldClim2_Annual_Mean_Temperature,
         "Annual precipitation (mm)" = WorldClim2_Annual_Precipitation,
         "Human development (%)" = Human_Development_Percentage,
         "Aridity Index" = Aridity_Index,
         "Cation exchange capacity (cmolc/kg)" = SG_CEC_015cm,
         "Sand content (%)" = SG_Sand_Content_015cm,
         "Net Primary Productivity" = Npp,
         "Enhanced Vegeation Index" = EVI) %>% 
  pivot_longer(cols = everything(),
               names_to = "Variable",
               values_to = "Value") %>% 
  na.omit() 


plot_overall <- ggplot(data = plot_data, aes(Value)) +
  geom_density(position = "stack", 
               aes(fill = "red",
                   alpha = 0.25,
                   y = ..scaled..)) +
  facet_wrap(vars(Variable), 
             scales = "free_x", 
             nrow = 2, 
             strip.position = "bottom") +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(),
    axis.ticks = element_line(),
    # plot.title = element_text(hjust = 0.5, size = 10),
    legend.position = "none",
    strip.text = element_text(size = 8),
    strip.placement = "outside",
    aspect.ratio = 1) +
  ylab("Relative frequency") +
  xlab("")

plot_overall

save_plot("output/sampled_range.pdf", plot = plot_overall, base_height = 5, base_width = 12)

```