2024/03/04 誤字修正、03/05 引用追加、文章校正
このブログでこうゆう話を書くのは珍しいのですが、今日は自分も良く使っているANIについてなるべく分かりやすく説明します。
２つの菌のゲノムDNA間を比較するAverage Nucleotide Identity (ANI) 比較は、wet実験のDNA-DNAハイブリダイゼーションのin silico代替としてdDDHと共に広まり、現在ではルーティンに行われるようになりました。このDNA-DNAハイブリダイゼーション（DDH）とは、比較したい２つの菌の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は繰り返し配列のフィルター）；残りのパラメーターはデフォルト設定で使用した。これらの設定は、より遠縁のゲノムを比較する場合、デフォルトの設定よりも感度が良くなる。
（１パラグラフ省略）
クエリゲノムとリファレンスゲノムの間の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より翻訳して引用
このように書かれています。整理すると、１）クエリゲノムをおよそ1kbに断片化し、BLASTnプログラムを使ってターゲットゲノムのベストヒットを探索して記録する。２）クエリとターゲットを入れ替えて逆検索する。３）順方向と逆方向の検索によって計算されたベストマッチの相互の平均値を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つの菌が同種なのか新種相当の距離があるか１つの証拠、あるいは傍証を得る事ができます。比較したい菌が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...となります（断片化せずに２ゲノム間の配列同一性を計算すると）。表の数値はこの値と桁があっていて、上手くいってそうです。
同様に小さな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軸がリファレンスゲノムです。リファレンスゲノム配列の７箇所で、リファレンスにはなくて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比較は万能なペアゲノム間の距離比較メトリクスではないということです。むしろ構造変異に対しては脆弱です。この話を現実的なシナリオで例えると、例えばクローンな２つの菌の片方だけがプラスミドを獲得した時、ANI比較ではプラスミドを獲得した菌としていない方の菌でANI比較しても100%のANIを示す可能性があるということです（コア遺伝子の系統推定も同様；片方にしかない遺伝子は比較対象にならない）。ANI値が100％を示す時、２つの菌は完全にクローンな菌である可能性もあるわけですが、断定するのは危険な可能性もあると言えます。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）検索された場合にのみ、その同一性値がその後の計算に含まれる。
