Formatting phylogeny to map projection (`phylo.to.plot`, or alternate method) in R

Question

I am hoping someone can help me with the formating from phylo.to.plot() or suggest another method that can produce a similar output.

I have followed tutorial(s) here to produce an output but it seems difficult to alter the resulting figures.

Briefly these are my questions. I will expand further below.

How to plot a subregion of a "WorldHires" map, not entire region?
Change the shape of the points on the map, but maintain the colour?
Add gradient of continuous variable to map

Reproducible example:

Here is a very basic tree with some randomly assigned geographic locations

myTree <- ape::read.tree(text='((A, B), ((C, D), (E, F)));')
plot(myTree)

# It needs to be rooted for `phylo.to.map()` to work
myTree$branch.length = NULL
rooted_cladogram = ape::compute.brlen(myTree)

# Sample information
Sample <- c("A","B","C","D","E","F")
coords <- matrix(c(56.001966,57.069417,50.70228, 51.836213, 54.678997, 54.67831,-5.636926,-2.47805,-3.8975018, -2.235444,-3.4392211, -1.751833), nrow=6, ncol=2)

rownames(coords) <- Sample  
head(coords)

## Plot phylo.to.map
obj<-phylo.to.map(rooted_cladogram,coords,database="worldHires", regions="UK",plot=FALSE,xlim=c(-11,3), ylim=c(49,59),direction="rightwards")
plot(obj,direction="rightwards",fsize=0.5,cex.points=c(0,1), lwd=c(3,1),ftype="i")

Plot output here:

Question 1: How do I plot a subregion of a "WorldHires" map, not the entire region?

I would like to only have mainland Britain which is a subregion of the "UK" in the WorldHires database. To access it normally I would do:

map1 <- ggplot2::map_data(map = "worldHires", region = c("UK"),xlim=c(-11,3), ylim=c(49,59))
GB <- subset(map1, subregion=="Great Britain")
# Plot
GB_plot<- ggplot(GB )+
  geom_polygon(aes(x = long, y = lat, group = group), fill = "white", colour = "black")+
theme_classic()+
  theme(axis.line=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        axis.title=element_blank(),
        panel.border = element_blank())

Which looks like this:

I have tried but it ignore the subregion argument.

obj<-phylo.to.map(ttree,coords,database="worldHires", regions="UK", subregion="Great Britain",plot=FALSE,xlim=c(-11,3), ylim=c(49,59),direction="rightwards")

Is there a way to provide it directly with a map instead of using WorldHires?

Question 2: How do I change the shape of the points on the map but keep maintain the colour?

I want to use shapes on the map to indicate the 3 major clade on my tree geographically. However, when I add a pch argument in, it correctly changes the shapes but the points then become black instead of following the colour that they were before. The lines from the tree to the map maintain the colour, it is just the points themselves that seem to turn black.

This is how I have tried to change the shape of the points:

# Original code - points
cols <-setNames(colorRampPalette(RColorBrewer::brewer.pal(n=6, name="Dark2"))(Ntip(myTree)),myTree$tip.label)

obj<-phylo.to.map(rooted_cladogram,coords,database="worldHires", regions="UK",plot=FALSE,xlim=c(-11,3), ylim=c(49,59),direction="rightwards")
plot(obj,direction="rightwards",fsize=0.5,cex.points=c(0,1), colors=cols,lwd=c(3,1),ftype="i")

Point and lines are coloured. I would like to change the shape of points

# Code to change points = but points are no longer coloured
shapes <- c(rep(2,2),rep(1,2),rep(0,2))

obj<-phylo.to.map(rooted_cladogram,coords,database="worldHires", regions="UK",plot=FALSE,xlim=c(-11,3), ylim=c(49,59),direction="rightwards")
plot(obj,direction="rightwards",fsize=0.5,cex.points=c(0,1), colors=cols,pch=shapes,lwd=c(3,1),ftype="i")

Output: The shapes are changed but they are no longer coloured in:

Question 3: How do I add a gradient to the map?

Given this fake dataset, how to I create a smoothed gradient of the value variable?

Any help and advice on this would be very much appreciated. It would also be useful to know how to change the size of points

Thank you very much in advance, Eve

Q1: it doesn't appear that the package author added `subregion=` as an argument in `phylo.to.map()`. One very hacky approach to dealing with this is to subset the coordinate pairs you want to plot. For instance, `obj$map$x <- obj$map$x[1450:20000]` and `obj$map$y <- obj$map$y[1450:20000]` — Skaqqs, Jul 02 '21 at 20:22
Q2: set `cex.points = 0` in `plot()`, then add `points(x = coords[,2], y = coords[,1], pch = shapes, col = cols, cex = 1) after plot()`. Can also change the size to whatever you want with `cex`. I'm not sure what you mean by 'add a gradient', but I may be able to come up with a half-baked solution if you elaborate :) — Skaqqs, Jul 02 '21 at 20:36

Skaqqs · Accepted Answer · 2021-07-22T20:25:47.590

I improved (somewhat) on my comments by using the map you made in your question. Here's the code:

library(mapdata)
library(phytools)
library(ggplot2)

myTree <- ape::read.tree(text='((A, B), ((C, D), (E, F)));')
plot(myTree)

# It needs to be rooted for `phylo.to.map()` to work
myTree$branch.length = NULL
rooted_cladogram = ape::compute.brlen(myTree)

# Sample information
Sample <- c("A","B","C","D","E","F")
coords <- matrix(
  c(56.001966,
    57.069417,
    50.70228,
    51.836213,
    54.678997,
    54.67831,
    -5.636926,
    -2.47805,
    -3.8975018,
    -2.235444,
    -3.4392211,
    -1.751833),
  nrow=6,
  ncol=2)

rownames(coords) <- Sample  
head(coords)

obj <- phylo.to.map(
  rooted_cladogram,
  coords,
  database="worldHires",
  regions="UK",
  plot=FALSE,
  xlim=c(-11,3),
  ylim=c(49,59),
  direction="rightwards")

# Disable default map
obj2 <- obj
obj2$map$x <- obj$map$x[1]
obj2$map$y <- obj$map$y[1]

# Set plot parameters
cols <- setNames(
  colorRampPalette(
    RColorBrewer::brewer.pal(n=6, name="Dark2"))(Ntip(myTree)),myTree$tip.label)
shapes <- c(rep(2,2),rep(1,2),rep(0,2))
sizes <- c(1, 2, 3, 4, 5, 6)

# Plot phylomap
plot(
  obj2,
  direction="rightwards",
  fsize=0.5,
  cex.points=0,
  colors=cols,
  pch=shapes,
  lwd=c(3,1),
  ftype="i")

# Plot new map area that only includes GB
uk <- map_data(
  map = "worldHires",
  region = "UK")
gb <- uk[uk$subregion == "Great Britain",]
points(x = gb$long,
       y = gb$lat,
       cex = 0.001)

# Plot points on map
points(
  x = coords[,2],
  y = coords[,1],
  pch = shapes,
  col = cols,
  cex = sizes)

e: Use sf object instead of points to illustrate GB. It is tough to provide more advice beyond this on how to add symbology for your spatially varying variable, but sf is popular and very well documented, e.g. https://r-spatial.github.io/sf/articles/sf5.html. Let me know if you have any other questions!

ee: Added lines to plot name and symbol on tips.

eee: Added gradient dataset to map.

library(phytools)
library(mapdata)
library(ggplot2)
library(sf)

myTree <- ape::read.tree(text='((A, B), ((C, D), (E, F)));')
plot(myTree)

# It needs to be rooted for `phylo.to.map()` to work
myTree$branch.length = NULL
rooted_cladogram = ape::compute.brlen(myTree)

# Sample information
Sample <- c("A","B","C","D","E","F")
coords <- matrix(c(56.001966,57.069417,50.70228, 51.836213, 54.678997, 54.67831,-5.636926,-2.47805,-3.8975018, -2.235444,-3.4392211, -1.751833), nrow=6, ncol=2)
rownames(coords) <- Sample  
head(coords)

obj <- phylo.to.map(
  rooted_cladogram,
  coords,
  database="worldHires",
  regions="UK",
  plot=FALSE,
  xlim=c(-11,3),
  ylim=c(49,59),
  direction="rightwards")

# Disable default map
obj2 <- obj
obj2$map$x <- obj$map$x[1]
obj2$map$y <- obj$map$y[1]

## Plot tree portion of map

# Set plot parameters
cols <- setNames(
  colorRampPalette(
    RColorBrewer::brewer.pal(n=6, name="Dark2"))(Ntip(myTree)),myTree$tip.label)
shapes <- c(rep(2,2),rep(1,2),rep(0,2))
sizes <- c(1, 2, 3, 4, 5, 6)

# Plot phylomap
plot(
  obj2,
  direction="rightwards",
  fsize=0.5,
  cex.points=0,
  colors=cols,
  pch=shapes,
  lwd=c(3,1),
  ftype="i")

tiplabels(pch=shapes, col=cols, cex=0.7, offset = 0.2)
tiplabels(text=myTree$tip.label, col=cols, cex=0.7, bg = NA, frame = NA, offset = 0.2)

## Plot GB portion of map

# Plot new map area that only includes GB
uk <- map_data(map = "worldHires", region = "UK")
gb <- uk[uk$subregion == "Great Britain",]

# Convert GB to sf object
gb_sf <- st_as_sf(gb, coords = c("long", "lat"))

# Covert to polygon
gb_poly <- st_sf(
  aggregate(
    x = gb_sf$geometry,
    by = list(gb_sf$region),
    FUN = function(x){st_cast(st_combine(x), "POLYGON")}))

# Add polygon to map
plot(gb_poly, col = NA, add = TRUE)

## Load and format gradient data as sf object

