Extracting posterior ages of nodes from beast outputs, to graph distributions

[Jo Wolfe]

I would like to create density graphs of the posterior/marginal prior resulting from divergence time analysis. Similar to Fig S13 in my previous paper: https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frspb.2019.0079&file=rspb20190079supp1.pdf
but the previous analysis was done with Phylobayes, which outputs a datedist file that contains the age of every node for every tree. Obviously this is formatted in a more complicated way in Beast, within nexus files.

I've found this, https://raw.githubusercontent.com/taming-the-beast/FBD-tutorial/master/scripts/plot_geoscaled_tree.R
but it plots the node bars and the phylogeny, none of which I'm trying to do. So far in R, I can import the posterior distribution of trees from Beast, and remove the burnin. I can get the $tip.label for every tree, but after that am a bit unsure. I realize there is getMRCA in ape and findMRCA in phytools but I am not sure how to get these correctly for all nodes? My attempt to use getMRCA returns the same age for every tree (for just one node, the root), which is definitely not right. And then I'm not sure how to do it for other nodes...my Phylobayes script points at each node of choice to make the graph.

Here is how it looked for the Phylobayes script:

MyTrees = read.tree("filename.datedist")

tempTree = MyTrees[[1]]

#this gets the node number
tempTree$node.label = seq(length(tempTree$node.label)) 

#creates an empty matrix to put the node dates but you need the node numbers 

nodeDates = matrix(ncol = length(MyTrees[[1]]$node.label),nrow=1000000) 

#this loop extracts the dates from the node labels in the datedist file

n=0
for (i in MyTrees){
  dates = as.numeric(i$node.label)
  n=n+1
  nodeDates[n,]=dates
}

I was wondering if anyone on this list has ever done this or published an example a bit closer to what I am trying to do?

Thanks!