ggcor替代包：linkET，相关图，mantel test可视化

本文首发于公众号：医学和生信笔记

“
医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

在几年前出现了一个ggcor包，可以用来可视化mantel test的结果，最开始还可以通过cran安装，不过后来也不行了，而且这个包由于一些原因已经停止维护了，最近的更新是2年前了！

但是那张图却一直很风靡。。。其实原作者已经开发了新的包用于可视化mantel test，名字叫linkET，只是由于缺少宣传，大家知道的比较少。

善于搜索一搜就能搜到，我在之前的可能是最适合初学者的R包安装教程，视频中提到了这个linkET，但是大家不愿意看，一个劲的问封面图是怎么画的，我真是服了。。。

所以，今天专门介绍用于mantel test可视化的linkET包。

安装

首先是安装R包，这个包只能通过github安装，或者下载到本地安装，使用install.packages()100%报错，使用BiocManager::install()也是100%报错！

很多初学者最大的拦路虎绝对是R包安装，有些人宁愿花钱找tb，也不愿意自己学习一下，搞不懂！

# install.packages("devtools")
devtools::install_github("Hy4m/linkET", force = TRUE)
1

使用

一般的相关性分析是用于两列数据之间的，而mantel test 是用于两个矩阵的相关性检验，我在工作中用的很少，做微生物的小伙伴应该用的比较多，做肠道菌群、宏基因组的应该也会用到。

加载R包和数据：

library(linkET)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)

data("varechem", package = "vegan")
data("varespec", package = "vegan")
1

看看这两个数据长什么样，因为数就是图，图就是数！

class(varechem)
## [1] "data.frame"
class(varespec)
## [1] "data.frame"


glimpse(varechem)
## Rows: 24
## Columns: 14
## $ N         19.8, 13.4, 20.2, 20.6, 23.8, 22.8, 26.6, 24.2, 29.8, 28.1, 2…
## $ P         42.1, 39.1, 67.7, 60.8, 54.5, 40.9, 36.7, 31.0, 73.5, 40.5, 3…
## $ K         139.9, 167.3, 207.1, 233.7, 180.6, 171.4, 171.4, 138.2, 260.0…
## $ Ca        519.4, 356.7, 973.3, 834.0, 777.0, 691.8, 738.6, 394.6, 748.6…
## $ Mg        90.0, 70.7, 209.1, 127.2, 125.8, 151.4, 94.9, 45.3, 105.3, 11…
## $ S         32.3, 35.2, 58.1, 40.7, 39.5, 40.8, 33.8, 27.1, 42.5, 60.2, 3…
## $ Al        39.0, 88.1, 138.0, 15.4, 24.2, 104.8, 20.7, 74.2, 17.9, 329.7…
## $ Fe        40.9, 39.0, 35.4, 4.4, 3.0, 17.6, 2.5, 9.8, 2.4, 109.9, 4.6, …
## $ Mn        58.1, 52.4, 32.1, 132.0, 50.1, 43.6, 77.6, 24.4, 106.6, 61.7,…
## $ Zn        4.5, 5.4, 16.8, 10.7, 6.6, 9.1, 7.4, 5.2, 9.3, 9.1, 8.1, 10.2…
## $ Mo        0.30, 0.30, 0.80, 0.20, 0.30, 0.40, 0.30, 0.30, 0.30, 0.50, 0…
## $ Baresoil  43.90, 23.60, 21.20, 18.70, 46.00, 40.50, 23.00, 29.80, 17.60…
## $ Humdepth  2.2, 2.2, 2.0, 2.9, 3.0, 3.8, 2.8, 2.0, 3.0, 2.2, 2.7, 2.5, 2…
## $ pH        2.7, 2.8, 3.0, 2.8, 2.7, 2.7, 2.8, 2.8, 2.8, 2.8, 2.7, 2.9, 2…

glimpse(varespec)
## Rows: 24
## Columns: 44
## $ Callvulg  0.55, 0.67, 0.10, 0.00, 0.00, 0.00, 4.73, 4.47, 0.00, 24.13, …
## $ Empenigr  11.13, 0.17, 1.55, 15.13, 12.68, 8.92, 5.12, 7.33, 1.63, 1.90…
## $ Rhodtome  0.00, 0.00, 0.00, 2.42, 0.00, 0.00, 1.55, 0.00, 0.35, 0.07, 0…
## $ Vaccmyrt  0.00, 0.35, 0.00, 5.92, 0.00, 2.42, 6.05, 2.15, 18.27, 0.22, …
## $ Vaccviti  17.80, 12.13, 13.47, 15.97, 23.73, 10.28, 12.40, 4.33, 7.13, …
## $ Pinusylv  0.07, 0.12, 0.25, 0.00, 0.03, 0.12, 0.10, 0.10, 0.05, 0.12, 0…
## $ Descflex  0.00, 0.00, 0.00, 3.70, 0.00, 0.02, 0.78, 0.00, 0.40, 0.00, 0…
## $ Betupube  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.02, 0.00, 0.00, 0.00, 0…
## $ Vacculig  1.60, 0.00, 0.00, 1.12, 0.00, 0.00, 2.00, 0.00, 0.20, 0.00, 0…
## $ Diphcomp  2.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0…
## $ Dicrsp    0.00, 0.33, 23.43, 0.00, 0.00, 0.00, 0.03, 1.02, 0.30, 0.02, …
## $ Dicrfusc  1.62, 10.92, 0.00, 3.63, 3.42, 0.32, 37.07, 25.80, 0.52, 2.50…
## $ Dicrpoly  0.00, 0.02, 1.68, 0.00, 0.02, 0.02, 0.00, 0.23, 0.20, 0.00, 0…
## $ Hylosple  0.00, 0.00, 0.00, 6.70, 0.00, 0.00, 0.00, 0.00, 9.97, 0.00, 0…
## $ Pleuschr  4.67, 37.75, 32.92, 58.07, 19.42, 21.03, 26.38, 18.98, 70.03,…
## $ Polypili  0.02, 0.02, 0.00, 0.00, 0.02, 0.02, 0.00, 0.00, 0.00, 0.00, 0…
## $ Polyjuni  0.13, 0.23, 0.23, 0.00, 2.12, 1.58, 0.00, 0.02, 0.08, 0.02, 0…
## $ Polycomm  0.00, 0.00, 0.00, 0.13, 0.00, 0.18, 0.00, 0.00, 0.00, 0.00, 0…
## $ Pohlnuta  0.13, 0.03, 0.32, 0.02, 0.17, 0.07, 0.10, 0.13, 0.07, 0.03, 0…
## $ Ptilcili  0.12, 0.02, 0.03, 0.08, 1.80, 0.27, 0.03, 0.10, 0.03, 0.25, 0…
## $ Barbhatc  0.00, 0.00, 0.00, 0.08, 0.02, 0.02, 0.00, 0.00, 0.00, 0.07, 0…
## $ Cladarbu  21.73, 12.05, 3.58, 1.42, 9.08, 7.23, 6.10, 7.13, 0.17, 23.07…
## $ Cladrang  21.47, 8.13, 5.52, 7.63, 9.22, 4.95, 3.60, 14.03, 0.87, 23.67…
## $ Cladstel  3.50, 0.18, 0.07, 2.55, 0.05, 22.08, 0.23, 0.02, 0.00, 11.90,…
## $ Cladunci  0.30, 2.65, 8.93, 0.15, 0.73, 0.25, 2.38, 0.82, 0.05, 0.95, 2…
## $ Cladcocc  0.18, 0.13, 0.00, 0.00, 0.08, 0.10, 0.17, 0.15, 0.02, 0.17, 0…
## $ Cladcorn  0.23, 0.18, 0.20, 0.38, 1.42, 0.25, 0.13, 0.05, 0.03, 0.05, 0…
## $ Cladgrac  0.25, 0.23, 0.48, 0.12, 0.50, 0.18, 0.18, 0.22, 0.07, 0.23, 0…
## $ Cladfimb  0.25, 0.25, 0.00, 0.10, 0.17, 0.10, 0.20, 0.22, 0.10, 0.18, 0…
## $ Cladcris  0.23, 1.23, 0.07, 0.03, 1.78, 0.12, 0.20, 0.17, 0.02, 0.57, 0…
## $ Cladchlo  0.00, 0.00, 0.10, 0.00, 0.05, 0.05, 0.02, 0.00, 0.00, 0.02, 0…
## $ Cladbotr  0.00, 0.00, 0.02, 0.02, 0.05, 0.02, 0.00, 0.00, 0.02, 0.07, 0…
## $ Cladamau  0.08, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
## $ Cladsp    0.02, 0.00, 0.00, 0.02, 0.00, 0.00, 0.02, 0.02, 0.00, 0.07, 0…
## $ Cetreric  0.02, 0.15, 0.78, 0.00, 0.00, 0.00, 0.02, 0.18, 0.00, 0.18, 0…
## $ Cetrisla  0.00, 0.03, 0.12, 0.00, 0.00, 0.00, 0.00, 0.08, 0.02, 0.02, 0…
## $ Flavniva  0.12, 0.00, 0.00, 0.00, 0.02, 0.02, 0.00, 0.00, 0.00, 0.00, 0…
## $ Nepharct  0.02, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
## $ Stersp    0.62, 0.85, 0.03, 0.00, 1.58, 0.28, 0.00, 0.03, 0.02, 0.03, 0…
## $ Peltapht  0.02, 0.00, 0.00, 0.07, 0.33, 0.00, 0.00, 0.00, 0.00, 0.02, 0…
## $ Icmaeric  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.00, 0.00, 0…
## $ Cladcerv  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
## $ Claddefo  0.25, 1.00, 0.33, 0.15, 1.97, 0.37, 0.15, 0.67, 0.08, 0.47, 1…
## $ Cladphyl  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
1

接下来先进行mantel test检验，然后画图即可。

# 计算
mantel <- mantel_test(varespec, varechem,
                      spec_select = list(Spec01 = 1:7,
                                         Spec02 = 8:18,
                                         Spec03 = 19:37,
                                         Spec04 = 38:44)) %>% 
  mutate(rd = cut(r, breaks = c(-Inf, 0.2, 0.4, Inf), # 对相关系数进行分割，便于映射大小
                  labels = c("< 0.2", "0.2 - 0.4", ">= 0.4")),
         pd = cut(p, breaks = c(-Inf, 0.01, 0.05, Inf), # 对P值进行分割，便于映射颜色
                  labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))
## `mantel_test()` using 'bray' dist method for 'spec'.
## `mantel_test()` using 'euclidean' dist method for 'env'.

# 画图
qcorrplot(correlate(varechem), type = "lower", diag = FALSE) +
  geom_square() +
  geom_couple(aes(colour = pd, size = rd), # 这行代码是关键
              data = mantel, 
              curvature = nice_curvature()) +
  
  # 下面就是各种颜色和名称设置
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu")) +
  scale_size_manual(values = c(0.5, 1, 2)) +
  scale_colour_manual(values = color_pal(3)) +
  guides(size = guide_legend(title = "Mantel's r",
                             override.aes = list(colour = "grey35"), 
                             order = 2),
         colour = guide_legend(title = "Mantel's p", 
                               override.aes = list(size = 3), 
                               order = 1),
         fill = guide_colorbar(title = "Pearson's r", order = 3))
1

简单方便，快捷画图，出图效果也很好。

当然各种细节都是可以用ggplot2语法修改的，喜欢折腾的可以自己尝试下~

这个图看起来很复杂，但是展示的信息确实很多，不过既然你需要这张图，那我相信你应该理解这张图的意思~

最后，我还是要强调下基础知识的重要性！如果要经常用R，最好还是系统学习下，吃饭的工具怎么能马马虎虎呢？

本文首发于公众号：医学和生信笔记

“
医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

本文由 mdnice 多平台发布

相关阅读:
栈的基本操作
 Lit（三）：样式
 如何使用 Python 在终端中创建令人惊叹的图形
 以自主技术跃进的综合冲压的顶级制造商
 Anaconda一文入门笔记
 docker-network网络
 会议OA之我的审批(查询&签字）
公众号爆文写作怎么做？或许这些领域才适合你！
内存管理篇——物理内存的管理
 mongodb
原文地址：https://blog.csdn.net/Ayue0616/article/details/126869631