JP5985040B2 - データ解析装置、及びその方法 - Google Patents

Description

本発明はDNA配列のデータ解析装置に係り、特に超並列型のDNAシーケンサ装置から得られるDNA配列データの解析技術に関するものである。

癌や生活習慣病や遺伝病などに対して、いわゆる個別化医療として、患者個人に適した治療法を選択したり、予後の予測を行うために、患者個人の遺伝的背景を調べることが必要になる。そのため、ゲノムやトランスクリプトーム（転写産物）などのDNA(deoxyribonucleic acid) 配列解析が行われる。そのとき用いられるDNAシーケンサ装置では断片化された短いDNA配列しか得られない。そのため、長大な参照ゲノム配列と比較して、得られた断片配列がゲノムのどの部分に由来するかを調べ、さらに、そこに含まれる一塩基変異（SNP， Single Nucleotide Polymorphism）や挿入・欠失などの変異を調べるデータ処理が必要となる。そのようなデータ処理は一般にマッピング処理と呼ばれる。

いわゆる次世代型DNAシーケンサとよばれる超並列型のDNAシーケンサでは、1回の計測で比較的短い100塩基程度の長さの断片配列（リード）が数億本以上も得られる。また、ヒトの場合、参照ゲノム配列の長さは約３ギガ塩基（30億塩基）にも及ぶ。マッピング処理では、これらのリード配列を1本ずつ参照ゲノム配列と比較して、対応する位置を特定し、そこに含まれる変異を同定する。これには非常に大きな計算コストが必要となるため、専用の効率的なアルゴリズムが開発され利用されている。代表的な方法は、参照ゲノム配列をBurrows-Wheeler変換（BWT, Burrows-Wheeler Transformation）（非特許文献１）によりデータベース化しておき、リード配列内の短い塩基配列を検索キーとして検索し、マッチした領域の前後で、シーケンシング・エラーや変異の可能性を考慮して、アラインメントを行う（非特許文献２）。

一般に、次世代型DNAシーケンサでは1%程度の読み取りエラーが生じ、また、長大なゲノム領域の中には、類似した配列が多数散在している。そのため、リード1本ごとのマッピング結果には誤りが生じる可能性がある。例えば、あるリード配列に対して、完全一致する領域は参照ゲノム配列内に無いが、少数のシーケンシング・エラーを仮定すれば対応するゲノム領域が複数か所見つかる場合がある。このような場合、どの領域を選択するかには任意性があり、その判断はマッピング処理のヒューリスティクスに任されている。そこで、変異解析を正確に行うためには、後続の処理、すなわち下流の処理において、多数のリードのマッピング結果を比較して多数決をとる、再マッピング処理が行われる（（非特許文献３）。そのため、全ゲノム解析を行う場合は、通常、ゲノム全体を数十倍カバーできるほどの配列量（数十ギガ塩基以上）がシーケンシングされる。また、マッピング先に任意性があるときにマッピング処理に依存したバイアスがかかることが懸念されるため、複数種類のマッピング・ツールによる結果を比較して、そのようなバイアスが生じていないかを確認することも行われる。なお、以上の技術に関連する特許文献としては、例えば特許文献１がある。

特開２００３−３３０９３４号公報

M. Burrows and D. Wheeler: A block-sorting lossless data compression algorithm. Technical Report 124, Digital Equipment Corporation, 1994. Li H. and Durbin R. (2009) Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 25:1754-60. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA (2010). The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20:1297-303. Mantaci, S., Restivo, A. ; Sciortino, M: "An extension of the Burrows Wheeler transform to k words." Data Compression Conference, 2005. Proceedings. DCC 2005. Markus J. Bauer, Anthony J. Cox, Giovanna Rosone: "Light-weight BWT Construction for Very Large String Collections," Combinatorial Pattern Matching, Lecture Notes in Computer ScienceVolume 6661, 2011, pp 219-231 Paolo Ferragina, Travis Gagie, Giovanni Manzini: "Light-weight Data Indexing and Compression in External Memory," Algorithmica, July 2012, Volume 63, Issue 3, pp 707-730. Kimura K, Suzuki Y, Sugano S, Koike A: "Computation of rank and select functions on hierarchical binary string and its application to genome mapping problems for short-read DNA sequences," J Comput Biol. 2009 Nov;16(11):1601-13. Ge Nong; Sen Zhang; Wai Hong Chan; , "Linear Suffix Array Construction by Almost Pure Induced-Sorting," Data Compression Conference, 2009. DCC '09. , vol., no., pp.193-202, 16-18 March 2009

上述した下流の処理では高い精度を得るために大きな計算コストが必要となるため、DNAシーケンサから得られた全てのリードを一括して取り扱うことは出来ない。そこで、効率的なアルゴリズムを採用したマッピング処理の結果を利用して、注目する遺伝子領域に、領域に由来する可能性が高いと思われるリード配列を選び出し、それらに対して下流の解析が行われる。

また、他方では、比較的リード長が長く（500塩基程度以上）リード本数が比較的少ない（100本程度）のキャピラリー型DNAシーケンサで多数回解析して得られたリード配列データをデータベース化し、ゲノム内の注目した遺伝子領域の配列をクエリーとしてホモロジー検索を行い、得られたリード配列をマルチプル・アラインメントしてバクテリアの同定を行う方法が知られている（特許文献１）。しかしながら、ヒトゲノムを次世代型DNAシーケンサで解析したときのデータ量は数十ギガ塩基以上に達するほど巨大であるため、実用に耐える計算時間でホモロジー検索することは不可能となる。

ヒトの全ゲノム解析を行う場合、総シーケンス量が数十ギガ塩基にも達するほどリード配列データ量が多いため、効率的なアルゴリズムを採用しているにも関わらず、マッピング処理の計算コストは大きく、この計算コストを下げることが課題となる。

また、シーケンシング・エラーの扱い方によりマッピング先に任意性がある場合、マッピング処理でヒューリスティクスを用いてその中からマッピング先を選ぶことは、マッピング処理に依存したバイアスが生じることを意味する。このようなヒューリスティクな判断を避けて、どのマッピング先の可能性も等価に扱われるような、中立的な処理方法を提供することが課題となる。

現在、ヒトの参照ゲノム配列は、ヒト白血球型抗原（Human Leukocyte Antigen： HLA）領域などの特殊な領域を除けば、一通りしかない。しかし、参照ゲノム配列を複数通り用意しておけば、患者の属する人種グループに適したものをその中から選択して、より精密に変異解析を行える可能性がある。マッピング処理は、全リード配列と参照ゲノム配列との組合せに対して行われるため、参照ゲノムを変更することは、全てのマッピング処理をやり直すことを意味する。そこで、全リード配列と参照ゲノム配列をそれぞれ独立に処理し、その組合せを変えて解析する時の計算コストの増加を抑えることができるようにすることが課題となる。

本願発明の目的は、上記の課題の少なくとも一つを解決し、マッピング処理の計算コストを下げることが可能な、或いは中立的な処理が可能なデータ解析装置、及びその方法を提供することにある。

上記の目的を達成するため、本発明においては、処理部と記憶部を備えたデータ解析装置であって、記憶部は、ゲノム配列データがデータベース化されたゲノム配列データベースと、リード配列データがデータベース化されたリード配列データベースを記憶し、処理部は、特定された解析対象のゲノム領域の配列に基づいて、検索用の塩基配列であるキー配列を選択し、リード配列データベースの中にあるキー配列の深度を求め、リード配列データベースの中にあるキー配列を含むリード配列データを抽出し、抽出されたリード配列データとゲノム領域の配列とを比較してデータ解析を行う構成のデータ解析装置を提供する。

また、上記の目的を達成するため、本発明においては、データ解析装置の処理部によるデータ解析方法であって、ゲノム配列データが検索可能な形式にデータベース化されたゲノム配列データベースと、リード配列データが検索可能な形式にデータベース化されたリード配列データベースを用い、特定された解析対象のゲノム領域の配列に基づいて検索用の塩基配列であるキー配列を選択し、リード配列データベースの中にあるキー配列の深度を求め、リード配列データベースの中にあるキー配列を含むリード配列データを抽出し、抽出されたリード配列データとゲノム領域の配列とを比較してデータ解析を行うデータ解析方法を提供する。

本発明の解析装置、及びその方法により、計算コストを抑えることができる。また、処理法に依存するバイアスは発生せず、中立的な処理が可能である。

実施例１に係る、変異解析を行うための処理手順を示すフローチャート図である。 実施例１に係る、配列比較により変異の有無を判定する方法を説明するための説明図である。 実施例１に係る、MLUと深度に基づき変異の有無を推定する方法を示すフローチャート図である。 実施例１に係る、長さが不揃いの多数のリード配列からなるデータに対して、一般化されたBurrows-Wheeler変換（BWT）を定義する説明図である。 実施例１に係る、SLCP（sorted list of cyclic permutations）を用いて文字列wの深度D(w)を求める方法を示す説明図である。 実施例１に係る、BWTを利用する際に用いる補助関数を説明する説明図である。 実施例１に係る、検索キー配列に対する深度を計算するためのフローチャート図である。 実施例１に係る、キー配列を含む全てのリード配列を、キー配列の左方にある塩基により分類して該当するリード配列の本数を求める方法を示したフローチャート図である。 実施例１に係る、BWTの計算方法を示したフローチャート図である。 実施例１に係る、旧文字列から新文字列への転記方法を示す説明図である。 実施例１に係る、ゲノム配列DBとリード配列DBの構成を示す説明図である。 実施例１に係る、ゲノム座標xにおけるMLUの値L(x)の計算方法を示すフローチャート図である。 実施例２に係る、注目する遺伝子領域の中から、スプライシングが生じている可能性が高い位置をMLUと深度に基づいて推定し、さらに、推定された位置において、配列比較に基づいてスプライシングの有無を判定する方法を示した説明図である。 実施例２に係る、順方向に探索してスプライシング有りと推定された位置において、配列比較によりスプライシングの有無を判定する処理を説明するための説明図である。 実施例２に係る、リード配列の共通配列Sとゲノム配列を比較してスプライシングの有無を判定する方法を示すフローチャート図である。 実施例２に係る、順方向に探索してゲノム座標xにおいて、MLUと深度に基づきスプライシングの有無を推定する方法を示すフローチャート図である。 実施例に係る、解析装置の内部の一構成例を示すブロック図である。 実施例のデータ解析方法の処理に係る、式１〜式３を示す図である。 実施例のデータ解析方法の処理に係る、式４〜式７を示す図である。

以下、本発明の種々の実施の形態を図面に従い説明するが、それに先立ち本発明の好適な態様の概要を説明する。本発明のデータ解析装置、及びその方法にあっては、全てのリード配列とそれら全ての巡回置換（cyclic permutation）、或いは接尾辞（suffix）を辞書式順序（lexicographic order）にソートしてデータベース化し、任意の短い塩基配列をキーに検索して、キー配列を含むリードの本数を即答でき、また、全リード配列の中からキー配列を含む全てのリード配列を即座に抽出できるようにする。

そして、参照ゲノム配列側では、各塩基位置からスタートする部分配列が何塩基の長さに達すれば参照ゲノム配列内で、相補鎖も考慮して一意的となるかを事前に調べてデータベース化し、任意の塩基位置でそのような一意性を保証できる最小長（MLU, minimum length for uniqueness）を即答できるようにする。このような計算は参照ゲノム配列データ単独で行える。そのため、参照ゲノム配列に対して一度計算しておけば、どのようなリード配列データに対しても再利用できる。

更に、下流の変異解析処理では、注目する遺伝子領域内部を１塩基ずつスキャンして、各塩基位置でMLUを参照ゲノムデータベースに問合せ、さらに長さMLUのゲノム部分配列を参照ゲノムデータベースに問合せて取得してから、これを検索キーとしてリード配列データベースに問合せて、キー配列を含むリード配列の本数（深度, depth）を得る。周辺と比較して深度の値が大きく低下する塩基位置が見つかれば、そこに変異が含まれる可能性が高いと推定する。

変異が含まれる可能性と高いと推定された場所では、再度、その近辺で深度が低下していない別の位置でMLUを参照ゲノムデータベースに問合せ、さらに長さMLUのゲノム塩基配列を参照ゲノムデータベースに問合せて取得してから、これを検索キーとしてリード配列データベースに問合せて、キー配列を含む全てのリード配列を抽出する。これにより、その近辺に由来する可能性が高いリード配列を集めることができる。これらを対象にして、詳細な変異解析処理を行う。

実施例１は、注目する遺伝子領域の中から、変異が含まれる可能性が高い塩基位置をMLUと深度に基づいて推定し、更に推定された位置において、配列比較に基づいて変異の有無を判定する解析装置、及び方法の実施例である。

図１は、本実施例により変異解析を行うための処理手順を示すフローチャートである。なお、各実施例に係る解析を実現する解析装置は、通常の計算機の構成を有するサーバ等のコンピュータで実現される。

図１７に、本実施例を含め全ての実施例に係る解析装置の一構成例を示した。同図において、解析装置1700は、処理部である中央処理部（CPU：Central Processing Unit）1701、プログラム等が記憶される記憶部であるメモリ1702、操作のためのGUI(Graphical User Interface)や、解析結果等を表示する表示部1703、データベース(DB)等を記憶する記憶部として機能するハードディスクドライブ（HDD）1704、パラメータ入力などを行うキーボード等の入力部1705、インターネット等に接続するためのネットワークインタフェース (NIF)1706がバス1707に接続された構成を備えている。HDD1704に記憶されるデータベース（DB）は、解析装置1700の外部に設置された記憶装置に記憶することもできるし、ネットワークを介してデータセンタ等に記憶することも可能である。以下の実施例で説明される各種のフローチャートは、CPU1701のプログラム実行などで実現される。

さて、図１に示したフローチャートにおいて、超並列型DNAシーケンサ102によりDNAサンプル101を解析して、多数の短い塩基配列からなるリード配列データ103を得る。これらに対してリード配列データベース(DB)化処理104を行い、リード配列DB105を得る。リード配列DB化処理104では、全てのリード配列とそれら全ての巡回置換（cyclic permutation）、或いは接尾辞（suffix）を辞書式順序にソートする。このDB化により、任意の短い塩基配列をキーに検索して、キー配列を含むリードの本数（深度）を即答でき、また、全リード配列の中からキー配列を含む全てのリード配列を即座に抽出できるようにする。

参照ゲノム配列データ106は、ゲノム配列DB化処理107によりデータベース化して、ゲノム配列DB108を構築する。参照ゲノムDB化処理では、参照ゲノム配列の各塩基位置で、そこからスタートする部分配列が何塩基の長さに達すれば参照ゲノム配列内で、相補鎖も考慮して一意となるかを調べて記憶する。DB化しておくことにより、任意の塩基位置でそのような一意性を保証できる最小長（MLU）を即答できるようにする。また、塩基配列をそのまま座標順に記憶し、任意に指定された座標範囲の塩基配列を即答できるようにする。

座標xにおけるMLU、即ち、xからスタートする部分配列がゲノム配列内で一意となる最小の長さをL(x)と表す。また、座標xからスタートする長さL(x)のゲノム部分配列をキー配列としたときのリード配列の深度をxにおける深度とよびD(x)と表す。以下、図１８の式１で示した記法を導入する。図１８の式１の記法を用いると、ゲノム座標が増大する方向（順方向、forward）に探索（スキャン）を行う場合、MLUと深度は、図１８の式２のように計算できる。または、ゲノム座標が減少する方向（逆方向, backward）に探索を行う場合、MLUと深度は、図１８の式３のように計算できる。

図１のフローチャートにおいて、注目する遺伝子領域の範囲を規定する座標(Xmin, Xmax)を入力する(109)。Xmin以上Xmax以下の各xに対して、ゲノム配列DBに問合せてMLU、即ち、L(x)の値を得る(110)。また、リード配列DBに問合せて、各xにおける深度D(x)を得る (111)。

そして、変異の有無推定基準となるパラメータを入力し(112)、xを注目する領域の左端Xminとして(113)、以下の繰り返し処理を開始する。MLUと深度に基づき、xにおける変異の有無を推定する(114)。変異無と推定されるか、または、推定不能となった場合は、直ちにxの値をx＋1に更新する(120)。そうでない場合は、xの近くの別の塩基位置yを深度に基づいて一つ選択し、yからスタートして長さL(y)のゲノム部分配列をゲノム配列DBに問合せて取得する(115)。これをキー配列として、これを含む全てのリード配列を、リード配列DBに問合せて取得する(116)。これらのリード配列と参照ゲノム配列を詳細に比較して変異解析を行い(117)、変異が見つかった場合は、その結果を端末(122)や記憶装置(123)に出力する(119)。その後、xをx＋1に更新する(120)。xが領域の右端を超えていなければ、処理を繰り返す(121)。そうでなければ、処理を終了する。

図２は、上述した実施例１の解析法において、順方向に探索して、変異有りと推定された位置xにおいて、配列比較により変異の有無を判定する方法を説明するための説明図である。逆方向に探索する場合も同様である。図２中のグラフ203は、横軸201にゲノム位置座標、縦軸202にカウント数をとり、位置座標xのおける深度D(x)をプロットしたものである。SNPがゲノムの位置座標zにある場合、204の範囲に示すように、位置xからスタートする長さL(x)のゲノム部分配列がzを通るとき、深度D(x)の値は周辺に比べて大きく低下する。深度が低下する位置は、順方向探索では、変異がある位置の逆方向側にある。そこで、深度が低下した位置xの順方向側の近隣の位置yで深度D(y)が低下していないような任意の位置yを一つとる。yからスタートして長さL(y)のゲノム部分配列をゲノム配列DBに問合せて取得し、その配列をキー配列としてリード配列DBに問合せて、キー配列を含むリード配列を全て集める。

206は集められたリード配列を表し、キー配列に該当する部分を下線で示した。205はx周辺のゲノム配列を表し、キー配列に該当する部分を下線で示した。x周辺のゲノム配列205は、ゲノム配列DBに問合せて得られる。集められたリード配列206は大量にある全リード配列データのごく一部なので、これらをx周辺のゲノム配列205と比較する計算コストは抑えられる。これらの配列をキー配列を基準に揃えてアラインメントし、yの逆方向側の近辺で変異を探す。変異が見つかった場合（例ではzの位置で太字で示すようなGからAへの一塩基変異がある）、変異有りと判定する。そうでない場合は、変異無しと判定する。

図３は、本実施例で、塩基位置xにおいて、MLUと深度に基づき変異の有無を推定する方法を示すフローチャート図である。x自身を含め、xの周辺の各塩基位置yでは、ゲノム配列DBとリード配列DBへの問合せの結果として、MLUと深度の値、即ち、L(y)と D(y)の値を取得済みとする。また、d1, d2, h1, h2,h3, l1は推定基準パラメータであり、これらの値は処理に先だって入力される（112）。位置xにおいて深度D(x)が十分大きければ変異無しと推定する(301)。xの周辺でMLUが大きくなることがあれば推定不能と推定する(302)。xの周辺で平均深度が十分大きくなければ、推定不能とする(303)。周辺と比較して、xで深度が大きく低下していれば変異有りと推定し、そうでなければ、変異無しと推定する(304)。

図４は、本実施例において、長さが不揃いの多数のリード配列からなるデータに対して、一般化されたBurrows-Wheeler変換（BWT）を定義する説明図である。本来のBWTは1本の文字列に対して定義されており（上述した非特許文献１参照）、複数の文字列に対して一般化した定義は幾通りか知られている（本実施例の説明の末尾に記載した非特許文献４、非特許文献５、非特許文献６参照）。ここでは、長さが不揃いの多数の文字列に対して、次のように一般化して定義する。図4では、簡単のため、2本のリード配列に対する定義を説明するが、リード配列の本数が何本に増えても同様に定義できる。

同図の401と402を、対象とする2本のリード配列とする。それらの末尾に区切り文字$を付加した文字列を、403と404とする。これらの文字列の全ての巡回置換（cyclic permutation）（巡回シフト、cyclic shift）として得られる文字列のリストを、405と406とする。これらのリストを統合したのち、辞書式順序（lexicographic order）にソートして得られる文字列のリストを407とする。但し、アルファベットの比較順位は、$＜A＜C＜G＜T＜Nとする。ここで、Nの順位をTよりも後ろにしたのは、NがA, C, G, Tの何れかの塩基を表す特殊な文字だからである。また、２つの文字列を先頭文字から順に比較する際、同じ文字位置に区切り文字$が現れた場合は、それ以降の比較は行わず、それらの順番は任意とする。これらのソート済みの文字列の末尾の文字をソート順に並べて得られる文字列を408とする。ここでは、対応関係を分かり易く示すために、文字列408を縦書きで示している。リード配列データ（401と402）に対して、それらの末尾に区切り文字$を付加し、それらの全ての巡回置換として得られる文字列のリストを辞書式順序にソートして得られる文字列リスト（407）の末尾の文字をソート順に並べて得られる文字列408を、リード配列データ（401と402）のBWTと定義する。また、このとき途中で得られる407をソート済み巡回置換文字列リスト（SLCP, sorted list of cyclic permutations）とよぶ。

図５は、本実施例において、SLCPを用いて文字列wの深度D(w)を求める方法を示す説明図である。501はSLCPを表す。501の各要素の末尾の文字を連結したものはBWTであるが、煩雑さをさけるため図中には示していない。SLCPは辞書式順序にソート済みゆえ、501の先頭の文字を連結したものは、図１９の式４の形の文字列となる。

SLCPは辞書式順序にソート済みなので、A, C, G, T, Nからなる任意の文字列wに対して、先頭の文字から順にwをSLCPの要素と比較すれば、SLCP内にwで始まる要素が現れる直前の位置R(w)と直後の位置S(w)を定めることができる。ここで、SLCPの最初の要素の直前の位置は0で表し、SLCPの最後の要素の直後の位置は、SLCPの全要素数nで表す。nはリード配列データ内の総塩基数とリード配列の本数の和に等しく、図１９の式５で与えられる。

L S CP内にwで始まる要素が１つでも有る場合はR(w) ＜ S(w)となり、その差S(w) - R(w) はリード配列データ全体の中で文字列wが出現する回数、即ち、wの深度D(w) に等しい。一方、そのような要素が１つも無い場合は、R(w) = S(w) となり、これは、辞書式順序を崩さずにSLCPにwを追加する際の挿入位置を示す。（R(w), S(w)）を文字列wの順位区間とよぶ。なお、文字列の比較において、文字の比較は先頭文字から順に初めて異なる文字が現れるかまたは区切り文字$が現れるまで行われるため、wの末尾文字が$である場合にも、wに対する順位区間（R(w), S(w)）が定まる。

図６は、本実施例でBWTを利用する際に用いる補助関数を説明する説明図である。BWTはSLCP(501)の各要素（文字列）の末尾文字を連結したものであり、これらの文字列の長さ一般に異なるが、図６では、これらの末尾文字を右寄せして揃え、文字Aが所々に現れることを例示している。

0以上n以下の任意の整数rとA, C, G, T, N, $の何れかの文字zに対して、文字列BWTの先頭文字からr番目の文字までの範囲内に、文字zが出現する回数をO(z, r)で表す。文字列の先頭から指定した位置rまでの範囲における特定の文字zの出現回数O(z, r) は、ランク関数とよばれ、それを効率的に計算する方法は公知である（非特許文献７）。なお、r = nの場合は、各z = $, A, C, G, T, Nに対して、O(z, n) = n(z) となり、これらの値は、リード配列データ全体を事前に1回だけスキャンすることにより、既に求められている。

SLCPは辞書順にソートされているため、その要素を順に抜き出したものも辞書順にソートされている。特に、A, C, G, T, N, $の何れかの文字zに対して、zで始まる要素を抜き出したもの、即ち、R(z)からS(z)までの範囲は辞書順にソートされている。また同様に、A, C, G, T, N, $の何れかの文字zに対して、zで終わる要素を抜き出したもの全体も、辞書順にソートされている。SLCPは巡回置換（巡回シフト）で生成される全ての文字列から構成されているため、zで始まる要素の全体とzで終わる要素の全体は巡回置換により一対一に対応する。特に、A, C, G, T, Nからなる任意の文字列wに対して、zwで始まる文字列の全体は、wで始まりzで終わる文字列の全体と巡回置換により一対一に対応する。zwで始まる文字列の全体は、SLCP内の順位区間(R(zw), S(zw))で与えられるが、wで始まりzで終わる文字列全体は、SLCP内の順位区間(R(w), S(w))の中でzで終わる要素の全体によって与えられ、これらは一般にはSLCPの順位区間(R(w), S(w))内の飛び飛びの順位を占める。

図７は、本実施例において、この一対一の対応関係を利用して、検索キー配列に対する深度を計算する方法を示すフローチャート図である。キー配列kを入力し(701)、wをkの末尾の一文字zからなる文字列（kの接尾辞, suffix）とする。w=zに対するR(w)とS(w)の値は、全ての文字uに対してn(u)の値が既知であるため、直ちに計算できる(702)。接尾辞wがキー配列kに一致するときは、R(k)とS(k)の値から深度D(k)の値を計算して処理を終了し(706)、そうでない場合は、以下の処理を繰り返す(703)。キー配列kにおいて、接尾wの直前の文字をzとして、上述の一対一の対応関係を利用して、R(zw)とS(zw)の値を計算する(704)。ここで、R(z)の値は、702と同様に、zより順位の小さい文字uに対するn(u)の総和として計算できる。zwを新たなwとして(705)、703に戻り処理を繰り返す。

図８は、本実施例において、キー配列を含む全てのリード配列を、キー配列の左方にある塩基により分類して該当するリード配列の本数を求める方法を示したフローチャート図である。キー配列kに対するS(k)とR(k)の値は、上述の方法により計算済みであるとする。A, C, G, T, N, $の各文字zに対して、以下の処理を繰り返す(803)。上述の巡回置換による一対一の対応関係を利用すると、文字列zkに対する深度D(zk)を計算することができる(804)。この値が正のとき(805)、zkを含む配列がD(zk)本存在するので報告する(806)。

図８でキー配列を各塩基を用いて左方に1塩基だけ延長してそれぞれの深度を求める方法を示したで、この方法を繰り返し行うことにより、キー配列を種々の塩基を用いて左方に複数塩基延長し、それぞれの深度を計算することができる。また、それらの延長された配列をそれぞれの深度の値が示す重複度で重複させれば、リード配列データ全体の中から、元のキー配列を含むリード配列を左方に複数塩基延長したものを重複度込みで全て抽出できたことになる。

図９は、本実施例において、BWTの計算方法を示したフローチャート図である。BWTの計算においても、上述の巡回置換による一対一の対応関係を利用する。

先ず、901において、SLCP内の$で始まる要素n($)個からなるサブリストに対応する処理を行う。ここに対応するBWTの先頭のn($)文字の部分文字列をQ($)で表す。Q($)の中には、リード配列データの中の塩基文字A, C, G, T, Nが全て含まれる。リード配列データの中に空のリード配列がもし含まれていれば、それらは予め取り除いておく。その結果、Q($)の中に$は含まれていない。ソートの際の文字列比較では$より後ろの文字は比較されないため、Q($)の中の文字の並び順は任意でよい。そこで、Q($)を図１９の式６に示す文字列とする。文字の繰り返しは×で、文字の連結は＋またはΣで表現することにより、Q($)は901内に示した式により計算される。

また、全てのリード配列rの先頭に区切り文字$を付加し、任意の順番に並べたリストpを作る。pは文字列$の順位区間に属するSLCPの要素からなるリストであり、pの各要素の末尾文字はBWTの部分文字列としてQ($)の中に登録されている。また、z = $、I($)を空の数値リストとし、各y = A, C, G, T, Nに対して、P(y)とP’(y)を空リストに、Q’(y)を空文字列に、I’(y)を空の数値リストに初期化して、902以下の繰り返し処理に進む。( ) は空リストを、“ ”は空文字列を表す。

902では、pを或る文字列w$の順位区間に属するSLCPの要素の$以降の文字列からなるリストとし、pの各要素の末尾文字は既にBWTの部分文字列として何れかのQ(y)の中に登録されているが、巡回置換を施してそれらの末尾文字が先頭文字に来るように一対一に対応させたSLCPの要素（以下、シフトした要素とよぶ）の末尾文字は未だ何れのQ(y)にも登録されていないとする。

リストpの各要素を末尾文字で分類する。各y = A, C, G, T, Nに対して、yで終わる要素からなるpのサブリストをp(y)とする(902)。

次に、各y = A, C, G, T, Nに対して、p(y)の各要素の文字列から末尾文字yを削除する (903)。これらは、シフトした要素の$以降の文字列のリストである。

各y = A, C, G, T, Nに対して、p(y)の要素のうち$一文字だけからなるものの数をq($,y)とし、それら以外のp(y)の要素を末尾の文字x = A, C, G, T, Nで分類して数え、それらの数をq(x, y)とする(904)。

905において、各y = A, C, G, T, Nに対して、903で求めたリストp(y)の要素の末尾文字をBWTの部分文字列として文字列Q’(z)の後ろに登録する。ソートの際の文字列比較では$より後ろの文字は比較されないため、それらの中での文字の並び順は任意でよい。そこで、図１９の式７の形の文字列を連結して得られる文字列を新たなQ’(z)とする（905内の式では、文字の繰り返しは×で、文字の連結は＋またはΣで表す）。

これらの中で最初のq($, y)個の要素（末尾に$をもつ）に関しては、シフトした要素の末尾文字は、既に901でQ($)内に登録されている。一方、残りの下式で示す、
q(A, y) ＋ q(C, y) ＋ q(G, y) ＋ q(T, y) ＋ q(N, y) ‐‐‐（式８）
個の要素（末尾にA, C, G, T, N何れかの文字をもつ）に関しては、シフトした要素の末尾文字は未だ何れのQ(y)にも登録されていない。そこで、後続の処理のために、各y = A, C, G, T, Nに対して、リストP’(y)の後ろにリストp(y)を連結したものを新たなリストP’(y)とする（なお、図９中の905内の式で、リストの連結は＋で表す）。

また、数値リストI’(y)の後ろに２要素：
q($, y) , - ( q(A, y) ＋ q(C, y) ＋ q(G, y) ＋ q(T, y) ＋ q(N, y))‐‐‐（式９）
を追加する。ここで、正の数は、登録が済んだ要素の数、負の数は後続の処理で登録が必要な要素の数の符号を反転させたものを表す。

次に数値リストI(z)が空でないか調べる(906)。

それが空でない場合は、リストの先頭要素を取り出して削除し、取り出した要素をiとする。iが負の場合は、文字列リストP(z)から先頭の (- i ) 個の要素を取り出して削除し、取り出した要素からなるリストを新たなリストpとする。905で述べたことにより、このpは902で仮定した条件を満たしている。そこで、902に戻って処理を繰り返す。一方、iが正の場合、Q(z)の次のi個の要素をシフトした要素の末尾文字は、既に何れかのQ(y) 内に登録されているので、それらを新しいQ'(y)内に転記して、その転記量を各I'(y)に登録する(910)。

また、906で数値リストI(z)が空になった場合は、zが$かNに等しいかどうか調べる(911)。zがその何れにも等しくない場合は、zを次の順位の文字Next(z)に置き換え(912)、906に戻って処理を継続する。ここで、
Next(A) = C, Next(C) = G, Next(G) = T, Next(T) = N‐‐‐（式１０）
である。

911においてzが$またはNに等しい場合は、全てのz = A, C, G, T, Nに対して、P(z)が空リストかどうかを調べる(913)。それらの中に空リストでものがある場合は、全てのz = A, C, G, T, Nに対して、P’(z), Q’(z), I’(z)をそれぞれ新たなP(z), Q(z), I(z)として(914)、また、zを最初の順位の文字Aとして(915)、906に戻り処理を繰り返す。

913において、全てのz = A, C, G, T, Nに対してP(z)が空リストとなった場合は、文字列Q($), Q(A), Q(C), Q(G), Q(T), Q(N)を連結してBWTを得て、これを出力する(916)。

図１０は、図９中の910において旧文字列Qから新文字列Q'への転記方法を示す説明図である。各z = A, C, G, T, Nに対して、Q(z)内の転記元の位置を示す転記元ポインタと、Q (z)内の転記先の位置を示す転記先ポインタがある。これらは処理開始時に全て0に初期化されており、910内で繰り返し更新され、914内で全て0にリセットされる。また、転記対象を選ぶ選択先ポインタがあり、これはQ(A), Q(C), Q(G), Q(T), Q(N)内の何れかの位置を示す。選択先ポインタは、処理の開始時にQ(A)の先頭にリセットされており、910内で繰り返し更新され、912で新しいzに対するQ(z)の先頭にリセットされ、915でQ(A)の先頭にリセットされる。旧文字列Qに対しては、選択先ポインタと転記元ポインタの2種類が使われるため、図１０では、混乱を避けるために、選択先ポインタが指し示すQと転記元ポインタが指し示すQを繰り返し表示しているが、これらは同一物を指している。

910でQからQ'への1個の要素を転記するには、次のように行う。先ず、選択先ポインタが示す文字列Q(z)内の文字yを読み取り、選択先ポインタを＋１進める。図１０の例では、z = A, y = Cである。Q(y)を選択し、その転記元ポインタが示す文字（図の例ではT）を、Q'(y)の転記先ポインタが示す位置に転記し、転記元ポインタ、転記先ポインタとも＋１進める。910でQからQ'へのi個の要素を転記するには、1個の要素の転記をi回繰り返す。その結果、各 y = A, C, G, T, Nに対して、Q'(y)の転記先ポインタが＋i(y)進められたとすると、数値リストI'(y)の末尾に、i(y)を追加する。

図１１は、ゲノム配列DBとリード配列DBの構成を示す説明図である。リード配列DB(105)は、リード配列(106)のBWT (408)と、BWT上でランク関数の高速計算に必要な補助テーブル(1101)とから構成される。ゲノム配列DB(108)は、座標順に並んだゲノムの塩基配列データ(1102)とMLUデータ(1103)とから構成される。塩基配列データ(1102)からは、任意に指定された範囲の塩基配列を迅速に取り出すことができる。MLUデータは、二進符号化データ(1104)と検索用補助テーブル(1105)から構成される。

MLUの二進符号化データ(1104)は、長さ2nのバイナリ文字列であり、以下のように構成される。先ず、全ての要素を0に初期化する。参照ゲノム配列の全ての塩基位置xに対して、MLUの値L(x)を計算し、
k(x) = 2 x ＋ L(x) - 1 ‐‐‐（式１１）
として、二進符号化データのk(x) 番目の要素を1にセットする。但し、先頭要素は0番目の要素と数える。

末尾以外の任意のゲノム座標xで整数 l = L(x) - 1 をとると、xにおけるMLUの定義から、xからスタートする長さlの配列と等しい配列が、xとは異なる別の位置yにある。このとき先頭の一文字を無視すれば、x ＋ 1 からスタートする長さ( l - 1 ) の配列と等しい配列が、x ＋ 1とは異なる別の位置 y ＋ 1 にあることになる。従って、x ＋ 1におけるMLUの定義から、L(x ＋ 1) は少なくとも( l - 1 ) よりは大きくなければならない。

