macでインフォマティクス

macでインフォマティクス

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

メタゲノムシークエンシングデータから微生物真核生物ゲノムを取り出すパイプライン Eukfinder

 

 微生物群集のホールゲノムショットガン(WGS)メタゲノムシークエンシングにより、多様な生態系に生息する微生物の原核生物や真核生物の機能、生理、進化の歴史を発見することができる。その重要性にもかかわらず、微生物真核生物のメタゲノム研究は、WGSデータから高品質の真核生物ゲノムを同定し、アセンブルすることが困難なため、原核生物の研究に比べて遅れている。この問題に対処するため、本著者らは、WGSメタゲノミクスデータから真核微生物の核ゲノムおよびミトコンドリアゲノムを復元し、アセンブルするバイオインフォマティクスパイプラインであるEukfinderを開発した。ワークフローの一環として、2つの特殊なデータベースを利用し、対象データセットや環境に合わせてカスタマイズ可能な分類法に基づいてリードを分類する。Eukfinderをヒト腸内細菌叢のWGSメタゲノミックシーケンスデータに適用し、ヒトや動物の消化管に非常に多く寄生する原生寄生生物Blastocystis sp.からゲノムを回収した。Eukfinderは、細菌リードと組み合わせたBlastocystisリードの数を変化させた一連のシミュレート腸内細菌叢データセットと、Blastocystisを含む実際のメタゲノム腸内サンプルの両方を用いてテストした。Eukfinderの結果を他の公開ワークフローと比較した。十分なリードが読まれている場合、Eukfinderは、リファレンスゲノムを用いることなく、メタゲノミックデータから多様なBlastocystisサブタイプから高品質でほぼ完全な核およびミトコンドリアゲノムを効率的にアセンブルした。さらに、十分な深さの配列サンプリングにより、Eukfinderはメタゲノムデータから真核生物のゲノムを復元するために使用される同様のツールよりも優れた性能を示した。Eukfinderは、環境メタゲノムシーケンスサンプルから、リファレンスに依存せず、培養なしで真核微生物ゲノムを研究するための有用なツールとなるだろう。

 

wiki

https://github.com/dzhao2019/eukfindertest/wiki

(マニュアルより)(a)ショートリードを用いるEukfinder_shortのランでは、まずショートリードを5つの異なる分類学的カテゴリー(Archaeal、Bacterial、Viral、Eukaryotic、Unknown)に分類する、 (b)Eukfinder_longのランでは、MAGのコンティグや、NanoporeやPacbioプラットフォームのロングリードシーケンスデータを使用する。Eukfinder_longは、EukaryoticおよびUnknownコンティグを選択するために1ラウンドの分類を行い、教師付きビニングの後、真核生物ゲノムおよびミトコンドリアゲノムを生成する。

 

インストール

依存

  • Python >= 3.7
  • ete3,numpy, pandas, joblib, pyqt, spades, seqkit, trimmomatic, bowtie2, centrifuge, acc2tax, plast

git clone https://github.com/RogerLab/Eukfinder.git
cd Eukfinder/
mamba create -n Eukfinder python=3.9 -y
conda activate Eukfinder
mamba install -c conda-forge -c bioconda -y ete3 numpy pandas joblib pyqt=5 spades seqkit trimmomatic bowtie2

#centrifuge(github)
git clone https://github.com/DaehwanKimLab/centrifuge
cd centrifuge
make -j20
sudo make install

PLAST (paper) とacc2tax(紹介)も必要

#PLAST
wget https://github.com/PLAST-software/plast-library/releases/download/v2.3.2/plastbinary_linux_v2.3.2.tar.gz
tar -zxf plastbinary_linux_v2.3.2.tar.gz
export PATH=$PWD/plastbinary_linux_v2.3.2/build/bin:$PATH

#acc2tax
git clone https://github.com/richardmleggett/acc2tax.git
cd acc2tax/
cc -o acc2tax acc2tax.c
export PATH=$PWD:$PATH
cd ../

 

データベース

4つダウンロードする。

https://perun.biochem.dal.ca/Metagenomics-Scavenger/

 

実行方法

ショートリードとロングリードを使用するコマンドが用意されている。

 

ショートリードを使うコマンドの1例

python3 eukfinder.py read_prep --r1 Illumina_1.fastq --r2 Illumina_2.fastq -n 20 --hcrop 10 -l 15 -t 15 --wsize 40 --qscore 25 --mlen 40 --hg GCF_000001405.39_GRCh38.p13_genomic.fna -o Eukfinder_output --cdb Centrifuge_NewDB_Sept2020/ -i TrueSeq2_NexteraSE-PE.fa --mhlen 40

 

  • EukfinderをBiocondaで利用できるようにするため更新が続いている。
  • 現在開発が続いており、まだヘルプなどは用意されていない。

引用

Eukfinder: a pipeline to retrieve microbial eukaryote genomes from metagenomic sequencing data.

Dandan Zhao, Dayana E. Salas-Leiva, Shelby K. Williams, Katherine A. Dunn,  Andrew J. Roger

bioRxiv, Posted December 28, 2023

 

 

(メタ)ゲノムのARGプロファイリングを行うSnakemakeパイプライン ARGprofiler

 

 メタゲノム解析は、抗菌薬耐性遺伝子(ARG)の機能や分布を理解する上で非常に有用である。しかし、研究の比較可能性を確保するために、標準化された再現可能なワークフローが必要である。現在の選択肢には、それぞれ特定の目的を念頭に設計された様々なツールや参照データベースが含まれているからである。本研究では、ARGの構成、分布、機能を研究するために、大量の生シーケンスリードを処理するワークフローARGprofilerを作成した。ARGprofilerは、14,078のユニークなARGからなるPanResデータベースを提供することで、どの参照データベースを使うべきかという課題に取り組んでいる。ARGprofilerのパイプラインは、遺伝子や微生物のアバンダンステーブルを作成するだけでなく、ARGextenderを用いてARGのフランキング領域を再構築するように設計されている。ARGextenderはKMAとSPAdesを組み合わせたバイオインフォマティックアプローチで、ターゲットとなるde novoアセンブリのためにリードをリクルートする。本目的はARGにあるが、パイプラインは、シーケンシングランを高速に検索・比較するためのMashスケッチも作成する。

 ARGprofilerパイプラインは、メタゲノムシーケンスデータの再利用をサポートするSnakemakeワークフローであり、簡単にインストールでき、https://github.com/genomicepidemiology/ARGprofilerから利用できる。

 

インストール

依存

  • Snakemake

Github

git clone https://github.com/genomicepidemiology/ARGprofiler.git
cd ARGprofiler/
mamba env create --name argprofiler --file rules/environment_argprofiler.yaml
coonda activate argprofiler

 

テストラン

レポジトリのルートにあるinput.jsonが認識される。

cd ARGprofiler/
snakemake --profile profile_argprofiler

input.json

run_accessionはリードシーケンスデータセットのENAのID。READ_TYPEはPAIREDまたはSINGLEのいずれかなので、2つ分のデータセット(片方はペアエンド、片方はシングルエンド)が解析対象になる。

ENAのデータがダウンロードされて解析される。

resultsに全ての結果が保存される。

results/

