アセンブリ時のハプロタイプ重複に対処するためのツール群 Mabs

　真核生物のゲノム解析は、ゲノム解読法の進歩にもかかわらず、未だエラーフリーには至っていない。真核生物のゲノムアセンブリの問題の中には、対立遺伝子が誤ってパラロガスとしてアセンブリされるケースとして現れる、いわゆる「ハプロタイプ重複」と呼ばれるものがある。ハプロタイプ重複は、遺伝子ファミリーの拡大という幻想をもたらし、ゲノムの進化や機能に関して誤った結論を導く可能性があるため、危険である。本発表では、著名なゲノムアセンブラであるHifiasmとFlyeのパラメータ最適化ツール群Mabsを紹介する。Mabsは、HifiasmとFlyeのパラメータを最適化することで、できるだけ正確に遺伝子が組み合わされたゲノムアセンブリを作成しようとするものである。6種類の真核生物ゲノムを用いたテストでは、6件中6件で、MabsはHifiasmとFlyeをデフォルトのパラメータで実行した場合よりも、より正確に遺伝子が組み立てられたアセンブリ配列を作成したことが確認された。また、Mabs、Hifiasm、Flyeのアセンブリをハプロタイプ重複除去ツールPurge_dupsで後処理したところ、6件中5件でMabsの方がより正確に遺伝子がアセンブリされた。MabsはPythonで書かれており、https://github.com/shelkmike/Mabs で入手できる。

Mabsは基本的にゲノムアセンブラのパラメータオプティマイザとして動作する。最適化器として動作し、AG（"Accurately assembled Genes ")を最大化するパラメータの値を見つけようとする。

インストール

Github

git clone https://github.com/shelkmike/Mabs.git
cd Mabs/
bash install.sh

> ./mabs-flye.py

Mabs-flye, a program for genome assembly.

Main options:

1) --nanopore_reads Path to Nanopore reads.

2) --pacbio_clr_reads Path to PacBio CLR reads, also known as

"old PacBio" reads.

3) --pacbio_hifi_reads Path to PacBio HiFi reads, also known as

CCS reads.

[any of the above files may be in FASTQ or FASTA, gzipped or not]

4) --download_busco_dataset Name of a file from

http://mikeshelk.site/Data/BUSCO_datasets/Latest/ . It should be the

most taxonomically narrow dataset for your species. For example, for a

human genome assembly, use "--download_busco_dataset

primates_odb10.2021-02-19.tar.gz" and for a drosophila genome assembly

use "--download_busco_dataset diptera_odb10.2020-08-05.tar.gz".

Mabs-flye will download the respective file. This option is mutually

exclusive with "--local_busco_dataset".

5) --threads Number of CPU threads to be used by Mabs-flye. The

default value is 10.

6) --output_folder Output folder for Mabs-flye results. The

default is "Mabs_results".

7) --number_of_busco_orthogroups How many BUSCO orthogroups

should Mabs-flye use. Should be either a positive integral value or

"all" to use all orthogroups. The default value is 1000.

8) --genome_size Haploid genome size. Should be either "auto"

for automatic estimation, or a number ending with "k", "m" or "g". For

example, 1.5g means 1.5 gigabases. The default value is "auto".

9) --max_intron_length Maximum allowed length of an intron.

Should be either "from_BUSCO" to use a value from a BUSCO dataset, or

a number, possibly ending with "k", "m" or "g". For example, 20k means

20 kilobases. The default is "from_BUSCO". Change --max_intron_length

if you assemble a genome with unusually long introns.

10) --local_busco_dataset Path to a local BUSCO dataset,

manually pre-downloaded from

http://mikeshelk.site/Data/BUSCO_datasets/Latest/ or

http://busco-data.ezlab.org/v5/data/lineages/. Example:

"--local_busco_dataset

/home/test/Data/primates_odb10.2021-02-19.tar.gz". May be a .tar.gz

file or a decompressed folder. This option is mutually exclusive with

"--download_busco_dataset".

Informational options:

11) --help Print this help.

12) --version Print the version of Mabs-flye.

13) --run_test Assemble a small dataset to test whether

Mabs-flye was installed properly. The assembly takes approximately 10

minutes. If a non-empty file

./Mabs_results/The_best_assembly/assembly.fasta appears, then

Mabs-flye was installed successfully.

Example 1:

mabs-flye.py --nanopore_reads nanopore_reads.fastq

--download_busco_dataset eudicots_odb10.2020-09-10.tar.gz --threads 40

Example 2:

mabs-flye.py --nanopore_reads nanopore_reads.fastq --pacbio_hifi_reads

pacbio_hifi_reads.fastq --download_busco_dataset

diptera_odb10.2020-08-05.tar.gz --threads 40

> ./mabs-hifiasm.py

Mabs-hifiasm, a program for genome assembly.

Main options:

1) --pacbio_hifi_reads Path to PacBio HiFi reads, also known as CCS reads.

2) --short_hi-c_reads_R1 Path to short Hi-C reads, first reads from a pair.

3) --short_hi-c_reads_R2 Path to short Hi-C reads, second reads from a pair.

[any of the above files may be in FASTQ or FASTA, gzipped or not]

4) --download_busco_dataset Name of a file from

http://mikeshelk.site/Data/BUSCO_datasets/Latest/ . It should be the

most taxonomically narrow dataset for your species. For example, for a

human genome assembly, use "--download_busco_dataset

primates_odb10.2021-02-19.tar.gz" and for a drosophila genome assembly

use "--download_busco_dataset diptera_odb10.2020-08-05.tar.gz".

Mabs-hifiasm will download the respective file. This option is

mutually exclusive with "--local_busco_dataset".

5) --threads Number of CPU threads to be used by Mabs-hifiasm.

The default value is 10.

6) --output_folder Output folder for Mabs-hifiasm results. The

default is "Mabs_results".

7) --number_of_busco_orthogroups How many BUSCO orthogroups

should Mabs-hifiasm use. Should be either a positive integral value or

"all" to use all orthogroups. The default value is 1000.

8) --ploidy Ploidy of the genome. The default value is 2.

9) --genome_size Haploid genome size. Should be either "auto"

for automatic estimation, or a number, possibly ending with "k", "m"

or "g". For example, 1.5g means 1.5 gigabases. The default value is

"auto".

10) --max_intron_length Maximum allowed length of an intron.

Should be either "from_BUSCO" to use a value from a BUSCO dataset, or

a number, possibly ending with "k", "m" or "g". For example, 20k means

20 kilobases. The default is "from_BUSCO". Change --max_intron_length

if you assemble a genome with unusually long introns.

11) --local_busco_dataset Path to a local BUSCO dataset,

manually pre-downloaded from

http://mikeshelk.site/Data/BUSCO_datasets/Latest/ or

http://busco-data.ezlab.org/v5/data/lineages/. Example:

"--local_busco_dataset

/home/test/Data/primates_odb10.2021-02-19.tar.gz". May be a .tar.gz

file or a decompressed folder. This option is mutually exclusive with

"--download_busco_dataset".

Informational options:

12) --help Print this help.

13) --version Print the version of Mabs-hifiasm.

14) --run_test Assemble a small dataset to test whether

Mabs-hifiasm was installed properly. The assembly takes approximately

10 minutes. If a non-empty file

./Mabs_results/The_best_assembly/assembly.fasta appears, then

Mabs-hifiasm was installed successfully.

Example 1:

mabs-hifiasm.py --pacbio_hifi_reads hifi_reads.fastq

--download_busco_dataset eudicots_odb10.2020-09-10.tar.gz --threads 40

Example 2:

mabs-hifiasm.py --pacbio_hifi_reads hifi_reads.fastq

--short_hi-c_reads_R1 hi-c_reads_trimmed_R1.fastq

--short_hi-c_reads_R2 hi-c_reads_trimmed_R2.fastq

--download_busco_dataset diptera_odb10.2020-08-05.tar.gz --threads 40

実行方法

Mabs-flyeはNanoporeリードとPacBio CLRリードを対象としている。Mabs-hifiasmは、PacBio HiFiを使用したアセンブリを対象としている。

hifiasm.pyのコマンド例

mabs-hifiasm.py --pacbio_hifi_reads hifi_reads.fastq --download_busco_dataset eudicots_odb10.2020-09-10.tar.gz --threads 40

出力例

mabs_results/

最良のアセンブリはThe_best_assembly/に保存される。繰り返すために時間はかかる（レポジトリによると通常のアセンブリの3倍の時間）。mabs-flye.pyも同様に使える。

calculate_AGツールはMabs-hifiasmとMabs-flyeの内部で使用されているが、ユーザーのアセンブリの品質の評価にも使用できる。

calculate_AG.py --assembly contigs.fasta --nanopore_reads nanopore_reads.fastq --local_busco_dataset /mnt/lustre/username/Datasets/eudicots_odb10 --threads 40

終了するとAG値のテキストとシングルコピーオーソログループ、真のマルチコピーオーソログループ、偽のマルチコピーオーソログループの遺伝子数が保存される。また、sinaplotと呼ばれる図が保存される。これはマルチコピー・オーソログループ遺伝子のカバレッジ分布が単峰性か二峰性かを確認するために利用できる。明確に二峰性なら品質はかなり悪いと推測される。

レポジトリより

異なる技術で作成された複数のリードデータセットがある場合、これらのオプションは同時に使用できる。
MabsはBUSCO遺伝子のアセンブリの品質を、"AG "（Accurately assembled Genes）と呼ばれる特別な指標を用いて評価する。AGの計算方法については、レポジトリのcalculate_AGで説明されている（当然だが論文中にはより詳しい説明がある）。
基本的にAGが大きいほど、アセンブリが良好であることを意味する。
"Mabs" は "Busco Score を最大化する Miniasm-based Assembler" という意味。

レポジトリのQ＆Aも役立ちます。ハプロタイプの重複があまり期待できない場合にMabsを使う意味があるのか、高精度（99％以上）のNanoporeリードではどうするか、Mabsの初期バージョンの話、AGの導入動機などなど説明されています。

引用

Mabs, a suite of tools for gene-informed genome assembly

Mikhail I. Schelkunov

bioRxiv, Posted February 13, 2023