大規模遺伝子ネットワーク推定ソフトウェア SiGN

人の体は，そのほとんどがタンパク質から出来ています．そのタンパク質の設計図が DNA という形で記録されている遺伝子です．生物がその生命活動を維持するためには，必要なタンパク質を適切なタイミングで作り出す必要があります．生命の最小単位は細胞ですが，この細胞内でタンパク質が作られます．細胞の中でタンパク質生成の制御を行っているのもまたタンパク質です．つまり細胞内にはどのような遺伝子を機能させるか，遺伝子同士が協調作用し，遺伝子の機能を制御するネットワークが出来ているのです．遺伝子の機能を制御するのは遺伝子だけではありませんが，制御の基本となるのは遺伝子同士の制御であり，この制御のネットワークを「遺伝子ネットワーク」と呼んでいます．

遺伝子のネットワークは生命現象の基本と言って良いものですので，したがって生命現象を理解するためには「遺伝子ネットワーク」の理解が欠かせません．病気の仕組みを理解し，その原因となる遺伝子を突き止めたり，副作用の少ない新しい薬を開発するためにも，遺伝子ネットワークを知る必要があります．しかし細胞内で行われている遺伝子間の制御は非常に複雑で，制御の関係が分かっている遺伝子は極々一部に限られます．人には 2 万を超える遺伝子があるとされており，それらすべての遺伝子が他のどの遺伝子に影響を与えるのかをしらみつぶしで調べるのは現実的ではありません．しかし現在主に行われている分子生物学における研究手法は，まさにこの「しらみつぶし」を行っている状態であり，一つ一つの遺伝子間の関係を解き明かすために非常に時間がかかり，また出来る範囲も非常に限られたものになります．

細胞内の全ての遺伝子の活動の様子（遺伝子発現）を一度に観測することが出来れば，どの遺伝子が活動しているときにどの遺伝子が活動している，といった情報から遺伝子間の関係を網羅的に知ることが出来ます．しかし現在の技術では，細胞内のすべての遺伝子に対応するタンパク質の量を一度に計測することは難しく実現されていません．DNA として記録されている遺伝子が実際にタンパクに「翻訳」される前に，細胞内ではメッセンジャー RNA (mRNA) という物質に変換されます．この mRNA のすべてがそのままタンパクに翻訳されるされるわけではありませんが，タンパクを計測する代わりに細胞内の遺伝子 1 つ 1 つの mRNA の量を「マイクロアレイ」と呼ばれるツールで計測することが出来ます．

1 枚のマイクロアレイによって計測された「遺伝子発現データ」は細胞内のある瞬間に発現していた遺伝子 (mRNA) の数値データです．これは言わば細胞内の遺伝子発現の「スナップショット」です．1 枚のマイクロアレイからでは意味のある情報を取り出すことはできません．しかし数多くのマイクロアレイを集めたり，あるいは刺激を細胞に与えた後で計測するなどすると，そこから意味のある情報を取り出すことが出来るようになります．たとえば，薬を細胞に投与し，1 時間毎にマイクロアレイを計測すれば，投与した薬に対してどのような遺伝子が反応したが，さらにそれらの遺伝子の活動が時間が進むにつれてどのように変化していったか，といった情報を得ることができます．また，遺伝子の発現に外乱を与えたマイクロアレイを，様々な遺伝子に対して用意することにより，細胞内の遺伝子ネットワークの全体像を浮かび上がらせることができるようになります．

マイクロアレイ等によって観測された「遺伝子発現データ」から「遺伝子ネットワーク」を数学的・統計学的な計算によって推定する技術があります．本研究プロジェクトで開発している「大規模遺伝子ネットワーク推定ソフトウェア SiGN （サイン）」もその一つです．遺伝子ネットワークの推定は観測された事象から内部の仕組みを解明・予測する方法であるため「生命のリバースエンジニアリング」と呼ばれることもあります．

SiGN は遺伝子ネットワークを数学的に表現するために「ベイジアンネットワーク」といわれるモデルを用います．ベイジアンネットワークでは，数学的には有向グラフと呼ばれるものでネットワークを表現します．簡単に言うと，遺伝子を点で表し，遺伝子間の依存関係を矢印で結んだものです．ベイジアンネットワークは遺伝子ネットワークに限らず，様々な現象の依存関係を解き明かすために広く用いられているモデルですが，データからベイジアンネットワークのネットワーク構造を予測することは非常に計算に時間のかかる問題であるため，大規模な事象の解析にはあまり使われてきませんでした．

本研究プロジェクトではこのベイジアンネットワークを使って，2 万以上あるとされるヒト全遺伝子を網羅可能な遺伝子ネットワーク推定ソフトウェアを研究・開発しています．現在開発中の次世代スーパーコンピュータによって実現される 10 ペタフロップススケールの計算により，これまで 1000 遺伝子程度が限界であったベイジアンネットワークによる遺伝子ネットワーク推定を一気に全遺伝子を網羅可能になります．大規模な遺伝子ネットワークでは推定されるネットワークの精度にやはり限界があります．従って，比較的小規模なネットワークを 10 ペタフロップススケールの計算により超高精度に推定する技術も同時に研究・開発しています．

