JP2014089690A - 塩基配列アラインメントシステム及び方法 - Google Patents

Abstract

【課題】マッピングの正確性を保障しマッピング時の複雑度を改善し処理速度を上げる。
【解決手段】リード配列から複数の断片（fragment）配列を生成する断片配列生成部と、生成された該複数の断片配列から候補断片配列集合を構成するフィルタリング部と、該候補断片配列それぞれの該参照配列でのマッピング位置の個数を計算し、計算された該マッピング位置の個数が設定された値を超える断片配列を選択し、該参照配列でのマッピング位置の個数が該設定された値以下になるまで選択された断片配列のサイズを拡張する断片配列拡張部と、該参照配列を複数の区間に分割し該各区間別に該候補断片配列の総マッピング長さを計算するマッピング長さ計算部と、計算された該総マッピング長さが基準値以上の区間を選択し、選択された区間に対して該リード配列に対する全域アラインメント（global alignment）を遂行するアラインメント部とを備える。
【選択図】図１

Description

本発明は、ゲノムの塩基配列を分析するための技術に関するものである。

低いコスト及び速いデータ生産によって大容量の短い配列を生産する次世代シーケンシング（ＮＧＳ：Next Generation Sequencing）が伝統的なサンガー（Sanger）シーケンシング方式を迅速に置き換えている。また、多様なＮＧＳ配列の組換えプログラムが正確度に焦点を合わせて開発されてきた。しかしながら、最近、次世代シーケンシング技術が発展することによって、断片配列を作成するコストが以前の半分以下となり、これにより用いることができるデータの量が多くなったため、大容量の短い配列を急速かつ正確に処理するための技術が必要となった。

配列組換えの第一の段階は、塩基配列アラインメント（alignment）アルゴリズムを通じてリードを参照配列の正確な位置にマッピング（mapping）するものである。これにおける問題点は、同一種の個体であるとしても多様な遺伝的変異によって、ゲノムの配列に差があり得るという点である。また、シーケンシング過程におけるエラーによっても塩基配列に差が生じ得る。従って、塩基配列アラインメントアルゴリズムは、このような差及び変異を効果的に考慮してマッピングの正確度を高めるべきである。

それ故に、ゲノム情報の分析を行うためには、可能な限り多くの正確な全体ゲノム情報データが必要である。なお、このためには、何よりも優れた正確性及び大量の処理量を有する塩基配列アラインメントアルゴリズムを開発することを先行して行わなければならない。しかしながら、従来の方法は、このような要求条件を満たすのに限界があった。

本発明は、マッピングの正確性を保障すると同時にマッピング時の複雑度を改善し、処理速度を上げることができる塩基配列アラインメント手段を提供することをその目的としている。

前記課題を解決するための、本発明の一実施形態による塩基配列アラインメントシステムは、リード（read）配列を参照配列にアラインメントするためのシステムであって、前記リード配列から複数の断片（fragment）配列を生成する断片配列生成部と、生成された前記複数の断片配列から候補断片配列集合を構成するフィルタリング部と、前記候補断片配列それぞれの前記参照配列でのマッピング位置の個数を計算し、計算された前記マッピング位置の個数が、設定された値を超える断片配列を選択し、前記参照配列でのマッピング位置の個数が前記設定された値以下になるまで選択された断片配列のサイズを拡張する断片配列拡張部と、前記参照配列を複数の区間に分割し、前記各区間別に前記候補断片配列の総マッピング長さを計算するマッピング長さ計算部と、計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行するアラインメント部とを備える。

一方、前記課題を解決するための本発明の一実施形態によるリード（read）配列を参照配列にアラインメントするための方法は、断片配列生成部において、前記リード配列から複数の断片（fragment）配列を生成する段階と、フィルタリング部において、生成された前記複数の断片配列から候補断片配列集合を構成する段階と、断片配列拡張部において、生成された前記候補断片配列それぞれの前記参照配列でのマッピング位置の個数を計算する段階と、前記断片配列拡張部において、計算された前記マッピング位置の個数が設定された値を超える断片配列を選択する段階と、前記断片配列拡張部において、前記参照配列でのマッピング位置の個数が前記設定された値以下になるまで選択された断片配列のサイズを拡張する段階と、マッピング長さ計算部において、前記参照配列を複数の区間に分割し、前記各区間別に前記候補断片配列の総マッピング長さを計算する段階と、アラインメント部において、計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行する段階とを備える。

一方、前記課題を解決するための本発明の他の実施形態による塩基配列アラインメントシステムは、リード（read）配列を参照配列にアラインメントするためのシステムであって、前記リード配列から複数の断片（fragment）配列を生成する断片配列生成部と、生成された前記複数の断片配列から候補断片配列集合を構成するフィルタリング部と、前記参照配列を複数の区間に分割し、前記各区間別に前記候補断片配列の総マッピング長さを計算するマッピング長さ計算部と、計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行するアラインメント部とを備える。

本発明の実施形態による場合、リード配列のアラインメント時にリード配列の特定領域のみを考慮するものではなく、リード全体を考慮してシード(断片配列)を選択するため、リードの一部分のみを考慮するアルゴリズムと比較して正確性を向上させることができる。

また、断片配列別に参照塩基配列内における繰り返し数を制限し、これを超えるシードに対しては、シードの長さを拡張することによって、マッピングの正確性を高めると同時に速度もまた向上させることができる効果がある。

