【原】Protein Cell：基于R語言的微生物組數(shù)據(jù)挖掘最佳流程（大眾評審截止14號晚6點）

宏基因組 2023-05-13 發(fā)布于北京

展開全文

Protein & Cell綜述：基于R語言的微生物組數(shù)據(jù)挖掘的最佳流程

近日，中國農(nóng)業(yè)科學院劉永鑫組聯(lián)合南京農(nóng)業(yè)大學袁軍組在國際期刊 Protein & Cell (IF = 15.3) 發(fā)表了題為”“The best practice for microbiome analysis using R”的綜述論文，建立了微生物組分析的R語言代碼庫EasyMicrobiomeR，助力微生物組學發(fā)展。

IF: 15.328
DOI: https:///10.1093/procel/pwad024
上線時間: 2023/05/02
第一作者：
文濤（Tao Wen），牛國慶（Guoqing Niu）
通訊作者：
袁軍Jun Yuan）（junyuan@njau.edu.cn)、劉永鑫（Yong-Xin Liu）（liuyongxin@caas.cn）
合作作者：
陳同（ Tong Chen），沈其榮（Qirong Shen）（shenqirong@njau.edu.cn)

大眾評審

同行評審通常只有2-4個審稿人，樣本量不足可能會遺漏一些問題。我們開展大眾評審模式，對同行評審的文章進行補充，希望進一步提高文章質(zhì)量。

文章目前處于提前在線發(fā)表階段，正式文檔正在排版校對中。歡迎廣大讀者下載文章Proof: http://www.ience/iMeta/temp/PROCEL_pwad024_wt.pdf ，提交修改建議至【金山文檔】文章大眾評審
https:///l/cnUufkyt5F5A 中幫忙文章質(zhì)量提高。

對于較重要的修改意見，將會在文章中致謝作者姓名。有興趣對本項目新增功能和貢獻代碼的用戶長期招募、歡迎加入，共同參與EasyMicrobiomeR下一版開發(fā)和共同發(fā)表，共同打造本領(lǐng)域基礎(chǔ)設(shè)施，推動微生物組學發(fā)展。

圖文摘要

本文詳細介紹了基于R語言的324個常用R包進行微生物組數(shù)據(jù)挖掘的過程。特別關(guān)注了六個微生物組分析集成R包（phyloseq, microbiome, MicrobiomeAnalystR, Animalcules, microeco, amplicon）。并提出了基于R語言進行微生物組分析的最佳流程，相關(guān)代碼都可以在 https://github.com/taowenmicro/EasyMicrobiomeR 獲取。如果覺得項目有用，歡迎點擊GitHub主頁右上角Star支持本項目。

摘要

隨著測序技術(shù)的逐步成熟，許多微生物組研究成果相繼發(fā)表，推動了相關(guān)分析工具的出現(xiàn)和發(fā)展。R語言是目前廣泛使用的用于微生物數(shù)據(jù)分析的平臺，具有強大的功能。然而，數(shù)以萬計的R包和無數(shù)類似的分析工具給許多研究人員挖掘微生物組數(shù)據(jù)帶來了重大挑戰(zhàn)。如何從眾多的R包中選擇合適、高效、方便、易學的工具已經(jīng)成為許多微生物組研究人員面臨的問題。我們整理了324個用于微生物組分析的常用R包，并根據(jù)應(yīng)用類別(多樣性、差異性、生物標志物、相關(guān)性和網(wǎng)絡(luò)、功能預測等)對其進行分類，以幫助研究人員快速找到用于微生物組分析的相關(guān)R包。此外，我們對微生物組分析的集成R包(Phyloseq、Microbiome、MicrobiomeAnalystR、Animalcules、Microeco和Amplicon)進行了系統(tǒng)的分類，并總結(jié)了它們的優(yōu)點和局限性，以幫助研究人員選擇合適的工具。最后，我們對用于微生物組分析的R包進行了全面的回顧，總結(jié)了微生物組中大多數(shù)常見的分析內(nèi)容，構(gòu)建最適合微生物組分析的流程。本文附帶了GitHub中數(shù)百個代碼的例子，可以幫助初學者學習，也可以幫助分析人員比較和測試不同的工具。本文對R在微生物組中的應(yīng)用進行了系統(tǒng)的梳理，為今后開發(fā)更好的微生物組工具提供了重要的理論依據(jù)和實踐參考。所有代碼都可以在GiHub：github.com/taowenmicro/EasyMicrobiomeR 上找到。

關(guān)鍵詞: R包，微生物組，數(shù)據(jù)分析，可視化，擴增子，宏基因組

文章亮點

本文按照六大功能類別（多樣性、差異性、生物標志物、相關(guān)性和網(wǎng)絡(luò)、功能預測和其他分析）對324個常用R包進行功能分類；
對微生物組分析的集成R包（phyloseq, microbiome, MicrobiomeAnalystR, Animalcules, microeco, and amplicon）進行了系統(tǒng)的介紹，并總結(jié)了其優(yōu)點和局限性；
本文總結(jié)了微生物組數(shù)據(jù)分析中常見的內(nèi)容，并提供了一套最適合微生物組數(shù)據(jù)挖掘的分析流程；
在GitHub上，分享了包含上萬行代碼的例子，這些例子不僅可以幫助初學者學習，還可以幫助專業(yè)人員比較和測試不同的工具。

引言

楊盛蝶翻譯，南京農(nóng)業(yè)大學，博士在讀；

宏基因組分析通過對微生物群落的DNA或RNA序列進行測序、定量、注釋和分析，來研究微生物的多樣性、結(jié)構(gòu)和功能。微生物組研究中常用的高通量測序技術(shù)主要有擴增子測序(amplicon sequencing)和隨機宏基因組測序(shotgun metagenomic sequencing)。擴增子測序以其成本低、分析體系成熟、分析過程簡單等優(yōu)點被廣泛應(yīng)用于微生物組研究中。隨機宏基因組測序提供了微生物的功能信息和更準確的微生物組成信息，測序成本較高，所需計算資源較大。我們在之前的綜述中系統(tǒng)地總結(jié)了這兩種測序的詳細流程(劉永鑫等，2021年，Protein & Cell)。微生物群落作為生物多樣性的重要組成部分，在生物學、生態(tài)學、生物技術(shù)、農(nóng)業(yè)和醫(yī)學等領(lǐng)域發(fā)揮著至關(guān)重要的作用。微生物群落分析需要多種生物信息學方法，主要包括三個部分：1)數(shù)據(jù)預處理，2)量化和注釋，3)統(tǒng)計和可視化(圖1A)。在預處理階段，對原始數(shù)據(jù)進行過濾和質(zhì)量控制，以確保數(shù)據(jù)質(zhì)量。在量化和注釋步驟中，使用工具和數(shù)據(jù)庫來識別微生物代表序列并注釋微生物分類和功能。微生物群落分析的前兩部分已經(jīng)進行了很好的討論，根據(jù)我們之前的論文(劉永鑫等，2023，iMeta)可以很好地完成這兩部分。最后，在統(tǒng)計和可視化步驟中，使用各種統(tǒng)計方法來探索微生物群落的多樣性、結(jié)構(gòu)和潛在功能。

圖1. 微生物群落數(shù)據(jù)分析工作流程及相關(guān)R軟件包。

(A) 微生物群落數(shù)據(jù)分析工作流程概述。核心文件是特征表（OTU）、分類表、樣品元數(shù)據(jù)（Metadata）、系統(tǒng)發(fā)育樹（Tree）和代表序列（Rep.fa）。(B) 詳細的微生物群落分析工作流程。首先，原始數(shù)據(jù)可以通過使用USEARCH/VSEARCH、QIIME2、DADA2軟件包進行處理。然后，將重要文件保存，用于RStudio軟件下的R語言環(huán)境開展下游分析。許多微生物分析方法都依賴于R語言開發(fā)的R包。詞云圖中的字體大小代表R包的引用次數(shù)。(C) 用于數(shù)據(jù)清理/處理和可視化的常用的R包。(D) 微生物群落分析R包劃分為六個類別。

隨著高通量測序技術(shù)的發(fā)展，許多研究利用擴增子技術(shù)(Thompson等，2017；Proctor等，2019)和隨機宏基因組測序技術(shù)(Carrión等，2019；Li等，2022；Paoli等，2022年)，這導致了微生物組分析方法學、軟件和流程的開發(fā)，例如QIIME(Caporaso等，2010年)、Mothur(Schlos等，2009年)、USEARCH(Edgar，2010)、VSEARCH (Rognes等，2016)、QIIME 2(Bolyen等，2019年)、Parallel‐Meta Suite (Chen等，2022)、EasyAmplicon(Liu等，2023)、Kraken(Wood和Salzberg等，2014)、MEGAN (Huson等，2007)、MetaPhlAn2(Truong等，2015)，HUMAnN2 (Franzosa et al.，2018)等。作為擴增測序數(shù)據(jù)分析最關(guān)鍵和最基本的步驟，OTU(操作分類單元)聚類方法在2015年前流行起來，而非聚類方法則是在最近幾年逐步發(fā)展和廣泛應(yīng)用的。目前，常見的非聚類方法包括DADA2(Callahan等，2016)、deblur(Amir等，2017)、unoise3(Edgar和FlyvbJerg，2015)。其中最具代表性的非聚類算法之一是DADA2，它是用R語言創(chuàng)建的。這使得R語言(Ihaka和Gentleman，1996)在擴增子測序的原始數(shù)據(jù)處理中占據(jù)了重要地位。與許多可用于微生物群落測序數(shù)據(jù)分析的上游步驟的軟件相比，下游的分析步驟嚴重依賴于R語言和各種軟件包。這些分析主要包括：1)多樣性分析；2)差異分析；3)相關(guān)性和網(wǎng)絡(luò)分析；4)生物標志物識別；5)功能預測；6)微生物群落與其他指標的整合分析(包括系統(tǒng)發(fā)育分析和環(huán)境因素分析等)。除了多元變量統(tǒng)計分析之外，R中還有各種數(shù)據(jù)清洗包允許在不同的分析之間轉(zhuǎn)換數(shù)據(jù)。

R是用于數(shù)據(jù)統(tǒng)計分析和可視化的自由、開源語言和環(huán)境，由新西蘭奧克蘭大學的Ross Ihaka和Robert Gentleman創(chuàng)建，現(xiàn)在由“R Development Core Team”負責。與SPSS、MINITAB, MATLAB等更適合處理標準化數(shù)據(jù)的分析工具相比，R語言既能處理已經(jīng)被修改的數(shù)據(jù)，也能處理原始數(shù)據(jù)。R可以很容易地實現(xiàn)幾乎所有的分析方法，許多最新的方法或算法都是最先在R上展示的。此外，R顯示了出色的數(shù)據(jù)可視化，特別是對于復雜數(shù)據(jù)。強大而靈活的交互分析也是R的一個優(yōu)勢，同時支持可視化的數(shù)據(jù)探索。R語言的功能在很大程度上依賴于數(shù)以千計的R包，它們提供了各種各樣的數(shù)據(jù)處理和分析的策略，幾乎允許在R中完成任何數(shù)據(jù)分析過程。CRAN上發(fā)布的R包總數(shù)為18,981個，Bioconductor上發(fā)布的R包總數(shù)為2,183個(截至2023年1月31日)。這些程序包展示了R強大的數(shù)據(jù)處理和分析性能。

近年來，在R平臺上開發(fā)了大量用于微生物組下游分析的R包，為相關(guān)研究領(lǐng)域做出了重要貢獻。然而，下游分析R包的數(shù)量已經(jīng)達到了令人眼花繚亂的水平(圖1B)。此外，包含大量微生物組分析的集成R包逐漸出現(xiàn)，如phyloseq(McMurdie和Holmes，2013)、microeco(Liu等，2020)和amplicon(Liu等，2023)。如此豐富的R包為微生物組分析人員提供了更多的選擇，但也使得在許多類似的分析工具中確定最合適的工具存在難度。此外，過多的R包使初學者很難走上一條組織良好的微生物組分析學習之路。因此，當務(wù)之急是比較相似的分析功能，并提取其異同點，選擇最適合微生物組分析的流程，幫助初學者更有效地學習。

本文試圖對324個常見的R包(附圖1)進行分類和運行。特別是用于微生物組分析的集成R包，并完成以下三個部分：1)根據(jù)微生物組分析的功能類別比較不同的R包分析流程，分析結(jié)果并總結(jié)范例代碼；2)根據(jù)微生物組分析的功能類別介紹六個集成R包的內(nèi)容，比較分析結(jié)果，生成范例代碼；3)根據(jù)所有R包，使用R語言選擇最優(yōu)的分析方法，并提供范例代碼供研究人員參考和學習。

微生物數(shù)據(jù)分析前準備工作

微生物組下游分析需要準備五個數(shù)據(jù)文件，包括特征表、特征注釋文件、樣本分組文件、系統(tǒng)發(fā)育樹和代表性序列。對于初學者來說，重要的是了解這些文件的格式和基本數(shù)據(jù)結(jié)構(gòu)，并學習如何將這些文件導入R語言。此外，不同的分析內(nèi)容往往對數(shù)據(jù)有不同的要求，需要學習一些數(shù)據(jù)操作技巧來滿足各種功能的需求。最后，有必要學習R出圖的基礎(chǔ)知識，以便于展示結(jié)果。

數(shù)據(jù)準備和清洗

在對序列進行預處理、量化和注釋后，需要對輸出的文件進行進一步的分析，包括導入這些文件、清理數(shù)據(jù)、轉(zhuǎn)換格式，這是后續(xù)R中微生物組分析所需的。在進行統(tǒng)計分析之前，必須掌握R語言的基本程序，以應(yīng)對不同R包的數(shù)據(jù)輸入要求。這一部分包括：數(shù)據(jù)的導入、整理、過濾、基本計算、轉(zhuǎn)換、均一化和修改。原始數(shù)據(jù)處理中經(jīng)常使用五種數(shù)據(jù)形式，包括特征表(文件格式為.csv/.txt/.xlsx/.biom，通常為分類和功能表，包括OTU/ASV/分類/基因/模塊/路徑表)、特征注釋(.csv/.txt/.xlsx/.biom)、樣本分組文件(.csv/.txt)、進化/系統(tǒng)發(fā)育樹(.nwk/.tree)、代表性序列(.fast a/.fas/.fa)。所有與數(shù)據(jù)清理相關(guān)的包如圖1C所示。微生物群落的表格數(shù)據(jù)輸入主要是使用utils包(代碼1A，GitHub中的腳本)中的read.table()、read.delim()和read.csv()等函數(shù)完成的。進化樹文件的讀取依賴于ape/ggtree/treeio包中的read.tree()或Phyloseq包中的read_tree()之類的函數(shù)。為了讀取微生物組中具有代表性的序列文件，通常使用Biostrings包(Pages等，2016)中的readDNAStringSet()。目前，微生物組的大數(shù)據(jù)集成已成為一種趨勢，并導致出現(xiàn)了用于整合多項研究數(shù)據(jù)的R包，如curatedMetagenomicData(Pasolli等，2017年)。只需要導入該R包，就可以重新分析篩選的數(shù)據(jù)，而不是輸入的原始數(shù)據(jù)。

數(shù)據(jù)整理的本質(zhì)可以概括為三個步驟：數(shù)據(jù)分割、使用函數(shù)處理以及將輸出結(jié)果合并為所需的格式。R中基礎(chǔ)包的函數(shù)組合在一起，可以滿足微生物組數(shù)據(jù)操作的大部分要求。例如，結(jié)合基本統(tǒng)計函數(shù)[sum()、mean()、sd()等]的“for循環(huán)”，可用于進行微生物相對豐度(代碼1B)的基本統(tǒng)計分析和數(shù)據(jù)轉(zhuǎn)換；base包提供apply系列函數(shù)，包括apply()、sapply()、lapply()、Tapply()、aggregate()等，可以快速完成數(shù)據(jù)處理的三個階段。apply家族函數(shù)提供了一個框架，起到替代for循環(huán)的作用，比R(代碼1B)中的基本“for loop”函數(shù)快得多。類似的purr(https://github.com/tidyverse/purrr))包也可以代替“for loop”來執(zhí)行高效的操作。

plyr包在base包的基礎(chǔ)上升級了數(shù)據(jù)框，列表等多種數(shù)據(jù)形式的數(shù)據(jù)處理過程。plyr包提供了在一個函數(shù)內(nèi)同時完成“Split - Apply - Combine”三個數(shù)據(jù)處理階段，并且，plyr包實現(xiàn)R類型（向量、列表和數(shù)據(jù)框）之間的分組變換，基本上可以取代base包中的apply家族函數(shù)。可以很方便的處理分組計算，例如，計算不同分類水平下的微生物豐度（示例二）。reshape2包在數(shù)據(jù)處理過程中提供了長寬格式轉(zhuǎn)化。由于ggplot2包大多數(shù)建模函數(shù)如lm()、glm()、gam()等會使用的是長數(shù)據(jù)格式，而微生物組數(shù)據(jù)一般都是寬數(shù)據(jù)格式，所以使用reshape2可以完成微生物組數(shù)據(jù)到繪圖的轉(zhuǎn)化。

dplyr包是tidyverse家族的一員，創(chuàng)新性地放棄了R中數(shù)據(jù)保存的常見形式，使用tibble格式(比data.Frame格式更強大)進行數(shù)據(jù)處理，可以更高效地完成行和列內(nèi)的數(shù)據(jù)框選擇、合并和統(tǒng)計，以及數(shù)據(jù)框長度和寬度格式的變化，配合%>%管道符號可以完成更復雜的數(shù)據(jù)處理過程。Tibble格式可以在分析和建模過程中存儲數(shù)據(jù)，這對數(shù)據(jù)分析非常重要。例如，我們演示了使用dplyr和流程來運行隨機森林建模以及重要變量的挑選過程(代碼1E)。

R語言中的可視化

牛國慶翻譯，南京農(nóng)業(yè)大學，博士在讀；

在大多數(shù)情況下，微生物組數(shù)據(jù)分析中，我們通常習慣于繪制標準圖形，例如alpha/beta多樣性和分類組成。圖1C中顯示是與可視化相關(guān)的R包。由于ggplot2的廣泛使用（代碼2A），出現(xiàn)了許多基于ggplot2擴展的R包，具有豐富的繪圖樣式、顏色和主題。這些R包主要包括用于繪制三元圖（Hamilton和Ferry，2018）的ggtern（代碼2B）、用于繪制網(wǎng)絡(luò)圖（Si等，2022）的ggraph（代碼2C）、用于繪制進化樹或物種分類圖（Xu等，2022）的ggtree（代碼2D）、ggalluvial包、ggVennDiagram包（代碼2E），用于繪制餅圖ggstatsplot包以及提供各種不同主題和顏色的輸出的ggpubr包。此外，基于grid繪圖系統(tǒng)的pheatmap和ComplexHeatmap軟件包（GU，2022）可以繪制不同樣本中特征的相對豐度（代碼2F），VennDiagram包（Chen和Boutros，2011）可以顯示不同樣本中特征的數(shù)量。UpSetR包（Conway等，2017）可以繪制了一種類似于Venn圖的新型圖形，稱為Upset視圖。基于base的繪圖系統(tǒng)雖然比較復雜并且難于學習，但是功能卻非常強大，對于繪制復雜圖形（如由微生物組成的和弦圖），是一個不錯的選擇，例如使用circlize（Gu等，2014）包（代碼2G）繪制微生物組成的弦圖。

此外，經(jīng)常有許多涉及圖形組合的微生物組繪圖工作。目前，R中有許多可以組合圖形的工具，如cowplot、patchwork和aplot。patchwork包具有最強大的功能，支持模塊化拼接圖形（代碼2H）。

微生物群落分析

我們在圖1D中將微生物組數(shù)據(jù)分析分為以下六個主要類型：多樣性分析、差異分析、生物標志物鑒定、相關(guān)性和網(wǎng)絡(luò)分析、功能預測以及其他微生物組分析（包括朔源分析、群落組裝過程以及微生物群落與環(huán)境因素之間的關(guān)聯(lián)分析）。然后，我們將會對所有相關(guān)的R包進行組織、比較和總結(jié)。

多樣性分析

微生物群落多樣性主要包括Alpha多樣性（Richness、Shannon、Simpson、Chao1、ACE等）、稀釋曲線、Beta多樣性（排序和聚類分析）、分類或功能組成。這里必須介紹一下vegan包（Oksanen等，2007），它是由包括來自芬蘭的Oksanen在內(nèi)的九位數(shù)量生態(tài)學家編寫的，是一款最初用于處理群落生態(tài)學數(shù)據(jù)的工具。該包提供了各種數(shù)據(jù)標準化和轉(zhuǎn)換方法。例如，用于Alpha多樣性分析的數(shù)據(jù)可以使用rrarefy()函數(shù)在相同測序深度下進行標準化，用于排序分析的數(shù)據(jù)可以使用decostand()函數(shù)進行標準化處理（代碼3A）。測序數(shù)據(jù)經(jīng)過樣本歸一化后，多樣性計算結(jié)果會更合理。此外，還可以使用adespatial（Dray等，2018）、ade4（Dray和Dufour，2007）和picante包（Kembel等，2010）來計算Alpha多樣性指標。例如，可以使用picante包中的pd()函數(shù)計算系統(tǒng)發(fā)育多樣性（代碼3A）。vegan不僅可以進行Alpha多樣性分析，還提供了一些其他功能，例如rda()函數(shù)用于主成分分析（PCA）和冗余分析（RDA），cca()函數(shù)用于對應(yīng)分析（CA）和典型對應(yīng)分析（CCA），decorana()函數(shù)用于決策曲線分析（DCA），metaMDS()函數(shù)用于行微生物組排序分析的非度量多維尺度分析（NMDS）（代碼3B）。stats包中的prcomp()函數(shù)可用于主成分分析（PCA），這是一種降維分析方法。MASS包中提供的mca()函數(shù)和FactoMineR包中提供的MCA()函數(shù)可用于多重對應(yīng)分析（代碼3B）；ape包提供的pcoa()函數(shù)可用于主坐標分析（PCoA）；MASS包提供的lda()函數(shù)可用于線性判別分析（LDA，代碼3C）。在運行多個排序操作之前，通常需要進行群落聚類。vegan包中的vegdist()函數(shù)可計算euclidean，manhattan, bray, canberra等多種生態(tài)距離（代碼3B）。此外，距離計算也可以使用stats包中的dist()函數(shù)完成。除了排序分析，距離矩陣還可用于聚類分析。stats包中的hclust()函數(shù)可用于聚類分析，factoextra包、kmeans包也可以實現(xiàn)類似的功能（代碼3D）。微生物組成分析主要用于顯示微生物的豐度，需要使用dplyr包對數(shù)據(jù)進行整理，然后使用ggplot2包進行展示。

差異分析

差異分析分為群落水平分析和特征水平（包括任何分類和功能層次）分析。群落水平的差異分析主要使用函數(shù)包括：vegan包中的adonis()、anosim()、mrpp()函數(shù)以及ape包中的mantel.test()函數(shù)（代碼4A）。特征水平的差異分析中組成差異可以利用stats包中的wilcox.test()函數(shù)（代碼4B）和t.test()函數(shù)（代碼4C）。隨后，針對測序數(shù)據(jù)專門開發(fā)了數(shù)據(jù)矯正算法，例如edgeR包中的：Upper Quartile(UQ)，Trimmed Mean of M-values (TMM)，Relative Log Expression (RLE)；DESeq2包中的Median（MED）以及metagenomeSeq（https://github.com/sirusb/metagenomeSeq）
包中Cumulative-Sum Scaling （CSS）算法（代碼4F）。此外，ALDEx2包提供了多項式模型，可用于推斷特征豐度并使用非參數(shù)檢驗、t檢驗或廣義線性模型計算特征差異（代碼4G）。ANCOM-BC包試圖通過對數(shù)線性模型矯正偏差來解決樣本異質(zhì)性問題。此外，用于微生物組數(shù)據(jù)矯正和差異檢驗的其他R包括limma（代碼4H）、DR、ANCOM（Lin和Peddada，2020）（代碼4I）、corncob（代碼4J）、Maaslin2（代碼4K）等。Nearing等（2022）表明，他們比較了這些差異分析方法，并提出ALDEx2和ANCOM-II（anchom_v2.1.R，代碼4L）是微生物群落差異分析中最好的方法。至于顯著性檢驗，不同的包使用不同的方法進行顯著性檢驗。例如，edgeR包使用Fisher檢驗，DESeq2和corncob包使用Wald檢驗，limma包使用t檢驗。還有其他的顯著性檢驗方法，如Wilcoxon秩和檢驗（ALDEx2和ANCOM-II），方差分析（ANOVA）（Maaslin2）等。

生物標志物判斷

特征微生物的尋找是為了解釋某些問題，如肥胖或高血壓人群的腸道生物標志物，或鐮刀菌枯萎病發(fā)病土壤的生物標志物等。通過差異分析挑選的微生物通常無法確定它們是否是關(guān)注的主要差異。因此，權(quán)重分析或機器學習方法被用來進一步區(qū)分特征微生物。

常用于權(quán)重分析的主要方法包括線性判別分析（LEfSe）、主成分分析（PCA）等（代碼5A）。LEfSe是專門針對微生物組數(shù)據(jù)開發(fā)的，其核心功能是用LDA（Fisher，1936年）和MASS（Ripley等，2013年）等包來實現(xiàn)。通過提取PCA排序的loading矩陣，可以找到對樣本變化影響最大的微生物作為生物標志物（代碼5B）。

在機器學習方面，目前微生物組分析中廣泛使用的隨機森林模型是通過使用randomforest包（Liaw和Wiener，2002年）來實現(xiàn)（代碼5C）。還有許多基于決策樹的機器學習模型，例如mboost（Hofner等，2014年）包提供了基于提升算法的方法，e1071（Dimitriadou等，2008年）包提供了支持向量機svm()（代碼5D），以及樸素貝葉斯（naiveBayes()）。xgboost包可以將多個樹模型集成在一起形成強大的分類器，還可以通過多種策略，包括正則化項、Shrinkage and Column Subsampling等，防止過擬合等問題。此外，使用pROC（Robin等，2011年）包用來繪制ROC曲線（代碼5D）以評估機器學習模型的效果。Caret包提供交叉驗證用來確定特征微生物數(shù)量（Kuhn，2008年）。目前，Wirbel等人（2021年）開發(fā)了一個開源的R包SIAMCAT，這是一個專門針對微生物組數(shù)據(jù)的特點進行定制，功能強大且友好的計算機器學習工具包。

功能預測

針對16S rDNA的功能預測，已開發(fā)了Tax4Fun（A?hauer等，2015年）包（代碼7A），使用擴增子數(shù)據(jù)更準確地預測微生物群落功能變化。該包已更新為Tax4Fun2（Wemheuer等，2020年）。Microeco可以實現(xiàn)對細菌/古菌的FAPROTAX（Louca等，2016年）預測和對真菌的FUNGuild（Nguyen等，2016年）預測，該預測是基于精選發(fā)表論文中分類功能描述數(shù)據(jù)庫。功能預測可以預測微生物群落功能并進行后續(xù)的統(tǒng)計分析。此外，vegan包可用于多樣性分析，而edgeR、DESeq2和limma包可用于差異分析。對于功能富集，clusterProfiler（代碼7B）包可進行GO、KEGG、GSEA和GSVA富集分析，它考慮了基因/通路的豐度，被推薦使用。此外，clusterProfiler包基于ggplot語法提供繪圖功能，可以簡單地繪制圖形?？梢允褂?strong>WGCNA進行基因/通路網(wǎng)絡(luò)分析，ggClusterNet進行網(wǎng)絡(luò)參數(shù)計算和可視化。然而，功能預測結(jié)果的可靠性，特別是對于環(huán)境樣品而言，目前存在爭議，因此通常需要進一步驗證分析結(jié)果。

其他微生物組分析

微生物群落形成過程的分析常用Stegen等人（2013）提出的框架，通常結(jié)合使用R包minpack.lm、picante、Hmisc、eulerr、FSA、ape、stats4等計算βNTI和RC-Bray指數(shù)（代碼8A），并推斷群落形成過程。Ning等人（2020年）通過基于系統(tǒng)發(fā)育分箱的零模型分析定量推斷群落構(gòu)建機制，并開發(fā)了R包iCAMP（代碼8B）。它可以定量評估不同生態(tài)過程（如同質(zhì)化選擇、異質(zhì)化選擇、擴散和漂移）對整個群落和每個系統(tǒng)發(fā)育分箱（通常由具有不同生態(tài)特征的單科或單目分類群組成）的相對重要性。此外，該包還提供中性理論模型，在群落和分類水平上進行系統(tǒng)發(fā)育和分類學零模型分析，計算類群間的生態(tài)位差異和系統(tǒng)發(fā)育距離，以及檢驗單個系統(tǒng)發(fā)育分箱內(nèi)的系統(tǒng)發(fā)育信號。

微生物群落經(jīng)常與環(huán)境指標進行相關(guān)性分分析，例如，使用vegan包提供的mantel.test()來檢驗微生物群落與環(huán)境指標的相關(guān)性，使用wascores()、mantel.correlog()來檢測微生物群落與環(huán)境因素之間的系統(tǒng)發(fā)育信號（代碼8C）。此外，ggClusterNet包可用于計算微生物/微生物群落與環(huán)境因子之間的共現(xiàn)關(guān)系，并生成可供發(fā)表的圖片（代碼8D）。

Knights等人（2011）提出了基于R語言的微生物組朔源工具source tracker。Metcalf等（2016）預測了死亡時間并追蹤了真實尸體上微生物群落的來源微生物，隨后微生物溯源分析逐漸普及。Shenhav等（2019）在R語言中提出了一種新的算法FEAST（代碼8E），使微生物溯源分析更加高效快捷，且誤報率低。

微生物組集成R包

謝鵬昊翻譯，南京農(nóng)業(yè)大學，博士在讀；

隨著微生物組測序的日益普及，專用于微生物組數(shù)據(jù)處理的R包也逐漸出現(xiàn)（圖2）。McMurdie和Holmes（2013）開發(fā)了phyloseq包：這是一個全面的微生物組數(shù)據(jù)處理工具（包括特征表、系統(tǒng)發(fā)育樹和特征注解）用于聚類，整合數(shù)據(jù)導入、存儲、分析和輸出。該包利用R中的許多工具進行生態(tài)和系統(tǒng)發(fā)育分析（如vegan, ade4, ape, 和picante），并使用ggplot2來輸出高標準的圖形。數(shù)據(jù)存儲結(jié)構(gòu)使用類似S4的存儲系統(tǒng)將所有相關(guān)數(shù)據(jù)存儲為單一的對象，從而使數(shù)據(jù)共享和分析復現(xiàn)更加容易。隨后，microbiome（https://github.com/microbiome/microbiome）
、MicrobiomeAnalystR（Chong等，2020）、microViz（Barnett等，2021）和micreobiomeSeq等包在此框架下出現(xiàn)。接下來，根據(jù)R6框架開發(fā)的microeco包提供了更多的分析功能。隨著數(shù)據(jù)交互分析的需求，Animalcules（Zhao等，2021）出現(xiàn)了。EasyMicroPlot（https://github.com/xielab2017/EasyMicroPlot）
也使用交互式界面進行微生物組數(shù)據(jù)的探索，允許快速進行微生物組的下游分析（圖3；表1）。

使用Phyloseq包進行微生物組數(shù)據(jù)分析

Phyloseq利用S4類對象，更適合面向?qū)ο蟮木幊蹋瑢ξ⑸锝M數(shù)據(jù)分析產(chǎn)生了巨大影響（圖2/3，圖S2A-G， Pipeline 1. phyloseq.Rmd）。通過S4類對象，phyloseq允許數(shù)據(jù)的五個部分（特征表、注釋信息、分組信息、代表性序列和進化樹）在同一框架下保持對應(yīng)，并提供了關(guān)于微生物特征和樣本的多種過濾功能，允許五部分數(shù)據(jù)在不考慮數(shù)據(jù)差異的情況下保持一致的過濾。它還通過微生物數(shù)據(jù)過濾和標準化，提供了微生物組多種分析，包括多樣性計算（圖S2A-B），微生物組成可視化（圖S2C-D），進化樹可視化，和網(wǎng)絡(luò)分析（圖S2E）。其Beta多樣性功能提供了超過30種距離算法，遠超過vegan等包提供的。其次，phyloseq包使用ggplot進行圖形可視化（圖S2F），這使得生成和修改圖形變得更加容易。此外，phyloseq可以在樹的枝和葉上整合進化樹和物種分類和豐度（圖S2G），使得樹信息豐富且美觀。