# Load data
g <- read.csv("gradient_data.txt", sep = " ", na.strings = c("NA", " "))
# Check for, then remove NAs
table(is.na(g))
g2 <- g[!is.na(g$Lng),]
# For demonstration purposes, make dataset easier to manage
# Delete this sampling line to use the full dataset
g2 <- g2[sample(1:nrow(g2), size = 1000),]
# Create sf point object
gpt <- st_as_sf(g2, coords = c("Lng", "Lat"))

## Set symbology and plot 

# Cut data into 5 groups based on "value"
groups <- cut(gpt$value,
              breaks = seq(min(gpt$value), max(gpt$value), len = 5),
              include.lowest = TRUE)
# Set colors
gpt$colors <- colorRampPalette(c("yellow", "red"))(5)[groups]
# Plot
plot(gpt$geometry, pch = 16, col = gpt$colors, add = TRUE)

## Optional legend for gradient data

# Order labels and colors for the legend
lev <- levels(groups)
# Used rev() here to make colors in correct order
fil <- rev(levels(as.factor(gpt$colors)))
legend("topright", legend = lev, fill = fil, add = TRUE)

## Plot sample points on GB

# Plot points on map
points(
  x = coords[,2],
  y = coords[,1],
  pch = shapes,
  col = cols,
  cex = sizes)

see here for more info on gradient symbology and legends: R: Gradient plot on a shapefile

This is superbly helpful - thank you! Do you know how I would fill the shapes with the colours as well? Thanks again — QPaps, Jul 03 '21 at 16:14
In terms of the gradient - for example if I have a continuous variable for the whole of mainland Britain - how would I add this to the map background? — QPaps, Jul 03 '21 at 16:21
You are welcome! You can change your points from open to closed (filled) with `pch`, see http://www.sthda.com/english/wiki/r-plot-pch-symbols-the-different-point-shapes-available-in-r. Adding color based on a variable is certainly possible, but we might have to change a couple things (like using an sf object instead of points to illustrate GB). Would you be able to share those data as part of your original question? — Skaqqs, Jul 03 '21 at 18:19
Hi @Skaqqs - thank you for the shape info. Unfortunately the dataset is rather large and I do not think I am allowed to share it. If possible could you show me how to make this but using an sf objcet so I can attempt the gradient myself? Thank you so much — QPaps, Jul 05 '21 at 10:47
Just FYI - I have just used your code on my full dataset and the map is made up of lots of little dots, so I altered the `points(x = gb$long,y = gb$lat,cex = 0.001))` to `lines()` instead. — QPaps, Jul 05 '21 at 10:48
I updated by original answer to show how you could use an `sf` object in your map. Good note on `lines()`; that seems much better! — Skaqqs, Jul 05 '21 at 14:51
This is amazing, thank you so much for all your help - I would upvote a lot more times if only it were possible! :D — QPaps, Jul 05 '21 at 17:28
Sorry to ask you for more help. You wouldn't know how to add the corresponding symbol to the end of the tree? So it would look like tree + symbol + name? I have looked at `tiplabels(pch=shapes, col=cols, cex=0.7,offset=-0.25)` but if I would ratehr it was positioned just before the writing and after the tree. I can't work out how to shift the labels from the `phylo.to.map`. Thanks again — QPaps, Jul 06 '21 at 09:42
Sure! Just add `tiplabels(pch=shapes, col=cols, cex=0.7, offset = 0.2); tiplabels(text=myTree$tip.label, col=cols, cex=0.7,offset=0, bg = NA, frame = NA)` after you plot your tree. — Skaqqs, Jul 06 '21 at 12:10
The labels are already printed with the `plot(obj...)` so anything on the end is just written over the top? Is there a way to not print the labels in the first instance? — QPaps, Jul 06 '21 at 12:25
I've never used phylomap before this, so I can't say for sure how the labelling works. But from what I can tell, adding `tiplabels()` after plotting `obj` overwrites labels in the output. I updated my answer and posted the result. The colors and shapes are a bit jumbled for some reason, but that shouldn't be hard to organize. — Skaqqs, Jul 06 '21 at 12:36
Thank you so much - you have no idea how many headaches you have saved me. Thank you again :D — QPaps, Jul 06 '21 at 12:45
Hi @Skaqqs - my apologies to get back to you but I cant seem to get the graident to work and I do not know who else to ask. I have added a dropbox transfer link to the datapoints with a made up variable. Would you mind helping me plot a smotthed gradient of this in the map on the plot? — QPaps, Jul 21 '21 at 10:49
hello! Just writing to say that I see your comment and will have time to take a look at this next week. If you solve this before then, let me know. Or if you need help sooner, feel free to post as another question... someone may have a quick answer! — Skaqqs, Jul 21 '21 at 19:01
Hi @QPaps! I updated my answer to show how you could illustrate your gradient data as points. It is not the prettiest, but I hope it gets you on the right track. Let me know if you have any questions! — Skaqqs, Jul 22 '21 at 20:25
Thank you so much, this helps alot @Skaqqs! I will have a go at smoothing the gradient thanks again — QPaps, Jul 24 '21 at 16:53

Formatting phylogeny to map projection (`phylo.to.plot`, or alternate method) in R

1 Answers1