多次元データモデリングソフトウェアの開発

荒川 正幹, 船津 公人


Return

1 はじめに

近年化学の分野においては、コンビナトリアルケミストリやハイスループットスクリーニング(HTS)の発達にともない、多くの情報が日々蓄積されている。これらの情報の中からいかにして意味のある知識を抽出できるかが、効率的な創薬の鍵を握っているといえる。また、構造活性相関(QSAR)、構造物性相関(QSPR)をはじめとするケモメトリックスに関する研究も盛んに行われており、数値データの有効活用の重要性が高まっている。
通常これらの数値データは多次元/多変量であり、既存の統計ソフトウェアでは取り扱うことが困難である。そこで本研究では、このような場面において多次元データの操作、解析、グラフ描画などを行うことを目指し、多次元データモデリングソフトウェアの開発を行った。本システムでは、通常の表計算ソフトウェアで取り扱われるような行列データに加え、多次元配列データを扱うことが可能である。また、様々な形式でグラフを描画することによって、このようなデータの構造を視覚的に捉えることが可能となっている。さらに、これらのデータに対し、Multi-way PLS法による回帰分析、Kohonenニューラルネットワークの学習などを行い、より詳細な解析を行うことも可能である。
これらのケモメトリックス手法を化学データへと応用した研究として、文献[1, 2]を挙げることができる。
文献[1]は、CoMFA法による3D-QSAR解析において、適切な配座および重ね合わせルールを推定するために、Multi-way PLS法を用いたものである。はじめにCoMFAフィールド行列を、異なる配座および重ね合わせルールを用いて複数計算する。次にこれらの行列を結合し、説明変数となる配列を作成する。1つのCoMFAフィールドが行列、つまり2次元の配列であるため、説明変数配列はこれに配座、重ね合わせルールの次元を加えた4次元の配列となる。そして、この4次元配列を説明変数、生物学的活性値を目的変数としてMulti-way PLS法による回帰分析を行う。その結果得られる回帰係数を解析することにより、活性を説明するのに最も適した配座および重ね合わせルールを推定することができる。最後に選択された配座と重ね合わせルールを用いて通常のCoMFAを行い、3D-QSARモデルを作成する。文献[1]では、この解析手順をcholecystokinin-B (CCK-B)の拮抗剤であるbenzodiazepine誘導体へと適用し、良好な結果を得ている。
文献[2]は、KohonenニューラルネットワークおよびMulti-way PLS法を応用することにより、分子表面上の静電、疎水ポテンシャルを用いた3D-QSAR解析を行った例である。分子表面における静電、疎水ポテンシャルはその化合物の生物学的な活性を説明するために極めて重要であると考えられる。しかし、これらのポテンシャルと生物学的活性値との間でQSARモデルを構築しようとする場合、各分子構造によって分子表面の形状が異なることが問題となる。この手法では、Kohonenニューラルネットワークを利用することによって、この問題点を克服している。まずサンプル化合物のなかから1つの分子構造をテンプレート分子として選択し、その分子表面の形状をKohonenニューラルネットワークにより2次元マップへと写像する。つまり、分子表面上に一定間隔で設定されたサンプリング点の3次元座標を入力ベクトルとし、Kohonenニューラルネットワークによる非線形写像を行う。この操作により分子表面はそのトポロジカルな関係を維持しつつ、2次元のマップへと写像されることとなる。そしてこの発火情報をもとにしてKohonenマップと同じサイズの行列に、分子表面上の静電、疎水ポテンシャルをそれぞれ写像する。こうして作成されたマップを1つの配列に結合したものを説明変数として用いることにより、3D-QSAR解析が可能となる。説明変数配列は4次元の配列となるため、モデリングにはMulti-way PLS法を使用している。そして、回帰モデルの回帰係数値を分子表面上に逆射影することによって、活性にとって重要であると思われる領域を特定することが可能となる。この手法を36のエストロゲン拮抗剤を用いた3D-QSAR解析に適用し、良好な結果を得ている[2]。
上記2つの研究においては、いずれも多次元配列を取り扱う必要があり、Kohonenニューラルネットワーク、Multi-way PLS法を用いることでこれらを可能にしている。1つ目の配座、重ね合わせルール推定手法では、Multi-way PLS法を用いることによって、配座と重ね合わせルールを同時に推定することが可能である。これらは互いに密接に関係しているため、同時に取り扱うことが重要であり、そのためにはMulti-way PLS法などの多次元配列を扱う手法が不可欠である。2つ目の分子表面上の静電、疎水ポテンシャルを用いた3D-QSAR解析においても、Multi-way PLS法を用いることで、静電ポテンシャルと疎水ポテンシャルの両方を考慮したモデルを作成することに成功している。また、分子表面の形状の違いを取り除くために、Kohonenニューラルネットワークによる非線形写像を行っている。これらのケモメトリックス手法は非常に有益なものであるが、まだ一般的な統計ソフトウェアには実装されていないのが現状である。そこで本研究では、このような解析をGUIによって簡便に行うことができるシステムの開発を行った。

