之前详细介绍了宏基因组二、三代混合组装软件OPERA-MS的Nature Biotechnology文章。详见下文：

• NBT：宏基因组二、三代混合组装软件OPERA-MS

今天带来软件的使用经验，主要参考如下官方文档：

https://github.com/CSB5/OPERA-MS

同时使用自己的数据进行经验分享。

简介

OPERA-MS是一种混合式宏基因组组装软件，结合了短读长和长读长技术的优点，可提供高质量的组装，解决了短读组装的低连续性和长读组装的低碱基对质量的问题。 OPERA-MS已在使用不同的长读技术（包括牛津纳米孔，PacBio和Illumina合成长读长）测序的合成群落和真实群落中进行了广泛的测试，并且对读长数据中的噪声特别稳健。

OPERA-MS采用分段组装策略，该策略旨在利用低覆盖率的长读取数据来改善基因组组装。它从构建一个短读长的宏基因组拼接集（默认值：MEGAHIT）开始，该程序集可以很好地表示宏基因组中的基础序列，但可能会被片段化。然后，将长读取和短读取映射到拼接集，以识别重叠群之间的连通性并计算读取的覆盖率信息。这是OPERA-MS算法核心的基础，该算法将利用基于贝叶斯模型的方法利用覆盖率和连接性信息将重叠群准确地聚类到基因组中。 OPERA-MS的另一个重要优点是，它可以对宏基因组中的菌株进行去卷积，可以选择使用参考基因组中的信息来对此进行支持。对于从容易出错的长读的拼接开始的流程来说，这从根本上是挑战。聚类后，OPERA-LG使用准确的支架(scaffold)对各个基因组进行进一步的延伸和缺口填充。

OPERA-MS可以从宏基因组数据集中组装几乎完整的基因组，而其长读长覆盖率仅为9倍。软件的设计是保守的，避免进行激进的组装，这是许多现代组装程序倾向于报告高连续性统计信息的策略。 OPERA-MS用于人类肠道微生物组数据，可提供数百个高质量的草图基因组，其中大多数N50 > 100kbp。我们观察到完整质粒的装配，其中许多是新颖的，并且包含以前看不见的抗性基因组合。此外，即使在复杂的宏基因组中存在多种物种的菌株，OPERA-MS仍可以非常准确地组装基因组，从而使我们能够使用纵向数据将质粒与宿主基因组关联起来。

使用

软件安装

# 下载软件并编译
git clone https://github.com/CSB5/OPERA-MS.git
cd OPERA-MS
make
# 检查依赖关系
perl OPERA-MS.pl check-dependency

# 缺perl模块的安装
perl tools_opera_ms/install_perl_module.pl
export PERL5LIB="/home/$USER/perl5/lib/perl5${PERL5LIB:+:${PERL5LIB}}";

内置测试数据运行

输入文件：

• 短序列：—short-read1/2
• 长序列：—long-read
• (可选)预组装结果：—contig-file
• 输出目录：—out-dir
• 关闭有参聚类：—no-ref-clustering

输入预组装的重叠群，短和长读长序列。关闭有参聚类，24线程加速，输出至RESULTS目录，耗时1m31s。结果N50有39K

cd test_files
time perl ../OPERA-MS.pl \
    --contig-file contigs.fasta \
    --short-read1 R1.fastq.gz \
    --short-read2 R2.fastq.gz \
    --long-read long_read.fastq \
    --no-ref-clustering --num-processors 24 \
    --out-dir RESULTS 2> log.err
cat RESULTS/assembly.stats

计算共分为9步。

1. MEGAHIT/Metaspades短序列组装，指定—contig-file则跳过
2. 上读长比对
3. 短读长比对
4. 层级聚类
5. 有参聚类
6. 菌株聚类
7. 其他聚类组装
8. 补洞
9. 纠错

不指定contig文件时，结果下降非常大；N50由271k下降到18k

输出结果

可以在指定的输出目录（即RESULTS）中找到以下输出文件。 文件contigs.polished.fasta（如果未纠错组装的则为contigs.fasta）包含组装的contig，assembly.stats提供总体组装统计信息（例如组装统计，N50，最长的重叠群等），而contig_info.txt提供组装后的重叠群的详细概述。

最后，可以在目录RESULTS / opera_ms_clusters / all中找到OPERA-MS菌株水平聚类（每个菌株一个fasta文件），而cluster_info.txt提供了这些群落的组装统计的详细概述。 请注意，这些簇是为产出高质量组装而构造的，因此是保守的。 重叠群可以使用诸如MaxBin2或MetaBAT2之类的方法进一步进行分箱。

• assembly.stats # 组装统计
• contig_info.txt # 重叠群信息，如长度
• contigs.fasta # 组装结果，重叠群
• contigs.polished.fasta # 纠错后的重叠群
• intermediate_files/ # 中间文件
• opera-ms-utils.config # 参数文件

常用参数

显示帮助

perl ../OPERA-MS.pl --help

• 参考基因组目录：—genome-db，使用GTDB的基因组数据库和分类信息，详见数据库部分
• 关闭参考基因组聚类：—no-ref-clustering
• 关闭菌株水平聚类:—no-strain-clustering
• 关闭纠错：—no-polishing，默认使用短读长软件Pilon进行纠错，结果为contigs.polished.fasta。样本高覆盖度和/或高度复杂时需要大量计算时间。
• 长序列比对软件选择：—long-read-mapper，默认为 blasr，可选minimap2，仅测试了v2.11-r79
• 短序列组装：—short-read-assembler，默认为 megahit ，可选spades
• 禁用补洞：—no-gap-filling
• kmer大小：—kmer-size，默认60
• 重叠群长度阈值：—contig-len-thr，默认500，小于会被去除
• 重叠末端截断长度：—contig-edge-len，计算重叠群覆盖度时，在每个重叠群末端过滤一定长度以避免低比对效率，默认为80
• 重叠 —contig-window-len，用于覆盖度估计的窗口大小，默认为340。 推荐使用最小长度 — 2 X 末端截断长度(500 - 2 X 80 = 340)
• 使用预组装的重叠群：—contig-file，使用会跳过短序列组装步骤。
• 线程数：—num-processors，默认为2，可以使用10-20等适合的计算资源加速

自定义参数分析

关闭参考聚类，和菌株层面聚类，关闭纠错(可显著减少计算时间)，使用megahit组装，使用24个线程。耗时1min，N50 仅为 18k

perl ../OPERA-MS.pl \
    --short-read1 R1.fastq.gz \
    --short-read2 R2.fastq.gz \
    --long-read long_read.fastq \
    --short-read-assembler megahit \
    --num-processors 24 \
    --no-ref-clustering --no-strain-clustering --no-polishing \
    --out-dir RESULTS1 2> log.err
cat RESULTS1/assembly.stats

关闭参考聚类，和菌株层面聚类，组装方法改为spades(需要预安装，默认没有)，使用24个线程。耗时4min，N50 40k

perl ../OPERA-MS.pl \
    --short-read1 R1.fastq.gz \
    --short-read2 R2.fastq.gz \
    --long-read long_read.fastq \
    --short-read-assembler spades \
    --num-processors 24 \
    --no-ref-clustering --no-strain-clustering --no-polishing \
    --out-dir RESULTS2 2> log.err
cat RESULTS2/assembly.stats

OPERA-MS旨在与深度的短读测序一起使用，但就长读测序可以在较低的覆盖范围内工作。 实际上，建议短读覆盖率> 15倍，而OPERA-MS可以使用低至9倍的长读覆盖率来提高装配的连续性。 基于此，我们建议至少9Gbp的短读长数据和3Gbp的长读长数据，以允许细菌基因组在宏基因组中相对丰度为1％时进行组装。

测试数据结果评估

测试数据默认提供了contigs.fasta预组装的结果(pre)，还有使用默认参数下利用三代提高的结果(default)。同时我还使用测试数据分别使用megahit和metaspades从头组装的结果。使用quast进行评估(安装方法为conda install quast)：

quast.py --label "pre,default,megahit,metapasdes" \
    contigs.fasta \
    RESULTS/contigs.polished.fasta \
    RESULTS1/contigs.polished.fasta \
    RESULTS2/contigs.fasta \
    -o quast

可以看到基于预组装并结合三代修正的结果最高，重叠群最大997kb，N50高达391kb。然而奇怪的是我从头组装的结果N50却只有20k和28k，而作者预组装的有61k。是我的组装参数不对，还是作者使用了更大的测试数据集吗？

此外，我还用我自己的二、三代数据测试了metaspades和OPERA-MS，结果如下：

整体上在总长度方面metaspades更好。而N50则是opera-ms更好，但优势不明显。opera-ms及polished的结果基本一致，在长度上也重合。

辅助工具

OPERA-MS-UTILS脚本对组装进行后处理。 实现简化的分析工具，可以计算短读序列和长读序列之间代表类群的一致性，从而对分箱、重叠群进行评估，从而可以评估分箱质量并从宏基因组中鉴定新物种的基因组。

这些工具的完整说明详见以下链接：

https://github.com/CSB5/OPERA-MS/wiki/OPERA-MS-Utilities

OPERA-MS-UTILS提供了不同的实用程序，对OPERA-MS组装，质控不同技术读长，并方便地执行必要的组装后宏基因组分析：

• opera-ms-db：生成用于OPERA-MS基于参考聚类的自定义数据库
• read-concordance：计算长和短读长测序数据之间的丰度分布相关性
• binning：使用MetaBAT2或MaxBin2分箱OPERA-MS的重叠群
• bin-evaluation：使用CheckM评估bin质量
• novel-species：OPERA-MS宏基因组组装基因组（MAG）的新物种分析，以鉴定密切相关的物种和新物种的MAG

数据库下载

方法1. 程序命令下载数据库，国内dropbox可能无法下载

perl OPERA-MS.pl install-db

方法2. 手动构建，输入为gtdb的taxonomy和基因组文件

https://gtdb.ecogenomic.org/ —— Downloads —— GTDB Data Files —— lastest (2020年6月发行的95版)

mkdir -p gtdb/95/genomic_files_reps/ && cd gtdb/95/genomic_files_reps/
# 基因组32G
wget -c https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/genomic_files_reps/gtdb_genomes_reps.tar.gz
# 解压
tar xvzf gtdb_genomes_reps_r95.tar.gz
rename 's/_genomic.fna.gz//' gtdb_genomes_reps_r95/*
# 获得基因组ID
ls gtdb_genomes_reps_r95/GC*|cut -f 2 -d '/'|sed 's/_genomic.fna.gz//' > gtdb_genomes_reps_r95.id
cd ..

# 物种注释
# 古菌3k
wget -c https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/ar122_taxonomy.tsv
# 细菌，190k
wget -c https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/bac120_taxonomy.tsv
# 合并，并修改与基因组文件名一致
cat ar122_taxonomy.tsv bac120_taxonomy.tsv | sed 's/^GB_//'> genome_taxonomy.tsv
# 筛选与解压一致的基因组
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$0} NR>FNR{print a[$1]}' genome_taxonomy.tsv genomic_files_reps/gtdb_genomes_reps_r95.id > gtdb_genomes_reps_r95_taxonomy.tsv

# 建索引
db=/home/meta/db/gtdb/95
python OPERA-MS-UTILS.py opera-ms-db \
   --genomes-dir ${db}/genomic_files_reps/gtdb_genomes_reps_r95 \
   --taxonomy ${db}/gtdb_genomes_reps_r95_taxonomy.tsv \
   --db-name db

复制文件结束后报错

Traceback (most recent call last):
  File "OPERA-MS-UTILS.py", line 287, in <module>
    main(args)
  File "OPERA-MS-UTILS.py", line 200, in main
    opera_ms_db(args.genomes_dir, args.taxonomy, args.db_name, args.thread)
  File "OPERA-MS-UTILS.py", line 121, in opera_ms_db
    read_taxonomy_file(genomes_dir, taxonomy, genome_db, genome_list, genome_size)
  File "OPERA-MS-UTILS.py", line 71, in read_taxonomy_file
    tax_info = line_list[1].split(";")
IndexError: list index out of range