開発ソフトウェアの記載フォーマット

• ヒトの全遺伝子・転写産物を対象したネットワーク解析を可能にする大規模遺伝子ネットワーク推定ソフトウェア。
• PPIチャレンジ：１０００×１０００の超大規模計算を可能にする網羅的タンパク質間相互作用推定ソフトウェア（タンパク質ドッキング解析プログラム）。
• 1人あたり５０万SNPを用い、疾患と遺伝子の関連を発見し（GWAS）、薬物反応性を予測するプログラム。
• 「個」のデータを「一般」のモデルに合理的にフィットさせる生命システムのためのデータ同化を柱としたプログラム群。
• 以上を、融合し、統合的に利用活用するためのソフトウェア環境。

人間は一人一人の人がそれぞれ自分だけの個性的な遺伝子を持っていて 人によって病気のかかかりやすさやかかりにくさが違います。このことを逆に利用して、ある病気、たとえば糖尿病や乳がんなどの 病気にかかった人を何千人・何万人と集め、病気にならなかった人をまた何千人・何万人とあつめて、その遺伝子を比較すると、 病気にかかった人たちが共通に持つ遺伝子が見えてきます。このような遺伝子を疾患関連遺伝子と言います。 このような遺伝子は、病気を引き起こす原因である可能性が高くそれを調べることで、病気の原因の解明が進み、さらには 病気の治療法や治療薬の開発へとつなげることができます。このような研究を関連解析と言います。 関連解析により、糖尿病、乳がん、リューマチ、肥満、高血圧など次々と疾患関連遺伝子が見つかってきています。

全ゲノム関連解析（Genome-wide association study; GWAS）では、患者集団とコントロール集団の間で極めて多数（数十万～数百万） のSNPの頻度比較のための検定を行います。

従来は、この数十万～数百万の検定を独立と考えて、Bonferroniの補正という方法で、有意水準を決めていました。 しかし、問題は、この数十万の検定が互いに独立ではなく、連鎖不平衡という関係があることです。 そのため、Bonferroniの補正では、検定が保守的になりすぎて、本当に疾患に関連している遺伝子でも、 見落としてしまうと言うことが起こりえます。このように、疾患に関連している遺伝子を見落としてしまう エラーをtype I errorといいます。 我々は、連鎖不平衡がある条件下で、type I errorの確率を求め Haplotype頻度から多項分布を利用し、計算する理論式およびアルゴリズムを論文を2008年に発表しました。 連鎖不平衡は連続した複数座位におけるハプロタイプ頻度から計算できます。 このことにより、従来の方法でも見落としていた疾患関連遺伝子が見つけることができるようになり、 新たな治療法開発へとつながることが期待されます。

Misawa K, Fujii S, Yamazaki T, Takahashi A, Takasaki J, Yanagisawa M, Ohnishi Y, Nakamura Y, Kamatani N. New correction algorithms for multiple comparisons in case-control multilocus association studies based on haplotypes and diplotype configurations. J Hum Genet. 53:789-801, 2008.

プログラムパッケージparaHaploは、この理論式とアルゴリズムを組みこんだソフトウエアです。 全ゲノム関連解析は、今後データが揃えば、数万人×数百万SNPもの大量データに関して解析を行うので、時間がかかります。 特に、連鎖不平衡がある場合のtype I errorの確率の計算は、従来法に比べてさらに時間がかかります。 しかし、患者さん達は 一刻も早い治療を必要としています。 そこでparaHaploは、並列コンピューティング技術を使い、PCクラスタ上で高速に計算できるように開発されました。 将来のスーパーコンピュータへの移植も視野に入れたアルゴリズムを採用しています。

Misawa K, Kamatani N. 2009. ParaHaplo: A program package for haplotype-based whole-genome association study using parallel computing. Source Code Biol Med 4:7. 。

入力されるデータ形式としては、HapMap Database、胞状危胎を用いたハプロタイプデータベース、数千人のゲノムワイドSNP遺伝子型データより得られるハプロタイプデータに対応しています。

ハプロタイプを基礎とした関連解析ソフトウェアでは、ハプロタイプ推定に時間がかかることがわかったため、ハプロタイプ推定も並列化して高速に行うプログラムも実装しました。これをとりいれたプログラムパッケージを、paraHaplo version 2.0として公開しました。 Misawa K, Kamatani N. 2009. ParaHaplo 2.0: a program package for haplotype-estimation and haplotype-based whole-genome association study using parallel computing. Source Code Biol Med 5:5.

プログラムパッケージがparaHaploとしてsourceForgeから公開されています。 http://sourceforge.jp/projects/parallelgwas/releases/?package_id=9706

全ゲノム関連解析では単純な検定を行う場合もあるが、一つのSNP座位に2つ以上の検定を行い、低い方のP値を採用することも多い。 また、ハプロタイプを基礎にした検定でもどのハプロタイプに注目するかにより複雑な検定となる。そのように複雑な検定を用いた場合、検定統計量の分布は単純な分布を取らず、正確な検定ができないことになる。 Permutationは症例とコントロールの表現型の個人の表現型を入れ替えることにより検定統計量の分布を経験的に求める方法である。 SNP数と個体数が多いとPermutationの数は膨大となり、ペタフロップスレベルの計算機が必要となります。

そのプログラムには、我々のアルゴリズムの他に、通常のPermutation、さらには Kimmel and Shamir (2006)のRAT法を使ったハプロタイプ解析のアルゴリズムも組み込みました。 このプログラムはparaHaploとしてsourceForgeから公開されています。 http://sourceforge.jp/projects/parallelgwas/releases/?package_id=9706

入力されるデータ形式としては、HapMap Database、胞状危胎を用いたハプロタイプデータベース、数千人のゲノムワイドSNP遺伝子型データより得られるハプロタイプデータに対応しています。

ベイジアンネットワークによる大規模遺伝子ネットワーク推定ソフトウェア

SiGN（サイン）は DNA チップ・マイクロアレイによる遺伝子発現データなどから遺伝子間の発現の依存関係を表す遺伝子ネットワークを推定するソフトウェアです．遺伝子ネットワーク推定技術により，これまで研究者の長年の知識や経験を頼りにしらみつぶしに行ってきた遺伝子機能の同定や，遺伝子間の発現の依存関係の同定を網羅的に予測可能になります．観測されたデータに基づいて行われる「数学的に裏付けのある予測」により，生物学的な発見や新規薬剤標的遺伝子の同定，薬剤作用機序同定に繋がるまでの時間を大幅に短縮します． 10 ペタフロップススケールの計算により従来では不可能だった遺伝子ネットワーク推定が可能になり生物・医薬分野での適用・応用範囲が大幅に広がります．

より詳しい解説は大規模遺伝子ネットワーク推定ソフトウェア SiGN のページをご覧ください．

以下の特徴を持っています．

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

L1正則化法による生体内分子の大規模ネットワーク推定プログラム (L1-regularization for Gene Network inference)

• 数万から数十万トランスクリプトームに関する遺伝子発現データから、生体内分子間相互作用を表す大規模遺伝子ネットワークを推定するプログラム
• グラフィカルガウシアンモデル、ベクトル自己回帰モデル、同時方程式モデルに対するパラメータ推定、構造学習アルゴリズムとして以下を実装
  • Weighted Lasso アルゴリズム: 2009/12/06現在、数千ノードのスーパーコンピュータを使用して、数千から数万遺伝子を対象とした遺伝子ネットワークの推定が可能。
  • Recursive Elastic Net アルゴリズム: 2009/12/06現在、数千ノードのスーパーコンピュータを使用して、数千から数万遺伝子を対象とした遺伝子ネットワークの推定が可能。Weighted Lassoアルゴリズムよりも計算量がかかるが、より高精度なネットワークを推定することができる。
  • Relevance-Weighted Recursive Elastic Net アルゴリズム: 2009/12/06現在、複数の実験条件下における数百遺伝子を対象とした遺伝子ネットワークを推定し、実験条件間のシステム的な違いを抽出することが可能。次世代スーパーコンピュータの使用を前提として、数千から数万遺伝子ネットワークを対象とした遺伝子ネットワークの推定・比較に拡張可能。
• 情報量規準によるモデルの自動選択
• Sun Grid Engine により 統計解析システム R のバッチ処理を並列化
• 遺伝子発現データ解析プラットフォームSBiP (旧名称: CSML Pipeline) のコンポーネントとして使用可能

図: L1GNによって推定された4実験条件の遺伝子ネットワークを、Cell Illustrator Onlineで表示した図。

状態空間モデル(State Space Model)による時系列データからの遺伝子ネットワーク推定