遺伝子ネットワーク推定技術により，これまで研究者の長年の知識や経験を頼りにしらみつぶしに行ってきた遺伝子機能の同定や，遺伝子間の発現の依存関係の同定を網羅的に予測可能になります．観測されたデータに基づいて行われる「数学的に裏付けのある予測」により，生物学的な発見に繋がるまでの時間を大幅に短縮します．

たとえば，薬を投与した細胞の発現データから推定した遺伝子ネットワークを用いることにより，完全には分かっていない薬の詳しい仕組みを解明するための大幅な近道になります．個人毎の副作用のありなしや効きの違いなど，その原因・仕組みが分かっていない薬も多くあります．遺伝子ネットワーク推定技術を用いることにより，そいういった薬の仕組みの解明のための研究が高効率化されます．

病気の仕組みや既存の薬の仕組みを解明することにより，より効果の高い副作用の少ない薬の開発に繋がります．特に副作用の予測は遺伝子ネットワークなどの細胞全体を網羅するモデルなしには不可能であるため，より安全な薬の開発に役立つことが期待されています．

より研究の現場に近い部分での効果についても同様です．たとえばこれまで行った来た遺伝子ネットワークによる解析では，推定することができる遺伝子ネットワークのサイズが数千遺伝子程度に限られていたため，あらかじめ旧来の統計学的手法などにより遺伝子を選別しておく必要がありました．本研究プロジェクトで開発されるソフトウェアを用いれば事前知識なしに全遺伝子を網羅的にモデル化可能であるため，そのような選別を行う必要が無があります．これにより，本来は重要な遺伝子であるのに，旧来の手法によってうまく選別されず，遺伝子ネットワークによる解析対象そしてそれによる新しい発見の対象にならないといった問題を克服することが出来ます．

前述の通り，SiGN はベイジアンネットワークを遺伝子ネットワークモデルとして用いる遺伝子ネットワーク推定ソフトウェアです．一般的に知られるベイジアンネットワークは離散的データ，つまり On/Off などの二値や 0, 1, 2 などの整数値，あるいは +/-/none などのカテゴリカルなデータをモデリングの対象にしています．マイクロアレイによる遺伝子発現データは連続値（実数値）であるため，これを離散化すると本来持っている情報が失われてしまいます．したがって SiGN では情報を離散化せずに実数値のまま扱います．具体的には B-Spline によるノンパラメトリック回帰という方法を用いてベイジアンネットワークにおけるノード（遺伝子）の親子間のモデル化をします．ただし，この方法は離散モデルに比べると 100 倍程度計算に時間がかかります．

ベイジアンネットワークの推定で問題になるのはネットワーク構造の決定です．ベイジアンネットワークではデータに最適な非巡回有向グラフ (DAG: Directed Acyclic Graph) を探索する必要がありますが，解候補となるこの DAG の数が膨大にあります．たとえば 10 ノードで 4.17 × 10¹⁸ もの DAG の総数つまり解候補があります．20 ノードでは 2.34 × 10⁷² とさらに膨れあがり天文学的な値になります．アルゴリズム（計算方法）を工夫し効率的な探索を行うことにより 30 ノード程度まででしたら，現在のスーパーコンピュータを用いてこの膨大の解候補のなかから最適な構造を探すことが可能です．しかしそれ以上のサイズの大規模なネットワークでは最適解を探索することは事実上不可能です．従って，このような大規模な遺伝子ネットワークの構造探索にはいわゆる発見的アルゴリズムと呼ばれる局所解探索を行います．１つの発見的アルゴリズムで適用可能な遺伝子ネットワークのサイズが限られるため，本研究プロジェクトでは，推定したい遺伝子ネットワークの大きさにあわせて複数のアルゴリズムを研究開発しています．

NNSR (Neighbor Node Sampling & Repeat) アルゴリズムは本研究プロジェクトで開発された全く新しい遺伝子ネットワーク推定アルゴリズムです．これはスカラ型スーパーコンピュータ専用に設計された超並列アルゴリズムであり，従来 1,000 ノード（遺伝子）程度が限界であった遺伝子ネットワーク推定の限界を大幅に超え， 2 テラフロップス程度の計算資源で全ヒト遺伝子に匹敵する 20,000 ノードの遺伝子ネットワークを 2 時間で推定可能になりました．より計算資源を投入することで，高精度化が可能な他，遺伝子発現の個人差などを考慮した 100,000 ノード規模の遺伝子ネットワークを推定する場合，現在の 1,000 倍以上，すなわち 10 ペタフロップススケールの計算資源が必要になります．

1,000 遺伝子程度までの大規模遺伝子ネットワークの構造推定は発見的アルゴリズムの一種である Greedy Hill-Climbing (HC) アルゴリズムを用います．この規模の遺伝子ネットワークの構造推定では，発見的アルゴリズムを用いるという理由もあり，１回の構造推定では得られた遺伝子ネットワークに信頼性の問題が生じます．これを克服するためにブートストラップ法と呼ばれる繰り返し計算を行います．

ブートストラップ法は繰り返し計算により推定されるモデルの確率的な確からしさを評価できる方法です．この計算には，いわゆるデータ並列性があるため超並列化が可能でスカラ型スーパーコンピュータの性能を最大限に享受することができます． 従来のテラフロップス級のスパコンで 1,000 回程度の繰り返し計算が可能であったものが 10 ペタフロップススケールの計算により 1,000,000 回以上の繰り返し計算が可能なり，推定される遺伝子ネットワークの確率的な安定性・確度および精度を飛躍的に高めることが可能になります．

すでに上で説明したように，観測したデータに最も適したネットワーク構造を，すべての可能な解候補の中から探索することは非常に困難です．しかし当研究グループでは過去に動的計画を用いた非常に効率の良い最適解探索 (OS) アルゴリズムを考案し，20 遺伝子程度と非常に小規模な遺伝子ネットワークであれば最適ネットワーク構造の探索が可能になりました．OS アルゴリズムは大量のメモリを使用する必要もあって大規模なメモリを搭載する共有メモリ型の並列計算機では並列動作が可能ですが，次世代スパコンで採用される分散メモリ型スーパーコンピュータでは並列動作は出来ません．本研究プロジェクトではこの OS アルゴリズムをベースに，分散メモリ型スーパーコンピュータ上で動作可能なアルゴリズム Para-OS の開発に成功しました．

この新しい並列アルゴリズムを用い，2 テラフロップス程度の計算資源を 199 時間 (約 8 日間) 用いて 31 ノードの遺伝子ネットワークの最適構造探索に成功しました．この時，使用したメモリは総計で 464.3 GB でした．これまで論文などで報告されているベイジアンネットワークの最適解構造探索は 29 ノードが最大でしたので，この記録は世界記録になります．現在は次世代スーパーコンピュータ上で超並列動作をするための改良を行っています．ベイジアンネットワークの構造探索はノード数に対して指数時間で増加していくため，10 ペタフロップス，つまり今の 5,000 倍の計算資源を理想的に使用できたと仮定しても 10 ノード程度適用可能なネットワークサイズが増えるだけです．しかし局所解ではない最適解が求めることが出来ますので，10 ペタフロップススケールの計算により，データに対して最適な遺伝子ネットワークの推定をより大きなサイズで行うことが出来ることになり，最適遺伝子ネットワーク推定アルゴリズムの応用範囲が大幅に広がります．

遺伝子ネットワークは細胞内の遺伝子間の制御ネットワークをモデル化するものであるため，生物学的な知識や他の情報をうまく活用することにより探索すべき解候補を絞ることが出来ます．この考えを用いることにより，解候補となるネットワーク構造に制約を与え，その上で最適ネットワーク構造を探索するアルゴリズム COS の開発に成功しています．現在アルゴリズムのプログラム化が完了している段階で，今後次世代スーパーコンピュータで動作可能な並列版アルゴリズムの研究開発を進める予定でいます．

• ノンパラメトリック回帰によるベイジアンネットワーク推定プログラム
• マイクロアレイなどの遺伝子発現データから遺伝子間の依存関係を表す遺伝子ネットワークを予測・推定
• 次世代スーパーコンピュータによる 10 ペタフロップススケールの計算に対応
• ベイジアンネットワーク構造探索アルゴリズムとして以下のものを実装
  • Neighbor Node Sampling & Repeat (NNSR) アルゴリズム：20,000 遺伝子以上の超大規模遺伝子ネットワークが推定可能な新アルゴリズム．超並列・分散メモリ型スーパーコンピュータが前提
  • Greedy Hill-climbing アルゴリズム：〜 1,000 遺伝子程度の大規模遺伝子ネットワーク推定アルゴリズム
  • 並列版 Optimal Search (OS) アルゴリズム：30 遺伝子前後まで対応可能な動的計画法を用いた最適解探索アルゴリズム
  • COS アルゴリズム：制約付き最適解アルゴリズムで 70 遺伝子前後まで対応可能
• それぞれのアルゴリズムでブートストラップ計算が可能
• MPI による並列実行に対応．ブートストラップ計算は SGE を利用しグリッド計算にも対応．
• 2008 年 5 月から開発をスタートし，旧来の計算ルーチンの刷新と MPI 対応化による超並列動作が可能に
  • 行列演算ライブラリの使用 (LAPACK)，無駄な計算の省略：従来比 2〜3 倍の高速化
• 開発言語： C99
• ソースコード：計 47,869 行 (09/10/29 現在)
• CSML による出力に対応し，CI online や CSML Pipeline での直接表示・解析が可能