A cool (in my opinion) phytools plotting function that I've yet to see applied
in a published study is markChanges
, which can be used to visualize
the posterior density of changes (as opposed to trait values) from a
sample of stochastically mapped trees.
Here are a couple of different examples - one with lots of uncertainty about the location of trait changes on the tree (feeding mode in elopomorph eel evolution; from Collar et al. 2014); and another in which there should be much less (piscivory vs. non-piscivory in Centrarchidae).
First, the elopomorph dataset:
library(phytools)
eel.tree
##
## Phylogenetic tree with 61 tips and 60 internal nodes.
##
## Tip labels:
## Moringua_edwardsi, Kaupichthys_nuchalis, Gorgasia_taiwanensis, Heteroconger_hassi, Venefica_proboscidea, Anguilla_rostrata, ...
##
## Rooted; includes branch lengths.
feed.mode
## Albula_vulpes Anguilla_anguilla
## suction suction
## Anguilla_bicolor Anguilla_japonica
## suction suction
## Anguilla_rostrata Ariosoma_anago
## suction suction
## Ariosoma_balearicum Ariosoma_shiroanago
## suction suction
## Bathyuroconger_vicinus Brachysomophis_crocodilinus
## bite bite
## Conger_japonicus Conger_myriaster
## suction suction
## Conger_verreaxi Conger_wilsoni
## suction suction
## Congresox_talabonoides Cynoponticus_ferox
## suction bite
## Dysomma_anguillare Elops_saurus
## bite suction
## Facciolella_gilbertii Gavialiceps_taeniola
## bite bite
## Gnathophis_longicauda Gorgasia_taiwanensis
## suction suction
## Gymnothorax_castaneus Gymnothorax_flavimarginatus
## bite bite
## Gymnothorax_kidako Gymnothorax_moringa
## bite bite
## Gymnothorax_pseudothyrsoideus Gymnothorax_reticularis
## bite bite
## Heteroconger_hassi Ichthyapus_ophioneus
## suction bite
## Kaupichthys_hyoproroides Kaupichthys_nuchalis
## bite bite
## Megalops_cyprinoides Moringua_edwardsi
## suction bite
## Moringua_javanica Muraenesox_bagio
## bite bite
## Muraenesox_cinereus Myrichthys_breviceps
## bite suction
## Myrichthys_maculosus Myrichthys_magnificus
## suction suction
## Myrophis_vafer Nemichthys_scolopaceus
## bite bite
## Nettastoma_melanurum Ophichthus_serpentinus
## bite suction
## Ophichthus_zophochir Oxyconger_leptognathus
## suction bite
## Parabathymyrus_macrophthalmus Paraconger_notialis
## suction suction
## Pisodonophis_cancrivorus Poeciloconger_kapala
## bite suction
## Rhinomuraena_quaesita Rhynchoconger_flavus
## bite suction
## Saurenchelys_fierasfer Scolecenchelys_breviceps
## bite suction
## Scuticaria_tigrina Serrivomer_beanii
## bite bite
## Serrivomer_sector Simenchelys_parasitica
## bite suction
## Uroconger_lepturus Uropterygius_micropterus
## suction bite
## Venefica_proboscidea
## bite
## Levels: bite suction
mtrees<-make.simmap(eel.tree,feed.mode,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## bite suction
## bite -0.01582783 0.01582783
## suction 0.01582783 -0.01582783
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## bite suction
## 0.5 0.5
## Done.
cols<-setNames(c("red","blue"),c("bite","suction"))
dotTree(eel.tree,feed.mode,colors=cols,fsize=0.7,ftype="i",
legend=FALSE)
add.simmap.legend(x=0,y=-4,colors=cols,prompt=FALSE,
vertical=FALSE,shape="circle")
nulo<-sapply(mtrees,markChanges,sapply(cols,
make.transparent,0.1))
add.simmap.legend(colors=sapply(setNames(cols[2:1],
c("bite->suction","suction->bite")),
make.transparent,0.1),prompt=FALSE,x=50,y=-4,
vertical=FALSE)
Just to clarify the interpretation of this plot, the hash marks show the entire posterior sample of changes on the tree - across 100 samples from the posterior distribution. Any individual tree will have fewer (but still a lot of changes):
pd<-summary(mtrees)
pd
## 100 trees with a mapped discrete character with states:
## bite, suction
##
## trees have 30.95 changes between states on average
##
## changes are of the following types:
## bite,suction suction,bite
## x->y 17.58 13.37
##
## mean total time spent in each state is:
## bite suction total
## raw 1155.3150660 828.7861757 1984.101
## prop 0.5822863 0.4177137 1.000
Now here is the second dataset. This one piscivory vs. non-piscovory in Centrarchidae:
cent.tree
##
## Phylogenetic tree with 28 tips and 27 internal nodes.
##
## Tip labels:
## Acantharchus_pomotis, Lepomis_gibbosus, Lepomis_microlophus, Lepomis_punctatus, Lepomis_miniatus, Lepomis_auritus, ...
##
## Rooted; includes branch lengths.
piscivory
## Acantharchus_pomotis Lepomis_gibbosus Lepomis_microlophus
## pisc non non
## Lepomis_punctatus Lepomis_miniatus Lepomis_auritus
## non non non
## Lepomis_marginatus Lepomis_megalotis Lepomis_humilis
## non non non
## Lepomis_macrochirus Lepomis_gulosus Lepomis_symmetricus
## non pisc non
## Lepomis_cyanellus Micropterus_coosae Micropterus_notius
## pisc pisc pisc
## Micropterus_treculi Micropterus_salmoides Micropterus_floridanus
## pisc pisc pisc
## Micropterus_punctulatus Micropterus_dolomieu Centrarchus_macropterus
## pisc pisc non
## Enneacantus_obesus Pomoxis_annularis Pomoxis_nigromaculatus
## non pisc pisc
## Archolites_interruptus Ambloplites_ariommus Ambloplites_rupestris
## pisc pisc pisc
## Ambloplites_cavifrons
## pisc
## Levels: non pisc
mtrees<-make.simmap(cent.tree,piscivory,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## non pisc
## non -4.228856 4.228856
## pisc 4.228856 -4.228856
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## non pisc
## 0.5 0.5
## Done.
cols<-setNames(c("red","blue"),c("pisc","non"))
dotTree(cent.tree,piscivory,colors=cols,ftype="i")
nulo<-sapply(mtrees,markChanges,sapply(cols,
make.transparent,0.2))
add.simmap.legend(colors=sapply(setNames(cols[2:1],
c("non->pisc","pisc->non")),
make.transparent,0.2),prompt=FALSE,
x=par()$usr[1]+0.1*max(nodeHeights(cent.tree)),
y=-1.2,vertical=FALSE)
For comparison to our previous analysis:
summary(mtrees)
## 100 trees with a mapped discrete character with states:
## non, pisc
##
## trees have 6.95 changes between states on average
##
## changes are of the following types:
## non,pisc pisc,non
## x->y 3.27 3.68
##
## mean total time spent in each state is:
## non pisc total
## raw 0.6532209 1.041858 1.695079
## prop 0.3853630 0.614637 1.000000
That's it.
Quick question, I have phytools on R, but its tells me eel.tree is not present. How do you get it?
ReplyDelete