Science组合图表解读
陈威，水利部中国科学院水工程生态研究所，藻类生态学方向，R语言爱好者。
关于此图的讨论已经有一段时间了。我发现一个事实，对此图教程表现出强烈渴望的小伙伴名字后面都有“生态”二字。不管是土壤生态、草地生态还是水生态。非生态的大佬及吃瓜群众也被图形的美学及提供的丰富信息量所吸引。R小白的我也尝试着去还原文中的美图，但是一直进展缓慢。这几天，擂台赛似的相继出来了几种画法：“坐标法”，“python法”（原谅我也不知道用的什么法），“拼接法”，原图的效果大致都出来了：
(R学习教程看这里->R语言学习系列教程汇总 ）
但是有几位大佬好像把mantel test的mantel拼错了。于是乎，我决定以文献原文为基础，尝试结合 corrplot和 mantel test讲一个小故事。先结合图表简单介绍一下原文。题目：基于浮游动物群落揭示氨氮的生态阈值
图1，左下角展示了13种环境因子之间的相关性，右上角展示了浮游动物zooplankton（包括枝角类cladocera，桡足类copepod和轮虫rotiffer，其实还应该有原生动物，不知道这里为什么没有提），以及3个subsample枝角类，桡足类和轮虫群落与每一个环境因子之间的mantel 相关性。
回到图中左下部分，可以看出，总氮TN与硝氮NO3（，氨氮TAN（均为总氮的一部分）相关性非常高，化学需氧量CODMn和生物需氧量BOD5相关性非常高，还有其他。其实，首先应该试探性将所有环境因子作为解释变量进行初步分析，查看其方差膨胀因子(variance inflation factor，VIF)大小， 一般认为 VIF＞10时，因子间共线性明显，需对其缩减，原则为pearson相关系数大于0.75的所有因子只保留其中一个，具体选择过程就不赘述了。
回到图中右上部分，可以看出，与浮游动物总群落相关性较高的有TN，PO4，TAN（这三者的自相关也非常高，可以认为是一个因子吗？）。TAN与浮游动物总群落和桡足类的相关性大于0.3，与枝角类的相关性在 0.2~0.3，是所有环境因子里面最大的，我想这就是这篇文章主题的来源吧。
图2A，物种矩阵与环境矩阵的冗余分析（ RDA），揭示环境因子对物种群落的影响。可以看出TAN的箭头处往第一轴上做垂线时，是最长的。此处与图1的结果吻合。图2B则单独拿出TAN含量与其样本在RDA第一轴上的得分做回归，相关系数当然高。
可以说，以上的所有内容只为找到这个 TAN，剩下的内容就是拿TAN来做文章，此处就不说了。
回头看另一篇更早的文章，Structure and function of theglobal ocean microbiome，找出的关键因素是Temperature，套路几乎一模一样。
弄清楚了套路，接下来谈谈数据是怎么分析的，图是怎么画的吧。
1、首先得有两个矩阵，一个是物种矩阵，另一个是影响物种组成的环境因子矩阵，两个矩阵有相同的行名称(如果有的话)及行数量，且物种矩阵每一行的和不能为0，暂且分别命名为otu和env。
2、计算env矩阵的相关系数，并以图形的方式展示出来。
otu<-read.csv(“otu.csv”,head=T，row.names=1) #如果没有行名，则不需row.names=1
env<-read.csv(“env.csv”,head=T，row.names=1) #如果没有行名，则不需row.names=1
library(corrplot)
rdf<-cor(env)
#得到相关矩阵图
corrplot(rdf,method = "ellipse" ,type = "upper",addrect = 1,insig="blank",rect.col = "blue",rect.lwd = 2)
以上因子的排序是按照env中的原始排序。考虑到后面的操作，我们更愿意将相关性高的一类因子放在一起，因此可以加入参数order="AOE"，另外"FPC"，"hclust"也有类似的效果。但是需要特别注意的是，后续mantel test 结果表中的因子顺序需按此重新排序，以免发生错配。
corrplot(cor(env),method = "ellipse",type = "upper",order="AOE",addrect = 1,insig="blank",rect.col = "blue",rect.lwd = 2)
为方便，本文还是以默认排序进行演示。
3、计算总类群及分类群与各个环境因子的mantel相关系数及显著性。
这一部分的内容感谢群里的李金明老师提供的方案，我做得有限。首先看封装的完整函数：
mt<-function(otu,env){
+ library(vegan)
+ library(dplyr)
+ vars <- colnames(env)
+ models<-list()
+ for (i in seq_along(vars)){
+ otu_bray<-vegdist(otu,method = "bray")
+ env_dis<-vegdist(env[vars[i]],method = "euclidean")
+ model <- mantel(otu_bray,env_dis, permutations=999)
+ name <- vars[i]
+ statistic <- model$statistic
+ signif <- model$signif
+ models[[i]] <- data.frame(name = name, statistic = statistic, signif = signif, row.names = NULL)
+ }
+ models %>% bind_rows()
+ }函数的功能大概就是将env矩阵中的每一个环境因子（已通过筛选）与otu进行mantel test，并从返回的model中将相关系数statistic和p值signif提取出来，并按顺序返回到一个新的dataframe中。此函数大大降低了工作量，只需作者整理好完整物种矩阵及各个subsample矩阵。当然也可以是多个独立的，但是都与同一环境因子矩阵相关联的物种矩阵，但是会损失一些可比信息。导入函数后，我们直接运行：
mantRpTotal<-mt(otu,env)
#不同subsample最好不要重合
otu_sub1<-otu[,1:18]
otu_sub2<-otu[,-(1:18)]
#同样必须保证otu_sub1和otu_sub2中每一行的和不为0
mantRpTotal<-mt(otu,env)
mantRpsub1<-mt(otu_sub1,env)
mantRpsub2<-mt(otu_sub2,env)
#得到三组相关矩阵> mantRpTotal
name statistic signif
1 T 0.07099868 0.277
2 Cond 0.21553934 0.051
3 DO -0.03173941 0.540
4 pH 0.30217734 0.017
5 ORP -0.09269039 0.796
6 Chla 0.10956145 0.196
> mantRpsub1
name statistic signif
1 T 0.08162417 0.257
2 Cond 0.21331172 0.041
3 DO -0.08547460 0.700
4 pH 0.26544007 0.014
5 ORP -0.07304401 0.744
6 Chla 0.15479697 0.101
> mantRpsub2
name statistic signif
1 T 0.007865918 0.446
2 Cond 0.177458987 0.107
3 DO 0.012635374 0.407
4 pH 0.291917464 0.021
5 ORP -0.116943736 0.853
6 Chla 0.044119183 0.332利用此表做出文献图1中右上部分的线段图还不够，还需要另外输入坐标信息。此处参考厚蕴老师的方法。
set.seed(123)
n <- ncol(env)
grp <- c('Total', 'Sub1', 'Sub2') # 分组名称
subx <- c(-2, -1, 0) # 分组的X坐标
suby <- c(4.5, 2.5, 1) # 分组的Y坐标
df <- data.frame(
grp = rep(grp, each = n), # 分组名称，每个重复n次
subx = rep(subx, each = n), # 组X坐标，每个重复n次
suby = rep(suby, each = n), # 组Y坐标，每个重复n次
x = rep(0:(n - 1) - 0.5, 3), # 变量连接点X坐标
y = rep(n:1, 3), # 变量连接点Y坐标
)
> df
grp subx suby x y
1 Total -2 5.0 -0.5 6
2 Total -2 5.0 0.5 5
3 Total -2 5.0 1.5 4
4 Total -2 5.0 2.5 3
5 Total -2 5.0 3.5 2
6 Total -2 5.0 4.5 1
7 Sub1 -1 2.5 -0.5 6
8 Sub1 -1 2.5 0.5 5
9 Sub1 -1 2.5 1.5 4
10 Sub1 -1 2.5 2.5 3
11 Sub1 -1 2.5 3.5 2
12 Sub1 -1 2.5 4.5 1
13 Sub2 0 1.0 -0.5 6
14 Sub2 0 1.0 0.5 5
15 Sub2 0 1.0 1.5 4
16 Sub2 0 1.0 2.5 3
17 Sub2 0 1.0 3.5 2
18 Sub2 0 1.0 4.5 1
df2 <-rbind(mantRpTotal, mantRpsub1, mantRpsub2)
df_segment<-cbind(df,df2)
df_segment
name statistic signif grp subx suby x y
1 T 0.070998676 0.277 Total -2.0 4.5 -0.5 6
2 Cond 0.215539337 0.051 Total -2.0 4.5 0.5 5
3 DO -0.031739414 0.540 Total -2.0 4.5 1.5 4
4 pH 0.302177339 0.017 Total -2.0 4.5 2.5 3
5 ORP -0.092690387 0.796 Total -2.0 4.5 3.5 2
6 Chla 0.109561448 0.196 Total -2.0 4.5 4.5 1
7 T 0.081624169 0.257 Sub1 -1.0 2.5 -0.5 6
8 Cond 0.213311719 0.041 Sub1 -1.0 2.5 0.5 5
9 DO -0.085474598 0.700 Sub1 -1.0 2.5 1.5 4
10 pH 0.265440072 0.014 Sub1 -1.0 2.5 2.5 3
11 ORP -0.073044012 0.744 Sub1 -1.0 2.5 3.5 2
12 Chla 0.154796969 0.101 Sub1 -1.0 2.5 4.5 1
13 T 0.007865918 0.446 Sub2 0 1.0 -0.5 6
14 Cond 0.177458987 0.107 Sub2 0 1.0 0.5 5
15 DO 0.012635374 0.407 Sub2 0 1.0 1.5 4
16 pH 0.291917464 0.021 Sub2 0 1.0 2.5 3
17 ORP -0.116943736 0.853 Sub2 0 1.0 3.5 2
18 Chla 0.044119183 0.332 Sub2 0 1.0 4.5 14、此时需考虑线条美学的映射，按原文图的表示，并不是按数值大小完全映射，而是划分范围后映射，此处对此时的我来说是知识盲区，又一次参考厚蕴老师的案例。
df_segment <- df_segment %>%
mutate(
lcol = ifelse(signif <= 0.001, '#1B9E77', NA),
# p值小于0.001时，颜色为绿色，下面依次类推
lcol = ifelse(signif > 0.001 & signif <= 0.01, '#88419D', lcol),
lcol = ifelse(signif > 0.01 & signif <= 0.05, '#A6D854', lcol),
lcol = ifelse(signif > 0.05, '#B3B3B3', lcol),
lwd = ifelse(statistic >= 0.3,6, NA),
# statistic >= 0.3 时，线性宽度为6，下面依次类推
lwd = ifelse(statistic >= 0.15 & statistic < 0.3, 3, lwd),
lwd = ifelse(statistic < 0.15, 1, lwd)
)5、一切就绪，开始拼图。
corrplot(cor(env),method = "ellipse",type = "upper", addrect = 1,insig="blank",rect.col = "blue",rect.lwd = 2)
segments(df_segment$subx, df_segment$suby, df_segment$x, df_segment$y, lty = 'solid', lwd = df_segment$lwd, col = df_segment$lcol, xpd = TRUE)
结果已经很明显了，pH就是我们要找的那个目标。下面继续完善图形内容。
points(subx, suby, pch = 24, col = 'blue', bg = 'blue', cex = 2, xpd = TRUE)
points(x,y,pch=21,col="red",bg="red",cex=1,xpd=TRUE)
text(subx - 0.5, suby, labels = grp, adj = c(0.8, 0.5), cex = 1.2, xpd = TRUE)
图例的内容我就不搬了。
6、下面继续挖掘并聚焦已经找到的目标，pH。
library(vegan)
otu_rda<- rda(otu,env)
summary(otu_rda)
plot(otu_rda)
#获取各样本pH在RDA第一轴的坐标,然后与env$pH做散点图
pH_sca<- otu_rda$CCA$QR$qr$pH具体画法很简单，不赘述（不班门弄斧了）。
至此，上述图形就都完成了，小故事也讲完了，有不对地方敬请指出。
参考资料：
[1]. Yang, J., et al., Ecogenomics of ZooplanktonCommunity Reveals Ecological Threshold of Ammonia Nitrogen. Environ SciTechnol, 2017. 51(5): p. 3057-3064.
[2]Sunagawa, S., Coelho, L. P., Chaffron, S., Kultima, J.R., Labadie, K., Salazar, G., … & Bork, P. (2015). Structure and function of theglobal ocean microbiome. Science, 348(6237), 1261359-1261359.
[3]厚缊,诹图——照虎画猫site:http://houyun.xyz/post/20190421