I am trying to specify the colours used in my ancestral state reconstruction analysis so that I can then use the same colours in other graphs the same data. However, when I used nodelabels
from ape the colours I specify do not line up with what is displayed. Does anyone know why this is? And how to work around it? See example of the problem below:
rm(list = ls())
library(phytools)
library(ape)
library(RColorBrewer)
set.seed(1237)
tree = rtree(50)
plot(tree)
variable = rTraitDisc(tree, k = 5)
names(variable) = tree$tip.label
cols = data.frame(type = levels(variable),
color = I(brewer.pal(nlevels(variable), name = 'Set1')))
cols_vector = cols$color
names(cols_vector) = cols$type
fit = ace(x = variable, phy = tree, type = 'discrete')
nodelabels(node=1:tree$Nnode+Ntip(tree),
pie=fit$lik.anc, piecol=cols_vector)
tiplabels(pie=to.matrix(variable,sort(unique(variable))),piecol=cols_vector)
tiplabels(text = variable)
Which gives the image:
Most of the colours line up ok, but E is displayed as purple but has the code #FF7F00 (which is orange). Which is exemplified by this code:
t = table(variable)
barplot(t, col = cols_vector[match(names(t), names(cols_vector))])
Any advice would be appreciated.
The problem is in command to.matrix(variable,sort(unique(variable)))
. Given your seed, no tips have value D and thus if you construct a matrix with columns representing unique values, only the first four colours from your palette will be used.
head(to.matrix(variable,sort(unique(variable))))
A B C E
t47 1 0 0 0
t14 1 0 0 0
t2 1 0 0 0
t19 1 0 0 0
t29 1 0 0 0
t5 1 0 0 0
Replacing the seq =
argument with the vector of all your possible trait values should solve the problem with non-matching colours.
tiplabels(pie = to.matrix(x = variable, seq = LETTERS[1:5]), piecol = cols_vector)