最適アラインメント計算装置及びプログラム
【課題】リード長に制限を設けずに、計算機のワード演算のビット並列性を利用して、動的計画法（DP法）によるアラインメント処理を１プロセッサで並列化し高速化する。
【解決手段】DP行列を計算する範囲を、計算機のワード長に対応する幅をもつ対角線周辺領域に制限し、その領域内の情報を対角線と直行する方向のビットの並びに分割して表現し、それらの情報に対する計算をビット並列化する。
【解決手段】DP行列を計算する範囲を、計算機のワード長に対応する幅をもつ対角線周辺領域に制限し、その領域内の情報を対角線と直行する方向のビットの並びに分割して表現し、それらの情報に対する計算をビット並列化する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、文字列、DNA塩基配列等の文字列データ（文字列データとして取り扱い可能な形式に加工された任意のデータを含む。）を、検索文字列と比較照合する処理を高速化するための技術に関する。
【背景技術】
【０００２】
近年、従来型のキャピラリー型ＤＮＡシーケンサとは異なる原理に基づき、一回のランで数千万本以上にも及ぶ大量の配列を読み取ることができる超並列ＤＮＡシーケンサが出現し、シーケンシングに要するコストを劇的に低減させた（非特許文献１）。そのため、大量の配列データを高速に解析する技術に対する需要が高まりつつある。中でも、ゲノム・マッピング解析は、最も基本的な配列データ解析の一つであり、シーケンサにより新たに読まれた塩基配列（リード）を、データベース中の参照ゲノム配列と比較して、ゲノム内の位置を特定し、さらに、そこに含まれる変異（塩基の置換・挿入・欠失）を同定する。
【０００３】
初期の超並列ＤＮＡシーケンサのリード長は30塩基程度と比較的短かったため、ハッシュ法やBurrows-Wheeler変換などを用いて高速・高精度なゲノム・マッピング解析がなされた（非特許文献２）。
【０００４】
しかし、リード長が100塩基以上まで伸びてくると、従来型のＤＮＡシーケンサ向けのゲノム・マッピング解析で用いられてきたseed-and-extend strategyのような方法が必要となってくる（非特許文献３）。この方法では、第一ステップで、リードに含まれる短い配列（シード）がゲノム内に出現する位置を高速に検索し、第二ステップで、シードの出現位置周辺でリード全体がゲノム配列に類似しているか（即ち、シードとゲノムの対応関係を、多少の変異を許容しながら、リード全体にまで拡張できるか）を調べる。この第一ステップでは、シード長が短いため、ハッシュ法やBurrows-Wheeler変換などを用いた高速な検索方法が利用でき、これにより、長大なゲノム配列全体の中からリードのマッピング候補位置を高速に絞り込むことができる。第二ステップでは、第一ステップで絞り込まれた候補位置にあるゲノム部分配列（ターゲット配列）に対して、動的計画法（DP法）などを用いてリード（クエリー配列）とのアラインメントを行い、リード長が長い場合であっても、塩基変異を高精度に調べることができる（非特許文献８）。
【０００５】
しかし、一般に、動的計画法では、配列長の2乗に比例する大きな計算コストが必要となる。そこで、動的計画法によるアラインメント処理を高速化することは、ゲノム・マッピング解析処理全体を高速化する上で有益となる。
【０００６】
動的計画法によるアラインメント処理を高速化するために、複数のプロセッサを利用して処理を並列化する方法が知られている（非特許文献４）。このとき、プロセッサ間で受け渡されるデータの依存関係を考慮し、また、プロセッサ間通信コストを抑えるような負荷分散方式を採用する必要がある。それにより、多数のプロセッサを使えば、高い並列性が得られる。
【０００７】
ところで、現在主流のデジタル計算機では、8バイト整数（ワード）に対する演算が基本演算となっているが、これは、8バイト即ち64ビットの並列演算と見ることもできる。G. Myersは、この整数演算のビット並列性を利用して、動的計画法によるアラインメント処理を１プロセッサで並列化し、配列長に比例する時間で計算する方法を示した（非特許文献５）。
【０００８】
G. Myersの方法では、長さmのクエリー配列（リード）と長さnのターゲット配列（ゲノム部分配列）のアラインメントを求めるために、サイズm×nのDP行列の情報を長さmの列方向の情報のn個の並びに分解して表現し、新たに読み込んだターゲット配列の１文字とクエリー配列の全m文字との比較に関わる処理を、mビット並列に行い、行方向（ターゲット配列の方向）に1塩基ずつ移動しながらこの処理を繰り返す。
【０００９】
式を用いて詳細に説明すると、次のようになる。１塩基の置換・挿入・欠失のコストを全て１として、DP行列の第 (i, j) 成分C[i, j] （クエリー配列 q = q1 q2 ... qm の長さ iの部分文字列q1 q2 ... qi とターゲット配列 t = t1 t2 ... tn の長さ jの部分文字列 t1 t2 ... tjとの最適なアラインメントのコスト）に対して、
等式C[i, j]−C[i, j−1] = 1が成立するか否かを表わすブール変数をPh[i, j]、
等式C[i, j]−C[i, j−1] = −1が成立するか否かを表わすブール変数をMh[i, j]、
等式C[i, j]−C[i−1, j] = 1が成立するか否かを表わすブール変数をPv[i, j]、
等式C[i, j]−C[i−1, j] = −1が成立するか否かを表わすブール変数をMv[i, j]、
クエリー配列のインデクス iの文字（i番目の文字）qiとターゲット配列のインデクスjの文字tjとが一致するか否かを表わすブール変数をEq(i, j)、
と表わし、これらの 1ビット（偽または真、即ち、0または1）で表わされるブール変数をi = 1, 2, ..., m について下位ビットより昇順に並べてできるワード変数を、
（数１)
Ph(j) = [ Ph[m, j], ・・・ , Ph[2, j], Ph[1, j] ]
Mh(j) = [ Mh[m, j], ・・・ , Mh[2, j], Mh[1, j] ]
Pv(j) = [ Pv[m, j], ・・・ , Pv[2, j], Pv[1, j] ]
Mv(j) = [ Mv[m, j], ・・・ , Mv[2, j], Mv[1, j] ]
Eq(j) = [ Eq[m, j], ・・・ , Eq[2, j], Eq[1, j] ]
と表わし、さらに記号の簡略化のためにjを一つの値に固定して、Ph(j), Mh(j), Pv(j), Mv(j), Eq(j) の値をそれぞれ、Ph, Mh, Pv, Mv, Eqと略記し、Xv、Xhを補助ワード変数として、
（数２)
Xh = ((Eq & Pv) + Pv) ^ Pv | Eq
Ph = Mv | ~ (Xh | Pv)
Mh = Pv & Xh
Ph <<= 1
Mh <<= 1
Xv = Eq | Mv
Pv = Mh | ~(Xv | Ph)
Mv = Ph & Xv
によりビット並列な整数演算を行い、Ph, Mh, Pv, Mvの値を上書きする。ここで、ワード変数a、bに対して、a|bはビットごとのOR演算、a&bはビットごとのAND演算、a ^ bはビットごとの排他的論理和演算、a + bは整数としての加算、~aはビットごとのNOT演算、a<<=1は1ビットの左シフト演算を表わす。これにより、Ph, Mh, Pv, Mvの値は、Ph(j+1), Mh(j+1), Pv(j+1), Mv(j+1) の値に更新される。これが、アラインメント処理をビット並列化する中核部分であり、これにより処理が高速化される。
【００１０】
しかし、この方法を直接的に用いることができるのは、リード長（クエリー配列長）m が計算機のワード長（64）以下の場合に限られる。この方法をリード長が計算機のワード長を越える場合に拡張するためには、リード（クエリー配列）をワード長以下の長さの複数の断片配列に分割する必要がある（非特許文献５、非特許文献６）。このとき、これらの断片配列に対する計算には互いに依存関係があり、並列化する際に制約が生じ、処理効率も低下する。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】Rajendrani Mukhopadhyay. DNA sequencers: the next generation. Analytical Chemistry, Vol. 81, No. 5, March 1, 2009
【非特許文献２】Kouichi Kimura, Yutaka Suzuki, Sumio Sugano, Asako Koike. Computation of Rank and Select Functions on Hierarchical Binary String and Its Application to Genome Mapping Problems for Short-Read DNA Sequences. Journal of Computational Biology. November 2009, 16(11): 1601-1613.
【非特許文献３】Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005 May 1;21(9):1859-75.
【非特許文献４】Martins WS, Del Cuvillo JB, Useche FJ, Theobald KB, Gao GR. A multithreaded parallel implementation of a dynamic programming algorithm for sequence comparison. Pac Symp Biocomput. 2001:311-22.
【非特許文献５】G. Myers. A fast bit-vector algorithm for approximate string matching based on dynamic programming. Journal of the ACM, 46(3):395-415, 1999.
【非特許文献６】HEIKKI HYYRO. BIT-PARALLEL COMPUTATION OF LOCAL SIMILARITY SCORE MATRICES WITH UNITARY WEIGHTS. International Journal of Foundations of Computer Science, 17(6): 1325-1344, 2006.
【非特許文献７】Kevin Judd McKernan, Heather E. Peckham, Gina L. Costa, et al. Sequence and structural variation in a human genome uncovered by short-read, massively parallel ligation sequencing using two-base encoding. Genome Res. 2009 19: 1527-1541.
【非特許文献８】Richard Durbin, Sean R. Eddy, Andrew Krogh, Graeme Mitchison, Biological Seuqnece Analysis, probabilistic models of proteins and nucleic acids, Cambridge University Press 1988（邦訳、阿久津達也、浅井潔、矢田哲士訳、バイオインフォマティクス―確率モデルによる遺伝子配列解析―、医学出版、2001）
【発明の概要】
【発明が解決しようとする課題】
【００１２】
本発明が解決しようとする課題は、リード長に対する制限を設けることなく、また、リード配列を分割するなどの処理効率を低下させる方法もとらず、計算機のワード演算のビット並列性を利用して、動的計画法によるアラインメント処理を1プロセッサで並列化し高速化する方法を提供することにある。
【課題を解決するための手段】
【００１３】
本発明による方法では、長さmのクエリー配列（リード）と長さnのターゲット配列（ゲノム部分配列）のアラインメントを求めるために、サイズm×nのDP行列の計算を対角線（行のインデクスiと列のインデクスjが等しいところ）の周辺領域に制限し、その領域内の情報を対角線に直行する方向（反対角線方向）の情報のビットの並びに分解して表現し、計算機のワード演算のビット並列性を利用して、反対角線上でのクエリー配列とターゲット配列の文字比較に関わる処理を1プロセッサで並列化し、対角線方向に移動しながらこの処理を繰り返す。
【００１４】
式を用いて詳細に説明すると、次のようになる。C[i, j], Ph[i, j], Mh[i, j], Pv[i, j], Mv[i, j], Eq[i, j] を上記の通りとして、これらを対角線周辺で反対角線方向に並べてできるワード変数を、
（数３）
Ph(k) = [ Ph[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Mh(k) = [ Mh[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Pv(k) = [ Pv[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Mv(k) = [ Mv[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Eq(k) = [ Eq[i, j] | i + j = k, −m ＜ i−j ≦ m ]
とする。ここで、ビット変数をワード変数に詰め込む順番は、jについて下位ビットから昇順とする。さらに、記号の簡略化のためにkを一つの値に固定して、Ph(k), Mh(k), Pv(k), Mv(k), Eq(k) の値をそれぞれ、Ph, Mh, Pv, Mv, Eqと略記し、Xv、Xhを補助ワード変数として、
kが偶数のときは、
（数４）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
kが奇数のときは、
（数５）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' << = 1; Mv' <<= 1;
に従ってビット並列な整数演算を行い、Ph', Mh', Pv,' Mv'の値を計算する。これらは、Ph(k+1), Mh(k+1), Pv(k+1), Mv(k+1) の値を与える。これが、アラインメント処理をビット並列化する中核部分である。
【発明の効果】
【００１５】
G. Myersの方法と同様、アラインメント処理が１プロセッサでビット並列化され高速化される。G. Myersの方法と異なるのは、動的計画法のビット並列計算を、DP行列の対角線周辺領域に限定して対角線方向に移動させるため、リード長に関する制限が無くなり、従って、リード長がワード長を越えた場合にリードを分割する必要がなくなる点である。また、類似性の高い配列の対する動的計画法の計算は、DP行列の対角線の周辺を行えば十分であるため、対角線から離れた領域での不要な計算を避けることができ、処理効率が向上する。また、行と列に関して対称な形で計算を行えるような新規なビット並列化方式であるため、（数４）及び（数５）において補助変数Xhの計算式はXvの計算式と対称な形に簡易化され、これはG. Myersによる（数２）における補助変数Xhのやや複雑な計算式よりも少ないマシン命令数で計算可能となる。
【図面の簡単な説明】
【００１６】
【図１】本発明による最適アラインメント計算方法において、ビット並列演算の適用箇所を示した説明図である（実施例１）。
【図２】次世代型DNAシーケンサの配列データと参照ゲノム配列データベースとの高速比較照合を行うシステムの構成例を示す図である（実施例１）。
【図３】DP行列内のセル周辺における入力側と出力側のブール変数の位置関係を示す説明図である。
【図４】和が k となる2重インデクス [i, j] とワード変数のインデクスとの対応関係を例示する説明図である。
【図５】一定の k に対するブール変数Ph, Mh, Pv, Mv のインデクスとワード変数のインデクスの対応関係を例示する説明図である。
【図６】セル周辺のブール変数Ph, Mh, Pv, Mvのインデクスが、k が偶数か奇数かにより違いが生じることを示す説明図である。
【図７】ブール変数 Eq をビット並列に計算する方法を示す説明図である。
【図８】対角線周辺領域の外部で事前に設定されるブール変数値を示す説明図である。
【図９】DPコストの絶対値を求めるためのDPコストの差分値の累積法を示す説明図である。
【発明を実施するための形態】
【００１７】
［実施例１］
以下、本発明の実施例を図面を用いて詳細に説明する。
本実施例では、リード（クエリー配列）と参照ゲノム配列データが与えられたとき、リードに類似した部分配列が存在するゲノム配列内の位置（マッピング位置）を求め、クエリー配列とマッピング位置周辺のゲノム部分配列（ターゲット配列）との最適アラインメントを求めるための方法を説明する。通常、DNAシーケンサから得られるリード配列は塩基配列であるため、以下、リード配列は塩基配列であると仮定する。特殊な場合として、２ベース・エンコーディング方式を採用したDNAシーケンサでは、リード配列は、隣り合う２塩基を4色に符号化したカラー配列となる（非特許文献７）。このような場合は、参照ゲノム配列をカラー配列に変換して、本実施例で述べる方法に類似の方法を適用することは、容易な類推により可能である。また、アミノ酸配列をクエリー配列として、参照タンパク配列データベースを検索する課題に対しても、容易な類推により、本実施例で述べる方法に類似の方法が適用可能である。また、文字列や画像を含む文書その他のデータベースを検索する課題に対しても、容易な類推により、本実施例で述べる方法に類似の方法が適用可能である。
【００１８】
図2は、本実施例のシステム構成例を示す図である。201は、このシステムが稼働している計算機であり、外部記憶装置202に格納された参照ゲノム配列データ203の情報を主記憶209内に読み込み、またＤＮＡシーケンサ204を用いてＤＮＡサンプル205を解析して得られた塩基配列データ（クエリー配列、リード配列）206を、入力装置207を介して主記憶209内に読み込み、また、ユーザにより指定される処理パラメータ（類似配列のアラインメント・コストの上限値M）を入力装置208を介して主記憶209内に読み込む。
【００１９】
計算機201は、主記憶209 に蓄えられたリード配列と参照ゲノム配列に対して、シード検索処理210を行い、リードに含まれる短い配列（シード配列）がゲノム配列内に出現する位置を求める。シード配列としては、リード配列の中から長さ12〜20塩基程度の短い部分配列を切り出したものを使えばよく、また、シード配列の出現位置を求めるためには、Burrows-Wheeler変換など用いた公知の高速な検索方法を用いればよい（非特許文献２）。
【００２０】
詳細比較処理部211は、シード配列がゲノム配列内に出現する位置の周辺で、最適アラインメント計算212を行い、そのアラインメントのコスト値がユーザにより指定された上限値M以下ならば、そこにリード配列に類似した配列があると判断し、出力処理部213にゲノム配列内の出現位置情報とアラインメント情報を伝える。出力処理部213は、出力装置214に、検索結果としてゲノム配列内の出現位置情報（ゲノム・マッピング位置情報）215と、最適アラインメント結果216を出力する。
【００２１】
アラインメントのコストは、アラインメントに内に含まれる１塩基の置換、挿入、欠失のコストを１として、その総和として計算する。クエリー配列のp塩基目の位置から始まる長さsのシード配列が、ゲノム配列中の座標区間 [x, x + s−1] に出現している場合、リード配列の出現位置は、座標区間 [x−p+1−M, x−p+m+M] に含まれる。ここで、mはクエリー配列長、Mはアラインメント・コストの上限値を表わす。リード配列の出現位置の範囲が±Mだけ広がるのは、シード配列の出現位置の前後に、最大M塩基の挿入・欠失が含まれる可能性があるからである。最適アラインメント計算212では、この座標区間 [x−p+1−M, x−p+m+M] のゲノム部分配列をターゲット配列として、クエリー配列とターゲット配列との間で、セミ・グローバルなアラインメントを行う。即ち、クエリー配列の全体と、ターゲット配列の部分配列との間に、挿入・欠失を許しながら対応関係を定める。その計算は、具体的には、DP行列を計算することによってなされる（非特許文献８）。
【００２２】
クエリー配列長をm、ターゲット配列長をnとする。１≦ i ≦m，１≦ j ≦nに対して、クエリー配列 q = q1 q2 ... qm の長さiの部分文字列 q1 q2 ... qi とターゲット配列 t = t1 t2 ... tn の長さjの部分文字列 t1 t2 ... tj との最適なアラインメントのコストを表わすDP行列の第 (i, j) 成分を C[i, j] として、G. Myersに従い（非特許文献５）、
等式C[i, j]−C[i, j−1] = 1が成立するか否かを表わすブール変数をPh[i, j]、
等式C[i, j]−C[i, j−1] = −1が成立するか否かを表わすブール変数をMh[i, j]、
等式C[i, j]−C[i−1, j] = 1が成立するか否かを表わすブール変数をPv[i, j]、
等式C[i, j]−C[i−1, j] = −1が成立するか否かを表わすブール変数をMv[i, j]、
クエリー配列のインデクスiの文字（i番目の文字）qi とターゲット配列のインデクスjの文字（j番目の文字） tj が一致するか否かを表わすブール変数をEq(i, j)、
と表わす。また、セミ・グローバルなアラインメントに対応する境界条件として、C[i, 0] = i, C[0, j] = 0と定める。任意の１≦ i ≦m，１≦ j ≦nに対して、
（数６）
Mh = Mh[i−1, j]; Ph = Ph[i−1, j];
Mv = Mh[i, j−1]; Pv = Pv[i, j−1];
Mh' = Mh[i, j]; Ph' = Ph[i, j];
Mv' = Mv[i, j]; Pv' = Pv[i, j];
Eq = Eq[i, j];
と略記すると、これらブール変数の間には、
（数７）
Xv = Eq or Mv; Pv' = Mh or not(Xv or Ph); Mv' = Ph and Xv;
Xh = Eq or Mh; Ph' = Mv or not(Xh or Pv); Mh' = Pv and Xh;
の関係式が成立する（非特許文献５）。ここで、XvとXhは補助のブール変数を表わす。この関係式を反復して利用することによりC[i, j]の値を全て計算することができ、DP行列の下辺でC[m, j]（j = 1, 2, ... , n）の最小値を与える点 (m, j) からトレースバックすることにより最適なアラインメントを求めることができる（非特許文献８）。
【００２３】
ここで、複数組のPh, Mh, Pv, Mvから複数組のPh', Mh', Pv', Mv'を同時にビット並列に計算して、処理を高速化する方法を考える。図３に、これらのブール変数のDP行列内部における位置関係を示す。これらのブール変数は、DP行列内のセル109の周りにおいて、水平方向または垂直方向に隣接するDP行列の格子点301におけるコスト値 C[i, j] の差分に関する情報を表現している。（数６）より、入力側302のブール変数Ph, Mh, Pv, Mvの2重インデクスは [i−1, j] または [i, j−1] でその和は i + j−1 であり、また、出力側303のブール変数Ph', Mh', Pv', Mv'の2重インデクスは [i, j] でその和は i + j である。従って、2重インデクス [i, j] の和 k が小さいものから順に計算することによって、全てのブール変数Ph, Mh, Pv, Mv の値を計算することができる。さらに、クエリー長に関する制限を設けることなく、これらの計算をビット並列化できる方法があることを以下に説明する。
【００２４】
ところで、和が k に等しい2重インデクス [i, j] の個数は、 k と共に増大するが、類似配列の最適アラインメントの計算でトレースバックする際に参照すべきブール変数の2重インデクスは対角線 i = j の周辺の範囲に限られる。そこで、この範囲に制限し、（数７）の関係式をビット並列に計算する方法を考える。
【００２５】
計算機の基本演算に用いられる整数ワードのビット長をwとする。wは一般に２のべき乗であり、現在主流の計算機では、w = 64 となっている。また、半ワード長をw' = w／2とする。対角線の周辺の範囲として、−w < i−j≦w の範囲を考えると、和が k となる2重インデクス [i, j] は、
k が偶数 k = 2k' の場合、
（数８）
[i, j] = [k' + w', k'−w'], ..., [k'+1, k'−1], [k', k'], [k'−1, k'+1],
..., [k'−w' + 1, k' + w'−1]
k が奇数 k = 2k' + 1 の場合、
（数９）
[i, j] = [k' + w', k'−w' + 1], ..., [k'+1, k'], [k', k'+1],
..., [k'−w' + 1, k' + w']
となる。そこで、これらの2重インデクス（格子点）に、（1重）インデクスα= k'−i を対応させれば、kが偶数であるか奇数であるかに関わらず、これらは、
α = −w', −w' + 1, ..., −1, 0, 1, ..., w'−1
に対応する。そこで、wビットのワード変数の各ビットを、下位ビットより順にこれらでインデクス付けして対応させる。以下、説明を簡単にするために w = 4 として、図を用いてこれらの対応関係を具体的に説明する。
【００２６】
図４は、w = 4 の場合に、対角線の周辺範囲 −w < i−j ≦w にある和が k となる2重インデクス [i, j] と4ビットのワード変数のインデクスとの対応関係を例示したものである。101は、DP行列の2重インデクス [i, j] の全体を表わす正方格子である。102と103に表示したインデクスi、jは、各々、クエリー配列とターゲット配列のインデクス（文字位置）である。104で示した対角線 i = j の周辺領域−w < i−j ≦wは、105に示す2本の破線に挟まれた範囲106である。この範囲106の中で、2重インデクスの和がk = 5（奇数）となるものは、407に示す破線で囲まれた中にある4個の白丸（408）の位置にある。同様に、この範囲106の中で、2重インデクスの和がk = 6（偶数）となるものは、409に示す破線で囲まれた中にある4個の黒丸（410）の位置にある。k = 7（奇数）の場合も同様である。これらに対応するインデクスα = −2, −1, 0, 1 の値を傍らに付した。また、これらのインデクスに対応する4ビット整数のビット位置を411に示した。
【００２７】
和が k となる2重インデクス [i, j] をもつブール変数Ph, Mh, Pv, Mv に対しても、同様にワード変数の（1重）インデクスαを対応させる。ブール変数Ph, Mhは水平方向に隣接する格子点のコストの差分値の情報を表わしているが、それらには左方の格子点のインデクスαを対応させる（図5の501参照）。即ち、Ph[i, j], Mh[i, j] には [i, j−1] に対応するインデクスαを対応させて、Ph[α], Mh[α] と表わす。同様に、ブール変数Pv, Mvは垂直方向に隣接する格子点のコストの差分値の情報を表わしているが、それらには下方の格子点のインデクスαを対応させる（図5の502参照）。即ち、Pv[i, j], Mv[i, j] には [i, j] に対応するインデクスαを対応させて、Pv[α], Mv[α] と表わす。図5に、w = 4 の場合の例を具体的に示す。h[α]はブール変数Ph[α]またはMh[α]を表わし、また、v[α]はブール変数Pv[α]またはMv[α]を表わす。
【００２８】
対角線周辺領域106内に、一定の k に対する Ph[α], Mh[α], Pv[α], Mv[α] は各々w個ずつあるので、対応するインデクスαに従いそれらのビット値503を1ワードに詰め込むことができ、そのようにして得られるワード変数504, 505をPh, Mh, Pv, Mv と表わす。図5はw = 4 の場合を示す。ここで、hはPhまたはMhを、vはPvまたはMvを表わす。ブール変数Eq = Eq[i, j]についても、同様に、2重インデクス [i, j] に対応するインデクスαを対応させて Eq[α] と表わし、それに従ってビット値をワード変数に詰め込んだものをEqと書く。
【００２９】
これらのワード変数は、 k について昇順に計算することができる。図1は、w = 4 の場合に、その計算過程におけるワード変数のデータの依存関係と、その時に計算されるブール変数のDP行列内における位置を示したものである。対角線周辺領域106（−w < i−j ≦w）の中で、kが奇数のとき、ブール変数Ph[α], Mh[α], Pv[α], Mv[α]に対応する水平または垂直方向の隣接格子点を結ぶ辺を太い実線107で、また、kが偶数のとき、対応する同様の辺を太い点線108で表わした。また、辺で囲まれた各セル109には、ブール変数Eq[α]が対応している。これらのブール変数をワード変数に詰め込んだものがPh, Mh, Pv, Mv, Eqであり、これらは110に示すワード・データを構成する。これらは、以下に説明するビット並列演算111によりkについて昇順に計算される。
【００３０】
ビット並列演算111を説明するために、図5に戻って、各セル109の周りで計算すべき式を確認する。但し、セルの左辺と上辺に対応するkの値が偶数か奇数かによって事情が異なる。先ず、kが偶数の場合は、図6左に示すように、kに対応する左辺と上辺の入力側ブール変数601はPv[α], Mv[α] , Ph[α], Mh[α] であり、また、k+1に対応する下辺と右辺の出力側ブール変数602はPh[α], Mh[α], Pv[α], Mv[α]であり、これらは全て同じインデクスαをもっている。従って、この場合は、ワード変数 Ph, Mh, Pv, Mv に対して、kに対応するものをPh, Mh, Pv, Mv、また、k+1に対応するものをPh', Mh', Pv', Mv'と書いて区別すれば、（数７）と同様の関係式、
（数４）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
により、ビット並列に計算できる。ここで、ワード変数a、bに対して、a|bはビットごとのOR演算、a&bはビットごとのAND演算、~aはビットごとのNOT演算を表わす。一方、kが奇数の場合は、図6右に示すように、kに対応する左辺と上辺の入力側ブール変数603はPv[α], Mv[α] , Ph[α+1], Mh[α+1] であり、また、k+1に対応する下辺と右辺の出力側ブール変数604はPh[α], Mh[α], Pv[α+1], Mv[α+1] であり、インデクスαに+1のずれが生じている。そこで、この場合は、上記のワード変数 Ph, Mh, Pv, Mv 及び Ph', Mh', Pv', Mv' に対して、前後にビットシフト演算を行えば、（数７）と同様な関係式が適用可能となる。即ち、この場合は、
（数５）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' <<= 1; Mv' <<= 1;
により、ビット並列に計算できる。ここで、ワード変数aに対して、a<<=1 は1ビットの左シフト演算、a>>=1 は1ビットの右シフト演算を表わす。
【００３１】
一方、ワード変数 Eq については、
k が偶数 k = 2k' の場合、
（数10）
Eq(2k') = [ Eq[−w'], ..., Eq[−1], Eq[0], Eq[1], ..., Eq[w'−1 ] ]
= [ Eq[k'−w' + 1, k' + w'−1], ...,
Eq[k'−1, k'+1], Eq[k', k'], Eq[k'+1, k'−1],
..., Eq[k' + w', k'−w'] ]
k が奇数 k = 2k' + 1 の場合、
（数11）
Eq(2k'+1) = [ Eq[−w'], ..., Eq[−1], Eq[0], Eq[1], ..., Eq[w'−1 ] ]
= [ Eq[k'−w' + 1, k' + w'], ...,
Eq[k', k'+1], Eq[k'+1, k'],
..., Eq[k' + w', k'−w' + 1] ]
となる。ここで、外側の[ , , ... , ] はビットを並べてできる整数値を表わし（右端が最下位ビット）、Eq[i, j] は、クエリー配列 q = q1 q2 ... qm のi番目の文字（塩基）qiとターゲット配列 t = t1 t2 ... tn のj番目の文字（塩基）tjが一致するか否かを表わすブール変数を表わす。4種類の塩基A, C, G, Tの2ビット符号化 Aa00, Ca01, Ga10, Ta11 における下位ビットをl、上位ビットをu と表わし、l(A) = l(G) = 0, l(G) = l(T) = 1, u(A) = u(C) = 0, u(G) = u(T) = 1 とすると、
（数12）
Eq[i, j] = not( ( l(qi) xor l(tj) ) or ( u(qi) xor u(tj) ))
が成り立つ。ここで、xorはブール変数の排他的論理和を表わす。そこで、ブール変数l(tj)、u(tj)を最下位ビットから昇順に並べてできるワード変数をTl, Tu、ブール変数l(qi)、u(qi)を最上位ビットから降順に並べてできるワード変数をQl, Quとすれば、
ワード変数Eqは、
（数13）
Eq = ~( ( Ql ^ Tl ) | ( Qu ^ Tu ))
のようにビット並列に計算できる。ここで、ワード変数a, bに対して、a ^ bはビットごとの排他的論理和演算を表わす。一方、Tl, Tu, Ql, Quは、k について帰納的にビットシフト演算を用いて計算できる。但し、ビットシフトによって空いたビットには、クエリー配列またはターゲット配列の対応する文字の符号化ビットをセットする。図7にw = 4 の場合のk = 4, 5, 6に対応するTl, Tu, Ql, Quの計算方法を示す。一般の場合に、kが増加するときのTl, Tu, Ql, Quの値の更新方法を具体的に式で表わすと、
kが奇数（k = 2k' + 1）のとき、
（数14）
Ql <<= 1; Ql |= l(qi);
Qu <<= 1; Qu |= u(qi);
kが偶数（k = 2k'）のとき、
（数15）
Tl >>= 1; Tl |= (l(tj) << (w−1));
Tu >>= 1; Tu |= (u(tj) << (w−1));
となる。ここで、qi, tjはクエリー配列とターゲット配列から新たに読み込まれる文字を表わし、i = w' + k' + 1、j = w' + k' である。また、ワード変数a, b に対して、a = a | b を a |= b と略記する。
【００３２】
上述の計算において、対角線周辺領域106の境界105 （ i−j = ±w ）に内部から接する辺に対応するブール変数Pv[α], Mv[α] , Ph[α], Mh[α] を計算する際、この境界に外部から接する辺に対応するブール変数Pv[α], Mv[α] , Ph[α], Mh[α] の値を参照する必要が生じる。後者の値は計算されないので、境界条件として、それらに適切な値を事前に定めておく必要がある。そこで、対角線周辺領域106の外部では、クエリー配列とターゲット配列の文字比較結果は常に不一致と仮定する。これは、偽の類似配列が求められることを避けるための（悲観的）仮定である。このことは、対角線周辺領域106の外部では C[i, j]−C[i−1, j−1] = 1 が常に成り立つことを意味する。そのためには、図8に示すように、対角線周辺領域の外部では常に Pv[α] = 1, Mv[α] = Ph[α] = Mh[α] = Eq[α] = 0 が成り立つと仮定すればよい。そこで、これを境界条件とする。このとき、（数７）が成立することは、これらの値を代入して直接確かめることができる。
【００３３】
また、k の値がワード長 w 以下のとき、ワード変数Ph, Mh, Pv, Mv は、DP行列外部の i < 0，j < 0 となる領域に対応するブール変数Pv[α], Mv[α] , Ph[α], Mh[α], Eq[α] を参照するが、これらに対しても同様に、初期条件としてPv[α] = 1, Mv[α] = Ph[α] = Mh[α] = Eq[α] = 0 が常に成り立つと仮定することにより、偽の類似配列が求められることを避けることができる。
【００３４】
最適アラインメントを求めるトレースバック処理を開始するためには、DP行列の下辺上のDPコスト値 C[m, j]（j = 1, 2, ... , n）の最小値を求める必要がある。類似配列に対しては、この最小値は対角線との交点の近くに現れる。そこで、最小値を探す範囲を、DP行列の下辺上で、対角線周辺領域106に含まれる範囲902に制限する。対角線周辺領域106の内部では、隣接する格子点間のDPコストの差分値を上述の方法で計算しているので、それらの値を図9内の矢印901の列で示すように、C[0, 0] = 0から出発して、先ず対角線に沿って加算し、次に、DP行列の下辺で、対角線との交点から左右両方向に加算することにより、対角線周辺領域106内に含まれるDP行列下辺902上のDPコスト値を計算し、その最小値を求めることができる。また、類似配列の最適アラインメントを求めるトレースバック処理で参照されるのは、対角線周辺領域内の隣接する格子点間のDPコストの差分値なので（非特許文献８）、上述の計算過程で各 k に対するPh, Mh, Pv, Mv の値を計算機主記憶内に保持しておき、トレースバック処理でこれらの値を参照できる。これらの計算は、配列長に比例する少ない計算コストで計算できる。
【００３５】
［実施例２］
実施例１では、DP行列の対角線周辺領域として −w < i−j≦w を満たす2重インデクス (i, j) の範囲を採用した。ここで、wは計算機のワード長を表わす。この場合、実施例１の方法で最適なアラインメントが求まるためには、挿入長、欠失長の累積値が w 未満でなければならない。従って、より多くの挿入・欠失を含むような最適なアラインメントを求める必要がある場合には、DP行列の対角線周辺領域の幅を広げる必要がある。本実施例では、DP行列の対角線周辺領域の幅を、実施例1の2倍に広げるために、実施例1の方法を修正して適用する方法を説明する。DP行列の対角線周辺領域の幅を3倍以上に広げる場合も、容易な類推により、同様な方法が可能である。
【００３６】
DP行列の対角線周辺領域の幅を、実施例1の2倍に広げると、−2 w < i−j ≦2 w を満たす2重インデクス (i, j) の範囲が計算の対象となる。そこで、この領域を対角線を挟む２つの部分領域、−2 w < i−j ≦0 と0 < i−j ≦2 w に分割する。それぞれの部分領域は、幅wであるため、その内部では実施例１と類似の方法によるビット並列なDP行列の計算が可能となる。但し、これら２つの部分領域が互いに接する（対角線に接する）位置では、互いに他の部分領域で計算されたブール変数Pv[α], Mv[α] , Ph[α], Mh[α], Eq[α] の値を境界条件として使用しなければならない。その他は、実施例1と全く同様に実施できる。
【符号の説明】
【００３７】
101 DP行列の2重インデクス [i, j] の全体を表わす正方格子
102 クエリー配列のインデクス（文字位置）
103 ターゲット配列のインデクス（文字位置）
104 DP行列の対角線：i = j
105 DP行列内の対角線周辺領域の境界線：i−j = ±w
106 DP行列内の対角線周辺領域：−w < i−j ≦w
107 奇数の k に対応する水平または垂直方向の隣接格子点を結ぶ辺
108 偶数の k に対応する水平または垂直方向の隣接格子点を結ぶ辺
109 DP行列内の隣接格子点を結ぶ辺で囲まれたセル
110 ワード・データ：Ph, Mh, Pv, Mv, Eq
111 ワード・データに対するビット並列演算
201 計算機
202 外部記憶装置
203 ゲノム配列データベース
204 DNAシーケンサ
205 DNAサンプル
206 リード配列データ（クエリー配列データ）
207 データ入力装置
208 パラメータ入力装置
209 主記憶装置
210 シード検索処理部
211 詳細比較処理部
212 最適アラインメント計算処理部
213 出力処理部
214 出力装置
215 検索結果（ゲノム・マッピング位置情報）
216 最適アラインメント結果（塩基変異の情報）
301 DP行列の格子点
302 入力側ブール変数
303 出力側ブール変数
407 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）の集合
408 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）
409 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）の集合
410 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）
411 ワード変数のインデクス
501 水平方向の隣接辺に対応するブール変数への（１重）インデクスαの対応付け
502 垂直方向の隣接辺に対応するブール変数への（１重）インデクスαの対応付け
503 対角線周辺領域内で、一定のkに対応する（1重）インデクス付けされた（水平または垂直方向の隣接辺に対応する）ブール変数の集合
504 対角線周辺領域内で、一定のkに対応する（1重）インデクス付けされた（水平方向の隣接辺に対応する）ブール変数の集合を詰め込んだワード変数
505 対角線周辺領域内で、一定のkに対応する（1重）インデクス付けされた（垂直方向の隣接辺に対応する）ブール変数の集合を詰め込んだワード変数
601 kが偶数のときの入力側ブール変数
602 k + 1が奇数のときの出力側ブール変数
603 kが奇数のときの入力側ブール変数
604 k + 1が偶数のときの出力側ブール変数
901 DPコスト値を得るために累積されるDPコストの差分値
902 DP行列の下辺上で、対角線周辺領域106に含まれる範囲（DPコスト値の最小値を探す範囲）
【技術分野】
【０００１】
本発明は、文字列、DNA塩基配列等の文字列データ（文字列データとして取り扱い可能な形式に加工された任意のデータを含む。）を、検索文字列と比較照合する処理を高速化するための技術に関する。
【背景技術】
【０００２】
近年、従来型のキャピラリー型ＤＮＡシーケンサとは異なる原理に基づき、一回のランで数千万本以上にも及ぶ大量の配列を読み取ることができる超並列ＤＮＡシーケンサが出現し、シーケンシングに要するコストを劇的に低減させた（非特許文献１）。そのため、大量の配列データを高速に解析する技術に対する需要が高まりつつある。中でも、ゲノム・マッピング解析は、最も基本的な配列データ解析の一つであり、シーケンサにより新たに読まれた塩基配列（リード）を、データベース中の参照ゲノム配列と比較して、ゲノム内の位置を特定し、さらに、そこに含まれる変異（塩基の置換・挿入・欠失）を同定する。
【０００３】
初期の超並列ＤＮＡシーケンサのリード長は30塩基程度と比較的短かったため、ハッシュ法やBurrows-Wheeler変換などを用いて高速・高精度なゲノム・マッピング解析がなされた（非特許文献２）。
【０００４】
しかし、リード長が100塩基以上まで伸びてくると、従来型のＤＮＡシーケンサ向けのゲノム・マッピング解析で用いられてきたseed-and-extend strategyのような方法が必要となってくる（非特許文献３）。この方法では、第一ステップで、リードに含まれる短い配列（シード）がゲノム内に出現する位置を高速に検索し、第二ステップで、シードの出現位置周辺でリード全体がゲノム配列に類似しているか（即ち、シードとゲノムの対応関係を、多少の変異を許容しながら、リード全体にまで拡張できるか）を調べる。この第一ステップでは、シード長が短いため、ハッシュ法やBurrows-Wheeler変換などを用いた高速な検索方法が利用でき、これにより、長大なゲノム配列全体の中からリードのマッピング候補位置を高速に絞り込むことができる。第二ステップでは、第一ステップで絞り込まれた候補位置にあるゲノム部分配列（ターゲット配列）に対して、動的計画法（DP法）などを用いてリード（クエリー配列）とのアラインメントを行い、リード長が長い場合であっても、塩基変異を高精度に調べることができる（非特許文献８）。
【０００５】
しかし、一般に、動的計画法では、配列長の2乗に比例する大きな計算コストが必要となる。そこで、動的計画法によるアラインメント処理を高速化することは、ゲノム・マッピング解析処理全体を高速化する上で有益となる。
【０００６】
動的計画法によるアラインメント処理を高速化するために、複数のプロセッサを利用して処理を並列化する方法が知られている（非特許文献４）。このとき、プロセッサ間で受け渡されるデータの依存関係を考慮し、また、プロセッサ間通信コストを抑えるような負荷分散方式を採用する必要がある。それにより、多数のプロセッサを使えば、高い並列性が得られる。
【０００７】
ところで、現在主流のデジタル計算機では、8バイト整数（ワード）に対する演算が基本演算となっているが、これは、8バイト即ち64ビットの並列演算と見ることもできる。G. Myersは、この整数演算のビット並列性を利用して、動的計画法によるアラインメント処理を１プロセッサで並列化し、配列長に比例する時間で計算する方法を示した（非特許文献５）。
【０００８】
G. Myersの方法では、長さmのクエリー配列（リード）と長さnのターゲット配列（ゲノム部分配列）のアラインメントを求めるために、サイズm×nのDP行列の情報を長さmの列方向の情報のn個の並びに分解して表現し、新たに読み込んだターゲット配列の１文字とクエリー配列の全m文字との比較に関わる処理を、mビット並列に行い、行方向（ターゲット配列の方向）に1塩基ずつ移動しながらこの処理を繰り返す。
【０００９】
式を用いて詳細に説明すると、次のようになる。１塩基の置換・挿入・欠失のコストを全て１として、DP行列の第 (i, j) 成分C[i, j] （クエリー配列 q = q1 q2 ... qm の長さ iの部分文字列q1 q2 ... qi とターゲット配列 t = t1 t2 ... tn の長さ jの部分文字列 t1 t2 ... tjとの最適なアラインメントのコスト）に対して、
等式C[i, j]−C[i, j−1] = 1が成立するか否かを表わすブール変数をPh[i, j]、
等式C[i, j]−C[i, j−1] = −1が成立するか否かを表わすブール変数をMh[i, j]、
等式C[i, j]−C[i−1, j] = 1が成立するか否かを表わすブール変数をPv[i, j]、
等式C[i, j]−C[i−1, j] = −1が成立するか否かを表わすブール変数をMv[i, j]、
クエリー配列のインデクス iの文字（i番目の文字）qiとターゲット配列のインデクスjの文字tjとが一致するか否かを表わすブール変数をEq(i, j)、
と表わし、これらの 1ビット（偽または真、即ち、0または1）で表わされるブール変数をi = 1, 2, ..., m について下位ビットより昇順に並べてできるワード変数を、
（数１)
Ph(j) = [ Ph[m, j], ・・・ , Ph[2, j], Ph[1, j] ]
Mh(j) = [ Mh[m, j], ・・・ , Mh[2, j], Mh[1, j] ]
Pv(j) = [ Pv[m, j], ・・・ , Pv[2, j], Pv[1, j] ]
Mv(j) = [ Mv[m, j], ・・・ , Mv[2, j], Mv[1, j] ]
Eq(j) = [ Eq[m, j], ・・・ , Eq[2, j], Eq[1, j] ]
と表わし、さらに記号の簡略化のためにjを一つの値に固定して、Ph(j), Mh(j), Pv(j), Mv(j), Eq(j) の値をそれぞれ、Ph, Mh, Pv, Mv, Eqと略記し、Xv、Xhを補助ワード変数として、
（数２)
Xh = ((Eq & Pv) + Pv) ^ Pv | Eq
Ph = Mv | ~ (Xh | Pv)
Mh = Pv & Xh
Ph <<= 1
Mh <<= 1
Xv = Eq | Mv
Pv = Mh | ~(Xv | Ph)
Mv = Ph & Xv
によりビット並列な整数演算を行い、Ph, Mh, Pv, Mvの値を上書きする。ここで、ワード変数a、bに対して、a|bはビットごとのOR演算、a&bはビットごとのAND演算、a ^ bはビットごとの排他的論理和演算、a + bは整数としての加算、~aはビットごとのNOT演算、a<<=1は1ビットの左シフト演算を表わす。これにより、Ph, Mh, Pv, Mvの値は、Ph(j+1), Mh(j+1), Pv(j+1), Mv(j+1) の値に更新される。これが、アラインメント処理をビット並列化する中核部分であり、これにより処理が高速化される。
【００１０】
しかし、この方法を直接的に用いることができるのは、リード長（クエリー配列長）m が計算機のワード長（64）以下の場合に限られる。この方法をリード長が計算機のワード長を越える場合に拡張するためには、リード（クエリー配列）をワード長以下の長さの複数の断片配列に分割する必要がある（非特許文献５、非特許文献６）。このとき、これらの断片配列に対する計算には互いに依存関係があり、並列化する際に制約が生じ、処理効率も低下する。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】Rajendrani Mukhopadhyay. DNA sequencers: the next generation. Analytical Chemistry, Vol. 81, No. 5, March 1, 2009
【非特許文献２】Kouichi Kimura, Yutaka Suzuki, Sumio Sugano, Asako Koike. Computation of Rank and Select Functions on Hierarchical Binary String and Its Application to Genome Mapping Problems for Short-Read DNA Sequences. Journal of Computational Biology. November 2009, 16(11): 1601-1613.
【非特許文献３】Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005 May 1;21(9):1859-75.
【非特許文献４】Martins WS, Del Cuvillo JB, Useche FJ, Theobald KB, Gao GR. A multithreaded parallel implementation of a dynamic programming algorithm for sequence comparison. Pac Symp Biocomput. 2001:311-22.
【非特許文献５】G. Myers. A fast bit-vector algorithm for approximate string matching based on dynamic programming. Journal of the ACM, 46(3):395-415, 1999.
【非特許文献６】HEIKKI HYYRO. BIT-PARALLEL COMPUTATION OF LOCAL SIMILARITY SCORE MATRICES WITH UNITARY WEIGHTS. International Journal of Foundations of Computer Science, 17(6): 1325-1344, 2006.
【非特許文献７】Kevin Judd McKernan, Heather E. Peckham, Gina L. Costa, et al. Sequence and structural variation in a human genome uncovered by short-read, massively parallel ligation sequencing using two-base encoding. Genome Res. 2009 19: 1527-1541.
【非特許文献８】Richard Durbin, Sean R. Eddy, Andrew Krogh, Graeme Mitchison, Biological Seuqnece Analysis, probabilistic models of proteins and nucleic acids, Cambridge University Press 1988（邦訳、阿久津達也、浅井潔、矢田哲士訳、バイオインフォマティクス―確率モデルによる遺伝子配列解析―、医学出版、2001）
【発明の概要】
【発明が解決しようとする課題】
【００１２】
本発明が解決しようとする課題は、リード長に対する制限を設けることなく、また、リード配列を分割するなどの処理効率を低下させる方法もとらず、計算機のワード演算のビット並列性を利用して、動的計画法によるアラインメント処理を1プロセッサで並列化し高速化する方法を提供することにある。
【課題を解決するための手段】
【００１３】
本発明による方法では、長さmのクエリー配列（リード）と長さnのターゲット配列（ゲノム部分配列）のアラインメントを求めるために、サイズm×nのDP行列の計算を対角線（行のインデクスiと列のインデクスjが等しいところ）の周辺領域に制限し、その領域内の情報を対角線に直行する方向（反対角線方向）の情報のビットの並びに分解して表現し、計算機のワード演算のビット並列性を利用して、反対角線上でのクエリー配列とターゲット配列の文字比較に関わる処理を1プロセッサで並列化し、対角線方向に移動しながらこの処理を繰り返す。
【００１４】
式を用いて詳細に説明すると、次のようになる。C[i, j], Ph[i, j], Mh[i, j], Pv[i, j], Mv[i, j], Eq[i, j] を上記の通りとして、これらを対角線周辺で反対角線方向に並べてできるワード変数を、
（数３）
Ph(k) = [ Ph[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Mh(k) = [ Mh[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Pv(k) = [ Pv[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Mv(k) = [ Mv[i, j] | i + j = k, −m ＜ i−j ≦ m ]
Eq(k) = [ Eq[i, j] | i + j = k, −m ＜ i−j ≦ m ]
とする。ここで、ビット変数をワード変数に詰め込む順番は、jについて下位ビットから昇順とする。さらに、記号の簡略化のためにkを一つの値に固定して、Ph(k), Mh(k), Pv(k), Mv(k), Eq(k) の値をそれぞれ、Ph, Mh, Pv, Mv, Eqと略記し、Xv、Xhを補助ワード変数として、
kが偶数のときは、
（数４）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
kが奇数のときは、
（数５）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' << = 1; Mv' <<= 1;
に従ってビット並列な整数演算を行い、Ph', Mh', Pv,' Mv'の値を計算する。これらは、Ph(k+1), Mh(k+1), Pv(k+1), Mv(k+1) の値を与える。これが、アラインメント処理をビット並列化する中核部分である。
【発明の効果】
【００１５】
G. Myersの方法と同様、アラインメント処理が１プロセッサでビット並列化され高速化される。G. Myersの方法と異なるのは、動的計画法のビット並列計算を、DP行列の対角線周辺領域に限定して対角線方向に移動させるため、リード長に関する制限が無くなり、従って、リード長がワード長を越えた場合にリードを分割する必要がなくなる点である。また、類似性の高い配列の対する動的計画法の計算は、DP行列の対角線の周辺を行えば十分であるため、対角線から離れた領域での不要な計算を避けることができ、処理効率が向上する。また、行と列に関して対称な形で計算を行えるような新規なビット並列化方式であるため、（数４）及び（数５）において補助変数Xhの計算式はXvの計算式と対称な形に簡易化され、これはG. Myersによる（数２）における補助変数Xhのやや複雑な計算式よりも少ないマシン命令数で計算可能となる。
【図面の簡単な説明】
【００１６】
【図１】本発明による最適アラインメント計算方法において、ビット並列演算の適用箇所を示した説明図である（実施例１）。
【図２】次世代型DNAシーケンサの配列データと参照ゲノム配列データベースとの高速比較照合を行うシステムの構成例を示す図である（実施例１）。
【図３】DP行列内のセル周辺における入力側と出力側のブール変数の位置関係を示す説明図である。
【図４】和が k となる2重インデクス [i, j] とワード変数のインデクスとの対応関係を例示する説明図である。
【図５】一定の k に対するブール変数Ph, Mh, Pv, Mv のインデクスとワード変数のインデクスの対応関係を例示する説明図である。
【図６】セル周辺のブール変数Ph, Mh, Pv, Mvのインデクスが、k が偶数か奇数かにより違いが生じることを示す説明図である。
【図７】ブール変数 Eq をビット並列に計算する方法を示す説明図である。
【図８】対角線周辺領域の外部で事前に設定されるブール変数値を示す説明図である。
【図９】DPコストの絶対値を求めるためのDPコストの差分値の累積法を示す説明図である。
【発明を実施するための形態】
【００１７】
［実施例１］
以下、本発明の実施例を図面を用いて詳細に説明する。
本実施例では、リード（クエリー配列）と参照ゲノム配列データが与えられたとき、リードに類似した部分配列が存在するゲノム配列内の位置（マッピング位置）を求め、クエリー配列とマッピング位置周辺のゲノム部分配列（ターゲット配列）との最適アラインメントを求めるための方法を説明する。通常、DNAシーケンサから得られるリード配列は塩基配列であるため、以下、リード配列は塩基配列であると仮定する。特殊な場合として、２ベース・エンコーディング方式を採用したDNAシーケンサでは、リード配列は、隣り合う２塩基を4色に符号化したカラー配列となる（非特許文献７）。このような場合は、参照ゲノム配列をカラー配列に変換して、本実施例で述べる方法に類似の方法を適用することは、容易な類推により可能である。また、アミノ酸配列をクエリー配列として、参照タンパク配列データベースを検索する課題に対しても、容易な類推により、本実施例で述べる方法に類似の方法が適用可能である。また、文字列や画像を含む文書その他のデータベースを検索する課題に対しても、容易な類推により、本実施例で述べる方法に類似の方法が適用可能である。
【００１８】
図2は、本実施例のシステム構成例を示す図である。201は、このシステムが稼働している計算機であり、外部記憶装置202に格納された参照ゲノム配列データ203の情報を主記憶209内に読み込み、またＤＮＡシーケンサ204を用いてＤＮＡサンプル205を解析して得られた塩基配列データ（クエリー配列、リード配列）206を、入力装置207を介して主記憶209内に読み込み、また、ユーザにより指定される処理パラメータ（類似配列のアラインメント・コストの上限値M）を入力装置208を介して主記憶209内に読み込む。
【００１９】
計算機201は、主記憶209 に蓄えられたリード配列と参照ゲノム配列に対して、シード検索処理210を行い、リードに含まれる短い配列（シード配列）がゲノム配列内に出現する位置を求める。シード配列としては、リード配列の中から長さ12〜20塩基程度の短い部分配列を切り出したものを使えばよく、また、シード配列の出現位置を求めるためには、Burrows-Wheeler変換など用いた公知の高速な検索方法を用いればよい（非特許文献２）。
【００２０】
詳細比較処理部211は、シード配列がゲノム配列内に出現する位置の周辺で、最適アラインメント計算212を行い、そのアラインメントのコスト値がユーザにより指定された上限値M以下ならば、そこにリード配列に類似した配列があると判断し、出力処理部213にゲノム配列内の出現位置情報とアラインメント情報を伝える。出力処理部213は、出力装置214に、検索結果としてゲノム配列内の出現位置情報（ゲノム・マッピング位置情報）215と、最適アラインメント結果216を出力する。
【００２１】
アラインメントのコストは、アラインメントに内に含まれる１塩基の置換、挿入、欠失のコストを１として、その総和として計算する。クエリー配列のp塩基目の位置から始まる長さsのシード配列が、ゲノム配列中の座標区間 [x, x + s−1] に出現している場合、リード配列の出現位置は、座標区間 [x−p+1−M, x−p+m+M] に含まれる。ここで、mはクエリー配列長、Mはアラインメント・コストの上限値を表わす。リード配列の出現位置の範囲が±Mだけ広がるのは、シード配列の出現位置の前後に、最大M塩基の挿入・欠失が含まれる可能性があるからである。最適アラインメント計算212では、この座標区間 [x−p+1−M, x−p+m+M] のゲノム部分配列をターゲット配列として、クエリー配列とターゲット配列との間で、セミ・グローバルなアラインメントを行う。即ち、クエリー配列の全体と、ターゲット配列の部分配列との間に、挿入・欠失を許しながら対応関係を定める。その計算は、具体的には、DP行列を計算することによってなされる（非特許文献８）。
【００２２】
クエリー配列長をm、ターゲット配列長をnとする。１≦ i ≦m，１≦ j ≦nに対して、クエリー配列 q = q1 q2 ... qm の長さiの部分文字列 q1 q2 ... qi とターゲット配列 t = t1 t2 ... tn の長さjの部分文字列 t1 t2 ... tj との最適なアラインメントのコストを表わすDP行列の第 (i, j) 成分を C[i, j] として、G. Myersに従い（非特許文献５）、
等式C[i, j]−C[i, j−1] = 1が成立するか否かを表わすブール変数をPh[i, j]、
等式C[i, j]−C[i, j−1] = −1が成立するか否かを表わすブール変数をMh[i, j]、
等式C[i, j]−C[i−1, j] = 1が成立するか否かを表わすブール変数をPv[i, j]、
等式C[i, j]−C[i−1, j] = −1が成立するか否かを表わすブール変数をMv[i, j]、
クエリー配列のインデクスiの文字（i番目の文字）qi とターゲット配列のインデクスjの文字（j番目の文字） tj が一致するか否かを表わすブール変数をEq(i, j)、
と表わす。また、セミ・グローバルなアラインメントに対応する境界条件として、C[i, 0] = i, C[0, j] = 0と定める。任意の１≦ i ≦m，１≦ j ≦nに対して、
（数６）
Mh = Mh[i−1, j]; Ph = Ph[i−1, j];
Mv = Mh[i, j−1]; Pv = Pv[i, j−1];
Mh' = Mh[i, j]; Ph' = Ph[i, j];
Mv' = Mv[i, j]; Pv' = Pv[i, j];
Eq = Eq[i, j];
と略記すると、これらブール変数の間には、
（数７）
Xv = Eq or Mv; Pv' = Mh or not(Xv or Ph); Mv' = Ph and Xv;
Xh = Eq or Mh; Ph' = Mv or not(Xh or Pv); Mh' = Pv and Xh;
の関係式が成立する（非特許文献５）。ここで、XvとXhは補助のブール変数を表わす。この関係式を反復して利用することによりC[i, j]の値を全て計算することができ、DP行列の下辺でC[m, j]（j = 1, 2, ... , n）の最小値を与える点 (m, j) からトレースバックすることにより最適なアラインメントを求めることができる（非特許文献８）。
【００２３】
ここで、複数組のPh, Mh, Pv, Mvから複数組のPh', Mh', Pv', Mv'を同時にビット並列に計算して、処理を高速化する方法を考える。図３に、これらのブール変数のDP行列内部における位置関係を示す。これらのブール変数は、DP行列内のセル109の周りにおいて、水平方向または垂直方向に隣接するDP行列の格子点301におけるコスト値 C[i, j] の差分に関する情報を表現している。（数６）より、入力側302のブール変数Ph, Mh, Pv, Mvの2重インデクスは [i−1, j] または [i, j−1] でその和は i + j−1 であり、また、出力側303のブール変数Ph', Mh', Pv', Mv'の2重インデクスは [i, j] でその和は i + j である。従って、2重インデクス [i, j] の和 k が小さいものから順に計算することによって、全てのブール変数Ph, Mh, Pv, Mv の値を計算することができる。さらに、クエリー長に関する制限を設けることなく、これらの計算をビット並列化できる方法があることを以下に説明する。
【００２４】
ところで、和が k に等しい2重インデクス [i, j] の個数は、 k と共に増大するが、類似配列の最適アラインメントの計算でトレースバックする際に参照すべきブール変数の2重インデクスは対角線 i = j の周辺の範囲に限られる。そこで、この範囲に制限し、（数７）の関係式をビット並列に計算する方法を考える。
【００２５】
計算機の基本演算に用いられる整数ワードのビット長をwとする。wは一般に２のべき乗であり、現在主流の計算機では、w = 64 となっている。また、半ワード長をw' = w／2とする。対角線の周辺の範囲として、−w < i−j≦w の範囲を考えると、和が k となる2重インデクス [i, j] は、
k が偶数 k = 2k' の場合、
（数８）
[i, j] = [k' + w', k'−w'], ..., [k'+1, k'−1], [k', k'], [k'−1, k'+1],
..., [k'−w' + 1, k' + w'−1]
k が奇数 k = 2k' + 1 の場合、
（数９）
[i, j] = [k' + w', k'−w' + 1], ..., [k'+1, k'], [k', k'+1],
..., [k'−w' + 1, k' + w']
となる。そこで、これらの2重インデクス（格子点）に、（1重）インデクスα= k'−i を対応させれば、kが偶数であるか奇数であるかに関わらず、これらは、
α = −w', −w' + 1, ..., −1, 0, 1, ..., w'−1
に対応する。そこで、wビットのワード変数の各ビットを、下位ビットより順にこれらでインデクス付けして対応させる。以下、説明を簡単にするために w = 4 として、図を用いてこれらの対応関係を具体的に説明する。
【００２６】
図４は、w = 4 の場合に、対角線の周辺範囲 −w < i−j ≦w にある和が k となる2重インデクス [i, j] と4ビットのワード変数のインデクスとの対応関係を例示したものである。101は、DP行列の2重インデクス [i, j] の全体を表わす正方格子である。102と103に表示したインデクスi、jは、各々、クエリー配列とターゲット配列のインデクス（文字位置）である。104で示した対角線 i = j の周辺領域−w < i−j ≦wは、105に示す2本の破線に挟まれた範囲106である。この範囲106の中で、2重インデクスの和がk = 5（奇数）となるものは、407に示す破線で囲まれた中にある4個の白丸（408）の位置にある。同様に、この範囲106の中で、2重インデクスの和がk = 6（偶数）となるものは、409に示す破線で囲まれた中にある4個の黒丸（410）の位置にある。k = 7（奇数）の場合も同様である。これらに対応するインデクスα = −2, −1, 0, 1 の値を傍らに付した。また、これらのインデクスに対応する4ビット整数のビット位置を411に示した。
【００２７】
和が k となる2重インデクス [i, j] をもつブール変数Ph, Mh, Pv, Mv に対しても、同様にワード変数の（1重）インデクスαを対応させる。ブール変数Ph, Mhは水平方向に隣接する格子点のコストの差分値の情報を表わしているが、それらには左方の格子点のインデクスαを対応させる（図5の501参照）。即ち、Ph[i, j], Mh[i, j] には [i, j−1] に対応するインデクスαを対応させて、Ph[α], Mh[α] と表わす。同様に、ブール変数Pv, Mvは垂直方向に隣接する格子点のコストの差分値の情報を表わしているが、それらには下方の格子点のインデクスαを対応させる（図5の502参照）。即ち、Pv[i, j], Mv[i, j] には [i, j] に対応するインデクスαを対応させて、Pv[α], Mv[α] と表わす。図5に、w = 4 の場合の例を具体的に示す。h[α]はブール変数Ph[α]またはMh[α]を表わし、また、v[α]はブール変数Pv[α]またはMv[α]を表わす。
【００２８】
対角線周辺領域106内に、一定の k に対する Ph[α], Mh[α], Pv[α], Mv[α] は各々w個ずつあるので、対応するインデクスαに従いそれらのビット値503を1ワードに詰め込むことができ、そのようにして得られるワード変数504, 505をPh, Mh, Pv, Mv と表わす。図5はw = 4 の場合を示す。ここで、hはPhまたはMhを、vはPvまたはMvを表わす。ブール変数Eq = Eq[i, j]についても、同様に、2重インデクス [i, j] に対応するインデクスαを対応させて Eq[α] と表わし、それに従ってビット値をワード変数に詰め込んだものをEqと書く。
【００２９】
これらのワード変数は、 k について昇順に計算することができる。図1は、w = 4 の場合に、その計算過程におけるワード変数のデータの依存関係と、その時に計算されるブール変数のDP行列内における位置を示したものである。対角線周辺領域106（−w < i−j ≦w）の中で、kが奇数のとき、ブール変数Ph[α], Mh[α], Pv[α], Mv[α]に対応する水平または垂直方向の隣接格子点を結ぶ辺を太い実線107で、また、kが偶数のとき、対応する同様の辺を太い点線108で表わした。また、辺で囲まれた各セル109には、ブール変数Eq[α]が対応している。これらのブール変数をワード変数に詰め込んだものがPh, Mh, Pv, Mv, Eqであり、これらは110に示すワード・データを構成する。これらは、以下に説明するビット並列演算111によりkについて昇順に計算される。
【００３０】
ビット並列演算111を説明するために、図5に戻って、各セル109の周りで計算すべき式を確認する。但し、セルの左辺と上辺に対応するkの値が偶数か奇数かによって事情が異なる。先ず、kが偶数の場合は、図6左に示すように、kに対応する左辺と上辺の入力側ブール変数601はPv[α], Mv[α] , Ph[α], Mh[α] であり、また、k+1に対応する下辺と右辺の出力側ブール変数602はPh[α], Mh[α], Pv[α], Mv[α]であり、これらは全て同じインデクスαをもっている。従って、この場合は、ワード変数 Ph, Mh, Pv, Mv に対して、kに対応するものをPh, Mh, Pv, Mv、また、k+1に対応するものをPh', Mh', Pv', Mv'と書いて区別すれば、（数７）と同様の関係式、
（数４）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
により、ビット並列に計算できる。ここで、ワード変数a、bに対して、a|bはビットごとのOR演算、a&bはビットごとのAND演算、~aはビットごとのNOT演算を表わす。一方、kが奇数の場合は、図6右に示すように、kに対応する左辺と上辺の入力側ブール変数603はPv[α], Mv[α] , Ph[α+1], Mh[α+1] であり、また、k+1に対応する下辺と右辺の出力側ブール変数604はPh[α], Mh[α], Pv[α+1], Mv[α+1] であり、インデクスαに+1のずれが生じている。そこで、この場合は、上記のワード変数 Ph, Mh, Pv, Mv 及び Ph', Mh', Pv', Mv' に対して、前後にビットシフト演算を行えば、（数７）と同様な関係式が適用可能となる。即ち、この場合は、
（数５）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' <<= 1; Mv' <<= 1;
により、ビット並列に計算できる。ここで、ワード変数aに対して、a<<=1 は1ビットの左シフト演算、a>>=1 は1ビットの右シフト演算を表わす。
【００３１】
一方、ワード変数 Eq については、
k が偶数 k = 2k' の場合、
（数10）
Eq(2k') = [ Eq[−w'], ..., Eq[−1], Eq[0], Eq[1], ..., Eq[w'−1 ] ]
= [ Eq[k'−w' + 1, k' + w'−1], ...,
Eq[k'−1, k'+1], Eq[k', k'], Eq[k'+1, k'−1],
..., Eq[k' + w', k'−w'] ]
k が奇数 k = 2k' + 1 の場合、
（数11）
Eq(2k'+1) = [ Eq[−w'], ..., Eq[−1], Eq[0], Eq[1], ..., Eq[w'−1 ] ]
= [ Eq[k'−w' + 1, k' + w'], ...,
Eq[k', k'+1], Eq[k'+1, k'],
..., Eq[k' + w', k'−w' + 1] ]
となる。ここで、外側の[ , , ... , ] はビットを並べてできる整数値を表わし（右端が最下位ビット）、Eq[i, j] は、クエリー配列 q = q1 q2 ... qm のi番目の文字（塩基）qiとターゲット配列 t = t1 t2 ... tn のj番目の文字（塩基）tjが一致するか否かを表わすブール変数を表わす。4種類の塩基A, C, G, Tの2ビット符号化 Aa00, Ca01, Ga10, Ta11 における下位ビットをl、上位ビットをu と表わし、l(A) = l(G) = 0, l(G) = l(T) = 1, u(A) = u(C) = 0, u(G) = u(T) = 1 とすると、
（数12）
Eq[i, j] = not( ( l(qi) xor l(tj) ) or ( u(qi) xor u(tj) ))
が成り立つ。ここで、xorはブール変数の排他的論理和を表わす。そこで、ブール変数l(tj)、u(tj)を最下位ビットから昇順に並べてできるワード変数をTl, Tu、ブール変数l(qi)、u(qi)を最上位ビットから降順に並べてできるワード変数をQl, Quとすれば、
ワード変数Eqは、
（数13）
Eq = ~( ( Ql ^ Tl ) | ( Qu ^ Tu ))
のようにビット並列に計算できる。ここで、ワード変数a, bに対して、a ^ bはビットごとの排他的論理和演算を表わす。一方、Tl, Tu, Ql, Quは、k について帰納的にビットシフト演算を用いて計算できる。但し、ビットシフトによって空いたビットには、クエリー配列またはターゲット配列の対応する文字の符号化ビットをセットする。図7にw = 4 の場合のk = 4, 5, 6に対応するTl, Tu, Ql, Quの計算方法を示す。一般の場合に、kが増加するときのTl, Tu, Ql, Quの値の更新方法を具体的に式で表わすと、
kが奇数（k = 2k' + 1）のとき、
（数14）
Ql <<= 1; Ql |= l(qi);
Qu <<= 1; Qu |= u(qi);
kが偶数（k = 2k'）のとき、
（数15）
Tl >>= 1; Tl |= (l(tj) << (w−1));
Tu >>= 1; Tu |= (u(tj) << (w−1));
となる。ここで、qi, tjはクエリー配列とターゲット配列から新たに読み込まれる文字を表わし、i = w' + k' + 1、j = w' + k' である。また、ワード変数a, b に対して、a = a | b を a |= b と略記する。
【００３２】
上述の計算において、対角線周辺領域106の境界105 （ i−j = ±w ）に内部から接する辺に対応するブール変数Pv[α], Mv[α] , Ph[α], Mh[α] を計算する際、この境界に外部から接する辺に対応するブール変数Pv[α], Mv[α] , Ph[α], Mh[α] の値を参照する必要が生じる。後者の値は計算されないので、境界条件として、それらに適切な値を事前に定めておく必要がある。そこで、対角線周辺領域106の外部では、クエリー配列とターゲット配列の文字比較結果は常に不一致と仮定する。これは、偽の類似配列が求められることを避けるための（悲観的）仮定である。このことは、対角線周辺領域106の外部では C[i, j]−C[i−1, j−1] = 1 が常に成り立つことを意味する。そのためには、図8に示すように、対角線周辺領域の外部では常に Pv[α] = 1, Mv[α] = Ph[α] = Mh[α] = Eq[α] = 0 が成り立つと仮定すればよい。そこで、これを境界条件とする。このとき、（数７）が成立することは、これらの値を代入して直接確かめることができる。
【００３３】
また、k の値がワード長 w 以下のとき、ワード変数Ph, Mh, Pv, Mv は、DP行列外部の i < 0，j < 0 となる領域に対応するブール変数Pv[α], Mv[α] , Ph[α], Mh[α], Eq[α] を参照するが、これらに対しても同様に、初期条件としてPv[α] = 1, Mv[α] = Ph[α] = Mh[α] = Eq[α] = 0 が常に成り立つと仮定することにより、偽の類似配列が求められることを避けることができる。
【００３４】
最適アラインメントを求めるトレースバック処理を開始するためには、DP行列の下辺上のDPコスト値 C[m, j]（j = 1, 2, ... , n）の最小値を求める必要がある。類似配列に対しては、この最小値は対角線との交点の近くに現れる。そこで、最小値を探す範囲を、DP行列の下辺上で、対角線周辺領域106に含まれる範囲902に制限する。対角線周辺領域106の内部では、隣接する格子点間のDPコストの差分値を上述の方法で計算しているので、それらの値を図9内の矢印901の列で示すように、C[0, 0] = 0から出発して、先ず対角線に沿って加算し、次に、DP行列の下辺で、対角線との交点から左右両方向に加算することにより、対角線周辺領域106内に含まれるDP行列下辺902上のDPコスト値を計算し、その最小値を求めることができる。また、類似配列の最適アラインメントを求めるトレースバック処理で参照されるのは、対角線周辺領域内の隣接する格子点間のDPコストの差分値なので（非特許文献８）、上述の計算過程で各 k に対するPh, Mh, Pv, Mv の値を計算機主記憶内に保持しておき、トレースバック処理でこれらの値を参照できる。これらの計算は、配列長に比例する少ない計算コストで計算できる。
【００３５】
［実施例２］
実施例１では、DP行列の対角線周辺領域として −w < i−j≦w を満たす2重インデクス (i, j) の範囲を採用した。ここで、wは計算機のワード長を表わす。この場合、実施例１の方法で最適なアラインメントが求まるためには、挿入長、欠失長の累積値が w 未満でなければならない。従って、より多くの挿入・欠失を含むような最適なアラインメントを求める必要がある場合には、DP行列の対角線周辺領域の幅を広げる必要がある。本実施例では、DP行列の対角線周辺領域の幅を、実施例1の2倍に広げるために、実施例1の方法を修正して適用する方法を説明する。DP行列の対角線周辺領域の幅を3倍以上に広げる場合も、容易な類推により、同様な方法が可能である。
【００３６】
DP行列の対角線周辺領域の幅を、実施例1の2倍に広げると、−2 w < i−j ≦2 w を満たす2重インデクス (i, j) の範囲が計算の対象となる。そこで、この領域を対角線を挟む２つの部分領域、−2 w < i−j ≦0 と0 < i−j ≦2 w に分割する。それぞれの部分領域は、幅wであるため、その内部では実施例１と類似の方法によるビット並列なDP行列の計算が可能となる。但し、これら２つの部分領域が互いに接する（対角線に接する）位置では、互いに他の部分領域で計算されたブール変数Pv[α], Mv[α] , Ph[α], Mh[α], Eq[α] の値を境界条件として使用しなければならない。その他は、実施例1と全く同様に実施できる。
【符号の説明】
【００３７】
101 DP行列の2重インデクス [i, j] の全体を表わす正方格子
102 クエリー配列のインデクス（文字位置）
103 ターゲット配列のインデクス（文字位置）
104 DP行列の対角線：i = j
105 DP行列内の対角線周辺領域の境界線：i−j = ±w
106 DP行列内の対角線周辺領域：−w < i−j ≦w
107 奇数の k に対応する水平または垂直方向の隣接格子点を結ぶ辺
108 偶数の k に対応する水平または垂直方向の隣接格子点を結ぶ辺
109 DP行列内の隣接格子点を結ぶ辺で囲まれたセル
110 ワード・データ：Ph, Mh, Pv, Mv, Eq
111 ワード・データに対するビット並列演算
201 計算機
202 外部記憶装置
203 ゲノム配列データベース
204 DNAシーケンサ
205 DNAサンプル
206 リード配列データ（クエリー配列データ）
207 データ入力装置
208 パラメータ入力装置
209 主記憶装置
210 シード検索処理部
211 詳細比較処理部
212 最適アラインメント計算処理部
213 出力処理部
214 出力装置
215 検索結果（ゲノム・マッピング位置情報）
216 最適アラインメント結果（塩基変異の情報）
301 DP行列の格子点
302 入力側ブール変数
303 出力側ブール変数
407 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）の集合
408 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）
409 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）の集合
410 対角線周辺領域に含まれ、和がk = 5（奇数）となる2重インデクス [ i, j ] の位置（格子点）
411 ワード変数のインデクス
501 水平方向の隣接辺に対応するブール変数への（１重）インデクスαの対応付け
502 垂直方向の隣接辺に対応するブール変数への（１重）インデクスαの対応付け
503 対角線周辺領域内で、一定のkに対応する（1重）インデクス付けされた（水平または垂直方向の隣接辺に対応する）ブール変数の集合
504 対角線周辺領域内で、一定のkに対応する（1重）インデクス付けされた（水平方向の隣接辺に対応する）ブール変数の集合を詰め込んだワード変数
505 対角線周辺領域内で、一定のkに対応する（1重）インデクス付けされた（垂直方向の隣接辺に対応する）ブール変数の集合を詰め込んだワード変数
601 kが偶数のときの入力側ブール変数
602 k + 1が奇数のときの出力側ブール変数
603 kが奇数のときの入力側ブール変数
604 k + 1が偶数のときの出力側ブール変数
901 DPコスト値を得るために累積されるDPコストの差分値
902 DP行列の下辺上で、対角線周辺領域106に含まれる範囲（DPコスト値の最小値を探す範囲）
【特許請求の範囲】
【請求項１】
文字列データベースを記憶する記憶装置と、
検索文字列を入力する入力装置と、
二つの文字列が類似していると判定されるためのアラインメント・コストの上限値を入力する入力装置と、
前記検索文字列と前記文字列データベースを比較して、両者に共有される部分配列（シード配列）の文字列データベース中の出現位置を検索するシード検索処理部と、
前記シード検索処理部により得られる前記文字列データベース内の前記シード配列の出現位置の各々に対して、前記文字列データベース内におけるその近傍の部分文字列（ターゲット配列）と前記検索文字列全体とを比較して最適なアラインメントを算出し、そのコストを前記上限値と比較することにより、その位置に前記検索文字列全体に類似した配列が含まれているか否かを判定する詳細比較処理部と、
前記詳細比較処理部の前記判定結果に基づき、前記検索文字列に類似した配列が現れる前記文字列データベースの出現位置、及び、前記最適なアラインメントを出力する出力装置を有し、
前記詳細比較処理部は、前記検索文字列と前記ターゲット配列の二つの文字列データに対する前記最適なアラインメントを算出するために、部分配列に対する最適なアラインメントを順次伸長させる処理（伸長処理）を繰り返す動的計画法に基づく処理を実行し、
前記伸長処理において、それぞれの文字列データから新たな文字を読み込んで比較を行い最適なアラインメントを伸長させる際、比較される文字のインデクス（文字列先頭から数えた文字位置）の差の絶対値が計算機のワード長以下になるものに処理を制限することにより、それらのインデクスの和が一定になる伸長処理の数をワード長以下に抑えて、計算機のワードに対する演算のビット並列性を利用して、それらの伸長処理を１プロセッサで並列に行う
ことを特徴とする最適アラインメント計算装置。
【請求項２】
検索文字列とターゲット配列を入力する入力装置と、
二つの文字列が類似していると判定されるためのアラインメント・コストの上限値を入力する入力装置と、
前記ターゲット配列と前記検索文字列全体とを比較して最適なアラインメントを算出し、そのコストを前記上限値と比較することにより、前記ターゲット配列内に前記検索文字列全体に類似した配列が含まれているか否かを判定する詳細比較処理部と、
前記詳細比較処理部の前記判定結果、及び、前記最適なアラインメントを出力する出力装置を有し、
前記詳細比較処理部は、前記検索文字列と前記ターゲット配列の二つの文字列データに対する前記最適なアラインメントを算出するために、部分配列に対する最適なアラインメントを順次伸長させる処理（伸長処理）を繰り返す動的計画法に基づく処理を実行し、
前記伸長処理において、それぞれの文字列データから新たな文字を読み込んで比較を行い最適なアラインメントを伸長させる際、比較される文字のインデクス（文字列先頭から数えた文字位置）の差の絶対値が計算機のワード長以下になるものに処理を制限することにより、それらのインデクスの和が一定になる伸長処理の数をワード長以下に抑えて、計算機のワードに対する演算のビット並列性を利用して、それらの伸長処理を１プロセッサで並列に行う
ことを特徴とする最適アラインメント計算装置。
【請求項３】
請求項１に記載の最適アラインメント計算装置において、
アラインメント内での１文字の置換、挿入、欠失のコストを１として、
前記検索文字列の先頭からインデクスiまでの部分文字列と前記近傍文字列の先頭からインデクスjまでの部分文字列との最適なアラインメントのコストをC[i, j]で表わし
等式C[i, j]−C[i, j−1] = 1が成立するか否かをブール変数をPh[i, j]で表わし、
等式C[i, j]−C[i, j−1] = −1が成立するか否かをブール変数Mh[i, j]で表わし、
等式C[i, j]−C[i−1, j] = 1が成立するか否かをブール変数Pv[i, j]で表わし、
等式C[i, j]−C[i−1, j] = −1が成立するか否かをブール変数Mv[i, j]で表わし、
前記検索文字列のインデクスiの文字と前記近傍文字列のインデクスjの文字が一致するか否かをブール変数Eq(i, j)で表わし、
計算機のワード長をwで表わし、
iとjの差の絶対値がw未満で和がkとなるiとjの組に対して、0または1の1ビットで表わされる前記ブール変数Ph[i, j]をjについて最下位ビットより昇順に並べてできるワードをPh(k)と表わし、また、同様な方法で前記ブール変数Mh[i, j]、 Pv[i, j]、 Mv[i, j]、 Eq[i, j]を並べてできるワード変数をMh(k)、Pv(k)、Mv(k)、Eq(k)と表し、
計算機によるワード変数a、bに対する演算として、ビットごとのOR演算をa|b、ビットごとのAND演算をa&b、ビットごとのNOT演算を~a、1ビットの左シフト演算を a<<=1、1ビットの右シフト演算を a>>=1と表わし、
以下、記号の簡略化のためkを一つの値に固定して、前記ワード変数 Ph(k)、Mh(k)、Pv(k)、Mv(k)、Eq(k) を、Ph、Mh、Pv、Mv、Eq と表わし、
Xv、Xhを補助ワード変数として、
kが偶数のとき、
（数１）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
に従い、
kが奇数のとき、
（数２）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' <<= 1; Mv' <<= 1;
に従い、Ph'、Mh'、Pv'、Mv' を計算してワード変数 Ph(k+1)、Mh(k+1)、Pv(k+1)、Mv(k+1) の値を得ることにより伸長処理の計算を行う
ことを特徴とする最適アラインメント計算装置。
【請求項４】
請求項２に記載の最適アラインメント計算装置において、
アラインメント内での１文字の置換、挿入、欠失のコストを１として、
前記検索文字列の先頭からインデクスiまでの部分文字列と前記近傍文字列の先頭からインデクスjまでの部分文字列との最適なアラインメントのコストをC[i, j]で表わし
等式C[i, j]−C[i, j−1] = 1が成立するか否かをブール変数をPh[i, j]で表わし、
等式C[i, j]−C[i, j−1] = −1が成立するか否かをブール変数Mh[i, j]で表わし、
等式C[i, j]−C[i−1, j] = 1が成立するか否かをブール変数Pv[i, j]で表わし、
等式C[i, j]−C[i−1, j] = −1が成立するか否かをブール変数Mv[i, j]で表わし、
前記検索文字列のインデクスiの文字と前記近傍文字列のインデクスjの文字が一致するか否かをブール変数Eq(i, j)で表わし、
計算機のワード長をwで表わし、
iとjの差の絶対値がw未満で和がkとなるiとjの組に対して、0または1の1ビットで表わされる前記ブール変数Ph[i, j]をjについて最下位ビットより昇順に並べてできるワードをPh(k)と表わし、また、同様な方法で前記ブール変数Mh[i, j]、 Pv[i, j]、 Mv[i, j]、 Eq[i, j]を並べてできるワード変数をMh(k)、Pv(k)、Mv(k)、Eq(k)と表し、
計算機によるワード変数a、bに対する演算として、ビットごとのOR演算をa|b、ビットごとのAND演算をa&b、ビットごとのNOT演算を~a、1ビットの左シフト演算を a<<=1、1ビットの右シフト演算を a>>=1と表わし、
以下、記号の簡略化のためkを一つの値に固定して、前記ワード変数 Ph(k)、Mh(k)、Pv(k)、Mv(k)、Eq(k) を、Ph、Mh、Pv、Mv、Eq と表わし、
Xv、Xhを補助ワード変数として、
kが偶数のとき、
（数１）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
に従い、
kが奇数のとき、
（数２）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' <<= 1; Mv' <<= 1;
に従い、Ph'、Mh'、Pv'、Mv' を計算してワード変数 Ph(k+1)、Mh(k+1)、Pv(k+1)、Mv(k+1) の値を得ることにより伸長処理の計算を行う
ことを特徴とする最適アラインメント計算装置。
【請求項５】
コンピュータに、
検索文字列とターゲット配列の入力を受付ける処理と、
二つの文字列が類似していると判定されるためのアラインメント・コストの上限値の入力を受付ける処理と、
前記ターゲット配列と前記検索文字列全体とを比較して最適なアラインメントを算出し、そのコストを前記上限値と比較することにより、前記ターゲット配列内に前記検索文字列全体に類似した配列が含まれているか否かを判定する詳細比較処理と、
前記詳細比較処理の前記判定結果、及び、前記最適なアラインメントを出力装置に出力する処理とを実行させるプログラムであり、
前記詳細比較処理は、前記検索文字列と前記ターゲット配列の二つの文字列データに対する前記最適なアラインメントを算出するために、部分配列に対する最適なアラインメントを順次伸長させる処理（伸長処理）を繰り返す動的計画法に基づく処理を実行し、
前記伸長処理において、それぞれの文字列データから新たな文字を読み込んで比較を行い最適なアラインメントを伸長させる際、比較される文字のインデクス（文字列先頭から数えた文字位置）の差の絶対値が計算機のワード長以下になるものに処理を制限することにより、それらのインデクスの和が一定になる伸長処理の数をワード長以下に抑えて、計算機のワードに対する演算のビット並列性を利用して、それらの伸長処理を１プロセッサで並列に行う
ことを特徴とするプログラム。
【請求項１】
文字列データベースを記憶する記憶装置と、
検索文字列を入力する入力装置と、
二つの文字列が類似していると判定されるためのアラインメント・コストの上限値を入力する入力装置と、
前記検索文字列と前記文字列データベースを比較して、両者に共有される部分配列（シード配列）の文字列データベース中の出現位置を検索するシード検索処理部と、
前記シード検索処理部により得られる前記文字列データベース内の前記シード配列の出現位置の各々に対して、前記文字列データベース内におけるその近傍の部分文字列（ターゲット配列）と前記検索文字列全体とを比較して最適なアラインメントを算出し、そのコストを前記上限値と比較することにより、その位置に前記検索文字列全体に類似した配列が含まれているか否かを判定する詳細比較処理部と、
前記詳細比較処理部の前記判定結果に基づき、前記検索文字列に類似した配列が現れる前記文字列データベースの出現位置、及び、前記最適なアラインメントを出力する出力装置を有し、
前記詳細比較処理部は、前記検索文字列と前記ターゲット配列の二つの文字列データに対する前記最適なアラインメントを算出するために、部分配列に対する最適なアラインメントを順次伸長させる処理（伸長処理）を繰り返す動的計画法に基づく処理を実行し、
前記伸長処理において、それぞれの文字列データから新たな文字を読み込んで比較を行い最適なアラインメントを伸長させる際、比較される文字のインデクス（文字列先頭から数えた文字位置）の差の絶対値が計算機のワード長以下になるものに処理を制限することにより、それらのインデクスの和が一定になる伸長処理の数をワード長以下に抑えて、計算機のワードに対する演算のビット並列性を利用して、それらの伸長処理を１プロセッサで並列に行う
ことを特徴とする最適アラインメント計算装置。
【請求項２】
検索文字列とターゲット配列を入力する入力装置と、
二つの文字列が類似していると判定されるためのアラインメント・コストの上限値を入力する入力装置と、
前記ターゲット配列と前記検索文字列全体とを比較して最適なアラインメントを算出し、そのコストを前記上限値と比較することにより、前記ターゲット配列内に前記検索文字列全体に類似した配列が含まれているか否かを判定する詳細比較処理部と、
前記詳細比較処理部の前記判定結果、及び、前記最適なアラインメントを出力する出力装置を有し、
前記詳細比較処理部は、前記検索文字列と前記ターゲット配列の二つの文字列データに対する前記最適なアラインメントを算出するために、部分配列に対する最適なアラインメントを順次伸長させる処理（伸長処理）を繰り返す動的計画法に基づく処理を実行し、
前記伸長処理において、それぞれの文字列データから新たな文字を読み込んで比較を行い最適なアラインメントを伸長させる際、比較される文字のインデクス（文字列先頭から数えた文字位置）の差の絶対値が計算機のワード長以下になるものに処理を制限することにより、それらのインデクスの和が一定になる伸長処理の数をワード長以下に抑えて、計算機のワードに対する演算のビット並列性を利用して、それらの伸長処理を１プロセッサで並列に行う
ことを特徴とする最適アラインメント計算装置。
【請求項３】
請求項１に記載の最適アラインメント計算装置において、
アラインメント内での１文字の置換、挿入、欠失のコストを１として、
前記検索文字列の先頭からインデクスiまでの部分文字列と前記近傍文字列の先頭からインデクスjまでの部分文字列との最適なアラインメントのコストをC[i, j]で表わし
等式C[i, j]−C[i, j−1] = 1が成立するか否かをブール変数をPh[i, j]で表わし、
等式C[i, j]−C[i, j−1] = −1が成立するか否かをブール変数Mh[i, j]で表わし、
等式C[i, j]−C[i−1, j] = 1が成立するか否かをブール変数Pv[i, j]で表わし、
等式C[i, j]−C[i−1, j] = −1が成立するか否かをブール変数Mv[i, j]で表わし、
前記検索文字列のインデクスiの文字と前記近傍文字列のインデクスjの文字が一致するか否かをブール変数Eq(i, j)で表わし、
計算機のワード長をwで表わし、
iとjの差の絶対値がw未満で和がkとなるiとjの組に対して、0または1の1ビットで表わされる前記ブール変数Ph[i, j]をjについて最下位ビットより昇順に並べてできるワードをPh(k)と表わし、また、同様な方法で前記ブール変数Mh[i, j]、 Pv[i, j]、 Mv[i, j]、 Eq[i, j]を並べてできるワード変数をMh(k)、Pv(k)、Mv(k)、Eq(k)と表し、
計算機によるワード変数a、bに対する演算として、ビットごとのOR演算をa|b、ビットごとのAND演算をa&b、ビットごとのNOT演算を~a、1ビットの左シフト演算を a<<=1、1ビットの右シフト演算を a>>=1と表わし、
以下、記号の簡略化のためkを一つの値に固定して、前記ワード変数 Ph(k)、Mh(k)、Pv(k)、Mv(k)、Eq(k) を、Ph、Mh、Pv、Mv、Eq と表わし、
Xv、Xhを補助ワード変数として、
kが偶数のとき、
（数１）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
に従い、
kが奇数のとき、
（数２）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' <<= 1; Mv' <<= 1;
に従い、Ph'、Mh'、Pv'、Mv' を計算してワード変数 Ph(k+1)、Mh(k+1)、Pv(k+1)、Mv(k+1) の値を得ることにより伸長処理の計算を行う
ことを特徴とする最適アラインメント計算装置。
【請求項４】
請求項２に記載の最適アラインメント計算装置において、
アラインメント内での１文字の置換、挿入、欠失のコストを１として、
前記検索文字列の先頭からインデクスiまでの部分文字列と前記近傍文字列の先頭からインデクスjまでの部分文字列との最適なアラインメントのコストをC[i, j]で表わし
等式C[i, j]−C[i, j−1] = 1が成立するか否かをブール変数をPh[i, j]で表わし、
等式C[i, j]−C[i, j−1] = −1が成立するか否かをブール変数Mh[i, j]で表わし、
等式C[i, j]−C[i−1, j] = 1が成立するか否かをブール変数Pv[i, j]で表わし、
等式C[i, j]−C[i−1, j] = −1が成立するか否かをブール変数Mv[i, j]で表わし、
前記検索文字列のインデクスiの文字と前記近傍文字列のインデクスjの文字が一致するか否かをブール変数Eq(i, j)で表わし、
計算機のワード長をwで表わし、
iとjの差の絶対値がw未満で和がkとなるiとjの組に対して、0または1の1ビットで表わされる前記ブール変数Ph[i, j]をjについて最下位ビットより昇順に並べてできるワードをPh(k)と表わし、また、同様な方法で前記ブール変数Mh[i, j]、 Pv[i, j]、 Mv[i, j]、 Eq[i, j]を並べてできるワード変数をMh(k)、Pv(k)、Mv(k)、Eq(k)と表し、
計算機によるワード変数a、bに対する演算として、ビットごとのOR演算をa|b、ビットごとのAND演算をa&b、ビットごとのNOT演算を~a、1ビットの左シフト演算を a<<=1、1ビットの右シフト演算を a>>=1と表わし、
以下、記号の簡略化のためkを一つの値に固定して、前記ワード変数 Ph(k)、Mh(k)、Pv(k)、Mv(k)、Eq(k) を、Ph、Mh、Pv、Mv、Eq と表わし、
Xv、Xhを補助ワード変数として、
kが偶数のとき、
（数１）
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
に従い、
kが奇数のとき、
（数２）
Ph >>= 1; Mh >>= 1;
Xv = Eq | Mv; Pv' = Mh | ~(Xv | Ph); Mv' = Ph & Xv;
Xh = Eq | Mh; Ph' = Mv | ~(Xh | Pv); Mh' = Pv & Xh;
Pv' <<= 1; Mv' <<= 1;
に従い、Ph'、Mh'、Pv'、Mv' を計算してワード変数 Ph(k+1)、Mh(k+1)、Pv(k+1)、Mv(k+1) の値を得ることにより伸長処理の計算を行う
ことを特徴とする最適アラインメント計算装置。
【請求項５】
コンピュータに、
検索文字列とターゲット配列の入力を受付ける処理と、
二つの文字列が類似していると判定されるためのアラインメント・コストの上限値の入力を受付ける処理と、
前記ターゲット配列と前記検索文字列全体とを比較して最適なアラインメントを算出し、そのコストを前記上限値と比較することにより、前記ターゲット配列内に前記検索文字列全体に類似した配列が含まれているか否かを判定する詳細比較処理と、
前記詳細比較処理の前記判定結果、及び、前記最適なアラインメントを出力装置に出力する処理とを実行させるプログラムであり、
前記詳細比較処理は、前記検索文字列と前記ターゲット配列の二つの文字列データに対する前記最適なアラインメントを算出するために、部分配列に対する最適なアラインメントを順次伸長させる処理（伸長処理）を繰り返す動的計画法に基づく処理を実行し、
前記伸長処理において、それぞれの文字列データから新たな文字を読み込んで比較を行い最適なアラインメントを伸長させる際、比較される文字のインデクス（文字列先頭から数えた文字位置）の差の絶対値が計算機のワード長以下になるものに処理を制限することにより、それらのインデクスの和が一定になる伸長処理の数をワード長以下に抑えて、計算機のワードに対する演算のビット並列性を利用して、それらの伸長処理を１プロセッサで並列に行う
ことを特徴とするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１２−３２９７５（Ｐ２０１２−３２９７５Ａ）
【公開日】平成２４年２月１６日（２０１２．２．１６）
【国際特許分類】
【出願番号】特願２０１０−１７１４３２（Ｐ２０１０−１７１４３２）
【出願日】平成２２年７月３０日（２０１０．７．３０）
【出願人】（０００００５１０８）株式会社日立製作所 (27,607)
【Ｆターム（参考）】
【公開日】平成２４年２月１６日（２０１２．２．１６）
【国際特許分類】
【出願日】平成２２年７月３０日（２０１０．７．３０）
【出願人】（０００００５１０８）株式会社日立製作所 (27,607)
【Ｆターム（参考）】
[ Back to top ]