[go: up one dir, main page]

WO2024185097A1 - 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム - Google Patents

近似文字列照合方法及び該方法を実現するためのコンピュータプログラム Download PDF

Info

Publication number
WO2024185097A1
WO2024185097A1 PCT/JP2023/008899 JP2023008899W WO2024185097A1 WO 2024185097 A1 WO2024185097 A1 WO 2024185097A1 JP 2023008899 W JP2023008899 W JP 2023008899W WO 2024185097 A1 WO2024185097 A1 WO 2024185097A1
Authority
WO
WIPO (PCT)
Prior art keywords
key
string
predetermined
computer program
occurrences
Prior art date
Application number
PCT/JP2023/008899
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 PCT/JP2023/008899 priority Critical patent/WO2024185097A1/ja
Publication of WO2024185097A1 publication Critical patent/WO2024185097A1/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/30—Information retrieval; Database structures therefor; File system structures therefor of unstructured textual data
    • G06F16/31—Indexing; Data structures therefor; Storage structures

Definitions

  • the present invention relates to approximate string matching technology, and in particular to an approximate string matching device and method that can be used to analyze the human genome, as well as a computer program for implementing the method.
  • the human genome is the set of genetic information that a person possesses, and the substance that carries this information is DNA (deoxyribonucleic acid), which is made up of approximately 3 billion pairs of bases.
  • the bases are adenine (A), guanine (G), cytosine (C), and thymine (T).
  • A adenine
  • G guanine
  • C cytosine
  • T thymine
  • a person's genetic information is determined by the arrangement (sequence) of these bases.
  • a device called a sequencer is used to read the human genome.
  • the sequencer reads a sample human genome, chops it up into a specified upper limit (about several hundred base pairs), amplifies it, and outputs it as a huge data sequence made up of data fragments.
  • Current sequencers can read one person's human genome in about an hour.
  • the individual pieces of data read by the sequencer are compared with the human genome reference sequence, which has been established as the standard genome sequence for humans, and are then reconstructed into a human genome sequence of the original length for analysis. For example, the position of each piece of data is examined (mapped) in comparison with the human genome reference sequence, and analysis is also performed to determine what mutations are present. Since the data read from the sequencer normally contains errors, statistical processing is performed using, for example, 30 times the amount of redundant data per sample (the equivalent of one person's human genome sequence) to reduce the errors. Therefore, analysis of the human genome sequence requires an enormous amount of calculation, and so typically high-performance computers such as supercomputers, cluster computers, and FPGA (Field-Programmable Gate Array)-based computers are used.
  • high-performance computers such as supercomputers, cluster computers, and FPGA (Field-Programmable Gate Array)-based computers are used.
  • a program tool such as BWA (Burrow-Wheeler Aligner) is used to map data fragments.
  • BWA is composed of three algorithms, BWA-backtrack, BWA-SW, and BWA-MEM.
  • BWA-MEM is widely used as a high-speed algorithm that supports Indels (insertion deletions).
  • BWA-MEM is an algorithm that creates an index using a suffix sequence for the portion of a query string based on the read data fragment that appears repeatedly in the human genome reference sequence, and performs mapping (Non-Patent Document 1).
  • the data fragments are then compared to partial sequences of the human genome reference sequence obtained by mapping, and any mutations present are analyzed.
  • an alignment algorithm is used to compare the partial sequences of the human genome reference sequence with the data fragments.
  • approximate character strings are identified while calculating the degree of mutation (similarity) for each element of the sequence, known as an alignment table, according to dynamic programming (Non-Patent Document 2).
  • BWA-MEM is widely used as an algorithm that allows for gaps and high-speed mapping, but the size of the index created is very large, which causes the problem of requiring a large amount of memory resources.
  • memory needs to be accessed a number of times according to the number of characters (i.e., the number of bases) in the query string, which causes frequent random access to memory, resulting in the problem that the processing speed of the processor is limited by the time it takes to access the memory.
  • BWA-MEM in order to identify the position where a substring (sometimes called a "key") in a query string appears in the human genome reference sequence, random access to the main memory is performed the same number of times as the total number of characters in the query string. This is because the amount of data used in the analysis of the human genome is so large that the necessary data sequence cannot be accommodated all at once in a cache memory that allows faster access.
  • the latency of one random access to the main memory (the time from when an access occurs until the data is actually obtained) is about 1 ⁇ s (0.000001 sec) per bank in current computers.
  • random access to the main memory alone takes about 90,000 seconds (about 25 hours). Therefore, even if a computer equipped with a high-performance processor is used, there is a limit to how much the mapping time can be reduced due to the constraint of memory access time.
  • mapping is completed within the time it takes for the sequencer to read the genome, so a balance between the two is generally maintained.
  • next-generation sequencers are becoming cheaper and are said to have the potential to further shorten read times, which makes it desirable to shorten analysis times accordingly.
  • One approach to shorten analysis times would be to further improve the performance of computers, but the high cost of improving computer performance makes it difficult to put this to practical use.
  • the present invention aims to provide a new technology that can quickly and/or efficiently analyze given query data using reference data.
  • one object of the present invention is to provide an approximate string matching device and an approximate string matching method using the same that can quickly and/or efficiently perform approximate string matching between a given query string and a reference string.
  • Another object of the present invention is to provide an approximate string matching device and an approximate string matching method using the same that can quickly and/or efficiently perform analysis using a human genome reference sequence based on data fragments read by a sequencer.
  • Another object of the present invention is to provide a technique for creating a hierarchical index based on reference strings that are suitable for approximate string matching.
  • the present invention which aims to solve the above problems, comprises the following invention-specific matters and technical features.
  • the present invention is a computer program for causing a computing device to realize a method for searching for an approximate string in a reference string based on a query string.
  • the method includes: creating a hierarchical index based on the reference string; mapping the query string to the reference string by referring to the hierarchical index to identify a substring in the reference string that matches at least a part of the query string; and deriving a string that approximates the matched string as the approximate string based on a matched string based on the reference string and a matched string based on the query string, which are created based on the substring identified by the mapping.
  • creating the hierarchical index includes extracting each first key of a predetermined length from the reference string; creating a first key array in which a hash value calculated based on the first key by a predetermined hash function is assigned to each of the extracted first keys; updating the created first key array; and outputting the updated first key array as the hierarchical index.
  • updating the first key arrangement includes, for each of the first keys in the first key arrangement, identifying the number of occurrences of the first key in the reference string, and updating the first key arrangement by adding a first additional key consisting of at least one character following the first key in the reference string to the first key according to the number of occurrences of the identified first key.
  • Creating the first key array may include sorting each of the first keys in the first key array according to the hash value.
  • updating the first key sequence may include determining whether the identified number of occurrences exceeds a predetermined tolerance; if it is determined that the identified number of occurrences exceeds the predetermined tolerance, creating a new first key by adding to the first key the first additional key consisting of at least one character following the first key in the reference string; and identifying the number of occurrences of the new first key in the reference string for the new first key.
  • updating the first key array may further include sorting the new first keys in the first key array according to the first additional key.
  • updating the first key array may include creating a new first key by sequentially adding a new first additional key to the current first key until it is determined that the identified number of occurrences does not exceed a predetermined tolerance.
  • outputting the first key arrangement as the hierarchical index may include outputting the current first key arrangement as the hierarchical index if it is determined that the identified number of occurrences does not exceed a predetermined allowable value.
  • Performing the mapping may include extracting each second key of a predetermined length from the query string, creating a second key array in which each second key extracted from the query string is assigned a hash value calculated based on the second key using the predetermined hash function, and, for each second key, referencing the hierarchical index at a predetermined sampling interval according to the hash value to identify the start position of appearance and the number of appearances of the second key.
  • Identifying the start position and the number of occurrences of the second key may include determining whether the number of occurrences of the second key exceeds the predetermined tolerance; if it is determined that the number of occurrences of the second key exceeds the predetermined tolerance, creating a new second key by adding a second additional key consisting of at least one character following the second key in the query string to the second key; if it is determined that the number of occurrences of the second key does not exceed the predetermined tolerance, outputting the identified current second key as a matching string and outputting the start position of the occurrences of the second key.
  • identifying the start position of occurrence and the number of occurrences of the second key may further include sequentially adding the new second additional key to the current second key to create the new second key until it is determined that the identified number of occurrences of the second key does not exceed the predetermined tolerance.
  • the predetermined sampling interval of the second key may be changed to be larger.
  • the deriving of the approximate string may include receiving a string pair consisting of the matched string and the match string based on the substring identified by the mapping, performing a predetermined alignment process to derive at least one approximate string based on the string pair, and outputting the derived at least one approximate string.
  • Performing the predetermined alignment process may include creating a predetermined alignment table based on the matched string and the match string, setting a calculation area having a width m centered on a diagonal element of the alignment table, calculating the degree of variation for each element in the set calculation area, determining a maximum degree of variation based on the calculated degree of variation, and deriving the at least one approximate string based on the determined maximum degree of variation.
  • executing the predetermined alignment process may include: comparing the maximum mutation degree with a predetermined lower limit value to determine whether the maximum mutation degree exceeds the predetermined lower limit value; widening the width m of the calculation area to set a new calculation area when it is determined that the maximum mutation degree does not exceed the predetermined lower limit value; and deriving the at least one approximate character string based on an element having the maximum mutation degree when it is determined that the maximum mutation degree exceeds the predetermined lower limit value.
  • the method may be configured to repeat the process of setting a new calculation area by widening the calculation area and calculating the mutation degree until it is determined that the maximum mutation degree exceeds the lower limit value.
  • the specified lower limit value may be the value of the degree of variation when it is assumed that there are m consecutive gaps in a specified element string, and the remaining parts are a perfect or substantial match.
  • the method may further include creating the matched string by appending a corresponding predetermined character string in the reference string to the substring, and creating the match string by appending a corresponding predetermined character string in the query string to the substring.
  • the invention is also a computer program product for causing a computing device to execute a method for creating a hierarchical index for searching for a reference string based on a query string.
  • the method includes extracting each first key of a predetermined length from the reference string, creating a first key array to which each extracted first key is assigned a hash value calculated based on the first key using a predetermined hash function, updating the created first key array, and outputting the updated first key array as the hierarchical index.
  • updating the first key arrangement includes: for each of the first keys in the first key arrangement, identifying a start position of appearance and a number of appearances of the first key in the reference string; and updating the first key arrangement by adding, to the first key, a first additional key consisting of at least one character following the first key in the reference string according to the start position of appearance and the number of appearances of the identified first key.
  • the present invention is a computer program for causing a computing device to execute a method for mapping a query string to a reference string.
  • the method includes: reading a hierarchical index consisting of each first key extracted from the reference string and a first key array in which a hash value calculated based on the first key by a predetermined hash function is assigned to each of the first keys, the first key array having an additional key added to the first key according to the appearance start position and the number of appearances of the first key in the reference string; extracting each second key having a predetermined key length from the query string; creating a second key array in which each of the second keys extracted from the query string is assigned a hash value calculated based on the second key by a predetermined hash function; and calculating the hash value for each of the second keys.
  • the method includes: referring to the hierarchical index at a predetermined sampling interval according to the above, identifying the appearance start position and the number of appearances of the second key by comparing with the first key in the first key array; determining whether the identified number of appearances exceeds a predetermined threshold; if it is determined that the identified number of appearances exceeds a predetermined allowable value, creating a new second key by adding an additional key consisting of at least one or more characters following the second key in the query string to the second key; and if it is determined that the identified number of appearances does not exceed the predetermined threshold, outputting the appearance start position and the identified current second key.
  • identifying the appearance start position and the number of appearances of the second key includes sequentially adding the new additional key to the current second key until it is determined that the identified number of appearances does not exceed the predetermined threshold to create the new second key.
  • the predetermined sampling interval of the second key may be changed to be larger.
  • the present invention is a computer program for causing a computing device to realize a method for identifying mutations between a substring in a reference string and a query string by a predetermined alignment process.
  • the method includes receiving a string pair consisting of a matched string based on the reference string and a match string based on the query string, which are created based on the substring identified by mapping, performing a predetermined alignment process to derive at least one approximate string that is approximate to the matched string based on the string pair, and outputting the derived at least one approximate string.
  • the performing of the predetermined alignment process includes creating a predetermined alignment table based on the matched string and the match string, setting a calculation area having a width m centered on an element on the diagonal of the alignment table, calculating a degree of mutation for each element in the set calculation area, determining a maximum degree of mutation based on the calculated degree of mutation, and deriving the at least one approximate string based on the determined maximum degree of mutation.
  • executing the predetermined alignment process may include: comparing the maximum mutation degree with a predetermined lower limit value to determine whether the maximum mutation degree exceeds the predetermined lower limit value; widening the width m of the calculation area to set a new calculation area when it is determined that the maximum mutation degree does not exceed the predetermined lower limit value; and deriving the at least one approximate character string based on an element having the maximum mutation degree when it is determined that the maximum mutation degree exceeds the predetermined lower limit value. Then, the process of setting a new calculation area by widening the calculation area and calculating the mutation degree may be repeated until it is determined that the maximum mutation degree exceeds the lower limit value.
  • performing the predetermined alignment process may further include setting the value of the degree of variation when it is assumed that there are m consecutive gaps in the predetermined element string and that the remaining portions are completely or substantially identical as the predetermined lower limit value.
  • the present invention also applies to a computer-readable recording medium that non-temporarily records the computer program.
  • the present invention also applies to an apparatus that is made up of hardware and/or firmware configured to execute the method.
  • An approximate string matching device is one aspect of the present invention.
  • a means or unit does not simply mean a physical means, but also includes cases where the functions of that means or unit are realized by software. Furthermore, the functions of one means or unit may be realized by two or more physical means, and the functions of two or more means or units may be realized by one physical means. Furthermore, a "system” refers to a logical collection of multiple devices (or functional modules that realize a specific function), and it does not matter whether each device or functional module is contained within a single housing.
  • the present invention allows fast and/or efficient analysis of given query data using reference data.
  • the present invention allows fast and/or efficient approximate string matching between given query strings and reference strings.
  • the present invention makes it possible to perform fast and/or efficient analysis using a human genome reference sequence based on the data fragments read by the sequencer.
  • the present invention provides a hierarchical index based on reference strings suitable for approximate string matching.
  • FIG. 1 is a block diagram showing an example of a schematic configuration of a computer system according to an embodiment of the present invention.
  • FIG. 2 is a flowchart outlining an example of approximate string matching processing by a computer system according to an embodiment of the present invention.
  • FIG. 3 is a flowchart illustrating an example of an index creation process performed by a computer system according to an embodiment of the present invention.
  • FIG. 4 is a diagram showing an example of a reference string used in a computer system according to an embodiment of the present invention.
  • FIG. 4A is a diagram for explaining an example of a data array structure in the process of creating a hierarchical index based on the reference character string shown in FIG. FIG.
  • FIG. 4B is a diagram for explaining an example of a data array structure in the process of creating a hierarchical index based on the reference character string shown in FIG.
  • FIG. 4C is a diagram for explaining an example of a data array structure in the process of creating a hierarchical index based on the reference character string shown in FIG.
  • FIG. 5 is a diagram showing an example of a table structure indicating key start positions and occurrence counts in the process of creating a hierarchical index based on the reference character string shown in FIG.
  • FIG. 6 is a diagram showing an example of a reference string used in a computer system according to an embodiment of the present invention.
  • FIG. 7 is a flowchart illustrating an example of a process for searching for a reference string based on a query string by a computer system according to an embodiment of the present invention.
  • FIG. 8A is a diagram showing an example of an alignment table for explaining the Needleman-Wunsch algorithm.
  • FIG. 8B is a diagram showing an example of an alignment table for explaining the Needleman-Wunsch algorithm.
  • FIG. 8C is a diagram showing an example of an alignment table for explaining the Needleman-Wunsch algorithm.
  • FIG. 8D is a diagram showing an example of an alignment table for explaining the Needleman-Wunsch algorithm.
  • FIG. 9 is a flowchart illustrating an example of alignment processing using dynamic programming by a computer system according to an embodiment of the present invention.
  • FIG. 9 is a flowchart illustrating an example of alignment processing using dynamic programming by a computer system according to an embodiment of the present invention.
  • FIG. 10A is a diagram showing an example of an alignment table for explaining alignment using dynamic programming according to one embodiment of the present invention.
  • FIG. 10B is a diagram showing an example of an alignment table for explaining alignment using dynamic programming according to one embodiment of the present invention.
  • FIG. 10C is a diagram showing an example of an alignment table for explaining alignment using dynamic programming according to one embodiment of the present invention.
  • FIG. 11A is a diagram showing an example of an approximate character string obtained by alignment using dynamic programming according to an embodiment of the present invention.
  • FIG. 11B is a diagram showing an example of an approximate character string obtained by alignment using dynamic programming according to an embodiment of the present invention.
  • FIG. 12 is a diagram showing an example of the hardware configuration of a computing device in a system according to an embodiment of the present invention.
  • the computer system 1 includes, for example, a host computer 10, a plurality of lower computers 20(1)-20(n), and a database 30.
  • the host computer 10 and the lower computers 20(1)-20(n) are connected to each other so as to be able to communicate with each other via a predetermined interface.
  • the computer system 1 is configured as a centralized computer system in which the host computer 10 comprehensively controls the plurality of lower computers 20(1)-20(n), but is not limited to this and may be configured as, for example, a distributed computer system.
  • each computer can operate cooperatively by parallel distributed processing, but for a specific process, only one representative computer may execute the process.
  • the hardware configuration of the host computer 10 and the lower computer 20 is exemplified in FIG. 12, but detailed description thereof will be omitted since it is obvious to those skilled in the art.
  • the lower computers 20 process tasks in parallel under the control of the host computer 10. For example, each of the multiple lower-level computers 20 performs an analysis process based on a given individual query string in comparison with a reference string.
  • the lower-level computers 20(1) to 20(n) may be referred to simply as "lower-level computers 20" unless there is a need to particularly distinguish between them.
  • the database 30 stores various data, such as reference character strings, under the control of the host computer 10.
  • the database 30 stores a human genome reference sequence.
  • the human genome reference sequence is data showing a base pair sequence defined as a standard human genome sequence.
  • the database 30 also stores an index created based on the reference character string.
  • the index is a type of data array structure used to search for the reference character string. Note that in the figure, the database 30 is configured to be connected so as to be accessible only to the host computer 10, but this is not limiting, and the database 30 may also be configured to be connected so as to be accessible to the lower computer 20 as well.
  • FIG. 2 is a flow chart outlining an example of approximate string matching processing by a computer system according to an embodiment of the present invention. This processing is realized, for example, by a host computer 10 and a number of lower level computers 20, which execute an approximate string matching program under the control of a processor, in cooperation with other hardware components.
  • the host computer 10 reads out the reference string stored in the database 30, and creates an index based on the read reference string (S201). Specifically, the host computer 10 creates an index having a predetermined data array structure by expanding and/or sorting and arranging substrings in the reference string according to predetermined conditions.
  • the data array structure is a hierarchical tree structure. In this disclosure, such an index is referred to as a hierarchical index.
  • the host computer 10 stores the created index based on the reference string in the database 30. Details of the process of creating a hierarchical index based on a reference string will be described later. Note that if a hierarchical index based on a reference string has already been created and stored in the database 30, process S201 may be omitted.
  • the host computer 10 receives a query string from an external device (not shown) and performs mapping based on the received query string by referring to the hierarchical index (S202).
  • the external device is a sequencer that reads out the human genome
  • the query string is a data fragment of, for example, about 150 to 200 base pairs output from the sequencer.
  • the lower computer 20 also performs mapping based on the assigned query string by referring to the hierarchical index under the control of the higher computer 10.
  • Mapping is a process that identifies substrings (matching strings) in a reference string that match at least a portion of a query string. In other words, mapping identifies the start position and length (number of characters) of a substring in the reference string that matches at least a portion of the query string.
  • a hierarchical index is searched or explored for each key in the query string, thereby identifying each key of the query string and its start position of appearance in the reference string. Note that, to speed up processing, the created hierarchical index can be expanded in a high-speed memory such as main memory or cache memory.
  • the sampling interval of the start position of the key appearance in the query string is adjusted to be somewhat larger than that of the conventional technology.
  • the interval of the start position of the key appearance is set to 3 to 5 characters, but in the search according to the present disclosure, it can be set to, for example, 10 characters or more.
  • the sampling interval of the start position of the key appearance is small during a search, the number of searches increases and the efficiency decreases, whereas if the sampling interval of the start position of the key appearance is always large, the probability of overlooking increases in exchange for faster search. Therefore, in the present disclosure, the sampling interval of the start position of the key appearance is increased while the length of the key is increased, if it exceeds a certain length, thereby preventing key matches from being overlooked and improving search efficiency.
  • the upper computer 10 and the lower computer 20 each perform mapping according to the assigned query string, but this is not limited to the above.
  • only the upper computer 10 may perform a search by referring to the hierarchical index based on the assigned query string.
  • the string pair consists of a matched string in the reference string and a matching string in the query string, which includes a matching string identified by mapping.
  • the upper computer 10 and the lower computer 20 create a matched string by adding a corresponding string of a predetermined length in the reference string to the beginning and end of the matching string identified by mapping.
  • the matched string is a string in the vicinity of the matching string identified in the reference string, including the matching string.
  • the reference string is "CCGATCTGTATAACCCTACGA” and the matching string is "TACC”
  • the string "TATACCCT” obtained by adding two characters to the beginning and end becomes the matched string.
  • the length of bases added to the beginning and end of the matching string for the reference string can be, for example, about 50 bases each.
  • the host computer 10 also creates a comparison string by adding a corresponding string of a predetermined length in the query string to the end of the match string.
  • the length of the string to be added to the end of the match string for the query string is, for example, about 50 bases.
  • a string of a predetermined length is added to the beginning and end of the matching string, but this is not limited to the above, and for example, a string of a predetermined length may be added only to either the beginning or the end. Also, for a query string, a string of a predetermined length may be added to the beginning and end of the matching string, or the matching string itself may be treated as a matching string without adding any string.
  • the upper computer 10 and the lower computer 20 derive at least one approximate string based on a string pair consisting of a matched string based on the reference string and a match string based on the query string (S204).
  • a predetermined alignment using dynamic programming is applied to derive the approximate string.
  • Alignment is a method of comparing two sequences (strings) while allowing substitution, insertion, and deletion between their elements, and calculating the degree of variation (similarity) according to a defined score/penalty.
  • Known algorithms for achieving alignment include the Smith-Waterman algorithm and the Needleman-Wunsch algorithm.
  • Dynamic programming is a method of calculating the optimal solution for the next stage based on the optimal solution obtained at a certain stage. In other words, in calculating the degree of variation, a score is calculated for each element of the array-like alignment table according to dynamic programming, and the element with the maximum score is determined, thereby deriving at least one approximate string.
  • a hierarchical index is created based on a reference string, and then a matching string (and its length) is identified by searching the hierarchical index according to a given query string.
  • An approximate string matching is then performed on a string pair consisting of a matched string and a matching string based on the identified matching string, thereby deriving an approximate string.
  • FIG. 3 is a flowchart illustrating an example of an index creation process by a computer system according to an embodiment of the present invention. That is, FIG. 3 shows details of the hierarchical index creation process (S201) shown in FIG. 2. Also, FIG. 4 is a diagram showing an example of a reference string for creating a hierarchical index, and FIGS. 4A to 4C are diagrams showing an example of a data array structure in the process of creating a hierarchical index based on a reference string. Furthermore, FIG. 5 is a diagram showing an example of a table structure indicating the start position and number of occurrences of a key in the process of creating a hierarchical index, and FIG. 6 is a diagram for explaining a hierarchical index.
  • the host computer 10 receives a reference string for creating a hierarchical index (S301).
  • a reference string for creating a hierarchical index For example, in the case of a human genome reference sequence, the host computer 10 accesses the database 30 and reads out the stored human genome reference sequence.
  • the reference string "CCGATCTGTATAACCCTACGA” consisting of the four characters "A”, "C", "G", and "T” (see FIG. 4).
  • the host computer 10 extracts each substring (i.e., key) from the received reference string according to a specified key length to create a key arrangement (S302). For example, if the key length is "2", the key arrangement will be as shown in Figure 4A (a). In the figure, the number on the left is the array number. Also, to simplify the explanation, the key beginning with the 19th "A", which is the end of the reference string, has been omitted here.
  • the host computer 10 calculates a hash value for each key in the created key arrangement using a predetermined hash function and assigns this to the key (S303). This results in a key arrangement as shown in FIG. 4A(b).
  • each key array after sorting may include the array number before sorting.
  • the 0th key "AC” in the sorted key array holds the (original) array number "11" before sorting
  • the 1st key "AC” after sorting holds the array number "16" before sorting.
  • the calculated hash value is assigned to each extracted key, and the sorting is not limited to this.
  • a key array may be prepared in which hash values are calculated and assigned based on the combination of all characters appearing in the reference string, and the appearance start position corresponding to the extracted key may be assigned. That is, in a reference string consisting of four characters "A”, “C”, “G”, and "T", for example, if the length of the key is two characters, a key array having 42 elements is first created, and a hash value is further assigned to each of them. Next, the extracted key is assigned to the corresponding element (element of the same key) in the created key array together with the appearance start position in the reference string, thereby obtaining a key array as shown in FIG. 4A(c).
  • the host computer 10 identifies the appearance start position and the number of appearances of each key in the current key arrangement (S305).
  • Figure 5(a) shows the appearance start position and the number of appearances of each key in the key arrangement shown in Figure 4A(c). For example, it is shown that the key “AC” appears twice in the key arrangement, starting from the appearance start position "0" (sequence number "0"). Also, it is shown that the key “CC” appears three times, starting from the appearance start position "4".
  • the appearance start positions and the number of appearances for the key patterns consisting of all combinations of each character are shown, and for example, the key "AA" which is not included in the illustrated key arrangement is shown as the appearance start position "-" and the number of appearances "0".
  • the host computer 10 determines whether the number of occurrences of each key exceeds a predetermined tolerance (S306).
  • the predetermined tolerance is a value that allows for the same key to exist multiple times in the hierarchical index.
  • the predetermined tolerance is "1". In other words, if the predetermined tolerance is "1", each key will be unique in the hierarchical index. Furthermore, the larger the predetermined tolerance, the higher the possibility that duplicate keys will exist in the hierarchical index, while the faster the creation of the hierarchical index. If the host computer 10 determines that there is a key whose number of occurrences exceeds the predetermined tolerance (Yes in S306), it adds an additional key to that key (S308).
  • An additional key is one or more characters following the key in the original string.
  • the additional key is one character.
  • the substring obtained by adding the additional key is considered to be a new key.
  • the new key to which the additional key has been added is called an "extended key" to distinguish it from the original key, and the array is sometimes called an extended key array.
  • FIG. 4B (a) shows an example of an extended key array consisting of an extended key in which one additional key has been added to an original key. For the extended key "GA” with the array number "12”, since there is no character following "A” in the original string, for example, "$" is assigned as the termination string. For the keys with the array numbers "13", "17”, and "18”, since they appear only once, no additional key is added.
  • the host computer 10 further sorts each key (i.e., the extended keys) according to the additional key, for example in ascending order, as shown in FIG. 4B (b) (S308).
  • the sort order of the original keys takes precedence.
  • the host computer 10 then returns to step S306 and repeats the above process until the number of occurrences of all keys does not exceed the predetermined allowable value.
  • the host computer 10 identifies the start position of appearance and the number of appearances of each key in the key arrangement (S305), and if it determines that there is no key whose number of appearances exceeds a predetermined allowable value (No in S306), the desired hierarchical index has been created and the process ends.
  • FIG. 5 is a diagram for explaining the appearance start position and the number of appearances of a key in the process of creating a hierarchical index based on the reference string shown in FIG. 4.
  • an additional key is added to the key with an appearance count of 2 or more (FIG. 4B(a)), and each key (extended key) is sorted according to the additional key (FIG. 4B(b)).
  • FIG. 5(b) the appearance start position and the number of appearances of each key in the current key arrangement are identified.
  • the appearance count of the extended keys "CGA" and "TAC" is 2.
  • the host computer 10 creates a hierarchical index based on, for example, the human genome reference sequence. Because such a hierarchical index is sorted according to the hash value, the partial data array structure of the hierarchical index associated with a specific key (substring) is expanded in a concentrated manner in a specific address area in the main memory (data is made sequential), which makes it possible to significantly reduce the number of random accesses to the main memory.
  • an extremely short character string has been used as an example, but when dealing with the human genome, for example, it is preferable to set the key length of the substring to about 10-20, the key length of the additional key to 2-8, and the predetermined tolerance to 5-40, and it is even more preferable to set the key length of the substring to about 10-15, the key length of the additional key to 2-4, and the predetermined tolerance to 10-20.
  • FIG. 7 is a flowchart illustrating an example of a search process for a reference string based on a query string by a computer system according to one embodiment of the present invention. That is, FIG. 7 shows details of the search process (S202) for a reference string based on a query string shown in FIG. 2. This search process maps a specific key in the query string to the reference string, and identifies a matching string and its appearance start position in the reference string. Note that although the search process by the upper computer 10 is described below, the search process by the lower computer 20 operating in parallel is similar.
  • the host computer 10 reads the hierarchical index from the database 30, expands it in memory, and stores it (S701). From the perspective of increasing speed, the host computer 10 expands the data that constitutes the hierarchical index continuously in main memory and stores it. Here, expanding continuously includes arranging the data at consecutive memory addresses in the same bank.
  • the host computer 10 receives a query string for mapping to the reference string (S702).
  • the query string is a data fragment of base pairs read from a sequencer.
  • TACC the query string
  • the host computer 10 successively extracts each key from the received query string according to a predetermined key length (S703).
  • the initial value of the predetermined key length is "2". Therefore, the extracted keys are "TA”, “AC”, “CC”, and "C$".
  • the initial value of the sampling interval for each key is 1.
  • the sampling interval for a key is the interval between the start positions for sequentially referencing the hierarchical index according to the key. Note that in the case of human genome analysis, the initial value of the sampling interval can be, for example, around 3 to 5. As will be described later, the sampling interval is extended according to the length of the key.
  • a hash function is defined as a function that outputs a value expressed in base 4 using the numerical values "0" to "3" assigned to the four characters "A", “C”, “G”, and “T”, respectively.
  • the host computer 10 refers to the hierarchical index according to the calculated hash value (S705) and identifies the appearance start position and the number of appearances of each key in the reference string (S706). For example, for the key "TA”, the appearance start position in the hierarchical index is identified to be 14 and the number of appearances is identified to be 3 according to the hash value (see FIG. 6(a)).
  • the host computer 10 determines whether the number of occurrences of each key exceeds a predetermined threshold value (S707).
  • the predetermined tolerance is 1. If the host computer 10 determines that the number of occurrences of the key exceeds the predetermined tolerance value (Yes in S707), it then determines whether the current key length exceeds a predetermined upper limit value (S708).
  • the predetermined upper limit value may be around 20 in the case of human genome analysis, for example, but is not limited to this.
  • the host computer 10 determines that the current key length for each key does not exceed the predetermined upper limit (No in S708), it creates an extended key by adding an additional key to that key (S709) and returns to processing in S706.
  • the host computer 10 adds an additional key "C” to the key "TA” and checks the number of occurrences (see FIG. 6(b)), and then adds an additional key "C” to the key (extended key) to which the additional key has been added and checks the number of occurrences (see FIG. 6(c)), repeating this process until the number of occurrences becomes 1.
  • the host computer 10 determines that the current key length for each key exceeds the predetermined upper limit (Yes in S708), it extends (increases) the sampling interval for that key (S710) and returns to the process of S705. Note that the extension of the sampling interval may be performed only a predetermined number of times (for example, once) for each key.
  • the host computer 10 determines that the number of occurrences of the current key (i.e., the extended key) does not exceed a predetermined allowable value (No in S707), it outputs the key (matching string) in the identified reference string and its start position of occurrence (S711).
  • the host computer 10 can use the hierarchical index to search for at least a part of the query string from within the reference string.
  • a search according to the query string is performed on a partial, continuous data array structure in the hierarchical index, making it possible to significantly reduce the number of random accesses to the main memory.
  • the sampling interval for the start position of appearance is extended, which effectively reduces the number of searches and speeds up the search.
  • extending the sampling interval may increase the probability of missing a key that should have been identified, but since the sampling interval is lengthened when the key length exceeds a certain length, the occurrence of such oversights is effectively suppressed.
  • the host computer 10 executes alignment processing using dynamic programming to estimate the degree of variation (degree of variation) of the query string (matching string) with respect to the string (matched string) near the identified appearance start position in the reference string.
  • Alignment is a method of comparing two sequences (character strings) while allowing substitutions, insertions, and deletions between their elements, and calculating the degree of variation (similarity) according to a defined score/penalty. For example, in a comparison between character string X (i.e., a matched character string based on a reference character string) and character string Y (i.e., a matched character string based on a query character string), if some characters in character string Y are replaced, inserted, or deleted, character string Y is determined not to match character string X.
  • character string X i.e., a matched character string based on a reference character string
  • character string Y i.e., a matched character string based on a query character string
  • the relationship between characters at each position in character strings X and Y can be said to be either a match (match), a mismatch (mismatch), or one character not being present (gap).
  • the score for a match is defined as, for example, "+2”
  • the score for a mismatch as, for example, "-1”
  • the score for a gap as, for example, "-2”.
  • each element (character) of character string X and character string Y is compared, and the degree of variation of each element is calculated using these scores.
  • an improved alignment algorithm based on the Needleman-Wunsch algorithm, which is an example of global alignment, is used to calculate the degree of variation. Below, we explain the basic concept of the Needleman-Wunsch algorithm, and then explain the improved alignment algorithm applied to the present invention.
  • the degree of variation F(i,j) between a character x i and a character y j is defined as follows.
  • max is a function that outputs the maximum value among the values of a given formula
  • FIG. 8A is an alignment table that arranges the strings to be compared in this example.
  • " ⁇ " is a null character
  • F(0,0) 0.
  • the argument for the mutation degree F(i,j) is a negative value, it is outside the range and the calculation is omitted (the output is Null).
  • the alignment table is expanded in memory as a type of data structure and made available for use by the processor.
  • an alignment table as shown in FIG. 8D is created.
  • the mutation degree F(6,7) of element position (6,7) has the maximum value "7". Therefore, as shown in the figure, backtracking is performed from element position (6,7) to element position (1,1).
  • an arrow pointing diagonally downward to the right indicates a match of characters
  • an arrow pointing to the right indicates a missing character
  • an arrow pointing downward indicates an insertion.
  • this disclosure proposes the following improved alignment algorithm, which reduces the amount of calculation and speeds up processing.
  • the improved alignment algorithm applied to the present invention involves, broadly speaking, defining a calculation region of a prescribed width centered on an element located on the diagonal of the alignment table, calculating the degree of variation only for elements within the calculation region to determine its maximum value, and deriving an approximate string from the element based on the maximum value if the maximum value satisfies a prescribed condition. If the maximum value does not satisfy the prescribed condition, the calculation region is expanded, and similarly, the degree of variation is calculated to determine the maximum value, and this is repeated until the maximum value satisfies the prescribed condition.
  • FIG. 9 is a flowchart illustrating an example of alignment processing using dynamic programming by a computer system according to one embodiment of the present invention. That is, FIG. 9 shows details of the alignment processing (S204) shown in FIG. 2. Note that, although the alignment processing based on one query string by the upper computer 10 is described below, the same applies to alignment processing based on other query strings by the lower computer 20 operating in parallel.
  • the host computer 10 creates an alignment table based on string pairs consisting of a matched string and a matching string (S1001). As described above, the alignment table is expanded in memory as a type of data structure.
  • the host computer 10 sets the width m (where m is a positive number) of the calculation area to an initial value (S1002).
  • the initial value of width m can be, for example, the length of the matching string obtained by mapping, but is not limited to this.
  • width m is set to an odd number value since it is centered on the diagonal element position of the alignment table, but is not limited to this. This determines a calculation area of width m centered on the diagonal element position of the alignment table.
  • the host computer 10 sets a predetermined dummy value to each element that defines the boundary of the calculation domain (S1003).
  • the dummy value is selected to be a value that is sufficiently small as the value of the degree of mutation F.
  • the dummy value is a negative value, for example, about 2 to 3 times the width m of the initial value, and can be set arbitrarily.
  • the host computer 10 calculates the mutation degree F for each element in the calculation domain according to Equation 1 (S1004). Then, the host computer 10 determines the maximum mutation degree Fmax that has the maximum value in the calculation domain and identifies the position of the element (S1005). Note that the element having the maximum mutation degree Fmax is not necessarily one.
  • the host computer 10 also calculates a lower limit F Low of the degree of variation F in a row or column of the alignment table (S1006).
  • the host computer 10 compares the maximum mutation degree Fmax with the lower limit value FLow to determine whether the maximum mutation degree Fmax exceeds the lower limit value FLow (S1007). If the host computer 10 determines that the maximum mutation degree Fmax does not exceed the lower limit value FLow (No in S1007), it expands the width m by a predetermined size ⁇ (S1008). For example, ⁇ is set to m+1 (however, m does not exceed the number of characters in the string).
  • the host computer 10 performs the same process on the expanded calculation region of width m, and repeats the above process until the maximum value F max exceeds the lower limit value F low .
  • the host computer 10 determines that the maximum value Fmax exceeds the lower limit Flow (Yes in S1007), it performs backtracking from the position of the element having the maximum value to determine an approximate character string (S1009).Then, the host computer 10 outputs the determined approximate character string (S1010).
  • Equation 1 is used to calculate the degree of mutation F. Note that the length of the matched string and the length of the matching string do not need to be the same, and typically the length of the matched string is longer than the length of the matching string.
  • the host computer 10 prepares an alignment table in which the strings to be compared are arranged, and sets an initial value for width m.
  • the initial value of width m(0) is set to 7.
  • the host computer 10 also sets dummy values to elements at the boundary of the calculation area defined according to width m.
  • the dummy value is set to -20.
  • the host computer 10 calculates the mutation degree F of each element in the calculation region R(0) according to Equation 1.
  • FIG. 10B shows the state in which the mutation degree F of each element in the calculation region has been calculated.
  • the host computer 10 determines the maximum mutation degree F max from the calculated mutation degrees F. In this example, the maximum mutation degree F max is 17. In the figure, elements with a maximum mutation degree F max of 17 are shown hatched.
  • the host computer 10 calculates the lower limit F Low for a row or column of the alignment table.
  • the host computer 10 compares the maximum mutation degree Fmax with the lower limit value FLow , and determines that the maximum mutation degree Fmax does not exceed the lower limit value FLow , and therefore expands the width m by ⁇ .
  • the expanded width m(1) is expanded to 15.
  • the calculation region of the expanded width m(1) is defined as R(1).
  • the host computer 10 calculates the mutation degree F of each element in the calculation region R(1) according to Equation 1.
  • FIG. 10C shows the state in which the mutation degree F of each element in the calculation region has been calculated.
  • the region that has been added to the calculation region R(1) by widening is hatched.
  • the maximum mutation degree F max becomes 19.
  • the lower limit value F low at this time becomes 18.
  • the improved alignment algorithm does not calculate the degree of mutation for all elements in the alignment table, but instead defines a calculation region of a specified width and calculates the degree of mutation while expanding the calculation region as necessary, thereby reducing the amount of calculation and thereby speeding up processing.
  • a hierarchical index based on a reference string is created, and then a matching string (and its length) is identified by searching the hierarchical index according to a given query string.
  • An approximate string match is then performed on a string pair consisting of a matched string and a matching string based on the identified matching string, thereby deriving an approximate string.
  • steps, operations, or functions may be performed in parallel or in a different order, provided that the results are not inconsistent.
  • the steps, operations, and functions described are provided merely as examples, and some of the steps, operations, and functions may be omitted or combined into one, or other steps, operations, or functions may be added, without departing from the spirit of the invention.

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

本発明は、クエリ文字列に基づいて参照文字列における近似文字列を検索する方法である。前記方法は、前記参照文字列に基づいて階層的インデックスを作成することと、前記クエリ文字列の少なくとも一部と一致する前記参照文字列における部分文字列を同定するために、前記階層的インデックスを参照して、前記参照文字列に対するクエリ文字列のマッピングを行うことと、前記マッピングにより同定される少なくとも１以上の前記部分文字列に基づいて、前記近似文字列を導出することを含む。前記階層的インデックスは、前記参照文字列から切り出される各キーに対して、その出現回数に従い追加キーを追加しソートすることを繰り返しながら作成される。

Description

近似文字列照合方法及び該方法を実現するためのコンピュータプログラム
　本発明は、近似文字列照合技術に関し、特に、ヒトゲノムの解析に利用可能な近似文字列照合装置及び近似文字列照合方法並びに該方法を実現するためのコンピュータプログラムに関する。
　ヒトゲノムは、人が持つ遺伝情報のセットであり、これを担っている物質が、約３０億対の塩基が連なったＤＮＡ（デオキシリボ核酸）である。塩基は、アデニン（Ａ）、グアニン（Ｇ）、シトシン（Ｃ）、及びチミン（Ｔ）がある。すなわち、人の遺伝情報は、これらの塩基の並び（配列）によって決定される。
　ヒトゲノムの読み取りにはシークエンサと称される装置が用いられる。シークエンサは、サンプルとなるヒトゲノムを読み取って、これを所定の上限値（数百塩基対程度）に細断して増幅し、データ片からなる膨大なデータ配列として出力する。現行のシークエンサは、一人分のヒトゲノムを１時間ほどで読み出すことができる。
　シークエンサにより読み出されるばらばらのデータ片は、人の標準的なゲノム配列として定められたヒトゲノム参照配列と比較されることによって、元の長さのヒトゲノム配列に再構築され解析される。例えば、各データ片が、ヒトゲノム参照配列との比較において、どの位置にあるかが調べられ（マッピング）、また、どのような変異があるかといった解析がなされる。通常、シーケンサから読み出されるデータには誤差が含まれるため、１サンプルあたりヒトゲノム配列一人分の例えば３０倍の冗長データを用いて統計処理を行うことにより誤差を小さくしている。したがって、ヒトゲノム配列の解析には膨大な計算量が必要とされるため、典型的には、スーパーコンピュータやクラスタコンピュータ、ＦＰＧＡ（Field-Programmable Gate Array）ベースのコンピュータといった高性能なコンピュータが用いられる。
　データ片のマッピングには、例えばＢＷＡ（Burrow-Wheeler Aligner）といったプログラムツールが用いられる。ＢＷＡは、３つのアルゴリズム、ＢＷＡ－ｂａｃｋｔｒａｃｋ、ＢＷＡ－ＳＷ、及びＢＷＡ－ＭＥＭから構成される。このうち、ＢＷＡ－ＭＥＭは、Indel（挿入欠失）に対応した高速アルゴリズムとして広く利用されている。ＢＷＡ－ＭＥＭは、読み出したデータ片に基づくクエリ文字列のうち、ヒトゲノム参照配列に繰り返し現れる部分に対して、接尾辞配列を用いてインデックスを作成し、マッピングを行うアルゴリズムである（非特許文献１）。
　また、データ片は、マッピングにより得られるヒトゲノム参照配列の部分配列と照合され、どのような変異があるかが解析される。ヒトゲノム参照配列の部分配列とデータ片との照合には、例えば、アラインメントアルゴリズムが用いられる。アラインメントアルゴリズムでは、動的計画法に従って、アラインメント表と称される配列の各要素について変異度（類似度）を算出しながら、近似文字列を同定する（非特許文献２）。
Heng Li, "Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM", May 26, 2013, arXiv:1303.3997 (q-bio.GN) 内山将夫　他，「近似文字列照合による全文検索のための接尾辞配列の高速走査法」，２００２年９月１５日，情報処理学会，Vol. 43 No. SIG 9(TOD 15)
　上述したＢＷＡ－ＭＥＭは、ギャップを許容しかつ高速マッピングが可能なアルゴリズムとして広く利用されているが、作成されるインデックスのサイズが非常に大きくなるため、大量のメモリリソースを必要とするという問題があった。また、ＢＷＡ－ＭＥＭでは、クエリ文字列の文字数（すなわち、塩基数）に応じた回数だけメモリへのアクセスが必要となるため、メモリへのランダムアクセスが頻繁に発生し、これにより、プロセッサの処理速度がメモリへのアクセス時間により律速されてしまうという問題があった。
　より具体的には、ＢＷＡ－ＭＥＭでは、クエリ文字列における部分文字列（これを「キー」と称することがある。）がヒトゲノム参照配列内に出現する位置を同定するために、メインメモリへのランダムアクセスがクエリ文字列の総文字数分の回数だけ実行される。これは、ヒトゲノムの解析に用いるデータの量が非常に大きいため、必要なデータ配列をより高速アクセス可能なキャッシュメモリに一度に収容しきれないからである。ここで、メインメモリへの１回のランダムアクセスの待ち時間（アクセスが発生してから実際にデータが得られるまでの時間）は、現在のコンピュータでは、１バンクあたり約１μ秒（０．０００００１秒）かかっている。したがって、３０億塩基対のヒトゲノムの解析の場合、冗長性を３０とすると、総文字数に対する総アクセス時間Ｔは、
　　Ｔ＝３，０００，０００，０００×０．０００００１×３０
　　　＝９０，０００（秒）
となる。つまり、ヒトゲノムの解析において、メインメモリへのランダムアクセスだけで、約９万秒（約２５時間）かかることになる。このため、たとえ、高性能なプロセッサを備えたコンピュータを用いたとしても、メモリアクセス時間が制約となって、マッピング時間を短縮するには限界があった。
　また、上述したように、シークエンサにより一人分のヒトゲノムを読み出すためには、現状、１時間ほど要している。一方で、高性能コンピュータを用いてＢＷＡ－ＭＥＭを実行した場合、シークエンサの読み出し時間以内にマッピングが完了し、両者のバランスは概ね保たれている。
　一方で、次世代型シークエンサは、低コスト化が進み、また、読み出し時間を更に短縮し得ると言われており、これに伴って、解析時間もまた短縮することが望まれる。解析時間を短縮するための一つのアプローチとして、コンピュータの更なる高性能化が考えられるが、コンピュータの高性能化のためには非常にコストが嵩むため実用化へのハードルが高い。
　更に、マッピングされたデータ片の解析に用いられる従前のアラインメントアルゴリズムでは、アラインメント表の全ての要素について変異度を算出するため、計算量が多くなり、時間がかかるという問題がある。
　そこで、本発明は、参照データを用いて、与えられるクエリデータの解析を高速及び／又は効率的に行うことができる新たな技術を提供することを目的とする。
　より具体的には、本発明の一つの目的は、与えられるクエリ文字列と参照文字列との間の近似文字列照合を高速及び／又は効率的に行うことができる近似文字列照合装置及びこれを用いた近似文字列照合方法を提供することである。
　また、本発明の一つの目的は、シークエンサによって読み出されたデータ片に基づくヒトゲノム参照配列を用いた解析を高速及び／又は効率的に行うことができる近似文字列照合装置及びこれを用いた近似文字列照合方法を提供することである。
　また、本発明の一つの目的は、前記近似文字列照合に適合した参照文字列に基づく階層的インデックスを作成する技術を提供することである。
　上記課題を解決するための本発明は、以下に示す発明特定事項乃至は技術的特徴を含んで構成される。
　ある観点に従う本発明は、コンピューティングデバイスに、クエリ文字列に基づいて参照文字列における近似文字列を検索するための方法を実現させるためのコンピュータプログラムである。前記方法は、前記参照文字列に基づいて階層的インデックスを作成することと、前記クエリ文字列の少なくとも一部と一致する前記参照文字列における部分文字列を同定するために、前記階層的インデックスを参照して、前記参照文字列に対する前記クエリ文字列のマッピングを行うことと、前記のマッピングにより同定される前記部分文字列に基づいて作成される、前記参照文字列に基づく被照合文字列と前記クエリ文字列に基づく照合文字列とに基づいて、前記被照合文字列に近似する文字列を前記近似文字列として導出することと、を含む。ここで、前記階層的インデックスを作成することは、前記参照文字列から所定長の各第１のキーを切り出すことと、切り出された前記各第１のキーについて、所定のハッシュ関数により該第１のキーに基づいて算出されるハッシュ値を割り当てた第１のキー配列を作成することと、作成された前記第１のキー配列を更新することと、更新された前記第１のキー配列を前記階層的インデックスとして出力することと、を含む。また、前記第１のキー配列を更新することは、前記第１のキー配列における前記各第１のキーについて、前記参照文字列における該第１のキーの出現回数を同定することと、同定された前記第１のキーの前記出現回数に従って、該第１のキーに対して前記参照文字列における該第１のキーに続く少なくとも１以上の文字からなる第１の追加キーを追加することにより前記第１のキー配列を更新することと、を含む。
　前記第１のキー配列を作成することは、前記ハッシュ値に従って前記第１のキー配列における前記各第１のキーをソートすることを含み得る。
　また、前記第１のキー配列を更新することは、前記同定された出現回数が所定の許容値を超えているか否かを判断することと、前記同定された出現回数が前記所定の許容値を超えていると判断される場合に、前記第１のキーに対して前記参照文字列における該第１のキーに続く少なくとも１以上の文字からなる前記第１の追加キーを追加することにより新たな第１のキーを作成することと、前記新たな第１のキーについて、前記参照文字列における該新たな第１のキーの出現回数を同定することと、を含み得る。
　また、前記第１のキー配列を更新することは、前記第１の追加キーに従って前記第１のキー配列における前記新たな第１のキーをソートすることを更に含み得る。
　また、前記第１のキー配列を更新することは、前記同定された前記出現回数が所定の許容値を超えていないと判断されるまで、現在の前記第１のキーに新たな前記第１の追加キーを順次に追加することにより新たな前記第１のキーを作成することを含み得る。
　また、前記第１のキー配列を前記階層的インデックスとして出力することは、前記同定された前記出現回数が所定の許容値を超えていないと判断される場合に、現在の前記第１のキー配列を前記階層的インデックスとして出力することを含み得る。
　前記マッピングを行うことは、前記クエリ文字列から所定長の各第２のキーを切り出すことと、前記クエリ文字列から切り出された前記各第２のキーについて、前記所定のハッシュ関数により該第２のキーに基づいて算出されるハッシュ値を割り当てた第２のキー配列を作成することと、前記各第２のキーについて、前記ハッシュ値に従って、所定のサンプリング間隔で、前記階層的インデックスを参照し、該第２のキーの出現開始位置及び出現回数を同定することと、を含み得る。
　前記第２のキーの前記出現開始位置及び前記出現回数を同定することは、前記第２のキーの前記出現回数が前記所定の許容値を超えているか否かを判断することと、前記第２のキーの前記出現回数が前記所定の許容値を超えていると判断される場合に、前記第２のキーに対して前記クエリ文字列における該第２のキーに続く少なくとも１以上の文字からなる第２の追加キーを追加することにより新たな第２のキーを作成することと、前記第２のキーの前記出現回数が前記所定の許容値を超えていないと判断される場合に、同定された現在の前記第２のキーを一致文字列として出力するとともに該第２のキーの前記出現開始位置を出力することと、を含み得る。
　また、前記第２のキーの前記出現開始位置及び前記出現回数を同定することは、前記第２のキーの前記同定された前記出現回数が前記所定の許容値を超えていないと判断されるまで、現在の前記第２のキーに新たな前記第２の追加キーを順次に追加して、前記新たな第２のキーを作成することを更に含み得る。
　また、前記第２のキーの前記出現回数が前記所定の許容値を超えていると判断される場合に、該第２のキーの前記所定のサンプリング間隔が大きくなるように変更され得る。
　前記近似文字列として導出することは、前記マッピングにより同定された前記部分文字列に基づく、前記被照合文字列と前記照合文字列とからなる文字列ペアを受信することと、前記文字列ペアに基づいて少なくとも１つの前記近似文字列を導出するために、所定のアラインメント処理を実行することと、導出された前記少なくとも１つの近似文字列を出力することと、を含み得る。
　前記所定のアラインメント処理を実行することは、前記被照合文字列と前記照合文字列とに基づいて所定のアラインメント表を作成することと、前記アラインメント表の対角線上の要素を中心にした幅ｍを有する計算領域を設定することと、設定された前記計算領域における各要素について、変異度を算出することと、算出された前記変異度に基づいて、最大変異度を決定することと、決定された前記最大変異度に基づいて、前記少なくとも１つの近似文字列を導出することを含み得る。
　また、前記所定のアラインメント処理を実行することは、前記最大変異度と所定の下限値とを比較して、前記最大変異度が前記所定の下限値を超えているかを判断することと、前記最大変異度が前記所定の下限値を超えていないと判断される場合に、新たな計算領域を設定するために、前記計算領域の前記幅ｍを拡幅することと、前記最大変異度が前記所定の下限値を超えていると判断される場合に、前記最大変異度を有する要素に基づいて、前記少なくとも１つの近似文字列を導出することと、を含み得る。そして、前記方法は、前記最大変異度が前記下限値を超えると判断されるまで、前記計算領域を拡幅することにより新たな計算領域を設定して前記変異度を算出することを繰り返すように構成され得る。
　前記所定の下限値は、所定の要素列にｍ個の連続したギャップがあり、それ以外の部分は完全又は実質的に一致したと仮定した場合の変異度の値であり得る。
　また、前記方法は、前記部分文字列に対して前記参照文字列における対応する所定の文字列を追加することにより前記被照合文字列を作成することと、前記部分文字列に対して前記クエリ文字列における対応する所定の文字列を追加することにより前記照合文字列を作成することと、を更に含み得る。
　また、ある観点に従う本発明は、コンピューティングデバイスに、クエリ文字列に基づいて参照文字列を探索するための階層的インデックスを作成する方法を実行させるためのコンピュータプログラムである。
　前記方法は、前記参照文字列から所定長の各第１のキーを切り出すことと、切り出された前記各第１のキーについて、所定のハッシュ関数により該第１のキーに基づいて算出されるハッシュ値を割り当てた第１のキー配列を作成することと、作成された前記第１のキー配列を更新することと、更新された前記第１のキー配列を前記階層的インデックスとして出力することと、を含む。
　ここで、前記第１のキー配列を更新することは、前記第１のキー配列における前記各第１のキーについて、前記参照文字列における該第１のキーの出現開始位置及び出現回数を同定することと、同定された前記第１のキーの前記出現開始位置及び前記出現回数に従って、該第１のキーに対して前記参照文字列における該第１のキーに続く少なくとも１以上の文字からなる第１の追加キーを追加することにより前記第１のキー配列を更新することと、を含む。
　また、ある観点に従う本発明は、コンピューティングデバイスに、参照文字列に対してクエリ文字列のマッピングを行う方法を実行させるためのコンピュータプログラムである。前記方法は、前記参照文字列から切り出される各第１のキー、及び前記各第１のキーについて、該第１のキーに基づいて所定のハッシュ関数により算出されるハッシュ値が割り当てられた第１のキー配列であって、前記参照文字列における前記第１のキーの出現開始位置及び出現回数に従って前記第１のキーに追加キーが追加された該第１のキー配列からなる階層的インデックスを読み出すことと、前記クエリ文字列から所定のキー長を有する各第２のキーを切り出すことと、前記クエリ文字列から切り出された前記各第２のキーについて、所定のハッシュ関数により該第２のキーに基づいて算出されるハッシュ値を割り当てた第２のキー配列を作成することと、前記各第２のキーについて、前記ハッシュ値に従って、所定のサンプリング間隔で、前記階層的インデックスを参照し、前記第１のキー配列における前記第１のキーとの比較により、該第２のキーの出現開始位置及び出現回数を同定することと、前記同定した出現回数が所定のしきい値を超えているか否かを判断することと、前記同定された前記出現回数が所定の許容値を超えていると判断される場合に、前記第２のキーに対して前記クエリ文字列における該第２のキーに続く少なくとも１以上の文字からなる追加キーを追加することにより新たな第２のキーを作成することと、前記同定された前記出現回数が所定のしきい値を超えていないと判断される場合に、同定された現在の前記第２のキーの出現開始位置及び該キーを出力することと、を含む。そして、前記第２のキーの前記出現開始位置及び前記出現回数を同定することは、前記同定された前記出現回数が所定のしきい値を超えていないと判断されるまで、現在の前記第２のキーに新たな前記追加キーを順次に追加して、前記新たな第２のキーを作成することを含む。
　ここで、前記同定された前記出現回数が所定の許容値を超えていると判断される場合に、前記第２のキーの前記所定のサンプリング間隔が大きくなるように変更され得る。
　また、ある観点に従う本発明は、コンピューティングデバイスに、参照文字列における部分文字列とクエリ文字列との間の変異を所定のアラインメント処理により同定する方法を実現させるためのコンピュータプログラムである。前記方法は、マッピングにより同定される前記部分文字列に基づいて作成される、前記参照文字列に基づく被照合文字列と前記クエリ文字列に基づく照合文字列とからなる文字列ペアを受信することと、前記文字列ペアに基づいて、前記被照合文字列に近似する文字列を少なくとも１つの近似文字列として導出するために、所定のアラインメント処理を実行することと、導出された前記少なくとも１つの近似文字列を出力することと、を含む。そして、前記所定のアラインメント処理を実行することは、前記被照合文字列と前記照合文字列とに基づいて所定のアラインメント表を作成することと、前記アラインメント表の対角線上の要素を中心にした幅ｍを有する計算領域を設定することと、設定された前記計算領域における各要素について、変異度を算出することと、算出された前記変異度に基づいて、最大変異度を決定することと、決定された前記最大変異度に基づいて、前記少なくとも１つの近似文字列を導出することを含む。
　また、前記所定のアラインメント処理を実行することは、前記最大変異度と所定の下限値とを比較して、前記最大変異度が前記所定の下限値を超えているかを判断することと、前記最大変異度が前記所定の下限値を超えていないと判断される場合に、新たな計算領域を設定するために、前記計算領域の前記幅ｍを拡幅することと、前記最大変異度が前記所定の下限値を超えていると判断される場合に、前記最大変異度を有する要素に基づいて、前記少なくとも１つの近似文字列を導出することと、を含み得る。そして、前記最大変異度が前記下限値を超えると判断されるまで、前記計算領域を拡幅することにより新たな計算領域を設定して前記変異度を算出することが繰り返され得る。
　また、前記所定のアラインメント処理を実行することは、所定の要素列にｍ個の連続したギャップがあり、それ以外の部分は完全に又は実質的に一致したと仮定した場合の変異度の値を前記所定の下限値として設定することを更に含み得る。
　なお、本発明は、前記コンピュータプログラムを非一時的に記録するコンピュータ可読記録媒体としても成立する。また、本発明は、前記方法を実行するように構成されたハードウェア及び／又はファームウェアからなる装置としても成立する。近似文字列照合装置は、本発明の一形態である。
　なお、本明細書等において、手段又は部（unit）とは、単に物理的手段を意味するものではなく、その手段又は部が有する機能をソフトウェアによって実現する場合も含む。また、１つの手段又は部が有する機能が２つ以上の物理的手段により実現されても、２つ以上の手段又は部の機能が１つの物理的手段により実現されても良い。また、「システム」とは、複数の装置（又は特定の機能を実現する機能モジュール）が論理的に集合した物のことをいい、各装置や機能モジュールが単一の筐体内にあるか否かは特に問わない。
　本発明によれば、与えられるクエリデータに対する参照データを用いた解析を高速及び／又は効率的に行うことができる。とりわけ、本発明によれば、与えられるクエリ文字列と参照文字列との間の近似文字列照合を高速及び／又は効率的に行うことができる。
　また、本発明によれば、シークエンサによって読み出されたデータ片に基づくヒトゲノム参照配列を用いた解析を高速及び／又は効率的に行うことができる。
　更に、本発明によれば、近似文字列照合に適合した参照文字列に基づく階層的インデックスを提供することができる。
　本発明の他の技術的特徴、目的、及び作用効果乃至は利点は、添付した図面を参照して説明される以下の実施形態により明らかにされる。本明細書に記載された効果はあくまで例示であって限定されるものではなく、また他の効果があっても良い。
図１は、本発明の一実施形態に係るコンピュータシステムの概略的構成の一例を示すブロックダイアグラムである。 図２は、本発明の一実施形態に係るコンピュータシステムによる近似文字列照合処理の概略の一例を説明するフローチャートである。 図３は、本発明の一実施形態に係るコンピュータシステムによるインデックス作成処理の一例を説明するフローチャートである。 図４は、本発明の一実施形態に係るコンピュータシステムにおいて用いられる参照文字列の一例を示す図である。 図４Ａは、図４に示される参照文字列に基づく階層的インデックスの作成過程におけるデータ配列構造の一例を説明するための図である。 図４Ｂは、図４に示される参照文字列に基づく階層的インデックスの作成過程におけるデータ配列構造の一例を説明するための図である。 図４Ｃは、図４に示される参照文字列に基づく階層的インデックスの作成過程におけるデータ配列構造の一例を説明するための図である。 図５は、図４に示される参照文字列に基づく階層的インデックスの作成過程におけるキー開始位置及び出現回数を示すテーブル構造の一例を示す図である。 図６は、本発明の一実施形態に係るコンピュータシステムにおいて用いられる参照文字列の一例を示す図である。 図７は、本発明の一実施形態に係るコンピュータシステムによるクエリ文字列に基づく参照文字列の探索処理の一例を説明するフローチャートである。 図８Ａは、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムを説明するためのアラインメント表の一例を示す図である。 図８Ｂは、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムを説明するためのアラインメント表の一例を示す図である。 図８Ｃは、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムを説明するためのアラインメント表の一例を示す図である。 図８Ｄは、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムを説明するためのアラインメント表の一例を示す図である。 図９は、本発明の一実施形態に係るコンピュータシステムによる動的計画法を用いたアラインメント処理の一例を説明するフローチャートである。 図１０Ａは、本発明の一実施形態に係る動的計画法を用いたアラインメントを説明するためのアラインメント表の一例を示す図である。 図１０Ｂは、本発明の一実施形態に係る動的計画法を用いたアラインメントを説明するためのアラインメント表の一例を示す図である。 図１０Ｃは、本発明の一実施形態に係る動的計画法を用いたアラインメントを説明するためのアラインメント表の一例を示す図である。 図１１Ａは、本発明の一実施形態に係る動的計画法を用いたアラインメントにより得られる近似文字列の一例を示す図である。 図１１Ｂは、本発明の一実施形態に係る動的計画法を用いたアラインメントにより得られる近似文字列の一例を示す図である。 図１２は、本発明の一実施形態に係るシステムにおけるコンピューティングデバイスのハードウェア構成の一例を示す図である。
　以下、図面を参照して本発明の実施の形態を説明する。ただし、以下に説明する実施形態は、あくまでも例示であり、以下に明示しない種々の変形や技術の適用を排除する意図はない。本発明は、その趣旨を逸脱しない範囲で種々変形（例えば各実施形態を組み合わせる等）して実施することができる。また、以下の図面の記載において、同一又は類似の部分には同一又は類似の符号を付して表している。図面は模式的なものであり、必ずしも実際の寸法や比率等とは一致しない。図面相互間においても互いの寸法の関係や比率が異なる部分が含まれていることがある。
　図１は、本発明の一実施形態に係るコンピュータシステムの概略的構成の一例を示すブロックダイアグラムである。同図に示すように、コンピュータシステム１は、例えば、上位コンピュータ１０と、複数の下位コンピュータ２０（１）～２０（ｎ）と、データベース３０と含み構成される。上位コンピュータ１０と下位コンピュータ２０（１）～２０（ｎ）とは、所定のインターフェースを介して通信可能に接続される。本開示では、コンピュータシステム１は、上位コンピュータ１０が複数の下位コンピュータ２０（１）～２０（ｎ）を統括的に制御する中央集権型コンピュータシステムとして構成されるが、これに限られず、例えば分散型コンピュータシステムとして構成されても良い。分散型コンピュータシステムにおいては、各コンピュータが並列分散処理により協調的に動作し得るが、特定の処理に関しては、代表する一のコンピュータのみが該処理を実行する場合があっても良い。上位コンピュータ１０及び下位コンピュータ２０のハードウェア構成は、図１２に例示されるが、当業者にとって自明であるため、その詳細な説明は省略する。本開示では、下位コンピュータ２０は、上位コンピュータ１０の制御の下、並列的にタスクを処理する。例えば、複数の下位コンピュータ２０のそれぞれは、与えられた個々のクエリ文字列に基づいて参照文字列との比較において解析処理を行う。以下では、下位コンピュータ２０（１）～２０（ｎ）について、それらを特に区別する必要がない限り、単に、「下位コンピュータ２０」と表記することがある。
　データベース３０は、上位コンピュータ１０の制御の下、各種のデータ、例えば参照文字列を格納する。一例として、データベース３０は、ヒトゲノム参照配列を格納する。ヒトゲノム参照配列は、人の標準的なゲノム配列として定められた塩基対の配列を示すデータである。また、データベース３０は、参照文字列に基づいて作成されたインデックスを格納する。インデックスは、参照文字列を探索するために用いられるある種のデータ配列構造である。なお、図中、データベース３０は、上位コンピュータ１０にのみアクセス可能に接続される構成となっているが、これに限られず、下位コンピュータ２０もまたアクセス可能に接続される構成であっても良い。
　図２は、本発明の一実施形態に係るコンピュータシステムによる近似文字列照合処理の概略の一例を説明するフローチャートである。かかる処理は、例えば、上位コンピュータ１０及び複数の下位コンピュータ２０が、プロセッサの制御の下、近似文字列照合プログラムを実行することにより他のハードウェアコンポーネントと協働し、実現される。
　すなわち、同図に示すように、まず、上位コンピュータ１０は、データベース３０に格納された参照文字列を読み出し、読み出した参照文字列に基づいて、インデックスを作成する（Ｓ２０１）。具体的には、上位コンピュータ１０は、参照文字列における部分文字列を所定の条件に従って拡張し及び／又はソートして配列化することにより、所定のデータ配列構造を有するインデックスを作成する。データ配列構造は、階層的なツリー構造からなる。本開示では、このようなインデックスを階層的インデックスと称するものとする。上位コンピュータ１０は、作成した参照文字列に基づくインデックスをデータベース３０に格納する。参照文字列に基づく階層的インデックスの作成処理の詳細は後述される。なお、参照文字列に基づく階層的インデックスが既に作成されデータベース３０に格納されている場合には、処理Ｓ２０１は省略され得る。
　続いて、上位コンピュータ１０は、図示しない外部装置からクエリ文字列を受信し、受信したクエリ文字列に基づいて、階層的インデックスを参照して、マッピングを行う（Ｓ２０２）。一例では、外部装置は、ヒトゲノムを読み出すシークエンサであり、クエリ文字列は、シークエンサから出力される例えば１５０～２００塩基程度の塩基対のデータ片である。本開示では、下位コンピュータ２０もまた、上位コンピュータ１０の制御の下、割り当てられたクエリ文字列に基づいて、階層的インデックスを参照して、マッピングを行う。
　マッピングは、クエリ文字列の少なくとも一部と一致する参照文字列における部分文字列（一致文字列）を同定する処理である。つまり、マッピングにより、参照文字列におけるクエリ文字列の少なくとも一部と一致する部分文字列の出現開始位置及び長さ（文字数）が同定される。本開示に係るマッピングでは、クエリ文字列における各キーについて、階層的インデックスを検索ないしは探索することにより、参照文字列におけるクエリ文字列の各キー及びその出現開始位置が同定される。なお、処理の高速化のため、作成された階層的インデックスは、メインメモリ又はキャッシュメモリ等の高速メモリに展開され得る。
　また、本開示に係るマッピングでは、階層的インデックスの探索に際して、クエリ文字列におけるキーの出現開始位置のサンプリング間隔が従来技術に比較してある程度大きくなるように調整される。例えばＢＬＡＳＴ（Basic Local Alignment Search Tool）と称されるＤＮＡの塩基配列やタンパク質のアミノ酸配列のシーケンスアラインメントを行うためのアルゴリズムでは、キーの出現開始位置の間隔は３～５文字に設定されるが、本開示に係る探索では、例えば１０文字又はそれ以上に設定され得る。探索に際して、キーの出現開始位置のサンプリング間隔が小さいと、探索の回数が増え、効率が低下してしまうのに対して、キーの出現開始位置のサンプリング間隔を常に大きくしてしまうと、探索の高速化と引き換えに、見落としの確率が上昇してしまう。そこで、本開示では、キーの長さを長くしながら、それが一定長を超える場合には、出現開始位置のサンプリング間隔を大きくすることにより、キーの一致を見落とすことを防ぐとともに、探索効率の向上を図っている。
　なお、上記の例では、上位コンピュータ１０及び下位コンピュータ２０が、それぞれ、割り当てられたクエリ文字列に従ってマッピングを行っているが、これに限られず、例えば、上位コンピュータ１０のみが、割り当てられたクエリ文字列に基づいて、階層的インデックスを参照して、探索を行っても良い。
　次に、上位コンピュータ１０及び下位コンピュータ２０は、それぞれ、後述する近似文字列照合のための文字列ペアを作成する（Ｓ２０３）。文字列ペアは、マッピングにより同定される一致文字列を含む、参照文字列における被照合文字列とクエリ文字列における照合文字列とからなる。
　具体的には、上位コンピュータ１０及び下位コンピュータ２０は、マッピングにより同定された一致文字列の先頭及び末尾のそれぞれに、参照文字列における対応する所定長の文字列を追加することにより、被照合文字列を作成する。つまり、被照合文字列は、参照文字列において同定された一致文字列を含む該一致文字列近傍の文字列である。例えば、参照文字列が「ＣＣＧＡＴＣＴＧＴＡＴＡＣＣＣＴＡＣＧＡ」であって、一致文字列が「ＴＡＣＣ」である場合に、例えば前後２文字ずつ追加した文字列「ＴＡＴＡＣＣＣＴ」が被照合文字列となる。ヒトゲノムの解析の場合、参照文字列について、一致文字列の先頭及び末尾に追加する塩基の長さは、それぞれ、例えば５０塩基程度であり得る。
　また、上位コンピュータ１０は、一致文字列の末尾にクエリ文字列における対応する所定長の文字列を追加することにより、照合文字列を作成する。ヒトゲノムの解析の場合、クエリ文字列について、一致文字列の末尾に追加する文字列の長さは、例えば５０塩基程度である。
　なお、本開示では、参照文字列について、一致文字列の先頭及び末尾に所定長の文字列が追加されるものとしたが、これに限られず、例えば、先頭又は末尾の一方にのみ所定長の文字列が追加されても良い。また、クエリ文字列について、一致文字列の先頭及び末尾に、それぞれ、所定長の文字列を追加するようにしても良いし、或いは、文字列を追加せずに、一致文字列そのものを照合文字列として扱っても良い。
　次に、上位コンピュータ１０及び下位コンピュータ２０は、参照文字列に基づく被照合文字列とクエリ文字列に基づく照合文字列とからなる文字列ペアに基づいて少なくとも１つの近似文字列を導出する（Ｓ２０４）。近似文字列の導出には、動的計画法を用いた所定のアラインメントが適用される。アラインメントは、２つの配列（文字列）をその要素どうしで置換、挿入及び欠損を許容しつつ比較して、定義されたスコア／ペナルティに従って変異度（類似度）を算出する手法である。アラインメントを実現するアルゴリズムとしては、Ｓｍｉｔｈ－ＷａｔｅｒｍａｎアルゴリズムやＮｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムが知られている。また、動的計画法とは、ある段階で得られた最適解に基づいて次の段階の最適解を算出する手法である。つまり、変異度の算出では、動的計画法に従って、配列状のアラインメント表の各要素に対してスコアが算出され、最大スコアを持つ要素が決定され、これにより、少なくとも１つ以上の近似文字列が導出される。
　以上のように、本実施形態の近似文字列照合では、参照文字列に基づく階層的インデックスが作成された後、与えられたクエリ文字列に従って、該階層的インデックスを探索することにより一致文字列（及びその長さ）が同定され、同定された一致文字列に基づく被照合文字列と照合文字列とからなる文字列ペアに対して近似文字列照合がなされることにより、近似文字列が導出される。
　図３は、本発明の一実施形態に係るコンピュータシステムによるインデックス作成処理の一例を説明するフローチャートである。すなわち、図３は、図２に示した階層的インデックスの作成処理（Ｓ２０１）の詳細を示している。また、図４は、階層的インデックスを作成するための参照文字列の一例を示す図であり、図４Ａ～４Ｃは、参照文字列に基づく階層的インデックスの作成過程におけるデータ配列構造の一例を示す図である。更に、図５は、階層的インデックスの作成過程におけるキーの出現開始位置及び出現回数を示すテーブル構造の一例を示す図であり、図６は、階層的インデックスを説明するための図である。
　まず、図３に示すように、上位コンピュータ１０は、階層的インデックスを作成するための参照文字列を受信する（Ｓ３０１）。例えば、ヒトゲノム参照配列であれば、上位コンピュータ１０は、データベース３０にアクセスし、格納されているヒトゲノム参照配列を読み出す。以下では、理解容易のため、４種類の文字「Ａ」、「Ｃ」、「Ｇ」、及び「Ｔ」から構成される参照文字列「ＣＣＧＡＴＣＴＧＴＡＴＡＣＣＣＴＡＣＧＡ」を例にして説明する（図４参照）。
　上位コンピュータ１０は、受信した参照文字列について、所定のキー長に従って各部分文字列（すなわち、キー）を切り出して、キー配列を作成する（Ｓ３０２）。例えば、キー長が「２」である場合、キー配列は、図４Ａ（ａ）のようになる。同図中、左端の番号は、配列番号である。また、参照文字列の末端である１９番目の「Ａ」で始まるキーは、説明の簡略化のため、ここでは省略している。
　次に、上位コンピュータ１０は、作成したキー配列の各キーについて、所定のハッシュ関数を用いてハッシュ値を算出し、これを該キーに割り当てる（Ｓ３０３）。これにより、キー配列は、図４Ａ（ｂ）のようになる。本開示におけるハッシュ関数は、４種類の文字「Ａ」、「Ｃ」、「Ｇ」、及び「Ｔ」にそれぞれ割り当てた「０」～「３」の数値により、４進数で表現した値を出力する関数として定義されるが、これに限られない。例えば、１番目のキー「ＣＧ」については、ハッシュ値は、ｈ（ＣＧ）＝１×４＋２＝６となり、また、１１番目のキー「ＡＣ」については、ハッシュ値は、ｈ（ＡＣ）＝０×４＋１＝１となる。
　次に、上位コンピュータ１０は、キー配列の各キーを、割り当てたハッシュ値に従って、例えば昇順にソートする（Ｓ３０４）。これにより、キー配列は、図４Ａ（ｃ）のようになる。本例では、ソート後の各キー配列は、ソート前の配列番号を含み得る。例えば、図４Ａ（ｃ）中、ソート後のキー配列における０番目のキー「ＡＣ」は、ソート前の（元の）配列番号「１１」を保持し、また、ソート後の１番目のキー「ＡＣ」は、ソート前の配列番号「１６」を保持している。
　なお、上記の例では、切り出された各キーについて、算出したハッシュ値を割り当てて、ソートするものとしているが、これに限られない。例えば、参照文字列から切り出されるキーに拘わらず、参照文字列に現れる全ての文字の組み合わせに基づいてハッシュ値を算出して割り当てたキー配列を用意し、切り出されるキーに対応する出現開始位置を割り当てても良い。すなわち、４種類の文字「Ａ」、「Ｃ」、「Ｇ」、及び「Ｔ」から構成される参照文字列において、例えば、キーの長さが２文字であれば、４２個の要素を有するキー配列がまず作成され、更に、ハッシュ値がそれぞれ割り当てられる。続いて、切り出されたキーは、参照文字列における出現開始位置とともに、作成されたキー配列における対応する要素（同じキーの要素）に割り当てられることにより、図４Ａ（ｃ）に示すようなキー配列が得られる。
　次に、上位コンピュータ１０は、現在のキー配列における各キーの出現開始位置及び出現回数を同定する（Ｓ３０５）。図５（ａ）は、図４Ａ（ｃ）に示されるキー配列における各キーの出現開始位置及び出現回数を示している。例えば、キー「ＡＣ」は、キー配列において、出現開始位置「０」（配列番号「０」）を基点にして２回出現することが示されている。また、キー「ＣＣ」は、出現開始位置「４」を基点にして３回出現することが示されている。なお、ここでは、各文字どうしの全ての組み合わせからなるキーのパターン（すなわち、１６パターン）に対するその出現開始位置及び出現回数が示されており、例示したキー配列に含まれていない例えばキー「ＡＡ」については、出現開始位置「－」及び出現回数「０」のように示されている。
　次に、上位コンピュータ１０は、各キーについて、その出現回数が所定の許容値を超えているか否かを判断する（Ｓ３０６）。本開示において、所定の許容値は、階層的インデックスにおいて同じキーが重複して存在し得ることを許容する値である。本例では、所定の許容値は「１」としている。つまり、所定の許容値が「１」であれば、階層的インデックスにおいて各キーは唯一の存在となる。また、所定の許容値が大きいほど、階層的インデックスにおいて重複したキーが存在する可能性が高くなる一方、階層的インデックスの作成は高速化される。上位コンピュータ１０は、出現回数が所定の許容値を超えているキーがあると判断する場合（Ｓ３０６のＹｅｓ）、そのキーに対して追加キーを追加する（Ｓ３０８）。
　追加キーは、元の文字列における該キーに続く１以上の文字である。本例では、追加キーは１文字としている。追加キーの追加により得られる部分文字列は、新たなキーとみなされる。以下では、追加キーが追加された新たなキーを元のキーと区別するために「拡張キー」と称し、その配列を拡張キー配列と称する場合がある。追加キーが追加されることにより、各キーどうしが異なるものとして識別されることになる。図４Ｂ（ａ）は、元のキーに１個の追加キーが追加された拡張キーからなる拡張キー配列の一例を示している。なお、配列番号「１２」の拡張キー「ＧＡ」については、元の文字列において「Ａ」に続く文字がないため、終端文字列として例えば「＄」を割り当てている。また、配列番号「１３」、「１７」及び「１８」のキーについては、その出現回数が１回であるため、追加キーは追加されない。
　次に、上位コンピュータ１０は、各キー（すなわち、拡張キー）を、図４Ｂ（ｂ）に示すように、該追加キーに従って例えば昇順に更にソートする（Ｓ３０８）。この場合、元のキーのソート順が優先される。同図中、例えば、配列番号「２」及び「３」のキー配列「ＡＴ」については、追加キーによるソートで、その順序が入れ替わっていることがわかる。続いて、上位コンピュータ１０は、処理Ｓ３０６に戻り、全てのキーの出現回数が所定の許容値を超えなくなるまで上記の処理を繰り返す。
　すなわち、上位コンピュータ１０は、キー配列における各キーの出現開始位置及び出現回数を同定し（Ｓ３０５）、出現回数が所定の許容値を超えているキーがないと判断する場合（Ｓ３０６のＮｏ）、所望の階層的インデックスが作成されたため、処理を終了する。
　図５は、図４に示される参照文字列に基づく階層的インデックスの作成過程におけるキーの出現開始位置及び出現回数を説明するための図である。例えば、図５（ａ）において出現回数が２以上であるキーには、追加キーが追加され（図４Ｂ（ａ））、各キー（拡張キー）は、追加キーに従ってソートされる（図４Ｂ（ｂ））。これにより、図５（ｂ）に示されるように、現在のキー配列における各キーの出現開始位置及び出現回数が同定される。同図に示す例では、拡張キー「ＣＧＡ」及び「ＴＡＣ」の出現回数が２となっている。したがって、これらの拡張キーのそれぞれについて、同様に、追加キーが追加されソートされる（図４Ｃ（ａ）及び（ｂ））。なお、元の配列番号「１７」（ハッシュ値でソート後の配列番号「８」）のキー「ＣＧ」については、キー「Ａ」に続く文字がないため、終端文字列として例えば「＄」が割り当てられている。これにより、図５（ｃ）に示されるように、該キーの出現開始位置及び出現回数が同定される。以上により、全ての拡張キーは、その出現回数が「１」となったため、インデックスの作成処理が終了する。
　一例として、キー「ＴＡ」について考える。キー「ＴＡ」は、配列番号「１４」～「１６」にあることから（図４Ａ（ｃ）参照）、図６（ａ）に示すように、その出現開始位置は「１４」、出現回数は「３」となる。次に、キー「ＴＡ」に追加キーが追加されソートされることより（図４Ｂ（ｂ）参照）、図６（ｂ）に示すように、追加キー「Ｃ」を含む拡張キー「ＴＡＣ」については、その出現開始位置は「１４」、出現回数は「２」となる一方、追加キー「Ｔ」を含む拡張キー「ＴＡＴ」については、その出現開始位置は「１６」、出現回数は「１」となる。したがって、出現回数が「２」のキー「ＴＡＣ」について、更に追加キー「Ｃ」及び「Ｇ」がそれぞれ追加され、これにより、図６（ｃ）に示すように、キー「ＴＡＣＧ」の出現回数は「１」となる。このようにして、拡張キー配列は、階層的なツリー構造として把握される。
　以上のようにして、上位コンピュータ１０は、例えばヒトゲノム参照配列に基づいて、階層的インデックスを作成する。このような階層的インデックスは、ハッシュ値に従ってソートされているため、特定のキー（部分文字列）に関連する階層的インデックスの部分的なデータ配列構造は、メインメモリにおける特定のアドレス領域に集約的に展開され（データのシーケンシャル化）、これにより、メインメモリに対するランダムアクセスの回数を大幅に減らすことができるようになる。
　なお、上記では、簡単化のため、極めて短い文字列を例にして説明したが、例えば、ヒトゲノムを扱う場合には、部分文字列のキー長を１０～２０程度、追加キーのキー長を２～８、所定の許容値を５～４０とすることが好ましく、部分文字列のキー長を１０～１５程度、追加キーのキー長を２～４、所定の許容値を１０～２０とすることがより好ましい。
　図７は、本発明の一実施形態に係るコンピュータシステムによるクエリ文字列に基づく参照文字列の探索処理の一例を説明するフローチャートである。すなわち、図７は、図２に示したクエリ文字列に基づく参照文字列の探索処理（Ｓ２０２）の詳細を示している。かかる探索処理により、クエリ文字列における所定のキーが参照文字列にマッピングされ、参照文字列における一致文字列及びその出現開始位置が同定される。なお、以下では、上位コンピュータ１０による探索処理が説明されるが、並列的に動作する下位コンピュータ２０による探索処理も同様である。
　同図に示すように、上位コンピュータ１０は、データベース３０から階層的インデックスを読み出して、メモリ上に展開し、記憶する（Ｓ７０１）。上位コンピュータ１０は、高速化の観点から、階層的インデックスを構成するデータをメインメモリ上に連続的に展開し、記憶する。ここで、連続的に展開とは、データが同一バンクにおける連続的なメモリアドレスに配置されることを含む。
　次に、上位コンピュータ１０は、参照文字列に対してマッピングを行うためのクエリ文字列を受信する（Ｓ７０２）。例えば、参照文字列がヒトゲノム参照配列であれば、クエリ文字列は、シークエンサから読み出される塩基対のデータ片である。以下では、理解容易のため、クエリ文字列は「ＴＡＣＣ」であるものとして説明する。
　次に、上位コンピュータ１０は、受信したクエリ文字列について、所定のキー長に従って各キーを連続的に切り出す（Ｓ７０３）。本例では、所定のキー長の初期値は「２」であるものとする。したがって、切り出される各キーは、「ＴＡ」、「ＡＣ」、「ＣＣ」、及び「Ｃ＄」となる。また、本例では、各キーのサンプリング間隔の初期値は１であるものとする。キーのサンプリング間隔とは、該キーに従って階層的インデックスを順に参照する開始位置の間隔である。なお、ヒトゲノムの解析であれば、サンプリング間隔の初期値は、例えば３～５程度であり得る。後述するように、サンプリング間隔は、キーの長さに応じて伸長される。
　続いて、上位コンピュータ１０は、切り出した各キーについて、所定のハッシュ関数を用いてハッシュ値を算出する（Ｓ７０４）。上述したように、ハッシュ関数は、４種類の文字「Ａ」、「Ｃ」、「Ｇ」、及び「Ｔ」にそれぞれ割り当てた「０」～「３」の数値により、４進数で表現した値を出力する関数として定義される。例えば、キー「ＴＡ」について、ハッシュ値は、ｈ（ＴＡ）＝３×４＋０＝１２となる。
　次に、上位コンピュータ１０は、算出したハッシュ値に従って、階層的インデックスを参照し（Ｓ７０５）、参照文字列における各キーの出現開始位置及び出現回数を同定する（Ｓ７０６）。例えば、キー「ＴＡ」であれば、ハッシュ値に従って、階層的インデックスにおける出現開始位置は１４、出現回数は３であることが同定される（図６（ａ）参照）。
　次に、上位コンピュータ１０は、各キーについて、その出現回数が所定の回数しきい値を超えているか否かを判断する（Ｓ７０７）。上述したように、本例では、所定の許容値は１としている。上位コンピュータ１０は、キーの出現回数が所定の許容値を超えていると判断する場合（Ｓ７０７のＹｅｓ）、続いて、現在のキー長が所定の上限値を超えているか否かを判断する（Ｓ７０８）。所定の上限値は、例えば、ヒトゲノムの解析であれば、２０程度であり得るが、これに限られない。
　上位コンピュータ１０は、各キーについて、現在のキー長が所定の上限値を超えていないと判断する場合（Ｓ７０８のＮｏ）、該キーに対して追加キーを追加することにより、拡張キーを作成し（Ｓ７０９）、Ｓ７０６の処理に戻る。
　例えば、上位コンピュータ１０は、キー「ＴＡ」に対して追加キー「Ｃ」を追加して、その出現回数を調べ（図６（ｂ）参照）、更に、追加キーが追加されたキー（拡張キー）に対して追加キー「Ｃ」を追加して、その出現回数を調べ（図６（ｃ）参照）、出現回数が１になるまで繰り返す。
　一方、上位コンピュータ１０は、各キーについて、現在のキー長が所定の上限値を超えていると判断する場合（Ｓ７０８のＹｅｓ）、該キーのサンプリング間隔を伸長（大きく）し（Ｓ７１０）、Ｓ７０５の処理に戻る。なお、サンプリング間隔の伸長は、キーについて、所定回数（例えば一回）だけ行われるようにしても良い。
　また、上位コンピュータ１０は、現在のキー（すなわち、拡張キー）について、その出現回数が所定の許容値を超えていないと判断する場合（Ｓ７０７のＮｏ）、同定された参照文字列におけるキー（一致文字列）及びその出現開始位置を出力する（Ｓ７１１）。
　以上のようにして、上位コンピュータ１０は、階層的インデックスを用いて、参照文字列の中からクエリ文字列の少なくとも一部を探し出すことができる。とりわけ、本実施形態によれば、階層的インデックスの各キーはハッシュ値に従ってソートされているため、クエリ文字列に従った探索は、階層的インデックスにおける部分的・連続的なデータ配列構造に対して行われことになり、メインメモリに対するランダムアクセスの回数を大幅に減らすことができるようになる。
　また、階層的インデックスの探索において、現在のキー（拡張キー）の長さが一定長以上になった場合に出現開始位置のサンプリング間隔を伸長するので、探索の回数が効率的に削減され、探索の高速化を図ることができる。一方で、サンプリング間隔の伸長は、本来同定されるべきキーを見逃す確率が高くなる可能性があるが、キーの長さが一定長以上になった場合にサンプリング間隔を長くしているので、このような見逃しの発生を効果的に抑制している。
　次に、動的計画法を用いたアラインメントについて説明する。すなわち、上位コンピュータ１０は、参照文字列における同定した出現開始位置近傍の文字列（被照合文字列）に対してクエリ文字列（照合文字列）がどれくらい変異しているか（変異度）を推定するために、動的計画法を用いたアラインメント処理を実行する。
　アラインメントとは、２つの配列（文字列）をその要素どうしで置換、挿入及び欠損を許容しつつ比較して、定義されたスコア／ペナルティに従って変異度（類似度）を算出する手法である。例えば、文字列Ｘ（すなわち、参照文字列に基づく被照合文字列）と文字列Ｙ（すなわち、クエリ文字列に基づく照合文字列）との比較において、文字列Ｙの一部の文字が置換され、挿入され、又は欠損する場合に、文字列Ｙは文字列Ｘに一致しないと判断される。つまり、文字列Ｘ及びＹのそれぞれにおける各位置の文字どうしの関係は、一致する場合（一致）、一致しない場合（不一致）、及び一方の文字が存在しない場合（ギャップ）のいずれかであるといえる。ここで、一致する場合のスコアを例えば「＋２」、不一致の場合のスコアを例えば「－１」、及びギャップの場合のスコアを例えば「－２」と定義する。そして、文字列Ｘ及び文字列Ｙの各要素（文字）どうしを比較して、これらのスコアを用いて各要素の変異度が算出される。本開示では、変異度の算出のために、グローバルアラインメントの一例であるＮｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムをベースにした改良アラインメントアルゴリズムが用いられる。以下、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムの基本的な考え方を説明し、更に、本発明に適用される改良アラインメントアルゴリムを説明する。
　例えば、文字列Ｘ＝ｘ１ｘ２.. ｘiと文字列Ｙ＝ｙ１ｙ２.. ｙｊとの比較において、文字ｘiと文字ｙｊとの変異度Ｆ（ｉ，ｊ）は以下のように定義される。
Figure JPOXMLDOC01-appb-M000001
 　ここで、ｍａｘは与えられた式の値の中から最大値を出力する関数、ｓはスコア関数（一致：ｓ＝２、不一致：ｓ＝－１、ギャップ：ｓ＝－２）、ｄはギャップによるペナルティ（ｄ＝２）である。
　以下では、理解容易のため、被照合文字列「ＧＣＣＴＣＧＣＴ」と照合文字列「ＧＣＣＡＴＴＣＡ」との間での動的計画法を用いたアラインメントを説明する。
　図８Ａは、本例における比較対象の文字列どうしを配列したアラインメント表である。表中、「φ」は空文字であり、Ｆ（０，０）＝０とする。また、式１において、変異度Ｆ（ｉ，ｊ）の引数が負の値になる場合、範囲外であるため、計算は省略される（出力をＮｕｌｌとする。）。アラインメント表は、ある種のデータ構造としてメモリ上に展開され、プロセッサの利用に供される。
　まず、ｊ＝０の場合において、Ｆ（１，０）については、式１より、ｍａｘ関数内のそれぞれは、
　　Ｆ（０，－１）＋ｓ（ｘ１，ｙ０）＝Ｎｕｌｌ
　　Ｆ（０，０）－ｄ＝０－２＝－２
　　Ｆ（１，－１）－ｄ＝Ｎｕｌｌ
であるから、
　　Ｆ（１，０）＝－２
となる。
　次に、Ｆ（２，０）は、
　　Ｆ（２，０）＝Ｆ（１，０）－ｄ＝－４
となる。同様にして、
　　Ｆ（ｎ，０）＝－ｎｄ
　　Ｆ（０，ｎ）＝－ｎｄ
となるため、アラインメント表は図８Ｂに示すようになる。
　次に、ｊ＝１の場合においても、同様に算出される。Ｆ（１，１）については、ｍａｘ関数内のそれぞれは、
　　Ｆ（０，０）＋ｓ（ｘ１，ｙ１）＝０＋２＝２
　　Ｆ（０，１）－ｄ＝－２－２＝－４
　　Ｆ（１，０）－ｄ＝－２－２＝－４
であり、これにより、
　　Ｆ（１，１）＝２
となり、アラインメント表は図８Ｃに示すようになる。
　以上の計算を同様に繰り返すことにより、図８Ｄに示すようなアラインメント表が作成されることになる。作成されたアラインメント表において、要素位置（６，７）の変異度Ｆ（６，７）が最大値「７」を有している。したがって、同図に示すように、要素位置（６，７）から要素位置（１，１）までバックトラックがなされる。表中、右下斜め方向への矢印は文字の一致を示し、右方向への矢印は欠損を示し、下方向への矢印は挿入を示している。これにより、近似文字列「ＧＣＣ＜Ａ＞Ｔ［Ｔ］Ｇ」が導出されることになる。ただし、記号「＜　＞」は、文字間への挿入を表し、記号「［　］」は置換を表すものとする。つまり、参照文字列「ＧＣＣＴＣＧ」であるところ、クエリ文字列は、参照文字列の「Ｃ」と「Ｔ」の間に「Ａ」が挿入され、参照文字列の「ＴＣＧ」の「Ｃ」が「Ｔ」に置換されていることがわかる。
　以上のように、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムに従って、アラインメント表における変異度が最大値である要素位置を特定することにより、そこから近似文字列を導出することができる。しかしながら、Ｎｅｅｄｌｅｍａｎ－Ｗｕｎｓｃｈアルゴリズムではアラインメント表の全ての要素の変異度を算出するため、計算量が膨大となり、時間がかかっていた。そこで、本開示では、以下のような改良アラインメントアルゴリズムを提案し、これにより、計算量を削減し、処理の高速化を図っている。
　すなわち、本発明に適用される改良アラインメントアルゴリズムは、概略的には、アラインメント表における対角線上に位置する要素を中心とする所定の幅を有する計算領域を定め、該計算領域内の要素についてのみ変異度を算出することによりその最大値を決定し、該最大値が所定の条件を満たす場合に、該最大値に基づく要素から近似文字列を導出することを含む。最大値が所定の条件を満たさない場合には、計算領域が拡大され、同様に、変異度が算出されることによりその最大値を決定し、該最大値が所定の条件を満たすまで繰り返される。
　図９は、本発明の一実施形態に係るコンピュータシステムによる動的計画法を用いたアラインメント処理の一例を説明するフローチャートである。すなわち、図９は、図２に示したアラインメント処理の（Ｓ２０４）の詳細を示している。なお、以下では、上位コンピュータ１０による一のクエリ文字列に基づくアラインメント処理が説明されるが、並列的に動作する下位コンピュータ２０による他のクエリ文字列に基づくアラインメント処理も同様である。
　同図に示すように、上位コンピュータ１０は、被照合文字列と照合文字列とからなる文字列ペアに基づいてアラインメント表を作成する（Ｓ１００１）。上述したように、アラインメント表は、ある種のデータ構造としてメモリ上に展開される。
　次に、上位コンピュータ１０は、計算領域の幅ｍ（ただし、ｍは正数）を初期値に設定する（Ｓ１００２）。幅ｍの初期値は、例えば、マッピングにより得られた一致文字列の長さであり得るが、これに限られない。また、幅ｍは、アラインメント表の対角線上の要素位置を中心とすることから、奇数の値に設定されるが、これに限られるものではない。これにより、アラインメント表の対角線上の要素位置を中心とする幅ｍの計算領域が決定される。
　続いて、上位コンピュータ１０は、計算領域の境界を画定する各要素に所定のダミー値を設定する（Ｓ１００３）。ダミー値は、変異度Ｆの値として十分に小さい値が選択される。例えば、ダミー値は、初期値の幅ｍの例えば２～３倍程度の負の値であり、任意に設定することができる。
　次に、上位コンピュータ１０は、計算領域内の各要素について、式１に従って変異度Fを算出する（Ｓ１００４）。続いて、上位コンピュータ１０は、計算領域において最大値を有する最大変異度Ｆｍａｘを決定し、その要素の位置を特定する（Ｓ１００５）。なお、最大変異度Ｆｍａｘを持つ要素は、１つであるとは限らない。
　また、上位コンピュータ１０は、アライメント表の行又は列における変異度Ｆの下限値ＦＬｏｗを算出する（Ｓ１００６）。下限値ＦＬｏｗは、該行又は列において、参照文字列とクエリ文字列とを比較して、ｍ個の連続したギャップがあり、それ以外の部分は完全に又は実質的に一致したと仮定した場合の変異度Ｆの値である。すなわち、下限値ＦＬｏｗは、
　　ＦＬｏｗ＝（文字列の長さ－ｍ）×ｓ　　…式２
　ただし、ｓ＝２である。
で算出される。
　次に、上位コンピュータ１０は、最大変異度Ｆｍａｘと下限値ＦＬｏｗとを比較して、最大変異度Ｆｍａｘが下限値ＦＬｏｗを超えているか否かを判断する（Ｓ１００７）。上位コンピュータ１０は、最大変異度Ｆｍａｘが下限値ＦＬｏｗを超えていないと判断する場合（Ｓ１００７のＮｏ）、幅ｍを所定の大きさδだけ拡幅する（Ｓ１００８）。例えば、δは、ｍ＋１とする（ただし、ｍは文字列の文字数を超えないものとする。）。
　上位コンピュータ１０は、拡幅された幅ｍの計算領域に対して、同様に処理を行い、最大値Ｆｍａｘが下限値ＦＬｏｗを超えるまで、上記処理を繰り返す。
　上位コンピュータ１０は、最大値Ｆｍａｘが下限値ＦＬｏｗを超えていると判断する場合（Ｓ１００７のＹｅｓ）、該最大値を持つ要素の位置からバックトラックを行って、近似文字列を決定する（Ｓ１００９）。そして、上位コンピュータ１０は、決定した近似文字列を出力する（Ｓ１０１０）。
　例えば、被照合文字列「ＧＧＧＡＴＣＣＧＡＴＡＡＴＣＧＧＴＣＣＣＣＴＡＧＧ」（２４文字）に対して照合文字列「ＧＧＧＣＡＴＴＣＡＡＣＡＴＡＡＧＴＣＧＧＣＣＴＧ」（２４文字）との間での、本発明に係るアラインメント法による変異度の算出例を説明する。なお、変異度Ｆの算出に式１を用いる点は、上述した例と同様である。なお、被照合文字列の長さと照合文字列の長さとは一致する必要はなく、典型的には、被照合文字列の長さの方が照合文字列の長さよりも長い。
　まず、上位コンピュータ１０は、比較対象の文字列どうしを配列したアラインメント表を用意し、幅ｍの初期値を設定する。本例では、幅ｍ（０）の初期値は７であるものとする。また、上位コンピュータ１０は、幅ｍに従って規定される計算領域の境界部分の要素にダミー値を設定する。本例では、ダミー値は－２０であるものとする。図１０Ａは、ダミー値が設定されたアラインメント表を示している。表中、ハッチングが描かれている要素が幅ｍ（０）＝７での計算領域Ｒ（０）である。
　上位コンピュータ１０は、計算領域Ｒ（０）内の各要素の変異度Ｆを式１に従って計算する。図１０Ｂは、計算領域内の各要素の変異度Ｆが算出された状態を示している。上位コンピュータ１０は、算出された変異度Ｆの中から、最大変異度Ｆｍａｘを決定する。本例では、最大変異度Ｆｍａｘは１７である。図中、最大変異度Ｆｍａｘが１７である要素にはハッチングが示されている。
　続いて、上位コンピュータ１０は、アラインメント表の行又は列における下限値ＦＬｏｗを算出する。本例では、下限値ＦＬｏｗは、
　　ＦＬｏｗ＝（２４－ｍ）×ｓ
　　　　　＝１７×２
　　　　　＝３４
となる。
　続いて、上位コンピュータ１０は、最大変異度Ｆｍａｘと下限値ＦＬｏｗとを比較し、これにより、最大変異度Ｆｍａｘが下限値ＦＬｏｗを超えていないと判断するため、幅ｍをδだけ拡幅する。本例では、拡幅された幅ｍ（１）を１５に拡幅する。また、拡幅された幅ｍ（１）の計算領域をＲ（１）とする。
　上位コンピュータ１０は、同様にして、計算領域Ｒ（１）内の各要素の変異度Ｆを式１に従って算出する。図１０Ｃは、計算領域内の各要素の変異度Ｆが算出された状態を示している。図中、計算領域Ｒ（１）に対して拡幅により追加された領域にハッチングが描かれている。これにより、最大変異度Ｆｍａｘは１９となる。また、このときの下限値ＦＬｏｗは１８となる。
　したがって、上位コンピュータ１０は、最大変異度Ｆｍａｘが下限値ＦＬｏｗを超えていると判断するため、最大変異度Ｆｍａｘ＝１９である要素からバックトラックし、これにより得られるパスに従って近似文字列を特定する。
　すなわち、図１０Ｃに示す例では、上位コンピュータ１０は、最大変異度Ｆｍａｘ＝１９である要素（１８，２２）を起点（現在の要素位置）として、要素（０，０）方向に向けて隣接する要素のうち変異度Ｆの値が最も大きい要素を同定し、そこに遷移する。したがって、変異度Ｆ＝１７を持つ要素（１７，２１）が現在の要素位置となる。このような遷移を要素（０，０）まで繰り返すことにより、最終的に、近似文字列が導出される。なお、バックトラックにより得られるパスは、１つとは限られず、複数である場合がある。
　より具体的には、図１０Ｃに示すアラインメント表において、バックトラックにより得られるパスは６通りあり、各パスに従う近似文字列は、以下のとおりとなる。
（ａ）第１のパス：ＧＧＧ＜Ｃ＞Ａ＜Ｔ＞ＴＣ＜ＡＡ＞Ｃ－ＡＴＡＡ＜Ｇ＞ＴＣＧ［Ｇ］ＣＣ
（ｂ）第２のパス：ＧＧＧ＜Ｃ＞Ａ＜Ｔ＞ＴＣ＜Ａ＞［Ａ］［Ｃ］ＡＴＡＡ＜Ｇ＞ＴＣＧ［Ｇ］ＣＣ
（ｃ）第３のパス：ＧＧＧ＜Ｃ＞ＡＴ［Ｔ］Ｃ＜Ａ＞［Ａ］［Ｃ］ＡＴＡＡ＜Ｇ＞ＴＣＧ［Ｇ］ＣＣ
（ｄ）第４のパス：ＧＧＧ＜Ｃ＞ＡＴ［Ｔ］Ｃ［Ａ］＜ＡＣ＞ＡＴＡＡ＜Ｇ＞ＴＣＧ［Ｇ］ＣＣ
（ｅ）第５のパス：ＧＧＧ＜Ｃ＞Ａ＜Ｔ＞ＴＣ＜ＡＡ＞Ｃ－ＡＴＡＡ＜Ｇ＞ＴＣＧ［Ｇ］ＣＣ
（ｆ）第６のパス：ＧＧＧ＜Ｃ＞ＡＴ＜Ｔ＞Ｃ＜Ａ＞［Ａ］［Ｃ］ＡＴＡＡ＜Ｇ＞ＴＣＧ［Ｇ］ＣＣ
　ただし、記号「＜　＞」は、文字間への挿入を表し、記号「［　］」は文字の置換を表し、記号「－」は文字の欠損を表すものとする。
　なお、理解容易のため、上記の各パスに従う近似文字列を参照文字列との対比において図１１Ａ及び１１Ｂに示している。
　このように、改良されたアラインメントアルゴリズムでは、アラインメント表の全ての要素について変異度を算出するのではなく、所定の幅を有する計算領域を定め、必要に応じた範囲で計算領域を拡大させながら変異度を算出していくので、計算量を削減し、これにより、処理の高速化を図ることができるようになる。
　以上のように、本実施形態によれば、参照文字列に基づく階層的インデックスが作成された後、与えられたクエリ文字列に従って、該階層的インデックスを探索することにより一致文字列（及びその長さ）が同定され、同定された一致文字列に基づく被照合文字列と照合文字列とからなる文字列ペアに対して近似文字列照合がなされることにより、近似文字列が導出される。
　上記各実施形態は、本発明を説明するための例示であり、本発明をこれらの実施形態にのみ限定する趣旨ではない。本発明は、その要旨を逸脱しない限り、さまざまな形態で実施することができる。
　例えば、本明細書に開示される方法においては、その結果に矛盾が生じない限り、ステップ、動作又は機能を並行して又は異なる順に実施しても良い。説明されたステップ、動作及び機能は、単なる例として提供されており、ステップ、動作及び機能のうちのいくつかは、発明の要旨を逸脱しない範囲で、省略でき、また、互いに結合させることで一つのものとしてもよく、また、他のステップ、動作又は機能を追加してもよい。
　また、本明細書では、さまざまな実施形態が開示されているが、一の実施形態における特定のフィーチャ（技術的事項）を、適宜改良しながら、他の実施形態に追加し、又は該他の実施形態における特定のフィーチャと置換することができ、そのような形態も本発明の要旨に含まれる。
１…コンピュータシステム
１０…上位コンピュータ
２０…下位コンピュータ
３０…データベース

Claims (21)

  1. 　コンピューティングデバイスに、クエリ文字列に基づいて参照文字列における近似文字列を検索するための方法を実行させるためのコンピュータプログラムであって、
    　前記方法は、
    　前記参照文字列に基づいて階層的インデックスを作成することと、
    　前記クエリ文字列の少なくとも一部と一致する前記参照文字列における部分文字列を同定するために、前記階層的インデックスを参照して、前記参照文字列に対する前記クエリ文字列のマッピングを行うことと、
    　前記のマッピングにより同定される前記部分文字列に基づいて作成される、前記参照文字列に基づく被照合文字列と前記クエリ文字列に基づく照合文字列とに基づいて、前記被照合文字列に近似する文字列を前記近似文字列として導出することと、を含み、
    　前記階層的インデックスを作成することは、
    　前記参照文字列から所定長の各第１のキーを切り出すことと、
    　切り出された前記各第１のキーについて、所定のハッシュ関数により該第１のキーに基づいて算出されるハッシュ値を割り当てた第１のキー配列を作成することと、
    　作成された前記第１のキー配列を更新することと、
    　更新された前記第１のキー配列を前記階層的インデックスとして出力することと、を含み、
    　前記第１のキー配列を更新することは、
    　前記第１のキー配列における前記各第１のキーについて、前記参照文字列における該第１のキーの出現回数を同定することと、
    　同定された前記第１のキーの前記出現回数に従って、該第１のキーに対して前記参照文字列における該第１のキーに続く少なくとも１以上の文字からなる第１の追加キーを追加することにより前記第１のキー配列を更新することと、を含む、
    コンピュータプログラム。
  2. 　前記第１のキー配列を作成することは、前記ハッシュ値に従って前記第１のキー配列における前記各第１のキーをソートすることを含む、
    請求項１に記載のコンピュータプログラム。
  3. 　前記第１のキー配列を更新することは、
    　前記同定された出現回数が所定の許容値を超えているか否かを判断することと、
    　前記同定された出現回数が前記所定の許容値を超えていると判断される場合に、前記第１のキーに対して前記参照文字列における該第１のキーに続く少なくとも１以上の文字からなる前記第１の追加キーを追加することにより新たな第１のキーを作成することと、
    　前記新たな第１のキーについて、前記参照文字列における該新たな第１のキーの出現回数を同定することと、を含む、
    請求項１又は２に記載のコンピュータプログラム。
  4. 　前記第１のキー配列を更新することは、前記第１の追加キーに従って前記第１のキー配列における前記新たな第１のキーをソートすることを更に含む、
    請求項１から３のいずれか一項に記載のコンピュータプログラム。
  5. 　前記第１のキー配列を更新することは、前記同定された前記出現回数が所定の許容値を超えていないと判断されるまで、現在の前記第１のキーに新たな前記第１の追加キーを順次に追加することにより新たな前記第１のキーを作成することを含む、
    請求項３又は４に記載のコンピュータプログラム。
  6. 　前記第１のキー配列を前記階層的インデックスとして出力することは、
    　前記同定された前記出現回数が所定の許容値を超えていないと判断される場合に、現在の前記第１のキー配列を前記階層的インデックスとして出力することを含む、
    請求項３から５のいずれか一項に記載のコンピュータプログラム。
  7. 　前記マッピングを行うことは、
    　前記クエリ文字列から所定長の各第２のキーを切り出すことと、
    　前記クエリ文字列から切り出された前記各第２のキーについて、前記所定のハッシュ関数により該第２のキーに基づいて算出されるハッシュ値を割り当てた第２のキー配列を作成することと、
    　前記各第２のキーについて、前記ハッシュ値に従って、所定のサンプリング間隔で、前記階層的インデックスを参照し、該第２のキーの出現開始位置及び出現回数を同定することと、を含む、
    請求項１から６のいずれか一項に記載のコンピュータプログラム。
  8. 　前記第２のキーの前記出現開始位置及び前記出現回数を同定することは、
    　前記第２のキーの前記出現回数が前記所定の許容値を超えているか否かを判断することと、
    　前記第２のキーの前記出現回数が前記所定の許容値を超えていると判断される場合に、前記第２のキーに対して前記クエリ文字列における該第２のキーに続く少なくとも１以上の文字からなる第２の追加キーを追加することにより新たな第２のキーを作成することと、
    　前記第２のキーの前記出現回数が前記所定の許容値を超えていないと判断される場合に、同定された現在の前記第２のキーを前記部分文字列として出力するとともに該第２のキーの前記出現開始位置を出力することと、を含む、
    請求項７に記載のコンピュータプログラム。
  9. 　前記第２のキーの前記出現開始位置及び前記出現回数を同定することは、前記第２のキーの前記同定された前記出現回数が前記所定の許容値を超えていないと判断されるまで、現在の前記第２のキーに新たな前記第２の追加キーを順次に追加して、前記新たな第２のキーを作成することを更に含む、
    請求項８に記載のコンピュータプログラム。
  10. 　前記第２のキーの前記出現回数が前記所定の許容値を超えていると判断される場合に、該第２のキーの前記所定のサンプリング間隔を大きくする、
    請求項８又は９に記載のコンピュータプログラム。
  11. 　前記近似文字列として導出することは、
    　前記マッピングにより同定された前記部分文字列に基づく、前記被照合文字列と前記照合文字列とからなる文字列ペアを受信することと、
    　前記文字列ペアに基づいて少なくとも１つの前記近似文字列を導出するために、所定のアラインメント処理を実行することと、
    　導出された前記少なくとも１つの近似文字列を出力することと、を含む、
    請求項８から１０のいずれか一項に記載のコンピュータプログラム。
  12. 　前記所定のアラインメント処理を実行することは、
    　前記被照合文字列と前記照合文字列とに基づいて所定のアラインメント表を作成することと、
    　前記アラインメント表の対角線上の要素を中心にした幅ｍを有する計算領域を設定することと、
    　設定された前記計算領域における各要素について、変異度を算出することと、
    　算出された前記変異度に基づいて、最大変異度を決定することと、
    　決定された前記最大変異度に基づいて、前記少なくとも１つの近似文字列を導出することを含む、
    請求項１１に記載のコンピュータプログラム。
  13. 　前記所定のアラインメント処理を実行することは、
    　前記最大変異度と所定の下限値とを比較して、前記最大変異度が前記所定の下限値を超えているかを判断することと、
    　前記最大変異度が前記所定の下限値を超えていないと判断される場合に、新たな計算領域を設定するために、前記計算領域の前記幅ｍを拡幅することと、
    　前記最大変異度が前記所定の下限値を超えていると判断される場合に、前記最大変異度を有する要素に基づいて、前記少なくとも１つの近似文字列を導出することと、を含み、
    　前記最大変異度が前記下限値を超えると判断されるまで、前記計算領域を拡幅することにより新たな計算領域を設定して前記変異度を算出することを繰り返す、
    請求項１２に記載のコンピュータプログラム。
  14. 　前記所定の下限値は、所定の要素列にｍ個の連続したギャップがあり、それ以外の部分は一致したと仮定した場合の変異度の値である、
    請求項１３に記載のコンピュータプログラム。
  15. 　前記部分文字列に対して前記参照文字列における対応する所定の文字列を追加することにより前記被照合文字列を作成することと、
    　前記部分文字列に対して前記クエリ文字列における対応する所定の文字列を追加することにより前記照合文字列を作成することと、を更に含む、
    請求項１１から１４のいずれか一項に記載のコンピュータプログラム。
  16. 　コンピューティングデバイスに、クエリ文字列に基づいて参照文字列を探索するための階層的インデックスを作成する方法を実行させるためのコンピュータプログラムであって、
    　前記方法は、
    　前記参照文字列から所定長の各第１のキーを切り出すことと、
    　切り出された前記各第１のキーについて、所定のハッシュ関数により該第１のキーに基づいて算出されるハッシュ値を割り当てた第１のキー配列を作成することと、
    　作成された前記第１のキー配列を更新することと、
    　更新された前記第１のキー配列を前記階層的インデックスとして出力することと、を含み、
    　前記第１のキー配列を更新することは、
    　前記第１のキー配列における前記各第１のキーについて、前記参照文字列における該第１のキーの出現開始位置及び出現回数を同定することと、
    　同定された前記第１のキーの前記出現開始位置及び前記出現回数に従って、該第１のキーに対して前記参照文字列における該第１のキーに続く少なくとも１以上の文字からなる第１の追加キーを追加することにより前記第１のキー配列を更新することと、を含む、
    コンピュータプログラム。
  17. 　コンピューティングデバイスに、参照文字列に対してクエリ文字列のマッピングを行う方法を実行させるためのコンピュータプログラムであって、
    　前記方法は、
    　前記参照文字列から切り出される各第１のキー、及び前記各第１のキーについて、該第１のキーに基づいて所定のハッシュ関数により算出されるハッシュ値が割り当てられた第１のキー配列であって、前記参照文字列における前記第１のキーの出現開始位置及び出現回数に従って前記第１のキーに追加キーが追加された該第１のキー配列からなる階層的インデックスを読み出すことと、
    　前記クエリ文字列から所定のキー長を有する各第２のキーを切り出すことと、
    　前記クエリ文字列から切り出された前記各第２のキーについて、所定のハッシュ関数により該第２のキーに基づいて算出されるハッシュ値を割り当てた第２のキー配列を作成することと、
    　前記各第２のキーについて、前記ハッシュ値に従って、所定のサンプリング間隔で、前記階層的インデックスを参照し、前記第１のキー配列における前記第１のキーとの比較により、該第２のキーの出現開始位置及び出現回数を同定することと、
    　前記同定した出現回数が所定のしきい値を超えているか否かを判断することと、
    　前記同定された前記出現回数が所定の許容値を超えていると判断される場合に、前記第２のキーに対して前記クエリ文字列における該第２のキーに続く少なくとも１以上の文字からなる追加キーを追加することにより新たな第２のキーを作成することと、
    　前記同定された前記出現回数が所定のしきい値を超えていないと判断される場合に、同定された現在の前記第２のキーの出現開始位置及び該キーを出力することと、を含み、
    　前記第２のキーの前記出現開始位置及び前記出現回数を同定することは、前記同定された前記出現回数が所定のしきい値を超えていないと判断されるまで、現在の前記第２のキーに新たな前記追加キーを順次に追加して、前記新たな第２のキーを作成することを含む、
    コンピュータプログラム。
  18. 　前記方法は、前記同定された前記出現回数が所定の許容値を超えていると判断される場合に、前記第２のキーの前記所定のサンプリング間隔を大きくするように構成される、
    請求項１７に記載のコンピュータプログラム。
  19. 　コンピューティングデバイスに、参照文字列における部分文字列とクエリ文字列との間の変異を所定のアラインメント処理により同定する方法を実行させるためのコンピュータプログラムであって、
    　前記方法は、
    　マッピングにより同定される前記部分文字列に基づいて作成される、前記参照文字列に基づく被照合文字列と前記クエリ文字列に基づく照合文字列とからなる文字列ペアを受信することと、
    　前記文字列ペアに基づいて、前記被照合文字列に近似する文字列を少なくとも１つの近似文字列として導出するために、所定のアラインメント処理を実行することと、
    　導出された前記少なくとも１つの近似文字列を出力することと、を含み、
    　前記所定のアラインメント処理を実行することは、
    　前記被照合文字列と前記照合文字列とに基づいて所定のアラインメント表を作成することと、
    　前記アラインメント表の対角線上の要素を中心にした幅ｍを有する計算領域を設定することと、
    　設定された前記計算領域における各要素について、変異度を算出することと、
    　算出された前記変異度に基づいて、最大変異度を決定することと、
    　決定された前記最大変異度に基づいて、前記少なくとも１つの近似文字列を導出することを含む、
    コンピュータプログラム。
  20. 　前記所定のアラインメント処理を実行することは、
    　前記最大変異度と所定の下限値とを比較して、前記最大変異度が前記所定の下限値を超えているかを判断することと、
    　前記最大変異度が前記所定の下限値を超えていないと判断される場合に、新たな計算領域を設定するために、前記計算領域の前記幅ｍを拡幅することと、
    　前記最大変異度が前記所定の下限値を超えていると判断される場合に、前記最大変異度を有する要素に基づいて、前記少なくとも１つの近似文字列を導出することと、を含み、
    　前記最大変異度が前記下限値を超えると判断されるまで、前記計算領域を拡幅することにより新たな計算領域を設定して前記変異度を算出することを繰り返す、
    請求項１９に記載のコンピュータプログラム。
  21. 　前記所定のアラインメント処理を実行することは、所定の要素列にｍ個の連続したギャップがあり、それ以外の部分は一致したと仮定した場合の変異度の値を前記所定の下限値として設定することを更に含む、
    請求項２０に記載のコンピュータプログラム。
PCT/JP2023/008899 2023-03-08 2023-03-08 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム WO2024185097A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/JP2023/008899 WO2024185097A1 (ja) 2023-03-08 2023-03-08 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2023/008899 WO2024185097A1 (ja) 2023-03-08 2023-03-08 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム

Publications (1)

Publication Number Publication Date
WO2024185097A1 true WO2024185097A1 (ja) 2024-09-12

Family

ID=92674215

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2023/008899 WO2024185097A1 (ja) 2023-03-08 2023-03-08 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム

Country Status (1)

Country Link
WO (1) WO2024185097A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN119322917A (zh) * 2024-12-19 2025-01-17 安康大数据运营有限公司 一种应用于基层平台的多源数据管理方法及系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120239706A1 (en) * 2011-03-18 2012-09-20 Los Alamos National Security, Llc Computer-facilitated parallel information alignment and analysis
US20210201163A1 (en) * 2019-12-28 2021-07-01 Intel Corporation Genome Sequence Alignment System and Method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120239706A1 (en) * 2011-03-18 2012-09-20 Los Alamos National Security, Llc Computer-facilitated parallel information alignment and analysis
US20210201163A1 (en) * 2019-12-28 2021-07-01 Intel Corporation Genome Sequence Alignment System and Method

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN119322917A (zh) * 2024-12-19 2025-01-17 安康大数据运营有限公司 一种应用于基层平台的多源数据管理方法及系统

Similar Documents

Publication Publication Date Title
Kim et al. Geniehd: Efficient dna pattern matching accelerator using hyperdimensional computing
JP3672242B2 (ja) パターン検索方法、パターン検索装置、コンピュータプログラム及び記憶媒体
Drew et al. Polymorphic malware detection using sequence classification methods
US10521441B2 (en) System and method for approximate searching very large data
US20140344195A1 (en) System and method for machine learning and classifying data
US11941534B2 (en) Genome sequence alignment system and method
JP2002278761A (ja) 否定項を含む相関ルール抽出方法およびシステム
CN112735528A (zh) 一种基因序列比对方法及系统
US12125559B2 (en) Parallelizable sequence alignment systems and methods
JP2009116559A (ja) 大量配列の一括検索方法及び検索システム
US20150006577A1 (en) Method and system for searching and storing data
Loukides et al. Bidirectional string anchors for improved text indexing and top-$ k $ similarity search
WO2024185097A1 (ja) 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム
US20240404642A1 (en) Genome graph analysis method, device and medium based on in-memory computing
JP7422367B2 (ja) 近似文字列照合方法及び該方法を実現するためのコンピュータプログラム
US10867134B2 (en) Method for generating text string dictionary, method for searching text string dictionary, and system for processing text string dictionary
CN110892401B (zh) 生成用于k个不匹配搜索的过滤器的系统和方法
Luo et al. A collective approach to scholar name disambiguation
Yang et al. Improving regular-expression matching on strings using negative factors
de Armas et al. K-mer Mapping and de Bruijn graphs: the case for velvet fragment assembly
Quedenfeld et al. Variant tolerant read mapping using min-hashing
Lavenier et al. A reconfigurable index FLASH memory tailored to seed-based genomic sequence comparison algorithms
JP2024169825A (ja) インデックス作成用コンピュータプログラム、マッピング用コンピュータプログラム、インデックス作成用コンピュータ装置、マッピング用コンピュータ装置、階層的インデックスのデータ構造。
EP3014482B1 (en) Method and system for searching and storing data
Alatabbi et al. On the repetitive collection indexing problem

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: 23926307

Country of ref document: EP

Kind code of ref document: A1