System and method for cleaning noisy genetic data and determining chromosome copy number

[0001] This application is a continuation of U.S. application Ser. No. 12/076,348, filed Mar. 17, 2008, now U.S. Pat. No. 8,515,679; which is a continuation-in-part of U.S. application Ser. No. 11/496,982, filed Jul. 31, 2006, now abandoned; a continuation-in-part of U.S. application Ser. No. 11/603,406, filed Nov. 22, 2006, now U.S. Pat. No. 8,532,930; and a continuation-in-part of U.S. application Ser. No. 11/634,550, filed Dec. 6, 2006, now abandoned; which claims the benefit of U.S. 60/846,610, filed Sep. 22, 2006; U.S. 60/817,741, filed Jun. 30, 2006; U.S. 60/789,506, filed Apr. 4, 2006; U.S. 60/774,976, filed Feb. 21, 2006; U.S. 60/754,396, filed Dec. 29, 2005; U.S. 60/742,305, filed Dec. 6, 2005; the disclosures thereof are incorporated by reference herein in their entirety. U.S. patent application Ser. No. 12/076,348 also claims the benefit of priority to U.S. 60/918,292, filed Mar. 16, 2007; U.S. 60/926,198, filed Apr. 25, 2007; U.S. 60/932,456, filed May 31, 2007; U.S. 60/934,440, filed Jun. 11, 2007; U.S. 61/003,101, filed Nov. 12, 2007; and U.S. 61/008,637, filed Dec. 21, 2007; the disclosures thereof are incorporated by reference herein in their entirety. U.S. application Ser. No. 11/496,982 claims the benefit of priority to U.S. 60/703,415, filed Jul. 29, 2005; and U.S. application Ser. No. 11/603,406 claims the benefit of priority to U.S. 60/739,882, filed Nov. 26, 2005; the disclosures thereof are incorporated by reference herein in their entirety.

[0002] This invention was made with government support under Grant No. R44HD054958-02A2 awarded by the National Institutes of Health. The government has certain rights in the invention.

[0003] The invention relates generally to the field of acquiring, manipulating and using genetic data for medically predictive purposes, and specifically to a system in which imperfectly measured genetic data of a target individual are made more accurate by using known genetic data of genetically related individuals, thereby allowing more effective identification of genetic variations, specifically aneuploidy and disease linked genes, that could result in various phenotypic outcomes.

[0004] In 2006, across the globe, roughly 800,000 in vitro fertilization (IVF) cycles were run. Of the roughly 150,000 cycles run in the US, about 10,000 involved pre-implantation genetic diagnosis (PGD). Current PGD techniques are unregulated, expensive and highly unreliable: error rates for screening disease-linked loci or aneuploidy are on the order of 10%, each screening test costs roughly $5,000, and a couple is forced to choose between testing aneuploidy, which afflicts roughly 50% of IVF embryos, or screening for disease-linked loci on the single cell. There is a great need for an affordable technology that can reliably determine genetic data from a single cell in order to screen in parallel for aneuploidy, monogenic diseases such as Cystic Fibrosis, and susceptibility to complex disease phenotypes for which the multiple genetic markers are known through whole-genome association studies.

[0005] Most PGD today focuses on high-level chromosomal abnormalities such as aneuploidy and balanced translocations with the primary outcomes being successful implantation and a take-home baby. The other main focus of PGD is for genetic disease screening, with the primary outcome being a healthy baby not afflicted with a genetically heritable disease for which one or both parents are carriers. In both cases, the likelihood of the desired outcome is enhanced by excluding genetically suboptimal embryos from transfer and implantation in the mother.

[0006] The process of PGD during IVF currently involves extracting a single cell from the roughly eight cells of an early-stage embryo for analysis. Isolation of single cells from human embryos, while highly technical, is now routine in IVF clinics. Both polar bodies and blastomeres have been isolated with success. The most common technique is to remove single blastomeres from day 3 embryos (6 or 8 cell stage). Embryos are transferred to a special cell culture medium (standard culture medium lacking calcium and magnesium), and a hole is introduced into the zona pellucida using an acidic solution, laser, or mechanical techniques. The technician then uses a biopsy pipette to remove a single blastomere with a visible nucleus. Features of the DNA of the single (or occasionally multiple) blastomere are measured using a variety of techniques. Since only a single copy of the DNA is available from one cell, direct measurements of the DNA are highly error-prone, or noisy. There is a great need for a technique that can correct, or make more accurate, these noisy genetic measurements.

[0007] Normal humans have two sets of 23 chromosomes in every diploid cell, with one copy coming from each parent. Aneuploidy, the state of a cell with extra or missing chromosome(s), and uniparental disomy, the state of a cell with two of a given chromosome both of which originate from one parent, are believed to be responsible for a large percentage of failed implantations and miscarriages, and some genetic diseases. When only certain cells in an individual are aneuploid, the individual is said to exhibit mosaicism. Detection of chromosomal abnormalities can identify individuals or embryos with conditions such as Down syndrome, Klinefelter's syndrome, and Turner syndrome, among others, in addition to increasing the chances of a successful pregnancy. Testing for chromosomal abnormalities is especially important as the age of a potential mother increases: between the ages of 35 and 40 it is estimated that between 40% and 50% of the embryos are abnormal, and above the age of 40, more than half of the embryos are like to be abnormal. The main cause of aneuploidy is nondisjunction during meiosis. Maternal nondisjunction constitutes 88% of all nondisjunction of which 65% occurs in meiosis 1 and 23% in meiosis II. Common types of human aneuploidy include trisomy from meiosis 1 nondisjunction, monosomy, and uniparental disomy. In a particular type of trisomy that arises in meiosis II nondisjunction, or M2 trisomy, an extra chromosome is identical to one of the two normal chromosomes. M2 trisomy is particularly difficult to detect. There is a great need for a better method that can detect for many or all types of aneuploidy at most or all of the chromosomes efficiently and with high accuracy.

[0008] Karyotyping, the traditional method used for the prediction of aneuploidy and mosaicism is giving way to other more high-throughput, more cost effective methods such as Flow Cytometry (FC) and fluorescent in situ hybridization (FISH). Currently, the vast majority of prenatal diagnoses use FISH, which can determine large chromosomal aberrations and PCR/electrophoresis, and which can determine a handful of SNPs or other allele calls. One advantage of FISH is that it is less expensive than karyotyping, but the technique is complex and expensive enough that generally a small selection of chromosomes are tested (usually chromosomes 13, 18, 21, X, Y; also sometimes 8, 9, 15, 16, 17, 22); in addition, FISH has a low level of specificity. Roughly seventy-five percent of PGD today measures high-level chromosomal abnormalities such as aneuploidy using FISH with error rates on the order of 10-15%. There is a great demand for an aneuploidy screening method that has a higher throughput, lower cost, and greater accuracy.

[0009] The number of known disease associated genetic alleles is currently at 389 according to OMIM and steadily climbing. Consequently, it is becoming increasingly relevant to analyze multiple positions on the embryonic DNA, or loci, that are associated with particular phenotypes. A clear advantage of pre-implantation genetic diagnosis over prenatal diagnosis is that it avoids some of the ethical issues regarding possible choices of action once undesirable phenotypes have been detected. A need exists for a method for more extensive genotyping of embryos at the pre-implantation stage.

[0010] There are a number of advanced technologies that enable the diagnosis of genetic aberrations at one or a few loci at the single-cell level. These include interphase chromosome conversion, comparative genomic hybridization, fluorescent PCR, mini-sequencing and whole genome amplification. The reliability of the data generated by all of these techniques relies on the quality of the DNA preparation. Better methods for the preparation of single-cell DNA for amplification and PGD are therefore needed and are under study. All genotyping techniques, when used on single cells, small numbers of cells, or fragments of DNA, suffer from integrity issues, most notably allele drop out (ADO). This is exacerbated in the context of in-vitro fertilization since the efficiency of the hybridization reaction is low, and the technique must operate quickly in order to genotype the embryo within the time period of maximal embryo viability. There exists a great need for a method that alleviates the problem of a high ADO rate when measuring genetic data from one or a small number of cells, especially when time constraints exist.

[0011] Listed here is a set of prior art which is related to the field of the current invention. None of this prior art contains or in any way refers to the novel elements of the current invention. In U.S. Pat. No. 6,489,135 Parrott et al. provide methods for determining various biological characteristics of in vitro fertilized embryos, including overall embryo health, implantability, and increased likelihood of developing successfully to term by analyzing media specimens of in vitro fertilization cultures for levels of bioactive lipids in order to determine these characteristics. In US Patent Application 20040033596 Threadgill et al. describe a method for preparing homozygous cellular libraries useful for in vitro phenotyping and gene mapping involving site-specific mitotic recombination in a plurality of isolated parent cells. In U.S. Pat. No. 5,635,366 Cooke et al. provide a method for predicting the outcome of IVF by determining the level of 11β-hydroxysteroid dehydrogenase (11β-HSD) in a biological sample from a female patient. In U.S. Pat. No. 7,058,517 Denton et al. describe a method wherein an individual's haplotypes are compared to a known database of haplotypes in the general population to predict clinical response to a treatment. In U.S. Pat. No. 7,035,739 Schadt at al. describe a method is described wherein a genetic marker map is constructed and the individual genes and traits are analyzed to give a gene-trait locus data, which are then clustered as a way to identify genetically interacting pathways, which are validated using multivariate analysis. In US Patent Application US 2004/0137470 A1, Dhallan et al. describe using primers especially selected so as to improve the amplification rate, and detection of, a large number of pertinent disease related loci, and a method of more efficiently quantitating the absence, presence and/or amount of each of those genes. In World Patent Application WO 03/031646, Findlay et al. describe a method to use an improved selection of genetic markers such that amplification of the limited amount of genetic material will give more uniformly amplified material, and it can be genotyped with higher fidelity.

[0012] Current methods of prenatal diagnosis can alert physicians and parents to abnormalities in growing fetuses. Without prenatal diagnosis, one in 50 babies is born with serious physical or mental handicap, and as many as one in 30 will have some form of congenital malformation. Unfortunately, standard methods require invasive testing and carry a roughly 1 percent risk of miscarriage. These methods include amniocentesis, chorion villus biopsy and fetal blood sampling. Of these, amniocentesis is the most common procedure; in 2003, it was performed in approximately 3% of all pregnancies, though its frequency of use has been decreasing over the past decade and a half. A major drawback of prenatal diagnosis is that given the limited courses of action once an abnormality has been detected, it is only valuable and ethical to test for very serious defects. As result, prenatal diagnosis is typically only attempted in cases of high-risk pregnancies, where the elevated chance of a defect combined with the seriousness of the potential abnormality outweighs the risks. A need exists for a method of prenatal diagnosis that mitigates these risks. It has recently been discovered that cell-free fetal DNA and intact fetal cells can enter maternal blood circulation. Consequently, analysis of these cells can allow early Non-Invasive Prenatal Genetic Diagnosis (NIPGD). A key challenge in using NIPGD is the task of identifying and extracting fetal cells or nucleic acids from the mother's blood. The fetal cell concentration in maternal blood depends on the stage of pregnancy and the condition of the fetus, but estimates range from one to forty fetal cells in every milliliter of maternal blood, or less than one fetal cell per 100,000 maternal nucleated cells. Current techniques are able to isolate small quantities of fetal cells from the mother's blood, although it is very difficult to enrich the fetal cells to purity in any quantity. The most effective technique in this context involves the use of monoclonal antibodies, but other techniques used to isolate fetal cells include density centrifugation, selective lysis of adult erythrocytes, and FACS. Fetal DNA isolation has been demonstrated using PCR amplification using primers with fetal-specific DNA sequences. Since only tens of molecules of each embryonic SNP are available through these techniques, the genotyping of the fetal tissue with high fidelity is not currently possible.

[0013] Much research has been done towards the use of pre-implantation genetic diagnosis (PGD) as an alternative to classical prenatal diagnosis of inherited disease. Most PGD today focuses on high-level chromosomal abnormalities such as aneuploidy and balanced translocations with the primary outcomes being successful implantation and a take-home baby. A need exists for a method for more extensive genotyping of embryos at the pre-implantation stage. The number of known disease associated genetic alleles is currently at 389 according to OMIM and steadily climbing. Consequently, it is becoming increasingly relevant to analyze multiple embryonic SNPs that are associated with disease phenotypes. A clear advantage of pre-implantation genetic diagnosis over prenatal diagnosis is that it avoids some of the ethical issues regarding possible choices of action once undesirable phenotypes have been detected.

