Calculating site-specific evolutionary rates at the amino-acid or codon level yields similar rate estimates

Generation of simulated alignments

Our simulation approach was similar to the one employed by Spielman, Wan & Wilke (2016). In brief, we first generated a set of balanced, binary trees with different branch lengths and numbers of taxa, using the R package ape (Paradis, Claude & Strimmer, 2004). We then simulated sequence evolution along these trees using the python library pyvolve (Spielman & Wilke, 2015a).

We generated a total of 40 trees, using all pairwise combinations of five different branch lengths and eight different numbers of taxa. The branch lengths we used were 0.0025, 0.01, 0.04, 0.16, and 0.64. These numbers indicate the divergence in mutations per site between two nodes in a tree. The numbers of taxa we used were 16, 32, 64, 128, 256, 512, 1,024, and 2,048.

To generate alignments with site-specific dN∕dS values, we simulated sequences with 100 codon sites using a site-specific Muse–Gaut model (Muse & Gaut, 1994). To simulate sequences with constant dS, we set dS = 1 at all sites and set dN at each site to a different value randomly drawn from a uniform distribution between 0.1 and 1.6. To simulate sequences with variable dS, we assigned each site a distinct dN and dS value, by first choosing a randomly drawn dN∕dS value, then choosing a randomly drawn dS value, and then setting dN = dS × (dN∕dS). The dN∕dS values were drawn from a uniform distribution between 0.1 and 1.6, and the dS values were drawn from a uniform distribution between 0.5 and 2. We generated 50 replicate sequence alignments for each combination of branch length, number of taxa (128–2,048), and choice of dS (constant or variable), for a total of 2,500 sequence alignments.

We also generated alignments with gamma-distributed site-specific dN∕dS values. We simulated sequences with 100 codon sites using a site-specific Muse–Gaut model. We set dN∕dS to a randomly drawn value from a gamma distribution. We ran simulations for six distinct gamma distributions, using shape parameters α and rate parameters β previously estimated for six HIV-1 proteins (Table 2 in Meyer & Wilke, 2015b). For each protein (i.e., distinct gamma distribution), we generated 50 replicate sequence alignments for each combination of branch length and number of taxa (16–256), for a total of 1,500 sequence alignments per protein.

For sequences simulated according to MutSel models, we used sequence alignments previously published in Spielman, Wan & Wilke (2016), specifically the alignments simulated with unequal nucleotide frequencies. These sequences were simulated using the Halpern and Bruno model (HB98) (Halpern & Bruno, 1998), and we had alignments for the same tree parameters, dS variation (constant/variable), and replicate numbers as our simulations of the dN∕dS model, for a total of 2,500 sequence alignments.

Rate inference

To acquire the Rate4Site scores, the simulated sequences were translated into amino acids using biopython. The translated sequences were inputted into Rate4Site along with their corresponding trees. We ran Rate4Site with the following options:

 
 
rate4site -s aln_file -t tree_file -o norm_rates_file \ 
         -y orig_rates_file    

Here, aln_file is the input fasta file with aligned sequences. The file tree_file contains the phylogenetic tree. The file norm_rates_file is the output file into which Rate4Site writes z-normalized rate scores, and orig_rates_file is the output file into which Rate4Site writes original rate scores. The option -y causes Rate4Site to output original scores. (By default, Rate4Site only outputs z-transformed scores.) In our analysis we used only the original scores, renormalized such that they had a mean of 1.

We inferred site-specific dN∕dS using the one-parameter fixed-effects likelihood method (FEL1) implemented in HyPhy (Kosakovsky Pond, Frost & Muse, 2005). We ran HyPhy using the FEL1 script provided in Spielman, Wan & Wilke (2016). After running the dN∕dS inference, we explicitly set dN∕dS = 0 at all sites that did not experience any amino-acid changes. We did so because the FEL1 method assigns a site-wise dN∕dS of 1 to completely conserved sites that contain no synonymous and no non-synoymous mutations. However, dN∕dS should equal 0 at such sites in a one-parameter model, which implicitly assumes that dS is the same at all sites and hence will be non-zero even at completely conserved sites.

Analysis of empirical datasets

