1. 程式人生 > >QIIME 2使用者文件. 7差異丰度分析gneiss(2018.11)

QIIME 2使用者文件. 7差異丰度分析gneiss(2018.11)

文章目錄

前情提要

QIIME 2使用者文件. 7差異丰度分析gneiss

Differential abundance analysis with gneiss

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

此例項需要一些基礎知識,要求完成本系列文章前兩篇內容:《1簡介和安裝》《4人體各部位微生物組分析》

在本教程中,您將學習如何使用gneiss中的balances來進行差異丰度分析。我們將關注的主要問題是如何以組成連貫的方式識別差異豐富的分類群。

詳者注:balances批比值中的分子,對應中文中的平衡、均衡、天平或餘額都有相似,但又不完全一致。目前沒有標準翻譯,文中還是使用原詞 balances

組成性指的是處理比例的問題。為了考慮測序深度的差異,微生物的丰度通常被標準化為比例(例如,相對丰度)。正因為如此,準確地推斷出哪些微生物正在發生變化就變得具有挑戰性——因為比例加和為1,單個微生物的變化也會改變其餘微生物的比例。

考慮以下示例:

image
圖1. 相對丰度變化示意圖

在左邊,我們有十個物種的真實丰度,第一個物種在時間點1和時間點2之間加倍。當我們將這些按比例歸一化時,似乎所有的物種在這兩個時間點之間都發生了變化。單從比例上看,我們永遠不會意識到這個問題,實際上我們不能僅僅根據比例精確地確定哪些物種正在發生變化。

雖然我們不能確切地解決鑑別不同物種的問題,但我們可以放寬這個問題,並詢問微生物的哪些分割槽/部分(partitions)正在改變。在上面的例子中,如果我們計算第一種和第二種之間的比率,對於原始丰度和比例,這個比率在時間點1是1:1,在時間點2是2:1。這是balances

試圖解決的問題。我們可以把重點放在分類群(或分類組)之間的比率上,而不是放在單個分類群上,因為這些比率是由所觀察物種的真實丰度和所觀察的比例構成的。我們通常對這些比率進行對數轉換,以優化視覺化效果(“對數比”)。計算多個物種的balances(或ratios)的概念可以擴充套件到樹,如下面的示例所示。

image

圖2. 樹形圖展示原始值、相對丰度、比例間的關係

在左邊,我們定義一棵樹,其中每個樹枝尖對應一個分類單元,下面是第一個樣本中每個分類單元的比例。內部節點(即balances)定義了底下分類群之間的對數比率。右邊是同一棵樹,下面是不同樣本中每個分類群的比例。只有一個分類群數量變化,正如我們之前所觀察到的,所有分類群的比例都會改變,但是看看balances,只有包含紫色分類群的平衡才會改變。在這種情況下,balances不會改變,因為它只考慮紅色與分類群之間的比率。通過觀察balances而不是比例,我們可以通過限制觀察只關注給定balances內的分類群來消除一些差異。

這裡突出的問題是,我們如何構造balances樹來控制這種變異,並識別出有趣的分類群差異豐富的分割槽?在gneiss中,這有三種主要的方法:

  • 相關聚類(Correlation clustering)。如果我們沒有關於如何將生物體聚集的相關先驗資訊,我們可以根據它們彼此共生的頻率來將生物體分組在一起。這可以在correlation-clustering命令中獲得,併為ilr-hierarchical建立樹輸入檔案。
  • 梯度聚類(Gradient clustering)。使用元資料類別對在相似樣本型別中發現的分類群進行聚類。例如,如果我們要評估pH是否是驅動因子,我們可以根據觀察到的分類群的pH進行聚類,並觀察低pH生物與高pH生物的比例是否隨著pH的變化而變化。這可以在gradient-clustering命令中獲得,併為ilr-hierarchical建立樹輸入檔案。
  • 系統發育分析(Phylogenetic analysis)。也可以使用gneiss以外建立的系統發育樹(例如,rooted-tree.qza)。在這種情況下,您可以使用您的系統發育樹作為ilr-phylogenetic的輸入。

一旦我們有了一棵樹,我們就可以使用下面的等式來計算balances

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

image

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

在計算出balances後,可以進行標準統計過程,如方差分析和線性迴歸。我們將使用慢性疲勞綜合症(Chronic Fatigue Syndrome )資料集演示執行這些過程。

建立balances

Creating balances

在Giloteaux等人(2016)發表的慢性疲勞綜合徵資料集中,有87人,48名患者和39名健康對照組。本教程使用的資料是在Illumina MiSeq上使用地球微生物組專案高變區4 (V4) 16S rRNA基因測序方法產出的。

對於上文提到了兩種常用安裝方法,我們每次在分析資料前,需要開啟工作環境,根據情況選擇對應的開啟方式。

啟動工作環境並建立工作目錄

# 建立qiime2學習目錄並進入
mkdir -p qiime2
cd qiime2

# Miniconda安裝的請執行如下命令載入工作環境
source activate qiime2-2018.11

# 如果是docker安裝的請執行如下命令,預設載入當前目錄至/data目錄
# docker run --rm -v $(pwd):/data --name=qiime -it  qiime2/core:2018.11

# 建立本節學習目錄
mkdir qiime2-chronic-fatigue-syndrome-tutorial
cd qiime2-chronic-fatigue-syndrome-tutorial

實驗資料下載

注意:QIIME 2 官方測試資料部分儲存在Google伺服器上,國內下載比較困難。可使用代理伺服器(如藍燈)下載,或公眾號後臺回覆"qiime2"獲取測試資料批量下載連結,你還可以跳過以後的wget步驟

下載來源Google文件的實驗設計

wget \
  -O "sample-metadata.tsv" \
  "https://data.qiime2.org/2018.11/tutorials/gneiss/sample-metadata.tsv"

下載特徵表和物種註釋

wget \
  -O "table.qza" \
  "https://data.qiime2.org/2018.11/tutorials/gneiss/table.qza"
wget \
  -O "taxa.qza" \
  "https://data.qiime2.org/2018.11/tutorials/gneiss/taxa.qza"

首先,我們將定義我們想要構建balances的微生物的分割槽。同樣,有多種可能的方法來構建一個樹(即層級結構),它定義了我們要為其構建balances的微生物balances的分割槽。我們將在這個資料集中展示相關聚類( correlation-clustering)和梯度聚類(gradient-clustering)。

請注意,我們將執行的差異丰度技術將使用對數比轉換(log ratio transforms)。由於不允許取零的對數,下面的兩種聚類方法都包含一個預設的偽計數引數。這將用1替換表中的所有零,這樣我們就可以在轉換後的表上應用對數。
輸入表是原始計數表(FeatureTable[Frequency])。

選項1:相關性聚類

Option 1: Correlation-clustering

這個選項應該是您的預設選項。我們將通過Ward的分層聚類來採用無監督聚類,以獲得Principal Balances。從本質上講,這將使用Ward層次聚類來定義通常彼此共存的微生物的分割槽,這是由以下度量標準定義的。

D(x,y)=V[Ln(x/y)]

其中x和y代表兩種微生物在所有樣品中的比例。如果兩種微生物高度相關,那麼這個數量將縮小到接近零。然後,Ward層次叢集將使用這個距離度量來迭代地將相互關聯的微生物群聚集在一起。最後,我們得到的樹將突出顯示高層結構,並識別資料中的任何塊。

qiime gneiss correlation-clustering \
  --i-table table.qza \
  --o-clustering hierarchy.qza

輸出物件:

  • table.qza: 特徵表
  • taxa.qza: 物種註釋
  • hierarchy.qza: 層級結構

選項2:梯度聚類

Option 2: Gradient-clustering

相關聚類的另一種選擇是基於數字元資料類別建立樹。通過梯度聚類,我們可以對出現在元資料類別類似範圍內的分類群進行分組。在本例中,我們將使用元資料類別年齡建立樹(層次結構)。請注意,元資料類別不能缺少變數,並且必須是數字。

qiime gneiss gradient-clustering \
  --i-table table.qza \
  --m-gradient-file sample-metadata.tsv \
  --m-gradient-column Age \
  --o-clustering gradient-hierarchy.qza

輸出物件:

  • gradient-hierarchy.qza: 梯度層級

下游分析的一個重要考慮因素是過度擬合問題。當使用梯度聚類時,您正在建立一個樹來最好地突出所選元資料類別中的成分差異,並且可能會得到假陽性結果。小心使用梯度聚類。

用平衡建立線性模型

Building linear models using balances

現在我們有了一個定義分割槽的樹,我們可以執行等軸測對數比率(ILR)轉換。ILR轉換計算樹中每個節點上組之間的對數比率。

qiime gneiss ilr-hierarchical \
  --i-table table.qza \
  --i-tree hierarchy.qza \
  --o-balances balances.qza

輸出物件:

  • balances.qza 對數比結果

既然我們有了樹中每個節點的對數比率,我們就可以對balances進行線性迴歸。我們將要執行的線性迴歸稱為多變數響應線性迴歸(multivariate response linear regression),其歸根結底是對每個balances分別執行線性迴歸。

我們可以用它來嘗試將微生物丰度與環境變數聯絡起來。與標準單變量回歸相比,執行此模型有多方面的優勢,因為它可以避免與過度擬合相關的許多問題,並且可以讓我們根據環境引數瞭解群落範圍內的擾動。

由於微生物丰度可以直接對映到balances上,我們可以直接對balances進行多元迴歸。我們將要建立的模型表示如下

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

image

其中y是代表待預測的balances矩陣,βi→表示協變數i的係數向量,Xi→表示協變數i的測量。

記住,方差分析(ANOVA)是線性迴歸的一種特殊情況——方差分析可以解決的每個問題都可以重新表述為線性迴歸。詳情請參閱本帖。所以我們仍然可以用這種方法回答同樣的差異丰度問題,我們可以控制多個不同的潛在混淆變數和互動效應以獲得更精確的答案。

為了更好的理解以下公式,我們先看一下實驗設計的樣式

less -S sample-metadata.tsv
#SampleID       Sample_Name_s   BarcodeSequence LinkerPrimerSequence    Subject Sex     Age     Pittsburgh      Bell    BMI     sCD14ugml       LBPugml LPSpgml IFABPpgml       Physical_functioning    Role_physical   Role_emotional  Energy_fatigue  Emotional_well_being    Social_functioning      Pain    General_health  Description
ERR1331797      LR17    AAGGGACAAGTG    TATGGTAATTGT    Patient Female  67      2       40.0    32.89   1.65    15.01   279.3   145.0   40.0    0.0     0.0     35.0    52.0    25.0    23.0    30.0    
ERR1331823      IC19    TCCACCCTCTAT    TATGGTAATTGT    Control Female  43                      26.15   1.36    18.76   100.54  184.1                                                                   
ERR1331841      LR49    CTGTAGCTTGGC    TATGGTAATTGT    Control Female  60                      30.89   1.23    12.43   97.36   234.4   

線性迴歸分析:

qiime gneiss ols-regression \
  --p-formula "Subject+Sex+Age+BMI+sCD14ugml+LBPugml+LPSpgml" \
  --i-table balances.qza \
  --i-tree hierarchy.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization regression_summary.qzv

輸出物件:

  • regression_summary.qzv 迴歸摘要

image

圖5. 迴歸結果總結。說明見下方段落。

現在我們總結了迴歸模型。具體來說,我們想看看哪些協變數(covariates)對模型影響最大,哪些balances有意義,以及有多少潛在的過度擬合正在發生。

在迴歸摘要中,有幾點需要注意。摘要中有一個R2(Rsquared),它給出了有關回歸模型解釋群落中有多少方差的資訊。從我們看到的,迴歸可以解釋大約11%的群落結構變異(Rsquared: 0.1108)。這是典型的人類腸道微生物群,因為由於遺傳、飲食、環境等因素,有非常高的混雜變異。

為了評價一個單一協變數的解釋模型,採用了一個遺漏變數的方法。遺漏一個變數,計算r2的變化(R2diff)。變化越大,協變數的作用越強。在這種情況下,課題(Subject)是最大的解釋因素,解釋了2.41%的變化。

為了確保我們不過度擬合,進行了10倍交叉驗證。這將把資料分成10個分割槽,在其中9個分割槽上構建模型,並使用剩餘的分割槽來度量預測的準確性。此過程重複10次,每輪交叉驗證一次。每輪交叉驗證報告模型內誤差(within model error, mse)、R2和預測精度(R2 and the prediction accuracy, pred_mse)。在這裡,預測精度低於模型內誤差,表明不會發生過擬合