[0014] Many techniques exist for isolating single cells. The FACS machine has a variety of applications; one important application is to discriminate between cells based on size, shape and overall DNA content. The FACS machine can be set to sort single cells into any desired container. Many different groups have used single cell DNA analysis for a number of applications, including prenatal genetic diagnosis, recombination studies, and analysis of chromosomal imbalances. Single-sperm genotyping has been used previously for forensic analysis of sperm samples (to decrease problems arising from mixed samples) and for single-cell recombination studies.

[0015] Isolation of single cells from human embryos, while highly technical, is now routine in in vitro fertilization clinics. To date, the vast majority of prenatal diagnoses have used fluorescent in situ hybridization (FISH), which can determine large chromosomal aberrations (such as Down syndrome, or trisomy 21) and PCR/electrophoresis, which can determine a handful of SNPs or other allele calls. Both polar bodies and blastomeres have been isolated with success. It is critical to isolate single blastomeres without compromising embryonic integrity. The most common technique is to remove single blastomeres from day 3 embryos (6 or 8 cell stage). Embryos are transferred to a special cell culture medium (standard culture medium lacking calcium and magnesium), and a hole is introduced into the zona pellucida using an acidic solution, laser, or mechanical drilling. The technician then uses a biopsy pipette to remove a single visible nucleus. Clinical studies have demonstrated that this process does not decrease implantation success, since at this stage embryonic cells are undifferentiated.

[0016] There are three major methods available for whole genome amplification (WGA): ligation-mediated PCR (LM-PCR), degenerate oligonucleotide primer PCR (DOP-PCR), and multiple displacement amplification (MDA). In LM-PCR, short DNA sequences called adapters are ligated to blunt ends of DNA. These adapters contain universal amplification sequences, which are used to amplify the DNA by PCR. In DOP-PCR, random primers that also contain universal amplification sequences are used in a first round of annealing and PCR. Then, a second round of PCR is used to amplify the sequences further with the universal primer sequences. Finally, MDA uses the phi-29 polymerase, which is a highly processive and non-specific enzyme that replicates DNA and has been used for single-cell analysis. Of the three methods, DOP-PCR reliably produces large quantities of DNA from small quantities of DNA, including single copies of chromosomes. On the other hand, MDA is the fastest method, producing hundred-fold amplification of DNA in a few hours. The major limitations to amplification material from a single cells are (1) necessity of using extremely dilute DNA concentrations or extremely small volume of reaction mixture, and (2) difficulty of reliably dissociating DNA from proteins across the whole genome. Regardless, single-cell whole genome amplification has been used successfully for a variety of applications for a number of years.

[0017] There are numerous difficulties in using DNA amplification in these contexts. Amplification of single-cell DNA (or DNA from a small number of cells, or from smaller amounts of DNA) by PCR can fail completely, as reported in 5-10% of the cases. This is often due to contamination of the DNA, the loss of the cell, its DNA, or accessibility of the DNA during the PCR reaction. Other sources of error that may arise in measuring the embryonic DNA by amplification and microarray analysis include transcription errors introduced by the DNA polymerase where a particular nucleotide is incorrectly copied during PCR, and microarray reading errors due to imperfect hybridization on the array. The biggest problem, however, remains allele drop-out (ADO) defined as the failure to amplify one of the two alleles in a heterozygous cell. ADO can affect up to more than 40% of amplifications and has already caused PGD misdiagnoses. ADO becomes a health issue especially in the case of a dominant disease, where the failure to amplify can lead to implantation of an affected embryo. The need for more than one set of primers per each marker (in heterozygotes) complicate the PCR process. Therefore, more reliable PCR assays are being developed based on understanding the ADO origin. Reaction conditions for single-cell amplifications are under study. The amplicon size, the amount of DNA degradation, freezing and thawing, and the PCR program and conditions can each influence the rate of ADO.

[0018] All those techniques, however, depend on the minute DNA amount available for amplification in the single cell. This process is often accompanied by contamination. Proper sterile conditions and microsatellite sizing can exclude the chance of contaminant DNA as microsatellite analysis detected only in parental alleles exclude contamination. Studies to reliably transfer molecular diagnostic protocols to the single-cell level have been recently pursued using first-round multiplex PCR of microsatellite markers, followed by real-time PCR and microsatellite sizing to exclude chance contamination. Multiplex PCR allows for the amplification of multiple fragments in a single reaction, a crucial requirement in the single-cell DNA analysis. Although conventional PCR was the first method used in PGD, fluorescence in situ hybridization (FISH) is now common. It is a delicate visual assay that allows the detection of nucleic acid within undisturbed cellular and tissue architecture. It relies firstly on the fixation of the cells to be analyzed. Consequently, optimization of the fixation and storage condition of the sample is needed, especially for single-cell suspensions.

[0019] Advanced technologies that enable the diagnosis of a number of diseases at the single-cell level include interphase chromosome conversion, comparative genomic hybridization (CGH), fluorescent PCR, and whole genome amplification. The reliability of the data generated by all of these techniques rely on the quality of the DNA preparation. PGD is also costly, consequently there is a need for less expensive approaches, such as mini-sequencing. Unlike most mutation-detection techniques, mini-sequencing permits analysis of very small DNA fragments with low ADO rate. Better methods for the preparation of single-cell DNA for amplification and PGD are therefore needed and are under study. The more novel microarrays and comparative genomic hybridization techniques, still ultimately rely on the quality of the DNA under analysis.

[0020] Several techniques are in development to measure multiple SNPs on the DNA of a small number of cells, a single cell (for example, a blastomere), a small number of chromosomes, or from fragments of DNA. There are techniques that use Polymerase Chain Reaction (PCR), followed by microarray genotyping analysis. Some PCR-based techniques include whole genome amplification (WGA) techniques such as multiple displacement amplification (MDA), and Molecular Inversion Probes (MIPS) that perform genotyping using multiple tagged oligonucleotides that may then be amplified using PCR with a singe pair of primers. An example of a non-PCR based technique is fluorescence in situ hybridization (FISH). It is apparent that the techniques will be severely error-prone due to the limited amount of genetic material which will exacerbate the impact of effects such as allele drop-outs, imperfect hybridization, and contamination.

[0021] Many techniques exist which provide genotyping data. Taqman is a unique genotyping technology produced and distributed by Applied Biosystems. Taqman uses polymerase chain reaction (PCR) to amplify sequences of interest. During PCR cycling, an allele specific minor groove binder (MGB) probe hybridizes to amplified sequences. Strand synthesis by the polymerase enzymes releases reporter dyes linked to the MGB probes, and then the Taqman optical readers detect the dyes. In this manner, Taqman achieves quantitative allelic discrimination. Compared with array based genotyping technologies, Taqman is quite expensive per reaction (˜$0.40/reaction), and throughput is relatively low (384 genotypes per run). While only 1 ng of DNA per reaction is necessary, thousands of genotypes by Taqman requires microgram quantities of DNA, so Taqman does not necessarily use less DNA than microarrays. However, with respect to the IVF genotyping workflow, Taqman is the most readily applicable technology. This is due to the high reliability of the assays and, most importantly, the speed and ease of the assay (˜3 hours per run and minimal molecular biological steps). Also unlike many array technologies (such as 500 k Affymetrix arrays), Taqman is highly customizable, which is important for the IVF market. Further, Taqman is highly quantitative, so anueploidies could be detected with this technology alone.

[0022] Illumina has recently emerged as a leader in high-throughput genotyping. Unlike Affymetrix, Illumina genotyping arrays do not rely exclusively on hybridization. Instead, Illumina technology uses an allele-specific DNA extension step, which is much more sensitive and specific than hybridization alone, for the original sequence detection. Then, all of these alleles are amplified in multiplex by PCR, and then these products hybridized to bead arrays. The beads on these arrays contain unique “address” tags, not native sequence, so this hybridization is highly specific and sensitive. Alleles are then called by quantitative scanning of the bead arrays. The Illlumina Golden Gate assay system genotypes up to 1536 loci concurrently, so the throughput is better than Taqman but not as high as Affymetrix 500 k arrays. The cost of Illumina genotypes is lower than Taqman, but higher than Affymetrix arrays. Also, the Illumina platform takes as long to complete as the 500 k Affymetrix arrays (up to 72 hours), which is problematic for IVF genotyping. However, Illumina has a much better call rate, and the assay is quantitative, so anueploidies are detectable with this technology. Illumina technology is much more flexible in choice of SNPs than 500 k Affymetrix arrays.

[0023] One of the highest throughput techniques, which allows for the measurement of up to 250,000 SNPs at a time, is the Affymetrix GeneChip 500K genotyping array. This technique also uses PCR, followed by analysis by hybridization and detection of the amplified DNA sequences to DNA probes, chemically synthesized at different locations on a quartz surface. A disadvantage of these arrays are the low flexibility and the lower sensitivity. There are modified approaches that can increase selectivity, such as the “perfect match” and “mismatch probe” approaches, but these do so at the cost of the number of SNPs calls per array.

[0024] Pyrosequencing, or sequencing by synthesis, can also be used for genotyping and SNP analysis. The main advantages to pyrosequencing include an extremely fast turnaround and unambiguous SNP calls, however, the assay is not currently conducive to high-throughput parallel analysis. PCR followed by gel electrophoresis is an exceedingly simple technique that has met the most success in preimplantation diagnosis. In this technique, researchers use nested PCR to amplify short sequences of interest. Then, they run these DNA samples on a special gel to visualize the PCR products. Different bases have different molecular weights, so one can determine base content based on how fast the product runs in the gel. This technique is low-throughput and requires subjective analyses by scientists using current technologies, but has the advantage of speed (1-2 hours of PCR, 1 hour of gel electrophoresis). For this reason, it has been used previously for prenatal genotyping for a myriad of diseases, including: thalassaemia, neurofibromatosis type 2, leukocyte adhesion deficiency type I, Hallopeau-Siemens disease, sickle-cell anemia, retinoblastoma, Pelizaeus-Merzbacher disease, Duchenne muscular dystrophy, and Currarino syndrome.

[0025] Another promising technique that has been developed for genotyping small quantities of genetic material with very high fidelity is Molecular Inversion Probes (MIPs), such as Affymetrix's Genflex Arrays. This technique has the capability to measure multiple SNPs in parallel: more than 10,000 SNPS measured in parallel have been verified. For small quantities of genetic material, call rates for this technique have been established at roughly 95%, and accuracy of the calls made has been established to be above 99%. So far, the technique has been implemented for quantities of genomic data as small as 150 molecules for a given SNP. However, the technique has not been verified for genomic data from a single cell, or a single strand of DNA, as would be required for pre-implantation genetic diagnosis.

[0026] The MIP technique makes use of padlock probes which are linear oligonucleotides whose two ends can be joined by ligation when they hybridize to immediately adjacent target sequences of DNA. After the probes have hybridized to the genomic DNA, a gap-fill enzyme is added to the assay which can add one of the four nucleotides to the gap. If the added nucleotide (A,C,T,G) is complementary to the SNP under measurement, then it will hybridize to the DNA, and join the ends of the padlock probe by ligation. The circular products, or closed padlock probes, are then differentiated from linear probes by exonucleolysis. The exonuclease, by breaking down the linear probes and leaving the circular probes, will change the relative concentrations of the closed vs. the unclosed probes by a factor of 1000 or more. The probes that remain are then opened at a cleavage site by another enzyme, removed from the DNA, and amplified by PCR. Each probe is tagged with a different tag sequence consisting of 20 base tags (16,000 have been generated), and can be detected, for example, by the Affymetrix GenFlex Tag Array. The presence of the tagged probe from a reaction in which a particular gap-fill enzyme was added indicates the presence of the complimentary amino acid on the relevant SNP.

