feat(phylo): add --iqtree-min-freq to filter rare nucleotide states

Introduces --iqtree-min-freq (default 0.001) to treat low-frequency nucleotide states as missing data during IQ-TREE alignment generation when --free-loss is active. This triggers a recoding pass that folds rare states into the missing symbol, followed by non-informative site removal and alphabet recomputation to maintain output consistency. The change also adds Sankoff model configuration files and updates related tests and documentation.

src/obikmer/src/cmd/phylo/args.rs

@@ -267,6 +267,25 @@ pub struct PhyloArgs {
     #[arg(long)]
     pub iqtree: bool,
 
+    /// Under `--iqtree --free-loss`, also recode to `?` (the same
+    /// missing-data treatment as `-`) any state whose empirical frequency
+    /// in the alignment falls below this threshold — not just genuinely
+    /// absent calls. States encoding 3 or 4 simultaneously-observed
+    /// central bases (IUPAC `V`/`H`/`K`.../`N` for 3, `N` for 4) are rare
+    /// by construction and, on real data, land in exactly this low-frequency
+    /// range — more likely assembly/detection noise than a genuine,
+    /// widely-preserved multi-way polymorphism, the same "sampling failure,
+    /// not true signal" reasoning `--free-loss` already applies to absence.
+    /// Confirmed on real data to matter: `iqtree3`'s own "Numerical
+    /// underflow for lh-derivative" warnings and the exact-zero
+    /// exchangeability rows this was meant to fix (see
+    /// `DevDocMD/theory/evolutionary_distances.md`) both trace back to
+    /// states this thin. No effect without `--free-loss` (there is no
+    /// missing-data symbol to recode to otherwise). `<prefix>_iqtree_states.csv`
+    /// reports the frequency actually used to decide.
+    #[arg(long, default_value = "0.001")]
+    pub iqtree_min_freq: f64,
+
     /// Scale factor applied before rounding real-valued costs to the
     /// integers both `--tnt`'s smatrix/cost commands and `--phyg`'s `tcm:`
     /// matrix require. Keep this small: the total tree score is this scale

src/obikmer/src/cmd/phylo/iqtree.rs

@@ -107,6 +107,19 @@ fn drop_ascertainment_noninformative(alignment: &SnpAlignment) -> SnpAlignment {
     SnpAlignment { sequences }
 }
 
+/// Recode every occurrence of a byte in `symbols` to `-` — the same
+/// "absent" byte `drop_ascertainment_noninformative`/`compact_alphabet`
+/// already treat specially under `--free-loss` (recoded to `?` further
+/// downstream). Used by `--iqtree-min-freq` to fold rare, likely-noisy
+/// states into the missing-data treatment before a second
+/// `compact_alphabet` pass, without duplicating that treatment's logic.
+fn recode_symbols_as_absent(alignment: &SnpAlignment, symbols: &[u8]) -> SnpAlignment {
+    let sequences = alignment.sequences.iter()
+        .map(|seq| seq.iter().map(|&b| if symbols.contains(&b) { b'-' } else { b }).collect())
+        .collect();
+    SnpAlignment { sequences }
+}
+
 fn compact_alphabet(alignment: &SnpAlignment, free_loss: bool) -> CompactAlphabet {
     let iupac_to_state = state_index_table();
 
@@ -240,6 +253,7 @@ pub(super) fn write_iqtree(
     labels: &[String],
     output: &Option<PathBuf>,
     free_loss: bool,
+    min_freq: f64,
 ) {
     let filtered;
     let alignment = if free_loss {
@@ -257,7 +271,41 @@ pub(super) fn write_iqtree(
         alignment
     };
 
-    let alphabet = compact_alphabet(alignment, free_loss);
+    let mut alphabet = compact_alphabet(alignment, free_loss);
+
+    // `--iqtree-min-freq`: fold rare (likely-noisy) states into the same
+    // missing-data treatment `-` already gets under `--free-loss`, then
+    // recompute the alphabet on the further-filtered alignment — see
+    // `args.rs`'s docs on `--iqtree-min-freq` and
+    // `DevDocMD/theory/evolutionary_distances.md` for why (real data:
+    // states encoding 3-4 simultaneous central bases land in exactly this
+    // low-frequency range and correlate with `iqtree3`'s own numerical
+    // instability warnings).
+    let refiltered;
+    let alignment = if free_loss {
+        let low_freq_symbols: Vec<u8> = alphabet.compact_to_old.iter().zip(alphabet.freq.iter())
+            .filter(|&(_, &f)| f < min_freq)
+            .map(|(&old, _)| STATE_SYMBOL[old as usize] as u8)
+            .collect();
+        if low_freq_symbols.is_empty() {
+            alignment
+        } else {
+            let recoded = recode_symbols_as_absent(alignment, &low_freq_symbols);
+            let before = recoded.sequences.first().map(|s| s.len()).unwrap_or(0);
+            refiltered = drop_ascertainment_noninformative(&recoded);
+            let after = refiltered.sequences.first().map(|s| s.len()).unwrap_or(0);
+            info!(
+                "--iqtree-min-freq {min_freq}: {} rare state(s) ({}) recoded as missing, {before} → {after} sites",
+                low_freq_symbols.len(),
+                low_freq_symbols.iter().map(|&b| b as char).collect::<String>(),
+            );
+            alphabet = compact_alphabet(&refiltered, free_loss);
+            &refiltered
+        }
+    } else {
+        alignment
+    };
+
     let states_path = write_iqtree_states_csv(&alphabet, output);
     let model_path = write_iqtree_model(matrix, &alphabet, output);
     let (fasta_path, n_sites) = write_iqtree_alignment(alignment, labels, &alphabet, output, free_loss);
@@ -349,4 +397,39 @@ mod tests {
         assert_eq!(lines.next(), Some("1,G,0.5"));
         assert!(lines.next().is_none());
     }
+
+    #[test]
+    fn iqtree_min_freq_folds_rare_states_into_missing() {
+        // 20 common A/C sites (60 calls total across 3 genomes) plus one
+        // site where genome 0 carries the rare ambiguity state `M` and
+        // genome 1 carries `A` (kept informative by the first
+        // ascertainment filter: two distinct non-`-` calls) — `M` ends up
+        // at 1/62, well below the 0.05 threshold used here.
+        let mut sequences: Vec<Vec<u8>> = vec![Vec::new(); 3];
+        for i in 0..20 {
+            let (a, b, c) = if i % 2 == 0 { (b'A', b'C', b'A') } else { (b'C', b'A', b'C') };
+            sequences[0].push(a);
+            sequences[1].push(b);
+            sequences[2].push(c);
+        }
+        sequences[0].push(b'M');
+        sequences[1].push(b'A');
+        sequences[2].push(b'-');
+        let alignment = SnpAlignment { sequences };
+        let labels = vec!["g1".to_string(), "g2".to_string(), "g3".to_string()];
+        let matrix = [[0.0f64; 16]; 16];
+        let prefix = std::env::temp_dir().join(format!("obikmer_test_iqtree_minfreq_{}", std::process::id()));
+        let output = Some(prefix.clone());
+
+        write_iqtree(&matrix, &alignment, &labels, &output, true, 0.05);
+
+        let states_path = format!("{}_iqtree_states.csv", prefix.display());
+        let csv = std::fs::read_to_string(&states_path).unwrap();
+        assert!(!csv.contains(",M,"), "M (freq ~1/62) must be folded into missing under --iqtree-min-freq 0.05, got:\n{csv}");
+        assert!(csv.contains(",A,") && csv.contains(",C,"), "A/C must survive (well above threshold), got:\n{csv}");
+
+        for suffix in ["_iqtree_states.csv", "_iqtree.model", "_iqtree.fasta"] {
+            std::fs::remove_file(format!("{}{suffix}", prefix.display())).ok();
+        }
+    }
 }

src/obikmer/src/cmd/phylo/mod.rs

@@ -249,7 +249,7 @@ pub fn run(args: PhyloArgs) {
             write_sankoff_phyg(&matrix, &args.output, args.sankoff_cost_scale);
         }
         if args.iqtree {
-            write_iqtree(&matrix, &alignment, &kept_labels, &args.output, args.free_loss);
+            write_iqtree(&matrix, &alignment, &kept_labels, &args.output, args.free_loss, args.iqtree_min_freq);
         }
     }
 

src/obikphylo/src/siblings/tests.rs

@@ -15,7 +15,7 @@ use super::cardinality::CardinalityExt;
 use super::distance::DistanceExt;
 use super::entropy::ShannonEntropyExt;
 use super::entropy_annex::{EntropyAnnex, ENTROPY_ANNEX_FILE_NAME};
-use super::helpers::is_minorant;
+use super::helpers::{central_base, is_minorant};
 use super::sankoff_bundle::SankoffBundleExt;
 use super::stats::SiblingStatsExt;
 use super::subsample::EntropyBias;
@@ -457,6 +457,44 @@ fn sankoff_bundle_matches_old_separate_calls() {
     assert_eq!(bundle.alignment.sequences, expected_alignment.sequences);
 }
 
+#[test]
+fn base_pair_tally_accumulates_base_a_diagnostic() {
+    // Diagnostic for a user-reported observation on real data:
+    // `composition_transitions`'s row/column for base 'A' (index 0) was
+    // *entirely* zero, including the diagonal (`same[0]`, "A stays A"),
+    // despite 'A' appearing at ~6% frequency in the pseudo-alignment
+    // itself. g1/g2 are identical (both central base 'A'), g3 differs
+    // (central base 'C') at the same family — makes the family variable
+    // (family_size=2) while g1/g2 form an eligible (ratio 0), ordinary
+    // pair both resolving `single_form` to `Some(0)` = 'A'.
+    let dir = tempdir().unwrap();
+    let g1 = build_single_genome_index(dir.path(), "g1", b"AAAAAAAAAAA");
+    let g2 = build_single_genome_index(dir.path(), "g2", b"AAAAAAAAAAA");
+    let g3 = build_single_genome_index(dir.path(), "g3", b"AAAAACAAAAA");
+    let mut rep = Reporter::new();
+    let merged = KmerIndex::merge(
+        &dir.path().join("merged.idx"),
+        &[&g1, &g2, &g3],
+        MergeMode::Presence,
+        false,
+        false,
+        1.0,
+        &mut rep,
+    ).expect("merge");
+    merged.build_sibling_annex().expect("build_sibling_annex");
+
+    let n_genomes = merged.meta().genomes.len();
+    let exclude_mask = vec![false; n_genomes];
+    let bundle = merged.sankoff_bundle(None, None, 0.5, &exclude_mask).expect("sankoff_bundle");
+
+    assert!(
+        bundle.base_pair_tally.same[0] >= 1,
+        "base 'A' (index 0) must accumulate in BasePairTally.same for an identical, included g1/g2 \
+         pair both resolving to 'A' — got same={:?}, counts={:?}",
+        bundle.base_pair_tally.same, bundle.base_pair_tally.counts,
+    );
+}
+
 #[test]
 fn entropy_annex_builds_on_demand_and_biases_selection() {
     // Same fixture, same single non-monomorphic family, known entropy15 =
@@ -537,6 +575,132 @@ fn diag_real_index_layer_distribution() {
     println!("top10 counts: {:?}", &sorted_counts[..10.min(sorted_counts.len())]);
 }
 
+#[test]
+#[ignore]
+fn diag_real_index_base_a_representation() {
+    // Diagnostic for a user-reported observation, reproduced on two
+    // unrelated real datasets (a plant and this bacterial index):
+    // `composition_transitions`'s row/column for base 'A' is entirely
+    // zero, even though 'A' appears at real (~6-22%) frequency in the
+    // pseudo-alignment itself. A minimal synthetic reproduction
+    // (`base_pair_tally_accumulates_base_a_diagnostic`) showed the
+    // pairwise tally mechanism itself works correctly for base A in
+    // isolation — so this checks one level lower, purely structural
+    // (annex bits only, no cross-partition/genome-level resolution): does
+    // base A even show up as *present* in minorant families' own
+    // `FamilyMask`s at all, at a rate proportional to the other 3 bases?
+    // If yes, the anomaly is specific to the pairwise/`single_form`
+    // resolution stage; if `has_base[0]` is itself near-zero here, the
+    // anomaly originates earlier, in `build_sibling_annex`/`central_base`.
+    let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
+        .expect("open real index");
+    let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
+    let mut has_base = [0u64; 4];
+    let mut minorant_total = 0u64;
+    let mut minorant_central_a = 0u64;
+    for layer_dir in &layer_dirs {
+        let annex = SiblingAnnex::open(&layer_dir.join(ANNEX_FILE_NAME)).unwrap();
+        for slot in 0..annex.len() {
+            let Some(mask) = annex.get(slot) else { continue };
+            if !mask.is_minorant() {
+                continue;
+            }
+            minorant_total += 1;
+            for b in 0..4u8 {
+                if mask.has(b) {
+                    has_base[b as usize] += 1;
+                }
+            }
+        }
+        // Cross-check: for a sample of this layer's minorants, is the
+        // slot's *own* central base ever actually 'A' (bit 0)? — reads
+        // the real k-mer via the layer's MPHF, not just the mask.
+        let meta = PartitionMeta::load(layer_dir.parent().unwrap()).unwrap();
+        let mphf = MphfLayer::open(layer_dir, &meta.mode).unwrap();
+        for (order, kmer) in mphf.enumerate_kmers().take(2_000_000) {
+            let Some(mask) = annex.get(order) else { continue };
+            if mask.is_minorant() && central_base(kmer, idx.kmer_size()) == 0 {
+                minorant_central_a += 1;
+            }
+        }
+    }
+    println!(
+        "minorant_total={minorant_total} has_base(A,C,G,T)={has_base:?} minorant_central_a_sampled={minorant_central_a}"
+    );
+}
+
+#[test]
+#[ignore]
+fn diag_real_index_sankoff_bundle_base_a() {
+    // Structural check (`diag_real_index_base_a_representation`) shows
+    // base A fully, proportionally represented at the family level (~23M
+    // minorants, the *highest* of the four bases) — so the anomaly must
+    // be downstream, in `sankoff_bundle`'s actual pairwise resolution.
+    // Reproduces through the real cross-partition code path (not the
+    // minimal synthetic fixture, which showed the mechanism working in
+    // isolation), bounded by `--subsample` to stay fast.
+    let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
+        .expect("open real index");
+    let n_genomes = idx.meta().genomes.len();
+    let exclude_mask = vec![false; n_genomes];
+    let bundle = idx.sankoff_bundle(Some(1_000_000), None, 0.5, &exclude_mask).expect("sankoff_bundle");
+    println!("base_pair_tally.same={:?}", bundle.base_pair_tally.same);
+    println!("base_pair_tally.counts={:?}", bundle.base_pair_tally.counts);
+    let alignment_a_count: usize = bundle.alignment.sequences.iter()
+        .map(|seq| seq.iter().filter(|&&b| b == b'A').count())
+        .sum();
+    println!("alignment raw 'A' byte count={alignment_a_count}");
+}
+
+#[test]
+#[ignore]
+fn diag_real_index_genome_mask_for_base_a_families() {
+    // Directly inspects the raw `genome_mask` array `sankoff_bundle`'s
+    // closures see, for the first few variable families where base A is
+    // present — to check whether A ever co-occurs as `single_form` in two
+    // *different* genomes at the same family at all (the precondition for
+    // `bp_same`/`bp_counts` to ever increment for base A), rather than
+    // reasoning about it further.
+    use std::sync::Arc;
+    let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
+        .expect("open real index");
+    let n_parts = idx.n_partitions();
+    let n_genomes = idx.meta().genomes.len();
+    let with_counts = idx.meta().config.with_counts;
+    let k = idx.kmer_size();
+    let n_bits = n_parts.trailing_zeros() as usize;
+    let partition = obikpartitionner::KmerPartition::open_with_config(
+        idx.root_path(), idx.kmer_size(), idx.minimizer_size(), n_bits,
+    ).unwrap();
+    let cache = Arc::new(super::cache::PartitionCache::build(&partition, n_parts, with_counts).unwrap());
+    let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).unwrap();
+
+    let mut printed = 0usize;
+    for layer_dir in &layer_dirs {
+        super::family_scan::scan_layer_families(
+            layer_dir, n_parts, n_genomes, with_counts, k, &cache, &super::family_scan::Selection::All,
+            |_family_idx, mask, genome_mask| {
+                if printed >= 15 || !mask.has(0) || mask.family_size() < 2 {
+                    return;
+                }
+                let single_a: Vec<usize> = (0..n_genomes).filter(|&g| genome_mask[g] == 1).collect();
+                let others: Vec<(usize, u8)> = (0..n_genomes)
+                    .filter(|&g| genome_mask[g] != 0 && genome_mask[g] != 1)
+                    .map(|g| (g, genome_mask[g]))
+                    .collect();
+                println!(
+                    "family bits={:#06b} size={} single_A_genomes={:?} other_nonzero_genomes={:?}",
+                    mask.bits(), mask.family_size(), single_a, others,
+                );
+                printed += 1;
+            },
+        ).unwrap();
+        if printed >= 15 {
+            break;
+        }
+    }
+}
+
 #[test]
 #[ignore]
 fn diag_plant_index_presence_matrix_sparsity() {