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)

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

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	25	0.0001169	57	58	GGT	GGC	G	G	TRUE	#16A35E
0	32	0.0556395	57	53	GGT	GAT	G	D	FALSE	#E84B2A
0	34	0.0399654	53	57	GAT	GGT	D	G	FALSE	#E84B2A
0	35	0.0092495	57	53	GGT	GAT	G	D	FALSE	#E84B2A
0	35	0.1793139	53	57	GAT	GGT	D	G	FALSE	#E84B2A
0	47	0.0263713	57	49	GGT	GCT	G	A	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