Nonparametric Disequilibrium Mapping of Functional Sites Using Haplotypes of Multiple Tightly Linked Single-Nucleotide Polymorphism Markers

Corresponding author: Ming D. Li, Mail Code 7792, 7703 Floyd Curl Dr., San Antonio, TX 78229., lim2{at}uthscsa.edu (E-mail)

Communicating editor: Z-B. ZENG

	ABSTRACT

TOP ABSTRACT METHODS RESULTS DISCUSSION LITERATURE CITED

As the speed and efficiency of genotyping single-nucleotide polymorphisms (SNPs) increase, using the SNP map, it becomes possible to evaluate the extent to which a common haplotype contributes to the risk of disease. In this study we propose a new procedure for mapping functional sites or regions of a candidate gene of interest using multiple linked SNPs. Based on a case-parent trio family design, we use expectation-maximization (EM) algorithm-derived haplotype frequency estimates of multiple tightly linked SNPs from both unambiguous and ambiguous families to construct a contingency statistic S for linkage disequilibrium (LD) analysis. In the procedure, a moving-window scan for functional SNP sites or regions can cover an unlimited number of loci except for the limitation of computer storage. Within a window, all possible widths of haplotypes are utilized to find the maximum statistic S* for each site (or locus). Furthermore, this method can be applied to regional or genome-wide scanning for determining linkage disequilibrium using SNPs. The sensitivity of the proposed procedure was examined on the simulated data set from the Genetic Analysis Workshop (GAW) 12. Compared with the conventional and generalized TDT methods, our procedure is more flexible and powerful.

MOST human disorders of interest likely result from the cumulative effect of alleles at multiple susceptibility loci, none of which on its own is either necessary or sufficient to cause the disease. Because of this, the classical strategies of analyzing monogenic disorders have been unsuccessful. Alternative approaches such as genome-wide linkage and association analysis have been proposed and utilized in many studies (e.g., see RISCH 2000 ). Large-scale association studies would employ a dense single-nucleotide polymorphism (SNP) map to detect association between markers and disease. In the causal hypothesis, most genetic determinants of a disease are SNPs that are likely to be selected as markers. In the proximity hypothesis, most disease determinants will not be included among markers but may be detected through linkage disequilibrium with other SNPs. Simulation studies based on monotonic population expansion suggest that useful association does not usually extend beyond 3 kb along the genome (COLLINS et al. 1999 ; KRUGLYAK 1999 ). Recent advances in the identification of SNP haplotype blocks in the human genome (PATIL et al. 2001 ; GABRIEL et al. 2002 ; ZHANG et al. 2002A ) showed that using selected representative SNPs by haplotype block partition algorithms can dramatically reduce the time and effort for genotyping without losing much haplotype information. Therefore, association studies using SNP markers represent an important tool in identifying susceptibility loci for human disorders because of their ability to identify unique chromosomal segments likely to harbor disease-predisposing genes.

Many haplotype analysis methods in the literature require phase information inferred from genotype data. However, as the number of loci increases, the information loss due to haplotype ambiguity would increase rapidly (HODGE et al. 1999 ). Several strategies that make use of the expectation-maximization (EM) algorithm (DEMPSTER et al. 1977 ) have been proposed to overcome the problem of missing phase information for estimating haplotype frequencies from unphased diploid genotype data (EXCOFFIER and SLATKIN 1995 , EXCOFFIER and SLATKIN 1998 ; HAWLEY and KIDD 1995 ; LONG et al. 1995 ; SLATKIN and EXCOFFIER 1996 ; CHIANO and CLAYTON 1998 ; CLAYTON 1999 ). Recently, several research groups have discussed and compared the accuracy of haplotype frequency estimation using the EM algorithm (FALLIN and SCHORK 2000 ; TISHKOFF et al. 2000 ; STEPHENS et al. 2001 ; NIU et al. 2002 ; QIN et al. 2002 ). It was suggested that even in the worst case, individual haplotype frequency estimates obtained by the EM algorithm from a sample size of 100 or more would not deviate >5% from their true values among the sampled individuals (FALLIN and SCHORK 2000 ). Furthermore, FALLIN and SCHORK 2000 suggested that one should be more concerned about the quality of sampling than the possibility of estimation errors when assessing haplotypes among unphased individuals. Therefore, EM estimation of haplotype frequencies for multiple diallelic genotypes may represent a viable alternative to the recruitment of additional family members or intensive laboratory haplotyping for haplotype-based genetic studies.

Haplotype frequency estimates from tightly linked multilocus genotyping data have been used for linkage-disequilibrium (LD) analysis (CLAYTON 1999 ; TOIVONEN et al. 2000 ; ZHAO et al. 2000 ; FALLIN et al. 2001 ; SCHAID et al. 2002 ). For example, ZHAO et al. 2000 proposed a new statistical method for the transmission/disequilibrium test (TDT) using EM-derived haplotype frequency estimates of multiple tightly linked markers from both unambiguous and ambiguous families. Their results suggest that the power of the method was higher than that of the conventional TDT method (SPIELMAN et al. 1993 ) or other TDT methods, such as using estimated haplotype frequencies of multiple markers from only unambiguous families (WILSON 1997 ; CLAYTON and JONES 1999 ). CLAYTON 1999 also proposed a new TDT statistic for the situation in which some haplotype phases are unknown or some parental genotypes are missing. DUDBRIDGE et al. 2000 showed that, in some families for which haplotypes are known, a potentially serious bias is introduced into the TDT if the loss of information from families with ambiguous haplotypes is not taken into consideration. The methods proposed by CLAYTON 1999 and ZHAO et al. 2000 can avoid loss of information from families with ambiguous haplotypes, but both require one to preassign a window width prior to analyzing multisite parental transmission data under this given fixed width.

SACHIDANANDAM et al. 2001 pointed out that in practice, when a gene has been implicated in causing the disease of interest as a result of linkage analysis, known biological function(s), or an expression pattern, it is desirable to survey allelic variation within it exhaustively for any potential association to the disease. Using SNP markers, it becomes possible to evaluate the extent to which common haplotypes contribute to disease risk. In this study, we apply EM-derived haplotype frequency estimates of multiple tightly linked SNP markers from both unambiguous and ambiguous families to construct a contingency table statistic, namely S(a, b), for LD analysis. Our procedure assumes no or very rare recombination among multiple tightly linked markers. In our procedure, a moving-window scan for functional sites or regions can cover any number of loci without limitation except for that imposed by computer storage. Within a given maximum window width, all possible widths of haplotypes are utilized to find the maximum statistic S* for each site (or locus). This represents a significant difference between the procedure reported here and the published methods in the literature (CLAYTON 1999 ; ZHAO et al. 2000 ). We apply the proposed method to scan the simulated sequence data from the Genetic Analysis Workshop (GAW) 12 (ALMASY et al. 2001 ) to identify the functional sites of candidate genes with tightly linked SNP markers. Comparing the results from our procedure with those from the conventional and generalized TDT methods, it appears that our method is more flexible and powerful.

	METHODS

TOP ABSTRACT METHODS RESULTS DISCUSSION LITERATURE CITED

We describe a nonparametric method of detecting functional SNP site(s) of haplotypes associated with the disease of interest in a LD study. The method naturally incorporates tightly linked multilocus genotype data, in the sense that multiallelic loci are utilized and the existing disequilibrium among markers (in both haplotypes that are transmitted and those that are not transmitted to affected offspring) is built into the test. The approach searches for evidence of ancestral haplotypes that are shared more often among chromosomes transmitted than among chromosomes not transmitted to affected offspring, followed by comparison of the observed data to the distribution of data expected under the null hypothesis.

EM algorithm estimation:
When a haplotype segment with certain marker(s) occurs in a chromosome transmitted to affected offspring with a frequency higher than that in the chromosome not transmitted, there exists an association of the haplotype with the disease. This finding would be considered striking and consistent with the hypothesis that the marker positions are near a disease gene or located at a disease locus (e.g., a functional site lies within the sequence of a disease gene), explaining the haplotype as being identical by descent (IBD) with that from a common founding ancestor. We consider N randomly sampled trios, each containing an affected offspring along with two parents, with each member genotyped for multiple tightly linked markers. In some trios, it may be impossible to unambiguously reconstruct haplotypes from the given genotypes. In such a design with data from an affected singleton offspring and both parents, if the haplotypes can be reconstructed from the known genotypes unambiguously, we let H be the total number of possible haplotypes and n_j and n^'_j be numbers of a particular haplotype h_j (j = 1, 2, ... , H) transmitted and not transmitted to the affected offspring, respectively (see Table 1). If the haplotypes cannot be reconstructed unambiguously, then, to fully use the information of the multiple tightly linked markers, the EM algorithm is employed to estimate the transmitted and not transmitted haplotype frequencies f_j and f^'_j (j = 1, 2, ... , H). In Table 1, we use 2N x f_j and 2N x f^'_j to represent directly counted numbers n_j and n^'_j, respectively, even when some of the actual numbers are ambiguous.

View this table: In this window In a new window   Table 1. H x 2 table for calculating the contingency table statistic S (S*(x) = S(a, b))

Suppose G_i (= g_i, g_i,m, g_i,f) is the trio for the ith family, where g_i, g_i,m, and g_i,f represent the genotypes of the child, the mother, and the father, respectively. Let

• h_u be paternal haplotype transmitted to the affected offspring

• h_v be paternal haplotype not transmitted to the affected offspring

• h_l be maternal haplotype transmitted to the affected offspring

• h_m be maternal haplotype not transmitted to the affected offspring.

Conditional on G_i, the probability or weight to a transmitted haplotype h_j from the ith family is

where (h_u, h_v, h_l, h_m) G_i denotes the haplotype group (h_u, h_v, h_l, h_m) that is compatible with the genotype group of the ith family and the factor c^j_ul(i) (= 0, 1, or 2) depends on the counts of haplotype frequency f_j occurring in the pair of haplotype frequencies (f_u, f_l). The summation in the denominator is over all haplotype groups that are compatible with G_i. The w⁽ⁱ⁾_j are defined as the estimated weight of the ith family when the parent(s) has haplotype (h_u=j, h_v) and/or (h_l=j, h_m), and h_j is transmitted. Similarly,

is the weight for a not-transmitted haplotype h_j.

Given the genotypes, the likelihood function for the data is

where

So the log-likelihood is

Next we consider estimation of the haplotype frequencies f_j and f^'_j (j = 1, 2, ... , H). The EM algorithm (see, e.g., EXCOFFIER and SLATKIN 1995 , EXCOFFIER and SLATKIN 1998 ; HAWLEY and KIDD 1995 ; LONG et al. 1995 ; SLATKIN and EXCOFFIER 1996 ; CHIANO and CLAYTON 1998 ; CLAYTON 1999 ) is a way of attempting to find the f_j and f^'_j (j = 1, 2, ... , H) that maximize the log-likelihood logL(data|G₁, G₂, ... , G_N). Similar to the study reported by ZHAO et al. 2000 , we assume that all parents represent a random sample of unrelated individuals from a population in Hardy-Weinberg equilibrium. The entire set of weights, w⁽ⁱ⁾_j and w^'(i)_j (i = 1, 2, ... , N; j = 1, 2, ... , H) can thus be computed in a single pass through the H possible haplotypes as indicated below:

In principle, the maximum-likelihood (ML) estimates of the haplotype frequencies could be found analytically by solving a set of equations with the Lagrange multipliers ₁ and ₂:

and the additional partial derivatives with respect to ₁ and ₂. We obtain the ML estimators of f_j and f^'_j (j = 1, 2, ... , H) in the (p + 1)st EM iteration,

Initially, all haplotypes are set equally frequent, so that all possible complementary haplotype pairs are equally likely. In the (p + 1)st expectation step (E-step), the weights of the transmitted and not-transmitted haplotypes, w^(i)(p+1)_j and w^'(i)(p+1)_j, can be obtained from the given genotype information and current estimates of the haplotype frequencies, ^(p)_j and ^'(p)_j. The (p + 1)st maximization step (M-step) gives the maximum-likelihood estimates of ^(p+1)_j and ^'(p+1)_j. When the difference between the values of the haplotype estimates in the previous M-step and current one becomes less than a predetermined quantity (e.g., 10^-7), the iteration is stopped and the final estimates are obtained. The estimated numbers _j (=2N x _j) and ^'_j (=2N x ^'_j) of transmitted and not-transmitted haplotypes are then used for further computing the statistic S and estimating the empirical 100(1 - ) percentiles and P value by simulation permutation (CHURCHILL and DOERGE 1994 ).

Statistic S:
The key feature of this method is encompassed by defining a statistic S(a, b) (GAO and WRIGHT 1999 ) that is computed on the set of haplotypes beginning at marker position a and ending at position b. The association evidence at position x is

The statistic S is chosen to reflect a striking association of all haplotypes with the disease. In other words, for an arbitrary statistic S and at position x, we search over all haplotype widths (containing x) to find the haplotype in most striking association with the disease. Under the null hypothesis, this haplotype has no association with the disease. Under the assumption that haplotypes are sampled independently, we compute a P value for the ² contingency table test of haplotype vs. disease status and refer to the statistic S(a, b) = -log₁₀ (P value) as S. The P value corresponds to the k x 2 standard ² contingency table test of k (1 < k H) unique haplotypes (with k - 1 d.f.) that begin at position a and end at b vs. disease status, under the assumption of independence. The number k of unique haplotypes will depend on the choice of a and b.

Permutation for empirical P value:
To evaluate type I error accurately, an empirical P value that appropriately corrects for the testing of multiple marker locations is obtained from the EM estimates of haplotype numbers and frequencies (_j, ^'_j, _j, and ^'_j; j = 1, 2, ... , H), by permuting the transmitted status of individual haplotypes, and computing S*(x) for each permutation (ZHAO et al. 2000 ). Briefly, we first generate haplotypes conditional on each family's genotypes. If multiple possible compatible haplotype groups are conditional on a given family's genotypes, first a compatible haplotype group is randomly selected. For each family, haplotypes are then randomly permuted within each transmitted/nontransmitted haplotype pair for each parent to form a new set of family genotypes. The transmitted and nontransmitted haplotype frequencies are then reestimated using the EM algorithm for each random selection and permutation of the data. By using this procedure, we obtain the null distribution of the maximum of S* from 10,000 randomized permutations. Then it is easy to obtain the 100(1 - ) percentile and appropriate P value for the maximum of S* over the region with multiple tightly linked markers. A computer program for the proposed method is available upon request.

	RESULTS

TOP ABSTRACT METHODS RESULTS DISCUSSION LITERATURE CITED

We applied the proposed procedure to the simulated sequence data set of GAW 12. This data set contains computer-simulated sequence data with multiple SNP markers for seven candidate genes in 23 extended pedigrees with a total of 1497 individuals for two populations: a general and an isolated population. There were 50 replicates for each population model. We randomly sampled 10 replicates from the 50 replicates of the general population and separated them into two groups each with five replicates. Group 1 contained replicates 1, 10, 18, 42, and 48, while group 2 contained replicates 21, 23, 33, 34, and 38. For each extended pedigree in a replicate, we randomly sampled a trio with one affected offspring and two parents (affected or unaffected), for whom genotyping information was available. In total, 95 and 97 trios were in groups 1 and 2, respectively. Last, we pooled the two samples together to form sample 3, with a total of 192 trios. All the SNP variants present in a sample were counted, and only those with a frequency >1% were chosen as possible SNP sites (BRODER and VENTER 2000 ).

Theoretically, the EM algorithm can be applied to an unlimited number of loci with any number of alleles. However, in practice, implementation of this algorithm is limited by the need to store the estimated haplotype frequencies for every possible haplotype contained in the sample. These storage requirements increase exponentially with the number of loci under investigation. For example, the numbers of SNPs for candidate gene 1 are 155, 158, and 152 in samples 1, 2, and 3, respectively. If any individual is heterozygous at 150 loci, then the number of possible haplotypes in that sample is 2¹⁵⁰. Also, as the number of markers increases, there will be an increased variance in the estimates. In this study, we first selected SNP markers at 1000-bp intervals and set the window width at 5 (i.e., m = 5; five tightly linked markers) to scan every selected SNP marker for all candidate gene sequences. The procedure of window scanning implemented in the program can be described as follows: a pattern is set up where the first window consists of markers 1–5 (i.e., 1 a x b 5) and the second window of markers 2–6 (i.e., 2 a x b 6), and the shifting of the window continues until all of the SNPs of a gene or genomic fragment have been scanned. The chromosome locations, lengths of sequence data, and number of SNPs for these candidate genes are shown in Table 2.

Linkage disequilibrium analysis was performed for all selected SNPs of these seven candidate genes. The distributions of the statistic S* for these candidate genes are shown in Fig 1. The highest peak locations of S* and their corresponding P values for all candidate genes are summarized in Table 3. As described in METHODS, for an arbitrary statistic S and at the position x, we searched over all haplotype widths. If we set m = 5 and x = 5, the varied pairs of (a, b) containing the x will be (1, 5), (2, 6), (3, 7), (4, 8), (5, 9), (2, 5), (3, 6), (4, 7), (5, 8), (3, 5), (4, 6), (5, 7), (4, 5), (5, 6), and (5, 5), and the statistic S calculated from the last pair of (5, 5) is equivalent to a conventional single-marker test. Some of the peak locations of S* in Table 3 contained two linked markers with an equal maximum S*, while others might contain three, four, or five linked markers. One functional SNP was identified for candidate genes 1 and 6. No functional SNP site was detected at the 0.01 significance level for the other five candidate genes (Table 3). In fact, for candidate gene 2 multiple simulated functional sites directly affected a quantitative trait, but no single functional site directly affected disease status. For the other four candidate genes, no functional sites were simulated in the original model (ALMASY et al. 2001 ). For candidate gene 1, the functional site at the nucleotide 557 position was detected in samples 2 and 3, as expected from the underlying assumption of the simulated data set. For candidate gene 6, a region from nucleotides 5782 to 7332 (i.e., the maximum statistic S* obtained using the haplotype data from two tightly linked loci), which contains the simulated functional site at nucleotide 5782, was detected in both samples 1 and 3. For candidate gene 1 in sample 1 and candidate gene 6 in sample 2, the statistic S* peaks were detected at the same locations as in the other samples, except that the corresponding P values did not reach the 0.01 level. It should be noted that, after pooling samples 1 and 2 to form sample 3, the statistic S* values at the peak locations increased as the sample size became larger for both candidate genes 1 and 6, but did not increase, in fact even decreased, for the other five candidate genes at some peak locations (Fig 1).

View larger version (24K): In this window In a new window Download PPT slide   Figure 1. The chromosomal distribution of statistic S* for seven candidate genes. The lines with open circles, triangles, and squares represent the results from samples 1, 2, and 3, respectively. Average SNP interval for all candidate genes is 1000 bp.

Next, we compared our procedure to the generalized TDT (CLAYTON 1999 ) and the conventional TDT (SPIELMAN et al. 1993 ) methods for candidate genes 1 and 6 in samples 1–3. For the conventional TDT method, each SNP site was studied separately, leading to a P value for each SNP site. For the generalized TDT method implemented in the TRANSMIT program (CLAYTON 1999 ), we fixed the window widths at two, three, four, or five linked SNPs followed by choosing the smallest P value for each SNP site from the multiple-locus and single-locus analysis results under different windows. As shown in Fig 2, for candidate gene 1, the highest peaks of the statistic S* from our procedure are greater than (or at least equal to) those obtained from the generalized or conventional TDT methods for all samples. For candidate gene 6, the largest statistics obtained from our procedure are all greater than those from the TRANSMIT or TDT procedures. The highest peaks of statistic S* from the TRANSMIT program are greater than or equal to those from the conventional TDT method.

View larger version (34K): In this window In a new window Download PPT slide   Figure 2. Comparison of linkage disequilibrium analyses among our procedure (nonparametric disequilibrium, ND), the generalized TDT (using the program TRANSMIT; CLAYTON 1999 ), and the conventional TDT method, for candidate genes 1 and 6 on samples 1, 2, and 3. The lines with open circles, triangles, and squares indicate the results from ND, TRANSMIT, and TDT, respectively. The average SNP interval for all candidate genes is 1000 bp.

Fig 3 compares the results obtained from our procedure with the generalized and conventional TDT methods for candidate gene 1 in sample 3 (N = 192 trios) under different SNP densities. A total of 41, 23, 8, and 5 SNPs were selected from the original SNPs for candidate gene 1, which gave an average SNP interval of 500, 1000, 3000, and 5000 bp, respectively. The simulated functional site (i.e., 557 bp) was included for only the first two high densities. As shown in Fig 3, the highest peaks of statistic S* detected by these methods are all at or around the 557-bp functional site. Compared to what is found with high densities, the highest peaks of statistic S*'s detected at low densities are smaller. Similarly, we compared the results from these three methods for candidate gene 6 under three different SNP densities. As shown in Fig 4, the functional site is detected by all three methods for candidate gene 6 under the SNP density of 1000- and 2000-bp intervals, with the highest peak at functional site 5782 bp for 1000-bp intervals and at SNP site 7332 bp (which is 1550 bp away from the functional 5782-bp site) for 2000-bp intervals. Maximum statistics S* obtained from our procedure are always greater than those obtained from the generalized or conventional TDT (see Fig 4A and Fig B). When 5 SNPs at 4000-bp intervals were used, no linkage disequilibrium was detected by any method, but our procedure and the generalized TDT method still produced smaller P values than those from the conventional TDT (Fig 4C).

View larger version (16K): In this window In a new window Download PPT slide   Figure 3. Comparison of linkage disequilibrium analyses of our procedure with the generalized TDT (TRANSMIT; CLAYTON 1999 ) and the conventional TDT for candidate gene 1 under SNP densities of 500-, 1000-, 3000-, and 5000-bp intervals. The lines with open circles, triangles, and squares represent the results from the ND procedure, the generalized TDT, and the conventional TDT method, respectively.

View larger version (18K): In this window In a new window Download PPT slide   Figure 4. Comparison of linkage disequilibrium analyses of the ND procedure with the generalized TDT and conventional TDT methods for candidate gene 6 under SNP densities of (A) 1000-, (B) 2000-, and (C) 4000-bp intervals, respectively. The lines with open circles, triangles, and squares represent the results from our procedure, the generalized TDT, and the conventional TDT, respectively.

	DISCUSSION

TOP ABSTRACT METHODS RESULTS DISCUSSION LITERATURE CITED

To date, 2.1 million SNPs have been identified in the human genome (VENTER et al. 2001 ). This rich resource should allow the initiation of genome-wide linkage disequilibrium mapping of disease genes in the human population (LANDER et al. 2001 ). In the past, many LD analysis methods have been proposed (e.g., FALK and RUBINSTEIN 1987 ; OTT 1989 ; THOMSON et al. 1989 ; TERWILLIGER and OTT 1992 ; SPIELMAN et al. 1993 ; SHAM and CURTISS 1995; THOMSON 1995 ; SCHAID 1996 ; SPIELMAN and EWENS 1996 , SPIELMAN and EWENS 1998 ; MARTIN et al. 1997 , MARTIN et al. 2000 ; SHAM 1997 ; XIONG and GUO 1997 ; LAZZERONI and LANGE 1998 ; CORDELL and ELSTON 1999 ; SUN et al. 1999 ; TENG and RISCH 1999 ; ZHENG and ELSTON 1999 ; ZHU and ELSTON 2000 , ZHU and ELSTON 2001 ) and proved to be powerful approaches to identify susceptibility genes for the disease of interest. However, to apply most, if not all, of these methods, one needs to either discard families with ambiguous haplotypes or analyze the markers separately, resulting in potential loss of power.

In this study, we propose a new procedure for nonparametric disequilibrium mapping at a functional site (or region) within a candidate gene by using multiple tightly linked SNPs. For an arbitrary statistic S at position x, we search over all haplotype widths (containing x) and find that some of the peak locations of the maximum statistic S* contain two linked markers while others might contain three, four, or five linked markers (see Table 3). With our procedure, researchers do not need to know prior to analysis what window width of haplotype is more appropriate for detecting a functional site or region. In contrast, the methods proposed by CLAYTON 1999 and ZHAO et al. 2000 require assigning a window width first and then analyzing multisite parental transmission data under this fixed width. Therefore, our procedure is more flexible, and it can find the maximum statistic S* with an appropriate window width of haplotype transmission data within a preset range (m) spanning a single, two, or more multiple linked markers (m). Although it is generally thought that a multiple-locus approach may be more powerful than single-locus analysis, this does not mean that this is true in all cases. For example, for candidate gene 1 in the simulated GAW 12 data set, it seems that the single-marker approach is better than that of multiple markers. In contrast, for candidate gene 6, we found that using the haplotype transmission data of two linked markers appears to be better than using those of a single marker, or three, four, or five linked markers.

The simulation model implemented in the GAW 12 data set was rather complex (ALMASY et al. 2001 ). Seven major genes (MG) were assumed to influence one or more of five quantitative traits (Q1–Q5), disease liability, and age at onset. The disease liability was constructed as a function of Q1–Q5, MG6, and a mitochondrial component. MG6 provides a baseline liability value, which is adjusted for Q1–Q5. MG6 and the mitochondrial component have direct effects on liability and account for 22 and 10% of its variance, respectively, conditional on other covariates like gender, age, and environmental factors. MG1–MG5 affect liability only through the quantitative traits. MG1 accounts for 4.6% of the variance in liability, MG2 for 4.4%, MG3 for 8%, MG4 for 0.3%, and MG5 for 9%. Our procedure could detect the functional site for candidate gene 1 that directly affects liability and the functional site for candidate gene 6, where only one functional site directly contributes to Q1 and Q2 and indirectly affects the affection status. But we could not detect the multiple functional sites for candidate gene 2, a quantitative trait gene with multiple functional sites. Therefore, for quantitative traits controlled by multiple genes with relatively small effects or multiple functional sites each with a small genetic effect, such as candidate gene 2, other methods of using haplotype information should be considered. Several methods focusing on this type of data have been proposed (CLAYTON 1999 ; LI et al. 2001 ; SCHAID et al. 2002 ).

Using the procedure reported in this study, linkage disequilibrium due to the two simulated functional sites of candidate genes 1 and 6 (defined as MG6 and MG1, respectively, in the original model) were detected. Multiple alleles were in candidate gene 2; all changes in regulatory elements or in the first or second base pair of a codon leading to amino acid substitutions were functional. However, this was not the case for candidate gene 6. Furthermore, candidate gene 2 directly contributes to Q5 from multiple functional sites and then indirectly affects the affection status.

On the basis of the mapping results from our method and the generalized and conventional TDT methods on the GAW 12 simulated data set, it appears that our procedure is as sensitive and powerful as the TDT methods for candidate gene 1, but is more powerful than the TDT methods for candidate gene 6. Additionally, we investigated how SNP density may affect the performance of our procedure. A functional SNP site for candidate gene 1 was detected by all three methods, regardless of SNP density and whether the functional SNP site was included or not. In contrast, a functional SNP site for candidate gene 6 was detected only at high SNP densities, not at low density. To explain this difference, we examined selected SNPs used under each density for both genes and found that a SNP at the 189-bp position of gene 1, which is only 368 bp away from the functional SNP site, was used for all analyses. However, for candidate gene 6 at the 4000-bp intervals, the closest SNP sites used in the analysis were at the 4848- and 9952-bp positions, which are either loosely linked or too far away from the functional site (5782 bp). Regardless of which SNP density was used, overall, we found that our procedure and the generalized TDT method are less dependent on the SNP density than the conventional TDT method is. Moreover, our program is more user friendly and less time consuming than the TRANSMIT program with respect to file preparation before and after mapping analysis. For example, for candidate gene 1 under a density of 500-bp intervals, only 1 file was needed for our procedure, while for TRANSMIT we had to prepare 40, 39, 38, and 37 files for window widths of two, three, four, and five linked SNPs, respectively. Although these files can be prepared using a shell script for each selected window width, we still have to summarize a significant number of output files for comparing the results.

Theoretically, our search procedure can cover as many loci as one wishes without limitation, except for the computer storage requirement. As needed, our method can be used for regional or genome-wide scanning to study LD using SNP markers. This represents another advantage of our program compared to other currently available procedures. It should be pointed out that if genome-wide high-density SNP data were available, then the size of the data file would become significantly large. The newly proposed haplotype block partition algorithms (PATIL et al. 2001 ; ZHANG et al. 2002A ) can be employed to find appropriate blocks and select representative SNPs. To overcome the difficulty of extreme time cost and the limitation of computer storage in estimating large haplotypes of SNPs, NIU et al. 2002 and QIN et al. 2002 proposed a partition-ligation (PL) strategy together with Gibbs sampling and an EM algorithm. Their PL-EM algorithm (QIN et al. 2002 ) can be applied to haplotype estimation within a window for a large number of SNPs (e.g., >20). These multimarker approaches can effectively extract information from a case-parent trio as well as case-control data and increase the efficiency and detection power in regional or genome-wide LD analysis.

An alternative approach to EM in inferring haplotypes has been suggested by CLARK 1990 and used in several recent studies (REICH et al. 2001 ; STEPHENS et al. 2001 ). However, as discussed by CLARK 1990 , several problems may arise with his procedure, e.g., the possibility of never being able to start the iterative algorithm because of the absence of any unambiguous individuals in the sample. Instead of finding only a list of possible haplotypes by Clark's algorithm, the EM approach has the advantage of estimating haplotype frequencies, which in turn can be used directly for further ² testing in LD analysis. STEPHENS et al. 2001 proposed a new statistical method, applicable for genotype data at linked loci, for haplotype reconstruction from population data, based on an iterative sampling algorithm using a pseudo-Gibbs sampler. Their simulated results showed that error rates were reduced by >50% relative to its nearest competitor (EM-derived haplotypes). However, ZHANG et al. 2001 compared the two methods and pointed out that one limitation of Stephens et al.'s study is the fact that their simulations were based on coalescent models, which may not be good approximations of human population evolutionary histories. Furthermore, their method cannot be applied to estimate the transmission haplotype frequencies. The PL-Gibbs sampler strategy proposed by NIU et al. 2002 not only can infer the haplotype phases with a large number of SNPs, but is also supposed to run faster than the method proposed by STEPHENS et al. 2001 . It is possible to extend their method to estimate the affected (case) and unaffected (control) haplotype frequencies from population data and to estimate transmission haplotype frequencies for independent case-parent trio data as done in ZHAO et al. 2000 , together with the procedure reported herein. Indeed, this idea has been extended for case-control association analysis with multiple tightly linked SNP markers by using EM-derived haplotype frequencies (FALLIN et al. 2001 ). For example, if the haplotype frequencies or numbers of multiple SNPs are estimated by the PL-EM algorithm (QIN et al. 2002 ) or the algorithm proposed by NIU et al. 2002 , we can calculate the maximum statistic S* for each site with a wider moving window to find associated functional sites or regions, using simulation permutation to obtain P values as indicated here.

In the current version of our program, the missing genotype of a parent and/or the affected offspring at any SNP site is assigned a code that differs from the codes assigned to the two possible SNPs at the position. As expected, this approach will yield rare haplotypes, which may eventually affect the power of the test to detect functional sites.

The chi-square statistic is highly sensitive to small cell counts. Recently, SELTMAN et al. 2001 pointed out that although haplotypes are more informative than alleles of individual SNPs in the transmission disequilibrium-based approach, the larger number of haplotypes relative to alleles at individual loci might decrease the statistical power of haplotype-based TDT. The reason for this may be the larger number of degrees of freedom required. As indicated earlier, our program automatically clusters haplotypes into fewer groups, and eventually a maximum S(x) with the smallest P value is chosen as S*(x) at locus x. By applying the method to the simulated data sets, we found that the shorter the window width, the smaller the number of degrees of freedom. This suggests that the clustering of haplotypes into fewer groups may, in fact, improve the power of the test.

The number of possible haplotypes increases rapidly with the size of the searching window. To avoid very rare haplotypes, we compared sample sizes from 100 to 1000 and found that the maximum searching window width should be set in the range approximately four to eight tightly linked SNPs (equivalent to 16–256 possible haplotypes). By empirical study, NIU et al. 2002 and QIN et al. 2002 also suggested that five to eight SNPs appeared to be a good choice for the atomistic unit size before the ligation step. Specifically, the maximum width was set at 5 in all the LD analyses reported in this study. If the sample size is large enough, it is sufficient to use a single initial condition for the haplotype frequencies (e.g., all haplotypes are set equally frequent) to start the EM estimation (HAWLEY and KIDD 1995 ; SLATKIN and EXCOFFIER 1996 ). Given the sample size in the GAW 12 data set, we used this single initial condition. On the other hand, if the sample size is not large enough, then multiple initial conditions should be used for EM estimation to avoid going to a local optimum.

An example has been given for the first five selected SNP loci of candidate gene 1 to illustrate our procedure. Table 4 shows the estimated numbers of haplotypes found by using the EM algorithm for the first five selected SNPs of candidate gene 1 and P values for the maximum statistic S* obtained from 10,000 randomized permutations in all samples. The nonancestral and ancestral sequence variants of candidate gene 1 were coded as alleles 1 and 2, respectively. On the basis of the estimated haplotype numbers, it is easy to see evidence of ancestral haplotypes that are shared more often among the haplotypes transmitted than among those not transmitted to affected offspring. In this case, allele 2 at base pair 557 is associated with the disease.

Recently, ZHANG et al. 2002B applied the haplotype block partition and "tag SNPs," distinguishing algorithms to association studies. After haplotype block partition, tag SNPs were further selected within each haplotype block and used for association analysis. Their simulation results indicated that the genotyping efforts can be significantly reduced by using the tag SNPs without much loss of power, and haplotype-based LD analysis can be much more powerful than marker-by-marker analysis. It seems that a more effective and powerful searching strategy should combine the procedures of identifying haplotype block structure and the corresponding tag SNPs with the PL-EM (QIN et al. 2002 ) and/or PL-Gibbs sampling (NIU et al. 2002 ) algorithms as well as our moving-window statistic S(x) and permutation methods.

The procedure proposed in this study is nonparametric and robust. For a complex disease controlled by multiple genes, the effect of a functional site accounting for as little as 4.6% of the variance in liability (e.g., candidate gene 6) could be detected with the sample size studied. Covariate factors such as gender, age, environment, the disease mode (dominant, recessive, additive, or multiplicative), and penetrance, were not considered in the proposed method. Also, the power of a haplotype-based test study depends on the data structure, the trait of interest, the polymorphism information content of marker loci, the degree of LD among the markers, the allele frequencies, the amount of population stratification, and the marker density. Our approach makes the implicit assumption that the underlying population is homogeneous. In a non-admixed population, if linkage disequilibrium exists across a region, then recombination must be quite infrequent and probably can be safely ignored. If population migration, admixture, or stratification is present, this should affect the estimates of allele or haplotype frequencies and decrease the power of LD detection. A model-based procedure implementing some appropriate parameters may then improve the power of the method. Also, it remains to be determined in a future study whether our idea can be extended to other data structures, such as nuclear families with multiple sibs and/or one or both parents missing, or to extended pedigrees.

	FOOTNOTES

¹ Present address: Columbia Genome Center, Columbia University, 1150 S. Nicholas Ave., New York, NY 10032.
² Present address: Department of Biostatistics, University of North Carolina, Chapel Hill, NC 27599.

	ACKNOWLEDGMENTS

The authors thank the editor, anonymous reviewers, and Dr. Hongyu Zhao for their valuable comments and suggestions on the manuscript. Also, we thank Dr. David Clayton for allowing us to use his TRANSMIT program and providing advice on using it. The use of the GAW 12 simulated data was permitted by the Southwest Foundation for Biomedical Research, which was supported by a GAW grant, GM31575 from National Institute of General Medical Sciences. We are grateful to Dr. Jean MacCluer and Ms. Vanessa Olmo for handling our request for the GAW 12 data. This research was supported in part by National Institutes of Health (NIH) grant DA12844 (to M.D.L.); GM28356, RR03655, and DK57292 (to R.C.E.); GM58934 (to F.A.W.); National Science Foundation grant DMS-9971770 (to S.L.); and a General Clinical Research Center grant from the NIH (M01-RR00211) awarded to the University of Tennessee Health Science Center.

Manuscript received October 10, 2002; Accepted for publication March 24, 2003.

	LITERATURE CITED

TOP ABSTRACT METHODS RESULTS DISCUSSION LITERATURE CITED

ALMASY, L., J. D. TERWILLIGER, D. NIELSEN, T. D. DYER, and D. ZAYKIN et al., 2001  GAW12: simulated genome scan, sequence, and family data for a common disease. Genet. Epidemiol. 21(Suppl. 1):S332-S338.

BRODER, S. and J. C. VENTER, 2000  Sequencing the entire genomes of free-living organisms: the foundation of pharmacology in the new millennium. Annu. Rev. Pharmacol. Toxicol. 40:97-132.[Medline]

CHIANO, M. and D. CLAYTON, 1998  Fine genetic mapping using haplotypes and the missing data problem. Ann. Hum. Genet. 62:55-60.[Medline]

CHURCHILL, G. A. and R. W. DOERGE, 1994  Empirical threshold values for quantitative trait mapping. Genetics 138:963-971.[Abstract]

CLARK, A. G., 1990  Inference of haplotypes from PCR-amplified samples of diploid populations. Mol. Biol. Evol. 7:111-122.[Abstract]

CLAYTON, D. G., 1999  A generalization of the transmission/disequilibrium test for uncertain-haplotype transmission. Am. J. Hum. Genet. 65:1170-1177.[Medline]

CLAYTON, D. G. and H. JONES, 1999  Transmission/disequilibrium tests for extended marker haplotypes. Am. J. Hum. Genet. 65:1161-1169.[Medline]

COLLINS, A., C. LONJOU, and N. E. MORTON, 1999  Genetic epidemiology of single-nucleotide polymorphisms. Proc. Natl. Acad. Sci. USA 96:15173-15177.[Abstract/Free Full Text]

CORDELL, H. J. and R. C. ELSTON, 1999  Fieller's theorem and linkage disequilibrium mapping. Genet. Epidemiol. 17:237-252.[Medline]

DEMPSTER, A., N. LAIRD, and D. RUBIN, 1977  Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. B 39:1-38.

DUDBRIDGE, F., B. P. C. KOELEMAN, J. A. TODD, and D. G. CLAYTON, 2000  Unbiased application of the transmission/disequilibrium test to multilocus haplotypes. Am. J. Hum. Genet. 66:2009-2012.[Medline]

EXCOFFIER, L. and M. SLATKIN, 1995  Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population. Mol. Biol. Evol. 12:921-927.[Abstract]

EXCOFFIER, L. and M. SLATKIN, 1998  Incorporating genotypes of relatives into a test of linkage disequilibrium. Am. J. Hum. Genet. 62:171-180.[Medline]

FALK, C. T. and P. RUBINSTEIN, 1987  Haplotype relative risks: an easy reliable way to construct a proper control sample for risk calculations. Ann. Hum. Genet. 51:227-233.[Medline]

FALLIN, D. and N. J. SCHORK, 2000  Accuracy of haplotype frequency estimation for biallelic loci, via the expectation-maximization algorithm for unphased diploid genotype data. Am. J. Hum. Genet. 67:947-959.[Medline]

FALLIN, D., A. COHEN, L. ESSIOUX, I. CHUMAKOV, and M. BLUMENFELD et al., 2001  Genetic analysis of case/control data using estimated haplotype frequencies: application to APOE locus variation and Alzheimer's disease. Genome Res. 11:143-151.[Abstract/Free Full Text]

GABRIEL, S. B., S. F. SCHAFFNER, H. NGUYEN, J. M. MOORE, and J. ROY et al., 2002  The structure of haplotype blocks in the human genome. Science 296:2225-2229.[Abstract/Free Full Text]

GAO, X. and F. A. WRIGHT, 1999  Nonparametric disequilibrium mapping when haplotypes are available. Am. J. Hum. Genet. 65(Suppl.):A250.

HAWLEY, M. E. and K. K. KIDD, 1995  HAPLO: a program using the EM algorithm to estimate frequencies of multi-site haplotypes. J. Hered. 86:409-411.[Free Full Text]

HODGE, S. E., M. BOEHNKE, and M. A. SPENCE, 1999  Loss of information due to ambiguous haplotyping SNPs. Nat. Genet. 21:360-361.[Medline]

KRUGLYAK, L., 1999  Prospects for whole-genome linkage disequilibrium mapping of common disease genes. Nat. Genet. 22:139-144.[Medline]

LANDER, E. S., L. M. LINTON, B. BIRREN, C. NUSBAUM, and M. C. ZODY et al., 2001  Initial sequencing and analysis of the human genome. Nature 409:860-921.[Medline]

LAZZERONI, L. C. and K. LANGE, 1998  A conditional inference framework for extending the transmission/disequilibrium test. Hum. Hered. 48:67-81.[Medline]

LI, J., D. WANG, J. DONG, R. JIANG, and K. ZHANG et al., 2001  The power of transmission disequilibrium tests for quantitative traits. Genet. Epidemiol. 21(Suppl. 1):S632-637.

LONG, J. C., R. C. WILLIAMS, and M. URBANEK, 1995  An E-M algorithm and testing strategy for multiple-locus haplotypes. Am. J. Hum. Genet. 56:799-810.[Medline]

MARTIN, E. R., N. L. KAPLAN, and B. S. WEIR, 1997  Tests for linkage and association in nuclear families. Am. J. Hum. Genet. 61:439-448.[Medline]

MARTIN, E. R., S. A. MONKS, L. L. WARREN, and N. L. KAPLAN, 2000  A test for linkage and association in general pedigrees: the pedigree disequilibrium test. Am. J. Hum. Genet. 67:146-154.[Medline]

NIU, T., Z. S. QIN, X. XU, and J. S. LIU, 2002  Bayesian haplotype inference for multiple linked single-nucleotide polymorphisms. Am. J. Hum. Genet. 70:157-169.[Medline]

OTT, J., 1989  Statistical properties of the haplotype relative risk. Genet. Epidemiol. 6:127-130.[Medline]

PATIL, N., A. J. BERNO, D. A. HINDS, W. A. BARRETT, and J. M. DOSHI et al., 2001  Blocks of limited haplotype diversity revealed by high-resolution scanning of human chromosome 21. Science 294:1719-1723.[Abstract/Free Full Text]

QIN, Z. S., T. NIU, and J. S. LIU, 2002  Partition-ligation-expectation-maximization algorithm for haplotype inference with single-nucleotide polymorphisms. Am. J. Hum. Genet. 71:1242-1247.[Medline]

REICH, D. E., M. CARGILL, S. BOLK, J. IRELAND, and P. C. SABETI et al., 2001  Linkage disequilibrium in the human genome. Nature 411:199-204.[Medline]

RISCH, N. J., 2000  Searching for genetic determinants in the new millennium. Nature 405:847-856.[Medline]

SACHIDANANDAM, R., D. WEISSMAN, S. C. SCHMIDT, J. M. KAKOL, and L. D. STEIN et al., 2001  A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature 409:928-933.[Medline]

SELTMAN, H., K. ROEDER, and B. DEVLIN, 2001  Transmission/disequilibrium test meets measured haplotype analysis: family-based association analysis guided by evolution of haplotypes. Am. J. Hum. Genet. 68:1250-1263.[Medline]

SCHAID, D. J., 1996  General score tests for associations of genetic markers with disease using cases and their parents. Genet. Epidemiol. 13:423-449.[Medline]

SCHAID, D. J., C. M. ROWLAND, D. E. TINES, R. M. JACOBSON, and G. A. POLAND, 2002  Score tests for association between traits and haplotypes when linkage phase is ambiguous. Am. J. Hum. Genet. 70:425-434.[Medline]

SHAM, P., 1997  The transmission/disequilibrium tests for multiallelic loci. Am. J. Hum. Genet. 61:774-778.[Medline]

SHAM, P. C. and D. CURTIS, 1995  An extended transmission/disequilibrium test (TDT) for multi-allele marker loci. Ann. Hum. Genet. 59:323-336.[Medline]

SLATKIN, M. and L. EXCOFFIER, 1996  Testing for linkage disequilibrium in genotypic data using the EM algorithm. Heredity 76:377-383.

SPIELMAN, R. S. and W. J. EWENS, 1996  The TDT and other family-based tests for linkage disequilibrium and association. Am. J. Hum. Genet. 59:983-989.[Medline]

SPIELMAN, R. S. and W. J. EWENS, 1998  A sibship test for linkage in the presence of association: the sib transmission/disequilibrium test. Am. J. Hum. Genet. 62:450-458.[Medline]

SPIELMAN, R. S., R. E. MCGINNIS, and W. J. EWENS, 1993  The transmission test for linkage disequilibrium: the insulin gene and insulin-dependent diabetes mellitus (IDDM). Am. J. Hum. Genet. 52:506-516.[Medline]

STEPHENS, M., N. J. SMITH, and P. DONNELLY, 2001  A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet. 68:978-989.[Medline]

SUN, F., W. D. FLANDERS, Q. YANG, and M. J. KHOURY, 1999  Transmission disequilibrium test (TDT) when only one parent is available: the 1-TDT. Am. J. Epidemiol. 150:97-104.[Abstract/Free Full Text]

TENG, J. and N. RISCH, 1999  The relative power of family-based and case-control designs for linkage disequilibrium studies of complex human diseases. II. Individual genotyping. Genome Res. 9:234-241.[Abstract/Free Full Text]

TERWILLIGER, J. D. and J. OTT, 1992  A haplotype-based "haplotype relative risk" approach to detecting allelic associations. Hum. Hered. 42:337-346.[Medline]

THOMSON, G., 1995  Mapping disease genes: family-based association studies. Am. J. Hum. Genet. 57:487-498.[Medline]

THOMSON, G., W. P. ROBINSON, M. K. KUHNER, S. JOE, and W. KLITZ, 1989  HLA, insulin gene, and Gm associations with IDDM. Genet. Epidemiol. 6:155-160.[Medline]

TISHKOFF, S. A., A. J. PAKSTIS, G. RUANO, and K. K. KIDD, 2000  The accuracy of statistical methods for estimation of haplotype frequencies: an example from the CD4 locus. Am. J. Hum. Genet. 67:518-522.[Medline]

TOIVONEN, H. T., P. ONKAMO, K. VASKO, V. OLLIKAINEN, and P. SEVON et al., 2000  Data mining applied to linkage disequilibrium mapping. Am. J. Hum. Genet. 67:133-145.[Medline]

VENTER, J. C., M. D. ADAMS, E. W. MYERS, P. W. LI, and R. J. MURAL et al., 2001  The sequence of the Human Genome. Science 291:1304-1351.[Abstract/Free Full Text]

WILSON, S. R., 1997  On extending the transmission/disequilibrium test (TDT). Ann. Hum. Genet. 61:151-161.[Medline]

XIONG, M. and S. W. GUO, 1997  Fine-scale genetic mapping based on linkage disequilibrium: theory and applications. Am. J. Hum. Genet. 60:1513-1531.[Medline]

ZHANG, K., M. DENG, T. CHEN, M. S. WATERMAN, and F. SUN, 2002a  A dynamic programming algorithm for haplotype block partitioning. Proc. Natl. Acad. Sci. USA 99:7335-7339.[Abstract/Free Full Text]

ZHANG, K., P. CALABRESE, M. NORDBORG, and F. SUN, 2002b  Haplotype block structure and its applications to association studies: power and study designs. Am. J. Hum. Genet. 71:1386-1394.[Medline]

ZHANG, S. L., A. J. PAKSTIS, K. K. KIDD, and H. Y. ZHAO, 2001  Comparisons of two methods for haplotype reconstruction and haplotype frequency estimation from population data. Am. J. Hum. Genet. 69:906-911.[Medline]

ZHAO, H. Y., S. L. ZHANG, K. R. MERIKANGAS, M. TRIXLER, and D. B. WILDENAUER et al., 2000  Transmission/disequilibrium tests using multiple tightly linked markers. Am. J. Hum. Genet. 67:936-946.[Medline]

ZHENG, C. and R. C. ELSTON, 1999  Multipoint linkage disequilibrium mapping with particular reference to the African-American population. Genet. Epidemiol. 17:79-101.[Medline]

ZHU, X. and R. C. ELSTON, 2000  Power comparison of regression methods to test quantitative traits for association and linkage. Genet. Epidemiol. 18:322-330.[Medline]

ZHU, X. and R. C. ELSTON, 2001  Transmission/disequilibrium tests for quantitative traits. Genet. Epidemiol. 20:57-74.[Medline]

This article has been cited by other articles:

				Y. V. Sun, D. M. Jacobsen, and S. L. R. Kardia ChromoScan: a scan statistic application for identifying chromosomal regions in genomic studies Bioinformatics, December 1, 2006; 22(23): 2945 - 2947. [Abstract] [Full Text] [PDF]

				M. S. Orloff, S. K. Iyengar, C. A. Winkler, K. A. B. Goddard, R. A. Dart, T. S. Ahuja, M. Mokrzycki, W. A. Briggs, S. M. Korbet, P. L. Kimmel, et al. Variants in the Wilms' tumor gene are associated with focal segmental glomerulosclerosis in the African American population Physiol Genomics, April 14, 2005; 21(2): 212 - 221. [Abstract] [Full Text] [PDF]

				S. Lin, A. Chakravarti, and D. J. Cutler Haplotype and Missing Data Inference in Nuclear Families Genome Res., August 1, 2004; 14(8): 1624 - 1632. [Abstract] [Full Text] [PDF]