From 135ffd2e41c2d340ae7761ac478a73fc82b6c9ac Mon Sep 17 00:00:00 2001 From: YPARK Date: Thu, 1 Oct 2026 18:14:13 -0700 Subject: [PATCH] Locus matching: tagged loci, composite rows, exact preference - The gene rule strips a feature-type tag from a coordinate and keeps the rest whole: chr1:0-100_Peaks matches chr1:0-100 again; Mixed keys it as 1:0-100. - GeneIndex finds a row's locus part (whole name, '/'-core, or the locus of an id_name join such as chr1:100-200_GENE1). A locus query matches the same name first, then locus key plus rest, then a bare key; a whole-locus row wins its key over one that only starts with it. - h5ad feature types are read from var/feature_type or var/feature_types (scanpy, muon). - The non-colon warning fires at a tenth of the axis, and also when an explicit overlap kind finds no colon loci; an untyped list that only partly reads as intervals is logged; the import count is rows. Requires legume-genomic-types 0.5.3. Release 0.7.2. --- Cargo.lock | 6 +- Cargo.toml | 4 +- src/aux/feature_names.rs | 57 +++++++++---- src/handlers/builders/from_h5ad.rs | 8 +- src/utilities/name_matching.rs | 121 ++++++++++++++++++++------- src/utilities/name_matching/tests.rs | 24 ++++++ 6 files changed, 166 insertions(+), 54 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index b07de6e..f9375e9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -930,7 +930,7 @@ dependencies = [ [[package]] name = "data-beans" -version = "0.7.1" +version = "0.7.2" dependencies = [ "anyhow", "approx", @@ -2058,9 +2058,9 @@ checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" [[package]] name = "legume-genomic-types" -version = "0.5.2" +version = "0.5.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fa51741a411361db709efda8c16f3dea273f19b3afce10241247a3234f2b3c97" +checksum = "a1abfb0da9e371a7c9cb679d5d062f6784bf3115824c76613a95a0a12a319df0" dependencies = [ "anyhow", "dashmap", diff --git a/Cargo.toml b/Cargo.toml index 2affd33..526862e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,7 +9,7 @@ edition = "2021" readme = "README.md" # Own version (decoupled from workspace.version): data-beans evolves its CLI / # QC API independently of the shared utility crates. -version = "0.7.1" +version = "0.7.2" rust-version = "1.91" # Ship only the sources, README and license: nothing else in the working tree # (data stores, saved name lists, scratch files) can end up on crates.io. @@ -49,7 +49,7 @@ hdf5 = ["dep:hdf5", "ndarray"] [dependencies] legume-numeric = { version = "0.8.14", default-features = false, features = ["matrix", "param"] } -legume-genomic-types = "0.5.2" +legume-genomic-types = "0.5.3" anyhow = "1.0" clap = { version = "4.5.20", features = ["derive"] } diff --git a/src/aux/feature_names.rs b/src/aux/feature_names.rs index 312c1c0..e688d72 100644 --- a/src/aux/feature_names.rs +++ b/src/aux/feature_names.rs @@ -117,19 +117,7 @@ impl FeatureNameKind { let pct_locus = n_locus as f32 / n as f32; let pct_gene = n_gene_like as f32 / n as f32; if pct_locus < 0.50 { - let n_spelled = names - .iter() - .filter(|name| { - !coordinates::is_locus(name) && coordinates::import_interval(name).is_some() - }) - .count(); - if n_spelled * 2 >= n { - log::warn!( - "{n_spelled} of {n} row names read as intervals only in a non-colon \ - spelling (e.g. `chr1-100-200`); they are not loci here. Re-import them so \ - peaks are named `chr:start-end`." - ); - } + warn_non_colon_loci(names); } if pct_locus >= 0.10 && pct_gene >= 0.10 { Self::Mixed @@ -207,7 +195,8 @@ pub use genomic_data::coordinates::locus_key; /// Per-name rule of [`FeatureNameKind::Mixed`]: locus key, else the gene /// rule for gene-style names, else the name unchanged. fn mixed_canonicalize(name: &str) -> Box { - locus_key(name).unwrap_or_else(|| gene_canonicalize(name, '_')) + let untagged = strip_feature_type_suffix(name, '_'); + locus_key(untagged).unwrap_or_else(|| gene_canonicalize(name, '_')) } /// Build the overlap-merge canonical map from a flat list of row names @@ -294,9 +283,29 @@ pub fn build_locus_overlap_canonical_map(names: &[Box]) -> HashMap out.insert(names[i].clone(), cluster.locus_key()); } } + if out.is_empty() { + warn_non_colon_loci(names); + } out } +/// Warn when many of `names` read as intervals only in a non-colon spelling +/// (`chr1-100-200`): they are not loci, so a locus rule passes them through. +fn warn_non_colon_loci(names: &[Box]) { + let n_spelled = names + .iter() + .filter(|name| !coordinates::is_locus(name) && coordinates::import_interval(name).is_some()) + .count(); + if n_spelled > 0 && n_spelled * 10 >= names.len() { + log::warn!( + "{n_spelled} of {} row names read as intervals only in a non-colon spelling \ + (e.g. `chr1-100-200`); they are not loci here. Re-import them so peaks are named \ + `chr:start-end`.", + names.len() + ); + } +} + /// Build a `RowNameCanonicalizer` for [`FeatureNameKind::LocusOverlap`]. /// `names` should be the concatenation of every input backend's row /// names (in any order). The returned canonicalizer does @@ -341,7 +350,15 @@ pub fn build_mixed_kind_canonicalizer(names: &[Box]) -> RowNameCanonicalize /// first so the actual symbol becomes the rsplit target. fn gene_canonicalize(name: &str, delim: char) -> Box { let (head, rest) = gene_part(name); - match gene_symbol(head, delim) { + // A coordinate keeps itself, without a feature-type tag + // (`chr1:0-100_Peaks` matches `chr1:0-100`). + let untagged = strip_feature_type_suffix(head, delim); + let symbol = if coordinates::is_region(untagged) { + Some(untagged) + } else { + gene_symbol(head, delim) + }; + match symbol { Some(symbol) if rest.is_empty() => symbol.into(), Some(symbol) => format!("{symbol}{rest}").into_boxed_str(), None => name.into(), @@ -745,11 +762,19 @@ mod tests { gene.canonicalize("ENSG000_GENE1/count/spliced").as_ref(), "GENE1/count/spliced" ); + // A tagged locus loses its tag and nothing else. + assert_eq!(gene.canonicalize("chr1:0-100_Peaks").as_ref(), "chr1:0-100"); + assert_eq!( + FeatureNameKind::Mixed + .canonicalize("chr1:0-100_Peaks") + .as_ref(), + "1:0-100" + ); // A locus is never cut at `_`, even when its contig name has one or // it carries a feature-type tag. assert_eq!( gene.canonicalize("chrUn_CTG1v1:0-100_Peaks").as_ref(), - "chrUn_CTG1v1:0-100_Peaks" + "chrUn_CTG1v1:0-100" ); assert_eq!( gene.canonicalize("chrUn_CTG1v1:0-100").as_ref(), diff --git a/src/handlers/builders/from_h5ad.rs b/src/handlers/builders/from_h5ad.rs index c2e4101..ad52706 100644 --- a/src/handlers/builders/from_h5ad.rs +++ b/src/handlers/builders/from_h5ad.rs @@ -211,8 +211,11 @@ pub fn run_build_from_h5ad(args: &FromH5adArgs) -> anyhow::Result<()> { assert_eq!(nrows, row_ids.len()); assert_eq!(nrows, row_names.len()); - // Feature types from var/feature_type (often categorical) - let typed = read_h5ad_column(&var_group, "feature_type").ok(); + // Feature types from var/feature_type or var/feature_types (scanpy, + // muon), often categorical + let typed = read_h5ad_column(&var_group, "feature_type") + .or_else(|_| read_h5ad_column(&var_group, "feature_types")) + .ok(); let has_types = typed.is_some(); let mut row_types: Vec> = typed.unwrap_or_else(|| vec![Box::from(""); nrows]); if nrows < row_types.len() { @@ -221,7 +224,6 @@ pub fn run_build_from_h5ad(args: &FromH5adArgs) -> anyhow::Result<()> { assert_eq!(nrows, row_types.len()); // Peak rows in chr:start-end form, then composite row names: id_name - let mut row_ids = row_ids; let n_peaks = colon_peak_names( &mut row_ids, &mut row_names, diff --git a/src/utilities/name_matching.rs b/src/utilities/name_matching.rs index b64f580..621677b 100644 --- a/src/utilities/name_matching.rs +++ b/src/utilities/name_matching.rs @@ -1,5 +1,5 @@ use crate::sparse_io::ROW_SEP; -use genomic_data::coordinates::{import_interval, locus_key}; +use genomic_data::coordinates::{import_interval, is_locus, locus_key}; use rayon::prelude::*; use rustc_hash::FxHashMap as HashMap; @@ -66,21 +66,36 @@ pub fn colon_peak_names( } None => true, }; - if types.is_none() && !ids.iter().all(|id| import_interval(id).is_some()) { - return 0; + if types.is_none() { + let n_read = ids + .iter() + .filter(|id| import_interval(id).is_some()) + .count(); + if n_read < ids.len() { + if n_read * 2 >= ids.len() { + log::info!( + "{n_read} of {} untyped rows read as intervals, but not all; \ + the names are kept as written", + ids.len() + ); + } + return 0; + } } let mut n = 0; for (i, (id, name)) in ids.iter_mut().zip(names.iter_mut()).enumerate() { if !is_peak(i) { continue; } + let mut changed = false; for s in [id, name] { if let Some(l) = import_interval(s) { let colon = l.to_string().into_boxed_str(); - n += usize::from(*s != colon); + changed |= *s != colon; *s = colon; } } + n += usize::from(changed); } n } @@ -340,43 +355,84 @@ pub struct GeneIndex { exact: HashMap, symbol: HashMap, ensg: HashMap, - /// Locus rows by their locus key. Loci match only here, case kept. + /// Rows with a locus part ([`locus_part`]), by their name as written, + /// then by locus key plus the rest of the name, then (for rows with a + /// rest) by the bare locus key. Loci match only here, case kept. + locus_raw: HashMap, usize>, locus: HashMap, usize>, } +/// A name's locus part and the rest after it: the whole name when it is a +/// locus, its `/`-core (`chr1:1-2/count/spliced`), or the locus before an +/// `id{ROW_SEP}name` join (`chr1:1-2_GENE1`). `None` when no part is a locus. +fn locus_part(name: &str) -> Option<(&str, &str)> { + if is_locus(name) { + return Some((name, "")); + } + let core = name.split('/').next().unwrap_or(name); + let cut = if is_locus(core) { + core.len() + } else { + let colon = core.find(':')?; + colon + 1 + core[colon + 1..].find(ROW_SEP)? + }; + let (part, rest) = name.split_at(cut); + is_locus(part).then_some((part, rest)) +} + +/// The key a name with a locus part is matched by: the locus key, then the +/// rest of the name as written. +fn locus_match_key(part: &str, rest: &str) -> Option> { + let key = locus_key(part)?; + Some(if rest.is_empty() { + key + } else { + format!("{key}{rest}").into_boxed_str() + }) +} + #[allow(dead_code)] // consumed by downstream crates (geu, senna), not the data-beans bin impl GeneIndex { /// Build the index from the dictionary's gene-name order. The first row /// wins on duplicate keys (matching positional-scan semantics). #[must_use] pub fn build(gene_names: &[Box]) -> Self { - // A locus row gets its key and an empty lowered name, which keeps it - // out of every gene tier, the fallback scan included. A row whose - // `/`-core is a locus (`chr1:1-2/count/spliced`) also gets its core - // key, and stays in the gene tiers for its full name. - let (lowered, keys): (Vec, Vec>>) = gene_names + // A whole-locus row gets an empty lowered name, which keeps it out + // of every gene tier, the fallback scan included. A row with a + // locus part plus a rest stays in the gene tiers for its full name. + let lowered: Vec = gene_names .par_iter() - .map(|g| match locus_key(g) { - Some(key) => (String::new(), Some(key)), - None => { - let core = g.split('/').next().unwrap_or(g); - let key = if core.len() < g.len() { - locus_key(core) - } else { - None - }; - (g.to_lowercase(), key) + .map(|g| { + if is_locus(g) { + String::new() + } else { + g.to_lowercase() } }) - .unzip(); - let mut exact: HashMap = HashMap::default(); - let mut symbol: HashMap = HashMap::default(); - let mut ensg: HashMap = HashMap::default(); + .collect(); + let mut locus_raw: HashMap, usize> = HashMap::default(); let mut locus: HashMap, usize> = HashMap::default(); - for (i, (low, key)) in lowered.iter().zip(keys).enumerate() { - if let Some(key) = key { + let mut bare: Vec<(Box, usize)> = Vec::new(); + for (i, g) in gene_names.iter().enumerate() { + let Some((part, rest)) = locus_part(g) else { + continue; + }; + locus_raw.entry(g.clone()).or_insert(i); + if let Some(key) = locus_match_key(part, rest) { locus.entry(key).or_insert(i); } + if !rest.is_empty() { + bare.extend(locus_key(part).map(|k| (k, i))); + } + } + // A whole-locus row wins its key over a row that only starts with it. + for (key, i) in bare { + locus.entry(key).or_insert(i); + } + let mut exact: HashMap = HashMap::default(); + let mut symbol: HashMap = HashMap::default(); + let mut ensg: HashMap = HashMap::default(); + for (i, low) in lowered.iter().enumerate() { if low.is_empty() { continue; } @@ -401,6 +457,7 @@ impl GeneIndex { exact, symbol, ensg, + locus_raw, locus, } } @@ -408,10 +465,14 @@ impl GeneIndex { /// Row index for `gene`, or `None` if unmatched (tiers above). #[must_use] pub fn match_gene(&self, gene: &str) -> Option { - // A locus matches a locus row by key, strictly: no case folding, - // aliasing or prefix fallback. - if let Some(key) = locus_key(gene) { - return self.locus.get(&key).copied(); + // A name with a locus part matches only a row with one: the same + // name as written first, then by locus key plus the rest. Strictly: + // no case folding, aliasing or prefix fallback. + if let Some((part, rest)) = locus_part(gene) { + if let Some(&i) = self.locus_raw.get(gene) { + return Some(i); + } + return locus_match_key(part, rest).and_then(|k| self.locus.get(&k).copied()); } let gl = gene.to_lowercase(); if let Some(&i) = self.exact.get(&gl) { diff --git a/src/utilities/name_matching/tests.rs b/src/utilities/name_matching/tests.rs index bce9777..c4e0ae4 100644 --- a/src/utilities/name_matching/tests.rs +++ b/src/utilities/name_matching/tests.rs @@ -120,3 +120,27 @@ fn an_untyped_list_converts_only_when_every_row_is_an_interval() { assert_eq!(colon_peak_names(&mut ids, &mut nm, None), 0); assert_eq!(ids, names(&["chr1-100-200", "GENE1"])); } + +#[test] +fn locus_queries_prefer_the_exact_row_and_reach_composite_rows() { + let idx = GeneIndex::build(&names(&["chr2:5-9/count/spliced", "chr2:5-9"])); + assert_eq!(idx.match_gene("chr2:5-9"), Some(1), "exact name first"); + assert_eq!( + idx.match_gene("2:5-9"), + Some(1), + "a whole-locus row wins its key" + ); + assert_eq!(idx.match_gene("chr2:5-9/count/spliced"), Some(0)); + assert_eq!( + idx.match_gene("2:5-9/count/spliced"), + Some(0), + "chr prefix free" + ); + let idx = GeneIndex::build(&names(&["chrX:0-100", "X:0-100"])); + assert_eq!(idx.match_gene("X:0-100"), Some(1)); + assert_eq!(idx.match_gene("chrX:0-100"), Some(0)); + // An id_name composite peak row is reached by its locus. + let idx = GeneIndex::build(&names(&["GENE9", "chr1:100-200_GENE1"])); + assert_eq!(idx.match_gene("chr1:100-200"), Some(1)); + assert_eq!(idx.match_gene("x:0-100/count/spliced"), None); +}