Semiparametric functional mapping of quantitative trait loci governing long-term HIV dynamics
Department of Statistics, University of Florida, Gainesville, FL 32611, USA
*To whom correspondence should be addressed.
| ABSTRACT |
|---|
|
|
|---|
Motivation: Functional mapping has proven to be powerful for characterizing quantitative trait loci (QTL) that control complex dynamic traits. More recently, functional mapping has been extended to identify the host QTL responsible for HIV dynamics by incorporating a parametric bi-exponential function for earlier stages of viral load trajectories. However, existing functional mapping cannot be used to map long-term HIV dynamics because no mathematical functions are available for later stages of HIV dynamic changes.
Results: We derived a statistical model for functional mapping of dynamic QTL through characterizing HIV load trajectories during a long-term period semiparametrically. The new model was constructed within the maximum likelihood framework and implemented with the EM-simplex algorithm. It allows for the test of differences in the genetic control of short- and long-term HIV dynamics and the characterization of the effects of viral-host genome interaction. Extensive simulation studies have been performed to test the statistical behavior of this model. The new model will provide an important tool for genetic and genomic studies of human complex diseases like HIV/AIDS and their pathological progression.
Availability: Available on request from the corresponding author.
Contact: rwu{at}stat.ufl.edu
| 1 INTRODUCTION |
|---|
|
|
|---|
The invasion of AIDS into human society is one of the most devastating events of our time. Up to 2003, the World Health Organization (WHO) has estimated 25 million person deaths due to AIDS, and 42 million person infections by the human immunodeficiency virus (HIV). However, until today is an efficient treatment to AIDS still lacking because our understanding of the genetic basis of HIV infections is very poor. Although traditional genetic research has been shown to be useful, it is too slow to decipher the complete picture of the genetic control of HIV dynamics. Part of the reason for this is due to the fact that HIV/AIDS progression is complex, in which a web of genes are involved, many displaying strong sensitivities to developmental and environmental signals in a complicated manner (Ho et al., 1995; O'Brien and Nelson, 2004).
Although AIDS itself is not a genetic disease, it is tightly associated with a host's genetic background because the host supplies HIV virus with all necessary resources to complete its life cycle and flourish. Therefore, if the underlying genetic factors can be identified and well understood, it would greatly help to build up new strategies to combat with AIDS. Currently, there are several types of treatment towards AIDS, including reverse transcriptase and protease inhibitors. It has been observed that HIV-infected individuals are heterogenous to respond to these anti-retroviral treatments. For example, some individuals respond drug treatment rapidly, while others respond much more slowly or even have no response. It is clear that such across-individual variation is at least partially determined by variants in host genes (O'Brien and Nelson, 2004). Therefore, the identification of these genes would facilitate applying the medications more efficiently and effectively, and realizing personalized drugs ultimately.
With this great promise, a natural issue arises about how to identify these underlying genes. The past decade has witnessed the rapid development of genetic mapping techniques to characterize individual genes, known as quantitative trait loci (QTL), that control variation in a complex quantitative trait (Lander and Botstein, 1989). More importantly, a new mapping strategy, called functional mapping, which is capable of characterizing QTL that regulate the pattern and form of a dynamic process, has been invented by R. Wu and colleagues (Ma et al., 2002; Wu and Lin, 2006). Functional mapping integrates the mathematical aspects of biological processes into a statistical mapping framework for QTL mapping. It attempts to estimate the genotype-specific mathematical parameters that define a curve. Functional mapping can not only provide a quantitative and testable framework for studying the interface between gene actions and development, but also increase statistical power for gene detection by reducing the number of unknown parameters to be estimated. In particular, functional mapping capitalizes on statistical modeling of the covariance structure based on autoregressive (AR) and structured antedependence (SAD) processes and, thus, makes a longitudinal high-dimensional problem computationally tractable. Because of these favorable features, functional mapping is suitable to map host QTL that affect HIV dynamics in natural human populations after anti-retroviral treatment (Wang and Wu, 2004; Wang et al., 2005; Wu et al., 2006).
The mathematical model used in HIV-related QTL mapping was a bi-exponential function derived for earlier stages of HIV viral load trajectories from fundamental biological and biochemical principles of host-viral interactions. However, from a practical viewpoint, it is more interesting to detect how genes control the dynamic change of viral loads during a long-term period. This is because genetic information about long-term HIV dynamics helps to design personalized medications for the control of HIV/AIDS progression. However, no specific parametric functions are currently derived to specify the long-term changes of HIV viral loads. In this article, we develop a statistical model for functional mapping of long-term HIV dynamics by incorporating the parametric function of earlier stages and the non-parametric function of later stages for HIV progression. This so-called semiparametric functional mapping has power to study the genetic architecture of HIV viral load trajectories in a long-term period.
| 2 METHOD |
|---|
|
|
|---|
2.1 Semiparametric mathematic model for HIV dynamics
Based on the kinetic theory of HIV replication (Ho et al., 1995; Perelson et al., 1997), Wu and Ding (1999) derived a bi-exponential model to approximate short-term ( < 12 weeks) dynamic changes of HIV virion copies in AIDS patients after initiation of highly active anti-retroviral therapy (HAART). This model is expressed as
|
| (1) |
where V(t) is the plasma HIV load at time t,
1 and
2 are two different viral decay rates in the first and second phase, representing the minimum turnover rates of two compartments, i.e. productively infected cells and latently/long-lived infected cells, respectively and P1 and P2 are related to the baseline viral loads for the two above compartments when the treatment is initiated. By estimating both viral decay rates
1 and
2, the shape of the viral load trajectory can be described and, furthermore, the antiviral effect and antiviral treatment can be quantified and assessed in clinical studies.
Although this parsimonious model enjoys both advantages of a simple computational form and easy biological interpretation, it suffers from the lack of incorporating the characteristics of long-term HIV viral load changes. That is, V(t) is a monotonically decreasing function as time increases, but during the late stage of treatment ( > 24 weeks), the HIV virion load may rebound back to high levels as observed from an actual AIDS clinical study, ACTG protocol 315 (Wu and Zhang, 2002).
In general, after the treatment is applied to AIDS patients, there are two phases of viral load decay, the early rapid decay and the late slow decay corresponding to the cleaning of free and latent viruses, respectively. It has been found that the first-phase decay, which is mainly controlled by
1 in model (1), is completed in a very short time period [usually within one or two weeks (Ho et al., 1995; Perelson et al., 1997)]. Thus,
1 is unlikely to change with the long-term treatment, which would have altered the treatment environment dramatically. This is because active free viruses are killed rather rapidly by the drug. Meanwhile, due to the latent nature of the latently/long-lived infected cells, the effect of the second phase decay, which is mainly controlled by
2, can last almost through the whole treatment period. In an extreme case, if the patient has become totally resistant to the anti-retroviral drug over a long period, the latently/long-lived infected cells may transform to productively infected cells and bring virus loads back to a high level, which leads to an extremely small or even negative value of
2. Therefore, it is not sensible to assume constant
2 over a long-term treatment period. To accommodate this characteristic of long-term HIV viral trajectories, a natural extension for the original model is to allow
2 to change as a function of time, denoted as
2(t) (Wu and Zhang, 2002). The modified model is expressed as
|
| (2) |
Equation (2) is a semiparametric model which not only preserves the interpretability of the original parametric model (1), but also provides the flexibility of a non-parametric model.
2.2 Natural cubic spline
Many non-parametric approaches can be used to estimate
2(t). As a commonly used one, natural cubic spline displays many favorable properties. For example, it can be constructed easily, has a good flexibility to fit the underlying curves of various shapes, and is the best twice-continuously differentiable interpolate for a twice-continuously differentiable function.
A cubic spline is a piecewise third-order polynomial function that passes through a set of control points. Let
(
0,
1, ... ,
m) be a non-decreasing sequence of control points, where
0 and
m are the endpoints and
1, ... ,
are the interior points. Also, let the coordinates of the control points be (
i, yi) for
. When a point
,
), define
with
. The i-th piece of the spline is then represented as
|
|
The spline function is a summation of these piecewise functions, i.e.
|
|
To ensure the second-order smoothness through all the interior points, the following conditions should be satisfied:
|
|
for
. The above constraints yield
equations but there are 4m unknown parameters. To obtain two more equations, usually some boundary conditions need to be specified; e.g. set the second derivative of each polynomial to zero at the endpoints, i.e.
|
|
This particular boundary conditions produce the so-called natural cubic spline and lead to a simple tridiagonal system that can be solved easily to give the coefficients of the polynomials (Bartels et al., 1987).
The natural cubic spline
can also be constructed via B-spline bases of degree 3 (Knott, 1999). Terminologically, the control points above are now called knots, with boundary knots of
0 and
m, and interior knots of
1, ... ,
. In general, for a B-spline of degree d with m – 1 interior knots,
basis functions are needed to span the linear space formed by all the B-spline functions. Since natural cubic spline has two additional constraints at the endpoints, the dimension of the linear space formed by natural cubic B-spline functions is
. Additionally, because a B-spline of degree d globally has (d – 1)-continuous derivatives except at the endpoints, the smoothness conditions specified on the interior points for a cubic spine are automatically satisfied by a B-spline representation. To describe a curve, another important issue for B-spline is the location of the knots, since different locations of knots yield different shapes of the spline functions. Usually, the knots are placed uniformly across the interval, which would generate a simple function representation. An alternative way is to place the knots unevenly, which could fit the curve more flexibly but with a more complex expression. In our study, the knots are placed unevenly according the sample percentiles of the design time, i.e. more knots are placed where more design points are available. In this way, the knot placement is design oriented and, thus, is more practical. For computational convenience, we will simply use R to generate the B-spline basis matrix for natural cubic spline.
2.3 The quantitative genetic model
2.3.1 Population Structure
Suppose there is a random sample of size n drawn from a natural human population at Hardy–Weinberg equilibrium. In this sample, multiple polymorphic sites, e.g. single nucleotide polymorphism (SNP), are genotyped, aimed at the identification of QTL affecting HIV dynamics. Assume that a segregating QTL with alleles A and
in the human population affects long-term viral load trajectories. The allele frequencies of A and
are expressed as q and
, respectively. For a particular genotype j (j = 0 for
, 1 for
and 2 for
) of this QTL, the parameters describing its HIV dynamics by model (2) are denoted as
, where
represents the control parameter to model
non-parametrically. The comparisons of these parameters between the three different genotypes can determine whether and how this putative QTL affects HIV dynamics.
Suppose this QTL is genetically associated with a codominant SNP marker with three genotypes
,
and
. Let p and
be the allele frequencies of M and m, respectively, and D be the coefficient of (gametic) linkage disequilibrium between the marker and QTL. According to the linkage disequilibrium-based mapping theory, the detection of significant linkage disequilibrium between the marker and QTL implies that the QTL may be linked with and, therefore, can be genetically manipulated by the marker. The marker and QTL together form four haplotypes,
, with respective frequencies expressed as
. These four haplotypes randomly unite to generate 16 combinations, some of which are collapsed to form nine distinguishable genotypes with frequencies and observed sample sizes tabulated in Table 1. If we could observe the two diplotypes for heterozygote
(
and
, summing to n5), haplotype frequencies could be estimated from these counts explicitly. However, because the two diplotypes are not observable, the EM algorithm will be used to estimate haplotype frequencies.
|
2.3.2 Linear model linking genetic and residual effects
In practice, viral loads are measured for each subject i at a finite set (Ti) of time points. We recognize that measurement times can be irregularly spaced, with different measurement schedules for all subjects. Let
|
|
where the observation vector at different time points,
, the mean vector for QTL genotype j,
and the
within-subject covariance matrix,
, are all subject-specific. Also notice that the values in
come from model (2) described by parameters
. Then, the relationship between the observation vector and genotypic mean vector can be described by
|
|
where
is an indicator variable for the QTL genotype of subject i defined as 1 if a particular QTL genotype j is indicated and 0 otherwise, and
is the residual effect vector of subject i, including the aggregate effect of polygenes and error effect, and distributed as
, where
can be fitted by some widely used variance-covariance models, such as first-order AR(1) model or SAD model [14].
2.4 Irregularly spaced structured Antedependence Model
Previously, Zhao et al. (2005) compared the SAD with AR(1) model to characterize the structure of
, and found that the former is generally more flexible than the later. For this reason, we will implement the SAD(1) model to approximate
for the accommodation of non-stationary variances and correlations. In this section, we show the structure of the irregularly spaced structured antedependence (ISSAD) model modified from the SAD(1) model.
For the SAD model, it is assumed that the residual error at a particular time t depends on the previous ones, with the degree of dependence (defined by antedependence coefficients) decaying with time lag, and the one that is formed at the current time point. The error formed at the current time point, denoted by
, is called the innovation error with
. If an error at time t is independent of all errors before
, this antedependence model is said to be of order r expressed by SAD(r). For a simple SAD(1) model, the errors at a particular time point can be expressed as
. Thus, the relationships among the residual errors, antedependence parameters and innovation errors occurring at irregularly spaced time points for subject i in an HIV/AIDS study can be expressed in matrix form as:
|
|
|
|
|
|
If the innovative variance is assumed as a constant
, the residual variance–covariance matrix of the longitudinal trait is then expressed as
|
|
where
is the innovation variance–covariance matrix expressed as
|
|
Closed forms for the inverse and determinant of matrix
can be obtained, which facilitates the estimation of the parameters,
, that model the matrix
|
|
Notice that even if innovation variances are assumed to be constant over time points, the residual variance can still change with time. Furthermore, the correlation function in ISSAD(1) is design-dependent because the correlation between the residuals of two measurements does not depend on the order of the measurements, but rather depends on the actual time interval between them.
2.5 The likelihood and estimation
Assuming the independence of n sampled subjects, the likelihood of longitudinal HIV trajectories (y) and molecular markers (M) at the underlying QTL can be constructed through a multivariate mixture model, expressed as
|
| (3) |
, and
is the vector of population genetic parameters that specify allele frequencies of the marker and QTL and their linkage disequilibrium,
, the vector for the quantitative genetic parameters that define genotype-specific HIV dynamics and
is the vector for matrix-structuring parameters with the ISSAD(1) model.
The likelihood function of Equation (3) provides a model for obtaining the maximum likelihood estimates (MLEs) of the unknown parameters (
's). This can be done by solving the likelihood equations obtained by differentiating the log-likelihood with respect to any element
in the unknown vector space
and setting the derivatives equal zero, i.e.
|
|
It is very difficult to solve the these likelihood equations in an explicit form. Hence, we implement an EM-simplex algorithm, which has been shown to be very efficient for the parameter estimation. The details for the algorithm has been given in the previous literature (Ma al., 2002; Zhao et al., 2004, 2005). In essence, for the curve parameters
and residual variance parameters
, whose derivatives are difficult to derive, the simplex algorithm is implemented. For the population genetic parameters
, the EM algorithm provides closed-form solutions and is integrated with the iterative procedure of the simplex algorithm.
| 3 HYPOTHESIS TESTING |
|---|
|
|
|---|
As shown in Wang and Wu (2004) and Wu et al. (2004), functional mapping for HIV dynamics allows for the tests of a number of biologically and clinically meaningful hypotheses. These hypothesis tests can be a global test for the existence of significant QTL, a regional test for the overall effect of QTL on a particular period of HIV dynamics, a local test for the genetic effect on viral load at a particular time point or interaction tests for the change of QTL expression across time. Below, we introduce some fundamental tests.
3.1 Global test
It tests the existence of QTL for the whole period of HIV dynamics. The hypothesis is formulated as
|
|
The
states that there is no QTL affecting HIV dynamics (the reduced model), whereas H1 proposes that such a QTL does exist (the full model). The test statistics for testing the hypotheses is the log-likelihood ratio (LR) of the reduced model to the full model. An empirical approach for determining the critical threshold is based on permutation test. The global effects of different genetic components, additive and dominant, on the shapes of entire growth curves can also be tested.
3.2 Regional test
Instead of testing the whole dynamic curve, we can additionally test how the detected QTL trigger their effects on the difference of viral load trajectories during a particular period
. This test can be based on the area under curve with the null hypothesis expressed as
|
|
3.3 Local test
It is also important to know how the detected QTL affect viral load at a given time point (t*) of interest. The null hypothesis for this test is constructed as
|
|
By moving the time point from 1 to T, this so-called local test allows for the characterization of the timing at which the QTL are switched on or off. The effects of different genetic components for a particular time point can be tested accordingly.
3.4 Interaction test
The effects of QTL may change with time, which suggests the occurrence of QTL x time interaction effects on the HIV dynamics. The differentiation of µj(t) with respect to time t represents a slope of the viral load curve (decay rate). If the slopes at a particular time point t* are different between the curves of different QTL genotypes, this means that significant QTL x time interaction occurs between this time point and next. The test for QTL x time interaction can be formulated with the null hypothesis:
|
|
The effect of QTL x time interaction on HIV pathogenesis can be examined during entire viral load trajectories.
3.5 Pleiotropic test
It tests whether a QTL pleiotropically affects the earlier and later stages of HIV dynamics. The effect of the QTL on the earlier and later stage can be tested by formulating the following null hypotheses, respectively
|
| (4) |
|
| (5) |
When both null hypotheses (4) and (5) are rejected, this means that the QTL exerts a pleiotropic effect on the two stages of HIV dynamics.
| 4 MONTE CARLO SIMULATION |
|---|
|
|
|---|
Monte Carlo simulation experiments were performed to examine the statistical properties of the model proposed for the functional mapping of long-term HIV dynamics. The experimental design for the simulation is composed of two parts, sampling strategies for patients drawn from a population and clinical trials in which patients are measured for viral load on a regular scheme. We randomly choose n = 100 or 400 subjects from a human population at Hardy–Weinberg equilibrium, and suppose one marker of two alleles M and m with allele frequency p = 0.6 for allele M. This maker is used to infer a QTL of two alleles A and a for long-term HIV dynamics based on the non-random association between marker and QTL. The allele frequency for A is assumed as q = 0.6, and the linkage disequilibrium between alleles M and A is assumed as
|
The measurements of viral load require the patients to visit the clinic on a dynamic basis. We consider a currently used design in AIDS clinical trials with the number of visits at days 0, 2, 7, 10, 14, 21, 28, 56, 84, 168 and 336. To model HIV dynamic curves with a natural cubic function, an interior knot is placed at day 21 as suggested by Wu and Zhang (2002). Therefore, we have three control parameters to describe the function of
j(t). Based on the above given parameters and design, we simulated the phenotypic and marker data by assuming the heritability (H2) of 0.1 or 0.4 at the middle measurement point.
The simulated data sets are analyzed by the proposed model. The means and mean square errors of parameter estimates were calculated from 200 simulation replicates. The results suggest that the QTL responsible for long-term viral dynamics can be detected using the marker in association with the QTL. The accuracy of curve estimates is consistently better at a higher heritability (Fig. 1a
Fig. 1b, Fig. 1c and Fig. 1d), (Fig. 1e and Fig. 1f). At each heritability level, the accuracy of curve estimates increase consistently with a sample size from 100 to 200 to 400. The curve estimation accuracy is also sensitive to the level of heritability. For example, a small sample size n = 100 with a large heritability H2 = 0.4 yields estimations as good as a large sample size n = 400 with a small heritability H2 = 0.1. This implies that in practice we need to put more efforts on the experiment design to precisely measure HIV loads, instead of simply increasing the sample size.
|
The parameters that define the curves of HIV dynamics and the structure of the covariance matrix can be well estimated with greater precision for a higher heritability (0.4) and for a larger sample size (400), as shown by the sampling errors of the MLEs (Table 2). All the population genetic parameters can be precisely estimated for the simulated data with given heritabilities and sample sizes. The estimation precision of marker allele frequencies only depends on the sample size used. The estimation precision of any parameters related to QTL segregation, QTL allele frequency and linkage disequilibrium relies on both sample size and heritability (Table 2). As expected, the estimates of the marker parameters display better precision than those of QTL-related parameters, but the latter is very responsive to sample size and heritability.
| 5 DISCUSSION |
|---|
|
|
|---|
Genetics has been thought to play an importance role in shaping the temporal change of HIV loads (O'Brien and Nelson, 2004). The underlying genetic factors may reside on the genome of the patient, like drug resistance genes, or on the genome of HIV virus, like rapidly evolved viral genes under strong selection pressure causing the HIV immune to the anti-retroviral drug, or even on both. The observation that AIDS patients respond differentially to the anti-retroviral drug will make it possible to unravel these genetic factors. In this article, we have focused on the host genome, assuming the homogeneity of HIV viral strains. This assumption has been used in several molecular genetic studies, but can be relaxed by incorporating more sophisticated mathematical and statistical models for characterizing host-virus interactions (Wang et al., 2005).
Different previous studies of parametric functional mapping for short-term HIV dynamics (Wang and Wu, 2004; Wang et al., 2005), the model derived in this article is founded on a semiparametric method that models the earlier stage of HIV dynamics by a parametric function and the later stage by a non-parametric function. It is therefore more general and useful to study the detailed genetic architecture of HIV/AIDS progression in a long-term period. Its biological advantages lie in the formulation of a number of clinically meaningful hypotheses at the interplay between gene actions and the temporal pattern of HIV load trajectories.
In this article, we used natural cubic spline to estimate the decay rate in the later stage (
2(t)). This idea has also been used by Wu and Zhang (2002), although they have never considered any genetic control of HIV dynamics. This spline function possesses several favorable properties in terms of its formulation, flexibility and interpolation. Monte Carol simulation studies were performed to investigate the statistical behavior of the semiparametric functional mapping. The results from simulation studies indicate that the new model is robust for a sample size as small as n = 100. The model is able to detect the underlying QTL for long-term HIV dynamics with reasonable precision.
We have presented a new statistical design and model for functional mapping of long-term HIV dynamics by a semiparametric approach. We assumed that one QTL is affecting HIV dynamics and associated with a marker. Although these assumptions may not be realistic, the extension of our model to include more QTL and more markers is statistically straightforward. In such a multi-locus analysis, the estimation of linkage disequilibria at different orders can gain better insights into the evolutionary pattern and process of human populations infected by HIV/AIDS. In statistics, a class of non-linear mixed-effects models can be incorporated into our semiparametric model that will make the functional mapping model more flexible. These extended models will provide an important tool for genetic and genomic studies of HIV/AIDS and its pathological progression.
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
We thank four anonymous referees for their constructive comments. The preparation of this manuscript is supported by NSF grant (0540745) to R.W.
Conflict of Interest: none declared.
| REFERENCES |
|---|
|
|
|---|
Bartels RH, et al. An Introduction to Splines for Use in Computer Graphics and Geometric Modelling, ( (1987) ) San Francisco, CA: Morgan Kaufmann, Morgan Kaufmann Publishers..
Ho DD, et al. Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature, ( (1995) ) 373, : 123–126.[CrossRef][Medline].
Knott GD. Interpolating Cubic Splines, ( (1999) ) Boston, MA: Birkhäuser..
Lander ES, Botstein D. Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps. Genetics, ( (1989) ) 121, : 185–199.
Ma CX, et al. Functional mapping of quantitative trait loci underlying the character process: a theoretical framework. Genetics, ( (2002) ) 161, : 1751–1762.
O'Brien SJ, Nelson GW, et al. Human genes that limit AIDS. Nat. Genet, ( (2004) ) 13, : 565–574..
Perelson AS, et al. Decay characteristics of HIV-1-infected compartments during combination therapy. Nature, ( (1997) ) 387, : 188–191.[CrossRef][Medline].
Wang ZH, Wu RL. A statistical model for high-resolution mapping of quantitative trait loci determining HIV dynamics. Stat. Med, ( (2004) ) 23, : 3033–3051.[CrossRef][ISI][Medline].
Wang ZH, et al. A statistical model to analyze quantitative trait locus interactions for HIV dynamics from the virus and human genomes. Stat. Med, ( (2005) ) 25, : 495–511.[CrossRef][ISI].
Wu H, Ding A. Population HIV-1 dynamics in vivo: applicable models and inferential tools for virological data from AIDS clinical trials. Biometrics, ( (1999) ) 55, : 410–418.[CrossRef][ISI][Medline].
Wu H, Zhang JT. The study of long-term HIV dynamics using semiparametric nonlinear mixed-effects models. Stat. Med, ( (2002) ) 21, : 3655–3675.[CrossRef][ISI][Medline].
Wu RL, Lin M. Functional mapping – how to map and study the genetic architecture of dynamic complex traits. Nat. Rev. Genet, ( (2006) ) 7, : 229–237.[CrossRef][ISI][Medline].
Wu RL. Functional mapping of growth quantitative trait loci using a transform-both-sides logistic model. Biometrics, ( (2004) ) 60, : 729–738.[CrossRef][ISI][Medline].
Wu S. Multilocus linkage disequilibrium mapping of epistatic quantitative trait loci that regulate HIV dynamics: a simulation approach. Stat. Med, ( (2006) ) 25, : 3826–3849.[CrossRef][ISI][Medline].
Zhao W, et al. A fast algorithm for functional mapping of complex traits. Genetics, ( (2004) ) 167, : 2133–2137.
Zhao W, et al. A nonstationary model for functional mapping of complex traits. Bioinformatics, ( (2005) ) 21, : 2469–2477.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||







