Mutation-Selection

Mutation Selection (MS) model (1, 2) is the core model implemented in PalantiR. Mutation Selection is a codon substitution model that describes the interplay between mutation rate, selection strength, and population size.

For a substitution from codon \(i\) to codon \(j\): \[ s_{ij} = \psi_j - \psi_i \\ Q_{ij} = 2N_e \mu_{ij} \frac{1-e^{s_{ij}}}{1-e^{-2N_es_{ij}}} \] Only codons differing in one nucleotide positions are considered. Rates of transitions for multiple codon differences are \(0\).

\(N_e\) is the effective population size
\(\mu_{ij}\) is the mutation rate between a pair of nucleotides. This rate is taken from a nucleotide substitution model, such as Hasegawa Kishino Yano [1].
\(\psi_i\) is the absolute fitness of amino acid for codon \(i\).
\(s_{ij}\) is selection coefficient of a substitution. For synonymous substitutions, the second term of the product will be equal to \(1\).

The substitution rates are also scaled to assign particular meaning to branch length of a phylogeny. For more detail, see Rate Scaling.

Example

Setup

Load model parameters:

use_genetic_code("Standard nuclear")

# Load phylogeny
mammals <- Phylogeny("../inst/extdata/mammals.newick")

# Load nucleotide equilibrium distribution
# Nucleotide order is T,C,A,G
nuc_pi <- c(.3, .2, .25, .25)

# Load amino acid fitness vector
# Amino acid order is:
# A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y
aa_psi <- as.matrix(read.csv("../inst/extdata/amino_acid_fitness_N_1000.csv"))

Create model:

# Create nucleotide substitution model
hky <- HasegawaKishinoYano(nuc_pi)

# Create codon substitution model
ms <- MutationSelection(
    population_size = 1000,
    mutation_rate = 1e-8,
    nucleotide_model = hky,
    fitness = aa_psi[1,],
    scaling_type = "synonymous")

Here "synonymous" selects the dS gauge, so the tree lengths are interpreted as neutral nucleotide substitutions per site. Use "synonymous-per-codon" only when one branch-length unit should instead mean one expected synonymous event per codon at stationarity.

Sample sequence from model equilibrium distribution to start at root:

s <- sample_sequence(model = ms, length = 100)

Simulate:

sim <- simulate_over_phylogeny(
    phylogeny = mammals,
    model = ms,
    sequence = s)

Substitution history

Plot substitution history:

plot(sim, sites = 0:10)

Synonymous substitutions are shown in green, non-synonymous in red.

Check substitution history table:

head(sim$substitutions)

node	time	from	to	codon_from	codon_to	amino_acid_from	amino_acid_to	synonymous	color
22	0.0011622	34	36	ACC	ACG	T	T	TRUE	#16A35E
23	0.0106934	34	35	ACC	ACA	T	T	TRUE	#16A35E
24	0.0060985	34	5	ACC	TCC	T	S	FALSE	#E84B2A
24	0.0332322	5	4	TCC	TCT	S	S	TRUE	#16A35E
25	0.0661563	4	7	TCT	TCG	S	S	TRUE	#16A35E
26	0.0146581	4	0	TCT	TTT	S	F	FALSE	#E84B2A

Alignment

To export FASTA sequence:

as.fasta(sim$alignment, file = "palantir_sim.fa")

Show alignment in a scrollable viewer:

plot(sim$alignment)

References

Golding, B. and Felsenstein, J. (1990). A maximum likelihood approach to the detection of selection from a phylogeny. Journal of Molecular Evolution, 31(6):511–523.
Yang, Z. and Nielsen, R. (2008). Mutation-Selection Models of Codon Substitution and Their Use to Estimate Selective Strengths on Codon Usage. Molecular Biology and Evolution. 25(3):568-579.
Hasegawa, M., Kishino, H. and Yano, T. (1985). Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution 22:160-174.