Model

The Goldman–Yang 1994 (GY94) model describes codon substitutions using a transition/transversion rate ratio, kappa, and a non-synonymous/synonymous rate ratio, omega. For distinct sense codons (i) and (j), PalantiR first constructs the unscaled generator

\[ \widetilde q_{ij} = \begin{cases} 0, & \text{if }i\text{ and }j\text{ differ at more than one nucleotide},\\ \pi_j, & \text{synonymous transversion},\\ \kappa\pi_j, & \text{synonymous transition},\\ \omega\pi_j, & \text{non-synonymous transversion},\\ \omega\kappa\pi_j, & \text{non-synonymous transition}. \end{cases} \]

Here (pi_j) is the stationary frequency of target codon (j). Diagonal entries are the negative row sums. Both omega and kappa must be finite and strictly positive. The familiar interpretation is omega < 1 for purifying selection, omega = 1 for neutral amino-acid evolution, and omega > 1 for an excess of non-synonymous change, although the constructor does not impose those biological categories.

use_genetic_code("Standard nuclear")

gy <- GoldmanYang94(
    omega = 0.2,
    kappa = 2,
    frequencies = c(T = 0.30, C = 0.20, A = 0.25, G = 0.25),
    frequency_model = "F1x4",
    scaling_type = "standard")

gy$omega

## [1] 0.2

gy$scaling_type

## [1] "substitution"

GY94() is a short alias for GoldmanYang94().

Stationary codon frequencies

All frequency adapters operate on the sense codons of the active genetic code. They remove stop codons and normalize the resulting sense-codon probabilities to sum to one.

FEqual() assigns equal probability to every sense codon.
F1x4() accepts four nucleotide frequencies shared by all three codon positions.
F3x4() accepts a 3 by 4 matrix (positions by nucleotides), its 4 by 3 transpose, or a length-12 vector in three position blocks.
F61() accepts one positive value per sense codon. The name is conventional: the actual length follows the active genetic code. Named input is safest and must contain every sense codon exactly once.

The nucleotide order for unnamed input is T, C, A, G. Named input may be in any order and is rearranged internally.

use_genetic_code("Standard nuclear")

pi_equal <- FEqual()
pi_1x4 <- F1x4(c(T = 0.30, C = 0.20, A = 0.25, G = 0.25))

position_frequencies <- rbind(
    c(T = 0.30, C = 0.20, A = 0.25, G = 0.25),
    c(T = 0.20, C = 0.30, A = 0.30, G = 0.20),
    c(T = 0.25, C = 0.25, A = 0.20, G = 0.30))
pi_3x4 <- F3x4(position_frequencies)

f61_input <- setNames(rep(1, length(GeneticCode())), names(GeneticCode()))
pi_61 <- F61(f61_input)

sum(pi_3x4)

## [1] 1

head(pi_3x4)

##        TTT        TTC        TTA        TTG        TCT        TCC 
## 0.01590668 0.01590668 0.01272534 0.01908802 0.02386002 0.02386002

codon_frequencies() performs the same adaptation and can infer F1x4, F3x4, or F61 from input length. An explicit frequency_model is recommended in reproducible scripts. A CodonFrequencies object records the active genetic code; PalantiR rejects it if the genetic code is subsequently changed, because the codon state ordering may no longer agree.

Branch-length scaling

Let $\mathcal A$ be the set of ordered codon transitions that determines scaling. For one GY94 process, the corresponding stationary event rate is

\[ S_{\mathcal A} = \sum_i \pi_i \sum_{j:(i,j)\in\mathcal A}\widetilde q_{ij}. \]

The available choices for $\mathcal A$ are:

"standard" or "substitution": all codon substitutions;
"non-synonymous": amino-acid-changing substitutions only;
"synonymous": synonymous substitutions only;
"none": no normalization.

"standard" is an input alias for "substitution"; the model stores the canonical value "substitution". A standalone model uses $Q=\widetilde Q/S_{\mathcal A}$, so one branch-length unit represents one expected event from $\mathcal A$ per codon at stationarity.

For a discrete mixture with class proportions $p_k$, PalantiR follows the mixture-wide convention used for codon site models in PAML:

\[ \overline S_{\mathcal A}=\sum_k p_k S_{\mathcal A,k}, \qquad Q_k=\widetilde Q_k/\overline S_{\mathcal A}. \]

Every component therefore shares one denominator, while its stationary rate $S_{\mathcal A,k}/\overline S_{\mathcal A}$ is retained. The mixture-weighted mean of those rates is one; individual component rates need not be one. Under synonymous or non-synonymous scaling, independent normalization would erase omega-dependent rate differences among classes.

For a branch-site model, PalantiR computes the weighted denominator separately for each branch type, because the collection of active class-specific generators differs between background and foreground branches:

\[ \overline S_{\mathcal A,b}=\sum_k p_k S_{\mathcal A,k,b}. \]

