Sequence data analyzer, DNA analysis system and sequence data analysis method

Abstract

Sequence data analyzer (1) comprising: a read directory preparation unit (21) configured to generate a concatenation series for each of the sample DNA fragments, obtained by joining a pair of a left sequence (11a) and a right sequence (11b) by a concatenation character, both obtained by sequencing each of the sample DNA fragments from the two ends of each of the sample DNA fragments, and to prepare a read sequence directory (14) based on a series obtained by joining the concatenation series of the sample DNA fragments with end characters; a query recovery unit (23) configured to start searching the read sequence directory (14) for a hit position (17a), which is a base character coordinate, at which a query sequence generated from a base series of a genome DNA appears;a sample reconstruction unit (25) configured to start searching the read sequence directory (14) for the sequence from the hit position (17a) to the end character located on both sides of the hit position (17a) in order to recover the sequence as a sample sequence; to start searching the sample sequence for the concatenation character located on both sides of the hit position (17a) from the hit position (17a) and to extract as a pairing sequence from the retrieved concatenation character to the end character either on the left sequence (11a) or the right sequence (11b) where the hit position (17a) does not exist; and a mapping unit (26) configured to search the genome DNA for a base character coordinate of the base sequence in which the pairing sequence appears.

Description

CROSS-REFERENCE TO RELATED REGISTRATIONS

[0001] The present application claims priority pursuant to 35 USC §119 of Japanese patent application No. 2014-77278, filed on April 3, 2014, entitled “SEQUENCE DATA ANALYZER, DNA ANALYSIS SYSTEM AND SEQUENCE DATA ANALYSIS METHOD”. The contents of that application are incorporated herein in full by reference. BACKGROUND OF THE INVENTION Area of ​​the invention

[0002] The present invention relates to a sequence data analyzer, a DNA analysis system and a sequence data analysis method. Discussion of the background

[0003] The base sequence of a genome DNA fragment has already been sequenced in its entirety, and the base sequences are primarily exposed on servers connected to the internet. Researchers perform a comparison (genome mapping) of base positions between sample DNA fragments from a subject, sequenced by an automated sequencing system, and the genome DNA as a model (reference data) to detect mutations in the sample DNA. A mutation is, for example, a different part, such as a single nucleotide polymorphism (SNP) or a structural mutation (SV), between a base sequence of a genome DNA fragment and a base sequence of a sample DNA fragment. A sample DNA fragment is obtained by performing fragmentation processing on a piece of sample DNA in an automated sequencing system.

[0004] Since the sequencing machine can only sequence a limited length of a DNA sequence, a paired-ends method is known for handling two pieces of (paired) read sequences. This method involves sequencing a limited length (only a portion of a sample DNA fragment) from each of the two ends of sample DNA fragments that are essentially aligned at the same length. In other words, a paired-ends method leaves a section in the middle of a sample DNA fragment that does not belong to any of the paired read sequences and is not to be sequenced.In genome mapping of a paired read sequence, which is sequenced by a method of paired ends at a base position in a genome DNA, it is therefore possible to map a piece of a sample DNA fragment in a genome DNA with high accuracy by position mapping together with the other read sequence instead of by position mapping only one read sequence of a pair.

[0005] Since the amount of data from base sequences of a genome DNA and sample DNA fragments is enormous when two pieces of (paired) read sequences in a genome DNA are mapped using a paired-end method, as described above, processing efficiency and data compression are required. Consequently, various technologies for improving genome mapping efficiency have been proposed, as described below.

[0006] The diTag of paired ends, in which the sequence direction from the 5' end to the 3' end is maintained on a full-length cDNA, is sequenced in a genome-wide transcriptome profiling and a method for efficiently detecting diTag sequences adjacent to a defined spacer sequence is known (Patent Literature 1).

[0007] In order to quickly and accurately align the enormous number of short readings generated by a new DNA sequencing machine to a reference genome sequence, the technology of performing a Burrows-Wheeler transformation on a genome sequence and rapidly identifying the position in a reference genome where the sequence corresponding to the leading several dozen germ sequences of a short reading appears is widely and extensively used (for example, Non-Patent Document 1).

[0008] To reduce ambiguity arising from the existence of multiple possible positions that should be aligned in this case, a paired short reading map was also extensively implemented (e.g., Non-Patent Literature 2). Normally, the same name (e.g., QNAME in Non-Patent Literature 1) is given to two short readings to indicate that they are paired.

[0009] A method for using the large amount of data obtained by mapping a large number of short reads onto a reference genome is generally and extensively used to perform mutation analysis (e.g., non-patent literature 3). To perform mutation analysis without the computationally expensive mapping calculation, a further possible method exists for using BW transformations on a reference genome sequence and short read sequences (patent literature 2, non-patent literature 4). Non-patent literature Patent literature 1: Japanese unexamined patent disclosure no. 2008-547080 Patent literature 2: Japanese unexamined patent disclosure no. 2013-038919 Non-patent literature 1: Li H., Handsaker B., Wysoker A., ​​Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup: The Sequence alignment / map (SAM) format and SAMtools. Bioinformatics, 25, 2078–9 (2009). Non-patent literature 2: Li H. and Durbin R. Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 25:1754–60 (2009). Non-patent literature 3: MA DePristo et al., A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43(5): 491–498 (2011). Non-patent literature 4: K. Kimura and A. Koike: A new approach to DNA sequence variation analysis using Burrows-Wheeler transform of massive short-read data, in Proceedings of Advances in Genome Biology & Technology Conference (AGBT) 2013, Marco Island, p. 202. SUMMARY OF THE INVENTION Problems to be solved by the invention

[0010] To handle the paired information in a computer for the large amount of read sequence data sequenced using a paired-ends method, identifiers (pointer information) are required to link one read sequence of a pair to another read sequence of the same pair. These identifiers require multiple bits to uniquely identify a pairing partner. This increases the size of the identifiers, leading to a decrease in computational efficiency. For example, if there are 8 billion short reads of 100 bases each, 4 billion identifiers are needed. Since each identifier requires at least 4 bytes, the total identifier size for all reads can reach up to 4 GB (gigabytes).

[0011] The present invention is intended to efficiently locate one read sequence from the other read sequence, both forming a pair, which is sequenced from both ends of a sample DNA fragment. Solution to the problem

[0012] To solve the preceding problem, a sequence data analyzer of the present invention comprises a read directory preparation unit configured to generate a concatenation sequence for each of the sample DNA fragments, obtained by joining a pair of left-hand sequences and right-hand sequences with a concatenation character, both obtained by sequencing each of the sample DNA fragments from the two ends of each of the sample DNA fragments, and to prepare a read sequence directory based on a sequence obtained by joining the concatenation sequences of the sample DNA fragments with end characters; a query recovery unit configured to start searching the read sequence directory for a hit position that is a base character coordinate at which a query sequence generated from a base row of a genome DNA appears; a sample reconstruction unit configured to start searching the read sequence directory for the sequence from the hit position to the end character located on both sides of the hit position in order to recover the sequence as a pattern sequence, to start searching the sample sequence from the hit position for the concatenation character located on both sides of the hit position, and to extract as a pairing sequence from the recovered concatenation character to the end character either on the left or the right sequence where the hit position does not exist; and an imaging unit configured to search the genome DNA for a base sign coordinate of the base series in which the mating sequence appears.

[0013] The other measures are described below. Advantageous effects of the invention

[0014] According to the present invention, one read sequence can be efficiently located from the other read sequence, both of which form a pair to be sequenced from both ends of a sample DNA fragment. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] The present invention is described in detail below with reference to the accompanying drawings, in which:

[0016] Fig. 1 is a configuration showing a DNA analysis system according to an embodiment of the present invention.

[0017] Fig. 2 is a configuration showing a sequence data analyzer according to an embodiment of the present invention.

[0018] Fig. 3 is a flowchart showing a process for generating each directory according to an embodiment of the present one.

[0019] Fig. Four explanatory drawings are shown, illustrating a processing operation for generating a BW series during the processing operation for generating a directory according to an embodiment of the present invention. Fig. 4A represents a process for generating concatenation sequences from read sequences. Fig. 4B represents a process for generating BW series from concatenation series.

[0020] Fig. The 5 explanatory drawings describe a conversion process for converting into a wavelet tree form during processing for generating directories according to an embodiment of the present invention. Fig. 5A represents the binary tree that is referenced for conversion to wavelet tree form. Fig. 5B represents the conversion processing for converting a BW series into the wavelet tree form.

[0021] Fig. 6 is a flowchart showing an SNP analysis processing according to an embodiment of the present invention.

[0022] Fig. Seven explanatory drawings are shown, illustrating all data obtained during SNP analysis processing according to one embodiment of the present invention. Fig. 7A displays SNP information. Fig. Section 7B presents query sequences generated by the SNP information and a genome sequence directory. Fig. 7C represents the state in which the query sequence is recovered from the read sequence directory. Fig. 7D shows the analysis results of the SNP.

[0023] Fig. Figure 8 is a flowchart showing the analysis processing of the structural mutation according to an embodiment of the present invention.

[0024] Fig. The 9 explanatory drawings show data in the analysis processing of a structural mutation according to an embodiment of the present invention. Fig. 9A shows structural mutation information. Fig. 9B presents the query sequences generated from the structural mutation information and the genome sequence directory. Fig. 9C is the state in which the query sequences are recovered from the read sequence directory. Fig. Figure 9D shows the analysis results of the structural mutations. DESCRIPTION OF THE EXECUTION FORMS

[0025] One embodiment of the present invention will now be described in detail with reference to the drawings.

[0026] Fig. 1 is a configuration that shows a DNA analysis system.

[0027] A sequence data analyzer 1 It is implemented using a computer, a server, or similar, with a configuration similar to a conventional computer.

[0028] The sequence data analyzer 1 includes a CPU (central processing unit) 201 , a storage 202 , which is a storage unit in which programs are mainly stored, a GUI (graphical user interface) for operation, a display unit 203 To display analysis results, an HDD (hard disk drive) is used. 204 , which acts as a storage unit, the sequence directories (a read sequence directory) 14 and a genome sequence directory 15 in Fig. 2) and the like, stores an input unit 205 a keyboard or similar device primarily for entering mutation information such as SNPs and parameters, and a network interface (NIF) 206mainly to connect to the internet, all via a bus 207 are connected.

[0029] The HDD 204 Stored sequence directories can be located either in a storage device located outside the sequence data analyzer. 1 is installed, or stored in a data center or the like via a network.

[0030] The various flowcharts described below are mainly implemented by executing programs on the CPU. 201 implemented.

[0031] The NIF 206 of the sequence data analyzer 1 is with a genome server 8 and a sequencing machine 9 connected by a network.

[0032] The sequencing machine 9 sequences a pair of sequences at both ends (read sequence at the 5' end and read sequence at the 3' end) of each sample DNA fragment to the sequence data analyzer 1to provide the result.

[0033] Furthermore, it is a common way of writing reading sequences (base sequences) to describe the base characters at the 5' end on the left and the base characters at the 3' end on the right; consequently, the 5' end is represented as 'left' below and the 3' end is represented as 'right'.

