在 R 中格式化系统发育以映射投影(`phylo.to.plot` 或替代方法)

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

我希望有人可以帮助我从 phylo.to.plot() 格式化或建议另一种可以产生类似输出的方法。

我已经按照教程 here 生成了输出,但似乎很难更改结果数字。

简而言之,这些是我的问题。我将在下面进一步展开。

  1. 如何绘制“WorldHires”地图的一个子区域,而不是整个区域?
  2. 更改地图上点的形状,但保持颜色?
  3. 将连续变量的梯度添加到地图

可重现的例子:

这是一棵非常基本的树,其中包含一些随机分配的地理位置

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")

在此处绘制输出:

问题 1:如何绘制“WorldHires”地图的一个子区域,而不是整个区域?

我只想在 WorldHires 数据库中包含英国大陆,它是“英国”的一个子区域。要正常访问它,我会这样做:

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())

看起来像这样:

我试过了,但它忽略了子区域参数。

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")

有没有办法直接提供地图而不是使用 WorldHires?

问题2:如何改变地图上点的形状但保持颜色不变?

我想使用地图上的形状在地理上指示我的树上的 3 个主要进化枝。然而,当我添加一个 pch 参数时,它正确地改变了形状,但这些点随后变成了黑色,而不是遵循它们之前的颜色。从树到地图的线保持颜色,只是点本身似乎变黑了。

这就是我尝试改变点形状的方法:

# 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")

点和线是彩色的。我想改变点的形状

# 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")

输出:形状已更改但不再着色:

问题3:如何给地图添加渐变?

鉴于此 fake dataset,如何创建值变量的平滑梯度?

如有任何帮助和建议,我们将不胜感激。 知道如何更改点的大小也很有用

非常感谢您, 夏娃

通过使用您在问题中制作的地图,我(在某种程度上)改进了我的评论。这是代码:

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: 用sf对象代替点来说明GB。很难就如何为空间变化变量添加符号系统提供更多建议,但 sf 很受欢迎并且有据可查,例如https://r-spatial.github.io/sf/articles/sf5.html。如果您还有其他问题,请告诉我!

ee:在提示上添加了绘图名称和符号的线条。

eee:向地图添加了梯度数据集。

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)

有关渐变符号系统和图例的更多信息,请参阅此处:R: Gradient plot on a shapefile