Find species occurrences within a region

December 30, 2017January 6, 2018
by brian_ward

Over the past few weeks I have been negligent in my duties to this website; I attribute my lapse in activity (of all kinds) to the holidays. However, I haven’t been totally sedentary. I was just home in New Hampshire for a few days, where we received 6 inches of fresh snow on Christmas day, on top of the foot or so that was already on the ground. I went out for a snowshoe with my parents and we were truly in a winter wonderland:

winter hunting cabin — The 2017-2018 winter in New Hampshire, off to a good start

However, the downside is that temperatures dropped precipitously after that, eventually staying in the vicinity of 0° F (-18° C) for the remainder of my trip. The moral of the story is that while it’s nice to visit New Hampshire in the winter, it’s also nice to be able to leave.

Anyways, I guess that’s enough about me. On to the primary topic of this post:

Problem

A coworker recently posed the following problem that she needed to solve in an extensible way:

Given a country X, find all species of weeds that have been reported there, and compare this list against a preexisting list containing the names of weed species of concern

In order to solve this problem, the first resource I would turn to is the Global Biodiversity Information Facility (GBIF). GBIF hosts a massive database that collects observations of species (referred to as “occurrences”) from all over the world; as of this writing the database contains roughly 1 billion species occurrences total.

Outline

We can break this problem down into a few simple steps:

Clean our preexisting list of species names to match GBIF formatting. In the example I’ll use, we will have to get rid of some extra information included with the species names. One example is naming authorities – so for our purposes the gray wolf, first named by Linnaeus in 1758, should be written as Canis lupus instead of Canis lupus L. In addition, we will remove information regarding subspecies, varieties, and cultivars.
Identify an appropriate dataset to download. We already know that we need to download species occurrences for a particular country, but downloading all occurrences from any given country could be far more information than we need. Therefore, we need to restrict the size of the download as much as is feasible.
Initiate the download from GBIF servers
Clean species names within the downloaded GBIF data. Generally the species names in GBIF datasets are formatted well, though it pays to also clean this data to ensure reliable matches with our own list.
Compare our list of species names with the GBIF data
[Optional] write out the downloaded occurrence data to a .csv file

R example

The script that this code is copied from is available here.

For the example, I will be using the excellent rOpenSci package rgbif, which, you may have guessed, provides several R functions for accessing GBIF data. If you don’t have rgbif installed, just open an R terminal and type:

install.packages("rgbif")

Input variables and data

Now that that is out of the way, let’s start by loading required packages, defining some variables, and simulating some data. We’ll use the country of Chad as an example. For our list of species, I am just going to create a data frame containing some various plants along with two animals with scientific names that I find humorous. Normally, our list of species would probably come from a preexisting .csv file that we read in. Finally, we need supply a model or type species, and a desired taxonomic rank of interest (more on that below).

library(rgbif)

## Country of interest
country <- "Chad"

## Species of interest - would normally obtain this from a pre-existing .csv file
species <- data.frame("sci_name" = c("Pica pica L.", 
                                     "Tuta absoluta (Meyrick)", 
                                     "Arachis hypogaea var. Georgia-12Y", 
                                     "Vigna vexillata (L.) A. Rich", 
                                     "ZEA MAYS"),
                      "common_name" = c("Eurasian magpie",
                                        "tomato leafminer",
                                        "peanut variety Georgia-12Y",
                                        "wild cowpea",
                                        "maize"))

## Model or type species, and taxonomic rank of interest
mod_sp <- "Arabidopsis thaliana"
rank <- "phylum"

## Select whether to output occurrence data (TRUE/FALSE)
## And specify output file name (.csv extension required)
## File name irrelevant if output_occs not set to TRUE
output_occs <- FALSE
out_file_name <- "Chad_vascular_plant_occurences.csv"

Define Functions

Now we need to create a couple of functions. The first will read in a vector of species names, and will output a cleaned-up version of this vector. Specifically, it will only keep the genus and species names, and will ensure that the genus name is capitalized (and that the species name is not). It will also get rid of any additional information such as naming authority, variety, subspecies, etc.

