Evaluation of methods for detecting recombination from DNA sequences: Computer simulations

The Coalescent with Recombination.

Multiple genealogies for samples of n = 10 sequences were simulated under the coalescent with recombination (18–26). In the recombination model implemented here, there are n sequences with l sites, and the population consists of N diploid individuals. A continuous time approximation is used, and the time is scaled in units of 2N generations. The recombination rate (ρ) is defined as ρ = 4Nrl, where r is the rate of recombination per site per generation.

The coalescent is built backwards in time. It is constructed by waiting for recombination or coalescent events until all ancestral sites in the n sequences have found a most recent common ancestor (not necessarily the same ancestor for all sites). The waiting times to a coalescent event are exponentially distributed with mean k(k − 1)/2 (k is the number of sequences at a given generation). The waiting times to a recombination event are also exponentially distributed, but with mean 2NrG. G is the number of potential locations where recombination can happen (20). This quantity is a number between 0 and (l − 1)k, and it depends on the outcome of the previous events. G can be written as: 1 where g_i is the number of ways a sequence can be a recombinant descendant of two sequences both ancestral to the sample. As an example, the values of g_i associated with the 4-tuples (a sequence with four sites) (1, 0, 0, 1) (1, 0, 1, 0), (1, 1, 0, 0), and (1, 0, 0, 0) are 3, 2, 1, and 0, respectively (20) (where 1 denotes that a site is ancestral and 0 that it is nonancestral). A site between two segments of ancestral material is called “trapped site;” in the tuple (1, 0, 1, 0), the first 0 is trapped, whereas the second is not (21).

Because coalescent and recombination events are independent, the time to one of these events happening is exponentially distributed with mean k(k − 1)/2 + 2NrG. The probabilities that a given event is either a coalescence or a recombination are: 2 3 Simulation of a single realization of this process is performed by starting with k = n sequences at time 0, and determining when the first event (coalescence or recombination) happens by drawing a random number u from the uniform distribution: 4 A decision is made whether this event is a coalescence or a recombination based in their relative probabilities, and then the number of sequences k is updated. If the event is a coalescence, two sequences are chosen uniformly, and their material is coalesced. The number of sequences decreases by one (k = k − 1). In the case of a recombination event, one sequence is chosen at random by assigning probabilities to the sequences based on their relative g_i values. A recombination breakpoint is then chosen uniformly over the ancestral material and the nonancestral material that is trapped between two blocks of ancestral material. When a recombination event happens, the number of sequences increases by one (k = k + 1). After each event, the value of G is updated [initially, G = (l − 1)k]. The process is continued until each site in the extant sequences has found a most recent common ancestor (MRCA). With recombination, different parts of the alignment are likely to have different coalescent trees and different times to the MRCA.

The number of recombinational events in the history of a sample of size n, R(n), has the expectation 5 (19). Not all recombination events, R(n) in total, are detectable. Regarding their effect on the genealogy, there are three types of recombination events: events that do not change the branch lengths, events that do change the branch lengths but do not change the topology, and events that change the topology. Wiuf et al. (15) provide the expectations for these three types of events.

Sequence Evolution.

Sequences were evolved on the simulated (potentially) recombinant genealogies. Several models of nucleotide substitution were used (Table 4) to study the effect of base frequency and transition/transversion ratio on the detection of recombination. Different mutation rates were used to obtain alignments with different levels of divergence (Tables 2 and 3). The expected average pairwise sequence divergence (p) depends on the parameters of the model of nucleotide substitution. A rough approximation for models with no rate variation would be 6 where θ = 4Nμl, and μ is the mutation rate per site per generation. For example, a value of θ = 100 would indicate that a randomly chosen pair of sequences is expected to differ in 0.1/(1 + 0.1)≈9% of their sites (given that sequence length, l = 1,000).

Performance Evaluation.

The recombination detection algorithms were applied to the simulated data sets, and the number of times a method inferred the presence of recombination of the 100 replicates was recorded. This number approximates the probability of detecting recombination for each method and therefore is a convenient indicator of performance. Although some methods provide a qualitative answer for the presence of recombination (yes or no), most methods calculate a P value. In the later case, recombination was inferred when the provided P value was smaller than 0.05.

Methods for Detecting Recombination.

We evaluated 14 methods for the detection of recombination (Table 5). A detailed description of these methods is published as supporting information on the PNAS web site, www.pnas.org. In general, we can tentatively classify these methods as:

(i) Distance Methods.

Distance methods look for inversions of distance patterns among the sequences (27). In general, they use a sliding window approach and the estimation of some statistic based on genetic distances among the sequences. Because the phylogeny does not need to be known, these are normally fast methods.

(ii) Phylogenetic Methods.

Several methods infer recombination when phylogenies from different parts of the genome result in discordant topologies or when orthologous genes from different species are clustered. When comparisons of adjacent sequences yield different branching patterns, there is reason to suspect the involvement of recombinational events. If the consequence of such changes results in reconciling different sequence phylogenies to a single phylogeny, then the existence of such events becomes a reasonable hypothesis (28–33). These are the methods most extensively used in the literature.

