Branch Heterogeneity

Branch heterogeneous simulations specify different models (called modes) for each branch on a phylogeny. This is done by providing a newick tree with model indexes instead of branch length, and a list of models.

For more details on time-heterogeneous simulation, see Temporal Heterogeneity Rescaling.

Here is an example file, which represents a model switch on the primate lineage in the mammalian phylogeny:

scan("../inst/extdata/mammals_switch.newick", what = character(), quiet = T)

## [1] "((((((((((((((Human:1,Pan_troglodytes:1):1,Gorilla_gorilla:1):1,Pongo_abelii:1):1,Nomascus_leucogenys:1):1,(Macaca_mulatta:1,Papio_hamadryas:1):1):1,Callithrix_jacchus:1):1,Tarsius_syrichta:1):1,(Microcebus_murinus:1,Otolemur_garnettii:1):1):1,Tupaia_belangeri:1):1,(((((Mus_musculus:0,Rattus_norvegicus:0):0,Dipodomys_ordii:0):0,Cavia_porcellus:0):0,Ictidomys_tridecemlineatus:0):0,(Oryctolagus_cuniculus:0,Ochotona_princeps:0):0):0):0,((((Vicugna_pacos:0,(Tursiops_truncatus:0,(Bos_taurus:0,Ovis_aries:0):0):0):0,Sus_scrofa:0.):0,((Equus_caballus:0,(Felis_catus:0,((Ailuropoda_melanoleuca:0,Mustela_putorius_furo:0):0,Canis_familiaris:0):0):0):0,(Myotis_lucifugus:0,Pteropus_vampyrus:0):0):0):0,(Erinaceus_europaeus:0,Sorex_araneus:0):0):0):0,(((Loxodonta_africana:0,Procavia_capensis:0):0,Echinops_telfairi:0):0,(Dasypus_novemcinctus:0,Choloepus_hoffmanni:0):0):0):0,(Monodelphis_domestica:0,(Macropus_eugenii:0,Sarcophilus_harrisii:0):0):0):0,Ornithorhynchus_anatinus:0);"

This phylogeny can be loaded with the Phylogeny constructor, with an additional type = "mode" argument.

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

# Load branch modes
mammal_switches <- Phylogeny("../inst/extdata/mammals_switch.newick", type = "mode")

We can view the phylogeny and the assigned modes:

plot(mammals, mammal_switches)

Example

In the example below, we simulate on a mammalian phylogeny with a population size expansion on the primate lineage.

Setup

Load model parameters:

use_genetic_code("Standard nuclear")

# 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, ], 5000)

Create models:

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

# Create codon models

# Smaller population size
ms_5k <- MutationSelection(
    population_size = 5000,
    mutation_rate = 1e-8,
    nucleotide_model = hky,
    fitness = psi)

# Larger population size
ms_30k <- MutationSelection(
    population_size = 30000,
    mutation_rate = 1e-8,
    nucleotide_model = hky,
    fitness = psi)

Sample sequence from model equilibrium distribution to start at root:

s <- sample_sequence(ms_5k, 100)

Simulate:

sim <- simulate_over_interval_phylogeny(
    phylogeny = mammals,
    mode_phylogeny = mammal_switches,
    models = list(ms_5k, ms_30k),
    sequence = s,
    start_mode = 0)

Note that we have to specify which mode the root starts with using start_mode argument.

Substitution history

Plot substitution history:

plot(sim, sites = 0:24, circle_opacity = 0.6)

Synonymous substitutions are shown in green, non-synonymous in red, and branch modes as colored branches. An increase in population size is expected to produce a transient excess of non-synonymous substitutions, as selection becomes strong enough to replace mildly deleterious variants that fixed under the smaller population size. The shifted clade is short, so all of its branches fall within this transient period.

The two regimes differ more in the direction of amino-acid change than in its rate. At stationarity, substitutions that increase fitness and substitutions that decrease it occur at equal rates, so the mean fitness change per substitution is zero. After the expansion, fitness-increasing substitutions predominate until the clade reaches its new equilibrium. The intervals table assigns each branch to its model, which allows the shifted clade to be compared with the remainder of the tree:

aa <- c("A", "C", "D", "E", "F", "G", "H", "I", "K", "L",
        "M", "N", "P", "Q", "R", "S", "T", "V", "W", "Y")
nsyn <- subset(sim$substitutions, !synonymous)
nsyn$group <- ifelse(sim$intervals$state[match(nsyn$node, sim$intervals$node)] == 1,
                     "shifted clade", "rest of tree")
d_psi <- psi[match(nsyn$amino_acid_to, aa)] - psi[match(nsyn$amino_acid_from, aa)]
s_scaled <- 2 * 5000 * d_psi   # scaled by the ancestral population size

counts <- rbind(gain = tapply(d_psi > 0, nsyn$group, sum),
                loss = tapply(d_psi < 0, nsyn$group, sum))
mean_s <- tapply(s_scaled, nsyn$group, mean)
se_s   <- tapply(s_scaled, nsyn$group, function(x) sd(x) / sqrt(length(x)))

par(mfrow = c(1, 2), mar = c(4, 4.5, 2, 1))
bars <- barplot(counts, beside = TRUE, col = c("#2E8B57", "#C0392B"),
                ylab = "Non-synonymous substitutions", ylim = c(0, max(counts) * 1.15),
                legend.text = c("fitness gain", "fitness loss"),
                args.legend = list(bty = "n", x = "topright"))
text(bars, counts, labels = counts, pos = 3, cex = 0.9)
bars2 <- barplot(mean_s, col = "#8DA0CB", ylab = "Mean scaled fitness change (2Ns)",
                 ylim = range(0, mean_s + sign(mean_s) * (se_s + 0.3)))
arrows(bars2, mean_s - se_s, bars2, mean_s + se_s, angle = 90, code = 3, length = 0.05)
abline(h = 0)

In the rest of the tree, gains and losses occur in nearly equal numbers and the mean population-scaled fitness change is close to zero, as expected at stationarity. In the shifted clade, gains outnumber losses 104 to 5 and the mean scaled change per substitution is +1.4 (error bars are standard errors of the mean).

Check substitution history table:

head(sim$substitutions)

node	time	from	to	codon_from	codon_to	amino_acid_from	amino_acid_to	synonymous	color
16	0.0075900	57	60	GGT	GGG	G	G	TRUE	#16A35E
59	0.0016043	57	59	GGT	GGA	G	G	TRUE	#16A35E
63	0.0112522	57	58	GGT	GGC	G	G	TRUE	#16A35E
63	0.0621408	58	57	GGC	GGT	G	G	TRUE	#16A35E
63	0.0716124	57	41	GGT	AGT	G	S	FALSE	#E84B2A
66	0.0309177	57	58	GGT	GGC	G	G	TRUE	#16A35E

Alignment

To export FASTA sequence:

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

Show alignment in a scrollable viewer:

plot(sim$alignment)