また、マッピングヒストグラムを利用して、参照塩基配列内にリードがマッピングされる可能性が高い特定領域を選択し、該当領域内でのみ全域アラインメント（global alignment）を遂行することによって、全域アラインメントの時間を大幅に減少させることができる。

さらに、リードから導き出された断片配列のマッピング位置と組み合わせを検索する複雑な過程の代わりに、組み合わせを構成する可能性が高い断片配列に対して直ちに全域アラインメントを遂行することによって、全域アラインメントの速度をさらに高めることができ、全域アラインメントの位置を記憶しておき、該当位置の周辺では繰り返して全域アラインメントを遂行しないように構成することによって、不要な全域アラインメントの回数を減少させることができるといった長所がある。

本発明の一実施形態による塩基配列アラインメント方法１００を説明するための図である。 本発明の一実施形態による塩基配列アラインメント方法１００における段階１０８でのｍＥＢの計算過程を例示するための図である。 本発明の一実施形態による塩基配列アラインメント方法１００における段階１１２での断片配列の生成過程を説明するための図である。 本発明の一実施形態によるマッピングヒストグラムを利用したマッピング対象区間の選択過程を例示した図である。 本発明の一実施形態による全域アラインメント時に不要な全域アラインメント回数を減少させるための方法を説明するための例示図である。 本発明の一実施形態による塩基配列アラインメントシステム６００を示したブロック図である。

以下、図面を参照して、本発明の具体的な実施形態を説明する。しかしながら、これらは例示に過ぎず、本発明はこれらに限定されるものではない。

本発明を説明するにあたり、本発明に係る公知技術に対する具体的な説明が、本発明の要旨を不明瞭にすると判断される場合には、その詳細な説明を省略する。そして、後述する用語は、本発明における機能を考慮して定義された用語であって、これは使用者、運用者の意図又は慣例などにより変わり得る。従って、その定義は本明細書の全般にわたる内容に基づいて下されなければならない。

本発明の技術的思想は特許請求の範囲によって定められ、以下の実施形態は、本発明の技術的思想を本発明が属する技術分野において通常の知識を有する者に効率的に説明するための一手段に過ぎない。

本発明の実施形態を詳しく説明するに先立ち、まず、本発明において用いられる用語に関して説明すると以下のとおりである。

先ず、「リード（read）配列」(又は、略して「リード」と称する)とは、ゲノムシーケンサ（genome sequencer）から出力される短い長さの塩基配列データである。リード配列の長さは、ゲノムシーケンサの種類によって、一般的に３５〜５００ｂｐ（base pair）程度で多様に構成され、一般的にＤＮＡ塩基の場合、Ａ、Ｃ、Ｇ、Ｔのアルファベット文字で表現される。

「参照塩基配列」とは、前記リード配列から全体塩基配列を生成するのに参照となる塩基配列（reference sequence）を意味する。塩基配列の分析では、ゲノムシーケンサから出力される多量のリードを、参照塩基配列を参照してマッピングすることによって、全体塩基配列を完成するようになる。本発明において、前記参照塩基配列は、塩基配列の分析時に既に設定された配列(例えば、ヒトの全体塩基配列等)であり得、或いはゲノムシーケンサから出力された塩基配列を参照塩基配列として用いることもできる。

「塩基（base）」は、参照塩基配列及びリードを構成する最小単位である。前述したように、ＤＮＡ塩基の場合、Ａ、Ｃ、Ｇ、及びＴの４種類のアルファベット文字で構成され得、これらそれぞれを塩基と表現する。換言すれば、ＤＮＡ塩基の場合、４種の塩基で表現され、これはリード配列もまた同様である。

「断片（fragment）配列」（又は、シード（seed））とは、リード配列のマッピングのためにリード配列と参照塩基配列とを比較する時の単位となるシーケンスである。理論的に、リードを参照塩基配列にマッピングするためには、リードの全体を参照塩基配列の最初の部分から順次比較していきながら、リードのマッピング位置を計算すべきである。しかしながら、このような方法の場合、一つのリードをマッピングするのに非常に長い時間及び計算能力（コンピュータパワー）が求められるため、実際ではリードの一部分で構成された小片である断片配列を先に参照塩基配列にマッピングすることによって、全体リード配列のマッピング候補位置を検出し、該当候補位置に全体リード配列をマッピング（global alignment）するようになる。

図１は、本発明の一実施形態による塩基配列アラインメント方法１００を説明するための図である。本発明の実施形態において、塩基配列アラインメント方法１００は、ゲノムシーケンサ（genome sequencer）から出力されるリード配列を参照塩基配列と比較してリード配列の前記参照配列でのマッピング(又は、アラインメント)位置を決定する一連の過程を意味する。

先ず、ゲノムシーケンサ（genome sequencer）からリード配列が入力されると（段階１０２）、リード配列全体と前記参照塩基配列との一致・整合（exact matching）を試みる（段階１０４）。仮に、前記試みの結果、リード全体に対する一致・整合が成功した場合には、以後のアラインメント段階を遂行せずにアラインメントに成功したものと判断する（段階１０６）。

ヒトの塩基配列を対象とした実験の結果、ゲノムシーケンサから出力される１００万個のリード配列をヒトの塩基配列に一致・整合する場合、総２００万回のアラインメントのうち(順方向シーケンス１００万回、逆相補（reverse complement）方向シーケンス１００万回)２３１,５６４回の一致・整合が生じるものと示された。従って、前記段階１０４の遂行結果、約１１．６％程度のアラインメント所要量を減少させることができた。

