系統核心實做與架構

本研究主要的目標是探測基因序列中的調控模組，亦即搜尋出一組調控模組 G.RM 使得 Score(G.RM) 為最大值。我們利用模擬退火法【Kirkpatrick et. al.

1983】為核心並且參考 Gibbs Sampler 【Lawrence et. al. 1993】的流程，撰構出 我 們 的 主 程 式 核 心 程 式 ，SAMLA （ Simulated Annealing for Multiple Local Sequences Alignment）。（圖3.5.1）為虛擬碼。原始的模擬退火法中，每一次迭 代只需要隨機的選取新的狀態，並且以波茲曼概率分布來決定是否改變目前狀 態，進而逼近最佳解。系統若以此方式來進行預測，勢必需要花費較長的時間來 收尋，方能收斂於最佳解。因此我們參考 Gibbs Sampler 的流程【Lawrence et. al.

1993】【Liu,J.S. 1994】【Liu et. al. 1995】，改進原先隨機選取新狀態的過程，亦 即，我們將會針對所有可能的調控序列位置進行比較，選取分數高者與目前的狀 態做比較，來決定下一迭代的新狀態。

程式之初，我們從 G.C（所有 PRM 的集合）中隨意選出起始的模組

RMcurrent。每一次的迴圈中，我們針對每條輸入之基因序列 gi 執行三種不同的操

作來嘗試改變 RMcurrent 的狀態，分別為：加入新的 PRM（ADD-Operator），移 除一組已存在的 PRM（DELETE-Operator），將一組已存在的 PRM 與其他代 換（SWAP-Operator）。除此之外，由於模擬退火法的原始設計為尋找能量函數 E(*) 的最小值，但是我們現在的問題卻是必須尋找 Score(*) 的最大值，因此我 們修改了狀態變動時的機率計算方式：原先是計算 exp(-ΔE/t) 值，我們改成計 算 exp(ΔE/t) 來符合尋求最大值的要求。

（圖3.5.1）SAMLA的虛擬碼

z ADD-Operator：走訪基因序列 gi 中所有的 PRM，c∈gi.C^[-G.RM]，嘗試 加入新的 PRM 於現有的模組。計算 Score(G.RM_current∪c)。最後回傳 分數最高的新調控模組 G.RM^+c。時間複雜度為O(l_i^m)，m 表示調控序 列的各數，li 則是目前正在操作的序列 gi 的長度。

For all c∈gi.C^[-G.RM], DO {

Calculate Score(G.RMcurrent∪c);

}

Return the best scored G.RM^+c; ADD-Operator

(1) Score(*) ← energy function (definition in 3.2);

(2) T ← annealing temperature schedule;

(3) G.RM_current ← randomly initialize from G.C;

(4) While( !( Score(G.RM_current) no changed || System “freezes” || iteration > 200) ) {

(a) For each g_i∈G, DO {

i. G.RM_next ← max{ ADD-Operator(g_i,G.RM_current), DELETE-Operator(g_i,G.RM_current),

SWAP-Operator(g_i,G.RM_current) };

ii. ΔE = Score(G.RM_next) - Score(G.RM_current);

iii. Draw a Uniform(0,1) random variable U;

iv. If( U ≦ min{1, exp(ΔE/t)} ), then G.RM_current ←G.RM_next;

}

(b) G.RM_current ← SHIFT-Operator(G.RM_current);

(c) t ← T_n; }

(5) Return G.RMcurrent as G.RM;

z DELETE-Operator：此操作會移除目前 G.RM_current 中屬於基因序列 gi

的一個 PRM，c∈gi.RMcurrent。計算 Score(G.RMcurrent/c)。最後回傳 Score(G.RM_current/c) 分數最高的新調控模組 G.RM^-c。時間複雜度為 O(max_sites_per_seq)，max_sites_per_seq 為使用者輸入的參數，代表每 個輸入的基因序列中最多允許模組的個數。

z SWAP-Operator：將目前 g_i.RMcurrent 中的所有 PRM 與 gi.C^[-G.RM] 中的 PRM 作置換動作。此操作會針對所有的 cRM∈gi.RMcurrent 以及 cSWAP

∈G.C^[-G.RM] 來嘗試置換的結果，時間複雜度為O(l_i^m)，m 表示調控序

列的各數，l_i 則是目前正在操作的序列 gi 的長度。此操作回傳置換之 後分數最高的新調控模組 G.RM ^SWAP。

SWAP-Operator

For all cRM∈g_i.RMcurrent, DO {

S ←G.RM_current / c_RM ;

For all c_SWAP∈g_i.C^[-G.RM], DO {

Calculate Score(S∪cSWAP);

} }

Return the best scored G.RM^SWAP; DELETE-Operator

For all c∈gi.RMcurrent, DO {

Calculate Score(G.RMcurrent / c);

}

Return the best scored G.RM^-c;

SHIFT-Operator 來避開此種困境【Lawrence et. al. 1993】。

z SHIFT-Operator：為了避免程式落入全域最佳位移解，我們在每一次迭 代之後，將目前 G.RMcurrent 中每一個調控序列 pi 向左以及向右移動一 個核甘酸位置。最後回傳分數最高的移動方式。SHIFT-Operator 的時 間複雜度為O(m)。m 表示調控序列的個數。

SHIFT-Operator

RM_shift ← NULL;

For all p_k∈G.RM_current, DO {

(1) A_left ← A_current;

(2) If( ai,j == k ), then set ai,j = 0 and ai,j-1 = 1, for 1≦i≦n, 1≦j

≦li, and ai,j∈A_left;

(3) Form one new RM_left by A_left and D;

(4) A_right ← A_current;

(5) If( ai,j == k ), then set ai,j = 0 and ai,j+1 = 1, for 1≦i≦n, 1≦j

≦li, and ai,j∈A_left;

(6) Form one new RM_right by A_right and D;

(7) If( Score(RM_right) > Score(RM_left) ), then RM_shift ←RM_right; Otherwise, RM_shift ←RM_left;

}

If( Score(G.RM_current) > Score(RM_shift) ), then return G.RM_current; Otherwise, return RMshift;

（圖3.5.2）黑色粗線為輸入的序列，G = { gi | 1≦i≦n }。假設真正的調控模組 G.RM 由 p1，p2 所 組成。而又有另一調控模組 G.RM’ 由 x1，p2 所組成。程式在迭代搜尋 G.RM 時，會因為 x1 與

p1 重疊而無法找出比 Score(G.RM’) 更高分數的 G.RM’’ 而停止程式。因此在每次迭代搜尋時 必須測試位移調控序列 pi，來解決落入區域最佳位移解的困境。

g₁

g2

gn

…

: 調控序列 p1

: 調控序列 p2

: 假調控序列 x1