Selection in Overlapping Genes Figure 3. Simulation results in same-strand (SS) and opposite-strand (OP) overlaps. Estimations of the ratios of nonsynonymous to synonymous rates in the two genes (v1 and v2) by simultaneous estimation (blue line) and by independent estimation (red line) are plotted against the true value (X axis, dashed green line) for five types of overlap. The simulated value of v1 was set to 0.2 and v2 was varied between 0.2 and 2. k was set to 1 and t was set to 0.35. Each data point is the median of 100 replications. Vertical lines mark the lower and upper quartiles. Top: estimation of v1. Bottom: estimation of v2. Dotted black lines (X = 1 and Y = 1) illustrate the range of parameters that result in false inference of positive selection by independent estimation, i.e., when simulated v2v1 and estimated v2w1. doi:10.1371/journal.pone.0003996.g003 ‘ðtÞ~ X X nAij log pAipAij ðtÞ i j z X X nBij log pBipBij ðtÞ i j 9 ð Þ Here, pAi and pBi are the equilibrium frequency of codons in gene A and gene B respectively, estimated from the non- A B overlapping regions of the genes. n ij and n ij are the number of sites in the alignment consist of codons i and j for gene A and gene B, respectively. The method was implemented in Matlab and is available at http://nsmn1.uh.edu/,dgraur/Software.html. Running time is ,7 seconds for a pair of aligned sequences of length 1000 codons. Similar to the single-gene model, this method can be extended to deal with multiple sequences in a phylogenetic context and to test hypotheses concerning variable selection pressures among lineages and sites [40–42]. Results Simulation studies We tested the performance of our new method for simultaneous estimation of selection intensities in comparison to the indepen- dent estimation that does not account for gene overlap (as described in equation 1). We examined the effects of nonsynon- ymous/synonymous rate ratio in each gene (v1 and v2), transition/transversion rate ratio (k), and sequence divergence (t). In all of the methods, we used the F364 model [38] to estimate codon equilibrium frequencies. For each set of parameters, we generated 100 replications of random overlapping gene pairs (each gene was 2000 codons in length with 1000 codons in the overlap) by sampling codons from a uniform distribution of sense codons. To simulate the evolution along a branch of length t, we divided the sequence of the overlapping gene pair into three regions: non- overlapping region of gene one, non-overlapping region of gene two, and overlapping region. For the non-overlapping regions, we calculated the transition-probability matrixes based on the non- overlapping model in equation 1. For the overlapping region, we calculated the transition-probability matrixes (based on the overlapping models in equations 5 and 6). Using the three probability matrixes, we simulated nucleotide substitutions at each codon independently [38]. Different selection pressures To examine the effect of different selection pressures, we initially set k = 1 and t = 0.35, which result in a sequence divergence of ,10%. We set v1~0:2 and varied v2 between 0.2 and 2. In Figure 3, we compare the simultaneous estimation of v1 and v2 (blue line) and the independent estimation (red line) to the true simulated value (X axis, dashed green line) in the five types of overlaps. Each data point is the median of 100 replications. We use the median rather than mean since ratios are not normally distributed. In all overlap types, the estimation of our method is in near-perfect match to the simulated value (blue and green lines, Figure 3) and the bias in the independent estimation of v2 is greater than that of v1. As expected, we found a similar pattern of bias in all overlap types except opposite-strand phase 2. In all of these overlap types (same-strand phase 1, same-strand phase 2, opposite-strand phase 0, and opposite-strand phase 1), the independent estimation of v1 is overestimated for v2v1 and underestimated for v2w1. The independent estimation of v2 is overestimated throughout the range of the simulation resulting in the false inference of positive selection in gene 2, while in reality this gene is under weak purifying selection. For example, the independent estimation of v2 in same-strand phase 1 is greater than one (apparent positive selection) for simulated values of v2 between 0.5 and one. The bias in opposite-strand phase 2 differs from the other overlap types because this overlap contains positions that are synonymous in both genes (Figure 1). Because of this factor, the independent estimation of v1 is underestimated for v2v1 and overestimated for v2w1. The independent estimation of v2 is underestimated throughout the range of the simulation, resulting PLoS ONE | www.plosone.org 4 December 2008 | Volume 3 | Issue 12 | e3996