WO2010119783A1 - ペア文字列検索システム - Google Patents

ペア文字列検索システム Download PDF

Info

Publication number
WO2010119783A1
WO2010119783A1 PCT/JP2010/056146 JP2010056146W WO2010119783A1 WO 2010119783 A1 WO2010119783 A1 WO 2010119783A1 JP 2010056146 W JP2010056146 W JP 2010056146W WO 2010119783 A1 WO2010119783 A1 WO 2010119783A1
Authority
WO
WIPO (PCT)
Prior art keywords
character string
index
pair
processing unit
sequence
Prior art date
Application number
PCT/JP2010/056146
Other languages
English (en)
French (fr)
Inventor
木村　宏一
Original Assignee
株式会社日立製作所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社日立製作所 filed Critical 株式会社日立製作所
Priority to US13/264,037 priority Critical patent/US8788522B2/en
Priority to JP2011509261A priority patent/JP5279897B2/ja
Publication of WO2010119783A1 publication Critical patent/WO2010119783A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/90—Details of database functions independent of the retrieved data types
    • G06F16/903—Querying
    • G06F16/90335—Query processing
    • G06F16/90344—Query processing by using string matching techniques
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search

Definitions

  • the present invention relates to a method for quickly searching two character strings appearing in proximity from large-scale genome sequence data or a large amount of general document data using a computer.
  • the present invention relates to a method for specifying a corresponding position on the genome at a high speed with respect to a large amount of DNA sequence data of a massively parallel DNA sequencer obtained by using the pair-end method.
  • Non-patent Document 1 a massively parallel DNA sequencer based on a new principle, which is completely different from the conventional capillary DNA sequencer, has appeared.
  • a single run of a massively parallel DNA sequencer can read a large number of sequences, reaching tens of millions at a time.
  • the pair-end method is used to compensate for the disadvantages (Non-Patent Document 3).
  • genomes are sequenced by sequencing several tens of bases from both ends of a large number of genomic sequence fragments that are aligned to a substantially constant length (a few hundred bases to several kilobases). Information on pairs of base sequences (pair-end sequences) that are separated by a substantially constant interval can be obtained.
  • the genome mapping calculation for the reference genome sequence is performed on the paired-end sequence data obtained by analyzing the genomic DNA sample by the sequencer, that is, the sequence obtained by the sequencing is located at any position on the reference genome sequence.
  • the expected interval such as insertion / deletion between the sample genome and the reference genome Structural polymorphism can be detected.
  • the distance between mapping positions is larger than the expected distance, it is considered that there is a deletion between these sequence pairs on the sample / genome side, and conversely, the distance between mapping positions is expected. If it is smaller than the interval, the insertion between the sequence pairs is considered to occur on the sample genome side (Non-patent Document 3).
  • mapping position is not uniquely determined in the genome mapping calculation, and many positions may be listed as candidates. Even in such a case, by using the mapping position interval of the paired arrays as a constraint, that is, by using the constraint that these mapping positions are close to each other, can the mapping position as the paired array be uniquely determined? Alternatively, it can be expected that the number of candidates can be narrowed down.
  • introns may enter between mapping positions of paired sequences, and the mapping position interval may be increased by the length of the intron. is there. Even in such a case, the restriction that the interval between mapping positions does not exceed the length of the gene region on the genome can be used in mapping calculation as a pair sequence.
  • Non-Patent Document 2 When performing genome mapping calculation of a large amount of sequence sequence data (query sequence data) on a large-scale reference genome sequence, the reference genome sequence data is usually indexed in advance. The search for the query sequence can be accelerated by using this index. As an indexing method, a suffix array (Non-Patent Document 2) can be used.
  • MAQ is known as software capable of mapping a large amount of pair sequence data of a massively parallel DNA sequencer analyzed by the pair-end method onto a reference genome sequence with high accuracy and high speed (Non-patent Document 6).
  • index information is created in advance for a large amount of pair sequence data, and mapping position candidates are searched by referring to the index information while scanning the reference genome sequence.
  • the candidate for the position to be mapped on the plus strand of the genome sequence is assigned to each query sequence in the scanning process.
  • mapping position is obtained at high speed using an index independently for each sequence
  • a combination satisfying the space constraint may be selected from the combinations of the mapping positions.
  • the sequence length is relatively short
  • mapping positions are obtained independently for two pairs of sequences without using space constraints
  • a large number of candidate positions are obtained first, after which A small number satisfying the space constraint is selected from among the above, and the processing speed may be reduced.
  • an oversight may occur depending on their relative positional relationship.
  • An object of the present invention is that when a large-scale reference genome sequence data and a large amount of pair-end sequence sequence data (query sequence data) are given, the mapping position interval between the pairs is equal to or less than a predetermined value. It is to obtain a combination of mapping positions of such pair arrangements at high speed.
  • the entire genome sequence has information equivalent to that of a normal suffix array. Furthermore, even for partial regions of the genome restricted at various scales, index information restricted to that region (local index information) Index of the genome sequence data is performed so that Further, a method is used in which the index information for a certain region of the genome is divided into index information for each partial region when the region is divided into two. For paired sequences, first, the entire genome region is searched independently for each sequence to obtain index information, and then, by dividing into two until the length of the region becomes 1, those sequences are obtained. The position coordinates on the genome for the index are gradually determined in detail. In the process, the distance restriction between pairs is taken into consideration, and index information that does not satisfy the restriction is removed.
  • a pair character string search system includes a storage device that stores a reference character string, a first character string and a second character string that form a pair, and a distance between the first character string and the second character string.
  • An input device for inputting an upper limit value, a suffix array construction processing unit for constructing a suffix array of reference character strings, and an index in the suffix array independently of the first character string and the second character string
  • An index area search processing unit that calculates a region, an LSA construction processing unit that constructs index information LSA that can be localized, an index region localization processing unit, a distance evaluation processing unit between index regions, and a distance within an upper limit value
  • An output device that outputs a position of the paired first character string and the second character string on the reference character string, and the index area localization processing unit Distribution, and the results of evaluation by inter index range distance evaluation processing unit, removes the information of the index region which can not make a pair at a distance of more than the upper limit.
  • the LSA construction processing unit uses the suffix array of the reference character string as an initial block, and divides the suffix array into two blocks corresponding to the first and second partial regions of the reference character string according to whether the most significant bit is 0 or 1 Then, the arrangement of the most significant bits referred to at that time is used as the first column of the LSA, and the obtained two blocks are further divided into two blocks according to whether the second bit is 0 or 1. As a result, the entire reference genome sequence is divided into four blocks corresponding to the partial region to be divided into four, the second-order bits referenced at that time are arranged as the second column of LSA, and so on. The process is repeated until the length of 1 becomes 1, and in the process, each column of the LSA is calculated to construct the LSA.
  • the index area localization processing unit uses the index area in the suffix array for each of the first character string and the second character string as the index area in the initial block, and corresponds to the first half and second half partial areas of the reference character string.
  • the index area in the initial block is divided into the index areas in each block using the information in the first column of the LSA, and the entire reference character string is divided into four
  • the index area in the two blocks is divided into the index areas in the four blocks using the information in the second column of the LSA corresponding to the division into the four blocks corresponding to, and the same division process is repeated thereafter.
  • a calculation is performed to localize the index area in the initial block to the index area in the partial area of the reference character string.
  • the index area localization processing unit applies a rank function to the blocks in each column of the LSA, and calculates the position of the index area in the two blocks obtained by dividing the block.
  • the inter-index-area-distance evaluation processing unit assumes that the position coordinates on the reference character string of the index area localized in the partial area leave an indeterminate range that is only the length of the partial area. The upper limit of the distance between the index areas of the converted first character string and second character string is evaluated.
  • the high-speed operation is performed while considering the restriction of the mapping position interval.
  • Index search mapping calculation
  • Example 1 a method for obtaining a mapping position close to a specified interval or less when a pair of reference genome sequence data and a query sequence is given will be described.
  • the mapping between the query sequence and the genome sequence is considered to be an exact match, but application to the incomplete match that allows a mismatch of several bases is possible by simple analogy.
  • only the upper limit of the interval between the mapping positions of the pair array is considered, but simultaneously, the lower limit of the interval between the mapping positions of the pair array can be considered by simple analogy.
  • FIG. 1 is a diagram illustrating a configuration example of a system that performs paired array mapping calculation according to the present embodiment.
  • Reference numeral 101 denotes a computer in which this system is operating, which reads the reference genome sequence data 103 stored in the external storage device 102 and analyzes the DNA sample 115 using the DNA sequencer 116, and the pair sequence data obtained A (query sequence pair) 107 is read via the input device 106, and a parameter 114 (upper limit value D between pairs) specified by the user is read via the input device 118.
  • the computer 101 uses the suffix array construction processing unit 105 to connect the base sequences in the reference genome sequence data 103 via the delimiter $, calculates the suffix array, and stores it in the main memory 104.
  • the LSA construction processing unit 108 constructs index information LSA (localized suffix array) based on the suffix array.
  • the index area search processing unit 109 searches the suffix array independently for each of the paired query sequences, and obtains index area information.
  • the index area localization processing unit 110 localizes (subdivides) an index area that satisfies the restriction of the mapping position interval between pairs based on the index area information of the query sequence and the LSA of the reference genome sequence. Then, the position of the index region on the genome is refined from a state where it spreads throughout the initial genome to a state where it is finally determined in detail at the base level.
  • mapping position information (search result 113) to the output device 112.
  • the mapping position information (search result 113) of the pair sequence includes information on the chromosome name, its direction (+/- which strand) and coordinates indicating the appearance position.
  • Non-patent Document 2 In order to construct a suffix array from a reference genome sequence, a known method such as Manber-Myers may be used (Non-patent Document 2). Further, in order to obtain the index information of the query sequence, that is, the index range in which the query sequence appears as a suffix prefix in the suffix array, a general binary search method may be used.
  • each suffix is represented by an integer not less than 0 and less than 2 to the power of h.
  • An empty string is the final suffix.
  • the index of the suffix corresponding to the first half of the genome is represented by an integer with the most significant bit (h bit from the least significant bit) being 0 in binary notation, and the index of the suffix corresponding to the second half of the genome is 2 When expressed in hexadecimal, the most significant bit (hth bit from the least significant bit) is represented by an integer of 1.
  • the suffix array B for the whole genome region is divided into two by the most significant bit being 0 or 1 while maintaining the internal relative order, the suffix array subsequence B0 corresponding to the first half of the genome, A subsequence B1 of the suffix array corresponding to the second half of the genome is obtained.
  • the same sequence as the query sequence appears scattered on the first half 205 of the genome sequence and the second half 207 of the genome sequence as shown by 206 and 208, respectively.
  • 210 and 212 appear together in B0 and B1, respectively.
  • the same bisection can be performed by checking whether the value of the second bit from the top (h-1 bit from the bottom) of the index is 0 or 1. , Thereby dividing the whole genome into four parts. Subsequently, 8 divisions, 16 divisions, and the like are possible, and the division can be repeated recursively until the length of the area finally reaches 1.
  • Reference numeral 301 represents the entire LSA
  • 302 represents the most significant bit of the index in B that was referenced when the entire genome region was divided into the first half and the latter half in the first two divisions.
  • the next four divisions that is, in the two divisions of the above two areas, the second bit from the higher order of the indexes referenced in B0 and B1 is arranged.
  • 304 is a table in which the third bit from the top of the index referenced from each area after the four divisions is arranged vertically.
  • the fourth and subsequent columns are similarly configured for 16 divisions and thereafter.
  • FIG. 4 is a diagram for explaining a method for evaluating the distance on the genome between index regions included in different blocks in the inter-index region distance evaluation processing unit 117. As will be described later, this evaluation method is used when configuring the index region localization processing unit 110.
  • 401 represents a block. Assume that index areas 402 and 403 for two query sequences p and q are included in different blocks 404 and 405. The position coordinates on the genome of these index regions are fixed, leaving uncertain ranges as indicated by 406 and 407. The block length is calculated based on the total length of the genome and the number of divisions. Therefore, the upper limit 409 and the lower limit 408 of the distance between these index regions can be calculated by counting the number of blocks between the block 404 and the block 405.
  • FIG. 5 shows that when index regions in the whole genome for two query sequences p and q are given, the position coordinates on the genome sequence as a pair sequence recursively divide the genome region in the index region localization processing unit 110. It is a figure explaining how it decides while reducing uncertainty by repeating automatically.
  • Reference numeral 501 denotes a suffix array for the entire genome region (initial block).
  • Reference numeral 502 denotes a partial array of the suffix array corresponding to the first and second half blocks when the whole genome is divided into two.
  • Reference numeral 503 denotes a partial array of suffix arrays corresponding to each block when the whole genome is divided into four.
  • reference numerals 504 and 505 are obtained by vertically arranging the partial columns of the suffix array corresponding to each block when the whole genome is divided into 8 and 16, respectively.
  • Reference numeral 506 denotes a case where the block length finally reaches 1. Solid arrows connecting between 501 to 506 indicate that the index area for the query sequence p is detailed by recursive localization processing (block segmentation processing).
  • a broken-line arrow connecting between 501 to 506 indicates that the index area for the query sequence q is detailed by the localization process.
  • processing that does not satisfy the constraint on the distance between paired sequences is aborted.
  • those that do not satisfy the constraint on the distance between the paired sequences there are those in which the order of the paired sequences is reversed, such as 509.
  • 510 and 511 remaining up to 506 give the mapping position of the pair sequence.
  • Fig. 6 shows the LSA construction processing procedure.
  • block B is set as a suffix array for the entire reference genome sequence calculated by the suffix array construction processing unit 105 and stored in the main memory 104.
  • B is an array in which indexes (h-bit integers) are arranged to the power of 2h.
  • Step 604 is the end determination process.
  • the processing of steps 606 and 607 is performed for each block represented by Ba with a character string a of length (i-1) consisting of 0 and 1.
  • step 606 for each index (h-bit integer) included in the block Ba, a bit string formed by arranging the i-th bit from the higher order is defined as La.
  • each index included in Ba is classified into block Ba0 or Ba1, depending on whether those bits are 0 or 1. At that time, the relative order between the indexes is maintained.
  • the value of i is incremented by 1 in step 608 and the process returns to immediately before step 604. If the end condition is satisfied in step 605, LSA LSA is constructed on the main memory 104 in step 609.
  • FIG. 7 is an explanatory diagram of a method for constructing an LSA from La.
  • the LSA 701 is an array of bits in 2 h rows and h columns. Its first column 702 is equal to L.
  • the second column 703 is obtained by connecting L0 and L1.
  • the third column 704 is obtained by connecting L00, L01, L10, and L11.
  • the fourth column 705 is obtained by connecting L000, L001, L010, L011, L100, L101, L110, and L111. The same applies to the fifth and subsequent columns.
  • Reference numeral 801 denotes a list (including an empty character string) of all suffixes of “abracad” and their indexes.
  • Reference numeral 802 denotes a suffix array in which they are sorted in lexicographic order.
  • Reference numeral 803 denotes a suffix array index value in binary notation, and B corresponding to the entire genome.
  • Reference numeral 804 denotes an LSA constructed from B by the method described in FIG. 6, and corresponds to FIG.
  • FIG. 9 is an explanatory diagram for explaining a process of constructing an LSA from B corresponding to the entire genome according to the method described in FIG. 6 when the reference genome sequence is “abracad”.
  • 901 is B corresponding to the whole genome.
  • Reference numeral 902 is a list in which B0 and B1 corresponding to the first half and the second half of the genome are separated by a broken line. When dividing from B to B0 and B1, the most significant bit (bold) of B is referred to. The arrangement of these bits is L0, which is the first column of the LSA indicated by 905.
  • reference numeral 903 denotes B00, B01, B10, and B11 corresponding to the four divisions of the genome, which are separated by a broken line.
  • L0 and L1 are a sequence of these bits, which is the second column of LSA indicated by 905.
  • reference numeral 904 shows B000, B001, B010, B011, B100, B101, B110, and B111 corresponding to 8 divisions of the genome, which are separated by a broken line.
  • L00, L01, L10, and L11 are a sequence of these bits, which is the third column of the LSA indicated by 905. At 904, the block length becomes 1, and the recursive division ends here. Therefore, the 8-by-3 bit matrix indicated by 905 is LSA.
  • FIG. 10 is a diagram for explaining the internal processing flow of the index area localization process 110.
  • the upper limit value D of the distance is read from the user via the input device 106.
  • the appearance range (s p , t p ) of the query sequence p in the suffix array B of the entire genome sequence and the query sequence q The appearance range (s q , t q ) is read from the main memory 104.
  • B has already been constructed by the suffix array construction processing unit 105, and the values of (s p , t p ) and (s q , t q ) have been calculated by the index area search unit 109, and these are stored in the main memory 104. Is held in.
  • information List (p) indicating a candidate list of mapping positions of query sequence p
  • information List (q) indicating a candidate list of mapping positions of query sequence q are initialized as shown in the following equation.
  • Step 1004 is an end determination process of this iterative process. When the block length becomes 1, it is determined to end.
  • step 1005 according to the method described below with reference to FIG. 11, each element belonging to the list List (p) is subdivided to update the list List (p), and similarly each element belonging to List (q) is subdivided. And list List (q) is updated.
  • the mapping positions of the index areas of each query sequence are determined with the uncertainty of only the block length, and based on this, the distance between the mapping positions of the index areas belonging to different blocks is determined.
  • the upper and lower limits can be evaluated.
  • the elements of the lists List (p) and List (q) are evaluated in combination, and the lists List (p) and List (p) such that a pair whose upper limit of the distance between the pairs is not more than D cannot be formed. Delete the element of q).
  • the lists List (p) and List (q) are returned to the main memory 104 as answers in step 1007.
  • FIG. 11 is a diagram for explaining the processing procedure of the list segmentation processing 1005 for updating the information List_in indicating the candidate list of the mapping position of the query sequence and replying with a new list List_out.
  • This processing procedure is the same as the method described with reference to FIG.
  • a list List_in is input, and the list List_out is set as an empty list.
  • List_in is List (p) or List (q) before being subdivided in step 1005.
  • the following processing is performed on each element (Ba, s, t) of the list List_in.
  • a represents a character string consisting of 0 and 1
  • s and t are integers indicating the index area in the block Ba of the query sequence.
  • integers0, t0, s1, and t1 indicating the index areas in the subdivided blocks Ba0 and Ba1 are calculated according to the following equation.
  • rank (La, s, 0) rank (La, t, 0)
  • t0 rank (La, t, 0)
  • s1 rank (La, s, 1)
  • t1 rank (La, t, 1)
  • rank (La, s, 0) is a rank function and represents the number of 0 appearing before the position of the s-th row in the bit string La in the LSA
  • rank (La, s, 1) is also a rank.
  • the function represents the number of 1 appearing before the position of the s-th row in the bit string La in the LSA.
  • the rank function can be efficiently calculated by a known method (Non-Patent Document 5). The same applies to rank (La, t, 0) and rank (La, t, 1).
  • step 1104 depending on whether each bit in the corresponding bit string La in LSA is 0 or 1, (Ba, s, t) is divided into two blocks (Ba0, s0, t0) and (Ba1, s1, t1). And add them to the list List_out.
  • step 1105 the obtained list List_out is returned in the main memory 104 as an answer. This becomes List (p) or List (q) after the subdivision in step 1005.
  • mapping position candidates for the query sequence are compactly expressed as index areas (two integers s and t) in each block. Further, by repeatedly performing the division of the block and the accompanying division of the index area, the mapping position of the index area can be finally determined in detail up to the base level while gradually reducing the uncertainty range. At that time, the upper limit of the distance between the mapping positions leaving the indeterminate range is evaluated, and the index information for the mapping positions that cannot satisfy the mapping position interval restriction (distance D or less) specified between the pairs is obtained. Remove and avoid useless calculations. Also, by referring to the LSA, the index area division accompanying the block division can be calculated efficiently. A procedure for efficiently building LSA was also given.
  • Example 2 In the above example, the method for calculating the genomic mapping of the query sequence obtained by the paired-end method was shown, but the query obtained by sequencing the fragment sequence of the cDNA sequence made from the spliced mRNA.
  • the present invention can also be implemented for the problem of calculating genome mapping of sequences.
  • FIG. 12 is a diagram illustrating a configuration example of a system that performs mapping calculation including splicing according to the present embodiment.
  • Reference numeral 1201 denotes a computer in which this system is operated, and reads the reference genome sequence data 1203 stored in the external storage device 1202 and also obtains the query sequence obtained by analyzing the cDNA sample 1216 using the DNA sequencer 1217.
  • the data 1207 is read via the input device 1206, and the parameter 1219 (intron length upper limit value D) designated by the user is read via the input device 1218.
  • the computer 1201 uses the suffix array construction processing unit 1205 to connect the base sequences in the reference genome sequence data 1203 via the delimiter $, and calculates the suffix array and stores it in the main memory 1204. Similarly, calculation results of other processing units are also stored in the main memory 1204, and the contents are referred to by each processing unit as necessary.
  • the LSA construction processing unit 1208 constructs new index information LSA (localized “suffix” array) based on the suffix array.
  • the pair array generation processing unit 1209 divides the query array into two at various positions and creates various pair arrays.
  • the index area search processing unit 1210 searches the suffix array independently for each pair array to obtain index information.
  • coordinates corresponding to an index satisfying the constraint that the mapping position interval between pairs is equal to or less than the upper limit of the intron length based on the index information of the pair sequence and the LSA of the reference genome sequence. Is determined while gradually localizing from the entire genome. Thereby, the mapping position as a pair arrangement
  • the consensus sequence inspection processing unit 1212 examines the genome sequence around the splice site and removes mapping results that do not satisfy the consensus.
  • the output processing unit 1213 outputs the mapping position information (search result 1215) to the output device 1214.
  • the mapping information (search result 1215) of the spliced sequence includes information on the chromosome name, its direction (+/ ⁇ which strand) and coordinates indicating two separated appearance positions.
  • the processing contents of the suffix array construction processing unit 1205, the LSA construction processing unit 1208, the index region search processing unit 1210, the index region localization processing unit 1211, and the like are the same as those in the first embodiment.
  • the pair array generation processing unit 1209 specifies in advance a lower limit value of the array length of the pair array obtained by dividing, and generates all pair arrays that satisfy the condition. They are processed as paired sequence data by the index region search processing unit 1210 and the index region localization processing unit 1211, but are interpreted and checked as a spliced mapping by the consensus sequence inspection processing unit 1212.
  • the part sandwiched between the genome mapping positions as a pair sequence is interpreted as an intron sequence, and the condition that the intron sequence starts with GT and ends with AG is imposed. Mapping results of pair sequences that do not satisfy this condition are removed.
  • mapping position interval between pairs is equal to or less than a predetermined value.
  • Example 3 In the above embodiment, the embodiment for the base sequence data composed of A, C, G, T, and N has been described.
  • the present invention is also applicable to general character string data that is not limited to base sequence data.
  • amino acid fragment sequence data will be described.
  • the amino acid sequence is a sequence composed of about 20 kinds of amino acids.
  • FIG. 13 is a diagram showing a configuration example of a system that performs calculation for identifying a protein from the amino acid fragment sequence of this example.
  • Reference numeral 1301 denotes a computer in which this system is operating.
  • the computer 1301 reads the reference protein sequence data 1303 stored in the external storage device 1302 and analyzes the protein sample 1306 using the mass spectrometry system 1305.
  • the query sequence data 1307 is read via the input device 1308, and the parameter 1318 (upper limit value D of the protein sequence length) designated by the user is read via the input device 1319.
  • the computer 1301 uses the suffix array construction processing unit 1310 to connect the amino acid sequences in the reference protein sequence data 1303 via the delimiter $, calculate the suffix array, and store it in the main memory 1309.
  • calculation results of other processing units are also stored in the main memory 1309, and the contents thereof are referred to by each processing unit as necessary.
  • Processing contents in the suffix array construction processing unit 1310, the LSA construction processing unit 1311, the index region search processing unit 1312, the index region localization processing unit 1313, and the like are the same as those in the above embodiment.
  • FIG. 14 is a diagram illustrating a configuration example of a system that performs calculation for identifying a document in which two search character strings appear close to each other in general document data according to the present embodiment.
  • Reference numeral 1401 denotes a computer in which this system is operating, which reads the reference document data 1403 stored in the external storage device 1402, and reads a search character string pair (query) 1407 via the input device 1406. Also, the parameter 1411 (upper limit value D of the distance between pairs) designated by the user is read via the input device 1416.
  • the reference document data 1403 is a collection of documents consisting of character string bodies with document names. However, even for general structured documents having a hierarchical structure, simple reference is made. The present invention can be applied.
  • the calculator 1401 uses a suffix array construction processing unit 1405 to connect character strings in the document data 1403 via identifiable delimiters, and calculates the suffix array and stores it in the main memory 1404. Similarly, calculation results of other processing units are also stored in the main memory 1404, and the contents thereof are referred to by each processing unit as necessary.
  • the LSA construction processing unit 1408 constructs index information LSA (localized suffix array) based on the suffix array.
  • the index area search processing unit 1409 searches the suffix array independently for each search character string to obtain index information. In the index area localization processing unit 1410, based on the index information of the search character string and the LSA of the document data, the appearance position in the document with respect to the index that appears close to the document data 1403 is gradually increased from the entire document.
  • the document identification processing unit 1412 checks whether two search character strings appear in the same document in the document data.
  • the output processing unit 1413 outputs to the output device 1414 a search result 1415 including the document name of the document in which the two search character strings appear close to each other and information on the appearance position of each search character string.
  • the occurrence positions of character string pairs such that the interval between the appearance positions between the pairs is equal to or less than a predetermined value shown in the above-described embodiment.

Abstract

　計算機上でペア文字列検索を高速に行うためのインデクス情報のデータ構造を提供し、また、そのインデクス情報を用いて文書中に近接して現れるペア文字列を高速に検索する方法を提供する。　参照文書データのサフィックス・アレイのビットを並べ替えて、局所化可能な、即ち、文書の部分領域に対するインデクスとしても利用可能な、インデクス情報ＬＳＡを作成する。これを利用して、文書全体を初期領域とする領域の二分割処理を繰り返して、クエリー文字列に対するインデクス情報の参照文書データ中の位置を段階的に詳細化し、ペア間の距離を評価して候補を絞り、最終的にペア文字列が近接して出現する位置を確定する。

Description

ペア文字列検索システム
　本発明は、大規模なゲノム配列データ、又は大量の一般の文書データの中から、近接して現れる二つの文字列を、計算機を利用して高速に検索するための方法に関する。特に、ペア・エンド法を用いて得られた超並列ＤＮＡシーケンサの大量ＤＮＡ配列データに対して、対応するゲノム上の位置を高速に特定するための方法に関する。
　近年、従来型のキャピラリー型ＤＮＡシーケンサとは全く異なる、新しい原理に基づいた超並列ＤＮＡシーケンサが出現した（非特許文献１）。超並列ＤＮＡシーケンサの一回のランでは、一度に数千万本にも及ぶ大量の配列を読み取ることができる。しかし、読取り配列長が数十塩基長程度と比較的短いという短所がある。そこで、その短所を補うためにペア・エンド法が用いられる（非特許文献３）。ペア・エンド法では、ほぼ一定の長さ（百数十塩基から数キロ塩基程度）に揃えられた多数のゲノム配列断片に対して、その両端から数十塩基の部分をシーケンスすることにより、ゲノム上でほぼ一定の間隔だけ離れた塩基配列のペアの情報（ペア・エンド配列）を得ることができる。
　ゲノムＤＮＡサンプルをシーケンサにより解析して得られたペア・エンド配列データに対して参照ゲノム配列に対するゲノム・マッピング計算を行い、即ち、それらシーケンスして得られた配列が参照ゲノム配列上のどの位置に現れるかの計算を行い、ペアをなす配列の参照ゲノム上のマッピング位置間の距離が期待される間隔だけ離れているかを調べることにより、サンプル・ゲノムと参照ゲノムの間の挿入・欠失などの構造多型を検出できる。即ち、マッピング位置間の距離が期待される間隔より大きければ、サンプル・ゲノム側にはそれらの配列ペア間に欠失が生じていると考えられ、また、逆に、マッピング位置間の距離が期待される間隔より小さければ、サンプル・ゲノム側にはそれらの配列ペア間に挿入が生じていると考えられる（非特許文献３）。
　シーケンス配列長が短い場合には、ゲノム・マッピング計算において、マッピング位置が一意に定まらず、多数の位置が候補として挙がることがある。そのような場合でも、ペアをなす配列のマッピング位置間隔を制約として用いることにより、即ち、それらのマッピング位置は近接しているという制約を用いることにより、ペア配列としてのマッピング位置を一意に決めるか、或いは、少数の候補に絞り込むことができることが期待できる。　また、転写産物サンプルをシーケンサにより解析して得られたペア・エンド配列データに関しては、ペアをなす配列のマッピング位置の間にイントロンが入り込み、マッピング位置の間隔がイントロン長の分だけ長くなることもある。その場合でも、マッピング位置間の間隔は、ゲノム上の遺伝子領域の長さを超えないという制約を、ペア配列としてのマッピング計算の際に利用することができる。
　大規模な参照ゲノム配列に対して、大量のシーケンス配列データ（クエリー配列データ）のゲノム・マッピング計算を行う際には、通常、参照ゲノム配列データを予めインデクス化しておく。クエリー配列の検索は、このインデクスを利用することにより、高速化できる。インデクス化の方法として、サフィックス・アレイ（非特許文献２）を利用できる。
　ペア・エンド法で解析された超並列ＤＮＡシーケンサの大量のペア配列データを参照ゲノム配列上に精度良く高速にマッピングすることが可能なソフトウェアとして、ＭＡＱが知られている（非特許文献６）。ＭＡＱでは、大量のペア配列データに対して予めインデクス情報を作成しておき、参照ゲノム配列上をスキャンしながらインデクス情報を参照してマッピング位置の候補を検索する。ペアをなす２つの配列が距離の制約を満たして近接してマッピングされるという条件を考慮するためには、スキャンの過程では、ゲノム配列のプラス鎖上にマッピングされる位置の候補を各クエリー配列当り最新の２か所まで計算機の記憶領域内に保持し、同時にマイナス鎖上にマッピングされる位置の候補を検索して、見つかったマイナス鎖上のマッピング位置が保持されたプラス鎖上のマッピング位置と距離の制約を満たすペアを作れるか否かを検査する。これにより、ペアをなす２つの配列のマッピング位置の候補の組み合わせの評価を効率的に行っている。
Service RF. "Gene sequencing. The race for the $1000 genome." Science. 2006 Mar 17;311(5767):1544-6. Manber, U. and Myers, G.: Suffix arrays: A new method for on-line string searches, in 1st ACM-SIAM, Symposium on Discrete Algorithms, pp. 319-327 (1990). Jan O. Korbel, Alexander Eckehart Urban, Jason P. Affourtit, Brian Godwin, Fabian Grubert, Jan Fredrik Simons, Philip M. Kim, Dean Palejev, Nicholas J. Carriero, Lei Du, Bruce E. Taillon, Zhoutao Chen, Andrea Tanzer, A. C. Eugenia Saunders, Jianxiang Chi, Fengtang Yang, Nigel P. Carter, Matthew E. Hurles, Sherman M. Weissman, Timothy T. Harkins, Mark B. Gerstein, Michael Egholm, and Michael Snyder, Paired-End Mapping Reveals Extensive Structural Variation in the Human Genome, Science 19 October 2007: 420-426. Ross Lippert, Space-efficient whole genome comparisons with Burrows-Wheeler transforms, Journal of Computational Biology, 12(4), pp. 407-415, 2005. R. Gonzalez, S. Grabowski, V. Makinen, and G. Navarro. Practical Implementation of Rank and Select Queries. In Proc. WEA'05, pages 27-38, 2005. Heng Li, Jue Ruan and Richard Durbin, Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 2008 18: 1851-1858.
　ペアをなす２つの配列に対して、指定した間隔以下に近接したマッピング位置を高速に検索する場合は、先ず、各々の配列に対して独立に、インデクスを利用して高速にマッピング位置を求め、次に、それらのマッピング位置の組み合わせの中から、間隔の制約を満たすものを選び出せばよい。しかしながら、配列長が比較的短い場合は、ペアをなす２つの配列に対して、間隔の制約を用いずにそれぞれ独立にマッピング位置を求めると、最初に多数の候補位置が得られ、その後でそれらの中から間隔の制約を満たす少数のものを選び出すことになり、処理速度の低下を招くことがある。一方、サフィックス・アレイなどのインデクスを用いる通常の高速検索方法では、間隔の制約を考慮しながら検索を行うことは困難である。また、ＭＡＱでは、ペアをなす２つの配列のマッピング候補位置が多数あるとき、それらの相対的な位置関係によっては、見落としが生ずる場合もでてくる。
　本発明の課題は、大規模な参照ゲノム配列データと大量のペア・エンドのシーケンス配列データ（クエリー配列データ）が与えられた時、ペア間のマッピング位置の間隔が予め指定された値以下となるようなペア配列のマッピング位置の組み合わせを高速に求めることである。
　ゲノム配列全体に対しては通常のサフィックス・アレイと等価な情報をもち、さらに、種々のスケールで制限されたゲノムの部分領域に対しても、その領域に制限したインデクス情報（局所的なインデクス情報）が得られるような、ゲノム配列データのインデクス化を行う。また、ゲノムの或る領域に対するインデクス情報を、その領域を二分したとき、各々の部分領域に対するインデクス情報に分割する方法を用いる。ペア配列に対しては、先ず、各々の配列に対して独立にゲノム全体領域を検索してインデクス情報を求め、次に、領域の長さが１になるまで二分割を繰り返すことにより、それらのインデクスに対するゲノム上の位置座標を徐々に詳細に決めていく。その過程でペア間の距離の制約を考慮し、制約を満たさないインデクス情報を除去する。
　本発明によるペア文字列検索システムは、参照文字列を記憶する記憶装置と、ペアをなす第１の文字列と第２の文字列及び第１の文字列と第２の文字列間の距離の上限値を入力する入力装置と、参照文字列のサフィックス・アレイを構築するサフィックス・アレイ構築処理部と、第１の文字列と第２の文字列に対して独立に、サフィックス・アレイ内のインデクス領域を計算するインデクス領域検索処理部と、局所化可能なインデクス情報ＬＳＡを構築するＬＳＡ構築処理部と、インデクス領域局所化処理部と、インデクス領域間距離評価処理部と、上限値以内の距離でペアをなす第１の文字列と第２の文字列の参照文字列上の位置を出力する出力装置とを有し、インデクス領域局所化処理部は、計算の過程で、空となるインデクス領域の情報、及び、インデクス領域間距離評価処理部による評価の結果、上限値以下の距離でペアを作れないようなインデクス領域の情報を除去する。
　ＬＳＡ構築処理部は、参照文字列のサフィックス・アレイを初期ブロックとして、その最上位ビットが０か１かに従ってサフィックス・アレイを参照文字列の前半と後半の部分領域に対応する２つのブロックに分割し、その際に参照した最上位ビットを並べたものをＬＳＡの第一列とし、得られた２つのブロックに対して第二位ビットが０か１かに従って更に２つずつのブロックに分割することにより、参照ゲノム配列全体を４分割する部分領域に対応する４つのブロックに分割し、その際に参照した第二位ビットを並べたものをＬＳＡの第二列とし、以下同様にして部分領域の長さが１となるまで処理を繰り返し、その過程でＬＳＡの各列を計算してＬＳＡを構築する。
　インデクス領域局所化処理部は、第１の文字列と第２の文字列の各々に対するサフィックス・アレイ内のインデクス領域を初期ブロック内のインデクス領域として、参照文字列の前半と後半の部分領域に対応する２つのブロックへの分割に対応して、初期ブロック内のインデクス領域をＬＳＡの第一列の情報を用いてそれぞれのブロック内のインデクス領域に分割し、参照文字列全体を４分割する部分領域に対応する４つのブロックへの分割に対応して２つのブロック内のインデクス領域をＬＳＡの第二列の情報を用いて４つのブロック内のインデクス領域に分割し、以下同様な分割処理を繰り返し、初期ブロック内のインデクス領域を、参照文字列の部分領域内のインデクス領域に局所化する計算を行う。インデクス領域局所化処理部は、ＬＳＡの各列のブロックにランク関数を適用して、当該ブロックが分割された２つのブロックでのインデクス領域の位置を計算する。
　インデクス領域間距離評価処理部は、部分領域に局所化されたインデクス領域の参照文字列上の位置座標は部分領域の長さだけの不確定範囲を残すものとして、参照文字列の部分領域に局所化された第１の文字列と第２の文字列のインデクス領域どうしの距離の上限を評価する。
　本発明では、部分領域に対するインデクス情報は、領域分割により位置座標を詳細化する過程で、ペア間の距離制約を満たさないものが除去されていくため、マッピング位置の間隔の制約を考慮しながら高速にインデクス検索（マッピング計算）を行うことが可能になる。
ペア配列のゲノム・マッピング計算を行うシステムの構成例を示す図である。 領域分割に伴いインデクス情報を分割する方法を説明するための図である。 ＬＳＡ（localized suffix array）の説明図である。 インデクス領域どうしの距離の評価方法の説明図である。 インデクス領域の局所化の過程を説明する図である。 ＬＳＡの構築方法を説明するためのフローチャートである。 ＬＳＡの構築方法の説明図である。 ＬＳＡの構築の具体例を説明する図である。 ＬＳＡの構築の具体例をビットレベルで詳細に説明するための図である。 インデクス領域局所化処理の処理手順を説明するフローチャートである。 インデクス領域のリストの更新方法を説明するフローチャートである。 スプライシングを含むゲノム・マッピング計算を行うシステムの構成例を示す図である。 アミノ酸配列断片から蛋白質を同定するシステムの構成例を示す図である。 一般の文書データから二つの文字列が近接して現れるような文書を同定する計算を行うシステムの構成例を示す図である。
　以下、本発明の実施例を、図面を用いて詳細に説明する。
［実施例１］
　本実施例では、参照ゲノム配列データと、クエリー配列のペアが与えられた時、指定された間隔以下に近接したマッピング位置を求めるための方法を説明する。本実施例では、クエリー配列とゲノム配列の間のマッピングは完全一致の場合を考えるが、数塩基のミスマッチを許容するような不完全一致の場合への応用は、容易な類推により可能である。また、本実施例では、ペア配列のマッピング位置の間隔の上限のみの制約を考えるが、同時にペア配列のマッピング位置の間隔の下限の制約を考えることも、容易な類推により可能である。
　図１は、本実施例のペア配列のマッピング計算を行うシステムの構成例を示す図である。１０１は、このシステムが稼働している計算機であり、外部記憶装置１０２に格納された参照ゲノム配列データ１０３を読み込み、またＤＮＡシーケンサ１１６を用いてＤＮＡサンプル１１５を解析して得られたペア配列データ（クエリー配列のペア）１０７を、入力装置１０６を介して読み込み、また、ユーザにより指定されたパラメータ１１４（ペア間距離の上限値Ｄ）を、入力装置１１８を介して読み込む。計算機１０１は、サフィックス・アレイ構築処理部１０５で、参照ゲノム配列データ１０３内の塩基配列を区切り文字＄を介して連結して、そのサフィックス・アレイを計算し主記憶１０４内に記憶する。主記憶１０４内には、同様に、他の処理部の計算結果も記憶され、その内容は必要に応じて各処理部から参照される。ＬＳＡ構築処理部１０８は、サフィックス・アレイを基に、インデクス情報ＬＳＡ（localized suffix array）を構築する。インデクス領域検索処理部１０９は、ペアをなすクエリー配列のそれぞれに対して独立に、サフィックス・アレイを検索して、インデクス領域の情報を得る。インデクス領域局所化処理部１１０では、クエリー配列のインデクス領域の情報と参照ゲノム配列のＬＳＡをもとにして、ペア間のマッピング位置間隔の制約を満たすようなインデクス領域を局所化（細分化）していき、インデクス領域のゲノム上の位置を、初期のゲノム全体に広がった状態から最終的に塩基レベルで詳細に確定した状態まで詳細化していく。
　これにより最終的に、ペア配列のマッピング位置が得られる。ペア間のマッピング位置間隔の制約の成否を調べるために、異なるインデクス領域間のゲノム上の位置の距離が、インデクス領域間距離評価処理部１１７で評価される。出力処理部１１１は、ペア配列のマッピング位置情報（検索結果１１３）を出力装置１１２に出力する。ペア配列のマッピング位置情報（検索結果１１３）は、染色体名とその向き（＋／－何れのストランドか）と出現位置を示す座標の情報を含んでいる。
　参照ゲノム配列からサフィックス・アレイを構築するには、Manber-Myersなどの公知の方法を利用すればよい（非特許文献２）。また、クエリー配列のインデクス情報、即ち、サフィックス・アレイ内においてクエリー配列がサフィックスのプレフィックスとして現れるインデクスの範囲を求めるためには、一般的な二分探索法を用いればよい。
　図２は、領域の分割に際して、インデクス情報を分割する方法を説明するための説明図である。この方法は、後述するようにインデクス領域局所化処理部１１０を構成するために用いられる。ここでは、説明を容易にするため、ゲノムの長さは２のべき乗（２のｈ乗）マイナス１であると仮定する。ゲノムの長さがそれ以外の一般の場合への拡張は、容易な類推により可能である。クエリー配列と同一の配列が、２０１で示すゲノム配列上に、２０２で示すように多数散在して出現しているとする。このとき、ゲノム配列全体の領域に対するサフィックス・アレイＢ（２０３）上では、サフィックス・アレイの性質により、それら散在した多数の出現位置に対応するインデクスが、２０４に示すように一か所にまとまって現れる。ゲノムの長さを２のｈ乗マイナス１と仮定したので、各々のサフィックスの位置は０以上、２のｈ乗未満の整数で表わされる。空文字列は最終サフィックスとする。ゲノムの前半に対応するサフィックスのインデクスは、二進表記したとき最上位ビット（最下位からｈビット目）が０の整数で表わされ、また、ゲノムの後半に対応するサフィックスのインデクスは、二進表記したとき最上位ビット（最下位からｈビット目）が１の整数で表わされる。そこで、ゲノム全体の領域に対するサフィックス・アレイＢを、内部の相対的な順序を保ったまま、最上位ビットが０か１かにより二分すると、ゲノム前半に対応するサフィックス・アレイの部分列Ｂ０と、ゲノム後半に対応するサフィックス・アレイの部分列Ｂ１が得られる。ゲノム配列全体を二分することにより、クエリー配列と同一の配列は、ゲノム配列前半２０５上とゲノム配列後半２０７上に、それぞれ２０６と２０８で示すように多数散在して出現するが、それらのインデクスは、２１０と２１２で示すように、それぞれＢ０とＢ１内に一か所にまとまって現れる。
　ゲノム全体の領域に対するサフィックス・アレイＢにおいて、先頭行を０行目として、クエリー配列に対応するインデクス領域２０４の先頭行をｓ行目、最終行の次の行をｔ行目とする。同様に、ゲノム前半の領域に対するサフィックス・アレイＢ０において、クエリー配列に対応するインデクス領域２１０の先頭行をｓ０行目、最終行の次の行をｔ０行目とし、また、ゲノム後半の領域に対するサフィックス・アレイＢ１において、クエリー配列に対応するインデクス領域２１２の先頭行をｓ１行目、最終行の次の行をｔ１行目とする。ゲノム全体の領域に対するサフィックス・アレイＢの最上位ビット列Ｌにおいて、ｓ行目より前にある０の個数はランク関数rank(Ｌ，ｓ，０)で表わされ、また、ｓ行目より前にある１の個数はランク関数rank(Ｌ，ｓ，１)で表わされ、
　　　　　ｓ０＝rank(Ｌ，ｓ，０)
　　　　　ｓ１＝rank(Ｌ，ｓ，１)
となる。ランク関数は公知の方法で高速に計算することができる（非特許文献５）。同様に、
　　　　　ｔ０＝rank(Ｌ，ｔ，０)
　　　　　ｔ１＝rank(Ｌ，ｔ，１)
となる。従って、インデクスの最上位ビットによるゲノム全体領域の二分割に際して、クエリー配列に対応するインデクスの範囲（ｓ，ｔ）のゲノム前半とゲノム後半に対応する部分への分割（ｓ０，ｔ０），（ｓ１，ｔ１）は高速に計算される。なお、分割されたうちの一方は空になることもある。
　ゲノムの前半と後半のそれぞれに対しても、インデクスの上位から２ビット目（最下位からｈ－１ビット目）の値が０か１かを調べることにより、同様な二分割を行うことができ、それによりゲノム全体は四分割される。続けて、８分割、１６分割なども可能で、最終的には領域の長さが１に達するまで、再帰的に分割を繰り返すことができる。
　このようなゲノム領域の再帰的な分割に伴うインデクス領域の分割の計算に必要な、サフィックスのインデクスのビットを並べたものをＬＳＡ（localized suffix array）とよぶ。図３はｈ＝３のときのＬＳＡを説明するための図である。３０１はＬＳＡ全体を表し、３０２は、最初の二分割でゲノム全体の領域を前半と後半の２つの領域に分割する際に参照したＢ内のインデクスの最上位ビットを並べたものである。３０３は、次の４分割において、すなわち、上の２つの領域の二分割において、Ｂ０、Ｂ１内で参照されたインデクスの上位から２ビット目のビットを並べたものである。同様に、３０４は次の８分割において、４分割後の各領域から参照されたインデクスの上位から３ビット目のビットを上下に並べたものである。ｈ＞３のときは、１６分割以降に対しても、同様に４列目以降が構成される。
　以上のようなゲノム全体の領域の再帰的な分割の過程で得られる部分領域をブロックとよぶことにする。ブロックの中には、クエリー配列に対応するインデクスの範囲を有するものも有しないものもある。ブロック内のインデクス領域に含まれるインデクスのゲノム上の位置はそのブロック内に限定されるため、インデクス領域のゲノム上の位置座標はブロックの長さに対応する不確定さを残して確定する。従って、そのようなインデクス領域どうしのゲノム上での距離（間隔）を、不確定さを残しつつ評価することが可能になる。
　図４は、インデクス領域間距離評価処理部１１７において、異なるブロック内に含まれるインデクス領域どうしのゲノム上の距離の評価方法を説明する図である。後述するように、この評価方法は、インデクス領域局所化処理部１１０を構成する際に用いられる。４０１はブロックを表す。二つのクエリー配列ｐ、ｑに対するインデクス領域４０２，４０３が異なるブロック４０４，４０５に含まれているとする。それらのインデクス領域のゲノム上の位置座標は、４０６，４０７に示すような不確定範囲を残して確定している。また、ブロックの長さは、ゲノムの全長と分割回数に基づき計算される。従って、それらのインデクス領域の間の距離の上限４０９と下限４０８は、ブロック４０４とブロック４０５の間にあるブロックの数を数えることにより計算できる。
　　（距離の下限）＝（間にあるブロックの数）×（ブロック長）＋１
　　（距離の上限）＝（（間にあるブロックの数）＋２）×（ブロック長）－１
　最終の分割でブロック長が１になると、これらの値は一致する。即ち、不確定さを残さずに、インデクスのゲノム上の位置座標が確定する。
　図５は、２つのクエリー配列ｐ、ｑに対するゲノム全体でのインデクス領域が与えられとき、インデクス領域局所化処理部１１０において、ペア配列としてのゲノム配列上の位置座標が、ゲノム領域の分割を再帰的に繰り返すことにより、不確定さを減少させながら、確定していくようすを説明する図である。５０１はゲノム全体の領域（初期ブロック）に対するサフィックス・アレイである。５０２は、ゲノム全体を２分割したときの前半と後半のブロックに対応するサフィックス・アレイの部分列である。５０３は、ゲノム全体を４分割したときの各ブロックに対応するサフィックス・アレイの部分列を上下に並べたものである。同様に、５０４，５０５は、ゲノム全体をそれぞれ８分割、１６分割したときの、各ブロックに対応するサフィックス・アレイの部分列を上下に並べたものである。５０６は、最終的にブロック長が１に達したときのものである。５０１～５０６の間をつなぐ実線の矢印は、クエリー配列ｐに対するインデクス領域が再帰的な局所化処理（ブロックの細分化処理）により詳細化される様子を示す。
　同様に、５０１～５０６の間をつなぐ破線の矢印は、クエリー配列ｑに対するインデクス領域が局所化処理により詳細化される様子を示す。その過程で、５０７や５０８，５０９に示すように、ペア配列間の距離の制約を満たさないものは、処理を打ち切られる。ペア配列間の距離の制約を満たさないものの中には、５０９のようにペア配列の順序が逆になったものも含まれる。最終的に５０６上まで残った５１０や５１１がペア配列のマッピング位置を与える。
　図６に、ＬＳＡの構築処理手順を示す。ステップ６０１では、ブロックＢを、サフィックス・アレイ構築処理部１０５で計算され、主記憶１０４内に記憶されている、参照ゲノム配列全体に対するサフィックス・アレイとする。これまで同様、説明を容易にするため、参照ゲノム配列長は２のｈ乗マイナス１であると仮定する。Ｂはインデクス（ｈビットの整数）を２のｈ乗個並べた配列である。ステップ６０３でｉ＝１として以下の処理を繰り返す。ステップ６０４はその終了判定処理である。ステップ６０５では、０と１からなる長さ（ｉ－１）の文字列ａにより、Ｂａと表わされる各ブロックに対して、ステップ６０６，６０７の処理を行う。ステップ６０６では、ブロックＢａ内に含まれる各インデクス（ｈビットの整数）に対して、それらの上位からｉビット目を並べてできるビット列をＬａとする。ステップ６０７では、それらのビットが０であるか１であるかに従い、Ｂａ内に含まれる各インデクスを、ブロックＢａ０又はＢａ１に分類する。その際、インデクス間の相対順序は保つ。そのような全てのＢａに対して処理が完了したら、ステップ６０８でｉの値を１だけ増やしてステップ６０４の直前に戻る。ステップ６０５で終了条件が成立したら、ステップ６０９でＬａからＬＳＡを主記憶１０４上に構築する。
　図７は、ＬａからＬＳＡを構築する方法の説明図である。ＬＳＡ７０１は、２のｈ乗行、ｈ列にビットを並べたものである。その第１列７０２はＬに等しい。その第２列７０３は、Ｌ０とＬ１をつないだものである。その第３列７０４は、Ｌ００，Ｌ０１，Ｌ１０，Ｌ１１をつないだものである。第４列７０５は、Ｌ０００，Ｌ００１，Ｌ０１０，Ｌ０１１，Ｌ１００，Ｌ１０１，Ｌ１１０，Ｌ１１１をつないだものである。第５列以降も同様である。
　ゲノム配列は実際には塩基を表す文字Ａ，Ｃ，Ｇ，Ｔ，Ｎからなるが、本発明による方法は一般のアルファベットからなる文字列に対して適用可能である。図８は、ｈ＝３で、参照ゲノム配列を“ａｂｒａｃａｄ”としたときにＬＳＡを構築する例を説明する図である。８０１は、“ａｂｒａｃａｄ”の全てのサフィックスとそのインデクスのリスト（空文字列を含む）である。８０２は、それらを辞書式順番にソートしたもので、サフィックス・アレイである。８０３は、サフィックス・アレイのインデクスの値を２進表記したものであり、ゲノム全体に対応するＢである。８０４は、図６で説明した方法でＢから構築したＬＳＡであり、図７に対応する。
　図９は、参照ゲノム配列を“ａｂｒａｃａｄ”としたときに、図６で説明した方法に従って、ゲノム全体に対応するＢからＬＳＡを構築する過程を説明するための説明図である。９０１は、ゲノム全体に対応するＢである。９０２は、ゲノムの前半と後半に対応するＢ０とＢ１を破線で区切って並べたものである。ＢからＢ０とＢ１への分割に際しては、Ｂの最上位ビット（太字）を参照する。これらのビットを並べたものがＬ０であり、９０５で示すＬＳＡの第１列目となっている。同様に、９０３は、ゲノムの４分割に対応するＢ００，Ｂ０１，Ｂ１０，Ｂ１１を破線で区切って並べたものである。Ｂ０とＢ１からＢ００，Ｂ０１，Ｂ１０，Ｂ１１への分割に際しては、Ｂ０とＢ１の上から２番目のビット（太字）を参照する。それらのビットを並べたものがＬ０とＬ１であり、９０５で示すＬＳＡの第２列目となっている。同様に、９０４は、ゲノムの８分割に対応するＢ０００，Ｂ００１，Ｂ０１０，Ｂ０１１，Ｂ１００，Ｂ１０１，Ｂ１１０，Ｂ１１１を破線で区切って並べたものである。Ｂ００，Ｂ０１，Ｂ１０，Ｂ１１からこれら８個への分割に際しては、Ｂ００，Ｂ０１，Ｂ１０，Ｂ１１の上から３番目のビット（太字）を参照する。それらのビットを並べたものがＬ００，Ｌ０１，Ｌ１０，Ｌ１１であり、９０５で示すＬＳＡの第３列目となっている。９０４においてブロック長は１となり、ここで再帰的な分割は終了する。従って、９０５に示す８行３列のビット行列がＬＳＡである。
　図１０は、インデクス領域局所化処理１１０の内部の処理の流れを説明する図である。ステップ１００１でユーザから入力装置１０６を介して距離の上限値Ｄを読み込み、ステップ１００２でゲノム配列全体のサフィックス・アレイＢ内におけるクエリー配列ｐの出現範囲（ｓp，ｔp）とクエリー配列ｑの出現範囲（ｓq，ｔq）を主記憶１０４内から読み込む。Ｂはサフィックス・アレイ構築処理部１０５で既に構築済みであり、（ｓp，ｔp）と（ｓq，ｔq）の値はインデクス領域検索部１０９で計算済みであり、これらは主記憶１０４内に保持されている。ステップ１００３では、クエリー配列ｐのマッピング位置の候補リスト示す情報List(ｐ)と、クエリー配列ｑのマッピング位置の候補リスト示す情報List(ｑ)を、次式に示すように初期化する。
　　　　　List(ｐ)＝｛（Ｂ，ｓp，ｔp）｝，
　　　　　List(ｑ)＝｛（Ｂ，ｓq，ｔq）｝
　これらのリストは、マッピング位置座標を詳細に決めていくために、以下の反復処理により更新される。ステップ１００４は、この反復処理の終了判定処理であり、ブロック長が１になったら終了と判定する。ステップ１００５では、図１１を用いて以下に説明する方法に従い、リストList(ｐ)に属する各要素を細分化してリストList(ｐ)を更新し、同様にList(ｑ)に属する各要素を細分化してリストList(ｑ)を更新する。
　図４を用いて説明したように、各クエリー配列のインデクス領域のマッピング位置は、ブロック長だけの不確定さを残して確定し、それに基づいて、異なるブロックに属するインデクス領域のマッピング位置どうしの距離の上限や下限を評価できる。ステップ１００６では、リストList(ｐ)とList(ｑ)の要素を組み合わせて評価し、ペア間の距離の上限がＤ以下となるようなペアを作れないような、リストList(ｐ)とList(ｑ)の要素を削除する。リストList(ｐ)とList(ｑ)の要素の組み合わせを考える際には、それらの要素をマッピング位置の順にソートしておくと効率的である。反復処理が終了したとき、ステップ１００７においてリストList(ｐ)とList(ｑ)を回答として、主記憶１０４内に返す。
　図１１は、クエリー配列のマッピング位置の候補リスト示す情報List_inを更新して新しいリストList_outを回答するリストの細分化処理１００５の処理手順を説明する図である。この処理手順は、図２を用いて説明した方法と同様な方法である。ステップ１１０１で、リストList_inを入力し、リストList_outを空リストとする。List_inは、ステップ１００５の細分化される前のList(ｐ)又はList(ｑ)である。ステップ１１０２では、リストList_inの各要素（Ｂａ，ｓ，ｔ）に対して以下の処理を行う。ここで、ａは０と１からなる文字列を表し、ｓとｔはクエリー配列のブロックＢａ内のインデクス領域を示す整数である。ステップ１１０３では、次式に従い、細分化されたブロックＢａ０とＢａ１内のインデクス領域を示す整数であるｓ０，ｔ０，ｓ１，ｔ１を計算する。
　　　　　ｓ０＝rank(Ｌａ，ｓ，０)
　　　　　ｔ０＝rank(Ｌａ，ｔ，０)
　　　　　ｓ１＝rank(Ｌａ，ｓ，１)
　　　　　ｔ１＝rank(Ｌａ，ｔ，１)
　ここで、rank(Ｌａ，ｓ，０)はランク関数で、ＬＳＡ内のビット列Ｌａにおいて、ｓ行目の位置より前に現れる０の個数を表し、また、rank(Ｌａ，ｓ，１)もランク関数で、ＬＳＡ内のビット列Ｌａにおいて、ｓ行目の位置より前に現れる１の個数を表す。ランク関数は公知の方法で効率的に計算できる（非特許文献５）。rank(Ｌａ，ｔ，０)とrank(Ｌａ，ｔ，１)についても同様である。ステップ１１０４では、ＬＳＡ内の対応するビット列Ｌａ内の各ビットが０か１かに応じて、（Ｂａ，ｓ，ｔ）を２つのブロック（Ｂａ０，ｓ０，ｔ０）と（Ｂａ１，ｓ１，ｔ１）に分割して、これらをリストList_outに追加する。ステップ１１０５では、得られたリストList_outを回答として、主記憶１０４内に返す。これは、ステップ１００５の細分化された後のList(ｐ)又はList(ｑ)となる。
　本発明では、クエリー配列に対する多数のマッピング位置候補が各ブロック内のインデクス領域（二つの整数ｓとｔ）としてコンパクトに表現される。また、ブロックの分割とそれに伴うインデクス領域の分割を繰り返し行うことにより、インデクス領域のマッピング位置を、不確定範囲を徐々に減少させながら、最終的に塩基レベルまで詳細に決められる。その際、不確定範囲を残したマッピング位置間の距離の上限を評価して、ペア間で指定されたマッピング位置の間隔の制約（距離Ｄ以下）を満たし得ないようなマッピング位置に対するインデクス情報を除去して、無駄な計算を避ける。また、ＬＳＡを参照することにより、ブロック分割に伴うインデクス領域の分割は効率的に計算される。ＬＳＡを効率的に構築する処理手順も与えた。
　以上により、大規模な参照ゲノム配列データと大量のペア・エンドのシーケンス配列データ（クエリー配列データ）が与えられた時、ペア間のマッピング位置の間隔が予め指定された値以下となるようなペア配列のマッピング位置の組み合わせを高速に求めることが可能となる。
［実施例２］
　前記実施例では、ペア・エンド法で得られたクエリー配列のゲノム・マッピングの計算法を示したが、スプライシングを受けたｍＲＮＡから作られたｃＤＮＡ配列の断片配列をシーケンシングして得られたクエリー配列のゲノム・マッピング計算する問題に対しても本発明を実施することができる。
　図１２は、本実施例のスプライシングを含むマッピング計算を行うシステムの構成例を示す図である。１２０１は、このシステムが稼働している計算機であり、外部記憶装置１２０２に格納された参照ゲノム配列データ１２０３を読み込み、また、ＤＮＡシーケンサ１２１７を用いてｃＤＮＡサンプル１２１６を解析して得られたクエリー配列データ１２０７を、入力装置１２０６を介して読み込み、また、ユーザにより指定されたパラメータ１２１９（イントロン長の上限値Ｄ）を、入力装置１２１８を介して読み込む。計算機１２０１は、サフィックス・アレイ構築処理部１２０５で、参照ゲノム配列データ１２０３内の塩基配列を区切り文字＄を介して連結して、そのサフィックス・アレイを計算し主記憶１２０４内に記憶する。主記憶１２０４内には、同様に、他の処理部の計算結果も記憶され、その内容は必要に応じて各処理部から参照される。
　ＬＳＡ構築処理部１２０８は、サフィックス・アレイを基に、新たなインデクス情報ＬＳＡ（localized suffix array）を構築する。ペア配列生成処理部１２０９は、クエリー配列を種々の位置で二分割して種々のペア配列を作成する。インデクス領域検索処理部１２１０は、ペア配列のそれぞれに対して独立に、サフィックス・アレイを検索して、インデクス情報を得る。インデクス領域局所化処理部１２１１では、ペア配列のインデクス情報と参照ゲノム配列のＬＳＡをもとにして、ペア間のマッピング位置間隔がイントロン長上限値以下となる制約を満たすようなインデクスに対応する座標を、ゲノム全体から徐々に局所化しながら決めていく。これにより、ペア配列としてのマッピング位置が得られる。
　コンセンサス配列検査処理部１２１２では、スプライス・サイト周辺のゲノム配列を調べて、コンセンサスを満たさないようなマッピング結果を除去する。出力処理部１２１３は、マッピング位置情報（検索結果１２１５）を出力装置１２１４に出力する。スプライシングを受けた配列のマッピング情報（検索結果１２１５）は、染色体名とその向き（＋／－何れのストランドか）と分断された二つの出現位置を示す座標の情報を含んでいる。
　サフィックス・アレイ構築処理部１２０５、ＬＳＡ構築処理部１２０８、インデクス領域検索処理部１２１０、インデクス領域局所化処理部１２１１等の処理内容は、実施例１と同様である。ペア配列生成処理部１２０９では、分割して得られるペア配列の配列長の下限値を予め指定しておいて、その条件を満たすペア配列を全て生成する。それらはペア配列データとして、インデクス領域検索処理部１２１０、インデクス領域局所化処理部１２１１で処理されるが、コンセンサス配列検査処理部１２１２で、スプライシングを受けたマッピングとして解釈され検査される。スプライス・サイト周辺のゲノム配列のコンセンサスとしては、ペア配列としてのゲノム・マッピング位置にはさまれた部分をイントロン配列と解釈して、イントロン配列がＧＴで始まりＡＧで終わるという条件を課す。この条件を満たさないペア配列のマッピング結果は除去される。
　大規模な参照ゲノム配列データとｃＤＮＡをシーケンスした大量のクエリー配列データが与えられた時、前記実施例で示した、ペア間のマッピング位置の間隔が予め指定された値以下となるようなペア配列のマッピング位置の組み合わせを高速に求める方法を上記のようにして利用することにより、妥当な長さのイントロン配列長（例えば、数百キロ塩基以下）をもつスプライシングを含むゲノム・マッピングを高速に計算することが可能となる。
［実施例３］
　前記実施例では、Ａ，Ｃ，Ｇ，Ｔ，Ｎからなる塩基配列データに対する実施例を説明した。本発明は、塩基配列データとは限らない一般の文字列データに対しても適用可能である。ここでは、アミノ酸断片配列データに対する実施例を説明する。アミノ酸配列は、２０種類程のアミノ酸からなる配列である。
　図１３に、本実施例のアミノ酸断片配列から蛋白質を同定する計算を行うシステムの構成例を示す図である。１３０１は、このシステムが稼働している計算機であり、外部記憶装置１３０２に格納された参照蛋白質配列データ１３０３を読み込み、また、質量分析システム１３０５を用いて蛋白質サンプル１３０６を解析して得られた複数のクエリー配列データ１３０７を、入力装置１３０８を介して読み込み、また、ユーザにより指定されたパラメータ１３１８（蛋白質配列長の上限値Ｄ）を、入力装置１３１９を介して読み込む。計算機１３０１は、サフィックス・アレイ構築処理部１３１０で、参照蛋白質配列データ１３０３内のアミノ酸配列を区切り文字＄を介して連結して、そのサフィックス・アレイを計算し主記憶１３０９内に記憶する。主記憶１３０９内には、同様に、他の処理部の計算結果も記憶され、その内容は必要に応じて各処理部から参照される。
　ＬＳＡ構築処理部１３１１は、サフィックス・アレイを基にインデクス情報ＬＳＡ（localized suffix array）を構築する。インデクス領域検索処理部１３１２は、クエリー配列のそれぞれに対して独立に、サフィックス・アレイを検索して、インデクス情報を得る。インデクス領域局所化処理部１３１３では、クエリー配列のインデクス情報と参照ゲノム配列のＬＳＡをもとにして、参照蛋白質配列データ１３０３内に近接して現れるようなインデクスの出現位置を、参照蛋白質配列データ全体から徐々に局所化しながら決めていく。これにより、クエリー配列の集団としての参照蛋白質配列データ中の出現位置が得られる。蛋白質同定処理部１３１４では、クエリー配列の参照蛋白質配列データ中の出現位置に基づきクエリー配列を含む蛋白質を同定し、複数のクエリー配列を含むような蛋白質を同定する。出力処理部１３１５では、同定された蛋白質の名称とクエリー配列の出現位置（断片位置）の情報を含む検索結果１３１７を出力装置１３１６に出力する。
　サフィックス・アレイ構築処理部１３１０、ＬＳＡ構築処理部１３１１、インデクス領域検索処理部１３１２、インデクス領域局所化処理部１３１３等における処理内容は、前記実施例と同様である。
　大規模な参照蛋白質配列データと蛋白質サンプルを質量分析システムで解析して得られた大量のクエリー配列データ（アミノ酸配列データ）が与えられた時、前記実施例で示した、ペア間のマッピング位置の間隔が予め指定された値以下となるようなペア配列のマッピング位置の組み合わせを高速に求める方法を上記のようにして利用することにより、互いに近接して現れるようなクエリー配列の参照蛋白質配列データ中の出現位置が求まり、さらに、それらのクエリー配列の出現位置を複数含むような信頼性の高い蛋白質を高速に同定することが可能となる。
［実施例４］
　前記実施例では、塩基配列やアミノ酸配列などの生物学的な配列データを扱ってきたが、本発明は、一般の文書データへの適用も可能である。
　図１４に、本実施例の一般の文書データの中から二つの検索文字列が近接して現れるような文書を同定する計算を行うシステムの構成例を示す図である。１４０１は、このシステムが稼働している計算機であり、外部記憶装置１４０２に格納された参照文書データ１４０３を読み込み、また、検索文字列のペア（クエリー）１４０７を、入力装置１４０６を介して読み込み、また、ユーザにより指定されたパラメータ１４１１（ペア間距離の上限値Ｄ）を、入力装置１４１６を介して読み込む。本実施例では、参照文書データ１４０３は、文書名を付けられた文字列本文からなる文書の集まりと仮定するが、階層的な構造をもつ一般の構造化文書に対しても、容易な類推により本発明を適用することができる。
　計算機１４０１は、サフィックス・アレイ構築処理部１４０５で、文書データ１４０３内の文字列を識別可能な区切り文字を介して連結して、そのサフィックス・アレイを計算し主記憶１４０４内に記憶する。主記憶１４０４内には、同様に、他の処理部の計算結果も記憶され、その内容は必要に応じて各処理部から参照される。ＬＳＡ構築処理部１４０８は、サフィックス・アレイを基にインデクス情報ＬＳＡ（localized suffix array）を構築する。インデクス領域検索処理部１４０９は、検索文字列のそれぞれに対して独立に、サフィックス・アレイを検索して、インデクス情報を得る。インデクス領域局所化処理部１４１０では、検索文字列のインデクス情報と文書データのＬＳＡをもとにして、文書データ１４０３内に近接して現れるようなインデクスに対する文書中の出現位置を、文書全体から徐々に局所化しながら決めていく。これにより、検索文字列ペアとしての文書中の出現位置が得られる。文書同定処理部１４１２では、二つの検索文字列が文書データ内の同一の文書に出現しているかを検査する。出力処理部１４１３では、二つの検索文字列が近接して出現するような文書の文書名と各々の検索文字列の出現位置の情報を含む検索結果１４１５を出力装置１４１４に出力する。
　サフィックス・アレイ構築処理部１４０５、ＬＳＡ構築処理部１４０８、インデクス領域検索部１４０９、インデクス領域局所化処理部１４１０等の処理内容は、前記実施例と同様である。
　大規模な文書データと大量の検索文字列ペアが与えられた時、前記実施例で示した、ペア間の出現位置の間隔が予め指定された値以下となるような文字列ペアの出現位置の組み合わせを高速に求める方法を上記のようにして利用することにより、互いに近接して現れるような検索文字列の文書データ中の出現位置が求まり、それらが出現する文書を高速に同定することが可能となる。
２０１　参照ゲノム配列
２０２　クエリー配列と同一の配列が参照ゲノム配列上に出現する位置
２０４　参照ゲノム全体に対するサフィックス・アレイ内で、クエリー配列に対応するインデクス領域
２０５　参照ゲノムの前半の配列
２０６　クエリー配列と同一の配列が参照ゲノムの前半の配列上に出現する位置
２０７　参照ゲノムの後半の配列
２０８　クエリー配列と同一の配列が参照ゲノムの後半の配列上に出現する位置
２１０　参照ゲノムの前半に対するサフィックス・アレイ内で、クエリー配列に対応するインデクス領域
２１２　参照ゲノムの後半に対するサフィックス・アレイ内で、クエリー配列に対応するインデクス領域

Claims (6)

  1. 　参照文字列を記憶する記憶装置と、
    　ペアをなす第１の文字列と第２の文字列及び前記第１の文字列と第２の文字列間の距離の上限値を入力する入力装置と、
    　前記参照文字列のサフィックス・アレイを構築するサフィックス・アレイ構築処理部と、
    　前記第１の文字列と第２の文字列に対して独立に、前記サフィックス・アレイ内のインデクス領域を計算するインデクス領域検索処理部と、
    　前記サフィックス・アレイを初期ブロックとして、その最上位ビットが０か１かに従って前記サフィックス・アレイを前記参照文字列の前半と後半の部分領域に対応する２つのブロックに分割し、その際に参照した前記最上位ビットを並べたものをＬＳＡの第一列とし、得られた前記２つのブロックに対して第二位ビットが０か１かに従って更に２つずつのブロックに分割することにより、前記参照ゲノム配列全体を４分割する部分領域に対応する４つのブロックに分割し、その際に参照した前記第二位ビットを並べたものをＬＳＡの第二列とし、以下同様にして前記部分領域の長さが１となるまで処理を繰り返し、その過程でＬＳＡの各列を計算して、局所化可能なインデクス情報ＬＳＡを構築するＬＳＡ構築処理部と、
    　前記第１の文字列と第２の文字列の各々に対する前記サフィックス・アレイ内のインデクス領域を前記初期ブロック内のインデクス領域として、前記参照文字列の前半と後半の部分領域に対応する２つのブロックへの分割に対応して、前記初期ブロック内のインデクス領域を前記ＬＳＡの第一列の情報を用いてそれぞれのブロック内のインデクス領域に分割し、前記参照文字列全体を４分割する部分領域に対応する４つのブロックへの分割に対応して前記２つのブロック内のインデクス領域を前記ＬＳＡの第二列の情報を用いて前記４つのブロック内のインデクス領域に分割し、以下同様な分割処理を繰り返し、前記初期ブロック内のインデクス領域を、前記参照文字列の部分領域内のインデクス領域に局所化する計算を行うインデクス領域局所化処理部と、
    　前記部分領域に局所化された前記インデクス領域の前記参照文字列上の位置座標は前記部分領域の長さだけの不確定範囲を残すものとして、前記参照文字列の部分領域に局所化された前記第１の文字列と第２の文字列の前記インデクス領域どうしの距離の上限を評価するインデクス領域間距離評価処理部と、
    　前記上限値以内の距離でペアをなす前記第１の文字列と第２の文字列の前記参照文字列上の位置を出力する出力装置とを有し、
    　前記インデクス領域局所化処理部は、前記計算の過程で、空となるインデクス領域の情報、及び、前記インデクス領域間距離評価処理部による評価の結果、前記上限値以下の距離でペアを作れないようなインデクス領域の情報を除去することを特徴とするペア文字列検索システム。
  2. 　請求項１に記載のペア文字列検索システムにおいて、前記インデクス領域局所化処理部は、前記ＬＳＡの各列のブロックにランク関数を適用して、当該ブロックが分割された２つのブロックでのインデクス領域の位置を計算することを特徴とするペア文字列検索システム。
  3. 　請求項１又は２に記載のペア文字列検索システムにおいて、
    　前記参照文字列は参照ゲノム配列データであり、
    　前記ペアをなす第１の文字列と第２の文字列はＤＮＡシーケンサから得られたＤＮＡの第１の塩基配列と第２の塩基配列であることを特徴とするペア文字列検索システム。
  4. 　請求項１又は２に記載のペア文字列検索システムにおいて、
    　前記参照文字列は参照ゲノム配列データであり、
    　前記ペアをなす第１の文字列と第２の文字列は、ＤＮＡシーケンサから得られたｃＤＮＡの塩基配列を２分割して生成された第１の塩基配列と第２の塩基配列のペアであり、
    　前記上限値はイントロン長の上限値であり、
    　前記参照ゲノム配列データ内のスプライス・サイト周辺でコンセンサス配列の検査を行うコンセンサス検査処理部を有し、
    　前記コンセンサス検査処理部は、前記第１の塩基配列と第２の塩基配列に挟まれた部分をイントロン配列と解釈し、スプライス・サイト周辺にコンセンサス配列が現れないような検索結果を除去することを特徴とするペア文字列検索システム。
  5. 　請求項１又は２に記載のペア文字列検索システムにおいて、
    　蛋白質同定処理部を有し、
    　前記参照文字列は複数の参照蛋白質配列データであり、
    　前記ペアをなす第１の文字列と第２の文字列は、質量分析システムから得られた断片化された蛋白質の２つのアミノ酸配列であり、
    　前記サフィックス・アレイ構築処理部は、前記複数の参照蛋白質配列データ内のアミノ酸配列を区切り文字を介して連結してそのサフィックス・アレイを計算し、
    　前記蛋白質同定処理部は前記ペアをなす２つのアミノ酸配列を含む蛋白質を同定し、
    　前記出力装置は、同定された蛋白質の名称と前記ペアをなす２つのアミノ酸配列の出現位置の情報を出力することを特徴とするペア文字列検索システム。
  6. 　請求項１又は２に記載のペア文字列検索システムにおいて、
    　文書同定処理部を有し、
    　前記参照文字列は複数の文書データを識別可能な区切り文字を介して連結した参照文書データであり、
    　前記文書同定処理部は、前記ペアをなす第１の文字列と第２の文字列が同一の文書に出現しているかを検査し、
    　前記出力装置は、前記ペアをなす第１の文字列と第２の文字列が含まれる文書の文書名と前記第１の文字列及び第２の文字列の出現位置の情報を出力することを特徴とするペア文字列検索システム。
PCT/JP2010/056146 2009-04-13 2010-04-05 ペア文字列検索システム WO2010119783A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US13/264,037 US8788522B2 (en) 2009-04-13 2010-04-05 Pair character string retrieval system
JP2011509261A JP5279897B2 (ja) 2009-04-13 2010-04-05 ペア文字列検索システム

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2009-097270 2009-04-13
JP2009097270 2009-04-13

Publications (1)

Publication Number Publication Date
WO2010119783A1 true WO2010119783A1 (ja) 2010-10-21

Family

ID=42982445

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2010/056146 WO2010119783A1 (ja) 2009-04-13 2010-04-05 ペア文字列検索システム

Country Status (3)

Country Link
US (1) US8788522B2 (ja)
JP (1) JP5279897B2 (ja)
WO (1) WO2010119783A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012203456A (ja) * 2011-03-23 2012-10-22 Hitachi Ltd 文書検索システム、文書検索方法、及びプログラム
WO2014132497A1 (ja) * 2013-02-28 2014-09-04 株式会社日立ハイテクノロジーズ データ解析装置、及びその方法
JP2015219856A (ja) * 2014-05-21 2015-12-07 株式会社日立製作所 解析装置、データベース作成方法、およびシステム
CN108920483A (zh) * 2018-04-28 2018-11-30 南京搜文信息技术有限公司 基于后缀数组的字符串快速匹配方法

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9201916B2 (en) * 2012-06-13 2015-12-01 Infosys Limited Method, system, and computer-readable medium for providing a scalable bio-informatics sequence search on cloud
WO2014145503A2 (en) * 2013-03-15 2014-09-18 Lieber Institute For Brain Development Sequence alignment using divide and conquer maximum oligonucleotide mapping (dcmom), apparatus, system and method related thereto
US11150871B2 (en) * 2017-08-18 2021-10-19 Colossio, Inc. Information density of documents
CN114943021B (zh) * 2022-07-20 2022-11-08 之江实验室 一种tb级增量数据筛选方法和装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000259646A (ja) * 1999-03-05 2000-09-22 Ricoh Co Ltd 情報索引装置
JP2001243238A (ja) * 2000-02-29 2001-09-07 Ricoh Co Ltd 索引作成装置および文字列照合装置、索引作成方法および文字列照合方法ならびに記憶媒体

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000259646A (ja) * 1999-03-05 2000-09-22 Ricoh Co Ltd 情報索引装置
JP2001243238A (ja) * 2000-02-29 2001-09-07 Ricoh Co Ltd 索引作成装置および文字列照合装置、索引作成方法および文字列照合方法ならびに記憶媒体

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HENG LI ET AL.: "Mapping short DNA sequencing reads and calling variants using mapping quality scores", GENOME RESEARCH, vol. 18, November 2008 (2008-11-01), pages 1851 - 1858 *
KOICHI KIMURA: "Cho Heiretsu DNA Sequencer Muke Kosoku?Koseido Genome Mapping Shuho", SHOROKU CD SEIKAGAKU, 4 December 2008 (2008-12-04) *
KOUICHI KIMURA ET AL.: "LOCALIZED SUFFIX ARRAY AND ITS APPLICATION TO GENOME MAPPING PROBLEMS FOR PAIRED-END SHORT READS", GENOME INFORMATICS 2009, vol. 23, 16 December 2009 (2009-12-16), pages 60 - 71 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012203456A (ja) * 2011-03-23 2012-10-22 Hitachi Ltd 文書検索システム、文書検索方法、及びプログラム
WO2014132497A1 (ja) * 2013-02-28 2014-09-04 株式会社日立ハイテクノロジーズ データ解析装置、及びその方法
JP5985040B2 (ja) * 2013-02-28 2016-09-06 株式会社日立ハイテクノロジーズ データ解析装置、及びその方法
JP2015219856A (ja) * 2014-05-21 2015-12-07 株式会社日立製作所 解析装置、データベース作成方法、およびシステム
CN108920483A (zh) * 2018-04-28 2018-11-30 南京搜文信息技术有限公司 基于后缀数组的字符串快速匹配方法

Also Published As

Publication number Publication date
JP5279897B2 (ja) 2013-09-04
JPWO2010119783A1 (ja) 2012-10-22
US20120041977A1 (en) 2012-02-16
US8788522B2 (en) 2014-07-22

Similar Documents

Publication Publication Date Title
JP5279897B2 (ja) ペア文字列検索システム
Baeza-Yates and G. Navarro Faster approximate string matching
Lam et al. Compressed indexing and local alignment of DNA
Eskin et al. Finding composite regulatory patterns in DNA sequences
WO2017120128A1 (en) Systems and methods for adaptive local alignment for graph genomes
CA2424031C (en) System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map
JP5183155B2 (ja) 大量配列の一括検索方法及び検索システム
US20050227278A1 (en) Recursive categorical sequence assembly
WO2006130947A1 (en) A method of syntactic pattern recognition of sequences
US20060218138A1 (en) System and method for improving search relevance
WO2015123269A1 (en) System and methods for analyzing sequence data
US7467047B2 (en) Method of discovering patterns in symbol sequences
EP2963575B1 (en) Data analysis device and method therefor
Kumar et al. Fast and memory efficient approach for mapping NGS reads to a reference genome
Sagot et al. Identifying satellites and periodic repetitions in biological sequences
US7085651B2 (en) Method and device for assembling nucleic acid base sequences
Goad Computational analysis of genetic sequences
CN115662523A (zh) 面向群体基因组索引表示与构建的方法及设备
KR100538451B1 (ko) 분산 컴퓨팅 환경에서의 유전자 및 단백질 유사서열 검색시스템 및 그 방법
He Using suffix tree to discover complex repetitive patterns in DNA sequences
Šrámek et al. On-line Viterbi algorithm for analysis of long biological sequences
JP2003256433A (ja) 遺伝子構造解析方法およびその装置
Lecroq et al. Sequence indexing
CN115168661B (zh) 原生图数据处理方法、装置、设备及存储介质
CN111916153B (zh) 一种并行多重序列比对方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 10764368

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2011509261

Country of ref document: JP

WWE Wipo information: entry into national phase

Ref document number: 13264037

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 10764368

Country of ref document: EP

Kind code of ref document: A1