接下來,我們採用熱圖視覺化每個balances的所有協變數的顯著性水平(p值)。熱圖的列代表balances,熱圖的行代表協變數。熱圖由p值的負對數著色,突出潛在的顯著p值。使用滑鼠懸停功能可以獲得特定的係數值及其對應的p值,並啟用縮放可以探索感興趣的協變數和balances。

接下來是預測(prediction)和殘差(residual)圖。這裡,只繪製了前兩個balances,模型中的預測殘差被投影到這兩個balances上。從這些圖中我們可以看出,預測點與原始群落位於同一區域內。然而,我們可以看到,殘差與預測的方差大致相同。這有點讓人不安——但請注意,我們只能解釋10%的群落差異,所以這些計算很正常。

視覺化樹中的分支長度也可以通過模型中可解釋的平方和(sum of squares)進行標準化。最長的分支長度對應於最有用的balances。這可以讓我們對模型中最重要的balances進行高層次的概述。從這個圖和上面的熱圖中,我們可以看到balances y0很重要。這些balances不僅具有很小的p值(p<0.05),能用於區分受試者,而且它們在樹形圖中具有最大的分支長度。這表明微生物的這種劃分可以區分CFS患者和對照組。

我們可以在一個熱圖上視覺化這些balances,以瞭解它們代表的分類群。預設情況下,特徵表中的值是按對數標準化的,意味著樣本值以零為中心。

qiime gneiss dendrogram-heatmap \
  --i-table table.qza \
  --i-tree hierarchy.qza \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column Subject \
  --p-color-map seismic \
  --o-visualization heatmap.qzv

輸出視覺化結果:

  • heatmap.qzv 熱圖

image

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

具體來說,我們將繪製一個箱線圖,並確定能夠解釋對照組和患者組之間差異的分類群。

qiime gneiss balance-taxonomy \
  --i-table table.qza \
  --i-tree hierarchy.qza \
  --i-taxonomy taxa.qza \
  --p-taxa-level 2 \
  --p-balance-name 'y0' \
  --m-metadata-file sample-metadata.tsv \
  --m-metadata-column Subject \
  --o-visualization y0_taxa_summary.qzv

輸出視覺化結果:

  • y0_taxa_summary.qzv 熱圖

image

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

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

  1. 患者組與健康對照組y0分子的平均分類數增加。
  2. 患者組與健康對照組間y0分母的分類平均下降;
  3. 上述兩種情況的結合;
  4. y0分子和y0分母的類群丰度均增加,但y0分子的類群丰度比y0分子增加更多。
  5. y0分子和y0分母的類群丰度均下降,而y0分子比y0分母的類群丰度相對增加較多。

為了進一步縮小這些假設,需要生物先驗知識或實驗驗證。

Reference

  1. 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://doi.org/10.7287/peerj.preprints.27295v2

譯者簡介

劉永鑫,博士。2008年畢業於東北農大微生物學專業。2014年中科院遺傳發育所獲生物資訊學博士學位,2016年博士後出站留所工作,任巨集基因組學實驗室工程師,目前主要研究方向為巨集基因組學、資料分析與可重複計算和植物微生物組、QIIME 2專案參與人。發於論文12篇,SCI收錄9篇。2017年7月創辦“巨集基因組”公眾號,目前分享巨集基因組、擴增子原創文章400+篇,代表博文有《擴增子圖表解讀、分析流程和統計繪圖三部曲》,關注人數3.2萬+,累計閱讀500萬+。

猜你喜歡

寫在後面

為鼓勵讀者交流、快速解決科研困難,我們建立了“巨集基因組”專業討論群,目前己有國內外3000+ 一線科研人員加入。參與討論,獲得專業解答,歡迎分享此文至朋友圈,並掃碼加主編好友帶你入群,務必備註“姓名-單位-研究方向-職稱/年級”。技術問題尋求幫助,首先閱讀《如何優雅的提問》學習解決問題思路,仍末解決群內討論,問題不私聊,幫助同行。
image

學習擴增子、巨集基因組科研思路和分析實戰,關注“巨集基因組”
image

image

點選閱讀原文,跳轉最新文章目錄閱讀
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA