macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

コアゲノム推定にメタゲノムアセンブルゲノムを活用するためのロバストなベイズアプローチ mOTUpan

 

 近年のシーケンサーバイオインフォマティクスの進歩により、メタゲノムアセンブルゲノム(MAG)やシングルセルアセンブルゲノム(SAG)を通じて、環境に関連する未培養クレードのゲノムを提供し、生命の系譜を拡大している。このような多様性の拡大により、微生物の集団構造に関する新しい知見が得られる一方で、コアゲノム推定に利用できるほとんどのツールは、ゲノムの完全性に過敏に反応する。その結果、環境ゲノム解析によって明らかになった膨大な系統的多様性の大部分が、このような解析から除外されたままになっている。本論文では、多様なゲノムアセンブリに対して、コアゲノムを計算するための新しい反復ベイズ法を紹介する。本手法は、多様なゲノムアセンブリに対して、各遺伝子クラスターがコアゲノムに属する可能性、あるいはアクセサリゲノムに属する可能性を、ゲノムアセンブリにおけるその存在/非存在パターンの確率を計算することにより推定する。コアゲノム予測は計算効率が高く、数千ゲノムまでスケールアップ可能である。高品質なゲノムに対しては、最新ツールであるRoaryやPPanGGoLiNと同等の推定値を示し、より低い完全性の閾値ではそれらを凌駕する。mOTUpanは、特定のコアゲノム予測の品質を推定するためのブートストラップ手順をラップしている。

各ゲノムは、まず、ゲノムにコードされた形質の集合として記述される。ここでは遺伝子クラスターを用いるが、mOTUpanはその形質にとらわれず、遺伝子、COG、機能アノテーション、その他ゲノムがコードするあらゆる形質を用いることが可能である。2つの仮説(コアーアクセサリー形質)それぞれについて、各ゲノムについて推定されたゲノム完全性事前分布を使用して確率が計算される(ゲノムの完全性はCheckMを使用して計算することもできるし、固定値を使用することもできる)。その結果、最も可能性の高い形質カテゴリー(コアまたはアクセサリ)が、その形質のクラスとして選ばれる。この新しい分類を用いて、事後的な完全性推定値を計算し、これを2回目の反復のための事前分布として用いることができる。そして、収束するまでこのプロセス全体を繰り返す。

Roaryのような利用可能なツールのほとんどは、コアゲノムを定義するためのハードな存在/不在の閾値に依存している。この制限により、不完全で断片化したMAGやSAGを扱う場合、そのようなツールはほとんど使用できない。PPanGGOLiNはネットワークの性質上、ある程度の不完全性に対応することができるが、シンテニーのパターンを探してゲノムの残存率を決定しているため、(MAGやSAGによく見られる)あまりに断片化するとゲノムの残存率の計算で問題が発生する可能性がある。MAGやSAGに代表される種の場合、PPanGGOLiNによって「持続性」と分類された遺伝子はコアの一部である可能性が非常に高いが、このアプローチではいくつかのコア遺伝子を見落としてしまう可能性がある。一方、mOTUpanは、不完全性と断片化の両方の制限を回避し、不完全で断片化したMAGとSAGアセンブリに対して、コアゲノムとパンゲノムを強固に推定することができる。mOTUpanはまた、コアゲノム/パンゲノム分割のためのブートストラップ偽発見率および感度を計算する。ビニングツールが誤って残りのMAGと一緒にクラスタリングすることがあるので、そのゲノムの真の一部でないかもしれないコンティグによってMAGが汚染されるという広くかつ妥当な懸念が存在する。しかし、mOTUpanは、Coreと注釈された遺伝子は汚染される可能性が非常に低く、ゲノムの品質予測に利用できるため、この既知の問題に別の方法で対処することができる。このように、mOTUpanは、CheckMが推定した完全性の値をユーザーが改良することを可能にする。また、他のゲノムセット、ウイルス、プラスミドについても、mOTUpanは完全性の推定が可能である。

 

 

mOTUlizer: Bayesian approach to leverage metagenomic bins for pan-genome analysis

 

Githubより

mOTUsという、リードから直接OTUテーブルを作成するツールがありますが、それとは異なるツールです(mOTUs2)。この mOTUpanは、多少なりとも出所の怪しい、密接に関連したMAG/ゲノム/bin/SUBのグループを解析するためのユーティリティです。現在、いくつかのプログラムで構成されています。

 

インストール

ubuntu18にてcondaで環境を作って導入した。

依存

  • fastANI
  • mmseqs2

Github

