salmon_model/
libdetect.rs

1//! Automatic library-type detection.
2//!
3//! Port of salmon's `LibraryTypeDetector`
4//! (`include/.../model/LibraryTypeDetector.hpp`). During the first reads of a
5//! run, the observed [`LibraryFormat`] of each confidently mapped fragment is
6//! tallied; once enough samples are seen, the most likely orientation and
7//! strandedness are inferred from the count ratios using salmon's 30%/70%
8//! thresholds.
9
10use salmon_core::{LibraryFormat, ReadOrientation, ReadStrandedness, ReadType};
11use std::sync::atomic::{AtomicBool, AtomicI64, AtomicU64, Ordering};
12
13/// Default number of samples to collect before guessing (matches salmon).
14pub const DEFAULT_SAMPLES_NEEDED: i64 = 50_000;
15
16/// Accumulates observed library formats and infers the most likely type.
17#[derive(Debug)]
18pub struct LibraryTypeDetector {
19    active: AtomicBool,
20    read_type: ReadType,
21    samples_needed: AtomicI64,
22    counts: Vec<AtomicU64>,
23}
24
25impl LibraryTypeDetector {
26    pub fn new(read_type: ReadType) -> Self {
27        let counts = (0..=LibraryFormat::MAX_FORMAT_ID)
28            .map(|_| AtomicU64::new(0))
29            .collect();
30        Self {
31            active: AtomicBool::new(true),
32            read_type,
33            samples_needed: AtomicI64::new(DEFAULT_SAMPLES_NEEDED),
34            counts,
35        }
36    }
37
38    pub fn is_active(&self) -> bool {
39        self.active.load(Ordering::Relaxed)
40    }
41
42    /// True once enough samples have been collected to guess.
43    pub fn can_guess(&self) -> bool {
44        self.samples_needed.load(Ordering::Relaxed) <= 0
45    }
46
47    /// Record one confidently mapped fragment's observed format. Only formats
48    /// matching the detector's read type are counted, and only until the sample
49    /// budget is exhausted. Thread-safe.
50    pub fn add_sample(&self, f: LibraryFormat) {
51        if f.read_type == self.read_type && self.samples_needed.load(Ordering::Relaxed) >= 0 {
52            self.counts[f.format_id() as usize].fetch_add(1, Ordering::Relaxed);
53            self.samples_needed.fetch_sub(1, Ordering::Relaxed);
54        }
55    }
56
57    /// Infer the most likely library format from the accumulated counts. Marks
58    /// the detector inactive (so this is effectively done once). Returns `None`
59    /// if the detector was already inactive.
60    pub fn most_likely_type(&self) -> Option<LibraryFormat> {
61        if !self.active.swap(false, Ordering::AcqRel) {
62            return None;
63        }
64
65        let count = |id: u8| self.counts[id as usize].load(Ordering::Relaxed);
66
67        let fmt = match self.read_type {
68            ReadType::SingleEnd => {
69                let mut nf = 0u64;
70                let mut nr = 0u64;
71                for id in 0..=LibraryFormat::MAX_FORMAT_ID {
72                    let f = LibraryFormat::from_format_id(id);
73                    let c = count(id);
74                    nf += if f.strandedness == ReadStrandedness::S {
75                        c
76                    } else {
77                        0
78                    };
79                    nr += if f.strandedness == ReadStrandedness::A {
80                        c
81                    } else {
82                        0
83                    };
84                }
85                let strandedness = if nf + nr == 0 {
86                    ReadStrandedness::U
87                } else {
88                    // Single-end uses the matching (S/A) encoding, like a
89                    // paired "same"-orientation library.
90                    strandedness_from_fw_ratio(nf as f64 / (nf + nr) as f64, true)
91                };
92                LibraryFormat::new(ReadType::SingleEnd, ReadOrientation::None, strandedness)
93            }
94            ReadType::PairedEnd => {
95                let (mut nsf, mut nsr) = (0u64, 0u64);
96                let (mut ninward, mut noutward, mut nsame) = (0u64, 0u64, 0u64);
97                for id in 0..=LibraryFormat::MAX_FORMAT_ID {
98                    let f = LibraryFormat::from_format_id(id);
99                    let c = count(id);
100                    nsf += matches!(f.strandedness, ReadStrandedness::S | ReadStrandedness::SA)
101                        .then_some(c)
102                        .unwrap_or(0);
103                    nsr += matches!(f.strandedness, ReadStrandedness::A | ReadStrandedness::AS)
104                        .then_some(c)
105                        .unwrap_or(0);
106                    match f.orientation {
107                        ReadOrientation::Toward => ninward += c,
108                        ReadOrientation::Away => noutward += c,
109                        ReadOrientation::Same => nsame += c,
110                        ReadOrientation::None => {}
111                    }
112                }
113
114                let num_orient = ninward + noutward + nsame;
115                if num_orient > 0 && (nsf + nsr) > 0 {
116                    let ratio_in = ninward as f64 / num_orient as f64;
117                    let ratio_out = noutward as f64 / num_orient as f64;
118                    let ratio_same = nsame as f64 / num_orient as f64;
119
120                    let (orientation, same) = if ratio_in >= ratio_out && ratio_in >= ratio_same {
121                        (ReadOrientation::Toward, false)
122                    } else if ratio_out >= ratio_in && ratio_out >= ratio_same {
123                        (ReadOrientation::Away, false)
124                    } else {
125                        (ReadOrientation::Same, true)
126                    };
127
128                    let ratio_fw = nsf as f64 / (nsf + nsr) as f64;
129                    let strandedness = strandedness_from_fw_ratio(ratio_fw, same);
130                    LibraryFormat::new(ReadType::PairedEnd, orientation, strandedness)
131                } else {
132                    LibraryFormat::new(
133                        ReadType::PairedEnd,
134                        ReadOrientation::Toward,
135                        ReadStrandedness::U,
136                    )
137                }
138            }
139        };
140
141        Some(fmt)
142    }
143}
144
145/// Map a forward-strand fraction to a strandedness using salmon's 30%/70%
146/// thresholds. `same` selects between the matching (S/A) and opposite (SA/AS)
147/// stranded encodings for paired-end "same"-orientation libraries; for
148/// single-end pass `false`.
149fn strandedness_from_fw_ratio(ratio_fw: f64, same: bool) -> ReadStrandedness {
150    if ratio_fw < 0.3 {
151        if same {
152            ReadStrandedness::A
153        } else {
154            ReadStrandedness::AS
155        }
156    } else if ratio_fw < 0.7 {
157        ReadStrandedness::U
158    } else if same {
159        ReadStrandedness::S
160    } else {
161        ReadStrandedness::SA
162    }
163}
164
165#[cfg(test)]
166mod tests {
167    use super::*;
168
169    #[test]
170    fn single_end_detects_sense() {
171        let d = LibraryTypeDetector::new(ReadType::SingleEnd);
172        let sf = LibraryFormat::parse("SF").unwrap();
173        let sr = LibraryFormat::parse("SR").unwrap();
174        for _ in 0..90 {
175            d.add_sample(sf);
176        }
177        for _ in 0..10 {
178            d.add_sample(sr);
179        }
180        let f = d.most_likely_type().unwrap();
181        assert_eq!(f.canonical(), "SF");
182        // detector becomes inactive after guessing
183        assert!(d.most_likely_type().is_none());
184    }
185
186    #[test]
187    fn single_end_detects_unstranded() {
188        let d = LibraryTypeDetector::new(ReadType::SingleEnd);
189        let sf = LibraryFormat::parse("SF").unwrap();
190        let sr = LibraryFormat::parse("SR").unwrap();
191        for _ in 0..50 {
192            d.add_sample(sf);
193            d.add_sample(sr);
194        }
195        assert_eq!(d.most_likely_type().unwrap().canonical(), "U");
196    }
197
198    #[test]
199    fn paired_end_detects_isr() {
200        let d = LibraryTypeDetector::new(ReadType::PairedEnd);
201        let isr = LibraryFormat::parse("ISR").unwrap();
202        for _ in 0..100 {
203            d.add_sample(isr);
204        }
205        // ISR: inward + antisense -> toward + AS
206        assert_eq!(d.most_likely_type().unwrap().canonical(), "ISR");
207    }
208
209    #[test]
210    fn paired_end_detects_iu() {
211        let d = LibraryTypeDetector::new(ReadType::PairedEnd);
212        let isf = LibraryFormat::parse("ISF").unwrap();
213        let isr = LibraryFormat::parse("ISR").unwrap();
214        for _ in 0..50 {
215            d.add_sample(isf);
216            d.add_sample(isr);
217        }
218        // balanced strandedness -> unstranded, inward -> IU
219        assert_eq!(d.most_likely_type().unwrap().canonical(), "IU");
220    }
221
222    #[test]
223    fn sample_budget_is_respected() {
224        let d = LibraryTypeDetector::new(ReadType::SingleEnd);
225        assert!(!d.can_guess());
226        let sf = LibraryFormat::parse("SF").unwrap();
227        // exhaust the budget
228        let mut n = DEFAULT_SAMPLES_NEEDED + 5;
229        while n > 0 {
230            d.add_sample(sf);
231            n -= 1;
232        }
233        assert!(d.can_guess());
234    }
235}
salmon_model/libdetect.rs

salmon_model/
libdetect.rs
