Bioinformatics Advance Access originally published online on May 16, 2008
Bioinformatics 2008 24(14):1596-1602; doi:10.1093/bioinformatics/btn236
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Penalized estimation of haplotype frequencies
1Department of Biomathematics, 2Department of Human Genetics and 3Department of Statistics, University of California, Los Angeles, CA 90095, USA
*To whom correspondence should be addressed.
| ABSTRACT |
|---|
|
|
|---|
Motivation: Low haplotype diversity and linkage disequilibrium are the rule in short genomic segments. This fact suggests that parsimony should be enforced in estimation of haplotype frequencies. The current article introduces a diversity penalty that automatically discards potential haplotypes with low explanatory power. The standard EM algorithm for haplotype frequency estimation can accommodate the penalty if one passes over to a more general minorize–maximize (MM) scheme for estimation.
Results: Our new MM algorithm converges in fewer iterations, eliminates marginal haplotypes from further consideration and reduces the computational complexity of each iteration. Estimation by the MM algorithm also improves haplotyping and genotype imputation compared to naive application of the EM algorithm. Thus, the MM algorithm is a useful substitute for the EM algorithm. Compared to the most sophisticated current methods of haplotyping and genotype imputation, the MM algorithm is slightly less accurate but at least an order of magnitude faster.
Availability: Our software will be made available in the next release the program Mendel at http://www.genetics.ucla.edu/software/.
Contact: kayers{at}ucla.edu
| 1 INTRODUCTION |
|---|
|
|
|---|
Estimation of haplotype frequencies serves a variety of purposes. For example, good estimates help distinguish ethnic groups, quantify the extent of linkage disequilibrium and guide imputation of missing genotypes. In mapping Mendelian disease genes, haplotype signatures provide evidence of unique mutation events. In association studies with common diseases, these signatures can offer more definitive predictors than single marker alleles (Akey et al., 2001; Ayers et al., 2007). With the advent of large-scale genome association studies and massive single nucleotide polymorphism (SNP) genotyping, haplotyping and associated tasks have taken on greater urgency. Fortunately, the enormous energy expended by geneticists in improving haplotyping is beginning to pay dividends in faster and more accurate software. Halperin and Eskin (2004) Scheet and Stephens (2006) and Marchini et al. (2006) summarize and compare the recent computational approaches.
The EM algorithm lying at the heart of many of these methods relies on a classical gene counting argument (Excoffier and Slatkin, 1995; Hawley and Kidd, 1995; Long et al., 1995; Qin et al., 2002). The algorithm operates on population data by filling in missing phase information based on current haplotype frequencies. Given reconstructed phases, the EM algorithm equates haplotype frequencies to imputed haplotype proportions. This iterative process of imputation and re-estimation is natural and effective. One of its strengths is that it accommodates a Dirichlet prior on haplotype frequencies (Lange, 2002). In this Bayesian context, the EM algorithm simply adds fixed pseudo-counts to imputed counts before forming its new haplotype proportions. The drawback of a Dirichlet prior is that it can only encourage the inclusion of rare haplotypes. If we want to discourage the inclusion of rare haplotypes with low explanatory power, we must turn elsewhere. In this article we propose a haplotype diversity penalty that has the desired opposite effect. Simple modification of the EM algorithm yields a novel algorithm that maximizes the penalized likelihood.
Our algorithm is example of an minorize–maximize (MM) algorithm. All EM algorithms are MM algorithms but not vice versa. Many MM algorithms, ours included, dispense with the missing data structures required by EM algorithms. In their stead one must construct a surrogate function that is optimized at each iteration. Derivation of surrogate functions requires manipulation of mathematical inequalities. Compared to the traditional EM algorithm for haplotype frequency estimation, our new MM algorithm converges in fewer iterations, eliminates marginal haplotypes from further consideration and reduces the computational complexity of each iteration. Imposition of the diversity penalty also improves haplotyping and genotype imputation. Compared to more sophisticated methods of haplotyping such as PHASE (Marchini et al., 2006; Stephens et al., 2001) and fastPHASE (Scheet and Stephens, 2006; Stephens and Scheet, 2005), the MM algorithm is slightly less accurate but considerably faster.
| 2 METHODS |
|---|
|
|
|---|
An MM algorithm for maximization involves minorizing an objective function f(p) by a surrogate function g(p|pn) anchored at the current iterate pn of a search (De Leeuw and Heiser, 1977; Groenen, 1993; Hunter and Lange, 2004; Lange, 2004). Minorization is defined by the two properties
|
| (1) |
|
| (2) |
g(p|pn) lies below the surface p
f(p) and is tangent to it at the point p=pn. Construction of the minorizing function g(p|pn) constitutes the first M of the MM algorithm.
The second M of the algorithm maximizes the surrogate g(p|pn) rather than f(p). If pn+1 denotes the maximizer of g(p|pn), then this action forces the ascent property f(pn+1)
f(pn). The proof of the property follows from the inequalities
|
|
It is instructive to derive the traditional EM algorithm for haplotype frequency estimation from the MM perspective. Let Hi be the set of maternal–paternal haplotype pairs consistent with the observed genotype of person i at each marker. If pj is the frequency of haplotype j, then the likelihood of i's observed multi-marker genotype is
|
|
|
|
is a positive constant that depends on the previous parameter vector pn but not on the current parameter vector p, and j ranges over all haplotypes consistent with at least one multi-marker genotype. A brief calculation shows that |
|
jc
lnpj+c
subject to the linear equality constraint
jpj=1 and the lower bounds pj
0. Note that the surrogate function separates parameters. This desirable feature carries over to the penalized loglikelihood.
Our diversity penalty is modeled on the lasso penalty, which was introduced in regression analysis to perform continuous model selection and enforce sparse solutions in underdetermined problems (Chen et al., 1998; Claerbout and Muir, 1973; Santosa and Symes, 1986; Taylor et al., 1979; Tibshirani, 1996). Unfortunately, the lasso penalty 
j|pj|=
is worthless in the haplotype setting because it simply reduces to the tuning parameter
. A more sensible penalty is linear for small haplotype frequencies and levels off thereafter. We therefore suggest the penalty 
jf(pj), where
|
|
. This choice of the penalty still discourages small positive estimates. The optimal value of the tuning constant
can be determined by numerical experimentation.
|
The overall minorization
|
|
, we minorize –f(q) by –q. When qn
, we minorize –f(q) by the constant –
. One can easily draw a simple graph illustrating how the tangency conditions are met in each case; see Figure 1. If Sn is the haplotype set (j: p
<
), then the surrogate function minorizing the penalized loglikelihood is |
|
In maximizing the surrogate function, the bound pj
0 can be ignored because the term c
lnpj tends to
as pj tends to 0. The equality constraint
jpj=1 must be faced, however. This is done by introducing a Lagrange multiplier
and looking for a stationary point of the Lagrangian
|
|
|
| (3) |
j pj=1 requires |
|
k
Sn pk. Substituting this result in Equation (3) produces
|
| (4) |
|
|
j
Snc
and e=
j
Snc
, then we can recast this condition as |
|
t)(d+e –
t+
) and rearranging terms leads to the quadratic
|
| (5) |
|
|
Our computer implementation of the MM algorithm in the genetic analysis program Mendel(Lange et al., 2001) simultaneously conducts haplotype frequency estimation, haplotyping and genotype imputation. Mendel uses a haplotype window surrounding a central marker flanked by f markers on the left and f markers on the right. The central marker is the object of phase and genotype imputation. The value of f is determined by the user. Mendel's current default of 9 gives a window of length 19. Estimation of haplotype frequencies commences with a defined list of haplotypes much shorter than the full list of available haplotypes. We will comment in a moment on how this list is generated and how windows at the ends of chromosomes are handled. If an individual is untyped at a marker, then during haplotype frequency estimation, all genotypes at the marker are assumed possible for the individual. Mendel takes initial haplotype frequencies to be uniform and iterates via the MM algorithm until the
1 distance
k|p
–p
| between successive iterations n and n+1 drops below 10–4 or the number of iteration exceeds 100. Once convergence is declared, Mendel discards all haplotypes with estimated frequencies below 10–8.
Given haplotype frequency estimates, Mendel imputes phase at the central marker for a given person by finding the ordered genotype at the central marker with the highest posterior probability over all consistent haplotype pairs. In the absence of pedigree data, this discovery by itself does not pin phase down. However, if we imagine sliding the haplotype window from left to right across a chromosome, then imputed ordered genotypes to the left of the central marker will be available. We can therefore assign phase to any consistent haplotype pair. If a consistent haplotype pair that agrees with the already imputed phases is not found, Mendel will search for haplotypes pairs that disagree at only one position. To the left of the central marker, a mismatch can involve either one phase switch or one allele mismatch. Allele mismatches are not allowed at the current central marker. To the right of the central marker, a mismatch can involve only allele mismatches because phase has not yet been imputed. In the rare case that a consistent haplotype pair is still not found, the most common genotype is used to fill in the missing ordered genotype.
Imputation of missing genotypes in the absence of haplotyping is handled a little differently. We now divide the consistent haplotype pairs into groups depending on the unordered genotype at the central marker. Mendel assigns a probability to each group by summing the product probabilities of its haplotype pairs. If no consistent haplotype pairs are found, then Mendel will allow for one allele mismatch. The group with highest probability determines the missing genotype at the central marker. In other words, Mendel selects the unordered genotype with highest posterior probability.
When we slide the haplotype window one marker to the right, we must construct a new abbreviated list of possible haplotypes. As we mentioned, we discard haplotypes from the existing list with estimated frequencies below 10–8. For each remaining haplotype, we crop its leftmost allele and add on its right one of the possible alleles at the new marker. If the new marker has m alleles, this action propagates each cropped haplotype into m different haplotypes in the new list. For example, the current SNP haplotype 1-2-2-1-2 is cropped to 2-2-1-2 and expanded to the two new haplotypes 2-2-1-2-1 and 2-2-1-2-2. Our retention–propagation strategy keeps all pertinent haplotypes in play. Penalization weeds outs many of the haplotypes in the new list and keeps the list from growing geometrically.
This description omits initialization of the haplotype list. Since computation times scale as the square of list length, it is imperative to adopt a strategy that minimizes list length. Thus, at the leftmost marker, we start with a window of length 1 and extend it as just described, except for imputation and cropping, until it hits length f+1. At that point, we commence haplotype and genotype imputation at the leftmost marker but still omit cropping. When the window reaches full length 2f+1, then we begin haplotype cropping. At the right end of the chromosome, haplotype propagation is omitted as soon as new markers are exhausted. Haplotype and genotype imputation continue until the rightmost marker is processed. These tactical adjustments entail more book keeping and shift the focus away from the center of the window. In compensation, they successfully keep all haplotype lists short.
In some regions little or no linkage disequilibrium exists, and the number of haplotypes can balloon out of control. Many individuals will have unique haplotypes; other individuals will be consistent with many haplotype pairs. The result is a large list of haplotypes, with many haplotypes having equal frequency estimates. In these regions, genotype imputation is already poor. To decrease computation time, we limit the number of haplotypes in a window to hmax. We order the haplotypes by decreasing frequency and find the frequency of haplotype hmax+1. Any haplotype with frequency less than or equal to this amount is dropped before moving to the next window. In very rare cases, this tactic deletes too many haplotypes, so we impose a lower limit hmin on the number of haplotypes retained. When hmax=hmin, all hmax haplotypes are kept. Setting a rigorous bound on retained haplotypes is also central to other haplotyping programs such as SNPHAP (http://www-gene.cimr.cam.ac.uk/clayton/software/).
| 3 RESULTS |
|---|
|
|
|---|
3.1 Haplotype frequency estimation
To compare the MM and the EM algorithms in haplotype frequency estimation, we randomly generated multilocus autosomal genotypes from male X chromosome haplotypes. The fathers in the 30 European (CEU) parent–offspring trios of the HapMap project are a convenient source of data (http://www.hapmap.org). We chose groups of fully typed consecutive markers outside the pseudoautosomal region with 8–11 markers per group. We then simulated 100 sets of 50, 100 and 500 genotypes from each group, sampling haplotypes with replacement. For this analysis, we began estimation with the list of all consistent haplotypes. Table 1 gives the results of applying the MM algorithm as a function of
and
for 100 genotypes. Results for 50 and 500 genotypes were similar (data not shown). The EM algorithm correspond to the choice
=0. The error column gives the average value of the
1 error
. As a rule of thumb, we suggest choosing
between 100 and 1000, with larger values for larger sample sizes. Error rates are also relatively flat in
. We recommend the choice
=0.005 for haplotype frequency estimation in a small window. To achieve the highly accurate estimates in this test, we departed from Mendel's defaults and chose the more stringent
1 convergence criterion of 10–5 and the more liberal value of 150 for the maximum number of iterations.
|
To test how the EM and MM algorithms perform in conjunction with our specific haplotype extension strategy, we simulated 100 autosomal genotypes using a longer stretch of the same HapMap X chromosome data. All 30 European haplotypes in this region of 110 consecutive markers are distinct. We initiated estimation at the first marker and extended haplotypes by adding one marker at a time. At each extension step, we computed haplotype frequencies and dropped those haplotypes with frequencies below 10–8. At most 100 haplotypes were retained at each step. Table 2 records our average results over 100 random replicates for
=.005 and
=100 and
=1000. Column 1 gives window length, Column 2 the number of iterations until convergence, Column 3 the time in seconds until convergence, Column 4 the
1 error, Column 5 the maximum number of haplotypes encountered and Column 6 the actual number of haplotypes in the sample. Sampling was done with replacement, and time is cumulative. For testing convergence, we used Mendel's default criteria.
Some interesting conclusions emerge from the table. First, standard EM does surprisingly well when paired with our extension–elimination strategy. Nonetheless, MM takes about half as many iterations and about half the time. For this increase in speed, MM pays a small but manageable price in
1 error. Error rates stabilize because there are only 30 generating haplotypes. Once all of these are included in the model, accuracy remains almost constant as new markers are added.
3.2 Genotype imputation
To compare the performance of the MM and EM algorithms in genotype imputation, we analyzed the X chromosome HapMap data on all 54 males of the African population (Yoruban). Wefirst removed pseudoautosomal markers and markers with missing genotypes. After the datawere cleaned, we constructed 30 genotypes by sampling haplotypes with replacement from the first 10 000 markers. We then randomly deleted 1% of the constructed genotypes. These steps generated 3008 missing genotypes and positioned us to evaluate the accuracy of the MM and EM algorithms in genotype imputation. In our experience, imputation by posterior probability is more accurate than imputation by most likely haplotype pair. Table 3, therefore records counts of imputation errors using posterior probabilities. Since error rates depend on the number of flanking markers, the table lists results in the range of 6–10 flanking markers. In the table, Cm denotes the number of incorrectly imputed genotypes, and Am denotes the number of incorrectly imputed alleles. These numbers differ slightly because a few imputed genotypes incorrectly specify both alleles.
|
|
The EM algorithm is overwhelmed by the sheer number of haplotypes when the number of flanking markers reaches eight. The MM algorithm discards most haplotypes and can attack much longer segments. Introducing strict limits on the number of haplotypes within a window allows the EM algorithm to recover. Haplotype frequency estimation was performed under Mendel's default convergence criterion and an upper limit of hmax=hmin=100 haplotypes per window. Compared to more stringent criteria, these choices greatly reduce computing times with virtually no effect on error rates.
Inspection of Table 3 shows that both error rates reach their approximate minima for nine flanking markers and the value
=1000. Recall that
=0 corresponds to the EM algorithm. Experiments not displayed in the table suggests that the choice
=.005 performs almost as well as our current choice
=.01. At the bottom of the table, we list the more accurate but far slower results of fastPHASE. For fastPHASE, we invoked the options –H–4 and –K10. The –H option shuts off haplotype estimation, and the –K10 options sets the number of haplotype clusters to 10. Both of these choices promote faster computation times at the expense of a slight increase in error rates.
We also compared results on two other populations with 60 individuals each from the SeattleSNPs resequencing project (http://pga.gs.washington.edu). Our initial findings on 50 different genes (data not shown) are similar to the HapMap findings. The SeattleSNPs analysis was also done in both PHASE v2.1 and fastPHASE. The software for these programs were downloaded at (http://www.stat.washington.edu/stephens/software.html). We found PHASE to have similar error rates to fastPHASE, varying by a fraction of a percent, but much longer computation times. Because of PHASE's inability to handle large numbers of markers simultaneously, we abandoned PHASE on a comparison involving 10 000 markers. In this larger dataset, the simple default of filling in missing genotypes with the most common genotype in the population results in 873 mistakes for an error rate of 29%. Mendel reduces this error rate to 4.6%, and fastPHASE reduces it further to 2.5%. In timing comparisons Mendel is about 50 times faster than fastPHASE.
3.3 Haplotyping
To compare the MM and EM algorithms on large-scale haplotyping, we reverted to the simulated data constructed from the African HapMap X chromosome data. Again we elected Mendel's default convergence criteria and set hmax=hmin=100. In this case, we filled in missing phases and missing genotypes using the ordered genotypes rather than the unordered genotypes with the highest posterior probabilities. Table 4 records the number Cm of incorrectly imputed genotypes and the number Cs of phase switch errors under this strategy. Markers with missing genotypes were not included in the switch error because an imputed genotype can differ from the true genotype. The bottom line of the table displays fastPHASE's result on the same data under the –K10 option.
|
In this comparison, Mendel's best genotype imputation error rate of 4.9% is nearly double fastPHASE's error rate of 2.6%. Mendel's best phase switch error rate of 5.4% is also about double fastPHASE's error rate of 2.8%. Increasing the number of flanking markers continues to improve Mendel's phase switch error rate, but it eventually increases the genotype imputation error rate. Haplotyping also decreases computation times because it takes advantage of the phase information to the left of the central marker. For almost every value of
, the MM algorithm outperforms the EM algorithm in phasing, genotype imputation and speed. Although fastPHASE makes fewer mistakes than Mendel, it is slower by one to two orders of magnitude, depending on parameter settings. | 4 DISCUSSION |
|---|
|
|
|---|
The EM algorithm has long served as a computational engine in haplotyping schemes. Our analysis demonstrates that penalization improves haplotype frequency estimation, genotype imputation, and haplotyping. In essence, penalization captures the parsimony nature imposes. The MM implementation of penalized estimation converges in fewer iterations than EM. Combining the MM algorithm with haplotype extension–elimination along a sliding window of markers makes it possible to handle hundreds of thousands of markers efficiently. Overall computational complexity scales linearly in the number of markers. The software described here will be made available to the public in the next release of Mendel.
Although our combination of methods does not lead to the lowest error rates in imputing missing genotypes, it is not clear that this is a serious handicap. If we accept 1% missing data as reasonable on the best genotyping platforms, then Mendel's overall error rate of 1/2000 should lead to very few incorrect inferences in association studies. The vast majority of errors committed still get one of the two alleles correct, and errors are less damaging in association studies than they are in linkage studies. This optimistic attitude should not be equated with complacency. Every source of error should be attacked.
The alternative to haplotyping via linkage disequilibrium is haplotyping via Mendelian inferences in pedigrees. When pedigree data are available, the two methods can be combined. The obvious tactic is to apply genotype elimination first marker by marker (Lange and Goradia, 1987). The partial phase information gleaned can then guide haplotype frequency estimation and genotype imputation, treating the genotyped pedigree members as unrelated. One anticipated problem with this approach is that the first stage may uncover genetic inconsistencies. In this rare circumstance, we suggest ignoring stage one and proceeding directly to haplotype frequency estimation. Neither haplotype frequency estimation nor Mendelian inferences depend on allele frequencies or map distances.
At the expense of more complex programming, there are several options for improving the speed and accuracy of the MM algorithm. For instance, one can chose window widths to reflect the local extent of linkage disequilibrium. Long windows are more compatible with strong linkage disequilibrium. We have used a fairly strict convergence criterion. Relaxing it cuts the number of iterations until convergence. There are obvious tradeoffs between speed and accuracy. Since there is no guarantee that the objective function is concave in penalized estimation, one can safeguard estimation by trying several different starting points. This tactic obviously increases computational times.
Although penalized estimation by itself does not lead to the most accurate haplotyping, it is important to stress its potential in more sophisticated schemes of haplotyping. It has much to offer in the more complex algorithms incorporated in fastPHASE. Bayesian estimation based on Markov chain Monte Carlo is less likely to be a beneficiary. Finally, it is worth mentioning that penalized estimation is apt to pay dividends in other areas of population genetics. High-dimensional estimation problems are here to stay in genetics, and one of our first reflexes in solving them should be to consider parameter regularization.
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
Funding: This research supported in part by National Institutes of Health (HG02536 to K.L.A.); United States Public Health Service (GM53275 to K.L., MH59490 to K.L.).
Conflict of Interest: none declared.
| FOOTNOTES |
|---|
Associate Editor: Martin Bishop
Received on March 29, 2008; revised on May 13, 2008; accepted on May 14, 2008
| REFERENCES |
|---|
|
|
|---|
Akey J, et al. Haplotypes vs single marker linkage disequilibrium tests: what do we gain? Eur. J. Hum. Genet (2001) 9:291–300.[CrossRef][Web of Science][Medline]
Ayers KL, et al. A dictionary model for haplotyping, genotype calling, and association testing. Genet. Epi (2007) 31:672–683.[CrossRef]
Chen SS, et al. Atomic decomposition by basis pursuit. SIAM J. Sci. Comput (1998) 20:33–61.[CrossRef]
Claerbout JF, Muir F. Robust modeling with erratic data. Geophysics (1973) 38:826–844.[CrossRef][Web of Science]
De Leeuw J, Heiser WJ. Geometric Representations of Relational Data (1977) Ann Arbor, MI: Mathesis Press.
Excoffier L, Slatkin M. Maximum likelihood estimation of molecular haplotype frequencies in a diploid population. Mol. Biol. Evol (1995) 12:921–927.[Abstract]
Groenen PJ. The Majorization Approach to Multidimensional Scaling: Some Problems and Extensions (1993) Leiden, The Netherlands: DSWO Press.
Halperin E, Eskin E. Haplotype reconstruction from genotype data using imperfect phylogeny. Bioinformatics (2004) 20:1842–1849.
Hawley ME, Kidd KK. Haplo: a program using the EM algorithm to estimate the frequencies of multi-site haplotypes. J. Hered (1995) 86:409–411.
Hunter DR, Lange K. A tutorial on MM algorithms. Am. Stat (2004) 58:30–37.
Lange K. Mathematical and Statistical Methods for Genetic Analysis (2002) New York: Spring-Verlag.
Lange K. Optimization (2004) New York: Spring-Verlag.
Lange K, Goradia TM. An algorithm for automatic genotype elimination. Am. J. Hum. Genet (1987) 40:250–256.[Web of Science][Medline]
Lange K, et al. Mendel version 4.0: a complete package for the exact genetic analysis of discrete traits in pedigree and population data sets. Am. J. Hum. Genet (2001) 69((Suppl. 4)):A1886.
Long JC, et al. An E-M algorithm and testing strategy for multiple-locus haplotypes. Am. J. Hum. Genet (1995) 56:225–232.
Marchini J, et al. A comparison of phasing algorithms for trios and unrelated individuals. Am. J. Hum. Genet (2006) 78:437–450.[CrossRef][Web of Science][Medline]
Qin ZS, et al. Partition-ligation-expectation-maximization algorithm for haplotype inference with single-nucleotide polymorphisms. Am. J. Hum. Genet (2002) 71:1242–1247.[CrossRef][Web of Science][Medline]
Santosa F, Symes WW. Linear inversion of band-limited reflection seismograms. SIAM J. Sci. Stat. Comput (1986) 7:1307–1330.[CrossRef]
Scheet P, Stephens M. A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. Am. J. Hum. Genet (2006) 78:629–644.[CrossRef][Web of Science][Medline]
Stephens M, Scheet P. Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. Am. J. Hum. Genet (2005) 76:449–462.[CrossRef][Web of Science][Medline]
Stephens M, et al. A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet (2001) 68:978–989.[CrossRef][Web of Science][Medline]
Taylor HL, et al. Deconvolution with the
1norm. Geophysics (1979) 44:39–52.[CrossRef][Web of Science]
Tibshirani R. Regression shrinkage and selection via the Lasso. JRSS-B (1996) 58:267–288.
This article has been cited by other articles:
![]() |
T. T. Wu, Y. F. Chen, T. Hastie, E. Sobel, and K. Lange Genome-wide association analysis by lasso penalized logistic regression Bioinformatics, March 15, 2009; 25(6): 714 - 721. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||





