JP2009116559A - 大量配列の一括検索方法及び検索システム - Google Patents

Abstract

【課題】 大規模な文字配列データベースに対して、相互に類似した多数の文字配列に対する検索を一括して効率良く行う方法を提供する。特に、超並列DNAシーケンサを用いて生体内でのmRNAの発現解析を行う際、大量で冗長なmRNA断片配列データに対して、既知のゲノム配列（または既知の遺伝子配列）データベースに対する検索を、計算機を用いて一括して効率良く高速に行うための方法を提供する。
【解決手段】 ターゲット配列（ゲノム配列など）をサフィックス・アレイ等のようなデータ構造で表現するとともに、クエリー配列（mRNA断片配列など）もサフィックス・アレイ等のようなデータ構造で表現し、これらのデータ構造を比較することにより、サフィックスを共有する多数のクエリー配列に対する検索処理を一度に纏めて行い、多数のクエリー配列に対する検索を一括して効率良く行う。
【選択図】 図１

Description

本発明は、大規模な文字配列データベースに対して、相互に類似した多数の文字配列を一括して効率良く検索するための方法に関するものである。

生体内で発現している遺伝子を調べるために、生体サンプル内からmRNAを抽出してcDNAライブラリを作成し、PCR増幅してクローン選択を行い、キャピラリー型DNAシーケンサを用いてシーケンシングを行い、得られた配列と一致するようなゲノム中の遺伝子領域を検索する方法は、従来より一般的に広く行われている。通常、生体内では数万種類以上のmRNAが発現していると考えられているが、キャピラリー型DNAシーケンサで一度に解析できる配列本数は100本程度である。各mRNAの長さは数千塩基程度であることが多いが、キャピラリー型DNAシーケンサでは各配列を一度に500〜800塩基程度の長さまで読むことが出来る。

ヒトのゲノム配列は長大で30億塩基対にも及ぶため、効率的な検索方法が必要である。そこで、検索ターゲットとなるゲノム配列を予めインデクス付けしておくことにより、各mRNA断片配列と一致する配列をゲノム配列中から高速に検索することができる。また、全く無関係な２つの配列が十分な長さを持っていれば、それらが偶然一致するような確率は極めて低いため、検索処理の過程において検出された一致が生じる確率を評価して、統計的に有意でない（偶然の一致の可能性が疑われるような）情報をその過程で捨てることにより、高速化を図ることも行われている。

また、ゲノム配列中には、広範囲に非常に多数回繰り返し現れる反復配列が相当の割合を占めていることが知られている。このような反復配列を検索すると膨大なヒットが検出されることになり、検索処理の効率の低下を招く原因となる。そこで、一般的には、検索に先立って、反復配列が存在するゲノム領域を予めマスクしておくことが行われている。但し、反復配列に関する情報は完全には知られていないため、未知の反復配列による検索処理の効率低下、または、誤った反復配列情報に起因する検出漏れも生じうる。

しかしながら、現実的には多くの場合、前述の方法を用いることにより、具体的には、例えばBLAST（非特許文献１）とRepeatMasker（非特許文献２）のような専用のプログラムを用いることにより、キャピラリー型DNAシーケンサの一回のランで得られるmRNA断片配列データに対して、ゲノム配列データベースに対する検索を、検出漏れを抑えて高速に行うことができる。一般的に、検索に要する計算機処理時間はシーケンシングに要する時間よりも短いほどである。

一方、検出漏れを起こさないために統計的な評価に頼らず、また、不完全な反復配列情報にも頼らない、より精密な検索方法も知られている。その代表的な方法は、ゲノム配列をサフィックス・ツリーやサフィックス・アレイ等のようなフル・テキスト・インデックス化したデータ構造で表現し、そのデータ構造を利用して高速に配列検索を行う方法である（非特許文献３）。そのようなデータ構造は非常に大規模になるため、効率的に構築することが容易ではなく、そのメモリ使用量を削減するための方法や、その構築処理を高速化する方法が、多々提案されてきている。しかしながら、このようなデータ構造は一度構築してしまえば、一回一回の検索は非常に高速に行うことができる。

さらに、比較ゲノムの分野では、インデクス化されたデータ構造どうしを比較する方法が用いられている（非特許文献４）。そこでは、生物のゲノム配列の進化を解明することを目的として、2種類の生物種間で対応するゲノム配列の部分領域であって、進化的に保存されたものを求める問題が検討されている。そのような保存領域は、進化の過程で受けた変異の蓄積により、完全に一致する訳ではないため、直接探し出すことは困難である。そこで、保存領域を探す手がかりとして、短い配列であって、2種類の生物種のゲノム配列中に共有され、かつ、それぞれの生物種のゲノム配列中には1度しか現れないような配列を求める問題が検討されている。このような短い配列は、2種類の生物種のゲノム配列間の対応関係の候補をピンポイントで示していると考えられる。この問題の解決のために、それぞれの生物種のゲノム配列（数本から数十本程度の巨大な染色体配列からなる）をサフィックス・アレイ等のようなデータ構造で表現し、それらのデータ構造を比較することにより、その回答を得ようとする方法が行われている。

一方、近年、キャピラリー型DNAシーケンサとは全く異なる、新しい原理に基づいた超並列DNAシーケンサが出現した（非特許文献５）。超並列DNAシーケンサは、一度に数千万本にも及ぶ配列を読み取ることができる。このような圧倒的な並列性を活かすことにより、一細胞内で数万種類以上も発現しているmRNAを網羅的に解析することが可能になる。

しかしながら、超並列DNAシーケンサで得られる配列データでは、配列1本当たりの塩基長が30塩基長程度と短い。そのため、類似した遺伝子やスプライス・バリアントや選択的プロモータの違いを区別することが困難な場合も生ずる。それでもなお、大多数の場合は、発現している遺伝子を同定することが可能である。

このような超並列DNAシーケンサの出現により、配列検索処理に求められる性能も変化してきている。一回のランで一度に数千万本の配列が得られるため、一度に数千万本の配列に対して既知配列データベース（ゲノム配列または遺伝子配列）に対する検索を行う必要がでてきた。そのために要する計算機処理コストも膨大となり、大型の並列計算機が必要となってきた。

また、超並列DNAシーケンサから得られる配列の配列長が短いため、配列一致の統計的有意性の評価を利用する従来の高速配列検索方法では、単なる偶然の一致として評価されるような場合が生じ易く、それが原因となって検出漏れが多々生じる可能性が高い。そこで、サフィックス・アレイなどを用いた、より精密な検索方法を用いる必要が出てくる。このような精密な検索方法では、一旦インデクス化されたデータ構造を構築してしまえば高速な検索が可能ではあるが、一度に数千万本もの配列を検索するとなると、やはり大型の並列計算機が必要となってくる。

多数回の検索は、それぞれ独立に一回ごとの検索の繰り返しとして行われる。並列計算機ではその独立性を利用して、各計算ノードは相互に通信し合うことなく、それぞれ独立に検索処理を平行して行っている。

Altschul， S.F.， Madden， T.L.， Schaeffer， A.A.， Zhang， J.， Zhang， Z.， Miller， W. & Lipman, D.J. (1997) "Gapped BLAST and PSI-BLAST: a new generation of protein database search programs." Nucleic Acids Res. 25:3389-3402. http://www.repeatmasker.org/webrepeatmaskerhelp.html Sadakane K， Shibuya T. "Indexing huge genome sequences for solving various problems." Genome Inform. 2001;12:175-83. Ross Lippert, Space-efficient whole genome comparisons with Burrows-Wheeler transforms, Journal of Computational Biology, 12(4), pp. 407-415, 2005. Service RF. "Gene sequencing. The race for the $1000 genome." Science. 2006 Mar 17;311(5767):1544-6.

超並列DNAシーケンサの一回のランで得られるような、数千万本以上にも及ぶ多数の配列に対して、一度に配列データベース検索を行う際、それらの多数の検索配列の中に同じものまたは類似したものが多数含まれていると、それぞれの検索を独立に行う従来の方法では、同じような検索処理を何度も繰り返すという無駄が生じる。特に発現解析データにおいては、発現量の大きな遺伝子に由来するmRNA配列が検索配列データの中に何度も繰り返し冗長に現れるため、そのような状況が生じ易い。

そこで、本発明の目的は、多数の検索配列の相互の類似性を組織的に調べて利用し、無駄な処理の回避を図り、より効率的に多数の検索を一括して行う方法を提供することにある。

本発明では、多数回の検索を一括して効率的に行う方法として、ターゲット側の既知配列をサフィックス・アレイ等のようなデータ構造で表現するとともに、数千万本にも及ぶ多数の検索配列データもサフィックス・アレイ等のようなデータ構造で表現し、それらのデータ構造の比較により、多数の検索配列に対する処理を一度に纏めて効率良く行う。

多数の類似した検索配列が、サフィックス・アレイ等のようなデータ構造上で一箇所に纏まって表現されるため、それらに対する検索処理を一度に纏めて行うことが出来る。それにより、特に、発現量の大きな遺伝子に由来するmRNA配列に対して同じ検索が何度も繰り返される無駄が生じる問題が解決される。

（実施例１）
以下、本発明の実施例を図面を用いて詳細に説明する。

本実施例では、ターゲット配列は既知のゲノム配列（数十本程度の染色体配列からなる）であり、塩基種別を表す４種類の文字A, C, G, Tからなる文字列であるとする。また、クエリー配列は大量のmRNA断片配列であり、同じく４種類の文字A, C, G, Tからなる文字列であるとする。本実施例では、これらのmRNA断片配列全てに対して、一括してゲノム配列への検索を行い、ゲノム配列中に完全一致する部分配列として出現する位置（マッピング位置）を報告するものとする。但し、マッピング位置の数が、予め定められた規定値を超えるときは、そのようなmRNA断片配列はリピート配列であると報告し、マッピング位置の報告は行わない。また、マッピング位置を持たないようなmRNA断片配列については、マッピング可能な（ゲノム配列中に完全一致する部分配列として出現する）最長のサフィックスを報告するものとする。

図1は、本発明の特徴を最も良く現している大量配列の一括検索方法の説明図である。この図では、計算機主記憶上に構築されるデータ構造とそれら相互の論理的関係、並びに、説明用の補助情報を表現している。

現在検索処理中の塩基配列データ（文字列）１は、多数の配列の一括検索において、一般には複数のクエリー配列の末尾に共通して現れる配列である。クエリー側サフィックス・アレイ２は、全てのクエリー配列を区切り文字を介して連結して得られる一本の巨大配列に対して、その全てのサフィックスを求め、それらを辞書式順番にソートしたときの、サフィックス位置を表す整数を並べた配列データである。ここで、文字列のサフィックスとは、文字列中の任意の位置から文字列末尾までからなる部分文字列を意味する。図１では、クエリー側サフィックス・アレイ２の中央に、実際に計算機主記憶上に記憶される整数配列データ（サフィックス位置を格納している）を示し、その左側にはアレイ内の順位を示し、またその右側には、説明用の補助情報として、対応するサフィックス内の塩基の並び（文字列）を表示している。実際には、このような補助情報やアレイ内の順位は計算機主記憶上に記憶する必要はない。クエリー側サフィックス・アレイ２の中で、２つの整数i0とi1は、現在検索中の配列１に対応するサフィックスが、ソートされたインデクス上でi0番目からi1番目までの順位の範囲を占めており、サフィックス位置がp0，…,p1で与えられことを表現している。