clean_sp <- function(s) {
  
  s <- gsub("_", " ", s, fixed = TRUE)
  s <- tolower(s)
  l <- strsplit(s, " ", fixed = TRUE)
  l <- lapply(l, function(x) {paste(x[1], x[2], sep = " ")})
  cleaned <- unlist(l)
  
  substr(cleaned, 1, 1) <- toupper(substr(cleaned, 1, 1))
  
  return(cleaned)
}

The second function will take a species that we choose as input, as well a taxonomic rank of our choosing (either kingdom, phylum, class, order, family, or genus). The GBIF functions that we use below require a “key” (which is an integer) for identifying a specific taxonomic rank. We can get these keys using the rgbif function name_backbone(). For example, if we want to find all occurrences of vascular plants (phylum Tracheophyta) within a country, the easiest way to get the GBIF key for Tracheophyta is to look up an example plant species that we know is in this phylum using the name_backbone() function, and then extract the key from it. The phylum key for Tracheophyta happens to be 7707728 – it doesn’t exactly roll off the tongue. The function tax_key() below is a wrapper for name_backbone() which allows us to easily look up a rank key for any given species.

tax_key <- function(species, rank) {
  
  if(!rank %in% c("kingdom", "phylum", "class", "order", "family", "genus")) {
    stop("rank argument requires value of 'kingdom', 'phylum', 'class', 'order', 'family', or 'genus'")
  }
  
  tax <- name_backbone(species)
  key <- unlist(tax[grep(paste(rank, "Key", sep = ""), names(tax))])
  print(paste("Using GBIF key", key, "for", rank, tax[[rank]]))
  
  return(key)
}

Executable example

Now we’re finally able to run some code. First we’ll clean up our species names, obtain the key value for phylum Tracheophyta (using the species Arabidopsis thaliana as an example), and obtain the two-letter ISO 3166-1 country code for Chad:

## Clean up input vector of species names
species$clean_name <- clean_sp(species$sci_name)

## Obtain integer key for given species/rank combination
key <- tax_key(species = mod_sp, rank = rank)

## Obtain ISO 3166-1 two-letter country code
country_code <- rgb_country_codes(country)$code

The rgbif package includes three functions that can download data from GBIF servers:

occ_download() – performs exactly the same operation as a manual data download. Requires GBIF username and password (which can be obtained for free through the GBIF website)
occ_search() – allows for finer control of download queries. Can download data, metadata, or hierarchical trees
occ_data() – faster version of occ_search() optimized for downloading data only

For the rest of this example, we’ll focus on the later two functions, though using occ_download() is a good idea if you will be grabbing a large amount of data.

Before we download any occurrence data, it is always a good idea to check on how many results there are for our desired search query. We can do this using the occ_search() function, by telling it to return metadata, and setting the limit value to 0. You’ll notice below that I’m calling occ_search() in an unusual fashion, by using the function do.call(). do.call() allows us to specify a function, and then supply its arguments as a list. The reason I’m doing this has to do with our key value for phylum Tracheophyta. Normally if I called occ_search() and wanted to only select this phylum, I would need to supply the argument manually:

occ_search(phylumKey = 7707728, …).

However, if on another occasion I wanted to restrict my search to the family Brassicaceae (key 3112), I would need to use the argument “familyKey”:

occ_search(familyKey = 3112, …)

This means that I would have to manually change the function call every time I search by a different taxonomic rank. The tax_rank() function that we defined earlier returns a named integer for the key, so we can supply its name (“phylumKey” in this example) in the list supplied to do.call() to avoid having to ever manually change the occ_search() call.

So, in the block below I define a list, gbif_args, which contains some arguments that we will use repeatedly. Notice “key” is listed without any assignment (=), because it contains both a name and a value. When I actually call the occ_search() function below, I append on some additional arguments to return metadata, and to limit the number of occurrences to 0. Note that I also include a prompt so the user can select whether or not to proceed with the download in the next step. The script must be sourced (i.e. not run interactively) for it to kick in:

## Construct list of arguments for GBIF functions
gbif_args <- list(country = country_code,
                  key,
                  hasCoordinate = TRUE,
                  hasGeospatialIssue = FALSE,
                  curlopts = list(progress = TRUE))

## Find number of occurences using specified arguments
n_occs <- do.call(occ_search, c(gbif_args, list(return = "meta", limit = 0)))$count
print(paste("Number of occurences returned by specified parameters:", n_occs))
y_n <- readline(prompt = "Proceed with occurrences download? [y/n]")
if (y_n != "y") {
  stop("User declined to download occurrences, or did not supply a y/n response")
}

We can now finally download our species occurrence data for real. We have our choice of either the occ_search() or occ_data() functions. Of the two, occ_data() tends to be faster. Running the code below on my computer/internet connection, occ_search() downloaded 3,225 occurrences in 0.64 minutes, while occ_data() completed the same task in 0.24 minutes.

## Download occurences into data frame using occ_search() function
start_time <- Sys.time()
occs <- do.call(occ_search, c(gbif_args, list(return = "data", limit = n_occs)))
end_time <- Sys.time()

print(paste("Elapsed time to download", nrow(occs), "occurences:", 
            round(difftime(end_time, start_time, units = "mins"), 2), "minutes"))

## Download occurences into data frame using occ_data() function
## Note "return" argument is not used by occ_data()
start_time <- Sys.time()
occs <- do.call(occ_data, c(gbif_args, list(limit = n_occs)))$data
end_time <- Sys.time()

print(paste("Elapsed time to download", nrow(occs), "occurences:", 
            round(difftime(end_time, start_time, units = "mins"), 2), "minutes"))

Finally, we can compare the GBIF species occurrences with our own list of species. Some GBIF occurrences are only identified to the genus level, so I’ll remove these by deleting any occurrences with scientific names that don’t contain a space. While the species names in the GBIF data tend to be pretty clean, we might as well run them through our cleaning function to be thorough. While we’re at it, we might as well save the occurrence data to a file, though this step is optional.

occs <- occs[grep(" ", occs$name, fixed = TRUE), ]
occs$clean_name <- clean_sp(occs$name)

intersect(species$clean_name, occs$clean_name)

if (output_occs) {
 write.csv(occs, file = out_file_name, row.names = FALSE)
}

When we run intersect(), we find that the three plant species in our original dataframe (which were not at all cherry-picked) are indeed among the list of vascular plants that have been observed in Chad. Critically, the Eurasian magpie and tomato leafminer are not shown here, since they are not vascular plants (or plants of any kind for that matter).

Going further

The rgbif package allows for many customized operations. For instance, we can isolate occurrences from a specific state or province. In this example, I count the number of occurrences of vascular plants in the U.S. state of Delaware:

## Count all occurences in US state of Delaware
gbif_args <- list(country = "US",
                  stateProvince = "Delaware",
                  key,
                  hasCoordinate = TRUE,
                  hasGeospatialIssue = FALSE,
                  limit = 0,
                  return = "meta",
                  curlopts = list(progress = TRUE))
n_occs <- do.call(occ_search, gbif_args)$count

If country borders seem a little too arbitrary, you can also define a custom polygon in well known text (WKT) format. I recommend this site for constructing the polygon (the site only seems to work with MS Internet Explorer or MS Edge). Below we search for vascular plants within a rectangle covering portions of Chad, Cameroon, and Nigeria:

## Find all occurences within a polygon
library(wicket)
wkt <- "POLYGON((17.02 11.52, 10.47 11.03, 10.78 7.78, 17.20 8.57, 17.02 11.52))"
wkt <- wkt_correct(wkt)
gbif_args <- list(geometry = wkt,
                  hasCoordinate = TRUE,
                  hasGeospatialIssue = FALSE,
                  limit = 100,
                  curlopts = list(progress = TRUE))
occs <- do.call(occ_data, gbif_args)$data

Cautionary notes

This is a pretty powerful method for finding species occurrences within a specific region. However, there are a few rules you should keep in mind:

Reduce, reduce, reduce. You want your queries to be as specific as possible while still being relevant. Otherwise you’ll end up wasting time downloading lots of irrelevant data. I recommend using the lowest taxonomic rank and smallest geographic region possible. In this example, we downloaded vascular plant occurrences from Chad. If I had searched for all vascular plant occurrences in China, or the United States, the download would take a long time, or would just fail outright. That brings us to the next rule:
Be nice to GBIF. They are providing a very valuable free service. Please don’t use scripts like this to spam their servers with multiple huge download requests.

Anyways, that’s about it. Now you can search out species occurrence data to your heart’s content.

R coding

Multi-Panel Figures in R

November 9, 2017November 12, 2017
by brian_ward

Background

This post will detail how to make some nice, (nearly) publication-ready, multi-panel figures in R. I want to stress that this isn’t a post intended for R beginners. I will be mixing some standard R code with tidyverse syntax, using the latter’s almighty pipe operator (%>%) and some left-to-right assignment (->).

As the title suggests, this post will be about that most dreaded part of preparing a manuscript for publication – making multi-panel figures. These are often an intimidating part of writing any article; if you’re too much of a perfectionist like me, it can be hard enough to make a single graph that looks just right, let alone making multiple graphs and fitting them all together.

When I first started using R, I was at a loss when it came time to put multiple graphs together. At first I would just create them individually, pasting them together afterwards in image editing software. Unfortunately this technique often leads to things that don’t quite look right – panels that don’t exactly line up, lines/points/labels that are different sizes, etc.

Eventually I learned about packages such as lattice and ggplot, which allow you to take the data stored within a dataframe, split it up into pieces based upon a factor variable, and then plot each piece. So, here’s an example using lattice and some car engine data from the famous Motor Trend dataset:

library(lattice)
data(mtcars)
xyplot(hp ~ disp | factor(cyl), pch = 16, data = mtcars)

The xyplot() function is plotting engine displacement on the x-axis against horsepower on the y-axis, and splitting the data up by number of cylinders. The result:

Plotting engine displacement, horsepower, and cylinders with the lattice package

Here’s how to make approximately the same thing using ggplot, and the output:

library(ggplot2)
data(mtcars)
ggplot(mtcars, aes(x = disp, y = hp)) +
  geom_point(color = "blue", fill = "blue") +
  facet_grid(~ cyl) +
  theme_bw()

Plotting engine displacement, horsepower, and cylinders with the ggplot package

Not too bad – the ggplot code is a little more involved, but offers nearly endless options for customization. I want to emphasize nearly, because there are certain things that ggplot and lattice can do with multiple panels, and certain things they can’t. If you’re ever using ggplot and you think to yourself “this graph would be great if I could just change the scale on only the left pane” you’ll quickly learn that it’s relatively difficult to do. This is by design – “customized” multi-pane figures are one way in which people commonly try to cook up some significant looking findings.

In the example below, I’ll show you how to mix and match different plots together. I want to emphasize that this is often a bad idea, but just like everything in life, there are exceptions. Ultimately it is up to you to decide when it is appropriate to use this technique. However, even if you are working with graphs that aren’t drastically different, this can be a good method to let you format and customize each one individually to get everything looking good. I’ll first go through the example, and at the end give my justification for why I think the practice is okay in this situation.

For the example, I will create a dataset consisting of 5 sunflower genotypes tested across 4 environments (props if you can figure out which state I pulled these town names from). I will simulate two traits for all genotype/environment combinations – plant height, and seed oil content. Then I’ll make genotype-by-environment interaction plots for both traits. The code will use ggplot in conjunction with the cowplots package to arrange plots into a nice, professional-looking grid.

Data Simulation

The first step is to create our simulated data. I’ll define a function, dat_sim(), and then use it to create a dataset “gxe”, consisting of the environment name, genotype name, and then the oil and height data. The dat_sim() function will create random data defined by both the standard deviation between environments, and the standard deviation between genotypes within environments. So, I’ll give the oil content trait some messy genotype-by-environment interaction – the trait variation will be much greater across environments than across genotypes. Conversely, I’ll give the height trait highly stable interaction, with trait variation across environments much less than the variation across genotypes:

