微生物多樣性分析美吉生物講義

1、微生物多樣性標準分析微生物生信開發(fā)部朱娟中級生信分析工程師meta1234群落組成分析- Diversity排序分析聚類分析群落組成分析第一章群落組成分析概念群落 (Biocoenosis)或稱為“生基本原則一個群落環(huán)境的小環(huán)境越是多樣，那么其中的群落就越是物種豐富。（熱帶雨林）一個群落環(huán)境的非生物因素和全球平均值相差越物群落”。指的是學中，在一個群落生境內(nèi)相互之間具有直接或間接關(guān)系的所有生物，即生物的總合。大，其物種也會越少，但數(shù)目卻會越多。這個區(qū)域的群落通常是高特異性的。（鹽湖，深海)一個群落環(huán)境的生存環(huán)境越是緩慢連續(xù)的交替改變，群落也會越是豐富。（珊瑚礁)OTUs的數(shù)量+豐度=群落展示方

2、法群落otu2otu3otu1otu rep species群落組成分析方法一：柱狀圖輸出文件輸入文件R & barplotPhylum表格柱狀圖TaxonABCD.Actinobacteria12335214239.Bacteroidetes.Cyanobacteria1439230871077010328.Firmicutes297367554034525.Fusobacteria1021.Gemmatimonadetes1030.Parcubacteria91115.Proteobacteria23941624533864702.Saccharibacteria81412.Ten

3、ericutes0001.群落組成分析方法二：餅圖輸入文件輸出文件R & piePhylum表格餅圖TaxonABCD.Actinobacteria12335214239.Bacteroidetes.Cyanobacteria1439230871077010328.Firmicutes297367554034525.Fusobacteria1021.Gemmatimonadetes1030.Parcubacteria91115.Proteobacteria23941624533864702.Saccharibacteria81412.Tenericutes0001.群落組成分析方法三：

4、Venn輸出文件輸入文件R & plotVennVennOTU表格TaxonABCD.Actinobacteria12335214239.Bacteroidetes.Cyanobacteria1439230871077010328.Firmicutes297367554034525.Fusobacteria1021.Gemmatimonadetes1030.Parcubacteria91115.Proteobacteria23941624533864702.Saccharibacteria81412.Tenericutes0001.群落組成分析方法四：熱圖輸入文件輸出文件R &

5、 heatmapPhylum表格熱圖TaxonABCD.Actinobacteria12335214239.Bacteroidetes.Cyanobacteria1439230871077010328.Firmicutes297367554034525.Fusobacteria1021.Gemmatimonadetes1030.Parcubacteria91115.Proteobacteria23941624533864702.Saccharibacteria81412.Tenericutes0001.- Diversity第二章- Diversity概 念u 在地區(qū)尺度上，物種組成沿著某個梯

6、度方向從一個群落到另一個群落的變化率。即沿著某一環(huán)境梯度，物種替代的速率、物種周轉(zhuǎn)率、生物變化速率等。u 不同群落或某環(huán)境梯度上不同點之間的共有種越少，多樣性越大。它可以指示生境被物種的程度；多樣性的測定值可以用來比較不同地段的生境多樣性；多樣性與多樣性一起生物異質(zhì)性。了總體多樣性或一定地段的Ecological distance = Dissimilarity- DiversitySimilarity- 物種組成越相近，相似性越高- 指數(shù)范圍： 0 - 1 .Dissimilarity-dissimilarity = 1 - similarity- 通過特定的距離計算方法 (Euclidea

7、n).- Diversity Beta diversity metrics- pairwise sample dissimilarityu Non-phylogenetic(Component-based)“Distance”u Phylogenetic(UniFrac)S1S2S3S4S5S100.19840.2388830.2227630.259351S20.198400.1273240.1537680.160351S30.2388830.12732400.1943980.147496S40.2227630.1537680.19439800.20812S50.2593510.1603510

8、.1474960.208120- Diversity Non-phylogenetic metrics - Component-basedDistanceMethods:EuclideanManhattan Bray Curtis JaccardHellinger Pearson.TaxonS1S2S3Acidobacteria150Actinobacteria30801357642Bacteroidetes113Chlorobi010Chloroflexi51055Cyanobacteria0419Fibrobacteres001Firmicutes119224682732Fusobacte

9、ria031315OTU nameS1S2S3OTU1333150672OTU278015630OTU34168367OTU481367162OTU52549353OTU6471228OTU713239180- Diversity Phylogenetic metrics- Diversity Non-phylogenetic metrics距離名稱距離特征Bray_curtis同時考慮物種有無和物種豐度，不考慮各物種之間的進化關(guān)系或關(guān)聯(lián)信息。Bray_curtis 和bray_curtis_faith是標準化的曼哈頓距離；而bray_curtis_magurran是定量的sorensen距離

10、Bary_curtis_faithBray_curtis_magurranAbund_jaccardJaccard距離只考慮物種有無，而abund_jaccard添加了物種豐度Manhattan考慮物種豐度，但沒有考慮各物種之間的相關(guān)性。歐氏距離對（大的）異常值較為敏感，而曼哈頓距離卻對異常值不太敏感。若 樣本區(qū)分不明顯，選擇歐氏距離，相當于放大異常值，從而放大區(qū)分EuclideanCanberra考慮物種豐度，但沒有考慮各物種之間的相關(guān)性。受異常值的影響較小，使其適合應(yīng)用于具有高度偏斜或含異常值的數(shù)據(jù)chisq考慮物種豐度。利用列聯(lián)表分析的方法得到一個卡方統(tǒng)計量來衡量兩之間的差異性?？ǚ浇y(tǒng)計

11、量越大，表明的選擇對變量的取值有顯著影響，這也意味著兩之間的差異越大。稀有物種會對距離產(chǎn)生明顯影響，應(yīng)先剔除稀有物種后再進行分析gower考慮物種豐度，對數(shù)據(jù)進行了標準化，減小異常值的影響hellinger通過物種在樣本中的相對豐度進行距離計算，比絕對豐度更精確kulczynski考慮到兩樣本間物種豐度，及共有和獨有物種的豐度- Diversity Non-phylogenetic metrics距離名稱距離特征morisita_horn考慮物種相對豐度soergel考慮物種豐度specprof考慮物種相對豐度pearson相關(guān)系數(shù)是衡量隨 量X與Y相關(guān)程度的 法，相關(guān)系數(shù)的取值范圍是-1,1

12、。相關(guān)系數(shù)的絕對值越大，則表明X與Y相關(guān)度越高。當X與Y線性相關(guān)時，相關(guān)系數(shù)取值為1（正線性相關(guān)）或-1（負線性相關(guān)）?？紤]物種豐度，適合符合正態(tài)分布的樣本Spearman_approx利用兩變量的秩次大小作線性相關(guān)分析，具體是將兩要素的樣本值按數(shù)據(jù)的大小順序排列位次，以各要素樣本值的位次代替實際數(shù)據(jù)而求得的一種統(tǒng)計量。Spearman對原始變量的分布不作要求，屬于非參數(shù)統(tǒng)計方法，適用范圍要廣些Binary_距離名稱算法名稱前帶有“binary-”的算法為先將OTU表中的數(shù)值轉(zhuǎn)換為二進制布爾類型，再進行計算。例如“binary_euclidean”算法，先將OTU表中的數(shù)值為零的不變，大于零的

