導(dǎo)讀

多樣性分析是16S測(cè)序分析中常見的分析方法。本文旨在向初學(xué)者介紹多樣性分析中alpha多樣性和beta多樣性的由來(lái)、概念、計(jì)算，以及展示“如何用R語(yǔ)言實(shí)現(xiàn)alpha多樣性和beta多樣性的計(jì)算和可視化”的工作流程。

一、概念介紹

Alpha多樣性、beta多樣性和gamma多樣性的概念由美國(guó)杰出的生態(tài)學(xué)家Robert Harding Whittaker提出。Whittaker將一個(gè)環(huán)境中總的物種多樣性命名為gamma多樣性。Gamma多樣性由alpha多樣性和beta多樣性共同決定。Alpha多樣性指在環(huán)境中的一個(gè)區(qū)域內(nèi)的平均物種多樣性。Beta多樣性指環(huán)境中不同區(qū)域之間的差異。

二、Alpha多樣性指數(shù)

1. Shannon指數(shù)：

Shannon指數(shù)最初是由Claude Shannon提出用來(lái)計(jì)算字符串文本信息熵的指標(biāo)，后來(lái)逐漸發(fā)展成生態(tài)學(xué)研究中最常用的多樣性指數(shù)。Shannon指數(shù)不只考慮物種豐富度（Richness，樣本中物種數(shù)），而且同時(shí)考慮物種的均勻度（Evenness，平均程度），所以它是反應(yīng)群落結(jié)構(gòu)的綜合指標(biāo)。
計(jì)算方法：

公式

Pi：第i個(gè)物種的個(gè)體數(shù)占樣本中總個(gè)體數(shù)的比例
R：樣本中的總物種數(shù)
樣本的Shannon指數(shù)越高，則其物種多樣性也越高。

2. Simpson指數(shù)：

Simpson指數(shù)是由Edward H. Simpson在1949年提出評(píng)價(jià)生物多樣性的另一種常用指標(biāo)。它也是既考慮樣本中的物種豐富度，又考慮物種均勻度的綜合指標(biāo)。
計(jì)算方法：

公式

Pi：第i個(gè)物種的個(gè)體數(shù)占樣本中總個(gè)體數(shù)的比例
R：樣本中的總物種數(shù)
Simpson指數(shù)的取值范圍在0和1之間。樣本的Simpson指數(shù)越高，則其物種多樣性也越高。

三、Beta多樣性的計(jì)算

Beta多樣性專用于不同樣本間的比較，它不能直接通過某一個(gè)樣本的物種豐富度和均勻度計(jì)算出該樣本的多樣性度量值。Beta多樣性是利用不同樣本間的豐度變化或進(jìn)化關(guān)系來(lái)計(jì)算樣本間距離，從而反映樣本間是否具有顯著的微生物群落差異。
計(jì)算beta多樣性的方法有很多：有最為常用的bray curtis距離、Jaccard距離還有歐式距離，他們考慮的是樣本間物種豐度（有無(wú)）和均度（相對(duì)豐度）；另外還有Unifrac距離法，它是根據(jù)系統(tǒng)發(fā)生樹進(jìn)行比較。

四、準(zhǔn)備工作

軟件準(zhǔn)備（window環(huán)境）：

1. R
地址：https://www.r-project.org/
版本：3.4.1

2. R包：openxlsx
功能：打開Excel文件
版本：4.1

3. R包：vegan
功能：包內(nèi)含有多種生態(tài)學(xué)分析必備的函數(shù)
版本：2.4

4. Excel
功能：數(shù)據(jù)處理

數(shù)據(jù)準(zhǔn)備：

標(biāo)準(zhǔn)化菌屬相對(duì)豐度表。獲取方法請(qǐng)參考16S測(cè)序分析系列（一）菌屬豐度表獲取。

五、Alpha多樣性的計(jì)算和可視化

install.packages("openxlsx")  ##  下載安裝openxlsx包
install.packages("vegan")  ##  下載安裝vegan包

setwd("C:/mywd")  ##  設(shè)置工作目錄為C盤mywd文件夾
getwd()  ##  進(jìn)入工作目錄

library(openxlsx)  ##  調(diào)用openxlsx包
library(vegan)  ##  調(diào)用vegan包

workbook <- "C:/mywd/data.xlsx"
##  讀取并存儲(chǔ)data.xlsx的內(nèi)容到workbook

說明：data.xlsx必須是標(biāo)準(zhǔn)化后的菌屬相對(duì)豐度表，格式如下：第一行是菌屬名稱，其余每行代表一個(gè)樣本，每列代表一類菌屬的在每一個(gè)樣本中的相對(duì)豐度。

圖片.png

mydataframe <- read.xlsx(workbook, 1)
##  讀取并存儲(chǔ)workbook中的第一個(gè)sheet的內(nèi)容到mydataframe

