JP6653628B2 - DNA sequence analyzer, DNA sequence analysis method, and DNA sequence analysis system - Google Patents

Description

  The present invention relates to a DNA sequence analysis device, a DNA sequence analysis method, and a DNA sequence analysis system.

  The DNA sequencing technique (Next Generation Sequencing, NGS) using a next-generation sequencer can determine the genome sequence at a dramatically lower cost than the Sanger method (conventional). For example, the cost of genome sequence determination, which was about $ 100 million in 2001, has fallen to $ 1,245 in 2015 with the advent of next-generation sequencers ("DNA Sequencing Costs, http: //www.genome. gov / sequencingcosts / ”). The next-generation sequencer can not only reduce the cost but also obtain a huge amount of sequence data in a short period of time. For example, next-generation sequencers capable of generating a huge amount of sequence data exceeding one trillion bases at a time have been commercialized. These techniques have made it possible to determine the individual genomic sequence of a large number of subjects.

  When analyzing a human genome sequence using a next-generation sequencer, a sequence obtained by the next-generation sequencer (hereinafter referred to as a “read sequence”) must be compared with a reference genome sequence, which is a standard human genome sequence. Is required. The calculation process for comparing the read sequence with the reference genome sequence, specifying the position of the reference genome sequence corresponding to the read sequence, and clarifying the difference from the reference genome sequence is called “mapping”. The mapping reveals sequences that are unique to the subject's genome and reveals the characteristics of the gene that the subject has. The gene information thus obtained is extremely useful in predicting the disease risk and drug responsiveness of the subject.

  Since mapping is indispensable for processing sequence data in a next-generation sequencer, when the sequence data is published on a database of a public institution, the result mapped to the genome may be published. On the other hand, in recent years, the amount of sequence data handled by the next-generation sequencer has become enormous, which is a problem. Therefore, a technique for recording only the difference from the reference genome sequence using the result of mapping the read sequence to the reference genome sequence and compressing the read sequence has also been developed (Fritz et al., Genome Res. 2011, 21). (5): 734-40).

  In addition, since the amount of read sequences mapped to a reference genome sequence is enormous, the process of mapping read sequences is a major bottleneck in analyzing sequence data on a computer. Therefore, many techniques have been developed for efficiently executing the analysis processing on a computer (for example, Non-Patent Document 1). These techniques employ a technique of directly comparing a given read sequence to a reference genomic sequence.

Li and Durbin 2009, Bioinformatics 2009; 25 (14): 1754-1760.

  As described above, the results of mapping the read sequence to the reference genomic sequence are provided, but the reference genomic sequence is not the only fixed sequence. In fact, there are parts of the human genomic sequence that are difficult to sequence (eg, very long repetitive sequences). For this reason, it is not easy to provide a complete sequence as a reference genome sequence, and the improvement of the reference genome sequence is still ongoing.

  Further, as described above, recently, genome sequences of many individuals can be analyzed by the next-generation sequencer, but the genome sequence of each individual is different from the reference genome sequence. For this reason, it is expected that various reference genome sequences will have to be handled in the future. Along with this, it is expected that a read sequence mapped to one reference genome sequence will need to be re-mapped to another reference genome sequence, and a mechanism for streamlining the mapping process will also be required. Conceivable.

  However, in the above-mentioned Non-Patent Document 1, when a read sequence is already mapped to a certain reference genome sequence (hereinafter referred to as “old genome sequence”), another reference genome sequence ( There is no mechanism provided for efficient mapping to the "new genome sequence."

  As a matter of course, in a region where the two old and new reference genome sequences match, it is conceivable that an existing mapping result is used as it is as a mapping result of the new reference genome sequence. However, when a sequence identical to the sequence matching between the two new and old reference genome sequences is newly inserted into the new genome sequence, the mapping destination of the read sequence is set to the sequence existing in the old genome sequence. It is not possible to uniquely specify whether the sequence should be a newly inserted sequence. That is, there is a problem that the original mapping result cannot be used as it is.

  Accordingly, the present invention provides a technique for efficiently mapping a huge number of read sequences mapped to an old genome sequence to a new genome sequence different from the old genome sequence by using the mapping result to the old genome sequence. I will provide a.

  In order to solve the above-described problems, the present invention employs, for example, a configuration described in the claims. The present specification includes a plurality of means for solving the above-mentioned problems. For example, "(1) a first genomic sequence called an old genome sequence and a second genomic sequence called a new genome sequence" A large number of DNA sequences called a read sequence, a storage device for storing the result of mapping the read sequence to the old genomic sequence, and (2) a search for an arbitrary character string existing in the new genomic sequence. A new genome index constructing unit for constructing an index, and (3) a query which records a combination of mutations with respect to the old genomic sequence present in the read sequence, the mutations being located near a processing position on the old genomic sequence. A data structure called a table, a query table update unit that is constructed based on the result of mapping, and (4) constructed by applying a combination of mutations stored in the query table to the old genome sequence A mapping destination search unit that compares a sequence with the new genome sequence and comprehensively outputs a portion where the number of K (one or more natural numbers) bases or more completely matches the new genome sequence in the constructed sequence. DNA sequence analyzer. "

  According to the present invention, an enormous number of read sequences mapped to an old genome sequence can be efficiently mapped to a new genome sequence different from the old genome sequence by utilizing the result of mapping to the old genome sequence. . Problems, configurations, and effects other than those described above will be apparent from the following description of the embodiments.

FIG. 7 is a diagram for explaining a conventional technique for mapping a read sequence. FIG. 4 is a diagram for explaining an intermediate stage of the lead sequence mapping method according to the embodiment. FIG. 4 is a view for explaining a lead array mapping method according to the embodiment. FIG. 1 is a diagram illustrating a configuration example of a DNA sequence analyzer according to an embodiment. 5 is a flowchart illustrating an outline of a DNA sequence analysis process according to the embodiment. The figure which shows the sequence of the DNA sequence analysis processing which concerns on an Example. 5 is a flowchart illustrating an example of a processing operation performed by a new genome index construction unit. The figure which shows an example of a query table. 9 is a flowchart illustrating an example of a processing operation executed by a query table update unit. 9 is a flowchart for explaining an existing query allele update process (details of S210) executed by the query table update unit. The figure which shows the query table after updating the query table of FIG. 9 is a flowchart illustrating an example of a processing operation performed by a mapping destination search unit. The figure which shows an example of the query table in the area | region where there is no variation in the vicinity. 11 is a flowchart for explaining processing (details of S306) when a matching area is lost, which is executed by the mapping destination search unit. 9 is a flowchart for explaining processing (details of S307) performed when a matching area is extended by a mapping destination search unit. 9 is a flowchart for explaining processing (details of S309) executed by a processing unit used for searching for a new matching area. The figure which shows the example of visualization of a query table.

  Hereinafter, embodiments of the present invention will be described with reference to the drawings. The embodiments of the present invention are not limited to the examples described later, and various modifications are possible within the scope of the technical idea.

(1) Basic Concept In order to efficiently execute the process of remapping a large amount of read sequences mapped to the old genome sequence 108 to the new genome sequence 114, the read sequence 107 used in the conventional technique must be used. It is preferable to avoid the individual remapping process (FIG. 1). This process is inefficient because it does not use known information and requires a very long processing time.

  Therefore, in the embodiment described below, first, a sequence obtained by applying the mutation 401 contained in the read sequence 107 to the old genome sequence 108 (in the example of FIG. 2, the region surrounded by the broken line with respect to the old genome sequence 108) Is compared with the new genome sequence 114 to clarify the corresponding positions in the old genome sequence 108 and the new genome sequence 114. In the example of FIG. 2, of the two similar regions, the region on the right side in the new genome sequence 114 has a higher degree of similarity, and thus this region is the mapping destination. Next, in the embodiment described below, the read sequence 107 mapped to the old genome sequence 108 is collectively moved to the corresponding position in the new genome sequence 114 as shown in FIG. Enhance.

(2) Terms First, definitions of terms used in this specification will be described.
-"Alphabet" is a set of characters constituting a character string to be analyzed. For DNA sequences, the alphabet is {A, C, G, T}.
"Character string" is a string obtained by connecting characters included in the alphabet. For example, "ATTG" is a character string.
-The length of the character string s is expressed as | s |. When s = ATTG, | s | = 4.
-The i-th character of the character string s is described as s [i]. When s = ATTG, s [1] = A. In this specification, the first character is the first character and not the zeroth character.
-A partial character string starting from the beginning of a given character string is called a prefix. “ATTG”, “ATT”, “AT”, “A” are prefixes of the character string “ATTG”.
A substring that ends at the end of a given string is called a suffix. “ATTG”, “TTG”, “TG”, and “G” are suffixes of the character string “ATTG”.
When the suffixes of the given character string s are arranged in dictionary order, an integer array that records the order in which the suffixes start in the character string s is defined as a suffix array (Suffix Array: SA). (See "Gusfield, Algorithms on Strings, Trees and Sequences, Cambridge University Press, 1997"). In the case of s = ATTG, if the suffixes are arranged in dictionary order, they are "ATTG", "G", "TG", and "TTG". Therefore, the suffix array is SA = [1, 4, 3, 2]. . The i-th element of the suffix array SA is denoted as SA [i].
"LCP (Longest common prefix)" is the longest common prefix when two character strings are compared.
The “LCP array” is an integer array obtained by calculating the length of the LCP for all suffixes adjacent in the suffix array. In the above example of the suffix array, the LCP array is [0, 0, 1]. For the calculation of the LCP sequence, a known method (see “Kasai et al., Proc. CPM pp. 181-192, 2001”) can be used.
The "Burrows-Wheeler transform (BWT)" of the given character string s is a character string having the same length as | s |, and the ith character is s [SA [i] -1]. It is a character string (see “Navarro & Makkinen, ACM Computing Surveys 39 (1): Article No. 2, 2007”). However, when SA [i] = 1, the i-th character of the BWT is “＄”. "@" Is a new character that is not in the alphabet, called terminal symbol.
"FM-index" is a data structure obtained by extending BWT so that it can be used as a search index for searching for an arbitrary character string included in the original character string ("Ferragina & Manzini Proc.FOCS pp.390- 398,2000 ").
"Bit vector" is a sequence obtained by arranging 0s and 1s. For example, “0100110” is a bit vector.
"Function rank (S, a, i)" returns the number of times a appears up to the i-th character string or bit vector S. For example, rank (ATTG, T, 3) = 2.
"Function select (S, a, i)" returns the position where the i-th "a" appears from the beginning in the character string or bit vector S. For example, select (ATTG, T, 2) = 3.
The "dictionary" of a bit vector is a data structure for high-speed calculation of a rank function and a select function for a bit vector ("Navarro & Makkinen, ACM Computing Surveys 39 (1): Article No. 2 , 2007 ”).
"Mutation" is a difference in DNA sequence derived from individual differences or the like.
-"Allyls" are different DNA sequences that are generated by mutation. For example, at the position of a certain mutation in the genomic sequence, the allele of the old genomic sequence may be A and the allele of the read sequence may be T.

(3) Example 1
Next, a mechanism for efficiently mapping the read sequence 107 to the new genome sequence 114 using the mapping result for the old genome sequence 108 will be described.

(3-1) Prerequisites of Embodiment For the sake of simplicity, the following assumption is adopted in Embodiment 1. First, in Example 1, it is assumed that the new genome sequence 114 is one character string. When dealing with multiple chromosomes, this assumption does not hold, but for example, adding a new character "#" to the alphabet and concatenating multiple sequences of the new genome sequence across "#" into one character string In a way, this constraint can be avoided. This method is a general method when handling a plurality of character strings in a suffix array.

  It is also assumed that the mutation detected during processing is a single base mutation. In the case of treating a mutation of a plurality of bases, the treatment may be performed assuming that a plurality of single base mutations are present continuously.

  Generally, the process of aligning the read sequence 107 with the reference genomic sequence includes a process of calculating a position on the reference genomic sequence to which the read sequence 107 can be mapped (hereinafter, referred to as a “seed”), and thereafter, a process of calculating the seed. (Seed) A comparison is made between whether the surrounding character string matches the read sequence 107 while allowing a slight sequence difference, and it is determined whether mapping is possible (hereinafter referred to as "extension"). It consists of.

  In the following embodiment, a procedure for calculating a seed in mapping is provided. A known method such as “Smith-Waterman method (Smith & Waterman, journal of Molecular Biology 147: 195-197, 1981)” can be used for the extension process.

  Prior to the start of the analysis process, two parameters are given. The first parameter L is the length of the read array 107. The length of the read sequence 107 was about 30 bases in the early next-generation sequencer, but in recent years may exceed 200 bases. L is set so as to match the read array length of the data to be processed. The second parameter is the length K (a natural number of 1 or more) of a perfect match required when searching for a mapping destination candidate. For example, in the case of L = 100 bases, if an error of up to 5% is allowed, at least one place has a perfect match of 15 bases, so that K = 15 and all perfect matches having a length of K bases or more are found. If so, it is possible to exhaustively enumerate the locations that can be mapped with an error of up to 5%.

  In the following description, it is assumed that the read sequence 107 includes the result of mapping the read sequence 107 to the old genome sequence 108. Furthermore, in this specification, the two terms array and character string are used without distinction. Similarly, the two words base and letter are used interchangeably.

(3-2) System Configuration FIG. 4 shows a configuration example of a DNA sequence analysis apparatus 100 according to the present embodiment and a DNA sequence analysis system 115 configured using the apparatus. The DNA sequence analyzer 100 includes a CPU (Central Processing Unit) 101, a main storage device (memory) 102, and an auxiliary storage device 103, and uses a removable medium 104 as necessary.

  The main storage device 102 is a memory such as a RAM (Random Access Memory) in which various programs executed by the CPU 101 and various data necessary for executing the programs by the CPU 101 are stored. The main storage device 102 stores a read sequence (1) 107-1, an old genome sequence (1) 108-1, a new genome index (1) 109-1 constructed based on a new genome, and a query table 110. The main storage device 102 also stores programs to be executed by the CPU 101. Through the execution of the program, the functions of a new genome index constructing unit 111, a query table updating unit 112, and a mapping destination searching unit 113, which will be described later, are provided.

  The auxiliary storage device 103 is a storage device such as an HDD (Hard Disk Drive), and records a read sequence (2) 107-2, an old genome sequence (2) 108-2, a new genome index (2) 109-2, and the like. Is done. The auxiliary storage device 103 may further record the new genome sequence (1) 114-1. The removable medium 104 is a storage medium such as a CD and a DVD that can be attached to and detached from the DNA sequence analyzer 100, and includes a read sequence (3) 107-3, an old genome sequence (3) 108-3, and a new genome index (3 ) 109-3, new genome sequence (2) 114-2, etc. are recorded.

  The DNA sequence analyzer 100 may be connected to a user interface unit 106, or may be connected to a storage device via a network 105 such as a LAN (Local Area Network). Here, the network 105 is incompatible with a LAN (Local Area Network) or the Internet. The network 105 does not need to be limited to a wired connection but may be a wireless connection. The user interface unit 106 is an input / output device that provides a user interface, and includes, for example, a keyboard, a mouse, and a display. The user interface unit 106 may be integrated with the DNA sequence analyzer 100.

  The storage devices connected via the network 105 include a read sequence (4) 107-4, an old genome sequence (4) 108-4, a new genome index (4) 109-4, and a new genome sequence (3). 114-3 can be recorded. In this embodiment, a configuration in which the user interface unit 106 and these storage devices are connected to the DNA sequence analysis device 100 is referred to as a DNA sequence analysis system 115. However, this division is not fixed, and for example, a system that executes a part of a program in cooperation with another computer connected via the network 105 may be referred to as a DNA sequence analysis system. Further, the user may have only the user interface unit 106 and employ a mechanism for confirming the result of processing executed on the server connected via the network. The distinction between the DNA analysis device 100 and the DNA analysis system 115 is for convenience, and the device configuration may be the same.

  In FIG. 4, symbols (1) to (4) representing respective data are given according to the storage destination of the data even for the same type of data. 107, an old genome sequence 108, a new genome index 109, and a new genome sequence 114.

  These data may be stored in the main storage device 102 as necessary, for example, when reading and writing by the CPU 101, or when the power of the DNA sequence analysis device 100 is turned off or the free space of the main storage device 102 is exhausted. In this case, the data may be copied from the main storage device 102 to another storage medium or storage device (including a storage device). Before the analysis is started, data is stored in a removable medium 104 or an external storage device from another device (not shown), and the stored data is stored in the DNA sequence analysis device 100 when the DNA sequence analysis device 100 is started or analysis is started. The data may be copied to the main storage device 102 or the auxiliary storage device 103 in the analyzer 100.

(3-3) Processing Operation Hereinafter, a processing operation performed by the DNA sequence analyzer 100 will be described with reference to FIGS. 5 and 6. FIG. 5 shows an outline of processing until a candidate for mapping the read sequence 107 to the new genome sequence 114 is output. FIG. 6 is a processing sequence for explaining the cooperative operation of each processing unit. The individual processes will be described with reference to these drawings.

(3-3-1) Processing of the New Genome Index Construction Unit 111 The new genome index construction unit 111 searches for a search index (ie, a new genome index) for facilitating the search for an arbitrary partial sequence included in the new genome sequence 114. 109). The new genome index 109 is composed of two data. One is FM-index, which is a general character string search index. The remaining one is a bit vector based on the LCP sequence, and is data for completely searching for a match between the old genome index and the new genome index 109.

FIG. 7 shows a processing operation executed by the new genome index construction unit 111. First, the new genome index construction unit 111 constructs an FM-index for the new genome sequence 114 (S101). Next, the new genome index construction unit 111 calculates an LCP sequence for the new genome sequence 114 (S102). Thereafter, the new genome index construction unit 111 constructs a bit vector _BL having the same length as the LCP array calculated in S102 (S103). Here, the i-th numerical value of B _L is “1” when a numerical value equal to or more than K is written at the i-th value of the LCP array, and is “0” otherwise. Subsequently, the new genome index constructing unit 111 constructs a dictionary for calculating rank and select at high speed with respect to the bit vector _BL (S104).

(3-3-2) Process of Query Table Update Unit 112 When the new genome index constructing unit 111 constructs a dictionary, the process of the query table update unit 112 starts. The query table updating unit 112 executes a process of reading the old genome sequence 108 one base at a time and recording a mutation position in the query table 110. By updating (preparing) the query table 110 here, a plurality of read sequences 107 mapped to the same location of the old genome sequence 108 can be collectively mapped to the new genome sequence 114, and the mapping efficiency is improved. It is greatly improved.

  However, the read sequence 107 does not always completely match the mapped old genome sequence 108, and a difference may occur due to mutation. Therefore, in the present embodiment, a sequence constructed by reflecting the allele of the mutation contained in the read sequence 107 in the old genome sequence 108, instead of the old genome sequence 108 itself where the read sequence 107 is mapped, is used for the search. Adopt the method. In order to realize this method, the query table updating unit 112 records the mutation found in the read sequence 107 in the query table 110, and manages information on the mutation to be reflected in the old genome sequence 108.

  The query table updating unit 112 is used to update the query table 110 when the query table 110 is constructed (initial state) and when searching. FIG. 8 shows an example of the query table 110. FIG. 8 illustrates not the data of the query table 110 itself but a table format for explanation. The query table stores combinations of alleles of mutations existing within L bases from the processing position on the old genome sequence 108. This combination is hereinafter referred to as query allele 1101. A sequence obtained by modifying the old genome sequence 108 by reflecting the query allele 1101 is called a query sequence. For each query allele 1101, the length of the query sequence corresponding to each query allele 1101 (column 1102), the range in which the query sequence exists in the new genome index (column 1103), and the allele at each mutation position are recorded in column 1104. I do. A list 1105 of alleles appearing in the read sequence 107 may be recorded for each mutation position.

  FIG. 9 shows a processing operation executed by the query table update unit 112. Here, the query table updating unit 112 updates the query table 110 so that the mutation of the read sequence 107 mapped near the processing position (on the old genome sequence 108) of the mapping destination searching unit 113 is recorded in the query table 110. Is updated as appropriate.

  First, the query table updating unit 112 removes a mutation that is far from the vicinity (S201 to S203). Specifically, in S201, the query table updating unit 112 determines whether or not the query table 110 has a mutation at a position that is at least L (a natural number equal to or greater than 1) bases from the position being processed. If there is such a mutation, the query table update unit 112 proceeds to S202, and deletes the column in which the mutation is recorded from the query table 110. Next, the query table updating unit 112 proceeds to S203, and integrates the same query allele 1101. At this time, the query table updating unit 112 compares the individual query alleles 1101 recorded in the query table 110, and when there is a query allele 1101 composed of the same combination of alleles having the same length value, Only one query allele 1101 is left, and the other query allele 1101 is deleted.

  Subsequently, the query table updating unit 112 adds the newly discovered mutation to the query table 110 (S204 to S205). Specifically, in S204, the query table update unit 112 adds a new column of the query table 110. In S205, the query table updating unit 112 records each allele of the old genome sequence 108 and the read sequence 107 so as to be associated with a number indicating the position on the old genome sequence 108 determined to have a mutation.

  When these processes are completed, the query table updating unit 112 executes a process (S208 to S209) for adding a query allele when the query table 110 before the process is empty, and when the query table 110 is not empty, The processing of updating the query allele (S208, S210) is executed. Specifically, in S208, the query table updating unit 112 determines whether one or more query alleles 1101 already exist in the query table 110.

  If not present (in S209), the query table updating unit 112 adds a new query allele 1101 having the allele of the read sequence 107 due to the found mutation to the query table 110. However, when there is no mutation at the position being processed, as shown in FIG. 13, a row having no corresponding mutation is added instead of the query allele. Here, the “range of the new genome index” in the query table is the entire range ([1, new genome sequence length]), and the “length” is 1. On the other hand, when it is determined in S208 that the existing query allele 1101 exists (in S210), the query table update unit 112 executes an update process (FIG. 10) of each existing query allele 1101.

  FIG. 10 shows a detailed operation of the update process executed in S210. First, the query table updating unit 112 sequentially processes all of the existing query alleles 1101 (S2101, S2109). Specifically, in S2101, the query table update unit 112 selects one unprocessed query allele 1101 from the query table 110. In step S2109, the query table update unit 112 determines whether an unprocessed query allele 1101 remains, and returns to step S2101 if it remains. By this loop processing, the processing of S2102 to S2108 described later is repeated for each query allele 1101 one by one.

  Next, the query table updating unit 112 sequentially processes all the alleles of the newly added mutation (S2102, S2105). Specifically, in S2102, the query table update unit 112 selects one unprocessed allele from the allele of the new mutation. In S2105, the query table update unit 112 determines whether or not unprocessed alleles of the new mutation remain, and returns to S2102 if alleles remain. By this loop processing, the processing of S2103 to S2104 described later is repeated for all new mutation alleles one by one.

  Subsequently, the query table updating unit 112 determines whether there is a read sequence 107 including both the new mutation allele and the existing query allele 1101 (S2103). When there is such a read sequence, the query table updating unit 112 adds a new query allele 1101 obtained by combining the new mutation allele and the existing query allele 1101 to the query table 110 (S2104). On the other hand, when a negative result is obtained in S2103 (when there is no read sequence 107 including any of the newly added alleles of the mutation), the query table update unit 112 skips S2104 and proceeds to S2105.

  When a positive result is obtained in S2105, the query table updating unit 112 determines whether the additional process of S2104 has not been performed on any of the new alleles of the existing query allele 1101 being processed (S2106). . Here, when a positive result is obtained, the query table updating unit 112 executes a process at the time of loss of the matching area to output the matching area found so far (S2107). On the other hand, when a negative result is obtained in S2106, or when the process of S2107 is completed, the query table updating unit 112 deletes the existing query allele 1101 that becomes unnecessary by performing the above-described process (S2108).

  With reference to FIG. 11, how the query table 110 is changed by the update processing by the query table update unit 112 will be described. FIG. 11 shows the query table 110 after updating the query table 110 shown in FIG. FIG. 8 assumes the state after processing the “1104” base of the old genome sequence, and FIG. 11 assumes the state after processing the “1103” base. Also assume that L = 100 and K = 15.

  In the case of the query table 110 shown in FIG. 8, there is no base that is at least L = 100 bases away from the “1103” -th position, and thus there is no column deleted in S201 to S203. Next, the mutation found at the “1103” th position of the old genome sequence 108 is added to the query table 110 of FIG. 11 as a new column (S204 to S205). Further, the query alleles 1101 written in the query table 110 are updated one by one by the query table update processing (S208, S210) (S2101, S2109).

  First, consider the case where G and A described in the first row of the query allele 1101 in the query table 110 in FIG. 8 are processed. In this example, since the new mutation found at the “1103” th position is T or C, two of “T, G, A” and “C, G, A” are included in the query allele 1101 with these added at the beginning. is there. If these allele combinations do not exist in the read sequence 107, there is no need to record them in the query table 110 and process them. Therefore, the presence or absence of the read sequence 107 including these is determined (S2103), and is added to the query table 110 only when “exist” is present (S2104). In FIG. 11, it is assumed that both the read sequences 107 including “T, G, A” and “C, G, A” exist, and both are added to the query table 110. The query alleles 1101 of only G and A among the query alleles 1101 described in the query table 110 of FIG. 8 are unnecessary and are deleted (S2108).

  Next, consider the case where C and A described in the second row of the query allele 1101 in the query table 110 of FIG. 11 are processed. When combined with the new mutation found at the "1103" position, there are two possibilities: "T, C, A" and "C, C, A". However, it is determined in S2103 that none of these exist in the read sequence 107, and they are not left in FIG. Instead, in S2106, it is determined that there is no read sequence coexisting with the new mutation, and the process at the time of loss of the matching region (S2107) is executed.

  In the query table 1101 of FIG. 11, only “C, G, T” from G and T described in the third row of the query allele 1101 is described from C and T described in the fourth row. Since only “T, C, T” exists in the read sequence 107, it is added to the query table 110 as shown in FIG. A query allele 1101 of only C and a query allele 1101 of only T are added to the query table 110 of FIG. This is added because it was found in S3074 (FIG. 15) described later that the prefix of length K = 15 coincides with a new location on the new genome sequence.

(3-3-3) Processing of the Mapping Destination Searching Unit 113 The mapping destination searching unit 113 converts the character string obtained by reflecting the mutation described in the query table 110 into the old genome sequence 108 into the new genome sequence 114. Is searched using the new genome index 109. By the way, the mapping destination search unit 113 processes the bases one by one from the right end to the left end of the old genome sequence 108 and, while referring to the query table 110, as long as the same sequence as the query sequence exists in the new genome sequence 114, Extend the sequence one base at a time.

  FIG. 12 shows a processing operation executed by the mapping destination search unit 113. First, the mapping destination search unit 113 performs an initialization process (S301). Specifically, the variable i is initialized to an integer value equal to the length of the old genome sequence 108. Next, the mapping destination search unit 113 executes a process of updating the query table 110 (S302). Specifically, the processing of the query table updating unit 112 is executed for the i-th position of the old genome sequence 108.

  For each query allele 1101 of the query table 110, the mapping destination search unit 113 executes a process of extending the matching position in the new genome sequence 114 by one base (S303 to S308). Specifically, in S303, the mapping destination search unit 113 selects one unprocessed query allele 1101 from among the query alleles 1101 recorded in the query table 110, and uses the query allele 1101 to generate i of the old genome sequence 108. A query sequence is generated by modifying the length sequence recorded in the “length” column 1102 of the query table 110, starting from the base number.

  However, when there is no mutation near the position being processed, there is a case where there is no mutation as in the query table 110 shown in FIG. In this case, the “length” sequence recorded in column 1102 of the query table 110, starting from the i-th base of the old genome sequence 108, is used as it is as the query sequence. In the present specification, such a row indicating a query sequence that has no mutation and matches the old genome sequence 108 is also treated in the same manner as the query allele 1101.

Next, in S304, the mapping destination search unit 113 searches the first character of the query sequence including the query allele 1101 and having a length obtained by adding 1 to the value recorded in the “length” column 1102. , And the suffix starting from the second character is x, and the range of the new genome index 109 where the character string bx exists is calculated. This range can be calculated as follows using LF-mapping (Navarro & Makkinen, ACM Computing Surveys 39 (1): Article No. 2, 2007).
beg (bx) = C [b] + rank (BWT, b, beg (x) -1) +1
end (bx) = C [b] + rank (BWT, b, end (x))

  Here, C [b] is the sum of the total number of characters smaller than b in the dictionary order in the new genome sequence. Beg (x) and end (x) are the start position and end position described in the range (column 1103) of the new genome index 109, respectively.

  In the next step S305, the mapping destination search unit 113 determines whether or not the calculation result is beg (bx)> end (bx). If beg (bx)> end (bx), it is known from the nature of LF-mapping that the character string bx does not exist in the new genome sequence 114. In this case, the mapping destination search unit 113 sets the matching region The process at the time of loss is executed (S306).

  FIG. 14 shows the detailed operation of S306. In this process, the mapping destination search unit 113 executes a process of judging and outputting a necessary portion from the matching portion between the query sequence and the new genome sequence 114 which has been extended until immediately before the loss of the matching region. First, in S3061, the mapping destination searching unit 113 determines whether or not | x | ≧ K. If | x | <K, nothing is output. In S3062, the mapping destination search unit 113 determines whether or not the output range is already output when processing another query allele or the like. If the data has already been output, the mapping destination search unit 113 does not output the data again. In S3063, the mapping destination search unit 113 outputs a range that has not been output for the value of i + 1 among the “range of the new genome index” as a matching area for i + 1.

  Returning to the description of FIG. In step S305, if beg (bx) ≦ end (bx), the character string bx is found to exist in the new genome sequence from the property of LF-mapping. In this case, the mapping destination search unit 113 The process at the time of extension of the matching area is executed (S307).

  FIG. 15 shows the detailed operation of S307. In this process, the mapping destination search unit 113 executes approximately three processes. The first process is a process of updating the range in which the character string bx obtained by extending the query sequence x corresponding to each query allele recorded in the query table information by one character to the left appears in the new genome. (S3071). The second process is a process of outputting a matching area that has reached the read length L (S3072 to S3073). The third process is a process of increasing the number of matching points by shortening the query sequence bx (S3074 to S3076), and calculating a new position on the new genome sequence 114 where the prefix of the length K of bx matches. I do. Since a match with a length less than K is not an output target, no processing is required, and a match with a length of K + 1 or more can be processed by extending a match having a length of K or more already searched.

  Here, in S3071, the mapping destination search unit 113 replaces “the range of the new genome index” of the query allele 1101 to be processed with (beg (bx), end (bx)). Further, the mapping destination search unit 113 adds 1 to the “length” of the query allele 1101. In step S3072, the mapping destination searching unit 113 determines whether the “length” has become L or more.

If the “length” is equal to or greater than L, the mapping destination search unit 113 proceeds to S3073, and if there is a read sequence mapped at position i of the old genome sequence, i and “the range of the new genome index after replacement” , The range not yet output for the same value of i is output, and the "length" is returned to L-1. When the “length” is less than L or after S3073 has been executed, the mapping destination search unit 113 proceeds to S3074, where the prefix y of the length K of bx is other than the prefix y of bx in the new genome sequence. Check if it exists at the location. For this purpose, the mapping destination search unit 113 calculates the following two values beg (y) and end (y) using the bit vector _BL .
beg (y) = select (B _L , 0, rank (B _L , 0, beg (bx) -1)) + 1
end (y) = select (B _L , 0, rank (B _L , 0, end (bx) -1) +1)

Then, it is determined whether or not beg (y) ≠ beg (bx) or end (y) ≠ end (bx) holds. If this condition is satisfied, it is determined that y exists at a position other than the bx prefix in the new genome sequence. In the above equation, subtracting 1 from beg (x) and end (bx) means that the bet (x) -th and end (x) -th elements of the BWT are the beg (x) -first of the _BL , end (x) -1 to correspond to the first element.

  In step S3075, the mapping destination search unit 113 determines whether the condition in step S3074 is satisfied. If the condition is satisfied, the process advances to step S3076 to add a new query allele combination to the query table 110. In this query allele, all the alleles within K bases are copied in the above query allele, the “range on the new genome index” is [beg (y), end (y)], and the “length” is K. If there is no mutation within K bases, a row having no mutation is added to the query table as in the example of FIG. If the condition of S3074 is not satisfied, the mapping destination search unit 113 ends the processing of S307.

  Returning to the description of FIG. When the above process ends, S308 is executed. In S308, the mapping destination search unit 113 determines whether or not the unprocessed query allele 1101 remains in the query table 110, and if it remains, repeats the processes of S303 to S308 described above. If no unprocessed query allele 1101 remains, the mapping destination detection unit 113 executes a process of searching for a new matching area, which is not an extension of the existing matching area, in S309.

  FIG. 16 shows the detailed operation of S309. In this processing, the mapping destination search unit 113 examines combinations of alleles that are not present in the existing query allele 1101 (S3091 to S3093, S3098), reflects the combination in the K base region of the old genome sequence (S3094), and Is searched for K bases or more (S3095). If there is a matching part, the mapping destination search unit 113 adds it to the query table 110 (S3096 to S3097).

  Here, in S3091, the mapping destination searching unit 113 lists all combinations of all mutations within K bases from the position on the old genome sequence 108 indicated by the variable i. When there is no mutation within K bases, the same processing as when there is an allele completely matching the old genome sequence is performed. In the next step S3092, the mapping destination search unit 113 selects one of the combinations of alleles listed in step S3091. In the next step S3093, the mapping destination searching unit 113 determines whether or not there is a read sequence containing the combination of alleles selected in S3092. If not, the mapping destination search unit 113 moves to S3098 and proceeds to the next combination of alleles.

  On the other hand, if a positive result is obtained in S3093, the mapping destination search unit 113 proceeds to S3094 and adds the selected allele combination to the K base old genome sequence 108 from the position on the old genome sequence 108 indicated by the variable i. Generate the applied array. In the next step S3095, the mapping destination search unit 113 determines whether or not the sequence in S3094 exists in the new genome sequence 114 using the new genome index 109.

  If it is determined that the combination does not exist (No in S3096), the mapping destination search unit 113 moves to S3098 and proceeds to the next combination of alleles. On the other hand, if it is determined that the combination exists (Yes in S3096), the mapping destination search unit 113 moves to S3097, and if the combination of alleles processed after S3092 does not yet exist in the query table, queries as a new query allele 1101. Add to Table 110. If it is determined in S3091 that there is no mutation within K bases and that there is a completely identical allele, a row having no mutation as illustrated in FIG. 13 is added to the query table 110. Then, in S3098, if there is any unprocessed combination of alleles, the mapping destination search unit 113 returns to S3092 to process another combination of alleles.

  Returning to the description of FIG. When the above-described processing is completed, the mapping destination search unit 113 moves one base at a time from the right end to the left end of the old genome sequence 108 (S310 to S311). Specifically, in S310, the mapping destination search unit 113 subtracts 1 from the variable i. Further, in S311, if the value of the variable i is 1 or more, the process returns to S302, and the above-described processing is repeated.

(3-4) Effect By executing the above-described processing, in S3073, for any sequence within the length L obtained by applying the mutation found in the read sequence 107 to the old genome sequence 108, the position that appears in the new genome sequence 114 is It can output comprehensively. Therefore, according to the present embodiment, candidates for the position in the new genome sequence 114 to which the read sequence 107 mapped to the same position in the old genome sequence 108 is to be mapped can be efficiently calculated in a lump. .

(4) Example 2
In order to evaluate the validity of the mapping result, it is desirable that the user can inspect the contents of the query table 110 occupying an important part in the DNA analyzer 100. Therefore, in the present embodiment, a mechanism for providing a graphical interface for the user to check the contents of the query table 110 will be described.

  FIG. 17 shows an example of the interface screen 1700. The interface screen 1700 is displayed on a display screen constituting the user interface unit 106. The interface screen 1700 has a display area 1701 for displaying the contents of the query table 110. In the lower part of the display area 1701, a display area 1710 of a query sequence corresponding to the query table 110 is displayed.

  In the display area 1710, a straight line 1703 representing a query sequence is displayed in parallel with a straight line (horizontal axis) 1702 representing the old genome sequence 108. On the straight line 1703, for each position 1704 of each mutation in the query table 110, the allele 1705 of the query sequence at that position is displayed, and the number 1706 of the read sequences 107 that match each query sequence (for example, “3 ×” is 3) This indicates a book.) Is displayed. The range on the old genome sequence 108 to be displayed on the interface screen 1700 can be arbitrarily selected by the user using an interface such as the direction button 1707, for example. The query table 110 records alleles of mutations occurring in the neighboring (less than L) region. That is, since the query table 110 is irrelevant to an area separated by L or more, the query table of the area that the user wants to observe can be efficiently reconstructed.

(5) Other Embodiments The above embodiments show preferred examples, and are not intended to limit the present invention to a specific configuration. Various changes can be made without departing from the gist of the above description. For example, it is not always necessary to provide all the components of the above-described embodiment. Further, a part of one embodiment can be replaced with the configuration of another embodiment. Further, the configuration of one embodiment can be added to the configuration of another embodiment. Further, for a part of the configuration of each embodiment, a part of the configuration of another embodiment can be added, deleted, or replaced.

  In addition, each of the above-described configurations, functions, processing units, processing means, and the like may be partially or entirely realized by hardware, for example, by designing an integrated circuit. Further, each of the above-described configurations, functions, and the like may be realized by a processor interpreting and executing a program that realizes each function (that is, as software). Information such as a program, a table, and a file for realizing each function can be stored in a storage device such as a memory, a hard disk, an SSD (Solid State Drive), or a storage medium such as an IC card, an SD card, and a DVD. Further, the control lines and the information lines indicate those considered to be necessary for the description, and do not represent all the control lines and the information lines necessary for the product. In fact, it can be considered that almost all components are interconnected.

100: DNA sequence analyzer 101: CPU
102 main storage device 103 auxiliary storage device 104 removable medium 105 network 106 user interface unit 107 mapping result of read sequence and its old genome sequence 108 old genome sequence 109 new genome index 110 query table 111 ... New genome index constructing unit 112 Query table updating unit 113 Mapping destination searching unit 114 New genome sequence 115 DNA sequence analysis system 401 Lead sequence mutation 1101 Query alleles 1102 stored in the query table 1102 Length 1103 of the corresponding query sequence ... Range in the new genome index corresponding to each query allele 1104 ... Allele at each mutation position corresponding to each query allele 1105 ... List of alleles found in the read sequence at each mutation position 17 Reference numeral 1 denotes a display area for displaying the contents of the query table 1702 ... a straight line 1703 representing the old genome sequence 1703 ... a straight line 1704 representing the query sequence 1705 of the position of each mutation in the query table 1705 ... alleles of the query sequence at the position 1706 ... Number of matching lead sequences 1707: Direction buttons for changing the range to be displayed

Claims (9)

A storage for storing a first genomic sequence called an old genomic sequence, a second genomic sequence called a new genomic sequence, a large number of DNA sequences called a read sequence, and a result of mapping the read sequence to the old genomic sequence. Equipment and
A new genome index constructing unit that constructs an index for searching for any character string present in the new genome sequence,
A data structure called a query table that records a combination of mutations in the read sequence that is a mutation with respect to the old genome sequence and is located near the processing position on the old genome sequence, based on the result of the mapping. A query table update unit to be constructed,
A sequence constructed by applying the combination of mutations stored in the query table to the old genomic sequence is compared with the new genomic sequence, and K (1 or more natural number) bases or more in the constructed sequence are compared with the new genomic sequence. A DNA sequence analysis device comprising: a mapping destination search unit that comprehensively outputs a part that completely matches a genome sequence. The DNA sequence analyzer according to claim 1,
The new genome index construction unit,
As the index,
FM-index,
When comparing the suffixes corresponding to the adjacent numbers in the suffix array, construct a bit vector in which the length of the longest prefix that matches is 1 at a position exceeding the K and 0 at other positions. A DNA sequence analyzer characterized by the above-mentioned. The DNA sequence analyzer according to claim 1,
The query table update unit,
A combination of allele combinations of mutations within L (one or more natural numbers) bases, wherein allele combinations present in both the read sequence and the new genome sequence are stored in the query table. Sequence analyzer. The DNA sequence analyzer according to claim 2,
The mapping destination search unit,
Using the bit vector, searching for a position on the new genomic sequence where the sequence reflecting the allele of the read sequence in the old genomic sequence and the new genomic sequence completely match at least K bases. DNA sequence analyzer. A storage for storing a first genomic sequence called an old genomic sequence, a second genomic sequence called a new genomic sequence, a large number of DNA sequences called a read sequence, and a result of mapping the read sequence to the old genomic sequence. Equipment and
A new genome index constructing unit that constructs an index for searching for any character string present in the new genome sequence,
A data structure called a query table that records a combination of mutations in the read sequence that is a mutation with respect to the old genome sequence and is located near the processing position on the old genome sequence, based on the result of the mapping. A query table update unit to be constructed,
A sequence constructed by applying the combination of mutations stored in the query table to the old genomic sequence is compared with the new genomic sequence, and K (1 or more natural number) bases or more in the constructed sequence are compared with the new genomic sequence. A mapping destination search unit that comprehensively outputs a part that completely matches the genome sequence,
A first region for displaying the content and position of the mutation appearing in the read sequence in association with each other, the combination of the mutation, the length of the corresponding query sequence on the new genome sequence, and the corresponding query sequence And a user interface unit for displaying a confirmation screen including a second area for displaying a range present in the index in association with the second area. The DNA sequence analysis system according to claim 5,
The user interface unit displays one or a plurality of straight lines corresponding to the query sequence in parallel with a first axis corresponding to the old genome sequence, and is displayed in the query table in the straight line. A DNA sequence analysis system, wherein the content of the mutation in the query sequence is displayed for each mutation position. DNA sequence analyzer,
According to the processing contents, a first genomic sequence called an old genomic sequence, a second genomic sequence called a new genomic sequence, a large number of DNA sequences called a read sequence, and the read sequence were mapped to the old genomic sequence. Reading a part or all of the result from the storage device;
A process of constructing an index for searching for any character string present in the new genome sequence,
A data structure called a query table that records a combination of mutations in the read sequence that is a mutation with respect to the old genome sequence and is located near the processing position on the old genome sequence, based on the result of the mapping. The process to build,
A sequence constructed by applying the combination of mutations stored in the query table to the old genomic sequence is compared with the new genomic sequence, and K (1 or more natural number) bases or more in the constructed sequence are compared with the new genomic sequence. And a process of comprehensively outputting a part that completely matches the genome sequence. The DNA sequence analysis method according to claim 7,
The DNA sequence analyzer,
A first region for displaying the content and position of the mutation appearing in the read sequence in association with each other, the combination of the mutation, the length of the corresponding query sequence on the new genome sequence, and the corresponding query sequence Further comprising the step of: displaying a confirmation screen including a second area for displaying a range associated with the index in the second area. The DNA sequence analysis method according to claim 8,
The DNA sequence analyzer,
One or more straight lines corresponding to the query sequence are displayed in parallel to a first axis corresponding to the old genome sequence, and for each mutation position displayed in the query table on the straight line A DNA sequence analysis method, further comprising the step of displaying the contents of the mutation in the query sequence.

Images