diff --git a/.DS_Store b/.DS_Store index 4373600f..eeb93388 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/.gitignore b/.gitignore index 5bb83c6d..b19afb41 100644 --- a/.gitignore +++ b/.gitignore @@ -12,10 +12,12 @@ data-stress *.csv Betula_exilis--IGA-24-33 benchmark/genomes +benchmark/genomes_orig benchmark/simulated_data benchmark/specimen_index_presence benchmark/specimen_index_count benchmark/global_index_presence +benchmark/global_index_presence_orig benchmark/global_index_presence_sav benchmark/all_specific benchmark/global_index_count diff --git a/Le_bug_des_A.md b/Le_bug_des_A.md new file mode 100644 index 00000000..3f25b9c2 --- /dev/null +++ b/Le_bug_des_A.md @@ -0,0 +1,43 @@ +Voici la version corrigée : + +--- + +**Bug** : dans `base_pair_tally`, toutes les transitions/comptes depuis/vers A valent 0 dans `_sankoff_params.yaml`, alors que C/G/T sont corrects. + +**Contexte** : obikmer, pipeline phylogénétique `--sankoff`. L’index est construit sur 20 génomes bactériens. Même symptôme sur un jeu de 100 génomes de plantes : A est toujours à 0. + +**Fichier clé** : `src/obikphylo/src/siblings/sankoff_bundle.rs` (Pass A + Pass B). + +**Ce qui a été vérifié** : +- Le fichier de sortie `_sankoff_params.yaml` montre bien `composition_transitions` avec A à 0 partout. +- L’index contient bien des familles avec A (`mask.has(0) == true`), et même des familles où A co-existe avec d’autres bases (`mask == 0b0011` par ex.). +- Un k-mer propriétaire de famille avec `mask == 0b0001` (A seul) a été identifié : forward `GAACAAGAGATCTCGATCTTGTCTACAAGGA`, revcomp `TCCTTGTAGACAAGATCGAGATCTCTTGTTC`. +- Le diagnostic CLI sur l’index réel donne : + - Pass A : `a_pairs=623342 a_snp=623342 a_shared=0 a_both_a=0` + - Pass B : `families_with_a=22965521 a_single_form_genomes=22913238 a_included_pairs=0 a_same_incremented=0 bp_same=[0, 96389, 222720, 277909] bp_counts[0]=[0, 0, 0, 0]` + +**Interprétation** : A est fréquemment en `single_form` (mask == 1) chez certains génomes, mais **jamais simultanément** chez deux génomes différents dans la même famille. Donc toutes les paires “avec A” sont 100% SNP → ratio = 1.0 > `ratio_ceiling=0.5` → toutes exclues par le filtre `included`. C’est pourquoi `bp_same[0]` et `bp_counts[0][*]` restent à 0. + +**Point crucial** : le bug n’apparaît **que sur l’index compacté sparse**. Sur le même index avant compaction (matrice dense `matrix.pbmx`), `--sankoff` produit des tallies corrects pour A. Dès qu’on compacte avec `pack --sparse`, A disparaît. + +**Vérifications supplémentaires (diagnostic sparse)** : +- La compaction `pack --sparse` produit une matrice `PersistentSparseBitMatrix` dont le contenu est **strictement identique** à la matrice dense d'origine : vérification exhaustive coordonnée par coordonnée sur **1 804 774 880 cellules** (512 partitions × 2 layers), **zéro différence**. +- `fill_row` et `fill_sub_matrix` (les deux chemins de lecture utilisés par le pipeline phylogénétique) restituent les mêmes bits sur dense et sparse. +- **Conclusion** : le bug n'est **pas** dans la compaction sparse elle-même, ni dans les chemins de lecture individuels. La structure stocke correctement A, C, G, T. + +**Conséquence logique** : +Si les matrices sont identiques mais que le résultat final diffère, le bug se situe dans l'**intersection** des informations — c'est-à-dire dans le code qui **combine** les lectures des deux matrices (ou qui transforme les résultats bruts en tallies). Deux endroits possibles : +1. **Le scan `sankoff_bundle`** (`family_scan.rs` + `sankoff_bundle.rs`) : la boucle qui lit les matrices, construit `genome_mask`, et accumule `bp_counts` / `same`. C'est l'étape d'intersection proprement dite. +2. **La conversion des tallies en YAML** (`obikmer/src/cmd/phylo/sankoff.rs`) : moins probable, mais possible si quelque chose sélectionne/filtre les transitions avant écriture. + +**Hypothèse la plus probable** : bug dans la résolution cross-partition lors de la construction de l'annex sibling (`build_sibling_annex`). A (bit 0) serait systématiquement manquant ou mal résolu quand on interroge les variants d'une famille depuis une partition différente. À vérifier dans `src/obikphylo/src/siblings/build.rs` et `src/obikphylo/src/siblings/cache.rs` (`PartitionCache::find` / `find_presence_batch`). + +**Prochaine étape logique** : +1. Inspecter `build_sibling_annex` pour voir si les variants avec base A sont bien générés et bien recherchés dans `cache.find`. +2. Vérifier `PartitionCache::find` et `resolve_layer_hits` pour un éventuel biais contre le bit 0. +3. Si besoin, ajouter un diagnostic ciblé (compteurs par base) **uniquement** dans `cache.rs` ou `build.rs`, pas dans `sankoff_bundle.rs` qui est déjà propre. + +**Contraintes** : +- Ne pas modifier `sankoff_bundle.rs` davantage. +- Ne pas toucher à git. +- Faire des diagnostics minimaux et ciblés. diff --git a/benchmark/test_transition_sankoff_params.yaml b/benchmark/test_transition_sankoff_params.yaml new file mode 100644 index 00000000..721e01d7 --- /dev/null +++ b/benchmark/test_transition_sankoff_params.yaml @@ -0,0 +1,167 @@ +ratio_ceiling: 0.5 +cardinality_transitions: +- from: 0 + to: 0 + count: 115955295 + probability: 0.8297007290571883 +- from: 0 + to: 1 + count: 12901663 + probability: 0.0923159153460836 +- from: 0 + to: 2 + count: 10519367 + probability: 0.07526975347801175 +- from: 0 + to: 3 + count: 339782 + probability: 0.0024312591600108434 +- from: 0 + to: 4 + count: 39459 + probability: 0.00028234295870548727 +- from: 1 + to: 0 + count: 12901663 + probability: 0.9007741539906029 +- from: 1 + to: 1 + count: 973491 + probability: 0.0679676357956696 +- from: 1 + to: 2 + count: 444884 + probability: 0.031061112720426456 +- from: 1 + to: 3 + count: 2739 + probability: 0.00019123274323474897 +- from: 1 + to: 4 + count: 84 + probability: 5.864750066344985e-6 +- from: 2 + to: 0 + count: 10519367 + probability: 0.9590597257562327 +- from: 2 + to: 1 + count: 444884 + probability: 0.04056045644508228 +- from: 2 + to: 2 + count: 3796 + probability: 0.00034608458084699006 +- from: 2 + to: 3 + count: 327 + probability: 0.000029812870900149037 +- from: 2 + to: 4 + count: 43 + probability: 3.920346937940087e-6 +- from: 3 + to: 0 + count: 339782 + probability: 0.9908029486551426 +- from: 3 + to: 1 + count: 2739 + probability: 0.007986913010007698 +- from: 3 + to: 2 + count: 327 + probability: 0.0009535306879417734 +- from: 3 + to: 3 + count: 58 + probability: 0.00016912776728019222 +- from: 3 + to: 4 + count: 30 + probability: 0.00008747987962768563 +- from: 4 + to: 0 + count: 39459 + probability: 0.9957604663487016 +- from: 4 + to: 1 + count: 84 + probability: 0.0021197668256491787 +- from: 4 + to: 2 + count: 43 + probability: 0.0010851187321775557 +- from: 4 + to: 3 + count: 30 + probability: 0.0007570595805889924 +- from: 4 + to: 4 + count: 11 + probability: 0.0002775885128826305 +composition_transitions: +- from: 'A' + to: 'A' + count: 169043 + probability: 0.49442957633191476 +- from: 'A' + to: 'C' + count: 27672 + probability: 0.08093712982055894 +- from: 'A' + to: 'G' + count: 124160 + probability: 0.36315242983957063 +- from: 'A' + to: 'T' + count: 21020 + probability: 0.06148086400795566 +- from: 'C' + to: 'A' + count: 27672 + probability: 0.08651690662664728 +- from: 'C' + to: 'C' + count: 145183 + probability: 0.4539167409213838 +- from: 'C' + to: 'G' + count: 23425 + probability: 0.07323859994684925 +- from: 'C' + to: 'T' + count: 123565 + probability: 0.3863277525051197 +- from: 'G' + to: 'A' + count: 124160 + probability: 0.3846082361177367 +- from: 'G' + to: 'C' + count: 23425 + probability: 0.07256320820761906 +- from: 'G' + to: 'G' + count: 148125 + probability: 0.4588441927749658 +- from: 'G' + to: 'T' + count: 27112 + probability: 0.08398436289967846 +- from: 'T' + to: 'A' + count: 21020 + probability: 0.06258131551760583 +- from: 'T' + to: 'C' + count: 123565 + probability: 0.3678810776371534 +- from: 'T' + to: 'G' + count: 27112 + probability: 0.08071858355439246 +- from: 'T' + to: 'T' + count: 164186 + probability: 0.4888190232908483 diff --git a/benchmark/transition.awk b/benchmark/transition.awk new file mode 100755 index 00000000..93fcbd79 --- /dev/null +++ b/benchmark/transition.awk @@ -0,0 +1,18 @@ +#!/usr/bin/awk -f + +/^>/ { + print + next +} + +{ + gsub(/[Aa]/, "a") + gsub(/[Gg]/, "A") + gsub(/a/, "G") + + gsub(/[Cc]/, "c") + gsub(/[Tt]/, "C") + gsub(/c/, "T") + + print +} diff --git a/src/Cargo.lock b/src/Cargo.lock index 1c8798ad..d3a61d54 100644 --- a/src/Cargo.lock +++ b/src/Cargo.lock @@ -493,6 +493,16 @@ dependencies = [ "impl-tools", ] +[[package]] +name = "compare_sparse" +version = "0.1.0" +dependencies = [ + "anyhow", + "fastrand", + "obicompactvec", + "obikindex", +] + [[package]] name = "console" version = "0.15.11" diff --git a/src/Cargo.toml b/src/Cargo.toml index f34f6918..84935806 100644 --- a/src/Cargo.toml +++ b/src/Cargo.toml @@ -1,5 +1,5 @@ [workspace] resolver = "3" -members = ["obikseq", "obiread", "obiskbuilder", "obifastwrite", "obikmer","obikrope","obipipeline", "obikpartitionner","obiskio","obidebruinj","obilayeredmap", "obicompactvec", "obisys", "obikindex", "obitaxonomy", "obikentropy", "obikphylo"] +members = ["obikseq", "obiread", "obiskbuilder", "obifastwrite", "obikmer","obikrope","obipipeline", "obikpartitionner","obiskio","obidebruinj","obilayeredmap", "obicompactvec", "obisys", "obikindex", "obitaxonomy", "obikentropy", "obikphylo", "compare_sparse"] [profile.release] debug = 1 diff --git a/src/bin/compare_sparse.rs b/src/bin/compare_sparse.rs new file mode 100644 index 00000000..372116b0 --- /dev/null +++ b/src/bin/compare_sparse.rs @@ -0,0 +1,130 @@ +use std::path::PathBuf; + +use obicompactvec::{BinaryMatrix, PersistentBitMatrix, PersistentSparseBitMatrix}; +use obikindex::KmerIndex; + +fn main() -> anyhow::Result<()> { + let orig_root = std::env::args().nth(1).expect("usage: compare_sparse "); + let sparse_root = std::env::args().nth(2).expect("usage: compare_sparse "); + + let orig = KmerIndex::open(&orig_root)?; + let sparse = KmerIndex::open(&sparse_root)?; + + let n_parts = orig.n_partitions(); + let n_layers = orig.n_layers_per_partition()?; + let n_genomes = orig.meta().genomes.len(); + + println!("index orig: {} partitions, {} layers, {} genomes", n_parts, n_layers, n_genomes); + println!("index sparse:{} partitions, {} layers, {} genomes", sparse.n_partitions(), sparse.n_layers_per_partition()?, sparse.meta().genomes.len()); + + let mut total_layers = 0usize; + let mut mismatches = 0usize; + let mut slot_checked = 0usize; + + for part in 0..n_parts { + let part_dir_orig = orig.partition().part_dir(part); + let part_dir_sparse = sparse.partition().part_dir(part); + + if !part_dir_orig.join("index").exists() || !part_dir_sparse.join("index").exists() { + continue; + } + + for layer in 0..n_layers { + let layer_dir_orig = part_dir_orig.join("index").join(format!("layer_{layer}")); + let layer_dir_sparse = part_dir_sparse.join("index").join(format!("layer_{layer}")); + + if !layer_dir_orig.exists() || !layer_dir_sparse.exists() { + continue; + } + + let presence_orig = layer_dir_orig.join("presence"); + let presence_sparse = layer_dir_sparse.join("presence"); + + if !presence_orig.exists() || !presence_sparse.exists() { + continue; + } + + // Ouvrir la matrice dense (PackedBitMatrix) depuis l'original + let dense = match PersistentBitMatrix::open(&presence_orig) { + Ok(m) => m, + Err(e) => { + eprintln!("ERREUR ouverture dense part={part} layer={layer}: {e}"); + continue; + } + }; + + // Ouvrir la matrice sparse + let sp = match PersistentSparseBitMatrix::open(&presence_sparse) { + Ok(m) => m, + Err(e) => { + eprintln!("ERREUR ouverture sparse part={part} layer={layer}: {e}"); + continue; + } + }; + + if dense.n_cols() != sp.n_cols() { + eprintln!("ERREUR n_cols différent part={part} layer={layer}: dense={} sparse={}", dense.n_cols(), sp.n_cols()); + continue; + } + + let n = dense.n(); + if n != sp.n() { + eprintln!("ERREUR n différent part={part} layer={layer}: dense={} sparse={}", n, sp.n()); + continue; + } + + if n == 0 { + continue; + } + + // Échantillon de 100 slots aléatoires par layer + let mut rng = fastrand::Rng::with_seed(part as u64 * 1000 + layer as u64); + let sample_size = 100.min(n); + let mut checked_any = false; + + for _ in 0..sample_size { + let slot = rng.usize(0..n); + let mut dense_row = vec![0u32; dense.n_cols()]; + let mut sparse_row = vec![0u32; sp.n_cols()]; + + dense.fill_row(slot, &mut dense_row); + sp.fill_row(slot, &mut sparse_row); + + if dense_row != sparse_row { + // Premier mismatch pour ce layer : diagnostiquer la colonne A + if !checked_any { + let a_col = 0; // A est la colonne 0 + let mut a_dense = 0u32; + let mut a_sparse = 0u32; + // scanner tous les slots pour compter les différences sur la colonne A + let mut diff_a = 0usize; + for s in 0..n { + let mut d = vec![0u32; dense.n_cols()]; + let mut s2 = vec![0u32; sp.n_cols()]; + dense.fill_row(s, &mut d); + sp.fill_row(s, &mut s2); + if d[a_col] != s2[a_col] { + diff_a += 1; + } + } + eprintln!("MISMATCH part={part} layer={layer} slot={slot}: colonne A différente sur {diff_a}/{n} slots"); + mismatches += 1; + checked_any = true; + } + } + slot_checked += 1; + } + + total_layers += 1; + } + } + + println!("Vérifié {} layers, {} slots échantillonnés", total_layers, slot_checked); + if mismatches == 0 { + println!("OK : aucune différence détectée."); + } else { + println!("MISMATCHES détectés dans {} layers", mismatches); + } + + Ok(()) +} diff --git a/src/compare_sparse/Cargo.toml b/src/compare_sparse/Cargo.toml new file mode 100644 index 00000000..f11f19d0 --- /dev/null +++ b/src/compare_sparse/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "compare_sparse" +version = "0.1.0" +edition = "2021" + +[[bin]] +name = "compare_sparse" +path = "src/main.rs" + +[dependencies] +obikindex = { path = "../obikindex" } +obicompactvec = { path = "../obicompactvec" } +anyhow = "1" +fastrand = "2" diff --git a/src/compare_sparse/src/main.rs b/src/compare_sparse/src/main.rs new file mode 100644 index 00000000..37b291eb --- /dev/null +++ b/src/compare_sparse/src/main.rs @@ -0,0 +1,100 @@ +use obicompactvec::{PersistentBitMatrix, PersistentSparseBitMatrix}; +use obikindex::KmerIndex; + +fn main() -> anyhow::Result<()> { + let sparse_root = std::env::args().nth(1).expect("usage: compare_full "); + let dense_root = std::env::args().nth(2).expect("usage: compare_full "); + + let sparse = KmerIndex::open(&sparse_root)?; + let dense = KmerIndex::open(&dense_root)?; + + let n_parts = sparse.n_partitions(); + let n_layers = sparse.n_layers_per_partition()?; + + println!("index sparse: {} partitions, {} layers", sparse.n_partitions(), n_layers); + println!("index dense: {} partitions, {} layers", dense.n_partitions(), dense.n_layers_per_partition()?); + + let mut total_layers = 0usize; + let mut total_cells = 0usize; + let mut mismatched_cells = 0usize; + let mut first_mismatch = None; + + for part in 0..n_parts { + let part_dir_sparse = sparse.partition().part_dir(part); + let part_dir_dense = dense.partition().part_dir(part); + + if !part_dir_sparse.join("index").exists() || !part_dir_dense.join("index").exists() { + continue; + } + + for layer in 0..n_layers { + let layer_dir_sparse = part_dir_sparse.join("index").join(format!("layer_{layer}")); + let layer_dir_dense = part_dir_dense.join("index").join(format!("layer_{layer}")); + + if !layer_dir_sparse.exists() || !layer_dir_dense.exists() { + continue; + } + + let presence_sparse = layer_dir_sparse.join("presence"); + let presence_dense = layer_dir_dense.join("presence"); + + if !presence_sparse.exists() || !presence_dense.exists() { + continue; + } + + let sp = match PersistentSparseBitMatrix::open(&presence_sparse) { + Ok(m) => m, + Err(e) => { + eprintln!("ERREUR ouverture sparse part={part} layer={layer}: {e}"); + continue; + } + }; + + let dn = match PersistentBitMatrix::open(&layer_dir_dense) { + Ok(m) => m, + Err(e) => { + eprintln!("ERREUR ouverture dense part={part} layer={layer}: {e}"); + continue; + } + }; + + if dn.n_cols() != sp.n_cols() || dn.n() != sp.n() || dn.n() == 0 { + continue; + } + + let n = dn.n(); + let n_cols = dn.n_cols(); + + // Vérification exhaustive : tous les slots, toutes les colonnes + let mut dense_row = vec![0u32; n_cols]; + let mut sparse_row = vec![0u32; n_cols]; + + for slot in 0..n { + dn.fill_row(slot, &mut dense_row); + sp.fill_row(slot, &mut sparse_row); + + for c in 0..n_cols { + total_cells += 1; + if dense_row[c] != sparse_row[c] { + mismatched_cells += 1; + if first_mismatch.is_none() { + first_mismatch = Some((part, layer, slot, c, dense_row[c], sparse_row[c])); + } + } + } + } + + total_layers += 1; + } + } + + println!("Vérifié {} layers, {} cellules totales", total_layers, total_cells); + if let Some((part, layer, slot, col, d, s)) = first_mismatch { + eprintln!("PREMIER MISMATCH : part={part} layer={layer} slot={slot} col={col} dense={d} sparse={s}"); + println!("MISMATCHES : {} cellules différentes", mismatched_cells); + } else { + println!("OK : toutes les {} cellules sont identiques entre dense et sparse.", total_cells); + } + + Ok(()) +} diff --git a/src/obikphylo/src/siblings/cache.rs b/src/obikphylo/src/siblings/cache.rs index 0cadf0e6..3c729b18 100644 --- a/src/obikphylo/src/siblings/cache.rs +++ b/src/obikphylo/src/siblings/cache.rs @@ -1,10 +1,12 @@ use rayon::prelude::*; +use std::path::Path; + use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix, PersistentSparseBitMatrix}; use obikpartitionner::KmerPartition; use obikseq::CanonicalKmer; -use obilayeredmap::Layer; -use obilayeredmap::meta::PartitionMeta; +use obilayeredmap::{Layer, OLMResult}; +use obilayeredmap::meta::{IndexMode, PartitionMeta}; use obisys::progress_bar; use obikindex::OKIResult; @@ -50,6 +52,31 @@ pub(super) enum Mat { } impl Mat { + /// Open one layer's matrix, auto-detecting count vs. dense-presence vs. + /// sparse-presence from what is actually on disk — the single source of + /// truth for *every* caller that opens a layer's own matrix (this + /// module's cross-partition [`PartitionCache::build`] and + /// `family_scan::scan_layer_families`'s own-layer lookup alike), so the + /// two can never disagree about which format a layer's `presence/` was + /// packed to. Before this existed, `scan_layer_families` open-coded its + /// own (non-sparse-aware) copy of this decision: on a `pack --sparse`d + /// layer, `presence/matrix.pbmx` no longer exists (removed by + /// `pack_sparse_bit_matrix`), so a caller that only checks for it and + /// otherwise unconditionally opens `PersistentBitMatrix` doesn't error — + /// `PersistentBitMatrix::open`'s own auto-detection falls all the way + /// through to the `Implicit` (mono-genome, all-present) case, silently + /// corrupting every read. + pub(super) fn open(layer_dir: &Path, mode: &IndexMode, with_counts: bool) -> OLMResult { + if with_counts && layer_dir.join("counts").exists() { + return Layer::::open(layer_dir, mode).map(Mat::Count); + } + if layer_dir.join("presence").join("is_multi.prsb").exists() { + Layer::::open(layer_dir, mode).map(Mat::SparsePresence) + } else { + Layer::::open(layer_dir, mode).map(Mat::Presence) + } + } + fn find_slot(&self, kmer: CanonicalKmer) -> Option { match self { Mat::Count(l) => l.find_slot(kmer), @@ -156,21 +183,7 @@ impl PartitionCache { let mut mats = Vec::with_capacity(meta.n_layers); for l in 0..meta.n_layers { let layer_dir = index_dir.join(format!("layer_{l}")); - let use_counts = with_counts && layer_dir.join("counts").exists(); - // `pack --sparse` converts a layer's `presence/` directory - // in place, leaving `is_multi.prsb` as the on-disk marker - // that distinguishes the sparse format - // (`obicompactvec::PersistentSparseBitMatrix`) from the - // dense one — see `bitmatrix/sparse.rs`'s module docs. - let is_sparse = layer_dir.join("presence").join("is_multi.prsb").exists(); - let mat = if use_counts { - Layer::::open(&layer_dir, &meta.mode).ok().map(Mat::Count) - } else if is_sparse { - Layer::::open(&layer_dir, &meta.mode).ok().map(Mat::SparsePresence) - } else { - Layer::::open(&layer_dir, &meta.mode).ok().map(Mat::Presence) - }; - let Some(mat) = mat else { continue }; + let Ok(mat) = Mat::open(&layer_dir, &meta.mode, with_counts) else { continue }; mats.push(mat); } pb.inc(1); diff --git a/src/obikphylo/src/siblings/family_scan.rs b/src/obikphylo/src/siblings/family_scan.rs index d76cfacc..709b28c2 100644 --- a/src/obikphylo/src/siblings/family_scan.rs +++ b/src/obikphylo/src/siblings/family_scan.rs @@ -53,9 +53,7 @@ use std::sync::atomic::{AtomicU8, Ordering}; use rayon::prelude::*; -use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix}; use obikseq::CanonicalKmer; -use obilayeredmap::Layer; use obilayeredmap::meta::PartitionMeta; use obipipeline::{ThrottleGuard, throttle}; @@ -195,12 +193,7 @@ pub(super) fn scan_layer_families( let meta = PartitionMeta::load(index_dir).map_err(olm_to_ok)?; let annex = Arc::new(SiblingAnnex::open(&layer_dir.join(ANNEX_FILE_NAME))?); - let use_counts = with_counts && layer_dir.join("counts").exists(); - let mat = if use_counts { - Mat::Count(Layer::::open(layer_dir, &meta.mode).map_err(olm_to_ok)?) - } else { - Mat::Presence(Layer::::open(layer_dir, &meta.mode).map_err(olm_to_ok)?) - }; + let mat = Mat::open(layer_dir, &meta.mode, with_counts).map_err(olm_to_ok)?; let n_cols = mat.n_cols().min(n_genomes); let ctx = Arc::new(LayerCtx { mat, n_parts, n_genomes, n_cols, k }); diff --git a/src/obikphylo/src/siblings/tests.rs b/src/obikphylo/src/siblings/tests.rs index f984a376..b21b21c2 100644 --- a/src/obikphylo/src/siblings/tests.rs +++ b/src/obikphylo/src/siblings/tests.rs @@ -1089,3 +1089,96 @@ fn bench_real_persistent_sparse_bit_matrix_on_disk_size() { let _ = std::fs::remove_dir_all(&out_dir); } + +#[test] +#[ignore] +fn diag_real_index_verify_kmer_resolution() { + 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()); + + // The k-mer with A+other (mask 0b0011 = A+C) + let kmer_str = "AGCTAGCTATTGAGCCTGGTCCGTATGAAAC"; + let kmer = CanonicalKmer::from_str(kmer_str, k).unwrap(); + let rc = kmer.revcomp(); + let rc_str = kmer_to_string(&rc, k); + + println!("kmer_forward={} kmer_revcomp={}", kmer_str, rc_str); + println!("kmer partition={}", kmer.partition(n_parts)); + + // Check if this k-mer is found in its own partition + let dest_part = kmer.partition(n_parts); + if let Some(layer) = cache.find(dest_part, kmer) { + println!("FOUND in partition {} layer {}", dest_part, layer); + } else { + println!("NOT FOUND in partition {}", dest_part); + } + + // Now check the variant with A at central position + // The mask is 0b0011 = A+C, so there should be a variant with A + // Let's generate the canonical neighbors and check which one has A + let mut found_variants = Vec::new(); + for variant in kmer.central_canonical_neighbors() { + if variant == kmer { + continue; + } + let base = central_base(variant, k); + if base == 0 { // A + let v_part = variant.partition(n_parts); + println!("variant with A: partition={} found={}", v_part, cache.find(v_part, variant).is_some()); + found_variants.push(variant); + } + } + + // Also check the C variant + for variant in kmer.central_canonical_neighbors() { + if variant == kmer { + continue; + } + let base = central_base(variant, k); + if base == 1 { // C + let v_part = variant.partition(n_parts); + println!("variant with C: partition={} found={}", v_part, cache.find(v_part, variant).is_some()); + } + } +} + +fn kmer_to_string(kmer: &CanonicalKmer, k: usize) -> String { + let mut s = String::with_capacity(k); + for i in 0..k { + let n = kmer.nucleotide(i); + s.push(match n { + 0 => 'A', + 1 => 'C', + 2 => 'G', + 3 => 'T', + _ => unreachable!(), + }); + } + s +} + +fn kmer_revcomp_to_string(kmer: &CanonicalKmer, k: usize) -> String { + let rc = kmer.revcomp(); + let mut s = String::with_capacity(k); + for i in 0..k { + let n = rc.nucleotide(i); + s.push(match n { + 0 => 'A', + 1 => 'C', + 2 => 'G', + 3 => 'T', + _ => unreachable!(), + }); + } + s +} diff --git a/src/sankoff_params.yaml b/src/sankoff_params.yaml new file mode 100644 index 00000000..1e488fab --- /dev/null +++ b/src/sankoff_params.yaml @@ -0,0 +1,167 @@ +ratio_ceiling: 0.5 +cardinality_transitions: +- from: 0 + to: 0 + count: 122014779 + probability: 0.870049812090934 +- from: 0 + to: 1 + count: 17966272 + probability: 0.12811195254940888 +- from: 0 + to: 2 + count: 233683 + probability: 0.001666321505518981 +- from: 0 + to: 3 + count: 24109 + probability: 0.00017191385413811494 +- from: 0 + to: 4 + count: 0 + probability: 0.0 +- from: 1 + to: 0 + count: 17966272 + probability: 0.967023782579572 +- from: 1 + to: 1 + count: 602798 + probability: 0.03244523972983381 +- from: 1 + to: 2 + count: 9562 + probability: 0.0005146688978673966 +- from: 1 + to: 3 + count: 303 + probability: 0.00001630879272681669 +- from: 1 + to: 4 + count: 0 + probability: 0.0 +- from: 2 + to: 0 + count: 233683 + probability: 0.9585500516842502 +- from: 2 + to: 1 + count: 9562 + probability: 0.039222603245442765 +- from: 2 + to: 2 + count: 438 + probability: 0.0017966429848885097 +- from: 2 + to: 3 + count: 105 + probability: 0.00043070208541847837 +- from: 2 + to: 4 + count: 0 + probability: 0.0 +- from: 3 + to: 0 + count: 24109 + probability: 0.9827171564831044 +- from: 3 + to: 1 + count: 303 + probability: 0.012350711286838137 +- from: 3 + to: 2 + count: 105 + probability: 0.004279949455834998 +- from: 3 + to: 3 + count: 16 + probability: 0.0006521827742224758 +- from: 3 + to: 4 + count: 0 + probability: 0.0 +- from: 4 + to: 0 + count: 0 + probability: 0.0 +- from: 4 + to: 1 + count: 0 + probability: 0.0 +- from: 4 + to: 2 + count: 0 + probability: 0.0 +- from: 4 + to: 3 + count: 0 + probability: 0.0 +- from: 4 + to: 4 + count: 0 + probability: 0.0 +composition_transitions: +- from: 'A' + to: 'A' + count: 0 + probability: 0.0 +- from: 'A' + to: 'C' + count: 0 + probability: 0.0 +- from: 'A' + to: 'G' + count: 0 + probability: 0.0 +- from: 'A' + to: 'T' + count: 0 + probability: 0.0 +- from: 'C' + to: 'A' + count: 0 + probability: 0.0 +- from: 'C' + to: 'C' + count: 96389 + probability: 0.9663832688335907 +- from: 'C' + to: 'G' + count: 1290 + probability: 0.012933368089671353 +- from: 'C' + to: 'T' + count: 2063 + probability: 0.020683363076737984 +- from: 'G' + to: 'A' + count: 0 + probability: 0.0 +- from: 'G' + to: 'C' + count: 1290 + probability: 0.005696948820201646 +- from: 'G' + to: 'G' + count: 222720 + probability: 0.9835848381669073 +- from: 'G' + to: 'T' + count: 2427 + probability: 0.010718213012891003 +- from: 'T' + to: 'A' + count: 0 + probability: 0.0 +- from: 'T' + to: 'C' + count: 2063 + probability: 0.007305266661709142 +- from: 'T' + to: 'G' + count: 2427 + probability: 0.008594223067362136 +- from: 'T' + to: 'T' + count: 277909 + probability: 0.9841005102709287