WO2006059770A1 - 確率算出方法 - Google Patents

確率算出方法 Download PDF

Info

Publication number
WO2006059770A1
WO2006059770A1 PCT/JP2005/022317 JP2005022317W WO2006059770A1 WO 2006059770 A1 WO2006059770 A1 WO 2006059770A1 JP 2005022317 W JP2005022317 W JP 2005022317W WO 2006059770 A1 WO2006059770 A1 WO 2006059770A1
Authority
WO
WIPO (PCT)
Prior art keywords
sample
haplotype
frequency
diplotype
locus
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
PCT/JP2005/022317
Other languages
English (en)
French (fr)
Inventor
Naoyuki Kamatani
Toshimasa Yamazaki
Masao Yanagisawa
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
RIKEN
Original Assignee
RIKEN
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 RIKEN filed Critical RIKEN
Priority to JP2006546685A priority Critical patent/JP4755601B2/ja
Publication of WO2006059770A1 publication Critical patent/WO2006059770A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/40—Population genetics; Linkage disequilibrium

Definitions

  • the present invention relates to a probability calculation method for calculating the probability (significance level) of the first type of error in an independence test using case-contro 1 correlation analysis using SNP (Single Nucleotide Polymorphism). .
  • case-control correlation analysis using SNP as risk factors, (1) allele frequency model (AZa), (2) dominant genetic model (AA + Aa / aa), (3) recessive genetic model (AA / Aa + aa), (4) Genotype model (AAZAaZaa) is assumed.
  • A represents a major allele
  • a represents a minor allele.
  • a 2 X 2 or 2 X 3 contingency table is created and an independence test is performed.
  • Bonferroni's correction method the value obtained by multiplying the p-value obtained by the independence test by the number of independence tests (4 times in the case of the four models described above) is regarded as the multiple test. This excludes the effects of multiple testing. This method interprets the test results in the most conservative (strict) way.
  • the method described in Non-Patent Document 1 is known.
  • the linkage disequilibrium index value ( ⁇ ) between each two SNPs is calculated for a small number of SNPs, and approximation is performed by performing eigenvalue decomposition of the matrix containing the index values as components. Corrections are made.
  • the method described in Non-Patent Document 1 is a method for correcting the p-value in a multiple test focusing on linkage disequilibrium.
  • Patent Document 1 Dale R. Nyholt, A simple correction for multiple test ing for single- nucleotide polymorphisms in linkage disequilibrium with each other, Am. J. Hum. Genet., Vol. 74, pp. 765- 769 Disclosure of the Invention
  • the Bonferroni correction method has a problem that there is a risk of overlooking a disease-related SNP because the corrected p-value is likely to be larger than the true p-value.
  • the p-value of the entire multiple test is set to a commonly used value such as 0.01 or 0.05, in order to allow the set p-value. Since the p-value of each SNP is very small, that is, it is a very severe condition, there is also a problem that there is a risk of overlooking the disease-related SNP as well. For example, if 50000 SNP data are targeted and the p-value of the entire multiple test is set to 0.01, the p-value of each SNP is always 0.0000002 (0.0.1 + 50000) / J, Value.
  • Non-Patent Document 1 has a problem that haplotype frequency is not taken into account because it focuses on pairwise linkage disequilibrium. Furthermore, there was a problem that the direct p-value was not calculated.
  • the present invention has been made in view of the above problems.
  • the P value in the independence test performed by case-control correlation analysis using SNP that is, the probability of the first type of error is considered, and the haplotype frequency is considered. It is an object of the present invention to provide a probability calculation method that can be directly calculated and, as a result, can prevent oversight of disease-related SNPs.
  • the probability calculation method calculates the probability of the first type error in the independence test performed by case-control correlation analysis using SNP.
  • a sample acquisition step for acquiring a first sample and a second sample including a plurality of haplotype information related to haplotypes composed of a plurality of SNPs, and a first sample acquired in the sample acquisition step.
  • a haplotype frequency calculation step for calculating the haplotype frequency for each haplotype that can exist, and each locus of each haplotype is a minor allele.
  • a matrix definition step for defining a matrix having the predetermined value as a component, a haplotype frequency of each haplotype in the first sample, and a haplotype frequency of each haplotype in the second sample are set.
  • a haplotype frequency setting step that gives a combination of each haplotype frequency of the first sample and each haplotype frequency of the second sample, a plurality of haplotype frequencies corresponding to the first sample set in the haplotype frequency setting step, and the matrix definition step
  • the allele frequency at each locus in the first sample is calculated based on the matrix defined in, and the multiple haplotype frequencies corresponding to the set second sample and the defined matrix at each locus in the second sample are calculated.
  • a contingency table creation step for creating a contingency table corresponding to each locus based on the allele frequency at each locus in the first sample and the allele frequency at each locus in the second sample; Either the independence test execution step for executing the independence test for each locus or the test result of the independence test executed in the independence test execution step based on the contingency table created in the generation step is If significant, based on the number of samples in the first sample, the number of samples in the second sample, the total number of haplotypes that may be present, the haplotype frequency in the first sample, the haplotype frequency in the second sample, and the no-protype frequency The probability corresponding to the combination of each haplotype frequency of the first sample and each haplotype frequency of the second sample set in the haplotype frequency setting step is a haplo of the first sample and the second sample set in advance.
  • the probability calculation method according to claim 2 is a probability calculation method for calculating the probability of the first type of error in the independence test performed in case-contro relation analysis using SNP.
  • the contingency table creation step for creating the contingency table corresponding to each locus and the contingency table created in the contingency table creation step for each locus
  • one of the independence test execution step that executes the independence test and the test result of the independence test executed in the independence test execution step is significant
  • the probability corresponding to the combination of each diplotype form frequency of the first sample set in the diplotype form frequency setting step and each diplotype form frequency of the second sample is set between the first sample and the second sample set in advance.
  • the probability calculation method according to claim 3, which is useful for the present invention, is a probability calculation method for calculating the probability of the first type of error in the independence test performed in case-contro relation analysis using SNP.
  • the haplotype frequency calculation step for calculating the haplotype frequency for each haplotype that can exist, and whether each locus of each haplotype has a minor allele or major allele is predetermined.
  • a matrix definition step for defining a matrix having the predetermined value as a component and the first sample
  • the diplotype form frequency setting step that gives the combination, the plurality of diplotype form frequencies corresponding to the first sample set in the diplotype form frequency setting step, and the matrix defined in the matrix definition step
  • Genotypes that calculate genotype frequencies at each locus and calculate genotype frequencies at each locus in the second sample based on multiple diplotype frequencies corresponding to the set second sample and a defined matrix The genotype frequency at each locus in the first sample calculated in the frequency calculation step and the genotype frequency calculation step.
  • the contingency table creation step for creating the contingency table corresponding to each locus, and the contingency table created in the contingency table creation step ,
  • the independence test execution step that executes the independence test for each locus and the test result of the independence test executed in the independence test execution step!
  • the first sample Based on the number of samples, the number of samples in the second sample, the total number of possible haplotypes, the number of diplotypes in the first sample, the number of diplotypes in the second sample, and the frequency of the diplotype, Probability corresponding to the combination of each diplotype power of the first sample and each diplotype power of the second sample set in the setting step is related to the diplotype power of the first and second samples set in advance.
  • It includes a probability calculation step of calculating using the interleaf distribution, and the diplotype power setting step, the genotype of The number calculation step, the contingency table creation step, the independence test execution step, and the probability calculation step are performed for all combinations of each diplotype form frequency of the first sample and each of the divotype form frequencies of the second sample.
  • the probability of the first type error is finally calculated by adding the probabilities calculated each time repeatedly.
  • the probability calculation method according to claim 4 which is useful for the present invention, approximates the probability of the first type of error in the independence test performed by case-control correlation analysis using SNP using the MCMC method.
  • the first haplotype selection step for selecting the first haplotype and the haplotype power of the first haplotype in the sample selected in the sample selection step
  • the second haplotype selection step for selecting a second haplotype different from the first haplotype, and the haplotype of the first haplotype based on the preset haplotype frequency of the first haplotype and the haplotype frequency of the second haplotype.
  • a haplotype frequency update step for updating the haplotype frequency and the haplotype frequency of the second haplotype, a plurality of haplotype frequencies corresponding to the first sample updated in the haplotype frequency update step, and a matrix defined in the matrix definition step.
  • An allele frequency calculating step to calculate, and the allele Create a contingency table for each locus based on the allele frequency at each locus in the first sample and the allele frequency at each locus in the second sample calculated in the number calculation step For the step and the contingency table created in the contingency table creation step, the independence test execution step for performing the independence test for each locus and the independence test test result performed in the independence test execution step.
  • a cumulative number updating step of updating the cumulative number of times in the case where any one of the two is significant the sample selection step, the first haptic A mouth type selection step, the second haplotype selection step, the haplotype frequency update step, the allele frequency calculation step, the contingency table creation step, the independence detection execution step, and the cumulative number update step are repeated a predetermined number of times, The ratio of the final cumulative number of times to the predetermined number of times is calculated as the first type error probability.
  • the probability calculation method includes: (1) acquiring a first sample and a second sample including a plurality of haplotype information related to a haplotype composed of a plurality of SNPs; and (2) acquiring the acquired first Based on the haplotype information contained in the sample and the haplotype information contained in the acquired second sample, the haplotype frequency is calculated for each haplotype that can exist. (3) Each locus of each haplotype is a minor allele (with the same locus). By setting a predetermined value according to whether it has a less frequent allele or a major allele (the more frequent allele at the same sitting position), a matrix having the predetermined value as a component is defined.
  • each haplotype of the first sample is set.
  • the haplotype frequencies of the second sample are combined, and (5) the haplotype frequencies corresponding to the set first sample and the matrix frequencies defined in the first sample based on the defined matrix and the defined matrix.
  • the contingency table corresponding to each locus (note that the contingency table is contingency t able, translated into contingency table in Japanese) (7) Independence test is performed for each locus on the generated contingency table.
  • Sample number of the first sample, second sample The number of haplotypes in the first sample and the number of haplotypes in the second sample, based on the number of haplotypes in the first sample, the haplotype frequency in the first sample, the haplotype frequency in the second sample, and the frequency of the haplotypes The probability corresponding to the combination of and is calculated using a preset distribution of the haplotype frequencies of the first and second samples. And (4), ( 5), (6), (7), and (8) are repeated for all combinations of each haplotype frequency of the first sample and each haplotype frequency of the second sample. By adding the probabilities, the probability of type 1 error is finally calculated. As a result, it is possible to directly calculate the probability of the first type error in consideration of the haplotype frequency, and as a result, it is possible to prevent the oversight of the disease-related SNP.
  • the probability calculation method includes: (1) acquiring a first sample and a second sample including a plurality of haplotype information related to a haplotype composed of a plurality of SNPs, and (2) acquiring the acquired first Based on the haplotype information contained in the specimen and the haplotype information contained in the acquired second specimen, the haplotype frequency is calculated for each haplotype that can exist, and (3) each locus of each haplotype is identified as a minor or major allele.
  • a matrix having the predetermined value as a component is defined.
  • the diplotype frequency of each diplotype in the first sample and each matrix in the second sample By setting the diplotype power of the diplotype, the set of each diplotype power of the first sample and each diplotype power of the second sample (5) Calculate the allele frequency at each locus in the first sample based on the defined diplotype frequencies corresponding to the set first sample and the defined matrix, and correspond to the set second sample Calculate the allele frequency at each locus in the second sample based on the multiple diplotype frequencies and the defined matrix, and (6) the calculated allele frequency at each locus in the first sample and the second sample Based on the allele frequency at each sitting position, create an even table corresponding to each sitting position. (7) Perform an independence test for each sitting position.
  • the number of samples in the first sample, the number of samples in the second sample, the total number of haplotypes that can exist, the diplotype form frequency in the first sample, the number of samples in the second sample Diplotype power and Noplota Based on the frequency, the probability corresponding to the combination of each diplotype form frequency of the first sample and each diplotype form frequency of the second sample is determined between the preset first sample and second sample. Calculated using the bond distribution for the diplotype form frequency.
  • the probability calculation method includes (1) acquiring a first sample and a second sample including a plurality of haplotype information related to haplotypes configured by a plurality of SNPs, and (2) acquiring Based on the haplotype information included in the first sample and the haplotype information included in the acquired second sample, the haplotype frequency is calculated for each haplotype that can exist.
  • Each locus of each haplotype is a minor allele or major By setting a predetermined value according to whether alleles are present, a matrix with the predetermined value as a component is defined, and (4) the diplotype frequency of each diplotype in the first sample and the second sample By setting the diplotype power of each diplotype in, each diplotype power of the first sample and each diplotype power of the second sample (5) Based on the multiple diplotype frequencies corresponding to the set first sample and the defined matrix, genotype frequencies at each locus in the first sample are calculated, and the set second sample is calculated. The genotype frequency at each locus in the second sample is calculated based on the corresponding multiple diplotype frequencies and the defined matrix.
  • the genotype frequency and the Based on the genotype frequency at each locus in 2 samples create an even table corresponding to each locus, (7) Run an independent test for each locus on the created even table, (8) If one of the results of the independence test performed is significant, the number of samples in the first sample, the number of samples in the second sample, the total number of haplotypes that can exist, the diplotype form factor of the first sample, The diplotype form factor of the second sample.
  • the probability corresponding to the combination of each diplotype form frequency of the set first sample and each diplotype form frequency of the second sample is determined based on the frequency of the first and second samples. Calculated using the joint distribution for the diplotype frequency with the second sample.
  • the probability calculation method includes (1) setting a predetermined value that is predetermined according to whether each locus of each haplotype has a minor allele or a major allele. Define a matrix whose components are predetermined values, (2) generate a combination of haplotype frequencies for the first and second samples set in advance, and (3) either the first sample or the second sample (4) Select the first haplotype, (5) If the haplotype frequency of the first haplotype in the selected sample is not 0, select a second haplotype that is different from the first haplotype, (6 ) Based on the haplotype frequency of the first haplotype and the haplotype frequency of the second haplotype set in advance, the haplotype frequency of the first haplotype and the haplotype frequency of the second haplotype (7) Calculate the allele frequencies at each locus in the first sample based on multiple haplotype frequencies corresponding to the updated first sample and the defined matrix, and update the second sample The allele frequencies at each locus in the second sample are
  • FIG. 1 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type of error using an allele frequency model.
  • Figure 2 shows the probability calculation that accurately calculates the probability of type 1 error in the dominant / recessive genetic model. It is a flowchart which shows an example of the exit method.
  • FIG. 3 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type error in the genotype model.
  • FIG. 4 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type of error when an accurate haplotype frequency is not known.
  • FIG. 5 is a flowchart showing an example of a probability calculation method for approximately calculating the probability of the first type of error using the MCMC method in the allele frequency model.
  • FIG. 6 is a flowchart showing an example of a probability calculation method for calculating the probability of the first type of error using the MCMC method in a dominant / recessive genetic model.
  • FIG. 7 is a flowchart showing an example of a probability calculation method for approximately calculating the probability of the first type of error using the MCMC method in the genotype model.
  • FIG. 8 is a diagram illustrating an example of information included in the first sample and the second sample.
  • FIG. 9 is a diagram showing an example of information on all haplotypes that can exist.
  • FIG. 10 is a diagram showing an example of haplotype frequencies for all haplotypes that can exist.
  • FIG. 11 is a diagram illustrating an example of a defined matrix.
  • FIG. 12 is a diagram showing an example of information related to haplotype frequency.
  • FIG. 13 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • FIG. 14 is a diagram showing the relationship between diplotype model numbers and haplotype numbers.
  • FIG. 15 is a diagram showing an example of information related to diplotype form frequencies.
  • FIG. 16 is a diagram showing an example of a contingency table created in the case of a recessive genetic model.
  • FIG. 17 is a diagram showing an example of a contingency table created in the case of a dominant genetic model.
  • FIG. 18 is a diagram showing an example of a contingency table created in the case of a genotype model.
  • FIG. 19 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • FIG. 20 is a diagram showing an example of an even table created in the case of the allele frequency model.
  • FIG. 21 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • FIG. 22 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • FIG. 23 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • FIG. 24 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • FIG. 25 is a diagram showing an example of an even table created in the case of an allele frequency model.
  • the probability calculation method for the present invention which accurately calculates the probability of the first type of error, is as follows: (1 1) allele frequency model, (1 2) dominant. Recessive genetic model, (1 3) genotype The model, (1-4) If the exact haplotype frequency is weak, explain in detail with reference to the figure.
  • FIG. 1 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type of error using the allele frequency model.
  • a first sample and a second sample containing a plurality of haplotype information related to a haplotype composed of a plurality of SNPs are acquired (step SA-1). Specifically, as shown in Fig. 8, N pieces of haplotype information consisting of L (L is a positive integer) chained SNPs are included.
  • the haplotype frequency is calculated for each haplotype that may exist (Ste SA—2).
  • the total number of haplotypes that can exist is specifically 2 or less and M may be described as shown in FIG. is there.
  • the haplotype frequency h (where j is a number that identifies the haplotype and satisfies l ⁇ j ⁇ M)
  • the haplotype frequency h satisfies the following formula 1.
  • a predetermined value is set according to the force at which each locus of each haplotype has a minor allele or major allele, thereby defining a matrix having the predetermined value as a component (step SA-3). Specifically, if the i-th position in the haplotype of the haplotype number j (i is an integer that identifies the locus and satisfies l ⁇ i ⁇ L) is the minor allele (a), 1 is set as the major allele. In the case of (A), 0 is set as a predetermined value.
  • the matrix to be defined can be expressed as shown in Fig. 11, and a value of 0 or 1 for each combination of the locus number i (i is an integer satisfying l ⁇ i ⁇ L) and the haplotype number j. Is stored.
  • the haplotype frequency X of haplotype number j is a random variable.
  • the actual value X of the haplotype frequency X of the haplotype number j in the first sample is set for all haplotype numbers j, and the second
  • Haplotype frequency of sample with haplotype number j Realized value of X
  • step SA-4 based on the multiple haplotype frequencies corresponding to the first sample set in step SA-4 and the matrix defined in step SA-3, the allele frequency at each locus in the first sample is calculated.
  • the key at each locus in the second sample is Calculate the real frequency (step SA—5).
  • the minor allele at the locus number i in the first sample and the minor allele at the locus number i in the second sample is calculated.
  • the allele frequency Y is a random variable, and is expressed by the following Equation 3. And below
  • ⁇ is the matrix defined in step SA-3.
  • Step SA-5 based on the allele frequency at each locus in the first sample and the allele frequency at each locus in the second sample calculated in step SA-5, the occurrence corresponding to each locus.
  • Step SA—6 Specifically, 2 x 2 contingency tables shown in Fig. 13 are created for the number of sitting positions (L).
  • step SA-7 an independence test is performed for each locus on the contingency table created in step SA-6 (step SA-7).
  • the independence test will be briefly described.
  • X defined by Equation 4 below is calculated for an even table, and the calculated X value has a degree of freedom (1, 1) and a significance level (for example, 0.01 or 0.05).
  • Equation 4 i is a number that identifies the specimen.
  • j is a number for identifying the allele.
  • rij is the allele frequency.
  • n is the minor allele frequency (n) and measure for sample number i
  • n is the allele frequency of the case group corresponding to allele number j (
  • n is the number of cases and lj 2j
  • n ⁇ 2 ⁇ 2 n .
  • step SA—8 Yes
  • the number of samples in the first sample, the second sample Sample haplotypes, the total number of haplotypes that can be present, the haplotype frequency of the first sample, the haplotype frequency of the second sample, and the frequency of the negative protypes, and then each haplo of the first sample set in step SA-4.
  • the probability corresponding to the combination of the type frequency and each haplotype frequency of the second sample is calculated using the preset distribution of the haplotype frequencies of the first and second samples (step SA-9). Specifically, the number of samples of the first sample, N, the number of samples of the second sample, N,
  • Equation 5 k is a number that identifies the sample.
  • step SA—10 Yes
  • step SA—4 to Step SA—9
  • All the probabilities that have been issued are added, and the added probability is finally set as the first type error probability (step SA-11).
  • step SA-10 No
  • step SA-4 executes step SA-4 to step SA-9 and repeat until all combinations are completed.
  • the probability calculation method (1) obtains a first sample and a second sample including a plurality of haplotype information related to a haplotype composed of a plurality of SNPs, and (2) obtains Based on the haplotype information contained in the first sample and the haplotype information contained in the acquired second sample, calculate the haplotype frequency for each haplotype that may exist. (3) Each locus of each haplotype is minor By setting a predetermined value according to whether the allele or major allele is present, a matrix with the predetermined value as a component is defined.
  • the haplotype frequency of each haplotype in the first sample and the second sample By setting the haplotype frequency of each haplotype, a combination of each haplotype frequency of the first sample and each haplotype frequency of the second sample is given, (5) The allele frequencies at each locus in the first sample are calculated based on the defined haplotype frequencies corresponding to the first sample and the defined matrix, and the haplotype frequencies and the corresponding haplotype frequencies corresponding to the set second sample are calculated. Based on the defined matrix, the allele frequency at each locus in the second sample is calculated.
  • step SA-5 the allele frequency Y of the minor allele at the locus number i in the first sample shown in Equation 3 and the minor allele at the locus number i in the second sample are shown.
  • the allele frequency Y of haplota is related to the allele frequency p of the minor allele at locus number i.
  • step SA-9 the combined distribution related to the haplotype frequency between the first sample and the second sample shown in Equation 5 is the combined distribution related to the haplotype frequency of the first sample shown in Equation 7 below and Equation 8 below. It is given by the product of the second sample shown and the joint distribution of haplotype frequencies.
  • two samples are obtained from the same population as the null hypothesis (H).
  • haplotype frequencies of haplotype number j in these two samples are random variables X and X, respectively, the two samplings are
  • the 2j probabilities are expressed by the following formula 7 and the following formula 8, respectively. Note that both samples are assumed to be from a population with the same haplotype frequency (s) (null hypothesis).
  • X and X lj 2j are real values of the random variables X and X, respectively.
  • the difference between the concept of haplotype copy and haplotype in the experiment of extracting haplotypes is as follows. For example, if one individual has two identical wild-types (homozygotes), this individual is considered to possess one haplotype and two haplotype copies. Also, the first species calculated in step SA-11 The probability of error can be expressed by Equation 9 below. The probability calculated by Equation 9 below is the probability under the null hypothesis, and means the probability that any one of the test results of the independence test performed in Step SA-7 will be significant.
  • Equation 9 Z is 1 if one of the test results of the independence test performed in Step SA-7 is significant, and 0 if it is not. It is a random variable defined to take The variable Z is a function of the random variable Y described above, and
  • the random variable Y is a function of the random variable X described above
  • the variable Z is a function f (x, ⁇ ⁇ ki kj 11
  • ⁇ a ⁇ is an s XL matrix.
  • the function Is 21 22) is defined as follows. ⁇ dl,
  • b is the proportion of js + l of s-Lth haplotype haplotype copies of the jth sample that is a minor allele at the kth position. This is a value that varies depending on the composition of the X actual haplotypes of the s + 1st haplotype of the result, and the run js + l
  • the independence test is performed for each sitting position on the created contingency table. Then define a function f that is 1 if any of the results of the independence test is significant, and 0 otherwise.
  • a and b are non-negative integers, and m ⁇ a and n ⁇ b.
  • Equation 12 is negative when a ⁇ bmZn and positive when a> bmZn. It should be noted that m ⁇ a and n ⁇ b in Equation 12 and Equation 13 are not 0 at the same time. Therefore, T is monotonically decreasing at a bm / n and monotonically increasing at a ⁇ bmZn. When b is fixed, in the range of a ⁇ a ⁇ a
  • Equation 13 is negative when b ⁇ anZm and positive when b> anZm.
  • T is monotonically decreasing when b ⁇ anZm, and monotonically increasing when b> anZm. Therefore, T is b ⁇ b ⁇ b when a is fixed
  • Equation 11 is expressed as (a 2 a, b 2 b), (a
  • the conservative probability calculation method for the major haplotypes is described in detail again.
  • the allele at the k-th position of the s + 1st haplotype is all minor alleles or all in the first and second samples.
  • the test is performed as a major allele, and if any test is significant, the k-position is determined to be significant.
  • s (the cumulative frequency of major haplotypes) should be as large as possible. Then the confidence in the estimated frequency of major haplotypes will be less reliable.
  • the reliability of the haplotype frequency at the s + l, ..., L locus may be low.
  • the minor allele frequency at the kth locus in a population is often known. If the minor allele frequency at the k locus in this population is q, then k
  • the minor allele frequency of the kth locus in other haplotypes can be considered as q— ⁇ s «h.
  • the number of copies of other haplotypes in the jth sample is X
  • the calculation method may be a method such as linear regression or curve regression.
  • FIG. 2 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type of error in the dominant / recessive genetic model.
  • a first sample and a second sample containing a plurality of haplotype information related to a haplotype composed of a plurality of SNPs are acquired (step SB-1). Note that the explanation regarding the acquisition of the first specimen and the second specimen is the same as the above-mentioned (11) allele frequency model.
  • the haplotype frequency is calculated for each haplotype that may exist. (Step SB—2).
  • the haplotype frequency is the same as the (1-1) allele frequency model described above.
  • a predetermined value is set according to the force at which each locus of each haplotype has a minor allele or major allele, thereby defining a matrix having the predetermined value as a component (step SB-3).
  • the explanation about the definition of the matrix is the same as the above (11) array frequency model.
  • each diplotype power and the first diplotype power in the second sample are set.
  • each diplotype form factor of the two samples step SB—4.
  • Mouth type form frequency X is a random variable, and satisfies the following formula 14 respectively.
  • the actual value X of the diplopivo frequency of diplotype model number jj 'in the first sample is set for all diplotype model numbers jj', and the second sample
  • the allele frequencies at each locus in the first sample And calculate the allele frequencies at each position in the second sample based on the multiple probabilistic frequencies corresponding to the second sample set in step SB-4 and the matrix defined in step SB-3. (Step SB—5).
  • the allele frequency Y of the minor allele at sitting position i in the first sample and the minor allele at sitting position i in the second sample is calculated.
  • the allele frequency Y is a random variable, and is represented by the following formula 15. And below Based on Formula 15, the actual value y of the minor allele allele frequency at locus i in the first sample and the minor allele frequency y at locus i in the second sample y li
  • the allele frequency of the minor allele (aa + Aa) is calculated using Formula 15, but in the case of the dominant genetic model, the minor allele (aa) is calculated using Formula 16 below. Calculate the allele frequency.
  • the contingency table corresponding to each locus is calculated. Create (Step SB—6). Specifically, in the case of the recessive genetic model, the 2 x 2 contingency table shown in Figure 16 for the allele frequency calculated using Equation 15 is created for the number of loci (L), and the dominant genetic model Then, the 2 x 2 contingency table for the allele frequency calculated using Equation 16 is created for the number of sitting positions (S) shown in Fig. 17.
  • step SB-7 an independence test is performed for each sitting position on the contingency table created in step SB-6 (step SB-7).
  • step SB—8 Yes
  • the number of samples in the first sample, the second sample Each diplo of the first sample set in step SB-4 based on the number of haplotypes, the total number of haplotypes that can be present, the diplotype form frequency of the first sample, the diploma form frequency of the second sample, and the haplotype frequency.
  • the probability is calculated using the preset distribution of the diplotype frequency of the first and second samples (step SB-9). Specifically, the number of samples of the first sample N, the second
  • each diplotype power of the first sample and each diplotype power of the second sample (X, ⁇ ⁇ , X, X, ⁇ ⁇ ⁇ ⁇ , X)
  • the corresponding probability P is calculated using the joint distribution of the diplotype frequency between the first and second samples shown in Equation 17 below.
  • Equation 17 k is a number that identifies the sample.
  • step SB-10 it is determined whether or not all combinations of the diplotype powers of the first sample and the diplotype powers of the second sample have been completed, and if all the combinations have been completed (step SB-10: Yes), each time step SB-4 to step SB-9 are repeated, the probability calculated in step SB-9 is added, and the added probability is finally set as the probability of type 1 error ( Step SB—11). On the other hand, if all combinations have not been completed (step SB-10: No), execute steps SB-4 to SB-9 and repeat until all combinations are completed.
  • the probability calculation method (1) obtains a first sample and a second sample including a plurality of haplotype information related to haplotypes configured by a plurality of SNPs, and (2 ) Calculate the haplotype frequency for each haplotype that can exist based on the haplotype information contained in the acquired first sample and the haplotype information contained in the acquired second sample, and (3) each haplotype Sitting position minor allele or major allele By defining a predetermined value according to whether or not it has, a matrix with the predetermined value as a component is defined.
  • the diplotype form frequency of each diplotype in the first sample and each diplotype in the second sample By setting the type diplotype power, the combination of each diplotype power of the first sample and each diplotype power of the second sample is given, and (5) corresponds to the set first sample.
  • the allele frequency at each locus in the first sample is calculated based on the multiple diplotype powers and the defined matrix, and the second diplotype powers corresponding to the set second sample and the defined matrix are calculated. Calculate the allele frequencies at each locus in the two samples, and (6) calculate the allele frequencies at each locus in the first sample and the allele frequencies at each locus in the second sample.
  • the predetermined diplotype of the first sample and the second sample of the probability corresponding to the combination of each diplotype form frequency of the first sample and each diplotype form factor of the second sample Calculation is performed using a joint distribution relating to the form frequency.
  • step SB-9 the joint distribution of the first and second samples in terms of the diplotype frequency shown in Equation 17 is the joint distribution of X for all j and j ', that is,
  • the dominant genetic model test is a major allele.
  • the test method that compares the ratio of individuals with those that do not, and the ratio of individuals in two groups.
  • the test with the recessive genetic model is a test that compares the ratio of individuals with minor alleles and those that do not in two groups. Let's go.
  • j j 'may be sufficient. If the total number of haplotypes is M, the total number of permutations diplotype is M 2. If Hardy-Weinberg equilibrium is established for haplotypes, the frequency of D in the population is hh. Now, N people as the first specimen and N people as the second specimen
  • Equation 18 It is a multinomial distribution, and the probability is given by Equation 18 and Equation 19 below for each of the first and second samples.
  • the probability of type 1 error calculated in step SB-11 can be expressed by the following equation 20.
  • the probability calculated by Equation 20 below means the probability that any one of the results of the independence test performed in Step SB-7 will be significant.
  • the probability of the first error calculated by this test coincides with the probability of the first error calculated by the above-described (1 1) allele frequency model.
  • the probability of the first type of error due to the method of performing the three tests of the dominant genetic model, the recessive genetic model, and the allele frequency model and determining that it is significant if any of them is significant should be accurately calculated. Is possible. Specifically, for the contingency tables in Fig. 13, Fig. 16 and Fig. 17, an independence test is performed at each locus, and if the deviation is significant, the function f is set to 1, and Otherwise, it can be 0.
  • FIG. 3 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type of error in the genotype model.
  • a first sample and a second sample containing a plurality of haplotype information related to haplotypes composed of multiple SNPs are acquired (step SC-1).
  • the explanation for obtaining the first specimen and the second specimen is the same as the (1 1) allele frequency model and (1 2) dominant 'recessive genetic model described above.
  • the haplotype frequency is calculated for each haplotype that may exist ( Step SC—2).
  • the explanation for the calculation of the haplotype frequency is the same as the (1 1) allele frequency model and (1 2) dominant / recessive genetic model described above.
  • a predetermined value is set according to the force at which each locus of each haplotype has a minor allele or major allele, thereby defining a matrix having the predetermined value as a component (step SC-3).
  • the description of the matrix definition is the same as the (1-1) array frequency model and (1-2) dominant / recessive genetic model described above.
  • each diplotype power and the first diplotype power in the second sample are set.
  • a combination of two specimens with each diplotype power is given (step SC—4).
  • the explanation regarding the setting of the diplotype form frequency is the same as the above-mentioned (12) dominant / recessive genetic model.
  • the genes at each locus in the first sample based on the multiple diplotype diopter corresponding to the first sample set in step SC-4 and the matrix defined in step SC-3.
  • the genotype at each locus in the second sample is calculated based on the multiple dive mouth type frequencies corresponding to the second sample set in step SC-4 and the matrix defined in step SC-3. Calculate the frequency (step SC—5).
  • the major homozygote (AA) genotype frequency Y at locus i in the first sample and at locus i in the second sample Major homozygous (AA) genotype frequency Y, minor homozygous at locus i in the first sample (a)
  • the genotype frequency Y "of hetero (Aa) at locus i in 2i li and the second sample is
  • step SC-5 based on the genotype frequency at each locus in the first sample and the genotype frequency at each locus in the second sample calculated in step SC-5, the even number corresponding to each locus is calculated.
  • Create a current table step SC—6). Specifically, the 2 ⁇ 3 contingency table shown in FIG. 18 relating to the genotype frequency calculated using Equation 22 is created for the number of loci (L).
  • step SC-7 an independence test is performed for each locus on the contingency table created in step SC-6 (step SC-7).
  • step SC—8 Yes
  • the number of samples in the first sample, the second sample Each diplo of the first sample set in step SC-4 based on the number of samples, the total number of haplotypes that can be present, the diplotype form frequency of the first sample, the diploma form frequency of the second sample, and the haplotype frequency.
  • the probability corresponding to the combination of the type frequency and each diplotype frequency of the second sample is calculated using the preset distribution of the diplotype frequency of the first and second samples (step SC). — 9). Note that the explanation regarding the calculation of the probability is described above. (1 2) Same as dominant / recessive genetic model.
  • step SC-10 it is determined whether or not all combinations of each diplotype power of the first sample and each diplotype power of the second sample have been completed, and if all the combinations have been completed (step SC-10: Yes), each time step SC-4 to step SC-9 are repeated, the probability calculated in step SC-9 is added, and the added probability is finally set as the probability of type 1 error ( Step SC—11). On the other hand, if all combinations have not been completed (step SC-10: No), execute steps SC-4 to SC-9 and repeat until all combinations are completed.
  • the probability calculation method (1) obtains a first sample and a second sample including a plurality of haplotype information related to haplotypes composed of a plurality of SNPs, and (2 ) Calculate the haplotype frequency for each haplotype that can exist based on the haplotype information contained in the acquired first sample and the haplotype information contained in the acquired second sample, and (3) each haplotype By setting a predetermined value according to whether the sitting position has a minor allele or a major allele, a matrix with the predetermined value as a component is defined.
  • each diplotype power for the first sample and each diplotype for the second sample (5) Calculate and set genotypic frequencies at each locus in the first sample based on multiple diplotype shape frequencies corresponding to the set first sample and the defined matrix The genotype frequency at each locus in the second sample is calculated based on a plurality of diplotype frequencies corresponding to the second sample and the defined matrix. (6) Genes at each locus in the calculated first sample Based on the type frequency and the genotype frequency at each locus in the second sample, create an even table corresponding to each locus, and (7) perform an independence test for each locus on the created even table.
  • the number of samples in the first sample, the number of samples in the second sample, the total number of haplotypes that can exist, and the diplotype form of the first sample Frequency, diplotype shape of the second sample Based on the protype frequency, the probability corresponding to the combination of each diplotype form frequency of the first sample and each dive mouth type form frequency of the second sample is set to the preset first and second samples. With specimen Calculation is made using a bond distribution relating to diplotype power.
  • FIG. 4 is a flowchart showing an example of a probability calculation method for accurately calculating the probability of the first type of error when the exact haplotype frequency is inefficient.
  • step SD-1 obtain the first and second samples that contain multiple pieces of haplotype information related to haplotypes composed of multiple SNPs.
  • the explanation for obtaining the first specimen and the second specimen is the same as the above-mentioned (1 1) allele frequency model, (1 2) dominant / recessive genetic model, and (1-3) genotype model.
  • a predetermined value is set according to the force at which each locus of each haplotype has a minor allele or major allele, thereby defining a matrix having the predetermined value as a component (step SD-2).
  • the explanation for the definition of the matrix is the same as the above (1 1) allelic frequency model, (1 2) dominant / recessive genetic model, and (1 3) genotype model.
  • the haplotype frequency X of haplotype number j is a random variable.
  • the actual value X of the haplotype frequency X of the haplotype number j in the first sample is set for all haplotype numbers j.
  • haplotype frequency of haplotype number j in 2 samples Realized value of X Set for type number j.
  • step SD -Four the explanation regarding the calculation of the allele frequency is the same as the above-described (1-1) allele frequency model.
  • step SD-4 based on the allele frequency at each locus in the first sample and the allele frequency at each locus in the second sample calculated in step SD-4, the contingency table corresponding to each locus is calculated. Create (Step SD—5).
  • step SD-6 the independence test is performed for each locus based on the contingency table created in step SD-5 (step SD-6).
  • step SD—7 Yes
  • the number of samples in the first sample, the number of samples in the second sample Based on the total number of haplotypes that can exist, the haplotype frequency of the first sample, and the haplotype frequency of the second sample, each haplotype frequency of the first sample and each haplotype frequency of the second sample set in step SD—3
  • the probability corresponding to the combination is calculated using the preset distribution of the haplotype frequencies of the first and second samples (step SD-8). Specifically, sample number N of the first sample, sample number N of the second sample, haplotypes that may exist
  • each haplo of the first sample set in step SD-3 Based on the actual value X of the type frequency, each haplo of the first sample set in step SD-3
  • the probability P corresponding to 11 1 21 2 is calculated using the joint distribution related to the haplotype frequency between the first and second samples shown in Equation 24 below. Note that the first and second samples shown in Equation 24 below
  • the joint distribution regarding the haplotype frequency is hypergeometric distribution.
  • step SD—9 it is determined whether all combinations of each haplotype frequency of the first sample and each haplotype frequency of the second sample have been completed, and all combinations have been completed.
  • step SD—9 Yes
  • step SD-3 to step SD-8 all the probabilities calculated in step SD-8 are added, and the added probability is finally set as the first type error probability (step SD - Ten).
  • Step SD-9 No
  • Step SD-3 to Step SD-8 execute Step SD-3 to Step SD-8 and repeat until all combinations are completed.
  • the probability calculation method (1) obtains a first sample and a second sample including a plurality of haplotype information related to a haplotype composed of a plurality of SNPs, and (2 ) By defining a predetermined value according to whether each locus of each haplotype has a minor allele or a major allele, a matrix whose component is the predetermined value is defined.
  • the probability corresponding to the combination of each haplotype frequency of the set first sample and each haplotype frequency of the second sample is set in advance. Calculated using the joint distribution of the haplotype frequency between the first and second samples. Then, repeat (3), (4), (5), (6), and (7) for all combinations of each haplotype frequency in the first sample and each haplotype frequency in the second sample. By adding the probabilities calculated in (7), we finally calculate the probability of type 1 error. As a result, it is possible to directly calculate the probability of type 1 error, and as a result, it is possible to prevent oversight of disease-related SNPs.
  • Equation 25 The probability of type 1 error calculated in step SD-10 can be expressed by Equation 25 below.
  • the probability calculated by Equation 25 below means the probability that any one of the results of the independence test performed in step SD-6 will be significant.
  • Equation 2 5 Equation 25
  • Z is 1 if one of the test results of the independence test performed in step SD-6 is significant, and 0 if it is not. It is a random variable defined to take The variable Z is a function of the random variable Y described above, and
  • the random variable Y is a function of the random variable X described above. Also, the random variable X
  • Equation 23 the variable Z can be expressed by the function f (x, ..., X)
  • the probability calculation method according to the present invention for approximately calculating the probability of the first type of error using the MCMC method is as follows: (2-1) Allele frequency model, (2-2) (2) 3) Genotype model, (2) 4)
  • Allele frequency model (2-2) (2) 3) Genotype model, (2) 4)
  • FIG. 5 is a flowchart showing an example of a probability calculation method for approximately calculating the first type error probability using the MCMC method V in the allele frequency model. It is assumed that all haplotype frequencies are known.
  • a predetermined value is set according to whether each locus of each haplotype has a minor allele or a major allele, thereby defining a matrix having the predetermined value as a component (step SE-1).
  • the description regarding the definition of the matrix is the same as the (1-1) allele frequency model in the first embodiment described above.
  • a combination of haplotype frequencies is generated for the first and second samples set in advance (step SE-2). Specifically, the realized value X of the haplotype frequency X of haplotype number j in the first sample and the haplotype j of haplotype number j in the second sample
  • one of the samples of the first sample and the second sample is selected (step SE-3).
  • step SE—6 If the haplotype frequency of the haplotype number j in the selected sample number k is X ⁇ 0 (or X> 0), in step SE-3
  • the haplotype frequency of the first haplotype and the haplotype frequency of the second haplotype are updated based on the haplotype frequency of the first haplotype and the haplotype frequency of the second haplotype set in advance (step SE-7).
  • the value of c defined by Equation 26 below is calculated.
  • h and h are the haplotype frequency of node protype number j and the haplotype frequency of haplotype number j ′, respectively, and their values are preset.
  • kj Keep the values of kj and X.
  • step SE-7 based on a plurality of haplotype frequencies corresponding to the first sample updated in step SE-7 and the matrix defined in step SE-1, the key at each locus in the first sample is determined.
  • the allele at each locus in the second sample is calculated based on the multiple haptic type frequencies corresponding to the second sample after the real frequency is calculated and updated in step SE-7 and the matrix defined in step SE-1 Calculate the frequency (step SE—8).
  • the explanation for the calculation of the allele frequency is the same as the (1-1) allele frequency model of the first embodiment described above.
  • Step SE-9 based on the allele frequencies at each locus in the first sample and the allele frequencies at each locus in the second sample calculated in step SE-8, the contingency table corresponding to each locus is created. (Step SE-9).
  • Step SE-10 an independence test is performed for each locus based on the contingency table created in step SE-9 (step SE-10).
  • Step SE-11 Yes
  • the cumulative number in that case is updated (Step SE). — 12). Specifically, 1 is added to the cumulative number.
  • step SE-13 it is determined whether or not the force has reached the predetermined number of times, and when the predetermined number of times has ended (step SE-13: Yes), the ratio of the final cumulative number to the predetermined number of times (cumulative number ⁇ predetermined number of times) Is calculated as the probability of type 1 error (step SE-14). On the other hand, if the predetermined number of times is over (Step SE-13: No), execute Step SE-3 to Step SE-12 and repeat until the predetermined number of times is completed.
  • the probability calculation method that works for the present invention generates haplotype samples according to the multinomial distribution expected for each sample by the Metropolis-Hastings method, and assigns each sample to each sample. ! /, Perform a test at each position, and find the percentage that is significant at any position.
  • the state of the Markov chain is represented by a contingency table of 2 X M (total number of haplotypes that can exist) based on the actual value of the haplotype frequency.
  • the test is performed by sampling a conservative test method or other haplotype copy number.
  • the probability of type 1 error can be calculated by calculating the proportion of steps determined to be significant.
  • the first and second samples are independently generated by the Monte Carlo sampler that follows the multinomial distribution, and the above-described calculation of the function f is performed. A method is conceivable.
  • Figure 6 shows a dominant / recessive genetic model that approximates the probability of type 1 error using the MCMC method. It is a flowchart which shows an example of the probability calculation method calculated in step. It is assumed that all haplotype frequencies are known.
  • a predetermined value is set according to whether each locus of each haplotype has a minor allele or a major allele, thereby defining a matrix having the predetermined value as a component (step SF-1).
  • the description of the matrix definition is the same as the (1-2) dominant / recessive genetic model in the first embodiment described above.
  • a combination of diplotype powers is generated for the first and second samples set in advance (step SF-2). Specifically, the actual value X of the diplotype form frequency of diplotype form number ⁇ 'in the first sample is set to ⁇ ⁇ ' for all diplotype form numbers ⁇ '.
  • the realization value X of the diplotype form frequency of the diplotype form number jj 'in the second sample is generated so as to satisfy the above-mentioned equation 14. (X, ..., X, X
  • step SF—3 either the first sample or the second sample is selected.
  • Step SF—4 select the first diplotype. Specifically, the diplotype model number '(1 ⁇ ordered natural numbers i and j ⁇ M ⁇ S 1 ⁇ ); L is the number of sitting positions) is selected.
  • the diplotype frequency of the first diplotype form in the specimen selected in step SF-3 is 0 (SF-5: Yes)
  • the first diplotype form is different from the first diplotype form. 2
  • Select the diplotype type (Step SF—6). Specifically, the diplotype shape number j j 'of the selected sample number k sample x half 0 (or X
  • the value of the diplotype power X is maintained as kjljl 'kjljl'.
  • the diplotype of the first diplotype Degree and second diplotype Update the diplotype form factor for (step SF—7). Specifically, first, the value of c defined by Equation 27 below is calculated. In Equation 27 below, h, h, h and h
  • jl jl '] 2 is the haplotype frequency of haplotype number j and the haplotype of j'
  • each locus in the first sample is Based on multiple diplotype powers corresponding to the second sample after calculation of the allele frequency and updated in step SF—7 and the matrix defined in step SF— !, and each locus in the second sample Calculate the allele frequency at (step SF—8).
  • the explanation for calculating the allele frequency is the same as the (12) dominant 'recessive genetic model of the first embodiment described above.
  • Step SF-9 based on the allele frequencies at each locus in the first sample and the allele frequencies at each locus in the second sample calculated in step SF-8, the contingency table corresponding to each locus is created. (Step SF-9).
  • step SF-10 the independence test is performed for each locus based on the contingency table created in step SF-9 (step SF-10).
  • step SF—11 Yes
  • step SF— 12 update the cumulative number in that case. Specifically, 1 is added to the cumulative number.
  • Step SF—13 Yes
  • Step SF-13 No
  • Step SF-3 Step SF-12 and repeat until the predetermined number of times ends.
  • the probability calculation method according to the present invention generates a diplotype sample according to the multinomial distribution expected for each sample by the Metropolis-Hastings method, and each locus for each sample. Perform a test at, and find the proportion that is significant at any of these loci.
  • FIG. 7 is a flowchart showing an example of a probability calculation method for approximately calculating the first type error probability using the MCMC method in the genotype model. It is assumed that all haplotype frequencies are known.
  • a matrix having the predetermined value as a component is defined by setting a predetermined value according to whether each locus of each haplotype has a minor allele or a major allele (step SG-1).
  • the description regarding the definition of the matrix is the same as the (1-3) genotype model in the first embodiment described above.
  • step SG-2 a combination of diplotype form frequencies is generated for the first and second samples set in advance (step SG-2).
  • the explanation for the generation of the diplotype form frequency is the same as in (2-2) Dominant / recessive genetic model described above.
  • step SG-3 either the first sample or the second sample is selected.
  • step SG-3 if the diplotype frequency of the first diplotype form in the sample selected in step SG-3 is 0 (SG-5: Yes), the second diplotype form is different from the first diplotype form.
  • step SG—6 The explanation for the selection of the second diplotype is the same as in (2-2) Dominant / recessive genetic model described above.
  • the diplotype of the first diplotype is updated (step SG7 ).
  • the explanation for updating the diplotype form frequency is the same as the (2-2) dominant / recessive genetic model described above.
  • genotypes at each locus in the first sample based on the multiple diplotypes corresponding to the first sample after updating in step SG-7 and the matrix defined in step SG-1
  • the genotype at each locus in the second sample based on the multiple diplotype frequencies corresponding to the second sample after updating the frequency in step SG-7 and the matrix defined in step SG-1 Calculate the frequency (step SG—8).
  • the description regarding the calculation of the genotype frequency is the same as the (1-3) genotype model of the first embodiment described above.
  • step SG-8 based on the genotype frequency at each locus in the first sample and the genotype frequency at each locus in the second sample calculated in step SG-8, the contingency table corresponding to each locus. (Step SG—9).
  • step SG-10 the independence test is performed for each sitting position using the contingency table created in step SG-9 (step SG-10).
  • step SG—11 Yes
  • step SG— 12 the cumulative number in that case is updated. Specifically, 1 is added to the cumulative number.
  • step SG—13: Yes the ratio of the final cumulative number to the predetermined number (cumulative number ⁇ predetermined number of times) Is calculated as the probability of type 1 error.
  • step SG-13: No the ratio of the final cumulative number to the predetermined number (cumulative number ⁇ predetermined number of times) Is calculated as the probability of type 1 error.
  • the probability calculation method that works on the present invention is Metropolis- Hastings.
  • a diplotype sample that follows the expected multinomial distribution for each specimen is generated by the method, and a test is performed at each locus for each sample, and a ratio that is significant at any one of the loci is obtained.
  • the true haplotype frequency is often unknown. Therefore, even if the haplotype frequency is not weak, if the diplotype shape in the specimen is weak, it can be calculated as follows.
  • Permutation haplotype copy numbers follow a multinomial distribution given the haplotype frequency.
  • probability under specific observational data rather than the probability under certain conditions of the sample space. For example, consider the probability under the condition of a set of results that satisfy the restriction as shown in Equation 28 or 29 below.
  • the above-described dominant genetic model, recessive genetic model, allele frequency model, and genotype model can all be tested if the combined diplotype copy number is known.
  • ⁇ j ⁇ i follows hypergeometric distribution.
  • the probabilities are expressed by the following formula 30 and the following formula 31, respectively.
  • Equation 31 (A-2) Generate a sample according to the hypergeometric distribution shown in Equation 31 using the MCMC method. However, ⁇ y ⁇ satisfies Equation 32 below.
  • Equation 32 (A-2)
  • A-3 The generated sample is tested at each locus using a dominant genetic model, a recessive genetic model, an allele frequency model, and a genotype model.
  • the number of sitting positions is given under certain constraints.
  • the state space of the Markov chain (Markov-chain) is a plurality of element y forces that take different values under the predetermined constraint.
  • N 1 A given (given) fixed non-negative integer value.
  • N is included in sample number k
  • the diplotype model number uv is determined by selecting two ordered integers (u, V) with equal probability.
  • U is ;
  • L satisfies the number of sitting positions, and
  • V satisfies "l ⁇ v ⁇ u”.
  • the diplotype shape number uv of the selected sample number k sample uv has two orders different from the selected (u, V) when y
  • the diplotype model number WS is further determined by further selecting (W, S), which is an integer.
  • the power described for the independence test using one of the three models (dominant genetic model, recessive genetic model, and genotype model).
  • the probability calculation method according to the present invention includes three methods. It can be easily extended to an independence test using the entire model. Therefore, as described in (B-6) above, the independence test is performed at each sitting position using three different contingency tables. Independence test If the constant test result is significant, the test result of the overall (all loci) independence test is significant.
  • the probability calculation method according to the present invention can directly calculate the probability of the first type of error in consideration of the haplotype frequency, and as a result, the disease-related SNP can be calculated. Oversight can be prevented. Therefore, the probability calculation method according to the present invention is extremely useful in fields such as medicine and drug discovery.

Landscapes

  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Genetics & Genomics (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Chemical & Material Sciences (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Analytical Chemistry (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

　第１種の過誤の確率を、ハプロタイプ頻度を考慮して直接的に算出することができ、その結果、疾患関連ＳＮＰの見落としを防ぐことができる確率算出方法を提供することを課題とする。本発明にかかる確率算出方法は、（１）標本を取得し、（２）ハプロタイプ頻度を算出し、（３）各アレルに対応する所定値を成分とする行列を定義し、（４）ハプロタイプ度数の組合せを与え、（５）各座位でのアレル度数を算出し、（６）各座位に対応する偶現表を作成し、（７）偶現表につき、座位ごとに独立性検定を実行し、（８）独立性検定の検定結果のいずれか１つが有意であった場合、ハプロタイプ度数の組合せに対応する確率を結合分布を用いて算出する。そして、（４）、（５）、（６）、（７）および（８）を、ハプロタイプ度数の全ての組合せについて繰り返し、繰り返すごとに（８）で算出される確率を加算することで、第１種の過誤の確率を算出する。

Description

明 細 書
確率算出方法
技術分野
[0001] 本発明は、 SNP (Single Nucleotide Polymorphism)を用いた case— contro 1相関解析で行う独立性検定における第 1種の過誤の確率 (有意水準)を算出する確 率算出方法に関するものである。
背景技術
[0002] case— control研究では、ある疾患に罹った者またはある病態にある者（case)とそ うでない者（control)とがあらかじめどのような危険因子にさらされたかを調べる。 SN Pを用いた case— control相関解析では、当該危険因子として、（1)アレル頻度モデ ル (AZa)、 (2)優性遺伝モデル (AA+Aa/aa)、（3)劣性遺伝モデル (AA/Aa + aa)、 （4)遺伝子型モデル (AAZAaZaa)、を想定する。なお、 Aはメジャーァレ ル（major allele)、 aはマイナーアレル（minor allele)を示す。そして、各モデルに っ 、て、 2 X 2または 2 X 3の偶現表を作成し、独立性検定を行う。
[0003] し力し、上述したモデル間で共通したデータを使っているので、多重検定の偶然の 影響を除いた上で依然として検定の統計的有意差が認められるかが課題となる。
[0004] そこで、多重検定の影響を除外するための方法として、まず、 Bonferroniの補正の 方法が知られている。 Bonferroniの補正の方法では、独立性検定で得られた p値に 独立性検定の回数 (上述した 4つのモデルの場合では 4回）を乗じた値を多重検定 全体力 得られた P値とみなすことで、多重検定の影響を除外している。この方法は、 検定結果を最も保守的 (厳格)に解釈するものである。また、多重検定の影響を除外 するための方法として、つぎに、非特許文献 1に記載の方法が知られている。非特許 文献 1に記載の方法では、少数個の SNPを対象として各 2つの SNP間の連鎖不平 衡指標値（ Δ )を計算し当該指標値を成分とする行列の固有値分解を行うことで近似 的な補正を施している。つまり、非特許文献 1に記載の方法は連鎖不平衡に着目し た多重検定における p値の補正の方法である。
[0005] 特許文献 1 : Dale R. Nyholt, A simple correction for multiple test ing for single― nucleotide polymorphisms in linkage disequilibrium with each other, Am. J. Hum. Genet. , Vol. 74, pp. 765— 769 発明の開示
発明が解決しょうとする課題
[0006] しかしながら、 Bonferroniの補正の方法では、補正された p値が真の p値より大きい 可能性が高いので、疾患関連 SNPを見落とす危険性がある、という問題点があった 。また、大量の SNPデータを対象とした場合、多重検定全体の p値を一般的に用いら れている例えば 0. 01や 0. 05などに設定すると、設定した p値を許容するためには 各 SNPの p値が非常に小さな値、つまり非常に厳しい条件となるので、同様に疾患 関連 SNPを見落とす危険性がある、という問題点があった。例えば、 50000個の SN Pデータを対象として、多重検定全体の p値を 0. 01に設定した場合、各 SNPの p値 は 0. 0000002 (0. 01 + 50000)と 常に/ J、さな値となる。
[0007] また、非特許文献 1に記載の方法では、ペアワイズ (pairwise)な連鎖不平衡に着 目しているので、ハプロタイプ頻度が考慮されていない、という問題点があった。さら に、直接的な p値を算出していない、という問題点があった。
[0008] 本発明は上記問題点に鑑みてなされたもので、 SNPを用いた case— control相関 解析で行う独立性検定における P値、つまり第 1種の過誤の確率を、ハプロタイプ頻 度を考慮して直接的に算出することができ、その結果、疾患関連 SNPの見落としを 防ぐことができる確率算出方法を提供することを目的とする。
課題を解決するための手段
[0009] 上記目的を達成するために、本発明にかかる請求項 1に記載の確率算出方法は、 SNPを用いた case— control相関解析で行う独立性検定における第 1種の過誤の 確率を算出する確率算出方法において、複数の SNPで構成されるハプロタイプに関 するハプロタイプ情報を複数含む第 1標本および第 2標本を取得する標本取得ステ ップと、前記標本取得ステップで取得した第 1標本に含まれるハプロタイプ情報およ び取得した第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイ プごとにハプロタイプ頻度を算出するハプロタイプ頻度算出ステップと、各ハプロタイ プの各座位がマイナーアレルまたはメジャーアレルを有するかに従って予め定めた 所定値を設定することで、当該所定値を成分とする行列を定義する行列定義ステツ プと、第 1標本における各ハプロタイプのハプロタイプ度数および第 2標本における 各ハプロタイプのハプロタイプ度数を設定することで、第 1標本の各ハプロタイプ度数 と第 2標本の各ハプロタイプ度数との組合せを与えるハプロタイプ度数設定ステップ と、前記ハプロタイプ度数設定ステップで設定した第 1標本に対応する複数のハプロ タイプ度数および前記行列定義ステップで定義した行列に基づいて第 1標本におけ る各座位でのアレル度数を算出し、設定した第 2標本に対応する複数のハプロタイプ 度数および定義した行列に基づいて第 2標本における各座位でのアレル度数を算 出するアレル度数算出ステップと、前記アレル度数算出ステップで算出した第 1標本 における各座位でのアレル度数および第 2標本における各座位でのアレル度数に基 づいて、各座位に対応する偶現表を作成する偶現表作成ステップと、前記偶現表作 成ステップで作成した偶現表にっ ヽて、座位ごとに独立性検定を実行する独立性検 定実行ステップと、前記独立性検定実行ステップで実行した独立性検定の検定結果 のいずれか 1つが有意であった場合、第 1標本の標本数、第 2標本の標本数、存在し 得るハプロタイプの総数、第 1標本のハプロタイプ度数、第 2標本のハプロタイプ度数 およびノヽプロタイプ頻度に基づ ヽて、前記ハプロタイプ度数設定ステップで設定した 第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイプ度数との組合せに対応す る確率を、予め設定した第 1標本と第 2標本とのハプロタイプ度数に関する結合分布 を用いて算出する確率算出ステップと、を含み、前記ハプロタイプ度数設定ステップ 、前記アレル度数算出ステップ、前記偶現表作成ステップ、前記独立性検定実行ス テツプおよび前記確率算出ステップを、第 1標本の各ハプロタイプ度数と第 2標本の 各ハプロタイプ度数との全ての組合せについて繰り返し、繰り返すごとに算出される 前記確率を加算することで、最終的に第 1種の過誤の確率を算出すること、を特徴と する。
また、本発明に力かる請求項 2に記載の確率算出方法は、 SNPを用いた case— co ntro湘関解析で行う独立性検定における第 1種の過誤の確率を算出する確率算出 方法にお 、て、複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を 複数含む第 1標本および第 2標本を取得する標本取得ステップと、前記標本取得ス テツプで取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に含 まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロタイプ頻度 を算出するハプロタイプ頻度算出ステップと、各ハプロタイプの各座位がマイナーァ レルまたはメジャーアレルを有するかに従って予め定めた所定値を設定することで、 当該所定値を成分とする行列を定義する行列定義ステップと、第 1標本における各 ディプロタイプ形のディプロタイプ形度数および第 2標本における各ディプロタイプ形 のディプロタイプ形度数を設定することで、第 1標本の各ディプロタイプ形度数と第 2 標本の各ディプロタイプ形度数との組合せを与えるディプロタイプ形度数設定ステツ プと、前記ディプロタイプ形度数設定ステップで設定した第 1標本に対応する複数の ディプロタイプ形度数および前記行列定義ステップで定義した行列に基づいて第 1 標本における各座位でのアレル度数を算出し、設定した第 2標本に対応する複数の ディプロタイプ形度数および定義した行列に基づ ヽて第 2標本における各座位での アレル度数を算出するアレル度数算出ステップと、前記アレル度数算出ステップで算 出した第 1標本における各座位でのアレル度数および第 2標本における各座位での アレル度数に基づ ヽて、各座位に対応する偶現表を作成する偶現表作成ステップと 、前記偶現表作成ステップで作成した偶現表について、座位ごとに独立性検定を実 行する独立性検定実行ステップと、前記独立性検定実行ステップで実行した独立性 検定の検定結果のいずれか 1つが有意であった場合、第 1標本の標本数、第 2標本 の標本数、存在し得るハプロタイプの総数、第 1標本のディプロタイプ形度数、第 2標 本のディプロタイプ形度数およびノヽプロタイプ頻度に基づ ヽて、前記ディプロタイプ 形度数設定ステップで設定した第 1標本の各ディプロタイプ形度数と第 2標本の各デ ィプロタイプ形度数との組合せに対応する確率を、予め設定した第 1標本と第 2標本 とのディプロタイプ形度数に関する結合分布を用いて算出する確率算出ステップと、 を含み、前記ディプロタイプ形度数設定ステップ、前記アレル度数算出ステップ、前 記偶現表作成ステップ、前記独立性検定実行ステップおよび前記確率算出ステップ を、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との全て の組合せについて繰り返し、繰り返すごとに算出される前記確率を加算することで、 最終的に第 1種の過誤の確率を算出すること、を特徴とする。 また、本発明に力かる請求項 3に記載の確率算出方法は、 SNPを用いた case— co ntro湘関解析で行う独立性検定における第 1種の過誤の確率を算出する確率算出 方法にお 、て、複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を 複数含む第 1標本および第 2標本を取得する標本取得ステップと、前記標本取得ス テツプで取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に含 まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロタイプ頻度 を算出するハプロタイプ頻度算出ステップと、各ハプロタイプの各座位がマイナーァ レルまたはメジャーアレルを有するかに従って予め定めた所定値を設定することで、 当該所定値を成分とする行列を定義する行列定義ステップと、第 1標本における各 ディプロタイプ形のディプロタイプ形度数および第 2標本における各ディプロタイプ形 のディプロタイプ形度数を設定することで、第 1標本の各ディプロタイプ形度数と第 2 標本の各ディプロタイプ形度数との組合せを与えるディプロタイプ形度数設定ステツ プと、前記ディプロタイプ形度数設定ステップで設定した第 1標本に対応する複数の ディプロタイプ形度数および前記行列定義ステップで定義した行列に基づいて第 1 標本における各座位での遺伝子型度数を算出し、設定した第 2標本に対応する複数 のディプロタイプ形度数および定義した行列に基づ 、て第 2標本における各座位で の遺伝子型度数を算出する遺伝子型度数算出ステップと、前記遺伝子型度数算出 ステップで算出した第 1標本における各座位での遺伝子型度数および第 2標本にお ける各座位での遺伝子型度数に基づ ヽて、各座位に対応する偶現表を作成する偶 現表作成ステップと、前記偶現表作成ステップで作成した偶現表について、座位ごと に独立性検定を実行する独立性検定実行ステップと、前記独立性検定実行ステップ で実行した独立性検定の検定結果の!/、ずれか 1つが有意であった場合、第 1標本の 標本数、第 2標本の標本数、存在し得るハプロタイプの総数、第 1標本のディプロタイ プ形度数、第 2標本のディプロタイプ形度数およびノヽプロタイプ頻度に基づいて、前 記ディプロタイプ形度数設定ステップで設定した第 1標本の各ディプロタイプ形度数 と第 2標本の各ディプロタイプ形度数との組合せに対応する確率を、予め設定した第 1標本と第 2標本とのディプロタイプ形度数に関する結合分布を用いて算出する確率 算出ステップと、を含み、前記ディプロタイプ形度数設定ステップ、前記遺伝子型度 数算出ステップ、前記偶現表作成ステップ、前記独立性検定実行ステップおよび前 記確率算出ステップを、第 1標本の各ディプロタイプ形度数と第 2標本の各デイブロタ イブ形度数との全ての組合せについて繰り返し、繰り返すごとに算出される前記確率 を加算することで、最終的に第 1種の過誤の確率を算出すること、を特徴とする。 また、本発明に力かる請求項 4に記載の確率算出方法は、 SNPを用いた case— co ntrol相関解析で行う独立性検定における第 1種の過誤の確率を MCMC法を用い て近似的に算出する確率算出方法であって、各ハプロタイプの各座位がマイナーァ レルまたはメジャーアレルを有するかに従って予め定めた所定値を設定することで、 当該所定値を成分とする行列を定義する行列定義ステップと、ハプロタイプ度数の組 合せを予め設定した第 1標本および第 2標本に対して生成するハプロタイプ度数生 成ステップと、第 1標本または第 2標本の 、ずれかの標本を選択する標本選択ステツ プと、第 1ハプロタイプを選択する第 1ハプロタイプ選択ステップと、前記標本選択ス テツプで選択した標本における前記第 1ハプロタイプのハプロタイプ度数力^でない 場合、当該第 1ハプロタイプとは異なる第 2ハプロタイプを選択する第 2ハプロタイプ 選択ステップと、予め設定した前記第 1ハプロタイプのハプロタイプ頻度および前記 第 2ハプロタイプのハプロタイプ頻度に基づいて、前記第 1ハプロタイプのハプロタイ プ度数および前記第 2ハプロタイプのハプロタイプ度数を更新するハプロタイプ度数 更新ステップと、前記ハプロタイプ度数更新ステップで更新した後の第 1標本に対応 する複数のハプロタイプ度数および前記行列定義ステップで定義した行列に基づ ヽ て第 1標本における各座位でのアレル度数を算出し、更新した後の第 2標本に対応 する複数のハプロタイプ度数および定義した行列に基づいて第 2標本における各座 位でのアレル度数を算出するアレル度数算出ステップと、前記アレル度数算出ステツ プで算出した第 1標本における各座位でのアレル度数および第 2標本における各座 位でのアレル度数に基づ ヽて、各座位に対応する偶現表を作成する偶現表作成ス テツプと、前記偶現表作成ステップで作成した偶現表について、座位ごとに独立性検 定を実行する独立性検定実行ステップと、前記独立性検定実行ステップで実行した 独立性検定の検定結果のいずれか 1つが有意であった場合、当該場合の累積回数 を更新する累積回数更新ステップと、を含み、前記標本選択ステップ、前記第 1ハプ 口タイプ選択ステップ、前記第 2ハプロタイプ選択ステップ、前記ハプロタイプ度数更 新ステップ、前記アレル度数算出ステップ、前記偶現表作成ステップ、前記独立性検 定実行ステップおよび前記累積回数更新ステップを所定回数繰り返し、前記所定回 数に対する最終的な累積回数の割合を第 1種の過誤の確率として算出すること、を 特徴とする。
発明の効果
本発明にかかる請求項 1に記載の確率算出方法は、（1)複数の SNPで構成される ハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を取得し 、 (2)取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に含ま れるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロタイプ頻度を 算出し、（3)各ハプロタイプの各座位がマイナーアレル（同じ座位で頻度の低い方の アレル）またはメジャーアレル（同じ座位で頻度の高い方のアレル）を有するかに従つ て予め定めた所定値を設定することで、当該所定値を成分とする行列を定義し、（4) 第 1標本における各ハプロタイプのハプロタイプ度数および第 2標本における各ハプ 口タイプのハプロタイプ度数を設定することで、第 1標本の各ハプロタイプ度数と第 2 標本の各ハプロタイプ度数との組合せを与え、（5)設定した第 1標本に対応する複数 のハプロタイプ度数および定義した行列に基づいて第 1標本における各座位でのァ レル度数を算出し、設定した第 2標本に対応する複数のハプロタイプ度数および定 義した行列に基づいて第 2標本における各座位でのアレル度数を算出し、（6)算出 した第 1標本における各座位でのアレル度数および第 2標本における各座位でのァ レル度数に基づいて、各座位に対応する偶現表（なお、偶現表とは contingency t ableのことであり、 日本語では分割表と訳すこともある。）を作成し、（7)作成した偶現 表について、座位ごとに独立性検定を実行し、（8)実行した独立性検定の検定結果 のいずれか 1つが有意であった場合、第 1標本の標本数、第 2標本の標本数、存在し 得るハプロタイプの総数、第 1標本のハプロタイプ度数、第 2標本のハプロタイプ度数 およびノヽプロタイプ頻度に基づいて、設定した第 1標本の各ハプロタイプ度数と第 2 標本の各ハプロタイプ度数との組合せに対応する確率を、予め設定した第 1標本と 第 2標本とのハプロタイプ度数に関する結合分布を用いて算出する。そして、（4)、 ( 5)、（6)、（7)および (8)を、第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイ プ度数との全ての組合せについて繰り返し、繰り返すごとに（8)で算出される確率を 加算することで、最終的に第 1種の過誤の確率を算出する。これにより、第 1種の過 誤の確率を、ハプロタイプ頻度を考慮して直接的に算出することができ、その結果、 疾患関連 SNPの見落としを防ぐことができる、という効果を奏する。
本発明にかかる請求項 2に記載の確率算出方法は、（1)複数の SNPで構成される ハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を取得し 、 (2)取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に含ま れるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロタイプ頻度を 算出し、 (3)各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有す るかに従って予め定めた所定値を設定することで、当該所定値を成分とする行列を 定義し、（4)第 1標本における各ディプロタイプ形のディプロタイプ形度数および第 2 標本における各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1標本 の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与え、 ( 5)設定した第 1標本に対応する複数のディプロタイプ形度数および定義した行列に 基づいて第 1標本における各座位でのアレル度数を算出し、設定した第 2標本に対 応する複数のディプロタイプ形度数および定義した行列に基づいて第 2標本におけ る各座位でのアレル度数を算出し、（6)算出した第 1標本における各座位でのアレル 度数および第 2標本における各座位でのアレル度数に基づ 、て、各座位に対応する 偶現表を作成し、（7)作成した偶現表について、座位ごとに独立性検定を実行し、 ( 8)実行した独立性検定の検定結果のいずれか 1つが有意であった場合、第 1標本 の標本数、第 2標本の標本数、存在し得るハプロタイプの総数、第 1標本のディプロ タイプ形度数、第 2標本のディプロタイプ形度数およびノヽプロタイプ頻度に基づ 、て 、設定した第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数と の組合せに対応する確率を、予め設定した第 1標本と第 2標本とのディプロタイプ形 度数に関する結合分布を用いて算出する。そして、（4)、（5)、（6)、（7)および (8)を 、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との全ての 組合せについて繰り返し、繰り返すごとに（8)で算出される確率を加算することで、最 終的に第 1種の過誤の確率を算出する。これにより、第 1種の過誤の確率を、ハプロ タイプ頻度を考慮して直接的に算出することができ、その結果、疾患関連 SNPの見 落としを防ぐことができる、という効果を奏する。
また、本発明にかかる請求項 3に記載の確率算出方法は、（1)複数の SNPで構成 されるハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を 取得し、 (2)取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本 に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロタイプ 頻度を算出し、 (3)各ハプロタイプの各座位がマイナーアレルまたはメジャーアレル を有するかに従って予め定めた所定値を設定することで、当該所定値を成分とする 行列を定義し、（4)第 1標本における各ディプロタイプ形のディプロタイプ形度数およ び第 2標本における各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを 与え、（5)設定した第 1標本に対応する複数のディプロタイプ形度数および定義した 行列に基づいて第 1標本における各座位での遺伝子型度数を算出し、設定した第 2 標本に対応する複数のディプロタイプ形度数および定義した行列に基づいて第 2標 本における各座位での遺伝子型度数を算出し、（6)算出した第 1標本における各座 位での遺伝子型度数および第 2標本における各座位での遺伝子型度数に基づいて 、各座位に対応する偶現表を作成し、（7)作成した偶現表について、座位ごとに独 立性検定を実行し、 (8)実行した独立性検定の検定結果の 、ずれか 1つが有意であ つた場合、第 1標本の標本数、第 2標本の標本数、存在し得るハプロタイプの総数、 第 1標本のディプロタイプ形度数、第 2標本のディプロタイプ形度数およびノヽプロタイ プ頻度に基づ!、て、設定した第 1標本の各ディプロタイプ形度数と第 2標本の各ディ プロタイプ形度数との組合せに対応する確率を、予め設定した第 1標本と第 2標本と のディプロタイプ形度数に関する結合分布を用いて算出する。そして、（4)、 （5)、 (6 )、 （7)および (8)を、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイ プ形度数との全ての組合せについて繰り返し、繰り返すごとに（8)で算出される確率 を加算することで、最終的に第 1種の過誤の確率を算出する。これにより、第 1種の過 誤の確率を、ハプロタイプ頻度を考慮して直接的に算出することができ、その結果、 疾患関連 SNPの見落としを防ぐことができる、という効果を奏する。
[0016] また、本発明にかかる請求項 4に記載の確率算出方法は、（1)各ハプロタイプの各 座位がマイナーアレルまたはメジャーアレルを有するかに従って予め定めた所定値 を設定することで、当該所定値を成分とする行列を定義し、（2)ハプロタイプ度数の 組合せを予め設定した第 1標本および第 2標本に対して生成し、（3)第 1標本または 第 2標本のいずれかの標本を選択し、（4)第 1ハプロタイプを選択し、（5)選択した標 本における第 1ハプロタイプのハプロタイプ度数が 0でない場合、当該第 1ハプロタイ プとは異なる第 2ハプロタイプを選択し、 (6)予め設定した第 1ハプロタイプのハプロ タイプ頻度および第 2ハプロタイプのハプロタイプ頻度に基づ 、て、第 1ハプロタイプ のハプロタイプ度数および第 2ハプロタイプのハプロタイプ度数を更新し、 (7)更新し た後の第 1標本に対応する複数のハプロタイプ度数および定義した行列に基づいて 第 1標本における各座位でのアレル度数を算出し、更新した後の第 2標本に対応す る複数のハプロタイプ度数および定義した行列に基づいて第 2標本における各座位 でのアレル度数を算出し、（8)算出した第 1標本における各座位でのアレル度数およ び第 2標本における各座位でのアレル度数に基づ 、て、各座位に対応する偶現表を 作成し、（9)作成した偶現表について、座位ごとに独立性検定を実行し、（10)実行 した独立性検定の検定結果のいずれか 1つが有意であった場合、当該場合の累積 回数を更新する。そして、（3)、（4)、 (5)、 (6)、 (7)、（8)、 (9)および（10)を所定回 数繰り返し、所定回数に対する最終的な累積回数の割合を第 1種の過誤の確率とし て算出する。これにより、第 1種の過誤の確率を、ハプロタイプ頻度を考慮して直接的 に算出することができ、その結果、疾患関連 SNPの見落としを防ぐことができる、とい う効果を奏する。また、本確率算出方法では第 1種の過誤の確率を MCMC法で近 似的に算出しているので、計算時間を短縮することができ、実用性がより高まる、とい う効果を奏する。
図面の簡単な説明
[0017] [図 1]図 1は、アレル頻度モデルで第 1種の過誤の確率を正確に算出する確率算出 方法の一例を示すフローチャートである。
[図 2]図 2は、優性 ·劣性遺伝モデルで第 1種の過誤の確率を正確に算出する確率算 出方法の一例を示すフローチャートである。
[図 3]図 3は、遺伝子型モデルで第 1種の過誤の確率を正確に算出する確率算出方 法の一例を示すフローチャートである。
[図 4]図 4は、正確なハプロタイプ頻度がわ力 ない場合に第 1種の過誤の確率を正 確に算出する確率算出方法の一例を示すフローチャートである。
[図 5]図 5は、アレル頻度モデルで第 1種の過誤の確率を MCMC法を用いて近似的 に算出する確率算出方法の一例を示すフローチャートである。
[図 6]図 6は、優性.劣性遺伝モデルで第 1種の過誤の確率を MCMC法を用いて近 似的に算出する確率算出方法の一例を示すフローチャートである。
[図 7]図 7は、遺伝子型モデルで第 1種の過誤の確率を MCMC法を用いて近似的に 算出する確率算出方法の一例を示すフローチャートである。
圆 8]図 8は、第 1標本および第 2標本に含まれる情報の一例を示す図である。
[図 9]図 9は、存在し得る全てのハプロタイプに関する情報の一例を示す図である。
[図 10]図 10は、存在し得る全てのハプロタイプに関するハプロタイプ頻度の一例を示 す図である。
[図 11]図 11は、定義した行列の一例を示す図である。
[図 12]図 12は、ハプロタイプ度数に関する情報の一例を示す図である。
[図 13]図 13は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。
[図 14]図 14は、ディプロタイプ形番号とハプロタイプ番号との関係を示す図である。
[図 15]図 15は、ディプロタイプ形度数に関する情報の一例を示す図である。
[図 16]図 16は、劣性遺伝モデルの場合に作成する偶現表の一例を示す図である。
[図 17]図 17は、優性遺伝モデルの場合に作成する偶現表の一例を示す図である。
[図 18]図 18は、遺伝子型モデルの場合に作成する偶現表の一例を示す図である。
[図 19]図 19は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。
[図 20]図 20は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。
[図 21]図 21は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。
[図 22]図 22は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。
[図 23]図 23は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。 [図 24]図 24は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。
[図 25]図 25は、アレル頻度モデルの場合に作成する偶現表の一例を示す図である。 発明を実施するための最良の形態
[0018] 以下に、本発明にかかる確率算出方法の実施の形態を図面に基づいて詳細に説 明する。なお、この実施の形態によりこの発明が限定されるものではない。
[0019] [第 1実施形態]
ここでは、第 1種の過誤の確率を正確に算出する本発明に力かる確率算出方法に ついて、（1 1)アレル頻度モデル、 （1 2)優性.劣性遺伝モデル、 （1 3)遺伝子 型モデル、 （1—4)正確なハプロタイプ頻度がわ力もない場合、の順に図を参照して 詳細に説明する。
[0020] (1 - 1)アレル頻度モデル
図 1は、アレル頻度モデルで第 1種の過誤の確率を正確に算出する確率算出方法 の一例を示すフローチャートである。
まず、複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む 第 1標本および第 2標本を取得する (ステップ SA— 1)。具体的には、図 8に示すよう に、 L (Lは正の整数)個の連鎖した SNPで構成されるハプロタイプ情報を N個含む
1 第 1標本および L個の連鎖した SNPで構成されるハプロタイプ情報を N個含む第 2
2
標本を取得する。
[0021] っ 、で、ステップ SA— 1で取得した第 1標本に含まれるハプロタイプ情報および第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロ タイプ頻度を算出する（ステップ SA— 2)。ここで、 SNPの場合には 2つのアレルが存 在するので、存在し得るハプロタイプの総数は、具体的には、図 9に示すように 2 以 下、 Mと記載することがある）個である。そして、図 10に示すように、存在し得る M個 のハプロタイプのそれぞれに対して、ステップ SA— 1で取得した第 1標本に含まれる
Figure imgf000014_0001
、て、ハプロタ イブ頻度 h (jはハプロタイプを識別する番号であり、 l≤j≤Mを満たす整数)を算出
J
する。なお、ハプロタイプ頻度 hは下記数式 1を満たす。
[数 1] M
· · · (数式 1 )
7=1
[0022] ついで、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有する 力に従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定 義する（ステップ SA—3)。具体的には、ハプロタイプ番号 jのハプロタイプにおける i 番目（iは座位を識別する番号であり、 l≤i≤Lを満たす整数)の座位がマイナーァレ ル (a)である場合は 1を、メジャーアレル (A)である場合は 0を所定値として定める。そ して、定義する行列は、図 11のように表すことができ、座位番号 i (iは l≤i≤Lを満た す整数）とハプロタイプ番号 jとの組合せごとに、 0または 1の値を格納する。
[0023] ついで、第 1標本における各ハプロタイプのハプロタイプ度数および第 2標本にお ける各ハプロタイプのハプロタイプ度数を設定することで、第 1標本の各ハプロタイプ 度数と第 2標本の各ハプロタイプ度数との組合せを与える (ステップ SA—4)。ここで 、第 1標本におけるハプロタイプ番号 jのハプロタイプ度数 Xおよび第 2標本における
1]
ハプロタイプ番号 jのハプロタイプ度数 Xは、ランダム変数であり、それぞれ下記数式
2j
2を満たす。そして、下記数式 2を満たすように、第 1標本におけるハプロタイプ番号 j のハプロタイプ度数 Xの実現値 Xを全てのハプロタイプ番号 jに対して設定し、第 2
1] 1]
標本におけるハプロタイプ番号 jのハプロタイプ度数 Xの実現値 Xを全てのハプロタ
2j 2j
イブ番号 jに対して設定する（図 12参照)。これにより、（X , · · · , X , X , · · · , X )
11 1 21 2 の組合せが与えられた。
[数 2]
M M
∑Χυ = 2Ν, , ∑X2j = 2N2 · · · (数式 2 ) ゾ=1 ゾ =1
[0024] ついで、ステップ SA— 4で設定した第 1標本に対応する複数のハプロタイプ度数お よびステップ SA— 3で定義した行列に基づ 、て第 1標本における各座位でのアレル 度数を算出し、ステップ SA— 4で設定した第 2標本に対応する複数のハプロタイプ度 数およびステップ SA— 3で定義した行列に基づ 、て第 2標本における各座位でのァ レル度数を算出する (ステップ SA— 5)。ここで、第 1標本における座位番号 iでのマイ ナーアレルのアレル度数 Yおよび第 2標本における座位番号 iでのマイナーアレル
li
のアレル度数 Yは、ランダム変数であり、それぞれ下記数式 3で表される。そして、下
2i
記数式 3に基づいて、第 1標本における座位番号 iでのマイナーアレルのアレル度数 Yの実現値 yおよび第 2標本における座位番号 iでのマイナーアレルのアレル度数 li li
Yの y
2i 実現値 2iを算出する。
[数 3]
M M
Υ2ί =∑α, χ Χ 2j 数式 3
/=1 ゾ = 数式 3において、 αはステップ SA— 3で定義した行列である。
[0025] っ 、で、ステップ SA— 5で算出した第 1標本における各座位でのアレル度数およ び第 2標本における各座位でのアレル度数に基づ 、て、各座位に対応する偶現表を 作成する (ステップ SA— 6)。具体的には、図 13に示す 2 X 2の偶現表を座位数分（ L個）だけ作成する。
[0026] ついで、ステップ SA— 6で作成した偶現表について、座位ごとに独立性検定を実 行する (ステップ SA— 7)。ここで、独立性検定について簡単に説明する。一般に、偶 現表に対して、下記数式 4で定義される Xを算出し、算出した Xの値が、自由度（1, 1)で有意水準 (例えば 0. 01や 0. 05など）の X (カイ）二乗分布値 χ a 2以上な
1,1
らば有意とする。
[数 4]
(数式 4 )
Figure imgf000016_0001
η 数式 4において、 iは標本を識別する番号である。ここでは、 i= lの場合はケース群 を示し、 i= 2の場合はコントロール群を示す。 jはアレルを識別する番号である。ここ では、 j = 1の場合はマイナーアレルを示し、 j = 2の場合はメジャーアレルを示す。 rij はアレル度数である。 nは標本番号 iにおけるマイナーアレル度数 (n )およびメジャ
i. il
一アレル度数 (n )の和である。 nはアレル番号 jに対応するケース群のアレル度数（
i2 .j
n )およびコントロール群のアレル度数（n )の和である。 nはケース群のサンプル数と lj 2j
コントロール群のサンプル数の和である。具体的には、 n=∑ 2∑ 2 nである。
i=l ij
[0027] っ 、で、ステップ SA— 7で実行した独立性検定の検定結果の 、ずれか 1つが有意 であった場合 (ステップ SA— 8 :Yes)、第 1標本の標本数、第 2標本の標本数、存在 し得るハプロタイプの総数、第 1標本のハプロタイプ度数、第 2標本のハプロタイプ度 数およびノヽプロタイプ頻度に基づ 、て、ステップ SA— 4で設定した第 1標本の各ハ プロタイプ度数と第 2標本の各ハプロタイプ度数との組合せに対応する確率を、予め 設定した第 1標本と第 2標本とのハプロタイプ度数に関する結合分布を用いて算出す る (ステップ SA— 9)。具体的には、第 1標本の標本数 N、第 2標本の標本数 N、存
1 2 在し得るハプロタイプの総数 M (2L)、第 1標本のハプロタイプ度数の実現値 X、第 2
1] 標本のハプロタイプ度数の実現値 Xおよびノ、プロタイプ頻度 hに基づいて、ステップ
2j j
SA— 4で設定した第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイプ度数と の組合せ (X , · · · , X , X , · · · , X )に対応する確率 Pを、下記数式 5に示す第 1
11 1 21 2
標本と第 2標本とのハプロタイプ度数に関する結合分布を用いて算出する。
[数 5]
± ( 1 1，ズ 12，… , Χ\Μ , X21， Χ22 ,…， Χ2Μ ) · · · (数式 5 )
Figure imgf000017_0001
数式 5において、 kは標本を識別する番号である。ここでは、 k= lの場合は第 1標 本を示し、 k= 2の場合は第 2標本を示す。
[0028] ついで、第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイプ度数との全て の組合せが終了したカゝ否かを判定し、全ての組合せが終了した場合 (ステップ SA— 10 : Yes)、ステップ SA— 4〜ステップ SA— 9を繰り返すごとにステップ SA— 9で算 出された確率を全て加算し、加算した確率を最終的に第 1種の過誤の確率とする (ス テツプ SA—11)。一方、全ての組合せが終了していない場合 (ステップ SA— 10 : N o)、ステップ SA— 4〜ステップ SA— 9を実行し、全ての組合せが終了するまで繰り 返す。
以上、説明したように、本発明にかかる確率算出方法は、（1)複数の SNPで構成さ れるハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を取 得し、 (2)取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に 含まれるハプロタイプ情報に基づ!/、て、存在し得るハプロタイプごとにハプロタイプ頻 度を算出し、 (3)各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを 有するかに従って予め定めた所定値を設定することで、当該所定値を成分とする行 列を定義し、（4)第 1標本における各ハプロタイプのハプロタイプ度数および第 2標 本における各ハプロタイプのハプロタイプ度数を設定することで、第 1標本の各ハプ 口タイプ度数と第 2標本の各ハプロタイプ度数との組合せを与え、（5)設定した第 1標 本に対応する複数のハプロタイプ度数および定義した行列に基づいて第 1標本にお ける各座位でのアレル度数を算出し、設定した第 2標本に対応する複数のハプロタイ プ度数および定義した行列に基づいて第 2標本における各座位でのアレル度数を算 出し、（6)算出した第 1標本における各座位でのアレル度数および第 2標本における 各座位でのアレル度数に基づいて、各座位に対応する偶現表を作成し、（7)作成し た偶現表について、座位ごとに独立性検定を実行し、（8)実行した独立性検定の検 定結果のいずれか 1つが有意であった場合、第 1標本の標本数、第 2標本の標本数 、存在し得るハプロタイプの総数、第 1標本のハプロタイプ度数、第 2標本のハプロタ イブ度数およびノヽプロタイプ頻度に基づ 、て、設定した第 1標本の各ハプロタイプ度 数と第 2標本の各ハプロタイプ度数との組合せに対応する確率を、予め設定した第 1 標本と第 2標本とのハプロタイプ度数に関する結合分布を用いて算出する。そして、 ( 4)、 （5)、 （6)、 （7)および (8)を、第 1標本の各ハプロタイプ度数と第 2標本の各ハプ 口タイプ度数との全ての組合せについて繰り返し、繰り返すごとに（8)で算出される確 率を加算することで、最終的に第 1種の過誤の確率を算出する。これにより、第 1種の 過誤の確率を、ハプロタイプ頻度を考慮して直接的に算出することができ、その結果 、疾患関連 SNPの見落としを防ぐことができる。
[0030] ここで、ステップ SA— 5において、数式 3に示す第 1標本における座位番号 iでのマ イナ一アレルのアレル度数 Yおよび第 2標本における座位番号 iでのマイナーアレル
li
のアレル度数 Y は、座位番号 iでのマイナーアレルのアレル頻度 pに関するハプロタ
2i i
イブ頻度 hの関数（下記数式 6参照）に基づ 、て定められる。
[数 6]
M
Pi = α·. χ /ζ ,. · · · (数式 6 )
[0031] また、ステップ SA— 9において、数式 5に示す第 1標本と第 2標本とのハプロタイプ 度数に関する結合分布は、下記数式 7に示す第 1標本のハプロタイプ度数に関する 結合分布と下記数式 8に示す第 2標本のハプロタイプ度数に関する結合分布との積 で与えられる。ここで、 case— control association研究などから得られた観察デー タを用いて検定を行う場合、帰無仮説 (H )として 2つの標本が同一の集団から得ら
0
れたことを仮定する。この 2つの標本 (第 1標本および第 2標本）におけるハプロタイプ 番号 jのハプロタイプ度数をそれぞれランダム変数 X、Xとすると、 2つの標本抽出は
lj 2j
独立なので、ランダム変数 Xおよび Xはそれぞれ二項分布に従う。よって、 X
lj 2j 11、 X
12
、 · · ·、 X の結合分布 (joint distribution)は多項分布であり、 X
1 21、 X
22、 · · ·、 X
2 の結合分布はそれとは独立な多項分布である。従って、ランダム変数 X
ljおよび Xの
2j 確率は、それぞれ下記数式 7および下記数式 8で示される。ここで、 2つの標本はい ずれも、同じハプロタイプ頻度 (複数)を持つ集団からのものであると仮定して 、ること に注意する（帰無仮説)。また、下記数式 7および下記数式 8において、 Xおよび X lj 2j は、それぞれランダム変数 Xおよび Xの実現値である。
lj 2j
[数 7]
Figure imgf000020_0001
Figure imgf000020_0002
また、集団から、第 1標本として 2N個、第 2標本として 2N個のハプロタイプコピー
1 2
を抽出する実験において、ハプロタイプコピーとハプロタイプとの概念の違いは次の 通りである。例えば、 1つの個体に 2つの同じノヽプロタイプ（ホモ接合体）が存在する ならば、この個体は 1つのハプロタイプと 2つのハプロタイプコピーを保有すると考える また、ステップ SA— 11で算出した第 1種の過誤の確率は、下記数式 9で表現する ことができる。そして、下記数式 9で求められる確率は、帰無仮説のもとでの確率であ り、ステップ SA— 7で実行した独立性検定の検定結果のいずれ力 1つが有意となる 確率を意味する。
[数 9]
Figure imgf000021_0001
ゾ =1 k=\
(数式 9 ) 数式 9において、 Zは、ステップ SA— 7で実行した独立性検定の検定結果のいず れカ 1つが有意であった場合には 1を、そうでな力つた場合には 0をとるように定義し たランダム変数である。変数 Zは上述したランダム変数 Yの関数であり、さらに当該ラ ki
ンダム変数 Yは上述したランダム変数 Xの関数であるので、変数 Zは関数 f (x , · · ki kj 11
•，x ，x ， · · · , X )で表すことができる。
1 21 2
[0034] ここで、アレル頻度モデルで第 1種の過誤の確率を算出する確率算出方法におい て、主要なハプロタイプによる保守的な確率算出方法について説明する。
上述では M (2^;Lは SNP座位数)個全てのハプロタイプについて計算を行ったが 、 M (2L)個の全てのハプロタイプが存在するとは限らない。その場合、ある程度以上 の累積頻度 (例えば 95%)に達するまで頻度の高いハプロタイプを数え、その個数が sだとすると、残りのハプロタイプ（以下では、「その他ハプロタイプ」と記載する場合が ある）を全て s + 1番目のハプロタイプとして取り扱うことが適切であると思われる。ただ し、その他ハプロタイプとして算出される第 1種の過誤の確率の誤差を考慮するべき である。
[0035] 具体的には、主要なハプロタイプによる保守的な確率算出方法は以下の手順で行 われる,
(1)累積頻度の高いハプロタイプから加え、 s個のハプロタイプを選択する。ここで、 累積度数を Pとする。その他ハプロタイプの頻度は 1—pである。また、 i番目のハプロ タイプ頻度を h (i= l s)とし、 h を 1—∑ とする。
i 〜 s+1 i=l i (2) 1〜5番目のハプロタイプについて、 i番目のハプロタイプの k番目の座位がマイ ナーアレルであれば 1、メジャーアレルであれば 0となる変数 a (i= l〜s, k= 1〜L
ik
)を定義する。 { a }は s XLの行列となる。
(3)f( X , X • · ·, X , X , X ', χ
Is 21 22 )なる関数を以下のように定義する。†dl、
2s
し、 x.ま iu = ）番目のハプロタイプの第 j (j = l, 2)標本におけるハプロタイプコ ピーの数に対応するランダム変数 Xの実現値である。
まず、任意の k(k=l〜L)について、図 19に示す偶現表を作成する。ここで、図 19 において、 bは、第 j標本の第 s- L番目のハプロタイプのハプロタイプコピー X 個の js+l うち、 k番座位についてマイナーアレルである割合である。これは、結果の第 s +1番 目のハプロタイプの X 個の実際のパプ口タイプの構成により変化する値であり、ラン js+l
ダム変数 Bの実現値と考えられる。
jk
ついで、作成した偶現表について、座位ごとに独立性検定を実行する。 ついで、独立性検定の結果のいずれかが有意になった場合は 1、それ以外の場合 は 0となる関数 fを定義する。
[0036] ここで、下記数式 10に示す標本ォッズ比は、 b =1および b =0の時に最大で、 b
lk 2k ]
=0および b =1の時に最小である。なぜなら、下記数式 10の標本ォッズ比の分子 k 2k
は b について単調増カロ、分母は b について単調減少だ力もである c
lk 2一-k
[数 10]
(数式 1 o )
Figure imgf000022_0001
[0037] ここで、表 1に示す偶現表があるとする。
[表 1]
(表 1)
a m— a
b 但し、表 1において、 a、 bは非負の整数で、 m≥a、 n≥bとする。
[0038] このとき、 Pearsonの χ 2統計量は下記数式 11のようになる。また、 Τを aで偏微分す ると下記数式 12となり、 bで偏微分すると下記数式 13となる。
[ w
(数式 1 1
Figure imgf000023_0001
[数 12]
(数式 1 2 )
Figure imgf000023_0002
[数 13]
Figure imgf000023_0003
[0039] 数式 12は、 a<bmZnで負、 a>bmZnで正である。なお、数式 12および数式 13 の m— aおよび n—bは同時に 0になることはないものとする。従って、 Tは aく bm/n で単調減少、 a≥bmZnで単調増加である。 bが固定された時、 a≤a≤aの範囲で
0 1
Tの最大値を考えると、 bに関わらず a = aまたは a = aで最大値をとる。同様にして、
0 1
数式 13は、 b< anZmで負、 b〉anZmで正である。すなわち、 Tは b< anZmで単 調減少、 b〉anZmで単調増加である。従って、 Tは aが固定された時、 b≤b≤bの
0 1 範囲では b =bまたは b =bで最大値をとる。
0 1
[0040] こうして、数式 11は、 a≤a≤aおよび b≤b≤bの範囲では、（a二 a ， b二 b )、（a
0 1 0 1 0 0
=a , b二 b )、（a二 a , b二 b )、（a二 a ， b二 b )のいずれかの点で最大値をとる。な
0 1 1 0 1 1
ぜなら、それ以外の点力 に最大値を与えるとすると、それは aが aでも aでもないか
0 1
、 bが bでも bでもないことを意味する力 前者の場合は aを aまたは aにすることによ
0 1 0 1 り更に大きい τを得ることに矛盾し、後者の場合は bを bまたは bにすることにより更に
0 1
大きい Tを得ることに矛盾する力もである。ただし、 a≤a≤aおよび b≤b≤bで、 a =b = 0でなく、 a=m, b =nでな!/、とする。従って、図 19の偶現表では、 (b =0, b lk 2k
= 0)、（b =0, b = 1)、（b = 1, b =0)、（b = 1, b = 1)のいずれかの時に最 lk 2k lk 2k lk 2k
大の％ 2統計量を示す。
[0041] 改めて主要なハプロタイプによる保守的な確率算出方法を具体的に述べると、 s + 1番目のハプロタイプの k番座位のアレルが第 1標本、第 2標本のそれぞれにおいて 、全てマイナーアレルまたは全てメジャーアレルとして検定を行い、いずれかの検定 で有意となれば k番座位について有意と判定する方法である。第 1標本、第 2標本の それぞれについて、 s + 1番ハプロタイプの全ての k番座位をマイナーアレルまたはメ ジャーアレルであると考えることにより、図 20から図 23に示す 4種類の偶現表ができ る。この内、 1つでも有意と判定される場合、図 20から図 23に示す 4つの偶現表のい ずれかが有意と判定される。すなわち、図 20から図 23に示す 4つの偶現表のいずれ かが有意であれば、関数 f= lとする。
[0042] つぎに、アレル頻度モデルで第 1種の過誤の確率を算出する確率算出方法におい て、上述したその他ハプロタイプのそれぞれの座位のアレル数のかわりに、頻度に応 じてサンプルして第 1種の過誤の確率を算出する確率算出方法について説明する。 上述した保守的な確率計算方法において、 p
s (主要ハプロタイプの累積頻度）は出 来る限り大きくするべきである力 そうすると主要ハプロタイプの頻度の推定値の信頼 性は低くなるであろう。ここで述べる確率算出方法は、第 k座位に対して、第 1標本、 第 2標本のそれぞれにつ!/、て、その他ハプロタイプのマイナーアレル数の代わりに、 頻度によりサンプルした数を用いる方法である。全てのハプロタイプの頻度が既知で あるとすれば、主要ハプロタイプの第 k座位でのマイナーアレル累積頻度は∑ s a i=l ik hであり、その他ハプロタイプの第 k座位でのマイナーアレル累積頻度は∑ L a h i i=s+l ik i である。しかし、一般に、 s + l、 · · ·、 L座位でのハプロタイプ頻度の信頼性は低い可 能性がある。一方、一般に、集団における第 k座位でのマイナーアレル頻度は既知 であることが多い。この集団における k座位でのマイナーアレル頻度を qとすると、そ k
の他ハプロタイプにおける第 k座位のマイナーアレル頻度は、 q—∑ s « hと考えら k i=l ik i れる。第 j標本のその他ハプロタイプのコピー数は X
js+lである。ゆえに、 B (X , q - js+1 k
∑ s a h )の二項分布に従い、標本ごと、座位ごとに第 k座位のマイナーアレル数を i=l ik i サンプルする。その実現値を第 k座位、第 j標本についてランダム変数 Zとすると、図
jk
24に示す偶現表を得る。そして、全ての座位について図 24に示す偶現表を作成し、 作成した偶現表について、座位ごとに独立性検定を行った後、そのいずれかが有意 であれば関数 f= lとし、それ以外は 0とする。
[0043] これまで述べてきた確率算出方法では、各 SNP座位での有意水準 (第 1種の過誤 の確率）が与えられた時、全体の第 1種の過誤の確率を算出した。そこで、ここでは、 全体の有意水準 (例えば P = 0. 05)が与えられた時に、望ましい有意水準を各 SNP 座位での有意水準を算出する確率算出方法について説明する。具体的には、以下 の手順で行う。
(1)各 SNPでの検定のための有意水準を、例えば P =0. 001、 0. 002、 0. 00
locus
3、 · · ·、 0. 010のように少しずつ増やしていき、上述した確率算出方法または後述 する近似的に確率を計算する方法により全体の第 1種の過誤の確率 P
allを算出する。
(2) P と P との関係により、望ましい P の値 (例えば 0. 05)に対応する P を計 locus all all locus 算する。なお、計算方法は直線回帰や曲線回帰などの方法を用いればよい。
(3)計算した各 SNPの検定における望ま U、有意水準 P を用いて各 SNPの検
locus
定を行う。これにより、連鎖不平衡を考慮した多数の SNP座位での検定が可能にな る。
[0044] 以上、アレル頻度モデルで第 1種の過誤の確率を算出する確率算出方法に関する 説明を終了する。
[0045] (1 2)優性'劣性遺伝モデル
図 2は、優性 ·劣性遺伝モデルで第 1種の過誤の確率を正確に算出する確率算出 方法の一例を示すフローチャートである。
まず、複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む 第 1標本および第 2標本を取得する (ステップ SB— 1)。なお、第 1標本および第 2標 本の取得に関する説明は、上述した（1 1)アレル頻度モデルと同様である。
[0046] っ 、で、ステップ SB— 1で取得した第 1標本に含まれるハプロタイプ情報および取 得した第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごと にハプロタイプ頻度を算出する (ステップ SB— 2)。なお、ハプロタイプ頻度の算出に 関する説明は、上述した（1— 1)アレル頻度モデルと同様である。
[0047] ついで、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有する 力に従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定 義する (ステップ SB— 3)。なお、行列の定義に関する説明は、上述した（1 1)ァレ ル頻度モデルと同様である。
[0048] ついで、第 1標本における各ディプロタイプ形のディプロタイプ形度数および第 2標 本における各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1標本の 各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与える (ス テツプ SB— 4)。ここで、第 1標本におけるディプロタイプ形番号 jj' (j'はハプロタイプ 番号であり、 l≤j'≤Mを満たす整数、ディプロタイプ形番号については図 14参照） のディプロタイプ形度数 X および第 2標本におけるディプロタイプ形番号 jj 'のデイブ ΐϋ'
口タイプ形度数 X は、ランダム変数であり、それぞれ下記数式 14を満たす。そして、
2jj'
下記数式 14を満たすように、第 1標本におけるディプロタイプ形番号 jj'のデイブロタ イブ形度数の実現値 X を全てのディプロタイプ形番号 jj 'に対して設定し、第 2標本
ljj
におけるディプロタイプ形番号 jj 'のディプロタイプ形度数の実現値 X を全てのデイブ
2jj，
口タイプ形番号 jj 'に対して設定する（図 15参照)。これにより、（X , · · · , X , χ ,
111 1 211
· · · , χ )の組合せが与えられた。
2
[数 14]
M M M M
∑∑^jf = Nx , ∑∑X2jf = N2 · · · (数式 1 4 ) ゾ =1 /=1 ゾ =1 /=1
[0049] ついで、ステップ SB— 4で設定した第 1標本に対応する複数のディプロタイプ形度 数およびステップ SB— 3で定義した行列に基づ 、て第 1標本における各座位でのァ レル度数を算出し、ステップ SB— 4で設定した第 2標本に対応する複数のデイブロタ イブ形度数およびステップ SB— 3で定義した行列に基づ 、て第 2標本における各座 位でのアレル度数を算出する (ステップ SB— 5)。ここで、第 1標本における座位 iでの マイナーアレルのアレル度数 Yおよび第 2標本における座位 iでのマイナーアレルの
li
アレル度数 Yは、ランダム変数であり、それぞれ下記数式 15で表される。そして、下 記数式 15に基づいて、第 1標本における座位 iでのマイナーアレルのアレル度数の 実現値 yおよび第 2標本における座位 iでのマイナーアレルのアレル度数の実現値 y li
2iを算出する。
[数 15]
M M M M
=∑∑" ' X = "リ' X · · · （数式 1 5 )
1 /=\ j=\ f=\ 数式 15において、 および はステップ SB— 3で定義した行列である。
[0050] ここで、劣性遺伝モデルの場合では、数式 15を用いてマイナーアレル（aa+ Aa)の アレル度数を算出するが、優性遺伝モデルの場合では下記数式 16を用いてマイナ アレル（aa)のアレル度数を算出する。
[数 16]
M M M M
=∑∑(1一 ) ( ') χ υ/ =∑∑(1一 ) (卜""') X '
=1 =1 =1 =1
• · · (数式 1 6 )
[0051] っ 、で、ステップ SB— 5で算出した第 1標本における各座位でのアレル度数および 第 2標本における各座位でのアレル度数に基づ ヽて、各座位に対応する偶現表を作 成する (ステップ SB— 6)。具体的には、劣性遺伝モデルの場合では数式 15を用い て算出したアレル度数に関する図 16に示す 2 X 2の偶現表を座位数分 (L個）だけ作 成し、優性遺伝モデルの場合では数式 16を用いて算出したアレル度数に関する図 1 7に示す 2 X 2の偶現表を座位数分 ( S)だけ作成する。
[0052] ついで、ステップ SB— 6で作成した偶現表について、座位ごとに独立性検定を実 行する (ステップ SB— 7)。
[0053] っ 、で、ステップ SB— 7で実行した独立性検定の検定結果の 、ずれか 1つが有意 であった場合 (ステップ SB— 8 :Yes)、第 1標本の標本数、第 2標本の標本数、存在 し得るハプロタイプの総数、第 1標本のディプロタイプ形度数、第 2標本のデイブロタ イブ形度数およびハプロタイプ頻度に基づ 、て、ステップ SB— 4で設定した第 1標本 の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せに対応す る確率を、予め設定した第 1標本と第 2標本とのディプロタイプ形度数に関する結合 分布を用いて算出する (ステップ SB— 9)。具体的には、第 1標本の標本数 N、第 2
1 標本の標本数 N、存在し得るハプロタイプの総数 M (2^)、第 1標本のディプロタイプ
2
形度数の実現値 X 、第 2標本のディプロタイプ形度数の実現値 X およびハプロタイ ltf 2jj'
プ頻度 hに基づいて、ステップ SB— 4で設定した第 1標本の各ディプロタイプ形度数 と第 2標本の各ディプロタイプ形度数との組合せ (X , · · · , X , X , · · · , X )に
111 1 211 2 対応する確率 Pを、下記数式 17に示す第 1標本と第 2標本とのディプロタイプ形度数 に関する結合分布を用いて算出する。
[数 17]
X\MM -> X2W Χ212，…， Χ2ΜΜノ · · · (数式 1 7 )
Figure imgf000028_0001
数式 17において、 kは標本を識別する番号である。ここでは、 k= lの場合は第 1標 本を示し、 k= 2の場合は第 2標本を示す。
[0054] ついで、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数と の全ての組合せが終了したカゝ否かを判定し、全ての組合せが終了した場合 (ステツ プ SB - 10 : Yes)、ステップ SB - 4〜ステップ SB - 9を繰り返すごとにステップ SB - 9で算出される確率を加算し、加算した確率を最終的に第 1種の過誤の確率とする（ ステップ SB— 11)。一方、全ての組合せが終了していない場合 (ステップ SB— 10 : No)、ステップ SB— 4〜ステップ SB— 9を実行し、全ての組合せが終了するまで繰り 返す。
[0055] 以上、説明したように、本発明にかかる確率算出方法は、（1)複数の SNPで構成さ れるハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を取 得し、 (2)取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に 含まれるハプロタイプ情報に基づ!/、て、存在し得るハプロタイプごとにハプロタイプ頻 度を算出し、 (3)各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを 有するかに従って予め定めた所定値を設定することで、当該所定値を成分とする行 列を定義し、（4)第 1標本における各ディプロタイプ形のディプロタイプ形度数および 第 2標本における各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1 標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与 え、 (5)設定した第 1標本に対応する複数のディプロタイプ形度数および定義した行 列に基づいて第 1標本における各座位でのアレル度数を算出し、設定した第 2標本 に対応する複数のディプロタイプ形度数および定義した行列に基づいて第 2標本に おける各座位でのアレル度数を算出し、（6)算出した第 1標本における各座位でのァ レル度数および第 2標本における各座位でのアレル度数に基づ 、て、各座位に対応 する偶現表を作成し、（7)作成した偶現表について、座位ごとに独立性検定を実行 し、（8)実行した独立性検定の検定結果のいずれ力 1つが有意であった場合、第 1 標本の標本数、第 2標本の標本数、存在し得るハプロタイプの総数、第 1標本のディ プロタイプ形度数、第 2標本のディプロタイプ形度数およびノヽプロタイプ頻度に基づ V、て、設定した第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度 数との組合せに対応する確率を、予め設定した第 1標本と第 2標本とのディプロタイ プ形度数に関する結合分布を用いて算出する。そして、（4)、 （5)、 （6)、 （7)および（ 8)を、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との全 ての組合せについて繰り返し、繰り返すごとに（8)で算出される確率を加算すること で、最終的に第 1種の過誤の確率を算出する。これにより、第 1種の過誤の確率を、 ハプロタイプ頻度を考慮して直接的に算出することができ、疾患関連 SNPの見落とし を防ぐことができる。
また、ステップ SB— 9において、数式 17に示す第 1標本と第 2標本とのディプロタイ プ形度数に関する結合分布は、全ての j、 j'について X の結合分布、すなわちテ
kjj'
ンソル {X }の分布は多項分布であるので、下記数式 18に示す第 1標本のディプロ
kjj'
タイプ形度数に関する結合分布と下記数式 19に示す第 2標本のディプロタイプ形度 数に関する結合分布との積で与えられる。ここで、優性遺伝モデルまたは劣性遺伝 モデルによる検定について第 1の過誤の確率を計算するためにはディプロタイプ形 について考える必要がある。なお、優性遺伝モデルによる検定とはメジャーアレルを 保有する個体とそうでな 、個体の比率を 2群で比較する検定法を 、、劣性遺伝モ デルによる検定とはマイナーアレルを保有する個体とそうでない個体の比率を 2群で 比較する検定法を 、う。 2つのハプロタイプコピーの順列を順列ディプロタイプ形と ヽ うことにし、ハプロタイプ番号 jのハプロタイプとハプロタイプ番号 j'のハプロタイプとに よる（この順番の)順列ディプロタイプ形を Dとする。一般に、 D =D は成り立たず、
JJ JJ J J
j=j 'であってもよい。ハプロタイプの総数を Mとすれば、順列ディプロタイプ形の総数 は M2である。ハプロタイプに関する Hardy— Weinberg平衡が成り立つているとする と、集団内の Dの頻度は hhである。今、第 1標本として N人、第 2標本として N人を
Ώ ί ϊ 1 2 ランダムに抽出すると、第 k標本について Dを持つ人の総数は二項分布 B(N , hh
JJ k j j
)に従うランダム変数である。第 k標本について、このランダム変数を X とする。特定
kjj'
の kについて、全ての X の結合分布、すなわち行列 {X }および {X }の分布は、
kjj' ljj' 2jj'
多項分布であり、その確率は、第 1標本および第 2標本のそれぞれについて、下記数 式 18および下記数式 19で与えられる。
[数 18]
^ ^\\\ ― XU ->^\ 2 ― Χ\\2·>'"·>^-ΙΜΜ = Χ\ΜΜ)
Μ Μ 、
YlYl{hjhfY" · · · (数式 18)
ゾ=1 =1
Figure imgf000030_0001
[数 19] (^211 = 211, A 212 ~ Χ2\2·> "' ·>^2ΜΜ = Χ2ΜΜ .
Μ Μ
Ν2、.
Μ Μ Ylfl{ jhfY"' · · · (数式 19)
1
Figure imgf000030_0002
また、ステップ SB— 11で算出した第 1種の過誤の確率は、下記数式 20で表現する ことができる。そして、下記数式 20で求められる確率は、ステップ SB— 7で実行した 独立性検定の検定結果のいずれか 1つが有意となる確率を意味する。
[数 20] p[z = 1]
Figure imgf000031_0001
=0 =0
∑ ·· ∑
Figure imgf000031_0002
Figure imgf000031_0003
：数式 2 0 ) 数式 20において、 Zは、ステップ SB— 7で実行した独立性検定の検定結果のいず れカ 1つが有意であった場合は 1、そうでな力つた場合は 0と定義したランダム変数で ある。すなわち、変数 Zは上述したランダム変数 Yの関数であり、さらに当該ランダム 変数 Yは上述したランダム変数 X の関数であるので、変数 Zは関数 f(x , · · ·, X
, x ··· χ )と表すことができる。また、数式 20において、 s (k=l, 2;j, j'=
1 M)は、下記数式 21で定義される。
[数 21]
Figure imgf000032_0001
Sk2\ ~ Sk\M ― Xk\M , Sk21 ― Sk2\ ― Xk2\ ' " * ' Sk2M一 Sk1M - \ ― ^k2M-\
SkM\ ― SkM~\ M ~ kM-\M
SkM2一 SkM\ ― XkM\
SkMM-\ ― Sk履- 2 ― 謝 - 2
• · · (数式 2 1 )
[0058] このように、優性遺伝モデルや劣性遺伝モデルで検定を行う場合、ハプロタイプコ ピーではなぐディプロタイプ形コピーの数を数える必要がある。しかし、ディプロタイ プ形コピーの数カもハプロタイプコピーの数を数えることは可能である。このことは重 要であり、ランダム変数 X の実現値 X を用いて、優性遺伝モデルや劣性遺伝モデ
kjj kjj
ルでの検定だけでなぐアレル頻度モデルの検定を行うことが可能だ力もである。そ して、連鎖した多数の SNP座位について、優性遺伝モデル、劣性遺伝モデルおよび アレル頻度モデルで検定を行い、そのいずれかが有意の場合は、有意とした場合の 全体の第 1種の過誤の確率を正確に計算することができる。具体的には、順列デイブ 口タイプ形コピーの数であるランダム変数 X の実現値 X を用いて、図 25に示す偶
kjj kjj
現表を作成すればよい。そして、全ての座位について独立性検定を行い、いずれか が有意の場合は、上述した関数 f (x , · · · , X , X , · · · , X )を 1とし、それ以外
111 1 211 2
の場合は 0とすればよい。この検定で算出した第 1の過誤の確率は、上述した（1 1 )アレル頻度モデルで算出した第 1の過誤の確率と一致する。
[0059] さらに、優性遺伝モデル、劣性遺伝モデルおよびアレル頻度モデルの 3つの検定 を行い、いずれかで有意の場合に有意と判定する方法による第 1種の過誤の確率も 、正確に計算することが可能である。具体的には、図 13、図 16および図 17の偶現表 につ 、て各座位のそれぞれで独立性検定を行 、、その 、ずれかで有意の場合は関 数 fを 1とし、それ以外では 0とすればよい。
[0060] 以上、優性'劣性遺伝モデルで第 1種の過誤の確率を算出する確率算出方法に関 する説明を終了する。
[0061] (1 3)遺伝子型モデル
図 3は、遺伝子型モデルで第 1種の過誤の確率を正確に算出する確率算出方法の 一例を示すフローチャートである。
まず、複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む 第 1標本および第 2標本を取得する (ステップ SC— 1)。なお、第 1標本および第 2標 本の取得に関する説明は、上述した（1 1)アレル頻度モデル、（1 2)優性'劣性 遺伝モデルと同様である。
[0062] っ 、で、ステップ SC— 1で取得した第 1標本に含まれるハプロタイプ情報および第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとにハプロ タイプ頻度を算出する (ステップ SC— 2)。なお、ハプロタイプ頻度の算出に関する説 明は、上述した（1 1)アレル頻度モデル、（1 2)優性 ·劣性遺伝モデルと同様で ある。
[0063] ついで、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有する 力に従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定 義する (ステップ SC— 3)。なお、行列の定義に関する説明は、上述した（1— 1)ァレ ル頻度モデル、 (1 - 2)優性 ·劣性遺伝モデルと同様である。
[0064] ついで、第 1標本における各ディプロタイプ形のディプロタイプ形度数および第 2標 本における各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1標本の 各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与える (ス テツプ SC— 4)。なお、ディプロタイプ形度数の設定に関する説明は、上述した（1 2)優性 ·劣性遺伝モデルと同様である。
[0065] ついで、ステップ SC— 4で設定した第 1標本に対応する複数のディプロタイプ形度 数およびステップ SC— 3で定義した行列に基づ 、て第 1標本における各座位での遺 伝子型度数を算出し、ステップ SC— 4で設定した第 2標本に対応する複数のデイブ 口タイプ形度数およびステップ SC— 3で定義した行列に基づ 、て第 2標本における 各座位での遺伝子型度数を算出する (ステップ SC— 5)。ここで、第 1標本における 座位 iでのメジャーホモ (AA)の遺伝子型度数 Yおよび第 2標本における座位 iでの メジャーホモ (AA)の遺伝子型度数 Y、第 1標本における座位 iでのマイナーホモ（a
2i
a)の遺伝子型度数 Y 'および第 2標本における座位 iでのマイナーホモ (aa)の遺伝
li
子型度数 Y 'ならびに第 1標本における座位 iでのへテロ (Aa)の遺伝子型度数 Y "
2i li および第 2標本における座位 iでのへテロ (Aa)の遺伝子型度数 Y "は、それぞれラ
2i
ンダム変数であり、それぞれ下記数式 22で表される。そして、下記数式 22に基づい て、遺伝子型度数 Yの
li 実現値 y
liおよび遺伝子型度数 Yの
2i 実現値 y
2i、遺伝子型度 数 γ 'の実現値 y 'および遺伝子型度数 Y 'の実現値 y 'ならびに遺伝子型度数 Y " li li 2i 2i li の実現値 y
li "および遺伝子型度数 Y "の
2i 実現値 y
2i "を算出する。
[数 22]
YM
Figure imgf000034_0001
, M M r M M
Y =∑∑ ( " ( " ， 4 =∑∑(l -aIJ )(l - alf ) X2j/ Y ' - ^ - (Υ, + Υ, ' ) ， 4" = N2— (
• . · (数式 2 2 )
[0066] っ 、で、ステップ SC— 5で算出した第 1標本における各座位での遺伝子型度数お よび第 2標本における各座位での遺伝子型度数に基づ ヽて、各座位に対応する偶 現表を作成する (ステップ SC— 6)。具体的には、数式 22を用いて算出した遺伝子 型度数に関する図 18に示す 2 X 3の偶現表を座位数分 (L個）だけ作成する。
[0067] ついで、ステップ SC— 6で作成した偶現表について、座位ごとに独立性検定を実 行する (ステップ SC— 7)。
[0068] っ 、で、ステップ SC— 7で実行した独立性検定の検定結果の 、ずれか 1つが有意 であった場合 (ステップ SC— 8 :Yes)、第 1標本の標本数、第 2標本の標本数、存在 し得るハプロタイプの総数、第 1標本のディプロタイプ形度数、第 2標本のデイブロタ イブ形度数およびハプロタイプ頻度に基づ 、て、ステップ SC— 4で設定した第 1標本 の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せに対応す る確率を、予め設定した第 1標本と第 2標本とのディプロタイプ形度数に関する結合 分布を用いて算出する (ステップ SC— 9)。なお、確率の算出に関する説明は、上述 した（ 1 2)優性 ·劣性遺伝モデルと同様である。
[0069] ついで、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数と の全ての組合せが終了したカゝ否かを判定し、全ての組合せが終了した場合 (ステツ プ SC - 10 : Yes)、ステップ SC - 4〜ステップ SC - 9を繰り返すごとにステップ SC - 9で算出される確率を加算し、加算した確率を最終的に第 1種の過誤の確率とする（ ステップ SC— 11)。一方、全ての組合せが終了していない場合 (ステップ SC— 10 : No)、ステップ SC— 4〜ステップ SC— 9を実行し、全ての組合せが終了するまで繰り 返す。
[0070] 以上、説明したように、本発明にかかる確率算出方法は、（1)複数の SNPで構成さ れるハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を取 得し、 (2)取得した第 1標本に含まれるハプロタイプ情報および取得した第 2標本に 含まれるハプロタイプ情報に基づ!/、て、存在し得るハプロタイプごとにハプロタイプ頻 度を算出し、 (3)各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを 有するかに従って予め定めた所定値を設定することで、当該所定値を成分とする行 列を定義し、（4)第 1標本における各ディプロタイプ形のディプロタイプ形度数および 第 2標本における各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1 標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与 え、 (5)設定した第 1標本に対応する複数のディプロタイプ形度数および定義した行 列に基づいて第 1標本における各座位での遺伝子型度数を算出し、設定した第 2標 本に対応する複数のディプロタイプ形度数および定義した行列に基づいて第 2標本 における各座位での遺伝子型度数を算出し、（6)算出した第 1標本における各座位 での遺伝子型度数および第 2標本における各座位での遺伝子型度数に基づいて、 各座位に対応する偶現表を作成し、（7)作成した偶現表について、座位ごとに独立 性検定を実行し、 (8)実行した独立性検定の検定結果のいずれか 1つが有意であつ た場合、第 1標本の標本数、第 2標本の標本数、存在し得るハプロタイプの総数、第 1標本のディプロタイプ形度数、第 2標本のディプロタイプ形度数およびノ、プロタイプ 頻度に基づいて、設定した第 1標本の各ディプロタイプ形度数と第 2標本の各デイブ 口タイプ形度数との組合せに対応する確率を、予め設定した第 1標本と第 2標本との ディプロタイプ形度数に関する結合分布を用いて算出する。そして、（4)、（5)、（6)、 (7)および (8)を、第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ 形度数との全ての組合せについて繰り返し、繰り返すごとに（8)で算出される確率を 加算することで、最終的に第 1種の過誤の確率を算出する。これにより、第 1種の過 誤の確率を、ハプロタイプ頻度を考慮して直接的に算出することができ、その結果、 疾患関連 SNPの見落としを防ぐことができる。
[0071] 以上、遺伝子型モデルで第 1種の過誤の確率を算出する確率算出方法に関する 説明を終了する。
[0072] (1 -4)正確なハプロタイプ頻度がわ力 な 、場合
図 4は、正確なハプロタイプ頻度がわ力 な 、場合に第 1種の過誤の確率を正確に 算出する確率算出方法の一例を示すフローチャートである。
まず、複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む 第 1標本および第 2標本を取得する (ステップ SD— 1)。なお、第 1標本および第 2標 本の取得に関する説明は、上述した（1 1)アレル頻度モデル、（1 2)優性'劣性 遺伝モデル、 (1 - 3)遺伝子型モデルと同様である。
[0073] ついで、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有する 力に従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定 義する (ステップ SD— 2)。なお、行列の定義に関する説明は、上述した（1 1)ァレ ル頻度モデル、（1 2)優性.劣性遺伝モデル、（1 3)遺伝子型モデルと同様であ る。
[0074] ついで、第 1標本における各ハプロタイプのハプロタイプ度数および第 2標本にお ける各ハプロタイプのハプロタイプ度数を設定することで、第 1標本の各ハプロタイプ 度数と第 2標本の各ハプロタイプ度数との組合せを与える (ステップ SD— 3)。ここで 、第 1標本におけるハプロタイプ番号 jのハプロタイプ度数 Xおよび第 2標本における
1]
ハプロタイプ番号 jのハプロタイプ度数 Xは、ランダム変数であり、それぞれ下記数式
2j
23を満たす。そして、下記数式 23を満たすように、第 1標本におけるハプロタイプ番 号 jのハプロタイプ度数 Xの実現値 Xを全てのハプロタイプ番号 jに対して設定し、第
1] 1]
2標本におけるハプロタイプ番号 jのハプロタイプ度数 Xの実現値 Xを全てのハプロ タイプ番号 jに対して設定する。これにより、（X , · · · , X , X , · · · , X )の組合せが
11 1 21 2
与えられた。
[数 23]
2 M M
∑Xkj = X, ， ∑X
Figure imgf000037_0001
, ∑X2j = N2 · · · (数式 2 3 ) k=\ 7=1 ゾ =1
[0075] っ 、で、ステップ SD— 3で設定した第 1標本に対応する複数のハプロタイプ度数お よびステップ SD— 2で定義した行列に基づ 、て第 1標本における各座位でのアレル 度数を算出し、ステップ SD— 3で設定した第 2標本に対応する複数のハプロタイプ 度数およびステップ SD— 2で定義した行列に基づ ヽて第 2標本における各座位での アレル度数を算出する (ステップ SD—4)。なお、アレル度数の算出に関する説明は 、上述した（1— 1)アレル頻度モデルと同様である。
[0076] ついで、ステップ SD— 4で算出した第 1標本における各座位でのアレル度数およ び第 2標本における各座位でのアレル度数に基づ 、て、各座位に対応する偶現表を 作成する (ステップ SD— 5)。
[0077] っ 、で、ステップ SD— 5で作成した偶現表にっ 、て、座位ごとに独立性検定を実 行する（ステップ SD— 6)。
[0078] ついで、ステップ SD— 6で実行した独立性検定の検定結果のいずれか 1つが有意 であった場合 (ステップ SD— 7 : Yes)、第 1標本の標本数、第 2標本の標本数、存在 し得るハプロタイプの総数、第 1標本のハプロタイプ度数および第 2標本のハプロタイ プ度数に基づ 、て、ステップ SD— 3で設定した第 1標本の各ハプロタイプ度数と第 2 標本の各ハプロタイプ度数との組合せに対応する確率を、予め設定した第 1標本と 第 2標本とのハプロタイプ度数に関する結合分布を用いて算出する (ステップ SD— 8 )。具体的には、第 1標本の標本数 N、第 2標本の標本数 N、存在し得るハプロタイ
1 2
プの総数 M (2L)、第 1標本のハプロタイプ度数の実現値 Xおよび第 2標本のハプロ
1]
タイプ度数の実現値 X に基づいて、ステップ SD— 3で設定した第 1標本の各ハプロ
2j
タイプ度数と第 2標本の各ハプロタイプ度数との組合せ (X , · · · , X , X , · · · , X )
11 1 21 2 に対応する確率 Pを、下記数式 24に示す第 1標本と第 2標本とのハプロタイプ度数に 関する結合分布を用いて算出する。なお、下記数式 24に示す第 1標本と第 2標本と のハプロタイプ度数に関する結合分布は、超幾何分布である。
[数 24]
1，2，·'·，Μ
(数式 2 4 )
Figure imgf000038_0001
=1 k=l 数式 24において、 kは標本を識別する番号である。ここでは、 k= lの場合は第 1標 本を示し、 k = 2の場合は第 2標本を示す。
[0079] ついで、第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイプ度数との全て の組合せが終了したカゝ否かを判定し、全ての組合せが終了した場合 (ステップ SD— 9 : Yes)、ステップ SD— 3〜ステップ SD— 8を繰り返すごとにステップ SD— 8で算出 された確率を全て加算し、加算した確率を最終的に第 1種の過誤の確率とする (ステ ップ SD— 10)。一方、全ての組合せが終了していない場合 (ステップ SD— 9 : No)、 ステップ SD— 3〜ステップ SD— 8を実行して、全ての組合せが終了するまで繰り返 す。
[0080] 以上、説明したように、本発明にかかる確率算出方法は、（1)複数の SNPで構成さ れるハプロタイプに関するハプロタイプ情報を複数含む第 1標本および第 2標本を取 得し、 (2)各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するか に従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定義 し、（3)第 1標本における各ハプロタイプのハプロタイプ度数および第 2標本における 各ハプロタイプのハプロタイプ度数を設定することで、第 1標本の各ハプロタイプ度数 と第 2標本の各ハプロタイプ度数との組合せを与え、（4)設定した第 1標本に対応す る複数のハプロタイプ度数および定義した行列に基づいて第 1標本における各座位 でのアレル度数を算出し、設定した第 2標本に対応する複数のハプロタイプ度数およ び定義した行列に基づ 、て第 2標本における各座位でのアレル度数を算出し、 (5) 算出した第 1標本における各座位でのアレル度数および第 2標本における各座位で のアレル度数に基づいて、各座位に対応する偶現表を作成し、（6)作成した偶現表 について、座位ごとに独立性検定を実行し、（7)実行した独立性検定の検定結果の いずれか 1つが有意であった場合、第 1標本の標本数、第 2標本の標本数、存在し得 るハプロタイプの総数、第 1標本のハプロタイプ度数および第 2標本のハプロタイプ 度数に基づいて、設定した第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイ プ度数との組合せに対応する確率を、予め設定した第 1標本と第 2標本とのハプロタ イブ度数に関する結合分布を用いて算出する。そして、（3)、（4)、（5)、（6)および（ 7)を、第 1標本の各ハプロタイプ度数と第 2標本の各ハプロタイプ度数との全ての組 合せについて繰り返し、繰り返すごとに（7)で算出される確率を加算することで、最終 的に第 1種の過誤の確率を算出する。これにより、第 1種の過誤の確率を、直接的に 算出することができ、その結果、疾患関連 SNPの見落としを防ぐことができる。
また、ステップ SD— 10で算出した第 1種の過誤の確率は、下記数式 25で表現する ことができる。そして、下記数式 25で求められる確率は、ステップ SD— 6で実行した 独立性検定の検定結果のいずれか 1つが有意となる確率を意味する。
[数 25]
Figure imgf000039_0001
数式 2 5 ) 数式 25において、 Zは、ステップ SD— 6で実行した独立性検定の検定結果のいず れカ 1つが有意であった場合には 1を、そうでな力つた場合には 0をとるように定義し たランダム変数である。変数 Zは上述したランダム変数 Yの関数であり、さらに当該ラ
ki
ンダム変数 Yは上述したランダム変数 Xの関数である。また、ランダム変数 Xにつ
ki kj kj いては数式 23の制限がある。よって、変数 Zは関数 f (x , · · · , X )で表すことができ
11 1
る。なお、その他の Xはその他の変数の関数として表される。
kj
[0082] 以上、正確なハプロタイプ頻度がわ力 ない場合に第 1種の過誤の確率を正確に 算出する確率算出方法に関する説明を終了する。
[0083] [第 2実施形態]
ここでは、第 1種の過誤の確率を MCMC法（マルコフ連鎖モンテカルロ法）を用い て近似的に算出する本発明にかかる確率算出方法について、（2—1)アレル頻度モ デル、（2— 2)優性 ·劣性遺伝モデル、 （2— 3)遺伝子型モデル、 （2— 4)正確なハプ 口タイプ頻度がわ力もない場合、の順に図を参照して詳細に説明する。
[0084] (2— 1)アレル頻度モデル
図 5は、アレル頻度モデルで第 1種の過誤の確率を MCMC法を用 V、て近似的に 算出する確率算出方法の一例を示すフローチャートである。なお、全てのハプロタイ プ頻度が既知であることを前提とする。
まず、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに 従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定義す る (ステップ SE— 1)。なお、行列の定義に関する説明は、上述した第 1実施形態に おける（1— 1)アレル頻度モデルと同様である。
[0085] ついで、ハプロタイプ度数の組合せを予め設定した第 1標本および第 2標本に対し て生成する (ステップ SE— 2)。具体的には、第 1標本におけるハプロタイプ番号 jの ハプロタイプ度数 Xの実現値 Xおよび第 2標本におけるハプロタイプ番号 jのハプロ
1] 1]
タイプ度数 Xの実現値 Xを、上述した数式 2を満たすように生成する。これにより、（X
2j 2j
, · · · , χ , χ , · · · , X )の組合せを生成する。
11 1 21 2
[0086] っ 、で、第 1標本または第 2標本の 、ずれかの標本を選択する (ステップ SE— 3)。
具体的には、標本番号 k(k= l, 2)の標本を等確率で選択する。
[0087] つ!、で、第 1ハプロタイプを選択する（ステップ SE— 4)。具体的には、ハプロタイプ 番号 j (l≤j≤M ( = 2L)； Lは座位数)を選択する。
[0088] つ!、で、ステップ SE— 3で選択した標本における第 1ハプロタイプのハプロタイプ度 数が 0でない場合（SE— 5 : Yes)、当該第 1ハプロタイプとは異なる第 2ハプロタイプ を選択する (ステップ SE— 6)。具体的には、選択した標本番号 kの標本におけるハ プロタイプ番号 jのハプロタイプ度数 X ≠0 (または X >0)の場合、ステップ SE— 3で
kj kj
選択したハプロタイプ番号 jとは異なるハプロタイプ番号 j ' ( 1≤ j '≤ M)を選択する。な お、ハプロタイプ度数 X =0の場合、ハプロタイプ度数 Xの値を維持する。
kj kj
[0089] ついで、予め設定した第 1ハプロタイプのハプロタイプ頻度および第 2ハプロタイプ のハプロタイプ頻度に基づいて、第 1ハプロタイプのハプロタイプ度数および第 2ハプ 口タイプのハプロタイプ度数を更新する (ステップ SE— 7)。具体的には、まず、下記 数式 26で定義される cの値を算出する。なお、下記数式 26において、 hおよび h,は、 それぞれノヽプロタイプ番号 jのハプロタイプ頻度およびハプロタイプ番号 j'のハプロタ イブ頻度であり、その値は予め設定されている。
[数 26] c = hf X" · · · (数式 2 6 )
/2バ¾, + 1) ついで、 c≥lの場合、 Xの値から 1を減算し、 X の値に 1をカ卩算する。また、 c < lの
kj kj'
場合、 cの確率で Xの値から 1を減算し、 X の値に 1を加算し、 1— cの確率で Xおよ
kj kj び X の値を維持する。
kj'
[0090] っ 、で、ステップ SE— 7で更新した後の第 1標本に対応する複数のハプロタイプ度 数およびステップ SE— 1で定義した行列に基づ 、て第 1標本における各座位でのァ レル度数を算出し、ステップ SE— 7で更新した後の第 2標本に対応する複数のハプ 口タイプ度数およびステップ SE— 1で定義した行列に基づ 、て第 2標本における各 座位でのアレル度数を算出する (ステップ SE— 8)。なお、アレル度数の算出に関す る説明は、上述した第 1実施形態の（1— 1)アレル頻度モデルと同様である。
[0091] ついで、ステップ SE— 8で算出した第 1標本における各座位でのアレル度数および 第 2標本における各座位でのアレル度数に基づ ヽて、各座位に対応する偶現表を作 成する（ステップ SE— 9)。
[0092] っ 、で、ステップ SE— 9で作成した偶現表にっ 、て、座位ごとに独立性検定を実 行する（ステップ SE— 10)。 [0093] ついで、ステップ SE— 10で実行した独立性検定の検定結果のいずれか 1つが有 意であった場合 (ステップ SE— 11： Yes)、当該場合の累積回数を更新する (ステツ プ SE— 12)。具体的には、累積回数に 1を加算する。
[0094] ついで、所定回数が終了した力否かを判定し、所定回数が終了した場合 (ステップ SE—13 :Yes)、所定回数に対する最終的な累積回数の割合 (累積回数 ÷所定回 数)を第 1種の過誤の確率として算出する (ステップ SE— 14)。一方、所定回数が終 了して ヽな 、場合 (ステップ SE— 13： No)、ステップ SE - 3〜ステップ SE - 12を実 行し、所定回数が終了するまで繰り返す。
[0095] 以上、説明したように、本発明に力かる確率算出方法は、 Metropolis— Hastings 法により、標本ごとに期待される多項分布に従うハプロタイプのサンプルを生成し、そ れぞれのサンプルにつ!/、て各座位で検定を行 、、その 、ずれかの座位で有意となる 割合を求める。なお、 Markov— chainの状態は、ハプロタイプ度数の実現値による 2 X M (存在し得るハプロタイプの総数)の偶現表により代表される。
[0096] ここで、主要ハプロタイプを用いた計算では、主要ハプロタイプに加え、その他ハプ 口タイプを設ける。従って、マルコフ連鎖（Markov— chain)の状態は X (k = l , 2 ;j kj
= 1, 2, · · · , s)の 2標本の主要ハプロタイプコピーの数とその他ハプロタイプコピー の数 X により代表される。従って、 X (k= l , 2 ;j = l , 2, · · · , s+ 1)の全てのハプ ks+l kj
口タイプを用いた場合、 x (k= l, 2 ;j = l, 2, · · · , M ( = 2L) )と同様に MCMC法を
kj
実行すればよい。なお、検定については、上述したように保守的な検定法またはその 他ハプロタイプコピー数をサンプルすることにより行う。第 1種の過誤の確率は、同様 に有意と判定されたステップの割合を算出すればよい。
[0097] また、多項分布に従う Monte Carloサンブラにより第 1標本、第 2標本を独立に生 成し、上述した関数 fの計算を行い、これが 1となる割合を第 1種の過誤の確率とする 方法が考えられる。
[0098] 以上、アレル頻度モデルで第 1種の過誤の確率を MCMC法を用いて近似的に算 出する確率算出方法に関する説明を終了する。
[0099] (2— 2)優性 ·劣性遺伝モデル
図 6は、優性 ·劣性遺伝モデルで第 1種の過誤の確率を MCMC法を用いて近似的 に算出する確率算出方法の一例を示すフローチャートである。なお、全てのハプロタ イブ頻度が既知であることを前提とする。
まず、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに 従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定義す る (ステップ SF— 1)。なお、行列の定義に関する説明は、上述した第 1実施形態に おける（1— 2)優性 ·劣性遺伝モデルと同様である。
[0100] っ 、で、ディプロタイプ形度数の組合せを予め設定した第 1標本および第 2標本に 対して生成する (ステップ SF— 2)。具体的には、第 1標本におけるディプロタイプ形 番号 ϋ 'のディプロタイプ形度数の実現値 X を全てのディプロタイプ形番号 ϋ 'に対し ΐϋ'
て設定し、第 2標本におけるディプロタイプ形番号 jj'のディプロタイプ形度数の実現 値 X を、上述した数式 14を満たすように生成する。これにより、（X , · · · , X , X
111 1 211
, · · · , X )の組合せを生成する。
2
[0101] ついで、第 1標本または第 2標本のいずれかの標本を選択する (ステップ SF— 3)。
具体的には、標本番号 k(k= l, 2)の標本を等確率で選択する。
[0102] っ 、で、第 1ディプロタイプ形を選択する (ステップ SF—4)。具体的には、ディプロ タイプ形番号 '（1≤順序つきの自然数 iおよび j^M ^ S1^)； Lは座位数)を選択 する。
[0103] っ 、で、ステップ SF— 3で選択した標本における第 1ディプロタイプ形のディプロタ イブ形度数が 0でな 、場合 (SF- 5 :Yes)、当該第 1ディプロタイプ形とは異なる第 2 ディプロタイプ形を選択する (ステップ SF— 6)。具体的には、選択した標本番号 kの 標本におけるディプロタイプ形番号 j j 'のディプロタイプ形度数 x 半 0 (または X
1 1 kjljl' kjljl'
>0)の場合、ステップ SF— 3で選択したディプロタイプ形番号 j j 'とは異なるディプロ
1 1
タイプ形番号 j j ' (1≤順序つきの自然数 jおよび j '≤M ( = 2L)； Lは座位数)を選択
2 2 2 2
する。なお、ディプロタイプ形度数 X =0の場合、ディプロタイプ形度数 X の値を kjljl' kjljl' 維持する。
[0104] ついで、予め設定した第 1ディプロタイプ形を構成する各ハプロタイプのハプロタイ プ頻度および第 2ディプロタイプ形を構成する各ハプロタイプのハプロタイプ頻度に 基づ 、て、第 1ディプロタイプ形のディプロタイプ形度数および第 2ディプロタイプ形 のディプロタイプ形度数を更新する (ステップ SF— 7)。具体的には、まず、下記数式 27で定義される cの値を算出する。なお、下記数式 27において、 h、 h 、 hおよび h
jl jl' ]2 は、それぞれハプロタイプ番号 jのハプロタイプ頻度、ハプロタイプ番号 j 'のハプロ
)2' 1 1 タイプ頻度、ハプロタイプ番号 jのハプロタイプ頻度およびハプロタイプ番号 j 'のハプ
2 2 口タイプ頻度であり、その値は予め設定されている。
[数 27]
c = ゾ 2 · · · (数式 2 7 )
hj、h八 X , + 1)
A J\ 2ゾ 2 ついで、 c≥lの場合、 X の値から 1を減算し、 X の値に 1を加算する。また、 c <
kjljl' kj2j2'
1の場合、 cの確率で x の値から 1を減算し、 X の値に 1をカ卩算し、 l—cの確率で
kjljl' kj2j2'
X および X の値を維持する。
kjljl' kj2j2'
[0105] っ 、で、ステップ SF— 7で更新した後の第 1標本に対応する複数のディプロタイプ 形度数およびステップ SF— 1で定義した行列に基づ 、て第 1標本における各座位で のアレル度数を算出し、ステップ SF— 7で更新した後の第 2標本に対応する複数の ディプロタイプ形度数およびステップ SF— 1で定義した行列に基づ!/、て第 2標本に おける各座位でのアレル度数を算出する (ステップ SF— 8)。なお、アレル度数の算 出に関する説明は、上述した第 1実施形態の（1 2)優性'劣性遺伝モデルと同様で ある。
[0106] ついで、ステップ SF— 8で算出した第 1標本における各座位でのアレル度数および 第 2標本における各座位でのアレル度数に基づ ヽて、各座位に対応する偶現表を作 成する（ステップ SF— 9)。
[0107] っ 、で、ステップ SF— 9で作成した偶現表にっ 、て、座位ごとに独立性検定を実 行する（ステップ SF— 10)。
[0108] ついで、ステップ SF— 10で実行した独立性検定の検定結果のいずれ力 1つが有 意であった場合 (ステップ SF— 11： Yes)、当該場合の累積回数を更新する (ステツ プ SF— 12)。具体的には、累積回数に 1を加算する。
[0109] ついで、所定回数が終了した力否かを判定し、所定回数が終了した場合 (ステップ SF—13 :Yes)、所定回数に対する最終的な累積回数の割合 (累積回数 ÷所定回 数)を第 1種の過誤の確率として算出する (ステップ SF— 14)。一方、所定回数が終 了して ヽな 、場合 (ステップ SF— 13： No)、ステップ SF - 3〜ステップ SF— 12を実 行し、所定回数が終了するまで繰り返す。
[0110] 以上、説明したように、本発明に力かる確率算出方法は、 Metropolis— Hastings 法により、標本ごとに期待される多項分布に従うディプロタイプ形のサンプルを生成し 、それぞれのサンプルについて各座位で検定を行い、そのいずれかの座位で有意と なる割合を求める。
[0111] 以上、優性'劣性遺伝モデルで第 1種の過誤の確率を MCMC法を用いて近似的 に算出する確率算出方法に関する説明を終了する。
[0112] (2— 3)遺伝子型モデル
図 7は、遺伝子型モデルで第 1種の過誤の確率を MCMC法を用いて近似的に算 出する確率算出方法の一例を示すフローチャートである。なお、全てのハプロタイプ 頻度が既知であることを前提とする。
まず、各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに 従って予め定めた所定値を設定することで、当該所定値を成分とする行列を定義す る (ステップ SG— 1)。なお、行列の定義に関する説明は、上述した第 1実施形態に おける（1— 3)遺伝子型モデルと同様である。
[0113] ついで、ディプロタイプ形度数の組合せを予め設定した第 1標本および第 2標本に 対して生成する (ステップ SG— 2)。なお、ディプロタイプ形度数の生成に関する説明 は、上述した（2— 2)優性 ·劣性遺伝モデルと同様である。
[0114] ついで、第 1標本または第 2標本のいずれかの標本を選択する (ステップ SG— 3)。
具体的には、標本番号 k(k= l, 2)の標本を等確率で選択する。
[0115] ついで、第 1ディプロタイプ形を選択する (ステップ SG—4)。具体的には、ディプロ タイプ形番号 j j ' (1≤順序つきの自然数 jおよび j '≤M ( = 2L)； Lは座位数)を選択
1 1 1 1
する。
[0116] ついで、ステップ SG— 3で選択した標本における第 1ディプロタイプ形のディプロタ イブ形度数が 0でな 、場合 (SG- 5 :Yes)、当該第 1ディプロタイプ形とは異なる第 2 ディプロタイプ形を選択する (ステップ SG— 6)。なお、第 2ディプロタイプ形の選択に 関する説明は、上述した (2— 2)優性 ·劣性遺伝モデルと同様である。
[0117] ついで、予め設定した第 1ディプロタイプ形を構成する各ハプロタイプのハプロタイ プ頻度および第 2ディプロタイプ形を構成する各ハプロタイプのハプロタイプ頻度に 基づ 、て、第 1ディプロタイプ形のディプロタイプ形度数および第 2ディプロタイプ形 のディプロタイプ形度数を更新する (ステップ SG 7)。なお、ディプロタイプ形度数 の更新に関する説明は、上述した (2— 2)優性 ·劣性遺伝モデルと同様である。
[0118] ついで、ステップ SG— 7で更新した後の第 1標本に対応する複数のディプロタイプ 形度数およびステップ SG— 1で定義した行列に基づ 、て第 1標本における各座位で の遺伝子型度数を算出し、ステップ SG— 7で更新した後の第 2標本に対応する複数 のディプロタイプ形度数およびステップ SG— 1で定義した行列に基づ 、て第 2標本 における各座位での遺伝子型度数を算出する (ステップ SG— 8)。なお、遺伝子型度 数の算出に関する説明は、上述した第 1実施形態の（1— 3)遺伝子型モデルと同様 である。
[0119] ついで、ステップ SG— 8で算出した第 1標本における各座位での遺伝子型度数お よび第 2標本における各座位での遺伝子型度数に基づ ヽて、各座位に対応する偶 現表を作成する (ステップ SG— 9)。
[0120] っ 、で、ステップ SG— 9で作成した偶現表にっ 、て、座位ごとに独立性検定を実 行する（ステップ SG— 10)。
[0121] ついで、ステップ SG— 10で実行した独立性検定の検定結果のいずれか 1つが有 意であった場合 (ステップ SG— 11： Yes)、当該場合の累積回数を更新する (ステツ プ SG— 12)。具体的には、累積回数に 1を加算する。
[0122] ついで、所定回数が終了した力否かを判定し、所定回数が終了した場合 (ステップ SG—13 :Yes)、所定回数に対する最終的な累積回数の割合 (累積回数 ÷所定回 数)を第 1種の過誤の確率として算出する (ステップ SG— 14)。一方、所定回数が終 了して 、な 、場合 (ステップ SG— 13： No)、ステップ SG— 3〜ステップ SG - 12を実 行し、所定回数が終了するまで繰り返す。
[0123] 以上、説明したように、本発明に力かる確率算出方法は、 Metropolis— Hastings 法により、標本ごとに期待される多項分布に従うディプロタイプ形のサンプルを生成し 、それぞれのサンプルについて各座位で検定を行い、そのいずれかの座位で有意と なる割合を求める。
[0124] 以上、遺伝子型モデルで第 1種の過誤の確率を MCMC法を用いて近似的に算出 する確率算出方法に関する説明を終了する。
[0125] (2—4)正確なハプロタイプ頻度がわ力もない場合
正確なハプロタイプ頻度がわ力 ない場合に第 1種の過誤の確率を MCMC法を用 いて近似的に算出する確率算出方法では、上述した数式 24の超幾何分布に従うノ、 プロタイプのサンプルを生成し、それぞれのサンプルにつ ヽて各座位での検定を行 い、そのいずれかの座位で有意となる割合、つまり上述した関数 fが 1となる割合を求 める。
[0126] ここで、一般に、真のハプロタイプ頻度がわからな 、ことも多 、。そこで、ハプロタイ プ頻度がわ力もなくても、標本内のディプロタイプ形がわ力 ている場合には、以下 のように計算することが可能である。
[0127] 順列ハプロタイプコピー数は、ハプロタイプ頻度が与えられれば多項分布に従う。こ こで、標本空間の確実な出来事の条件の下での確率ではなぐ特定の観察データの 下での確率を考える。例えば、下記数式 28または下記数式 29のように制限をカロえ、 それを満足する結果の集合という条件の下での確率を考える。
[数 28]
Figure imgf000047_0001
· · · (数式 2 8 )
ん=1
[数 29]
Figure imgf000048_0001
(数式 29
k=l
[0128] ここで、数式 29において、 Y は k番目の標本における順列ディプロタイプ形コピー
kij
数ではなぐ組合せディプロタイプ形コピー数をランダム変数で表したものである。一 般に、順列ディプロタイプ形を知ることは困難であり、 Y を用いる方が妥当であろう。
kij
また、一般に上述した優性遺伝モデル、劣性遺伝モデル、アレル頻度モデルおよび 遺伝子型モデルでの検定は、全て組合せディプロタイプ形コピー数が分かれば可能 である。
[0129] 上述した条件の下では、ランダム変数テンソル {X }(k=l, 2;i=l, 2, ···, M(
kij
= 2L) ;j = l, 2, ···, M( = 2L))または {Y }(k=l, 2;i=l, 2， ···， M( = 2L);1
kij
≤j≤i)は超幾何分布に従う。その確率はそれぞれ下記数式 30および下記数式 31 で示される。
[数 30]
P (数式 30)
Figure imgf000048_0002
'=l j=\ k=\
[数 31]
M M y.ij - M M 2 (数式 3 1)
Figure imgf000048_0003
)!]!!!!! !
r=l j=\ k=\ なお、これらの超幾何分布に従って第 1種の過誤の確率を正確に算出する確率算 出方法は、先に述べたように、多項分布に従って第 1種の過誤の確率を算出する確 率算出方法と同等である。
[0131] つぎに、冗長な計算を避けるため、第 1種の過誤の確率を MCMC法を用いて近似 的に算出する確率算出方法について説明する。
[0132] 実際の手順は以下のとおりである。
(A— 1)観察データ {Y }(k=l, 2;i=l, 2, ···, M( = 2L);l≤j≤i)が与えられ
kij
ている。
(A— 2) MCMC法により数式 31に示す超幾何分布に従うサンプルを生成する。伹 し、 {y }は下記数式 32を満たす。ここで、 MCMC法では、第 1標本の 1つの組合せ kij
ディプロタイプ形と第 2標本の 1つの組合せディプロタイプ形をランダムに選択し、 y
kij に 1つ加える、または減じるなどの方法で結果的に超幾何分布になるようにすればよ い。
[数 32] y, =∑ykJJ · · · (数式 3 2)
=1
(A— 3)生成されたサンプルについて、各座位のそれぞれで、優性遺伝モデル、劣 性遺伝モデル、アレル頻度モデルおよび遺伝子型モデルで検定を行う。
(A— 4)検定結果のいずれかで有意である場合に f(y ) =1となり、それ以外で 0と
kij
なる関数 fを定義する。
(A- 5) MCMC法で生成したサンプルのうち、 f = 1となる割合を経験的な第 1種の 過誤の確率とする。
[0133] さらに詳細に述べると、実際の手順は以下のとおりである。
(B— 1)任意の整数値である y (k=l, 2;i=l, 2, ···, M( = 2L)； l≤j≤i;Lは
kij
座位数） 、所定の制約のもとに与えられる。ここで、マルコフ連鎖（Markov— chain )の状態空間は、当該所定の制約のもと、互いに異なる値をとる複数の要素 y力 な
kij る集合 {y }である。なお、所定の制約とは、 "任意の iおよび j(i=l, 2, ···, M( = 2L
kij
);l≤j≤i;Lは座位数）に対して「∑ =y」を満たす"および"任意の k(k=l, 2
)に対して「∑
Figure imgf000049_0001
=n」を満たす"というものである。ただし、 yは、予め観測さ
=1 れた (与えられた）固定の非負整数値である。また、 nは、標本番号 kの標本に含まれ
k
る標本数である。
(B- 2)標本番号 k (k= 1 , 2)を等確率で選択し、選択した標本番号 kの値に応じ て標本番号 k' (k' = l, 2)をさらに選択する。具体的には、標本番号 kとして 1を選択 した場合には、標本番号 k'として 2を選択し、標本番号 kとして 2を選択した場合には 標本番号 k'として 1を選択する。
(B 3) 2つの順序つき整数である (u, V)を等確率で選択することにより、ディプロ タイプ形番号 uvを決定する。なお、 uは
Figure imgf000050_0001
； Lは座位数」を満たし、 V は「l≤v≤u」を満たす。
(B— 4)選択した標本番号 kの標本におけるディプロタイプ形番号 uvのディプロタイ プ形度数 y が「y >Ojを満たす場合には、選択した (u, V)とは異なる 2つの順序
kuv
つき整数である (W, S)をさらに選択することによりディプロタイプ形番号 WSをさらに決 定する。ここで、 wは「l≤w≤M ( = 2L)； Lは座位数」を満たし、 sは「l≤s≤w」を満 たす。一方、ディプロタイプ形度数 y 力^ y =0」を満たす場合には、ディプロタイプ
kuv kuv
形度数 y の値を維持して、後述する（B— 6)へ進み、そして上述した (B— 2)へ戻る kuv さらに、選択した標本番号 k'の標本におけるディプロタイプ形番号 wsのディプロタ イブ形度数 y , 力^ y , =0」を満たす場合には、ディプロタイプ形度数 y , の値を維
k ws k ws k ws
持して、後述する（B— 6)へ進み、そして上述した (B— 2)へ戻る。
(B—5)新たなディプロタイプ形度数 (y* 、 y 、 y , 、y , )を算出する。具体的
kuv kws k uv k ws
には、まず、下記数式 33で定義される cの値を算出する。なお、下記数式 33におい て、 y* は「y* =y 1」を満たし、 y* は「y* =y + 1
kws 」を満たし、 y*
k uvは「y*， kuv kuv kuv kws k uv
=y , + 1」を満たし、 y* =y , — 1
k uv k， wsは「y*，
k ws k ws 」を満たす。
[数 33]
C―
Figure imgf000050_0002
• · · (数式 3 3 ) ついで、 c≥lの場合、新たなディプロタイプ形度数として、 y 、y 、y , および y ,
kuv kws k uv k ws の代わりに y* 、y 、y , および y* , を用い、後述する（B— 6)へ進み、そして上述
kuv kws k uv k ws
した (B— 2)へ戻る。一方、 c< lの場合、 cの確率で、新たなディプロタイプ形度数と して、 y 、y 、y 、y* 、y*
kuv kws k uvおよび y の
k ws 代わりに y*
kuv kws k uvおよび y*
k wsを用い、 l—cの 確率で y 、y 、y , および y , の値を維持し、後述する（B— 6)へ進み、そして上述 kuv kws κ uv k ws
した (B— 2)へ戻る。
(B— 6)劣性遺伝モデル、優性遺伝モデル、アレル頻度モデル、遺伝子型モデル の各モデルで、 3つの偶現表のどれかを用いて、表現型と遺伝子型との間の独立性 検定を座位ごとに実行する。ここで、劣性遺伝モデルの場合には図 16の偶現表を用 い、優性遺伝モデルの場合には図 17の偶現表を用いる。ただし、図 16において、数 式 15の X および X の代わりにそれぞれ、 y および y を用いる。また、図 17におい
ltf 2jj' ljj' 2jj'
て、数式 16の x および X の代わりにそれぞれ、 y および y を用いる。また、遺伝
ltf 2jj' ltf 2jj'
子型モデルの場合には図 18の偶現表を用い、アレル頻度モデルの場合には図 25 の偶現表を用いる。ただし、図 18において、数式 22の X および X の代わりにそれ
ltf 2jj'
ぞれ、 y および y を用いる。また、図 25において、 x および x の代わりにそれぞ ltf ¾j' ltf 2jj'
れ、 y および y を用いる。
ltf ¾j'
(B- 7)実行した独立性検定の検定結果のいずれか 1つが有意であった場合、実 行した全体 (全ての座位にっ 、て）の独立性検定の検定結果が有意であると判定す る。
(B— 8)これまでの MCMC法による処理を所定回数分繰り返し実行し、（B— 7)で 有意であると判定された回数 (有意判定回数)の割合 (有意判定回数 ÷所定回数)を 、全体 (全ての座位について）の独立性検定における経験的な第 1種の過誤の確率 として算出する。
ここで、上述した手順では 3つのモデル (優性遺伝モデル、劣性遺伝モデル、遺伝 子型モデル)のどれかを用いた独立性検定について述べた力 本発明にカゝかる確率 算出方法は、 3つのモデル全てを用いた独立性検定に容易に拡張することができる 。したがって、上述した (B— 6)で述べたように、独立性検定は 3つの異なる偶現表を 用いて各座位にっ 、て行われ、 、ずれかのモデル · 、ずれかの座位での独立性検 定の検定結果が有意であるなら、全体 (全ての座位にっ 、て）の独立性検定の検定 結果を有意であるとする。
[0134] 以上、正確なハプロタイプ頻度がわ力 ない場合に第 1種の過誤の確率を MCMC 法を用いて近似的に算出する確率算出方法に関する説明を終了する。
産業上の利用可能性
[0135] 以上のように、本発明に力かる確率算出方法は、第 1種の過誤の確率を、ハプロタ イブ頻度を考慮して直接的に算出することができ、その結果、疾患関連 SNPの見落 としを防ぐことができる。そのため、本発明にかかる確率算出方法は、医療や創薬な どの分野において極めて有用である。

Claims

請求の範囲
SNPを用いた case— control相関解析で行う独立性検定における第 1種の過誤の 確率を算出する確率算出方法にぉ 、て、
複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む第 1標 本および第 2標本を取得する標本取得ステップと、
前記標本取得ステップで取得した第 1標本に含まれるハプロタイプ情報および取得 した第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとに ハプロタイプ頻度を算出するハプロタイプ頻度算出ステップと、
各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに従つ て予め定めた所定値を設定することで、当該所定値を成分とする行列を定義する行 列定義ステップと、
第 1標本における各ハプロタイプのハプロタイプ度数および第 2標本における各ハ プロタイプのハプロタイプ度数を設定することで、第 1標本の各ハプロタイプ度数と第
2標本の各ハプロタイプ度数との組合せを与えるハプロタイプ度数設定ステップと、 前記ハプロタイプ度数設定ステップで設定した第 1標本に対応する複数のハプロタ イブ度数および前記行列定義ステップで定義した行列に基づいて第 1標本における 各座位でのアレル度数を算出し、設定した第 2標本に対応する複数のハプロタイプ 度数および定義した行列に基づいて第 2標本における各座位でのアレル度数を算 出するアレル度数算出ステップと、
前記アレル度数算出ステップで算出した第 1標本における各座位でのアレル度数 および第 2標本における各座位でのアレル度数に基づ ヽて、各座位に対応する偶現 表を作成する偶現表作成ステップと、
前記偶現表作成ステップで作成した偶現表につ！ヽて、座位ごとに独立性検定を実 行する独立性検定実行ステップと、
前記独立性検定実行ステップで実行した独立性検定の検定結果のいずれか 1つ が有意であった場合、第 1標本の標本数、第 2標本の標本数、存在し得るハプロタイ プの総数、第 1標本のハプロタイプ度数、第 2標本のハプロタイプ度数およびノ、プロ タイプ頻度に基づ 、て、前記ハプロタイプ度数設定ステップで設定した第 1標本の各 ハプロタイプ度数と第 2標本の各ハプロタイプ度数との組合せに対応する確率を、予 め設定した第 1標本と第 2標本とのハプロタイプ度数に関する結合分布を用いて算出 する確率算出ステップと、
を含み、
前記ハプロタイプ度数設定ステップ、前記アレル度数算出ステップ、前記偶現表作 成ステップ、前記独立性検定実行ステップおよび前記確率算出ステップを、第 1標本 の各ハプロタイプ度数と第 2標本の各ハプロタイプ度数との全ての組合せについて繰 り返し、繰り返すごとに算出される前記確率を加算することで、最終的に第 1種の過 誤の確率を算出すること、
を特徴とする確率算出方法。
SNPを用いた case— control相関解析で行う独立性検定における第 1種の過誤の 確率を算出する確率算出方法にぉ 、て、
複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む第 1標 本および第 2標本を取得する標本取得ステップと、
前記標本取得ステップで取得した第 1標本に含まれるハプロタイプ情報および取得 した第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとに ハプロタイプ頻度を算出するハプロタイプ頻度算出ステップと、
各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに従つ て予め定めた所定値を設定することで、当該所定値を成分とする行列を定義する行 列定義ステップと、
第 1標本における各ディプロタイプ形のディプロタイプ形度数および第 2標本にお ける各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1標本の各ディ プロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与えるデイブロタ イブ形度数設定ステップと、
前記ディプロタイプ形度数設定ステップで設定した第 1標本に対応する複数のディ プロタイプ形度数および前記行列定義ステップで定義した行列に基づいて第 1標本 における各座位でのアレル度数を算出し、設定した第 2標本に対応する複数のディ プロタイプ形度数および定義した行列に基づいて第 2標本における各座位でのァレ ル度数を算出するアレル度数算出ステップと、
前記アレル度数算出ステップで算出した第 1標本における各座位でのアレル度数 および第 2標本における各座位でのアレル度数に基づ ヽて、各座位に対応する偶現 表を作成する偶現表作成ステップと、
前記偶現表作成ステップで作成した偶現表につ！ヽて、座位ごとに独立性検定を実 行する独立性検定実行ステップと、
前記独立性検定実行ステップで実行した独立性検定の検定結果のいずれか 1つ が有意であった場合、第 1標本の標本数、第 2標本の標本数、存在し得るハプロタイ プの総数、第 1標本のディプロタイプ形度数、第 2標本のディプロタイプ形度数および ハプロタイプ頻度に基づ ヽて、前記ディプロタイプ形度数設定ステップで設定した第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せに 対応する確率を、予め設定した第 1標本と第 2標本とのディプロタイプ形度数に関す る結合分布を用いて算出する確率算出ステップと、
を含み、
前記ディプロタイプ形度数設定ステップ、前記アレル度数算出ステップ、前記偶現 表作成ステップ、前記独立性検定実行ステップおよび前記確率算出ステップを、第 1 標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との全ての組合 せについて繰り返し、繰り返すごとに算出される前記確率を加算することで、最終的 に第 1種の過誤の確率を算出すること、
を特徴とする確率算出方法。
SNPを用いた case— control相関解析で行う独立性検定における第 1種の過誤の 確率を算出する確率算出方法にぉ 、て、
複数の SNPで構成されるハプロタイプに関するハプロタイプ情報を複数含む第 1標 本および第 2標本を取得する標本取得ステップと、
前記標本取得ステップで取得した第 1標本に含まれるハプロタイプ情報および取得 した第 2標本に含まれるハプロタイプ情報に基づ 、て、存在し得るハプロタイプごとに ハプロタイプ頻度を算出するハプロタイプ頻度算出ステップと、
各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに従つ て予め定めた所定値を設定することで、当該所定値を成分とする行列を定義する行 列定義ステップと、
第 1標本における各ディプロタイプ形のディプロタイプ形度数および第 2標本にお ける各ディプロタイプ形のディプロタイプ形度数を設定することで、第 1標本の各ディ プロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せを与えるデイブロタ イブ形度数設定ステップと、
前記ディプロタイプ形度数設定ステップで設定した第 1標本に対応する複数のディ プロタイプ形度数および前記行列定義ステップで定義した行列に基づいて第 1標本 における各座位での遺伝子型度数を算出し、設定した第 2標本に対応する複数のデ ィプロタイプ形度数および定義した行列に基づいて第 2標本における各座位での遺 伝子型度数を算出する遺伝子型度数算出ステップと、
前記遺伝子型度数算出ステップで算出した第 1標本における各座位での遺伝子型 度数および第 2標本における各座位での遺伝子型度数に基づ 、て、各座位に対応 する偶現表を作成する偶現表作成ステップと、
前記偶現表作成ステップで作成した偶現表につ！ヽて、座位ごとに独立性検定を実 行する独立性検定実行ステップと、
前記独立性検定実行ステップで実行した独立性検定の検定結果のいずれか 1つ が有意であった場合、第 1標本の標本数、第 2標本の標本数、存在し得るハプロタイ プの総数、第 1標本のディプロタイプ形度数、第 2標本のディプロタイプ形度数および ハプロタイプ頻度に基づ ヽて、前記ディプロタイプ形度数設定ステップで設定した第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との組合せに 対応する確率を、予め設定した第 1標本と第 2標本とのディプロタイプ形度数に関す る結合分布を用いて算出する確率算出ステップと、
を含み、
前記ディプロタイプ形度数設定ステップ、前記遺伝子型度数算出ステップ、前記偶 現表作成ステップ、前記独立性検定実行ステップおよび前記確率算出ステップを、 第 1標本の各ディプロタイプ形度数と第 2標本の各ディプロタイプ形度数との全ての 組合せについて繰り返し、繰り返すごとに算出される前記確率を加算することで、最 終的に第 1種の過誤の確率を算出すること、
を特徴とする確率算出方法。
SNPを用いた case— control相関解析で行う独立性検定における第 1種の過誤の 確率を MCMC法を用 、て近似的に算出する確率算出方法であって、
各ハプロタイプの各座位がマイナーアレルまたはメジャーアレルを有するかに従つ て予め定めた所定値を設定することで、当該所定値を成分とする行列を定義する行 列定義ステップと、
ハプロタイプ度数の組合せを予め設定した第 1標本および第 2標本に対して生成 するハプロタイプ度数生成ステップと、
第 1標本または第 2標本のいずれかの標本を選択する標本選択ステップと、 第 1ハプロタイプを選択する第 1ハプロタイプ選択ステップと、
前記標本選択ステップで選択した標本における前記第 1ハプロタイプのハプロタイ プ度数力^でない場合、当該第 1ハプロタイプとは異なる第 2ハプロタイプを選択する 第 2ハプロタイプ選択ステップと、
予め設定した前記第 1ハプロタイプのハプロタイプ頻度および前記第 2ハプロタイプ のハプロタイプ頻度に基づ 、て、前記第 1ハプロタイプのハプロタイプ度数および前 記第 2ハプロタイプのハプロタイプ度数を更新するハプロタイプ度数更新ステップと、 前記ハプロタイプ度数更新ステップで更新した後の第 1標本に対応する複数のハ プロタイプ度数および前記行列定義ステップで定義した行列に基づいて第 1標本に おける各座位でのアレル度数を算出し、更新した後の第 2標本に対応する複数のハ プロタイプ度数および定義した行列に基づいて第 2標本における各座位でのアレル 度数を算出するアレル度数算出ステップと、
前記アレル度数算出ステップで算出した第 1標本における各座位でのアレル度数 および第 2標本における各座位でのアレル度数に基づ ヽて、各座位に対応する偶現 表を作成する偶現表作成ステップと、
前記偶現表作成ステップで作成した偶現表につ！ヽて、座位ごとに独立性検定を実 行する独立性検定実行ステップと、
前記独立性検定実行ステップで実行した独立性検定の検定結果のいずれか 1つ が有意であった場合、当該場合の累積回数を更新する累積回数更新ステップと、 を含み、
前記標本選択ステップ、前記第 1ハプロタイプ選択ステップ、前記第 2ハプロタイプ 選択ステップ、前記ハプロタイプ度数更新ステップ、前記アレル度数算出ステップ、 前記偶現表作成ステップ、前記独立性検定実行ステップおよび前記累積回数更新 ステップを所定回数繰り返し、前記所定回数に対する最終的な累積回数の割合を第 1種の過誤の確率として算出すること、
を特徴とする確率算出方法。
PCT/JP2005/022317 2004-12-03 2005-12-05 確率算出方法 Ceased WO2006059770A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2006546685A JP4755601B2 (ja) 2004-12-03 2005-12-05 確率算出方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2004-351839 2004-12-03
JP2004351839 2004-12-03

Publications (1)

Publication Number Publication Date
WO2006059770A1 true WO2006059770A1 (ja) 2006-06-08

Family

ID=36565195

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2005/022317 Ceased WO2006059770A1 (ja) 2004-12-03 2005-12-05 確率算出方法

Country Status (2)

Country Link
JP (1) JP4755601B2 (ja)
WO (1) WO2006059770A1 (ja)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004192018A (ja) * 2002-10-16 2004-07-08 Japan Biological Informatics Consortium Ｄｎａプールによるハプロタイプ頻度推定方法
JP2004229511A (ja) * 2003-01-28 2004-08-19 Ntt Data Corp ハプロタイプ解析装置、および、ハプロタイプ解析方法をコンピュータに実行させることを特徴とするプログラム

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004192018A (ja) * 2002-10-16 2004-07-08 Japan Biological Informatics Consortium Ｄｎａプールによるハプロタイプ頻度推定方法
JP2004229511A (ja) * 2003-01-28 2004-08-19 Ntt Data Corp ハプロタイプ解析装置、および、ハプロタイプ解析方法をコンピュータに実行させることを特徴とするプログラム

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
KAMATANI N. ET AL: "Haplotype o Kiso ni shita Human Genome Takei to Hyogengata no Kanren no Kaiseki Shuho", JAPANESE FEDERATION OF STATISTICAL SCIENCE ASSOCIATIONS, vol. 2004, 6 September 2004 (2004-09-06), pages 188 - 189, XP002999646 *
NYHOLT D.R.: "A simple correction for multiple testing for single-nucleotide polymorphisms in linkage disequilibrium with each other", AM.J.HUM.GENET., vol. 74, 2004, pages 765 - 769, XP002999648 *
TOMITA M. ET AL: "Rensa Fuheiko Kaiseki ni okeru Tazai no Kankei", JAPANESE FEDERATION OF STATISTICAL SCIENCE ASSOCIATIONS, vol. 2003, 5 September 2003 (2003-09-05), pages 227 - 228, XP002999647 *
XIONG M. ET AL: "The haplotype linkage disequilibrium test for genome-wide screens: its power and study design", PACIFIC SYMPOSIUM ON BIOCOMPUTING, 2000, pages 672 - 681, XP002999649 *

Also Published As

Publication number Publication date
JP4755601B2 (ja) 2011-08-24
JPWO2006059770A1 (ja) 2008-06-05

Similar Documents

Publication Publication Date Title
Mathers et al. Chromosome-scale genome assemblies of aphids reveal extensively rearranged autosomes and long-term conservation of the X chromosome
Di Rienzo et al. Heterogeneity of microsatellite mutations within and between loci, and implications for human demographic histories
Liu et al. Ancient and modern genomes unravel the evolutionary history of the rhinoceros family
Burgess et al. Estimation of hominoid ancestral population sizes under Bayesian coalescent models incorporating mutation rate variation and sequencing errors
Vakhrusheva et al. Genomic signatures of recombination in a natural population of the bdelloid rotifer Adineta vaga
Fernández-Mazuecos et al. Resolving recent plant radiations: power and robustness of genotyping-by-sequencing
Thomas et al. Coding single-nucleotide polymorphisms associated with complex vs. Mendelian disease: evolutionary evidence for differences in molecular effects
Tsai et al. Population genomics of the wild yeast Saccharomyces paradoxus: quantifying the life cycle
Schaffner et al. Calibrating a coalescent simulation of human genome sequence variation
Blumenstiel et al. An age-of-allele test of neutrality for transposable element insertions
Morgan et al. Informatics resources for the Collaborative Cross and related mouse populations
Ding et al. Genome structure-based Juglandaceae phylogenies contradict alignment-based phylogenies and substitution rates vary with DNA repair genes
Whelan et al. Estimating the frequency of events that cause multiple-nucleotide changes
Long et al. Low base-substitution mutation rate in the germline genome of the ciliate Tetrahymena thermophila
Kim Allele frequency distribution under recurrent selective sweeps
Collet et al. Rapid evolution of the intersexual genetic correlation for fitness in Drosophila melanogaster
Menardo et al. Reconstructing the evolutionary history of powdery mildew lineages (Blumeria graminis) at different evolutionary time scales with NGS data
Barton et al. New methods for inferring the distribution of fitness effects for INDELs and SNPs
Mugal et al. Polymorphism data assist estimation of the nonsynonymous over synonymous fixation rate ratio ω for closely related species
Riebler et al. Bayesian variable selection for detecting adaptive genomic differences among populations
Ebler et al. Pangenome-based genome inference
Maruki et al. Evolutionary genomics of a subdivided species
WO2015006668A1 (en) Methods for identification of individuals
Fazalova et al. Low spontaneous mutation rate and pleistocene radiation of pea aphids
Cohen et al. The social supergene dates back to the speciation time of two Solenopsis fire ant species

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A1

Designated state(s): AE AG AL AM AT AU AZ BA BB BG BR BW BY BZ CA CH CN CO CR CU CZ DE DK DM DZ EC EE EG ES FI GB GD GE GH GM HR HU ID IL IN IS JP KE KG KM KN KP KR KZ LC LK LR LS LT LU LV LY MA MD MG MK MN MW MX MZ NA NG NI NO NZ OM PG PH PL PT RO RU SC SD SE SG SK SL SM SY TJ TM TN TR TT TZ UA UG US UZ VC VN YU ZA ZM ZW

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): GM KE LS MW MZ NA SD SL SZ TZ UG ZM ZW AM AZ BY KG KZ MD RU TJ TM AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HU IE IS IT LT LU LV MC NL PL PT RO SE SI SK TR BF BJ CF CG CI CM GA GN GQ GW ML MR NE SN TD TG

121 Ep: the epo has been informed by wipo that ep was designated in this application
WWE Wipo information: entry into national phase

Ref document number: 2006546685

Country of ref document: JP

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 05811804

Country of ref document: EP

Kind code of ref document: A1

WWW Wipo information: withdrawn in national office

Ref document number: 5811804

Country of ref document: EP