The challenge problem was as follows.
Use the following files containing (1) a set of trees from the Bayesian posterior distribution for eels & (2) phenotypic data for two characters for a subset of the species in (1).
Conduct stochastic mapping of feeding mode on this set of trees. Summarize the results using phytools summary methods, and plotting your findings.
First, let's load packages:
library(phytools)
Now, we can read the datasets from file:
eel.trees<-read.nexus("anguila.nex")
eel.data<-read.csv("elopomorph.csv",row.names=1)
This is not obligatory - in fact, we can conduct stochastic mapping in which we have incomplete data - however, for this exercise let's just prune the taxa from each of our trees that are absent from the data. Since our trees in eel.trees
come from a Bayesian analysis, we will just assume that they all contain the same set of taxa.
library(geiger)
chk<-name.check(eel.trees[[1]],eel.data)
chk
## $tree_not_data
## [1] "Ahlia_egmontis" "Anguilla_australis"
## [3] "Anguilla_bengalensis" "Anguilla_celebesensis"
## [5] "Anguilla_dieffenbachi" "Anguilla_interioris"
## [7] "Anguilla_malgumora" "Anguilla_marmorata"
## [9] "Anguilla_megastoma" "Anguilla_mossambica"
## [11] "Anguilla_nebulosa" "Anguilla_obscura"
## [13] "Anguilla_reinhardti" "Avocettina_infans"
## [15] "Bassanago_albescens" "Bassanago_bulbiceps"
## [17] "Bassanago_hirsutus" "Bathycongrus_retrotinctus"
## [19] "Chilorhinus_platyrhynchus" "Chilorhinus_suensonii"
## [21] "Coloconger_cadenati" "Derichthys_serpentinus"
## [23] "Diastobranchus_capensis" "E.pelecanoides"
## [25] "Echelus_myrus" "Echidna_xanthospilos"
## [27] "Elops_affinis" "Elops_sp"
## [29] "Eurypharynx_pelecanoides" "Gnathophis_nystromi"
## [31] "Gorgasia_sp" "Histiobranchus_bathybius"
## [33] "Hoplunnis_tenuis" "Ilyophis_brunneus"
## [35] "Macrocephenchelys_soela" "Meadia_abyssalis"
## [37] "Meadia_roseni" "Megalops_atlanticus"
## [39] "Muraenichthys_sp" "Myrophis_platyrhynchus"
## [41] "Myrophis_punctatus" "Nessorhamphus_danae"
## [43] "Nessorhamphus_ingolfianus" "Nettastoma_parviceps"
## [45] "Ophichthus_evermanni" "Ophichthus_gomesii"
## [47] "Ophichthus_puncticeps" "Ophisurus_macrorhynchos"
## [49] "Saccopharynx_ampullaceus" "Serrivomer_lanceolatoides"
## [51] "Serrivomer_samoensis" "Stemonidium_hypomelas"
## [53] "Synaphobranchus_affinis" "Synaphobranchus_kaupii"
## [55] "Thalassenchelys_coheni" "Uroconger_syringinus"
## [57] "Xenomystax_congroides" "Xyrias_revulsus"
##
## $data_not_tree
## character(0)
## prune trees:
eel.trees<-lapply(eel.trees,drop.tip,chk$tree_not_data)
class(eel.trees)<-"multiPhylo"
Now, let's do stochastic mapping:
x<-setNames(eel.data[,1],rownames(eel.data))
x
## 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
mapped.trees<-make.simmap(eel.trees,x,model="ARD",message=FALSE)
mapped.trees
## 100 phylogenetic trees with mapped discrete characters
This has just generated one stochastic map per tree - but normally we would use more (10 or 100). In this case, doing so would result in a computationally burdensome analysis time.
Now, we are going to use the S3 summary
method for objects
of class "multiSimmap"
. When trees are incongruent, the method
will compute a consensus tree; however instead we are going to supply it
with one computed using the phytools method consensus.edges
.
## compute consensus to supply as a reference tree
cons.tree<-consensus.edges(eel.trees,method="least.squares")
## perform summary
obj<-summary(mapped.trees,ref.tree=cons.tree)
obj
## 100 trees with a mapped discrete character with states:
## bite, suction
##
## trees have 30.84 changes between states on average
##
## changes are of the following types:
## bite,suction suction,bite
## x->y 17.18 13.66
##
## mean total time spent in each state is:
## bite suction total
## raw 15.0791819 11.1153447 26.19453
## prop 0.5781189 0.4218811 1.00000
## plot summary
colors<-setNames(c("red","blue"),c("bite","suction"))
plot(obj,ftype="i",fsize=0.7,colors=colors)
add.simmap.legend(x=0,y=5,prompt=FALSE,colors=colors,shape="circle")
That's it.
Written by Liam J. Revell. Last updated 12 Mar. 2016.