1use std::sync::Arc;
31
32use crate::sparse_io_vector::RowNameCanonicalizer;
33use genomic_data::coordinates::{self, chr_stripped, PeakCoord};
34use rustc_hash::FxHashMap as HashMap;
35
36#[derive(Clone, Debug, Default, PartialEq, Eq)]
41pub enum FeatureNameKind {
42 #[default]
44 Exact,
45 Gene { delim: char },
49 Locus { merge_overlapping: bool },
56 Mixed,
62}
63
64impl FeatureNameKind {
65 pub fn canonicalize(&self, name: &str) -> Box<str> {
72 match self {
73 FeatureNameKind::Exact => name.into(),
74 FeatureNameKind::Gene { delim } => gene_canonicalize(name, *delim),
75 FeatureNameKind::Locus { .. } => locus_key(name).unwrap_or_else(|| name.into()),
76 FeatureNameKind::Mixed => mixed_canonicalize(name),
77 }
78 }
79
80 pub fn is_exact(&self) -> bool {
82 matches!(self, FeatureNameKind::Exact)
83 }
84
85 pub fn needs_global_pass(&self) -> bool {
89 matches!(
90 self,
91 FeatureNameKind::Locus {
92 merge_overlapping: true
93 } | FeatureNameKind::Mixed
94 )
95 }
96
97 pub fn auto_detect(names: &[Box<str>]) -> Self {
104 let n = names.len();
105 if n == 0 {
106 return Self::Exact;
107 }
108 let mut n_locus = 0usize;
109 let mut n_gene_like = 0usize;
110 for name in names {
111 if coordinates::is_locus(name) {
112 n_locus += 1;
113 } else if is_gene_like(name, '_') {
114 n_gene_like += 1;
115 }
116 }
117 let pct_locus = n_locus as f32 / n as f32;
118 let pct_gene = n_gene_like as f32 / n as f32;
119 if pct_locus < 0.50 {
120 let n_spelled = names
121 .iter()
122 .filter(|name| {
123 !coordinates::is_locus(name) && coordinates::import_interval(name).is_some()
124 })
125 .count();
126 if n_spelled * 2 >= n {
127 log::warn!(
128 "{n_spelled} of {n} row names read as intervals only in a non-colon \
129 spelling (e.g. `chr1-100-200`); they are not loci here. Re-import them so \
130 peaks are named `chr:start-end`."
131 );
132 }
133 }
134 if pct_locus >= 0.10 && pct_gene >= 0.10 {
135 Self::Mixed
136 } else if pct_locus >= 0.50 {
137 Self::Locus {
138 merge_overlapping: true,
139 }
140 } else if pct_gene >= 0.50 {
141 Self::Gene { delim: '_' }
142 } else {
143 Self::Exact
144 }
145 }
146
147 #[must_use]
162 pub fn reconcile(kinds: &[FeatureNameKind]) -> FeatureNameKind {
163 if kinds.iter().any(|k| matches!(k, FeatureNameKind::Mixed)) {
164 return FeatureNameKind::Mixed;
165 }
166 let gene = kinds
167 .iter()
168 .find(|k| matches!(k, FeatureNameKind::Gene { .. }));
169 let locus = kinds
170 .iter()
171 .find(|k| matches!(k, FeatureNameKind::Locus { .. }));
172 match (gene, locus) {
173 (Some(_), Some(_)) => FeatureNameKind::Mixed,
174 _ => gene.or(locus).cloned().unwrap_or(FeatureNameKind::Exact),
175 }
176 }
177
178 pub fn into_canonicalizer(self) -> Option<RowNameCanonicalizer> {
184 if self.is_exact() {
185 return None;
186 }
187 Some(Arc::new(move |name: &str| self.canonicalize(name)))
188 }
189}
190
191pub fn parse_locus(name: &str) -> Option<(Box<str>, u64, u64)> {
198 let (chr, start, end) = coordinates::split_interval(name)?;
199 Some((chr_stripped(chr).into(), start as u64, end as u64))
200}
201
202pub use genomic_data::coordinates::locus_key;
206
207fn mixed_canonicalize(name: &str) -> Box<str> {
210 locus_key(name)
211 .or_else(|| gene_symbol(name, '_').map(Into::into))
212 .unwrap_or_else(|| name.into())
213}
214
215pub fn build_locus_overlap_canonical_map(names: &[Box<str>]) -> HashMap<Box<str>, Box<str>> {
226 let n = names.len();
227 let parsed: Vec<Option<PeakCoord>> = names
228 .iter()
229 .map(|n| coordinates::parse_interval(n))
230 .collect();
231
232 let mut by_chr: HashMap<&str, Vec<usize>> = HashMap::default();
234 for (i, p) in parsed.iter().enumerate() {
235 if let Some(p) = p {
236 by_chr.entry(chr_stripped(&p.chr)).or_default().push(i);
237 }
238 }
239
240 let mut parent: Vec<usize> = (0..n).collect();
242 fn find(p: &mut [usize], mut x: usize) -> usize {
243 while p[x] != x {
244 let g = p[p[x]];
245 p[x] = g;
246 x = g;
247 }
248 x
249 }
250
251 let mut cluster_extent: HashMap<usize, (i64, i64)> = HashMap::default();
253 for (_, mut idxs) in by_chr {
254 idxs.sort_by_key(|&i| parsed[i].as_ref().map_or(0, |p| p.start));
255 let mut current_root: Option<usize> = None;
256 let mut current_min_start: i64 = 0;
257 let mut current_max_end: i64 = 0;
258 for i in idxs {
259 let PeakCoord {
260 start: s, end: e, ..
261 } = parsed[i].as_ref().unwrap();
262 match current_root {
263 Some(root) if *s < current_max_end => {
264 let ra = find(&mut parent, root);
265 let rb = find(&mut parent, i);
266 if ra != rb {
267 parent[rb] = ra;
268 }
269 current_max_end = current_max_end.max(*e);
270 cluster_extent
271 .insert(find(&mut parent, i), (current_min_start, current_max_end));
272 }
273 _ => {
274 current_root = Some(i);
275 current_min_start = *s;
276 current_max_end = *e;
277 cluster_extent.insert(i, (*s, *e));
278 }
279 }
280 }
281 }
282
283 let mut out: HashMap<Box<str>, Box<str>> = HashMap::default();
285 for (i, p) in parsed.iter().enumerate() {
286 if let Some(p) = p {
287 let root = find(&mut parent, i);
288 let (start, end) = cluster_extent.get(&root).copied().unwrap_or((0, 0));
289 let cluster = PeakCoord {
292 chr: p.chr.clone(),
293 start,
294 end,
295 };
296 out.insert(names[i].clone(), cluster.locus_key());
297 }
298 }
299 out
300}
301
302pub fn build_locus_overlap_canonicalizer(names: &[Box<str>]) -> RowNameCanonicalizer {
308 let map = Arc::new(build_locus_overlap_canonical_map(names));
309 Arc::new(move |name: &str| {
310 map.get(name)
311 .cloned()
312 .or_else(|| locus_key(name))
313 .unwrap_or_else(|| name.into())
314 })
315}
316
317pub fn build_mixed_kind_canonicalizer(names: &[Box<str>]) -> RowNameCanonicalizer {
328 let map = Arc::new(build_locus_overlap_canonical_map(names));
329 Arc::new(move |name: &str| {
330 map.get(name)
331 .cloned()
332 .unwrap_or_else(|| mixed_canonicalize(name))
333 })
334}
335
336fn gene_canonicalize(name: &str, delim: char) -> Box<str> {
345 gene_symbol(name, delim).unwrap_or(name).into()
346}
347
348fn is_gene_like(name: &str, delim: char) -> bool {
350 gene_symbol(name, delim).is_some()
351}
352
353fn gene_symbol(name: &str, delim: char) -> Option<&str> {
360 if !name.contains(delim) {
361 return None;
362 }
363 let stripped = strip_feature_type_suffix(name, delim);
364 if coordinates::is_locus(stripped) {
365 return None;
366 }
367 let symbol = stripped.rsplit(delim).next().unwrap_or(stripped);
368 (!symbol.bytes().all(|b| b.is_ascii_digit())).then_some(symbol)
369}
370
371fn strip_feature_type_suffix(name: &str, delim: char) -> &str {
377 const TAGS: &[&str] = &[
382 "Gene_Expression",
383 "Gene",
384 "Antibody_Capture",
385 "CRISPR_Guide_Capture",
386 "Multiplexing_Capture",
387 "Custom",
388 "Peaks",
389 ];
390 for tag in TAGS {
391 if let Some(rest) = name.strip_suffix(tag).and_then(|r| r.strip_suffix(delim)) {
394 return rest;
395 }
396 }
397 name
398}
399
400#[derive(clap::ValueEnum, Clone, Debug, Default, serde::Serialize, serde::Deserialize)]
407#[serde(rename_all = "kebab-case")]
408pub enum FeatureNameKindArg {
409 #[default]
410 Auto,
411 Exact,
412 Gene,
413 Locus,
414 LocusOverlap,
415 Mixed,
416}
417
418impl FeatureNameKindArg {
419 pub fn resolve_or_gene(&self) -> FeatureNameKind {
423 Option::<FeatureNameKind>::from(self.clone())
424 .unwrap_or(FeatureNameKind::Gene { delim: '_' })
425 }
426}
427
428impl From<FeatureNameKindArg> for Option<FeatureNameKind> {
429 fn from(arg: FeatureNameKindArg) -> Self {
430 match arg {
431 FeatureNameKindArg::Auto => None,
432 FeatureNameKindArg::Exact => Some(FeatureNameKind::Exact),
433 FeatureNameKindArg::Gene => Some(FeatureNameKind::Gene { delim: '_' }),
434 FeatureNameKindArg::Locus => Some(FeatureNameKind::Locus {
435 merge_overlapping: false,
436 }),
437 FeatureNameKindArg::LocusOverlap => Some(FeatureNameKind::Locus {
438 merge_overlapping: true,
439 }),
440 FeatureNameKindArg::Mixed => Some(FeatureNameKind::Mixed),
441 }
442 }
443}
444
445#[cfg(test)]
446#[path = "feature_names_tests.rs"]
447mod feature_names_tests;
448
449#[cfg(test)]
450mod tests {
451 use super::*;
452
453 #[test]
454 fn exact_passthrough() {
455 let k = FeatureNameKind::Exact;
456 assert_eq!(
457 k.canonicalize("ENSG00000000003_TSPAN6").as_ref(),
458 "ENSG00000000003_TSPAN6"
459 );
460 assert!(k.is_exact());
461 assert!(k.into_canonicalizer().is_none());
462 }
463
464 #[test]
465 fn gene_takes_last_underscore_component() {
466 let k = FeatureNameKind::Gene { delim: '_' };
467 assert_eq!(k.canonicalize("ENSG00000000003_TSPAN6").as_ref(), "TSPAN6");
468 assert_eq!(k.canonicalize("TSPAN6").as_ref(), "TSPAN6");
470 assert_eq!(k.canonicalize("A_B_C").as_ref(), "C");
473 assert!(!k.is_exact());
474 assert!(k.into_canonicalizer().is_some());
475 }
476
477 #[test]
478 fn gene_strips_cell_ranger_feature_type_suffix() {
479 let k = FeatureNameKind::Gene { delim: '_' };
480 assert_eq!(
483 k.canonicalize("ENSG00000187634_SAMD11_Gene").as_ref(),
484 "SAMD11"
485 );
486 assert_eq!(
488 k.canonicalize("ENSG00000187634_SAMD11_Gene_Expression")
489 .as_ref(),
490 "SAMD11"
491 );
492 assert_eq!(k.canonicalize("FakeGene").as_ref(), "FakeGene");
495 }
496
497 #[test]
498 fn locus_strips_chr_and_leaves_other_spellings_alone() {
499 let k = FeatureNameKind::Locus {
500 merge_overlapping: false,
501 };
502 assert_eq!(k.canonicalize("chr1:1000-2000").as_ref(), "1:1000-2000");
503 assert_eq!(k.canonicalize("ChrX:5000-6000").as_ref(), "X:5000-6000");
504 assert_eq!(k.canonicalize("1_1000_2000").as_ref(), "1_1000_2000");
506 }
507
508 #[test]
511 fn parse_locus_accepts_common_formats() {
512 assert_eq!(
514 parse_locus("chr1:1000-2000"),
515 Some(("1".into(), 1000, 2000))
516 );
517 assert_eq!(parse_locus("1:1000-2000"), Some(("1".into(), 1000, 2000)));
518 assert_eq!(
519 parse_locus("CHR1:1000-2000"),
520 Some(("1".into(), 1000, 2000))
521 );
522 assert_eq!(
523 parse_locus("chrX:5000-6000"),
524 Some(("X".into(), 5000, 6000))
525 );
526 assert_eq!(parse_locus("chrMT:1-100"), Some(("MT".into(), 1, 100)));
527 }
528
529 #[test]
530 fn parse_locus_keeps_contig_names_with_separators() {
531 assert_eq!(
532 parse_locus("chrUn_CTG1v1:0-100"),
533 Some(("Un_CTG1v1".into(), 0, 100))
534 );
535 assert_eq!(
537 parse_locus("Un_CTG1v1:0-100"),
538 Some(("Un_CTG1v1".into(), 0, 100))
539 );
540 }
541
542 #[test]
543 fn contig_peaks_stay_loci_on_a_mixed_axis() {
544 let names: Vec<Box<str>> = vec![
545 "chr1_CTG1v1_random:5-10".into(),
546 "chr4_CTG2v2_random:5-10".into(),
547 "ENSG000_GENE1".into(),
548 ];
549 let canon = build_mixed_kind_canonicalizer(&names);
550 assert_eq!(canon(&names[0]).as_ref(), "1_CTG1v1_random:5-10");
551 assert_eq!(canon(&names[1]).as_ref(), "4_CTG2v2_random:5-10");
552 assert_eq!(canon(&names[2]).as_ref(), "GENE1");
553 }
554
555 #[test]
556 fn every_locus_path_gives_one_key() {
557 let names: Vec<Box<str>> = vec!["chrChr1:0-100".into(), "chr1:0-100".into()];
558 let map = build_locus_overlap_canonical_map(&names);
559 let k = FeatureNameKind::Locus {
560 merge_overlapping: false,
561 };
562 for name in &names {
563 assert_eq!(map.get(name).unwrap(), &k.canonicalize(name));
564 }
565 }
566
567 #[test]
568 fn locus_canonical_keeps_case_on_every_path() {
569 let names: Vec<Box<str>> = vec!["chrX:0-100".into(), "chr1:0-100".into()];
570 let map = build_locus_overlap_canonical_map(&names);
571 assert_eq!(map.get(&names[0]).unwrap().as_ref(), "X:0-100");
572 assert_eq!(map.get(&names[1]).unwrap().as_ref(), "1:0-100");
573 let canon = build_locus_overlap_canonicalizer(&names);
575 assert_eq!(canon("chrX:200-300").as_ref(), "X:200-300");
576 let mixed = build_mixed_kind_canonicalizer(&names);
577 assert_eq!(mixed("chrX:0-100").as_ref(), "X:0-100");
578 assert_eq!(mixed("chrM:200-300").as_ref(), "M:200-300");
579 let k = FeatureNameKind::Locus {
580 merge_overlapping: false,
581 };
582 assert_eq!(k.canonicalize("chrM:0-100").as_ref(), "M:0-100");
583 }
584
585 #[test]
586 fn parse_locus_rejects_non_loci() {
587 assert!(parse_locus("TGFB1").is_none()); assert!(parse_locus("ENSG00000105329").is_none()); assert!(parse_locus("chr1:bad-2000").is_none()); assert!(parse_locus("chr1:1000").is_none()); assert!(parse_locus("chr1:2000-1000").is_none()); assert!(parse_locus("").is_none()); assert!(parse_locus("chr1").is_none()); assert!(parse_locus("chr:1-2").is_none()); assert!(parse_locus("ENSG000_GENE1").is_none()); assert!(parse_locus("GENE1-AS1").is_none()); assert!(parse_locus("chr1_1000_2000").is_none()); assert!(parse_locus("chr1-1000-2000").is_none()); }
600
601 #[test]
602 fn overlap_map_merges_two_overlapping_intervals() {
603 let names = vec![
605 "chr1:1-20".to_string().into_boxed_str(),
606 "chr1:15-30".to_string().into_boxed_str(),
607 ];
608 let map = build_locus_overlap_canonical_map(&names);
609 let c0 = map.get(&names[0]).unwrap();
610 let c1 = map.get(&names[1]).unwrap();
611 assert_eq!(c0, c1, "both inputs should map to the same canonical");
612 assert_eq!(c0.as_ref(), "1:1-30"); }
614
615 #[test]
616 fn overlap_map_keeps_non_overlapping_separate() {
617 let names = vec![
618 "chr1:1-20".to_string().into_boxed_str(),
619 "chr1:100-200".to_string().into_boxed_str(),
620 "chr2:1-20".to_string().into_boxed_str(),
621 ];
622 let map = build_locus_overlap_canonical_map(&names);
623 assert_eq!(map.get(&names[0]).unwrap().as_ref(), "1:1-20");
624 assert_eq!(map.get(&names[1]).unwrap().as_ref(), "1:100-200");
625 assert_eq!(map.get(&names[2]).unwrap().as_ref(), "2:1-20");
627 }
628
629 #[test]
630 fn overlap_map_handles_transitive_chain() {
631 let names = vec![
635 "chr1:1-20".to_string().into_boxed_str(),
636 "chr1:15-30".to_string().into_boxed_str(),
637 "chr1:25-40".to_string().into_boxed_str(),
638 ];
639 let map = build_locus_overlap_canonical_map(&names);
640 let c0 = map.get(&names[0]).unwrap();
641 let c1 = map.get(&names[1]).unwrap();
642 let c2 = map.get(&names[2]).unwrap();
643 assert_eq!(c0, c1);
644 assert_eq!(c1, c2);
645 assert_eq!(c0.as_ref(), "1:1-40"); }
647
648 #[test]
649 fn overlap_map_handles_full_containment() {
650 let names = vec![
652 "chr1:1-100".to_string().into_boxed_str(),
653 "chr1:30-50".to_string().into_boxed_str(),
654 ];
655 let map = build_locus_overlap_canonical_map(&names);
656 let c0 = map.get(&names[0]).unwrap();
657 let c1 = map.get(&names[1]).unwrap();
658 assert_eq!(c0, c1);
659 assert_eq!(c0.as_ref(), "1:1-100");
660 }
661
662 #[test]
663 fn overlap_map_treats_adjacent_as_separate() {
664 let names = vec![
667 "chr1:1-20".to_string().into_boxed_str(),
668 "chr1:20-30".to_string().into_boxed_str(),
669 ];
670 let map = build_locus_overlap_canonical_map(&names);
671 assert_ne!(map.get(&names[0]).unwrap(), map.get(&names[1]).unwrap());
672 }
673
674 #[test]
675 fn overlap_map_normalizes_chr_prefix_within_cluster() {
676 let names = vec![
679 "chr1:1-20".to_string().into_boxed_str(),
680 "1:15-30".to_string().into_boxed_str(),
681 ];
682 let map = build_locus_overlap_canonical_map(&names);
683 let c0 = map.get(&names[0]).unwrap();
684 let c1 = map.get(&names[1]).unwrap();
685 assert_eq!(c0, c1);
686 assert_eq!(c0.as_ref(), "1:1-30");
687 }
688
689 #[test]
690 fn overlap_map_leaves_non_colon_spellings_out() {
691 let names = vec![
693 "chr1:1-20".to_string().into_boxed_str(),
694 "chr1_15_30".to_string().into_boxed_str(),
695 ];
696 let map = build_locus_overlap_canonical_map(&names);
697 assert_eq!(map.get(&names[0]).unwrap().as_ref(), "1:1-20");
698 assert!(!map.contains_key(&names[1]));
699 }
700
701 #[test]
702 fn underscore_peaks_are_not_collapsed_by_the_gene_rule() {
703 let names: Vec<Box<str>> = (1..=20)
706 .map(|i| format!("chr{i}_100_200").into_boxed_str())
707 .collect();
708 assert_eq!(FeatureNameKind::auto_detect(&names), FeatureNameKind::Exact);
709 let mixed = build_mixed_kind_canonicalizer(&names);
710 assert_eq!(mixed("chr2_100_200").as_ref(), "chr2_100_200");
711 let gene = FeatureNameKind::Gene { delim: '_' };
712 assert_eq!(gene.canonicalize("chr2_100_200").as_ref(), "chr2_100_200");
713 assert_eq!(gene.canonicalize("ENSG000_GENE1").as_ref(), "GENE1");
714 assert_eq!(gene.canonicalize("GENE1_Gene").as_ref(), "GENE1");
715 assert_eq!(
716 gene.canonicalize("chr2_100_200_Peaks").as_ref(),
717 "chr2_100_200_Peaks"
718 );
719 assert_eq!(
722 gene.canonicalize("chrUn_CTG1v1:0-100_Peaks").as_ref(),
723 "chrUn_CTG1v1:0-100_Peaks"
724 );
725 assert_eq!(
726 gene.canonicalize("chrUn_CTG1v1:0-100").as_ref(),
727 "chrUn_CTG1v1:0-100"
728 );
729 }
730
731 #[test]
732 fn overlap_map_ignores_non_locus_names() {
733 let names = vec![
736 "TGFB1".to_string().into_boxed_str(),
737 "chr1:1-20".to_string().into_boxed_str(),
738 ];
739 let map = build_locus_overlap_canonical_map(&names);
740 assert!(!map.contains_key(&names[0]));
741 assert!(map.contains_key(&names[1]));
742 }
743
744 #[test]
745 fn overlap_map_skips_empty_intervals() {
746 let names = vec!["chr1:1000-1000".to_string().into_boxed_str()];
748 let map = build_locus_overlap_canonical_map(&names);
749 assert!(map.is_empty());
750 }
751
752 #[test]
753 fn overlap_canonicalizer_falls_back_for_unmatched() {
754 let names = vec!["chr1:1-20".to_string().into_boxed_str()];
755 let canon = build_locus_overlap_canonicalizer(&names);
756 assert_eq!(canon("chr1:1-20").as_ref(), "1:1-20");
758 assert_eq!(canon("chr2:500-600").as_ref(), "2:500-600");
760 assert_eq!(canon("GENE1").as_ref(), "GENE1");
762 }
763
764 #[test]
767 fn auto_detect_pure_locus_axis() {
768 let names: Vec<Box<str>> = (0..100)
769 .map(|i| format!("chr1:{}-{}", i * 100, i * 100 + 50).into_boxed_str())
770 .collect();
771 assert!(matches!(
772 FeatureNameKind::auto_detect(&names),
773 FeatureNameKind::Locus {
774 merge_overlapping: true
775 }
776 ));
777 }
778
779 #[test]
780 fn auto_detect_pure_gene_axis() {
781 let names: Vec<Box<str>> = (0..100)
782 .map(|i| format!("ENSG000_GENE{}", i).into_boxed_str())
783 .collect();
784 assert!(matches!(
785 FeatureNameKind::auto_detect(&names),
786 FeatureNameKind::Gene { delim: '_' }
787 ));
788 }
789
790 #[test]
791 fn auto_detect_mixed_axis() {
792 let mut names: Vec<Box<str>> = (0..80)
794 .map(|i| format!("chr1:{}-{}", i * 1000, i * 1000 + 500).into_boxed_str())
795 .collect();
796 names.extend((0..20).map(|i| format!("ENSG000_GENE{}", i).into_boxed_str()));
797 assert!(matches!(
798 FeatureNameKind::auto_detect(&names),
799 FeatureNameKind::Mixed
800 ));
801 }
802
803 #[test]
804 fn auto_detect_empty_or_exact() {
805 assert!(matches!(
806 FeatureNameKind::auto_detect(&[]),
807 FeatureNameKind::Exact
808 ));
809 let names = vec!["TGFB1".into(), "CD4".into(), "IL2".into(), "GAPDH".into()];
810 assert!(matches!(
811 FeatureNameKind::auto_detect(&names),
812 FeatureNameKind::Exact
813 ));
814 }
815
816 #[test]
817 fn mixed_dispatcher_canonicalizes_each_name_by_kind() {
818 let names: Vec<Box<str>> = vec![
819 "chr1:1-20".into(), "chr1:15-30".into(), "ENSG000_TGFB1".into(), "CD4".into(), ];
824 let canon = build_mixed_kind_canonicalizer(&names);
825 assert_eq!(canon("chr1:1-20").as_ref(), "1:1-30");
826 assert_eq!(canon("chr1:15-30").as_ref(), "1:1-30");
827 assert_eq!(canon("ENSG000_TGFB1").as_ref(), "TGFB1");
828 assert_eq!(canon("CD4").as_ref(), "CD4");
829 }
830}