Я пытаюсь построить данные о признаках филогении, используя пакет 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
является вектором названный (!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)
Вы можете применить ту же процедуру для вашего большего дерева.
Я предлагаю вам использовать функцию cleand.data
из пакета dispRity
, чтобы сопоставить ваше дерево и ваш набор данных:
## 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)
Аааа спасибо большое. Я проголосую и приму, как только у меня будет время запустить его на своем компьютере, чтобы убедиться, что он работает.