## Abstract

Recombination is a powerful evolutionary force that merges historically distinct genotypes. But the extent of recombination within many organisms is unknown, and even determining its presence within a set of homologous sequences is a difficult question. Here we develop a new statistic, Φ_{w}, that can be used to test for recombination. We show through simulation that our test can discriminate effectively between the presence and absence of recombination, even in diverse situations such as exponential growth (star-like topologies) and patterns of substitution rate correlation. A number of other tests, Max χ^{2}, NSS, a coalescent-based likelihood permutation test (from LDHat), and correlation of linkage disequilibrium (both *r*^{2} and |*D*′|) with distance, all tend to underestimate the presence of recombination under strong population growth. Moreover, both Max χ^{2} and NSS falsely infer the presence of recombination under a simple model of mutation rate correlation. Results on empirical data show that our test can be used to detect recombination between closely as well as distantly related samples, regardless of the suspected rate of recombination. The results suggest that Φ_{w} is one of the best approaches to distinguish recurrent mutation from recombination in a wide variety of circumstances.

RECOMBINATION is a fundamental biological process that can, for example, increase viral or bacterial pathogenicity by diffusing genetic material throughout populations (Awadalla 2003). The biological mechanisms of recombination differ across organisms, but in broad terms recombination results in the creation of mosaic sequences where the evolutionary history at each site may be different. Violating this tree-like assumption of evolution can lead to serious consequences when performing phylogenetic analyses for a set of sequences. Indeed, as the evolution of the sequences cannot be described by a single tree, this can lead to overestimation or underestimation of branch lengths among other problems (Schierup and Hein 2000a,b; Posada 2001; Posada and Crandall 2002). Thus, an important question for a given set of aligned sequences is to determine whether or not recombination is likely to have occurred.

The ability of a large number of general methods to detect recombination has recently been evaluated empirically and through simulation (Crandall and Templeton 1999; Brown *et al.* 2001; Posada and Crandall 2001; Wiuf *et al.* 2001; Posada 2002). These studies have established that methods such as Geneconv (Sawyer 1989), Max χ^{2} (Maynard Smith 1992), RDP (Martin and Rybicki 2000), Phypro (Weiller 1998), RecPars (Hein 1990, 1993), and neighbor similarity score (NSS) (Jakobsen and Easteal 1996) efficiently detect recombination in a wide range of circumstances (Brown *et al.* 2001; Posada and Crandall 2001; Wiuf *et al.* 2001; Posada 2002). These tests infer the presence of recombination either directly through sequence comparisons or indirectly through phylogenetic means. As no underlying assumptions are made concerning the origin of the sequences, these tests can be applied to detect recombination within any set of aligned homologous sequences. Indeed, these techniques can be used to detect recombination within either closely or distantly related genotypes (Posada 2002). Moreover, these methods can be termed general since no specific assumptions concerning sample history (beyond sequence homology) are made.

In contrast to general methods for inferring recombination, there are also population-specific methods for detecting recombination, where the samples consist of genotypes from closely related individuals. Within a single population, recombination can be tested for using nonparametric approaches such as permutation tests based on summary statistics like the correlation of linkage disequilibrium with distance (Miyashita and Langley 1988; Schaeffer and Miller 1993; Awadalla *et al.* 1999). Linkage disequilibrium is typically measured using the statistics *r*^{2} and |*D*′| (Lewontin 1964; Hill and Robertson 1968).

Recently, coalescent (Kingman 1982) methods have been developed that can specifically detect (Brown *et al.* 2001; McVean *et al.* 2002) or characterize the rate of recombination (Griffiths and Marjoram 1996; Hey and Wakeley 1997; Kuhner *et al.* 2000; Nielsen 2000; Wall 2000; Fearnhead and Donnelly 2001; Hudson 2001; McVean *et al.* 2002) for a set of samples within a single population. Recombination can be modeled under either a basic crossing-over model (Hudson 1983) or a more complex model of gene conversion (Wiuf and Hein 2000). Only a few methods (Kuhner *et al.* 2000; Fearnhead and Donnelly 2001; McVean *et al.* 2002) relax the infinite-sites model (Kimura 1969) under which a site can undergo at most a single mutation. Relaxing the infinite-sites model is important for many bacterial and viral data sets, since under the infinite-sites model, high levels of recurrent mutation can cause patterns consistent with recombination (McVean *et al.* 2002).

The basic coalescent operates under several assumptions that include constant population size, no selection, random mating, and no population structure (Hein *et al.* 2005). Whereas these assumptions can be relaxed using additional parameters such as a term for population growth (Slatkin and Hudson 1991), these additional parameters are presently not accounted for in current methods that characterize and detect recombination (Kuhner *et al.* 2000; Fearnhead and Donnelly 2001; McVean *et al.* 2002). Importantly, the influence of population structure and demographic history may adversely affect the ability of coalescent methods to correctly infer the rate of recombination (McVean *et al.* 2002; Haydon *et al.* 2004).

The myriad of methods available to detect, characterize, and find recombinant sequences is somewhat bewildering. Traditionally, general approaches have been used for recombination analysis between distantly related genotypes, whereas population genetic-based approaches have been used for recombination analysis between closely related genotypes. However, in many cases the line between the approaches is blurred, and both approaches have been used to infer the presence of recombination in bacteria, viral, and animal mitochondrial data sets (McVean *et al.* 2002; Posada 2002; Piganeau *et al.* 2004).

Often, one of the primary questions for any data analysis is to determine whether recombination is likely to be present within a set of sequences at all (Awadalla *et al.* 1999; Maynard Smith and Smith 2002; McVean *et al.* 2002; Posada 2002; Piganeau *et al.* 2004; Tsaousis *et al.* 2005). Indeed, there are still open questions with regard to the extent of recombination in animal mitochondrial DNA (Maynard Smith and Smith 2002; Piganeau *et al.* 2004; Tsaousis *et al.* 2005). Moreover, if the sequences are obtained from closely related, yet distinct, organisms or from many different populations, it is inappropriate to analyze the sequences in a framework that assumes a single population, such as linkage disequilibrium or coalescent approaches (Tsaousis *et al.* 2005). But determining whether recombination has occurred in such circumstances is an important question that cannot be easily answered in a parametric framework. A robust nonparametric test for recombination can help distinguish between the presence and absence of recombination in such cases.

Testing for recombination can statistically validate visual evidence of recombination obtained using, for instance, phylogenetic network approaches (*e.g.*, Huson and Bryant 2006) or independently verify the presence of recombination if a positive estimate of the rate of recombination is inferred (*e.g.*, McVean *et al.* 2002). Moreover, it is often difficult to distinguish between rate heterogeneity and recombination in many circumstances (Grassly and Holmes 1997; McGuire and Wright 2000) and thus regions that exhibit phylogenetic inconsistencies can be individually tested for recombination. Additionally, testing for recombination can be used as a prior probability for the presence of recombination when inferring the points at which infrequent recombination may have occurred (Minin *et al.* 2005). In this sense, testing for recombination can be used in conjunction with other methods.

Ideally, a single test could correctly determine whether recombination is present within any given set of aligned sequences, regardless of population history, demographic history, recombination rate, or mutation rate. Preferably, such a test would also minimize the production of false positives. Here we develop a new test that is powerful under many of these different situations and produces few false positives. Through simulation and empirical data analysis we characterize the performance of our test under various rates of recombination, rates of mutation, demographic histories, and sample sizes. We also show through simulation that a simple model of substitution rate autocorrelation (consistent with mutational “hot spots”) gives rise to a signal similar to recombination for two different general tests, Max χ^{2} and NSS, but not for our method.

## METHODS

Tests for recombination based on the principle of compatibility have proved to be among the most powerful (Brown *et al.* 2001; Posada and Crandall 2001; Wiuf *et al.* 2001; Posada 2002). The traditional binary notion of compatibility (Le Quesne 1969) is well suited for sites with at most two alleles, but can be directly extended into a broader notion (Penny and Hendy 1986) that we term here as refined incompatibility. We then develop a new statistic to test for recombination, the Φ_{w}- (or pairwise homoplasy index, PHI) statistic that uses this notion of refined incompatibility.

#### Compatibility and incompatibility:

It is not obvious how to determine the genealogical history of a single site. As such, the pattern of mutation present at multiple sites must be used to infer the genealogy of the sample as a whole. One possibility is to use the observed patterns at pairs of sites, in particular the notion of compatibility (Le Quesne 1969) or the “four-gametes” test (Hudson and Kaplan 1985). Two sites *i* and *j* are compatible if and only if there is a genealogical history that can be inferred parsimoniously that does not involve any recurrent or convergent mutations (known as homoplasies as in Figure 1b). If the two sites are not compatible, they are termed incompatible. Under an infinite-sites model (Kimura 1969) of sequence evolution, the possibility of a homoplasy does not exist, and so incompatibility for a pair of sites implies that at least one recombination event must have occurred, as in Figure 1a. This can be used to estimate the minimum number of recombination events present in the sample as a whole (Hudson and Kaplan 1985; Song and Hein 1999; Myers and Griffiths 2003). Testing for compatibility can be accomplished by checking if all four combinations of {00, 01, 10, 11} are present among the sequences (Le Quesne 1969).

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

The traditional, binary notion of either compatibility or incompatibility treats a single homoplasy the same as many homoplasies. That is, although in some situations more than one homoplasy can be parsimoniously inferred for a pair of sites (Camin and Sokal 1965; Penny and Hendy 1986), this information is disregarded. Consider two sites *i* and *j*, with |χ_{i}| and |χ_{j}| representing the number of observed states (alleles) at each site. Let *l*(χ_{i}, χ_{j}) denote the minimum number of mutations required by *any tree* used to represent the genealogical history of both sites. Thus *l*(χ_{i}, χ_{j}) represents the maximum parsimony score for these two characters over all trees. Note that *l*(χ_{i}, χ_{j}) ≥ (|χ_{i}| − 1) + (|χ_{j}| − 1) as each state (except the ancestral state) must arise at least once in the tree. Define the refined incompatibility score of sites *i* and *j* asThe refined incompatibility score relates to the traditional notion of compatibility in the following way: two sites are compatible if and only if *i*(χ_{i}, χ_{j}) = 0; if *i*(χ_{i}, χ_{j}) > 0 the two sites are incompatible. There are also two interpretations of this refined incompatibility score: in the absence of recombination, this score represents the minimum number of homoplasies that have occurred in the history of the samples for these two sites (Penny and Hendy 1986); in the absence of recurrent or convergent mutations, this score represents the minimum number of recombinations that have occurred between the two sites (T. Bruen and D. Bryant, unpublished data). This latter result depends on viewing recombinations as unrooted subtree-prune and regraft operations (see Hein *et al.* 2005). Importantly, this score can be calculated quickly [linear time in the number of sequences (Bruen and Bryant 2006)], which allows alignments with large numbers of sequences to be evaluated rapidly.

A parsimony informative site has at least two different alleles that are represented by at least two different sequences each (there must be at least four sequences at a site for the site to be parsimony informative) (Felsenstein 2004). A compatibility matrix (Sneath *et al.* 1975; Jakobsen and Easteal 1996) is traditionally used to represent compatibility between all pairs of parsimony informative sites. This matrix can also easily be extended into a refined incompatibility matrix by setting each entry (*i*, *j*) equal to the refined incompatibility score between any two sites *i* and *j*.

Sites that have the same history will tend to be more compatible than sites that have different histories (Sneath *et al.* 1975; Jakobsen and Easteal 1996; Drouin *et al.* 1999). One way to measure the extent of “clustering” in the matrix is to consider the proportion of neighboring cells in the matrix that are either compatible or incompatible. The resulting statistic is termed the NSS and has been used as a powerful test for recombination (Jakobsen and Easteal 1996; Brown *et al.* 2001; Posada and Crandall 2001; Wiuf *et al.* 2001; Posada 2002). However, simulations suggest that the NSS produces an excess of “false positives” in certain situations (see results and discussion) and so we have developed an alternative statistic.

#### Test statistic (Φ_{w}):

The degree of genealogical correlation between neighboring sites is negatively correlated with the rate of recombination (Hudson and Kaplan 1985). In the case of finite levels of recombination, the genealogical correlation of sites is partially reflected by a tendency of closely linked sites to have greater compatibility than distant sites (Hagenblad and Nordborg 2002; Innan and Nordborg 2002).

To measure the similarity between closely linked sites, we propose calculating a new statistic, the pairwise homoplasy index (PHI). The idea is to calculate the mean refined incompatibility score from nearby sites by using the first *k* off-diagonal rows of a refined incompatibility matrix (see Figure 2). Let *w* denote a fixed width (measured in bases) and choose *k* so that it is proportional to *w*. Specifically, let *q* denote the proportion of parsimony informative sites within the alignment and set *k* = *wq*. The statistic thus measures the mean refined incompatibility score of sites up to (approximately) *w* bases apart. We can now formally define the Φ or PHI statistic asThe term “pairwise homoplasy index” refers to the fact that the refined incompatibility score can be interpreted as the minimum number of convergent or recurrent mutations (homoplasies) necessarily present on any tree describing the history of any two sites *i* and *j*. The term *k*(2*n* − *k* − 1)/2 is a normalizing factor.

Clearly *w* should be somewhat less than the total number of sites but large enough that a number of comparisons are made. For all simulated and empirical analyses *w* was set to 100 and *k* chosen according to the above formula. Other choices of *w* were also considered (*w* = 50 and *w* = 150), but simulations (across different sequence lengths) suggested that *w* = 100 was slightly better than the other two choices (results not shown).

#### Significance:

Significance of the observed Φ_{w}-statistic can be obtained by using a permutation test. Under the null hypothesis of no recombination, the genealogical correlation of adjacent sites is invariant to permutations of the sites as all sites have the same history. But in the case of finite levels of recombination, the order of the sites is important, as distant sites will tend to have less genealogical correlation than adjacent sites. Let denote the observed value of the Φ_{w}-statistic on the original alignment and let *Z*_{0} denote the value of the Φ_{w}-statistic for a random permutation of the sites. Hence *Z*_{0} is distributed according to the null hypothesis of no recombination. To determine the significance of the observed value , a Monte Carlo *P*-value can be directly estimated by permuting the alignment many times and counting the proportion of times the Φ_{w}-statistic on a permuted alignment is less than or equal to . However, computation of *P*-values based on permutations of the alignment is time consuming. One way to circumvent this problem is to determine the distribution of the test statistic under permutations of the alignment. The expectation (*E*_{0}(Φ_{w}) = μ′) and variance (Var_{0}(Φ_{w}) = σ^{2}) of Φ_{w} can be calculated analytically (see appendix a for details). Moreover, initial simulations indicated that the distribution of Φ_{w} under permutations of the alignment is approximately normal (results not shown). Using these assumptions, the value of can be calculated aswhere *n*(τ | μ′, σ^{2}) denotes a normal probability distribution function with mean μ′ and variance σ^{2}. This alternative to the permutation test has the advantage that it can be obtained quickly and gives a more precise *P*-value under an assumption of normality.

The normality of the distribution of the test statistic can be explained by noting that for a large refined incompatibility matrix, calculating the Φ_{w}-statistic amounts to taking the mean of a small sample of values from the matrix. The simplest version of the central limit theorem then suggests that taking the mean of a small sample within a “large” matrix has a limiting normal distribution, if the terms are independent and identically distributed (Casella and Berger 2001). However, in this case the central limit theorem provides a guide rather than a formal equivalence.

For every data set examined (both simulated and empirical) the significance of the observed Φ_{w}-statistic was calculated using the permutation test directly as well as the normal alternative. The *P*-values obtained by using the permutation test are written as *P*_{P}(Φ_{w}) whereas the *P*-values obtained by using the normal alternative are written as *P*_{N}(Φ_{w}).

#### Simulation study:

We repeated many of the same simulations that had been performed in other studies (Posada and Crandall 2001; Wiuf *et al.* 2001) but expanded the parameter search space and considered the Φ_{w}-statistic as well as additional tests. The protocol followed was based on simulations from the neutral coalescent model (Kingman 1982) with recombination (Hudson 1983).

The coalescent model provides a natural foundation for simulation (Crandall and Templeton 1999; Brown *et al.* 2001; Posada and Crandall 2001; Wiuf *et al.* 2001). Simulations were almost all conducted using the program Treevolve (Grassly *et al.* 1999). For very high rates of recombination (ρ = 128), simulations were performed using the program Hudson (Schierup and Hein 2000a,b) since the program Treevolve did not run at such high rates of recombination. Mutations were added according to a Jukes–Cantor model (Jukes and Cantor 1969). Other methods of sequence evolution were also examined, including the addition of extreme rate heterogeneity (α = 0.1), which resulted in a moderate decrease in power for all methods (results not shown). For each parameter setting, 1000 replicate data sets were created, with each replicate consisting of an alignment of length 1000 (see appendix b for further details). Significance was set at the 0.05 level.

In addition to the Φ_{w}-statistic, four of the best nonparametric tests were computed for each parameter setting, namely the Max χ^{2}-statistic (Maynard Smith 1992), the NSS (Jakobsen and Easteal 1996), and two measures of correlation of linkage disequilibrium (*r*^{2} and |*D*′|) with distance (Lewontin 1964; Hill and Robertson 1968; Miyashita and Langley 1988; Schaeffer and Miller 1993). Furthermore, results obtained from a coalescent-based likelihood permutation test (LPT) from LDHat (McVean *et al.* 2002) are reported as well. The Max χ^{2}-statistic has been found to be the best general test for detecting recombination in a recent empirical study (Posada 2002), and the NSS statistic has been found to be very efficient as well (Brown *et al.* 2001; Posada and Crandall 2001; Wiuf *et al.* 2001; Posada 2002). Correlation of linkage disequilibrium with distance using *r*^{2} has been found to be the strongest nonparametric approach for detecting recombination within populations (McVean *et al.* 2002). Recently, the likelihood permutation test was introduced as a powerful alternative to methods based on linkage disequilibrium (McVean *et al.* 2002). For the Max χ^{2}-statistic a fixed window size of the number of polymorphic sites divided by 1.5 was used following a previously described protocol (Posada and Crandall 2001; Posada 2002). For both measures of correlation of *r*^{2} and *D*′ with distance, only sites with two alleles segregating and minor allele frequencies of at least 0.1 were used, as this approach tends to maximize power (Weir and Hill 1986; McVean *et al.* 2002). For the likelihood permutation test, precomputed likelihood files were used on the basis of 101 grid points with a value of θ per site of either 0.001 or 0.1. For each replicate, if the expected mean sequence diversity was <10%, then a likelihood file with a θ per site value of 0.001 was used; otherwise a likelihood file with a θ per site value of 0.1 was used (under a constant-size population the expected mean sequence diversity of 10% corresponds to an expected value of θ per site of ∼0.12). The significance for each of the statistics was obtained using a permutation test. For the power determination, 1000 permutations were performed, whereas for the false positives, 200 permutations were performed.

##### Power:

To determine power in the presence of recombination, the recombination rate ρ (under population growth ρ^{†}) varied among 0, 1, 2, 4, 8, 16, and 128; the expected nucleotide diversity *p* between any two sequences varied among 1, 5, 10, 15, and 25%; and the growth rate of the population β varied between 0 (constant-size populations) and 5000. The sample size *m* varied among 5, 10, 15, 25, and 50. For ρ = 128 simulations with β = 5000 were not performed since this option was not available with the program Hudson. More details explaining the protocol can be found in appendix b and elsewhere (Wiuf *et al.* 2001).

##### False positives:

Substitution rate heterogeneity across sites on a genealogy was modeled here using a Γ-distribution (Uzzell and Corbin 1971; Yang 1993). In this case, the substitution rate at each site *i*, *Z _{i}*, is drawn from a Γ-distribution with shape parameter α and scale parameter 1/α (Yang 1993).

Autocorrelation among substitution rates was modeled assuming Markov dependence among rates (Yang 1995). To achieve this, two random variables *Y _{i}* and

*Y*

_{i}_{+1}were drawn from a bivariate normal distribution with correlation ρ

_{N}and transformed into two marginally distributed gamma random variables

*Z*and

_{i}*Z*

_{i}_{+1}with correlation ρ

_{G}(Yang 1995). Using the bivariate normal distribution of

*Y*and

_{i}*Y*

_{i}_{+1}(including correlation ρ

_{N}), the probability distribution function of random variable

*Y*

_{i}_{+1}was obtained conditional on the random variable

*Y*, allowing Markov-dependent substitution rates to be drawn. The substitution rates

_{i}*Z*and

_{i}*Z*

_{i}_{+1}then represent draws from a bivariate Γ-distribution with correlation ρ

_{G}. The value of ρ

_{G}is positively correlated with the value ρ

_{N}but not identical (Yang 1995).

Data sets were simulated using a modified version of Treevolve (Grassly *et al.* 1999) with a number of the sampling functions taken from PAML (Yang 1997). The correlation parameter ρ_{N} varied among 0 (no correlation), 0.3, 0.6, and 0.9; the expected nucleotide diversity *p* between any two sequences varied among 1, 5, 10, 15, and 25%; the value of α for the Γ-distribution varied among 0.1, 1.0, and ∞; and the growth rate of the population β varied between 0 (constant-size populations) and 5000. The sample size *m* varied among 5, 10, 15, 25, and 50.

#### Empirical data:

A number of population and species level data sets were examined. The presence of recombination in each of these data sets was debated, unknown, or suspected. The rate of recombination in these data sets ranged from rare to very frequent. In general, data sets with at least a few hundred sites were chosen.

Tests for recombination were performed using the Φ_{w}-statistic as well as the Max χ^{2}-statistic (Maynard Smith 1992) and the NSS statistic (Jakobsen and Easteal 1996). As in the simulation studies, *w* was set to 100 for all analyses. One thousand permutations were performed to obtain significance. Additional results are reported for the population level data sets, using permutation tests based on *r*^{2} and |*D*′| (Lewontin 1964; Hill and Robertson 1968; Miyashita and Langley 1988; Schaeffer and Miller 1993) as well as a coalescent-based LPT with LDHat (McVean *et al.* 2002). Furthermore, an estimate of the rate of recombination was also obtained in LDHat using a model of crossing over rather than gene conversion. The maximum value of ρ was set to 100 and 100 grid points were used in LDHat. The value of Tajima's *D*-statistic is also reported, as it can be an indicator of population growth or selective pressure (Tajima 1989). Table 1 summarizes the data sets used. The data sets include sequences from bacteria, viruses, and fungi. Two of the data sets were from animal mitochondrial DNA (mtDNA).

For the Boletales data set additional analysis was performed by first estimating a neighbor-joining tree (Saitou and Nei 1987) using PAUP* (Swofford 1998). Branch lengths for the tree, a transition/transversion ratio, codon frequencies, a value of α for the substitution rate heterogeneity (Yang 1993), as well as the degree of substitution rate autocorrelation (estimated using the autodiscrete gamma model) (Yang 1995), were then estimated using a codon model in PAML (Yang 1997). A parametric bootstrap of 1000 replicates was then performed under the estimated parameters using a modified version of PAML that allowed autocorrelated substitution rates. For each replicate, a test for recombination was performed using the Max χ^{2}-statistic, the NSS statistic, and the Φ_{w}-statistic (with 1000 permutations). Significance was set at 0.05.

## RESULTS AND DISCUSSION

#### Simulation studies:

##### Analytical calculation of P-values:

Table 2 shows the proportion of times that recombination was inferred using Φ_{w}, when the rate of recombination ρ was set to 0 and there was no population growth (β = 0). Since the significance level was set to 0.05, the Φ_{w}-test is too conservative when the mean sequence diversity is ∼1% or when there are few samples (*e.g.*, *m* = 5). This is partly due to the fact that there are very few informative sites or incompatibilities produced in these situations (results not shown). Table 2 also indicates that when the sequence diversity and sample size are small, obtaining significance using the permutation test (*P*_{P}(Φ_{w})) is even more conservative than obtaining significance using the normal distribution (*P*_{N}(Φ_{w})). On the other hand, Figure 3 shows that both methods for obtaining significance give very similar answers for higher amounts of sequence diversity (at least 10%), with at least 15 samples. These results suggest that it is sufficient to obtain significance for Φ_{w} using the normal distribution. For all subsequent simulations, the results quickly obtained with the Φ_{w}-statistic using the normal distribution are reported.

##### Time:

The time to calculate Φ_{w} is much faster than other population genetic methods especially for moderate numbers of sites and sequences. For instance, several simulated alignments of 25 samples with 5000 sites with moderate sequence diversity (10%), corresponding to viral genomic samples, were analyzed on a Mac G4 desktop computer. The time taken to analyze each alignment was ∼20 sec using Φ_{w} without the permutation test, 30 sec using Φ_{w} with the permutation test, 7 min with the linkage disequilibrium methods (using LDHat), and 8 hr using the likelihood permutation test of LDHat (using a precomputed likelihood file). For longer alignments, however, the permutation test becomes impractical even for Φ_{w} and in these cases analytical *P*-values are the only way to practically test for recombination. It is worth noting that since the power to detect recombination increases as a function of sequence length (Wiuf *et al.* 2001), this constitutes an important advantage for the Φ_{w}-test, since faint recombinant signals may be detectable using only very long sequences.

##### Power:

Figure 4 shows the power to detect recombination for Φ_{w}, Max χ^{2}, NSS, the LPT in LDHat, and two measures of correlation of linkage disequilibrium with distance (*r*^{2} and |*D*′|), when the rate of recombination ρ is greater than zero, for two different sample sizes (*m* = 10 and *m* = 50). Two principal types of genealogies were created: with and without population growth. If there is population growth, the genealogies created will be more star-like with long branches at the leaves (Griffiths and Tavaré 1998; Wiuf *et al.* 2001). If there is no population growth, there are short branches at the tip but long branches at the root. When genealogies are more star-like, recurrent mutations will tend to mask the initial recombination, and the recombination events are best considered to be “ancestral.”

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

The top rows of Figure 4, a and b, show that without population growth (β = 0), all six methods performed similarly, although overall Φ_{w} is the most powerful method with a large number of samples. Without population growth, the power to detect recombination of all six methods generally increases as a function of both sequence diversity and the rate of recombination, similar to earlier observations (Posada and Crandall 2001; Wiuf *et al.* 2001). A notable exception is the LPT for which there is a slight decline in power when the mean sequence diversity reaches 10%. At this point, a likelihood file with a value of θ per site of 0.1 was used rather than a likelihood file with a value of θ per site of 0.001. However, when the sequence diversity reaches 10%, the expected value of θ per site is ∼0.12, suggesting that a value of θ per site of 0.1 is a better choice. Nonetheless, more power may be obtained by using a gross underestimate of θ, although previous work has demonstrated a relative insensitivity of the LPT to a specific estimate of θ (McVean *et al.* 2002).

The top rows of Figure 4, a and b, suggest that the Φ_{w} method performs similarly to the linkage disequilibrium approaches when there is very little sequence diversity (*e.g.*, *p* = 1%), despite the fact that the test is too conservative in these circumstances (Table 2). For very little sequence diversity (*i.e.*, *p* = 1%), the coalescent-based method LPT is the most powerful method in constant-size populations, but has about the same power as Φ_{w} for growing populations. However, the results suggest that all methods may underestimate the presence of recombination if few sequences are present with very little divergence, especially in an expanding population (or “star-like” genealogy).

By comparing the bottom rows of Figure 4, a and b, to the top rows of Figure 4, a and b, it is evident that detecting the presence of recombination under population growth (β = 5000) is a more difficult task than detecting the presence of recombination without population growth (β = 0). Of all six methods, the bottom rows of Figure 4, a and b, suggest that Φ_{w} is much better at detecting recombination under population growth than Max χ^{2}, NSS, the coalescent-based LPT, or the linkage disequilibrium approaches. For the coalescent-based LPT, it is worth noting that population growth could be incorporated in the method in the future, possibly increasing power. The decline of linkage disequilibrium in expanding populations using *r*^{2} is consistent with previous observations (Slatkin 1994; McVean 2002), but the results suggest that the performance of the |*D*′| statistic is similar. The results for the Φ_{w}-test suggest that subsequent mutations do not “mask” the recombinant signal for this method. Interestingly, this is similar behavior to the RECPARS method (Hein 1993; Wiuf *et al.* 2001) and may be of particular importance when trying to determine ancestral recombination between diverged genotypes. The results also suggest that the Φ_{w}-statistic can be used to distinguish between star-like genealogies due to population growth and star-like genealogies due to recombination (Schierup and Hein 2000b).

A comparison of the top row of Figure 4a to the top row of Figure 4b reveals that an increase in sample size from *m* = 10 to *m* = 50 causes an increase in the ability of all six methods to infer recombination when there is no population growth (β = 0). For population growth (the bottom rows of Figure 4, a and b), the power to detect recombination for the NSS statistic for actually decreases sharply from *m* = 10 to *m* = 50. But for the other five tests, the power to detect recombination generally increases when moving from *m* = 10 to *m* = 50 even under population growth. These results expand upon some previous observations (Wiuf *et al.* 2001).

Under a neutral coalescent model with recombination, it is possible to use a likelihood-ratio test to determine whether the hypothesis of no recombination (ρ = 0) should be rejected at a given significance level (Kuhner *et al.* 2000; Brown *et al.* 2001). However, even when data are simulated according to the neutral coalescent with low levels of recombination, the hypothesis ρ = 0 is rejected only a limited proportion of the time (Brown *et al.* 2001). However, such a simulation represents an ideal situation, where the likelihood-ratio test is guaranteed to be the most powerful (Brown *et al.* 2001) and the model used to infer ρ is identical to the model used to generate samples. This suggests that it might be difficult for any test to correctly infer the presence of recombination for very low recombination rates. Additionally, a theoretical analysis shows that generating small sets of samples using a low rate of recombination produces only a limited number of incompatibilities (Wiuf *et al.* 2001). It is thus possible that full-likelihood approaches (Kuhner *et al.* 2000; Fearnhead and Donnelly 2001) or a phylogenetic network (Huson and Bryant 2006) approach could be particularly useful to determine whether there is any possibility of recombination when only a weak recombinant signal exists.

Table 3 demonstrates that Φ_{w} can detect recombination even under extremely high recombination rates (ρ = 128). Except for low sequence diversity (*p* = 1%), the presence of recombination is correctly inferred each time. But even for low sequence diversity, the presence of recombination can be inferred nearly every time by increasing the sample size from *m* = 10 to *m* = 50.

It is worth noting that the Φ_{w}-statistic can also be calculated without the refined incompatibility score, but using only the traditional notion of compatibility. For cases without population growth (β = 0), the results are almost identical (results not shown). On the other hand, with population growth (β = 5000), there is an increase in power using the refined incompatibility score when the number of samples is large (*e.g.*, *m* = 50) and there is some recurrent mutation. For a rate of recombination of ρ = 1, a sample size of 50, and exponential growth, the gains in power using the refined incompatibility score rather than the compatibility score were 2, 5, and 12% for mean pairwise sequence divergences of 10, 15, and 25%, respectively. Similar results are obtained for ρ = 2 but not for higher rates of recombination (results not shown). This suggests that the refined incompatibility score is a useful extension to the traditional notion of compatibility especially for large sample sizes with sites that experience recurrent mutations.

For no population growth, the Φ_{w}-test and the linkage disequilibrium approaches perform similarly, although Φ_{w} is more powerful for a large number of samples. However, Φ_{w} is applicable even if the samples are from different species or different populations, whereas the linkage disequilibrium and coalescent approaches are not (Tsaousis *et al.* 2005). Under population growth, however (β = 5000), only Φ_{w} continues to consistently infer the presence of recombination as the power of the other five methods suffers sharp declines. This suggests that, of all six methods, Φ_{w} has the greatest flexibility in detecting recombination in the different circumstances studied.

##### False positives:

Of particular concern for any test for recombination is the effect of confounding processes such as substitution rate heterogeneity and autocorrelated substitution rates. Autocorrelation of substitution rates implies that the rate of substitution of one site is not independent of the rate of substitution of a neighboring site and can create “mutational hot spots” within a sequence. This can potentially create the same patterns as recombination.

Figure 5 shows the proportion of false positives for Max χ^{2} and NSS when there is no recombination (ρ = 0) but “mosaic” sequences are artificially induced by using a range of autocorrelated substitution rates. Figure 5 shows that both Max χ^{2} and NSS falsely infer the presence of recombination >50% of the time in certain cases. The results for the linkage disequilibrium, likelihood permutation test, and Φ_{w} are omitted from Figure 5 since these methods did not falsely infer recombination >7% of the time, although Table 4 shows this information for Φ_{w}. Table 4 shows that the Φ_{w}-statistic did not infer recombination >6% of the time when recombination was falsely inferred >50% of the time using both Max χ^{2} and NSS. Although the global model of substitution rate autocorrelation employed by this study is quite simple since it ignores codon positions and substitution rate correlation within local patterns of substitution (McVean 2001), it nonetheless provides a guide to the effect of autocorrelated substitution rates.

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

- Download figure
- Open in new tab
- Download powerpoint
- Download figure

The problem of false positives in NSS and Max χ^{2} is most severe for large sample sizes (*e.g.*, *m* = 50), both under constant-size populations (Figure 5b) and under population growth (Figure 5c). Although the problem is in general greater for higher substitution heterogeneity (Figure 5, top rows) it is also a problem with lower substitution rate heterogeneity (Figure 5, bottom rows).

The level of false positives of both NSS and Max χ^{2} suggests caution in interpreting evidence for recombination, especially when autocorrelated rates are an issue. For instance, inferring the presence of recombination in mitochondrial DNA should be done cautiously as substitution rate correlation is known (Yang 1995; Nielsen 1997).

The results using Φ_{w} contrast strongly with the results using the NSS (which is also compatibility based). This is likely due to the difference in the statistics themselves. The Φ_{w}-statistic uses compatibility between closely linked sites directly whereas the NSS statistic measures clustering within a compatibility matrix. As the clustering can be caused by substitution rate correlation, and not only by recombination, this might explain the difference between the two statistics. For Max χ^{2} the problem is possibly due to pairs of sequences that differ greatly on one side of a site (due to high mutation) but share a great degree of similarity on the other side of a site (due to low mutation). Local “bursts” of mutation (McVean 2001) likely exacerbate the problem, especially for linkage disequilibrium approaches that are based on allele frequencies at different sites.

#### Empirical data:

The general information concerning the empirical data sets is summarized in Table 1. Tables 5 and 6 show the results of tests for recombination on all the empirical data sets. In addition to the results obtained using the Φ_{w}-statistic, results using Max χ^{2} (Maynard Smith 1992), NSS (Jakobsen and Easteal 1996), correlation of *r*^{2} and |*D*′| with distance (Lewontin 1964; Hill and Robertson 1968), and a LPT (McVean *et al.* 2002) are shown. The estimates of ρ for the population level data sets were obtained using LDHat (McVean *et al.* 2002). Tests for recombination within populations (*i.e.*, *r*^{2}, |*D*′|, and LPT) were not applied to data sets that contained individuals from different species.

##### Recombinant examples:

Table 5 shows that the null hypothesis of no recombination is rejected by all tests for most of the suspected recombinant data sets, including the Candida example that had very little sequence diversity (0.7%). Whereas a lack of sequence diversity in the simulations made recombination harder to detect, this may be partially overcome by using longer alignments, such as that for the Candida example, which had 2553 sites. Interestingly, the null hypothesis of no recombination was not universally rejected for two of the bacterial data sets: Cowdria and *Helicobacter pylori*. For these two bacterial examples, evidence for recombination was found using the Φ_{w}-statistic as well as the coalescent-based likelihood permutation test. However, recombination was detected in the Cowdria example using the correlation of distance with *r*^{2} only after sites with minor alleles were removed. Moreover, in the *H. pylori* data set neither NSS nor Max χ^{2} found significant evidence for recombination. This could be due to the high suspected rate of recombination in the *H. pylori* example, which has conditions approaching linkage equilibrium (Suerbaum *et al.* 1998). The linkage disequilibrium methods seem to be highly sensitive to sites with low allele frequencies and consistent results are obtained only after the removal of these sites.

##### Possibly recombinant examples:

The results obtained from the data sets for which the status of recombination is debated are quite interesting (Table 6). For the Norovirus example, evidence of recombination is found using Φ_{w}, Max χ^{2}, and the LPT. There is some evidence of recombination found with *r*^{2}, but after sites with minor allele frequencies <0.1 are removed no further evidence is found by the linkage disequilibrium methods. Since the samples came from a number of different cities, it could be that evidence of recent recombination is weakened by removing these sites. However, the LPT finds evidence of recombination regardless of whether or not these sites are removed.

For the bacterial symbiont nematode Wolbachia, there is little prior reason to suspect recombination (Jiggins 2002). Nonetheless, evidence for recombination is found using correlation of *r*^{2} with distance and marginal evidence for recombination is found by using the likelihood permutation test when sites with minor alleles frequencies <0.1 are removed. The results obtained using the Φ_{w}-statistic also suggest that there is marginal evidence for recombination with Wolbachia. The possible presence of recombination in Wolbachia should be tested further using more data.

Recombination in the animal mitochondrial DNA of Apodemus was first proposed (Ladoukakis and Zouros 2001) and then disputed (Maynard Smith and Smith 2002). Tests for recombination using Φ_{w} and Max χ^{2} indicate that there is little evidence for recombination, although the NSS statistic does find evidence for recombination. The evidence for recombination within Apodemus using the Max χ^{2}-test is even weaker here than in previous studies (Maynard Smith and Smith 2002), possibly due to the fact that this implementation of the Max χ^{2}-test uses a “fixed window size.” Given the high level of false positives of NSS, the results suggest that evidence for recombination within Apodemus is lacking.

For the fungal Boletales, results using the Φ_{w}-statistic are quite distinct from the results obtained using both the NSS and the Max χ^{2}-statistic. The Φ_{w}-based tests find no evidence for recombination whereas both other tests find strong evidence for recombination. Interestingly, although most other methods for detecting recombination find evidence for recombination within this data set, Geneconv (Sawyer 1989), another powerful sequence-based test for recombination, does not (Posada 2002).

One possibility for the Boletales data set is that the Φ_{w}-statistic is too conservative and produced a type II error (“false negative”). The Boletales data set is a saturated data set with a strong A + T bias (Kretzer and Bruns 1999). The strong A + T bias results in an estimated transition/transversion ratio of 0.4. Simulations show, however, that even under such conditions, there is reason to believe that recombination will still create distinct patterns of compatibility and incompatibility that should be detectable using the Φ_{w}-statistic (results not shown). Moreover, simulations indicate that the Φ_{w}-statistic appears to be more powerful than the NSS statistic (which is also compatibility based), suggesting that a type II error for the Φ_{w}-statistic, but not for the NSS statistic, is unlikely.

Another possibility for the Boletales example is that both Max χ^{2} and the NSS statistic are producing type I errors, which, according to the simulations, autocorrelated substitution rates might induce. To test this, a parametric bootstrap with 1000 replicates simulating codons (with no recombination) was performed using a substitution rate heterogeneity of 1.31 and global substitution rate correlation ρ_{G} = 0.35 as estimated from the data set. Figure 6 shows the distribution of estimated *P*-values obtained on the 1000 replicates using the Max χ^{2}-statistic, NSS statistic, and the Φ_{w}-statistic. Recombination was inferred 5.7% of the time using the Φ_{w}-statistic, 8.5% of the time with the Max χ^{2}-statistic, and 37.5% of the time using the NSS statistic. Since none of the replicates contained recombination, the *P*-values for each of the three methods should follow a uniform distribution. Figure 6 shows that the parametric bootstrap creates conditions similar to recombination for both Max χ^{2} and NSS [a one-sided Kolmogorov–Smirnov test (Massey 1951) rejects the uniform distribution at a significance level of 10^{−7} for both Max χ^{2} and NSS but fails to find any evidence to reject the uniform distribution for Φ_{w}]. Whereas the results for Max χ^{2} are less striking than those for NSS, the parametric bootstrap fails to account for local patterns of mutation (Hey 2000; McVean 2001; McVean *et al.* 2002), which are likely to exacerbate the observed bias. These results suggest that there is reason to doubt the validity of the inferences of Max χ^{2} and NSS concerning the presence of recombination in the Boletales data set.

#### Conclusion:

We have presented a simple, powerful test for detecting recombination that can be used regardless of sample history. The approach is very general (*e.g.*, does not assume a single population) and aims to determine simply whether there is a recombinant signal present within the sequences. In contrast to two other general tests, Max χ^{2} and NSS, our test does not falsely infer the presence of recombination because of mutation rate correlation (which is present in some mitochondrial DNA). Interestingly, our approach performs very well even in the presence of population growth, in contrast to methods based on linkage disequilibrium (*r*^{2} and |*D*′|), a coalescent-based likelihood permutation test (from LDHat), Max χ^{2}, and NSS. Our method can be used by itself, or to validate the visual presence of recombination from a phylogenetic network approach, or to independently verify the presence of recombination if a positive estimate of the rate of recombination is obtained. The approach may be particularly useful in distinguishing recurrent mutation from recombination when assumptions such as a single, randomly mating, and constant-size population are not met. The test can be used easily when many sequences and sites are present because of its computational efficiency and indeed is more powerful in such circumstances. A program implementing our test as well as both Max χ^{2} and NSS is available as a stand-alone program at the following address: http://www.mcb.mcgill.ca/^{∼}trevor. The test is also implemented in SplitsTree 4.2, available at http://www.splitstree.org.

## APPENDIX A

The normal approximation to the permutation test requires calculation of the expectation and variance of the Φ_{w}-statistic under permutations of the alignment. This section contains derivations for both the mean and the variance and outlines how to compute both values efficiently. Again, assume that the proportion of informative sites is *q* and let *w* be a fixed width (in bases). Throughout this section, let *k* = *wq*.

Let *M* = (*M _{i}*

_{,j}) be a given

*n*×

*n*refined incompatibility matrix. Note that

*M*is symmetric. Let

*I*= {1, …,

*n*} be an index set. Let σ be any permutation of the index set, and define a permutation of the matrix as σ(

*M*) = (

*M*

_{σ(i),σ(j)}).

Define the sample space Ω by Ω = {σ(*M*): σ ∈ *S _{n}*}. Assume that every permutation σ is equally likely. Define an

*n*×

*n*random matrix by

*X*= σ(

*M*). Note that

*X*is symmetric, a fact that is used throughout without further mention.

Define for all 1 ≤ *i* ≤ *n*: and .

Also define .

Lemma 1. *Let X be a random matrix. Then for any arbitrary but distinct* {*i*, *j*, *k*, *l*}

*Proof*. Note that a permutation σ of *I* can be viewed as mapping to . Denote the value of σ(*i*) by σ_{i}. The total number of permutations is then *n*!. The number of permutations that have *m* distinct elements fixed in some mapping is (*n* − *m*)! (*e.g*., σ(*a*_{1}) = *b*_{1}, σ(*a*_{2}) = *b*_{2},…, σ(*a _{m}*) =

*b*). Since every permutation is equally likely the probability of such a permutation isNote that every distinct pair (

_{m}*i*,

*j*),

*i*≠

*j*can be mapped to any distinct pair (

*a*,

*b*),

*a*≠

*b*, by some σ. Note also that Pr[

*X*

_{i}_{,j}=

*M*

_{a}_{,b}] = Pr[σ

_{a}=

*i*∧ σ

_{b}=

*j*]. Finally, for notational convenience the summation is written as . Hence,▪

Consider the statistic Φ_{w} defined on a random matrix *X* asDefine (for 1 ≤ *a*, *b* ≤ *n*)Note thatThen

Theorem 1. *The expectation and variance of* Φ* _{w} can be written as*(

*for n*≥ 2

*k*),

*where*

*Moreover*,

*both E*[Φ

_{w}]

*and*Var[Φ

_{w}]

*can be calculated in O*(

*n*

^{2})

*time*.

*Proof.* The expectation is straightforward:The variance is a little more involved,whereand ≺ denotes standard lexicographical ordering.

Note that *Q _{k}* can be partitioned into two disjoint sets

*Q*

_{k}_{,0}and

*Q*

_{k}_{,1}, where

*Q*

_{k}_{,m}= {((

*a*,

*b*), (

*c*,

*d*)) ∈

*Q*: |{

_{k}*a*,

*b*} ∩ {

*c*,

*d*}| =

*m*} [by definition

*Q*does not contain pairs of the type ((

_{k}*a*,

*b*), (

*a*,

*b*))]. One way to determine

*Q*

_{k}_{,1}is to set up a recurrence.

Note thatso thatHence |*Q*_{1,1}| = (*n* − 2).

Next let ((*a*_{1}, *a*_{2}), (*a*_{3}, *a*_{4})) ∈ *Q _{k}* −

*Q*

_{k}_{−1}. Then at least one (

*a*

_{1},

*a*

_{2}) = (

*a*,

*a*+

*k*) or (

*a*

_{3},

*a*

_{4}) = (

*a*,

*a*+

*k*) must be true. Consider the four subcases:

*Case 1:* ((*a*, *b*), (*a*, *a* + *k*)), where 1 ≤ *a* ≤ *n* − *k* and *a* < *b* < *a* + *k.* There are precisely (*n* − *k*)(*k* − 1) terms of this type.

*Case 2:* ((*a*, *a* + *k*), (*b*, *a* + *k*)), where 1 ≤ *a* ≤ *n* − *k* and *a* < *b* < *a* + *k.* Again, there are precisely (*n* − *k*)(*k* − 1) terms of this type.

*Case 3:* ((*a*, *a* + *k*), (*a* + *k*, *b*)), where 1 ≤ *a* ≤ *n* − *k* and *a* + *k* < *b* ≤ min(*a* + 2*k*, *n*). For *n* ≥ 2*k* there are (*k*)((*n* − *k*) − *k*) + (*k*)(*k* − 1)/2 such terms.

*Case 4:* ((*b*, *a*), (*a*, *a* + *k*)), where 1 ≤ *a* ≤ *n* − *k* and max(1, *a* − *k*) ≤ *b* < *a.* For *n* ≥ 2*k* there are again (*k*)((*n* − *k*) − *k*) + (*k*)(*k* − 1)/2 such terms.

Cases 3 and 4 can coincide for *n* ≥ 2*k* when |*a* − *b*| = *k.* All other combinations of cases are disjoint. There are precisely (*n* − *k*) − *k* such coincidences. This gives the following recurrence for *Q _{k}*

_{,1}:The recurrence can be solved by standard techniques resulting inNote that . Since

*Q*is the disjoint union of

_{k}*Q*

_{k}_{,0}and

*Q*

_{k}_{,1}, thenThe variance of Φ

_{w}can then be written asNoting that Cov[

*X*

_{a}_{,b}

*X*

_{c}_{,d}] =

*E*[

*X*

_{a}_{,b}

*X*

_{c}_{,d}] −

*E*[

*X*

_{a}_{,b}]

*E*[

*X*

_{c}_{,d}] and Var[

*X*

_{a}_{,b}] =

*E*[

*X*

_{a}_{,b}

^{2}] −

*E*[

*X*

_{a}_{,b}]

^{2}, the constants

*c*

_{1},

*c*

_{2}, and

*c*

_{3}can be solved for using the relations from the previous lemma. Since the quantities

*u*,

*v*, and

*w*can be computed in

*O*(

*n*

^{2}) time, so can the variance and expectation. ▪

## APPENDIX B

The rate of recombination is here referred to as ρ = 4*Nrt*, where *r* is the per base recombination rate and *t* is the sequence length. Here *N* was set to 1000 (diploid population), *t* was set to 1000 as well, and *r* solved for accordingly.

For population growth ρ^{†} was obtained so that the expected number of recombinations was equal under scenarios (*i.e*., *E*_{β=5000}[*R*(*m*)] = *E*_{β=0}[*R*(*m*)]), where *R*(*m*) is the number of recombinations for a sample of size *m* (Wiuf *et al.* 2001), and β = *Nb*, where *b* is the population growth rate per generation (Wiuf *et al.* 2001). The expected number of recombinations for β = 0 can be found by the following formula (Hudson and Kaplan 1985):Table B1 shows the values used for ρ = 1 (when β = 0). For values of ρ > 1 (*e.g*., ρ = 2) one can simply double the values in the table.

Similarly, the rate of mutation is here referred to as θ = 4*N*μ*t*, where μ is the per base mutation rate and *t* is the sequence length. Under a Jukes–Cantor model if β = 0 then(Wiuf *et al.* 2001). This allows θ to be found for a fixed amount of sequence diversity *p*. For β = 5000 the appropriate value of θ was found by simulation. The values used are shown in Table B2.

## Acknowledgments

T.B. thanks Kirk and Rachel Bevan, Scott Bunnell, Daniel Huson, and Russell Steele, as well as the two anonymous referees for a number of helpful suggestions that greatly improved the manuscript. T.B. is supported by the National Science Engineering and Research Council (NSERC) (postgraduate scholarship B) and by Le Fonds québécois de la recherche sur la nature et les technologies (FQRNT grant 2003-NC-81840). D.B. is supported in part by NSERC (grant 238975-01). H.P. acknowledges Génome Québec.

## Footnotes

Communicating editor: M. Veuille

- Received July 30, 2005.
- Accepted February 3, 2006.

- Copyright © 2006 by the Genetics Society of America