A couple of days ago I
posted
an example showing how to overplot a "contMap" graph on a "phylo.to.map"
projection of a phylogeny to a geographic map.
The example I showed, however, was for a right-facing tree plotted to the left of a vertical (portrait) geographic map.
This is also possible to do for a downward oriented tree plotted above a landscape orientation geographic map.
Here my phylogenetic tree, lat/long information, and trait data are in
World.tree, World.latlong, and World.trait, respectively, but I'm imagining
that I have multiple geographic observations per species, but only one phenotypic trait
value. (Simulation code is
here.)
## load packages
library(phytools)
packageVersion("phytools") ## remember, you need this version!
## [1] '1.0.4'
Data:
World.tree
##
## Phylogenetic tree with 26 tips and 25 internal nodes.
##
## Tip labels:
## B.mjdns, U.ira, P.igtzpylhm, K.ubhoryjdva, R.gozlnst, Y.btopuyzmq, ...
##
## Rooted; includes branch lengths.
World.latlong
## lat long
## B.mjdns 16.919556 113.863811
## U.ira 12.147042 117.430738
## P.igtzpylhm -45.751652 -75.329094
## K.ubhoryjdva -45.998366 -28.283127
## R.gozlnst 47.641988 -42.318668
## Y.btopuyzmq 16.757773 -41.207613
## X.zytplwfvum 37.170279 -36.422596
## G.hzfjatxd 26.896980 -16.049381
## D.sczwfvx 14.566583 6.157820
## Z.xvlswt 28.904822 33.803170
## Q.voyj 50.286729 56.479024
## L.pibqh 13.173950 -69.101749
## G.hbkwlzex 65.249831 19.752462
## G.qxfragjpdl 60.363030 -83.871754
## Q.vhxg 52.227775 -30.332072
## K.hscaoe 54.415468 -83.886786
## F.gkdz 41.718283 -45.708856
## O.tvpuxgo 50.265279 23.357039
## Q.ifmlj -29.188547 -48.539389
## X.dbcl 39.347408 17.642643
## Q.lzqxyasg 40.843105 -12.175260
## C.xiyogqzmk 42.591035 -20.614369
## Z.phgimc -40.716567 104.626179
## N.pfwn 7.269377 43.531950
## N.nigtewk 9.608382 74.669080
## Y.yiomvf -15.202712 15.029529
## B.mjdns 10.480000 106.958868
## B.mjdns 17.037866 111.604051
## U.ira 14.592380 114.101974
## U.ira 13.460313 120.932577
## P.igtzpylhm -51.634712 -77.239777
## P.igtzpylhm -49.239105 -69.240627
## K.ubhoryjdva -51.357916 -31.586794
## K.ubhoryjdva -38.005277 -32.274457
## R.gozlnst 47.933340 -41.674420
## R.gozlnst 45.228145 -46.494576
## Y.btopuyzmq 25.343239 -39.459538
## X.zytplwfvum 41.796516 -32.244525
## G.hzfjatxd 33.081089 -16.118202
## G.hzfjatxd 34.067357 -11.421670
## D.sczwfvx 13.763942 7.515051
## D.sczwfvx 14.128540 6.998809
## Z.xvlswt 24.014239 32.823913
## Z.xvlswt 25.507820 35.315387
## Q.voyj 47.178446 61.363431
## L.pibqh 20.456394 -75.428737
## L.pibqh 9.072784 -68.180213
## G.hbkwlzex 64.992838 15.199572
## G.hbkwlzex 62.840271 19.476022
## G.qxfragjpdl 62.514989 -87.234281
## G.qxfragjpdl 52.789052 -86.182006
## Q.vhxg 49.949079 -33.518558
## K.hscaoe 57.536418 -84.353937
## F.gkdz 38.713806 -48.498001
## F.gkdz 40.125053 -46.138437
## O.tvpuxgo 55.737767 23.339125
## Q.ifmlj -28.456218 -46.215204
## Q.ifmlj -32.896343 -46.065622
## X.dbcl 31.452855 12.216477
## Q.lzqxyasg 41.807206 -17.495484
## C.xiyogqzmk 42.976644 -23.147773
## C.xiyogqzmk 38.301593 -13.177576
## Z.phgimc -47.393028 94.014434
## Z.phgimc -38.529515 102.938858
## N.pfwn 13.275540 47.454354
## N.nigtewk 3.238530 71.269795
## Y.yiomvf -14.415705 11.733332
## Y.yiomvf -20.686126 9.026657
World.trait
## B.mjdns U.ira P.igtzpylhm K.ubhoryjdva R.gozlnst Y.btopuyzmq
## -11.157621 -12.184928 -10.300926 -9.098500 -1.911623 -1.796410
## X.zytplwfvum G.hzfjatxd D.sczwfvx Z.xvlswt Q.voyj L.pibqh
## -10.591982 -2.242069 -13.570787 -11.901156 -5.563962 -11.873892
## G.hbkwlzex G.qxfragjpdl Q.vhxg K.hscaoe F.gkdz O.tvpuxgo
## -11.291176 -12.162927 -4.779088 -2.447187 -1.363460 -8.349076
## Q.ifmlj X.dbcl Q.lzqxyasg C.xiyogqzmk Z.phgimc N.pfwn
## -12.051234 -17.144910 -18.559916 -18.251086 -3.390340 -10.642816
## N.nigtewk Y.yiomvf
## -10.049277 -5.857286
Next, compute our "phylo.to.map" object:
obj<-phylo.to.map(World.tree,World.latlong,plot=FALSE)
## objective: 202
## objective: 182
## objective: 182
## objective: 182
## objective: 180
## objective: 180
## objective: 158
## objective: 158
## objective: 156
## objective: 152
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
## objective: 150
obj
## Object of class "phylo.to.map" containing:
##
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
##
## (2) A geographic map with range:
## -85.19N, 83.6N
## -180W, 180W.
##
## (3) A table containing 68 geographic coordinates (may include
## more than one set per species).
##
## If optimized, tree nodes have been rotated to maximize alignment
## with the map when the tree is plotted in a downwards direction.
Set colors for plotting:
cols<-setNames(sample(rainbow(n=Ntip(World.tree))),
World.tree$tip.label)
cols
## B.mjdns U.ira P.igtzpylhm K.ubhoryjdva R.gozlnst Y.btopuyzmq
## "#FF00B1" "#00FFC4" "#FFB100" "#00FFFF" "#0014FF" "#00C4FF"
## X.zytplwfvum G.hzfjatxd D.sczwfvx Z.xvlswt Q.voyj L.pibqh
## "#0089FF" "#FF7600" "#FF0000" "#FF0076" "#9D00FF" "#FFEB00"
## G.hbkwlzex G.qxfragjpdl Q.vhxg K.hscaoe F.gkdz O.tvpuxgo
## "#004EFF" "#2700FF" "#FF003B" "#62FF00" "#00FF4E" "#27FF00"
## Q.ifmlj X.dbcl Q.lzqxyasg C.xiyogqzmk Z.phgimc N.pfwn
## "#D8FF00" "#6200FF" "#D800FF" "#00FF89" "#FF3B00" "#FF00EB"
## N.nigtewk Y.yiomvf
## "#9DFF00" "#00FF14"
We can even plot a test version without our "contMap" phylogeny:
plot(obj,colors=cols,ftype="i",fsize=0.6,cex.points=c(0.7,1.2))
Now compute our "contMap" and update the color palette:
cMap<-contMap(obj$tree,World.trait,plot=FALSE)
cMap
## Object of class "contMap" containing:
##
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
##
## (2) A mapped continuous trait on the range (-18.559916, -1.36346).
cMap<-setMap(cMap,c("lightblue","blue","purple","red","darkred"))
Finally, we're reading for our full graph. The one big difference from before is that
to make things easier, after plotting our geographic map & base tree, we'll flip the
vertical (y) axis before adding our "contMap" tree to make things a tiny bit easier.
## plot our phylogeny & map
plot(obj,colors=cols,ftype="i",fsize=0.6,cex.points=c(0.7,1.2),
pts=FALSE)
## get "last_plot.phylo" object
lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## rescale our "contMap" tree
tt<-cMap$tree
tt<-rescaleSimmap(tt,depth=diff(range(lastPP$yy)))
## graph the contMap tree
## in this case we'll flip
plot(tt,cMap$cols,add=TRUE,lwd=4,ftype="off",
ylim=lastPP$y.lim[2]-lastPP$y.lim,
xlim=lastPP$x.lim,
tips=lastPP$xx[1:Ntip(cMap$tree)],
mar=rep(0.01,4),asp=1,direction="upwards")
## add a color gradient legend
add.color.bar(70,cMap$cols,"phenotype",lims=cMap$lims,
digits=1,prompt=FALSE,lwd=6,outline=TRUE,x=-170,y=255,
subtitle="",fsize=0.8)
Cool.
This comment has been removed by the author.
ReplyDeleteHey, awesome blog. I have taken inspiration from this for my PhD project.
ReplyDeleteQuick question. How would I go about retaining the pts? as you remove them:
plot(obj,colors=cols,ftype="i",fsize=0.6,cex.points=c(0.7,1.2),
pts=FALSE)
Say I wanted to keep them in the final plot? In their current state, if pts=T then it becomes covered by the new plot's overlay.
Thanks again.