WO2022054178A1 - 個体ゲノムの構造変異検出方法及び装置 - Google Patents

個体ゲノムの構造変異検出方法及び装置 Download PDF

Info

Publication number
WO2022054178A1
WO2022054178A1 PCT/JP2020/034166 JP2020034166W WO2022054178A1 WO 2022054178 A1 WO2022054178 A1 WO 2022054178A1 JP 2020034166 W JP2020034166 W JP 2020034166W WO 2022054178 A1 WO2022054178 A1 WO 2022054178A1
Authority
WO
WIPO (PCT)
Prior art keywords
sequence data
mapping
processor
sequence
read
Prior art date
Application number
PCT/JP2020/034166
Other languages
English (en)
French (fr)
Inventor
宏一 木村
Original Assignee
株式会社日立ハイテク
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社日立ハイテク filed Critical 株式会社日立ハイテク
Priority to PCT/JP2020/034166 priority Critical patent/WO2022054178A1/ja
Publication of WO2022054178A1 publication Critical patent/WO2022054178A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search

Definitions

  • the present disclosure relates to a method and an apparatus for detecting structural variation of an individual genome.
  • the most typical method for detecting structural variation is the method by pair-end sequencing.
  • the personal genome is usually fragmented into a large number of sequences (called inserts) having a length of about several hundred bases, the sequences at both ends of each insert are read by a sequencer, and a read sequence having a length of about 100 bases is read. Get a lot of pairs.
  • the position (mapping position) in the standard genome corresponding to the read sequence can be uniquely specified (Non-Patent Document 2).
  • the mapping positions of the paired read sequences are separated from each other by a distance approximately equal to the mean insert length on the standard genome sequence coordinates. Such a pair is called an accordant pair.
  • the distance between the mapping positions of the paired read sequences deviates from the average value of the insert length, it is presumed that some structural variation has occurred. Such pairs are called discordant pairs.
  • the distance between the mapping positions of the paired read sequences is smaller than the average value of the insert length, it is presumed that the decrease reflects the length of the insertion mutation that occurred in the individual genome. Therefore, an insertion mutation with a length exceeding the insert length cannot be detected. If the mapping of one or both of the paired read sequences fails, the data of that pair will not be used.
  • the present disclosure provides a technique for detecting an insertion mutation having a length exceeding the insert length by using the sequence data of a short-read type next-generation DNA sequencer by pair-end sequencing.
  • the method for detecting structural mutations in an individual genome of the present disclosure is a method for detecting structural mutations in an individual genome executed by a computer processor, and the standard genome sequence data is received by the processor.
  • the processor calculates the probability that the mapping will fail at each position of the standard genome sequence data, assuming that there is no mutation, and the processor receives the read sequence data of the individual genome. , The processor failed to map each end of the read sequence data onto the standard genomic sequence data, and the processor failed to map the end of the read sequence data at each position of the standard genomic sequence data.
  • the frequency is calculated, and the processor determines whether the frequency of failure of the mapping is significantly higher than the probability of failure of the mapping at each position of the standard genomic sequence data.
  • the processor extracts a region on the standard genome sequence data in which a position determined to have a significantly high frequency of failure of the mapping continues for a predetermined length or longer as a candidate for a mutation region, and the processor extracts the region. It includes outputting the mutation region candidate as a detection result of a structural mutation.
  • the figure explaining the method of calculating the read frequency with each base coordinate of a standard genome sequence The figure explaining the method of calculating the read frequency with each base coordinate of a standard genome sequence.
  • the figure explaining the method of extracting the candidate of a mutation region on the standard genome sequence The figure explaining the method of extracting the candidate of a mutation region on the standard genome sequence.
  • the figure explaining the method of detecting a breakpoint on a standard genome sequence The functional block diagram of the structural variation detection system which concerns on 2nd Embodiment.
  • the flowchart which shows the whole processing of the structural variation detection system which concerns on 2nd Embodiment.
  • FIG. 1 is a hardware configuration diagram of a structural variation detection system 1 for an individual genome (individual genome) according to the first embodiment.
  • the structural variation detection system 1 includes a computer 100, a database for storing human standard genome sequence data 111, and a database for storing personal genome read sequence data 112.
  • the computer 100 is a device such as a server having a normal computer configuration.
  • the computer 100 includes a CPU 101 (processor), a memory 102, a storage device 103, a network interface (NIF) 104, an input device 105, a display / output device 106, and a bus 107. Each component of the computer 100 is connected to each other by a bus 107.
  • the CPU (Central Processing Unit) 101 reads out a program temporarily stored in the memory 102 and various data, and executes processing necessary for detecting structural variation of the personal genome.
  • another processing device such as an MPU (MicroProcessingUnit) may be used.
  • the storage device 103 stores the standard genome sequence dictionary data 121, the standard genome MLRU (Minimum Length for Robust Uniqueness) data 122, and the personal genome read sequence dictionary data 123 generated by the processing of the CPU 101.
  • a hard disk drive (HDD), a solid state drive (SSD), a magnetic disk, an optical disk, or the like can be used as the storage device 103.
  • the standard genome sequence dictionary data 121, the standard genome MLRU data 122, and the personal genome read sequence dictionary data 123 may be stored in a storage device externally connected to the computer 100, or may be connected to the computer 100 via a network. It may be stored in a data center or the like.
  • the network interface 104 communicates with an external device of the computer 100 via a network such as a LAN (Local Area Network) and the Internet.
  • the CPU 101 can access and download the standard genome sequence data 111 and the personal genome read sequence data 112 stored in an external database via the network interface 104. Each data acquired from the outside is stored in the storage device 103.
  • the standard genome sequence data 111 is data of a human genome sequence defined as a reference, such as an international reference genome or a Japanese reference genome sequence.
  • the personal genome read sequence data 112 is, for example, a set of sequence data of fragments of the personal genome read by a short-read next-generation DNA sequencer by pair-end sequencing.
  • the input device 105 is, for example, a mouse, a keyboard, a touch panel, a camera, a microphone, or the like.
  • the display / output device 106 is, for example, a display, a touch panel, a printer, a speaker, or the like.
  • the display / output device 106 displays a GUI (Graphical User Interface) 108 for operation by the user, an analysis result, and the like on the display.
  • the user of the computer 100 can input information such as commands and parameters via the GUI 108 by operating the input device 105.
  • the input commands and parameters are stored in the memory 102 or the storage device 103.
  • FIG. 2 is a functional block diagram of the structural variation detection system.
  • the structural variation detection system includes a standard genome sequence dictionary creation unit 201, an MLRU calculation unit 202, an individual genome read sequence dictionary creation unit 203, an input unit 205, a display / output unit 206, and a pair end mapping unit 211.
  • Pair classification unit 212 mutation region candidate extraction unit 213, read sequence extraction unit 214, alignment unit 215, breakpoint (BP) extraction unit 216, mapping failure frequency evaluation unit 221 and mapping failure probability evaluation unit 222.
  • BP breakpoint
  • the input unit 205 stores the processing parameter 218 used in each unit of the CPU 101 designated (input) by the user via the GUI 108, and various input / output data (standard genome sequence data 111, personal genome read) described later.
  • the storage destinations of the sequence data 112, the standard genome sequence dictionary data 121, the standard genome MLRU data 122, the personal genome read sequence dictionary data 123, and the structural mutation detection result 113) are read and stored in the memory 102. Further, the input unit 205 acquires the standard genome sequence data 111 and the personal genome read sequence data 112 from an external database via the network interface 104, and stores them in the storage device 103 according to the designated storage destination.
  • the standard genome sequence dictionary creation unit 201 receives the input of the standard genome sequence data 111 from the storage device 103 and generates the standard genome sequence dictionary data 121.
  • the MLRU calculation unit 202 receives the input of the standard genome sequence data 111 from the storage device 103 and generates the standard genome MLRU data 122.
  • the personal genome read sequence dictionary creation unit 203 receives the input of the personal genome read sequence data 112 from the storage device 103, and generates the personal genome read sequence dictionary data 123. These generated data are stored in the storage device 103 according to the designated storage destination.
  • the pair end mapping unit 211 receives the input of the individual genome read sequence data 112, the standard genome sequence dictionary data 121, and the standard genome MLRU data 122, performs the pair end mapping process, and outputs the result to the pair classification unit 212.
  • the pair classification unit 212 accepts the input of the processing parameter 218 input by the user using the GUI 108 and the input device 105, sets the classification standard, and classifies the pair based on the result of the pair end mapping unit 211.
  • the pair classification unit 212 outputs the pair classification result to the mapping failure frequency evaluation unit 221.
  • the mapping failure frequency evaluation unit 221 calculates (evaluates) the frequency of actual mapping failures based on the pair classification result.
  • the mapping failure frequency evaluation unit 221 outputs the calculation result to the mutation region candidate extraction unit 213.
  • the mapping failure probability evaluation unit 222 calculates (evaluates) the probability of mapping failure based on the base reading error rate, which is one of the processing parameters 218, and the probability model using the standard genome MLRU data 122.
  • the mapping failure probability evaluation unit 222 outputs the calculation result to the mutation region candidate extraction unit 213.
  • the mutation region candidate extraction unit 213 receives the calculation result (frequency of mapping failure) of the mapping failure frequency evaluation unit 221 and also receives the calculation result (probability of mapping failure) of the mapping failure probability evaluation unit 222.
  • the mutation region candidate extraction unit 213 extracts a region in which the frequency of mapping failure is significantly higher than the probability of mapping failure under the condition of error rate, which is one of the processing parameters 218, as a mutation region candidate.
  • the result is output to the read sequence extraction unit 214.
  • the read sequence extraction unit 214 receives the input of the standard genome MLRU data 122 from the storage device 103, and determines the length of the seed sequence, which is a partial sequence used for extracting the read sequence. Next, the read sequence extraction unit 214 receives the input of the standard genome sequence data 111 and extracts the seed sequence (partial sequence) from the input. Next, the read sequence extraction unit 214 receives the input of the personal genome read sequence dictionary data 123 from the storage device 103, extracts all the read sequences including the seed sequence, and outputs the result to the alignment unit 215.
  • the alignment unit 215 receives the input of the standard genome sequence data 111 from the storage device 103, and aligns the read sequence extracted by the read sequence extraction unit 214 with the seed sequence commonly included in the standard genome sequence. Next, the alignment unit 215 extends the alignment to the region adjacent to the seed array, and outputs the result to the breakpoint extraction unit 216.
  • the breakpoint extraction unit 216 accepts the input of the processing parameter 218 and sets the alignment criterion. Next, the breakpoint extraction unit 216 extracts the boundary between the seed array and the adjacent region as a breakpoint when the alignment of the adjacent region of the seed array obtained by the alignment unit 215 does not meet the matching criterion, and displays the result. / Output to the output device 106.
  • the display / output unit 206 displays the breakpoint extraction result (structural variation detection result 113) obtained by the breakpoint extraction unit 216 as the GUI 108 on the display / output device 106. Further, the display / output unit 206 outputs the structural variation detection result 113 to the external device via the network interface 104 (not shown in FIG. 2).
  • FIG. 3 is a flowchart showing the entire process of the structural variation detection system. Details of some of the steps shown in FIG. 3 will be described later with reference to FIGS. 4-9.
  • the input unit 205 stores the processing parameter 218 used in each unit of the CPU 101 designated (input) by the user via the GUI 108, and various input / output data (standard genome sequence data 111, personal genome read sequence data).
  • the input of the storage destination of 112, the standard genome sequence dictionary data 121, the standard genome MLRU data 122, the personal genome read sequence dictionary data 123, and the structural mutation detection result 113) is received and stored in the memory 102.
  • Step 302 The input unit 205 acquires (receives) the standard genome sequence data 111 and the personal genome read sequence data 112 from the external database via the network interface 104, and stores them in the storage device 103 according to the storage destination specified in step 301. do.
  • the standard genome sequence data 111 is composed of a base sequence of a plus chain of each chromosome (a character string consisting of letters A, C, G, and T representing a base).
  • the standard genome sequence dictionary creation unit 201 reads the standard genome sequence data 111 from the storage device 103, concatenates these sequences and their reverse complementary sequences with the delimiter $, and forms a single standard genome sequence G. do.
  • the position of each character in the standard genome sequence G is specified by the coordinate X of an integer with the left end of the standard genome sequence G as 1.
  • the coordinate X specifies a character in the plus chain
  • the coordinate with a plus sign + X is used, and the base position in the complementary chain sequence (minus chain) corresponding to X is used.
  • Negatively signed coordinates-X are used to specify.
  • the standard genome sequence dictionary creating unit 201 calculates the Burrows-Wheeler transformation (FM index) of the standard genome sequence G using a known method (Non-Patent Document 5), and outputs the result as the standard genome sequence dictionary data 121. do.
  • the standard genome sequence dictionary creation unit 201 stores the generated standard genome sequence dictionary data 121 in the storage device 103.
  • the personal genome read sequence dictionary creation unit 203 reads the personal genome read sequence data 112 from the storage device 103, and uses known methods (Patent Document 1, Non-Patent Document 6) to Burrows-Wheeler of the personal genome read sequence data 112. The conversion (FM index) is calculated, and the result is output as the personal genome read sequence dictionary data 123.
  • the personal genome read sequence dictionary creation unit 203 stores the generated personal genome read sequence dictionary data 123 in the storage device 103.
  • the MLRU calculation unit 202 reads the standard genome sequence data 111 from the storage device 103 and calculates the standard genome MLRU data 122.
  • the MLRU calculation unit 202 stores the generated standard genome MLRU data 122 in the storage device 103.
  • the pair end mapping unit 211 reads the personal genome read sequence data 112, the standard genome sequence dictionary data 121, and the standard genome MLRU data 122 from the storage device 103, performs the pair end mapping process, and obtains the mapping positions at both ends of each pair.
  • the pair end mapping unit 211 outputs the mapping positions at both ends of each obtained pair to the pair classification unit 212.
  • the pair classification unit 212 sets each pair into a matching pair AP (acccordant pair), an inconsistent pair DP (discordant pair), and an incompletely mapped pair IP (incompletely mapped pair) based on the mapping position obtained by the pair end mapping unit 211. And, it is classified into one of the unmapping pair UP (unmappable pair).
  • the pair classification unit 212 outputs the pair classification result to the mapping failure frequency evaluation unit 221.
  • the mapping failure frequency evaluation unit 221 obtains the frequency of the matching pair AP and the frequency of the incomplete map pair IP at each base position of the standard genome sequence G based on the classification result of the pair.
  • the mapping failure frequency evaluation unit 221 calculates the frequency at which mapping actually fails as the relative frequency of the incomplete map pair IP to the total frequency of the matching pair AP and the incomplete map pair IP. Further, the mapping failure probability evaluation unit 222 calculates the probability of mapping failure based on the probability model of the binomial distribution using the base reading error rate, which is one of the processing parameters 218, and the standard genome MLRU data 122. do.
  • the mutation region candidate extraction unit 213 receives input of the frequency of actual mapping failure from the mapping failure frequency evaluation unit 221 and input of the probability of mapping failure from the mapping failure probability evaluation unit 222. Then, the mutation region candidate extraction unit 213 scans on the standard genome sequence G, and one region as a mutation region candidate is a region in which the frequency of actual mapping failure is significantly higher than the probability of mapping failure. Extract one by one.
  • Step 309 The mutation region candidate extraction unit 213 determines whether or not the mutation region candidate has been extracted. If there is no mutation region candidate (No), the entire process is terminated. At this time, the mutation region candidate extraction unit 213 may output a result indicating that there is no mutation region candidate to the display / output unit 206, and the display / output unit 206 indicates that there is no mutation region candidate. The result may be displayed as a GUI 108 on the display / output device 106 (display). If a mutation region candidate is obtained (Yes), the process proceeds to step 310.
  • Step 310 The mutation region candidate extraction unit 213 repeats the following processing with the coordinate X as the X coordinate at the right end of the extracted mutation region candidate.
  • Step 311 The mutation region candidate extraction unit 213 determines whether or not the coordinate X is within the mutation region candidate. If X is not among the mutation region candidates (No), the process returns to step 308 (mutation region candidate extraction process) to search for the next candidate. If the coordinate X is within the mutation region candidate (Yes), the process proceeds to step 312.
  • the read sequence extraction unit 214 obtains the seed sequence S having the coordinate X at the left end, and extracts all the read sequences including the seed sequence S from the personal genome read sequence dictionary data 123.
  • the seed sequence S having the coordinate X as the left end is a partial sequence of length MLRU (+ X) extracted from the standard genome sequence data 111 starting from the element at the position of the coordinate X.
  • a known method LF mapping, last-to-first column mapping
  • Step 313 The read sequence extraction unit 214 determines whether or not the read sequence extracted in step 312 exists. If there is no read sequence (No), the process proceeds to step 318. When there is a read sequence (Yes), the read sequence extraction unit 214 outputs the extracted read sequence to the alignment unit 215, and the process proceeds to step 314.
  • Step 3128 The mutation region candidate extraction unit 213 subtracts 1 from the coordinate X and returns to step 311.
  • the alignment unit 215 aligns the extracted read sequence and the standard genome sequence G at the left end X of the seed sequence S commonly contained therein, and extends the alignment to the left of X.
  • the alignment unit 215 outputs the result of the extension of the alignment to the breakpoint extraction unit 216.
  • Step 315) The breakpoint extraction unit 216 determines whether or not base mismatch occurs frequently in the alignment extended to the left of the coordinate X. If the discrepancy does not occur frequently (No), the process proceeds to step 318, and the mutation region candidate extraction unit 213 subtracts 1 from the coordinate X and returns to step 311. If the mismatch occurs frequently (Yes), the process proceeds to step 316.
  • Step 316 The breakpoint extraction unit 216 determines that there is a breakpoint at the coordinate X, and outputs the determination result to the display / output unit 206.
  • Step 317) The display / output unit 206 outputs the coordinates X of the breakpoint to the display / output device 106.
  • the display / output unit 206 displays, for example, the alignment result of the coordinate X of the breakpoint and its surroundings, and the fluctuation of the frequency of the matching pair AP and the frequency of the incomplete map pair IP around the breakpoint as described later as GUI 108. / Can be displayed on the output device 106 (display) or output as a file.
  • step 309 After that, the process returns to step 308, and the same process is repeated.
  • FIG. 4 is a flowchart showing a method in which the MLRU calculation unit 202 calculates the standard genome MLRU data 122 from the standard genome sequence data 111 in step 304.
  • Step 401 The MLRU calculation unit 202 reads out the standard genome sequence data 111 from the storage device 103.
  • Step 402 The MLRU calculation unit 202 reads out the standard genome sequence data 111 from the storage device 103, and concatenates those sequences and their reverse complementary sequences with the delimiter $ to form one standard genome sequence G.
  • the MLRU calculation unit 202 calculates the PLCP (Permutated Longest Common Prefix Array) of the standard genome sequence G by a known method (Non-Patent Document 4).
  • MLRU (X) which represents minimum length for robust uniqueness, is a robust comparison condition that allows a difference of 1 base or less in the standard genome sequence G from the partial sequence of the standard genome sequence G starting from the coordinate X.
  • the minimum required to be unique below ie, its subsequence does not match any other subsequence of the standard genomic sequence G originating from a location other than X, except for a difference of at most one base). It is defined as representing the length of.
  • L MLU (X) +1 is set as the initial value of L. Used (described later in step 408).
  • the partial sequence of the length of MLU When determining the coordinates X in the standard genome sequence G by comparing the partial sequences, if the partial sequence of the length of MLU is used, it is affected by SNPs (Single Nucleotide Polymorphism) contained in a large number in the individual genome. The position of the coordinate X may be wrong. On the other hand, when a partial sequence of MLRU length is used, the conditions for sequence comparison are relaxed so as to allow a difference of 1 base or less, so that the coordinates X are robustly unaffected by SNP. Can be correctly determined.
  • SNPs Single Nucleotide Polymorphism
  • Step 405 The MLRU calculation unit 202 reads out the standard genome sequence dictionary data 121 from the storage device 103.
  • Step 406 The MLRU calculation unit 202 substitutes 1 which is the leftmost coordinate of the standard genome sequence G into the coordinate X.
  • Step 408 The MLRU calculation unit 202 sets MLU (X) + 1 as the initial value of the value L of the MLRU (X) at the coordinates X of the standard genome sequence G.
  • the MLRU calculation unit 202 uses S as a partial sequence of length L starting from the coordinate X (left end) in the standard genome sequence G.
  • the MLRU calculation unit 202 sets Q as a sequence in which a single base substitution, insertion or deletion is inserted in the partial sequence S.
  • the MLRU calculation unit 202 uses the standard genome sequence dictionary data 121 for each sequence Q by a known method (Non-Patent Document 2), and the number of times the sequence Q appears in the standard genome sequence G. Obtain Occ (G, Q). Then, the MLRU calculation unit 202 determines whether or not Occ (G, Q)> 0. If the value of Occ (G, Q) is positive (Yes), the process proceeds to step 412. On the other hand, if the value of Occ (G, Q) is 0 (No), the process proceeds to step 413.
  • Step 412 The MLRU calculation unit 202 increases the value of L by 1 and returns to step 409 for determining the partial array S.
  • Step 413 The MLRU calculation unit 202 determines whether or not the tests of all the sequences Q in the 8L manner described above have been completed. When the test of all the sequences Q is completed (Yes), the process proceeds to step 414. If the test of all the sequences Q has not been completed (No), the process returns to step 410 for defining the sequence Q, and the test process of step 411 is performed on the other sequences Q.
  • the MLRU calculation unit 202 defines MLRU (X), which is the value of MLRU at the coordinate X, as L, outputs it as standard genome MLRU data 122, and stores it in the storage device 103.
  • Step 415) The MLRU calculation unit 202 increments X by 1 and returns to the comparison process (step 407) with the coordinate Xmax at the right end of the standard genome sequence G.
  • X> Xmax in step 407 the process ends.
  • FIG. 5 is a flowchart showing a method in which the pair end mapping unit 211 obtains the mapping positions at both ends of each pair in step 305.
  • the pair end mapping unit 211 reads the standard genome sequence dictionary data 121 and the standard genome MLRU data 122 from the storage device 103.
  • Step 502 The pair end mapping unit 211 reads out the personal genome read sequence data 112 from the storage device 103.
  • the personal genome read sequence data 112 is a set of paired read sequences, and each pair consists of two read sequences. Both ends of the pair mean the 5'ends of each pair of read sequences.
  • the pair end mapping unit 211 determines whether or not there is an unprocessed read sequence. Processing is performed on all read sequences, and when there are no unprocessed read sequences (No), the processing ends. If there is an unprocessed read sequence (Yes), the process proceeds to step 504.
  • Step 504 The pair end mapping unit 211 sets the unprocessed read sequence as R and sets it as the next processing target.
  • Step 505 The pair end mapping unit 211 initializes the query length L to 1.
  • the pair end mapping unit 211 takes a partial sequence of length L from the 5'end of the read sequence R and uses it as the query sequence Q.
  • the pair end mapping unit 211 uses the standard genome sequence dictionary data 121 by a known method (Non-Patent Document 2 and Non-Patent Document 3) to obtain Occ (G) the number of times the query sequence Q appears in the standard genome sequence G. , Q), and let the value be F.
  • Step 508 The pair end mapping unit 211 determines whether or not F> 1. If F> 1, (Yes), the process proceeds to step 509. If F> 1 is not (No), the process proceeds to step 510.
  • Step 509 The pair end mapping unit 211 increases L by 1 and returns to step 506 for determining the query sequence Q.
  • Step 511 The pair end mapping unit 211 determines that the mapping of the read sequence R has failed, and outputs the determination result to the mutation region candidate extraction unit 213. After that, the process returns to step 503.
  • Step 512 The pair end mapping unit 211 obtains the coordinates X of the only appearance position of the query sequence Q by a known method (Non-Patent Document 2), and updates the value of L to MLRU (X).
  • Step 513 The pair end mapping unit 211 takes a partial sequence of length L from the 5'end of the read sequence R and uses it as the query sequence Q.
  • the pair end mapping unit 211 uses the standard genome sequence dictionary data 121 by a known method (Non-Patent Document 2 and Non-Patent Document 3) to obtain Occ (G) the number of times the query sequence Q appears in the standard genome sequence G. , Q), and let the value be F.
  • Step 516 The pair end mapping unit 211 determines that the mapping of the read sequence R is successful, and outputs the determination result and the mapping coordinate X (R) of the read sequence R to the mutation region candidate extraction unit 213. After that, the process returns to step 503, and when it is determined that there is no unprocessed read sequence, the process ends.
  • ⁇ Pair classification method> 6A-6F are diagrams for explaining how the pair classification unit 212 classifies each pair into one of the matched pair AP, the inconsistent pair DP, the incomplete map pair IP, and the unmapping pair UP in step 306. be. This classification is performed independently for each pair.
  • FIGS. 6A-F shows a typical case where a pair is classified as either AP, DP, IP or UP.
  • the horizontal axis 601 is a coordinate axis indicating a position (coordinates) on the standard genome sequence G.
  • the mapping coordinates X (R1) are abbreviated as X1.
  • the mapping coordinate X (R2) is abbreviated as X2.
  • the read sequence R1 indicated by the right-pointing arrow means that the end of the read sequence is mapped to the positive strand of the standard genomic sequence G, and the read sequence R1 indicated by the left-pointing arrow has the end of the read sequence. It means that it is mapped to the minus chain of the standard genome sequence G. The same applies to the read sequence R2 indicated by the right-pointing or left-pointing arrow.
  • the x mark 606 (FIG. 6E, FIG. 6F) represents a case where the mapping of the end of the read sequence R2 fails in the pair end mapping unit 211 and the mapping position cannot be determined.
  • the x mark 607 (FIG. 6F) represents a case where the mapping at the end of the read sequence R1 fails and the mapping position cannot be determined.
  • the average value (or representative value) of the insert length is set to M
  • the threshold value that can be tolerated when the magnitude of the deviation from the average value of the insert length is within the normal range is set to V.
  • the average value M and the threshold value V of the insert length are a part of the processing parameter 218, and are acquired by the input unit 205 through the GUI 108 and stored in the memory 102 or the storage device 103.
  • FIG. 6A shows a case where it is classified into a matching pair AP.
  • the mapping of the ends of the read sequences R1 and R2 was successful, the mapping destinations of the ends of the read sequences R1 and R2 were on the same chromosome, the read sequence R1 was mapped to the plus chain, and the read sequence was used.
  • mapping of the ends of the read sequences R1 and R2 is successful, if they cannot be classified as a matching pair AP, they are all classified as an inconsistent pair DP. For example, if the mapping of the ends of the read sequences R1 and R2 is successful and both mapping destinations are on different chromosomes, this pair is classified as an unmatched pair DP.
  • FIG. 6B shows the case of being classified into an inconsistent pair DP.
  • the mapping destinations of both are on the same chromosome, and both are mapped to the plus strand, this pair is converted to the unmatched pair DP. Classify.
  • the ends of the read sequences R1 and R2 are successfully mapped, both are mapped to the same chromosome, and both are mapped to the minus chain, the pair is classified as an unmatched pair DP. ..
  • FIG. 6C shows the case of being classified as an inconsistent pair DP.
  • the mapping of the ends of the read sequences R1 and R2 was successful, both mapping destinations are on the same chromosome, the read sequence R1 is mapped to the plus strand, and the read sequence R2 is mapped to the minus strand. And if the coordinate X1 is located to the right of the coordinate X2, this pair is classified as an inconsistent pair DP. Further, when the roles of the read sequences R1 and R2 are exchanged in FIG. 6C, they are also classified into the inconsistent pair DP.
  • FIG. 6D shows the case of being classified as an inconsistent pair DP.
  • the mapping of the ends of the read sequences R1 and R2 was successful, both mapping destinations are on the same chromosome, the read sequence R1 is mapped to the plus chain, and the read sequence R2 is mapped to the minus chain.
  • FIG. 6E shows the case of being classified as an incomplete map pair IP.
  • the pair is classified as an incomplete map pair IP.
  • the coordinate X1 is called the mapping position of the incomplete map pair IP. The same applies when the roles of the read sequences R1 and R2 are exchanged in FIG. 6E.
  • FIG. 6F shows a case where it is classified into a non-mapping pair UP. As shown in FIG. 6F, if the mapping of the ends of the read sequences R1 and R2 both fails, this pair is classified as an unmapping pair UP.
  • ⁇ Calculation method of read sequence frequency> 7A and 7B are diagrams illustrating a method of calculating the read sequence frequency (read frequency) at each coordinate X of the standard genome sequence G.
  • FIG. 7A shows the case where the frequency of the read sequence is determined on the plus strand of the standard genome sequence G.
  • the 3'end is mapped to XR + L-1.
  • L represents the length of the read sequence.
  • the read frequency at the coordinates X of the standard genome sequence G is given by the number P (X) of the read sequences satisfying XR ⁇ X ⁇ XR + L-1.
  • P (X) the number of the read sequences satisfying XR ⁇ X ⁇ XR + L-1.
  • all read sequences R are sorted in ascending order of the coordinates XR of the mapping position at the 5'end, as is usually done, and XR in ascending order for all the coordinates X.
  • the number of read sequences satisfying ⁇ X ⁇ XR + L-1 may be counted.
  • FIG. 7B shows the case where the frequency of the read sequence is determined on the minus chain of the standard genome sequence G.
  • the 3'end is mapped to XR-L + 1.
  • L represents the length of the read sequence.
  • the read frequency at the coordinates X of the standard genome sequence G is given by the number of read sequences M (X) satisfying XR-L + 1 ⁇ X ⁇ XR.
  • M (X) In order to calculate M (X) efficiently, all read sequences R are sorted in descending order of the coordinates XR of the mapping position at the 5'end, as is usually done, and XRs in descending order for all coordinates X.
  • the number of read sequences satisfying ⁇ L + 1 ⁇ X ⁇ XR may be counted.
  • the frequency of the matching pair AP at each base position of the standard genome sequence G is the read frequency limited to the read sequence R belonging to the pair classified into the matching pair AP.
  • the frequency of the incomplete map pair IP at each base position of the standard genome sequence G is the read frequency limited to the read sequence R belonging to the pair classified into the incomplete map pair IP.
  • the mutation region candidate extraction unit 213 obtains a read sequence belonging to a pair classified into a matching pair AP and a read sequence belonging to a pair classified into an incomplete map pair IP based on the classification result of the pair classification unit 212. Using the read frequency calculation method described using 7A and 7B, the frequency of matched pair APs and the frequency of incomplete map pair IPs at each base position of the standard genome sequence G are determined.
  • ⁇ Extraction method of mutation region candidates> 8A and 8B are diagrams illustrating a method of extracting a candidate of a region containing a mutation by the mutation region candidate extraction unit 213 in step 308.
  • FIG. 8A shows a situation in which candidates for regions that may contain homozygous mutations are selected.
  • the mutation region candidate extraction unit 213 has the coordinates X of the standard genome sequence G on the horizontal axis and the read frequency on the vertical axis. Plot the frequency 805.
  • the frequency of the incomplete map pair IP, the frequency of the matched pair AP, and the frequency of the inconsistent pair DP at the coordinate X are expressed as IP (X), AP (X), and DP (X), respectively.
  • the mutation region candidate extraction unit 213 scans on the X-axis and determines that either the relative frequency of the incomplete map pair IP or the relative frequency of the inconsistent pair DP is significantly higher based on the probability model described later. A region in which the position continues for a specified length (about the read sequence length L) or more is extracted, and the extracted region is used as a candidate for a region containing a mutation.
  • the region 813 in FIG. 8A is a candidate for a mutation region extracted when the relative frequency of the incomplete map pair IP exceeds a threshold value (for example, about 0.25). In the case of homozygous mutations, RIP (X) takes a value close to 1 at the maximum.
  • FIG. 8A does not show an example of a candidate for a mutation region in which the relative frequency of the unmatched pair DP exceeds the threshold value, but in such a region, it is not obtained by a known method (Non-Patent Document 1). Mutations can be detected using matched pairs.
  • FIG. 8B shows a situation in which candidates for regions that may contain heterozygous mutations are selected.
  • the horizontal axis is the coordinate X of the standard genome sequence G
  • the vertical axis is the read frequency
  • the frequency 823 of the incomplete map pair IP, the frequency 824 of the matched pair AP, and the frequency 825 of the unmatched pair DP are plotted.
  • the mutation region candidate extraction unit 213 scans on the X-axis and determines that either the relative frequency of the incomplete map pair IP or the relative frequency of the inconsistent pair DP is significantly higher based on the probability model described later. A region in which the position continues for a specified length (about the read sequence length L) or more is extracted, and the extracted region is used as a candidate for a region containing a mutation.
  • the region 833 in FIG. 8B is a candidate for a mutation region extracted when the relative frequency of the incomplete map pair IP exceeds a threshold value (for example, about 0.25). In the case of heterozygous mutations, RIP (X) takes up to a value close to 0.5.
  • FIG. 8B does not show an example of a candidate for a mutation region in which the relative frequency of the unmatched pair DP exceeds the threshold value, but in such a region, it is not obtained by a known method (Non-Patent Document 1). Mutations can be detected using matched pairs.
  • the above method for extracting candidates for regions that can contain homozygous or heterozygous mutations can be applied to the analysis of individual genomes of organisms that distinguish between males and females.
  • the frequency of pairs is plotted for two purposes.
  • the entire standard genomic sequence G is plotted and scanned on the X-axis to relative the incomplete map pair IP.
  • the frequency and the relative frequency of the inconsistent pair DP are compared with the threshold value and used to extract candidates for regions containing mutations.
  • the plot around the candidate of the extracted mutation region is displayed as GUI 108 on the display / output device 106, and is used to show the user the situation where the candidate is extracted.
  • FIG. 9 is a diagram illustrating a method of extracting a breakpoint (BP) associated with a structural variation from the inside of a mutation region candidate by the alignment unit 215 and the breakpoint extraction unit 216 in steps 314, 315 and 316.
  • FIG. 9 describes a method for investigating the positive chain. When investigating the minus chain, flip the left and right and use the same method. On the plus chain, the inside of the mutation region candidate is scanned one base at a time from left to right, and the following processing is repeated.
  • BP breakpoint
  • Step 901 The alignment unit 215 uses X as the coordinate being scanned.
  • Step 902 The alignment unit 215 acquires the MLRU value MLRU (+ X) at the coordinate X on the plus chain from the standard genome MLRU data 122.
  • the alignment unit 215 extracts a partial sequence having the leftmost coordinate X and a length MLRU (X) on the plus strand from the standard genome sequence data 111, and uses it as the seed sequence S.
  • Step 904 In the alignment unit 215, the coordinate adjacent to the left of the coordinate X during scanning is set to X-1.
  • Step 905 The alignment unit 215 acquires the value MLRU (-(X-1)) of MLRU at X-1 on the minus chain from the standard genome MLRU data 122.
  • the alignment unit 215 extracts a partial sequence having the rightmost coordinate X-1 and a length MLRU (-(X-1)) on the minus chain from the standard genome sequence data 111, and sets it as the standard extension sequence F0. ..
  • the alignment unit 215 obtains all the read sequences including the seed sequence S from the personal genome read sequence dictionary data 123 by a known method (Non-Patent Document 2 and Non-Patent Document 3).
  • the alignment unit 215 obtains all the partial sequences having a length of MLRU (-(X-1)) or more adjacent to the left side of the seed sequence S in the obtained read sequence, and uses these as individual extension sequences (FIG. 9). In the example of F1, F2, F3, F4 and F5).
  • the alignment unit 215 aligns the standard extension sequence F0 and the individual extension sequences F1 to F5 by a known method such as dynamic programming, and outputs the alignment result to the breakpoint extraction unit 216.
  • the breakpoint extraction unit 216 obtains the edit distance between each of the individual extension sequences F1 to F5 and the standard extension sequence F0.
  • the breakpoint extraction unit 216 states that the read sequence in which the ratio of the editing distance to the length MLRU (-(X-1)) of the standard extension sequence F0 is a certain value (for example, about 0.5) or more is derived from a structural variation. judge.
  • the breakpoint extraction unit 216 obtains the ratio determined to be derived from the structural variation in the individual extension sequences F1 to F5, and uses this as the mutation rate.
  • the read sequence 921 containing F1, F2 and F3 is determined to be derived from a structural variation
  • the read sequence 922 containing F4 and F5 is determined not to be derived from a structural variation.
  • the rate is 0.6.
  • the breakpoint extraction unit 216 compares the mutation rate with 1 (in the case of homozygous mutation) or 0.5 (in the case of heterozygous mutation) to determine the presence or absence of homozygous or heterozygous mutation.
  • the breakpoint extraction unit 216 has an individual extension sequence (F1, F2 and F3) contained in the read sequence determined to be derived from the structural variation of the individual genome. It is determined that it is a part of the insertion sequence generated by the structural variation. Further, the breakpoint extraction unit 216 determines that the breakpoint of the structural variation is at the boundary position between the seed sequence S and the standard extension sequence F0 on the standard genome sequence G.
  • the breakpoint extraction unit 216 outputs these determination results as the structural variation detection result 113, stores them in the storage device 103, and outputs them to the display / output unit 206.
  • the display / output unit 206 presents the structural variation detection result 113 and the alignment result of the seed sequence S and the extension sequence used for the determination to the user by displaying the alignment result as the GUI 108 on the display / output device 106.
  • mapping failure probability evaluation unit 222 calculates the probability of mapping failure based on the probability model, and the probability that the mutation region candidate extraction unit 213 fails in mapping is compared with the frequency at which mapping actually fails. A method of determining whether or not the frequency of actual mapping failures is significantly high will be described.
  • both the pair end mapping unit 211 and the read sequence extraction unit 214 use a partial sequence having a length equal to MLRU in the read sequence.
  • the typical length of the MLRU of the human genome is about 20-30.
  • the condition for successful mapping is that the entire subsequence exactly matches the subsequence in the standard genomic sequence, as described with reference to FIG. Therefore, if the read sequence contains a base reading error during DNA sequencing, the mapping will fail.
  • the number of base reading errors that occur in a partial sequence can be approximated by a binomial distribution, and the error rate per base is usually about 0.01 (1%). Therefore, if MLRU is m and the error rate per base is e, the probability p that one read sequence fails to map in the absence of a mutation is given by the following [Equation 1].
  • the mutation region candidate extraction unit 213 sets ⁇ as an error rate (one of the processing parameters 218 specified by the user through the GUI 108, usually about 1 to 5%), and the probability E (r / n) of failure in mapping. ) Is less than the error rate ⁇ , it is determined that the frequency of mapping failures is significantly higher than the probability of mapping failure (calculated based on the probability model). As a result, the probability that a candidate for an erroneous mutation region is extracted by determining that the frequency of mapping failure is significantly higher than the probability of mapping failure when the mutation does not actually exist is determined. It becomes ⁇ or less.
  • the method for detecting structural mutations in the individual genome (individual genome) in the first embodiment is executed by the CPU 101 (processor) of the computer 100, and the input unit 205 receives the standard genome sequence data 111.
  • the mapping failure probability evaluation unit 222 evaluates the probability that mapping will fail at each position (coordinate X) of the standard genome sequence data 111 (standard genome sequence G) when it is assumed that there is no mutation, and the input unit.
  • the 205 receives the read sequence data 112 of the individual genome, the paired end mapping unit 211 maps each end of the read sequence data 112 onto the standard genome sequence data 111, and the mapping failure frequency evaluation unit 221.
  • the probability of mapping failure is evaluated assuming that there is no mutation, the frequency of mapping failure, which is information that has not been used in the past, is investigated, and when there are too many failures to match the above probability.
  • the read sequence extraction unit 214 extracts a seed sequence S, which is a partial sequence of standard genomic sequence data, from a candidate for a mutation region, and reads including the seed sequence S.
  • the sequence data is extracted, the standard genomic sequence data and the read sequence data are aligned to the adjacent extended sequence by aligning the common seed sequence S by the alignment unit 215, and the extended sequence is aligned by the breakpoint extraction unit 216.
  • the breakpoint extraction unit 216 is determined when the ratio of the editing distance of the read sequence data and the standard genomic sequence data to the length of the extended sequence is equal to or more than a predetermined value, it is determined that the read sequence data is derived from the structural variation, and the breakpoint extraction unit 216 is determined.
  • the ratio (variation rate) of the read sequence data determined to be derived from the structural variation is compared with 1 in the case of homotype mutation and 0.5 in the case of heterotype mutation. Includes determining the presence or absence of heterozygous mutations.
  • FIG. 10 is a functional block diagram of the structural variation detection system 2 according to the second embodiment.
  • the structural variation detection system 2 of the present embodiment is obtained by extracting necessary components from the configuration (FIG. 2) of the structural variation detection system 1 of the first embodiment, and has the same configuration as that of the first embodiment.
  • the reference numeral is attached, and the description thereof will be omitted.
  • the mutation region candidate extraction unit 1013 is different from the mutation region candidate extraction unit 213 (extracting mutation region candidates one by one) of the first embodiment in that all the mutation region candidates are extracted at once.
  • FIG. 11 is a flowchart showing the entire process of the structural variation detection system 2 according to the second embodiment.
  • the same processing as in the first embodiment is designated by the same reference numerals, and the description thereof will be omitted.
  • Steps 301 and 302 The input unit 205 executes steps 301 and 302 as in the first embodiment.
  • the standard genome sequence dictionary creating unit 201 calculates the Burrows-Wheeler transformation (FM index) of the standard genome sequence G using a known method (Non-Patent Document 5), and uses the result as the standard genome sequence dictionary data 121. ..
  • the standard genome sequence dictionary creation unit 201 stores the generated standard genome sequence dictionary data 121 in the storage device 103.
  • Steps 304-307 are performed as in the first embodiment.
  • the mapping failure frequency evaluation unit 221 calculates the frequency at which mapping actually fails as the relative frequency of the incomplete map pair IP to the total frequency of the matching pair AP and the incomplete map pair IP. Further, the mapping failure probability evaluation unit 222 calculates the probability of mapping failure based on the probability model of the binomial distribution using the base reading error rate, which is one of the processing parameters 218, and the standard genome MLRU data 122. do.
  • the mutation region candidate extraction unit 1013 receives input of the frequency of actual mapping failure from the mapping failure frequency evaluation unit 221 and input of the probability of mapping failure from the mapping failure probability evaluation unit 222. Then, the mutation region candidate extraction unit 1013 scans on the standard genome sequence G and extracts all regions in which the frequency of actual mapping failure is significantly higher than the probability of mapping failure as mutation region candidates. do.
  • the mutation region candidate extraction unit 1013 outputs the extracted mutation region candidate as the structural variation detection result 113, and ends the process.
  • the mutation region candidate extraction unit 1013 may output the structural variation detection result 113 to the display / output unit 206 and display it as GUI 108.
  • the output result of the mutation region candidate can be used as follows.
  • a long-read DNA sequencer throughput is Structural variation analysis is performed by directly reading the structural variation sequence using a sequencer that can directly read a long sequence of bases of several kilobases to several tens of kilobases, although it is low.

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biophysics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

コンピュータのプロセッサにより実行される個体ゲノムの構造変異検出方法であって、標準ゲノム配列データを受信することと、標準ゲノム配列データの各位置において、変異がないと仮定した場合にマッピングに失敗する確率を計算することと、個体ゲノムのリード配列データを受信することと、リード配列データのそれぞれの末端を標準ゲノム配列データ上にマッピングすることと、標準ゲノム配列データの各位置でリード配列データの末端のマッピングに失敗した頻度を計算することと、標準ゲノム配列データの各位置において、マッピングに失敗した頻度がマッピングに失敗する確率と比較して有意に大きいか否かを判定することと、標準ゲノム配列データ上で、マッピングに失敗した頻度が有意に大きいと判定される位置が所定の長さ以上続く領域を、変異領域候補として抽出することと、変異領域候補を構造変異の検出結果として出力することと、を含む。

Description

個体ゲノムの構造変異検出方法及び装置
　本開示は、個体ゲノムの構造変異検出方法及び装置に関する。
　ショートリード型の次世代型ＤＮＡシーケンシング技術の進歩により、ＤＮＡシーケンサで読み取った個人のゲノム配列を癌や生活習慣病などの予防及び治療に役立てることが可能になってきている。その際、個人のゲノム配列を標準ゲノム配列と比較して、標準ゲノム内のどの位置に、個人ゲノム側でどのような変異が生じているかを調べることが必要になる。変異には、１塩基～数塩基のスケールで生じるものから数キロ～数十キロの大きなスケールで生じるものまで様々ある。通常、ショートリード型の次世代型ＤＮＡシーケンサで直接読み取れる塩基長は１００～１５０塩基程度に限られるため、それを超える大きなスケールで生じる変異（構造変異とよぶ）は間接的な方法を用いて検出される。
　構造変異を検出するための最も代表的な方法が、ペアエンド・シーケンシングによる方法である。この方法では、通常、個人ゲノムを数百塩基程度の長さの多数の配列（インサートとよぶ）に断片化して、各インサートの両端の配列をシーケンサで読み取って、長さ１００塩基程度のリード配列のペアを多数得る。大多数の場合、リード配列に対応する標準ゲノム内の位置（マッピング位置）は一意に特定することができる（非特許文献２）。個人ゲノムに構造変異が生じていない場合、ペアをなすリード配列のマッピング位置は、標準ゲノム配列座標上でインサート長の平均値にほぼ等しい距離だけ互いに離れている。このようなペアは整合ペア（accordant pair）とよばれる。それに対して、ペアをなすリード配列のマッピング位置間の間隔がインサート長の平均値から外れている場合は、何らかの構造変異が生じていると推定される。このようなペアは非整合ペア（discordant pair）とよばれる。
米国特許第８７９８９３６号
Chen, K., Wallis, J., McLellan, M. et al. BreakDancer: an algorithm for high-resolution mapping of genomic structural variation. Nat Methods 6, 677-681 (2009). Li H. and Durbin R. (2009) Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 25:1754-60. Ferragina P, Manzini G. Proceedings of the 41st Symposium on Foundations of Computer Science (FOCS 2000). Los Alamitos, CA, USA: IEEE Computer Society; 2000. Opportunistic data structures with applications; p. 390-398. Karkkainen J., Manzini G., Puglisi S.J. (2009) Permuted Longest-Common-Prefix Array. In: Kucherov G., Ukkonen E. (eds) Combinatorial Pattern Matching. CPM 2009. Lecture Notes in Computer Science, vol 5577. Springer, Berlin, Heidelberg G. Nong, S. Zhang and W. H. Chan, "Linear Suffix Array Construction by Almost Pure Induced-Sorting," 2009 Data Compression Conference, Snowbird, UT, 2009, pp. 193-202, doi: 10.1109/DCC.2009.42. Heng Li, Fast construction of FM-index for long sequence reads, Bioinformatics, Volume 30, Issue 22, 15 November 2014, Pages 3274-3275
　ペアをなすリード配列のマッピング位置間の間隔がインサート長の平均値よりも小さい場合は、その減少分は個人ゲノムに生じた挿入変異の長さを反映していると推定される。そのため、インサート長を超える長さの挿入変異は検出することができない。また、ペアをなすリード配列の一方又は両方のマッピングに失敗した場合は、そのペアのデータは利用されない。
　そこで、本開示は、ペアエンド・シーケンシングによるショートリード型の次世代型ＤＮＡシーケンサの配列データを用いて、インサート長を超える長さの挿入変異を検出する技術を提供する。
　上記課題を解決するために、本開示の個体ゲノムの構造変異検出方法は、コンピュータのプロセッサにより実行される個体ゲノムの構造変異検出方法であって、前記プロセッサにより、標準ゲノム配列データを受信することと、前記プロセッサにより、前記標準ゲノム配列データの各位置において、変異がないと仮定した場合にマッピングに失敗する確率を計算することと、前記プロセッサにより、個体ゲノムのリード配列データを受信することと、前記プロセッサにより、前記リード配列データのそれぞれの末端を前記標準ゲノム配列データ上にマッピングすることと、前記プロセッサにより、前記標準ゲノム配列データの各位置で前記リード配列データの末端のマッピングに失敗した頻度を計算することと、前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記マッピングに失敗した頻度が前記マッピングに失敗する確率と比較して有意に大きいか否かを判定することと、前記プロセッサにより、前記標準ゲノム配列データ上で、前記マッピングに失敗した頻度が有意に大きいと判定される位置が所定の長さ以上続く領域を、変異領域候補として抽出することと、前記プロセッサにより、前記変異領域候補を構造変異の検出結果として出力することと、を含む。
　本開示に関連する更なる特徴は、本明細書の記述、添付図面から明らかになるものである。また、本開示の態様は、要素及び多様な要素の組み合わせ及び以降の詳細な記述と添付される特許請求の範囲の様態により達成され実現される。
　本明細書の記述は典型的な例示に過ぎず、本開示の特許請求の範囲又は適用例を如何なる意味に於いても限定するものではない。
　本開示の技術によれば、インサート長を超える長さの挿入変異を検出することができる。
　上記以外の課題、構成及び効果は、以下の実施の形態の説明により明らかにされる。
第１の実施形態に係る構造変異検出システムのハードウェア構成図。 第１の実施形態に係る構造変異検出システムの機能ブロック図。 第１の実施形態に係る構造変異検出システムの処理全体を示すフローチャート。 標準ゲノム配列データから標準ゲノムＭＬＲＵデータを計算する方法を示すフローチャート。 ペア末端を標準ゲノム配列上にマッピングする方法を示すフローチャート。 ペアを分類する方法を示す図。 ペアを分類する方法を示す図。 ペアを分類する方法を示す図。 ペアを分類する方法を示す図。 ペアを分類する方法を示す図。 ペアを分類する方法を示す図。 標準ゲノム配列の各塩基座標でリード頻度を計算する方法を説明する図。 標準ゲノム配列の各塩基座標でリード頻度を計算する方法を説明する図。 標準ゲノム配列上で変異領域の候補を抽出する方法を説明する図。 標準ゲノム配列上で変異領域の候補を抽出する方法を説明する図。 標準ゲノム配列上でブレークポイントを検出する方法を説明する図。 第２の実施形態に係る構造変異検出システムの機能ブロック図。 第２の実施形態に係る構造変異検出システムの処理全体を示すフローチャート。
　以下、図面を参照して本開示の実施形態を説明する。各実施形態において、ヒトのゲノム（個人ゲノム）を構造変異の検出対象として説明するが、本開示の技術は、他の生物のゲノム（個体ゲノム）の解析にも適用することができる。
［第１の実施形態］
＜構造変異検出システムの構成例＞
　図１は、第１の実施形態に係る個人ゲノム（個体ゲノム）の構造変異検出システム１のハードウェア構成図である。構造変異検出システム１は、コンピュータ１００、ヒトの標準ゲノム配列データ１１１を格納するデータベース及び個人ゲノムリード配列データ１１２を格納するデータベースを備える。
　コンピュータ１００（構造変異検出装置）は、通常の計算機の構成を有するサーバ等の装置である。コンピュータ１００は、ＣＰＵ１０１（プロセッサ）、メモリ１０２、記憶装置１０３、ネットワークインタフェース（ＮＩＦ）１０４、入力装置１０５、表示／出力装置１０６及びバス１０７を備える。コンピュータ１００の各構成要素は、バス１０７により互いに接続される。
　ＣＰＵ（Central Processing Unit）１０１は、メモリ１０２に一時的に記憶されたプログラム及び種々のデータを読み出して、個人ゲノムの構造変異検出に必要な処理を実行する。なお、ＣＰＵの代わりに、ＭＰＵ（Micro Processing Unit）など他の処理装置（プロセッサ）を用いてもよい。
　記憶装置１０３は、ＣＰＵ１０１の処理により生成された標準ゲノム配列辞書データ１２１、標準ゲノムＭＬＲＵ（Minimum Length for Robust Uniqueness）データ１２２及び個人ゲノムリード配列辞書データ１２３を格納する。記憶装置１０３としては、例えばハードディスクドライブ（ＨＤＤ）、ソリッドステートドライブ（ＳＳＤ）、磁気ディスク又は光学ディスクなどを用いることができる。なお、標準ゲノム配列辞書データ１２１、標準ゲノムＭＬＲＵデータ１２２及び個人ゲノムリード配列辞書データ１２３は、コンピュータ１００に外部接続された記憶装置に記憶されてもよいし、ネットワークを介してコンピュータ１００に接続されたデータセンタなどに記憶されていてもよい。
　ネットワークインタフェース１０４は、ＬＡＮ（Local Area Network）及びインターネットなどのネットワークを介してコンピュータ１００の外部装置と通信する。ＣＰＵ１０１は、ネットワークインタフェース１０４を介して外部のデータベースに格納される標準ゲノム配列データ１１１及び個人ゲノムリード配列データ１１２にアクセスし、これらをダウンロードすることができる。外部から取得された各データは、記憶装置１０３に格納される。標準ゲノム配列データ１１１は、例えば、国際基準ゲノム又は日本人基準ゲノム配列など、基準と定められたヒトのゲノム配列のデータである。個人ゲノムリード配列データ１１２は、例えば、ペアエンド・シーケンシングによるショートリード型の次世代型ＤＮＡシーケンサにより読み取られた個人ゲノムの断片の配列データの集合である。
　入力装置１０５は、例えばマウス、キーボード、タッチパネル、カメラ、マイクなどである。表示／出力装置１０６は、例えばディスプレイ、タッチパネル、プリンタ、スピーカなどである。表示／出力装置１０６は、ユーザによる操作のためのＧＵＩ（Graphical User Interface）１０８や解析結果などをディスプレイに表示する。コンピュータ１００のユーザは、入力装置１０５を操作することにより、ＧＵＩ１０８を介してコマンドやパラメータなどの情報を入力することができる。入力されたコマンドやパラメータは、メモリ１０２又は記憶装置１０３に格納される。
　図２は、構造変異検出システムの機能ブロック図である。図２に示すように、構造変異検出システムは、標準ゲノム配列辞書作成部２０１、ＭＬＲＵ計算部２０２、個人ゲノムリード配列辞書作成部２０３、入力部２０５、表示／出力部２０６、ペア末端マッピング部２１１、ペア分類部２１２、変異領域候補抽出部２１３、リード配列抽出部２１４、アラインメント部２１５、ブレークポイント（ＢＰ）抽出部２１６、マッピング失敗頻度評価部２２１及びマッピング失敗確率評価部２２２を有する。以上の各部の処理は、ＣＰＵ１０１がメモリ１０２に記憶されたプログラムを実行することにより実現される。
　入力部２０５は、ユーザによりＧＵＩ１０８を介して指定（入力）された、ＣＰＵ１０１の各部で用いられる処理パラメータ２１８の格納先、及び、後述する各種の入出力データ（標準ゲノム配列データ１１１、個人ゲノムリード配列データ１１２、標準ゲノム配列辞書データ１２１、標準ゲノムＭＬＲＵデータ１２２、個人ゲノムリード配列辞書データ１２３及び構造変異検出結果１１３）の格納先を読み込み、メモリ１０２に記憶させる。また、入力部２０５は、ネットワークインタフェース１０４を介して外部のデータベースから標準ゲノム配列データ１１１及び個人ゲノムリード配列データ１１２を取得し、指定された格納先に従って記憶装置１０３に格納する。
　標準ゲノム配列辞書作成部２０１は、記憶装置１０３から標準ゲノム配列データ１１１の入力を受け付け、標準ゲノム配列辞書データ１２１を生成する。ＭＬＲＵ計算部２０２は、記憶装置１０３から標準ゲノム配列データ１１１の入力を受け付け、標準ゲノムＭＬＲＵデータ１２２を生成する。個人ゲノムリード配列辞書作成部２０３は、記憶装置１０３から個人ゲノムリード配列データ１１２の入力を受け付け、個人ゲノムリード配列辞書データ１２３を生成する。これらの生成されたデータは、指定された格納先に従って記憶装置１０３に格納される。
　ペア末端マッピング部２１１は、個人ゲノムリード配列データ１１２、標準ゲノム配列辞書データ１２１及び標準ゲノムＭＬＲＵデータ１２２の入力を受け付けてペア末端マッピング処理を行い、その結果をペア分類部２１２に出力する。ペア分類部２１２は、ユーザがＧＵＩ１０８及び入力装置１０５を用いて入力した処理パラメータ２１８の入力を受け付けて分類基準を設定し、ペア末端マッピング部２１１の結果に基づいてペアの分類を行う。ペア分類部２１２は、ペアの分類結果をマッピング失敗頻度評価部２２１に出力する。
　マッピング失敗頻度評価部２２１は、ペアの分類結果に基づいて、実際にマッピングに失敗した頻度を計算する（評価する）。マッピング失敗頻度評価部２２１は、計算結果を変異領域候補抽出部２１３に出力する。
　マッピング失敗確率評価部２２２は、処理パラメータ２１８の一つである塩基読み取りエラー率と標準ゲノムＭＬＲＵデータ１２２とを用いた確率モデルとに基づいて、マッピングに失敗する確率を計算する（評価する）。マッピング失敗確率評価部２２２は、計算結果を変異領域候補抽出部２１３に出力する。
　変異領域候補抽出部２１３は、マッピング失敗頻度評価部２２１の計算結果（マッピングに失敗した頻度）を受け取るとともに、マッピング失敗確率評価部２２２の計算結果（マッピングに失敗する確率）も受け取る。変異領域候補抽出部２１３は、処理パラメータ２１８の一つである誤り率の条件の下で、マッピングに失敗した頻度がマッピングに失敗する確率よりも有意に大きい領域を変異領域の候補として抽出し、その結果をリード配列抽出部２１４に出力する。
　リード配列抽出部２１４は、記憶装置１０３から標準ゲノムＭＬＲＵデータ１２２の入力を受け付け、リード配列を抽出するために用いる部分配列であるシード配列の長さを決める。次に、リード配列抽出部２１４は、標準ゲノム配列データ１１１の入力を受け付け、その中からシード配列（部分配列）を取り出す。次に、リード配列抽出部２１４は、記憶装置１０３から個人ゲノムリード配列辞書データ１２３の入力を受け付け、シード配列を含む全てのリード配列を抽出し、その結果をアラインメント部２１５に出力する。
　アラインメント部２１５は、記憶装置１０３から標準ゲノム配列データ１１１の入力を受け付け、リード配列抽出部２１４で抽出されたリード配列と標準ゲノム配列に共通に含まれるシード配列を揃えてアラインメントを行う。次に、アラインメント部２１５は、アラインメントをシード配列に隣接する領域まで延長して、その結果をブレークポイント抽出部２１６に出力する。
　ブレークポイント抽出部２１６は、処理パラメータ２１８の入力を受け付けてアラインメントの一致基準を設定する。次に、ブレークポイント抽出部２１６は、アラインメント部２１５で得られたシード配列の隣接領域のアラインメントが一致基準を満たさない場合、シード配列と隣接領域の境界をブレークポイントとして抽出し、その結果を表示／出力装置１０６に出力する。
　表示／出力部２０６は、ブレークポイント抽出部２１６で得られたブレークポイントの抽出結果（構造変異検出結果１１３）をＧＵＩ１０８として表示／出力装置１０６に表示する。また、表示／出力部２０６は、ネットワークインタフェース１０４（図２には不図示）を介して構造変異検出結果１１３を外部装置に出力する。
＜構造変異検出方法＞
　図３は、構造変異検出システムの処理全体を示すフローチャートである。図３に示されるいくつかのステップの詳細については図４～図９を用いて後述する。
（ステップ３０１）
　入力部２０５は、ユーザによりＧＵＩ１０８を介して指定（入力）された、ＣＰＵ１０１の各部で用いられる処理パラメータ２１８の格納先、及び、各種の入出力データ（標準ゲノム配列データ１１１、個人ゲノムリード配列データ１１２、標準ゲノム配列辞書データ１２１、標準ゲノムＭＬＲＵデータ１２２、個人ゲノムリード配列辞書データ１２３及び構造変異検出結果１１３）の格納先の入力を受け付け、メモリ１０２に記憶させる。
（ステップ３０２）
　入力部２０５は、ネットワークインタフェース１０４を介して、外部データベースから標準ゲノム配列データ１１１及び個人ゲノムリード配列データ１１２をそれぞれ取得（受信）し、ステップ３０１で指定された格納先に従って、記憶装置１０３に格納する。
（ステップ３０３）
　ヒトゲノムＤＮＡは複数の染色体に分かれ、その各々はプラス鎖とマイナス鎖からなる２重らせん構造をもつ。標準ゲノム配列データ１１１は、各染色体のプラス鎖の塩基配列（塩基を表すＡ，Ｃ，Ｇ，Ｔの文字からなる文字列）からなる。標準ゲノム配列辞書作成部２０１は、記憶装置１０３から標準ゲノム配列データ１１１を読み出し、これらの配列とその相補鎖配列（reverse complementary sequences）を区切り文字＄で連結して一本の標準ゲノム配列Ｇとする。標準ゲノム配列Ｇ内の各文字の位置は、標準ゲノム配列Ｇの左端を１として整数の座標Ｘで指定される。また、説明上の必要に応じて、座標Ｘがプラス鎖内の文字を指定する場合、正符号付の座標＋Ｘを用い、また、Ｘに対応する相補鎖配列（マイナス鎖）内の塩基位置を指定するために負符号付の座標－Ｘを用いる。
　標準ゲノム配列辞書作成部２０１は、公知の方法（非特許文献５）を利用して標準ゲノム配列ＧのＢｕｒｒｏｗｓ－Ｗｈｅｅｌｅｒ変換（ＦＭインデクス）を計算し、その結果を標準ゲノム配列辞書データ１２１として出力する。標準ゲノム配列辞書作成部２０１は、生成した標準ゲノム配列辞書データ１２１を記憶装置１０３に格納する。
　個人ゲノムリード配列辞書作成部２０３は、記憶装置１０３から個人ゲノムリード配列データ１１２を読み出し、公知の方法（特許文献１、非特許文献６）を利用して個人ゲノムリード配列データ１１２のＢｕｒｒｏｗｓ－Ｗｈｅｅｌｅｒ変換（ＦＭインデクス）を計算し、その結果を個人ゲノムリード配列辞書データ１２３として出力する。個人ゲノムリード配列辞書作成部２０３は、生成した個人ゲノムリード配列辞書データ１２３を記憶装置１０３に格納する。
（ステップ３０４）
　ＭＬＲＵ計算部２０２は、記憶装置１０３から標準ゲノム配列データ１１１を読み出し、標準ゲノムＭＬＲＵデータ１２２を計算する。ＭＬＲＵ計算部２０２は、生成した標準ゲノムＭＬＲＵデータ１２２を記憶装置１０３に格納する。
（ステップ３０５）
　ペア末端マッピング部２１１は、記憶装置１０３から個人ゲノムリード配列データ１１２、標準ゲノム配列辞書データ１２１及び標準ゲノムＭＬＲＵデータ１２２を読み出し、ペア末端マッピング処理を行い、各ペアの両端のマッピング位置を求める。ペア末端マッピング部２１１は、求めた各ペアの両端のマッピング位置をペア分類部２１２に出力する。
（ステップ３０６）
　ペア分類部２１２は、ペア末端マッピング部２１１が求めたマッピング位置に基づいて、各ペアを整合ペアＡＰ（acccordant pair）、非整合ペアＤＰ（discordant pair）、不完全マップペアＩＰ（incompletely mapped pair）及びマッピング不能ペアＵＰ（unmappable pair）のいずれかに分類する。ペア分類部２１２は、ペアの分類結果をマッピング失敗頻度評価部２２１に出力する。
（ステップ３０７）
　マッピング失敗頻度評価部２２１は、ペアの分類結果に基づいて、標準ゲノム配列Ｇの各塩基位置で整合ペアＡＰの頻度と不完全マップペアＩＰの頻度を求める。
（ステップ３０８）
　マッピング失敗頻度評価部２２１は、実際にマッピングに失敗した頻度を、整合ペアＡＰと不完全マップペアＩＰの頻度の合計に対する不完全マップペアＩＰの相対頻度として計算する。また、マッピング失敗確率評価部２２２は、マッピングに失敗する確率を、処理パラメータ２１８の一つである塩基読み取りエラー率と標準ゲノムＭＬＲＵデータ１２２とを用いて、二項分布の確率モデルに基づいて計算する。
　変異領域候補抽出部２１３は、マッピング失敗頻度評価部２２１から実際にマッピングに失敗した頻度の入力を受け付け、マッピング失敗確率評価部２２２からマッピングに失敗する確率の入力を受け付ける。そして、変異領域候補抽出部２１３は、標準ゲノム配列Ｇ上を走査（スキャン）して、実際にマッピングに失敗した頻度がマッピングに失敗する確率よりも有意に高い領域を変異領域の候補として一つずつ抽出する。
（ステップ３０９）
　変異領域候補抽出部２１３は、変異領域候補が抽出されたか否かを判定する。変異領域候補がなかった場合（Ｎｏ）は、全処理を終了する。このとき、変異領域候補抽出部２１３は、変異領域候補がなかったことを示す結果を表示／出力部２０６に出力してもよく、表示／出力部２０６は、変異領域候補がなかったことを示す結果をＧＵＩ１０８として表示／出力装置１０６（ディスプレイ）に表示させてもよい。変異領域候補が得られた場合（Ｙｅｓ）は、処理はステップ３１０に移行する。
（ステップ３１０）
　変異領域候補抽出部２１３は、座標Ｘを抽出された変異領域候補の右端のＸ座標として、以下の処理を繰り返す。
（ステップ３１１）
　変異領域候補抽出部２１３は、座標Ｘがその変異領域候補内にあるか否かを判定する。Ｘがその変異領域候補内にない場合（Ｎｏ）は、ステップ３０８（変異領域候補の抽出処理）に戻って次の候補を探す。座標Ｘがその変異領域候補内にある場合（Ｙｅｓ）は、処理はステップ３１２に移行する。
（ステップ３１２）
　リード配列抽出部２１４は、座標Ｘを左端とするシード配列Ｓを求め、個人ゲノムリード配列辞書データ１２３から、シード配列Ｓを含むリード配列を全て抽出する。ここで、座標Ｘを左端とするシード配列Ｓとは、標準ゲノム配列データ１１１から、座標Ｘの位置にある要素を始点として長さＭＬＲＵ（＋Ｘ）の部分配列を取り出したものである。シード配列Ｓを含む全てのリード配列を抽出するためには、個人ゲノムリード配列辞書データ１２３に対して公知の方法（ＬＦマッピング、last-to-first column mapping）を適用すればよい。
（ステップ３１３）
　リード配列抽出部２１４は、ステップ３１２において抽出されたリード配列が存在するか否かを判定する。リード配列がない場合（Ｎｏ）は、処理はステップ３１８に移行する。リード配列がある場合（Ｙｅｓ）は、リード配列抽出部２１４は、抽出されたリード配列をアラインメント部２１５に出力し、処理はステップ３１４に移行する。
（ステップ３１８）
　変異領域候補抽出部２１３は、座標Ｘから１を引いて、ステップ３１１に戻る。
（ステップ３１４）
　アラインメント部２１５は、抽出されたリード配列と標準ゲノム配列Ｇとを、それらに共通に含まれるシード配列Ｓの左端Ｘで揃えてアラインメントし、そのアラインメントをＸの左方にまで延長する。アラインメント部２１５は、アラインメントの延長した結果をブレークポイント抽出部２１６に出力する。
（ステップ３１５）
　ブレークポイント抽出部２１６は、座標Ｘの左方に延長されたアラインメントにおいて塩基の不一致が頻出するか否かを判定する。不一致が頻出しない場合（Ｎｏ）には、処理はステップ３１８に移行し、変異領域候補抽出部２１３は、座標Ｘから１を引いて、ステップ３１１に戻る。不一致が頻出する場合（Ｙｅｓ）には、処理はステップ３１６に移行する。
（ステップ３１６）
　ブレークポイント抽出部２１６は、座標Ｘにブレークポイントがあると判定し、判定結果を表示／出力部２０６に出力する。
（ステップ３１７）
　表示／出力部２０６は、ブレークポイントの座標Ｘを表示／出力装置１０６に出力する。表示／出力部２０６は、例えば、ブレークポイントの座標Ｘ及びその周辺のアラインメント結果や、後述するようなブレークポイント周辺における整合ペアＡＰの頻度と不完全マップペアＩＰの頻度の変動をＧＵＩ１０８として、表示／出力装置１０６（ディスプレイ）に表示若しくはファイル出力することができる。
　その後、処理はステップ３０８に戻り、同様の処理が繰り返される。ステップ３０９において新たな変異領域候補がないと判定されたときに、全処理が終了する。
＜ＭＬＲＵの計算方法＞
　図４は、ステップ３０４においてＭＬＲＵ計算部２０２が標準ゲノム配列データ１１１から標準ゲノムＭＬＲＵデータ１２２を計算する方法を示すフローチャートである。
（ステップ４０１）
　ＭＬＲＵ計算部２０２は、記憶装置１０３から標準ゲノム配列データ１１１を読み出す。
（ステップ４０２）
　ＭＬＲＵ計算部２０２は、記憶装置１０３から標準ゲノム配列データ１１１を読み出し、それらの配列とその相補鎖配列（reverse complementary sequences）を区切り文字＄で連結して一本の標準ゲノム配列Ｇとする。
（ステップ４０３）
　ＭＬＲＵ計算部２０２は、公知の方法（非特許文献４）により、標準ゲノム配列ＧのＰＬＣＰ（Permutated Longest Common Prefix Array）を計算する。
（ステップ４０４）
　ＭＬＲＵ計算部２０２は、標準ゲノム配列Ｇ上の各座標Ｘについて、ＭＬＵ（Ｘ）＝ＰＣＬＰ（Ｘ）＋１と定める。ＰＬＣＰの性質により、ＭＬＵ（minimum length for uniqueness）は、座標Ｘを起点とする標準ゲノム配列Ｇの部分配列が標準ゲノム配列Ｇ内で完全一致の条件の下で一意的となる（即ち、その部分配列がＸ以外の箇所を起点とする標準ゲノム配列Ｇの他の部分配列と完全一致することはない）ために必要な最小の長さを与える。
　それに対して、minimum length for robust uniquenessを表すＭＬＲＵ（Ｘ）は、座標Ｘを起点とする標準ゲノム配列Ｇの部分配列が標準ゲノム配列Ｇ内で１塩基以下の違いを許容するロバストな比較条件の下で一意的となる（即ち、その部分配列がＸ以外の箇所を起点とする標準ゲノム配列Ｇの他の部分配列と高々１塩基の違いを除いて一致することはない）ために必要な最小の長さを表すものと定める。ＭＬＲＵとＭＬＵとを比較すると、部分配列を比較する条件を緩和したため、ＭＬＲＵ（Ｘ）の値はＭＬＵ（Ｘ）の値よりも常に大きい。そこで、ＭＬＲＵ計算部２０２は、以下に示す手順で標準ゲノム配列Ｇの各座標ＸでのＭＬＲＵ（Ｘ）の値Ｌを計算する際には、Ｌの初期値としてＬ＝ＭＬＵ（Ｘ）＋１を用いる（ステップ４０８で後述）。
　部分配列の比較により標準ゲノム配列Ｇ内の座標Ｘを定める際、ＭＬＵの長さの部分配列を用いると、個人ゲノムの中に多数含まれているＳＮＰ（Single Nucleotide Polymorphism）の影響を受けて、座標Ｘの位置を誤ることが起きる。それに対して、ＭＬＲＵの長さの部分配列を用いた場合には、配列比較の条件が１塩基以下の違いを許容するように緩和されているため、ＳＮＰの影響を受けずにロバストに座標Ｘを正しく定めることができる。
（ステップ４０５）
　ＭＬＲＵ計算部２０２は、記憶装置１０３から標準ゲノム配列辞書データ１２１を読み出す。
（ステップ４０６）
　ＭＬＲＵ計算部２０２は、座標Ｘに標準ゲノム配列Ｇの左端の座標である１を代入する。
（ステップ４０７）
　ＭＬＲＵ計算部２０２は、座標Ｘを標準ゲノム配列Ｇの右端の座標Ｘｍａｘ（Ｘｍａｘ＝標準ゲノム配列Ｇの長さ）と比較する。Ｘ＞Ｘｍａｘであるならば（Ｙｅｓ）、標準ゲノムＭＬＲＵデータの計算処理を終了する。Ｘ＞Ｘｍａｘでない場合（Ｎｏ）は、ステップ４０８に移行する。
（ステップ４０８）
　ＭＬＲＵ計算部２０２は、標準ゲノム配列Ｇの座標ＸでのＭＬＲＵ（Ｘ）の値Ｌの初期値として、ＭＬＵ（Ｘ）＋１をセットする。
（ステップ４０９）
　ＭＬＲＵ計算部２０２は、標準ゲノム配列Ｇ内で座標Ｘを起点（左端）とする長さＬの部分配列をＳとする。
（ステップ４１０）
　ＭＬＲＵ計算部２０２は、部分配列Ｓに１塩基の置換、挿入又は欠失を入れた配列をＱとする。長さＬの部分配列Ｓに１塩基置換を入れた配列Ｑは３Ｌ通りある。何故ならば、Ｌ箇所ある部分配列Ｓの塩基位置のそれぞれにおいて、元の塩基（Ａ，Ｃ，Ｇ，Ｔのいずれか）をそれ以外の塩基に置き換える方法が３通りあるからである。同様に、長さＬの部分配列Ｓに１塩基挿入を入れた配列Ｑは４Ｌ通りある。また、長さＬの部分配列Ｓに１塩基欠失を入れた配列ＱはＬ通りある。したがって、長さＬの部分配列Ｓに１塩基の置換、挿入又は欠失を入れた配列Ｑは全部で８Ｌ通りある（ただし、部分配列Ｓの中に同じ塩基が連続して並ぶ場合には、８Ｌ通り数え方に重複が含まれる）。
（ステップ４１１）
　ＭＬＲＵ計算部２０２は、一つ一つの配列Ｑに対して、公知の方法（非特許文献２）により、標準ゲノム配列辞書データ１２１を利用して、標準ゲノム配列Ｇ内に配列Ｑが出現する回数Ｏｃｃ（Ｇ，Ｑ）を求める。そして、ＭＬＲＵ計算部２０２は、Ｏｃｃ（Ｇ，Ｑ）＞０であるか否かを判定する。Ｏｃｃ（Ｇ，Ｑ）の値が正である（Ｙｅｓ）ならば、ステップ４１２に移行する。一方、Ｏｃｃ（Ｇ，Ｑ）の値が０である（Ｎｏ）ならば、ステップ４１３に移行する。
（ステップ４１２）
　ＭＬＲＵ計算部２０２は、Ｌの値を１だけ増やして、部分配列Ｓを定めるステップ４０９に戻る。
（ステップ４１３）
　ＭＬＲＵ計算部２０２は、前述した８Ｌ通りの全ての配列Ｑのテストが完了したか否かを判定する。全ての配列Ｑのテストが完了した場合（Ｙｅｓ）には、ステップ４１４に移行する。全ての配列Ｑのテストが完了していなかった場合（Ｎｏ）には、配列Ｑを定めるステップ４１０に戻って、他の配列Ｑに対して、ステップ４１１のテスト処理を行う。
（ステップ４１４）
　ＭＬＲＵ計算部２０２は、座標ＸにおけるＭＬＲＵの値であるＭＬＲＵ（Ｘ）をＬと定めて、標準ゲノムＭＬＲＵデータ１２２として出力し、記憶装置１０３に格納する。
（ステップ４１５）
　ＭＬＲＵ計算部２０２は、Ｘを１だけ増やして、標準ゲノム配列Ｇの右端の座標Ｘｍａｘとの比較処理（ステップ４０７）に戻る。ステップ４０７においてＸ＞Ｘｍａｘとなったとき、処理を終了する。
＜マッピング位置を求める方法＞
　図５は、ステップ３０５においてペア末端マッピング部２１１が各ペアの両端のマッピング位置を求める方法を示すフローチャートである。
（ステップ５０１）
　ペア末端マッピング部２１１は、記憶装置１０３から標準ゲノム配列辞書データ１２１及び標準ゲノムＭＬＲＵデータ１２２を読み出す。
（ステップ５０２）
　ペア末端マッピング部２１１は、記憶装置１０３から個人ゲノムリード配列データ１１２を読み出す。
（ステップ５０３）
　個人ゲノムリード配列データ１１２はペアをなすリード配列の集合であり、各ペアは２本のリード配列からなる。ペアの両端とは、ペアをなすリード配列それぞれの５’末端を意味する。ペア末端マッピング部２１１は、未処理のリード配列があるか否かを判定する。全てのリード配列に対する処理を行い、未処理のリード配列がなくなったとき（Ｎｏ）、処理は終了する。未処理のリード配列がある場合（Ｙｅｓ）、ステップ５０４に移行する。
（ステップ５０４）
　ペア末端マッピング部２１１は、未処理のリード配列をＲとし、次の処理対象とする。
（ステップ５０５）
　ペア末端マッピング部２１１は、クエリー長Ｌを１に初期化する。
（ステップ５０６）
　ペア末端マッピング部２１１は、リード配列Ｒの５’末端から長さＬの部分配列をとり、それをクエリー配列Ｑとする。
（ステップ５０７）
　ペア末端マッピング部２１１は、公知の方法（非特許文献２、非特許文献３）により、標準ゲノム配列辞書データ１２１を利用して、標準ゲノム配列Ｇ内にクエリー配列Ｑが出現する回数Ｏｃｃ（Ｇ，Ｑ）を求め、その値をＦとする。
（ステップ５０８）
　ペア末端マッピング部２１１は、Ｆ＞１であるか否かを判定する。Ｆ＞１ならば（Ｙｅｓ）、ステップ５０９に移行する。Ｆ＞１でないならば（Ｎｏ）、ステップ５１０に移行する。
（ステップ５０９）
　ペア末端マッピング部２１１は、Ｌを１だけ増やして、クエリー配列Ｑを定めるステップ５０６に戻る。
（ステップ５１０）
　ペア末端マッピング部２１１は、Ｆ＝０であるか否かを判定する。Ｆ＝０ならば（Ｙｅｓ）、ステップ５１１に移行する。Ｆ＝０でない場合（Ｎｏ）、即ち、Ｆ＝１の場合は、ステップ５１２に移行する。
（ステップ５１１）
　ペア末端マッピング部２１１は、リード配列Ｒのマッピングは失敗したと判定してその判定結果を変異領域候補抽出部２１３に出力する。その後、処理はステップ５０３に戻る。
（ステップ５１２）
　ペア末端マッピング部２１１は、クエリー配列Ｑの唯一の出現位置の座標Ｘを公知の方法（非特許文献２）により求めて、Ｌの値をＭＬＲＵ（Ｘ）に更新する。
（ステップ５１３）
　ペア末端マッピング部２１１は、リード配列Ｒの５’末端から長さＬの部分配列をとり、それをクエリー配列Ｑとする。
（ステップ５１４）
　ペア末端マッピング部２１１は、公知の方法（非特許文献２、非特許文献３）により、標準ゲノム配列辞書データ１２１を利用して、標準ゲノム配列Ｇ内にクエリー配列Ｑが出現する回数Ｏｃｃ（Ｇ，Ｑ）を求め、その値をＦとする。
（ステップ５１５）
　ペア末端マッピング部２１１は、Ｆを１と比較し、Ｆ＝１であるか否かを判定する。Ｆ＝１の場合（Ｙｅｓ）は、ステップ５１６に移行する。Ｆ＝１でない場合（Ｎｏ）は、Ｆ＝０であり、ステップ５１１に移行する。その後、ステップ５０３に戻り、未処理のリード配列がないと判定されたとき、処理は終了する。
（ステップ５１６）
　ペア末端マッピング部２１１は、リード配列Ｒのマッピングが成功したと判定して、その判定結果とリード配列Ｒのマッピング座標Ｘ（Ｒ）がＸであることを変異領域候補抽出部２１３に出力する。その後、ステップ５０３に戻り、未処理のリード配列がないと判定されたとき、処理は終了する。
＜ペアの分類方法＞
　図６Ａ～６Ｆは、ステップ３０６においてペア分類部２１２が各ペアを整合ペアＡＰ、非整合ペアＤＰ、不完全マップペアＩＰ又はマッピング不能ペアＵＰのいずれかに分類する方法を説明するための図である。この分類は、ペアごとに、それぞれ独立に行う。図６Ａ～Ｆの各々は、ペアがＡＰ、ＤＰ、ＩＰ又はＵＰのいずれかに分類される代表的な場合を示している。
　図６Ａ～Ｆの各々において、横軸６０１は、標準ゲノム配列Ｇ上の位置（座標）を示す座標軸である。ペアをなす２本のリード配列をそれぞれＲ１及びＲ２として、ペア末端マッピング部２１１によりリード配列Ｒ１の末端のマッピングが成功した場合、そのマッピング座標Ｘ（Ｒ１）を略してＸ１と表す。同様に、リード配列Ｒ２の末端のマッピングが成功した場合、そのマッピング座標Ｘ（Ｒ２）を略してＸ２と表す。
　右向きの矢印で示されたリード配列Ｒ１は、リード配列の末端が標準ゲノム配列Ｇのプラス鎖にマッピングされることを意味し、左向きの矢印で示されたリード配列Ｒ１は、リード配列の末端が標準ゲノム配列Ｇのマイナス鎖にマッピングされることを意味する。右向き又は左向きの矢印で示されたリード配列Ｒ２についても同様である。
　×印６０６（図６Ｅ、図６Ｆ）は、ペア末端マッピング部２１１においてリード配列Ｒ２の末端のマッピングが失敗してマッピング位置が定まらなかった場合を表す。同様に、×印６０７（図６Ｆ）は、リード配列Ｒ１の末端のマッピングが失敗してマッピング位置が定まらなかった場合を表す。
　また、インサート長の平均値（又は代表値）をＭとし、インサート長の平均値からのずれの大きさが正常の範囲にあると許容できる閾値をＶとする。インサート長の平均値Ｍ及び閾値Ｖは、処理パラメータ２１８の一部であり、入力部２０５がＧＵＩ１０８を通じて取得し、メモリ１０２又は記憶装置１０３に記憶するものである。
　図６Ａは、整合ペアＡＰに分類される場合を示している。図６Ａに示すように、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも成功し、リード配列Ｒ１及びＲ２の末端のマッピング先は同じ染色体にあり、リード配列Ｒ１はプラス鎖にマッピングされ、リード配列Ｒ２はマイナス鎖にマッピングされ、座標Ｘ１は座標Ｘ２の左方に位置し、座標Ｘ１とＸ２の間の距離Ｄ＝Ｘ２－Ｘ１＋１とインサート長の平均値Ｍとの差が閾値Ｖの範囲内に収まっている（｜Ｄ－Ｍ｜≦Ｖを満たす）場合、このペアは整合ペアＡＰに分類される。
　また、図６Ａにおけるリード配列Ｒ１及びＲ２の役割を交換した場合も整合ペアＡＰに分類することができる。即ち、リード配列Ｒ２及びＲ１の末端のマッピングが両者とも成功し、両者のマッピング先は同じ染色体にあり、リード配列Ｒ２はプラス鎖にマッピングされ、リード配列Ｒ１はマイナス鎖にマッピングされ、座標Ｘ２は座標Ｘ１の左方に位置し、座標Ｘ２とＸ１との間の距離Ｄ＝Ｘ１－Ｘ２＋１とインサート長の平均値Ｍとの差が閾値Ｖの範囲内に収まっている場合、このペアは整合ペアＡＰに分類される。
　一方、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも成功しても、整合ペアＡＰに分類できなかった場合は全て、非整合ペアＤＰに分類する。例えば、リード配列Ｒ１とＲ２の末端のマッピングが両者とも成功し、両者のマッピング先が別の染色体にある場合、このペアは非整合ペアＤＰに分類される。
　図６Ｂは、非整合ペアＤＰに分類される場合を示している。図６Ｂに示すように、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも成功し、両者のマッピング先は同じ染色体にあり、両者ともプラス鎖にマッピングされる場合、このペアを非整合ペアＤＰに分類する。同様に、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも成功し、両者のマッピング先は同じ染色体にあり、両者ともマイナス鎖にマッピングされる場合も、このペアは非整合ペアＤＰに分類される。
　図６Ｃは、非整合ペアＤＰに分類される場合を示している。図６Ｃに示すように、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも成功し、両者のマッピング先は同じ染色体にあり、リード配列Ｒ１はプラス鎖にマッピングされ、リード配列Ｒ２はマイナス鎖にマッピングされ、座標Ｘ１が座標Ｘ２の右方に位置する場合、このペアは非整合ペアＤＰに分類される。また、図６Ｃでリード配列Ｒ１及びＲ２の役割を交換した場合も同様に、非整合ペアＤＰに分類される。
　図６Ｄは、非整合ペアＤＰに分類される場合を示している。図６Ｄに示すように、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも成功し、両者のマッピング先は同じ染色体にあり、リード配列Ｒ１はプラス鎖にマッピングされ、リード配列Ｒ２はマイナス鎖にマッピングされ、座標Ｘ１が座標Ｘ２の左方に位置し、座標Ｘ１とＸ２との間の距離Ｄ＝Ｘ２－Ｘ１＋１とインサート長の平均値Ｍとの差が閾値Ｖを超えた（｜Ｄ－Ｍ｜＞Ｖを満たす）場合、このペアは非整合ペアＤＰに分類される。また、図６Ｄでリード配列Ｒ１及びＲ２の役割を交換した場合も同様に、非整合ペアＤＰに分類される。
　図６Ｅは、不完全マップペアＩＰに分類される場合を示している。図６Ｅに示すように、リード配列Ｒ１の末端のマッピングに成功し、リード配列Ｒ２の末端のマッピングに失敗した場合、このペアは不完全マップペアＩＰに分類される。座標Ｘ１を不完全マップペアＩＰのマッピング位置と呼ぶ。また、図６Ｅでリード配列Ｒ１及びＲ２の役割を交換した場合も同様である。
　図６Ｆは、マッピング不能ペアＵＰに分類される場合を示している。図６Ｆに示すように、リード配列Ｒ１及びＲ２の末端のマッピングが両者とも失敗した場合、このペアはマッピング不能ペアＵＰに分類される。
＜リード配列の頻度の計算方法＞
　図７Ａ及び７Ｂは、標準ゲノム配列Ｇの各座標Ｘでリード配列の頻度（リード頻度）を計算する方法を説明する図である。
　図７Ａは、標準ゲノム配列Ｇのプラス鎖上でリード配列の頻度を求める場合を示している。プラス鎖上では、リード配列Ｒの５’末端が座標ＸＲの位置にマッピングされるとき、３’末端はＸＲ＋Ｌ－１にマッピングされる。ここで、Ｌはリード配列の長さを表す。標準ゲノム配列Ｇの座標Ｘでリード頻度は、ＸＲ≦Ｘ≦ＸＲ＋Ｌ－１を満たすリード配列の本数Ｐ（Ｘ）によって与えられる。Ｐ（Ｘ）を効率的に計算するためには、通常よく行われるように、全リード配列Ｒを５’末端のマッピング位置の座標ＸＲの昇順にソートして、全ての座標Ｘについて昇順にＸＲ≦Ｘ≦ＸＲ＋Ｌ－１を満たすリード配列の本数を数えればよい。
　図７Ｂは、標準ゲノム配列Ｇのマイナス鎖上でリード配列の頻度を求める場合を示している。マイナス鎖上では、リード配列Ｒの５’末端が座標ＸＲの位置にマッピングされるとき、３’末端はＸＲ－Ｌ＋１にマッピングされる。ここで、Ｌはリード配列の長さを表す。標準ゲノム配列Ｇの座標Ｘでリード頻度は、ＸＲ－Ｌ＋１≦Ｘ≦ＸＲを満たすリード配列の本数Ｍ（Ｘ）によって与えられる。Ｍ（Ｘ）を効率的に計算するためには、通常よく行われるように、全リード配列Ｒを５’末端のマッピング位置の座標ＸＲの降順にソートして、全ての座標Ｘについて降順にＸＲ－Ｌ＋１≦Ｘ≦ＸＲを満たすリード配列の本数を数えればよい。
　標準ゲノム配列Ｇの各塩基位置での整合ペアＡＰの頻度とは、整合ペアＡＰに分類されたペアに属するリード配列Ｒに限定したリード頻度である。同様に、標準ゲノム配列Ｇの各塩基位置での不完全マップペアＩＰの頻度とは、不完全マップペアＩＰに分類されたペアに属するリード配列Ｒに限定したリード頻度である。変異領域候補抽出部２１３は、ペア分類部２１２の分類結果に基づいて、整合ペアＡＰに分類されたペアに属するリード配列と不完全マップペアＩＰに分類されたペアに属するリード配列を求め、図７Ａ及び７Ｂを用いて説明したリード頻度の計算方法を用いて、標準ゲノム配列Ｇの各塩基位置での整合ペアＡＰの頻度と不完全マップペアＩＰの頻度を求める。
＜変異領域候補の抽出方法＞
　図８Ａ及び８Ｂは、ステップ３０８において変異領域候補抽出部２１３により変異が含まれる領域の候補を抽出する方法を説明する図である。
　図８Ａは、ホモ型（homozygous）の変異が含まれる可能性がある領域の候補が選び出される状況を示している。変異領域候補抽出部２１３は、横軸に標準ゲノム配列Ｇの座標Ｘをとり、縦軸にリード頻度をとり、不完全マップペアＩＰの頻度８０３、整合ペアＡＰの頻度８０４及び非整合ペアＤＰの頻度８０５をプロットする。
　座標Ｘにおける不完全マップペアＩＰの頻度、整合ペアＡＰの頻度及び非整合ペアＤＰの頻度をそれぞれＩＰ（Ｘ）、ＡＰ（Ｘ）及びＤＰ（Ｘ）と表す。不完全マップペアＩＰの相対頻度を
　ＲＩＰ（Ｘ）＝ＩＰ（Ｘ）／（ＡＰ（Ｘ）＋ＩＰ（Ｘ））
と定めて、これを「マッピングに失敗した頻度」ともよぶ。また、非整合ペアＤＰの相対頻度を
　ＲＤＰ（Ｘ）＝ＤＰ（Ｘ）／（ＡＰ（Ｘ）＋ＤＰ（Ｘ））
と定める。
　変異領域候補抽出部２１３は、Ｘ軸上を走査して、不完全マップペアＩＰの相対頻度及び非整合ペアＤＰの相対頻度のいずれかが、後述の確率モデルに基づいて有意に大きいと判定される位置が指定された長さ（リード配列長Ｌ程度）以上続く領域を抽出し、抽出した領域を変異が含まれる領域の候補とする。図８Ａ中の領域８１３は、不完全マップペアＩＰの相対頻度が閾値（例えば０．２５程度）を超えて抽出された変異領域の候補である。ホモ型の変異の場合、ＲＩＰ（Ｘ）は最大で１に近い値をとる。一方、領域８１４及び８１５では、不完全マップペアＩＰの相対頻度は閾値以下であり、領域８１４及び８１５の不完全マップペアＩＰは偶発的原因で生じたノイズであると考えられる。また、図８Ａには、非整合ペアＤＰの相対頻度が閾値を超えて抽出された変異領域の候補の例は示していないが、そのような領域では公知の方法（非特許文献１）により非整合ペアを利用して変異を検出できる。
　図８Ｂは、ヘテロ型（heterozygous）の変異が含まれる可能性がある領域の候補が選び出される状況を示している。横軸に標準ゲノム配列Ｇの座標Ｘをとり、縦軸にリード頻度をとり、不完全マップペアＩＰの頻度８２３、整合ペアＡＰの頻度８２４及び非整合ペアＤＰの頻度８２５をプロットする。
　変異領域候補抽出部２１３は、Ｘ軸上を走査して、不完全マップペアＩＰの相対頻度及び非整合ペアＤＰの相対頻度のいずれかが、後述の確率モデルに基づいて有意に大きいと判定される位置が指定された長さ（リード配列長Ｌ程度）以上続く領域を抽出し、抽出した領域を変異が含まれる領域の候補とする。図８Ｂ中の領域８３３は、不完全マップペアＩＰの相対頻度が閾値（例えば０．２５程度）を超えて抽出された変異領域の候補である。ヘテロ型の変異の場合、ＲＩＰ（Ｘ）は最大で０．５に近い値をとる。一方、領域８３４及び８３５では、不完全マップペアＩＰの相対頻度は閾値以下であり、領域８３４及び８３５の不完全マップペアＩＰは偶発的原因で生じたノイズであると考えられる。また、図８Ｂには、非整合ペアＤＰの相対頻度が閾値を超えて抽出された変異領域の候補の例は示していないが、そのような領域では公知の方法（非特許文献１）により非整合ペアを利用して変異を検出できる。
　以上のようなホモ型又はヘテロ型の変異が含まれ得る領域の候補を抽出する方法は、雌雄の区別のある生物の個体ゲノムの解析に適用することができる。
　なお、ペアの頻度のプロットは２つの目的で行う。第一の目的では、上述のように、コンピュータ１００の内部（変異領域候補抽出部２１３）において、標準ゲノム配列Ｇの全体にわたるプロットを行い、Ｘ軸上を走査して不完全マップペアＩＰの相対頻度及び非整合ペアＤＰの相対頻度を閾値と比較して変異が含まれる領域の候補を抽出するために利用する。第２の目的では、抽出された変異領域の候補の周辺のプロットをＧＵＩ１０８として表示／出力装置１０６に表示して、候補が抽出された状況をユーザに示すために利用する。
＜ブレークポイントの抽出方法＞
　図９は、ステップ３１４、３１５及び３１６においてアラインメント部２１５及びブレークポイント抽出部２１６により変異領域候補の内部から構造変異に伴うブレークポイント（ＢＰ）を抽出する方法を説明する図である。図９では、プラス鎖について調べる方法を説明する。マイナス鎖について調べる場合は、左右を反転して同様の方法を用いる。プラス鎖上では、変異領域候補の内部を左から右に順に１塩基ずつ走査して以下の処理を繰り返す。
（ステップ９０１）
　アラインメント部２１５は、走査中の座標をＸとする。
（ステップ９０２）
　アラインメント部２１５は、標準ゲノムＭＬＲＵデータ１２２から、プラス鎖上の座標ＸにおけるＭＬＲＵの値ＭＬＲＵ（＋Ｘ）を取得する。
（ステップ９０３）
　アラインメント部２１５は、標準ゲノム配列データ１１１から、プラス鎖上で左端座標がＸであり長さがＭＬＲＵ（Ｘ）の部分配列を取り出し、それをシード配列Ｓとする。
（ステップ９０４）
　アラインメント部２１５は、走査中の座標Ｘの左に隣接する座標をＸ－１とする。
（ステップ９０５）
　アラインメント部２１５は、標準ゲノムＭＬＲＵデータ１２２から、マイナス鎖上のＸ－１におけるＭＬＲＵの値ＭＬＲＵ（－（Ｘ－１））を取得する。
（ステップ９０６）
　アラインメント部２１５は、標準ゲノム配列データ１１１から、マイナス鎖上で右端座標がＸ－１であり長さがＭＬＲＵ（－（Ｘ－１））の部分配列を取り出し、それを標準延長配列Ｆ０とする。
（ステップ９０７）
　アラインメント部２１５は、公知の方法（非特許文献２、非特許文献３）により、個人ゲノムリード配列辞書データ１２３の中からシード配列Ｓを含む全てのリード配列を求める。
（ステップ９０８）
　アラインメント部２１５は、求めたリード配列内で、シード配列Ｓの左方に隣接する長さＭＬＲＵ（－（Ｘ－１））以上の部分配列を全て求め、これらを個人延長配列とする（図９の例ではＦ１、Ｆ２、Ｆ３、Ｆ４及びＦ５）。アラインメント部２１５は、動的計画法などの公知の方法により、標準延長配列Ｆ０と個人延長配列Ｆ１～Ｆ５のアラインメントを行い、アラインメントの結果をブレークポイント抽出部２１６に出力する。
　ブレークポイント抽出部２１６は、個人延長配列Ｆ１～Ｆ５のそれぞれと標準延長配列Ｆ０との編集距離（edit distance）を求める。ブレークポイント抽出部２１６は、編集距離と標準延長配列Ｆ０の長さＭＬＲＵ（－（Ｘ－１））との比が一定値（例えば０．５程度）以上であるリード配列は構造変異に由来すると判定する。
　また、ブレークポイント抽出部２１６は、個人延長配列Ｆ１～Ｆ５の中で構造変異に由来すると判定された割合を求め、これを変異率とする。図９の例では、Ｆ１～Ｆ５のうち、Ｆ１、Ｆ２及びＦ３を含むリード配列９２１は構造変異に由来すると判定され、Ｆ４及びＦ５を含むリード配列９２２は構造変異に由来しないと判定され、変異率は０．６である。
　ブレークポイント抽出部２１６は、変異率を１（ホモ型の変異の場合）又は０．５（ヘテロ型の変異の場合）と比較して、ホモ型又はヘテロ型の変異の有無を判定する。ホモ型又はヘテロ型の変異ありと判定された場合には、ブレークポイント抽出部２１６は、構造変異に由来すると判定されたリード配列に含まれる個人延長配列（Ｆ１、Ｆ２及びＦ３）が個人ゲノムの構造変異で生じた挿入配列の一部であると判定する。また、ブレークポイント抽出部２１６は、その構造変異のブレークポイントが標準ゲノム配列Ｇ上のシード配列Ｓと標準延長配列Ｆ０の境界の位置にあると判定する。
　ブレークポイント抽出部２１６は、これらの判定結果を構造変異検出結果１１３として出力し、記憶装置１０３に格納するとともに、表示／出力部２０６に出力する。表示／出力部２０６は、構造変異検出結果１１３と、判定に利用したシード配列Ｓ及び延長配列のアラインメント結果をＧＵＩ１０８として表示／出力装置１０６に表示させることで、ユーザに提示する。
＜マッピングに失敗する確率と実際に失敗した頻度について＞
　マッピング失敗確率評価部２２２が確率モデルに基づいてマッピングに失敗する確率を計算する方法、及び、変異領域候補抽出部２１３がマッピングに失敗する確率と実際にマッピングに失敗した頻度とを比較して、実際にマッピングに失敗した頻度が有意に大きいか否かを判定する方法について説明する。
　上述のように、ペア末端マッピング部２１１においても、リード配列抽出部２１４においても、リード配列内のＭＬＲＵに等しい長さの部分配列が用いられる。ヒトゲノムのＭＬＲＵの典型的な長さは２０～３０程度である。ペア末端マッピング部２１１においては、マッピングが成功するための条件は、図５を用いて説明したように、部分配列全体が標準ゲノム配列内の部分配列と完全一致することである。したがって、リード配列内にＤＮＡシーケンシングの際の塩基読み取りエラーが含まれていれば、マッピングに失敗する。部分配列内に生じる塩基読み取りエラーの数は二項分布で近似でき、通常、１塩基当たりのエラー率は０．０１（１％）程度である。したがって、ＭＬＲＵをｍとし、１塩基当たりのエラー率をｅとすれば、変異が存在しない場合に１本のリード配列がマッピングに失敗する確率ｐは、以下の［式１］で与えられる。
Figure JPOXMLDOC01-appb-M000001
　したがって、変異が存在しない場合に、ｎ本のリード配列中ｒ本以上がマッピングに失敗する確率Ｅ（ｒ／ｎ）は、以下の［式２］で与えられる。
Figure JPOXMLDOC01-appb-M000002
　マッピング失敗頻度評価部２２１により、標準ゲノム配列の座標Ｘの位置で整合ペアＡＰの頻度ＡＰ（Ｘ）と不完全マップペアＩＰの頻度ＩＰ（Ｘ）が得られたとき、ｎ＝ＡＰ（Ｘ）＋ＩＰ（Ｘ）とし、ｒ＝ＩＰ（Ｘ）とすると、実際にマッピングに失敗した頻度ＲＩＰ（Ｘ）はｒ／ｎに等しい。このとき、［式２］で計算されたマッピングに失敗する確率Ｅ（ｒ／ｎ）は、実際には変異がないにもかかわらず、ｎ本中ｒ本以上のリード配列のマッピングが失敗する確率を表す。そこで、変異領域候補抽出部２１３は、αを誤り率（ＧＵＩ１０８を通じてユーザにより指定される処理パラメータ２１８の一つであり、通常１～５％程度）として、マッピングに失敗する確率Ｅ（ｒ／ｎ）が誤り率α未満になるとき、マッピングに失敗した頻度がマッピングに失敗する確率（確率モデルに基づき計算される）と比較して有意に大きいと判定する。これにより、実際には変異が存在しない場合に、マッピングに失敗した頻度がマッピングに失敗する確率と比較して有意に大きいと判定されることにより誤った変異領域の候補が抽出される確率は、α以下となる。
＜第１の実施形態のまとめ＞
　以上のように、第１の実施形態における個人ゲノム（個体ゲノム）の構造変異検出方法は、コンピュータ１００のＣＰＵ１０１（プロセッサ）により実行され、入力部２０５により、標準ゲノム配列データ１１１を受信することと、マッピング失敗確率評価部２２２により、標準ゲノム配列データ１１１（標準ゲノム配列Ｇ）の各位置（座標Ｘ）において、変異がないと仮定した場合にマッピングに失敗する確率を評価することと、入力部２０５により、個人ゲノムのリード配列データ１１２を受信することと、ペア末端マッピング部２１１により、リード配列データ１１２のそれぞれの末端を標準ゲノム配列データ１１１上にマッピングすることと、マッピング失敗頻度評価部２２１により、標準ゲノム配列データ１１１の各位置でリード配列データ１１２の末端のマッピングに失敗した頻度を評価することと、変異領域候補抽出部２１３により、標準ゲノム配列データ１１１の各位置において、マッピングに失敗した頻度がマッピングに失敗する確率と比較して有意に大きいか否かを判定することと、変異領域候補抽出部２１３により、標準ゲノム配列データ１１１上で、マッピングに失敗した頻度が有意に大きいと判定される位置が所定の長さ以上続く領域を、変異領域候補として抽出することと、表示／出力部２０６により、変異領域候補を構造変異の検出結果として出力することと、を含む。
　このように、変異がないと仮定してマッピングに失敗する確率を評価し、従来は用いられていない情報であるマッピングに失敗した頻度を調べ、上記確率に合わないほど多数の失敗があった場合に、何らかの異常が生じたことを示す証拠だと考えて、その領域に構造変異が生じていると判定する。これにより、インサート長を超える長さの挿入変異がある可能性がある領域を抽出することができる。
　また、第１の実施形態に係る個人ゲノムの構造変異検出方法は、リード配列抽出部２１４により、変異領域候補から標準ゲノム配列データの部分配列であるシード配列Ｓを取り出し、シード配列Ｓを含むリード配列データを抽出することと、アラインメント部２１５により、標準ゲノム配列データとリード配列データとを共通のシード配列Ｓを揃えて隣接する延長配列までアラインメントすることと、ブレークポイント抽出部２１６により、延長配列において、リード配列データ及び標準ゲノム配列データの編集距離と延長配列の長さとの比が所定の値以上である場合、そのリード配列データは構造変異に由来すると判定することと、ブレークポイント抽出部２１６により、構造変異に由来すると判定されたリード配列データの割合（変異率）を、ホモ型の変異の場合は１と比較し、ヘテロ型の変異の場合は０．５と比較してホモ型又はヘテロ型の変異の有無を判定することと、含む。
　これにより、挿入配列長がインサート長より小さいか否かに関わらず、挿入変異に伴うブレークポイントを検出でき、また、挿入配列の一部を取得することができる。構造変異などの異常が生じていなければ、個人ゲノムの配列データを用いて計算したペアの片方でマッピングに失敗した相対頻度は、標準ゲノム配列のみから計算したマッピングに失敗する確率にほぼ等しいはずであるが、これらの不一致を検出することにより、変異候補領域を絞り込める。それにより、アラインメントを用いてブレークポイント検出処理を行うゲノム領域を限定でき、計算コストを抑える効果が得られる。
［第２の実施形態］
　第１の実施形態では、ショートリード型のＤＮＡシーケンサで得られた個人ゲノムの配列データを用いて個人ゲノムの構造変異検出を行い、ブレークポイントの位置を塩基単位で詳細に決定する方法を示した。本開示の第２の実施形態では、ショートリード型のＤＮＡシーケンサで得られた配列データを用いて変異領域候補を求め、その結果を出力するにとどめる方法を提案する。
＜構造変異検出システム＞
　図１０は、第２の実施形態に係る構造変異検出システム２の機能ブロック図である。本実施形態の構造変異検出システム２は、第１の実施形態の構造変異検出システム１の構成（図２）から必要なものを抜き出したものであり、第１の実施形態と同じ構成には同様の符号を付しており、その説明を省略する。標準ゲノム配列辞書作成部２０１、ＭＬＲＵ計算部２０２、入力部２０５、表示／出力部２０６、ペア末端マッピング部２１１、ペア分類部２１２、変異領域候補抽出部１０１３、マッピング失敗頻度評価部２２１及びマッピング失敗確率評価部２２２を有する。変異領域候補抽出部１０１３は、変異領域の候補を一度に全て抽出する点で、第１の実施形態の変異領域候補抽出部２１３（変異領域の候補を一つずつ抽出する）と異なっている。
＜構造変異検出方法＞
　図１１は、第２の実施形態に係る構造変異検出システム２の処理全体を示すフローチャートである。第１の実施形態と同様の処理には同様の符号を付しており、その説明を省略する。
（ステップ３０１及び３０２）
　入力部２０５は、第１の実施形態と同様にステップ３０１及び３０２を実行する。
（ステップ１１０３）
　標準ゲノム配列辞書作成部２０１は、公知の方法（非特許文献５）を利用して標準ゲノム配列ＧのＢｕｒｒｏｗｓ－Ｗｈｅｅｌｅｒ変換（ＦＭインデクス）を計算し、その結果を標準ゲノム配列辞書データ１２１とする。標準ゲノム配列辞書作成部２０１は、生成した標準ゲノム配列辞書データ１２１を記憶装置１０３に格納する。
（ステップ３０４～３０７）
　第１の実施形態と同様にステップ３０４～３０７が実行される。
（ステップ１１０８）
　マッピング失敗頻度評価部２２１は、実際にマッピングに失敗した頻度を、整合ペアＡＰと不完全マップペアＩＰの頻度の合計に対する不完全マップペアＩＰの相対頻度として計算する。また、マッピング失敗確率評価部２２２は、マッピングに失敗する確率を、処理パラメータ２１８の一つである塩基読み取りエラー率と標準ゲノムＭＬＲＵデータ１２２とを用いて、二項分布の確率モデルに基づいて計算する。
　変異領域候補抽出部１０１３は、マッピング失敗頻度評価部２２１から実際にマッピングに失敗した頻度の入力を受け付け、マッピング失敗確率評価部２２２からマッピングに失敗する確率の入力を受け付ける。そして、変異領域候補抽出部１０１３は、標準ゲノム配列Ｇ上を走査（スキャン）して、実際にマッピングに失敗した頻度がマッピングに失敗する確率よりも有意に高い領域を変異領域の候補として全て抽出する。
　変異領域候補抽出部１０１３は、抽出した変異領域の候補を構造変異検出結果１１３として出力し、処理を終了する。変異領域候補抽出部１０１３は、構造変異検出結果１１３を表示／出力部２０６に出力して、ＧＵＩ１０８として表示してもよい。
　本実施形態では、変異領域候補の出力結果は次のように利用することができる。公知の方法（ＰＣＲ法やハイブリダイゼーション法）を利用して、個人ゲノムのＤＮＡサンプルの中から変異領域候補を含むＤＮＡ断片を増幅又は濃縮したサンプルを調整し、ロングリード型のＤＮＡシーケンサ（スループットは低いが数キロ塩基から数十キロ塩基程度の長い配列の塩基を直接読み取ることができるシーケンサ）を用いて、構造変異配列を直接読み取ることにより、構造変異解析を行う。
［変形例］
　本開示は、上述した実施形態に限定されるものでなく、様々な変形例を含んでいる。例えば、上述した実施形態は、本開示を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備える必要はない。また、ある実施形態の一部を他の実施形態の構成に置き換えることができる。また、ある実施形態の構成に他の実施形態の構成を加えることもできる。また、各実施形態の構成の一部について、他の実施形態の構成の一部を追加、削除又は置換することもできる。
１００…コンピュータ
１０１…ＣＰＵ
１０２…メモリ
１０３…記憶装置
１０４…ネットワークインタフェース
１０５…入力装置
１０６…表示／出力装置

Claims (14)

  1. 　コンピュータのプロセッサにより実行される個体ゲノムの構造変異検出方法であって、
    　前記プロセッサにより、標準ゲノム配列データを受信することと、
    　前記プロセッサにより、前記標準ゲノム配列データの各位置において、変異がないと仮定した場合にマッピングに失敗する確率を計算することと、
    　前記プロセッサにより、個体ゲノムのリード配列データを受信することと、
    　前記プロセッサにより、前記リード配列データのそれぞれの末端を前記標準ゲノム配列データ上にマッピングすることと、
    　前記プロセッサにより、前記標準ゲノム配列データの各位置で前記リード配列データの末端のマッピングに失敗した頻度を計算することと、
    　前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記マッピングに失敗した頻度が前記マッピングに失敗する確率と比較して有意に大きいか否かを判定することと、
    　前記プロセッサにより、前記標準ゲノム配列データ上で、前記マッピングに失敗した頻度が有意に大きいと判定される位置が所定の長さ以上続く領域を、変異領域候補として抽出することと、
    　前記プロセッサにより、前記変異領域候補を構造変異の検出結果として出力することと、を含む個体ゲノムの構造変異検出方法。
  2. 　前記プロセッサにより、前記変異領域候補から前記標準ゲノム配列データの部分配列であるシード配列を取り出し、前記シード配列を含む前記リード配列データを抽出することと、
    　前記プロセッサにより、前記標準ゲノム配列データと前記リード配列データとを共通の前記シード配列を揃えて隣接する延長配列までアラインメントすることと、
    　前記プロセッサにより、前記延長配列において、前記リード配列データ及び前記標準ゲノム配列データの編集距離と前記延長配列の長さとの比が所定の値以上である場合、そのリード配列データは構造変異に由来すると判定することと、
    　前記プロセッサにより、前記構造変異に由来すると判定された前記リード配列データの割合を、ホモ型の変異の場合は１と比較し、ヘテロ型の変異の場合は０．５と比較してホモ型又はヘテロ型の変異の有無を判定することと、をさらに含む請求項１に記載の個体ゲノムの構造変異検出方法。
  3. 　前記プロセッサは、前記リード配列データのＢｕｒｒｏｗｓ－Ｗｈｅｅｌｅｒ変換を利用して、前記シード配列を含む前記リード配列データを抽出する請求項２に記載の個体ゲノムの構造変異検出方法。
  4. 　前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記各位置を始点とする部分配列が前記標準ゲノム配列データ内で一意的になる最小の長さを求めることをさらに含み、
    　前記プロセッサは、前記マッピングに失敗する確率を、前記最小の長さの中に塩基読み取りエラーが１個以上生じる確率として計算する請求項１に記載の個体ゲノムの構造変異検出方法。
  5. 　前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記各位置を始点とする部分配列が前記標準ゲノム配列データ内で１個以下のミスマッチを許す条件下において一意的になる最小の長さを求めることをさらに含み、
    　前記プロセッサは、前記マッピングに失敗する確率を、前記最小の長さの中に塩基読み取りエラーが１個以上生じる確率として計算する請求項１に記載の個体ゲノムの構造変異検出方法。
  6. 　前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記各位置を始点とする部分配列が前記標準ゲノム配列データ内で一意的になる最小の長さ、又は、１個以下のミスマッチを許す条件下において一意的になる最小の長さを求めることをさらに含み、
    　前記プロセッサは、前記リード配列データの末端が前記マッピングに成功するための条件を、
    　前記リード配列データの末端から切り出した部分配列が、前記最小の長さの前記標準ゲノム配列データの部分配列と完全一致することとする請求項２に記載の個体ゲノムの構造変異検出方法。
  7. 　前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記各位置を始点とする部分配列が前記標準ゲノム配列データ内で一意的になる最小の長さ、又は、１個以下のミスマッチを許す条件下において一意的になる最小の長さを求めることをさらに含み、
    　前記プロセッサは、前記シード配列と前記延長配列の長さを前記最小の長さとする請求項２に記載の個体ゲノムの構造変異検出方法。
  8. 　前記ホモ型又はヘテロ型の変異があると判定された場合に、
    　前記プロセッサにより、前記構造変異に由来すると判定された前記リード配列データに含まれる前記延長配列が前記個体ゲノムの前記構造変異で生じた挿入配列の一部であると判定することと、
    　前記プロセッサにより、前記構造変異のブレークポイントが前記標準ゲノム配列データ上で前記シード配列と前記延長配列との境界の位置にあると判定することと、
    　前記プロセッサにより、前記挿入配列の一部であるとの判定結果及び前記ブレークポイントの位置を出力することと、をさらに含む請求項２に記載の個体ゲノムの構造変異検出方法。
  9. 　前記プロセッサにより、表示画面として、前記変異領域候補の周辺で、前記標準ゲノム配列データの座標軸上に沿って、前記リード配列データの末端の前記マッピングに成功した頻度と失敗した頻度とをプロット表示することをさらに含む請求項１に記載の個体ゲノムの構造変異検出方法。
  10. 　前記マッピングに失敗した頻度を計算することは、
    　前記プロセッサにより、前記リード配列データの末端のペアを整合ペア、非整合ペア、不完全マップペア及びマッピング不能ペアのうちいずれかに分類することを含み、
    　前記不完全マップペアは、前記ペアの一方のみで前記マッピングに成功したペアであり、
    　前記プロセッサは、前記マッピングに失敗した頻度を、前記整合ペアの頻度と前記不完全マップペアの頻度の合計に対する前記不完全マップペアの相対頻度として計算する請求項１に記載の個体ゲノムの構造変異検出方法。
  11. 　前記プロセッサにより、塩基読み取りエラー率の入力を受け付けることをさらに含み、
    　前記プロセッサは、前記塩基読み取りエラー率の条件下で前記マッピングに失敗する確率を計算する請求項１に記載の個体ゲノムの構造変異検出方法。
  12. 　前記プロセッサにより、前記標準ゲノム配列データの各位置において、前記各位置を始点とする部分配列が前記標準ゲノム配列データ内で１個以下のミスマッチを許す条件下において一意的になる最小の長さを求めることをさらに含み、
    　前記プロセッサは、前記塩基読み取りエラー率と前記最小の長さとを用いた確率モデルに基づいて、前記マッピングに失敗する確率を計算する請求項１１に記載の個体ゲノムの構造変異検出方法。
  13. 　前記プロセッサにより、前記マッピングに失敗する確率についての所定の誤り率の入力を受け付けることをさらに含み、
    　前記プロセッサは、前記マッピングに失敗する確率が前記誤り率未満になるとき、前記マッピングに失敗した頻度が前記マッピングに失敗する確率と比較して有意に大きいと判定する請求項１に記載の個体ゲノムの構造変異検出方法。
  14. 　プロセッサと、前記プロセッサにより実行されるプログラムを格納するメモリと、を備える個体ゲノムの構造変異検出装置であって、
    　前記プロセッサは、
    　標準ゲノム配列データを受信する処理と、
    　前記標準ゲノム配列データの各位置において、変異がないと仮定した場合にマッピングに失敗する確率を計算する処理と、
    　個体ゲノムのリード配列データを受信する処理と、
    　前記リード配列データのそれぞれの末端を前記標準ゲノム配列データ上にマッピングする処理と、
    　前記標準ゲノム配列データの各位置で前記リード配列データの末端のマッピングに失敗した頻度を計算する処理と、
    　前記標準ゲノム配列データの各位置において、前記マッピングに失敗した頻度が前記マッピングに失敗する確率と比較して有意に大きいか否かを判定する処理と、
    　前記標準ゲノム配列データ上で、前記マッピングに失敗した頻度が有意に大きいと判定される位置が所定の長さ以上続く領域を、変異領域候補として抽出する処理と、
    　前記変異領域候補を構造変異の検出結果として出力する処理と、を実行する、個体ゲノムの構造変異検出装置。
PCT/JP2020/034166 2020-09-09 2020-09-09 個体ゲノムの構造変異検出方法及び装置 WO2022054178A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/JP2020/034166 WO2022054178A1 (ja) 2020-09-09 2020-09-09 個体ゲノムの構造変異検出方法及び装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2020/034166 WO2022054178A1 (ja) 2020-09-09 2020-09-09 個体ゲノムの構造変異検出方法及び装置

Publications (1)

Publication Number Publication Date
WO2022054178A1 true WO2022054178A1 (ja) 2022-03-17

Family

ID=80631406

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/034166 WO2022054178A1 (ja) 2020-09-09 2020-09-09 個体ゲノムの構造変異検出方法及び装置

Country Status (1)

Country Link
WO (1) WO2022054178A1 (ja)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015197899A (ja) * 2014-04-03 2015-11-09 株式会社日立ハイテクノロジーズ 配列データ解析装置、ｄｎａ解析システムおよび配列データ解析方法
JP2018533143A (ja) * 2015-08-06 2018-11-08 エイアールシー バイオ リミテッド ライアビリティ カンパニー ゲノム分析のためのシステムおよび方法
WO2020067603A1 (ko) * 2018-09-28 2020-04-02 한양대학교 산학협력단 다중 참조 유전체에 기반한 유전체 구조변이 검출 방법 및 구조변이 검출 장치
CN111326212A (zh) * 2020-02-18 2020-06-23 福建和瑞基因科技有限公司 一种结构变异的检测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015197899A (ja) * 2014-04-03 2015-11-09 株式会社日立ハイテクノロジーズ 配列データ解析装置、ｄｎａ解析システムおよび配列データ解析方法
JP2018533143A (ja) * 2015-08-06 2018-11-08 エイアールシー バイオ リミテッド ライアビリティ カンパニー ゲノム分析のためのシステムおよび方法
WO2020067603A1 (ko) * 2018-09-28 2020-04-02 한양대학교 산학협력단 다중 참조 유전체에 기반한 유전체 구조변이 검출 방법 및 구조변이 검출 장치
CN111326212A (zh) * 2020-02-18 2020-06-23 福建和瑞基因科技有限公司 一种结构变异的检测方法

Similar Documents

Publication Publication Date Title
US10192026B2 (en) Systems and methods for genomic pattern analysis
US10783984B2 (en) De novo diploid genome assembly and haplotype sequence reconstruction
US11756652B2 (en) Systems and methods for analyzing sequence data
US20200399719A1 (en) Systems and methods for analyzing viral nucleic acids
Hasan et al. Performance evaluation of indel calling tools using real short-read data
CN106068330B (zh) 将已知等位基因用于读数映射中的系统和方法
US9165109B2 (en) Sequence assembly and consensus sequence determination
US20130138358A1 (en) Algorithms for sequence determination
US11646101B2 (en) Continuous wavelet-based dynamic time warping method and system
US10810239B2 (en) Sequence data analyzer, DNA analysis system and sequence data analysis method
US20150169823A1 (en) String graph assembly for polyploid genomes
US20150286775A1 (en) String graph assembly for polyploid genomes
US20220254444A1 (en) Systems and methods for detecting recombination
CN106529211A (zh) 变异位点的获取方法及装置
EP1859268A2 (en) System, method and computer program for non-binary sequence comparison
US20210350876A1 (en) System and method for direct subsequence searching and mapping in nanopore raw signal
WO2022054178A1 (ja) 個体ゲノムの構造変異検出方法及び装置
US20220157401A1 (en) Method and system for mapping read sequences using a pangenome reference
WO2016205767A1 (en) String graph assembly for polyploid genomes
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
US20170024514A1 (en) Distance maps using multiple alignment consensus construction
Oja et al. Clustering of human endogenous retrovirus sequences with median self-organizing map
Jiang et al. Long-read based novel sequence insertion detection with rCANID
CN110570908A (zh) 测序序列多态识别方法及装置、存储介质、电子设备
Matroud Nested tandem repeat computation and analysis

Legal Events

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

Ref document number: 20953244

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20953244

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP