E5-2 アラインメントされた配列集合からモチーフを 福本翔平 抽出する方法 北上始 森康真 広島市立大学情報科学部知能工学科 広島市立大学大学院情報科学研究科知能工学専攻 731-3194 広島市安佐南大塚東 3 丁目 4 番 1 号 E-mail: s20160@edu.ipc.hiroshima-cu.ac.jp {kitakami, mori}@hiroshima-cu.ac.jp あらまし配列データベースから類似部分の多い部分配列, すなわち, モチーフを抽出す る方法は, 数多く提案されている. その中でも, 文字の出現頻度を考慮した GS( ギブスサ ンプリング ) 法は, 最も有名な抽出方法として知られている. しかし, その方法は必ずし も正しいとは限らない. その理由の一つとして,GS 法は, 配列データごとに, 類似部分配 列の位置を確率的に計算している為, それらの位置が常に変動するという問題がある. 本 稿では, その問題を解決する為に, 配列データベースを予めアラインメントした上で, ギ ャップを考慮した新しいプロファイル計算法を用いて高い出現頻度を持つ類似部分配列を 抽出する方法を提案する. また, 提案手法と従来手法 (GS 法 ) の比較実験を行ったので, その実験結果について報告する. キーワードデータマイニング, テキストマイニング, バイオインフォマティクス, 科学 データ管理 科学データベース 1. はじめに 配列データベースから類似部分が存在するパターンを抽出する方法は, テキストデータの規則的な共通部分を取り出すだけでなく, アミノ酸などの分子配列データからモチーフを抽出する方法として幅広く使用されている. アミノ酸は 20 種類存在し, それぞれにはアルファベット1 文字を対応させて表現している. モチーフとは, アミノ酸配列で生物学的に重要な機能を果たす特徴的な類似パターンのことである. 自然界には, さまざまなモチーフが存在するため, それらは,PROSITE [2] や Pfam などのデータベースで管理されている. 同じ機能を持つモチーフを集めてみると, それらのアミノ酸配列は, お互いに類似しており, 完全に一致しない場合が多い. すなわち, モチーフ の表現には曖昧性が含まれている. このためモチーフの曖昧性は正規表現を用いて表現されている. 正規表現により, モチーフに含まれるアミノ酸配列の類似性を簡単に把握することが出来る. 正規表現されたモチーフを導出するためには, 類似する部分配列を配列データベースから抽出する方法が大変重要になるが, 様々な抽出方法が存在する. その中でも最も有名な抽出方法として知られているのが,GS( ギブスサンプリング ) 法である. GS 法は文字の出現頻度を考慮した確率的最適化アルゴリズムであり, 抽出する配列の長さ k をユーザが指定する事で, その長さ k を持つ類似部分配列 ( 以後,k- 類似部分配列と呼ぶ ) を出力する事が出来る. しかし,GS 法は, 計算結果の精度が
毎回変動し, 必ずしも良い結果を出力するとは限らない. 理由の一つとして,GS 法は配列データごとに, ランダムに与えた初期値に基づいて類似部分配列の位置を確率的に更新している為, 計算途中でそれらの位置が常に変動し, 結果が安定しないという問題が発生する. 本稿では, この問題を解決する為に, 配列データベースを予めアラインメントした上で, 高い出現頻度を持つ類似部分配列を抽出する GS-Align 法を提案する. 具体的には, 先ず, マルチプルアラインメント [1] により各配列の最適な位置にギャップを挿入する. これにより, 各配列の長さを整え, 配列間の類似部分を同じ位置に配置したデータベースを作成できる. 次に, ギャップ文字ができるだけ類似部分配列に含まれないようにするために, ギャップを考慮した新しいプロファイル計算法を導入し, プロファイルを計算する. 最後に,GS 法で用いられている評価方法を取扱い, 類似部分がより多く存在するパターンを抽出する. これにより GS 法を使用する際に常に変動していた結果を可能な限り抑え, ある程度まとまった結果を抽出可能にする. 以下, 本稿の構成を示す.2 章では類似部分配列抽出に関する関連研究について述べる.3 章は従来の手法である GS 法の詳細について,4 章では提案手法である GS-Align 法の詳細について述べる.5 章では GS 法と GS-Align 法の結果を比較して評価を行い,6 章では本稿のまとめと今後の課題について述べる. の多い文字列に変更するという基本操作がある. この基本操作を繰り返し実行することで, モチーフの構成要素となる k- 類似部分配列を探索する. この一連の動作は統計値を推定するアルゴリズムではなく, 解を繰り返し求め直す事で近似の解を求めていく焼きなまし法と言う手法である. 正確には焼きなまし法の一種である確率的最適化アルゴリズムと見做すことができる. GS 法は k- 類似部分配列を抽出するにあたって非常に有効な手段となっている. しかしながら, 出力されたデータは, 必ずしも正確なモチーフを抽出しているとは限らない. 配列データベース内の k- 部分文字列を取り出す際は, 配列ごとに確率的に k- 類似部分配列を見つけ出すという計算を行っているので,k- 類似部分配列の存在位置の探索が不安定であり, 結果が本来抽出すべきモチーフから外れる可能性があるという問題を抱えている. この問題を解決するために, 本稿では, 以下の処理手順から成る GS-align 法を提案している. (1) 予め収集した配列データベースをマルチプルアラインメントする. その結果をマルチプルアラインメント済み配列集合と呼ぶ. (2) マルチプルアラインメント済み配列集合に対して新しいプロファイル計算法を適用し, プロファイルから相対エントロピーが最大となるクラスタを1つだけ選択する. ただし, クラスタとは, マルチプルアラインメント済み配列集合の同じ列から選択される k- 部分文字列の集合をさす. 2. 関連研究配列データベースから k- 類似部分配列を抽出する方法で, 最も有効な手法は,Lawrence らが提案した GS( ギブスサンプリング ) 法 [3][4] である. GS 法では, 配列データベース内からランダムに取り出した複数の文字列の出現頻度を計算して頻度 3. 従来の抽出手法本章では, 従来の手法である GS 法を用いた類似部分配列の抽出方法と, 抽出した類似部分配列 (k- 部分配列集合 ) の評価方法について説明する. k- 部分配列集合とは,GS 処理において, 配列データベースの各配列からランダムに取り出される長
さ k の部分文字列集合の事である.k の値はユーザー側が任意で与え, それにより k の値分の部分文字列を取り出す仕組みとなっている. 以下に配列データベースと k - 部分配列集合の詳細を示す. 3.1 ギブスサンプリング配列データベース DB は n 種類の文字からなる文字集合 = {a 1, a 2,, a n } で定義されているとする. また, 全配列数はDB = {s 1, s 2,, s N } の式から N 本として見なす. る生起確率を計算する. それにより n k 個の結果が算出され, それらを n k 行列で表現したものをプロファイル (p n,k ) と呼ぶ. プロファイルは出現頻度を求める際に必要となる. (2) 出現頻度候補となる k- 部分文字列集合に存在する1 つの文字列 x = a 1, a 2,, a k の出現頻度 P x を計算する. 計算方法は P x = p 11 p 22 p nk と定められている. ちなみに x の部 位に存在する文字 a i が行列の i 行目に対応するならば,p ij は j 列目における文字 a i の生起確率と見なす. これにより, 文字列 x の確率が高ければ k- 部分配列集合の総意に類似し, 低ければ類似しない事が意味される. 図 1 配列データベースと k- 部分配列集合 GS 法の主な目的は図 1 のように, 文字列集合 である DB = {s 1, s 2,, s N } から, ユーザが定めた k 値分の k - 部分配列を取り出し, お互いにできるだ 図 2 k- 部分配列集合とプロファイル け類似した部分配列集合となるように様々な計算 方法を行って変更していくものである. その計算を行うためには, プロファイル, 出現頻度, 背景頻度の三種類の方法が挙げられる.GS 法を実行した際に存在する k - 部分配列集合に対して, これら三種類の計算方法は以下のように定められている. (1) プロファイルプロファイルの初期値は, 各配列データからランダムに選択された k- 部分配列の集合を用いて計算される. その後は, 新しく計 (3) 背景頻度解候補である k - 部分配列集合以外の部位 BS に存在する各文字 = {a 1, a 2,, a n } の出現確率を背景頻度としている. 文字 a i の背景頻度 b ai は BS に存在する文字 a i の生起確率と見なす. これにより k- 部分文字列集合に存在する一つの文字列 x = a 1, a 2,, a k の背景頻度 Q x は,b a1 b a2 b ak と計算する事で求められる. 算された出現頻度や背景頻度を用いて, 再 計算される. 取り出した k- 部分配列集合に おいて,k 個の列ごとに n 個の各文字に対す
ばれる評価関数を用いている. その計算を行うた めには先ず, ベイズ統計解析を考慮したプロファ イル E ij を式 (1) のように定義する. E ij = (C ij + b i ) ((n 1) + B) (1) C ij とは, プロファイル p nk の i 行目に該当する 図 3 k- 部分文字列群に対する BS 法 GS 法は,DB = {s 1, s 2,, s N } からランダムに選択された行列 Z を用いる事で, 出現頻度が高くかつ背景頻度の低い k- 部分文字列集合を抽出する処理を行っており, そのアルゴリズムを図 4に示す. 文字が j 列目に現れる数である.n は配列総数,B は (n) 1 2と定め, プロファイルの i 行目に該当する文字の全配列に対する相対出現頻度をf i とする. また, プロファイルの i 行目に該当する文字の疑似度数 b i はf i Bとしており, 分子のゼロ除算を回避するた めに扱われている. 1 DB の各配列に対して,k- 部分配列の開始点 st i をランダムに選び, それらを行列順に並べた k- 部分文字配列 S = {st 1, st 2,, st N } を初期値とする. 2 DB からランダムに一つの配列 Z を選択する. 3 Z 以外である N-1 個の配列データベース DB -Z から図 2 のようなプロファイル (p n,k ) を算出する. 4 配列 Z の長さを L と見なす.Z 内に存在する l i k + 1 (i = 1,, L) 個の k- 部分配列 x について, 出現頻度 P x および背景頻度 Q x を計算し, この計算によって算出されたプロファイル E ij による相対エントロピー F は以下の式 (2) となる. また,DB 内に存在する文字の種類は, 扱う DB の 種類によって変動するので, ここでは 20 種類と仮 定する. k 20 F = C ij log i=1 j=1 ( E ij b i ) (2) この式を k- 部分文字配列に当てはめる事によ って, 得られた値が 0 に近ければ類似部分配列と して近似しており,0 よりもマイナス側に遠ざかれ ば類似していないものとして判断する事が出来る. 双方の比である R x = P x /Q x を算出する. 5 {R 1, R 2,, R li k+1} (i = 1,, L) となった各値から, 比例した確率でランダムにE r を選択し, E r に対応する k- 部分配列を新たな開始点 st Z として更新する. 6 結果が収束するまで2~6を繰り返す. 繰り返し回数は多いほど良い結果が出力されるが, その分実行時間が大幅に伸びる. 図 4 GS 法のアルゴリズム 4. 提案する類似部分配列抽出法本章では, 配列データベース DB から比較的安定した類似部分配列を抽出させる GS-Align 法を提案する. そのために, 先ず,DB に含まれる配列データの長さを統一するために利用されるマルチプルアラインメント操作について説明する. 次に, マルチプルアラインメントが行われた DB から同じ長さ k の類似部分配列の集合を抽出する GS-Align 法について述べる. GS-Align 法と GS 法 3.2 k- 部分配列集合の評価法 配列データベース DB から抽出する類似部分配 列を評価する方法として, 相対エントロピーと呼 を利用するに当たって, どちらも, ユーザ側が k- 部分文字配列の k 値を予め設定しなければならな い.
4.1 マルチプルアラインメントマルチプルアラインメントとは生物学などで扱われる手法の一つであり,DNA やアミノ酸等といった配列を類似した部分で特定できるように並べ替えたものである. 以下では, これを単にアラインメントと呼ぶ. 図 5にアラインメントの例を示す. この手法によって, 配列データの各文字が他の配列データのどの文字に対応するのかを決めることができる. このため, 類似部分配列を一目で見つけやすくなる. アラインメント結果には, ギャップと呼ばれる記号 (-) が存在しているが, これは類似部分を整列化させる為に組み込まれた記号である. アラインメントする為に用いられる表現方法であり,1 行目はシーケンスデータの詳細,2 行目以降は実際のデータの文字列で構成されている. 本稿でもその形式を採用する必要があるが, PROSITE 内のアミノ酸データは FASTA 形式で既に記述されている為, 抽出を行うのみで良く, 余計な変換をする必要はない. これによってアラインメントされたデータを文字列集合 DB として扱い, 提案手法であるプログラム GS-Align 法に与え, 類似部分を抽出していく. ただし, 問題となる部分もあり,ClustalX によるアラインメントはノイズも少なからず関与しているので, 完全にアラインメントされた結果が出 力される訳ではない. 4.2 提案手法 GS-Align 法は,GS 法の問題点である出力結果の変動や精度を改善することを意図して提案された手法である.GS-Align 法で入力するデータはアラインメントした DB 以外に, 抽出する配列の長さ 図 5 アラインメントの例 k を設定する必要がある.k の値は抽出するモチー フの長さから決定し,DB はギャップも含まれて 配列データベース DB に含まれる各配列データの長さは不統一であるが, ギャップを組み込むことで同じ長さに統一している.GS 法は DB の長さが不統一であっても動作は可能であるが, 解を安定的に取得できない. 我々が提案する手法は, アラインメントを予め実施するので, 類似部分の整列化と長さの統一化が図られ, 解の安定的な取得が期待できる. そのアラインメントを行うプログラムとして, 本稿では ClustalX [1] と呼ばれている系統解析用のプログラムを使用した. 扱うデータに関しては PROSITE [2] から抽出した特定のアミノ酸データを用いる. 動作を行うためには先ず,ClustalX に読み いる為その分の長さも考慮して決定する. GS-Align 法ではギャップを考慮した新しいプロファイル計算法により算出されるプロファイル G nk と,GS 法で行った k- 部分配列集合の評価関数である相対エントロピーの式 (1)(2) を用いる. プロファイルG nk の計算方法として, 先ず GS 法で用いたプロファイルp nk をG nk とする. その後 k 列ごとの生起確率の合計値 SUM を求め, 列ごとに存在するギャップの生起確率を計算する 最後にそれらを n 個の各文字に割り当てるように1 文字ごとの確率 R を加算する. プロファイルG nk を算出するための, 新しいプロファイル計算方法は以下の式 (3)(4) の通りである. 込ませるデータを FASTA 形式に変換する必要があ る.FASTA 形式とは, 塩基配列やアミノ酸配列を R = 1 SUM n (3)
G nk = G nk + R (4) 評価する部分は GS 法と同じく予測される k- 部分 タを抽出し, 類似部分配列として表現する. 図 7 に提案するアルゴリズムを示す. 配列集合である t k だが, 参照する範囲は異なる. 先ずアラインメントによって統一された DB 全体の長さを L とすると,1つの k- 部分配列集合の長さが k である事から,L k + 1 個の k- 部分配列集合が作られる. 以下では,DB を矩形の N L 行列とみなし, 各集合 (k L 行列 ) を特にクラスタと呼ぶ ( 図 6).L k + 1 個のクラスタを GS 法の評価法に基づいて計算し, 相対エントロピーの値が最も高いクラスタを類似部分配列集合と見なして抽出する. 1 ClustalX にてアラインメントされた DB から, 長さ k の値分の部分配列集合 (t k) を1 クラスタと見なす. 2 プロファイルp nk から1クラスタ分のプロファイルG nk を算出する. 3 プロファイルG nk を基に相対エントロピーの評価関数を算出する. 4 L-k+1 個分のクラスタが終了するまで,3~ 4を繰り返す. L 文字 (L-k+1) 番目 5 L-k+1 個分の相対エントロピーから最も値の 大きい k- 部分配列集合を呼び出し,k- 類似部分 N 個 1 番目のクラスタ i 番目のクラスタ クラスタ 配列集合として出力する. 図 7. 提案するアルゴリズム k 文字 k 文字 k 文字図 6 DB 内に存在するクラスタ 5. 評価実験 本章では, 閾値を 0.1 として, 提案手法の評価 実験を行う. 性能評価のために使用した配列デー また, あるクラスタにおける文字列集合の始点と終点にギャップが存在している場合は, その数をカウントし,{ 一クラスタにおけるカウント数 t 閾値 } を満たす場合に限り, 解の候補とするため評価関数を計算する. 理由の1つとして, モチーフを表現する正規表現において, 始点や終点にギャップは存在しないことが挙げられる. よって, そのような k- 部分配列集合を解の候補から除外するための処置として行っている. しかし, このような操作を行っても, 始点と終点にギャップが含まれている文字列は完全には除去できないので, そのようなギャップ入りの部分文字列が存在する場合はその k- 部分文字列だけを削除する操作を行っている. GS-Align はアラインメントされた DB を用いて タベースは,PROSITE 内に登録されているアミノ 酸データセットを5つ用いた. 扱う5つのデータ は, タンパク質の中でも有効な働きをする有名な アミノ酸配列データベースである. 詳細は以下に 示す. モチーフの長さに関しては, ギャップ込み で表 1 および表 2 に示している. 表 1 PROSITE のデータセット 番号 モチーフ名 登録番号 長さ 件数 1 Kringle PS00021 24 95 2 Homeobox PS00027 129 1316 3 PTS_EIIA PS00372 22 51 4 HTH_ASNC PS00519 37 43 5 HTH_DEOR PS00894 35 81 複数のクラスタから評価関数が最も大きいクラス
表 3 提案手法と GS 法との精度結果比較 番号 モチーフ名 提案手法 (%) 従来手法 (%) 1 Kringle 79.70 64.44 2 Homeobox 79.71 87.75 3 PTS_EIIA 54.90 49.18 4 HTH_ASNC 59.27 41.76 5 HTH_DEOR 20.07 17.06 3 [DENQ]-x(6)-[LIVMF]-[GA]-x(7)-[LIVM]-A 表 2 各モチーフの正規表現番号正規表現 ( ギャップ有 ) 1 [FY]-C-[RH]-[NS]-x(18)-[WY]-C 2 [LIVMFYG]-[ASLVR]-x(2)-[LIVMSTACN]- x-[livm]-{y}-x(32)-{l}-[liv]-[rknqest AIY]-[LIVFSTNKH]-W-[FYVC]-x-[NDQTA H]-x(80)-[RKNAIMW] -[LIVM]-P-H-[GAC] 4 [GSTAP]-x(2)-[DNEQA]-[LIVM]-[GSA]-x(2 )-[LIVMFYT]-[GAN]-[LIVMST]-[ST]-x(6)- R-[LIVT]-x(2)-[LIVM]-x(13)-G 5 R-{G}-x(2)-[LIVM]-x(3)-[LIVM]-x(17)-[ST A]-x(2)-T-[LIVMA]-[RH]-[KRNAQ]-D-[LIV MF] 比較結果を見る限りでは, 提案方式が従来手法よりも確実に優れている訳ではない. そのため, 提案手法に対する改良の余地が未だに残されている. 改良する部分としては, ギャップの少ない DB に対する計算方法の変更にある.4 章で述べたように, アラインメントされた DB には少なからずノイズが発生する. そのノイズの割合や挿入され るギャップの頻度によって, 提案手法で抽出され 従来手法 (GS 法 ) の実行において, 処理の繰 る類似部分配列集合の結果が変わることがある. り返し回数は, データセットごとの DB 内の文字数 とした. 繰り返し回数をそれ以上に増加させたが, 相対エントロピーが既に収束しているので, 計算 結果が変化することがなかった. 表 4 データセットにおけるモチーフ内のノイズ 番号モチーフ名ノイズの割合 1 Kringle 0.17894 2 Homeobox 0.19680 従来手法と提案手法との性能を比較するために, 以下で定義される精度の式 (5) を利用する. 精度 (%)= B 100 (5) B + C ただし, そして検索で合致した範囲を B として, 検索されたノイズの部分を C とする. この式は本来抽出すべき範囲と抽出した範囲がどの程度合致しているかを数値化したものであり, 百分率 (%) で表され, 数値が高い程一致している部分が多いと見なされる. 3 PTS_EIIA 0.45098 4 HTH_ASNC 0.04651 5 HTH_DEOR 0.02500 表 3 および表 4 より, ノイズが少なくギャップ の割合も少ない HTH_DEOR は, モチーフを抽出し にくい傾向となっている. すなわち, 元の DB をア ラインメントする際に, モチーフ ( 正解 ) ではな い他の部分も類似した文字集合として整列されて しまう可能性があるので, そういった部分を類似 部分配列として探索してしまう事が問題となって 表 3 は, 従来手法と提案手法の精度を比較した いる. 表である.
1 16 31 46 61 76 91 106 121 136 151 相対エントロピー値 0-1000 -2000-3000 -4000-5000 -6000 文字列の先頭からの位置 図 8 提案手法を HTH_DEOR で実行した結果 配列パターンを含むクラスタは, 正解になるかどうかの検討が重要である. (2) 従来手法の図 3 にて行われている背景頻度の取り扱いや DB 内の文字 = {a 1, a 2 a n } の関係性を数値化した Blosum62 など, クラスタ内のプロファイル以外から類似部分を探索する新たな操作が必要である. 完全な解を探すことは難しいが, 今後, 新たな手 法を導入することにより, 高精度な近似解を見つ 図 8 において, 縦線が引かれている x 軸のパターンが抽出すべきモチーフ ( 正解 ) であり, 点で示されている部分は, 提案手法で抽出されたクラスタである. ただし, クラスタの位置とは, クラスタの 1 列目が DB 上に存在する位置をさす. この図からもわかるように, モチーフパターン ( 正解 ) 以外のクラスタにも相対エントロピーの最大値に近いピークがいくつか存在する. すなわち, モチーフ ( 正解 ) 以外に,k- 類似部分配列集合と見なされる可能性のあるクラスタが複数存在するという新たな問題がある. 6. まとめ本稿では,DB をアラインメントする事によって k- 類似部分配列であるモチーフの抽出を行った. 5 件のデータセットを用いた評価実験では,1 件のデータセット (Homeobox) 除いて, 提案手法は従来手法よりも精度が向上し, 安定した結果を出力した.Homeobox のデータセットについては, 精度が向上しなかった原因を調査し, 改良を加える必要がある. 今後の課題は, 以下のとおりである. (1) 相対エントロピーに関して, 複数のピークから正しいピークを選択する方法として, 最小汎化集合の要素を支持数でランキングする方法も考えられる. 支持数が最大となる汎化 けることができるのではないかと考えられる. 謝辞本研究の一部は, 日本学術振興会 科学研究費補助金 ( 基盤研究 (C), 課題番号 :20500137) の支援により行われた. 参考文献 [1] M.A.Larkin, G.Blackshields, N.P.Brown, R.Chenna, P.A.McGettigan, H.McWilliam, F.Valentin, I.M.Wallace, A.Wilm, R.Lopez, J.D.Thompson, T.J.Gibson and D.G.Higgins : Clustal W and Clustal X version 2.0, Bioinformatics, Applications Note, Vol.23 No.21, pp.2947-2948, 2007. (ClustalX:http://www.clustal.org/) [2] PROSITE:http://prosite.expasy.org/ [3] Lawrence C. E., ALtschul, S. F., Bogushi, M. S., Liu, J. S., Neuwald, A. N. and Wotton, J.: Detecting subtle sequence signals: A Gibbs Sampling Strategy for Multiple Alignment, Science, 263, pp.208-214, 1993. [4] Liu,J.S., Neuwald,A.N. and Lawrence,C.E.: Bayesian Model for Multiple Local Sequence Alignment and Gibbs Sampling Strategies, JASA, 90, pp.1156-1170, 1995.