GRUMP Workshop | Explore Global Ocean Plankton Data in R

01 · Prepare

Install and load the packages.

Install the latest version of R from CRAN first, then install the latest free version of RStudio Desktop from Posit. This workshop uses the tidyverse for data wrangling and plotting, plus rnaturalearth and sf for the map.

install.packages("tidyverse")
install.packages("rnaturalearth")
install.packages("rnaturalearthdata")
install.packages("sf")

library("tidyverse")
library("rnaturalearth")
library("rnaturalearthdata")
library("sf")

02 · Load data

Set the working directory.

Download the GRUMP data from CMAP or Zenodo and place it on your Desktop in GRUMP-Workshop/02-Data. Name the file Grump_data.csv.

setwd("~/Desktop/GRUMP-Workshop/02-Data/")
grump_data <- read_csv("Grump_data.csv")

str(grump_data)

03 · Browse taxonomy

Find your plankton.

First, run this code in R to build a smaller taxonomy table. It contains one row for each ASV and makes it easier to inspect the organisms represented in GRUMP.

grump_taxonomy <- grump_data %>%
  select(
    Domain, Supergroup, Division, Phylum, Class,
    Order, Family, Genus, Species, ProPortal_ASV_Ecotype,
    Sequence_Type, ASV_hash, ASV
  ) %>%
  group_by(ASV_hash, Domain) %>%
  distinct(ASV_hash, .keep_all = TRUE)

Search all GRUMP taxonomy levels

Type a taxon such as Vibrio, Prochlorococcus, or Dinoflagellata. Results show the matching taxonomic level. Select the result you want and the R code below will update automatically.

Start typing to search the GRUMP taxonomy.

Your taxonomy chunk

Genus · Vibrio

Copy this personalized chunk into your R script immediately after the taxonomy-table code above, then run it. You can replace it at any time by searching for and selecting another taxon.

my_planktons_taxonomy <- grump_taxonomy %>%
  filter(Genus %in% c("Vibrio"))

my_plankton_data <- grump_data %>%
  filter(Genus %in% c("Vibrio"))

04 · Calculate abundance

Sum ASVs within each sample.

Many organisms have multiple ASVs in one sample. This step calculates total relative abundance for your selected taxon and keeps one row per sample. The taxonomic level in this code also updates when you select a search result above.

my_plankton_data <- my_plankton_data %>%
  group_by(SampleID, Genus) %>%
  mutate(
    Total_Relative_Abundance = sum(Relative_Abundance)
  ) %>%
  ungroup() %>%
  distinct(SampleID, .keep_all = TRUE)

05 · Make a map

Map global distribution.

Create a Robinson-projection world map, add bubbles scaled by relative abundance, and save the result as a PDF in the workshop figures folder.

world <- ne_countries(
  scale = "medium",
  returnclass = "sf"
) %>%
  st_set_crs(4326)

robinson <- paste0(
  "+proj=robin +lon_0=-130 +x_0=0 +y_0=0 ",
  "+ellps=WGS84 +datum=WGS84 +units=m +no_defs"
)

world_robinson <- world %>%
  st_break_antimeridian(lon_0 = -130) %>%
  st_transform(crs = robinson)

my_plankton_sf <- st_as_sf(
  my_plankton_data,
  coords = c("Longitude", "Latitude"),
  crs = 4326
) %>%
  st_transform(crs = robinson)

my_plankton_map <- ggplot() +
  geom_sf(
    data = world_robinson,
    fill = "grey90",
    color = NA
  ) +
  geom_sf(
    data = my_plankton_sf,
    aes(size = Total_Relative_Abundance),
    color = "#00b2f6",
    alpha = 0.7
  ) +
  scale_size_continuous(
    name = "Relative Abundance",
    range = c(2, 16)
  ) +
  ggtitle("My plankton in the ocean") +
  coord_sf(crs = robinson) +
  theme_minimal(base_size = 18) +
  theme(
    plot.title = element_text(
      size = 20, face = "bold", hjust = 0.5,
      margin = margin(b = 10)
    ),
    axis.title = element_blank(),
    axis.text = element_text(size = 14),
    legend.position = "bottom",
    legend.title = element_text(size = 16, face = "bold"),
    legend.text = element_text(size = 14),
    panel.grid = element_line(color = "grey85", size = 0.2)
  )

