Search code examples
rtreephylogeny

Plotting traits on phylogeny with phytools


I am trying to plot trait data on a phylogeny using the phytools package. I'm sure this should be simple but I'm getting an unhelpful error message and I don't know what to try.

Here is my code including data download.



# General
library(dplyr)

# Phylogenetic libraries.
library(caper)

library(phytools)


#+ data_read
p <- read.table(file = 'http://esapubs.org/archive/ecol/E090/184/PanTHERIA_1-0_WR05_Aug2008.txt',
  header = TRUE, sep = "\t", na.strings = c("-999", "-999.00"))




## Some data cleaning

# Remove NAs in response and response where litter size is less than one (doesn't make sense).
p <- p %>% 
       filter(!is.na(X15.1_LitterSize)) %>% 
       filter(X15.1_LitterSize >= 1) %>% 
       mutate(y = log1p(X15.1_LitterSize)) %>% 
       dplyr::select(-X15.1_LitterSize, -References, -X24.1_TeatNumber)


## Get phylogeny data.

### read in phylogeny data.

# Read in trees
tree <- read.nexus('https://onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1111%2Fj.1461-0248.2009.01307.x&file=ELE_1307_sm_SA1.tre')

# Select best supported tree
tree <- tree[[1]]
tree$tip.label <- gsub('_', ' ', tree$tip.label)

# Check if species are available.
mean(p$MSW05_Binomial %in% tree$tip.label)
in_phylo <- p$MSW05_Binomial %in% tree$tip.label

# Remove data that is not in the phylogeny.

p <- p %>% filter(in_phylo) 

# Try just vulpes.

unneededTips <- tree$tip.label[!grepl('Vulpes', tree$tip.label) | !(tree$tip.label %in% p$MSW05_Binomial)]

# Prune tree down to only needed tips.
pruneTree <- drop.tip(tree, unneededTips)

dotTree(pruneTree, p$y[grepl('Vulpes', p$MSW05_Binomial)])




# Try all species

unneededTips <- tree$tip.label[!(tree$tip.label %in% p$MSW05_Binomial)]

# Prune tree down to only needed tips.
pruneTree <- drop.tip(tree, unneededTips)

dotTree(pruneTree, p$y)

I have tried plotting a smaller subset of the tree and the full tree but in both cases I get the error:

Error in if (k <= 0.8 && any(rr > (strwidth("W") * fsize/2))) rr <- rr/max(rr) * : missing value where TRUE/FALSE needed


Solution

  • for dotTree and similar functions in phytools (e.g. contMap), your trait value must be a named vector with the names corresponding to the tips in your tree. In your example you need to make sure p$y is a named vector (!is.null(names(p$y)) should be TRUE):

    ## Prune down the non Vulpes tips
    vulpes_tree <- drop.tip(tree, tree$tip.label[-grep("Vulpes", tree$tip.label)])
    
    ## Naming the variables in p$y
    all_vulpes <- grepl('Vulpes', p$MSW05_Binomial)
    traits_to_plot <- p$y[all_vulpes]
    names(traits_to_plot) <- p$MSW05_Binomial[all_vulpes]
    
    ## Plotting the Vulpes and the traits
    dotTree(vulpes_tree, traits_to_plot)
    

    You can apply the same procedure for your bigger tree. I suggest you use the function cleand.data from the dispRity package to match your tree and your dataset:

    ## Matching the tree and the data (using the dispRity package)
    library(dispRity)
    ## Attributing rownames to the dataset
    rownames(p) <- p$MSW05_Binomial
    ## Cleaning both the data and the tree
    cleaned_data <- dispRity::clean.data(p, tree)
    ## Extracting the cleaned dataset and the cleaned tree
    clean_p <- cleaned_data$data
    clean_tree <- cleaned_data$tree
    
    ## Same for the complete tree
    all_traits <- clean_p$y
    names(all_traits) <- clean_p$MSW05_Binomial
    
    ## Plotting all species and their traits
    dotTree(clean_tree, all_traits)