Update documentation for CLI command renames and new commands
Renames the reindex command to convert across all usage guides and navigation menus. Adds documentation for the newly introduced name-tree command. Updates flag references, such as replacing --metric with --distance in pack and phylo modules, and adjusts MkDocs configuration to match the revised structure. All changes are strictly limited to documentation and configuration files.
This commit is contained in:
Eric Coissac
@@ -20,7 +20,7 @@ An index directory is organized as `KmerIndex → partitions → layers`, with a
|
||||
|
||||
## Parallel execution and NUMA awareness
|
||||
|
||||
Partition-level work (index construction, `merge`, `filter`, `reindex`, `select`, `phylo`'s sibling-annex/Sankoff computations) is dispatched by a partition runner that adapts to the machine's memory topology, detected automatically at startup via hwloc:
|
||||
Partition-level work (index construction, `merge`, `filter`, `convert`, `select`, `phylo`'s sibling-annex/Sankoff computations) is dispatched by a partition runner that adapts to the machine's memory topology, detected automatically at startup via hwloc:
|
||||
|
||||
- On a multi-socket / multi-NUMA-node machine, one thread pool is pinned per NUMA node, and each partition is processed entirely by threads pinned to one node — keeping the memory a partition touches local to that node's DRAM. This matters because touching kmer data across NUMA nodes without pinning can degrade throughput by an order of magnitude or more on large multi-socket machines.
|
||||
- On a single-socket machine, Apple Silicon, or if hwloc cannot report NUMA topology, all cores are treated as one node with no pinning and negligible overhead — this is the default behavior on macOS.
|
||||
|
||||
@@ -19,7 +19,7 @@ Each partition's surviving kmers are mapped to a dense range of integer slots by
|
||||
|
||||
## Evidence: exact vs. approximate
|
||||
|
||||
Two verification modes are available, selected at build time (`index --approx`) and convertible afterwards ([`reindex`](../usage/reindex.md)):
|
||||
Two verification modes are available, selected at build time (`index --approx`) and convertible afterwards ([`convert`](../usage/convert.md)):
|
||||
|
||||
- **Exact** (default): the hashed slot stores a pointer back into the partition's unitig data. At query time the kmer is reconstructed from that location and compared directly to the query. Zero false positives, at the cost of one extra random read per lookup.
|
||||
- **Approximate** (`--approx`): the slot stores a short fingerprint (`--evidence-bits` bits) instead of a pointer; verification is a single fingerprint comparison. This trades a small, bounded false-positive rate ($1/2^b$ per kmer, reduced further to about $1/2^{b \cdot z}$ for a read requiring $z$ consecutive matching kmers via the `-z`/`--findere-z` parameter) for lower memory and disk usage, since no reconstruction index is needed. See [`estimate`](../usage/estimate.md) to explore this trade-off before building.
|
||||
@@ -51,6 +51,6 @@ Two verification modes are available, selected at build time (`index --approx`)
|
||||
|
||||
A **layer** corresponds to one increment of kmer content added to a partition — most commonly, one [`merge`](../usage/merge.md) operation that introduces kmers not already present in the index. Genomes already present in the index simply gain new columns in the existing layers' count/presence data; only genuinely new kmer content is assembled into a new layer. Because of this, merging cost scales with the novel kmer content being added, not with the accumulated size of the index. A query against an index with several layers checks each layer's MPHF in turn.
|
||||
|
||||
Sources merged together must share the same kmer size, minimizer size, partition count, and evidence mode (including matching approximate-mode parameters); mismatches are rejected rather than silently reconciled — [`reindex`](../usage/reindex.md) one of the sources first if needed.
|
||||
Sources merged together must share the same kmer size, minimizer size, partition count, and evidence mode (including matching approximate-mode parameters); mismatches are rejected rather than silently reconciled — [`convert`](../usage/convert.md) one of the sources first if needed.
|
||||
|
||||
`obikmer pack` consolidates a partition's per-column files (counts/presence) into a single file, reducing the number of file opens needed at query time.
|
||||
|
||||
+2
-3
@@ -24,11 +24,10 @@ All functionality is exposed through a single binary, `obikmer`, organized as su
|
||||
| [`query`](usage/query.md) | Query an index with sequences and annotate matches |
|
||||
| [`dump`](usage/dump.md) | Dump indexed kmers as CSV |
|
||||
| [`annotate`](usage/annotate.md) | Add, update, or dump genome metadata |
|
||||
| [`phylo`](usage/phylo.md) | Compute pairwise evolutionary-distance proxies, trees, and phylogenetic exports |
|
||||
| [`name-tree`](usage/name-tree.md) | Translate a TNT/PhyG numeric-label tree export back to real taxon names |
|
||||
| [`phylo`](usage/phylo.md) | Compute pairwise genome distances, trees, and phylogenetic exports |
|
||||
| [`unitig`](usage/unitig.md) | Dump the unitigs of an index as FASTA |
|
||||
| [`estimate`](usage/estimate.md) | Estimate approximate-index parameters before indexing |
|
||||
| [`reindex`](usage/reindex.md) | Convert an index's evidence representation (exact ↔ approximate) |
|
||||
| [`convert`](usage/convert.md) | Convert an index's evidence representation (exact/approximate/hybrid), in place |
|
||||
| [`utils`](usage/utils.md) | Miscellaneous index maintenance and inspection utilities |
|
||||
| [`pack`](usage/pack.md) | Pack per-column matrix files into a single-file format |
|
||||
|
||||
|
||||
@@ -0,0 +1,29 @@
|
||||
# convert
|
||||
|
||||
Convert an existing index's evidence representation in place, between exact, approximate, and hybrid.
|
||||
|
||||
```bash
|
||||
obikmer convert INDEX (--exact-evidence | --approx-evidence BITS | --hybrid-evidence) [OPTIONS]
|
||||
```
|
||||
|
||||
## Arguments
|
||||
|
||||
| Argument | Description |
|
||||
|---|---|
|
||||
| `INDEX` | Index directory to convert (modified in place) |
|
||||
|
||||
## Options
|
||||
|
||||
Exactly one of the first three is required:
|
||||
|
||||
| Option | Description |
|
||||
|---|---|
|
||||
| `--exact-evidence` | Convert to exact evidence (zero false positives) |
|
||||
| `--approx-evidence BITS` | Convert to approximate (fingerprint-only) evidence; `BITS` = fingerprint bits per slot (b) |
|
||||
| `--hybrid-evidence` | Convert to hybrid evidence (both exact and approximate bundles kept) |
|
||||
| `--evidence-bits BITS` | Fingerprint bits per slot (b) — required with `--hybrid-evidence` when the source index is currently exact; rejected otherwise (the source already fixes `b`) |
|
||||
| `-z, --findere-z Z` | Findere z parameter: number of consecutive stored kmers that must all match to confirm a hit. This does not shorten the indexed kmer length (fixed forever at `index` build time) — it extends the effective match window: on a k=31 index, `z=2` requires 32 consecutive matching bases, not 30 |
|
||||
| `--fp FP` | Target false-positive rate per z-window (e.g. `0.01`); derives `b` or `z` when one of them isn't given directly |
|
||||
| `--block-size N` | Block size for exact evidence's on-disk index (unitigs per block). Ignored when converting to pure approximate evidence. Default `1` |
|
||||
|
||||
See [`index`](index_command.md#exact-vs-approximate-evidence) for the exact/approximate trade-off and the underlying false-positive model, and [`estimate`](estimate.md) to explore parameters beforehand. The index directory is locked for exclusive access during conversion.
|
||||
@@ -15,4 +15,4 @@ obikmer estimate [OPTIONS]
|
||||
| `--evidence-bits` | none | Fingerprint bits per slot (b) |
|
||||
| `--fp` | none | Target false-positive rate per z-window |
|
||||
|
||||
Any two of `-z`, `--evidence-bits`, `--fp` may be given; the third is derived using the same model as `index --approx` and `reindex --approx` ($FP = 1 / 2^{b \cdot z}$). The report printed to stdout includes: query $k$, effective indexed $k$ ($k-z+1$), $z$, evidence bits, per-kmer false-positive rate, and per-z-window false-positive rate.
|
||||
Any two of `-z`, `--evidence-bits`, `--fp` may be given; the third is derived using the same model as `index --approx` and `convert --approx-evidence` ($FP = 1 / 2^{b \cdot z}$). The report printed to stdout includes: query $k$, effective indexed $k$ ($k-z+1$), $z$, evidence bits, per-kmer false-positive rate, and per-z-window false-positive rate.
|
||||
|
||||
@@ -45,6 +45,6 @@ With `--approx`, evidence is stored as a compact **fingerprint** instead, tradin
|
||||
|
||||
$$FP = \frac{1}{2^{b \cdot z}}$$
|
||||
|
||||
where $b$ is `--evidence-bits` and $z$ is `--findere-z`. Any two of `-z`, `--evidence-bits`, `--fp` can be given and the third is derived; if none are given, defaults are $b=8$, $z=1$ ($FP \approx 1/256$). See [`estimate`](estimate.md) to explore this trade-off before building an index, and [`reindex`](reindex.md) to convert an existing index between the two representations.
|
||||
where $b$ is `--evidence-bits` and $z$ is `--findere-z`. Any two of `-z`, `--evidence-bits`, `--fp` can be given and the third is derived; if none are given, defaults are $b=8$, $z=1$ ($FP \approx 1/256$). See [`estimate`](estimate.md) to explore this trade-off before building an index, and [`convert`](convert.md) to change an existing index's representation afterwards.
|
||||
|
||||
`z` must be strictly less than k: the effective indexed kmer length under approximate evidence is k−z+1.
|
||||
|
||||
@@ -1,21 +0,0 @@
|
||||
# name-tree
|
||||
|
||||
Translate a numerically-labelled tree export (TNT, PhyG, or any plain Newick file with bare `1`, `2`, `3`, … leaf labels) back to real taxon names, reading the label order from the FASTA that produced it.
|
||||
|
||||
```bash
|
||||
obikmer name-tree TREE --fasta FASTA -o OUTPUT
|
||||
```
|
||||
|
||||
## Arguments
|
||||
|
||||
| Argument | Description |
|
||||
|---|---|
|
||||
| `TREE` | Tree file to translate — a TNT-style NEXUS export (`tree NAME = [&U] ...;`) or a plain Newick file |
|
||||
| `--fasta` | FASTA file whose record order gives the numeric taxon labels (1-based) — typically the `_sankoff.fasta`/`_snp.fasta` used to produce `TREE` |
|
||||
| `-o, --output` | Output NEXUS file path |
|
||||
|
||||
## Output
|
||||
|
||||
A NEXUS file with a `taxa` block, a `translate` table (numeric label → taxon name, from `--fasta`'s header order), and every tree found in `TREE`, topology unchanged — readable directly in FigTree, PearTree, `ape` (R), etc.
|
||||
|
||||
`--tnt`'s and `--phyg`'s exports (see [phylo](phylo.md)) both number taxa `1..N` in the same order as the pseudo-alignment FASTA they were built from (`<prefix>_sankoff.fasta`), so pass that same file as `--fasta` here.
|
||||
@@ -24,6 +24,6 @@ The index directory is locked for exclusive access while packing.
|
||||
|
||||
Matrix data (which genomes carry each kmer, or with what count) is often mostly empty — most kmers are present in only a handful of genomes out of the whole collection. The default (dense) packed format stores one entry per genome for every kmer regardless of how many genomes actually carry it; `--sparse` instead stores each kmer's genome list directly. For presence/absence matrices, identical genome lists shared by many kmers are also deduplicated (common in real data, since kmers from the same conserved region tend to be carried by the same genomes); for count matrices, the genome list is deduplicated the same way but each kmer's actual counts are kept per-kmer, since two kmers sharing the same genome list rarely carry the same counts.
|
||||
|
||||
On real genome collections this has measured at roughly 7x smaller on disk than the dense format for presence/absence, and single-kmer lookups (the shape `phylo`'s sibling-annex/entropy/Sankoff computations use) are typically faster too, since the smaller files mean less data to read from disk. The trade-off: reading a whole genome column at once (used by `--metric` distance-matrix computations) is much slower on the sparse format than on the dense one, since there is no native column layout to read sequentially — prefer the dense format (the default, no `--sparse`) for indexes you mainly query with `phylo`'s plain `--metric` distance matrices.
|
||||
On real genome collections this has measured at roughly 7x smaller on disk than the dense format for presence/absence, and single-kmer lookups (the shape `phylo`'s sibling-annex/entropy/Sankoff computations use) are typically faster too, since the smaller files mean less data to read from disk. The trade-off: reading a whole genome column at once (used by `--distance` matrix computations) is much slower on the sparse format than on the dense one, since there is no native column layout to read sequentially — prefer the dense format (the default, no `--sparse`) for indexes you mainly query with `phylo`'s `--distance` matrices.
|
||||
|
||||
`--sparse` applies to both presence/absence and count matrices — a count index (`--metric` distance matrices included) is packed sparse the same as a presence index.
|
||||
`--sparse` applies to both presence/absence and count matrices — a count index (`--distance` matrix computations included) is packed sparse the same as a presence index.
|
||||
|
||||
+185
-125
@@ -1,6 +1,6 @@
|
||||
# phylo
|
||||
|
||||
Compute pairwise evolutionary-distance proxies between the genomes stored in an index — a plain distance matrix, optionally trees (NJ/UPGMA), and optionally a central-position SNP model with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).
|
||||
Compute pairwise distances between the genomes stored in an index, optionally build trees (NJ/UPGMA) from them, and optionally calibrate a 16-state parsimony model for a central-position SNP character with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).
|
||||
|
||||
```bash
|
||||
obikmer phylo INDEX [OPTIONS]
|
||||
@@ -12,133 +12,217 @@ obikmer phylo INDEX [OPTIONS]
|
||||
|---|---|
|
||||
| `INDEX` | Index directory |
|
||||
|
||||
## Distance matrix
|
||||
## Distance matrix (`--distance`)
|
||||
|
||||
| Option | Default | Description |
|
||||
|---|---|---|
|
||||
| `--metric` | `jaccard` | One of `jaccard`, `mash`, `hamming`, `bray-curtis`, `relfreq-bray-curtis`, `euclidean`, `relfreq-euclidean`, `hellinger`, `hellinger-euclidean` |
|
||||
| `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for Jaccard/Mash on a count index |
|
||||
| `--shared-kmers` | off | Also write the shared-kmer count matrix |
|
||||
| `--distance` | `jaccard` | See the two tables below for the full list of accepted values |
|
||||
| `--gamma-shape ALPHA` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). No effect on the other values; rejected if given together with a value that doesn't support it |
|
||||
| `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for `jaccard`/`mash` on a count index |
|
||||
| `--csv` | off | Write the matrix as plain CSV instead of the default relaxed-PHYLIP format |
|
||||
| `--shared-kmers` | off | Also write the shared-kmer count matrix. Only valid with a whole-index metric, not a `snp-*` value |
|
||||
| `--nj` | off | Compute and write a Neighbor-Joining tree (Newick) |
|
||||
| `--upgma` | off | Compute and write a UPGMA tree (Newick) |
|
||||
| `-o, --output` | none (stdout) | Output file prefix; without it, the distance matrix is printed to stdout as CSV |
|
||||
| `-o, --output` | none (stdout) | Output file prefix |
|
||||
|
||||
`hamming` requires a presence/absence index. All other metrics work on either index type; on a presence index, `jaccard`/`mash`/`hamming` are the only ones available.
|
||||
Every value routes to one of two independent computations:
|
||||
|
||||
### Metric definitions
|
||||
### Whole-index metrics
|
||||
|
||||
- **jaccard**: $D = 1 - \dfrac{|A \cap B|}{|A \cup B|}$ over the sets of kmers present in each genome.
|
||||
- **mash**: derived from the Jaccard distance via $D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)$ where $J = 1 - D_{\text{jaccard}}$ and $k$ is the index's kmer size; clamped to 1.0 when $J \le 0$.
|
||||
- **hamming**: number of kmer positions where presence differs between the two genomes (presence index only, not normalized): $D = \sum_i \mathbb{1}[a_i \ne b_i]$.
|
||||
- **bray-curtis**: $D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}$ on raw per-kmer counts.
|
||||
- **relfreq-bray-curtis**: the same formula computed on per-genome relative frequencies $p_i = c_i / \sum_j c_j$ instead of raw counts.
|
||||
- **euclidean**: $D = \sqrt{\sum_i (c_i^A - c_i^B)^2}$ on raw counts.
|
||||
- **relfreq-euclidean**: the same formula on relative frequencies.
|
||||
- **hellinger**: $D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}$ on relative frequencies, bounded in $[0, 1]$.
|
||||
- **hellinger-euclidean**: the unnormalized variant, $D = \sqrt{2} \times D_{\text{hellinger}}$.
|
||||
| Value | Definition |
|
||||
|---|---|
|
||||
| `jaccard` | $D = 1 - \dfrac{\lvert A \cap B \rvert}{\lvert A \cup B \rvert}$ over the sets of kmers present in each genome |
|
||||
| `mash` | derived from the Jaccard distance via $D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)$ where $J = 1 - D_{\text{jaccard}}$ and $k$ is the index's kmer size; clamped to 1.0 when $J \le 0$ |
|
||||
| `hamming` | number of kmer positions where presence differs between the two genomes (presence index only, not normalized): $D = \sum_i \mathbb{1}[a_i \ne b_i]$ |
|
||||
| `bray-curtis` | $D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}$ on raw per-kmer counts |
|
||||
| `relfreq-bray-curtis` | the same formula computed on per-genome relative frequencies $p_i = c_i / \sum_j c_j$ instead of raw counts |
|
||||
| `euclidean` | $D = \sqrt{\sum_i (c_i^A - c_i^B)^2}$ on raw counts |
|
||||
| `relfreq-euclidean` | the same formula on relative frequencies |
|
||||
| `hellinger` | $D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}$ on relative frequencies, bounded in $[0, 1]$ |
|
||||
| `hellinger-euclidean` | the unnormalized variant, $D = \sqrt{2} \times D_{\text{hellinger}}$ |
|
||||
|
||||
`hamming` requires a presence/absence index; the others work on either index type.
|
||||
|
||||
### `snp-*` corrections
|
||||
|
||||
Computed from the central-position SNP model (see "Central-position SNP model" below): a family is the set of up to 4 kmers sharing identical flanking sequence and differing only at the central base. These values require the sibling annex (`--sibling-annex`, below) and are, by default, computed exhaustively over every non-monomorphic family in the index; add `--subsample N` to bound the computation to approximately `N` families instead (see "Sampling at scale" below — the same flag `--pseudo-alignment`/`--sankoff` use, but optional here).
|
||||
|
||||
For a genome pair, let $L$ be its total number of eligible loci (both genomes single-copy at that family), $p$ the raw proportion of substitutions among those loci, $P$/$Q$ the transition/transversion proportions, $Q_1$/$Q_2$ Kimura's two transversion categories (A↔C & G↔T vs. A↔T & C↔G), $P_1$/$P_2$ the purine (A↔G) / pyrimidine (C↔T) transition proportions, and $\pi_A,\pi_C,\pi_G,\pi_T$ the pair's pooled base frequencies.
|
||||
|
||||
**`snp-raw`**
|
||||
|
||||
$$d = p$$
|
||||
|
||||
**`snp-jc`**
|
||||
|
||||
$$d = -\frac{3}{4}\ln\!\left(1-\frac{4p}{3}\right)$$
|
||||
|
||||
**`snp-k2p`**
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
a_1 &= 1-2P-Q \\
|
||||
a_2 &= 1-2Q \\
|
||||
d &= -\frac{1}{2}\ln a_1-\frac{1}{4}\ln a_2
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
**`snp-k81`**
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
a_1 &= 1-2P-2Q_1 \\
|
||||
a_2 &= 1-2P-2Q_2 \\
|
||||
a_3 &= 1-2Q_1-2Q_2 \\
|
||||
d &= -\frac{1}{4}\left(\ln a_1+\ln a_2+\ln a_3\right)
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
**`snp-f81`**
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
E &= 1-\left(\pi_A^2+\pi_C^2+\pi_G^2+\pi_T^2\right) \\
|
||||
d &= -E\ln\!\left(1-\frac{p}{E}\right)
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
**`snp-t92`**
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
g &= \pi_C+\pi_G \\
|
||||
w &= 2g(1-g) \\
|
||||
a_1 &= 1-\frac{P}{w}-Q \\
|
||||
a_2 &= 1-2Q \\
|
||||
d &= -w\ln a_1-\frac{1}{2}(1-w)\ln a_2
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
**`snp-tn93`**
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
g_R &= \pi_A+\pi_G \\
|
||||
g_Y &= \pi_C+\pi_T \\
|
||||
k_1 &= \frac{2\pi_A\pi_G}{g_R} \\
|
||||
k_2 &= \frac{2\pi_C\pi_T}{g_Y} \\
|
||||
k_3 &= 2\left(g_Rg_Y-\frac{\pi_A\pi_G\,g_Y}{g_R}-\frac{\pi_C\pi_T\,g_R}{g_Y}\right) \\
|
||||
w_1 &= 1-\frac{P_1}{k_1}-\frac{Q}{2g_R} \\
|
||||
w_2 &= 1-\frac{P_2}{k_2}-\frac{Q}{2g_Y} \\
|
||||
w_3 &= 1-\frac{Q}{2g_Rg_Y} \\
|
||||
d &= -k_1\ln w_1-k_2\ln w_2-k_3\ln w_3
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
**`snp-tv`** — transversions only, deliberately uncorrected:
|
||||
|
||||
$$d = Q$$
|
||||
|
||||
`--gamma-shape ALPHA` applies to every value above except `snp-raw` and `snp-tv`: each $-\ln(x)$ term in the formulas above is replaced by $\alpha\left(x^{-1/\alpha}-1\right)$ (the same weight, same $x$).
|
||||
|
||||
### Output
|
||||
|
||||
Without `-o`, the matrix goes to stdout in relaxed-PHYLIP format (`n` on the first line, then one `label<TAB>value...` row per genome). With `--csv`, the format is instead a header row `genome,<label1>,<label2>,...` followed by one `<label>,<value1>,<value2>,...` row per genome, 6 decimals. Both formats are symmetric with a zero diagonal, except where noted below.
|
||||
|
||||
## `--exclude-genome`, `--min-shared-family`
|
||||
|
||||
| Option | Description |
|
||||
|---|---|
|
||||
| `--exclude-genome LABEL` | Exclude a genome (repeatable). Drops its row/column from the distance/shared-kmer matrix output, and removes it from the sampling used by `--pseudo-alignment`/`--sankoff`/a `snp-*` `--distance` value. Does not change the value computed for any remaining pair |
|
||||
| `--min-shared-family N` | Auto-exclude, on top of `--exclude-genome`, any genome whose mean shared-family count against every other genome (see "Family Overlap" below) falls below `N`. Applies only to `--pseudo-alignment`/`--sankoff`/`snp-*` `--distance` — never to the whole-index metrics or their matrix/NJ/UPGMA output |
|
||||
|
||||
Neighbor-Joining and UPGMA trees (`--nj`/`--upgma`) are always built from every genome in the index, regardless of `--exclude-genome`/`--min-shared-family`.
|
||||
|
||||
## Central-position SNP model
|
||||
|
||||
This is a separate operation from the distance-matrix computation above: if any option below is used, no `--metric` matrix is computed in the same invocation.
|
||||
Requires the sibling annex, built once per index:
|
||||
|
||||
A **family** is the set of up to 4 kmers that share identical flanking sequence and differ only at the exact central base. Because $k$ is odd, the central position is well defined and maps to itself under reverse complementation. All computations below first require building the **sibling annex**, an index-wide record of which of the 4 possible central bases are observed at each family, across every genome.
|
||||
| Option | Description |
|
||||
|---|---|
|
||||
| `--sibling-annex` | Build (or rebuild) the sibling-count/minorant annex — prerequisite for every option in this section, and for a `snp-*` `--distance` value |
|
||||
| `--sibling-stats` | Write `<prefix>_siblings.csv`: the family-size distribution, per genome and globally |
|
||||
| `--sibling-hist` | Print the global family-size histogram (1-4 members) only |
|
||||
| `--family-overlap` | Write `<prefix>_family_overlap.csv`: for every genome pair, how many variable families both genomes carry a call for |
|
||||
| `--shannon` | Write `<prefix>_entropy.csv`: per-family Shannon entropy, one row per family, full unsampled scan |
|
||||
| `--pseudo-alignment` | Write `<prefix>_alignment.fasta`: a SNP-only pseudo-alignment. Requires `--subsample N` |
|
||||
|
||||
A family is eligible for a genome pair $(i,j)$ only if both genomes carry exactly one of its observed forms (single-copy, unambiguous).
|
||||
|
||||
### `--sibling-stats`
|
||||
|
||||
`<prefix>_siblings.csv` — family size = number of distinct central bases observed at a family (1-4).
|
||||
|
||||
| Column | Meaning |
|
||||
|---|---|
|
||||
| `genome` | genome label, or the literal `global` for the last row |
|
||||
| `1`, `2`, `3`, `4` | for a genome row: number of families of that size where the genome carries ≥ 1 member. For the `global` row: the actual deduplicated family-size histogram — not the sum of the rows above |
|
||||
|
||||
### Family Overlap
|
||||
|
||||
`--family-overlap` writes `<prefix>_family_overlap.csv`: header `genome,<label1>,<label2>,...`, one row per genome, cell `[i][j]` = number of variable families (family size ≥ 2) where both genome `i` and genome `j` carry a call. The diagonal is always `0`. Every genome is written, unfiltered by `--exclude-genome`/`--min-shared-family`.
|
||||
|
||||
`--min-shared-family N` uses the mean of each genome's own row (excluding the diagonal) against this same matrix as its exclusion statistic. There is no universal value for `N` — inspect `--family-overlap`'s own output to find where the real gap sits in a given genome collection before choosing a threshold.
|
||||
|
||||
### `--pseudo-alignment`
|
||||
|
||||
`<prefix>_alignment.fasta` — one record per non-excluded genome, one column per variable family (family size ≥ 2). Each site is IUPAC-coded from the genome's presence mask at that family: a single observed form → the plain base; several forms → the matching IUPAC ambiguity code; no form → `-`.
|
||||
|
||||
| Option | Default | Description |
|
||||
|---|---|---|
|
||||
| `--sibling-annex` | off | Build the sibling presence-mask annex (prerequisite for every option below) |
|
||||
| `--exclude-genome LABEL` | none | Exclude a genome (repeatable) from every SNP/Sankoff/export computation below |
|
||||
| `--min-shared-family N` | none | Auto-exclude any genome whose mean shared-family count against every other genome (see `--family-overlap`) falls below `N` — same exclusion as `--exclude-genome`, applied on top of it |
|
||||
| `--sibling-stats` | off | Write the family-size (sibling count) distribution, per genome and globally |
|
||||
| `--raw-snp-distance` | off | Write the single-copy central-SNP p-distance matrix |
|
||||
| `--raw-snp-counts` | off | Write per-pair diagnostic counts (n_snp, n_shared, n_eligible) instead of a matrix |
|
||||
| `--snp` | off | Write a SNP-only pseudo-alignment in FASTA, IUPAC-coded |
|
||||
| `--family-overlap` | off | Write an NxN matrix of, for each genome pair, how many variable families both genomes actually carry a call for; the diagonal holds each genome's own total family count |
|
||||
|
||||
### Locus eligibility
|
||||
|
||||
A family is eligible for a genome pair $(i, j)$ only if genome $i$ carries exactly one of the family's observed forms (single-copy, unambiguous) and genome $j$ also carries exactly one. A genome carrying more than one form at a locus makes that locus ineligible for any pair involving it.
|
||||
|
||||
`--raw-snp-distance` tallies, over every eligible locus of every genome pair, $n_{\text{snp}}$ (the two genomes' single forms differ) versus $n_{\text{shared}}$ (they agree — this includes invariant families). The output ratio is $\hat{p} = \dfrac{n_{\text{snp}}}{n_{\text{snp}} + n_{\text{shared}}}$.
|
||||
|
||||
`--snp` restricts itself to *variable* families (family size $\ge 2$) and writes one FASTA record per genome, one column per family, IUPAC-coded from each genome's presence mask at that family (a single form → the plain base; several forms → the matching IUPAC ambiguity code; no form → `-`).
|
||||
|
||||
`--exclude-genome` removes a genome from these computations, re-checking column variability among the remaining genomes so that a column made monomorphic by the exclusion is dropped rather than kept artificially. It does not affect the `--metric` distance-matrix path.
|
||||
|
||||
### Family overlap and low-coverage genomes
|
||||
|
||||
`--family-overlap` writes, for every genome pair, how many variable families both genomes actually carry a call for (neither is absent) — a direct measure of how much informative content two genomes actually share. On genome-skim or otherwise incomplete-coverage collections, a genome with very little overlap with everything else has almost nothing left to constrain its position in a tree, and tends to end up placed unstably (near-zero branch length, grafted inside an unrelated clade) by `--tnt`/`--iqtree`.
|
||||
|
||||
`--min-shared-family N` automates the fix: it excludes, before any computation, every genome whose mean shared-family count against all other genomes (the same statistic, averaged per row of the `--family-overlap` matrix) falls below `N`. There is no universal value for `N` — it depends on how divergent and how completely covered the genome collection is; inspect `--family-overlap`'s own output to find where the real gap sits before choosing a threshold.
|
||||
|
||||
## Sampling at scale: `--subsample`, `--shannon`, `--entropy`
|
||||
|
||||
On a large index (billions of families), building a full pseudo-alignment or fully calibrating the Sankoff model is not just slow — it may not fit in the time you have. `--subsample` bounds the work to a fixed number of families; `--shannon` reports how informative each family is; `--entropy`/`--entropy-sd` bias which families get kept toward the informative ones instead of choosing uniformly at random.
|
||||
|
||||
| Option | Default | Description |
|
||||
|---|---|---|
|
||||
| `--subsample N` | none (keep everything) | Cap the number of variable families (family size ≥ 2) retained, to approximately `N` |
|
||||
| `--shannon` | off | Write `<prefix>_shannon.csv`: per-family Shannon entropy, one row per family |
|
||||
| `--subsample N` | none (mandatory here) | Target number of families to sample |
|
||||
| `--free-loss` | off | Treat a genome carrying none of a family's observed members as missing data (`?`) instead of `-` |
|
||||
| `--no-ambiguity` | off | Treat a genome carrying more than one member of a family as missing data (`?`) instead of an IUPAC ambiguity code |
|
||||
| `--entropy MU` | off (`1.0` if only `--entropy-sd` is given) | Center of the entropy band to favor when sampling |
|
||||
| `--entropy-sd SIGMA` | off (`0.5` if only `--entropy` is given) | Width of that band |
|
||||
|
||||
`--subsample`/`--entropy`/`--entropy-sd` affect every option that scans variable families: `--snp`, `--family-overlap`, `--shannon`, and the whole Sankoff pipeline (`--sankoff`/`--tnt`/`--phyg`/`--iqtree`, next section) — all of them draw from the *same* selection of families in one invocation, so the Sankoff calibration and the pseudo-alignment it calibrates always describe the same sites, and `--family-overlap`'s counts stay consistent with `--snp`'s columns. `--raw-snp-distance`/`--raw-snp-counts` are not affected — they always scan every family, since their p-distance estimate is a whole-index statistic, not something that benefits from being restricted to a sample.
|
||||
## Sampling at scale: `--subsample`, `--shannon`, `--entropy`
|
||||
|
||||
### `--subsample N`
|
||||
`--subsample`, `--free-loss`, `--no-ambiguity`, `--entropy`/`--entropy-sd` are shared by `--pseudo-alignment`, `--sankoff` (and everything it implies: `--tnt`/`--phyg`/`--iqtree`), and a `snp-*` `--distance` value — one draw feeds all of them in a single invocation. `--subsample` is mandatory for `--pseudo-alignment`/`--sankoff`; for a `snp-*` `--distance` value it is optional (omitted means every non-monomorphic family in the index, not an approximation).
|
||||
|
||||
Without `--subsample`, every variable family (family size ≥ 2, i.e. every family where at least one genome differs from the rest) is used. With `--subsample N`, roughly `N` families are kept instead, chosen at random but in proportion to how many candidate families each part of the index actually holds — so the sample stays representative of the whole index, not skewed toward whichever part happens to be scanned first. If the index has fewer than `N` candidate families in the first place, `--subsample` has no effect: everything is kept.
|
||||
|
||||
`--subsample` trades completeness for speed: `--snp`'s alignment gets fewer columns, `--sankoff`'s calibration is based on fewer observations, but the resolution work (the expensive part of a `phylo` run on a large index) scales with `N` instead of with the index's true size. Pick `N` as large as your time budget allows — a few hundred thousand to a few million families is usually enough for the transition-probability estimates in `--sankoff`'s calibration to stabilize; a smaller `N` speeds up exploratory runs.
|
||||
Without `--subsample`, every variable family (family size ≥ 2) is used. With `--subsample N`, roughly `N` families are kept instead, drawn in proportion to how many candidate families each part of the index actually holds, so the sample stays representative of the whole index. If the index has fewer than `N` candidate families, `--subsample` has no effect.
|
||||
|
||||
### `--shannon`: measuring how informative a family is
|
||||
|
||||
Not every variable family is equally useful for a tree: a family that differs in only one genome out of a thousand carries very little signal, and one where the pattern looks essentially random across genomes may be too saturated (multiple substitutions have overwritten the original signal) to carry real information either. `--shannon` quantifies this with the Shannon entropy (in bits) of each family's states across the genomes that carry it — low entropy means "almost everyone agrees" (an invariant or near-invariant family, phylogenetically shallow), while entropy near the ceiling for a 4-state character means "close to a random draw between the possible bases" (saturated).
|
||||
|
||||
`<prefix>_shannon.csv` (or `shannon.csv` without `-o`) has one row per family visited:
|
||||
`<prefix>_entropy.csv` has one row per family visited:
|
||||
|
||||
| Column | Meaning |
|
||||
|---|---|
|
||||
| `layer` | an internal index-layer identifier — stable within one run, not meaningful across indexes |
|
||||
| `family_idx` | the family's position within that layer |
|
||||
| `entropy15` | Shannon entropy (bits) over the 16 possible states (the 15 non-empty subsets of `{A,C,G,T}` — the same alphabet `--sankoff`'s 16-state model uses), genomes absent from the family excluded from the count |
|
||||
| `entropy4` | Shannon entropy (bits) reduced to the 4 plain bases, kept alongside `entropy15` for comparison — a genome carrying more than one base at once counts once per base, so this can differ from `entropy15` |
|
||||
| `family_size` | number of distinct central bases observed anywhere in the index for this family (2–4, since monomorphic families aren't visited) |
|
||||
| `entropy15` | Shannon entropy (bits) over the 16 possible states (the 15 non-empty subsets of `{A,C,G,T}`), genomes absent from the family excluded from the count |
|
||||
| `entropy4` | Shannon entropy (bits) reduced to the 4 plain bases, kept alongside `entropy15` for comparison |
|
||||
| `family_size` | number of distinct central bases observed anywhere in the index for this family (2-4) |
|
||||
| `n_genomes_present` | how many genomes the entropy was computed over |
|
||||
|
||||
Run with `--subsample N --shannon` to get a bounded diagnostic sample instead of a full-index pass — useful to inspect the entropy distribution and decide reasonable `--entropy`/`--entropy-sd` values (see below) before committing to a full run.
|
||||
Run with `--subsample N --shannon` to get a bounded diagnostic sample instead of a full-index pass — useful for choosing `--entropy`/`--entropy-sd` values before a full run.
|
||||
|
||||
### `--entropy MU` / `--entropy-sd SIGMA`: biasing the sample toward informative families
|
||||
|
||||
By default, `--subsample` draws families uniformly — every candidate family has the same chance of being kept, regardless of how informative it actually is. `--entropy`/`--entropy-sd` change that: instead of a uniform draw, each family's chance of being kept is weighted by how close its own entropy (the `entropy15` value `--shannon` reports) is to `MU`, using a bell-shaped (Gaussian) curve of width `SIGMA` — a family with entropy exactly `MU` is the most likely to be kept, and the chance falls off smoothly the further its entropy is from `MU`, with no hard cutoff (a few families outside the target band can still get in, just less often).
|
||||
By default, `--subsample` draws families uniformly. With `--entropy`/`--entropy-sd`, each family's chance of being kept is instead weighted by how close its own entropy (`entropy15`) is to `MU`, using a Gaussian curve of width `SIGMA` — no hard cutoff. The filter activates as soon as either flag is given; the other defaults to `1.0`/`0.5`. Combine with `--subsample N` (expect somewhat fewer than `N` families kept in practice) or use alone (a soft filter over the whole index, no size target).
|
||||
|
||||
The filter activates as soon as *either* `--entropy` or `--entropy-sd` is given; whichever one you don't set defaults to `1.0`/`0.5`. It can be combined with `--subsample N` (the target count is still approximately `N`, now biased toward the entropy band instead of uniform — expect somewhat *fewer* than `N` families in practice, since low-weight families are dropped rather than replaced) or used alone (`--entropy` without `--subsample`: a soft entropy filter over the whole index, no size target).
|
||||
|
||||
The first `phylo` run on a given index that uses `--entropy`/`--entropy-sd` pays a one-time extra cost (every candidate family's entropy has to be computed once, up front, and is then saved alongside the index). Every following run — even with different `MU`/`SIGMA` values — reuses that saved data and stays fast.
|
||||
The first `phylo` run on a given index that uses `--entropy`/`--entropy-sd` pays a one-time extra cost (every candidate family's entropy is computed once and saved alongside the index); later runs, even with different `MU`/`SIGMA`, reuse that saved data.
|
||||
|
||||
## Sankoff calibration and phylogenetic exports
|
||||
|
||||
| Option | Default | Description |
|
||||
|---|---|---|
|
||||
| `--sankoff` | off | Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment |
|
||||
| `--sankoff-ratio-ceiling` | `0.5` | Exclude genome pairs whose raw SNP ratio exceeds this value from the calibration |
|
||||
| `--free-loss` | off | Recode a family's non-detection as the `?` missing-data symbol instead of an ordinary, costed state, in `--sankoff`'s pseudo-alignment and every export built from it |
|
||||
| `--sankoff` | off | Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment. Requires `--subsample N` |
|
||||
| `--sankoff-ratio-ceiling` | `0.5` | Exclude genome pairs whose raw SNP ratio exceeds this value from the base-composition part of the calibration |
|
||||
| `--free-loss` | off | Recode a family's non-detection as the `?` missing-data symbol instead of an ordinary, costed state, throughout `--sankoff` and every export built from it |
|
||||
| `--tnt` | off | Also write a TNT script (implies `--sankoff`) |
|
||||
| `--phyg` | off | Also write PhyG input files (implies `--sankoff`) |
|
||||
| `--iqtree` | off | Also write an IQ-TREE custom model and alignment (implies `--sankoff`) |
|
||||
| `--iqtree-min-freq` | `0.001` | With `--iqtree --free-loss`: also treat as missing any state rarer than this in the alignment |
|
||||
| `--sankoff-cost-scale` | `100` | Integer scaling factor applied to costs before rounding (required by TNT/PhyG's integer-only cost commands) |
|
||||
| `--sankoff-cost-scale` | `100` | Integer scaling factor applied to costs before rounding, for TNT/PhyG's integer-only cost commands |
|
||||
|
||||
### The 16-state model
|
||||
|
||||
Each family is treated as a character with 16 possible states: one per subset of the 4 possible central bases actually observed (including the empty subset). Calibration combines two tallies, both restricted to genome pairs at or below `--sankoff-ratio-ceiling`:
|
||||
Each family is a character with 16 possible states: one per subset of the 4 possible central bases (including the empty subset). Calibration combines a $5 \times 5$ transition matrix over family cardinality (0-4 observed forms) and a $4 \times 4$ base-substitution matrix from unambiguous single-copy loci, both restricted to genome pairs at or below `--sankoff-ratio-ceiling`, into a row-normalized $16 \times 16$ transition probability matrix $P$, converted to a symmetric cost matrix via $\text{cost}(a,b) = -\ln P(a,b)$.
|
||||
|
||||
- a $5 \times 5$ transition matrix over family cardinality (0–4 observed forms) between paired genomes, and
|
||||
- a $4 \times 4$ base-substitution transition matrix from unambiguous single-copy loci,
|
||||
`--sankoff` alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol `0` (never a gap character). It does not run any external tool.
|
||||
|
||||
which are combined into a row-normalized $16 \times 16$ transition probability matrix $P$, converted to a symmetric cost matrix via $\text{cost}(a,b) = -\ln P(a,b)$.
|
||||
|
||||
`--sankoff` alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol `0` (never a gap character, to avoid ambiguity with external tools' own gap semantics). It does not run any external tool.
|
||||
|
||||
With `--free-loss`, the empty state is recoded to `?` instead — TNT/PhyG/IQ-TREE's own missing-data symbol — rather than an ordinary, costed 16th state. This matters for genome-skim or otherwise incomplete-coverage collections, where non-detection of a family is dominated by sampling failure rather than true evolutionary loss: scoring it as a real state risks grouping genomes by shared undersampling instead of shared ancestry. `?` rather than `-` because `-` still carries gap/indel semantics in these tools, and a non-detected family is not an observed deletion. `--free-loss` also drops the cardinality-transition cost between any two states, not just to/from the empty one: whether a genome shows 1 vs. 2 (etc.) detected members of a family it does carry is exactly as vulnerable to sampling failure as whether the family was detected at all, so gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family. Combine with `--min-shared-family`/`--family-overlap` above: `--free-loss` removes the false signal from non-detection, but a genome left with too little real overlap with everything else will still be placed unstably — excluding it is the other half of the fix.
|
||||
With `--free-loss`, the empty state is recoded to `?` (TNT/PhyG/IQ-TREE's own missing-data symbol) instead of an ordinary, costed 16th state — `?` rather than `-`, since `-` still carries gap/indel semantics in these tools. `--free-loss` also zeroes the cardinality-transition cost between any two states, not just to/from the empty one: gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family.
|
||||
|
||||
### Exports
|
||||
|
||||
@@ -146,58 +230,34 @@ All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment,
|
||||
|
||||
- **`--tnt`**: a self-contained TNT script (alignment recoded to TNT's fixed 16-symbol alphabet, integer-scaled cost matrix re-closed to a metric, a default search block).
|
||||
- **`--phyg`**: a custom cost-matrix file plus a PhyG script reusing the `--sankoff` alignment directly.
|
||||
- **`--iqtree`**: a custom substitution-model file (exchangeability matrix recovered as $R(a,b) = e^{-\text{cost}(a,b)}$, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths (unlike the parsimony step-counts from TNT/PhyG). Only states actually occurring in the alignment are kept and compactly renumbered.
|
||||
- **`--iqtree`**: a custom substitution-model file (exchangeability matrix recovered as $R(a,b) = e^{-\text{cost}(a,b)}$, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths. Only states actually occurring in the alignment are kept and compactly renumbered.
|
||||
|
||||
TNT and PhyG both write trees with bare numeric leaf labels (`1`, `2`, …, in the same order as `<prefix>_sankoff.fasta`). Use [`name-tree`](name-tree.md) on the tool's own tree output plus that same FASTA to get a NEXUS file with real taxon names.
|
||||
TNT and PhyG both write trees with bare numeric leaf labels (`1`, `2`, ..., in the order the genomes appear in `<prefix>_sankoff.fasta`).
|
||||
|
||||
## Output files
|
||||
|
||||
With `-o/--output PREFIX`, the relevant subset of the files below is written. Without `-o`, only the plain `--metric` distance matrix is produced, on stdout. All matrices use genome labels (from the index metadata) as row/column headers, in index order; all CSVs are comma-separated with a header row.
|
||||
With `-o/--output PREFIX`, the relevant subset of the files below is written. Without `-o`, only the distance matrix is produced, on stdout. All matrices use genome labels as row/column headers, in index order.
|
||||
|
||||
### Distance matrix
|
||||
|
||||
| File | Written by | Format | Content |
|
||||
|---|---|---|---|
|
||||
| `<prefix>_dist.csv` | always | CSV matrix | the `--metric` distance, 6 decimals, symmetric, diagonal 0 |
|
||||
| `<prefix>_dist.phy` | always, unless `--csv` | relaxed PHYLIP | the `--distance` matrix |
|
||||
| `<prefix>_dist.csv` | `--csv` | CSV matrix | the `--distance` matrix, 6 decimals |
|
||||
| `<prefix>_shared.csv` | `--shared-kmers` | CSV matrix | shared-kmer count per genome pair (integers) |
|
||||
| `<prefix>_nj.nwk` | `--nj` | Newick | Neighbor-Joining tree, branch lengths from the `--metric` matrix |
|
||||
| `<prefix>_upgma.nwk` | `--upgma` | Newick | UPGMA tree, same matrix |
|
||||
| `<prefix>_nj.nwk` | `--nj` | Newick | Neighbor-Joining tree |
|
||||
| `<prefix>_upgma.nwk` | `--upgma` | Newick | UPGMA tree |
|
||||
|
||||
Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below): header `genome,<label1>,<label2>,...`, one data row per genome, `<label>,<value1>,<value2>,...`.
|
||||
CSV matrix layout (`_dist.csv`, `_shared.csv`, `_family_overlap.csv`): header `genome,<label1>,<label2>,...`, one data row per genome, `<label>,<value1>,<value2>,...`.
|
||||
|
||||
### Central-position SNP model
|
||||
|
||||
| File | Written by | Format | Content |
|
||||
|---|---|---|---|
|
||||
| `<prefix>_siblings.csv` | `--sibling-stats` | CSV table | family-size distribution, per genome and global |
|
||||
| `<prefix>_rawsnp.csv` | `--raw-snp-distance` | CSV matrix | single-copy central-SNP p-distance ($\hat p$), or `NA` |
|
||||
| `<prefix>_rawsnp_counts.csv` | `--raw-snp-counts` | CSV table | per-pair diagnostic counts behind `_rawsnp.csv` |
|
||||
| `<prefix>_snp.fasta` | `--snp` | FASTA | SNP-only pseudo-alignment, IUPAC-coded |
|
||||
| `<prefix>_family_overlap.csv` | `--family-overlap` | CSV matrix | variable families both genomes of a pair carry a call for |
|
||||
| `<prefix>_shannon.csv` | `--shannon` | CSV table | per-family Shannon entropy, see "Sampling at scale" above |
|
||||
|
||||
**`_siblings.csv`** — family size = number of distinct central bases observed at a family (1–4), not "sibling count" (0–3).
|
||||
|
||||
| Column | Meaning |
|
||||
|---|---|
|
||||
| `genome` | genome label, or the literal `global` for the last row |
|
||||
| `1`, `2`, `3`, `4` | for a genome row: number of families of that size where the genome carries ≥ 1 member. For the `global` row: the actual deduplicated family-size histogram — **not** the sum of the rows above (a family shared by several genomes would otherwise be counted once per genome) |
|
||||
|
||||
**`_rawsnp.csv`** — same matrix layout as `_dist.csv`; each cell is $\hat p = n_{\text{snp}}/(n_{\text{snp}}+n_{\text{shared}})$, 6 decimals, or `NA` when the pair has zero eligible loci (distinguishes "identical everywhere eligible" from "nothing eligible at all").
|
||||
|
||||
**`_rawsnp_counts.csv`** — one row per unordered genome pair (not a matrix), the counts `_rawsnp.csv`'s ratio is computed from:
|
||||
|
||||
| Column | Meaning |
|
||||
|---|---|
|
||||
| `genome_a`, `genome_b` | the pair |
|
||||
| `n_snp` | eligible loci where the two genomes' single forms differ |
|
||||
| `n_shared` | eligible loci where they agree (includes invariant families) |
|
||||
| `n_eligible` | `n_snp + n_shared` |
|
||||
| `ratio` | $\hat p$ = `n_snp / n_eligible`, or `NA` if `n_eligible = 0` |
|
||||
|
||||
**`_snp.fasta`** — one record per non-excluded genome, one column per variable family (family size ≥ 2), header carries an `n_sites` annotation. Each site is IUPAC-coded from the genome's presence mask at that family: single observed form → plain base; several forms → matching IUPAC ambiguity code; no form → `-`.
|
||||
|
||||
**`_family_overlap.csv`** — same matrix layout as `_dist.csv`; cell `[i][j]` = number of `_snp.fasta` columns where both genome `i` and `j` carry a call (neither is `-`). Diagonal `[i][i]` is kept (not skipped): it holds genome `i`'s own total variable-family count.
|
||||
| `<prefix>_entropy.csv` | `--shannon` | CSV table | per-family Shannon entropy, see "Sampling at scale" above |
|
||||
| `<prefix>_alignment.fasta` | `--pseudo-alignment` | FASTA | SNP-only pseudo-alignment, IUPAC-coded |
|
||||
|
||||
### Sankoff calibration and exports
|
||||
|
||||
@@ -205,7 +265,7 @@ Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below):
|
||||
|---|---|---|---|
|
||||
| `<prefix>_sankoff_matrix.csv` | `--sankoff`/`--tnt`/`--phyg`/`--iqtree` | CSV matrix | calibrated 16×16 cost matrix |
|
||||
| `<prefix>_sankoff_params.yaml` | same flags | YAML | calibration report (raw tallies + derived probabilities) |
|
||||
| `<prefix>_sankoff.fasta` | same flags | FASTA | Sankoff-recoded pseudo-alignment |
|
||||
| `<prefix>_sankoff.fasta` | same flags | FASTA | Sankoff-recoded pseudo-alignment, header carries an `n_sites` annotation |
|
||||
| `<prefix>_sankoff.tnt` | `--tnt` | TNT script | ready-to-run parsimony search |
|
||||
| `<prefix>_sankoff.tcm` | `--phyg` | PhyG TCM | cost matrix in PhyG's own format |
|
||||
| `<prefix>_sankoff.pg` | `--phyg` | PhyG script | ready-to-run parsimony search |
|
||||
@@ -213,33 +273,33 @@ Matrix layout (`_dist.csv`, `_shared.csv`, and every other "CSV matrix" below):
|
||||
| `<prefix>_iqtree.fasta` | `--iqtree` | FASTA | alignment recoded for that model |
|
||||
| `<prefix>_iqtree_states.csv` | `--iqtree` | CSV table | maps `_iqtree.model`/`_iqtree.fasta`'s compact state symbols back to `_sankoff_matrix.csv`'s alphabet |
|
||||
|
||||
**`_sankoff_matrix.csv`** — header `state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N`: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), `0` standing for the empty/absent state (not `-`, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost $-\ln P(a,b)$, 4 decimals.
|
||||
**`_sankoff_matrix.csv`** — header `state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N`: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), `0` standing for the empty/absent state. One row per source state, one value per destination state, cost $-\ln P(a,b)$, 4 decimals.
|
||||
|
||||
**`_sankoff_params.yaml`** — everything the calibration estimated, structured so it can be reloaded rather than re-parsed:
|
||||
**`_sankoff_params.yaml`**:
|
||||
|
||||
| Key | Meaning |
|
||||
|---|---|
|
||||
| `ratio_ceiling` | the `--sankoff-ratio-ceiling` value used |
|
||||
| `cardinality_transitions` | 5×5 list of `{from, to, count, probability}`, family cardinality (0–4 observed forms) |
|
||||
| `cardinality_transitions` | 5×5 list of `{from, to, count, probability}`, family cardinality (0-4 observed forms) |
|
||||
| `composition_transitions` | 4×4 list of `{from, to, count, probability}`, base letters `A/C/G/T`, single-copy substitutions |
|
||||
|
||||
**`_sankoff.fasta`** — same sites as `_snp.fasta`, recoded to match `_sankoff_matrix.csv`'s alphabet: absent state is `0` (or `?` under `--free-loss`). Excluded genomes dropped; columns left monomorphic by that exclusion are re-checked and dropped too.
|
||||
**`_sankoff.fasta`** — recoded to match `_sankoff_matrix.csv`'s alphabet: absent state is `0` (or `?` under `--free-loss`). Excluded genomes dropped; columns left monomorphic by that exclusion are re-checked and dropped too.
|
||||
|
||||
**`_sankoff.tnt`** (`--tnt`) — self-contained TNT script: `xread` block (alignment recoded to TNT's fixed `0-9A-F` alphabet), an integer-scaled (`--sankoff-cost-scale`) and metric-closed `smatrix`, a default `hold 20; mult; export` search. Run with `printf 'proc <path>;\nquit;\n' | tnt`. Produces `<prefix>_sankoff.tre` (bare numeric leaf labels, order matching `_sankoff.fasta`) — feed both into [`name-tree`](name-tree.md) to recover taxon names.
|
||||
**`_sankoff.tnt`** (`--tnt`) — `xread` block (alignment recoded to TNT's fixed `0-9A-F` alphabet), an integer-scaled (`--sankoff-cost-scale`) and metric-closed `smatrix`, a default `hold 20; mult; export` search. Run with `printf 'proc <path>;\nquit;\n' | tnt`. Produces `<prefix>_sankoff.tre` (bare numeric leaf labels, order matching `_sankoff.fasta`).
|
||||
|
||||
**`_sankoff.tcm`** (`--phyg`) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state `0`, since it's never actually triggered).
|
||||
**`_sankoff.tcm`** (`--phyg`) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state `0`).
|
||||
|
||||
**`_sankoff.pg`** (`--phyg`) — script: `read(prefasta:..., tcm:...)` against `_sankoff.fasta`/`_sankoff.tcm`, a default 300s/4-instance `search`, `report(...)` writing `<prefix>_sankoff.tre` (bare numeric labels, as for `--tnt`). Run with `phyg` from the output directory (the script uses relative file names). Feed the tree plus `_sankoff.fasta` into [`name-tree`](name-tree.md) for taxon names.
|
||||
**`_sankoff.pg`** (`--phyg`) — script: `read(prefasta:..., tcm:...)` against `_sankoff.fasta`/`_sankoff.tcm`, a default 300s/4-instance `search`, `report(...)` writing `<prefix>_sankoff.tre`. Run with `phyg` from the output directory (the script uses relative file names).
|
||||
|
||||
**`_iqtree.model`** (`--iqtree`) — lower-triangular exchangeability matrix $R(a,b) = e^{-\text{cost}(a,b)}$ (one row of increasing length per state, whitespace-separated, PAML order), followed by one line of empirical state frequencies. Only states actually occurring in the alignment are kept, compactly renumbered `0..k-1`.
|
||||
|
||||
**`_iqtree.fasta`** (`--iqtree`) — alignment recoded to that same compact `0..k-1` alphabet (symbols `0-9A-F`). Under `--free-loss`, non-detection becomes `?` and columns left non-informative once missing calls are ignored are dropped first (required for `+ASC`); with `--iqtree-min-freq` also set (the default), any state rarer than that threshold is folded into the same `?` treatment, and non-informative columns are re-checked and dropped again after that. Run with:
|
||||
**`_iqtree.fasta`** (`--iqtree`) — alignment recoded to that same compact `0..k-1` alphabet (symbols `0-9A-F`). Under `--free-loss`, non-detection becomes `?` and columns left non-informative once missing calls are ignored are dropped first (required for `+ASC`); with `--iqtree-min-freq` also set (the default), any state rarer than that threshold is folded into the same `?` treatment, and non-informative columns are re-checked and dropped again. Run with:
|
||||
```
|
||||
iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
|
||||
```
|
||||
|
||||
**`_iqtree_states.csv`** (`--iqtree`) — one row per state actually kept in `_iqtree.model`/`_iqtree.fasta` (header `iqtree_symbol,canonical_symbol,frequency`): `iqtree_symbol` is the compact `0-9A-F` symbol as written in those two files, `canonical_symbol` is the matching `_sankoff_matrix.csv` state, `frequency` is that state's empirical frequency at full precision (`_iqtree.model`'s own frequency line is rounded to 6 decimals). Under `--free-loss`, absent (`0`/`?`) is never a kept state, so it never appears here — nor does any state `--iqtree-min-freq` folded away for being too rare. Use this file to identify which real state a given row/column of `_iqtree.model`'s matrix corresponds to — e.g. to check whether a state showing zero exchangeability with everything else is expected (a state combination that never co-occurs with anything else in this data) or worth investigating further.
|
||||
**`_iqtree_states.csv`** (`--iqtree`) — one row per state actually kept in `_iqtree.model`/`_iqtree.fasta` (header `iqtree_symbol,canonical_symbol,frequency`): `iqtree_symbol` is the compact `0-9A-F` symbol as written in those two files, `canonical_symbol` is the matching `_sankoff_matrix.csv` state, `frequency` is that state's empirical frequency at full precision. Under `--free-loss`, absent (`0`/`?`) is never a kept state, so it never appears here — nor does any state `--iqtree-min-freq` folded away for being too rare.
|
||||
|
||||
### Rare states and `--iqtree-min-freq`
|
||||
|
||||
States that combine 3 or 4 central bases at once (IUPAC `V`/`H`/`K`.../`N`) are inherently rare — and, on real data, rare enough that they can make `iqtree3` itself numerically unstable ("Numerical underflow for lh-derivative" warnings, near-degenerate likelihood optimization). They're also more likely to be assembly/detection noise than genuine, widely-shared multi-way polymorphism, the same "sampling failure, not true signal" reasoning `--free-loss` already applies to non-detection. With `--free-loss` set, `--iqtree-min-freq` (default `0.001`, i.e. one in a thousand) extends that same missing-data treatment to any state below this frequency, not just absence. Check `_iqtree_states.csv` to see exactly which states survived and at what frequency; set `--iqtree-min-freq 0` to disable this and keep every state that occurs at all (the old behavior). Has no effect without `--free-loss` — there is no missing-data symbol to fold rare states into otherwise.
|
||||
States that combine 3 or 4 central bases at once (IUPAC `V`/`H`/`K`.../`N`) are inherently rare, and can make `iqtree3` itself numerically unstable ("Numerical underflow for lh-derivative" warnings). With `--free-loss` set, `--iqtree-min-freq` (default `0.001`, one in a thousand) extends the missing-data treatment to any state below this frequency, not just absence. Check `_iqtree_states.csv` to see exactly which states survived and at what frequency; set `--iqtree-min-freq 0` to keep every state that occurs at all. Has no effect without `--free-loss`.
|
||||
|
||||
@@ -1,25 +0,0 @@
|
||||
# reindex
|
||||
|
||||
Convert an existing index's evidence representation in place, between exact and approximate.
|
||||
|
||||
```bash
|
||||
obikmer reindex INDEX [OPTIONS]
|
||||
```
|
||||
|
||||
## Arguments
|
||||
|
||||
| Argument | Description |
|
||||
|---|---|
|
||||
| `INDEX` | Index directory to convert (modified in place) |
|
||||
|
||||
## Options
|
||||
|
||||
| Option | Default | Description |
|
||||
|---|---|---|
|
||||
| `--approx` | off | Convert to approximate evidence (default direction is approximate → exact); requires `-z`/`--evidence-bits`/`--fp` |
|
||||
| `-z, --findere-z` | none | Findere z parameter (≥ 1) |
|
||||
| `--evidence-bits` | none | Fingerprint bits per slot (b) |
|
||||
| `--fp` | none | Target false-positive rate per z-window |
|
||||
| `--block-size` | `1` | Block size for the exact on-disk index (ignored when converting to approximate) |
|
||||
|
||||
See [`index`](index_command.md#exact-vs-approximate-evidence) for the exact/approximate trade-off and the underlying false-positive model, and [`estimate`](estimate.md) to explore parameters beforehand. The index directory is locked for exclusive access during conversion.
|
||||
+2
-30
@@ -693,34 +693,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="/usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="/usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -778,14 +750,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="/usage/reindex/" class="md-nav__link">
|
||||
<a href="/usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -700,34 +700,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -785,14 +757,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../usage/reindex/" class="md-nav__link">
|
||||
<a href="../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1201,7 +1173,7 @@
|
||||
<li>Genome columns (count or presence data) are kept at a consistent width across every layer and partition after a merge, which is what allows whole-index aggregate distances (Jaccard, Bray-Curtis, Euclidean, Hellinger, …) to be computed as a two-pass cascade (local partial sums per partition, then a global combination) with no double counting.</li>
|
||||
</ul>
|
||||
<h2 id="parallel-execution-and-numa-awareness">Parallel execution and NUMA awareness</h2>
|
||||
<p>Partition-level work (index construction, <code>merge</code>, <code>filter</code>, <code>reindex</code>, <code>select</code>, <code>phylo</code>'s sibling-annex/Sankoff computations) is dispatched by a partition runner that adapts to the machine's memory topology, detected automatically at startup via hwloc:</p>
|
||||
<p>Partition-level work (index construction, <code>merge</code>, <code>filter</code>, <code>convert</code>, <code>select</code>, <code>phylo</code>'s sibling-annex/Sankoff computations) is dispatched by a partition runner that adapts to the machine's memory topology, detected automatically at startup via hwloc:</p>
|
||||
<ul>
|
||||
<li>On a multi-socket / multi-NUMA-node machine, one thread pool is pinned per NUMA node, and each partition is processed entirely by threads pinned to one node — keeping the memory a partition touches local to that node's DRAM. This matters because touching kmer data across NUMA nodes without pinning can degrade throughput by an order of magnitude or more on large multi-socket machines.</li>
|
||||
<li>On a single-socket machine, Apple Silicon, or if hwloc cannot report NUMA topology, all cores are treated as one node with no pinning and negligible overhead — this is the default behavior on macOS.</li>
|
||||
|
||||
@@ -702,34 +702,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -787,14 +759,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/reindex/" class="md-nav__link">
|
||||
<a href="../../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1203,7 +1175,7 @@
|
||||
<h2 id="minimal-perfect-hash-function-mphf">Minimal perfect hash function (MPHF)</h2>
|
||||
<p>Each partition's surviving kmers are mapped to a dense range of integer slots by a minimal perfect hash function: no collisions, near-optimal space (a few bits per key), O(1) lookup. Because an MPHF maps <em>any</em> input to some slot — including kmers that were never indexed — a lookup alone cannot distinguish a genuinely indexed kmer from an arbitrary one; every lookup is followed by an evidence check.</p>
|
||||
<h2 id="evidence-exact-vs-approximate">Evidence: exact vs. approximate</h2>
|
||||
<p>Two verification modes are available, selected at build time (<code>index --approx</code>) and convertible afterwards (<a href="../../usage/reindex/"><code>reindex</code></a>):</p>
|
||||
<p>Two verification modes are available, selected at build time (<code>index --approx</code>) and convertible afterwards (<a href="../../usage/convert/"><code>convert</code></a>):</p>
|
||||
<ul>
|
||||
<li><strong>Exact</strong> (default): the hashed slot stores a pointer back into the partition's unitig data. At query time the kmer is reconstructed from that location and compared directly to the query. Zero false positives, at the cost of one extra random read per lookup.</li>
|
||||
<li><strong>Approximate</strong> (<code>--approx</code>): the slot stores a short fingerprint (<code>--evidence-bits</code> bits) instead of a pointer; verification is a single fingerprint comparison. This trades a small, bounded false-positive rate (<span class="arithmatex">\(1/2^b\)</span> per kmer, reduced further to about <span class="arithmatex">\(1/2^{b \cdot z}\)</span> for a read requiring <span class="arithmatex">\(z\)</span> consecutive matching kmers via the <code>-z</code>/<code>--findere-z</code> parameter) for lower memory and disk usage, since no reconstruction index is needed. See <a href="../../usage/estimate/"><code>estimate</code></a> to explore this trade-off before building.</li>
|
||||
@@ -1230,7 +1202,7 @@
|
||||
</code></pre></div>
|
||||
<p><code>unitigs.bin</code> is the only file from which the indexed kmer content can be fully recovered; it is always retained. Every other file (MPHF, evidence, counts) is derived from it.</p>
|
||||
<p>A <strong>layer</strong> corresponds to one increment of kmer content added to a partition — most commonly, one <a href="../../usage/merge/"><code>merge</code></a> operation that introduces kmers not already present in the index. Genomes already present in the index simply gain new columns in the existing layers' count/presence data; only genuinely new kmer content is assembled into a new layer. Because of this, merging cost scales with the novel kmer content being added, not with the accumulated size of the index. A query against an index with several layers checks each layer's MPHF in turn.</p>
|
||||
<p>Sources merged together must share the same kmer size, minimizer size, partition count, and evidence mode (including matching approximate-mode parameters); mismatches are rejected rather than silently reconciled — <a href="../../usage/reindex/"><code>reindex</code></a> one of the sources first if needed.</p>
|
||||
<p>Sources merged together must share the same kmer size, minimizer size, partition count, and evidence mode (including matching approximate-mode parameters); mismatches are rejected rather than silently reconciled — <a href="../../usage/convert/"><code>convert</code></a> one of the sources first if needed.</p>
|
||||
<p><code>obikmer pack</code> consolidates a partition's per-column files (counts/presence) into a single file, reducing the number of file opens needed at query time.</p>
|
||||
|
||||
|
||||
|
||||
+5
-37
@@ -813,34 +813,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -898,14 +870,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="usage/reindex/" class="md-nav__link">
|
||||
<a href="usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1282,11 +1254,7 @@
|
||||
</tr>
|
||||
<tr>
|
||||
<td><a href="usage/phylo/"><code>phylo</code></a></td>
|
||||
<td>Compute pairwise evolutionary-distance proxies, trees, and phylogenetic exports</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><a href="usage/name-tree/"><code>name-tree</code></a></td>
|
||||
<td>Translate a TNT/PhyG numeric-label tree export back to real taxon names</td>
|
||||
<td>Compute pairwise genome distances, trees, and phylogenetic exports</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><a href="usage/unitig/"><code>unitig</code></a></td>
|
||||
@@ -1297,8 +1265,8 @@
|
||||
<td>Estimate approximate-index parameters before indexing</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><a href="usage/reindex/"><code>reindex</code></a></td>
|
||||
<td>Convert an index's evidence representation (exact ↔ approximate)</td>
|
||||
<td><a href="usage/convert/"><code>convert</code></a></td>
|
||||
<td>Convert an index's evidence representation (exact/approximate/hybrid), in place</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><a href="usage/utils/"><code>utils</code></a></td>
|
||||
|
||||
@@ -888,34 +888,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -973,14 +945,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../usage/reindex/" class="md-nav__link">
|
||||
<a href="../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
Binary file not shown.
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/reindex/" class="md-nav__link">
|
||||
<a href="../../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -806,34 +806,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -891,14 +863,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/reindex/" class="md-nav__link">
|
||||
<a href="../../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/reindex/" class="md-nav__link">
|
||||
<a href="../../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -790,34 +790,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -875,14 +847,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/reindex/" class="md-nav__link">
|
||||
<a href="../../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -801,34 +801,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/unitig/" class="md-nav__link">
|
||||
|
||||
@@ -886,14 +858,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../../usage/reindex/" class="md-nav__link">
|
||||
<a href="../../usage/convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -773,34 +773,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -858,14 +830,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -23,7 +23,7 @@
|
||||
|
||||
|
||||
|
||||
<title>reindex - obikmer — User Guide</title>
|
||||
<title>convert - obikmer — User Guide</title>
|
||||
|
||||
|
||||
|
||||
@@ -64,7 +64,7 @@
|
||||
<div data-md-component="skip">
|
||||
|
||||
|
||||
<a href="#reindex" class="md-skip">
|
||||
<a href="#convert" class="md-skip">
|
||||
Skip to content
|
||||
</a>
|
||||
|
||||
@@ -100,7 +100,7 @@
|
||||
<div class="md-header__topic" data-md-component="header-topic">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
</span>
|
||||
</div>
|
||||
@@ -704,34 +704,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -805,7 +777,7 @@
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -823,7 +795,7 @@
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1144,9 +1116,9 @@
|
||||
|
||||
|
||||
|
||||
<h1 id="reindex">reindex</h1>
|
||||
<p>Convert an existing index's evidence representation in place, between exact and approximate.</p>
|
||||
<div class="highlight"><pre><span></span><code>obikmer<span class="w"> </span>reindex<span class="w"> </span>INDEX<span class="w"> </span><span class="o">[</span>OPTIONS<span class="o">]</span>
|
||||
<h1 id="convert">convert</h1>
|
||||
<p>Convert an existing index's evidence representation in place, between exact, approximate, and hybrid.</p>
|
||||
<div class="highlight"><pre><span></span><code>obikmer<span class="w"> </span>convert<span class="w"> </span>INDEX<span class="w"> </span><span class="o">(</span>--exact-evidence<span class="w"> </span><span class="p">|</span><span class="w"> </span>--approx-evidence<span class="w"> </span>BITS<span class="w"> </span><span class="p">|</span><span class="w"> </span>--hybrid-evidence<span class="o">)</span><span class="w"> </span><span class="o">[</span>OPTIONS<span class="o">]</span>
|
||||
</code></pre></div>
|
||||
<h2 id="arguments">Arguments</h2>
|
||||
<table>
|
||||
@@ -1164,39 +1136,42 @@
|
||||
</tbody>
|
||||
</table>
|
||||
<h2 id="options">Options</h2>
|
||||
<p>Exactly one of the first three is required:</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Option</th>
|
||||
<th>Default</th>
|
||||
<th>Description</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>--approx</code></td>
|
||||
<td>off</td>
|
||||
<td>Convert to approximate evidence (default direction is approximate → exact); requires <code>-z</code>/<code>--evidence-bits</code>/<code>--fp</code></td>
|
||||
<td><code>--exact-evidence</code></td>
|
||||
<td>Convert to exact evidence (zero false positives)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>-z, --findere-z</code></td>
|
||||
<td>none</td>
|
||||
<td>Findere z parameter (≥ 1)</td>
|
||||
<td><code>--approx-evidence BITS</code></td>
|
||||
<td>Convert to approximate (fingerprint-only) evidence; <code>BITS</code> = fingerprint bits per slot (b)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--evidence-bits</code></td>
|
||||
<td>none</td>
|
||||
<td>Fingerprint bits per slot (b)</td>
|
||||
<td><code>--hybrid-evidence</code></td>
|
||||
<td>Convert to hybrid evidence (both exact and approximate bundles kept)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--fp</code></td>
|
||||
<td>none</td>
|
||||
<td>Target false-positive rate per z-window</td>
|
||||
<td><code>--evidence-bits BITS</code></td>
|
||||
<td>Fingerprint bits per slot (b) — required with <code>--hybrid-evidence</code> when the source index is currently exact; rejected otherwise (the source already fixes <code>b</code>)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--block-size</code></td>
|
||||
<td><code>1</code></td>
|
||||
<td>Block size for the exact on-disk index (ignored when converting to approximate)</td>
|
||||
<td><code>-z, --findere-z Z</code></td>
|
||||
<td>Findere z parameter: number of consecutive stored kmers that must all match to confirm a hit. This does not shorten the indexed kmer length (fixed forever at <code>index</code> build time) — it extends the effective match window: on a k=31 index, <code>z=2</code> requires 32 consecutive matching bases, not 30</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--fp FP</code></td>
|
||||
<td>Target false-positive rate per z-window (e.g. <code>0.01</code>); derives <code>b</code> or <code>z</code> when one of them isn't given directly</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--block-size N</code></td>
|
||||
<td>Block size for exact evidence's on-disk index (unitigs per block). Ignored when converting to pure approximate evidence. Default <code>1</code></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
@@ -773,34 +773,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -858,14 +830,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -12,7 +12,7 @@
|
||||
<link rel="prev" href="../unitig/">
|
||||
|
||||
|
||||
<link rel="next" href="../reindex/">
|
||||
<link rel="next" href="../convert/">
|
||||
|
||||
|
||||
|
||||
@@ -704,34 +704,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -847,14 +819,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1158,7 +1130,7 @@
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p>Any two of <code>-z</code>, <code>--evidence-bits</code>, <code>--fp</code> may be given; the third is derived using the same model as <code>index --approx</code> and <code>reindex --approx</code> (<span class="arithmatex">\(FP = 1 / 2^{b \cdot z}\)</span>). The report printed to stdout includes: query <span class="arithmatex">\(k\)</span>, effective indexed <span class="arithmatex">\(k\)</span> (<span class="arithmatex">\(k-z+1\)</span>), <span class="arithmatex">\(z\)</span>, evidence bits, per-kmer false-positive rate, and per-z-window false-positive rate.</p>
|
||||
<p>Any two of <code>-z</code>, <code>--evidence-bits</code>, <code>--fp</code> may be given; the third is derived using the same model as <code>index --approx</code> and <code>convert --approx-evidence</code> (<span class="arithmatex">\(FP = 1 / 2^{b \cdot z}\)</span>). The report printed to stdout includes: query <span class="arithmatex">\(k\)</span>, effective indexed <span class="arithmatex">\(k\)</span> (<span class="arithmatex">\(k-z+1\)</span>), <span class="arithmatex">\(z\)</span>, evidence bits, per-kmer false-positive rate, and per-z-window false-positive rate.</p>
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1301,7 +1273,7 @@
|
||||
<p>By default, an index stores <strong>exact</strong> evidence: a kmer is either present or absent (or has an exact count with <code>--with-counts</code>), with no false positives.</p>
|
||||
<p>With <code>--approx</code>, evidence is stored as a compact <strong>fingerprint</strong> instead, trading a small, tunable false-positive rate for reduced memory/disk usage. The false-positive model is:</p>
|
||||
<div class="arithmatex">\[FP = \frac{1}{2^{b \cdot z}}\]</div>
|
||||
<p>where <span class="arithmatex">\(b\)</span> is <code>--evidence-bits</code> and <span class="arithmatex">\(z\)</span> is <code>--findere-z</code>. Any two of <code>-z</code>, <code>--evidence-bits</code>, <code>--fp</code> can be given and the third is derived; if none are given, defaults are <span class="arithmatex">\(b=8\)</span>, <span class="arithmatex">\(z=1\)</span> (<span class="arithmatex">\(FP \approx 1/256\)</span>). See <a href="../estimate/"><code>estimate</code></a> to explore this trade-off before building an index, and <a href="../reindex/"><code>reindex</code></a> to convert an existing index between the two representations.</p>
|
||||
<p>where <span class="arithmatex">\(b\)</span> is <code>--evidence-bits</code> and <span class="arithmatex">\(z\)</span> is <code>--findere-z</code>. Any two of <code>-z</code>, <code>--evidence-bits</code>, <code>--fp</code> can be given and the third is derived; if none are given, defaults are <span class="arithmatex">\(b=8\)</span>, <span class="arithmatex">\(z=1\)</span> (<span class="arithmatex">\(FP \approx 1/256\)</span>). See <a href="../estimate/"><code>estimate</code></a> to explore this trade-off before building an index, and <a href="../convert/"><code>convert</code></a> to change an existing index's representation afterwards.</p>
|
||||
<p><code>z</code> must be strictly less than k: the effective indexed kmer length under approximate evidence is k−z+1.</p>
|
||||
|
||||
|
||||
|
||||
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -704,34 +704,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -789,14 +761,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1198,15 +1170,15 @@
|
||||
<tr>
|
||||
<td><code>--sparse</code></td>
|
||||
<td>off</td>
|
||||
<td>Pack presence/absence matrices into a sparse, deduplicated format instead of the dense one</td>
|
||||
<td>Pack presence/absence and count matrices into a sparse, deduplicated format instead of the dense one</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p>The index directory is locked for exclusive access while packing.</p>
|
||||
<h2 id="-sparse"><code>--sparse</code></h2>
|
||||
<p>Presence/absence data (which genomes carry each kmer) is often mostly empty — most kmers are present in only a handful of genomes out of the whole collection. The default (dense) packed format stores one bit per genome for every kmer regardless of how many genomes actually carry it; <code>--sparse</code> instead stores each kmer's genome list directly, and deduplicates identical lists shared by many kmers (common in real data, since kmers from the same conserved region tend to be carried by the same genomes).</p>
|
||||
<p>On real genome collections this has measured at roughly 7x smaller on disk than the dense format, and single-kmer lookups (the shape <code>phylo</code>'s sibling-annex/entropy/Sankoff computations use) are typically faster too, since the smaller files mean less data to read from disk. The trade-off: reading a whole genome column at once (used by <code>--metric</code> distance-matrix computations) is much slower on the sparse format than on the dense one, since there is no native column layout to read sequentially — prefer the dense format (the default, no <code>--sparse</code>) for indexes you mainly query with <code>phylo</code>'s plain <code>--metric</code> distance matrices.</p>
|
||||
<p>Count matrices (<code>--metric</code> on a count index) are not affected by <code>--sparse</code> — only presence/absence matrices are.</p>
|
||||
<p>Matrix data (which genomes carry each kmer, or with what count) is often mostly empty — most kmers are present in only a handful of genomes out of the whole collection. The default (dense) packed format stores one entry per genome for every kmer regardless of how many genomes actually carry it; <code>--sparse</code> instead stores each kmer's genome list directly. For presence/absence matrices, identical genome lists shared by many kmers are also deduplicated (common in real data, since kmers from the same conserved region tend to be carried by the same genomes); for count matrices, the genome list is deduplicated the same way but each kmer's actual counts are kept per-kmer, since two kmers sharing the same genome list rarely carry the same counts.</p>
|
||||
<p>On real genome collections this has measured at roughly 7x smaller on disk than the dense format for presence/absence, and single-kmer lookups (the shape <code>phylo</code>'s sibling-annex/entropy/Sankoff computations use) are typically faster too, since the smaller files mean less data to read from disk. The trade-off: reading a whole genome column at once (used by <code>--distance</code> matrix computations) is much slower on the sparse format than on the dense one, since there is no native column layout to read sequentially — prefer the dense format (the default, no <code>--sparse</code>) for indexes you mainly query with <code>phylo</code>'s <code>--distance</code> matrices.</p>
|
||||
<p><code>--sparse</code> applies to both presence/absence and count matrices — a count index (<code>--distance</code> matrix computations included) is packed sparse the same as a presence index.</p>
|
||||
|
||||
|
||||
|
||||
|
||||
+354
-254
@@ -12,7 +12,7 @@
|
||||
<link rel="prev" href="../annotate/">
|
||||
|
||||
|
||||
<link rel="next" href="../name-tree/">
|
||||
<link rel="next" href="../unitig/">
|
||||
|
||||
|
||||
|
||||
@@ -748,22 +748,44 @@
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#distance-matrix" class="md-nav__link">
|
||||
<a href="#distance-matrix-distance" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Distance matrix
|
||||
Distance matrix (--distance)
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
<nav class="md-nav" aria-label="Distance matrix">
|
||||
<nav class="md-nav" aria-label="Distance matrix (--distance)">
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#metric-definitions" class="md-nav__link">
|
||||
<a href="#whole-index-metrics" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Metric definitions
|
||||
Whole-index metrics
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#snp-corrections" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
snp-* corrections
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#output" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Output
|
||||
|
||||
</span>
|
||||
</a>
|
||||
@@ -773,6 +795,17 @@
|
||||
</ul>
|
||||
</nav>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-exclude-genome-min-shared-family" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
--exclude-genome, --min-shared-family
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
@@ -788,10 +821,10 @@
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#locus-eligibility" class="md-nav__link">
|
||||
<a href="#-sibling-stats" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Locus eligibility
|
||||
--sibling-stats
|
||||
|
||||
</span>
|
||||
</a>
|
||||
@@ -799,10 +832,21 @@
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#family-overlap-and-low-coverage-genomes" class="md-nav__link">
|
||||
<a href="#family-overlap" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Family overlap and low-coverage genomes
|
||||
Family Overlap
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-pseudo-alignment" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
--pseudo-alignment
|
||||
|
||||
</span>
|
||||
</a>
|
||||
@@ -826,17 +870,6 @@
|
||||
<nav class="md-nav" aria-label="Sampling at scale: --subsample, --shannon, --entropy">
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-subsample-n" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
--subsample N
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-shannon-measuring-how-informative-a-family-is" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
@@ -916,7 +949,7 @@
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#distance-matrix_1" class="md-nav__link">
|
||||
<a href="#distance-matrix" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Distance matrix
|
||||
@@ -979,34 +1012,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -1064,14 +1069,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
@@ -1323,22 +1328,44 @@
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#distance-matrix" class="md-nav__link">
|
||||
<a href="#distance-matrix-distance" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Distance matrix
|
||||
Distance matrix (--distance)
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
<nav class="md-nav" aria-label="Distance matrix">
|
||||
<nav class="md-nav" aria-label="Distance matrix (--distance)">
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#metric-definitions" class="md-nav__link">
|
||||
<a href="#whole-index-metrics" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Metric definitions
|
||||
Whole-index metrics
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#snp-corrections" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
snp-* corrections
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#output" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Output
|
||||
|
||||
</span>
|
||||
</a>
|
||||
@@ -1348,6 +1375,17 @@
|
||||
</ul>
|
||||
</nav>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-exclude-genome-min-shared-family" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
--exclude-genome, --min-shared-family
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
@@ -1363,10 +1401,10 @@
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#locus-eligibility" class="md-nav__link">
|
||||
<a href="#-sibling-stats" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Locus eligibility
|
||||
--sibling-stats
|
||||
|
||||
</span>
|
||||
</a>
|
||||
@@ -1374,10 +1412,21 @@
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#family-overlap-and-low-coverage-genomes" class="md-nav__link">
|
||||
<a href="#family-overlap" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Family overlap and low-coverage genomes
|
||||
Family Overlap
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-pseudo-alignment" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
--pseudo-alignment
|
||||
|
||||
</span>
|
||||
</a>
|
||||
@@ -1401,17 +1450,6 @@
|
||||
<nav class="md-nav" aria-label="Sampling at scale: --subsample, --shannon, --entropy">
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-subsample-n" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
--subsample N
|
||||
|
||||
</span>
|
||||
</a>
|
||||
|
||||
</li>
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#-shannon-measuring-how-informative-a-family-is" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
@@ -1491,7 +1529,7 @@
|
||||
<ul class="md-nav__list">
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="#distance-matrix_1" class="md-nav__link">
|
||||
<a href="#distance-matrix" class="md-nav__link">
|
||||
<span class="md-ellipsis">
|
||||
|
||||
Distance matrix
|
||||
@@ -1557,7 +1595,7 @@
|
||||
|
||||
|
||||
<h1 id="phylo">phylo</h1>
|
||||
<p>Compute pairwise evolutionary-distance proxies between the genomes stored in an index — a plain distance matrix, optionally trees (NJ/UPGMA), and optionally a central-position SNP model with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).</p>
|
||||
<p>Compute pairwise distances between the genomes stored in an index, optionally build trees (NJ/UPGMA) from them, and optionally calibrate a 16-state parsimony model for a central-position SNP character with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).</p>
|
||||
<div class="highlight"><pre><span></span><code>obikmer<span class="w"> </span>phylo<span class="w"> </span>INDEX<span class="w"> </span><span class="o">[</span>OPTIONS<span class="o">]</span>
|
||||
</code></pre></div>
|
||||
<h2 id="arguments">Arguments</h2>
|
||||
@@ -1575,7 +1613,7 @@
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<h2 id="distance-matrix">Distance matrix</h2>
|
||||
<h2 id="distance-matrix-distance">Distance matrix (<code>--distance</code>)</h2>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
@@ -1586,19 +1624,29 @@
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>--metric</code></td>
|
||||
<td><code>--distance</code></td>
|
||||
<td><code>jaccard</code></td>
|
||||
<td>One of <code>jaccard</code>, <code>mash</code>, <code>hamming</code>, <code>bray-curtis</code>, <code>relfreq-bray-curtis</code>, <code>euclidean</code>, <code>relfreq-euclidean</code>, <code>hellinger</code>, <code>hellinger-euclidean</code></td>
|
||||
<td>See the two tables below for the full list of accepted values</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--gamma-shape ALPHA</code></td>
|
||||
<td>none</td>
|
||||
<td>Rate-heterogeneity correction, for <code>snp-*</code> values that support it (see below). No effect on the other values; rejected if given together with a value that doesn't support it</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--presence-threshold</code></td>
|
||||
<td><code>1</code></td>
|
||||
<td>Minimum count for a kmer to be considered present, for Jaccard/Mash on a count index</td>
|
||||
<td>Minimum count for a kmer to be considered present, for <code>jaccard</code>/<code>mash</code> on a count index</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--csv</code></td>
|
||||
<td>off</td>
|
||||
<td>Write the matrix as plain CSV instead of the default relaxed-PHYLIP format</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--shared-kmers</code></td>
|
||||
<td>off</td>
|
||||
<td>Also write the shared-kmer count matrix</td>
|
||||
<td>Also write the shared-kmer count matrix. Only valid with a whole-index metric, not a <code>snp-*</code> value</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--nj</code></td>
|
||||
@@ -1613,87 +1661,201 @@
|
||||
<tr>
|
||||
<td><code>-o, --output</code></td>
|
||||
<td>none (stdout)</td>
|
||||
<td>Output file prefix; without it, the distance matrix is printed to stdout as CSV</td>
|
||||
<td>Output file prefix</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><code>hamming</code> requires a presence/absence index. All other metrics work on either index type; on a presence index, <code>jaccard</code>/<code>mash</code>/<code>hamming</code> are the only ones available.</p>
|
||||
<h3 id="metric-definitions">Metric definitions</h3>
|
||||
<ul>
|
||||
<li><strong>jaccard</strong>: <span class="arithmatex">\(D = 1 - \dfrac{|A \cap B|}{|A \cup B|}\)</span> over the sets of kmers present in each genome.</li>
|
||||
<li><strong>mash</strong>: derived from the Jaccard distance via <span class="arithmatex">\(D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)\)</span> where <span class="arithmatex">\(J = 1 - D_{\text{jaccard}}\)</span> and <span class="arithmatex">\(k\)</span> is the index's kmer size; clamped to 1.0 when <span class="arithmatex">\(J \le 0\)</span>.</li>
|
||||
<li><strong>hamming</strong>: number of kmer positions where presence differs between the two genomes (presence index only, not normalized): <span class="arithmatex">\(D = \sum_i \mathbb{1}[a_i \ne b_i]\)</span>.</li>
|
||||
<li><strong>bray-curtis</strong>: <span class="arithmatex">\(D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}\)</span> on raw per-kmer counts.</li>
|
||||
<li><strong>relfreq-bray-curtis</strong>: the same formula computed on per-genome relative frequencies <span class="arithmatex">\(p_i = c_i / \sum_j c_j\)</span> instead of raw counts.</li>
|
||||
<li><strong>euclidean</strong>: <span class="arithmatex">\(D = \sqrt{\sum_i (c_i^A - c_i^B)^2}\)</span> on raw counts.</li>
|
||||
<li><strong>relfreq-euclidean</strong>: the same formula on relative frequencies.</li>
|
||||
<li><strong>hellinger</strong>: <span class="arithmatex">\(D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}\)</span> on relative frequencies, bounded in <span class="arithmatex">\([0, 1]\)</span>.</li>
|
||||
<li><strong>hellinger-euclidean</strong>: the unnormalized variant, <span class="arithmatex">\(D = \sqrt{2} \times D_{\text{hellinger}}\)</span>.</li>
|
||||
</ul>
|
||||
<p>Every value routes to one of two independent computations:</p>
|
||||
<h3 id="whole-index-metrics">Whole-index metrics</h3>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Value</th>
|
||||
<th>Definition</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>jaccard</code></td>
|
||||
<td><span class="arithmatex">\(D = 1 - \dfrac{\lvert A \cap B \rvert}{\lvert A \cup B \rvert}\)</span> over the sets of kmers present in each genome</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>mash</code></td>
|
||||
<td>derived from the Jaccard distance via <span class="arithmatex">\(D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)\)</span> where <span class="arithmatex">\(J = 1 - D_{\text{jaccard}}\)</span> and <span class="arithmatex">\(k\)</span> is the index's kmer size; clamped to 1.0 when <span class="arithmatex">\(J \le 0\)</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>hamming</code></td>
|
||||
<td>number of kmer positions where presence differs between the two genomes (presence index only, not normalized): <span class="arithmatex">\(D = \sum_i \mathbb{1}[a_i \ne b_i]\)</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>bray-curtis</code></td>
|
||||
<td><span class="arithmatex">\(D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}\)</span> on raw per-kmer counts</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>relfreq-bray-curtis</code></td>
|
||||
<td>the same formula computed on per-genome relative frequencies <span class="arithmatex">\(p_i = c_i / \sum_j c_j\)</span> instead of raw counts</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>euclidean</code></td>
|
||||
<td><span class="arithmatex">\(D = \sqrt{\sum_i (c_i^A - c_i^B)^2}\)</span> on raw counts</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>relfreq-euclidean</code></td>
|
||||
<td>the same formula on relative frequencies</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>hellinger</code></td>
|
||||
<td><span class="arithmatex">\(D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}\)</span> on relative frequencies, bounded in <span class="arithmatex">\([0, 1]\)</span></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>hellinger-euclidean</code></td>
|
||||
<td>the unnormalized variant, <span class="arithmatex">\(D = \sqrt{2} \times D_{\text{hellinger}}\)</span></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><code>hamming</code> requires a presence/absence index; the others work on either index type.</p>
|
||||
<h3 id="snp-corrections"><code>snp-*</code> corrections</h3>
|
||||
<p>Computed from the central-position SNP model (see "Central-position SNP model" below): a family is the set of up to 4 kmers sharing identical flanking sequence and differing only at the central base. These values require the sibling annex (<code>--sibling-annex</code>, below) and are, by default, computed exhaustively over every non-monomorphic family in the index; add <code>--subsample N</code> to bound the computation to approximately <code>N</code> families instead (see "Sampling at scale" below — the same flag <code>--pseudo-alignment</code>/<code>--sankoff</code> use, but optional here).</p>
|
||||
<p>For a genome pair, let <span class="arithmatex">\(L\)</span> be its total number of eligible loci (both genomes single-copy at that family), <span class="arithmatex">\(p\)</span> the raw proportion of substitutions among those loci, <span class="arithmatex">\(P\)</span>/<span class="arithmatex">\(Q\)</span> the transition/transversion proportions, <span class="arithmatex">\(Q_1\)</span>/<span class="arithmatex">\(Q_2\)</span> Kimura's two transversion categories (A↔C & G↔T vs. A↔T & C↔G), <span class="arithmatex">\(P_1\)</span>/<span class="arithmatex">\(P_2\)</span> the purine (A↔G) / pyrimidine (C↔T) transition proportions, and <span class="arithmatex">\(\pi_A,\pi_C,\pi_G,\pi_T\)</span> the pair's pooled base frequencies.</p>
|
||||
<p><strong><code>snp-raw</code></strong></p>
|
||||
<div class="arithmatex">\[d = p\]</div>
|
||||
<p><strong><code>snp-jc</code></strong></p>
|
||||
<div class="arithmatex">\[d = -\frac{3}{4}\ln\!\left(1-\frac{4p}{3}\right)\]</div>
|
||||
<p><strong><code>snp-k2p</code></strong></p>
|
||||
<div class="arithmatex">\[
|
||||
\begin{aligned}
|
||||
a_1 &= 1-2P-Q \\
|
||||
a_2 &= 1-2Q \\
|
||||
d &= -\frac{1}{2}\ln a_1-\frac{1}{4}\ln a_2
|
||||
\end{aligned}
|
||||
\]</div>
|
||||
<p><strong><code>snp-k81</code></strong></p>
|
||||
<div class="arithmatex">\[
|
||||
\begin{aligned}
|
||||
a_1 &= 1-2P-2Q_1 \\
|
||||
a_2 &= 1-2P-2Q_2 \\
|
||||
a_3 &= 1-2Q_1-2Q_2 \\
|
||||
d &= -\frac{1}{4}\left(\ln a_1+\ln a_2+\ln a_3\right)
|
||||
\end{aligned}
|
||||
\]</div>
|
||||
<p><strong><code>snp-f81</code></strong></p>
|
||||
<div class="arithmatex">\[
|
||||
\begin{aligned}
|
||||
E &= 1-\left(\pi_A^2+\pi_C^2+\pi_G^2+\pi_T^2\right) \\
|
||||
d &= -E\ln\!\left(1-\frac{p}{E}\right)
|
||||
\end{aligned}
|
||||
\]</div>
|
||||
<p><strong><code>snp-t92</code></strong></p>
|
||||
<div class="arithmatex">\[
|
||||
\begin{aligned}
|
||||
g &= \pi_C+\pi_G \\
|
||||
w &= 2g(1-g) \\
|
||||
a_1 &= 1-\frac{P}{w}-Q \\
|
||||
a_2 &= 1-2Q \\
|
||||
d &= -w\ln a_1-\frac{1}{2}(1-w)\ln a_2
|
||||
\end{aligned}
|
||||
\]</div>
|
||||
<p><strong><code>snp-tn93</code></strong></p>
|
||||
<div class="arithmatex">\[
|
||||
\begin{aligned}
|
||||
g_R &= \pi_A+\pi_G \\
|
||||
g_Y &= \pi_C+\pi_T \\
|
||||
k_1 &= \frac{2\pi_A\pi_G}{g_R} \\
|
||||
k_2 &= \frac{2\pi_C\pi_T}{g_Y} \\
|
||||
k_3 &= 2\left(g_Rg_Y-\frac{\pi_A\pi_G\,g_Y}{g_R}-\frac{\pi_C\pi_T\,g_R}{g_Y}\right) \\
|
||||
w_1 &= 1-\frac{P_1}{k_1}-\frac{Q}{2g_R} \\
|
||||
w_2 &= 1-\frac{P_2}{k_2}-\frac{Q}{2g_Y} \\
|
||||
w_3 &= 1-\frac{Q}{2g_Rg_Y} \\
|
||||
d &= -k_1\ln w_1-k_2\ln w_2-k_3\ln w_3
|
||||
\end{aligned}
|
||||
\]</div>
|
||||
<p><strong><code>snp-tv</code></strong> — transversions only, deliberately uncorrected:</p>
|
||||
<div class="arithmatex">\[d = Q\]</div>
|
||||
<p><code>--gamma-shape ALPHA</code> applies to every value above except <code>snp-raw</code> and <code>snp-tv</code>: each <span class="arithmatex">\(-\ln(x)\)</span> term in the formulas above is replaced by <span class="arithmatex">\(\alpha\left(x^{-1/\alpha}-1\right)\)</span> (the same weight, same <span class="arithmatex">\(x\)</span>).</p>
|
||||
<h3 id="output">Output</h3>
|
||||
<p>Without <code>-o</code>, the matrix goes to stdout in relaxed-PHYLIP format (<code>n</code> on the first line, then one <code>label<TAB>value...</code> row per genome). With <code>--csv</code>, the format is instead a header row <code>genome,<label1>,<label2>,...</code> followed by one <code><label>,<value1>,<value2>,...</code> row per genome, 6 decimals. Both formats are symmetric with a zero diagonal, except where noted below.</p>
|
||||
<h2 id="-exclude-genome-min-shared-family"><code>--exclude-genome</code>, <code>--min-shared-family</code></h2>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Option</th>
|
||||
<th>Description</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>--exclude-genome LABEL</code></td>
|
||||
<td>Exclude a genome (repeatable). Drops its row/column from the distance/shared-kmer matrix output, and removes it from the sampling used by <code>--pseudo-alignment</code>/<code>--sankoff</code>/a <code>snp-*</code> <code>--distance</code> value. Does not change the value computed for any remaining pair</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--min-shared-family N</code></td>
|
||||
<td>Auto-exclude, on top of <code>--exclude-genome</code>, any genome whose mean shared-family count against every other genome (see "Family Overlap" below) falls below <code>N</code>. Applies only to <code>--pseudo-alignment</code>/<code>--sankoff</code>/<code>snp-*</code> <code>--distance</code> — never to the whole-index metrics or their matrix/NJ/UPGMA output</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p>Neighbor-Joining and UPGMA trees (<code>--nj</code>/<code>--upgma</code>) are always built from every genome in the index, regardless of <code>--exclude-genome</code>/<code>--min-shared-family</code>.</p>
|
||||
<h2 id="central-position-snp-model">Central-position SNP model</h2>
|
||||
<p>This is a separate operation from the distance-matrix computation above: if any option below is used, no <code>--metric</code> matrix is computed in the same invocation.</p>
|
||||
<p>A <strong>family</strong> is the set of up to 4 kmers that share identical flanking sequence and differ only at the exact central base. Because <span class="arithmatex">\(k\)</span> is odd, the central position is well defined and maps to itself under reverse complementation. All computations below first require building the <strong>sibling annex</strong>, an index-wide record of which of the 4 possible central bases are observed at each family, across every genome.</p>
|
||||
<p>Requires the sibling annex, built once per index:</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Option</th>
|
||||
<th>Default</th>
|
||||
<th>Description</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>--sibling-annex</code></td>
|
||||
<td>off</td>
|
||||
<td>Build the sibling presence-mask annex (prerequisite for every option below)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--exclude-genome LABEL</code></td>
|
||||
<td>none</td>
|
||||
<td>Exclude a genome (repeatable) from every SNP/Sankoff/export computation below</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--min-shared-family N</code></td>
|
||||
<td>none</td>
|
||||
<td>Auto-exclude any genome whose mean shared-family count against every other genome (see <code>--family-overlap</code>) falls below <code>N</code> — same exclusion as <code>--exclude-genome</code>, applied on top of it</td>
|
||||
<td>Build (or rebuild) the sibling-count/minorant annex — prerequisite for every option in this section, and for a <code>snp-*</code> <code>--distance</code> value</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--sibling-stats</code></td>
|
||||
<td>off</td>
|
||||
<td>Write the family-size (sibling count) distribution, per genome and globally</td>
|
||||
<td>Write <code><prefix>_siblings.csv</code>: the family-size distribution, per genome and globally</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--raw-snp-distance</code></td>
|
||||
<td>off</td>
|
||||
<td>Write the single-copy central-SNP p-distance matrix</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--raw-snp-counts</code></td>
|
||||
<td>off</td>
|
||||
<td>Write per-pair diagnostic counts (n_snp, n_shared, n_eligible) instead of a matrix</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--snp</code></td>
|
||||
<td>off</td>
|
||||
<td>Write a SNP-only pseudo-alignment in FASTA, IUPAC-coded</td>
|
||||
<td><code>--sibling-hist</code></td>
|
||||
<td>Print the global family-size histogram (1-4 members) only</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--family-overlap</code></td>
|
||||
<td>off</td>
|
||||
<td>Write an NxN matrix of, for each genome pair, how many variable families both genomes actually carry a call for; the diagonal holds each genome's own total family count</td>
|
||||
<td>Write <code><prefix>_family_overlap.csv</code>: for every genome pair, how many variable families both genomes carry a call for</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--shannon</code></td>
|
||||
<td>Write <code><prefix>_entropy.csv</code>: per-family Shannon entropy, one row per family, full unsampled scan</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--pseudo-alignment</code></td>
|
||||
<td>Write <code><prefix>_alignment.fasta</code>: a SNP-only pseudo-alignment. Requires <code>--subsample N</code></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<h3 id="locus-eligibility">Locus eligibility</h3>
|
||||
<p>A family is eligible for a genome pair <span class="arithmatex">\((i, j)\)</span> only if genome <span class="arithmatex">\(i\)</span> carries exactly one of the family's observed forms (single-copy, unambiguous) and genome <span class="arithmatex">\(j\)</span> also carries exactly one. A genome carrying more than one form at a locus makes that locus ineligible for any pair involving it.</p>
|
||||
<p><code>--raw-snp-distance</code> tallies, over every eligible locus of every genome pair, <span class="arithmatex">\(n_{\text{snp}}\)</span> (the two genomes' single forms differ) versus <span class="arithmatex">\(n_{\text{shared}}\)</span> (they agree — this includes invariant families). The output ratio is <span class="arithmatex">\(\hat{p} = \dfrac{n_{\text{snp}}}{n_{\text{snp}} + n_{\text{shared}}}\)</span>.</p>
|
||||
<p><code>--snp</code> restricts itself to <em>variable</em> families (family size <span class="arithmatex">\(\ge 2\)</span>) and writes one FASTA record per genome, one column per family, IUPAC-coded from each genome's presence mask at that family (a single form → the plain base; several forms → the matching IUPAC ambiguity code; no form → <code>-</code>).</p>
|
||||
<p><code>--exclude-genome</code> removes a genome from these computations, re-checking column variability among the remaining genomes so that a column made monomorphic by the exclusion is dropped rather than kept artificially. It does not affect the <code>--metric</code> distance-matrix path.</p>
|
||||
<h3 id="family-overlap-and-low-coverage-genomes">Family overlap and low-coverage genomes</h3>
|
||||
<p><code>--family-overlap</code> writes, for every genome pair, how many variable families both genomes actually carry a call for (neither is absent) — a direct measure of how much informative content two genomes actually share. On genome-skim or otherwise incomplete-coverage collections, a genome with very little overlap with everything else has almost nothing left to constrain its position in a tree, and tends to end up placed unstably (near-zero branch length, grafted inside an unrelated clade) by <code>--tnt</code>/<code>--iqtree</code>.</p>
|
||||
<p><code>--min-shared-family N</code> automates the fix: it excludes, before any computation, every genome whose mean shared-family count against all other genomes (the same statistic, averaged per row of the <code>--family-overlap</code> matrix) falls below <code>N</code>. There is no universal value for <code>N</code> — it depends on how divergent and how completely covered the genome collection is; inspect <code>--family-overlap</code>'s own output to find where the real gap sits before choosing a threshold.</p>
|
||||
<h2 id="sampling-at-scale-subsample-shannon-entropy">Sampling at scale: <code>--subsample</code>, <code>--shannon</code>, <code>--entropy</code></h2>
|
||||
<p>On a large index (billions of families), building a full pseudo-alignment or fully calibrating the Sankoff model is not just slow — it may not fit in the time you have. <code>--subsample</code> bounds the work to a fixed number of families; <code>--shannon</code> reports how informative each family is; <code>--entropy</code>/<code>--entropy-sd</code> bias which families get kept toward the informative ones instead of choosing uniformly at random.</p>
|
||||
<p>A family is eligible for a genome pair <span class="arithmatex">\((i,j)\)</span> only if both genomes carry exactly one of its observed forms (single-copy, unambiguous).</p>
|
||||
<h3 id="-sibling-stats"><code>--sibling-stats</code></h3>
|
||||
<p><code><prefix>_siblings.csv</code> — family size = number of distinct central bases observed at a family (1-4).</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Column</th>
|
||||
<th>Meaning</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>genome</code></td>
|
||||
<td>genome label, or the literal <code>global</code> for the last row</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>1</code>, <code>2</code>, <code>3</code>, <code>4</code></td>
|
||||
<td>for a genome row: number of families of that size where the genome carries ≥ 1 member. For the <code>global</code> row: the actual deduplicated family-size histogram — not the sum of the rows above</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<h3 id="family-overlap">Family Overlap</h3>
|
||||
<p><code>--family-overlap</code> writes <code><prefix>_family_overlap.csv</code>: header <code>genome,<label1>,<label2>,...</code>, one row per genome, cell <code>[i][j]</code> = number of variable families (family size ≥ 2) where both genome <code>i</code> and genome <code>j</code> carry a call. The diagonal is always <code>0</code>. Every genome is written, unfiltered by <code>--exclude-genome</code>/<code>--min-shared-family</code>.</p>
|
||||
<p><code>--min-shared-family N</code> uses the mean of each genome's own row (excluding the diagonal) against this same matrix as its exclusion statistic. There is no universal value for <code>N</code> — inspect <code>--family-overlap</code>'s own output to find where the real gap sits in a given genome collection before choosing a threshold.</p>
|
||||
<h3 id="-pseudo-alignment"><code>--pseudo-alignment</code></h3>
|
||||
<p><code><prefix>_alignment.fasta</code> — one record per non-excluded genome, one column per variable family (family size ≥ 2). Each site is IUPAC-coded from the genome's presence mask at that family: a single observed form → the plain base; several forms → the matching IUPAC ambiguity code; no form → <code>-</code>.</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
@@ -1705,13 +1867,18 @@
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>--subsample N</code></td>
|
||||
<td>none (keep everything)</td>
|
||||
<td>Cap the number of variable families (family size ≥ 2) retained, to approximately <code>N</code></td>
|
||||
<td>none (mandatory here)</td>
|
||||
<td>Target number of families to sample</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--shannon</code></td>
|
||||
<td><code>--free-loss</code></td>
|
||||
<td>off</td>
|
||||
<td>Write <code><prefix>_shannon.csv</code>: per-family Shannon entropy, one row per family</td>
|
||||
<td>Treat a genome carrying none of a family's observed members as missing data (<code>?</code>) instead of <code>-</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--no-ambiguity</code></td>
|
||||
<td>off</td>
|
||||
<td>Treat a genome carrying more than one member of a family as missing data (<code>?</code>) instead of an IUPAC ambiguity code</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--entropy MU</code></td>
|
||||
@@ -1725,13 +1892,11 @@
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><code>--subsample</code>/<code>--entropy</code>/<code>--entropy-sd</code> affect every option that scans variable families: <code>--snp</code>, <code>--family-overlap</code>, <code>--shannon</code>, and the whole Sankoff pipeline (<code>--sankoff</code>/<code>--tnt</code>/<code>--phyg</code>/<code>--iqtree</code>, next section) — all of them draw from the <em>same</em> selection of families in one invocation, so the Sankoff calibration and the pseudo-alignment it calibrates always describe the same sites, and <code>--family-overlap</code>'s counts stay consistent with <code>--snp</code>'s columns. <code>--raw-snp-distance</code>/<code>--raw-snp-counts</code> are not affected — they always scan every family, since their p-distance estimate is a whole-index statistic, not something that benefits from being restricted to a sample.</p>
|
||||
<h3 id="-subsample-n"><code>--subsample N</code></h3>
|
||||
<p>Without <code>--subsample</code>, every variable family (family size ≥ 2, i.e. every family where at least one genome differs from the rest) is used. With <code>--subsample N</code>, roughly <code>N</code> families are kept instead, chosen at random but in proportion to how many candidate families each part of the index actually holds — so the sample stays representative of the whole index, not skewed toward whichever part happens to be scanned first. If the index has fewer than <code>N</code> candidate families in the first place, <code>--subsample</code> has no effect: everything is kept.</p>
|
||||
<p><code>--subsample</code> trades completeness for speed: <code>--snp</code>'s alignment gets fewer columns, <code>--sankoff</code>'s calibration is based on fewer observations, but the resolution work (the expensive part of a <code>phylo</code> run on a large index) scales with <code>N</code> instead of with the index's true size. Pick <code>N</code> as large as your time budget allows — a few hundred thousand to a few million families is usually enough for the transition-probability estimates in <code>--sankoff</code>'s calibration to stabilize; a smaller <code>N</code> speeds up exploratory runs.</p>
|
||||
<h2 id="sampling-at-scale-subsample-shannon-entropy">Sampling at scale: <code>--subsample</code>, <code>--shannon</code>, <code>--entropy</code></h2>
|
||||
<p><code>--subsample</code>, <code>--free-loss</code>, <code>--no-ambiguity</code>, <code>--entropy</code>/<code>--entropy-sd</code> are shared by <code>--pseudo-alignment</code>, <code>--sankoff</code> (and everything it implies: <code>--tnt</code>/<code>--phyg</code>/<code>--iqtree</code>), and a <code>snp-*</code> <code>--distance</code> value — one draw feeds all of them in a single invocation. <code>--subsample</code> is mandatory for <code>--pseudo-alignment</code>/<code>--sankoff</code>; for a <code>snp-*</code> <code>--distance</code> value it is optional (omitted means every non-monomorphic family in the index, not an approximation).</p>
|
||||
<p>Without <code>--subsample</code>, every variable family (family size ≥ 2) is used. With <code>--subsample N</code>, roughly <code>N</code> families are kept instead, drawn in proportion to how many candidate families each part of the index actually holds, so the sample stays representative of the whole index. If the index has fewer than <code>N</code> candidate families, <code>--subsample</code> has no effect.</p>
|
||||
<h3 id="-shannon-measuring-how-informative-a-family-is"><code>--shannon</code>: measuring how informative a family is</h3>
|
||||
<p>Not every variable family is equally useful for a tree: a family that differs in only one genome out of a thousand carries very little signal, and one where the pattern looks essentially random across genomes may be too saturated (multiple substitutions have overwritten the original signal) to carry real information either. <code>--shannon</code> quantifies this with the Shannon entropy (in bits) of each family's states across the genomes that carry it — low entropy means "almost everyone agrees" (an invariant or near-invariant family, phylogenetically shallow), while entropy near the ceiling for a 4-state character means "close to a random draw between the possible bases" (saturated).</p>
|
||||
<p><code><prefix>_shannon.csv</code> (or <code>shannon.csv</code> without <code>-o</code>) has one row per family visited:</p>
|
||||
<p><code><prefix>_entropy.csv</code> has one row per family visited:</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
@@ -1750,15 +1915,15 @@
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>entropy15</code></td>
|
||||
<td>Shannon entropy (bits) over the 16 possible states (the 15 non-empty subsets of <code>{A,C,G,T}</code> — the same alphabet <code>--sankoff</code>'s 16-state model uses), genomes absent from the family excluded from the count</td>
|
||||
<td>Shannon entropy (bits) over the 16 possible states (the 15 non-empty subsets of <code>{A,C,G,T}</code>), genomes absent from the family excluded from the count</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>entropy4</code></td>
|
||||
<td>Shannon entropy (bits) reduced to the 4 plain bases, kept alongside <code>entropy15</code> for comparison — a genome carrying more than one base at once counts once per base, so this can differ from <code>entropy15</code></td>
|
||||
<td>Shannon entropy (bits) reduced to the 4 plain bases, kept alongside <code>entropy15</code> for comparison</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>family_size</code></td>
|
||||
<td>number of distinct central bases observed anywhere in the index for this family (2–4, since monomorphic families aren't visited)</td>
|
||||
<td>number of distinct central bases observed anywhere in the index for this family (2-4)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>n_genomes_present</code></td>
|
||||
@@ -1766,11 +1931,10 @@
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p>Run with <code>--subsample N --shannon</code> to get a bounded diagnostic sample instead of a full-index pass — useful to inspect the entropy distribution and decide reasonable <code>--entropy</code>/<code>--entropy-sd</code> values (see below) before committing to a full run.</p>
|
||||
<p>Run with <code>--subsample N --shannon</code> to get a bounded diagnostic sample instead of a full-index pass — useful for choosing <code>--entropy</code>/<code>--entropy-sd</code> values before a full run.</p>
|
||||
<h3 id="-entropy-mu-entropy-sd-sigma-biasing-the-sample-toward-informative-families"><code>--entropy MU</code> / <code>--entropy-sd SIGMA</code>: biasing the sample toward informative families</h3>
|
||||
<p>By default, <code>--subsample</code> draws families uniformly — every candidate family has the same chance of being kept, regardless of how informative it actually is. <code>--entropy</code>/<code>--entropy-sd</code> change that: instead of a uniform draw, each family's chance of being kept is weighted by how close its own entropy (the <code>entropy15</code> value <code>--shannon</code> reports) is to <code>MU</code>, using a bell-shaped (Gaussian) curve of width <code>SIGMA</code> — a family with entropy exactly <code>MU</code> is the most likely to be kept, and the chance falls off smoothly the further its entropy is from <code>MU</code>, with no hard cutoff (a few families outside the target band can still get in, just less often).</p>
|
||||
<p>The filter activates as soon as <em>either</em> <code>--entropy</code> or <code>--entropy-sd</code> is given; whichever one you don't set defaults to <code>1.0</code>/<code>0.5</code>. It can be combined with <code>--subsample N</code> (the target count is still approximately <code>N</code>, now biased toward the entropy band instead of uniform — expect somewhat <em>fewer</em> than <code>N</code> families in practice, since low-weight families are dropped rather than replaced) or used alone (<code>--entropy</code> without <code>--subsample</code>: a soft entropy filter over the whole index, no size target).</p>
|
||||
<p>The first <code>phylo</code> run on a given index that uses <code>--entropy</code>/<code>--entropy-sd</code> pays a one-time extra cost (every candidate family's entropy has to be computed once, up front, and is then saved alongside the index). Every following run — even with different <code>MU</code>/<code>SIGMA</code> values — reuses that saved data and stays fast.</p>
|
||||
<p>By default, <code>--subsample</code> draws families uniformly. With <code>--entropy</code>/<code>--entropy-sd</code>, each family's chance of being kept is instead weighted by how close its own entropy (<code>entropy15</code>) is to <code>MU</code>, using a Gaussian curve of width <code>SIGMA</code> — no hard cutoff. The filter activates as soon as either flag is given; the other defaults to <code>1.0</code>/<code>0.5</code>. Combine with <code>--subsample N</code> (expect somewhat fewer than <code>N</code> families kept in practice) or use alone (a soft filter over the whole index, no size target).</p>
|
||||
<p>The first <code>phylo</code> run on a given index that uses <code>--entropy</code>/<code>--entropy-sd</code> pays a one-time extra cost (every candidate family's entropy is computed once and saved alongside the index); later runs, even with different <code>MU</code>/<code>SIGMA</code>, reuse that saved data.</p>
|
||||
<h2 id="sankoff-calibration-and-phylogenetic-exports">Sankoff calibration and phylogenetic exports</h2>
|
||||
<table>
|
||||
<thead>
|
||||
@@ -1784,17 +1948,17 @@
|
||||
<tr>
|
||||
<td><code>--sankoff</code></td>
|
||||
<td>off</td>
|
||||
<td>Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment</td>
|
||||
<td>Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment. Requires <code>--subsample N</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--sankoff-ratio-ceiling</code></td>
|
||||
<td><code>0.5</code></td>
|
||||
<td>Exclude genome pairs whose raw SNP ratio exceeds this value from the calibration</td>
|
||||
<td>Exclude genome pairs whose raw SNP ratio exceeds this value from the base-composition part of the calibration</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--free-loss</code></td>
|
||||
<td>off</td>
|
||||
<td>Recode a family's non-detection as the <code>?</code> missing-data symbol instead of an ordinary, costed state, in <code>--sankoff</code>'s pseudo-alignment and every export built from it</td>
|
||||
<td>Recode a family's non-detection as the <code>?</code> missing-data symbol instead of an ordinary, costed state, throughout <code>--sankoff</code> and every export built from it</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>--tnt</code></td>
|
||||
@@ -1819,30 +1983,25 @@
|
||||
<tr>
|
||||
<td><code>--sankoff-cost-scale</code></td>
|
||||
<td><code>100</code></td>
|
||||
<td>Integer scaling factor applied to costs before rounding (required by TNT/PhyG's integer-only cost commands)</td>
|
||||
<td>Integer scaling factor applied to costs before rounding, for TNT/PhyG's integer-only cost commands</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<h3 id="the-16-state-model">The 16-state model</h3>
|
||||
<p>Each family is treated as a character with 16 possible states: one per subset of the 4 possible central bases actually observed (including the empty subset). Calibration combines two tallies, both restricted to genome pairs at or below <code>--sankoff-ratio-ceiling</code>:</p>
|
||||
<ul>
|
||||
<li>a <span class="arithmatex">\(5 \times 5\)</span> transition matrix over family cardinality (0–4 observed forms) between paired genomes, and</li>
|
||||
<li>a <span class="arithmatex">\(4 \times 4\)</span> base-substitution transition matrix from unambiguous single-copy loci,</li>
|
||||
</ul>
|
||||
<p>which are combined into a row-normalized <span class="arithmatex">\(16 \times 16\)</span> transition probability matrix <span class="arithmatex">\(P\)</span>, converted to a symmetric cost matrix via <span class="arithmatex">\(\text{cost}(a,b) = -\ln P(a,b)\)</span>.</p>
|
||||
<p><code>--sankoff</code> alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol <code>0</code> (never a gap character, to avoid ambiguity with external tools' own gap semantics). It does not run any external tool.</p>
|
||||
<p>With <code>--free-loss</code>, the empty state is recoded to <code>?</code> instead — TNT/PhyG/IQ-TREE's own missing-data symbol — rather than an ordinary, costed 16th state. This matters for genome-skim or otherwise incomplete-coverage collections, where non-detection of a family is dominated by sampling failure rather than true evolutionary loss: scoring it as a real state risks grouping genomes by shared undersampling instead of shared ancestry. <code>?</code> rather than <code>-</code> because <code>-</code> still carries gap/indel semantics in these tools, and a non-detected family is not an observed deletion. <code>--free-loss</code> also drops the cardinality-transition cost between any two states, not just to/from the empty one: whether a genome shows 1 vs. 2 (etc.) detected members of a family it does carry is exactly as vulnerable to sampling failure as whether the family was detected at all, so gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family. Combine with <code>--min-shared-family</code>/<code>--family-overlap</code> above: <code>--free-loss</code> removes the false signal from non-detection, but a genome left with too little real overlap with everything else will still be placed unstably — excluding it is the other half of the fix.</p>
|
||||
<p>Each family is a character with 16 possible states: one per subset of the 4 possible central bases (including the empty subset). Calibration combines a <span class="arithmatex">\(5 \times 5\)</span> transition matrix over family cardinality (0-4 observed forms) and a <span class="arithmatex">\(4 \times 4\)</span> base-substitution matrix from unambiguous single-copy loci, both restricted to genome pairs at or below <code>--sankoff-ratio-ceiling</code>, into a row-normalized <span class="arithmatex">\(16 \times 16\)</span> transition probability matrix <span class="arithmatex">\(P\)</span>, converted to a symmetric cost matrix via <span class="arithmatex">\(\text{cost}(a,b) = -\ln P(a,b)\)</span>.</p>
|
||||
<p><code>--sankoff</code> alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol <code>0</code> (never a gap character). It does not run any external tool.</p>
|
||||
<p>With <code>--free-loss</code>, the empty state is recoded to <code>?</code> (TNT/PhyG/IQ-TREE's own missing-data symbol) instead of an ordinary, costed 16th state — <code>?</code> rather than <code>-</code>, since <code>-</code> still carries gap/indel semantics in these tools. <code>--free-loss</code> also zeroes the cardinality-transition cost between any two states, not just to/from the empty one: gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family.</p>
|
||||
<h3 id="exports">Exports</h3>
|
||||
<p>All three exports reuse the <code>--sankoff</code> calibrated matrix and pseudo-alignment, recoded for the target tool:</p>
|
||||
<ul>
|
||||
<li><strong><code>--tnt</code></strong>: a self-contained TNT script (alignment recoded to TNT's fixed 16-symbol alphabet, integer-scaled cost matrix re-closed to a metric, a default search block).</li>
|
||||
<li><strong><code>--phyg</code></strong>: a custom cost-matrix file plus a PhyG script reusing the <code>--sankoff</code> alignment directly.</li>
|
||||
<li><strong><code>--iqtree</code></strong>: a custom substitution-model file (exchangeability matrix recovered as <span class="arithmatex">\(R(a,b) = e^{-\text{cost}(a,b)}\)</span>, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths (unlike the parsimony step-counts from TNT/PhyG). Only states actually occurring in the alignment are kept and compactly renumbered.</li>
|
||||
<li><strong><code>--iqtree</code></strong>: a custom substitution-model file (exchangeability matrix recovered as <span class="arithmatex">\(R(a,b) = e^{-\text{cost}(a,b)}\)</span>, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths. Only states actually occurring in the alignment are kept and compactly renumbered.</li>
|
||||
</ul>
|
||||
<p>TNT and PhyG both write trees with bare numeric leaf labels (<code>1</code>, <code>2</code>, …, in the same order as <code><prefix>_sankoff.fasta</code>). Use <a href="../name-tree/"><code>name-tree</code></a> on the tool's own tree output plus that same FASTA to get a NEXUS file with real taxon names.</p>
|
||||
<p>TNT and PhyG both write trees with bare numeric leaf labels (<code>1</code>, <code>2</code>, ..., in the order the genomes appear in <code><prefix>_sankoff.fasta</code>).</p>
|
||||
<h2 id="output-files">Output files</h2>
|
||||
<p>With <code>-o/--output PREFIX</code>, the relevant subset of the files below is written. Without <code>-o</code>, only the plain <code>--metric</code> distance matrix is produced, on stdout. All matrices use genome labels (from the index metadata) as row/column headers, in index order; all CSVs are comma-separated with a header row.</p>
|
||||
<h3 id="distance-matrix_1">Distance matrix</h3>
|
||||
<p>With <code>-o/--output PREFIX</code>, the relevant subset of the files below is written. Without <code>-o</code>, only the distance matrix is produced, on stdout. All matrices use genome labels as row/column headers, in index order.</p>
|
||||
<h3 id="distance-matrix">Distance matrix</h3>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
@@ -1854,10 +2013,16 @@
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code><prefix>_dist.phy</code></td>
|
||||
<td>always, unless <code>--csv</code></td>
|
||||
<td>relaxed PHYLIP</td>
|
||||
<td>the <code>--distance</code> matrix</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_dist.csv</code></td>
|
||||
<td>always</td>
|
||||
<td><code>--csv</code></td>
|
||||
<td>CSV matrix</td>
|
||||
<td>the <code>--metric</code> distance, 6 decimals, symmetric, diagonal 0</td>
|
||||
<td>the <code>--distance</code> matrix, 6 decimals</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_shared.csv</code></td>
|
||||
@@ -1869,17 +2034,17 @@
|
||||
<td><code><prefix>_nj.nwk</code></td>
|
||||
<td><code>--nj</code></td>
|
||||
<td>Newick</td>
|
||||
<td>Neighbor-Joining tree, branch lengths from the <code>--metric</code> matrix</td>
|
||||
<td>Neighbor-Joining tree</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_upgma.nwk</code></td>
|
||||
<td><code>--upgma</code></td>
|
||||
<td>Newick</td>
|
||||
<td>UPGMA tree, same matrix</td>
|
||||
<td>UPGMA tree</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p>Matrix layout (<code>_dist.csv</code>, <code>_shared.csv</code>, and every other "CSV matrix" below): header <code>genome,<label1>,<label2>,...</code>, one data row per genome, <code><label>,<value1>,<value2>,...</code>.</p>
|
||||
<p>CSV matrix layout (<code>_dist.csv</code>, <code>_shared.csv</code>, <code>_family_overlap.csv</code>): header <code>genome,<label1>,<label2>,...</code>, one data row per genome, <code><label>,<value1>,<value2>,...</code>.</p>
|
||||
<h3 id="central-position-snp-model_1">Central-position SNP model</h3>
|
||||
<table>
|
||||
<thead>
|
||||
@@ -1898,90 +2063,25 @@
|
||||
<td>family-size distribution, per genome and global</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_rawsnp.csv</code></td>
|
||||
<td><code>--raw-snp-distance</code></td>
|
||||
<td>CSV matrix</td>
|
||||
<td>single-copy central-SNP p-distance (<span class="arithmatex">\(\hat p\)</span>), or <code>NA</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_rawsnp_counts.csv</code></td>
|
||||
<td><code>--raw-snp-counts</code></td>
|
||||
<td>CSV table</td>
|
||||
<td>per-pair diagnostic counts behind <code>_rawsnp.csv</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_snp.fasta</code></td>
|
||||
<td><code>--snp</code></td>
|
||||
<td>FASTA</td>
|
||||
<td>SNP-only pseudo-alignment, IUPAC-coded</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_family_overlap.csv</code></td>
|
||||
<td><code>--family-overlap</code></td>
|
||||
<td>CSV matrix</td>
|
||||
<td>variable families both genomes of a pair carry a call for</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_shannon.csv</code></td>
|
||||
<td><code><prefix>_entropy.csv</code></td>
|
||||
<td><code>--shannon</code></td>
|
||||
<td>CSV table</td>
|
||||
<td>per-family Shannon entropy, see "Sampling at scale" above</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><strong><code>_siblings.csv</code></strong> — family size = number of distinct central bases observed at a family (1–4), not "sibling count" (0–3).</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Column</th>
|
||||
<th>Meaning</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>genome</code></td>
|
||||
<td>genome label, or the literal <code>global</code> for the last row</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>1</code>, <code>2</code>, <code>3</code>, <code>4</code></td>
|
||||
<td>for a genome row: number of families of that size where the genome carries ≥ 1 member. For the <code>global</code> row: the actual deduplicated family-size histogram — <strong>not</strong> the sum of the rows above (a family shared by several genomes would otherwise be counted once per genome)</td>
|
||||
<td><code><prefix>_alignment.fasta</code></td>
|
||||
<td><code>--pseudo-alignment</code></td>
|
||||
<td>FASTA</td>
|
||||
<td>SNP-only pseudo-alignment, IUPAC-coded</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><strong><code>_rawsnp.csv</code></strong> — same matrix layout as <code>_dist.csv</code>; each cell is <span class="arithmatex">\(\hat p = n_{\text{snp}}/(n_{\text{snp}}+n_{\text{shared}})\)</span>, 6 decimals, or <code>NA</code> when the pair has zero eligible loci (distinguishes "identical everywhere eligible" from "nothing eligible at all").</p>
|
||||
<p><strong><code>_rawsnp_counts.csv</code></strong> — one row per unordered genome pair (not a matrix), the counts <code>_rawsnp.csv</code>'s ratio is computed from:</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
<th>Column</th>
|
||||
<th>Meaning</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<td><code>genome_a</code>, <code>genome_b</code></td>
|
||||
<td>the pair</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>n_snp</code></td>
|
||||
<td>eligible loci where the two genomes' single forms differ</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>n_shared</code></td>
|
||||
<td>eligible loci where they agree (includes invariant families)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>n_eligible</code></td>
|
||||
<td><code>n_snp + n_shared</code></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>ratio</code></td>
|
||||
<td><span class="arithmatex">\(\hat p\)</span> = <code>n_snp / n_eligible</code>, or <code>NA</code> if <code>n_eligible = 0</code></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><strong><code>_snp.fasta</code></strong> — one record per non-excluded genome, one column per variable family (family size ≥ 2), header carries an <code>n_sites</code> annotation. Each site is IUPAC-coded from the genome's presence mask at that family: single observed form → plain base; several forms → matching IUPAC ambiguity code; no form → <code>-</code>.</p>
|
||||
<p><strong><code>_family_overlap.csv</code></strong> — same matrix layout as <code>_dist.csv</code>; cell <code>[i][j]</code> = number of <code>_snp.fasta</code> columns where both genome <code>i</code> and <code>j</code> carry a call (neither is <code>-</code>). Diagonal <code>[i][i]</code> is kept (not skipped): it holds genome <code>i</code>'s own total variable-family count.</p>
|
||||
<h3 id="sankoff-calibration-and-exports">Sankoff calibration and exports</h3>
|
||||
<table>
|
||||
<thead>
|
||||
@@ -2009,7 +2109,7 @@
|
||||
<td><code><prefix>_sankoff.fasta</code></td>
|
||||
<td>same flags</td>
|
||||
<td>FASTA</td>
|
||||
<td>Sankoff-recoded pseudo-alignment</td>
|
||||
<td>Sankoff-recoded pseudo-alignment, header carries an <code>n_sites</code> annotation</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code><prefix>_sankoff.tnt</code></td>
|
||||
@@ -2049,8 +2149,8 @@
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><strong><code>_sankoff_matrix.csv</code></strong> — header <code>state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N</code>: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), <code>0</code> standing for the empty/absent state (not <code>-</code>, to avoid colliding with external tools' own gap syntax). One row per source state, one value per destination state, cost <span class="arithmatex">\(-\ln P(a,b)\)</span>, 4 decimals.</p>
|
||||
<p><strong><code>_sankoff_params.yaml</code></strong> — everything the calibration estimated, structured so it can be reloaded rather than re-parsed:</p>
|
||||
<p><strong><code>_sankoff_matrix.csv</code></strong> — header <code>state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N</code>: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), <code>0</code> standing for the empty/absent state. One row per source state, one value per destination state, cost <span class="arithmatex">\(-\ln P(a,b)\)</span>, 4 decimals.</p>
|
||||
<p><strong><code>_sankoff_params.yaml</code></strong>:</p>
|
||||
<table>
|
||||
<thead>
|
||||
<tr>
|
||||
@@ -2065,7 +2165,7 @@
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>cardinality_transitions</code></td>
|
||||
<td>5×5 list of <code>{from, to, count, probability}</code>, family cardinality (0–4 observed forms)</td>
|
||||
<td>5×5 list of <code>{from, to, count, probability}</code>, family cardinality (0-4 observed forms)</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td><code>composition_transitions</code></td>
|
||||
@@ -2073,17 +2173,17 @@
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
<p><strong><code>_sankoff.fasta</code></strong> — same sites as <code>_snp.fasta</code>, recoded to match <code>_sankoff_matrix.csv</code>'s alphabet: absent state is <code>0</code> (or <code>?</code> under <code>--free-loss</code>). Excluded genomes dropped; columns left monomorphic by that exclusion are re-checked and dropped too.</p>
|
||||
<p><strong><code>_sankoff.tnt</code></strong> (<code>--tnt</code>) — self-contained TNT script: <code>xread</code> block (alignment recoded to TNT's fixed <code>0-9A-F</code> alphabet), an integer-scaled (<code>--sankoff-cost-scale</code>) and metric-closed <code>smatrix</code>, a default <code>hold 20; mult; export</code> search. Run with <code>printf 'proc <path>;\nquit;\n' | tnt</code>. Produces <code><prefix>_sankoff.tre</code> (bare numeric leaf labels, order matching <code>_sankoff.fasta</code>) — feed both into <a href="../name-tree/"><code>name-tree</code></a> to recover taxon names.</p>
|
||||
<p><strong><code>_sankoff.tcm</code></strong> (<code>--phyg</code>) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state <code>0</code>, since it's never actually triggered).</p>
|
||||
<p><strong><code>_sankoff.pg</code></strong> (<code>--phyg</code>) — script: <code>read(prefasta:..., tcm:...)</code> against <code>_sankoff.fasta</code>/<code>_sankoff.tcm</code>, a default 300s/4-instance <code>search</code>, <code>report(...)</code> writing <code><prefix>_sankoff.tre</code> (bare numeric labels, as for <code>--tnt</code>). Run with <code>phyg</code> from the output directory (the script uses relative file names). Feed the tree plus <code>_sankoff.fasta</code> into <a href="../name-tree/"><code>name-tree</code></a> for taxon names.</p>
|
||||
<p><strong><code>_sankoff.fasta</code></strong> — recoded to match <code>_sankoff_matrix.csv</code>'s alphabet: absent state is <code>0</code> (or <code>?</code> under <code>--free-loss</code>). Excluded genomes dropped; columns left monomorphic by that exclusion are re-checked and dropped too.</p>
|
||||
<p><strong><code>_sankoff.tnt</code></strong> (<code>--tnt</code>) — <code>xread</code> block (alignment recoded to TNT's fixed <code>0-9A-F</code> alphabet), an integer-scaled (<code>--sankoff-cost-scale</code>) and metric-closed <code>smatrix</code>, a default <code>hold 20; mult; export</code> search. Run with <code>printf 'proc <path>;\nquit;\n' | tnt</code>. Produces <code><prefix>_sankoff.tre</code> (bare numeric leaf labels, order matching <code>_sankoff.fasta</code>).</p>
|
||||
<p><strong><code>_sankoff.tcm</code></strong> (<code>--phyg</code>) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state <code>0</code>).</p>
|
||||
<p><strong><code>_sankoff.pg</code></strong> (<code>--phyg</code>) — script: <code>read(prefasta:..., tcm:...)</code> against <code>_sankoff.fasta</code>/<code>_sankoff.tcm</code>, a default 300s/4-instance <code>search</code>, <code>report(...)</code> writing <code><prefix>_sankoff.tre</code>. Run with <code>phyg</code> from the output directory (the script uses relative file names).</p>
|
||||
<p><strong><code>_iqtree.model</code></strong> (<code>--iqtree</code>) — lower-triangular exchangeability matrix <span class="arithmatex">\(R(a,b) = e^{-\text{cost}(a,b)}\)</span> (one row of increasing length per state, whitespace-separated, PAML order), followed by one line of empirical state frequencies. Only states actually occurring in the alignment are kept, compactly renumbered <code>0..k-1</code>.</p>
|
||||
<p><strong><code>_iqtree.fasta</code></strong> (<code>--iqtree</code>) — alignment recoded to that same compact <code>0..k-1</code> alphabet (symbols <code>0-9A-F</code>). Under <code>--free-loss</code>, non-detection becomes <code>?</code> and columns left non-informative once missing calls are ignored are dropped first (required for <code>+ASC</code>); with <code>--iqtree-min-freq</code> also set (the default), any state rarer than that threshold is folded into the same <code>?</code> treatment, and non-informative columns are re-checked and dropped again after that. Run with:
|
||||
<p><strong><code>_iqtree.fasta</code></strong> (<code>--iqtree</code>) — alignment recoded to that same compact <code>0..k-1</code> alphabet (symbols <code>0-9A-F</code>). Under <code>--free-loss</code>, non-detection becomes <code>?</code> and columns left non-informative once missing calls are ignored are dropped first (required for <code>+ASC</code>); with <code>--iqtree-min-freq</code> also set (the default), any state rarer than that threshold is folded into the same <code>?</code> treatment, and non-informative columns are re-checked and dropped again. Run with:
|
||||
<div class="highlight"><pre><span></span><code>iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
|
||||
</code></pre></div></p>
|
||||
<p><strong><code>_iqtree_states.csv</code></strong> (<code>--iqtree</code>) — one row per state actually kept in <code>_iqtree.model</code>/<code>_iqtree.fasta</code> (header <code>iqtree_symbol,canonical_symbol,frequency</code>): <code>iqtree_symbol</code> is the compact <code>0-9A-F</code> symbol as written in those two files, <code>canonical_symbol</code> is the matching <code>_sankoff_matrix.csv</code> state, <code>frequency</code> is that state's empirical frequency at full precision (<code>_iqtree.model</code>'s own frequency line is rounded to 6 decimals). Under <code>--free-loss</code>, absent (<code>0</code>/<code>?</code>) is never a kept state, so it never appears here — nor does any state <code>--iqtree-min-freq</code> folded away for being too rare. Use this file to identify which real state a given row/column of <code>_iqtree.model</code>'s matrix corresponds to — e.g. to check whether a state showing zero exchangeability with everything else is expected (a state combination that never co-occurs with anything else in this data) or worth investigating further.</p>
|
||||
<p><strong><code>_iqtree_states.csv</code></strong> (<code>--iqtree</code>) — one row per state actually kept in <code>_iqtree.model</code>/<code>_iqtree.fasta</code> (header <code>iqtree_symbol,canonical_symbol,frequency</code>): <code>iqtree_symbol</code> is the compact <code>0-9A-F</code> symbol as written in those two files, <code>canonical_symbol</code> is the matching <code>_sankoff_matrix.csv</code> state, <code>frequency</code> is that state's empirical frequency at full precision. Under <code>--free-loss</code>, absent (<code>0</code>/<code>?</code>) is never a kept state, so it never appears here — nor does any state <code>--iqtree-min-freq</code> folded away for being too rare.</p>
|
||||
<h3 id="rare-states-and-iqtree-min-freq">Rare states and <code>--iqtree-min-freq</code></h3>
|
||||
<p>States that combine 3 or 4 central bases at once (IUPAC <code>V</code>/<code>H</code>/<code>K</code>.../<code>N</code>) are inherently rare — and, on real data, rare enough that they can make <code>iqtree3</code> itself numerically unstable ("Numerical underflow for lh-derivative" warnings, near-degenerate likelihood optimization). They're also more likely to be assembly/detection noise than genuine, widely-shared multi-way polymorphism, the same "sampling failure, not true signal" reasoning <code>--free-loss</code> already applies to non-detection. With <code>--free-loss</code> set, <code>--iqtree-min-freq</code> (default <code>0.001</code>, i.e. one in a thousand) extends that same missing-data treatment to any state below this frequency, not just absence. Check <code>_iqtree_states.csv</code> to see exactly which states survived and at what frequency; set <code>--iqtree-min-freq 0</code> to disable this and keep every state that occurs at all (the old behavior). Has no effect without <code>--free-loss</code> — there is no missing-data symbol to fold rare states into otherwise.</p>
|
||||
<p>States that combine 3 or 4 central bases at once (IUPAC <code>V</code>/<code>H</code>/<code>K</code>.../<code>N</code>) are inherently rare, and can make <code>iqtree3</code> itself numerically unstable ("Numerical underflow for lh-derivative" warnings). With <code>--free-loss</code> set, <code>--iqtree-min-freq</code> (default <code>0.001</code>, one in a thousand) extends the missing-data treatment to any state below this frequency, not just absence. Check <code>_iqtree_states.csv</code> to see exactly which states survived and at what frequency; set <code>--iqtree-min-freq 0</code> to keep every state that occurs at all. Has no effect without <code>--free-loss</code>.</p>
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -704,34 +704,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -789,14 +761,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -784,34 +784,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -869,14 +841,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -773,34 +773,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -858,14 +830,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -9,7 +9,7 @@
|
||||
|
||||
|
||||
|
||||
<link rel="prev" href="../name-tree/">
|
||||
<link rel="prev" href="../phylo/">
|
||||
|
||||
|
||||
<link rel="next" href="../estimate/">
|
||||
@@ -702,34 +702,6 @@
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -858,14 +830,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -9,7 +9,7 @@
|
||||
|
||||
|
||||
|
||||
<link rel="prev" href="../reindex/">
|
||||
<link rel="prev" href="../convert/">
|
||||
|
||||
|
||||
<link rel="next" href="../pack/">
|
||||
@@ -704,34 +704,6 @@
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../name-tree/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
name-tree
|
||||
|
||||
|
||||
|
||||
</span>
|
||||
|
||||
|
||||
|
||||
</a>
|
||||
</li>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../unitig/" class="md-nav__link">
|
||||
|
||||
@@ -789,14 +761,14 @@
|
||||
|
||||
|
||||
<li class="md-nav__item">
|
||||
<a href="../reindex/" class="md-nav__link">
|
||||
<a href="../convert/" class="md-nav__link">
|
||||
|
||||
|
||||
|
||||
<span class="md-ellipsis">
|
||||
|
||||
|
||||
reindex
|
||||
convert
|
||||
|
||||
|
||||
|
||||
|
||||
+1
-2
@@ -46,10 +46,9 @@ nav:
|
||||
- dump: usage/dump.md
|
||||
- annotate: usage/annotate.md
|
||||
- phylo: usage/phylo.md
|
||||
- name-tree: usage/name-tree.md
|
||||
- unitig: usage/unitig.md
|
||||
- estimate: usage/estimate.md
|
||||
- reindex: usage/reindex.md
|
||||
- convert: usage/convert.md
|
||||
- utils: usage/utils.md
|
||||
- pack: usage/pack.md
|
||||
- Predicates and taxonomy paths: usage/predicates.md
|
||||
|
||||
Reference in New Issue
Block a user