13、變?yōu)?，再進行“euclidean”距離分析。所有距離只要換算成二進制算法，均不再考慮物種豐度- Diversity Phylogenetic metrics距離名稱距離特征Unweighted_unifrac算法名稱中含有“unifrac”的算法需要各個物種分類單元（如OTU、屬等）的系統(tǒng)進化樹（程序會自動搜索項目中存在的完整進化樹），通過計算進化樹各物種的系統(tǒng)發(fā)育進化關(guān)系，從而計算樣本間距離，其中unweighted UniFrac距離算法沒有計入不同環(huán)境樣本的序列相對豐度，而weighted UniFrac算法在計算樹枝長度時將序列的豐度信息進行 計算，因此unweighted UniF

14、rac可以檢測樣本間變化的存在，而weighted UniFrac可以更進一步定量的檢測樣本間不同譜系上發(fā)生的變異。如果未添加full_tree信息，默認選擇豐度前50的OTU進行構(gòu)建進化樹，而添加full_tree則表明利用所有OTU構(gòu)建進化樹。Unweighted_unifrac_full_treeWeighted_normalized_unifracWeighted_unifrac排序分析第三章排序分析應(yīng)用群落學數(shù)據(jù)一般是數(shù)據(jù)，例如物種屬性或環(huán)境因子。多元統(tǒng)計分析是群落學常用的分析方法，而排序（ordination）是多元統(tǒng)計最常用的方法。多樣性分析中常用的排序分析方法：PCA、PCoA

15、、NMDS、RDA、db-RDA必須強調(diào)的是：排序?qū)ο?樣方，物種和環(huán)境因子)在排序空間的坐標值（ scores）一般來講，只有相對意義，沒有絕對意義。排序分析何為排序何為排序？排序分析概念u排序的過程是將樣方或植物種排列在一定的空間，使得排序軸能夠反映一定的梯度，從而，能夠解釋植被或植物種的分布與環(huán)境因子間的關(guān)系，也就是說排序是為了揭示植被-環(huán)境間的關(guān)系。u排序也叫梯度分析（gradient analysis）圖： 美國大煙山植被梯度分析一例排序分析分類u 間接梯度排序(indirect gradient analysis)：非約束性排序；尋求潛在的或在間接的環(huán)境梯度來解釋物種數(shù)據(jù)的變化；

16、如PCA、PCoA；u 直接梯度排序(direct gradient analysis)：約束性排序；在特定的梯度上（環(huán)境軸）上探討物種的變化情況； 如：RDA、db-RDA；排序分析物種響應(yīng)環(huán)境梯度模型：u 線性模型：包括直線和曲線線性關(guān)系某個植物種隨著某一環(huán)境因子的變化而呈線性變化或叫線性反應(yīng)（linear response）u 非線性模型：二次曲線模型：高斯模型（Gaussian m），正態(tài)曲線某個植物種的數(shù)隨某個環(huán)境因子值的增加而增加，當環(huán)境因子增加到某一值時，植物種的數(shù)達到最大值，此時的環(huán)境因子值稱為該種的最適值（optimum）；隨后當環(huán)境因子值繼續(xù)增加時，種的演變：單峰模型數(shù)逐漸

17、下降，最后消失排序分析物種響應(yīng)環(huán)境梯度模型：注：左邊圖表示兩個種對環(huán)境梯度反應(yīng)的不同模型， 右邊圖表示模型所對應(yīng)的關(guān)系，ad為線形模型，ef為單峰模型圖兩個種x、y在某一環(huán)境梯度上的關(guān)系類型排序分析-PCA 主成分分析Ø PCA：Principal components analysis是一種簡化數(shù)據(jù)集的技術(shù)，通過線性變化把數(shù)據(jù)變換到一個新的坐標系中， 使得任何數(shù)據(jù)投影的第一大方差在第一個坐標(稱為第一主成分)上，第二大方差在第二個坐標(第二主成分)上，依次類推。原理和方法：原理是設(shè)法將原來變量重新組一組新的相互無關(guān)的幾個綜合變量，同時根據(jù)實際需要從中可以取出幾個較少的總和變量盡可能