2 手法

本システムでは、多次元配列データの操作や各種グラフ描画による多次元情報の可視化の他、PLS法、Multi-way PLS法、Kohonenニューラルネットワークなどのケモメトリックス手法を用いたデータ解析が可能である。以下では本システムに実装されている解析手法および本システムで用いる多次元配列のためのファイル形式について説明する。

2. 1 PLS法、Multi-way PLS法

PLS法はケモメトリックスの分野などで広く用いられている線形重回帰分析手法である。説明変数X、目的変数Yから、それぞれTUと呼ばれるスコアを抽出しモデリングを行う。その基本式を以下に示す。

ここでPQはローディング、EFは残差である。PLS法は、通常の最小2乗法による線形重回帰と比較し、頑健で予測的なモデルが得られるという特徴を持っている。また、サンプル数よりも記述子数の方が多い状況であってもモデルが構築できるという利点もある。線形重回帰分析手法であるため、上記の式を変形することにより、以下の回帰式を得ることが可能である。

ここでBは回帰係数、Eはモデル残差である。PLS法においては目的変数の数は複数であってもよいため、YおよびBはベクトルではなく行列となる。目的変数が1つ、つまりYが1次元の列ベクトルである場合はBも列ベクトルとなる。
通常のPLS法による回帰分析においては、説明変数Xおよび目的変数Yは2次元の配列、つまり行列で表現される。しかし、扱うべきデータはときとして3次元以上の多次元配列となる。このような場面において安定した回帰モデルを得るための手法として提案されたのが、Multi-way PLS法[3]である。これはPLS法におけるXYを3次元以上の多次元配列へと拡張した手法であり、データが本質的に多次元の性質を持つ場合には、通常のPLS法に比べ頑健なモデルが得られることが知られている。以下にMulti-way PLS法の基本式を示す。

Xは説明変数、TX から抽出されたスコア、Wはウェイト、E が誤差である。YUCF についても同様である。XY などのアンダーラインは、これらが行列ではなく、多次元の配列であることを示している。Multi-way PLS法は通常のPLS法を一般化したものとみなすことができ、XおよびYが行列である場合には、Multi-way PLS法と通常のPLS法は全く同じ結果を与える。回帰式についても通常のPLS法と同様に、以下の式によって表現される。

回帰係数B も同様に多次元配列となる。例えば説明変数が4次元の配列、目的変数が1次元の列ベクトルである場合、回帰係数は3次元配列となる。

2. 2 Kohonenニューラルネットワーク

Kohonenニューラルネットワーク[4]はKohonen によって1982年に提案された、教師なし学習を行うニューラルネットワークの一種である。多次元空間の情報を、トポロジーを保持したまま2次元のマップ上へと非線形写像する能力を持つことから、多変量データの次元縮約やクラスタリングを行うためのツールとして応用されている。Figure 1にKohonenニューラルネットワークの概要図を示す。ネットワークは入力層と出力層からなっており、出力層の2次元平面上にはニューロンが並べられている。まずサンプルが入力層を通してネットワークへと入力され、勝者ニューロンの探索が行われる。出力層上の各ニューロンの重みベクトルと、入力ベクトルとの間で次式によって距離が計算され、最も近いニューロンが勝者ニューロンとして選択される。


Figure 1. Structure of KNN.


ここでxiは入力ベクトルのi番目の値、wijはニューロンjの重みベクトルのi番目の値、mは記述子の数、djはニューロンjと入力ベクトルとの間の距離を示す。そして選択された勝者ニューロンとその近傍のニューロンの重みベクトルを、入力ベクトルへと近づく方向へと修正する。

ここで、j*は勝者ニューロン、tは学習回数、aは学習係数を示す。すべての入力ベクトルについてこの操作を行うことで1世代の学習が完了する。そして、aの値を減少させながら十分収束するまでこれを繰り返すことにより、自己組織化学習が行われる。その結果、入力空間において近い距離にあるサンプルは出力層上においても近い位置へ、遠くにあるサンプルは遠い位置へと写像される。つまり、トポロジカル情報を保持したまま、入力空間のサンプルが2次元マップへと非線形写像される。

2. 3 多次元配列のためのファイル形式

統計ソフトウェアなどで一般に用いられているcsvなどのファイル形式では、多次元配列の情報を取り扱うことが困難である。そこで、本システム独自のファイル形式(mdd形式)を定義した。mddファイルは複数の多次元配列を格納することが可能なテキストファイルである。ファイルの先頭に、格納されている配列の個数が記述され、それに続いて各配列のデータが記述される。各配列の情報は、配列の名前、次元数、各次元の要素数、各次元の名前、各次元の各要素の名前、数値データの順に記述される。mddファイル形式を用いることにより、本システムの入出力となる複数の多次元配列を、1つのファイルへと格納することが可能となる。

3 実行例

以下に本システムの実行例を示す。1つ目はKohonenニューラルネットワークの学習を行った例、2つ目はMulti-way PLS法を用いた配座、重ね合わせルール選択手法[1]の計算例である。
Figure 2に本システムのメインウィンドウを示す。図左下のリストボックスは、配列一覧部であり、現在システムに読み込まれている配列の名前が一覧表示される。ここでは"配列_1"と、それを用いてKohonenニューラルネットワークの学習を行った結果である"発火位置"、"重みベクトル_X"の3つの配列が読み込まれている。この配列一覧部において選択されている配列の次元情報が左上のコンボボックス群によって示され、そこで任意の2つの次元を指定することによって配列内の値がウィンドウ右側の行列表示部に表示される。


Figure 2. Main window of the software.

ここでは、配列"重みベクトル_X"を例として説明を行う。この配列は10*10*10の3次元配列であり、各次元の名前は、マップ_縦、マップ_横、記述子となっている。それぞれの次元の各要素には以下のラベルがつけられている。
マップ_縦 : 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
マップ_横 : 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
記述子 : Fe, Ti, Ba, Ca, K, Mn, Rb, Sr, Y, Zr
ウィンドウ左上の次元情報を格納するためのコンボボックスは5つ用意されており、5次元までの配列を取り扱うことが可能である。ここでは、配列が3次元であるため、上の3つのみが使われており、残りの2つについてはNoneと表示されている。各コンボボックスには先頭に次元名が格納され、それに続いて各要素のラベルが格納される。例えば1つ目のコンボボックスの内容は以下のようになる。
マップ_縦, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
これらのコンボボックスを用いて多次元配列の任意のスライスを指定することによって、その値をウィンドウ右側の行列表示部に表示する。固定する次元についてはその要素ラベルを、その他の次元については次元の名前をコンボボックスで指定する。この例では、3つ目の次元をFeに固定しており、右側の行列では、行がマップ_縦に、列がマップ_横に対応している。このような仕組みによって、多次元配列の任意のスライスを指定し、行列状のデータとして画面上で確認することが可能となる。また表示されたスライスをcsvファイルとして書き出すことも可能である。
また、これらのデータを用いてグラフを描画することも可能である。Figure 3に本システムによって描画したグラフの例を示す。本システムでは、棒グラフ、折れ線グラフ、2次元/3次元散布図、曲面グラフの描画が可能である。各グラフはマウスにより回転、拡大、縮小、移動などの操作が可能である。


Figure 3. Graph examples.

Figure 4にKohonenニューラルネットワークの設定ダイアログを示す。このダイアログにおいて、計算条件を設定し、計算開始ボタンを押すことによって、Kohonenニューラルネットワークの学習を行うことができる。設定可能な計算条件は、マップサイズ、学習回数、学習範囲、学習率の値と、トーラスマッピング、オートスケーリング、データシャッフルのオン/オフである。


Figure 4. Configuration dialog of Kohonen neural network.