[0027] The molecular biological advantages of MIPS include: (1) multiplexed genotyping in a single reaction, (2) the genotype “call” occurs by gap fill and ligation, not hybridization, and (3) hybridization to an array of universal tags decreases false positives inherent to most array hybridizations. In traditional 500K, TaqMan and other genotyping arrays, the entire genomic sample is hybridized to the array, which contains a variety of perfect match and mismatch probes, and an algorithm calls likely genotypes based on the intensities of the mismatch and perfect match probes. Hybridization, however, is inherently noisy, because of the complexities of the DNA sample and the huge number of probes on the arrays. MIPs, on the other hand, uses multiplex probes (i.e., not on an array) that are longer and therefore more specific, and then uses a robust ligation step to circularize the probe. Background is exceedingly low in this assay (due to specificity), though allele dropout may be high (due to poor performing probes).

[0028] When this technique is used on genomic data from a single cell (or small numbers of cells) it will—like PCR based approaches—suffer from integrity issues. For example, the inability of the padlock probe to hybridize to the genomic DNA will cause allele dropouts. This will be exacerbated in the context of in-vitro fertilization since the efficiency of the hybridization reaction is low, and it needs to proceed relatively quickly in order to genotype the embryo in a limited time period. Note that the hybridization reaction can be reduced well below vendor-recommended levels, and micro-fluidic techniques may also be used to accelerate the hybridization reaction.

[0029] These approaches to reducing the time for the hybridization reaction will result in reduced data quality.

[0030] Once the genetic data has been measured, the next step is to use the data for predictive purposes. Much research has been done in predictive genomics, which tries to understand the precise functions of proteins, RNA and DNA so that phenotypic predictions can be made based on genotype. Canonical techniques focus on the function of Single-Nucleotide Polymorphisms (SNP); but more advanced methods are being brought to bear on multi-factorial phenotypic features. These methods include techniques, such as linear regression and nonlinear neural networks, which attempt to determine a mathematical relationship between a set of genetic and phenotypic predictors and a set of measured outcomes. There is also a set of regression analysis techniques, such as Ridge regression, log regression and stepwise selection, that are designed to accommodate sparse data sets where there are many potential predictors relative to the number of outcomes, as is typical of genetic data, and which apply additional constraints on the regression parameters so that a meaningful set of parameters can be resolved even when the data is underdetermined. Other techniques apply principal component analysis to extract information from undetermined data sets. Other techniques, such as decision trees and contingency tables, use strategies for subdividing subjects based on their independent variables in order to place subjects in categories or bins for which the phenotypic outcomes are similar. A recent technique, termed logical regression, describes a method to search for different logical interrelationships between categorical independent variables in order to model a variable that depends on interactions between multiple independent variables related to genetic data. Regardless of the method used, the quality of the prediction is naturally highly dependant on the quality of the genetic data used to make the prediction.

[0031] Normal humans have two sets of 23 chromosomes in every diploid cell, with one copy coming from each parent. Aneuploidy, a cell with an extra or missing chromosomes, and uniparental disomy, a cell with two of a given chromosome that originate from one parent, are belived to be responsible for a large percentage of failed implantations, miscarriages, and genetic diseases. When only certain cells in an individual are aneuploid, the individual is said to exhibit mosaicism. Detection of chromosomal abnormalities can identify individuals or embryos with conditions such as Down syndrome, Klinefelters syndrome, and Turner syndrome, among others, in addition to increasing the chances of a successful pregnancy. Testing for chromosomal abnormalities is especially important as mothers age: between the ages of 35 and 40 it is estimated that between 40% and 50% of the embryos are abnormal, and above the age of 40, more than half of the embryos are abnormal.

[0032] Karyotyping, the traditional method used for the prediction of aneuploides and mosaicism is giving way to other more high throughput, more cost effective methods. One method that has attracted much attention recently is Flow cytometry (FC) and fluorescence in situ hybridization (FISH) which can be used to detect aneuploidy in any phase of the cell cycle. One advantage of this method is that it is less expensive than karyotyping, but the cost is significant enough that generally a small selection of chromosomes are tested (usually chromosomes 13, 18, 21, X, Y; also sometimes 8, 9, 15, 16, 17, 22); in addition, FISH has a low level of specificity. Using FISH to analyze 15 cells, one can detect mosaicism of 19% with 95% confidence. The reliability of the test becomes much lower as the level of mosaicism gets lower, and as the number of cells to analyze decreases. The test is estimated to have a false negative rate as high as 15% when a single cell is analysed. There is a great demand for a method that has a higher throughput, lower cost, and greater accuracy.

[0033] Listed here is a set of prior art which is related to the field of the current invention. None of this prior art contains or in any way refers to the novel elements of the current invention. In U.S. Pat. No. 6,720,140, Hartley et al describe a recombinational cloning method for moving or exchanging segments of DNA molecules using engineered recombination sites and recombination proteins. In U.S. Pat. No. 6,489,135 Parrott et al. provide methods for determining various biological characteristics of in vitro fertilized embryos, including overall embryo health, implantability, and increased likelihood of developing successfully to term by analyzing media specimens of in vitro fertilization cultures for levels of bioactive lipids in order to determine these characteristics. In US Patent Application 20040033596 Threadgill et al. describe a method for preparing homozygous cellular libraries useful for in vitro phenotyping and gene mapping involving site-specific mitotic recombination in a plurality of isolated parent cells. In U.S. Pat. No. 5,994,148 Stewart et al. describe a method of determining the probability of an in vitro fertilization (IVF) being successful by measuring Relaxin directly in the serum or indirectly by culturing granulosa lutein cells extracted from the patient as part of an IVF/ET procedure. In U.S. Pat. No. 5,635,366 Cooke et al. provide a method for predicting the outcome of IVF by determining the level of 11β-hydroxysteroid dehydrogenase (11β-HSD) in a biological sample from a female patient. In U.S. Pat. No. 7,058,616 Larder et al. describe a method for using a neural network to predict the resistance of a disease to a therapeutic agent. In U.S. Pat. No. 6,958,211 Vingerhoets et al. describe a method wherein the integrase genotype of a given HIV strain is simply compared to a known database of HIV integrase genotype with associated phenotypes to find a matching genotype. In U.S. Pat. No. 7,058,517 Denton et al. describe a method wherein an individual's haplotypes are compared to a known database of haplotypes in the general population to predict clinical response to a treatment. In U.S. Pat. No. 7,035,739 Schadt at al. describe a method is described wherein a genetic marker map is constructed and the individual genes and traits are analyzed to give a gene-trait locus data, which are then clustered as a way to identify genetically interacting pathways, which are validated using multivariate analysis. In U.S. Pat. No. 6,025,128 Veltri et al. describe a method involving the use of a neural network utilizing a collection of biomarkers as parameters to evaluate risk of prostate cancer recurrence.

[0034] The cost of DNA sequencing is dropping rapidly, and in the near future individual genomic sequencing for personal benefit will become more common. Knowledge of personal genetic data will allow for extensive phenotypic predictions to be made for the individual. In order to make accurate phenotypic predictions high quality genetic data is critical, whatever the context. In the case of prenatal or pre-implantation genetic diagnoses a complicating factor is the relative paucity of genetic material available. Given the inherently noisy nature of the measured genetic data in cases where limited genetic material is used for genotyping, there is a great need for a method which can increase the fidelity of, or clean, the primary data.

[0035] The system disclosed enables the cleaning of incomplete or noisy genetic data using secondary genetic data as a source of information, and also the determination of chromosome copy number using said genetic data. While the disclosure focuses on genetic data from human subjects, and more specifically on as-yet not implanted embryos or developing fetuses, as well as related individuals, it should be noted that the methods disclosed apply to the genetic data of a range of organisms, in a range of contexts. The techniques described for cleaning genetic data are most relevant in the context of pre-implantation diagnosis during in-vitro fertilization, prenatal diagnosis in conjunction with amniocentesis, chorion villus biopsy, fetal tissue sampling, and non-invasive prenatal diagnosis, where a small quantity of fetal genetic material is isolated from maternal blood. The use of this method may facilitate diagnoses focusing on inheritable diseases, chromosome copy number predictions, increased likelihoods of defects or abnormalities, as well as making predictions of susceptibility to various disease- and non-disease phenotypes for individuals to enhance clinical and lifestyle decisions. The invention addresses the shortcomings of prior art that are discussed above.

[0036] In one aspect of the invention, methods make use of knowledge of the genetic data of the mother and the father such as diploid tissue samples, sperm from the father, haploid samples from the mother or other embryos derived from the mother's and father's gametes, together with the knowledge of the mechanism of meiosis and the imperfect measurement of the embryonic DNA, in order to reconstruct, in silico, the embryonic DNA at the location of key loci with a high degree of confidence. In one aspect of the invention, genetic data derived from other related individuals, such as other embryos, brothers and sisters, grandparents or other relatives can also be used to increase the fidelity of the reconstructed embryonic DNA. It is important to note that the parental and other secondary genetic data allows the reconstruction not only of SNPs that were measured poorly, but also of insertions, deletions, and of SNPs or whole regions of DNA that were not measured at all.

[0037] In one aspect of the invention, the fetal or embryonic genomic data which has been reconstructed, with or without the use of genetic data from related individuals, can be used to detect if the cell is aneuploid, that is, where fewer or more than two of a particular chromosome is present in a cell. The reconstructed data can also be used to detect for uniparental disomy, a condition in which two of a given chromosome are present, both of which originate from one parent. This is done by creating a set of hypotheses about the potential states of the DNA, and testing to see which hypothesis has the highest probability of being true given the measured data. Note that the use of high throughput genotyping data for screening for aneuploidy enables a single blastomere from each embryo to be used both to measure multiple disease-linked loci as well as to screen for aneuploidy.

[0038] In another aspect of the invention, the direct measurements of the amount of genetic material, amplified or unamplified, present at a plurality of loci, can be used to detect for monosomy, uniparental disomy, trisomy and other aneuploidy states. The idea behind this method is that measuring the amount of genetic material at multiple loci will give a statistically significant result.

[0039] In another aspect of the invention, the measurements, direct or indirect, of a particular subset of SNPs, namely those loci where the parents are both homozygous but with different allele values, can be used to detect for chromosomal abnormalities by looking at the ratios of maternally versus paternally miscalled homozygous loci on the embryo. The idea behind this method is that those loci where each parent is homozygous, but have different alleles, by definition result in a heterozygous loci on the embryo. Allele drop outs at those loci are random, and a shift in the ratio of loci miscalled as homozygous can only be due to incorrect chromosome number.

[0040] It will be recognized by a person of ordinary skill in the art, given the benefit of this disclosure, that various aspects and embodiments of this disclosure may implemented in combination or separately.

Conceptual Overview of the System

[0061] The goal of the disclosed system is to provide highly accurate genomic data for the purpose of genetic diagnoses. In cases where the genetic data of an individual contains a significant amount of noise, or errors, the disclosed system makes use of the expected similarities between the genetic data of the target individual and the genetic data of related individuals, to clean the noise in the target genome. This is done by determining which segments of chromosomes of related individuals were involved in gamete formation and, when necessary where crossovers may have occurred during meiosis, and therefore which segments of the genomes of related individuals are expected to be nearly identical to sections of the target genome. In certain situations this method can be used to clean noisy base pair measurements on the target individual, but it also can be used to infer the identity of individual base pairs or whole regions of DNA that were not measured. It can also be used to determine the number of copies of a given chromosome segment in the target individual. In addition, a confidence may be computed for each call made. A highly simplified explanation is presented first, making unrealistic assumptions in order to illustrate the concept of the invention. A detailed statistical approach that can be applied to the technology of today is presented afterward.

[0062] In one aspect of the invention, the target individual is an embryo, and the purpose of applying the disclosed method to the genetic data of the embryo is to allow a doctor or other agent to make an informed choice of which embryo(s) should be implanted during IVF. In another aspect of the invention, the target individual is a fetus, and the purpose of applying the disclosed method to genetic data of the fetus is to allow a doctor or other agent to make an informed choice about possible clinical decisions or other actions to be taken with respect to the fetus.

• [0063] SNP (Single Nucleotide Polymorphism): a single nucleotide that may differ between the genomes of two members of the same species. In our usage of the term, we do not set any limit on the frequency with which each variant occurs.
• [0064] To call a SNP: to make a decision about the true state of a particular base pair, taking into account the direct and indirect evidence.
• [0065] Locus: a particular region of interest on the DNA of an individual, which may refer to a SNP, the site of a possible insertion or deletion, or the site of some other relevant genetic variation. Disease-linked SNPs may also refer to disease-linked loci.
• [0066] To call an allele: to determine the state of a particular locus of DNA. This may involve calling a SNP, or determining whether or not an insertion or deletion is present at that locus, or determining the number of insertions that may be present at that locus, or determining whether some other genetic variant is present at that locus.
• [0067] Correct allele call: An allele call that correctly reflects the true state of the actual genetic material of an individual.
• [0068] To clean genetic data: to take imperfect genetic data and correct some or all of the errors or fill in missing data at one or more loci. In the context of this disclosure, this involves using genetic data of related individuals and the method described herein.
• [0069] To increase the fidelity of allele calls: to clean genetic data.
• [0070] Imperfect genetic data: genetic data with any of the following: allele dropouts, uncertain base pair measurements, incorrect base pair measurements, missing base pair measurements, uncertain measurements of insertions or deletions, uncertain measurements of chromosome segment copy numbers, spurious signals, missing measurements, other errors, or combinations thereof.
• [0071] Noisy genetic data: imperfect genetic data, also called incomplete genetic data.
• [0072] Uncleaned genetic data: genetic data as measured, that is, where no method has been used to correct for the presence of noise or errors in the raw genetic data; also called crude genetic data.
• [0073] Confidence: the statistical likelihood that the called SNP, allele, set of alleles, or determined number of chromosome segment copies correctly represents the real genetic state of the individual.
• [0074] Parental Support (PS): a name sometimes used for the any of the methods disclosed herein, where the genetic information of related individuals is used to determine the genetic state of target individuals. In some cases, it refers specifically to the allele calling method, sometimes to the method used for cleaning genetic data, sometimes to the method to determine the number of copies of a segment of a chromosome, and sometimes to some or all of these methods used in combination.
• [0075] Copy Number Calling (CNC): the name given to the method described in this disclosure used to determine the number of chromosome segments in a cell.
• [0076] Qualitative CNC (also qCNC): the name given to the method in this disclosure used to determine chromosome copy number in a cell that makes use of qualitative measured genetic data of the target individual and of related individuals.
• [0077] Multigenic: affected by multiple genes, or alleles.
• [0078] Direct relation: mother, father, son, or daughter.
• [0079] Chromosomal Region: a segment of a chromosome, or a full chromosome.
• [0080] Segment of a Chromosome: a section of a chromosome that can range in size from one base pair to the entire chromosome.
• [0081] Section: a section of a chromosome. Section and segment can be used interchangeably.
• [0082] Chromosome: may refer to either a full chromosome, or also a segment or section of a chromosome.
• [0083] Copies: the number of copies of a chromosome segment may refer to identical copies, or it may refer to non-identical copies of a chromosome segment wherein the different copies of the chromosome segment contain a substantially similar set of loci, and where one or more of the alleles are different. Note that in some cases of aneuploidy, such as the M2 copy error, it is possible to have some copies of the given chromosome segment that are identical as well as some copies of the same chromosome segment that are not identical.
• [0084] Haplotypic Data: also called ‘phased data’ or ‘ordered genetic data;’ data from a single chromosome in a diploid or polyploid genome, i.e., either the segregated maternal or paternal copy of a chromosome in a diploid genome.
• [0085] Unordered Genetic Data: pooled data derived from measurements on two or more chromosomes in a diploid or polyploid genome, i.e., both the maternal and paternal copies of a chromosome in a diploid genome.
• [0086] Genetic data ‘in’, ‘of’, ‘at’ or ‘on’ an individual: These phrases all refer to the data describing aspects of the genome of an individual. It may refer to one or a set of loci, partial or entire sequences, partial or entire chromosomes, or the entire genome.
• [0087] Hypothesis: a set of possible copy numbers of a given set of chromosomes, or a set of possible genotypes at a given set of loci. The set of possibilities may contain one or more elements.
• [0088] Target Individual: the individual whose genetic data is being determined. Typically, only a limited amount of DNA is available from the target individual. In one context, the target individual is an embryo or a fetus.
• [0089] Related Individual: any individual who is genetically related, and thus shares haplotype blocks, with the target individual.
• [0090] Platform response: a mathematical characterization of the input/output characteristics of a genetic measurement platform, such as TAQMAN or INFINIUM. The input to the channel is the true underlying genotypes of the genetic loci being measured. The channel output could be allele calls (qualitative) or raw numerical measurements (quantitative), depending on the context. For example, in the case in which the platform's raw numeric output is reduced to qualitative genotype calls, the platform response consists of an error transition matrix that describes the conditional probability of seeing a particular output genotype call given a particular true genotype input. In the case in which the platform's output is left as raw numeric measurements, the platform response is a conditional probability density function that describes the probability of the numeric outputs given a particular true genotype input.
• [0091] Copy number hypothesis: a hypothesis about how many copies of a particular chromosome segment are in the embryo. In a preferred embodiment, this hypothesis consists of a set of sub-hypotheses about how many copies of this chromosome segment were contributed by each related individual to the target individual.

A Allele Calling: Preferred Method

[0092] Assume here the goal is to estimate the genetic data of an embryo as accurately as possible, and where the estimate is derived from measurements taken from the embryo, father, and mother across the same set of n SNPs. Note that where this description refers to SNPs, it may also refer to a locus where any genetic variation, such as a point mutation, insertion or deletion may be present. This allele calling method is part of the Parental Support (PS) system. One way to increase the fidelity of allele calls in the genetic data of a target individual for the purposes of making clinically actionable predictions is described here. It should be obvious to one skilled in the art how to modify the method for use in contexts where the target individual is not an embryo, where genetic data from only one parent is available, where neither, one or both of the parental haplotypes are known, or where genetic data from other related individuals is known and can be incorporated.

[0093] For the purposes of this discussion, only consider SNPs that admit two allele values; without loss of generality it is possible to assume that the allele values on all SNPs belong to the alphabet A={A,C}. It is also assumed that the errors on the measurements of each of the SNPs are independent. This assumption is reasonable when the SNPs being measured are from sufficiently distant genic regions. Note that one could incorporate information about haplotype blocks or other techniques to model correlation between measurement errors on SNPs without changing the fundamental concepts of this invention.

[0094] Let e=(e₁,e₂) be the true, unknown, ordered SNP information on the embryo, e₁,e₂∈Aⁿ. Define e₁ to be the genetic haploid information inherited from the father and e₂ to be the genetic haploid information inherited from the mother. Also use e_i=(e_1i,e_2i) to denote the ordered pair of alleles at the i-th position of e. In similar fashion, let f=(f₁,f₂) and m=(m₁,m₂) be the true, unknown, ordered SNP information on the father and mother respectively. In addition, let g₁ be the true, unknown, haploid information on a single sperm from the father. (One can think of the letter g as standing for gamete. There is no g₂. The subscript is used to remind the reader that the information is haploid, in the same way that f₁ and f₂ are haploid.) It is also convenient to define r=(f,m), so that there is a symbol to represent the complete set of diploid parent information from which e inherits, and also write r_i=(f_i, m_i)=((f_1i,f_2i),(m_1i,m_2i)) to denote the complete set of ordered information on father and mother at the i-th SNP. Finally, let ê=(ê₁, ê₂) be the estimate of e that is sought, ê₁, ê₂ϵAⁿ.

[0095] By a crossover map, it is meant an n-tuple θϵ{1,2}ⁿ that specifies how a haploid pair such as (f₁,f₂) recombines to form a gamete such as e₁. Treating θ as a function whose output is a haploid sequence, define θ(f)_i=θ(f₁,f₂)_i=f_θi,i. To make this idea more concrete, let f₁=ACAAACCC, let f₂=CAACCACA, and let θ=11111222. Then θ(f₁,f₂)=ACAAAACA. In this example, the crossover map θ implicitly indicates that a crossover occurred between SNPs i=5 and i=6.

[0096] Formally, let θ be the true, unknown crossover map that determines e₁ from f, let ϕ be the true, unknown crossover map that determines e₂ from m, and let ψ be the true, unknown crossover map that determines g₁ from f. That is, e₁=θ(f), e₂=ϕ(m), g₁=ψ(f). It is also convenient to define X=(θ,ϕ,ψ) so that there is a symbol to represent the complete set of crossover information associated with the problem. For simplicity sake, write e=X(r) as shorthand for e=(θ(f),ϕ(m)); also write e_i=X(r_i) as shorthand for e_i=X(r)_i

[0097] In reality, when chromosomes combine, at most a few crossovers occur, making most of the 2ⁿ theoretically possible crossover maps distinctly improbable. In practice, these very low probability crossover maps will be treated as though they had probability zero, considering only crossover maps belonging to a comparatively small set Ω. For example, if Ω is defined to be the set of crossover maps that derive from at most one crossover, then |Ω|=2n.

[0098] It is convenient to have an alphabet that can be used to describe unordered diploid measurements. To that end, let B={A,B,C,X}. Here A and C represent their respective homozygous locus states and B represents a heterozygous but unordered locus state. Note: this section is the only section of the document that uses the symbol B to stand for a heterozygous but unordered locus state. Most other sections of the document use the symbols A and B to stand for the two different allele values that can occur at a locus. X represents an unmeasured locus, i.e., a locus drop-out. To make this idea more concrete, let f₁=ACAAACCC, and let f₂=CAACCACA. Then a noiseless unordered diploid measurement of {tilde over (f)} would yield f=BBABBBCB.

[0099] In the problem at hand, it is only possible to take unordered diploid measurements of e, f, and m, although there may be ordered haploid measurements on g₁. This results in noisy measured sequences that are denoted {tilde over (e)}ϵBⁿ, {tilde over (f)}ϵBⁿ, {tilde over (m)}ϵBⁿ, and {tilde over (g)}₁ϵAⁿ respectively. It will be convenient to define {tilde over (r)}=({tilde over (f)}, {tilde over (m)}) so that there is a symbol that represents the noisy measurements on the parent data. It will also be convenient to define {tilde over (D)}=({tilde over (r)}, {tilde over (e)}, {tilde over (g)}₁) so that there is a symbol to represent the complete set of noisy measurements associated with the problem, and to write {tilde over (D)}_i=({tilde over (r)}_i, {tilde over (e)}_i, {tilde over (g)}_1i)=({tilde over (f)}_i, {tilde over (m)}_i, {tilde over (e)}_i, {tilde over (g)}_1i) to denote the complete set of measurements on the i-th SNP. (Please note that, while f_i is an ordered pair such as (A,C), {tilde over (f)}_i is a single letter such as B.)

[0100] Because the diploid measurements are unordered, nothing in the data can distinguish the state (f₁, f₂) from (f₂, f₁) or the state (m₁, m₂) from (m₂, m₁). These indistinguishable symmetric states give rise to multiple optimal solutions of the estimation problem. To eliminate the symmetries, and without loss of generality, assign θ₁=ϕ₁=1.

[0101] In summary, then, the problem is defined by a true but unknown underlying set of information {r, e, g₁, X}, with e=X(r). Only noisy measurements {tilde over (D)}=({tilde over (r)}, {tilde over (e)}, {tilde over (g)}₁) are available. The goal is to come up with an estimate ê of e, based on {tilde over (D)}.

[0102] Note that this method implicitly assumes euploidy on the embryo. It should be obvious to one skilled in the art how this method could be used in conjunction with the aneuploidy calling methods described elsewhere in this patent. For example, the aneuploidy calling method could be first employed to ensure that the embryo is indeed euploid and only then would the allele calling method be employed, or the aneuploidy calling method could be used to determine how many chromosome copies were derived from each parent and only then would the allele calling method be employed. It should also be obvious to one skilled in the art how this method could be modified in the case of a sex chromosome where there is only one copy of a chromosome present.

Solution Via Maximum a Posteriori Estimation

[0103] In one embodiment of the invention, it is possible, for each of the n SNP positions, to use a maximum a posteriori (MAP) estimation to determine the most probable ordered allele pair at that position. The derivation that follows uses a common shorthand notation for probability expressions. For example, P(e′_i,{tilde over (D)}|X′) is written to denote the probability that random variable e_i takes on value e′_i and the random variable {tilde over (D)} takes on its observed value, conditional on the event that the random variable X takes on the value X′. Using MAP estimation, then, the i-th component of ê, denoted ê_i=(ê_1i, ê_2i) is given by

[0104] $\begin{array}{cc}{\hat{e}}_i& =\underset{e_i^{\prime }}{\arg \max }P\left(e_i^{\prime }❘\tilde{D}\right)\\ & =\underset{e_i^{\prime }}{\arg \max }P\left(e_i^{\prime },\tilde{D}\right)\\ & =\underset{e_i^{\prime }}{\arg \max }\sum _{X^{\prime }\in {\Omega }^3}P\left(X^{\prime }\right)P\left(e_i^{\prime },\tilde{D}❘X^{\prime }\right)\\ \left(a\right)& =\underset{e_i^{\prime }}{\arg \max }\sum _{X^{\prime }\in {\Omega }^3:{\theta }_1^{\prime }={\varphi }_1^{\prime }=1}P\left(X^{\prime }\right)P\left(e_i^{\prime },{\tilde{D}}_i❘X^{\prime }\right)\prod _{j\ne i}P\left({\tilde{D}}_j❘X^{\prime }\right)\\ \left(b\right)& =\underset{e_i^{\prime }}{\arg \max }\sum _{X^{\prime }\in {\Omega }^3:{\theta }_1^{\prime }={\varphi }_1^{\prime }=1}P\left(X^{\prime }\right)\sum _{r_i^{\prime }\in A^4}P\left(r_i^{\prime }\right)P\left(e_i^{\prime },{\tilde{D}}_i❘X^{\prime },r_i^{\prime }\right)\prod _{j\ne i}\sum _{r_j^{\prime }\in A^4}P\left(r_j^{\prime }\right)P\left({\tilde{D}}_j❘X^{\prime },r_j^{\prime }\right)\\ \left(c\right)& =\underset{e_i^{\prime }}{\arg \max }\sum _{X^{\prime }\in {\Omega }^3:{\theta }_1^{\prime }={\varphi }_1^{\prime }=1}P\left(X^{\prime }\right)\sum _{r_i^{\prime }\in A^4}P\left(r_i^{\prime }\right)P\left(e_i^{\prime }❘X^{\prime },r_i^{\prime }\right)P\left({\tilde{D}}_i❘X^{\prime },r_i^{\prime }\right)\prod _{j\ne i}\sum _{r_j^{\prime }\in A^4}P\left(r_j^{\prime }\right)P\left({\tilde{D}}_j❘X^{\prime },r_j^{\prime }\right)\\ \mathrm{(*})& =\underset{e_i^{\prime }}{\arg \max }\sum _{X^{\prime }\in {\Omega }^3:{\theta }_1^{\prime }={\varphi }_1^{\prime }=1}P\left(X^{\prime }\right)\prod _j\sum _{r_j^{\prime }\in A^4}1\left(i\ne j\mathrm{or}{X}^{\prime }\left(r_j^{\prime }\right)=e_i^{\prime }\right)P\left(r_j^{\prime }\right)P\left({\tilde{D}}_j❘X^{\prime },r_j^{\prime }\right)\end{array}$
In the preceding set of equations, (a) holds because the assumption of SNP independence means that all of the random variables associated with SNP i are conditionally independent of all of the random variables associated with SNP j, given X; (b) holds because r is independent of X; (c) holds because e_i and {tilde over (D)}_i are conditionally independent given r_i and X (in particular, e_i=X(r_i)); and (*) holds, again, because e_i=X(r_i), which means that P(e′_i|X′,r′_i) evaluates to either one or zero and hence effectively filters r′_i to just those values that are consistent with e′_i and X′.

[0105] The final expression (*) above contains three probability expressions: P(X′), P(r′_j), and P({tilde over (D)}_j|X′,r′_j). The computation of each of these quantities is discussed in the following three sections.

Crossover Map Probabilities

[0106] Recent research has enabled the modeling of the probability of recombination between any two SNP loci. Observations from sperm studies and patterns of genetic variation show that recombination rates vary extensively on kilobase scales and that much recombination occurs in recombination hotspots. The NCBI data about recombination rates on the Human Genome is publicly available through the UCSC Genome Annotation Database.

[0107] One may use the data set from the Hapmap Project or the Perlegen Human Haplotype Project. The latter is higher density; the former is higher quality. These rates can be estimated using various techniques known to those skilled in the art, such as the reversible-jump Markov Chain Monte Carlo (MCMC) method that is available in the package LDHat.

[0108] In one embodiment of the invention, it is possible to calculate the probability of any crossover map given the probability of crossover between any two SNPs. For example, P(θ=11111222) is one half the probability that a crossover occurred between SNPs five and six. The reason it is only half the probability is that a particular crossover pattern has two crossover maps associated with it: one for each gamete. In this case, the other crossover map is θ=22222111.

[0109] Recall that X=(θ,ϕ,ψ), where e₁=θ(f), e₂=ϕ(m), g₁=ψ(f). Obviously θ, ϕ, and ψ result from independent physical events, so P(X)=P(θ)P(ϕ)P(ψ). Further assume that P_θ(^●)=P_ϕ(^●)=P_ψ(^●), where the actual distribution P_θ(^●) is determined in the obvious way from the Hapmap data.

Allele Probabilities

[0110] It is possible to determine P(r_i)=P(f_i)P(m_i)=P(f_i1)P(f_i2)P(m_i1)P(m_i2) using population frequency information from databases such as dbSNP. Also, as mentioned previously, choose SNPs for which the assumption of intra-haploid independence is a reasonable one. That is, assume that

[0111] $P\left(r\right)=\prod _{i}P\left(r_i\right)$
Measurement Errors

[0112] Conditional on whether a locus is heterozygous or homozygous, measurement errors may be modeled as independent and identically distributed across all similarly typed loci. Thus:

[0113] $\begin{array}{c}P\left(\tilde{D}❘X,r\right)=\prod _{i}P\left({\tilde{D}}_i❘X,r_i\right)\\ =\prod _{i}P\left({\tilde{f}}_i,{\tilde{m}}_i,{\tilde{e}}_i,{\tilde{g}}_{1i}❘X,f_i,m_i\right)\\ =\prod _{i}P\left({\tilde{f}}_i❘f_i\right)P\left({\tilde{m}}_i❘m_i\right)P\left({\tilde{e}}_i❘\theta \left(f_i\right),\varphi \left(m_i\right)\right)P\left({\tilde{g}}_{1i}❘\Psi \left(f_i\right)\right)\end{array}$
where each of the four conditional probability distributions in the final expression is determined empirically, and where the additional assumption is made that the first two distributions are identical. For example, for unordered diploid measurements on a blastomere, empirical values p_d=0.5 and p_a=0.02 are obtained, which lead to the conditional probability distribution for P({tilde over (e)}_i|e_i) shown in Table 1.

[0114] Note that the conditional probability distributions mentioned above, P({tilde over (f)}_i|f_i), P({tilde over (m)}_i|m_i), P({tilde over (e)}_i|e_i), can vary widely from experiment to experiment, depending on various factors in the lab such as variations in the quality of genetic samples, or variations in the efficiency of whole genome amplification, or small variations in protocols used. Therefore, in a preferred embodiment, these conditional probability distributions are estimated on a per-experiment basis. We focus in later sections of this disclosure on estimating P({tilde over (e)}_i|e_i), but it will be clear to one skilled in the art after reading this disclosure how similar techniques can be applied to estimating P({tilde over (f)}_i|f_i) and P({tilde over (m)}_i|m_i). The distributions can each be modeled as belonging to a parametric family of distributions whose particular parameter values vary from experiment to experiment. As one example among many, it is possible to implicitly model the conditional probability distribution P({tilde over (e)}_i|e_i) as being parameterized by an allele dropout parameter p_d and an allele dropin parameter p_a. The values of these parameters might vary widely from experiment to experiment, and it is possible to use standard techniques such as maximum likelihood estimation, MAP estimation, or Bayesian inference, whose application is illustrated at various places in this document, to estimate the values that these parameters take on in any individual experiment. Regardless of the precise method one uses, the key is to find the set of parameter values that maximizes the joint probability of the parameters and the data, by considering all possible tuples of parameter values within a region of interest in the parameter space. As described elsewhere in the document, this approach can be implemented when one knows the chromosome copy number of the target genome, or when one doesn't know the copy number call but is exploring different hypotheses. In the latter case, one searches for the combination of parameters and hypotheses that best match the data are found, as is described elsewhere in this disclosure.

[0115] Note that one can also determine the conditional probability distributions as a function of particular parameters derived from the measurements, such as the magnitude of quantitative genotyping measurements, in order to increase accuracy of the method. This would not change the fundamental concept of the invention.

[0116] It is also possible to use non-parametric methods to estimate the above conditional probability distributions on a per-experiment basis. Nearest neighbor methods, smoothing kernels, and similar non-parametric methods familiar to those skilled in the art are some possibilities. Although this disclosure focuses parametric estimation methods, use of non-parametric methods to estimate these conditional probability distributions would not change the fundamental concept of the invention. The usual caveats apply: parametric methods may suffer from model bias, but have lower variance. Non-parametric methods tend to be unbiased, but will have higher variance.

[0117] Note that it should be obvious to one skilled in the art, after reading this disclosure, how one could use quantitative information instead of explicit allele calls, in order to apply the PS method to making reliable allele calls, and this would not change the essential concepts of the disclosure.

B Factoring the Allele Calling Equation

[0118] In a preferred embodiment of the invention, the algorithm for allele calling can be structured so that it can be executed in a more computationally efficient fashion. In this section the equations are re-derived for allele-calling via the MAP method, this time reformulating the equations so that they reflect such a computationally efficient method of calculating the result.

Notation

[0119] X*,Y*,Z*ϵ{A,C}^nx2 are the true ordered values on the mother, father, and embryo respectively.

[0120] H*ϵ{A,C}^nxh are true values on h sperm samples.

[0121] B*ϵ{A,C}^nxbx2 are true ordered values on b blastomeres.

[0122] D={x,y,z,B,H} is the set of unordered measurement data on father, mother, embryo, b blastomeres and h sperm samples. D_i={x_i,y_i,z_i,H_i,B_i,} is the data set restricted to the i-th SNP.

[0123] rϵ{A,C}⁴ represents a candidate 4-tuple of ordered values on both the mother and father at a particular locus.

[0124] {circumflex over (Z)}_i ϵ{A,C}² is the estimated ordered embryo value at SNP i.

[0125] Q=(2+2b+h) is the effective number of haploid chromosomes being measured, excluding the parents. Any hypothesis about the parental origin of all measured data (excluding the parents themselves) requires that Q crossover maps be specified.

[0126] χϵ{1,2}^n×Q is a crossover map matrix, representing a hypothesis about the parental origin of all measured data, excluding the parents. Note that there are 2^nQ different crossover matrices. χ_i χ_i, is the matrix restricted to the i-th row. Note that there are 2^Q vector values that the i-th row can take on, from the set χϵ{1,2}^Q.

[0127] f(x; y, z) is a function of (x, y, z) that is being treated as a function of just x. The values behind the semi-colon are constants in the context in which the function is being evaluated.

PS Equation Factorization

[0128] $\begin{array}{c}{\hat{Z}}_i=\underset{z_i}{\arg \max }P\left(Z_i,D\right)\\ =\underset{z_i}{\arg \max }\sum _{\chi }P\left(\chi \right)P\left(Z_i,D|\chi \right)\\ =\underset{z_i}{\arg \max }\sum _{\chi }P\left({\chi }_1\right)P\left({\chi }_2|{\chi }_1\right)\dots P\left({\chi }_n,{\chi }_{n-1}\right)\\ \left(\sum _{r\in {\left(A,C\right)}^4}P\left(r\right)P\left(Z_i,D_i|X_i,r\right)\right)\prod _{j\ne 1}\left(\sum _{r\in {\left(A,C\right)}^4}P\left(r\right)P\left(D_j|X_j,r\right)\right)\\ =\underset{z_i}{\arg \max }\sum _{\chi }P\left({\chi }_1\right)P\left({\chi }_2|{\chi }_1\right)\dots P\left({\chi }_n,{\chi }_{n-1}\right)f_1\left({\chi }_i;Z_i,D_i\right)\end{array}$ $\begin{array}{c}\prod _{j=1}{f}_2\left({\chi }_j;D_j\right)=\underset{z_i}{\arg \max }\sum _{{\chi }_1\in {\left\{1,2\right\}}^Q}\dots \sum _{{\chi }_2\in {\left\{1,2\right\}}^Q}P\left({\chi }_1\right)P\left({\chi }_2|{\chi }_1\right)\dots \\ P\left({\chi }_n,{\chi }_{n-1}\right)\\ f_1\left({\chi }_i;Z_i,D_i\right)\prod _{j=1}{f}_2\left({\chi }_j;D_j\right)\\ =\underset{z_i}{\arg \max }\sum _{{\chi }_1\in {\left\{1,2\right\}}^Q}P\left({\chi }_1\right)f_2\left({\chi }_1;D_1\right)\times \\ \sum _{{\chi }_2\in {\left\{1,2\right\}}^Q}P\left({\chi }_2|{\chi }_1\right)\\ f_2\left({\chi }_2;D_2\right)\times \dots \sum _{{\chi }_1\in {\left\{1,2\right\}}^Q}P\left({\chi }_i|X_{i-1}\right)\\ f_1\left({\chi }_i;Z_i,D_i\right)\times \sum _{{\chi }_n\in {\left\{1,2\right\}}^Q}P\left({\chi }_n|X_{n-1}\right)\\ f_2\left({\chi }_n;D_n\right)\end{array}$
The number of different crossover matrices χ is 2^nQ. Thus, a brute-force application of the first line above is U(n2^nQ). By exploiting structure via the factorization of P(χ) and P(z_i,D|χ), and invoking the previous result, final line gives an expression that can be computed in O(n2^2Q).
C Quantitative Detection of Aneuploidy

[0129] In one embodiment of the invention, aneuploidy can be detected using the quantitative data output from the PS method discussed in this patent. Disclosed herein are multiple methods that make use of the same concept; these methods are termed Copy Number Calling (CNC). The statement of the problem is to determine the copy number of each of 23 chromosome-types in a single cell. The cell is first pre-amplified using a technique such as whole genome amplification using the MDA method. Then the resulting genetic material is selectively amplified with a technique such as PCR at a set of n chosen SNPs at each of m=23 chromosome types.

[0130] This yields a data set [t_ij], i=1 . . . n, j=1 . . . m of regularized ct (ct, or CT, is the point during the cycle time of the amplification at which dye measurement exceeds a given threshold) values obtained at SNP i, chromosome j. A regularized ct value implies that, for a given (i,j), the pair of raw ct values on channels FAM and VIC (these are arbitrary channel names denoting different dyes) obtained at that locus are combined to yield a ct value that accurately reflects the ct value that would have been obtained had the locus been homozygous. Thus, rather than having two ct values per locus, there is just one regularized ct value per locus.

[0131] The goal is to determine the set {n_j} of copy numbers on each chromosome. If the cell is euploid, then n_j=² for all j; one exception is the case of the male X chromosome. If n_j≠2 for at least one j, then the cell is aneuploid; excepting the case of male X.

Biochemical Model

[0132] The relationship between ct values and chromosomal copy number is modeled as follows: α_ijn_jQ2^βijtij-Q_T In this expression, n_j is the copy number of chromosome j. Q is an abstract quantity representing a baseline amount of pre-amplified genetic material from which the actual amount of pre-amplified genetic material at SNP i, chromosome j can be calculated as α_ijn_jQ. α_ij is a preferential amplification factor that specifies how much more SNP i on chromosome j will be pre-amplified via MDA than SNP 1 on chromosome 1. By definition, the preferential amplification factors are relative to

[0133] ${\alpha }_{11}\stackrel{\Delta }{=}1.$

[0134] β_ij is the doubling rate for SNP i chromosome j under PCR. t_ij is the ct value. Q_T is the amount of genetic material at which the ct value is determined. T is a symbol, not an index, and merely stands for threshold.

[0135] It is important to realize that α_ij, β_ij, and Q_T are constants of the model that do not change from experiment to experiment. By contrast, n_j and Q are variables that change from experiment to experiment. Q is the amount of material there would be at SNP 1 of chromosome 1, if chromosome 1 were monosomic.

[0136] The original equation above does not contain a noise term. This can be included by rewriting it as follows:

[0137] $\mathrm{(*}){\beta }_{\mathrm{ij}}{t}_{\mathrm{ij}}=\log \frac{Q_T}{{\alpha }_{\mathrm{ij}}}-\log n_j-\log Q+Z_{\mathrm{ij}}$
The above equation indicates that the ct value is corrupted by additive Gaussian noise Z_ij. Let the variance of this noise term be σ_ij².
Maximum Likelihood (ML) Estimation of Copy Number

[0138] In one embodiment of the method, the maximum likelihood estimation is used, with respect to the model described above, to determine n_j. The parameter Q makes this difficult unless another constraint is added:

[0139] $\frac{1}{m}\sum _j\log n_j=1$
This indicates that the average copy number is 2, or, equivalently, that the average log copy number is 1. With this additional constraint one can now solve the following ML problem:

[0140] $\begin{array}{c}\hat{Q},{\hat{n}}_j=\underset{Q,n_j}{\arg \max }\prod _{ij}{f}_z\left(\begin{array}{c}\log n_j+\log Q-\\ \left(\log \frac{Q_T}{{\alpha }_{ij}}-{\beta }_{ij}{t}_{ij}\right)\end{array}\right)s.t.\frac{1}{m}\sum _j\log n_j=1\\ =\underset{Q,n_j}{\arg \min }\sum _{ij}\frac{1}{{\sigma }_{ij}^2}{\left(\begin{array}{c}\log n_j+\log Q-\\ \left(\log \frac{Q_T}{{\alpha }_{ij}}-{\beta }_{ij}{t}_{ij}\right)\end{array}\right)}^{2}s.t.\frac{1}{m}\sum _j\log n_j=1\end{array}$
The last line above is linear in the variables log n_j and log Q, and is a simple weighted least squares problem with an equality constraint. The solution can be obtained in closed form by forming the Lagrangian

[0141] $L\left(\log n_j,\log Q\right)=\sum _{ij}\frac{1}{{\sigma }_{ij}^2}{\left(\begin{array}{c}\log n_j+\log Q-\\ \left(\log \frac{Q_T}{{\alpha }_{ij}}-{\beta }_{ij}{t}_{ij}\right)\end{array}\right)}^2+\lambda \sum _j\log n_j$
and taking partial derivatives.
Solution when Noise Variance is Constant

[0142] To avoid unnecessarily complicating the exposition, set σ_ij²=1. This assumption will remain unless explicitly stated otherwise. (In the general case in which each σ_ij² is different, the solutions will be weighted averages instead of simple averages, or weighted least squares solutions instead of simple least squares solutions.) In that case, the above linear system has the solution:

[0143] $\log Q_j\stackrel{\Delta }{=}\frac{1}{n}\sum _i\left(\log \frac{Q_T}{{\alpha }_{ij}}-{\beta }_{ij}{t}_{ij}\right)$ $\log Q=\frac{1}{m}\sum _j\log Q_j-1$ $\log n_j=\log Q_j-\log Q=\log \frac{Q_j}{Q}$
The first equation can be interpreted as a log estimate of the quantity of chromosome j. The second equation can be interpreted as saying that the average of the Q_j is the average of a diploid quantity; subtracting one from its log gives the desired monosome quantity. The third equation can be interpreted as saying that the copy number is just the ratio

[0144] $\frac{Q_j}{Q}.$
Note that n_j is a ‘double difference’, since it is a difference of Q-values, each of which is itself a difference of values.
Simple Solution

[0145] The above equations also reveal the solution under simpler modeling assumptions: for example, when making the assumption α_ij=1 for all i and j and/or when making the assumption that β_ij=β for all i and j. In the simplest case, when both α_ij=1 and β_ij=β, the solution reduces to

[0146] $\mathrm{(*}\mathrm{*)}\log n_j=1+\beta \left(\frac{1}{mn}\sum _{ij}{t}_{ij}-\frac{1}{n}\sum _i{t}_{ij}\right)$
The Double Differencing Method

[0147] In one embodiment of the invention, it is possible to detect monosomy using double differencing. It should be obvious to one skilled in the art how to modify this method for detecting other aneuploidy states. Let {t_ij} be the regularized ct values obtained from MDA pre-amplification followed by PCR on the genetic sample. As always, t_ij is the ct value on the i-th SNP of the j-th chromosome. Denote by tj the vector of ct values associated with the j-th chromosome. Make the following definitions:

[0148] $\stackrel{\_}{t}\stackrel{\Delta }{=}\frac{1}{mn}\sum _{i,j}{t}_{ij}$ ${\tilde{t}}_j\stackrel{\Delta }{=}{t}_j-\stackrel{\_}{t}1$
Classify chromosome j as monosomic if and only if f^T {tilde over (t)}_i is higher than a certain threshold value, where f is a vector that represents a monosomy signature. f is the matched filter, whose construction is described next.

[0149] The matched filter f is constructed as a double difference of values obtained from two controlled experiments. Begin with known quantities of euploid male genetic data and euploid female genetic material. Assume there are large quantities of this material, and pre-amplification can be omitted. On both the male and female material, use PCR to sequence n SNPs on both the X chromosome (chromosome 23), and chromosome 7. Let {t_ij^X}, i=1 . . . n, jϵ{7, 23} denote the measurements on the female, and let {t_ij^Y} similarly denote the measurements on the male. Given this, it is possible to construct the matched filter f from the resulting data as follows:

[0150] ${\stackrel{\_}{t}}_7^X\stackrel{\Delta }{=}\frac{1}{n}\sum _i^{}{t}_{i,7}^X$ ${\stackrel{\_}{t}}_7^Y\stackrel{\Delta }{=}\frac{1}{n}\sum _i^{}{t}_{i,7}^Y$ ${\Delta }^X\stackrel{\Delta }{=}{t}_{23}^X-{\stackrel{\_}{t}}_7^{X}1$ ${\Delta }^Y\stackrel{\Delta }{=}{t}_{23}^Y-{\stackrel{\_}{t}}_7^{Y}1$ $f\stackrel{\Delta }{=}{\Delta }^Y-{\Delta }^X$
In the above, equations, t₇^X and t₇^Y are scalars, while Δ^X and Δ^Y are vectors. Note that the superscripts X and Y are just symbolic labels, not indices, denoting female and male respectively. Do not to confuse the superscript X with measurements on the X chromosome. The X chromosome measurements are the ones with subscript 23.

[0151] The next step is to take noise into account and to see what remnants of noise survive in the construction of the matched filter f as well as in the construction of {tilde over (t)}_j. In this section, consider the simplest possible modeling assumption: that β_ij=β for all i and j, and that α_ij=1 for all i and j. Under these assumptions, from (*) above: βt_ij-log Q_T-log n_j-log Q+Z_ij Which can be rewritten as:

[0152] $t_{ij}=\frac{1}{\beta }\log Q_T-\frac{1}{\beta }\log n_j-\frac{1}{\beta }\log Q+Z_{ij}$
In that case, the i-th component of the matched filter f is given by:

[0153] $f_i\stackrel{\Delta }{=}{\Delta }_i^Y-{\Delta }_i^X=\left\{t_{i,23}^Y-{\stackrel{\_}{t}}_7^Y\right\}-\text{{}\text{(}{t}_{i,23}^X-{\stackrel{\_}{t}}_7^X\text{}}=\left\{\left(\frac{1}{\beta }\log Q_T-\frac{1}{{\beta }_{i,23}}\log n_{23}^Y-\frac{1}{\beta }\log Q^Y+Z_{i,23}^Y\right)-\frac{1}{n}\sum _i^{}\left(\frac{1}{\beta }\log Q_T-\frac{1}{\beta }\log n_7^Y-\frac{1}{\beta }\log Q^Y+Z_{i,7}^Y\right)\right\}-\left\{\left(\frac{1}{\beta }\log Q_T-\frac{1}{\beta }\log n_{23}^X-1\frac{1}{\beta }\log Q^X+Z_{i,23}^X\right)-\frac{1}{n}\sum _i^{}\left(\frac{1}{\beta }\log Q_T-\frac{1}{\beta }\log n_7^X-\frac{1}{\beta }\log Q^X+Z_{i,7}^X\right)\right\}-\left\{\left(\frac{1}{\beta }+Z_{i,23}^Y\right)-\frac{1}{n}\sum _i^{}{Z}_{i,7}^Y\right\}-\left\{Z_{i,23}^X-\frac{1}{n}\sum _i^{}{Z}_{i,7}^X\right\}$
Note that the above equations take advantage of the fact that all the copy number variables are known, for example, n₂₃^Y=1 and that n₂₃^X=2.

[0154] Given that all the noise terms are zero mean, the ideal matched filter is 1/β1. Further, since scaling the filter vector doesn't really change things, the vector 1 can be used as the matched filter. This is equivalent to simply taking the average of the components of {tilde over (t)}_j. In other words, the matched filter paradigm is not necessary if the underlying biochemistry follows the simple model. In addition, one may omit the noise terms above, which can only serve to lower the accuracy of the method. Accordingly, this gives:

[0155] ${\tilde{t}}_{ij}\stackrel{\Delta }{=}{t}_j-\stackrel{\_}{t}=\left\{\frac{1}{\beta }\log Q_T-\frac{1}{\beta }\log n_j-\frac{1}{\beta }\log Q+Z_{ij}\right\}-\frac{1}{mn}\sum _{i,j}\left\{\frac{1}{\beta }\log Q_T-\frac{1}{\beta }\log n_j-\frac{1}{\beta }\log Q+Z_{ij}\right\}=\frac{1}{\beta }\left(1-\log n_j\right)+Z_{ij}-\frac{1}{mn}\sum _{i,j}{Z}_{ij}$
In the above, it is assumed that

[0156] $\frac{1}{mn}\sum _{i,j}^{}\log n_j=1.$
that is, that the average copy number is 2.

[0157] Each element of the vector is an independent measurement of the log copy number (scaled by 1/β), and then corrupted by noise. The noise term Z_ij cannot be gotten rid of: it is inherent in the measurement. The second noise term probably cannot be gotten rid of either, since subtracting out t is necessary to remove the nuisance term

[0158] $\frac{1}{\beta }\log Q.$
Again, note that, given the observation that each element of {tilde over (t)}_j is an independent measurement of

[0159] $\frac{1}{\beta }\left(1-\log n_j\right),$
it is clear that a UMVU (uniform minimum variance unbiased) estimate of

[0160] $\frac{1}{\beta }\left(1-\log n_j\right)$
is just the average of the elements of {tilde over (t)}_j. (In the case in which each σ_ij² is different, it will be a weighted average.) Thus, performing a little bit of algebra, the UMVU estimator for log n_j is given by:

[0161] $\frac{1}{n}\sum _i^{}{\tilde{t}}_{ij}\approx \frac{1}{\beta }\left(1-\log n_j\right)\Rightarrow \log n_j\approx 1-\beta ·\frac{1}{n}\sum _{ij}^{}{t}_{ij}=1-\beta \left(\frac{1}{n}\sum _i^{}{t}_{ij}-\frac{1}{mn}\sum _{i,j}^{}{t}_{ij}\right)$
Analysis Under the Complicated Model Now repeat the preceding analysis with respect to a biochemical model in which each β_ij and α_ij is different. Again, take noise into account and to see what remnants of noise survive in the construction of the matched filter f as well as in the construction of {tilde over (t)}_j. Under the complicated model, from (*) above:

[0162] ${\beta }_{ij}{t}_{ij}=\log \frac{Q_T}{{\alpha }_{ij}}-\log n_j-\log Q+Z_{ij}$
Which can be rewritten as:

[0163] $\mathrm{(*}\mathrm{**)}{t}_{\mathrm{ij}}=\frac{1}{{\beta }_{\mathrm{ij}}}\log \frac{Q^T}{{\alpha }_{\mathrm{ij}}}-\frac{1}{{\beta }_{\mathrm{ij}}}\log n_j-\frac{1}{{\beta }_{\mathrm{ij}}}\log Q+Z_{\mathrm{ij}}$
The i-th component of the matched filter f is given by:

[0164] $f_i\stackrel{\Delta }{=}{\Delta }_i^Y-{\Delta }_i^X=\left\{t_{i,23}^Y-{\stackrel{\_}{t}}_7^Y\right\}-\text{{}\left(t_{i,23}^X-{\stackrel{\_}{t}}_7^X\right\}=\left\{\left(\frac{1}{{\beta }_{i,23}}\log \frac{Q_T}{{\alpha }_{i,23}}-\frac{1}{{\beta }_{i,23}}\log n_{23}^Y-\frac{1}{{\beta }_{i,23}}\log Q^Y+Z_{i,23}^Y\right)-\frac{1}{n}\sum _i^{}\left(\frac{1}{{\beta }_{i,7}}\log \frac{Q_T}{{\alpha }_{i,7}}-\frac{1}{{\beta }_{i,7}}\log n_7^Y-\frac{1}{{\beta }_{i,7}}\log Q^Y+Z_{i,7}^Y\right)\right\}-\left\{\left(\frac{1}{{\beta }_{i,23}}\log \frac{Q_T}{{\alpha }_{i,23}}-\frac{1}{{\beta }_{i,23}}\log n_{23}^X-1\frac{1}{{\beta }_{i,23}}\log Q^X+Z_{i,23}^X\right)-\frac{1}{n}\sum _i^{}\left(\frac{1}{{\beta }_{i,7}}\log \frac{Q_T}{{\alpha }_{i,7}}-\frac{1}{{\beta }_{i,7}}\log n_7^X-\frac{1}{{\beta }_{i,7}}\log Q^X+Z_{i,7}^X\right)\right\}=\frac{1}{{\beta }_{i,23}}+\left(\frac{1}{{\beta }_{i,23}}-\left(\frac{1}{n}\sum _i^{}\frac{1}{{\beta }_{i,7}}\right)\right)\log \frac{Q^Y}{Q^X}+\left\{Z_{i,23}^Y-Z_{i,23}^X+\frac{1}{n}\sum _i^{}{Z}_{i,7}^X-\frac{1}{n}\sum _i^{}{Z}_{i,7}^Y\right\}$
Under the complicated model, this gives:

[0165] ${\tilde{t}}_{ij}\stackrel{\Delta }{=}{t}_j-\stackrel{\_}{t}=\left\{\frac{1}{{\beta }_{ij}}\log \frac{Q_T}{{\alpha }_{ij}}-\frac{1}{{\beta }_{ij}}\log n_j-\frac{1}{{\beta }_{ij}}\log Q+Z_{ij}\right\}-\frac{1}{mn}\sum _{ij}^{}\left\{\frac{1}{{\beta }_{ij}}\log \frac{Q_T}{{\alpha }_{ij}}-\frac{1}{{\beta }_{ij}}\log n_j-\frac{1}{{\beta }_{ij}}\log Q+Z_{ij}\right\}$
An Alternate Way to Regularize CT Values

[0166] In another embodiment of the method, one can average the CT values rather than transforming to exponential scale and then taking logs, as this distorts the noise so that it is no longer zero mean. First, start with known Q and solve for betas. Then do multiple experiments with known n_j to solve for alphas. Since aneuploidy is a whole set of hypotheses, it is convenient to use ML to determine the most likely n_j and Q values, and then use this as a basis for calculating the most likely aneuploid state, e.g., by taking the n_j value that is most off from 1 and pushing it to its nearest aneuploid neighbor.

Estimation of the Error Rates in the Embryonic Measurements.

[0167] In one embodiment of the invention, it is possible to determine the conditional probabilities of particular embryonic measurements given specific underlying true states in embryonic DNA. In certain contexts, the given data consists of (i) the data about the parental SNP states, measured with a high degree of accuracy, and (ii) measurements on all of the SNPs in a specific blastomere, measured poorly.

[0168] Use the following notation: U—is any specific homozygote, Ū is the other homozygote at that SNP, H is the heterozygote. The goal is to determine the probabilities (p_ij) shown in Table 2. For instance p₁₁ is the probability of the embryonic DNA being U and the readout being U as well. There are three conditions that these probabilities have to satisfy:
p₁₁+p₁₂+p₁₃+p₁₄=1  (1)
p₂₁+p₂₂+p₂₃+p₂₄=1  (2)
p₂₁=p₂₃  (3)
The first two are obvious, and the third is the statement of symmetry of heterozygote dropouts (H should give the same dropout rate on average to either U or Ū).

[0169] There are 4 possible types of matings: U×U, U×Ū, UxH, HxH. Split all of the SNPs into these 4 categories depending on the specific mating type. Table 3 shows the matings, expected embryonic states, and then probabilities of specific readings (p_ij). Note that the first two rows of this table are the same as the two rows of the Table 2 and the notation (p_ij) remains the same as in Table 2.

[0170] Probabilities p_3i and p_4i can be written out in terms of p_1i and p_2i.
p₃₁=½[p₁₁+p₂₁]  (4)
p₃₂=½[p₁₂+p₂₂]  (5)
p₃₃=½[p₁₃+p₂₃]  (6)
p₃₄=½[p₁₄+p₂₄]  (7)
p₄₁=¼[p₁₁+²p₂₁+p₁₃]  (8)
p₄₂=½[p₁₂+p₂₂]  (9)
p₄₃=¼[p₁₁+2p₂₃+p₁₃]  (10)
p₄₄=½[p₁₄+p₂₄]  (11)
These can be thought of as a set of 8 linear constraints to add to the constraints (1), (2), and (3) listed above. If a vector P=[p₁₁, p₁₂, p₁₃, p₁₄, p₂₁ . . . , p₄₄]^T (16×1 dimension) is defined, then the matrix A (11×16) and a vector C can be defined such that the constraints can be represented as:
AP=C  (12)
C=[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]^T. Specifically, A is shown in Table 4, where empty cells have zeroes.

[0171] The problem can now be framed as that of finding P that would maximize the likelihood of the observations and that is subject to a set of linear constraints (AP=C). The observations come in the same 16 types as p_ij. These are shown in Table 5. The likelihood of making a set of these 16 n_ij observations is defined by a multinomial distribution with the probabilities p_ij and is proportional to:

[0172] $\begin{array}{cc}L\left(P,n_{ij}\right)\propto \prod _{ij}^{}{p}_{ij}^{n_{ij}}& \left(13\right)\end{array}$
Note that the full likelihood function contains multinomial coefficients that are not written out given that these coefficients do not depend on P and thus do not change the values within P at which L is maximized. The problem is then to find:

[0173] $\begin{array}{cc}\underset{P}{\max }\left[L\left(P,n_{ij}\right)\right]=\underset{P}{\max }\left[\ln \left(L\left(P,n_{ij}\right)\right)\right]=\underset{P}{\max }\left(\sum _{ij}^{}{n}_{ij}\ln \left(p_{ij}\right)\right)& \left(14\right)\end{array}$
subject to the constraints AP=C.
Note that in (14) taking the ln of L makes the problem more tractable (to deal with a sum instead of products). This is standard given that value of x such that f(x) is maximized is the same for which ln(f(x)) is maximized. P(n_j,Q,D)=P(n_j)P(Q)P(D_j|Q,n_j)P(D_k≠j|Q).
D MAP Detection of Aneuploidy without Parents

[0174] In one embodiment of the invention, the PS method can be applied to determine the number of copies of a given chromosome segment in a target without using parental genetic information. In this section, a maximum a-posteriori (MAP) method is described that enables the classification of genetic allele information as aneuploid or euploid. The method does not require parental data, though when parental data are available the classification power is enhanced. The method does not require regularization of channel values. One way to determine the number of copies of a chromosome segment in the genome of a target individual by incorporating the genetic data of the target individual and related individual(s) into a hypothesis, and calculating the most likely hypothesis is described here. In this description, the method will be applied to ct values from TAQMAN measurements; it should be obvious to one skilled in the art how to apply this method to any kind of measurement from any platform. The description will focus on the case in which there are measurements on just chromosomes X and 7; again, it should be obvious to one skilled in the art how to apply the method to any number of chromosomes and sections of chromosomes.

Setup of the Problem

[0175] The given measurements are from triploid blastomeres, on chromosomes X and 7, and the goal is to successfully make aneuploidy calls on these. The only “truth” known about these blastomeres is that there must be three copies of chromosome 7. The number of copies of chromosome X is not known.

[0176] The strategy here is to use MAP estimation to classify the copy number N₇ of chromosome 7 from among the choices {1,2,3} given the measurements D. Formally that looks like this:

[0177] ${\hat{n}}_7=\underset{n_7\in \left(1,2,3\right)}{\arg \max }P\left(n_7,D\right)$
Unfortunately, it is not possible to calculate this probability, because the probability depends on the unknown quantity Q. If the distribution f on Q were known, then it would be possible to solve the following:

[0178] ${\hat{n}}_7=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)P\left(n_7,D|Q\right)dQ$
In practice, a continuous distribution on Q is not known. However, identifying Q to within a power of two is sufficient, and in practice a probability mass function (pmf) on Q that is uniform on say {2¹,2² . . . , 2⁴⁰} can be used. In the development that follows, the integral sign will be used as though a probability distribution function (pdf) on Q were known, even though in practice a uniform pmf on a handful of exponential values of Q will be substituted.

[0179] This discussion will use the following notation and definitions: N₇ is the copy number of chromosome seven. It is a random variable. n₇ denotes a potential value for N₇. N_X is the copy number of chromosome X. n_X denotes a potential value for N_X. N_j is the copy number of chromosome-j, where for the purposes here jϵ{7,X}. n_j denotes a potential value for N_j. D is the set of all measurements. In one case, these are TAQMAN measurements on chromosomes X and 7, so this gives D={D₇,D_X}, where D_j={t_ij^A,t_ij^C} is the set of TAQMAN measurements on this chromosome. t_ij^A is the ct value on channel-A of locus i of chromosome-j. Similarly, t_ij^C is the value on channel-C of locus i of chromosome-j. (A is just a logical name and denotes the major allele value at the locus, while C denotes the minor allele value at the locus.) Q represents a unit-amount of genetic material such that, if the copy number of chromosome-j is n_j, then the total amount of genetic material at any locus of chromosome-j is n_jQ. For example, under trisomy, if a locus were AAC, then the amount of A-material at this locus would be 2Q, the amount of C-material at this locus is Q, and the total combined amount of genetic material at this locus is 3Q. (n^A,n^C) denotes an unordered allele patterns at a locus when the copy number for the associate chromosome is n. n^A is the number of times allele A appears on the locus and n^C is the number of times allele C appears on the locus. Each can take on values in 0, . . . , n, and it must be the case that n^A+n^C=n. For example, under trisomy, the set of allele patterns is {(0,3), (1,2), (2,1), (3,0)}. The allele pattern (2,1) for example corresponds to a locus value of A^2C, i.e., that two chromosomes have allele value A and the third has an allele value of C at the locus. Under disomy, the set of allele patterns is {(0,2), (1,1), (2,0)}. Under monosomy, the set of allele patterns is {(0,1), (1,0)}.

Q_T is the (known) threshold value from the fundamental TAQMAN equation Q₀2^βr=Q_T.

β is the (known) doubling-rate from the fundamental TAQMAN equation Q₀2^βt=Q_T.

⊥ (pronounced “bottom”) is the ct value that is interpreted as meaning “no signal”.

f_Z(χ) is the standard normal Gaussian pdf evaluated at χ.

σ is the (known) standard deviation of the noise on TAQMAN ct values.

MAP Solution

In the solution below, the following assumptions have been made:

[0180] N₇ and N_χ are independent.

[0181] Allele values on neighboring loci are independent.

[0182] The goal is to classify the copy number of a designated chromosome. In this case, the description will focus on chromosome 7. The MAP solution is given by:

[0183] ${\hat{n}}_7=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)P\left(n_7,D|Q\right)dQ=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)\sum _{n_X\in \left\{1,2,3\right\}}P\left(n_7,n_X,D|Q\right)dQ=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)\sum _{n_X\in \left\{1,2,3\right\}}P\left(n_7\right)P\left(n_X\right)P\left(D_7|Q,n_7\right)P\left(D_X|Q,n_X\right)dQ=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)\left(P\left(n_7\right)P\left(D_7|Q,n_7\right)\right)\left(\sum _{n_X\in \left\{1,2,3\right\}}P\left(n_X\right)P\left(D_X|Q,n_X\right)\right)dQ=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)\left(P\left(n_7\right)\prod _{i}P\left(t_{i,7}^A,t_{i,7}^C|Q,n_7\right)\right)\left(\sum _{n_X\in \left\{1,2,3\right\}}P\left(n_X\right)\prod _{i}P\left(t_{i,X}^A,t_{i,X}^C|Q,n_X\right)\right)dQ=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }\int f\left(Q\right)\left(P\left(n_7\right)\prod _i\sum _{n^A+n^C=n_7}P\left(n^A,n^C|n_7,i\right)P\left(t_{i,7}^A|Q,n^A\right)P\left(t_{i,7}^C|Q,n^C\right)\right)\times \left(\sum _{n_X\in \left\{1,2,3\right\}}P\left(n_X\right)\prod _i\sum _{n^A+n^C=n_X}P\left(n^A,n^C|n_X,i\right)P\left(t_{i,X}^A|Q,n^A\right)P\left(t_{i,X}^C|Q,n^C\right)\right)dQ$
Allele Distribution Model

[0184] Equation (*) depends on being able to calculate values for P(n^A,n^C ln₇,i) and P(n^A,n^C ln_X,i). These values may be calculated by assuming that the allele pattern (n^A,n^C) is drawn i.i.d (independent and identically distributed) according to the allele frequencies for its letters at locus i. An example should suffice to illustrate this. Calculate P((2,1)ln₇=3) under the assumption that the allele frequency for A is 60%, and the minor allele frequency for C is 40%. (As an aside, note that P((2,1)ln₇-2)-0, since in this case the pair must sum to 2.) This probability is given by

[0185] $P\left(\left(2,1\right)|n_7=3\right)=\left(\begin{array}{c}3\\ 2\end{array}\right){\left(\mathrm{.60}\right)}^2\left(\mathrm{.40}\right)$
The general equation is

[0186] $P\left(n^A,n^C|n_j,i\right)=\left(\begin{array}{c}n\\ n^A\end{array}\right){\left(1-p_{\mathrm{ij}}\right)}^{n^A}{\left(p_{\mathrm{ij}}\right)}^{n^C}$
Where p_i,j is the minor allele frequency at locus i of chromosome j.
Error Model
Equation (*) depends on being able to calculate values for P(t^A|Q,n^A) and P(t^C|Q,n^C). For this an error model is needed. One may use the following error model:

[0187] $P\left(t^A|Q,n^A\right)=\left\{\begin{array}{cc}{p}_d& t^A=\perp \mathrm{and}{n}^A>0\\ \left(1-p_d\right)f_Z\left(\frac{1}{\sigma }\left(t^A-\frac{1}{\beta }\log \frac{Q_T}{n^{A}Q}\right)\right)& t^A\bumpeq ̸\perp \mathrm{and}{n}^A>0\\ 1-p_a& t^A=\perp \mathrm{and}{n}^A=0\\ 2p_a{f}_Z\left(\frac{1}{\sigma }\left(t^A-\perp \right)\right)& t^A\bumpeq ̸\perp \mathrm{and}{n}^A-0\end{array}\right\}$

[0188] Each of the four cases mentioned above is described here. In the first case, no signal is received, even though there was A-material on the locus. That is a dropout, and its probability is therefore p_d. In the second case, a signal is received, as expected since there was A-material on the locus. The probability of this is the probability that a dropout does not occur, multiplied by the pdf for the distribution on the ct value when there is no dropout. (Note that, to be rigorous, one should divide through by that portion of the probability mass on the Gaussian curve that lies below ⊥, but this is practically one, and will be ignored here.) In the third case, no signal was received and there was no signal to receive. This is the probability that no drop-in occurred, 1-p_a. In the final case, a signal is received even through there was no A-material on the locus. This is the probability of a drop-in multiplied by the pdf for the distribution on the ct value when there is a drop-in. Note that the ‘2’ at the beginning of the equation occurs because the Gaussian distribution in the case of a drop-in is modeled as being centered at 1. Thus, only half of the probability mass lies below ⊥ in the case of a drop-in, and when the equation is normalized by dividing through by one-half, it is equivalent to multiplying by 2. The error model for P(t^C|Q,n^C) by symmetry is the same as for P(t^A|Q,n^A) above. It should be obvious to one skilled in the art how different error models can be applied to a range of different genotyping platforms, for example the ILLUMINA INFINIUM genotyping platform.

Computational Considerations

[0189] In one embodiment of the invention, the MAP estimation mathematics can be carried out by brute-force as specified in the final MAP equation, except for the integration over Q. Since doubling Q only results in a difference in ct value of 1/β, the equations are sensitive to Q only on the log scale. Therefore to do the integration it should be sufficient to try a handful of Q-values at different powers of two and to assume a uniform distribution on these values. For example, one could start at Q=Q_T2^−20β, which is the quantity of material that would result in a ct value of 20, and then halve it in succession twenty times, yielding a final Q value that would result in a ct value of 40.

[0190] What follows is a re-derivation of a derivation described elsewhere in this disclosure, with slightly difference emphasis, for elucidating the programming of the math. Note that the variable D below is not really a variable. It is always a constant set to the value of the data set actually in question, so it does not introduce another array dimension when representing in MATLAB. However, the variables D_j do introduce an array dimension, due to the presence of the index j.

[0191] ${\hat{n}}_7=\underset{n_7\in \left\{1,2,3\right\}}{\arg \max }P\left(n_7,D\right)$ $P\left(n_7,D\right)=\sum _{Q}P\left(n_7,Q,D\right)$ $P\left(n_7,Q,D\right)=P\left(n_7\right)P\left(Q\right)P\left(D_7|Q,n_7\right)P\left(D_X|Q\right)$ $P\left(D_j|Q\right)=\sum _{n_j\in \left\{1,2,3\right\}}P\left(D_j,n_j|Q\right)$ $P\left(D_j,n_j|D\right)=P\left(n_j\right)P\left(D_j|Q,n_j\right)$ $P\left(D_j|Q,n_j\right)=\prod _{i}P\left(D_{ij}|Q,n_j\right)$ $P\left(D_{ij}|Q,n_j\right)=\sum _{n^A+n^C=n_j}P\left(D_{ij},n^A,n^C|Q,n_j\right)$ $P\left(\begin{array}{c}{D}_{ij},n^A,\\ n^C|Q,n_j\end{array}\right)=P\left(n^A,n^C|n_j,i\right)P\left(t_{ij}^A|Q,n^A\right)P\left(t_{ij}^C|Q,n^C\right)$ $P\left(n^A,n^C|n_j,i\right)=\left(\begin{array}{c}n\\ n^A\end{array}\right){\left(1-p_{ij}\right)}^{n^A}{\left(p_{ij}\right)}^{n^C}$ $P\left(\begin{array}{c}{t}_{ij}^A|Q,\\ n^A\end{array}\right)=\left\{\begin{array}{cc}{p}_d& t^A=\perp \mathrm{and}{n}^A>0\\ \left(1-p_d\right)f_Z\left(\frac{1}{\sigma }\left(t^A-\frac{1}{\beta }\log \frac{Q_T}{n^{A}Q}\right)\right)& t^A\bumpeq ̸\perp \mathrm{and}{n}^A>0\\ 1-p_a& t^A=\perp \mathrm{and}{n}^A=0\\ 2p_d{f}_Z\left(\frac{1}{\sigma }\left(t^A-\perp \right)\right)& t^A\bumpeq ̸\perp \mathrm{and}{n}^A=0\end{array}\right\}$
E MAP Detection of Aneuploidy with Parental Info