従って、
l - 1 ＜ L(x ＋ 1) ∴ L(x) - 1 ＝ l ≦ L(x ＋ 1) ‐‐‐（式１２）
となり、その結果
k(x) = 2 x＋L(x) -1≦2 x ＋L(x ＋1)＜ 2 x＋1＋L(x ＋ 1) = k(x ＋1)‐（式１３）
となる。

即ち、各塩基位置xに対してk(x) は相異なる値をとり、それらは二進符号化データの異なる要素を指示する。また、明らかに、x ＜ y ならば、k(x) ＜ k(y) である。従って、二進符号化データが得られていれば、任意のゲノム座標xにおけるL(x)の値を求めることができる。即ち、二進符号化データの中でx番目に1が現れる位置k(x)を求め、次のように計算すれば良い。

L(x) = k(x) - 2 x ＋ 1 ‐‐‐（式１４）
二進符号化データの中でx番目に1が現れる位置を求める関数k(x) = select(x) はセレクト関数と呼ばれており、補助テーブルを用いて効率的に計算する方法が知られている（非特許文献７）。1105は二進符号化データ(1104)上でセレクト関数を高速に計算する際に用いる補助テーブルである。

図１２は、本実施例の解析法における、ゲノム座標xにおけるMLUの値L(x)の計算方法を示すフローチャート図である。参照ゲノム配列データ(106)を入力し(1201)、両鎖を繋いだ参照ゲノム配列Gを作成し(1202)、GのサフィックスアレイSAを計算する(1203)。SAは、Gの全ての接尾辞(suffix)を辞書式順序にソートしたとき、接尾辞のスタート位置を示す整数をソート順に並べた整数配列である。ヒトゲノムの場合は、ゲノム・サイズが3ギガ塩基程度であるため、Gの長さは6ギガ程度となる。この程度の大きさのGに対しては、公知の方法(非特許文献８)を用いて、効率的にSAを計算することができる。SAは接尾辞のソート順位をスタート位置に変換する対応表であるが、その逆変換の対応表である逆サフィックスアレイ（ISA, inverse suffix array）を作成する(1204)。

また、Gの最長共通接頭辞長配列（LCP, longest common prefix length array）を計算する(1205)。Gの接尾辞を辞書式順序にソートしたときのr番目の要素をs(r)と表すと、LCPは整数配列であって、そのr番目の要素はs(r)と s(r - 1 )の最長の共通接頭辞の長さとして定義される。LCPは公知の方法（引用―LCP）を用いて効率的に計算することができる。s(r) と直前の要素s(r - 1)と先頭文字から順に比較すると第LCP(r) ＋ 1文字目の文字が異なる。同様に、s(r)を直後の要素s(r＋1)と先頭文字から順に比較すると第LCP(r＋ 1) ＋ 1文字目の文字が異なる。従って、下式１５で与えられる長さをもつs(r)の接頭辞はゲノム配列G内で一意的となる。

max ( LCP(r) ＋ 1, LCP(r＋1) ＋ 1 ) ‐‐‐（式１５）
ここで、maxは最大値をとることを意味する。従って、一意性を保証できる配列長（MLU, minimum length for uniqueness）は、1206内に示した式によって計算できる。

本実施例におけるDNAサンプルとしては、全ゲノム解析のサンプル、全エクソーム解析のサンプル、あるいは、注目するターゲット領域のDNA断片を濃縮したサンプル等を用いることができる。

実施例２として、注目する遺伝子領域の中から、スプライシングが生じている可能性が高い位置をMLUと深度に基づいて推定し、さらに、推定された位置において、配列比較に基づいてスプライシングの有無を判定する解析装置、及び方法の実施例を説明する。

図１３は、実施例２で、トランスクリプトームのスプライシング解析を行うための処理手順を示すフローチャート図である。

超並列型DNAシーケンサ(102)によりcDNAサンプル(1301)を解析して、多数の短い塩基配列からなるリード配列データ(103)を得る。以下、実施例１と同様に、リード配列DB化処理(104)を行い、リード配列DB(105)を得る。

参照ゲノム配列データ(106)に対しても、実施例１と同様に、ゲノム配列DB化処理(107)によりデータベース化して、ゲノム配列DB(108)を構築する。

注目する遺伝子領域の範囲を規定する座標(Xmin, Xmax)を入力する(109)。Xmin以上Xmax以下の各xに対して、実施例１と同様に、ゲノム配列DBに問合せてMLU、即ち、L(x)の値を得る(110)。また、リード配列DBに問合せて、各xにおける深度D(x)を得る(111)。

スプライシングの有無推定基準となるパラメータを入力し(1312)、xを注目する領域の左端Xminとして(113)、以下の繰り返し処理を開始する。MLUと深度に基づき、xにおけるスプライシングの有無を推定する(1314)。スプライシング無しと推定されるか、または、推定不能となった場合は、直ちにxの値をx＋1に更新する(120)。そうでない場合は、xの近くの別の塩基位置yを深度に基づいて一つ選択し、yからスタートして長さL(y)のゲノム部分配列をゲノム配列DBに問合せて取得する(115)。これをキー配列として、これを含む全てのリード配列を、リード配列DBに問合せて取得する(116)。これらのリード配列と参照ゲノム配列を比較することによりスプライシングの有無の判定を行う(1317)。スプライシング有りと判定された場合は、その結果を端末(122)や記憶装置(123)に出力する(119)。その後、xをx＋1に更新する(120)。xが領域の右端を超えていなければ、処理を繰り返す(121)。そうでなければ、処理を終了する。

図１４は、本実施例で順方向（ゲノム座標が増大する方向）に探索してスプライシング有りと推定された位置xにおいて、配列比較によりスプライシングの有無を判定する処理を説明するための説明図である。逆方向に探索する場合も同様である。横軸1401にゲノム位置座標、縦軸1402にカウント数をとり、グラフ1403は位置座標xにおける深度D(x)をプロットしたものである。ゲノム座標がzの位置で、その逆方向側（ゲノム座標が減少する方向）にzを末端とするイントロンがあり、その順方向側にエクソンがあるようなスプライシングが生じているとする。1404の範囲に示すように、位置xからスタートする長さL(x)のゲノム部分配列がイントロン内部に含まれるか、または、zを通るとき、深度D(x)の値はzの順方向側の周辺に比べて大きく低下する。即ち、zを境にしてzの逆方向側ではzの順方向側と比較して深度が大きく低下している。

そこで、順方向側と比較して逆方向側で深度が大きく低下する境界となるような位置座標xに対して、xの順方向側の近隣で深度が低下していない任意の位置yを一つとる。yからスタートして長さL(y)のゲノム部分配列をゲノム配列DBに問合せて取得し、その配列をキー配列としてリード配列DBに問合せて、キー配列を含むリード配列を全て集める。1405はx周辺のゲノム配列を表し、キー配列に該当する部分を下線で示した。1406は集められたリード配列を表し、キー配列に該当する部分を下線で示した。x周辺のゲノム配列1405は、ゲノム配列DBに問合せて得られる。

