生物信息學習的正確姿勢

NGS系列文章包括NGS基礎、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。

222Beta多樣性限制性排序CPCoA/CCA/RDA/LDA

本節作者：文濤，南京農業大學；劉永鑫，中科院遺傳發育所

版本1.0.5，更新日期：2020年8月12日

本項目永久地址：https://github.com/YongxinLiu/MicrobiomeStatPlot ，本節目錄 222CPCoA，包含R markdown(*.Rmd)、Word(*.docx)文檔、測試數據和結果圖表，歡迎廣大同行幫忙審核校對、并提修改意見。提交反饋的三種方式：1. 公眾號文章下方留言；2. 下載Word文檔使用審閱模式修改和批注后，發送至微信(meta-genomics)或郵件(metagenome@126.com)；3. 在Github中的Rmd文檔直接修改并提交Issue。審稿人請在創作者登記表 https://www.kdocs.cn/l/c7CGfv9Xc 中記錄個人信息、時間和貢獻，以免專著發表時遺漏。

基本概念

在本章節中以下概念可以通用：

限制性=約束性=有監督=典范

非限制性=非約束性=無監督

上一節我們主要講了非限制性排序分析(Unconstrained ordination / indirect gradient analyis)《221.Beta多樣性PCoA和NMDS排序》。今天重點講限制性(constrained ordination / direct gradient analysis)排序分析。

限制性排序，又稱為限制性排序，不同于非約束排序。因為限制性排序從一開始就考慮到了解釋變量，并進行運算，然后提取出來和解釋變量相關的數據結構。其實非限制性排序也可以加入解釋變量，只是在計算之后進行的被動添加，也就是說非限制性排序并不受外部變量的影響。

從表達的意義上來講限制性排序可以探索解釋變量矩陣和響應變量矩陣之間的關系，而非限制性排序只是一種描述響應變量矩陣的，探索性的分析方法。

限制性排序是兩個矩陣之間關系的尋找，這就牽扯到不同的限制性排序方法。如果僅僅關注于微生物生態學領域，經常用的有兩種，分別是冗余分析（RDA）和典范對應分析（CCA）。這兩種方法作為典型的生態方向的應用都是使用多元回歸尋找解釋變量矩陣和響應變量矩陣之間的關系，分別結合PCA排序或者CA排序展示結果。兩個矩陣之間的關系往往可以做顯著性分析，這里經常使用置換檢驗。

這里額外提及LDA和CPCoA，他們也是限制性排序，只是將連續型變量替換為分類型變量(分組)，目的是尋找能夠解釋已知分組對象的定量解釋變量的組合，也就是分組可以解釋響應變量矩陣的程度。

總結一下常見排序分析，包括主成分分析（Principal Component Analysis，PCA）、主坐標分析（Principal Coordinate Analysis，PCoA）、對應分析（Correspondence Analysis，CA）、去趨勢對應分析（Detrended Correspondence Analysis，DCA）和非度量多維尺度分析（Non-metric Multi-Dimensional Scaling，NMDS)、冗余分析（Redundancy Analysis，RDA）、典范對應分析（Canonical Correspondence Analysis，CCA）都屬于降維排序分析方法，其區別總結成一個表格如下：

表1. 常見限制性和約束制排序分類

類型原始數據(線性模型)原始數據(單峰模型)距離矩陣(多種距離)

非限制性排序(間接梯度排序)	PCA	CA，DCA	PCoA，NMDS
限制性排序(直接梯度排序)	RDA	CCA	CPCoA，db-RDA

冗余分析(RDA）

主成分分析（Principal Component Analysis，PCA）結合多元回歸分析，即為冗余分析（Redundancy analysis，RDA）。作為限制性排序的代表方法之一，在微生物生態學領域經常用于尋找環境變量和微生物群落之間的關系等。首先計算解釋變量同微生物群落之間的多元多重線性回歸，其次使用回歸擬合值矩陣進行PCA排序。如果回歸擬合值矩陣使用對應分析（Correspondence Analysis，CA）排序擬合，那么就是典范對應分析（CCA）。

如果有解釋變量作為協變量，從此衍生出來偏RDA分析，表示在控制協變量后解釋變量可以解釋的響應變量矩陣。

從上一節《221.Beta多樣性PCoA和NMDS排序》我們也逐漸認識到，隨著生態學的發展，以歐氏距離為代表的傳統分析方法并不能很好解決微生物組數據，由此衍生出來的Bray-Curtis距離等諸多適合于生態學領域的距離。這些距離是沒有辦法使用RDA分析的，隨之而來就有了基于距離的冗余分析（Distance-based redundancy analysis，db-RDA）。當然不僅僅是生態學的距離，任何距離都可以使用db-RDA。隨著研究進一步發展，我們了解到微生物和解釋變量之間會有大量的非線性關系，比如單峰響應，此時可以使用非線性關系的CCA分析。

無論是哪種方法，發展到如今，都可以通過R包Vegan進行很好的處理生態學數據，下面我們就一個公開發表文章中的案例進行結果解讀，而后進行簡單的示范性演示。

RDA分析結果得到排序，總結了響應矩陣中主要的變化模式，可以通過解釋變量矩陣來解釋。可以用于討論選擇適當的可視化尺度并解釋此排序。在RDA的結果中經常會顯示有關多個約束軸（RDA軸）和非約束軸（PCA軸）的信息。每個RDA軸都有一個與之關聯的特征值。由于總方差等于所有特征值的總和（受約束的），因此每個軸解釋的方差比例只是給定特征值與解的總方差的商。有時，殘差的排序和/或殘差之間的相關性可能比特征明確的殘差更具有生態學意義。通過排序和相關性檢查RDA解決方案的非規范（不受約束）向量，可以部分解釋殘差。或者，可以在對一組響應變量執行多元線性回歸（multiple linear regression，MLR）之后，對殘差矩陣執行PCA分析，探索生態學意義。RDA可以將PCA軸與RDA軸一起顯示。PCA軸匯總了不受約束的（殘差）方差。

閱讀RDA排序的雙序圖和三序圖

RDA結果可以表示為雙序圖或三序圖（圖1）。這些圖的解釋取決于所選擇的標尺。通常，如果對象之間的距離具有特定值，或者大多數解釋性變量是邏輯值，則考慮I型標準化(type I scaling，又常被翻譯為歸一化、縮放)。如果對變量之間的相關關系更感興趣，請考慮使用II型標準化。進一步的解釋在下面討論。

圖1. RDA雙序圖(a)和RDA三序圖(b)的示意圖。

a）RDA雙序圖將對象作為點，將響應或解釋變量作為向量（紅色箭頭）。標稱變量(Levels of nominal variables)的級別繪制為點（紅色）。

b）在一個三序圖中，對象被標為點（藍色），而響應變量和解釋變量（分別為紅色和綠色箭頭）被繪制為矢量。物種繪制為點（綠色）。

• I型標準化-距離圖（以對象為中心）

目標點之間的距離近似于歐幾里得距離。因此，可以期望將靠近在一起的對象具有相似的變量值。這并非總是如此，因為RDA僅展示了數據集中的部分變化（圖2）。

代表響應變量的矢量之間的角度是沒有意義的。

代表響應變量的點和代表解釋變量的矢量之間的角度反映了它們的（線性）相關性。

• II型標準化-相關圖（關注響應變量）

對象點之間的距離不是近似歐幾里得距離。

對象在表示響應變量的點上的直角投影近似于給定對象的變量值。

所有矢量之間的角度反映了它們的（線性）相關性。相關性等于矢量之間的角度的余弦值。

圖2. RDA圖坐標對象的投影和角度示例。
a）坐標對象在矢量上的投影，b）矢量之間的角度。

坐標點在變量向量上的投影（如圖2a中的點i所示）近似于為該對象實現的變量值。因此，目測與大多數其他對象相比，可以預期對象i受變量1影響較大。但是，可以預期對象ii相對于其他對象受的變量1影響較小。請注意，虛線通常不會在雙序圖中顯示，此處為演示而顯示。使用II型尺度縮放時，向量之間的夾角余弦近似于它們表示的變量之間的相關性。在這種情況下，∠a接近90，這表明變量“ 1”和“ 2”顯示出很小的相關性（即它們幾乎是正交的，就像獨立的軸一樣）。∠b小于90°，表明變量“ 2”和“ 3”之間呈正相關，而∠c接近180°，表明變量“ 2”和“ 4”之間具有強烈的負相關（即變量“增加”的方向） 2”和“ 4”彼此相對）。變量5是非數值解釋變量，由點表示。在變量4上的直角投影表明兩者是正相關的。

典范對應分析(CCA）

同RDA有許多相似之處，作為限制性排序典型算法的CCA，是基于單峰模型，同CA分析一樣的基礎算法，但在排序迭代過程中將解釋變量同響應變量矩陣之間做了加權回歸分析。這里重點關注CCA使用的條件：1. 物種在整個生態梯度環境中都有采樣；2. 物種在生態梯度上呈現單峰分布響應。

限制性主坐標分析(CPCoA）

限制性主坐標分析(Constrained PCoA，CPCoA），也稱限制性主坐標軸分析。是在PCoA的的分析中即添加了分組信息，以尋找在自定分組條件下可最大解釋組間差異的平面。結果展示為在指定分組條件下的解析率，要求至少3組及以上才可開展此類分析。

library(amplicon) (p1=beta_pcoa(beta_bray_curtis, metadata, groupID="Group")) ggsave("b5.pcoa.pdf", p1, width=91, height=59, units="mm") ggsave("b5.pcoa.png", p1, width=91, height=59, units="mm") (p2=beta_cpcoa_dis(beta_bray_curtis, metadata, groupID="Group")) ggsave("b6.cpcoa.pdf", p2, width=91, height=59, units="mm") ggsave("b6.cpcoa.png", p2, width=91, height=59, units="mm") library(cowplot) (p0 = plot_grid(p1, p2, labels = c("A", "B"), ncol = 2)) ggsave("b3_pcoa_cpcoa.pdf", p0, width=91*2, height=59, units="mm") ggsave("b3_pcoa_cpcoa.png", p0, width=91*2, height=59, units="mm")

圖3. PCoA和CPCoA比較。

表2. PCoA和CPCoA差異比較

特征PCoACPCoA

解析率	整體是100%	整體僅為0%~100%間，見左上角或圖注
坐標軸	PCo1/2，解析率為整體的百分比	CPCo1/2，解析率為本圖的百分比
組間分開	不明顯	明顯
顯著性	需要單獨統計	整體anova.cca統計值見左上角
分析數據量	>= 3個樣本	>= 3個實驗組

線性判別分析(LDA）

線性判別分析(Linear Discriminant Analysis,LDA)在模式識別領域（比如人臉識別，艦艇識別等圖形圖像識別領域）中有非常廣泛的應用，在生物學大數據研究中同樣也有廣泛應用，比如Sicence封面文章哈扎人菌群研究就使了此方法(詳見例4)。

LDA是一種監督學習的降維技術，也就是說它的數據集的每個樣本是有類別輸出的。這點和PCA不同。PCA是不考慮樣本類別輸出的無監督降維技術。LDA的思想可以用一句話概括，就是“投影后類內方差最小，類間方差最大”。更多信息詳見下文。

• LDA分析、作圖及添加置信-ggord

實例解讀

例1. CPCoA組間整體差異

本文是2016年德國馬普Paul Schulze-Lefert組發表于PNAS的文章 ，分析了百脈根根瘤的微生物群落組成，同時在根瘤缺失突變體中發現根和根際微生物組均有較大差異的變化。

圖2. 散點圖展示限制性主坐標軸分析取材部位和基因型間的差異。

(A) 基于野生型樣本的Bray-Curtis距離采用限制性主坐標分析不同取樣部分樣本間的差異。該條件下基于94個樣本最大可解析整體19.97%的變異，且組間至少存在一個兩組間的顯著差異 (P<0.001)。

(B) 以基因型為條件分析最大解釋基因型組間差異的空間平面，可解釋9.82%的變異，并且基因型間存在顯著差異，其中作者按形狀標出了各基因型；同時作者還按取樣部分進行著色，在這一平面上取樣部仍能很好的分開。

Fig. 2. (A) Constrained PCoA plot of Bray–Curtis distances between samples including only the WT constrained by compartment (19.97% of variance, P > 0.001; n=94). (B) Constrained PCoA plot of Bray–Curtis distances constrained by genotype (9.82% of variance explained, P < 0.001; n=164). Each point corresponds to a different sample colored by compartment, and each host genotype is represented by a different shape. The percentage of variation indicated in each axis corresponds to the fraction of the total variance explained by the projection.

圖表結果：主要結果包括取樣位置可解釋19.97%差異，且區分明顯；突變體與野生型可以區分，但區分不大(占9.82%變異中的17.75%的縱軸上可區分)；各突變體間很難區分，完全混在一起；在基因型最大解釋平面上，取樣位置仍能非常好的在第一軸上區分。

圖2在原文中共被引用了11次。讓我們跟著作者學習兩個小圖如何寫成幾大段文字的。

值得注意的是，我們發現與四個共生突變體中的每一個的根和根際相關的群落彼此相似，但與野生型植物的群落顯著不同（圖2B）。采用典范主坐標分析整個群落，宿主基因型對細菌群落存在顯著影響，解釋了9.82％的變異（圖2B）。根瘤是根上的錨定結構，但兩個生境仍具有獨特的細菌組合（圖2A）。

Remarkably, we found that communities associated with the roots and rhizospheres of each of the four symbiosis mutants were similar to each other, but significantly different from the communities of WT plants (Fig. 2B). Furthermore, CAP performed on the entire dataset revealed a prominent effect of the host genotype on bacterial communities, explaining 9.82% of the variance (Fig. 2B). Nodules are root-derived and -anchored structures, and yet the two organs host distinctive bacterial assemblages (Fig. 2A).

由于測試的共生突變體之間的菌群沒有顯著差異（圖2B，詳者注：沒有顯著差異也是發現），我們使用來自所有土壤批次的所有突變體基因型的組合樣本進行了分析。接下來，我們直接比較了野生型和突變體，以確定根或根際差異豐富的OTU（附錄，圖S10和S11，詳者注：除主圖外還需要有大量的附圖和表對同一結論的支撐）。我們發現，宿主基因型的使用群落結構轉變（圖2B）很大程度上是由大量OTU引起的，這些OTU為突變體根樣品相對于野生型根中特異富集（n=45）或消減（n=15）至少屬于14個細菌綱（附錄，圖2B，圖S10A，S11 A和B）。

Due to the fact that the bacterial assemblages of the tested symbiosis mutants do not significantly differ among each other (Fig. 2B), we performed our analyses using the combined samples from all mutant genotypes across all soil batches. Next, we compared directly the WT and mutants to identify OTUs differentially abundant in the root or rhizosphere (SI Appendix, Figs. S10 and S11). We found that the community shift that separates host genotypes (Fig. 2B) is largely caused by numerous OTUs that are specifically enriched (n=45) or depleted (n=15) in WT roots with respect to mutant root samples (SI Appendix, Fig. S10A) belonging to at least 14 bacterial orders (SI Appendix, Fig. S11 A and B).

我們的統計分析揭示了宿主基因型對百脈根的根系微生物群的重大影響，這解釋了數據中9.82％的差異（圖2B）。這種基因型效應幾乎比野生型和嚴重免疫受損的擬南芥突變體的根相關群落（討論：與2015的Science結果比較）或擬南芥與三種3千5百萬年分化的十字花科物種之間鑒定的效應大兩倍（討論：與2014的PNAS結果比較）。有趣的是，在nfr5，nin和lhk1突變體根中檢測到的群落轉變具有高度類似（圖2B）。

Our statistical analyses revealed a major impact of the host genotype on the Lotus root microbiota, which explains 9.82% of the variance in the data (Fig. 2B). This genotype effect is almost twofold larger than the effect identified when comparing the root associated communities of WT and severely immunocompromised A. thaliana mutants (52) or between A. thaliana and three Brassicaceae species that diverged up to 35 Mya (29). Interestingly, the community shifts detected in nfr5, nin, and lhk1 mutant roots are highly comparable (Fig. 2B).

為了評估不同取樣部分對細菌群落組裝的影響，我們使用Bray–Curtis距離比較了β-多樣性（樣本間多樣性），并進行了典范主坐標分析。該分析揭示了根、根際、根瘤和土壤不同取樣部分的明顯區別，這解釋了高達數據總方差的19.97％（圖2A；P < 0.001），而歸因于 土壤批次相對較小（方差的8.01％，P < 0.001）。

To assess the effect of the different compartments on the assembly of bacterial communities, we compared the β-diversity (between- samples diversity) using Bray–Curtis distances and performed a canonical analysis of principal coordinates (CAP) (30) (Materials and Methods). This analysis revealed a clear differentiation of samples belonging to the root, rhizosphere, nodule, and soil compartments that explains as much as 19.97% of the overall variance of the data (Fig. 2A; P < 0.001), whereas the effect attributable to the soil batch was comparatively small (8.01% of the variance, P < 0.001).

我們利用了科隆土壤中存在的百脈根和根瘤菌之間的相容共生關系，并對在這種土壤中生長的野生型植物的附生植物枯竭的功能性根瘤細菌群落進行了分析（材料和方法）。我們發現，與根部和根際中的根瘤相比，根瘤中有一個獨特的細菌群落（圖2A）。?重要的是，根瘤和根的群落與根際相似（從第二主坐標軸分離；圖2A）。

We took advantage of the compatible symbiotic association between Lotus and rhizobia present in Cologne soil and performed an analysis of the bacterial community of epiphyte-depleted, functional nodules of WT plants grown in this soil (Materials and Methods). We found that nodules were inhabited by a distinctive bacterial community compared with those present in the root and rhizosphere (Fig. 2A). Importantly, the nodule and root communities were similarly divergent from the rhizosphere (separation by second component; Fig. 2A).

圖表結論或規律：取樣部位對微生物組成影響較大，基因型其次；不同根瘤突變體差異極小。

圖片優點：配色選擇各組區分較好，不同圖配色方案一致；圖片使用矢量圖線條和文字清楚(大部分PNAS和Microbiome文章中是位圖，經過PDF的壓縮，文字非常模糊)。個人建議，只要不是照片，畫的圖都用矢量，無極縮放不失真，一般體積還小，而且方便編輯修改。

例2. CAP細菌/真菌群落和雙標圖biplot

這篇文章由瑞士伯爾尼大學的Matthias Erb和Klaus Schlaeppi團隊合作于2018年發表在Nature Communications期刊發表關于根系分泌代謝物通過塑造根際微生物群來驅動植物土壤對生長和防御的反饋的研究，揭示了植物根際微生物區系組成、植物性能和植物-草食相互作用的機制。詳細中文解讀

圖7. 苯并惡嗪類化合物(BX)調節可在第二代植物的根和根際中塑造土壤微生物群。

a，b對先前曾受野生型（BX +）或苯并惡嗪類（BX）缺乏的突變植物作用的土壤中生長的野生型（WT）植物（n=10）的根際和根細菌（a）和真菌（b）根際主坐標的偏典范分析（CAP）。使用Bray-Curtis距離的CAP分析樣本類型和土壤條件的限制。在圖上方標出了模型解釋的總方差和模型顯著性。軸報告受約束的軸解釋的總變化的比例。c，d用CAP的細菌操作分類單位（bOTU，c）和真菌OTU（fOTU，d）的雙標圖，顯示各個OTU對不同處理分離的貢獻。

Fig. 7 Benzoxazinoid-conditioning shapes soil microbiota in the roots and rhizosphere of the second plant generation. a, b Partial Canonical analysis of Principal Coordinates (CAP) of rhizosphere and root bacterial (a) and fungal (b) communities of wild type (WT) plants grown in soils previously conditioned by WT (BX+) or benzoxazinoid (BX)-deficient bx1 mutant plants (BX?) (n=10). CAP ordinations using Bray–Curtis distance were constrained for the factors sample type and soil conditioning. Model, explained fraction of total variance and model significance are indicated above the plots. Axes report the proportions of total variation explained by the constrained axes. c, d Biplot with the bacterial Operational Taxonomic Unit (bOTU, c) and fungal OTU (fOTU, d) scores of the CAP showing the contribution of individual OTUs to the separation of the different treatments.

結果描述：

為了了解BX土壤處理是否改變了響應植物的根相關菌群，研究人員分析了生長在BX+和BX?土壤中野生植物的根和根際樣品中的真菌和細菌群落。非限制性和限制性的縱坐標顯示，細菌和真菌在根際和根之間以及在BX+和BX?土壤之間有明顯的分離（圖7a, b，補充圖9a, b）。多變量統計證實了這些效應的顯著性（補充表7和8；參考補充數據3來深入分析BX土壤處理對微生物多樣性的影響）。放線菌門OTUs和屬于子囊菌綱和球囊菌門的OTUs的亞群對根和根際樣品的分離貢獻最大（圖7c, d）

To understand whether BX soil conditioning changes the rootassociated microbiota of the responding plants, we profiled fungal and bacterial communities in root and rhizosphere samples of WT plants grown in BX+ and BX? soils. Both unconstrained and constrained ordinations revealed a clear separation of bacteria and fungi between rhizosphere and roots and between BX+ and BX? soils (Fig. 7a, b, Supplementary Fig. 9a, b). Multivariate statistics confirmed the significance of these effects (Supplementary Tables 7 and 8; we refer to the Supplementary Data 3 for in depth analysis of BX soil conditioning effects on microbial diversity). Actinobacteria OTUs and a subset of OTUs belonging to Ascomycota and Glomeromycota contributed most strongly to the separation of root and rhizosphere samples (Fig. 7c, d).

• 總結

研究中同時對細菌和真菌測序分析，可以使用同類圖形并行展示，增加結果的全面和豐富性；

圖注中對圖形的描述比較清楚，在上文中已經加粗，供寫作參考；

同一個圖除了CAP還可采用雙標圖(Biplot)展示，增加結果樣本豐度度和看問題的角度。

• 結果：

在RDA分析中，2014年采集的土壤樣品中20個最高豐度門中的16個更高，而2015年采集的土壤樣品中只有3個更高，這與之前的結果一致。如圖7b所示，第一軸解釋方差為82.1％（p=.001；置換檢驗）和
第二軸解釋方差為0.9％（p=.992；置換檢驗），表明第一軸基于硝化微生物群落將所有環境變量和土壤樣品顯著分開。

In the RDA analysis, the relative abundance of 16 out of the 20 top phyla were higher in the soil samples collected in 2014, whereas only three of these phyla were higher in the soil samples collected in 2015, which was consistent with the previous results. As shown in Figure 7b, 82.1% (p=.001; permutation test) and 0.9% (p=.992; permutation test) of the cumulative variance of the relationship were explained by the first and second axes of RDA, respectively, which indicated that the first axis significantly separated all of the environmental variables and soil samples based on nitrifying microbial communities.

例3. RDA分析環境因子對微生物群落影響

本文由北京師范大學環境學院黃來斌(Laibin Huang)等人于2020年發表在Global Change Biology上，在人類活動強度高于全球沿海生態系統平均水平的情況下，分析了在大尺度生態系統解除脅迫因子后，微生物群落整體也能在短時間內恢復到與原始狀態相似的狀態，這意味著即使在強烈的人為干擾下，土壤微生物群落也具有很強的恢復能力。點擊訪問中文解讀

圖7. RDA分析對珠江三角洲沿岸微生物群落變化的環境解釋。（a） 前20門變化的環境解釋；（b）與硝化過程有關的屬變化的環境解釋。USF2014（粉色圓圈）和USF_2015（藍色圓圈）分別代表2014年和2015年上游的原位土壤；DSF（粉色鉆石）和DSF_（藍色鉆石）分別代表2014年和2015年下游的原位土壤（無限制置換檢驗，n=999）。

FIGURE7 Environmental explanation of the changes in microbial community along the Pearl River Estuary by RDA analysis. (a) Environmental explanation of the changes in top 20 phylum; (b) environmental explanation of the changes in genus related to nitrification process. USF_2014 (pink circle) and USF_2015 (blue circle) represent field soils from upstream in 2014 and 2015, respectively; DSF_2014 (pink diamond) and DSF_2015 (blue diamond) represent field soils from downstream in 2014 and 2015, respectively. .05; .01; <.001 (Unrestricted permutation test, n=999)

• 結果：

在RDA分析中，2014年采集的土壤樣品中20個最高豐度門中的16個更高，而2015年采集的土壤樣品中只有3個更高，這與之前的結果一致。如圖7b所示，第一軸解釋方差為82.1％（p=.001；置換檢驗）和
第二軸解釋方差為0.9％（p=.992；置換檢驗），表明第一軸基于硝化微生物群落將所有環境變量和土壤樣品顯著分開。

In the RDA analysis, the relative abundance of 16 out of the 20 top phyla were higher in the soil samples collected in 2014, whereas only three of these phyla were higher in the soil samples collected in 2015, which was consistent with the previous results. As shown in Figure 7b, 82.1% (p=.001; permutation test) and 0.9% (p=.992; permutation test) of the cumulative variance of the relationship were explained by the first and second axes of RDA, respectively, which indicated that the first axis significantly separated all of the environmental variables and soil samples based on nitrifying microbial communities.

例4. LDA展示腸道菌群四季的變化

本文為2017年8月25日發表在Sciences雜志上的封面文章為例，介紹了采集狩獵的非洲哈扎人腸道菌群存在一年四季的動態變化。詳細介紹見- 3分和30分文章差距在哪里？

圖 1. 哈扎人腸道菌群隨季節變化。(E) 線性判別分析 ?(LDA，一種限制性排序分析)展示按季節分組最大貢獻的OTUs；箭頭的長度和方向指示每個特征（OTUs）的標準化比例。

Fig. 1. Hadza gut microbial community compositions are cyclic and can be differentiated by season. (E) Linear discriminant analysis, a supervised learning approach that utilizes a linear combination of features to maximize the separation of classes, successfully separates the subseasons, except for the dry seasons.The length and direction of the arrows indicate the normalized scalings for each of the features (OTUs).

• 結果：

專門嘗試通過整合OTU的線性組合來區分群體的監督學習方法無法區分連續年份中的同一季節（旱季），從而支持了菌群周期性重新配置的規律（圖1E，圖S3A和表S2）。

A supervised learning approach that specifically attempts to distinguish groups by integrating a linear combination of OTUs was unable to differentiate the same season (dry) in sequential years, supporting the cyclic nature of the reconfiguration (Fig. 1E, fig. S3A, and table S2).

例5. CPCoA展示基因型對菌群的調控

本文是2019年5月中科院遺傳發育所白洋組與JIC的Anne Osbourn組合作在擬南芥代謝物調控根系微生物組領域取得重大突破，成果以長文形式發表于Science雜志，海南大學羅杰教授對本工作的意義進行點評。詳細報導請點擊下方鏈接：- Science：擬南芥三萜化合物特異調控根系微生物組* 專家點評。

圖4. 三萜通路的突變體特異調控的根系細菌類群。(A) 基于Bray-Curtis距離的有監督主坐標軸分析展示基因型對菌群調控的效應。用于分析的每組植物數量：Col-0（n = 12），thas-ko1（n = 12），thas-ko2（n = 9），thah-ko（n = 14），thao-ko（n = 13），thaa2-ko（n = 9）和thaa2-crispr（n = 12）。來自兩個獨立實驗（實驗1和2）的生物學重復（單個植物）分別用點和三角形表示。橢圓包括每種基因型的68％樣品。

Fig. 4 Modulation of specific root bacterial taxa in triterpene pathway mutants. (A) Constrained principal coordinate analysis (CPCoA) of Bray-Curtis dissimilarity showing triterpene mutant effects. Total number of individual plants used for analyses: Col-0 (n = 12), thas-ko1 (n = 12), thas-ko2 (n = 9), thah-ko (n = 14), thao-ko (n = 13), thaa2-ko (n = 9), and thaa2-crispr (n = 12). Biological replicates (individual plants) from two independent experiments (experiment 1 and 2) are indicated with dots and triangles, respectively. Ellipses include 68% of samples from each genotype.

• 結果：

有監督的主坐標軸分析(CPCoA)展示了野生型與突變體根系微生物組的差異(可以解析16.5%的總體變異，多元變量置換檢驗 P < 0.001，圖4A，附表18-22)。存在兩個獨立來源的突變體THAS (thas-ko1和thas-ko2) 和THAA2 (thaa2-ko和thaa2-cripsr)表現出相似的代謝缺陷（附圖18-22）也同樣擁有相似的微生物組成和多樣性（圖4A，附圖23，表23-26）。

Constrained principal coordinate analysis (CPCoA) revealed differences in root microbiota between the wild-type (Col-0) and mutant lines (16.5% of total variance was explained by the plant genotypes, P < 0.001, permutational multivariate analysis of variance) (Fig. 4A and tables S18 to S22). Pairs of independent mutant lines for THAS (thas-ko1 and thas-ko2) and THAA2 (thaa2-ko and thaa2-cripsr) each display similar metabolite defects (figs. S18 to S22) and have similar microbiota profiles and microbial diversity (Fig. 4A, fig. S23, and tables S23 to S26).

繪圖實戰

安裝和加載R包：amplicon

# 基于github安裝包，需要devtools，檢測是否存在，不存在則安裝 if (!requireNamespace("devtools", quietly=TRUE))install.packages("devtools") # 注：提示安裝源碼包的最新版時，推薦選否，加速安裝預編譯的穩定版。原碼包編譯時間長且容易出錯 # 第一次運行，會自動在C:\Users\User\Documents\R\win-library\4.0目錄中安裝75個包 # 加載github包安裝工具 library(devtools)# 檢測amplicon包是否安裝，沒有從源碼安裝 if (!requireNamespace("amplicon", quietly=TRUE))install_github("microbiota/amplicon") # 提示升級，選擇3 None不升級；升級會容易出現報錯 # library加載包，suppress不顯示消息和警告信息 suppressWarnings(suppressMessages(library(amplicon)))

輸入文件和基本參數

設置基本參數：讀取實驗設計、設置分組列名、圖片寬高

# 讀取元數據 metadata=read.table("metadata.txt", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F) # 設置實驗分組列名 group="Group" # 圖片16：10，即半版89 mm x 56 mm; 183 mm x 114 mm w=89 h=59

限制性PCoA

限制性PCoA(Constrained PCoA)，要求至少分組數量 >= 3時才可用，否則請跳過此步分析。

方法1. 基于特征OTU/ASV表使用vegan包計算矩陣矩陣并進行限制性主坐標分析。

library(amplicon) # 讀取特征ASV表，多樣性分析輸入抽平標準化的表 otutab=read.table("otutab_rare.txt", header=T, row.names=1, sep="\t", comment.char="") # 基于特征表、原數據、距離類型、分組列名、是否添加置信橢圓，是否添加樣本標簽 (p = beta_cpcoa(otutab, metadata, dis="bray", groupID="Group", ellipse=T, label=T)) ggsave("p1.cpcoa_otutab.pdf", p, width=w, height=h, units="mm") ggsave("p1.cpcoa_otutab.png", p, width=w, height=h, units="mm", dpi=300)

圖1. 基于Bray-Curtis距離的限制性主坐標分析不同基因型分組

我們輸入文件為更常用為各種類型的距離矩陣，可以基于OTU表使用QIIME的beta_diversity.py、USEARCH的-beta_div計算。

常用的矩陣類型常用Bray-Curtis、Jaccard、Unifrac和Unweighted unifrac等。下面以Unifrac距離矩陣為例繪制

distance_mat=read.table("unifrac.txt", header=T, row.names=1, sep="\t", comment.char="") distance_mat[1:3, 1:3] (p=beta_cpcoa_dis(distance_mat, metadata, groupID=group)) ggsave("p2.cpcoa_unifrac_distance.pdf", p, width=w, height=h, units="mm") ggsave("p2.cpcoa_unifrac_distance.png", p, width=w, height=h, units="mm", dpi=300)

圖2. 基于Unifrac距離的限制性主坐標分析不同基因型分組

RDA/CCA

加載代碼包和輸入文件

otu=read.delim("otutab_rare.txt",row.names=1) # 導入環境因子文件 env=read.delim("env.txt",row.names=1)

RDA 或CCA 模型的選擇原則：RDA是基于線性模型，CCA是基于單峰模型。一般會選擇CCA來做直接梯度分析。首先函數先做DCA(detrended correspondence analysis)分析，并提取結果中Lengths of gradient第一軸的大小。如果大于4.0，就應該選CCA；如果3.0-4.0 之間，選RDA 和CCA均可；如果小于3.0，RDA 的結果要好于CCA。

所以這里會有同時做RDA和CCA的情況。

result = RDA_CCA(otutab,metadata,ps=NULL,env=env,group="Group") #--提取圖 p=result[[1]] p ggsave("p3.RDA.pdf", p, width=w*1.5, height=h*1.5, units="mm") ggsave("p3.RDA.png", p, width=w*1.5, height=h*1.5, units="mm", dpi=300)

圖3. RDA分析結果圖。點代表樣本，箭頭代表環境因子。

#--提取作圖數據 plotdata=result[[2]] plotdata[1:3,1:3]# 提取帶有樣本標簽的圖形 p2=result[[3]] p2

檢測解釋變量對響應變量的解釋度和顯著性。

# 提取環境因子同群落的差異檢測 table=result[[4]] head(table)

Lengths ofgradient 的第一軸的大小在3.0-4.0 之間，會同時計算RDA和CCA，供選擇。

# 如果兩種模型都選擇，則列表中5，6就被激活了# result[[5]] result[[6]]

LDA

LDA部分的分析集成到BetaDiv函數中。

# 安裝Bioconductor的R包phyloseq if (!requireNamespace("BiocManager", quietly=TRUE))install.packages("BiocManager") suppressWarnings(suppressMessages(library(BiocManager))) if (!requireNamespace("phyloseq", quietly=TRUE))BiocManager::install("phyloseq") library(phyloseq)# 輸入抽平標準化的特征表、元數據、分組列名、距離類型、降維和統計方法 result=BetaDiv(otu=otutab_rare, map=metadata, group="Group",dist="bray", method="LDA") # 返回結果列表：標準圖，數據，標簽圖，成對比較結果，整體結果#提取排序散點圖(結果列表中的1) (p=result[[1]]) ggsave(paste0("p4.LDA.jpg"), p, width=w, height=h, units="mm") ggsave(paste0("p4.LDA.pdf"), p, width=w, height=h, units="mm")

圖4. LDA分析樣本微生物群落結構。本數據在LDA分析的結果還沒有PCoA按組分開明顯，更不如CPCoA。不同方法只是不同看問題的角度，適合的結果是不同方法嘗試結合背景知識綜合思考后選擇的。

# 提取出圖坐標 plotdata=result[[2]] plotdata[1:3,1:3] # 提取帶標簽排序散點圖 (p=result[[3]])

參考文獻

• Rafal Zgadzaj, Ruben Garrido-Oter, Dorthe Bodker Jensen, Anna Koprivova, Paul Schulze-Lefert & Simona Radutoiu. (2016). Root nodule symbiosis in Lotus japonicus drives the establishment of distinctive rhizosphere, root, and nodule bacterial communities. Proceedings of the National Academy of Sciences of the United States of America 113, E7996-E8005, doi: https://doi.org/10.1073/pnas.1616564113

• Lingfei Hu, Christelle A. M. Robert, Selma Cadot, Xi Zhang, Meng Ye, Beibei Li, Daniele Manzo, Noemie Chervet, Thomas Steinger, Marcel G. A. van der Heijden, Klaus Schlaeppi & Matthias Erb. (2018). Root exudate metabolites drive plant-soil feedbacks on growth and defense by shaping the rhizosphere microbiota. Nature Communications 9, 2738, doi: https://doi.org/10.1038/s41467-018-05122-7

• Laibin Huang, Junhong Bai, Xiaojun Wen, Guangliang Zhang, Chengdong Zhang, Baoshan Cui & Xinhui Liu. (2020). Microbial resistance and resilience in response to environmental changes under the higher intensity of human activities than global average level. Global Change Biology 26, 2377-2389, doi: https://doi.org/10.1111/gcb.14995

• Samuel A. Smits, Jeff Leach, Erica D. Sonnenburg, Carlos G. Gonzalez, Joshua S. Lichtman, Gregor Reid, Rob Knight, Alphaxard Manjurano, John Changalucha, Joshua E. Elias, Maria Gloria Dominguez-Bello & Justin L. Sonnenburg. (2017). Seasonal cycling in the gut microbiome of the Hadza hunter-gatherers of Tanzania. Science 357, 802-806, doi: https://doi.org/10.1126/science.aan4834

• Ancheng C. Huang, Ting Jiang, Yong-Xin Liu, Yue-Chen Bai, James Reed, Baoyuan Qu, Alain Goossens, Hans-Wilhelm Nützmann, Yang Bai & Anne Osbourn. (2019). A specialized metabolic network selectively modulates Arabidopsis root microbiota. Science 364, eaau6389, doi: https://doi.org/10.1126/science.aau6389

• GUSTA ME http://mb3is.megx.net/gustame/constrained-analyses/rda

• 排序分析PCA、PCoA、CA、NMDS、RDA、CCA等區別與聯系 https://www.omicsclass.com/article/148

責編：劉永鑫 中科院遺傳發育所

版本更新歷史

1.0.0，2020/6/27，文濤，初稿

1.0.1，2020/7/13，劉永鑫，LDA/CCA/CPCoA背景介紹

1.0.2，2020/7/17，席嬌，全文校對

1.0.3，2020/8/8，劉永鑫，全文修改

1.0.4，2020/8/12，文濤，添加LDA分析，修改RDA分析不依賴taxonomy文件

1.0.5，2020/8/12，劉永鑫，添加兩篇Science實例解讀LDA和CPCoA

送書

在上周的留言送書活動中，恭喜下面這位讀者獲得書籍“R語言數據高效處理指南”，請及時與生信寶典編輯（shengxinbaodian）聯系。

高級教程和生信基礎知識在生信寶典往期推文中都有，贈送的書籍可以作為一類延伸閱讀擴展知識面。

本期的留言主題是：請提出你想生信寶典分享的知識并說出為什么？歡迎轉發朋友圈并留言評論。

留言得贊前三者中選取最有緣者獲得下面由北京大學出版社贊助的書籍（聯系小編時請附上分享截圖），結果在下一期送書活動中公布：

《Python人工智能開發從入門到精通》主要介紹了Python進行人工智能開發所需的技術、基礎設施、核心理念、實施方法與流程，以及實戰操作應用。

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的送书 | 222Beta多样性限制性排序CPCoA/CCA/RDA/LDA的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。