1use std::fmt;
4use std::io::Write;
5use std::ops::Range;
6
7use crate::error::{Error, Result};
8use crate::format::Format;
9use crate::qual::{self, PHRED33};
10use crate::seq::{self, Alphabet, BaseCounts};
11
12#[derive(Debug, Clone, Default, PartialEq, Eq, Hash)]
30pub struct Sequence {
31 pub id: String,
33 pub description: Option<String>,
35 pub seq: Vec<u8>,
37 pub quality: Option<Vec<u8>>,
39}
40
41impl Sequence {
42 pub fn fasta<I: Into<String>>(id: I, seq: impl Into<Vec<u8>>) -> Sequence {
44 Sequence {
45 id: id.into(),
46 description: None,
47 seq: seq.into(),
48 quality: None,
49 }
50 }
51
52 pub fn fastq<I: Into<String>>(
54 id: I,
55 seq: impl Into<Vec<u8>>,
56 quality: impl Into<Vec<u8>>,
57 ) -> Result<Sequence> {
58 let id = id.into();
59 let seq = seq.into();
60 let quality = quality.into();
61 if seq.len() != quality.len() {
62 return Err(Error::LengthMismatch {
63 id,
64 seq: seq.len(),
65 quality: quality.len(),
66 });
67 }
68 Ok(Sequence {
69 id,
70 description: None,
71 seq,
72 quality: Some(quality),
73 })
74 }
75
76 pub fn with_description<D: Into<String>>(mut self, description: D) -> Sequence {
78 let description = description.into();
79 self.description = if description.is_empty() {
80 None
81 } else {
82 Some(description)
83 };
84 self
85 }
86
87 pub fn len(&self) -> usize {
89 self.seq.len()
90 }
91
92 pub fn is_empty(&self) -> bool {
94 self.seq.is_empty()
95 }
96
97 pub fn has_quality(&self) -> bool {
99 self.quality.is_some()
100 }
101
102 pub fn format(&self) -> Format {
104 if self.has_quality() {
105 Format::Fastq
106 } else {
107 Format::Fasta
108 }
109 }
110
111 pub fn header(&self) -> String {
119 match &self.description {
120 Some(d) => format!("{} {}", self.id, d),
121 None => self.id.clone(),
122 }
123 }
124
125 pub fn seq_str(&self) -> Result<&str> {
127 std::str::from_utf8(&self.seq).map_err(|e| Error::InvalidByte {
128 id: self.id.clone(),
129 pos: e.valid_up_to(),
130 byte: self.seq[e.valid_up_to()],
131 })
132 }
133
134 pub fn clear(&mut self) {
138 self.id.clear();
139 self.description = None;
140 self.seq.clear();
141 if let Some(q) = self.quality.as_mut() {
142 q.clear();
143 }
144 }
145
146 pub fn base_counts(&self) -> BaseCounts {
148 BaseCounts::of(&self.seq)
149 }
150
151 pub fn gc_content(&self) -> Option<f64> {
153 seq::gc_content(&self.seq)
154 }
155
156 pub fn ambiguous_count(&self) -> u64 {
158 self.base_counts().n
159 }
160
161 pub fn reverse_complement(&self) -> Sequence {
171 Sequence {
172 id: self.id.clone(),
173 description: self.description.clone(),
174 seq: seq::reverse_complement(&self.seq),
175 quality: self
176 .quality
177 .as_ref()
178 .map(|q| q.iter().rev().copied().collect()),
179 }
180 }
181
182 pub fn reverse_complement_in_place(&mut self) {
184 seq::reverse_complement_in_place(&mut self.seq);
185 if let Some(q) = self.quality.as_mut() {
186 q.reverse();
187 }
188 }
189
190 pub fn make_uppercase(&mut self) {
192 self.seq.make_ascii_uppercase();
193 }
194
195 pub fn subseq(&self, range: Range<usize>) -> Result<Sequence> {
199 if range.start > range.end || range.end > self.seq.len() {
200 return Err(Error::OutOfBounds {
201 id: self.id.clone(),
202 start: range.start as u64,
203 end: range.end as u64,
204 length: self.seq.len() as u64,
205 });
206 }
207 Ok(Sequence {
208 id: self.id.clone(),
209 description: self.description.clone(),
210 seq: self.seq[range.clone()].to_vec(),
211 quality: self.quality.as_ref().map(|q| q[range].to_vec()),
212 })
213 }
214
215 pub fn trim_to(&mut self, range: Range<usize>) -> Result<()> {
217 if range.start > range.end || range.end > self.seq.len() {
218 return Err(Error::OutOfBounds {
219 id: self.id.clone(),
220 start: range.start as u64,
221 end: range.end as u64,
222 length: self.seq.len() as u64,
223 });
224 }
225 self.seq.truncate(range.end);
226 self.seq.drain(..range.start);
227 if let Some(q) = self.quality.as_mut() {
228 q.truncate(range.end);
229 q.drain(..range.start);
230 }
231 Ok(())
232 }
233
234 pub fn translate(&self, frame: usize, stop_at_stop: bool) -> Sequence {
236 Sequence {
237 id: self.id.clone(),
238 description: self.description.clone(),
239 seq: seq::translate(&self.seq, frame, stop_at_stop),
240 quality: None,
241 }
242 }
243
244 pub fn kmers(&self, k: usize) -> impl Iterator<Item = &[u8]> {
246 seq::kmers(&self.seq, k)
247 }
248
249 pub fn quality_scores(&self) -> Option<Vec<u8>> {
251 self.quality.as_ref().map(|q| qual::scores(q, PHRED33))
252 }
253
254 pub fn quality_scores_with(&self, offset: u8) -> Option<Vec<u8>> {
256 self.quality.as_ref().map(|q| qual::scores(q, offset))
257 }
258
259 pub fn mean_quality(&self) -> Option<f64> {
261 qual::mean_quality(self.quality.as_deref()?, PHRED33)
262 }
263
264 pub fn expected_errors(&self) -> Option<f64> {
266 Some(qual::expected_errors(self.quality.as_deref()?, PHRED33))
267 }
268
269 pub fn convert_quality_offset(&mut self, from: u8, to: u8) {
271 if let Some(q) = self.quality.as_mut() {
272 for c in q.iter_mut() {
273 *c = qual::encode(qual::score(*c, from), to);
274 }
275 }
276 }
277
278 pub fn into_fasta(mut self) -> Sequence {
280 self.quality = None;
281 self
282 }
283
284 pub fn validate(&self, alphabet: Alphabet) -> Result<()> {
292 if self.id.is_empty() {
293 return Err(Error::Parse {
294 line: 0,
295 kind: crate::error::ParseError::EmptyId,
296 });
297 }
298 if let Some(pos) = self.id.bytes().position(|b| b.is_ascii_whitespace()) {
299 return Err(Error::InvalidByte {
300 id: self.id.clone(),
301 pos,
302 byte: self.id.as_bytes()[pos],
303 });
304 }
305 if let Some(d) = &self.description {
306 if let Some(pos) = d.bytes().position(|b| b == b'\n' || b == b'\r') {
307 return Err(Error::InvalidByte {
308 id: self.id.clone(),
309 pos,
310 byte: d.as_bytes()[pos],
311 });
312 }
313 }
314 alphabet.validate_named(&self.seq, &self.id)?;
315 if let Some(q) = &self.quality {
316 if q.len() != self.seq.len() {
317 return Err(Error::LengthMismatch {
318 id: self.id.clone(),
319 seq: self.seq.len(),
320 quality: q.len(),
321 });
322 }
323 if let Some(pos) = q.iter().position(|&b| !(33..=126).contains(&b)) {
324 return Err(Error::InvalidByte {
325 id: self.id.clone(),
326 pos,
327 byte: q[pos],
328 });
329 }
330 }
331 Ok(())
332 }
333
334 pub fn write_fasta<W: Write>(&self, out: &mut W, line_width: Option<usize>) -> Result<()> {
342 check_writable_residues(&self.seq, &self.id)?;
343 out.write_all(b">")?;
344 self.write_header(out)?;
345 match line_width.filter(|w| *w > 0) {
346 None => {
347 out.write_all(&self.seq)?;
348 out.write_all(b"\n")?;
349 }
350 Some(width) => {
351 for chunk in self.seq.chunks(width) {
352 out.write_all(chunk)?;
353 out.write_all(b"\n")?;
354 }
355 if self.seq.is_empty() {
356 out.write_all(b"\n")?;
357 }
358 }
359 }
360 Ok(())
361 }
362
363 pub fn write_fastq<W: Write>(&self, out: &mut W) -> Result<()> {
366 let quality = self.quality.as_ref().ok_or_else(|| Error::MissingQuality {
367 id: self.id.clone(),
368 })?;
369 if quality.len() != self.seq.len() {
370 return Err(Error::LengthMismatch {
371 id: self.id.clone(),
372 seq: self.seq.len(),
373 quality: quality.len(),
374 });
375 }
376 check_writable_residues(&self.seq, &self.id)?;
377 check_writable_fastq_sequence(&self.seq, &self.id)?;
378 check_writable_quality(quality, &self.id)?;
379 out.write_all(b"@")?;
380 self.write_header(out)?;
381 out.write_all(&self.seq)?;
382 out.write_all(b"\n+\n")?;
383 out.write_all(quality)?;
384 out.write_all(b"\n")?;
385 Ok(())
386 }
387
388 fn write_header<W: Write>(&self, out: &mut W) -> Result<()> {
389 out.write_all(self.id.as_bytes())?;
390 if let Some(d) = &self.description {
391 out.write_all(b" ")?;
392 out.write_all(d.as_bytes())?;
393 }
394 out.write_all(b"\n")?;
395 Ok(())
396 }
397
398 pub fn to_string_in(&self, format: Format) -> Result<String> {
400 let mut buf = Vec::with_capacity(self.seq.len() * 2 + 64);
401 match format {
402 Format::Fasta => self.write_fasta(&mut buf, Some(60))?,
403 Format::Fastq => self.write_fastq(&mut buf)?,
404 }
405 Ok(String::from_utf8_lossy(&buf).into_owned())
406 }
407
408 pub(crate) fn set_header_reusing(&mut self, header: &[u8], spare: &mut String) {
425 let (id, description) = split_header(header);
426 push_lossy(&mut self.id, id);
427 match description {
428 None => self.description = None,
429 Some(description) => {
430 let mut buffer = std::mem::take(spare);
431 buffer.clear();
432 push_lossy(&mut buffer, description);
433 self.description = Some(buffer);
434 }
435 }
436 }
437}
438
439fn push_lossy(target: &mut String, bytes: &[u8]) {
444 match std::str::from_utf8(bytes) {
445 Ok(text) => target.push_str(text),
446 Err(_) => target.push_str(&String::from_utf8_lossy(bytes)),
447 }
448}
449
450pub(crate) fn split_header(header: &[u8]) -> (&[u8], Option<&[u8]>) {
456 let header = trim_ascii_end(header);
457 let id = header_id(header);
458 let description = trim_ascii_start(&header[id.len()..]);
459 (
460 id,
461 if description.is_empty() {
462 None
463 } else {
464 Some(description)
465 },
466 )
467}
468
469pub(crate) fn header_id(header: &[u8]) -> &[u8] {
477 let end = header
478 .iter()
479 .position(|b| b.is_ascii_whitespace())
480 .unwrap_or(header.len());
481 &header[..end]
482}
483
484pub(crate) fn check_writable_residues(seq: &[u8], id: &str) -> Result<()> {
499 match memchr::memchr3(b'\n', b'\r', b'>', seq) {
500 None => Ok(()),
501 Some(pos) => Err(Error::InvalidByte {
502 id: id.to_string(),
503 pos,
504 byte: seq[pos],
505 }),
506 }
507}
508
509pub(crate) fn check_writable_fastq_sequence(seq: &[u8], id: &str) -> Result<()> {
511 if seq.first() == Some(&b'+') {
512 return Err(Error::InvalidByte {
513 id: id.to_string(),
514 pos: 0,
515 byte: b'+',
516 });
517 }
518 Ok(())
519}
520
521pub(crate) fn check_writable_quality(quality: &[u8], id: &str) -> Result<()> {
523 match memchr::memchr2(b'\n', b'\r', quality) {
524 None => Ok(()),
525 Some(pos) => Err(Error::InvalidByte {
526 id: id.to_string(),
527 pos,
528 byte: quality[pos],
529 }),
530 }
531}
532
533fn trim_ascii_start(mut bytes: &[u8]) -> &[u8] {
535 while let [first, rest @ ..] = bytes {
536 if first.is_ascii_whitespace() {
537 bytes = rest;
538 } else {
539 break;
540 }
541 }
542 bytes
543}
544
545fn trim_ascii_end(mut bytes: &[u8]) -> &[u8] {
546 while let [rest @ .., last] = bytes {
547 if last.is_ascii_whitespace() {
548 bytes = rest;
549 } else {
550 break;
551 }
552 }
553 bytes
554}
555
556impl fmt::Display for Sequence {
557 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
559 let mut buf = Vec::new();
560 let rendered = match self.format() {
561 Format::Fasta => self.write_fasta(&mut buf, Some(60)),
562 Format::Fastq => self.write_fastq(&mut buf),
563 };
564 rendered.map_err(|_| fmt::Error)?;
565 f.write_str(&String::from_utf8_lossy(&buf))
566 }
567}
568
569#[cfg(test)]
570mod tests {
571 use super::*;
572
573 #[test]
574 fn fastq_requires_matching_lengths() {
575 assert!(Sequence::fastq("r", b"ACGT", b"III").is_err());
576 assert!(Sequence::fastq("r", b"ACGT", b"IIII").is_ok());
577 }
578
579 #[test]
580 fn header_splitting() {
581 let mut s = Sequence::default();
582 s.set_header_reusing(b"chr1 human chromosome 1", &mut String::new());
583 assert_eq!(s.id, "chr1");
584 assert_eq!(s.description.as_deref(), Some("human chromosome 1"));
585
586 let mut s = Sequence::default();
587 s.set_header_reusing(b"lonely", &mut String::new());
588 assert_eq!(s.id, "lonely");
589 assert_eq!(s.description, None);
590
591 let mut s = Sequence::default();
592 s.set_header_reusing(b"tabbed\tdesc with spaces ", &mut String::new());
593 assert_eq!(s.id, "tabbed");
594 assert_eq!(s.description.as_deref(), Some("desc with spaces"));
595 }
596
597 #[test]
598 fn header_splitting_handles_non_ascii() {
599 let mut s = Sequence::default();
602 s.set_header_reusing("chr\u{a0}1 description".as_bytes(), &mut String::new());
603 assert_eq!(s.id, "chr\u{a0}1");
604 assert_eq!(s.description.as_deref(), Some("description"));
605
606 let mut s = Sequence::default();
608 s.set_header_reusing(&[b'i', b'd', 0xff, b' ', b'd'], &mut String::new());
609 assert!(s.id.starts_with("id"));
610 assert_eq!(s.description.as_deref(), Some("d"));
611
612 let mut s = Sequence::default();
615 s.set_header_reusing(b" \t ", &mut String::new());
616 assert!(s.id.is_empty());
617 }
618
619 #[test]
620 fn fasta_wrapping() {
621 let s = Sequence::fasta("x", b"AAAAACCCCC".to_vec());
622 let mut out = Vec::new();
623 s.write_fasta(&mut out, Some(5)).unwrap();
624 assert_eq!(out, b">x\nAAAAA\nCCCCC\n");
625
626 let mut out = Vec::new();
627 s.write_fasta(&mut out, None).unwrap();
628 assert_eq!(out, b">x\nAAAAACCCCC\n");
629
630 let mut out = Vec::new();
632 Sequence::fasta("empty", Vec::new())
633 .write_fasta(&mut out, Some(60))
634 .unwrap();
635 assert_eq!(out, b">empty\n\n");
636 }
637
638 #[test]
639 fn fastq_output_and_missing_quality() {
640 let r = Sequence::fastq("r", b"ACGT", b"IIII")
641 .unwrap()
642 .with_description("d");
643 let mut out = Vec::new();
644 r.write_fastq(&mut out).unwrap();
645 assert_eq!(out, b"@r d\nACGT\n+\nIIII\n");
646
647 let mut out = Vec::new();
648 assert!(matches!(
649 Sequence::fasta("r", b"ACGT".to_vec()).write_fastq(&mut out),
650 Err(Error::MissingQuality { .. })
651 ));
652 }
653
654 #[test]
655 fn trimming_keeps_quality_aligned() {
656 let mut r = Sequence::fastq("r", b"AACCGGTT", b"01234567").unwrap();
657 r.trim_to(2..6).unwrap();
658 assert_eq!(r.seq, b"CCGG");
659 assert_eq!(r.quality.as_deref(), Some(&b"2345"[..]));
660 assert!(r.trim_to(0..99).is_err());
661 }
662
663 #[test]
664 fn subseq_bounds() {
665 let r = Sequence::fastq("r", b"AACCGG", b"012345").unwrap();
666 let sub = r.subseq(1..3).unwrap();
667 assert_eq!(sub.seq, b"AC");
668 assert_eq!(sub.quality.as_deref(), Some(&b"12"[..]));
669 #[allow(clippy::reversed_empty_ranges)] {
671 assert!(r.subseq(4..2).is_err());
672 }
673 assert!(r.subseq(0..7).is_err());
674 }
675
676 #[test]
677 fn validation_catches_bad_records() {
678 let mut r = Sequence::fasta("ok", b"ACGT".to_vec());
679 assert!(r.validate(Alphabet::Dna).is_ok());
680 r.id = "has space".into();
681 assert!(r.validate(Alphabet::Dna).is_err());
682 r.id = String::new();
683 assert!(r.validate(Alphabet::Dna).is_err());
684
685 let mut r = Sequence::fastq("r", b"ACGT", b"IIII").unwrap();
686 r.quality = Some(b"II".to_vec());
687 assert!(matches!(
688 r.validate(Alphabet::Dna),
689 Err(Error::LengthMismatch { .. })
690 ));
691 r.quality = Some(b"II\nI".to_vec());
692 assert!(matches!(
693 r.validate(Alphabet::Dna),
694 Err(Error::InvalidByte { .. })
695 ));
696 }
697
698 #[test]
699 fn refuses_to_write_unrepresentable_residues() {
700 let mut out = Vec::new();
704 for bad in [&b"AC\r"[..], b"AC\nGT", b"AC>GT", b">AC"] {
705 let record = Sequence::fasta("x", bad.to_vec());
706 assert!(
707 matches!(
708 record.write_fasta(&mut out, Some(60)),
709 Err(Error::InvalidByte { .. })
710 ),
711 "{:?} should be rejected",
712 String::from_utf8_lossy(bad)
713 );
714 }
715 for good in [&b"ACGTN"[..], b"acgt-n.", b"MEEPQSDPSV*"] {
717 let record = Sequence::fasta("x", good.to_vec());
718 assert!(record.write_fasta(&mut out, Some(60)).is_ok());
719 }
720 }
721
722 #[test]
723 fn refuses_to_write_unrepresentable_fastq() {
724 let mut out = Vec::new();
725
726 let record = Sequence::fastq("x", b"+CGT", b"IIII").unwrap();
728 assert!(matches!(
729 record.write_fastq(&mut out),
730 Err(Error::InvalidByte {
731 pos: 0,
732 byte: b'+',
733 ..
734 })
735 ));
736 assert!(Sequence::fastq("x", b"A+GT", b"IIII")
738 .unwrap()
739 .write_fastq(&mut out)
740 .is_ok());
741
742 let record = Sequence::fastq("x", b"ACGT", b"II\nI").unwrap();
744 assert!(matches!(
745 record.write_fastq(&mut out),
746 Err(Error::InvalidByte { byte: b'\n', .. })
747 ));
748 assert!(Sequence::fastq("x", b"ACGT", b"@@@@")
750 .unwrap()
751 .write_fastq(&mut out)
752 .is_ok());
753 }
754
755 #[test]
756 fn quality_offset_conversion() {
757 let mut r = Sequence::fastq("r", b"ACGT", b"hhhh").unwrap();
758 r.convert_quality_offset(crate::qual::PHRED64, PHRED33);
759 assert_eq!(r.quality.as_deref(), Some(&b"IIII"[..]));
760 }
761
762 #[test]
763 fn clear_keeps_capacity() {
764 let mut r = Sequence::fastq("r", b"ACGT", b"IIII").unwrap();
765 let cap = r.seq.capacity();
766 r.clear();
767 assert!(r.id.is_empty() && r.seq.is_empty());
768 assert_eq!(r.quality.as_deref(), Some(&b""[..]));
769 assert_eq!(r.seq.capacity(), cap);
770 }
771}