BUSTED[S]: Detecting Positive Selection in the Presence of SRV

We adapted the existing BUSTED test of positive selection (Murrell et al. 2015 (link)) to account for the presence of SRV and call the new method BUSTED[S]. To explore the generality of our findings about FPRs in the presence of SRV we also investigated a second existing test of selection, the M1a versus M2a comparison from Wong et al. (2004) (link), modified slightly to employ MG94 substitution models.
BUSTED[S] is a straightforward extension of BUSTED (Murrell et al. 2015 (link)). The nucleotide substitution process is modeled using the standard finite state continuous time Markov process approach of Muse and Gaut (1994) (link), with entries of the instantaneous rate matrix Q corresponding to substitutions between sense codons i and j denoted as
$q_{ij}=\{\begin{array}{cc}{\alpha }^s{\theta }_{ij}{\pi }_j^p\hfill & 1-\text{step} \text{synonymous} \text{change},\hfill \\ {\alpha }^s{\omega }^{bs}{\theta }_{ij}{\pi }_j^p\hfill & 1-\text{step} \text{nonsynonymous} \text{change},\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$
The θ_ij (with ${\theta }_{ij}={\theta }_{ji}$ ) are parameters governing nucleotide substitution biases. For example, ${\theta }_{\text{ACT},\text{AGT}}={\theta }_{\text{CG}}$ and because we incorporate the standard nucleotide GTR model there are five identifiable θ_ij parameters: ${\theta }_{\text{AC}},\hspace{0.17em}{\theta }_{\text{AT}},\hspace{0.17em}{\theta }_{\text{CG}},\hspace{0.17em}{\theta }_{\text{CT}}$ , and ${\theta }_{\text{GT}}$ , with ${\theta }_{\text{AG}}\equiv 1$ . The position-specific equilibrium frequency of the target nucleotide of a substitution is ${\pi }_j^p$ ; for example, it is ${\pi }_{\mathrm{G}}^2$ for the second-position change associated with $q_{\text{ACT},\text{AGT}}$ . The ${\pi }_j^p$ and the stationary frequencies of codons under this model are estimated using the CF3 × 4 procedure (Kosakovsky Pond et al. 2010 (link)), adding nine parameters to the model. The ratio of nonsynonymous to synonymous substitution rates for site s along branch b is ω^bs, and this ratio is modeled using a 3-bin general discrete distribution (GDD) with five estimated hyperparameters: $0\le {\omega }_1\le {\omega }_2\le 1\le {\omega }_3,\hspace{0.17em}{p}_1=P({\omega }^{bs}={\omega }_1)$ , and $p_2=P({\omega }^{bs}={\omega }_2)$ . The procedure for efficient computation of the phylogenetic likelihood function for these models was described in Kosakovsky Pond et al. (2011) (link). The quantity α^s is a site-specific synonymous substitution rate (no branch-to-branch variation is modeled) drawn from a separate 3-bin GDD. The mean of this distribution is constrained equal to one to maintain statistical identifiability, resulting in four estimated hyperparameters: $0\le c{\alpha }_1<{\alpha }_2=c\le c{\alpha }_3,\hspace{0.17em}{f}_1=P({\alpha }^s={\alpha }_1)$ , and $f_2=P({\alpha }^s={\alpha }_2)$ , with c chosen to ensure that $E\left\{{\alpha }^s\right\}=1$ . Typical implementations, including ours, allow the number of α and ω rate categories to be separately adjusted by the user, for example, to minimize AIC_c or to optimize some other measure of model fit. The default setting of three categories generally provides a good balance between fit and performance when using this GDD approach for modeling. Our HyPhy implementation of BUSTED[S] will warn the user if there is evidence of model overfitting, such as the appearance of rate categories with very similar estimated rate values or very low frequencies.
The BUSTED[S] procedure for identifying positive selection is the likelihood ratio test comparing the full model described above to the constrained model formed when ω₃ is set equal to 1 (i.e., no positively selected sites). Critical values of the test are derived from a $50:50$ mixture distribution of ${\chi }_0^{\mathrm{2}}$ and ${\chi }_2^{\mathrm{2}}$ . Note that this asymptotic statistic differs from the 3-component mixture used by Murrell et al. (2015) (link); the simulation studies performed in the current study suggest that this less conservative mixture is sufficient to maintain nominal Type I errors. Both BUSTED[S] and BUSTED analyses in the current work use the same $50:50$ mixture test statistic. BUSTED[S] reduces to BUSTED by setting ${\alpha }^s=1$ , that is, by placing all the mass of the synonymous rate heterogeneity distribution at α = 1. The method is implemented as a part of HyPhy (version 2.5.1 or later). BUSTED[S] is available for free public use on the Datamonkey webserver (Weaver et al. 2018 (link)) at https://www.datamonkey.org/BUSTED (last accessed February 24, 2020).