Shared Substitution History

Shared substitution history is a switching simulation type similar to Markov-Modulated Mutation-Selection. However, instead of simulating each site with independent switching in the parent model, the switches are shared across sites. In essence, this is a conditioned simulation, given identities and times of parent model switches.

Unlike Shared Model Time History, this simulation shares times and identities.

This type of simulation is useful to emulate population size switches, such that population size change is common for all sites in the alignment.

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 equilibrium frequencies
# These will be converted to population-specific fitness vectors
# (the raw file carries a spurious 21st column from a trailing tab,
# so keep the twenty amino-acid columns)
aa_pi <- as.matrix(read.table("../inst/extdata/psi-c50-1.txt"))[, 1:20]

psi <- equilibrium_to_fitness(aa_pi[1, ], 1000)

Create models:

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

# Parent switching model
gtr <- GeneralTimeReversible(
    equilibrium = rep(1/2, 2),
    exchangeability = matrix(rep(1, 4), nrow = 2))

# Create codon models

# Small population size
ms_1k <- MutationSelection(
    population_size = 1000, 
    mutation_rate = 1e-8, 
    nucleotide_model = hky,
    fitness = psi,
    scaling_type = "synonymous")

# Large population size
ms_20k <- MutationSelection(
    population_size = 20000,
    mutation_rate = 1e-8, 
    nucleotide_model = hky,
    fitness = psi,
    scaling_type = "synonymous")

Both substitution models use the dS gauge, so branch lengths are neutral nucleotide substitutions per site even as the shared process switches between population sizes.

Sample sequence from model equilibrium distribution to start at root:

s <- sample_sequence(ms_1k, 100)

Simulate:

sim <- simulate_with_shared_substitution_heterogeneity(
    phylogeny = mammals,
    switching_model = gtr,
    switching_rate = 10,
    substitution_models = list(ms_1k, ms_20k),
    sequence = s,
    start_mode = 0)

Note that we have to specify which mode the root starts with using start_mode argument. We can additionally accelerate switching model substitution rate with switching_rate.

Substitution history

Plot substitution history:

plot(sim, circle_opacity = 0.1)

Synonymous substitutions are shown in green, non-synonymous in red. Branch modes are shown as colored branches. Each mode has its own color. After every model switch, there is a burst of non-synonymous substitutions.

Check substitution history table:

head(sim$substitutions)

site	node	time	from	to	codon_from	codon_to	amino_acid_from	amino_acid_to	synonymous	color
0	2	0.0373791	57	60	GGT	GGG	G	G	TRUE	#16A35E
0	24	0.0253312	60	59	GGG	GGA	G	G	TRUE	#16A35E
0	65	0.0300350	60	57	GGG	GGT	G	G	TRUE	#16A35E
0	80	0.0593585	57	58	GGT	GGC	G	G	TRUE	#16A35E
0	82	0.0145228	57	60	GGT	GGG	G	G	TRUE	#16A35E
0	82	0.1180168	60	28	GGG	CGG	G	R	FALSE	#E84B2A

Check model switching history table. Each row corresponds to a tree node (branch). state is the index of the model, from and to are branch time points where the model applies:

head(sim$intervals)

node	state	from	to
1	0	0.0000000	0.0561826
1	1	0.0561826	0.0716640
2	1	0.0000000	0.0368998
2	0	0.0368998	0.0948974
2	1	0.0948974	0.1493364
2	0	0.1493364	0.1784608

Alignment

To export FASTA sequence:

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

Show alignment in a scrollable viewer:

plot(sim$alignment)

de Koning Lab, University of Calgary · 2016–2026 · lab.jasondk.io · GitHub