1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
//! Removing overlapping hits with [`dedup`].
//!
//! Hits are grouped by sequence, motif and strand, though motif and strand can
//! be ignored. Within a group, overlapping hits form clusters; in each cluster
//! hits are taken in order of priority, with ties going to the earlier hit, and
//! kept unless they overlap a hit already kept.
use crate::motif::Strand;
/// A hit as seen by [`dedup`]. Positions are half-open, and a higher
/// priority is better.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Interval {
/// Index of the sequence.
pub sequence: usize,
/// Index of the motif.
pub motif: usize,
/// The strand of the hit, with [`Strand::Both`] for hits on neither.
/// Hits only overlap others on the same strand, unless strands are
/// ignored.
pub strand: Strand,
/// First position, counting from 0.
pub start: i64,
/// Position just past the hit, so that it covers `start..end`.
pub end: i64,
/// How good the hit is; higher is better.
pub priority: f64,
}
/// Which hits to keep, in input order.
pub fn dedup(hits: &[Interval], ignore_strand: bool, ignore_motif: bool) -> Vec<bool> {
let n = hits.len();
let mut keep = vec![false; n];
let key = |i: usize| {
let h = &hits[i];
(
h.sequence,
if ignore_motif { 0 } else { h.motif },
if ignore_strand { 0 } else { h.strand as u8 },
)
};
let mut ord: Vec<usize> = (0..n).collect();
ord.sort_by_key(|&i| (key(i), hits[i].start, i));
let mut flush = |cluster: &[usize]| {
if let [only] = cluster {
keep[*only] = true;
return;
}
let mut by_priority = cluster.to_vec();
by_priority.sort_by(|&a, &b| {
hits[b]
.priority
.total_cmp(&hits[a].priority)
.then(a.cmp(&b))
});
let mut kept: Vec<(i64, i64)> = Vec::new();
for i in by_priority {
let (s, e) = (hits[i].start, hits[i].end);
if !kept.iter().any(|&(ks, ke)| s < ke && ks < e) {
keep[i] = true;
kept.push((s, e));
}
}
};
let mut g = 0;
while g < n {
let mut g_end = g + 1;
while g_end < n && key(ord[g_end]) == key(ord[g]) {
g_end += 1;
}
let mut lo = g;
let mut max_end = hits[ord[g]].end;
for k in g + 1..g_end {
let h = &hits[ord[k]];
if h.start >= max_end {
flush(&ord[lo..k]);
lo = k;
max_end = h.end;
} else if h.end > max_end {
max_end = h.end;
}
}
flush(&ord[lo..g_end]);
g = g_end;
}
keep
}