CODEWITHSUNDEEP
rtreebioinformaticsphylogenyape
j91
j91

Reputation: 501

Color the tips of a phylogenetic tree by individual's species

I'm trying to build a phylogenetic tree with bootstrapping using the methodology defined here. Plotting the name of each individual (named 1,2,3,4,... for the sake of this examples) is not informative since I have a HUGE data set. Therefore, I want the tips to be coloured to reflect the individual's species. The individual species information is stored in a separated matrix called species:

species = data.frame(Ind=c(1 ,2  ,3  ,4  ,5 ),
                     Spe=c("s1","se1","se2","se2","se3"))

Nevertheless, whenever I tried to replace the names of the individuals with its species in fit$tree$label, the only results is empting the vector. Probably this is because fit is an object of class pml, which my generate conflict with my strategy, but I'm not sure. Is there a way to colour the tips of the tree based on the individual species? My code looks like this, for now:

library(phangorn)
#Create a tree from data in fasta format
dat = read.phyDat(file = "myalignment.fasta", format ="fasta")
tree <- pratchet(dat)          # parsimony tree
mt <- modelTest(dat, tree=tree, multicore=TRUE)
mt[order(mt$AICc),]
bestmodel <- mt$Model[which.min(mt$AICc)]
env = attr(mt, "env")
fitStart = eval(get("GTR+G+I", env), env)
fit = optim.pml(fitStart, rearrangement = "stochastic", optGamma=TRUE, optInv=TRUE, model="GTR")
bs = bootstrap.pml(fit, bs=100, optNni=TRUE, multicore=TRUE)

#Replace the names with the species...
fit$tree$tip.label <- species[which(species[,1] == fit$tree$tip.label),2]
#If I print fit$tree$tip.label here, the output is factor(0)

#...and create the tree with colored tips
plotBS(midpoint(fit$tree), bs, p = 50, type="p", show.tip.label = FALSE)
tiplabels(pch=19, col = as.factor(fit$tree$tip.label), adj = 2.5, cex = 2)

REPRODUCIBLE EXAMPLE

Save the following with the name "myalignment.fasta" and run the code above. It should create a toy example to play with:

>1
AACCAGGAGAAAATTAA
>2
AAAAA---GAAAATTAA
>3
ACACAGGAGAAAATTAA
>4
AACCTTGAGAAAATTAT
>5
CCTGAGGAGAAAATTAA

Upvotes: 0

Views: 1939

Answers (1)

MrFlick
MrFlick

Reputation: 206253

So if your problem is list plotting points as colors rather than labels, you can do

tiplabels(pch=19, cex=2, 
    col = species$Spe[match(fit$tree$tip.label, species$Ind)])

I do a bit of gynmastics just to make sure the labels match up to the species (though they were in the same order in this example).

enter image description here

Upvotes: 2

Related Questions