前情提要

• 文章導(dǎo)讀：QIIME 2可重復(fù)、交互和擴(kuò)展的微生物組數(shù)據(jù)分析流程

• 1簡介和安裝

• 2插件工作流程概述

• 3老司機(jī)上路指南

• 4人體各部位微生物組分析

• 5糞菌移植分析練習(xí)

• 5糞菌移植分析練習(xí)

• 6沙漠土壤分析Atacama soil

QIIME 2用戶文檔. 7差異豐度分析gneiss

Differential abundance analysis with gneiss

原文地址：https://docs./2018.11/tutorials/gneiss/

此實(shí)例需要一些基礎(chǔ)知識(shí)，要求完成本系列文章前兩篇內(nèi)容：《1簡介和安裝》和《4人體各部位微生物組分析》。

在本教程中，您將學(xué)習(xí)如何使用gneiss中的balances來進(jìn)行差異豐度分析。我們將關(guān)注的主要問題是如何以組成連貫的方式識(shí)別差異豐富的分類群。

詳者注：balances批比值中的分子，對(duì)應(yīng)中文中的平衡、均衡、天平或余額都有相似，但又不完全一致。目前沒有標(biāo)準(zhǔn)翻譯，文中還是使用原詞 balances。

組成性指的是處理比例的問題。為了考慮測序深度的差異，微生物的豐度通常被標(biāo)準(zhǔn)化為比例（例如，相對(duì)豐度）。正因?yàn)槿绱耍瑴?zhǔn)確地推斷出哪些微生物正在發(fā)生變化就變得具有挑戰(zhàn)性——因?yàn)楸壤雍蜑?，單個(gè)微生物的變化也會(huì)改變其余微生物的比例。

考慮以下示例：

圖1. 相對(duì)豐度變化示意圖

在左邊，我們有十個(gè)物種的真實(shí)豐度，第一個(gè)物種在時(shí)間點(diǎn)1和時(shí)間點(diǎn)2之間加倍。當(dāng)我們將這些按比例歸一化時(shí)，似乎所有的物種在這兩個(gè)時(shí)間點(diǎn)之間都發(fā)生了變化。單從比例上看，我們永遠(yuǎn)不會(huì)意識(shí)到這個(gè)問題，實(shí)際上我們不能僅僅根據(jù)比例精確地確定哪些物種正在發(fā)生變化。

雖然我們不能確切地解決鑒別不同物種的問題，但我們可以放寬這個(gè)問題，并詢問微生物的哪些分區(qū)/部分(partitions)正在改變。在上面的例子中，如果我們計(jì)算第一種和第二種之間的比率，對(duì)于原始豐度和比例，這個(gè)比率在時(shí)間點(diǎn)1是1:1，在時(shí)間點(diǎn)2是2:1。這是balances試圖解決的問題。我們可以把重點(diǎn)放在分類群（或分類組）之間的比率上，而不是放在單個(gè)分類群上，因?yàn)檫@些比率是由所觀察物種的真實(shí)豐度和所觀察的比例構(gòu)成的。我們通常對(duì)這些比率進(jìn)行對(duì)數(shù)轉(zhuǎn)換，以優(yōu)化可視化效果（“對(duì)數(shù)比”）。計(jì)算多個(gè)物種的balances（或ratios）的概念可以擴(kuò)展到樹，如下面的示例所示。

圖2. 樹形圖展示原始值、相對(duì)豐度、比例間的關(guān)系

在左邊，我們定義一棵樹，其中每個(gè)樹枝尖對(duì)應(yīng)一個(gè)分類單元，下面是第一個(gè)樣本中每個(gè)分類單元的比例。內(nèi)部節(jié)點(diǎn)（即balances）定義了底下分類群之間的對(duì)數(shù)比率。右邊是同一棵樹，下面是不同樣本中每個(gè)分類群的比例。只有一個(gè)分類群數(shù)量變化，正如我們之前所觀察到的，所有分類群的比例都會(huì)改變，但是看看balances，只有包含紫色分類群的平衡才會(huì)改變。在這種情況下，balances不會(huì)改變，因?yàn)樗豢紤]紅色與分類群之間的比率。通過觀察balances而不是比例，我們可以通過限制觀察只關(guān)注給定balances內(nèi)的分類群來消除一些差異。