raw_reads/にはダウンロードしたシーケンスデータセットが格納され、trimmed_reads/にはトリミングされたシーケンスデータが格納される。kma_mOTUs/にはmOTUsデータベースを含む全てのアライメント結果ファイルが含まれ、kma_panres/にはPanResデータベースとのアライメント結果ファイルがすべて格納されている。Mash/には各シーケンスデータセットのマッシュスケッチが含まれている。画像では失敗しているが、argextenderディレクトリには、ARG周辺のゲノムフランキング領域を抽出したファイルが格納される。

 

どのサブディレクトリもシングルリードとペアリードの結果に分けて保存されている。

raw_reads/

 

論文とレポジトリより

  • 公開されていないfastqを使用したい場合、local_readsという名前のディレクトリにfastqファイルを配置する。
  • clusterのフラグを付けることでHPC環境でも実行できる。
  • 論文で書かれているが、いくつかの方法が比較検討され、最終的に、fastq-dlにとるENAからのリードのダウンロード、fastpによる前処理、KMA(紹介)による選択したmOTUs3(紹介)やPanResなどの参照配列データベースへのマッピング、この論文のために開発された目的の遺伝子周辺のゲノムフランキング領域を拡張するアセンブルツール:ARGextenderによる抗菌薬耐性遺伝子のターゲットフランキング領域のアセンブリ、Mash(紹介)によるシーケンシングリードのスケッチというフローが行われる。
  • 抗生物質、重金属、殺生物剤に対する耐性をコードする細菌遺伝子はこれまでに同定され、いくつかのデータベースにまとめられている。本著者らは、興味のあるこれらの遺伝子をPanResと名付けた単一のユニークなコレクションに集めた: ResFinder、ResFinderFG、CARD、MegaRes、AMRFinderPlus、ARGANNOT。これらのコレクションに加えて、環境および臨床サンプルからクローン化され、機能決定されたARGのセット(市販されていない抗生物質に対する耐性が確認されたもの)を使用している(CsabaPalタグ)。重金属はしばしば抗生物質耐性を共選択するため、BacMet v1.1からも手動でキュレーションされてコレクションされている。すべての検索配列はクラスタリングされてユニークな配列となっている。
  • メタゲノムde novoアセンブリは計算要求が高すぎて日常的な使用には適さないことが多い。より低い計算負荷でアセンブリを行うためにARGextenderが開発された。ARGextenderではKMA とSPAdesを再帰的に使用する。

 

コメント

r-optparseのインストールに失敗し、一部のプロセスが完了しなかった。改善出来たら追記します。

引用

ARGprofiler—a pipeline for large-scale analysis of antimicrobial resistance genes and their flanking regions in metagenomic datasets 
Hannah-Marie Martiny, Nikiforos Pyrounakis, Thomas N Petersen, Oksana Lukjančenko, Frank M Aarestrup, Philip T L C Clausen, Patrick Munk
Bioinformatics, Volume 40, Issue 3, March 2024

 

関連

 

 

 

 

 

複数のラージゲノム間のシンテニーを高速に検出する ntsynt

 

 近年、リファレンスグレードのゲノムアセンブリは大幅に多様化している。このような豊富なデータにより、ゲノム間の配列保存に関する情報を提供し、種の進化に関する重要な知見に貢献するゲノムのシンテニーの検出を含む、スケーラブルな複数種の比較ゲノム解析のための堅牢なツールの開発が急務となっている。ntSyntは、最小化グラフベースのアプローチにより、大規模な複数ゲノムのシンテニーブロックを計算するスケーラブルなユーティリティである。ntSyntは、約3Gbpの複数のゲノムを用いて、34GBのメモリを用い、最大2時間で79-100%のカバー率を持つシンテニーブロックを生成する。既存の最先端手法と比較して、ntSyntは多様な入力ゲノム配列とシンテニーブロックの粒度に柔軟に対応できる。我々(著者ら)は、ntSyntが促進するマクロシンテニックゲノム解析が、生命のツリー全体にわたる種内および種間の重要な進化的洞察を生み出す上で、広範な有用性を持つことを期待している。

 

wiki

https://github.com/bcgsc/ntSynt/wiki/de-novo-statistics-summary

visualization

https://github.com/bcgsc/ntSynt/tree/main/visualization_scripts

 

インストール

依存

  • python 3.9+ with modules:
  • intervaltree
  • pybedtools
  • ncls
  • python-igraph

Github

mamba create -n ntsynt python=3.9 -y
conda activate ntsynt
mamba install -c bioconda -c conda-forge ntsynt -y

git clone https://github.com/bcgsc/ntSynt.git
cd ntSynt/
./run_ntSynt_demo.sh 

> ntSynt -h

usage: ntSynt [-h] -d DIVERGENCE [-p PREFIX] [-k K] [-w W] [-t T] [--fpr FPR] [-b BLOCK_SIZE] [--merge MERGE] [--w_rounds W_ROUNDS [W_ROUNDS ...]] [--indel INDEL] [-n] [--benchmark] [-f] [--dev] [-v] fastas [fastas ...]

 

ntSynt: Multi-genome synteny detection using minimizer graphs

 

positional arguments:

  fastas                Input genome fasta files

 

optional arguments:

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

  -d DIVERGENCE, --divergence DIVERGENCE

                        Approx. maximum percent sequence divergence between input genomes (Ex. -d 1 for 1% divergence).

                        This will be used to set --indel, --merge, --w_rounds, --block_size

                        See below for set values - You can also set any of those parameters yourself, which will override these settings.

  -p PREFIX, --prefix PREFIX

                        Prefix for ntSynt output files [ntSynt.k<k>.w<w>]

  -k K                  Minimizer k-mer size [24]

  -w W                  Minimizer window size [1000]

  -t T                  Number of threads [12]

  --fpr FPR             False positive rate for Bloom filter creation [0.025]

  -b BLOCK_SIZE, --block_size BLOCK_SIZE

                        Minimum synteny block size (bp)

  --merge MERGE         Maximum distance between collinear synteny blocks for merging (bp). 

                        Can also specify a multiple of the window size (ex. 3w)

  --w_rounds W_ROUNDS [W_ROUNDS ...]

                        List of decreasing window sizes for synteny block refinement

  --indel INDEL         Threshold for indel detection (bp)

  -n, --dry-run         Print out the commands that will be executed

  --benchmark           Store benchmarks for each step of the ntSynt pipeline

  -f, --force           Run all ntSynt steps, regardless of existing output files

  --dev                 Run in developer mode to retain intermediate files, log verbose output

  -v, --version         show program's version number and exit

 

Default parameter settings for divergence values:

< 1% divergence:    --block_size 500 --indel 10000 --merge 10000 --w_rounds 100 10

1% - 10% divergence:    --block_size 1000 --indel 50000 --merge 100000 --w_rounds 250 100

> 10% divergence:    --block_size 10000 --indel 100000 --merge 1000000 --w_rounds 500 250

If any of these parameters are set manually, those values will override the above.

 

If you have any questions about ntSynt, please create a GitHub issue: https://github.com/bcgsc/ntSynt

 

 

テスト

git clone https://github.com/bcgsc/ntSynt.git
cd ntSynt/
./run_ntSynt_demo.sh 

テストには線虫の染色体の一部が使われている。

 

実行方法

ゲノムを指定する。生のfastaファイルを受け付ける(gzip圧縮されているとエラーを起こす)。

ntSynt -d 5 assembly1.fa assembly2.fa assembly3.fa
  • -d    Approx. maximum percent sequence divergence between input genomes 

 

出力例

 

メインの出力ファイルは<prefix>.synteny_blocks.tsvで、このTSVファイルに計算されたシンテニーブロックが記録される。TSVの列は以下の通り。

  • synteny block ID - 同じIDを持つ行は同じsynteny blockの一部
  • ゲノムファイル名
  • ゲノム染色体/コンティグ
  • ゲノム開始座標
  • ゲノム終了座標
  • 染色体/コンティグ鎖
  • このシンテニブロックにマップされたミニマライザーの数
  • 前のシンテニーブロックとの不連続の理由

(マニュアルより)

 

その他

引用

Multi-genome synteny detection using minimizer graph mappings

Lauren Coombe,  Parham Kazemi,  Johnathan Wong,  Inanc Birol,  René L. Warren

bioRxiv, Posted February 13, 2024.



 

ブルームフィルタを用いて低メモリ使用量且つ高速にsamの重複マークを行う streammd

 

重複テンプレートの同定は、バルクシークエンシング解析における一般的な前処理ステップである。streammdは、Picard MarkDuplicatesの出力を忠実に再現しながら、大幅に高速化し、SAMBLASTERよりはるかに少ないメモリで動作する。streammdは、GitHub https://github.com/delocalizer/streammd からMITライセンスの下で入手可能なC++プログラムである。

 

特徴(レポジトリより)

  • 高速 - デフォルトの設定で、streammdはPicard MarkDuplicatesより5倍高速
  • 大規模なライブラリでも少ないメモリフットプリント - デフォルトの設定で
  • streammd が 1B のテンプレートを処理するのに必要なメモリはわずか 4G 
  • Picard MarkDuplicatesのメトリックスとの高い一致
  • ソフトクリップされたリードを正しく処理
  • ターゲットの偽陽性率を調整可能
  • ストリーミング入出力

インストール

依存

  • Requires c++17 with <charconv>
  • gcc >= 8.1 or clang >= 7 should work.

Github

#リリースよりstreammd-4.3.0.tar.gzをダウロードする
https://github.com/delocalizer/streammd/releases

cd streammd/
./configure
make -j8
sudo make install

$ streammd --help

Usage: streammd [-h] [--input INPUT] [--output OUTPUT] [--fp-rate FP_RATE] [--mem MEM] [--allow-overcapacity] [--metrics METRICS_FILE] [--remove-duplicates] [--show-capacity] [--single] [--strip-previous]

 

Read a SAM file from STDIN, mark duplicates in a single pass and stream processed records to STDOUT. Input must begin with a valid SAM header followed by qname-grouped records. Default log level is 'info' — set to something else (e.g. 'debug') via SPDLOG_LEVEL environment variable.

 

Optional arguments:

  -h, --help                shows help message and exits 

  -v, --version             prints version information and exits 

  --input INPUT             Input file. [default: STDIN] 

  --output OUTPUT           Output file. [default: STDOUT] 

  -p, --fp-rate FP_RATE     The maximum acceptable marginal false-positive rate. [default: 1e-06]

  -m, --mem MEM             Memory allowance for the Bloom filter, e.g "4GiB". Both binary (kiB|MiB|GiB) and decimal (kB|MB|GB) formats are understood. As a result of implementation details, a value that is an exact power of 2 (512MiB, 1GiB, 2GiB etc) gives a modest processing speed advantage (~5%) over neighbouring values. [default: "4GiB"]

  --allow-overcapacity      Warn instead of error when Bloom filter capacity is exceeded. [default: false] 

  --metrics METRICS_FILE    Output metrics file. [default: "streammd-metrics.json"]

  --remove-duplicates       Omit detected duplicates from the output. 

  --show-capacity           Do no work, just print the capacity of the Bloom filter that would be constructed with the given --fp-rate and --mem values. 

  --single                  Accept single-ended reads as input. [default: paired-end] 

  --strip-previous          Unset duplicate flag for any reads that have it set and are no longer considered duplicate. Only ever required if records have previously been through a duplicate marking step. [default: false] 

 

 

テスト

run unit tests

make check

10分くらいかかる。

 

実行方法

STDINからSAMファイルを読み込み、1回のパスで重複をマークし、処理したレコードをSTDOUTにストリームする。

#bwa memと組み合わせる例
bwa mem ref.fa r1.fq r2.fq|streammd > out.sam

 

レポジトリより

  • 60倍のヒトWGS 2x150bpペアエンドシーケンスはn≈6.00E+08テンプレートで構成され、デフォルトの偽陽性率1.00E-06でこれを処理するには、デフォルトのメモリ設定4GiBで十分
  • 入力は有効な SAMのヘッダで始まり、その後に qnameでグループ化されたレコードが続く必要がある。
  • シングルパス処理の性質上、streammdは最初に出会ったテンプレートをオリジナルとして保持し、それ以降に出会ったコピーを光学的要因による複製、あるいはPCRによる複製(optical duplicactes or PCR duplicactes )として、2つを区別せずマークする。
  • 現在の実装ではSAM形式の入力のみを扱う。

引用

streammd: fast low-memory duplicate marking using a Bloom filter 
Conrad Leonard
Bioinformatics, Volume 39, Issue 4, April 2023

 

関連

https://kazumaxneo.hatenablog.com/entry/2017/08/20/203834

 

細菌の近傍に存在するタンパク質ファミリーを調べる ProFaNA

 

 機能的に関連する遺伝子は、特に原核生物において、ゲノム上でしばしば近傍にグループ化されることがよく知られている。この現象が起こる進化的メカニズムは様々であるが、未知の遺伝子の機能を予測するのに利用できる。ここでは、現在利用可能な膨大なゲノムデータを活用した、シンプルで頑健な統計的アプローチを提供する。タンパク質のドメインを機能単位とみなし、問い合わせたドメインのゲノム近傍に有意に多く存在する他の機能単位(ドメイン)を探索することができる。この解析は異なる分類学的レベルにわたって行うことができる。また、非常に近縁な多数の株(例えば病原性株)に焦点を当てることが多いゲノム配列決定プロジェクトによる分類学的空間の不均一なサンプリングを補正するための規定を設けることもできる。この目的のために、亜分類群内の出現頻度を平均化するオプションの手順が用意されている。

