文章目錄

• MetaWRAP—a flexible pipeline for genome-resolved metagenomic data analysis
• 分析實戰
• 0.下載腸道巨集基因組資料
• 1.read_qc質控和去宿主
• 2. 組裝
• 3. kraken物種註釋(可選)
• 4. 三種方法Bin
• 5.Bin提純
• 6. Blobology視覺化bin
• 7. Bin定量
• 8.重組裝
• 9.bin物種註釋
• 10.Bin功能註釋
• 我們還能做什麼
• 常見錯誤
• bin提純報錯
• Reference
• 猜你喜歡
• 寫在後面

MetaWRAP—a flexible pipeline for genome-resolved metagenomic data analysis

MetaWRAP——靈活的單基因組精度巨集基因組分析流程

關於巨集基因組Binning，有無數的軟體和資料庫，大家分析費時費力，結果也差別很大。現在有了MetaWRAP，一個軟體就夠了，整合3個主流分箱工具、並增長了提純和重組裝環節，保證結果最約。

關於軟體評價，文章導讀，請閱讀

• Microbiome：巨集基因組分箱流程MetaWRAP 簡介

關於軟體的安裝和資料庫的部署，請閱讀

• Microbiome：巨集基因組分箱流程MetaWRAP 安裝和資料庫部署

今天帶來本軟體最後一節，實戰和結果解讀。

分析實戰

https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md

0.下載腸道巨集基因組資料

# 建立工作目錄
cd ~
mkdir metawrap && cd metawrap

## 下載3個樣品，共6G，解壓後15G；包括7G的序列，我下載了12小時
wget ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011347/ERR011347_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011347/ERR011347_2.fastq.gz

wget ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011348/ERR011348_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011348/ERR011348_2.fastq.gz

wget ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011349/ERR011349_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR011349/ERR011349_2.fastq.gz

# 解壓
gunzip *.gz

# 移動至子目錄
mkdir RAW_READS
mv *fastq RAW_READS
 

1.read_qc質控和去宿主

預設使用人類基因組的hg38的bmt索引去宿主，需要在軟體安裝和資料庫配置中提前佈置好。對於環境樣本，也可以不去宿主，使用 --skip-bmtagger 跳過。

我們對每個樣本進行處理

mkdir READ_QC
# 24執行緒程質控5GB資料，大約29分鐘(m)單個樣本，但軟體好像不支援多執行緒
time metawrap read_qc -1 RAW_READS/ERR011347_1.fastq -2 RAW_READS/ERR011347_2.fastq -t 24 -o READ_QC/ERR011347
metawrap read_qc -1 RAW_READS/ERR011348_1.fastq -2 RAW_READS/ERR011348_2.fastq -t 24 -o READ_QC/ERR011348
metawrap read_qc -1 RAW_READS/ERR011349_1.fastq -2 RAW_READS/ERR011349_2.fastq -t 24 -o READ_QC/ERR011349

我們可以檢視每個樣品目錄中的報告，來比較質控前後的變化，如READ_QC/ERR011347/pre-QC_report/ERR011347_2_fastqc.html包含質控前的質量評估結果。

READ_QC/ERR011347/post-QC_report/final_pure_reads_2_fastqc.html包含質控前的質量評估結果。

如果確定質量合格(質量不符合你自己要求，可檢視幫助調整引數重新質控)，可移動結果到新目錄，質控後15G變為9G。

mkdir CLEAN_READS
for i in READ_QC/*; do 
	b=${i#*/}
	mv ${i}/final_pure_reads_1.fastq CLEAN_READS/${b}_1.fastq
	mv ${i}/final_pure_reads_2.fastq CLEAN_READS/${b}_2.fastq
done

2. 組裝

檔案合併，方便拼接簡化輸入檔案、組裝獲得統一的參考contig，如果檔案過大，也可能需要分批處理。

cat CLEAN_READS/ERR*_1.fastq > CLEAN_READS/ALL_READS_1.fastq
cat CLEAN_READS/ERR*_2.fastq > CLEAN_READS/ALL_READS_2.fastq

組裝，獲得疊連群(contigs)，是Binning的基本單元

# 執行1h，仍然報錯，去掉--metaspades後8m完成
time metawrap assembly -1 CLEAN_READS/ALL_READS_1.fastq -2 CLEAN_READS/ALL_READS_2.fastq -m 200 -t 96  -o ASSEMBLY # --metaspades

組裝結果為ASSEMBLY/final_assembly.fasta，統計報告見ASSEMBLY/assembly_report.html

檢視最長的3個contig資訊

grep ">" ASSEMBLY/final_assembly.fasta | head -n3

>NODE_1_length_196124_cov_2.427049
>NODE_2_length_176373_cov_3.889994
>NODE_3_length_163601_cov_3.070200

看到最長的contig將近200kb，Top 3大於150kb，我們只使用了測試的7G資料，通常使用30 - 100 GB，拼接結果會更好。

3. kraken物種註釋(可選)

kraken在reads和contig層面進行物種註釋。當然kraken這麼全能的工具，還可以在基因和bin層面

# 7G演示資料，8執行緒耗時3m17s
metawrap kraken -o KRAKEN -t 96 -s 1000000 CLEAN_READS/ERR*fastq ASSEMBLY/final_assembly.fasta

結果資料夾中有註釋結果檔案*.kraken(每個reads的註釋結果)、*.krona(註釋結果分類彙總)，和視覺化的Krona網頁圖表kronagram.html

---

如果你己獲得了contigs，可以從下面開始！！！

4. 三種方法Bin

基於contig和三種bin方法，需要拼接結果和原始序列，用時34m

metawrap binning -o INITIAL_BINNING -t 8 -a ASSEMBLY/final_assembly.fasta --metabat2 --maxbin2 --concoct CLEAN_READS/ERR*fastq

結果目錄中檔案或目錄說明

insert_sizes.txt 樣本量統計和估計的建庫時文庫片段大小

concoct_bins maxbin2_bins metabat2_bins 三個目錄為三種bin的結果

work_files有三種bin分析所需要的檔案，如不同格式的bin覆蓋度或丰度資訊。

concoct, metabat2和maxbin2三個軟體分別獲得了47, 29 和20個bins，但我們並不知道它們質量如何，可以新增--run-checkm引數獲得質量評估，但我們更喜歡下面獨立分析！

沒有比較就沒有傷害！只是為了突出本流程的優秀。

5.Bin提純

三種主流bin結果各有優缺點，我們將對結果進行整合和優化，以獲得更好的結果

如果你有超過3種結果，也可以合併，如5種結果，先合併1，2，3；再合併4，5；最後將兩類合併。

推薦使用完整度70，汙染率5的閾值；但本演示測序資料較小，僅使用50和10級別的閾值，保證有足夠的結果用於演示和評估。要求越高，bin越少，請根據個人需要調整。

主要引數有-o輸出目錄；-t執行緒數；-A/B/C三種Bin結果；-c完整度閾值；-x汙染率閾值；如下8執行緒，耗時67m

metawrap bin_refinement -o BIN_REFINEMENT -t 8 -A INITIAL_BINNING/metabat2_bins/ -B INITIAL_BINNING/maxbin2_bins/ -C INITIAL_BINNING/concoct_bins/ -c 50 -x 10

結果目錄中有三個原始bin的結果與統計。metaWRAP目錄包含最終結果。如果想檢視中間檔案見Binning_refiner目錄。

.stat檔案包含每個bin的統計：完整性、汙染率、GC含量、物種、N50、大小和來源。

如cat BIN_REFINEMENT/metaWRAP_bins.stats

bin     completeness    contamination   GC      lineage N50     size    binner
bin.1   98.92   0.866   0.401   Clostridiales   13297   2373299 binsBC
bin.8   93.73   0.884   0.312   Euryarchaeota   5058    1485694 binsC
bin.2   85.57   0.223   0.370   Clostridiales   5890    2030996 binsA

提純的結果如何看，詳見BIN_REFINEMENT/figures/目錄中的圖片，有eps向量圖和Png點陣圖供檢視和修改發表使用，作者太貼心了！

6. Blobology視覺化bin

我們使用Blobology模組，將Contig的GC含量和豐度進行散點圖視覺化，並按物種或Bin進行著色，來觀察結果。

視覺化Bin的contig丰度和GC含量，耗時31m

metawrap blobology -a ASSEMBLY/final_assembly.fasta -t 8 -o BLOBOLOGY --bins BIN_REFINEMENT/metaWRAP_bins CLEAN_READS/ERR*fastq

存在如下圖所 需的原始資料.binned.blobplot.檔案，包括GC含量、丰度、平均丰度等，可使用ggplot2方便視覺化如下圖，難點是顏色分配。

參考指令碼見程式目中
/conda2/envs/metawrap/bin/metawrap-scripts/blobology/makeblobplot_with_bins.R final_assembly.binned.blobplot 1 taxlevel_phylum有三個引數，指令碼較老，需要自己修改一下。

按門水平著色的GC含量vs丰度散點圖

按Bin著色的GC含量vs丰度散點圖

7. Bin定量

我們想知道這些單菌基因組草圖(bin)在每個樣品中的丰度。quant_bin模組幫你實現，它也依賴salmon來實現定量，估計每一個scaffold的丰度，以及bin的平均丰度。

實際時間1m39s，計算時間45m6s，使用呼叫了我的所有執行緒。此處可設定執行緒，但salmon仍儘可能多呼叫資源。

metawrap quant_bins -b BIN_REFINEMENT/metaWRAP_bins -t 8 -o QUANT_BINS -a ASSEMBLY/final_assembly.fasta CLEAN_READS/ERR*fastq

結果包括bin丰度熱圖QUANT_BINS/genome_abundance_heatmap.png。

如果想自己畫圖，原始資料位於QUANT_BINS/abundance_table.tab

Genomic bins	ERR011349	ERR011348	ERR011347
bin.9	0.113912116828	35.851964987	39.8440491514
bin.10	0.273774684856	9.52869077293	39.988244574

8.重組裝

提純的bin還可以通過再組裝進一步改善結果。基於原始reads對結果優化，只有結果更優的情況，才對結果進行更新。

先呼叫bwa比對原始序列到bin，使用SPAdes不同引數下(寬鬆和嚴謹)精細組裝，去除小於1500bp的疊連群，checkM評估，結果統計，繪圖。

用時41m, 執行緒總用時178m

metawrap reassemble_bins -o BIN_REASSEMBLY -1 CLEAN_READS/ALL_READS_1.fastq -2 CLEAN_READS/ALL_READS_2.fastq -t 8 -m 800 -c 50 -x 10 -b BIN_REFINEMENT/metaWRAP_bins

結果統計見BIN_REASSEMBLY/reassembled_bins.stats，發現3個bin通過嚴格模式下得到優化，6個通過寬鬆模式得到改進，4個沒有變化。

bin	completeness	contamination	GC	lineage	N50	size	binner
bin.10.orig	65.78	0.0	0.263	Bacteria	3045	1159966	NA
bin.7.strict	54.94	0.671	0.501	Clostridiales	3947	1474089	NA
bin.4.permissive	99.32	1.342	0.408	Clostridiales	72135	2088821	NA

改進的情況如何呢？要檢視結果BIN_REASSEMBLY/reassembly_results.png，比對重組裝前後的變化，N50、這完整度和汙染率均有改進。

BIN_REASSEMBLY/reassembled_bins.png展示CheckM對bin評估結果的視覺化。

綠色為單拷貝基因，灰色為缺失基因；藍色梯度代表不同的雜合率；紅色梯度代表汙染率。

9.bin物種註釋

其實Bin提純和重組裝中，在checkM的stat檔案中，就有物種的註釋結果。但軟體和資料庫都不完善。

基於NCBI_nt和NCBI_tax資料庫，我們使用 Taxator-tk 進行每條contig物種註釋，再估計bin整體的物種。註釋結果的準確性由參考資料庫決定(如參考庫中有錯誤，我們無法識別，只能將錯就錯)。

此步8執行緒用時12分鐘。

# real 12m, user 84m
metawrap classify_bins -b BIN_REASSEMBLY/reassembled_bins -o BIN_CLASSIFICATION -t 8

註釋結果見BIN_CLASSIFICATION/bin_taxonomy.tab

bin.1.orig.fa	Bacteria;Firmicutes;Clostridia;Clostridiales
bin.5.orig.fa	Archaea;Euryarchaeota;Methanobacteria;Methanobacteriales;Methanobacteriaceae;Methanobrevibacter;Methanobrevibacter smithii
bin.11.orig.fa	Bacteria;Bacteroidetes;Bacteroidia;Bacteroidales;Bacteroidaceae;Bacteroides

注意物種註釋層級不同，因為有些物種在資料庫中沒有近親。

10.Bin功能註釋

此步基於PROKKA進行基因註釋。

# 4m5s
metaWRAP annotate_bins -o FUNCT_ANNOT -t 8 -b BIN_REASSEMBLY/reassembled_bins/

每個bin基因註釋的gff檔案見 FUNCT_ANNOT/bin_funct_annotations

head FUNCT_ANNOT/bin_funct_annotations/bin.1.orig.gff

NODE_75_length_31799_cov_0.983871	Prodigal:2.6	CDS	2866	3645	.	-	0	ID=HMOHEJHL_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=HMOHEJHL_00001;product=hypothetical protein
NODE_75_length_31799_cov_0.983871	Prodigal:2.6	CDS	3642	4478	.	-	0	ID=HMOHEJHL_00002;inference=ab initio prediction:Prodigal:2.6;locus_tag=HMOHEJHL_00002;product=hypothetical protein
NODE_75_length_31799_cov_0.983871	Prodigal:2.6	CDS	4606	5859	.	-	0	ID=HMOHEJHL_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=HMOHEJHL_00003;product=hypothetical protein

基因蛋白和核酸序列見FUNCT_ANNOT/bin_translated_genes和 FUNCT_ANNOT/bin_untranslated_genes目錄。PROKKA輸出結果見 FUNCT_ANNOT/prokka_out.

我們還能做什麼

現在的Binning分析己非常全面了。接下來我們還能做些什麼呢？Binning的結果主要是細菌基因組的草圖，相當於單菌基因組的研究領域，我們可以結合功能註釋、進化和泛基因組學研究手段進行研究。

比如Binning分析的經典高分文章，2個樣本發NC。
請參考如下相關文章：

• Nat Commun：巨集基因組學提示曙古菌門的代謝和進化(中大李文均組)

比如提取同源基因進化分析

• 進化樹 一文讀懂 Evolview基礎 Evolview進階 iTOL美化
• iTOL快速繪製顏值最高的進化樹！
• 線上RaxML構建系統發育樹

比如在Bin中代謝物和基因簇的挖掘，

• antiSMASH：微生物次生代謝物基因簇預測

常見錯誤

bin提純報錯

[Errno 13] Permission denied: '/conda2/envs/metawrap/lib/python2.7/site-packages/checkm/DATA_CONFIG’

沒有許可權，請新增此檔案許可權，checkm才能寫入正確的資料庫位置

[Error] No bins found. Check the extension (-x) used to identify bins.

可能原因：資料太小，沒有bin，需要加大測序量，或輸入資料量

Reference

主頁和軟體安裝教程：https://github.com/bxlab/metaWRAP

資料庫佈署：https://github.com/bxlab/metaWRAP/blob/master/installation/database_installation.md

使用教程：https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md

猜你喜歡

• 10000+: 菌群分析
  寶寶與貓狗 提DNA發Nature 實驗分析誰對結果影響大 Cell微生物專刊 腸道指揮大腦
• 系列教程：微生物組入門 Biostar 微生物組 巨集基因組
• 專業技能：生信寶典 學術圖表 高分文章 不可或缺的人
• 一文讀懂：巨集基因組 寄生蟲益處 進化樹
• 必備技能：提問 搜尋 Endnote
• 文獻閱讀 熱心腸 SemanticScholar Geenmedical
• 擴增子分析：圖表解讀 分析流程 統計繪圖
• 16S功能預測 PICRUSt FAPROTAX Bugbase Tax4Fun
• 線上工具：16S預測培養基 生信繪圖
• 科研經驗：雲筆記 雲協作 公眾號
• 程式設計模板: Shell R Perl
• 生物科普: 腸道細菌 人體上的生命 生命大躍進 細胞暗戰 人體奧祕

寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了“巨集基因組”專業討論群，目前己有國內外2300+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註“姓名-單位-研究方向-職稱/年級”。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。

學習擴增子、巨集基因組科研思路和分析實戰，關注“巨集基因組”

點選閱讀原文，跳轉最新文章目錄閱讀
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA