R/03_data-processing-phylogeny.R

In phytomosaic/pnw: Tree functional strategies in the US Pacific Northwest

######################################################################
#
#  Get phylogeny for PNW trees
#
#  Rob Smith, phytomosaic@gmail.com, 02 April 2020
#
##  CC-BY-SA 4.0 License (Creative Commons Attribution-ShareAlike 4.0)


### preamble
rm(list=ls())
require(ecole)
require(viridis)
require(V.PhyloMaker)
require(brranching)
load('./data/pnw_tree.rda', verbose=T)
d <- pnw_tree  ;  rm(pnw_tree)

### get family / genus / species
(tax <- sort(unique(d$spp)))
tax  <- brranching::phylomatic_names(tax, db='apg')
tax[gsub('\\/.*','',tax) == 'NA']  # expect 0
tax  <- data.frame(do.call(rbind, strsplit(tax,'/')))
tax  <- data.frame(species = tax[,2],
                   genus = sub('\\_.*', '', tax[,2]),
                   family = tax[,1])
anyNA(tax)
tax

### query V.PhyloMaker ! ! ! TIMEWARN ! ! !  57 taxa ~30 sec
system.time(p <-V.PhyloMaker::phylo.maker(sp.list=tax,scenarios='S3'))
phy <- p$scenario.3
phy$tip.label <- tolower(phy$tip.label)
table(p$species.list$status)
if (all(tax$species %in% phy$tip.label) &&
    all(phy$tip.label %in% tax$species) ) {
     cat('all taxa appear in `phy`\n')
     rm(tax, p)
}
pnw_phy <- phy
save(pnw_phy, file='./data/pnw_phy.rda')
rm(pnw_phy)

### plot phylogeny, colored by frequency
frq  <- aggregate(d$spp, list(d$spp), length)
(frq <- frq[match(phy$tip.label, frq[,1]),]) # order to match phy
cc <- ecole::colvec(log10(frq$x), pal=viridis::viridis(99), alpha=1)
png('./fig/fig_00_phylogeny.png', wid=6.5, hei=4.5, unit='in',
    bg='transparent', res=900)
plot(ape::ladderize(phy), no.margin=TRUE, cex=0.7, tip.col=cc,
     label.offset=10, edge.width = 2)
tiplabels(pch=19, col=cc, cex=0.8)
dev.off()

# ### Ancestral trait reconstruction
# ### continuous trait reconstruction _across branches_
# trm <- contMap(phy, x=tr, lims=c(16.5,33.5), plot=FALSE)
# trm$cols[1:length(trm$cols)] <- viridis::inferno(length(trm$cols))
# ### continuous trait reconstruction _at each node_
# fa <- fastAnc(phy, tr)
# plot(trm, fsize=0.4, lwd=3, outline=F, leg.txt='Trait values')
# nodelabels(pch=19, col = colvec(fa), cex = 0.9)
# tiplabels(pch=19, col = colvec(tr,alpha=1), cex = 0.8)

####    END    ####