Parametric inference of recombination in HIV genomes

Recombination is an important event in the evolution of HIV. It affects the global spread of the pandemic as well as evolutionary escape from host immune response and from drug therapy within single patients. Comprehensive computational methods are needed for detecting recombinant sequences in large databases, and for inferring the parental sequences. We present a hidden Markov model to annotate a query sequence as a recombinant of a given set of aligned sequences. Parametric inference is used to determine all optimal annotations for all parameters of the model. We show that the inferred annotations recover most features of established hand-curated annotations. Thus, parametric analysis of the hidden Markov model is feasible for HIV full-length genomes, and it improves the detection and annotation of recombinant forms. All computational results, reference alignments, and C++ source code are available at http://bio.math.berkeley.edu/recombination/.


Introduction
Retroviral recombination is a significant contributor to genetic variation in HIV-1 genomes [29]. When an individual is infected by two or more different strains of HIV-1, recombination can yield new forms of the virus that are mosaics of the original strains. Given a particular viral genome, we would like to determine whether it was formed by recombination and, if so, what parental strains recombined to form it. This determination is important in studying the geographic epidemiology of HIV-1 [21], and in the intra-patient evolution of immune escape and of drug resistance in response to therapy [15,36]. In this paper, we present a method for inferring the parental strains of a recombinant genome. This method relies on a particular hidden Markov model (HMM) and involves statistical inference for all choices of model parameters.
Many models have been suggested for parental inference, i.e., for the identification and characterization of recombinant sequences. All rely, at some point, on a set of parameters. However, small differences in the parameter values may result in substantially different predictions by the model. This is also true of the model we present, but we will use the technique of parametric inference to ensure that we have a complete understanding of how choices for the parameter values affect predictions.
Given a multiple alignment of a recombinant viral genome with its possible parental genomes, we use a probabilistic HMM to predict, for each position in the recombinant genome, the parental strain that gave rise to it. Our HMM is motivated by the copy-choice model [3], which gives one possible biological mechanism for viral recombination. This model is based on the fact that recombination in retroviruses results from two RNA molecules being packaged in one virion [10]. In multiply infected cells, two distinct strains may be packed into a single virion [13]. In the copy-choice model, a mosaic DNA molecule results from reverse transcriptase jumping between two different RNA templates during reverse transcription in the subsequently infected cell.
We analyze our HMM over all choices of parameters using the parametric inference technique, which is a general method for evaluating graphical models. For example, parametric analysis has been used successfully for pairwise sequence alignment [8,23,35,37]. In this paper, we describe the methods of parametric inference as they are applied to our recombination HMM. To evaluate our parametric method, we use it to infer the parental subtypes of HIV-1 circulating recombinant forms (CRFs) obtained from the Los Alamos HIV Sequence Database.
We show that a simple HMM, combined with the ability to evaluate the model over its entire parameter space, is effective in predicting parental subtypes for recombinant HIV-1 genomes. We identify the range of parameter values that maximize the accuracy of our predictions on the test set, where we measure accuracy based on concurrence with hand-curated annotations. We demonstrate that parental subtype inference over all parameter values is feasible for HIV full-length genomes, and that it is much more informative than restricting to a particular choice of parameters.

Related work
There are many methods for the identification and characterization of recombinant sequences (for a current list, see http://bioinf.man.ac.uk/robertson/ recombination/). Most of these methods take as input a multiple alignment of a putative recombinant with a collection of parental sequences. These methods either output a list of parental sequences giving rise to each query position or simply state whether or not the query is a recombinant.
One method that attempts to identify the parental subtype for each position within a putative recombinant sequence was introduced in [9]. This method models each column of the input alignment by a hidden state that represents the tree topology of the column. Given fixed costs for substitutions and recombination events (changes in tree topology), the algorithm gives a most parsimonious explanation for the query sequence. This model was later formulated probabilistically as an HMM [22]. Parameter estimation and extensions for this model have been studied in depth [11,12]. These HMMs are similar to multiple alignment methods that have incorporated the concept of recombination [14,18]. Indeed, our model can be regarded as a specialization of "jumping alignments" that were introduced for remote homology detection [33].
Several methods that have been applied to HIV recombinant forms use a sliding window technique [5,20,26,32,31]. In this approach, subsequences of a certain fixed length are considered along the genome and a test statistic is calculated from each segment. These methods are local in the sense that they do not solve a global optimization problem, but detect changes of locally computed characteristics along the genome.
Methods like those just mentioned and others require many parameters, which must be estimated. Although many of these methods can be effective, their performance is highly dependent on their parameter values and the rates of evolutionary events that formed the input sequences [27]. In addition, the nucleotide substitution models for these methods are usually fixed, even though different models may be better for different data sets [28]. The framework of Bayesian statistics offers a way to deal with uncertainty in model parameters, albeit at considerable computational cost due to Markov Chain Monte Carlo (MCMC) methods that are needed for estimating the posterior probabilities [34].
We remedy parameter uncertainty by determining all solutions for all parameter choices. For the presented model, a specific parametric analysis was first conducted in [19] for the detection of recombinant sequences. Employing cost minimization, the parametric solutions were analyzed in order to identify features that indicate recombination. Here, we take a probabilistic point of view and present parametric inference methods that generalize readily to more complex probabilistic graphical models. Even though the algorithms we present are quite general and return complete information about dependence on parameter choices, they can master full-length HIV genomes of about 10,000 base pairs.

Hidden Markov model and dynamic programming
We consider sequences over the 5-letter alphabet Σ ′ = {A, C, G, T, -} of nucleotides supplemented by the gap character. For a given multiple alignment A of n columns and N + 1 DNA sequences y, s (1) , s (2) , . . . , s (N ) ∈ Σ ′n , the task is to explain the distinguished sequence y as a recombinant of the remaining sequences s (i) , i ∈ {1, . . . , N }. An annotation of y = (y 1 , . . . , y n ) is a sequence x = (x 1 , . . . , x n ) over the alphabet Σ = {1, 2, . . . , N }, where x j = i signifies that the jth character of y originates from the jth character of s (i) . We introduce an HMM for inferring unobserved annotations from observed sequences. An explanation is an annotation that maximizes the a posteriori probability of the data, given fixed model parameters. Since the "true" values for the parameters are not known with any certainty, we present a parametric analysis that yields all maximum a posteriori (MAP) annotations for all parameter values.
GFED @ABC X 1 Figure 1: Graph of the recombination hidden Markov model for n = 4 alignment columns. For column j of a multiple alignment, the observed random variable, Y j , is the character of the recombinant sequence for that column, and the hidden random variable, X j , represents the parental sequence from which that character was derived.
The HMM has hidden random variables X = (X 1 , . . . , X n ), encoding annotations in Σ n , and observed random variables Y = (Y 1 , . . . , Y n ), encoding sequences in Σ ′n . The underlying graph of the model is depicted in Figure 1.
The dynamics of the model are given by the N ×N transition matrix θ, which is the same for all transitions, and the N × 5 emission matrices θ ′ j , for 1 ≤ j ≤ n. The joint probability of an observed sequence Y = y and an annotation X = x is The likelihood of an observed sequence Y = y is obtained by marginalization, The sum can be evaluated efficiently due to the factorization which yields the forward algorithm. The factorization is valid exactly because multiplication distributes over addition, a property characteristic of semirings.
We have the following generalization.
Dynamic Programming Principle. Let S be any set equipped with two binary operations ⊕ and ⊙ such that (S, ⊕, ⊙) forms a semiring. Let π, θ, and θ ′ j be matrices over S. Then the element can be computed by the factorization For variations and generalizations of this principle, see e.g. [1,2,4,6,7,17,25]. We will revisit the Dynamic Programming Principle several times, with different semirings.
The explanations of y are exactly the annotations x satisfying This condition on x remains unchanged when passing to the logarithms of the parameters π i , θ i,i ′ , θ ′ i,j,k . Then the maximum (3) can be computed efficiently using the factorization (2) in the (max, +) algebra. This procedure is exactly the classical Viterbi algorithm. The (max, +) algebra is also known as the tropical semiring, denoted (R ∪ {−∞}, max, +), and solving (3) is equivalent to evaluating the polynomial (1) in tropical arithmetic [24].

Model specification
Motivated by the copy-choice model of recombination, we customize the general HMM by specializing the matrices θ and θ ′ j . We assume that sequences change over time by two mechanisms: mutation and recombination. Furthermore, these evolutionary events occur uniformly over the sequence with unknown probabilities µ and ρ for mutation and recombination, respectively. For simplicity, we neglect insertions and deletions in our model and encode alignment gaps with the additional character "-". Given the multiple alignment A and parameters µ and ρ, we define for i, i ′ ∈ Σ, k ∈ Σ ′ , and j = 1, . . . , n. MAP estimation in this HMM is equivalent to optimizing the following intuitive scoring scheme. Given an annotation x, we identify two salient features. The first of these is r, the number of indices j (1 ≤ j ≤ n − 1) such that x j = x j+1 , that is, the number of recombination events. The second is m, the number of j such that y j = s (xj) j , that is, the number of mutations that must have occurred given that y j originated from the sequence s (xj ) . We include in m instances where one of y j or s (xj ) i is a gap character and the other is a nucleotide. Given choices for two real-valued parameters R and M , corresponding to recombination and mutation, respectively, the annotation score of x is R · r + M · m. We seek to find the annotation x that maximizes this score. The biological meaningful solutions correspond to negative parameters R and M .
We now show that this scoring scheme is equivalent to our hidden Markov model. Given values for R and M , let Note that 1−(N −1)ρ = 1 1+(N −1)e R and 1−4µ = 1 1+4e M are the other transition and emission probabilities appearing in θ and θ ′ j . Given an annotation x with r recombination events and m mutation events, one can check that Since the denominator is constant over all annotations x, and since exponentiation is an increasing function, Prob(X = x, Y = y) is maximized exactly when the score R · r + M · m is maximized.
Conversely, given ρ and µ, let R = log(ρ) − log (1 − (N − 1)ρ) and M = log(µ) − log(1 − 4µ). Since these are simply equations (4) solved for R and M , the most probable annotation given by the HMM is exactly the annotation maximizing the score R·r+M ·m. Therefore our probabilistic HMM is equivalent to the scoring scheme described above, and results can easily be translated back and forth. For the following parametric analysis, the scoring scheme formulation is more natural.

Concurrence of explanations
We want to compare inferred explanations x with the "true" biological annotation a, which we take to be the hand-curated Los Alamos annotation that is based on phylogenetic analyses [30]. This comparison can be done for all parameter values of the model and hence allows for determining optimal parameters.
Let y be a sequence with true annotation a ∈ Σ n . When comparing inferred annotations x to the true annotation a, we consider a rating scheme that encapsulates the most important aspect of the annotation x, namely, whether it correctly states the set of sequences that have recombined to form y. We define the parental set of x as I x = {i | x j = i for some j}, i.e., the set of recombining sequences that x indicates. The concurrence of a parental set I to the true annotation a is defined as c a (I) = |I a ∩ I| |I a ∪ I| , and the concurrence of an annotation x to the true annotation is c a (x) = c a (I x ). Thus, if x correctly names the parental sequences, then I x = I a and c a (x) = 1. If x and a have no recombining sequence in common, then c a (x) = 0. In general, a fixed choice of parameters will yield several optimal parental sets. For a collection X of annotations, we average the concurrence c a over all sets I that appear as a parental set I x , for x ∈ X . To be precise, define I X = {I x | x ∈ X }, and let c a (X ) be the average of c a (I) over all I ∈ I X . Note that we use an unweighted average: the number of occurrences of a set I as a parental set I x does not affect it.
In order to rate all optimal annotations, we need to compute their parental sets. This can be achieved by another instance of the Dynamic Programming Principle. Let I X denote the collection of parental sets that corresponds to a collection X of explanations. We consider ordered pairs (φ, I) ∈ S = R ∪ {−∞} × 2 2 Σ , consisting of a number φ and a collection I of subsets of Σ. For φ 1 , φ 2 ∈ R ∪ {−∞} and I 1 , I 2 ⊂ 2 2 Σ , we define the operations which make (S, ⊕, ⊙) a semiring. If we define the matrices π, θ, and θ ′ over this semiring by setting then the object defined by expression (1) is exactly the pair (φ, I X ) with φ the optimal score and I X the corresponding collection of optimal parental sets. Thus, by virtue of the factorization (2), these collections can be computed efficiently.

Data set
We considered all of the HIV-1 full-length genomes in the Los Alamos Sequence Database. A multiple DNA sequence alignment of length 12,635 was obtained from the 2003 HIV and SIV Alignments web site (http://hiv-web.lanl.gov/). We have omitted the CRF 01 AE, because of the lacking putative parental E strain. We further excluded all recombinants that involve a CRF as one of the recombining sequences. The resulting set of 341 genomes comprises 11 different subtypes and 11 different CRFs (Table 1).
Pure subtype sequences were used to build subtype consensus sequences. The resulting consensus alignment was trimmed such that all initial and terminal gaps were removed. The subalignment covers positions 1,626 to 11,199 relative to the original alignment. These positions correspond to 1,174 and 9,144, respectively, in the HXB2 numbering scheme [16]. By trimming, we effectively excluded the flanking LTR regions, and small portions of the gag (at the 3' end) and nef (at the 5' end) genes.
Thus, the data for each of the 341 recombination inference problems consists of a multiple alignment of length 9,574 of one of the sequences and the eleven consensus sequences.

Results
We have computed, for all 341 HIV-1 genomes, the annotation polygons and all concurrence ratings for all parameter values of the HMM. The Dynamic Programming Principle has been implemented in C++, and its various occurrences have been realized using templates and operator overloading, i.e., implementation of the respective semirings. The annotation polygons can be computed in  Figure 2, the concurrence (y-axis) of the optimal parental sets with the true parental set is plotted. MAP estimates remain constant over each parameter region (numbered according to Figure 2), which gives rise to a step function. O(N 2 n 5/3 ) time [23], and computing a typical HIV annotation polygon takes less than one minute on a 2.6 Ghz Linux workstation. The complete computational results and the source code are available at http://bio.math.berkeley.edu/ recombination/.

Explanations
An annotation polygon represents the set of all possible annotations of a sequence as a recombinant of the 11 consensus sequences. For example, Figure 2 shows the annotation polygon for sequence CY.94.CY032 AF049337, which is believed to be a CRF 04 comprising subtypes A1, G, H, and K in the Los Alamos database.
The annotation polygon induces a subdivision of the parameter space into regions such that the same annotations are optimal for all parameter values in a given region. Figure 2 shows the biologically meaningful part of the polygon (middle) and the corresponding subdivision of the parameter space (right) for all non-positive M and R (i.e., where mutation and recombination is penalized). Which annotations are optimal is uniquely determined by the angle α = tan −1 R/M , and each parameter region is a cone enclosed by a pair of rays. We refer to these parameter regions by (intervals of) the angle α ∈ [0, π/2].
For our analysis of the HIV-1 genomes, an annotation is a string over Σ = {1, . . . , 11} of length n = 9,574. In Figure 4 we have picked one explanation for each vertex in Figure 2. From bottom to top, the explanations correspond to parameter choices of increasing α, that is, of increasing ratio of the probability of a mutation to the probability of a recombination. Explanations 25 and 26 use exactly the set of recombining sequences A1, G, H, and K, although some small sequence segments of the Los Alamos hand-curated annotation are missing. However, most of these small segments appear in explanations that correspond to smaller angles α.

Performance
Computing the sets of recombining sequences for each vertex of the polygon yields the concurrences for all parameter values in the corresponding normal cones ( Figure 3). Concurrences are plotted against the angle α. The resulting function takes values in [0, 1], where 0 means that the parental sets of the predicted annotations and the true annotation have no sequence in common, and 1 means that both sets coincide. Since every parameter value in a region yields the same set of explanations, the function is a step function, constant within each region. For example, for sequence CY.94.CY032 AF049337, the maximum score of 1 is attained for the two subsequent parameter regions denoted 25 and 26, which correspond to the two explanations labeled as such in Figure 4. We have computed the concurrences for all 341 genomes in order to analyze the performance of the HMM in recovering the established annotation from inferred explanations. In Figure 5 we first analyze the 271 pure subtypes.
For α = π/2 ≈ 1.57 all scores are 1, since this angle corresponds to ρ = 0, i.e., recombination is impossible. In fact, all α ≥ 1.556 yield the correct annotations for all subtypes and so have a score of 1. The average scores per subtype are displayed in the lower rightmost plot of Figure 5. For example, for α ≥ 1.472, the average of these averages is still ≥ 0.95. This average can be regarded as the expected score when assuming uniform distribution of the subtypes.
For the CRFs, concurrence is no longer an increasing function, because ρ = 0 cannot yield the true annotation ( Figure 6). For many but not all CRFs, a score of 1, i.e., perfect annotation, is reached. The maximum expected score under the assumption of uniform CRF distribution is 0.874 and it is attained for α ∈ [1.481, 1.485]. This unique local and global maximum is rather sharp indicating that the optimal value does not vary a lot between different CRFs. Scores ≥ 0.85 are expected for α ∈ [1.472, 1.494], a parameter region in which the expected score for the pure subtypes is ≥ 0.95.

Discussion
We have presented an HMM for the detection and annotation of recombinant sequences. The unifying Dynamic Programming Principle was applied for MAP estimation, parametric inference, and validation via concurrence rating. As a probabilistic model, the parameters ρ and µ of the HMM can be interpreted as the probability of observing a recombination event and a mutation event, respectively, in the data set. The HMM also allows one to compute the posterior probability of each putative parent at each site of the genome. The model can be extended in many directions, while maintaining the same conceptual and computational framework. The emission matrices θ ′ j may be parametrized to implement more flexible substitution models than the Jukes-Cantor model (which we implicitly assume), e.g., by accounting for transition/transversion bias. The model might be expanded by explicit modeling of insertions and deletions, thereby simultaneously searching for recombination and aligning the query sequence to the reference alignment. This problem can also be cast in an HMM. The corresponding higher dimensional polytopes and concurrence ratings for any of these models can be computed efficiently using the Dynamic Programming Principle. In addition, other more refined rating systems could be used to determine the accuracy of the explanations returned by the model. On the other hand, the power of our model to predict recombination, as measured by concurrence ratings, is already quite high, and we need only two parameters.
Parametric inference is very efficient for the presented model. We have used it to analyze 341 HIV genomes of size approximately 10,000 characters. Furthermore, the 2-parameter model has a concise representation of its set of explanations (Figure 4), and hence the parametric inference solution can easily be inspected. Indeed, the case of HIV-1 genomes has shown that the set of all solutions for all parameters can be much more informative than a single estimate. For example, as discussed in Section 3.1, certain parameter values may be appropriate for determining the parental set of a sequence, whereas other values may pick up shorter parental subsequences. All of this information is inherent in the annotation polygon.
Solving the parametric inference problem should be regarded as an offline precomputing step. Once its solution is available, statistical inference for fixed parameters becomes very cheap. In fact, MAP estimation in the HMM corresponds to solving a linear program on the annotation polygon. In the absence of labeled training data, i.e., annotated recombinants, the model parameters may be estimated by maximum likelihood using the Expectation Maximization algorithm. In those situations the annotation polygon can be used to assess the sensitivity of MAP estimation with respect to the uncertainty inherent in parameter estimation. The relative position of the parameter estimate to the boundaries of the discrete parameter regions determines the robustness of explanations.
Parametric analysis also provided the basis for evaluating the HMM. We used the parental sets as the feature of an explanation to be compared to the true annotation. Unlike the collection of explanations for a given choice of parameters, the collection of parental sets can be efficiently computed. For example, there are 28,704 explanations for sequence CY.94.CY032 AF049337 in parameter region 25 (Figures 2,3,4), but all share the same single parental set {A1, G, H, K}. Thus, restricting to parental sets and ignoring recombination breakpoints is a compromise between computational feasibility and accuracy. The analysis of these concurrence ratings has identified optimal parameter regions. In practice, prior knowledge about the mosaic structure of a sequence will affect the choice of the parameter region to investigate.
In summary, annotation polygons provide a powerful tool for inferring the recombinant structure of HIV genomes.