Neotoma Part II: extinct mammal distributions

library(ggplot2)
library(tidyverse)
library(skimr)
library(taxize)

## Read in previous data.
all_fauna_data <- read.csv("../data/all_fauna_data.csv")

## subset only mammals
mammals <- all_fauna_data %>%
  filter(taxon.group == "Mammals")

skim(all_fauna_data) 

# I have a love / hate relationship with this first species
humans <- mammals %>% 
  filter(taxon.name == "Homo sapiens" & long > -160) #ignoring samples out in the middle of the ocean

# Start with world
world <- map_data("world") 

ggplot() + 
  geom_polygon(data = world, 
               fill = "grey38", 
               aes(x = long, y = lat, group = group)) + 
  coord_fixed(1.3) +
  geom_point(data = humans, 
             aes(x = long, 
                 y = lat,
                 color = age.older),
             alpha = .7, size = .7) +
  scale_color_continuous(guide = guide_legend(title = "Age (years ago)"))  +
    guides(colour = guide_legend(override.aes = list(size = 3))) +
  ggtitle("Geographic Distribution of Homo sapien Samples")

usa <- map_data("state") 

ggplot() + 
  geom_polygon(data = usa, 
               fill = "grey38", 
               aes(x = long, y = lat, group = group)) + 
  coord_fixed(1.3) +
  geom_point(data = humans, 
             aes(x = long, 
                 y = lat,
                 color = age.older),
             alpha = .7, size = .4) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
   ggtitle("Geographic Distribution of Homo sapien Samples") 

skim(mammals$age.older)

skim(mammals$age.younger)

smilodon <- mammals[grepl("Smilodon ", mammals$taxon.name),]

unique(smilodon$taxon.name) # looks like 5 unique categories.

usa <- map_data("state") 

ggplot() + 
  geom_polygon(data = usa, fill = "grey30", 
               aes(x = long, y = lat, group = group)) + 
  coord_fixed(1.3) +
  geom_point(data = smilodon, 
             aes(x = long, y = lat, color = taxon.name),
             alpha = .9, size = .7) +
    scale_color_brewer(palette = "Pastel1") +
    theme(legend.key = element_rect(fill = "white")) +
    guides(colour = guide_legend(override.aes = list(size = 3))) +
    ggtitle("Geographic Distribution of Smilodon Samples")

## Mammuthus is 
mammoth_test <- mammals[grepl("Mammuthus", mammals$taxon.name),]
unique(mammoth_test$taxon.name) # 16 unique Mammoth taxon names

## Take a min to run
classification("Mammuthus", db = 'tol')

## Does not like how my taxa.names are written
whole_mammoth_test <- classification(unique(mammoth_test$taxon.name), db = 'tol') 

## Show which species did not work
as.data.frame(is.na(whole_mammoth_test))

## From Simon 
## "This doesn't catch things like "sensu stricto" and others."

mammoth_test$taxon.name <- gsub("(\\?|\\-type|cf\\.\\s|aff\\.|\\sundiff\\.)", "", mammoth_test$taxon.name, perl = TRUE)

table(mammoth_test$taxon.name) # use table() to both view the new and improved Mammoth taxon.names, and to see how many samples fall under each name.

clean_mammoth_test <- classification(unique(mammoth_test$taxon.name), db = 'tol') 

## Yay! Much better. No NAs.
as.data.frame(is.na(clean_mammoth_test))

get_phylogeny <- function(phylo_list) {
  #phylo list is the tree of life mammoth info data
  #taxon.name <- names(phylo_list) # the current taxa

  order <- phylo_list %>% 
              filter(rank == "order") %>% 
              {.$name}
  family <- phylo_list %>% 
              filter(rank == "family") %>%
              {.$name}
  genus <- phylo_list %>% 
              filter(rank == "genus") %>%
              {.$name}
  tmp_list <- list("order" = order,"family"=family,"genus"=genus)
  
  x <- phylo_list %>% 
    filter(rank %in% c("order", "family", "genus")) %>%
    select(-id) %>%
    spread(rank, name)
  
  if(dim(x)[2] != 3) { #probably need a better check here.
    missing <- setdiff(c("order", "family", "genus"), x$rank)
    add <- data.frame(name = rep(NA, length(missing)), rank=missing)
    x <- rbind(x, add)
   }
  tmp_list <- list("order"= x$order, "family"= x$family,"genus"= x$genus)
  return(tmp_list)
}

tol_classification_list <- sapply(clean_mammoth_test, get_phylogeny)
tol_classif <- as.data.frame(t(tol_classification_list))
tol_classif$taxon.name <- rownames(tol_classif)
rownames(tol_classif) <- NULL

merged <- merge(mammoth_test, tol_classif, by = "taxon.name")

## Lets take a peak of the dataset
merged %>%
  select(taxon.name, order, family, genus, lat, long, ) %>%
  head()

neotoma_taxa <- neotoma::get_table("taxa")
extinct_taxa <- neotoma_taxa %>% 
  filter(Extinct == "TRUE")

## Find the intersect of extinct taxa and our mammal dataset
extinct_species_in_mammals <- intersect(extinct_taxa$TaxonName, mammals$taxon.name)

## Great! 383 of these extinct species are found in our
## mammal dataset.
skim(extinct_species_in_mammals)

extinct_mammals <- mammals %>% 
                      filter(taxon.name %in% extinct_species_in_mammals)

## Check
## skim(extinct_mammals) 
nrow(extinct_mammals)

