zombi2.sequences¶
Level 3: a sequence evolving inside a gene, along its gene tree. The level takes a whole genome run, not bare gene trees — a sequence sees the species tree only through its gene tree, but its rate is set by the species branch the gene sits on.
simulate_sequences(genomes, *, model: SubstitutionModel | None = None, length: int | None = None, partitions=None, profiles=None, families=None, intergene_model: SubstitutionModel | None = None, intergene_speed=3.0, insertion=0.0, deletion=0.0, insertion_extent=3.0, deletion_extent=3.0, substitution=None, divergence=None, seed=None, parallel=False, stream_to=None, outputs=None, flat: bool = False, genes=None, joint: bool = False, record: bool = False, progress=False) -> 'SequencesResult | StreamedSequences'
Evolve one sequence down each family's gene tree under a substitution model.
genomes is a genome run — the FamilyGenomesResult that
genomes.simulate_genomes_family(...) returned. Its gene_trees are what the sequences
evolve along and its complete_tree is the species tree the lineage clock rides; bare gene
trees are rejected (they would run, but with no clock and no species phylogram — a silent
degradation). Each family's complete gene tree is evolved, so the true history is complete and
ancestral sequences exist for extinct/lost lineages too; the observable alignments are the
extant tips.
model is a substitution model from the menu (substitution_models) — nucleotide
jc69 · k80 · hky85 · gtr, or protein poisson · jtt · dayhoff ·
wag · lg; its alphabet is what the sequences are written in (ACGT or the 20 amino
acids). length is the number of sites. substitution is the per-site substitution
rate (default 1.0): a branch of Δt time accrues substitution · Δt substitutions/site.
The founding sequence of each family is drawn from the model's stationary frequencies. Deterministic
given seed.
insertion and deletion give the sequences indels, so an alignment gains gaps: a
lineage loses sites it had, or gains sites the others never did. They are the family and ordered
resolutions' indels — a nucleotide genome owns its own, because there a base pair has a
position, and passing these alongside one is refused rather than merged. The rates are
relative to substitution: deletion=0.05 is five deletions for every hundred substitutions
a site expects, so a lineage's clock reaches its indels too and the number means the same on a
tree of any height. insertion_extent / deletion_extent are how many sites an event takes,
a bare number being the mean of a geometric draw; any distribution works, so Fixed(1) is a
single-site indel and a power law is the shape indel lengths really have.
The columns are the union of every site any lineage ever had, and a lineage that carries none of
one shows a gap there — the same alignment a nucleotide run writes, arrived at by the level that
knows the sites. Refused beside partitions or profiles, which are written against a site
count an indel changes.
substitution may carry a lineage clock — one factor per species branch, shared across
families, computed once before evolving, rescaling each gene-tree branch by the clock of the species
branch it sits on: PerSite(1.0).varying_among('lineages', LogNormal(0.0, σ)) is the
uncorrelated clock (each branch drawn
i.i.d.), and PerSite(1.0).varying_among('lineages', Drift(LogNormal(0.0, σ))) is the
autocorrelated clock (the factor drifts
parent→child down the species tree).
It may also carry a driver — PerSite(1.0).scaled_by(habitat, {"cave": 0.5, "surface": 1.0}),
where habitat is a trait grown first (the ~zombi2.traits.TraitsResult, or the path to the
trait_events.tsv it wrote). That is conditioning, not a joint run: SPEC §3 allows the pair
Traits–Sequences to be conditioned and never joined, so naming a live level ("trait") raises
rather than starting one. A clock and a driver compose — the factors multiply (SPEC §5), so a
lineage's dealt tempo and its state both count — and several drivers on one rate multiply too.
A discrete driver switches mid-branch, and the branch length is the driver integrated across
the branch rather than one sample of it (clock), so the phylograms are the trees the alignments
were actually drawn along. Any other modifier (the Markov clock, a draw among
families), a second lineage clock, or a non-PerSite scope raises.
record=True keeps this level's own history — every substitution, and every site gained
or lost — in result.events, written as sequence_events.tsv. It is off by default and is
the only level's log that is, because it is the one log bigger than the output it explains: three
hundred sites over thirty time units at rate 1.0 is nine thousand rows for one family. Recording
also changes the sampler: an ordinary run draws each branch's end from exp(Q·bl) and never
simulates the path between the two ends, while a recorded run walks that path site by site. The
two are the same process and the same distribution at a branch's end, so a recorded run is a
valid run — but it is a different realisation for the same seed, which is the price of asking
what happened rather than only where it ended. A site is named by an id rather than a position,
because a position moves with every insertion above it; see zombi2.sequences._record.
families restricts the run to named families — families=["chaperone"] — instead of
evolving every one. The names are the genome run's, declared there with
families=[family("chaperone")]. Two uses: a run of one family, whose pooled composition is
then that family's composition, and a pair of runs that let one gene's sequence drive another's
rate. That pair needs it: without it the second run re-evolves the driver alongside its target,
so the rate reads one history of the driver while the file on disk holds a different one.
Rate variation across sites rides on model, not on substitution:
model=hky85(2.0).across_sites(gamma_shape=0.5, invariant=0.1) sorts the sites into a
discretised-Gamma set of rate classes plus a class that never changes. The two axes are
orthogonal and compose — the clock says which lineages run fast, the model which sites do.
partitions splits a family's sites into blocks each under its own model, in place of
model= and length=::
partitions=[(hky85(kappa=2.0), 600), (jc69(), 400)]
— a 1000-site gene whose first 600 sites evolve under HKY85 and whose last 400 evolve under
JC69, concatenated in that order into one sequence per gene copy. Each partition's model may
carry its own across_sites classes, which is how the field usually spells a codon-position
split. Giving model= or length= alongside is refused rather than merged: the partitions
already answer both, and the length is their site counts summed. Every partition must be over the
same alphabet, because they concatenate into one sequence.
All the partitions share the run's one substitution rate, and so the family keeps its one
phylogram — exactly, not approximately: every model is normalised to one expected substitution
per site per unit branch length, and every set of rate classes to a mean of 1, so a branch of
Δt accrues rate · Δt substitutions per site in each partition alike. There is
deliberately no per-partition speed; it would make one phylogram a weighted average of the trees
the partitions were really drawn along. Family and ordered runs only — a nucleotide run's
blocks already carry their own lengths and their own gene/spacer models, so it refuses
partitions. Experimental (SPEC §9): Python API, no CLI flag yet.
profiles gives chosen families a site profile — an (L, K) array whose row i is the
equilibrium frequencies at position i, in the model's own alphabet order::
profiles={4: array_of_shape_300_by_20}
Every model on the menu gives a gene one set of frequencies shared by all its sites; a profile
says which residues belong at each position instead, which is what a buried hydrophobic site
and the loop beside it differ by. The base model's exchangeabilities are kept — which pairs
interchange easily is chemistry, and the profile's business is only where each residue belongs —
so the site's matrix is the run's model over the site's own frequencies. Families you leave out
evolve under model untouched, and a family's profile changes nothing about any other family.
A row of exact zeros is refused: it makes that site's matrix degenerate, and a real profile says
a residue is unlikely at a position rather than impossible — add a pseudocount. A flat
profile, every row the model's own frequencies, is the model without one; it is statistically
identical rather than byte-identical, because L single-site models consume the random stream
differently from one L-site model. Profiles compose with across_sites — a profile says which
residues, +Γ says how fast — and are refused alongside partitions (both decide a
family's per-site models) and alongside parallel (which ships one shared partition set to
every worker).
An amino-acid profile needs a protein model and so belongs to a family or ordered run: a nucleotide genome is measured in base pairs and read on either strand, so it refuses protein models outright. Profiles still apply there, over the four bases — a row per base pair — and the row count must equal that block's length, which the genome run already fixed. Experimental (SPEC §9): Python API, no CLI flag yet.
On a nucleotide genome run every root block is evolved — spacer as well as genes — each at its
own length in bp, so length does not apply and is rejected. model evolves the genes and
intergene_model (default jc69) the spacer, at intergene_speed times the rate (default
3.0). Each carries its own across-site variation: decorating model with
across_sites does not reach the spacer, whose default jc69() stays flat — the spacer's
job is to be the unconstrained null, and silently giving it the genes' Gamma would make it
something else. Give intergene_model a decorated model to vary the spacer too. Because the whole genome is covered, the run also puts the genomes back together:
.node_genomes holds every node's chromosomes, blocks concatenated in physical order — the complete
tree, reconstructed — and .initial_genome the one the run started with.
The result carries the phylograms the sequences were drawn along — each gene tree and the
species tree, with branch lengths converted from time to substitutions/site by the same
base × clock × Δt.
parallel opts into evolving the gene trees concurrently — one gene tree per worker
process, which is where a run's time goes when the sequences are long or a nucleotide genome has
thousands of blocks. False (the default) runs the serial engine above. True uses every
core; a positive int sets the worker count. It is a separate engine: each family draws
from its own RNG stream (spawned from seed), so every worker count returns the same bytes,
but that realisation differs from the serial one for a given seed — parallel is a speed choice,
made once, not a drop-in for the default. Threads would not help here (numpy releases the GIL too
little for these array sizes), so this is process-backed and only pays off above a work threshold.
Because it spawns processes, a script that calls it with parallel set must guard its entry with
if __name__ == "__main__": (the standard multiprocessing requirement); the zombi2 CLI
already does, so --parallel there needs nothing extra.
stream_to=DIR writes each family's files as it finishes and keeps nothing, returning a
StreamedSequences handle instead of a SequencesResult. This is the level where a run's memory
actually goes — every alignment and every ancestral sequence live at once — so it is the dial for
a run whose result would not fit. The files are the same ones .write(DIR) would leave, so a
streamed run and an in-memory one at the same seed are the same dataset; outputs picks which,
exactly as .write does, and flat is passed through the same way. It composes with
parallel. A nucleotide run cannot stream: it puts whole genomes back together, and that
needs every block's sequence at once, which is the opposite of keeping nothing.
Source code in zombi2/sequences/__init__.py
1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112 2113 2114 2115 2116 2117 2118 2119 2120 2121 2122 2123 2124 2125 2126 2127 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 2140 2141 2142 2143 2144 2145 2146 | |
dataclass
¶
SequencesResult(alignments: 'Mapping[int, dict[str, str]]', ancestral: 'Mapping[int, dict[str, str]]', founding: dict[int, str], phylograms: dict[int, dict[str, str | None]], species_phylogram: dict[str, str | None], seed: int | None, node_genomes: 'Mapping[str, dict[int, str]]' = dict(), initial_genome: dict[int, str] = dict(), unit: str = 'family', alphabet: str = '', extant_tips: tuple[str, ...] = (), families: tuple[str, ...] = (), events: list = list(), _insertions: dict = dict(), _raw_rows: dict = dict(), _raw_ancestral: dict = dict())
What simulate_sequences() returns.
alignments—{family: {n<species>_g<copy>: sequence}}: the observable gene alignment, one entry per extant gene-tree tip, keyed by its (unique, per-segment) gene id — the same labels as the gene tree's / phylogram's Newick leaves. Empty for a family with no surviving copy.ancestral—{family: {n<species>_g<copy>: sequence}}: the true sequence at every node that is not an extant tip — internal nodes (the family's root gene included) and the dead tips, where a copy was lost or its species went extinct. Withalignmentsit accounts for every node of the tree exactly once, so every label in the complete phylogram names a sequence.founding—{family: sequence}: the sequence the family started with, at its origination — the state the phylogram's root branch leads from. It is drawn from the model's stationary frequencies and then evolves across the stem into the root gene's sequence, so it is not the same string asancestral[family]["g<root copy>"]unless the stem is empty. Kept out ofancestralon purpose: those keys pair one-to-one with phylogram nodes, and the origination is a point on a branch, not a node.phylograms—{family: {"complete": newick, "extant": newick | None}}: each gene tree with branch lengths in substitutions/site (base × lineage-clock × Δt) — the ground-truth tree behind each alignment. Every node is labelledn<species>_g<copy>, so the tips match thealignmentskeys and the internal nodes match theancestralkeys (the phylogram pairs one-to-one with the sequences)."extant"isNonefor a family with no survivor.species_phylogram—{"complete": newick, "extant": newick | None}: the species tree with branch lengths in substitutions/site — the molecular clock made visible (which lineages ran hot / cold). Always present: a run always comes from a genome run, which carries its tree.node_genomes—{lineage: {chromosome id: sequence}}: every node's assembled genome, its blocks concatenated in physical order (reverse-complemented where the genome carries them inverted) — extant tips, ancestors and the lineages that went extinct alike. The same coverage, and the same name, as the genome level'snode_genomes.genomesis the extant tips alone — the observed genomes — exactly as it is one level down. Only a nucleotide genome run has any: a family or ordered run has gene families, not coordinates, so there is no genome to lay out and both are empty.initial_genome—{chromosome id: sequence}: the genome the run started with, at the root lineage's origination. Not ingenomes, because it belongs to no node: the root branch is real simulated time, so the root node's genome is this one plus whatever happened along the stem. It stands togenomesasfoundingstands toancestral.seed— the run's seed.unit— what the integer key ofalignments/ancestral/founding/phylogramsnames:"family"(a gene family id) on a family or ordered run,"block"(an index into the genome run'sroot_blocks) on a nucleotide one, where every block evolves and spacer has no family. They are different numbering schemes over the same ints, so a gene family id is not a key here on a nucleotide run — go throughblock_of(). It is also what the filenames say.alphabet— what the sequences are written in:"ACGT", or the 20 amino acids. One per run (every partition of a family shares it), and whatcompositionchecks the letters it is asked to count against.
property
¶
The observed genomes — the assembled genome of each extant tip, keyed by tip name
(n5): what you would have sequenced. Empty unless the run came from a nucleotide
genome run, which is the only kind that has coordinates to lay out.
node_genomes is the whole record — every node, ancestors and extinct lineages included.
The pair is the genome level's ~zombi2.genomes.FamilyGenomesResult.genomes /
node_genomes, and it reads here the way alignments reads against ancestral:
what you observe, and the truth behind it.
The share of a lineage's sequence that is one of letters, at every instant, as a
conditioning driver (~zombi2.sequences._composition.Composition)::
proteins = simulate_sequences(g, model=lg(), length=300, seed=1)
simulate_discrete(tree, states=["mesophile", "thermophile"], start="mesophile", seed=2,
switch=PerLineage(0.2).scaled_by(proteins.composition("KR"),
Curve(lambda x: 40.0 ** (x - 0.1))))
This is how an amino-acid frequency is asked for: one residue ("K") or a set of them
("KR", "AVLIMFWP"). The letters must be in this run's alphabet; gc is the same
driver over "GC", named because it is the one people ask for by name.
A number, so it takes a ~zombi2.params.mapping.Curve or a ~zombi2.params.mapping.Scalar, and
it drives what comes after a sequence — a trait, or a further sequence run — never the
genome the gene trees came from.
absent is what a branch reads where the run has no sequence at all on that lineage,
and it is what makes one family's composition a usable driver: restrict the run with
families=["chaperone"] and this statistic is that family's, but the family is missing on
some branches and a driver has to answer for every branch the target walks::
chaperone = simulate_sequences(g, families=["chaperone"], model=lg(), length=300, seed=3)
chaperone.composition("KR", absent=0.08)
On a restricted run with a gap and no absent, resolving the driver raises rather than
carrying the parent's value forward, which would be a different model from the one asked
for.
Source code in zombi2/sequences/__init__.py
This run's GC content as a conditioning driver: composition over "GC", the
fraction of a lineage's DNA that is G or C, pooled over every family the run evolved::
seqs = simulate_sequences(g, model=hky85(2.0), length=300, seed=1)
simulate_continuous(tree, rate=PerLineage(1.0).scaled_by(seqs.gc(),
Curve(lambda x: 4.0 * x)),
seed=2)
Nucleotide runs only, because G and C are also glycine and cysteine: on a protein run the
call is ambiguous rather than wrong, so it is refused and composition asked for instead.
family is here only to refuse one — one family's GC is asked for by restricting the run
to it with families=, and absent then says what a branch without it reads
(composition).
Source code in zombi2/sequences/__init__.py
What this run produced, as a plain dict — the payload of sequences_summary.json.
The identity is the number worth having: it is what says whether the alignments carry signal, and it depends on the height of the tree the run went down, which no flag shows you. The run prints it and warns when it is near the floor; this is the machine-readable copy of that.
Source code in zombi2/sequences/__init__.py
write(directory, outputs=('alignments', 'phylograms', 'species_phylogram', 'genomes', 'initial_genome', 'summary'), *, flat: bool = False) -> None
Write chosen outputs to directory (created if needed). <u> below is
fam<family> on a family or ordered run and block<index> on a nucleotide one — the
integer keys mean different things, so the files say which (see unit):
"summary"→sequences_summary.json: what came out (summary); the ancestral count is dropped when"ancestral"was not asked for."alignments"→<u>.fastaunderalignments/(skipped for empty families)."ancestral"→sequences_ancestral_<u>.fastaunderancestral/."founding"→sequences_founding.fasta, one record<u>apiece: the sequence each family originated with, before its stem."phylograms"→phylogram_<u>_{complete,extant}.nwk(subs/site) underphylograms/."species_phylogram"→clock_species_tree_{complete,extant}.nwk: the species tree with its branches in substitutions/site — the molecular clock made visible."genomes"→genome_<lineage>.fastaundergenomes/, one file per node — extant, extinct and ancestral alike — with one record per chromosome. Nucleotide runs only; nothing is written otherwise. The big one: a real genome times every node in the tree."initial_genome"→genome_initial.fasta, ingenomes/with the rest: it is a whole-genome FASTA like they are, and it belongs beside them.
Everything that is one file per family or per node gets a subdirectory, or the two trees and
the one founding FASTA would be lost among thousands; flat=True writes everything into
directory instead. Nothing is created for an output this run has none of, so a family run
leaves no empty genomes/ behind.
Source code in zombi2/sequences/__init__.py
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 | |
dataclass
¶
StreamedSequences(directory: str, seed: int | None, n_families: int, n_sequences: int, outputs: tuple, identity: 'float | None' = None, sites: 'int | None' = None, n_ancestral: int = 0, unit: str = 'family')
A sequence run written straight to disk, family by family — what stream_to= returns.
Thin on purpose: the outputs are the files, so this carries where they are and how big the run
was, not the run itself. The sequence level is where a run's memory actually goes — every family's
alignment plus its ancestral sequences, all live at once — so this is the handle for the size at
which a SequencesResult would not fit.
The same payload SequencesResult.summary() builds, from what the sink counted as it wrote.
A streamed run has to produce the same files as an in-memory one at the same seed — that is the whole contract — so it produces this one too. Two fields it cannot know: the ancestral count and the assembled genomes, neither of which a streamed family run writes.
Source code in zombi2/sequences/__init__.py
Mean identity over every within-family sequence pair, or None when no family holds two
sequences to compare.
Exhaustive, and the same quantity however the run was written: it used to be a bounded random
sample, and the streamed path sampled differently, so one flag that only chooses how much memory
a run uses moved a number in the report while leaving every alignment byte-identical. Nothing
said the number was an estimate, so it read as a fingerprint two matching runs could disagree on.
_identity_counts makes counting them all linear in the alignment, so there is no longer a reason
to estimate.
Source code in zombi2/sequences/__init__.py
The substitution-model menu¶
A substitution model is a K×K rate matrix Q, normalised to one expected substitution per site
per unit branch length, and its stationary frequencies. Different models are genuinely different
matrices, so they stay a menu of constructors rather than one grammar.
Substitution models — the menu.
A substitution model is the chemistry of a sequence: a K×K rate matrix Q (normalised to
one expected substitution per site per unit branch length) and its stationary frequencies π.
Different models are genuinely different matrices — Jukes–Cantor, K80, HKY85, GTR differ in their
transition/transversion structure and base composition — so, unlike the clock, they do not
collapse to one grammar: they stay a menu of constructors, each taking its own physical parameters.
Two alphabets are on the menu: the four nucleotide models (4 states, ACGT — jc69() ·
k80() · hky85() · gtr()) and the five protein models (20 states,
AMINO_ACIDS — poisson() · jtt() · dayhoff() · wag() · lg()).
The protein models are empirical: their exchangeabilities and frequencies were estimated once from
large alignments and are read off the published matrices (_aa_matrices), so they take no
free parameters — you pick one, you do not tune it. Codon models are not in the menu; adding one is
a pure extension of it, no refactor.
A matrix of your own goes through reversible(): a symmetric exchangeability matrix S and
frequencies π, normalised the same way, over any alphabet. It is the door the menu itself uses,
so gtr() is its four-state case and lg() its twenty-state one. It takes S and π rather
than Q because the engine requires a time-reversible model and S·π cannot be anything
else; a general Q is refused rather than evaluated wrongly (see SubstitutionModel).
Across-site rate variation decorates any model on the menu rather than adding entries to it:
SubstitutionModel.across_sites() returns the same chemistry with its sites sorted into
rate classes — a discretised Gamma (+Γ), a class that never changes (+I), or both. That is
where the field puts it (HKY85+I+G4) and, more to the point, it is not a rate modifier:
SPEC §5's modifiers multiply one rate by a context factor, while this splits the sites of one
rate into classes. Putting it on the model keeps the rate grammar to one job and lets the spacer of
a nucleotide run carry its own (or no) variation.
Every model here is time-reversible, so the transition matrix over a branch of length t (in
substitutions/site), P(t) = exp(Q·t), is computed by eigendecomposition of the symmetric
matrix B = diag(√π)·Q·diag(1/√π) (numpy only, no scipy):
P(t) = diag(1/√π)·V·exp(Λt)·Vᵀ·diag(√π).
dataclass
¶
SubstitutionModel(name: str, Q: ndarray, stationary: ndarray, alphabet: str = BASES, site_rates: tuple[float, ...] = (1.0,), site_shares: tuple[float, ...] = (1.0,))
A K-state reversible model: a normalised K×K rate matrix Q, its stationary
frequencies, and the ordered alphabet whose order Q / stationary follow.
Built through the menu constructors (jc69(), k80(), hky85(), gtr()) or, for a matrix of
your own, through reversible() — which takes a symmetric exchangeability matrix and the
frequencies rather than Q itself, so a non-reversible model cannot be spelled. The reversible
eigendecomposition behind p_matrix() is precomputed once in __post_init__, and so is the
check that the matrix it is being applied to is reversible at all: the transform is only a
similarity transform of Q under detailed balance, so a matrix that violates it would give
transition probabilities that are wrong without looking wrong. That check lives here rather than
only in the constructors because this class is public, and every route into it has to end up in a
state the engine can honour.
site_rates and site_shares carry across-site rate variation: the sites of a sequence
are sorted into classes, site_rates[c] multiplying the branch length for the sites in class
c and site_shares[c] giving the proportion of sites that fall there. A plain model has
the one class (1.0,) / (1.0,) and every site evolves at the same speed; across_sites()
builds the rest. The two are dimensionless — a class scales a branch length, it is not a rate in
SPEC §5's time⁻¹ sense — and they obey one invariant:
``sum(rate × share) == 1``
which is what makes a branch length the mean substitutions per site. Every phylogram in the result is drawn with that meaning, so the invariant is checked here rather than trusted.
Transition probabilities over branch length t (substitutions/site).
P(t) = exp(Qt) via the reversible eigendecomposition; clipped to [0, ∞) to scrub tiny
negative round-off so every row is a valid probability distribution.
Source code in zombi2/sequences/substitution_models.py
across_sites(*, gamma_shape: float | None = None, invariant: float = 0.0, rate_categories: int = 4) -> 'SubstitutionModel'
The same model with its sites sorted into rate classes — +Γ, +I, or both.
Every site of a gene evolving at exactly one speed is a model no real gene obeys: some positions are held nearly fixed by function, others race ahead. Two standard knobs say so, and they compose:
gamma_shape(α) draws each site's rate from a Gamma with mean 1, discretised intorate_categoriesequal-probability classes (Yang 1994; see_site_rates). A small shape is strong variation — 0.5 gives a few fast sites among many slow ones — and a large one is nearly flat. There is no default: asking for+Γmeans choosing how unequal the sites are.invariant(0 by default) is the proportion of sites that never change: one more class, at rate 0. Real alignments have columns that are constant because the site cannot change, not because it happened not to, and a Gamma alone fits those badly.
A site's class is drawn once and holds for the whole tree — that is what across-site
variation means, as against a rate that varies along a branch. Sites stay independent, so
this is a knob on the existing engine and not a new one (SPEC §9).
The classes are normalised to mean 1, invariant sites included: the Gamma classes are
scaled by 1 / (1 - invariant) to make room for the sites that never change. So a branch
length stays the mean substitutions per site, a run with variation and one without have
the same phylogram at the same rate, and only the spread of change across columns differs.
This is not a rate modifier (SPEC §5), and simulate_sequences will still refuse one
on substitution: a modifier multiplies one rate by a context factor, while this splits
one rate's sites into classes. It belongs to the model, which is where the field puts it —
the decorated name reads HKY85+I+G4, and it is what the run reports.
Returns a new model; the original is unchanged (they are frozen).
Source code in zombi2/sequences/substitution_models.py
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 | |
Jukes–Cantor (1969): equal rates, equal base frequencies. No free parameters.
Kimura 2-parameter (1980): transition/transversion ratio kappa, equal frequencies.
HKY85 (Hasegawa–Kishino–Yano 1985): transition bias kappa with unequal base frequencies (A,C,G,T).
Source code in zombi2/sequences/substitution_models.py
gtr(exchangeabilities=(1, 1, 1, 1, 1, 1), frequencies=(0.25, 0.25, 0.25, 0.25)) -> SubstitutionModel
General time-reversible: 6 exchangeabilities [AC,AG,AT,CG,CT,GT] and frequencies
(A,C,G,T). This is reversible() restricted to four states, and the two name the numbers the
same thing on purpose — they are the same numbers.
Source code in zombi2/sequences/substitution_models.py
reversible(exchangeabilities, frequencies, *, name: str = 'Custom', alphabet: str = BASES) -> SubstitutionModel
A model of your own matrix: a symmetric K×K exchangeability matrix and K
stationary frequencies, over any alphabet.
Q_ij = S_ij · π_j for i ≠ j, the diagonal set from the row sums, and the whole matrix
scaled so that -Σ π_i Q_ii = 1 — one expected substitution per site per unit branch length,
exactly as every model on the menu is scaled. That normalisation is what keeps a branch length
meaning substitutions per site whatever matrix produced it, so a phylogram from a model built
here is comparable with one from hky85() without a conversion.
This is the same door the menu already goes through, so the menu is its special cases: gtr() is
this with K = 4 and the six nucleotide exchangeabilities spelled out, and lg() is this with
K = 20 and a published triangle. Give alphabet=AMINO_ACIDS for a 20-state protein matrix
of your own, or any K distinct ASCII characters for an alphabet that is neither — though the
rest of the level only makes sense for those two: a nucleotide genome run reads its blocks on
both strands and so refuses anything but ACGT, and the FASTA it writes is whatever letters
you chose.
You give S and π rather than Q, and that is the point. The engine's p_matrix() computes
exp(Qt) by eigendecomposing diag(√π)·Q·diag(1/√π), which is a symmetric matrix similar to
Q only when π_i Q_ij = π_j Q_ji — detailed balance. A symmetric S times π obeys
that by construction, so there is no way to spell a model here that the engine would evaluate
wrongly. A general Q handed straight to SubstitutionModel could be, which is why the
dataclass checks the same identity and refuses. Non-reversible models (UNREST and its relatives)
are not implemented — a statement about this code, not about the models (SPEC §5).
S must be symmetric, non-negative, and zero on the diagonal; π strictly positive and
summing to 1. name is what the run prints and what error messages quote, so it is worth
setting to whatever the matrix actually is.
There is deliberately no command-line flag for this. SPEC §5 asks for one written form
everywhere, and a K×K matrix has no written form a command line can carry: at K = 20 it
is 190 numbers, which is a file format, not an argument. It stays a Python constructor; the CLI's
--model offers the menu, which is what a flag can honestly hold.
Experimental (SPEC §9): Python API only, and the spelling may still move.
Source code in zombi2/sequences/substitution_models.py
Poisson: equal exchangeabilities, equal frequencies — the JC69 of proteins. No free parameters.
Source code in zombi2/sequences/substitution_models.py
JTT (Jones, Taylor & Thornton 1992): the empirical matrix from close protein homologues.
Dayhoff (Dayhoff, Schwartz & Orcutt 1978): the original PAM matrix, in PAML's values.
WAG (Whelan & Goldman 2001): estimated by maximum likelihood over a wide protein database.
LG (Le & Gascuel 2008): WAG's successor, fitted with across-site rate variation — the default protein model of modern phylogenetics.
Map an array of integer states back to a string over alphabet — ACGT by default, or
AMINO_ACIDS for a protein model (callers pass model.alphabet).
states are indices into alphabet, so the whole array is one numpy gather into an ASCII
lookup table — lut[states] — read out in a single .tobytes().decode() rather than one
Python step per site. This is called once per node of every gene tree, so the per-site loop it
replaces was the dominant cost of a sequence run; the result is byte-for-byte the same string. The
lookup table itself is cached per alphabet (_DECODE_LUT), built once instead of per call.
states may be multi-dimensional: a 2-D (rows, length) array decodes to the rows' strings
concatenated back to back (row-major), which lets a caller decode a whole gene tree's nodes in one
gather + one ASCII decode and slice the fixed-length rows out — see _split().
Source code in zombi2/sequences/substitution_models.py
The inverse of decode(): a string over alphabet to its integer states. Used to found a
run's blocks from a real fasta= — the supplied DNA becomes a block's founding states. A character not in
alphabet raises (the FASTA reader already rejects non-ACGT, so this is a second guard).