いくつかの例から、このアプローチにより、未解決の遺伝子ファミリーに有用な機能予測を提供できること、また、この情報を他のアプローチとどのように組み合わせればよいかがわかる。この手法はウェブサーバー(http://bioinfo.sggw.edu.pl/neighborhood_analysis)として公開されている。

 

Documentation

http://bioinfo.sggw.edu.pl/news/

Usage

http://bioinfo.sggw.edu.pl/news/1/

 

webサービス

http://bioinfo.sggw.edu.pl/neighborhood_analysis/

解析するには最初にアカウントを作ってログインする必要がある。右上から登録できる。

 

解析には3つの方法がある。

1,遺伝子リストを指定 - 1つ目は、クエリは遺伝子リスト(IMG-JGI識別子、1行に1つ、例えば2264882819)を指定する。指定した遺伝子リストの全遺伝子の近傍領域が解析される。

遺伝子リスト、ジョブ名を指定して実行する(画像ではログイン前なのでLogin firstとなっている)。任意で解析するDNA鎖の設定、多重比較検定やその検定テストの閾値を設定できる。

 

example遺伝子リスト(先頭のみ)。

サブミット後は右上のTaskから管理する。ラボの他のツールと一緒に管理できるようになっている。一番上がProFaNA。

example遺伝子リストで試したが、エラーとなった。

 

2,Pfamドメイン名を指定 - クエリはPfamドメインで指定(Pfam ID、例えばPF02696) 

結果を属や科の分類学的レベルで平均化するか、平均化しない(all)で3つから選ぶ。

PFAM Domain as query (genera level)を選択した。

 

Pfam IDと属を指定する。他の項目は1と同様。

解析が完了すると、結果へのリンクがメールで通知される。

 

結果は表として提示される。GeneOntologyのアノテーションと、一致する有意に過剰発現しているドメインが示されている。このクエリでは、Uncharacterized ACR, YdiU/UPF0061 family、Fatty acid desaturase、Transcriptional regulatory protein, C terminalのP値がゼロとなっている。機能未知遺伝子についてはGene Ontologyの列は空白となっている。

CSV形式でダウンロードしてスプレッドシートアプリケーションで閲覧できる、

 

3,タンパク質配列を指定 - 相同なタンパク質の検索にはDIAMONDが使用される。選択されたカットオフ値で得られたヒットが集められ、その近傍が解析される。

タンパク質配列を入力する。

exampleタンパク質配列で試したが、ジョブがキューに登録されなかった。

 

論文より

  • 細菌では、局所的な再配列のためにオペロンとして分類することはできないが、共発現や共機能性を保持していることが多い。酵母、植物、哺乳類などの真核生物ゲノムでも、類似または関連した機能を持つグループへの遺伝子クラスタリングが起こる。
  • ゲノム中の共発現遺伝子の位置はランダムではない。例えば、タンパク質複合体のメンバーやシグナル伝達・代謝経路のエレメントなど、共通の生物学的機能に関与する遺伝子がゲノム内で隣接する傾向があることはよく知られている。したがって、原核生物ゲノム内の遺伝子groupingを調べることは、未特定遺伝子の機能を予測したり、タンパク質相互作用を予測したり、ゲノムが進化してきた事象を解析したりするための出発点となり得る。
  • ゲノム近傍の統計的有意性と生物学的関連性を評価する方法を開発することが課題である。ProFaNAは、クエリドメインをコードする遺伝子のゲノム近傍において、どのタンパク質構造/機能ドメインが有意に多くコードされているかを評価する。
  • ProFaNAはオペロンのような関係を予測することなく、物理的なゲノムの近接性にのみ着目し、偏りのない方法でゲノム近傍を解析している。
  • JGIのIMG/Mの88,754の細菌ゲノムコレクションが使用されている(https://img.jgi.doe.gov/)。
  • 近傍領域のサイズを各方向に5,000bp(ゲノム配列の10kbp)とすると、近傍領域全体で約10個の遺伝子に相当する。ほとんどのオペロン様機能単位は10遺伝子を超えないと予想されるが、このパラメータはユーザが変更することができる。
  • ProFaNAツールは、著者らの知る限り、数千のゲノムを含む原核生物のゲノム近傍を統計的評価を適用しながら大規模に解析できる唯一のツールである。

 

コメント

論文では機能未知の良く見つかるドメインだったり、機能既知のタンパク質へのProFaNAの適用例も示しています。読んでみて下さい。

引用

Protein family neighborhood analyzer—ProFaNA

Bartosz Baranowski and Krzysztof Pawłowski

PeerJ. 2023; 11: e15715. Published online 2023 Jul 21

 

関連

https://kazumaxneo.hatenablog.com/entry/2020/02/27/073000

 

https://kazumaxneo.hatenablog.com/entry/2024/01/18/003341

 

https://kazumaxneo.hatenablog.com/entry/2022/06/19/023905

 

https://kazumaxneo.hatenablog.com/entry/2021/04/07/073000

 

ゲノムからメタコミュニティの幅広いデータに対応したロバストな機能アノテーションを行うツール MetaCerberus

2024/03/5 更新

 

 MetaCerberusは、超並列、高速、低メモリ、スケーラブルなアノテーションツールであり、ゲノムからメタコミュニティにわたる遺伝子機能を推論する。MetaCerberusは、HMM/HMMERベースのツールを低メモリで高速に提供する。KEGG(KO)、COGs、CAZy、FOAM、そしてVOGsやPHROGsなどのウイルスに特化したデータベースを含む主要な公開データベースに対して、シングルゲノムからメタコミュニティまで、スケーラブルな遺伝子解明を提供する。MetaCerberusはシングルノードでeggNOG-mapper v2より1.3倍高速で、HMM/HMMERモードのみを使用した場合、5倍少ないメモリで実行できた。直接比較すると、MetaCerberusはDRAM、Prokka、InterProScanよりもウイルス、ファージ、古細菌ウイルスのアノテーションが優れている。MetaCerberusは、DRAMと比較して、186倍小さいデータベースと63倍少ないメモリで、ドメイン間でより多くのKOをアノテーションする。MetaCerberusは、差分統計ツール(DESeq2やedgeRなど)、パスウェイエンリッチメント(GAGE R)、pathview Rを使用した統計とパスウェイの自動解析のために完全に統合されている。
 MetaCerberusはPythonで書かれており、BSD-3ライセンスの下で配布されている。MetaCerberusのソースコードPython 3と互換性があり、Mac OS XLinuxの両方で動作する::https://github.com/raw-lab/metacerberus。MetaCerberusはbiocondaを使って簡単にインストールすることもできる: mamba create -n metacerberus -c bioconda -c conda-forge metacerberus。

 

 

 

レポジトリより

MetaCerberusは、生のショットガン・メタオミクスシーケンスデータをナリッジに変換する。メタゲノムFunctional Ontology Assignments for Metagenomes(FOAM)である。KEGG、CAZy、VOG/pVOG、PHROG、COGデータベースを生態系全体のメタボローム解析のための隠れマルコフモデル(HMM)を介して汎用的に解析するためのスタートからゴールまでのPythonコードとして機能する。MetaCerberusは、DESeq2/EdgeRを使用した自動統計差分解析、GAGEを使用したパスウェイエンリッチメント解析、Pathview Rを使用したパスウェイの可視化も提供する。

インストール

Apple siliconのMacへの導入もサポートされている(レポジトリ参照)。ここではcondaを使ってubuntu20.04 LTSに導入した。

依存

  • python >= 3.8
  • fastqc

  • fastp

  • porechop

  • bbmap

  • prodigal

  • HMMER

Github

mamba create -n metacerberus python=3.10 -y
conda activate metacerberus
mamba install -c conda-forge -c bioconda metacerberus -y

> metacerberus.py -h

usage: metacerberus.py [-c CONFIG] [--prodigal PRODIGAL] [--fraggenescan FRAGGENESCAN] [--super SUPER] [--prodigalgv PRODIGALGV] [--phanotate PHANOTATE] [--protein PROTEIN] [--rollup ROLLUP] [--illumina | --nanopore | --pacbio]

                       [--setup] [--uninstall] [--dir_out DIR_OUT] [--meta] [--scaffolds] [--minscore MINSCORE] [--evalue EVALUE] [--skip_decon] [--skip_pca] [--cpus CPUS] [--chunker CHUNKER] [--grouped] [--replace] [--keep]

                       [--hmm HMM] [--class CLASS] [--tmpdir TMPDIR] [--version] [-h] [--adapters ADAPTERS] [--qc_seq QC_SEQ]

 

options:

  --illumina            Specifies that the given FASTQ files are from Illumina

  --nanopore            Specifies that the given FASTQ files are from Nanopore

  --pacbio              Specifies that the given FASTQ files are from PacBio

 

Required arguments

At least one sequence is required.

<accepted formats {.fastq .fasta .faa .fna .ffn}>

Example:

> metacerberus.py --prodigal file1.fasta

> metacerberus.py --config file.config

*Note: If a sequence is given in .fastq format, one of --nanopore, --illumina, or --pacbio is required.:

  -c CONFIG, --config CONFIG

                        Path to config file, command line takes priority

  --prodigal PRODIGAL   Prokaryote nucleotide sequence (includes microbes, bacteriophage)

  --fraggenescan FRAGGENESCAN

                        Eukaryote nucleotide sequence (includes other viruses, works all around for everything)

  --super SUPER         Run sequence in both --prodigal and --fraggenescan modes

  --prodigalgv PRODIGALGV

                        Giant virus nucleotide sequence

  --phanotate PHANOTATE

                        Phage sequence

  --protein PROTEIN, --amino PROTEIN

                        Protein Amino Acid sequence

  --rollup ROLLUP       Rolled up annotations from HMMER

 

optional arguments:

  --setup               Set this flag to ensure dependencies are setup [False]

  --uninstall           Set this flag to remove downloaded databases and FragGeneScan+ [False]

  --dir_out DIR_OUT     path to output directory, creates "pipeline" folder. Defaults to current directory. [./results-metacerberus]

  --meta                Metagenomic nucleotide sequences (for prodigal) [False]

  --scaffolds           Sequences are treated as scaffolds [False]

  --minscore MINSCORE   Score cutoff for parsing HMMER results [60]

  --evalue EVALUE       E-value cutoff for parsing HMMER results [1e-09]

  --skip_decon          Skip decontamination step. [False]

  --skip_pca            Skip PCA. [False]

  --cpus CPUS           Number of CPUs to use per task. System will try to detect available CPUs if not specified [Auto Detect]

  --chunker CHUNKER     Split files into smaller chunks, in Megabytes [Disabled by default]

  --grouped             Group multiple fasta files into a single file before processing. When used with chunker can improve speed

  --replace             Flag to replace existing files. [False]

  --keep                Flag to keep temporary files. [False]

  --hmm HMM             Specify a coma seperated list of databases for HMMER. Use quotes around the list, or avoid spaces. (KOFam_all, KOFam_eukaryote, KOFam_prokaryote, COG, CAZy, PHROG, COG) [KOFam_all]

  --class CLASS         path to a tsv file which has class information for the samples. If this file is included scripts will be included to run Pathview in R

  --tmpdir TMPDIR       temp directory for RAY [system tmp dir]

  --version, -v         show the version number and exit

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

 

  --adapters ADAPTERS   FASTA File containing adapter sequences for trimming

  --qc_seq QC_SEQ       FASTA File containing control sequences for decontamination

 

Args that start with '--' can also be set in a config file (specified via -c). Config file syntax allows: key=value, flag=true, stuff=[a,b,c] (for details, see syntax at https://goo.gl/R74nmi). In general, command-line values override

config file values which override defaults.

 

 

データベース

OSFで管理されている。--setupオプションでダウンロード出来る。

https://osf.io/3uz2j/

metacerberus.py --setup
  • --setup    Set this flag to ensure dependencies are setup [False]
  • --uninstall    Set this flag to remove downloaded databases and FragGeneScan+ [False]

ダウンロードには40分ほどかかった。

 

実行方法

ゲノムアセンブリfastaファイル、プロテオームのfastaファイル、Illumina, PacBio, Oxford Nanoporeなどの主要なHTSのfastqフォーマットに対応している。

 

1、fastaファイルを使用 : ゲノムアセンブリを指定してゲノムの機能的アノテーションを行う(.fasta, .fa, .fna, .ffn)。遺伝子予測して得たprotein fastaにも対応。

バクテリアやバクテリオファージなら”--prodigal”でfastaを指定する。真核生物やバクテリオファージ以外のウィルスなら”--fraggenescan”でfastaを指定する。”--hmm”でデータベースを指定する。

#bacteriaの例

#全データベース指定
metacerberus.py --prodigal lambda.fna --hmm "KOFam_all, COG, VOG, PHROG, CAZy" --dir_out lambda_dir

#KOFam_all DBのみ
metacerberus.py --prodigal lambda.fna --hmm "KOFam_all" --dir_out lambda_ko-only_dir

#KOFam_prokaryote DBのみ
metacerberus.py --prodigal ecoli.fna --hmm "KOFam_prokaryote" --dir_out ecoli_ko-only_dir

#KOFam_eukaryote DBのみ
metacerberus.py --fraggenescan human.fna --hmm "KOFam_eukaryote" --dir_out human_ko-only_dir

#Viral/Phage DBのみ
metacerberus.py --prodigal lambda.fna --hmm "VOG, PHROG" --dir_out lambda_vir-only_dir

#複数fasta、全データベース指定(4つ以上あるなら)
metacerberus.py --prodigal fasta_dir/ --hmm "KOFam_all, COG, VOG, PHROG, CAZy" --dir_out lambda_dir
  • --hmm     Specify a coma seperated list of databases for HMMER. Use quotes around the list, or avoid spaces. (KOFam_all, KOFam_eukaryote, KOFam_prokaryote, COG, CAZy, PHROG, COG) [KOFam_all] 
  • --dir_out    path to output directory, creates "pipeline" folder. Defaults to current directory. [./results-metacerberus] 
  • --prodigal    Prokaryote nucleotide sequence (includes microbes, bacteriophage)

 

2、fastqを使用 : fastqを指定して(メタ)コミュニティの機能的プロファイリングを行う(.fastq format)。

リードを指定した場合、リード全ての機能的アノテーションが行われる。

#illlumina fastq

#バクテリアアーキア、バクテリオファージのmetagenomes/metatranscriptomes
metacerberus.py --prodigal [input_folder] --illumina --meta --dir_out [out_folder]

#真核生物とウィルスのmetagenomes/metatranscriptomes
metacerberus.py --fraggenescan [input_folder] --illumina --meta --dir_out [out_folder]

(レポジトリではPacBioとNanoporeのリードを使う例も書かれている。)

 

DRAMより遥かに高速化されているが、それでもHmmerのサーチとフィルタリングにはそれなりの時間がかかる。細菌ゲノム1つを使ってランすると30~40分ほどかかり、細菌ゲノム50個ほど使うと3~4時間かかった(結果をパースしたりsummaryを作るために時間がかかる)。さらに100MBx2のfastqを1組だけ使うと30分ほどかかった(いずれも3990X使用、デフォルトで全CPUが使用される)。デフォルトでクラスタに対応しているので、よりリッチな計算機環境でランすることが推奨される。

 

出力例

--dir_outで出力ディレクトリを指定しなかった場合、results-metacerberus/が作成され、ステップごとにサブディスプレイに結果が保存される。

results-metacerberus/

step_08-hmmer/

step_09-parse/

step_10-visualizeData/prodigal~

sunburst_KEGG.html

クリックする事でインタラクティブに操作可能。

(COGとFOAMもKEGGと同様にプロットされる、metatranscriptomeやmetagenomeでもこのプロットは得られるので、コミュニティの中でエンリッチされてる機能的アノテーションを推測できる)

 

step_10-visualizeData/prodigal~/annotation_summary.tsv

 

step_10-visualizeData/combine/

stats.html

アノテーションされたタンパク質数とタンパク質長、アノテーションデータベースごとのアノテーション数などのstats。

 

4つ以上のサンプルを使ってランした場合、各アノーテーションごとのサンプルの類似性を調べるためのPCA plotなども出力される。step_10-visualizeData/にはサンプルごとのサブディレクトリができ、combineのファイル構造も変化する。

step_10-visualizeData/combine/

 

 

stats.html - サンプルごとの統計。画像はゲノムそれぞれのデータベースごとのアノテーションカウント値。

 

アノテーションそれぞれについて、インタラクティブなPCA plot.htmlも出力される。

 

その他

  • HMMファイルを処理した後、MetaCerberusはKEGG/FOAMからKO(KEGG Orthology)カウントテーブルを計算し、GAGEとPathViewで処理する。パスウェイエンリッチメント解析にはGAGEを、代謝パスウェイの可視化にはPathViewが推奨される。この解析を実行するには、--class オプションで "class "ファイルが必要である。また、このタイプの解析には少なくとも4つのサンプルを使用する必要がある。
  • MetaCerberusはRAYによるマルチプロセッシング/マルチコンピューティングに対応している。ジョブはlocal Ray instanceとしてwebブラウザ上で管理できる。
  • 真核生物では、珪藻でテストされている。しかし、他の真核生物についてもベータテスターが必要。
  • マルチコアの計算機で複数サンプルを提供してランした場合、計算が終わったサンプルからstepごとにサブディレクトリに保存されていく。全部のサンプルの終了までに時間差が生じる。

コメント

これまでのメタコミュニティの機能的アノテーションを行うツールは、インストールでエラーが発生しやすい(snakemakeやnextflowの環境が作れないことも含む)、データベースのダウンロードエラーが起きやすい( NCBIなどがデータのURLを変更するため)、などでトラブルが高確率で起きるのがネックで、使えるようになるまで1日溶けることも珍しくありませんでしたが、MetaCerberusは一度もエラーの発生が無く最後までランできました。非常に使いやすい印象です。注意点として、テストのために1つのワークステーションで並行して2つのジョブを実行すると、後から実行した方のHMMプロファイリングがスキップされる現象が発生しました。ジョブ管理されているはずですが、1つずつランした方が良いかもしれません。

引用

MetaCerberus: distributed highly parallelized HMM-based processing for robust functional annotation across the tree of life 
Jose L Figueroa, Eliza Dhungel, Madeline Bellanger, Cory Brouwer, Richard Allen White, III
Bioinformatics, Published: 29 February 2024

 

関連

https://kazumaxneo.hatenablog.com/entry/2022/02/10/000615

ANIについて

2024/03/04 誤字修正、03/05 引用追加、文章校正

 

このブログでこうゆう話を書くのは珍しいのですが、今日は自分も良く使っているANIについてなるべく分かりやすく説明します。

 

2つの菌のゲノムDNA間を比較するAverage Nucleotide Identity (ANI) 比較は、wet実験のDNA-DNAハイブリダイゼーションのin silico代替としてdDDHと共に広まり、現在ではルーティンに行われるようになりました。このDNA-DNAハイブリダイゼーション(DDH)とは、比較したい2つの菌のDNAを抽出し、800塩基対程度に切断後、一本鎖にしてハイブリダイゼーションさせ、一本鎖に戻る温度を評価してDNA全体がどのくらい似ているかどうか評価することで、2つの菌が同種かどうか判定するために使用される細菌分類のための実験手法です(wikiより)。古くから知見があり確立された実験手法ですが(ref.1)、実験までの準備が大変だったり、再現性良く実験するには熟練の分子生物学の研究者でないと難しかったりと、簡単に使えるものではないようです(ref.2、まず16Sの比較を行い、はっきりしなかった時に最も近縁な既知の菌と未知の株のDNAを使って実行する)。ANIはそのDDHのin silicoの代替としてよく使われています。

遡ると、DDHが報告された半世紀以上の昔は、細菌もウィルスもゲノム配列の決定は夢のような時代で、分類学に巨大な計算機を使う意味も必要性も無かった時代です。しかし商用の全自動DNAシークエンサーが登場することでゲノム解読プロジェクトが始動し、早くも1990年代には最初の完全長細菌ゲノムが報告され、2000年代になると数百の細菌ゲノム配列が決定されていました(ref.3)。2004年あるいは2005年に商用のRoche GS20が市場に登場して(ref.4,  5)サンガー法とは桁違いなハイスループットのシークエンシングが可能になり、計算機の性能向上とin silicoのアセンブル手法の発展とも相まり、細菌の(ドラフト)ゲノムを短い時間で決定することが可能になりました。そうすると、煩雑なウエットの実験にとって代わり、DNAの塩基配列を計算機上で比較して同種かどうか素早く識別できないかと考える研究者が出てくるのは自然な事です。2005年にKonstantinos T. KonstantinidisとJames M. Tiedjeが発表した論文では(ref.6)、70%のDDHが93-94%のANIに対応することから、ANI値94%以上が従来の細菌分類における同種に相当するだろうと報告しています。2007年のJohan Gorisらの論文では(ref.7)、ゲノム配列が得られた28株について総当たりでDDHの値を決定し、種の識別に推奨されるDDH 70 %カットオフポイントはANI値95%に相当すると報告しています。その後のより大規模な調査では、同種か別種かの境界はANI95-96%ANI付近と報告しています(ref.8)。16S rRNAによる種の閾値と似ていて、比較対象が大規模になるにつれて閾値は少しずつ厳しくなっている傾向がありますが、現在の細菌の種という分類が完全に均一な進化距離での分類の最小単位ではない以上(例;大腸菌赤痢菌)、研究に使用した種によって閾値は揺らぐ可能性があり、大まかには95%がカットオフと言って誤解は無いかと思います。

このANI計算では、比較したい2つの細菌のDNA配列を用意し、それぞれバラバラに切って、ベストマッチの断片間の塩基配列同一性を計算します。実際の方法はref.6の論文では詳細が書かれておらず、ref.7のJohan Gorisらの研究のMethodsに書かれています。以下の通りです。

--------------------------------------------------------------------------------------------------

配列に基づく比較 - すべてのペアワイズ、全ゲノム配列比較は以下のように行った。ペアの片方のゲノム配列「クエリ」を連続した1020 ntの断片に切断した。1020ntのカットオフは、DDH実験中にゲノムDNAが約1kbの断片に断片化するのに対応するように使用した。異なるカットオフ(例えば、より小さい断片)を用いても、結果は特に変化しなかった(データは示さず)。次に、1020 ntの断片を、blastnアルゴリズム(Altschul et al., 1997)を用いて、ペアのもう一方のゲノム配列(「リファレンス」)全体に対して検索した。blastnアルゴリズムは以下の設定で実行した: X=150(Xはギャップアライメントのドロップオフ値)、q=-1(qはヌクレオチドミスマッチのペナルティ)、F=F(Fは繰り返し配列のフィルター);残りのパラメーターはデフォルト設定で使用した。これらの設定は、より遠縁のゲノムを比較する場合、デフォルトの設定よりも感度が良くなる。

(1パラグラフ省略)

クエリゲノムとリファレンスゲノムの間のANIは、その長さの少なくとも70%のアライメント可能な領域にわたって、30%以上の全配列同一性(全配列に沿った同一性に再計算)を示したすべてのblastnマッチの平均同一性として計算された。このカットオフは、アライメントされた配列間の類似度が低いため、相同性の推測に誤差が生じやすい類似性検索の「トワイライトゾーン」を超えている(Rost, 1999; Sander & Schneider, 1991)。したがって、計算では相同なDNA断片のみを考慮したと仮定できる。

逆検索、すなわちリファレンスゲノムをクエリとして使用する検索も、相互値を提供するために実行した。全ゲノム配列ファイルからの1020 nt断片の抽出、Blast検索用のデータベースのフォーマット、Blast出力の自動解析にはPerlスクリプトを使用した。”

--------------------------------------------------------------------------------------------------

https://www.microbiologyresearch.org/content/journal/ijsem/10.1099/ijs.0.64483-0#tab2より翻訳して引用

 

このように書かれています。整理すると、1)クエリゲノムをおよそ1kbに断片化し、BLASTnプログラムを使ってターゲットゲノムのベストヒットを探索して記録する。2)クエリとターゲットを入れ替えて逆検索する。3)順方向と逆方向の検索によって計算されたベストマッチの相互の平均値をANI値とする、となります。両方向から検索することでパラログへのアラインメントの影響を軽減出来ます。

