用phytools绘制系统发育特征

Plotting traits on phylogeny with phytools

我正在尝试使用 phytools 包绘制系统发育的性状数据。我确信这应该很简单,但我收到了一条无用的错误消息,我不知道该尝试什么。

这是我的代码,包括数据下载。



# 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)

我尝试绘制树的较小子集和完整树,但在这两种情况下我都收到错误:

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

对于 dotTreephytools 中的类似函数(例如 contMap),您的特征值必须是一个命名向量,其名称与树中的提示相对应。 在您的示例中,您需要确保 p$y 是一个 named 向量(!is.null(names(p$y)) 应该是 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)

您可以对更大的树应用相同的程序。 我建议您使用 dispRity 包中的函数 cleand.data 来匹配您的树和数据集:

## 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)