マップサイズは出力層の一辺のニューロン数、学習回数は学習を繰り返す数である。学習範囲は重みベクトルの修正を行う際、勝者ニューロンの周囲のどれだけの範囲のニューロンを修正するかを指定する値である。例えば学習範囲が5であれば、各学習において第5近傍までのニューロンが入力ベクトルへ近づく方向へと修正される。学習率は式(6)におけるaの初期値であり、学習回数とともに減少していく。設定可能な値の範囲は0から1である。大きな値を指定すれば少ない学習回数で収束するが、大きすぎると適切な写像が行われない可能性がある。トーラスマッピングの設定をオンにすると、出力層としてトーラス状の平面が用いられる。トーラスマッピングでは、マップの上辺と下辺、右辺と左辺がそれぞれ隣接しているものとして、近傍範囲の計算が行われる。これにより、多次元空間から2次元平面への写像が無理なく行われることが多いため、Kohonenニューラルネットワークでは一般的にトーラスマッピングが用いられる。オートスケーリングの設定をオンにすると、入力サンプルの各記述子が平均0、分散1となるようにスケーリングされる。これにより、各記述子の平均値や分散の違いに影響を受けないモデリングが可能となる。記述子間で単位が異なる場合など、その平均値や分散に特に意味がない場合には、通常前処理としてオートスケーリングを行う。データシャッフルの設定をオンにすると、入力サンプルの学習の順番をランダムにすることができる。Kohonenニューラルネットワークの学習では、先に入力されたサンプルよりも、後に入力されたサンプルの影響が大きくなることが知られているが、この設定を用いることにより、学習の順番による影響を少なくすることが可能である。
学習条件の設定を行った後、計算開始ボタンを押すことによりKohonenニューラルネットワークの学習が行われ、計算結果が自動的にシステムに読込まれる。計算結果は、"発火位置"、"重みベクトル"という2つの配列として得られる。"発火位置"は2次元の配列であり、それぞれの入力サンプルがマップ上のどの位置に発火するのかを示したものである。各サンプルについて、それが発火するニューロンのX座標とY座標およびそのときの重みベクトルと入力ベクトルとの誤差の値が格納されている。"重みベクトル"は3次元の配列であり、各記述子に対する重みの値が格納されている。1つのマップが2次元であるため、これに記述子の次元が加わり3次元の配列となる。得られた重みベクトルを用いてFigure 5に示すような曲面グラフを描くことによって、マップを視覚化することができる。重みの値が面の高さによって表現されている。この図を見ることによってKohonenマップの様子を視覚的に捉えることが可能である。また、この図はマウスにより、移動、拡大、縮小などの操作が可能である。さらに、色やフォント、座標軸などについての様々な設定を行うことも可能である。


Figure 5. Example of surface graph.

2つ目の例として、Multi-way PLS法を用いた配座、重ね合わせルール選択手法[1]の計算例を示す。CCK-Bの拮抗剤であるbenzodiazepine誘導体16化合物について、14種類の配座および3種類の重ね合わせルールを用いてCoMFAフィールド値の計算を行い、説明変数Xを作成した。説明変数は4次元の配列であり、サイズは16 (化合物数) * 5096 (記述子数) * 14 (配座数) * 3 (重ね合わせルール数) である。この計算は、本研究室で別途開発されたソフトウェアによって行い、mddファイルを介して、システムへと入力された。次に生物学的活性値を目的変数として、Multi-way PLS法による回帰分析を行った。 Figure 6にMulti-way PLS法の設定ダイアログを示す。


Figure 6. Configuration dialog of Multi-way PLS.

説明変数、目的変数のコンボボックスには、現在システムに読み込まれている配列の一覧が表示される。そのなかから説明変数、目的変数として用いる配列を選択し、計算を行う最大成分数、センタリング、スケーリング、クロスバリデーションの有無を指定した後、計算開始ボタンを押すことによりMulti-way PLS法による回帰分析が行われる。ここでは、センタリング、クロスバリデーションを行い、3成分までの計算を行った。計算結果として、"R2乗"、"回帰係数"、"スコア_T"、"Q2乗"、"予測値"の配列が得られ、メインウィンドウの配列一覧部へ追加される。
次に回帰係数の2乗和の比率を計算し、活性配座および適切な重ね合わせルールの推定を行った。説明変数が4次元の配列であるため、得られる回帰係数は3次元の配列であり、そのサイズは5096 (記述子数) * 14 (配座数) * 3 (重ね合わせルール数) となる。この配列を用いて文献[1]の定義に従い、下記の式によってWCWAの2つの指標を求めた。

この値が最大となるものが、活性配座および適切な重ね合わせルールであると考えられる。Figure 7WCの計算例を示す。Multi-way PLS法によって得られた回帰係数の配列を選択した状態で、メニューから[2乗和の比率を計算]を選択することにより(a)のダイアログが表示される。そして、基準とする次元を指定し計算を行うことにより、(b)に示すような結果が得られる。


Figure 7. Calculation of WC.

4 まとめ

多次元/多変量の数値データを効果的に処理することを目標に、多次元データモデリングソフトウェアの開発を行った。数値データをもとにしてグラフを描くことにより、的確にデータの様子を把握することが可能である。また、Multi-way PLS法、Kohonenニューラルネットワークなどのモデリング手法を利用することによって、詳細な解析を行うことができる。本システムは今後Web上で一般に公開することを予定している。

参考文献

[ 1] K. Hasegawa, M. Arakawa, K. Funatsu, Comput. Biol. Chem., 27, 211-216 (2003).
[ 2] K. Hasegawa, M. Arakawa, K. Funatsu, Comput. Biol. Chem., 27, 381-386 (2003).
[ 3] R. Bro, J. Chemom., 10, 47-61 (1996).
[ 4] T. Kohonen, Self-Organizing Maps, 3rd Ed., Springer (2000).


Return