2009年にはANI計算をスイス・アーミーナイフ的に使いやすくしたJavaアプリケーション - JSpeciesが公開され(download  link)(ref.9)、GUIアプリからANI計算が出来るようになりました。2015年にはJSpeciesのweb版であるJSpeciesWSも公開され、web上でANI計算が出来るようになりました(ref.10)。2016年にはreciprocal blastnを使うOrthoANIが発表され(ref.11)(*1)、2019年には、MASHのMinHash技術を応用したMashmapを使うことで、近似であるものの2桁あるいは3桁以上高速なFastANI(ref.12)が報告されるなど、近年はよりユーザーフレンドリーに、そしてより大規模な比較でもスケールできるように発展が続いています。dRepやGTDB-tkのように種の判定にANIを利用するツールも増えています。

ANI比較によって2つの菌が同種なのか新種相当の距離があるか1つの証拠、あるいは傍証を得る事ができます。比較したい菌が3つ以上ある時は全てのペア間の総当たりでANI計算を行い、比較した菌間で同種に相当するグループが存在するのか調べる事もできます。Pyaniは、ANI比較を総当たりで行い、階層的クラスタリングしたヒートマップで視覚化する人気のpython実装です(pyani、star数175)。Pyaniを使用したと思われる論文は良く見ます。

 

このANI値ですが、クローンな株間でANI比較をしたとしても、塩基置換が1か所でも発生していればANI値は塩基置換の数だけ敏感に低下します。indel変異が発生していても敏感に低下します。これを確認するため、1995年に完全長ゲノムが解読されたHaemophilus influenzae Rd株(ref.13)を使用して、ランダムに塩基置換を発生させてみます(FTP)。

for i in {1..10}
do
msbar -sequence GCA_000027305.1.fna -outseq "output_${i}.fna" -count 10 -point 0 -block 0 -codon 4
done

emboss msbar(紹介)(ref.14)を使って塩基置換を10回ランダムに発生させ、それを10回ループ実行しています。

得られた10個のゲノムと元のゲノム1個について、pyaniで総当たりのANI計算をします。

average_nucleotide_identity.py -m ANIb -g -i ./ -o ANIb_outdir

結果の表です。太字が元のゲノムを含む列と行です。pyaniは逆方向の検索をしないのでクエリとターゲットがどちらかで値が微妙に変化します。

H. influenzae Rd株株のゲノムサイズは1,830,138塩基対なので、ゲノムに沿って10か所で塩基置換が発生すると、0.999994536...となります(断片化せずに2ゲノム間の配列同一性を計算すると)。表の数値はこの値と桁があっていて、上手くいってそうです。

 

