Dev bqtools 0.5.10

[noamteyssier]

No description provided.

[gemini-code-assist]

Code Review

This pull request introduces a new FastQC-inspired quality control command (bqtools qc) to the CLI, implementing parallelized modules for base/sequence quality, GC content, sequence length, duplication levels, and overrepresented sequences. It also updates some internal collections to use hashbrown and refactors grep pattern matching. The review feedback identifies several key issues: a bug in duplication sampling when using --span due to absolute record indexing, potential OOM risks in sequence length tracking using dense vectors, unidiomatic Rust patterns (such as .div() and bool::then for side effects), a potential panic from slicing non-UTF-8 character boundaries, and redundant computations of overrepresented sequences.

[gemini-code-assist]

When processing a specific range of records (e.g., using the --span CLI option), record.index() returns the absolute 0-based index of the record in the file, not the relative index within the processed span.

If the span starts at or after self.sample_size (which defaults to 100,000), record.index() as usize >= self.sample_size will be true for all records in the span, causing the duplication and overrepresented sequence modules to sample zero records and produce empty reports.

To fix this, we should pass the starting index of the processed range/span to the duplication module (e.g., via QcConfig::build_qc_modules(start_index)) and check against start_index + sample_size.

[gemini-code-assist]

Using a dense Vec<usize> for SeqLenHistogram can lead to extremely high memory usage or Out-Of-Memory (OOM) crashes if the input BINSEQ file contains very long sequences (e.g., assembly scaffolds or long reads) or if a record has a corrupted/huge length.

Using std::collections::BTreeMap<usize, usize> is much more robust, memory-efficient for sparse/large lengths, and naturally keeps the lengths sorted for serialization without needing to filter out zero counts.

#[derive(Debug, Clone, Default)]
pub struct SeqLenHistogram { 
    /// Map from sequence length to count
    inner: std::collections::BTreeMap<usize, usize>,
}
impl SeqLenHistogram {
    fn is_empty(&self) -> bool {
        self.inner.is_empty()
    }
    /// Track a single read's length
    fn push(&mut self, len: usize) {
        if len == 0 {
            return;
        }
        *self.inner.entry(len).or_insert(0) += 1;
    }
    fn ingest(&mut self, other: &mut Self) {
        for (&len, &count) in &other.inner {
            *self.inner.entry(len).or_insert(0) += count;
        }
        other.inner.clear();
    }
    fn serialize_to<W: Write>(&self, wtr: &mut W) -> Result<()> {
        if self.is_empty() {
            return Ok(());
        }

        let mut ser = csv::WriterBuilder::default()
            .delimiter(b'\t')
            .has_headers(true)
            .from_writer(wtr);

        self.inner
            .iter()
            .try_for_each(|(&len, &count)| -> Result<()> {
                ser.serialize(&SeqLenRecord { len, count })
                    .map_err(Into::into)
            })?;

        ser.flush().map_err(Into::into)
    }

    fn total(&self) -> usize {
        self.inner.values().sum()
    }

    fn min_len(&self) -> Option<usize> {
        self.inner.keys().next().copied()
    }

    fn max_len(&self) -> Option<usize> {
        self.inner.keys().next_back().copied()
    }

    fn mean(&self) -> f64 {
        let total = self.total();
        if total == 0 {
            0.0
        } else {
            let sum: usize = self.inner.iter().map(|(&len, &c)| len * c).sum();
            sum as f64 / total as f64
        }
    }

    fn mode(&self) -> usize {
        self.inner
            .iter()
            .max_by_key(|&(_, &c)| c)
            .map_or(0, |(&len, _)| len)
    }

    fn summary_table(&self) -> Option<String> {
        if self.is_empty() {
            return None;
        }
        Some(table(
            &["Metric", "Value"],
            &[
                vec!["Reads".into(), self.total().to_string()],
                vec!["Min Length".into(), self.min_len().unwrap_or(0).to_string()],
                vec!["Max Length".into(), self.max_len().unwrap_or(0).to_string()],
                vec!["Mean Length".into(), format!("{:.2}", self.mean())],
                vec!["Mode Length".into(), self.mode().to_string()],
            ],
        ))
    } 
}

[gemini-code-assist]

Using the .div() method from std::ops::Div is unidiomatic in Rust when standard division / can be used directly.

Suggested change

	let binned_mean = (total as f64).div(&(qual.len() as f64)).round() as usize;
	let binned_mean = (total as f64 / qual.len() as f64).round() as usize;

[gemini-code-assist]

Remove the unused ops::Div import since standard division / is preferred.

Suggested change

	use std::{io::Write, ops::Div, path::Path, sync::Arc};
	use std::{io::Write, path::Path, sync::Arc};

[gemini-code-assist]

Using bool::then purely for side effects (calling add_module) is unidiomatic in Rust. Standard if statements are much more readable and avoid creating unnecessary closures.

        if self.per_base_qual {
            add_module(QcModuleType::new_base_quality());
        }
        if self.per_seq_qual {
            add_module(QcModuleType::new_seq_quality());
        }
        if self.per_base_content {
            add_module(QcModuleType::new_base_content());
        }
        if self.per_seq_gc {
            add_module(QcModuleType::new_gc_content());
        }
        if self.seq_length {
            add_module(QcModuleType::new_seq_length());
        }

[gemini-code-assist]

Slicing a str directly by a byte index (&seq[..MAX_LEN]) can panic if the index does not fall on a UTF-8 character boundary. While nucleotide sequences are typically ASCII, if there are any invalid UTF-8 bytes, String::from_utf8_lossy will insert the 3-byte replacement character ``, which could cause a panic if sliced in the middle.

It is safer to find the nearest valid character boundary or use character iterators.

Suggested change

	fn truncate_sequence(seq: &[u8]) -> String {
	const MAX_LEN: usize = 40;
	let seq = String::from_utf8_lossy(seq);
	if seq.len() > MAX_LEN {
	format!("{}...", &seq[..MAX_LEN])
	} else {
	seq.into_owned()
	}
	}
	fn truncate_sequence(seq: &[u8]) -> String {
	const MAX_LEN: usize = 40;
	let seq = String::from_utf8_lossy(seq);
	if seq.len() > MAX_LEN {
	let mut end = MAX_LEN;
	while end > 0 && !seq.is_char_boundary(end) {
	end -= 1;
	}
	format!("{}...", &seq[..end])
	} else {
	seq.into_owned()
	}
	}

[gemini-code-assist]

counter.overrepresented(...) is an expensive operation because it iterates over the entire HashMap, collects the elements, and sorts them. Calling it twice (once to check if it's empty, and once inside serialize_overrepresented_to) is inefficient.

We should compute it once and reuse the result.

            if self.emit_overrepresented {
                let overrepresented = counter.overrepresented(self.overrepresented_threshold);
                if overrepresented.is_empty() {
                    trace!(
                        "No {label} sequences met the overrepresented threshold ({}%)",
                        self.overrepresented_threshold
                    );
                } else {
                    let mut handle = match_output(Some(outdir.as_ref().join(overrep_path)))?;
                    let mut ser = csv::WriterBuilder::default()
                        .delimiter(b'\t')
                        .has_headers(true)
                        .from_writer(&mut handle);

                    overrepresented
                        .into_iter()
                        .try_for_each(|(seq, count, pct)| -> Result<()> {
                            ser.serialize(&OverrepresentedRecord {
                                sequence: String::from_utf8_lossy(seq).into_owned(),
                                count,
                                pct,
                            })
                            .map_err(Into::into)
                        })?;

                    ser.flush()?;
                }
            }