本文介紹兩個(gè)分析微生物群落和環(huán)境因子相關(guān)性的工具。

Bioenv

Bioenv比較常用，是vegan包中的函數(shù)。 功能：找到與群落差異有最大(秩)相關(guān)性的環(huán)境變量的最佳子集

 1library(vegan)
 2#說(shuō)明
 3bioenv(comm, env, method = "spearman", index = "bray",
 4       upto = ncol(env), trace = FALSE, partial = NULL, 
 5       metric = c("euclidean", "mahalanobis", "manhattan", "gower"),
 6       parallel = getOption("mc.cores"), ...)
 7comm:群落
 8env:環(huán)境因子
 9method:相關(guān)性方法，同cor函數(shù)
10index：群落不相似性矩陣的構(gòu)建方法。同vegdist函數(shù)
11upto：子集中最多元素個(gè)數(shù)
12partial：是否做partial分析
13metric：環(huán)境因子矩陣構(gòu)建方法，有四種方法。

對(duì)于p個(gè)環(huán)境因子，存在2^p-1個(gè)子集，因此環(huán)境因子多的時(shí)候會(huì)非常非常非常慢。20個(gè)變量包含>1e6組合。

 1#例子
 2data(varespec)
 3data(varechem)
 4sol <- bioenv(wisconsin(varespec) ~ log(N) + P + K + Ca + pH + Al, varechem)
 5#結(jié)果給出了不同數(shù)量環(huán)境因子條件下相關(guān)性最高的組合
 6summary(sol)
 7                    size correlation
 8P                      1      0.2513
 9P Al                   2      0.4004
10P Ca Al                3      0.4005
11P Ca pH Al             4      0.3619
12log(N) P Ca pH Al      5      0.3216
13log(N) P K Ca pH Al    6      0.2822

bvStep

bvStep在sinkr包中。功能和Bioenv一樣。在環(huán)境因子很多的時(shí)候表現(xiàn)更好。

 1# install.packages("remotes")
 2# remotes::install_github("marchtaylor/sinkr")
 3library('sinkr')
 4#說(shuō)明
 5bvStep(
 6  fix.mat,
 7  var.mat,
 8  fix.dist.method = "bray",
 9  var.dist.method = "euclidean",
10  scale.fix = FALSE,
11  scale.var = TRUE,
12  max.rho = 0.95,
13  min.delta.rho = 0.001,
14  random.selection = TRUE,
15  prop.selected.var = 0.2,
16  num.restarts = 10,
17  var.always.include = NULL,
18  var.exclude = NULL,
19  output.best = 10
20)
21fix.mat:固定的群落或環(huán)境因子數(shù)據(jù)(一般是群落)
22var.mat:可變的群落或環(huán)境因子數(shù)據(jù)（一般是環(huán)境因子）
23fix.dist.method    :計(jì)算固定數(shù)據(jù)的不相似矩陣（默認(rèn)bray）
24var.dist.metho:計(jì)算可變數(shù)據(jù)的不相似矩陣(默認(rèn)歐氏距離)
25scale.fix:固定數(shù)據(jù)是否中心化和標(biāo)準(zhǔn)化（默認(rèn)否）
26scale.va：可變數(shù)據(jù)是否中心化和標(biāo)準(zhǔn)化（默認(rèn)是）
27max.rho:最大Spearman 相關(guān)性，用于停止計(jì)算，默認(rèn)0.95
28min.delta.rho：Spearman 最小的變化
29random.selection：是否隨機(jī)開(kāi)始
30prop.selected.var：每次重新開(kāi)始包含的變量比例
31num.restarts:重新開(kāi)始次數(shù)，默認(rèn)50
32var.always.include：每次重新開(kāi)始包含的環(huán)境因子
33output.best：結(jié)果返回多少種組合，默認(rèn)10

