導讀
多樣性分析是16S測序分析中常見的分析方法。本文旨在向初學者介紹多樣性分析中alpha多樣性和beta多樣性的由來、概念、計算，以及展示“如何用R語言實現(xiàn)alpha多樣性和beta多樣性的計算和可視化”的工作流程。
一、概念介紹
Alpha多樣性、beta多樣性和gamma多樣性的概念由美國杰出的生態(tài)學家Robert Harding Whittaker提出。Whittaker將一個環(huán)境中總的物種多樣性命名為gamma多樣性。Gamma多樣性由alpha多樣性和beta多樣性共同決定。Alpha多樣性指在環(huán)境中的一個區(qū)域內(nèi)的平均物種多樣性。Beta多樣性指環(huán)境中不同區(qū)域之間的差異。
二、Alpha多樣性指數(shù)
1. Shannon指數(shù)：
Shannon指數(shù)最初是由Claude Shannon提出用來計算字符串文本信息熵的指標，后來逐漸發(fā)展成生態(tài)學研究中最常用的多樣性指數(shù)。Shannon指數(shù)不只考慮物種豐富度（Richness，樣本中物種數(shù)），而且同時考慮物種的均勻度（Evenness，平均程度），所以它是反應群落結(jié)構(gòu)的綜合指標。
計算方法：
Pi：第i個物種的個體數(shù)占樣本中總個體數(shù)的比例
R：樣本中的總物種數(shù)
樣本的Shannon指數(shù)越高，則其物種多樣性也越高。
2. Simpson指數(shù)：
Simpson指數(shù)是由Edward H. Simpson在1949年提出評價生物多樣性的另一種常用指標。它也是既考慮樣本中的物種豐富度，又考慮物種均勻度的綜合指標。
計算方法：
Pi：第i個物種的個體數(shù)占樣本中總個體數(shù)的比例
R：樣本中的總物種數(shù)
Simpson指數(shù)的取值范圍在0和1之間。樣本的Simpson指數(shù)越高，則其物種多樣性也越高。
三、Beta多樣性的計算
Beta多樣性專用于不同樣本間的比較，它不能直接通過某一個樣本的物種豐富度和均勻度計算出該樣本的多樣性度量值。Beta多樣性是利用不同樣本間的豐度變化或進化關(guān)系來計算樣本間距離，從而反映樣本間是否具有顯著的微生物群落差異。
計算beta多樣性的方法有很多：有最為常用的bray curtis距離、Jaccard距離還有歐式距離，他們考慮的是樣本間物種豐度（有無）和均度（相對豐度）；另外還有Unifrac距離法，它是根據(jù)系統(tǒng)發(fā)生樹進行比較。
四、準備工作
軟件準備（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)學分析必備的函數(shù)
版本：2.4
4. Excel
功能：數(shù)據(jù)處理
數(shù)據(jù)準備：
標準化菌屬相對豐度表。獲取方法請參考16S測序分析系列（一）菌屬豐度表獲取。
五、Alpha多樣性的計算和可視化
install.packages("openxlsx") ## 下載安裝openxlsx包
install.packages("vegan") ## 下載安裝vegan包
setwd("C:/mywd") ## 設置工作目錄為C盤mywd文件夾
getwd() ## 進入工作目錄
library(openxlsx) ## 調(diào)用openxlsx包
library(vegan) ## 調(diào)用vegan包
workbook <- "C:/mywd/data.xlsx"
## 讀取并存儲data.xlsx的內(nèi)容到workbook
說明：data.xlsx必須是標準化后的菌屬相對豐度表，格式如下：第一行是菌屬名稱，其余每行代表一個樣本，每列代表一類菌屬的在每一個樣本中的相對豐度。
mydataframe <- read.xlsx(workbook, 1)
## 讀取并存儲workbook中的第一個sheet的內(nèi)容到mydataframe
alpha_shannon <- diversity(mydataframe, "shannon")
## 用vegan包中的diversity函數(shù)計算shannon多樣性
write.xlsx(alpha_shannon, "data_diversity.xlsx")
## 獲得shannon多樣性指標的計算結(jié)果，即data_diversity.xlsx
## 保存alpha_shannon到當前的data_diversity.xlsx中（可自動創(chuàng)建）
數(shù)據(jù)處理：打開data_diversity.xlsx文件如“1”。其中第一個樣本的shannon index是shannon index 1，以此類推。然后我們把這些shannon index加上變量名和分組信息（將其分為A，B兩組），得到“2”。利用R或者Excel可以對“2”中的數(shù)據(jù)進行分組T檢驗，方法很簡單，所以這里不再加以展示。接下來我將利用“2”制作boxplot來更直觀的觀察兩組的shannon index是否有明顯的差異。
workbook2 <- "C:/mywd/data_diversity.xlsx"
## 讀取并存儲data_diversity.xlsx的內(nèi)容到workbook2
mydataframe2 <- read.xlsx(workbook2, 1)
## 讀取并存儲workbook2中的第一個sheet的內(nèi)容到mydataframe2
boxplot(mydataframe2$Diversity~ mydataframe2$Group, main="Shannon index", col=c("green", "red"))
## mian：命名；col：上色
## 可視化結(jié)果如下。
由圖可見，A組的多樣性明顯高于B組。
六、Beta多樣性的計算和可視化
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中保存的是每個樣本的二維坐標信息，它是進行此分析的最關(guān)鍵的文件。用windows打開data_ord.xlsx文件，格式如“1”。第一行是變量名，其余每行是每個樣本的坐標信息。添加分組信息，如下圖“2”。為了方便作圖，我這里直接使用顏色”Color”代替組別”Group”，即red代替Group A，blue代替Group B。接下來利用“2”中的樣本坐標信息進行繪圖。
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為橫坐標，V2為縱坐標
## xlab：設置坐標名稱
## main：設置圖片名稱
## pch=19：用點表示樣本
## col=mydataframe3$Color：根據(jù)Color中的信息對點上色
legend("topright", pch=c(19, 19), col=c("red", "blue"), legend=c("Group A", "Group B"))
## "topright"：設置圖注在右上角
## pch=c(19, 19)：繪制兩個點代表兩個組
## col：給點上色
## legend：給點命名
結(jié)束語
到這里就完成了beta多樣性的計算和可視化的所有工作。從圖中可以看出A和B兩個組的樣本的腸道菌群有明顯的差異，所以可以推測“分組變量”很可能和腸道菌群有聯(lián)系。要想知道具體是“哪些菌”對此聯(lián)系做了貢獻，還需要進行更深層次的分析。本文到此就結(jié)束了，如有疑問歡迎留言交流。下次將為大家推出“如何從腸道菌群中尋找你感興趣的細菌？”來介紹尋找“哪些菌”的方法。
相關(guān)閱讀：
16S測序分析（一）菌屬豐度表獲取
16S測序分析（二）菌群多樣性分析
16S測序分析（三）用LEfSe尋找組間差異細菌
16S測序分析（四）用MaAsLin尋找組間差異細菌
16S測序分析（五）用RandomForest尋找關(guān)鍵細菌
16S測序分析（六）用PICRUSt預測菌群KEGG代謝通路