18、多地反映原來變量的信息；方法主要是通過對協(xié)方差矩陣進行特征分解，以得出數(shù)據(jù)的主成分（即特征向量）與它們的權(quán)值（即特征值）。排序分析-PCAM=0,那么F=1； M=1,那么F=0M或F去掉實際上沒有任何信息的損失，因為只要保留一列就可以完全還原另一列。IDMF101210310401510601710排序分析-PCA每條可以被表示為一個10維向量：(22347, 5867, 15051, 183, 95, 158, 1.62%, 28707, 2925583, 421)T降維當然意味著信息的丟失，不過鑒于實際數(shù)據(jù)本身常常存在的相關(guān)性，我 們可以想辦法在降維的同時將信息的損失盡量降低?！盀g覽量”

19、和“訪客數(shù)”往往具有較強的相關(guān)關(guān)系 “訂單數(shù)”和“成交數(shù)”往往也具有較強的相關(guān)關(guān)系時段流量時段付款時段上架時段瀏覽量(PV)訪客數(shù)(UV)銷售額銷售量訂單數(shù)單價轉(zhuǎn)化率成交量總銷售額高質(zhì)寶貝數(shù)1日22347586715051183951581.62%2870729255834212日111542638508156331541.25%1092511320451863日51031230234532141681.14%5015526410764日3194737945951890.68%1602181455445日2252555768861281.08%56268447166日51600.95%6277

20、4155137日30369041158971650.77%1013129333168日65642177385342321201.47%5003508392519日24257785810239105621650.79%91581137090211排序分析-PCA計算過程第一步，分別求每組的平均值，然后對于所有的樣例，都減去對應(yīng)的均值。第二步，求特征協(xié)方差矩陣第三步，求協(xié)方差的特征值和特征向量第四步，將特征值按照從大到小的順序排序，選擇其中最大的k個，然后將其對應(yīng)的k個特征向量分別作為列向量組成特征向量矩陣第五步，將樣本點投影到選取的特征向量上 PCA (用歐氏距離方程計算點間距離)PCA本質(zhì)上是

21、將方差最大的方向作為主要特征，并且在各個正交方向上將數(shù)據(jù)“離相關(guān)”，也就是讓它們在不同正交方向上沒有相關(guān)性。PCA也存在一些限制，例如它可以很好的解除線性相關(guān)，但是對于高階相關(guān)性就沒有辦法了。PCA假設(shè)數(shù)據(jù)各主特征是分布在正交方向上，如果在非正交方向上存在幾個方差較大的方向，PCA的效果就大打折扣了。排序分析-PCA實現(xiàn)方法stats package: prcomp參考文獻：Legendre P, Legendre L. Numerical Ecology. 2nd ed. Amsterdam: Elsevier, 1998. ISBN 978-0444892508.Jackson DA (1

22、993) StoppingRules in Principal Components Analysis: A Comparison of Heuristical and Statistical Approaches. Ecology. 74(8): 2204-2214. Ramette A (2007) Multivariate analyses in microbial ecology. FEMS Microbiol Ecol. 62(2): 142160.ter Braak CJF (1994) Canonical commuordination. Part I: basic theory

23、 and linear methods. Ecoscience. 1: 127140.排序分析-PCoA主坐標分析u PCoA (Principal coordinates analysis)：主坐標分析u 根據(jù)樣方之間的差異即距離的數(shù)據(jù)進行分析u N個樣本，在N-1維坐標系中，點間的歐氏距離的平方正好等于原來的差異數(shù)據(jù)排序分析-PCoA計算過程1)計算樣方間的距離系數(shù)，N×N 距離矩陣D協(xié)方差矩陣S(即離差矩陣)2)3)求S矩陣的特征根和相對應(yīng)的特征向量4)求排序坐標排序分析-PCoA實現(xiàn)方法stats package: cmdscale參考文獻：Gower JC (1966) S

