alpha多樣性是進行生物群落研究的最基礎方法,關於alpha多樣性指數的介紹和計算在網上有茫茫多的資料,本文將不會花費大量篇幅介紹這些基礎知識,而是重點帶來對alpha多樣性進行統計學分析並直接生成圖像的方法。
Alpha多樣性指數Alpha多樣性用於分析樣品內(Within-community)的微生物群落多樣性,可以反映樣品內的微生物群落的豐富度和多樣性。
alpha多樣性指數包括豐富度、多樣性、均一性等。
豐富度指數只是簡單的估算群落中物種的數量,而不考慮群落中每個物種的豐度情況,例如Chao1、ACE等。
多樣性指數同時考慮的物種在群落中是否存在即其所佔的豐度,例如Shannon、Simpson等。
均一性指數用於評估群落中各個微生物在群落中的均勻程度,所有微生物均具有相同的豐度則均一性最高,例如Pielou等。
Alpha多樣性指數計算目前可以用來計算alpha多樣性指數的工具非常的多,可以說基本上所有的微生物群落分析工具或相關的統計學軟體都包含了alpha多樣性指數的計算功能,最不濟我們還可以根據每個指數的公式進行手動計算。
這裡推薦兩種alpha多樣性指數的計算工具:
再給大家安利一下PAST這個軟體,真是可以說是最輕便、學習成本最低、功能最多樣的統計學分析軟體了,之前做了一套PAST的完整使用教程,大家可以看一下,可以說是瞬間就可以學會。
Alpha多樣性差異檢驗在微生物群落的alpha多樣性指數分析中,最常用的就是利用統計學分析檢驗不同組樣本間微生物群落alpha多樣性指數的差異顯著性。
兩組樣本分析當研究的樣本只有兩組時,一般使用t-test檢驗組間差異。
數據格式為兩列,第一列為各樣本的alpha多樣性指數,第二列為樣本的分組信息。
首先導入分析數據:
alpha <- read.table("alpha_diversity.txt",header = TRUE,row.names = 1,sep = "\t")手動定義下分組名稱的順序,以確定畫圖時各組的位置,默認畫圖時時按照組名的字母順序排列。
alpha$Group <- factor(alpha$Group,levels = c("Water","Sediment"))使用t.test()函數進行t-test,並將結果保存。
t <- t.test(Shannon~Group,data = alpha)write.csv(t$p.value,file = "t.test.shannon.txt")進行結果圖的繪製,結果以箱須圖展示並根據t-test的p-value標註組間差異性。
png(filename = "Shannon.png",width = 5400,height = 7200,res = 600)par(mar = c(5,8,2,2))par(xpd = TRUE)box <- boxplot(Shannon~Group,data = alpha,col = "grey80",lwd = 4,range = 0,axes = FALSE,ylim = c(5,11),xlab = "",ylab = "")axis(side = 1,at = c(1,2),labels = c("",""),lwd = 4)segments(x0 = 0.42,y0 = 4.76,x1 = 2.5,y1 = 4.76,lwd = 5)text(x = c(1,2),y = 4.52,labels = c("water","Sediment"),cex = 3,font = 2)axis(side = 2,at = c(5,6,7,8,9,10),labels = c("","","","","",""),lwd = 4)text(x = 0.37,y = c(5,6,7,8,9,10),labels = c("5","6","7","8","9","10"),adj = c(1,0.5),cex = 2.5)segments(x0 = 0.42,y0 = 4.76,x1 = 0.42,y1 = 11,lwd = 5)mtext(side = 2,text = "Shannon index",cex = 3.5,line = 4.5)segments(x0 = 1,y0 = 7.3,x1 = 1,y1 = 10.6,lwd = 5)segments(x0 = 1,y0 = 10.6,x1 = 2,y1 = 10.6,lwd = 5)segments(x0 = 2,y0 = 10.6,x1 = 2,y1 = 10,lwd = 5)text(x = 1.5,y = 10.8,labels = ifelse(t$p.value <= 0.001, "***",ifelse(t$p.value <= 0.01, "**", ifelse(t$p.value <= 0.05 ,"*"," "))),cex = 6)jitter <- jitter(c(rep(1,36),rep(2,36)))points(x = jitter,y = alpha$Shannon,col = "black",bg = "white",cex = 3,lwd= 3,pch = 21)dev.off()代碼中各參數的意義詳見之前的推文:每天學習一點R:19.箱須圖的繪製!!
多組樣本分析大多數情況下,我們分析的樣本可能並不是兩組而是多組,此時需要使用ANVOA來檢驗樣本間多樣性指數的差異,但是ANOVA只是用於評估所有樣本整體是否具有差異,各組樣本間的兩兩比較還需要進一步使用Turkey HSD test來進行檢驗。
當樣本組數目較多時,如果使用兩兩之間連線+顯著性標註的方法來展示樣本簡單的差異會比較亂,因而通常是使用相同的字母來代表組間沒有差異的樣本,而常規的Turkey test檢驗結果只會給出p-value而並不會標出各組之間的分類結果,下面就來介紹一個基於R語言自動進行統計學檢驗並以字母標出組間差異的方法。
分析的數據格式與兩組樣本一致,均為兩列,一列為alpha多樣性指數,另一列為分組信息。
首先還是要導入數據,並手動指定分組因子的順序。
Sediment <- read.table("Sediment.txt",header = T,sep = "\t",row.names = 1)Sediment$Group <- factor(Sediment$Group,levels = c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))接下來進行Turkey檢驗並保存結果。
fit <- aov(chao1~Group,data = Sediment)Tukey <- TukeyHSD(fit)write.csv(Tukey$Group,file = "Tukey.test.S.chao1.csv")接來下進行圖像繪製,繪圖需要依賴multcomp包。
library(multcomp)png(filename = "S.chao1.png",width = 5400,height = 5400,res = 700)par(mar=c(4,6,6,2)) tuk<-glht(fit,linfct=mcp(Group="Tukey"))p1 <- plot(cld(tuk,level=0.05),col="lightgrey",cex = 4)dev.off()當然由於這項分析的主要目的是給出統計學結果,因而對於結果可視化上就沒有進行優化,結果圖像中的各項參數相對比較難調整,建議應用該方法得到樣本組分類對應的字母編號,之後在繪製合適的箱須圖後手動添加這些字母標註。
關注公眾號「紅皇后學術」,後臺回復「alpha」獲取本文介紹的繪圖示例文件和完整代碼。