使用microbiome包進行微生物組數(shù)據(jù)分析

microbiome包也像phyloseq一樣使用S4類對象，也可以執(zhí)行大部分微生物組的分析（圖2/3，圖S3A-G， Pipeline 2. Microbiome.Rmd）。它包括微生物多樣性分析（圖S3A-E），以及差異分析（圖S3F-G）。與phyloseq相比，microbiome包在alpha多樣性指標上更為豐富，提供了超過30種alpha多樣性指標。其次，它提供了核心微生物計算和可視化功能?？偟膩碚f，它可以作為對phyloseq的補充，或者與之聯(lián)合使用。

使用MicrobiomeAnalystR包進行微生物組數(shù)據(jù)分析

MicrobiomeAnalystR是MicrobiomeAnalyst在線版本的R包（圖2/3，圖S4A-J， Pipeline 3. MicrobiomeAnalystR.Rmd）。包括多樣性分析（圖S4A-F）、差異分析（圖S4G）、生物標記識別（圖S4H-I）以及樣本測序庫大小概覽（圖S4J），這些功能比前兩個包更強大。在可視化方面，它結(jié)合了基礎(chǔ)包、ggplot繪圖以及交互式繪圖。在網(wǎng)絡(luò)分析方面，它提供了計算和繪制更適合微生物組數(shù)據(jù)的SparCC網(wǎng)絡(luò)的過程。然而，這個包依賴于許多來自CRAN、Bioconductor和GitHub的R包，因此完全安裝MicrobiomeAnalystR需要付出大量的努力。

使用Animalcules包進行微生物組數(shù)據(jù)分析

Animalcules包是在交互式平臺上進行分析的一種方式（圖2/3，圖S5A-J， Pipeline 4. Animalcules.Rmd）。它可以對樣本進行基礎(chǔ)信息統(tǒng)計與繪圖（圖S5A）或交互式餅圖（圖S5B），計算和可視化alpha多樣性（圖S5C），微生物組成或功能組成的熱圖和堆積柱狀圖（圖S5D-E），特征豐度箱線圖（圖S5F），屬水平的bray距離熱圖（圖S5G），排序分析（圖S5H-I），使用隨機森林、邏輯回歸選擇生物標記（圖S5J）等其他分析。這些分析的結(jié)果通?？梢酝ㄟ^交互式修改參數(shù)進行再分析，圖片可以進行交互式的放大和縮小，點擊查看詳細信息，以及由鼠標進行的其他操作以更好地修改。然而，結(jié)果不能導出為矢量格式，這并不滿足出版的要求。其次，分析內(nèi)容太少，特別是微生物組網(wǎng)絡(luò)分析，微生物組與其他指標之間的相關(guān)性分析。

使用microeco包進行微生物組數(shù)據(jù)分析

microeco包非常強大，使用R6類數(shù)據(jù)結(jié)構(gòu)（圖2/3，圖S6A-L， Pipeline 5. microeco.Rmd）。它包括微生物多樣性（圖S6A/B）物種組成（圖S6C-E），差異分析（圖S6F-H），生物標記物尋找（圖S6I-J），網(wǎng)絡(luò)分析（圖S6K），以環(huán)境因子探索（圖S6L），以及系統(tǒng)發(fā)育多樣性分析等。它幾乎可以完成目前所有的微生物組分析內(nèi)容。然而，它并不適合初學者使用，因為使用S6類對象有一定的門檻。另外，由于功能過多，對輸入數(shù)據(jù)的要求不同，導致一些功能難以使用。

使用amplicon包進行微生物組數(shù)據(jù)分析

amplicon包是微生物組分析工具包EasyMicrobiome（Liu等，2023）中的一個分析和繪圖工具（圖2/3，圖S7A-I，Pipeline 6. Amplicon.Rmd）。它能夠進行各種多樣性分析，包括alpha多樣性（圖S7A），稀化曲線（圖S7B），聚類距離熱圖（圖S7C）和PCoA（圖S7D），NMDS，LDA和PCA，物種組成（圖S7E/F），差異分析（圖S7G/H）。它可以輕松生成高質(zhì)量的圖像，如用于多樣性分析的箱線圖，散點圖，用于分類或功能組成的堆疊條形圖，環(huán)形圖和maptree圖。它的一個顯著特點是能夠精細調(diào)整圖像的呈現(xiàn)，從而生成可以發(fā)布的圖像。此外，amplicon包內(nèi)的幾個工具可用于微生物組數(shù)據(jù)的轉(zhuǎn)換，有助于使用LEfSe和STAMP等工具進行后續(xù)分析。然而，在當前版本，amplicon包并未提供網(wǎng)絡(luò)分析，微生物組-環(huán)境相互作用分析，以及群落形成過程分析的一些功能。作者在EasyAmplicon流程中提供了一些腳本來實現(xiàn)這些功能，并在發(fā)表的論文中提到計劃在未來完成這些功能。。

圖2. 微生物分析集成R包的功能介紹

微生物群落分析可分為多樣性分析、差異分析、生物標志物鑒定、相關(guān)和網(wǎng)絡(luò)分析、功能預測、以及其他微生物群落分析(群落建立/構(gòu)建過程，與其他指標的關(guān)聯(lián)分析）。

圖3. 微生物群落分析集成R包的典型結(jié)果和類似結(jié)果的比較

根據(jù)微生物群落分析功能的主要類別，對多個集成R包的分析結(jié)果進行分組。樹狀圖中的每個主要分支代表微生物群落分析的一種類型，共有10個主要分支：特征多樣性分析，包括（i）alpha多樣性分析；(ii)beta多樣性分析；(iii)特征（群落分類或功能）組成分析；(iv)進化或分類樹分析；（v）差異分析；（vi）生物標志物鑒定；（vii）相關(guān)和網(wǎng)絡(luò)分析；（viii）功能預測；（ix）群落形成/構(gòu)建過程分析；（x）其他分析，如與其他指標的關(guān)聯(lián)分析。每片葉子（圓圈）代表分析結(jié)果的一種風格樣式，葉子外面的圓圈號碼代表分析結(jié)果所來自的集成R包編號。

基于R語言的微生物組數(shù)據(jù)挖掘最佳流程

謝鵬昊，牛國慶翻譯，南京農(nóng)業(yè)大學，博士在讀；

過多的R包可能會阻礙微生物組研究人員有效地選擇適合微生物組相關(guān)分析的R包。因此，我們在六種分析中選擇了高效、常用、用戶友好的函數(shù)（圖S8）：1) 多樣性分析（圖S9A-I; 圖S10A-E），2) 差異分析（圖S10F-I; 圖S11A-B），3) 生物標志物識別（圖S11C-D），4) 相關(guān)性和網(wǎng)絡(luò)分析（圖S11E-I），5) 功能預測，6) 其他微生物組分析（圖S12A-I）。所有腳本都可以在Pipeline.BestPractice.Rmd文件中找到。

在這個流程中，我們使用amplicon包進行Alpha多樣性稀疏曲線（圖4A；圖S9A）和PCoA分析（圖4B；圖S9B），使用ggplot2包可視化微生物群落組成，使用ggClusterNet構(gòu)建Venn網(wǎng)絡(luò)（Chen等，2021年）（圖4C），使用ggtree和ggtrextre構(gòu)建進化樹（圖4D），使用LEfSe生成 cladograms（圖4E）。我們使用stst4、ggplot2和cowplot包進行差異分析并生成STAMP圖（圖4F），使用edgeR進行差異分析并在曼哈頓圖中進行可視化（圖4G），并應(yīng)用DESeq2進行差異分析并生成多組火山圖（圖4H）。我們還使用el071、caret、randomforest和ROC包進行各種機器學習分析并生成特征微生物重要性加權(quán)圖（圖4I）。此外，我們還使用ggClusterNet進行微生物網(wǎng)絡(luò)分析（圖4J），構(gòu)建網(wǎng)絡(luò)圖和繪制組合圖以探究環(huán)境因素和微生物群落之間的關(guān)系（圖4K）。最后，我們使用FEAST包進行群落溯源分析并構(gòu)建餅圖（圖4L）。其他分析包括微生物群落組成的堆積柱狀圖（圖S9E/H）、和弦圖（圖S10A）、Venn圖（圖S10C）、Upset圖（圖S10D）、差異分析火山圖（圖S10F）、功能預測等。