しかしながら、これとは異なり、前記段階１０６で該当リード配列が、一致・整合しないものと判断される場合には、該当リード配列を前記参照配列にアラインメントした時に生じ得るエラーの最小個数(ｍＥＢ：minimum error bound)を計算する（段階１０８)。

図２は、前記段階１０８でのｍＥＢ計算過程を例示するための図である。先ず、図２の（ａ）に示すように、最初ｍＥＢを０に設定して、リード配列の一番目の塩基から右側に一塩基ずつ移動しながら一致・整合を試みる。この時、（ｂ）に示すように、リード配列の特定塩基(図において二番目のＴで示した部分)からこれ以上一致・整合が不可能であると仮定する。この場合は、リード配列の整合の始めの位置から現在の位置の間の区間のどこかでエラーが生じたことを意味する。従って、この場合には、ｍＥＢ値を１ほど増加させ(ｍＥＢ＝１)、次の位置で新たに一致・整合を開始する(図において、（ｃ）に示す)。以後、再び一致・整合が不可能であると判断される場合には、一致・整合を新たに開始した位置から現在位置の間の区間のどこかでまたエラーが生じたものであるため、ｍＥＢ値を再び１ほど増加させ(ｍＥＢ＝２)、次の位置で新たに一致・整合を開始する(図において、（ｄ）に示す)。このような過程を通じて、リードの最後まで到逹した場合のｍＥＢ値が該当リードのｍＥＢ値となる。

前記のような過程を通じて、リード配列のｍＥＢ値を計算すると、計算されたｍＥＢ値が既設定された最大エラー許容値(max error)を超えるか否かを判断し（段階１１０)、超える場合には、該当リード配列に対するアラインメントが失敗したものと判断してアラインメントを終了する。

前述したヒトの塩基配列を対象とした実験で、最大エラー許容値(max error)を３とし、残りのリードのｍＥＢを計算した結果、総８４４,８９１回に該当するリードが前記最大エラー許容値を超えるものと示された。すなわち、前記段階１０８の遂行結果、約４２．２％程度のアラインメント所要量を減少させることができた。

しかしながら、これとは異なり、前記段階１１０における判断の結果、計算されたｍＥＢ値が前記最大エラー許容値以下である場合には、次のような過程を通じて、該当リード配列に対するアラインメントを遂行する。

先ず、前記リード配列から複数の断片（fragment）配列を生成し（段階１１２)、生成された前記複数の断片配列のうち、前記参照配列とマッチングする断片配列のみを含む候補断片配列集合を構成する(段階１１４)。以後、前記参照配列を複数の区間に分割し、前記各区間別に前記候補断片配列の総マッピング長さを計算することによって、マッピングヒストグラムを生成し（段階１１６)、生成されたマッピングヒストグラムから前記総マッピング長さが基準値（histogram cut）以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行するようになる（段階１１８）。この時、前記全域アラインメントの結果、リードのエラー個数が既設定された最大エラー許容値（max error）を超える場合にはアラインメント失敗と判断され、そうではない場合にはアラインメントに成功したものと判断される（段階１２０）。

以下では、前記段階１１２乃至段階１１８の具体的な過程を詳しく説明する。

リード配列から複数の断片配列生成（段階１１２）
本段階は、本格的にリード配列のアラインメントを遂行するためにリード配列から複数の小片である断片配列を生成する段階である。本段階では、前記リード配列の最初の塩基から最後の塩基まで設定された間隔（shift size）ほど移動しながら設定されたサイズ（fragment size）ほどリード配列の値を読み取ることによって、前記断片配列を生成するようになる。

図３は、前記段階１１２での断片配列の生成過程を説明するための図である。図３は、リード配列の長さが７５ｂｐ（base pair）、リードの最大エラー許容許容値が３ｂｐ、断片配列のサイズ（fragment size）が１５ｂｐ、移動間隔（shift size）が４ｂｐである場合の実施形態を示したものである。すなわち、リード配列の最初の塩基から４ｂｐずつ右側に移動しながら断片配列を生成するようになる。ただし、示した実施形態は、単なる例示に過ぎず、例えば、前記移動間隔、断片配列のサイズ等は、リード配列の長さ、リードの最大エラー許容値等の値を考慮して適切に定められ得る。換言すれば、本発明の権利範囲は、特定の断片配列のサイズ及び移動間隔に限定されるものではないことに留意されたい。

生成された断片配列のフィルタリング及び拡張（段階１１４）
前記のような過程を通じて断片配列が生成されると、次に生成された断片配列のうち参照配列とマッチングしない断片配列を除くフィルタリング過程を通じて、候補断片配列集合（sub-candidate）を構成する。すなわち、生成された断片配列と前記参照配列との一致・整合（exact matching）を試み、その結果、不一致となる塩基の数が既設定された許容値以下である断片配列(候補断片配列)で前記候補断片配列集合を構成するようになる。この時、前記許容値が０である場合、前記候補断片配列集合には、前記参照配列と一致・整合する断片配列のみが含まれる。

