WO2016068245A1 - 遺伝子型推定装置、方法、及びプログラム - Google Patents

遺伝子型推定装置、方法、及びプログラム Download PDF

Info

Publication number
WO2016068245A1
WO2016068245A1 PCT/JP2015/080573 JP2015080573W WO2016068245A1 WO 2016068245 A1 WO2016068245 A1 WO 2016068245A1 JP 2015080573 W JP2015080573 W JP 2015080573W WO 2016068245 A1 WO2016068245 A1 WO 2016068245A1
Authority
WO
WIPO (PCT)
Prior art keywords
genotype
estimation
data
threshold
snp
Prior art date
Application number
PCT/JP2015/080573
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 JP2016556630A priority Critical patent/JP6280657B2/ja
Publication of WO2016068245A1 publication Critical patent/WO2016068245A1/ja
Priority to US15/419,903 priority patent/US11355219B2/en

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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • 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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • 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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/40—Population genetics; Linkage disequilibrium

Definitions

  • Embodiments of the present invention relate to a genotype estimation apparatus, method, and program.
  • a gene has genetic information that greatly affects the appearance and constitution of the individual with the gene. What holds the information of the gene in the living body is called DNA.
  • the DNA of many organisms has a double helix structure excellent in self-replication for growth and structural strength, and retains genetic information as double-stranded DNA.
  • a single DNA has four types of units called bases and has genetic information as a base sequence. These bases include A (adenine), T (thymine), C (cysin), and G (guanine). The combination of these four patterns represents differences in organisms. In double-stranded DNA, A (adenine) in one place based on one DNA, T (thymine) in another DNA, G (guanine) in C (sicin), etc. In the other DNA, the corresponding bases are determined in a complementary manner. In addition, the number of bases (base-pair) constituting a living body varies depending on the organism, and baboon DNA is composed of about 3 billion base pairs.
  • SNP single nucleotide polymorphism
  • SNP refers to a point where a specific base of a base sequence varies depending on an individual, and is a type of mutation.
  • a combination pattern of a plurality of base sequences (alleles) at a certain SNP locus is called a genotype.
  • the genotype of SNP is related to the individual's appearance and constitution, the degree of risk of a specific disease, how the drug works, and how the side effects of the drug appear.
  • GWAS Gene-Wide Association Study
  • DNA microarray technology is a genotyping technology that can determine genotypes of tens of thousands to hundreds of thousands of SNPs (tag SNPs) that almost cover the whole genome information of a plurality of people at a high throughput at a time. .
  • a known base sequence around a certain SNP is used as a probe, an unknown base sequence of a sample is hybridized to measure the fluorescence intensity, the fluorescence intensity of a plurality of specimens is mapped to a cluster space, and clustering is performed. Determine the genotype. This utilizes the fact that each SNP is a cluster in which two types of fluorescence intensity values reflecting genotype information are grouped for each genotype in space.
  • the clustering strength is also calculated as the clustering reliability.
  • the clustering strength is high when the genotypes of the specimen group are well separated and there is a grouping of each genotype. However, it is difficult to accurately assign genotypes to specimens with fluorescence intensities farther from the cluster, so do not assign genotypes to unreliable specimens by setting a threshold for cluster intensity. There is also.
  • genotype estimation apparatus capable of accurately estimating a genotype.
  • the genotype estimation apparatus includes an acquisition unit and an estimation unit.
  • the acquisition unit acquires clustering strength of genotype data of a plurality of samples including an unknown sample whose genotype is unknown and a known sample whose genotype is known.
  • the estimation unit estimates the genotype of the unknown specimen based on the genotype data and outputs an estimation result.
  • the figure which shows an example of the operation screen of GUI of the genotype estimation apparatus of FIG. The figure which shows the other example of the operation screen of GUI of the genotype estimation apparatus of FIG.
  • FIG. 1 is a block diagram illustrating a functional configuration of the estimation apparatus according to the present embodiment.
  • the estimation apparatus includes a specimen data storage unit 1, a reference data storage unit 2, a clustering strength acquisition unit 3, a clustering strength determination unit 4, a genotype estimation unit 5, and an estimation result display. Part 6.
  • the specimen data storage unit 1 stores data (specimen data) relating to a specimen whose genotype has been determined by the DNA microarray technology.
  • the sample data includes, for example, genotype data, signal intensity data, and clustering data, but is not limited thereto.
  • Genotype data is data indicating the result of genotype determination by DNA microarray technology.
  • the genotype data includes the genotype determined for each specimen and for each SNP.
  • FIG. 2 is a diagram showing an example of genotype data.
  • the genotype data in FIG. 2 includes the genotype determination results of SNPrs000001 to rs9999999 of specimens 01 to N.
  • the genotype of SNPrs000001 of specimen 01 is “CG”.
  • “-” indicates that the genotype could not be determined by the DNA microarray technique.
  • the estimation apparatus estimates such a genotype that could not be determined by the DNA microarray technology.
  • a specimen whose genotype is known in a certain SNP that is, a specimen whose genotype is determined by the DNA microarray technology is referred to as a known specimen.
  • a specimen whose genotype is unknown that is, a specimen whose genotype could not be determined by the DNA microarray technology is referred to as an unknown specimen.
  • the sample 01 is a known sample and the sample 02 is an unknown sample.
  • the signal strength data is data indicating the measurement result of the signal strength by the DNA microarray technology.
  • the signal intensity is, for example, measured values such as fluorescence intensity, current, and voltage, but is not limited thereto. Further, the signal intensity may be an arbitrary parameter calculated from the above measured value.
  • the signal intensity data includes the value of each signal intensity at each SNP of each specimen.
  • FIG. 3 is a diagram showing an example of signal strength data.
  • the signal strength data in FIG. 3 includes the values of the signal strengths x1 to xn of SNPrs000001 to rs9999999 of the specimens 01 to N.
  • the value of the signal intensity x1 of SNPrs000001 of the specimen 01 is 0.8.
  • Clustering data is data indicating the result of clustering performed for each SNP when determining the genotype by DNA microarray technology.
  • the clustering by the DNA microarray technique may be hierarchical clustering or non-hierarchical clustering. In the following, it is assumed that the clustering by the DNA microarray technology is non-hierarchical clustering.
  • the clustering data includes, for example, cluster coordinate data and clustering strength data, but is not limited thereto.
  • each cluster corresponds to each genotype in the SNP, the number of genotypes is generated. For example, when there are three SNP genotypes, CC, CT, and TT, three clusters are generated on the cluster space.
  • the cluster coordinates are, for example, the coordinates of the center of gravity of the cluster.
  • FIG. 4 is a diagram showing an example of cluster coordinate data.
  • the cluster coordinate data in FIG. 4 includes the coordinates of clusters 1 to m of SNPrs000001 to rs999999.
  • the coordinates of each cluster are indicated by two axes v1 and v2 in the cluster space.
  • the cluster coordinate data may include not only the coordinates of each cluster, but also the coordinates of each specimen in the cluster space.
  • vn xn may be sufficient.
  • the cluster space is an n-dimensional space with signal intensities x1 to xn.
  • the clustering strength data is data indicating the clustering strength for each SNP.
  • the clustering strength is an index indicating the reliability of clustering. The greater the clustering strength, the higher the clustering reliability.
  • FIG. 5 is a diagram showing an example of clustering strength data.
  • the clustering strength data in FIG. 5 includes clustering strengths of SNPrs000001 to rs9999999.
  • the clustering strength of SNPrs000001 is 0.95.
  • the clustering strength for example, an average value of distances between the clusters can be used.
  • the clustering strength can be obtained from the cluster coordinate data by the following formula.
  • the inter-cluster distance ij is the Euclidean distance between any two clusters i and j
  • (vi1, vi2) are the coordinates of the center of gravity of cluster i
  • (vj1, vj2) are the coordinates of the center of gravity of cluster j. is there.
  • m is the number of clusters.
  • FIG. 6 is a diagram showing the clustering strength data generated from the cluster coordinate data of FIG.
  • the inter-cluster distances i and j and the clustering strength are standardized so that the clustering strength is a value of 0 or more and 1 or less.
  • the reference data storage unit 2 stores known data (reference data) related to genes.
  • Reference data includes, but is not limited to, linkage disequilibrium statistical data, reference haplotype data, and reference genotype frequency data, for example.
  • reference data large-scale project data such as the international HapMap project or the 1000-person genome project can be used.
  • Linkage disequilibrium statistical data (hereinafter referred to as “LD data”) is data indicating the correlation between SNPs.
  • FIG. 7 is a diagram illustrating an example of LD data.
  • the LD data in FIG. 7 includes SNPrs125678 and SNPrs129688 scores, SNPrs125678 and SNPrs986675 scores, and SNPrs129688 and SNPrs986675 scores.
  • the score is an index indicating the strength of correlation between SNPs.
  • the LD data in FIG. 7 includes a linkage disequilibrium score (D ′), a correlation coefficient (r 2 ), and a logarithm of odds ratio (LOD) as scores.
  • D ′ linkage disequilibrium score
  • r 2 correlation coefficient
  • LOD logarithm of odds ratio
  • Reference haplotype data is data indicating any combination of alleles (bases) of SNPs statistically related on the same chromosome. That is, each reference haplotype data indicates a highly probable base combination in some SNPs.
  • the SNP included in the reference haplotype data is selected based on, for example, LD data.
  • FIG. 8 is a diagram showing an example of reference haplotype data.
  • the reference haplotype data in FIG. 8 includes reference haplotype data refHTD1 to refHTD6.
  • Each reference haplotype data includes alleles of SNPrs123456, rs623456, rs987456, rs987123, rs598456, rs387456, rs912346, rs778456, rs873456, rs98709.
  • the allele of SNPrs123456 in the reference haplotype data refHTD1 is A.
  • Reference genotype frequency data is data indicating the genotype frequency (reference genotype frequency) of each SNP in a certain population.
  • FIG. 9 is a diagram showing an example of reference genotype frequency data.
  • the reference genotype frequency data in FIG. 9 includes the frequencies of SNPrs125678 genotypes CC, CT, and TT.
  • the allele of SNPrs125678 is C or T.
  • the total frequency of each genotype is 1.
  • the frequency of the genotype CC of SNPrs125678 is 0.42.
  • Clustering strength acquisition unit 3 acquires the clustering strength for each SNP.
  • the acquisition unit 3 acquires the clustering intensity data from the sample data storage unit 1.
  • the acquisition unit 3 acquires the cluster coordinate data from the sample data storage unit 1, and calculates the clustering strength of each SNP based on the cluster coordinate data. May be.
  • the calculation method of the clustering strength is as described above.
  • the acquisition unit 3 obtains genotype data from the sample data storage unit 1.
  • the reference genotype frequency data may be acquired from the reference data storage unit 2, and the clustering strength of each SNP may be calculated based on the genotype data and the reference genotype frequency data.
  • the calculation method of the clustering strength is as follows.
  • the acquisition unit 3 calculates the frequency of each genotype (DNA microarray genotype frequency) of each SNP based on the genotype data.
  • the DNA microarray genotype frequency is the frequency of the genotype determined by the DNA microarray technology.
  • the acquisition unit 3 calculates the clustering strength of each SNP by the following formula based on the DNA microarray genotype frequency and the reference genotype frequency.
  • m is the number of genotypes
  • f i is the reference genotype frequency of genotype i
  • f i, D are the DNA microarray genotype frequencies of genotype i. Since the genotype frequency indicates the probability for each genotype , the sum of fi , r and the sum of fi, r are both 1.
  • the genotype of a certain SNP is CC, CT, TT
  • the reference genotype frequencies are 0.5, 0.3, 0.2, respectively
  • the DNA microarray genotype frequencies are 0.4, 0.4, respectively.
  • the clustering strength determination unit 4 obtains the clustering strength for each SNP acquired by the acquisition unit 3, the threshold value ⁇ 1 (first threshold value), and the threshold value ⁇ 2 (second threshold value). Compare.
  • the threshold values ⁇ 1, ⁇ 2 ( ⁇ 1 ⁇ ⁇ 2) are values set in advance to determine the reliability of clustering by the DNA microarray technology.
  • the determination unit 4 determines that the clustering reliability is high when the clustering strength is greater than the threshold ⁇ 1, and determines that the reliability is low when the clustering strength is less than the threshold ⁇ 2, and when the clustering strength is greater than or equal to the threshold ⁇ 2 and less than or equal to the threshold ⁇ 1.
  • the clustering reliability is determined to be medium.
  • the threshold values ⁇ 1 and ⁇ 2 depend on the clustering strength, and when the clustering strength is a value in the range of 0 to 1, the value is in the range of 0 to 1. For example, when the clustering strength is 0 or more and 1 or less, the threshold values ⁇ 1 and ⁇ 2 are set to 0.8 and 0.4, respectively.
  • the determination unit 4 determines that the reliability is high when the clustering strength is greater than the threshold ⁇ 1, and determines that the reliability is low when the clustering strength is equal to or less than the threshold ⁇ 1.
  • the genotype estimation unit 5 (hereinafter referred to as “estimation unit 5”) estimates the genotype of an unknown sample in each SNP of the genotype data. For example, the estimation unit 5 estimates the genotype of SNPrs000003 of the sample 01 and the genotype of SNPrs000002 of the sample 02 in the genotype data of FIG.
  • the estimation unit 5 selects an estimation method based on the determination result of the clustering reliability by the determination unit 4. For example, when the clustering strength is larger than the threshold ⁇ 1, that is, when the reliability of clustering by the DNA microarray technology is high, the estimation unit 5 estimates the genotype by the k-nearest neighbor method based on the genotype data. Further, when the clustering strength is smaller than the threshold value ⁇ 2, that is, when the reliability of clustering by the DNA microarray technology is low, the estimation unit 5 estimates the genotype by the imputation method based on the genotype data and the reference data.
  • the estimation unit 5 estimates the genotype by using the k-nearest neighbor method and the imputation method together when the clustering strength is not less than the threshold ⁇ 2 and not more than the threshold ⁇ 1, that is, the reliability of clustering by the DNA microarray technology is medium. . And the estimation part 5 outputs an estimation result. Specific examples of the genotype estimation method will be described later in detail.
  • the estimation result display unit 6 (hereinafter referred to as “display unit 6”) displays the estimation result by the estimation unit 5.
  • the display unit 6 may display genotype data and various types of information used in the estimation together with the estimation result.
  • the estimation apparatus includes a computer 100.
  • the computer 100 includes a CPU (central processing unit) 101, an input device 102, a display device 103, a communication device 104, and a storage device 105, which are connected to each other via a bus 106.
  • CPU central processing unit
  • the CPU 101 is a control device and a calculation device of the computer 100.
  • the CPU 101 performs arithmetic processing based on data or a program input from each device (for example, the input device 102, the communication device 104, and the storage device 105) connected via the bus 106, and outputs the calculation result and the control signal.
  • the data is output to each device (for example, the display device 103, the communication device 104, and the storage device 105) connected via the bus 106.
  • the CPU 101 executes an OS (operating system) of the computer 100, a genotype estimation program (hereinafter referred to as “estimation program”), and the like, and controls each device constituting the computer 100.
  • the estimation program is a program that causes the computer 100 to realize the above-described functional configurations of the estimation device. When the CPU 101 executes the estimation program, the computer 100 functions as an estimation device.
  • the input device 102 is a device for inputting information to the computer 100.
  • the input device 102 is, for example, a keyboard, a mouse, and a touch panel, but is not limited thereto.
  • the user can input information such as threshold values ⁇ ⁇ b> 1 and ⁇ ⁇ b> 2 by using the input device 102.
  • the display device 103 is a device for displaying an image, a video, or the like based on data output from the CPU 101 or the like.
  • the display device 103 is, for example, an LCD (liquid crystal display), a CRT (CRT), and a PDP (plasma display), but is not limited thereto.
  • the display unit 6 can be configured using the display device 103.
  • the communication device 104 is a device for the computer 100 to communicate with an external device wirelessly or by wire.
  • the communication device 104 is, for example, a modem, a hub, and a router, but is not limited thereto.
  • Information such as sample data and reference data can be input by being received from an external device via the communication device 104.
  • data such as calculation results output from the CPU 101 can be transmitted to an external device.
  • the storage device 105 is a storage medium that stores the OS of the computer 100, an estimation program, data necessary for execution of the estimation program, data generated and output by execution of the estimation program by the CPU 101, and the like.
  • the storage device 105 includes a main storage device and an external storage device.
  • the main storage device is, for example, a RAM, a DRAM, or an SRAM, but is not limited thereto.
  • the external storage device is a hard disk, an optical disk, a flash memory, and a magnetic tape, but is not limited thereto.
  • the sample data storage unit 1 and the reference data storage unit 2 can be configured using the storage device 105.
  • the computer 100 may include one or more CPUs 101, input devices 102, display devices 103, communication devices 104, and storage devices 105, or may be connected to peripheral devices such as printers and scanners. .
  • the estimation device may be configured by a single computer 100 or may be configured as a system including a plurality of computers 100 connected to each other.
  • the estimation program may be stored in advance in the storage device 105 of the computer 100, may be stored in a storage medium such as a CD-ROM, or may be uploaded on the Internet.
  • the estimation apparatus can be configured by installing the estimation program in the computer 100 and executing it.
  • FIG. 11 is a flowchart showing an outline of the operation of the estimation apparatus according to the present embodiment.
  • step S1 the acquisition unit 3 acquires genotype data from the specimen data storage unit 1.
  • step S2 the acquisition unit 3 extracts the genotype data of the target SNP from the genotype data.
  • the target SNP is a SNP including an unknown sample.
  • the acquisition unit 3 extracts genotype data of SNPrs000002 and rs000003.
  • step S3 the acquisition unit 3 acquires the clustering strength CS of each target SNP.
  • the acquisition unit 3 can acquire the clustering strength CS based on genotype data, clustering data, reference genotype frequency data, and the like.
  • step S4 the determination unit 4 acquires the clustering strength CS of each target SNP from the acquisition unit 3 and compares it with the threshold value ⁇ 1.
  • the threshold value ⁇ 1 may be the same for each target SNP or may be different.
  • step S4 the determination unit 4 determines that the clustering reliability is high (YES in step S4), and the process proceeds to step S5.
  • step S5 the estimation unit 5 estimates the genotype of the unknown specimen by the k-nearest neighbor method based on the genotype data.
  • a genotype estimation method by the k-nearest neighbor method will be described later.
  • step S4 If CS ⁇ ⁇ 1 (NO in step S4), the process proceeds to step S6.
  • step S6 the determination unit 4 compares the clustering strength CS of each target SNP acquired from the acquisition unit 3 with the threshold value ⁇ 2.
  • the threshold value ⁇ 2 may be the same for each target SNP or may be different.
  • step S6 the determination unit 4 determines that the clustering reliability is low (YES in step S6), and the process proceeds to step S7.
  • step S7 the estimation unit 5 estimates the genotype of the unknown specimen by the imputation method based on the genotype data and the reference data.
  • the genotype estimation method by the imputation method will be described later.
  • step S6 the determination unit 4 determines that the clustering reliability is medium, and the process proceeds to step S8.
  • step S8 the estimation unit 5 estimates the genotype of the unknown specimen by using both the k-nearest neighbor method and the imputation method.
  • a genotype estimation method using the k-nearest neighbor method and the imputation method will be described later.
  • the display unit 6 displays the estimation result by the estimation unit 5 in step S9.
  • the k-nearest neighbor method is a method for estimating the genotype of an unknown specimen based on the genotypes of k nearest samples.
  • a method using a known specimen and a cluster line as samples will be described.
  • FIG. 12 is a flowchart showing an example of a genotype estimation method by the k-nearest neighbor method.
  • a known specimen is used as a sample.
  • step S511 the estimation unit 5 acquires genotype data and signal intensity data of the unknown sample S of the target SNP from the sample data storage unit 1. In the following, it is assumed that there is one unknown sample S. However, when there are a plurality of unknown samples S, the following processing is performed for each unknown sample S.
  • step S512 the estimation unit 5 acquires the genotype data and signal intensity data of the known sample group ST of the target SNP from the sample data storage unit 1.
  • the known sample group ST is a set of known samples Si included in the target SNP.
  • the estimation unit 5 calculates a distance di for each known sample Si included in the known sample group ST.
  • the distance di is a distance between the unknown sample S and the known sample Si. For example, when the signal intensity data of the unknown specimen S is (x1, x2,..., Xn) and the signal intensity data of the known specimen Si is (xi1, xi2,. Calculated by the formula.
  • step S514 the estimation unit 5 selects the k known samples Si nearest to the known sample group ST, that is, k known samples Si in ascending order of the distance di.
  • the parameter k is an arbitrary natural number set in advance. A method for setting the parameter k will be described later.
  • FIG. 13 is a diagram for explaining a method of selecting a known specimen Si.
  • FIG. 13 shows an example of the clustering map of the target SNP.
  • a genotype CC a genotype.
  • a circle is a known specimen with a genotype CC
  • a triangle is a genotype.
  • three known specimens with a genotype of CC and two known specimens with a genotype of CG are selected in ascending order of the distance di.
  • step S515 the estimation unit 5 estimates the genotype of the unknown sample S based on the selected k known sample Si genotypes.
  • the estimation unit 5 estimates the genotype of the unknown specimen S using, for example, a majority algorithm. That is, the estimation unit 5 estimates the genotype having the largest number of samples (number of votes) among the k known genotypes of the known sample Si as the genotype of the unknown sample S.
  • FIG. 14 is a diagram for explaining a genotype estimation method using the majority algorithm.
  • the number of votes for each of the genotypes AG, GG, and AA is 3, 1, and 1, respectively, the genotype of the unknown sample S is estimated to be AG having the largest number of votes.
  • the estimation unit 5 may estimate the genotype of the unknown specimen S using a weighted majority algorithm. In this case, the estimation unit 5 first calculates the weight of each selected known specimen Si. As the weight of the known specimen Si, the ratio of SNPs whose genotypes are determined in the known specimen Si can be used. For example, the weight of the known specimen Si in which 150,000 genotypes are determined out of 200,000 SNPs by the DNA microarray technology is 0.75.
  • the estimation unit 5 uses the weight of each known specimen Si as the number of votes, and estimates the genotype with the largest number of votes as the genotype of the unknown specimen S.
  • FIG. 15 is a diagram for explaining a genotype estimation method using a weighted majority algorithm.
  • the genotype of the unknown sample S is estimated to be the AG with the largest number of votes.
  • FIG. 16 is a flowchart showing another example of the genotype estimation method by the k-nearest neighbor method.
  • cluster lines are used as samples.
  • step S521 the estimation unit 5 acquires the genotype data and signal intensity data of the unknown sample S of the target SNP from the sample data storage unit 1.
  • Step S521 is the same as step S511 described above.
  • step S522 the estimation unit 5 acquires the genotype data and signal intensity data of the known sample group ST of the target SNP from the sample data storage unit 1.
  • Step S522 is the same as step S512 described above.
  • the estimation unit 5 creates a cluster line Ci based on the signal intensity data of the known sample group ST.
  • the cluster line Ci is obtained by approximating a known sample included in each cluster (each genotype) on the cluster space by a straight line or a curve.
  • the cluster line Ci can be created by performing regression analysis on the coordinates of a known specimen in the cluster space. The regression analysis may be linear regression or non-linear regression.
  • FIG. 17 is a diagram for explaining a method of creating the cluster line Ci.
  • FIG. 17 shows an example of the clustering map of the target SNP.
  • the cluster line Ci is a straight line and is created one for each cluster.
  • the cluster line Ci is represented by the following equation.
  • m i and c i are constants and are obtained by regression analysis.
  • the estimation unit 5 can create the cluster line Ci as described above by performing linear regression analysis for each cluster.
  • a plurality of cluster lines Ci may be created for each cluster.
  • the estimation unit 5 may divide each cluster into a plurality of sub-clusters and create a cluster line Ci for each sub-cluster as shown in Expression (7).
  • cluster line Ci is not limited to a straight line, but may be a curved line. Further, the number of cluster lines Ci can be arbitrarily selected.
  • step S524 the estimation unit 5 calculates a distance Di for each cluster line Ci.
  • the distance Di is a distance between the unknown sample S and the cluster line Ci.
  • step S525 the estimation unit 5 selects the k cluster lines Ci nearest to each other, that is, the k cluster lines Ci in ascending order of the distance Di from the plurality of created cluster lines Ci.
  • step S526 the estimation unit 5 estimates the genotype of the unknown sample S based on the genotypes of the selected k cluster lines Ci.
  • the estimation unit 5 estimates the genotype of the unknown specimen S using, for example, a majority algorithm. That is, the estimation unit 5 estimates the genotype having the largest number of cluster lines (number of votes) as the genotype of the unknown specimen S among the g types of the selected k cluster lines Ci.
  • FIG. 19 is a diagram for explaining a genotype estimation method using a majority vote algorithm.
  • the number of votes for each of the genotypes AG, GG, and AA is 3, 1, and 1, respectively, the genotype of the unknown sample S is estimated to be AG having the largest number of votes.
  • the estimation unit 5 may estimate the genotype of the unknown specimen S using a weighted majority algorithm. In this case, the estimation unit 5 first calculates the weight of each selected cluster line Ci. As the weight of the cluster line Ci, an average value for each cluster of the ratio of SNPs whose genotypes are determined in the known specimen Si can be used. The estimation unit 5 estimates the genotype with the largest number of votes as the genotype of the unknown sample S using the weight of each cluster line Ci as the number of votes.
  • FIG. 20 is a flowchart showing a method for setting the parameter k.
  • the estimation unit 5 sets the parameter k by cross verification.
  • step S531 the estimation unit 5 acquires genotype data and signal strength data of one or more verification SNPs from the specimen data storage unit 1.
  • the verification SNP is a SNP having a large clustering strength CS and all samples are known samples.
  • the verification SNP is, for example, an SNP in which the clustering strength CS is larger than the threshold value ⁇ 1.
  • FIG. 21 is a diagram for explaining the verification SNP.
  • the samples of SNPrs00001 and rs000003 are all known samples.
  • the estimation unit 5 extracts SNPrs00001 and rs000003 as verification SNPs, and acquires these genotype data and signal strength data.
  • step S532 the estimation unit 5 selects an evaluation sample and a learning sample.
  • An evaluation sample is a sample handled as an unknown sample.
  • the learning sample is a sample handled as a known sample.
  • the genotype of the sample selected as the evaluation sample is used as correct data for cross-validation.
  • FIG. 22 is a diagram illustrating an example of an evaluation sample and a learning sample.
  • samples 01 to 10 are selected as samples for evaluation, and samples 11 to N are selected as samples for learning. Note that the sample for evaluation and the sample for learning can be arbitrarily selected.
  • step S533 the estimation unit 5 sets a plurality of parameter k candidates k ′.
  • the estimation unit 5 can set an arbitrary natural number as the parameter k candidate k ′.
  • step S534 the estimation unit 5 estimates the genotype of each evaluation sample based on the genotype data and signal intensity data of the learning sample. At this time, the estimation unit 5 estimates the genotype of the sample for evaluation by the k-nearest neighbor method using each candidate k ′ as the parameter k.
  • step S535 the estimation unit 5 calculates the estimation accuracy of each candidate k ′ by cross verification. That is, the estimation unit 5 compares the estimation result of the genotype of the evaluation sample with the known genotype of the evaluation sample, and calculates the ratio of the evaluation sample with the genotype correctly estimated.
  • FIG. 23 is a diagram illustrating an example of a calculation result of the estimation accuracy.
  • the estimation accuracy is calculated for each candidate k ′ of each verification SNP.
  • the estimation unit 5 may calculate an average value (average estimation accuracy) of the estimation accuracy of each candidate k ′.
  • step S536 the estimation unit 5 sets the candidate k ′ having the maximum estimation accuracy as the parameter k. Further, when the cross-validation is performed for a plurality of verification SNPs, the estimation unit 5 may set a candidate k ′ having the maximum average estimation accuracy as the parameter k. For example, in the case of FIG. 23, the parameter k is set to 5 at which the average estimation accuracy is maximized.
  • FIG. 24 is a flowchart showing a genotype estimation method by the imputation method.
  • step S71 the estimation unit 5 refers to the reference data storage unit 2 and confirms whether there is LD data of the target SNP in the reference data storage unit 2. Since the imputation method uses the LD data of the target SNP, if there is no LD data of the target SNP (NO in step S71), the process proceeds to step S5, and the genotype of the unknown sample S is estimated by the k-nearest neighbor method. To do.
  • the estimation method by the k-nearest neighbor method is as described above.
  • step S71 when there is LD data of the target SNP (YES in step S71), the process proceeds to step S72.
  • step S72 the estimation unit 5 acquires the LD data of the target SNP from the reference data storage unit 2.
  • step S73 the estimation unit 5 refers to the LD data of the target SNP and extracts L or more SNPs having a high score.
  • FIG. 25 is a flowchart showing the SNP extraction method in step S73.
  • step S731 the estimation unit 5 extracts the LD data of the estimation SNP from the LD data of the target SNP.
  • the estimation SNP is a SNP having a large clustering strength CS and all samples are known samples.
  • the estimation SNP is, for example, an SNP in which the clustering strength CS is larger than the threshold value ⁇ 1.
  • the LD data of the target SNP includes scores of other SNPs with respect to the target SNP.
  • the estimation unit 5 refers to the genotype data and clustering data of each SNP whose score is included in the LD data of the target SNP, and extracts the LD data of the estimation SNP.
  • step S732 the estimation unit 5 sets the score threshold LD ⁇ to the highest score included in the LD data of the estimation SNP.
  • the score can be arbitrarily selected from among linkage disequilibrium score, correlation coefficient, logarithm of odds ratio, and the like included in the LD data.
  • the estimation unit 5 refers to the LD data of the estimation SNP and extracts an SNP having a score equal to or higher than the threshold LD ⁇ from the estimation SNP.
  • step S734 the estimation unit 5 determines whether the number of extracted SNPs is equal to or greater than a predetermined value L.
  • L can be set arbitrarily. If the number of extracted SNPs is smaller than L (NO in step S734), the process proceeds to step S735.
  • the reduction amount ⁇ of the threshold value LD ⁇ is, for example, 0.01.
  • the estimation unit 5 repeats the processing of steps S733 to S735 until the number of extracted SNPs becomes L or more.
  • the estimation unit 5 can extract L or more estimation SNPs having a score equal to or greater than the threshold value LD ⁇ .
  • step S74 the estimation unit 5 acquires the reference haplotype data of the SNP extracted in step S73 from the reference data storage unit 2.
  • step S75 the estimation unit 5 acquires the genotype data of the unknown sample S from the sample data storage unit 1.
  • the estimation unit 5 creates haplotype data of the unknown sample S from the genotype data of the unknown sample S.
  • the haplotype data of the unknown specimen S can be created by extracting the SNP genotype from the genotype data using a fading (phasing) algorithm and determining the sequence of alleles present in each chromosome.
  • a fading algorithm for example, BEAGLE, fastPHASE, IMPUTEv2, MACH, and ShapeIT can be used.
  • FIG. 26 is a diagram showing an example of haplotype data of the unknown sample S.
  • two haplotype data HTD1 and HTD2 are created from the genotype data of the unknown specimen S by the fading algorithm.
  • the allelic pair of the SNP whose genotype is unknown is unknown, so it is indicated by “ ⁇ ”.
  • step S77 the estimation unit 5 selects reference haplotype data most similar to the two haplotype data HTD1 and HTD2 of the unknown sample S from the reference haplotype data acquired in step S74.
  • the similarity between the haplotype data HTD1 and HTD2 and the reference haplotype data means that the sequences of alleles in SNPs other than the SNP whose genotype is unknown are similar.
  • FIG. 27 is a diagram illustrating a method for selecting reference haplotype data.
  • the estimation unit 5 selects the reference haplotype data refHTD5 as the reference haplotype data most similar to the haplotype data HTD1, and is most similar to the haplotype data HTD2.
  • Reference haplotype data refHTD3 is selected as reference haplotype data. Details of the method for selecting the reference haplotype data will be described later.
  • step S78 the estimation unit 5 estimates the genotype of the unknown specimen S based on the alleles of the target SNP in the two selected reference haplotype data. For example, when reference haplotype data is selected as shown in FIG. 27, the estimation unit 5 estimates one allele of SNPrs987709 as G and estimates the other allele as A. And the estimation part 5 estimates the genotype of SNPrs987007 of the unknown sample S based on these alleles as AG.
  • FIG. 28 is a flowchart showing a method for selecting reference haplotype data in step S77.
  • step S771 the estimation unit 5 replaces the alleles of the extracted reference haplotype data and haplotype data of the unknown sample S with numerical values.
  • FIG. 29 is a diagram illustrating an example of reference haplotype data in which alleles are replaced with numerical values and haplotype data of an unknown sample S.
  • alleles A, C, G, and T are replaced with numerical values 1, 2, 3, and 4, respectively.
  • step S772 the estimation unit 5 calculates the distance dhi.
  • the distance dhi is a distance between each haplotype data of the unknown specimen S and each reference haplotype data.
  • the distance dhi is calculated by the following equation, for example.
  • step S773 the estimation unit 5 selects the reference haplotype data having the smallest distance dhi as the most similar reference haplotype data for each haplotype data of the unknown sample S.
  • FIG. 30 is a flowchart showing a genotype estimation method using both the k-nearest neighbor method and the imputation method.
  • step S81 the estimation unit 5 estimates the genotype of the unknown specimen S by the k-nearest neighbor method, and obtains a genotype group GT1 composed of one or a plurality of genotype candidates.
  • the genotype group GT1 selects, for example, ⁇ genotypes in descending order of votes as genotype candidates, or ⁇ It can be obtained by estimating the genotype using the parameter k.
  • the estimation unit 5 estimates the genotype of the unknown specimen S by the imputation method, and obtains a genotype group GT2 composed of one or a plurality of genotype candidates.
  • the genotype group GT2 selects, for example, ⁇ reference haplotype data that is similar for each haplotype data of the unknown sample S, and Can be obtained by estimating.
  • step S83 the estimation unit 5 estimates the genotype of the unknown specimen S from the genotype candidates included in the genotype groups GT1 and GT2 using a majority algorithm. As the number of votes of the majority algorithm, the number of genotypes included in the genotype groups GT1 and GT2 can be used.
  • the estimation apparatus and method according to the present embodiment estimate a genotype that could not be determined by the DNA microarray technology by a method according to the reliability of clustering by the DNA microarray technology. That is, when the reliability is low, the estimation is performed by the imputation method using the reference data, and when the reliability is high, the estimation is performed by the k-nearest neighbor method using the genotype data determined by the DNA microarray technology. . Thereby, the estimation apparatus and method which concern on this embodiment can estimate a genotype accurately.
  • FIG. 31 is a diagram illustrating an example of a GUI operation screen displayed by the display unit 6.
  • FIG. 31 shows a GUI when the genotype is estimated by the k-nearest neighbor method using a known specimen as a sample.
  • the GUI includes an SNP selection unit G1, a sample selection unit G2, a k value selection unit G3, a sample list display unit G4, a selection result display unit G5, and a genotype display unit G6. .
  • the SNP selection unit G1 is a drop-down list for the user to select a target SNP.
  • the drop-down list of the SNP selection unit G1 includes IDs of all SNPs including unknown samples.
  • the SNP selection unit G1 displays the ID of the target SNP selected by the user.
  • the specimen selection unit G2 is a drop-down list for the user to select an unknown specimen S whose genotype is estimated.
  • the drop-down list of the sample selection unit G2 includes IDs of all unknown samples S included in the genotype data of the target SNP.
  • the content of the drop-down list of the sample selection unit G2 changes according to the target SNP selected by the user.
  • the ID of the unknown sample S selected by the user is displayed.
  • the k value selection unit G3 is a drop-down list for the user to set the parameter k.
  • the drop-down list of the k value selection unit G3 includes a plurality of parameter k value candidates.
  • the k value selection unit G3 displays the value of the parameter k set by the user. In FIG. 31, the parameter k is set to 5. In the k value selection unit G3, it is preferable that the value of the parameter k having the highest estimation accuracy is set as a default value.
  • the sample list display unit G4 is a list of IDs of known samples (known samples Si included in the known sample group ST) included in the genotype data of the target SNP selected by the user, each known sample Si, and the user's selection.
  • the distance di between the unknown sample S and the unknown sample S is displayed.
  • the distance di displayed in FIG. 31 is a distance calculated by the above equation (6).
  • the selection result display unit G5 displays the ID, distance di, genotype, and weight of k known samples Si having a small distance di selected from the known samples Si displayed on the sample list display unit G4. To do. In FIG. 31, since the parameter k is 5, five known specimens Si are displayed.
  • the weight is a weight used in the majority algorithm, and 1.0 is set as a default value. When the weight is 1.0, a majority algorithm without weight is obtained. When the weighted algorithm is used, the weight of each known specimen Si calculated by the above method is displayed as the weight.
  • the genotype display unit G6 displays the genotype estimation result by the k-nearest neighbor method.
  • the estimated genotype is AG.
  • FIG. 32 is a diagram illustrating another example of a GUI operation screen of the estimation apparatus.
  • the estimation apparatus provided with the GUI of FIG. 32 stores information indicating the relationship between the SNP and the disease in the reference data storage unit 2.
  • the GUI includes a patient selection unit G7, a disease selection unit G8, a SNP information display unit G9, and a genotype information display unit G10.
  • the patient selection unit G7 is a drop-down list for the user to select a patient.
  • the patient here corresponds to a specimen whose genotype has been determined by the DNA microarray technology.
  • the drop-down list of the patient selection unit G7 includes IDs of a plurality of patients (specimens).
  • the patient selection unit G7 displays the ID of the patient selected by the user.
  • the disease selection unit G8 is a drop-down list for the user to select a disease.
  • the drop-down list of the disease selection unit G8 includes names of a plurality of diseases stored in the reference data storage unit 2.
  • the name of the disease selected by the user is displayed.
  • the SNP information display unit G9 displays SNP information related to the patient and the disease selected by the user.
  • SNP information includes the type of chromosome, SNP ID, locus, disease-related SNP, odds ratio (OR), Addr information, base sequence information, and the like.
  • the odds ratio here is a scale used as a method of showing the results of medical clinical trials, and is a statistical scale showing the susceptibility to disease in two groups. Addr information and base sequence information will be described later.
  • the SNP information is stored in the reference data storage unit 2.
  • the SNP information display part G9 displays the 1st to 5th, 10th to 12th, and XY chromosomes.
  • a hatched portion on each chromosome indicates a SNP whose genotype is known among the disease-related SNPs, and a dot portion indicates a SNP whose genotype is unknown among the disease-related SNPs.
  • the SNP on each chromosome is a command button, and when the user selects (clicks), Addr information on the SNP and base sequences around the SNP are displayed.
  • the Addr information includes the number of the chromosome (Chromosome) to which the SNP belongs, the locus on the chromosome (Position), the name of the gene to which the SNP belongs (Gene), and the ID of the SNP.
  • the Addr information of the SNP selected by the user with the command button is displayed in the Addr column of the SNP information display unit G9.
  • Base sequence information is base sequence data of a gene locus excluding SNP.
  • the base sequence around the selected SNP is extracted from the base sequence information, and the genotype data of the SNP included in the extracted base sequence range is extracted from the specimen data storage unit 1 And displayed in the base sequence column of the SNP information display part G9.
  • alleles A and B of SNP are displayed in the format [A / B].
  • the alleles of SNPrs547794 are G and T.
  • the genotype information display unit G10 displays genotype information related to the SNP selected by the user.
  • the genotype information is generated from various data stored in the specimen data storage unit 1.
  • the genotype information display unit G10 When the genotype of the SNP selected by the user is known, the genotype information display unit G10, as shown in FIG. 32, the clustering map (Genotype Clustering) of the selected SNP, the genotype ratio by HapMap, etc. Is displayed.
  • the clustering map Geneotype Clustering
  • the clustering map may be stored in advance in the sample data storage unit 1 or may be generated from signal intensity data stored in the sample data storage unit 1 or the like.
  • the genotype ratio by HapMap is the ratio of the genotype of the selected SNP in the ethnic group to which the patient belongs.
  • the genotype ratio by HapMap can be extracted from the reference genotype frequency data.
  • the genotype information display unit G10 estimates the genotype estimation result or genotype by the estimation unit 5 as shown in FIG. Displays data indicating the process to be performed.
  • Imputation displayed on the genotype information display unit G10 is a display label, which indicates a genotype estimation method by the estimation unit 5.
  • the genotype display unit G10 displays the reference haplotype data, the haplotype data of the patient (specimen), the haplotype data after fading, and the post-imputation data. Display haplotype data.
  • the haplotype data after imputation includes the estimated genotype.
  • the estimated genotype (allele) is displayed in the base sequence column of the SNP information display part G9. Further, the genotype information display unit G10 may display the genotype ratio by HapMap as in FIG.
  • the genotype information display unit G10 displays the sample list display unit G4, the selection result display unit G5, the genotype display unit G6, etc. in FIG.
  • the displayed information may be displayed.
  • FIG. 34 to 58 A second embodiment will be described with reference to FIGS. 34 to 58.
  • FIG. 34 to 58 a genotype estimation method using a threshold method will be described.
  • the threshold method is a genotype estimation method in which a correspondence relationship between a signal intensity interval and a genotype is learned, and a genotype of each specimen is estimated based on the learned correspondence relationship.
  • Each section of the signal strength is defined by a signal strength threshold. Details of the threshold method will be described later.
  • FIG. 34 is a block diagram showing a functional configuration of the estimation apparatus according to this embodiment.
  • the estimation apparatus includes a sample data storage unit 1, an estimation unit 5, and a display unit 6.
  • the first embodiment differences from the first embodiment will be described.
  • the specimen data storage unit 1 stores genotype data and signal intensity data as specimen data, and does not store clustering data. Further, the estimation device does not include the reference data storage unit 2, the acquisition unit 3, and the determination unit 4.
  • the threshold method does not use clustering data, reference data, and clustering strength.
  • clustering data is stored in the specimen data storage unit 1 and reference data is stored in the estimation device.
  • the unit 2, the acquisition unit 3, and the determination unit 4 may be provided.
  • the estimation unit 5 includes a threshold learning unit 51 and a threshold method estimation unit 52.
  • the threshold learning unit 51 learns the correspondence between the signal strength interval used in the threshold method and the genotype based on the signal strength of the Fullcall SNP. Specifically, the learning unit 51 learns the correspondence relationship between the threshold value that defines the signal intensity interval and the genotype.
  • the full call SNP is an SNP in which the genotypes of all the samples are determined by the DNA microarray technology, that is, an SNP in which all the samples are known samples.
  • an SNP in which at least one specimen genotype has not been determined by the DNA microarray technology, that is, an SNP including at least one unknown specimen is referred to as a Nocall SNP.
  • FIG. 35 is a diagram illustrating an example of genotype data stored in the specimen data storage unit 1.
  • the genotypes of all specimens are determined for SNPrs00001 and rs999999. Therefore, SNPrs00001, rs999999 are full call SNPs.
  • SNPrs000002 and rs000003 specimens 02 and 01 are unknown specimens, respectively. Therefore, SNPrs000002 and rs000003 are no call SNPs.
  • the learning unit 51 can grasp the full call SNP and the no call SNP by referring to the genotype data.
  • the learning unit 51 extracts the full call SNP from the genotype data and the signal intensity of each sample of the full call SNP from the signal intensity data in order to learn the threshold value.
  • the learning unit 51 may extract any one type of signal strength to be learned.
  • the type of signal intensity with which the learning unit 51 learns the threshold value can be arbitrarily set by the user of the estimation device.
  • the learning unit 51 extracts the signal strength x1 and learns the threshold value of the signal strength x1 will be described as an example.
  • the number of thresholds learned by the learning unit 51 can be arbitrarily set by the user of the estimation device, and may be one or plural.
  • the number of thresholds is preferably set according to the type of genotype included in each SNP.
  • the number of signal strength sections is specified by one more than the threshold number. Therefore, when the maximum number of types of genotypes included in each SNP is X, the learning unit 51 may learn, for example, X ⁇ 1 threshold values.
  • a threshold value x 1 first threshold value
  • a threshold value x r second threshold value
  • the threshold method estimation unit 52 determines the genotype of each sample of the no-call SNP based on the correspondence between the signal intensity interval (threshold) learned by the learning unit 51 and the genotype. presume.
  • the no-call SNP includes an unknown sample and a known sample. Therefore, in this embodiment, not only estimation of the genotype of an unknown sample but also estimation (re-determination) of the genotype of a known sample is performed.
  • SNPrs000002 estimates the genotype of the sample 02 which is an unknown sample, and also estimates the genotypes of the samples 01 and N which are known samples.
  • the genotype estimation method using the threshold method will be described in detail later.
  • the hardware configuration of the estimation apparatus according to the present embodiment is the same as that of the first embodiment. That is, when the computer 100 executes the estimation program, each functional configuration described above of the estimation device is realized.
  • FIG. 36 is a flowchart showing an outline of a threshold learning method. Details of each step will be described later.
  • the learning unit 51 generates a threshold combination list.
  • the threshold combination list is a list including a plurality of threshold combinations.
  • the threshold combination is a combination of threshold candidates.
  • the threshold combination is a combination of a threshold x 1 candidate and a threshold x r candidate.
  • step S11 the learning unit 51 calculates a genotype frequency for evaluating each threshold combination included in the threshold combination list.
  • step S12 the learning unit 51 calculates an evaluation value for each threshold combination based on the threshold candidates included in each threshold combination and the genotype frequency.
  • step S13 the learning unit 51 selects the threshold combination having the maximum evaluation value from the threshold combinations included in the threshold combination list.
  • Each threshold candidate included in the selected threshold combination is employed as a threshold for estimating the genotype by the threshold method.
  • FIG. 37 is a flowchart illustrating an example of a method for generating a threshold combination list in step S10.
  • step S101 the learning unit 51 acquires the genotype data of all SNPs and the signal strength data of the signal strength x1 of all SNPs from the specimen data storage unit 1.
  • FIG. 38 is a diagram showing an example of acquired genotype data and signal intensity data.
  • genotype data and signal intensity data of SNPrs000001 to rs9999999 are acquired.
  • step S102 the learning unit 51 replaces each genotype included in the genotype data acquired in step S101 with a cluster number.
  • the cluster number is a value assigned according to the relative position of each cluster on the clustering map.
  • the learning unit 51 first assigns a cluster number to each cluster of each SNP.
  • FIG. 39 is a diagram illustrating an example of a cluster number assignment method.
  • cluster numbers 0, 1, and 2 are assigned in order from the cluster located on the right side of the clustering map. This corresponds to assigning cluster numbers 0, 1, and 2 in order from the cluster having the highest signal intensity x1 at the center of gravity.
  • the learning unit 51 can calculate the signal intensity x1 of the center of gravity of each cluster from the genotype data and the signal intensity data acquired in step S101, and can assign a cluster number.
  • cluster coordinate data is stored in the sample data storage unit 1
  • the learning unit 51 acquires cluster coordinate data from the sample data storage unit 1, and refers to the acquired cluster coordinate data to obtain a cluster number. May be assigned.
  • Cluster number is assigned to each SNP by a common method. Therefore, in some SNPs, cluster number 0 may be assigned to a cluster of genotype AA, and in other SNPs, cluster number 0 may be assigned to a cluster of genotype CC.
  • the learning unit 51 assigns a cluster number to each cluster, and then replaces the genotype of each specimen with the cluster number assigned to the cluster containing the specimen. For example, when the cluster number 0 is assigned to a certain cluster, the learning unit 51 replaces the genotype of each sample included in the cluster with 0.
  • FIG. 40 is a diagram showing an example of genotype data after the genotype is replaced with the cluster number.
  • the genotype data in FIG. 40 corresponds to the genotype data in FIG.
  • 0, 1, and 2 are cluster numbers corresponding to genotypes, respectively, and -1 is a cluster number corresponding to a genotype not being determined.
  • the genotype CG is replaced with cluster number 1, and the genotype CC is replaced with cluster number 2.
  • the genotype AT is replaced with cluster number 1
  • the genotype TT is replaced with cluster number 2. This indicates that the relative position of the genotype CG cluster in SNPrs000001 is equal to the relative position of the genotype AT in SNPrs000002.
  • the cluster numbers are assigned to each cluster in descending order of the signal strength x1, but may be assigned in ascending order of the signal strength x1, or may be assigned in descending or ascending order of the signal strength x2. May be.
  • genotype of each specimen is expressed using a cluster number.
  • the learning unit 51 refers to the genotype data after replacement and extracts a full call SNP.
  • SNPrs000001, rs999998, and rs999999 are extracted as full call SNPs.
  • step S104 the learning unit 51 extracts the signal strength data of the full call SNP extracted in step S103 from the signal strength data acquired in step S101.
  • FIG. 41 is a diagram illustrating an example of the extracted signal strength data of the full call SNP.
  • the signal strength data of FIG. 41 is obtained by extracting the signal strength data of the full call SNP from the signal strength data of FIG.
  • step S105 the learning unit 51 divides the signal strength data of the full call SNP extracted in step S104 into signal strength data of one cluster SNP and signal strength data of a plurality of cluster SNPs.
  • 1-cluster SNP is an SNP that contains only one type of genotype as a determination result. That is, one cluster SNP is a SNP in which all specimens are determined to have the same genotype.
  • a multi-cluster SNP is a SNP that includes multiple types of genotypes as a determination result. As a determination result, all SNPs including two or more genotypes are included in the multiple cluster SNP. However, the above-mentioned genotype-1 is not included in the two or more genotypes mentioned here.
  • the learning unit 51 counts the types of genotypes included in the genotype data of each SNP in order to divide the signal strength data.
  • the type of genotype included in the genotype data of a certain SNP is one type (for example, genotype 1)
  • the learning unit 51 determines that the SNP is one cluster SNP.
  • the learning unit 51 determines that the SNP is a multi-cluster SNP.
  • the learning unit 51 divides the signal strength data of the full call SNP based on the SNP determination result thus obtained.
  • FIG. 42 is a diagram illustrating an example of signal strength data of one class SNP and signal strength data of a plurality of cluster SNPs.
  • the signal strength data in FIG. 42 is obtained by dividing the signal strength data in FIG.
  • SNPrs9999998 is a one-cluster SNP that includes only one type of genotype 0, and SNPrs000001 and rs999999 are multi-cluster SNPs that include three types of genotypes 0,1,2. Therefore, as shown in FIG. 42, the signal strength data of one cluster SNP includes the signal strength data of SNPrs999998, and the signal strength data of the plurality of cluster SNPs includes the signal strength data of SNPrs000001 and 999999. Yes.
  • the signal strength data of the one cluster SNP and the plurality of cluster SNPs are divided because the difference between the cluster distribution in the one cluster SNP and the cluster distribution in the plurality of cluster SNPs is large.
  • the accuracy of genotype estimation can be improved by performing signal intensity threshold learning and genotype estimation using the learned threshold for each of the one cluster SNP and the plurality of cluster SNPs.
  • the estimation device can also process one cluster SNP and a plurality of cluster SNPs together.
  • the learning unit 51 does not divide the signal strength data in step S105, and may process the one cluster SNP and the plurality of clusters SNP collectively for the subsequent processing.
  • the learning unit 51 refers to the signal strength data of one cluster SNP divided in step S105, and calculates a signal strength statistical value for each genotype. Further, the learning unit 51 refers to the signal strength data of the plurality of cluster SNPs divided in step S105, and calculates a signal strength statistical value for each genotype.
  • Statistic value of signal strength includes minimum value, average value, maximum value, and standard deviation value.
  • the learning unit 51 refers to the genotype data after replacement, extracts the signal intensity of genotype 0 from the signal intensity data of one cluster SNP (or multiple clusters SNP) after division, and the statistical value of the extracted signal intensity , The statistical value of the signal intensity of genotype 0 of one cluster SNP (or multiple clusters SNP) can be calculated. Other genotypes are calculated in the same manner.
  • FIG. 43 is a diagram illustrating an example of a statistical value of signal strength of one cluster SNP and a statistical value of signal strength of a plurality of cluster SNPs.
  • the minimum value of the signal intensity of the genotype 2 of one cluster SNP is ⁇ 6.29.
  • step S107 the learning unit 51 associates the SNP signal intensity interval of one cluster SNP with the genotype based on the average value of the signal intensity of each genotype of one cluster SNP. Further, the learning unit 51 associates the signal intensity intervals of the multiple clusters SPN with the genotypes based on the average value of the signal intensity of each genotype of the multiple clusters SNP.
  • the learning unit 51 learns two signal intensities xl and xr, three signal intensity sections are formed.
  • the learning unit 51 associates genotypes having a small average value of signal strength in order from a section having a low signal strength.
  • FIG. 44 is a diagram illustrating an example of a correspondence relationship between signal intensity intervals and genotypes.
  • the correspondence relationship of FIG. 44 is based on the average value of the signal intensity of FIG. In the example of FIG. 43, the average value of the signal intensity is smaller in the order of genotype 2, 1, 0. Therefore, in the example of FIG. 44, genotypes 2, 1, and 0 are associated in order from the section with the lowest signal intensity.
  • the signal strength is genotype 2 is associated with the period of less than x l
  • the signal strength is genotype 1 is associated with the following section x l or x r
  • the signal strength is x r greater than interval Is associated with genotype 0.
  • step S108 the learning unit 51 generates a threshold candidate list of one cluster SNP based on the statistical value of one cluster SNP calculated in step S106.
  • the learning unit 51 generates a threshold candidate list for the multi-cluster SNP based on the statistical value of the multi-cluster SNP calculated in step S106.
  • the threshold candidate list is a list including a plurality of threshold candidates.
  • the threshold candidate is a candidate of signal intensity thresholds x 1 and x r .
  • the threshold candidates for one cluster SNP are, for example, the minimum value, the average value, the maximum value, the average value + N ⁇ standard deviation of the signal intensity of each genotype of one cluster SNP (or a plurality of cluster SNPs).
  • N is an integer), but is not limited thereto.
  • FIG. 45 is a diagram illustrating an example of a threshold candidate list for one cluster SNP and a threshold candidate list for a plurality of cluster SNPs.
  • the threshold candidate list in FIG. 45 corresponds to the statistical values in FIG. 43 and includes nine threshold candidates (minimum value, average value, maximum value of each genotype).
  • the threshold candidate ⁇ 6.29 included in the one-cluster SNP threshold candidate list is the minimum value of the signal intensity of the genotype 2 of the one-cluster SNP.
  • each threshold candidate x i is calculated by the following equation.
  • n is the number of threshold candidates included in the threshold candidate list
  • x min is the minimum statistical value
  • x max is the maximum statistical value
  • d is the threshold candidate interval.
  • the minimum value x min of the statistical value of one cluster SNP corresponds to ⁇ 6.29 which is the minimum value of the signal strength of genotype 2
  • the maximum value x max is the signal strength of genotype 0.
  • the maximum value of 7.46 corresponds to the maximum value of 7.46.
  • threshold candidates included in the threshold candidate list are not limited to those described above, and can be generated by an arbitrary method from the statistical values of signal strength. Further, the threshold candidate list may include an arbitrary value set in advance as a threshold candidate.
  • step S109 the learning unit 51 refers to the threshold candidate list for one cluster SNP generated in step S108 and generates a threshold combination list for one cluster SNP. Further, the learning unit 51 refers to the threshold candidate list for the multiple cluster SNP generated in step S108 and generates a threshold combination list for the multiple cluster SNP.
  • the threshold combination list is a list including a plurality of threshold combinations.
  • the learning unit 51 generates a threshold combination by combining the threshold candidates included in the threshold candidate list, and generates a threshold combination list including a plurality of threshold combinations.
  • n threshold candidates are included in the threshold candidate list and r threshold candidates are included in the threshold combination
  • n! / (Nr)! r! Threshold combinations are generated. Therefore, when nine threshold candidates are included in the threshold candidate list of one cluster SNP and two threshold candidates x 1 and x r are included in the threshold combination, a maximum of 36 threshold combinations are generated.
  • FIG. 46 is a diagram illustrating an example of a threshold combination list of one cluster SNP.
  • the threshold combination list in FIG. 46 corresponds to the threshold candidate list in FIG.
  • the threshold combination list in FIG. 46 includes 36 threshold combinations.
  • a threshold combination list as shown in FIG. 46 is also generated for a plurality of cluster SNPs.
  • the learning unit 51 calculates a genotype frequency of one cluster SNP and a genotype frequency of a plurality of cluster SNPs.
  • the genotype frequency is the number of samples having a signal intensity included in a signal intensity interval defined by two adjacent threshold candidates. The genotype frequency is calculated for each genotype.
  • the learning unit 51 refers to the signal intensity data divided in step S105, the genotype data after replacement, and the threshold candidate list generated in step S108, so that each section of the one-cluster SNP and the plurality of cluster SNPs The genotype frequency of each genotype can be calculated.
  • the signal strength of the sample 01 of SNPrs9999998 is 0.3. It can be seen that this signal strength is included in a section between the threshold candidate 0.69 and the threshold candidate 2.11. Then, referring to the genotype data of FIG. 40, the genotype of the sample 01 of SNPrs999998 is 2. As a result, the genotype frequency of genotype 2 in the section between 0.69 and 2.11.
  • the learning unit 51 refers to each signal intensity included in the signal intensity data of one cluster SNP, and adds the genotype frequencies of each genotype in each section as described above, and the genotype frequency of one cluster SNP. Is calculated. The genotype frequency is also calculated for the multiple cluster SNPs by the same method.
  • FIG. 47 is a diagram showing an example of the genotype frequency of one cluster SNP.
  • the genotype frequency in FIG. 47 corresponds to the threshold candidate list in FIG. In the example of FIG. 47, the genotype frequency of genotype 2 is high, the genotype frequency of genotype 1 is low, and the genotype 0 is low in the interval between threshold candidate-2.11 and threshold candidate-1.79. Genotype frequency is zero.
  • the learning unit 51 calculates an evaluation value of each threshold combination of one cluster SNP, and selects a threshold combination of one cluster SNP based on the calculated evaluation value. Further, the learning unit 51 calculates an evaluation value of each threshold combination of the multiple cluster SNP, and selects a threshold combination of the multiple cluster SNP based on the calculated evaluation value.
  • the learning unit 51 re-determines the genotype of one cluster SNP in the full call SNP based on each threshold combination.
  • two re-determination methods will be described.
  • the learning unit 51 selects one threshold combination from the threshold combination list of one cluster SNP, and extracts the signal intensity and genotype of the specimen of one cluster SNP. Then, the learning unit 51, the signal strength is again determined genotype of x l is less than the sample and genotype 2, the signal strength is the genotype of the following sample x l or x r re determined to genotype 1, signal strength is again determined that the genotype 0 genotype x r larger specimens.
  • FIG. 48 is a diagram illustrating an example of a clustering map before and after redetermination by the first redetermination method.
  • a circle indicates each specimen, and a value in the circle indicates a genotype value.
  • the upper diagram in FIG. 48 shows the genotype before redetermination, and the lower diagram shows the genotype after redetermination.
  • the genotype of the sample whose signal intensity is included in the section between the threshold candidate xl and the threshold candidate x r is re-determined as all genotypes 1.
  • step S11 can be omitted when the first redetermination method is employed.
  • the genotype frequency is used in the second re-determination method.
  • the signal strength genotype following sample x l or x r is re-determined by majority algorithm using genotype frequencies. That is, a specimen having a signal intensity of xl or more and xr or less is re-determined to have a genotype having the highest genotype frequency in the section including the signal intensity of the specimen.
  • genotype 2 which has the highest genotype frequency in the interval between ⁇ 2.11 and ⁇ 1.79.
  • FIG. 49 is a diagram illustrating an example of a clustering map before and after redetermination by the second redetermination method. As shown in FIG. 49, in the second re-determination method, the specimen whose signal intensity is included in the section between the threshold candidate xl and the threshold candidate x r is not necessarily re-determined as genotype 1. .
  • the learning unit 51 re-determines the genotype of one cluster SNP by the re-determination method as described above, and then calculates a match rate between the genotype before re-determination and the genotype after re-determination. For example, in the example of FIG. 48, the genotypes of 12 specimens out of 16 specimens match before and after redetermination, so the coincidence rate is 0.75. In the example of FIG. 49, since the genotypes of 13 samples out of 16 samples match before and after redetermination, the matching rate is 0.81.
  • the coincidence rate calculated in this way is used as the evaluation value of the threshold combination used for redetermination.
  • the learning unit 51 calculates the evaluation value of each threshold combination included in the threshold combination list of one cluster SNP by the above method.
  • the learning unit 51 also calculates an evaluation value of each threshold combination included in the threshold combination list of the multiple cluster SNP by the same method.
  • FIG. 50 is a diagram illustrating an example of a threshold combination list including evaluation values.
  • the evaluation value of threshold combination 1 is 0.80
  • the evaluation value of threshold combination 24 is 0.97.
  • the learning unit 51 selects the threshold combination having the maximum evaluation value from the threshold combinations included in the threshold combination list of one cluster SNP. Each threshold candidate included in the selected threshold combination is learned as a threshold for estimating the genotype of one cluster SNP by the threshold method.
  • the learning unit 51 selects the threshold combination having the maximum evaluation value from the threshold combinations included in the threshold combination list of the multiple cluster SNP.
  • Each threshold candidate included in the selected threshold combination is learned as a threshold for estimating the genotypes of the multiple cluster SNPs by the threshold method.
  • the threshold combination with the maximum evaluation value in the threshold combination list is the threshold combination 24, and thus the threshold candidates included in the threshold combination 24 are learned as thresholds for estimating the genotype. Is done. That is, when the genotype is estimated by the threshold method, ⁇ 0.80 is used as the threshold x 1 and 2.11 is used as the threshold x r .
  • the threshold value By learning the threshold value as described above, it is possible to improve the accuracy of genotype estimation by the threshold method. This is because the higher the evaluation value (matching rate) of the threshold combination is, the closer the genotype estimation result by the threshold combination and the genotype estimation result in the full call SNP are. In general, since the genotype estimation accuracy in the full call SNP is high, the estimation accuracy becomes higher as the threshold combination can obtain an estimation result closer to the estimation result in the full call SNP.
  • FIG. 51 is a flowchart showing an outline of a genotype estimation method by the threshold method.
  • step S14 the estimation unit 52 acquires genotype data of all SNPs and signal strength data of signal strengths x1 of all SNPs from the specimen data storage unit 1.
  • step S15 the estimation unit 52 replaces each genotype included in the genotype data acquired in step S14 with a cluster number.
  • the genotype replacement method is as described in step S102.
  • the estimation unit 52 refers to the genotype data after replacement and extracts a no-call SNP.
  • SNPrs000002 and rs000003 are extracted as no-call SNPs.
  • the estimation part 52 selects object SNP from the no call SNP extracted by step S16.
  • the target SNP is a SNP that is a target whose genotype is estimated by the threshold method.
  • threshold learning is performed using all full-call SNPs, but genotype estimation is performed for each individual no-call SNP.
  • the method for selecting the target SNP is arbitrary.
  • step S18 the estimation unit 52 estimates the genotype of each specimen of the target SNP selected in step SS17 by the threshold method. Step S18 will be described in detail later.
  • step S19 If there is an unselected no-call SNP as the target SNP (YES in step S19), the estimation unit 52 selects the next target SNP from the unselected no-call SNPs (step S17). Thereafter, steps S17 to S19 are repeated until there is no unselected no-call SNP.
  • the estimation unit 52 ends the genotype estimation process.
  • the estimation unit 52 estimates the genotypes of all samples of the target SNP using only the signal intensity thresholds x 1 and x r .
  • FIG. 52 is a flowchart showing the first estimation method.
  • step S1801 the estimation unit 52 extracts the genotype data and signal strength data of the target SNP from the genotype data and signal strength data after replacement.
  • step S1802 the estimation unit 52 refers to the genotype data of the target SNP and acquires the number of clusters of the target SNP.
  • the method for obtaining the number of clusters is as described above. That is, the estimation unit 52 acquires the number of clusters of the target SNP by counting the types of genotypes included in the genotype data of the target SNP.
  • step S1802 may be omitted.
  • the estimation part 52 should just collect the 1 cluster SNP and multiple cluster SNP, and may perform the subsequent processes.
  • step S ⁇ b> 1803 the estimation unit 52 acquires the threshold values x l and x r of the signal strength according to the number of clusters of the target SNP from the learning unit 51.
  • Estimating unit 52 when the number of clusters target SNP is 1, the threshold x l of 1 cluster SNP, acquires x r, if the number of clusters target SNP is 2 or more, the threshold x l multiple clusters SNP, x r To get.
  • the estimation unit 52 estimates the genotype of each sample of the target SNP based on the threshold values x 1 and x r acquired in step S1803. Specifically, the estimation unit 52 estimates the genotype of the specimen whose signal intensity x i is smaller than x l (x i ⁇ x l ) as genotype 2, and the signal intensity is x l or more and x r or less (x l ⁇ x i ⁇ x r) the genotype of the sample was estimated to genotype 1, the signal strength is estimated that x r larger than (x r ⁇ x i) genotype the genotype of the sample 0.
  • FIG. 53 is a diagram showing an example of a clustering map before and after estimation by the first estimation method.
  • the upper diagram of FIG. 53 shows the genotype before estimation by the first estimation method, and the lower diagram shows the genotype after estimation by the first estimation method.
  • the genotypes of all the samples including the unknown sample and the known sample are estimated based on the threshold values x 1 and x r .
  • the estimation unit 52 ends the genotype estimation process of the target SNP.
  • the estimation unit 52 selects the next target SNP (step S17).
  • the estimation unit 52 estimates the genotype of the known sample of the target SNP using the signal intensity thresholds x l and x r, and then uses the k-nearest neighbor method to determine the genotype of the unknown sample.
  • FIG. 54 is a flowchart showing the second estimation method. Steps S1801 to S1803 in FIG. 54 are the same as those in the first estimation method, and thus description thereof is omitted.
  • step S1805 based on the threshold values x 1 and x r acquired in step S1803, the estimation unit 52 generates a genotype of a sample having the maximum signal intensity and a sample having the minimum signal intensity among the samples of the target SNP. Is estimated. Specifically, the estimation unit 52, the signal strength is a genotype x l is less than the sample estimates that genotype 2, the signal strength is estimated that the original genotype The genotype of the following sample x l or x r and the signal strength of the genotype x r larger sample estimates genotype 0.
  • FIG. 55 is a diagram showing an example of the clustering map before and after the estimation in step S1805.
  • the upper diagram in FIG. 55 shows the genotype before estimation, and the lower diagram shows the genotype after estimation.
  • the genotype of the sample with the lowest signal intensity (leftmost sample) is genotype 2
  • the genotype of the sample with the highest signal intensity (rightmost sample) is the gene.
  • Type 1 The signal intensities of these two analytes are all smaller than the threshold value x l.
  • the genotypes of the two specimens are both estimated to be genotype 2.
  • the genotype of the right end of the specimen is estimated genotype 1 (original genotype)
  • the signal strength of the right end of the sample of Figure 55 is less than x l or x r
  • step S1806 the estimation unit 52 determines whether the genotypes of the two samples estimated in step S1805 are the same. If the genotypes of the two specimens are the same (YES in step S1806), the process proceeds to step S1807.
  • step S1807 the estimation unit 52 estimates the genotypes of all samples as the same genotypes as the genotypes of the two samples estimated in step S1805. This is because when the genotypes of the specimens having the maximum and minimum signal intensity are both estimated to be the same genotype X, the genotypes of all the specimens are estimated to be the same genotype X in the subsequent processing. .
  • steps S1805 to S1807 only the genotypes of two specimens are estimated first, and the genotypes of all specimens are estimated based on the estimation results, thereby reducing the calculation amount of the estimation apparatus.
  • FIG. 56 is a diagram showing an example of the clustering map before and after the estimation in step S1807.
  • the upper diagram in FIG. 56 shows the genotype before estimation, and the lower diagram shows the genotype after estimation.
  • the upper diagram in FIG. 56 corresponds to the lower diagram in FIG.
  • the genotypes of the left end and right end specimens are both estimated to be genotype 2
  • the genotypes of all specimens are estimated to be genotype 2 in step S1807.
  • step S1807 not only the known specimen but also the genotype of the unknown specimen is estimated, so the k-nearest neighbor method is not used.
  • the estimation unit 52 ends the genotype estimation process of the target SNP.
  • the estimation unit 52 selects the next target SNP (step S17).
  • step S1805 determines whether the genotypes of the two samples estimated in step S1805 are different (NO in step S1806). If the genotypes of the two samples estimated in step S1805 are different (NO in step S1806), the process proceeds to step S1808.
  • step S1808 the estimation unit 52 estimates the genotypes of all samples of the target SNP based on the threshold values x 1 and x r acquired in step S1803.
  • the estimation method is the same as that in step S1805. That is, the estimation unit 52, the signal strength estimates the genotype x l is less than the sample and genotype 2, the signal strength estimates the original genotype The genotype of the following sample x l or x r, signal strength genotypes x r larger sample estimates genotype 0.
  • FIG. 57 is a diagram showing an example of the clustering map before and after the estimation in step S1808.
  • the upper diagram in FIG. 57 shows the genotype before estimation, and the lower diagram shows the genotype after estimation.
  • step S1809 the estimation unit 52 calculates a match rate between the genotype before estimation and the genotype after estimation in step S1808, and determines whether the calculated match rate is greater than a match rate threshold value.
  • the matching rate threshold can be arbitrarily set. In the example of FIG. 57, since the genotypes of 12 samples out of 16 samples match before and after the estimation, the matching rate is 0.75.
  • step S1810 If the match rate is greater than the match rate threshold (YES in step S1809), the process proceeds to step S1810.
  • step S1810 the estimation unit 52 determines whether there is an unknown sample among the samples of the target SNP.
  • the unknown sample corresponds to a sample having a genotype of -1. Therefore, the estimation unit 52 can determine whether there is an unknown sample by referring to the estimation result in step S1808 and confirming whether there is a genotype-1 sample.
  • the estimation unit 52 ends the genotype estimation process of the target SNP.
  • the estimation unit 52 selects the next target SNP (step S17).
  • step S1810 when there is an unknown sample in the sample of the target SNP (YES in step S1810), the process proceeds to step S1811.
  • step S1811 the estimation unit 52 extracts an unknown sample from the sample of the target SNP, and estimates the genotype of each unknown sample by the k-nearest neighbor method.
  • the estimation result in step S1808 is used as the genotype of the known specimen.
  • the genotype estimation method by the k-nearest neighbor method is as described in the first embodiment.
  • the estimation unit 52 extracts, for example, k specimens (or cluster lines) closest to the unknown specimen (the Euclidean distance on the clustering map is closest), and out of the genotypes of the extracted specimens (or cluster lines) The most common genotype may be estimated as the genotype of the unknown specimen.
  • the estimation unit 52 ends the genotype estimation process of the target SNP.
  • the estimation unit 52 selects the next target SNP (step S17).
  • step S1812 the match rate is equal to or less than the match rate threshold (NO in step S1809), the process proceeds to step S1812.
  • step S1812 the estimation unit 52 estimates a genotype by a majority method for each cluster.
  • the estimation unit 52 classifies the known samples of the target SNP for each genotype before estimation in step S1811, and generates a cluster corresponding to each genotype.
  • Each cluster includes a specimen having a genotype corresponding to the cluster.
  • the estimation unit 52 refers to the estimation result in step S1808, and estimates the most common genotype among the genotypes after estimation of the samples included in each cluster as the genotype of the cluster. Then, the estimation unit 52 estimates the genotypes of all the samples included in each cluster as the genotypes of the clusters including the sample.
  • FIG. 58 is a diagram showing an example of the clustering map before and after the estimation in steps S1808 and S1812.
  • the upper diagram in FIG. 58 shows the genotype before estimation in step S1808, the middle diagram shows the genotype after estimation in step S1808, and the lower diagram shows the genotype after estimation in step S1812.
  • the upper diagram and the middle diagram in FIG. 58 correspond to FIG.
  • the estimation unit 52 first generates a cluster 1 including four genotype 1 specimens and a cluster 2 including 11 genotype 2 specimens (see the above figure).
  • the estimation unit 52 refers to the estimation result in step S1808, and estimates the most common genotype among the genotypes after estimation of the samples included in the cluster 1 as the genotype of the cluster 1 (see the middle diagram). ).
  • cluster 1 includes three genotype 0 samples and one genotype 1 sample. Therefore, the estimation unit 52 estimates the genotype of cluster 1 as genotype 1.
  • the estimation unit 52 estimates the genotypes of all the samples included in the cluster 1 as the genotype 0 that is the genotype of the cluster 1 (see the following diagram).
  • the estimation unit 52 also estimates the genotype of each specimen included in the cluster 2 by the same method. As a result, the genotypes of all the specimens included in cluster 2 are estimated as genotype 2.
  • step S1810 The processing after step S1810 is as described above, and in step S1811, the genotype of the unknown sample is estimated by the k-nearest neighbor method.
  • the estimation apparatus estimates the genotype of a no-call SNP by a threshold method using a signal intensity threshold. Since the threshold is learned using the signal intensity data of the full call SNP whose genotype has been determined with high accuracy, the estimation device can accurately estimate the genotype of the no call SNP.
  • the estimation method according to the present embodiment can be used even when the reference data is not sufficiently obtained.
  • the value of the parameter k used in the k-nearest neighbor method may be optimally set by cross verification.
  • the method for setting the parameter k by cross verification is as described in the first embodiment.
  • the estimation apparatus according to the present embodiment may be able to execute the estimation method according to the first embodiment.
  • the estimation method can be selected by the user through the GUI.
  • the estimation device may execute the estimation method selected by the user.
  • FIG. 59 to 61 A third embodiment will be described with reference to FIGS. 59 to 61.
  • FIG. 61 a modified example of the genotype estimation method based on the k-nearest neighbor method used in the estimation methods according to the first embodiment and the second embodiment will be described.
  • FIG. 59 is a diagram illustrating an example of a clustering map in which a plurality of unknown specimens are gathered at relatively close positions.
  • specimens s 1 to s 3 are unknown specimens and are gathered at positions close to each other.
  • the three samples closest to the sample s 1 are the samples s2 and s3 with the genotype -1, and the sample with the genotype 0.
  • One is selected.
  • the genotype of the specimen s 1 is estimated to be ⁇ 1 by the majority algorithm. That is, the genotype of the sample s 1 can not be estimated. For the same reason, the genotypes of the specimens s 2 and s 3 cannot be estimated.
  • FIG. 60 is a flowchart showing the k-nearest neighbor method according to this embodiment.
  • step S20 the estimation unit 5 generates an unknown sample list.
  • the unknown sample list is a list including all unknown samples of the target SNP.
  • step S21 the estimation unit 5 determines whether the unknown sample list generated in step S20 is empty. If the unknown sample list is empty (YES in step S21), that is, if there is no unknown sample in the target SNP, the estimation unit 5 ends the genotype estimation process of the unknown sample by the k-nearest neighbor method.
  • step S21 if the unknown sample list is not empty (NO in step S21), that is, if there is an unknown sample in the target SNP, the process proceeds to step S22.
  • step S22 the estimation unit 5 estimates the genotype of each unknown sample included in the unknown sample list by the k-nearest neighbor method using all samples.
  • the genotype estimation by the k-nearest neighbor method executed in step S22 is as described in the first embodiment.
  • step S23 the estimation unit 5 deletes the unknown sample whose genotype was estimated in step S22 from the unknown sample list.
  • the unknown sample deleted here is treated as a known sample in the subsequent processing.
  • step S24 When one or more unknown samples are deleted from the unknown sample list in step S23 (YES in step S24), the process returns to step S21. This is because the genotype of an unknown sample that could not be genotyped in step S22 may be estimated when one or more unknown samples become new known samples.
  • steps S21 to S24 are repeated until one or more unknown specimens are not deleted from the unknown specimen list (until the genotype of one or more unknown specimens is no longer estimated in step S22).
  • step S25 the process proceeds to step S25.
  • the genotype of an unknown sample cannot be estimated by the k-nearest neighbor method using all samples. That is, as described above, it means that a plurality of unknown specimens are gathered at relatively close positions.
  • step S25 the estimation unit 5 selects a target sample from unknown samples included in the unknown sample list.
  • the target specimen here is an unknown specimen that is a genotype estimation target.
  • the estimation unit 5 may select the target specimen at random or may select it by the following method.
  • the estimation unit 5 calculates the average distance between each unknown sample included in the unknown sample list and k nearest known samples. Then, the estimation unit 5 selects an unknown sample having the smallest average distance from k known samples as a target sample.
  • step S26 the estimation unit 5 estimates the genotype of the target sample selected in step S25 by the k-nearest neighbor method using only known samples. As a result, k known samples are selected as the k samples closest to the target sample. Therefore, the genotype of the target sample can be estimated based on the genotype of the selected known sample.
  • step S27 the estimation unit 5 deletes the target sample from the unknown sample list. Thereafter, the process returns to step S21. Thereafter, the processes of steps S21 to S27 are repeated until the unknown sample list becomes empty.
  • FIG. 61 is a diagram showing a clustering map corresponding to a process in which the genotype of an unknown specimen is estimated by the k-nearest neighbor method according to the present embodiment.
  • the upper left diagram in FIG. 61 is the same as FIG. To simplify the explanation, it is assumed that the samples s 1 to s 3 cannot be genotyped by the k-nearest neighbor method using all the samples. Samples s 1 to s 3 are selected as target samples in this order. Further, the parameter k is assumed to be 3.
  • the estimating unit 5 selects the sample s 1 as the target analyte (step S25). Then, the estimation unit 5 estimates the genotype of the sample s 1 by the k-nearest neighbor method using a known sample. In the example of FIG. 61, as the three known specimens nearest sample s 1, three known sample genotype 0 is selected. Therefore, the estimation unit 5 estimates the genotype of the sample s 1 as genotype 0 as shown in the lower left diagram (step S26). Then, the estimation unit 5 deletes the sample s 1 from the unknown sample list (step S27). Thereafter, the sample s 1 becomes a known sample of genotype 0.
  • the estimation unit 5 selects a sample s 2 as the target analyte (step S25).
  • the estimation unit 5, the k-nearest neighbor method using a known analyte, to estimate the genotype of the sample s 2.
  • two known samples of genotype 1 and one known sample (specimen s 1 ) of genotype 0 are selected as the three known samples nearest to sample s 2.
  • the estimation unit 5 as shown in the upper right diagram, the genotype of the sample s 2 estimates the genotype 1 (step S26).
  • the estimation unit 5 deletes the sample s 2 from an unknown sample list (step S27). Thereafter, the sample s 2 becomes a known sample of genotype 1.
  • the estimating unit 5 selects a sample s 3 as the target analyte (step S25). Then, the estimation unit 5 estimates the genotype of the sample s 3 by the k-nearest neighbor method using the known sample.
  • the estimation unit 5 estimates the genotype of the specimen s 3 as genotype 1 as shown in the lower right diagram (step S26). Then, the estimation unit 5 deletes the sample s 3 from the unknown sample list (step S27). Thereby, the sample s 3 becomes a known sample of the genotype 1.
  • the genotype of an unknown specimen can be estimated even when a plurality of unknown specimens are gathered at relatively close positions.
  • the k-nearest neighbor method according to the present embodiment is applicable to both the first embodiment and the second embodiment.
  • the present invention is not limited to the above-described embodiments as they are, and can be embodied by modifying the constituent elements without departing from the scope of the invention in the implementation stage.
  • various inventions can be formed by appropriately combining a plurality of constituent elements disclosed in the above embodiments. Further, for example, a configuration in which some components are deleted from all the components shown in each embodiment is also conceivable. Furthermore, you may combine suitably the component described in different embodiment.