圖4. R語言中微生物群落分析的最佳實踐結(jié)果實例

所選結(jié)果包括稀疏曲線（A）、
主坐標分析散點圖（B），維恩網(wǎng)絡(luò)圖（C），進化樹（D），LEfSe圖（E），STAMP的差異分析誤差條圖（F），差異分析曼哈頓圖（G），差異分析多組火山圖（H），生物標志物選擇環(huán)列圖 (I），網(wǎng)絡(luò)圖（J），組合圖（K），溯源分析餅圖（L）

總結(jié)與展望

在過去的十年中，R語言和眾多R包在微生物組數(shù)據(jù)分析中發(fā)揮了重要作用。R語言易于使用和入門，吸引了許多研究人員學習使用。然而，在微生物組數(shù)據(jù)分析中，供給和需求之間仍存在一些矛盾。例如，在Windows系統(tǒng)下支持多線程往往很困難；其次，許多R包的運行速度相對較慢；第三，R包在微生物組應(yīng)用中仍需要進一步發(fā)展。例如，缺乏允許探索基于時間序列的微生物組成的軟件包，以及更強大的交互式軟件包來分析復雜的微生物數(shù)據(jù)。此外，ggplot2缺乏創(chuàng)建復雜和組合圖形的能力，無法滿足多個復雜指標與微生物群落數(shù)據(jù)之間關(guān)系的可視化要求。因此，有必要開發(fā)新的R包以更適合繪制復雜圖形和組合圖形。

隨著測序技術(shù)的發(fā)展，數(shù)據(jù)分析方法也隨著R包的發(fā)展而在微生物組領(lǐng)域取得了進步。這些R軟件包從經(jīng)典的R包如vegan，到集成的R包如phyloseq（包含了許多功能，并設(shè)置了統(tǒng)一的數(shù)據(jù)處理框架，已經(jīng)能夠?qū)崿F(xiàn)微生物組分析的大多數(shù)功能，包括微生物多樣性、差異、生物標志物鑒定、相關(guān)性和網(wǎng)絡(luò)、系統(tǒng)發(fā)育分析等）。然而，這些R包存在一些冗余的功能；例如，phyloseq、microbiome等可以進行微生物多樣性分析，差異僅在于可視化方法和方案。這種類似的情況在微生物組分析的R軟件包中一直存在，因此我們希望在未來的開發(fā)中，盡量避免重復使用相同部分或類似內(nèi)容，以突出R包的優(yōu)勢。

雖然這些R包可以實現(xiàn)很多功能，但在一些具體的分析中可能存在不足，比如在Alpha和Beta多樣性分析中，圖形往往沒有展示差異檢驗結(jié)果，難以直觀地從圖形中觀察出差異。此外，還有一些內(nèi)容需要進一步開發(fā)，例如將更多的機器學習方法應(yīng)用于微生物組數(shù)據(jù)，包括學習方法、模型和重要變量的評估等。其次，宏基因組的應(yīng)用越來越廣泛，基于Kraken（Wood和Salzberg，2014）、MEGAN（Huson等，2007）、MetaPhlAn2（Truong等，2015）、HUMAnN2（Franzosa等，2018）、eggNOG-mapper（Huerta-Cepas等，2017）等物種和功能注釋的結(jié)果從兆字節(jié)（M）升級到千兆字節(jié)（G）。因此，應(yīng)使用更快速的數(shù)據(jù)處理R包來進行微生物組數(shù)據(jù)分析，例如data.table、fst、tidyfst等。

使用適當?shù)臄?shù)據(jù)結(jié)構(gòu)可以加速微生物組數(shù)據(jù)分析。起初，我們使用S4類對象進行微生物組數(shù)據(jù)封裝，可以全面有效地完成各種分析。R6類對象和其他對象的出現(xiàn)極大地影響了微生物組數(shù)據(jù)的處理，并在很大程度上促進了數(shù)據(jù)的處理。隨著R語言tidy家族的發(fā)展，最近出現(xiàn)了基于tidy的數(shù)據(jù)結(jié)構(gòu)，用于微生物組數(shù)據(jù)挖掘。例如，MicrobiotaProcess包（Xu等，2023）。這種結(jié)構(gòu)更適用于微生物組數(shù)據(jù)挖掘、機器學習建模等分析，在分析中可以更容易提取實驗設(shè)計、時間、空間等因素對微生物組數(shù)據(jù)的影響，發(fā)現(xiàn)深層次的規(guī)律。我們期待R語言能使微生物組分析更加高效，幫助大家發(fā)現(xiàn)更多關(guān)于微生物在人類、動物、植物和環(huán)境中的作用，并利用微生物為我們造福，使世界變得更加美好。

翻譯：楊盛蝶/牛國慶/謝鵬昊，南京農(nóng)業(yè)大學，博士在讀
校對：文濤，南京農(nóng)業(yè)大學，博士后
審核：劉永鑫，中國農(nóng)業(yè)科學院深圳農(nóng)業(yè)基因組研究所，研究員

引文格式

Tao Wen, Guoqing Niu, Tong Chen, Qirong Shen, Jun Yuan, Yong-Xin Liu. 2023. The best practice for microbiome analysis using R. Protein & Cell pwad024. https:///10.1093/procel/pwad024