alpha_shannon <- diversity(mydataframe, "shannon") 
##  用vegan包中的diversity函數(shù)計(jì)算shannon多樣性

write.xlsx(alpha_shannon, "data_diversity.xlsx")
## 獲得shannon多樣性指標(biāo)的計(jì)算結(jié)果，即data_diversity.xlsx
## 保存alpha_shannon到當(dāng)前的data_diversity.xlsx中（可自動(dòng)創(chuàng)建）

數(shù)據(jù)處理：打開data_diversity.xlsx文件如“1”。其中第一個(gè)樣本的shannon index是shannon index 1，以此類推。然后我們把這些shannon index加上變量名和分組信息（將其分為A，B兩組），得到“2”。利用R或者Excel可以對(duì)“2”中的數(shù)據(jù)進(jìn)行分組T檢驗(yàn)，方法很簡(jiǎn)單，所以這里不再加以展示。接下來(lái)我將利用“2”制作boxplot來(lái)更直觀的觀察兩組的shannon index是否有明顯的差異。

圖片.png

workbook2 <- "C:/mywd/data_diversity.xlsx"
##  讀取并存儲(chǔ)data_diversity.xlsx的內(nèi)容到workbook2

mydataframe2 <- read.xlsx(workbook2, 1)
##  讀取并存儲(chǔ)workbook2中的第一個(gè)sheet的內(nèi)容到mydataframe2

boxplot(mydataframe2$Diversity~ mydataframe2$Group, main="Shannon index", col=c("green", "red"))
## mian：命名；col：上色
## 可視化結(jié)果如下。

圖片.png

由圖可見，A組的多樣性明顯高于B組。

六、Beta多樣性的計(jì)算和可視化

setwd("C:/mywd")
getwd()

library(openxlsx)
library(vegan)

workbook <- "C:/mywd/data.xlsx"
mydataframe <- read.xlsx(workbook, 1)
ord <- cmdscale(vegdist(mydataframe, method="bray"))

write.xlsx(ord, "data_ord.xlsx")

數(shù)據(jù)處理：data_ord.xlsx中保存的是每個(gè)樣本的二維坐標(biāo)信息，它是進(jìn)行此分析的最關(guān)鍵的文件。用windows打開data_ord.xlsx文件，格式如“1”。第一行是變量名，其余每行是每個(gè)樣本的坐標(biāo)信息。添加分組信息，如下圖“2”。為了方便作圖，我這里直接使用顏色”Color”代替組別”Group”，即red代替Group A，blue代替Group B。接下來(lái)利用“2”中的樣本坐標(biāo)信息進(jìn)行繪圖。

圖片.png

workbook3 <- "C:/mywd/data_ord.xlsx"
mydataframe3 <- read.xlsx(workbook3, 1)

plot(x=mydataframe3$V1,y=mydataframe3$V2,xlab="mds1",ylab="mds2",main="Beta diversity ( Bray Curtis )" , pch=19, col=mydataframe3$Color)
##  x=mydataframe3$V1,y=mydataframe3$V2：以V1為橫坐標(biāo)，V2為縱坐標(biāo)
##  xlab：設(shè)置坐標(biāo)名稱
##  main：設(shè)置圖片名稱
##  pch=19：用點(diǎn)表示樣本
##  col=mydataframe3$Color：根據(jù)Color中的信息對(duì)點(diǎn)上色

legend("topright", pch=c(19, 19), col=c("red", "blue"), legend=c("Group A", "Group B"))
##  "topright"：設(shè)置圖注在右上角
##  pch=c(19, 19)：繪制兩個(gè)點(diǎn)代表兩個(gè)組
##  col：給點(diǎn)上色
##  legend：給點(diǎn)命名

圖片.png

結(jié)束語(yǔ)

到這里就完成了beta多樣性的計(jì)算和可視化的所有工作。從圖中可以看出A和B兩個(gè)組的樣本的腸道菌群有明顯的差異，所以可以推測(cè)“分組變量”很可能和腸道菌群有聯(lián)系。要想知道具體是“哪些菌”對(duì)此聯(lián)系做了貢獻(xiàn)，還需要進(jìn)行更深層次的分析。本文到此就結(jié)束了，如有疑問歡迎留言交流。下次將為大家推出“如何從腸道菌群中尋找你感興趣的細(xì)菌？”來(lái)介紹尋找“哪些菌”的方法。

相關(guān)閱讀：
16S測(cè)序分析（一）菌屬豐度表獲取
16S測(cè)序分析（二）菌群多樣性分析
16S測(cè)序分析（三）用LEfSe尋找組間差異細(xì)菌
16S測(cè)序分析（四）用MaAsLin尋找組間差異細(xì)菌
16S測(cè)序分析（五）用RandomForest尋找關(guān)鍵細(xì)菌
16S測(cè)序分析（六）用PICRUSt預(yù)測(cè)菌群KEGG代謝通路