print(my_plankton_map)

setwd("~/Desktop/GRUMP-Workshop/03-Figures/")
ggsave(
  "my_plankton_map.pdf",
  plot = my_plankton_map,
  width = 11.7,
  height = 7
)

06 · Plot depth

Explore the upper 200 metres.

This function lets you select one or more ocean basins while retaining depth resolution and using a consistent abundance scale.

my_plankton_data_depth_filtered <- my_plankton_data %>%
  filter(depth <= 200)

global_max_abundance <- max(
  my_plankton_data_depth_filtered$Total_Relative_Abundance,
  na.rm = TRUE
)

plot_depth_profile <- function(
  ...,
  plot_title = "Depth Profile"
) {
  selected_cruises <- c(...)

  plot_data <- my_plankton_data_depth_filtered %>%
    filter(Ocean_Basin %in% selected_cruises)

  ggplot(
    plot_data,
    aes(
      x = lat,
      y = -depth,
      size = Total_Relative_Abundance,
      color = Ocean_Basin
    )
  ) +
    geom_point(alpha = 1) +
    scale_size_continuous(
      range = c(1, 20),
      name = "Relative Abundance",
      limits = c(0, global_max_abundance)
    ) +
    theme_minimal() +
    labs(
      title = plot_title,
      x = "Latitude",
      y = "Depth (m)",
      color = "Ocean_Basin"
    ) +
    theme(
      text = element_text(size = 16),
      axis.text.x = element_text(
        angle = 45, hjust = 1, size = 14
      ),
      axis.text.y = element_text(size = 14),
      axis.title.x = element_text(size = 18, face = "bold"),
      axis.title.y = element_text(size = 18, face = "bold"),
      plot.title = element_text(
        size = 20, face = "bold", hjust = 0.5
      ),
      legend.position = "bottom",
      legend.text = element_text(size = 14),
      legend.title = element_text(size = 16, face = "bold")
    )
}

all_ocean_basins_depth_profile <- plot_depth_profile(
  "Atlantic.Ocean",
  "Southern.Ocean",
  "Indian.Ocean",
  "Arctic.Ocean",
  "Pacific.Ocean",
  plot_title = "Depth Profile"
)

print(all_ocean_basins_depth_profile)

setwd("~/Desktop/GRUMP-Workshop/03-Figures/")
ggsave(
  "all_ocean_basins_depth_profile.pdf",
  plot = all_ocean_basins_depth_profile,
  width = 11.7,
  height = 7
)

07 · Compare provinces

Plot Longhurst biogeography.

Compare abundance across Longhurst provinces with a log-scaled box plot and save the final figure.

Longhurst_Box_Plot <- ggplot(
  my_plankton_data,
  aes(
    x = factor(Longhurst_Long),
    y = Total_Relative_Abundance,
    fill = Longhurst_Long
  )
) +
  geom_boxplot(alpha = 1, outlier.shape = NA) +
  geom_jitter(width = 0.2, size = 1, alpha = 0.5) +
  scale_y_log10() +
  theme_minimal() +
  labs(
    title = "My Plankton Abundance by Longhurst Province",
    x = "Ocean Basin",
    y = "Total Relative Abundance",
    fill = "Longhurst Province"
  ) +
  theme(
    text = element_text(size = 16),
    axis.text.x = element_text(
      angle = 90, hjust = 1, size = 14
    ),
    axis.text.y = element_text(size = 14),
    axis.title.x = element_text(size = 18, face = "bold"),
    axis.title.y = element_text(size = 18, face = "bold"),
    plot.title = element_text(
      size = 20, face = "bold", hjust = 0.5
    ),
    legend.position = "none"
  )

Longhurst_Box_Plot

setwd("~/Desktop/GRUMP-Workshop/03-Figures/")
ggsave(
  "My_Plankton_Longhurst_Box_Plot.pdf",
  plot = Longhurst_Box_Plot,
  width = 11.69,
  height = 8.27,
  units = "in"
)

Explore plankton in R.