For analysis of empirical datasets, we used data from Spielman & Wilke (2013) and Meyer & Wilke (2015b). From Spielman & Wilke (2013), we acquired unaligned sequence data for six arbitrarily chosen membrane proteins: Mannose-6-phosphate receptor M6PR (Ensembl transcript ID: ENST00000000412), CD74 (Ensembl transcript ID: ENST00000009530), CD4 (Ensembl transcript ID: ENST00000011653), G protein-coupled receptor class C, GPRC5A (Ensembl transcript ID: ENST00000014914), Gamma-aminobutyric acid type A receptor, GABRA1 (Ensembl transcript ID: ENST00000023897), and TNF receptor superfamily member 17, TNFRSF17 (Ensembl transcript ID: ENST00000053243). We aligned the amino-acid sequences using MAFFT 7.305b (Multiple Alignment using Fast Fourier Transform) (Katoh & Standley, 2013). We ran MAFFT using default options with:

 
 
mafft input_fasta_file > output_fasta_file    

Here, input_fasta_file is the fasta file containing sequences to be aligned and output_fasta_file is the output file into which the alignment is written. The aligned amino-acid sequences were subsequently back-translated to codon sequences using the unaligned nucleotide sequences.

For sequence data from Meyer & Wilke (2015b), we used the aligned amino and nucleotide sequence files for all of the proteins used in the paper.

We inferred phylogenetic trees from amino-acid sequences using RAxML (Stamatakis, 2014). We ran RAxML with the following options:

 
 
raxmlHPC-PTHREADS-SSE3 -T 48 -s fasta_file -w output_directory \ 
        -n tree_name -m PROTCATLG -p 1    

Here, fasta_file is the input file containing the aligned sequences. RAxML outputs all output files into the directory indicated by output_directory, and tree_name is the name for the output tree files. The option -m PROTCATLG picks a protein CAT model with LG matrix for the tree inference, and the option -p 1 generates a random number seed for parsimony inference.

Finally, for all empirical datasets, we inferred Rate4Site scores and per-site dN∕dS values as described in the subsection “Rate inference.”

Relationship between Rate4Site scores and true dN∕dS

To determine the relationship between Rate4Site and dN∕dS models, we began by simulating sequence evolution with known, site-specific dN∕dS values and then comparing these true dN∕dS values to Rate4Site scores inferred from the simulated alignments (Fig. 1). We first considered the case of constant dS among all sites. Thus, for each site in each alignment, we randomly drew a dN from a uniform distribution ranging from 0.1 to 1.6. We set dS = 1 for all sites, such that the dN∕dS ratios similarly varied from 0.1 to 1.6. We ran simulations along a set of 25 balanced trees with different branch lengths and numbers of taxa, as used previously in a study of dN∕dS inference (Spielman, Wan & Wilke, 2016). Simulated sequences were 100 codon sites long, and we generated 50 replicate simulations for each simulation condition.

We calculated the correlation between each site’s true dN∕dS and its inferred Rate4Site score to assess how well the Rate4Site scores agreed with the simulated rates. We then plotted the mean correlation strengths in replicate simulations against the simulations’ branch lengths and number of taxa. We found that correlation strengths systematically increased with both increasing branch lengths and number of taxa (Fig. 2A). While correlations were low to moderate for the least-diverged and smallest alignments, for larger and/or more diverged alignments correlations approached values ranging from 0.8 to 1.0.

We also performed a comparison of the magnitude of Rate4Site scores and dN∕dS scores, by calculating root-mean-square deviations (RMSD) between these scores. Because these two types of scores are not measured in the same units, this comparison may not seem meaningful. However, we can convert both types of scores into normalized, relative scores by dividing them by their mean score. These normalized scores have comparable interpretations and RMSDs between them are meaningful quantities.

We found that RMSD values were generally moderate, between 0.1 and 0.6 (Fig. 2B). They declined with both increasing number of taxa and increasing sequence divergence. However, overall RMSD depended more strongly on branch length than on the number of taxa. Visual inspection of normalized Rate4Site scores plotted against normalized dN∕dS scores revealed no major systematic differences between these scores (Fig. 3). Differences seemed to be driven primarily by the sampling noise inherent in estimating site-specific evolutionary rates.

