Hasegawa-Kishino-Yano

Hasegawa Kishino Yano (HKY) [1] is a nucleotide substitution model which accounts for nucleotide equilibrium distribution rates of transitions (\(\alpha\)) and transversions(\(\beta\)).

\[ Q_{ij} = \begin{bmatrix} . & \alpha \pi_C & \beta \pi_A & \beta \pi_G \\ \alpha \pi_T & . & \beta \pi_A & \beta \pi_G \\ \beta \pi_T & \beta \pi_C & . & \alpha \pi_G \\ \beta \pi_T & \beta \pi_C & \alpha \pi_A & . \\ \end{bmatrix} \] In PalantiR The HKY model is most commonly used to define nucleotide transition rates for Mutation Selection model, but it can also be used as a stand-alone simulation model.

Example

Load model parameters:

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

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

Create model:

hky <- HasegawaKishinoYano(
    equilibrium = pi,
    transition_rate = 2,
    transversion_rate = 1)

We first need to provide a sequence that to start off the simulation at the root. We can sample a sequence from the equilibrium distribution of the model we are simulating:

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

Simulate:

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

Substitution history

Plot substitution history:

plot(sim, sites = 0:10)

Transition substitutions are shown in green, transversions in red.

Check substitution history table:

head(sim$substitutions)

node	time	from	to	nucleotide_from	nucleotide_to	transition	color
22	0.0026142	1	2	C	A	FALSE	#E84B2A
23	0.0240532	1	0	C	T	TRUE	#16A35E
24	0.0137177	1	0	C	T	TRUE	#16A35E
27	0.1724469	1	0	C	T	TRUE	#16A35E
31	0.0057174	1	2	C	A	FALSE	#E84B2A
33	0.0566507	2	0	A	T	FALSE	#E84B2A

We can change the colors on the plot by changing columns in the substitutions dataframe:

transition_substitutions <- sim$substitutions$transition == TRUE
sim$substitutions$color[transition_substitutions] <- "orange"
sim$substitutions$color[!transition_substitutions] <- "blue"
plot(sim, sites = 0:10)

Alignment

To export FASTA sequence:

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

Show alignment in a scrollable viewer:

plot(sim$alignment)

References

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.