symcurve/
fasta.rs

1//! Functions for working with FASTA files.
2
3use std::sync::Arc;
4
5use noodles_core::Position;
6use noodles_fasta::Record;
7use noodles_fasta::record::Sequence;
8
9/// One Record will be split into multiple RecordPieces.
10/// The original Record is kept as an Arc so that each of the
11/// RecordPieces can share the same ownership. Arc rather than Rc because the
12/// pieces are scored in parallel, and Rc is not Send.
13pub struct RecordPiece {
14    pub record: Arc<Record>,
15    pub start: Position,
16    pub end: Position,
17}
18
19impl RecordPiece {
20    fn new(record: Arc<Record>, start: Position, end: Position) -> Self {
21        Self { record, start, end }
22    }
23
24    /// Get the sequence of the RecordPiece by slicing into the original Record.
25    ///
26    /// This copies the slice out. Prefer [`RecordPiece::bases`] when the bytes are only
27    /// going to be read, which is the case on the scoring path.
28    pub fn sequence(&self) -> Sequence {
29        self.record.sequence().slice(self.start..=self.end).unwrap()
30    }
31
32    /// Borrow this piece's bases directly out of the shared Record.
33    ///
34    /// Unlike [`RecordPiece::sequence`] this allocates nothing, so scoring a piece does
35    /// not begin by copying it. `start` and `end` are 1-based and inclusive.
36    pub fn bases(&self) -> &[u8] {
37        let start = usize::from(self.start) - 1;
38        let end = usize::from(self.end);
39        &self.record.sequence().as_ref()[start..end]
40    }
41}
42
43/// Returns true for any base that cannot be scored and must therefore break the sequence.
44///
45/// A, C, G and T are scoreable in either case: RepeatMasker lowercases repetitive regions,
46/// but soft-masking is an annotation rather than missing data, so `acgt` is ordinary
47/// sequence. Everything else -- `N` and the IUPAC ambiguity codes (R, Y, S, W, K, M, B,
48/// D, H, V) -- represents a base that is genuinely unknown.
49fn is_gap(base: u8) -> bool {
50    !matches!(base.to_ascii_uppercase(), b'A' | b'C' | b'G' | b'T')
51}
52
53#[allow(dead_code)]
54/// Given a record, split the sequence at runs of unscoreable bases.
55///
56/// Returns a vector of pieces, each covering a stretch of sequence that contains only
57/// A, C, G and T (in either case). Runs of `N` or IUPAC ambiguity codes are dropped, and
58/// the sequence is split there. Each piece records its own start-end position in the
59/// original record, the positions being 1-based.
60///
61/// Splitting matters because curvature is computed from a sliding window over a running
62/// sum of coordinates. Merely deleting the unknown bases would let a window span the gap
63/// and derive a value from bases that are far apart in the real sequence.
64///
65/// Input:
66/// ```text
67/// >chr42
68/// ATGCATGC
69/// NNNNATGC
70/// A
71/// ```
72///
73/// Output:
74/// ```text
75/// >chr42 1-8
76/// ATGCATGC
77/// >chr42 13-17
78/// ATGCA
79/// ```
80pub fn split_seq_by_gaps(record: Record) -> Vec<RecordPiece> {
81    // Move the record into a single Arc up front. Every piece then clones this
82    // one handle, so they all point at the same allocation. Calling Arc::new
83    // per piece would instead allocate a fresh box holding a full copy of the
84    // sequence, which is what the shared ownership here is meant to avoid.
85    let record = Arc::new(record);
86    let mut records = Vec::new();
87    let n = record.sequence().len();
88    let seq = record.sequence().as_ref();
89    let mut pos = 0;
90    // classic two-pointer approach is tried-and-true
91    // but might not be the most idiomatic Rust
92    while pos < n {
93        while (pos < n) && is_gap(seq[pos]) {
94            pos += 1;
95        }
96        let left = pos;
97        while (pos < n) && !is_gap(seq[pos]) {
98            pos += 1;
99        }
100        let right = pos;
101        if left < right {
102            // Position is 1-based so add 1 to left
103            let start = Position::try_from(left + 1).unwrap();
104            let end = Position::try_from(right).unwrap();
105            let piece = RecordPiece::new(Arc::clone(&record), start, end);
106            records.push(piece);
107        }
108    }
109    records
110}
111
112#[cfg(test)]
113mod tests {
114    use super::*;
115    use crate::curve::matrix;
116    use approx::assert_relative_eq;
117
118    #[test]
119    fn test_read_fasta() {
120        // two sequences
121        let src = b">sq0\nACGT\n>sq1\nN\n";
122        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
123        let first_rec = reader.records().next().unwrap().unwrap();
124        let second_rec = reader.records().next().unwrap().unwrap();
125        assert_eq!(first_rec.name(), b"sq0");
126        assert_eq!(second_rec.name(), b"sq1");
127        let start = Position::try_from(2).unwrap();
128        let end = Position::try_from(3).unwrap();
129        assert_eq!(
130            first_rec.sequence().slice(start..=end).unwrap().as_ref(),
131            b"CG".to_vec()
132        );
133    }
134
135    #[test]
136    fn test_windows() {
137        let seq = b"ACGTACGTACGTACGTACGT";
138        let mut seq_bytes = seq.windows(3);
139        assert_eq!(seq_bytes.next().unwrap(), b"ACG");
140        assert_eq!(seq_bytes.next().unwrap(), b"CGT");
141        assert_eq!(seq_bytes.next().unwrap(), b"GTA");
142    }
143
144    #[test]
145    fn test_splitting() {
146        let src = b">chr42\nATGCATGCNNNNATGCA\n";
147        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
148        let split_records: Vec<_> = reader
149            .records()
150            .flat_map(|rec| split_seq_by_gaps(rec.unwrap()))
151            .collect();
152        assert_eq!(split_records.len(), 2);
153        assert_eq!(split_records[0].sequence().as_ref(), b"ATGCATGC".to_vec());
154        assert_eq!(split_records[1].sequence().as_ref(), b"ATGCA".to_vec());
155        assert_eq!(split_records[1].sequence().as_ref(), b"ATGCA".to_vec());
156        assert_eq!(usize::from(split_records[1].start), 13);
157        assert_eq!(usize::from(split_records[1].end), 17);
158    }
159
160    #[test]
161    fn test_pieces_share_one_record() {
162        // Four pieces separated by runs of Ns. All of them must point at the
163        // same underlying Record rather than each holding their own copy.
164        let src = b">chr42\nACGTNNACGTNNACGTNNACGT\n";
165        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
166        let record = reader.records().next().unwrap().unwrap();
167        let pieces = split_seq_by_gaps(record);
168        assert_eq!(pieces.len(), 4);
169        // One allocation, one handle per piece.
170        assert_eq!(Arc::strong_count(&pieces[0].record), 4);
171        for piece in &pieces[1..] {
172            assert!(Arc::ptr_eq(&pieces[0].record, &piece.record));
173        }
174        // Sharing the record must not disturb the slicing.
175        assert_eq!(pieces[0].sequence().as_ref(), b"ACGT".to_vec());
176        assert_eq!(pieces[3].sequence().as_ref(), b"ACGT".to_vec());
177        assert_eq!(usize::from(pieces[3].start), 19);
178        assert_eq!(usize::from(pieces[3].end), 22);
179    }
180
181    #[test]
182    fn test_softmasked_bases_are_kept() {
183        // Lowercase acgt is soft-masked repeat sequence, not missing data, so it must
184        // flow through as ordinary sequence rather than splitting the record.
185        let src = b">chr42\nACGTacgtACGT\n";
186        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
187        let record = reader.records().next().unwrap().unwrap();
188        let pieces = split_seq_by_gaps(record);
189        assert_eq!(pieces.len(), 1);
190        assert_eq!(pieces[0].sequence().as_ref(), b"ACGTacgtACGT".to_vec());
191        assert_eq!(usize::from(pieces[0].start), 1);
192        assert_eq!(usize::from(pieces[0].end), 12);
193    }
194
195    #[test]
196    fn test_softmasked_bases_score_as_their_uppercase_form() {
197        // The lookup must agree across cases, or soft-masked regions would yield
198        // different curvature than the same sequence unmasked.
199        for (upper, lower) in [(b"ACG", b"acg"), (b"TTT", b"ttt"), (b"CCA", b"cca")] {
200            assert_relative_eq!(
201                matrix::matrix_lookup(upper, &matrix::ROLL_SIMPLE).unwrap(),
202                matrix::matrix_lookup(lower, &matrix::ROLL_SIMPLE).unwrap(),
203                epsilon = 1e-12
204            );
205        }
206    }
207
208    #[test]
209    fn test_ambiguity_codes_split_like_n() {
210        // R and Y are genuinely unknown bases. Before this they reached the matrix
211        // lookup and panicked; now they gap-split exactly as N does.
212        let src = b">chr42\nACGTRYACGT\n";
213        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
214        let record = reader.records().next().unwrap().unwrap();
215        let pieces = split_seq_by_gaps(record);
216        assert_eq!(pieces.len(), 2);
217        assert_eq!(pieces[0].sequence().as_ref(), b"ACGT".to_vec());
218        assert_eq!(pieces[1].sequence().as_ref(), b"ACGT".to_vec());
219        assert_eq!(usize::from(pieces[1].start), 7);
220        assert_eq!(usize::from(pieces[1].end), 10);
221    }
222
223    #[test]
224    fn test_every_iupac_code_is_a_gap() {
225        for &code in b"NRYSWKMBDHVnryswkmbdhv" {
226            assert!(is_gap(code), "{:?} should be a gap", code as char);
227        }
228        for &code in b"ACGTacgt" {
229            assert!(!is_gap(code), "{:?} should not be a gap", code as char);
230        }
231    }
232
233    #[test]
234    fn test_mixed_gaps_and_softmasking() {
235        // A realistic shape: soft-masked repeat, an assembly gap, an ambiguity code.
236        let src = b">chr42\nACGTacgtNNNNacgtRACGT\n";
237        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
238        let record = reader.records().next().unwrap().unwrap();
239        let pieces = split_seq_by_gaps(record);
240        assert_eq!(pieces.len(), 3);
241        assert_eq!(pieces[0].sequence().as_ref(), b"ACGTacgt".to_vec());
242        assert_eq!(pieces[1].sequence().as_ref(), b"acgt".to_vec());
243        assert_eq!(pieces[2].sequence().as_ref(), b"ACGT".to_vec());
244    }
245
246    #[test]
247    fn test_lookup_rejects_unknown_base_distinctly() {
248        // The two failure modes used to be conflated under "must be of length 3".
249        let bad_base = matrix::matrix_lookup(b"AAN", &matrix::ROLL_SIMPLE).unwrap_err();
250        assert!(
251            bad_base.to_string().contains("unrecognized nucleotide"),
252            "got: {bad_base}"
253        );
254        let bad_len = matrix::matrix_lookup(b"AA", &matrix::ROLL_SIMPLE).unwrap_err();
255        assert!(bad_len.to_string().contains("length 3"), "got: {bad_len}");
256    }
257
258    #[test]
259    fn test_splitting_empty() {
260        let src = b">chr42\n\n";
261        let mut reader = noodles_fasta::io::Reader::new(&src[..]);
262        let split_records: Vec<_> = reader
263            .records()
264            .flat_map(|rec| split_seq_by_gaps(rec.unwrap()))
265            .collect();
266        assert_eq!(split_records.len(), 0);
267    }
268}
symcurve/fasta.rs

symcurve/
fasta.rs