[0034] The sequencing machine 9 is configured as a superparallel type DNA sequencing machine (so-called next generation type) to sequence multiple (for example, 1 billion pieces of) sample DNA fragments in parallel.

[0035] The genome server also provides 8 the genome sequence obtained by sequencing a genome DNA with the sequence data analyzer 1 will be received.

[0036] Furthermore, it is possible to detect splice variants by using mRNA sequences for a reference sequence and analyzing the sequencing data of paired ends of cDNA samples, instead of using the genome sequence obtained by sequencing a genome DNA for a reference sequence and analyzing sample DNA fragments. This is because the loss of exons due to splice variants corresponds to the "deletion" of a structural mutation, and the integration of new exons corresponds to the "insertion" of a structural mutation.

[0037] Fig. 2 is a configuration that uses the sequence data analyzer 1 represents the outline of each component. Fig. 2 described and the details of each component are revealed below.

[0038] The sequence data analyzer 1 receives the inputs of a read sequence group (a left sequence) 11a and a right sequence 11b) from the sequencing machine 9 .

[0039] The reading sequence group 11 is a group of reading sequences (of the left sequence) 11a and the right sequence 11b ), which are obtained from each sample DNA fragment that is processed by the sequencing machine 9 is sequenced.

[0040] The left sequence 11a is a read sequence obtained by sequencing towards the 3' end from the endpoint at the 5' end as the base point of a sample DNA fragment.

[0041] The right sequence 11b is a read sequence obtained by sequencing towards the 5' end with the endpoint at the 3' end as the base point of a sample DNA fragment.

[0042] Here are the lengths of the left sequence 11a and the right sequence 11bFor example, approximately 100 bases, and assuming that the length of a sample DNA fragment is approximately 300 bases, there are approximately 100 bases in the middle that are not in the left sequence. 11a still in the right sequence 11b The fragment contains parts that are excluded from sequencing. Alternatively, a sample DNA fragment approximately 20,000 bases in length includes approximately 19,800 bases at positions excluded from sequencing.

[0043] A reading list preparation unit 21 prepares the reading sequence directory 14 from the reading sequence group 11 for preparation prior to DNA analysis processing. A genome inventory preparation unit. 22 prepares the genome sequence directory 15 from a reference genome sequence 12 before.

[0044] Furthermore, the reference genome sequence 12, which is determined in advance for each species to be analyzed, a group of sequences in a full length of each of the chromosomes.

[0045] The sequence analyzer 1 receives the inputs of SNP information 13a or structural mutation information 13b This information serves as analytical data for specifying the mutation to be analyzed this time. A SNP means that the base at a specific position in a sample DNA fragment differs from the base at the same position in a genome DNA sequence. A structural mutation involves the insertion and deletion of several continuous bases.

[0046] A query generation unit 23 generates a query sequence 16 with the analysis information 13 specified mutation with reference to the genome sequence directory.

[0047] The query recovery unit 24 wins a scoring position17a again, which is a base position coordinate that refers to the query sequence 16 (the query sequence 16 (appears) from the reading sequence directory 14 to be depicted.

[0048] A sample reconstruction unit 25 reconstructs reading sequences (sample sequences) 17 ) of a sample DNA fragment with the hit position 17a Here is the hit location 17a either in the left sequence 11a or in the right sequence 11b contained, consequently one of the reading sequences in which the hit position 17a not included, as a mating sequence 17b designated.

[0049] A mapping unit 26 refers to the genome sequence directory 15 , to specify a position in a genome DNA (genome mapping) from which the mating sequence 17b derived from the sample DNA fragment.

[0050] A sample determination unit 27 determines whether the sample DNA fragment matches the mating sequence 17a the analysis information 13 The specified mutation includes, according to success or failure in mapping the mating sequence. 17b Furthermore, the sample determination unit provides 27 the result of the analysis information 13 out of.

[0051] Fig. Figure 3 shows flowcharts that describe the creation process of each directory. First, the read directory creation steps (S101–S105) are shown, in which the read directory preparation unit 21 the reading sequence directory 15 from the reading sequence group 11 prepared, described.

[0052] In S101 the sequence data analyzer receives 1 the inputs of the reading sequence group 11 (the left sequence 11a and the right sequence 11b ) from the sequencing machine9 .

[0053] To facilitate understanding with Fig. To create 4A, two pieces of sample DNA fragments are needed. 301 and 305 depicted, and in the sample DNA fragment 301 “GA” will be added to the left sequence 11a counted as 2 bases and “T” at the right sequence 11b is counted as 1 base, while in the sample DNA fragment 305 “C” on the left sequence 11a is counted as 1 base and “TA” is at the right sequence 11b is counted as 2 bases.

[0054] The sequencing machine 9 notifies the sequence data analyzer 1 about the reading sequence group 11 , which are converted into characters by sequencing in FASTQ format 361 , that the left sequence 11a indicates, and signs 362 , that the right sequence 11b The reading sequences of two sample DNA fragments are indicated. 361 or signs362 listed, in which 4 lines of information about 1 reading sequence are described.

[0055] The first line “@seq1, @seq2” in the FASTQ format represents an identifier (Id) of a sample DNA fragment, and the second line “GA, T, C, TA” represents a read sequence. The correspondence between the first line “@seq1” of the character 361 and the first line “seq1” of the character 362 for example, it proves to be a pair that can be identified by sequencing the same DNA fragment 301 will be received.

[0056] In S102 in Fig. 3 links the reading list preparation unit 21 Using a concatenation symbol “&”, the base symbols of a pair of the left and right sequences are used to create a concatenation sequence that represents a segment of a sample DNA fragment. Additionally, the base symbols “A”, C, G, T” each represent four types of bases, and “N” represents an unknown base.

[0057] In Fig. 4A is a pair of the left sequence 11a “GA” and the right sequence 11b “T” of the same sample DNA fragment 301 with a concatenation sign “&” 302 concatenated and terminated with a "$" symbol 303 added to the end, so that a concatenation series 304 “A&T$” will be received. Likewise, a chaining series will be obtained. 306 “C&TA$” from a pair of the sample DNA fragment 305 receive.

[0058] Furthermore, several different types of concatenation symbols can be used according to the lengths of the sample DNA fragments. For example, "&" can be used for a concatenation symbol in the concatenation sequence generated from approximately 300 bases of a sample DNA fragment, while "#" can be used for a concatenation symbol in the concatenation sequence generated from approximately 20,000 bases of a sample DNA fragment. This makes it possible to obtain not only the mating partner of a sample DNA fragment but also the length of a sample DNA fragment from concatenation symbols.

[0059] In S103 in Fig. 3 introduces the reading list preparation unit 21 a BW transformation on the chaining series 304 and 306 through to a BW series 311 to produce.

[0060] In Fig. 4B will be the BW series 311generated, for example, by the following procedure. If, in this calculation process, a comparison is made between a character from the concatenation series 304 and a sign from the concatenation series 306 is performed sequentially from the head when a comparison between “$” from the concatenation series 304 and “$” from the concatenation series 306 A comparison is performed, and if a comparison between “&” from the concatenation series is carried out, the comparison is terminated. 304 and “&” from the concatenation series 306 When comparing them, the comparison continues.

[0061] (Process 1) A sequence list 307 is achieved by performing a cyclic shift of the chain series 304 received and a list of orders 308 is achieved by performing a cyclic shift of the chain series 306 receive.

[0062] (Process 2) A combined list 309is achieved by combining two lists of 307 and 308 receive.

[0063] (Process 3) A sorted list 310 sorting the combined list 309 The characters are preserved in alphabetical order. The sorting rank of characters is, for example, “$ < # < & < A < C < G ​​< T < N” in this case.

[0064] (Process 4) A BW series 311 is achieved by concatenating the end characters of each line of the sorted list 310 educated.

[0065] Since the BW series thus obtained 311 If the data has already been sorted, the same characters often continue. However, the amount of data can be reduced by performing run-length encoding on the BW series. 311 be compressed.

[0066] In S105 in Fig. 3 sets the reading list preparation unit 21 the BW series 311 to the wavelet tree format to convert the read sequence directory 14to prepare for efficient recovery.

[0067] Fig. 5A represents a binary tree 320 that, on which to implement the BW series 311 This refers to the wavelet tree format. This binary tree 320 denotes all characters ($, &, A, C, G, T, N) 321 , which are used for the series, as root nodes.

[0068] The binary tree 320 describes the subdivision procedure of recursive subdivision of all characters ($, &, A, C, G, T, N) 321 , which are used for the series, into 2 groups, so that at most 2 types of characters are included at the end of the subdivision.

[0069] At the root node of the binary tree 320 All characters will be 321 in A, T (W) 24 and the others (S) 325 subdivided. The characters ($, &, C, G, N) 3250Characters that are subdivided into S are likewise subdivided into M and K. This is repeated recursively until at most two types of characters remain at the end of the subdivision. The character “&” 302 , which indicates a pair, and the closing symbol “$” 303 However, they are subdivided in such a way that they are at one end 331 appear together in the subdivision.

[0070] A wavelet tree 340 in Fig. 5B is a binary tree that represents the BW series. 311 indicates the root node of the wavelet tree. 340 is achieved by transforming the BW series 311 into a binary series 341 according to two groups of W and S in the binary tree 320 received. The reading list preparation unit 21 generates a subseries 342 , which is formed by extracting the characters that are subdivided into W, and a subseries 343 , which is formed by extracting the characters that are subdivided into S.

[0071] The reading list preparation unit 21 continues the sub-series 342 (including only 2 types of characters A and T) into a binary series 344 um, which consists of 0 and 1.

[0072] The reading list preparation unit 21 The sub-series also sets 343 into a binary series 345 according to two groups of M and K formed by symbols subdivided into S (by the binary tree) 320 (specified).

[0073] The reading list preparation unit 21 generates a subseries 346 , which is formed by extracting the characters that are subdivided into M, and a subseries 347 , which is formed by extracting the characters that are subdivided into K.

[0074] By repeating the processing of the generation of each sub-series, as described above, 6 pieces of binary series are formed, as shown. At one end 351 The subdivision is a series 350 , which are formed from the character “&” 302 , which indicates a pair, and the closing symbol “$” 303 consists of a binary series 351 The length (number of bits) of the binary sequence is represented. 351 is equal to the total number of readings (2 times the total number of pairs).

[0075] Consequently, the wavelet tree 340 by performing an invertible transformation on the BW series 311 received, therefore the BW series can be obtained 311 from the wavelet tree 340 be reconstructed.

[0076] In S106 in Fig. 3 provides the reading list preparation unit 21 the reading sequence directory 14 , which comes from the wavelet tree 340The generated output is sent to a storage unit. Here, the read directory preparation unit can be used. 21 to the reading sequence directory 14 next to the wavelet tree 340 small (approximately 3.5% data set relative to the BW series) 311 ) Auxiliary data for calculating a rank function and a selection function of the wavelet tree 340 add.

[0077] Rank (p, c) is a function that returns a number of occurrences of a character “c” from sequence elements 0 and p.

[0078] Select (i, c) is a function that returns a sequence position where the (i + 1)th character “c” appears.

