MPB：南土所褚海燕组-​​利用种分布模型绘制微生物分布图谱

为进一步提高《微生物组实验手册》稿件质量，本项目新增大众评审环节。文章在通过同行评审后，采用公众号推送方式分享全文，任何人均可在线提交修改意见。公众号格式显示略有问题，建议电脑端点击文末阅读原文下载PDF审稿。在线文档(https://kdocs.cn/l/cL8RRqHIL)大众评审页面登记姓名、单位和行号索引的修改建议。修改意见的征集截止时间为推文发布后的72小时，文章将会结合有建设性的修改意见进一步修改后获得DOI在线发表，同时根据贡献程度列为审稿人或致谢。感谢广大同行提出宝贵意见。

利用种分布模型绘制微生物分布图谱

Illustrating Microbial Distribution Map Utilizing Species Distribution Modeling 

李云涛¹，褚海燕^{1, 2,} *

¹中国科学院南京土壤研究所，土壤与农业可持续发展国家重点实验室，南京，江苏省，210008；

²中国科学院大学，北京，100049

*通讯作者邮箱：[e-mail]

摘要：当以点向面进行变量预测时，常使用插值法来实现。但是，环境微生物是异质性极高的复杂有机群体，使用插值法对其多样性或群落结构的预测效果通常很差。另一方面，环境微生物极易受到环境及气候因子的影响，这一特性使得我们可以通过环境或气候因素来对微生物的地理分布进行高效而准确的预测。本文以中国华北平原麦田土壤细菌的地理分布研究为示例，基于环境因子对微生物群落的驱动作用，利用种分布模型，在R软件中绘制微生物的分布图谱，给出了微生物分布图谱的标准化绘制及验证流程，为环境微生物的生物地理学研究提供了新的思路。

关键词：微生物群落，种分布模型，分布图谱

仪器设备

1.普通Windows系统个人电脑，内存8 G，需求硬盘空间 (含软件)约500 M

软件

1.R (v3.5.1)，所需依赖包：sp、raster、rgdal、dismo、ggplot2和ggthemes

注：本教程是基于已经在个人电脑上安装好的相关软件和依赖包进行的。如果安装出现问题，请参考以下链接：

https://cran.r-project.org/web/packages/sp/index.html

https://cran.r-project.org/web/packages/raster/index.html

https://cran.r-project.org/web/packages/rgdal/index.html

https://cran.r-project.org/web/packages/dismo/index.html

https://cran.r-project.org/web/packages/ggplot2/index.html

https://cran.r-project.org/web/packages/ggthemes/index.html

实验步骤

1.数据准备

本分析中需要用到三个数据，即微生物群落数据 (如多样性指数)，土壤理化性状数据 (与采样点一一对应)，以及待分析区域内的背景土壤多边形数据（一般为shapefile格式）。

注：背景土壤理化性状多边形数据可能需要联系相关专业单位或课题组获得。本文中所使用的数据来源于国家土壤信息服务平台（http://www.soil.csdb.cn/）。

2.数据导入及土壤理化性状筛选

2.1导入观测的微生物多样性数据及土壤理化性状数据（图1）

obs = read.table ("observation_data.txt", header = T)

图1. 微生物群落指标及土壤理化性状数据集

2.2使用广义线性模型筛选出与微生物多样性显著相关的土壤理化性状，以PD值为例

m = glm (PD~., data = obs)

summary (m)

3.读取背景土壤多边形数据

shpmap = readOGR (dsn = "7Province", layer = "bg_soil_data")

注：其中dsn指代包含所有背景地图数据图层文件的文件夹，layer指代文件夹中图层文件的名字 (不含扩展名)。

如有必要，将地图投影转换为常用的WGS84坐标系：

shpmap84 = spTransform (shpmap, CRS ("+proj = longlat +ellps = WGS84"))

4.提取需要的土壤背景栅格数据

4.1创建一个空白栅格对象，其尺度等于背景地图的尺度

r1 = raster (extent(shpmap84))

4.2设置空白栅格的分辨率，其单位与背景地图保持一致，数值可自行调整

res(r1) = c (1/50, 1/50)

4.3将背景多边形数据的其中一个土壤理化性状提取到新创建的空白栅格对象中（图2）

r1 = rasterize (shpmap84, r1, field = "pH")

plot (r1)

图2. 栅格对象r1中的土壤背景pH分布

注：选择的土壤理化性状须与步骤2中所筛选的因子完全一致，如有多个性状，则应创建多个空白栅格 (r1, r2, r3…)并重复本步骤。

5.种分布模型构建

5.1合并步骤4中的所有栅格对象

predictors = stack (r1, r2...)

5.2查看合并后每个土壤理化性状对应的名字

Names (predictors)

注：合并后，根据stack ()函数中每个栅格对象输入的顺序，predictors对象中土壤理化性状的名字将会变为layer.1、layer.2...，因此还需将观测数据中的对应土壤理化性状也更改为layer.1、layer.2...，再次建立步骤2中的广义线性模型，以保持变量名的统一。

5.3使用合并后的土壤栅格数据以及广义线性模型对象来进行微生物多样性的种分布模型构建

p = predict (predictors, m)

6.导出微生物多样性预测值的栅格数据

将步骤5中的种分布模型栅格对象输出为文本格式（图3）

p.xy = as.data.frame (p, xy = TRUE)

write.table (p.xy, file = "PD_prediction.txt", sep = "\t", col.names = TRUE)

注：输出的文本文件包含每个栅格的经纬度坐标以及预测的微生物多样性数值（x表示经度，y表示纬度），这里可以保存并过滤少量异常值点。

图3. 种分布模型预测结果展示

7.将种分布模型绘制成图

7.1输入文本格式的栅格数据

p.xy = read.table ("PD_prediction.txt", header = T)

7.2使用ggplot2包生成图像（图4）

g = ggplot (p.xy, aes(x, y))

g + geom_tile (aes (fill = PD)) + scale_fill_gradientn (colours = c ("blue", "green", "red")) + theme_few () + xlab (label = "Latitude") + ylab (label = "Longitude")

图4. ggolot2包绘制微生物群落多样性的分布图谱

8.交叉验证种分布模型的效力

注：在交叉验证时，首先将观测数据 (示例中为243个)随机分为两部分，其中2/3 (162个)用于广义线性模型构建，称为模型数据集，另外1/3 (81个)用于验证模型的效力，称为验证数据集。

8.1对观测数据随机取样并建立模型数据集及验证数据集

w = sample(1:243,size = 162)

mod_p = obs[which(rownames(obs) %in% w),]

write.table (mod_p, file = "modeling_part.txt ", sep = "\t")

val_p = obs[-which(rownames(obs) %in% w),]

write.table (val_p, file = "validation_part.txt", sep = "\t")

8.2读取模型数据集并建立广义线性模型

mod = read.table ("modeling_part.txt", header = T)

m1 = glm (PD ~ pH + ..., data = mod)

注：这里选择的土壤理化性状须与步骤2中保持一致。

8.3将验证数据集中的土壤理化性状数据代入上述方程，生成对应的多样性指数预测值 (该步骤不再做展示)然后将验证数据集中的预测值和观测值一同读取

crsval = read.table ("PD_obs_pre.txt", header = T)

8.4使用线性回归方程考察预测值与观测值的回归关系，R²即为交叉验证的效力，方程的斜率和截距可用于画图

m.val = lm (obs_PD ~ pre_PD, data = crsval)

summary (m.val)

8.5使用ggplot2包生成交叉验证图（图5）

g2=ggplot (crsval, aes (x = obs_PD,y = pre_PD))

g2 + geom_point (size = 3) + theme_few () + geom_abline (intercept = 124.73, slope = 0.6108) + theme (text = element_text (size = 20)) + labs (x = "Observed phylogenetic diversity", y = "Predicted phylogenetic diversity")

图5. ggplot2包绘制种分布模型效力的交叉验证结果

注意事项

种分布模型的优势和局限

优势：只需有限数量的调查样点，就能够高效、准确地预测目标范围内的各项微生物指标，包括多样性、群落结构、相对丰度等。

局限：研究范围中的生境差异性不宜过大，不同生境中的驱动因子存在差异，整体绘制微生物分布图谱时会使预测效力降低。此外，在进行种分布模型预测前，建议先分析微生物群落构建过程。如果随机性过程主导了微生物群落组装过程，那么通过该方法进行微生物地理分布预测的可靠性会降低。

参考文献

1.Hijmans, R., Phillips, S., Leathwick, J. & Elith, J.. (2014). Dismo: species distribution modeling.

2.Shi, Y., Li, Y., Yuan, M., Adams, J. M., & Chu, H.. (2019). A biogeographic map of soil bacterial communities in wheats field of the north china plain. Soil Ecology Letters, 1(1), 50-58.

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文下载PDF审稿，或浏览器直接访问下载链接：http://210.75.224.110/github/MicrobiomeProtocol/04Review/210324/2003667HaiyanChu939029/Protocol2003667.pdf