這里突出的問題是，我們?nèi)绾螛?gòu)造balances樹來控制這種變異，并識(shí)別出有趣的分類群差異豐富的分區(qū)？在gneiss中，這有三種主要的方法：

• 相關(guān)聚類(Correlation clustering)。如果我們沒有關(guān)于如何將生物體聚集的相關(guān)先驗(yàn)信息，我們可以根據(jù)它們彼此共生的頻率來將生物體分組在一起。這可以在correlation-clustering命令中獲得，并為ilr-hierarchical創(chuàng)建樹輸入文件。

• 梯度聚類(Gradient clustering)。使用元數(shù)據(jù)類別對(duì)在相似樣本類型中發(fā)現(xiàn)的分類群進(jìn)行聚類。例如，如果我們要評(píng)估pH是否是驅(qū)動(dòng)因子，我們可以根據(jù)觀察到的分類群的pH進(jìn)行聚類，并觀察低pH生物與高pH生物的比例是否隨著pH的變化而變化。這可以在gradient-clustering命令中獲得，并為ilr-hierarchical創(chuàng)建樹輸入文件。

• 系統(tǒng)發(fā)育分析(Phylogenetic analysis)。也可以使用gneiss以外創(chuàng)建的系統(tǒng)發(fā)育樹（例如，rooted-tree.qza）。在這種情況下，您可以使用您的系統(tǒng)發(fā)育樹作為ilr-phylogenetic的輸入。

一旦我們有了一棵樹，我們就可以使用下面的等式來計(jì)算balances：

bi=rsr+s_logg(xr)g(xs)

其中i表示樹中的第i個(gè)內(nèi)部節(jié)點(diǎn)，g(x)表示集合x內(nèi)的幾何平均值，x r表示balances分子中的分類群豐富度集合，x s表示balances分母中的分類群豐富度集合，r和s分別表示x r和x s內(nèi)的分類群數(shù)量。

在計(jì)算出balances后，可以進(jìn)行標(biāo)準(zhǔn)統(tǒng)計(jì)過程，如方差分析和線性回歸。我們將使用慢性疲勞綜合癥(Chronic Fatigue Syndrome )數(shù)據(jù)集演示運(yùn)行這些過程。

創(chuàng)建balances

Creating balances

在Giloteaux等人(2016)發(fā)表的慢性疲勞綜合征數(shù)據(jù)集中，有87人，48名患者和39名健康對(duì)照組。本教程使用的數(shù)據(jù)是在Illumina MiSeq上使用地球微生物組項(xiàng)目高變區(qū)4 (V4) 16S rRNA基因測序方法產(chǎn)出的。

對(duì)于上文提到了兩種常用安裝方法，我們每次在分析數(shù)據(jù)前，需要打開工作環(huán)境，根據(jù)情況選擇對(duì)應(yīng)的打開方式。

啟動(dòng)工作環(huán)境并創(chuàng)建工作目錄

實(shí)驗(yàn)數(shù)據(jù)下載

注意：QIIME 2 官方測試數(shù)據(jù)部分保存在Google服務(wù)器上，國內(nèi)下載比較困難?？墒褂么矸?wù)器(如藍(lán)燈)下載，或公眾號(hào)后臺(tái)回復(fù)”qiime2”獲取測試數(shù)據(jù)批量下載鏈接，你還可以跳過以后的wget步驟。

下載來源Google文檔的實(shí)驗(yàn)設(shè)計(jì)

下載特征表和物種注釋

首先，我們將定義我們想要構(gòu)建balances的微生物的分區(qū)。同樣，有多種可能的方法來構(gòu)建一個(gè)樹（即層級(jí)結(jié)構(gòu)），它定義了我們要為其構(gòu)建balances的微生物balances的分區(qū)。我們將在這個(gè)數(shù)據(jù)集中展示相關(guān)聚類( correlation-clustering)和梯度聚類(gradient-clustering)。

請(qǐng)注意，我們將運(yùn)行的差異豐度技術(shù)將使用對(duì)數(shù)比轉(zhuǎn)換(log ratio transforms)。由于不允許取零的對(duì)數(shù)，下面的兩種聚類方法都包含一個(gè)默認(rèn)的偽計(jì)數(shù)參數(shù)。這將用1替換表中的所有零，這樣我們就可以在轉(zhuǎn)換后的表上應(yīng)用對(duì)數(shù)。
輸入表是原始計(jì)數(shù)表（FeatureTable[Frequency]）。

