前兩天看中科院微生物所高程老師的文章：Fungal community assembly in drought-stressed sorghum shows stochasticity, selection, and universal ecological dynamics，里面有這樣一個圖：
仔細看了材料與方法部分才知道這是臨界指示物種分析法（Threshold Indicator Taxa ANalysis肯腕，TITAN）枢冤。這種方法經(jīng)常用于分析群落在時間梯度或者是環(huán)境梯度上變化的閾值。如上圖表示的是隨著時間變化黍瞧，不同真菌群落的閾值變化躬翁。該方法于 2010 年被開發(fā)出來：
加載示例數(shù)據(jù)
先加載 R 包和示例數(shù)據(jù)：
rm(list = ls())
library(TITAN2)
data("glades.taxa")
data("glades.env")
這個方法需要的數(shù)據(jù)有兩個：
-
群落數(shù)據(jù)：示例數(shù)據(jù)是通過
log10(x+1)進行轉(zhuǎn)換過的昧旨，開發(fā)者認為通常不用進行轉(zhuǎn)換拜轨。進行轉(zhuǎn)換的目的是避免稀有種對群落的影響抽减。群落數(shù)據(jù)的行為樣本，列為群落豐度數(shù)據(jù)橄碾。
image 環(huán)境數(shù)據(jù)：可以是時間梯度卵沉，也可以是環(huán)境梯度。高程老師文章使用的是時間梯度堪嫂，在示例數(shù)據(jù)中使用的是總磷含量作為環(huán)境梯度：
參數(shù)一覽
默認參數(shù)：
res = titan(glades.env,glades.taxa)
完整參數(shù)：
res = titan(glades.env,
glades.taxa,
minSplt = 5,
numPerm = 999,
boot = TRUE,
nBoot = 500,
imax = FALSE,
ivTot = FALSE,
pur.cut = 0.95,
rel.cut = 0.95,
ncpus = 1,
memory = FALSE)
參數(shù)解釋：
-
glades.env：環(huán)境數(shù)據(jù)偎箫； -
glades.taxa：群落數(shù)據(jù)木柬； -
minSplt：計算在突變點任一側(cè)相關(guān)統(tǒng)計量的最小觀察數(shù)皆串，不能小于 3； -
numPerm：隨機置換檢驗的次數(shù)眉枕，不能小于 250恶复，通常是 500 或 1000怜森； -
boot：是否自舉重抽樣； -
nBoot：自舉重抽樣的得到的新數(shù)據(jù)集的數(shù)據(jù)量谤牡，小于 500 能夠有效減少計算資源的消耗副硅，但通常是 500 或 1000； -
imax：是否用IndVal最大值或者 z-score 來決定每個群落的閾值翅萤； -
ivTot：是否通過平均相對豐度來計算IndVal恐疲； -
pur.cut：純響應(yīng)方向的閾值； -
rel.cut：可靠響應(yīng)幅度的閾值套么； -
ncpu：自舉重抽樣部分的處理的核數(shù)培己； -
memory：是否在自舉重抽樣的時候建立臨時文件夾存放文件。
運行耗時
樣品數(shù)量多少胚泌、置換檢驗次數(shù)省咨、自舉重抽樣次數(shù)等幾個因素決定了運行時長，從幾分鐘到幾個小時不等玷室。示例數(shù)據(jù)有 126 個樣本零蓉，164 個群落數(shù)據(jù)，使用下方的代碼計算運行時長：
system.time(titan(glades.env,
glades.taxa,
minSplt = 5,
numPerm = 250,
boot = TRUE,
nBoot = 500,
imax = FALSE,
ivTot = FALSE,
pur.cut = 0.95,
rel.cut = 0.95,
ncpus = 1,
memory = FALSE)
)
可以看到穷缤，耗時在 21 分鐘左右敌蜂。
結(jié)果可視化
老版本的可視化是這樣的（圖中的點代表的是臨界點（突變點），點的大小代表響應(yīng)強度的大猩鹣睢）：
代碼如下：
plot_taxa(glades.titan, xlab = 'Surface Water TP (ug/l)')
新的可視化方法是這樣的：
- 總的：
plot_taxa_ridges(glades.titan,
xlab = expression(paste('Surface water total phosphorus ('*mu*'g/l)')))
-
Decreasers：
plot_taxa_ridges(glades.titan,
z2 = FALSE,
xlab = expression(paste('Surface water total phosphorus ('*mu*'g/l)')))
-
Increasers：
plot_taxa_ridges(glades.titan,
z1 = FALSE,
xlab = expression(paste('Surface water total phosphorus ('*mu*'g/l)')))
