近年的研究熱點集中于環(huán)境和生物體相互作用的微生物群體，而大量復(fù)雜的微生物群體存在培養(yǎng)困難，構(gòu)成復(fù)雜（包括細(xì)菌,、古菌,、真菌、原生生物,、病毒甚至小型真核生物）,。因此如何用高通量精準(zhǔn)的了解這些群體的構(gòu)成，基因功能分布以及具體的表達(dá)活性和代謝狀況成為首要問題,。

高通量測序技術(shù)的發(fā)展,，讓我們可以不經(jīng)過培養(yǎng)，一次性了解微生物群落構(gòu)成甚至基因代謝組成,。

隨著技術(shù)的進(jìn)步,，檢測方法也逐漸豐富，對應(yīng)的分析手段和軟件算法也逐步完善,，使我們可以根據(jù)研究需要選擇不同的檢測和分析策略來獲得海量的數(shù)據(jù)并進(jìn)行相應(yīng)的研究分析,。

01

簡 介

免于培養(yǎng)的微生物學(xué)研究方法主要基于測序，高通量測序使我們一次可以獲得整個微生物群體的數(shù)據(jù)信息,，簡單來說包括兩種策略：

1,、基于特定標(biāo)記基因的擴(kuò)增測序方案（常見的16s，ITs,，18s或特定功能基因）

2,、對整個群落DNA進(jìn)行測序，獲取全部微生物基因組進(jìn)而進(jìn)行分類和功能分析的策略（鳥槍法宏基因組測序shotgun metagenomics）,。

基于16s基因的分析方法

由于其極低的成本,，對于樣本DNA的低要求非常適合于大規(guī)模群體樣本的調(diào)查和分析，隨著DADA2等分析方法的改進(jìn),，物種分類精度和準(zhǔn)確度也有所提升,，加上PICRUST等功能預(yù)測方法一定程度上彌補(bǔ)了基因信息的缺失，因此16s這類基于基因的微生物研究方法仍然是不可或缺的方案,。

下表列了16s常見的分析軟件,，目前QIIME2作為整合包使用最為方便，VSEARCH也作為UPARSE的開源版本使用也非常廣泛,。

16s測序的分析流程如下圖,，獲得序列經(jīng)過聚類后獲得OTU或ASV，并得到相對豐度,。

經(jīng)過PICRUSt可以得到預(yù)測的基因分類豐度,，進(jìn)而進(jìn)行alpha多樣性和Beta多樣性以及組間差異和相關(guān)性分析。

PICRSt的工作原理如下圖,，將OTU表內(nèi)16s序列進(jìn)行對應(yīng)物種16s拷貝數(shù)標(biāo)準(zhǔn)化后,，將物種豐度乘以已經(jīng)整理好的物種的基因注釋數(shù)表就獲得基因的預(yù)測豐度,。

02

淺宏基因組

淺宏基因組測序方案是去年knights-lab在msystems上發(fā)表的針對16s分辨率和宏基因組高成本之間的一個折中方案，通過降低測序深度,，每個樣本50萬reads,，但是物種的分辨率并沒有低于一般宏基因組（普遍5~10G數(shù)據(jù)量）。

不通過拼接組裝,，直接基于kraken2等kmer,，或MetaPhlAn2等標(biāo)記基因的參考基因組方法進(jìn)行種屬豐度分類。結(jié)合其到菌株的物種分類和豐度數(shù)據(jù)可較16s方案下的PICRUST更加準(zhǔn)確的預(yù)測基因構(gòu)成,。

Hillmann B, Al-Ghalith GA, Shields-Cutler RR, Zhu Q, Gohl DM, Beckman KB, Knight R, Knights D. 2018. Evaluating the information content of shallow shotgun metagenomics. mSystems 3:e00069-18. https:///10.1128/mSystems.00069-18.

我們發(fā)現(xiàn)有些小伙伴的需求是：

想要獲得更全和更精細(xì)分類精度同時不需要獲得完整基因組序列和重建菌群基因的,。

那么這時候，我們提供的淺宏基因組測序就可以成為很好的選擇,，其成本低（快要接近16s測序分析的價格了,，文末有福利），分析簡便快速,，同樣能獲得宏基因組的基本豐度數(shù)據(jù)。不過淺宏基因組也有其適用范圍,，根據(jù)樣品類型的不同,，一些樣品可能包含 >99％的人類宿主DNA，這不僅增加了序列成本,，而且給測量帶來了不確定性,。

在許多研究中也會采取在進(jìn)行宏基因組測序文庫的準(zhǔn)備之前去除宿主DNA的方法。但是,，在去除宿主DNA后,，可能沒有足夠的微生物基因組DNA用于宏基因組測序，這通常需要最少50ng的輸入,。因此淺宏基因組較適合于宿主DNA含量較低的樣本,，如人類糞便、水體,、土壤等,；而如口腔唾液、肺泡灌洗液,、血液等人體體液類樣本就不太適合,。

下圖是宏基因組測序數(shù)據(jù)中比對到人類基因組的序列比例，根據(jù)樣本類型不同而不同,。

我們可以免費提供針對糞便及環(huán)境樣本助力臨床/科研取樣,。

人體口腔、痰液,、腹水,、腦脊液,、尿液、皮膚,、陰道分泌物等高寄主細(xì)胞含量樣本可根據(jù)我們的處理方案簡單處理后大幅降低宿主DNA比例,。

處理方案如下：

可向上滑動

本處理方案以后宿主DNA可以降低8%以下。

03

宏基因組

說起宏基因組,，對于熟悉宏基因組或者打算做宏基因組的同學(xué)可能已經(jīng)迫不及待想知道這個怎么分析啊,，怎么看結(jié)果啊之類的問題... 但在這之前，首先你應(yīng)該了解的是宏基因組是什么,，做宏基因組你能得到什么,。

此外，對于缺乏深度研究和高質(zhì)量參考基因組的樣本,，如土壤和特殊環(huán)境下的樣本,，宏基因組獲得的較為完整的基因組不僅可以豐富參考基因組數(shù)據(jù)庫，同時可以提供更加準(zhǔn)確的物種分類,。

至此,，我們了解了16s,、淺宏基因組、宏基因組三種方式,，我們將它們各自的特點總結(jié)如下表,，便于你更直觀地去了解（文末有福利~）。

宏基因組報告中有哪些分析內(nèi)容,？

測序數(shù)據(jù)需要經(jīng)過質(zhì)檢,，去除接頭和低質(zhì)量序列，一般還會進(jìn)行一步過濾人的基因組序列,，然后分為兩個路徑,，使用參考數(shù)據(jù)的比對方法和從頭組裝的方法，下圖是一個完整的宏基因組分析流程：

看完上圖,，可以對宏基因組測序的基本流程有個大致了解,。

對于宏基因組測序而言，最重要的就是獲得微生物群準(zhǔn)確的物種構(gòu)成及其豐度,。

首先你需要了解的是無論16S測序還是宏基因組測序獲得的均是相對豐度,，即每種菌占所有菌屬的比例。

要獲得絕對的豐度需要在取樣時做好取樣量的計量,，并在提取和建庫中加入已知絕對量的參照DNA,。

宏基因組測序獲得物種構(gòu)成及其豐度有以下兩條路可以走：

我們先講其中之一：直接比對 。

直接比對是基于參考數(shù)據(jù)的,，那么基于參考數(shù)據(jù)的物種構(gòu)成分析主要有兩類方法：

一類是基于Kmer和LCA比對特征來分析對應(yīng)物種豐度,，如kraken2等。

另一類是基于特征標(biāo)記基因進(jìn)行分析的,，如MetaPhlAn2等,。

基于參考基因組的分析工具如下表：

除了上面表中列出來的,，另外還有

Centrifuge：比kraken2慢2x,內(nèi)存使用少很多

Sourmash：類似CLARK,可以使用整個refseq作為數(shù)據(jù)庫,。

主流的kraken2——快速、準(zhǔn)確度高,、內(nèi)存要求高

目前主要使用kraken2為主,，因為快速，準(zhǔn)確度也相當(dāng)不錯,。不過,，對于內(nèi)存的要求較高，另外受數(shù)據(jù)庫本身質(zhì)量影響較大,，默認(rèn)kraken2的參考數(shù)據(jù)庫只包括了細(xì)菌,、古菌、病毒和人,，還需要添加其他域的參考基因組,。但涵蓋的測序參考種仍然有限，對于菌株水平的鑒定受一定影響,。后續(xù)使用Bracken可以針對kraken2的比對結(jié)果進(jìn)行計算相對豐度,。

MetaPhlAn2——物種跨度大,、實用

MetaPhlAn2首先從全基因組數(shù)據(jù)庫中找出clade-specific marker genes，然后利用這個marker genes的數(shù)據(jù)庫對高通量測序得到的shotgun序列進(jìn)行注釋,，目前主要用于后面直接使用reads獲得基因和代謝通路豐度的HUMANn2的流程中,，其物種跨度較大，速度也可以接受,。

以上我們了解了直接使用reads獲得豐度,。

如果有足夠測序深度和樣本數(shù)量還可以通過組裝出參考基因組來鑒定獲得。該部分我們在下面的組裝和分箱流程部分詳細(xì)講,。

使用Kraken2對其中的微生物進(jìn)行物種注釋。我們的Kraken2使用的數(shù)據(jù)庫是由Refseq（2020.04.20）細(xì)菌,，古細(xì)菌,、真菌、原生動物和病毒庫以及GRCh38人類基因組構(gòu)建的,。

通過查詢數(shù)據(jù)庫序列中的每個k-mer,，然后使用所得的LCA分類單元集確定序列的適當(dāng)標(biāo)簽，對序列進(jìn)行分類,。數(shù)據(jù)庫中沒有k-mers的序列不會被Kraken2分類,。這里我們是在使用k-mer=35的條件下進(jìn)行物種注釋。

使用Bracken對物種注釋結(jié)果計算相對豐度,。Bracken是一種高度精確的統(tǒng)計方法,，可從宏基因組學(xué)樣本計算DNA序列中物種的豐度。Braken使用Kraken2分配的分類標(biāo)簽來估計源自樣本中每種物種的讀數(shù)數(shù)量,。

對物種注釋結(jié)果使用 KRONA 進(jìn)行可視化展示,。

注：圓圈從內(nèi)到外依次代表不同的分類級別（界門綱目科屬種），扇形的大小代表不同注釋結(jié)果的相對比例,。

上面的是使用KRONA對單個樣本的構(gòu)成圖形化,，所有樣本合并使用柱狀圖就可以了解具體的樣本構(gòu)成豐度，從門-綱-目-科-屬-種-甚至菌株每個層次都可以進(jìn)行顯示（下面是截取我們報告中的相關(guān)圖）,。

如果嫌柱狀圖的展示方式單一,，當(dāng)然也可以有別的選擇。比如說以Circos的環(huán)圖形式展現(xiàn)：

也可以進(jìn)行聚類分析：

有了這些數(shù)據(jù)我們就可以進(jìn)行alpha多樣性（指每個樣本內(nèi)部菌群多樣性）的分析了,。

各樣本和多組之間也可以進(jìn)行Beta多樣性的比較分析：

計算樣本之間的菌屬構(gòu)成相似度：

組間的差異分析：尋找差異或代表性菌屬,，如下：

Trukey多組間檢驗

LefSe分析

其中LEfSe基于線性判別分析（Linear discriminant analysis，LDA）的分析方法,，篩選組與組之間生物標(biāo)記物Biomarker（基因,、通路和分類單元等），即組間差異顯著物種或基因,。當(dāng)分組較多時較難獲得每個分組獨特的Biomarker,。

以上是關(guān)于物種組成部分,，但是有些小伙伴會有這樣一些疑惑：物種構(gòu)成變化很大怎么辦？個體差異也很大,？之類的諸多疑問,。

下圖可以幫你更好地理解上面這段話,。從圖中我們可以看到,，舌背樣本和糞便樣本雖然在種屬上有很大差異，但它們在基因功能層面卻有著相似的構(gòu)成,。

與物種構(gòu)成豐度的分析類似,，基因功能構(gòu)成分析也同樣可以包括兩種方法：

方法一、通過直接基于reads的參考數(shù)據(jù)庫方法獲得

方法二,、通過組裝后預(yù)測注釋基因并得到豐度

在具體展開方法之前,，我們需要先了解關(guān)于基因功能的基本概念。

我們先來看方法一,，具體是如何操作的,？

主流的HUMAnN2——獲得基因和代謝通路豐度的同時可直接進(jìn)行下游分析

基于測序原始序列直接獲得基因構(gòu)成豐度的軟件目前最主要的是HUMAnN2,，其首先使用MetaPhlAn2進(jìn)行物種分類（關(guān)于這個軟件我們在前面物種組成部分已經(jīng)講過），并提取相應(yīng)物種參考基因組用于比對,，未比對上的用于進(jìn)一步和uniref數(shù)據(jù)庫進(jìn)行蛋白質(zhì)序列比對,。原理見下圖：

下圖是基于HUMAnN2的不同代謝通路的菌貢獻(xiàn)比例圖：

在我們的宏基因組報告中獲得的是這樣的：

而另外一種方法是通過組裝獲得，我們在前面物種構(gòu)成小節(jié)也已經(jīng)提到過組裝分析,，那么這里我們就組裝拼接分析這部分展開講解一下,。

什么樣的條件下可以進(jìn)行組裝分析？

當(dāng)測序深度足夠的情況下,，目前illumina二代和Pacbio以及Nanopore等長片段測序技術(shù)已經(jīng)足以組裝出高質(zhì)量的細(xì)菌基因組草圖,，結(jié)合Binning方法可以一次性獲得大量物種的高質(zhì)量接近完成基因組。此外還有Hi-C等手段可以進(jìn)一步完成基因組以及對應(yīng)質(zhì)粒的完整拼接,。

組裝的流程是什么樣的,？

來看一下整個基于組裝的流程：

① 提取、測序

首先從樣本中提取基因組DNA,，進(jìn)行測序,，可以使用Illumina的段片段深度測序也可以輔助三代長片段測序。

② 獲得contig序列

接著對序列經(jīng)過質(zhì)檢過濾處理后直接使用序列進(jìn)行拼接,，獲得contig序列,，這時通常每個菌的基因組會有幾十到數(shù)千個contig片段，由于構(gòu)成復(fù)雜,，很多近緣菌之間的基因組存在大量相似序列,，以及每種菌豐度都不一致，所以contig階段的片段仍然較多,。

③ Binning分析 

基于序列構(gòu)成特征如GC含量,、核苷酸多態(tài)性、覆蓋度以及基因的物種相似度等多種數(shù)據(jù),，如果有多個樣本或梯度可以同時結(jié)合樣本豐度變化來進(jìn)行分箱也就是Binning分析,，將具有相同特征和變化的contig聚類歸為同一個來源的箱，每個bins通常來自單一菌也就是一個菌株的基因組（我們的數(shù)據(jù)分析中包含這部分分析內(nèi)容）,。

④ 進(jìn)一步質(zhì)檢評估

之后會進(jìn)行進(jìn)一步的質(zhì)檢,，如checkM等評估每個Bin的完整度（核心基因以及rRNA等的完整性）和污染比例（如錯誤拼接，不同物種來源等）,。一般要求50%以上的完整度以及10%以下的污染,，當(dāng)然樣本數(shù)量越多，測序深度越高，測序讀長越長理論上binning的質(zhì)量也會更好,，能獲得更多高質(zhì)量的單一菌完整基因組,。

借用一張分箱的說明PPT：

⑤ 注釋

最后使用PROKKA進(jìn)行基因預(yù)測,，獲得的編碼序列我們經(jīng)過進(jìn)一步CD-Hit聚類去冗余，然后使用eggNOG-mapper對其進(jìn)行進(jìn)一步的功能比對注釋,。使用salmon完成基因的定量,，這樣我們就得到基于組裝注釋的基因豐度數(shù)據(jù)了。之后就可以進(jìn)行基因和功能層面的多樣性,、構(gòu)成以及樣本和組間差異分析,。

我們獲得的最基礎(chǔ)的uniref，eggNOG,，KEGG和GO等注釋如下：

KEGG

COG

eggNOG

組間差異分析,，如KEGG途徑：

除此之外，還可以使用其他的功能基因數(shù)據(jù)庫來進(jìn)行進(jìn)一步的基因注釋和分析,。比如：

CAZy：

VFDB毒力因子注釋：

抗性基因注釋：

TCDB數(shù)據(jù)庫注釋：

PHI數(shù)據(jù)庫注釋：

BCGs分析：

以及基于antiSMASH和BiG-SCAPE來對代謝物的合成生物基因簇BCGs進(jìn)行分析,。

固定代謝能力評估：

或更聚焦于特定代謝的如下圖中的氮、磷,、硫和碳固定代謝能力和水平的評估：

當(dāng)有了大量樣本的菌群構(gòu)成豐度信息,，以及各種基因和代謝豐度數(shù)據(jù)后，我們需要根據(jù)樣本的meta信息,，基于不同分組,，時間或環(huán)境因子等數(shù)據(jù)進(jìn)行統(tǒng)計分析和檢驗，進(jìn)而發(fā)現(xiàn)和探索可能的關(guān)聯(lián)以及背后的生物學(xué)意義,。

那么在面對宏基因組這類數(shù)據(jù)時在進(jìn)行統(tǒng)計檢驗分析時需要注意什么呢,，應(yīng)該采用哪些分析，并如何解讀這些結(jié)果呢,？

首先,，微生物組數(shù)據(jù)分析分為四大類：

在對所有數(shù)據(jù)進(jìn)行統(tǒng)計檢驗前一般建議對數(shù)據(jù)進(jìn)行基本的質(zhì)量過濾,。一類是去除絕大部分樣本都不存在的物種和基因,，如Prevalence in samples (20%)，還有一類是去除變異度過小的Percentage to remove (10%)基于Inter-quantile range。

為什么可以過濾這兩類,？

宏基因組數(shù)據(jù)具有一些獨特的特征,，例如測序深度的巨大差異，稀疏性（包含許多零）和分布的巨大差異（過度分散）,。在進(jìn)行后續(xù)的統(tǒng)計檢驗之前建議針對不同的分析方法進(jìn)行相應(yīng)匹配的標(biāo)準(zhǔn)化處理,。標(biāo)準(zhǔn)化包括：

Rarefaction和縮放方法：這些方法通過將樣本放到相同的比例進(jìn)行比較來處理不均勻的測序深度；

轉(zhuǎn)換方法：包括處理稀疏性,，組成性和數(shù)據(jù)中較大變化的方法,。

那么各種標(biāo)準(zhǔn)化方法是什么，應(yīng)該選擇哪種方法,？

參考MicrobiomeAnalyst網(wǎng)站提供的信息,，以下是一個簡短的介紹：

請注意，數(shù)據(jù)標(biāo)準(zhǔn)化主要用于可視數(shù)據(jù)探索,，例如beta多樣性和聚類分析,。有時候不使用標(biāo)準(zhǔn)化也能獲得最佳結(jié)果，比如：單變量統(tǒng)計和LEfSe,。

同時,，其他比較分析將使用其自己的特定標(biāo)準(zhǔn)化方法。例如,，對metagenomeSeq使用累積總和定標(biāo)（CSS）標(biāo)準(zhǔn)化,，對edgeR應(yīng)用M值的修剪均值（TMM）。

經(jīng)常有小伙伴問,，這個數(shù)據(jù)是用的什么標(biāo)準(zhǔn)化,？沒有做標(biāo)準(zhǔn)化怎么辦？這類問題,。

目前,，尚無關(guān)于應(yīng)使用標(biāo)準(zhǔn)化的共識性指南。建議大家可以探索不同的方法,，然后目視檢查分離模式（即PCoA圖）以評估不同標(biāo)準(zhǔn)化程序?qū)嶒灄l件或其他感興趣的宏基因組數(shù)據(jù)的影響,。

有關(guān)這些方法的詳細(xì)討論，請參考使用者最近發(fā)表的兩篇論文 

①  Paul J. McMurdie等

（https:///10.1371/journal.pcbi.1003531） 

② Jonathan Thorsen等

（ http:///10.1186/s40168-016-0208-8  ）

以上是關(guān)于標(biāo)準(zhǔn)化的這部分內(nèi)容需要了解的知識,，接下來我們來看具體如何操作,，怎么得到那些圖表？它們分別代表什么,？

一般我們需要先進(jìn)行探索性分析,，也就是不設(shè)預(yù)訂的假設(shè)，首先從主成分分析結(jié)果中了解樣本的菌屬和基因的大概分布。

主成分分析是根據(jù)不同距離算法計算樣本之間的距離矩陣,，然后進(jìn)行降維,，最終形成一個三維的空間分布。樣本之間在空間上分隔越遠(yuǎn)表明樣本之間的差異越大,。

比如我們報告中的下圖,，疾病和正常樣本可以較好的區(qū)分，一般此處我們還會進(jìn)行一個統(tǒng)計檢驗,，來判別PC1和PC2這幾個維度上兩組之間是否真的存在統(tǒng)計差異,。

基于豐度圖來評估各樣本和分組的基本構(gòu)成，如：

之后我們可以針對不同分組或處理之間的樣本進(jìn)行統(tǒng)計檢驗,，可以使用的檢驗方法包括兩組間的非參數(shù)統(tǒng)計檢驗T-test/ANOVA,，3組以上組間統(tǒng)計檢驗可以使用Tukey?test，其直接生成各組將的統(tǒng)計差異,，并提供字母標(biāo)注,，直觀簡便，如：

具體的統(tǒng)計方法選擇可以參考下表：

除了常規(guī)的非參數(shù)檢驗外,，包括metagenomeSeq和DEseq以及edgeR等統(tǒng)計方法包可以很好的分析組間差異特征,。LEfSe則一般用于尋找特征標(biāo)志物。

那么有了大量的差異特征菌屬或基因之后,，我們是否能基于這些差異菌屬有效的區(qū)分不同的分組呢,，或構(gòu)件一個模型來預(yù)測或分類呢？

這時候可以使用隨機(jī)森林（Random Forest）一類的決策樹機(jī)器學(xué)習(xí)模型,，來利用這些差異特征構(gòu)建分類模型,，并使用AUC等指標(biāo)來評估基于這些模型的預(yù)測有效性和準(zhǔn)確度（我們報告中如下圖）。

當(dāng)然也可以使用其他更復(fù)雜的如深度學(xué)習(xí)等方法來構(gòu)建分類模型,。

除了性別,、疾病、地點等分類差異之外,，我們通常還有很多元數(shù)據(jù),，包括臨床指標(biāo)或環(huán)境因子等信息，這些數(shù)據(jù)通常是連續(xù)型數(shù)值,，對于這類數(shù)據(jù)我們可以進(jìn)行相關(guān)性分析,。

當(dāng)然反過來，將菌群特征作為表型也可以和如基因組的基因型或SNP構(gòu)成來進(jìn)行相關(guān)性分析,。

對于菌群數(shù)據(jù)的相關(guān)性分析比較推薦：

SparCC方法,，可以構(gòu)建菌種或菌屬之間的相關(guān)性網(wǎng)絡(luò)，相對穩(wěn)定,。

對于與疾病或環(huán)境變量進(jìn)行相關(guān)性分析可以使用：

Sperman秩相關(guān)分析,。

另外RDA/CCA分析也可以有效的反映菌屬與環(huán)境因子等指標(biāo)直接的關(guān)系（我們報告中如下圖）,。

Mentel檢驗也可以用于判斷菌群構(gòu)成特征與單個或一組環(huán)境因子之間是否存在顯著相關(guān)。

以上看完后，你應(yīng)該對宏基因組的數(shù)據(jù)分析流程有了整體的認(rèn)識,，也學(xué)會了相應(yīng)的一些操作,，但是不一定能直接從自己的這些數(shù)據(jù)、圖表中真正探索到和實際生物學(xué)相關(guān)的有價值的研究成果,。

所以,，我們又選取了一些已發(fā)表的研究作為案例，結(jié)合實際問題來具體分析,，從實驗設(shè)計到具體分析流程方法和圖表的展示,，再到相應(yīng)的結(jié)論，掌握這類文章的總體思路,。

之后無論是剛開始的實驗設(shè)計,，還是后面的分析，都會更加得心應(yīng)手,。

建議想好整個實驗思路再開始（或者也可以咨詢我們,，我們專業(yè)的數(shù)據(jù)分析團(tuán)隊會為你提供切實可行的項目方案）。

04

案例解析

  < 案例一 > 肥胖患者的腸道微生物組

第一項研究是關(guān)于肥胖患者減肥手術(shù)后的宏基因組和代謝數(shù)據(jù)的分析研究,。

文獻(xiàn)來源：Aron-Wisnewsky J, Prifti E, Belda E, et al. Major microbiota dysbiosis in severe obesity: fate after bariatric surgery. Gut. 2019;68(1):70–82. doi:10.1136/gutjnl-2018-316103

研究納入了61名嚴(yán)重肥胖的受試者,，他們是可調(diào)節(jié)胃束帶術(shù)（AGB，n = 20）或Roux-en-Y胃旁路術(shù)（RYGB,，n = 41）的候選人,。減肥手術(shù)后1、3和12個月隨訪24名受試者,。使用宏基因組學(xué)測序和LC-MS分析腸道菌群和血清代謝組,。另外納入了10人和147人分別作為宏基因組和代謝檢測的驗證集。

研究思路

這樣的設(shè)計分別有什么作用,？

從下圖可以看到研究針對樣本的總基因多樣性水平與生理指標(biāo)和疾病狀態(tài)進(jìn)行相關(guān)性分析和組間差異分析,，圖中給出了顯著相關(guān)和差異的指標(biāo)。

使用的統(tǒng)計檢驗方法是pearson和sperman相關(guān)和t-test以及Kruskal-Wallis檢驗,。

下圖是研究將MAGs與各項生理和代謝值進(jìn)行相關(guān)性分析后的熱力圖,。該研究由于測序較早，并未獨立拼接,，而是直接使用了之前一項人類腸道菌群研究獲得組裝基因組參考序列,。

進(jìn)一步研究分析了術(shù)后特定變化模式的MAGs以及它們與代謝生理指標(biāo)的相關(guān)性，見下圖：

上圖的研究可以通過pattern search的方法尋找特定變化模式的菌種,。

研究的主要結(jié)論發(fā)現(xiàn)是低基因豐富度（LGC）存在于75％的患者中,，并且與軀干脂肪質(zhì)量和合并癥（2型糖尿病，高血壓和嚴(yán)重程度）增加相關(guān),。LGC改變了78種宏基因組種（MGS）,，其中50％與不良的身體成分和代謝表型有關(guān)。九種血清代謝產(chǎn)物（包括谷氨酸鹽,，3-甲氧基苯基乙酸和L-組氨酸）和含有參與其代謝的蛋白質(zhì)家族的功能模塊與低MGR密切相關(guān),。術(shù)后一年，BS會增加MGR,，但盡管RYGB患者的代謝改善比AGB患者大,，但術(shù)后一年的MGR仍然很低。

  < 案例二 >  食物與人類腸道微生物組

第二項研究是Dan Knights實驗室發(fā)表在Cell Host & Microbe,2019的一篇針對34個人17天每日飲食和菌群變化的相關(guān)研究,，試圖揭示日常食物選擇與人類腸道微生物組組成之間的精細(xì)關(guān)系,。

文獻(xiàn)來源：Johnson Abigail J,Vangay Pajau,Al-Ghalith Gabriel A et al. Daily Sampling Reveals Personalized Diet-Microbiome Associations in Humans.[J] .Cell Host Microbe, 2019, 25: 789-802.e5.

可以看到，研究同時記錄了糞便樣本的菌群宏基因組和每日的飲食記錄,。研究的核心在于將每日飲食的食物通過營養(yǎng)構(gòu)成進(jìn)行量化,，并構(gòu)建類似物種進(jìn)化樹的食物物候樹,。

此外由于有每日的數(shù)據(jù)，可以通過前一日的食物與第二日的菌群數(shù)據(jù)進(jìn)行時間序列分析,，構(gòu)建食物與菌之間的關(guān)聯(lián)以及時間相關(guān)性,。

最后基于菌群數(shù)據(jù)和前一日飲食來構(gòu)建模型預(yù)測判斷后一日的菌群狀態(tài)，幫助我們了解食物對于個體菌群的影響因素并實現(xiàn)定量和預(yù)測,。

研究中對數(shù)據(jù)的處理過濾標(biāo)準(zhǔn)如下：刪除所有具有低讀取計數(shù)（每個樣品<23,500個讀?。┑臉悠贰Ｎ锓N級別的分類表僅限于研究對象中至少存在25％的研究日,，且在10%的研究樣本對象中發(fā)現(xiàn)的那些物種,。

最后,，相對豐度<0.01％的稀有物種被丟棄,，將物種數(shù)量限制為290個注釋。將得到的分類表匯總到較高的分類級別（即屬,，科,，門等），以進(jìn)行下游分析,。

菌群和飲食以及營養(yǎng)構(gòu)成的堆疊圖很好展現(xiàn)了變化和對應(yīng),。

下面這張圖很好的顯示了飲食食物的變化與菌群變化之間的時間變化關(guān)系：

下面這張圖通過對每個人單獨進(jìn)行菌屬與食物的Spearman相關(guān)，展現(xiàn)了菌與食物之間的關(guān)聯(lián)的個體化差異,，在特定菌屬對應(yīng)相同食物不同人會出現(xiàn)完全不同方向的變化,，這也正是這項研究所揭示的，這種關(guān)聯(lián)關(guān)系的復(fù)雜性,。

  < 案例三 >  類風(fēng)濕關(guān)節(jié)炎的人群腸道微生物組 

接下來的一個研究是比較典型的宏基因組組裝并與疾病進(jìn)行關(guān)聯(lián)分析的案例,，研究的是日本人群類風(fēng)濕關(guān)節(jié)炎的腸道微生物組的全基因組關(guān)聯(lián)研究。

文獻(xiàn)來源：Kishikawa Toshihiro, Maeda Yuichi,Nii Takuro et al. Metagenome-wide association study of gut microbiome revealed novel aetiology of rheumatoid arthritis in the Japanese population.[J] .Ann. Rheum. Dis., 2020, 79: 103-111.

研究使用較高深度的宏基因組shotgun測序（每個樣本平均13 Gb）對日本人群（病例 = 82,，對照 = 42）進(jìn)行了RA腸道微生物組的MWAS分析,。MWAS由三個主要的生物信息學(xué)分析渠道（系統(tǒng)發(fā)育分析、功能基因分析,、途徑分析）組成,。

使用了之前研究中6139個完成拼接日本人腸道宏基因組作為參考序列以及其他幾項研究的參考基因組，在過濾部分種過多的基因組之后,，最后一共使用了7881個參考基因組,。

將QC后的序列直接比對到參考基因組,，并根據(jù)基因組長度計算對應(yīng)物種的相對豐度。

基因方面選擇denovo組裝,，使用MegaHIT,，然后再contig上完成ORF預(yù)測和CD-HIT聚類去冗余，最后與UniRef和KEGG數(shù)據(jù)庫比對,。

最后使用bowtie2將測序序列比對到注釋后的unigene序列上獲得基因豐度,，經(jīng)過KEGG注釋得到代謝途徑的豐度。研究的數(shù)據(jù)處理流程圖如下：

在數(shù)據(jù)分析流程和方案選擇上人體腸道菌群由于研究眾多,，以及有多個深度測序拼接完成的Binning參考基因組數(shù)據(jù)集,，確實可以直接使用參考基因組直接比對。

對于其他一些環(huán)境或來源的樣本這個深度的數(shù)據(jù)量可以考慮獨立拼接和分箱,。該研究中使用已有參考基因組,，大概88%的序列能比對到參考基因組，如果直接組裝這個比例應(yīng)該可以更高一些,。另外在獲得基因豐度是可以考慮使用Salmon,，比對獲得基因豐度更為方便。

獲得相應(yīng)數(shù)據(jù)后對相對豐度,，該研究使用Box-Cox transformation對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化,，并過濾了一些低豐度的菌屬。Case-control的相關(guān)性分析使用的R的glm2模塊,，將年齡,、性別和測序上機(jī)分組作為協(xié)變量。

對于菌屬的關(guān)聯(lián)分析,，最終將顯著性結(jié)果以火山圖和GraPhlAn圖的形式展現(xiàn)如下：

上面其中D圖使用每個菌的豐度進(jìn)行UMAP分析,，并映射關(guān)聯(lián)效應(yīng)的展示比較有意思。

不過在基因?qū)用?/strong>上并未找到相應(yīng)的關(guān)聯(lián),，可以看到下圖UniRef90的NMDS分布圖兩組之間無法有效區(qū)分,，多樣性也沒有顯著差異。

  < 案例四 ＞ 永凍土中參與有機(jī)物降級的關(guān)鍵菌群

最后這項研究是對來自永凍土融化梯度的214個樣品的宏基因組測序組裝了1,529個基因組，揭示了參與有機(jī)物降解的關(guān)鍵種群，包括其基因組編碼先前未描述的木糖降解真菌途徑的細(xì)菌,。

文獻(xiàn)來源：Woodcroft Ben J,Singleton Caitlin M,Boyd Joel A et al. Genome-centric view of carbon processing in thawing permafrost.[J] .Nature, 2018, 560: 49-54.

通過宏基因組denovo組裝和分箱Binning,，最終獲得了1529個永凍土菌群基因組?；谶@些數(shù)據(jù)描繪了永凍土融化梯度下的菌群構(gòu)成特征,，如下圖。

論文是2018年發(fā)表的,，測序是在2011和12年測的,，使用的是CLC Genomics Workbench 較早的4.4版分單樣本組裝，然后使用MetaBAT進(jìn)行分箱,，最后的標(biāo)準(zhǔn)是70%完成度和低于10%的污染,。

其中糖苷水解酶基因使用dbCAN數(shù)據(jù)庫的HMM進(jìn)行預(yù)測，碳代謝使用KEGG數(shù)據(jù),。

研究還同時選擇了部分樣本進(jìn)行了宏轉(zhuǎn)錄組和宏蛋白組的測序,，對碳代謝同時結(jié)合轉(zhuǎn)錄組和蛋白組的數(shù)據(jù)，展現(xiàn)了特定通路下不同永凍土的菌群構(gòu)成和表達(dá)豐度差異,。

基因組拼接的分布情況,，以及不同地域樣本分布和菌屬豐度情況如下：

木糖降解途徑在每個樣本中的分布和維恩圖,，另外詳細(xì)的展現(xiàn)了主要門對每個代謝途徑的貢獻(xiàn)和基因表達(dá)豐度,，如下圖：

這張圖分析了特定菌與地理位置和CO2以及甲烷的濃度的關(guān)聯(lián)性，如下圖：

對關(guān)鍵物種的CH 4 :CO2濃度比相關(guān)代謝途徑的重建,，以及相應(yīng)基因的基因家族分析,。

p.s. 以上展示的圖表,，我們都可以幫你實現(xiàn)~ 

05

工具分享

 一、 MicrobiomeAnalyst 

網(wǎng)址：https://www./,，只需要biom文件或豐度表就可以進(jìn)行絕大部分統(tǒng)計檢驗分析,，而且生成圖表完整,，可以直接使用。偶爾會有服務(wù)器不穩(wěn)定,，上傳提示錯誤的情況,。

特別推薦其中的Taxon Set Enrichment Analysis模塊，直接提交物種列表（一般是找到的差異物種列表）,，可以直接在各種已有的相關(guān)性（人體基因-菌屬相關(guān),，生活方式-菌屬相關(guān)，疾病-菌屬相關(guān)）中進(jìn)行富集分析,，能很好的幫助判斷和提供差異菌群的具體關(guān)聯(lián)和證據(jù)支持,。

完整的支持分析包括：

可以直接生成下面的圖：

基本上常見的分析和圖都能在線實現(xiàn)。

 二,、 gcMeta 

另一個是https://gcmeta./,，是中科院微生物研究所搞的平臺，里面包括了宏基因組的樣本數(shù)據(jù)和在線分析平臺,，可以直接上傳原始數(shù)據(jù),，直接使用工具進(jìn)行在線分析，大部分常見工具都有,，也有一些流程,。

對于缺乏計算資源或想自己動手分析的朋友挺友好的，非常推薦試試看,。

最后,，幫大家整理了宏基因組可投稿的期刊，具體研究方向和影響因子見下表：