腸道菌群:宏基因組測序分析流程解讀(上)

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

- 概述

- 宏基因組測序實驗流程

- 宏基因組分析流程

- 文章分析思路

- 常見問題

- 相關名詞解釋

- 參考文獻

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

概述

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

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

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

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

高通量測序平臺

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

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

細菌16S rDNA基因

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

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

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

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

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

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

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

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

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

宏基因組測序實驗流程

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

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

實驗流程圖DNA樣品檢測

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

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

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

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

文庫構建及質檢

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

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

上機測序

建庫質檢合格後，把不同文庫按照有效濃度及目標下機數據量的需求，混合後進行Illumina測序。

宏基因組分析流程

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

信息分析流程圖測序數據預處理

•常規數據質控工具

ØFASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/commandline.html)

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

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

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

•去除宿主汙染

ØSoapAligner (http://soap.genomics.org.cn/soapaligner.html)

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

ØBowtie (http://sourceforge.net/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 測序平臺測序獲得的原始數據(Raw Data)存在一定比例低質量數據，裡面含有帶接頭的、重複的，以及測序質量很低的reads，這些 reads 會影響組裝和後續分析，為了保證後續分析的結果準確可靠，需要對原始的測序數據進行預處理，獲取用於後續分析的有效數據(Clean Data)。

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

例如具體處理步驟如下：

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

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

3) 如果樣品存在宿主汙染，需與宿主資料庫進行比對，過濾掉可能來源於宿主（可以採用 SoapAligner 軟體，參數設置: identity ≥ 90%, -l 30, -v 7, -M 4,-m 200,-x 400）的 reads；

數據處理信息統計表

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

物種注釋metaphlan2

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

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

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

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

精確的分類群分配

準確估計物種的相對豐度

種水平精度

株鑑定與追蹤

超快的分析速度

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)

實驗室主頁: http://segatalab.cibio.unitn.it

看看軟體主要能幹什麼

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

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

使用實例定量物種譜

一條命令就可以從原始數據得到物種的相對豐度！！

最常用的是使用雙端壓縮fastq文件，參數--nproc調用多個線程，並輸出比較結果

metaphlan2.py --bowtie2out metagenome.bowtie2.bz2 --nproc 20 --input_type fastq <(zcat metagenome_1.fq.gz metagenome_2.fq.gz) > profiled_metagenome1.txt

例如人類腸道的數據，一般幾個小時內就出結果了。

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

#SampleID Metaphlan2_Analysis

k__Bacteria 100.0

k__Bacteria|p__Bacteroidetes 58.71768

k__Bacteria|p__Firmicutes 27.73866

k__Bacteria|p__Verrucomicrobia 9.51729

k__Bacteria|p__Proteobacteria 3.20978

k__Bacteria|p__Actinobacteria 0.81659

k__Bacteria|p__Bacteroidetes|c__Bacteroidia 58.71768

k__Bacteria|p__Firmicutes|c__Bacilli 17.31386

k__Bacteria|p__Firmicutes|c__Clostridia 10.26799

k__Bacteria|p__Verrucomicrobia|c__Verrucomicrobiae 9.51729

k__Bacteria|p__Proteobacteria|c__Deltaproteobacteria 2.43682

k__Bacteria|p__Actinobacteria|c__Actinobacteria 0.81659

k__Bacteria|p__Proteobacteria|c__Betaproteobacteria 0.62445

k__Bacteria|p__Firmicutes|c__Erysipelotrichia 0.15682

k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria 0.14851

k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales 58.71768

k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales 17.31386

批量處理多個樣品

假如我們有1-20一共20個樣本，調用for循環並轉後臺並行處理

for i in `seq 1 20`;do

metaphlan2.py --nproc 5 --input_type fastq <(zcat Sample_${i}.1.fq.gz Sample_${i}.2.fq.gz) > metaphlan2_${i}.txt &

done

樣品物種譜合併

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

merge_metaphlan_tables.py metaphlan2*.txt | sed 's/metaphlan2_//g' > merged.txt

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

種相對豐度概況

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

門水平菌落結構柱狀圖基於GraPhlAn的分類學組成信息可視化

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

功能注釋humann2

Species-level functional profiling of metagenomes and metatranscriptomes

Nature Methods, Article, 2018-10-30

DOI: http://dx.doi.org/10.1038/s41592-018-0176-y

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

通訊作者：Curtis Huttenhower

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

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

humann2 --threads 20 --input test.fastq --output humann2_out/

humann2_renorm_table --input abundance.tsv --output abundance_relab.tsv --units relab

humann2_join_tables --input humann2_out --output humann2_pathabundance.tsv --file_name pathabundance_relab

merge_metaphlan_tables.py *tsv > metaphlan2_merged.tsv

Metagenome 組裝和注釋Metagenome 組裝

兩種常見策略：

1 基於序列overlap關係進行拼接

2 基於de Bruijn圖進行組裝

由於現階段的主流測序方法是二代短片段測序，序列短而且數目龐大，如果利用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.

如下是三種主流軟體分析，運行所消耗時間、內存比較:

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

目前MEGAHIT在現有組裝軟體中，資源消耗基本上是最低的，因此很適合宏基因組中的複雜環境樣品。

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

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

從各樣品質控後的Clean reads出發，組裝主要分為3步：單樣本組裝，多樣本組裝結果合併和豐度過濾。

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

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

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

4) 將各樣品質控後的CleanData採用SoapAligner軟體比對至各樣品組裝後的Scaftigs上，獲取未被利用上的PE reads；比對參數：-u,-2,-m200

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

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

7) 對於單樣品和混合組裝生成的Scaftigs，過濾掉500bp以下的片段，並進行統計分析和後續基因預測；

各樣品組裝結果基本信息統計

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

樣品的scaftigs長度分布統計如下：

單個樣品scaftigs長度分布統計

各樣品scaftigs數目統計：

所有樣品scaftigs數目統計圖組裝結果質量評估工具 MetaQUAST

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

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

2) 對各樣品及混合組裝的 ORF 預測結果，採用CD-HIT軟體進行去冗餘，以獲得非冗餘的初始 genecatalogue（此處，操作上，將非冗餘的連續基因編碼的核酸序列稱之為 genes），默認以 identity 95%,coverage 90% 進行聚類，並選取最長的序列為代表性序列；採用參數：-c 0.95, -G 0, -aS 0.9, -g 1, -d 0

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

4）過濾掉在各個樣品中支持 reads 數目<=2 的基因，獲得最終用於後續分析的 genecatalogue（Unigenes）

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

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

Gene catalogue 基本信息

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

gene catalogue 長度分布統計core-pan 基因分析

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

core-pan 基因稀釋曲線組間基因數目差異分析

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

組間基因數目差異小提琴圖

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

基因數目韋恩圖分析

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

基因數目韋恩圖（花瓣圖）分析

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

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

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

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

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

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

PCA圖

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

2) NMDS分析

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

NMDS圖組間差異物種的LEfSe分析

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

常用功能資料庫注釋

功能注釋主要是基於功能同源的序列往往具有序列相似性的原理，將去冗餘後的基因蛋白序列，與不同的蛋白功能資料庫進行序列比對， 然後用比對到的序列的功能作為目標序列的功能。

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

目前常用的功能資料庫主要有：

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

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

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

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

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

功能注釋基本步驟

1）使用DIAMOND軟體將 Unigenes 與各功能資料庫進行比對（blastp，evalue ≤ 1e-5）；

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

3）從比對結果出發，統計不同功能層級的相對豐度（各功能層級的相對豐度等於注釋為該功能層級的基因的相對豐度之和），其中，KEGG 資料庫劃分為 5 個層級，eggNOG 資料庫劃分為 3 個層級，CAZy 資料庫劃分為 3 個層級，各資料庫的詳細劃分層級如下所示：

資料庫名稱劃分層級該層級的描述KEGGLevel 1KEGG 代謝通路第一層級 6 大代謝通路；KEGGLevel 2KEGG 代謝通路第二層級 43 種子 pathway；KEGGLevel 3KEGG pathway id（例：ko00010）；KEGGkoKEGG ortholog group (例：K00010)；eggNOGLevel 124大功能類eggNOGLevel 2ortholog group descriptioneggNOGogortholog group IDCAZyLevel 16大功能類CAZyLevel 2CAZy familyCAZyLevel 3EC number

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

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

注釋基因數目統計

從 Unigenes 注釋結果出發，繪製各個資料庫的注釋基因數目統計圖，展示結果如下圖所示：

各資料庫注釋基因數目統計圖

說明：從上至下依次為CAZy，eggNOG 的 Unigenes 注釋數目統計圖。橫坐標軸是各資料庫中 level1 各功能類的代碼，代碼的解釋見對應的圖例說明。

功能相對豐度概況

從各資料庫 level1 的相對豐度表出發，繪製出各個資料庫中，各樣品對應的在 level1 層級上的豐度統計圖。

KEGG:

Eggnog:

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

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

功能相對豐度聚類分析

根據所有樣品在各個資料庫中的功能注釋及豐度信息，選取豐度排名前 35 的功能及它們在每個樣品中的豐度信息繪製熱圖，並從功能差異層面進行聚類。

Kegg:

Eggnog:

功能豐度聚類熱圖

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

基於功能豐度的降維分析

基於不同資料庫在各個分類層級的功能豐度進行 PCA 和 NMDS 降維 分析，如果樣品的功能組成越相似，則它們在降維圖中的距離越接近。基於KEGG的KO功能豐度進行PCA和NMDS 分析的結果展示如下：

抗生素抗性基因注釋

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

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

抗性基因豐度概況

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

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

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

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

根據抗性基因類型的豐度表，繪製聚類熱圖：

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

組間抗性基因數目差異分析

註：橫軸為分組，縱軸為抗性基因類型數目

物種間相關性網絡分析

機器學習算法

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

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

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

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

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

文章分析思路

文章分析思路：整體概述——物種組成——功能組成——關聯分析——因果關係驗證

文章主要結果

樣本背景、數據評估、物種、功能描述

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

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

環境因子與微生物組的關係

菌群或單菌移植到無菌小鼠體內，驗證因果關係

常見問題16S擴增子和宏基因組分析結果為什麼存在差別？

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

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

2、兩者採用的物種注釋方法及資料庫都存在著一定差別：

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

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

宏基因組分析策略

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

關於宏基因組的測序量，這個要根據樣品的複雜度來看：

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

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

抗性基因和毒力基因分析所用資料庫是什麼？

隨著人們對抗性基因相關研究的廣泛關注，我們宏基因組的標準分析中推出了抗性基因的相關分析。並且，由於自2009年ARDB資料庫再無更新，因此我們目前所用的抗性基因資料庫為CARD資料庫。

VFDB（Virulence Factors Database）資料庫用來注釋毒力基因。

環境因子與微生物組的關係

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

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

如何進行腸型分析？

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

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

2）3類腸型中的優勢菌，在不同樣本中是不同的。

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

腸型分類器的在線使用

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

宏基因組測序能否鑑定到菌株水平？CAG (co-abundance gene group)

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

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

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

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

兩款鑑定到菌株的軟體相關名詞解釋

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

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

Paired-End測序 雙末端測序，對插入片段兩端進行測序，產生具有Paired-End關係的reads。

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

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

Raw data（reads）即下機數據，原始的reads。

Clean data（reads）即下機數據經過去接頭汙染、過濾低質量reads、去duplication（針對大片段數據）等之後，實際用於組裝及分析的數據。

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

Scaffold序列 又稱super Contig，通過reads的Paired-End關係將Contig連接，並用N補充其中的內洞而形成。

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

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

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

物種分類學注釋：將基因集序列與參考資料庫比對，並根據分類學譜系系統得到該序列的物種分類地位，從而分析得到整個基因集的物種信息；

COG功能注釋：將基因集序列與COG參考資料庫比對，獲得基因集的同源聚類蛋白群的相應功能注釋信息；

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

CAZy資料庫注釋：將基因集序列與特殊資料庫比對，如 CAZy資料庫（Carbohydrate-Active enZYmes Database，碳水化合物活性酶）

CARD資料庫（Comprehensive Antibiotic Resistance Database，抗生素抗性基因）， 獲得基因集特殊的功能信息

參考文獻

[1] Luo et al.: SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. GigaScience 2012 1:18.

[2] Li J, Jia H, Cai X, et al. An integrated catalog of reference genes in the human gut microbiome[J]. Nature biotechnology, 2014, 32(8): 834-841.

[3] Oh J, Byrd A L, Deming C, et al. Biogeography and individuality shape function in the human skin metagenome[J]. Nature, 2014, 514(7520): 59-64.

[4] Zhu, Wenhan, Alexandre Lomsadze, and Mark Borodovsky. "Ab initio gene identification in metagenomic sequences." Nucleic acids research 38.12 (2010): e132-e132

[5] Duy Tin Truong, Eric A Franzosa,et al. MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nature Methods 12, 902-903 (2015)

[6 Kanehisa M, Goto S, Hattori M, Aoki-Kinoshita KF, Itoh M, Kawashima S, et al. (2006). From genomics to chemical genomics: new developments in KEGG. Nucleic Acids Res 34(Database issue): D354–7.

[7 Kanehisa, M., Goto, S., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M.; Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 42, D199–D205 (2014).

[8] Powell S, Forslund K, Szklarczyk D, et al. eggNOG v4. 0: nested orthology inference across 3686 organisms[J]. Nucleic acids research, 2013: gkt1253.

[9] Liu B, Pop M. ARDB—Antibiotic Resistance Genes Database[J]. Nucleic Acids Research, 2009, 37(Database issue):D443-7.

[10] Li W, Godzik A: Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 2006, 22(13):1658-1659.

[11] Fu L, Niu B, Zhu Z, Wu S, Li W: CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 2012, 28(23):3150-3152.

這次的內容就介紹到這了，有問題的小夥伴們歡迎留言討論！感謝您的支持，歡迎點擊「在看"和轉發！！長按下方二維碼，即可關注 「基因的生物信息學分析」。

相關閱讀：

腸道菌群：16S測序分析流程解讀