• マイクロアレイなどの遺伝子発現時系列データから生体システムに関する以下の情報を推定するプログラム
  • 遺伝子ネットワーク：遺伝子間の依存（制御）関係
  • 転写モジュールネットワーク：共発現する遺伝子群の依存（制御）関係
  • 発現パターン予測モデル：生体システムの抽象的シミュレーションモデル
• 特徴
  • 次元縮約による超高次元短時系列データからの構造パラメータの推定
    • 遺伝子発現時系列データは、超高次元短時系列データ。計測される遺伝子数(数千から数万）に比べて時点数が極端に少ない（通常15点弱）。通常このようなデータから、上記の情報を推定することは困難。本手法ではデータに含まれる代表的変動パターンを生成する低次元の動的システムを推定することで困難を克服
  • 制約付きEMアルゴリズムによる、一意性を保証したパラメータ推定
  • 情報量規準による最適モジュール数の決定
  • 繰り返し計測(replicate)データの利用による推定精度の向上
  • 動的予測によるシステム間の差異の探索

生化学反応シミュレーターのデータ同化プログラム

生命現象をシステムとして理解する上で，細胞内外を張る巨大な生化学反応ネットワークの動的特性を明らかにすることが必要不可欠です．LiSDASは，計算機内で仮想生化学反応系を大量に作り出し，パラメータとネットワーク構造の大規模ショットガンサーチを行うことで，高品質のイン・シリコモデルを効率的に設計するための生命科学アプリケーションです．従来モデラーが手作業で行ってきたシミュレーターの生化学反応パラメータや変数の初期条件の決定，あるいは生化学反応系のネットワーク構造決定に要する膨大な作業を，計算機上の統計的学習機能で代替し，効率化することが目的です．

開発プログラムの機能

• 生化学反応系イン・シリコモデルのパラメータ探索
• ネットワーク構造の最適化及びロバスト設計
• ネットワーク構造の予測機能

生化学反応ネットワークとシステム生物学 細胞内で起こるタンパク質の発現およびその生化学機能の決定には，タンパク質や核酸など，生体内分子同士の相互作用メカニズムが大きな役割を果たしてます．生化学反応系において最も基本的な役割を担う分子は，DNA（deoxyribonucleic acid）とRNA（ribonucleic acid），並びにタンパク質です．DNAのコード領域に記されたATGCから成る塩基配列の情報は，転写と呼ばれる過程を経てmRNA （messenger ribonucleic acid）に変換されます．次に，mRNAの塩基配列に則してアミノ酸が重合され，ポリペプチド鎖が合成されます．ポリペプチド（タンパク質）の合成過程は 翻訳と呼ばれます．細胞の分化や増殖，抗ウイルス免疫応答など，生命現象の多くはこのような段階的に起こる生化学反応の連鎖を介して制御されています．タ ンパク質やRNA分子の相互作用の有無を模式的に表すと，細胞内外を張る巨大な因果ダイアグラム，すなわちネットワーク（パスウェイ）が形成されます．シ ステム生物学のマイルストーンは，ネットワークの全体象の把握とその動的特性を明らかにすることです．

生化学反応系シミュレーション，そしてデータ同化へ システム生物学におけるシミュレーションの役割は，生化学反応系に内在する動的特性を計算機実験によって再現することです．地球物理や生命科学，ナノサイ エンスなど，巨大かつ複雑なシステムを理解する上で，シミュレーションに基づく計算機実験は，現代科学の方法論として揺るぎない地位を確立しました．シ ミュレーション科学の本質は，(1) 現象のシミュレート，(2) 実験データやアプリオリな知識との整合性検証 (3) モデルの再構築，というルーティンから成る帰納推論です．データ同化の概念は，シミュレーションモデルと「現実」，すなわち，実験データや科学的知識を陽 に結び付けるための「統計科学」に由来します．本来，データ同化という言葉は，気象予測などに代表される地球物理のシミュレーションから派生したもので す．しかしながら，その概念はシミュレーション科学全般に適用されうるものです．生命科学においてデータ同化の概念を体現するために開発されたツールボッ クスがLiSDASです．従来モデラーが手作業で行ってきたシミュレーターの生化学反応パラメータや変数の初期条件の決定，あるいは生化学反応系のネットワーク構造決定に要する膨大な作業を，計算機上の統計的学習機能で代替し，効率化することが目的です．

バイオ実験技術の革新 現代の分子生物学において，パラメータ推定や構造学習に比較的低コストで利用できる情報源は，mRNAやタンパク質の発現量です．DNAマイクロアレイや 質量分析計を利用してmRNAやタンパク質の発現量を網羅的に測定する技術は，現在では世界中の実験系ラボで標準的に利用できるようになりつつあります．実験技術の誕生以来，産学の垣根を超えた品質改善努力が実を結び，今や実用上十分な精度で発現量の網羅的測定値が入手できるようになり，分子生物学や医学の研究 形態は大きく変容を遂げました．とりわけ近年では，発現量の時間変化を観測した時系列データから，生化学反応経路の動的特性を理解しようという試みがなされています．観測データから効率的に情報を抽出しながら高性能のモデルを機能的に構築していく過程において，背後の知見を得るという推論方法は，まさに統 計科学といえるでしょう．

大規模モンテカルロ計算，データ同化，次世代スーパーコンピュータ ここで，データ同化の技術的な部分について概説します．われわれが提唱するデータ同化技術は，シミュレーションモデルの不完全性を前提にします．モデルは あくまで科学的仮説の表現です．したがって，誤りや未知の物理メカニズムの取りこぼし，モデルのパラメータや初期条件の不確実性など，様々な要因が積み重 なった結果，現実とは合いません．データ同化では，不完全性に結びつくあらゆるシナリオを，「モンテカルロサンプル」（粒子）で表現した上で，大量のシ ミュレーションを同時に走らせます．その過程において，各シミュレーションと観測データとの整合性を適切に評価しながら，現実を良く説明できるシナリオを 選択的に選び出します．あるいは，合わないものを淘汰します．これは計算機上で組織的に実現するための統計計算です．統計科学では，この種の計算は「逐次モンテカルロ法」（Sequential Monte Carlo Method）と呼ばれています．

逐次モンテカルロ法は，統計数理研究所の北川源四郎教授（現所長）によって，1990年代中盤に提唱され，現在に至るまで，統計的方法論として大いに発展を遂げ，多岐に渡る応用分野に波及してきました．しかしながら，逐次モンテカルロ法をデータ同化に適用する際，重大な問題が残されています．データ同化では，モデルの不完全性に結びつく「あらゆるシナリオ」を粒子で表現すると先述しましたが，それを実現するためには数億，数兆，あるいはそれ以上の粒子を計算機上で生成する必要があります．とりわけ，生命科学が対象とする系（生化学反応経路など）では，不確定性が支配する部分がかなり大きく，必要粒子数の増加は顕著です．現在の計算機性能では，数億粒子の生成は大きな壁です．

統数研データ同化グループも，同様の問題意識の下で，研究を推進しています．２００８年９月の時点で，一億粒子の実装に成功し，その成果を生物情報学の国際会議PSB２００９(Pacific Symposium on Biocomputing)において，研究速報として論文を発表しました．一億粒子の逐次モンテカルロ法の実装は，当時世界初の試みであり，現時点でも未踏の領域です．実際，遺伝子の転写制御ネットワークのインシリコモデルに一億粒子の逐次モンテカルロを適用し，その威力を実証しました．ここで示した一億粒子のアルゴリズムを実行するには，単コアによる計算時間が実時間で８日程度かかりますが（Opteron 2200, およそ5ギガフロップス），ペタコンではこの処理が数分程度で完了することになります．

ペタコン時代の到来は，シミュレーションモデル自体をも大きく変容させる可能性があります．現在われわれが用いている生化学反応系モデルは，細胞内外で起こっている真の分子メカニズムをかなり大胆に抽象化したものです．タンパク質の細胞局在など，空間的な情報も実際の生化学反応系では本質的な役割を担うことが知られています．DNAのクロマチン構造の修飾による発現のオン・オフ変換や膜の内側と外側と物質の流出入に関するマルチスケール/マルチフィジクス などを取り込んだ，これまでと比較にならないような高機能・統合的なモデルがシステム生物学において支配的になっていくことも予想されます．

大規模タンパク質間相互作用ネットワーク予測のためのプログラム

本ソフトウェアは、タンパク質の立体構造データをもとに、タンパク質間相互作用ネットワーク予測を行う。 システム生物学の対象となるような大規模な生命情報解析に、これまで大規模な解析にはあまり利用できなかった タンパク質立体構造データを活用することを可能にする。

タンパク質は生命活動の基本となる重要な要素であり、これらのタンパク質が互いに相互作用するネットワークが、シグナル伝達系など生命の重要な機能を担っている。どのタンパク質とどのタンパク質が実際の細胞内で相互作用しているのか（タンパク質間相互作用ネットワーク) をつきとめることは、生命科学における重要な課題の一つである。

一方、タンパク質の立体構造は、その機能や相互作用機構を調べるための物理化学的基礎であり、本質的な生命情報である。 2009年12月現在、公開データベースPDBには6万余の構造データが蓄積されているが、 これまでのタンパク質間相互作用ネットワーク予測にはこれらのデータが十分に活用されてこなかった。

本ソフトウェアは、形状相補性に基づくタンパク質のドッキングシミュレーションとそのプロファイルの解析により、 タンパク質間相互作用ネットワークを予測するものである。 例えば、ある疾患に関連する数千の遺伝子が発現するタンパク質の立体構造からその相互作用ネットワークを予測し、 未知相互作用の可能性を提示することにより、疾患理解への貢献が望まれる。 実際に、本ソフトウェアのプロトタイプを利用して、肺がんにかかわるシグナル伝達系に関連するとされる 遺伝子群を対象とした解析を始めている。

目標とする解析規模

• 1000 の候補タンパク質(1000 x 1000 = 1,000,000 規模の組み合わせ）の入力から相互作用ネットワークを予測する。
• 一つの対象タンパク質に対して、相互作用可能性のあるタンパク質を全PDBから予測する。

ドッキング計算を高速に行うための工夫

• タンパク質構造データに予めFFTをかけてデータベース化し再利用する。同じ内容の計算の繰り返しを避けられるため、対象のタンパク質数が多くなっても高速に計算が可能。
• ドッキングのスコア関数について、物理化学的相互作用を考慮した計算手法を導入しても、計算時間があまり増加しないような新たなモデル（rPSC法）を考案。

解析アルゴリズムの開発

• 予測構造のクラスタリングによるドッキングプロファイル解析手法を開発し、実装する。
• 相互作用予測の信頼性を評価するための手法を開発し、予測結果をランキング形式で表示する。

遺伝子相互作用は異なった遺伝子の表現型の効果の非独立のことを言う。単純なSNPごとの表現型関連遺伝子検索の作業の後は、遺伝子相互作用の検索が重要な研究分野となると考えられている。例えば、転写調節因子の遺伝子とターゲットの遺伝子間や、物理的に結合して効果を発する蛋白質をコードする遺伝子間にはこのような相互作用が存在する可能性がある。しかし、全ゲノムを対照とした場合、可能のある相互作用の数は膨大となり、ペタフロップスレベルの計算機が必要である。

特定の疾患や薬物反応性に関連する遺伝的要因とそれ以外の要因（男女、年齢や生活習慣などの要因）が多数見つかった場合、それらを用いて表現型の予測を行うことが重要となる。それにより、治療方針や生活習慣の改善の方針が決まるからである。出来る限り正しい表現型の予測のためには極めて膨大な計算が必要である。具体的には線形回帰、ロジスティック回帰、ベイズ法などを用いるが、対象とする個体数が膨大なためペタフロップスレベルの計算機が必要である。

データ解析融合プラットフォーム

既存ソフトウェア Cell Illustrator Online (CIO) をベースとして、機能モジュール追加により、本プロジェクトで開発されたソフトウェアを実装した外部計算サーバ（京など）及び BIOBASE を初めとする既存の商用・非商用データベースを統合的に利用可能にするソフトウェア環境です。この解析フロー設計ツールを用いることで、ローカルでのジョブ実行とサーバ側での実行をシームレスに実行できる環境を構築することを目指しています。また、さまざまなカスタム解析フローをGUI上で設計できるため、研究開発において必要となる解析フローをオンタイムで提供することができます。

このアプリケーション上で現在、300 以上の解析コンポーネント群が利用できるようになっています。また、本他プロジェクトで研究開発が進んでいる

• ベイジアンネットワークによる大規模遺伝子ネットワーク推定ソフトウェア (JSiGN)、
• 開発コード：L1GN L1正則化法による生体内分子の大規模ネットワーク推定プログラム (L1-regularization for Gene Network inference)
• 開発コード：SSM 状態空間モデル(State Space Model)による時系列データからの遺伝子ネットワーク推定
• 開発コード：PF データ同化ソフトウェア

は、このカスタム解析フロー設計ソフトウェアのコンポーネントとして実装されており簡単に利用できるようになっています。今後は、 300 以上のコンポーネントのバグ修正と改良、また、解析に必要となるコンポーネントの追加、ドキュメントの整備、本他プロジェクトで開発されるアプリケーションを追加できるようにする予定です。