ERROR CORRECTION IN ANCESTRY CLASSIFICATION
Error correction in ancestry classification includes obtaining, from a classifier, initial ancestry classifications associated with portions of two phased haplotypes of a chromosome pair of an individual; performing error correction on an initial ancestry classification, including detecting a phasing error in the initial ancestry classifications; and outputting a corrected ancestry classification in which the phasing error is corrected.
This application is a continuation of and claims priority to U.S. application Ser. No. 13/801,653, titled “ERROR CORRECTION IN ANCESTRY CLASSIFICATION,” filed Mar. 13, 2013, which claims priority to U.S. Provisional Patent Application No. 61/724,228, titled ANCESTRY PAINTING WITH LOCAL ANCESTRY INFERENCE, filed Nov. 8, 2012, each of which is incorporated herein by reference in its entirety for all purposes.
BACKGROUND OF THE INVENTIONAncestry deconvolution refers to identifying the ancestral origin of chromosomal segments in individuals. Ancestry deconvolution in admixed individuals (i.e., individuals whose ancestors such as grandparents are from different regions) is straightforward when the ancestral populations considered are sufficiently distinct (e.g., one grandparent is from Europe and another from Asia). To date, however, existing approaches are typically ineffective at distinguishing between closely related populations (e.g., within Europe). Moreover, due to their computational complexity, most existing methods for ancestry deconvolution are unsuitable for application in large-scale settings, where the reference panels used contain thousands of individuals.
Various embodiments of the invention are disclosed in the following detailed description and the accompanying drawings.
The invention can be implemented in numerous ways, including as a process; an apparatus; a system; a composition of matter; a computer program product embodied on a computer readable storage medium; and/or a processor, such as a processor configured to execute instructions stored on and/or provided by a memory coupled to the processor. In this specification, these implementations, or any other form that the invention may take, may be referred to as techniques. In general, the order of the steps of disclosed processes may be altered within the scope of the invention. Unless stated otherwise, a component such as a processor or a memory described as being configured to perform a task may be implemented as a general component that is temporarily configured to perform the task at a given time or a specific component that is manufactured to perform the task. As used herein, the term ‘processor’ refers to one or more devices, circuits, and/or processing cores configured to process data, such as computer program instructions.
A detailed description of one or more embodiments of the invention is provided below along with accompanying figures that illustrate the principles of the invention. The invention is described in connection with such embodiments, but the invention is not limited to any embodiment. The scope of the invention is limited only by the claims and the invention encompasses numerous alternatives, modifications and equivalents. Numerous specific details are set forth in the following description in order to provide a thorough understanding of the invention. These details are provided for the purpose of example and the invention may be practiced according to the claims without some or all of these specific details. For the purpose of clarity, technical material that is known in the technical fields related to the invention has not been described in detail so that the invention is not unnecessarily obscured.
A pipelined ancestry deconvolution process to predict an individual's ancestry based on genetic information is disclosed. Unphased genotype data associated with the individual's chromosomes is received and phased to generate phased haplotype data. In some embodiments, dynamic programming that does not require the unphased genotype data to be included in the reference data is implemented to facilitate phasing. The phased data is divided into segments, which are classified as being associated with specific ancestries. The classification is performed using a learning machine in some embodiments. The classification output undergoes an error correction process to reduce noise and correct for any phasing errors (also referred to as switch errors) and/or correlated classification errors. The error corrected output is optionally recalibrated, and ancestry labels are optionally clustered according to a geographical hierarchy to be displayed to the user.
In some embodiments, genotype data comprising gene sequences and/or genetic markers is used to represent an individual's genome. Examples of such genetic markers include Single Nucleotide Polymorphisms (SNPs), which are points along the genome, each corresponding to two or more common variations; Short Tandem Repeats (STRs), which are repeated patterns of two or more repeated nucleotide sequences adjacent to each other; and Copy-Number Variants (CNVs), which include longer sequences of deoxyribonucleic acid (DNA) that could be present in varying numbers in different individuals. Although SNP-based genotype data is described extensively below for purposes of illustration, the technique is also applicable to other types of genotype data such as STRs and CNVs. As used herein, a haplotype refers to DNA on a single chromosome of a chromosome pair. Haplotype data representing a haplotype can be expressed as a set of markers (e.g., SNPs, STRs, CNVs, etc.) or a full DNA sequence set.
Processor 102 is coupled bi-directionally with memory 110, which can include a first primary storage, typically a random access memory (RAM), and a second primary storage area, typically a read-only memory (ROM). As is well known in the art, primary storage can be used as a general storage area and as scratch-pad memory, and can also be used to store input data and processed data. Primary storage can also store programming instructions and data, in the form of data objects and text objects, in addition to other data and instructions for processes operating on processor 102. Also as is well known in the art, primary storage typically includes basic operating instructions, program code, data, and objects used by the processor 102 to perform its functions (e.g., programmed instructions). For example, memory 110 can include any suitable computer-readable storage media, described below, depending on whether, for example, data access needs to be bi-directional or uni-directional. For example, processor 102 can also directly and very rapidly retrieve and store frequently needed data in a cache memory (not shown).
A removable mass storage device 112 provides additional data storage capacity for the computer system 100, and is coupled either bi-directionally (read/write) or uni-directionally (read only) to processor 102. For example, storage 112 can also include computer-readable media such as magnetic tape, flash memory, PC-CARDS, portable mass storage devices, holographic storage devices, and other storage devices. A fixed mass storage 120 can also, for example, provide additional data storage capacity. The most common example of mass storage 120 is a hard disk drive. Mass storage 112, 120 generally store additional programming instructions, data, and the like that typically are not in active use by the processor 102. It will be appreciated that the information retained within mass storage 112 and 120 can be incorporated, if needed, in standard fashion as part of memory 110 (e.g., RAM) as virtual memory.
In addition to providing processor 102 access to storage subsystems, bus 114 can also be used to provide access to other subsystems and devices. As shown, these can include a display monitor 118, a network interface 116, a keyboard 104, and a pointing device 106, as well as an auxiliary input/output device interface, a sound card, speakers, and other subsystems as needed. For example, the pointing device 106 can be a mouse, stylus, track ball, or tablet, and is useful for interacting with a graphical user interface.
The network interface 116 allows processor 102 to be coupled to another computer, computer network, or telecommunications network using a network connection as shown. For example, through the network interface 116, the processor 102 can receive information (e.g., data objects or program instructions) from another network or output information to another network in the course of performing method/process steps. Information, often represented as a sequence of instructions to be executed on a processor, can be received from and outputted to another network. An interface card or similar device and appropriate software implemented by (e.g., executed/performed on) processor 102 can be used to connect the computer system 100 to an external network and transfer data according to standard protocols. For example, various process embodiments disclosed herein can be executed on processor 102, or can be performed across a network such as the Internet, intranet networks, or local area networks, in conjunction with a remote processor that shares a portion of the processing. Additional mass storage devices (not shown) can also be connected to processor 102 through network interface 116.
An auxiliary I/O device interface (not shown) can be used in conjunction with computer system 100. The auxiliary I/O device interface can include general and customized interfaces that allow the processor 102 to send and, more typically, receive data from other devices such as microphones, touch-sensitive displays, transducer card readers, tape readers, voice or handwriting recognizers, biometrics readers, cameras, portable mass storage devices, and other computers.
In addition, various embodiments disclosed herein further relate to computer storage products with a computer readable medium that includes program code for performing various computer-implemented operations. The computer-readable medium is any data storage device that can store data which can thereafter be read by a computer system. Examples of computer-readable media include, but are not limited to, all the media mentioned above: magnetic media such as hard disks, floppy disks, and magnetic tape; optical media such as CD-ROM disks; magneto-optical media such as optical disks; and specially configured hardware devices such as application-specific integrated circuits (ASICs), programmable logic devices (PLDs), and ROM and RAM devices. Examples of program code include both machine code, as produced, for example, by a compiler, or files containing higher level code (e.g., script) that can be executed using an interpreter.
The computer system shown in
In some embodiments, DNA samples (e.g., saliva, blood, etc.) are collected from genotyped individuals and analyzed using DNA microarray or other appropriate techniques. The genotype information is obtained (e.g., from genotyping chips directly or from genotyping services that provide assayed results) and stored in database 208 and is used by system 206 to make ancestry predictions. Reference data, including genotype data of unadmixed individuals (e.g., individuals whose ancestors came from the same region), simulated data (e.g., results of machine-based processes that simulate biological processes such as recombination of parents' DNA), pre-computed data (e.g., a precomputed reference haplotype graph used in out-of-sample phasing) and the like can also be stored in database 208 or any other appropriate storage unit.
The input to phasing module 302 comprises unphased genotype data, and the output of the phasing module comprises phased genotype data (e.g., two sets of haplotype data). In some embodiments, phasing module 302 performs out-of-sample phasing where the unphased genotype data being phased is not included in the reference data used to perform phasing. The phased genotype data is input into local classification module 304, which outputs predicted ancestry information associated with the phased genotype data. In some embodiments, the phased genotype data is segmented, and the predicted ancestry information includes one or more ancestry predictions associated with the segments. The posterior probabilities associated with the predictions are also optionally output. The predicted ancestry information is sent to error correction module 306, which averages out noise in the predicted ancestry information and corrects for phasing errors introduced by the phasing module and/or correlated prediction errors introduced by the local classification module. The output of the error correction module can be presented to the user (e.g., via an appropriate user interface). Optionally, the error correction module sends its output (e.g., error corrected posterior probabilities) to a recalibration module 308, which recalibrates the output to establish confidence levels based on the error corrected posterior probabilities. Also optionally, the calibrated confidence levels are further sent to label clustering module 310 to identify appropriate ancestry assignments that meet a confidence level requirement.
The modules described above can be implemented as software components executing on one or more processors, as hardware such as programmable logic devices and/or Application Specific Integrated Circuits designed to perform certain functions or a combination thereof. In some embodiments, the modules can be embodied by a form of software products which can be stored in a nonvolatile storage medium (such as optical disk, flash storage device, mobile hard disk, etc.), including a number of instructions for making a computer device (such as personal computers, servers, network equipment, etc.) implement the methods described in the embodiments of the present application. The modules may be implemented on a single device or distributed across multiple devices. The functions of the modules may be merged into one another or further split into multiple sub-modules.
In addition to being a part of the pipelined ancestry prediction process, the modules and their outputs can be used in other applications. For example, the output of the phasing module can be used to identify familial relatives of individuals in the reference database.
Details of the modules and their operations are described below.
PhasingAt a given gene locus on a pair of autosomal chromosomes, a diploid organism (e.g., a human being) inherits one allele of the gene from the mother and another allele of the gene from the father. At a heterozygous gene locus, two parents contribute different alleles (e.g., one A and one C). Without additional processing, it is impossible to tell which parent contributed which allele. Such genotype data that is not attributed to a particular parent is referred to as unphased genotype data. Typically, initial genotype readings obtained from genotyping chips manufactured by companies such as Illumina® are in an unphased form.
Phasing is often done using statistical techniques. Such techniques are also referred to as population-based phasing because genotype data from a reference collection of a population of individuals (e.g., a few hundred to a thousand) is analyzed. BEAGLE is a commonly used population-based phasing technique. It makes statistical determinations based on the assumption that certain blocks of haplotypes are inherited in blocks and therefore shared amongst individuals. For example, if the genotype data of a sample population comprising many individuals shows a common pattern of “?A ?C ?G ?T” (where “?” can be any other allele), then the block “ACGT” is likely to be a common block of haplotypes that is present in these individuals. The population-based phasing technique would therefore identify the block “ACGT” as coming from one parent whenever “?A ?C ?G ?T” is present in the genotype data. Because BEAGLE requires that the genotype data being analyzed be included in the reference collection, the technique is referred to as in-sample phasing.
In-sample phasing is often computationally inefficient. Phasing of a large database of a user's genome (e.g., 100,000 or more) can take many days, and it can take just as long whenever a new user has to be added to the database since the technique would recompute the full set of data (including the new user's data). There can also be mistakes during in-sample phasing. One type of mistake, referred to as phasing errors or switch errors, occurs where a section of the chromosome is in fact attributed to one parent but is misidentified as attributed to another parent. Switch errors can occur when a stretch of genotype data is not common in the reference population. For example, suppose that a parent actually contributed the haplotype of “ACCC” and another parent actually contributed the haplotype of “AAGT” to genotype 502. Because the block “ACGT” is common in the reference collection and “ACCC” has never appeared in the reference collection, the technique attributes “ACGT” and “AACC” to two parents respectively, resulting in a switch error.
Embodiments of the phasing technique described below permit out-of-sample population-based phasing. In out-of-sample phasing, when genotype data of a new individual needs to be phased, the genotype data is not necessarily immediately combined with the reference collection to obtain phasing for this individual. Instead, a precomputed data structure such as a predetermined reference haplotype graph is used to facilitate a dynamic programming based process that quickly phases the genotype data. For example, given the haplotype graph and unphased data, the likely sequence of genotype data can be solved using the Viterbi algorithm. This way, on a platform with a large number of users forming a large reference collection (e.g., at least 100,000 individuals), when a new individual signs up with the service and provides his/her genotype data, the platform is able to quickly phase the genotype data without having to recompute the common haplotypes of the existing users plus the new individual.
At 602, unphased genotype data of the individual is obtained. In some embodiments, the unphased genotype data such as sequence data 502 is received from a database, a genotyping service, or as an upload by a user of a platform such as 100.
At 604, the unphased genotype data is processed using dynamic programming to determine phased data, i.e., sets of likely haplotypes. The processing requires a reference population and is therefore referred to as population-based phasing. In some embodiments, the dynamic programming relies on a predetermined reference haplotype graph. The predetermined haplotype graph is precomputed without referencing the unphased genotype data of the individual. Thus, the unphased genotype data is said to be out-of-sample with respect to a collection of reference genotype data used to compute the predetermined reference haplotype graph. In other words, if the unphased genotype data is from a new user whose genotype data is not already included in the reference genotype data and therefore is not incorporated into the predetermined reference haplotype graph, it is not necessary to include the unphased genotype data from the new user in the reference genotype data and recompute the reference haplotype graph. Details of dynamic programming and the predetermined reference haplotype graph are described below.
At 606, trio-based phasing is optionally performed to improve upon the results from population-based phasing. As used herein, trio-based phasing refers to phasing by accounting for the genotyping data of one or more biological parents of the individual.
At 608, the likely haplotype data is output to be stored to a database and/or processed further. In some embodiments, the likely haplotype data is further processed by a local classifier as shown in
The likely haplotype data can also be used in other applications, such as being compared with haplotype data of other individuals in a database to identify the amount of DNA shared among individuals, thereby determining people who are related to each other and/or people belonging to the same population groups.
In some embodiments, the dynamic programming process performed in step 604 uses a predetermined reference haplotype graph to examine possible sequences of haplotypes that could be combined to generate the unphased genotype data, and determine the most likely sequences of haplotypes. Given a collection of binary strings of length L, a haplotype graph is a probabilistic deterministic finite automaton (DFA) defined over a directed acyclic graph. The nodes of the multigraph are organized into L+1 levels (numbered from 0 to L), such that level 0 has a single node representing the source (i.e., initial state) of the DFA and level L has a single node representing the sink (i.e., accepting state) of the DFA. Every directed edge in the multigraph connects a node from some level i to a node in level (i+1) and is labeled with either 0 or 1. Every node is reachable from the source and has a directed path to the sink. For each path through the haplotype graph from the source to the sink, the concatenation of the labels on the edges traversed by the path is a binary string of length L. Semantically, paths through the graph represent haplotypes over a genomic region comprising L biallelic markers (assuming an arbitrary binary encoding of the alleles at each site). A probability distribution over the set of haplotypes included in a haplotype graph can be defined by associating a conditional probability with each edge (such that the sum of the probabilities of the outgoing edges for each node is equal to 1), and generated by starting from the initial state at level 0, and choosing successor states by following random outgoing edges according to their assigned conditional probabilities.
Each edge is associated with a probability computed based on the reference collection of genotype data. In this example, a collection of genotype data is comprised of genotype data from 1000 individuals, of which 400 have the “A” allele at the first locus, and 600 have the “G” allele at the first locus. Accordingly, the probability associated with edge 708 is 400/1000 and the probability associated with edge 710 is 600/1000. All of the first 400 individuals have the “C” allele at the second locus, giving edge 712 a probability of 400/400. All of the next 600 individuals who had the “G” allele at the first locus have the “T” allele at the second locus, giving edge 714 a probability of 600/600, and so on. The probabilities associated with the respective edges are labeled in the diagram. The probability associated with a specific path is expressed as the product of the probabilities associated with the edges included in the path. For example, the probability associated with path 702 is computed as:
The dynamic programming process searches the haplotype graph for possible paths, selecting two paths h1 and h2 for which the product of their associated probabilities is maximized, subject to the constraint that when the two paths are combined, the alleles at each locus must match the corresponding alleles in the unphased genotype data (g). The following expression is used in some cases to characterize the process:
maximize P(h1)P(h2), subject to h1+h2=g
For out-of-sample phasing, the reference haplotype graph is built once and reused to identify possible haplotype paths that correspond to the unphased genotype data of a new individual (a process also referred to as “threading” the new individual's haplotype along the graph). The individual's genotype data sometimes does not correspond to any existing path in the graph (e.g., the individual has genotype sequences that are unique and not included in the reference population), and therefore cannot be successfully threaded based on existing paths of the reference haplotype graph. To cope with the possibility of a non-existent path, several modifications are made to the reference haplotype graph to facilitate the out-of-sample phasing process.
Recombination is one reason to extend graph 700 for out-of-sample phasing. As used herein, recombination refers to the switching of a haplotype along one path to a different path. Recombination can happen when segments of parental chromosomes cross over during meiosis. In some embodiments, reference haplotype graph 700 is extended to account for the possibility of recombination/path switching. Recombination events are modeled by allowing a new haplotype state to be selected (independent of the previous haplotype state) with probability τ at each level of the haplotype graph. By default, τ≈0.00448, which is an estimate of the probability of recombination between adjacent sites, assuming 500,000 uniformly spaced markers, a genome length of 37.5 Morgans, and 30 generations since admixture. Referring to the example of
Genotyping error is another reason to extend graph 700 for out-of-sample phasing. Genotyping errors can occur because the genotyping technology is imperfect and can make false readings. The rate of genotyping error for a given technology (e.g., a particular genotyping chip) can be obtained from the manufacturer. In some embodiments, when the search for possible paths for a new individual cannot be done according to the existing reference graph, the existing reference haplotype graph is extended to account for the possibility of genotyping errors. For example, suppose the new individual's unphased genotype data is “AG, CT, GG, CT, GG, CG,” (SEQ ID NO: 2) which cannot be split into two haplotypes by threading along existing paths in graph 700. Referring to
The example graphs shown include a small number of nodes and edges, and thus represent short sequences of genotype data. In practice, the begin state node corresponds to the first locus on the chromosome and the accepting state node the last locus on the chromosome, and the number of edges in a path corresponds to the number of SNPs in a chromosome (L), which can be on the order of 50,000 in some embodiments. The thickest portion of the graph (i.e., a locus with the greatest number of possible paths), which depends at least in part on the DNA sequences of individuals used to construct the graph (K), can be on the order of 5,000 in some embodiments. A large number of computations would be needed (O(LK4) in the worst case) for a naïve implementation of a dynamic programming solution based on the Viterbi algorithm.
In some embodiments, the paths are pruned at each state of the graph to further improve performance. In other words, only likely paths are kept in the modified graph and unlikely paths are discarded. In some embodiments, after i markers (e.g., 3 markers), paths with probabilities below a certain threshold ϵ (e.g., less than 0.0001%) are discarded. For example, a haplotype along a new path that accounts for both recombination and switching error would have very low probability of being formed, and thus can be discarded. As another example, in the case of unphased genotype data of “AG, CT, GG, CT, GG, CG,” (SEQ ID NO: 2) a new haplotype accounting for recombination can be forged by switching paths several times along the graph (additional edges would need to be added but are not shown in the diagram). Given the low probability associated with each switch, however, the formation of such a haplotype is very unlikely and would be pruned from the resulting graph, while the path that includes the genotyping error 825 has sufficiently high probability, and is kept in the graph and used to thread the unphased genotype data into phased genotype data. By pruning unlikely paths from the modified graph, the dynamic programming-based phasing process is prevented from exploring very unlikely paths in the graph when threading a new haplotype along it. The choice of e determines the trade-off between the efficiency of the algorithm (in both time and space) and the risk of prematurely excluding the best Viterbi path. Computation savings provided by pruning can be significant. In some cases, phasing using a naïve implementation can require 15 days per person while phasing with pruning only requires several minutes per person.
In some embodiments, the nodes and edges of the haplography can be represented as follows:
Even with a pruned haplotype graph, the number of nodes and edges can be large and using the above data structures to represent the graph would require a vast amount of memory (on the order of several gigabytes in some cases). In some embodiments, the graph is represented in a compressed form, using segments. The term “segment” used herein refers to the data structure used to represent the graph in a compressed form and is different from the DNA segments used elsewhere in the specification. Each segment corresponds to a contiguous set of edges in the graph, with the following constraints: the end of the segment has up to 1 branch (0 branches are permitted), and no segment points to the middle of another segment. In some embodiments, the data structure of a segment is represented as follows:
On a system such as the personal genomics services platform provided by 23andMe®, DNA sequence information of one or both parents of the individual is sometimes available and can be used to further refine phasing. With the exception of sites where all three individuals are heterozygous, the parental origin of each allele can be determined unambiguously. For ambiguous sites, knowledge of patterns of local linkage disequilibrium can be used to statistically estimate the most likely phase. In some embodiments, a refinement process that accounts for parental DNA sequence information, referred to as trio-based phasing, is optionally performed following the population-based phasing process to correct any errors in the output of the population-based phasing process and improve phasing accuracy. In some embodiments, the trio-based phasing technique is a post-processing step to be applied to sequences for which a previous population-based linkage-disequilibrium phasing approach has already been applied. The trio-based phasing technique can be used in combination with any existing phasing process to improve phasing quality, provided that an estimate of the switch error rate (also referred to as the phasing error rate) is available.
In some embodiments, trio-based phasing receives as inputs a set of preliminary phased haplotype data (e.g., output of an out-of-sample population-based phasing technique described above), and employs a probabilistic graphic model (also referred to as a dynamic Bayesian network) that models the observed alleles, hidden states, and relationships of the parental and child haplotypes. The input includes the set of preliminary phased haplotype data as well as the phased haplotype data of at least one parent. The genotype data at a particular site (e.g., the i-th SNP on a chromosome) for each individual in the trio (i.e., mom, dad, or child (i.e. the individual whose genetic data is being phased)) are represented by the following variables:
G0*,i, G1*,i∈{0, 1}: the observed alleles for haplotypes 0 and 1, provided as input data. For the child, the input data can be obtained from the output of the population-based phasing process (e.g., the preliminary haplotype data). For the parent, the input data can be the output of the population-based phasing process or the final output of a refined process.
Hm*,i, Hp*,i∈{0, 1}: the hidden true alleles of the individual's maternal (m) and paternal (p) haplotypes.
P*,i ∈{m, p}: a hidden binary phase indicator variable that is set to m whenever G0*,i corresponds to Hm*,i and set to p whenever G0*,i corresponds to Hp*,i.
The relationship between parental and child haplotypes are encoded by two additional variables, Tmom,i, Tdad,i∈{a, b}, where a indicates transmission of the parent's maternal haplotype to the child and b indicates transmission of the parent's paternal haplotype to the child. In some embodiments, a=0 and b=1.
The following assumptions are made about the model:
1. The hidden true alleles for each parent at each position (i.e., H*(mom,dad),i), the initial phase for each individual (i.e., P*,1), and the initial transmission for each parent (i.e., T*,1) are independently drawn from uniform Bernoulli priors.
2. The phase indicator variables for each individual and the transmission indicator variables for each parent are each sampled according to independent first order Markov processes. Specifically,
where s is the estimated switch error probability between consecutive sites in the input haplotypes and r is the estimated recombination probability between sites in a single meiosis. In some embodiments, s is set to a default value of 0.02 and r is set to a default value of
3. The hidden true alleles for the child at each position (i.e., H*kid,i) are deterministically set on the parents' true hidden haplotypes (i.e., neglecting the possibility of private mutations) and their respective transmission variables.
4. The observed alleles are sampled conditionally on the true alleles and the phase variables with genotyping error, according to the following model:
according to the estimated genotyping error rate.
The following expression is used to characterize the trio-based phasing process:
maximize Pr(Hmkid,Hpkid,Hmmom,Hpmom,Hmdad,Hpdad)
given Hmi+Hpi=G0i+G1i∀i∈{kid,mom,dad}
Trio-based phasing includes using the probabilistic model to estimate the most probable setting of all unobserved variables, conditioned on the observed alleles. In some embodiments, the most probable H variables are determined using a standard dynamic programming-based technique (e.g., Viterbi). One can visualize the model as plates corresponding to i∈{1, 2, . . . , L} being stacked in sequential order, and the paths are formed by the interconnections of nodes on the same plate, as well as nodes across plates.
At 1102, emission probabilities are precomputed for each plate of model 1000. In some embodiments, the emission probabilities, which correspond to the most likely setting for the H variables given the G, P, and T variables, are found using a dynamic programming (e.g., Viterbi) based process. Referring to
At 1104, transition probabilities are computed based at least in part on the values of transition probabilities from the previous position. Referring to
At 1106, based on the computed probabilities, the settings of transition variables T and P across the entire chromosome sequence (i.e., for i=1, . . . , L) are searched to determine the settings that would most likely result in the observed values. In some embodiments, the determination is made using a dynamic programming technique such as Viterbi, and 25*25*L states are searched.
At 1108, the setting of H variables is looked up across the entire sequence to determine the settings that would most likely result in the given G, P, and T variables. This requires L table lookups.
The trio-based phasing solves the most likely settings for the H variables (the hidden true alleles for the individual's maternal and paternal haplotypes at a given location). The solution is useful for phasing the child's DNA sequence information as well as for phasing a parent's DNA sequence information (if the parent's DNA sequence information is unphased initially). In the event that only one parent's DNA sequence information is available, the other parent's DNA sequence information can be partially determined based on the DNA sequence information of the known parent and the child (e.g., if the child's alleles at a particular location is “AC” and the mother's alleles at the same location are “CC”, then one of the father's alleles would be “A” and the other one is unknown). The partial information can be marked (e.g., represented using a special notation) and input to the model. The quality of trio-based phasing based on only one parent's information is still higher than population-based phasing without using the trio-based method.
In addition to improved haplotypes data, the result of trio-based phasing also indicates whether a specific allele is deemed to be inherited from the mother or the father. This information is stored and can be presented to the user in some embodiments.
Local ClassificationLocal classification refers to the classification of DNA segments as originating from an ancestry associated with a specific geographical region (e.g., Eastern Asia, Scandinavia, etc.) or ethnicity (e.g., Ashkenazi Jew).
Local classification is based on the premise that, T generations ago, all the ancestors of an individual were unadmixed (i.e., originating from the same geographical region). Starting at T generation, ancestors from different geographical regions produced admixed offspring. Genetic recombination breaks chromosomes and recombines them at each generation. After T generations, 2T meiosis occurred. As a result, the expected length of a recombination-free segment is expressed as:
where F corresponds to a fixed length segment. In some embodiments, the expected length L is determined to be 100 SNPs. It is used as the segment size (also referred to as the window size) used in local classification.
Initially, at 1202, a set of K local ancestries is obtained. In some embodiments, the specification of the local ancestries depends on the ancestries of unadmixed individuals whose DNA sequence information is used as reference data. For example, the set of local ancestries can be pre-specified to include the following: African, Native American, Ashkenazi, Eastern Asian, Southern Asian, Balkan, Eastern European, Western European, Middle Eastern, British Isles, Scandinavian, Finnish, Oceanian, Iberian, Greek, Sardinian, Italian, and Arabic. Many other specifications are possible; for example, in some embodiments the set of local ancestries correspond to individual countries such as the UK, Ireland, France, Germany, Finland, China, India, etc.
At 1204, a classifier is trained using reference data. In this example, the reference data includes DNA sequence information of unadmixed individuals, such as individuals who are self-identified or identified by the system as having four grandparents of the same ancestry (i.e., from the same region), DNA sequence information obtained from public databases such as 1KG, HGDP-CEPH, HapMap, etc. The DNA sequence information and their corresponding ancestry origins are input into the classifier, which learns the corresponding relationships between the DNA sequence information (e.g., DNA sequence segments) and the corresponding ancestry origins. In some embodiments, the classifier is implemented using a known machine learning technique such as a support vector machine (SVM), a neural network, etc. A SVM-based implementation is discussed below for purposes of illustration.
At 1206, phased DNA sequence information of a chromosome of the individual is divided into segments. In some embodiments, phased data is obtained using the improved phasing technique described above. Phased data can also be obtained using other phasing techniques such as BEAGLE. The length of the segments can be a predetermined fixed value, such as 100 SNPs. It is assumed that each segment corresponds to a single ancestry.
At 1208, the DNA sequence segments are input into the trained classifier to obtain corresponding predicted ancestries. In some embodiments, the classifier determines probabilities associated with the set of local ancestries (i.e., how likely a segment is from a particular local ancestry), and the ancestry associated with the highest probability is selected as the predicted ancestry for a particular segment.
In some embodiments, one or more SVMs are used to implement the classifier. An SVM is a known type of non-probabilistic binary classifier. It constructs a hyper plane that maximizes the distance to the closest training data point of each class (in this case, a class corresponds to a specific ancestry). A SVM can be expressed using the following general expression:
where w is the normal vector to the hyper plane, C is a penalty term (fixed), the ξ are slack variables, xi represents the features of the data point i to be classified, and yi is the class of data point i.
Since a SVM is a binary classifier and there are K (e.g., 18) classes of local ancestries to be classified, the classification can be decomposed into a set of binary problems (e.g., should the sequences be classified as African or Native American, African or Ashkenazi, Native American or Ashkenazi, etc.). One approach is the “one vs. one” technique where a total of
classifiers are trained and combined to form a single local ancestry classifier. Specifically, there is one classifier configured to determine the likelihood that a sequence is African or Native American, another to determine African or Ashkenazi, another to determine Native American or Ashkenazi, etc. During the training process, reference data of DNA sequences and their corresponding ancestries is fed to the SVM for machine learning. When an ancestry prediction for a DNA sequence segment is to be made, each trained SVM makes a determination about which one of the ancestry pair the DNA sequence segment more likely corresponds to, and the results are combined to determine which ancestry is most likely. Specifically, the ancestry that wins the highest number of determinations is chosen as the predicted ancestry. Another approach is the “one vs. all” technique where K classifiers are trained.
Several refinements can be made to improve the SVM. For example, the number of unadmixed reference individuals can vary greatly per ancestral origin. If 700 samples are from Western Europe but only 200 samples are from South Asia, the imbalance in the number of samples can cause the Western European-South Asian SVM to “favor” the larger class. Thus, the larger class is penalized to compensate for the imbalance according to the following:
where w is the normal vector to the hyper plane, CG is a penalty term for class G, the ξ are slack variables, x; represents the features of the data point i to be classified, and yi, is the class of data point i.
Another refinement is to encode strings of SNPs according to the presence or absence of features. One approach is to encode one feature at each SNP according to the presence or absence of the minor allele. Another approach is to take substrings of length 2 which have 4 features per position and which can be encoded based on their presence or absence as 00, 01, 10, and 11. A more general approach is to use a window of length L, and encode (L−k+1)·2k features of length k according to the presence or absence of the features.
The general approach is not always feasible for practical implementation, given that there are
features in a window of length L. With L=100, this number is approximately 1030, too large for most memory systems. Thus, in some embodiments, a modified kernel is used. In some embodiments, a specialized string kernel is used that computes the similarity between any two given windows as the total number of substrings they share. This approach takes into account that even very similar windows contain sites that have mutated, resulting in common subsequences along with deleted, inserted, or replaced symbols. Therefore, the specialized string kernel is a more relevant way of comparing the similarity between two 100 SNP windows, and achieves much higher accuracy than the standard linear kernel.
Another refinement is to use supervised learning. Supervised learning refers to the task of training (or learning) a classifier using a pre-labeled data, also referred to as the training set. Specifically, an SVM classifier is trained (or learned) using a training set of customers whose ancestry was known (e.g., self-reported ancestries). Parameters of the SVM classifier are adjusted during the process. The trained classifier is then used to predict a label (ancestry) for any new unlabeled data.
Error CorrectionThe results of the local classifier can contain errors.
In some embodiments, error correction is implemented using a Hidden Markov Model (HMM), which is a statistical model in which the input data is being modeled as a Markov process with unobserved (hidden) states. In an HMM, the observed signal (the input data) is being generated by a hidden process in a sequential manner. A standard HMM assumes that an observation, given the hidden state that generated it, is independent of all previous observations. The hidden state at any given position only depends on the hidden state at the previous position. In some embodiments, the input (observed data) to the HMM includes the predicted ancestries of DNA sequence segments (e.g., the ancestries as predicted by the local classifier for segments that are 100 SNPs in length). The hidden state corresponds to the true ancestries of the segments. The output of the HMM forms a set of smoothed ancestry origins for the segments.
The basic model averages out the output of the learning machine and generates a smoother and less noisy output, but does not correct many of the errors in the output. For example, in
Also, in
The graph defines a probabilistic model as follows:
Pr(H1,H2,H3, . . . ,O1,O2,O3, . . . )=Pr(H1)Pr(O1|H1)Pr(H2|H1)Pr(O2|H2,O1)Pr(H3|H2)Pr(O3|H3,O2) . . .
where probabilities P(Oj|Hj, Qj-1) are referred to as the emission parameters, and probabilities Pr(Hi|Hi-1) are referred to as the transition parameters.
The model outputs probabilities associated with ancestry assignments of the most probable sequence. Training is required to estimate the emission parameters and the transition parameters. In some embodiments, an expectation maximization method is used to estimate the parameters.
The emission parameters characterize how well the local classifier predicts the ancestry. Specifically, given the underlying true state of a segment, what is the probability that the local classifier will output the true state.
The transition parameters correspond to the probability of a particular hidden state of the model given the previous hidden state. They represent the statistical likelihood of observing certain sequences of true ancestries in the population, and therefore need to be determined based on admixed data. However, it is not possible to obtain fully transitioned and accurately labeled genomes from actual admixed individuals. Thus, to determine the transition parameters, an iterative approach is used. Initially, the transition parameters are arbitrarily chosen to establish an initial model. The initial model is used to perform error correction. Based on the error corrected results, the model is updated by applying an expectation maximization method. The process can be repeated until a final convergent model is achieved.
Once the emission parameters and the transition parameters are established, the model is fully specified. Thus, the most likely sequence of hidden variables can be determined based on the observed states using conventional HMM techniques. For example, a probabilistic scoring scheme is used to determine the most likely sequence in some embodiments. All the possibilities associated with the hidden states are listed, and a set of scoring rules are specified to reward or penalize certain events by adding or subtracting a score associated with a sequence. For example, a change in adjacent haplotypes is likely an error; therefore, whenever two adjacent haplotypes are different, the score is reduced according to the rules. A mismatched observed state/hidden state pair also indicates likely error; therefore, whenever there is a mismatch of predicted ancestry and the underlying ancestry, the score is reduced. The most likely sequence of hidden states can be determined by computing scores associated with all possible combinations of observed states and hidden states, and selecting the sequence that leads to the highest score. In some embodiments, more efficient techniques for selecting the most likely sequence such as dynamic programming are employed to break the problem down into subproblems and reduce the amount of computation required. For example, the problem can be reduced to recursively determine the best ancestry assignment for everything to the left or the right of a particular position.
As described above, training is required to obtain parameters for the PHMM (or APHMM). In some embodiments, an ensemble technique is used where reference population is grouped into distinct subsets to serve as different types of training data resulting in different types of models. For example, different types of reference individuals that tend to have similar ancestries are identified and grouped into subsets. Such subsets can be formed from admixed individuals (e.g., Latinos, Eurasians, etc.). as well as unadmixed individuals (e.g., East Asians, Northern Europeans, etc.) Data from a subset is used to determine the parameters of the model for that subset. The resulting model is a model specific to the reference group (e.g., a Latino-specific model, a Eurasian-specific model, an East Asian specific-model, a Northern European-specific model, etc.). In some embodiments, the error correction process applies its input to all available models, and the results are weighted based on confidence measures and then combined according to a Bayesian model averaging technique.
The error correction module outputs the most probable sequence of ancestry assignments for a pair of haplotypes, and posterior probabilities associated with the corresponding assignments. The posterior probabilities are recalibrated to establish confidence measures associated with the ancestry assignments. A well calibrated prediction with a probability of P should be correct P of the times. How well the posterior probability of the output is calibrated can be determined based on reference data of actual unadmixed individuals and/or simulated admixed individuals. For example, if it is known that in the reference data, 10% of the haplotype segments correspond to East European ancestry, but the output predicts with 80% posterior probability that 20% of all the haplotype segments correspond to East European ancestry, then the posterior probability is overly confident. By tallying the percentage of the reference data that corresponds to a specific ancestry, and applying the reference data to the predictive engine to obtain the posterior probability, a reliability plot of accuracy vs. posterior probability can be determined for each reference population corresponding to a specific ancestry.
In some embodiments, Platt's recalibration technique is used to recalibrate the posterior probabilities. Logistics regression is applied to posterior probabilities. A feature matrix X (e.g., 2nd degree polynomials) is defined, and a fit is determined based on the following:
K-class recalibration is then performed (K being the number of local ancestries). In some embodiments, K logistic curves are fit and renormalized. In some embodiments, K logistic curves are fit and multinomial logistic regression (i.e., softmax) is performed according to the following:
In some embodiments, an isotonic regression technique (e.g., the Zadrozny and Elkan method) is used to recalibrate the posterior probabilities, where recalibrated probabilities are estimated as percentages of well classified training examples falling in each bin.
Given the input of (yi, pi)i=1, . . . , n, the input is sorted in increasing order of pi. ϕi that monotonically increases with pi but close to yi are found. In some embodiments, a pool-adjacent-violators (PAV) algorithm is used to solve:
where yi is the label predicted for individual i, pi is the uncalibrated probability associated with the prediction and ϕi is the recalibrated probability.
In some embodiments, modified isotonic regression techniques are used. For example, pi can be bracketed into bins, and weights proportional to the bin sizes are introduced to reduce computational cost. As another example, regularization terms can be introduced to ensure smoothness of the output curves as follows:
In some embodiments, separate calibration regimes are used for individuals with different amounts of effective switch error. Specifically, separate calibration curves are fitted for unadmixed individuals (who have a low rate of effective switching error) or admixed individuals (who have a high rate of effective switching error).
Label ClusteringIn some embodiments, the recalibrated results are required to meet a threshold level of confidence before they are presented to the user. If the threshold level is unmet, the assignments are clustered and repeated as necessary until a total confidence level meets the threshold level.
The output of the label cluster outputs the predicted ancestry for each haplotype segment. In some embodiments, the information is stored in a database and/or sent to an application to be displayed.
Display of Ancestry Composition InformationIn some embodiments, once the ancestries associated with the individual's chromosomes are determined, the results can be presented via various user interfaces upon user requests. The user interfaces can also present ancestry information obtained using other techniques so long as the data being presented includes requisite information such as the specific ancestries and proportions of the individual's genotype data that corresponds to the specific ancestries.
At 2402, a request to display ancestry composition of an individual is received. In some embodiments, the request is received from an application that allows users to interact with their genetic information and/or ancestry data. Referring to
Returning to
At 2406, the ancestry composition information is presented to be displayed via a user interface.
In some embodiments, the ancestry composition information is initially displayed according to geographical regions and proportions of ancestries deemed to correspond to those geographical regions. Subsequently, the user can request different data to be displayed via user interfaces provided by the application. A further user request is optionally received at 2408. At 2409, the type of request is determined. In response to the further user request and depending on the type of request, different information is displayed. As shown, if the user request is a request to display subregions of a specific ancestry, subregions and proportions of the individual's ancestries corresponding to the subregions are displayed (or caused to be displayed on a display device by processors) at 2410 in response. If the user request is a request to display ancestries inherited from one or more parents, such information is displayed (or caused to be displayed) at 2412 if available. If the user request is a request to display ancestry composition information for a specific chromosome, the proportions of ancestries associated with the specific chromosome is displayed (or caused to be displayed) at 2414. Other types of requests/displays are possible.
The user is provided with the ability to expand the regions and view more detailed information pertaining to subregions.
The subregions can be further expanded.
In some embodiments, the ancestry composition information that is obtained includes inheritance information, including proportions of the individual's ancestries that are deemed to be inherited from either the father or the mother. In other words, the inheritance information pertains to how much of the individual's DNA corresponding to a specific ancestry is inherited from each parent. For example, the trio-based phasing result can indicate that for chromosome 1, haplotype 0, segments 1-20 correspond to Scandinavian ancestry and are inherited from the mother, segments 21-45 correspond to Italian ancestry and are inherited from the mother also, segments 46-73 correspond to Greek ancestry and are inherited from the father, and so on. The segments from either the mother or the father and the corresponding ancestries of the segments are tallied, and proportions of the ancestries attributed to each parent are computed. The inheritance information computation can be done following trio-based phasing, at the time the request to display inheritance from parents is made, or at some other appropriate time. Ancestry composition information of how much of the individual's DNA corresponding to a specific ancestry is inherited from each parent is displayed.
In some embodiments, since the ancestry deconvolution process is applied to individual chromosomes and the results are stored on the basis of individual chromosomes, the user has the option to select a specific autosomal chromosome or an X-chromosome to view its ancestral composition.
A pipelined ancestry deconvolution process and display of results have been described. The accuracy of ancestry predictions is greatly improved over existing techniques, and the results are presented in an informative and user-friendly fashion.
Although the foregoing embodiments have been described in some detail for purposes of clarity of understanding, the invention is not limited to the details provided. There are many alternative ways of implementing the invention. The disclosed embodiments are illustrative and not restrictive.
Claims
1. A method of error correction in ancestry classification, comprising:
- obtaining, from a classifier, initial ancestry classifications associated with portions of two phased haplotypes of a chromosome pair of an individual;
- performing error correction on an initial ancestry classification using one or more computer processors, including detecting a phasing error in the initial ancestry classifications; and
- outputting a corrected ancestry classification in which the phasing error is corrected.
2. The method of claim 1, wherein performing error correction on the initial ancestry classification further includes detecting a correlated prediction error that was caused by the classifier.
3. The method of claim 1, wherein performing error correction on the initial ancestry classification further includes detecting a correlated prediction error that was caused by the classifier, and the correlated prediction error is corrected in the corrected ancestry classification.
4. The method of claim 1, wherein performing error correction on the initial ancestry classification includes obtaining a Pair Hidden Markov Model (PHMM) in which an observed state corresponds to an initial ancestry classification associated with a portion of one of the two haplotypes, and a hidden state corresponds to ancestries associated with a portion of the two haplotypes.
5. The method of claim 4, wherein performing error correction on the initial ancestry classification further includes determining a likely sequence of hidden states given the initial ancestry classifications.
6. The method of claim 5, wherein the determination includes performing dynamic programming based on the PHMM.
7. The method of claim 4, wherein the PHMM is an Autoregressive Pair Hidden Markov Model (APHMM) in which the observed state is dependent on its corresponding hidden state and on a previous observed state.
8. The method of claim 1, wherein performing error correction on the initial ancestry classification includes obtaining a plurality of Pair Hidden Markov Models (PHMMs), and wherein each PHMM corresponds to a distinct reference population.
9. The method of claim 8, wherein performing error correction on the initial ancestry classification further includes applying the initial ancestry classifications associated with portions of two phased haplotypes to the plurality of PHMMs, and performing a Bayesian model averaging on outputs of the plurality of PHMMs.
10. A system of error correction in ancestry classification, comprising:
- one or more processors configured to: obtain, from a classifier, initial ancestry classifications associated with portions of two phased haplotypes of a chromosome pair of an individual; perform error correction on an initial ancestry classification, including to detect a phasing error in the initial ancestry classifications; and output a corrected ancestry classification in which the phasing error is corrected;
- one or more memories coupled to the one or more processors, configured to provide the one or more processors with instructions.
11. The system of claim 10, wherein to perform error correction on the initial ancestry classification further includes to detect a correlated prediction error that was caused by the classifier.
12. The system of claim 10, wherein to perform error correction on the initial ancestry classification further includes to detect a correlated prediction error that was caused by the classifier, and the correlated prediction error is corrected in the corrected ancestry classification.
13. The system of claim 10, wherein to perform error correction on the initial ancestry classification includes to obtain a Pair Hidden Markov Model (PHMM) in which an observed state corresponds to an initial ancestry classification associated with a portion of one of the two haplotypes, and a hidden state corresponds to ancestries associated with a portion of the two haplotypes.
14. The system of claim 13, wherein to perform error correction on the initial ancestry classification further includes to determine a likely sequence of hidden states given the initial ancestry classifications.
15. The system of claim 14, wherein the determination includes performing dynamic programming based on the PHMM.
16. The system of claim 13, wherein the PHMM is an Autoregressive Pair Hidden Markov Model (APHMM) in which the observed state is dependent on its corresponding hidden state and on a previous observed state.
17. The system of claim 10, wherein to perform error correction on the initial ancestry classification includes to obtain a plurality of Pair Hidden Markov Models (PHMMs), and wherein each PHMM corresponds to a distinct reference population.
18. The system of claim 17, wherein to perform error correction on the initial ancestry classification further includes to apply the initial ancestry classifications associated with portions of two phased haplotypes to the plurality of PHMMs, and to perform a Bayesian model averaging on outputs of the plurality of PHMMs.
19. A computer program product for of error correction in ancestry classification, the computer program product being embodied in a tangible computer readable storage medium and comprising computer instructions for:
- obtaining, from a classifier, initial ancestry classifications associated with portions of two phased haplotypes of a chromosome pair of an individual;
- performing error correction on an initial ancestry classification, including detecting a phasing error in the initial ancestry classifications; and
- outputting a corrected ancestry classification in which the phasing error is corrected.
Type: Application
Filed: Apr 10, 2018
Publication Date: Apr 18, 2019
Inventors: Chuong Do (Mountain View, CA), Eric Durand (Sceaux), John Michael Macpherson (Mountain View, CA)
Application Number: 15/950,023