提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  其他

当前栏目

跟着NBT学做图：样本地理信息图

2023-03-07 09:47:08 时间

今天我们来学习一下刘永鑫老师2019年发表在Nature Biotechnology上的文章NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice中的代码。

其代码和数据都已经在其github[1]免费分享，按下图操作即可全部打包。下载缓慢的朋友也可在公众号回复「20220902」获得压缩包。

https://github.com/microbiota/Zhang2019NBT

图1 a Diagram of original collection sites (44 countries) of indica (red) and japonica (blue) rice

图1 a主要介绍了本研究中样本籼稻（indica）和粳稻（japonica）的位置。

源代码

在作图之前我们需要准备「包含样本地理位置信息的表」，至少包含品种和经纬度。

首先我放上经过我详细注释后的代码。

library(dplyr)
library(maptools)
library(ggthemes)
library(ggplot2)
library(maps)

geotable = read.table("varieties_geo.txt", header = T, sep = "\t") # 读取地理信息
worldmap = map_data("world") # 引入世界地图数据
# 通过?map_data查看，仅能提供county,france,italy,nz,state,usa,world,world2等几款地图

fig1 = ggplot(geotable, aes(Longitude, Latitude, color = Subspecies)) +
  geom_polygon(data = worldmap, aes(x = long, y = lat, group = group, fill = NA), 
               color = "grey70", size = 0.25)+ # 绘制地图
  geom_point(size = 2, alpha = 0.5)+ # 绘制具体的位置点
  scale_colour_brewer(palette = "Set1") + # 设置主题颜色
  scale_fill_brewer(palette = "Set1")+  # 设置主题颜色
  scale_y_continuous(breaks = (-3:3)*30)+ # 将Y轴的刻度限制为-90～90
  scale_x_continuous(breaks = (-6:6)*30)+ # 将X轴的刻度限制为-180～180
  labs(x="Longitude", y="Latitude", colour = "Subspecies" ) + # 修改X轴、Y轴及图例
  theme_tufte() # 需要library(ggthemes)
fig1
ggsave(paste0(output_dir, "minicore-worldmap.pdf", sep=""), fig1, width = 9, height = 5)

源代码出的图

函数详解

map_data()

map_data()[2]是ggplot2的一个函数，用于将map包中的数据转换为适合ggplot2绘图的框架。

「主要用法」

map_data(map, region = ".", exact = FALSE, ...)

「主要参数」

  • map：可以提供的地图，包括county(United States County Map), france,italy, nz(New Zealand Basic Map),state(United States State Boundaries Map), usa(United States Coast Map), world(Low (mid) resolution World Map), world2(Pacific Centric Low resolution World Map)。
  • region：包含子区域的名称，默认包含所有子区域。
  • exact：将region视为正则表达式(FALSE)还是固定字符串(TRUE)。

「示例」

library(ggplot2)
states <- map_data("state") # 引入美国地图
arrests <- USArrests
names(arrests) <- tolower(names(arrests)) # 将arrests的列名改为小写
arrests$region <- tolower(rownames(USArrests)) # 将USArrests的行名转换为小写并作为新增的一列保存为region，用于后续合并

choro <- merge(states, arrests, sort = FALSE, by = "region") # 通过region合并states和arrests
ggplot(choro, aes(long, lat)) + # 根据经纬度绘制地图
  geom_polygon(aes(group = group, fill = assault)) +
  coord_map("albers", lat0 = 45.5, lat1 = 29.5) 
# coord_map() 将地球的一部分近似球面投影到一个平面的2D平面上
# albers(lat0,lat1)：equal-area, true scale on lat0 and lat1

USArrests[3]包含了1973年美国50个州每10万居民中因袭击(assault)、谋杀(murder)和强奸(rap)而被捕的人数，还给出了生活在城市地区的人口百分比(urbanpop，Percent urban population)。

arrests

不添加coord_map("albers", lat0 = 45.5, lat1 = 29.5)

geom_polygon()

geom_polygon()[4]非常类似于由geom_path()绘制的路径，只不过起始点和结束点相连，内部可以填充。group映射决定了哪些case被连接在一起成为一个多边形。

「主要用法」

geom_polygon(
  mapping = NULL,
  data = NULL,
  stat = "identity",
  position = "identity",
  rule = "evenodd",
  ...,
  na.rm = FALSE,
  show.legend = NA,
  inherit.aes = TRUE)

「主要参数」

  • mapping：同ggplot2
  • data：同ggplot2
  • stat：使用的统计转换
  • position：位置调整
  • rule：可选evenodd或winding，如果要绘制带孔的多边形，这个参数定义了如何解释孔的坐标，示例[5]。
  • na.rm：默认情况下（False），缺失值会被移除并警告；选择True，缺失值会被悄无声息的移除。
  • show.legend：默认不展示。
  • inherit.aes：如果选择False，将会覆盖默认的映射（aesthetics），而不是将它们结合。

「注意⚠️」

在源代码中，关于geom_polygon的参数group，有一点需要注意。

geom_polygon(data = worldmap, aes(x = long, y = lat, group = group, fill = NA), 
               color = "grey70", size = 0.25)

worldmap

在这里group=group，而观察后我们发现，每个region对应了一个group，比如Aruba对应group 1，Afghanistan对应group 2，enclave对应group 1627，既然这样，能否group=region呢？

group=region

此时你应该明白，地图数据中的group是有意义的，决定了连线的先后顺序，在其他地图包中也有该顺序。

New Zealand Basic Map

美化

我对代码进行了微调，包括地图背景颜色、字/点的大小、透明度，把图限制在了一定区间范围，同时修改了主题。

library(dplyr)
library(maptools)
library(ggthemes)
library(ggplot2)
library(maps)
library(ggsci)
geotable = read.table("varieties_geo.txt", header = T, sep = "\t") # 读取地理信息
worldmap = map_data("world") # 引入世界地图数据
# 通过?map_data查看，仅能提供county,france,italy,nz,state,usa,world,world2等几款地图

fig = ggplot(geotable, aes(Longitude, Latitude, color = Subspecies)) +
  geom_polygon(data = worldmap, aes(long, lat, group = group, fill = NA), 
               color = "grey70",fill = "#F5F5DC", size = 0.25) + # 绘制地图
  geom_point(size = 1.7, alpha = 0.8)+ # 绘制具体的位置点
  labs(x="Longitude", y="Latitude", color = "Species" ) + # 修改X轴、Y轴及图例
  coord_cartesian(xlim = c(-120,150), ylim = c(-40,70)) +
  theme_tufte() + # 需要library(ggthemes) 
  theme(axis.ticks = element_line(linetype = "blank"), # 去掉坐标轴
  axis.text = element_text(colour = NA), # 去掉坐标轴的数字标记
  panel.background = element_rect(fill = '#F0FFFF'), # 背景色
  axis.title = element_text(size = 13),
  legend.text = element_text(size = 10),
  legend.title = element_text(size = 12),
  legend.position = c(0.06, 0.075)) + # 修改图例位置
  scale_color_aaas() # 修改主题 
 
fig

如果你用上了他们的代码，请务必按照下列方法引用：

Jingying Zhang#, Yong-Xin Liu#, Na Zhang#, Bin Hu#, Tao Jin#, Haoran Xu, Yuan Qin, Pengxu Yan, Xiaoning Zhang, Xiaoxuan Guo, Jing Hui, Shouyun Cao, Xin Wang, Chao Wang, Hui Wang, Baoyuan Qu, Guangyi Fan, Lixing Yuan, Ruben Garrido-Oter, Chengcai Chu* & Yang Bai*. NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology. 2019. doi:10.1038/s41587-019-0104-4

参考资料

[1]

github链接: https://github.com/microbiota/Zhang2019NBT

[2]

map_data(): https://ggplot2.tidyverse.org/reference/map_data.html

[3]

USArrests: https://www.rdocumentation.org/packages/datasets/versions/3.6.2/topics/USArrests

[4]

geom_polygon(): https://ggplot2.tidyverse.org/reference/geom_polygon.html

[5]

geom_polygon() rule参数示例: https://rdrr.io/r/grid/grid.path.html

相关文章