用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
对于 dotTree
和 phytools
中的类似函数(例如 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)
我正在尝试使用 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
对于 dotTree
和 phytools
中的类似函数(例如 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)