Locus matching: tagged loci, composite rows, exact preference

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"] }

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<str> {
-    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<str>]) -> HashMap<Box<str>
             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<str>]) {
+    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<str>]) -> RowNameCanonicalize
 /// first so the actual symbol becomes the rsplit target.
 fn gene_canonicalize(name: &str, delim: char) -> Box<str> {
     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(),

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<Box<str>> = 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,

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<String, usize>,
     symbol: HashMap<String, usize>,
     ensg: HashMap<String, usize>,
-    /// 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<Box<str>, usize>,
     locus: HashMap<Box<str>, 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<Box<str>> {
+    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<str>]) -> 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<String>, Vec<Option<Box<str>>>) = 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<String> = 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<String, usize> = HashMap::default();
-        let mut symbol: HashMap<String, usize> = HashMap::default();
-        let mut ensg: HashMap<String, usize> = HashMap::default();
+            .collect();
+        let mut locus_raw: HashMap<Box<str>, usize> = HashMap::default();
         let mut locus: HashMap<Box<str>, usize> = HashMap::default();
-        for (i, (low, key)) in lowered.iter().zip(keys).enumerate() {
-            if let Some(key) = key {
+        let mut bare: Vec<(Box<str>, 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<String, usize> = HashMap::default();
+        let mut symbol: HashMap<String, usize> = HashMap::default();
+        let mut ensg: HashMap<String, usize> = HashMap::default();
+        for (i, low) in lowered.iter().enumerate() {
             if low.is_empty() {
                 continue;
             }
@@ -401,17 +457,22 @@ impl GeneIndex {
             exact,
             symbol,
             ensg,
+            locus_raw,
             locus,
         }
     }
 
     /// Row index for `gene`, or `None` if unmatched (tiers above).
     #[must_use]
     pub fn match_gene(&self, gene: &str) -> Option<usize> {
-        // 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) {

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);
+}