Thursday, April 28, 2022

Combining contMap and phylo.to.map plots part II: Downward facing trees

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

plot of chunk unnamed-chunk-6

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)

plot of chunk unnamed-chunk-8

Cool.

2 comments:

  1. This comment has been removed by the author.

    ReplyDelete
  2. Hey, awesome blog. I have taken inspiration from this for my PhD project.

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

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.