集められたリード配列1406は、大量にある全リード配列データのごく一部なので、以下の処理の計算コストは抑えられる。1411はリードの共通配列Sを表す。これは、1406に属するリード配列をキー配列の位置で合わせてアラインメントを行い、各塩基位置での最も高頻度で現れる塩基を並べて得られる配列である。共通配列S上では、キー配列（下線部分）で対応するゲノム配列Gの塩基位置座標を延長した座標系(1412)を用いる。共通配列Sを解析して、スプライシングが検出された場合、このスプライシングを報告する。

ここで、共通配列Sの部分文字列と、ゲノム配列データおよびリード配列データにおけるその深度を表すために、次の式１６に示す記法と用語を導入する。
S[z0, z1] : 位置座標が z0 以上 z1 以下の S の部分配列
Occ(s, G) : 塩基配列 s が G の部分配列として現れる回数（出現数）
（ゲノム配列データにおける s の深度）
Occ(s, R) : 塩基配列 s が R の部分配列として現れる回数（出現数）
（リード配列データにおける s の深度）
Loc(s, G) ： 塩基配列 s が G の部分配列として一意的に現れる
ときの出現位置座標 （Occ(s, G) = 1のとき）‐‐‐（式１６）

ここで、リード配列データにおけるsの深度は、実施例１で図7を用いて説明した方法により効率的に計算できる。また、ゲノム配列データにおけるsの深度も参照ゲノム配列GのBWTを用いて同様に計算できる。GのBWTはGのサフィックスアレイから直ちに計算することができ、実施例１で述べたように参照ゲノム配列Gのサフィックスアレイは公知の方法を用いて効率的に計算できる（非特許文献８）。また、塩基配列sがGの部分配列として現れる場所が一意のとき、その位置座標はGのBWTを用いて効率的に計算できることも公知である（非特許文献７）。

図１５は、本実施例で、リード配列の共通配列Sとx周辺のゲノム配列を比較して、xの位置でスプライシングの有無を判定する方法を示すフローチャート図である。ここで、p0, p1, p2, p3はユーザにより指示される判定基準パラメータである。先ず、変数x1をスプライシング有りと推定された位置座標xに初期設定する(1501)。x1とxを比較し(1502)、両者がp0以上離れている場合は、スプライシング無しと判定して処理を終了する(1517)。そうでない場合は、x0の値を更新して(1503)、Sの部分配列s を定める(1504)。ゲノム配列データにおけるsの深度 Occ(s, G) を計算して(1505)それが1より大の場合は、x0の更新処理1506に進む。それ以外の場合で、Occ(s, G)が１に等しい場合(1509)は、リード配列データにおけるsの深度Occ(s, R)がp1より大きい場合(1510)は、1511に進む。一方、Occ(s, G)が0の場合(1509)やOcc(s, R)がp1以下の場合(1510)はx1の更新処理1508に進み、x1とxの比較判定処理1502に戻って処理を継続する。また、x0の更新処理(1506)の後は、共通配列Sがx0まで定まっているか確認し、そうである場合は1504に戻ってsを更新し処理を継続する。

そうでない場合は、x1の更新処理1508に進んで処理を継続する。1511に進んだ場合はsの一意的な出現位置vを求め、xからvまでの距離がp2以下であれば、p2以下の長さの短い欠失を検出したと判定し(1513)、スプライシング無しと判定し(1517)、処理を終了する。そうでない場合は、vがxの左方にあり、かつ、xからvまでの距離がs3以下であれば(1514)、スプライシング有りと判定して処理を終了する(1515)。そうでない場合は、キメラ遺伝子（融合遺伝子）が検出されたと判定し(1516)、スプライシング有りと判定して(1515)、処理を終了する。

図１６は、本実施例において、順方向に探索して塩基位置xでスプライシングの有無を、MLUと深度に基づいて推定する方法を示すフローチャート図である。x自身を含め、xの周辺の各塩基位置yでは、ゲノム配列DBとリード配列DBへの問合せの結果として、MLUと深度の値、即ち、L(y)と D(y)の値を取得済みとする。また、d1, d2, h1, h2, h3, h4, l1は推定基準パラメータであり、これらの値は処理に先だって入力される。位置xにおいて深度D(x)が十分大きければスプライシング無しと推定する(1601)。xの周辺でMLUが大きくなることがあれば推定不能と推定する(1602)。xの順方向側の周辺で平均深度が十分大きくなければ、推定不能とする(1603)。Xの順方向側の周辺と比較して、xの逆方向側の周辺で深度が大きく低下していればスプライシング有りと推定し、そうでなければ、スプライシング無しと推定する(1604)。

以上説明した本発明では、各リード配列に対して、シーケンシング・エラーの種々の可能性を考慮してそのマッピング先を調べるような処理は一切行わない。リード配列データに対しては、（リード配列自身も含め）その全ての接尾辞を辞書式順序にソートするだけである。このようなソート処理は単純で、何らの任意性も残されていない。従って、処理法に依存するバイアスは発生せず、中立的な処理が可能である。また、処理の単純さゆえ、マッピング処理と比較して、計算コストも抑えることができる。

また、参照ゲノム側でのMLUの計算と、リード配列側での全ての接尾辞のソート処理は全く独立に行われる。そのため、複数の参照ゲノム配列が用意された場合、参照ゲノムとリード配列の組合せに依存する処理を行う必要はない。

さらに、下流の変異解析処理では、解析したい遺伝子領域内で、参照ゲノム配列データベースへの問合せ（MLU、または、MLUの長さの部分配列）とリード配列データベースへの問合せ（キー配列に対する深度、または、キー配列を含むリード配列）を行うことにより、変異が含まれる可能性が高い領域を推定し、全リード配列の中から対象とするリード配列を絞り込んで詳細な解析を行うことができる。従って、従来のマッピング処理を行わずに、本発明により対象を絞り込んで、変異解析を効率的に行うことが可能になる。

なお、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。例えば、上記した実施例は本発明のより良い理解のために詳細に説明したのであり、必ずしも説明の全ての構成を備えるものに限定されものではない。また、一の実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、一の実施例の構成に他の実施例の構成を加えることが可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。

更に、上述した各構成、機能、処理部等は、それらの一部又は全部を実現するプログラムを作成する例を説明したが、それらの一部又は全部を例えば集積回路で設計する等によりハードウェアで実現しても良いことは言うまでもない。

100、122 ユーザ端末
101 DNAサンプル
102 超並列型DNAシーケンサ
103 リード配列データ
104 リード配列データベース(DB)化処理
105 リード配列データベース（DB）
106 参照ゲノム配列データ
107 ゲノム配列データベース(DB)化処理
108 ゲノム配列データベース(DB)
123 ディスク
407 ソート済み巡回置換文字列リスト（SLCP）
408 リード配列データのBWT
501 ソート済み巡回置換文字列リスト（SLCP）
1104 ゲノム配列のMLU（minimum length for uniqueness）の二進符号化データ
1700 解析装置
1701 処理部（CPU）
1702 メモリ
1703 表示部
1704 記憶装置（HDD）
1705 入力部
1706 ネットワークインタフェース(NIF)
1707 バス