同様に小さなindelが発生した時も調べます。まずランダムに1塩基挿入(-point 2)を10か所で発生させます。10回行います。

#insertion
for i in {1..10}
do
msbar -sequence GCA_000027305.1.fna -outseq "In_${i}.fna" -count 10 -point 2 -block 0 -codon 0
done

pyaniの結果です。

今度は全てのペアのANI値は完全に一致していて、0.99998853...となりました。

 

次はランダムに1塩基欠失(-point 3)を10か所で発生させます。

#deletion
for i in {1..10}
do
msbar -sequence GCA_000027305.1.fna -outseq "output_${i}.fna" -count 10 -point 3 -block 0 -codon 0
done

pyaniの結果です。

1塩基挿入の時と同じく、全てのペアのANI値は完全に一致して0.99998853...となりました。ANIはsmall indelの発生数に応じて値が減少します。

 

大きな挿入として、サイズが1,000から10,000塩基の範囲でランダム塩基の挿入を、ゲノムのランダムな10か所で発生させます。10回行います。

#insertion
for i in {1..10}
do
msbar -sequence GCA_000027305.1.fna -outseq "In_${i}.fna" -count 10 -point 0 -block 2 -codon 0 -minimum 1000 -maximum 10000
done

pyaniの結果です。

セルフヒットを除いてANIの全ての平均は

0.99999693(99.9996%)でした。1塩基挿入を発生させた時のANI値の全平均の

0.99998953(99.9989%)より高い値です。違和感を感じないでしょうか?

 

より挿入のサイズを大きくしてみます。10kbから100kbの範囲で10回。出力はLarge_In_Mutant配列と呼びます。

msbar -sequence CA_000027305.1.fna -outseq Large_In_Mutant.fna -count 10 -point 0 -block 2 -codon 0 -minimum 10000 -maximum 100000

リファレンス配列とLarge_In_Mutant配列のANI値をorthoaniで計算させると0.9995440157480315となりました。大きな挿入が組み込まれているのですがANI値は高いままですね。

dot plotでもリファレンス配列とLarge_In_Mutant配列を比較してみます。D-Genies(ref.15)を使用しています。

y軸がリファレンスゲノムです。リファレンスゲノム配列の7箇所で、リファレンスにはなくてLarge_In_Mutant配列にだけ存在する領域があることが分かります。

 

以下はartemis synteny plot(ref.16)(紹介)を使ったシンテニープロット図です。上がリファレンスゲノム、下がLarge_In_Mutant配列です。Large_In_Mutant配列にリファレンスゲノムには存在しない長い配列が存在している事がより分かりやすく視覚化されています。

0.9995..というANI値はこの挿入がある無しのゲノムを比較していることを考えると、高すぎると思うでしょう。しかしANIの定義的に数値は間違っていません。後で説明します。

 

最後にmsbarを使ってサイズが1,000から10,000塩基の範囲でランダムなサイズの欠失をゲノムのランダムな10か所で発生させました(10回独立に実行)。ANI比較結果です。

セルフヒットを除いてANIの全ての全平均は0.99993924(99.994%)でした。これは大きな挿入変異の導入時の平均よりも低く、塩基置換で例えればゲノムに沿って111か所で塩基置換が発生させた時のANI値に相当する値です。1kbから10kbの大きな欠失をゲノムの10か所から発生させているにも関わらず、感覚よりもずっと高いANI値と感じませんでしょうか?

 

なぜこうなるかというと、発生させた大きな挿入配列が元のゲノムにはないランダムな配列だからです。この挿入を含むDNA断片は挿入がない方のゲノムのどのDNA断片ともマッチせず(ANIは70%以下は考慮しない)、ANI値には反映されません。

大きな挿入の例として、系統的に近い菌株間で起きやすいとされる遺伝子水平伝播(HGT)を例にイラストを書いています。ref.BがHGTで獲得したDNA断片は、ref.Aのゲノムのどの断片にもヒットしません(長さが違うように書いているのは誤解を招くかもしれませんが)。

 

続いて欠失です。ref.Bのfragment 2に相当するDNA領域が欠失すると、ref.Aのfragment 2はref.BのどのDNA断片とも一致しません。ref.Bをクエリとする場合は、fragment 2のDNA断片はそもそも存在しません。結果、ANI値の計算では欠失は全く数値に反映されない可能性があります。

 

欠失と挿入の説明から分かるように、ANIではサイズが大きくなると変異を考慮して値を出すことが難しいということです。これは、ほかの種類の構造変異ではより悪化する可能性があります。

 

下の図は逆位の例です。fragment3の逆位のブレイクポイントが断片のサイズに一致した場合、ANI値は100%のまま変化しないことがあり得ます。

 

転座の例です。ゲノム上のfragment2と3の位置がシャッフルされたとしても、ANI値は100%のまま変化しないことがあり得ます。

 

まとめ

この話の教訓は、ANI比較は万能なペアゲノム間の距離比較メトリクスではないということです。むしろ構造変異に対しては脆弱です。この話を現実的なシナリオで例えると、例えばクローンな2つの菌の片方だけがプラスミドを獲得した時、ANI比較ではプラスミドを獲得した菌としていない方の菌でANI比較しても100%のANIを示す可能性があるということです(コア遺伝子の系統推定も同様;片方にしかない遺伝子は比較対象にならない)。ANI値が100%を示す時、2つの菌は完全にクローンな菌である可能性もあるわけですが、断定するのは危険な可能性もあると言えます。ref.6の論文では、ヒトとチンパンジー間のANIは98.7%とありますが、これも一部の比較可能なDNA領域だけでANI値が算出されているからなのかもしれません(このような可能性があるので、JGIのANI計算ではアラインメントfractionも考慮する)。ゲノムの大きな変化を調べるには、シンテニー比較やドットプロット解析が有効です。また、完全長ゲノムが得られている菌ならゲノムサイズを比較するだけでクローンかどうか簡易判定できるでしょう(同じになるはず)。

ANIは同種か別種かを判定する基準として使用されています。このような文脈で、ANI比較も参考にしつつ、多相的なアプローチによって同種か別種か総合的に解釈するのは良いANIの使い方だと思います。一方で、ANIはユークリッド距離での進化距離に相当すると考えて拡大解釈するとゲノムの重要な変化を見逃す可能性があり、これはANIの落とし穴と言えます。特にANIがblastnやmummerに依存している以上、これらのプログラムではギャップとして考慮できない大きなサイズの変化に脆弱です。ANIに限らず万能なゲノム比較メトリクスは無いということに注意してください。

 

コメント

web上にANIについてまとめた情報が無かったので簡単に説明してみました。間違いがあれば教えて下さい。

 

引用および参考文献

1 Chapter15 - DNA–DNA Hybridization. Ramon Rosselló-Móra, Mercedes Urdiain, Arantxa López-López. Methods in Microbiology Volume 38, 2011, Pages 325-347

2 細菌の系統分類と同定方法. 河村好章. 日本細菌学雑誌55 (3): 545-584, 2000

3 Insights from 20 years of bacterial genome sequencing
Miriam Land, Loren Hauser, Se-Ran Jun, Intawat Nookaew, Michael R. Leuze, Tae-Hyuk Ahn, Tatiana Karpinets, Ole Lund, Guruprased Kora, Trudy Wassenaar, Suresh Poudel, and David W. Usser. Funct Integr Genomics. 2015; 15(2): 141–161.

4 Generations of Sequencing Technologies: From First to Next Generation. Mehdi Kchouk, Jean-François Gibrat and Mourad Elloumi.  Biol Med (Aligarh) 2017, 9:3

5 次世代シーケンス技術の現状と今後―2020. 中村昇太. 生物工学会誌第99巻第5号 242–245. 2021

6 Genomic insights that advance the species definition for prokaryotes. Konstantinos T. Konstantinidis and James M. Tiedje. Proc Natl Acad Sci U S A. 2005 Feb 15; 102(7): 2567–2572.

7 DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Johan Goris, Konstantinos T Konstantinidis, Joel A Klappenbach, Tom Coenye, Peter Vandamme, James M Tiedje. Int J Syst Evol Microbiol
. 2007 Jan;57(Pt 1):81-91.

8 Towards a taxonomic coherence between average nucleotide identity and 16S rRNA gene sequence similarity for species demarcation of prokaryotes. Mincheol Kim, Hyun-Seok Oh, Sang-Cheol Park and Jongsik Chun. Int J Syst Evol Microbiol
. 2014 Feb;64(Pt 2):346-351.

9 Shifting the genomic gold standard for the prokaryotic species definition. Michael Richter, Ramon Rosselló-Móra. Proc Natl Acad Sci USA. 2009 Nov 10;106(45):19126-31.

10 JSpeciesWS: a web server for prokaryotic species circumscription based on pairwise genome comparison . Michael Richter, Ramon Rosselló-Móra, Frank Oliver Glöckner, Jörg Peplies. Bioinformatics, Volume 32, Issue 6, March 2016, Pages 929–931.

11 OrthoANI: An improved algorithm and software for calculating average nucleotide identity Free. Imchang Lee​, Yeong Ouk Kim​, Sang-Cheol Park,​ and Jongsik Chun. Int J Syst Evol Microbiol. 2016 Feb;66(2):1100-1103.

12 High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Chirag Jain, Luis M Rodriguez-R, Adam M Phillippy, Konstantinos T Konstantinidis, Srinivas Aluru. Nat Commun. 2018 Nov 30;9(1):5114.

13 Whole-Genome Random Sequencing and Assembly of Haemophilus influenzae Rd. R D Fleischmann  1 , M D Adams, O White, R A Clayton, E F Kirkness, A R Kerlavage, C J Bult, J F Tomb, B A Dougherty, J M Merrick, et al(35名). Science. 1995 Jul 28;269(5223):496-512.

14 EMBOSS: The European Molecular Biology Open Software Suite. Peter Rice a, Ian Longden a, Alan Bleasby. Volume 16, Issue 6, 1 June 2000, Pages 276-277.

15 D-GENIES : Dot plot large GENomes in an interactive, efficient and simple way. Floréal Cabanettes, Christophe Klopp​. PeerJ. 2018; 6: e4958.

16 ACT: the Artemis Comparison Tool. Carver TJ, Rutherford KM, Berriman M, Rajandream MA, Barrell BG, Parkhill J. Bioinformatics. 2005 Aug 15;21(16):3422-3. Epub 2005 Jun 23.

 

*1

オリジナルのANIとOrthoANIの主な違いは論文で以下の通りと説明されている; 1) OrthoANIでは、両方のゲノムがin silicoで断片化される、2) OrthoANIでは1020bp未満の断片は使用しない、3) OrthoANIでは、BLASTnプログラムを用いて2つの断片がベストヒットとして相互に(reciprocally)検索された場合にのみ、その同一性値がその後の計算に含まれる。