可以先經(jīng)過(guò)一輪初步篩選，得到相關(guān)性比較高的因子后，再進(jìn)行二次篩選。

 1#例子
 2library(vegan)
 3data(varespec)
 4data(varechem)
 5set.seed(1)
 6res.biobio1 <- bvStep(wisconsin(varespec), wisconsin(varespec), 
 7                      fix.dist.method="bray", var.dist.method="bray",
 8                      scale.fix=FALSE, scale.var=FALSE, 
 9                      max.rho=0.95, min.delta.rho=0.001,
10                      random.selection=TRUE,
11                      prop.selected.var=0.3,
12                      num.restarts=50,
13                      output.best=10,
14                      var.always.include=NULL
15)
16res.biobio1 
17#結(jié)果中包含相關(guān)性最高的組合，和不同因子數(shù)相關(guān)性最高的結(jié)果。
18$order.by.best
19                            var.incl n.var       rho
201     3,4,13,15,24,28,29,30,35,37,39    11 0.8404258
212  3,4,13,15,24,28,29,30,35,37,39,44    12 0.8385896
223       4,13,15,24,28,29,30,35,37,39    10 0.8346131
234       4,15,17,24,26,29,34,35,37,43    10 0.8302690
245     3,4,15,17,24,26,29,34,35,37,43    11 0.8296674
256            3,4,5,15,22,23,29,38,42     9 0.8294197
267              4,5,15,19,25,26,29,38     8 0.8286999
278         3,4,5,15,17,22,23,29,38,42    10 0.8275624
289                 4,5,15,19,26,29,38     7 0.8260670
2910          4,5,14,15,19,25,26,29,38     9 0.8249346
30
31$order.by.i.comb
32                             var.incl n.var       rho
331                                  26     1 0.5186498
342                               23,26     2 0.6100455
353                            15,26,29     3 0.6873777
364                         15,19,26,29     4 0.7428844
375                       4,15,19,26,29     5 0.7760409
386                    4,15,19,26,29,38     6 0.8038093
397                  4,5,15,19,26,29,38     7 0.8260670
408               4,5,15,19,25,26,29,38     8 0.8286999
419             3,4,5,15,22,23,29,38,42     9 0.8294197
4210       4,13,15,24,28,29,30,35,37,39    10 0.8346131
4311     3,4,13,15,24,28,29,30,35,37,39    11 0.8404258
4412  3,4,13,15,24,28,29,30,35,37,39,44    12 0.8385896
4513 2,3,6,7,11,12,13,14,16,35,42,43,44    13 0.7157158
46
47$best.model.vars
48[1] "Rhodtome,Vaccmyrt,Dicrpoly,Pleuschr,Cladstel,Cladgrac,Cladfimb,Cladcris,Cetreric,Flavniva,Stersp"
49
50$best.model.rho
51[1] 0.8404258
52
53$var.always.include
54NULL
55
56$var.exclude
57[1]  8 20 21 32
58
59#23，26，29一直存在，因此第二次篩選將這幾個(gè)樣本包含進(jìn)去。
60set.seed(1)
61res.biobio2  <- bvStep(wisconsin(varespec), wisconsin(varespec), 
62                       fix.dist.method="bray", var.dist.method="bray",
63                       scale.fix=FALSE, scale.var=FALSE, 
64                       max.rho=0.95, min.delta.rho=0.001,
65                       random.selection=TRUE,
66                       prop.selected.var=0.3,
67                       num.restarts=50,
68                       output.best=10,
69                       var.always.include=c(23,26,29)
70)
71res.biobio2 # Best rho equals 0.895 (15 of 44 variables)
72
73#最優(yōu)的組合包含15個(gè)環(huán)境因子
74$best.model.vars
75[1] "Vaccmyrt,Vaccviti,Polypili,Pohlnuta,Cladarbu,Cladrang,Cladcocc,Cladcorn,Cladgrac,Cladfimb,Cladcris,Cladamau,Cetreric,Nepharct,Cladcerv"

sinkr包重點(diǎn)包含了多變量方法和地理數(shù)據(jù)的處理。

sinkr也有做bioenv的函數(shù)為bioEnv。

本文分享自微信公眾號(hào) - Listenlii（gh_1a9e56035563），作者：水岸風(fēng)堤

原文出處及轉(zhuǎn)載信息見(jiàn)文內(nèi)詳細(xì)說(shuō)明，如有侵權(quán)，請(qǐng)聯(lián)系 yunjia_community@tencent.com 刪除。

原始發(fā)表時(shí)間：2020-11-15

本文參與騰訊云自媒體分享計(jì)劃，歡迎正在閱讀的你也加入，一起分享。