Anvi'o 宏基因组分析流程快速指南

写在前面：Anvi'o 是一个功能强大且可扩展的组学数据分析和可视化平台，可以轻松应用于泛基因组分析以及宏基因组分析。本文整理自 Anvi'o 官网的宏基因组数据分析教程。

软件流程

大致的流程简单来说就是首先把经过质控去宿主后的fq文件组装成contig，再把contig与fq文件进行比对，生成bam文件。接着使用 Anvi'o 将contig的fa文件构建contig数据库并识别 ORF 区域，使用hmmer 扫描contig数据库中的单拷贝基因，对数据库进行菌群物种分类学注释和功能注释，单个样品信息( bam )导入到数据库，再将多个数据库进行合并，同时 Binning，最后进行结果可视化。

软件安装

推荐大家使用docker进行安装，方便省事，一行代码搞定：

docker pull meren/anvio

docker 是什么：参见 我的 docker 学习笔记。

安装成功后就可以使用下面的代码进入docker容器了：

docker run --rm -it -v `pwd`:`pwd` -w `pwd` -p 8080:8080 meren/anvio:latest

-v参数可以在容器启动的时候挂载宿主机的一个目录，冒号 " : " 前面的目录是宿主机目录，后面的目录是容器内目录。

-w参数则会进入该目录。

-p参数指端口映射，会显式将一个或者一组端口从容器里绑定到宿主机上。

除了用docker安装之外，还可以用conda来安装，这里就不赘述了，参见：http://merenlab.org/2016/06/26/installation-v2/

快速使用指南

在使用 Anvi'o 之前，我们需要准备一个包含宏基因组数据contig的fa文件和将contig.fa比对到宏基因组fq数据的bam文件。

1. 数据质控和过滤宿主基因组

这里我使用了KneadData对宏基因组测序数据进行质量控制，当然这步大家也可使用其他软件。KneadData会先调用 Trimmomatic 对数据进行质控，接着调用 Bowtie2 将数据比对到宿主参考基因组上。

kneaddata -i 00_RAW/Sample_01-R1.fastq -i 00_RAW/Sample_01-R2.fastq -o 01_QC \
-db /home/shared/bowtiedb/GRCh38_PhiX --trimmomatic /usr/local/prg/Trimmomatic-0.36/ \
-t 4 --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
--bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output

关于KneadData的具体安装及使用参见：https://bitbucket.org/biobakery/kneaddata/wiki/Home

2. contig 组装拼接

这一步使用MEGAHIT进行 contig 组装拼接，为了同时处理多个样本，可以指定两个环境变量：

# py27
R1s=`ls 01_QC/*_kneaddata_paired_1.fastq | python -c 'import sys; print ",".join([x.strip() for x in sys.stdin.readlines()])'`
R2s=`ls 01_QC/*_kneaddata_paired_2.fastq | python -c 'import sys; print ",".join([x.strip() for x in sys.stdin.readlines()])'`

看看这两个变量：

echo $R1s
01_QC/Sample_01-QUALITY_PASSED_R1.fastq,01_QC/Sample_02-QUALITY_PASSED_R1.fastq
echo $R2s
01_QC/Sample_01-QUALITY_PASSED_R2.fastq,01_QC/Sample_02-QUALITY_PASSED_R2.fastq

运行MEGAHIT：

megahit -1 $R1s -2 $R2s --min-contig-len 1000  -m 0.85 -o 02_ASSEMBLY/ -t 50

关于MEGAHIT的具体安装及使用参见：https://github.com/voutcn/megahit

3. Mapping

使用Bowtie2将组装的contig比对到fq文件上：

mkdir 04_MAPPING
# build an index for our contigs
bowtie2-build 03_CONTIGS/contigs.fa 04_MAPPING/contigs

生成一个样本的bam：

bowtie2 --threads 50 -x 04_MAPPING/contigs -1 01_QC/Sample_01-QUALITY_PASSED_R1.fastq -2 01_QC/Sample_01-QUALITY_PASSED_R2.fastq -S 04_MAPPING/Sample_01.sam
samtools view -F 4 -bS 04_MAPPING/Sample_01.sam > 04_MAPPING/Sample_01-RAW.bam
anvi-init-bam 04_MAPPING/Sample_01-RAW.bam -o 04_MAPPING/Sample_01.bam
rm 04_MAPPING/Sample_01.sam 04_MAPPING/Sample_01-RAW.bam

现在，我们有了 Anvi'o 所需要的fa文件和bam文件，就可以进入 Anvi'o 环境进行下游分析啦。

4. 使用 Anvi'o 构建 contig 数据库

在命令行中输入anvi-敲两下TAB键就可以看到所有 Anvi'o 的 108 个命令，这些命令都提供了一个比较详细的文档，建议大家可以多用用-h。

进入 Anvi'o 环境：

docker run --rm -it -v `pwd`:`pwd` -w `pwd` -p 8080:8080 meren/anvio:latest

对contigs.fa文件进一步过滤和序列重命名：

mkdir 03_CONTIGS
anvi-script-reformat-fasta 02_ASSEMBLY/final.contigs.fa -o 03_CONTIGS/contigs.fa --min-len 2500 --simplify-names --report name_conversions.txt

Input ........................................: 02_ASSEMBLY/final.contigs.fa
Output .......................................: 03_CONTIGS/contigs.fa
Minimum length ...............................: 2,500
Total num contigs ............................: 74,426
Total num nucleotides ........................: 100,206,801
Contigs removed ..............................: 68940 (92.63% of all)
Nucleotides removed ..........................: 48177561 (48.08% of all)
Deflines simplified ..........................: True

使用 Anvi'o 构建 contig 数据库，同时也会调用 Prodigal 来预测 ORF：

anvi-gen-contigs-database -f contigs.fa -o contigs.db -n 'An example contigs datbase'

Input FASTA file .............................: /data1/zwbao/metagenome/19.4.26/anvio_tmp/03_CONTIGS/contigs.fa
Name .........................................: An example contigs datbase
Description ..................................: No description is given
Split Length .................................: 20,000
K-mer size ...................................: 4
Skip gene calling? ...........................: False
External gene calls provided? ................: None
Ignoring internal stop codons? ...............: False
Splitting pays attention to gene calls? ......: True

Finding ORFs in contigs
===============================================
Genes ........................................: /tmp/tmpjoviq0ti/contigs.genes
Amino acid sequences .........................: /tmp/tmpjoviq0ti/contigs.amino_acid_sequences
Log file .....................................: /tmp/tmpjoviq0ti/00_log.txt
Result .......................................: Prodigal (v2.6.3) has identified 52431 genes.

Contigs with at least one gene call ..........: 5486 of 5486 (100.0%)
Contigs database .............................: A new database, contigs.db, has been created.
Number of contigs ............................: 5,486
Number of splits .............................: 6,132
Total number of nucleotides ..................: 52,029,240
Gene calling step skipped ....................: False
Splits broke genes (non-mindful mode) ........: False
Desired split length (what the user wanted) ..: 20,000
Average split length (wnat anvi'o gave back) .: 22,313

5. 使用 hmmer 识别 contig 数据库中的单拷贝基因

anvi-run-hmms -c contigs.db -T 50

6. 物种分类学注释和功能注释

anvi-run-ncbi-cogs

Anvi'o 提供了一个非常简便的方法来进行 COG 功能注释：

# 下载数据库
anvi-setup-ncbi-cogs -T 50 --just-do-it
# 注释
anvi-run-ncbi-cogs -c contigs.db  -T 50

anvi-import-functions

Anvi'o 还可以导入其他软件生成的基因功能注释，例如 eggnog-mapper，Prokka 等。详情参见：http://merenlab.org/2016/06/18/importing-functions/

anvi-import-taxonomy

Anvi'o 通过导入其他软件生成的物种分类学注释来实现菌群注释的功能，例如 Kaiju，Centrifuge 等。详情参见：http://merenlab.org/2016/06/18/importing-taxonomy/

7. 将单个样品信息导入 contig 数据库

首先需要将 bam文件排序并建立索引，Anvi'o 也十分贴心地提供了一行命令：

anvi-init-bam SAMPLE-01-RAW.bam -o SAMPLE-01.bam

如何对多个样本批量操作呢？可以创建一个文本文件SAMPLE_IDs，包含所有的样本 ID，再用for循环读取即可。

$ cat SAMPLE_IDs
SAMPLE-01
SAMPLE-02
SAMPLE-03

for sample in `cat SAMPLE_IDs`; do anvi-init-bam $sample-RAW.bam -o $sample.bam; done

Anvi'o 的 profile 数据库存储了有关样本 contig 的特定信息。使用anvi-profile命令对样本的bam文件进行概要分析创建 profile 数据库，该概要文件会根据比对结果报告每个样本中的每个 contig 的属性。它是这个宏基因组工作流程中最关键（也是最复杂，计算要求最高）的步骤之一。

anvi-profile -i SAMPLE-01.bam -c contigs.db -T 50

8. 合并 profile 数据库

anvi-merge SAMPLE-01/PROFILE.db SAMPLE-02/PROFILE.db SAMPLE-03/PROFILE.db -o SAMPLES-MERGED -c contigs.db

或直接使用下面这个命令：

anvi-merge */PROFILE.db -o SAMPLES-MERGED -c contigs.db

9. 结果可视化

Anvi'o 的交互界面是流程中最复杂的部分之一。在宏基因组工作流程中，交互式界面可以让你更直观地浏览多维度的数据和可视化的分箱结果，并进行优化。关于 Anvi'o 的可视化界面详见：http://merenlab.org/2016/02/27/the-anvio-interactive-interface/

anvi-interactive -p SAMPLES-MERGED/PROFILE.db -c contigs.db

Reference

• http://merenlab.org/tutorials/assembly-based-metagenomics/

• http://merenlab.org/2016/06/22/anvio-tutorial-v2/

猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。