(iii) Compatibility Methods.

Compatibility methods test for partition phylogenetic incongruence in a site-by-site basis and do not require the phylogeny of the sequences analyzed to be known (34, 35).

(iv) Substitution Distribution.

Nucleotide substitution distribution methods examine the sequences for a significant clustering of substitutions or fit to an expected statistical distribution (S. A. Sawyer, http://www.math.wustl.edu/∼sawyer/geneconv/index.html; refs. 36–43).

Implementation of Methods for Detecting Recombination.

Unless otherwise noted, the number of permutations was 1,000, and the family significance level used was 0.05. Permuted alignments were obtained by randomizing the position of the columns in the alignment.

The windows program simplot (33) was generously modified by Stuart Ray (simplot's author) to implement the bootscanning (44) of every sequence in the alignment against the rest. We used a sliding window size of 200 base pairs and a step size of 10 nucleotides. Neighbor-joining trees were estimated by using F84 distances (45, 46), and bootstrap values were obtained from 100 replicates. Several bootstrapping thresholds for assignment of parenthood were explored (70, 90, and 95%), but only the 95% threshold provided reasonable false positive rates. To implement the method of Sawyer (http://www.math.wustl.edu/∼sawyer), we used a modified version of geneconv 1.81. The global permutation P values based on blast-like global scores, obtained from 10,000 replicates smaller than 0.05, were considered evidence of recombination. A multiple comparison correction is already built into these P values, so there was no need for further correction. The parameter gscale, which scales the mismatch penalty, was set to 0. To implement the homoplasy test (42), two qbasic programs written by Maynard Smith were translated into a single c program, which was benchmarked against the original implementation. Because an outgroup was not used, and to be conservative, the number of effective sites, Se, was taken to be 0.6 × the total number of sites.

The program pist (A. Rambaut and M. Worobey, http://evolve.zoo.ox.ac.uk/) was modified for simulations to implement the informative sites test (43). pist takes as input and alignment a tree and the parameter values of a model of evolution. For each data set, a maximum likelihood (ML) tree was estimated under the Hasegawa–Kinshino–Yano (HKY) + Γ model. At the same time, we obtained ML estimates of the parameters in the HKY + Γ model (π, κ, and α). The trees and model parameter estimates for each data set were used in the pist analysis. For the parametric simulation of the null distribution of the statistic (see supporting information, www.pnas.org), 100 replicates were used. A computer program was written in c implementing a modification of Maynard Smith's maximum χ² method (15, 41) by using only variable sites. The statistic is the maximum χ² in the original alignment. The P value equals the number of times the original statistic is smaller than the statistic from permuted alignments divided by the number of permutations. For all calculations, a sliding window was used, with the width of the windows set to the number of polymorphic sites divided by 1.5. This window moved in steps of one nucleotide at a time. A previous implementation calculated the P values by calculating the value of the statistic in the permuted data sets exactly at the same position (breakpoint and sequences) where the original maximum was found. This strategy resulted in many false positives and was discarded. The computer program chimaera was written in c, implementing the maximum mismatch χ² method (D.P., unpublished work) (see supporting information). The rest of the implementation is the same as in the maximum χ² method. A computer program was written in c implementing an extension of the Phylogenetic Profiles (phypro) method (27), which in its original form does not provide statistical significance. Only variable sites were used. The statistic is the minimum distance vector correlation in the original alignment. The rest of the implementation was the same as in the maximum χ² and chimaera methods.

The program plato (30) was also modified for simulations. For each data set, a maximum likelihood tree was estimated, with parameter estimates then obtained under the Hasegawa–Kinshino–Yano + Γ model of evolution. Those values were used in the calculation of the likelihoods in plato. A null distribution was simulated by 100 Monte Carlo replicates. The default window settings were used (minimum size, 5; step, 1). The windows program rdp (32) was generously modified by Darren Martin for the simulations. After exploring different conditions, the best settings were using internal and external reference sequences and a window size of 10 nucleotides. Turning off the multiple significance correction improved performance. The c program recpars (K. Fisker, ftp://ftp.daimi.aau.dk/pub/empl/kfisker/programs/RecPars) was modified for the simulations to implement the recombination parsimony method (28, 29). A recombination cost of three times the maximum substitution cost (d = 3 × s) was found to perform the best with the statistic below. The substitution costs were the true values of the parameters of the model. The statistic used was the number of histories recovered. If more than one history was recovered, recombination was inferred.

The c program reticulate (34) was modified for the simulations. The statistic used was the neighbor similarity score. A c program was written implementing the runs test (40). The program was benchmarked against results in Takahata's paper (40). Sneath's method (38, 47) implemented in qbasic was translated into c and benchmarked against the original implementation. Because P values are calculated for each pairwise comparison, P values were Bonferroni-corrected. To evaluate Stephens' method (36), an improved implementation written in fortran by Mary Kuhner (48) was translated into c and benchmarked against the original implementation. Because many tests are made, the Bonferroni correction was applied with a family α level of 0.01 (a 0.05 level gave an excess of false positives).