<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>
<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>
<td><spanclass="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 <spanclass="arithmatex">\(D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)\)</span> where <spanclass="arithmatex">\(J = 1 - D_{\text{jaccard}}\)</span> and <spanclass="arithmatex">\(k\)</span> is the index's kmer size; clamped to 1.0 when <spanclass="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): <spanclass="arithmatex">\(D = \sum_i \mathbb{1}[a_i \ne b_i]\)</span></td>
</tr>
<tr>
<td><code>bray-curtis</code></td>
<td><spanclass="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 <spanclass="arithmatex">\(p_i = c_i / \sum_j c_j\)</span> instead of raw counts</td>
</tr>
<tr>
<td><code>euclidean</code></td>
<td><spanclass="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><spanclass="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 <spanclass="arithmatex">\([0, 1]\)</span></td>
<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 <spanclass="arithmatex">\(L\)</span> be its total number of eligible loci (both genomes single-copy at that family), <spanclass="arithmatex">\(p\)</span> the raw proportion of substitutions among those loci, <spanclass="arithmatex">\(P\)</span>/<spanclass="arithmatex">\(Q\)</span> the transition/transversion proportions, <spanclass="arithmatex">\(Q_1\)</span>/<spanclass="arithmatex">\(Q_2\)</span> Kimura's two transversion categories (A↔C & G↔T vs. A↔T & C↔G), <spanclass="arithmatex">\(P_1\)</span>/<spanclass="arithmatex">\(P_2\)</span> the purine (A↔G) / pyrimidine (C↔T) transition proportions, and <spanclass="arithmatex">\(\pi_A,\pi_C,\pi_G,\pi_T\)</span> the pair's pooled base frequencies.</p>
<p><code>--gamma-shape ALPHA</code> applies to every value above except <code>snp-raw</code> and <code>snp-tv</code>: each <spanclass="arithmatex">\(-\ln(x)\)</span> term in the formulas above is replaced by <spanclass="arithmatex">\(\alpha\left(x^{-1/\alpha}-1\right)\)</span> (the same weight, same <spanclass="arithmatex">\(x\)</span>).</p>
<h3id="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>
<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>
<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>
<p>A family is eligible for a genome pair <spanclass="arithmatex">\((i,j)\)</span> only if both genomes carry exactly one of its observed forms (single-copy, unambiguous).</p>
<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>
<h3id="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>
<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>
<h2id="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>
<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>
<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>
<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>
<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>
<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 <spanclass="arithmatex">\(5 \times 5\)</span> transition matrix over family cardinality (0-4 observed forms) and a <spanclass="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 <spanclass="arithmatex">\(16 \times 16\)</span> transition probability matrix <spanclass="arithmatex">\(P\)</span>, converted to a symmetric cost matrix via <spanclass="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>
<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 <spanclass="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>
<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>
<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>
<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>
<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 <spanclass="arithmatex">\(-\ln P(a,b)\)</span>, 4 decimals.</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 <spanclass="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. Run with:
<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>
<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>
<scriptid="__config"type="application/json">{"annotate":null,"base":"../..","features":[],"search":"../../assets/javascripts/workers/search.2c215733.min.js","tags":null,"translations":{"clipboard.copied":"Copied to clipboard","clipboard.copy":"Copy to clipboard","search.result.more.one":"1 more on this page","search.result.more.other":"# more on this page","search.result.none":"No matching documents","search.result.one":"1 matching document","search.result.other":"# matching documents","search.result.placeholder":"Type to start searching","search.result.term.missing":"Missing","select.version":"Select version"},"version":null}</script>
