JP2017016665A - Method for selecting variation information from sequence data, system, and computer program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To acquire variation information having high reliability from lead information derived from a DNA sample.SOLUTION: A threshold of lead depth is determined by finding that coincidence between genome variation information based on lead information acquired from a next generation sequencer and genome variation information acquired by another means such as an SNP array is closely associated with the depth of lead information (lead depth) and using the coincidence as an index, and it is found that genome variation information having high reliability can be acquired by filtering the genome variation information on the basis of the lead information.SELECTED DRAWING: Figure 1

Description

  The present invention relates to an organism mutation information selection method using a computer, a system for performing the method, and a computer program.

  Prior to the introduction of so-called next-generation sequencers (NGS), microarrays were mainly used for comprehensive genome analysis. However, when microarrays are used, it is difficult to completely detect unknown genomic mutation information. This is because the probe must be designed in advance.

  Therefore, establishment of means for detecting genomic mutation information by using NGS is expected in the industry.

  However, with respect to the lead information provided by NGS, there is a risk such as an error in mapping.

  Therefore, it has been a problem to find highly reliable genome mutation information by eliminating such risks from the read information.

  The present inventor found that the degree of coincidence between genome mutation information based on read information obtained from NGS or the like and genome mutation information obtained by other means such as a SNP array is the depth of the read information (read depth). I found it closely related. Then, the inventor determines a read depth threshold by using the degree of coincidence as an index, and obtains genome mutation information with high reliability by filtering genome mutation information based on the read information. The present invention has been completed by finding out what can be done. Furthermore, the present inventor has found that the reliability of genome mutation information from the read information can be improved by using other filtering means in addition to the basic filtering means.

The present invention relates to a method for selecting gene mutation information by a computer, wherein the following steps (1) to (3) are performed. Method).
(1) A step of performing the following for each DNA sample.
(1-1) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism is grouped into a subset having a predetermined range of the read depth. Divided,
(1-2) Mutation information between the primary mutation set information and the secondary mutation set information obtained from each DNA sample by means other than mapping, based on the measured value of the genotype coincidence of the mutation information of each subset. By comparing genotypes of
(1-3) determining an appropriate range of read depths by combining the range of read depths of subsets having a genotype match measurement value higher than a predetermined threshold;
(2) For each DNA sample, from the primary mutation set information in step (1-1), removing mutation information having a read depth outside the appropriate range of read depth to obtain remaining mutation information for each DNA sample .
(3) As a result of performing the above steps (1) and (2) for at least one DNA sample among all DNA samples, the remaining mutation information is extracted together with the locus to which the mutation information belongs. And selecting as target mutation information.

  In the present invention, a “DNA sample” is a sample containing DNA of a test subject (target organism), and the DNA includes genomic DNA, cDNA and the like. The organism that is the source of the DNA sample is not particularly limited as long as the organism has a genome, and includes mammals such as humans, reptiles, insects, plants, bacteria, and viruses. Of these, humans are preferred.

  The term “mutation” means a variation in the nucleotide sequence of a genome, and specifically includes a single nucleotide variation (SNV), insertion, or deletion. “Mutation information” is digital information representing these.

  The term “read information” is digital information about a nucleotide sequence obtained by processing a genome with a sequencer. The lead information is divided into two types of information. The one is single read information composed of single sequence information, and the other is paired end read information including sequence information of two adjacent regions on the genome. In the present invention, the lead information includes both types of lead information. Sequencers that provide read information are “high-performance sequencers” including a so-called “next-generation sequencer” that can provide a large amount of information in a relatively short time. Examples of next-generation sequencers at the time of the present invention include Genome Sequence FLX (Roche (454)); Genome Analyzer IIx, HiSeq2000, HiSeq2500, and MiSeq (Illumina); SOLiD (Applied Biosystem); Examples include, but are not limited to, PacBio RS II (Pacific Biosciences) and include all sequencers provided at the present time and in the future.

  The term “mapping algorithm” used in generating primary mutation set information is an algorithm for specifying a position in a reference sequence—from which a read sequence is read—for example, BWA-MEM, and Bowtie2 is included. In the present specification, claims, and drawings, “mapping” is synonymous with “alignment”. The “mutation call algorithm” is an algorithm for extracting lead information in which mutation information is found from the mapped lead information, and includes, for example, Bcftools and Genome Analysis Toolkit (GATK). By applying this mutation call algorithm to the mapped read information, the mutation information is extracted, and “reference homo”, “hetero”, and “non-reference homo” are classified. “Reference homo” is a form in which the same base as the standard genomic sequence (or reference genomic sequence) exists in a predetermined site in a homozygous genotype, and “non-reference homo” is different from the standard genomic sequence. In this form, the base exists in a homozygous genotype. “Hetero” is a form in which bases different from the same base as the standard genome sequence exist in a heterozygous genotype.

  The term “read depth” is synonymous with coverage, meaning the number of corresponding reads for each mutation site (locus), calculated during mapping of the DNA sample, where duplicate site readings are tested This is an index indicating how much the target mutation site has been performed. In this specification, it may be simply described as “depth”.

  “Means other than mapping” for obtaining “secondary mutation set information” are not particularly limited, and include, for example, Sanger sequence, mass array (registered trademark), and SNP array, and to obtain primary mutation set information. It has at least the same accuracy as the sequencer used, preferably higher accuracy. Among these means, the SNP array is preferable because it is inexpensive and can determine a large number of mutation sites with high accuracy at a time. In this case, “many with high accuracy at one time” preferably corresponds to about 300,000 to about 3,000,000 SNPs.

  In the “secondary mutation set information”, the above-mentioned classification of “reference homo”, “hetero”, and “non-reference homo” is assigned by a method according to a specific “means other than mapping”. . For example, in the case of an SNP array, a clustering algorithm called “genotype call” is applied based on the color of the probe dye or the fluorescence of the fluorescent dye, so that the above “reference homo”, “hetero”, and “ “Non-reference homo” is classified. In the case of the Sanger sequence, the “reference homo”, “hetero”, and “non-reference homo” classifications are analyzed by analyzing the waveform of the chart based on the color of the probe dye and the fluorescence of the fluorescent dye. Is known to be made.

  The “measured value of genotype coincidence” in the step (2-2) is preferably (1) precision and / or (2) detection rate (power). That is, both the precision and the detection rate, or one of them can be used as the measurement value. Specifically, a predetermined threshold value is set for these measurement index values, and the measured value is equal to or higher than the threshold value of the accuracy rate or the detection rate, or the threshold values of both the accuracy rate and the detection rate. The above case can be recognized as “the genotype matching degree is higher than a predetermined threshold”. As a preferred aspect, if either the precision ratio or the detection ratio is used, it is preferable to prioritize the precision ratio as a threshold index, and it is more preferable that both the precision ratio and the detection ratio be threshold indices. It is. Here, these measurement index values will be described.

<Indicator values (about precision and detection rate)>
These index values are values that can be calculated by the “mutation call” performed in step (1). By executing the mutation call algorithm, the non-matching information corresponding to the “non-matching information whose lead information does not match the reference sequence information”, which is recognized when mapping the lead information or using the SNP array performed in advance, is performed. Is a mismatch due to an error or a mismatch due to a mutation that the DNA sample truly has. Here, if the non-matching information is not “mismatch due to an error” and it is determined by executing the mutation call algorithm, it is determined that there is a mutation in the portion where the non-matching has occurred. Furthermore, in diploid organisms such as humans, when all the read information of the site where the mutation is identified is non-reference type, it is “non-reference homozygous”, which means that there is a mutation in DNA from both parents It is regarded as having done. In addition, when the reference type and the non-reference type are mixed in the read information of the site where the mutation is identified, it is “hetero”, and it is considered that the mutation exists in one of the parents. Further, even when the mutation call algorithm does not determine that the mutation is present, it is possible to regard the site where all the read information is of the reference type as being in the reference homo state.

  In summary, for example, when the target organism is a human, the non-matching information when the mutation call algorithm is executed includes (a) reference homo, (b) hetero, and (c) non-reference homo. Classification is possible, and when the primary mutation set information is combined with the secondary mutation set information, nine cases are recognized.

  In Table 1 below, these nine patterns are illustrated, and the precision and detection rate calculated based on these are described. The result of the mutation call based on NGS generates primary mutation set information, and the result of the mutation call based on SNP array generates secondary mutation set information. The alphabets (ai) in the table are the number of genotypes detected as mutation information in each subset. For example, “d” in Table 1 below is the number of genotypes determined to be “hetero” by NGS and “reference homo” by the SNP array.

  Here, the relevance rate is “primary mutation set information (NSG in this example) and determined as“ hetero or non-reference homo ”, and the same determination is also made in secondary mutation set information (SNP array in this example). Is a ratio represented by “(e + i) / (d + e + f + h + i)” in Table 1. The detection rate is the ratio of “secondary mutation set information (SNP array in this example) determined as“ hetero or non-reference homo ”, and the same determination was made in primary mutation set information (NSG in this example)” In Table 1, the value is represented by “(e + i) / (b + c + e + f + h + i)”. The display format of both the precision rate and the detection rate is not limited, and for example, direct ratio value display or percentage display may be used.

  The threshold value of the precision in the step (2-3) is preferably 0.995 or more, and more preferably 0.998 or more. Similarly, the threshold of the detection rate is preferably 0.93 or more, and more preferably 0.95 or more.

Further, in the selection method of the present invention, the above-mentioned steps (1) to (4) are carried out to obtain the first specific mutation information, which is excluded in one or more steps selected from the following (A) to (D). One or more mutation information defined as power mutation information is removed from the first specific mutation information to obtain the second specific mutation information, and the second specific mutation information is selected as the target mutation information. It is also possible to adopt an aspect.
(A) When the first specific mutation information is obtained by using two or more DNA samples, the ratio of the specific mutation information having a read depth outside the appropriate range of the read depth is the step (3). If more than 5% to 20% in the mutation site in the, the specific mutation information is defined as the mutation information to be excluded above;
(B) All the mutation information existing in the area assigned to the low complexity area is selected from the primary mutation set information, and the measured value of the genotype coincidence is the all mutation information and each low complexity area. All mutation information in the low-complexity region that is calculated between the mutation information in the second mutation set information and has a genotype-matched measurement value below a predetermined threshold is defined as the mutation information to be excluded. ;
(C) Step (1) or steps (1) to (4) are performed by using a different mapping algorithm and different mutation call algorithm from the first step (1), and are different for each DNA sample. Mutation information to obtain primary mutation set information or different mutation information and is included in the first mutation information but not in different primary mutation set information or different mutation information is defined as mutation information to be excluded;
(D) An algorithm related to the Hardy-Weinberg equilibrium test is performed, and mutation information deviating from the Hardy-Weinberg equilibrium is defined as mutation information to be removed.

  These additional steps (A) to (D) are four types shown as steps (A), (B), (C) and (D) on the premise that the selection method of the present invention described above is performed. One of the other filtering means is further performed, or two to four of them are combined and performed. By adding these additional steps (A) to (D) to the selection method of the present invention, it is possible to obtain mutation information having higher reliability derived from lead information. The “specific mutation information” in the additional step (A) means, for example, mutation information such as individual SNV.

  Furthermore, the present invention provides a computer system for performing the above-described selection method of the present invention, a computer program including an algorithm for causing a computer to execute the selection method of the present invention, and a computer reading the computer program of the present invention. Provide digital media, stored as possible.

  The present invention also provides mutation information obtained by the selection method of the present invention.

Furthermore, the present invention provides a method for determining an appropriate range of the lead depth by a computer, which includes the following steps (1) to (3) (hereinafter also referred to as the depth determination method of the present invention).
(1) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, the primary mutation set information, and the predetermined range of the read depth Group into subsets with
(2) The gene of the mutation information between the primary mutation set information and the secondary mutation set information obtained from each DNA sample by means other than mapping, using the measured value of the genotype coincidence of the mutation information of each subset. Calculate by comparing types,
(3) The appropriate range of the read depth is determined by combining the read depth ranges of the subsets having the genotype matching degree higher than the predetermined threshold.

  The individual definitions of the depth determination method of the present invention are as described in the selection method of the present invention. For example, the “measured value of genotype coincidence” is preferably the above-described “precision and / or detection rate (power)”.

  Furthermore, the present invention provides a computer system or a computer apparatus for performing the depth determination method of the present invention, and provides a computer program including an algorithm for causing a computer to execute the depth determination method of the present invention. Provide digital media stored in a computer readable state.

  According to the present invention, highly reliable mutation information can be obtained from lead information derived from a DNA sample.

It is drawing which shows the whole image of the SNV filter of this invention.

  Specifically, the above-described selection method and depth determination method of the present invention are executed in the computer system of the present invention, and are executed by causing the computer program algorithm of the present invention to be realized in a computer. With these forms, highly reliable mutation information can be easily obtained from lead information derived from a DNA sample.

[Computer System of the Present Invention]
The computer system or computer apparatus is mainly a computer system for performing each step of the selection method of the present invention in a computer. A computer system is substantially synonymous with a computer device. A computer system or apparatus usually includes a recording unit (HD), an arithmetic processing unit (CPU), an operation unit (keyboard), and a display unit (display). In the recording unit, electronic information forming a basis for performing the selection method of the present invention is recorded, and the processing of the necessary steps by the arithmetic processing unit is performed on the electronic information read out to the temporary recording unit (memory). Done, and the result is also recorded in the recorder. In the arithmetic processing unit, image data for promoting the operation of the operation unit or displaying the processing result is created, and the image data is displayed on the display unit via the video RAM. In this way, the process of advancing each process is performed, and finally necessary mutation information is derived. The same applies to the computer system or apparatus for performing the depth determination method of the present invention, and the processing of each step for determining the appropriate lead depth is performed in the computer, whereby a desired appropriate lead depth is derived.

A main aspect of the computer system or apparatus of the present invention is a computer system that selects specific mutation information from base sequence data, and includes at least a recording unit and an arithmetic processing unit, and includes the following processes (a) to (f): Is executed by a computer system or apparatus.
(A) Based on the mapping algorithm and the mutation call algorithm, the recording unit maps the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and mapping from each DNA sample Secondary mutation set information obtained by means other than is recorded,
(B) In the arithmetic processing unit, a process of grouping primary mutation set information read from the recording unit into a subset having a predetermined range of read depth is executed,
(C) A genotype comparison process of the mutation information between the primary mutation set information and the secondary mutation set information read from the recording unit as the measured value of the genotype matching degree of the mutation information of each subset Is calculated by executing
(D) The determination process of the appropriate range of the read depth is performed by combining the range of the read depth of the subset having the genotype matching degree higher than the predetermined threshold into one,
(E) From the primary mutation set information for each DNA sample, mutation information having a read depth outside the appropriate range of the read depth determined in (d) above is executed, and the remaining mutation information for each DNA sample Is extracted,
(F) As a result of performing the above processing (a) to (e) for at least one DNA sample among all DNA samples, the remaining mutation information is extracted together with the gene locus to which the mutation information belongs. , A process for identifying the genetic information as target specific mutation information.

  The said process (a)-(f) is a process corresponding to each process (1)-(3) in the selection method of this invention. That is, process (a) and (b) correspond to process (1-1), process (c) corresponds to process (1-2), and process (d) corresponds to process (1-3). The process (e) corresponds to the process (2), and the process (f) corresponds to the process (3).

  As described above, the “primary mutation set information” in the process (a) is mutation information to which a read depth is assigned for each DNA sample of one or more individuals of the target organism based on the mapping algorithm and the mutation call algorithm. is there. The “mapping algorithm” and the “mutation call algorithm” are as described above. The means for securing the “primary mutation set information” is not limited. For example, in the same computer terminal, from the lead information, those generated based on the “mapping algorithm” and the “mutation call algorithm” may be used, “Primary mutation set information” separately generated elsewhere may be used for the computer system of the present invention. Furthermore, as for “secondary mutation set information” obtained by means other than mapping such as an SNP array, the securing means is not limited at all.

  The process (a) is a process corresponding to “selection procedure 1” in an embodiment described later.

  The grouping process (process (b)) into subsets for each read depth range, and the genotype comparison process (process (c)) of the mutation information between the primary mutation set information and the secondary mutation set information, This is processing corresponding to “selection procedure 2 (the middle graph of“ procedure 2 ”in FIG. 1)” in an embodiment described later. As described above, the “measured value of the genotype coincidence of mutation information of each subset” in the processing (c) is preferably a precision and / or detection rate (power), more preferably As the precision, it is most preferably used as the precision and the detection rate. The threshold value of the precision is preferably 0.995 or more, and more preferably 0.998 or more. The detection rate threshold is preferably 0.93 or more, and more preferably 0.95 or more.

  Next, the determination process of the appropriate range of the read depth (process (d)) and the exclusion process of the mutation information having the read depth outside the appropriate range of the read depth (process (e)) are performed in the “selection” example described later. This is a process corresponding to the procedure 2 ("procedure 2" in FIG. 1).

  Furthermore, an example of the result of the target specific mutation information specifying process (process (f)) is shown in “Procedure 2” in Table 2 of Examples described later.

In the computer system of the present invention, the specific mutation information obtained by the execution of the computer system should be excluded from one or more processes selected from the following (α) to (δ) as the first specific mutation information. One or more mutation information defined as mutation information is removed from the first specific mutation information to obtain second specific mutation information, and the second specific mutation information is used as target mutation information. It is preferable that the computer system executes the extracted process.
(Α) When the first specific mutation information is obtained by using two or more DNA samples, the ratio of the specific mutation information having a read depth outside the appropriate range of the read depth is determined by the above processing (f The processing in which the specific mutation information is defined as the mutation information to be excluded when it exceeds 5% to 20% of the mutation information obtained by (1).
(Β) All the mutation information existing in the area assigned to the low complexity area is extracted from the primary mutation set information, and the measured value of the genotype coincidence is the all mutation information and each low complexity area. All mutation information in the low-complexity region that is calculated and processed with the mutation information in the second mutation set information regarding and has a genotype-matched measurement value below a predetermined threshold is defined as the mutation information to be excluded. Processing.
(Γ) processing different primary mutation set information or different mutation information for each DNA sample generated by using a different mapping algorithm and different mutation call algorithm from the generation of the first used primary mutation set information ( The process in which the mutation information that is included in the mutation information extracted first by f) but not included in the different primary mutation set information or the different mutation information is defined as the mutation information to be removed.
(Δ) A process in which an algorithm related to the Hardy-Weinberg equilibrium test is performed, and mutation information deviating from the Hardy-Weinberg equilibrium is defined as mutation information to be removed.

  The additional process α relating to “removal of mutation information based on the read depth for each mutation site” is a process corresponding to the additional step (a) of the selection method of the present invention, and “selection procedure 3” (see FIG. 1 corresponds to “Procedure 3”, particularly in the middle), and an example of the result is described in “Procedure 3” in Table 2.

  The additional process β related to “removal of mutation information based on genome complexity” is a process corresponding to the additional step (b) of the selection method of the present invention, and “selection procedure 4” (see FIG. This corresponds to “Procedure 4”, in particular, the table on the left side and the right figure.

  The additional process γ related to “removal of mutation information related to bias of analysis method” is a process corresponding to the additional step (c) of the selection method of the present invention, and “selection procedure 5” (see “ An example of the result is described in “Procedure 5” in Table 2.

  The additional process (δ) relating to “removal of population genetic variation information” is a process corresponding to the additional step (d) of the selection method of the present invention. Specifically, a Hardy-Weinberg equilibrium test is performed. The mutation information deviating from the equilibrium is removed. This additional processing (δ) corresponds to “selection procedure 6” (“procedure 6” in FIG. 1) of an embodiment described later, and an example of the result is described in “procedure 6” of Table 2.

The computer system or apparatus of the present invention can be an embodiment for executing not only the selection method of the present invention but also the depth determination method of the present invention. Specifically, the computer system selects a specific read depth range from base sequence data, and includes at least a recording unit and an arithmetic processing unit, and performs the following processes (a) to (d). A featured computer system or apparatus.
(A) Based on the mapping algorithm and the mutation call algorithm, the recording unit maps the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and mapping from each DNA sample Secondary mutation set information obtained by means other than is recorded,
(B) In the arithmetic processing unit, a process of grouping primary mutation set information read from the recording unit into a subset having a predetermined range of read depth is executed,
(C) A genotype comparison process of the mutation information between the primary mutation set information and the secondary mutation set information read from the recording unit as the measured value of the genotype matching degree of the mutation information of each subset Is calculated by executing
(D) The process of determining the appropriate range of the read depth is executed by combining the range of the read depth of the subset having the genotype matching degree higher than the predetermined threshold.

The processes (a) to (d) are substantially the same as (a) to (d) among the processes (a) to (f) in the computer system or apparatus for executing the selection method of the present invention. is there.
[Computer program of the present invention]
A first aspect of the computer program of the present invention is a computer program including an algorithm for executing the selection method of the present invention. Specifically, the computer program selects specific mutation information from base sequence data. The computer program includes an algorithm for causing a computer to realize the following functions (a) to (f).
(A) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and obtained by means other than mapping from each DNA sample A function to read out the information from the recording section where the secondary mutation set information is recorded.
(B) A function of performing a process of grouping primary mutation set information read from the recording unit into subsets having a predetermined range of read depth.
(C) A process for comparing genotypes of mutation information between primary mutation set information and the secondary mutation set information read from the recording unit, based on a measurement value of the genotype matching degree of mutation information of each subset A function that calculates by executing
(D) a function of executing a process of determining an appropriate range of read depth by combining the range of read depths of a subset having a measurement value of genotype matching higher than a predetermined threshold;
(E) From the primary mutation set information for each DNA sample, mutation information having a read depth outside the appropriate range of the read depth determined in (d) above is executed, and the remaining mutation information for each DNA sample The function to extract.
(F) As a result of executing the above processes (a) to (e) for at least one DNA sample among all DNA samples, the remaining mutation information is extracted together with the locus to which the mutation information belongs. The function of identifying the genetic information as the specific mutation information of interest.

  The functions (a) to (f) correspond directly to the “processing (a) to (f)” in the computer system or apparatus of the present invention. That is, the computer program of the present invention includes an algorithm that allows the computer to perform the processes (a) to (f) by realizing the functions (a) to (f) in the computer. .

  The algorithm for realizing the function (a) is as described as the function (a).

  As an algorithm for realizing the function (b), for example, an algorithm for classifying each read depth value assigned to the primary mutation set information can be cited, and the degree of the classification is not particularly limited, and for each “1”. It is possible to do this, and it is also possible to do it with more values. In the example described later, the classification was performed for each “1”.

  The algorithm for realizing the function (c) is, for example, the above-described individual mutation information for each subset of the mutation information classified by read depth based on the primary mutation set information obtained by the function (b). Classification processing into “reference homo”, “hetero”, and “non-reference homo”, and secondary mutation set information also into “reference homo”, “hetero”, and “non-reference homo” An algorithm for performing classification processing and calculating / summing (averaging) the above-mentioned “adaptation rate” and / or “detection rate” is given.

  The algorithm for realizing the function (d) is, for example, calculated with respect to a predetermined threshold (as described above) for “accuracy rate” and / or “detection rate” and a subset of variation information classified by read depth. Comparing the summed (averaged) values of “Precision rate” and / or “Detection rate”, extracting a subset for each lead depth that has a value exceeding the threshold, and integrating the read depth of the subset Then, an algorithm for generating an appropriate range of the lead depth can be mentioned.

  The algorithm for realizing the function (e) is as described as the function (e).

  The algorithm for realizing the function (f) is as described as the function (f).

In addition to the algorithm for realizing the functions (a) to (f) in the computer program, the computer program of the present invention further includes one or more functions selected from the following (α) to (δ). An algorithm that realizes a function of obtaining new second specific mutation information by removing one or more pieces of mutation information defined as the mutation information to be added from the existing first specific mutation information can be included.
(Α) When the first specific mutation information is obtained by using two or more DNA samples, the ratio of specific mutation information having a read depth outside the appropriate range of the read depth is the function (f The function of defining the specific mutation information as the above-described mutation information to be excluded when exceeding 5% to 20% of the mutation information obtained by (1).
(Β) All the mutation information existing in the area assigned to the low complexity area is extracted from the primary mutation set information, and the genotype match value is measured as the all mutation information and each low complexity area. All mutation information in the low-complexity region having a genotype-matched measurement value below a predetermined threshold is defined as the mutation information to be excluded. function.
(Γ) Different primary mutation set information or different mutation information for each DNA sample generated by using a different mapping algorithm and different mutation call algorithm from the generation of the first used primary mutation set information. A function of defining, as the mutation information to be excluded, the mutation information that is included in the mutation information extracted first in (f) but is not included in the different primary mutation set information or the different mutation information.
(Δ) A function for defining the mutation information deviating from the Hardy-Weinberg equilibrium as the mutation information to be removed by performing an algorithm related to the Hardy-Weinberg equilibrium test.

  The additional functions (α) to (δ) correspond to “addition processing (α) to (δ)” as they are in the computer system or apparatus of the present invention. That is, the computer program of the present invention includes an algorithm that allows the computer to perform the processes (α) to (δ) by realizing the additional functions (α) to (δ) in the computer. Yes.

  The algorithm for realizing the additional function (α) is, for example, extracting the read depth of specific mutation information (such as SNV information) in the first specific mutation information for each DNA sample, and the read depth is the above appropriate Calculate the ratio of the number of DNA samples outside the range to the total number of DNA samples, and determine whether the calculated value exceeds a specific threshold (any of 5 to 20% as a lower limit), In the case of exceeding, there is an algorithm for extracting the specific mutation information as “mutation information to be removed”.

  As an algorithm for realizing the additional function (β), for example, mutation information in a predetermined low complexity region is extracted in a primary mutation set and a secondary mutation set, and “matching rate” and / or “detection rate” is determined. Compared with the predetermined threshold (as described above) and the value of the “compliance rate” and / or “detection rate” calculated and summed (averaged) for the subset of mutation information divided by read depth, There is an algorithm for extracting variation information of “matching rate” and / or “detection rate” below the threshold as “mutation information to be removed” regardless of the read depth.

  The algorithm for realizing the additional function (γ) is, for example, “for each DNA sample generated by using a different mapping algorithm and a different mutation call algorithm from the generation of the first used primary mutation set information. Match or mismatch between “different primary mutation set information or different mutation information” and “mutation information first extracted by the function (f)” is extracted, and the mismatch mutation information is first extracted by the above function (f). Algorithm to extract “mutation information to be removed from the mutation information that has been performed”.

  The algorithm for realizing the additional function (δ) is as described as the additional function (δ).

The computer program of the present invention can be an algorithm for realizing not only the selection method of the present invention but also the depth determination method of the present invention in a computer. Specifically, it is a computer program for calculating an appropriate range of read depth from base sequence data, and includes an algorithm that causes the computer to realize the following functions (a) to (d): Is a computer program.
(A) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and obtained by means other than mapping from each DNA sample A function to read out the information from the recording section where the secondary mutation set information is recorded,
(B) a function of performing a process of grouping primary mutation set information read from the recording unit into a subset having a predetermined range of read depths;
(C) A genotype comparison process of the mutation information between the primary mutation set information and the secondary mutation set information read from the recording unit as the measured value of the genotype matching degree of the mutation information of each subset Function to calculate by executing
(D) a function of executing a process of determining an appropriate range of read depths by combining the range of read depths of a subset having a genotype matching degree measurement value higher than a predetermined threshold;

  The functions (a) to (d) are substantially the same as the functions (a) to (d) among the functions (a) to (f) in the computer program for executing the selection method of the present invention.

  The computer program of the present invention can be written not only in a high-level computer language such as C language, Java (registered trademark), Perl, or Python, but also in binary numbers or in assembly language. The computer program is recorded in the recording unit as needed or in advance or in external hardware, and arithmetic processing according to the recorded algorithm is performed in the arithmetic processing unit as necessary.

Examples of the present invention will be described below.
[ Sequencing of whole genome ]
In order to prevent mixing of samples, genomic DNA specimens were handled using a 96-well plate for library construction. Genomic DNA was prepared to a concentration of 20 ng / μL using a laboratory automation system (Biomek NXP, Beckman) and averaged to 550 bp using a 96-well plate DNA sonicator (Covaris LE220, Covaris). Fragmented. The DNA thus obtained was converted into a library using a TruSeq DNA PCR-Free HT sample prep kit (Illumina) using a Bravo liquid handling apparatus (Agilent Technologies). Finally, the prepared library was transferred to a 1.5 mL test tube identified by a barcode label and subjected to denaturation treatment for library quality control.

  Library quantification and quality control were performed with the newly developed quantitative MiSeq method (qMiSeq). In this procedure, 8 μL or 10 μL of the prepared library is first denatured with an equal amount of 0.1N NaOH for 5 minutes and then diluted with 49 times volume of ice-cooled “Illumina HT1 buffer”. Next, 96 50 μL diluted libraries containing 3 target samples previously analyzed with HiSeq are mixed. Of this, 60 μL is diluted with 540 μL of ice-cold Illumina HT1 buffer and analyzed using the MiSeq 25 bp paired end protocol. The index ratio determined by the analysis of the MiSeq sequencer is used as a relative concentration for determining the execution condition when executing HiSeq. Details of the qMiSeq method are already known as method papers [Katsuoka F, et al. An efficient quantitation method of next-generation sequencing libraries by using MiSeq sequencer. Anal. Biochem. 466c, 27-29 (2014)]. In addition to the MiSeq QC, Fragment Analyzer software (Advanced Analytical) was also used for quality control of the DNA library.

The aforementioned DNA library was used for analysis with a HiSeq 2500 sequencer according to the manufacturer's protocol. The analysis was performed with a Rapid-Run Mode of 162 bp pair end using TruSeq Rapid PE Cluster Kit (Illumina) and one-and-a-half TruSeq Rapid SBS Kit (200 cycles, Illumina). Based on the result of qMiSeq, the library was adjusted to an appropriate concentration, and on-board cluster formation (Illumina) was performed. The cluster density reported first was confirmed, and the operation of determining whether to continue the analysis depending on the density (about 550 to 650 K / mm ² ) was repeated.

[ Selection Procedure 1 ] Alignment and Mutation Call of Each DNA Sample The sequence read of each genomic DNA sample (hereinafter also referred to as a DNA sample) output from the HiSeq sequencer is used as a decoy sequence (hs37d5) as a reference genomic sequence (GRCh37 / hg19). ). Alignment was performed with Bowtie 2 (version 2.1.0) using the “-X 2000” option and BWA-MEM (ver. 0.7.5a-r405) using the standard option. Mutation calls were made using “Bcftools software” (ver. 0.1.17-dev) and “the Genome Analysis Toolkit” (GATK version 2.5-2) using these alignment results (“FIG. 1“ GATK version 2.5-2 ”). Procedure 1 "). For use in subsequent filtering, the read depth of all SNV sites of each DNA sample was calculated (the bottom row in “Procedure 1” in FIG. 1). Here, it is confirmed that the lead depth is the number of sequence leads that are aligned with a mapping quality of 5 or more at each SNV position.

  Next, in order to evaluate the precision and detection rate (= recall) of the sites designed in the SNP array, the same sample as the SNV call was analyzed, “HumanOmni2.5 A comparison was made with the genotype call of “-8 BeadChip”. The SNV (Bowtie2 + Bcftools) match rate obtained using the Bcftools mutation caller for the Bowtie2 alignment is always SNV (BWA-MEM + GATK) obtained by fitting the GATK mutation caller to the BWA-MEM alignment. ) Was higher. On the other hand, the detection rate of BWA-MEM + GATK was higher than that of Bowtie2 + Bcftools. This result is consistent with previous studies [Li H. Toward better understanding of artifacts in variant calling from high-coverage samples. Bioinformatics 30, 2843-2851 (2014)]. From these results, it can be considered that the SNV obtained by “BWA-MEM + GATK” is a result aimed at high sensitivity, that is, finding as many SNV candidates as possible. On the other hand, the SNV obtained by “Bowtie2 + Bcftools” has high reliability, and is used as a highly reliable SNV candidate in the following filtering step. Unless otherwise specified, “SNV” refers to a variant that has been mutated and called by “Bowtie2 + Bcftools”.

[ Selection Procedure 2 ] SNV Filtering Based on Read Depth for Each DNA Sample Genotypes with extremely high or low read depth were filtered for each DNA sample. For each DNA sample, the range of the read depth that leaves the genotype was determined based on the precision and detection rate calculated from the comparison with the genotyping results of “HumanOmni2.5-8 BeadChip”. SNVs were grouped according to the lead depth of NGS, and the precision and detection rate were calculated for each group (middle graph of “Procedure 2” in FIG. 1). At that time, only reads with a mapping quality of 5 or higher were considered in order to eliminate the influence of reads aligned to multiple chromosomes. As a result, a genotype in the range of read depths of SNVs belonging to a group having a fitness greater than 0.998 and a detection rate greater than 0.95 is left as a reliable SNV (see “Procedure 2” in FIG. 1). "). In this selection procedure 2, 2% of the separation sites were removed ("Procedure 2" in Table 2).

  The “known mutation” in Table 2 is the SNV reported in “dbSNP build 138”.

[ Selection Procedure 3 ] SNV Filtering Based on Read Depth for Each SNV Site As a result of applying the filtering based on the read depth in the selection procedure 2, an SNV in which genotypes of 10% or more of all DNA samples do not exist The whole was removed. For example, sites that lead to extremely high or low read depths as a result of aligning reads to repetitive sequences may be subject to this filtering. The percentage of SNV sites remaining by this filtering was calculated (middle of “Procedure 3” in FIG. 1). In this selection procedure 3, 4.99% of SNV sites were removed ("Procedure 3" in Table 2 above).

[ Selection Procedure 4 ] Filter based on genome complexity It is generally difficult to accurately call mutations in regions with low sequence complexity such as tandem repeats. The RepeatMasker program [Smit, AFA, Hubley, R & Green, P. RepeatMasker Open-4.0. 2013-2015 <http://www.repeatmasker.org>.] Has 56% of the genomic region with such low complexity. Judged to be an area. For each group of repeats determined by Repeat Masker, genotypic fitness with SNP array genotyping of the same sample was calculated (the table on the left side of “Procedure 4” in FIG. 1). The SNV call fitness of each genomic region was calculated based on the result of genotype call obtained in the analysis of “HumanOmni2.5-8 BeadChip”. 1) (the right figure of “Procedure 4”) was removed. Low complexity genomic regions include Alu, ERVK, Low_complexity, Satelite, Simple_repeat, TcMar-Mariner, CR1, DNA, Deu, Dong-R4, ERV, ERV1, ERVK, ERVL, ERVL-MaLR, GENCODEHel, Gyps, Elps , L1, L2, LTR, MIR, Merlin, MuDR, Penelope, PiggyBac, RNA, RTE, RTE-BovB, SINE, TcMar, TcMar-Tc2, TcMar-Tigger, acro, centr, hAT, hAT-BlackAThack HAT-Tip100, intelligent, rRNA, scRNA, snRNA, srpRNA, tRNA, elo and the like. This procedure removed 14.22% of the SNV (“Procedure 4” in Table 2 above).

[ Selection procedure 5 ] Filter against bias in analysis method In order to control the bias in the analysis method in the selection procedure 1, only SNV sites found in different analysis methods were left. Specifically, as described above, “BWA-MEM + GATK” in which high-sensitivity detection is performed is used as “another method”, and SNVs that are not found in the “other method” are additionally excluded. For example, in the upper diagram of “Procedure 5” in FIG. 1, the leftmost SNV was not found by other methods (“BWA-MEM + GATK”), so this SNV was excluded from the subsequent analysis. This procedure removed 0.57% of the SNV (“Procedure 5” in Table 2 above).

[ Selection Procedure 6 ] Population Genetic Filtering SNVs having a p value of less than 10 ⁻⁵ were removed by Hardy-Weinberg equilibrium test, which was calculated from the SNV genotype frequency obtained in the above selection procedure 5. This filter is intended to remove SNVs whose genotype frequencies deviate from Hardy Weinberg equilibrium. In many cases removed here, it is assumed that the reference sequence itself is incomplete or is an artifact due to a systematic alignment error. This procedure removed 1.03% of the SNV, and finally 77.20% of the original SNV remained after filtering (“Procedure 6” in FIG. 1 and “Procedure 6” in Table 2 above).

Claims (23)

A method for selecting gene mutation information by a computer, wherein the following steps (1) to (3) are executed.
(1) The following steps for each DNA sample:
(1-1) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism is grouped into a subset having a predetermined range of the read depth. Divided,
(1-2) Genes of mutations between the primary mutation set information and the secondary mutation set information obtained from each DNA sample by means other than mapping, based on the measured value of the genotype coincidence of each subset of mutations Calculate by comparing types,
(1-3) determining an appropriate range of read depths by combining the range of read depths of subsets having a genotype match measurement value higher than a predetermined threshold;
(2) For each DNA sample, from the primary mutation set information in step (1-1), removing mutations having a read depth outside the appropriate range of read depths to obtain remaining mutation information for each DNA sample;
(3) As a result of performing the above steps (1) and (2) for at least one DNA sample among all DNA samples, the remaining mutation information is extracted together with the locus to which the mutation information belongs. And selecting as target mutation information.   The method according to claim 1, wherein the gene mutation information includes gene mutation information for single nucleotide mutation (SNV), insertion, or deletion.   The method according to claim 1, wherein the target organism is a human.   The method according to any one of claims 1 to 3, wherein means other than mapping for obtaining secondary mutation set information is a SNP array.   The method according to any one of claims 1 to 4, wherein the measured value of genotype coincidence is precision and / or detection rate (power).   The method according to any one of claims 1 to 5, wherein the predetermined threshold value in the step (1-3) is a matching rate of 0.995 or more.   The method according to any one of claims 1 to 5, wherein the predetermined threshold value in the step (1-3) is a matching rate of 0.995 or more and a detection rate of 0.93 or more. . 1 defined as mutation information to be removed in one or more steps selected from the following (A) to (D) by performing the above steps (1) to (3) to obtain first specific mutation information. The above mutation information is removed from the first specific mutation information to obtain second specific mutation information, and the second specific mutation information is selected as target mutation information. The method in any one of 1-7.
(A) When the first specific mutation information is obtained by using two or more DNA samples, the ratio of the specific mutation information having a read depth outside the appropriate range of the read depth is the step (3). If more than 5% to 20% in the mutation site in the, the specific mutation information is defined as the mutation information to be excluded above;
(B) All the mutation information existing in the area assigned to the low complexity area is selected from the primary mutation set information, and the measured value of the genotype coincidence is the all mutation information and each low complexity area. All mutation information in the low-complexity region that is calculated between the mutation information in the second mutation set information and has a genotype-matched measurement value below a predetermined threshold is defined as the mutation information to be excluded. ;
(C) By using a different mapping algorithm and different mutation call algorithm from the generation of the primary mutation set information used in the first step (1), step (1) or steps (1) to (3) Different primary mutation set information or different mutation information generated is obtained for each DNA sample, and mutation information that is included in the first mutation information but not in the different primary mutation set information or different mutation information is excluded. Defined as power variation information;
(D) An algorithm related to the Hardy-Weinberg equilibrium test is performed, and mutation information deviating from the Hardy-Weinberg equilibrium is defined as mutation information to be removed.   The method according to claim 8, wherein the ratio of mutation information having a lead depth outside the appropriate range of the lead depth exceeds 7% to 15% in the step (A).   The method according to claim 8 or 9, wherein the secondary mutation set information in the step (B) is obtained by using a SNP array.   The method according to any one of claims 8 to 10, wherein the measured value of the genotype coincidence in the step (B) is a precision and / or a detection rate.   12. The method according to claim 11, wherein the predetermined threshold value in the step (B) is a relevance ratio of 0.995 or more.   The method according to claim 11, wherein the predetermined threshold value in the step (B) is a matching rate of 0.995 or higher and a detection rate of 0.93 or higher. A computer system for selecting specific mutation information from base sequence data, comprising at least a recording unit and an arithmetic processing unit, and the following processes (a) to (f):
(A) Based on the mapping algorithm and the mutation call algorithm, the recording unit maps the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and mapping from each DNA sample Secondary mutation set information obtained by means other than is recorded,
(B) In the arithmetic processing unit, a process of grouping primary mutation set information read from the recording unit into a subset having a predetermined range of read depth is executed,
(C) A genotype comparison process of the mutation information between the primary mutation set information and the secondary mutation set information read from the recording unit as the measured value of the genotype matching degree of the mutation information of each subset Is calculated by executing
(D) The determination process of the appropriate range of the read depth is performed by combining the range of the read depth of the subset having the genotype matching degree higher than the predetermined threshold into one,
(E) From the primary mutation set information for each DNA sample, mutation information having a read depth outside the appropriate range of the read depth determined in (d) above is executed, and the remaining mutation information for each DNA sample Is extracted,
(F) As a result of performing the above processing (a) to (e) for at least one DNA sample among all DNA samples, the remaining mutation information is extracted together with the gene locus to which the mutation information belongs. , A process for identifying the genetic information as target specific mutation information;
A computer system characterized in that is executed. The specific mutation information obtained by the execution of the computer system is defined as the first specific mutation information, and is further defined as mutation information to be removed from one or more processes selected from the following (α) to (δ) The above-described mutation information is removed from the first specific mutation information to obtain the second specific mutation information, and the second specific mutation information is extracted as target mutation information. The computer system according to claim 14, wherein:
(Α) When the first specific mutation information is obtained by using two or more DNA samples, the ratio of the specific mutation information having a read depth outside the appropriate range of the read depth is When 5% to 20% of the mutation information obtained by the process (f) is exceeded, the specific mutation information is defined as the mutation information to be removed;
(Β) All the mutation information existing in the area assigned to the low complexity area is extracted from the primary mutation set information, and the measured value of the genotype coincidence is the all mutation information and each low complexity area. All mutation information in the low-complexity region that is calculated and processed with the mutation information in the second mutation set information regarding and has a genotype-matched measurement value below a predetermined threshold is defined as the mutation information to be excluded. Processing;
(Γ) different primary mutation set information or different mutation information for each DNA sample generated by using a different mapping algorithm and different mutation call algorithm from the generation of the first used primary mutation set information. A process in which the mutation information that is included in the mutation information extracted first by the process (f) of 14 but is not included in the different primary mutation set information or the different mutation information is defined as the mutation information to be removed;
(Δ) A process in which an algorithm related to the Hardy-Weinberg equilibrium test is performed, and mutation information deviating from the Hardy-Weinberg equilibrium is defined as mutation information to be removed. A computer program for selecting specific mutation information from base sequence data, and the following functions (a) to (f):
(A) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and obtained by means other than mapping from each DNA sample A function to read out the information from the recording section where the secondary mutation set information is recorded,
(B) a function of performing a process of grouping primary mutation set information read from the recording unit into a subset having a predetermined range of read depths;
(C) A genotype comparison process of the mutation information between the primary mutation set information and the secondary mutation set information read from the recording unit as the measured value of the genotype matching degree of the mutation information of each subset Function to calculate by executing
(D) a function of executing a process of determining an appropriate range of read depth by combining the range of read depths of a subset having a measurement value of genotype matching higher than a predetermined threshold;
(E) From the primary mutation set information for each DNA sample, mutation information having a read depth outside the appropriate range of the read depth determined in (d) above is executed, and the remaining mutation information for each DNA sample The ability to extract,
(F) As a result of executing the above processes (a) to (e) for at least one DNA sample among all DNA samples, the remaining mutation information is extracted together with the locus to which the mutation information belongs. , The function of identifying the genetic information as the specific mutation information of interest;
The computer program characterized by including the algorithm which implement | achieves. In addition to the algorithm for realizing the functions (a) to (f) in the computer program, it is further defined as mutation information to be removed by realizing one or more functions selected from the following (α) to (δ). The one or more mutation information is removed from the existing first specific mutation information, and an algorithm for realizing a function of obtaining new second specific mutation information is included. Computer program.
(Α) When the first specific mutation information is obtained by using two or more DNA samples, the ratio of the specific mutation information having a read depth outside the appropriate range of the read depth is the function (f) A function that defines the specific mutation information as the above-described mutation information to be excluded when exceeding 5% to 20% of the mutation information obtained by
(Β) All the mutation information existing in the area assigned to the low complexity area is extracted from the primary mutation set information, and the genotype match value is measured as the all mutation information and each low complexity area. All mutation information in the low-complexity region having a genotype-matched measurement value below a predetermined threshold is defined as the mutation information to be excluded. function;
(Γ) Different primary mutation set information or different mutation information for each DNA sample generated by using a different mapping algorithm and different mutation call algorithm from the generation of the first used primary mutation set information. a function of defining, as the above-described mutation information to be excluded, the mutation information that is included in the mutation information extracted first by f) but not included in the different primary mutation set information or the different mutation information;
(Δ) A function for defining the mutation information deviating from the Hardy-Weinberg equilibrium as the mutation information to be removed by performing an algorithm related to the Hardy-Weinberg equilibrium test.   A digital program in which the computer program according to claim 16 or 17 is stored in a computer-readable state.   Mutation information obtained by the method according to claim 1. A method for determining an appropriate range of lead depth by a computer, comprising the following steps (1) to (3).
(1) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, the primary mutation set information, and the predetermined range of the read depth Group into subsets with
(2) The gene of the mutation information between the primary mutation set information and the secondary mutation set information obtained from each DNA sample by means other than mapping, using the measured value of the genotype coincidence of the mutation information of each subset. Calculate by comparing types,
(3) The appropriate range of the read depth is determined by combining the read depth ranges of the subsets having the genotype matching degree higher than the predetermined threshold.   21. The method according to claim 20, wherein the measured value of genotype coincidence is a precision and / or detection rate. A computer program for calculating an appropriate range of read depth from base sequence data, and the following functions (a) to (d):
(A) Based on the mapping algorithm and the mutation call algorithm, the primary mutation set information to which the read depth is assigned for each DNA sample of one or more individuals of the target organism, and obtained by means other than mapping from each DNA sample A function to read out the information from the recording section where the secondary mutation set information is recorded,
(B) a function of performing a process of grouping primary mutation set information read from the recording unit into a subset having a predetermined range of read depths;
(C) A genotype comparison process of the mutation information between the primary mutation set information and the secondary mutation set information read from the recording unit as the measured value of the genotype matching degree of the mutation information of each subset Function to calculate by executing
(D) a function of executing a process of determining an appropriate range of read depths by combining the range of read depths of a subset having a genotype matching degree measurement value higher than a predetermined threshold;
The computer program characterized by including the algorithm which implement | achieves.   23. A digital medium in which the computer program according to claim 22 is stored in a computer-readable state.

Images