上回給大家講述了16S測序分析流程，應各位小伙伴們的要求，本期的宏基因組測序分析來啦~

　　- 概述

　　- 宏基因組測序?qū)嶒灹鞒?/p>

　　- 宏基因組分析流程

　　- 文章分析思路

　　- 常見問題

　　- 相關名詞解釋

　　- 參考文獻

　　16S rDNA擴增子測序（16S rDNA amplicons sequencing）和鳥槍法宏基因組測序（whole metageome shotgun sequencing）都是研究微生物組的重要方法，那么問題來了：這兩者到底有什么區(qū)別呢？什么情況下做16S測序？什么情況下做宏基因組測序？

　　概述

　　先要來說說什么是宏基因組測序：

　　宏基因組 ( Metagenome) (也稱元基因組)。是由 Handelsman 等 1998 年提出的新名詞, 其定義為“the genomes of the total microbiota found in nature” , 即生境中全部微小生物遺傳物質(zhì)的總和。它包含了可培養(yǎng)的和未可培養(yǎng)的微生物的基因, 目前主要指環(huán)境樣品中的細菌和真菌的基因組總和。

　　宏基因組學 (或元基因組學, metagenomics)就是一種以環(huán)境樣品中的微生物群體基因組為研究對象, 以功能基因篩選和/或測序分析為研究手段, 以微生物多樣性、 種群結構、 進化關系、 功能活性、 相互協(xié)作關系及與環(huán)境之間的關系為研究目的的新的微生物研究方法。

　　宏基因組測序（Metagenomics Sequencing）是對環(huán)境樣品中全部微生物的總DNA進行高通量測序，主要研究微生物種群結構、基因功能、微生物之間的相互協(xié)作關系以及微生物與環(huán)境之間的關系。宏基因組測序研究擺脫了微生物分離純培養(yǎng)的限制，擴展了微生物資源的利用空間，為環(huán)境微生物群落的研究提供了有效工具。

　　高通量測序平臺

　　16S擴增子測序和宏基因組測序的主要區(qū)別如下：

　　測序原理不同

　　16S rDNA基因存在于所有細菌的基因組中，具有高度的保守性。該序列包含9個高變區(qū)和10個保守區(qū)（如下圖），通過對某一段高變區(qū)序列（V4區(qū)或V3-V4區(qū)）進行PCR擴增后進行測序，得到300-500bp左右的序列。

　　細菌16S rDNA基因

　　宏基因組測序則是將微生物基因組DNA隨機打斷成500bp的小片段，然后在片段兩端加入通用引物進行PCR擴增測序，再通過組裝的方式，將小片段拼接成較長的序列。

　　研究目的不同

　　16S測序主要研究群落的物種組成、物種間的進化關系以及群落的多樣性。

　　宏基因組測序在16S測序分析的基礎上還可以進行基因和功能層面的深入研究，宏基因組測序可以回答這樣的問題"who is there?"和 "what are they doing?"。

　　物種鑒定分辨率不同

　　16S測序得到的序列很多注釋不到種水平，而宏基因組測序則能鑒定微生物到種水平甚至菌株水平。

　　對于16S測序而言，任何一個高變區(qū)或幾個高變區(qū)，盡管具有很高的特異性，但是某些物種（尤其是分類水平較低的種水平）在這些高變區(qū)可能非常相近，能夠區(qū)分它們的特異性片段可能不在擴增區(qū)域內(nèi)。

　　宏基因組測序通過對微生物基因組隨機打斷，并通過組裝將小片段拼接成較長的序列。因此，在物種鑒定過程中，宏基因組測序具有較高的優(yōu)勢。

　　Tips：通常情況下，建議同時結合宏基因組測序和16S測序兩種技術手段，可以更高效、更準確地研究微生物群落組成結構、多樣性以及功能情況。

　　如果樣本污染宿主DNA比較嚴重，例如腸道粘膜樣本，直接宏基因組測序會產(chǎn)生大量的宿主污染，為了降低實驗成本，可以使用16S測序。

　　如果想快速鑒定未知病原感染，直接通過metagenome測序可以鑒定是細菌、真菌或者是病毒感染。

　　宏基因組測序?qū)嶒灹鞒?/p>

　　從環(huán)境（如土壤、海洋、淡水、腸道等）中采集實驗樣本，將原始采樣樣本或已提取的 DNA 樣本低溫運輸（0℃以下），對樣品進行樣品檢測。檢測合格的 DNA 樣品，進行文庫構建以及文庫檢測，檢測合格的文庫將采用 Illumina 高通量測序平臺進行測序，測序得到的下機數(shù)據(jù)(Raw Data)將用于后期信息分析。

　　為保證測序數(shù)據(jù)的準確性、可靠性，對樣品檢測、建庫、測序每一個生產(chǎn)步驟都嚴格把控，從根本上確保高質(zhì)量數(shù)據(jù)的產(chǎn)出，具體的實驗流程圖如下：

　　實驗流程圖

　　DNA樣品檢測

　　對 DNA 樣品的檢測主要包括 3 種方法：

　　1) 瓊脂糖凝膠電泳（AGE）分析 DNA 的純度和完整性；

　　2) Nanodrop 檢測 DNA 的純度（OD260/280 比值）；

　　3) Qubit 對 DNA 濃度進行精確定量；

　　文庫構建及質(zhì)檢

　　檢測合格的 DNA 樣品用超聲波破碎儀隨機打斷成長度約為350bp的片段，經(jīng)末端修復、3’端加A、加測序接頭、純化、片段選擇、PCR 擴增等步驟完成整個文庫制備。

　　文庫構建完成后，先用電泳及Nanodrop進行初步定量，對濃度>=15ng/ul的文庫進行Qubit定量，用毛細管電泳對文庫的插入片段大小進行檢測，插入片段大小符合預期后，使用qPCR方法對文庫的有效濃度進行準確定量（文庫有效濃度＞3nM），以保證文庫上機質(zhì)量。

　　上機測序

　　建庫質(zhì)檢合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求，混合后進行Illumina測序。

　　宏基因組分析流程

　　下面先放一張宏基因組分析流程圖，供小伙伴們快速了解一下。

　　信息分析流程圖

　　測序數(shù)據(jù)預處理

　　?常規(guī)數(shù)據(jù)質(zhì)控工具

　　?FASTX-Toolkit (http://hannonlab./fastx_toolkit/commandline.html)

　　?Cutadapt (https://github.com/marcelm/cutadapt)

　　?Trimmomatic (http://www./cms/?page=trimmomatic)

　　?Sickle (https://github.com/najoshi/sickle)

　　?去除宿主污染

　　A member of the SOAP (Short Oligonucleotide Analysis Package). It is an updated version of SOAP software for short oligonucleotide alignment.

　　?Bowtie (http:///projects/bowtie-bio/files/)

　　An ultrafast, memory-efficient short read aligner, which aligns short DNA sequences (reads) to the human genome at a rate of over 25 million 35-bp reads per hour.

　　采用 Illumina 測序平臺測序獲得的原始數(shù)據(jù)(Raw Data)存在一定比例低質(zhì)量數(shù)據(jù)，里面含有帶接頭的、重復的，以及測序質(zhì)量很低的reads，這些 reads 會影響組裝和后續(xù)分析，為了保證后續(xù)分析的結果準確可靠，需要對原始的測序數(shù)據(jù)進行預處理，獲取用于后續(xù)分析的有效數(shù)據(jù)(Clean Data)。

　　可以使用過濾軟件 ，可以從任意一段切除低質(zhì)量的堿基，同時還支持滑窗過濾，根據(jù)情況設定滑窗的大小和閾值，當滑窗內(nèi)的堿基質(zhì)量與設定的閾值進行比較，如果數(shù)值低于閾值則切除整個滑窗的堿基。高通量測序一般會包括接頭序列以及引物片段，可以使用 來去除這些序列。

　　例如具體處理步驟如下：

　　1) 去除所含低質(zhì)量堿基（質(zhì)量值≤38）超過一定比例（默認設為 40bp）的 reads；

　　2) 去除 N 堿基達到一定比例的 reads（默認設為10bp）；

　　數(shù)據(jù)處理信息統(tǒng)計表

　　注：RawData 表示下機原始數(shù)據(jù)；CleanData 表示過濾得到的有效數(shù)據(jù)；CleanQ20，CleanQ30 表示 CleanData 中測序錯誤率小于0.01(質(zhì)量值大于 20)和 0.001(質(zhì)量值大于 30)的堿基數(shù)目的百分比；Clean_GC(%) 表示 CleanData 中堿基的 GC 含量；Effective(%) 表示有效數(shù)據(jù)( CleanData )與原始數(shù)據(jù)( RawData )的百分比。

　　物種注釋metaphlan2

　　主頁：http://segatalab.cibio./tools/metaphlan2/

　　從Clean read出發(fā)，使用metaphlan2軟件分析，獲得不同分類層級的物種豐度表。

　　MetaPhlAn2是分析微生物群落(細菌、古菌、真核生物和病毒)組成的工具，它在宏基因組研究中非常有用，只需一條命,即可獲得微生物的物種豐度信息。同時配合自帶的腳本可進一步統(tǒng)計和可視化。

　　MetaPhlAn2整理了超過17000個參考基因組，包括13500個細菌和古菌，3500個病毒和110種真核生物，匯編整理了>1百萬類群特異的標記基因，可以實現(xiàn)：

　　精確的分類群分配

　　準確估計物種的相對豐度

　　種水平精度

　　株鑒定與追蹤

　　超快的分析速度

　　Duy Tin Truong, Eric Franzosa, Timothy L Tickle, Matthias Scholz, George Weingart, Edoardo Pasolli, Adrian Tett, Curtis Huttenhower, and Nicola SegataMetaPhlAn2 for enhanced metagenomic taxonomic profilingNature Methods 12, 902–903 (2015)

　　看看軟件主要能干什么

　　分析人類微生物組(HMP)數(shù)據(jù)，在樣本間進行層級聚類

　　人類微生物組數(shù)據(jù)多時間點普氏（Prevotella copri ）菌株水平的指紋

　　使用實例

　　定量物種譜

　　一條命令就可以從原始數(shù)據(jù)得到物種的相對豐度??！

　　最常用的是使用雙端壓縮fastq文件，參數(shù)--nproc調(diào)用多個線程，并輸出比較結果

　　例如人類腸道的數(shù)據(jù)，一般幾個小時內(nèi)就出結果了。

　　輸出結果為各層級物種相對豐度值，可以直接作為lefse的輸入文件。

　　批量處理多個樣品

　　假如我們有1-20一共20個樣本，調(diào)用for循環(huán)并轉(zhuǎn)后臺并行處理

　　樣品物種譜合并

　　mergemetaphlantables.py命令可以將每個樣品結果表合并，程序位于程序的utils目錄中，請自行添加環(huán)境變量或使用絕對路徑。合并時支持輸入文件多個文件空格分隔，或使用通配符(如下)。可結合sed刪除樣本名中共有部分，精簡樣品名方便可視化簡潔美觀

　　獲得了矩陣表，下面可以進行各種統(tǒng)計分析與可視化啦！

　　種相對豐度概況

　　從不同分類層級的相對豐度表出發(fā)，選取出在各樣品中的最大相對豐度排名前 15 的物種，繪制出各樣品對應的物種注釋結果在不同分類層級上的相對豐度柱形圖。

　　門水平菌落結構柱狀圖

　　基于GraPhlAn的分類學組成信息可視化

　　注：分類等級樹展示了樣本總體中，從門到屬（從內(nèi)圈到外圈依次排列）所有分類單元（以節(jié)點表示）的等級關系，節(jié)點大小對應于該分類單元的平均相對豐度，字母上的陰影顏色同對應節(jié)點顏色一致。

　　功能注釋humann2

　　Species-level functional profiling of metagenomes and metatranscriptomes

　　Nature Methods, Article, 2018-10-30

　　第一作者： Eric A. Franzosa, Lauren J. McIver

　　通訊作者：Curtis Huttenhower

　　主要單位：哈佛醫(yī)學院統(tǒng)計系；哈佛和麻省理工博德(Broad)研究所

　　humman2是一套分析流程，它包括調(diào)用metaphlan流程來分析物種組成，和自身分析功能基因和代謝通路組成。HUMAnN2 Workflow

　　Metagenome 組裝和注釋

　　Metagenome 組裝

　　兩種常見策略：

　　1 基于序列overlap關系進行拼接

　　2 基于de Bruijn圖進行組裝

　　由于現(xiàn)階段的主流測序方法是二代短片段測序，序列短而且數(shù)目龐大，如果利用overlap關系直接進行組裝，這要求每對reads之間都進行一次序列比較，這會很耗費時間，而且結果并不可靠。為迎合二代測序的特點，一種基于k-mer的de Bruijn組裝策略則成為更有效的解決方法。

　　能用于宏基因組組裝的軟件有很多，如metavelvet，SOAPdenovo，megahit，ibda-ud，metaspades等。

　　1、SOAPdenovo

　　SOAPdenovo is a novel short-read assembly method that can build a de novo draft assembly for the human-sized genomes， which specially designed to assemble Illumina GA short reads.

　　2、MEGAHIT

　　MEGAHIT is a single node assembler for large and complex metagenomics NGS reads, such as soil. Compare to SOAPdenovo, it generates longer contigs and consumes less memory.

　　3、IDBA-UD

　　IDBA-UD is a iterative De Bruijn Graph De Novo Assembler for Short Reads Sequencing data with Highly Uneven Sequencing Depth.

　　4、SPAdes和metaSPAdes

　　An assembly toolkit containing various assembly pipelines, which works with Illumina or IonTorrent reads and is capable of providing hybrid assemblies using PacBio, Oxford Nanopore and Sanger reads.

　　如下是三種主流軟件分析，運行所消耗時間、內(nèi)存比較:

　　https://2017-cicese-metagenomics./en/latest/assemble.html

　　目前MEGAHIT在現(xiàn)有組裝軟件中，資源消耗基本上是最低的，因此很適合宏基因組中的復雜環(huán)境樣品。

　　還有上面介紹過的軟件SPAdes，無論單菌、宏基因組還是宏病毒組都表現(xiàn)不錯 metaSPAdes: a new versatile metagenomics assembler）中顯示即使在復雜環(huán)境（土壤），組裝效果也大大優(yōu)于megahit、IDBA-UD等，但是沒有megahit資源消耗低。

　　以下使用 SOAPdenovo[1]進行Metagenome 組裝為例。

　　從各樣品質(zhì)控后的Clean reads出發(fā)，組裝主要分為3步：單樣本組裝，多樣本組裝結果合并和豐度過濾。

　　1) 經(jīng)過預處理后得到Clean Data，使用SOAPdenovo[1]組裝軟件進行組裝分析(Assembly Analysis)；

　　2) 對于單個樣品，首先選取一個K-mer(默認選取55)進行組裝,得到該樣品的組裝結果；組裝參數(shù)：-d1,-M3,-R,-u,-F

　　3) 將組裝得到的Scaffolds從N連接處打斷，得到不含N的序列片段，稱為Scaftigs(i.e.,continuous sequences within scaffolds)；

　　4) 將各樣品質(zhì)控后的CleanData采用SoapAligner軟件比對至各樣品組裝后的Scaftigs上，獲取未被利用上的PE reads；比對參數(shù)：-u,-2,-m200

　　5) 將各樣品未被利用上的reads放在一起，進行混合組裝，考慮到計算消耗和時間消耗，只選取一個kmer進行組裝(默認-K55)，其他組裝參數(shù)與單樣品組裝參數(shù)相同；

　　6) 將混合組裝的Scaffolds從N連接處打斷，得到不含N的Scaftigs序列；

　　7) 對于單樣品和混合組裝生成的Scaftigs，過濾掉500bp以下的片段，并進行統(tǒng)計分析和后續(xù)基因預測；

　　各樣品組裝結果基本信息統(tǒng)計

　　說明：SampleID為樣品名稱；Total Length表示組裝得 Scaftigs 的總長； Scaftigs 表示組裝得到的Scaftigs 總條數(shù)； N50，N90 表示將 Scaftigs 按照長度進行排序，然后由長到短加和，當加和值達到 Scaftigs 總長的 50%，90%時的 Scaftigs 的長度值；Max len 表示組裝得到的最長 Scaftigs 的長度值；Min len表示組裝得到的最長 Scaftigs 的長度值。

　　樣品的scaftigs長度分布統(tǒng)計如下：

　　單個樣品scaftigs長度分布統(tǒng)計

　　各樣品scaftigs數(shù)目統(tǒng)計：

　　所有樣品scaftigs數(shù)目統(tǒng)計圖

　　組裝結果質(zhì)量評估工具 MetaQUAST

　　基因預測及豐度分析

　　基因預測及豐度分析基本步驟

　　1) 從各樣品及混合組裝的 Scaftigs（>=500bp）出發(fā)，采用MetaGeneMark進行ORF (OpenReading Frame) 預測，并從預測結果出發(fā)，過濾掉長度小于 100nt的信息；預測參數(shù)：采用默認參數(shù)進行

　　3) 采用SoapAligner，將各樣品的 Clean Data 比對至初始 gene catalogue，計算得到基因在各樣品中比對上的reads 數(shù)目；比對參數(shù)：-m 200, -x 400, identity ≥ 95%

　　4）過濾掉在各個樣品中支持 reads 數(shù)目

　　5）從比對上的 reads 數(shù)目及基因長度出發(fā)，計算得到各基因在各樣品中的豐度信息

　　6）基于 gene catalogue 中各基因在各樣品中的豐度信息，進行基本信息統(tǒng)計，core-pan 基因分析，樣品間相關性分析，及基因數(shù)目韋恩圖分析。

　　Gene catalogue 基本信息

　　說明：ORFs NO. 表示 gene catalogue 中基因的數(shù)目；integrity:start 表示只含有起始密碼子的基因數(shù)目及百分比；integrity:end 表示只含有終止密碼子的基因數(shù)目及百分比；integrity:none 表示沒有起始密碼子也沒有終止密碼子的基因數(shù)目及百分比；integrity:all 表示完整基因（既有起始密碼子也有終止密碼子）數(shù)目的百分比；Total Len.(Mbp) 表示 gene catalogue 中基因的總長，單位是百萬；Average Len. 表示 gene catalogue 中基因的平均長度；GC Percent 表示預測的 gene catalogue 中基因的整體 GC 含量值。

　　gene catalogue 長度分布統(tǒng)計

　　core-pan 基因分析

　　從基因在各樣品中的豐度表出發(fā)，可以獲得各樣品的基因數(shù)目信息，通過隨機抽取不同數(shù)目的樣品，可以獲得不同數(shù)目樣品組合間的基因數(shù)目，由此我們構建和繪制了 Core 和 Pan 基因的稀釋曲線，圖片展示如下:

　　core-pan 基因稀釋曲線

　　組間基因數(shù)目差異分析

　　為了考察組與組間的基因數(shù)目差異情況，繪制了組間基因數(shù)目差異箱圖，展示結果如下：

　　組間基因數(shù)目差異小提琴圖

　　說明：圖中，橫坐標為各個分組信息；縱坐標為基因數(shù)目。

　　基因數(shù)目韋恩圖分析

　　為了考察指定樣品（組）間的基因數(shù)目分布情況，分析不同樣品（組）之間的基因共有、特有信息，繪制了韋恩圖(Venn Graph)或花瓣圖，展示結果如下：

　　基因數(shù)目韋恩圖（花瓣圖）分析

　　說明：圖中，每個圈代表一個樣品；圈和圈重疊部分的數(shù)字代表樣品之間共有的基因個數(shù)；沒有重疊部分的數(shù)字代表樣品的特有基因個數(shù)。

　　基于基因豐度的樣品間相關性分析

　　樣品間基因豐度相關性是檢驗實驗可靠性和樣本選擇是否合理性的重要指標。相關系數(shù)越接近1，表明樣品之間基因豐度模式的相似度越高。

　　說明：圖中不同顏色代表相關系數(shù)的高低，相關系數(shù)與顏色間的關系見有圖圖例說明。

　　基于物種豐度的降維分析1) PCA分析

　　主成分分析PCA（Principal component analysis）是一種研究數(shù)據(jù)相似性或差異性的可視化方法，通過一系列的特征值和特征向量進行排序后，選擇主要的前幾位特征值，采取降維的思想，PCA 可以找到距離矩陣中最主要的坐標，結果是數(shù)據(jù)矩陣的一個旋轉(zhuǎn)，它沒有改變樣品點之間的相互位置關系，只是改變了坐標系統(tǒng)。PCA 可以觀察個體或群體間的差異。下圖每一個點代表一個樣本，相同顏色的點來自同一個分組，兩點之間距離越近表明兩者的群落構成差異越小。

　　PCA圖

　　說明：種水平PCA分析，橫坐標表示第一主成分，百分比則表示第一主成分對樣品差異的貢獻值；縱坐標表示第二主成分，百分比表示第二主成分對樣品差異的貢獻值；圖中的每個點表示一個樣品，同一個組的樣品使用同一種顏色表示。

　　2) NMDS分析

　　NMDS（Non-metric Multidimensional scaling）非度量多維尺度分析是一種將多維空間的研究對象（樣本或變量）簡化到低維空間進行定位、分析和歸類，同時又保留對象間原始關系的數(shù)據(jù)分析方法。適用于無法獲得研究對象間精確的相似性或相異性數(shù)據(jù)，僅能得到他們之間等級關系數(shù)據(jù)的情形。其基本特征是將對象間的相似性或相異性數(shù)據(jù)看成點間距離的單調(diào)函數(shù)，在保持原始數(shù)據(jù)次序關系的基礎上，用新的相同次序的數(shù)據(jù)列替換原始數(shù)據(jù)進行度量型多維尺度分析。其特點是根據(jù)樣品中包含的物種信息，以點的形式反映在多維空間上，而對不同樣品間的差異程度，則是通過點與點間的距離體現(xiàn)的，最終獲得樣品的空間定位點圖。

　　NMDS圖

　　組間差異物種的LEfSe分析

　　為了篩選組間具有顯著差異的物種Biomarker，首先通過秩和檢驗的方法檢測不同分組間的差異物種并通過LDA（線性判別分析）實現(xiàn)降維并評估差異物種的影響大小，即得到LDA score；組間差異物種的LEfSe分析結果包括三部分，分別是LDA值分布柱狀圖，進化分支圖（系統(tǒng)發(fā)育分布）和組間具有統(tǒng)計學差異的Biomarker在不同組中豐度比較圖。差異物種的LDA值分布圖和進化分支圖如下：

　　常用功能數(shù)據(jù)庫注釋

　　功能注釋主要是基于功能同源的序列往往具有序列相似性的原理，將去冗余后的基因蛋白序列，與不同的蛋白功能數(shù)據(jù)庫進行序列比對， 然后用比對到的序列的功能作為目標序列的功能。

　　從 gene catalogue 出發(fā)，進行代謝通路（KEGG）[6,7]，同源基因簇（eggNOG）[8]，碳水化合物酶（CAZy）的功能注釋和豐度分析。基于物種豐度表和功能豐度表，可以進行豐度聚類分析，PCA和NMDS 降維分析，Anosim分析，樣品聚類分析；當有分組信息時，可以進行Metastat和LEfSe多元統(tǒng)計分析以及代謝通路比較分析，挖掘樣品之間的物種組成和功能組成差異。

　　目前常用的功能數(shù)據(jù)庫主要有：

　　KEGG，全稱Kyoto Encyclopedia of Genes and Genomes，是一個關于基因功能注釋方面的綜合性數(shù)據(jù)庫，包括基因的功能，分類，代謝通路(KEGG Pathway數(shù)據(jù)庫， 是KEGG最核心的數(shù)據(jù)庫)等諸多方面的信息。

　　KEGG 數(shù)據(jù)庫于 1995 年由 Kanehisa Laboratories 推出 0.1 版，目前發(fā)展為一個綜合性數(shù)據(jù)庫，其中最核心的為KEGG PATHWAY 和 KEGG ORTHOLOGY 數(shù)據(jù)庫。在 KEGG ORTHOLOGY 數(shù)據(jù)庫中，將行使相同功能的基因聚在一起，稱為Ortholog Groups (KO entries)，每個 KO 包含多個基因信息，并在一至多個 pathway 中發(fā)揮作用。

　　KEGG Pathway數(shù)據(jù)庫將生物代謝通路劃分為6大類(A級分類)，分別為：新陳代謝（Metabolism）、遺傳信息處理（Genetic Information Processing）、 環(huán)境信息處理（Environmental Information Processing）、細胞過程（Cellular Processes）、生物體系統(tǒng)（Organismal Systems）、 人類疾病（Human Diseases），其中每大類又被系統(tǒng)分類為B、C、D3個級別。 其中B級分類目前包括有 43 種子功能；C級分類即為代謝通路圖；D級分類為每個代謝通路圖的具體注釋信息。

　　eggNOG 數(shù)據(jù)庫是利用 Smith-Waterman 比對算法對構建的基因直系同源簇 (Orthologous Groups) 進行功能注釋， eggNOG（v4.5）目前涵蓋了2031個物種，構建了包含25個大類，約19萬個直系同源簇。

　　CAZy 數(shù)據(jù)庫是研究碳水化合物酶的專業(yè)級數(shù)據(jù)庫，主要涵蓋 6 大功能類：糖苷水解酶（Glycoside Hydrolases ，GHs）， 糖基轉(zhuǎn)移酶（Glycosyl Transferases，GTs），多糖裂合酶（Polysaccharide Lyases，PLs），碳水化合物酯酶（Carbohydrate Esterases，CEs）， 輔助氧化還原酶(Auxiliary Activities , AAs)和碳水化合物結合模塊（Carbohydrate-Binding Modules， CBMs）。 其中每一個大類有可以分類很多小的家族，比如CE1，CE2等等，注釋結果中的CE0表示沒有小家族分類的結果。

　　功能注釋基本步驟

　　1）使用DIAMOND軟件將 Unigenes 與各功能數(shù)據(jù)庫進行比對（blastp，evalue ≤ 1e-5）；

　　2）比對結果過濾：對于每一條序列的 比對結果，選取 score 最高的比對結果（one HSP > 60 bits）進行后續(xù)分析；

　　3）從比對結果出發(fā)，統(tǒng)計不同功能層級的相對豐度（各功能層級的相對豐度等于注釋為該功能層級的基因的相對豐度之和），其中，KEGG 數(shù)據(jù)庫劃分為 5 個層級，eggNOG 數(shù)據(jù)庫劃分為 3 個層級，CAZy 數(shù)據(jù)庫劃分為 3 個層級，各數(shù)據(jù)庫的詳細劃分層級如下所示：

　　4）從功能注釋結果及基因豐度表出發(fā)，獲得各個樣品在各個分類層級上的基因數(shù)目表，對于某個功能在某個樣品中的基因數(shù)目，等于在注釋為該功能的基因中，豐度不為 0 的基因數(shù)目；

　　5）從各個分類層級上的豐度表出發(fā)，進行注釋基因數(shù)目統(tǒng)計，相對豐度概況展示，豐度聚類熱圖展示，PCA和NMDS降維分析，基于功能豐度的Anosim組間（內(nèi)）差異分析，代謝通路比較分析，組間功能差異的Metastat和LEfSe分析。

　　注釋基因數(shù)目統(tǒng)計

　　從 Unigenes 注釋結果出發(fā)，繪制各個數(shù)據(jù)庫的注釋基因數(shù)目統(tǒng)計圖，展示結果如下圖所示：

　　各數(shù)據(jù)庫注釋基因數(shù)目統(tǒng)計圖

　　說明：從上至下依次為CAZy，eggNOG 的 Unigenes 注釋數(shù)目統(tǒng)計圖。橫坐標軸是各數(shù)據(jù)庫中 level1 各功能類的代碼，代碼的解釋見對應的圖例說明。

　　功能相對豐度概況

　　從各數(shù)據(jù)庫 level1 的相對豐度表出發(fā)，繪制出各個數(shù)據(jù)庫中，各樣品對應的在 level1 層級上的豐度統(tǒng)計圖。

　　KEGG:

　　Eggnog:

　　功能注釋在 level1 上的相對豐度柱形圖

　　說明：從上至下依次為 Kegg,eggNOG 的結果展示。縱軸表示注釋到某功能類的相對比例；橫軸表示樣品名稱；各顏色區(qū)塊對應的功能類別見右側(cè)圖例。

　　功能相對豐度聚類分析

　　根據(jù)所有樣品在各個數(shù)據(jù)庫中的功能注釋及豐度信息，選取豐度排名前 35 的功能及它們在每個樣品中的豐度信息繪制熱圖，并從功能差異層面進行聚類。

　　Kegg:

　　Eggnog:

　　功能豐度聚類熱圖

　　說明：從上至下依次為 KEGG, eggNOG 的結果展示。橫向為樣品信息；縱向為功能注釋信息；

　　基于功能豐度的降維分析

　　基于不同數(shù)據(jù)庫在各個分類層級的功能豐度進行 PCA 和 NMDS 降維 分析，如果樣品的功能組成越相似，則它們在降維圖中的距離越接近?；贙EGG的KO功能豐度進行PCA和NMDS 分析的結果展示如下：

　　抗生素抗性基因注釋

　　抗生素的濫用導致人體和環(huán)境中微生物群落發(fā)生不可逆的變化，對人體健康和生態(tài)環(huán)境造成風險，因此抗性基因的相關研究受到了研究者的廣泛關注。

　　目前，抗生素抗性基因的研究主要利用ARDB（Antibiotic Resistance Genes Database）數(shù)據(jù)庫和CARD（Comprehensive Antibiotic Resistance Database）數(shù)據(jù)庫來注釋抗生素抗性基因。通過該數(shù)據(jù)庫的注釋，可以找到抗生素抗性基因（Antibiotic Resistance Genes ，ARGs）及其抗性類型（Antibiotic Resistance Class）以及這些基因所耐受的抗生素種類（ Antibiotic）等信息。

　　抗性基因豐度概況

　　依據(jù)耐受的抗生素種類豐度結果，提取豐度前15抗生素種類繪制柱狀圖如下：

　　不同抗生素抗性基因在各樣品中的豐度柱形圖

　　注：橫軸為樣品，縱軸為抗生素抗性基因

　　抗性基因類型豐度聚類分析

　　根據(jù)抗性基因類型的豐度表，繪制聚類熱圖：

　　注：橫軸為樣品，縱軸為抗性基因類型，不同顏色達標圖例中對應的豐度范圍。

　　組間抗性基因數(shù)目差異分析

　　注：橫軸為分組，縱軸為抗性基因類型數(shù)目

　　物種間相關性網(wǎng)絡分析

　　機器學習算法

　　分類器是通過特征工程中的特征選擇的方法挑選出若干不同級別的bio-markers，并根據(jù)bio-markers來構建的一種分類模型，它采用分類評價指標來評估模型效果的好壞，目的是用于識別疾病與健康人群。

　　分類器是來自機器學習中的分類算法，包括Logistictic回歸、Support Vector Machine (SVM)、Random Forest、人工神經(jīng)網(wǎng)絡、樸素貝葉斯等。除了分類算法，特征選擇對分類性能也會有直接的影響，特征選擇就是在我們獲得的給定特征集中選擇出與分類相關的特征子集的過程。

　　分類評價指標是用來評價分類器分類性能好壞的一種指標，常見的評價指標有正確率、錯誤率、靈敏度、特異度、精度、Matthews correlation coefficient (MCC)、Receiver operating characteristic (ROC)。其中MCC和ROC不受分組樣本數(shù)量影響；ROC是根據(jù)一系列不同的二分類方式，以真陽性率(靈敏度)為縱坐標，假陽性率(1-特異度)為橫坐標繪制的曲線，直觀描繪了分類器在TP和FP間的trade-off。

　　AUC (area under curve)是ROC曲線下方面積之和，它以數(shù)值的形式來評估分類器的好壞，AUC的取值范圍在(0.5, 1.0)范圍內(nèi)，AUC的值越大，分類的性能就越好。若分類評價指標最終顯示分類器性能不佳，則需重新構建分類器直至分類性能達到穩(wěn)定為止。

　　分類器的引入對于構建用于臨床診斷的基于腸道微生物的檢測方法具有重要意義。

　　文章分析思路

　　文章分析思路：整體概述——物種組成——功能組成——關聯(lián)分析——因果關系驗證

　　文章主要結果

　　樣本背景、數(shù)據(jù)評估、物種、功能描述

　　物種組成整體差異和分組具體差異

　　功能組成整體差異和分組具體差異

　　環(huán)境因子與微生物組的關系

　　菌群或單菌移植到無菌小鼠體內(nèi)，驗證因果關系

　　常見問題

　　16S擴增子和宏基因組分析結果為什么存在差別？

　　1、兩者的實驗方法存在較大差異：

　　16S是先擴增后測序，而且不同物種DNA的擴增有偏好；在宏基因組測序中，得到的相對物種豐度的差異與測序深度、DNA提取以及測序的方法都密切相關；

　　2、兩者采用的物種注釋方法及數(shù)據(jù)庫都存在著一定差別：

　　16S采用的是將16S rDNA與Greengenes/Silva數(shù)據(jù)庫進行比對注釋，只能注釋到細菌或古菌；而宏基因組則是將預測得到的基因與NR數(shù)據(jù)庫或特定的marker基因數(shù)據(jù)庫比對從而進行注釋，宏基因組注釋得到的物種信息更為全面，不僅包括細菌，還包括真菌、古菌以及病毒等。

　　很多文章在物種豐度水平更多使用16s的結果，而功能注釋則使用metagenome的結果。

　　宏基因組分析策略

　　宏基因組研究主要有兩種分析策略：

　　第一種有參考基因組，直接將reads比對到參考基因組上進行研究。

　　第二種為無參考基因組，需要通過序列組裝、基因預測、再進行物種功能注釋，是一種基于de novo組裝的研究。

　　我們平時所做的宏基因組研究大多數(shù)是把兩種研究方法結合起來。

　　為什么基因組的組裝不完整？

　　宏基因組組裝的效果主要跟以下幾個因素有關：樣本的測序數(shù)據(jù)量，物種多樣且豐度分布不均勻等，這些因素都會造成宏基因組組裝比細菌等單物種的組裝更加困難，這也是目前宏基因組研究中有待突破的重點。

　　隨著測序讀長不斷增加，測序質(zhì)量的不斷提高，將三代測序數(shù)據(jù)應用到宏基因組中將是未來研究的一大方向。

　　為什么不把所有樣本合并在一起進行組裝？

　　不同樣本中物種豐度差異很大，如果把所有樣本都混合在一起，對服務器的要求很高（例如需要大內(nèi)存服務器），將會大大增加數(shù)據(jù)的復雜度，組裝效果可能會更差。

　　宏基因組測序的測序技術選擇，二代測序還是三代？

　　二代測序的短讀長限制了contig的長度，基因之間的物理位置很難準確確定。

　　隨著三代測序技術成熟，基于三代PacBio的SMRT單分子測序和Nanopore長讀長特性可以跨越大的重復區(qū)域，使得宏基因組拼接組裝效果有了很大程度的提高，能組裝出更長的contig，增加了完整基因的數(shù)量，可以對復雜的重復區(qū)域進行分析。

　　宏基因組測序深度的選擇？

　　1）由于受到測序深度及測序成本的影響，在現(xiàn)在的宏基因組文章中，測序數(shù)據(jù)量一般選擇大于5G，可以測出樣品中絕大多數(shù)的微生物，但是對于一些低豐度的物種，因為測序深度的原因，確實很有可能會組裝不出來；

　　2）在宏基因組分析中，也一般多關注的是較高豐度物種的組成情況，如果要對低豐度物種進行特殊分析，一般需要加大測序數(shù)據(jù)量，或者在前期提取過程中經(jīng)過一些特殊的處理，盡可能的富集出多的低豐度物種，再進行測序分析。

　　關于宏基因組的測序量，這個要根據(jù)樣品的復雜度來看：

　　如果diversity比較低，例如像人腸道之類的， 每個樣品測個5-10Gb的數(shù)據(jù)就可以了，對于復雜樣品，例如土壤，致使測序幾十到幾百G，也可能也會深度不足。

　　如果不知道樣品的diversity，一般建議先對樣品做一個 16S rRNA測序的survey，來看一下你的樣品里面大致有多少種微生物。

　　抗性基因和毒力基因分析所用數(shù)據(jù)庫是什么？

　　隨著人們對抗性基因相關研究的廣泛關注，我們宏基因組的標準分析中推出了抗性基因的相關分析。并且，由于自2009年ARDB數(shù)據(jù)庫再無更新，因此我們目前所用的抗性基因數(shù)據(jù)庫為CARD數(shù)據(jù)庫。

　　VFDB（Virulence Factors Database）數(shù)據(jù)庫用來注釋毒力基因。

　　環(huán)境因子與微生物組的關系

　　RDA 或者CCA定義是基于對應分析發(fā)展而來的一種排序方法，將對應分析與多元回歸分析相結合。此分析是主要用來反映菌群與環(huán)境因子之間關系。RDA是基于線性模型，CCA 是基于單峰模型。分析可以檢測環(huán)境因子、樣品、菌群三者之間的關系或者兩兩之間的關系。

　　注：圖中數(shù)字表示樣品名，不同顏色或形狀表示不同環(huán)境或條件下的樣本組；箭頭表示環(huán)境因子；圖中藍色上三角表示不同屬的細菌；物種與環(huán)境因子之間的夾角代表物種與環(huán)境因子間的正、負相關關系（銳角：正相關；鈍角：負相關；直角：無相關性）；由不同的樣本向各環(huán)境因子做垂線，投影點越相近說明樣本間該環(huán)境因子屬性值越相似，即環(huán)境因子對樣品的影響程度相當。

　　如何進行腸型分析？

　　2011年腸型（Enterotypes）概念第一次被提出，腸型被分為3種Bacteroides (enterotype 1), Prevotella (enterotype 2) and Ruminococcus (enterotype 3)，文中指出：

　　1）在不同的腸道樣本中，均有明顯的類別存在；

　　2）3類腸型中的優(yōu)勢菌，在不同樣本中是不同的。

　　影響腸型結果的各種因素：

　　1）不同聚類方法、距離計算方式最終得到的腸型會有不同

　　2）不同OTU的層級得到的腸型不同

　　3）不同測序區(qū)域得到的腸型不同

　　4）OTU的計算策略得到的腸型不同

　　5）測序策略得到的腸型不同

　　腸型分類器的在線使用

　　腸型在線分類器主頁: http:///基于MetaHIT數(shù)據(jù)集，建立了一個基于屬水平的在線分類器(基于HMP1和中國二型糖尿病 )

　　宏基因組測序能否鑒定到菌株水平？

　　CAG (co-abundance gene group)

　　CAG方法是2014年在《自然?生物技術》雜志上報道出來的人體腸道宏基因組數(shù)據(jù)分析方法。這種方法的思想和MGS的思想有一定相同之處，都利用了這樣一種思想或理論：在同一個株或種中的基因，這些基因的豐度在群體中具有豐度一致性。為了解決使用MGS方法時僅關注基因標記的局限性，canopy算法被用來鑒定共豐富基因簇(CAG)，CAG分析基于基因豐度鑒定物種，并且每個CAG可以被看作是部分微生物或完整微生物。

　　當基因豐度表包含足夠數(shù)量的樣本時，可以應用canopy算法來識別CAG。canopy算法使用Pearson相關系數(shù)和Spearman相關系數(shù)作為閾值，在對只有1個基因的canopies進行第一輪聚類和過濾之后，根據(jù)簇的平均豐度，使用Pearson相關系數(shù)作為閾值再次被應用到canopy算法中，第二輪簇可以包含重疊基因。

　　因此，對于出現(xiàn)在多于1個簇中的基因，基因及其相關簇之間的距離能夠被確定，對于每個重疊基因，選擇最近的簇。最后，選擇含有超過700個基因的簇最可能是含有細菌基因組部分的CAG。

　　基于CAG分析的結果，可以在人腸道微生物群中鑒定出可能顯著影響宿主的未分類或難以理解的微生物(也稱為生物暗物質(zhì))。在過去的研究中，生物暗物質(zhì)一般被忽略，因此，許多疾病相關的腸道微生物可能未被檢測到。CAG中的重疊群和基因可以產(chǎn)生無法從當前的細菌基因組數(shù)據(jù)庫獲得的大量新信息。

　　兩款鑒定到菌株的軟件

　　PanPhlAn：PanPhlAn - Strain-level resolution in Metagenomics（一種基于菌株特異性泛基因組學的分析方法）

　　ConStrains: ConStrains identifies microbial strains in metagenomic datasets

　　相關名詞解釋

　　De novo測序從頭測序，無需任何參考序列，直接對一個物種進行測序，然后進行拼接、組裝成該物種的基因組序列圖譜。

　　建庫在基因組隨機打斷的片段上機前，為保證在測序時有足夠的數(shù)據(jù)強度支持，所進行的基于PCR反應的片段擴增，區(qū)別于通常說的基于Fosmid/BAC等載體的建立文庫。

　　Paired-End測序雙末端測序，對插入片段兩端進行測序，產(chǎn)生具有Paired-End關系的reads。

　　測序策略sequences strategyPE100：即Paired-End (100,100)，采用雙末端測序法對插入片段兩端進行讀長為100 bp的高通量測序。

　　讀長測序儀所能獲取的實際長度，即reads的長度，例如：90bp、125bp、300bp。

　　Raw data（reads）即下機數(shù)據(jù)，原始的reads。

　　Clean data（reads）即下機數(shù)據(jù)經(jīng)過去接頭污染、過濾低質(zhì)量reads、去duplication（針對大片段數(shù)據(jù)）等之后，實際用于組裝及分析的數(shù)據(jù)。

　　Contig序列由來源于同一基因組，具有overlapping關系的reads拼接而成的片段。

　　Scaffold序列又稱super Contig，通過reads的Paired-End關系將Contig連接，并用N補充其中的內(nèi)洞而形成。

　　N50/N90將組裝所得片段（Scaffold/Contig）按照從長到短排序并累加求和，累加值達到基因組總長度一半時的片段長度即是該組裝結果的N50值，通常用來衡量組裝情況；N90與之類似，即累加長度為基因組總長90%時，該片段的長度。

　　基因集（Gene Catalog）由所有樣本中檢測到的基因構成的集合。

　　非冗余基因集（Non-redundant Gene Catalog）將所有樣本的基因以一定的閾值進行聚類，每個類別取最長的基因作為代表序列。

　　物種分類學注釋：將基因集序列與參考數(shù)據(jù)庫比對，并根據(jù)分類學譜系系統(tǒng)得到該序列的物種分類地位，從而分析得到整個基因集的物種信息；

　　COG功能注釋：將基因集序列與COG參考數(shù)據(jù)庫比對，獲得基因集的同源聚類蛋白群的相應功能注釋信息；

　　KEGG功能注釋：將基因集序列與KEGG的基因組數(shù)據(jù)庫比對，獲得基因集在代謝、酶、反應和功能模塊方面的注釋信息。

　　CAZy數(shù)據(jù)庫注釋：將基因集序列與特殊數(shù)據(jù)庫比對，如 CAZy數(shù)據(jù)庫（Carbohydrate-Active enZYmes Database，碳水化合物活性酶）

　　CARD數(shù)據(jù)庫（Comprehensive Antibiotic Resistance Database，抗生素抗性基因）， 獲得基因集特殊的功能信息