MetaWRAP分箱流程实战和结果解读

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.gzwget 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.gzwget 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/*; dob=${i#*/}mv ${i}/final_pure_reads_1.fastq CLEAN_READS/${b}_1.fastqmv ${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

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读