2024/08/29 推敲, 8/30追記, 9/3 誤字修正、補足*7追加, 9/10 画像追加, 9/20追記
今日はパンゲノム（pan-genome）解析について簡単に紹介します。パンゲノム解析とは、解析対象の種内の全ての株に存在する遺伝子や、一部の株にしか存在しない遺伝子のレパートリーを調べ、その共通や差分から生物種や特定の株を特徴づける研究手法です。パンゲノムとはTettelin らの提案した造語で、特定の分類群に属する全ての株が持つ遺伝子の全体を指します。遺伝子プールとか、全ての遺伝子が入った袋ともみなされます（Tim Kahlke et al, 2012）。しかし、パンゲノムとは通常のゲノムとは何が違うのでしょうか？同じ種であれば、系統が違う株間でゲノムの多少の差異はあったとしても、遺伝子セットはほぼ同一ではないのでしょうか？
パンゲノムとは
（原核生物の）パンゲノムの先駆的な研究として引用されることが多いのは、Welchらが行った大腸菌の3株のゲノム比較研究です（Welch et al, 2002）。Welchらは、尿路病原性大腸菌CFT073株のゲノムを決定し、既に解読されている2株の大腸菌、腸管出血性大腸菌EDL933株、実験室株MG1655株（K-12）のゲノムにコードされた遺伝子と比較して、人への病原性に関わってる遺伝子座を探索しました。その中で、3株すべてのゲノムに共通する遺伝子は全体の39％しかなかったことを報告しました。つまり、同じ種でも、株ごとに遺伝子レパートリーにはかなりの差が見られたということです。しかしゲノムが少ないため限られた情報しか得られませんでした。
3年後、Herve Tettelinらは、Streptococcus agalactiae（B群連鎖球菌）の異なる血清型を対象に、完全長のゲノム配列（1株）とドラフトゲノム配列（5株）を解読し、既存の2株のゲノムも加えた合計8株のゲノムを比較することで、S. agalactiaeの遺伝子レパートリーを調査しました。目的はゲノム全体をモデル化することで種を説明できる遺伝子セットを明らかにすることです（Hervé Tettelin et al, 2005）。Tettelinらの論文中には、驚くべき発見と強調して次の記述があります。
"group B Streptococcus （GBS）のコアゲノムは最小1,806遺伝子（95%信頼区間=1,750-1,841）に達し、さらに多くのゲノムが追加されても、比較的一定に保たれることが示された。〜驚くべきことに、8つのゲノムの塩基配列を決定した後もユニークな遺伝子が検出され（論文図３）、数学的外挿すると、さらに多くの株の塩基配列を決定した後も新しい遺伝子が発見されると予測される"
つまり、各株が持つユニークな遺伝子の数に上限が見当たらないことを意味します。株ごとに新しい遺伝子が次々と見つかることから、研究者たちは細菌種内の遺伝的可変性が従来考えられていた以上に高い可能性を指摘しています。この研究から、GBSの遺伝子レパートリーはcore genome（全株に共通する遺伝子）とshared genome（部分的な株に存在する遺伝子）、specific genome（特定の株にのみ存在する遺伝子）で構成され、これらすべてを含めたパンゲノムで種を記述できると述べています（can be described）。次のように書かれています。
"Comparative analysis of the six newly sequenced genomes and the two genomes already available in the databases suggests that a bacterial species can be described by its “pan-genome” (pan, from the Greek word π αν, meaning whole), which includes a core genome containing genes present in all strains and a dispensable genome composed of genes absent from one or more strains and genes that are unique to each strain. "
”pan”は古代ギリシャ語の全てを意味する接頭辞で（Wiktionary）、パンデミックのパンも同じ意味だそうです。パンゲノムを直訳すれば”全てのゲノム”でしょうか。
2010年、Oksana Lukjancenkoらは公開されている61の大腸菌と赤痢菌のパンゲノムを調べ、予測されたパンゲノムは15,741の遺伝子ファミリーから構成され、コアゲノムを構成するすべてのゲノムで表現されているファミリーは993（6％）に過ぎないと報告しました（Oksana Lukjancenko et al, 2010）。2015年のMiriam Landらの総説では、2,000以上の大腸菌ゲノムを比較した結果、大腸菌のコアゲノムは約3,100の遺伝子ファミリーを持ち、合計約89,000の異なる遺伝子ファミリーとなったと報告しました（Miriam Land et al, 2015）。
最近では、Erwin Tantosoが1324個の完全長ゲノムから大腸菌パンゲノムを構築し、パンゲノムのサイズは約25,000遺伝子ファミリーで、コアゲノムはゲノムの数が増えるにつれて減少するがソフトコアゲノム（株の95％以上）はゲノムの総数にかかわらず〜3000遺伝子ファミリーで安定していると報告しています（Erwin Tantoso et al, 2022）。
なぜ株の数が増えるとパンゲノムサイズも収束することなく増加するのでしょうか？多くの株に共有される遺伝子レパートリーには限りがあるので、共有されている遺伝子ではパンゲノムサイズの永続的な増分は説明できません。パンゲノムサイズの永続的な増分は、各株が固有の遺伝子を持つと仮定したときのみ説明できます。このようなパンゲノムはオープン・パンゲノムと呼ばれ、複数の環境にコロニーを形成し環境中から遺伝子を取り込む方法を持つような種が典型的と考えられています（ヒープス則でαが１より低い、しかし１に近い時は判断が難しい可能性がある（Lars Snipen & Kristian Hovde Liland, 2015）)。
一方で、株特異的な遺伝子が少ない分類群では、パンゲノムサイズは株の数を増やしても新しい遺伝子が見つかりにくくなり、最終的には一定の値に収束（漸近）します。これはクローズド・パンゲノムと呼ばれます（（ヒープス則でαが１より高い）。クローズド・パンゲノムは、遺伝子プールへのアクセスが限られた環境に生息する種に典型的で（Duccio Medini et al, 2005）、例えば病病原菌（Xiangyu Deng et al, 2010）や、ビフィズス菌の1種（A. O’Callaghan et al, 2015）で報告があります。Staphylococcus lugdunensisは、遺伝子の水平伝播を防ぐ仕組みを持つことでクローズド・パンゲノムになっていると説明されています（Xavier Argemi et al, 2018）。
このように、細菌種は、それぞれの株特異的な遺伝子を一定の数獲得している系統と、株特異的な遺伝子数が少ないかほぼゼロの系統があるということが分かってきています。その後の研究で、core genomeはハウスキーピング遺伝子や必須の細胞機能に関わる遺伝子が多く存在し、ある種のゲノムのサブセットにのみ存在するアクセサリ遺伝子は、株特異的な病原性因子、抗生物質耐性、ニッチ適応など、生息環境特異的な機能に関係する遺伝子が多いことが報告されています（Cynthia Maria Chibani et al, 2020）。さらに遺伝子の差異比較の考え方を取り入れると、ある細菌種に特有の表現型形質は、その細菌種独自のコアゲノム、すなわち他の近縁生物と比較してこの細菌群に特有のコア遺伝子によって反映されると考えることができます（Tim Kahlke et al, 2012）。ちなみに、株特異的な遺伝子からは短い機能未知遺伝子をコードするORFやISがよく見つかるようです *1, 7。
パンゲノム解析のためのツール
パンゲノムについて簡単に説明しましたが、実際にパンゲノム解析を行うにはどうすればよいのでしょうか。パンゲノム解析でもっとも基本的なプロセスは、解析対象の2つ以上の原核生物ゲノム間における相同性の探索です。パンゲノム解析の主要な焦点は遺伝子にあるので、ゲノム間のホモログ（相同遺伝子；homologous genes *2）を探索します。相同遺伝子はオルソログとパラログに分けられ、通常はパラログと区別してオルソログを探索・報告します*3（wiki）。何らかの方法で相同遺伝子を推定し、その結果から新しい発見があるか考えることになります。
難しいのは、どのような配列探索方法と閾値設定で相同遺伝子を探索するかです。閾値の設定次第で、解析対象の分類群のコア遺伝子やアクセサリ遺伝子、株特異的遺伝子の個数は敏感に変化する可能性があります（Savio Souza Costa et al, 2020）。従って、パンゲノムソフトウェアの選択と、配列探索パラメータは注意深く設定する必要があります。
例として、Roary（Andrew J. Page et al, 2015）（紹介）というパンゲノムソフトウェアの探索手順を確認します。パンゲノム解析は遺伝子セット間の全対全比較が必要なためにNP困難の性質を持ち（Vincenzo Bonnici et al, 2018）、BLAST総当たりの配列探索の実行時間は入力データの大きさに応じてほぼ二次関数的に増大してしまいますが、Roaryは事前の配列クラスタリングをすることで計算量を大幅に軽減し、大規模なデータセットで計算不可能にならないようにしています*4。
Roaryの解析手順は論文で説明されています。Roaryの入力はProkka（Torsten Seemann, 2014）形式のGFF3フォーマットです（GFF3末尾にゲノム配列が含まれる）。GFF3ファイルから抽出されたコード領域は、全対全の比較前にCD-HIT (Limin Fu, 2012)の高い閾値設定で反復的にクラスター化されます。論文には書かれていませんが、PIRATES（Sion C Bayliss et al, 2019）と同様、各CD-HITクラスターから最長の配列が配列類似性検索の代表として使用されると理解しています。CD-HITのクラスタリングによりタンパク質配列のセットは大幅に削減されます。削減された配列に対して、ユーザー定義の配列同一性パーセンテージ（デフォルト95%）でBLASTPによる全対全比較が実行され、この全配列間の類似性の行列はMCL（ Enright A.J et al, 2002）を用いてクラスタ化されます *5。最後にCD-HITによるクラスタリング前の結果がMCLの結果と一緒にマージされます。
Roaryの主要な出力は、１）遺伝子の有無を示したgene presence/absenceテーブルとそのバイナリ形式テーブル、２）コア遺伝子数などの要約レポート、３）パンゲノムの各遺伝子の代表配列を連結したfastaファイル、４）定義されたパラメータでのコア遺伝子の代表配列だけを連結したfastaファイル（"-e" or "-e --mafft" *6）など（順不同）。さらに５）アクセサリ遺伝子の有無のバイナリファイルから作成されたNewickツリーファイル（これはアクセサリーゲノムを基に分離株を”大まかに”グループ分けするための簡単な（レポジトリでは"dirty"な）ツリー）、６）クラスタ内の配列リスト（clustered_proteins）、７）Gephiなどの大きなグラフを扱えるグラフソフトウェアで扱えるcore_accessory_graph.dotファイル；パンゲノムのコンティグレベルで遺伝子がどのように連結されているかを示したもの、も出力されます。"-e"で作成される４のmulti-FASTAアライメント”core_gene_alignment.aln”は、全てのコア遺伝子の多重整列後の連結配列となるので、必ず株間で同じ長さになっています。これを出発として、さらにトリミングなどを行い、パンゲノムの株のコア遺伝子の系統推定を行うこともできます。Roaryはこのような流れで、コア遺伝子、アクセサリ遺伝子、株特異的遺伝子は定義されます。類似のパンゲノムソフトウェアも似た出力を出してくれることが多いです。
実装はパンゲノム解析ソフトによっても異なりますが、大規模にスケールできると主張するものは、類似の計算量を低減する仕組みを持つものが多いです。最近の発表だと、PanTA(Duc Quang Le et al, 2024)（紹介）は、ゲノム数に対してリニアにスケールし、ワークステーション一台で数千のゲノム配列を扱うことができます。
パンゲノム解析ツールを使うことで、関心のある分類群における遺伝子レパートリーを調べることができます。パンゲノムは最初１つの細菌種の全てのゲノムという定義で登場しましたが、現在では、１つの細菌種だけでなく、より高次の分類階級でパンゲノム解析 *8を行った論文も出ています（Tim Kahlke et al, 2012）。また、真菌（Charley G P McCarthy and David A Fitzpatrick, 2019）、植物（Wei Li et al, 2022）や動物（Ran Li et al, 2023）、昆虫（Xiaoling Tong et al, 2022）へ適用した研究も増えています 。パンゲノム解析結果と形質情報を組み合わせて形質を説明する遺伝子を統計的に調べる「pan-GWAS解析」を行うScoary（紹介、Scoary2も紹介予定）のような実装、パンゲノム解析とmulti-omicsやバリアントプロファイルと組み合わせる包括的な解析（Ting Wang et al, 2023）、解析プラットフォーム（Xiaobo Cui et al, 2023）、Pantranscriptome（link）*9 など、さまざまなアプローチが報告されるようになっています。
コメント
今日はパンゲノムについて簡単に説明しました。ただし自分の説明が間違っている可能性もあるので、パンゲノム解析について学ぶならレビュー論文を読まれることをお勧めします。もう一回だけパンゲノム解析に関連したブログ記事を投稿する予定です。
２回目
https://kazumaxneo.hatenablog.com/entry/2024/09/25/180049
注釈
*1 株ごとに多くの特異的遺伝子が見つかる細菌種があることは、同種の別株が本当に同じ表現型であるか、そして種が分類の基本単位というのは本当なのかという疑問を生じさせる。W. Ford DoolittleとOlga Zhaxybayevaは、どのような種の定義を採用するにしても、そして今後より頑健な種の定義が提案されたとしても、その方法とカットオフ値を合理化できるかという問題が残ることから、原核生物種の存在そのものを疑い、"株"が原核生物の分類学の唯一の生物学的に意味のある基本単位であると仮定をしている（W. Ford Doolittle and Olga Zhaxybayeva）。種内のnatural ANI gapを報告した2023年の研究（Luis M. Rodriguez-R et al, 2023）や、それに続いてANI提案者のKonstantinos T氏が報告した種内natural ANI gapの単著（Konstantinos T. Konstantinidis, 2023）を読むと、この考え方はあながち間違っていないと感じる。
*2 https://ultrabem-branch3.com/informatics/bioinformatics/homolog
進化の観点から、遺伝子は相同（homologous）または類似（analogous ）に分類される。homologous遺伝子とは共通の祖先から生まれた遺伝子で、 類似遺伝子とは収斂進化によって独立に進化した遺伝子のこと。パンゲノム解析では相同遺伝子：ホモログを探し、可能ならオルソログとパラログに区別して分類する。
*3 オルソログ（オーソロガス遺伝子）とは、2つ以上の細菌に共通し、同等の生物学的機能を持つ遺伝子のことだが、パラログ（パラロガス遺伝子）は遺伝子重複によって分岐したもので、変異によって生物学的機能が変化している。オルソログはパラログよりも保存されている傾向があるとされる（オーソログ予想 (ortholog conjecture) - Nathan L. Nehrt et al, 2011; Xiaoshu Chen and Jianzhi Zhang, 2012）。自分の経験則では、これはよく当てはまるように思います。
*4 Roaryは具体的には1000個の単離株のパンゲノム探索を１つのCPUだけを使い4.5時間と13GBのRAM使用量で探索を完了するとされる（Andrew J. Page et al, 2015）。グラフィカルなレポートと豊富な下流解析サポートスクリプトも利用でき、roaryを中心としたパンゲノム解析ツールの連携する環境ができていることから人気のソフトウエアとなっており*2、PubMed調べでは2460回も引用されている（link）。ただし、Roaryは既に保守終了している点に注意（bugも修正されない）。またRoaryは同種内を対象（#514）に開発されているので、どうしても別種間で使う場合、近縁な別種でもデフォルトの感度では低すぎる可能性がある（コア遺伝子数が減り、非コア遺伝子や株特異的遺伝子が増える恐れが高い）。
*5 MCLのインフレーション値は1.5がデフォルトで、最適な結果を出すとされる。インフレーション値は通常１－２の範囲で、値を増やすとクラスタはより緻密に分割され、値を下げるとより大まかなクラスタが形成される。
*6 roaryに-eパラメータを渡すと、全てのコア遺伝子のmulti-FASTAアライメントが作成される。デフォルトではコドンを意識したアライメントを行うPRANKが使用される。遅いが正確とされる。-eに加えて"--mafft"も付けると、ヌクレオチドアライメントを行うMAFFTが使用される。非常に速いが精度は低い(Roaryマニュアルより)。
*7 短いORFは、sORFとして重要な機能を担っている可能性もある（注；sORFはおそらく断片化した遺伝子にも対応できるprodigalでもアノテーションが難しい。baktaが有効な可能性がある）。しかし短いORFは信頼性が低く、何もコードしていない可能性もある。例えばmicropanのレポジトリでは、ランダムに発生させたゲノムからprodigalを使って遺伝子予測し、その時のスコアと実際の遺伝子のスコアを比較している；ランダムに発生させたゲノム配列からのORFはほとんど短く、しかしORFとして予測されてしまっている。スコアは低いものの、実際のゲノムでの遺伝子予測のスコア分布とわずかに重複しているように見える。従って実際の解析ではランダムなゲノムから予測された時のスコアと被るORFを予めフィルタリングしてからパンゲノム解析することを提案している。これによって短いhypoのタンパク質は減るので、特に株特異的遺伝子は減る可能性がある。
*8 遠縁な系統間でパンゲノム解析を試みる場合、精度という観点からチャレンジングになる。パンゲノム解析のソフトウエアは初期から同種の株間の遺伝子を比較するために設計されていて保守的な閾値を持つのが理由。実際に同種内向けに設計されたRoaryを同属だが別種の株間で使うと、デフォルトの感度（blastp identity threshold "0.95"）では１つのオルソログを2個以上の遺伝子ファミリーに分割してしまうことが頻繁に起こる。下の図は同種（>95% ANI）のroaryの結果である（グラフに載せるため、敢えて"-i 60"で実行している）。
図1. BLASTP idenntityの分布。この１つの種の数十の株の総当たりBLASTP identity %の主要な分布は95-96%付近にピークを持つ比較的狭い分布をしている。しかしながら、カットオフが60%でもヒットするペアは見られる（パラログの可能性もある。そうやって見ると弱いピークがidentityが65、70, 80%付近にあるようにも見え、パラログの分布かもしれない）。roaryのレポートとして出力されるRplots.pdfの１スライドより。
図より、roaryデフォルトの95%カットオフでは同種の株間比較でもかなりのホモログペアを切ってしまう（デフォルト値はもっとクローン性の株に適したパラメータと推測される）。このように閾値を一義的に決めるのは種レベルでも難しく、より最近のツールでは、パラメータも異なっている。例えばPPanGGOLiNは、種レベルでの配列のアラインメントとクラスタリングに経験的に効果的なパラメータとして、--identity (default: 0.8) と--coverage (default: 0.8) がデフォルト値となっている。get-homologueではbidirectional best-hits (BDBH)が使われている。
遺伝子は正の自然選択と純化選択の影響で異なる速度で変異を獲得し（Erick Denamur and Ivan Matic, 2006）。そのため、リラックスした閾値は関連遺伝子ファミリーの過剰なクラスタリングを引き起こす危険性があり、保守的な閾値は同一遺伝子を複数のクラスターに誤分類する危険性があるとされる。Sion C Baylissらは、特に広範な分類群をパンゲノム解析の対象にする場合、単一の同一性閾値で遺伝子が同じファミリーに属するかどうか定義することは難しいと述べている（Sion C Bayliss et al, 2019）。
閾値問題は、パンゲノムではなく真核生物の種で、Lineage-specific genesやorphan遺伝子の文脈でより深く議論されている印象。2023年の遺伝子ファミリー創始者推定フレームワークGenEraの論文（Josué Barrera-Redondo et al, 2023）（紹介）では、以下のように書かれている。
”ゲノム系統分類の最大の注意点は、相同性検出の失敗、すなわち、ペアワイズローカルアライナーが中立的な配列発散によってのみしか遠縁のホモログをさかのぼることができないために、短くて進化の速い遺伝子がしばしば若い遺伝子として誤ってアノテーションされ、その結果、偽のtaxonomically restricted genes (TRGs)と見做されてしまうことである。"
この論文で発表されたGenEraは、v.1.1からDIAMONDだけでなくFoldseekを使った構造類似性探索もできるようになっている。構造類似性（つまりタンパク質の折り畳みの類似性）は一次配列相同性より進化的に保存されている傾向があり、有名なものだと、グロビン（HANS-ERIK G et al, 1994）のようなタンパク質なら、たとえkingdomを跨っていても3次構造はよく似ている（アミノ酸配列同一性だと10％程度しかない時がある）。そのため、structuralアライナーなら有意に類似した構造として検出できる可能性がある。
このように、一次配列類似性で遠縁な系統間のパンゲノム解析を行なっても、BLASTPや類似ツールの感度によって相同遺伝子かランダムマッチかの区別が難しくなる。BLASTは、クエリ配列からユーザーが定義した長さのwordを取り出し、データベース内のエントリーの中からこれらのwordに高得点でマッチするものを検索することによりホモログ検索を開始するので、ワードサイズを短くすることで感度を上げられる可能性はあるが、根本的な解決にはなっていない。
では精度の低下が懸念される場合どうするかだが、例えばTim Kahlke et al, 2012の研究では、"ソフトウェアのパラメータに基づくクラスタリングの変化を最小限に抑えるため、パラメータ値を変化させて複数回のOrthoMCLを実行し、すべての条件で安定しないクラスタを除外した"という記述がある。このように、難しい系でパンゲノム解析を行なったなら、結果が妥当かどうか検証するプロセスが必要かもしれない。
追記；あるいは、結果は精度と感度のバランスが取れたよい推定結果を提供するが、誤差も存在している可能性があると考えて、個別の遺伝子ファミリーに対してさらにセミマニュアルのキュレーションを行うステップが必要かもしれない（難しい系の時）。
(追記: 有効な追加解析として､もし特定の遺伝子にだけ興味があるなら､系統解析を推奨する｡興味ある表現型の株だけ特定のクレードを形成しないかを観察する｡つまり感度が高すぎて1つのオルソクラスターに複数のオルソログが混ざっているなら､系統推定でさらに細分化される)
*9 パントランスクリプトーム
"パントランスクリプトームは、ゲノミクスの分野で新たに生まれつつある 「パンゲノミクス 」の概念に基づいている。通常、個人のゲノムデータの変異を評価する場合、科学者は個人のゲノムをDNA塩基の一本鎖からなるリファレンスと比較する。パンゲノムを用いることで、研究者は個人のゲノムを、生物地理学的に多様な祖先を持つ個体から得られた、遺伝学的に多様な参照配列のコホートと一度に比較することができる。~ 遺伝子発現を理解するためのRNAシーケンスデータのマッピングは、RNA配列が細胞機構によってスプライシングされるため困難である。つまり、1セットのRNAデータはゲノムの非連結領域から得られる可能性があり、それらをリファレンスに正しくアライメントすることは困難である。これらのスプライシング部位は、ヒト集団全体で一様ではなく、個体によって異なる。また、RNAがどのハプロタイプに由来するのかを知ることも難しい。つまり、その遺伝子群が、その人の母親から受け継いだ染色体群に由来するのか、父親から受け継いだ染色体群に由来するのかを知ることも難しい。オープンソースツールの新しいパイプラインによって、研究者は個人のRNAのスプライスセグメントを取り出し、パンゲノム上のどこに並んでいるかをマッピングし、そのデータがどのハプロタイプに属するかを特定し、遺伝子発現を解析することができる。"
関連
RIBAP；属レベルなどでパンゲノム解析を行うために、遺伝子の近接性情報をより積極的に使用する。また、roaryの閾値を変えたときのクラスタリング結果を視覚化したインタラクティブなレポートを提供し、単一閾値でのクラスタリング結果が妥当かどうかの判断をサポートする。ランタイムがかなり長めなので、まずは少数のゲノムでテストすることを推奨。