[0079] Auxiliary data is a “hierarchical binary string,” as described, for example, in the reference work “Kouichi Kimura, Yutaka Suzuki, Sumio Sugano, and Asako Koike, Computation of Rank and Select Functions and a Hierarchical Binary String and Its Application to Genome Mapping Problems for Short-Read DNA Sequences,” Journal of Computational Biology, November 2009, 16(11): 1601–1613.” This auxiliary data is used to efficiently perform the recovery of primarily all read sequence fragments corresponding to a base sequence arbitrarily chosen from the BW series. 311 is selected.

[0080] Since the auxiliary data used for the efficient calculation of a rank function and a selection function are provided at the sequence obtained by BW transformation, it is possible to efficiently find all read sequence fragments corresponding to a sequence s chosen arbitrarily, and furthermore to perform a reconstruction calculation on left and right sequences (each extending up to “$”) obtained from each read sequence fragment, using the properties of the BW transformation that two types of separators $ and & are used.

[0081] The reading list preparation steps (S101 to S105) were described above with reference to Fig. 3 to Fig. 5 described. Meanwhile, genome inventory preparation steps (S101b to S105b in Fig. 3) similar to the reading directory generation steps (S101 to S105).

[0082] In S101b, the genome inventory preparation unit receives... 22 the inputs from the reference genome sequence 12 instead of the reading sequence group 11 in S101.

[0083] In S102b, the genome inventory preparation unit is linked 22 several chromosome sequences (base rows) of a genome DNA, which are replaced by the reference genome sequence 12 They are represented as they are, with the ending symbol “$” to create a segment of a series. Since the reference genome sequence 12 Since the pair is not present here in paired form, pair concatenation processing with the concatenation symbol “&”, as shown in S102, is not required.

[0084] In S105b, the genome inventory preparation unit 22 the genome sequence directory 15 instead of the reading sequence directory 14 exit at S105.

[0085] Fig. 6 is a flowchart that represents the analysis processing of the SNP.

[0086] In S121 the sequence data analyzer receives 1 the SNP information inputs 13a as analysis information 13 , which represent the mutation to be analyzed this time.

[0087] Fig. 7A presents a table 400 that show the SNP information 13a shows. As in every row of the table 400 As shown, each SNP contains the information of a chromosome name, a base position coordinate on a chromosome, a type of base in the reference genome sequence (normal base), and a type of base that appears as a SNP (mutated base).

[0088] The first row of the table 400 shows that the SNP is located in the 123456th base on a chromosome 7 is arranged and a base “A” of the reference genome is mutated into a base “G”.

[0089] In S122 in Fig. 6 generates the query generation unit 23 the query sequence 16with SNP, which is determined by the SNP information 13a is indicated with reference to the genome sequence directory 15 .

[0090] A description field 420 in Fig. 7B represents S122, in which a horizontal axis 421 specifies a base position coordinate on a chromosome. First, the query generation unit refers to... 23 to the genome sequence directory 15 , to create a base sequence 422 in the series of (for example, approximately 10 bases to the left and to the right) an SNP position 424 to find the one in the first row of the table 400 The query generation unit is shown. 23 mutates the bases of the base sequence 422 in the SNP position 424 , to create a sequence 423 to generate a query sequence 16 is referred to as. Alternatively, the query generation unit is referred to as... 23the base sequence that does not include any mutation, as a query sequence 16 instead of the base sequence 423 with a mutation to detect a normal base located in the SNP position 424 appears.

[0091] In S123 in Fig. 6 wins the query recovery unit 24 the hit position 17a , which is a base position coordinate, within the read sequence directory 14 again, which refers to the query sequence 16 to be depicted (the orientations of the base rows are accordingly).

[0092] In the calculation process based on the read sequence directory 14 , which is implemented through BW transformation when a comparison is made between a character from the read sequence directory 14 and a character from the query sequence 16 is performed sequentially from the heads if “$” is from the read sequence directory. 14and “$” from the query sequence 16 Furthermore, if a comparison is made with each other, a comparison is terminated, and if “&” is from the read sequence directory 14 and “&” from the query sequence 16 If they are compared with each other, the comparison will continue.

[0093] In S131 to S139, the sequence data analyzer performs 1 a loop processing at each of the hit positions 17a from those found in S123.

[0094] In S132, the sample reconstruction unit reconstructs 25 the reading sequence (the sample sequence) 17 ) of the sample DNA fragment with the hit position 17a , which is currently selected in the loop of S131. In this reconstruction process, the read sequence directory is 14 , which is generated by BW transformation, by extending from the hit position 17aup to the end character “$” using a rank function and a selection function, so that the sample sequence 17 between the left and right end sign “$”.

[0095] One such method for using a rank function and a selection function is described, for example, in the reference “Ferragina, P. and Manzini, G., Opportunistic data structures with applications, In 41 st described in IEEE Symposium on Foundations of Computer Science, FOCS, pages 390–398.

[0096] The sample sequence 17 includes the concatenation symbol “&” and the end symbol “$” as concatenation series 304 and 306 This sample sequence 17 is separated by the concatenation symbol “&”, so that a pair of two reading sequences (the left sequence) 11a , the right sequence 11b ) without using an identifier for each sample DNA fragment.

[0097] In S133, the sample reconstruction unit probes (examines). 25 the sample sequence 17 up to the concatenation symbol “&”, as described below, to form the pairing sequence 17b to obtain (S134).

[0098] Since in the case of [1] in Fig. 7C the concatenation symbol “&” appears in the sequence that differs from the sequence 423 in the hit position 17a Extending to the right, is the obtained mating sequence 17b the right sequence 11b , which extends to the right side of the concatenation sign “&”.

[0099] Since in the case of [2] in Fig. C the concatenation symbol “&” does not appear in any sequence that differs from the sequence 423 in the hit position 17a extending to the right or to the left, is the obtained mating sequence 17b the left sequence 17b , which extends on the left side of the concatenation sign “&”.

[0100] If, in the case of [3] in Fig. C, the concatenation sign “&” is not extended to the right or to the left of the sequence 423 in the hit position 17a The sample sequence appears 17 An independent reading sequence without partners forming a pair. Such an independent reading sequence can be considered unreliable and is therefore ignored.

[0101] Alternatively, only if it is confirmed that there is only one position in which the sequence 423 appears in a genome, by querying the genome directory 15 after the sequence 423 Before introducing a mutation, it can be determined that an SNP has been detected (although this is unreliable).

[0102] In S135, the mapping unit refers to 26 to the genome sequence directory 15 , to indicate in the genome DNA the position (genome mapping) from which the mating sequence 17b, which is obtained in S134, is derived. The mapping unit 26 For example, it cuts a short subsequence (e.g., approximately 20 bases) from the pairing sequence. 17b and queries the genome directory 15 depending on whether the partial sequence appears only in one position in the genome.

[0103] If the subsequence does not appear in any position, it is assumed that the subsequence contains a sequence error or polymorphism; consequently, the mapping unit queries 26 the genome directory 15 after another sub-sequence.

[0104] If the subsequence appears in multiple positions, the mapping unit queries 26 Furthermore, the subsequence is repeated at an increased length or as a different subsequence. If the position of the short subsequence is thus determined on a genome, the position of the pairing sequence can be determined. 17bThe short subsequence must also be determined at the genome level (in other words, successful mapping). If, on the other hand, the position cannot be determined, the mapping unit fails. 26 the illustration.

[0105] If the determination of S135 is yes, the loop continues to S136, while if the determination of S135 is no, this loop ends (S139) to return to S131 to select the next hit position.

[0106] In S136, the sample determination unit is determined 27 , whether the sample DNA fragment of the sequence 17a S132 contains an SNP mutation, which is determined by the analysis information 13 The result is given depending on whether the distance is normal (correct) or not, as described below. The sample determination unit 27 For example, consider the distance between the image position of the mating sequence. 17b, which was successful in the mapping, and the SNP position determined by the SNP information 13a (the hit position 17a ) is specified, approximately equal to the length of the sample DNA fragment, that the left sequence 11a and the right sequence 11b , which the sample sequence 17 form, accordingly, so that they are determined as “SNP detected” (in other words, according to the SNP information) 13a ).

[0107] If the determination of S136 is yes, the loop continues to S137, while if the determination of S136 is no, the loop continues to S138.

[0108] In S137, the sample determination unit increases 27 the detection count counter value of the SNP information 13a , which are determined to be “SNP detected”, by one.

[0109] In S138, the sample determination unit increases 27not the detection count counter value due to the determination as “SNP not detected”.

[0110] During the loop from S131 to S139, the sample determination unit is given in S141. 27 (reports to a user) the result of the analysis information 13 from where the SNP information 13a are assigned to the detection count counter value (the number of sample DNA fragments in which SNP is detected).

[0111] A table 460 in Fig. 7D is an example of information output in S141. It represents the number of read fragments (detection count of mutated bases) in which a mutated base (SNP) is found at the SNP position. 424 is detected and the number of read fragments (number of normal bases detected) in which a normal base is detected for each SNP, which is listed in the table 400 in Fig. The values ​​shown in 7A are read from the detection count counter value and entered into the table. 460 written.

[0112] Fig. 8 is the flowchart that specifies the analysis processing of the structural mutation. Fig. Section 8 below focuses on the difference to the analysis processing of the SNP in Fig. 6 described.

[0113] In S121b, the sequence data analyzer receives 1 the input of the structure mutation information 13b instead of the SNP information 13a as analysis information 13 , which specify the mutation to be analyzed this time.

[0114] The structural mutation information 13b in Fig. 9A includes a chromosome name, a base position coordinate on a chromosome, a mutation type (insertion or deletion), and information about a mutation length, which is entered in each row of a table. 600 shown.

[0115] The first row of the table 600 indicates that a structural mutation is located at the position of base 654321 on the third chromosome and that the deletion, in which the orientation of 600 The loss of continuous bases from the normal genome is caused by this.

[0116] In S122b in Fig. 8 creates a query generation unit 23 the query sequence 16 in the series of a structural mutation, which is determined by the structural mutation information 13b is specified with reference to the genome sequence directory 15 .

[0117] A description field 620 in Fig. 9B describes the generation method of the query sequence 16 The horizontal axis is shown in S122b and is described. 421 Specifies a base position coordinate on a chromosome. First, with reference to the genome sequence directory. 15 short (for example, approximately 20 bases) base sequences 622and 623 in the row (for example, in position ten bases to the left and right) of the position 624 , in which the structural mutation is caused, found so that they can each be used as a query sequence 16 be defined. The query sequence 16 is a short sequence in the series (to the left or to the right) of the position 624 , in which a structural mutation is caused.

[0118] In S133b in Fig. 8 is executed after processing S133 (in Fig. (8 not shown) determines whether the chaining series of the sample sequence reconstructed from the hit position is excluded from the structural mutation analysis (mutation information is not retained) or not. If the determination of S133b is Yes, this loop (S139) ends, while if the determination of S133b is No, this loop continues to S134. The determination processing of S133b is described below with reference to the description field. 640 from Fig. 9C described.

[0119] As in the case of [1] in Fig. 9C, where the concatenation symbol “&” is in the sample sequence 17 appears, which is reconstructed in such a way that the query sequence 622 to the left of the position 624 , which determines the hit position 17a The query sequence extends to the right. 622 in the left sequence 11a contain.

[0120] Consequently, the sequence641 (the right sequence 11b ) to the right of the concatenation sign “&” as a pairing sequence 17b reconstructed. Furthermore, the position can be determined. 624 , in which a structural mutation is caused, between the query sequence 622 and the sequence 641 It should be included.

[0121] Although in Fig. 9C is not shown if the concatenation sign “&” appears in a series that is reconstructed in such a way that the query sequence 622 If it extends to the left, on the other hand there is no possibility of a structural mutation, consequently it is determined to be outside a target (mutation information is not preserved).

[0122] As in the case of [2] in Fig. 9(c) will make the similar determination on the query sequence 623 to the right of the position 624 , which determines the hit position 17a is, by horizontally rotating [1] into Fig. 9C performed. If the concatenation symbol “&” is in the sample sequence 17 appears, which is reconstructed in such a way that the query sequence 623 extending to the left, the query sequence 623 into the right sequence 11b recorded.

[0123] Consequently, the sequence 642 (the left sequence 11a ), which extends to the left from the concatenation sign “&”, as a pairing sequence 17b reconstructed. Furthermore, the position is 624 , in which a structural mutation is caused, probably between the query sequence 623 and the sequence 624 recorded.

[0124] Although in Fig. 9C is not shown if the concatenation character “&” appears in the sequence that is reconstructed in such a way that the query sequence 623If it extends to the right, on the other hand there is no possibility of a structural mutation, consequently it is determined to be outside a target (mutation information is not preserved).

[0125] The above describes the determination processing of S133b.

[0126] If in S136b, which is performed instead of S136, the distance between the pair of read sequences is between the pairing sequence 17b and the sample sequence 17 , from which the mating sequence 17b The sample determination unit evaluates the extractable DNA fragments to be almost equal to the length of the sample DNA fragments. 27 , that the distance is a normal value (S136b, Yes).

[0127] If the S136b determination is negative, a structural mutation is detected; consequently, the sample determination unit is increased. 27 the counter value of the “number of detected mutations” of the corresponding structural mutation information 13b at one.

[0128] Here the sample determination unit can be used. 27 Sum the number of detections separately for a mutation type "deletion" and a mutation type "insertion". If the distance between a pair is longer than predicted by the length of a sample DNA fragment, it can therefore be determined that a deletion (corresponding to a length of the difference) is caused; conversely, if the distance between a pair is shorter than predicted by the length of a sample DNA fragment, it can be determined that an insertion (corresponding to a length of the difference) is caused.

[0129] If the S136b determination is positive, no structural mutation is detected; consequently, the sample determination unit is increased. 27 the counter value “of the number of undetected mutations” of the corresponding structural mutation information 13b by one (S138b).

[0130] In S141b, the sequence data analyzer reports 1the counter values ​​of S137b and S138b as detection results of structural mutations, so that the values ​​shown in a table 660 from Fig. The information shown in 9D will be displayed.

[0131] The table 660 reports for each one in the table 600 The structural mutation shown represents the number of specific reading fragments as the existence of a mutation and the number of specific reading fragments as the non-existence of a mutation (in other words, corresponding to a normal genome) as the number of detections of the existence of mutations or the number of detections of the non-existence of mutations.

[0132] In the embodiment described above, the reading directory preparation unit connects 21 the left sequence 11a and the right sequence 11b, which form a pair of a sample DNA fragment, with the concatenation character and joins pairs of each of the sample DNA fragments with the end characters to form a concatenation chain. And the reading directory preparation unit 21 performs a BW transformation on the chaining sequence and then prepares the read sequence directory. 14 before, which is formatted as a wavelet tree.

[0133] The sample reconstruction unit 25 can, if they find the mating partner (the mating sequence) 17b ) of the sample DNA fragment with the hit position 17a from the point of impact 17a the query sequence 16 in the reading sequence directory 14 reconstructed the mating sequence 17b based on the reading sequence directory 14 Reconstruct the embedded concatenation marker. Consequently, the identifier is not required for each sample DNA fragment, so the pairing sequence17b , which forms a pair, can be efficiently computed for any reading sequence.

[0134] Furthermore, the amount of data and processing that can be done with the sequence data analyzer 1 The number of issues that need to be addressed can be reduced by not using an identifier for each sample DNA fragment, as described below.

[0135] As a comparative example, if an identifier for each sample DNA fragment such as “@seq1”, which is represented by a character 361 in Fig. 4 is used if the read sequence group 11 , which is read from the sample DNA fragment, from 1 billion of the left sequence 11a and the right sequence 11b This would require 5 billion types of identifiers. Assuming that the data size of an identifier is 4 Bytes is required, which determines the data size of all identifiers. 4Gigabytes. If the identifier is used for each sample DNA fragment, it is also necessary to determine the pairing sequence. 17b , which forms the pair, with the identifier for the hit read sequence as the recovery key.

[0136] On the other hand, in this embodiment the data size of the control characters (concatenation characters, end characters) that are in the read sequence directory 14 Since 1 bit is embedded for a sample DNA fragment, the total data size of all sample DNA fragments is 1 billion bits (equivalent to 0.125 gigabytes), which is approximately 32 times smaller than the data size when using identifiers. Furthermore, since the hit read sequence and the mating sequence 17b , which form the pair, adjacent via the concatenation symbol in the read sequence directory 14 are aligned, the recovery load of the mating sequence can 17b be suppressed.

[0137] Furthermore, the present invention is not limited to the embodiment described above, but includes various variations. The embodiment described above, for example, is a detailed explanation to facilitate the description of the present invention and is not necessarily limited to the embodiment with all the described configurations.

[0138] Furthermore, it is possible to exchange some configurations of one embodiment for configurations of another embodiment, and it is also possible to add configurations of another embodiment to configurations of one embodiment.

[0139] Furthermore, it is possible to add, delete, or replace other configurations for some configurations of each embodiment. Here, part or all of each configuration, function, processing unit, and processing means as described above can be designed, for example, by an integrated circuit for implementation in hardware.

[0140] Furthermore, each of the configurations and functions as above can be implemented in software in such a way that a processor interprets and executes the program that implements each function.

[0141] The information, such as program, table and file for executing each function, can be stored in recording devices such as a memory, a hard drive and an SSD (semiconductor drive) or recording media such as an IC card (integrated circuit board), an SD card, a DVD (digital versatile disk).

[0142] Furthermore, control and information lines necessary for explanation are indicated; consequently, not all control and information lines for the product are necessarily shown. It can be assumed that almost all configurations are actually interconnected. Reference symbol list 1 Sequence data analyzer 8 Genome servers 9 Sequencing machine 11 Reading sequence group 11a left sequence 11b right sequence 12 Reference genome sequence 13 Analysis Information 13a SNP Information 13b structural mutation information 14 Reading sequence directory 15 Genome sequence directory 16 query sequence 17 Sample sequence 17a Hit position 17b mating sequence 21 Reading List Preparation Unit 22 Genome Directory Preparation Unit 23 Query generation unit 24 Query Recovery Unit 25 Sample Reconstruction Unit 26 imaging units 27 Sample determination unit

Claims

[1] Sequence data analyzer, which includes: a read directory preparation unit configured to generate, for each of the sample DNA fragments, a concatenation series obtained by joining a pair of a left sequence and a right sequence by a concatenation character, both obtained by sequencing each of the sample DNA fragments from the two ends of each of the sample DNA fragments, and to prepare a read sequence directory based on a series obtained by joining the concatenation series of the sample DNA fragments with end characters; a query recovery unit configured to start searching the read sequence directory for a hit position that is a base character coordinate at which a query sequence generated from a base row of a genome DNA appears; a sample reconstruction unit configured to start searching the read sequence directory for the sequence from the hit position to the end character located on either side of the hit position in order to recover the sequence as a sample sequence, to start searching the sample sequence from the hit position for the concatenation character located on either side of the hit position, and to extract as a pairing sequence from the retrieved concatenation character to the end character either on the left sequence or the right sequence where the hit position does not exist; and an imaging unit configured to search the genome DNA for a base sign coordinate of the base series in which the mating sequence appears. [2] Sequence data analyzer according to claim 1, further comprising: a query generation unit and a sample determination unit, wherein The query generation unit mutates a base character at a predetermined position below the base row of the genome DNA into a different base character to generate the query sequence with the mutated base character and a base row of the genome DNA that is positioned at both ends of the mutated base character, and The sample determination unit, when a base character coordinate of the pairing sequence is specified by the imaging unit, detects SNPs (single nucleotide polymorphisms) in the sample DNA fragment based on the length of the sample DNA fragment, which is determined by the character interval of the base character coordinate between the hit position and the pairing sequence. [3] Sequence data analyzer according to claim 1, further comprising: a query generation unit; and a sample determination unit, wherein The query generation unit generates the query sequence from the base sequence of the genome DNA, which is arranged in the sequence at a predetermined position below the base sequence of the genome DNA, and The sample determination unit, when the base sign coordinate of the mating sequence is determined by the imaging unit, detects a structural mutation in the sample DNA fragment based on the length of the sample DNA fragment, which is determined by the interval sign of the base sign coordinate between the hit position and the mating sequence. [4] Sequence data analyzer according to any one of claims 1 to 3, wherein the read directory preparation unit generates a BW series obtained by performing a BW transformation (Burrows-Wheeler transformation) on the series obtained by joining the concatenation rows of each of the sample DNA fragments with the end characters, and converts the BW series into the wavelet tree format to prepare the read sequence directory. [5] Sequence data analyzer according to any one of claims 1 to 3, wherein the read directory preparation unit generates the concatenation sequence using several types of concatenation symbols corresponding to a length of each of the sample DNA fragments. [6] Sequence data analyzer according to any one of claims 1 to 3, wherein the sample reconstruction unit excludes the hit position from the targets to be processed by the imaging unit if the concatenation sign does not exist in the series of hit positions. [7] DNA analysis system that includes: a sequence data analyzer according to any one of claims 1 to 6; and a sequencing machine configured to sequence the sample DNA fragments from both ends and to transfer a pair of read sequences, consisting of the left sequence and the right sequence both obtained as a result of sequencing, to the sequence data analyzer. [8] Sequence data analysis methods, wherein the sequence data analyzer comprises a read directory preparation unit, a query recovery unit, a sample reconstruction unit and an imaging unit, wherein The read directory preparation unit generates a concatenation sequence for each of the sample DNA fragments, obtained by joining a pair of left-hand and right-hand sequences with a concatenation character, both obtained by sequencing the two ends of each of the sample DNA fragments, and prepares a read sequence directory based on a sequence obtained by joining the concatenation sequences of the sample DNA fragments with end characters. The query recovery unit starts searching the read sequence directory for a hit position that is a base character coordinate where a query sequence generated from a base series of a genome DNA appears, The sample reconstruction unit starts searching the read sequence directory for the sequence from the hit position to the end character located on both sides of the hit position in order to recover the sequence as a sample sequence; it then searches the sample sequence for the concatenation character located on both sides of the hit position, starting from the hit position, and extracts the pairing sequence from the recovered concatenation character to the end character either on the left sequence or the right sequence where the hit position does not exist. The imaging unit searches the genome DNA for a base sign coordinate of the base series in which the mating sequence appears.