同様に、ターゲット側サフィックス・アレイ３は、全てのターゲット配列を区切り文字を介して連結して得られる一本の巨大配列に対して、その全てのサフィックスを求め、それらを辞書式順番にソートしたときの、サフィックス位置を表す整数を並べた配列データである。図１では、ターゲット側サフィックス・アレイ３の中央に、実際に計算機主記憶上に記憶される整数配列データ（サフィックス位置を格納している）を示し、その左側にはアレイ内の順位を示し、またその右側には、説明用の補助情報として、対応するサフィックス配列の塩基の並びを表示している。実際には、このような補助情報やアレイ内の順位は計算機主記憶上に記憶する必要はない。ターゲット側サフィックス・アレイ３の中で、２つの整数j0とj1は、現在検索中の配列１に対応するサフィックスが、ソートされたインデクス上でj0番目からj1番目までの順位の範囲を占めており、サフィックス位置がq0,…, q1で与えられることを表現している。
データの参照関係４は、現在検索処理中の作業配列１に対して、クエリー側サフィックス・アレイ２内の範囲(i0,i1)とターゲット側サフィックス・アレイ３内の範囲(j0,j1)が対応していることを表している。

新しい検索配列５及び６は、現在検索中の配列を先頭から1文字だけ伸長させたものである。これらは、一般には文字の種類だけ複数の候補がある。クエリー側サフィックス・アレイ７は、計算機主記憶上ではクエリー側サフィックス２と同一のものであるが、新しい検索配列からその内部への参照先が変化している。新しい検索配列５からの参照先は(i00, i01)であり、また、新しい検索配列６からの参照先は(i10, i11)である。参照先の変化を表現する状態遷移５は、クエリー側サフィックス・アレイ内部への参照先が変化したことを表現している。

同様に、ターゲット側サフィックス・アレイ９は、計算機主記憶上ではターゲット側サフィックス３と同一のものであるが、新しい検索配列からその内部への参照先が変化している。新しい検索配列５からの参照先は(j00, j01)であり、また、新しい検索配列６からの参照先は(j10, j11)である。参照先の変化を表現する状態遷移10は、ターゲット側サフィックス・アレイ内部への参照先が変化したことを表現している。

データの参照関係11は、新しい作業配列５に対して、クエリー側サフィックス・アレイ７内の範囲(i00,i01)とターゲット側サフィックス・アレイ９内の範囲(j00,j01)が対応していることを表している。同様に、データの参照関係12は、新しい作業配列６に対して、クエリー側サフィックス・アレイ７内の範囲(i10,i11)とターゲット側サフィックス・アレイ９内の範囲(j10,j11)が対応していることを表している。

破線１３は、計算機主記憶上において、クエリー側の情報とターゲット側の情報を論理的に区分している説明用の分割線である。

本発明は、以上の構成を備えているので、多数の配列に対する検索を、それらの類似性を利用して、一括して効率的に行うことが出来る。図2は、本発明の動作を説明するための、計算機処理手順の概要を示すフローチャートである。

21は多数のクエリー配列からなる入力データを表す。クエリー配列は有限種類の文字からなる文字列である。本実施例では、クエリー配列はmRNA断片配列であり、塩基種別を表す４種類の文字A, C, G, Tからなる文字列である。

22では、多数のクエリー配列を区切り文字$を挟んで連結して一本の巨大な文字列とし、そのサフィックス・アレイを構築する。これは、その巨大な文字列の全てのサフィックス（途中位置から末尾までの部分文字列）を求め、それらを辞書式順番にソートしたときの、サフィックス位置を表す整数を並べた配列データである。サフィックス・アレイは、非特許文献６などの公知の技術を用いることにより、計算機上で効率的に構築することが出来る。これをクエリー側サフィックス・アレイとよぶ。

23はもう一つの入力データであるターゲット配列データを表す。ターゲット配列も、クエリー配列と同じ文字種からなる文字列である。本実施例では、ターゲット配列は既知のゲノム配列（染色体配列）であり、同じく４種類の文字A, C, G, Tからなる文字列であるとする。

上記と同様に、24では、複数のターゲット配列を区切り文字$を挟んで連結して一本の巨大な文字列とし、非特許文献６などの公知の技術を用いて、そのサフィックス・アレイを構築する。これをターゲット側サフィックス・アレイとよぶ。

サフィックス・アレイ上では、全てのサフィックスが辞書式順番にソートされているため、一つの文字列Sが与えられたとき、Sで始まるサフィックスが複数ある場合は、それらはサフィックス・アレイ上の連続した位置を占める。従って、その範囲はサフィックス・アレイ上の開始位置と終了位置のソート順位を表す整数のペア(i0,, i1)で表すことができる。これを、Sのサフィックス・アレイ内出現位置とよぶ。

検索処理の過程では、処理中の検索文字列Sの長さを1文字ずつ伸長しながら、その出現位置を絞り込んでゆく。そこで、処理中の検索文字列Sに対しては、そのクエリー側サフィックス・アレイ内出現位置(i0, i1)とそのターゲット側サフィックス・アレイ内出現位置(j0, j1)を常に関連付けながら処理を進める。

25では、そのための初期化を行う。すなわち、文字を全く含まない空文字列を初期の検索文字列Sとして、空文字列のクエリー側サフィックス・アレイ内出現位置をその全範囲(0, Nq)とし、また、空文字列のターゲット側サフィックス・アレイ内出現位置をその全範囲(0, Nt)とする。ここで、Nq及びNtは、それぞれ、クエリー側、ターゲット側の全配列を連結して出来た一本の巨大な文字列の長さである。

26では、以下に図3を用いて説明する「検索文字列の1文字伸長処理および棄却判断」を行う。1文字伸長の仕方は、一般に文字種に対応するだけの可能性があるため、それら複数の可能性を平行して調べる。図2に示されているように、26はループ処理の中に含まれているため、ループ処理を繰り返すたびに、平行して調べる場合の数は増えていく。しかし、以下に図3を用いて説明するように、「検索文字列の1文字伸長処理および棄却判断」26では、不要になった検索文字列が棄却されるため、実際には、平行して調べる場合の数は抑制される。特に、ループの回数は、クエリー配列データに含まれる最長のものの長さを超えることはない。そこで、逐次計算機上で複数の可能性を平行して調べるために、深さ優先探索を行う。

全長一致するようなターゲット配列の部分配列を持たないクエリー配列に対しては、図3を用いて以下に説明するように、26の内部で検索文字列の棄却の判断が行われ、そのようなクエリー配列のターゲット配列に対する最長一致サフィックスを報告する(27)。さらに、後述の31と同様な処理により、それらの最長一致サフィックスがターゲット配列中に出現する位置（マッピング位置）を報告することもできる。

28では、平行して調べている多数の検索文字列の中から、それと全長一致するクエリー配列をもつような検索文字列を抽出して、それらに対して29の繰り返し処理を行う。なお、平行して調べている全ての検索文字列に対して、それが28で抽出対象になったか否かに拘わらず、26に戻るループ処理は継続する。

28において、検索文字列Sに全長一致するクエリー配列があるか否かを調べるためには、区切り文字$を使って検索文字列を1文字伸長させて$Sとし、そのクエリー側サフィックス・アレイ内出現位置を調べればよい。そのような出現位置(i0, i1)が見つかれば、それに対応するi1-i0+1本のクエリー配列がSと全長一致する。また、そのような出現位置が見つからなければ、Sと全長一致するようなクエリー配列は存在しない。区切り文字$による検索文字列の1文字伸長処理と、それに伴う、クエリー側サフィックス・アレイ内出現位置の更新方法は、以下で図3を用いて説明する「検索文字列の1文字伸長処理および棄却判断」26と同様な方法で実現できる。

繰り返し処理29では、28で抽出した検索文字列Sごとに、30以下の判断と処理を行う。30の判断で必要となるSのターゲット配列内出現回数は、Sのターゲット側サフィックス・アレイ内出現位置(j0, j1)を用いてj1-j0+1と計算される。これが規定値以下に収まっていれば、31においてSと全長一致するクエリー配列に対して、それらのターゲット配列内の部分配列としての出現位置（マッピング位置）を報告する。そうでない場合は、32において、Sと全長一致するクエリー配列はリピート配列であると報告する。

ここで、Sと全長一致するクエリー配列を求めるためには、以下のようにすれば良い。Sのクエリー側サフィックス・アレイ内出現位置(i0, i1)の範囲に含まれるサフィックス位置のリストp0，…,p1を求める。これらのサフィックス位置はそれぞれ一本のクエリー配列に対応している。サフィックス位置pでの区切り文字$に関するrank関数の値ｎ=rank(p)を計算すれば、対応するクエリー配列がn+1番目のものであることが分かる。ここで、pにおけるrank関数の値がnであるとは、pより前方にn個の区切り文字$があることを意味する。このようなrank関数は、非特許文献７などの公知の方法を用いて効率的に計算できる。 Sと全長一致するターゲット配列の出現位置を求めるには、以下のようにすれば良い。本実施例では、ターゲット配列23は既知ゲノム配列であり、24では各染色体配列を染色体番号順に区切り文字$を介して連結して連結して一本の巨大配列を作り、ターゲット側サフィックス・アレイを構築したとする。先ず、Sのターゲット側サフィックス・アレイ内出現位置(j0, j1)の範囲に含まれるサフィックス位置のリストq0,…, q1を求める。これらのサフィックス位置はそれぞれ一本の染色体配列内の特定の位置に対応している。サフィックス位置qでの区切り文字$に関するrank関数の値ｎ=rank(q)を計算すれば、n+1番目の染色体が対応していることが分かる。次に区切り文字$に関するselect関数のnにおける値s=select(n)を計算すれば、サフィックス位置qがn+1番目の染色体のq−s塩基目の位置であることが分かる。ここで、区切り文字$に関するselect関数のnにおける値がsであるとは、区切り文字$がn番目に現れる位置が、染色体配列を連結した巨大配列内のs番目の文字位置であることを意味する。このようなselect関数は、非特許文献７などの公知の方法を用いて効率的に計算できる。

このようにして、1本の検索文字列Sを通じて、Sのクエリー側サフィックス・アレイ内出現位置(i0, i1)の範囲に対応するi1−i0+1本のクエリー配列に対して、一括して、それらのゲノム配列内の出現位置を求めることができる。クエリー配列データの冗長性が高いときは、クエリー配列の本数i1−i0+1は大きい場合もあり、クエリー配列ごとに独立に検索を行う従来法ではi1−i0+1回繰り返される検索処理が、本発明では1回に纏めて行われるため、処理の効率化が達成される。

次に図3と図1を用いて、「検索文字列の1文字伸長処理と棄却判断」26の内部処理を説明する。図3において41は図2の処理26の入力データであり、検索文字列Sとそのクエリー側サフィックス・アレイ内出現位置(i0, i1)とそのターゲット側サフィックス・アレイ内出現位置(j0, j1)からなる。繰り返し処理42では、塩基を表す文字A, G, C, Tごとに以下の処理を繰り返す。今、選んだ文字をxで表す（xはA, G, C, Tの何れかである）。43では、xSを現在処理中の検索文字列であるとする。44では、xSのクエリー側サフィックス・アレイ内の出現位置を求める。これは、Sのクエリー側サフィックス・アレイ内の出現位置(i0, i1)に基づいて、非特許文献４などの公知の技術を用いることにより、効率的に計算することができる。

45では、xSのクエリー側サフィックス・アレイ内の出現位置が空か否かを判断している。これが空になった場合は、xSはクエリー配列に含まれていないため、xSの検索を継続する必要性が無いので、46において検索配列xSを棄却する。

そうでない場合は、47において、xSのターゲット側サフィックス・アレイ内出現位置を求める。これも、Sのターゲット側サフィックス・アレイ内出現位置(j0, j1)に基づいて、非特許文献４などの公知の技術を用いることにより、効率的に計算することができる。

48では、xSのターゲット側サフィックス・アレイ内の出現位置が空か否かを判断している。これが空になった場合は、xSはターゲット配列に含まれず、xSをサフィックスにもつようなクエリー配列のターゲット配列に対する最長の部分一致サフィックスはSになることが分かる。特に、これらのクエリー配列が全長一致するような、ターゲット内部分配列は存在しない。そこで、これらのクエリー配列に対しては、49において最長の部分一致サフィックスがSであることを報告して、50において検索配列xSを棄却する。xSをサフィックスにもつようなクエリー配列を求めるためには、xSのクエリー側サフィックス・アレイ内の出現位置の範囲に含まれるサフィックス位置のリストを作り、前述した図2の28の内部処理と同様な方法により、区切り文字$に関するそれらのサフィックス位置でのrank関数の値を用いて計算すればよい。

48では、xSのターゲット側サフィックス・アレイ内の出現位置が空でないと判断された場合は、51のように新たな検索文字列xSとそのクエリー側サフィックス・アレイ内出現位置とそのターゲット側サフィックス・アレイ内出現位置を、「検索文字列の1文字伸長処理および棄却判断」26の出力とする。

図1の例では1に示すようにS=TCGであり、状態遷移8のそれぞれの矢印は、繰り返し処理42においてx=A, C, Gを選んだ場合に対応し、ATCG, CTCG, GTCGのクエリー側サフィックス・アレイ内の出現位置がそれぞれ、(i00, i01)、(i10, i11)、(i20, i21)であることを示している。また、x=Tを選んだ場合、TTCGがクエリー側サフィックス・アレイに出現しないとしため、新たな検索文字列TTCGは棄却される。

状態遷移10のそれぞれの矢印は、繰り返し処理42においてx=A, Cを選んだ場合に対応し、ATCG, CTCGのターゲット側サフィックス・アレイ内の出現位置がそれぞれ、(j00, j01)、(j10, j11)であることを示している。ここで、x=Tを選んだ場合の検索文字列TTCGは、46において棄却されているため、図に示す(j20, j21)が検索されることはない。また、x=Gを選んだ場合、GTCGがターゲット側サフィックス・アレイには出現しないとしため、GTCGをサフィックスにもつクエリー配列のターゲット配列に対する最長の部分一致サフィックスはTCGとなる。そこで、49においてこれらのクエリー配列に対する最長の部分一致サフィックスがTCGであることを報告して、新たな検索文字列GTCGは棄却される。

最終的に、ATCGとCTCGが1塩基伸長された新たな検索文字列であり、それらのクエリー側サフィックス・アレイ内出現位置とそのターゲット側サフィックス・アレイ内出現位置がそれぞれ、(i00, i01)と(j00, j01)、及び、(i10, i11)と(j10, j11)となる。

次に図6を用いて、本発明によりmRNA発現解析の計算を行うシステム全体の構成を説明する。91は、本発明による検索処理の計算を行う計算機である。92は超並列DNAシーケンサ装置でありcDNAサンプル90を解析して、cDNA配列データ93を生成する。また、94は外部記憶装置であり、ゲノム配列データ95と遺伝子領域座標データ96を保持している。計算機91は、超並列DNAシーケンサ装置92および外部記憶装置94とネットワーク97を介して接続されている。98は計算機91に繋がれた表示出力装置であり、計算機91で計算された発現解析結果99を表示する。

計算機91の内部で行われる処理の概要を図7に示す。入力処理部101では、cDNA配列データ93、ゲノム配列データ95、および、遺伝子領域座標データ96を入力する。マッピング位置計算処理部102では、cDNA配列データ93をクエリー配列データ21として、また、ゲノム配列データ95をターゲット配列データ23として、図2を用いて説明した検索処理を行い、各クエリー配列に対するマッピング位置31を計算して、それらの計算結果をcDNAマッピング位置データ103として計算機91内部のメイン・メモリ104に保持する。座標データ比較処理部105では、これらのcDNAマッピング位置データ103と各遺伝子の遺伝子領域座標データ96を比較し、遺伝子領域と同じ染色体上でその開始位置から終了位置までの範囲内にマッピング位置が含まれるようなクエリー配列を、その遺伝子に属するクエリー配列であると判定する。配列カウント処理部104では、各遺伝子に属するクエリー配列の数を数える。これが、細胞内での各遺伝子の発現の強さを表す指標（発現量）となる。出力処理部105では、そのカウント結果を発現解析結果99として出力する。
（実施例２）
前記実施例においては、mRNA断片配列のゲノム配列へのマッピング位置を求める条件として、mRNA断片配列がゲノム配列中の完全一致する部分配列として出現することとした。しかし、実際には、SNP（1塩基多型）などが知られているように、レファランスとなるゲノム配列とシーケンシングで得たmRNA断片配列が正しく対応している場合であっても、ゲノムやmRNAを採取した個人の違いにより、それらの塩基の並びに若干の相違が見られることがある。そこで、このような場合にも、正しい対応関係を求めるために、mRNA断片配列のゲノム配列へのマッピング位置を求める条件を緩和して、mRNA断片配列がゲノム配列中に近似的に一致する部分配列として出現することとする。ここで２つの配列が近似的に一致するとは、予め指定した許容数以下の箇所の文字の置換や挿入や欠失により、それら２つの配列が一致することである。このときの、文字の置換や挿入や欠失の数を不一致数とよぶ。

また、ゲノム配列やmRNA断片配列には、正確に塩基の種類を決めることができなかった場合の塩基位置を表す文字Nが現れることもある。このような場合は、Nは他の任意の塩基A, G, C, T, Nと対応するものとして扱う。但し、Nを多数含むような配列に含まれる情報は少ないので、mRNA断片配列のゲノム配列へのマッピング位置を求める際には、そこに現れるNの数は予め定めた規定値以下となるように制限する。

このようにマッピングの条件を緩和した場合でも、前記実施例と同様な方法により、大量のmRNA断片配列に対するマッピング位置を一括して効率良く求めることが出来る。そのためには、検索処理中の作業用の文字列Sに対して、ターゲット側サフィックス・アレイにおいて、Sと近似的に一致するサフィックスが出現する範囲を調べる必要があり、また、それらは一般に複数種類あるためリスト化して管理する必要がある。また、それらの不一致数やNの数が規定値以下に納まっていることを保障する必要がある。（なお、クエリー側サフィックス・アレイに対しては、検索文字列Sと完全一致するサフィックスが出現する範囲を調べればよいので、前記実施例からの変更の必要はない。）
以下、図4を用いてこのような、近似的に一致するターゲット配列内のサフィックスの管理方法を説明する。61は検索処理中の作業用の文字列Sであり、この例では、ACGである。62は、Sと近似的に一致するターゲット側のサフィックスを管理するリストである。各行がそれぞれSと近似的に一致するターゲット側のサフィックスQに対応する。この例ではQは3種類ある。62の第1列目と第2列目は、それらのサフィックスのターゲット側サフィックス・アレイ63内の出現範囲の開始と終了位置であり、ターゲット側サフィックス・アレイ63内でのソート順位を表す整数である。ターゲット側サフィックス・アレイ63内は、前記実施例の図１における7や9と同一のものであるが、その内部で参照される範囲が複数ある点が異なっている。また、62の第3列目は、SとQの不一致数（それらの間で、塩基の置換、挿入、欠失により差異が生じている箇所の数）とNの数の和である。不一致数とNの数は独立にカウントし、サフィックスを棄却するか否かの判定に用いることも出来るが、本実施例では、不一致数とNの数との和を用いることとし、これを簡単のため単に不一致数とよぶことにする。

処理の流れは、前記実施例の図2および図3と同様である。但し、図3の47において、前記実施例においてはSと完全一致するサフィックスがターゲット側サフィックス・アレイ内に出現範囲する範囲を一つ求めれば良かったのに対し、本実施例では、Sと近似的に一致する複数種のサフィックスQに対して、Qがターゲット側サフィックス・アレイ内の出現範囲をそれぞれ求める必要がある。そこで、図3の47に代わる処理の流れを、図５を用いて説明する。

今、作業用の検索文字列Sに対して、その1文字伸長処理を行って、xSを新たな検索文字列にしようとしているとする。ここで、xはA, C, G, T, Nの何れかの文字を表す。図4のようにリスト管理されているSと近似的に一致しているそれぞれのサフィックスQに対して、図5に示すような処理を行う。

71はこれから処理しようとするサフィックスQである。72においてyをA, C, G, T, Nの何れかの文字として、文字種ごとに以下の処理を繰り返す。Qの先頭に1文字yを追加して新たなサフィックスyQを作り(73)、そのターゲット側サフィックス・アレイ内の出現範囲を求める(74)。これは、前記実施例の47と同様にできる。その出現範囲が空でないか調べ(75)、もし空ならばyQを棄却する(76)。そうでない場合は、Sに対する伸長文字xとQに対する伸長文字yが不一致であるか、または、xかyの少なくとも一方がNであるかをしらべ(77)、もしそうであるならば、yQの不一致数を１増加する(78)。これは、サフィックスの1文字伸長処理に伴う、塩基の置換による不一致数の増加、または、Nの数の増加を調べている。こうして得られたサフィックスyQは86の判断に委ね、再び72に戻って次の文字種に対して上記の処理を繰り返す。

また、79においては、Qの不一致数を１だけ増加させ、直ちに86の判断に委ねる。これは、ターゲット配列側において欠失が新たに1箇所増える場合を評価している。

80においては、1塩基挿入が起きる場合の評価を目的として、yQに対してさらに1文字の伸長処理を行う処理を行うために、zをA, C, G, T, Nの何れかの文字として、文字種ごとに以下の処理を繰り返す。yQの先頭に1文字zを追加して新たなサフィックスzyQを作り(81)、そのターゲット側サフィックス・アレイ内の出現範囲を求める(82)。これは、前記実施例の47と同様にできる。その出現範囲が空でないか調べ(83)、もし空ならばzyQを棄却する(84)。そうでない場合は、zyQの不一致数を１増加する(78)。これは、ターゲット側において１塩基挿入が起きる場合の不一致数の増加である。こうして得られたサフィックスzyQは86の判断に委ね、再び80に戻って次の文字種に対して上記の処理を繰り返す。

このようにして新たに得られた全てのサフィックスに対して、その不一致数が規定値以下か検査し(86)、そうでない場合は、そのサフィックスを棄却する(87)。そうである場合は、xSと近似的に一致するサフィックス(88)が新たに求められたことになるので、これを図４に示したようなリストに追加し、古いサフィックスQに関する情報はそのリストから削除する。

以上のようにして、検索文字列の1文字伸長処理において、図4に示した検索文字列と近似的に一致するターゲット配列内サフィックスの出現範囲のリストを正しく管理することができる。
（実施例３）
前記実施例においては、ターゲット側とクエリー側のサフィックスの出現範囲を管理するために、図１の2や3のようなサフィックス・アレイを用いている。例えば、ヒト・ゲノム配列の場合、総塩基数は約30億にも達するため、サフィックス・アレイ内のソート順位を表す整数には少なくとも４バイトが必要となり、サフィックス・アレイのサイズは12ギガバイトにも達する。そこで、サフィックス・アレイと情報的に等価な内容で、よりコンパクトは表現ができるコンプレスト・サフィックス・アレイ（非特許文献８）やバローズ・ホィーラー変換（非特許文献４）などを、サフィックス・アレイの代わりに用いて、上記と類似の構成で上記の目的を達することができる。
（実施例４）
前記実施例においては、A, G, C, T, Nの5種類の文字からなる塩基配列データに対して、本発明を適用し、mRNAの発現解析を行う方法を説明した。本発明は、より多種類の文字からなる一般的な文書データに対しても適用することが可能である。図8を用いて、本発明を用いて、大量の文字句データが大量の蓄積文書データ内部で扱われている話題に対して、どのような引用を行っているか解析するシステム全体の構成を説明する。より具体的には、例えば、大量の蓄積文書データとしては大量の広告文書を想定し、大量の文字句データとしては不特定多数の人から得られた大量のアンケート調査結果を想定し、どの文書のどの話題からの引用が多いかを解析すること想定してもよい。

