微生物群落多樣性測序與功能分析

微生物群落多樣性測序與功能分析微生物群落測序是指對微生物群體進(jìn)行高通量測序，通過分析測序序列的構(gòu)成分析特定環(huán)境中微生物群體的構(gòu)成情況或基因的組成以及功能。借助不同環(huán)境下微生物群落的構(gòu)成差異分析我們可以分析微生物與環(huán)境因素或宿主之間的關(guān)系，尋找標(biāo)志性菌群或特定功能的基因。對微生物群落進(jìn)行測序包括兩類，一類是通過16s rDNA，18s rDNA，ITS區(qū)域進(jìn)行擴(kuò)增測序分析微生物的群體構(gòu)成和多樣性；還有一類是宏基因組測序，是不經(jīng)過分離培養(yǎng)微生物，而對所有微生物DNA進(jìn)行測序，從而分析微生物群落構(gòu)成，基因構(gòu)成，挖掘有應(yīng)用價值的基因資源。以16s rDNA擴(kuò)增進(jìn)行測序分析主要用于微生物群落多樣性和構(gòu)成的分析，目前的生物信息學(xué)分析也可以基于16s rDNA的測序?qū)ξ⑸锶郝涞幕驑?gòu)成和代謝途徑進(jìn)行預(yù)測分析，大大拓展了我們對于環(huán)境微生物的微生態(tài)認(rèn)知。目前我們根據(jù)16s的測序數(shù)據(jù)可以將微生物群落分類到種（species）（一般只能對部分菌進(jìn)行種的鑒定），甚至對亞種級別進(jìn)行分析，幾個概念：16S rDNA（或16S rRNA）：16S rRNA基因是編碼原核生物核糖體小亞基的基因，長度約為1542bp，其分子大小適中，突變率小，是細(xì)菌系統(tǒng)分類學(xué)研究中最常用和最有用的標(biāo)志。16S rRNA基因序列包括9個可變區(qū)和10個保守區(qū)，保守區(qū)序列反映了物種間的親緣關(guān)系，而可變區(qū)序列則能體現(xiàn)物種間的差異。16S rRNA基因測序以細(xì)菌16S rRNA基因測序為主,核心是研究樣品中的物種分類、物種豐度以及系統(tǒng)進(jìn)化。OTU：operational taxonomic units (OTUs)在微生物的免培養(yǎng)分析中經(jīng)常用到，通過提取樣品的總基因組DNA，利用16S rRNA或ITS的通用引物進(jìn)行PCR擴(kuò)增，通過測序以后就可以分析樣品中的微生物多樣性，那怎么區(qū)分這些不同的序列呢，這個時候就需要引入operational taxonomic units，一般情況下，如果序列之間，比如不同的 16S rRNA序列的相似性高于97%就可以把它定義為一個OTU，每個OTU對應(yīng)于一個不同的16S rRNA序列，也就是每個OTU對應(yīng)于一個不同的細(xì)菌（微生物）種。通過OTU分析，就可以知道樣品中的微生物多樣性和不同微生物的豐度。測序區(qū)段：由于16s rDNA較長（1.5kb），我們只能對其中經(jīng)常變化的區(qū)域也就是可變區(qū)進(jìn)行測序。16s rDNA包含有9個可變區(qū)，分別是v1-v9。一般我們對v3-v4雙可變區(qū)域進(jìn)行擴(kuò)增和測序，也有對v1-v3區(qū)進(jìn)行擴(kuò)增測序。工具/原料 16s rDNA測序首先需要提取環(huán)境樣品的DNA，這些DNA可以來自土壤、糞便、空氣或水體等任何來源。 提取DNA后需要經(jīng)過質(zhì)檢和純化，一般16s rDNA測序擴(kuò)增對DNA的總量要求并不高，總量大于100ng，濃度大于10ng/ul一般都可以滿足要求。如果是來自和寄主共生的環(huán)境如昆蟲的腸道微生物，提取時可能包括了寄主本身的大量DNA，對DNA的總量要求會提高。微生物菌群多樣性測序受DNA提取和擴(kuò)增影響很大，不同的擴(kuò)增區(qū)段和擴(kuò)增引物甚至PCR循環(huán)數(shù)的差異都會對結(jié)果有所影響。因而建議同一項目不同樣品的都采用相同的條件和測序方法，這樣相互之間才存在可比性。 完成PCR之后的產(chǎn)物一般可以直接上測序儀測序，在上機(jī)測序前我們需要對所有樣本進(jìn)行定量和均一化，通常要進(jìn)行熒光定量PCR。完成定量的樣品混合后就可以上機(jī)測序。 16s rDNA測序目前可以采用多種不同的測序儀進(jìn)行測序，包括羅氏的454，Illumina的MiSeq，Life的PGM或Pacbio的RSII三代測序儀。不同的儀器各有優(yōu)缺點，目前最主流的是Illumina公司的MiSeq，因為其在通量、長度和價格三者之間最為平衡。MiSeq測序儀可以產(chǎn)生2x300bp的測序讀長，一次可以產(chǎn)生15Gb的測序數(shù)據(jù)遠(yuǎn)遠(yuǎn)大于其他測序儀的測序通量。方法/步驟1. 116s rDNA分析基本流程：2. 2原始數(shù)據(jù)處理：原始測序數(shù)據(jù)需要去除接頭序列，并將雙端測序序列進(jìn)行拼接成單條序列。根據(jù)測序barcode序列區(qū)分不同的樣本序列。過濾低質(zhì)量序列和無法比對到16s rDNA數(shù)據(jù)庫的序列。3. 3OTU分類和統(tǒng)計：OTU(operational taxonomic units) 是在系統(tǒng)發(fā)生學(xué)研究或群體遺傳學(xué)研究中，為了便于進(jìn)行分析，人為給某一個分類單元（品系，種，屬，分組等）設(shè)置的同一標(biāo)志。通常按照 97% 的相似性閾值將序列劃分為不同的 OTU，每一個 OTU 通常被視為一個微生物物種。相似性小于97%就可以認(rèn)為屬于不同的種，相似性小于93%-95%，可以認(rèn)為屬于不同的屬。樣品中的微生物多樣性和不同微生物的豐度都是基于對OTU的分析。使用QIIME（version 1.8.0）工具包進(jìn)行統(tǒng)計注釋。使用QIIME（version 1.9.0, http:/bio.cug.edu.cn/qiime/）的ucluster方法根據(jù)97%的序列相似度將所有序列進(jìn)行同源比對并聚類成operational taxonomic units (OTUs)。然后與數(shù)據(jù)庫GreenGenes（version gg_13_8, http:/greengenes.lbl.gov/cgi-bin/JD_Tutorial/nph-16S.cgi）進(jìn)行比對，比對方法uclust，identity 0.9 。然后對每個OTUs進(jìn)行reads數(shù)目統(tǒng)計。下面的2個表，其中一個表是對每個樣本的測序數(shù)量和OTU數(shù)目進(jìn)行統(tǒng)計，并且在表栺中列出了測序覆蓋的完整度（顯示前10個樣本）。另一個表是對每個樣本在分類字水平上的數(shù)量進(jìn)行統(tǒng)計，并且在表栺中列出了在每個分類字水平上的物種數(shù)目（顯示前10個樣本）?？梢钥吹浇^大部分的OTU都分類到了屬（Genus），也有很多分類到了種（Species）。但是仍然有很多無法完全分類到種一級，這是由于環(huán)境微生物本身存在非常豐富的多樣性，還有大量的菌仍然沒有被測序和發(fā)現(xiàn)。測序數(shù)目統(tǒng)計表主要是對每個樣本的測序數(shù)量和OTU數(shù)目進(jìn)行統(tǒng)計，并且在表格中列出了測序覆蓋的完整度（顯示前10個樣本，如果樣本超過10個，請查看結(jié)果中otu_stat.txt文件）其中 SampleName表示樣本名稱；SampleSize表示樣本序列總數(shù)；OTUsNumber表示注釋上的OTU數(shù)目；OTUsSeq表示注釋上OTU的樣本序列總數(shù)。Coverage是指各樣品文庫的覆蓋率，其數(shù)值越高，則樣本中序列沒有被測出的概率越低。該指數(shù)實際反映了本次測序結(jié)果是否代表樣本的真實情況。計算公式為：C=1-n1/N 其中n1 = 只含有一條序列的OTU的數(shù)目； N = 抽樣中出現(xiàn)的總的序列數(shù)目。分類水平統(tǒng)計表主要是對每個樣本在分類學(xué)水平上的數(shù)量進(jìn)行統(tǒng)計，并且在表格中列出了在每個分類學(xué)水平上的物種數(shù)目（只顯示前10個樣本，如果樣本超過10個，請查看結(jié)果中taxon_all.txt文件）其中SampleName表示樣本名稱；Phylum表示分類到門的OTU數(shù)量；Class表示分類到綱的OTU數(shù)量；Order表示分類到目的OTU數(shù)量；Family表示分類到科的OTU數(shù)量；Genus表示分類到屬的OTU數(shù)量；Species表示分類到種的OTU數(shù)量。4. 4我們還可以對這些種屬的構(gòu)成進(jìn)行柱狀圖顯示：橫坐標(biāo)中每一個條形圖代表一個樣本，縱坐標(biāo)代表該分類層級的序列數(shù)目或比例。同一種顏色代表相同的分類級別。圖中的每根柱子中的顏色表示該樣本在不同級別（門、綱、目等）的序列數(shù)目，序列數(shù)目只計算級別最低的分類，例如在屬中計算過了，則在科中則不重復(fù)計算。Q: 為什么要選擇V3-V4區(qū)的測序長度？為什么有些文獻(xiàn)是V6區(qū)，有什么區(qū)別？A: 16S rRNA總長約1540 bp，包含 9個可變區(qū)。由于高通量測序的測序長度的限制，不可能將16S rRNA的9個可變區(qū)全部測序，所以在PCR擴(kuò)增時往往只能選擇1-3個可變區(qū)作為擴(kuò)增片段。Kozich 等評估了Miseq測序儀分析的不同16S rRNA可變區(qū)的準(zhǔn)確性發(fā)現(xiàn)，測定 V4 區(qū)效果最佳。根據(jù)我們的測序長度，v3-v4區(qū)是最佳選擇。5. 5我們還需要對樣本之間或分組之間的OTU進(jìn)行比較獲得韋恩圖：注意，韋恩圖目前一般最多只能顯示5個樣本或分組，過多的樣本無法無法進(jìn)行韋恩圖繪制6. 6樣品構(gòu)成豐度：稀釋曲線微生物多樣性分析中需要驗證測序數(shù)據(jù)量是否足以反映樣品中的物種多樣性，稀釋曲線（豐富度曲線）可以用來檢驗這一指標(biāo)。稀釋曲線是用來評價測序量是否足以覆蓋所有類群，并間接反映樣品中物種的豐富程度。稀釋曲線是利用已測得16S rDNA序列中已知的各種OTU的相對比例，來計算抽取n個（n小于測得reads序列總數(shù)）reads時出現(xiàn)OTU數(shù)量的期望值，然后根據(jù)一組n值（一般為一組小于總序列數(shù)的等差數(shù)列）與其相對應(yīng)的OTU數(shù)量的期望值做出曲線來。當(dāng)曲線趨于平緩或者達(dá)到平臺期時也就可以認(rèn)為測序深度已經(jīng)基本覆蓋到樣品中所有的物種；反之，則表示樣品中物種多樣性較高，還存在較多未被測序檢測到的物種。下圖中的稀釋曲線橫坐標(biāo)代表隨機(jī)抽取的序列數(shù)量；縱坐標(biāo)代表觀測到的OTU數(shù)量。樣本曲線的延伸終點的橫坐標(biāo)位置為該樣本的測序數(shù)量，如果曲線趨于平坦表明測序已趨于飽和，增加測序數(shù)據(jù)無法再找到更多的OTU；反之表明不飽和，增加數(shù)據(jù)量可以發(fā)現(xiàn)更多OTU。7. 7Shannon-Winner曲線Shannon-Wiener 曲線，是利用shannon指數(shù)來進(jìn)行繪制的，反映樣品中微生物多樣性的指數(shù)，利用各樣品的測序量在不同測序深度時的微生物多樣性指數(shù)構(gòu)建曲線，以此反映各樣本在不同測序數(shù)量時的微生物多樣性。 當(dāng)曲線趨向平坦時，說明測序數(shù)據(jù)量足夠大，可以反映樣品中絕大多數(shù)的微生物物種信息。與上圖一樣，橫坐標(biāo)代表隨機(jī)抽取的序列數(shù)量；縱坐標(biāo)代表的是反映物種多樣性的Shannon指數(shù)。樣本曲線的延伸終點的橫坐標(biāo)位置為該樣本的測序數(shù)量，如果曲線趨于平坦表明測序已趨于飽和，增加測序數(shù)據(jù)無法再找到更多的OTU；反之表明不飽和，增加數(shù)據(jù)量可以發(fā)現(xiàn)更多OTU。其中曲線的最高點也就是該樣本的Shannon指數(shù)，指數(shù)越高表明樣品的物種多樣性越高。Q: Shannon指數(shù)怎么算的？A: Shannon指數(shù)公式：其中，Sobs=實際測量出的OTU數(shù)目；ni=含有i 條序列的OTU數(shù)目；N=所有的序列數(shù)。8. 8Rank-Abundance曲線用于同時解釋樣品多樣性的兩個方面，即樣品所含物種的豐富程度和均勻程度。物種的豐富程度由曲線在橫軸上的長度來反映，曲線越寬，表示物種的組成越豐富；物種組成的均勻程度由曲線的形狀來反映，曲線越平坦，表示物種組成的均勻程度越高。一般超過20個樣本圖就會變得非常復(fù)雜而且不美觀，所以一般20個樣本以下會做該圖，圖片保存為結(jié)果目錄中rank.pdf。橫坐標(biāo)代表物種排序的數(shù)量；縱坐標(biāo)代表觀測到的相對豐度。樣本曲線的延伸終點的橫坐標(biāo)位置為該樣本的物種數(shù)量，如果曲線越平滑下降表明樣本的物種多樣性越高，而曲線快速陡然下降表明樣本中的優(yōu)勢菌群所占比例很高，多樣性較低。9. 9Alpha多樣性（樣本內(nèi)多樣性）Alpha多樣性是指一個特定區(qū)域或者生態(tài)系統(tǒng)內(nèi)的多樣性，常用的度量指標(biāo)有Chao1 豐富度估計量（Chao1 richness estimator） 、香農(nóng) - 威納多樣性指數(shù)（Shannon-wiener diversity index）、辛普森多樣性指數(shù)（Simpson diversity index）等。計算菌群豐度：Chao、ace；計算菌群多樣性：Shannon、Simpson。Simpson指數(shù)值越大，說明群落多樣性越高；Shannon指數(shù)越大，說明群落多樣性越高。表中顯示前10個樣本，如果樣本大于10個，詳見結(jié)果目錄中的alpha_div.txt。Q: 能不能解釋下每個指數(shù)（如chao1、shannon）？A: Chao1：是用chao1 算法估計群落中含OTU 數(shù)目的指數(shù)，chao1 在生態(tài)學(xué)中常用來估計物種總數(shù)，由Chao (1984) 最早提出。Chao1值越大代表物種總數(shù)越多。Schao1=Sobs+n1(n1-1)/2(n2+1)其中Schao1為估計的OTU數(shù)，Sobs為觀測到的OTU數(shù)，n1為只有一條序列的OTU數(shù)目，n2為只有兩條序列的OTU數(shù)目。Shannon：用來估算樣品中微生物的多樣性指數(shù)之一。它與 Simpson 多樣性指數(shù)均為常用的反映 alpha 多樣性的指數(shù)。Shannon值越大，說明群落多樣性越高。Ace：用來估計群落中含有OTU 數(shù)目的指數(shù)，由Chao 提出，是生態(tài)學(xué)中估計物種總數(shù)的常用指數(shù)之一，與Chao1 的算法不同。Simpson：用來估算樣品中微生物的多樣性指數(shù)之一，由Edward Hugh Simpson ( 1949) 提出，在生態(tài)學(xué)中常用來定量的描述一個區(qū)域的生物多樣性。Simpson 指數(shù)值越大，說明群落多樣性越高。辛普森多樣性指數(shù)=隨機(jī)取樣的兩個個體屬于不同種的概率=1-隨機(jī)取樣的兩個個體屬于同種的概率10. 10Beta多樣性分析（樣品間差異分析）Beta多樣性度量時空尺度上物種組成的變化,是生物多樣性的重要組成部分,與許多生態(tài)學(xué)和進(jìn)化生物學(xué)問題密切相關(guān),因此在最近10年間成為生物多樣性研究的熱點問題之一。PCoA分析PCoA（principal co-ordinates analysis）是一種研究數(shù)據(jù)相似性或差異性的可視化方法，通過一系列的特征值和特征向量進(jìn)行排序后，選擇主要排在前幾位的特征值，PCoA 可以找到距離矩陣中最主要的坐標(biāo)，結(jié)果是數(shù)據(jù)矩陣的一個旋轉(zhuǎn)，它沒有改變樣品點之間的相互位置關(guān)系，只是改變了坐標(biāo)系統(tǒng)。通過PCoA 可以觀察個體或群體間的差異。每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構(gòu)成差異越小。PCoA有多張圖，分別代表的PCoA1-2,2-3,3-1。11. 11NMDS分析（非度量多維尺度分析）NMDS（Nonmetric Multidimensional Scaling）常用于比對樣本組之間的差異，可以基于進(jìn)化關(guān)系或數(shù)量距離矩陣。橫軸和縱軸：表示基于進(jìn)化或者數(shù)量距離矩陣的數(shù)值 在二維表中成圖。與PCA分析的主要差異在于考量了進(jìn)化上的信息。每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構(gòu)成差異越小。12. 12PCA分析主成分分析PCA（Principal component analysis）是一種研究數(shù)據(jù)相似性或差異性的可視化方法，通過一系列的特征值和特征向量進(jìn)行排序后，選擇主要的前幾位特征值，采取降維的思想，PCA 可以找到距離矩陣中最主要的坐標(biāo)，結(jié)果是數(shù)據(jù)矩陣的一個旋轉(zhuǎn)，它沒有改變樣品點之間的相互位置關(guān)系，只是改變了坐標(biāo)系統(tǒng)。詳細(xì)關(guān)于主成分分析的解釋推薦大家看一篇文章，http:/blog.csdn.net/aywhehe/article/details/5736659 。通過PCA 可以觀察個體或群體間的差異。每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構(gòu)成差異越小。以上三個圖可能遇到的問題：1：PCA，PcoA，NMDS分析分別是基于什么數(shù)據(jù)畫的？回答：PCA，PcoA，NMDS分析均是基于OTU分類taxon數(shù)據(jù)所畫，用的是R語言Vegan包中的相關(guān)函數(shù)畫成，其中PcoA與NMDS還要基于樣本之間的距離矩陣才能畫成。2：PCA分析如果圖中大部分點集中在一起，少數(shù)點在很遠(yuǎn)的外圍，是什么原因造成的？回答：是因為樣本OTU分類時候，少數(shù)樣本某些菌含量特別高所造成，導(dǎo)致這些樣本偏離正常范圍，建議單獨拿出這些樣本觀察，看是否是實驗錯誤。3：PCA分析時，不是有PC1，PC2，PC3三個坐標(biāo)嗎？是給出三張圖嗎？還是三維立體圖？回答：PCA作圖時，會得出PC1，PC2，PC3三個坐標(biāo)，可以根據(jù)PC12,PC13,PC23分別作圖，一般給出的是PC12的圖，當(dāng)PC12圖質(zhì)量不好，看不出明顯的樣本分類效果時，可以看PC13或PC23的圖分類是否清晰，也可以用R語言rgl包做出PC123三維圖。QIIME本身結(jié)果中有提供PCA的三維圖結(jié)果，可以通過網(wǎng)頁打開。13. 13LDA差異貢獻(xiàn)分析PCA和LDA的差別在于，PCA，它所作的只是將整組數(shù)據(jù)整體映射到最方便表示這組數(shù)據(jù)的坐標(biāo)軸上，映射時沒有利用任何數(shù)據(jù)內(nèi)部的分類信息，是無監(jiān)督的，而LDA是由監(jiān)督的，增加了種屬之間的信息關(guān)系后，結(jié)合顯著性差異標(biāo)準(zhǔn)測試(克魯斯卡爾-沃利斯檢驗和兩兩Wilcoxon測試)和線性判別分析的方法進(jìn)行特征選擇。除了可以檢測重要特征，他還可以根據(jù)效應(yīng)值進(jìn)行功能特性排序，這些功能特性可以解釋頂部的大部分生物學(xué)差異。詳細(xì)說明可以參考這篇文章http:/blog.csdn.net/sunmenggmail/article/details/8071502 。不同顏色代表不同樣本或組之間的顯著差異物種。使用LefSe軟件分析獲得，其中顯著差異的logarithmic LDA score設(shè)為2。問題：LDA分析有什么用？回答：組間差異顯著物種又可以稱作生物標(biāo)記物（biomarkers），該分析主要是想找到組間在豐度上有顯著差異的物種。14. 物種進(jìn)化樹的樣本群落分布圖是將不同樣本的群落構(gòu)成及分布以物種分類樹的形式在一個環(huán)圖中展示。數(shù)據(jù)經(jīng)過分析后，將物種分類樹和分類豐度信息通過軟件GraPhlAn(http:/huttenhower.sph.harvard.edu/GraPhlAn )進(jìn)行繪制。其目的是將物種之間的進(jìn)化關(guān)系以及不同樣本的物種分布豐度和最高分布樣本的信息在一個視覺集中的環(huán)圖中一次展示，其提供的信息量較其他圖最為豐富。中間為物種進(jìn)化分類樹，不同顏色的分支代表不同的綱（具體的代表顏色見右上角的圖例），之后外圈的灰色標(biāo)示字母的環(huán)表示的是本次研究中比例最高的15個科（字母代表的科參見左上角的圖例）。之后的外圈提供的是熱力圖，如果樣本數(shù)<=10個則繪制樣本，如果樣本數(shù)超過10個則按照分組繪制，每一環(huán)為一個樣本，根據(jù)其豐度繪制的熱力圖。最外圈為柱狀圖，繪制的是該屬所占比例最高的樣本的豐度和樣本顏色（樣本顏色見環(huán)最下方的樣本名字的顏色）。其中熱力圖和柱狀圖取值均為原比例值x10000后進(jìn)行l(wèi)og2轉(zhuǎn)換后的值參考文獻(xiàn)：1. Vazquez-Baeza Y, Pirrung M, Gonzalez A, Knight R. 2013. Emperor: A tool for visualizing high-throughput microbial community data. Gigascience 2(1):16.2. Legendre, P. and Legendre, L. 1998. Numerical Ecology. Second English Edition. Developments in Environmental Modelling 20. Elsevier, Amsterdam.3. Segata N, Izard J, Waldron L, et al. Metagenomic biomarker discovery and explanationJ. Genome Biol, 2011, 12(6): R60.4. Langille MGI, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA et al. (2013). Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol 31: 814821.15. 物種相關(guān)性分析根據(jù)各個物種在各個樣品中的豐度以及變化情況，計算物種之間的相關(guān)性，包括正相關(guān)和負(fù)相關(guān)。相關(guān)性分析使用CCREPE算法，首先對原始16s測序數(shù)據(jù)的種屬數(shù)量進(jìn)行標(biāo)準(zhǔn)化，然后進(jìn)行Spearman和Pearson秩相關(guān)分析并進(jìn)行統(tǒng)計檢驗，計算出各個物種之間的相關(guān)性，之后在所有物種中根據(jù)simscore絕對值的大小，挑選出相關(guān)性最高的前100組數(shù)據(jù)，基于Cytoscap繪制共表達(dá)分析網(wǎng)絡(luò)圖，網(wǎng)絡(luò)圖采用兩種不同的形式表現(xiàn)出來。物種相關(guān)性網(wǎng)絡(luò)圖A：圖中每一個點代表一個物種，存在相關(guān)性的物種用連線連接，其中，紅色的連線代表負(fù)相關(guān)，綠色的先代表正相關(guān)，連線顏色的深淺代表相關(guān)性的高低。 物種相關(guān)性網(wǎng)絡(luò)圖B：圖中每一個點代表一個物種，點的大小表示與其他物種的關(guān)聯(lián)關(guān)系的多少，其中與之有相關(guān)性的物種數(shù)越多，點的半徑和字體越大，連線的粗細(xì)代表兩物種之間相關(guān)性的大小，連線越粗，相關(guān)性越高。參考文獻(xiàn)：Schwager E, Weingart G, Bielski C, et al. CCREPE: Compositionality Corrected by Permutation and RenormalizationJ. 2014.16. 聚類分析根據(jù)OUT數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理（1wlog10）之后，選取數(shù)目最多的前60個物種，基于R heatmap進(jìn)行作圖，熱圖中的每一個色塊代表一個樣品的一個屬的豐度，樣品橫向排列，屬縱向排列，兩個熱圖，差異是是否對樣品進(jìn)行聚類，從聚類中可以了解樣品之間的相似性以及屬水平上的群落構(gòu)成相似性。如果聚類結(jié)果中出現(xiàn)大面積的白或黑是因為大量的菌含量非常低，導(dǎo)致都沒有數(shù)值，可以在繪制之前進(jìn)行標(biāo)準(zhǔn)化操作，對每一類菌單獨自身進(jìn)行Z標(biāo)準(zhǔn)化。17. 群落功能差異分析通過對已有測序微生物基因組的基因功能的構(gòu)成進(jìn)行分析后，我們可以通過16s測序獲得的物種構(gòu)成推測樣本中的功能基因的構(gòu)成，從而分析不同樣本和分組之間在功能上的差異（PICRUSt Nature Biotechnology, 1-10. 8 2013）。通過對宏基因組測序數(shù)據(jù)功能分析和對應(yīng)16s預(yù)測功能分析結(jié)果的比較發(fā)現(xiàn)，此方法的準(zhǔn)確性在84%-95%，對腸道微生物菌群和土壤菌群的功能分析接近95%，能非常好的反映樣品中的功能基因構(gòu)成。為了能夠通過16s測序數(shù)據(jù)來準(zhǔn)確的預(yù)測出功能構(gòu)成，首先需要對原始16s測序數(shù)據(jù)的種屬數(shù)量進(jìn)行標(biāo)準(zhǔn)化，因為不同的種屬菌包含的16s拷貝數(shù)不相同。然后將16s的種屬構(gòu)成信息通過構(gòu)建好的已測序基因組的種屬功能基因構(gòu)成表映射獲得預(yù)測的功能結(jié)果。（根據(jù)屬這個水平，對不同樣本間的物種豐度進(jìn)行顯著性差異兩兩檢驗，我們這里的檢驗方法使用STAMP中的two-sample中T-TEST方法，Pvalue值過濾為0.05，作Extent error bar圖。）此處提供COG，KO基因預(yù)測以及KEGG代謝途徑預(yù)測。用戶也可自行使用我們提供的文件和軟件（STAMP）對不同層級以及不同分組之間進(jìn)行統(tǒng)計分析和制圖，以及選擇不同的統(tǒng)計方法和顯著性水平。參考文獻(xiàn)：Donovan H. Parks1 ,Gene W. Tyson,STAMP: statistical analysis of taxonomic and functional profiles, Bioinformatics(2014)30(21):3123-3124.doi:10.109318. COG構(gòu)成差異分析圖圖中不同顏色代表不同的分組，列出了COG構(gòu)成在組間存在顯著差異的功能分類以及在各組的比例，此外右側(cè)還給出了差異的比例和置信區(qū)間以及P-value。19. KEGG代謝途徑差異分析圖通過KEGG代謝途徑的預(yù)測差異分析，我們可以了解到不同分組的樣品之間在微生物群落的功能基因在代謝途徑上的差異，以及變化的高低。為我們了解群落樣本的環(huán)境適應(yīng)變化的代謝過程提供一種簡便快捷的方法。圖解讀：圖中不同顏色代表不同的分組，列出了在第三層級的構(gòu)成在組間存在顯著差異的KEGG代謝途徑第三層分類以及在各組的比例，此外右側(cè)還給出了差異的比例和置信區(qū)間以及P-value。本例圖所顯示的是第三層級的KEGG代謝途徑的差異分析，也可以針對第二或第一層的分級進(jìn)行分析。20. 基因的差異分析圖除了能對大的基因功能分類和代謝途徑進(jìn)行預(yù)測外，我們還能提供精細(xì)的功能基因的數(shù)量和構(gòu)成的預(yù)測，以及進(jìn)行樣本間以及組間的差異分析，并給出具有統(tǒng)計意義和置信區(qū)間的分析結(jié)果。這一分析將我們對于樣本群落的差異進(jìn)一步深入到了每一類基因的層面。圖解讀：圖中不同顏色代表不同的分組，列出了在組間/樣本間存在顯著差異的每一個功能基因（酶）以及在各組的比例，此外右側(cè)還給出了差異的比例和置信區(qū)間以及P-value。21. 在獲得標(biāo)準(zhǔn)報告后如果希望單獨修改分組或?qū)δ承┙M之間進(jìn)行顯著性差異分析，可以使用STAMP軟件在自己的電腦上進(jìn)行數(shù)據(jù)分析。STAMP提供了豐富的統(tǒng)計檢驗方法和圖形化結(jié)果的輸出。在使用STAMP之前需要首先準(zhǔn)備需要的spf格式文件和樣品分組信息表。在我們的報告中已經(jīng)將KEGG和KO以及COG的結(jié)果文件后經(jīng)過轉(zhuǎn)換生成了適用于STAMP軟件打開的spf格式文件，還有對應(yīng)的分組信息表文件groupfile.txt。以下是使用STAMP時的一些相關(guān)問題，詳細(xì)的STAMP使用教程可以參考我們提供的STAMP使用教程。1、stamp作圖用的原始數(shù)據(jù)的來源？STAMP 可以直接使用來自QIIME的biom文件和PICUST的KEGG和ko 文件，groupfile.txt文件的格式為tab-saperated value (tab鍵隔開的數(shù)據(jù))2、分組問題：導(dǎo)入數(shù)據(jù)之后，viewgroup legend ,在窗口右側(cè)會出現(xiàn)分組欄，根據(jù)需要進(jìn)行分組。3、Unclassiffied選項中，remain Unclassiffied reads、remove Unclassiffied reads、和use only for calculating frequency profiles 方法的區(qū)別？remain Unclassiffied reads和use only for calculating frequency profiles方法會保留所有的數(shù)據(jù)，而remove Unclassiffied reads僅僅保留有確定分組信息的數(shù)據(jù)。4、Statistical test 中，Welchs t-test、t-test、whites non-parametric t-test的區(qū)別，各自優(yōu)缺點？為了確保統(tǒng)計學(xué)意義和準(zhǔn)確度和精確性，需要足夠多的樣本數(shù)目，t-test檢驗可以在最少樣本數(shù)為4的時候確保高的準(zhǔn)確度和精確性。當(dāng)兩個樣本之間具有相同方差的時候，用t-test更為準(zhǔn)確，當(dāng)兩個樣本沒有相同方差，Welchs t-test更為準(zhǔn)確。當(dāng)樣本數(shù)目少于8的時候，可以使用whites non-parametric t-test，該計算時間較長，當(dāng)樣本數(shù)目過多的時候不宜使用該方法。5、Two-group 中type: one side 和 two side 的區(qū)別？One side 只會顯示前一個group與后一個group差異的比例，而two side 兩者之間的比例均會顯示。6、 STAMP在使用時首先打開了一個分析文件，如果新打開一個可能會導(dǎo)致顯示錯誤？目前版本的STAMP存在一些小問題，一次分析只能使用一個數(shù)據(jù)文件，如果要打開新的需要關(guān)閉軟件后再打開。22. 環(huán)境因子分析RDA分析CCA/RDA分析基于對應(yīng)分析發(fā)展的一種排序方法，將對應(yīng)分析與多元回歸分析相結(jié)合，每一步計算均與環(huán)境因子進(jìn)行回歸，又稱多元直接梯度分析。主要用來反映菌群與環(huán)境因子之間的關(guān)系。RDA 是基于線性模型，CCA是基于單峰模型。分析可以檢測環(huán)境因子、樣品、菌群三者之間的關(guān)系或者兩兩之間的關(guān)系。橫軸和縱軸：RDA 和CCA 分析，模型不同，橫縱坐標(biāo)上的刻度為每個樣品或者物種在與環(huán)境因子進(jìn)行回歸分析計算時產(chǎn)生的值，可以繪制于二維圖形中。圖解讀：冗余分析可以基于所有樣品的OTU作圖，也可以基于樣品中優(yōu)勢物種作圖；箭頭射線：箭頭分別代表不同的環(huán)境因子；夾角：環(huán)境因子之間的夾角為銳角時表示兩個環(huán)境因子之間呈正相關(guān)關(guān)系，鈍角時呈負(fù)相關(guān)關(guān)系。環(huán)境因子的射線越長，說明該影響因子的影響程度越大； 不同顏色的點表示不同組別的樣品或者同一組別不同時期的樣品，圖中的拉丁文代表物種名稱，可以將關(guān)注的優(yōu)勢物種也納入圖中； 環(huán)境因子數(shù)量要少于樣本數(shù)量，同時在分析時，需要提供環(huán)境因子的數(shù)據(jù)，比如 pH值，測定的溫度值等。23. 有其他問題可以聯(lián)系谷禾信息。END注意事項 假設(shè)您要對人體腸道微生物菌群進(jìn)行測序，一般可以選擇對糞便進(jìn)行取樣。一般建議采取新鮮樣品，因為微生物本身是活的群體，樣品長期保存或不在原環(huán)境下保存會改變原有菌群的構(gòu)成，最終導(dǎo)致我們得到的菌群構(gòu)成發(fā)生偏差。比如糞便樣品如果4度低溫保存了一段時間，則其中部分耐低溫的菌可能仍然在持續(xù)繁殖，這樣最終樣品的菌群會發(fā)生偏差。如果是送往公司測序也建議首先對樣品進(jìn)行DNA提取后再寄送，因為原始樣品的寄送過程也可能導(dǎo)致菌群變化。