GenomegaMap: Within-Species Genome-Wide dN/dS Estimation from over 10,000 Genomes

Abstract

The dN/dS ratio provides evidence of adaptation or functional constraint in protein-coding genes by quantifying the relative excess or deficit of amino acid-replacing versus silent nucleotide variation. Inexpensive sequencing promises a better understanding of parameters, such as dN/dS, but analyzing very large data sets poses a major statistical challenge. Here, I introduce genomegaMap for estimating within-species genome-wide variation in dN/dS, and I apply it to 3,979 genes across 10,209 tuberculosis genomes to characterize the selection pressures shaping this global pathogen. GenomegaMap is a phylogeny-free method that addresses two major problems with existing approaches: 1) It is fast no matter how large the sample size and 2) it is robust to recombination, which causes phylogenetic methods to report artefactual signals of adaptation. GenomegaMap uses population genetics theory to approximate the distribution of allele frequencies under general, parent-dependent mutation models. Coalescent simulations show that substitution parameters are well estimated even when genomegaMap's simplifying assumption of independence among sites is violated. I demonstrate the ability of genomegaMap to detect genuine signatures of selection at antimicrobial resistance-conferring substitutions in Mycobacterium tuberculosis and describe a novel signature of selection in the cold-shock DEAD-box protein A gene deaD/csdA. The genomegaMap approach helps accelerate the exploitation of big data for gaining new insights into evolution within species.

Keywords: adaptation; big data; dN/dS; natural selection; parent-dependent mutation; recombination.

Fig. 1.

Comparison of omegaMap and genomegaMap estimates of the $d_{\mathrm{N}}/d_{\mathrm{S}}$ ratio ω along the porB3 outer membrane protein gene of Neisseria meningitidis. Solid lines and shaded regions show the point estimates (posterior medians) and 95% credibility intervals, respectively, for omegaMap (in blue) and genomegaMap (in red). The genomegaMap runs were 4.9 times faster for these 23 sequences at 92 min each.

Fig. 2.

Performance of genomegaMap inference of ω, κ, and θ in simulations. In the Unlinked simulations (top row), every codon was simulated independently, favoring the genomegaMap assumption. In the Clonal simulations (bottom row), all codons were completely linked, disfavoring the genomegaMap assumption. Point estimates (posterior medians) and 95% credibility intervals are indicated by the circles and solid vertical lines, respectively, the latter colored red when they exclude the actual parameter. The number of simulations (out of 100) in which the 95% credibility intervals included the actual values of ω, κ, and θ were 98, 98, and 97 in the Unlinked simulations and 92, 92, and 88 in the Clonal simulations. The correlation between the point estimates and actual values of $\text{log}\hspace{0.17em}\omega ,\hspace{0.17em}\hspace{0.17em}\text{log}\hspace{0.17em}\kappa$, and $\text{log}\hspace{0.17em}\theta$ were 0.86, 0.69, and 0.92 in the Unlinked simulations and 0.82, 0.61, and 0.88 in the Clonal simulations.

Fig. 3.

The evidence for positive selection across 3,979 genes in 10,209 Mycobacterium tuberculosis genomes. Each column is a stacked bar chart showing the proportion of codons in one gene with a given strength of evidence for positive selection, indicated by color. Blue indicates weakest evidence, $\mathrm{Pr}(\omega >1)\approx 0$, whereas red indicates strongest evidence, $\mathrm{Pr}(\omega >1)\approx 1$. Genes are ordered left-to-right by the mean $\mathrm{Pr}(\omega >1)$ across codons, from highest to lowest. Notable genes containing codons with strong evidence of positive selection are labeled; these occur across the spectrum. The genes with predominantly sky blue color, scattered between pncA and katG, contained little information because they mapped poorly to the reference genome.

Fig. 4.

Evidence of positive selection in ten Mycobacterium tuberculosis genes across 10,209 genomes. Genes are ordered by the mean $\mathrm{Pr}(\omega >1)$ across codons, from highest (gidB) to lowest (gyrA). Point estimates (black points) and 95% credibility intervals (gray bars) for ω are shown across codons. Codons for which $\mathrm{Pr}(\omega >1)\ge 0.9$ are highlighted with yellow boxes. Stacked points indicate the number of alleles that are nonsynonymous (orange) or synonymous (green) relative to the commonest allele.