FASTA and FASTQ

use fasta.qweruse fasta_parse, fasta_write from fasta      // Quillo 0.6: only these

serde_json = { version = "1", features = ["preserve_order"] }

// A small sequence bench: FASTA and FASTQ read, GC content, the other strand, a protein, reading// frames, qualities and k-mers, on sequences kept here so it runs anywhere.////   quillo run fasta.qwe//   then open http://localhost:8080/fasta, /sequence?s=ATGGCC, /fastq, /orfs, /kmers and /refuseduse "fasta.qwer"let sample = ">sp|P01308 insulin fragment; a comment lineATGGCCCTGTGGATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCAGCC>controlACGTNNACGT"route GET /fasta {    let records = fasta_parse(sample)?    {        records: records,        back: fasta_write(records, 20)?,    }}route GET /sequence ? { s: str = "ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG" } {    {        gc: fasta_gc(s),        rc: fasta_reverse_complement(s)?,        rna: fasta_transcribe(s),        protein: fasta_translate(s, {})?,        to_stop: fasta_translate(s, { to_stop: true })?,        frame2: fasta_translate(s, { frame: 2 })?,        minus: fasta_translate(s, { frame: -1 })?,    }}// A stop in the standard code that codes tryptophan in mitochondria, and codons left open by N.route GET /codes = {    standard: fasta_translate("ATGTGAAGA", {})?,    mitochondrial: fasta_translate("ATGTGAAGA", { table: 2 })?,    open_codon: fasta_translate("GCNTTNATH", {})?,    iupac: fasta_reverse_complement("ACGTRYSWKMBDHVNacgt")?,    rna: fasta_reverse_complement("AUGGC")?,}let reads = "@read1 lane 1GATTTGGGGTTCAAAGCAGT+!''*((((***+))%%%++)@read2ACGT+read2IIII"route GET /fastq {    let f = fastq_parse(reads, {})?    {        encoding: f.encoding,        reads: f.records,        scores: fastq_scores("II5#", 33)?,        back: fastq_write([{ id: "x", sequence: "ACGT", scores: [40, 30, 20, 2] }])?,        trimmed: fastq_trim("ACGTACGTAC", "IIIIIII###", 20, 3)?,    }}route GET /orfs = fasta_orfs("CCATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAGCC", 3, 1)?route GET /kmers = {    plain: fasta_kmers("ACGTACGTAC", 4, { limit: 3 })?,    canonical: fasta_kmers("ACGTACGT", 3, { canonical: true })?,    skipped: fasta_kmers("ACGNACG", 3, {})?.total,}route GET /refused = {    headless: fasta_parse("ACGT\n").error,    bad_base: fasta_parse(">x\nAC GT5\n").error,    short_quality: fastq_parse("@r\nACGT\n+\nII\n", {}).error,    no_plus: fastq_parse("@r\nACGT\n", {}).error,    wrong_offset: fastq_parse("@r\nAC\n+\n!!\n", { encoding: "phred64" }).error,    table: fasta_translate("ATG", { table: 99 }).error,    k: fasta_kmers("ACGT", 0, {}).error,    letter: fasta_reverse_complement("ACGX").error,}

//! FASTA and FASTQ read and written: sequences with their ids, GC content, reverse complements with IUPAC codes, translation by genetic code, open reading frames, Phred quality scores and trimming, and k-mer counts.//! version: 1.0.0//!//! [dependencies]//! serde_json = { version = "1", features = ["preserve_order"] }use serde_json::{json, Map, Value};use std::collections::HashMap;const MAX: usize = 256 * 1024 * 1024;fn head(line: &str) -> (String, String) {    let t = line[1..].trim();    match t.split_once(|c: char| c.is_whitespace()) {        Some((id, d)) => (id.to_string(), d.trim().to_string()),        None => (t.to_string(), String::new()),    }}fn valid_seq(s: &str, line: usize) -> Result<(), String> {    if let Some(c) = s.chars().find(|c| !(c.is_ascii_alphabetic() || *c == '*' || *c == '-' || *c == '.')) {        return Err(format!("line {}: `{}` cannot be in a sequence", line, c));    }    Ok(())}/// FASTA text read: [{ id, description, sequence, length }]. Lines are joined; `;` comment lines are skipped; a sequence before the first `>` or a bad character is an error value naming the line.pub fn fasta_parse(text: &str) -> Result<Value, String> {    if text.len() > MAX {        return Err("that FASTA is over 256 MB".into());    }    let mut out: Vec<(String, String, String)> = vec![];    for (i, raw) in text.lines().enumerate() {        let line = raw.trim_end();        if line.is_empty() || line.starts_with(';') {            continue;        }        if line.starts_with('>') {            let (id, d) = head(line);            if id.is_empty() {                return Err(format!("line {}: a `>` with no id", i + 1));            }            out.push((id, d, String::new()));        } else {            let t = line.trim();            valid_seq(t, i + 1)?;            match out.last_mut() {                Some(r) => r.2.push_str(t),                None => return Err(format!("line {}: sequence before the first `>` header", i + 1)),            }        }    }    Ok(json!(out.iter().map(|(id, d, s)| json!({ "id": id, "description": d, "sequence": s, "length": s.len() })).collect::<Vec<_>>()))}/// Records as FASTA, sequence lines wrapped at `width` (0 for one long line; 60 and 70 are usual).pub fn fasta_write(records: serde_json::Value, width: i64) -> Result<String, String> {    let list = records.as_array().ok_or("give a list of { id, description, sequence }")?;    let mut o = String::new();    for (i, r) in list.iter().enumerate() {        let id = r.get("id").and_then(|v| v.as_str()).ok_or(format!("record {} has no id", i + 1))?;        let seq = r.get("sequence").and_then(|v| v.as_str()).ok_or(format!("record {} has no sequence", i + 1))?;        o.push('>');        o.push_str(id);        if let Some(d) = r.get("description").and_then(|v| v.as_str()).filter(|d| !d.is_empty()) {            o.push(' ');            o.push_str(d);        }        o.push('\n');        if width <= 0 {            o.push_str(seq);            o.push('\n');        } else {            for c in seq.as_bytes().chunks(width as usize) {                o.push_str(&String::from_utf8_lossy(c));                o.push('\n');            }        }    }    Ok(o)}// Phred+33 runs from `!`; Phred+64 (old Illumina) from `@`. Below `;` only +33 is possible, above// `J` or so only +64.fn guess_offset(lows: u8, highs: u8) -> i64 {    if lows < 59 { 33 } else if highs > 74 { 64 } else { 33 }}fn offset_of(options: &Value) -> Result<Option<i64>, String> {    match options.get("encoding").and_then(|v| v.as_str()) {        None | Some("auto") => Ok(None),        Some("phred33") | Some("sanger") | Some("illumina1.8") => Ok(Some(33)),        Some("phred64") | Some("illumina1.3") => Ok(Some(64)),        Some(o) => Err(format!("no encoding `{}`; there are phred33, phred64 and auto", o)),    }}/// FASTQ text read: { encoding, records: [{ id, description, sequence, quality, length, mean_quality }] }. Takes wrapped reads. A quality line of the wrong length names its record. `options.encoding` is phred33, phred64 or auto, guessed from the characters.pub fn fastq_parse(text: &str, options: serde_json::Value) -> Result<Value, String> {    let want = offset_of(&options)?;    let lines: Vec<&str> = text.lines().collect();    let mut i = 0;    let mut recs: Vec<(String, String, String, String)> = vec![];    let (mut lo, mut hi) = (255u8, 0u8);    while i < lines.len() {        let h = lines[i].trim_end();        if h.is_empty() {            i += 1;            continue;        }        if !h.starts_with('@') {            return Err(format!("line {}: a FASTQ record starts with `@`", i + 1));        }        let at = i + 1;        let (id, d) = head(h);        i += 1;        let mut seq = String::new();        while i < lines.len() && !lines[i].starts_with('+') {            let t = lines[i].trim();            valid_seq(t, i + 1)?;            seq.push_str(t);            i += 1;        }        if i >= lines.len() {            return Err(format!("line {}: the record `{}` has no `+` line", at, id));        }        let plus = lines[i].trim_end();        if plus.len() > 1 && plus[1..] != h[1..] {            return Err(format!("line {}: the `+` line repeats a different id", i + 1));        }        i += 1;        let mut q = String::new();        // The quality may be wrapped and may begin with `@`, so it is read by length.        while i < lines.len() && q.len() < seq.len() {            q.push_str(lines[i].trim_end());            i += 1;        }        if q.len() != seq.len() {            return Err(format!("record `{}`: {} quality scores for {} bases", id, q.len(), seq.len()));        }        for b in q.bytes() {            lo = lo.min(b);            hi = hi.max(b);            if !(33..=126).contains(&b) {                return Err(format!("record `{}`: byte {} is not a quality character", id, b));            }        }        recs.push((id, d, seq, q));    }    let off = want.unwrap_or_else(|| guess_offset(lo, hi));    if off == 64 && lo < 59 {        return Err("the qualities go below `;`, so they cannot be Phred+64".into());    }    let records: Vec<Value> = recs        .iter()        .map(|(id, d, s, q)| {            let scores: Vec<i64> = q.bytes().map(|b| b as i64 - off).collect();            json!({ "id": id, "description": d, "sequence": s, "quality": q, "length": s.len(), "mean_quality": mean_q(&scores) })        })        .collect();    Ok(json!({ "encoding": if off == 33 { "phred33" } else { "phred64" }, "records": records }))}// Average of error probabilities, as an average of Phred scores would flatter a read.fn mean_q(q: &[i64]) -> f64 {    if q.is_empty() {        return 0.0;    }    let p: f64 = q.iter().map(|x| 10f64.powf(-*x as f64 / 10.0)).sum::<f64>() / q.len() as f64;    (-10.0 * p.log10() * 100.0).round() / 100.0}/// Records as FASTQ, with quality text as read, or `scores` as a list of ints written in Phred+33.pub fn fastq_write(records: serde_json::Value) -> Result<String, String> {    let list = records.as_array().ok_or("give a list of { id, sequence, quality }")?;    let mut o = String::new();    for (i, r) in list.iter().enumerate() {        let id = r.get("id").and_then(|v| v.as_str()).ok_or(format!("record {} has no id", i + 1))?;        let seq = r.get("sequence").and_then(|v| v.as_str()).ok_or(format!("record {} has no sequence", i + 1))?;        let q = match (r.get("quality").and_then(|v| v.as_str()), r.get("scores").and_then(|v| v.as_array())) {            (Some(q), _) => q.to_string(),            (None, Some(s)) => s.iter().map(|x| char::from_u32((x.as_i64().unwrap_or(0).clamp(0, 93) + 33) as u32).unwrap_or('!')).collect(),            _ => return Err(format!("record `{}` has neither quality nor scores", id)),        };        if q.len() != seq.len() {            return Err(format!("record `{}`: {} quality scores for {} bases", id, q.len(), seq.len()));        }        let d = r.get("description").and_then(|v| v.as_str()).filter(|d| !d.is_empty()).map(|d| format!(" {}", d)).unwrap_or_default();        o.push_str(&format!("@{}{}\n{}\n+\n{}\n", id, d, seq, q));    }    Ok(o)}/// A quality string as scores: { scores, mean, error_rate, q20, q30 }, with q20 and q30 the share of bases at or above them. `offset` is 33 or 64.pub fn fastq_scores(quality: &str, offset: i64) -> Result<Value, String> {    if offset != 33 && offset != 64 {        return Err("offset is 33 or 64".into());    }    let q: Vec<i64> = quality.bytes().map(|b| b as i64 - offset).collect();    if let Some(bad) = q.iter().find(|x| **x < 0 || **x > 93) {        return Err(format!("{} is no Phred score: is the offset right?", bad));    }    let share = |t: i64| ((q.iter().filter(|x| **x >= t).count() as f64 / q.len().max(1) as f64) * 10000.0).round() / 10000.0;    let m = mean_q(&q);    Ok(json!({ "scores": q, "mean": m, "error_rate": (10f64.powf(-m / 10.0) * 1e6).round() / 1e6, "q20": share(20), "q30": share(30) }))}/// A read cut where it turns bad, as Trimmomatic's SLIDINGWINDOW does: a window of `window` bases whose mean falls under `threshold` ends the read, and low bases at the end are dropped. Gives { sequence, quality, removed }.pub fn fastq_trim(sequence: &str, quality: &str, threshold: i64, window: i64) -> Result<Value, String> {    if sequence.len() != quality.len() {        return Err(format!("{} quality scores for {} bases", quality.len(), sequence.len()));    }    if window < 1 {        return Err("the window is at least 1 base".into());    }    let w = window as usize;    let q: Vec<i64> = quality.bytes().map(|b| b as i64 - 33).collect();    let mut keep = q.len();    if q.len() >= w {        for i in 0..=q.len() - w {            let sum: i64 = q[i..i + w].iter().sum();            if sum < threshold * w as i64 {                keep = i;                break;            }        }    }    if keep < q.len() {        // The window's good start stays.        let mut end = keep;        while end < (keep + w).min(q.len()) && q[end] >= threshold {            end += 1;        }        keep = end;    }    Ok(json!({ "sequence": &sequence[..keep], "quality": &quality[..keep], "removed": q.len() - keep }))}/// GC content, N and gaps not counted: { gc, at, length, counted, counts, gc_skew }. RNA's U counts as T, and IUPAC S (G or C) as GC, W as AT.pub fn fasta_gc(sequence: &str) -> Value {    let mut c: HashMap<char, usize> = HashMap::new();    for ch in sequence.chars() {        *c.entry(ch.to_ascii_uppercase()).or_insert(0) += 1;    }    let n = |x: char| *c.get(&x).unwrap_or(&0);    let gc = n('G') + n('C') + n('S');    let at = n('A') + n('T') + n('U') + n('W');    let counted = gc + at;    let frac = |x: usize| if counted == 0 { 0.0 } else { (x as f64 / counted as f64 * 10000.0).round() / 10000.0 };    let skew = if n('G') + n('C') == 0 { 0.0 } else { ((n('G') as f64 - n('C') as f64) / (n('G') + n('C')) as f64 * 10000.0).round() / 10000.0 };    let mut keys: Vec<&char> = c.keys().collect();    keys.sort();    let counts: Map<String, Value> = keys.iter().map(|k| (k.to_string(), json!(c[k]))).collect();    json!({ "gc": frac(gc), "at": frac(at), "length": sequence.chars().count(), "counted": counted, "counts": counts, "gc_skew": skew })}fn comp(c: char, rna: bool) -> Result<char, String> {    let up = c.to_ascii_uppercase();    let o = match up {        'A' => if rna { 'U' } else { 'T' },        'T' | 'U' => 'A', 'G' => 'C', 'C' => 'G', 'R' => 'Y', 'Y' => 'R', 'S' => 'S', 'W' => 'W', 'K' => 'M', 'M' => 'K',        'B' => 'V', 'V' => 'B', 'D' => 'H', 'H' => 'D', 'N' => 'N', '-' => '-', '.' => '.',        _ => return Err(format!("`{}` is not a nucleotide code", c)),    };    Ok(if c.is_ascii_lowercase() { o.to_ascii_lowercase() } else { o })}fn revcomp(s: &str) -> Result<String, String> {    let rna = s.contains('U') || s.contains('u');    s.chars().rev().map(|c| comp(c, rna)).collect()}/// The reverse complement, with IUPAC ambiguity codes and case kept. A sequence with U comes back as RNA.pub fn fasta_reverse_complement(sequence: &str) -> Result<String, String> {    revcomp(sequence)}/// DNA to RNA: T becomes U, on the coding strand.pub fn fasta_transcribe(sequence: &str) -> String {    sequence.replace('T', "U").replace('t', "u")}const CODE1: &str = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG";const CODE2: &str = "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG";const CODE3: &str = "FFLLSSSSYY**CCWWTTTTPPPPHHQQRRRRIIMMTTTTNNKKSSRRVVVVAAAADDEEGGGG";const CODE4: &str = "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG";fn code_of(table: i64) -> Result<&'static [u8], String> {    Ok(match table {        1 | 11 => CODE1.as_bytes(),        2 => CODE2.as_bytes(),        3 => CODE3.as_bytes(),        4 => CODE4.as_bytes(),        o => return Err(format!("no genetic code {}; there are 1 (standard), 2 (vertebrate mitochondrial), 3 (yeast mitochondrial), 4 (mould, Mycoplasma) and 11 (bacterial)", o)),    })}fn bases(c: u8) -> &'static [usize] {    match c.to_ascii_uppercase() {        b'T' | b'U' => &[0], b'C' => &[1], b'A' => &[2], b'G' => &[3],        b'R' => &[2, 3], b'Y' => &[0, 1], b'S' => &[1, 3], b'W' => &[0, 2], b'K' => &[0, 3], b'M' => &[1, 2],        b'B' => &[0, 1, 3], b'D' => &[0, 2, 3], b'H' => &[0, 1, 2], b'V' => &[1, 2, 3], b'N' => &[0, 1, 2, 3],        _ => &[],    }}// An ambiguous codon still has one amino acid when all it may be agree: GCN is alanine.fn amino(codon: &[u8], code: &[u8]) -> u8 {    let mut found = 0u8;    let (a, b, c) = (bases(codon[0]), bases(codon[1]), bases(codon[2]));    if a.is_empty() || b.is_empty() || c.is_empty() {        return b'X';    }    for x in a {        for y in b {            for z in c {                let aa = code[x * 16 + y * 4 + z];                if found == 0 {                    found = aa;                } else if found != aa {                    return b'X';                }            }        }    }    found}/// A coding sequence as protein. `options`: { table (1, 2, 3, 4, 11), frame (1 to 3, -1 to -3), to_stop }. Stops are `*`; `to_stop` ends at the first. A codon its ambiguity leaves open is `X`.pub fn fasta_translate(sequence: &str, options: serde_json::Value) -> Result<String, String> {    let code = code_of(options.get("table").and_then(|v| v.as_i64()).unwrap_or(1))?;    let frame = options.get("frame").and_then(|v| v.as_i64()).unwrap_or(1);    if !matches!(frame, 1 | 2 | 3 | -1 | -2 | -3) {        return Err("frame is 1 to 3, or -1 to -3 for the other strand".into());    }    let clean: String = sequence.chars().filter(|c| !c.is_whitespace() && *c != '-').collect();    let strand = if frame < 0 { revcomp(&clean)? } else { clean };    let b = strand.as_bytes();    let start = (frame.unsigned_abs() as usize - 1).min(b.len());    let stop = options.get("to_stop").and_then(|v| v.as_bool()).unwrap_or(false);    let mut out = String::new();    for c in b[start..].chunks_exact(3) {        let a = amino(c, code);        if stop && a == b'*' {            break;        }        out.push(a as char);    }    Ok(out)}/// Open reading frames on both strands from ATG to a stop, at least `min_length` amino acids: [{ strand, frame, start, end, length, protein }], longest first. `start` and `end` count from 1 on the forward strand.pub fn fasta_orfs(sequence: &str, min_length: i64, table: i64) -> Result<Value, String> {    let code = code_of(table)?;    let fwd: String = sequence.to_uppercase().chars().filter(|c| !c.is_whitespace()).collect::<String>().replace('U', "T");    let rev = revcomp(&fwd)?;    let n = fwd.len();    let mut out: Vec<(usize, Value)> = vec![];    for (strand, s) in [(1i64, fwd.as_bytes()), (-1, rev.as_bytes())] {        for frame in 0..3usize {            let mut i = frame;            while i + 3 <= s.len() {                if &s[i..i + 3] == b"ATG" {                    let mut protein = String::new();                    let mut j = i;                    let mut closed = false;                    while j + 3 <= s.len() {                        let a = amino(&s[j..j + 3], code);                        j += 3;                        if a == b'*' {                            closed = true;                            break;                        }                        protein.push(a as char);                    }                    if closed && protein.len() as i64 >= min_length {                        let (st, en) = if strand == 1 { (i + 1, j) } else { (n - j + 1, n - i) };                        out.push((protein.len(), json!({ "strand": strand, "frame": frame as i64 + 1, "start": st, "end": en, "length": protein.len(), "protein": protein })));                    }                    i = j;                } else {                    i += 3;                }            }        }    }    out.sort_by(|a, b| b.0.cmp(&a.0));    Ok(json!(out.into_iter().map(|o| o.1).collect::<Vec<_>>()))}/// K-mer counts, commonest first, as [{ kmer, count }], capped at `limit`. `options`: { canonical (a k-mer and its reverse complement count as one), limit }. A k-mer with N is skipped. Gives { k, total, distinct, kmers }.pub fn fasta_kmers(sequence: &str, k: i64, options: serde_json::Value) -> Result<Value, String> {    if !(1..=64).contains(&k) {        return Err("k is 1 to 64".into());    }    let canonical = options.get("canonical").and_then(|v| v.as_bool()).unwrap_or(false);    let limit = options.get("limit").and_then(|v| v.as_i64()).unwrap_or(100).max(0) as usize;    let s = sequence.to_uppercase();    let b = s.as_bytes();    let mut c: HashMap<&[u8], u64> = HashMap::new();    let mut extra: HashMap<String, u64> = HashMap::new();    let mut total = 0u64;    for w in b.windows(k as usize) {        if !w.iter().all(|x| matches!(x, b'A' | b'C' | b'G' | b'T' | b'U')) {            continue;        }        total += 1;        if canonical {            let t = std::str::from_utf8(w).unwrap_or_default();            let r = revcomp(t).unwrap_or_default();            *extra.entry(if r.as_str() < t { r } else { t.to_string() }).or_insert(0) += 1;        } else {            *c.entry(w).or_insert(0) += 1;        }    }    let mut rows: Vec<(String, u64)> = if canonical { extra.into_iter().collect() } else { c.into_iter().map(|(k, v)| (String::from_utf8_lossy(k).into_owned(), v)).collect() };    rows.sort_by(|a, b| b.1.cmp(&a.1).then(a.0.cmp(&b.0)));    let distinct = rows.len();    rows.truncate(limit);    Ok(json!({ "k": k, "total": total, "distinct": distinct, "kmers": rows.iter().map(|(k, c)| json!({ "kmer": k, "count": c })).collect::<Vec<_>>() }))}