Geographical Heat Maps in R

I go crazy for fancy data visualizations in R, and a figure in a recent publication has had me wondering if there is an easy way to incorporate density distributions (or as in their case, a distribution of f4 statistics, that are often used to estimate genomic admixture), plotted as heat-maps overlaid on geographical maps. As it turns out, it’s a breeze to make these plots in R, using the ggplot2 package.

Denisovan ancestry in modern human populations, measured as an f4 statistic distribution across the world. Image courtesy: Figure 2 from Qin and Stoneking (2015) http://dx.doi.org/10.1093/molbev/msv141
Denisovan ancestry in modern human populations, measured as an f4 statistic distribution across the world. Image courtesy: Figure 2 from Qin and Stoneking (2015) http://dx.doi.org/10.1093/molbev/msv141

So I simulated some geographical data (which you should switch out for your real GPS coordinates), and obtained posterior density estimates for effective population sizes that I had previously used in another tutorial as an example. Importantly, you will have to edit (1) the density data, here represented by “hpars”, (2) the positions table – here I have simulated some data around coordinates for Philadelphia, PA, which you should replace with your own geographical coordinates for your data of interest from (1).

library(ggmap) #Load libraries
library(ggplot2)
hpars <- read.table("https://sites.google.com/site/arunsethuraman1/teaching/hpars.dat?revision=1") #Read in the density data
positions <- data.frame(lon=rnorm(20000, mean=-75.1803458, sd=0.05),
lat=rnorm(10000,mean=39.98352197, sd=0.05)) #Simulate some geographical coordinates #Switch out for your data that has real GPS coords
map <- get_map(location=c(lon=-75.1803458,
lat=39.98352197), zoom=11, maptype='roadmap', color='bw') #Get the map from Google Maps
ggmap(map, extent = "device") +
geom_density2d(data = positions, aes(x = lon, y = lat), size = 0.3) +
stat_density2d(data = positions,
                 aes(x = lon, y = lat, fill = ..level.., alpha = ..level..), size = 0.01,
                 bins = 16, geom = "polygon") + scale_fill_gradient(low = "green", high = "red") +
  scale_alpha(range = c(0, 0.3), guide = FALSE) #Plot

And voila! Do give it a try with your data, and leave your comments below.

Geographical heatmap produced using an effective population size density distribution, overlaid on a geographical map of Philadelphia.
Geographical heatmap produced using an effective population size density distribution, overlaid on a geographical map of Philadelphia.
This entry was posted in bioinformatics, howto, R, software and tagged , . Bookmark the permalink.