ヒトゲノムからの有意な散在反復配列の列挙
Enumeration of Interspersed Repeats from the Human Genome
中村篤祥
1∗1
北海道大学
1
Hokkaido University
Abstract: A half of human genome is composed of repeats and most of them are interspersed repeats. Found repeats already have been registered to databases but there are some repeats that have not been registered yet and different interpretations of repeats also exist. We are now studying a frequent pattern mining method that enumerates all the statistically significant interspersed repeats. We report our trial to systematically find all the statistically significant repeats from the whole human genome, whose size is 3000 Mega-basepairs.
1
はじめに
生物の DNA の塩基配列には多くの散在反復配列が 存在している．これは，レトロトランスポゾンと呼ば れる部分の転移によっておこるものと考えられ，それ らは遺伝子の発現制御に関係し，進化における機能発 現に密接に関連しているとされている [2]． 散在反復配列の発見法として主なものに，クラスタ リングを利用する手法と頻出パターンマイニングを利 用する手法がある．クラスタリングを利用する手法と して，シーケンスに対し頻度の高い k-mer(長さ k の文 字列) の切り出し，不一致やギャップを考慮しつなぎ合 わせたものを反復配列候補とし，それらをクラスタリ ングすることで類似する反復配列候補をまとめて，そ のクラスタ中心を反復配列とする手法が提案されてい る [3]．この手法では，k 長の切り出しをつなぎ合わせ る際に経験的な知見を利用した手法を使用することや， 発見される反復配列がクラスタリング性能に依存する こと，反復配列と判断された基準がはっきりしないな どの問題がある． 頻出パターンマイニングを利用する手法として，類 似性をハミング距離やアラインメントスコアを利用し て測ることで，反復配列パターンを抽出する手法が提 案されている [4, 6]．しかし，反復配列パターンを指定 された回数以上の出現で判断しているため，統計的に 意味があると考えられる長い反復配列パターンが少な い出現回数のために検出されないといった問題がある． 我々のグループでは，最近最小サポートをランダム シーケンスにおけるパターンの出現回数分布に依存し て決める方法を提案し，パターンの出現としてアライ ∗連絡先：北海道大学大学院情報科学研究科 〒 060-0814 北海道札幌市北区北 14 条西 9 丁目 E-mail: [email protected] メントスコアにおいて局所最適なもののみカウントす る頻出パターンマイニング手法 [4] に組み込む方法を提 案した [7]．実験ではまだ，ヒトゲノムの２１番染色体 という長さが 4500 万程度の規模の文字列にしか適用し ておらず，全ヒトゲノムから散在反復配列を抽出する には計算時間がネックになっていた． 本稿では，全ヒトゲノムの解析を可能にするための 高速化の取り組みについて報告する．2
問題設定
Σ を有限集合とし，Σ の要素を文字と呼ぶ．有限個 の文字の列を文字列という．s を構成する文字数を文 字列 s の長さと呼び|s| で表す．s の i 番目の文字列を s[i] で表し，長さ n の文字列 s は s[1]s[2]· · · s[n] とも 表現される．s の i 番目から j 番目までの文字からなる 部分列を s[i..j] と表す。ただし，s[i..i− 1] は空文字列 を表すものとする．Σ に含まれない文字 ‘-’ をギャップ と呼ぶ．ギャップ付文字列とは，Σ∪ {-} の要素からな る有限長の列である．非負整数 k 及び長さの差が高々 k 以下の２つの文字列 s1，s2に対し，以下の３条件を 満たすギャップ付文字列の組 (s′1, s′2) を文字列 s1と s2 の k-ギャップアライメントと呼ぶ． GA1k s′1，s′2はそれぞれ文字列 s1，s2にギャップを 合計高々k 個挿入することにより生成される． GA2 長さが等しい．(|s′1| = |s′2|) GA3 共にギャップの位置はない．(すべての i∈ {1, ... ,|s′1|} に対して s′1[i] = s′2[i] =‘-’ ではない．) 以下の３条件を満たす (Σ∪ {-}) × (Σ ∪ {-}) 上の実数 値関数 w(x, y) をスコア関数と呼ぶ． 人工知能学会研究会資料 SIG-FPAI-B506-09SF1 任意の x, y∈ Σ に対して，w(x, y) = w(y, x) SF2 任意の x∈ Σ に対して，
max{w(x, -), w(-, x)} < 0 < w(x, x)
SF3 x ̸= y を満たす任意の x, y ∈ Σ に対して，
w(x, -)+w(-, y)≤ w(x, y) ≤ min{w(x, x), w(y, y)}
文字列 s1と s2の k-ギャップアライメント (s′1, s′2) に対 するスコア SALk(s′1, s′2) を SALk(s′1, s′2) = |s′ 1| ∑ i=1 w(s′1[i], s′2[i]) と定義する．文字列 s1と s2の k-ギャップアライメン トの集合を ALk(s1, s2) と表記する．文字列 s1と s2の k-ギャップアライメント (s∗1, s∗2)∈ ALk(s1, s2) が以下 の条件を満たすとき，k-ギャップアライメント (s∗ 1, s∗2) は最適であるという． SALk(s∗1, s∗2) = max (s′1,s′2)∈ALk(s1,s2) SALk(s′1, s′2). max (s′1,s′2)∈ALk(s1,s2) SALk(s′1, s′2) を SOPTk(s1, s2) で表す． 文字列 s における文字列 p の局所最適 k-ギャップ出 現を以下のように定義する． 定義 1 文字列 s[1..n] の部分文字列 s[h..j] が以下の条 件を満たすとき， s[h..j] は文字列 p[1..m] の局所最適 k-ギャップ出現であるという． CS1k SOPTk(p[1..m], s[h..j]) > 0, CS2k 任意の 1≤ g ≤ m+1, 1 ≤ h′≤ j +1 に対して， SOPTk(p[1..m], s[h..j])≥ SOPTk(p[g..m], s[h′..j]) CS3k 任意の 0≤ m′ ≤ m, h − 1 ≤ j′ ≤ n に対して， SOPTk(p[1..m], s[h..j])≥ SOPTk(p[1..m′], s[h..j′]) 本質的に同じ部分を異なる出現とみなさないように 以下で定義される極小のものだけ考える． 定義 2 文字列 p[1..m] の局所最適 k-ギャップ出現であ る文字列 s[1..n] の部分文字列 s[h..j] が，以下の極小生 の条件をみたすとき s[h..j] は p[1..m] の極小局所最適 k-ギャップ出現という． CS2Mk (g, h′) = (1, h) を除く任意の 1 ≤ g ≤ m + 1, h≤ h′ ≤ j + 1 に対して， SOPTk(p[1..m], s[h..j]) > SOPTk(p[g..m], s[h′..j]) CS3Mk (m′, j′) = (m, j) を除く任意の 0 ≤ m′ ≤ m, h− 1 ≤ j′≤ j に対して， SOPTk(p[1..m], s[h..j]) > SOPTk(p[1..m′], s[h..j′]) 与えられた文字列 s, スコア関数 w に対し，s にお ける極小局所最適 k-ギャップ出現の数が多い文字列パ ターンを列挙する問題を考える．Σ の文字からなる全 ての文字列を対象とするのは計算時間がかかるので、s の部分文字列のみを文字列パターンとする．この場合， すべての文字列パターンは、極小局所最適 k-ギャップ 出現が最低でも１回は存在することになる． 頻出パターンマイニングでは，頻出であるかどうか の判断は固定された最小サポートで行うことが一般的 であるが，パターンの出現し易さがパターン毎に異な るため，ある文字列パターン p に対して、その極小局所 最適 k-ギャップ出現の数が多いかどうかの判断は，は その発生しやすさに依存して決められるべきである． そこで，ある小さな正数 0 < δ < 1 に対し、ある 情報源 S0で文字列 s が発生したという仮定の下に，p の極小局所最適 k-ギャップ出現の数を Xpとした場合， Xp≥ 1 という条件の下に Xp≥ σ を満たす確率が σ 以 下となるような最小な自然数 σ を σα(p) とし，σα(p) を p の最小サポートとして使う問題を考える． σα(p) = min{σ ∈ N | PS0{Xp≥ σ | Xp≥ 1} ≤ α} σα(p) の計算を簡単にするために，情報源 S0として 定常情報源のみ考えることにする．定常性を仮定する ことにより、どの位置からも文字列パターン p の極小 局所最適 k-ギャップ出現が起こる確率は変わらないも のとして計算できる．さらに以下の仮定をおき、σα(p) の近似値を計算するアプローチをとる． 仮定 1 任意の位置 i, j において，位置 i から文字列パ ターン p の極小局所最適 k-ギャップ出現が起こる 事象と位置 j から p の極小局所最適 k-ギャップ出 現が起こる事象は互いに独立である． 定常性とこの仮定により，任意の位置から文字列パター ン p[1..m] の極小局所最適 k-ギャップ出現が起こる確率 を β(p) とすれば，二項分布 B(n− m + 1, β(p)) で p の 極小局所最適 k-ギャップ出現の数 Xpの分布を表すこと が可能である．β(p) が非常に小さな場合には，二項分布 B(n−m+1, β(p)) はポアソン分布 Po((n−m+1)β(p)) で近似することが可能である．このとき，λ = (n−m+ 1)β(p) とすれば ˆ σα(p) = min { σ∈ N 1 1− e−λ ( 1− σ∑−1 i=0 λi i!e −λ ) ≤ α } は σα(p) の近似値となる． 以下のマイニング問題を考える． 問題 1 与えられた非負整数 k，文字列 s, スコア関数 w，小さい正の実数 0 < α < 1 および定常情報源 S0に 対し，極小局所最適 k-ギャップ出現の数が σα(s[h..j]) 以上である s の部分文字列 s[h..j] を列挙せよ． － 46 －
文字列パターン p の s における極小局所最適 k-ギャッ プ出現の数をF(p) とすれば，関数 F は逆単調性を満た さない．つまり，F(p[2..m]) < F(p[1..m]) や F(p[1..m− 1]) <F(p[1..m]) を満たす p は存在する． 出現数の関数が逆単調性を満たす場合、パターンを 少しでも拡張すると出現数が減るパターンを閉パター ン (closed pattern) という．極小局所最適 k-ギャップ 出現は逆単調を満たさないが定義の解釈を広げて，文 字列パターン p を連続部分文字列として含むすべての 文字列 p′に対し，F(p) > F(p′) が成り立つとき，p は 閉文字列パターンであるということにする．
3
無記憶情報源における
β(p)
計算
[7]
情報源 S0が無記憶でギャップを許さない場合に，任 意の文字列パターン p に対して任意の位置から極小局 所最適 0-ギャップ出現が起こる確率 β(p) は，スコア関 数 w が単純なものてある場合，動的計画法により計算 することができる．単純なスコア関数 w として以下の ものを使う． w(x, y) = { +a (x = y) −b (x ̸= y) (1) ただし，a, b は a < b を満たす自然数 (正の整数) と する． 今，簡単のために|Σ| = ℓ とし，情報源 S0において すべての文字は等確率 1/ℓ で発生するものとする．文字 列パターン p[1..m] に対する極小局所最適 0-ギャップ出 現 t[1..m] の不一致の位置の集合 F ≡ {i ∈ {1, ..., m} | p[i]̸= t[i]} の集合族を F(p) とする．すると，S0から 発生した文字列 t[1..m] が文字列パターン p[1..m] の極 小局所最適 0-ギャップ出現である確率 β(p) は， β(p) = m ∑ k=0 |{F ∈ F(p) | |F | = k}| ( 1 ℓ )m−k( 1−1 ℓ )k で計算できる．上式は|{F ∈ F(p) | |F | = k}| の値 を求めれば計算でき，p に依存せず長さ m の文字列パ ターンはすべて同じ値 βmをとる．|{F ∈ F(p) | |F | = k}| は以下のように求めることができる．今，時刻 0 に位置 0 から出発し，時刻 i + 1 に時刻 i の位置から w(p[i + 1], t[i + 1]) だけ移動するランダムウォークを考 えると，時刻 i の位置は S(p[1..i], t[1..i]) に一致する． これはスコア関数 w の定義から，確率 1/ℓ で +a，確 率 1−1/ℓ で −b だけ移動するランダムウォークとなる． 不一致の数が k のとき，最終的な類似スコアは (m− k)a− kb = ma − (a + b)k となる．今，γh,i,kを現在の 位置が i(≤ h) のとき，残り (h − i)/a 回の移動で，途 中で (0, h− (a + b)k) の範囲から出ずに (h − i)/a 回目 Step 1 1≤ h ≤ Ma，0 ≤ i ≤ h を満たす整数のペア (h, i) に対し γh,i,0を (h− i)%a = 0 であれば 1， そうでなければ 0 に設定する．k = 1 とする． Step 2 (a + b)k + 1 ≤ h ≤ Ma の各々の h に対し, i = h− (a + b)k − a を満たす i に対して γh,i,k= γh−(a+b),i−b,k−1に設定する． Step 3 (a + b)k + 1≤ h ≤ Ma の各々の h に対し, i = h−(a+b)k−a−1, h−(a+b)k−a−2, . . . , b+1 の 各々の i に対して大きい方から γh,i,k= γh,i+a,k+ γh−(a+b),i−b,k−1に設定する． Step 4 (a + b)k + 1 ≤ h ≤ Ma の各々の h に対し, i = b, b− 1, . . . , 0 の各々の i に対して大きい方か ら γh,i,k= γh,i+a,kに設定する． Step 5 k + 1 > M a a+b であれば停止．そうでなければ k = k + 1 として Step 2 へとぶ． 図 1: すべての 1 ≤ m ≤ M, 0 ≤ k ≤ ma a+b に対して γma,0,kの値を求める動的計画法のアルゴリズム に h− (a + b)k に到達する一致・不一致のシーケンス パターンの数を表すものとする. このとき γma,0,k =|{F ∈ F(p) | |F | = k}| が成り立つ．γh,i,kに関して，以下の漸化式が成り立つ． 漸化式 1 γh,i,k= k = 0 の場合 1 ((h− i)%a = 0) 0 ((h− i)%a ̸= 0) k≥ 1 の場合 γh,i+a,k (i≤ b) γh,i+a,k+ γh−(a+b),i−b,k−1 (b + 1≤ i ≤ h − (a + b)k − a − 1) γh−(a+b),i−b,k−1 (i≥ h − (a + b)k − a) ただし，% は剰余演算を表すものとする．漸化式 1 を 用いて，図 1 のような動的計画法によりすべての 1≤ m ≤ M, 0 ≤ k ≤ a+bma に対して O(M3) の計算量で αma,0,kの値を求めることが可能である． Σ の文字によって出現確率が異なる場合，最も発生 確率の高い文字の発生確率を qmaxとすれば，任意の文 字列 p[1..m] に対して β(p)≤ m ∑ k=0 |{F ∈ F(p) | |F | = k}|qm−k max (1− qmax)kであることが示されている．上式の右辺を β′mと置き， β(p) の代わりにその上界 βm′ を用いることにより，す べての p に対し，σα(p) 以上の最小サポート σα,mに設 定され，p の極小局所最小 0-ギャップ出現の数を Xpと した場合，PS0{Xp≥ σα,m| Xp≥ 1} ≤ α を満たす．
4
アルゴリズム
ESFLOO-kG.C[4]
問題 1 を解くためには，s の全ての部分文字列に対し て，その極小局所最適 k-ギャップ出現の数を数える必要 がある．与えられた固定した最小サポート σ を使って， 与えられた文字列 s に対し，極小局所最適 k-ギャップ出 現の数が σ 以上である s の部分文字列パターンを列挙す る，現在最も効率的なアリゴリズムは ESFLOO-kG.C である．ここでは，固定した σ ではなく文字列パター ン p に依存した σα(p) を使うように変更したバージョ ンを図 2 に示す． ESFLOO-kG.C では，s のすべての部分文字列を組 織的に生成するために，s の接尾辞木を作成し，根から 深さ優先の木のなぞりを行っている．左境界の局所最 適性 CS2kと右境界の局所最適性 CS3kの両方を満たす ためには，前からと後ろからのアライメントを行う必 要がある [1, 5] が，ギャップの数を合計 k 個に制限した 問題では，位置 h から始まる s の部分文字列が左右両 境界の局所最適性を満たすか否かのチェックを前から の計算のみで行うことが可能である．ESFLOO-kG.C は極小局所最適 k-ギャップ出現候補の左境界の位置 h 毎に小さいテーブル Dh, Hh, Fh, Mhを保持すること により省メモリー，高速化を実現している． 現在のパターンの右端位置 i に対し，−k ≤ ℓ ≤ k, 0≤ u + v ≤ k，u, v ≥ 0 の範囲でテーブル Dh[i, ℓ, u, v] を保持する．Dh[i, ℓ, u, v] は，1 ≤ g ≤ i + 1 および 1 ≤ h′ ≤ h + ℓ + i − 1 − u + v + 1 の範囲におけ る SOPTk−(u+v)(p[g..i], s[h′..h + ℓ + i− 1 − u + v]) の 最大値であり，位置 i で終わる p の部分文字列と位置 h + ℓ + i− 1 − u + v で終わる s の部分文字列の最適 な (k− (u + v))-ギャップアライメントのスコアである． 漸化式は以下のように表される． RF1kC (D h[i, ℓ, u, v] の漸化式) Dh[i, ℓ, u, v]← 0 if i = 0 or h+ℓ +i−1−u+v ̸∈ [1, n], and otherwise max 0 Dh[i− 1, ℓ, u, v + 1] + w(p[i], -) when u + v < k Dh[i, ℓ, u + 1, v] +w(-, s[h + ℓ + i− 1 − u + v]) when u + v < k Dh[i− 1, ℓ, u, v] +w(p[i], s[h + ℓ + i− 1 − u + v]) −k ≤ ℓ ≤ k の範囲の ℓ に対し Dhの定義から Dh[i, ℓ, 0, 0] は、1≤ g ≤ i+1 および 1 ≤ h′≤ h+ℓ+i の範囲におい て SOPTk(p[g..i], s[h′..h + ℓ + i−1]) の最大値となる．し たがって，s の部分文字列 s[h..h+ℓ+m−1] が条件 C2k を満たすための必要十分条件は SOPTk(p[1..m], s[h..h + ℓ + m− 1]) = Dh[m, ℓ, 0, 0] が成り Hh[i, ℓ, u, v] は，S OPTk−(u+v)(p[g..i], s[h′..h + ℓ + i− 1−u+v]) = Dh[i, ℓ, u, v] を達成する位置のペア (h′, g) の中で，辞書順において最大のペアである．Hh[i, ℓ, u, v] は，以下の漸化式を使って効率的に計算できる． RF2kC (Hh[i, ℓ, u, v] の漸化式) Hh[i, ℓ, u, v]← { (h + ℓ− u + v, 1) if i = 0 max Jh[i, ℓ, u, v] otherwise,ただし，
Jh[i, ℓ, u, v]
= {(h+ℓ+i−u+v, i+1) : Dh[i, ℓ, u, v] = 0} ∪{Hh[i− 1, ℓ, u, v + 1] : Dh[i, ℓ, u, v] = Dh[i−1, ℓ, u, v+1]+w(p[i], -)} ∪{Hh[i, ℓ, u + 1, v] : Dh[i, ℓ, u, v] = Dh[i, ℓ, u+ 1, v] +w(-, s[h+ℓ+i−1−u+v])} ∪{Hh[i− 1, ℓ, u, v] : Dh[i, ℓ, u, v] = Dh[i−1, ℓ, u, v]
+w(p[i], s[h+ℓ+i−1−u+v])} 条件 CS1k，CS2k，CS2Mkの３つの条件を満たす s の部分文字列 s[h..j] を (極小局所最小 k-ギャップ) 出現 候補という．出現候補の数に関しては逆単調を満たす． p[1..m] の極小局所最適 k ギャップ出現の可能な最大頻 度は，出現候補 s[h..j] の左端の位置 h の数である．パ ターン p[1..m] に対し，位置 h が出現候補 s[h..j] の開始 位置か否かは，SOPTk(p[1..m], s[h..h + m− 1 + ℓ]) > 0 と Hh[m, ℓ, 0, 0] = (h, 1) がある ℓ∈ [−k, k] に対して成 り立つことである． － 48 －
u, v ≥ 0，u + v ≤ k を満たす u, v に対し，Fh[i, u, v]
は以下のように定義する．ALu,v(p[1..i], s[h..h + i−1+ u− v]) を p[1..i] と s[h..h + i − 1 + u − v] のアライメ ント (p′, s′) で，ギャップ付文字列 p′と s′がそれぞれ u 個，v 個のギャップを含むようなものの集合とする． Fh[i, u, v] を AL u,v(p[1..i], s[h..h + i− 1 + u − v]) に属 するアライメントの最大スコアとする．つまり Fh[i, u, v] = max
(p′,s′)∈ALu,v(p[1..i],s[h..h+i−1+u−v])SALu+v(p
′, s′) と定義する．Fh[i, u, v] は次の漸化式を使って計算で きる． RF3k (Fh[i, u, v] の漸化式) Fh[i, u, v]←
0 if i = 0 and u = v, and otherwise
max −∞ Fh[i− 1, u, v − 1] + w(p[i], -) when i > 0, v > 0 Fh[i, u− 1, v] + w(-, s[h + i − 1 + u − v]) when u > 0, h≤ h + i − 1 + u − v ≤ n Fh[i− 1, u, v] + w(p[i], s[h + i − 1 + u − v]) when i > 0, h≤ h + i − 1 + u − v ≤ n
Fh[i] を，i と h を固定して h を動かしたときの p[1..i] と s[h..j] の k-ギャップアライメントの最大値とする．つ まり，
Fh[i] = max
h+i−1−k≤j<h+i−1+kSOPTk(p[1..i], s[h..j]).
と定義する．Fh[i] の値は Fh[i, u, v] の定義より，式 Fh[i] = max u,v≥0,0≤u+v≤kF h[i, u, v] で計算できる． Mh[m] を Fh[0], Fh[1], ..., Fh[m] の最大値とする．す ると s[h..j] が条件 CS3k満たすための必要十分条件は SOPTk(p[1..m], s[h..j]) = Mh[m] と書ける．Fh[i] が与 えられれば, Mh[i] は以下の漸化式を使って計算できる． RF4k (Mh[i] の漸化式)
Mh[i] = max{Mh[i− 1], Fh[i]}.
s[h..j] が条件 CS3kおよび条件 CS3Mk を満たすた めの必要十分条件は Mh[m] > Mh[m− 1] かつ j = min{h + m − 1 + u − v : Fh[m, u, v] = Fh[m], u, v ≥ 0, 0≤ u + v ≤ k} が成り立つことである． ESFLOO-kG.C は時間計算量が O(k3n3)、空間計算 量が O(k3n2) のアルゴリズムである1． 1ただし k = 0 の場合は時間計算量 O(n3)，空間計算量 O(n2) のアルゴリズムである． ESFLOO-kG.C(s) Input: s[1..n]: 文字列 1: s を大域変数にする． 2: 文字列 s の接尾辞木T を作る． 3: for h = 1 to n do 4: Occ[h]← h 5: for ℓ =−k to k do 6: for u, v≥ 0 s.t. u + v ≤ k do (Dh[0, ℓ, u, v], Hh[0, ℓ, u, v])← (0, (h + ℓ − u + v, 1)) 7: end for 8: for u = 0 to⌊k 2⌋ do F h[0, u, u]← 0 9: for 0≤ u < v ≤ k − u do Fh[0, u, v]← −∞ 10: for 0≤ v < u ≤ k − v do Fh[0, u, v]← Fh[0, u− 1, v] + w(-, s[h − 1 + u − v]) 11: Mh[0]← 0 12: end for 13: occ num← n 14: x← 接尾辞木 T の根節点
15: PatGen(x, 0,(Occ, occ num,{Dh}, {Hh}, {Fh}, {Mh}))
PatGen(x, ℓ,Tbls) Global: s[1..n]: 文字列 Input: x: 文字列 s の接尾辞木の節点 ℓ: 節点 x の深さ Update: Tbls: DPTblCal.Fk.C より更新れるテーブル 1: for 節点 x の子節点 y 各々 do 2: (f, t)← 辺 (x, y) のラベル
3: for e = f to t do DPTblCal.Fk.C(ℓ+e−f+1, s[f−ℓ..e],Tbls)
4: PatGen(y, ℓ + t− f + 1,Tbls)
5: end for
DPTblCal.Fk.C(i, p[1..i],(Occ, occ num,{Dh}, {Hh}, {Fh}, {Mh}))
Global: s[1..n]: Input: i: 現在のパターンの長さ p[1..i]: 現在のパターン Update: Occ: 出現候補の左端のリスト occ num: 出現候補数 {Dh}: Dh[i− 1, ·, ·, ·] に基づく Dh[i,·, ·, ·] の更新 {Hh}: H[i − 1, ·, ·, ·] に基づく Hh[i,·, ·, ·] の更新 {Fh}: Fh[i− 1, ·, ·] に基づく Fh[i,·, ·] の更新 {Mh}: Mh[i− 1] に基づく Mh[i] の更新 1: OC no← 0
2: for j = 1 to occ num do 3: h← Occ[j] 4: flag OC← false 5: for ℓ =−k to k do 6: for u, v≥ 0, 0 ≤ u + v ≤ k do 7: RF1kCと RF2kCを用いて Dh[i, ℓ, u, v] と Hh[i, ℓ, u, v] を計算 8: end for
9: if Hh[i, ℓ, 0, 0] = (h, 1) then flag OC← true
10: end for 11: Fh[i]← −∞
12: for u, v≥ 0, 0 ≤ u + v ≤ k do 13: RF3kを用いて Fh[i, u, v] を計算
14: if Fh[i, u, v] > Fh[i] then Fh[i]← Fh[i, u, v]
15: end for
16: if flag OC = true then 17: RF4kを用いて Mh[i] を計算
18: OC no← OC no + 1, Occ[OC no] ← h
19: end if 20: end for
21: if|{(h, j) : (h, 1) = Hh[i, j−(h+i−1), 0, 0], Mh[i] > Mh[i−1],
j = min{h + i − 1 + u − v : Fh[i, u, v] =
Fh[i]}}| ≥ σα(p[1..i]) then
22: Print p[1..i] 23: end if 24: occ num← OC no 図 2: Algorithm ESFLOO-kG.C
5
全ヒトゲノムを扱うための高速化
5.1
同じ部分文字列の出現の共有
現在，同じ部分文字列で異なる位置に出現する出現 候補を別々に計算している．これは特にパターン文字が短い場合に，同じ部分文字列が多く出現するため無 駄が多くなっている．そこで，接尾辞木を使った同じ 部分文字列の出現候補の共有を次のように実現する． • 接尾辞木の各節点に，その子孫の葉の数を記録 し，共有出現候補の出現回数をすぐに計算できる ようにする． • パターン長が m であれば，一つの出現候補の左 端位置 h で極小局所最適 k-ギャップ出現であるか 否かを判断するために見なければならない s の範 囲は，[h− 2k, h + m + 2K] となり，その長さは m + 4k + 1 である．接尾辞木において，根からこ の長さの文字列に対応する枝の位置の数だけ異な るものがあり，それだけのものに対して出現候補 であるか，また実際の出現であるかをチェックす ればよい．ただし，長さが m− k 以上であれば出 現である可能性があるので s の両端に関しては， 長さが [m− k, m + 4k] のものもチェックする． • ESFLOO-kG.C では，パターンを拡張しながら (出現候補は逆単調性を満たすので) コピーして出 現候補を引き継ぐことを行う．パターンの拡張に よりパターン長が m + 1 になると，出現候補と してみなければならない範囲も 1 増える．出現候 補のチェックすべき範囲がちょうど接尾辞木の節 点の位置であるとき、拡張すると分岐するのでそ のような場合は出現候補を分割する．
5.2
並列探索
s の接尾辞木をなぞりながらパターンを生成するが， 各節点の子節点は別パターンに対応し全く独立に頻出 チェックを行うことが可能である．また，見つかった頻 出パターンは，他の頻出パターンを見つけるのに全く 必要ない．そこで各節点ではデタッチ属性でスレッドを 立ち上げて並列処理する．スレッドの数を制限するた めにセマフォーを用い，スレッド生成に sem trywait で セマフォーの減算を行い，スレッド消滅時に sem post でセマフォーの加算を行ことによりスレッド数を指定 した数に抑える．また，結果の出力はスレッド毎別ファ イルを割り当て独立に行う．6
むすび
現在，第 5 節で述べた高速化を C++言語で EFLOO-0G.C に組み込み，Intel(R) Xeon(R) CPU E7-2850(2.00GHz) 80 コア (160 スレッド), 3TB メモリのマシンにおいて， 64 スレッドを使って，無記憶情報源 S0に対するパター ン長 m 依存の最小サポート σα,mを，α = 1.0× 10−20 に対して計算して用い，全ヒトゲノムから有意な散在 反復配列を抽出する実験を行っている．参考文献
[1] Bruce W. Erickson and Peter H. Sellers. Recogni-tion of patterns in genetic sequences. In David Sankoff and Joseph B. Kruskal, editors, Time
Warps, String Edits, and Macromolecules : The Theory and Practice of Sequence Comparison,
chapter 2, pages 55–91. Addison-Wesley, Reading, MA, 1983.
[2] G. J. Faulkner et al. The regulated retrotrans-poson transcriptome of mammalian cells. Nature
Genetics, 41:563–571, 2009.
[3] Ruiqiang Li, Jia Ye, Songgang Li, Jing Wang, Yu-jun Han, Chen Ye, Jian Wang, Huanming Yang, Jun Yu, Gane Ka-Shu Wong, and Jun Wang. Reas: Recovery of ancestral sequences for trans-posable elements from the unassembled reads of a whole genome shotgun. PLOS Computational
Biology, 1(4):e43, 2005.
[4] A. Nakamura, I. Takigawa, H. Tosaka, M. Kudo, and H. Mamitsuka. Mining approximate patterns with frequent locally optimal occurrences.
Dis-crete Applied Mathematics, 200:123–152, 2016.
[5] Jeanette P. Schmidt. All highest scoring paths in weighted grid graphs and their application to finding all approximate repeats in strings. SIAM
Journal on Computing, 27(4):972–992, 1998.
[6] Feida Zhu, Xifeng Yan, Jiawei Han, and Philip S. Yu. Efficient discovery of frequent approximate se-quential patterns. In Proceedings of the 7th IEEE
International Conference on Data Mining (ICDM 2007), pages 751–756, 2007.
[7] 中村武則, 中村篤祥, 工藤峰一. ランダム仮説下で の出現回数分布を利用した散在反復配列の発見. In
DEIM Forum 2016 予稿集, pages G5–3, 2016.
