阿越1229

V1

2022/09/15阅读:15主题:自定义主题1

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)

使用

一般的相关性分析是用于两列数据之间的,而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")

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

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


glimpse(varechem)
## Rows: 24
## Columns: 14
## $ N        <dbl> 19.8, 13.4, 20.2, 20.6, 23.8, 22.8, 26.6, 24.2, 29.8, 28.1, 2…
## $ P        <dbl> 42.1, 39.1, 67.7, 60.8, 54.5, 40.9, 36.7, 31.0, 73.5, 40.5, 3…
## $ K        <dbl> 139.9, 167.3, 207.1, 233.7, 180.6, 171.4, 171.4, 138.2, 260.0…
## $ Ca       <dbl> 519.4, 356.7, 973.3, 834.0, 777.0, 691.8, 738.6, 394.6, 748.6…
## $ Mg       <dbl> 90.0, 70.7, 209.1, 127.2, 125.8, 151.4, 94.9, 45.3, 105.3, 11…
## $ S        <dbl> 32.3, 35.2, 58.1, 40.7, 39.5, 40.8, 33.8, 27.1, 42.5, 60.2, 3…
## $ Al       <dbl> 39.0, 88.1, 138.0, 15.4, 24.2, 104.8, 20.7, 74.2, 17.9, 329.7…
## $ Fe       <dbl> 40.9, 39.0, 35.4, 4.4, 3.0, 17.6, 2.5, 9.8, 2.4, 109.9, 4.6, …
## $ Mn       <dbl> 58.1, 52.4, 32.1, 132.0, 50.1, 43.6, 77.6, 24.4, 106.6, 61.7,…
## $ Zn       <dbl> 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       <dbl> 0.30, 0.30, 0.80, 0.20, 0.30, 0.40, 0.30, 0.30, 0.30, 0.50, 0…
## $ Baresoil <dbl> 43.90, 23.60, 21.20, 18.70, 46.00, 40.50, 23.00, 29.80, 17.60…
## $ Humdepth <dbl> 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       <dbl> 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 <dbl> 0.55, 0.67, 0.10, 0.00, 0.00, 0.00, 4.73, 4.47, 0.00, 24.13, …
## $ Empenigr <dbl> 11.13, 0.17, 1.55, 15.13, 12.68, 8.92, 5.12, 7.33, 1.63, 1.90…
## $ Rhodtome <dbl> 0.00, 0.00, 0.00, 2.42, 0.00, 0.00, 1.55, 0.00, 0.35, 0.07, 0…
## $ Vaccmyrt <dbl> 0.00, 0.35, 0.00, 5.92, 0.00, 2.42, 6.05, 2.15, 18.27, 0.22, …
## $ Vaccviti <dbl> 17.80, 12.13, 13.47, 15.97, 23.73, 10.28, 12.40, 4.33, 7.13, …
## $ Pinusylv <dbl> 0.07, 0.12, 0.25, 0.00, 0.03, 0.12, 0.10, 0.10, 0.05, 0.12, 0…
## $ Descflex <dbl> 0.00, 0.00, 0.00, 3.70, 0.00, 0.02, 0.78, 0.00, 0.40, 0.00, 0…
## $ Betupube <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.02, 0.00, 0.00, 0.00, 0…
## $ Vacculig <dbl> 1.60, 0.00, 0.00, 1.12, 0.00, 0.00, 2.00, 0.00, 0.20, 0.00, 0…
## $ Diphcomp <dbl> 2.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0…
## $ Dicrsp   <dbl> 0.00, 0.33, 23.43, 0.00, 0.00, 0.00, 0.03, 1.02, 0.30, 0.02, …
## $ Dicrfusc <dbl> 1.62, 10.92, 0.00, 3.63, 3.42, 0.32, 37.07, 25.80, 0.52, 2.50…
## $ Dicrpoly <dbl> 0.00, 0.02, 1.68, 0.00, 0.02, 0.02, 0.00, 0.23, 0.20, 0.00, 0…
## $ Hylosple <dbl> 0.00, 0.00, 0.00, 6.70, 0.00, 0.00, 0.00, 0.00, 9.97, 0.00, 0…
## $ Pleuschr <dbl> 4.67, 37.75, 32.92, 58.07, 19.42, 21.03, 26.38, 18.98, 70.03,…
## $ Polypili <dbl> 0.02, 0.02, 0.00, 0.00, 0.02, 0.02, 0.00, 0.00, 0.00, 0.00, 0…
## $ Polyjuni <dbl> 0.13, 0.23, 0.23, 0.00, 2.12, 1.58, 0.00, 0.02, 0.08, 0.02, 0…
## $ Polycomm <dbl> 0.00, 0.00, 0.00, 0.13, 0.00, 0.18, 0.00, 0.00, 0.00, 0.00, 0…
## $ Pohlnuta <dbl> 0.13, 0.03, 0.32, 0.02, 0.17, 0.07, 0.10, 0.13, 0.07, 0.03, 0…
## $ Ptilcili <dbl> 0.12, 0.02, 0.03, 0.08, 1.80, 0.27, 0.03, 0.10, 0.03, 0.25, 0…
## $ Barbhatc <dbl> 0.00, 0.00, 0.00, 0.08, 0.02, 0.02, 0.00, 0.00, 0.00, 0.07, 0…
## $ Cladarbu <dbl> 21.73, 12.05, 3.58, 1.42, 9.08, 7.23, 6.10, 7.13, 0.17, 23.07…
## $ Cladrang <dbl> 21.47, 8.13, 5.52, 7.63, 9.22, 4.95, 3.60, 14.03, 0.87, 23.67…
## $ Cladstel <dbl> 3.50, 0.18, 0.07, 2.55, 0.05, 22.08, 0.23, 0.02, 0.00, 11.90,…
## $ Cladunci <dbl> 0.30, 2.65, 8.93, 0.15, 0.73, 0.25, 2.38, 0.82, 0.05, 0.95, 2…
## $ Cladcocc <dbl> 0.18, 0.13, 0.00, 0.00, 0.08, 0.10, 0.17, 0.15, 0.02, 0.17, 0…
## $ Cladcorn <dbl> 0.23, 0.18, 0.20, 0.38, 1.42, 0.25, 0.13, 0.05, 0.03, 0.05, 0…
## $ Cladgrac <dbl> 0.25, 0.23, 0.48, 0.12, 0.50, 0.18, 0.18, 0.22, 0.07, 0.23, 0…
## $ Cladfimb <dbl> 0.25, 0.25, 0.00, 0.10, 0.17, 0.10, 0.20, 0.22, 0.10, 0.18, 0…
## $ Cladcris <dbl> 0.23, 1.23, 0.07, 0.03, 1.78, 0.12, 0.20, 0.17, 0.02, 0.57, 0…
## $ Cladchlo <dbl> 0.00, 0.00, 0.10, 0.00, 0.05, 0.05, 0.02, 0.00, 0.00, 0.02, 0…
## $ Cladbotr <dbl> 0.00, 0.00, 0.02, 0.02, 0.05, 0.02, 0.00, 0.00, 0.02, 0.07, 0…
## $ Cladamau <dbl> 0.08, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
## $ Cladsp   <dbl> 0.02, 0.00, 0.00, 0.02, 0.00, 0.00, 0.02, 0.02, 0.00, 0.07, 0…
## $ Cetreric <dbl> 0.02, 0.15, 0.78, 0.00, 0.00, 0.00, 0.02, 0.18, 0.00, 0.18, 0…
## $ Cetrisla <dbl> 0.00, 0.03, 0.12, 0.00, 0.00, 0.00, 0.00, 0.08, 0.02, 0.02, 0…
## $ Flavniva <dbl> 0.12, 0.00, 0.00, 0.00, 0.02, 0.02, 0.00, 0.00, 0.00, 0.00, 0…
## $ Nepharct <dbl> 0.02, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
## $ Stersp   <dbl> 0.62, 0.85, 0.03, 0.00, 1.58, 0.28, 0.00, 0.03, 0.02, 0.03, 0…
## $ Peltapht <dbl> 0.02, 0.00, 0.00, 0.07, 0.33, 0.00, 0.00, 0.00, 0.00, 0.02, 0…
## $ Icmaeric <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.00, 0.00, 0…
## $ Cladcerv <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
## $ Claddefo <dbl> 0.25, 1.00, 0.33, 0.15, 1.97, 0.37, 0.15, 0.67, 0.08, 0.47, 1…
## $ Cladphyl <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…

接下来先进行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(-Inf0.20.4Inf), # 对相关系数进行分割,便于映射大小
                  labels = c("< 0.2""0.2 - 0.4"">= 0.4")),
         pd = cut(p, breaks = c(-Inf0.010.05Inf), # 对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.512)) +
  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))
mantel test可视化
mantel test可视化

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

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

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

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

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

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

分类:

其他

标签:

医学

作者介绍

阿越1229
V1

黄金矿工。