R包vegan的相似性分析(ANOSIM)判斷群落結構差異

R包vegan的相似性分析(ANOSIM)判斷群落結構差異
相似性分析(Analysis of similaritiesANOSIM),是一種用於分析高維度資料組間相似性的非引數檢驗方法。它首先透過變數計算物件間距離測度(或者相似度),然後計算關係排名,最後透過排名進行置換檢驗判斷物件組間差異是否顯著不同於組內差異。

ANOSIM示意圖,透過將每個類別內的距離與類別之間的平均距離進行比較,當類內距離較小且類間距離較大時,表明分類明顯。
例如在群落分析中,常使用ANOSIM結合排序分析一起,評估不同環境群落組成結構的相似/差異程度,組間群落差異是否顯著不同於組內差異。
本篇以某16S擴增子測序所得的細菌群落資料為例,展示Rvegan進行ANOSIM檢驗群落結構組成差異的一般過程。
示例資料及R程式碼的百度盤連結:
https://pan.baidu.com/s/11U3kECzMkPriCKl43Wk60w
示例資料概要

示例資料共涉及了4016S測序樣本(細菌群落),均來自土壤。這40個土壤樣本共來自5個不同的取樣地點,每個地點各包含8個樣本(生物學重複)。
此處希望透過ANOSIM,檢視不同取樣地點的土壤細菌群落組成結構是否具有顯著的不同。
檔案“otu_table.txt”為OTU水平的物種丰度表格。每一列為一個樣本,每一行為一種OTU,交叉區域為每種OTU在各樣本中的丰度。

檔案“group.txt”為樣本分組資訊。第一列(names)為各樣本名稱;第二列(site)為各樣本的分組資訊,即這些樣本所屬的取樣地點(s1,地點1s2,地點2s3……以此類推)。

檔案“bray.txt”為提前計算得到的距離矩陣檔案(此處為細菌群落間的Bray-curtis距離)。每一列為一個樣本,每一行為一個樣本,交叉區域為樣本間的Bray-curtis距離(取值範圍0-1,越接近於1表明樣本間細菌群落組成差異越大)。

ANOSIM檢驗群落結構差異
全域性檢驗
首先在所有分組水平上,使用ANOSIM檢驗整體差異,即檢視來自5個不同取樣地點所獲得的土壤細菌群落結構組成在整體上是否不具備一致性,各組間群落差異是否顯著不同於各組內的差異。
vegan中執行ANOSIM的函式為anosim()。輸入anosim()的資料即可以為變數矩陣+分組資訊,也可以為距離矩陣+分組資訊。
library(vegan)
##讀入檔案

#OTU 丰度表

otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

otu <- data.frame(t(otu))
#樣本分組檔案

group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)
##ANOSIM 分析(所有分組間比較,即整體差異)

#根據 group$site 這一列分組進行 ANOSIM 分析檢驗組間差異,基於 999 次置換,詳情 ?anosim
#(1)直接輸入 OTU 丰度表,在引數中指定距離型別

#使用 Bray-Curtis 距離測度

anosim_result <- anosim(otu, group$site, distance = 'bray', permutations = 999)
#(2)或者首先根據丰度表計算群落間距離,仍然使用 Bray-Curtis 距離測度,詳情 ?vegdist

dis1 <- vegdist(otu, method = 'bray')
#再將所的距離資料作為輸入

anosim_result <- anosim(dis1, group$site, permutations = 999)
#(3)若是已經提供好了距離矩陣,則直接使用現有的距離矩陣進行分析即可

#現有的距離矩陣,這裡為 Bray-Curtis 距離測度

dis <- read.delim('bray.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

dis1 <- as.dist(dis)

anosim_result <- anosim(dis1, group$site, permutations = 999)
特別注意,為了避免識別錯誤,不要使用純數字作為分組名稱。對於置換次數,視資料量而定。
以上展示了3種處理過程,由於所使用距離型別相同(均為Bray-Curtis距離),分組資訊、置換檢驗次數(均為999次)等也一致,若置換檢驗次數已達到飽和,則結果也是一致的。
#檢視結果,上述 3 條命令所計算的內容一致,以其中一個為例

summary(anosim_result)

#或者

names(anosim_result)

anosim_result$signif #p 值

anosim_result$statistic #R 值

summary()後螢幕輸出了包含執行命令列、置換檢驗次數、差異計算結果等資訊。
其中,可主要關注兩個重要統計值,R值(ANOSIM statistic R)和p值(Significance)。R值可以得出組間與組內比較的差異程度,其取值範圍(-11);R>0,說明組間差異大於組內差異,即組間差異顯著;R<0,說明組內差異大於組間差異;R值的絕對值越大表明相對差異越大。p值越低表明這種差異檢驗結果越顯著,一般以0.05為顯著性水平界限。

使用plot()函式能夠將計算結果視覺化。
#作圖展示

#pdf(paste('anosim.all.pdf', sep = ''), width = 10, height = 5)

png(paste('anosim.all.png', sep = ''), width = 800, height = 400)

plot(anosim_result, col = c('gray', 'red', 'green', 'blue', 'orange', 'purple'))

dev.off()

ANOSIM分析基於兩兩樣本之間的距離值排序獲得的秩,這樣任一兩兩組的比較可以獲得三個分類的資料。以箱線圖的形式展示組間與組內的秩的分佈,橫座標表示所有樣品Between)以及各分組(本示例為s1s2s3s4s5),縱座標表示距離(本示例使用Bray-Curtis距離)的秩。當Between組相對與其他每個分組的秩較高時,則表明組間差異大於組內差異。同時圖的上方標註了R值與p值兩個重要統計指標,便於我們直觀地對組間差異是否顯著不同於各組內的差異進行判斷。