library(tidyverse)
library(cowplot)
set.seed(99)

dat_sim <- function(n_genos, n_envs, grand_mean, env_sd, geno_sd) {
  vec <- abs(rnorm(n_envs, mean = grand_mean, sd = env_sd))
  out_list <- list()
  for (i in 1:length(vec)) {
    out_list[[i]] <- rnorm(n_genos, vec[i], sd = geno_sd)
  }
  return(do.call(c, out_list))
}

gxe <- data.frame("env" = sort(rep(c("Alpena", "Mellette", "Kranzburg", "Centerville"), 5)),
                  "geno" = rep(c("Zebulon", "Chocolate", "Sunbright", "Firecracker", "Greenburst"), 4),
                  "oil" = dat_sim(5, 4, 0.25, 0.08, 0.03),
                  "height" = dat_sim(5, 4, 180, 30, 10))

One thing you’ll notice is that the trait units here are completely different. The oil content is expressed as a fraction, and the height data is expressed in centimeters.

The Plotting Function

Now I’m going to create a quick-and-dirty plotting function that uses ggplot() internally. This function, plot_gen(), will have input arguments consisting of the gxe dataframe (df), the trait name (trait), and a label for the trait (ylabel). It will then perform a few steps:

Find the mean trait value of all genotypes within each environment
Sort the environments from the lowest mean to the highest mean
Sort the dataframe by environment to match this ordering
Plot a regression for each genotype across the environments

Keep in mind that we will be running this function separately for each trait.

plot_gen <- function(df, trait, ylabel) {
  
  ## Find the order of environment means
  group_by(df, env) %>%
    summarise_at(.vars = trait, .funs = mean) %>%
    arrange_(trait) %>%
    pull(env) %>%
    as.character() ->
    env_order

  ## This assigns the env variable levels in the order of
  ## their mean values before sorting
  gxe$env <- factor(gxe$env, levels = env_order)
  gxe <- arrange(gxe, desc(env), geno)

  ## Plot regressions
  ggplot(gxe, aes_string(x = "env", y = trait)) +
    geom_smooth(aes(group = geno, color = geno), size = 1, method = "lm", se = FALSE) +
    xlab("Environment") +
    ylab(ylabel) +
    theme_bw() ->
    plot_out
  
  return(plot_out)
}

IMPORTANT: Notice that when I sort the dataframe in line 14 above, I am first sorting by environment, and then by genotype. In this case, this is a critical step. The genotype column in our dataframe must be sorted the same way for each plot that is created, because ggplot will read the genotype names in the order that it encounters them. If the ordering of this column doesn’t stay consistent, the coloring of the genotypes will change between plots. It doesn’t matter what order the environments are in, as long as the genotypes are always listed in the same order within each environment.

The method that I’m using – sorting environments by their mean values, and then performing a regression for each genotype – is similar to a method for more formally assessing genotype-by-environment interaction, first described by Finlay and Wilkinson (1963).

The advantage to doing our plotting inside a function is that it makes it easy to reproduce the same plot for different traits. So, to make the plot for oil content (our low-heritability trait), we can just write:

plot_gen(gxe, "oil", "Oil Content")

Which gives us:

Sunflower oil content genotype-by-environment interaction plot

Nice. Most of the genotypes are behaving alike across environments. There’s just that one troublemaker -the genotype ‘Greenburst’ (which conveniently got colored green), which is exhibiting strong crossover interaction with all the other genotypes.

Assembling the Plot

Making our oil content plot is all well and good, but now we can get fancier by plotting the two traits together. This is where the cowplots package comes into play. It contains a function, plot_grid(), which allows us to easily smoosh together plots by supplying them as arguments. You can also specify the layout, and give each plot labels. So, to plot our two traits together, it’s as simple as:

plot_grid(plot_gen(gxe, "oil", "Oil Content"),
          plot_gen(gxe, "height", "Plant Height (cm)"),
          ncol = 2,
          labels = c("A", "B"))

The output:

So… it’s not great. The big problem is that the legend appears twice. This is where it comes in handy to have all our plotting wrapped up in our plot_gen() function. First, to show an example of how we could tweak the plots themselves, I’ll rotate the x-axis labels 45 degrees, just for the heck of it. This just requires adding an additional line to our call to ggplot() – line 22 in the box below, so that plot_gen() becomes:

plot_gen <- function(df, trait, ylabel) {
  
  ## Find the order of environment means
  group_by(df, env) %>%
    summarise_at(.vars = trait, .funs = mean) %>%
    arrange_(trait) %>%
    pull(env) %>%
    as.character() ->
    env_order

  ## This assigns the env variables levels in the order of
  ## their mean values before sorting
  gxe$env <- factor(gxe$env, levels = env_order)
  gxe <- arrange(gxe, desc(env), geno)

  ## Plot regressions
  ggplot(gxe, aes_string(x = "env", y = trait)) +
    geom_smooth(aes(group = geno, color = geno), size = 1, method = "lm", se = FALSE) +
    xlab("Environment") +
    ylab(ylabel) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) ->
    plot_out
  
  return(plot_out)
}

NOTE that our use of theme() in line 22 above comes after our use of theme_bw() on the line above it. If theme_bw() comes after, it will override any custom theme options that we have set – you can tell I know this by hard-won experience.

Now to take care of the real problem – the multiple legends. The cowplot package has another really handy function – get_legend(). As the name implies, this function pulls the legend out of a figure generated by ggplot. So, in the code block below, we will:

Generate a dummy figure (just one of the plots that will be in the final figure)
Pull the legend out of it using get_legend()
Change the title of the legend from “geno” to “Genotype”
Assemble the final multi-panel figure using our two traits and our extracted legend
Label the two panels “A” and “B”

Note that we can also assign the portion of the overall plotting space that is taken up by each panel. So below, I am dividing the space up into fifths, and assigning each plot to 2/5 of the overall space, and the legend to 1/5.

dummy_plot <- plot_gen(gxe, "oil", "Oil Content") +
  scale_color_discrete(name = "Genotype")
legend <- get_legend(dummy_plot)

plot_grid(plot_gen(gxe, "oil", "Oil Content") + theme(legend.position = "none"),
          plot_gen(gxe, "height", "Plant Height (cm)") + theme(legend.position = "none"),
          legend,
          ncol = 3,
          rel_widths = c(2, 2, 1),
          labels = c("A", "B", ""))

Our final multi-panel figure:

Genotype-by-environment interaction plots for sunflower oil content (A) and plant height (B) measured across four environments

Discussion

That gets us pretty close to what we want. You’ll notice that not only are the units different for each plot, but the ordering of environments on the x axes are different as well. If there are additional tweaks that we want, we can always easily modify the plot_gen() function.

I won’t go into the details here, but as a final thought, I always suggest saving images in svg format. This allows us to both scale the image to any size without loss in fidelity, and to further tweak things by hand in a program like Adobe Illustrator or Inkscape. A great example of an easy svg tweak appears in our final image above – you’ll notice that the x-axis is labeled with “Environment” on both plots. It would be nice to have that appear just once, centered below both plots. This is an easy fix with an svg image – we can just manually remove one of the “Environment” labels, and drag the other copy to where we want it.

The graphs look pretty close to what I imagined – our messy trait (oil content) has some crossover genotype-by-environment interaction, while our more stable trait (plant height) exhibits nice, parallel lines in its interaction plot. I’m sure any plant breeder would prefer to work with the latter trait.

Now, as promised, my explanation as to why we did this in the first place. My justification is that in this case, we are more interested in the patterns of interaction for each of these two traits, rather than specific data points. To be more specific, it’s a matter of the scope of our inference. I could imagine generating this sort of plot if we were using a model that treated both genotype and environment as random effects. In that case, we would be more interested in characterizing some pool of possible genotypes and environments that we sampled in our experiment, rather than characterizing the specific genotypes or environments that we happened to use. (I use the term pool instead of population here, because I find the latter gets confusing when we start talking about populations in the genetic sense vs. the statistical sense).

With that, you now you have everything you need to go forth and compile all manner of multi-panel figures in R – please just commit to using this power wisely.