Abstract
Background
DNA pooling is a technique to reduce genotyping effort while incurring only minor losses in accuracy of allele frequency estimates for single nucleotide polymorphism (SNP) markers.
Results
We present an algorithm for reconstructing haplotypes (alleles for multiple SNPs on same chromosome) from pools of two individual DNAs, in which HardyWeinberg equilibrium conditions or other assumptions are not required. The program outputs, in addition to inferred haplotypes, a minimal number of haplotypetagging SNPs that are identified after an exhaustive search procedure.
Conclusion
Our method and algorithms lead to a significant reduction in genotyping effort, for example, in casecontrol disease association studies while maintaining the possibility of reconstructing haplotypes under very general conditions.
Background
While SNPs in many ways are highly useful genetic markers, it may only be the joint effect of multiple SNPs in a gene that provides much information about association to disease because, taken together, the SNPs represent a much more polymorphic system than each SNP by itself. The effects of these SNPs are perhaps best represented by their haplotypes. Ideally, a researcher would want to obtain genotypes on all SNPs in a gene but this effort tends to be expensive. Thus, several authors have proposed that pools of DNA from n individuals be genotyped, which reduces genotyping costs by a factor of n [1]. In casecontrol association studies, an extreme approach is to form one pool for case and one pool for control individuals [2,3]. Originally, pooling was introduced as a method for efficiently screening for a rare disease [4] and thus for identifying individuals. On the other hand, grouping was later used to protect confidentiality [5].
Several methodological investigations of pooling efficiency have focused on marker allele frequency estimation [6,7] and on testing for association between markers and disease based on pooled data [8]. Extensions to two loci involve estimation of allele frequencies and the disequilibrium parameter based on pooled data [9].
Clearly, DNA pools of large numbers of individuals will only allow investigating SNP alleles and will be unable to look at haplotypes. A middle ground will be to form pools of small numbers of individuals. Such pools, for example, of n = 2 or n = 3 individuals are expected to lead to savings in genotyping costs while still allowing to recognize haplotypes. Here we present algorithms for inferring haplotypes from pools of n = 2 individuals. Our methods are applicable to any (reasonable) number of SNP markers (not limited to pairs of SNPs) and do not rely on the presence of HardyWeinberg equilibrium (HWE) as many other haplotyping approaches do, that is, they are largely parameterfree and datadriven.
Algorithm
Pool phenotypes
Consider a gene with s SNPs and a number m of pools (of n = 2 individuals). Denote the two alleles at a SNP by X and Y. Then five different phenotypes may be distinguished as shown in Table 1 (in addition to these five phenotypes there may be additional, more ambiguous situations, for example, that both X and Y alleles are known to be present but in unknown numbers).
Table 1. Pool phenotypes. "Assign?" indicates whether genotypes can be assigned to individuals in the pool, and "Homogeneous?" indicates whether the two individuals have the same genotype.
Of the five phenotypes, all but one (XXYY) allow an unambiguous assignment of SNP genotypes to the two individuals in a pool although, of course, the order of individuals is unknown. For example, phenotype XXXY clearly identifies one individual as XX and the other as XY. Among the four phenotypes with unambiguous assignment, one can distinguish "homogeneous" and "heterogeneous" pool phenotypes, where the former exhibit two identical and the latter two different SNP genotypes in the two individuals. For those SNPs that allow an unambiguous assignment of genotypes, genotypephase will be known (i.e., which individual has which genotypes at all SNPs) if at most one SNP is heterogeneous (the term genotypephase as we introduce it here is different from gametic phase, which refers to haplotypes).
Resolving ambiguities and missing information
The phenotype XXYY is ambiguous in the sense that the two individual phenotypes could be XX and YY, or XY each. Also, some phenotypes may be missing altogether. Before we can proceed with reconstructing haplotypes we need to "fill in such holes". We do this by imputing (partially) missing information and assume that the missing piece of information is missing at random, that is, it may be imputed based on the distribution of this type of information in the remainder of the data.
Missing phenotype
For a given SNP, a pool may be missing a phenotype completely. In this case we proceed in one of two ways. If among the remaining pools one phenotype is at least twice as frequent as the next frequent phenotype, then we impute the most frequent phenotype. Otherwise we impute one of the phenotypes at random.
Unknown number of X and Y alleles
A pool may be known to contain X and Y alleles for a given SNP but their numbers may be unknown. In other words, the phenotypes may be XYYY, XXYY, or XXXY. Here we count the respective numbers n_{X }and n_{Y }of X and Y alleles in the other pools and impute the genotype XYYY if n_{X }<n_{Y}, XXYY if n_{X }= n_{Y}, and XXXY if n_{X }>n_{Y}. Scoring of alleles in pooled DNA, for example, by mass spectrometry (fluorescence intensity) has a limited resolution so that occasionally the specific number of X or Y alleles is unclear.
Resolving genotypephase
As mentioned above, XXYY is the only phenotype not revealing genotypephase, that is, the individual genotypes cannot be recognized. For such a pool, we look at those other pools that do reveal genotypephase and count the number of individual genotypes as n_{XX}, n_{XY}, and n_{YY}. We assume genotypes (XX, YY) if n_{XX }+ n_{XY }>n_{YY }and assume genotypes (XY, XY) otherwise. If n_{XX }+ n_{XY }= n_{YY }then we make a random assignment based on equal probability of the two situations.
Because of the random element involved in some of these imputations they may to lead to different results at different times. Thus, it may be optimal to repeat these procedures a number of times and suitably combine results. Another option would be to allow for the missing observations by maximum likelihood (in the EM algorithm outlined below) but this would generally require consideration of an unsuitably large number of genotype vectors.
Estimating genotype vector frequencies
At this point, each pool unequivocally exhibits the SNP genotypes of each individual in the pool. We see the following six genotype pairs: XX/XX, YY/YY, XX/XY, XY/YY, XX/YY, and XY/XY. Among these genotype pairs, we distinguish two types: The two individuals in a pool have the same or different genotypes. We call the former homogeneous pools and the latter heterogeneous pools. The individual genotypes at different SNPs form a genotype vector that is analogous to a haplotype, which is a set of alleles at different loci. We recognize individual genotype vectors if among all SNPs at most one shows a heterogeneous phenotype.
Most of the time there will be multiple SNPs exhibiting heterogeneous phenotypes. Thus, we cannot generally count genotype vectors but must estimate them. We do this by a maximum likelihood method. By estimating frequencies of genotype vectors we do not need to make any assumptions on HardyWeinberg equilibrium. For s SNPs, there are a total of 3^{s }different genotype vectors. For example, s = 15 leads to 14,348,907 possible genotype vectors. In order to avoid having to enumerate all these, we focus on the subset compatible with the observed pools. For each pool, we generate a list of all genotype vectors that are compatible with the given pool. Initially, each of the different genotype vectors is assigned the same frequency, which is updated iteratively by an EMtype algorithm. At each iteration, we compute the sum of squared differences between current and previous genotype vector frequencies. Iteration stops when this sum of squares falls below 10^{11}.
Haplotype reconstruction
A number m of pools represent 2m genotype vectors (two for each individual in a pool). So, with given estimated frequencies, g_{i}, for the ith genotype vector, we prepare the following list: [2mg_{i}] is the estimated number of individuals with the ith genotype vector, where [r] represents the number r rounded to its nearest integer. For each genotype vector, its estimated frequency is converted into an estimated number of individuals with this genotype vector. From this list of individuals and their SNP genotypes, we now reconstruct individual haplotypes by Clark's method [10], which, again, does not assume HardyWeinberg equilibrium. At this point, of course, haplotypes may be reconstructed and their frequencies estimated with alternative approaches [1114] but most of these methods assume HardyWeinberg equilibrium. Although Clark's method does not necessarily lead to unique solutions, our main aim here was to continue analysis with a method not requiring the assumption of HardyWeinberg equilibrium. To evaluate the variability of the resulting haplotype frequency estimates it is recommended to carry out the whole procedure multiple times.
Implementation
As described above, we developed our approach in the following three steps: (1) For those pools not unequivocally exhibiting individual genotypes for a given SNP, we impute these genotypes based on known genotypes in the data. After this step we know the two genotypes at each SNP but not which person has which genotype. (2) For a given set of SNPs in a gene, the SNP genotypes of an individual form a genotype vector. We set up an EM algorithm to estimate the frequencies of all genotype vectors compatible with SNP genotypes and in this manner find the generally small number of genotype vectors with nonzero frequencies. (3) Based on the resulting SNP genotypes, Clark's algorithm [10] is applied to infer individual haplotypes and an exhaustive search of subsets of SNPs [15] is conducted to find the smallest number of such marker loci that uniquely identify ("tag") haplotypes (see add2 sample data, below), that is, that represent haplotypetagging SNPs (htSNPs). These algorithms have been implemented in a computer program package pools2 that is freely available to researchers. It contains the programs snp (for steps 1 and 2 above) and htSNP.py, and also the raw data for the add2 gene mentioned below. The package may be downloaded from the web site, http://linkage.rockefeller.edu/register webcite.
Results and Discussion
Efficiency of the snp program depends on the informativeness of the pool phenotypes. For example, a total of 25 SNPs and indels had originally been genotyped in the add2 gene (data file contained in the pools2 package mentioned above). But to evaluate all 25 or even only the first 20 of them resulted in over 100,000 genotype vectors compatible with pool phenotypes, which presumably would have required enormous computing times. Thus, analysis of the first 15 SNPs by the snp program is shown in the next paragraph. It required 31 iterations and an execution time of 13.7 minutes on a Dell PC (1.4 GHz clock speed, 786 MB of RAM, Windows 2000).
Analysis of add2 gene
For the first 15 SNPs in the add2 gene genotyped on 16 pools (two individuals each) of DNA we applied the procedure outlined above. Of the 6,252 genotype vectors compatible with the data (pool phenotypes), 12 have nonzero estimated frequencies. Haplotype reconstruction by the Clark method leads to 11 different haplotypes in the 2 × 16 = 32 individuals (Table 2).
Table 2. 15 SNPs and 11 haplotypes in the add2 gene. The two alleles at each SNP are labeled 0 and 1.
Some of the 15 SNPs are identical: s12 = s13 = s14 and s3 = s4. Thus, there are a total of 12 unique SNPs. Each haplotype forms a pattern of 0s and 1s at these SNPs. If for a subset of SNPs this pattern is different for each haplotype then that set of SNPs collectively tag the haplotypes and represent a set of htSNPs. To find the minimum number of haplotype tagging SNPs we proceed as follows [15]. Let h denote the given number of haplotypes and k the number of SNPs that tag these haplotypes. For example, two SNPs can generate at most four different allele patterns. Therefore, two SNPs cannot possibly tag more than four haplotypes. Generally, we must have 2^{k }≥ h or, equivalently,
k ≥ log(h)/log(2). (1)
To be certain that we find the smallest set of htSNPs we carry out an exhaustive search of all subsets of SNPs as previously proposed [15]. We begin with subsets of k SNPs, where k is the smallest integer satisfying equation (1). For each subset, the number t of different allele patterns is evaluated. If t = n we have found a set of htSNPs. We wrote a computer program, htSNP.py, to carry out this exhaustive search. For the 15 SNPs in Table 2, the program finds k = 6 as the smallest set of htSNPs and identifies 10 such sets. When we collapse s13 and s14 into s12, and s4 into s3, the following three sets of htSNPs remain: (s1, s3, s6, s7, s8, s15), (s1, s3, s6, s8, s9, s12), and (s1, s3, s6, s8, s9, s15). Thus, to identify the haplotypes in the add2 gene, genotyping effort can be reduced from genotyping 15 SNPs to genotyping only 6 SNPs. Which of these sets of htSNPs is optimal depends on the population haplotype diversity [15]; in the given sample, each set of htSNPs identifies all haplotypes perfectly but this is not necessarily the case in a new sample.
If we disregard the six haplotypes occurring only once we are left with five haplotypes comprising a total frequency of 90.6% (Table 2). For these common haplotypes, the only polymorphic and unique SNPs are s1, s3, and s11. They must be haplotype tagging because we need at least 3 SNPs to tag 5 haplotypes. Thus, to tag common haplotypes in the add2 gene it is sufficient to genotype those three SNPs.
Searching for htSNPs among all SNPs in a gene is somewhat related to looking for a cladistic structure among the haplotypes [16]. The most polymorphic SNPs are likely to be the oldest and tend to be included in the set of htSNPs. Thus, htSNPs tend to comprise those SNPs closest to an original disease mutation. A note of caution – while a possibly small number of htSNPs can identify all haplotypes, using only the htSNPs rather than all original SNPs in a gene does lead to reduced power in casecontrol disease association studies [17].
Uncertainties
Our procedures allow for some types of errors, for example, uncertain numbers of X and Y alleles in a pool, while assuming absence of other errors such as sample swaps [18]. The described method for imputing missing data may be suboptimal as it works with only one SNP at a time. We are currently working on improved approaches that take information from neighboring SNPs into account because SNPs in a gene tend to be highly correlated.
It should be pointed out that heterozygosity at multiple SNPs potentially introduces errors into the estimation of genotype vectors. Thus, some genotype vectors will be known with certainty and others only probabilistically. Also, the same pools causing genotypic ambiguity may also lead to haplotype ambiguity even given known genotypes. At this stage, we don't have any good answers to these observations except that increased sample size is expected to reduce such uncertainties.
Authors' contributions
JO and JH developed and implemented the algorithm and drafted the manuscript, MGL suggested the problem, FM and XP contributed single nucleotide polymorphism data, and DM carried out analyses.
Acknowledgements
References

Sham P, Bader JS, Craig I, O'Donovan M, Owen M: DNA pooling: A tool for largescale association studies.
Nature Reviews Genetics 2002, 3:862871. PubMed Abstract  Publisher Full Text

Arnheim N, Strange C, Erlich HH: Use of pooled DNA samples to detect linkage disequilibrium of polymorphic restriction fragments and human disease: Studies of HLA class II loci.
Proc Natl Acad Sci 1985, 82:69706974. PubMed Abstract

Plomin R, Hill L, Craig IW, McGuffin P, Purcell S, Sham P, Lubinski D, Thompson LA, Fisher PJ, Turic D, Owen MJ: A genomewide scan of 1842 DNA markers for allelic associations with general cognitive ability: a fivestage design using DNA pooling and extreme selected groups.
Behav Genet 2001, 31:497509. PubMed Abstract  Publisher Full Text

Dorfman R: The detection of defective members of large populations.

Gastwirth JL, Hammick PA: Estimation of the prevalence of a rare disease, preserving the anonymity of the subjects by group testing: application to estimating the prevalence of AIDS antibodies in blood.
J Stat Planning Inference 1989, 22:1527. Publisher Full Text

Jawaid A, Bader JS, Purcell S, Cherny SS, Sham P: Optimal selection strategies for QTL mapping using pooled DNA samples.
Eur J Hum Genet 2002, 10:125132. PubMed Abstract  Publisher Full Text

Shaw SH, Carrasquillo MH, Kashuk C, Puffenberger EG, Chakravarti A: Allele frequency distributions in pooled DNA samples: Applications to mapping complex disease genes.
Genome Res 1998, 8:111123. PubMed Abstract  Publisher Full Text

Risch N, Teng J: The relative power of familybased and casecontrol designs for linkage disequilibrium studies of complex human diseases. I. DNA pooling.
Genome Res 1998, 8:12731288. PubMed Abstract  Publisher Full Text

Pfeiffer RM, Rutter JL, Gail MH, Struewing J, Gastwirth JL: Efficiency of DNA pooling to estimate joint allele frequencies and measure linkage disequilibrium.
Genet Epidemiol 2002, 22:94102. PubMed Abstract  Publisher Full Text

Clark AG: Inference of haplotypes from PCRamplified samples of diploid populations.
Mol Biol Evol 1990, 7:111122. PubMed Abstract  Publisher Full Text

Chiano MN, Clayton DG: Fine genetic mapping using haplotype analysis and the missing data problem.
Ann Hum Genet 1998, 62:5560. PubMed Abstract  Publisher Full Text

Zhao JH, Curtis D, Sham PC: Modelfree analysis and permutation tests for allelic association.
Hum Hered 2000, 50:133139. PubMed Abstract  Publisher Full Text

Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data.
Am J Hum Genet 2001, 68:978989. PubMed Abstract  Publisher Full Text

Lin S, Cutler DJ, Zwick ME, Chakravarti A: Haplotype inference in random population samples.
Am J Hum Genet 2002, 71:11291137. PubMed Abstract  Publisher Full Text

Johnson GC, Esposito L, Barratt BJ, Smith AN, Heward J, Di Genova G, Ueda H, Cordell HJ, Eaves IA, Dudbridge F, et al.: Haplotype tagging for the identification of common disease genes.
Nat Genet 2001, 29:233237. PubMed Abstract  Publisher Full Text

Templeton AR, Weiss KM, Nickerson DA, Boerwinkle E, Sing CF: Cladistic structure within the human Lipoprotein lipase gene and its implications for phenotypic association studies.
Genetics 2000, 156:12591275. PubMed Abstract  Publisher Full Text

Zhang K, Calabrese P, Nordborg M, Sun F: Haplotype block structure and its applications to association studies: Power and study designs.

Kirk KM, Cardon LR: The impact of genotyping error on haplotype reconstruction and frequency estimation.
Eur J Hum Genet 2002, 10:616622. PubMed Abstract  Publisher Full Text