據此,我們可知5個取樣地點的土壤細菌群落結構在整體水平上是不一致的,各組間群落差異是否顯著不同於各組內的差異。
組間兩兩比較
可進一步使用ANOSIM評估兩兩環境(兩兩取樣地)之間的土壤細菌群落結構組成的差異,包括是否具有顯著性,以及組間和組內的變異程度等。
##ANOSIM 分析(使用迴圈處理,進行小分組間比較,如兩組間)

#推薦使用 OTU 丰度表作為輸入資料,每次篩選分組後重新計算樣本距離,避免由於樣本數減少可能導致的距離變動而造成誤差

group_name <- unique(group$site)
dir.create('anosim_two', recursive = TRUE)

anosim_result_two <- NULL

for (i in 1:(length(group_name) – 1)) {

for (j in (i + 1):length(group_name)) {

group_ij <- subset(group, site %in% c(group_name[i], group_name[j]))

otu_ij <- otu[group_ij$names, ]

anosim_result_otu_ij <- anosim(otu_ij, group_ij$site, permutations = 999, distance = 'bray') #Bray-Curtis 距離測度,基於 999 次置換

anosim_result_two <- rbind(anosim_result_two, c(paste(group_name[i], group_name[j], sep = '/'), 'Bray-Curtis', anosim_result_otu_ij$statistic, anosim_result_otu_ij$signif))
#每次迴圈輸出圖片

#pdf(paste('anosim_two/anosim.', group_name[i], '_', group_name[j], '.pdf', sep = ''), width = 7, height = 5)

png(paste('anosim_two/anosim.', group_name[i], '_', group_name[j], '.png', sep = ''), width = 600, height = 400)

plot(anosim_result_otu_ij, col = c('gray', 'red', 'blue'))

dev.off()

}

}
#帶 R 值和 p 值的表格

anosim_result_two <- data.frame(anosim_result_two, stringsAsFactors = FALSE)

names(anosim_result_two) <- c('group', 'distance', 'R', 'P_value')
#可選新增 p 值校正過程,例如 Benjamini 校正

anosim_result_two$P_value <- as.numeric(anosim_result_two$P_value)

anosim_result_two$P_adj_BH <- p.adjust(anosim_result_two$P_value, method = 'BH')
write.table(anosim_result_two, 'anosim_two/ANOSIM.result_two.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')

以上示例,在每一步迴圈中挑選出兩個分組,並進行兩兩分組間的ANOSIM。輸出的ANOSIM統計總表中,包含了計算所依據的距離型別(Distance)、ANOSIM分析的R值(R)、顯著性p值(P_value)以及BH校正後p值(P_adj_BH)等多項內容。

每步迴圈的ANOSIM箱線圖儲存至“anosim_two”中。

其餘的圖片結果,如“anosim.s1_s2.png”,透過該圖我們可知s1s2取樣地間的土壤細菌群落組成無較大差異,或者說組內差異大於了組間差異。具體表現在Between組相對與其他兩個分組的秩較低,以及p值不顯著(0.122,高於0.05水平)。

類似地,透過“anosim.s2_s4.png”,我們可知s2s4取樣地間的土壤細菌群落組成具有顯著差異,即兩組間群落差異大於了各組內的差異。具體表現在Between組的秩明顯高於其他兩個分組的秩,R值較高(0.586),且p值很顯著(0.001水平,不過這僅為校正前的p值,可再結合統計表檢視校正後的p值是否同樣顯著)。

即可獲知來自5個取樣地點的土壤細菌群落中,兩兩之間的相似/相異水平,以及組間和組內的變異程度,並幫助尋找在排序圖中肉眼觀察不明顯的組間差異。
連結

相關文章