#conda (link)
mamba create -n mOTUpan -y
conda activate mOTUpan
mamba install -c bioconda motulizer -y

#pip (pypi)
pip install mOTUlizer

>  mOTUlize.py -h

usage: mOTUlize [-h] [--output [OUTPUT]] [--force] [--checkm [CHECKM]]

                [--similarities [SIMILARITIES]] [--fnas [FNAS [FNAS ...]]]

                [--prefix [PREFIX]] [--MAG-completeness [MAG_COMPLETENESS]]

                [--MAG-contamination [MAG_CONTAMINATION]]

                [--SUB-completeness [SUB_COMPLETENESS]]

                [--SUB-contamination [SUB_CONTAMINATION]]

                [--similarity-cutoff [SIMILARITY_CUTOFF]] [--threads THREADS]

                [--keep-simi-file [KEEP_SIMI_FILE]] [--txt] [--long]

                [--version]

 

From a set of genomes, makes metagenomic Operational Taxonomic Units (mOTUs).

By default it makes a graph of 95% (reciprocal) ANI (with fastANI) connected

MAGs (with completeness > 40%, contamination < 5%). The mOTUs will be the

connected components of this graph, to which smaller "SUBs" with ANI > 95% are

recruited. If similarities provided, it should be a TAB-separated file with

columns as query, subject and similarity (in percent, e.g. [0-100]) if you

also provide fasta-files (for stats purpouses) query and names should

correspond to the fasta-files you provide. If the columns are file names, the

folders are removed (mainly so it can read fastANI output directly).

 

optional arguments:

  -h, --help            show this help message and exit

  --output [OUTPUT], -o [OUTPUT]

                        send output to this file

  --force, -f           force execution (overwritting existing files)

  --checkm [CHECKM], -k [CHECKM]

                        checkm file (or whatever you want to use as

                        completness), if not provided, all genomes are assumed

                        to be seed MAG (e.g complete enough)

  --similarities [SIMILARITIES], -I [SIMILARITIES]

                        file containing similarities between MAGs, if not

                        provided, will use fastANI to compute one

  --fnas [FNAS [FNAS ...]], -F [FNAS [FNAS ...]]

                        list of nucleotide fasta-files of MAGs or whatnot

  --prefix [PREFIX], -n [PREFIX]

                        prefix for the mOTU names, default : mOTU_

  --MAG-completeness [MAG_COMPLETENESS], --MC [MAG_COMPLETENESS], -M [MAG_COMPLETENESS]

                        completeness cutoff for seed MAGs, default : 40

  --MAG-contamination [MAG_CONTAMINATION], --Mc [MAG_CONTAMINATION], -m [MAG_CONTAMINATION]

                        contamination cutoff for seed MAGs, default : 5

  --SUB-completeness [SUB_COMPLETENESS], --SC [SUB_COMPLETENESS], -S [SUB_COMPLETENESS]

                        completeness cutoff for recruited SUBs, default : 0

  --SUB-contamination [SUB_CONTAMINATION], --Sc [SUB_CONTAMINATION], -s [SUB_CONTAMINATION]

                        contamination cutoff for recruited SUBs, default : 100

  --similarity-cutoff [SIMILARITY_CUTOFF], -i [SIMILARITY_CUTOFF]

                        distance cutoff for making the graph, default : 95

  --threads THREADS, -t THREADS

                        number of threads (default all, e.g. 32)

  --keep-simi-file [KEEP_SIMI_FILE], -K [KEEP_SIMI_FILE]

                        keep generated similarity file if '--similarities' is

                        not provided, does nothing if '--similarity' is

                        provided

  --txt, -T             the '--fnas' switch indicates a file with paths

  --long                longform output, a JSON-file with a lot more

                        information (might be cryptic...)

  --version, -v         get version

 

Let's do this

> mOTUpan.py -h

usage: mOTUpan.py [-h] [--output OUTPUT] [--force] [--checkm CHECKM]

                  [--seed [SEED]] [--length_seed] [--random_seed]

                  [--genome2gene_clusters_only] [--precluster]

                  [--faas [FAAS [FAAS ...]]] [--txt]

                  [--gene_clusters_file [GENE_CLUSTERS_FILE]] [--name NAME]

                  [--long] [--boots BOOTS] [--max_iter MAX_ITER]

                  [--threads THREADS] [--version]

 

From a buch of amino-acid sequences or COG-sets, computes concensus AA/COG

sets. Returns all to stdout by default.

 

