nuclease 0.2.0

Streaming FASTQ preprocessor with a focus on extensibility
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

//! Quality-based trimming transforms.

use crate::{
    plan::{
        BuildPlan, IntoExecutionStep, ReadTransform, TransformArena, TransformResult, TransformStep,
    },
    record::RecordView,
};

/// Read transform that applies 3' quality trimming using a Phred cutoff.
pub(crate) struct QualityTrim {
    cutoff: u8,
}

impl QualityTrim {
    /// Construct a new quality trimmer with the provided Phred cutoff.
    pub(crate) const fn new(cutoff: u8) -> Self {
        Self { cutoff }
    }

    /// Return the exclusive trim end selected by the 3' quality-trimming algorithm.
    fn quality_trim_end(&self, quality: &[u8]) -> usize {
        let mut best_sum = 0_i32;
        let mut running_sum = 0_i32;
        let mut trim_end = quality.len();

        for (idx, &ascii_q) in quality.iter().enumerate().rev() {
            let q = i32::from(ascii_q.saturating_sub(33));
            let score = i32::from(self.cutoff) - q;
            running_sum += score;

            if running_sum < 0 {
                running_sum = 0;
            } else if running_sum > best_sum {
                best_sum = running_sum;
                trim_end = idx;
            }
        }

        trim_end
    }
}

impl ReadTransform for QualityTrim {
    fn code(&self) -> &'static str {
        "quality_trim"
    }

    fn apply<'a>(&self, record: RecordView<'a>, _arena: &'a TransformArena) -> TransformResult<'a> {
        assert_eq!(
            record.sequence().len(),
            record.quality().len(),
            "quality trimming requires equal sequence and quality lengths"
        );

        let trim_end = self.quality_trim_end(record.quality());

        if trim_end == record.sequence().len() {
            return TransformResult {
                record,
                applied: false,
            };
        }

        TransformResult {
            record: record
                .with_sequence_and_quality(
                    &record.sequence()[..trim_end],
                    &record.quality()[..trim_end],
                )
                .expect("quality trimming should preserve equal sequence and quality lengths"),
            applied: true,
        }
    }
}

impl IntoExecutionStep for QualityTrim {
    fn into_execution_step(self) -> Box<dyn crate::plan::ExecutionStep> {
        Box::new(TransformStep(self))
    }
}

/// Fluent extension trait adding the `.quality_trim(...)` transform combinator to plans.
pub(crate) trait QualityTrimTransform: BuildPlan {
    /// Trim low-quality 3' sequence using the configured Phred cutoff.
    fn quality_trim(self, cutoff: u8) -> Self {
        self.step(QualityTrim::new(cutoff))
    }
}

impl<T> QualityTrimTransform for T where T: BuildPlan {}

#[cfg(test)]
mod tests {
    use crate::{
        plan::{ReadTransform, TransformArena},
        record::RecordView,
    };

    use super::QualityTrim;

    fn record(sequence: &'static [u8], quality: &'static [u8]) -> RecordView<'static> {
        RecordView::new(b"read1", sequence, quality)
    }

    #[test]
    fn quality_trim_leaves_high_quality_read_unchanged() {
        let trimmer = QualityTrim::new(20);
        let arena = TransformArena::new();
        let result = trimmer.apply(record(b"ACGT", b"IIII"), &arena);

        assert!(!result.applied);
        assert_eq!(result.record.sequence(), b"ACGT");
    }

    #[test]
    fn quality_trim_removes_low_quality_suffix() {
        let trimmer = QualityTrim::new(20);
        let arena = TransformArena::new();
        let result = trimmer.apply(record(b"ACGTAC", b"IIIII!"), &arena);

        assert!(result.applied);
        assert_eq!(result.record.sequence(), b"ACGTA");
        assert_eq!(result.record.quality(), b"IIIII");
    }

    #[test]
    fn quality_trim_cuts_at_low_quality_tail_start_even_if_tail_recovers() {
        let trimmer = QualityTrim::new(20);
        let arena = TransformArena::new();
        let result = trimmer.apply(record(b"ACGT", b"I!II"), &arena);

        assert!(result.applied);
        assert_eq!(result.record.sequence(), b"A");
        assert_eq!(result.record.quality(), b"I");
    }

    #[test]
    fn quality_trim_can_trim_entire_read_without_rejecting() {
        let trimmer = QualityTrim::new(20);
        let arena = TransformArena::new();
        let result = trimmer.apply(record(b"ACGT", b"!!!!"), &arena);

        assert!(result.applied);
        assert_eq!(result.record.sequence(), b"");
        assert_eq!(result.record.quality(), b"");
    }
}