24、ome distance properties of latent root and vector methods used in multivariate analysis. Biometrika. 53(3-4):325-338.Legendre P, Legendre L. Numerical Ecology. 2nd ed. Amsterdam: Elsevier, 1998. ISBN 978-0444892508.排序分析Differences between PCA and PCoA:u PCA是先給出N個點的坐標，去找出剛性旋轉(zhuǎn)的坐標，PCoA是只知其間的距離去重新建立各自的坐

25、標。u PCA is used for similarities and PCoA for dissimilaritties。u PCA僅僅比較的是OTU豐度的不同，PCoA分析可以加入進化距 離的內(nèi)容。排序分析-NMDS非度量尺度分析尺度法(NMDS)：非度量是一種將空間的研究對象（樣本或變量）簡化到低進行、分析和歸類，同時又保留對象間原始關(guān)系的數(shù)據(jù)分析方法。其基本特征是將對象間的相似性或相異性數(shù)據(jù)看成點間距離的單調(diào)函數(shù)，在保持原始數(shù)據(jù)次序關(guān)系的基礎(chǔ)上，用新的相同次序的數(shù)據(jù)列替換原始數(shù)據(jù)進行度量型尺度分析。Ø 間接梯度排序Ø 接受任何距離矩陣Ø 秩次排序方法&#

26、216; 計算量較大，對硬件要求高排序分析-NMDS計算過程Kruskal 算法 在給定的 Rk 空間里選擇一個擬合構(gòu)圖, 其 Stress(k) 值越小越好:1. 計算一種距離得到距離矩陣，確定需要降到的緯度，一般是2維；2. 隨機選擇一個2維的初始數(shù)據(jù)；3. 在這個二stress值；尋找一組點，滿足距離與原始距離的順序相同,并計算4. 尋找一個新的維度,使得stress減??；5. 重復 3-4 直至 Stress(k=2) 收斂（即Stress不再減?。?，即為最終結(jié)果。如果預先設(shè)定的排序軸數(shù)量比較少（例如K=2或3），在相同軸數(shù)的條件下，NMDS往往能夠獲得比PCoA更少失真的對象之間的關(guān)

27、系。但是NMDS計算需要不斷迭代，對計算機要求比較高。排序分析-NMDSstress:依據(jù)Kruskal應(yīng)力公式I計算,是排序空間內(nèi)對象結(jié)構(gòu)與原始距離矩陣之間的相異程度的度量。lstress<0.2時，可以用NMDS的二維點圖表示,其圖具有一定的解釋意義stress<0.1時，可以認為是一個好的排序stress<0.05時，具有很好的代表性llstress值擬合度0.2以上非常不好0.2不好0.1普通0.05很好0.025非常好0完美排序分析-NMDS實現(xiàn)方法Ø vegan packagemetaMDS: performs multiple NMDS runs an

