|
??專注R語言在??生物醫(yī)學(xué)中的使用 在幾年前出現(xiàn)了一個(gè)ggcor包，可以用來可視化mantel test的結(jié)果，最開始還可以通過cran安裝，不過后來也不行了，而且這個(gè)包由于一些原因已經(jīng)停止維護(hù)了，最近的更新是2年前了！ 但是那張圖卻一直很風(fēng)靡。。。其實(shí)原作者已經(jīng)開發(fā)了新的包用于可視化mantel test，名字叫l(wèi)inkET，只是由于缺少宣傳，大家知道的比較少。 善于搜索一搜就能搜到，我在之前的可能是最好用的R包安裝教程，視頻中提到了這個(gè)linkET，但是大家不愿意看，一個(gè)勁的問封面圖是怎么畫的，我真是服了。。。 所以，今天專門介紹用于mantel test可視化的linkET包。 安裝 首先是安裝R包，這個(gè)包只能通過github安裝，或者下載到本地安裝，使用install.packages()100%報(bào)錯(cuò)，使用BiocManager::install()也是100%報(bào)錯(cuò)！ 很多初學(xué)者最大的攔路虎絕對是R包安裝，有些人寧愿花錢找tb，也不愿意自己學(xué)習(xí)一下，搞不懂！# install.packages("devtools") devtools::install_github("Hy4m/linkET", force = TRUE) 使用 一般的相關(guān)性分析是用于兩列數(shù)據(jù)之間的，而mantel test 是用于兩個(gè)矩陣的相關(guān)性檢驗(yàn)，我在工作中用的很少，做微生物的小伙伴應(yīng)該用的比較多，做腸道菌群、宏基因組的應(yīng)該也會用到。 加載R包和數(shù)據(jù)：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") 看看這兩個(gè)數(shù)據(jù)長什么樣，因?yàn)閿?shù)就是圖，圖就是數(shù)！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… 接下來先進(jìn)行mantel test檢驗(yàn)，然后畫圖即可。# 計(jì)算 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), # 對相關(guān)系數(shù)進(jìn)行分割，便于映射大小 labels = c("< 0.2", "0.2 - 0.4", ">= 0.4")), pd = cut(p, breaks = c(-Inf, 0.01, 0.05, Inf), # 對P值進(jìn)行分割，便于映射顏色 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), # 這行代碼是關(guān)鍵 data = mantel, curvature = nice_curvature()) + # 下面就是各種顏色和名稱設(shè)置 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))  mantel test可視化 簡單方便，快捷畫圖，出圖效果也很好。 當(dāng)然各種細(xì)節(jié)都是可以用ggplot2語法修改的，喜歡折騰的可以自己嘗試下~ 這個(gè)圖看起來很復(fù)雜，但是展示的信息確實(shí)很多，不過既然你需要這張圖，那我相信你應(yīng)該理解這張圖的意思~ |
|
|