extinct_mammals$taxon.name <- gsub("(\\?|\\-type|cf\\.\\s|aff\\.|\\sundiff\\.)", "", extinct_mammals$taxon.name, perl = TRUE)

# this species throws an error later
# during classification step, remove
extinct_mammals <- extinct_mammals %>% 
  filter(taxon.name != 'Ammospermophilus leucurus') 

## remove white space at the beginning taxon names, 
## if it exists, using the base function trimws()
extinct_mammals$taxon.name <- trimws(extinct_mammals$taxon.name)

## Warning: This next line takes awhile to run. 
## and also requires manual interaction
## you can download data set below
## Manual selection key: 1, 0, 2, 0, 2, 1 
extinct_tol <- classification(unique(extinct_mammals$taxon.name), db = 'tol') 

## You end up with a list of dataframes. 
## One for each taxon.name input

length(names(extinct_tol[!is.na(extinct_tol)])) # 286 were classified
length(names((extinct_tol[is.na(extinct_tol)]))) #27 were not classified

## Keep only dataframes that are not NA
extinct_tol <- extinct_tol[!is.na(extinct_tol)]

extinct_tol_list <- list()
for (i in 1:length(names(extinct_tol))) {
  taxon.name <- names(extinct_tol)[i]

  x <- extinct_tol[[i]] %>% 
    filter(rank %in% c("order", "family", "genus")) %>%
    select(-id)
  
  if(dim(x)[1] != 3) { #probably need a better check here.
    missing <- setdiff(c("order", "family", "genus"), x$rank)
    add <- data.frame(name = rep(NA, length(missing)), rank=missing)
    x <- rbind(x, add)
  }
  
  y <- c(taxon.name, as.character(x$name))
  names(y) <- c("taxon.name", "order", "family", "genus")
  extinct_tol_list[[i]] <- y

}


tol_classification <- as.data.frame(do.call(rbind, extinct_tol_list))
tol_classification[,5] <- NULL #whoops repeated column

extinct_merged <- merge(extinct_mammals, tol_classification, by = "taxon.name", all.x = TRUE)
dim(extinct_merged)

## Write out
## write.csv(extinct_merged, "../data/part2_extinct_merged_27March2018_caryn.csv", row.names = FALSE)

## Read in data frame
extinct_merged <- read.csv("../data/part2_extinct_merged_27March2018.csv")

extinct_merged %>%
  group_by(order) %>%
  count() %>%
    ggplot(., aes(reorder(x = order, -n), y = n)) + 
      geom_bar(stat = "identity") +
      theme_bw() +
      theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
      xlab("Extinct Mammals") +
      ylab("# of samples") +
      ggtitle("Number of Extinct Mammal Samples in Each Animal Order")

orders_to_keep <- c("Rodentia", "Proboscidea", "Pilosa", "Perissodactyla", "Lagomorpha", "Didelphimorphia", "Cingulata", "Chiroptera", "Cetartiodactyla", "Carnivora", "Anguilliformes", "Primates")	

extinct_merged %>% 	
  filter(order %in% orders_to_keep) %>%	
  group_by(order) %>%	
  count %>%	
  ggplot(., aes(reorder(x = order, -n), y = n)) + 
    geom_bar(stat = "identity") +
    theme_bw() +	
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    xlab("Extinct Mammals") +
    ylab("# of samples") +
  ggtitle("Number of Extinct Mammal Samples in Each Animal Order")

usa <- map_data("state") 

extinct_merged %>%
  filter(order %in% orders_to_keep) %>%
  filter(lat < 50) %>%
    ggplot(.) + 
      geom_polygon(data = usa, 
              fill = "grey30", 
              aes(x = long, y = lat, group = group)) +
      coord_fixed(1.3) +
      geom_point( aes(x = long, y = lat, color = order),  
             alpha = .7, size = .5) +
  scale_color_brewer(palette = "Set3") +
  theme(legend.key = element_rect(fill = "white")) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  ggtitle("Geographic Distribution of Extinct Mammal Samples")

extinct_merged %>%
  filter(order %in% orders_to_keep) %>%
  filter(order == "Proboscidea") %>%
  filter(lat < 50) %>%
    ggplot(.) + 
      geom_polygon(data = usa, 
              fill = "grey30", 
              aes(x = long, y = lat, group = group)) +
      coord_fixed(1.3) +
      geom_point( aes(x = long, y = lat, color = genus ),  
             alpha = .7, size = .6) +
    scale_color_brewer(palette = "Pastel1") +
    theme(legend.key = element_rect(fill = "white")) +
    guides(colour = guide_legend(override.aes = list(size = 3))) +
    facet_grid(family~.) +
    ggtitle("Geographic Distribution of Proboscidea Samples")

extinct_merged %>%
  filter(order %in% orders_to_keep) %>%
  filter(order == "Proboscidea") %>%
  skim(age.older)

extinct_merged %>%
  filter(order %in% orders_to_keep) %>%
  filter(family == "Mammutidae") %>%
  filter(lat < 50) %>%
    ggplot(.) + 
      geom_polygon(data = usa, 
              fill = "grey30", 
              aes(x = long, y = lat, group = group)) +
      coord_fixed(1.3) +
      geom_point( aes(x = long, y = lat, color = taxon.name ),  
             alpha = .7, size = .7) +
    scale_color_brewer(palette="Pastel1") +
    theme(legend.key = element_rect(fill = "white")) +
    guides(colour = guide_legend(override.aes = list(size = 3))) +
    ggtitle("Geographic Distribution of Mammutidae Samples")