バイオインフォマティクス第１２回最適化シミュレーション山村雅幸（東工大総合理工）前回のまとめ自己組織化 K-means法によるクラスタリング自己組織化マップ進化系統樹の構築分子時計近隣結合法この辺まで最節約法確率モデル例：アラインメントパラメータによる異なる系統樹の構築 02/07/18 バイオインフォマティクス第12回 2 補足：Jukes & Cantorモデル(1969) 置換速度行列 A C G T A  - 3α   α  α   α  C G T α α α   - 3α α α  α - 3α α   α α - 3α  置換行列 A C G T 02/07/18 バイオインフォマティクス第12回 C st rt st st G st st rt st ( 1 1 + 3e − 4αt 4 1 st = 1 − e − 4αt 4 rt = ※無限時間後の平衡状態 A,T,C,Gが1/4ずつ A  rt   st s  t s  t ( T st   st  st   rt  ) ) 3 計算例 …A… …G… st=1/4 (1-exp(-4αt)) ｓt …C… …G… rt=1/4 (1+3exp(-4αｔ)) ｓt …C… …G… 葉ノードで同じ座位の塩基が不一致の場合 …A… ｓt …G… …G… st=1/4 (1-exp(-4αt)) …G… rt …G… rt=1/4 (1+3exp(-4αｔ)) …G… 葉ノードで同じ座位の塩基が一致する場合 02/07/18 バイオインフォマティクス第12回 4 親-葉一致 st 親-葉不一致 rt 全一致 st*st 親-片葉不一致 st*rt 全不一致 rt*rt 葉-葉不一致葉-葉一致 1vs1 1vs4 1vs16 1 0.009803 0.990197 9.61E-05 0.009707 0.980491 0.019605 0.980779 0.019228 0.018141 0.014372 0.01 2 0.019221 0.980779 0.000369 0.018851 0.961928 0.038442 0.963036 0.037021 0.033065 0.021042 0.01 4 0.036964 0.963036 0.001366 0.035598 0.927438 0.073928 0.931537 0.068867 0.055669 0.023769 0.01 8 0.068463 0.931537 0.004687 0.063776 0.867762 0.136925 0.881823 0.120744 0.082796 0.018306 0.01 16 0.118177 0.881823 0.013966 0.104211 0.777612 0.236354 0.819509 0.193694 0.106605 0.009782 0.01 32 0.180491 0.819509 0.032577 0.147914 0.671596 0.360981 0.769326 0.277712 0.126452 0.005436 0.01 64 0.230674 0.769326 0.05321 0.177463 0.591863 0.461348 0.751494 0.3467 0.14714 0.004773 0.01 128 0.248506 0.751494 0.061755 0.186751 0.564743 0.497012 0.750009 0.372763 0.157265 0.004982 0.01 256 0.249991 0.750009 0.062496 0.187496 0.562513 0.499982 0.75 0.374987 0.158197 0.005011 0.01 512 0.25 0.75 0.0625 0.1875 0.5625 0.5 0.75 0.375 0.158203 0.005011 0.01 1024 0.25 0.75 0.0625 0.1875 0.5625 0.5 0.75 0.375 0.158203 0.005011 進化速度時間 0.01 02/07/18 バイオインフォマティクス第12回 5 確率モデルの使い方確率モデル系統樹の評価＝モデルからデータを説明する尤度系統樹の作り方は別に考えるモデルは精密にできるかなりの計算量がかかる最適化問題：次回座位依存で進化速度が違ったら？ギャップがあったら？使い方予言力ｖｓ計算量、人間の手間はトレードオフ試行錯誤が必要 02/07/18 バイオインフォマティクス第12回 6 4.6 例：分子進化系統樹の構築 ① DNA配列のアラインメント処理進化距離 ② ①の結果より進化距離行列作成 A B ③ ②を用いての無根系統樹構築 C F E D 系統樹のトポロジー（④ ③に対して根の決定をする）アラインメントパラメータの依存性を知る 02/07/18 バイオインフォマティクス第12回 7 ①配列のアラインメントについて類似度の算出配列１ G A C C C T 配列２ G C C G T 配列１ G A C C C T 配列2 G － C C G T 変異、挿入・欠損の推定動的計画法アラインメント候補となるものを推定一致、不一致、ギャップ（欠損・挿入）を得点パラメータ →最大の得点和を持つものを結果として採用 02/07/18 バイオインフォマティクス第12回 8 計算機実験用いる生物真核生物（3種）・真正細菌（2種）・古細菌（2種）用いる遺伝子配列（GenomeNet - KEGG）ペプチド鎖伸張因子 EF-1α（Tu）アラインメントのパラメータをさまざまに変化近隣結合法により系統樹を構築妥当とされる系統樹ができるか 02/07/18 バイオインフォマティクス第12回 9 結果 1 パラメータ (1,-1,-4,-2） (1,-1,-3,-2) (1,-2,-4,-2）３パラメータ (1,-3,-4,-2) (2,-2,-4,-2) 系統樹のトポロジーの典型例２パラメータ (1,-2,-3,-2) (1,-3,-4,-3) (2,-2,-4,-3) ４パラメータ (1,-1,-5,-2) (2,-1,-3,-2) （一致、不一致、ギャップ、連続） 02/07/18 バイオインフォマティクス第12回 10 考察アラインメントパラメータの影響大 →系統樹の結果はアラインメントにかかっている置換ｖｓギャップ質の異なる情報を安易に混ぜてはまずいパラメータ強い依存性今回実問題 02/07/18 結果結果最適パラメータ判明事前の最適設定不可能最適パラメータ推定値バイオインフォマティクス第12回 11 系統樹構築のまとめ分子時計とUPGMA法加法性と近隣結合法置換数の最小化最尤法距離情報のみから、近隣節点の選択・結合最節約法近接配列のボトムアップなクラスタリング尤度の最大化例：アラインメントパラメータによるさまざまな系統樹の構築 02/07/18 バイオインフォマティクス第12回 12 本日のメニュー進化系統樹の構築最適化最適化とは関数最適化、組合せ最適化局所探索法：焼きなまし大域探索法：遺伝的アルゴリズム適用例：タンパク質の立体構造アラインメントシミュレーション 02/07/18 バイオインフォマティクス第12回 13 ５．最適化最適化とは局所探索法：焼きなまし大域探索法：遺伝的アルゴリズム例：立体構造アラインメント参考書自分でプログラムを組むときの Numerical Recipe in C （技術評論社) 脇本和昌「乱数の知識」(森北出版) 乱数の教科書はあまり多くない。これは手ごろ。茨木俊秀：「アルゴリズムとデータ構造」岩波講座情報科学「最適化」柳浦、茨木：「組合せ最適化」乱数・探索法を含む一通りのことを調べるのに便利最新。局所探索中心。ＧＡのあたりはちょっと不満 ※ GAについては、考え方がそれぞれ異なり、ミスリーディングなものが多いので紹介しない。 02/07/18 バイオインフォマティクス第12回 15 5.1 最適化とは optimization 最も適切な計画・設計・方策の作成・選択工学の基本的・中心的課題日常生活のありふれた活動数理的研究１９４０年代からオペレーションズリサーチ第二次世界大戦中イギリス、軍事作戦（オペレーション) システムズアプローチ 02/07/18 バイオインフォマティクス第12回 16 形式的定義解空間：Ｓ制約条件：Ｃ目的関数：ｆＸ Determine: 目的関数を最小化(最大化)する解候補形式化Ｃ Given: Ｓ最適化 min f ( x ) s.t x ∈ C ⊂ S 自分の抱えている問題は一度形式化してみると良い最適化対象による分類関数最適化、組合せ最適化 02/07/18 バイオインフォマティクス第12回 17 関数最適化重要な概念大域解（global minimum）局所解（local minima）単峰性近傍：ユークリッド空間なら十分小さい球凸集合＋凸関数＝単峰局所解は容易、大域解は困難関数最適化（数理計画）歴史的には凸計画の解法の提案から 02/07/18 バイオインフォマティクス第12回 18 凸計画の例：線形計画法例：栄養問題炭水化物は400g以上、タンパク質は300g以上、脂肪は600g以下。 2種の食品A, B 料理の価格を最小化 02/07/18 バイオインフォマティクス第12回 19 形式化 min Z = 1000 x + 1500 y s.t 400 x + 200 y ≥ 400 炭水化物 100 x + 300 y ≥ 300 タンパク質 150 x + 200 y ≤ 600 脂肪 x, y ≥ 0 min Z = 2 x + 3 y 約分 s.t 2 x + y ≥ 2 x + 3y ≥ 3 3 x + 4 y ≤ 12 x, y ≥ 0 02/07/18 バイオインフォマティクス第12回 20 ＬＰの探索空間凸空間に注意大域解は端点にあるシンプレクス法行列の掃き出しによる線形計画の解法 O(N3)程度の計算量数千変数でも実用に耐える性能 02/07/18 バイオインフォマティクス第12回 21 関数最適化のクラスと解法線形微分可能微分利用不能多峰性 02/07/18 線形計画最急降下法 Newton法共役勾配法直接探索法シンプレクス法（≠線形計画）近似探索法焼きなまし法遺伝的アルゴリズムだんだん効率が悪くなるバイオインフォマティクス第12回 22 組合せ最適化離散的＝手がかりがない最終的には全部数え上げるしかない。関数最適化の凸性に対応するものは？最適原理（Bellman 1950’s）動的計画法（dynamic programming:DP）により効率的に解決例：配列のアラインメント、確率モデルの同定 02/07/18 バイオインフォマティクス第12回 23 組合せ最適化のクラスと解法最適原理動的計画法数え上げ可能縦型・横型探索ヒューリスティクス探索分枝限定法近似探索法焼きなまし法遺伝的アルゴリズムだんだん効率が悪くなる列挙不可能（組合せ爆発) 02/07/18 バイオインフォマティクス第12回 24 5.2 局所探索法：SA local search 現在の探索点の近傍を調べる近傍適応度景観傾斜局所最適解ローカルミニマ 02/07/18 大域解バイオインフォマティクス第12回 25 焼き鈍し法（SA）エネルギー GENERATE(x) 探索点xはエネルギーE(x)を持つ現在の探索点xの近傍x’をランダム生成 ACCEPT(x) 1   −( E ( x )− E ( x ' )) T e if E ( x') ≤ E ( x ) otherwise T：温度パラメータ 02/07/18 バイオインフォマティクス第12回 26 02/07/18 バイオインフォマティクス第12回 27 SAのパラメータ初期温度ある温度での探索回数 REDUCE(Tk) 温度スケジュール、徐冷スケジュール 02/07/18 対数型(理論保証) 指数型（頻用) Tk +1 = T1 / log k Tk +1 = γTk バイオインフォマティクス第12回 (0.8 ≤ γ < 1) 28 SAの収束定理定理十分ゆっくり温度を下げてゆくならば，SAは最適解で収束する定理の意味に注意収束＝無限に時間をかけること実際にSAを用いる条件は異なる 02/07/18 全部の解候補を無限回訪問した末に最適解に収束解空間の一部しか訪問しない．だからこそ意味があるバイオインフォマティクス第12回 29 マルチスタート法局所最適解山登り，傾斜法．．→ 初期解がどこだったかマルチスタート初期解をランダム発生，局所探索を繰り返す系統的な発生方法格子状，前回となるべく離れた位置など ☆ 極めて強力 × 多次元（数十以上）になると通用しなくなる 02/07/18 バイオインフォマティクス第12回 30 局所探索のまとめ近傍山登り焼きなまし法局所解という難点はあるが、応用マルチスタート程度の工夫で相当に強力 02/07/18 バイオインフォマティクス第12回 31 5.3 大域探索法：GA simple (standard) GA 遺伝子（gene）遺伝子の列。通常固定長。固定数。例）０００１１００１００１０００１１１００個体（individual）最小構成要素。通常｛０、１｝染色体（chromosome）ミクロ／マクロな生命現象の模倣染色体の集合。通常は１染色体、１倍体集団（population）個体の集合。通常固定サイズ 02/07/18 バイオインフォマティクス第12回 32 遺伝オペレーション（ミクロ）突然変異（mutation）ある遺伝子座の遺伝子を対立遺伝子と置換例）０１１０１０１００１ Pm：突然変異率（mutation rate）１点交叉（one point crossover）２つの染色体を途中から組み換える例）０１１０１０１１００１１０００１１００１ Pc：交叉率（crossover rate） 02/07/18 バイオインフォマティクス第12回 33 遺伝オペレーション（マクロ）適応度関数（fitness function）個体の生存力・繁殖力・競争力を表す量発生・形態形成などの抽象化ルーレット選択（roulette wheele selection）用語注意 02/07/18 選択（selection）淘汰（selection）増殖（reproduction） 30.9% 14.4% 5.5% バイオインフォマティクス第12回 49.2% 34 単純GAの手順１２３４５６初期化適応度計算選択交叉突然変異停止条件 02/07/18 個体をランダムに生成各個体の適応度を計算両親をルーレット選択確率Pcで交叉して親と置換確率Pmで子を突然変異打ち切り世代で停止さもなければ２へバイオインフォマティクス第12回 35 動作例（初期世代）個体染色体表現型適応度選択率選択数１２３４０１１０１１３１１０００２４０１０００８１００１１１９１６９１４．４％１回５７６４９．２％２回６４５．５％０回３６１３０．９％１回バイナリ表現ルーレットの幅ｆ（ｘ）＝ｘ２（一例として） 02/07/18 バイオインフォマティクス第12回 36 動作例（第2世代）親染色体交叉点子染色体表現型適応度１２２４０１１０１１１０００１１０００１００１１４４２２０１１００１１００１１１０１０１００００１２２５２６１６１４４６２５６７６２５６ランダムに決定突然変異 02/07/18 バイオインフォマティクス第12回 37 遺伝的アルゴリズムの応用教科書・ウェブサイト公開ライブラリ不用意に使わないで！！単純GAベース最悪の部品の組合せ単純GA＝説明が簡単なだけ 1st choiceを紹介 02/07/18 バイオインフォマティクス第12回 38 関数最適化用交叉 α d2 d2 α d1 d1 親1 親2 α d1 α d2 実数コード化／BLX-α Eshelman 93 02/07/18 実数コード化／正規分布交叉小野95 バイオインフォマティクス第12回 39 汎用世代交代モデル first choice: MGG [佐藤 97] 02/07/18 バイオインフォマティクス第12回 40 5.4 例：立体構造アラインメント配列アラインメント配列が似ている → 構造が似ている → 機能が似ている MASAYUKIYAMAMURA YAMAMASA -mA-sayukiyAMA-murA y-Am-------AMAs---A * *** * masayukiYAMAMurA---------YAMAM--Asa ***** * 配列が似ている → × 構造が似ている → 機能が似ている立体構造アラインメント 02/07/18 バイオインフォマティクス第12回 41 既存方法の問題点研究の目的局所探索性：大域的探索： RMS-fitting DP-matching Genetic Algorithm タンパク質の形状に依存するタンパク質の形状に依存しない保存部位・活性部位が保存部位・活性部位を無視される：RMSd 重視する：評価関数不明瞭な精度評価：原子間の距離: RMSd とM(対応原子数) 精度記号 02/07/18 バイオインフォマティクス第12回 42 実数値ＧＡの適用 z s = ∑ exp( ε × dist c2i ) i =1 s + 1 .0 f = g + 1 .0 g =m−z δ = 5 .0 ε<0 評価関数：距離の近い対応の強調＝保存部位・機能部位の発見 02/07/18 バイオインフォマティクス第12回 43 比較実験：対応原子間の距離 YSAS : RMSd + DP-matching + Monte Carlo Stralign : RMSd + Graph-matching (T.Akutsu) δ=5.0, fragment length:15, max test times:100 GSA : Real-coded GA + UNDX + MGG Yale Structural Alignment Server(M.Gerstein) 集団サイズ：50, 交叉回数：100, 世代数：3000, UNDX(α=0.5, β=0.3), 突然変異率=0.01, δ=5.0, ε= － 0.8 精度記号 dist2 記号＞δ or GAP １なし２・＞1.5 and 異文字３：＞1.5 and 同文字４＊ ≦1.5 ５｜ ≦1.0 ６＃ 02/07/18 ≦0.5 記号4,5,6が多ければ良い保存部位・活性部位が表現可？バイオインフォマティクス第12回 44 2+ Ca -binding Proteins Active Sites : CD-hand , EF-hand Quotations from PDBsum 3icb.pdb 02/07/18 5cpv.pdb バイオインフォマティクス第12回 45 2+ Ca -binding Proteins NMR 02/07/18 バイオインフォマティクス第12回 46 ６．シミュレーションハイブリッドペトリネット力学シミュレーション 6.1 ハイブリッドペトリネットプレーストランジショントークンペトリネットもともとは離散事象のモデル分子＝トークン、反応＝トランジションハイブリッド←反応速度 02/07/18 バイオインフォマティクス第12回 48 6.2 力学モデリング Bioinformatics 生命現象の再現と予測ゲノム配列情報遺伝子ネットワークたんぱく質立体構造 etc. 個体器官細胞・・・・生物の個体をまるごと再現 Î 細胞レベルに注目 02/07/18 バイオインフォマティクス第12回たんぱく質分子 49 モデル生物：線虫Caenorhabditis elegans 発生細胞配置の重要性卵殻なし卵殻あり ABの娘細胞 ABp P2 ABa EMS 野生型 02/07/18 EMS P2 Lois G. Edgar et al.; Early transcription in Caenorhabditis elegans embryos, バイオインフォマティクス第12回 50 Development, 120, pp.448(1994) 改変研究の目的細胞配置シミュレーション卵殻と細胞分裂時の細胞伸長細胞レベルの力学モデリング細胞配置を再現力学的相互作用 ABp AB 野生型 P1 P2 ABa 2細胞期 02/07/18 3細胞期バイオインフォマティクス第12回 EMS 4細胞期 51 構成要素のモデリング細胞細胞膜・・・質点をメッシュ状に接続 Î ばね定数を距離の関数k(x)とする Fspr = −k ( x) ⋅ x ただし，k ( x) = k spr ⋅ x 2 半径は，質点間に働く力の安定距離を表すばね【細胞膜に働く拘束条件】質点ダンパー細胞の体積を一定にする力 FV = kV ⋅ dV dV:標準体積からの差分細胞膜を平らにする力細胞膜をランダムに振動細胞膜 Fpl = k pl ⋅θ dist θdist :隣り合う三角パッチ間の角度 02/07/18 バイオインフォマティクス第12回細胞 52 現象のモデリング細胞分裂過程細胞形状の伸長二つの中心体間に一定の大きさの斥力 FCEN − CEN = k CEN 中心体-細胞膜間にばねによる斥力中心体収縮環 FCEN − MEM = k ( x ) ⋅ x ただし，k ( x ) = k CEN − MEM ⋅ x 2 細胞収縮環の収縮収縮環の質点同士に引力 FCRNG = k ( x) ⋅ x 2 ただし，k ( x ) = k CRNG ⋅ x 中心体が離れられる最大距離細胞膜の追加･削除質点間の距離に応じて質点を追加・削除 02/07/18 バイオインフォマティクス第12回娘細胞 53 4細胞期の細胞配置の比較顕微鏡画像シミュレーション結果ダイアモンド型平行四辺形型 T型正方形型 ABp ABa P2 EMS 卵殻あり（野生型） EMS ABの娘細胞 P2 卵殻なし 02/07/18 バイオインフォマティクス第12回 54 CEN-MEM CEN-CEN max of CENs 細胞配置4パターンの分布 80 CEN-CEN 20 80 200000 80 200000 80 60 80 ABa-ABp 120 120 120 120 120 120 120 120 ABa-ABp 60 80 100000 CEN-CEN 200000 CEN-CEN 20 80 200000 80 200000 80 60 80 100 100 100 100 120 120 120 120 120 120 02/07/18 300000 100000 200000 300000 100000 200000 300000 100000 200000 300000 バイオインフォマティクス第12回 100 120 200000 300000 20 40 100 200000 80 40 100 100000 60 CEN-CEN 20 40 60 100000 300000 CEN-CEN 20 40 60 100000 300000 CEN-CEN 20 40 60 100000 300000 20 ABa-ABp 120 300000 300000 100000 200000 300000 100000 60 80 100 120 200000 300000 55 1.0 80 ABa-ABp 100 40 60 80 100 200000 120 200000 40 100 20 40 60 100 100000 100 40 100 300000 CEN-CEN 20 80 80 CEN-CEN 20 40 60 100000 300000 60 ABa-ABp 200000 80 200000 CEN-CEN 20 40 60 100000 ABa-ABp 100000 80 300000 ABa-ABp ABa-ABp 80 200000 CEN-CEN 20 40 60 100000 300000 CEN-CEN 20 40 60 200000 ABa-ABp CEN-CEN 20 100000 300000 ABa-ABp CEN-CEN 200000 ABa-ABp 300000 100000 300000 20 40 120 200000 300000 40 120 100000 120 200000 CEN-CEN 20 40 60 100000 300000 CEN-CEN 20 40 60 100000 300000 CEN-CEN 20 40 60 100000 300000 100 0.8 ABa-ABp CEN-CEN 200000 100 120 ABa-ABp 100000 300000 100 20 300000 200000 100 80 200000 100000 100 100 100000 120 100 300000 CEN-CEN 120 100 60 200000 120 60 40 100000 120 80 ABa-ABp ABa-ABp ABa-ABp 80 CEN-CEN 120 40 60 200000 100 60 0.6 100000 100 20 40 80 100 300000 CEN-CEN 20 40 100 ABa-ABp CEN-CEN 200000 60 ABa-ABp 100000 300000 80 20 40 100 ABa-ABp 200000 80 60 ABa-ABp 120 80 60 CEN-CEN 20 40 ABa-ABp 80 60 180 CEN-CEN 20 40 ABa-ABp 100 60 160 CEN-CEN 20 40 ABa-ABp 80 ABa-ABp ABa-ABp 60 CEN-CEN 20 40 40 100000 CEN-CEN 20 140 ABa-ABp CEN-CEN 20 120 ABa-ABp 100 ABa-ABp CEN-CEN 80 max of CENs 0.4 CEN-MEM 60 考察野生型と同じ細胞配置を再現するには中心体が離れられる最大距離（max of CENs）が大きく影響中心体間の力（CEN-CEN）と中心体-細胞膜間の力（CEN-MEM）に相対的な関係が必要＝細胞骨格の役割が重要 CEN-CEN CEN-MEM 細胞中心体 max of CENs 02/07/18 バイオインフォマティクス第12回 56 結論細胞レベルの力学に基づくモデル細胞，卵殻，細胞分裂 Î 野生型の細胞配置を再現可能細胞配置シミュレーションでのパラメータ解析野生型の細胞配置を再現するには卵殻と細胞伸長との力学的な相互作用細胞が十分に伸長細胞の伸長する力が大きいことが必要 02/07/18 バイオインフォマティクス第12回 57 発展野生型と異なる分裂パターンと比較検証【生物学実験】レーザー照射特定細胞にレーザーを当てて細胞分裂を止める突然変異体等分裂になる，分裂方向が異なる，など 02/07/18 バイオインフォマティクス第12回 58 レポート課題（確認）課題関心のある問題をひとつ取り上げ、講義で紹介したモデルのいずれかの上で形式化してみる。今回紹介した最適化でも良い。余裕があれば、予想される利点・欠点を述べよ。A4 １枚程度で。例未知種のゲノム上で、遺伝子領域を予測する → ウインドウで配列をスキャンするニューラルネットによりモデル化、近縁種の既知の遺伝子領域のデータを教師とする。シグナル伝達経路の微分方程式モデルの未知パラメータを推測する。解空間はｎ次元の数値パラメータ、制約条件はパラメータの範囲、評価関数は実験によって得たタイムコースデータとの誤差、これを最小化する。 02/07/18 バイオインフォマティクス第12回 59 総まとめ機械学習自己組織化クラスタリング、SOM ※ コンピュータは道具であるプロは道具を使い分ける進化系統樹構築道具は研がねば使えない ↓ 近隣結合法、確率モデル問題クラスとの適性に注意せよ最適化ニューラルネットワーク、決定木焼鈍し法、遺伝的アルゴリズムシミュレーション 02/07/18 バイオインフォマティクス第12回 60