Amid basting my tofurky, here’s a follow-up to my previous post on quick-fix overlays of admixture plots on geographical maps in R. I recently discovered a wonderful R package called “phytools” from Liam Revell, which makes really neat phylogenetic trees (with several new functions added to the “ape” package). I highly recommend Liam’s blog which has a wealth of examples and hacks to make cool figures. “phytools” also has a function called “phylo.to.map” which lets you make some neat geophylogenies.
I will use the same example from my previous post (you’ll need the same gps.data and admixprops.data files) to show you how. You will also need a NEWICK tree string for your phylogeny – I am assuming here that you have a population tree created using your method of choice. Importantly, the ID’s of your GIS coordinates for each locale (here described by the “ID” column) should be the same as the labels in your NEWICK string. Eg. (3:0.6,(4:0.1,(1:0.03,2:0.03):0.08):0.5);
library(phytools)
library(plotrix)
phy <- "(3:0.6,(4:0.1,(1:0.03,2:0.03):0.08):0.5);"
tree <- read.newick(text=phy)
gps <- read.csv(“gps.data”,header=TRUE) #Read input files
admix <- read.csv(“admixprops.data”,header=TRUE)
Now you should have your tree, GPS coordinates, and admixture proportions read. To plot them all in one figure:
obj <- phylo.to.map(tree,gps[,2:3], type="phylogram",database="state")
plot(obj,colors=c("blue","red"),ftype="i", fsize=1.5,asp=1.2,split=c(0.5,0.5))
for (x in 1:nrow(gps)){
floating.pie(gps$Lon[x],gps$Lat[x],
c(admix$K1[x],admix$K2[x]),radius=admix$Num[x]/8,
col=c("red","blue")
}
And tada!
There are a lot of way cooler tools that do the same thing, including GeoPhyloBuilder and GenGIS, but hey – if you can do it in four lines of code in R, why not? Happy Thanksgiving from all of us at The Molecular Ecologist! This year, like every year, I am very thankful for R.