111は、本発明による検索処理の計算を行う計算機である。112はテキストデータ入力装置であり、文字句データ113を取り込む。また、114は記憶装置であり、蓄積文書データ115と話題出現箇所データ116を保持している。計算機111は、テキストデータ入力装置112および記憶装置114とネットワーク117を介して接続されている118は計算機111に繋がれた表示出力装置であり、計算機111で計算された引用解析結果119を表示する。計算機111の内部で行われる処理の概要を図9に示す。入力処理部121では、文字句データ113、蓄積文書データ115、および、話題出現箇所データ116を入力する。引用箇所計算処理部122では、文字句データ113をクエリー配列データ21として、また、蓄積文書データ115をターゲット配列データ23として、図2を用いて説明した検索処理を行い、各クエリー配列に対するマッピング位置（本例では引用位置を意味する）31を計算して、それらの計算結果を文字句引用位置データ123として計算機111内部のメイン・メモリ124に保持する。図2を用いて説明した検索処理は、文字の種類数が増えても全く同様に行うことができる。出現箇所比較処理部125では、これらの文字句引用位置データ123と各話題の出現箇所データ116を比較し、文書内で話題が出現する箇所（開始位置から終了位置まで）の範囲内に引用位置が含まれるような文字句を、その話題に属する文字句であると判定する。引用カウント処理部124では、各話題に属する文字句の数を数える。これが、各話題への総引用数となる。出力処理部125では、そのカウント結果を引用解析結果119として出力する。

N. J. Larsson and K. Sadakane. Faster Suffix Sorting, Technical report LU-CS-TR:99-214, 1999. R. Gonzalez, S. Grabowski, V. Makinen, and G. Navarro. Practical Implementation of Rank and Select Queries. In Proc. WEA'05, pages 27-38, 2005. Roberto Grossi and Jeffrey Vitter, Compressed Suffix Arrays and Suffix Trees with Applications to Text Indexing and String Matching, SIAM Journal on Computing, 2005.

大量配列の一括検索において、計算機主記憶上に構築されるデータ構造とそれら相互の論理的関係、並びに、説明用の補助情報を表示した説明図である。 大量配列の一括検索において、本発明の動作の概略を示した説明図である。 大量配列の一括検索の部分処理である、検索文字列の1文字伸長処理の動作の概略を示した説明図である。 近似的に一致するターゲット配列内サフィックスの管理するためのデータ構造を示した説明図である。 近似的に一致するターゲット配列内サフィックスの1文字伸長処理の動作の概略を示した説明図である。 mRNA発現解析の計算を行うシステム全体の構成を示した説明図である。 mRNA発現解析の計算を行うシステムにおいて、計算機内部処理の手順を示した説明図である。 一般の文書データに対して引用解析を行うシステム全体の構成を示した説明図である。 一般の文書データに対して引用解析を行うシステムにおいて、計算機内部処理の手順を示した説明図である。

符号の説明

１ 現在処理中の検索文字列（塩基配列データ）
２ クエリー側サフィックス・アレイ
３ ターゲット側サフィックス・アレイ
４ 現在処理中の検索文字列（塩基配列データ）に対して、それが出現するクエリー側サフィックス・アレイ内のインデクス範囲、及び、ターゲット側サフィックス・アレイ内のインデクス範囲を示す、データ間の参照関係
５ Aにより１文字伸長された新たな検索文字列（塩基配列データ）
６ Cにより１文字伸長された新たな検索文字列（塩基配列データ）
７ データ参照関係が更新されたクエリー側サフィックス・アレイ
８ クエリー側サフィックス・アレイ内部への参照先の変化を表す状態遷移
９ データ参照関係が更新されたターゲット側サフィックス・アレイ
10 ターゲット側サフィックス・アレイ内部への参照先の変化を表す状態遷移
11 新たな検索文字列５に対して、それが出現するクエリー側サフィックス・アレイ内のインデクス範囲、及び、ターゲット側サフィックス・アレイ内のインデクス範囲を示す、データ間の参照関係
12 新たな検索文字列６に対して、それが出現するクエリー側サフィックス・アレイ内のインデクス範囲、及び、ターゲット側サフィックス・アレイ内のインデクス範囲を示す、データ間の参照関係
13 計算機主記憶上において、クエリー側の情報とターゲット側の情報を論理的に区分する説明用の分割線
21 クエリー配列データ（mRNAの塩基配列を表現する文字列データ）
22 クエリー側サフィックス・アレイの構築処理
23 ターゲット配列データ（ゲノムの塩基配列を表現する文字列データ）
24 ターゲット側サフィックス・アレイの構築処理
25 検索文字列とそのクエリー側サフィックス・アレイ内の出現範囲、および、ターゲット側サフィックス・アレイ内の出現範囲の初期化処理
26 検索文字列の1文字伸長処理および棄却判断
27 該当するクエリー配列の最長一致サフィックスの報告
28 全長一致するクエリー配列をもつ検索文字列の抽出処理
29 検索文字列ごとの繰り返し処理
30 検索文字列のターゲット配列内の出現回数の評価と判断
31 該当するクエリー配列のターゲット配列内のマッピング位置の報告
32 該当するクエリー配列をリピート配列として報告
41 「検索文字列の1文字伸長処理および棄却判断」26の入力データ
42 文字種（A, G, C, T）ごとの繰り返し処理
43 検索文字列の1文字伸長処理
44 新たな検索文字列に対するクエリー側サフィックス・アレイ内出現位置の更新
45 新たな検索文字列に対するクエリー側サフィックス・アレイ内の出現有無の判定
46 新たな検索文字列の棄却処理
47 新たな検索文字列に対するターゲット側サフィックス・アレイ内出現位置の更新
48 新たな検索文字列に対するターゲット側サフィックス・アレイ内の出現有無の判定
49 全長一致するターゲット内の部分配列をもたないようなクエリー配列に対する、最長一致サフィックスの報告
50 新たな検索文字列の棄却処理
51 「検索文字列の1文字伸長処理および棄却判断」26の出力データ
61 検索文字列S
62 Sと近似的に一致する複数種類のターゲット配列内サフィックスQの管理テーブル
63 ターゲット側サフィックス・アレイ
71 Sと近似的に一致するサフィックスQの1文字伸長処理の入力データ
72 文字種ごとの繰り返し処理
73 サフィックスQの伸長処理
74 新たな検索文字列yQに対するターゲット側サフィックス・アレイ内出現位置の更新
75 出現の有無の判定
76 yQの棄却処理
77 塩基置換による不一致数またはNの数の増加の判定
78 塩基置換またはNによるyQの不一致数の更新
79 塩基欠失によるQの不一致数の更新
80 塩基挿入のための文字種ごとの繰り返し処理
81 サフィックスyQの伸長処理
82 新たな検索文字列zyQに対するターゲット側サフィックス・アレイ内出現位置の更新
83 出現の有無の判定
84 zyQの棄却処理
85 塩基挿入によるQの不一致数の更新
86 不一致数の判定
87 サフィックスの棄却
88 Sと近似的に一致するサフィックスQの1文字伸長処理の出力データ。

Claims (4)

1. 複数の文字配列からなる集合をクエリーデータとし、夫々の前記文字配列より長い長さをもつ文字配列からなる集合をターゲットデータとして入力手段を介して入力し、
   前記クエリーデータに属する各文字配列が、前記ターゲットデータ中の配列に部分配列として出現するか否か、部分配列として出現する回数、部分配列として出現する位置を算出する複数文字配列の一括検索方法において、
   前記クエリーデータに属する全ての文字配列のサフィックスを辞書式順番にソートした第一インデクスデータと、前記ターゲットデータに属する全ての文字配列のサフィックスを辞書式順番にソートした第二インデクスデータとを記憶手段に記憶し、
   同一の作業配列が複数のクエリー配列のサフィックスとして共有される事例のリストを前記第一インデクスの範囲とし、前記同一の作業配列がターゲット配列中に出現する全ての位置のリストを前記第二インデクスの範囲として、該第一、第二のインデクス範囲から両インデクス間の対応関係求め、
   作業配列の長さを０から1文字ずつ伸長しながら両インデクス間の対応関係を求め、
   作業配列の長さがクエリーデータの配列の長さに達するまで前記対応関係を求める処理を繰り返し、作業配列と全長にわたり一致する全てのクエリー配列とそれらのターゲット配列中の全ての出現位置との対応関係を両インデクス間の対応関係として求めることを特徴とする複数文字配列の一括検索方法。
2. 請求項１に記載の複数文字配列の一括検索方法において、
   ２つの文字配列が予め指定した許容数以下の箇所の文字の置換や挿入や欠失により一致するとき、または他の任意の文字と一致する特殊文字が予め指定した許容回数以下だけ含まれるとき、それら２つの文字配列が近似的に一致するものとし、
   前記クエリーデータに属する各文字配列と近似的に一致する配列が、前記ターゲットデータ中の配列に部分配列として出現するか否か、部分配列として出現する回数、部分配列として出現する位置を算出することを特徴とする複数文字配列の一括検索方法。
3. 請求項１に記載の複数文字配列の一括検索方法において、
   前記クエリーデータに属する全ての配列のサフィックス、または、ターゲットデータに属する全ての配列のサフィックスに対して、それらを辞書式順番にソートしたインデクスデータを前記記憶手段に記憶する代わりに、
   それと情報論的に等価なバローズ・ホィーラー変換、サフィックス・トリー、または、コンプレスト・サフィックス・アレイを前記記憶手段に記憶することを特徴とする複数文字配列の一括検索方法。
4. cDNA配列データ、ゲノム配列データ、および、遺伝子領域座標データを入力する入力処理部と、
   入力された前記cDNA配列データをクエリー配列データとし、前記ゲノム配列データをターゲット配列データとし、前記クエリー配列データに属する全ての文字配列のサフィックスを辞書式順番にソートした第一インデクスデータと、前記ターゲット配列データに属する全ての文字配列のサフィックスを辞書式順番にソートした第二インデクスデータとを記憶装置に記憶し、
   同一の作業配列が複数のクエリー配列のサフィックスとして共有される事例のリストを前記第一インデクスの範囲とし、前記同一の作業配列がターゲット配列中に出現する全ての位置のリストを第二インデクスの範囲とし、
   これらのインデクス範囲の組により両インデクス間の対応関係を求め、作業配列の長さを０から1文字ずつ伸長しながら両インデクス間の対応関係を求め、
   作業配列の長さがクエリー配列データの配列の長さに達するまで前記対応関係を求める処理を繰り返し、作業配列と全長にわたり一致する全てのクエリー配列とそれらのターゲット配列中の全ての出現位置との対応関係を一括して両インデクス間の対応関係として求めることにより、前記クエリー配列データに属する各文字配列が、前記ターゲット配列データ中の配列に部分配列として出現するか否か、部分配列として出現する回数、部分配列として出現する位置を算出して検索処理を行い、各クエリー配列に対するマッピング位置を計算するマッピング位置計算処理部と、
   前記各クエリー配列のマッピング位置と前記遺伝子領域座標データを比較し、遺伝子領域と同じ染色体上でその開始位置から終了位置までの範囲内にマッピング位置が含まれるようなクエリー配列を、その遺伝子に属するクエリー配列であると判定する座標データ比較処理部と、
   前記遺伝子に属するクエリー配列の数を数える配列カウント処理部と、
   前記配列カウント処理部での処理結果を発現解析結果として出力する出力処理部とを有することを特徴とする検索システム。

Images