2024/09/26 推敲、9/27 追記
前回のパンゲノム解析の説明の続きとなります。
前回: https://kazumaxneo.hatenablog.com/entry/2024/08/28/163036
1、パンゲノムプロット
パンゲノム解析結果を視覚化する最も一般的な方法は、x軸にゲノム数を、縦軸にコア遺伝子（ファミリー）数もしくはパンゲノム遺伝子（ファミリー）数の累積数をとってグラフにプロットする方法です。この"パンゲノム曲線"を見れば、分析した株間でユニークな遺伝子（ファミリー）数がどのくらいあるのか視覚的に推定することができます。
図１
上の図はRoaryを使ってVibrio属のVibrio harveyi（n=51）のパンゲノムサイズを計算した時の結果です。Roaryのレポートファイル中のグラフをベクトルグラフィック編集ソフトでフォントだけ修正して示しています。図からパンゲノム遺伝子数がゲノムの増加に応じて顕著に増加していることが分かります。プロットに四分位範囲とヒゲがあるのは、RoaryのRscriptがゲノムをランダムにサンプリングして累積遺伝子数の計算が繰り返されているためです。このようなグラフは、どの程度のゲノムのサンプリングで観測対象の種の遺伝子数が飽和するのか、あるいはしないのかなど調べる目的があり、このような累積曲線はレアファクションカーブと呼ばれます。ユニーク遺伝子のみとコア遺伝子のみでも同じ形式の図をプロットすることができ、パンゲノム解析で最も一般的なグラフの１つとなっています。
２、株間の異質性を探るためのU-shape plot
もう１つパンゲノム解析でよく見かけるグラフは、ゲノム間で共有されている遺伝子数を示したヒストグラムのグラフです。横軸にはあるbinサイズでの共有されているゲノム数、縦軸には頻度をプロットします。どのようなグラフになるのか例を示すため、NCBIの分類でShigella flexneriとなっている完全長の15株、NCBIの分類で大腸菌となっている完全長の20株のゲノムをそれぞれダウンロードしてパンゲノム解析ソフトウェアのPPanGGOLiN：Partitioned PanGenome Graph Of Linked Neighbors（Guillaume Gautreau et al, 2020）（紹介）を使ってまとめて解析してみます。
PPanGGOLiN（紹介）のallワークフローでは、アノテーションから始まりパンゲノム解析の全フローを自動で実行してくれます。使うには、ゲノムのリストを指定します；１列目に名前（特殊文字やホワイトスペースなしの株名など）、２列目にfastaファイルのパス（カレントにあるならそのままファイル名）のTSVファイルです。このTSVを指定して実行します。
#30CPJU指定
ppanggolin all --fasta list.tsv -c 30 -o outdir
PPanGGOLiNは計算速度が素晴らしく高速です。今回の35株の解析ではアノテーションを含めてわずか数分で計算が終わりました（3990X、30CPU指定）*2。
出力ファイル中に遺伝子ファミリー数を示したヒストグラムが保存されます（図2）。例えば10株中全ての組み合わせの5株だけに共有されている遺伝子（ファミリー）が20個あれば、x=5の位置には20の高さの棒グラフがプロットされます。
Ushaped_plot.html
図２ S. flexneriの15株とE. coli 20株の遺伝子共有性を表す頻度プロット。横軸は株数、縦軸は共有されている遺伝子（ファミリー）数をそれぞれ表す。
図2の棒グラフは株間で共有されている割合によって３つの色：cloud、shell、persistentで色付けされています。cloudは1個もしくはごく少数の個体のゲノムでしか見つからない遺伝子ファミリーで、PPanGGOLiNのwikiではファージ関連遺伝子が多いと書かれています。shellは中程度の個体に共有されている遺伝子ファミリーで、しばしば遺伝子の水平伝播によって獲得された遺伝子、環境適応、病原性、病原性、あるいは二次代謝産物の生産に関連する機能がコードされていると書かれています。persistentはその種の個体ほぼ全てに共通に見られる遺伝子を持つ遺伝子ファミリーが含まれています。persistentの中でx軸の最大値はパンゲノム解析に供したゲノム数に等しく、従って全ての株が持つコア遺伝子の数がプロットされています。x=1には１株にしか見つからないユニークな遺伝子の数がプロットされています。縦の破線はソフトコアの閾値95%を表します。
さて、このヒストグラムですがU字の形状になっています（PPanGGOLiNではU-shape plotと呼ぶ）。このグラフの形状はたまたま生じたわけではなく、解析対象の群間の均質性と遺伝子の保存性が反映された再現性のあるものです。論文では遺伝子頻度スペクトラム（gene frequency spectrum）と呼ばれ、Marie Touchonらの大腸菌パンゲノム解析で遺伝子の共有性とコードしているタンパク質のタイプを議論するために用いられ（Marie Touchon et al, 2009）、2012年のFranz Baumdickerら（Franz Baumdicker et al, 2012）の研究とR. Eric Collinsらの研究（R. Eric Collins et al, 2012）では、この遺伝子頻度スペクトラムのグラフを提示し、パンゲノムをモデル化するデータの一部として使用されています。2013年にはEvgeny N. Gordienkoらが大腸菌、赤痢菌、サルモネラ菌の複数ゲノムを解析し、”よく知られたU字型プロットから、大腸菌ゲノムには存在しない赤痢菌特異的な遺伝子が事実上存在しない一方、 S. entericaと大腸菌と赤痢菌に限定された遺伝子のピークが存在する"、ことを報告しています（Evgeny N. Gordienko et al, 2013）。2018年にMikhail A. Moldovanらはこれらの結果を整理し、遺伝子頻度スペクトラムで観察される異質性の有無を（Mikhail A. Moldovan et al, 2018）客観的な種分類の指標の１つとして使うことを提案しています。Mikhail A. Moldovanらの論文の導入で以下の通り書かれています。
”パンゲノムを表現する便利な方法は遺伝子頻度スペクトラム関数G(k)を考えることであり、これはちょうどk個のゲノムからの遺伝子を含むorthologous gene groups (OGGs)の数として定義される(Baumdicker et al., 2012; Collins and Higgs, 2012)(図1A-C)。通常、同じ生物種に属する株のサンプルを考慮すると、パンゲノムのスペクトル関数は滑らかで、ノイズと区別されるような内側のピークを持たないUのような形状をしている（Gordienko et al）。 ～ しかし、少数でも複数の生物種から得られた混合サンプルの場合、スペクトル関数には内部ピークが生じる（図1C）。 U字型のスペクトル関数を持つゲノムの集合をhomogenous（均質）と呼び、内部ピークを持つ集合をnon-homogenous（不均質）と呼ぶことにする。理論的には、最初はhomogenousであったパンゲノムが進化の過程で非均質になるには、2つの基本的なシナリオがある。〜”
内部ピーク（スパイク）は均質ではない菌株間で遺伝子頻度スペクトラム関数を算出した時に生じます。その内部ピークに相当する遺伝子セットには、表現型の違いを説明する遺伝子セットが含まれている可能性があります。しかし図２をみると、大腸菌のx=20、S. flexneriのx=15の位置に顕著なスパイクは観察されません（PPanGGOLiNのデフォルト設定の結果）。従って、この結果からは２つの種に顕著な表現型の違いを引き起こすまとまった遺伝子セットは存在しないことが予想されます。Evgeny N. Gordienko et al, 2013の研究と再現性がある結果と言えます。
Mikhail A. Moldovanらは、同じ論文で、菌株セットの均質性に基づく細菌種を厳密に定義するための新しい基準を提案しています（Mikhail A. Moldovan et al, 2018）。
"単系統性の種は、(1)配列に基づく系統樹においてmonophyletic（単系統）でなければならず、(2)均質な菌株集合から構成されなければならず、(3)条件1と2を満たす菌株の最大集合でなければならない。弱い定義では、種はmonophyleticかpolymorphic（多系統）のいずれかであり、条件2を満たす系統の最大集合であることが必要である。”
と書かれています。解析対象のゲノム間でパンゲノム解析を実行し、U字状ヒストグラムで遺伝子頻度を視覚化することで、特定の株間にだけ強い異質性が存在するか客観的に評価することができます。提案された指標では、閾値の厳密さによって種やもっと狭い操作上の分類単位で分類することも可能にもなります。この客観的な細菌分類法は、これまで数十以上提案されてきた細菌分類の指標の１つとなっています。
３、パンゲノムグラフ
さて、細菌株間の異質性ですが、他の視覚化方法でも確認することができます。どのような方法が考えられるのでしょうか。１つはパンゲノムグラフを可視化して探索するやり方です。試してみます。
10kbの短い模擬ゲノム1、模擬ゲノム1のORFの１つを欠落させた模擬ゲノム2、模擬ゲノム1の3'末端近くに２つのORFをコードするゲノム領域を挿入させたゲノムを用意しました。３つのゲノムのORFのシンテニーをclinker（紹介）を使って視覚化すると以下のインタラクティブな図が出力されます。
#prodigalの実行。短いのでmetaをつける
prodigal -i strain1.fna -o strain1.gbk -p meta
prodigal -i strain2.fna -o strain2.gbk -p meta
prodigal -i strain3.fna -o strain3.gbk -p meta
#clinkerの実行
clinker *gbk -p putput.html
基準の模擬ゲノム1を真ん中に配置しています。同じ色は同じORFを表します。
今度はこの模擬ゲノム３つをPPanGGOLiNを使ってパンゲノム解析します。得られたpangenomeGraph.gexf *3をgephiで可視化します。ノード（遺伝子）にはclinkerと同じ色を手動で割り当てました。
拡大します。下の図では、３つの模擬ゲノム全てにコードされているORFは最も大きいノードでプロットされ、部分的なゲノムでのみコードされているORFはより小さなノードでプロットされています。辺（エッジ）は隣接したORF間であることを示し、３つのゲノム全てで隣接したORF（ノード）間のエッジは最も太く表現されています。コア遺伝子の順番が維持したまま保存されているゲノム領域と、株によってバリエーションが存在する領域が存在することが示唆されます。
グラフの見方が確認できたところで、さきほどのS. flexneriの15株とE. coli 20株のパンゲノムグラフを見てみます（読み込み手順はgephiで紹介してます）。
ノードは共有されているゲノムの割合でアサインしています。コア遺伝子は濃い赤色になります。
基本的にコア遺伝子はどの株でも連続してコードされていることが分かります。ですが、ところどころで完全には保存されていない遺伝子が存在することが分かります。一か所拡大してみます。
選択後、新しいワークスペースにコピーし、expansionレイアウトで少しノード間の距離を開けました。
辺の色をweightにして色合いも強調します。
赤色のノードはコア遺伝子、緑のノードは非コア遺伝子です。緑は非コアですが、ノードの大きさは中程度で、色も緑であることから、一定の数の株で保存されていることになります。また辺が太めで単純な接続であることから、これらのORF群は連続したORFとしてコードされていると予想されます。遺伝子の並びが強く保存されている場合、遺伝子クラスターとして特定の（２次）代謝系に関与しているか、オペロンとして共転写している可能性も示唆されます。
PPanGGOLiNの著者らは、Regions of Genome Plasticity (RGP) 、すなわちゲノムの可塑性領域を予測するPPanGGOLiNの実装をPPanGGOLiNの論文とは別に論文発表しています（Adelme Bazin et al, 2020）（紹介）*4。パンゲノム解析から遺伝子頻度スペクトラムやパンゲノムグラフの分析を経て集団の均質性を調べることで、特定の亜集団にのみ存在する遺伝子群を特徴づけることが可能になります*5。U字状ヒストグラムにプロットするのは異質性があるかどうかの傾向を確認するために適しており、一方グラフでも全体の傾向を見れますが、グラフは特に異質性が強いゲノム領域があるのかどうか探るためにも適していると思われます。
２回に分けてパンゲノムの解説をしました。何か間違ったことを書いていれば教えて下さい。
引用
PPanGGOLiN: Depicting microbial diversity via a partitioned pangenome graph
Guillaume Gautreau, Adelme Bazin, Mathieu Gachet, Rémi Planel, Laura Burlot, Mathieu Dubois, Amandine Perrin, Claudine Médigue, Alexandra Calteau, Stéphane Cruveiller, Catherine Matias, Christophe Ambroise, Eduardo P. C. Rocha, David Vallenet
PLoS Comput Biol. 2020 Mar; 16(3): e1007732
panRGP: a pangenome-based method to predict genomic islands and explore their diversity
Adelme Bazin, Guillaume Gautreau, Claudine Médigue, David Vallenet, Alexandra Calteau
Bioinformatics. 2020 Dec 30;36(Suppl_2):i651-i658
Clinker & clustermap.js: Automatic generation of gene cluster comparison figures
Cameron L M Gilchrist, Yit-Heng Chooi
Bioinformatics, Published: 18 January 2021
*1 PPanGGOLiNは種や亜種レベルでは問題なく動作するようにチューニングされているが、より多系統のゲノムでPPanGGOLiNを使用する場合、パラメータのチューニングが必要となるとPPanGGOLiNのWikiに書かれている。
*2 PPanGGOLiNは非常に効率的な実装になっており、wikiでは、最大20656ゲノムを扱い、その計算には約1日と120GのRAMを必要としたと書かれている。多くのパンゲノムソフトではこのようにスケールしない。
https://ppanggolin.readthedocs.io/en/latest/user/practicalInformation.html
*3 PPanGGOLiNは2種類のパンゲノムグラフをXML形式のグラフフォーマットであるGEXF形式で出力する（参考）。軽量なlight.gexfファイルのノードには、遺伝子ファミリーに属する遺伝子数、最も一般的な遺伝子名（アノテーションを提供した場合）、最も一般的なプロダクト名（GFFまたはGBFF入力ファイルにアノテーションを提供した場合）、属するパーティション、ヌクレオチド単位の平均サイズと中央値、この遺伝子ファミリーを持つゲノム数が含まれている。エッジにはパンゲノム中に存在する回数が含まれる。lightではない方の.gexfのファイルにはこれに加えて、各遺伝子ファミリーに属する遺伝子、その名前、産物文字列、サイズ、そしてエッジを通して記述される各遺伝子ペアの近傍関係に関するすべての情報が含まれる。
https://ppanggolin.readthedocs.io/en/latest/user/PangenomeAnalyses/pangenomeGraphOut.html
*4 PPanGGOLiN：Partitioned PanGenome Graph Of Linked Neighborsは遺伝子レパートリーの変化をグラフとして表現するため、グラフからバリエーションがある領域やその挿入スポット、保存されたモジュールにおけるセグメンテーションを同定することが可能になっている。しばらく前にメジャーアップデートと思われる（その前はv1.2.105）v.2が出ている（link）。重くなりがちだったHDF-5ファイルの大幅な軽量化が大きな改良点としてアナウンスされている。
*5 見つかった遺伝子群は単系統かもしれないし、病原性因子のような多系統な遺伝子群かもしれない。このブログではグラフのノードのサイズと色のどちらにも共有されている遺伝子数を反映させたが、どちらかをメタデータに反映させれば、例えば病原性の株とそうでない株の遺伝子で色付けすれば、表現型と関係する遺伝子群がより明確に分かるかもしれない（別記事で予定）。
*6 PATOではアクセサリ遺伝子を視覚化している。
https://github.com/irycisBioinfo/PATO/wiki/Manual
