seqtk_rs/
trim.rs

1use crate::io_utils::{FqReader, FxWriter};
2/// Trims low-quality bases from a FASTQ data based on a quality threshold `Q`.
3/// Outputs the FASTQ to [`std::io::stdout()`].
4///
5/// The algorithm:
6///
7/// (1) Scan the read from 5' to 3' to find the first base with `quality >= Q`.
8///     This is the starting position of the trimmed read.
9/// (2) Compute a running score from 5' to 3'
10///        `score(i) = score(i-1) + quality(i) - Q`
11///      where the initial score is `score(start) = quality(start) - Q`.
12/// (3) Find out the maximun score. This is the end of the trimmed read.
13///
14/// Note:
15/// (1) If all base quality are less than `Q`, the read is discarded.
16/// (2) If the trimmed read is shorter than `min_len`, we first extend it from the 3' end. If still too short, extend form the 5' end until `min_len` is reached or no more bases are available.
17///
18///
19/// # Arguments
20///
21/// * `path` - FASTQ path
22/// * `q_plus_ascii` - The sum of quality threshold and asciibase.
23/// * `minilen` - The minimum length of read.
24///
25/// # Errors
26///
27/// Return an error if the operation cannot be completed.
28///
29/// # Notes
30/// Quality trimming is no longer necessary in most modern sequencing pipelines.
31/// Its usefulness depends on the sequencing technology and the goals of your downstream analysis.
32///
33/// For example:
34/// * Illumina reads may benefit from light trimming low-quality tails at 3' end.
35/// * Long-read technologies are generally not necessary to do quality trimming, because it may remove informative regions and reduce read length unnecessarily.
36///
37pub fn trimfq(fq_path: &str, q_plus_ascii: u8, minlen: usize) -> Result<(), std::io::Error> {
38    let reader = FqReader::new(fq_path)?;
39    let mut writer = FxWriter::new(false);
40    let (mut start, mut end, mut size): (usize, usize, usize);
41    for record in reader.records() {
42        let read = record.unwrap();
43        (start, end) = trim_read_by_q(read.qual(), q_plus_ascii);
44
45        // discard read if start == end (Q < q_threshold for all bases)
46        if start < end {
47            size = read.qual().len();
48
49            if size < minlen {
50                // write full length is size < minlen
51                (start, end) = (0, size - 1);
52            } else if minlen > (end - start + 1) {
53                if size - start >= minlen {
54                    // append 3' seq to reach the minlen
55                    end = start + minlen - 1;
56                } else {
57                    // append 5' seq if the whole 3' end is not enough
58                    end = size - 1;
59                    start = size - minlen;
60                }
61            }
62            writer.write(
63                read.id(),
64                &read.seq()[start..(end + 1)],
65                read.desc(),
66                &read.qual()[start..(end + 1)],
67            )?;
68        }
69    }
70    Ok(())
71}
72
73/// return start, end
74fn trim_read_by_q(qual: &[u8], q_plus_ascii: u8) -> (usize, usize) {
75    let mut start: usize = 0;
76    let len = qual.len();
77    for &q in qual {
78        if q >= q_plus_ascii {
79            break;
80        } else {
81            start += 1;
82        }
83    }
84    if start == len {
85        return (len, len);
86    }
87
88    let qthd_add_ascii_i = q_plus_ascii as isize;
89    let mut err_sum: Vec<isize> = vec![0; len];
90
91    err_sum[start] = qual[start] as isize - qthd_add_ascii_i;
92    let (mut max_idx, mut max_val) = (start, err_sum[start]);
93
94    ((start + 1)..len).for_each(|i| {
95        err_sum[i] = err_sum[i - 1] + qual[i] as isize - qthd_add_ascii_i;
96        if err_sum[i] >= max_val {
97            max_idx = i;
98            max_val = err_sum[i];
99        }
100    });
101
102    // println!("[!] thres: {} sum: {:?}", qthd_add_ascii_i, &err_sum[start..]);
103
104    (start, max_idx)
105}
106#[cfg(test)]
107mod tests {
108    use super::*;
109
110    #[test]
111    fn test_trim_read() {
112        let qual = b"&&&&&...&...++.&.++...&...''..'.&..+..++++.++";
113        let (start, end) = trim_read_by_q(qual, 33 + 10);
114        assert_eq!(end - start + 1, 40);
115        let qual= b"?@<DDD;2A<><FIBHB?FCGAHHEBEHAFCBEFGGB@4CCA@?*?DH9?B?BDC/?81.8BFFEHG@>@G=DGCECEHEHF>>;?;?B=3@CC:(,(98',>5(82)<5>>223@4+4>@+5855>>>@:AA@>:43>@##########";
116
117        let (start, end) = trim_read_by_q(qual, 33 + 20);
118        assert_eq!(end - start + 1, 140);
119    }
120}
seqtk_rs/trim.rs

seqtk_rs/
trim.rs