All site classes active on branch type $b$ share $\overline S_{\mathcal A,b}$. Thus, under standard scaling, branch length is the expected number of codon substitutions per site averaged over the whole mixture. Synonymous and non-synonymous scaling give the analogous mean for the specified transition set. This is the convention implemented for site and branch-site models in the PAML codeml source, which notes that scaling applies to all class-specific matrices for the whole branch. All classes and branch modes within one mixture must use the same scaling type.

Branch heterogeneity

GY94BranchModel() uses the same mode-tree convention as other PalantiR time-heterogeneous models. Mode labels are zero-based indices into models, and start_mode specifies the process at the root. Any GY94 parameter may change between modes, including omega, kappa, and the stationary codon frequencies.

mammals <- Phylogeny("../inst/extdata/mammals.newick")
mammal_modes <- Phylogeny("../inst/extdata/mammals_switch.newick", type = "mode")

background <- GY94(omega = 0.2, kappa = 2, frequencies = pi_1x4)
foreground <- GY94(omega = 3.0, kappa = 2, frequencies = pi_3x4)

episodic <- GY94BranchModel(
    models = list(background, foreground),
    mode_phylogeny = mammal_modes,
    start_mode = 0)

When stationary frequencies change, descendants enter the new process with the state distribution inherited from the parent; they are not redrawn from the new equilibrium. PalantiR simulates the resulting transient and uses the requested temporal rescaling method to preserve each component’s rate relative to the common mixture denominator. The mixture-weighted rate therefore retains the chosen branch-length currency. The exact time change is the default; the segment approximation remains available through rescale_method = "segments".

Arbitrary discrete site classes

GY94SiteModel() accepts any collection of homogeneous GY94 models and time-heterogeneous GY94BranchModel specifications. Supply either n_sites, one count per class, or site_classes, one assignment per output position. Character assignments refer to class labels; numeric assignments are one-based indices into models. Branch-heterogeneous classes in one mixture must use the same mode phylogeny and number of branch models, because their generators are normalized together on each branch type. Output site and class_index fields are zero-based.

site_model <- GY94SiteModel(
    models = list(
        constrained = background,
        neutral = GY94(omega = 1, kappa = 2, frequencies = pi_1x4),
        episodic = episodic),
    n_sites = c(40, 30, 30),
    labels = c("constrained", "neutral", "episodic"))

For an explicit interleaved assignment, use for example site_classes = c("constrained", "episodic", "neutral", ...); the simulator restores that requested order after efficiently simulating sites in class blocks.

Zhang–Nielsen–Yang branch-site model

ZNYBranchSiteModel() builds the four standard branch-site classes as a special case of the general site-class representation:

Site class	Background omega	Foreground omega
0	`omega0`	`omega0`
1	1	1
2a	`omega0`	`omega2`
2b	1	`omega2`

The mode tree must label background branches 0 and foreground branches 1. site_counts gives counts in the order 0, 1, 2a, 2b. Foreground frequencies default to the background frequencies, but may be supplied separately when the stationary codon distribution also changes on foreground branches.

branch_site <- ZNYBranchSiteModel(
    mode_phylogeny = mammal_modes,
    site_counts = c(4, 4, 4, 4),
    omega0 = 0.2,
    omega2 = 3,
    kappa = 2,
    frequencies = c(T = 0.30, C = 0.20, A = 0.25, G = 0.25),
    frequency_model = "F1x4",
    scaling_type = "standard")

ZNYBranchSite() is a short alias.

Simulation and output metadata

The simulator samples a root state for each site from the equilibrium of that class’s root model unless a codon Sequence is supplied explicitly.

sim <- simulate_gy94_site_model(mammals, branch_site)

head(sim$site_classes)
head(sim$substitutions)
plot(sim$alignment)
plot(sim, sites = 0:9)

For the class counts above, sim$site_classes begins:

  site site_class class_index
1    0          0           0
2    1          0           0
3    2          0           0
4    3          0           0
5    4          1           1
6    5          1           1

The complete assignment table is stored in sim$site_classes, including sites with no substitutions. The same table is attached to sim$alignment; the interactive alignment header displays the class label. Every event in sim$substitutions carries site_class and class_index, which also appear in substitution-history tooltips. sim$site_model retains the general class specification and sim$site_class_models maps labels to their homogeneous or branch-heterogeneous model definitions. The generating GY94SiteModel also records class_weights, the mixture_scaling denominator for each branch mode, and stationary_class_rates relative to that denominator.

References

Goldman, N. and Yang, Z. (1994). A codon-based model of nucleotide substitution for protein-coding DNA sequences. Molecular Biology and Evolution 11, 725–736.
Yang, Z., Nielsen, R., Goldman, N. and Pedersen, A.-M. K. (2000). Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics 155, 431–449.
Yang, Z. and Nielsen, R. (2002). Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Molecular Biology and Evolution 19, 908–917.
Zhang, J., Nielsen, R. and Yang, Z. (2005). Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Molecular Biology and Evolution 22, 2472–2479.

Goldman–Yang (GY94) codon models