WO2017115741A1 - アリル特異的遺伝子発現頻度の推定方法、推定用コンピュータシステム及び推定用プログラム - Google Patents

アリル特異的遺伝子発現頻度の推定方法、推定用コンピュータシステム及び推定用プログラム Download PDF

Info

Publication number
WO2017115741A1
WO2017115741A1 PCT/JP2016/088640 JP2016088640W WO2017115741A1 WO 2017115741 A1 WO2017115741 A1 WO 2017115741A1 JP 2016088640 W JP2016088640 W JP 2016088640W WO 2017115741 A1 WO2017115741 A1 WO 2017115741A1
Authority
WO
WIPO (PCT)
Prior art keywords
allele
calculated
distribution
lead
value
Prior art date
Application number
PCT/JP2016/088640
Other languages
English (en)
French (fr)
Inventor
直樹 成相
正朗 長崎
要 小島
隆広 三森
洋介 河合
Original Assignee
国立大学法人東北大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立大学法人東北大学 filed Critical 国立大学法人東北大学
Publication of WO2017115741A1 publication Critical patent/WO2017115741A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/10—Gene or protein expression profiling; Expression-ratio estimation or normalisation

Definitions

  • the present invention relates to the field of genetic analysis based on nucleic acid information, and more specifically, based on RNA sequence (RNA-seq) analysis information derived from a high-performance sequencer, allele-specific gene expression frequency, in other words,
  • the invention relates to a means for estimating the gene expression level in each allele of a polyploid organism.
  • the polyploid organism is a normal diploid organism such as a mammal including a human
  • each allele is an allele derived from the paternal and an allele derived from the maternal.
  • Allele-specific expression corresponds to various expression modes of allele genes of polyploid organisms.
  • Non-patent Document 1 For example, in “genomic imprinting”, a phenomenon that depends on whether it is a paternal gene or a maternal gene, one of the two female X chromosomes has been studied. Inactivation of sex chromosomes that inactivate is also known to be based on ASE (Non-patent Document 1). Recent studies have shown that ASE is a relatively universal phenomenon (Non-Patent Document 2), and that many cis-acting mutations alter gene expression (Non-Patent Document 2). Patent Document 3). In some cases, differences in gene expression in the three alleles are associated with disease trends such as colorectal cancer, and importantly, the amount of gene transcription is a susceptibility locus for normal diseases such as diabetes and obesity. It has been reported that it can be used as a quantitative index for discovering (Non-patent Documents 4 and 5).
  • RNA sequence (RNA-seq) analysis is used as a method for examining allele-specific gene expression. This method uses read information of cDNA derived from a sample RNA fragmented using a high-performance sequencer, but several problems have been reported at present.
  • RNA-seq reads are mapped to a plurality of positions in a reference nucleotide sequence, causing a problem in quantitative expression of genes (Mortazavi et al., Nature Methods, 5, 621-628 (2008). )).
  • the specimen has a heterogeneous SNP (single base substitution) different from the reference base sequence, it allows the presence of bias during mapping (Degner et al., Bioinformatics, 25 3207-3212 (2009)).
  • a) the genome sequences of paternal and maternal alleles (alleles) in diploid organisms are determined, and b) the sequence portions (exons) that are expressed as mRNA from the genome sequences by splicing.
  • the present inventors have examined a means for estimating the allele-specific gene expression of polyploid organisms as a whole gene based on the nucleotide sequence of the RNA-seq lead of the specimen. went.
  • the nucleotide sequence of the lead of the RNA-seq is used as observation data, a specific latent element is virtually associated with the allele-specific gene expression frequency, and the frequency is statistically estimated and calculated by a computer.
  • the present invention was completed by finding that the allele-specific gene expression frequency of somatic organisms can be estimated very accurately.
  • the present invention provides the following steps (1) to (6) for data in which the read information of cDNA of K polyploid organisms (K is an integer of 2 or more) is mixed (hereinafter also referred to as “allele mixed data”). Is provided, and the method for estimating the allele-specific gene expression frequency by computer (hereinafter also referred to as “estimation method of the present invention”) is provided.
  • Step (3) The total expected mapping number calculated in step (2) is divided by the sum of the total expected mapping number for each isoform to calculate the isoform amount ratio, and (ii) each isoform In each allele, the allele content of the total expected mapping number of the isoform is divided by the total expected mapping number of the isoform to calculate the allele preference degree, (4) Step (1) is executed again for the new isoform ratio and allyl preference calculated in step (3), and a new expected mapping number is calculated, (5) Step (2) or (3) is executed again for the new expected number of mappings obtained in step (4), and the isoform ratio and the allyl preference level are newly calculated, (6) Steps (4) and (5) show that the difference between the number of expected mappings calculated in (i) step (4) and the number of expected mappings calculated in the previous step (4) is all Or (ii) isoform ratio and / or allyl preference calculated in step (5) and isoform ratio and / or ratio calculated in the previous step (5) Optimized data where the expected number of mappings
  • K polyploid organism refers to an organism having K alleles whose chromosomes are K (K is an integer of 2 or more), and includes both homoploids and heteroploids. All organisms that fit this definition are “K polyploid organisms” of the present invention. Even if the organism has a single-phase generation with one chromosome, a multiphase generation with two chromosomes also corresponds to the K polyploid organism in the present invention.
  • the K number of chromosomes of somatic cells of many organisms including humans is 2, that is, a diploid organism.
  • “Lead information” is the base sequence information of the gene fragment of the specimen obtained with a high-performance sequencer.
  • the high-performance sequencer is a sequencer that can be used to implement the present invention and can provide a large amount of read information in a relatively short period of time, and includes a so-called “next-generation sequencer”.
  • Next-generation sequencers at the time of carrying out the present invention include, for example, Genome Sequencer FLX (Roche (454)), Genome Analyzer IIx, HiSeq2000, HiSeq2500, MiSeq (both Illumina), SOLiD (Applied Biosystem), PacBio RS II ( Pacific Biosciences) and the like, but is not limited thereto, and includes all of the high-performance sequencers that are currently provided in the future.
  • two types of lead information are accepted: single-end lead information and pair-end lead information.
  • Single-ended read information is read information about a fixed length or variable length (approximately 50 to 300 bp) at one end of the base sequence of the DNA fragment corresponding to the read.
  • the content of the lead information is steadily changing according to the progress of technology, but in the present invention, the lead information provided now or in the future can be applied.
  • mapping in the present invention has the same meaning as “alignment”.
  • the above step (1) is a step in which the number of expected mappings is digitized.
  • the “expected number of mappings” is the number of mappings defined for each allele for each lead.
  • the step (2) is a step in which the total expected mapping number is calculated based on the expected mapping number.
  • the “total expected number of mappings” is the number of mappings defined for each isoform.
  • Step (3) is a step in which “isoform amount ratio” and “allyl preference”, which are target parameters of the estimation method of the present invention, are calculated.
  • the above steps (4) to (6) are steps related to loop processing in the estimation method of the present invention, and the objective parameters are converged.
  • the optimization method of the present invention when the number of expected mappings in steps (1) and (4) and the isoform amount ratio and allyl preference in steps (3) and (5) are calculated, respectively.
  • an Expectation-Maximization algorithm EM algorithm
  • a variational Bayesian method is used to perform Bayesian estimation.
  • other estimation methods in the estimation method of the present invention include the sparse Bayes method, Gibbsib sampling method, MCMC method, EP method, and PowerEP method.
  • estimation elements in the estimation method of the present invention observation data, objective parameters, and latent variables that link observation data and objective parameters can be cited.
  • the estimation method of the present invention is expressed as follows using these estimation elements.
  • the estimation method according to the present invention includes a step of obtaining the number of expected mappings for isoforms derived from each allele in each lead using the base sequence of the entire read in the allele mixed data as the observation data R, and the isoform amount ⁇ as the target parameter.
  • observation data R n is “the nucleotide sequence of the cDNA read n in the allele mixed data” as described above.
  • the parental allele mixed data is the reference sequence of the paternal cDNA base sequence estimated from the genomic sequence of the specimen and the maternal cDNA base sequence separately from the read of the cDNA by the specimen RNA sequence.
  • the reference sequence can also be obtained by estimation from the read sequence of the genomic DNA of the specimen.
  • the paternal and maternal cDNA sequences used here as references are assumed to be known, but on average, the paternal and maternal cDNA sequences are 99.9% or more of the same sequence, and the individual cDNAs to be mapped It is extremely difficult to identify whether the lead sequence belongs to the paternal or the mother by simply examining the sequence directly.
  • the allele-specific gene expression frequency of the specimen provider can be estimated from such allele mixed data.
  • [1] -1 RNA sequence data of specimen RNA-sequence is performed in a sequencer to obtain cDNA sequence information based on RNA obtained from the specimen.
  • Samples include the above K polyploid organisms, more specifically mammals such as humans, monkeys, dogs, cats, horses, cows, sheep, fish, reptiles, amphibians, insects, seed plants, gymnosperms. It is a specimen capable of extracting total RNA such as mosses and bryophytes.
  • the body fluid such as blood
  • the biological tissues such as organs, fish, reptiles, amphibians, and insects can be used in addition to body fluids and organs.
  • plants including bryophytes
  • they are plant bodies, seeds, spores, etc., and if necessary, appropriate preservation treatment is performed.
  • RNA sequencing As a typical process, total RNA is extracted from the specimen, and poly A-RNA obtained using poly (T) oligo beads is reverse-transcribed to obtain single-stranded cDNA. Usually, this is fragmented and further made into a double strand using a random primer, and then DNA amplification by PCR or the like is performed by a predetermined operation. Subsequently, the PCR amplification product is sequenced to obtain “the nucleotide sequence of the read obtained by RNA sequencing”. As described above, the nucleotide sequence may be derived from a “single-end read” sequenced from only one side of the fragment or from a “pair-end read” sequenced from both sides of the fragment. Applicable to both of these leads.
  • the cDNA read information is electronic information that can be used in a computer, specifically, “FASTA format” or “FASTQ format” in which a sequence sequence and a quality value indicating the reading accuracy of the sequence are described. Provided. This is the “sample RNA sequence data”.
  • [1] -2 cDNA sequence for each allele
  • the cDNA sequence for each allele used as a reference sequence is a sequence obtained by estimation from a K polyploid genome sequence.
  • the cDNA sequence for each allele means a paternal or maternal allele cDNA sequence for humans.
  • the cDNA sequence can be obtained by directly applying an estimation algorithm to the genome sequence when the K polyploid genome sequence of the specimen providing source is obtained. Even when the K polyploid genome sequence cannot be obtained, it can be obtained by executing the step of estimating the K polyploid genome sequence from the read information of normal genomic DNA. It is highly preferred that the source of the reference sequence is the same as the source of the RNA sequence data described above.
  • the specimen for obtaining genomic DNA is not particularly limited as long as it can extract genomic DNA. If it is a mammal including a human, not only body fluids such as blood, biological tissues such as organs, but also hair, nails, etc. Etc. are also possible.
  • eggs can be used in addition to body fluids and organs, and in the case of plants (including bryophytes), they are plants, seeds, spores and the like.
  • Genomic DNA is subjected to gene amplification such as PCR by a predetermined operation, and a DNA fragment of the gene amplification product is obtained.
  • whole genome sequence data is obtained as read information by a technique in which the PCR amplification product is sequenced or not.
  • the nucleotide sequence of the lead may be derived from a “single-ended read” sequenced from only one side of the fragment, or from a “paired-end read” sequenced from both sides of the fragment. Can be applied to both of these leads.
  • the whole genome read information includes electronic information that can be used in a computer, specifically, a sequence sequence and a quality value indicating the reading accuracy of the sequence, such as “FASTA format” or “FASTQ format”, etc. Provided in.
  • This whole genome read information is mapped to a reference sequence, and the mapping information and depth are provided in SAM (sequence alignment / map format) or BAM format that is binary (binary) of this.
  • SAM sequence alignment / map format
  • BAM format binary (binary) of this.
  • VCF Variariant Call Format
  • a desired putative diplotype genome sequence can be obtained on both the paternal and maternal sides by further applying a diplotype genome construction estimation algorithm to this VCF.
  • a desired putative diplotype genome sequence can be obtained on both the paternal and maternal sides by further applying a diplotype genome construction estimation algorithm to this VCF.
  • the allele mixed data is the RNA sequence data of the sample obtained by (1) -1 for each allele obtained from (1) -2 (for humans, paternal and maternal).
  • the mapping information is added by mapping all the alleles at the same time with respect to the cDNA sequence inferred in (Alleles), and is provided in the SAM format or the BAM format.
  • the mapping is preferably performed in a format including all possibilities, that is, performed in a format allowing a plurality of mappings for one lead information.
  • observation data R n is a nucleotide sequence in the n-th read data among N allele mixed data (N is a natural number: converted to the number of reads) as described above. It is assumed that there are N independent uniformly distributed read data. In the case of a single-ended read, one observation data R n is applied to one lead. In the case of a pair-end read, two observation data corresponding to nucleotide sequences at both ends for one fragment. Are applied as a pair.
  • the purpose parameter object parameter is a parameter estimate based on the observed data R n is performed.
  • the present invention involves two objective parameters.
  • T is the number of isoform types
  • Isoforms are transcripts that are transcribed from the same gene but have different mRNA sequences, and the number of isoforms is the total number of isoforms totaled for all genes. Means the base sequence of cDNA.
  • the difference between isoforms may be due to mutations in bases such as SNPs in the genomic DNA sequence and may be due to alternative splicing of mRNA. In the present invention, the latter is the latter.
  • ⁇ tk is a scalar variable representing the allyl preference .
  • the allele-specific gene expression frequency which is the object of the present invention, can be estimated.
  • ⁇ t1 is abbreviated as ⁇ t .
  • the isoform is paternal allele-specific if close to 1 and maternal allele-specific if close to 0. Presumed. Moreover, it can be estimated by estimating the desired parameters theta t above, together the amount of the paternal or maternal isoforms are expressed in reality.
  • Latent variable The latent variable in the present invention describes from which isoform the observed data R n was generated, from which allele, from which location of the allele-specific isoform. Therefore, it is an unobserved variable.
  • at least the above three latent variables [(a) T n related to isoform selection of lead n, (b) H n related to haplotype selection of lead n , and (c) related to the start position of lead n.
  • Latent variable H n as described above, it relates to a haplotype selected lead n, is a variable that depends on the phi tk and T n.
  • H n 0 means that lead n originates from the paternal allele
  • H n 1 means that lead n originates from the maternal allele. Can be set to mean.
  • the latent variable S n is a variable depending on the T n and H n regarding the start position of the lead n.
  • S n s occurs from lead n from position s (1 ⁇ s ⁇ l th ⁇ L + 1) (l th is the length of haplotype h isoform t and L is the lead length) It means that Here, the length l th of the isoform is generally several hundred bases to several thousand bases, and is generally longer than the base length of the lead.
  • a start position s of 1 means that a read has been read from the first base of the isoform.
  • S n In other words, the when mapping the lead reference cDNA, means that the start position in the reference.
  • FIG. 1 is a directed graphical model in a model using single-ended lead information, illustrating dependency relationships between parameters and variables in the present invention described above.
  • the complete likelihood of the data ( ⁇ 1 ,..., ⁇ T ) ′ is defined as a T-dimensional vector, and ⁇ is defined as a T ⁇ K matrix (elements are ⁇ tk ).
  • ⁇ is a T ⁇ 2 matrix
  • ⁇ t1 and ⁇ t2 are parameters for each t.
  • the parameters can be expressed as ⁇ t1 and 1 ⁇ t1 .
  • the phi t1 abbreviated as phi t, represents the phi t and 1-phi t as a parameter of allyl preference.
  • ⁇ ) is the probability that lead n occurs from isoform t, given ⁇ . This probability can be calculated as p (T n t
  • ⁇ ) ⁇ t ((1) (a)).
  • p (H n h
  • the probability, p (H n 0
  • p (S n s
  • Equation (2) (Where subst (,,) is a base quality score-dependent substitution probability function that takes a numerical value obtained by subtracting a sequence substitution error from 1, and r n [x] is a base character at position x of lead n Q n [x] is the base quality score at position x of lead n and c t [x] is the base letter at position x of the DNA sequence of isoform t of haplotype h) Can be calculated.
  • the base quality score substitution probability function, “subst (,,)” can be determined according to the Phred base quality score, and can be estimated from the best alignment of the lead to the reference DNA sequence from the DNA-Seq data.
  • the Phred base quality score is a score that indicates a base reading accuracy provided together with the base sequence information output from the high-performance sequencer in the FASTQ format, that is, a score indicating the error rate output by the sequencer (Phred quality score).
  • the above formula (2) is a calculation formula that takes into account sequence substitution errors, but a calculation formula that takes into account sequence insertion / deletion errors can also be easily derived (Nariai et al., Bioinformatics: 15; 29 (18)).
  • F n is the length of a base fragment (fragment) inferred from the mapping of a pair of paired end reads “R na and R nb ” to a reference sequence.
  • ⁇ ) is the probability that lead n occurs from isoform t, given ⁇ . This probability can be calculated as p (T n t
  • ⁇ ) ⁇ t ((3) (a)).
  • p (F n f
  • l th is the length of the reference sequence of haplotype h of isoform t.
  • the length distribution d F (x) of the base fragment is given as a normal distribution with an average ⁇ F and a standard deviation ⁇ F , for example.
  • the mean ⁇ F and standard deviation ⁇ F may be specified if the base fragment length distribution is experimentally known when the base fragment is created. May be estimated ((3) (b)).
  • p (H n h
  • subst (,,) is a base quality score-dependent substitution probability function that takes a value obtained by subtracting a sequence substitution error from 1
  • r na [x] is the first base sequence of the set of leads n
  • R na is the base character at position x
  • q na [x] is the base quality score of the first base sequence position x in the set of leads n
  • c tha [x] is the haplotype of isoform t
  • r nb [x] is the second base sequence R nb of the set of lead n
  • the base character at position x, q nb [x] is the base quality score of the first base sequence position x of the set of leads n
  • c thb [x] is the DNA of haplotype h of isoform t
  • the above formulas (3-1) and (3-2) are the basis for obtaining the posterior probabilities and posterior distributions of the latent variables according to the estimation method of the present invention using the paired end model.
  • the above formula (2), formulas (3-1), and (3-2) are calculation formulas that take into account sequence substitution errors, but formulas that take sequence insertion / deletion errors into account are also available. Similarly, it can be easily derived (Nariai et al., Bioinformatics: 15; 29 (18)).
  • the present invention provides a computer system for realizing the estimation method of the present invention, and computer software including an algorithm capable of realizing the estimation method of the present invention in a computer.
  • the specific method for performing the estimation method of the present invention is not limited, and examples thereof include a maximum likelihood estimation method and a Bayes estimation method.
  • the parameter estimation algorithm based on the maximum likelihood estimation method includes the EM method
  • the parameter estimation algorithm based on the Bayesian estimation method includes the sparse Bayes estimation method, the variational Bayes estimation method, the Gibbs sampling method, the MCMC. Method, EP method, PowerEP method and the like.
  • the EM method or the variational Bayes estimation method is preferable.
  • the target allele-specific gene expression can be accurately obtained even by using the EM method or the variational Bayes method, which does not impose a relatively heavy computational burden. It is possible to estimate.
  • EM expectation-maximization
  • the EM method is a method for obtaining a maximum likelihood estimate in a probability model in which a latent variable exists.
  • the following contents can be used.
  • the estimation method of the present invention by the EM method is exemplified as the estimation method of the present invention characterized by performing the following steps (a) to (f) (this method is also referred to as “EM method of the present invention”). ).
  • Steps (a) and (b) of the EM method of the present invention are steps for defining initial values of the method for performing repeated optimization).
  • the initial value of Z nths is defined as 1 to N for n, 1 to T for t, and 0 or 1 for diploid.
  • h is 1 to K.
  • l is the length of isoform t of haplotype h
  • L is the lead length.
  • ⁇ ( ⁇ 1 ,..., ⁇ T ) ′ is defined as a T-dimensional vector, and ⁇ is defined as a T ⁇ K matrix (elements are ⁇ tk ).
  • ⁇ t1 and ⁇ t2 are parameters for each t.
  • the parameters can be expressed as ⁇ t1 and 1 ⁇ t1 I can do it.
  • the phi t1 abbreviated as phi t, represents the phi t and 1-phi t as a parameter of allyl preference.
  • Step (e) and step (f) are as described above.
  • a convergence criterion a relative change of 10 ⁇ 3 for an isoform with an abundance parameter ⁇ t > 10 ⁇ 7 is used as the convergence criterion, although a different convergence criterion can be used.
  • step (e) is performed one or more times, but this repeating step is not performed once or sufficiently. It is also possible to finish the step in a state where it has not converged and estimate the allele-specific gene expression.
  • variational Bayes method is a method for estimating more accurately by estimating the posterior probability distribution of parameters in the Bayesian estimation method.
  • the present method is performed with the following contents.
  • ⁇ ( ⁇ ) is a gamma function) Is used.
  • a single hyperparameter ⁇ 0 is used for all isoforms based on the assumption that there is no prior knowledge of the relative differences in isoform amounts.
  • the single hyperparameter ⁇ 0 controls the complexity of the target parameter.
  • ⁇ 0 ⁇ 1 ⁇ 0 ⁇ 1 can be interpreted as a prior count of leads assigned to the isoform, and when ⁇ 0 ⁇ 1, the prior distribution tends to have some isoform amount zero. give. Then, ⁇ 0 that maximizes the lower limit of the log marginal likelihood is selected.
  • ⁇ ( ⁇ 1 ,..., ⁇ T ) ′ is defined as a T-dimensional vector, and ⁇ is defined as a T ⁇ K matrix (elements are ⁇ tk ).
  • ⁇ t1 and ⁇ t2 are parameters for each t.
  • the parameters can be expressed as ⁇ t1 and 1 ⁇ t1 .
  • the phi t1 abbreviated as phi t, represents the phi t and 1-phi t as a parameter of allyl preference.
  • ⁇ distribution (Where ⁇ t , 1 > 0 and ⁇ t , 2 > 0 are hyperparameters, and B (•, •) is a ⁇ function.) Is used.
  • ⁇ t , 1 and ⁇ t , 2 can be said to be pre-counts of leads assigned to the paternal and maternal alleles, respectively, in order to calculate the paternal / maternal ratio.
  • the prior information of these data is known Can also be used as the initial values of ⁇ t and ⁇ tk .
  • ⁇ 0 and ⁇ t , 1 and ⁇ t , 2 can be set according to the known information.
  • the lower limit of the logarithmic marginal likelihood is maximized iteratively by a variational Bayesian estimation algorithm.
  • the estimation method of the present invention by the variational Bayes method is exemplified as the estimation method of the present invention characterized by performing the following steps (a) to (e) (this method is referred to as “variation Bayes of the present invention. Law ").
  • step (E) A step of evaluating the convergence of the expected value of the posterior distribution of ⁇ t and ⁇ t obtained in step (d), and if convergence at the expected value is recognized, the expected convergence value is expressed as ⁇ t And ⁇ tk are estimated values. If convergence is not recognized, the repetition of steps (b), (c), and (d) is determined.
  • Steps (a) and (b) of the variational Bayes method of the present invention are steps for defining initial values of the method for performing repeated optimization. Specifically, for each isoform t, assuming that the expression level of each isoform is equal for ⁇ * t , which is a parameter of q * ( ⁇ ), Is preferably set.
  • Step (c) is a step of calculating E z [Z nths ] given the current estimates of q * ( ⁇ ) and q * ( ⁇ ).
  • the value of Z nths is each number from 1 to N for n, each number from 1 to T for t, 0 or 1 for h, and 1 ⁇ s ⁇ l for s. th ⁇ L + 1, l th is the length of isoform t of haplotype h, and L is the lead length.
  • E z [Z nths ] is based on the current estimate of q * ( ⁇ ) and q * ( ⁇ ), Is calculated as
  • Step (d) is a step of calculating E ⁇ [ ⁇ t ] and E ⁇ [ ⁇ t ] given the current estimate of q * (Z).
  • q * ( ⁇ ) is also a Dirichlet distribution
  • p ( ⁇ ) is a conjugate prior distribution
  • q * ( ⁇ t ) is also a ⁇ distribution
  • prior distribution p ( ⁇ t ) is a conjugate prior distribution
  • step (d) as a preferred example of a convergence criterion, a relative change 10 ⁇ 3 for an isoform with an abundance parameter E ⁇ [ ⁇ t ]> 10 ⁇ 7 is used as the convergence criterion, but a different convergence. It is also possible to use criteria.
  • step (d) is performed one or more times, but this repeating step is not performed once or sufficiently. It is also possible to estimate the allele-specific gene expression by completing the step without convergence.
  • the log marginal likelihood in the variational Bayes method of the present invention is Can be broken down as
  • L (q) the log marginal likelihood is L (q) Gives a lower bound. That is, since the Cullback-Liber distance is always 0 or more, L (q) constitutes the lower limit of the marginal likelihood. It is generally indicated that L (q) (the lower limit of the logarithmic marginal likelihood) is increased every time the above steps (b) and (c) are repeatedly updated (Bishop CM. Pattern Recognition and Machine Learning. Springer Science: Business Media, LLC, New York, NY, USA; 2006).
  • Computer System and Computer Program of the Present Invention [1] Computer System of the Present Invention
  • the computer system of the present invention is a system that serves as means for performing the above-described estimation method of the present invention, and the program of the present invention performs the estimation method of the present invention on the computer system of the present invention.
  • a computer program having an algorithm for making “Algorithm” means a formalized representation of a procedure for solving a problem, as in the general concept of the computer field.
  • the computer system of the present invention can include hardware related to a normal computer system.
  • a CPU Central Processing Unit
  • DSP Digital Signal Processor
  • RISK Reduced Instruction Set Computer
  • a RAM Random Access Memory
  • main memory main memory
  • keyboard a mouse
  • touch panel or the like
  • operation unit a display
  • Corresponding serial or parallel interface, etc. corresponding “output / input interface (IF) unit”, video memory and D / A converter unit, and output analog signal according to video system of display unit “communication interface (IF)” Part ".
  • the communication IF unit can exchange data with external information, particularly human genome information such as a human genome database.
  • the “arithmetic processing unit” is operated by operating the “operation unit” to acquire the data of the human genome database in particular through the “communication IF unit”, record it in the “recording unit”, and appropriately store the data from the “recording unit”. After reading into the “temporary storage unit” and performing a predetermined calculation process, the result is recorded in the “recording unit” again.
  • the “arithmetic processing unit” creates screen data for prompting operation of the “operation unit” and screen data for displaying the processing results, and displays these images on the “display unit” via the video RAM of the input IF unit. To do.
  • the program of the present invention is used or recorded in advance in the “recording unit”, or recorded in an external hardware resource, and in the “arithmetic processing unit” according to the algorithm described in the “arithmetic processing unit” as necessary. Is done.
  • a computer system for optimizing the expected number of mappings for all allele-derived isoforms in all leads for data containing mixed cDNA read information comprising a recording unit and an arithmetic processing unit; (1) In the recording unit, read information of cDNAs derived from all alleles of the test organism is recorded as observation data of lead sequences and lead mapping destinations, (2) The arithmetic processing unit alternatively executes the following initialization step (2) -1 or (2) -2 based on the read observation data, (2) -1: Calculation process of initial value of variable ⁇ t related to isoform amount, calculation process of initial value of variable ⁇ tk related to allyl preference , (2) -2: Three latent variables that mediate the base sequence of the above-mentioned variables ⁇ t and ⁇
  • the calculation processing unit performs a calculation process of the index variable Z nths , (4)
  • the first likelihood estimation values of the variables ⁇ t and ⁇ tk are calculated based on the index variable Z nths calculated in step (2) -2 or step (3).
  • Update value is calculated, (5)
  • said step (3) and step (4) are again performed based on the 1st update value of the maximum likelihood estimated value of the variables (theta) t and (phi) tk calculated at said step (4).
  • the loop process for calculating the second update values of the variables ⁇ t and ⁇ tk is repeatedly executed until the difference between the new update value and the previous update value is substantially not recognized. Then, recording is performed in the recording section as data of ⁇ t and ⁇ tk in which the converged variables ⁇ t and ⁇ tk are optimized; A computer system for optimizing read information of cDNAs derived from individual alleles.
  • a computer system for optimizing the expected number of mappings for all allele-derived isoforms in all leads for data containing mixed cDNA read information comprising a recording unit and an arithmetic processing unit; (1) In the said recording part, the read information of cDNA derived from each allele of the test organism is recorded as the observation data of the lead sequence and lead mapping destination, (2) The arithmetic processing unit alternatively executes the following initialization step (2) -1 or (2) -2 based on the read observation data, (2) -1: Calculation process of updated value of posterior distribution of ⁇ t based on initial value of hyper parameter ⁇ t indicating distribution of prior knowledge about abundance of isoform t, and lead assigned to each allele Posterior of ⁇ tk based on initial values of hyper parameters ⁇ t
  • Computer program of the present invention is a computer program having an algorithm for executing the above-described computer system of the present invention. Unless otherwise specified, the symbols and concepts are as described above. is there.
  • the third function; (4) In the arithmetic processing unit, based on the posterior probability of the indicator variable Z nths calculated in step (2) -2 or step (3), the maximum likelihood estimated values of the variables ⁇ t and ⁇ tk A fourth function for performing a calculation process of the first update value; (5) In the arithmetic processing unit, the steps (3) and (4) are performed again based on the first update value of the maximum likelihood estimated value of the variables ⁇ t and ⁇ tk calculated in the step (4). And repeatedly execute the loop process in which the second update values of the variables ⁇ t and ⁇ tk are calculated until the difference between the new update value and the previous update value is substantially not recognized.
  • the calculation processing unit Based on the distributions of the variables ⁇ t and ⁇ tk calculated in step (2) -1, the calculation processing unit performs a calculation process of the distribution of the index variable Z nths (the expected mapping number of lead n).
  • step (3) and step (4) are again performed based on the 1st update posterior distribution of the variable (theta) t and (phi) tk calculated at said step (4), ,
  • the loop process in which the second updated posterior distribution of the variables ⁇ t and ⁇ tk is calculated, the difference between the newly updated expected posterior distribution and the previously updated expected posterior distribution is A fifth function that repeatedly executes until it is not recognized, and records the converged expected values of ⁇ t and ⁇ tk as optimized ⁇ t and ⁇ tk data in the recording unit;
  • the computer program which optimizes the read information of cDNA derived from each allele characterized by including the algorithm which implement
  • the computer program of the present invention can be described in, for example, C language, Java (registered trademark), Perl, Python, or the like.
  • the present invention further provides a computer-readable recording medium or a recording medium connectable to a computer (hereinafter also referred to as a recording medium of the present invention), in which the program of the present invention is recorded.
  • recording media include magnetic media such as flexible disks, flash memories, and hard disks, optical media such as CD, DVD, and BD, magneto-optical media such as MO and MD, and the like. is not.
  • a typical computer system of the present invention is characterized by executing the program of the present invention.
  • FIG. 2 is a flow sheet showing a processing flow of the estimation method of the present invention for estimating the allele-specific gene expression frequency of a human specimen.
  • Each element of the processing flow is electronically performed in a computer based on instructions for implementing a computer software algorithm via a computer network or a database as necessary.
  • Box B1-1 indicates the presence of the read data of the genomic DNA of the specimen
  • Box B1-2 indicates the presence of the human genome reference sequence.
  • Read data (box B1-1) is a part or all of the base sequence of each genomic DNA fragment (read) obtained by processing the genomic DNA of a specimen with a high-performance sequencer. Or a pair-end format. Usually, information indicating the reading accuracy (quality score) is also added (FASTQ: a term derived from the FASTA format representing the base sequence of DNA).
  • FASTQ a term derived from the FASTA format representing the base sequence of DNA.
  • the “human genome reference sequence” (box B1-2) results in a single individual genome sequence information, and the information source is not particularly limited. It is also assumed that the race is different in the specimen and the genomic sequence information, and it contains heterogeneity in the isoform gene sequence.
  • human genome reference sequences include human genome reference sequences (hg19 and the like) managed by UCSC (University California Santa Cruz) and human genome reference sequences (GRCh37 and the like) managed by Genome Reference Consortium. It is not limited to.
  • Step S1 is a step of mapping genomic DNA, and the step of mapping the “read data” in box B1-1 to the “human genome reference sequence” in box B1-2.
  • this mapping it is possible for those skilled in the art to create a computer program including an algorithm capable of realizing the mapping by a computer, but it is also possible to use already provided software. Examples of existing mapping software include BWA-MEM (http: // bio-bwa. Source. Net /).
  • the above-mentioned human genome reference sequence preferably includes a “decoy sequence”.
  • “Decoy sequence” is a genomic sequence that does not exist in a reference sequence prepared in advance. Without a decoy sequence, reads that are not derived from sequences that are not registered in the human genome reference sequence are more likely to be mapped anywhere in the existing human genome reference sequence, which may reduce the accuracy of the mapping. There is.
  • hs37d5 ftp://ftp.1000genmes.ebi.ac.uk/voll/ftp/technical/phase2_reference_assembly_sequence/hs37d5.fa.gz
  • Box B2 is a box indicating the presence of the “genomic DNA mapping result” in S1, and the mapping result exists in a format such as SAM format or BAM format.
  • Step S2 is a step of performing mutation calling and haplotype fading on the “result of genomic DNA mapping” in box B2.
  • Mutation calls are means for specifying mutations such as SNPs
  • haplotype fading is means for estimating the base sequence of genomic DNA containing mutations called at the haplotype level. It is possible for a person skilled in the art to create a computer program including an algorithm that can implement mutation calls and fading separately or together on a computer, but it is also possible to use already provided software . Examples of existing mutation calls and fading software include HapMonster (http://nagasakilab.csml.org/en/hapmonster).
  • Box B3 indicates that the sequence information of the genomic DNA estimated as a haplotype obtained as a result of the mutation call and fading performed in step S2 exists in the VCF format.
  • Step S3 is a step of estimating and constructing the paternal haplotype genomic DNA sequence and the maternal haplotype genomic DNA base sequence by estimating the diplotype sequence information from the haplotype genomic DNA sequence information in box B3.
  • a person skilled in the art to create a computer program including an algorithm capable of realizing this estimation construction by a computer, it is also possible to use already provided software. Examples of the existing software for estimation construction include vcf2diploid (Rozowskyowet al., Mol Syst Biol. 2011 Aug 2; 7: 522.).
  • Box B4 (B4-1 and B4-2) indicates the presence of the haploid genome base sequences of the paternal and the maternal obtained in step S3.
  • Box B4-1 is the paternal genomic DNA base sequence
  • box B4-2 is the maternal genomic DNA sequence.
  • the diplotype genomic base sequence is directly used and the previous steps are used. The next step S4 can be performed from the stage of the box B4 without going through the above.
  • Step S4 is a step of estimating the paternal and maternal genomic DNA base sequences of Box B4 into cDNA base sequences.
  • the estimation in step S4 can be performed by creating a computer program including an algorithm capable of realizing this, but it is also possible to use already provided software.
  • rSeq Jiang H, Wong WH. Statistical inferences for isoform expression in RNA-Seq. Bioinformatics.2009; 25: 1026-1032
  • the existing software for estimation can be used as the existing software for estimation.
  • Box B5 (B5-1 and B5-2) indicates the presence of paternal and maternal cDNA base sequences obtained in step S4.
  • Box B5-1 is the paternal cDNA sequence
  • box B5-2 is the maternal cDNA sequence.
  • Step S5 is a step in which “sample RNA sequence data” (read data) shown in box B6 is mapped using a combination of the base sequences of paternal and maternal cDNAs shown in box B5 as a reference sequence. .
  • the read data may be single end data or pair end data.
  • Box B7 indicates the presence of electronic data based on the mapping obtained in step S5, that is, “parental allele mixed data” used as observation data R n .
  • step S6 is a step of optimizing the “parental allele mixed data”, and details thereof are shown in FIGS. 3 and 4 as a flow sheet partially overlapping with the above.
  • the mode shown in FIG. 3 is a step of performing maximum likelihood estimation by the EM algorithm in order to optimize the “parental allele mixed data” and estimate the allele-specific gene expression frequency.
  • the mode shown in FIG. 4 is a step of performing Bayesian estimation by the variational Bayes method in order to optimize the “parental allele mixed data” and estimate the allele-specific gene expression frequency.
  • Step S6-11 shows the contents of step S5 and box B7 shown in FIG. 2, and indicates the existence of “parental allele mixed data”.
  • N is the “total number of leads”
  • T is the “total number of isoform types”.
  • E [Z nths ] is initially set as 1 / M n (M n is the total number of alignment destinations of the observed lead n) on the premise of no information and uniform distribution in advance.
  • step S6-13 based on ⁇ t and ⁇ t that were updated immediately before (initially the initial values in step S6-12 above), the amount of allocation to position s of allyl h of isoform t of lead n
  • step S6-13 based on ⁇ t and ⁇ t that were updated immediately before (initially the initial values in step S6-12 above), the amount of allocation to position s of allyl h of isoform t of lead n
  • step S6-14 the M step of the EM method of the present invention, which calculates the maximum likelihood estimates of ⁇ t and ⁇ t based on E z [Z nths ] (r nths ) calculated in step S6-13.
  • This is a step corresponding to.
  • ⁇ t and ⁇ t that maximize the log likelihood are calculated.
  • Isoform t takes 1 to T (total number of isoform types), It is.
  • Step S6-15 is a step for selecting whether to perform loop processing that repeats steps S6-13 and S6-14 or to end the loop processing.
  • steps S6-13 and S6-14 executed in the previous loop.
  • the change between the measured “ ⁇ t ′ ” and the “ ⁇ t ′ ” is substantially not recognized, not only the “ ⁇ t ” but also “ ⁇ t ” is recognized as the “convergence value”, and a series of optimization processes are performed. finish. The loop process is executed until this convergence is achieved.
  • the convergence check performed here may be performed on “ ⁇ t ” or on both “ ⁇ t ” and “ ⁇ t ”, but the same result is obtained in either case.
  • Step S6-22 is a parameter E [Z nths ] (assignment to the position s of the allele h of the isoform t of the lead n) used in the optimization process by the execution of the computer program.
  • Z nths 1 expected value per minute
  • ⁇ t ratio of leads assigned to isoform t relative to total lead
  • ⁇ t ratio of expression level from father allele in isoform t
  • S is “start position of lead n”
  • l th is the length of h allyl of isoform t
  • L is the lead length.
  • E [Z nths ] is initially set as 1 / M n (M n is the total number of alignment destinations of the observed lead n) on the premise of no information and uniform distribution in advance.
  • E ⁇ [ ⁇ t ], which is an expected value of ⁇ t is given as ⁇ * t / ⁇ t ′ ⁇ * t ′ with the Dirichlet distribution as a prior distribution.
  • ⁇ * t ⁇ 0 + N / T.
  • E ⁇ [ ⁇ t ], which is an expected value of ⁇ t is given as ⁇ * t, 1 / ⁇ * t, 1 + ⁇ * t, 2 with a ⁇ distribution as a prior distribution.
  • ⁇ * t, 1 ⁇ t , 1 + N / 2T
  • ⁇ * t, 2 ⁇ t , 2 + N / 2T.
  • ⁇ t, 1 and ⁇ t, 2 are hyperparameters that give uniform distribution, assuming both father allele and mother allele are equally likely, and assuming no prior knowledge.
  • the estimated value of E z [Z nths ] can be calculated in the following manner.
  • Step S6-24 is a step of calculating E ⁇ [ ⁇ t ] and E ⁇ [ ⁇ t ] based on the estimated value of E z [Z nths ] calculated in step S6-23. This is a step corresponding to the M step of the Bayes method.
  • E ⁇ [ ⁇ t ] is based on the estimated value of E z [Z nths ] calculated in the immediately preceding S6-23. Is calculated as
  • E ⁇ [ ⁇ t ] is based on the estimated value of E z [Z nths ] calculated in immediately preceding S6-23. Is calculated as
  • Step S6-25 is a step for selecting whether to perform a loop process that repeats steps S6-23 and S6-24 or to end the loop process.
  • steps S6-23 and S6-24 executed in the previous loop in the estimated lead value E ⁇ [ ⁇ t ] assigned to the isoform t with respect to the total amount of leads calculated in step S6-24.
  • E ⁇ [ ⁇ t ] assigned to the isoform t with respect to the total amount of leads calculated in step S6-24.
  • the optimization process is executed by [1] maximum likelihood estimation by the EM algorithm and [2] Bayes estimation by variational Bayes.
  • ASE is a term that means allele-specific gene expression.
  • NA12878 diploid genomic sequences [these were constructed from hg19 and are available on the website (http://sv.gersteinlab.org / NA12878_diploid) based on synthetic RNA-Seq data (30 million reads, 100 bp ⁇ 2 (paired end), average fragment size 400 bp, standard deviation 40 bp) from the sample RNA sequence. Prepared as the obtained lead information.
  • Paternal and maternal cDNA sequences were generated from diploid genomic sequences based on the UCSC genome annotation file (refFlat.txt) with rSeq (version 0.2.1), then 10,000 isoforms were randomly selected Expression levels were assigned to follow a lognormal distribution. Subsequently, two sets of RNA-Seq data with 0.1% substitution, deletion, and insertion errors were generated based on the above RNA-Seq data, while one was generated based on the null hypothesis that there was no ASE. However, the other was created based on the alternative hypothesis that ASE exists in some of the isoforms.
  • the lead information is mapped to the haplotypes of both parents, and the “ ⁇ k” option is specified for Bowtie2, and all possible alignments are retained.
  • the BAM file was used as an input, and the lead alignment for the isoform was optimized between the paternal and maternal alleles by the modified Bayesian estimation algorithm shown in the flow sheet of FIG.
  • FIG. 5 shows the true distribution of paternal / maternal ratios for the null hypothesis and the ASE hypothesis, the distribution predicted based on the alignment of the software based on the above existing approach of the lead to the diploid genome, And the distribution estimated with the estimation method of this invention is shown. Isoforms with more than 10 assigned leads were compared. Regardless of the presence or absence of ASE, the predicted distribution with the estimation method of the present invention was more similar to the true distribution, especially in the region where the paternal / maternal ratio was close to 0 or 1. On the other hand, the predicted distributions with the above existing approaches show “peaks” in their extreme regions, but in fact they did not exist in the true distribution.
  • the preferred result of the estimation method of the present invention is that in Bayesian estimation, when the ⁇ distribution for haplotype selection is updated, a pre-count of 1 is added to each allele of the isoform to calculate the paternal / maternal ratio.
  • a smoothing effect referred to as Laplace smoothing or add-one smoothing. This feature may be particularly beneficial for isoforms with originally low expression levels or when there are not many heterozygous SNPs that can be used to distinguish alleles.
  • RNA-Seq data 100 bp ⁇ 2 generated from the lymphoblastoid cell line GM12878 (Marinov et al., Genome Research, 24, 496-510 (2014)) was used for the estimation method of the present invention.
  • the estimation method of the present invention can accurately estimate allele-specific gene expression based on read information obtained by RNA sequencing of a specimen.
  • This technology is expected to enable the creation of biomarkers for early detection of diseases such as cancer, drug discovery, and association analysis between genome mutations and disease-related gene functions.
  • diseases such as cancer, drug discovery, and association analysis between genome mutations and disease-related gene functions.
  • not only diploids but also triploids more than triploids are known, especially in plants and fish, and artificial polyploids are produced as a means for improving various varieties. Estimating the frequency of allele-specific gene expression in these polyploids, and obtaining clues to improve the yield to higher yields is considered to be extremely important for industrial policy.
  • Step B1-1 Box indicating presence of read data of genomic DNA B1-2 ... Box indicating presence of human genome reference sequence S1 ... Mapping of genomic DNA Step B2 ... Result of mapping of S1
  • the data of step S3... B3 indicating the presence of the mutation call and fading result of step B3... S2 for performing the mutated call and haplotype fading on the mapping result of step S2.
  • Step 4-1 for constructing a diploid genome sequence based on the above: Box B4-2 indicating presence of paternal haploid genome sequence B4-2: Box S4 indicating presence of maternal haploid genome sequence Estimating paternal and maternal haploid genome sequences into cDNA sequences Step B5-1 ...
  • Paternal cD Box B5-2 indicating the presence of the A sequence Box B6 indicating the presence of the maternal cDNA sequence Box B5 indicating the presence of the RNA sequence read data
  • Mapping of the parental cDNA sequence of the RNA sequence data is performed Step B7 ... Box S6 indicating the existence of the mapping result of S5 ... Step S6-11 for optimizing the parental allele mixed data Step S6-21 ... indicating the presence of the parental allele mixed data
  • Step S6-12 indicating the presence of allele mixed data
  • Step S6-22 for initializing parameters for performing optimization processing
  • Step S6-22 for initializing parameters for performing optimization processing 13...
  • Step S6-23 for executing the function of the E step for performing the optimization process.
  • Step S6-14 for executing the function of the E step for performing the optimization process
  • Step S6-24 for executing the function of the M step for performing the optimization process.
  • Step S6-15 for executing the function of M step for performing ... Loop for performing the optimization process
  • Step S6-25 for performing the convergence function ... Loop for performing the optimization process
  • Step B8 for executing the convergence function:

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Genetics & Genomics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本発明は、諸生物におけるアリル特異的遺伝子発現を的確に推定可能な統計確率的な手段を提供することを課題とし、本発明者らは、倍数体生物のｃＤＮＡのリード情報が混在したデータに対して、（１）～（４）の内容をコンピュータにおいて行うことにより、上記の課題を解決し得ることを見出した。 （１）各リードにおける各アイソフォーム及び当該アイソフォームの各アリルに対する期待マッピング数の数値化。 （２）（１）において数値化された期待マッピング数の各アイソフォーム毎における合算による合計期待マッピング数の算出と、当該合計期待マッピング数の各アリル分の算出。（３）（２）において算出された合計期待マッピング数に関する、（ｉ）アイソフォーム量比の算出と（ii）アリル好選度の算出。 （４）上記ステップのループ処理による収束値をアリル特異的遺伝子発現の最適データとしての認定。

Description

アリル特異的遺伝子発現頻度の推定方法、推定用コンピュータシステム及び推定用プログラム 関連出願
　本出願は、２０１５年１２月２８日に出願された日本国特許出願２０１５－２５６８６０号に基づく優先権を主張し、その開示の全てが参照によりここに取り込まれる。
　本発明は、核酸情報に基づく遺伝的解析の分野に関する発明であり、より詳しくは、高性能シークエンサから導出されるＲＮＡシークエンス（ＲＮＡ－ｓｅｑ）解析情報に基づき、アリル特異的遺伝子発現頻度、言い換えれば、倍数体生物の各アリルにおける遺伝子発現量を推定する手段に関する発明である。倍数体生物がヒトを含んだ哺乳類等の通常の２倍体生物の場合は、各アリルとは父方由来のアリルと母方由来のアリルである。
　アリル特異的遺伝子発現（Allele-specific expression（ＡＳＥ））は、倍数体生物の各アリルの遺伝子の多様な発現態様に対応するものである。
　例えば、父方由来の遺伝子か、又は、母方由来の遺伝子か、に依存する現象である「ゲノム刷り込み(genomic imprinting)」において、既に研究がなされており、女性のＸ染色体２本のうち１本が不活性化する性染色体の不活性化も、ＡＳＥに基づくものとして知られている(非特許文献１)。最近の研究では、ＡＳＥは比較的普遍的な現象であり(非特許文献２)、多くのシス作用性(cis-acting)配列変異により、遺伝子発現が変化することが明らかになっている(非特許文献３)。ある場合には、３つのアリルにおける遺伝子発現の違いが、結腸直腸癌等の疾患傾向に関連しており、重要なことに遺伝子の転写量は、糖尿病や肥満等の通常の疾患に対する感受性遺伝子座を発見するための量的な指標として用いられ得ることが報告されている(非特許文献４，５)。
　近年、このような背景から、父方由来遺伝子と母方由来遺伝子の各々からの発現、すなわち、アリル特異的遺伝子発現を検討することにより、ａ）癌等の疾患の早期発見のためのバイオマーカーの作出、ｂ）創薬、ｃ）ゲノム変異と疾患関連遺伝子機能との関連解析、が可能となると期待されている。
　さらに、特に植物及び魚類では、２倍体のみならず、３倍体以上の倍数体が知られている他、各種の品種改良の手段として人為的な倍数体の作出が行われている。これらの倍数体におけるアリル特異的遺伝子発現頻度を推定し、より収量の期待出来る品種へ改良する手がかりを得る事等が、産業政策上極めて重要と考えられる。
Lyon MF. Gene action in the X-chromosome of the mouse (Mus musculus L.). Nature. 1961;190:372-373 Knight JC. Allele-specific gene expression uncovered. Trends Genet. 2004;20:113-116 Buckland PR. Allele-specific gene expression differences in humans. Hum Mol Genet. 2004; 13; 2:R255-260. Schadt et al., Genetics of gene expression surveyed in maize, mouse and man. Nature. 2003; 422:297-302. Fairfax et al., Innate immune activity conditions the effect of regulatory variants upon monocyte gene expression. Science. 2014; 343:6175.
　現在、アリル特異的遺伝子発現の検討手法として、ＲＮＡシークエンス（ＲＮＡ－ｓｅｑ）解析が用いられている。この手法は、高性能シークエンサを用いて断片化された検体ＲＮＡに由来するｃＤＮＡのリード情報を用いる方法であるが、現在いくつかの問題点が報告されている。
　第１に、多くの場合、ＲＮＡ－ｓｅｑリードは参照塩基配列の複数位置にマッピングされてしまい、遺伝子発現の定量性に問題が生じている(Mortazavi et al., Nature Methods,5,621-628(2008))。
　第２に、検体が、参照塩基配列とは異なるヘテロのＳＮＰ（一塩基置換）を有する場合には、マッピングの際に偏向(bias)の存在を許してしまう(Degner et al., Bioinformatics,25,3207-3212(2009))。
　特に、この第２の問題点を解決するためのアリル特異的遺伝子発現の推定方法として、例えば、Rozowskyらが２０１１年に行った方法が挙げられる（Rozowsky et al., AlleleSeq:analysis of allele-specific expression and binding in a network. Molecular Systems Biology,7,522,2011）。
　この方法は、ａ）２倍体生物における父方由来と母方由来のアリル（対立遺伝子）のそれぞれのゲノム配列を決定し、ｂ）当該ゲノム配列からスプライシングによりｍＲＮＡとして発現する配列部分（エキソン）のつながりを反映するｃＤＮＡ配列を予測し、ｃ）当該予測ｃＤＮＡ配列に対して、検体ＲＮＡのｃＤＮＡのリード情報のマッピングを行い、当該ｃＤＮＡ配列中のアリル特異的なＳＮＰ、構造変異等を計測し、ｄ）当該マッピング情報から、父方由来と母方由来のアリル各々の発現割合を推定する、という方法である。
　しかしながら、この方法では、ｃＤＮＡ配列中のＳＮＰ、構造変異等が認められない遺伝子座にマッピングしているリードの存在を無視することになり、高精度なアリル特異的遺伝子発現の推定方法とはいえない。
　本発明者らは上記の課題の解決に向けて、検体のＲＮＡ－ｓｅｑのリードのヌクレオチド配列を基に、倍数体生物のアリル特異的遺伝子発現を遺伝子全体として推定するための手段についての検討を行った。
　その結果、当該ＲＮＡ－ｓｅｑのリードのヌクレオチド配列を観測データとして、特定の潜在要素をアリル特異的遺伝子発現頻度に関連付けて仮想して、当該頻度をコンピュータにおいて統計的に推定算出することにより、倍数体生物のアリル特異的遺伝子発現頻度を非常に正確に推定することが可能であることを見出して、本発明を完成した。
［本発明の推定方法］
　本発明は、Ｋ倍数体生物（Ｋは２以上の整数）のｃＤＮＡのリード情報が混在したデータ（以下、「アリル混在データ」ともいう）に対して、以下（１）～（６）のステップが実行される、コンピュータによるアリル特異的遺伝子発現頻度の推定方法（以下、「本発明の推定方法」ともいう）を提供する。
（１）　各リードにおける各アイソフォーム及び当該アイソフォームの各アリルに対する期待マッピング数の数値化が行われるステップ、
（２）　ステップ（１）において数値化された期待マッピング数が各アイソフォーム毎に合算されて合計期待マッピング数が算出され、かつ、当該合計期待マッピング数の各アリル分がそれぞれ算出されるステップ、
（３）　ステップ（２）において算出された合計期待マッピング数は、（ｉ）各アイソフォームに対する合計期待マッピング数の和で除されてアイソフォーム量比が算出され、かつ、(ii)各アイソフォームの各アリルにおいて、当該アイソフォームの合計期待マッピング数の当該アリル分が当該アイソフォームの合計期待マッピング数で除されてアリル好選度が算出されるステップ、
（４）　ステップ（３）において算出された、新たなアイソフォーム量比とアリル好選度に対して再びステップ（１）が実行されて、新たな期待マッピング数が算出されるステップ、
（５）　ステップ（４）において得られた新たな期待マッピング数に対して、再びステップ（２）又は（３）が実行され、アイソフォーム量比とアリル好選度が新たに算出されるステップ、
（６）　ステップ（４）と（５）が、（ｉ）ステップ（４）において算出される期待マッピング数と、前回のステップ（４）で算出された期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、(ii)ステップ（５）において算出されるアイソフォーム量比及び／又はアリル好選度と、前回のステップ（５）において算出されたアイソフォーム量比及び／又はアリル好選度の間における差が全てのアリルにおいて認められなくなるまで、繰り返し実行され、収束した期待マッピング数、あるいは、収束したアイソフォーム量比及び／又はアリル好選度が、最適化されたデータとして認定されるステップ。
　本発明における「Ｋ倍数体生物」とは、染色体がＫ本（Ｋは２以上の整数）のアリルを有する生物のことで、同質倍数体と異質倍数体の双方が含まれる。この定義に当て嵌まる全ての生物が、本発明の「Ｋ倍数体生物」である。染色体数が１本の単相世代を有する生物であっても、染色体数が２本となる複相世代も本発明におけるＫ倍数体生物に該当する。ヒトを含む多くの生物の体細胞の染色体のＫ数は２、すなわち、２倍体生物である。このような通常の２倍体生物の場合は、「雄若しくは父アリルと、雌若しくは母アリル」の概念で、本発明における「アリル特異的遺伝子発現」を捉えることが可能な場合が多く、ヒトを含む哺乳類、鳥類、は虫類、両生類、さらに多くの植物が該当する。以下、このような場合をヒトを対象として便宜上「父方と母方」と表現するが、これを「雄と雌」としてより広い生物種に適用可能である。Ｋ数が３以上の生物は、一部の植物や魚類等の自然界に認められるものの他、人工的に作出されたものが含まれる。
　「リード情報」は、高性能シークエンサで得られる検体の遺伝子断片の塩基配列情報である。高性能シークエンサは、本発明の実施に用いることができる、大量のリード情報を比較的短期間で提供することができるシークエンサであり、いわゆる「次世代シークエンサ」を含むものである。本発明を行った時点における次世代シークエンサとしては、例えば、Genome Sequencer FLX（Roche(454)社）、Genome Analyzer IIx、HiSeq2000、HiSeq2500、MiSeq (共にIllumina社)、SOLiD (Applied Biosystem社)、PacBio RS II (Pacific Biosciences 社)等が挙げられるが、これらに限定されるものではなく、現在、将来において提供される高性能シークエンサの全てを含むものである。リード情報には、概ね、シングルエンドのリード情報と、ペアエンドのリード情報の２種類の形式が認められる。シングルエンドのリード情報とは、リードに対応するＤＮＡ断片の塩基配列の片端の一定長又は可変長（概ね５０～３００ｂｐ程度）についてのリード情報であり、ペアエンドのリード情報とは、当該ＤＮＡ断片の両端の一定長又は可変長（概ね５０～３００ｂｐ程度）についてのリード情報である。技術の進歩に応じてリード情報の内容も日進月歩であるが、本発明においては現在又は将来提供されるリード情報を適用させることが可能である。
　本発明における「マッピング」は、「アライメント」と同意義である。
　上記ステップ（１）は、期待マッピング数の数値化が行われるステップである。「期待マッピング数」は、リード毎に各アリルに対して定義されるマッピング数である。
　上記ステップ（２）は、上記期待マッピング数に基づき合計期待マッピング数が算出されるステップである。「合計期待マッピング数」は、アイソフォーム毎に定義されるマッピング数である。
　上記ステップ（３）は、本発明の推定方法の目的パラメータである「アイソフォーム量比」と「アリル好選度」が算出されるステップである。
　上記ステップ（４）～（６）は、本発明の推定方法におけるループ処理に関するステップで、上記目的パラメータの収束化が行われる。　本発明の最適化方法においては、ステップ（１）と（４）における期待マッピング数と、ステップ（３）と（５）におけるアイソフォーム量比とアリル好選度、のそれぞれが算出される際に、最尤推定を行うためにExpectation Maximizationアルゴリズム（ＥＭアルゴリズム）が使用されることが好ましく、更に好ましくは、ベイズ推定を行うために変分ベイズ法が使用される。さらに本発明の推定方法における推定手法としては、他に、スパースベイズ法、Gibbs sampling法、ＭＣＭＣ法、ＥＰ法、ＰｏｗｅｒＥＰ法等が挙げられる。
［本発明の推定方法における推定要素］
　本発明の推定方法における推定要素として、観測データ、目的パラメータ、及び、観測データと目的パラメータを結びつける潜在変数が挙げられる。これらの推定要素を用いて本発明の推定方法を表現すると、下記のようになる。
　すなわち、本発明の推定方法は、「アリル混在データにおけるリード全体の塩基配列を観測データＲとして、各リードにおける各アリル由来のアイソフォームに対する期待マッピング数を求めるステップ、目的パラメータであるアイソフォーム量θt（ｔ＝１，…，Ｔ：Ｔはアイソフォームの種類数）及び当該アイソフォームｔのアリル選好度φtk（φtk：ｋ＝１，…，Ｋ）の推定値を求めるステップ、並びに、当該推定値に基づいて、アリル特異的遺伝子発現頻度の決定を行うステップ、を含むアリル特異的遺伝子発現頻度のコンピュータによる推定方法において、
　上記目的パラメータθtとφtk、及び、観測データＲを媒介する３種の潜在変数として、（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔn、（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈn、及び、（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳnについて、
　これらの３種の潜在変数を、リードｎの塩基配列を観測データＲnとして、観測データＲnからの目的パラメータθtとφtkの推測工程において、観測データＲnが依存するように繰り入れて、当該推定値を算出することを特徴とする、アリル特異的遺伝子発現頻度の推定方法。　」、として表現される。
　なお、本明細書、特許請求の範囲、及び、図面に記載されたパラメータや数等を示す記号は、発明の開示の便宜ための記号であり、本発明はそれらの記号の種類に全く限定されない。
［１］観測データ
　観測データＲnは、上記の通りに「アリル混在データにおける、ｃＤＮＡリードｎのヌクレオチド配列」である。
　ヒトを例にとると、父母アリル混在データは、検体のＲＮＡシークエンスによるｃＤＮＡのリードを、これとは別に、検体のゲノム配列から推定される父方のｃＤＮＡ塩基配列と母方のｃＤＮＡ塩基配列をリファレンス配列としてマッピングを行って得られたリード個別の情報の総和として提供されるデータである。当該リファレンス配列は、検体のゲノムＤＮＡのリードシークエンスからの推定により得ることも可能である。
　ここでリファレンスとして用いられる父方のｃＤＮＡ配列と母方のｃＤＮＡ配列は、既知と仮定されるが、父方と母方のｃＤＮＡ配列は平均して９９．９％以上同じ配列であり、マッピングを行う個々のｃＤＮＡのリード配列が、父方に属するのか、母方に属するのか、をその配列を直接検討するだけで識別することは極めて困難である。本発明では、このようなアリル混在データから、検体提供者のアリル特異的遺伝子発現頻度を推定することができる。
　［１］－１：検体のＲＮＡシークエンスデータ
　ＲＮＡ－シークエンスは、検体から得られたＲＮＡを基にｃＤＮＡの配列情報とするためにシークエンサにおいて行われる。
　検体は、上記のＫ倍数体生物、さらに具体的には、ヒト、サル、イヌ、ネコ、ウマ、ウシ、ヒツジ等の哺乳類をはじめ、魚類、は虫類、両生類、昆虫類、種子植物類、裸子植物類、蘚苔類等の全ＲＮＡを抽出可能な検体である。検体は、例えば、ヒトを含む哺乳類であれば、血液等の体液や、臓器等の生体組織等、魚類、は虫類、両生類、昆虫類であれば、体液や臓器の他、卵を用いることが可能であり、植物類（蘚苔類を含む）であれば、植物体、種子、胞子等であり、必要であれば適切な保存処理が行われる。典型的な工程として、当該検体から全ＲＮＡを抽出し、さらにポリ（Ｔ）オリゴビーズ等を用いて得られたポリＡ－ＲＮＡ、に対して逆転写を行って一本鎖のｃＤＮＡを得て、通常はこれを断片化して、さらにランダムプライマーを用いて二本鎖とし、これに対して所定の操作によりＰＣＲ等によるＤＮＡ増幅が行われる。次いで、当該ＰＣＲ増幅産物のシークエンスが行われ、「ＲＮＡシークエンスにより得られたリードのヌクレオチド配列」が得られる。上述したように、当該ヌクレオチド配列は、断片の片側だけからシークエンスがなされる「シングルエンドリード」由来の場合と、断片の両側からシークエンスがなされる「ペアエンドリード」由来の場合があり、本発明はこれら双方のリードに対して適用が可能である。これらのｃＤＮＡリード情報は、コンピュータにおいて用いることが可能な電子情報、具体的には、シークエンス配列と、シークエンスの読み取り精度を表すクオリティー値が記述された、「ＦＡＳＴＡ形式」又は「ＦＡＳＴＱ形式」等で提供される。これが「検体のＲＮＡシークエンスデータ」である。
　［１］－２：アリル毎のｃＤＮＡ配列
　リファレンス配列として用いられるアリル毎のｃＤＮＡ配列は、Ｋ倍数体ゲノム配列から推定されて得られる配列である。アリル毎のｃＤＮＡ配列は、ヒトに例を取ると、父方又は母方のアリルのｃＤＮＡ配列を意味している。
　当該ｃＤＮＡ配列は、検体提供源のＫ倍数体ゲノム配列が得られる場合は、当該ゲノム配列に対して直接推定アルゴリズムを適用することにより得ることができる。また、Ｋ倍数体のゲノム配列が得られない場合であっても、通常のゲノムＤＮＡのリード情報からＫ倍数体のゲノム配列を推定するステップを実行することにより得ることができる。リファレンス配列の源は、上記のＲＮＡシークエンスデータの源と同一であることが極めて好ましい。ゲノムＤＮＡを取得するための検体は、ゲノムＤＮＡを抽出可能なものであれば特に限定されず、ヒトを含む哺乳類であれば、血液等の体液、臓器等の生体組織のみならず、毛髪、爪等であっても可能である。魚類、は虫類、両生類であれば、体液や臓器の他、卵を用いることが可能であり、植物類（蘚苔類を含む）であれば、植物体、種子、胞子等である。ゲノムＤＮＡは、所定の操作によりＰＣＲ等の遺伝子増幅が行われ、当該遺伝子増幅産物のＤＮＡ断片が得られる。次いで、当該ＰＣＲ増幅産物のシークエンスが行われる手法もしくは行われない手法によって全ゲノムのシークエンスデータが、リード情報として得られる。上述したように、リードのヌクレオチド配列は、断片の片側だけからシークエンスがなされる「シングルエンドリード」由来の場合と、断片の両側からシークエンスがなされる「ペアエンドリード」由来の場合があり、本発明はこれら双方のリードに対して適用が可能である。これらの全ゲノムリード情報は、コンピュータにおいて用いることが可能な電子情報、具体的には、シークエンス配列と、シークエンスの読み取り精度を表すクオリティー値が記述された、「ＦＡＳＴＡ形式」又は「ＦＡＳＴＱ形式」等で提供される。この全ゲノムリード情報に対して、リファレンス配列とのマッピングが行われ、マッピング情報と深度等がＳＡＭ(sequence alignment/map format)、又は、これをバイナリー（２進数）化したＢＡＭ形式で提供される。これにさらに変異コールとフェージングを行うことでＶＣＦ(Variant Call Format)が得られる。ヒトの場合、さらにこのＶＣＦに対して、ディプロタイプのゲノム構築の推測アルゴリズムを適用することにより、所望の推定ディプロタイプゲノム配列を父方と母方、双方において得ることができる。ここで用いられるステップの一態様については、後述する。
　上記［１］－１と［１］－２のシークエンスは、上述した「高性能シークエンサ」において行うことができる。
　［１］－３：アリル混在データ
　アリル混在データは、（１）－１により得られた検体のＲＮＡシークエンスデータを、（１）－２より得られたアリル毎（ヒトであれば、父方と母方のアリル）において推測されたｃＤＮＡ配列に対して、全てのアリルに対して同時にマッピングを行うことにより、マッピング情報が付加されて、ＳＡＭ形式、又は、ＢＡＭ形式により提供される。マッピングは、全ての可能性を含める形式で行うこと、すなわち、一つのリード情報に対して複数のマッピングを許容する形式で行うこと、が好適である。
　［１］－４：観測データ
　観測データＲnは、上記の通りにＮ個（Ｎは自然数：リードの本数換算）のアリル混在データのうちのｎ番目のリードデータにおけるヌクレオチド配列である。Ｎ個の独立した一様に分布したリードデータであると仮定される。シングルエンドリードの場合は一本のリードに対して一つの観測データＲnが当て嵌められるが、ペアエンドリードの場合には、一本の断片に対して両端のヌクレオチド配列に対応した２つの観測データが組として当て嵌められる。すなわち、ペアエンドリードの場合は、例えばＲnaとＲnbの組を構成するが、これらは同じ断片から由来する塩基情報であるため、このリードの組を一つの単位として扱うことによって、シングルエンドリードの場合と同様の統計的なモデルのもとで扱うことが容易に可能である（Nariai et al., Bioinformatics:15;29(18) 2013）。具体的には、当該先行文献に基づいた後述の［ペアエンドモデルの完全尤度］の記載の通りに計算対応が可能である。
［２］　目的パラメータ
　目的パラメータは、上記の観測データＲnを基に推定がなされるパラメータである。本発明では、２つの目的パラメータを伴っている。
　［２］－１：目的パラメータθt
　目的パラメータθtは（ｔ＝１，…，Ｔ：Ｔはアイソフォームの種類数）、アイソフォーム t の量を表すスカラーである。アイソフォームとは、同じ遺伝子から転写されるが、異なるｍＲＮＡ配列を持つ転写産物であり、アイソフォームの種類数とは、アイソフォームを全ての遺伝子について総計した合計数であり、本発明では、対応するｃＤＮＡの塩基配列を意味する。アイソフォームの違いは、ゲノムＤＮＡ配列におけるＳＮＰ等の塩基の変異による場合と、ｍＲＮＡの選択的スプライシングによる場合があるが、本発明においては後者とする。各アイソフォームについての存在度の分数を、
Figure JPOXMLDOC01-appb-M000001
 　の制約の下で示すことができる。この場合、アイソフォームはＴ個（Ｔは自然数）存在すると仮定され、個々のアイソフォームはｔ（ｔ＝１，…，Ｔ）としてカウントされる。
　［２］－２：目的パラメータφtk
　目的パラメータφtkは、アリルの好選度を表すスカラー変数である。例えば、φtk（φtk：ｋ＝１，…，Ｋ）は、各アイソフォームｔについての各アリルの割合を、0 ≦φtk≦1、
Figure JPOXMLDOC01-appb-M000002
 　の制約の下で示すことができる。
　［２］－３：目的パラメータとアリル特異的遺伝子発現頻度
　上記の２種類の目的パラメータを推定することにより、本発明の目的であるアリル特異的遺伝子発現頻度の推定を行うことができる。例えば、ヒト等の２倍体生物の場合、パラメータは各ｔについてφt1及びφt2となるが、φt1＋φt2＝１の制約から、パラメータはφt1及び１－φt1と表すことができる。以下、φt1をφtと省略する。φt＝０．５よりも１に近いか、０に近いか、によってアイソフォームｔが父方アリル特異的か、母方アリル特異的か、を推定することができる。具体的には、上記の設定でφtを「父方アリルの割合」とした場合に、１に近ければ当該アイソフォームは父方アリル特異的であり、０に近ければ母方アリル特異的である、と推定される。また、上記の目的パラメータθtの推定により、現実に発現している当該父方又は母方のアイソフォームの量を併せて推定することができる。
［３］潜在変数
　本発明における潜在変数は、上記観測データＲnがどのアイソフォームから生成されたか、どちらのアリルから生成されたか、アリル特異的なアイソフォームのどの場所から生成されたかを記述するため繰り入れられる非観測変数である。本発明においては、少なくとも上記の３種の潜在変数［（ａ）リードｎのアイソフォーム選択に関するＴn、（ｂ）リードｎのハプロタイプ選択に関するＨn、及び、（ｃ）リードｎの開始位置に関するＳn］を繰り入れてパラメータθtとφtkを算出推定することで、的確にこれらの目的変数の推定を行い、さらにアリル特異的遺伝子発現頻度の推定を行うことができる。これらの３種の潜在変数を、上記観測データＲnからの目的パラメータθtとφtkの推測工程に、観測データＲnが依存するように繰り入れて、パラメータθtとφtkを算出推定することで、さらにアリル特異的遺伝子発現頻度の推定を的確に行うことができる。
　［３］－１：潜在変数Ｔn
　上記したように潜在変数Ｔnは、リードｎのアイソフォーム選択に関する、上記θtに依存する変数である。Ｔn＝ｔは、リードｎがアイソフォームｔから発生することを意味している。
　［３］－２：潜在変数Ｈn
　上記したように潜在変数Ｈnは、リードｎのハプロタイプ選択に関する、上記φtkとＴnに依存する変数である。ヒト等の２倍体生物の場合、Ｈn＝０は、リードｎが父方アリルから発生していることを意味し、一方Ｈn＝１は、リードｎが母方アリルから発生していることを意味するように設定することができる。
　［３］－３：潜在変数Ｓn
　上記したように潜在変数Ｓnは、リードｎの開始位置に関する、上記Ｔn及びＨnに依存する変数である。Ｓn＝ｓは、リードｎが、位置ｓ（１≦ｓ≦ｌth－Ｌ＋１）（ｌthは、ハプロタイプｈのアイソフォームｔの長さであり、Ｌはリード長である）から発生していることを意味している。ここで、一般的にアイソフォームの長さｌthは、数百塩基長から数千塩基長であり、リードの塩基長よりも長いことが一般的である。また、開始位置ｓが１とは、アイソフォームの最初の塩基からリードが読まれたことを意味する。言い換えればＳnは、リードをリファレンスｃＤＮＡにマッピングした際の、リファレンスにおける開始位置のことを意味している。
［シングルエンドモデルの場合］
　図１は、上述した本発明においてパラメータと変数同士の依存関係を図示した、シングルエンドのリード情報を用いたモデルにおける有向グラフィカルモデルである。
　データの完全尤度（事後同時分布）は、条件付き確率の積として分解される。具体的には、下記式（１）により表される。各記号は、特に断らない限り、前記した通りである。一般的には、θ＝（θ1，…，θT）’をＴ次元ベクトル、φをＴ×Ｋの行列（要素はφtk）と定義する。
　以下ではＫ＝２の場合について詳細を記述する。φはＴ×２の行列となるため、各ｔについてφt1及びφt2がパラメータとなるが、φt1＋φt2　＝１　の制約から、パラメータはφt1及び１－φt1と表すことが出来る。以下、φt1をφtと省略し、φt及び１－φt　をアリル選好度のパラメータとして表す。また、以下ではリードｎが父方ハプロタイプから生成された場合はＨn＝０、母方ハプロタイプから生成された場合はＨn＝１と定義する。
Figure JPOXMLDOC01-appb-M000003
　式（１）において、
　ｐ（Ｔn＝ｔ｜θ）は、θが所与のもと、リードｎがアイソフォームｔから発生する確率である。この確率は、ｐ（Ｔn＝ｔ｜θ）＝θtとして計算され得る（（１）（ａ））。
　ｐ（Ｈn＝ｈ｜Ｔn＝ｔ，φ）は、アイソフォーム選択とφが所与のもと、リードｎがハプロタイプｈ（父方又は母方のいずれか）から発生する確率である。この確率は、ｐ（Ｈn＝０｜Ｔn＝ｔ，φ）＝φt（リードｎが父方アリルから発生する場合）、又は、ｐ（Ｈn＝１｜Ｔn＝ｔ，φ）＝１－φt（リードｎが母方アリルから発生する場合）として計算され得る（（１）（ｂ））。
　ｐ（Ｓn＝ｓ｜Ｔn＝ｔ，Ｈn＝ｈ）は、アイソフォームｔとハプロタイプｈが所与のもと、リードｎが位置ｓから発生する確率である。この確率は、ｐ（Ｓn＝ｓ｜Ｔn＝ｔ，Ｈn＝ｈ）＝１／（ｌth　－Ｌ＋１）として計算され得る（（１）（ｃ））。
　ｐ（Ｒn｜Ｔn＝ｔ，Ｈn＝ｈ，Ｓn＝ｓ）は、アイソフォーム選択、ハプロタイプ選択、及び、リードｎの開始位置が所与のもと、リードｎのヌクレオチド配列を観測する確率である。ここで、Ｔn、Ｈn及びＳnを要約するために、指標変数Ｚnthsを導入することが好適である。Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合、１に等しく、さもなければゼロである。仮に、πnを、リードｎの可能なアライメントについての全（ｔ，ｈ，ｓ）組のセットとして、その時、各（ｔ，ｈ，ｓ）∈πnについて、下記式（２）：
Figure JPOXMLDOC01-appb-M000004
 （式中、subst ( , , )は、１からシークエンスの置換エラーを引いた数値を取るベースクオリティスコア依存置換確率関数であり、ｒn［ｘ］は、リードｎの位置ｘの塩基文字であり、ｑn［ｘ］は、リードｎの位置ｘのベースクオリティスコアであり、ｃt［ｘ］は、ハプロタイプｈのアイソフォームｔのＤＮＡ配列の位置ｘの塩基文字である）によってリード配列の確率を計算することができる。ベースクオリティスコア置換確率関数,「subst ( , , )」 は、Phredベースクオリティスコアにしたがって決定することも可能であり、DNA－Seqデータからリードの参照ＤＮＡ配列に対する最も良いアラインメントから見積もることもできる。なお、Phredベースクオリティスコアは、高性能シークエンサからＦＡＳＴＱフォーマットとして出力される塩基配列情報と共に提供される塩基読み取り精度の目安となるスコア、すなわち、シークエンサが出力するエラー率を示すスコアである(Phred quality score)。具体的には、当該スコアＱは、
　Ｑ＝－１０log10Ｙ（Ｙは、エラー率）、で表される。
　式（１）で示した、シングルエンドモデルにおける本発明の推定方法の完全尤度の数式は、以上のように解釈される。この数式（１）は、シングルエンドモデルにおける本発明の推定方法に係る潜在変数の事後確率や事後分布を求める基礎となるものである。
　また、上記式（２）は、シークエンスの置換エラーを考慮した計算式であるが、シークエンスの挿入・欠失エラーを考慮した計算式も、同様に容易に導出可能である（Nariai et al., Bioinformatics:15;29(18)）。
［ペアエンドモデルの場合］
　ペアエンドのリード情報を用いた場合の完全尤度（事後同時分布）は、例えば、下記式（３）により表される。各記号は、特に断らない限り、前記した通りである。
Figure JPOXMLDOC01-appb-M000005
　ここで、Ｆnはペアエンドリードの組「ＲnaとＲnb」のリファレンス配列へのマッピングから推測される塩基断片（フラグメント）の長さである。
　式（３）右辺において、
　ｐ（Ｔn＝ｔ｜θ）は、θが所与のもと、リードｎがアイソフォームｔから発生する確率である。この確率は、ｐ（Ｔn＝ｔ｜θ）＝θtとして計算され得る（（３）（ａ））。
　ｐ（Ｆn＝ｆ｜Ｔn＝ｔ，Ｈn＝ｈ）は、アイソフォームｔ、ハプロタイプｈが所与のもと、塩基断片の長さｆが発生する確率である。ｄF（ｘ）を、事前に与えられている塩基断片の長さの分布とすると、この確率は
Figure JPOXMLDOC01-appb-M000006
　ここでｌthは、アイソフォームｔのハプロタイプｈの参照配列の長さである。塩基断片の長さの分布ｄF（ｘ）は、例えば、平均μF、標準偏差σFの正規分布として与えられる。平均μF、標準偏差σFは、塩基断片を作成した際において塩基断片長の分布が実験的に分かっていればその値を指定しても良いが、ペアエンドリードのマッピングの結果からこれらのパラメータを推定しても良い（（３）（ｂ））。
　ｐ（Ｈn＝ｈ｜Ｔn＝ｔ，φ）は、アイソフォームｔとアリル好選度φが所与のもと、リードｎがハプロタイプｈから発生する確率である。
　例えば、ヒト等の２倍体生物の父母を考えると、ｐ（Ｈn＝０｜Ｔn＝ｔ，φ）＝φt（リードｎが父方アリルから発生する場合）、又は、ｐ（Ｈn＝１｜Ｔn＝ｔ，φ）＝１－φt（リードｎが母方アリルから発生する場合）として計算され得る（（３）（ｃ））。
　ｐ（Ｓn＝ｓ｜Ｔn＝ｔ，Ｈn＝ｈ，Ｆn＝ｆ）は、アイソフォームｔとハプロタイプｈが所与のもと、断片長ｆのリードｎの組が位置ｓから発生する確率である。この確率は、
　ｐ（Ｓn＝ｓ｜Ｔn＝ｔ，Ｈn＝ｈ，Ｆn＝ｆ）＝１／（ｌt　－ｆ＋１）として計算され得る。ｌtはアリルｔの参照配列の長さ、Ｌ はリード長、ｆは塩基断片の長さを表す（（３）（ｄ））。
　ｐ（Ｒna｜Ｔn＝ｔ，Ｈn＝ｈ，Ｓn＝ｓ）及びｐ（Ｒnb｜Ｔn＝ｔ，Ｈn＝ｈ，Ｓn＝ｓ，Ｆn＝ｆ）は、アイソフォーム選択、リードｎの組の開始位置、断片長が所与のもと、下記式（３－１，３－２）で計算される：
Figure JPOXMLDOC01-appb-M000007
　式中、subst( , , ) は、１からシークエンスの置換エラーを引いた数値を取るベースクオリティスコア依存置換確率関数であり、ｒna［ｘ］は、リードｎの組の一つ目の塩基配列Ｒnaの位置ｘの塩基文字であり、ｑna［ｘ］は、リードｎの組の一つ目の塩基配列位置ｘのベースクオリティスコアであり、ｃtha［ｘ］は、アイソフォームｔのハプロタイプｈのＤＮＡ配列について、リードｎの組の一つ目の塩基配列とマッピングされた位置ｘの塩基文字であり、ｒnb［ｘ］は、リードｎの組の二つ目の塩基配列Ｒnbの位置ｘの塩基文字であり、ｑnb［ｘ］は、リードｎの組の一つ目の塩基配列位置ｘのベースクオリティスコアであり、ｃthb［ｘ］は、アイソフォームｔのハプロタイプｈのＤＮＡ配列について、リードｎの組の二つ目の塩基配列とマッピングされた位置ｘの塩基文字である。
　上記式（３－１）、（３－２）は、ペアエンドモデルによる本発明の推定方法に係る潜在変数の事後確率や事後分布を求める基礎となるものである。
　また、上記式（２）、及び、式（３－１）、（３－２）は、シークエンスの置換エラーを考慮した計算式であるが、シークエンスの挿入・欠失エラーを考慮した計算式も、同様に容易に導出可能である（Nariai et al., Bioinformatics:15;29(18)）。後述のように、本発明は、本発明の推定方法の他、これを実現するためのコンピュータシステム、及び、本発明の推定方法をコンピュータにおいて実現可能なアルゴリズムが含まれるコンピュータソフトウエアを提供する。
　本発明により、ＲＮＡシークエンスのリードデータから、コンピュータによってアリル特異的遺伝子発現頻度を高精度に推定することが可能な方法が提供され、さらに、これを実現可能なコンピュータシステム、及び、コンピュータにおいて実現可能なアルゴリズムが含まれるコンピュータソフトウエア、が提供される。
シングルエンドモデルの本発明の推定方法を図示した有向グラフィカルモデルである。 ヒト検体のアリル特異的遺伝子発現頻度を推定するための本発明の推定方法の処理の流れを示すフローシートである。 ＥＭ法を用いた本発明の推定方法の処理の流れを示すフローシートである。 変分ベイズ法を用いた本発明の推定方法の処理の流れを示すフローシートである。 帰無仮説及びＡＳＥ仮説についての父系／母系比の真の分布（最上部左及び最上部右）、帰無仮説及びＡＳＥ仮説についての既存アプローチでの予測された分布（中央部左及び中央部右）、及び帰無仮説及びＡＳＥ仮説についての本発明の推定方法により予測された分布（最下部左及び最下部右）を示している。 常染色体遺伝子（上部左）、ｃｈｒ１上の遺伝子（上部右）、ｃｈｒ２２上の遺伝子（下部左）、及びｃｈｒＸ上の遺伝子（下部右）についての、本発明の推定方法で予測された分布を示している。 アリル特異的遺伝子発現（ＡＳＥ）を伴う常染色体アイソフォームの総存在量を、ＡＳＥを伴わない常染色体アイソフォームの総存在量と比較した結果を示している。
［本発明の推定方法］
　本発明の推定方法を行うための具体的な手法は限定されず、例えば、最尤推定法、ベイズ推定法が挙げられる。上述したように、最尤推定法に基づくパラメータ推定アリゴリズムとしてはＥＭ法等が挙げられ、ベイズ推定法に基づくパラメータ推定アルゴリズムとしては、スパースベイズ推定法、変分ベイズ推定法、Gibbs sampling法、ＭＣＭＣ法、ＥＰ法、ＰｏｗｅｒＥＰ法等が挙げられる。これらの中でも、ＥＭ法又は変分ベイズ推定法が好適である。言い換えれば、上述した有効グラフィカルモデルで表される本発明の推定法においては、比較的計算の負担が重くならないＥＭ法や変分ベイズ法を用いても、的確に目的のアリル特異的遺伝子発現を推定することが可能である。
　以下、このＥＭ法と変分ベイズ法を用いた本発明の推定方法の態様を説明する。用いられる記号は、特に断らない限り前述した通りである。
［１］ＥＭ（expectation-maximization）法
　ＥＭ法は、潜在変数が存在する確率モデルにおいて最尤推定値を求めるための手法である。本発明の推定方法においては、例えば、下記の内容で行うことができる。
　すなわちＥＭ法による本発明の推定方法は、下記（ａ）～（ｆ）のステップを行うことを特徴とする本発明の推定方法として例示される（この方法を「本発明のＥＭ法」ともいう）。
　（ａ）所与されたθtの初期値、及び、φtkの初期値に基づいて、Ｚnths＝１の事後確率の第１の更新値を算出し、さらに、当該Ｚnths＝１の事後確率の第１の更新値に基づいてθt、及び、φtkの最尤推定値の第１の更新値を算出するステップ、あるいは、（ｂ）所与されたＺnths＝１の事後確率の初期値に基づいて、θt、及び、φtkの最尤推定値の第１の更新値を算出するステップ、
　（ｃ）　直前のステップ（ｄ）（ただし、初回はステップ（ａ）又はステップ（ｂ）である）により算出されたθt、及び、φtkの最尤推定値の更新値に基づいて、新たにＺnths＝１の事後確率の更新値を算出するステップ、
　（ｄ）　ステップ（ｃ）において算出されたＺnths＝１の事後確率の更新値に基づいて、θt、及び、φtkの最尤推定値の更新値を新たに算出するステップ、
　（ｅ）　ステップ（ｄ）において算出された最尤推定値の更新値の対数尤度を算出するステップ、
　（ｆ）　（ｉ） ステップ（ｅ）において算出された対数尤度の収束性を評価するステップ、
　　　　　（ii） ステップ（ｃ）で算出されたＺnths＝１の事後確率の更新値の収束性を評価するステップ、
　　　　　（iii） ステップ（ｄ）で算出されたθt、及び、φtkの最尤推定値の更新値の収束性を評価するステップであって、
　収束が認められれば、当該収束値をθtとφtkの最終推定値として決定し、収束が認められなければ、ステップ（ｃ）、（ｄ）、（ｅ）、及び、（ｆ）の繰り返しを決定する。
　本発明のＥＭ法の（ａ）ステップと（ｂ）ステップは、繰り返し最適化を行う本法の初期値を定義するステップである。
　Ｚnthsの初期値は、ｎについては１からＮの各々の数字であり、ｔについては１からＴについての各々の数字であり、ｈについては２倍体の場合は０もしくは１と定義する。Ｋ倍数体の場合、ｈは１からＫである。ｓについては１≦ｓ≦ｌth－Ｌ＋１、ｌthは、ハプロタイプｈのアイソフォームｔの長さ、Ｌはリード長である。そして、あるリードが複数箇所にアライメントしている場合、同様に確からしいと仮定して、Ｚnthsの初期値は１／［リードｎがアライメントしている合計箇所］とすることが好ましい。θtの初期値は、各々のアイソフォームの発現量は等しいと仮定して、１／Ｔとすることが好ましく、φtkの初期値は、各アリル由来である確率が確からしいと仮定して１／Ｋとすることが好ましい。ヒトを含む２倍体生物（Ｋ＝２）の場合は、父アリル由来、母アリル由来である確率が同様に確からしいと仮定して１／２とすることが好ましい。なお、被験対象の生物の母集団が特定のカテゴリーに属する場合等（例えば、被験対象がヒトであり、かつ、特定の人種の集団である等）、これらのデータの事前情報が既知の場合には、それらを上記θtとφtkの初期値として用いることも可能である。
　（ｃ）ステップは、現在の推定値（最初は（ａ）で得られた初期化の値）に基づいて、Ｚnths＝１の事後確率ｒnthsを算出するステップであり、本発明のＥＭ法のＥステップに相当する。
　一般的には、θ＝（θ1，…，θT）’をＴ次元ベクトル、φをＴ×Ｋの行列（要素はφtk）と定義する。
　以下ではＫ＝２の場合について詳細を記述する。φはＴ×２の行列となるため、各ｔについてφt1及びφt2がパラメータとなるが、φt1＋φt2＝１　の制約条件の下、パラメータはφt1及び１－φt1と表すことが出来る。以下、φt1をφtと省略し、φt及び１－φtをアリル選好度のパラメータとして表す。また、以下ではリードｎが父方ハプロタイプから生成された場合はＨn＝０、母方ハプロタイプから生成された場合はＨn＝１　と定義する。
　すなわち、Ｐ（Ｚnths＝１｜θ*，φ*）を算出するステップであり（＊は更新されたことを表している。以下、同様。）、
Figure JPOXMLDOC01-appb-M000008
 として計算される。
　ここで、シングルエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000009
 （ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（１）（ａ）に、第２項は上記式（１）（ｂ）に、第３項は上記式（１）（ｃ）に、及び、第４項は上記式（２）について、対数計算することにより算出することができる）
　ペアエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000010
 （ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（３）（ａ）に、第２項は上記式（３）（ｂ）に、第３項は上記式（３）（ｃ）に、第４項は上記式（３）（ｄ）に、及び、第５項は上記式３－１と３－２について、対数計算することにより算出することができる）
　ここまでが本発明のＥＭ法のＥステップに相当する。
　（ｄ）ステップは、本発明のＥＭ法のＭステップに相当するステップであり、（ｃ）ステップで得られたＺnths＝１の事後確率ｒnthsに基づき、θtとφtの最尤推定値を求めるステップである。以下の通り、対数尤度を最大にするθtとφtを算出する。
Figure JPOXMLDOC01-appb-M000011
　（ｅ）ステップと（ｆ）ステップは、上記の通りである。収束基準の好適な一例として、存在量パラメータθt＞１０-7であるアイソフォームについての相対的変化１０-3が、収束基準として使用されるが、異なる収束基準を用いることも可能である。
　なお、通常は（ｅ）ステップの「ステップ（ｃ）、（ｄ）、及び、（ｅ）の繰り返し」が、１回以上は行われるが、この繰り返し工程を１回も行わない、又は十分に収束していない状態でステップを終了して、アリル特定遺伝子発現の推定を行うことも可能である。
［２］変分ベイズ法
　変分ベイズ法は、ベイズ推定法においてパラメータの事後確率分布を推定することによって、より的確な推測を行うための方法である。本発明においては、以下の内容で本法が行われる。
　式（１）と式（３）で示した、本発明の推定方法の完全尤度（事後同時分布）の数式は、全ての可能な潜在変数Ｚにわたる積分を伴い、解析的に解けない。そのため、当該完全尤度（事後同時分布）において、潜在的変数及びモデルパラメーターの因子分解を、
Figure JPOXMLDOC01-appb-M000012
と仮定することによって、近似値を求めるものである。
　θの事前分布について、本発明では、ディリクレ分布
Figure JPOXMLDOC01-appb-M000013
 （式中、αt ＞ 0は、ハイパーパラメータであり、
Figure JPOXMLDOC01-appb-M000014
そしてг(・)は、ガンマ関数である）
　を用いる。本発明において、アイソフォーム量の相対的差異について予備知識はないという仮定に基づいて、全てのアイソフォームについて単一のハイパーパラメータα0が用いられることが好適である。当該単一ハイパーパラメータα0は、目的パラメータの複雑さを制御する。α0 ≧１の時、α0－１は、アイソフォームに割り当てられるリードの事前カウントとして解釈され得、α0＜１の時、当該事前分布は、アイソフォーム量のいくつかがゼロの傾向を与える。次いで、対数周辺尤度の下限を最大化するα0が選択される。
　一般的には、θ＝（θ1，…，θT）’をＴ次元ベクトル、φをＴ×Ｋの行列（要素はφtk）と定義する。
　以下ではＫ＝２の場合について詳細を記述する。φはＴ×２の行列となるため、各ｔについてφt1及びφt2がパラメータとなるが、φt1＋φt2＝１　の制約の下、パラメータはφt1及び１－φt1と表すことができる。以下、φt1をφtと省略して、φt及び１－φtをアリル選好度のパラメータとして表す。また、以下ではリードｎが父方ハプロタイプから生成された場合はＨn＝０、母方ハプロタイプから生成された場合はＨn＝１と定義する。
　φの事前分布について、本発明では、β分布
Figure JPOXMLDOC01-appb-M000015
 （式中、βt,1 ＞０及びβt,2 ＞０はハイパーパラメータであり、Ｂ(・，・)は、β関数である。）
　を用いる。ここで、βt,1及びβt,2は、父方／母方比を計算するために、それぞれ父方及び母方アリルに割り当てられるリードの事前カウントともいえる。本発明では、全てのアイソフォームについて、無情報事前分布としてβt,1＝βt,2＝１とすることが好ましい。なお、被験対象の生物の母集団が特定のカテゴリーに属する場合等（例えば、被験対象がヒトであり、かつ、特定の人種の集団である等）、これらのデータの事前情報が既知の場合には、それらを上記θtとφtkの初期値として用いることも可能である。また、アイソフォーム量やアリル好選度の情報が既知の場合には、当該既知情報に応じてα0とβt,1、βt,2を設定することも可能である。さらにＫ＝３以上の場合は、例えば、上記β分布を多変量に拡張したディリクレ分布を用いて、上記Ｋ＝２でβ分布の場合と同様にハイパーパラメータβt,γ（γ＝１，…，Ｋ）を設定することができる。
　ハイパーパラメータα0、βt,1、βt,2が所与のもと、対数周辺尤度の下限は、変分ベイズ推定アルゴリズムによって繰り返しにより最大化される。
　すなわち変分ベイズ法による本発明の推定方法は、下記（ａ）～（ｅ）のステップを行うことを特徴とする本発明の推定方法として例示される（この方法を「本発明の変分ベイズ法」ともいう）。
　（ａ）　所与されたアイソフォームｔの存在量についての予備知識の分布を示すハイパーパラメータαtの初期値に基づくθtの事後分布の更新値、及び、所与された個々のアリルに割り当てられるリード比についての予備知識の分布を示すハイパーパラメータβt,1、βt,2の初期値に基づくφtの事後分布の更新値、に基づいてＺnths＝１の事後分布を算出し、さらに、当該Ｚnths＝１の事後分布に基づいてθt、及び、φtの第１の更新事後分布の更新値を算出するステップ、あるいは、（ｂ）所与されたＺnths＝１の初期分布に基づいてθt、及び、φtの第１の事後分布の更新値を算出するステップ、
　（ｃ）　直前のステップ（ｄ）（ただし、初回はステップ（ａ）又はステップ（ｂ）である）により算出されたθt、及び、φtの事後分布の更新値に基づいて、新たにＺnths＝１の事後分布を算出するステップ、
　（ｄ）　ステップ（ｃ）において算出された　Ｚnths＝１の事後分布を基にして、θtとφtの事後分布を新たに算出して更新するステップ、
　（ｅ）　ステップ（ｄ）において得られたθtとφtの事後分布の期待値の収束性を評価するステップであって、当該期待値における収束が認められれば、当該収束期待値をθtとφtkの推定値として決定し、収束が認められなければ、ステップ（ｂ）、（ｃ）、及び、（ｄ）の繰り返しを決定する。
　本発明の変分ベイズ法の（ａ）ステップと（ｂ）ステップは、繰り返し最適化を行う本法の初期値を定義するステップである。具体的には、各アイソフォームｔについて、ｑ*（θ）のパラメータであるα* tについて、各々のアイソフォームの発現量は等しいと仮定して、
Figure JPOXMLDOC01-appb-M000016
 を設定することが好ましい。
　また、ｑ*（φ）のパラメータであるβ*t,1、β*t,2について、父アリル由来、母アリル由来である確率が同様に確からしいと仮定して、
Figure JPOXMLDOC01-appb-M000017
 を設定することが好ましい。
　（ｃ）ステップは、ｑ*（θ）及びｑ*（φ）の現在の推定値が所与のもと、Ｅz［Ｚnths］を計算するステップである。Ｚnthsの値は、ｎについては１からＮの各々の数字であり、ｔについては１からＴについての各々の数字であり、ｈについては０もしくは１であり、ｓについては１≦ｓ≦ｌth－Ｌ＋１、ｌthは、ハプロタイプｈのアイソフォームｔの長さ、Ｌはリード長である。Ｅz［Ｚnths］は、ｑ*（θ）及ｑ*（φ）の現在の推定値に基づいて、
Figure JPOXMLDOC01-appb-M000018
 として計算される。
　ここで、シングルエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000019
 である。
（ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（１）（ａ）に、第２項は上記式（１）（ｂ）に、第３項は上記式（１）（ｃ）に、及び、第４項は上記式（２）について、対数計算することにより算出することができる）
　ペアエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000020
 （ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（３）（ａ）に、第２項は上記式（３）（ｂ）に、第３項は上記式（３）（ｃ）に、第４項は上記式（３）（ｄ）に、及び、第５項は上記式３－１と３－２について、対数計算することにより算出することができる）
　ここで、
Figure JPOXMLDOC01-appb-M000021
 （式中、ψ(・)はディガンマ関数である）
である。
　（ｄ）ステップは、ｑ*（Ｚ）の現在の推定値が所与のもと、Ｅθ［θt］及びＥφ［φt］を計算するステップである。
　Ｅθ［θt］は、ｑ*（Ｚ）の現在の推定値に基づいて、
Figure JPOXMLDOC01-appb-M000022
 として計算される。
　従って、ｑ*（θ）もまたディリクレ分布であり、事前分布ｐ（θ）は共役事前分布である。
　同様に、Ｅφ［φt］は、ｑ*（Ｚ）の現在の推定値に基づいて、
Figure JPOXMLDOC01-appb-M000023
 として計算される。
　従って、ｑ*（φt）もまたβ分布であり、事前分布ｐ（φt）は共役事前分布である。
　（ｄ）ステップにおいて、収束基準の好適な一例として、存在量パラメータＥθ［θt］＞１０-7であるアイソフォームについての相対的変化１０-3が、収束基準として使用されるが、異なる収束基準を用いることも可能である。
　なお、通常は（ｄ）ステップの「ステップ（ｂ）、（ｃ）、及び、（ｄ）の繰り返し」が、１回以上は行われるが、この繰り返し工程を１回も行わない、又は十分に収束していない状態でステップを終了して、アリル特異的遺伝子発現の推定を行う事も可能である。
　ここで変分の下限についての検討を行い、本発明の変分ベイズ法の更新により得られる収束値が、真の目的変数θとφを近似するものであることを示す。
　本発明の変分ベイズ法における対数周辺尤度は、
Figure JPOXMLDOC01-appb-M000024
 として分解することができる。
　ＫＬ（ｑ｜｜ｐ）は、ｑ（θ，φ，Ｚ）とｐ（θ，φ，Ｚ｜Ｒ）の間のカルバック・ライブラー距離であるので、対数周辺尤度は、Ｌ（ｑ）によって下限を与えられる。すなわち、カルバック・ライブラー距離は常に０以上であるため、Ｌ（ｑ）が当該周辺尤度の下限を構成する。上記の（ｂ）ステップと（ｃ）ステップを繰り返し更新する度に、Ｌ（ｑ）（対数周辺尤度の下限）を増加させることが一般的に示される（Bishop CM. Pattern Recognition and Machine Learning. Springer Science:Business Media, LLC, New York, NY, USA; 2006）。
　上記の因子分解仮定
Figure JPOXMLDOC01-appb-M000025
 を用いて、
Figure JPOXMLDOC01-appb-M000026
 である。
［本発明のコンピュータシステムとコンピュータプログラム］
〔１〕本発明のコンピュータシステム
　本発明のコンピュータシステムは、上述した本発明の推定方法を行う手段となるシステムであり、本発明のプログラムは、本発明のコンピュータシステムに本発明の推定方法を行わせるためのアルゴリズムを備えたコンピュータプログラムである。「アルゴリズム」とは、コンピュータ分野の一般的な概念と同じく、問題を解くための手順を定式化した形で表現したものを意味する。
　本発明のコンピュータシステムは、通常のコンピュータシステムに関わるハードウエアを備えることができる。すなわち、通常ハードディスクドライブが該当する「記録部」、ＣＰＵ（中央演算処理装置：Central Processing Unit）、ＤＳＰ(Digital Signal Processor)、ＲＩＳＫ（Reduced Instruction Set Computer）等が該当する「演算処理部」の他、例えば、メインメモリ等として用いられるＲＡＭ(Random Access Memory)が該当する「一時記憶部」、キーボード、マウス、タッチパネル等が該当する「操作部」、ディスプレイが該当する「表示部」、操作部に応じたシリアル又はパラレルインターフェース等が該当する「出入力インターフェース（ＩＦ）部」、ビデオメモリとＤ／Ａ変換部を備え、表示部のビデオ方式に応じたアナログ信号を出力する「通信インターフェース（ＩＦ）部」を備えている。当該通信ＩＦ部では、外部の情報、特に、ヒトゲノムデータベース等のヒトゲノム情報とデータ交換を行うことができる。
　以下においては特に断らない限り、本発明のコンピュータシステムの「演算処理部」が行う処理として説明する。「演算処理部」は、「操作部」が操作されて「通信ＩＦ部」を介して、特にヒトゲノムデータベースのデータを取得して「記録部」に記録し、適宜当該「記録部」からデータを「一時記憶部」に読み出し、所定の演算処理を行った後、その結果を再度「記録部」に記録する。当該「演算処理部」は、「操作部」の操作を促す画面データや処理結果を表示する画面データを作成し、入力ＩＦ部のビデオＲＡＭを介して、これらの画像を「表示部」に表示する。本発明のプログラムは、用時又は予め「記録部」に記録、あるいは、外部のハードウエア資源に記録されており、必要に応じて「演算処理部」において、記載されたアルゴリズムに従った演算処理が行われる。
〔１〕－１　ＥＭ法を用いたコンピュータシステム
　ＥＭ法を用いた本発明のコンピュータシステムは、下記の通りであり、各記号は、特に断らない限り、前記した通りである。
　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え；
（１）　当該記録部には、当該被験生物の全てのアリル由来のｃＤＮＡのリード情報が、リードの配列及びリードのマッピング先の観測データとして記録されており、
（２）　当該演算処理部では読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２が択一的に実行され、
　（２）－１：アイソフォーム量に関する変数θtの初期値の算出処理、及び、アリル好選度に関する変数φtkの初期値の算出処理、
　（２）－２：上記変数θtとφtk、及び、観測データである当該被験生物の全てのアリル由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数としての下記（ａ）、（ｂ）及び（ｃ）：
　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔn、
　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈn、
　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳn、
が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の初期値の算出処理、
（３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θtとφtkに基づき、当該指標変数Ｚnthsの算出処理がなされ、
（４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsに基づいて変数θtとφtkの最尤推定値の第１の更新値が算出され、
（５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの最尤推定値の第１の更新値に基づいて上記ステップ（３）とステップ（４）が再度実行され、さらに、変数θtとφtkの第２の更新値が算出されるループ処理が、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束した変数θtとφtkが最適化されたθtとφtkのデータとして、上記記録部に記録がなされる；
　ことを特徴とする、個々のアリル由来のｃＤＮＡのリード情報を最適化するコンピュータシステム。
〔１〕－２　変分ベイズ法を用いたコンピュータシステム
　変分ベイズ法を用いた本発明のコンピュータシステムは、下記の通りであり、各記号は、特に断らない限り、前記した通りである。
　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え；
（１）　当該記録部には、被験生物の各アリル由来のｃＤＮＡのリード情報が、リードの配列及びリードのマッピング先の観測データとして記録されており、
（２）　当該演算処理部では読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２が択一的に実行され、
　（２）－１：アイソフォームｔの存在量についての予備知識の分布を示すハイパーパラメータαtの初期値に基づくθtの事後分布の更新値の算出処理、及び、個々のアリルに割り当てられるリード比についての予備知識の分布を示すハイパーパラメータβt,γ（γは、１，…，Ｋの各整数であり、βt,γはＫ個存在する。）の初期値に基づくφtkの事後分布の更新値の算出処理、
　（２）－２：上記θtとφtkの分布、及び、観測データである被験者の父方由来と母方由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数分布としての下記（ａ）、（ｂ）及び（ｃ）：
　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔnの分布、
　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈnの分布、
　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳnの分布、
が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の分布の初期分布の算出処理、
（３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θとφの分布に基づき、当該指標変数Ｚnthsの事後分布の算出処理がなされ、
（４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsの事後分布の更新値に基づいて変数θtとφtkの第１更新事後分布の更新値が算出され、
（５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの第１の更新事後分布に基づいて上記ステップ（３）とステップ（４）が再度実行され、さらに、変数θtとφtkの第２の更新事後分布が算出されるループ処理が、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束したθtとφtkの期待値が最適化されたθtとφtkのデータとして、上記記録部に記録がなされる；
　ことを特徴とする、個々のアリル由来のｃＤＮＡのリード情報を最適化するコンピュータシステム。
〔２〕本発明のコンピュータプログラム
　本発明のコンピュータプログラムは、上記の本発明のコンピュータシステムを実行するためのアルゴリズムを備えたコンピュータプログラムであり、特に断らない限り、各記号や概念は前述した通りである。
〔２〕－１　ＥＭ法によるコンピュータプログラム　ＥＭ法を用いた本発明のコンピュータプログラムは下記の通りである。
　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに、
（１）　記録部に記録されている当該被験生物の各アリル由来のｃＤＮＡのリード情報を、リードの配列及びリードのマッピング先の観測データとして読み出す、第１の機能；
（２）　演算処理部において読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２を択一的に実行する、第２の機能；
　（２）－１：アイソフォーム量に関する変数θtの初期値の算出処理、及び、アリル好選度に関する変数φtkの初期値の算出処理、
　（２）－２：上記変数θtとφtk、及び、観測データである当該被験生物の各アリル由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数としての下記（ａ）、（ｂ）及び（ｃ）：
　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔn、
　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈn、
　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳn、
が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の事後確率（リードｎの期待マッピング数）の初期値の算出処理、
（３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θtとφtkに基づき、当該指標変数Ｚnthsの事後確率（リードｎの期待マッピング数）の算出処理を行う、第３の機能；
（４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsの事後確率に基づいて変数θtとφtkの最尤推定値の第１の更新値の算出処理を行う、第４の機能；
（５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの最尤推定値の第１の更新値に基づいて上記ステップ（３）とステップ（４）を再度実行し、さらに、変数θtとφtkの第２の更新値が算出されるループ処理を、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束した変数θtとφtkを最適化されたθtとφtkのデータとして、上記記録部への記録を行う、第５の機能；
　を実現させるアルゴリズムが含まれていることを特徴とする、各アリル由来のｃＤＮＡのリード情報を最適化するコンピュータプログラム。
〔２〕－２　変分ベイズ法によるコンピュータプログラム
　変分ベイズ法を用いた本発明のコンピュータプログラムは下記の通りである。
　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに、
（１）　記録部に記録されている当該被験生物の各アリル由来のｃＤＮＡのリード情報を、リードの配列及びリードのマッピング先の観測データとして読み出す、第１の機能；
（２）　演算処理部において読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２を択一的に実行する、第２の機能；
　（２）－１：アイソフォームｔの存在量についての予備知識の分布を示すハイパーパラメータαtの初期値に基づくθtの分布の期待値の算出処理、及び、各アリルに割り当てられるリード比についての予備知識の分布を示すハイパーパラメータβt,γ（γは、１，…，Ｋの各整数であり、βt,γはＫ個存在する。）の初期値に基づくφtkの分布の期待値の算出処理、
　（２）－２：上記θtとφtkの分布、及び、観測データである当該被験生物の全てのアリル由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数分布としての下記（ａ）、（ｂ）及び（ｃ）：
　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔnの分布、
　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈnの分布、
　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳnの分布、
が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の分布（リードｎの期待マッピング数）の初期分布の算出処理、
（３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θtとφtkの分布に基づき、当該指標変数Ｚnthsの分布（リードｎの期待マッピング数）の算出処理を行う、第３の機能；
（４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsの分布に基づいて変数θtとφtkの第１更新事後分布の算出処理を行う、第４の機能；
（５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの第１の更新事後分布に基づいて上記ステップ（３）とステップ（４）を再度実行し、さらに、変数θtとφtkの第２の更新事後分布が算出されるループ処理を、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束したθtとφtkの期待値を最適化されたθtとφtkのデータとして、上記記録部への記録を行う、第５の機能；
　を実現させるアルゴリズムが含まれていることを特徴とする、個々のアリル由来のｃＤＮＡのリード情報を最適化するコンピュータプログラム。
　本発明のコンピュータプログラムは、例えば、Ｃ言語、Java（登録商標）、Perl、Python等で記載することが可能である。
　本発明はさらに、本発明のプログラムが記録されていることを特徴とする、コンピュータにおいて読み取り可能な記録媒体又はコンピュータに接続し得る記録媒体（以下、本発明の記録媒体ともいう）を提供する。これらの記録媒体としては、フレキシブルディスク、フラッシュメモリ、ハードディスク等の磁気的媒体、ＣＤ、ＤＶＤ、ＢＤ等の光学的媒体、ＭＯ、ＭＤ等の磁気光学的媒体等が挙げられ、特に限定されるものではない。
　本発明のコンピュータシステムの典型は、本発明のプログラムを実行することを特徴とするものである。
〔本発明の推定方法のさらに具体的な形態〕
　ここでは、アリル特異的遺伝子発現頻度の最適化の対象を２倍体生物であるヒトとした態様を開示する。
　図２は、ヒト検体のアリル特異的遺伝子発現頻度を推定するための本発明の推定方法の処理の流れを示すフローシートである。この処理の流れの各々の要素は、全てコンピュータにおいて、必要に応じてコンピュータネットワークやデータベースを介して、コンピュータソフトウエアのアルゴリズムを実現する命令に基づき、電子的に行われるものである。
　ボックスＢ１－１は、検体のゲノムＤＮＡのリードデータの存在を示しており、ボックスＢ１－２は、ヒトゲノム参照配列の存在を示している。「リードデータ」（ボックスＢ１－１）は、検体のゲノムＤＮＡを高性能シークエンサで処理することにより得られる、個々のゲノムＤＮＡ断片（リード）の一部又は全部の塩基配列であり、シングルエンド形式であっても、ペアエンド形式であってもよい。通常は読み取り精度を示す情報（クオリティースコア）も付加されている（ＦＡＳＴＱ：ＤＮＡの塩基配列を表すＦＡＳＴＡフォーマットに由来する用語）。「ヒトゲノム参照配列」（ボックスＢ１－２）は、一個人のゲノム配列情報に帰着するものであり、その情報源は特に限定されない。検体と当該ゲノム配列情報の人種が異なっており、アイソフォームの遺伝子配列に関して異質性を含んでいることも想定される。参照配列のバージョンの変更に本発明の実質的な効果が依存するものではなく、常にバージョン変更に対応させてマッピングを行うことができる。ヒトゲノム参照配列としては、例えば、ＵＣＳＣ(University of California Santa Cruz)が管理しているヒトゲノム参照配列（ｈｇ１９等）、Genome Reference Consortiumが管理しているヒトゲノム参照配列（ＧＲＣｈ３７等）が挙げられるが、これらに限定されるものではない。
　ステップＳ１は、ゲノムＤＮＡのマッピングを行うステップであり、上記ボックスＢ１－１の「リードデータ」を、ボックスＢ１－２の「ヒトゲノム参照配列」に対してマッピングを行うステップである。このマッピングに関しては、コンピュータで当該マッピングを実現可能なアルゴリズムを含むコンピュータプログラムを作出して行うことも当業者において可能であるが、既に提供されているソフトウエアを用いることも可能である。既存のマッピング用のソフトウエアとしては、例えば、ＢＷＡ－ＭＥＭ(http://bio-bwa. Source. Net/)が挙げられる。
　なお、上記のヒトゲノム参照配列には「おとり配列」が含まれることが好適である。「おとり配列」とは、予め用意された参照配列には存在しないゲノム配列である。おとり配列を用いない場合、ヒトゲノム参照配列には登録されていない配列由来ではないリードが、既存のヒトゲノム参照配列のどこかにマッピングされてしまう危険性が高くなり、マッピングの精度が低下する可能性がある。おとり配列としては、ｈｓ３７ｄ５(ftp://ftp.1000genmes.ebi.ac.uk/voll/ftp/technical/phase2_reference_assembly_sequence/hs37d5.fa.gz)等が例示される。
　ボックスＢ２は、上記Ｓ１における「ゲノムＤＮＡのマッピング結果」の存在を示すボックスであり、マッピングの結果は、ＳＡＭ形式又はＢＡＭ形式等のフォーマットで存在している。
　ステップＳ２は、上記ボックスＢ２の「ゲノムＤＮＡのマッピングの結果」に対して、変異コールとハプロタイプフェージングを行うステップである。変異コールは、ＳＮＰ等の変異を特定する手段であり、ハプロタイプフェージングは、ハプロタイプレベルでコールされた変異を含んだゲノムＤＮＡの塩基配列を推定する手段である。コンピュータで変異コールとフェージングを、別個に又は一緒に実現可能なアルゴリズムを含むコンピュータプログラムを作出して行うことも当業者において可能であるが、既に提供されているソフトウエアを用いることも可能である。既存の変異コールとフェージング用のソフトウエアとしては、例えば、HapMonster（http://nagasakilab.csml.org/ja/hapmonster）等が挙げられる。
　ボックスＢ３は、ステップＳ２で行った変異コールとフェージングの結果得られた、ハプロタイプとして推定されたゲノムＤＮＡの配列情報がＶＣＦ形式で存在することを示している。
　ステップＳ３は、ボックスＢ３のハプロタイプのゲノムＤＮＡの配列情報から、ディプロタイプの配列情報を推定して父方のハプロタイプのゲノムＤＮＡ配列と母方のハプロタイプのゲノムＤＮＡ塩基配列をそれぞれ推定構築するステップである。コンピュータでこの推定構築を、これを実現可能なアルゴリズムを含むコンピュータプログラムを作出して行うことも当業者において可能であるが、既に提供されているソフトウエアを用いることも可能である。既存の当該推定構築用のソフトウエアとしては、例えば、vcf2diploid(Rozowsky et al., Mol Syst Biol. 2011 Aug 2;7:522.)が挙げられる。
　ボックスＢ４（Ｂ４－１とＢ４－２）は、ステップＳ３で得られた父方と母方それぞれのハプロイドゲノム塩基配列の存在を示している。ボックスＢ４－１が父方ゲノムＤＮＡ塩基配列で、ボックスＢ４－２が母方ゲノムＤＮＡ配列である。また、過去の実験結果等により、サンプル特異的な父方と母方が区別されたディプロタイプのゲノム塩基配列が既に得られている場合は、当該ディプロタイプゲノム塩基配列を直接用いて、それ以前のステップを経ないで、このボックスＢ４の段階から次のステップＳ４を行うことができる。
　ステップＳ４は、ボックスＢ４の父方と母方のゲノムＤＮＡ塩基配列を、ｃＤＮＡ塩基配列へと推定するステップである。このステップＳ４の推定は、これを実現可能なアルゴリズムを含むコンピュータプログラムを作出して行うことも可能であるが、既に提供されているソフトウエアを用いることも可能である。既存の当該推定用のソフトウエアとしては、例えば、ｒＳｅｑ（Jiang H, Wong WH. Statistical inferences for isoform expression in RNA-Seq. Bioinformatics.2009;25:1026-1032）等を用いることができる。
　ボックスＢ５（Ｂ５－１とＢ５－２）は、ステップＳ４で得られた父方と母方それぞれのｃＤＮＡ塩基配列の存在を示している。ボックスＢ５－１が父方ｃＤＮＡ塩基配列で、ボックスＢ５－２が母方ｃＤＮＡ配列である。
　ステップＳ５は、ボックスＢ６に示した「検体のＲＮＡシークエンスデータ」(リードデータ)を、ボックスＢ５に示した父方と母方のｃＤＮＡの塩基配列を一緒にしたものを参照配列としてマッピングを行うステップである。リードデータは、シングルエンドデータであっても、ペアエンドデータであってもよい。
　ボックスＢ７は、ステップＳ５で得られたマッピングに基づく電子データ、すなわち観測データＲnとして用いられる「父母アリル混在データ」の存在を示している。
　図２のうち、ここまでがマッピングプロセスを示す部分である。
　以下、ステップＳ６は、上記「父母アリル混在データ」を最適化するステップであり、その詳細を、上記と一部重複したフローシートとして図３及び図４に示す。図３に示す態様は、「父母アリル混在データ」を最適化してアリル特異的遺伝子発現頻度を推定するために、ＥＭアルゴリズムによる最尤推定を行う工程である。図４に示す態様は、「父母アリル混在データ」を最適化してアリル特異的遺伝子発現頻度を推定するために、変分ベイズ法によるベイズ推定を行う工程である。
〔１〕　ＥＭアルゴリズムによる最尤推定
　ステップＳ６－１１は、上記図２に示したステップＳ５及びボックスＢ７の内容を示しており、「父母アリル混在データ」の存在を示している。
　ステップＳ６－１２は、本コンピュータプログラムの実行による最適化処理に際して用いられるパラメータＥ［Ｚnths］（リードｎのアイソフォームｔのアリルｈの位置ｓへの割り当て分のＺnths＝１の期待値）と、θt（リード総量に対してアイソフォームｔに割り当てられたリードの割合）と、φt（アイソフォームｔにおける、父アリルからの発現量の割合）の初期化を行うためのステップである。ここでＮは「リード総数」であり、Ｔは「アイソフォームの種類の総数」である。ｈは、父アリルの場合ｈ＝０であり、母アリルの場合ｈ＝１である。Ｓは「リードｎの開始位置」であり、ｌthは、アイソフォームｔのｈアリルの長さであり、Ｌはリード長である。Ｅ［Ｚnths］は、事前の無情報・一様分布を前提として、１／Ｍn（Ｍnは、観測されたリードｎのアライメント先の総数）として初期設定がなされている。
　θtは、事前の無情報・一様分布を前提として、１／Ｔが初期値として与えられている。
　Φtは、ｒt,1／（ｒt,1＋ｒt,2）（父アリルからの発現量：ｒt,1＝Ｎ／２Ｔ、母アリルからの発現量：ｒt,2＝Ｎ／２Ｔ）が初期値として与えられており、例えば、事前の無情報・一様分布故に、父アリル由来、母アリル由来で同様に確からしいと仮定して１／２とすることが好適である。
　ステップＳ６－１３は、直前に更新された（初回は上記ステップＳ６－１２の初期化値）θtとφtに基づいて、リードｎのアイソフォームｔのアリルｈの位置ｓへの割り当て分のＺnths＝１の期待値：Ｅz［Ｚnths］）を再計算するステップであり、本発明のＥＭ法のＥステップに相当するステップである。
　再計算方法：Ｅz［Ｚnths］= Ｐ（Ｚnths＝１｜θ* t，φ* t）＝ｒnths
において前述のように、
Figure JPOXMLDOC01-appb-M000027
 として計算される。
　ここで、シングルエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000028
 である。
（ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（１）（ａ）に、第２項は上記式（１）（ｂ）に、第３項は上記式（１）（ｃ）に、及び、第４項は上記式（２）について、対数計算することにより算出することができる）
　ペアエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000029
 （ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（３）（ａ）に、第２項は上記式（３）（ｂ）に、第３項は上記式（３）（ｃ）に、第４項は上記式（３）（ｄ）に、及び、第５項は上記式３－１と３－２について、対数計算することにより算出することができる）
　ステップＳ６－１４は、上記ステップＳ６－１３において算出されたＥz［Ｚnths］（ｒnths）に基づき、θtとφtの最尤推定値を算出する、本発明のＥＭ法のＭステップに相当するステップである。以下の通りに、対数尤度を最大にするθtとφtを算出する。具体的には、
アイソフォームｔは、１からＴ（アイソフォームの種類の総数）をとり、
Figure JPOXMLDOC01-appb-M000030
 である。
ステップＳ６－１５は、上記ステップＳ６－１３とＳ６－１４を繰り返すループ処理を行うか、当該ループ処理を終了するかの選択を行うステップである。
　すなわち、上記ステップＳ６－１４において算出された、リード総量に対してアイソフォームｔに割り当てられたリードの割合「θt」における、前回のループで実行されたステップＳ６－１３とＳ６－１４において得られた「θt'」との間の変化が実質的に認められなくなった場合に、当該「θt」のみならず「φt」も「収束値」として認定され、一連の最適化処理が終了する。この収束がなされるまで、上記ループ処理は実行される。なお、ここで行われる収束のチェックは、「φt」に対して、あるいは、「θt」と「φt」の双方において行っても良いが、いずれの場合でも同様の結果が得られる。
〔２〕　変分ベイズ法によるベイズ推定
　ステップＳ６－２２は、本コンピュータプログラムの実行による最適化処理に際して用いられるパラメータＥ［Ｚnths］（リードｎのアイソフォームｔのアリルｈの位置ｓへの割り当て分のＺnths＝１の期待値）と、θt（リード総量に対してアイソフォームｔに割り当てられたリードの割合）と、φt（アイソフォームｔにおける、父アリルからの発現量の割合）の初期化を行うためのステップである。ここでＮは「リード総数」であり、Ｔは「アイソフォームの種類の総数」である。ｈは、父アリルの場合ｈ＝０であり、母アリルの場合ｈ＝１である。Ｓは「リードｎの開始位置」であり、ｌthはアイソフォームｔのｈアリルの長さであり、Ｌはリード長である。Ｅ［Ｚnths］は、事前の無情報・一様分布を前提として、１／Ｍn（Ｍnは、観測されたリードｎのアライメント先の総数）として初期設定がなされている。
　θtの期待値であるＥθ［θt］は、ディリクレ分布を事前分布として、α* t／Σt'α* t'として与えられる。ここでα* t＝α0＋Ｎ／Ｔである。前述したように、α0はハイパーパラメータであり、例えば、事前知識が何も無いと仮定して、一様分布を与えるα0= 1  が用いられる。
　φtの期待値であるＥφ［φt］は、β分布を事前分布として、β* t,1 ／　β* t,1＋β* t,2として与えられる。ここでβ* t,1＝βt,1 ＋Ｎ／２Ｔ、β* t,2＝βt,2 ＋Ｎ／２Ｔ、である。βt,1 、βt,2 はハイパーパラメータであり、例えば父アリル由来および母アリル由来がどちらも同様に確からしいとして、また事前知識が何も無いと仮定して、一様分布を与えるβt,1 ＝βt,2 ＝１ が用いられる。
　ステップＳ６－２３は、直前に更新された（初回は上記ステップＳ６－２２の初期化値）θtとφtに基づいて、リードｎのアイソフォームｔのアリルｈへの割り当て分のＺnths＝１の期待値：Ｅz［Ｚnths］）を再計算するステップであり、本発明の変分ベイズ法のＥステップに相当するステップである。Ｅz［Ｚnths］の推定値は、前記したように、以下の要領で算出できる。
　ｑ*（θ）及ｑ*（φ）の現在の推定値に基づいて、
Figure JPOXMLDOC01-appb-M000031
 として計算される。
　ここで、シングルエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000032
 である。
（ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（１）（ａ）に、第２項は上記式（１）（ｂ）に、第３項は上記式（１）（ｃ）に、及び、第４項は上記式（２）について、対数計算することにより算出することができる）
　ペアエンドリードの場合は、
Figure JPOXMLDOC01-appb-M000033
 （ｌｏｇρnthsを算出するための「ｈ＝０」と「それ以外」を示す上記２式における、右辺第１項は上記式（３）（ａ）に、第２項は上記式（３）（ｂ）に、第３項は上記式（３）（ｃ）に、第４項は上記式（３）（ｄ）に、及び、第５項は上記式３－１と３－２について、対数計算することにより算出することができる）
　ここで、
Figure JPOXMLDOC01-appb-M000034
 （式中、ψ(・)はディガンマ関数である）
である。
　ステップＳ６－２４は、上記ステップＳ６－２３において算出されたＥz［Ｚnths］の推定値に基づき、Ｅθ［θt］及びＥφ［φt］を計算するステップであり、本発明の変分ベイズ法のＭステップに相当するステップである。
　Ｅθ［θt］は、直前のＳ６－２３で算出されたＥz［Ｚnths］の推定値に基づいて、
Figure JPOXMLDOC01-appb-M000035
 として計算される。
　同様に、Ｅφ［φt］は、直前のＳ６－２３で算出されたＥz［Ｚnths］の推定値に基づいて、
Figure JPOXMLDOC01-appb-M000036
 として計算される。
　ステップＳ６－２５は、上記ステップＳ６－２３とＳ６－２４を繰り返すループ処理を行うか、当該ループ処理を終了するかの選択を行うステップである。
　すなわち、上記ステップＳ６－２４において算出された、リード総量に対してアイソフォームｔに割り当てられたリードの推定値Ｅθ［θt］における、前回のループで実行されたステップＳ６－２３とＳ６－２４において得られたＥθ［θt］との間の変化が実質的に認められなくなった場合に、当該Ｅθ［θt］のみならずＥφ［φt］も「収束値」として認定され、一連の最適化処理が終了する。この収束がなされるまで、上記ループ処理は実行される。なお、ここで行われる収束のチェックは、Ｅφ［φt］に対して、あるいは、Ｅθ［θt］とＥφ［φt］の双方において行っても良いが、いずれの場合でも同様の結果が得られる。
　以上のように、〔１〕ＥＭアルゴリズムによる最尤推定と、〔２〕変分ベイズによるベイズ推定、による最適化工程が実行される。
　これらの工程を経て得られた最適化されたθtとφtを推定値として、これらに基づき、アリル特異的遺伝子発現頻度の決定がなされる（ボックスＢ８）。
　以下に、対象をヒトとして本発明を用いたアリル特異的遺伝子発現について検証を行った実施例を開示する。前述の通りに、ＡＳＥはアリル特異的遺伝子発現を意味する用語である。
［１］　シミュレーションデータ分析
　本発明の推定方法の性能を評価するため、ＮＡ１２８７８の２倍体ゲノム配列［それらは、ｈｇ１９から構築されたものであり、ウェブサイト(http://sv.gersteinlab.org/NA12878_diploid)から公的に入手可能であった］に基づいて、合成ＲＮＡ－Ｓｅｑデータ（３０００万リード、１００ｂｐ×２（ペアエンド）、平均断片サイズ４００ｂｐ、標準偏差４０ｂｐ）を、検体のＲＮＡシークエンスから得られたリード情報として準備した。
　父系及び母系ｃＤＮＡ配列を、ｒＳｅｑ(version 0.2.1)で、ＵＣＳＣゲノム注釈ファイル(refFlat.txt)に基づいて２倍体ゲノム配列から作出し、次に、１０，０００アイソフォームをランダムに選択し、対数正規分布に従うように発現レベルを割り当てた。その後、上記のＲＮＡ－Ｓｅｑデータに基づいて０．１％の置換、欠失、及び挿入エラーを伴うＲＮＡ－Ｓｅｑデータを２セット作出し、一方は、ＡＳＥが存在しない帰無仮説に基づいて作出し、他方は、アイソフォームのいくつかにＡＳＥが存在する対立仮説に基づいて作出した。
　帰無仮説データセットについては、アイソフォームの１００％が、父系及び母系アリルを同程度に発現している（５０：５０の確率）と仮定した。他方、ＡＳＥデータセットについて、アイソフォームの１０％は、父系特異的発現（その中では、父系アリルは８０％の確率で選択されて発現すると設定したのに対し、母系アリルは２０％の確率で選択されて発現すると設定した）を有しており、アイソフォームの１０％は、母系特異的発現（その中では、母系アリルは８０％の確率で選択されて発現すると設定したのに対し、父系アリルは２０％の確率で選択されて発現すると設定した）を有しており、そして残りのアイソフォームはＡＳＥを有していないと仮定した。
　既存アプローチ（Rozowsky et al., AlleleSeq:analysis of allele-specific expression and binding in a network. Molecular Systems Biology,7,522,2011）として、上記リード情報を父系及び母系ハプロタイプの両方にマッピングし、リードのマッピングデータを得た。次に、各アイソフォームについて、ヘテロ接合ＳＮＰｓの数を数えて、父系／母系比を決定した。
　他方、本発明の推定方法では、父母両方のハプロタイプに上記リード情報をマッピングし、Ｂｏｗｔｉｅ２に「－ｋ」オプションを指定して、考えられるアラインメントを全て保持した。その後、ＢＡＭファイルをインプットにし、父系及び母系アリルの間で、さらに、前記図４のフローシートに示した変法ベイズ推定アルゴリズムにより、アイソフォームについてのリードアライメントを最適化した。ハイパーパラメータα0は０．１に設定し、βt,1＝βt,2＝１と設定した。この設定において、データの対数周辺尤度の変分下限を最大にした。
　図５は、帰無仮説及びＡＳＥ仮説についての父系／母系比の真の分布、２倍体ゲノムに対するリードの、上記既存アプローチに基づくソフトウエア（ｂｅｓｔ　ａｌｉｇｎｍｅｎｔ）のアライメントに基づいて予測された分布、及び、本発明の推定方法で予測された分布を示している。１０以上の割り当てられたリードを持つアイソフォームについて、比較検討された。ＡＳＥの有無に係わらず、本発明の推定方法での予測された分布は、真の分布に、特に、父系／母系比が０又は１に近い領域において、より類似していた。一方、上記の既存アプローチでの予測された分布は、それらの極端な領域において「ピーク」を示しているが、それは実際には、真の分布において存在していなかった。本発明の推定方法での好ましい結果は、ベイズ推定において、ハプロタイプ選択についてのβ分布が更新される際、１の事前カウントが父系／母系比を計算するためにアイソフォームの各アリルに加えられることによって平滑化効果が得られることによる（ラプラス・スムージング、又はアッド－ワン・スムージングと呼ばれる）。この特徴は、発現レベルが元々低いアイソフォームについて、あるいは、アリルを区別するのに使用することができるヘテロ接合ＳＮＰｓが多くない時に、特に有益であろう。
Figure JPOXMLDOC01-appb-T000037
［２］　実在のデータ分析
　本発明の推定方法を、リンパ芽球様細胞株ＧＭ１２８７８(Marinov et al., Genome Research, 24, 496-510(2014))から生成したＲＮＡ－Ｓｅｑデータ（１００ｂｐ×２（ペアエンド）の３６５０万リード：ＮＣＢＩ　ＳＲＡ受託番号ＳＲＸ２４５４３４に基づいて公的に入手可能である）をＨａｐＭａｐ　ＮＡ１２８７８個人の２倍体ゲノムにアライメントして得られたマッピング情報に対して当てはめた。この細胞株は、ＨａｐＭａｐ　ＮＡ１２８７８個人から得られ、その２倍体ゲノムは、シミュレーションデータ分析の際と同様に入手し使用された。
　これにより、父系又は母系アリルのいずれかに由来するＡＳＥを示したいくつかの常染色体遺伝子があることを見出された（図６中の上部左）。続く分析で、それらのアイソフォームの父系／母系比が、≧０．７５又は≦０．２５である場合、遺伝子をＡＳＥ遺伝子と見なした。遺伝子のどの機能別カテゴリーがアリル特異的方法で調節されているかを調べるために、ＤＡＶＩＤ（Dennis et al., Genome Biology,4,P3(2013)）を用いて、常染色体の１，２５１個のＡＳＥ遺伝子において多く含まれる機能別カテゴリーを同定した。多く含まれていた用語は、「多型（ｐｏｌｙｍｏｒｐｈｉｓｍ）」、「配列変異（ｓｅｑｕｅｎｃｅ　ｖａｒｉａｎｔ）」、及び「スプライシング変異（ｓｐｌｉｃｉｎｇ　ｖａｒｉａｎｔ）」を含んでおり（表１）、それらは、当該母集団の中のハプロタイプの間のゲノム変異によって説明されるかもしれない。例えば、「多型」注釈は、ヒトの中に少なくとも１つの変異があることを意味し、それは病気に直接関与していない(Boeckmann et al., Nucleic Acids Res,31(1),365-370,(2003))。しかし、Ｇｅｎｅ　Ｏｎｔｏｌｏｇｙ　Ｔｅｒｍｓ(The Gene Ontology Consortium,2000)の中のいずれの機能別カテゴリーも、この分析においてボンフェローニ調整ｐ値０．００１で有意であることが見出されなかった。さらに常染色体ＡＳＥアイソフォームの総存在量を、ＡＳＥを持たない常染色体アイソフォームのそれと比較したところ、前者が後者よりも小さい傾向があった（図７）。これは、ゲノム変異又は他の調節機構による１つのアリルからのより低い発現は、当該細胞株中のもう一方のアリルからの発現によって補われなかったことを示唆している。従って、当該細胞株においてＡＳＥを示す遺伝子は、一般に、ハウスキーピング遺伝子であるとは考えられなかった。
　驚くべきことに、各染色体上の発現したアイソフォームの父系／母系比を調べることによって、ＧＭ１２８７８の非対称な父系アリルのＸ－不活化が観察された（図４中の下部右）。この結果は、ＣＴＣＦ結合(McDaniell et al., Science,328,235-239(2010))及びＣｈＩＰ－ＳｅｑデータからのＲＮＡポリメラーゼII(Reddy et al., Genome Reserch,22,860-869(2012))の占有からＸ－染色体不活化におけるバイアスを示した以前の研究における知見と一致していた。
　本発明の推定方法、さらにそれに基づくコンピュータシステムとコンピュータプログラムは、検体のＲＮＡシークエンスにより得られたリード情報に基づいて、アリル特異的遺伝子発現を的確に推定することが可能である。本技術により、がん等の疾患の早期発見のためのバイオマーカーの作出、創薬、ゲノム変異と疾患関連遺伝子機能との関連解析、が可能となると期待されている。さらに、特に植物及び魚類では、２倍体のみならず、３倍体以上の倍数体が知られている他、各種の品種改良の手段として人為的な倍数体の作出が行われている。これらの倍数体におけるアリル特異的遺伝子発現頻度を推定し、より収量の期待出来る品種へ改良する手がかりを得ることが、産業政策上極めて重要と考えられる。
Ｂ１－１・・・ゲノムＤＮＡのリードデータの存在を示すボックス
Ｂ１－２・・・ヒトゲノム参照配列の存在を示すボックス
Ｓ１・・・ゲノムＤＮＡのマッピングを行うステップ
Ｂ２・・・Ｓ１のマッピングの結果の存在を示すステップ
Ｓ２・・・Ｂ２のマッピングの結果に対して、変異コールとハプロタイプフェージングを行うステップ
Ｂ３・・・Ｓ２の変異コールとフェージングの結果の存在を示すステップ
Ｓ３・・・Ｂ３のデータを基にディプロイドゲノム配列の構築を行うステップ
Ｂ４－１・・・父方のハプロイドゲノム配列の存在を示すボックス
Ｂ４－２・・・母方のハプロイドゲノム配列の存在を示すボックス
Ｓ４・・・Ｂ４の父方と母方のハプロイドゲノム配列をｃＤＮＡ配列へと推定するステップ
Ｂ５－１・・・父方のｃＤＮＡ配列の存在を示すボックス
Ｂ５－２・・・母方のｃＤＮＡ配列の存在を示すボックス
Ｂ６・・・ＲＮＡシークエンスリードデータの存在を示すボックス
Ｓ５・・・ＲＮＡシークエンスデータの父母ｃＤＮＡ配列のマッピングを行うステップ
Ｂ７・・・Ｓ５のマッピングの結果の存在を示すボックス
Ｓ６・・・父母アリル混在データを最適化するステップ
Ｓ６－１１・・・父母アリル混在データの存在を示すステップ
Ｓ６－２１・・・父母アリル混在データの存在を示すステップ
Ｓ６－１２・・・最適化処理を行うためのパラメータの初期化を行うステップ
Ｓ６－２２・・・最適化処理を行うためのパラメータの初期化を行うステップ
Ｓ６－１３・・・最適化処理を行うためのＥステップの機能を実行するためのステップ
Ｓ６－２３・・・最適化処理を行うためのＥステップの機能を実行するためのステップ
Ｓ６－１４・・・最適化処理を行うためのＭステップの機能を実行するためのステップ
Ｓ６－２４・・・最適化処理を行うためのＭステップの機能を実行するためのステップ
Ｓ６－１５・・・最適化処理を行うためのループ・収束機能を実行するためのステップ
Ｓ６－２５・・・最適化処理を行うためのループ・収束機能を実行するためのステップ
Ｂ８・・・アリル特異的遺伝子発現の推定結果の存在を示すステップ

Claims (19)

  1. 　Ｋ倍数体生物（Ｋは２以上の整数）のｃＤＮＡのリード情報が混在したデータに対して、以下のステップが実行される、コンピュータによるアリル特異的遺伝子発現頻度の推定方法。
    （１）　各リードにおける各アイソフォーム及び当該アイソフォームの各アリルに対する期待マッピング数の数値化が行われるステップ；
    （２）　ステップ（１）において数値化された期待マッピング数が各アイソフォーム毎に合算されて合計期待マッピング数が算出され、かつ、当該合計期待マッピング数の各アリル分がそれぞれ算出されるステップ；
    （３）　ステップ（２）において算出された合計期待マッピング数は、（ｉ）各アイソフォームに対する合計期待マッピング数の和で除されてアイソフォーム量比が算出され、かつ、(ii)各アイソフォームの各アリルにおいて、当該アイソフォームの合計期待マッピング数の当該アリル分が当該アイソフォームの合計期待マッピング数で除されてアリル好選度が算出されるステップ；
    （４）　ステップ（３）において算出された、新たなアイソフォーム量比とアリル好選度に対して再びステップ（１）が実行されて、新たな期待マッピング数が算出されるステップ；
    （５）　ステップ（４）において得られた新たな期待マッピング数に対して、再びステップ（２）又は（３）が実行され、アイソフォーム量比とアリル好選度が新たに算出されるステップ；
    （６）　ステップ（４）と（５）が、（ｉ）ステップ（４）において算出される期待マッピング数と、前回のステップ（４）で算出された期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、(ii)ステップ（５）において算出されるアイソフォーム量比及び／又はアリル好選度と、前回のステップ（５）において算出されたアイソフォーム量比及び／又はアリル好選度の間における差が全てのアリルにおいて認められなくなるまで、繰り返し実行され、収束した期待マッピング数、あるいは、収束したアイソフォーム量比及び／又はアリル好選度が、最適化されたデータとして認定されるステップ。
  2. 　Ｋ倍数体生物（Ｋは２以上の整数）のｃＤＮＡのリード情報が混在したデータにおけるリード全体の塩基配列を観測データＲとして、各リードにおける各アリル由来のアイソフォームに対する期待マッピング数を求めるステップ、目的パラメータであるアイソフォーム量θt（ｔ＝１，…，Ｔ：Ｔはアイソフォームの種類数）及び当該アイソフォームｔのアリル選好度φtk（φtk：ｋ＝１，…，Ｋ）の推定値を求めるステップ、並びに、当該推定値に基づいて、アリル特異的遺伝子発現頻度の決定を行うステップ、を含むアリル特異的遺伝子発現頻度のコンピュータによる推定方法において、
    　上記目的パラメータθtとφtk、及び、観測データＲを媒介する３種の潜在変数として、（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔn、（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈn、及び、（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳnについて、
    　リードｎの塩基配列を観測データＲnとして、観測データＲnからの目的パラメータθtとφtkの推測工程において、少なくともこれらの３種の潜在変数を、観測データＲnが依存するように繰り入れて、当該推定値を算出することを特徴とする、アリル特異的遺伝子発現頻度の推定方法。
  3. 　Ｋ倍数体生物（Ｋは２以上の整数）のｃＤＮＡのリード情報が混在したデータは、検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリルのｃＤＮＡの塩基配列に対してマッピングを行うことにより得られることを特徴とする、請求項１又は２に記載のアリル特異的遺伝子発現頻度の推定方法。
  4. 　Ｋ倍数体生物（Ｋは２以上の整数）の由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列は、検体のゲノムＤＮＡのリードシークエンスから推定されて得られることを特徴とする、請求項３に記載のアリル特異的遺伝子発現頻度の推定方法。
  5. 　最尤推定法、又は、ベイズ推定法に基づくステップの実行により、全てのリードにおける個々のアリル由来アイソフォームに対する期待マッピング数、目的パラメータであるアイソフォーム量、及び、アリル選好度、の推定値を算出することを特徴とする、請求項１～４のいずれか１項に記載のアリル特異的遺伝子発現頻度の推定方法。
  6. 　３種の潜在変数Ｔn、Ｈn及びＳnが要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）を潜在変数として用いることを特徴とする、請求項２～５のいずれか１項に記載のアリル特異的遺伝子発現頻度の推定方法。
  7. 　下記（１）～（６）のステップを行うことを特徴とする、請求項６に記載のアリル特異的遺伝子発現頻度の推定方法。
    　（１）所与されたθtの初期値、及び、φtkの初期値に基づいて、Ｚnths＝１の事後確率の第１の更新値を算出し、さらに、当該Ｚnths＝１の事後確率の第１の更新値に基づいてθt、及び、φtkの最尤推定値の第１の更新値を算出するステップ、あるいは、（２）所与されたＺnths＝１の事後確率の初期値に基づいて、θt、及び、φtkの最尤推定値の第１の更新値を算出するステップ、
    　（３）　直前のステップ（４）（ただし、初回はステップ（１）又はステップ（２）である）により算出されたθt、及び、φtkの最尤推定値の更新値に基づいて、新たにＺnths＝１の事後確率の更新値を算出するステップ、
    　（４）　ステップ（３）において算出されたＺnths＝１の事後確率の更新値に基づいて、θt、及び、φtkの最尤推定値の更新値を新たに算出するステップ、
    　（５）　ステップ（４）において算出された最尤推定値の更新値の対数尤度を算出するステップ、
    　（６）　ステップ（５）において算出された対数尤度の収束性を評価するステップ、もしくはステップ（３）で算出されたＺnths＝１の事後確率の更新値の収束性を評価するステップ、もしくはステップ（４）で算出されたθt、及び、φtkの最尤推定値の更新値の収束性を評価するステップであって、収束が認められれば、当該収束値をθtとφtkの最終推定値として決定し、収束が認められなければ、ステップ（３）、（４）、（５）、及び、（６）の繰り返しを決定する。
  8. 　下記（１）～（５）のステップを行うことを特徴とする、請求項６に記載のアリル特異的遺伝子発現頻度の推定方法。
    　（１）　所与されたアイソフォームｔの存在量についての予備知識の分布を示すハイパーパラメータαtの初期値に基づくθtの事後分布の更新値、及び、所与された個々のアリルに割り当てられるリード比についての予備知識の分布を示すハイパーパラメータβt,γ（γは、 １，…，Ｋの各整数であり、βt,γはＫ個存在する。）の初期値に基づくφtkの事後分布の更新値、に基づいてＺnths＝１の事後分布を算出し、さらに、当該Ｚnths＝１の事後分布に基づいてθt、及び、φtkの第１の更新事後分布の更新値を算出するステップ、あるいは、（２）所与されたＺnths＝１の初期分布に基づいてθt、及び、φtkの第１の事後分布の更新値を算出するステップ、
    　（３）　直前のステップ（４）（ただし、初回はステップ（１）又はステップ（２）である）により算出されたθt、及び、φtkの事後分布の更新値に基づいて、新たにＺnths＝１の事後分布を算出するステップ、
    　（４）　ステップ（３）において算出された　Ｚnths＝１の事後分布を基にして、θtとφtkの事後分布を新たに算出して更新するステップ、
    　（５）　ステップ（４）において得られたθtとφtkの事後分布の期待値の収束性を評価するステップであって、当該期待値における収束が認められれば、当該収束期待値をθtとφtkの推定値として決定し、収束が認められなければ、ステップ（２）、（３）、及び、（４）の繰り返しを決定する。
  9. 　Ｋ倍数体生物のアリル数を示すＫが２である、請求項１～８のいずれか１項に記載の推定方法。
  10. 　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え；
    （１）　当該記録部には、当該被験生物の全てのアリル由来のｃＤＮＡのリード情報が、リードの配列及びリードのマッピング先の観測データとして記録されており、
    （２）　当該演算処理部では読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２が択一的に実行され、
    　（２）－１：アイソフォーム量に関する変数θtの初期値の算出処理、及び、アリル好選度に関する変数φtkの初期値の算出処理、
    　（２）－２：上記変数θtとφtk、及び、観測データである当該被験生物の全てのアリル由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数としての下記（ａ）、（ｂ）及び（ｃ）：
    　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔn、
    　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈn、
    　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳn、
    が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の初期値の算出処理、
    （３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θtとφtkに基づき、当該指標変数Ｚnthsの算出処理がなされ、
    （４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsに基づいて変数θtとφtkの最尤推定値の第１の更新値が算出され、
    （５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの最尤推定値の第１の更新値に基づいて上記ステップ（３）とステップ（４）が再度実行され、さらに、変数θtとφtkの第２の更新値が算出されるループ処理が、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束した変数θtとφtkが最適化されたθtとφtkのデータとして、上記記録部に記録がなされる；
    　ことを特徴とする、個々のアリル由来のｃＤＮＡのリード情報を最適化するコンピュータシステム。
  11. 　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え；
    （１）　当該記録部には、被験生物の各アリル由来のｃＤＮＡのリード情報が、リードの配列及びリードのマッピング先の観測データとして記録されており、
    （２）　当該演算処理部では読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２が択一的に実行され、
    　（２）－１：アイソフォームｔの存在量についての予備知識の分布を示すハイパーパラメータαtの初期値に基づくθtの事後分布の更新値の算出処理、及び、個々のアリルに割り当てられるリード比についての予備知識の分布を示すハイパーパラメータβt,γ（γは、１，…，Ｋの各整数であり、βt,γはＫ個存在する。）の初期値に基づくφtkの事後分布の更新値の算出処理、
    　（２）－２：上記θtとφtkの分布、及び、観測データである被験者の父方由来と母方由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数分布としての下記（ａ）、（ｂ）及び（ｃ）：
    　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔnの分布、
    　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈnの分布、
    　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳnの分布、が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の分布の初期分布の算出処理、
    （３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θとφの分布に基づき、当該指標変数Ｚnthsの事後分布の算出処理がなされ、
    （４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsの事後分布の更新値に基づいて変数θtとφtkの第１更新事後分布の更新値が算出され、
    （５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの第１の更新事後分布に基づいて上記ステップ（３）とステップ（４）が再度実行され、さらに、変数θtとφtkの第２の更新事後分布が算出されるループ処理が、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束したθtとφtkの期待値が最適化されたθtとφtkのデータとして、上記記録部に記録がなされる；
    　ことを特徴とする、個々のアリル由来のｃＤＮＡのリード情報を最適化するコンピュータシステム。
  12. 　上記演算部において、被験生物の検体のＲＮＡ－シークエンスにより得られたリード情報の、当該検体の個々のアリル由来のｃＤＮＡの塩基配列のマッピング処理が同時に行われて、当該マッピング情報が、全てのアリル由来のｃＤＮＡのリード情報として上記記録部に記録されて以降の処理が実行されることを特徴とする、請求項１０又は１１に記載のコンピュータシステム。
  13. 　前記個々のアリル由来のｃＤＮＡの塩基配列は、検体のゲノムＤＮＡのリードシークエンスから推定されて得られることを特徴とする、請求項１２に記載のコンピュータシステム。
  14. 　Ｋ倍数体生物のアリル数を示すＫが２である、請求項９～１３のいずれか１項に記載のコンピュータシステム。
  15. 　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに、
    （１）　記録部に記録されている当該被験生物の各アリル由来のｃＤＮＡのリード情報を、リードの配列及びリードのマッピング先の観測データとして読み出す、第１の機能；
    （２）　演算処理部において読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２を択一的に実行する、第２の機能；
    　（２）－１：アイソフォーム量に関する変数θtの初期値の算出処理、及び、アリル好選度に関する変数φtkの初期値の算出処理、
    　（２）－２：上記変数θtとφtk、及び、観測データである当該被験生物の各アリル由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数としての下記（ａ）、（ｂ）及び（ｃ）：
    　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔn、
    　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈn、
    　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳn、
    が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の事後確率（リードｎの期待マッピング数）の初期値の算出処理、
    （３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θtとφtkに基づき、当該指標変数Ｚnthsの事後確率（リードｎの期待マッピング数）の算出処理を行う、第３の機能；
    （４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsの事後確率に基づいて変数θtとφtkの最尤推定値の第１の更新値の算出処理を行う、第４の機能；
    （５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの最尤推定値の第１の更新値に基づいて上記ステップ（３）とステップ（４）を再度実行し、さらに、変数θtとφtkの第２の更新値が算出されるループ処理を、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束した変数θtとφtkを最適化されたθtとφtkのデータとして、上記記録部への記録を行う、第５の機能；
    　を実現させるアルゴリズムが含まれていることを特徴とする、各アリル由来のｃＤＮＡのリード情報を最適化するコンピュータプログラム。
  16. 　Ｋ倍数体生物（Ｋは２以上の整数）の検体のＲＮＡ－シークエンスにより得られたリードを、当該検体の各アリル由来のｃＤＮＡの塩基配列に対して同時にマッピングを行うことにより得られる各アリル由来のｃＤＮＡのリード情報が混在したデータについて、全てのリードにおける各アリル由来のアイソフォームに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに、
    （１）　記録部に記録されている当該被験生物の各アリル由来のｃＤＮＡのリード情報を、リードの配列及びリードのマッピング先の観測データとして読み出す、第１の機能；
    （２）　演算処理部において読み出された前記観測データに基づき、下記の初期化ステップ（２）－１又は（２）－２を択一的に実行する、第２の機能；
    　（２）－１：アイソフォームｔの存在量についての予備知識の分布を示すハイパーパラメータαtの初期値に基づくθtの分布の期待値の算出処理、及び、各アリルに割り当てられるリード比についての予備知識の分布を示すハイパーパラメータβt,γ（γは、１，…，Ｋの各整数であり、βt,γはＫ個存在する。）の初期値に基づくφtkの分布の期待値の算出処理、
    　（２）－２：上記θtとφtkの分布、及び、観測データである当該被験生物の全てのアリル由来のｃＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する３種の潜在変数分布としての下記（ａ）、（ｂ）及び（ｃ）：
    　（ａ）リードｎのアイソフォーム選択に関する、θtに依存する変数Ｔnの分布、
    　（ｂ）リードｎのハプロタイプ選択に関する、φtkとＴnに依存する変数Ｈnの分布、
    　（ｃ）リードｎの開始位置に関する、Ｔn及びＨnに依存するＳnの分布、
    が要約された、指標変数Ｚnths（Ｚnthsは、（Ｔn，Ｈn，Ｓn）＝（ｔ，ｈ，ｓ）の場合１であり、それ以外は０である。）の分布（リードｎの期待マッピング数）の初期分布の算出処理、
    （３）　当該演算処理部において、上記ステップ（２）－１で算出された変数θtとφtkの分布に基づき、当該指標変数Ｚnthsの分布（リードｎの期待マッピング数）の算出処理を行う、第３の機能；
    （４）　当該演算処理部において、上記ステップ（２）－２、又は、ステップ（３）で算出された当該指標変数Ｚnthsの分布に基づいて変数θtとφtkの第１更新事後分布の算出処理を行う、第４の機能；
    （５）　当該演算処理部において、上記ステップ（４）で算出された変数θtとφtkの第１の更新事後分布に基づいて上記ステップ（３）とステップ（４）を再度実行し、さらに、変数θtとφtkの第２の更新事後分布が算出されるループ処理を、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束したθtとφtkの期待値を最適化されたθtとφtkのデータとして、上記記録部への記録を行う、第５の機能；
    　を実現させるアルゴリズムが含まれていることを特徴とする、個々のアリル由来のｃＤＮＡのリード情報を最適化するコンピュータプログラム。
  17. 　コンピュータの演算部において、検体のＲＮＡ－シークエンスにより得られたリード情報の、当該検体の各アリル由来のｃＤＮＡの塩基配列のマッピング処理を同時に行い、当該マッピング情報を、各アリル由来のｃＤＮＡのリード情報として上記記録部に記録する機能が、上記第１の機能の上流に付加されている特徴とする、請求項１５又は１６に記載のコンピュータプログラム。
  18. 　Ｋ倍数体生物のアリル数を示すＫが２である、請求項１５～１７のいずれか１項に記載のコンピュータプログラム。
  19. 　請求項１５～１８のいずれか１項に記載のコンピュータプログラムが記録されていることを特徴とする、コンピュータにおいて読み取り可能な記録媒体。
PCT/JP2016/088640 2015-12-28 2016-12-26 アリル特異的遺伝子発現頻度の推定方法、推定用コンピュータシステム及び推定用プログラム WO2017115741A1 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2015256860 2015-12-28
JP2015-256860 2015-12-28

Publications (1)

Publication Number Publication Date
WO2017115741A1 true WO2017115741A1 (ja) 2017-07-06

Family

ID=59224728

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2016/088640 WO2017115741A1 (ja) 2015-12-28 2016-12-26 アリル特異的遺伝子発現頻度の推定方法、推定用コンピュータシステム及び推定用プログラム

Country Status (1)

Country Link
WO (1) WO2017115741A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020227952A1 (zh) * 2019-05-15 2020-11-19 深圳华大基因股份有限公司 基于测序数据的碱基突变检测方法、装置及存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
NARIAI NAOKI: "TIGAR: transcript isoform abundance estimation method with gapped alignment of RNA-Seq data by variational Bayesian inference", BIOINFORMATICS, vol. 29, no. 18, 2013, pages 2292 - 2299, XP055396200, Retrieved from the Internet <URL:https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btt381> [retrieved on 20170130] *
NARIAI NAOKI: "TIGAR2: sensitive and accurate estimation of transcript isoform expression with longer RNA-Seq reads", BMC GENOMICS, vol. 15, no. 10, 2014, pages 55, XP021204616, Retrieved from the Internet <URL:http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-S10-S5> [retrieved on 20170130] *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020227952A1 (zh) * 2019-05-15 2020-11-19 深圳华大基因股份有限公司 基于测序数据的碱基突变检测方法、装置及存储介质

Similar Documents

Publication Publication Date Title
JP7368483B2 (ja) 相同組換え欠損を推定するための統合された機械学習フレームワーク
Naslavsky et al. Whole-genome sequencing of 1,171 elderly admixed individuals from Brazil
Oesper et al. Quantifying tumor heterogeneity in whole-genome and whole-exome sequencing data
Juric et al. The strength of selection against Neanderthal introgression
Zhu et al. Deconvolution of multiple infections in Plasmodium falciparum from high throughput sequencing data
Henn et al. Cryptic distant relatives are common in both isolated and cosmopolitan genetic samples
KR102514024B1 (ko) 유전적 변이의 비침습 평가를 위한 방법 및 프로세스
Wilson Sayres et al. Natural selection reduced diversity on human Y chromosomes
Prüfer snpAD: an ancient DNA genotype caller
US10777302B2 (en) Identifying variants of interest by imputation
Xie et al. H-PoP and H-PoPG: heuristic partitioning algorithms for single individual haplotyping of polyploids
US20190065670A1 (en) Predicting disease burden from genome variants
US20220130488A1 (en) Methods for detecting copy-number variations in next-generation sequencing
Tatsumoto et al. Direct estimation of de novo mutation rates in a chimpanzee parent-offspring trio by ultra-deep whole genome sequencing
Moltke et al. A method for detecting IBD regions simultaneously in multiple individuals—with applications to disease genetics
Rödelsperger et al. Identity-by-descent filtering of exome sequence data for disease–gene identification in autosomal recessive disorders
JP2006519440A (ja) 疾患の増大リスクの統計学的同定法
Kristmundsdóttir et al. popSTR: population-scale detection of STR variants
Cheung et al. A statistical framework to guide sequencing choices in pedigrees
Wu et al. Estimating error models for whole genome sequencing using mixtures of Dirichlet-multinomial distributions
Woodhams et al. Simulating and summarizing sources of gene tree incongruence
EP3239875B1 (en) Method for determining genotype of particular gene locus group or individual gene locus, determination computer system and determination program
WO2017115741A1 (ja) アリル特異的遺伝子発現頻度の推定方法、推定用コンピュータシステム及び推定用プログラム
Zou et al. Inferring parental genomic ancestries using pooled semi-Markov processes
Markus et al. Integration of SNP genotyping confidence scores in IBD inference

Legal Events

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

Ref document number: 16881716

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 16881716

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP