做過16s測序的小夥伴們都知道
測完之後會拿到一份結果報告
但這並不代表可以開始寫文章了
看似一大堆數據圖表卻不知如何下手
這是很多人頭疼的地方
那麼怎樣給報告中的數據賦予靈魂
讓它真正成為對你有幫助的分析呢？
今天我們來詳細解讀下。
沒有接觸過微生物16s測序的同學也不要慌，接下來的乾貨供（bi）你（mian）參（cai）考（keng），輕鬆搞定一篇文章。
【此處是新手上路的充電時刻，老司機請自行跳過前幾個問題往下看。】
16S rRNA 基因是編碼原核生物核糖體小亞基的基因，長度約為1542bp，其分子大小適中，突變率小，是細菌系統分類學研究中最常用和最有用的標誌。
16S rRNA基因序列包括9個可變區和10個保守區，保守區序列反映了物種間的親緣關係， 而可變區序列則能體現物種間的差異。
16S rRNA基因測序以細菌16S rRNA基因測序為主,核心是研究樣品中的物種分類、物種豐度以及系統進化。
目前二代測序是一個邊合成邊測序的過程，使用的是螢光可逆終止子。每個可逆終止子的鹼基3』端都有一個阻斷基團，而在側邊帶有一種螢光。由於有4種不同的鹼基（ATCG），因此也會有對應4種不同顏色的螢光。開始擴增每次結合上一個鹼基，DNA的擴增便會停止，此時能收到一種螢光信號。然後放試劑除去阻斷基團，進行下一個鹼基的結合，以此類推得到一連串的螢光信號組合序列。而根據螢光的顏色我們便可以確定每一個位點的基因型，即可以得到這一段DNA片段的序列。
在進行實驗設計前，這是有些小夥伴面臨的一個問題。環境樣本由於來源和條件不完全可控，每個樣品之間會存在很大的差異，即便是相同樣本的不同取樣時間和部位也會存在一定的差異。
基於高通量測序主要是為了了解樣品的菌群構成和功能分析，以及尋找不同環境之間的差異，包括菌和功能基因以及代謝。如果僅做單一樣本，很可能結論只能代表這個單一取樣樣本的信息，無法排除不同樣本重複之間的差異，也就可能得不到真正代表環境差異的結果。
所以環境樣品不僅要重複而且還應該以分組方式取儘量多的樣本以全面的代表一個環境條件下的各種變異情況。
確定做重複後，又面臨該怎麼選擇測序區段的問題。目前市面上有v1-v3區/v3-v4區/v4區等可供選擇。
16S rRNA編碼基因序列共有9個保守區和9個高可變區。其中，V4區其特異性好，資料庫信息全，我們通過大量的測序試驗證明用v4區擴增出菌群結果的可以很好的反應樣本的菌群結構用於後續的數據建模分析，是細菌多樣性分析注釋的最佳選擇。
基本確定好後，就要著手開始實驗，實驗完送樣又是個問題，以往給測序公司送樣往往是低溫運輸，且不說麻煩，還要提心弔膽怕運輸過程會不會有什麼問題。為此我們免費提供常溫保存取樣盒，就不用有這樣的顧慮，取樣及運輸全程都只需要常溫即可。
很多小夥伴有過這樣的經歷，在拿到公司出具的報告之後，仍然一頭霧水，幾十頁的報告內容看著豐富卻不知該怎麼運用。我們一起來理一下關鍵圖表的含義。
OTU是我們要搞清的一個重要概念，可以說是後續分析的基石。
OTU(operational taxonomic units) 是在系統發生學研究或群體遺傳學研究中，為了便於進行分析，人為給某一個分類單元（品系，種，屬，分組等）設置的同一標誌。通常按照 97% 的相似性閾值將序列劃分為不同的 OTU，每一個 OTU 通常被視為一個微生物物種。相似性小於97%就可以認為屬於不同的種，相似性小於93%-95%，可以認為屬於不同的屬。樣品中的微生物多樣性和不同微生物的豐度都是基於對OTU的分析。
有了OTU這個概念之後，就不難理解下表。對每個樣本的測序數量和OTU數目進行統計，並且在表栺中列出了測序覆蓋的完整度。
其中 SampleName表示樣本名稱；SampleSize表示樣本序列總數；OTUsNumber表示注釋上的OTU數目；OTUsSeq表示注釋上OTU的樣本序列總數。
Coverage是指各樣品文庫的覆蓋率，其數值越高，則樣本中序列沒有被測出的概率越低。該指數實際反映了本次測序結果是否代表樣本的真實情況。計算公式為：C=1-n1/N 其中n1 = 只含有一條序列的OTU的數目；N = 抽樣中出現的總的序列數目。
下表是對每個樣本在分類字水平上的數量進行統計，並且在表栺中列出了在每個分類字水平上的物種數目
其中SampleName表示樣本名稱；Phylum表示分類到門的OTU數量；Class表示分類到綱的OTU數量；Order表示分類到目的OTU數量；Family表示分類到科的OTU數量；Genus表示分類到屬的OTU數量；Species表示分類到種的OTU數量。
我們可以看到絕大部分的OTU都分類到了屬（Genus），也有很多分類到了種（Species）。但是仍然有很多無法完全分類到種一級，這是由於環境微生物本身存在非常豐富的多樣性，還有大量的菌仍然沒有被測序和發現。
橫坐標中每一個條形圖代表一個樣本，縱坐標代表該分類層級的序列數目或比例。同一種顏色代表相同的分類級別。圖中的每根柱子中的顏色表示該樣本在不同級別（門、綱、目等）的序列數目，序列數目只計算級別最低的分類，例如在屬中計算過了，則在科中則不重複計算。
我們還需要對樣本之間或分組之間的OTU進行比較獲得韋恩圖：
微生物多樣性分析中如何驗證測序數據量是否足以反映樣品中的物種多樣性？
稀釋曲線（豐富度曲線）可以派上用場。它是用來評價測序量是否足以覆蓋所有類群，並間接反映樣品中物種的豐富程度。
不免有同學有疑惑，稀釋曲線怎麼來的？
它是利用已測得16S rDNA序列中已知的各種OTU的相對比例，來計算抽取n個（n小於測得reads序列總數）reads時出現OTU數量的期望值，然後根據一組n值（一般為一組小於總序列數的等差數列）與其相對應的OTU數量的期望值做出曲線來。
至此，我們雖然知道了稀釋曲線的由來，那麼這個五彩繽紛的稀釋曲線該怎麼看呢？
當曲線趨於平緩或者達到平臺期時也就可以認為測序深度已經基本覆蓋到樣品中所有的物種，增加測序數據無法再找到更多的OTU；
反之，則表示樣品中物種多樣性較高，還存在較多未被測序檢測到的物種。
橫坐標代表隨機抽取的序列數量；縱坐標代表觀測到的OTU數量。樣本曲線的延伸終點的橫坐標位置為該樣本的測序數量。
Shannon-Wiener 曲線，是利用shannon指數來進行繪製的，反映樣品中微生物多樣性的指數，利用各樣品的測序量在不同測序深度時的微生物多樣性指數構建曲線，以此反映各樣本在不同測序數量時的微生物多樣性。
當曲線趨向平坦時，說明測序數據量足夠大，可以反映樣品中絕大多數的微生物物種信息。
橫坐標代表隨機抽取的序列數量；縱坐標代表的是反映物種多樣性的Shannon指數，樣本曲線的延伸終點的橫坐標位置為該樣本的測序數量。
其中曲線的最高點也就是該樣本的Shannon指數，指數越高表明樣品的物種多樣性越高。
其中，Sobs= 實際測量出的OTU數目；
ni= 含有i 條序列的OTU數目；N = 所有的序列數。
該曲線用於同時解釋樣品多樣性的兩個方面，即樣品所含物種的豐富程度和均勻程度。
橫坐標代表物種排序的數量；縱坐標代表觀測到的相對豐度。
樣本曲線的延伸終點的橫坐標位置為該樣本的物種數量
物種的豐富程度由曲線在橫軸上的長度來反映，曲線越寬，表示物種的組成越豐富；
物種組成的均勻程度由曲線的形狀來反映，曲線越平坦，表示物種組成的均勻程度越高。
如果曲線越平滑下降表明樣本的物種多樣性越高，而曲線快速陡然下降表明樣本中的優勢菌群所佔比例很高，多樣性較低。
但一般超過20個樣本圖就會變得非常複雜而且不美觀！所以假如沒超過20個樣可以考慮該圖哦~
Alpha多樣性是指一個特定區域或者生態系統內的多樣性，常用的度量指標有Chao1 豐富度估計量（Chao1 richness estimator） 、香農 - 威納多樣性指數（Shannon-wiener diversity index）、辛普森多樣性指數（Simpson diversity index）等。
計算菌群豐度：Chao、ace；
計算菌群多樣性：Shannon、Simpson。
Simpson指數值越大，說明群落多樣性越高；Shannon指數越大，說明群落多樣性越高。
看了那麼多指數，可能覺得有點暈，到底每個指數是什麼意思呢？
當然要解釋下咯：
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數目。
Shannon：用來估算樣品中微生物的多樣性指數之一。它與 Simpson 多樣性指數均為常用的反映 alpha 多樣性的指數。Shannon值越大，說明群落多樣性越高。
Ace：用來估計群落中含有OTU 數目的指數，由Chao 提出，是生態學中估計物種總數的常用指數之一，與Chao1 的算法不同。
Simpson：用來估算樣品中微生物的多樣性指數之一，由Edward Hugh Simpson ( 1949) 提出，在生態學中常用來定量的描述一個區域的生物多樣性。Simpson 指數值越大，說明群落多樣性越高。
分別對 Alpha diversity 的各個指數進行秩和檢驗分析（若兩組樣品比較則使用 R 中的wilcox.test 函數，若兩組以上的樣品比較則使用 R 中的 kruskal.test 函數），通過秩和檢驗篩選不同條件下的顯著差異的 Alpha Diversity指數。
也許我們有聽說Beta多樣性在最近10年間成為生物多樣性研究的熱點問題之一。具體解釋下：
Beta多樣性度量時空尺度上物種組成的變化, 是生物多樣性的重要組成部分, 與許多生態學和進化生物學問題密切相關！
PCoA（principal co-ordinates analysis）是一種研究數據相似性或差異性的可視化方法，通過一系列的特徵值和特徵向量進行排序後，選擇主要排在前幾位的特徵值，PCoA 可以找到距離矩陣中最主要的坐標，結果是數據矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關係，只是改變了坐標系統。
每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構成差異越小。
主成分分析（Principal component analysis）PCA 是一種研究數據相似性或差異性的可視化方法，通過一系列的特徵值和特徵向量進行排序後，選擇主要的前幾位特徵值，採取降維的思想，PCA 可以找到距離矩陣中最主要的坐標，結果是數據矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關係，只是改變了坐標系統。
詳細關於主成分分析的解釋推薦大家看一篇文章，http://blog.csdn.net/aywhehe/article/details/5736659
研究背景：全球塑料產量飛速增長，而且呈持續上升的趨勢，因此導致大量塑料廢物排放到環境中，從沿海河口到大洋環流，從東大西洋到南太平洋海域。塑料廢棄物具有化學穩定性和生物利用率低的特點，可長期存在於海洋中，從而影響海洋環境包括海洋生物的生存。
作為一個獨特的底物，塑料碎片可以吸附海洋中的微生物並形成個「塑性球」。以生物膜形式存在於塑料碎片上的微生物群落。許多研究表明，無論是在海洋還是淡水生態系統中，附著在塑料碎片上微生物群落的組成明顯不同於周圍環境（水和沉積物），而且易受位置、時間和塑料類型的影響。
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具有較高的致病性。
NMDS（Nonmetric Multidimensional Scaling）常用於比對樣本組之間的差異，可以基於進化關係或數量距離矩陣。
橫軸和縱軸：表示基於進化或者數量距離矩陣的數值在二維表中成圖。與PCA分析的主要差異在於考量了進化上的信息。
每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構成差異越小。
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)
在非限制性排序中，16S和宏基因組數據分析通常用到的是PCA分析和PCoA分析，兩者的區別在於：
PCA分析是基於原始的物種組成矩陣所做的排序分析，而PCoA分析則是基於由物種組成計算得到的距離矩陣得出的。
在PCoA分析中，計算距離矩陣的方法有很多種，包括如：Euclidean, Bray-Curtis, and Jaccard，以及(un)weighted Unifrac (利用各樣品序列間的進化信息來計算樣品間距離，其中weighted考慮物種的豐度，unweighted沒有對物種豐度進行加權處理)。
如果說 PCA，它所作的只是將整組數據整體映射到最方便表示這組數據的坐標軸上，映射時沒有利用任何數據內部的分類信息，是無監督的。
那麼LDA是有監督的，增加了種屬之間的信息關係後，結合顯著性差異標準測試(克魯斯卡爾-沃利斯檢驗、兩兩Wilcoxon測試)和線性判別分析的方法進行特徵選擇。
兩者相同點：
都可以對數據進行降維。
降維時都採用了矩陣特徵分解的思想。
差異：
1）LDA是有監督學習的降維方法，而PCA是無監督的降維方法。（註：監督學習是從標記的訓練數據來推斷一個功能的機器學習任務。）
2）LDA選擇分類性能最好的投影方向，而PCA選擇樣本點投影具有最大方差的方向。
除了可以檢測重要特徵，他還可以根據效應值進行功能特性排序，這些功能特性可以解釋大部分生物學差異。這部分希望能詳細了解的同學可以參考這篇文章http://blog.csdn.net/sunmenggmail/article/details/8071502 。
組間差異顯著物種又可以稱作生物標記物（biomarkers），這個LDA分析主要是想找到組間在豐度上有顯著差異的物種。
研究背景：研究表明遺傳和環境影響都在I型糖尿病的發展中起作用，增加的遺傳風險不足以引起疾病，環境因素也是需要的，而且起著至關重要的作用。腸道菌群也許就是這個重要的環境因素，腸道菌群在免疫系統的成熟中起重要作用，此外還影響自身免疫疾病發展。
不同遺傳風險分組中包含的常見菌屬，部分存在特定分組中
PCoA分析揭示不同遺傳風險兒童腸道菌群的在不同地域樣本中均存在顯著差異
點評：針對I型糖尿病疾病發生過程中遺傳HLA分型風險和對應腸道菌群菌的關聯分析，揭示了特定腸道菌群與宿主特定遺傳風險共同作用推進疾病發生。某些特定菌屬可能無法在遺傳高風險兒童腸道內定植，可能對疾病發生存在特定作用。此外對於其他遺傳風險的自身免疫疾病也具有重要提示意義，例如乳糜瀉和類風溼性關節炎。
這是另一款和LDA長得有點像的圖，當然功能可完全不一樣。它是將不同樣本的群落構成及分布以物種分類樹的形式在一個環圖中展示。數據經過分析後，將物種分類樹和分類豐度信息通過這款軟體GraPhlAn進行繪製
(http://huttenhower.sph.harvard.edu/GraPhlAn )。
其目的是將物種之間的進化關係以及不同樣本的物種分布豐度和最高分布樣本的信息在一個視覺集中的環圖中一次展示，其提供的信息量較其他圖最為豐富。
中間為物種進化分類樹
不同顏色的分支代表不同的綱（具體的代表顏色見右上角的圖例），
接著的外圈的灰色標示字母的環表示的是本次研究中比例最高的15個科（字母代表的科參見左上角的圖例）。
之後的外圈提供的是熱力圖，如果樣本數<=10個則繪製樣本，如果樣本數超過10個則按照分組繪製，每一環為一個樣本，根據其豐度繪製的熱力圖。
最外圈為柱狀圖，繪製的是該屬所佔比例最高的樣本的豐度和樣本顏色（樣本顏色見環最下方的樣本名字的顏色）。其中熱力圖和柱狀圖取值均為原比例值x10000後進行log2轉換後的值。
根據各個物種在各個樣品中的豐度以及變化情況，計算物種之間的相關性，包括正相關和負相關。
首先對原始16s測序數據的種屬數量進行標準化，然後進行Spearman和Pearson秩相關分析並進行統計檢驗，計算出各個物種之間的相關性，之後在所有物種中根據simscore絕對值的大小，挑選出相關性最高的前100組數據，基於Cytoscape繪製共表達分析網絡圖。
○ 圖中每一個點代表一個物種，存在相關性的物種用連線連接。
○ 紅色的連線代表負相關，綠色的先代表正相關。
○ 連線顏色的深淺代表相關性的高低。
○ 圖中每一個點代表一個物種
○ 點的大小表示與其他物種的關聯關係的多少
○ 其中與之有相關性的物種數越多，點的半徑和字體越大
○ 連線的粗細代表兩物種之間相關性的大小
連線越粗，相關性越高。
研究背景：氣候變化導致美國中部草原的降水模式發生變化，對土壤微生物群落構成及代謝影響很大。
研究希望明確土壤微生物群落對土壤水分變化的反應，並確定響應的特定代謝特徵。
研究結論：土壤乾燥導致土壤微生物組的組成和功能發生顯著變化。相反，潤溼後幾乎沒有變化。由於乾旱導致的土壤水分減少對土壤碳循環和土壤微生物組進行的其他關鍵生物地球化學循環的影響很大。導致滲透保護劑化合物產生的代謝途徑受到較大影響。
點評：
相對簡單的樣本和實驗設計，但是從多個維度探尋支持土壤微生物群落對溼潤和乾燥表型的反應。
與常見的環境採樣檢測不同，針對同一樣本在對照環境下進行環境控制孵化，然後比較菌群變化可以更為有效的控制背景差異。
根據OTU數據進行標準化處理（1wlog10）之後，選取數目最多的前60個物種，基於R heatmap進行作圖
○ 熱圖中的每一個色塊代表一個樣品的一個屬的豐度
○ 樣品橫向排列，屬縱向排列
○ 差異是是否對樣品進行聚類，從聚類中可以了解樣品之間的相似性以及屬水平上的群落構成相似性。
Tips：
如果聚類結果中出現大面積的白或黑是因為大量的菌含量非常低，導致都沒有數值，可以在繪製之前進行標準化操作，對每一類菌單獨自身進行Z標準化。
研究背景：妊娠期糖尿病（GDM）的患病率在全球範圍內迅速增加，構成一個重要的健康問題和產科實踐的重大挑戰（Ferrara，2007）。高脂血症是妊娠常見的合併症。在GDM患者中，血脂的生理變化可能導致懷孕期間潛在的代謝紊亂。腸道失調在宿主代謝異常中起著至關重要的作用，最近關於2型糖尿病(T2D)和肥胖的研究就證明了這一點。這些研究表明，妊娠期間腸道微生物ME的主要變化可能在GDM的發展中起著至關重要的作用。
GDM加高脂血症（M隊列）妊娠期間與顯著改變的脂質相關的腸道微生物群（屬）
研究結論：我們的結果表明，血脂水平可能反映了GDM發展過程中的一些異常變化。所鑑定的多種生物標誌物對GDM合併高脂血症的防治有一定的參考價值。
組間物種差異性盒形圖描述在不同分組之間具有差異顯著的某一物種做箱線圖，圖中以屬水平為例做物種差異性箱線圖，展示如下：
○ 圖中不同顏色代表不同的分組，更直觀顯示組間物種差異
○ 每一個箱線圖代表一個物種，圖上方是物種名。
Anosim分析是一種非參數檢驗，用來檢驗組間的差異是否顯著大於組內差異，從而判斷分組是否有意義
R-value介於（-1，1）之間，R-value大於0，說明組間差異顯著。
R-value小於0,說明組內差異大於組間差異。
統計分析的可信度用 P-value 表示，P< 0.05 表示統計具有顯著性。
對Anosim的分析結果，基於兩兩樣本之間的距離值排序獲得的秩（組間的為between，組內的為within），這樣任一兩兩組的比較可以獲得三個分類的數據，並進行箱線圖的展示（若兩個箱的凹槽互不重疊，則表明它們的中位數有顯著差異）
隨機森林是機器學習算法的一種，它可以被看作是一個包含多個決策樹的分類器。
其輸出的分類結果是由每棵決策樹「投票」的結果。由於每棵樹在構建過程中都採用了隨機變量和隨機抽樣的方法，因此隨機森林的分類結果具有較高的準確度，並且不需要「減枝」來減少過擬合現象。
物種重要性點圖。橫坐標為重要性水平，縱坐標為按照重要性排序後的物種名稱。上圖反映了分類器中對分類效果起主要作用的菌屬，按作用從大到小排列。
Error rate: 表示使用下方的特徵進行隨機森林方法預測分類的錯誤率，越高表示基於菌屬特徵分類準確度不高，可能分組之間菌屬特徵不明顯。圖中以所有水平為例，取前60個作圖。
ROC 曲線指受試者工作特徵曲線(receiver operating characteristic curve), 是反映敏感性和特異性連續變量的綜合指標，通過構圖法揭示敏感性和特異性的相互關係。
ROC 曲線將連續變量設定出多個不同的臨界值，從而計算出一系列敏感性和特異性，再以敏感性為縱坐標、（1-特異性）為橫坐標繪製成曲線。
曲線下面積越大，診斷準確性越高。展示如下：
FAPROTAX是一款在2016年發表在SCIENCE上的較新的基於16S測序的功能預測軟體。它整合了多個已發表的可培養菌文章的手動整理的原核功能資料庫，資料庫包含超過4600個物種的7600多個功能注釋信息，這些信息共分為80多個功能分組，其中包括如硝酸鹽呼吸、產甲烷、發酵、植物病原等。
如果說PICRUSt（後續會介紹）在腸道微生物研究更為適合，那麼FAPROTAX尤其適用於生態環境研究，特別是地球化學物質循環分析。
FAPROTAX適用於對環境樣本（如海洋、湖泊等）的生物地球化學循環過程（特別是碳、氫、氮、磷、硫等元素循環）進行功能注釋預測。因其基於已發表驗證的可培養菌文獻，其預測準確度可能較好，但相比於上述PICRUSt和Tax4Fun來說預測的覆蓋度可能會降低。
FAPROTAX可根據16S序列的分類注釋結果對微生物群落功能（特別是生物地化循環相關）進行注釋預測。
圖中橫坐標代表樣本，縱坐標表示包括碳、氫、氮、硫等元素循環相關及其他諸多功能分組。可快速用於評估樣品來源或特徵。
Bugbase也是16年所提供服務的一款免費在線16S功能預測工具，到今年才發表文章公布其軟體原理。該工具主要進行表型預測，其中表型類型包括革蘭氏陽性、革蘭氏陰性、生物膜形成、致病性、移動元件、氧需求，包括厭氧菌、好氧菌、兼性菌）及氧化脅迫耐受等7類。
通過對已有測序微生物基因組的基因功能的構成進行分析後，我們可以通過16s測序獲得的物種構成推測樣本中的功能基因的構成，從而分析不同樣本和分組之間在功能上的差異（PICRUSt Nature Biotechnology, 1-10. 8 2013）。
通過對宏基因組測序數據功能分析和對應16s預測功能分析結果的比較發現，此方法的準確性在84%-95%，對腸道微生物菌群和土壤菌群的功能分析接近95%，能非常好的反映樣品中的功能基因構成。
為了能夠通過16s測序數據來準確的預測出功能構成，首先需要對原始16s測序數據的種屬數量進行標準化，因為不同的種屬菌包含的16s拷貝數不相同。
然後將16s的種屬構成信息通過構建好的已測序基因組的種屬功能基因構成表映射獲得預測的功能結果。（根據屬這個水平，對不同樣本間的物種豐度進行顯著性差異兩兩檢驗，我們這裡的檢驗方法使用STAMP中的two-sample中T-TEST方法，Pvalue值過濾為0.05，作Extent error bar圖。）
此處提供COG，KO基因預測以及KEGG代謝途徑預測。當然，躍躍欲試的小夥伴也可自行使用我們提供的文件和軟體（STAMP）對不同層級以及不同分組之間進行統計分析和製圖，以及選擇不同的統計方法和顯著性水平。
這裡提到的STAMP有些小夥伴說不太了解，別急，後面會有更多介紹。
圖中不同顏色代表不同的分組，列出了COG構成在組間存在顯著差異的功能分類以及在各組的比例，此外右側還給出了差異的比例和置信區間以及P-value。
通過KEGG代謝途徑的預測差異分析，我們可以了解到不同分組的樣品之間在微生物群落的功能基因在代謝途徑上的差異，以及變化的高低。為我們了解群落樣本的環境適應變化的代謝過程提供一種簡便快捷的方法。
本例圖所顯示的是第三層級的KEGG代謝途徑的差異分析，也可以針對第二或第一層的分級進行分析。
圖中不同顏色代表不同的分組，列出了在第三層級的構成在組間存在顯著差異的KEGG代謝途徑第三層分類以及在各組的比例，此外右側還給出了差異的比例和置信區間以及P-value。
研究背景：儘管普遍認為腸道微生物組的生態多樣性和分類組成在肥胖和T2D中發生改變，但與單個微生物或微生物產物的關聯在研究之間不一致。缺乏大樣本群體研究，從而確定腸道微生物組，血漿代謝組，肥胖和糖尿病表型以及環境因素之間的幾種關聯。
按照肥胖和糖尿病對人群分為三組，同時進行了16S，代謝和宏基因組的檢測。
研究結論：確定了腸道微生物組，血漿代謝組，肥胖和糖尿病表型以及環境因素之間的幾種關聯。與腸道微生物組變異相關的主要是肥胖，不是2型糖尿病。存在與腸道微生物組變異相關的藥物和膳食補充劑。高鐵攝入量影響小鼠的腸道微生物組成。微生物組變異也反映在血清代謝物譜中。
點評：
相對大人群的隊列研究，同時涵蓋了菌群、代謝和疾病表型以及膳食補充調查的數據。
從結果看菌屬和血漿代謝存在關聯，但是貢獻度都較低，如果樣本數量不足很可能找不到顯著的聯繫，這也是這類大樣本隊列研究的意義。
本研究在人群分組時針對性的研究了肥胖-II型糖尿病和菌群的關聯，因而構建了三個主要分組人群，結果顯示肥胖與菌群的關聯度更大，解釋了大部分的菌群差異，而糖尿病的菌群變化較小。
本研究其中較為重要的是發現了不同膳食補充對菌群的影響，並在小鼠實驗中得到證實。
除了能對大的基因功能分類和代謝途徑進行預測外，我們還能提供精細的功能基因的數量和構成的預測，以及進行樣本間以及組間的差異分析，並給出具有統計意義和置信區間的分析結果。
這一分析將我們對於樣本群落的差異進一步深入到了每一類基因的層面。
圖中不同顏色代表不同的分組，列出了在組間/樣本間存在顯著差異的每一個功能基因（酶）以及在各組的比例，此外右側還給出了差異的比例和置信區間以及P-value。
在獲得標準報告後如果希望單獨修改分組或對某些組之間進行顯著性差異分析，可以使用STAMP軟體在自己的電腦上進行數據分析。STAMP提供了豐富的統計檢驗方法和圖形化結果的輸出。
在使用STAMP之前需要首先準備需要的spf格式文件和樣品分組信息表，但是如果數據不會處理，那也很不便。
而在我們的報告中已經將KEGG和KO以及COG的結果文件後經過轉換生成了適用於STAMP軟體打開的spf格式文件，還有對應的分組信息表文件groupfile.txt。
STAMP 可以直接使用來自QIIME的biom文件和PICUST的KEGG和ko 文件，groupfile.txt文件的格式為tab-saperated value (tab鍵隔開的數據)
導入數據之後，view group legend ,在窗口右側會出現分組欄，根據需要進行分組。
3、Unclassiffied選項中，remain Unclassiffied reads、remove Unclassiffied reads、和use only for calculating frequency profiles 方法的區別？
remain Unclassiffied reads和use only for calculating frequency profiles方法會保留所有的數據，而remove Unclassiffied reads僅僅保留有確定分組信息的數據。
4、Statistical test 中，Welch’s t-test、t-test、white’s non-parametric t-test的區別，各自優缺點？
為了確保統計學意義和準確度和精確性，需要足夠多的樣本數目，t-test檢驗可以在最少樣本數為4的時候確保高的準確度和精確性。
當兩個樣本之間具有相同方差的時候，用t-test更為準確，當兩個樣本沒有相同方差，Welch’s t-test更為準確。
當樣本數目少於8的時候，可以使用white’s non-parametric t-test，該計算時間較長，當樣本數目過多的時候不宜使用該方法。
5、Two-group 中 type: one side和two side的區別？
One side 只會顯示前一個group與後一個group差異的比例，而two side 兩者之間的比例均會顯示。
6、STAMP在使用時首先打開了一個分析文件，如果新打開一個可能會導致顯示錯誤？
目前版本的STAMP存在一些小問題，一次分析只能使用一個數據文件，如果要打開新的需要關閉軟體後再打開。
詳細的STAMP使用教程可以參考我們提供的STAMP使用教程。
冗餘分析（redundancy analysis, RDA）或者
典範對應分析（canonical correspondence analysis, CCA）都是基於對應分析發展的一種排序方法，將對應分析與多元回歸分析相結合，每一步計算均與環境因子進行回歸，又稱多元直接梯度分析。主要用來反映菌群與環境因子之間的關係。
RDA 是基於線性模型，CCA是基於單峰模型。分析可以檢測環境因子、樣品、菌群三者之間的關係或者兩兩之間的關係。
○ 冗餘分析可以基於所有樣品的OTU作圖，也可以基於樣品中優勢物種作圖；
○ 箭頭射線：箭頭分別代表不同的環境因子；
○ 夾角：環境因子之間的夾角為銳角時表示兩個環境因子之間呈正相關關係，鈍角時呈負相關關係。環境因子的射線越長，說明該影響因子的影響程度越大；
○ 不同顏色的點表示不同組別的樣品或者同一組別不同時期的樣品，圖中的拉丁文代表物種名稱，可以將關注的優勢物種也納入圖中；
○ 環境因子數量要少於樣本數量，同時在分析時，需要提供環境因子的數據，比如 pH值，測定的溫度值等。
除以上部分，還可以進行個性化圖表定製，像下面這樣：
看完以上內容，也許還有不明白的地方，沒關係，我們羅列了一些常見的問題。看看有沒有你想問的。
原始fastq格式是一個文本格式用於存儲生物序列（通常是核酸序列）和其測序對應的質量值。這些序列以及質量信息用ASCII字符標識。通常fastq文件中一個序列有4行信息：如
第一行：序列標識，以 @開頭。格式自由，允許添加描述信息，描述信息以空格分開。
第二行：序列信息，不允許出現空格或制表符。一般是明確的DNA或RNA字符，通常大寫
第三行：用於將序列信息和質量值分隔開。以 +開頭，後邊是描述信息或者不加。
第四行：質量值， 每個字符與第二行的鹼基一一對應，按照一定規則轉換為鹼基質量得分。進而反映該鹼基的錯誤率，因此字符數必須和第二行保持一致。
Fasta格式
fasta是一種基於文本用於表示核苷酸序列或胺基酸序列的格式。在這種格式中鹼基對或胺基酸用單個字母來編碼，且允許在序列前添加序列名及注釋。由兩部分信息組成：如
第一行：序列標記，以 >開頭，接序列的標識符，序列標識符以空格結束，後接描述信息。為保證分析軟體能區分每條序列，每個序列的標識必須具有唯一性。
第二行：序列信息，使用既定的核苷酸或胺基酸編碼符號。
原始數據（Raw data），常見的是illumina機器產生的fastq文件，這一類文件需要向NCBI的SRA資料庫進行提交，SRA是NCBI為了並行測序的高通量數據（massively parallel sequencing）提供的存儲平臺。完整提交SRA需要一些獨立項目的分步提交，包括BioProject、BioSample、Experiment、Run等，每一部分用以描述數據的不同屬性。
原始的Tags數據會經過質控、過濾、去嵌合體，最終得到有效數據（Effective Tags）。所以在判斷測序質量是否合格時應該從幾個方面去判斷。
打開文件01_sequence_statistic/sumOTUPerSample.txt
報告裡所有的txt打開如果格式不對的話，可以用excel表打開。
其中tags為經質量過濾後能正確overlap包含正確barcode和高質量序列的數據。
Singleton為非完全相同的序列，只要有1個鹼基的差異即為不同序列，該值的高低與OUT數量並無直接關係，OTU是以97%的相似度聚類，測序質量較低導致的鹼基錯誤、PCR擴增過程中的鹼基錯誤、菌種內部的多樣性以及OTU數量均會影響該數量。
Chimeras為通過與RDP等標準資料庫比對分析判斷可能由於PCR過程錯誤擴增導致的嵌合體比例，chimeras%為百分比，一般低於1。
首先判斷下機數據tags和有效數據 clean tags 的數據量是否滿足測序要求，一般下機數據量達到3萬條reads以上滿足測序需要，谷禾16s樣本的測序深度可以達到10萬條reads左右。如果數據量不夠則需要重新補測樣本。通過觀察嵌合體數chimras 和嵌合體所佔百分比chimeras%，可以反應出有效序列的轉化率，嵌合體的比例越小序列的利用轉化率就越高。
根據稀釋曲線可以判斷測序深度是否達到飽和，如圖中曲線都逐漸趨於平緩，就證明樣本的測序深度較好，測序深度基本覆蓋能測到的該樣本所有的物種，測序深度比較好。同時曲線趨於水平縱坐標的高低也能夠反映各樣本的微生物多樣性情況，曲線越高，證明測到的物種種類越多，樣本的微生物多樣性就越高。
而從該圖可以看出，個別樣本的曲線未趨於平緩，證明該樣本測序深度不夠，測序深度未能很好的反映出該樣本的完整菌群構成。如果測序數據量更大的的話會檢測到更多物種。
如何了解分組內部的多個樣本的重複性以及多樣性情況？
觀察分組內部多個樣本的重複性如何可以從以下幾個方面考慮。
首先在各分類水平的柱狀圖的菌屬構成來看
從構成圖來看，Flu組和ZW3.7組，組內樣本重複性較好。Ctrl組中Ctrl.2明顯區別於組內另外兩個樣本，可以去掉該樣本。而ZW3.8組內樣本間差異性較大。
比如人體腸道或小鼠腸道樣本本身個體差異性較大，菌群結構組成複雜，即便通過不同疾病的分類的樣本，但營養飲食、代謝以及環境的影響都會改變腸道菌群的構成，所以有可能組內樣本間差異性會比較大。而經過單因素處理的樣本組內差異會比較小。
所以在前期實驗設計時，儘量選擇同一批次相同處理的小鼠或其他樣本，避免組內差異的影響。並且要預留好多餘的樣本，比如組內只有3個樣本，如果去掉一個差異性較大的樣本，一個分組內只有2個樣本，會影響後續組間差異比較，組間差異性比較分析每組要至少要3個樣本。
通過beta多樣性分析PCA,PCoA,MNDS 也可以大致觀察組內樣本重複性情況，左圖組內樣本重複性較好，右圖組內樣本間差異性較大，兩組間的區割不是很明顯。
在加圈圖的beta多樣性分析中，右下角有給出PC1和PC2的P值，小於0.05則差異顯著。
Alpha多樣性是針對單個樣品中物種多樣性的分析，包括chao1指數、ace指數，shannon指數以及simpson指數等。前面4個指數越大，最後一個指數越小，說明樣品中的物種越豐富。
其中chao指數和ACE指數反映樣品中群落的豐富度（species richness），即簡單指群落中物種的數量，而不考慮群落中每個物種的豐度情況。指數對應的稀釋曲線還可以反映樣品測序量是否足夠。如果曲線趨於平緩或者達到平臺期時也就可以認為測序深度已經基本覆蓋到樣品中所有的物種；反之，則表示樣品中物種多樣性較高，還存在較多未被測序檢測到的物種。
而shannon指數以及simpson指數反映群落的多樣性（species diversity），受樣品群落中物種豐富度（species richness）和物種均勻度（species evenness）的影響。相同物種豐富度的情況下，群落中各物種具有越大的均勻度，則認為群落具有越大的多樣性。
稀釋曲線是利用已測得序列中已知的各種OTU的相對比例，來計算抽取n個（n小於測得Reads序列總數）Tags時各Alpha指數的期望值，然後根據一組n值（一般為一組小於總序列數的等差數列，本項目公差為500 ）與其相對應的Alpha指數的期望值繪製曲線。
不同的樣本之間差異大嗎？不同分組之間能否用菌群差異來區分？
觀察不同分組間差異的大小可以觀察隨機森林分類效果圖。
路徑在07_diff_analysis/RF
圖中以該分類水平下選取用於區分不同分組間的差異性起到關鍵性影響因素的物種作為標誌物作圖。標誌物按重要性從大到小排列，圖中隨機森林值error rate 表示用隨機森林方法預測分組之間的錯誤率，分值越高代表所選取的標誌物準確度不高，並不能很好的用於區分各分組，分組差異不顯著。分值越低證明分組效果比較好。
上圖中的隨機森林按照門和屬以及代謝途徑分別進行分析作圖，各自都有單獨文件，報告中僅給出了一個圖，其他文件需要到目錄中查看。可能存在門或屬區分效果不佳，但是代謝途徑區分效果較好。
隨機森林篩選出來的物種是用於區分所有分組的重要標誌。分值越高代表該物種用於區分所有組之間的重要性越大。
16s測序主要為了鑑定菌種，通常在做鑑定的時候區分標準是97%，區分亞種和菌株的時候相似度更高。
普通TAQ酶的複製錯誤率較高，可能在擴增過程中引入錯誤，這些錯配可能導致相似度下降從而分類錯誤。
一般我們不建議使用普通TAQ酶進行擴增，都選擇高保真酶。
16s rRNA長度為1.5k多，作為菌種鑑定一般選擇相似度97%的標準，相似度超過97%一般定義為同一種菌。
如果是sanger測序獲得16s全長的都可以鑑定到種，甚至能區分亞種。有些細菌並不只有1個16s序列，會包含有1-15拷貝的16s序列，所以單一的16s序列鑑定可能會出現偏差。
利用高通量如454或miseq測序一般由於讀長的緣故，通常只有300-500多個鹼基被測序，所以在物種鑑定上一般比較可靠的是能分類到屬，部分能分類到種。
根據我們的經驗，不同的樣品會有大約10-50的菌能分類到種。利用新的分析方法，我們現在也可以利用16s rRNA的群落多樣性高通測序數據進行亞種級別的分析。主要是利用16s中共同變化的SNP位點進行分型。這樣可以大大提高菌種的分類精度，尤其是在有些菌株之間表型差異巨大的時候。
16s序列能夠區分菌的種屬，但是並不包含這些菌的基因和代謝功能的信息。不過由於我們已經對大量的細菌基因組進行了測序，所以可以根據16s的菌種信息，利用這個菌屬已經測序的細菌基因組的基因信息和代謝功能信息來估計每類基因的上限和下限。
所以答案是可以利用16s序列測序來預測菌群的功能基因分布和代謝途徑分布情況。
目前主要使用的軟體是PICRUSt和新發表的Tax4Fun。
從我們實際分析和實驗結果來看，預測的準確性還是很高的，不過和樣品有很大關係。像腸道菌群和土壤以及一些致病菌的測序較多，所以預測的準確度較高可以到85-90%以上。一些海洋的菌由於測序的菌較少，預測準確性要差一些。目前發表的文獻基本都是用PICRUSt，新的軟體還有待驗證。
測16s rRNA能分到亞種嗎？不同菌株都有致病性差異光到種不解決問題啊！
16s rRNA如果是使用sanger測序可以細分到亞種甚至有些可以精確區分菌株，但是要看菌種。
如果是高通量測序，目前的常見分析一般以97%為標準，大部分情況只能到屬，少部分能區分到種。如果要進一步細分到亞種甚至更小的區分目前是有可能的，我們在使用oligotype一類的方法時可以將相同變化模式的SNP歸類，並對原來的OTU進行進一步細分，理論上可以區分到菌株。
不過這種區分不同菌屬差異很大，有些可以很理想的區分，主要用來了解在更細分化尺度上菌株構成的地理和時間變化。
僅通過16s高通量測序恐怕不能完全解決菌株致病性差異這種問題，但是通過對常見OTU的進一步深入分析可以提供可能的解釋或方向。如果明確了某一特定類型菌株的變化有關，可以採用比如毒力基因或菌株特異性標記等方法詳細了解不同菌株的比例和差異。
多項合作成果發表於Nature communications、PNAS、Plant biotechnology journal、DNA Research、Environmental Science & Technology、Plant、cell & environment、Science of The Total Environment 、Gut Microbes 、Frontiers in microbiologyt、Journal of environmental management 等國際著名學術期刊。
Li W, Zhang Y, Wu N, et al. Colonization characteristics of bacterial communities on plastic debris, influenced by environmental factors and polymer types in the Haihe Estuary of Bohai Bay, China[J]. Environmental Science & Technology, 2019.
Mao Z, Li Y, Dong T, et al. Exposure to Titanium Dioxide Nanoparticles During Pregnancy Changed Maternal Gut Microbiota and Increased Blood Glucose of Rat[J]. Nanoscale research letters, 2019, 14: 26.
He Z, Kong X, Shao T, et al. Alterations of the Gut Microbiome Associated with Promoting Efficacy of Prednisone by Bromofuranone in MRL/lpr Mice[J]. Frontiers in Microbiology, 2019, 10: 978.
Liu Y, Lu M, Zhang X, et al. Shift of the microbial communities from exposed sandstone rocks to forest soils during pedogenesis[J]. International Biodeterioration & Biodegradation, 2019, 140: 21-28.
Lu H, Wang W, Li F, et al. Mixed-surfactant-enhanced phytoremediation of PAHs in soil: Bioavailability of PAHs and responses of microbial community structure[J]. Science of The Total Environment, 2019, 653: 658-666.
Xu S, Wang W, Zhu L. Enhanced microbial degradation of benzo [a] pyrene by chemical oxidation[J]. Science of The Total Environment, 2019, 653: 1293-1300.
Ji C, Yan L, Chen Y, et al. Evaluation of the developmental toxicity of 2, 7-dibromocarbazole to zebrafish based on transcriptomics assay[J]. Journal of hazardous materials, 2019, 368: 514-522.
Ju M, Liu Y, Li M, et al. Baicalin improves intestinal microecology and abnormal metabolism induced by high-fat diet[J]. European journal of pharmacology, 2019: 172457.
Wu H, Li Y, Zhang W, et al. Bacterial community composition and function shift with the aggravation of water quality in a heavily polluted river[J]. Journal of environmental management, 2019, 237: 433-441.
Xia H, Wu Y, Chen X, et al. Effects of antibiotic residuals in dewatered sludge on the behavior of ammonia oxidizers during vermicomposting maturation process[J]. Chemosphere, 2019, 218: 810-817.
Sun Y, Geng W, Pan Y, et al. Supplementation with Lactobacillus kefiranofaciens ZW3 from Tibetan Kefir improves depression-like behavior in stressed mice by modulating the gut microbiota[J]. Food & function, 2019, 10: 925-937.
Wu S, Qiu J, Gao Q. QTL-BSA: A Bulked Segregant Analysis and Visualization Pipeline for QTL-seq[J]. Interdisciplinary Sciences: Computational Life Sciences, 2019: 1-8
Liu H, Pan L, Lv S, et al. Alterations of Gut Microbiota and Blood Lipidome in Gestational Diabetes Mellitus with Hyperlipidemia[J]. Frontiers in Physiology, 2019, 10: 1015.
我們除了做科研服務之外，還有對人體腸道菌群構成和個人健康狀態的檢測服務，目前已上線的產品如下：
如有合作意向或者各種疑問，歡迎打擾。
方式1：撥打商務經理電話（微信同號）：13336028502
方式2：撥打穀禾熱線：400-161-1580
無論何種方式，都有專門的工作人員與您交流探討