Abstract

【課題】遺伝子型を精度よく推定することができる遺伝子型推定装置及び方法を提供する。 【解決手段】一実施形態に係る遺伝子型推定装置は、取得部と、推定部と、を備える。取得部は、遺伝子型が未知の未知検体と、遺伝子型が既知の既知検体と、を含む複数の検体の遺伝子型データのクラスタリング強度を取得する。推定部は、クラスタリング強度が第１の閾値より大きい場合、遺伝子型データに基づいて、未知検体の遺伝子型を推定し、推定結果を出力する。

Description

遺伝子型推定装置、方法、及びプログラム
　本発明の実施形態は、遺伝子型推定装置、方法、及びプログラムに関する。
　遺伝子とは、その遺伝子を持つ個人の容姿や体質に非常に大きく影響を与える遺伝情報を持ったものである。その遺伝子の情報を生体内で保持しているものがＤＮＡと呼ばれるものである。多くの生物のDNAは増殖するための自己複製や構造上の強度を保つのに優れた二重らせん構造を持ち、二本鎖ＤＮＡとして遺伝情報を保持している。
　一本のＤＮＡは４種類の塩基（base）と呼ばれる単位を並べ、塩基配列として遺伝情報を保有している。この塩基にはＡ(アデニン)、Ｔ(チミン)、Ｃ(シ卜シン)、Ｇ(グアニン)があり、この４つのパターンの組み合わせで生物の違いを表している。二本鎖ＤＮＡでは、一本のＤＮＡを基準にしてある場所のＡ(アデニン)にはもう一本のＤＮＡではＴ（チミン）を、Ｃ(シ卜シン)にはＧ(グアニン)をというようにもう一本のＤＮＡでは対応する塩基が相補的に決まっている。また、生物によって生体を構成する塩基数（base pair）は異なり、ヒ卜のＤＮＡは約３０億塩基対から構成される。
　生物の種によって通常、塩基数は同数であり、塩基配列も約９９％は共通である。その配列は個人間で異なり、塩基配列の中に同じ生物の種、ＤＮＡの同位置であるが一つの塩基が異なる場所がある。この一塩基の違いを一塩基多型（ＳＮＰ：Single Nucleotide Polymorphism）と言う。ＳＮＰは塩基配列の特定の塩基一か所が個人により異なるところを指し、突然変異(variant)の一種である。あるＳＮＰの遺伝子座における複数の塩基配列（対立遺伝子）の組み合わせパターンを遺伝子型（Genotype）という。ＳＮＰの遺伝子型は個人の容姿や体質、特定の疾患のリスク度、薬の効き方、薬の副作用の現れ方に関連している。
　そこで、複数人の各々の遺伝子型を網羅的に解析し、疾病や薬に関して関連度の高いＳＮＰを同定するＧＷＡＳ(Genome-Wide Association Study)の研究が注目を集めている。例えば、ありふれた疾患のＤＮＡにおける共通変異(common variant)を求めて、その疾患の患者群と健常者群との遺伝子型の違いを全遺伝子領域で網羅的に解析し、関連度の高いＳＮＰを同定する。その中で、複数人の遺伝子型をハイスループットで決定することのできるジェノタイピングの技術は今後のＧＷＡＳ研究の発展において重要な鍵になる。
　従来、このようなジェノタイピング技術として、ＤＮＡマイクロアレイ技術が提案されている。ＤＮＡマイクロアレイ技術は、複数人の全ゲノム情報をほぼカバーするような数万から数十万までのＳＮＰ（タグＳＮＰ）の遺伝子型をハイスループットで一度に決定することのできる遺伝子型判定技術である。
　ＤＮＡマイクロアレイ技術では、あるＳＮＰ周辺の既知塩基配列をプローブとし、検体の未知塩基配列をハイブリダイゼーションさせ蛍光強度を測定、複数の検体の蛍光強度をクラスタ空間に写像し、クラスタリングすることでそのSNPの遺伝子型を判定する。これは、各ＳＮＰにおいて、遺伝子型の情報を反映した２種類の蛍光強度の値が空間上で遺伝子型毎にまとまっているクラスタであることを利用している。
　クラスタリングを行った後、クラスタリングの信頼度としてクラスタリング強度も計算を行う。クラスタリング強度は検体群の遺伝子型が良く分離し、各遺伝子型にまとまりがある場合は高くなる。しかし、まとまりより離れた蛍光強度を持つ検体は既存のクラスタリング技術では遺伝子型を精度よく割り当てることが困難であるため、クラスタ強度に閾値を設けることで信頼性の低い検体に遺伝子型を割り当てないこともある。
米国特許第８２００４４０号明細書 米国特許第７４６７１１７号明細書
　遺伝子型を精度よく推定することができる遺伝子型推定装置、方法、及びプログラムを提供する。
　一実施形態に係る遺伝子型推定装置は、取得部と、推定部と、を備える。取得部は、遺伝子型が未知の未知検体と、遺伝子型が既知の既知検体と、を含む複数の検体の遺伝子型データのクラスタリング強度を取得する。推定部は、クラスタリング強度が第１の閾値より大きい場合、遺伝子型データに基づいて、未知検体の遺伝子型を推定し、推定結果を出力する。
第１実施形態に係る遺伝子型推定装置の機能構成を示すブロック図。 遺伝子型データの一例を示す図。 信号強度データの一例を示す図。 クラスタ座標データの一例を示す図。 クラスタリング強度データの一例を示す図。 クラスタリング強度の算出方法を説明する図。 連鎖不平衡統計データの一例を示す図。 参照ハプロタイプデータの一例を示す図。 参照遺伝子型頻度データの一例を図。 図１の遺伝子型推定装置のハードウェア構成を示すブロック図。 図１の遺伝子型推定装置の動作の概要を示すフローチャート。 ｋ近傍法による遺伝子型の推定方法の一例を示すフローチャート。 ｋ個の既知検体の選択方法を説明する図。 多数決アルゴリズムを用いた遺伝子型の推定方法を説明する図。 重み付き多数決アルゴリズムを用いた遺伝子型の推定方法を説明する図。 ｋ近傍法による遺伝子型の推定方法の他の例を示すフローチャート。 クラスタ線の作成方法を説明する図。 クラスタ線の作成方法を説明する図。 多数決アルゴリズムを用いた遺伝子型の推定方法を説明する図。 パラメータｋの設定方法を示すフローチャート。 検証用ＳＮＰを説明する図。 評価用検体及び学習用検体の一例を示す図。 推定精度の算出結果の一例を示す図。 インピュテーション法による遺伝子型の推定方法を示すフローチャート。 ＳＮＰの抽出方法を示すフローチャート。 未知検体のハプロタイプデータの一例を示す図。 参照ハプロタイプデータの選択方法を説明する図。 参照ハプロタイプデータの選択方法を示すフローチャート。 参照ハプロタイプデータの選択方法を説明する図。 ｋ近傍法及びインピュテーション法を併用した遺伝子型の推定方法を示すフローチャート。 図１の遺伝子型推定装置のＧＵＩの操作画面の一例を示す図。 図１の遺伝子型推定装置のＧＵＩの操作画面の他の例を示す図。 図１の遺伝子型推定装置のＧＵＩの操作画面の他の例を示す図。 第２実施形態に係る遺伝子型推定装置の機能構成を示すブロック図。 遺伝子型データの一例を示す図。 信号強度の閾値の学習方法の概要を示すフローチャート。 閾値組合せリストの生成方法の一例を示すフローチャート。 遺伝子型データの一例を示す図。 クラスタ番号の割当て方法の一例を示す図。 置換後の遺伝子型データの一例を示す図。 フルコールＳＮＰの信号強度データの一例を示す図。 分割された信号強度データの一例を示す図。 信号強度の統計値の一例を示す図。 信号強度の区間と遺伝子型との対応関係の一例を示す図。 閾値候補リストの一例を示す図。 閾値組合せリストの一例を示す図。 １クラスタＳＮＰの遺伝子型頻度の一例を示す図。 第１の再判定方法による再判定の前後のクラスタリングマップの一例を示す図。 第２の再判定方法による再判定の前後のクラスタリングマップの一例を示す図。 評価値を含む閾値組合せリストの一例を示す図。 閾値法による遺伝子型の推定方法の概要を示すフローチャート。 第１の推定方法を示すフローチャート。 第１の推定方法による推定の前後のクラスタリングマップの一例を示す図。 第２の推定方法を示すフローチャート。 ステップＳ１８０５における推定の前後のクラスタリングマップの一例を示す図。 ステップＳ１８０７における推定の前後のクラスタリングマップの一例を示す図。 ステップＳ１８０８における推定の前後のクラスタリングマップの一例を示す図。 ステップＳ１８０８，Ｓ１８１２における推定の前後のクラスタリングマップの一例を示す図。 複数の未知検体が相対的に近い位置に集まったクラスタリングマップの一例を示す図。 第３実施形態に係るｋ近傍法を示すフローチャート。 未知検体の遺伝子型が推定される過程に対応するクラスタリングマップを示す図である。
　以下、本発明の実施形態について図面を参照して説明する。
（第１実施形態）
　第１実施形態について、図１～図３３を参照して説明する。まず、本実施形態に係る遺伝子型推定装置（以下、「推定装置」という）の機能構成について、図１～図９を参照して説明する。図１は、本実施形態に係る推定装置の機能構成を示すブロック図である。図１に示すように、この推定装置は、検体データ記憶部１と、参照データ記憶部２と、クラスタリング強度取得部３と、クラスタリング強度判定部４と、遺伝子型推定部５と、推定結果表示部６と、を備える。
　検体データ記憶部１は、ＤＮＡマイクロアレイ技術によって遺伝子型を判定された検体に関するデータ（検体データ）を記憶する。検体データは、例えば、遺伝子型データ、信号強度データ、及びクラスタリングデータを含むが、これに限られない。
　遺伝子型データは、ＤＮＡマイクロアレイ技術による遺伝子型の判定結果を示すデータである。遺伝子型データには、検体毎かつＳＮＰ毎に判定された遺伝子型が含まれる。
　図２は、遺伝子型データの一例を示す図である。図２の遺伝子型データには、検体０１～ＮのＳＮＰｒｓ０００００１～ｒｓ９９９９９９９の遺伝子型の判定結果が含まれる。例えば、図２において、検体０１のＳＮＰｒｓ０００００１の遺伝子型は「ＣＧ」である。
　図２において、「－」は、ＤＮＡマイクロアレイ技術によって遺伝子型を判定できなかったことを示している。推定装置は、このような、ＤＮＡマイクロアレイ技術によって判定できなかった遺伝子型を推定する。
　以下では、あるＳＮＰにおいて、遺伝子型が既知の検体、すなわち、ＤＮＡマイクロアレイ技術によって遺伝子型が判定された検体を、既知検体という。また、あるＳＮＰにおいて、遺伝子型が未知の検体、すなわち、ＤＮＡマイクロアレイ技術によって遺伝子型が判定できなかった検体を、未知検体という。例えば、図２のＳＮＰｒｓ０００００２において、検体０１は既知検体であり、検体０２は未知検体である。
　信号強度データは、ＤＮＡマイクロアレイ技術による信号強度の測定結果を示すデータである。信号強度は、例えば、蛍光強度、電流、及び電圧などの測定値であるが、これに限られない。また、信号強度は、上記の測定値から算出される任意のパラメータであってもよい。信号強度データには、各検体の各ＳＮＰにおける各信号強度の値が含まれる。
　図３は、信号強度データの一例を示す図である。図３の信号強度データには、検体０１～ＮのＳＮＰｒｓ０００００１～ｒｓ９９９９９９９の信号強度ｘ１～ｘｎの値が含まれる。例えば、図３において、検体０１のＳＮＰｒｓ０００００１の信号強度ｘ１の値は０．８である。
　図３の信号強度データには、ｎ種類の信号強度の値が含まれている。ｎは、任意に設定可能であるが、ほとんどの場合２である。ｎ＝２の場合、信号強度ｘ１，ｘ２として、２種類の蛍光強度の測定値Ａ，Ｂを用いることができる。また、測定値Ａ，Ｂから以下の式により算出されるパラメータを、信号強度ｘ１，ｘ２として用いてもよい。
Figure JPOXMLDOC01-appb-M000001
　蛍光強度の測定値Ａ，Ｂをこのように変換することにより、信号強度をクラスタ空間に写像しやすくすることができる。
　クラスタリングデータは、ＤＮＡマイクロアレイ技術により遺伝子型を判定する際に、ＳＮＰ毎に行われたクラスタリングの結果を示すデータである。ＤＮＡマイクロアレイ技術によるクラスタリングは、階層的クラスタリングであってもよいし、非階層的クラスタリングであってもよい。以下では、ＤＮＡマイクロアレイ技術によるクラスタリングは、非階層的クラスタリングであるものとする。クラスタリングデータは、例えば、クラスタ座標データ、及びクラスタリング強度データを含むが、これに限られない。
　各クラスタは、ＳＮＰにおける各遺伝子型と対応するため、遺伝子型の数だけ生成される。例えば、あるＳＮＰの遺伝子型が、ＣＣ，ＣＴ，ＴＴの３つである場合、クラスタ空間上には３個のクラスタが生成される。クラスタの座標は、例えば、クラスタの重心の座標である。
　図４は、クラスタ座標データの一例を示す図である。図４のクラスタ座標データには、ＳＮＰｒｓ０００００１～ｒｓ９９９９９９のクラスタ１～ｍの座標が含まれる。各クラスタの座標は、クラスタ空間における２つの軸ｖ１，ｖ２により示されている。例えば、ＳＮＰｒｓ０００００１におけるクラスタ１の座標は、（ｖ１，ｖ２）＝（１２，３２）である。なお、クラスタ座標データには、各クラスタの座標だけでなく、クラスタ空間上における各検体の座標が含まれてもよい。また、ｖｎ＝ｘｎであってもよい。この場合、クラスタ空間は、信号強度ｘ１～ｘｎのｎ次元空間となる。
　クラスタリング強度データは、ＳＮＰ毎のクラスタリング強度を示すデータである。クラスタリング強度とは、クラスタリングの信頼度を示す指標である。クラスタリング強度が大きいほど、クラスタリングの信頼度は高い。
　図５は、クラスタリング強度データの一例を示す図である。図５のクラスタリング強度データには、ＳＮＰｒｓ０００００１～ｒｓ９９９９９９９のクラスタリング強度が含まれる。例えば、図５において、ＳＮＰｒｓ０００００１のクラスタリング強度は０．９５である。
　クラスタリング強度として、例えば、各クラスタ間の距離の平均値を用いることができる。この場合、クラスタリング強度は、クラスタ座標データから、以下の式により求めることができる。
Figure JPOXMLDOC01-appb-M000002
　式（３）において、クラスタ間距離ｉｊは任意の２つのクラスタｉ，ｊ間のユークリッド距離、（ｖｉ１，ｖｉ２）はクラスタｉの重心の座標、（ｖｊ１，ｖｊ２）はクラスタｊの重心の座標である。また、式（４）において、ｍはクラスタの数である。
　ここで、図６は、図４のクラスタ座標データから生成したクラスタリング強度データを示す図である。図６において、クラスタ間距離ｉ，ｊ及びクラスタリング強度は、クラスタリング強度が０以上１以下の値となるように規格化されている。
　参照データ記憶部２は、遺伝子に関する既知のデータ（参照データ）を記憶する。参照データは、例えば、連鎖不平衡統計データ、参照ハプロタイプデータ、及び参照遺伝子型頻度データを含むが、これに限られない。参照データのソースとして、国際ＨａｐＭａｐプロジェクトや１０００人ゲノムプロジェクトなどの、大規模なプロジェクトデータを用いることができる。
　連鎖不平衡統計データ（以下、「ＬＤデータ」という）は、ＳＮＰ同士の相関を示すデータである。図７は、ＬＤデータの一例を示す図である。図７のＬＤデータには、ＳＮＰｒｓ１２５６７８及びＳＮＰｒｓ１２９６８８のスコアと、ＳＮＰｒｓ１２５６７８及びＳＮＰｒｓ９８６７５４のスコアと、ＳＮＰｒｓ１２９６８８及びＳＮＰｒｓ９８６７５４のスコアと、が含まれる。スコアは、ＳＮＰ同士の相関の強さを示す指標である。図７のＬＤデータには、スコアとして、連鎖不平衡スコア（Ｄ′）と、相関係数（ｒ２）と、オッズ比の対数（ＬＯＤ）とが含まれる。例えば、図７において、ＳＮＰｒｓ１２５６７８とＳＮＰｒｓ１２９６８８の連鎖不平衡スコアは０．９８、相関係数は０．９６、オッズ比の対数は１８．６９である。
　参照ハプロタイプデータは、同一染色体上で統計学的に関連のあるＳＮＰの、対立遺伝子（塩基）のいずれか一方の組合せを示すデータである。すなわち、各参照ハプロタイプデータは、一部のＳＮＰにおける、蓋然性の高い塩基の組合せを示す。参照ハプロタイプデータに含まれるＳＮＰは、例えば、ＬＤデータに基づいて選択される。
　図８は、参照ハプロタイプデータの一例を示す図である。図８の参照ハプロタイプデータには、参照ハプロタイプデータｒｅｆＨＴＤ１～ｒｅｆＨＴＤ６が含まれる。各参照ハプロタイプデータには、ＳＮＰｒｓ１２３４５６，ｒｓ６２３４５６，ｒｓ９８７４５６，ｒｓ９８７１２３，ｒｓ５９８４５６，ｒｓ３８７４５６，ｒｓ９１２３４６，ｒｓ７７８４５６，ｒｓ８７３４５６，ｒｓ９８７００９の対立遺伝子が含まれる。例えば、図８において、参照ハプロタイプデータｒｅｆＨＴＤ１のＳＮＰｒｓ１２３４５６の対立遺伝子は、Ａである。
　参照遺伝子型頻度データは、ある母集団における各ＳＮＰの遺伝子型の頻度（参照遺伝子型頻度）を示すデータである。図９は、参照遺伝子型頻度データの一例を示す図である。図９の参照遺伝子型頻度データには、ＳＮＰｒｓ１２５６７８の遺伝子型ＣＣ，ＣＴ，ＴＴの頻度が含まれている。図９において、ＳＮＰｒｓ１２５６７８の対立遺伝子はＣ又はＴである。また、各遺伝子型の頻度の合計は１となる。例えば、図９において、ＳＮＰｒｓ１２５６７８の遺伝子型ＣＣの頻度は、０．４２である。
　クラスタリング強度取得部３（以下、「取得部３」という）は、ＳＮＰ毎のクラスタリング強度を取得する。検体データに図５のようなクラスタリング強度データが含まれる場合、取得部３は、検体データ記憶部１からクラスタリング強度データを取得する。
　また、検体データに図４のようなクラスタ座標データが含まれる場合、取得部３は、検体データ記憶部１からクラスタ座標データを取得し、クラスタ座標データに基づいて各ＳＮＰのクラスタリング強度を算出してもよい。クラスタリング強度の算出方法は上述の通りである。
　さらに、検体データに図２のような遺伝子型データが含まれ、参照データに図９のような参照遺伝子型頻度データが含まれる場合、取得部３は、検体データ記憶部１から遺伝子型データを取得し、参照データ記憶部２から参照遺伝子型頻度データを取得し、遺伝子型データ及び参照遺伝子型頻度データに基づいて、各ＳＮＰのクラスタリング強度を算出してもよい。クラスタリング強度の算出方法は以下の通りである。
　まず、取得部３は、遺伝子型データに基づいて、各ＳＮＰの各遺伝子型の頻度（ＤＮＡマイクロアレイ遺伝子型頻度）を算出する。ＤＮＡマイクロアレイ遺伝子型頻度は、ＤＮＡマイクロアレイ技術により判定された遺伝子型の頻度である。
　次に、取得部３は、ＤＮＡマイクロアレイ遺伝子型頻度と、参照遺伝子型頻度と、に基づいて、以下の式により各ＳＮＰのクラスタリング強度を算出する。
Figure JPOXMLDOC01-appb-M000003
　式（５）において、ｍは遺伝子型の数、ｆｉ，ｒは遺伝子型ｉの参照遺伝子型頻度、ｆｉ，Ｄは遺伝子型ｉのＤＮＡマイクロアレイ遺伝子型頻度である。遺伝子型頻度は、遺伝子型毎の確率を示すため、ｆｉ，ｒの合計及びｆｉ，ｒの合計はいずれも１である。
　例えば、あるＳＮＰの遺伝子型がＣＣ，ＣＴ，ＴＴであり、参照遺伝子型頻度がそれぞれ０．５，０．３，０．２であり、ＤＮＡマイクロアレイ遺伝子型頻度がそれぞれ０．４，０．４，０．２である場合、このＳＮＰのクラスタリング強度は、式（５）により、０．９２（＝１－ｓｑｒｔ（（（０．５－０．４）２＋（０．３－０．４）２＋（０．２－０．２）２）／３））と算出される。
　クラスタリング強度判定部４（以下、「判定部４」という）は、取得部３が取得したＳＮＰ毎のクラスタリング強度と、閾値θ１（第１の閾値）及び閾値θ２（第２の閾値）と、を比較する。閾値θ１，θ２（θ１≧θ２）は、ＤＮＡマイクロアレイ技術によるクラスタリングの信頼度を判定するために予め設定された値である。判定部４は、クラスタリング強度が閾値θ１より大きい場合、クラスタリングの信頼度は高いと判定し、閾値θ２より小さい場合、信頼度は低いと判定し、クラスタリング強度が閾値θ２以上閾値θ１以下の場合、クラスタリングの信頼度は中程度と判定する。閾値θ１，θ２は、クラスタリング強度に依存し、クラスタリング強度が０以上１以下の範囲内の値である場合、０以上１以下の範囲内の値とされる。例えば、クラスタリング強度が０以上１以下である場合、閾値θ１，θ２は、それぞれ０．８，０．４に設定される。
　なお、以下では、θ１＞θ２の場合について説明するが、θ１＝θ２であってもよい。この場合、判定部４は、クラスタリング強度が閾値θ１より大きい場合、信頼度が高いと判定し、閾値θ１以下の場合、信頼度が低いと判定する。
　遺伝子型推定部５（以下、「推定部５」という）は、遺伝子型データの各ＳＮＰにおける未知検体の遺伝子型を推定する。例えば、推定部５は、図２の遺伝子型データにおける、検体０１のＳＮＰｒｓ０００００３の遺伝子型や、検体０２のＳＮＰｒｓ０００００２の遺伝子型を推定する。
　推定部５は、判定部４によるクラスタリングの信頼度の判定結果に基づいて、推定方法を選択する。例えば、推定部５は、クラスタリング強度が閾値θ１より大きい、すなわち、ＤＮＡマイクロアレイ技術によるクラスタリングの信頼度が高い場合、遺伝子型データに基づいてｋ近傍法により遺伝子型を推定する。また、推定部５は、クラスタリング強度が閾値θ２より小さい、すなわち、ＤＮＡマイクロアレイ技術によるクラスタリングの信頼度が低い場合、遺伝子型データ及び参照データに基づいてインピュテーション法により遺伝子型を推定する。さらに、推定部５は、クラスタリング強度が閾値θ２以上閾値θ１以下、すなわち、ＤＮＡマイクロアレイ技術によるクラスタリングの信頼度が中程度の場合、ｋ近傍法及びインピュテーション法を併用して遺伝子型を推定する。そして、推定部５は推定結果を出力する。遺伝子型の推定方法の具体例について、詳しくは後述する。
　推定結果表示部６（以下、「表示部６」という）は、推定部５による推定結果を表示する。表示部６は、推定結果とともに、遺伝子型データや、推定の際に用いられた各種の情報を表示してもよい。
　次に、本実施形態に係る推定装置のハードウェア構成について、図１０を参照して説明する。本実施形態に係る推定装置は、図１０に示すように、コンピュータ１００により構成される。コンピュータ１００は、ＣＰＵ（中央演算装置）１０１と、入力装置１０２と、表示装置１０３と、通信装置１０４と、記憶装置１０５と、とを備え、これらはバス１０６により相互に接続されている。
　ＣＰＵ１０１は、コンピュータ１００の制御装置及び演算装置である。ＣＰＵ１０１は、バス１０６を介して接続された各装置（例えば、入力装置１０２、通信装置１０４、記憶装置１０５）から入力されたデータやプログラムに基づいて演算処理を行い、演算結果や制御信号を、バス１０６を介して接続された各装置（例えば、表示装置１０３、通信装置１０４、記憶装置１０５）に出力する。ＣＰＵ１０１は、コンピュータ１００のＯＳ（オペレーティングシステム）や、遺伝子型推定プログラム（以下、「推定プログラム」という）などを実行し、コンピュータ１００を構成する各装置を制御する。推定プログラムとは、コンピュータ１００に、推定装置の上述の各機能構成を実現させるプログラムである。ＣＰＵ１０１が推定プログラムを実行することにより、コンピュータ１００が推定装置として機能する。
　入力装置１０２は、コンピュータ１００に情報を入力するための装置である。入力装置１０２は、例えば、キーボード、マウス、及びタッチパネルであるが、これに限られない。ユーザは、入力装置１０２を用いることにより、閾値θ１，θ２などの情報を入力することができる。
　表示装置１０３は、ＣＰＵ１０１から出力されたデータ等に基づき、画像や映像等を表示するための装置である。表示装置１０３は、例えば、ＬＣＤ（液晶ディスプレイ）、ＣＲＴ（ブラウン管）、及びＰＤＰ（プラズマディスプレイ）であるが、これに限られない。表示部６は、表示装置１０３を用いて構成することができる。
　通信装置１０４は、コンピュータ１００が外部装置と無線又は有線で通信するための装置である。通信装置１０４は、例えば、モデム、ハブ、及びルータであるが、これに限られない。検体データや参照データなどの情報は、通信装置１０４を介して外部装置から受信することにより入力することができる。また、ＣＰＵ１０１から出力された演算結果等のデータを、外部装置へ送信することもできる。
　記憶装置１０５は、コンピュータ１００のＯＳや、推定プログラム、推定プログラムの実行に必要なデータ、及びＣＰＵ１０１による推定プログラムの実行により生成し出力されたデータなどを記憶する記憶媒体である。記憶装置１０５には、主記憶装置と外部記憶装置とが含まれる。主記憶装置は、例えば、ＲＡＭ、ＤＲＡＭ、ＳＲＡＭであるが、これに限られない。また、外部記憶装置は、ハードディスク、光ディスク、フラッシュメモリ、及び磁気テープであるが、これに限られない。検体データ記憶部１や参照データ記憶部２は、記憶装置１０５を用いて構成することができる。
　なお、コンピュータ１００は、ＣＰＵ１０１、入力装置１０２、表示装置１０３、通信装置１０４、及び記憶装置１０５を、１つ又は複数備えてもよいし、プリンタやスキャナなどの周辺機器を接続されていてもよい。
　また、推定装置は、単一のコンピュータ１００により構成されてもよいし、相互に接続された複数のコンピュータ１００からなるシステムとして構成されてもよい。
　さらに、推定プログラムは、コンピュータ１００の記憶装置１０５に予め記憶されていてもよいし、ＣＤ－ＲＯＭなどの記憶媒体に記憶されていてもよいし、インターネット上にアップロードされていてもよい。いずれの場合も、推定プログラムをコンピュータ１００にインストールして実行することにより、推定装置を構成することができる。
　次に、本実施形態に係る推定装置の動作について、図１１～図３０を参照して説明する。図１１は、本実施形態に係る推定装置の動作の概要を示すフローチャートである。
　ステップＳ１において、取得部３は、検体データ記憶部１から遺伝子型データを取得する。
　ステップＳ２において、取得部３は、遺伝子型データから対象ＳＮＰの遺伝子型データを抽出する。対象ＳＮＰとは、未知検体を含むＳＮＰである。例えば、図２の遺伝子型データの場合、取得部３は、ＳＮＰｒｓ０００００２，ｒｓ０００００３の遺伝子型データを抽出する。
　ステップＳ３において、取得部３は、各対象ＳＮＰのクラスタリング強度ＣＳを取得する。上述の通り、取得部３は、遺伝子型データ、クラスタリングデータ、及び参照遺伝子型頻度データなどに基づいて、クラスタリング強度ＣＳを取得することができる。
　ステップＳ４において、判定部４は、取得部３から各対象ＳＮＰのクラスタリング強度ＣＳを取得し、閾値θ１と比較する。閾値θ１は、対象ＳＮＰ毎に同一であってもよいし、異なってもよい。
　ＣＳ＞θ１の場合、判定部４は、クラスタリングの信頼度は高いと判定し（ステップＳ４のＹＥＳ）、処理はステップＳ５に進む。
　ステップＳ５において、推定部５は、遺伝子型データに基づいて、ｋ近傍法により未知検体の遺伝子型を推定する。ｋ近傍法による遺伝子型の推定方法については後述する。
　ＣＳ≦θ１の場合（ステップＳ４のＮＯ）、処理はステップＳ６に進む。
　ステップＳ６において、判定部４は、取得部３から取得した各対象ＳＮＰのクラスタリング強度ＣＳと、閾値θ２と、を比較する。閾値θ２は、対象ＳＮＰ毎に同一であってもよいし、異なってもよい。
　ＣＳ＜θ２の場合、判定部４は、クラスタリングの信頼度は低いと判定し（ステップＳ６のＹＥＳ）、処理はステップＳ７に進む。
　ステップＳ７において、推定部５は、遺伝子型データ及び参照データに基づいて、インピュテーション法により未知検体の遺伝子型を推定する。インピュテーション法による遺伝子型の推定方法については後述する。
　ＣＳ≧θ２の場合（ステップＳ６のＮＯ）、判定部４は、クラスタリングの信頼度は中程度と判定し、処理はステップＳ８に進む。
　ステップＳ８において、推定部５は、ｋ近傍法とインピュテーション法とを併用して、未知検体の遺伝子型を推定する。ｋ近傍法とインピュテーション法とを併用した遺伝子型の推定方法については後述する。
　ステップＳ５，Ｓ７，Ｓ８において未知検体の遺伝子型が推定された後、ステップＳ９において、表示部６は、推定部５による推定結果を表示する。
　以下、ｋ近傍法、インピュテーション法、及びこれらを併用した方法による、遺伝子型の推定方法について詳細に説明する。
　まず、ステップＳ５におけるｋ近傍法による遺伝子型の推定方法について、図１２～図２３を参照して説明する。ここでいうｋ近傍法とは、ｋ個の最近傍のサンプルの遺伝子型に基づいて、未知検体の遺伝子型を推定する方法のことである。以下では、サンプルとして、既知検体及びクラスタ線を用いる方法について、それぞれ説明する。
　図１２は、ｋ近傍法による遺伝子型の推定方法の一例を示すフローチャートである。図１２の推定方法では、サンプルとして既知検体を用いる。
　ステップＳ５１１において、推定部５は、検体データ記憶部１から、対象ＳＮＰの未知検体Ｓの遺伝子型データ及び信号強度データを取得する。以下では、未知検体Ｓは１つであるものとするが、未知検体Ｓが複数ある場合には、各未知検体Ｓについて、以降の処理が行われる。
　ステップＳ５１２において、推定部５は、検体データ記憶部１から、対象ＳＮＰの既知検体群ＳＴの遺伝子型データ及び信号強度データを取得する。既知検体群ＳＴは、対象ＳＮＰに含まれる既知検体Ｓｉの集合のことである。
　ステップＳ５１３において、推定部５は、既知検体群ＳＴに含まれる各既知検体Ｓｉについて、距離ｄｉを算出する。距離ｄｉは、未知検体Ｓと、既知検体Ｓｉと、の距離である。距離ｄｉは、例えば、未知検体Ｓの信号強度データが（ｘ１，ｘ２，・・・，ｘｎ）、既知検体Ｓｉの信号強度データが（ｘｉ１，ｘｉ２，・・・，ｘｉｎ）の場合、以下の式により算出される。
Figure JPOXMLDOC01-appb-M000004
　ステップＳ５１４において、推定部５は、既知検体群ＳＴの中から、最近傍のｋ個の既知検体Ｓｉ、すなわち、距離ｄｉが小さい順にｋ個の既知検体Ｓｉを選択する。パラメータｋは、予め設定される任意の自然数である。パラメータｋの設定方法については、後述する。
　図１３は、既知検体Ｓｉの選択方法を説明する図である。図１３は、対象ＳＮＰのクラスタリングマップの一例を示している。図１３において、信号強度はｘ１，ｘ２の２種類（ｎ＝２）、パラメータｋは５（ｋ＝５）、星印は未知検体Ｓ、丸は遺伝子型がＣＣの既知検体、三角は遺伝子型がＣＧの既知検体、四角は遺伝子型がＧＧの既知検体である。図１３の場合、ステップＳ５１４において、距離ｄｉが小さい順に、遺伝子型がＣＣの既知検体３個と、遺伝子型がＣＧの既知検体が２個選択される。
　ステップＳ５１５において、推定部５は、選択したｋ個の既知検体Ｓｉの遺伝子型に基づいて、未知検体Ｓの遺伝子型を推定する。
　推定部５は、例えば、多数決アルゴリズムを用いて未知検体Ｓの遺伝子型を推定する。すなわち、推定部５は、選択したｋ個の既知検体Ｓｉの遺伝子型のうち、最も検体数（投票数）が多い遺伝子型を、未知検体Ｓの遺伝子型として推定する。
　図１４は、多数決アルゴリズムを用いた遺伝子型の推定方法を説明する図である。図１４では、５個の既知検体Ｓｉ（ｉ＝１～５）が選択されており、それぞれの遺伝子型は、ＡＧ，ＧＧ，ＡＧ，ＡＧ，ＡＡである。この場合、各遺伝子型ＡＧ，ＧＧ，ＡＡの投票数は、それぞれ３，１，１となるため、未知検体Ｓの遺伝子型は、投票数が最大であるＡＧと推定される。
　また、推定部５は、重み付き多数決アルゴリズムを用いて未知検体Ｓの遺伝子型を推定してもよい。この場合、推定部５は、まず、選択した各既知検体Ｓｉの重みを算出する。既知検体Ｓｉの重みとして、既知検体Ｓｉにおける遺伝子型を判定されたＳＮＰの割合を用いることができる。例えば、ＤＮＡマイクロアレイ技術により、２０万個のＳＮＰのうち１５万個の遺伝子型が判定された既知検体Ｓｉの重みは、０．７５となる。
　推定部５は、各既知検体Ｓｉの重みを投票数として利用して、投票数が最も多い遺伝子型を未知検体Ｓの遺伝子型として推定する。図１５は、重み付き多数決アルゴリズムを用いた遺伝子型の推定方法を説明する図である。図１５では、５個の既知検体Ｓｉ（ｉ＝１～５）が選択されており、それぞれの遺伝子型は、ＡＧ，ＧＧ，ＡＧ，ＡＧ，ＡＡであり、それぞれの重みは、０．６，０．４，０．９，０．７，０．５である。この場合、各遺伝子型ＡＧ，ＧＧ，ＡＡの投票数は、それぞれ２．２，０．４，０．５となるため、未知検体Ｓの遺伝子型は、投票数が最大であるＡＧと推定される。
　図１６は、ｋ近傍法による遺伝子型の推定方法の他の例を示すフローチャートである。図１６の推定方法では、サンプルとしてクラスタ線を用いる。
　ステップＳ５２１において、推定部５は、検体データ記憶部１から、対象ＳＮＰの未知検体Ｓの遺伝子型データ及び信号強度データを取得する。ステップＳ５２１は、上述のステップＳ５１１と同様である。
　ステップＳ５２２において、推定部５は、検体データ記憶部１から、対象ＳＮＰの既知検体群ＳＴの遺伝子型データ及び信号強度データを取得する。ステップＳ５２２は、上述のステップＳ５１２と同様である。
　ステップＳ５２３において、推定部５は、既知検体群ＳＴの信号強度データに基づいて、クラスタ線Ｃｉを作成する。クラスタ線Ｃｉとは、クラスタ空間上の各クラスタ（各遺伝子型）に含まれる既知検体を、直線や曲線により近似したものである。クラスタ線Ｃｉは、クラスタ空間上における既知検体の座標を回帰分析することにより作成することができる。回帰分析は、線形回帰であってもよいし、非線形回帰であってもよい。
　図１７は、クラスタ線Ｃｉの作成方法を説明する図である。図１７は、対象ＳＮＰのクラスタリングマップの一例を示している。図１７において、信号強度はｘ１，ｘ２の２種類（ｎ＝２）、パラメータｋは５（ｋ＝５）、星印は未知検体Ｓ、丸は遺伝子型がＣＣの既知検体、三角は遺伝子型がＣＧの既知検体、四角は遺伝子型がＧＧの既知検体である。また、クラスタ線Ｃｉは、直線であり、クラスタ毎に１本ずつ作成されている。図１７の場合、クラスタ線Ｃｉは、以下の式により表される。
Figure JPOXMLDOC01-appb-M000005
　式（７）において、ｍｉ，ｃｉは、定数であり、回帰分析により求められる。推定部５は、クラスタ毎に線形回帰分析を行うことにより、上記のようなクラスタ線Ｃｉを作成することができる。
　また、図１８に示すように、クラスタ線Ｃｉは、クラスタ毎にそれぞれ複数本ずつ作成されてもよい。この場合、推定部５は、各クラスタを複数のサブクラスタに分割し、それぞれのサブクラスタについて、式（７）のようにクラスタ線Ｃｉを作成すればよい。
　なお、クラスタ線Ｃｉは、直線に限られず、曲線であってもよい。また、クラスタ線Ｃｉの本数は任意に選択可能である。
　ステップＳ５２４において、推定部５は、各クラスタ線Ｃｉについて、距離Ｄｉを算出する。距離Ｄｉは、未知検体Ｓと、クラスタ線Ｃｉと、の距離である。距離Ｄｉは、例えば、未知検体Ｓの信号強度データが（ｘｓ１，ｘｓ２）、クラスタ線Ｃｉがｘ２＝ｍｉｘ１＋ｃｉの場合、以下の式により算出される。
Figure JPOXMLDOC01-appb-M000006
　ステップＳ５２５において、推定部５は、作成した複数のクラスタ線Ｃｉの中から、最近傍のｋ本のクラスタ線Ｃｉ、すなわち、距離Ｄｉが小さい順にｋ本のクラスタ線Ｃｉを選択する。
　例えば、図１７において、ｋ＝１の場合、距離Ｄｉが最も小さい遺伝子型ＣＣのクラスタ線Ｃｉが選択される。また、図１８において、ｋ＝３の場合、距離Ｄｉが小さい順に、遺伝子型がＣＣのクラスタ線２本と、遺伝子型がＣＧのクラスタ線１本が選択される。
　ステップＳ５２６において、推定部５は、選択したｋ本のクラスタ線Ｃｉの遺伝子型に基づいて、未知検体Ｓの遺伝子型を推定する。
　推定部５は、例えば、多数決アルゴリズムを用いて未知検体Ｓの遺伝子型を推定する。すなわち、推定部５は、選択したｋ本のクラスタ線Ｃｉの遺伝子型のうち、最もクラスタ線数（投票数）が多い遺伝子型を、未知検体Ｓの遺伝子型として推定する。
　図１９は、多数決アルゴリズムを用いた遺伝子型の推定方法を説明する図である。図１９では、５本のクラスタ線Ｃｉ（ｉ＝１～５）が選択されており、それぞれの遺伝子型は、ＡＧ，ＧＧ，ＡＧ，ＡＧ，ＡＡである。この場合、各遺伝子型ＡＧ，ＧＧ，ＡＡの投票数は、それぞれ３，１，１となるため、未知検体Ｓの遺伝子型は、投票数が最大であるＡＧと推定される。
　また、推定部５は、重み付き多数決アルゴリズムを用いて未知検体Ｓの遺伝子型を推定してもよい。この場合、推定部５は、まず、選択した各クラスタ線Ｃｉの重みを算出する。クラスタ線Ｃｉの重みとして、既知検体Ｓｉにおける遺伝子型を判定されたＳＮＰの割合の、クラスタごとの平均値を用いることができる。推定部５は、各クラスタ線Ｃｉの重みを投票数として利用して、投票数が最も多い遺伝子型を未知検体Ｓの遺伝子型として推定する。
　ここで、ｋ近傍法で用いるパラメータｋの設定方法について、図２０～図２３を参照して説明する。図２０は、パラメータｋの設定方法を示すフローチャートである。本実施形態において、推定部５は、クロス検証によりパラメータｋを設定する。
　ステップＳ５３１において、推定部５は、検体データ記憶部１から、１つ又は複数の検証用ＳＮＰの遺伝子型データ及び信号強度データを取得する。検証用ＳＮＰとは、クラスタリング強度ＣＳが大きく、かつ、全ての検体が既知検体であるＳＮＰのことである。検証用ＳＮＰは、例えば、クラスタリング強度ＣＳが閾値θ１より大きいＳＮＰである。
　図２１は、検証用ＳＮＰを説明する図である。図２１の遺伝子型データにおいて、ＳＮＰｒｓ００００１，ｒｓ０００００３の検体は、全て既知検体である。これらのＳＮＰのクラスタリング強度ＣＳが大きい場合、推定部５は、ＳＮＰｒｓ００００１，ｒｓ０００００３を検証用ＳＮＰとして抽出し、これらの遺伝子型データ及び信号強度データを取得する。
　ステップＳ５３２において、推定部５は、評価用検体及び学習用検体を選択する。評価用検体とは、未知検体として扱う検体である。学習用検体とは、既知検体として扱う検体である。評価用検体として選択した検体の遺伝子型は、クロス検証のための正解データとして利用される。
　図２２は、評価用検体及び学習用検体の一例を示す図である。図２２において、評価用検体として検体０１～１０が選択され、検体１１～Ｎが学習用検体として選択されている。なお、評価用検体及び学習用検体は、任意に選択可能である。
　ステップＳ５３３において、推定部５は、パラメータｋの候補ｋ′を複数設定する。推定部５は、パラメータｋの候補ｋ′として、任意の自然数を設定することができる。
　ステップＳ５３４において、推定部５は、学習用検体の遺伝子型データ及び信号強度データに基づいて、各評価用検体の遺伝子型を推定する。この際、推定部５は、パラメータｋとして各候補ｋ′を用いたｋ近傍法により、評価用検体の遺伝子型を推定する。
　ステップＳ５３５において、推定部５は、クロス検証により、各候補ｋ′の推定精度を算出する。すなわち、推定部５は、評価用検体の遺伝子型の推定結果と、評価用検体の既知の遺伝子型と、を比較し、遺伝子型を正しく推定された評価用検体の割合を算出する。
　図２３は、推定精度の算出結果の一例を示す図である。図２３に示すように、推定精度は、各検証用ＳＮＰの各候補ｋ′について算出される。例えば、図２３において、ＳＮＰｒｓ０００００１のｋ′＝１における推定精度は０．８である。また、図２３に示すように、検証用ＳＮＰが複数ある場合、推定部５は、各候補ｋ′の推定精度の平均値（平均推定精度）を算出してもよい。
　ステップＳ５３６において、推定部５は、推定精度が最大の候補ｋ′を、パラメータｋとして設定する。また、推定部５は、複数の検証用ＳＮＰについてクロス検証を行った場合、平均推定精度が最大の候補ｋ′を、パラメータｋとして設定してもよい。例えば、図２３の場合、パラメータｋは、平均推定精度が最大となる５に設定される。
　次に、ステップＳ７におけるインピュテーション法による遺伝子型の推定方法について、図２４～図２９を参照して説明する。図２４は、インピュテーション法による遺伝子型の推定方法を示すフローチャートである。
　ステップＳ７１において、推定部５は、参照データ記憶部２を参照して、参照データ記憶部２に対象ＳＮＰのＬＤデータがあるか否か確認する。インピュテーション法では、対象ＳＮＰのＬＤデータを利用するため、対象ＳＮＰのＬＤデータが無い場合（ステップＳ７１のＮＯ）、処理はステップＳ５へ進み、ｋ近傍法により未知検体Ｓの遺伝子型を推定する。ｋ近傍法による推定方法は上述の通りである。
　一方、対象ＳＮＰのＬＤデータがある場合（ステップＳ７１のＹＥＳ）、処理はステップＳ７２に進む。
　ステップＳ７２において、推定部５は、参照データ記憶部２から、対象ＳＮＰのＬＤデータを取得する。
　ステップＳ７３において、推定部５は、対象ＳＮＰのＬＤデータを参照して、スコアが高いＬ個以上のＳＮＰを抽出する。図２５は、ステップＳ７３におけるＳＮＰの抽出方法を示すフローチャートである。
　ステップＳ７３１において、推定部５は、対象ＳＮＰのＬＤデータから、推定用ＳＮＰのＬＤデータを抽出する。推定用ＳＮＰとは、クラスタリング強度ＣＳが大きく、かつ、全ての検体が既知検体であるＳＮＰのことである。推定用ＳＮＰは、例えば、クラスタリング強度ＣＳが閾値θ１より大きいＳＮＰである。
　対象ＳＮＰのＬＤデータには、対象ＳＮＰに対する他の複数のＳＮＰのスコアが含まれる。推定部５は、対象ＳＮＰのＬＤデータにスコアが含まれる各ＳＮＰの遺伝子型データ及びクラスタリングデータを参照して、推定用ＳＮＰのＬＤデータを抽出する。
　ステップＳ７３２において、推定部５は、スコアの閾値ＬＤθを、推定用ＳＮＰのＬＤデータに含まれるスコアの最高値に設定する。スコアは、ＬＤデータに含まれる連鎖不平衡スコア、相関係数、及びオッズ比の対数などの中から任意に選択可能である。
　ステップＳ７３３において、推定部５は、推定用ＳＮＰのＬＤデータを参照して、推定用ＳＮＰの中からスコアが閾値ＬＤθ以上のＳＮＰを抽出する。
　ステップＳ７３４において、推定部５は、抽出したＳＮＰの数が、所定値Ｌ以上であるか判定する。Ｌは任意に設定可能である。抽出したＳＮＰの数がＬより小さい場合（ステップＳ７３４のＮＯ）、処理はステップＳ７３５に進む。
　ステップＳ７３５において、推定部５は、閾値ＬＤθを低くする（ＬＤθ＝ＬＤθ－Δ）。閾値ＬＤθの減少量Δは、例えば、０．０１である。閾値ＬＤθを低くした後、処理は、ステップＳ７３３に進む。そして、推定部５は、抽出したＳＮＰの数がＬ以上になるまで、ステップＳ７３３～Ｓ７３５の処理を繰り返す。
　抽出したＳＮＰの数がＬ以上の場合（ステップＳ７３４のＹＥＳ）、ＳＮＰの抽出処理は終了し、処理はステップＳ７４に進む。以上の処理により、推定部５は、スコアが閾値ＬＤθ以上のＬ個以上の推定用ＳＮＰを抽出することができる。
　ステップＳ７４において、推定部５は、参照データ記憶部２から、ステップＳ７３で抽出したＳＮＰの参照ハプロタイプデータを取得する。
　ステップＳ７５において、推定部５は、検体データ記憶部１から、未知検体Ｓの遺伝子型データを取得する。
　ステップＳ７６において、推定部５は、未知検体Ｓの遺伝子型データから、未知検体Ｓのハプロタイプデータを作成する。未知検体Ｓのハプロタイプデータは、遺伝子型データからフェージング（相化）アルゴリズムを用いてＳＮＰの遺伝子型を抽出し、各染色体に存在する対立遺伝子の配列を決定することにより作成することができる。フェージングアルゴリズムとして、例えば、ＢＥＡＧＬＥ，ｆａｓｔＰＨＡＳＥ，ＩＭＰＵＴＥｖ２，ＭＡＣＨ，ＳｈａｐｅＩＴを用いることができる。
　図２６は、未知検体Ｓのハプロタイプデータの一例を示す図である。図２６に示すように、フェージングアルゴリズムにより、未知検体Ｓの遺伝子型データから、２つのハプロタイプデータＨＴＤ１，ＨＴＤ２が作成される。未知検体Ｓのハプロタイプデータにおいて、遺伝子型が未知のＳＮＰの対立遺伝子対は不明であるため、「－」で示されている。
　ステップＳ７７において、推定部５は、ステップＳ７４で取得した参照ハプロタイプデータの中から、未知検体Ｓの２つのハプロタイプデータＨＴＤ１，ＨＴＤ２に最も類似した参照ハプロタイプデータをそれぞれ選択する。ハプロタイプデータＨＴＤ１，ＨＴＤ２と参照ハプロタイプデータとが類似するとは、遺伝子型が未知のＳＮＰ以外のＳＮＰにおける対立遺伝子の配列が類似することをいう。
　図２７は、参照ハプロタイプデータの選択方法を説明する図である。例えば、ステップＳ７４において、図２７の参照ハプロタイプデータが抽出された場合、推定部５は、ハプロタイプデータＨＴＤ１に最も類似する参照ハプロタイプデータとして、参照ハプロタイプデータｒｅｆＨＴＤ５を選択し、ハプロタイプデータＨＴＤ２に最も類似する参照ハプロタイプデータとして、参照ハプロタイプデータｒｅｆＨＴＤ３を選択する。参照ハプロタイプデータの選択方法の詳細は後述する。
　ステップＳ７８において、推定部５は、選択した２つの参照ハプロタイプデータにおける対象ＳＮＰの対立遺伝子に基づいて、未知検体Ｓの遺伝子型を推定する。例えば、図２７のように参照ハプロタイプデータを選択した場合、推定部５は、ＳＮＰｒｓ９８７００９の一方の対立遺伝子をＧと推定し、他方の対立遺伝子をＡと推定する。そして、推定部５は、これらの対立遺伝子に基づいて、未知検体ＳのＳＮＰｒｓ９８７００９の遺伝子型をＡＧと推定する。
　図２８は、ステップＳ７７における、参照ハプロタイプデータの選択方法を示すフローチャートである。
　ステップＳ７７１において、推定部５は、抽出した参照ハプロタイプデータ及び未知検体Ｓのハプロタイプデータの対立遺伝子を数値に置換する。図２９は、対立遺伝子を数値に置換された参照ハプロタイプデータ及び未知検体Ｓのハプロタイプデータの一例を示す図である。図２９において、対立遺伝子Ａ，Ｃ，Ｇ，Ｔが、数値１，２，３，４にそれぞれ置換されている。
　ステップＳ７７２において、推定部５は、距離ｄｈｉを算出する。距離ｄｈｉは、未知検体Ｓの各ハプロタイプデータと、各参照ハプロタイプデータと、の間の距離である。距離ｄｈｉは、例えば、以下の式により算出される。
Figure JPOXMLDOC01-appb-M000007
　式（９）において、ｐは未知検体Ｓのハプロタイプデータに含まれるＳＮＰのうち、遺伝子型が未知のＳＮＰを除くＳＮＰの数、ｓｉｊ（ｊ＝１～ｐ）は参照ハプロタイプデータｉのｊ番目のＳＮＰの数値、ｓｊ（ｊ＝１～ｐ）は未知検体Ｓのハプロタイプデータのｊ番目のＳＮＰの数値である。
　例えば、図２９の場合、ハプロタイプデータＨＴＤ１と参照ハプロタイプデータｒｅｆＨＴＤ１との距離ｄｈｉは、０．３５（＝ｓｑｒｔ（（（１－１）２＋（４－４）２＋（３－３）２＋（４－１）２＋（２－２）２＋（４－４）２＋（２－３）２＋（１－１）２＋（２－２）２））／９）と算出される。
　ステップＳ７７３において、推定部５は、未知検体Ｓの各ハプロタイプデータについて、距離ｄｈｉが最小の参照ハプロタイプデータを、最も類似した参照ハプロタイプデータとして選択する。
　次に、ステップＳ８におけるｋ近傍法とインピュテーション法とを併用した遺伝子型の推定方法について、図３０を参照して説明する。図３０は、ｋ近傍法とインピュテーション法を併用した遺伝子型の推定方法を示すフローチャートである。
　ステップＳ８１において、推定部５は、ｋ近傍法により未知検体Ｓの遺伝子型を推定し、１つ又は複数の遺伝子型の候補からなる遺伝子型群ＧＴ１を取得する。遺伝子型群ＧＴ１に含まれる遺伝子型の候補の数をα個とすると、遺伝子型群ＧＴ１は、例えば、遺伝子型の候補として、投票数が大きい順にα個の遺伝子型を選択したり、α個のパラメータｋにより遺伝子型を推定したりすることにより取得できる。
　ステップＳ８２において、推定部５は、インピュテーション法により未知検体Ｓの遺伝子型を推定し、１つ又は複数の遺伝子型の候補からなる遺伝子型群ＧＴ２を取得する。遺伝子型群ＧＴ２に含まれる遺伝子型の候補の数をβ個とすると、遺伝子型群ＧＴ２は、例えば、未知検体Ｓのハプロタイプデータ毎に類似している参照ハプロタイプデータをβ個選択して遺伝子型を推定することにより取得できる。
　ステップＳ８３において、推定部５は、遺伝子型群ＧＴ１，ＧＴ２に含まれる遺伝子型の候補の中から、多数決アルゴリズムを用いて、未知検体Ｓの遺伝子型を推定する。多数決アルゴリズムの投票数として、遺伝子型群ＧＴ１，ＧＴ２に含まれる遺伝子型の数を用いることができる。
　以上説明した通り、本実施形態に係る推定装置及び方法は、ＤＮＡマイクロアレイ技術により判定できなかった遺伝子型を、ＤＮＡマイクロアレイ技術によるクラスタリングの信頼度に応じた方法により推定する。すなわち、信頼度が低い場合には、参照データを利用したインピュテーション法により推定し、信頼度が高い場合には、ＤＮＡマイクロアレイ技術により判定された遺伝子型データを利用したｋ近傍法により推定する。これにより、本実施形態に係る推定装置及び方法は、遺伝子型を精度よく推定することができる。
　なお、以上説明した本実施形態に係る推定装置は、ＧＵＩ（Graphical User Interface）により操作可能であるのが好ましい。図３１は、表示部６により表示されるＧＵＩの操作画面の一例を示す図である。図３１は、既知検体をサンプルとして用いたｋ近傍法により遺伝子型を推定する場合のＧＵＩである。図３１に示すように、このＧＵＩは、ＳＮＰ選択部Ｇ１と、検体選択部Ｇ２と、ｋ値選択部Ｇ３と、検体一覧表示部Ｇ４と、選択結果表示部Ｇ５と、遺伝子型表示部Ｇ６と、を備える。
　ＳＮＰ選択部Ｇ１は、ユーザが対象ＳＮＰを選択するためのドロップダウンリストである。ＳＮＰ選択部Ｇ１のドロップダウンリストには、未知検体を含む全てのＳＮＰのＩＤが含まれる。ＳＮＰ選択部Ｇ１には、ユーザにより選択された対象ＳＮＰのＩＤが表示される。
　検体選択部Ｇ２は、ユーザが、遺伝子型を推定する未知検体Ｓを選択するためのドロップダウンリストである。検体選択部Ｇ２のドロップダウンリストには、対象ＳＮＰの遺伝子型データに含まれる全ての未知検体ＳのＩＤが含まれる。検体選択部Ｇ２のドロップダウンリストの内容は、ユーザにより選択された対象ＳＮＰに応じて変化する。検体選択部Ｇ２には、ユーザにより選択された未知検体ＳのＩＤが表示される。
　ｋ値選択部Ｇ３は、ユーザがパラメータｋを設定するためのドロップダウンリストである。ｋ値選択部Ｇ３のドロップダウンリストには、パラメータｋの値の候補が複数含まれる。ｋ値選択部Ｇ３には、ユーザにより設定されたパラメータｋの値が表示される。図３１では、パラメータｋは５に設定されている。なお、ｋ値選択部Ｇ３には、推定精度が最も高いパラメータｋの値がデフォルト値として設定されているのが好ましい。
　検体一覧表示部Ｇ４は、ユーザにより選択された対象ＳＮＰの遺伝子型データに含まれる既知検体（既知検体群ＳＴに含まれる既知検体Ｓｉ）のＩＤの一覧と、各既知検体Ｓｉとユーザにより選択された未知検体Ｓとの間の距離ｄｉと、を表示する。図３１で表示された距離ｄｉは、上述の式（６）により算出された距離である。
　選択結果表示部Ｇ５は、検体一覧表示部Ｇ４に表示された既知検体Ｓｉの中から選択された、距離ｄｉが小さいｋ個の既知検体Ｓｉの、ＩＤ、距離ｄｉ、遺伝子型、及び重みを表示する。図３１において、パラメータｋは５のため、５つの既知検体Ｓｉが表示されている。重みは、多数決アルゴリズムに用いられる重みであり、デフォルト値として１．０が設定されている。重みが１．０の場合、重み無しの多数決アルゴリズムとなる。重み付きアルゴリズムを用いる場合には、重みとして、上述の方法で算出された各既知検体Ｓｉの重みが表示される。
　遺伝子型表示部Ｇ６は、ｋ近傍法による遺伝子型の推定結果を表示する。図３１において、推定された遺伝子型はＡＧである。
　図３２は、推定装置のＧＵＩの操作画面の他の例を示す図である。図３２のＧＵＩを備える推定装置は、参照データ記憶部２に、ＳＮＰと疾患との関連を示す情報を記憶している。図３２に示すように、このＧＵＩは、患者選択部Ｇ７と、疾患選択部Ｇ８と、ＳＮＰ情報表示部Ｇ９と、遺伝子型情報表示部Ｇ１０と、を備える。
　患者選択部Ｇ７は、ユーザが患者を選択するためのドロップダウンリストである。ここでいう患者は、ＤＮＡマイクロアレイ技術により遺伝子型を判定された検体に対応する。患者選択部Ｇ７のドロップダウンリストには、複数の患者（検体）のＩＤが含まれる。患者選択部Ｇ７には、ユーザにより選択された患者のＩＤが表示される。
　疾患選択部Ｇ８は、ユーザが疾患を選択するためのドロップダウンリストである。疾患選択部Ｇ８のドロップダウンリストには、参照データ記憶部２に記憶された複数の疾患の名称が含まれる。疾患選択部Ｇ８には、ユーザにより選択された疾患の名称が表示される。
　ＳＮＰ情報表示部Ｇ９は、ユーザにより選択された患者及び疾患に関連するＳＮＰ情報を表示する。ＳＮＰ情報には、染色体の種類、ＳＮＰのＩＤ、遺伝子座、疾患関連ＳＮＰ、オッズ比（ＯＲ）、Ａｄｄｒ情報、及び塩基配列情報などが含まれる。ここでいうオッズ比は、医学的な臨床試験の結果を示す方法として用いられる尺度であり、疾患へのかかりやすさを２つの群で比較して示す統計的な尺度である。また、Ａｄｄｒ情報及び塩基配列情報については後述する。ＳＮＰ情報は、参照データ記憶部２に記憶されている。
　図３２において、ＳＮＰ情報表示部Ｇ９は、１番～５番、１０番～１２番、及びＸＹ染色体を表示している。各染色体上の斜線部分は、疾患関連ＳＮＰのうち遺伝子型が既知のＳＮＰを示し、ドット部分は疾患関連ＳＮＰのうち遺伝子型が未知のＳＮＰを示している。図３２において、各染色体上のＳＮＰは、コマンドボタンとなっており、ユーザが選択（クリック）すると、そのＳＮＰに関するＡｄｄｒ情報や、ＳＮＰの周辺の塩基配列が表示される。
　Ａｄｄｒ情報には、ＳＮＰが所属する染色体（Chromosome）の番号、染色体上での遺伝子座（Position）、ＳＮＰが所属する遺伝子の名称（Gene）、及びＳＮＰのＩＤが含まれる。コマンドボタンによりユーザが選択したＳＮＰのＡｄｄｒ情報は、ＳＮＰ情報表示部Ｇ９のＡｄｄｒ欄に表示される。
　塩基配列情報は、ＳＮＰを除く遺伝子座の塩基配列データである。コマンドボタンによりユーザがＳＮＰを選択すると、選択されたＳＮＰの周辺の塩基配列が塩基配列情報から抽出され、抽出された塩基配列の範囲に含まれるＳＮＰの遺伝子型データが検体データ記憶部１から抽出され、ＳＮＰ情報表示部Ｇ９の塩基配列欄に表示される。図３２の塩基配列において、ＳＮＰの対立遺伝子Ａ，Ｂは、［Ａ／Ｂ］という形式で表示されている。例えば、図３２において、ＳＮＰｒｓ５４７９８４の対立遺伝子は、ＧとＴである。
　遺伝子型情報表示部Ｇ１０は、ユーザにより選択されたＳＮＰに関する遺伝子型情報を表示する。遺伝子型情報は、検体データ記憶部１に記憶された各種のデータから生成される。
　ユーザにより選択されたＳＮＰの遺伝子型が既知である場合、遺伝子型情報表示部Ｇ１０は、図３２に示すように、選択されたＳＮＰのクラスタリングマップ（Genotype Clustering）や、ＨａｐＭａｐによる遺伝子型の比率などを表示する。
　クラスタリングマップは、検体データ記憶部１に予め記憶されていてもよいし、検体データ記憶部１に記憶された信号強度データなどから生成されてもよい。また、ＨａｐＭａｐによる遺伝子型の比率とは、患者が所属する民族団体における、選択されたＳＮＰの遺伝子型の比率のことである。ＨａｐＭａｐによる遺伝子型の比率は、参照遺伝子型頻度データから抽出することができる。
　これに対して、ユーザにより選択されたＳＮＰの遺伝子型が未知である場合、遺伝子型情報表示部Ｇ１０は、図３３に示すように、推定部５による遺伝子型の推定結果や、遺伝子型を推定する過程を示すデータを表示する。
　図３３において、遺伝子型情報表示部Ｇ１０に表示されたImputationは、表示ラベルであり、推定部５による遺伝子型の推定方法を示している。インピュテーション法により遺伝子型が推定された場合、図３３に示すように、遺伝子型表示部Ｇ１０は、参照ハプロタイプデータ、患者（検体）のハプロタイプデータ、フェージング後のハプロタイプデータ、インピュテーション後のハプロタイプデータなどを表示する。インピュテーション後のハプロタイプデータには、推定された遺伝子型が含まれる。ＳＮＰ情報表示部Ｇ９の塩基配列欄には、推定された遺伝子型（対立遺伝子）が表示される。また、遺伝子型情報表示部Ｇ１０は、図３２と同様に、ＨａｐＭａｐによる遺伝子型の比率を表示してもよい。
　なお、推定部５がｋ近傍法により遺伝子型を推定した場合には、遺伝子型情報表示部Ｇ１０は、図３１における検体一覧表示部Ｇ４、選択結果表示部Ｇ５、及び遺伝子型表示部Ｇ６などに表示された情報を表示してもよい。
（第２実施形態）
　第２実施形態について、図３４～図５８を参照して説明する。本実施形態では、閾値法を用いた遺伝子型の推定方法について説明する。ここでいう閾値法とは、信号強度の区間と遺伝子型との対応関係を学習し、学習した対応関係に基づいて、各検体の遺伝子型を推定する、遺伝子型の推定方法のことである。信号強度の各区間は、信号強度の閾値により規定される。閾値法について、詳しくは後述する。
　まず、本実施形態に係る推定装置の機能構成について、図３４及び図３５を参照して説明する。図３４は、本実施形態に係る推定装置の機能構成を示すブロック図である。図３４に示すように、この推定装置は、検体データ記憶部１と、推定部５と、表示部６と、を備える。以下、第１実施形態との相違点について説明する。
　本実施形態において、検体データ記憶部１は、検体データとして、遺伝子型データと、信号強度データと、を記憶しており、クラスタリングデータを記憶していない。また、推定装置は、参照データ記憶部２と、取得部３と、判定部４と、を備えない。
　これは、閾値法では、クラスタリングデータ、参照データ、及びクラスタリング強度を使用しないためである。後述するように、本実施形態に係る推定方法と、第１実施形態に係る推定方法と、を併用する場合には、検体データ記憶部１にクラスタリングデータを記憶させるとともに、推定装置に参照データ記憶部２と、取得部３と、判定部４と、を設ければよい。
　また、推定部５は、閾値学習部５１と、閾値法推定部５２と、を備える。
　閾値学習部５１（以下、「学習部５１」という）は、フルコール（Fullcall）ＳＮＰの信号強度に基づいて、閾値法で用いる信号強度の区間と、遺伝子型と、の対応関係を学習する。具体的には、学習部５１は、信号強度の区間を規定する閾値と、遺伝子型と、の対応関係を学習する。
　フルコールＳＮＰとは、ＤＮＡマイクロアレイ技術によって全ての検体の遺伝子型が判定されたＳＮＰ、すなわち、全ての検体が既知検体であるＳＮＰのことである。これに対して、ＤＮＡマイクロアレイ技術によって少なくとも１つの検体の遺伝子型が判定されなかったＳＮＰ、すなわち、少なくとも１つの未知検体を含むＳＮＰを、ノーコール（Nocall）ＳＮＰと称する。
　ここで、フルコールＳＮＰ及びノーコールＳＮＰについて、図３５を参照して具体的に説明する。図３５は、検体データ記憶部１に記憶された遺伝子型データの一例を示す図である。図３５の例では、ＳＮＰｒｓ００００１，ｒｓ９９９９９９は、全ての検体の遺伝子型が判定されている。したがって、ＳＮＰｒｓ００００１，ｒｓ９９９９９９は、フルコールＳＮＰである。これに対して、ＳＮＰｒｓ０００００２，ｒｓ０００００３は、それぞれ検体０２，０１が未知検体である。したがって、ＳＮＰｒｓ０００００２，ｒｓ０００００３は、ノーコールＳＮＰである。このように、学習部５１は、遺伝子型データを参照することにより、フルコールＳＮＰ及びノーコールＳＮＰを把握することができる。
　学習部５１は、閾値を学習するために、遺伝子型データからフルコールＳＮＰを抽出し、信号強度データからフルコールＳＮＰの各検体の信号強度を抽出する。信号強度データにｎ種類の信号強度の値が含まれている場合、学習部５１は、学習する対象となるいずれか１種類の信号強度を抽出すればよい。学習部５１が閾値を学習する信号強度の種類は、推定装置のユーザが任意に設定可能である。以下では、学習部５１が信号強度ｘ１を抽出し、信号強度ｘ１の閾値を学習する場合を例に説明する。
　また、学習部５１が学習する閾値の数は、推定装置のユーザが任意に設定可能であり、１つでもよいし、複数でもよい。閾値の数は、各ＳＮＰに含まれる遺伝子型の種類に応じて設定されるのが好ましい。
　信号強度の区間は、閾値の数より１つ多く規定される。したがって、各ＳＮＰに含まれる遺伝子型の種類の最大値がＸ個の場合、学習部５１は、例えば、Ｘ－１個の閾値を学習することが考えられる。
　以下では、学習部５１が閾値ｘｌ（第１の閾値）と、ｘｌより大きい閾値ｘｒ（第２の閾値）と、の２つの閾値を学習する場合を例に説明する。これは、第１実施形態と同様に、各ＳＮＰに最大３種類の遺伝子型が含まれる場合を想定している。
　信号強度の閾値及びその学習方法について、詳しくは後述する。
　閾値法推定部５２（以下、「推定部５２」という）は、学習部５１が学習した信号強度の区間（閾値）と遺伝子型との対応関係に基づいて、ノーコールＳＮＰの各検体の遺伝子型を推定する。上述の通り、ノーコールＳＮＰには、未知検体及び既知検体が含まれる。したがって、本実施形態では、未知検体の遺伝子型の推定だけでなく、既知検体の遺伝子型の推定（再判定）も行われる。
　例えば、第１実施形態では、図３５のＳＮＰｒｓ０００００２は、未知検体である検体０２の遺伝子型のみが推定される。これに対して、本実施形態では、ＳＮＰｒｓ０００００２は、未知検体である検体０２の遺伝子型を推定されるとともに、既知検体である検体０１，Ｎの遺伝子型も推定される。閾値法を用いた遺伝子型の推定方法について、詳しくは後述する。
　なお、本実施形態に係る推定装置のハードウェア構成は、第１実施形態と同様である。すなわち、コンピュータ１００が推定プログラムを実行することにより、推定装置の上述の各機能構成が実現される。
　次に、本実施形態に係る推定装置の動作について、図３６～図５８を参照して具体的に説明する。以下では、学習部５１の動作と、推定部５２の動作と、について順に説明する。
　まず、学習部５１の動作について、図３６～図５０を参照して説明する。以下では、学習部５１が信号強度ｘ１の２つの閾値ｘｌ，ｘｒ（３つの区間）を学習する場合を例に説明する。図３６は、閾値の学習方法の概要を示すフローチャートである。各ステップについて、詳しくは後述する。
　まず、ステップＳ１０において、学習部５１は、閾値組合せリストを生成する。閾値組合せリストとは、複数の閾値組合せを含むリストである。閾値組合せとは、閾値候補の組合せのことである。２つの閾値ｘｌ，ｘｒ（ｘｌ＜ｘｒ）を学習する場合、閾値組合せは、閾値ｘｌの候補と、閾値ｘｒの候補と、の組合せとなる。
　次に、ステップＳ１１において、学習部５１は、閾値組合せリストに含まれる各閾値組合せを評価するための遺伝子型頻度を算出する。
　続いて、ステップＳ１２において、学習部５１は、各閾値組合せに含まれる閾値候補と、遺伝子型頻度と、に基づいて、各閾値組合せの評価値を算出する。
　そして、ステップＳ１３において、学習部５１は、閾値組合せリストに含まれる閾値組合せの中から、評価値が最大の閾値組合せを選択する。選択された閾値組合せに含まれる各閾値候補が、閾値法で遺伝子型を推定するための閾値として採用される。
　以下、ステップＳ１０～Ｓ１３について、詳細に説明する。図３７は、ステップＳ１０における閾値組合せリストの生成方法の一例を示すフローチャートである。
　ステップＳ１０１において、学習部５１は、検体データ記憶部１から、全てのＳＮＰの遺伝子型データと、全てのＳＮＰの信号強度ｘ１の信号強度データと、を取得する。
　図３８は、取得された遺伝子型データ及び信号強度データの一例を示す図である。図３８の例では、ＳＮＰｒｓ０００００１～ｒｓ９９９９９９９の遺伝子型データ及び信号強度データが取得されている。
　ステップＳ１０２において、学習部５１は、ステップＳ１０１で取得した遺伝子型データに含まれる各遺伝子型を、クラスタ番号に置換する。クラスタ番号とは、クラスタリングマップ上の各クラスタの相対位置に応じて割当てられた値である。学習部５１は、まず、各ＳＮＰの各クラスタに、クラスタ番号を割当てる。
　図３９は、クラスタ番号の割当て方法の一例を示す図である。図３９の例では、クラスタリングマップ上の右に位置するクラスタから順に、クラスタ番号０，１，２が割当てられている。これは、重心の信号強度ｘ１が大きいクラスタから順に、クラスタ番号０，１，２を割当てることに相当する。
　学習部５１は、ステップＳ１０１で取得した遺伝子型データ及び信号強度データから、各クラスタの重心の信号強度ｘ１を算出して、クラスタ番号を割当てることができる。また、検体データ記憶部１にクラスタ座標データが記憶されている場合には、学習部５１は、検体データ記憶部１からクラスタ座標データを取得し、取得したクラスタ座標データを参照して、クラスタ番号を割当ててもよい。
　クラスタ番号は、ＳＮＰ毎に、共通の方法で割当てられる。したがって、あるＳＮＰでは、遺伝子型ＡＡのクラスタにクラスタ番号０が割当てられ、他のＳＮＰでは、遺伝子型ＣＣのクラスタにクラスタ番号０が割当てられる、ということが起こりえる。
　学習部５１は、各クラスタにクラスタ番号を割当てた後、各検体の遺伝子型を、その検体が含まれるクラスタに割当てられたクラスタ番号に置換する。例えば、学習部５１は、あるクラスタに、クラスタ番号０が割当てられた場合、そのクラスタに含まれる各検体の遺伝子型を、０に置換する。
　図４０は、遺伝子型をクラスタ番号に置換後の遺伝子型データの一例を示す図である。図４０の遺伝子型データは、図３８の遺伝子型データに対応している。図４０において、０，１，２は、それぞれ遺伝子型に対応するクラスタ番号であり、－１は、遺伝子型が判定されていないことに対応するクラスタ番号である。
　例えば、ＳＮＰｒｓ０００００１において、遺伝子型ＣＧは、クラスタ番号１に置換され、遺伝子型ＣＣは、クラスタ番号２に置換されている。また、ＳＮＰｒｓ０００００２において、遺伝子型ＡＴは、クラスタ番号１に置換され、遺伝子型ＴＴはクラスタ番号２に置換されている。これは、ＳＮＰｒｓ０００００１における遺伝子型ＣＧのクラスタの相対位置と、ＳＮＰｒｓ０００００２における遺伝子型ＡＴの相対位置と、が等しいことを示している。
　なお、図３９の例では、クラスタ番号は、各クラスタに、信号強度ｘ１の降順で割当てられたが、信号強度ｘ１の昇順で割当てられてもよいし、信号強度ｘ２の降順又は昇順で割当てられてもよい。
　以下、各検体の遺伝子型を、クラスタ番号を用いて表す。例えば、図４０のＳＮＰｒｓ０００００１の検体０１の遺伝子型は、遺伝子型１となる。
　ステップＳ１０３において、学習部５１は、置換後の遺伝子型データを参照して、フルコールＳＮＰを抽出する。例えば、図４０の遺伝子型データを参照すると、フルコールＳＮＰとして、ＳＮＰｒｓ０００００１，ｒｓ９９９９９８，ｒｓ９９９９９９が抽出される。
　ステップＳ１０４において、学習部５１は、ステップＳ１０１で取得した信号強度データから、ステップＳ１０３で抽出したフルコールＳＮＰの、信号強度データを抽出する。
　図４１は、抽出されたフルコールＳＮＰの信号強度データの一例を示す図である。図４１の信号強度データは、図３８の信号強度データからフルコールＳＮＰの信号強度データを抽出したものである。
　ステップＳ１０５において、学習部５１は、ステップＳ１０４で抽出したフルコールＳＮＰの信号強度データを、１クラスタＳＮＰの信号強度データと、複数クラスタＳＮＰの信号強度データと、に分割する。
　１クラスタＳＮＰとは、判定結果として１種類の遺伝子型をしか含まないＳＮＰのことである。すなわち、１クラスタＳＮＰとは、全ての検体が同一の遺伝子型と判定されたＳＮＰのことである。これに対して、複数クラスタＳＮＰとは、判定結果として複数種類の遺伝子型を含むＳＮＰのことである。判定結果として、２種類以上の遺伝子型を含むＳＮＰは、全て複数クラスタＳＮＰに含まれる。ただし、ここでいう２種類以上の遺伝子型には、上述の遺伝子型－１は含まれない。
　学習部５１は、信号強度データを分割するために、各ＳＮＰの遺伝子型データに含まれる遺伝子型の種類を数える。学習部５１は、あるＳＮＰの遺伝子型データに含まれる遺伝子型の種類が１種類（例えば、遺伝子型１）である場合、そのＳＮＰを１クラスタＳＮＰと判定する。また、学習部５１は、あるＳＮＰの遺伝子型データに含まれる遺伝子型の種類が２種類以上である場合、そのＳＮＰを複数クラスタＳＮＰと判定する。学習部５１は、こうして得られたＳＮＰの判定結果に基づいて、フルコールＳＮＰの信号強度データを分割する。
　図４２は、１クラスＳＮＰの信号強度データ及び複数クラスタＳＮＰの信号強度データの一例を示す図である。図４２の信号強度データは、図４１の信号強度データを分割したものである。図４０からわかるように、ＳＮＰｒｓ９９９９９８は、１種類の遺伝子型０のみを含む１クラスタＳＮＰであり、ＳＮＰｒｓ０００００１，ｒｓ９９９９９９は、３種類の遺伝子型０，１，２を含む複数クラスタＳＮＰである。このため、図４２に示すように、１クラスタＳＮＰの信号強度データには、ＳＮＰｒｓ９９９９９８の信号強度データが含まれ、複数クラスタＳＮＰの信号強度データには、ＳＮＰｒｓ０００００１，９９９９９９の信号強度データが含まれている。
　このように、１クラスタＳＮＰ及び複数クラスタＳＮＰの信号強度データを分割するのは、１クラスタＳＮＰにおけるクラスタの分布と、複数クラスタＳＮＰにおけるクラスタの分布と、の間の相違が大きいためである。信号強度の閾値の学習と、学習した閾値を用いた遺伝子型の推定と、を１クラスタＳＮＰ及び複数クラスタＳＮＰのそれぞれについて行うことにより、遺伝子型の推定精度を向上させることができる。
　なお、本実施形態において、推定装置は、１クラスタＳＮＰ及び複数クラスタＳＮＰをまとめて処理することも可能である。この場合、学習部５１は、ステップＳ１０５における信号強度データの分割を行なわず、以降の処理についても、１クラスタＳＮＰ及び複数クラスタＳＮＰをまとめて処理すればよい。
　ステップＳ１０６において、学習部５１は、ステップＳ１０５で分割した１クラスタＳＮＰの信号強度データを参照して、信号強度の統計値を遺伝子型毎に算出する。また、学習部５１は、ステップＳ１０５で分割した複数クラスタＳＮＰの信号強度データを参照して、信号強度の統計値を遺伝子型毎に算出する。
　信号強度の統計値には、最小値、平均値、最大値、及び標準偏差値が含まれる。学習部５１は、置換後の遺伝子型データを参照して、分割後の１クラスタＳＮＰ（又は複数クラスタＳＮＰ）の信号強度データから遺伝子型０の信号強度を抽出し、抽出した信号強度の統計値を算出することにより、１クラスタＳＮＰ（又は複数クラスタＳＮＰ）の遺伝子型０の信号強度の統計値を算出することができる。他の遺伝子型も同様の方法で算出される。
　図４３は、１クラスタＳＮＰの信号強度の統計値と、複数クラスタＳＮＰの信号強度の統計値と、の一例を示す図である。図４３の例では、１クラスタＳＮＰの遺伝子型２の信号強度の最小値は、－６．２９である。
　ステップＳ１０７において、学習部５１は、１クラスタＳＮＰの各遺伝子型の信号強度の平均値に基づいて、１クラスタのＳＮＰの信号強度の区間と、遺伝子型と、を対応させる。また、学習部５１は、複数クラスタＳＮＰの各遺伝子型の信号強度の平均値に基づいて、複数クラスタＳＰＮの信号強度の区間と、遺伝子型と、を対応させる。
　学習部５１が、２つの信号強度ｘｌ、ｘｒを学習すると、信号強度の区間が３つ形成される。学習部５１は、信号強度が小さい区間から順に、信号強度の平均値が小さい遺伝子型を対応させる。
　図４４は、信号強度の区間と、遺伝子型と、の対応関係の一例を示す図である。図４４の対応関係は、図４３の信号強度の平均値に基づいている。図４３の例では、信号強度の平均値は、遺伝子型２，１，０の順に小さい。このため、図４４の例では、信号強度が小さい区間から順に、遺伝子型２，１，０が対応付けられている。具体的には、信号強度がｘｌ未満の区間に遺伝子型２が対応付けられ、信号強度がｘｌ以上ｘｒ以下の区間に遺伝子型１が対応付けられ、信号強度がｘｒより大きい区間に遺伝子型０が対応付けられている。
　ステップＳ１０８において、学習部５１は、ステップＳ１０６で算出した１クラスタＳＮＰの統計値に基づいて、１クラスタＳＮＰの閾値候補リストを生成する。また、学習部５１は、ステップＳ１０６で算出した複数クラスタＳＮＰの統計値に基づいて、複数クラスタＳＮＰの閾値候補リストを生成する。閾値候補リストとは、複数の閾値候補を含むリストのことである。閾値候補とは、信号強度の閾値ｘｌ、ｘｒの候補のことである。
　１クラスタＳＮＰ（又は複数クラスタＳＮＰ）の閾値候補は、例えば、１クラスタＳＮＰ（又は複数クラスタＳＮＰ）の各遺伝子型の信号強度の、最小値、平均値、最大値、平均値＋Ｎ×標準偏差（Ｎは整数）であるが、これに限られない。
　図４５は、１クラスタＳＮＰの閾値候補リストと、複数クラスタＳＮＰの閾値候補リストと、の一例を示す図である。図４５の閾値候補リストは、図４３の統計値に対応しており、９個の閾値候補（各遺伝子型の最小値、平均値、最大値）を含んでいる。例えば、１クラスタＳＮＰの閾値候補リストに含まれる閾値候補－６．２９は、１クラスタＳＮＰの遺伝子型２の信号強度の最小値である。
　また、１クラスタＳＮＰ（又は複数クラスタＳＮＰ）の閾値候補は、例えば、１クラスタＳＮＰ（又は複数クラスタＳＮＰ）の統計値の最大値と最小値との間を、等間隔に分割する値であってもよい。この場合、各閾値候補ｘｉは、以下の式で算出される。
Figure JPOXMLDOC01-appb-M000008
　式（１０），（１１）において、ｎは閾値候補リストに含まれる閾値候補の数、ｘｍｉｎは統計値の最小値、ｘｍａｘは統計値の最大値、ｄは閾値候補の間隔である。図４３の例では、１クラスタＳＮＰの統計値の最小値ｘｍｉｎは、遺伝子型２の信号強度の最小値である－６．２９に相当し、最大値ｘｍａｘは、遺伝子型０の信号強度の最大値である７．４６に相当する。
　なお、閾値候補リストに含まれる閾値候補は、上記のものに限られず、信号強度の統計値から任意の方法で生成可能である。また、閾値候補リストには、閾値候補として、予め設定された任意の値が含まれてもよい。
　ステップＳ１０９において、学習部５１は、ステップＳ１０８で生成した１クラスタＳＮＰの閾値候補リストを参照して、１クラスタＳＮＰの閾値組合せリストを生成する。また、学習部５１は、ステップＳ１０８で生成した複数クラスタＳＮＰの閾値候補リストを参照して、複数クラスタＳＮＰの閾値組合せリストを生成する。
　閾値組合せリストとは、上述の通り、複数の閾値組合せを含むリストである。学習部５１は、閾値候補リストに含まれる閾値候補を組み合わせて閾値組合せを生成し、複数の閾値組合せを含む閾値組合せリストを生成する。
　閾値候補リストに閾値候補がｎ個含まれ、閾値組合せに閾値候補がｒ個含まれる場合、最大でｎ！／（ｎ－ｒ）！ｒ！個の閾値組合せが生成される。したがって、１クラスタＳＮＰの閾値候補リストに９個の閾値候補が含まれ、閾値組合せに２つの閾値候補ｘｌ、ｘｒが含まれる場合、最大で３６個の閾値組合せが生成される。
　図４６は、１クラスタＳＮＰの閾値組合せリストの一例を示す図である。図４６の閾値組合せリストは、図４５の閾値候補リストに対応している。図４６の閾値組合せリストには、３６個の閾値組合せが含まれている。図４６の例では、例えば、閾値組合せ１は、（ｘｌ、ｘｒ）＝（－６．２９，７．４６）である。図４６のような閾値組合せリストが、複数クラスタＳＮＰについても生成される。
　次に、ステップＳ１１における遺伝子型頻度の算出方法について、詳細に説明する。閾値組合せリストを生成した後、学習部５１は、１クラスタＳＮＰの遺伝子型頻度と、複数クラスタＳＮＰの遺伝子型頻度と、をそれぞれ算出する。遺伝子型頻度とは、隣接する２つの閾値候補により規定される信号強度の区間に含まれる信号強度、を有する検体の数のことである。遺伝子型頻度は、遺伝子型毎に算出される。
　学習部５１は、ステップＳ１０５で分割した信号強度データと、置換後の遺伝子型データと、ステップＳ１０８で生成した閾値候補リストと、を参照することにより、１クラスタＳＮＰ及び複数クラスタＳＮＰの各区間の各遺伝子型の遺伝子型頻度を算出することができる。
　例えば、図４２の１クラスタＳＮＰの信号強度データを参照すると、ＳＮＰｒｓ９９９９９８の検体０１の信号強度は０．３である。この信号強度は、図４５の１クラスタＳＮＰの閾値候補リストを参照すると、閾値候補０．６９と閾値候補２．１１との間の区間に含まれることがわかる。そして、図４０の遺伝子型データを参照すると、ＳＮＰｒｓ９９９９９８の検体０１の遺伝子型は２である。結果として、１クラスタＳＮＰの０．６９と２．１１との間の区間の遺伝子型２の遺伝子型頻度が、１加算される。
　学習部５１は、１クラスタＳＮＰの信号強度データに含まれる各信号強度を参照して、上記のように各区間の各遺伝子型の遺伝子型頻度を加算していき、１クラスタＳＮＰの遺伝子型頻度を算出する。複数クラスタＳＮＰについても同様の方法で遺伝子型頻度が算出される。
　図４７は、１クラスタＳＮＰの遺伝子型頻度の一例を示す図である。図４７の遺伝子型頻度は、図４５の閾値候補リストに対応している。図４７の例では、閾値候補－２．１１と閾値候補－１．７９との間の区間において、遺伝子型２の遺伝子型頻度が大きく、遺伝子型１の遺伝子型頻度が小さく、遺伝子型０の遺伝子型頻度が０である。
　次に、ステップＳ１２における評価値の算出方法及びステップＳ１３における閾値組合せの選択方法について、詳細に説明する。遺伝子型頻度を算出した後、学習部５１は、１クラスタＳＮＰの各閾値組合せの評価値を算出し、算出した評価値に基づいて、１クラスタＳＮＰの閾値組合せを選択する。また、学習部５１は、複数クラスタＳＮＰの各閾値組合せの評価値を算出し、算出した評価値に基づいて、複数クラスタＳＮＰの閾値組合せを選択する。
　まず、学習部５１は、フルコールＳＮＰにおける１クラスタＳＮＰの遺伝子型を、各閾値組合せに基づいてそれぞれ再判定する。以下、２つの再判定方法について説明する。
　第１の再判定方法では、学習部５１は、１クラスタＳＮＰの閾値組合せリストから１つの閾値組合せを選択し、１クラスタＳＮＰの検体の信号強度及び遺伝子型を抽出する。そして、学習部５１は、信号強度がｘｌより小さい検体の遺伝子型を遺伝子型２と再判定し、信号強度がｘｌ以上ｘｒ以下の検体の遺伝子型を遺伝子型１と再判定し、信号強度がｘｒより大きい検体の遺伝子型を遺伝子型０と再判定する。
　図４８は、第１の再判定方法による再判定の前後のクラスタリングマップの一例を示す図である。図４８において、丸は各検体を示し、丸中の値は遺伝子型の値を示している。図４８の上図は、再判定前の遺伝子型を示し、下図は再判定後の遺伝子型を示している。図４８からわかるように、第１の再判定方法では、閾値候補ｘｌと閾値候補ｘｒとの間の区間に信号強度が含まれる検体の遺伝子型は、全ての遺伝子型１と再判定される。
　このように、第１の再判定方法では、閾値組合せのみが用いられ、遺伝子型頻度が用いられない。したがって、第１の再判定方法を採用する場合には、ステップＳ１１を省略することができる。
　これに対して、第２の再判定方法では、遺伝子型頻度が用いられる。具体的には、第２の再判定方法では、信号強度がｘｌ以上ｘｒ以下の検体の遺伝子型が、遺伝子型頻度を用いた多数決アルゴリズムにより再判定される。すなわち、信号強度がｘｌ以上ｘｒ以下の検体は、その検体の信号強度が含まれる区間における遺伝子型頻度が最大の遺伝子型に再判定される。
　例えば、（ｘｌ，ｘｒ）＝（－２．１１，２．３３）であり、ある検体の信号強度が－２．００であり、図４７に示す遺伝子型頻度が得られた場合、この検体の遺伝子型は、－２．１１と－１．７９との間の区間において遺伝子型頻度が最大である、遺伝子型２に再判定される。
　図４９は、第２の再判定方法による再判定の前後のクラスタリングマップの一例を示す図である。図４９に示すように、第２の再判定方法では、閾値候補ｘｌと閾値候補ｘｒとの間の区間に信号強度が含まれる検体が、必ずしも遺伝子型１と再判定されるわけではない。
　学習部５１は、上記のような再判定方法により、１クラスタＳＮＰの遺伝子型を再判定した後、再判定前の遺伝子型と、再判定後の遺伝子型と、の一致率を算出する。例えば、図４８の例では、再判定の前後で、１６個の検体のうち、１２個の検体の遺伝子型が一致するため、一致率は０．７５である。また、図４９の例では、再判定の前後で、１６個の検体のうち、１３個の検体の遺伝子型が一致するため、一致率は０．８１である。
　こうして算出された一致率が、再判定に使用された閾値組合せの評価値として用いられる。学習部５１は、以上の方法により、１クラスタＳＮＰの閾値組合せリストに含まれる各閾値組合せの評価値を算出する。また、学習部５１は、複数クラスタＳＮＰの閾値組合せリストに含まれる各閾値組合せの評価値も、同様の方法で算出する。
　図５０は、評価値を含む閾値組合せリストの一例を示す図である。図５０の例では、閾値組合せ１の評価値は０．８０、閾値組合せ２４の評価値は０．９７である。
　学習部５１は、１クラスタＳＮＰの閾値組合せリストに含まれる閾値組合せの中から、評価値が最大の閾値組合せを選択する。選択された閾値組合せに含まれる各閾値候補が、閾値法で１クラスタＳＮＰの遺伝子型を推定するための閾値として学習される。
　また、学習部５１は、複数クラスタＳＮＰの閾値組合せリストに含まれる閾値組合せの中から、評価値が最大の閾値組合せを選択する。選択された閾値組合せに含まれる各閾値候補が、閾値法で複数クラスタＳＮＰの遺伝子型を推定するための閾値として学習される。
　例えば、図５０の例では、閾値組合せリストの中で評価値が最大の閾値組合せは、閾値組合せ２４であるため、閾値組合せ２４に含まれる閾値候補が、遺伝子型を推定するための閾値として学習される。すなわち、閾値法により遺伝子型を推定する際、閾値ｘｌとして－０．８０が使用され、閾値ｘｒとして２．１１が使用される。
　以上のように閾値を学習することにより、閾値法による遺伝子型の推定精度を向上させることができる。これは、閾値組合せの評価値（一致率）が高いほど、その閾値組合せによる遺伝子型の推定結果と、フルコールＳＮＰにおける遺伝子型の推定結果と、が近くなるためである。一般に、フルコールＳＮＰにおける遺伝子型の推定精度は高いため、フルコールＳＮＰにおける推定結果と近い推定結果を得られる閾値組合せほど、推定精度が高くなる。
　次に、推定部５２による閾値を用いた遺伝子型の推定方法について、図５１～図５８を参照して説明する。以下では、信号強度の閾値ｘｌ、ｘｒは、学習済みであるものとする。図５１は、閾値法による遺伝子型の推定方法の概要を示すフローチャートである。
　まず、ステップＳ１４において、推定部５２は、検体データ記憶部１から、全てのＳＮＰの遺伝子型データと、全てのＳＮＰの信号強度ｘ１の信号強度データと、を取得する。
　次に、ステップＳ１５において、推定部５２は、ステップＳ１４で取得した遺伝子型データに含まれる各遺伝子型を、クラスタ番号に置換する。遺伝子型の置換方法は、ステップＳ１０２で説明した通りである。
　続いて、ステップＳ１６において、推定部５２は、置換後の遺伝子型データを参照して、ノーコールＳＮＰを抽出する。例えば、図４０の遺伝子型データを参照すると、ノーコールＳＮＰとして、ＳＮＰｒｓ０００００２，ｒｓ０００００３が抽出される。
　そして、ステップＳ１７において、推定部５２は、ステップＳ１６で抽出したノーコールＳＮＰの中から、対象ＳＮＰを選択する。対象ＳＮＰとは、閾値法により遺伝子型を推定する対象となるＳＮＰのことである。本実施形態では、上述の通り、閾値の学習は、全てのフルコールＳＮＰを用いて行われるが、遺伝子型の推定は、個々のノーコールＳＮＰごとに行われる。対象ＳＮＰの選択方法は任意である。
　その後、ステップＳ１８において、推定部５２は、ステップＳＳ１７で選択した対象ＳＮＰの各検体の遺伝子型を、閾値法により推定する。ステップＳ１８について、詳しくは後述する。
　対象ＳＮＰとして未選択のノーコールＳＮＰがある場合（ステップＳ１９のＹＥＳ）、推定部５２は、未選択のノーコールＳＮＰの中から、次の対象ＳＮＰを選択する（ステップＳ１７）。以降、未選択のノーコールＳＮＰがなくなるまで、ステップＳ１７～Ｓ１９が繰り返される。
　そして、対象ＳＮＰとして未選択のノーコールＳＮＰがなくなると（ステップＳ１９のＮＯ）、推定部５２は、遺伝子型の推定処理を終了する。
　ここで、ステップＳ１８における、遺伝子型の推定方法について、詳細に説明する。以下、２つの推定方法についてそれぞれ説明する。
　まず、第１の推定方法について説明する。第１の推定方法では、推定部５２は、信号強度の閾値ｘｌ，ｘｒのみを用いて、対象ＳＮＰの全検体の遺伝子型を推定する。図５２は、第１の推定方法を示すフローチャートである。
　ステップＳ１８０１において、推定部５２は、置換後の遺伝子型データ及び信号強度データの中から、対象ＳＮＰの遺伝子型データ及び信号強度データを抽出する。
　ステップＳ１８０２において、推定部５２は、対象ＳＮＰの遺伝子型データを参照して、対象ＳＮＰのクラスタ数を取得する。クラスタ数の取得方法は、上述の通りである。すなわち、推定部５２は、対象ＳＮＰの遺伝子型データに含まれる遺伝子型の種類を数えることにより、対象ＳＮＰのクラスタ数を取得する。
　なお、信号強度の閾値ｘｌ，ｘｒが、１クラスタＳＮＰ及び複数クラスタＳＮＰのそれぞれについて個別に学習されていない場合、ステップＳ１８０２は省略されてもよい。この場合、推定部５２は、１クラスタＳＮＰ及び複数クラスタＳＮＰをまとめて、以降の処理を実行すればよい。
　ステップＳ１８０３において、推定部５２は、学習部５１から、対象ＳＮＰのクラスタ数に応じた信号強度の閾値ｘｌ，ｘｒを取得する。推定部５２は、対象ＳＮＰのクラスタ数が１の場合、１クラスタＳＮＰの閾値ｘｌ，ｘｒを取得し、対象ＳＮＰのクラスタ数が２以上の場合、複数クラスタＳＮＰの閾値ｘｌ，ｘｒを取得する。
　ステップＳ１８０４において、推定部５２は、ステップＳ１８０３で取得した閾値ｘｌ，ｘｒに基づいて、対象ＳＮＰの各検体の遺伝子型を推定する。具体的には、推定部５２は、信号強度ｘｉがｘｌより小さい（ｘｉ＜ｘｌ）検体の遺伝子型を遺伝子型２と推定し、信号強度がｘｌ以上ｘｒ以下（ｘｌ≦ｘｉ≦ｘｒ）の検体の遺伝子型を遺伝子型１と推定し、信号強度がｘｒより大きい（ｘｒ＜ｘｉ）検体の遺伝子型を遺伝子型０と推定する。
　図５３は、第１の推定方法による推定の前後のクラスタリングマップの一例を示す図である。図５３の上図は、第１の推定方法による推定前の遺伝子型を示し、下図は第１の推定方法による推定後の遺伝子型を示している。図５３からわかるように、第１の推定方法では、未知検体及び既知検体を含む全ての検体の遺伝子型が、閾値ｘｌ，ｘｒに基づいて推定される。
　その後、推定部５２は、対象ＳＮＰの遺伝子型の推定処理を終了する。未選択のノーコールＳＮＰがある場合（ステップＳ１９のＹＥＳ）、推定部５２は、次の対象ＳＮＰを選択する（ステップＳ１７）。
　次に、第２の推定方法について説明する。第２の推定方法では、推定部５２は、信号強度の閾値ｘｌ，ｘｒを用いて、対象ＳＮＰの既知検体の遺伝子型を推定した後、ｋ近傍法を用いて未知検体の遺伝子型を推定する。図５４は、第２の推定方法を示すフローチャートである。図５４のステップＳ１８０１～Ｓ１８０３は、第１の推定方法と同様であるため、説明を省略する。
　ステップＳ１８０５において、推定部５２は、ステップＳ１８０３で取得した閾値ｘｌ，ｘｒに基づいて、対象ＳＮＰの検体のうち、信号強度が最大の検体と、信号強度が最小の検体と、の遺伝子型を推定する。具体的には、推定部５２は、信号強度がｘｌより小さい検体の遺伝子型を遺伝子型２と推定し、信号強度がｘｌ以上ｘｒ以下の検体の遺伝子型を元の遺伝子型と推定し、信号強度がｘｒより大きい検体の遺伝子型を遺伝子型０と推定する。
　図５５は、ステップＳ１８０５における推定の前後のクラスタリングマップの一例を示す図である。図５５の上図は、推定前の遺伝子型を示し、下図は推定後の遺伝子型を示している。図５５の例では、上図に示すように、信号強度が最小の検体（左端の検体）の遺伝子型は遺伝子型２であり、信号強度が最大の検体（右端の検体）の遺伝子型は遺伝子型１である。また、これら２つの検体の信号強度は、いずれも閾値ｘｌより小さい。この場合、下図に示すように、ステップＳ１８０５において、２つの検体の遺伝子型はいずれも遺伝子型２と推定される。これに対して、図５５の右端の検体の信号強度がｘｌ以上ｘｒ以下の場合には、ステップＳ１８０５において、右端の検体の遺伝子型は、遺伝子型１（元の遺伝子型）と推定される。
　ステップＳ１８０６において、推定部５２は、ステップＳ１８０５で推定した２つの検体の遺伝子型が同じか判定する。２つの検体の遺伝子型が同じ場合（ステップＳ１８０６のＹＥＳ）、処理はステップＳ１８０７に進む。
　ステップＳ１８０７において、推定部５２は、全検体の遺伝子型を、ステップＳ１８０５で推定された２つの検体の遺伝子型と同じ遺伝子型と推定する。これは、信号強度が最大及び最小の検体の遺伝子型がいずれも同じ遺伝子型Ｘと推定された場合、以降の処理で、全ての検体の遺伝子型が同じ遺伝子型Ｘと推定されるためである。
　ステップＳ１８０５～Ｓ１８０７において、２つの検体の遺伝子型だけ先に推定し、その推定結果に基づいて全ての検体の遺伝子型を推定することにより、推定装置の計算量を削減することができる。
　図５６は、ステップＳ１８０７における推定の前後のクラスタリングマップの一例を示す図である。図５６の上図は、推定前の遺伝子型を示し、下図は推定後の遺伝子型を示している。図５６の上図は、図５５の下図に相当する。図５６の例では、左端及び右端の検体の遺伝子型が、いずれも遺伝子型２と推定されているため、ステップＳ１８０７において、全ての検体の遺伝子型が遺伝子型２と推定されている。図５６からわかるように、ステップＳ１８０７では、既知検体だけでなく、未知検体の遺伝子型も推定されるため、ｋ近傍法は使用されない。
　その後、推定部５２は、対象ＳＮＰの遺伝子型の推定処理を終了する。未選択のノーコールＳＮＰがある場合（ステップＳ１９のＹＥＳ）、推定部５２は、次の対象ＳＮＰを選択する（ステップＳ１７）。
　一方、ステップＳ１８０５で推定した２つの検体の遺伝子型が異なる場合（ステップＳ１８０６のＮＯ）、処理はステップＳ１８０８に進む。
　ステップＳ１８０８において、推定部５２は、ステップＳ１８０３で取得した閾値ｘｌ，ｘｒに基づいて、対象ＳＮＰの全検体の遺伝子型を推定する。推定方法は、ステップＳ１８０５と同様である。すなわち、推定部５２は、信号強度がｘｌより小さい検体の遺伝子型を遺伝子型２と推定し、信号強度がｘｌ以上ｘｒ以下の検体の遺伝子型を元の遺伝子型と推定し、信号強度がｘｒより大きい検体の遺伝子型を遺伝子型０と推定する。
　図５７は、ステップＳ１８０８における推定の前後のクラスタリングマップの一例を示す図である。図５７の上図は、推定前の遺伝子型を示し、下図は推定後の遺伝子型を示している。図５７の例では、ステップＳ１８０８における推定の結果、信号強度がｘｒより大きい４つの検体の遺伝子型が、遺伝子型０と推定されている。
　ステップＳ１８０９において、推定部５２は、ステップＳ１８０８における推定前の遺伝子型と推定後の遺伝子型との一致率を算出し、算出した一致率が一致率閾値より大きいか判定する。一致率閾値は、任意に設定可能である。図５７の例では、１６個の検体のうち、推定の前後で１２個の検体の遺伝子型が一致しているため、一致率は０．７５となる。
　一致率が一致率閾値より大きい場合（ステップＳ１８０９のＹＥＳ）、処理はステップＳ１８１０に進む。
　ステップＳ１８１０において、推定部５２は、対象ＳＮＰの検体の中に、未知検体があるか判定する。上述の通り、未知検体は、遺伝子型が－１の検体に相当する。このため、推定部５２は、ステップＳ１８０８における推定結果を参照して、遺伝子型－１の検体があるか確認することにより、未知検体があるか判定することができる。
　対象ＳＮＰの検体に未知検体がなかった場合（ステップＳ１８１０のＮＯ）、推定部５２は、対象ＳＮＰの遺伝子型の推定処理を終了する。未選択のノーコールＳＮＰがある場合（ステップＳ１９のＹＥＳ）、推定部５２は、次の対象ＳＮＰを選択する（ステップＳ１７）。
　一方、図５７の例のように、対象ＳＮＰの検体に未知検体があった場合（ステップＳ１８１０のＹＥＳ）、処理はステップＳ１８１１に進む。
　ステップＳ１８１１において、推定部５２は、対象ＳＮＰの検体から未知検体を抽出し、ｋ近傍法により、各未知検体の遺伝子型を推定する。このとき、既知検体の遺伝子型として、ステップＳ１８０８における推定結果が用いられる。ｋ近傍法による遺伝子型の推定方法は、第１実施形態で説明した通りである。推定部５２は、例えば、未知検体に最近傍の（クラスタリングマップ上におけるユークリッド距離が最も近い）ｋ個の検体（又はクラスタ線）を抽出し、抽出した検体（又はクラスタ線）の遺伝子型のうち、最も多い遺伝子型を、未知検体の遺伝子型として推定すればよい。
　その後、推定部５２は、対象ＳＮＰの遺伝子型の推定処理を終了する。未選択のノーコールＳＮＰがある場合（ステップＳ１９のＹＥＳ）、推定部５２は、次の対象ＳＮＰを選択する（ステップＳ１７）。
　これに対して、一致率が一致率閾値以下の場合（ステップＳ１８０９のＮＯ）、処理はステップＳ１８１２に進む。
　ステップＳ１８１２において、推定部５２は、クラスタ毎に多数決法により遺伝子型を推定する。まず、推定部５２は、対象ＳＮＰの既知検体を、ステップＳ１８１１における推定前の遺伝子型毎に分類し、各遺伝子型に対応するクラスタを生成する。各クラスタには、そのクラスタに対応する遺伝子型を有する検体が含まれる。
　次に、推定部５２は、ステップＳ１８０８における推定結果を参照し、各クラスタに含まれる検体の推定後の遺伝子型のうち、最も多い遺伝子型を、そのクラスタの遺伝子型と推定する。そして、推定部５２は、各クラスタに含まれる全ての検体の遺伝子型を、その検体が含まれるクラスタの遺伝子型と推定する。
　図５８は、ステップＳ１８０８，Ｓ１８１２における推定の前後のクラスタリングマップの一例を示す図である。図５８の上図は、ステップＳ１８０８における推定前の遺伝子型を示し、中図は、ステップＳ１８０８における推定後の遺伝子型を示し、下図は、ステップＳ１８１２における推定後の遺伝子型を示している。図５８の上図及び中図は、図５７と対応している。
　図５８の例では、推定部５２は、まず、遺伝子型１の検体を４個含むクラスタ１と、遺伝子型２の検体を１１個含むクラスタ２と、を生成する（上図参照）。
　次に、推定部５２は、ステップＳ１８０８における推定結果を参照し、クラスタ１に含まれる検体の推定後の遺伝子型のうち、最も多い遺伝子型を、クラスタ１の遺伝子型と推定する（中図参照）。図５８の例では、クラスタ１には、遺伝子型０の検体が３個含まれ、遺伝子型１の検体が１個含まれる。したがって、推定部５２は、クラスタ１の遺伝子型を遺伝子型１と推定する。
　そして、推定部５２は、クラスタ１に含まれる全ての検体の遺伝子型を、クラスタ１の遺伝子型である遺伝子型０と推定する（下図参照）。推定部５２は、同様の方法で、クラスタ２に含まれる各検体の遺伝子型も推定する。結果として、クラスタ２に含まれる全ての検体の遺伝子型は、遺伝子型２と推定される。
　その後、処理はステップＳ１８１０に進む。ステップＳ１８１０以降の処理は、上述の通りであり、ステップＳ１８１１において、未知検体の遺伝子型がｋ近傍法により推定される。
　以上説明した通り、本実施形態に係る推定装置は、信号強度の閾値を用いた閾値法により、ノーコールＳＮＰの遺伝子型を推定する。閾値は、遺伝子型を精度よく判定されたフルコールＳＮＰの信号強度データを用いて学習されるため、推定装置は、ノーコールＳＮＰの遺伝子型を精度よく推定することができる。
　また、本実施形態に係る推定方法は、参照データを用いることなく遺伝子型を推定できるため、参照データが十分に得られない場合であっても、利用することができる。
　なお、本実施形態において、ｋ近傍法で使用するパラメータｋの値は、クロス検証により最適に設定されてもよい。クロス検証によるパラメータｋの設定方法は、第１実施形態で説明したとおりである。
　また、本実施形態に係る推定装置は、第１実施形態に係る推定方法を実行可能であってもよい。この場合、推定方法は、ＧＵＩによりユーザが選択可能であるのが好ましい。推定装置は、ユーザにより選択された推定方法を実行すればよい。
（第３実施形態）
　第３実施形態について、図５９～図６１を参照して説明する。本実施形態では、第１実施形態及び第２実施形態に係る推定方法で用いた、ｋ近傍法による遺伝子型の推定方法の変形例について説明する。
　上記の各実施形態では、ｋ近傍法により、未知検体の遺伝子型を推定可能であることが前提であった。しかしながら、複数の未知検体が、クラスタリングマップ上において、相対的に近い位置に集まっている場合、上記のｋ近傍法では、各未知検体の遺伝子型を推定できないことがあり得る。
　図５９は、複数の未知検体が相対的に近い位置に集まったクラスタリングマップの一例を示す図である。図５９において、検体ｓ１～ｓ３は、未知検体であり、互いに近い位置に集まっている。
　図５９の例では、ｋ近傍法のパラメータｋが３である場合、検体ｓ１の最近傍の３個の検体として、遺伝子型が－１の検体ｓ２，ｓ３検体と、遺伝子型が０の検体１つと、が選択される。結果として、多数決アルゴリズムにより、検体ｓ１の遺伝子型は、－１と推定される。すなわち、検体ｓ１の遺伝子型は推定できない。同様の理由で、検体ｓ２，ｓ３の遺伝子型も推定できない。
　本実施形態では、このような場合にも、未知検体の遺伝子型を推定可能なｋ近傍法について説明する。図６０は、本実施形態に係るｋ近傍法を示すフローチャートである。
　ステップＳ２０において、推定部５は、未知検体リストを生成する。未知検体リストとは、対象ＳＮＰの全ての未知検体を含むリストである。
　ステップＳ２１において、推定部５は、ステップＳ２０で生成した未知検体リストが空か判定する。未知検体リストが空の場合（ステップＳ２１のＹＥＳ）、すなわち、対象ＳＮＰに未知検体がない場合、推定部５は、ｋ近傍法による未知検体の遺伝子型の推定処理を終了する。
　一方、未知検体リストが空ではない場合（ステップＳ２１のＮＯ）、すなわち、対象ＳＮＰに未知検体がある場合、処理はステップＳ２２に進む。
　ステップＳ２２において、推定部５は、全ての検体を用いたｋ近傍法により、未知検体リストに含まれる各未知検体の遺伝子型を推定する。ステップＳ２２で実行されるｋ近傍法による遺伝子型の推定は、第１実施形態で説明した通りである。
　ステップＳ２３において、推定部５は、ステップＳ２２で遺伝子型が推定された未知検体を、未知検体リストから削除する。ここで削除された未知検体は、以降の処理では既知検体として扱われる。
　ステップＳ２３で１つ以上の未知検体が未知検体リストから削除された場合（ステップＳ２４のＹＥＳ）、処理はステップＳ２１に戻る。これは、１つ以上の未知検体が新たな既知検体となることにより、ステップＳ２２で遺伝子型を推定できなかった未知検体の遺伝子型が、推定できるようになる可能性があるためである。
　以降、１つ以上の未知検体が未知検体リストから削除されなくなるまで（ステップＳ２２で１つ以上の未知検体の遺伝子型が推定されなくなるまで）、ステップＳ２１～Ｓ２４の処理が繰り返される。
　一方、ステップＳ２３で１つ以上の未知検体が未知検体リストから削除されなかった場合（ステップＳ２４のＮＯ）、処理はステップＳ２５に進む。これは、全検体を用いたｋ近傍法では、未知検体の遺伝子型を推定できなくなった場合に相当する。すなわち、上述の通り、複数の未知検体が相対的に近い位置に集まっていることを意味する。
　ステップＳ２５において、推定部５は、未知検体リストに含まれる未知検体の中から、対象検体を選択する。ここでいう対象検体は、遺伝子型の推定対象となる未知検体のことである。推定部５は、対象検体をランダムに選択してもよいし、以下の方法により選択してもよい。
　まず、推定部５は、未知検体リストに含まれる各未知検体について、最近傍のｋ個の既知検体との間の平均距離を算出する。そして、推定部５は、ｋ個の既知検体との間の平均距離が最小の未知検体を、対象検体として選択する。
　ステップＳ２６において、推定部５は、ステップＳ２５で選択した対象検体の遺伝子型を、既知検体のみを用いたｋ近傍法により推定する。これにより、対象検体の最近傍のｋ個の検体として、ｋ個の既知検体が選択される。したがって、選択された既知検体の遺伝子型に基づいて、対象検体の遺伝子型を推定することができる。
　ステップＳ２７において、推定部５は、対象検体を未知検体リストから削除する。その後、処理はステップＳ２１に戻る。以降、未知検体リストが空になるまで、ステップＳ２１～Ｓ２７の処理が繰り返される。
　図６１は、本実施形態に係るｋ近傍法により、未知検体の遺伝子型が推定される過程に対応する、クラスタリングマップを示す図である。図６１の左上図は、図５９と同様である。説明を簡単にするために、検体ｓ１～ｓ３は、全検体を用いたｋ近傍法では遺伝子型を推定できないものとする。また、検体ｓ１～ｓ３は、この順番で対象検体として選択されるものとする。さらに、パラメータｋは３であるものとする。
　まず、１回目の反復処理において、推定部５は、検体ｓ１を対象検体として選択する（ステップＳ２５）。そして、推定部５は、既知検体を用いたｋ近傍法により、検体ｓ１の遺伝子型を推定する。図６１の例では、検体ｓ１の最近傍の３個の既知検体として、遺伝子型０の３個の既知検体が選択される。したがって、推定部５は、左下図に示すように、検体ｓ１の遺伝子型を遺伝子型０と推定する（ステップＳ２６）。その後、推定部５は、検体ｓ１を未知検体リストから削除する（ステップＳ２７）。以降、検体ｓ１は、遺伝子型０の既知検体となる。
　次に、２回目の反復処理において、推定部５は、検体ｓ２を対象検体として選択する（ステップＳ２５）。そして、推定部５は、既知検体を用いたｋ近傍法により、検体ｓ２の遺伝子型を推定する。図６１の例では、検体ｓ２の最近傍の３個の既知検体として、遺伝子型１の２個の既知検体と、遺伝子型０の１個の既知検体（検体ｓ１）と、が選択される。したがって、推定部５は、右上図に示すように、検体ｓ２の遺伝子型を遺伝子型１と推定する（ステップＳ２６）。その後、推定部５は、検体ｓ２を未知検体リストから削除する（ステップＳ２７）。以降、検体ｓ２は、遺伝子型１の既知検体となる。
　さらに、３回目の反復処理において、推定部５は、検体ｓ３を対象検体として選択する（ステップＳ２５）。そして、推定部５は、既知検体を用いたｋ近傍法により、検体ｓ３の遺伝子型を推定する。図６１の例では、検体ｓ３の最近傍の３個の既知検体として、遺伝子型１の２個の既知検体（検体ｓ２を含む）と、遺伝子型０の１個の既知検体（検体ｓ１）と、が選択される。したがって、推定部５は、右下図に示すように、検体ｓ３の遺伝子型を遺伝子型１と推定する（ステップＳ２６）。その後、推定部５は、検体ｓ３を未知検体リストから削除する（ステップＳ２７）。これにより、検体ｓ３は、遺伝子型１の既知検体となる。
　以上説明した通り、本実施形態に係るｋ近傍法によれば、複数の未知検体が相対的に近い位置に集まっている場合であっても、未知検体の遺伝子型を推定することができる。本実施形態に係るｋ近傍法は、第１実施形態及び第２実施形態のいずれにも適用可能である。
　なお、本発明は上記各実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記各実施形態に開示されている複数の構成要素を適宜組み合わせることによって種々の発明を形成できる。また例えば、各実施形態に示される全構成要素からいくつかの構成要素を削除した構成も考えられる。さらに、異なる実施形態に記載した構成要素を適宜組み合わせてもよい。
１：検体データ記憶部、２：参照データ記憶部、３：クラスタリング強度取得部、４：クラスタリング強度判定部、５：遺伝子型推定部、６：推定結果表示部、５１：閾値学習部、５２：閾値法推定部、１００：コンピュータ、１０１：ＣＰＵ、１０２：入力装置、１０３：表示装置、１０４：通信装置、１０５：記憶装置、Ｇ１：ＳＮＰ選択部、Ｇ２：検体選択部、Ｇ３：ｋ値選択部、Ｇ４：検体一覧表示部、Ｇ５：選択結果表示部、Ｇ６：遺伝子型表示部、Ｇ７：患者選択部、Ｇ８：疾患選択部、Ｇ９：ＳＮＰ情報表示部、Ｇ１０：遺伝子型情報表示部

Claims (38)

  1. 　遺伝子型が未知の未知検体と、前記遺伝子型が既知の既知検体と、を含む複数の検体の遺伝子型データのクラスタリング強度を取得する取得部と、
    　前記クラスタリング強度が第１の閾値より大きい場合、前記遺伝子型データに基づいて、前記未知検体の前記遺伝子型を推定し、推定結果を出力する推定部と、
    を備える遺伝子型推定装置。
  2. 　前記推定部は、前記クラスタリング強度が、前記第１の閾値より小さい第２の閾値より小さい場合、参照ハプロタイプデータに基づいて、前記未知検体の前記遺伝子型を推定する
    請求項１に記載の遺伝子型推定装置。
  3. 　前記推定部は、前記クラスタリング強度が、前記第１の閾値以下かつ前記第２の閾値以上の場合、前記遺伝子型データと前記参照ハプロタイプデータとに基づいて、前記未知検体の前記遺伝子型を推定する
    請求項１又は請求項２に記載の遺伝子型推定装置。
  4. 　前記推定部は、前記遺伝子型データに基づいて、ｋ近傍法により前記未知検体の前記遺伝子型を推定する
    請求項１乃至請求項３のいずれか１項に記載の遺伝子型推定装置。
  5. 　前記推定部は、前記未知検体に最近傍のｋ個のサンプルを選択し、選択した前記サンプルの前記遺伝子型に基づいて、前記未知検体の前記遺伝子型を推定する
    請求項１乃至請求項４のいずれか１項に記載の遺伝子型推定装置。
  6. 　前記サンプルは、前記既知検体又はクラスタ線である
    請求項１乃至請求項５のいずれか１項に記載の遺伝子型推定装置。
  7. 　前記推定部は、前記サンプルの数を投票数として用いる多数決アルゴリズムにより、前記未知検体の前記遺伝子型を推定する
    請求項５又は請求項６に記載の遺伝子型推定装置。
  8. 　前記推定部は、前記サンプルに付与した重みを投票数として用いる多数決アルゴリズムにより、前記未知検体の前記遺伝子型を推定する
    請求項５乃至請求項７のいずれか１項に記載の遺伝子型推定装置。
  9. 　前記推定部は、前記ｋ近傍法のパラメータｋの複数の候補について、クロス検証により推定精度を算出し、前記推定精度が最大の前記候補を前記パラメータｋに設定する
    請求項４乃至請求項８のいずれか１項に記載の遺伝子型推定装置。
  10. 　前記遺伝子型データと、信号強度データと、クラスタリングデータと、の少なくとも１つを含む検体データを記憶する検体データ記憶部を更に備える
    請求項１乃至請求項９のいずれか１項に記載の遺伝子型推定装置。
  11. 　前記信号強度データは、ＤＮＡマイクロアレイ技術により測定された蛍光強度の測定値及び前記測定値に基づいて算出されたパラメータの少なくとも一方を含む
    請求項１０に記載の遺伝子型推定装置。
  12. 　前記推定部は、前記参照ハプロタイプデータに基づいて、インピュテーション法により、前記未知検体の前記遺伝子型を推定する
    請求項２乃至請求項７のいずれか１項に記載の遺伝子型推定装置。
  13. 　連鎖不平衡統計データと、前記参照ハプロタイプデータと、参照遺伝子型頻度データと、の少なくとも１つを含む、参照データを記憶する参照データ記憶部を更に備える
    請求項１乃至請求項１２のいずれか１項に記載の遺伝子型推定装置。
  14. 　前記推定部は、前記連鎖不平衡統計データに含まれる少なくとも１つのスコアに基づいて、前記遺伝子型が未知のＳＮＰに関連する複数の推定用ＳＮＰを抽出し、
    　前記推定用ＳＮＰの前記参照ハプロタイプデータに基づいて、前記未知検体の前記遺伝子型を推定する
    請求項１３に記載の遺伝子型推定装置。
  15. 　前記推定部は、前記スコアが所定の閾値以上の前記ＳＮＰを前記推定用ＳＮＰとして抽出する
    請求項１４に記載の遺伝子型推定装置。
  16. 　前記推定部は、前記推定用ＳＮＰの前記参照ハプロタイプデータの中から、前記未知検体のハプロタイプデータに類似する少なくとも１つの前記参照ハプロタイプデータを選択し、選択した前記参照ハプロタイプデータに基づいて、前記未知検体の前記遺伝子型を推定する
    請求項１４又は請求項１５に記載の遺伝子型推定装置。
  17. 　前記推定部は、前記参照ハプロタイプデータ及び前記ハプロタイプデータに含まれる対立遺伝子を数値に置換し、置換した前記数値に基づいて、前記参照ハプロタイプデータと前記ハプロタイプデータとの距離を算出する
    請求項１６に記載の遺伝子型推定装置。
  18. 　前記取得部は、前記遺伝子型データ及び前記参照遺伝子型頻度データに基づいて、前記クラスタリング強度を算出する
    請求項１３乃至請求項１７のいずれか１項に記載の遺伝子型推定装置。
  19. 　前記ｋ近傍法により選択されたｋ個の前記サンプルと、前記未知検体の遺伝子型の推定結果と、を表示する表示部を更に備える
    請求項５乃至請求項１８のいずれか１項に記載の遺伝子型推定装置。
  20. 　前記遺伝子型データのクラスタリングマップと、ＨａｐＭａｐによる前記遺伝子型の比率と、の少なくとも１つを表示する表示部を更に備える
    請求項１乃至請求項１９のいずれか１項に記載の遺伝子型推定装置。
  21. 　前記未知検体の前記遺伝子型は、前記未知検体に最近傍のｋ個の前記既知検体の前記遺伝子型に基づいて推定される
    請求項１乃至請求項２０のいずれか１項に記載の遺伝子型推定装置。
  22. 　前記未知検体の遺伝子型の推定結果と、前記参照ハプロタイプデータと、前記未知検体の前記ハプロタイプデータと、フェージング後の前記ハプロタイプデータと、インピュテーション後の前記ハプロタイプデータと、ＨａｐＭａｐによる前記遺伝子型の比率と、の少なくとも１つを表示する表示部を更に備える
    請求項１３乃至請求項２１のいずれか１項に記載の遺伝子型推定装置。
  23. 　遺伝子型が未知の未知検体と、前記遺伝子型が既知の既知検体と、を含む複数の検体の遺伝子型データのクラスタリング強度を取得する工程と、
    　前記クラスタリング強度が第１の閾値より大きい場合、前記遺伝子型データに基づいて、前記未知検体の前記遺伝子型を推定し、推定結果を出力する工程と、
    を備える遺伝子型推定方法。
  24. 　遺伝子型が未知の未知検体と、前記遺伝子型が既知の既知検体と、を含む複数の検体の遺伝子型データのクラスタリング強度を取得する工程と、
    　前記クラスタリング強度が第１の閾値より大きい場合、前記遺伝子型データに基づいて、前記未知検体の前記遺伝子型を推定し、推定結果を出力する工程と、
    をコンピュータに実行させるプログラム。
  25. 　遺伝子型が既知の検体である既知検体の、ＤＮＡマイクロアレイ技術により測定された信号強度に基づいて、前記信号強度の閾値により規定される区間と、前記遺伝子型と、の対応関係を学習する閾値学習部と、
    　前記対応関係に基づいて、前記遺伝子型が未知の前記検体である未知検体の前記遺伝子型を推定する閾値法推定部と、
    を備える遺伝子型推定装置。
  26. 　前記閾値学習部は、全ての検体が前記既知検体であるフルコールＳＮＰの前記信号強度に基づいて、前記対応関係を学習する
    請求項２５に記載の遺伝子型推定装置。
  27. 　前記閾値学習部は、前記フルコールＳＮＰの前記信号強度の統計値を、前記遺伝子型毎の算出し、前記統計値に基づいて、前記閾値の候補である閾値候補を生成する
    請求項２６に記載の遺伝子型推定装置。
  28. 　前記閾値学習部は、前記フルコールＳＮＰの前記信号強度の平均値に基づいて、前記区間と前記遺伝子型と、を対応させる
    請求項２６又は請求項２７に記載の遺伝子型推定装置。
  29. 　前記閾値学習部は、前記閾値候補の組合せである閾値組合せを生成し、前記閾値組合せに基づいて前記フルコールＳＮＰの前記遺伝子型を推定し、推定前の前記遺伝子型と、推定後の前記遺伝子型と、の一致率を算出する
    請求項２７又は請求項２８に記載の遺伝子型推定装置。
  30. 　前記閾値学習部は、前記一致率が最大の前記閾値組合せに含まれる前記閾値により規定される前記区間と、前記遺伝子型と、の対応関係を学習する
    請求項２９に記載の遺伝子型推定装置。
  31. 　前記閾値学習部は、前記信号強度の第１の閾値と、前記第１の閾値より大きい第２の閾値と、により規定される３つ前記区間と、前記遺伝子型と、の対応関係を学習する
    請求項２５乃至請求項３０のいずれか１項に記載の遺伝子型推定装置。
  32. 　前記閾値法推定部は、前記対応関係に基づいて、少なくとも１つの前記検体が前記未知検体であるノーコールＳＮＰの前記遺伝子型を推定する
    請求項２５乃至請求項３１のいずれか１項に記載の遺伝子型推定装置。
  33. 　前記閾値法推定部は、前記ノーコールＳＮＰの前記検体の前記遺伝子型を、前記検体の前記信号強度が含まれる前記区間に対応する前記遺伝子型と推定する
    請求項３２に記載の遺伝子型推定装置。
  34. 　前記閾値法推定部は、前記信号強度が所定の前記区間に含まれる前記検体の遺伝子型を、ｋ近傍法により推定する
    請求項２５乃至請求項３３のいずれか１項に記載の遺伝子型推定装置。
  35. 　前記閾値法推定部は、前記信号強度が第１の閾値より小さい前記検体と、前記信号強度が第２の閾値より大きい前記検体と、の前記遺伝子型を、前記対応関係に基づいて推定し、
    　前記信号強度が、前記第１の閾値以上前記第２の閾値以下の前記検体の前記遺伝子型を、ｋ近傍法により推定する
    請求項２５乃至請求項３４のいずれか１項に記載の遺伝子型推定装置。
  36. 　前記未知検体の遺伝子型は、前記未知検体に最近傍のｋ個の前記既知検体の前記遺伝子型に基づいて推定される
    請求項２５乃至請求項３５のいずれか１項に記載の遺伝子型推定装置。
  37. 　遺伝子型が既知の検体である既知検体の、ＤＮＡマイクロアレイ技術により測定された信号強度に基づいて、前記信号強度の閾値により規定される区間と、前記遺伝子型と、の対応関係を学習する工程と、
    　前記対応関係に基づいて、前記遺伝子型が未知の前記検体である未知検体の前記遺伝子型を推定する工程と、
    を備える遺伝子型推定方法。
  38. 　遺伝子型が既知の検体である既知検体の、ＤＮＡマイクロアレイ技術により測定された信号強度に基づいて、前記信号強度の閾値により規定される区間と、前記遺伝子型と、の対応関係を学習する工程と、
    　前記対応関係に基づいて、前記遺伝子型が未知の前記検体である未知検体の前記遺伝子型を推定する工程と、
    をコンピュータに実行させるプログラム。
PCT/JP2015/080573 2014-10-30 2015-10-29 遺伝子型推定装置、方法、及びプログラム WO2016068245A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2016556630A JP6280657B2 (ja) 2014-10-30 2015-10-29 遺伝子型推定装置、方法、及びプログラム
US15/419,903 US11355219B2 (en) 2014-10-30 2017-01-30 Genotype estimation device, method, and program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2014-221614 2014-10-30
JP2014221614 2014-10-30

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US15/419,903 Continuation US11355219B2 (en) 2014-10-30 2017-01-30 Genotype estimation device, method, and program

Publications (1)

Publication Number Publication Date
WO2016068245A1 true WO2016068245A1 (ja) 2016-05-06

Family

ID=55857580

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2015/080573 WO2016068245A1 (ja) 2014-10-30 2015-10-29 遺伝子型推定装置、方法、及びプログラム

Country Status (4)

Country Link
US (1) US11355219B2 (ja)
JP (1) JP6280657B2 (ja)
TW (1) TWI584143B (ja)
WO (1) WO2016068245A1 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10777302B2 (en) * 2012-06-04 2020-09-15 23Andme, Inc. Identifying variants of interest by imputation

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002025489A1 (fr) * 2000-09-19 2002-03-28 Hitachi Software Engineering Co., Ltd. Technique d'affichage de donnees genetiques et support d'enregistrement a cet effet
JP2007334719A (ja) * 2006-06-16 2007-12-27 Hitachi Software Eng Co Ltd 遺伝子発現解析の欠損値補完システム
JP2010086142A (ja) * 2008-09-30 2010-04-15 National Institute Of Agrobiological Sciences 遺伝子クラスタリング装置およびプログラム
JP2010218150A (ja) * 2009-03-16 2010-09-30 Meiji Univ 発現プロファイル解析システム及びそのプログラム
US20140107933A1 (en) * 2006-11-30 2014-04-17 The Johns Hopkins University Gene expression barcode for normal and diseased tissue classification

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4307807B2 (ja) 2002-08-13 2009-08-05 株式会社エヌ・ティ・ティ・データ データ処理装置、方法及びコンピュータプログラム
US20060271300A1 (en) * 2003-07-30 2006-11-30 Welsh William J Systems and methods for microarray data analysis
US7035740B2 (en) 2004-03-24 2006-04-25 Illumina, Inc. Artificial intelligence and global normalization methods for genotyping
JP4777631B2 (ja) * 2004-09-27 2011-09-21 株式会社日立ハイテクノロジーズ 核酸増幅分析法および装置
JP2006107396A (ja) * 2004-10-08 2006-04-20 Institute Of Physical & Chemical Research Ｓｎｐ遺伝子型分類方法、ｓｎｐ遺伝子型分類装置およびｓｎｐ遺伝子型分類プログラム
US8200440B2 (en) 2007-05-18 2012-06-12 Affymetrix, Inc. System, method, and computer software product for genotype determination using probe array data
WO2009017204A1 (ja) 2007-08-01 2009-02-05 Olympus Corporation クラスタリング方法、プログラムおよび装置
US8429153B2 (en) * 2010-06-25 2013-04-23 The United States Of America As Represented By The Secretary Of The Army Method and apparatus for classifying known specimens and media using spectral properties and identifying unknown specimens and media
AU2011358564B9 (en) * 2011-02-09 2017-07-13 Natera, Inc Methods for non-invasive prenatal ploidy calling
TWI582425B (zh) * 2011-11-25 2017-05-11 凸版印刷股份有限公司 樣本分析晶片、樣本分析方法及基因解析方法
HUP1200622A2 (en) 2012-10-30 2014-05-28 Budapesti Mueszaki Es Gazdasagtudomanyi Egyetem Method and computer program product for genotype classification
CN104995314A (zh) * 2012-11-26 2015-10-21 财团法人国家卫生研究院 用于基因型鉴定结核分枝杆菌的引物、snp标记及方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002025489A1 (fr) * 2000-09-19 2002-03-28 Hitachi Software Engineering Co., Ltd. Technique d'affichage de donnees genetiques et support d'enregistrement a cet effet
JP2007334719A (ja) * 2006-06-16 2007-12-27 Hitachi Software Eng Co Ltd 遺伝子発現解析の欠損値補完システム
US20140107933A1 (en) * 2006-11-30 2014-04-17 The Johns Hopkins University Gene expression barcode for normal and diseased tissue classification
JP2010086142A (ja) * 2008-09-30 2010-04-15 National Institute Of Agrobiological Sciences 遺伝子クラスタリング装置およびプログラム
JP2010218150A (ja) * 2009-03-16 2010-09-30 Meiji Univ 発現プロファイル解析システム及びそのプログラム

Also Published As

Publication number Publication date
JPWO2016068245A1 (ja) 2017-06-08
US11355219B2 (en) 2022-06-07
TW201633195A (zh) 2016-09-16
TWI584143B (zh) 2017-05-21
JP6280657B2 (ja) 2018-02-14
US20170364631A1 (en) 2017-12-21

Similar Documents

Publication Publication Date Title
Stein-O’Brien et al. Decomposing cell identity for transfer learning across cellular measurements, platforms, tissues, and species
JP6854272B2 (ja) 遺伝子の変異の非侵襲的な評価のための方法および処理
JP7302081B2 (ja) 深層ニューラルネットワークに基づくバリアント分類器
Stranger et al. Patterns of cis regulatory variation in diverse human populations
Jordan et al. Human allelic variation: perspective from protein function, structure, and evolution
KR101542529B1 (ko) 대립유전자의 바이오마커 발굴방법
US8600718B1 (en) Computer systems and methods for identifying conserved cellular constituent clusters across datasets
US20130212125A1 (en) Bioinformatics search tool system for retrieving and summarizing genotypic and phenotypic data for diagnosing patients
JP2005527904A (ja) 複雑性疾患を構成疾患に細分するコンピュータ・システムおよび方法
CN116640847A (zh) 癌症进化检测和诊断
JP2003021630A (ja) 臨床診断サービスを提供するための方法
WO2003065282A1 (en) Computer systems and methods for identifying genes and determining pathways associated with traits
JP7041614B2 (ja) 生体データにおけるパターン認識のマルチレベルアーキテクチャ
KR20140061223A (ko) 차세대 시퀀싱 데이터의 질병변이마커 검출 방법
KR20220136462A (ko) 서열-특정 오류(sse)를 유발시키는 서열 패턴을 식별하기 위한 심층 학습-기반 프레임워크
JP2019514148A (ja) デジタルデータを解析するための方法
KR20150024232A (ko) 질병에 대한 약물 내성 유전체로부터 내성 원인 마커의 발굴 방법
JP6280657B2 (ja) 遺伝子型推定装置、方法、及びプログラム
Grieve et al. Genome-wide co-expression analysis in multiple tissues
Gu et al. Analysis of allele specific expression-A survey
Che et al. A meta-analysis strategy for gene prioritization using gene expression, SNP genotype, and eQTL data
WO2020109356A1 (en) Vector-based haplotype identification
Wang et al. Gene-Focused Networks Underlying Phenotypic Convergence in a Systematically Phenotyped Cohort With Heterogeneous Intellectual Disability
Ahuja et al. A Study and Analysis of Disease Identification using Genomic Sequence Processing Models: An Empirical Review
Althagafi Visualization and Simulation of Variants in Personal Genomes With an Application to Premarital Testing (VSIM)

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2016556630

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 15855363

Country of ref document: EP

Kind code of ref document: A1