Solution to Challenge Problem 5: Stochastic Mapping

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).

  1. anguila.nex.
  2. elopomorph.csv

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:

obj<-summary(mapped.trees,check.equal=TRUE)
## Note: Some trees are not equal.
## A "reference" tree will be computed if none was provided.
## 
## No reference tree provided & some trees are unequal.
## Computing majority-rule consensus tree.
obj
## 100 trees with a mapped discrete character with states:
##  bite, suction 
## 
## trees have 30.08 changes between states on average
## 
## changes are of the following types:
##      bite,suction suction,bite
## x->y        16.37        13.71
## 
## mean total time spent in each state is:
##            bite    suction    total
## raw  14.9253103 11.2692163 26.19453
## prop  0.5710969  0.4289031  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")

plot of chunk unnamed-chunk-5

In addition to computing a simple consensus topology, we can also supply a consensus. This I will do by first computing a tree using a method called consensus.edges.

## compute consensus to supply as a reference tree
cons.tree<-consensus.edges(eel.trees,method="least.squares")
## [1] "RSS: 1.19631614789753"
## 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.08 changes between states on average
## 
## changes are of the following types:
##      bite,suction suction,bite
## x->y        16.37        13.71
## 
## mean total time spent in each state is:
##            bite    suction    total
## raw  14.9253103 11.2692163 26.19453
## prop  0.5710969  0.4289031  1.00000
## plot summary
plot(obj,ftype="i",fsize=0.7,colors=colors)
add.simmap.legend(x=0,y=5,prompt=FALSE,colors=colors,shape="circle")

plot of chunk unnamed-chunk-6

That's it.

Written by Liam J. Revell. Last updated 2 August 2017.