We repeated the same analysis but now using simulations in which dS was allowed to vary among sites as well. The dN∕dS range was kept the same as before (0.1–1.6), but now each site had its own unique dS, randomly chosen from a uniform distribution ranging from 0.5 to 2. Overall, we found similar patterns in the variable dS case as we had seen for constant dS (compare Figs. 2C, 2D to Figs. 2A, 2B. However, correlations were generally somewhat weaker (Fig. 2C) and RMSD values somewhat higher (Fig. 2D) than what we had observed for constant dS. These results were to be expected, since Rate4Site as an amino-acid based metric does not take synonymous variation into account, and thus the dS variation acts simply as added random noise on the dN∕dS scores compared to Rate4Site scores.

Relationship between Rate4Site scores and scaled selection coefficients

The dN∕dS model is not a particularly realistic model of sequence evolution, because it does not have the notion of an underlying fitness landscape. A mutation increasing fitness should fix much more rapidly than the reverse mutation decreasing fitness. However, in a dN∕dS model, both mutations fix at the same rate. To increase realism in our analysis, we next investigated Rate4Site in the context of sequences simulated with mutation–selection (MutSel) models. MutSel models are specified by scaled selection coefficients, which describe the relative fitness of different amino acids (or codons) at each site in a sequence. We can derive expected dN∕dS values from these scaled selection coefficients (Spielman & Wilke, 2015b; Dos Reis, 2015) and hence we can ask how Rate4Site scores compare to the predicted dN∕dS values in MutSel models.

For this analysis, we employed previously published sequence alignments from Spielman, Wan & Wilke (2016). These alignments had been simulated using the Halpern and Bruno model (HB98) (Halpern & Bruno, 1998) along the same 25 phylogenetic trees we employed in our previous analysis (five branch lengths in all pairwise combinations with five datasets with different numbers of taxa). As before, there were 50 replicates per simulation condition, and we again had one dataset with constant dS and one with variable dS. In the dataset with constant dS, all synonymous codons have the same fitness (neutral synonymous codons). In the dataset with variable dS, there are fitness differences among synonymous codons (non-neutral synonymous codons). See Spielman, Wan & Wilke (2016) for details of parameter choices.

Our results for sequences simulated with MutSel models were broadly similar to our results for sequences simulated with dN∕dS models (Fig. 4). As before, correlations increased and approached 1 with increasing branch lengths and numbers of taxa, and RMSDs commensurately decreased. However, correlation strengths were consistently lower and RMSD values higher for the MutSel datasets than for the dN∕dS datasets at the same sequence divergence and number of taxa (compare Figs. 4 to 2). As before, differences between normalized Rate4Site scores and normalized true dN∕dS values seemed to be driven primarily by the sampling noise inherent in estimating site-specific evolutionary rates (Fig. 5). Finally, we found only minor differences between simulations with neutral synonymous codons and simulations with non-neutral synonymous codons. However, in the latter case, correlations were generally slightly lower and RMSD values somewhat higher (compare Figs. 4C and 4D to 4A and 4B).

Relationship between Rate4Site scores and inferred dN∕dS

The preceding analyses asked to what extent Rate4Site scores reflect the known underlying parameters used to generate the sequence alignments. An alternative question, possibly more applicable to practical sequence analysis, is to what extent Rate4Site scores mirror dN∕dS values inferred on the same sequence data. To address this second question, we inferred site-wise dN∕dS values for all sites in all alignments studied in the previous two subsections. The dN∕dS values were inferred using the one-rate fixed-effects likelihood method (FEL1) implemented in HyPhy (Kosakovsky Pond, Frost & Muse, 2005). The FEL1 method assigns one dN value per site and one dS value across all sites in the sequence (Spielman, Wan & Wilke, 2016). Therefore, the variation in the inferred dN∕dS values is captured entirely by dN.

We found that Rate4Site scores were very highly correlated to inferred dN∕dS across all datasets and simulation conditions (Figs. 6 and 7). For sequences simulated with the dN∕dS model, correlations for all branch lengths exceeded 0.8 (Figs. 6A and 6C). Correlations generally increased and approached 1 for sequence alignments with more divergent sequences or more taxa. RMSD values were large for the smallest and least-diverged alignments but declined rapidly as either branch length or number of taxa increased (Figs. 6B and 6D). There was little difference between alignments simulated with constant dS and with variable dS (Figs. 6A vs. 6C and 6B vs. 6D). Results for sequences simulated with MutSel models were similar (Fig. 7). The main difference was that correlation coefficients were more sensitive to the number of taxa. For the lowest number of taxa (128) correlation coefficients were systematically lower when sequences were simulated with MutSel models rather than with dN∕dS models. The pattern reversed for the highest number of taxa (2,048). RMSDs, however, were systematically higher for sequences simulated with MutSel models rather than with dN∕dS models. In all cases, there was little difference between sequences simulated with neutral synonymous codons and with non-neutral synonymous codons (Figs. 7A vs. 7D and 7B vs. 7D).

Analysis of more realistic parameter settings and of empirical datasets

The above simulations utilized rates that were uniformly distributed across sites. In empirical datasets, however, the majority of sites evolves slowly and only very few sites evolve rapidly, with a resulting distribution that is approximately gamma (see e.g., Meyer & Wilke, 2015b). Therefore, we wanted to assess if the observed relationships between Rate4Site scores and dN∕dS hold for simulations with gamma-distributed true rates. We simulated sequences with gamma-distributed true dN∕dS values, using the shape and rate parameters obtained by Meyer & Wilke (2015b) for six different HIV-1 proteins. For these simulations, we also included a range of smaller numbers of taxa (16, 32, and 64) than before, to increase realism in the simulations.

We found a wide range of correlations between Rate4Site and true dN∕dS (Fig. 8A and Parts A of Figs. S1–S5), but the overall pattern was similar to what we had previously observed for uniformly distributed rates: Correlations increased and approached 1 as either sequence divergence or the number of taxa increased. Also note that the lowest number of taxa now was 16 vs. 128 in Fig. 2, and thus at identical numbers of taxa and branch lengths, the correlations for gamma-distributed rates were generally higher than the correlations for uniformly distributed rates. On the other hand, the RMSD values between Rate4Site scores and true dN∕dS were higher than in our previous analysis (compare Figs. 8B to 2B).

We also inferred dN∕dS from the simulated sequences and compared them to the Rate4Site scores. Unlike in our previous analysis on inferred dN∕dS (Fig. 6), the correlations between Rate4Site and inferred dN∕dS now spanned a wide range of values (Figs. 8C and Parts C of Figs. S1–S5) but were generally higher than the corresponding correlations between Rate4Site and true dN∕dS. On the flip side, RMSD values were larger when calculated for inferred dN∕dS than when calculated for true dN∕dS, and there was a notable uptick in RMSD at the highest branch lengths (Figs. 8D and Parts D of Figs. S1–S5). These observations can be explained with imprecise point estimates produced by both Rate4Site and dN∕dS inference at both low and high sequence divergence: In general, both methods correctly assess which sites in an alignment evolve more rapidly than which other sites, and hence correlations are high, but the specific rate estimates assigned to individual sites are poor at both low and high sequence divergence, and hence RMSDs are high as well.

Finally, we asked to what extent the results found for simulated sequences carry over to empirical datasets. We inferred both Rate4Site scores and site-wise dN∕dS in two distinct datasets, one consisting of six membrane proteins in mammals (taken from Spielman & Wilke, 2013) and one consisting of six alignments of HIV-1 genes (taken from Meyer & Wilke, 2015b). For both datasets, we measured the correlation between the original (non-normalized) dN∕dS and normalized Rate4Sites scores. We found that Rate4Site scores and dN∕dS were highly correlated, with correlation coefficients exceeding 0.8 in all cases (Fig. 9). Thus, we can conclude that our simulation results carry over to empirical datasets, and that Rate4Site scores are generally comparable to inferred dN∕dS values.

Note that we plotted raw (non-normalized) dN∕dS in Fig. 9, and as a result points are not expected to fall onto the x = y line. We plotted the data in this way because dN∕dS values are not typically normalized to their mean, and we wanted to assess how similar in magnitude raw dN∕dS scores are to Rate4Site scores. We found that for proteins under strong positive selection (such as HIV-1 gp120), raw dN∕dS is virtually identical to Rate4Site. However, for more typical proteins that experience little positive selection, raw dN∕dS values tend to be up to a factor 10 smaller than the corresponding Rate4Site scores.