選項(xiàng)1：相關(guān)性聚類

Option 1: Correlation-clustering

這個(gè)選項(xiàng)應(yīng)該是您的默認(rèn)選項(xiàng)。我們將通過Ward的分層聚類來采用無監(jiān)督聚類，以獲得Principal Balances。從本質(zhì)上講，這將使用Ward層次聚類來定義通常彼此共存的微生物的分區(qū)，這是由以下度量標(biāo)準(zhǔn)定義的。

D（x，y）＝V[Ln(x/y)]

其中x和y代表兩種微生物在所有樣品中的比例。如果兩種微生物高度相關(guān)，那么這個(gè)數(shù)量將縮小到接近零。然后，Ward層次集群將使用這個(gè)距離度量來迭代地將相互關(guān)聯(lián)的微生物群聚集在一起。最后，我們得到的樹將突出顯示高層結(jié)構(gòu)，并識(shí)別數(shù)據(jù)中的任何塊。

輸出對(duì)象:

• table.qza: 特征表

• taxa.qza: 物種注釋

• hierarchy.qza: 層級(jí)結(jié)構(gòu)

選項(xiàng)2：梯度聚類

Option 2: Gradient-clustering

相關(guān)聚類的另一種選擇是基于數(shù)字元數(shù)據(jù)類別創(chuàng)建樹。通過梯度聚類，我們可以對(duì)出現(xiàn)在元數(shù)據(jù)類別類似范圍內(nèi)的分類群進(jìn)行分組。在本例中，我們將使用元數(shù)據(jù)類別年齡創(chuàng)建樹（層次結(jié)構(gòu)）。請(qǐng)注意，元數(shù)據(jù)類別不能缺少變量，并且必須是數(shù)字。

輸出對(duì)象:

• gradient-hierarchy.qza: 梯度層級(jí)

下游分析的一個(gè)重要考慮因素是過度擬合問題。當(dāng)使用梯度聚類時(shí)，您正在創(chuàng)建一個(gè)樹來最好地突出所選元數(shù)據(jù)類別中的成分差異，并且可能會(huì)得到假陽性結(jié)果。小心使用梯度聚類。

用平衡建立線性模型

Building linear models using balances

現(xiàn)在我們有了一個(gè)定義分區(qū)的樹，我們可以執(zhí)行等軸測對(duì)數(shù)比率（ILR）轉(zhuǎn)換。ILR轉(zhuǎn)換計(jì)算樹中每個(gè)節(jié)點(diǎn)上組之間的對(duì)數(shù)比率。

輸出對(duì)象:

• balances.qza 對(duì)數(shù)比結(jié)果

既然我們有了樹中每個(gè)節(jié)點(diǎn)的對(duì)數(shù)比率，我們就可以對(duì)balances進(jìn)行線性回歸。我們將要運(yùn)行的線性回歸稱為多變量響應(yīng)線性回歸(multivariate response linear regression)，其歸根結(jié)底是對(duì)每個(gè)balances分別執(zhí)行線性回歸。

我們可以用它來嘗試將微生物豐度與環(huán)境變量聯(lián)系起來。與標(biāo)準(zhǔn)單變量回歸相比，運(yùn)行此模型有多方面的優(yōu)勢，因?yàn)樗梢员苊馀c過度擬合相關(guān)的許多問題，并且可以讓我們根據(jù)環(huán)境參數(shù)了解群落范圍內(nèi)的擾動(dòng)。

由于微生物豐度可以直接映射到balances上，我們可以直接對(duì)balances進(jìn)行多元回歸。我們將要建立的模型表示如下

y→=β0→+βSubject→Xsubject→+βsex→Xsex→+βage→XAge→+βsCD14ugml→XsCD14ugml→+βLBPugml→XLBPugml→

其中y是代表待預(yù)測的balances矩陣，βi→表示協(xié)變量i的系數(shù)向量，Xi→表示協(xié)變量i的測量。

記住，方差分析(ANOVA)是線性回歸的一種特殊情況——方差分析可以解決的每個(gè)問題都可以重新表述為線性回歸。詳情請(qǐng)參閱本帖。所以我們?nèi)匀豢梢杂眠@種方法回答同樣的差異豐度問題，我們可以控制多個(gè)不同的潛在混淆變量和交互效應(yīng)以獲得更精確的答案。

為了更好的理解以下公式，我們先看一下實(shí)驗(yàn)設(shè)計(jì)的樣式

線性回歸分析：

輸出對(duì)象:

• regression_summary.qzv 回歸摘要

圖5. 回歸結(jié)果總結(jié)。說明見下方段落。

現(xiàn)在我們總結(jié)了回歸模型。具體來說，我們想看看哪些協(xié)變量(covariates)對(duì)模型影響最大，哪些balances有意義，以及有多少潛在的過度擬合正在發(fā)生。

在回歸摘要中，有幾點(diǎn)需要注意。摘要中有一個(gè)R2(Rsquared)，它給出了有關(guān)回歸模型解釋群落中有多少方差的信息。從我們看到的，回歸可以解釋大約11%的群落結(jié)構(gòu)變異(Rsquared:    0.1108)。這是典型的人類腸道微生物群，因?yàn)橛捎谶z傳、飲食、環(huán)境等因素，有非常高的混雜變異。

為了評(píng)價(jià)一個(gè)單一協(xié)變量的解釋模型，采用了一個(gè)遺漏變量的方法。遺漏一個(gè)變量，計(jì)算r2的變化(R2diff)。變化越大，協(xié)變量的作用越強(qiáng)。在這種情況下，課題(Subject)是最大的解釋因素，解釋了2.41%的變化。

為了確保我們不過度擬合，進(jìn)行了10倍交叉驗(yàn)證。這將把數(shù)據(jù)分成10個(gè)分區(qū)，在其中9個(gè)分區(qū)上構(gòu)建模型，并使用剩余的分區(qū)來度量預(yù)測的準(zhǔn)確性。此過程重復(fù)10次，每輪交叉驗(yàn)證一次。每輪交叉驗(yàn)證報(bào)告模型內(nèi)誤差（within model error, mse）、R2和預(yù)測精度（R2 and the prediction accuracy, pred_mse）。在這里，預(yù)測精度低于模型內(nèi)誤差，表明不會(huì)發(fā)生過擬合。

接下來，我們采用熱圖可視化每個(gè)balances的所有協(xié)變量的顯著性水平(p值)。熱圖的列代表balances，熱圖的行代表協(xié)變量。熱圖由p值的負(fù)對(duì)數(shù)著色，突出潛在的顯著p值。使用鼠標(biāo)懸停功能可以獲得特定的系數(shù)值及其對(duì)應(yīng)的p值，并啟用縮放可以探索感興趣的協(xié)變量和balances。

接下來是預(yù)測(prediction)和殘差(residual)圖。這里，只繪制了前兩個(gè)balances，模型中的預(yù)測殘差被投影到這兩個(gè)balances上。從這些圖中我們可以看出，預(yù)測點(diǎn)與原始群落位于同一區(qū)域內(nèi)。然而，我們可以看到，殘差與預(yù)測的方差大致相同。這有點(diǎn)讓人不安——但請(qǐng)注意，我們只能解釋10%的群落差異，所以這些計(jì)算很正常。

可視化樹中的分支長度也可以通過模型中可解釋的平方和(sum of squares)進(jìn)行標(biāo)準(zhǔn)化。最長的分支長度對(duì)應(yīng)于最有用的balances。這可以讓我們對(duì)模型中最重要的balances進(jìn)行高層次的概述。從這個(gè)圖和上面的熱圖中，我們可以看到balances y0很重要。這些balances不僅具有很小的p值（p<0.05），能用于區(qū)分受試者，而且它們在樹形圖中具有最大的分支長度。這表明微生物的這種劃分可以區(qū)分CFS患者和對(duì)照組。

我們可以在一個(gè)熱圖上可視化這些balances，以了解它們代表的分類群。默認(rèn)情況下，特征表中的值是按對(duì)數(shù)標(biāo)準(zhǔn)化的，意味著樣本值以零為中心。

輸出可視化結(jié)果:

• heatmap.qzv 熱圖

如圖例所示，每個(gè)balances的分子(numerators)以淺紅色突出顯示，而分母以深紅色突出顯示。從這里我們可以看到，y0的分母與y0的分子相比，幾乎沒有OTU。分母中的這些分類群可能重要，所以讓我們研究一下用balance_taxonomy構(gòu)成balance的分類。

具體來說，我們將繪制一個(gè)箱線圖，并確定能夠解釋對(duì)照組和患者組之間差異的分類群。

輸出可視化結(jié)果:

• y0_taxa_summary.qzv 熱圖

在這種特殊情況下，與對(duì)照組相比，患者組的對(duì)數(shù)比率更低。實(shí)質(zhì)上，這意味著與患者組相比，健康對(duì)照組y0分子中的分類群平均比y0分母中的分類群更豐富。

記住，根據(jù)本教程開頭給出的示例，不可能推斷出給定樣本中微生物的絕對(duì)變化。Balances將無法提供這類答案，但它可以限制可能出現(xiàn)情況的數(shù)量。具體來說，可能發(fā)生了以下五種情況之一。

1. 患者組與健康對(duì)照組y0分子的平均分類數(shù)增加。

2. 患者組與健康對(duì)照組間y0分母的分類平均下降；

3. 上述兩種情況的結(jié)合；

4. y0分子和y0分母的類群豐度均增加，但y0分子的類群豐度比y0分子增加更多。

5. y0分子和y0分母的類群豐度均下降，而y0分子比y0分母的類群豐度相對(duì)增加較多。

為了進(jìn)一步縮小這些假設(shè)，需要生物先驗(yàn)知識(shí)或?qū)嶒?yàn)驗(yàn)證。

Reference

Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet C, Al-Ghalith GA, Alexander H, Alm EJ, Arumugam M, Asnicar F, Bai Y, Bisanz JE, Bittinger K, Brejnrod A, Brislawn CJ, Brown CT, Callahan BJ, Caraballo-Rodríguez AM, Chase J, Cope E, Da Silva R, Dorrestein PC, Douglas GM, Durall DM, Duvallet C, Edwardson CF, Ernst M, Estaki M, Fouquier J, Gauglitz JM, Gibson DL, Gonzalez A, Gorlick K, Guo J, Hillmann B, Holmes S, Holste H, Huttenhower C, Huttley G, Janssen S, Jarmusch AK, Jiang L, Kaehler B, Kang KB, Keefe CR, Keim P, Kelley ST, Knights D, Koester I, Kosciolek T, Kreps J, Langille MG, Lee J, Ley R, Liu Y, Loftfield E, Lozupone C, Maher M, Marotz C, Martin BD, McDonald D, McIver LJ, Melnik AV, Metcalf JL, Morgan SC, Morton J, Naimey AT, Navas-Molina JA, Nothias LF, Orchanian SB, Pearson T, Peoples SL, Petras D, Preuss ML, Pruesse E, Rasmussen LB, Rivers A, Robeson, II MS, Rosenthal P, Segata N, Shaffer M, Shiffer A, Sinha R, Song SJ, Spear JR, Swafford AD, Thompson LR, Torres PJ, Trinh P, Tripathi A, Turnbaugh PJ, Ul-Hasan S, van der Hooft JJ, Vargas F, Vázquez-Baeza Y, Vogtmann E, von Hippel M, Walters W, Wan Y, Wang M, Warren J, Weber KC, Williamson CH, Willis AD, Xu ZZ, Zaneveld JR, Zhang Y, Zhu Q, Knight R, Caporaso JG. 2018. QIIME 2: Reproducible, interactive, scalable, and extensible microbiome data science. PeerJ Preprints 6:e27295v2 https:///10.7287/peerj.preprints.27295v2

譯者簡介

劉永鑫，博士。2008年畢業(yè)于東北農(nóng)大微生物學(xué)專業(yè)。2014年中科院遺傳發(fā)育所獲生物信息學(xué)博士學(xué)位，2016年博士后出站留所工作，任宏基因組學(xué)實(shí)驗(yàn)室工程師，目前主要研究方向?yàn)楹昊蚪M學(xué)、數(shù)據(jù)分析與可重復(fù)計(jì)算和植物微生物組、QIIME 2項(xiàng)目參與人。發(fā)于論文12篇，SCI收錄9篇。2017年7月創(chuàng)辦“宏基因組”公眾號(hào)，目前分享宏基因組、擴(kuò)增子原創(chuàng)文章400+篇，代表博文有《擴(kuò)增子圖表解讀、分析流程和統(tǒng)計(jì)繪圖三部曲》，關(guān)注人數(shù)3.2萬+，累計(jì)閱讀500萬+。