Introduce compare_sparse CLI tool to verify index consistency
Adds a new binary that validates bit-level consistency between dense and sparse K-mer index representations by sampling slots across partitions and reporting discrepancies. Refactors sibling cache matrix initialization into a centralized factory method to simplify control flow and standardize error handling. Introduces diagnostic configurations, benchmark scripts, and an ignored test case to support k-mer resolution analysis.
This commit is contained in:
Eric Coissac
@@ -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
|
||||
|
||||
@@ -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.
|
||||
@@ -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
|
||||
Executable
+18
@@ -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
|
||||
}
|
||||
Generated
+10
@@ -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"
|
||||
|
||||
+1
-1
@@ -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
|
||||
|
||||
@@ -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 <orig_dense> <sparse>");
|
||||
let sparse_root = std::env::args().nth(2).expect("usage: compare_sparse <orig_dense> <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(())
|
||||
}
|
||||
@@ -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"
|
||||
@@ -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 <sparse> <orig_dense>");
|
||||
let dense_root = std::env::args().nth(2).expect("usage: compare_full <sparse> <orig_dense>");
|
||||
|
||||
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(())
|
||||
}
|
||||
@@ -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<Self> {
|
||||
if with_counts && layer_dir.join("counts").exists() {
|
||||
return Layer::<PersistentCompactIntMatrix>::open(layer_dir, mode).map(Mat::Count);
|
||||
}
|
||||
if layer_dir.join("presence").join("is_multi.prsb").exists() {
|
||||
Layer::<PersistentSparseBitMatrix>::open(layer_dir, mode).map(Mat::SparsePresence)
|
||||
} else {
|
||||
Layer::<PersistentBitMatrix>::open(layer_dir, mode).map(Mat::Presence)
|
||||
}
|
||||
}
|
||||
|
||||
fn find_slot(&self, kmer: CanonicalKmer) -> Option<usize> {
|
||||
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::<PersistentCompactIntMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::Count)
|
||||
} else if is_sparse {
|
||||
Layer::<PersistentSparseBitMatrix>::open(&layer_dir, &meta.mode).ok().map(Mat::SparsePresence)
|
||||
} else {
|
||||
Layer::<PersistentBitMatrix>::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);
|
||||
|
||||
@@ -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::<PersistentCompactIntMatrix>::open(layer_dir, &meta.mode).map_err(olm_to_ok)?)
|
||||
} else {
|
||||
Mat::Presence(Layer::<PersistentBitMatrix>::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 });
|
||||
|
||||
@@ -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
|
||||
}
|
||||
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user