例えば、図３に示す実施形態では、前記リードの１５番目、３１番目、及び４７番目の部分でエラーが生じたものと仮定する(図において点線で示す)。この場合、前記エラーを含む断片配列(図において灰色で表示)の場合には、参照配列との一致・整合ができず、エラーの影響を受けない４９〜６３、５３〜６７、５７〜７１、及び６１〜７５の４個の断片配列のみが参照配列と一致・整合する。従って、この場合、前記候補断片配列集合には、上述した４個の断片配列のみが含まれる。

その一方で、一般的に参照塩基配列(例えば、ヒトのゲノム)は、多数の反復シーケンス（repeat sequence）を含む。このような反復シーケンスは、参照配列の様々の位置に分布し、同一の塩基配列を繰り返しで含むため、一部の断片配列の場合、参照配列とのマッピングの時、極めて多い位置で一致・整合が生じるようになる。しかしながら、このような反復シーケンスによって、一部の断片配列で極めて多数のマッピングが生成される場合、全体アラインメントアルゴリズムの複雑度及び正確度に悪影響を与えるため、この場合、適切な方法を利用してマッピングされる位置の反復数を減らす必要がある。

そのために、本段階では、候補断片配列のうち、前記参照配列でのマッピング位置が既設定された値(例えば、５０個)を超える場合、前記マッピング位置の個数が前記設定された値以下になるまで該当断片配列のサイズを拡張する段階をさらに含むことができる。

具体的に、本段階では、生成された前記候補断片配列それぞれの前記参照配列でのマッピング位置の個数を計算し、計算された前記マッピング位置の個数が設定された値を超える断片配列を選択した後、前記参照配列でのマッピング位置の個数が前記設定された値以下になるまで選択された断片配列のサイズを拡張するようになる。この時、前記選択された断片配列のサイズ拡張は、前記選択された断片配列の最初又は最後の部分に、該当位置に対応する前記リード配列での塩基を追加することによって、遂行することができる。

これを例を挙げて説明すると、次のとおりである。例えば、次のようにリード配列から断片配列が生成されると仮定する。

リード配列：Ａ Ｔ Ｔ Ｇ Ｃ Ｃ Ｔ Ｃ Ａ Ｇ Ｔ
断片配列：Ｔ Ｔ Ｇ Ｃ（リード配列で下線を引いた部分）

仮に、前記断片配列に対するマッピングの結果、参照配列でのマッピング位置の個数が基準値である５０個を超える６５個である場合には、下記のように前記マッピング位置の個数が基準値以下に低下するまで前記断片配列の長さを１ｂｐずつ拡張するようになる。

Ｔ Ｔ Ｇ Ｃ（６５マッピング位置）
Ｔ Ｔ Ｇ Ｃ Ｃ(５４マッピング位置)
Ｔ Ｔ Ｇ Ｃ Ｃ Ｔ(２７マッピング位置)

上記の例示の場合、リード配列を参照して２個の塩基を追加した場合、マッピング位置の個数が設定値以下に低下するため、最終断片配列は最初生成された値より２ｂｐ拡張されたＴ Ｔ Ｇ Ｃ Ｃ Ｔになる。一方、前述した他の例と同様に、前記設定値もまた参照配列、リード配列及び断片配列の特性等によって適切に定められ得る値であって、本発明は、特定の設定値にその権利範囲が限定されるものではないことに留意すべきである。

ヒトの塩基配列を対象とした実験で、１００万個のリードで断片配列の長さ１５ｂｐ、シフト間隔４ｂｐで断片配列を生成した後、生成された断片配列を参照配列にマッチングする場合、基準値を５０とする時に総１５，５４７，８５６個の断片配列のうち約７７％の断片配列が５０個以下のマッピングを有するものと示された。すなわち、実験の結果、基準値を５０とする場合、７７％の断片配列はこれをそのまま活用することができ、残りの２３％の断片配列は前述した方法による断片配列の拡張が必要なものと示された。

マッピングヒストグラム生成（段階１１６）
上述した過程を通じて、候補断片配列集合（sub-candidate）が構成されると、原則的にこれらの候補断片配列集合の前記参照配列内でのマッピング位置を利用してリード配列を参照配列にマッピングすることが可能である。しかしながら、この場合、候補断片配列の各マッピング位置の全ての組み合わせを考慮しなければならないため、リード配列のマッピングのための計算の複雑度が非常に高くなる。例えば、候補断片配列集合に含まれる候補断片配列が４個であり、それぞれの候補断片配列の参照配列でのマッピング位置の個数がそれぞれ３、６、２４、４９個である場合、２１，１６８（＝３＊６＊２４＊４９)個の組み合わせ全てを調査しなければならない。本発明では、このような計算の複雑度を低減するためにマッピングヒストグラムを利用することになる。

本発明において、マッピングヒストグラムは、一定のサイズを有する配列（integer array）で、配列の値は、参照配列を同一のサイズを有する複数の区間に分割する時の各区間に対応する。例えば、参照配列を６５５３６（＝２^１６）ｂｐサイズを有する区間に分割する場合、参照配列の０〜６５５３５ｂｐまでの区間は、マッピングヒストグラム（ｈ）の一番目の値であるｈ［０］に対応し、６５５３６〜１３１０７１までの区間は、マッピングヒストグラム（ｈ）の二番目の値であるｈ［１］に対応する。このような方式で参照配列の分割された各区間をマッピングヒストグラムに対応させることができる。

また、マッピングヒストグラムの各値（ｈ［ｉ］）には、対応する参照配列区間での前記候補断片配列の総マッピング長さ（Ａ）が保存され、さらに、該当参照配列区間にマッピングされる候補断片配列のリード配列での位置のうち最大値（Ｂ）が保存され得る。

例えば、図３に示す実施形態において、５３〜６７の断片配列がｈ［０］区間にマッピングされる場合、ｈ［０］のヒストグラム値は（１５、６７）となる(ここで、１５はｈ［０］区間にマッピングされる候補断片配列の総マッピング長さ、６７はマッピングされる候補断片配列のリード配列での最後の位置)。同一の方式で、４９〜６３の断片配列がｈ［１］区間にマッピングされる場合、ｈ［１］のヒストグラム値は（１５、６３）になる。以後、６１〜７５の断片配列がｈ［０］区間にマッピングされる場合、ｈ［０］のヒストグラム値は（２３、７５）にアップデートされ、その理由は次のとおりである。

一番目の値２３：最初にマッピングされた５３〜６７及び以後にマッピングされた６１〜７５のオーバーラップされる区間を考慮した総マッピング長さ

二番目の値７５：マッピングされる断片配列の最後の位置。すなわち、６１〜７５の最後の位置に該当

マッピング対象区間の選択及び全域アラインメント（global assignment）（１１８）
前記のような過程を通じてマッピングヒストグラムが生成されると、マッピングヒストグラムのヒストグラム値(該当区間での候補断片配列の総マッピング長さ)が既設定された基準値（ヒストグラムカット：histogram cut)以上の区間をマッピング対象区間で選択する。

図４は、本発明の一実施形態によるマッピングヒストグラムを利用したマッピング対象区間の選択過程を例示した図である。図４に示すように、参照配列がｈ［０］からｈ［３］まで４個の区間に分割され、各区間のヒストグラム値が下記のように計算されるものと仮定する。
ｈ［０］＝１５、ｈ［１］＝０、ｈ［２］＝２３、ｈ［３］＝１５

この時、前記ヒストグラムカット値が２２に設定される場合、ヒストグラムカット値より大きい区間はｈ［２］に該当する区間となり、本段階では、ｈ［２］に該当する区間をマッピング対象区間で選択するようになる。この時、仮にヒストグラム値がヒストグラムカットより大きい区間が複数の場合には、該当する全区間がマッピング対象区間となり、マッピング対象区間に含まれる複数の区間の全てにおいて全域アラインメントが遂行される。この場合、アラインメント速度を上げるためにマッピング対象区間に含まれる各区間のヒストグラム値を相互に比較し、ヒストグラム値が高い区間から順次に全域アラインメントを遂行することができる。ヒストグラム値が高いということは、マッピングされる断片配列の総長さがより長いという意味であるため、該当区間でリード配列がマッピングされる可能性が高いためである。また、この場合、各区間別のヒストグラム値が同一である場合には、各区間にマッピングされる候補断片配列の個数を計算し、断片配列の個数が多い区間から全域アラインメントを遂行することができる。

このようにマッピング対象区間が選択されると、次に候補断片配列（sub-candidate）のうち該当マッピング対象区間にマッピングされる候補断片配列を最終候補断片配列（candidate）で選定し、選定された最終候補断片配列それぞれのマッピング位置でリード配列に対する全域アラインメントを遂行することによって、リード配列に対するアラインメントを完了する。

例えば、前記図４の実施形態において、ｈ［２］区間にマッピングされる候補断片配列が４９〜６３、５３〜６７、６１〜７５の３個であると仮定すると、前記３個の候補断片配列が最終候補となり、これらの該当区間内のマッピング位置でリード配列の全域アラインメントを遂行するようになる。

一方、前記最終候補断片配列に対する全域アラインメント時には、全域アラインメントに所要する時間を減らすために、全域アラインメントを一回遂行した参照配列での位置を記憶し、それと近い位置では、繰り返し全域アラインメントが数回遂行されないようにする。具体的に、本段階では、前記マッピング対象区間を複数の小区間に分割した後、全域アラインメントが遂行された小区間の場合、これを記録しておくように構成される。以後、該当小区間に対する全域アラインメント時には、前記記録された情報を利用して該当小区間で全域アラインメントが既遂行されたか否かを判断し、前記判断の結果、全域アラインメントが既遂行されていない場合にのみ全域アラインメントを遂行するようになる。

これを例を挙げて説明すると、図５のとおりである。図５に示すように、マッピング対象区間が５個の小区間に分割され、前記３個の最終候補のうち４９〜６３及び５３〜６７は二番目の小区間に、６１〜７５は４番目の小区間にマッピングされるものと仮定する。この場合、仮に４９〜６３の断片配列に対して２番目の小区間で全域アラインメントが遂行される場合、その結果に関係なく同一の小区間に属する５３〜６７に対しては全域アラインメントが遂行されず、これはその反対の場合も同様である。従って、示された実施形態の場合、全域アラインメントは、４９〜６３／６１〜７５又は５３〜６７／６１〜７５の組み合わせに対してのみ実施するようになる。本発明のように参照配列の全体ではなくマッピング対象区間内でのみ全域アラインメントが遂行されるとしても、全域アラインメントのためには非常に時間が掛かるため、このような過程を通じる場合、全域アラインメントに所要する時間を減らすことができる。

ヒストグラムカット（histogram cut）の計算
前記実施形態において、ヒストグラムカットは、次のような方式で計算することができる。

先ず、ｆを断片配列のサイズ、ｓを断片配列を生成するためのリード配列内での移動間隔、Ｌをリード配列の長さ、ｅをリード配列で許容される最大エラーの個数、Ｈをヒストグラムカットとすると、リード配列でエラーの影響を受けない領域の長さＴは、下記の数式のように求めることができる。

［数式１］
Ｔ＝Ｌ−ｆ*ｅ−ｓ

この時、Ｌ及びｅは、本発明を実施する際、既に決定されている値であるため、ｆ及びｓの値によってＴが決定される。すなわち、ｆ及びｓの値を如何に変化させるかによってアルゴリズムの性能が変化するようになる。

先ず、Ｈ値を決定する時には、下記二つの条件を考慮する。このうち、必須条件は必ず満たさなければならず、追加条件は可能な場合に考慮する。

−必須條件：マッピングの基本単位が断片配列であるため、ヒストグラムカットがいくら小さいとしても少なくともオーバーラップされる２個以上の断片配列を含むことができるサイズでなければならない。仮に、図２のようにｆ＝１５、ｓ＝４である場合、オーバーラップされる２個の断片配列の最小長さは１５＋４＝１９となるため、少なくともＨ値は１９以上であるべきである。また、前記Ｈ値は、少なくとも２個の断片配列が含まれるように設定すべきであるため、最小限ｆ＋ｓよりは大きいか、或いは同一でなければならない。後述するように、ｆ値は最小限１５以上でなければならないため、ｓ値をその最小値である１と仮定する場合、Ｈは最小限１６（＝１５＋１）以上の値となる。

−追加条件：理想的な状況を仮定する時、Ｈ＝Ｔと定め、Ｔ以上のシーケンスがマッピングされたヒストグラムを検索すれば、与えられたエラーに対する全てのマッピングを検索することができる。しかしながら、前述したように、参照配列自体に反復が多い場合、状況によって断片配列の長さを拡張しなければならない場合が生じ得る。従って、これを考慮してＨ値を定める時には、Ｔより若干小さいＴ−ｓを用いることがマッピング率の側面で有利である。仮にＨ＝Ｔと仮定した場合、Ｈ＝Ｔ−ｆ＊ｅ−ｓになり、このうちｅを最小値である１と仮定した場合(ｅが０である場合は、参照配列と一致・整合する場合であるため、前述した段階１０４でマッピングが完了する)、Ｈ＝Ｔ−ｆ−ｓとなる。この値がヒストグラム値の最大値となる。仮に、Ｌ＝７５ｂｐ、ｆ＝１５ｂｐ、ｓ＝１と仮定した場合、Ｈの最大値は７５−１５−１＝５９となる。

つまり、前記Ｈ値は、次の範囲を満たさなければならない。

［数式２］
ｆ＋ｓ＜＝Ｈ＜＝Ｔ−（ｆ＋ｓ）

次に、ｆ値は、下記二つの条件を満たす値のうち大きい値を選択する。これもまた必須条件は、必ず満たさなければならず、追加条件は可能な場合に考慮する。

−必須條件：ｆは１５以上でなければならず、その理由として、断片配列の長さが１４以下である場合、参照配列内でのマッピング位置の個数が急激に増加するためである。

下記表１は、断片配列の長さによるヒトのゲノム内での断片配列の平均発現頻度を示したものである。

前記表からわかるように、断片配列の長さが１４以下である場合には、断片配列別の頻度が１０以上であるが、１５である場合には、３以下に減少する。すなわち、断片配列の長さを１５以上で構成する場合、１４以下で構成する場合と比較して断片配列の反復を大幅に減少させることができる。

−追加条件：ｆ≦Ｌ／（ｅ＋２）を満たさなければならず、これはＴの長さを断片配列２個のサイズ以上に保障するためである。

例えば、Ｌ＝１００、ｅ＝４である時、ｆは１６以下の値を有するべきである。

上記条件を併せて、ｆ、ｓ、及びＨを決定する方法をまとめたものを次に示す。
−ｓは４に固定した後、ｆとＨを決定する。
−１５≦ｆ≦Ｌ／（ｅ＋２）範囲内で最も大きい値をｆに決定する(ただし、必ずｆ≧１５)。
−Ｈは、下記式を利用して決定する。

［数式３］
Ｈ＝Ｌ−ｆ＊ｅ−２ｓ又はＨ＝ｆ＋ｓから計算される値のうち大きい値
(ここで、Ｈは基準値、Ｌはリード配列の長さ、ｆは断片配列の長さ、ｅはリード配列の最大エラー個数、ｓは各断片配列の移動間隔である)

例１）Ｌ＝７５、ｅ＝３である時、ｆ＝１５〜１５であるため１５、ｓ＝４、Ｈ＝７５−３＊１５−２＊４＝２２となる。
例２)Ｌ＝１００、ｅ＝４である時、ｆ＝１５〜１６であるため１６、ｓ＝４、Ｈ＝１００−４＊１６−２＊４＝３６−８＝２８となる。
例３)Ｌ＝７５、ｅ＝４である時、ｆ＝１２〜１５であるが、１５以上でなければならないため１５、ｓ＝４、Ｈ＝７５−４＊１５−２＊４＝１５−８＝７であるが、ｆ＋ｓ＝１９であるため、結果的にＨ＝１９となる。

