説明

配列解析装置、配列解析方法およびコンピュータプログラム

【課題】大量のショートリードの中から、高速に、編集距離が所定の範囲内にあるショートリードのペアを探索することができる配列解析装置を提供する。
【解決手段】配列解析装置10は、複数の塩基配列を入力する配列入力部11と、編集距離dとブロック分割数を決定する数値nまたは分割数(d+n)を入力する条件入力部12と、複数の塩基配列の各々を(d+n)個に分割してブロックを生成し、その中から選んだn個のすべてのブロックについて、読み出した部分配列、または、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致する条件を満たす塩基配列の集合を等価クラスとする等価クラス生成部13と、等価クラス生成部から、塩基配列のペアについて、編集距離を計算し、編集距離d以内であると計算されたペアを示すデータを出力する類似判定部14とを備えた。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、配列解析装置および配列解析方法に関し、特に、大量の塩基配列の中から、距離が所定範囲内にある塩基配列のペアを探索する配列解析装置および配列解析方法に関する。
【背景技術】
【0002】
DNA塩基配列の分析技術の発展により、ギガシーケンサと呼ばれるシーケンサを用いることで、大量のリード長の短い塩基配列(以下、「ショートリード」(short read)、または「SR」という)が得られる。また、昨今では、ショートリードを用いても有効な分析を行うことができる技術が提案されている。ゲノムワイドアライメント(genome-wide alignment)やデノボゲノムアセンブリ(de novo genome assembly)のような塩基配列分析では、まず、大量のショートリードの中から類似のペアを見つけ出すことが重要な作業となる。
【0003】
大量の(例えば、数千万の)ショートリードの中から類似のペアを見つけ出すために、すべてのショートリードのペアの距離を求めるとすると、膨大な量の計算を行わなければならなくなる。シーケンサからショートリードがn個出力された場合には、すべてのショートリードのペアの距離を求めるとすると、{n×(n−1)}/2回の距離計算を行わなければならない。
【0004】
距離が所定の範囲内にあるショートリードのペアを取りこぼすことなく、実際に距離計算をするショートリードのペアの数を削減する方法として、複合ソート法がある。図21は、複合ソート法を説明するための図である。複合ソート法では、複数のショートリードを対応する部分配列ごとにブロックとして分割し、その中から選んだ所定数のブロックにおいて、すべてのショートリードの部分配列を比較して、配列パターンが一致するショートリードによって等価クラス(equivalence class)を形成し、等価クラスに含まれるショートリードのペアについて距離計算を行う。この方法によれば、実際に距離を計算するショートリードのペアを削減でき、互いに類似するショートリードのペアを高速に探索できる。
【0005】
本出願に関連する先行技術文献として、以下の文献がある。
【先行技術文献】
【特許文献】
【0006】
【特許文献1】特開2009−116559号公報
【発明の概要】
【発明が解決しようとする課題】
【0007】
しかしながら、上記の複合ソート法は、ハミング距離(hamming distance)でショートリードの類似性を評価することを前提としている。ここで、ハミング距離とは、等しい長さの2つのショートリードにおいて、対応する位置にある異なった文字の個数であり、換言すれば、一方のショートリードを他方のショートリードに一致させるために必要な文字の置換回数を計測したものである。
【0008】
ハミング距離によるショットリードの類似性の評価では、配列の一部の塩基が欠損している場合には、配列の全体がずれるため、類似性を適切に判断できない場合がある。遺伝子を構成する一部の塩基が欠損したり、新しい塩基が挿入されたりすることは起こり得ることなので、このような配列どうしは類似すると判断した方が良い場合もある。
【0009】
配列の一部の欠損や挿入をも考慮して、配列どうしの距離を適切に測る方法として編集距離(edit distance)がある。編集距離は、ハミング距離と同様に、2つのショートリードの類似度を示す数値であるが、文字の挿入、削除、または置換によって、一方のショートリードを他方のショートリードに一致させるために必要な操作の最小回数である。編集距離でショートリードの類似性を評価する場合には、上記の複合ソート法をそのまま採用することはできない。
【0010】
本発明は、大量のショートリードの中から、高速に、編集距離または最大ギャップ(挿入および欠損の総和の閾値)が所定の範囲内にあるショートリードのペアを探索することができる配列解析装置および配列解析方法を提供することを目的とする。
【課題を解決するための手段】
【0011】
本発明の配列解析装置は、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置であって、探索対象の複数の塩基配列を入力する配列入力部と、編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力する条件入力部と、前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成部と、前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内であると計算されたペアを示すデータを出力する類似判定部とを備えた構成を有している。
【0012】
本発明者らは、膨大な量の塩基配列の中から所定の編集距離d以内のペアを探索する計算コストの低い方法を鋭意研究した結果、塩基配列を(d+n)個のブロックに分割すると、少なくともn個のブロックについては、所定の最大オフセットの範囲内でオフセットさせることにより部分配列が一致することを見出した。この知見に基づいて、最初に、前述のような条件(編集距離がd以内であることの必要条件)を満たす塩基配列の集合を等価クラスとして求め、等価クラスに含まれる塩基配列どうしの編集距離を計算する本発明を完成させた。本発明の構成によれば、等価クラス生成部から等価クラスとして報告された塩基配列の集合によってできるすべてのペアの中には、編集距離がd以下であるペアがもれなく含まれているので、計算コストを削減した上で、所定の編集距離以内にあるすべてのペアを見つけ出すことができる。
【0013】
また、本発明の別の態様の配列解析装置は、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置であって、探索対象の複数の塩基配列を入力する配列入力部と、編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)と、前記配列入力部に入力される前記複数の塩基配列の挿入および欠損の総和の閾値である最大ギャップg(gはd以下の自然数)を入力する条件入力部と、前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記最大ギャップgによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成部と、前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内かつ、挿入および欠損の総和がg以内であると計算されたペアを示すデータを出力する類似判定部とを備えた構成を有している。
【0014】
また、上記の配列解析装置において、前記等価クラス生成部は、前記n個のブロックのうちの1つのブロックについて前記窓枠を前記最大オフセット内でオフセットを変えて読み出した部分配列または当該ブロック内の部分配列が一致するという絞り込み条件を満たす塩基配列の集合を、前記等価クラスの候補となる候補クラスとして求め、当該候補クラスに含まれる塩基配列に対して、前記n個のブロックのうちの別のブロックについて前記窓枠を前記最大オフセット内でオフセットを変えて読み出した部分配列または当該別のブロック内の部分配列が前記絞り込み条件を満たす候補クラスを求めるという処理を、前記n個のブロックについて順次行い、n個のすべてのブロックを用いて絞り込まれた候補クラスを等価クラスとして生成する。このように候補クラスを順次絞り込んでいくことにより、n個のすべてのブロックについて、所定の最大オフセットの範囲内でオフセットさせることにより部分配列が一致するという条件を満たす塩基配列の集合を適切に求めることができる。
【0015】
また、上記の配列解析装置において、前記等価クラス生成部は、塩基配列に対して、絞り込みに用いられるブロックから前記最大オフセット内のあらゆるオフセット値で前記窓枠をオフセットさせて部分配列を読み出すオフセット部と、前記窓枠をオフセットさせて読み出された部分配列および前記窓枠をオフセットしないで読み出された部分配列を塩基の並び方に応じてソートして同じ配列パターンの部分配列が続く範囲をグループ化するソート部と、前記ソート部にて生成されたグループに複数の部分配列が含まれ、かつ、当該グループの中の少なくとも1つの部分配列のオフセット値が0であるときに、そのグループを候補クラスとして求める候補クラス生成部とを備える。少なくとも1つの部分配列のオフセット値が0でない場合には、少なくとも1つの部分配列のオフセット値が0である候補クラスが他に存在するので、候補クラスとしないことにより、重複する候補クラスの生成を回避できる。
【0016】
また、上記の配列解析装置において、前記候補クラス生成部は、候補クラスに含まれる塩基配列に関連付けて、候補クラスを絞り込んだブロックの位置、および、候補クラスに含まれる他の部分配列と一致するために用いたオフセット値をメモリに記憶し、前記n個のブロックのうちのm(m<n)個目のブロックを用いて候補クラスを求める際には、前記1つのグループの中で、オフセット値が0であると判定された塩基配列の中に、前のm−1個のブロックを用いた候補クラスの絞り込みにおいて、当該塩基配列の部分配列のオフセット値が連続して0であった塩基配列が存在しない場合には、当該グループを候補クラスとしない。
【0017】
本発明では、ブロックごとに部分配列のオフセット値が独立しているので、同じ塩基配列の中でも、あるブロックにおいてオフセット値を−1とした部分配列が一致し、別のブロックではオフセット値を+1とした部分配列が一致するということがあり得る。しかし、このような塩基配列は、全体として見れば、オフセットの与え方が一貫していないため、同じ候補クラス内にある他の塩基配列と所望の編集距離以内にはない可能性が高い。本発明の構成によれば、オフセット値が連続して0である塩基配列が一つも含まれていない場合には候補クラスとしないことにより、所望の編集距離のペアが存在しない等価クラスの生成をあらかじめ防止することができる。
【0018】
また、上記の配列解析装置において、前記候補クラス生成部は、1つのグループ内に、オフセット値の異なる同じ塩基配列の部分配列が存在する場合において、それらのオフセット値の中に0が含まれるときは0を、含まれないときはそれらのオフセット値の中から任意に選択した値を当該部分配列のオフセット値として記憶する。
【0019】
また、上記の配列解析装置において、前記類似判定部は、前記メモリから、前記等価クラスに含まれる塩基配列のブロックの位置、および、オフセット値を読み出し、等価クラスに含まれる塩基配列の任意のペアのうち、少なくとも一方の塩基配列のオフセット値がすべてのブロックで0であり、かつ、所定のブロック順位規則に基づく順位において前記ブロックの位置よりも前に、前記最大オフセット内でオフセットさせることにより一致する部分配列がないときに、当該ペアの塩基配列について、その編集距離を計算する。
【0020】
同じ塩基配列のペアが複数の異なる等価クラスに含まれる場合があるが、本発明の構成によれば、それらの中から代表となる等価クラスを一つ決定することができるため、代表となる等価クラスに含まれるペアについてのみ編集距離の計算を行うことが可能となり、重複した計算処理を回避できる。
【0021】
また、上記の配列解析装置において、前記等価クラス生成部は、前記候補クラス生成部によって複数の候補クラスが生成された場合には、そのうちの1つの候補クラスについて次のブロックを用いた候補クラスの絞り込みを行なう処理を、等価クラスが求まるか、または候補クラスがなくなるまで繰り返し行い、等価クラスが求まるか、または候補クラスがなくなった時点で1つ前のブロックを用いた処理で求まった他の候補クラスについて絞り込みを行う。
【0022】
また、上記の配列解析装置において、前記等価クラス生成部は、最初のブロックを用いて求められたすべての候補クラスについての絞り込みを終了したときに、組み合わせの異なるn個のブロックを選択し、新たに選択されたn個のブロックを用いて等価クラスを生成する。これにより、n個のブロックの組み合わせを順次変えて、編集距離d以内にある塩基配列のペアをもれなく含む複数の等価クラスを求めることができる。
【0023】
また、上記の配列解析装置において、前記等価クラス生成部は、前記新たに選択されたn個のブロックを用いて等価クラスを生成する場合において、前記候補クラス生成部において利用できる候補クラスが生成されているときは、当該既に生成されている候補クラスを利用することで、前記ソート部において同一の候補クラスが複数回整列されるのを回避する。これにより、前記ソート部において、同一の候補クラスが複数回整列されることが回避され、再帰的計算を用いてすでに得られた候補クラスを再利用できる。
【0024】
また、上記の配列解析装置において、前記等価クラス生成部は、前記候補クラスが生成されるごとに、当該候補クラスを利用するn個のブロックの組合せのすべてについて、等価クラスが求まるまで絞込みを行い、当該候補クラスを利用するすべてのn個のブロックの組合せについて等価クラスを求めた後に、当該候補クラスをメモリ上から解放する。これにより、前記ソート部において同一の候補クラスが複数回整列されるのを回避しつつ、作成された候補クラスを再利用するために候補クラスを記憶するために必要なメモリ容量を小さく抑えることができる。
【0025】
また、上記の配列解析装置において、前記配列入力部に入力される前記複数の塩基配列は所定の範囲D内で長さが異なっており、前記条件入力部に入力される編集距離dは、d≧Dを満たし、前記最大オフセットは、d/2以上の最小の整数である。
【0026】
また、上記の配列解析装置において、前記配列入力部に入力される前記複数の塩基配列は所定の範囲D内で長さが異なっており、前記条件入力部に入力される編集距離dは、d≧Dを満たし、前記最大オフセット量は、g/2以上の最小の整数である。
【0027】
また、上記の配列解析装置において、前記配列入力部に入力される前記複数の塩基配列は長さが互いに等しく、前記最大オフセットは、d/2以下の最大の整数である。
【0028】
また、上記の配列解析装置において、前記配列入力部に入力される前記複数の塩基配列は長さが互いに等しく、前記最大オフセットは、g/2以下の最大の整数である。
【0029】
また、上記の配列解析装置において、前記数値nまたは前記分割数(d+n)は、前記配列解析装置が所定のアルゴリズムに従って決定して前記条件入力部に入力される。
【0030】
また、本発明の別の態様は、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析方法であって、この方法は、探索対象の複数の塩基配列を入力する配列入力ステップと、編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力する条件入力ステップと、前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成ステップと、前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内であると計算されたペアを示すデータを出力する類似判定ステップとを含む構成を有している。
【0031】
また、本発明の別の態様の配列解析方法は、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析方法であって、この方法は、(1)探索対象の複数の塩基配列を入力するステップと、(2)編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力するステップと、(3)前記複数の塩基配列の各々を前記分割数(d+n)個のブロックに分割するとともに各ブロックにブロック番号を付与するステップと、(4)前記(d+n)個のブロックの中からn個のブロックをブロックセットとして選択するステップと、(5)作業数を1とし、前記選択されたn個のブロック中のブロック番号が最小のブロックを作業ブロックをとし、かつ、ステップ(1)で入力されたすべての塩基配列を候補クラスとするステップと、(6)候補クラス内のすべての塩基配列について、前記作業ブロック内の部分配列と、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該作業ブロックからオフセットさせて読み出した部分配列からなる部分配列の集合を生成するステップと、(7)前記部分配列の集合を塩基の並び方に応じてソートしてグループ集合を生成するステップと、(8)前記グループ集合を順にスキャンしていき、同じ配列パターンの部分配列が続く範囲を1つのグループとして、そのグループに複数の部分配列が含まれ、かつ、少なくとも1つの部分配列のオフセット値が0であるときに、作業数がnであるかを判断するステップと、(9)ステップ(8)において、作業数がnでないときに、そのグループを候補クラスとして報告するステップと、(10)ステップ(9)の後に、作業数をインクリメントし、前記作業ブロックを、ステップ(4)にて選択されたn個のブロックのうちの次のブロックに変更して、ステップ(6)に戻るステップと、(11)ステップ(8)において、作業数がnであるときに、そのグループを等価クラスとして報告するステップと、(12)ステップ(11)の後に、当該等価クラスに含まれる塩基配列の任意のペアのうちの少なくとも一方において、オフセット値の履歴にゼロ以外の値がなく、かつ所定の正準性判断に基づいて正準であると判断されたペアにつき、編集距離を計算するステップと、(13)ステップ(12)の後に、編集距離がd以内である場合に、当該ペアを示すデータを出力するステップと、(14)ステップ(13)の後に、前記グループ集合に未スキャン部分があるかを判断するステップと、(15)ステップ(14)において、未スキャン部分がある場合に、ステップ(8)に戻るステップと、(16)ステップ(14)において、未スキャン部分がない場合に、作業数をデクリメントして、作業ブロックを1つ前のブロックに戻すステップと、(17)ステップ(16)の後に作業数が1以上である場合に、ステップ(14)に戻るステップと、(18)ステップ(16)の後に作業数が0になった場合に、前記グループ集合を消去してまだ選択されていないブロックセットがあるかを判断するステップと、(19)ステップ(18)において、まだ選択されていないブロックセットがある場合に、ステップ(4)に戻るステップと、(20)ステップ(18)において、まだ選択されていないブロックセットがない場合に、処理を終了するステップとを含む構成を有している。
【0032】
また、上記の配列解析方法は、さらに、(21)ステップ(4)の後に、前記ブロックセットに含まれるブロックをブロック番号が小さいものから並べて、ブロック番号が小さい方から1または複数のブロックにより構成されるサブセットの候補クラスが記憶されているか否かを判断するステップと、(22)ステップ(21)で、候補クラスが記憶されていると判断された場合に、ステップ(5)に代えて、当該サブセットの候補クラスを読み出して、その候補クラスに対応する作業数および作業ブロックを設定するステップとを含んでいる。
【0033】
また、本発明のさらに別の態様の配列解析方法は、複数の塩基配列の中から、所定の範囲内にある塩基配列のペアを探索する配列解析方法であって、(51)探索対象の複数の塩基配列を入力するステップと、(52)編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力するステップと、(53)前記複数の塩基配列の各々を前記分割数(d+n)個のブロックに分割するとともに、各ブロックにブロック番号を付与するステップと、(54)作業数Wを1とし、作業ブロックbkをブロック1とし、作業開始ブロックSを1とするステップと、(55)絞込み対象候補クラスAを入力されたすべてのショートリードとするステップと、(56)絞込み対象候補クラスAの作業ブロックbkを展開・ソートして、グループ集合G(bk)を生成して記憶し、作業数Wの作業済み最右のブロックP(W)を作業ブロックbkとするステップと、(57)グループ集合G(bk)を上からスキャンして、候補クラスの条件を満たす最初のグループを選択してそのグループを絞り込み対象候補クラスAとし、かつ、作業数Wの絞込み対象候補クラスQ(W)をその絞り込み対象候補クラスAとするステップと、(58)ステップ(57)の後に、作業数Wがnであるか否かを判断するステップと、(59)ステップ(58)にて作業数Wがnでない場合に、作業数Wをインクリメントするステップと、(60)ステップ(59)の後に、bk+1がd+W以下であるか否かを判断するステップと、(61)ステップ(60)にて、bk+1がd+W以下である場合に、作業ブロックbkをインクリメントして、ステップ(56)に戻るステップと、(62)ステップ(58)にて作業数Wがnである場合に、絞込み対象候補クラスAを等価クラスとするステップと、(63)ステップ(62)の後に、等価クラスについて、類似判定をするステップと、(64)ステップ(63)の後に、グループ集合G(bk)に未スキャン部分があるか否かを判断し、未スキャン部分がある場合に、ステップ(57)に戻るステップと、(65)ステップ(64)にて、未スキャン部分がない場合に、グループ集合G(bk)をメモリ上から解放し、作業数Wをデクリメントするステップと、(66)ステップ(65)の後に、W=0であるか否かを判断するステップと、(67)ステップ(66)にて、W=0でない場合に、絞込み対象候補クラスAを作業数Wの絞込み対象候補クラスとするステップと、(68)ステップ(60)にて、bk+1がd+Wより大きい場合に、作業数Wをデクリメントし、作業ブロックbkをそのデクリメントされた作業数の作業済み最右のブロックP(W)とするステップと、(69)ステップ(66)にて、W=0である場合に、作業開始ブロックSが編集距離d以下であるか否かを判断し、作業開始ブロックSが編集距離dより大きい場合に、処理を終了するステップと、(70)ステップ(69)にて、作業開始ブロックSが編集距離d以下である場合に、作業開始ブロックSをインクリメントして、作業数Wを1にリセットして、ステップ(55)に戻るステップとを有する。
【0034】
本発明のさらに別の態様は、コンピュータで実行されることにより、コンピュータに上記の方法を実行させるコンピュータプログラムである。
【発明の効果】
【0035】
本発明の配列解析装置によれば、計算コストを削減して、高速に、所定範囲内にあるすべてのショートリードのペアを見つけ出すことができる。
【図面の簡単な説明】
【0036】
【図1】本発明の実施の形態における配列解析装置の入力と出力を説明する図
【図2】(a)ショートリードのペアを示す図 (b)編集距離を説明する図
【図3】(a)編集距離の定理を説明する図(挿入と削除の総数の上限が指定されていない場合) (b)編集距離の定理を説明する図(挿入と削除の総数の上限が指定されている場合)
【図4】(a)配列の例を示す図 (b)編集距離の定理を具体例で説明する図
【図5】本発明の実施の形態における配列解析装置のブロック図
【図6】本発明の実施の形態におけるブロック分割の例を示す図
【図7】本発明の実施の形態におけるオフセット部の動作を説明する図
【図8】本発明の実施の形態におけるオフセット履歴の情報が付加された部分配列の集合を示す図
【図9】本発明の実施の形態におけるソート部によってソートされた部分配列を示す図
【図10】本発明の実施の形態における候補クラス記憶部に記憶される情報を示す図
【図11】本発明の第1の実施の形態におけるブロックマスクによる距離計算対象の絞込みを説明する図
【図12】本発明の第1の実施の形態におけるブロック1、ブロック2、ブロック3がブロックセットとして選択された場合の絞込み手順を説明する図
【図13】本発明の第1の実施の形態におけるブロック1、ブロック2、ブロック4がブロックセットとして選択された場合の絞込み手順を説明する図
【図14】(a)本発明の第1の実施の形態における等価クラスの情報および等価クラスに含まれるショートリードの配列パターンの例を示す図 (b)本発明の第1の実施の形態における正準性判定の条件3の判定例を示す図 (c)本発明の第1の実施の形態におけるS22とS75のペアの例を示す図 (d)本発明の第1の実施の形態におけるS56とS75のペアの例を示す図
【図15】本発明の第1の実施の形態における配列解析装置の処理フロー図
【図16】本発明の第2の実施の形態における絞込み手順を説明する図
【図17】本発明の第2の実施の形態における絞込み手順を説明する図
【図18】本発明の第2の実施の形態における絞込み手順を説明する図
【図19】本発明の第2の実施の形態における絞込み手順を説明する図
【図20】本発明の第2の実施の形態における配列解析装置の処理フロー図
【図21】従来の複合ソート法を説明する図
【発明を実施するための形態】
【0037】
以下、本発明の実施の形態について説明する。本発明の実施の形態の配列解析装置は、ほぼ同じ長さの膨大な量のショートリードが与えられたときに、その中から編集距離d以内のペアを類似のペアとして、重複なく、完全に列挙することを目的としている。配列解析装置に入力されるショートリードには、それぞれ固有のIDが付与されている。本明細書では、ショートリードに付与されたIDを「インデクス」と称し、「S」に添え字を付す形式で表記する。ショートリードは、塩基配列であるので、「A」、「C」、「G」、「T」という文字の組合せからなる配列である。
【0038】
図1に示すように、本実施の形態の配列解析装置は、入力として複数のショートリード(S1,・・・,SN)が与えられたときに、編集距離d以内のペア{S1, S2}、{S1, S100}、 {S2, S10}、・・・を出力する装置である。編集距離としては、例えば、挿入と削除の総和が5以内で、編集距離が10以内というように、挿入と削除の総和の最大値を指定することもできる。
【0039】
配列解析装置が扱う複数のショートリードは、完全に長さが等しいものである必要はない。但し、入力される複数のショートリードの中で最も短いショートリードと最も長いショートリードとの長さの差が編集距離d以内であることが条件となる。
【0040】
以下、編集距離について詳しく説明する。編集距離とは、上述のように、文字の挿入、削除、または置換によって、一方のショートリードを他方のショートリードに一致させるために必要な操作の最小回数である。例えば、図2(a)のショートリードのペア(S1とS2)があった場合に、図2(b)に示すように、このS1の第2文字目を削除して、第6文字目と第7文字目の間にGを挿入し、さらに第11文字目のTをAに置換するとS2と一致する。よって、この例では、S1について、1回の削除、1回の挿入、1回の置換という3回の操作によってS2と一致することになるので、S1とS2の編集距離は3である。
【0041】
編集距離については、次の定理が成り立つ。いま、文字列S1と文字列S2の編集距離がdであるとすると、図3(a)に示すように、S1、S2を中心をそろえて並べ、d+n個(nは任意の自然数)のブロックに分割した場合、n個のブロックが、ずれ幅d/2以内で完全一致する。ここで、d/2が割り切れない時は、切り上げる(ただし、S1とS2の長さが同じ場合は、切り捨てる)。例えば、d=3のときは、一般的には、ずれ幅2以内となる(ただし、S1、S2が同じ長さのときはずれ幅1以内となる)。
【0042】
挿入と削除の総数の上限g(gはg≦dを満たす自然数)が指定されているときは、図3(b)に示すように、S1、S2を中心をそろえて並べ、d+n個(nは任意の自然数)のブロックに分割した場合、少なくともn個のブロックでは、ずれ幅g/2以内で完全一致する。ここで、g/2が割り切れない時は、切り上げる(ただし、S1とS2の長さが同じ場合は、切り捨てる)。例えば、g=3のときは、一般的には、ずれ幅2以内となる(ただし、S1、S2が同じ長さのときはずれ幅1以内となる)。
【0043】
上記の定理を具体例で説明する。編集距離d=3、n=2として、図4(a)に示す配列S1、S2が与えられたとする。この場合は、配列を(d+n=)5ブロックに分割する。この5つのブロックの中からn(=2)個のブロックを選択する組合せは、52=10通りある。上記の定理より、S1、S2が編集距離d(=3)以内であるのならば、10通りある組合せの中で、少なくとも1つの組合せにおいて、2つのブロックがずれ幅1以内で一致することになる。
【0044】
実際には、図4(b)に示すように、S1とS2は、ブロック1とブロック3、ブロック1とブロック4、ブロック3とブロック4の組合せで、ずれ幅1以内で部分配列が一致する。具体的には、読出し範囲を特定する窓枠をずらして部分配列を読み出し、部分配列の配列パターンを比較したときに、配列パターンが一致すると判断している。即ち、ブロック1については、窓枠のずれ量をゼロとしたときに、窓枠内の配列パターンが「AT」で一致し、ブロック3については、配列S1に対してS2の窓枠をずれ量−1だけずらしたときに、窓枠内の配列パターンが「AGC」で一致し、ブロック4については、窓枠のずれ量をゼロとしたときに、窓枠内の配列パターンが「GAT」で一致する。
【0045】
上記の定理を用いて、本実施の形態の配列解析装置では、複数のショートリードの中から編集距離d以下のショートリードのペアを探索するために、入力された複数のショートリードの各々を(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個(nは1以上の自然数)のすべてのブロックにおいて、窓枠を上記のように定義されるずれ幅(以下単に「最大オフセット」という)の範囲内で窓枠内の配列パターンが一致するという条件を満たすショートリードの集合を、等価クラスとして出力する。選択するn個のブロックの組合せを変えることで、次々に等価クラスが出力される。
【0046】
この等価クラス内のショートリードのすべてのペアが編集距離d以内にあるわけではないが、仮に、この等価クラス内のショートリードのすべてのペアについて編集距離d以内にあるかを実際に計算したとしても、シーケンサから与えられた複数のショートリードのすべてのペアに対して計算をする場合と比較すると格段に計算すべきペアの数は減少することになる。本実施の形態では、後述するように、異なる等価クラスに含まれる同一のショートリードのペアについては、重複して編集距離を計算・報告しないようにしている。
【0047】
以下、図面を参照して本実施の形態の配列解析装置を説明する。図5は、本実施の形態の配列解析装置の構成を示すブロック図である。図5において、配列解析装置10は、配列入力部11と、条件入力部12と、等価クラス生成部13と、類似計算部14とを備えている。なお、本実施の形態の配列解析装置10は、以下に説明する各処理を実行するモジュールを有するプログラムをコンピュータによって実行することによって実現することとしてもよい。このようなプログラムも本発明に含まれる。
【0048】
配列入力部11には、シーケンサから出力された複数のショートリードのデータ(以下、単に「ショートリード」ともいう)が入力される。配列入力部11は、入力されたショートリードを等価クラス生成部13に出力する。ショートリードの数量は、例えば数千万個である。また、ショートリードの長さは、例えば35〜50程度である。なお、本発明は、ショートリードの数量が多いほど計算コスト削減の効果が顕著になる。また、ショートリードの長さは、上記の長さに限られるものではなく、より短くてもよく、より長くてもよい(例えば200以上であってよい)。
【0049】
入力されるショートリードは、すべて同一の長さであってもよく、互いに長さが異なっていてもよい。但し、編集距離d以内にあるショートリードのペアを探索する場合には、入力されるショートリードの最長長さと最短長さとの差分Dはd以下でなければならない。
【0050】
シーケンサから配列解析装置10に与えられたショートリードが配列解析装置10内のハードディスク等の記憶装置に保存されてもよい。この場合は、配列入力部11は記憶装置からショートリードを読み出して、等価クラス生成部13に出力する。
【0051】
条件入力部12には、探索したいペアの編集距離の最大値である編集距離dの値、ショートリードをブロックに分割する際の分割数bをb=d+nとして決定するためのnの値が入力される。nは1以上の任意の自然数である。nはユーザが任意の値を入力してもよいし、ショートリードの長さや数量等に応じて所定のアルゴリズムで計算された値が入力されてもよい。編集距離dおよび最大ギャップgは、ユーザが所望する任意の値を入力する。配列解析装置10は、nを計算して求める構成を有して、計算により得られた結果を条件入力部12に入力するようにしてもよい。
【0052】
等価クラス生成部13は、配列入力部11に入力された複数のショートリードに対して、条件入力部12に入力されたdおよびnの値に基づいて、編集距離d以内にあるペアが存在し得る複数のショートリードの集合である等価クラスを出力する。条件入力部12に編集距離dとともに最大ギャップgが入力された場合には、等価クラス生成部13は、配列入力部11に入力された複数のショートリードに対して、条件入力部12に入力されたgおよびnの値に基づいて、最大ギャップg以内にあるペアが存在し得る複数のショートリードの集合である等価クラスを出力してもよい。以下の説明では、主に条件入力部12に編集距離dが入力され最大ギャップgは入力されていない場合について説明する。
【0053】
等価クラス生成部13は、ブロック分割部131と、ソート対象指定部132と、オフセット部133と、ソート部134と、候補クラス生成部135と、候補クラス記憶部136とを有している。ブロック分割部131は、配列入力部11より入力されたショートリードの各々を、対応する部分配列ごとに比較する処理で用いる処理単位をブロックとして、当該ブロックごとに分割した状態で記憶している。このときのブロック数は、条件入力部に入力されたdおよびnの値を用いて(d+n)とする。
【0054】
この結果、各ショートリードは、(d+n)個の部分配列に分割される。ブロックは左から順に順位が定義されており、最左のブロックをブロック番号1番として、右にいくにつれて番号が増加するように各ブロックにブロック番号が付与される。本明細書では、ブロック番号k番のブロックを「ブロックk」と表記する。各ショートリードの同一のブロック番号の部分配列は以下の処理においてブロックが対応する部分配列として扱われる。
【0055】
上述のように、一般的にはショートリードの長さはすべて等しいわけではないので、ブロック分割部131は、各ショートリードの中心をそろえて(d+n)個のブロックに分割する。ショートリードの長さが奇数の場合と偶数の場合とでは、中心は完全にそろわないが、その場合にも、すべての奇数のショートリードの中心がすべての偶数のショートリードの中心の例えば左側に来るようにして、全体として中心のずれが0.5文字になるようにした上で、同じ列に並ぶ文字が同じブロックに所属するように、すべてのショートリードを分割する。なお、各ブロックの幅、即ち部分配列の長さは必ずしも互いに等しくなくてもよい。
【0056】
図6は、ブロック分割の例を示す図である。図6の例は、S1〜S3のショートリードがあった場合において、S1〜S3の長さがそれぞれ15、14、13であったときのブロック分割の例を示している。図6に示すように、奇数長であるS1およびS3の中心は、それぞれ第8番目および第7番目の文字の中心であり、偶数長であるS2の中心は、第7番目の文字と第8番目の文字の間である。この場合に、奇数長のショートリードについては、偶数長のショートリードよりも左に半文字分ずらすことで、奇数長のショートリードおよび偶数長のショートリードを含むすべてのショートリードの中心がそろうことになる。
【0057】
ソート対象指定部132は、処理の過程においてどのショートリードのどのブロックをソート対象とするかを指定する。等価クラス生成部13は、入力された複数のショートリードをブロックに分割して、ブロックごとに評価することで、入力された複数のショートリードのあらゆる組合せのペアから編集距離がd以内にないペアを除外していく絞込みを行う。換言すれば、ソート対象指定部132は、この絞り込みに用いられるブロック内に含まれる絞り込み対象のショートリードの部分配列を指定する。絞込みの詳細は後述する。
【0058】
オフセット部133は、ソート対象指定部132により指定されたショートリードの集合に含まれる各ショートリードについて、ソート対象指定部132により指定されたブロックから窓枠をオフセットして読み出した部分配列を図示しないメモリ上に展開する。このために、オフセット部133は、まず、ソート対象指定部132により指定されたショートリードの集合に含まれる各ショートリードの、ソート対象指定部132により決定されたブロックの部分配列を読み出すとともに、当該ブロックを基準として最大オフセット以下の各オフセット値で窓枠をずらして、部分配列を読み出す。
【0059】
このとき、ブロック分割部131から読み出された各部分配列には、その部分配列が属するショートリードのインデクスの情報とオフセット値の情報が付加されている。以下、オフセット部133が上記のようにブロックから窓枠をずらして読み出した部分配列を「オフセット部分配列」ともいう。また、ブロック分割部131にて分割されて生成されたブロック内に収まる部分配列を「オフセット部分配列」と特に区別するときは、ブロック分割部131にて分割されて生成されたブロック内に収まる部分配列を「原部分配列」ともいう。両者を区別しないときは、単に「部分配列」という。
【0060】
最大オフセットは、上述のように、ショートリードの長さが均一であるか否か、および編集距離d、または挿入と削除の総数の上限が分かっている場合にはその総数の上限gを用いて決定される。例えば、ショートリードの長さが10文字以内でばらついている場合において、挿入と削除の総数がg=5以内である場合には、g/2=2.5となるが、ショートリードの長さがばらついているので、2.5を切り上げた3が最大オフセットとなる。この場合、オフセット部133は、ソート対象指定部132により指定されたショートリードの各々について、ソート対象指定部132により指定されたブロックから、窓枠を−3ずらして読み出したオフセット部分配列、窓枠を−2ずらして読み出したオフセット部分配列、窓枠を−1ずらして読み出したオフセット部分配列、窓枠をずらさずに読み出した原部分配列、窓枠を+1ずらして読み出したオフセット部分配列、窓枠を+2ずらして読み出したオフセット部分配列、および窓枠を+3ずらして読み出したオフセット部分配列を生成する。
【0061】
図7は、オフセット部133の動作を説明する図である。図7に示すように、ソート対象指定部132がショートリードの集合{S1、S24、S97、S257、……}の「ブロック2」をソート対象として指定すると、オフセット部133は、ブロック分割部131から該当する部分配列を読み出す。この例では、最大オフセットは1であるので、オフセット部133は、ブロック2の原部分配列とともに、ブロック2から−1ずれた位置の部分配列およびブロック2から+1ずれた位置の部分配列をそれぞれオフセット部分配列として読み出す。
【0062】
図7において、「SR」はショートリードのインデクスを表し、「p」はオフセット値を表す。図に示すように、オフセット部133によってブロック分割部131から読み出される原部分配列およびオフセット部分配列には、それぞれ、当該部分配列が所属するショートリードのインデクスとともに、当該部分配列が指定のブロックから窓枠をどれだけずらして読み出されたかを示すオフセット値が付加されている。なお、p=0はオフセット値がゼロ、即ちその部分配列が原部分配列であることを示している。
【0063】
オフセット部133は、各部分配列に対して、その部分配列を読み出したときのオフセット値に、その部分配列が所属するショートリードの過去のすべてのオフセット値を加えて、オフセット履歴として付加する。図8は、過去のすべてのオフセット値と当該部分配列を読み出したときのオフセット値とがオフセット履歴として付加された部分配列の集合を示している。図8の例では、過去のブロック1において、S1、S24、S97、S257は、それぞれオフセット0、1、0、0で候補クラス(後述する)として判断されたため、それらの過去のオフセット値が、各部分配列を読み出した際のオフセット値とともにオフセット履歴として、各部分配列に付加されている。
【0064】
ソート部134は、オフセット部133にて生成されたオフセット部分配列および原部分配列のすべてを塩基の並び順に応じてソートし、整列された部分配列の集合を形成する。このとき、ソート部134は、基数交換ソートによって部分配列をソートするので、グループ集合では同一の配列パターンを有する部分配列はグループ化されている。以下、ソートされて同一の配列パターンを有する部分配列ごとにグループ化された部分配列の集合を「グループ集合」という。具体的には、部分配列は上述のようにA、C、G、Tの組み合わせからなる配列であるので、ソート部134は、これを辞書順に並べ替える。このように、部分配列をソーティング(整列)するのは、原部分配列およびオフセット部分配列の中で一致する部分配列をまとめるためである。図9は、グループ集合の例を示す図である。
【0065】
候補クラス生成部135は、グループ集合を上から順にスキャンして、ソート部にて生成されたグループごとに注目し、注目するグループが候補クラスであるか否かを判断する。候補クラスであると判断する条件は、次の条件1および条件2である。
条件1:そのグループに複数の部分配列が含まれていること
条件2:少なくとも1つの部分配列のオフセット履歴に含まれるオフセット値がすべて0であること
これらの条件1および条件2をいずれも満足する場合には、当該グループは候補クラスと判断される。
【0066】
図9の例では、条件1および条件2を同時に満足するグループは、文字列「GCT」のグループおよび文字列「TGC」のグループである。従って、図9の例では、候補クラス生成部135は、この「GCT」のグループおよび文字列「TGC」のグループを候補グループであると判断する。なお、「TTG」のグループは、複数の部分配列が含まれてはいるが、オフセット履歴に含まれるオフセット値がすべて0であるという条件を満たす部分配列がないため、候補クラスとは判断されない。これに対して、「GCT」のグループでは、S1のオフセット履歴に含まれるオフセット値がすべて0であり、「TGC」のグループは、S257のオフセット履歴に含まれるオフセット値がすべて0であるので、これらのグループは候補グループであると判断される。
【0067】
候補クラス生成部135は、候補クラスであると判断したグループに含まれる部分配列のインデクスをソート対象指定部132に報告する。また、候補クラス生成部135は、候補クラスであると判断したグループに含まれる部分配列およびそのオフセット履歴に、その候補クラスが生成されるまでにソート対象となったブロックの履歴と関連付けて候補クラス記憶部136に保存する。
【0068】
なお、このとき、候補クラス生成部135は、同一のショートリードに由来する複数の部分配列が1つの候補クラスに重複して含まれないようにする。例えば、オフセット「0」の場合に部分配列が「AAA」、オフセット「1」の場合にも部分配列が「AAA」となるような塩基配列である。このように、候補クラス生成部135は、候補クラスと判断したグループの中に、同一のショートリードに由来する複数の部分配列が含まれている場合において、それらの部分配列の中にオフセットが0である部分配列があるときは、当該オフセットが0を、それらの部分配列の中にオフセットが0である部分配列がないときはそれらの部分配列のオフセット値のうちから任意に選択した値のみを当該同一のショートリードに由来する複数の部分配列のオフセット値として採用する。即ち、オフセット履歴中のオフセット値がすべてゼロのものがあれば、そのオフセット履歴を当該ショートリードのオフセット履歴とし、他のオフセット履歴については、重ねてソート対象指定部132に報告したり候補クラス記憶部136に保存したりしないようにする。
【0069】
候補クラス生成部135は、作業数(後述)が最大値(条件入力部12に入力されたnが作業数の最大値となる)のときに候補クラスである判断したショートリードの集合を等価クラスとして類似判定部14に出力(報告)する。具体的には、候補クラス生成部135は、等価クラスに含まれるショートリードのインデクス、その候補クラスが生成されるまでにソート対象となったブロックの履歴、およびオフセット履歴の情報を類似判定部14に出力する。
【0070】
候補クラス記憶部136は、ハードディスクで構成され、上述のように、候補クラス生成部135から報告を受けた候補クラスを記憶する。候補クラスの情報は、その候補クラスが生成されるまでにソート対象となったブロックの履歴と関連付けられて記憶される。このブロックの履歴は、即ち、その候補クラスに含まれるショートリードが最大オフセット内で一致するブロックを示す情報である。また、候補クラス記憶部136に記憶される候補クラスには、当該候補クラスに含まれる各ショートリードのインデクスのそれぞれに対して、当該ショートリードが過去に当該候補クラス内の他のショートリードと一致した各ブロックにおいてどれだけのオフセットで一致したのかを示すオフセット履歴の情報が付加されている。
【0071】
図10は、候補クラス記憶部136に記憶される情報を示す図である。図10において、具体的なデータは、図7〜図9の例に対応している。図10に示すように、一般的には、1つのブロック履歴には複数の候補クラスが対応している。図中、「r」は「root」の略記であり、複数のショートリードの全体を表している。
【0072】
次に、ソート対象指定部132によるソート対象の指定について説明する。編集距離がdであり、配列の分割数がb(b=d+n)であるとすると、上記で説明した定理により、編集距離がd以内にあるショートリードは、必ず、少なくともn個のブロック(b個のブロック中のn個のブロックの組合せは複数とおりあるがそのうちの少なくとも1つの組合せ)において、最大オフセット内で部分配列が一致するはずである。従って、少なくともn個のブロックにおいて最大オフセット内で部分配列が一致するショートリードのグループの中のペアについてのみ実際の編集距離を計算すればよい。そこで、等価クラス生成部13は、b個のブロックの中から選択するn個のブロック組合せを変更しながら、選択されたn個のブロックについて、ショートリードの部分配列が最大オフセット内で一致するか否かを判断するために、ブロックを1つずつ検討して、一致しなかったショートリードを対象から除外していくという、段階的な絞込みを行う。
【0073】
図11は、ブロックマスクによる距離計算対象の絞込みを説明する図である。本明細書では、距離計算対象の絞込みについて、2つの実施の形態を説明する。まず、図11〜図15を参照して、第1の実施の形態を説明する。図11は、編集距離がd=2であり、ブロック番号がb=d+n=5(n=3)である例を示している。図中、丸の中の数字は作業を行うブロック(作業ブロック)のブロック番号を示している。ソート対象指定部132は、5つのブロックのうちの2つをマスクして3つを選択する。この選択の方法は、53(=10)通りある。ソート対象指定部132は、図11に示すようなツリー構造で階層的にブロックを選択していきながら絞込みを行う。
【0074】
ソート対象指定部132は、まず、3つのブロックを選択する(2つのブロックをマスクする)。以下、選択されたブロックの組合せを「ブロックセット」ともいう。ソート対象指定部132は、ブロック番号の若いほうから順にブロックを選択してブロックセットを作る。最初は、ブロック1、ブロック2、ブロック3が選択される。そして、ソート対象指定部132は、入力されたすべてのショートリードについて、未作業であってかつブロック番号が最も若いブロック1をソート対象として指定する。そして、ソート対象指定部132にて指定された対象(すべてのショートリードのブロック1)について、オフセット部133でオフセット部分配列を作成し、ソート部134にて原部分配列と合わせてソートしてグループ集合を生成し、候補クラス生成部135にて候補クラスを生成する。
【0075】
候補クラス生成部135は、候補クラスに含まれる部分配列を含むショートリードのインデクスおよびその各部分配列のオフセット値を候補クラス記憶部136に記憶するとともに、ソート対象指定部132に報告する。ソート対象指定部132は、候補クラス生成部135から報告を受けた候補クラスに含まれるショートリードを次の段階の作業のソート対象とする。ソート対象指定部132は、既にブロック1、ブロック2、ブロック3というブロックセットを選択しているので、未作業であってかつブロック番号が最も若いブロック2を作業ブロックとして選択する。即ち、ソート対象指定部132は、候補クラス生成部135から報告を受けた候補クラスに含まれるショートリードのブロック2を次のソート対象として指定する。
【0076】
オフセット部133は、ソート対象指定部132にて指定されたソート対象について、原部分配列およびオフセット部分配列を含む部分配列を生成し、ソート部134はオフセット部134にて生成された部分配列をソートしてグループ集合を生成し、候補クラス135は候補クラスを生成する。そして、上記と同様にして、候補クラス生成部135は、候補クラスを候補クラス記憶部136に記憶するとともに、ソート対象指定部132に報告する。
【0077】
ソート対象指定部132は、報告を受けた候補クラスに含まれるショートリードの、既に選択したブロック1、ブロック2、ブロック3のうちの未作業かつ最も若いブロックであるブロック3をソート対象として指定する。そして、上記と同様にして、ソート対象指定部132にて指定されたソート対象について、原部分配列およびオフセット部分配列を含む部分配列を生成し、ソート部134はオフセット部134にて生成された部分配列をソートしてグループ集合を生成し、候補クラス生成部135は候補クラスを生成する。
【0078】
このように、すべてのショートリードのあるブロックをソート対象として候補クラスを生成する作業を作業数1の作業と称し、次のソート対象に対する作業を作業数2の作業と称し、このような作業が進んで、一般的に第k回目にソート対象指定部132にてソート対象として決定された対象に対する作業を作業数kの作業という。また、各作業数で作業の対象となるブロックを「作業ブロック」という。図11の例では、n=3であり、ソート対象指定部132は、3つのブロックを選択するので、作業数Wが3になった時点で、選択したブロックすべてについて作業が終了したことになる。
【0079】
候補クラス生成部135は、最終作業数での作業を終了したときに生成された候補クラスを等価クラスとして、その等価クラスに含まれるショートリードのインデクスおよびその等価クラスに含まれるショートリードの、過去の作業ブロックでのオフセット履歴の情報を類似判定部14に出力する。そして、候補クラス生成部135は、ソート対象指定部132に最終作業数で作業が終了したことを報告する。
【0080】
ソート対象指定部132は、選択したブロックセットについて、必要なすべての作業が終わると、選択するブロックを変更して新たなブロックセットを決定する。本実施の形態では、ソート対象指定部132は、ブロック1、ブロック2、ブロック3というブロックセットの次にブロック番号が小さい組合せであるブロック1、ブロック2、ブロック4を選択して新たなブロックセットとする。
【0081】
以下同様にして、この3つのブロックについて段階的な絞込みの作業を行い、10通りのすべての組合せのブロックセットについて、上記のような段階的な検索を行うことで、すべての等価クラスが報告される。しかしながら、次のブロック1、ブロック2、ブロック4についてみると、作業数1の作業として、入力されたすべてのショートリードのブロック1をソート対象とした作業を行い、次に、作業数2の作業として、作業数1の作業で生成された候補クラスのブロック2をソート対象とした作業を行って候補クラスを生成するという一連の作業は、ブロック1、ブロック2、ブロック3のブロックセットについて作業をした際に既に行っており、この一連の作業により作成された候補クラスは、候補クラス記憶部136に記憶されている。
【0082】
そこで、ソート対象指定部132は、ブロック1、ブロック2、ブロック4をブロックセットとして選択した際には、すでにブロック1からブロック2への絞込みで作成されている候補クラスを利用して作業数3の作業をブロック4に対して行うように、ソート対象を指定する。具体的には、ソート対象指定部132は、すでにブロック1からブロック2への絞込みで作成されている候補クラスに含まれるショートリードのブロック4をソート対象として指定する。
【0083】
一般的には、ソート対象指定部132は、既に作成されている候補クラスを利用できるか否かを次のようにして判断する。すなわち、ソート対象指定部132は、複数のブロックの組み合わせを新たなブロックセットとして選択したときに、この複数のブロックをブロック番号が小さいものから順に並べて一番最後のブロックを除いたブロックの組合せで既に作業が行われていないかを判断する。この判断のために、ソート対象指定部132は、ブロックの組合せのうち、作業済みのブロックの組合せを記憶しておいてもよいし、または、新たに選択したブロックの組合せをブロック番号が若いものから順に並べて一番最後のブロックを除いたブロックの組合せについて、候補クラス記憶部136に候補クラスが記憶されているかを直接確認してもよい。
【0084】
ソート対象指定部132は、選択した複数のブロックをブロック番号が若いものから順に並べて一番最後のブロックを除いたブロックの組合せが作業済みでなかった場合には、さらに、残ったブロックの組合せの中で最もブロック番号が大きいブロックを除いて、既に作業が済んでいないかを判断する。このようにして選択した複数のブロックをブロック番号が大きいものから順に除いていって、最後まで作業済みのブロックの組合せがない場合には、上記で説明したブロック1、ブロック2、ブロック3の例のようにして作業を行う。以下、ブロックセットのブロックをブロック番号が小さいものから並べて、ブロック番号が大きいものから順に除いていってできるブロックの組合せをブロックの「サブセット」という。
【0085】
図11の例では、ブロック1、ブロック3、ブロック4がブロックセットとして選択されたときは、ブロック4を除いたサブセット(ブロック1、ブロック3)の作業は未だ行われていないので、次にブロック3も除いて、作業数1でブロック1を対象とするという作業が行われていないかを判断する。作業数1でブロック1を対象とするという作業は既に行われているので、この作業を繰り返し行うことはせず、この作業の結果既に作成されている候補クラスを利用して、その候補クラスに含まれるショートリードのブロック3をソート対象とするところから始める。
【0086】
図11の例において、ブロック2、ブロック3、ブロック4が選択された際には、作業数1でブロック2を対象とするという作業は未だ行われていないので、上記で説明したブロック1、ブロック2、ブロック3の例と同様にして最初から作業を行う。このようにして、重複した処理を回避して作業数を減らすことで、処理負担を軽減し、かつ高速に等価クラスを報告できる。
【0087】
以下、ソート対象指定部132にてあるブロックセットが選択された場合の、そのブロックセットに対する、距離計算対象の絞込みの手順について、更に詳細に説明する。図12は、ブロック1、ブロック2、ブロック3がブロックセットとして選択された場合のこのブロックセットに対する絞込み手順を説明するための図である。
【0088】
ソート対象指定部132は、最初に、作業数1(W=1)の作業のソート対象として、入力されたすべてのショートリードのブロック1を指定する。そして、この対象について上記のようにして、オフセット部133でオフセット部分配列をメモリ上に展開し、ソート部134で原部分配列とオフセット部分配列をソートしてグループ集合G1を生成する。候補クラス生成部135は、グループ集合を上からスキャンしていき、候補クラス作成条件を満たすグループが発見されたら、その候補クラス(候補クラスCC(1)1)をソート対象指定部132に報告する。本実施の形態では、候補クラスCC(1)1を報告する前に、候補クラスCC(1)1の最下列の位置E1を記憶しておく。
【0089】
ソート対象指定部132は、報告を受けた候補クラスCC(1)1に含まれるショートリードのブロック2をソート対象として指定する。そして、ソート対象指定部132が指定した対象、即ち候補クラスCC(1)1に含まれるショートリードのブロック2について、オフセット部132がオフセット部分配列を生成し、ソート部133がソートを行いグループ集合G2を生成する(ステップS11)。以下、候補クラス生成部135から報告された候補クラスに含まれるショートリードの、ソート対象指定部132が指定した作業ブロックについて、オフセット部132がオフセット部分配列をメモリ上に展開し、ソート部133がソートを行って配列パターンが同一である部分配列同士をグループ化することを当該候補クラスの当該ブロックを「展開・ソートする」と表現する。
【0090】
候補クラス生成部135は、候補クラスCC(1)1のブロック1を展開・ソートして生成されたグループ集合G2を上から順にスキャンしていき、候補クラスの生成条件を満たすグループが存在すると、そのグループを候補クラス(候補クラスCC(1-2)11)としてソート対象指定部132に報告する。本実施の形態では、候補クラスCC(1-2)11を報告する前に、候補クラスCC(1-2)11の最下列の位置E2を記憶しておく。そして、ソート対象指定部132は、候補クラスCC(1-2)11のブロック3をソート対象として指定する。この指定に基づいて、オフセット部132およびソート部133は、候補クラスCC(1-2)11のブロック3を展開・ソートしてグループ集合G3を生成する(ステップS12)。
【0091】
候補クラス生成部135は、展開・ソートにより得られたグループ集合G3を上から順にスキャンしていき、候補クラスの条件を満たすグループを発見すると、そのグループを候補クラスであると判断する。そして、作業数が最終作業数3であるため、この候補クラスCC(1-2-3)111を等価クラスとして、そのクラスに含まれるショートリードのインデクスの情報および各ショートリードのオフセット履歴の情報を類似判定部14に出力する(ステップS13)。なお、本実施の形態では、候補クラスCC(1-2-3)111を報告する前に、候補クラスCC(1-2-3)111の最下列の位置E3を記憶しておく。
【0092】
類似判定部14は、候補クラス生成部135から報告された等価クラス(候補クラスCC(1-2-3)111)について、当該等価クラスに含まれるショートリードのペアの編集距離を計算する。類似判定部14については後述する。等価クラス(候補クラスCC(1-2-3)111)について類似判定部14による計算が終了すると、グループ集合G3に戻って、再び上から順に候補クラスの条件を満たすグループがないかを探索する。この際、最上の列からスキャンして、すでに候補クラスであると判断されたグループを再び候補クラスと判断しないようにして次の候補クラスを探してもよいが、本実施の形態では、候補クラスCC(1-2-3)111を報告する前に、候補クラスCC(1-2-3)111の最下列の位置を記憶してあるので、その列から下に向けて次の候補クラスを探索する。
【0093】
グループ集合G3のすべての候補クラスが等価クラスとして報告されると、すなわち、グループ集合G3の最後まで探索が終了すると、候補クラス生成部135は、グループ集合G2に戻る(ステップS14)。このとき、既に展開・ソートされているグループ集合G3はメモリ上から解放される。候補クラス生成部135は、記憶されている前回の探索終了位置E2から候補クラスの探索を再開し、次の候補クラスCC(1-2)12を生成してソート対象指定部132に報告する。ソート対象指定部132は、候補クラスCC(1-2)12のブロック3をソート対象と指定し、オフセット部133およびソート部134は、決定された候補クラスCC(1-2)12のブロック3を展開・ソートする(ステップS15)。
【0094】
候補クラス生成部135は、候補クラスCC(1-2)12のブロック3を展開・ソートして得られたグループ集合G4を上から順にスキャンしていって、候補クラスCC(1-2-3)121を発見すると、これを等価クラスとして類似判定部14に報告する(ステップS16)。そして、類似判定部14にてこの等価クラスについての編集距離の計算が終わると、グループ集合G4を引き続きスキャンしていって、候補クラスCC(1-2-3)122を発見すると、これを等価クラスとして類似判定部14に報告する(ステップS17)。類似判定部14にてこの等価クラスについての編集距離の計算が終わると、候補クラス生成部135は、グループ集合G4を引き続きスキャンする。候補クラス生成部135は、列の最後までスキャンすると、グループ集合G4をメモリ上から解放し、グループ集合G2に戻って(ステップS18)、次の候補クラスを探索する。
【0095】
候補クラス生成部135は、以下同様にして、候補クラスCC(1-2)12のブロック3を展開・ソートしてグループ集合G5を生成し(ステップS19)、このグループ集合G5から候補クラスCC(1-2-3)122を等価クラスとして出力し(ステップS20)、その等価クラスについての距離の計算が終わると、グループ集合G5をメモリ上から解放して、グループ集合G2に戻る(ステップS21)。そして、候補クラス生成部135は、グループ集合G2に候補クラスがなくなると、グループ集合G2をメモリ上から解放して、グループ集合G1に戻って(ステップS22)、前回の探索終了位置E1から、次の候補クラスを探索する。
【0096】
候補クラス生成部135は、以下同様にして、候補クラスCC(1)2をステップS23、ステップS24、・・・と展開していく。以上は、候補クラス記憶部136に既に記憶されている候補クラスを利用しないで、距離計算の対象を絞り込んでいく場合の手順である。
【0097】
次に、図13を参照して、候補クラス記憶部136に既に記憶されている候補クラスを利用して、距離計算の対象を絞り込んでいく場合の手順を説明する。図13の例は、ブロック1、ブロック2、ブロック3ののブロックセットについて距離計算が終了した後に、ソート対象指定部132がブロック1、ブロック2、ブロック4を新たなブロックセットとして選択した(ブロック3、ブロック5をマスクした)ときの手順を示している。
【0098】
ソート対象指定部132は、ブロック1、ブロック2、ブロック4を選択すると、ブロック1、ブロック2のサブセットについて候補クラス記憶部136に候補クラスが記憶されていないかを判断する。この場合は、既に以前のブロックセット(ブロック1、ブロック2、ブロック3)についての絞込みでブロック1、ブロック2のサブセットについての候補クラスが候補クラス記憶部136に記憶されているので、ソート対象指定部132はこれを読み出して、その候補クラスに含まれるショートリードのブロック4をソート対象として決定する。
【0099】
オフセット部133およびソート部134は、最初の候補クラスCC(1-2)11に含まれるショートリードのブロック4を展開・ソートしてグループ集合G6を生成する(ステップS31)。候補クラス生成部135は、このグループ集合G6を上から順にスキャンして候補クラスを探索する。候補クラス生成部135は、候補クラスCC(1-2-4)111を発見すると、これを等価クラスとして類似判定部14に報告し(ステップS32)、類似判定部14は、この等価クラスについて編集距離を計算し、編集距離がd以内にあるすべてのショートリードのペアを示す情報(当該ペアのインデクス)を出力する。類似判定部14による類似判定が終了すると、候補クラス生成部135は、引き続き候補クラスを探索する。候補クラス生成部135は、候補クラスCC(1-2-4)112を発見すると、これを等価クラスとして類似判定部14に報告し(ステップS33)、類似判定部14は、この等価クラスについて編集距離を計算し、編集距離がd以内にあるすべてのショートリードのペアをのインデクスを出力する。
【0100】
類似判定部14による類似判定が終了すると、候補クラス生成部135は、グループ集合G6にて引き続き候補クラスを探索する。候補クラス生成部135は、グループ集合G6の最後までスキャンすると、グループ集合G6をメモリ上から解放して、その旨をソート対象部132に報告する。ソート対象部132は、この報告を受けて、候補クラス記憶部136から読み出した次の候補クラスCC(1-2)12に含まれるショートリードのブロック4をソート対象として指定する(ステップS34)。
【0101】
以下同様にして、オフセット部133およびソート部134が候補クラスCC(1-2)12に含まれるショートリードのブロック4を展開してグループ集合G7を生成し(ステップS35)、候補クラス生成部135がグループ集合G7から候補クラスCC(1-2-4)121を等価クラスとして類似判定部14に報告する(ステップS36)。類似判定部14がこの等価クラスに含まれるショートリードのペアの類似判定を終了すると、候補クラス生成部135は、次に候補クラスCC(1-2-4)122を等価クラスとして類似判定部14に報告する(ステップS37)。
【0102】
類似判定部14がこの等価クラスに含まれるショートリードのペアの類似判定を終了すると、候補クラス生成部135は候補クラスがなくなったことをソート対象指定部132に報告し、グループ集合G7をメモリ上から解放する(ステップS38)。ソート対象指定部132は、候補クラス記憶部136から読み出した次の候補クラスCC(1-2)13に含まれるショートリードのブロック4をソート対象として決定し、これを受けて、オフセット部133およびソート部134が候補クラスCC(1-2)13に含まれるショートリードのブロック4を展開・ソートしてグループ集合G8を生成し(ステップS39)、以下同様にしてステップS40、ステップS41、ステップS42、・・・と作業を進めていく。
【0103】
ソート対象指定部132は、候補クラス記憶部136から読み出したすべての候補クラスについての展開が終了して、すべての等価クラスが報告されると、新たなブロックセットを選択する。
【0104】
以上のように、本実施の形態の配列解析装置10は、候補クラス生成部135が候補クラスを生成するたびにそれをソート対象指定部132に報告し、ソート対象指定部132が報告を受けた候補クラスの次のブロックをソート対象として決定して、オフセット部133およびソート部134がそのソート対象を展開するという処理を繰り返して、最終作業数まできたときに候補クラス生成部135が候補クラスを等価クラスとして類似判定部14に報告して、類似判定部14で直ちに等価クラス内のショートリードのペアの編集距離を計算して、編集距離がd以内であるショートリードのペアのインデクスを出力する。仮に、これとは異なり、ブロックのすべての組について、すべての候補クラスを展開して、すべての等価クラスを求めた後に、各等価クラスについてそこに含まれるショートリードのペアの編集距離を計算するという手順にすると、求めた等価クラスの情報をすべて記憶装置に記憶しておかなければならなくなる。これに対して、本実施の形態の配列解析装置10のように、候補クラスを取得するごとに作業数を増やしていき、かつ等価クラスを取得するごとに類似判定を行うようにすると、記憶すべき情報を少なくでき、必要な記憶装置の容量も小さくできる。
【0105】
必要な記憶装置の容量を小さくするという観点では、候補クラス記憶部136に記憶する候補クラスの情報についても、不要となった時点で消去していく構成が望ましい。図11の例では、例えば、ブロック1、ブロック2、ブロック5を選択して距離計算対象の絞込みおよび距離計算を行った後には、作業数1の作業ブロックとしてブロック1を選択して、作業数2の作業ブロックとしてブロック2を選択して絞込みを行うことで得られた候補クラス(図13の左列の候補クラス)の情報は、消去することができる。また、ブロック1、ブロック4、ブロック5を選択して距離計算対象の絞込みおよび距離計算を行った後には、作業数1の作業としてブロック1を選択して絞込みを行うことで得られた候補クラス(図12の最左の列の候補クラス)の情報は、消去することができる。
【0106】
次に、類似判定部14について説明する。類似判定部14は、正準性判定部141と、編集距離計算部142と、類似ペア出力部143とを有している。上述のように、候補クラス生成部135からは、等価クラスの情報が類似判定部14に出力される。この等価クラスの情報には、等価クラスに含まれるショートリードのインデクスのほか、その等価クラスが生成されるまでにソート対象となったブロックの履歴、およびブロック履歴に含まれる各ブロックにおいていくつのオフセット値で他の部分配列と一致したのかを示すオフセット履歴の情報が含まれている。
【0107】
正準性判定部141は、候補クラス生成部135から等価クラスが入力されると、その等価クラスに含まれる各ショートリードのブロック履歴およびオフセット履歴に基づいて、正準性を判断する。仮に、等価クラスに含まれる複数のショートリードのすべてのペアについて編集距離を計算すると、すべての等価クラスについて編集距離を計算する過程で、同一のペアについて重複して編集距離を計算することになる。正準性判定部141による正準性の判定は、このように同一のペアについての重複した編集処理の計算を避けるために行われる判定である。
【0108】
正準性判定部141は、ショートリードのペアが下記の条件3および条件4を同時に満足するときに、当該ペアを正準(canonical)であると判定する。
条件3:当該ペアの少なくとも一方のショートリードのオフセット履歴に含まれるオフセットがすべてゼロである
条件4:ブロック履歴に含まれる任意のブロックのブロック番号よりも小さいブロック番号のブロックにおいて、一方のショートリードの部分配列または窓枠を最大オフセット内でオフセットさせたオフセット部分配列が、他方のショートリードの当該ブロックの部分配列と一致しない
【0109】
以下、例を挙げて正準性判定を説明する。図14(a)の左側は等価クラス生成部13から出力される等価クラスの情報の例である。等価クラスの情報には、当該等価クラスに含まれるショートリードのインデクスの情報、当該等価クラスが生成される際に選択されたブロック(ブロック履歴)の情報、当該等価クラスが生成される過程で生成された候補クラスにおける各ショートリードの部分配列のオフセット値(オフセット履歴)の情報が含まれる。
【0110】
正準性判定部141は、正準性判定のために、この等価クラスに含まれるショートリードのあらゆるペアについて、条件3および条件4を満たすかを判断する。このために、正準判定部14は、等価クラスの情報中のインデクスが示すショートリードをブロック分割部131から読み出す。図14(a)の右側は、実際に読み出したショートリードの配列パターンである。
【0111】
正準性判定部141は、この等価クラスに含まれるショートリードのペアについて、まず条件3を満たすか否かを判断する。S75とS88のみが、オフセット履歴中のオフセット値がすべてゼロであるという条件を満たすので、S75、S88のいずれかを含むペアのみが条件3を満たすことになる。これを表にまとめたのが図14(b)である。図14(b)に示すように、本来は、5つのショートリードのうちの2つのショートリードの組合せは、52=10通りあるが、そのうちの少なくとも一方にS75、S88を含むペアのみが条件3を満たす。
【0112】
正準性判定部141は、例えばS22とS75については、S75のオフセット履歴がすべてゼロであるため、条件3を満たすと判断する。図14(c)、(d)は、条件4の判断例を説明する図である。図14(c)は、S22とS75のペア、図14(d)はS56とS75のペアを示している。なお、いずれのペアも、オフセット履歴中のオフセット値がすべてゼロであるS75を含むので、条件3を満たす。以下、このオフセット履歴中のオフセット値がすべてゼロであるショートリードを「基準ショートリード」といい、他方のショートリードを「オフセットショートリード」という。
【0113】
条件4の判断ではまず、ブロック履歴に含まれる任意のブロックのブロック番号よりも小さいブロック番号のブロックに注目する。図14(c)、(d)の例では、ブロック1はブロック履歴に含まれるブロック2(およびブロック3、ブロック5)よりもブロック番号が小さいブロックである。よって、正準性判定部141は、オフセットショートリードについて、このブロック1の部分配列およびそこから窓枠を最大オフセット内でオフセットさせた複数のオフセット部分配列を生成し、それらの部分配列を、基準ショートリードのブロック1の部分配列と比較する。
【0114】
図14(c)の例では、オフセットショートリードであるS22においてブロック1から窓枠を最大オフセット1内でオフセットさせても、基準ショートリードであるS75のブロック1とは一致しないため、このS22とS75のペアは条件4を満たす。一方、図14(d)の例では、オフセットショートリードであるS56においてブロック1から窓枠を最大オフセット1内でオフセットさせると、オフセット0で、基準ショートリードであるS75のブロック1と一致するため、このS56とS75のペアは条件4を満たさない。
【0115】
距離計算部142は、条件3および条件4を何れも満たすと判断されたペアについて、その間の編集距離を計算する。なお、2つのショートリードの間の編集距離は、例えば、Daniel Jurafsky and James H.Martin: Speech and Language Processing, pp.74, Prentice Hall, 2009等に記載された公知の方法により計算する。類似ペア出力部143は、距離計算部142にて計算された編集距離がd以下であるときは、当該ショートリードのペアを類似するペアと判定して、当該ショートリードのペアを示すインデクスのペアを出力する(図1参照)。
【0116】
以上のように、条件3および条件4によって正準性を判断し、正準であると判断したペアのみについて編集距離を計算することで、計算コストの高い編集距離を同一のペアについて重複して行うことによる計算コストの増加を避けることができる。
【0117】
本実施の形態の配列解析装置10によれば、等価クラス生成部13から等価クラスとして出力されたショートリードの集合によってできるすべてのペアの中には、編集距離がd以下であるペアがもれなく含まれており、かつ、上述のように同一のペアについての重複した編集距離の計算を回避できるので、計算コストを削減した上で、編集距離がd以内にあるすべてのペアを重複なく見つけ出すことができる。
【0118】
次に、本実施の形態の配列解析装置10によって実行される配列解析方法を図15を参照して説明する。図15は、配列解析装置10の処理フロー図である。まず、配列解析装置10は配列入力部11にて、シーケンサから出力された複数のショートリードを入力し、かつ条件入力部12にて配列解析の条件、即ち編集距離d、およびブロック分割数を(d+n)によって決定するための数値nを入力する(ステップS51)。
【0119】
次に、入力されたすべてのショートリードは、等価クラス生成部13のブロック分割部に記憶される。ブロック分割部131は、すべてのショートリードを(d+n)個のブロックに分割する(ステップS52)。ソート対象指定部132は、(d+n)個のブロックの中からn個のブロックを選択する(ステップS53)。
【0120】
続いて、ソート対象指定部132は、選択されたn個のブロックのうちのサブセットの候補クラスが候補クラス記憶部136に記憶されているかを判断する(ステップS54)。選択されたn個のブロックについてサブセットの候補クラスがない場合は(ステップS54でNO)、作業数Wを1とし、選択されたブロック中のブロック番号が最小のブロックを作業ブロックbkとし、かつ、候補クラスをすべてのショートリードとする(ステップS55)。
【0121】
続いて、オフセット部133およびソート部134にて候補クラスを展開して、グループ集合を生成する(ステップS56)。なお、ステップS54でサブセットの候補クラスがある場合には(ステップS54でYES)、当該サブセットの候補クラスを候補クラス記憶部136から読み出して、読み出した候補クラスに対応する作業数Wおよび作業ブロックbkを設定して、ステップS56に進む。この場合は、ステップS56では、候補クラス記憶部136から読み出した候補クラスを展開して、グループ集合を生成する。
【0122】
ステップS56の後に、生成したグループ集合を上から順にスキャンして(ステップS57)、同一の配列パターンのグループについて、候補クラスの条件(条件1および条件2)を満たすか否かを判断する(ステップS58)。そのグループが候補クラスの条件を満たす場合は(ステップS58でYES)、作業数Wがnであるか否かを判断する(ステップS59)。作業数がnでない場合には、候補クラス生成部135は、発見したグループを候補クラスとしてソート対象指定部132に報告し、かつ候補クラス記憶部136に保存する(ステップS60)。
【0123】
候補クラス記憶部136は、その候補クラスに含まれるショートリードのインデクスおよびオフセット履歴をブロック履歴に関連付けて記憶する。ソート対象指定部132は、候補クラス記憶部136から候補クラスの報告を受けると、作業数Wをインクリメントし、作業ブロックを、選択されている複数のブロックの次のブロックに変更し(ステップS61)、ステップS56に戻って、候補クラス生成部135から報告された候補クラスを展開して部分配列の列を生成する。
【0124】
このようにして、ステップS56〜S59の処理を繰り返して、ステップS59にて作業数Wがnになると、候補クラス生成部135は、その直前に生成した候補クラスを等価クラスとして類似判定部14に報告する(ステップS62)。類似判定部14は、等価クラスについて類似判定を行い、等価クラスに含まれるペアについて編集距離を計算する(ステップS63)。そして、候補クラス生成部135は、グループ集合に未スキャン部分があるか否かを判断する(ステップS64)。また、ステップS57にて、グループ集合をスキャンして同一の配列パターンのグループが候補クラスの条件を満たさなかった場合にも(ステップS58でNO)、ステップS64にて未スキャン部分があるか否かを判断する。候補クラス生成部135は、未スキャン部分がある場合には(ステップS64でNO)、ステップS57に戻って当該未スキャン部分をスキャンする。
【0125】
一方、未スキャン部分がなくなったとき、即ち、グループ集合をすべてスキャンし終わったときは、作業数Wをデクリメントして、ブロック番号をブロックセット中の1つ前の番号に変更し(ステップS65)、作業数Wがゼロであるかを判断する(ステップS66)。作業数Wがまだゼロになっていない場合には(ステップS66でNO)、ステップS64に戻って、ステップS65でデクリメントされた作業数で展開されたグループ集合に未スキャン部分があるか否かを判断する。
【0126】
ステップS66で作業数Wがゼロとなっているときは、そのとき選択されているブロックセットについては既にすべての等価クラスが報告されているということであるため、このブロックセットについて展開されたグループ集合をすべてクリアし(ステップS67)、ソート対象指定部132は、未選択のブロックの組合せがまだあるかを判断する(ステップS68)。未選択のブロックセットが未だある場合には(ステップS68でYES)、ステップS53に戻って新たなn個のブロックを選択する。一方、未選択のブロックセットがなくなった場合は(ステップS68でNO)、処理を終了する。
【0127】
以下、上記のフローが図12〜図13の例で実行されるときに、実際にどのようにフローが進むかを説明する。まず、ステップS51にて複数のショートリードおよび条件(d=2、n=3)が入力されて、ステップS52にて複数のショートリードがブロックに分割される。次に、ステップS53にてソート対象指定部132は、n個のブロックを選択してブロックセットとするが、このとき、ブロック1、ブロック2、ブロック3が選択される。
【0128】
次に、ステップS54にて、ブロック1、ブロック2、ブロック3のブロックセットについては、サブセットの候補クラスがまだないので、ステップS55に進み、作業数Wが1とされ、ブロック番号bkは、ブロック1、ブロック2、ブロック3のうちの一番小さいブロック1とされ、候補クラスはすべてのショートリードとされる。そして、ステップS56にて、候補クラスであるすべてのショートリードのブロック1が展開されて、グループ集合G1が生成される。そして、ステップS57にて、このグループ集合G1が上からスキャンされて、ステップS58では、候補クラスCC(1)1より上のグループについては、候補クラスの要件を満たさないため、ステップS64に進み、この時点ではまだグループ集合G1をすべてスキャンしていないので、ステップS57に戻る。
【0129】
このようにして、ステップS57にて、候補クラスCC(1)1が候補クラスの要件を満たすので、ステップS59にて作業数Wがn(=3)であるか否かを判断する。作業数Wは1であって、W=nは満たさないため、ステップS60に進み、この候補クラスCC(1)1が報告される。そして、ステップS61にて、作業数Wがインクリメントされて、W=2となり、作業ブロックが、選択されたブロック1の次のブロック、即ちブロック2とされる。その後、ステップS56に戻って、この候補クラスCC(1)1を展開して(図12のステップS11に対応する)、グループ集合G2を生成する。
【0130】
次に、ステップS57にて、グループ集合G2を上からスキャンして、ステップS58で候補クラスCC(1-2)11が見つかると、ステップS59では、W=2であって、W=n(=3)ではないので、ステップS60でこの候補クラスCC(1-2)11が報告される。そして、ステップS61で、作業数Wが2から3にインクリメントされ、作業ブロックは次のブロック3となる。
【0131】
そして、再びステップS56に戻って、候補クラスCC(1-2)11のブロック3が展開されてグループ集合G3が生成され(図12のステップS12に対応する)、ステップS57で、グループ集合G3がスキャンされて、ステップS58で、候補クラスCC(1-2-3)111が生成される。次のステップS59では、作業数がW=3であるので、W=n(=3)となり、ステップS62に移行して、候補クラス生成部135は、この候補クラスCC(1-2-3)111を等価クラスとして類似判定部14に報告する(図12のステップS13に対応する)。ステップS63では、類似判定部14がこの等価クラス候補クラスCC(1-2-3)111について、ショートリードのペアの編集距離を計算し、編集距離がd(=2)以内にあるショートリードのペアのインデクスを出力する。
【0132】
その後、グループ集合G3に未スキャン部分があるか判断される。未スキャン部分があるので、ステップS57に戻って、候補クラスCC(1-2-3)111より下の部分について、ステップS57、ステップS58、ステップS64を繰り返す。グループ集合G3の最後まで来ると、ステップS64にて、グループ集合G3に未スキャン部分がなくなるので、ステップS65で作業数を3から2にデクリメントして(図12のステップS14に対応する)、かつ作業ブロックをブロック3からブロック2に変更する。ステップS66では、W=2であって、W=0ではないので、ステップS64に戻って、グループ集合G2に未スキャン部分があるか否かが判断される。図12のステップS12の時点では、候補クラスCC(1-2)11までしかスキャンしていなかったので、ステップS64では、グループ集合G2に未スキャン部分があるということになり、ステップS57に戻る。
【0133】
そして、ステップS57で、グループ集合G2について、候補クラスCC(1-2)11の下からスキャンを再開する。ステップS58で、候補クラスCC(1-2)12を生成すると、ステップS59でW=nであるかを判断し、W=2であってW=d(=3)ではないので、ステップS60に移行して、この候補クラスCC(1-2)12を報告し、ステップS61で作業数を2から3にインクリメントして作業ブロックをブロック2から次のブロック3にして、ステップS56で、候補クラスCC(1-2)12のブロック3を展開して(図12のステップS15に対応する)、グループ集合G4を生成する。
【0134】
次に、ステップS57で、グループ集合G4がスキャンされて、ステップS58で、候補クラスCC(1-2-3)121が生成される。次のステップS59では、作業数がW=3であるので、W=n(=3)となり、ステップS62に移行して、候補クラス生成部135は、この候補クラスCC(1-2-3)121を等価クラスとして類似判定部14に報告する(図12のステップS16に対応する)。ステップS63では、類似判定部14がこの等価クラスCC(1-2-3)121について、ショートリードのペアの編集距離を計算し、編集距離がd(=2)以内にあるショートリードのペアのインデクスを出力する。
【0135】
その後のステップS62では、グループ集合G4にまだ未スキャン部分があるので、ステップS57に戻って、その未スキャン部分がスキャンされる。そして、ステップS58で、候補クラスCC(1-2-3)122が生成される。次のステップS59では、作業数がW=3であるので、W=n(=3)となり、ステップS62に移行して、候補クラス生成部135は、この候補クラスCC(1-2-3)122を等価クラスとして類似判定部14に報告する(図12のステップS17に対応する)。ステップS63では、類似判定部14がこの等価クラスCC(1-2-3)122について、ショートリードのペアの編集距離を計算し、編集距離がd(=2)以内にあるショートリードのペアのインデクスを出力する。
【0136】
その後、ステップS64からステップS57に戻って、ステップS57、ステップS58、ステップS64を繰り返し、グループ集合G4の最後まで来ると、ステップS65にて、作業数が3から2にデクリメントされて、作業ブロックが1つ前のブロック2に変更されて、グループ集合G2に戻る。作業数はまだゼロではないので、ステップS64に戻って、さらに、グループ集合G2に未スキャン部分があるのでステップS57に戻って、グループ集合G2の候補クラスCC(1-2)12の下からスキャンを開始する(図12のステップS18に対応する)。
【0137】
ステップS57、ステップS58、ステップS64を繰り返し、候補クラスCC(1-2)13まで来ると、ステップS58からステップS59に移行して、作業数2はW=d(=3)を満たさないので、ステップS60でこの候補クラスCC(1-2)13を報告して、作業数を2から3にインクリメントし、作業ブロックをブロック3とし、この候補クラスCC(1-2)13のブロック3を展開して(図12のステップS19に対応する)、グループ集合G5を生成する。
【0138】
ステップS57でグループ集合G5がスキャンされて、ステップS58で、候補クラスCC(1-2-3)131が生成される。次のステップS59では、作業数がW=3であるので、W=n(=3)となり、ステップS62に移行して、候補クラス生成部135は、この候補クラスCC(1-2-3)131を等価クラスとして類似判定部14に報告する(図12のステップS20に対応する)。ステップS63では、類似判定部14がこの等価クラス候補クラスCC(1-2-3)131について、ショートリードのペアの編集距離を計算し、編集距離がd(=2)以内にあるショートリードのペアのインデクスを出力する。
【0139】
その後、グループ集合G5に未スキャン部分があるか判断される。未スキャン部分があるので、ステップS57に戻って、候補クラスCC(1-2-3)131より下の部分について、ステップS57、ステップS58、ステップS64を繰り返す。グループ集合G5の最後まで来ると、ステップS64にて、グループ集合G5に未スキャン部分がなくなるので、ステップS65で作業数を3から2にデクリメントして、かつ作業ブロックをブロック3からブロック2に変更して(図12のステップS21に対応する)、ステップS66では、W=2であって、W=0ではないので、ステップS64に戻って、グループ集合G2に未スキャン部分があるか否かが判断される。図12のステップS21の時点では、候補クラスCC(1-2)13までしかスキャンしていなかったので、ステップS64では、未スキャン部分があるということになり、ステップS57に戻る。
【0140】
その後は、ステップS57、ステップS58、ステップS64を繰り返して、候補クラスが見つからずにグループ集合G2の最後の行まで来ると、ステップS64でNOとなり、ステップS65で、作業数を2から1にデクリメントして、作業ブロックも1つ前のブロック1に変更する。
【0141】
以下同様にして処理を進め、スキャンの対象がグループ集合G1の最終行までくると、ステップS64にてNOとなり、作業数1がデクリメントされて0になる。そうすると、ステップS67にて、これまでに展開したグループ集合をすべてクリアする。ステップS68では、まだブロック1、ブロック2、ブロック3という最初のブロックセットを処理しただけであるので、ステップS53に移行して、ソート対象指定部132は、ブロック1、ブロック2、ブロック4を新たなブロックセットとして選択する。
【0142】
そして、ステップS54では、サブセットの候補クラスが既にあるか否かが判断されるが、ここでは、いま選択されているブロックセット(ブロック1、ブロック2、ブロック4)のうちのブロック1、ブロック2というサブセットの候補クラスは、前のブロックセット(ブロック1、ブロック2、ブロック3)について処理を行った際に既に生成されているので、ステップS69に移行して、これが読み出される(図13の左列)。そして、この読み出した候補クラスCC(1-2)11、CC(1-2)12、CC(1-2)13、CC(1-2)21、CC(1-2)22、CC(1-2)31、CC(1-2)32、・・・に合わせて、作業数は3、作業ブロックはブロック4と設定される。
【0143】
そして、ステップS56で、読み出された最初の候補クラスCC(1-2)11に含まれるショートリードのブロック4を展開する(図13のステップS31に対応する)。以下同様にして、グループ集合G6から等価クラスCC(1-2-3)111、CC(1-2-3)112を順に報告した後、グループ集合G6の最後までスキャンして、ステップS64でNOとなった場合には、作業数を3から2にデクリメントして、ステップS66を経てステップS64に戻る。そして、ステップS64では、グループ集合(ここでは、候補クラス記憶部136より読み出した複数の候補クラス)に未スキャン部分(ここでは、まだ展開していない候補クラス)があるか否かを判断する。
【0144】
そして、候補クラス記憶部136より読み出した複数の候補クラスの中に、まだ展開していない候補クラスがあるときは、ステップS57に戻り、候補クラス記憶部136より読み出した候補クラスをスキャンする。ステップS58では、候補クラスCC(1-2)12についてYESに進み、Wはnまで達していないので、ステップS60で、ソート対象指定部132は、この候補クラスCC(1-2)12を次の候補クラスとして決定する。そして、ステップS61で、作業数が2から3にインクリメントされて、作業ブロックがブロック4とされて、次のステップS56では、候補クラスCC(1-2)12に含まれるショートリードのブロック4が展開される(図13のステップS35に対応する)。
【0145】
以下同様にして、図13のステップS36、S37、・・・、S42、・・・と処理が進み、候補クラス記憶部136から読み出された複数の候補クラスがすべて展開されて、ステップS63で、このブロックセットについての最後の等価クラスについて類似性が判定され、その後にステップS64でNOとなると、作業数が3から2にデクリメントされて、ステップS64に戻って、未展開の候補クラスがあるか否かが判断される。すべての候補クラスについて展開済みであるので、ステップS64ではNOとなり、作業数が更に2から1へデクリメントされる。W=0ではないのでステップS64に戻るが、このブロックセットについては、作業数1ではもともと作業が行われていないため、S64はNOになる。そして、ステップS65で作業数がさらに1から0にデクリメントされるので、ステップS66でYESとなり、ステップS67、S68を経て、ステップS53に戻って、次のブロックセットである、ブロック1、ブロック2、ブロック5が選択される。
【0146】
以上のように、本発明の実施の形態の配列解析装置によれば、大量のショートリードの中から類似のペアを完全に、かつ重複なく抽出することができる。また、上記の第1の実施の形態では、選択された複数のブロックをブロック番号が若い方から順に展開していくに当たって、あるブロックを展開して得られたすべての候補クラスを展開してから次のブロックを展開するのではなく、あるブロックを展開して得られたグループ集合にて候補クラスを発見するごとに当該候補クラスについて次のブロックを展開するという処理を繰り返し、処理の済んだグループ集合を逐次メモリ上から解放することで、展開したグループ集合を記憶するためのメモリの容量を削減することができる。
【0147】
第1の実施の形態では、まずブロックセットを選択して、当該選択したブロックセットについてのすべての処理が終了した後に新たなブロックセットを選択して処理を行った。このために、あるブロックセットで処理した際の候補クラスを候補クラス記憶部136に記憶し、新たに選択されたブロックセットで利用できる候補セットを候補クラス記憶部136から読み出す構成を採用した。このような構成は、候補クラス記憶部136から既に作成された候補クラスを読み出すのに時間がかかってしまう。一方で、既に作成されて後のブロックセットで再利用する候補クラスをすべてメモリに記憶しておくとすると、きわめて大きなメモリ容量が必要になる。
【0148】
そこで、以下では、第1の実施の形態と比較して、メモリ容量を増大させることなく、候補クラス記憶部136のようなハードディスクへの書き込みおよび読み出しを不要として処理速度を向上できる第2の実施の形態を説明する。第2の実施の形態の配列解析装置は、候補クラス記憶部136がないことを除き、図5に示した配列解析装置と同じである。
【0149】
図16〜図19は、本発明の第2の実施の形態における絞込み手順を説明する図である。なお、図16〜図19において、グループ集合のうち、候補クラスを形成しない部分配列のグループを省略している。また、図16〜図19では、編集距離d=1、数値n=3、分割数b(=d+n)=4の例を説明している。
【0150】
本実施の形態では、第1の実施の形態と同様に、ソート対象指定部132は、図11のツリー構造に従ってブロック番号の若いブロックから順番にブロックを選択する。但し、生成された候補クラスを記憶しておくことを避けるため、生成された候補クラスについて、当該候補クラスを含むグループ集合の他の候補クラスを展開・ソートする処理を行う前に、当該候補クラスを利用する処理をすべて行うという手順を採用する。
【0151】
具体的には、図16に示すように、ソート対象指定部132は、まず、作業開始ブロックをブロック1として、すべてのショートリードのブロック1を展開・ソートする。このときに、グループ集合G9には、候補クラス1−1、候補クラス1−2、候補クラス1−3、候補クラス1−4が含まれているとする。候補クラス生成部135は、グループ集合G9を上から順にスキャンしていって、候補クラス1−1を発見すると、ソート対象指定部132に報告する。ソート対象指定132は、候補クラス1−1に含まれるショートリードについて、ブロック2をソート対象として指定し、オフセット部133およびソート部134は、指定された対象を展開・ソートしてグループ集合G10を生成する。
【0152】
候補クラス生成部135は、このグループ集合G10をスキャンしていって候補クラス2−1を発見すると、その旨をソート対象指定部132に報告する。ソート対象指定部132は、候補クラス2−1に含まれるショートリードのブロック3をソート対象に指定する。オフセット部133およびソート部134は、指定された対象を展開・ソートしてグループ集合G11を生成する。候補クラス生成部135は、このグループ集合G11をスキャンしていって候補クラス3−1を発見すると、これを等価クラスとして類似判定部14に出力し、候補クラス3−1をメモリ上から解放する。候補クラス生成部135は、続いて、候補クラス3−2を発見すると、これを等価クラスとして類似判定部14に出力し、候補クラス3−2をメモリ上から解放する。ここまでの処理の手順は、第1の実施の形態と同様である。
【0153】
ソート対象指定部132は、次に、展開済みの候補クラスについて、まだ処理を行っていないブロックの組合せで処理を行う。図16の例では、ブロック1の候補クラス1−1をブロック2で展開・ソートして得られた候補クラス2−1については、ブロック1、2、4というブロックセットでは処理をしていないため、候補クラス2−1に含まれるショートリードについて、ブロック4を展開・ソートする。こうしてできたグループ集合G12について、等価クラス4−1、4−2、4−3を順に出力するとともに、出力し終わった等価クラスを順に消去していく。そうしてグループ集合G12について処理が終わると(このときグループ集合G12はすべてメモリ上から解放されている)、ソート済み部分配列集合G10の候補クラス2−1はすべて処理が終わったので、図17に示すように、この部分をメモリ上から解放する。
【0154】
その後、グループ集合G10を続けてスキャンし、候補クラス2−2を発見すると、この候補クラス2−2に含まれるショートリードについて、ブロック3を展開・ソートしてグループ集合G13を得て、等価クラス3−1、3−2を出力してグループ集合G13をメモリ上から解放し、その後に候補クラス2−2に含まれるショートリードについて、ブロック4を展開・ソートしてグループ集合G14を得て、等価クラス4−1、4−2、4−3、4−4を出力してグループ集合G14をメモリ上から解放する。
【0155】
ここまで終了すると、グループ集合G10の候補クラス2−2はもう使用しないため、メモリ上から解放する。その後、同様にして、候補クラス2−3に含まれるショートリードのブロック3、ブロック4について、順に展開・ソートして、等価クラスを出力する。このようにして、グループ集合G10に含まれるすべての候補クラスの処理が終了すると(この時点でグループ集合G10はすべてメモリ上から解放されている)、グループ集合G9の候補クラス1−1に含まれるショートリードのブロック3、ブロック4について、順に展開・ソートして、等価クラスを出力する。この時点で、グループ集合の候補クラス1−1のついてはすべての処理が終了したことになるので、図18に示すように、グループ集合G9の候補クラス1−1をメモリ上から解放し、候補クラス1−1の終わりからグループ集合G9をスキャンする。候補クラス1−2を発見すると、上記と同様にして、展開・ソート、候補クラスの生成を繰り返して、等価クラスを出力する。
【0156】
以下同様にして処理を進め、候補クラス1−4についてすべての処理が終了し、即ち、グループ集合G9についてすべての処理が終了すると、この時点でグループ集合G9はすべてメモリ上から解放されていることになる。このように、作業開始ブロックについてグループ集合がすべて消去されると、ソート対象指定部132は、図19に示すように、作業開始ブロックを次のブロック2として、上記と同様の処理を行う。
【0157】
このように、第2の実施の形態では、グループ集合を順にスキャンしていき、発見された候補クラスを利用するすべてのパターンのブロックセットについて等価クラスを求め、すべての等価クラスを出力したら当該候補クラスをメモリ上から解放して、そのグループ集合内の次の候補クラスを探すという処理を行う。これにより、第1の実施の形態のように、候補クラスを後のブロックセットで利用するために保存しておく必要がなくなり、候補クラス記憶部のようなハードディスクへの書き込みおよび読み出しを行う必要はなくなると同時に、必要なメモリ容量も第1の実施の形態と同程度に小さく抑えることができる。
【0158】
図20は、第2の実施の形態の配列解析装置の処理フロー図である。最初に複数のショートリードおよび条件(編集距離dおよび数値n)を入力し(ステップS71)、入力された複数のショートリードの各々をブロックに分割する(ステップS72)処理は、第1の実施の形態のステップS51、ステップS52(図15参照)と同じである。ブロック分割部131にブロックに分割された状態のショートリードが準備されると、次に、作業数Wを1とし、作業ブロックbkをブロック1とし、作業開始ブロックSを1とする(ステップS73)。そして、絞込みを行う対象となる候補クラス(以下、「絞込み対象候補クラス」という)Aを入力されたすべてのショートリードとする(ステップS74)。
【0159】
次に、絞込み対象候補クラスAの作業ブロックbkを展開・ソートして、グループ集合G(bk)を生成して記憶し、作業数Wの作業済みブロックであって、図16〜図19において最も右にあるブロックP(W)を作業ブロックbkとする(ステップS75)。そして、グループ集合G(bk)を上からスキャンして、候補クラスの条件を満たす最初のグループを選択してそのグループを絞り込み対象候補クラスAとし、かつ、作業数Wの絞込み対象候補クラスQ(W)をその絞り込み対象候補クラスAとする(ステップS76)。
【0160】
次に、作業数Wがnになっているかを確認し(ステップS77)、まだnに達していない場合は(ステップS77でNO)、作業数Wをインクリメントして(ステップS78)、作業ブロックbkをインクリメントさせた値(bk+1)がd+W以下であるか否かを判断する(ステップS79)。bk+1≦d+Wであるときは(ステップS79にてYES)、作業ブロックbkをインクリメントさせて(ステップS80)、ステップS75に戻る。
【0161】
ステップS75〜S80のループを繰り返して、ステップS77にて作業数Wがnに達すると(ステップS77にてYES)、そのときの絞込み対象候補クラスAを等価クラスとして(ステップS81)、当該等価クラスAについて、類似判定部14で類似判定を行う(ステップS82)。そして、グループ集合G(bk)に未スキャン部分があるか否かを判断し(ステップS83)、ある場合は(ステップS83にてYES)、ステップS76に戻ってスキャンを続ける。ステップS76、ステップS77、ステップS81、ステップS82、ステップS83を繰り返し、グループ集合G(bk)をすべてスキャンして、未スキャン部分がなくなった場合には(ステップS83にてNO)、当該グループ集合G(bk)をメモリ上から解放し、作業数Wをデクリメントする(ステップS84)。
【0162】
その後、作業数Wが0になっていないかを確認し(ステップS85)、Wが0になっていない場合は(ステップS85にてNO)、絞込み対象候補クラスAを作業数Wの絞込み対象候補クラスとし(ステップS86)、ステップS78に移行する。ステップS75〜S86を繰り返し、ステップS79にて、bk+1>d+Wとなった場合は(ステップS79でNO)、作業数Wをデクリメントし、作業ブロックbkをそのデクリメントされた作業数の作業済みの最右のブロックとして(ステップS87)、ステップS83に移行する。
【0163】
上記のようなステップS75〜S87を繰り返し、ステップS85にて作業数Wが0になった場合には(ステップS85でYES)、作業開始ブロックSが編集距離d以下であるかを判断し(ステップS88)、作業開始ブロックSが編集距離dを越えていない場合には(ステップS88でYES)、作業開始ブロックSをインクリメントして、作業数Wを1とし(ステップS89)、ステップS74に戻る。
【0164】
上記のようにして、ステップS74〜S89を繰り返し、ステップS88にて作業開始ブロックSが編集距離dより大きくなった場合には(ステップS88でNO)、処理を終了する。
【0165】
以下、上記のフローが図16〜図19の例で実行されるときに、実際にどのようにフローが進むかを説明する。まず、ステップS73で作業数Wを1とし、作業ブロックbkを1とし、作業開始ブロックSが1とされる。そして、ステップS74で、絞込み対象候補クラスAが入力されたすべてのショートリードとされる。次に、ステップS75で入力されたすべてのショートリードのブロック1が展開・ソートされてグループ集合G(1)(=G9)が生成され、作業数1の作業済み最右のブロックP(1)がブロック1とされる。
【0166】
次に、ステップS76にて、グループ集合G(1)(=G9)がスキャンされ、選択した候補クラス1−1が絞込み対象候補クラスAとされ、作業数1の絞込み対象候補クラスQ(1)が候補クラス1−1とされる。ステップS77では、作業数W(=1)はn(=3)ではないので、ステップS78に移行して、作業数Wが1から2にインクリメントされる。次に、ステップS79では、bk+1(=2)は、d+W(=3)以下であるので、ステップS80に移行して作業ブロックがブロック1からブロック2とされる。
【0167】
次に、ステップS75にて、いま絞込み対象候補クラスAとされている候補クラス1−1をブロック2について展開・ソートして、グループG(2)(=G10)が生成され、作業数2の作業済み最右のブロックP(2)をブロック2とする。続いて、ステップS76では、グループ集合G(2)(=G10)をスキャンして、候補クラス2−1を選択して、これを絞込み対象候補クラスAとし、作業数2の絞込み対象候補クラスQ(2)を候補クラス2−1とする。そして、ステップS77では、作業数W(=2)はn(=3)ではないので、ステップS78に移行して、作業数Wが2から3にインクリメントされる。次に、ステップS79では、bk+1(=3)は、d+W(=4)以下であるので、ステップS80に移行して作業ブロックがブロック2からブロック3とされる。
【0168】
以下、同様にして、ステップS75にて、いま絞込み対象候補クラスAとされている候補クラス2−1をブロック3について展開・ソートして、グループ集合G(3)(=G11)が生成され、作業数3の作業済み最右のブロックP(3)をブロック3とする。続いて、ステップS76では、グループ集合G(3)(=G11)をスキャンして、候補クラス3−1を選択して、これを絞込み対象候補クラスAとし、作業数3の絞込み対象候補クラスQ(3)を候補クラス3−1とする。そして、ステップS77では、W(=3)=n(=n)なので、ステップS81に移行して、絞込み対象候補クラス3−1を等価クラスとして出力し、ステップS82にて、類似判定部14にて等価クラス3−1について類似判定を行う。
【0169】
次に、ステップS83では、グループ集合G(3)(=G11)にまだ未スキャン部分があるので、ステップS76に戻る。ステップS76では、候補クラス3−2を選択し、それを絞込み対象候補クラスAとし、作業数3の絞込み対象候補クラスQ(3)を候補クラス3−2とする。ステップS77では、作業数W(=3)=n(=3)であるので、ステップS81に移行して、この候補クラス3−2を等価クラスとして出力し、ステップS82で等価クラス3−2について類似判定を行う。
【0170】
次に、ステップS83に移行すると、グループ集合G(3)(=G11)には未スキャン部分はないので、ステップS84に移行して、グループ集合G(3)(=G11)はメモリ上から解放され、作業数Wは3から2にデクリメントされる。そして、ステップS85では、作業数W(=2)は0ではないので、ステップS86に移行して、絞込み対象候補クラスAが、絞込み対象候補クラスQ(2)、即ち候補クラス2−1とされる。続いて、ステップS78にて、作業数Wが2から3にインクリメントされ、ステップS79では、bk+1(=4)は、d+W(=4)以下であるので、ステップS80に移行して、作業ブロックbkが3から4にインクリメントされる。
【0171】
ステップS75に戻り、いま絞込み対象候補クラスAとされている候補クラス2−1がブロック4について展開・ソートされ、グループ集合G(4)(=G12)が生成され、作業数3の作業済み最右のブロック(3)がブロック4とされる。続いて、ステップS76にて、グループ集合G(4)(=G12)がスキャンされて、最初の候補クラス4−1が選択され、これが絞込み対象候補クラスAとされ、作業数3の絞込み対象候補クラスQ(3)が候補クラス4−1とされる。ステップS77では、W(=3)=n(=3)であるので、ステップS81に移行して、候補クラス4−1が等価クラスとして出力され、ステップS82で、等価クラス4−1について類似判定が行われる。次に、ステップS83では、グループ集合G(4)(=G12)にはまだ未スキャン部分が残っているので、ステップS76に戻る。
【0172】
その後は、ステップS76、S77、S81、S82、S83を繰り返すことで、等価クラス4−2、等価クラス4−3を順に出力し、その都度類似判定を行う。等価クラス4−3まで類比判定が終わると、ステップS83では、グループ集合G(4)(=G12)に未スキャン部分がなくなり、ステップS84に移行して、グループ集合G(4)(=G12)がメモリ上から解放され、作業数Wは、3から2にデクリメントされる。ステップS85では、作業数W(=2)は0ではないので、ステップS86に移行して、作業数2の絞込み対象候補クラスQ(2)である2−1がひとまず絞込み対象候補クラスAとされ、ステップS78では、作業数Wが2から3にインクリメントされる。
【0173】
そして、ステップS79では、bk+1(=5)はd+W(=4)以下ではなくなるので、ステップS87に移行して、作業数Wが3から2にデクリメントされ、作業ブロックbkは、作業数2の作業済み最右のブロックP(2)、即ちブロック2とされ、ステップS83で、グループ集合G(2)(=G10)に未スキャン部分があるか判断される。グループ集合G(2)(=G10)には未スキャン部分があるので、ステップS76に戻り、グループ集合G(2)(=G10)のスキャンを再開して、候補クラス2−2を選択し、これを絞込み対象候補クラスAとして、作業数2の絞込み対象候補クラスQ(2)を候補クラス2−2とする。
【0174】
ステップS77では、作業数W(=2)はn(=3)と等しくないないので、ステップS78に移行して、作業数Wが2から3にインクリメントされる。ステップS79では、bk+1(=3)は、d+W(=4)以下であるので、ステップS80に移行して、作業ブロックbkがブロック2からブロック3に変更されて、ステップS75に戻る。ステップS76では、絞込み対象候補クラスAである候補クラス2−2が、作業ブロックbk(=ブロック3)について展開・ソートされて、グループ集合G(3)(=G13)が生成される。そして、このグループ集合G13について、ステップS76、S77、S81、S82、S83が繰り返されて、グループ集合G13の等価クラス3−1、3−2が出力された後に、ステップS83でグループ集合G13に未スキャン部分がないと判断されて、ステップS84に移行して、グループ集合G13がメモリ上から解放されて、作業数Wが3から2にデクリメントされる。
【0175】
ステップS85では、W(=2)は0ではないので、絞込み対象候補クラスAは作業数2の絞込み対象候補クラスQ(2)、即ち候補クラス2−2とされる。その後、ステップS78で、作業数Wが2から3にインクリメントされて、ステップS79では、bk+1(=4)は、d+W(=4)以下であるので、ステップS80に移行して作業ブロックbkがブロック3からブロック4に変更されて、ステップS75に戻る。ステップS75では、絞込み対象候補クラスAである候補クラス2−2が作業ブロックbk(=ブロック4)について展開・ソートされて、グループ集合G(4)(=G14)が生成され、作業数W(=3)の作業済み最右のブロックP(3)が作業ブロックbk(=ブロック4)とされる。
【0176】
その後、ステップS76、S77、S81、S82、S83を繰り返すことで、グループ集合G14から等価クラス4−1、4−2、4−3、4−4が出力された後、ステップS83にてグループ集合G14に未スキャン部分がなくなったと判断されると、ステップS84で、グループ集合G14がメモリ上から解放され、作業数Wは3から2にデクリメントされる。次に、ステップS85では、作業数W(=2)は0ではないので、ステップS86で絞込み対象候補クラスAは、一旦、作業数2の絞込み対象候補クラスQ(2)即ち候補クラス2−2とされる。そして、ステップS78で作業数Wが2から3にインクリメントされ、ステップS79では、bk+1(=5)がd+W(=4)以下ではなくなるので、ステップS87に移行して、作業数Wが3から2にデクリメントされて、作業ブロックは、作業数2の作業済み最右のブロックP(2)、即ちブロック2とされる。そして、ステップS83に移行して、グループ集合G(2)(=G10)に未スキャン部分があるので、ステップS76に移行する。
【0177】
以下同様にして、候補クラス2−3についてブロック3を展開して得られたすべての等価クラスが出力され、候補クラス2−3についてブロック4を展開して得られたすべての等価クラスが出力されると、ステップS82でブロック4の最後の等価クラスについて類似判定を行った後のステップS83では、グループ集合G(4)に未スキャン部分がなくなるので、ステップS84でグループ集合G(4)をメモリ上から解放して、作業数Wを3から2にデクリメントする。ステップS85では、W(=2)は0ではないので、ステップS86で、絞込み対象候補クラスAを作業数2の絞込み対象候補クラスQ(2)即ち、候補クラス2−3とする。そして、ステップS78で、作業数Wを2から3にデクリメントする。すると、bk+1(=5)はd+W(=4)以下ではなくなるので、ステップS87に移行して、作業数Wが3から2にデクリメントされて、作業ブロックbkは、作業数2の作業済み最右ブロックP(2)、即ちブロック2とされる。
【0178】
ステップS83では、グループ集合G(2)(=G10)について未スキャン部分の有無が判断されるが、このときは既に候補クラス2−3までスキャンが終わっているので、ステップS84に移行して、グループ集合G(2)(=G10)がメモリ上から解放されて、作業数Wが2から1にデクリメントされる。そして、ステップS86では、絞込み対象候補クラスAが作業数1の絞込み対象クラスQ(1)である候補クラス1−1とされる。そして、ステップS78で、作業数Wが1から2へとインクリメントされて、ステップS79では、bk+1(=3)はd+W(=3)以下であるので、ステップS80に移行して、作業ブロックbkがブロック2からブロック3に変更される。
【0179】
続くステップS75では、絞込み対象候補クラスAである候補クラス1−1をブロック数bk(=ブロック3)について展開・ソートして、グループ集合G(3)を生成し、作業数2の作業済み最右のブロックP(2)を作業ブロックbk(=ブロック3)とする。以下、上記と同様にして、ステップS76、S77、S78、S79、S80を経てステップS76に戻り、グループ集合G(4)の最初の候補ブロックが、作業数3の絞込み対象候補クラスQ(3)とされると、そのグループ集合G(4)について、ステップS76、S77、S81、S82、S83が繰り返されて、グループ集合G(4)に含まれる等価クラスの出力および出力された等価クラスについての類似判定が行われる。
【0180】
グループ集合G(4)の最後の等価クラスについての類比判定がステップS82で行われると、続くステップS83では、グループ集合G(4)に未スキャン部分がないと判断されて、ステップS84で、グループ集合G(4)がメモリ上から解放されて、作業数Wが3から2にデクリメントされる。続くステップS85では、作業数W(=2)は0ではないので、絞込み対象候補クラスAは作業数2の絞込み対象候補クラスQ(2)、即ちグループ集合G(3)の最初の候補クラスとされて、ステップS78で作業数Wが2から3にインクリメントされる。ステップS79では、bk+1(=5)はd+W(4)以下ではないので、ステップS87に移行して、作業数Wが3から2にデクリメントされて、作業ブロックbkが作業数2の作業済み最右のブロックP(2)、即ちブロック3とされる。
【0181】
そして、ステップS83で、グループ集合G(3)に未スキャン部分があるか判断される。グループ集合G(3)に未スキャン部分がある場合には、ステップS76に戻って、グループ集合G(3)の次の候補クラスを選択して、ブロック4について展開・ソートし、等価クラスを出力するという上述の処理を繰り返す。このようにして、グループ集合G(3)の最後の等価クラスについて展開したグループ集合G(4)の最後の等価クラスについて、ステップS82で類比判定を行うと、続くステップS83では、グループ集合G(4)に未スキャン部分がないことになり、ステップS84にてグループ集合G(4)をメモリ上から解放した後に、作業数Wが3から2にデクリメントされる。
【0182】
ステップS85では、作業数(=2)が0でないので、ステップS86に移行し、絞込み対象候補クラスAが作業数2の絞込み対象候補クラスQ(2)、即ちグループ集合G(3)の最後の候補ブロックとされる。そして、ステップS78で作業数Wが2から3にインクリメントされると、ステップS79では、bk+1(=5)がd+W(=4)以下ではないので、ステップS87に移行して、作業数Wが3から2にデクリメントされて、作業ブロックbkは、作業数2の作業済み最右ブロックP(2)、即ちブロック3とされて、ステップS83でグループ集合G(3)に未スキャン部分があるかが判断される。グループ集合G(3)にはすでに未スキャン部分はないので、ステップS84で、このグループ集合G(3)は、メモリ上から解放され、作業数Wは2から1にデクリメントされる。
【0183】
続いて、ステップS85では、作業数W(=1)は0ではないので、ステップS86に移行して、絞込み対象候補クラスAは、作業数1の絞込み対象候補クラスQ(1)、即ちグループ集合G(1)の最初の候補ブロックとされる。次に、ステップS78で、作業数Wが1から2にインクリメントされて、ステップS79では、bk+1(=4)がd+W(=2)以下ではないので、ステップS87に移行して、作業数Wが2から1にデクリメントされ、作業ブロックは、作業数1の作業済み最右のブロックP(1)、即ちブロック1とされる。そして、ステップS83では、グループ集合G(1)(=G9)については、まだ候補クラス1−1の処理が終わっただけで、まだ未スキャン部分があるので、ステップS76に戻る。
【0184】
グループ集合G9の候補クラス1−1についてした上述の処理を、候補クラス1−2、1−3、1−4についても行い、最後に候補クラス1−4をブロック3について展開してできたグループ集合G(3)の最後の候補クラスをブロック4について展開してできたグループ集合G(4)の最後の等価クラスについて、ステップS82で類比の判定を行うと、その後の処理は以下の通りである。ステップS83では、グループ集合G(4)に未スキャン部分がないので、ステップS84で、グループ集合G(4)がメモリから解放され、作業数Wが3から2にデクリメントされる。ステップS85では、作業数W(=2)は0ではないので、ステップS86で、絞込み対象候補クラスAが作業数2の絞込み対象候補クラスQ(2)、即ちグループ集合G(3)の最後の候補クラスとされる。
【0185】
ステップS78では、作業数Wが2から3にインクリメントされるが、ステップS79では、bk+1(=5)はd+W(=4)以下ではないので、ステップS87にて、作業数Wが3から2にデクリメントされて、作業ブロックbkは、作業数2の作業済み最右のブロックP(2)、即ち、ブロック3とされる。そして、ステップS83では、グループ集合G(3)について未スキャン部分がないので、ステップS84に移行して、グループ集合G(3)がメモリ上から解放され、作業数Wが2から1にデクリメントされる。
【0186】
ステップS85では、作業数W(=1)はまだ0ではないので、ステップS86に移行して、絞込み対象候補クラスAが作業数1の絞込み対象候補クラスQ(1)、即ちグループ集合G(1)の最後の候補クラスとされる。そして、ステップS78で、作業数Wが1から2にインクリメントされるが、ステップS79では、bk+1(=4)がd+W(=2)以下ではないので、ステップS87に移行して、作業数が2から1にデクリメントされて、作業ブロックbkは、作業数1の作業済み最右ブロックP(1)、即ちブロック1とされる。
【0187】
続くステップS83では、グループ集合G(1)(=G9)に既に未スキャン部分がないため、ステップS84に移行して、グループ集合G(1)がメモリ上から解放され、作業数Wは1から0にデクリメントされる。そして、ステップS85に移行すると、作業数Wは0であることから、ステップS88に移行する。ステップS88では、ステップS73で設定した作業開始ブロックSが1であるため、S(=1)≦d(=1)を満たし、ステップS89に進んで、作業開始ブロックSがブロック1からブロック2に変更され、作業数Wは1にリセットされる。この状態で、ステップS74に戻って、絞込み対象候補クラスAがすべてのショートリードとされて、ステップS75に移行し、その後は、上記で説明したのと同様にして、処理が行われる。
【0188】
そして、作業開始ブロックSをブロック2とした処理がすべて終了して、ステップS85でW=0となってステップS88に移行すると、ステップS88では、S(=2)がd(=1)以下ではないので、処理が終了する。
【0189】
以上のように、第2の実施の形態の配列解析装置によれば、生成された候補クラスについて、当該候補クラスを含むグループ集合の他の候補クラスを展開・ソートする処理を行う前に、当該候補クラスを利用するすべての処理を行うので、第1の実施の形態のように、候補クラスを後のブロックセットで再利用するためにハードディスクに記憶しておく必要がない。また、そのようにハードディスクを不要としても、必要なメモリ容量は第1の実施の形態とほぼ変わらない。
【産業上の利用可能性】
【0190】
以上のように、本発明にかかる配列解析装置は、計算コストを削減して、高速に、所定の編集距離以内にあるすべてのショートリードのペアを見つけ出すことができるという効果を有し、大量の塩基配列の中から、距離が所定範囲内にある塩基配列のペアを探索する配列解析装置等に適用できる。
【符号の説明】
【0191】
10 配列解析装置
11 配列入力部
12 条件入力部
13 等価クラス生成部
131 ブロック分割部
132 ソート対象指定部
133 オフセット部
134 ソート部
135 候補クラス生成部
136 候補クラス記憶部
14 類似判定部
141 正準性判定部
142 距離計算部
143 類似ペア出力部
G1〜G14 グループ集合

【特許請求の範囲】
【請求項1】
複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置であって、
探索対象の複数の塩基配列を入力する配列入力部と、
編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力する条件入力部と、
前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成部と、
前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内であると計算されたペアを示すデータを出力する類似判定部と、
を備えたことを特徴とする配列解析装置。
【請求項2】
複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置であって、
探索対象の複数の塩基配列を入力する配列入力部と、
編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)と、前記配列入力部に入力される前記複数の塩基配列の挿入および欠損の総和の閾値である最大ギャップg(gはd以下の自然数)を入力する条件入力部と、
前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記最大ギャップgによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成部と、
前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内かつ、挿入および欠損の総和がg以内であると計算されたペアを示すデータを出力する類似判定部と、
を備えたことを特徴とする配列解析装置。
【請求項3】
前記等価クラス生成部は、前記n個のブロックのうちの1つのブロックについて前記窓枠を前記最大オフセット内でオフセットを変えて読み出した部分配列または当該ブロック内の部分配列が一致するという絞り込み条件を満たす塩基配列の集合を、前記等価クラスの候補となる候補クラスとして求め、当該候補クラスに含まれる塩基配列に対して、前記n個のブロックのうちの別のブロックについて前記窓枠を前記最大オフセット内でオフセットを変えて読み出した部分配列または当該別のブロック内の部分配列が前記絞り込み条件を満たす候補クラスを求めるという処理を、前記n個のブロックについて順次行い、n個のすべてのブロックを用いて絞り込まれた候補クラスを等価クラスとして生成することを特徴とする請求項1または2に記載の配列解析装置。
【請求項4】
前記等価クラス生成部は、
塩基配列に対して、絞り込みに用いられるブロックから前記最大オフセット内のあらゆるオフセット値で前記窓枠をオフセットさせて部分配列を読み出すオフセット部と、
前記窓枠をオフセットさせて読み出された部分配列および前記窓枠をオフセットしないで読み出された部分配列を塩基の並び方に応じてソートして同じ配列パターンの部分配列が続く範囲をグループ化するソート部と、
前記ソート部にて生成されたグループに複数の部分配列が含まれ、かつ、当該グループの中の少なくとも1つの部分配列のオフセット値が0であるときに、そのグループを候補クラスとして求める候補クラス生成部と、
を備えることを特徴とする請求項3に記載の配列解析装置。
【請求項5】
前記候補クラス生成部は、
候補クラスに含まれる塩基配列に関連付けて、候補クラスを絞り込んだブロックの位置、および、候補クラスに含まれる他の部分配列と一致するために用いたオフセット値をメモリに記憶し、
前記n個のブロックのうちのm(m<n)個目のブロックを用いて候補クラスを求める際には、前記1つのグループの中で、オフセット値が0であると判定された塩基配列の中に、前のm−1個のブロックを用いた候補クラスの絞り込みにおいて、当該塩基配列の部分配列のオフセット値が連続して0であった塩基配列が存在しない場合には、当該グループを候補クラスとしない請求項4に記載の配列解析装置。
【請求項6】
前記候補クラス生成部は、
1つのグループ内に、オフセット値の異なる同じ塩基配列の部分配列が存在する場合において、それらのオフセット値の中に0が含まれるときは0を、含まれないときはそれらのオフセット値の中から任意に選択した値を当該部分配列のオフセット値として記憶する請求項5に記載の配列解析装置。
【請求項7】
前記類似判定部は、
前記メモリから、前記等価クラスに含まれる塩基配列のブロックの位置、および、オフセット値を読み出し、
等価クラスに含まれる塩基配列の任意のペアのうち、少なくとも一方の塩基配列のオフセット値がすべてのブロックで0であり、かつ、所定のブロック順位規則に基づく順位において前記ブロックの位置よりも前に、前記最大オフセット内でオフセットさせることにより一致する部分配列がないときに、当該ペアの塩基配列について、その編集距離を計算することを特徴とする請求項5または6に記載の配列解析装置。
【請求項8】
前記等価クラス生成部は、
前記候補クラス生成部によって複数の候補クラスが生成された場合には、そのうちの1つの候補クラスについて次のブロックを用いた候補クラスの絞り込みを行なう処理を、等価クラスが求まるか、または候補クラスがなくなるまで繰り返し行い、等価クラスが求まるか、または候補クラスがなくなった時点で1つ前のブロックを用いた処理で求まった他の候補クラスについて絞り込みを行うことを特徴とする請求項4ないし請求項7のいずれかに記載の配列解析装置。
【請求項9】
前記等価クラス生成部は、最初のブロックを用いて求められたすべての候補クラスについての絞り込みを終了したときに、組み合わせの異なるn個のブロックを選択し、新たに選択されたn個のブロックを用いて等価クラスを生成することを特徴とする請求項8に記載の配列解析装置。
【請求項10】
前記等価クラス生成部は、前記新たに選択されたn個のブロックを用いて等価クラスを生成する場合において、前記候補クラス生成部において利用できる候補クラスが生成されているときは、当該既に生成されている候補クラスを利用することで、前記ソート部において同一の候補クラスが複数回整列されるのを回避することを特徴とする請求項9に記載の配列解析装置。
【請求項11】
前記等価クラス生成部は、前記候補クラスが生成されるごとに、当該候補クラスを利用するn個のブロックの組合せのすべてについて、等価クラスが求まるまで絞込みを行い、当該候補クラスを利用するすべてのn個のブロックの組合せについて等価クラスを求めた後に、当該候補クラスをメモリ上から解放することを特徴とする請求項8に記載の配列解析装置。
【請求項12】
前記配列入力部に入力される前記複数の塩基配列は所定の範囲D内で長さが異なっており、
前記条件入力部に入力される編集距離dは、d≧Dを満たし、
前記最大オフセットは、d/2以上の最小の整数である
ことを特徴とする請求項1ないし請求項8のいずれかに記載の配列解析装置。
【請求項13】
前記配列入力部に入力される前記複数の塩基配列は所定の範囲D内で長さが異なっており、
前記条件入力部に入力される編集距離dは、d≧Dを満たし、
前記最大オフセット量は、g/2以上の最小の整数である
ことを特徴とする請求項2に記載の配列解析装置。
【請求項14】
前記配列入力部に入力される前記複数の塩基配列は長さが互いに等しく、
前記最大オフセットは、d/2以下の最大の整数である
ことを特徴とする請求項1ないし請求項8のいずれかに記載の配列解析装置。
【請求項15】
前記配列入力部に入力される前記複数の塩基配列は長さが互いに等しく、
前記最大オフセットは、g/2以下の最大の整数である
ことを特徴とする請求項2に記載の配列解析装置。
【請求項16】
前記数値nまたは前記分割数(d+n)は、前記配列解析装置が所定のアルゴリズムに従って決定して前記条件入力部に入力されることを特徴とする請求項1ないし請求項15のいずれかに記載の配列解析装置。
【請求項17】
複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析方法であって、
探索対象の複数の塩基配列を入力する配列入力ステップと、
編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力する条件入力ステップと、
前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成ステップと、
前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内であると計算されたペアを示すデータを出力する類似判定ステップと、
を含むことを特徴とする配列解析方法。
【請求項18】
複数の塩基配列の中から、所定の範囲内にある塩基配列のペアを探索する配列解析方法であって、
(1)探索対象の複数の塩基配列を入力するステップと、
(2)編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力するステップと、
(3)前記複数の塩基配列の各々を前記分割数(d+n)個のブロックに分割するとともに、各ブロックにブロック番号を付与するステップと、
(4)前記(d+n)個のブロックの中からn個のブロックをブロックセットとして選択するステップと、
(5)作業数を1とし、前記選択されたn個のブロック中のブロック番号が最小のブロックを作業ブロックとし、かつ、ステップ(1)で入力されたすべての塩基配列を候補クラスとするステップと、
(6)候補クラス内のすべての塩基配列について、前記作業ブロック内の部分配列と、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該作業ブロックからオフセットさせて読み出した部分配列からなる部分配列の集合を生成するステップと、
(7)前記部分配列の集合を塩基の並び方に応じてソートしてグループ集合を生成するステップと、
(8)前記グループ集合を順にスキャンしていき、同じ配列パターンの部分配列が続く範囲を1つのグループとして、そのグループに複数の部分配列が含まれ、かつ、少なくとも1つの部分配列のオフセット値が0であるときに、作業数がnであるかを判断するステップと、
(9)ステップ(8)において、作業数がnでないときに、そのグループを候補クラスとして報告するステップと、
(10)ステップ(9)の後に、作業数をインクリメントし、前記作業ブロックを、ステップ(4)にて選択されたn個のブロックのうちの次のブロックに変更して、ステップ(6)に戻るステップと、
(11)ステップ(8)において、作業数がnであるときに、そのグループを等価クラスとして報告するステップと、
(12)ステップ(11)の後に、当該等価クラスに含まれる塩基配列の任意のペアのうちの少なくとも一方において、オフセット値の履歴にゼロ以外の値がなく、かつ所定の正準性判断に基づいて正準であると判断されたペアにつき、編集距離を計算するステップと、
(13)ステップ(12)の後に、編集距離がd以内である場合に、当該ペアを示すデータを出力するステップと、
(14)ステップ(13)の後に、前記グループ集合に未スキャン部分があるかを判断するステップと、
(15)ステップ(14)において、未スキャン部分がある場合に、ステップ(8)に戻るステップと、
(16)ステップ(14)において、未スキャン部分がない場合に、作業数をデクリメントして、作業ブロックを1つ前のブロックに戻すステップと、
(17)ステップ(16)の後に作業数が1以上である場合に、ステップ(14)に戻るステップと、
(18)ステップ(16)の後に作業数が0になった場合に、前記グループ集合を消去してまだ選択されていないブロックセットがあるかを判断するステップと、
(19)ステップ(18)において、まだ選択されていないブロックセットがある場合に、ステップ(4)に戻るステップと、
(20)ステップ(18)において、まだ選択されていないブロックセットがない場合に、処理を終了するステップと、
を有することを特徴とする配列解析方法。
【請求項19】
さらに、
(21)ステップ(4)の後に、前記ブロックセットに含まれるブロックをブロック番号が小さいものから並べて、ブロック番号が小さい方から1または複数のブロックにより構成されるサブセットの候補クラスが記憶されているか否かを判断するステップと、
(22)ステップ(21)で、候補クラスが記憶されていると判断された場合に、ステップ(5)に代えて、当該サブセットの候補クラスを読み出して、その候補クラスに対応する作業数および作業ブロックを設定するステップと、
を有することを特徴とする請求項17に記載の配列解析方法。
【請求項20】
複数の塩基配列の中から、所定の範囲内にある塩基配列のペアを探索する配列解析方法であって、
(51)探索対象の複数の塩基配列を入力するステップと、
(52)編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力するステップと、
(53)前記複数の塩基配列の各々を前記分割数(d+n)個のブロックに分割するとともに、各ブロックにブロック番号を付与するステップと、
(54)作業数Wを1とし、作業ブロックbkをブロック1とし、作業開始ブロックSを1とするステップと、
(55)絞込み対象候補クラスAを入力されたすべてのショートリードとするステップと、
(56)絞込み対象候補クラスAの作業ブロックbkを展開・ソートして、グループ集合G(bk)を生成して記憶し、作業数Wの作業済み最右のブロックP(W)を作業ブロックbkとするステップと、
(57)グループ集合G(bk)を上からスキャンして、候補クラスの条件を満たす最初のグループを選択してそのグループを絞り込み対象候補クラスAとし、かつ、作業数Wの絞込み対象候補クラスQ(W)をその絞り込み対象候補クラスAとするステップと、
(58)ステップ(57)の後に、作業数Wがnであるか否かを判断するステップと、
(59)ステップ(58)にて作業数Wがnでない場合に、作業数Wをインクリメントするステップと、
(60)ステップ(59)の後に、bk+1がd+W以下であるか否かを判断するステップと、
(61)ステップ(60)にて、bk+1がd+W以下である場合に、作業ブロックbkをインクリメントして、ステップ(56)に戻るステップと、
(62)ステップ(58)にて作業数Wがnである場合に、絞込み対象候補クラスAを等価クラスとするステップと、
(63)ステップ(62)の後に、等価クラスについて、類似判定をするステップと、
(64)ステップ(63)の後に、グループ集合G(bk)に未スキャン部分があるか否かを判断し、未スキャン部分がある場合に、ステップ(57)に戻るステップと、
(65)ステップ(64)にて、未スキャン部分がない場合に、グループ集合G(bk)をメモリ上から解放し、作業数Wをデクリメントするステップと、
(66)ステップ(65)の後に、W=0であるか否かを判断するステップと、
(67)ステップ(66)にて、W=0でない場合に、絞込み対象候補クラスAを作業数Wの絞込み対象候補クラスとするステップと、
(68)ステップ(60)にて、bk+1がd+Wより大きい場合に、作業数Wをデクリメントし、作業ブロックbkをそのデクリメントされた作業数の作業済み最右のブロックP(W)とするステップと、
(69)ステップ(66)にて、W=0である場合に、作業開始ブロックSが編集距離d以下であるか否かを判断し、作業開始ブロックSが編集距離dより大きい場合に、処理を終了するステップと、
(70)ステップ(69)にて、作業開始ブロックSが編集距離d以下である場合に、作業開始ブロックSをインクリメントして、作業数Wを1にリセットして、ステップ(55)に戻るステップと、
を有することを特徴とする配列解析方法。
【請求項21】
コンピュータを、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置として機能させるコンピュータプログラムであって、該コンピュータプログラムは、コンピュータに、
探索対象の複数の塩基配列を入力する配列入力ステップと、
編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力する条件入力ステップと、
前記複数の塩基配列の各々を前記分割数(d+n)個に分割して(d+n)個のブロックを生成し、その中から選んだn個のすべてのブロックについて、当該ブロックから読み出した部分配列、または、前記編集距離dによって定まる最大オフセットの範囲内で、読出し範囲を特定する窓枠を当該ブロックからオフセットさせて読み出した部分配列が一致するという条件を満たす塩基配列の集合を、等価クラスとして報告する等価クラス生成ステップと、
前記等価クラス生成部から報告された等価クラス内の塩基配列のペアについて、編集距離を計算して、計算によって編集距離d以内であると計算されたペアを示すデータを出力する類似判定ステップと、
を実行させることを特徴とするコンピュータプログラム。
【請求項22】
コンピュータを、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置として機能させるコンピュータプログラムであって、該コンピュータプログラムは、コンピュータに、
(1)探索対象の複数の塩基配列を入力するステップと、
(2)編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力するステップと、
(3)前記複数の塩基配列の各々を前記分割数(d+n)個のブロックに分割するステップと、
(4)前記(d+n)個のブロックの中からn個のブロックをブロックセットとして選択するステップと、
(5)作業数を1とし、作業ブロックを前記選択されたn個のブロック中のブロック番号が最小のブロックとし、かつ、候補クラスをすべてのショートリードとするステップと、
(6)候補クラス内のすべての塩基配列について、前記作業ブロックに対して部分配列を読み出す窓枠を前記編集距離dによって定まる最大オフセット内のあらゆるオフセット値でオフセットさせて読み出した部分配列および前記作業ブロック内の部分配列からなる部分配列の集合を生成するステップと、
(7)前記部分配列の集合を塩基の並び方に応じてソートしてグループ集合を生成するステップと、
(8)前記グループ集合を順にスキャンしていき、同じ配列パターンの部分配列が続く範囲を1つのグループとして、そのグループに複数の部分配列が含まれ、かつ、少なくとも1つの部分配列のオフセット値が0であるときに、作業数がnであるかを判断するステップと、
(9)ステップ(8)において、作業数がnであるときに、そのグループを候補クラスとして報告するステップと、
(10)ステップ(9)の後に、作業数をインクリメントし、前記作業ブロックを、ステップ(4)にて選択されたn個のブロックの次のブロックに変更して、ステップ(6)に戻るステップと、
(11)ステップ(8)において、作業数がnであるときに、そのグループを等価クラスとして報告するステップと、
(12)ステップ(11)の後に、当該等価クラスに含まれる塩基配列の任意のペアのうちの少なくとも一方において、オフセット値の履歴にゼロ以外の値がなく、かつ所定の正準性判断に基づいて正準であると判断されたペアにつき、編集距離を計算するステップと、
(13)ステップ(12)の後に、編集距離がd以内である場合に、当該ペアを示すデータを出力するステップと、
(14)ステップ(13)の後に、前記グループ集合に未スキャン部分があるかを判断するステップと、
(15)ステップ(14)において、未スキャン部分がある場合に、ステップ(8)に戻るステップと、
(16)ステップ(14)において、未スキャン部分がない場合に、作業数をデクリメンと、して、作業ブロックを1つ前のブロックに戻すステップと、
(17)ステップ(16)の後に作業数が1以上である場合に、ステップ(14)に戻るステップと、
(18)ステップ(16)の後に作業数が0になった場合に、前記グループ集合を消去してまだ選択されていないブロックセットがあるかを判断するステップと、
(19)ステップ(18)において、まだ選択されていないブロックセットがある場合に、ステップ(4)に戻るステップと、
(20)ステップ(18)において、まだ選択されていないブロックセットがない場合に、処理を終了するステップと、
を実行させることを特徴とするコンピュータプログラム。
【請求項23】
コンピュータを、複数の塩基配列の中から、所定範囲内にある塩基配列のペアを探索する配列解析装置として機能させるコンピュータプログラムであって、該コンピュータプログラムは、コンピュータに、
(51)探索対象の複数の塩基配列を入力するステップと、
(52)編集距離dと、複数の塩基配列を対応する部分配列ごとに比較する処理で用いる処理単位をブロックとし、塩基配列を何ブロックに分割するかを示す分割数を決定する数値nまたは前記分割数(d+n)を入力するステップと、
(53)前記複数の塩基配列の各々を前記分割数(d+n)個のブロックに分割するとともに、各ブロックにブロック番号を付与するステップと、
(54)作業数Wを1とし、作業ブロックbkをブロック1とし、作業開始ブロックSを1とするステップと、
(55)絞込み対象候補クラスAを入力されたすべてのショートリードとするステップと、
(56)絞込み対象候補クラスAの作業ブロックbkを展開・ソートして、グループ集合G(bk)を生成して記憶し、作業数Wの作業済み最右のブロックP(W)を作業ブロックbkとするステップと、
(57)グループ集合G(bk)を上からスキャンして、候補クラスの条件を満たす最初のグループを選択してそのグループを絞り込み対象候補クラスAとし、かつ、作業数Wの絞込み対象候補クラスQ(W)をその絞り込み対象候補クラスAとするステップと、
(58)ステップ(57)の後に、作業数Wがnであるか否かを判断するステップと、
(59)ステップ(58)にて作業数Wがnでない場合に、作業数Wをインクリメントするステップと、
(60)ステップ(59)の後に、bk+1がd+W以下であるか否かを判断するステップと、
(61)ステップ(60)にて、bk+1がd+W以下である場合に、作業ブロックbkをインクリメントして、ステップ(56)に戻るステップと、
(62)ステップ(58)にて作業数Wがnである場合に、絞込み対象候補クラスAを等価クラスとするステップと、
(63)ステップ(62)の後に、等価クラスについて、類似判定をするステップと、
(64)ステップ(63)の後に、グループ集合G(bk)に未スキャン部分があるか否かを判断し、未スキャン部分がある場合に、ステップ(57)に戻るステップと、
(65)ステップ(64)にて、未スキャン部分がない場合に、グループ集合G(bk)をメモリ上から解放し、作業数Wをデクリメントするステップと、
(66)ステップ(65)の後に、W=0であるか否かを判断するステップと、
(67)ステップ(66)にて、W=0でない場合に、絞込み対象候補クラスAを作業数Wの絞込み対象候補クラスとするステップと、
(68)ステップ(60)にて、bk+1がd+Wより大きい場合に、作業数Wをデクリメントし、作業ブロックbkをそのデクリメントされた作業数の作業済み最右のブロックP(W)とするステップと、
(69)ステップ(66)にて、W=0である場合に、作業開始ブロックSが編集距離d以下であるか否かを判断し、作業開始ブロックSが編集距離dより大きい場合に、処理を終了するステップと、
(70)ステップ(69)にて、作業開始ブロックSが編集距離d以下である場合に、作業開始ブロックSをインクリメントして、作業数Wを1にリセットして、ステップ(55)に戻るステップと、
を実行させることを特徴とするコンピュータプログラム。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate

【図11】
image rotate

【図12】
image rotate

【図13】
image rotate

【図14】
image rotate

【図15】
image rotate

【図16】
image rotate

【図17】
image rotate

【図18】
image rotate

【図19】
image rotate

【図20】
image rotate

【図21】
image rotate


【公開番号】特開2012−18592(P2012−18592A)
【公開日】平成24年1月26日(2012.1.26)
【国際特許分類】
【出願番号】特願2010−156342(P2010−156342)
【出願日】平成22年7月9日(2010.7.9)
【出願人】(301021533)独立行政法人産業技術総合研究所 (6,529)
【Fターム(参考)】