use std::fmt;
use std::io::Write;
use std::ops::Range;
use crate::error::{Error, Result};
use crate::format::Format;
use crate::qual::{self, PHRED33};
use crate::seq::{self, Alphabet, BaseCounts};
#[derive(Debug, Clone, Default, PartialEq, Eq, Hash)]
pub struct Sequence {
pub id: String,
pub description: Option<String>,
pub seq: Vec<u8>,
pub quality: Option<Vec<u8>>,
}
impl Sequence {
pub fn fasta<I: Into<String>>(id: I, seq: impl Into<Vec<u8>>) -> Sequence {
Sequence {
id: id.into(),
description: None,
seq: seq.into(),
quality: None,
}
}
pub fn fastq<I: Into<String>>(
id: I,
seq: impl Into<Vec<u8>>,
quality: impl Into<Vec<u8>>,
) -> Result<Sequence> {
let id = id.into();
let seq = seq.into();
let quality = quality.into();
if seq.len() != quality.len() {
return Err(Error::LengthMismatch {
id,
seq: seq.len(),
quality: quality.len(),
});
}
Ok(Sequence {
id,
description: None,
seq,
quality: Some(quality),
})
}
pub fn with_description<D: Into<String>>(mut self, description: D) -> Sequence {
let description = description.into();
self.description = if description.is_empty() {
None
} else {
Some(description)
};
self
}
pub fn len(&self) -> usize {
self.seq.len()
}
pub fn is_empty(&self) -> bool {
self.seq.is_empty()
}
pub fn has_quality(&self) -> bool {
self.quality.is_some()
}
pub fn format(&self) -> Format {
if self.has_quality() {
Format::Fastq
} else {
Format::Fasta
}
}
pub fn header(&self) -> String {
match &self.description {
Some(d) => format!("{} {}", self.id, d),
None => self.id.clone(),
}
}
pub fn seq_str(&self) -> Result<&str> {
std::str::from_utf8(&self.seq).map_err(|e| Error::InvalidByte {
id: self.id.clone(),
pos: e.valid_up_to(),
byte: self.seq[e.valid_up_to()],
})
}
pub fn clear(&mut self) {
self.id.clear();
self.description = None;
self.seq.clear();
if let Some(q) = self.quality.as_mut() {
q.clear();
}
}
pub fn base_counts(&self) -> BaseCounts {
BaseCounts::of(&self.seq)
}
pub fn gc_content(&self) -> Option<f64> {
seq::gc_content(&self.seq)
}
pub fn ambiguous_count(&self) -> u64 {
self.base_counts().n
}
pub fn reverse_complement(&self) -> Sequence {
Sequence {
id: self.id.clone(),
description: self.description.clone(),
seq: seq::reverse_complement(&self.seq),
quality: self
.quality
.as_ref()
.map(|q| q.iter().rev().copied().collect()),
}
}
pub fn reverse_complement_in_place(&mut self) {
seq::reverse_complement_in_place(&mut self.seq);
if let Some(q) = self.quality.as_mut() {
q.reverse();
}
}
pub fn make_uppercase(&mut self) {
self.seq.make_ascii_uppercase();
}
pub fn subseq(&self, range: Range<usize>) -> Result<Sequence> {
if range.start > range.end || range.end > self.seq.len() {
return Err(Error::OutOfBounds {
id: self.id.clone(),
start: range.start as u64,
end: range.end as u64,
length: self.seq.len() as u64,
});
}
Ok(Sequence {
id: self.id.clone(),
description: self.description.clone(),
seq: self.seq[range.clone()].to_vec(),
quality: self.quality.as_ref().map(|q| q[range].to_vec()),
})
}
pub fn trim_to(&mut self, range: Range<usize>) -> Result<()> {
if range.start > range.end || range.end > self.seq.len() {
return Err(Error::OutOfBounds {
id: self.id.clone(),
start: range.start as u64,
end: range.end as u64,
length: self.seq.len() as u64,
});
}
self.seq.truncate(range.end);
self.seq.drain(..range.start);
if let Some(q) = self.quality.as_mut() {
q.truncate(range.end);
q.drain(..range.start);
}
Ok(())
}
pub fn translate(&self, frame: usize, stop_at_stop: bool) -> Sequence {
Sequence {
id: self.id.clone(),
description: self.description.clone(),
seq: seq::translate(&self.seq, frame, stop_at_stop),
quality: None,
}
}
pub fn kmers(&self, k: usize) -> impl Iterator<Item = &[u8]> {
seq::kmers(&self.seq, k)
}
pub fn quality_scores(&self) -> Option<Vec<u8>> {
self.quality.as_ref().map(|q| qual::scores(q, PHRED33))
}
pub fn quality_scores_with(&self, offset: u8) -> Option<Vec<u8>> {
self.quality.as_ref().map(|q| qual::scores(q, offset))
}
pub fn mean_quality(&self) -> Option<f64> {
qual::mean_quality(self.quality.as_deref()?, PHRED33)
}
pub fn expected_errors(&self) -> Option<f64> {
Some(qual::expected_errors(self.quality.as_deref()?, PHRED33))
}
pub fn convert_quality_offset(&mut self, from: u8, to: u8) {
if let Some(q) = self.quality.as_mut() {
for c in q.iter_mut() {
*c = qual::encode(qual::score(*c, from), to);
}
}
}
pub fn into_fasta(mut self) -> Sequence {
self.quality = None;
self
}
pub fn validate(&self, alphabet: Alphabet) -> Result<()> {
if self.id.is_empty() {
return Err(Error::Parse {
line: 0,
kind: crate::error::ParseError::EmptyId,
});
}
if let Some(pos) = self.id.bytes().position(|b| b.is_ascii_whitespace()) {
return Err(Error::InvalidByte {
id: self.id.clone(),
pos,
byte: self.id.as_bytes()[pos],
});
}
if let Some(d) = &self.description {
if let Some(pos) = d.bytes().position(|b| b == b'\n' || b == b'\r') {
return Err(Error::InvalidByte {
id: self.id.clone(),
pos,
byte: d.as_bytes()[pos],
});
}
}
alphabet.validate_named(&self.seq, &self.id)?;
if let Some(q) = &self.quality {
if q.len() != self.seq.len() {
return Err(Error::LengthMismatch {
id: self.id.clone(),
seq: self.seq.len(),
quality: q.len(),
});
}
if let Some(pos) = q.iter().position(|&b| !(33..=126).contains(&b)) {
return Err(Error::InvalidByte {
id: self.id.clone(),
pos,
byte: q[pos],
});
}
}
Ok(())
}
pub fn write_fasta<W: Write>(&self, out: &mut W, line_width: Option<usize>) -> Result<()> {
check_writable_residues(&self.seq, &self.id)?;
out.write_all(b">")?;
self.write_header(out)?;
match line_width.filter(|w| *w > 0) {
None => {
out.write_all(&self.seq)?;
out.write_all(b"\n")?;
}
Some(width) => {
for chunk in self.seq.chunks(width) {
out.write_all(chunk)?;
out.write_all(b"\n")?;
}
if self.seq.is_empty() {
out.write_all(b"\n")?;
}
}
}
Ok(())
}
pub fn write_fastq<W: Write>(&self, out: &mut W) -> Result<()> {
let quality = self.quality.as_ref().ok_or_else(|| Error::MissingQuality {
id: self.id.clone(),
})?;
if quality.len() != self.seq.len() {
return Err(Error::LengthMismatch {
id: self.id.clone(),
seq: self.seq.len(),
quality: quality.len(),
});
}
check_writable_residues(&self.seq, &self.id)?;
check_writable_fastq_sequence(&self.seq, &self.id)?;
check_writable_quality(quality, &self.id)?;
out.write_all(b"@")?;
self.write_header(out)?;
out.write_all(&self.seq)?;
out.write_all(b"\n+\n")?;
out.write_all(quality)?;
out.write_all(b"\n")?;
Ok(())
}
fn write_header<W: Write>(&self, out: &mut W) -> Result<()> {
out.write_all(self.id.as_bytes())?;
if let Some(d) = &self.description {
out.write_all(b" ")?;
out.write_all(d.as_bytes())?;
}
out.write_all(b"\n")?;
Ok(())
}
pub fn to_string_in(&self, format: Format) -> Result<String> {
let mut buf = Vec::with_capacity(self.seq.len() * 2 + 64);
match format {
Format::Fasta => self.write_fasta(&mut buf, Some(60))?,
Format::Fastq => self.write_fastq(&mut buf)?,
}
Ok(String::from_utf8_lossy(&buf).into_owned())
}
pub(crate) fn set_header_reusing(&mut self, header: &[u8], spare: &mut String) {
let (id, description) = split_header(header);
push_lossy(&mut self.id, id);
match description {
None => self.description = None,
Some(description) => {
let mut buffer = std::mem::take(spare);
buffer.clear();
push_lossy(&mut buffer, description);
self.description = Some(buffer);
}
}
}
}
fn push_lossy(target: &mut String, bytes: &[u8]) {
match std::str::from_utf8(bytes) {
Ok(text) => target.push_str(text),
Err(_) => target.push_str(&String::from_utf8_lossy(bytes)),
}
}
pub(crate) fn split_header(header: &[u8]) -> (&[u8], Option<&[u8]>) {
let header = trim_ascii_end(header);
let id = header_id(header);
let description = trim_ascii_start(&header[id.len()..]);
(
id,
if description.is_empty() {
None
} else {
Some(description)
},
)
}
pub(crate) fn header_id(header: &[u8]) -> &[u8] {
let end = header
.iter()
.position(|b| b.is_ascii_whitespace())
.unwrap_or(header.len());
&header[..end]
}
pub(crate) fn check_writable_residues(seq: &[u8], id: &str) -> Result<()> {
match memchr::memchr3(b'\n', b'\r', b'>', seq) {
None => Ok(()),
Some(pos) => Err(Error::InvalidByte {
id: id.to_string(),
pos,
byte: seq[pos],
}),
}
}
pub(crate) fn check_writable_fastq_sequence(seq: &[u8], id: &str) -> Result<()> {
if seq.first() == Some(&b'+') {
return Err(Error::InvalidByte {
id: id.to_string(),
pos: 0,
byte: b'+',
});
}
Ok(())
}
pub(crate) fn check_writable_quality(quality: &[u8], id: &str) -> Result<()> {
match memchr::memchr2(b'\n', b'\r', quality) {
None => Ok(()),
Some(pos) => Err(Error::InvalidByte {
id: id.to_string(),
pos,
byte: quality[pos],
}),
}
}
fn trim_ascii_start(mut bytes: &[u8]) -> &[u8] {
while let [first, rest @ ..] = bytes {
if first.is_ascii_whitespace() {
bytes = rest;
} else {
break;
}
}
bytes
}
fn trim_ascii_end(mut bytes: &[u8]) -> &[u8] {
while let [rest @ .., last] = bytes {
if last.is_ascii_whitespace() {
bytes = rest;
} else {
break;
}
}
bytes
}
impl fmt::Display for Sequence {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let mut buf = Vec::new();
let rendered = match self.format() {
Format::Fasta => self.write_fasta(&mut buf, Some(60)),
Format::Fastq => self.write_fastq(&mut buf),
};
rendered.map_err(|_| fmt::Error)?;
f.write_str(&String::from_utf8_lossy(&buf))
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn fastq_requires_matching_lengths() {
assert!(Sequence::fastq("r", b"ACGT", b"III").is_err());
assert!(Sequence::fastq("r", b"ACGT", b"IIII").is_ok());
}
#[test]
fn header_splitting() {
let mut s = Sequence::default();
s.set_header_reusing(b"chr1 human chromosome 1", &mut String::new());
assert_eq!(s.id, "chr1");
assert_eq!(s.description.as_deref(), Some("human chromosome 1"));
let mut s = Sequence::default();
s.set_header_reusing(b"lonely", &mut String::new());
assert_eq!(s.id, "lonely");
assert_eq!(s.description, None);
let mut s = Sequence::default();
s.set_header_reusing(b"tabbed\tdesc with spaces ", &mut String::new());
assert_eq!(s.id, "tabbed");
assert_eq!(s.description.as_deref(), Some("desc with spaces"));
}
#[test]
fn header_splitting_handles_non_ascii() {
let mut s = Sequence::default();
s.set_header_reusing("chr\u{a0}1 description".as_bytes(), &mut String::new());
assert_eq!(s.id, "chr\u{a0}1");
assert_eq!(s.description.as_deref(), Some("description"));
let mut s = Sequence::default();
s.set_header_reusing(&[b'i', b'd', 0xff, b' ', b'd'], &mut String::new());
assert!(s.id.starts_with("id"));
assert_eq!(s.description.as_deref(), Some("d"));
let mut s = Sequence::default();
s.set_header_reusing(b" \t ", &mut String::new());
assert!(s.id.is_empty());
}
#[test]
fn fasta_wrapping() {
let s = Sequence::fasta("x", b"AAAAACCCCC".to_vec());
let mut out = Vec::new();
s.write_fasta(&mut out, Some(5)).unwrap();
assert_eq!(out, b">x\nAAAAA\nCCCCC\n");
let mut out = Vec::new();
s.write_fasta(&mut out, None).unwrap();
assert_eq!(out, b">x\nAAAAACCCCC\n");
let mut out = Vec::new();
Sequence::fasta("empty", Vec::new())
.write_fasta(&mut out, Some(60))
.unwrap();
assert_eq!(out, b">empty\n\n");
}
#[test]
fn fastq_output_and_missing_quality() {
let r = Sequence::fastq("r", b"ACGT", b"IIII")
.unwrap()
.with_description("d");
let mut out = Vec::new();
r.write_fastq(&mut out).unwrap();
assert_eq!(out, b"@r d\nACGT\n+\nIIII\n");
let mut out = Vec::new();
assert!(matches!(
Sequence::fasta("r", b"ACGT".to_vec()).write_fastq(&mut out),
Err(Error::MissingQuality { .. })
));
}
#[test]
fn trimming_keeps_quality_aligned() {
let mut r = Sequence::fastq("r", b"AACCGGTT", b"01234567").unwrap();
r.trim_to(2..6).unwrap();
assert_eq!(r.seq, b"CCGG");
assert_eq!(r.quality.as_deref(), Some(&b"2345"[..]));
assert!(r.trim_to(0..99).is_err());
}
#[test]
fn subseq_bounds() {
let r = Sequence::fastq("r", b"AACCGG", b"012345").unwrap();
let sub = r.subseq(1..3).unwrap();
assert_eq!(sub.seq, b"AC");
assert_eq!(sub.quality.as_deref(), Some(&b"12"[..]));
#[allow(clippy::reversed_empty_ranges)] {
assert!(r.subseq(4..2).is_err());
}
assert!(r.subseq(0..7).is_err());
}
#[test]
fn validation_catches_bad_records() {
let mut r = Sequence::fasta("ok", b"ACGT".to_vec());
assert!(r.validate(Alphabet::Dna).is_ok());
r.id = "has space".into();
assert!(r.validate(Alphabet::Dna).is_err());
r.id = String::new();
assert!(r.validate(Alphabet::Dna).is_err());
let mut r = Sequence::fastq("r", b"ACGT", b"IIII").unwrap();
r.quality = Some(b"II".to_vec());
assert!(matches!(
r.validate(Alphabet::Dna),
Err(Error::LengthMismatch { .. })
));
r.quality = Some(b"II\nI".to_vec());
assert!(matches!(
r.validate(Alphabet::Dna),
Err(Error::InvalidByte { .. })
));
}
#[test]
fn refuses_to_write_unrepresentable_residues() {
let mut out = Vec::new();
for bad in [&b"AC\r"[..], b"AC\nGT", b"AC>GT", b">AC"] {
let record = Sequence::fasta("x", bad.to_vec());
assert!(
matches!(
record.write_fasta(&mut out, Some(60)),
Err(Error::InvalidByte { .. })
),
"{:?} should be rejected",
String::from_utf8_lossy(bad)
);
}
for good in [&b"ACGTN"[..], b"acgt-n.", b"MEEPQSDPSV*"] {
let record = Sequence::fasta("x", good.to_vec());
assert!(record.write_fasta(&mut out, Some(60)).is_ok());
}
}
#[test]
fn refuses_to_write_unrepresentable_fastq() {
let mut out = Vec::new();
let record = Sequence::fastq("x", b"+CGT", b"IIII").unwrap();
assert!(matches!(
record.write_fastq(&mut out),
Err(Error::InvalidByte {
pos: 0,
byte: b'+',
..
})
));
assert!(Sequence::fastq("x", b"A+GT", b"IIII")
.unwrap()
.write_fastq(&mut out)
.is_ok());
let record = Sequence::fastq("x", b"ACGT", b"II\nI").unwrap();
assert!(matches!(
record.write_fastq(&mut out),
Err(Error::InvalidByte { byte: b'\n', .. })
));
assert!(Sequence::fastq("x", b"ACGT", b"@@@@")
.unwrap()
.write_fastq(&mut out)
.is_ok());
}
#[test]
fn quality_offset_conversion() {
let mut r = Sequence::fastq("r", b"ACGT", b"hhhh").unwrap();
r.convert_quality_offset(crate::qual::PHRED64, PHRED33);
assert_eq!(r.quality.as_deref(), Some(&b"IIII"[..]));
}
#[test]
fn clear_keeps_capacity() {
let mut r = Sequence::fastq("r", b"ACGT", b"IIII").unwrap();
let cap = r.seq.capacity();
r.clear();
assert!(r.id.is_empty() && r.seq.is_empty());
assert_eq!(r.quality.as_deref(), Some(&b""[..]));
assert_eq!(r.seq.capacity(), cap);
}
}