optional arguments:

  -h, --help            show this help message and exit

  --output OUTPUT, -o OUTPUT

                        send output to this file

  --force, -f           force execution answering default answers

  --checkm CHECKM, -k CHECKM

                        checkm file if you want to seed completnesses with it,

                        accepts concatenations of multiple checkm-files, check

                        manual for more formating options

  --seed [SEED], -s [SEED]

                        seed completeness, advice a number around 90 (95

                        default), this is the default completeness prior

  --length_seed, --ls   seed completeness by length fraction [0-100]

  --random_seed, --rs   random seed completeness between 5 and 95 percent

  --genome2gene_clusters_only

                        returns genome2gene_clusters only

  --precluster          precluster proteomes with cd-hit, mainly for legacy

                        reasons, mmseqs2 is faaaaaast

  --faas [FAAS [FAAS ...]], -F [FAAS [FAAS ...]]

                        list of amino-acids faas of MAGs or whatnot, or a text

                        file with paths to the faas (with the --txt switch)

  --txt                 the '--faas' switch indicates a file with paths

  --gene_clusters_file [GENE_CLUSTERS_FILE], --gene_clusterss [GENE_CLUSTERS_FILE], -c [GENE_CLUSTERS_FILE]

                        file with COG-sets (see doc)

  --name NAME, -n NAME  if you want to name this bag of bins

  --long                longform output, a JSON-file with a lot more

                        information (might be cryptic...)

  --boots BOOTS, -b BOOTS

                        number of bootstraps for fpr and recall estimate

                        (default 0), careful, slows down program linearly

  --max_iter MAX_ITER, -m MAX_ITER

                        max number of iterations (set to one if you have only

                        few traits, e.g. re-estimation of completeness is

                        nonsensical)

  --threads THREADS, -t THREADS

                        number of threads (default all, e.g. 32), only gene

                        clustering is multithreaded right now, the rest is to

                        come

  --version, -v         get version

 

Let's do this

 

 

テストラン

OTUを作成し、統計情報を得る。ゲノムを指定する。

git clone https://github.com/moritzbuck/mOTUlizer.git
cd mOTUlizer/
mOTUlize.py --fnas example_files/fnas/*.fna -o output.tsv

mOTUlize.py は、ゲノム(ここでは、ゲノムという言葉を、同じ生物/集団に由来すると思われる塩基配列の集合の略語として使うが、不完全でも冗長でも汚染されていても構わない)を受け取り、それらをmOTUクラスタリングする。類似度スコア(デフォルトでは fastANI が計算した ANI、ユーザが他の類似度を指定可能)を用いて、類似度を特定の値(デフォルトでは 95%)で閾値処理することにより、(ユーザ定義の)より質の高いゲノム(歴史的理由により MAG と呼ぶ)を基にネットワークを構築する。このグラフの連結成分は mOTU である。さらに、低品質ゲノム(SUB, )は、類似度が閾値以上であれば、最も類似している MAG の mOTU に採用される(Githubより)。

output.tsv

 

 

protein.fastaを指定する。

mOTUpan.py --faas example_files/faas/*.faa -o output.tsv

mOTUpan.py は、ある形質があるゲノムアセンブリのコアゲノムに含まれる可能性など、あるゲノムアセンブリの全てに含まれると予想される遺伝子コード化された形質の尤度を計算する。基本的には、mOTUpan に、対象となるゲノムのプロテオームセット(例えば同じ mOTU や Genus から)と、これらのゲノムの完全性の事前評価値(例えば checkm 出力や固定値)を与えると、mmseqs2 を使って遺伝子クラスタが計算されるようになっている。これらの遺伝子クラスタそれぞれについて、コアである可能性とコアでない可能性を計算し、これらの尤度の比によって、形質がコアとみなされるかどうかが決定される。この新しいパーティショニングは、完全性事前分布を更新するために用いられ、収束するまで繰り返し計算される(Githubより)。

output.tsv

そのタンパク質の出現回数、コアであるかのどうかの尤度とコアかアクセサリかの判定結果、そのタンパク質を持っているゲノムファイル名と遺伝子名、1つのゲノムでの出現個数がまとめられている。

 

  • mOTUpanは、通常は全てのゲノムが95%の以上のANIで定義されたコンパクトなクラスタ内にあるゲノムを使用する。
  • ドラフトゲノム、コンプリートゲノム、MAG、SAGのいずれもあり得る(roaryのように対応していないツールが多い)。

 

引用

mOTUpan: a robust Bayesian approach to leverage metagenome assembled genomes for core-genome estimation
Moritz Buck,  Maliheh Mehrshad,  Stefan Bertilsson

bioRxiv, Posted June 25, 2021

 

関連