Follow-up on plotting facing "arc" or "fan" style trees

Yesterday, I posted about graphing two “arc” or “fan” style phylogenies in a facing (up/down) style.

In my example, the trees showed alternative ancestral state reconstructions of the same trait, in which the bottom tree was a rotated version of the upper. A phytools user pointed out, however, that it might be easier to read if the trees mirrored one another. I agree.

I like it but wonder if it would be easier to interpret if these were mirrored. Just a thought. https://t.co/omIjQ1QNCX

— A.N. Atomist (@Anatonomicon) December 10, 2023

Very fortunately, this is not too difficult to accomplish.

One option would be to use the code I posted yesterday, but wherein we simply rotated all the nodes of one tree or the other (using, for instance, phytools::rotateNodes).

Another option also occurred to me, however, which is to use the phytools tree plotting argument tips, which sets the tip order.

The following is an example using spawning mode (group vs. pair) and paternal care (present or absent) in bonyfishes from Benun Sutton and Wilson (2019).

## load phytools
library(phytools)

## load tree and data
data("bonyfish.tree")
bonyfish.tree

## 
## Phylogenetic tree with 90 tips and 89 internal nodes.
## 
## Tip labels:
##   Xenomystus_nigri, Chirocentrus_dorab, Talismania_bifurcata, Alepocephalus_tenebrosus, Misgurnus_bipartitus, Opsariichthys_bidens, ...
## 
## Rooted; includes branch lengths.

data("bonyfish.data")
head(bonyfish.data)

##                          spawning_mode paternal_care
## Xenomystus_nigri                  pair          male
## Chirocentrus_dorab               group          none
## Talismania_bifurcata             group          none
## Alepocephalus_tenebrosus         group          none
## Misgurnus_bipartitus              pair          none
## Opsariichthys_bidens              pair          none

spawning_mode<-bonyfish.data$spawning_mode
paternal_care<-bonyfish.data$paternal_care
names(spawning_mode)<-names(paternal_care)<-
  rownames(bonyfish.data)

## ancestral state reconsrunction of
## spawning mode
spawning_mode.anc<-ancr(fitMk(bonyfish.tree,
  spawning_mode,pi="fitzjohn"))
spawning_mode.anc

## Marginal ancestral state estimates:
##       group     pair
## 91 0.059480 0.940520
## 92 0.074462 0.925538
## 93 0.098419 0.901581
## 94 0.097894 0.902106
## 95 0.970123 0.029877
## 96 0.008840 0.991160
## ...
## 
## Log-likelihood = -29.51819

## ancestral state reconstruction of
## paternal care
paternal_care.anc<-ancr(fitMk(bonyfish.tree,
  paternal_care,pi="fitzjohn"))
paternal_care.anc

## Marginal ancestral state estimates:
##        male     none
## 91 0.011804 0.988196
## 92 0.002168 0.997832
## 93 0.001939 0.998061
## 94 0.002191 0.997809
## 95 0.003462 0.996538
## 96 0.004425 0.995575
## ...
## 
## Log-likelihood = -33.150827

par(mfrow=c(2,1))
cols<-setNames(viridisLite::viridis(n=2),
  levels(spawning_mode))
cex<-apply(spawning_mode.anc$ace,1,
  function(x) if(any(x>0.95)) 0.3 else 0.6)
plot(spawning_mode.anc,
  args.plotTree=list(type="arc",arc_height=0.5,
    fsize=0.5,mar=c(0.1,0.1,2.1,0.1)),
  args.nodelabels=list(piecol=cols,cex=cex),
  args.tiplabels=list(cex=0.2),
  legend=FALSE)
mtext("a) ancestral reconstruction of spawning mode",
  line=0,adj=0.05)
legend("topright",c("group","pair"),title="spawning mode",
  pch=16,col=cols,cex=1,pt.cex=2,bty="n")
xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
ylim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim
cols<-setNames(viridisLite::rocket(n=2),
  levels(paternal_care))
cex<-apply(paternal_care.anc$ace,1,
  function(x) if(any(x>0.95)) 0.2 else 0.5)
plot(paternal_care.anc,
  args.plotTree=list(type="arc",arc_height=0.5,
    xlim=xlim,ylim=-ylim[2:1],
    tips=Ntip(bonyfish.tree):1+Ntip(bonyfish.tree),
    part=1,fsize=0.5,mar=c(0.1,0.1,2.1,0.1)),
  args.nodelabels=list(piecol=cols,cex=cex),
  args.tiplabels=list(cex=0.2),
  legend=FALSE)
mtext("b) ancestral reconstruction of paternal care",
  line=0,adj=0.05)
legend("bottomright",c("present","absent"),
  title="paternal care",
  pch=16,col=cols,cex=1,pt.cex=2,bty="n")

It worked!