微生物擴增子測序圖表解讀(實例數據)

16s分析結果詳解

1.    OUT（是我們要搞清的一個重要概念，可以說是後續分析的基石）

OTU(operational taxonomic units) 是在系統發生學研究或群體遺傳學研究中，為了便於進行分析，人為給某一個分類單元（品系，種，屬，分組等）設置的同一標誌。通常按照 97% 的相似性閾值將序列劃分為不同的 OTU，每一個 OTU 通常被視為一個微生物物種。相似性小於97%就可以認為屬於不同的種，相似性小於93%-95%，可以認為屬於不同的屬。樣品中的微生物多樣性和不同微生物的豐度都是基於對OTU的分析。

有了OTU這個概念之後，就不難理解下表。對每個樣本的測序數量和OTU數目進行統計，並且在表栺中列出了測序覆蓋的完整度。

其中 SampleName表示樣本名稱；SampleSize表示樣本序列總數；OTUsNumber表示注釋上的OTU數目；OTUsSeq表示注釋上OTU的樣本序列總數。

2.    Coverage

Coverage是指各樣品文庫的覆蓋率，其數值越高，則樣本中序列沒有被測出的概率越低。該指數實際反映了本次測序結果是否代表樣本的真實情況。計算公式為：C=1-n1/N 其中n1 = 只含有一條序列的OTU的數目；N = 抽樣中出現的總的序列數目。

下表是對每個樣本在分類字水平上的數量進行統計，並且在表栺中列出了在每個分類字水平上的物種數目

其中SampleName表示樣本名稱；Phylum表示分類到門的OTU數量；Class表示分類到綱的OTU數量；Order表示分類到目的OTU數量；Family表示分類到科的OTU數量；Genus表示分類到屬的OTU數量；Species表示分類到種的OTU數量。

我們可以看到絕大部分的OTU都分類到了屬（Genus），也有很多分類到了種（Species）。但是仍然有很多無法完全分類到種一級，這是由於環境微生物本身存在非常豐富的多樣性，還有大量的菌仍然沒有被測序和發現。

當然，對這些種屬的構成還可以進行柱狀圖展示：

橫坐標中每一個條形圖代表一個樣本，縱坐標代表該分類層級的序列數目或比例。同一種顏色代表相同的分類級別。圖中的每根柱子中的顏色表示該樣本在不同級別（門、綱、目等）的序列數目，序列數目只計算級別最低的分類，例如在屬中計算過了，則在科中則不重複計算。

3.    韋恩圖

韋恩圖是對樣本之間或分組之間的OTU進行比較獲得。

4.    樣品構成豐度

4.1稀釋曲線

微生物多樣性分析中如何驗證測序數據量是否足以反映樣品中的物種多樣性？

稀釋曲線（豐富度曲線）可以派上用場。它是用來評價測序量是否足以覆蓋所有類群，並間接反映樣品中物種的豐富程度。

不免有同學有疑惑，稀釋曲線怎麼來的？

它是利用已測得16S rDNA序列中已知的各種OTU的相對比例，來計算抽取n個（n小於測得reads序列總數）reads時出現OTU數量的期望值，然後根據一組n值（一般為一組小於總序列數的等差數列）與其相對應的OTU數量的期望值做出曲線來。

至此，我們雖然知道了稀釋曲線的由來，那麼這個五彩繽紛的稀釋曲線該怎麼看呢？

當曲線趨於平緩或者達到平臺期時也就可以認為測序深度已經基本覆蓋到樣品中所有的物種，增加測序數據無法再找到更多的OTU；

反之，則表示樣品中物種多樣性較高，還存在較多未被測序檢測到的物種。

橫坐標代表隨機抽取的序列數量；縱坐標代表觀測到的OTU數量。樣本曲線的延伸終點的橫坐標位置為該樣本的測序數量。

4.2 Shannon-Winner曲線

Shannon-Wiener 曲線，是利用shannon指數來進行繪製的，反映樣品中微生物多樣性的指數，利用各樣品的測序量在不同測序深度時的微生物多樣性指數構建曲線，以此反映各樣本在不同測序數量時的微生物多樣性。

當曲線趨向平坦時，說明測序數據量足夠大，可以反映樣品中絕大多數的微生物物種信息。

橫坐標代表隨機抽取的序列數量；縱坐標代表的是反映物種多樣性的Shannon指數，樣本曲線的延伸終點的橫坐標位置為該樣本的測序數量。

其中曲線的最高點也就是該樣本的Shannon指數，指數越高表明樣品的物種多樣性越高。

好奇的同學又有疑問，Shannon指數怎麼算的？

這裡有Shannon指數的公式：

其中，Sobs= 實際測量出的OTU數目； ni= 含有i 條序列的OTU數目；N = 所有的序列數。

4.3 Rank-Abundance曲線

該曲線用於同時解釋樣品多樣性的兩個方面，即樣品所含物種的豐富程度和均勻程度。

橫坐標代表物種排序的數量；縱坐標代表觀測到的相對豐度。樣本曲線的延伸終點的橫坐標位置為該樣本的物種數量

物種的豐富程度由曲線在橫軸上的長度來反映，曲線越寬，表示物種的組成越豐富；物種組成的均勻程度由曲線的形狀來反映，曲線越平坦，表示物種組成的均勻程度越高。如果曲線越平滑下降表明樣本的物種多樣性越高，而曲線快速陡然下降表明樣本中的優勢菌群所佔比例很高，多樣性較低。

但一般超過20個樣本圖就會變得非常複雜而且不美觀！所以假如沒超過20個樣可以考慮該圖哦~

5．Alpha多樣性（樣本內多樣性）

Alpha多樣性是指一個特定區域或者生態系統內的多樣性，常用的度量指標有Chao1 豐富度估計量（Chao1 richness estimator） 、香農 - 威納多樣性指數（Shannon-wiener diversity index）、辛普森多樣性指數（Simpson diversity index）等。

計算菌群豐度：Chao、ace； 計算菌群多樣性：Shannon、Simpson。

Simpson指數值越大，說明群落多樣性越高；Shannon指數越大，說明群落多樣性越高。

看了那麼多指數，可能覺得有點暈，到底每個指數是什麼意思呢？

當然要解釋下咯：