28、d retains the best solution stressplot: used to visualise the Shepard stress plotexample_NMDS=metaMDS(commu_matrix, # Our commu-by-species matrixk=2) # The number of reduced dimensionsexample_NMDS$stressØisoMDS()：所使用的距離矩陣帶有isoMDS()函數(shù)運行NMDS分析值，這個時候可以用Ølabdsv packagebestnmds:減少僅達到局部最小化應(yīng)力函數(shù)風險

29、參考文獻：Legendre P, Legendre L. Numerical Ecology. 2nd ed. Amsterdam: Elsevier, 1998. ISBN 978-0444892508. Legendre P, Borcard D, Peres-Neto PR (2005) Analyzing beta diversity: Partitioning the spatial variation of commucomposition data. Ecol Monogr. 75:435450.Clarke KR (1993) Non-parametric multivaria

30、te analysis of changes in commustructure. Austral J Ecol 18:117-143.Kenkel NC, Orloci L (1986) Applying Metric and Nonmetric Multidimensional Scaling to Ecological Studies: Some New Results. Ecology 67:919928.排序分析-RDA冗余分析llllRedundancy analysis, RDA約束化的主成份分析直接梯度分析數(shù)據(jù)矩陣和環(huán)境矩陣Actually:RDA => constrai

31、ned PCARDA或CCA選擇原則：先進行DCA分析，看分析結(jié)果中Lengths of gradient 的第一軸的大小，如果大于4.0，就應(yīng)該選CCA，如果3.0-4.0之間，選RDA和CCA均可， 如果小于3.0, RDA的結(jié)果要好于CCA。排序分析-RDA計算過程可以看作是PCA的每一步加上與環(huán)境因子線性回歸。PCA的軸為物種的線性組成；RDA的軸為環(huán)境因子的線性組成，是實際的環(huán)境梯度，約束的軸不超過環(huán)境因子個數(shù)；在尋找最大解釋方差的軸時，樣方的坐標與環(huán)境因子做回歸：q+ åbkUkjk =1z j= b0計算樣方排序新值zj(j= 1,2, N)Z=Ub排序分析-RDARD

32、A圖解1、點與點之間距離越大，分布差異越大；2、環(huán)境因子、物種、因子與物種間夾角為銳角則成正相關(guān)， 為鈍角則成負相關(guān)；3、樣品在物種上投影點的相對位置代表的是該物種這些樣方中擬合的多度值(即能夠被排序模型所解釋的部分，非觀測值） 的排序情況；4. 從樣方的點向數(shù)量環(huán)境因子的箭頭做投影，投影點的位置可以近似表示該環(huán)境因子數(shù)值在這些樣方內(nèi)的排序；5. 環(huán)境因子箭頭的長短可以代表環(huán)境因子對于物種數(shù)據(jù)的影響程度（解釋量）的大小；排序分析-RDA實現(xiàn)方法vegan package: rda參考文獻：Legendre P, Legendre L. Numerical Ecology. 2nd ed. Am

33、sterdam: Elsevier, 1998. ISBN 978-0444892508. (2007) Multivariate analyses in microbial ecology. FEMS Microbiol Ecol. 62(2): 142160.Ramette Ater Braak CJF (1994) Canonical commuordination. Part I: basic theory and linear methods. Ecoscience. 1: 127140.排序分析-dbRDAdbRDADistance-based redundancy analysi

34、s,db-RDA可以看作是PCoA的每一步加上與環(huán)境因子線性回歸；步驟：1.Calculate your distance matrix 2.Run a PCoA3.Run an RDA on the eigenvaluesobtained from the PCoA排序分析-dbRDA實現(xiàn)方法vegan package:capscale參考文獻：Legendre P, Anderson MJ (1999) Distance- based redundancy analysis: testing experiments. Ecol Monogr.69:124.multispecies resp

35、onses in multifactorial ecological排序分析 總 結(jié) Ø PCA:Ø PCoA:Ø NMDS:Ø RDA:Principal components analysis Principal coordinates analysisNon-metric multidimensional scalingRedundancy analysisØ db-RDA: Distance based RDARaw data-baseDistance-baselinearunimodalunconstrainedPCACA,DCAPC

36、oA,NMDSconstrainedRDACCAdb-RDA聚類分析第四章聚類分析概念u 聚類是指按照或樣品(individuals, objects or subjects)的特征將它們分類，使同一類別內(nèi)的具有盡可能高的同質(zhì)性(homogeneity)，而類別之間則應(yīng)具有盡可能高的異質(zhì)性(heterogeneity)。u 為什么聚類？u 從數(shù)據(jù)的海洋中找到一種模式u 用于理解和分析數(shù)據(jù)聚類分析層次聚類層次聚類(Hierarchical cluster analysis):通過生成一系列嵌套的聚類樹來完成聚類。單點聚類處在樹的最底層，在樹的頂層有一個根節(jié)點聚類。根節(jié)點聚類覆蓋了全部的所有數(shù)據(jù)點。自下而上自上而下聚類分析層次聚類(Hierarch