図６は、本発明の一実施形態による塩基配列アラインメントシステム６００のブロック図である。本発明の一実施形態による塩基配列アラインメントシステム６００は前述した塩基配列アラインメント方法を遂行するための装置であって、断片配列生成部６０２と、フィルタリング部６０４と、マッピング長さ計算部６０６と、アラインメント部６０８と、断片配列拡張部６１０とを備える。

断片配列生成部６０２は、ゲノムシーケンサから得られたリード配列から複数の断片（fragment）配列を生成する。前述したように、断片配列生成部６０２は、前記リード配列の一番目の塩基から設定された間隔ほど移動しながら設定されたサイズほど前記リード配列の値を読み取ることによって、前記断片配列を生成するようになる。

フィルタリング部６０４は、生成された前記複数の断片配列のうち前記参照配列とマッチングする断片配列のみを含む候補断片配列集合を構成する。この時、前記参照配列とマッチングする断片配列は、前記参照配列との一致・整合（exact matching）の結果、不一致となる塩基の数が、設定された個数以下である断片配列を意味する。

マッピング長さ計算部６０６は、前記参照配列を複数の区間に分割し、前記各区間別に前記候補断片配列のマッピング位置及び各区間別に前記候補断片配列の総マッピング長さを計算する。

アラインメント部６０８は、マッピング長さ計算部６０６によって分割された区間のうち計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメントを遂行する。具体的に、アラインメント部６０８は、前記候補断片配列のうち前記選択された区間にマッピングされる候補断片配列の前記参照配列内でのマッピング位置に基づいて、前記リード配列に対する全域アラインメントを遂行する。

また、アラインメント部６０８は、前記選択された区間(マッピング対象区間)を複数の小区間に分割し、前記全域アラインメントを遂行しようとする前記参照配列内の位置が属する小区間で全域アラインメントが既遂行されたか否かを判断し、前記判断の結果、全域アラインメントが既遂行されていない場合にのみ前記全域アラインメントを遂行することによって、不要な全域アラインメントの回数を減少させるように構成され得る。

断片配列拡張部６１０は、前記フィルタリング部６０４で生成された前記候補断片配列それぞれの前記参照配列でのマッピング位置の個数を計算し、計算された前記マッピング位置の個数が設定された値を超える断片配列を選択し、前記参照配列でのマッピング位置の個数が前記設定された値以下になるまで選択された断片配列のサイズを拡張する。この時、断片配列拡張部６１０は、前記選択された断片配列の最初又は最後の部分に、該当位置に対応する前記リード配列での塩基を追加することによって、前記拡張を遂行する。

一方、本発明の実施形態は、本明細書において記述した方法をコンピュータ上において遂行するためのプログラムを含むコンピュータ判読可能記録媒体を備えることができる。前記コンピュータ判読可能記録媒体は、プログラム命令、ローカルデータファイル、ローカルデータ構造等を単独で又は組み合わせて備えることができる。前記媒体は、本発明のために特別に設計され、構成されたものであるか、或いはコンピュータソフトウェア分野において通常の知識を有する者に公知となって使用可能なものであり得る。コンピュータ判読可能記録媒体の例としては、ハードディスク、フロッピーディスク（登録商標）及び磁気テープのような磁気媒体、ＣＤ−ＲＯＭ、ＤＶＤのような光記録媒体、フロッピーディスクのような磁気−光媒体、及びＲＯＭ、ＲＡＭ、フラッシュメモリ等のようなプログラム命令を保存して遂行するように特別に構成されたハードウェア装置が含まれる。プログラム命令の例としては、コンパイラーによって作成されるような機械語コードのみではなくインタプリター等を用いてコンピュータによって実行され得る高級言語コードを含むことができる。

以上において、代表的な実施形態を通じて本発明について詳細に説明したが、本発明が属する技術分野において通常の知識を有する者は、詳述した実施形態について本発明の範疇から外れない範囲内で様々な変形が可能であることを理解しなければならない。

従って、本発明の権利範囲は説明された実施形態に限定して定めてはならず、後述する特許請求の範囲のみならず、この特許請求の範囲と均等なものなどによって定めなければならない。

６００：塩基配列アラインメントシステム
６０２：断片配列生成部
６０４：フィルタリング部
６０６：マッピング長さ計算部
６０８：アラインメント部
６１０：断片配列拡張部

Claims (19)

1. リード配列から複数の断片（fragment）配列を生成する断片配列生成部と、生成された前記複数の断片配列から候補断片配列集合を構成するフィルタリング部と、前記候補断片配列それぞれの参照配列でのマッピング位置の個数を計算し、計算された前記マッピング位置の個数が、設定された値を超える断片配列を選択し、前記参照配列でのマッピング位置の個数が前記設定された値以下になるまで選択された断片配列のサイズを拡張する断片配列拡張部と、前記参照配列を複数の区間に分割し、各区間別に前記候補断片配列の総マッピング長さを計算するマッピング長さ計算部と、計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行するアラインメント部とを備える塩基配列アラインメントシステム。
2. 前記断片配列生成部は、前記リード配列の一番目の塩基から設定された間隔ほど移動しながら設定されたサイズほど前記リード配列の値を読み取ることによって、前記断片配列を生成する、請求項１に記載の塩基配列アラインメントシステム。
3. 前記フィルタリング部は、前記参照配列との一致・整合（exact matching）の結果、不一致となる塩基の数が、設定された個数以下である断片配列を前記候補断片配列集合に含まれる、請求項１に記載の塩基配列アラインメントシステム。
4. 前記断片配列拡張部は、前記選択された断片配列の最初又は最後の部分に、該当位置に対応する前記リード配列での塩基を追加する、請求項１に記載の塩基配列アラインメントシステム。
5. 前記アラインメント部は、前記候補断片配列の中から前記選択された区間にマッピングされる候補断片配列を選択し、選択されたそれぞれの候補断片配列の前記参照配列内でのマッピング位置で前記リード配列に対する全域アラインメントを遂行する、請求項１に記載の塩基配列アラインメントシステム。
6. 前記アラインメント部は、前記選択された区間を複数の小区間に分割し、前記全域アラインメントを遂行しようとする前記参照配列内の位置が属する小区間で全域アラインメントが既遂行されたか否かを判断し、前記判断の結果、全域アラインメントが既遂行されていない場合にのみ前記全域アラインメントを遂行する、請求項５に記載の塩基配列アラインメントシステム。
7. 前記基準値は、次の数式
   Ｈ＝Ｌ−ｆ＊ｅ−２ｓ
   (ここで、Ｈは基準値、Ｌはリード配列の長さ、ｆは断片配列の長さ、ｅはリード配列の最大エラー個数、ｓは各断片配列の移動間隔である)、又は
   Ｈ＝ｆ＋ｓ
   より計算される値のうち大きい値である、請求項１に記載の塩基配列アラインメントシステム。
8. 前記基準値は、次の数式
   ｆ＋ｓ＜＝Ｈ＜＝Ｔ−（ｆ＋ｓ）
   を満たす、請求項７に記載の塩基配列アラインメントシステム。
9. 前記基準値は、１６以上５９以下である、請求項１に記載の塩基配列アラインメントシステム。
10. リード（read）配列を参照配列にアラインメントするための方法であって、断片配列生成部で、前記リード配列から複数の断片（fragment）配列を生成する段階と、フィルタリング部で、生成された前記複数の断片配列から候補断片配列集合を構成する段階と、断片配列拡張部で、生成された前記候補断片配列それぞれの前記参照配列でのマッピング位置の個数を計算する段階と、前記断片配列拡張部で、計算された前記マッピング位置の個数が設定された値を超える断片配列を選択する段階と、前記断片配列拡張部で、前記参照配列でのマッピング位置の個数が前記設定された値以下になるまで選択された断片配列のサイズを拡張する段階と、マッピング長さ計算部で、前記参照配列を複数の区間に分割し、各区間別に前記候補断片配列の総マッピング長さを計算する段階と、アラインメント部で、計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行する段階とを含む塩基配列アラインメント方法。
11. 前記断片配列を生成する段階は、前記リード配列の一番目の塩基から設定された間隔ほど移動しながら設定されたサイズほど前記リード配列の値を読み取ることによって、前記断片配列を生成する、請求項１０に記載の塩基配列アラインメント方法。
12. 前記候補断片配列集合を構成する段階は、前記参照配列との一致・整合（exact matching）の結果、不一致となる塩基の数が設定された個数以下である断片配列を前記候補断片配列集合に含ませる、請求項１０に記載の塩基配列アラインメント方法。
13. 前記選択された断片配列のサイズを拡張する段階は、前記選択された断片配列の最初又は最後の部分に、該当位置に対応する前記リード配列での塩基を追加するように構成される、請求項１０に記載の塩基配列アラインメント方法。
14. 前記全域アラインメントを遂行する段階は、前記候補断片配列の中から前記選択された区間にマッピングされる候補断片配列を選択し、選択されたそれぞれの候補断片配列の前記参照配列内でのマッピング位置で前記リード配列に対する全域アラインメントを遂行する、請求項１０に記載の塩基配列アラインメント方法。
15. 前記全域アラインメントを遂行する段階は、前記選択された区間を複数の小区間に分割する段階と、前記全域アラインメントを遂行しようとする前記参照配列内の位置が属する小区間で全域アラインメントが既遂行されたか否かを判断する段階をさらに含み、前記判断の結果、全域アラインメントが既遂行されていない場合にのみ前記全域アラインメントを遂行する、請求項１４に記載の塩基配列アラインメント方法。
16. 前記基準値は、次の数式
    Ｈ＝Ｌ−ｆ＊ｅ−２ｓ
    (ここで、Ｈは基準値、Ｌはリード配列の長さ、ｆは断片配列の長さ、ｅはリード配列の最大エラー個数、ｓは各断片配列の移動間隔である)、又は
    Ｈ＝ｆ＋ｓ
    より計算される値のうち大きい値である、請求項１０に記載の塩基配列アラインメント方法。
17. 前記基準値は、次の数式
    ｆ＋ｓ＜＝Ｈ＜＝Ｔ−（ｆ＋ｓ）
    を満たす、請求項１６に記載の塩基配列アラインメント方法。
18. 前記基準値は、１６以上５９以下である、請求項１０に記載の塩基配列アラインメント方法。
19. リード配列から複数の断片（fragment）配列を生成する断片配列生成部と、生成された前記複数の断片配列から候補断片配列集合を構成するフィルタリング部と、参照配列を複数の区間に分割し、各区間別に前記候補断片配列の総マッピング長さを計算するマッピング長さ計算部と、計算された前記総マッピング長さが基準値以上の区間を選択し、選択された区間に対して前記リード配列に対する全域アラインメント（global alignment）を遂行するアラインメント部とを備える塩基配列アラインメントシステム。