5.1 Chao1：是用chao1 算法計算群落中只檢測到1次和2次的OTU數估計群落中實際存在的物種數。Chao1 在生態學中常用來估計物種總數，由Chao (1984) 最早提出。Chao1值越大代表物種總數越多。

Schao1=Sobs+n1(n1-1)/2(n2+1)

其中Schao1為估計的OTU數，Sobs為觀測到的OTU數，

n1為只有一條序列的OTU數目，n2為只有兩條序列的OTU數目。

5.2 Shannon：用來估算樣品中微生物的多樣性指數之一。它與 Simpson 多樣性指數均為常用的反映 alpha 多樣性的指數。Shannon值越大，說明群落多樣性越高。

5.3 Ace：用來估計群落中含有OTU 數目的指數，由Chao 提出，是生態學中估計物種總數的常用指數之一，與Chao1 的算法不同。

5.4 Simpson：用來估算樣品中微生物的多樣性指數之一，由Edward Hugh Simpson ( 1949) 提出，在生態學中常用來定量的描述一個區域的生物多樣性。Simpson 指數值越大，說明群落多樣性越高。

6.    Alpha多樣性指數差異箱形圖

分別對 Alpha diversity 的各個指數進行秩和檢驗分析（若兩組樣品比較則使用 R 中的wilcox.test 函數，若兩組以上的樣品比較則使用 R 中的 kruskal.test 函數），通過秩和檢驗篩選不同條件下的顯著差異的 Alpha Diversity指數。

7.    Beta多樣性分析（樣品間差異分析）

也許我們有聽說Beta多樣性在最近10年間成為生物多樣性研究的熱點問題之一。

具體解釋下：

Beta多樣性度量時空尺度上物種組成的變化, 是生物多樣性的重要組成部分, 與許多生態學和進化生物學問題密切相關！

7.1 PCoA分析

PCoA（principal co-ordinates analysis）是一種研究數據相似性或差異性的可視化方法，通過一系列的特徵值和特徵向量進行排序後，選擇主要排在前幾位的特徵值，PCoA 可以找到距離矩陣中最主要的坐標，結果是數據矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關係，只是改變了坐標系統。

重要的是，它是可以用來觀察個體或群體間的差異的。

每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構成差異越小。

7.2 PCA分析

主成分分析（Principal component analysis）PCA 是一種研究數據相似性或差異性的可視化方法，通過一系列的特徵值和特徵向量進行排序後，選擇主要的前幾位特徵值，採取降維的思想，PCA 可以找到距離矩陣中最主要的坐標，結果是數據矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關係，只是改變了坐標系統。

詳細關於主成分分析的解釋推薦大家看一篇文章，http://blog.csdn.net/aywhehe/article/details/5736659

一起來看看包含PCoA研究的文章

案例解析

研究背景：全球塑料產量飛速增長，而且呈持續上升的趨勢，因此導致大量塑料廢物排放到環境中，從沿海河口到大洋環流，從東大西洋到南太平洋海域。塑料廢棄物具有化學穩定性和生物利用率低的特點，可長期存在於海洋中，從而影響海洋環境包括海洋生物的生存。

作為一個獨特的底物，塑料碎片可以吸附海洋中的微生物並形成個「塑性球」。以生物膜形式存在於塑料碎片上的微生物群落。許多研究表明，無論是在海洋還是淡水生態系統中，附著在塑料碎片上微生物群落的組成明顯不同於周圍環境（水和沉積物），而且易受位置、時間和塑料類型的影響。

主要圖表

兩兩群落差異指數的PCoA圖

PCoA 圖可以清楚地看到，SW區細菌群落的置信橢圓與pd和sd的置信橢圓有顯著的偏差(p<0.05)，而sd上細菌群落的置信橢圓幾乎覆蓋了pd的置信橢圓(p>0.05)，這表明pd和sd上的細菌群落有相似之處。

不同樣本和處理下的細菌群落（ 前 10 位）豐度分布

底物(SW、SD和Pd)上的主要屬為細菌和假互斥單胞菌，暴露兩周後，這些菌可能是分布廣泛和適應性強的三種底物(SW、SD和PD)。暴露4周後，弧菌相對豐度增加.此外，暴露6周後，自養細菌(如扁平菌和硝酸菌)的數量增加。這三種底物上個細菌群落的生長模式也與3.2的結果一致。圖5還顯示，在6個星期內，在429個原位點中，假單胞菌在pd上的相對豐度高於sw和sd(anova，p<0.05)。

研究結論：首先，營養物質 (TN 和 TP) 與生物膜的平均生長速率呈正相關，而鹽度與生物膜的平均生長速率呈負相關。鹽度是影響PD的個細菌多樣性的主要因素，而溫度、溶解氧和養分(TN和TP)在類似的鹽度條件下可能具有二次效應。儘管種聚合物類型對PD上的細菌群落的多樣性具有較少的影響，但是在細菌群落中的一些屬顯示對PD的聚合物類型的選擇性，並且傾向於將其優選的基質定殖。大的相對豐度SW、PD、SD間屬顯著差異。鹽度是改變河口地區Pd條件致病菌富集的主要因素。另外，在種病原物種豐富的基礎上，PD具有較高的致病性。

7.3 NMDS分析（非度量多維尺度分析）

NMDS（Nonmetric Multidimensional Scaling）常用於比對樣本組之間的差異，可以基於進化關係或數量距離矩陣。

每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構成差異越小。

7.4 排序分析

PCA，PcoA，NMDS分析都屬於排序分析（Ordination analysis）。

排序(ordination)的過程就是在一個可視化的低維空間或平面重新排列這些樣本。

目的：使得樣本之間的距離最大程度地反映出平面散點圖內樣本之間的關係信息。

排序又分兩種：非限制性排序和限制性排序。

1、非限制性排序(unconstrained ordination)

——只使用物種組成數據的排序

(1) 主成分分析(principal components analysis,PCA)

(2) 對應分析(correspondence analysis, CA)

(3) 去趨勢對應分析(Detrended correspondence analysis, DCA)

(4) 主坐標分析(principal coordinate analysis, PCoA)

(5) 非度量多維尺度分析(non-metric multi-dimensional scaling, NMDS)

2、限制性排序(constrained ordination)

——同時使用物種和環境因子組成數據的排序

(1) 冗餘分析(redundancy analysis,RDA)

(2) 典範對應分析(canonical correspondence analysis, CCA)

比較PCA和PCoA

在非限制性排序中，16S和宏基因組數據分析通常用到的是PCA分析和PCoA分析，兩者的區別在於：

PCA分析是基於原始的物種組成矩陣所做的排序分析，而PCoA分析則是基於由物種組成計算得到的距離矩陣得出的。

在PCoA分析中，計算距離矩陣的方法有很多種，包括如：Euclidean, Bray-Curtis, and Jaccard，以及(un)weighted Unifrac (利用各樣品序列間的進化信息來計算樣品間距離，其中weighted考慮物種的豐度，unweighted沒有對物種豐度進行加權處理)。

8.    組間菌群比較選取物種標誌物

8.1（屬水平）樣本-物種豐度關聯circos弦裝圖

樣本與物種的共線性關係circus 圖是一種描述樣本與物種之間對應關係的可視化圈圖，該圖不僅反映了每個樣本的優勢物種組成比例，同時也反映了各優勢物種在不同樣本之間的分布比例。

圖解讀：樣本與物種的共線性關係圖，左半邊表示樣本屬物種豐度情況。右半邊表示屬水平在不同樣本中的分布比例情況。在最內一圈：左邊不同顏色代表不同物種，寬度表示物種豐度，圈外數值表示物種豐度刻度值。一端連接右邊的樣本，不同顏色代表不同樣本，條帶端點寬度表示該樣本中對應物種的比例分布。最外兩圈：左邊不同顏色表示不同樣本在某一物種的比例，右邊不同顏色表示不同物種在某一樣本中的比例。

8.2 Ternary三元相圖

三元相圖是重心圖的一種，它有三個變量，在一個等邊三角形坐標系中，圖中某一點的位置代表三個變量間的比例關係。這裡表示三組樣本之間優勢物種的差異，通過三元圖可以展示出不同物種在分組中的比重關係。

圖解讀：三角分別代表三個或三組樣本，圖中的圓分別代表排名最高哦的屬水平的物種，三種顏色分別代表三組不同分組的優勢物種，圓圈大小代表物種的相對豐度，圓圈理哪個頂點接近，表示此物種在這個分組中的含量較高。該分析僅限三個樣本或三組樣本之間分析比較。

8.3 相關係數圖

通過R 軟體的corrplot 包繪製spearman 相關性熱圖，並通過該熱圖可以發現優勢物種/樣本之間重要的模式與關係。

圖解讀：藍色系的為正相關，紅色系的為負相關，×表示檢驗水平下無意義。越靠近顏色條兩頭，相關係數越大。所以說，我們可以通過實心圓的顏色和大小判斷相關的方向和相關係數的大小。

9.    LDA差異貢獻分析

如果說 PCA，它所作的只是將整組數據整體映射到最方便表示這組數據的坐標軸上，映射時沒有利用任何數據內部的分類信息，是無監督的。

那麼LDA是有監督的，增加了種屬之間的信息關係後，結合顯著性差異標準測試(克魯斯卡爾-沃利斯檢驗和兩兩Wilcoxon測試)和線性判別分析的方法進行特徵選擇。

兩者相同點：

都可以對數據進行降維，降維時都採用了矩陣特徵分解的思想。

差異：

1）LDA是有監督學習的降維方法，而PCA是無監督的降維方法。（註：監督學習是從標記的訓練數據來推斷一個功能的機器學習任務。）

2）LDA選擇分類性能最好的投影方向，而PCA選擇樣本點投影具有最大方差的方向。

除了可以檢測重要特徵，他還可以根據效應值進行功能特性排序，這些功能特性可以解釋大部分生物學差異。這部分希望能詳細了解的同學可以參考這篇文章http://blog.csdn.net/sunmenggmail/article/details/8071502 。

不同顏色代表不同樣本或組之間的顯著差異物種。

使用LefSe軟體分析獲得，其中顯著差異的logarithmic LDA score設為2。

LDA分析究竟能做什麼

組間差異顯著物種又可以稱作生物標記物（biomarkers），這個LDA分析主要是想找到組間在豐度上有顯著差異的物種。

·         案例解析

研究背景：研究表明遺傳和環境影響都在I型糖尿病的發展中起作用，增加的遺傳風險不足以引起疾病，環境因素也是需要的，而且起著至關重要的作用。腸道菌群也許就是這個重要的環境因素，腸道菌群在免疫系統的成熟中起重要作用，此外還影響自身免疫疾病發展。

不同遺傳風險兒童的LDA差異菌群

不同遺傳風險分組中包含的常見菌屬，部分存在特定分組中

PCoA分析揭示不同遺傳風險兒童腸道菌群的在不同地域樣本中均存在顯著差異

點評：針對I型糖尿病疾病發生過程中遺傳HLA分型風險和對應腸道菌群菌的關聯分析，揭示了特定腸道菌群與宿主特定遺傳風險共同作用推進疾病發生。某些特定菌屬可能無法在遺傳高風險兒童腸道內定植，可能對疾病發生存在特定作用。此外對於其他遺傳風險的自身免疫疾病也具有重要提示意義，例如乳糜瀉和類風溼性關節炎。

10. 物種進化樹的樣本群落分布圖

這是另一款和LDA長得有點像的圖，當然功能可完全不一樣。它是將不同樣本的群落構成及分布以物種分類樹的形式在一個環圖中展示。數據經過分析後，將物種分類樹和分類豐度信息通過這款軟體GraPhlAn進行繪製 (http://huttenhower.sph.harvard.edu/GraPhlAn )。

其目的是將物種之間的進化關係以及不同樣本的物種分布豐度和最高分布樣本的信息在一個視覺集中的環圖中一次展示，其提供的信息量較其他圖最為豐富。

·         中間為物種進化分類樹

·         不同顏色的分支代表不同的綱（具體的代表顏色見右上角的圖例），

·         接著的外圈的灰色標示字母的環表示的是本次研究中比例最高的15個科（字母代表的科參見左上角的圖例）。

·         之後的外圈提供的是熱力圖，如果樣本數<=10個則繪製樣本，如果樣本數超過10個則按照分組繪製，每一環為一個樣本，根據其豐度繪製的熱力圖。

·         最外圈為柱狀圖，繪製的是該屬所佔比例最高的樣本的豐度和樣本顏色（樣本顏色見環最下方的樣本名字的顏色）。其中熱力圖和柱狀圖取值均為原比例值x10000後進行log2轉換後的值。

11.  物種相關性分析

根據各個物種在各個樣品中的豐度以及變化情況，計算物種之間的相關性，包括正相關和負相關。相關性分析使用CCREPE算法。

怎麼畫的？

首先對原始16s測序數據的種屬數量進行標準化，然後進行Spearman和Pearson秩相關分析並進行統計檢驗，計算出各個物種之間的相關性，之後在所有物種中根據simscore絕對值的大小，挑選出相關性最高的前100組數據，基於Cytoscap繪製共表達分析網絡圖。

網絡圖一般有兩種表現方式：

物種相關性網絡圖A

○ 圖中每一個點代表一個物種，存在相關性的物種用連線連接。

○ 紅色的連線代表負相關，綠色的先代表正相關。

○ 連線顏色的深淺代表相關性的高低。

物種相關性網絡圖B

○ 圖中每一個點代表一個物種

○ 點的大小表示與其他物種的關聯關係的多少

○ 其中與之有相關性的物種數越多，點的半徑和字體越大

○ 連線的粗細代表兩物種之間相關性的大小

連線越粗，相關性越高。

·         案例解析

研究背景：氣候變化導致美國中部草原的降水模式發生變化，對土壤微生物群落構成及代謝影響很大。研究希望明確土壤微生物群落對土壤水分變化的反應，並確定響應的特定代謝特徵。

主要圖表

同一樣本在不同水分含量孵化處理下土壤菌群的變化

受到水分條件影響的土壤菌群代謝途徑和網絡分布

研究結論：土壤乾燥導致土壤微生物組的組成和功能發生顯著變化。相反，潤溼後幾乎沒有變化。由於乾旱導致的土壤水分減少對土壤碳循環和土壤微生物組進行的其他關鍵生物地球化學循環的影響很大。導致滲透保護劑化合物產生的代謝途徑受到較大影響。

點評：相對簡單的樣本和實驗設計，但是從多個維度探尋支持土壤微生物群落對溼潤和乾燥表型的反應。 與常見的環境採樣檢測不同，針對同一樣本在對照環境下進行環境控制孵化，然後比較菌群變化可以更為有效的控制背景差異。

12.  聚類分析

根據OTU數據進行標準化處理（1wlog10）之後，選取數目最多的前60個物種，基於R heatmap進行作圖

○ 熱圖中的每一個色塊代表一個樣品的一個屬的豐度○ 樣品橫向排列，屬縱向排列

○ 差異是是否對樣品進行聚類，從聚類中可以了解樣品之間的相似性以及屬水平上的群落構 成相似性。

Tips：

如果聚類結果中出現大面積的白或黑是因為大量的菌含量非常低，導致都沒有數值，可以在繪製之前進行標準化操作，對每一類菌單獨自身進行Z標準化。

·         案例解析

研究背景：妊娠期糖尿病（GDM）的患病率在全球範圍內迅速增加，構成一個重要的健康問題和產科實踐的重大挑戰（Ferrara，2007）。高脂血症是妊娠常見的合併症。在GDM患者中，血脂的生理變化可能導致懷孕期間潛在的代謝紊亂。腸道失調在宿主代謝異常中起著至關重要的作用，最近關於2型糖尿病(T2D)和肥胖的研究就證明了這一點。這些研究表明，妊娠期間腸道微生物ME的主要變化可能在GDM的發展中起著至關重要的作用。

GDM加高脂血症（M隊列）妊娠期間與顯著改變的脂質相關的腸道微生物群（屬）

研究結論：我們的結果表明，血脂水平可能反映了GDM發展過程中的一些異常變化。所鑑定的多種生物標誌物對GDM合併高脂血症的防治有一定的參考價值。

組間物種差異性箱形圖

組間物種差異性盒形圖描述在不同分組之間具有差異顯著的某一物種做盒形圖

圖中以屬水平為例做物種差異性盒形圖，展示如下：

○ 圖中不同顏色代表不同的分組，更直觀顯示組間物種差異

○ 每一個盒形圖代表一個物種，圖上方是物種名。

Anosim檢驗

Anosim分析是一種非參數檢驗，用來檢驗組間的差異是否顯著大於組內差異，從而判斷分組是否有意義

展示如下：

R-value介於（-1，1）之間，R-value大於0，說明組間差異顯著。

R-value小於0,說明組內差異大於組間差異。

統計分析的可信度用 P-value 表示，P< 0.05 表示統計具有顯著性。

對Anosim的分析結果，基於兩兩樣本之間的距離值排序獲得的秩（組間的為between，組內的為within），這樣任一兩兩組的比較可以獲得三個分類的數據，並進行箱線圖的展示（若兩個箱的凹槽互不重疊，則表明它們的中位數有顯著差異）

13.  隨機森林分類樹屬分類效果

隨機森林是機器學習算法的一種，它可以被看作是一個包含多個決策樹的分類器。其輸出的分類結果是由每棵決策樹「投票」的結果。由於每棵樹在構建過程中都採用了隨機變量和隨機抽樣的方法，因此隨機森林的分類結果具有較高的準確度，並且不需要「減枝」來減少過擬合現象。

隨機森林可以有效的對分組樣品進行分類和預測。

物種重要性點圖。橫坐標為重要性水平，縱坐標為按照重要性排序後的物種名稱。上圖反映了分類器中對分類效果起主要作用的菌屬，按作用從大到小排列。

Error rate: 表示使用下方的特徵進行隨機森林方法預測分類的錯誤率，越高表示基於菌屬特徵分類準確度不高，可能分組之間菌屬特徵不明顯。圖中以所有水平為例，取前60個作圖。

14. ROC曲線圖

ROC 曲線指受試者工作特徵曲線(receiver operating characteristic curve), 是反映敏感性和特異性連續變量的綜合指標，通過構圖法揭示敏感性和特異性的相互關係。

ROC 曲線將連續變量設定出多個不同的臨界值，從而計算出一系列敏感性和特異性，再以敏感性為縱坐標、（1-特異性）為橫坐標繪製成曲線。

曲線下面積越大，診斷準確性越高。展示如下：

15. FAPROTAX生態功能預測

FAPROTAX是一款在2016年發表在SCIENCE上的較新的基於16S測序的功能預測軟體。它整合了多個已發表的可培養菌文章的手動整理的原核功能資料庫，資料庫包含超過4600個物種的7600多個功能注釋信息，這些信息共分為80多個功能分組，其中包括如硝酸鹽呼吸、產甲烷、發酵、植物病原等。

FAPROTAX對環境樣本更友好

如果說PICRUSt（後續會介紹）在腸道微生物研究更為適合，那麼FAPROTAX尤其適用於生態環境研究，特別是地球化學物質循環分析。

FAPROTAX適用於對環境樣本（如海洋、湖泊等）的生物地球化學循環過程（特別是碳、氫、氮、磷、硫等元素循環）進行功能注釋預測。因其基於已發表驗證的可培養菌文獻，其預測準確度可能較好，但相比於上述PICRUSt和Tax4Fun來說預測的覆蓋度可能會降低。

FAPROTAX可根據16S序列的分類注釋結果對微生物群落功能（特別是生物地化循環相關）進行注釋預測。

圖中橫坐標代表樣本，縱坐標表示包括碳、氫、氮、硫等元素循環相關及其他諸多功能分組。可快速用於評估樣品來源或特徵。

17.基於BugBase的表型分類比較

Bugbase也是16年所提供服務的一款免費在線16S功能預測工具，到今年才發表文章公布其軟體原理。該工具主要進行表型預測，其中表型類型包括革蘭氏陽性、革蘭氏陰性、生物膜形成、致病性、移動元件、氧需求，包括厭氧菌、好氧菌、兼性菌）及氧化脅迫耐受等7類。

Gram Negative 革蘭氏陰性菌

18. Picrust群落功能差異分析

通過對已有測序微生物基因組的基因功能的構成進行分析後，我們可以通過16s測序獲得的物種構成推測樣本中的功能基因的構成，從而分析不同樣本和分組之間在功能上的差異（PICRUSt Nature Biotechnology, 1-10. 8 2013）。

Picrust對腸道菌群樣本更友好

通過對宏基因組測序數據功能分析和對應16s預測功能分析結果的比較發現，此方法的準確性在84%-95%，對腸道微生物菌群和土壤菌群的功能分析接近95%，能非常好的反映樣品中的功能基因構成。

怎麼做出來的？

為了能夠通過16s測序數據來準確的預測出功能構成，首先需要對原始16s測序數據的種屬數量進行標準化，因為不同的種屬菌包含的16s拷貝數不相同。

然後將16s的種屬構成信息通過構建好的已測序基因組的種屬功能基因構成表映射獲得預測的功能結果。（根據屬這個水平，對不同樣本間的物種豐度進行顯著性差異兩兩檢驗，我們這裡的檢驗方法使用STAMP中的two-sample中T-TEST方法，Pvalue值過濾為0.05，作Extent error bar圖。）

此處提供COG，KO基因預測以及KEGG代謝途徑預測。當然，躍躍欲試的小夥伴也可自行使用我們提供的文件和軟體（STAMP）對不同層級以及不同分組之間進行統計分析和製圖，以及選擇不同的統計方法和顯著性水平。

這裡提到的STAMP有些小夥伴說不太了解，別急，後面會有更多介紹。

18.1 COG構成差異分析圖

圖中不同顏色代表不同的分組，列出了COG構成在組間存在顯著差異的功能分類以及在各組的比例，此外右側還給出了差異的比例和置信區間以及P-value。

18.2 KEGG代謝途徑差異分析圖

通過KEGG代謝途徑的預測差異分析，我們可以了解到不同分組的樣品之間在微生物群落的功能基因在代謝途徑上的差異，以及變化的高低。為我們了解群落樣本的環境適應變化的代謝過程提供一種簡便快捷的方法。

本例圖所顯示的是第三層級的KEGG代謝途徑的差異分析，也可以針對第二或第一層的分級進行分析。

圖中不同顏色代表不同的分組，列出了在第三層級的構成在組間存在顯著差異的KEGG代謝途徑第三層分類以及在各組的比例，此外右側還給出了差異的比例和置信區間以及P-value。

·         案例解析

研究背景：儘管普遍認為腸道微生物組的生態多樣性和分類組成在肥胖和T2D中發生改變，但與單個微生物或微生物產物的關聯在研究之間不一致。缺乏大樣本群體研究，從而確定腸道微生物組，血漿代謝組，肥胖和糖尿病表型以及環境因素之間的幾種關聯。

主要圖表：

按照肥胖和糖尿病對人群分為三組，同時進行了16S，代謝和宏基因組的檢測。

與肥胖相關的菌屬以及代謝途徑

研究結論：確定了腸道微生物組，血漿代謝組，肥胖和糖尿病表型以及環境因素之間的幾種關聯。與腸道微生物組變異相關的主要是肥胖，不是2型糖尿病。存在與腸道微生物組變異相關的藥物和膳食補充劑。高鐵攝入量影響小鼠的腸道微生物組成。微生物組變異也反映在血清代謝物譜中。

點評：

相對大人群的隊列研究，同時涵蓋了菌群、代謝和疾病表型以及膳食補充調查的數據。 從結果看菌屬和血漿代謝存在關聯，但是貢獻度都較低，如果樣本數量不足很可能找不到顯著的聯繫，這也是這類大樣本隊列研究的意義。 本研究在人群分組時針對性的研究了肥胖-II型糖尿病和菌群的關聯，因而構建了三個主要分組人群，結果顯示肥胖與菌群的關聯度更大，解釋了大部分的菌群差異，而糖尿病的菌群變化較小。 本研究其中較為重要的是發現了不同膳食補充對菌群的影響，並在小鼠實驗中得到證實。

20. 基因的差異分析圖

除了能對大的基因功能分類和代謝途徑進行預測外，我們還能提供精細的功能基因的數量和構成的預測，以及進行樣本間以及組間的差異分析，並給出具有統計意義和置信區間的分析結果。

這一分析將我們對於樣本群落的差異進一步深入到了每一類基因的層面。

圖中不同顏色代表不同的分組，列出了在組間/樣本間存在顯著差異的每一個功能基因（酶）以及在各組的比例，此外右側還給出了差異的比例和置信區間以及P-value。

21. 貢獻圖

我們通過計算每個變量正常計數中值，進一步確定每個被選擇的OTU的特徵。如果某一變量的中位數數高於任何其他變量，則OTU被定義為對變量有貢獻。其中每個OTU條長度對應於多元模型中特徵的重要性(對於每個組件上的特定特徵，具有正號或負號的多元回歸係數)通過從底部開始降低重要性進行排序，並且顏色與貢獻變量相匹配。貢獻圖可以顯示任意指定級別的細菌分類。

圖解讀：加載在comp1組件和comp2組件上貢獻最大的OTU圖。顏色代表不同分組。條形圖越長說明對應OTU在此分組中貢獻最大。

22. 環境因子分析

冗餘分析（redundancy analysis, RDA）或者

典範對應分析（canonical correspondence analysis, CCA）都是基於對應分析發展的一種排序方法，將對應分析與多元回歸分析相結合，每一步計算均與環境因子進行回歸，又稱多元直接梯度分析。主要用來反映菌群與環境因子之間的關係。

RDA 是基於線性模型，CCA是基於單峰模型。分析可以檢測環境因子、樣品、菌群三者之間的關係或者兩兩之間的關係。

冗餘分析可以基於所有樣品的OTU作圖，也可以基於樣品中優勢物種作圖；

箭頭射線：箭頭分別代表不同的環境因子；

夾角：環境因子之間的夾角為銳角時表示兩個環境因子之間呈正相關關係，鈍角時呈負相關關係。環境因子的射線越長，說明該影響因子的影響程度越大；

不同顏色的點表示不同組別的樣品或者同一組別不同時期的樣品，圖中的拉丁文代表物種名稱，可以將關注的優勢物種也納入圖中；

○ 環境因子數量要少於樣本數量，同時在分析時，需要提供環境因子的數據，比如 pH值，測定的溫度值等。

23. 升級版的Alpha多樣性指數

基於Alpha多樣性指數中的Invsimpson(逆辛普森指數)繪製。它是均勻度一致的群落豐富程度的標誌。與其他Alpha多樣性指數相比，逆辛普森指數不容易受到抽樣誤差的影響。一個樣本有2個及以上觀測值即可生成該圖。但樣本過多就會讓人眼花繚亂，可以只挑選感興趣的樣本作圖。

圖解讀：不同顏色為不同分組，縱坐標為指數值。指數越高說明其多樣性與豐富程度越高。

個人更推薦將上圖的Invsimpson(逆辛普森指數)用於組間作圖（見下圖）。

樣本量比較大時可以用該圖，例如一個樣本有5個或以上觀測值時，只要是連續的值都可以。當組間群落多樣性都較為均勻時，可以通過添加Invsimpson(逆辛普森指數)來展示細微的組間差異。

圖解讀：用log10轉換的y軸顯示的腸道微生物組的細菌α-多樣性（通過辛普森指數反向測定）。通過圖中的連線可以明顯看出菌群多樣性在隨時間變化的趨勢。

24. 樣本變化軌跡圖

對於使用數據的要求：

1.一個樣本有多個觀測值，示例圖中一個樣本有三個觀測值。如果觀測值為多個連續的值，例如採樣時間或用藥時間點，這樣時間軌跡的變化會更明顯。

2.pcoa排序結果表明組間有顯著差異。而這種差異特徵與時間有相關性。因為點與點的距離是基於pcoa矩陣得來的，如果組間差異不明顯，不同組之間的點與點的距離較近，時間軌跡的連線就會很亂，這樣也表達不出相關的生物學意義。例如右圖是該項目的bray距離的pcoa圖，可以發現pcoa圖中組間差異不明顯。

圖解讀：不同顏色代表不同組，形狀不同的點表明了不同的時間點，點之間的連線表明時間軌跡，。兩點之間距離越近表明兩者的群落構成差異較小。橫軸表示儘可能最大解釋數據變化的主坐標成分，縱軸表示解釋餘下的變化度中佔比例最大的主坐標成分。

25. spls(稀疏偏最小二乘）回歸分析

sPLS回歸允許整合微生物群落數據矩陣和臨床變量矩陣以進行多元回歸。它可以處理數據中的共線性和噪聲，並且適合對多個響應變量進行建模。

這需要有大量的meta信息，例如一個樣本有幾十個臨床信息，你想知道這些信息與腸道菌群的相關性是怎樣的，我們將這些臨床信息利用adonis2檢驗它們與腸道菌群間是否有統計學意義。然後將具有統計學意義的信息利用spls按照它們之間的相關性從大到小排列。數據間的相關性越強越能很好的使用此分析。

 a

圖解讀：

a圖. 前兩個sPLS維度的相關圓圖顯示了> 0.2/< - 0.2的相關性。兩個灰色圓圈表示相關係數為0.5和1.0。OUT顯示為較小的圓點，根據所屬的cluster進行著色。表示變量的圓點附帶了標籤。距離較近的變量之間呈正相關，投影方向相反的變量之間呈負相關。彼此垂直放置的變量不相關。OTU解釋的方差在Component 1上為2.94%,在Component 2為8.77%.

b圖. 前兩個sPLS維度的聚類圖像映射，顯示了OTUs(右側)和臨床變量(底部)之間的兩兩相關。紅色和藍色分別表示正相關和負相關。在基於sPLS回歸模型的mixOmics cim()函數內進行層次聚類(聚類方法: complete linkage，距離法:Pearson相關)。

c圖. 分別在Component 1和Component 2上貢獻最大的OTU的荷載圖。長方形條狀是根據它們所屬的簇而著色的。各OTU的分類信息根據顏色著色(圖例見b圖)

可變區和測序選擇

目前針對擴增子測序可選擇的測序平臺和方案很多，不同平臺的讀長和適用的測序區段以及優勢各有不同。16s測序主要的測序區段包括V4、V3V4，V1V2，V6，此外還有全長等不同的區段選擇，不同可變區或全長由於引物的不同以及不同種屬相應區段內的變異多樣性差異，對菌屬的豐度評估會有一定的差異。

從長度來看，全長16S長度為1.5kb左右，單菌落的16S全長sanger一代測序仍然是菌種鑑定的主要手段，納米孔和Pacbio的三代測序可以高通量的獲得全長序列，對於希望更高解析度的分析菌種的研究有一定優勢。三代的測序準確度目前逐漸改進，直接測序準確度可以在90%以上，糾錯後可以提高到97~99%以上，已足夠提供高精度的分類。三代目前主要問題在於建庫成本相對較高，通過使用barcode可以降低部分但仍然偏高，此外普遍測序深度相對於二代測序要低許多。目前最主要的可變區選擇是V4區和V3V4區，V4區長度為256bp左右，加上兩側引物長度為290bp左右，使用雙端2x250bp或2x150bp可以測通，此外如454、life、Illumina Hiseq 4000的測序平臺讀長也可以主要涵蓋該區段讀長。例如採用Illumina Hiseq測序平臺對該項目進行雙端測序(Paired-end)，測序得到了fastq格式的原始數據（樣本對應一對序列S_1.fastq和S_2.fastq）。再配對拼接成單條序列。其引物通用性相對是所有可變區中最高的，大量的大規模菌群調查研究都採用V4區作為檢測區域，包括人體菌群研究如：HMP，腸道菌群如美國腸道計劃AGP，歐洲的FGFP等，以及全球土壤菌群調查，目前仍然是國際研究中使用最廣泛和認可的檢測區域。

Illumina的Miseq提供了長達2x300bp以及Hiseq2500和最近的NovoSeq提供有2x250bp的測序方案，為進一步利用讀長，目前有相當一部分研究選擇V3V4區，該區段長度在460bp左右，相較於V4度多出了V3區段約100bp左右的片段，在少部分菌屬中可以增加一定解析度。經過對比，V3V4區的檢測結果和V4區在絕大部分菌屬中的豐度一致，但由於引物不同，在少量菌屬中豐度會有不同偏向，V3V4從OTU層面上並未發現較V4區有明顯增加。引物的選擇和提取、儲存方法是影響菌群檢測豐度構成的主要因素，不同研究之間的比較需要考慮到實驗方案的一致，相同的方案可以直接比較。目前的高通量測序平臺可以較低成本的進行大規模的測序，從測序深度角度，土壤菌群的多樣性最高，一般需要5萬條以上序列可以達到飽和，腸道樣本在3萬條以上，水體和尿液等1萬條以上基本可以到達飽和。

同一批小鼠糞便樣本v4（10萬 clean reads）和 v3v4(5萬clean reads)測序數據比較:

原始序列數據：

V4

V3V4

以上兩表是對原始序列數據進行統計，表中可以看出有效序列tags、高質量序列clean_tags、otus數量  V4區都遠高於v3v4區。V4區測序獲得下機數據在13萬條左右，v4區測序獲得的下機數據在5萬條左右。

  Alpha多樣性指數比較：

V4

V3V4

以上兩個表分別是對Alpha多樣性指數計算的結果比較

Chao1 指數和ACE指數是用來評估樣本中所含OTU數目的指數，從Chao1 指數和ACE指數可以看出，用 v4測序獲得的結果要明顯大於v3v4的結果。這是因為v4測序通量更高，測序深度更好，每個樣下機的測序數據可以到10萬條以上，一般在13萬條左右，所以經過序列比對獲得的OTU數目更多，相比較用v3v4測序每個樣下機的數據大約在4到5萬條左右，經過序列比對獲得的OTU相對少一點。

Shannon指數和Simpson指數是用來評估菌群的豐富度和均一度 的。從Shannon指數和Simpson指數，用v4和v3v4測序指數相差不大，或v4比v3v4略高一點，證明兩種測序之間菌群的豐富度多樣性和均一度叫接近。

物種主要構成比較：

V4

V3V4

V3v4

屬水平前10個物種構成：Lactobacillus、Adlercreutzia、Flexispira、Allobaculum、Desulfovibrio、Prevotella、Odoribater、Oscillospira、[Prevotella]、Bacteroides

V4

屬水平前10個物種構成：Lactobacillus、Akkermansia、Helicobacter、Allobaculum、Desulfovibrio、Adlercreutzia、Odoribacter、Bacteroides、Prevotella、[Prevotella]

從前10個物種構成來看，有8個是相同的，物種的主要構成基本一致，測序的穩定性較好。從種類來看，v3v4測到的屬水平個數較多。

各分類水平鑑定到的物種種類比較：

V4

V3v4

以上兩張表代表了每個樣本在各分類水平上鑑定到的物種種類數。從整體上來看，分別用v4和v3v4測序得到的數據，在各分類水平上鑑定到的物種個數相對比較穩定和接近，（尤其在目水平和科水平上）用v3v4測序獲得的物種數比v4相對較多一點，單相差不大，在屬水平和種水平則不一定是這種規律，最終鑑定到的物種個數也跟該樣本的測序質量有關。

最後附幾篇頂級雜誌發表的16s v4區的文章

 Poyet, M., et al. "A library of human gut bacterial isolates paired with longitudinal multiomics data enables mechanistic microbiome research." Nature medicine 25.9 (2019): 1442-1452.

16S library preparation and sequencing. 16S rRNA gene libraries targeting the V4 region of the 16S rRNA gene were prepared by first normalizing template concentrations and determining optimal cycle number by way of qPCR. Two 25 µL reactions for each sample were amplified with 0.5 units of Phusion with 1X High Fidelity buffer, 200 μM of each dNTP, 0.3 μM of 515 F( 5′- AATGATACGGCGACCACCGAGATCTACACTATGGTAATTGTGTGCCAGCMGCCGCGGTAA-3′) and 806rcbc0 (5′- CAAGCAGAAGACGGCATACGAGATTCCCTTGTCTCCAGTCAGTCAGCCGGACTACHVGGGTWTCTAAT-3′).

Tito, Raul Y., et al. "Population-level analysis of Blastocystis subtype prevalence and variation in the human gut microbiota." Gut 68.7 (2019): 1180-1189.

We profiled stool samples from 616 healthy individuals from the FGFP cohort as well as 107 patients with IBD using amplicon sequencing targeting the V4 variable region of the 16S rRNA and 18S rRNA genes.

Call, Lee, et al. "Metabolomic signatures distinguish the impact of formula carbohydrates on disease outcome in a preterm piglet model of NEC." Microbiome 6.1 (2018): 111.

Gut contents and mucosal samples were collected and analyzed for microbial profiles by sequencing the V4 region of the 16S rRNA gene. Metabolomic profiles of cecal contents and plasma were analyzed by LC/GC mass spectrometry

Wang, Chao, et al. "High-salt diet has a certain impact on protein digestion and gut microbiota: a sequencing and proteome combined study." Frontiers in Microbiology 8 (2017): 1838.

In this study, C57BL/6J mice were fed low- or high-salt diets (0.25 vs. 3.15% NaCl) for 8 weeks, and then gut contents and feces were collected. Fecal microbiota was identified by sequencing the V4 region of 16S ribosomal RNA gene.

Bai, J., Y. Hu, and D. W. Bruner. "Composition of gut microbiota and its association with body mass index and lifestyle factors in a cohort of 7–18 years old children from the American Gut Project." Pediatric obesity 14.4 (2019): e12480.

AGP sequenced the V4 region of 16S rRNA gene

Luthold, Renata V., et al. "Gut microbiota interactions with the immunomodulatory role of vitamin D in normal individuals." Metabolism 69 (2017): 76-86.

The association between 25(OH)D and fecal microbiota (16S rRNA sequencing, V4 region) was tested by multiple linear regression.

Iszatt, Nina, et al. "Environmental toxicants in breast milk of Norwegian mothers and gut bacteria composition and metabolites in their infants at 1 month." Microbiome 7.1 (2019): 34.

Child fecal samples were characterized by 16S rRNA gene amplicon sequencing of the V4 region. We used Deblur, a novel sub-operational taxonomic-unit (sub-OTU) approach that provides a higher resolution than OTU-based analyses.

Vangay, Pajau, et al. "US immigration westernizes the human gut microbiome." Cell 175.4 (2018): 962-972.

We performed amplicon-based sequencing of the 16S rRNA gene V4 region on 550 stool samples (one sample per participant).

Suez, Jotham, et al. "Post-antibiotic gut mucosal microbiome reconstitution is impaired by probiotics and improved by autologous FMT." Cell 174.6 (2018): 1406-1423.

For 16S amplicon pyrosequencing, PCR amplification was performed spanning the V4 region using the primers 515F/806R of the 16S rRNA gene and subsequently sequenced using 2X250 bp paired-end sequencing (Illumina MiSeq).

Zmora, Niv, et al. "Personalized gut mucosal colonization resistance to empiric probiotics is associated with unique host and microbiome features." Cell 174.6 (2018): 1388-1405.

For 16S amplicon pyrosequencing, PCR amplification was performed spanning the V4 region using the primers 515F/806R of the 16S rRNA gene and subsequently sequenced using 2 × 250 bp paired-end sequencing (Illumina MiSeq).

Riquelme, Erick, et al. "Tumor microbiome diversity and composition influence pancreatic cancer outcomes." Cell 178.4 (2019): 795-806.

The 16S rDNA V4 region was amplified by PCR and sequenced in the MiSeq platform (Illumina) using the 2x250 bp paired-end protocol yielding pair-end reads that overlap almost completely. The primers used for amplification contain adapters for MiSeq sequencing and single-index barcodes so that the PCR products may be pooled and sequenced directly (Caporaso et al., 2012), targeting at least 10,000 reads per sample. 16S (variable region 4 [v4]) rRNA gene pipeline data incorporated phylogenetic and alignment based approaches to maximize data resolution.

Matson, Vyara, et al. "The commensal microbiome is associated with anti–PD-1 efficacy in metastatic melanoma patients." Science 359.6371 (2018): 104-108.

Specifically, the V4 region of the 16S rRNA gene (515F-806R) was PCR-amplified with region-specific primers that include sequencer adapter sequences used in the Illumina flowcell.

Raman, Arjun S., et al. "A sparse covarying unit that describes healthy and impaired human gut microbiota development." Science 365.6449 (2019): eaau4735.

Amplicons generated from variable region 4 (V4) of bacterial 16S rRNA genes present in these 2455 fecal samples were sequenced, and the resulting reads were assigned to operational taxonomic units with ≥97% nucleotide sequence identity (97%ID OTUs).

Gehrig, Jeanette L., et al. "Effects of microbiota-directed foods in gnotobiotic animals and undernourished children." Science365.6449 (2019): eaau4732.

Characterizing human fecal microbial communities Methods for V4-16S rRNA gene sequencing and data analysis, calculation of MAZ scores and functional microbiome maturity, and quantification of enteropathogen burden by means of multiplex quantitative polymerase chain reaction (qPCR) are described in the supplementary materials.

Lloyd-Price, Jason, et al. "Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases." Nature 569.7758 (2019): 655.

In brief, bacterial genomic DNA was extracted from the total mass of the biopsied specimens using the MoBIO PowerLyzer Tissue and Cells DNA isolation kit and sterile spatulas for tissue transfer. The 16S rDNA V4 region was amplified from the extracted DNA by PCR and sequenced in the MiSeq platform (Illumina) using the 2 × 250 bp paired-end protocol, yielding pair-end reads that overlapped almost completely.

Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases. Nature. 2019

In brief, bacterial genomic DNA was extracted from the total mass of the biopsied specimens using the MoBIO PowerLyzer Tissue and Cells DNA isolation kit and sterile spatulas for tissue transfer. The 16S rDNA V4 region was amplified from the extracted DNA by PCR and sequenced in the MiSeq platform (Illumina) using the 2 × 250 bp paired-end protocol, yielding pair-end reads that overlapped almost completely.

emporal development of the gut microbiome in early childhood from the TEDDY study. Nature. 2019

Bacterial DNA was extracted using the PowerMag Microbiome DNA isolation kit following the manufacturer’s instructions. The V4 region of the 16S rRNA gene was amplified by PCR and sequenced on the MiSeq platform (Illumina) using the 2 × 250 bp paired-end read protocol.

A communal catalogue reveals Earth’s multiscale microbial diversity. Nature. 2018

We surveyed bacterial and archaeal diversity using amplicon sequencing of the 16S rRNA gene, a common taxonomic marker for bacteria and archaea12 that remains a valuable tool for microbial ecology despite the introduction of whole-genome methods (e.g., metagenomics) that capture gene-level functional diversity13. We amplified the 16S rRNA gene (V4 region) using primers14 shown to recover sequences from most bacterial taxa and many archaea.

Root microbiota drive direct integration of phosphate stress and immunity. Nature. 2017.

For wild soil experiment 16S sequencing, we processed libraries according to Caporaso, et al.28. Three sets of index primers were used to amplify the V4 (515F-806R) region of the 16S rRNA gene of each sample. In each case, the reverse primer had a unique molecular barcode for each sample.