Claims (15)

1. 処理部と記憶部を備えたデータ解析装置であって、
   前記記憶部は、ゲノム配列データがデータベース化されたゲノム配列データベースと、リード配列データがデータベース化されたリード配列データベースを記憶し、
   前記処理部は、
   特定された解析対象のゲノム領域の配列に基づいて検索用の塩基配列であるキー配列を選択し、
   前記リード配列データベースの中にある前記キー配列の深度を求め、
   前記リード配列データベースの中にある前記キー配列を含むリード配列データを抽出し、抽出された前記リード配列データと前記ゲノム領域の配列とを比較してデータ解析を行う、
   ことを特徴とするデータ解析装置。
2. 請求項１に記載のデータ解析装置であって、
   前記ゲノム配列データベースは、問合せを受けた位置座標を始点とする部分配列が、ゲノム配列内で相補鎖も考慮して一意的となる最小の長さ（MLU, minimum length for uniqueness）を出力でき、
   前記処理部は、
   選択する前記キー配列として、前記解析対象のゲノム領域の部分配列であって、前記MLUの長さをもつものを選択する、
   ことを特徴とするデータ解析装置。
3. 請求項１に記載のデータ解析装置であって、
   前記処理部は、
   ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が局所的に低下する位置を変異がある可能性が高い位置として推定し、
   推定した前記位置で前記キー配列を含むリード配列データを前記リード配列データベースの中から抽出し、
   抽出された前記リード配列データと前記ゲノム領域の配列を比較して変異解析を行う、
   ことを特徴とするデータ解析装置。
4. 請求項１に記載のデータ解析装置であって、
   前記処理部は、
   ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が前方に比較して後方で局所的に低下する位置をスプライシングが生じている可能性が高い位置として推定し、
   前記スプライシングが生じている可能性が高いと推定された位置で、前記キー配列を含む前記リード配列データを前記リード配列データベースの中から抽出し、
   抽出された前記リード配列データとゲノム領域の配列を比較してスプライシング解析を行う、
   ことを特徴とするデータ解析装置。
5. 請求項１に記載のデータ解析装置であって、
   解析パラメータを入力する入力部を更に備え、
   前記処理部は、
   前記キー配列を含む前記リード配列データを前記リード配列データベースの中から抽出するか否かの判断を、前記リード配列データベースの中での前記キー配列の深度に基づいて行い、
   前記入力部から入力される解析パラメータにより、当該判断の基準を調整可能である、
   ことを特徴とするデータ解析装置。
6. 請求項５に記載のデータ解析装置であって、
   前記処理部は、
   ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が局所的に低下する位置を変異がある可能性が高い位置として推定し、
   推定した前記位置で前記キー配列を含むリード配列データを前記リード配列データベースの中から抽出し、
   抽出された前記リード配列データと前記ゲノム領域の配列を比較して変異解析を行い、
   前記解析パラメータにより、前記推定の基準を調整可能である、
   ことを特徴とするデータ解析装置。
7. 請求項５に記載のデータ解析装置であって、
   前記処理部は、
   ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が前方に比較して後方で局所的に低下する位置をスプライシングが生じている可能性が高い位置として推定し、
   前記スプライシングが生じている可能性が高いと推定された位置で、前記キー配列を含む前記リード配列データを前記リード配列データベースの中から抽出し、
   抽出された前記リード配列データとゲノム領域の配列を比較してスプライシング解析を行い、
   前記解析パラメータにより、前記推定の基準を調整可能である、
   ことを特徴とするデータ解析装置。
8. 請求項１に記載のデータ解析装置であって、
   前記処理部による前記データ解析の結果を表示する表示部を更に備える、
   ことを特徴とするデータ解析装置。
9. データ解析装置の処理部によるデータ解析方法であって、
   ゲノム配列データが検索可能な形式にデータベース化されたゲノム配列データベースと、リード配列データが検索可能な形式にデータベース化されたリード配列データベースを用い、
   特定された解析対象のゲノム領域の配列に基づいて検索用の塩基配列であるキー配列を選択し、
   前記リード配列データベースの中にある前記キー配列の深度を求め、
   前記リード配列データベースの中にある前記キー配列を含むリード配列データを抽出し、抽出された前記リード配列データと前記ゲノム領域の配列とを比較してデータ解析を行う、
   ことを特徴とするデータ解析方法。
10. 請求項９に記載のデータ解析方法であって、
    前記ゲノム配列データベースは、問合せを受けた位置座標を始点とする部分配列が、ゲノム配列内で相補鎖も考慮して一意的となる最小の長さ（MLU, minimum length for uniqueness）を出力可能であり、
    選択する前記キー配列として、前記解析対象のゲノム領域の部分配列であって、前記MLUの長さをもつものを選択する、
    ことを特徴とするデータ解析方法。
11. 請求項９に記載のデータ解析方法であって、
    ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が局所的に低下する位置を変異がある可能性が高い位置として推定し、
    推定した前記位置で前記キー配列を含むリード配列データを前記リード配列データベースの中から抽出し、
    抽出された前記リード配列データと前記ゲノム領域の配列を比較して変異解析を行う、
    ことを特徴とするデータ解析方法。
12. 請求項９に記載のデータ解析方法であって、
    ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が前方に比較して後方で局所的に低下する位置をスプライシングが生じている可能性が高い位置として推定し、
    前記スプライシングが生じている可能性が高いと推定された位置で、前記キー配列を含む前記リード配列データを前記リード配列データベースの中から抽出し、
    抽出された前記リード配列データとゲノム領域の配列を比較してスプライシング解析を行う、
    ことを特徴とするデータ解析方法。
13. 請求項９に記載のデータ解析方法であって、
    前記キー配列を含む前記リード配列データを前記リード配列データベースの中から抽出するか否かの判断を、前記リード配列データベースの中での前記キー配列の深度に基づいて行い、
    ユーザから指示される解析パラメータにより、当該判断の基準を調整する、
    ことを特徴とするデータ解析方法。
14. 請求項１３に記載のデータ解析方法であって、
    ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が局所的に低下する位置を変異がある可能性が高い位置として推定し、
    推定した前記位置で前記キー配列を含むリード配列データを前記リード配列データベースの中から抽出し、
    抽出された前記リード配列データと前記ゲノム領域の配列を比較して変異解析を行い、
    前記解析パラメータにより、前記推定の基準を調整する、
    ことを特徴とするデータ解析方法。
15. 請求項１３に記載のデータ解析方法であって、
    ゲノム領域内をスキャンして、前記リード配列データベースの中での前記キー配列の深度が前方に比較して後方で局所的に低下する位置をスプライシングが生じている可能性が高い位置として推定し、
    前記スプライシングが生じている可能性が高いと推定された位置で、前記キー配列を含む前記リード配列データを前記リード配列データベースの中から抽出し、
    抽出された前記リード配列データとゲノム領域の配列を比較してスプライシング解析を行い、
    前記解析パラメータにより、前記推定の基準を調整する、
    ことを特徴とするデータ解析方法。

Images