多快好省的宏基因組研究技巧

近年的研究熱點集中於環境和生物體相互作用的微生物群體，而大量複雜的微生物群體存在培養困難，構成複雜（包括細菌、古菌、真菌、原生生物、病毒甚至小型真核生物）。因此如何用高通量精準的了解這些群體的構成，基因功能分布以及具體的表達活性和代謝狀況成為首要問題。

高通量測序技術的發展，讓我們可以不經過培養，一次性了解微生物群落構成甚至基因代謝組成。

隨著技術的進步，檢測方法也逐漸豐富，對應的分析手段和軟體算法也逐步完善，使我們可以根據研究需要選擇不同的檢測和分析策略來獲得海量的數據並進行相應的研究分析。

01

簡 介

免於培養的微生物學研究方法主要基於測序，高通量測序使我們一次可以獲得整個微生物群體的數據信息，簡單來說包括兩種策略：

1、基於特定標記基因的擴增測序方案（常見的16s，ITs，18s或特定功能基因）

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

基於16s基因的分析方法

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

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

16s測序的分析流程如下圖，獲得序列經過聚類後獲得OTU或ASV，並得到相對豐度。

經過PICRUSt可以得到預測的基因分類豐度，進而進行alpha多樣性和Beta多樣性以及組間差異和相關性分析。

PICRSt的工作原理如下圖，將OTU表內16s序列進行對應物種16s拷貝數標準化後，將物種豐度乘以已經整理好的物種的基因注釋數表就獲得基因的預測豐度。

02

淺宏基因組

淺宏基因組測序方案是去年knights-lab在msystems上發表的針對16s解析度和宏基因組高成本之間的一個折中方案，通過降低測序深度，每個樣本50萬reads，但是物種的解析度並沒有低於一般宏基因組（普遍5~10G數據量）。

不通過拼接組裝，直接基於kraken2等kmer，或MetaPhlAn2等標記基因的參考基因組方法進行種屬豐度分類。結合其到菌株的物種分類和豐度數據可較16s方案下的PICRUST更加準確的預測基因構成。

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://doi.org/10.1128/mSystems.00069-18.

我們發現有些小夥伴的需求是：

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

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

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

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

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

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

處理方案如下：

• 高宿主含量DNA樣本（包括唾液、血液、肺泡灌洗液、腹水、陰道分泌物和黏膜類樣品）的取樣前處
  

• 將200微升唾液等體液樣本以10,000g離心8分鐘 

• 棄去上清液，通過移液將細胞沉澱重懸於200μl無菌水中，短暫渦旋，然後在室溫下靜置5分鐘，以滲透壓裂解哺乳動物細胞

• 添加終濃度為10μm的PMA（疊氮溴化丙錠）（向200μl樣品中添加10μl的0.2 mM PMA溶液），並將樣品短暫渦旋，然後在黑暗中於室溫溫育5分鐘

• 然後將樣品從標準臺式螢光燈放置在<20cm的冰上水平放置25分鐘，短暫離心並每5分鐘旋轉一次

• 完成後，可將樣品冷凍在−20°C或轉移到取樣管的儲存液中

Marotz CA, Sanders JG, Zuniga C, Zaramela LS, Knight R, Zengler K. Improving saliva shotgun metagenomics by chemical host DNA depletion. Microbiome. 2018;6(1):42. Published 2018 Feb 27. doi:10.1186/s40168-018-0426-3

可向上滑動

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

03

宏基因組

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

此外，對於缺乏深度研究和高質量參考基因組的樣本，如土壤和特殊環境下的樣本，宏基因組獲得的較為完整的基因組不僅可以豐富參考基因組資料庫，同時可以提供更加準確的物種分類。

因此，深度宏基因組測序是解析新環境樣本的核心方法，不過從單一樣本中重建出完整的菌株基因組有相當困難，一般需要較多樣本或設置梯度樣本從而利用更高深度和共同變化來獲取分箱信息，當然對應測序和分析成本會更高。

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

宏基因組報告中有哪些分析內容？

上圖可以快速預覽一下我們報告中的分析內容。

接下來，我們會詳細介紹這些內容是如何從原始數據開始一步步實現的，同時也會選取一些文章案例來給大家做詳細解讀，希望給大家帶來一些思路。

數據分析流程

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

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

對於宏基因組測序而言，最重要的就是獲得微生物群準確的物種構成及其豐度。

1.

物種構成

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

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

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

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

直接比對是基於參考數據的，那麼基於參考數據的物種構成分析主要有兩類方法：

一類是基於Kmer和LCA比對特徵來分析對應物種豐度，如kraken2等。

另一類是基於特徵標記基因進行分析的，如MetaPhlAn2等。

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

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

Centrifuge：比kraken2慢2x,內存使用少很多

Sourmash：類似CLARK,可以使用整個refseq作為資料庫。

主流的kraken2——快速、準確度高、內存要求高

目前主要使用kraken2為主，因為快速，準確度也相當不錯。不過，對於內存的要求較高，另外受資料庫本身質量影響較大，默認kraken2的參考資料庫只包括了細菌、古菌、病毒和人，還需要添加其他域的參考基因組。但涵蓋的測序參考種仍然有限，對於菌株水平的鑑定受一定影響。後續使用Bracken可以針對kraken2的比對結果進行計算相對豐度。

MetaPhlAn2——物種跨度大、實用

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

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

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

接下來，看一下我們報告中獲得的結果和圖：

使用Kraken2對其中的微生物進行物種注釋。我們的Kraken2使用的資料庫是由Refseq（2020.04.20）細菌，古細菌、真菌、原生動物和病毒庫以及GRCh38人類基因組構建的。

通過查詢資料庫序列中的每個k-mer，然後使用所得的LCA分類單元集確定序列的適當標籤，對序列進行分類。資料庫中沒有k-mers的序列不會被Kraken2分類。這裡我們是在使用k-mer=35的條件下進行物種注釋。

使用Bracken對物種注釋結果計算相對豐度。Bracken是一種高度精確的統計方法，可從宏基因組學樣本計算DNA序列中物種的豐度。Braken使用Kraken2分配的分類標籤來估計源自樣本中每種物種的讀數數量。

對物種注釋結果使用 KRONA 進行可視化展示。

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

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

如果嫌柱狀圖的展示方式單一，當然也可以有別的選擇。比如說以Circos的環圖形式展現：

也可以進行聚類分析：

有了這些數據我們就可以進行alpha多樣性（指每個樣本內部菌群多樣性）的分析了。

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

計算樣本之間的菌屬構成相似度：

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

Trukey多組間檢驗

LefSe分析

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

以上是關於物種組成部分，但是有些小夥伴會有這樣一些疑惑：物種構成變化很大怎麼辦？個體差異也很大？之類的諸多疑問。

是的，微生物群落一般對應特定的環境，其物種構成有時候變化迅速，而且個體或不同地點的構成差異極大。如人體的腸道菌群，個體之間的菌群構成差異很大，僅少量核心菌在絕大部分人的腸道內出現，個體特異性菌株也非常常見。那麼如此多樣性和複雜的構成如何應對相似的環境呢？

研究顯示不同的菌屬可能有著相似的基因或代謝能力，差異極大的種屬在基因功能層面可能有著相似的構成。因此，獲得微生物群的基因和功能代謝構成及分布對於解釋和了解微生物如何響應和適應環境就尤為重要。

2.

功能構成

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

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

方法一、通過直接基於reads的參考資料庫方法獲得

方法二、通過組裝後預測注釋基因並得到豐度

在具體展開方法之前，我們需要先了解關於基因功能的基本概念。

基因功能

每個菌的基因組中都包含大量的編碼基因（ORF）以及非編碼的RNA。這些基因之間又存在同源或序列相似性，達到一定相似程度的稱為同源基因（一般通過CD-hit聚類為unigene，gltA這類基因名稱，而資料庫中一般聚類為如uniref90，eggNOG_ortholog等不同相似度的非冗餘基因），這些同源基因除了序列相似同樣也有著相似的功能，基於其功能或具備的蛋白功能域可以進一步分類為基因家族（Pfam），酶（EC 1.4.1.13），代謝通路（ko:K00266），更進一步層層分類為GO或頂層代謝通路Metacyc或COG等。

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

主流的HUMAnN2——獲得基因和代謝通路豐度的同時可直接進行下遊分析

基於測序原始序列直接獲得基因構成豐度的軟體目前最主要的是HUMAnN2，其首先使用MetaPhlAn2進行物種分類（關於這個軟體我們在前面物種組成部分已經講過），並提取相應物種參考基因組用於比對，未比對上的用於進一步和uniref資料庫進行蛋白質序列比對。原理見下圖：

HUMAnN2的便利之處在於獲得基因和代謝通路豐度的同時可以直接進行下遊分析，將導出的表用於如LEFSE等差異分析，此外還可以反向給出不同樣本中每個基因或代謝通路裡的物種貢獻。

下圖是基於HUMAnN2的不同代謝通路的菌貢獻比例圖：

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

而另外一種方法是通過組裝獲得，我們在前面物種構成小節也已經提到過組裝分析，那麼這裡我們就組裝拼接分析這部分展開講解一下。

3.

基於組裝拼接的分析

什麼樣的條件下可以進行組裝分析？

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

組裝的流程是什麼樣的？

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

① 提取、測序

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

② 獲得contig序列

接著對序列經過質檢過濾處理後直接使用序列進行拼接，獲得contig序列，這時通常每個菌的基因組會有幾十到數千個contig片段，由於構成複雜，很多近緣菌之間的基因組存在大量相似序列，以及每種菌豐度都不一致，所以contig階段的片段仍然較多。

③ Binning分析 

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

④ 進一步質檢評估

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

借用一張分箱的說明PPT：

目前組裝contig方面比較好的軟體主要是SPAdes和MegaHIT。分箱方面MetaWRAP流程可以將整個組裝和分箱優化全部完成，包括前期質檢到組裝以及使用三種分箱方法concoct, metabat2和maxbin2，並最終進行合併提純優化，輸出最終的分箱。

同時還可以對每個分箱bins進行物種鑑定和定量，這樣我們就可以獲得基於拼接組裝後的物種豐度構成表，開展上述的物種多樣性和樣本差異統計分析。

⑤ 注釋

最後使用PROKKA進行基因預測，獲得的編碼序列我們經過進一步CD-Hit聚類去冗餘，然後使用eggNOG-mapper對其進行進一步的功能比對注釋。使用salmon完成基因的定量，這樣我們就得到基於組裝注釋的基因豐度數據了。之後就可以進行基因和功能層面的多樣性、構成以及樣本和組間差異分析。

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

KEGG

COG

eggNOG

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

除此之外，還可以使用其他的功能基因資料庫來進行進一步的基因注釋和分析。比如：

CAZy：

VFDB毒力因子注釋：

抗性基因注釋：

TCDB資料庫注釋：

PHI資料庫注釋：

BCGs分析：

以及基於antiSMASH和BiG-SCAPE來對代謝物的合成生物基因簇BCGs進行分析。

固定代謝能力評估：

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

當有了大量樣本的菌群構成豐度信息，以及各種基因和代謝豐度數據後，我們需要根據樣本的meta信息，基於不同分組，時間或環境因子等數據進行統計分析和檢驗，進而發現和探索可能的關聯以及背後的生物學意義。

4.

統計檢驗

那麼在面對宏基因組這類數據時在進行統計檢驗分析時需要注意什麼呢，應該採用哪些分析，並如何解讀這些結果呢？

首先，微生物組數據分析分為四大類：

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

為什麼可以過濾這兩類？

上述的兩類由於其攜帶的信息量和變化過小在進行組間比較統計檢驗的時候都建議過濾，因為要麼是汙染，要麼與差異無關。

宏基因組數據具有一些獨特的特徵，例如測序深度的巨大差異，稀疏性（包含許多零）和分布的巨大差異（過度分散）。在進行後續的統計檢驗之前建議針對不同的分析方法進行相應匹配的標準化處理。標準化包括：

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

轉換方法：包括處理稀疏性，組成性和數據中較大變化的方法。

那麼各種標準化方法是什麼，應該選擇哪種方法？

參考MicrobiomeAnalyst網站提供的信息，以下是一個簡短的介紹：

請注意，數據標準化主要用於可視數據探索，例如beta多樣性和聚類分析。有時候不使用標準化也能獲得最佳結果，比如：單變量統計和LEfSe。

同時，其他比較分析將使用其自己的特定標準化方法。例如，對metagenomeSeq使用累積總和定標（CSS）標準化，對edgeR應用M值的修剪均值（TMM）。

經常有小夥伴問，這個數據是用的什麼標準化？沒有做標準化怎麼辦？這類問題。

目前，尚無關於應使用標準化的共識性指南。建議大家可以探索不同的方法，然後目視檢查分離模式（即PCoA圖）以評估不同標準化程序對實驗條件或其他感興趣的宏基因組數據的影響。

有關這些方法的詳細討論，請參考使用者最近發表的兩篇論文 

①  Paul J. McMurdie等

（ 

② Jonathan Thorsen等

（   ）

以上是關於標準化的這部分內容需要了解的知識，接下來我們來看具體如何操作，怎麼得到那些圖表？它們分別代表什麼？

一般我們需要先進行探索性分析，也就是不設預訂的假設，首先從主成分分析結果中了解樣本的菌屬和基因的大概分布。

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

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

基於豐度圖來評估各樣本和分組的基本構成，如：

之後我們可以針對不同分組或處理之間的樣本進行統計檢驗，可以使用的檢驗方法包括兩組間的非參數統計檢驗T-test/ANOVA，3組以上組間統計檢驗可以使用Tukey test，其直接生成各組將的統計差異，並提供字母標註，直觀簡便，如：

具體的統計方法選擇可以參考下表：

除了常規的非參數檢驗外，包括metagenomeSeq和DEseq以及edgeR等統計方法包可以很好的分析組間差異特徵。LEfSe則一般用於尋找特徵標誌物。

那麼有了大量的差異特徵菌屬或基因之後，我們是否能基於這些差異菌屬有效的區分不同的分組呢，或構件一個模型來預測或分類呢？

這時候可以使用隨機森林（Random Forest）一類的決策樹機器學習模型，來利用這些差異特徵構建分類模型，並使用AUC等指標來評估基於這些模型的預測有效性和準確度（我們報告中如下圖）。

當然也可以使用其他更複雜的如深度學習等方法來構建分類模型。

除了性別、疾病、地點等分類差異之外，我們通常還有很多元數據，包括臨床指標或環境因子等信息，這些數據通常是連續型數值，對於這類數據我們可以進行相關性分析。

當然反過來，將菌群特徵作為表型也可以和如基因組的基因型或SNP構成來進行相關性分析。

對於菌群數據的相關性分析比較推薦：

SparCC方法，可以構建菌種或菌屬之間的相關性網絡，相對穩定。

對於與疾病或環境變量進行相關性分析可以使用：

Sperman秩相關分析。

另外RDA/CCA分析也可以有效的反映菌屬與環境因子等指標直接的關係（我們報告中如下圖）。

Mentel檢驗也可以用於判斷菌群構成特徵與單個或一組環境因子之間是否存在顯著相關。

要 點 

宏基因組從大量菌群和基因構成中尋找關聯是需要足夠的樣本量才能達到有效的統計效力，因為一次性獲得了大量的特徵數據，樣本量過少會帶來統計結論的無效，越是組內差異大，組間差異小的研究足夠大的樣本量才能得到可靠的結論。

一般動物樣本具有較好的背景可控，組內樣本數量建議至少6個，而人群研究由於背景複雜，個體多樣性高，一般建議組內50例以上較好。

以上看完後，你應該對宏基因組的數據分析流程有了整體的認識，也學會了相應的一些操作，但是不一定能直接從自己的這些數據、圖表中真正探索到和實際生物學相關的有價值的研究成果。

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

之後無論是剛開始的實驗設計，還是後面的分析，都會更加得心應手。

建議想好整個實驗思路再開始（或者也可以諮詢我們，我們專業的數據分析團隊會為你提供切實可行的項目方案）。

04

案例解析

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

第一項研究是關於肥胖患者減肥手術後的宏基因組和代謝數據的分析研究。

文獻來源：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名嚴重肥胖的受試者，他們是可調節胃束帶術（AGB，n = 20）或Roux-en-Y胃旁路術（RYGB，n = 41）的候選人。減肥手術後1、3和12個月隨訪24名受試者。使用宏基因組學測序和LC-MS分析腸道菌群和血清代謝組。另外納入了10人和147人分別作為宏基因組和代謝檢測的驗證集。

研究思路

這樣的設計分別有什麼作用？

第一點持續的動態採樣可以獲得持續變化情況，尤其是在一個特定變化後（減肥手術），持續的最終採樣有助於確認菌群的變化出現和特定事件或生理病理變化的前後，尤其是在確定因果中有重要幫助。

第二點獲得多維的數據有助於幫助我們全方位的了解菌群變化背後的帶來的生理和代謝變化以及之間的關聯。

第三點獨立驗證集的存在將大大增強研究的可信度，尤其是該研究納入的樣本量並不多，無法全面有效的控制無關因素，使得很多統計檢驗的效力無法顯現。這也導致該研究僅在基因總量和多樣性上獲得較好的重複效果，而更多的菌群精細特徵以及具體基因和代謝通路沒有得到深入分析。但是獨立驗證集保證了核心結論的可靠性和重複性，這點在宏基因組研究中非常重要。

從下圖可以看到研究針對樣本的總基因多樣性水平與生理指標和疾病狀態進行相關性分析和組間差異分析，圖中給出了顯著相關和差異的指標。

使用的統計檢驗方法是pearson和sperman相關和t-test以及Kruskal-Wallis檢驗。

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

進一步研究分析了術後特定變化模式的MAGs以及它們與代謝生理指標的相關性，見下圖：

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

研究的主要結論發現是低基因豐富度（LGC）存在於75％的患者中，並且與軀幹脂肪質量和合併症（2型糖尿病，高血壓和嚴重程度）增加相關。LGC改變了78種宏基因組種（MGS），其中50％與不良的身體成分和代謝表型有關。九種血清代謝產物（包括穀氨酸鹽，3-甲氧基苯基乙酸和L-組氨酸）和含有參與其代謝的蛋白質家族的功能模塊與低MGR密切相關。術後一年，BS會增加MGR，但儘管RYGB患者的代謝改善比AGB患者大，但術後一年的MGR仍然很低。

點 評 

總體而言該項研究可以使用淺宏基因組（在文章開頭第二部分詳細介紹過）來完成所有測序和分析，進一步擴大樣本數量，如果能同時獲得人的轉錄組數據甚至能更加明確的找到菌群變化與特定代謝通路的關聯關係。

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

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

文獻來源：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.

可以看到，研究同時記錄了糞便樣本的菌群宏基因組和每日的飲食記錄。研究的核心在於將每日飲食的食物通過營養構成進行量化，並構建類似物種進化樹的食物物候樹。

此外由於有每日的數據，可以通過前一日的食物與第二日的菌群數據進行時間序列分析，構建食物與菌之間的關聯以及時間相關性。

最後基於菌群數據和前一日飲食來構建模型預測判斷後一日的菌群狀態，幫助我們了解食物對於個體菌群的影響因素並實現定量和預測。

研究中對數據的處理過濾標準如下：刪除所有具有低讀取計數（每個樣品<23,500個讀取）的樣品。物種級別的分類表僅限於研究對象中至少存在25％的研究日，且在10%的研究樣本對象中發現的那些物種。

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

菌群和飲食以及營養構成的堆疊圖很好展現了變化和對應。

下面這張圖很好的顯示了飲食食物的變化與菌群變化之間的時間變化關係：

下面這張圖通過對每個人單獨進行菌屬與食物的Spearman相關，展現了菌與食物之間的關聯的個體化差異，在特定菌屬對應相同食物不同人會出現完全不同方向的變化，這也正是這項研究所揭示的，這種關聯關係的複雜性。

點 評 

本研究雖然有大量樣本，但並未進行組裝，而是直接使用了Refseq的細菌完成基因組序列作為參考。研究由於樣本數量眾多，測序深度也很有限，類似研究也可以使用淺宏基因組方式完成。

  < 案例三 >  類風溼關節炎的人群腸道微生物組 

接下來的一個研究是比較典型的宏基因組組裝並與疾病進行關聯分析的案例，研究的是日本人群類風溼關節炎的腸道微生物組的全基因組關聯研究。

文獻來源：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）進行了RA腸道微生物組的MWAS分析。MWAS由三個主要的生物信息學分析渠道（系統發育分析、功能基因分析、途徑分析）組成。

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

將QC後的序列直接比對到參考基因組，並根據基因組長度計算對應物種的相對豐度。

基因方面選擇denovo組裝，使用MegaHIT，然後再contig上完成ORF預測和CD-HIT聚類去冗餘，最後與UniRef和KEGG資料庫比對。

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

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

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

獲得相應數據後對相對豐度，該研究使用Box-Cox transformation對數據進行標準化，並過濾了一些低豐度的菌屬。Case-control的相關性分析使用的R的glm2模塊，將年齡、性別和測序上機分組作為協變量。

對於菌屬的關聯分析，最終將顯著性結果以火山圖和GraPhlAn圖的形式展現如下：

上面其中D圖使用每個菌的豐度進行UMAP分析，並映射關聯效應的展示比較有意思。

不過在基因層面上並未找到相應的關聯，可以看到下圖UniRef90的NMDS分布圖兩組之間無法有效區分，多樣性也沒有顯著差異。

點 評 

這項研究在菌層面發現了多個普雷沃氏菌屬的菌在日本人群中與類風溼性關節炎存在關聯，不過除此之外其他方面的發現並不多，僅找到一個基因存在顯著關聯，涉及的臨床調查也相對有限，且人群隊列數量不算多，並無獨立驗證集，因此亮點並不多。如果能納入免疫相應指標可能能研究的更細緻一些。

  < 案例四 ＞ 永凍土中參與有機物降級的關鍵菌群

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

文獻來源：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個永凍土菌群基因組。基於這些數據描繪了永凍土融化梯度下的菌群構成特徵，如下圖。

論文是2018年發表的，測序是在2011和12年測的，使用的是CLC Genomics Workbench 較早的4.4版分單樣本組裝，然後使用MetaBAT進行分箱，最後的標準是70%完成度和低於10%的汙染。

其中糖苷水解酶基因使用dbCAN資料庫的HMM進行預測，碳代謝使用KEGG數據。

研究還同時選擇了部分樣本進行了宏轉錄組和宏蛋白組的測序，對碳代謝同時結合轉錄組和蛋白組的數據，展現了特定通路下不同永凍土的菌群構成和表達豐度差異。

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

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

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

對關鍵物種的CH 4 :CO2濃度比相關代謝途徑的重建，以及相應基因的基因家族分析。

點 評 

總結一下這項研究，永凍土的菌群參考基因組數據缺乏，該研究從大量地點採集樣本重建了1500多個參考基因組。

首先從物種構成特徵上與永凍土融化階段特徵進行分析，並與重要環境因子進行分析，鎖定重要的特徵菌。

然後針對重要的代謝途徑和關鍵基因結合宏轉錄組和宏蛋白組全面解析代謝途徑的分布和差異變化。對關鍵的物種重建了相關代謝途徑並對其相關基因家族進行分析。

研究基本上從頭構建了一個生態環境下的菌群結構數據，並利用獲得的基因組深度解析特定代謝途徑和基因的構成和表達變化，應該說既寬又深。很多樣本採集和測序是2011年和 2012年就開展的，雖然測序技術遠不如現在成本低和成熟，但是其獨特的研究對象和全面深入的分析仍然使整項研究和目前的一些研究相比完成的更加出色。

p.s. 以上展示的圖表，我們都可以幫你實現~ 

05

工具分享

 一、 MicrobiomeAnalyst 

網址：https://www.microbiomeanalyst.ca/，只需要biom文件或豐度表就可以進行絕大部分統計檢驗分析，而且生成圖表完整，可以直接使用。偶爾會有伺服器不穩定，上傳提示錯誤的情況。

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

完整的支持分析包括：

可以直接生成下面的圖：

基本上常見的分析和圖都能在線實現。

 二、 gcMeta 

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

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

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