use crate::commands::command::Command;
use crate::commands::utils::{
BUFFER_SIZE, OwnedRecordIter, PairingRule, SplitNameCheck, aggregate_errors, check_at_most_two,
check_dash_at_most_once, check_distinct_inputs, default_dash, fmt_count, open_fastq_inputs,
pull_pair_interleaved, resolve_inputs, resolve_real_path, sniff_single_input,
};
use anyhow::{Result, anyhow};
use bgzf::{CompressionLevel, Compressor};
use chelae_lib::IUPAC_MASKS;
use chelae_lib::adapter_db::{KitAdapter, expand_kit_name};
use clap::{Parser, ValueEnum};
use crossbeam_channel::{Receiver, Sender, bounded};
use fgoxide::io::{DelimFile, Io};
use fgoxide::iter::IntoChunkedReadAheadIterator;
use log::{info, warn};
use read_structure::{ReadStructure, ReadStructureError, SegmentType, SkipHandling};
use seq_io::fastq::OwnedRecord;
use seq_io::fastq::Record;
use serde::{Deserialize, Serialize};
use std::fs::File;
use std::io::{BufRead, BufWriter, IsTerminal, Write};
use std::path::{Path, PathBuf};
use std::str::FromStr;
use std::sync::atomic::{AtomicBool, Ordering};
use std::thread;
use wide::{i8x16, u8x16, u8x32};
const LOG_EVERY: u64 = 5_000_000;
const DEFAULT_COMPRESSION_LEVEL: usize = 1;
const UMI_JOIN: u8 = b'-';
const UMI_ID_SEP: u8 = b':';
const MAX_READ_ID_FIELDS: usize = 8;
const INSERT_STATS_UPDATE_INTERVAL: u64 = 1000;
const INSERT_STATS_MIN_DETECTIONS: u64 = 64;
const ADAPTER_EVIDENCE_PROBE_LEN: usize = 16;
const ADAPTER_EVIDENCE_MAX_MM: usize = 5;
const NEAR_WALK_SHIFTS: isize = 4;
#[derive(Parser, Debug)]
#[command(version)]
#[clap(verbatim_doc_comment)]
pub(crate) struct Trim {
#[clap(long, short = 'i', num_args = 1..=2)]
inputs: Vec<PathBuf>,
#[clap(long, short = 'o', num_args = 1..=2)]
outputs: Vec<PathBuf>,
#[clap(long, value_enum, default_value_t = OutputCompression::Auto)]
output_compression: OutputCompression,
#[clap(long, short = 't', default_value = "4")]
threads: usize,
#[clap(long, short = 'c', default_value_t = DEFAULT_COMPRESSION_LEVEL)]
compression_level: usize,
#[clap(long, short = 'm')]
metrics: Option<PathBuf>,
#[clap(long, short = 'r', num_args = 1..=2)]
read_structures: Vec<ReadStructure>,
#[clap(long)]
discard_unsupported_segments: bool,
#[clap(long, short = 'a', num_args = 1..=2)]
adapter_sequence: Vec<String>,
#[clap(long, short = 'f')]
adapter_fasta: Option<PathBuf>,
#[clap(long, short = 'k', num_args = 1..)]
kit: Vec<String>,
#[clap(long)]
no_overlap_detection: bool,
#[clap(long, default_value = "30")]
overlap_min_length: usize,
#[clap(long, default_value = "0.10")]
overlap_max_mismatch_rate: f64,
#[clap(long, default_value = "64")]
overlap_diagnostic_length: usize,
#[clap(long, hide = true, default_value_t = 1e-4)]
overlap_trust_max_chance: f64,
#[clap(long, default_value = "6")]
adapter_min_length: usize,
#[clap(long, default_value = "0.125")]
adapter_mismatch_rate: f64,
#[clap(long, default_value = "10")]
trim_polyg: usize,
#[clap(long, num_args = 0..=1, default_missing_value = "10")]
trim_polyx: Option<usize>,
#[clap(long, num_args = 0..=1, default_missing_value = "8:20")]
quality_trim_3p: Option<QualityTrim>,
#[clap(long, num_args = 0..=1, default_missing_value = "8:20")]
quality_trim_5p: Option<QualityTrim>,
#[clap(long = "filter-length", short = 'l', default_value = "15")]
filter_length: LengthFilter,
#[clap(long = "filter-max-ns")]
filter_max_ns: Option<usize>,
#[clap(long = "filter-mean-qual")]
filter_mean_qual: Option<u8>,
#[clap(long = "filter-low-qual")]
filter_low_qual: Option<LowQualFilter>,
#[clap(long, short = 'j')]
json: Option<PathBuf>,
#[clap(long, hide = true, default_value = "1024")]
batch_size: usize,
#[clap(long)]
expected_insert_size: Option<usize>,
#[clap(long, default_value_t = false)]
insert_size_stats: bool,
}
impl Trim {
fn validate(&self) -> Result<()> {
let mut errors: Vec<String> = Vec::new();
let inputs = default_dash(&self.inputs);
let outputs = default_dash(&self.outputs);
if let Err(e) = resolve_inputs(&self.inputs, std::io::stdin().is_terminal()) {
errors.push(e.to_string());
}
check_dash_at_most_once(&inputs, "--inputs", &mut errors);
check_dash_at_most_once(&outputs, "--outputs", &mut errors);
check_at_most_two(&inputs, "--inputs", &mut errors);
check_at_most_two(&outputs, "--outputs", &mut errors);
check_distinct_inputs(&inputs, &mut errors);
for path in &inputs {
if path.as_os_str() != "-" && !path.exists() {
errors.push(format!("Input file {path:?} does not exist."));
}
}
for path in &outputs {
if path.as_os_str() == "-" {
continue;
}
if let Some(parent) = path.parent()
&& !parent.as_os_str().is_empty()
&& !parent.exists()
{
errors.push(format!(
"Output parent directory {parent:?} does not exist (for {path:?})."
));
}
}
self.check_no_output_overwrites_input(&inputs, &outputs, &mut errors);
self.check_no_duplicate_outputs(&outputs, &mut errors);
self.check_metrics_json_not_dash(&mut errors);
self.check_read_structures(&mut errors);
self.check_adapter_args(&mut errors);
self.check_filter_args(&mut errors);
if self.threads < 1 {
errors.push(format!("Threads must be at least 1, got {}.", self.threads));
}
if !(1..=12).contains(&self.compression_level) {
errors.push(format!(
"Compression level must be in 1..=12, got {}.",
self.compression_level
));
}
aggregate_errors(errors)
}
fn validate_post_detection(&self, num_mates: usize, num_outputs: usize) -> Result<()> {
let mut errors: Vec<String> = Vec::new();
if num_outputs > num_mates {
errors.push(format!(
"{num_outputs} output(s) given but the input was detected as single-end (1 mate \
per record); single-end input supports at most 1 output."
));
}
if !self.read_structures.is_empty() && self.read_structures.len() != num_mates {
errors.push(format!(
"Number of read-structures ({}) must be 0 or equal to the number of mates ({num_mates}).",
self.read_structures.len(),
));
}
if self.adapter_sequence.len() > num_mates {
errors.push(format!(
"{} --adapter-sequence values supplied but input has only {num_mates} mate(s); \
expected at most one per mate.",
self.adapter_sequence.len(),
));
}
aggregate_errors(errors)
}
fn check_read_structures(&self, errors: &mut Vec<String>) {
if self.read_structures.is_empty() || self.discard_unsupported_segments {
return;
}
for (idx, rs) in self.read_structures.iter().enumerate() {
for seg in rs.iter() {
match seg.kind {
SegmentType::SampleBarcode => errors.push(format!(
"Read-structure for mate {} ({rs}) contains a sample barcode (B) segment. \
Sample barcodes are not supported by `chelae trim` — run `fqtk demux` first \
to assign reads to samples, or pass `--discard-unsupported-segments` to \
treat B as skip.",
idx + 1
)),
SegmentType::CellularBarcode => errors.push(format!(
"Read-structure for mate {} ({rs}) contains a cellular barcode (C) \
segment. Cellular barcodes have no widely-adopted FASTQ convention; \
single-cell tools (CellRanger, STARsolo, alevin, kallisto-bustools) \
expect CB to remain in the R1 sequence. Leave C as template, or pass \
`--discard-unsupported-segments` to treat as skip.",
idx + 1
)),
SegmentType::Template | SegmentType::MolecularBarcode | SegmentType::Skip => {}
_ => {}
}
}
}
}
fn check_adapter_args(&self, errors: &mut Vec<String>) {
if !(0.0..=1.0).contains(&self.adapter_mismatch_rate) {
errors.push(format!(
"--adapter-mismatch-rate must be in 0.0..=1.0, got {}.",
self.adapter_mismatch_rate
));
}
if !(0.0..=1.0).contains(&self.overlap_max_mismatch_rate) {
errors.push(format!(
"--overlap-max-mismatch-rate must be in 0.0..=1.0, got {}.",
self.overlap_max_mismatch_rate
));
}
if self.adapter_min_length == 0 {
errors.push("--adapter-min-length must be at least 1.".to_string());
}
if self.overlap_min_length == 0 {
errors.push("--overlap-min-length must be at least 1.".to_string());
}
if self.overlap_diagnostic_length == 0 {
errors.push("--overlap-diagnostic-length must be at least 1.".to_string());
}
if !(0.0..=1.0).contains(&self.overlap_trust_max_chance) {
errors.push(format!(
"--overlap-trust-max-chance must be in [0, 1], got {}.",
self.overlap_trust_max_chance
));
}
for seq in &self.adapter_sequence {
if seq.is_empty() {
errors.push("--adapter-sequence values must not be empty.".to_string());
continue;
}
if let Err(msg) = validate_adapter_bases(seq.as_bytes()) {
errors.push(format!("--adapter-sequence {seq:?}: {msg}"));
}
}
for kit in &self.kit {
if expand_kit_name(kit).is_none() {
errors.push(format!(
"--kit {kit:?} is not a recognized preset. Known: truseq, nextera, small-rna, \
aviti, mgi (alias dnbseq), all."
));
}
}
if let Some(fasta) = &self.adapter_fasta
&& !fasta.exists()
{
errors.push(format!("--adapter-fasta path {fasta:?} does not exist."));
}
}
fn check_filter_args(&self, errors: &mut Vec<String>) {
if let Some(n) = self.trim_polyx
&& n == 0
{
errors.push("--trim-polyx min run must be >= 1.".to_string());
}
}
fn check_no_output_overwrites_input(
&self,
inputs: &[PathBuf],
outputs: &[PathBuf],
errors: &mut Vec<String>,
) {
let input_abs: Vec<PathBuf> = inputs
.iter()
.filter(|p| p.as_os_str() != "-" && p.exists())
.filter_map(|p| resolve_real_path(p))
.collect();
let mut check = |candidate: &Path, label: &str| {
if candidate.as_os_str() == "-" {
return;
}
let abs = match resolve_real_path(candidate) {
Some(a) => a,
None => return,
};
for input in &input_abs {
if &abs == input {
errors.push(format!(
"{label} path {candidate:?} resolves to an input file ({input:?}); \
refusing to overwrite."
));
}
}
};
for out in outputs {
check(out, "Output");
}
if let Some(m) = &self.metrics {
check(m, "Metrics");
}
if let Some(j) = &self.json {
check(j, "JSON");
}
}
fn check_no_duplicate_outputs(&self, outputs: &[PathBuf], errors: &mut Vec<String>) {
let mut labeled: Vec<(String, PathBuf)> = Vec::new();
for (i, out) in outputs.iter().enumerate() {
if out.as_os_str() == "-" {
continue;
}
if let Some(abs) = resolve_real_path(out) {
labeled.push((format!("Output {}", i + 1), abs));
}
}
if let Some(m) = &self.metrics
&& m.as_os_str() != "-"
&& let Some(abs) = resolve_real_path(m)
{
labeled.push(("--metrics".to_string(), abs));
}
if let Some(j) = &self.json
&& j.as_os_str() != "-"
&& let Some(abs) = resolve_real_path(j)
{
labeled.push(("--json".to_string(), abs));
}
for i in 0..labeled.len() {
for j in (i + 1)..labeled.len() {
if labeled[i].1 == labeled[j].1 {
errors.push(format!(
"{} and {} both resolve to {:?}; each output/metrics/json path must be \
distinct.",
labeled[i].0, labeled[j].0, labeled[i].1
));
}
}
}
}
fn check_metrics_json_not_dash(&self, errors: &mut Vec<String>) {
if let Some(m) = &self.metrics
&& m.as_os_str() == "-"
{
errors.push("--metrics does not support '-' (stdin/stdout); pass a file path.".into());
}
if let Some(j) = &self.json
&& j.as_os_str() == "-"
{
errors.push("--json does not support '-' (stdin/stdout); pass a file path.".into());
}
}
fn emit_summary(&self, metrics: &TrimMetrics, num_mates: usize, elapsed: std::time::Duration) {
let pass_pct = pct(metrics.reads_out, metrics.reads_in);
let filt = metrics.reads_filtered_length
+ metrics.reads_filtered_n
+ metrics.reads_filtered_quality
+ metrics.reads_filtered_low_qual;
let q20_before = metrics.q20_before_r1 + metrics.q20_before_r2;
let q30_before = metrics.q30_before_r1 + metrics.q30_before_r2;
let q20_after = metrics.q20_after_r1 + metrics.q20_after_r2;
let q30_after = metrics.q30_after_r1 + metrics.q30_after_r2;
let bases_before_qual = metrics.bases_in;
let bases_after_qual = metrics.total_bases_after_r1 + metrics.total_bases_after_r2;
let secs = elapsed.as_secs();
info!("chelae trim complete:");
info!(
" input: {} {}",
fmt_count(metrics.reads_in),
if num_mates == 2 { "pairs" } else { "reads" }
);
info!(
" output: {} ({pass_pct:.2}%) filtered: {} (length {}, n-base {}, quality {}, low-qual {})",
fmt_count(metrics.reads_out),
fmt_count(filt),
fmt_count(metrics.reads_filtered_length),
fmt_count(metrics.reads_filtered_n),
fmt_count(metrics.reads_filtered_quality),
fmt_count(metrics.reads_filtered_low_qual),
);
info!(
" bases: {} in / {} out",
fmt_count(metrics.bases_in),
fmt_count(metrics.bases_out)
);
info!(
" trimmed: read-structure {} adapter {} poly-G {} poly-X {} quality {}",
fmt_count(metrics.bases_trimmed_read_structure),
fmt_count(metrics.bases_trimmed_adapter),
fmt_count(metrics.bases_trimmed_polyg),
fmt_count(metrics.bases_trimmed_polyx),
fmt_count(metrics.bases_trimmed_quality),
);
info!(
" Q20 rate: {:.2}% \u{2192} {:.2}%",
pct(q20_before, bases_before_qual),
pct(q20_after, bases_after_qual),
);
info!(
" Q30 rate: {:.2}% \u{2192} {:.2}%",
pct(q30_before, bases_before_qual),
pct(q30_after, bases_after_qual),
);
info!(" elapsed: {}m {}s", secs / 60, secs % 60);
}
}
impl Command for Trim {
fn execute(&self) -> Result<()> {
let start = std::time::Instant::now();
self.validate()?;
let inputs = default_dash(&self.inputs);
let outputs = default_dash(&self.outputs);
info!("Trimming {} input file(s) to {} output file(s)", inputs.len(), outputs.len());
let mut sources = open_fastq_inputs(&inputs)?;
let batch_size = self.batch_size.max(1);
let read_ahead_chunk = batch_size.min(1024);
let read_ahead_buffer = 4usize;
let (num_mates, mut iters, interleaved_rule) = if inputs.len() == 2 {
let iters: Vec<_> = sources
.into_iter()
.map(|reader| {
OwnedRecordIter::new(reader).read_ahead(read_ahead_chunk, read_ahead_buffer)
})
.collect();
(2, iters, None)
} else {
let reader = sources.pop().expect("resolved to exactly one input");
let sniffed = sniff_single_input(reader)?;
let num_mates = if sniffed.is_empty {
outputs
.len()
.max(self.read_structures.len())
.max(self.adapter_sequence.len())
.clamp(1, 2)
} else if sniffed.interleaved {
2
} else {
1
};
let interleaved_rule =
(num_mates == 2).then(|| sniffed.pairing_rule.unwrap_or(PairingRule::CasavaOrBare));
let iters = vec![sniffed.records.read_ahead(read_ahead_chunk, read_ahead_buffer)];
(num_mates, iters, interleaved_rule)
};
self.validate_post_detection(num_mates, outputs.len())?;
let adapters =
build_adapter_set(&self.adapter_sequence, &self.adapter_fasta, &self.kit, num_mates)?;
let mut split_name_check = SplitNameCheck::Pending;
let first_batch = fill_batch_from_iters(
&mut iters,
batch_size,
num_mates,
interleaved_rule,
&mut split_name_check,
0,
)?;
let polyg_min_run: Option<usize> =
if self.trim_polyg == 0 { None } else { Some(self.trim_polyg) };
let overlap_adapter_library =
build_overlap_adapter_library(&self.adapter_sequence, &self.adapter_fasta, &adapters)?;
let output_index = if outputs.len() > 1 { [0, 1] } else { [0, 0] };
let output_encodings: Vec<OutputEncoding> =
outputs.iter().map(|p| resolve_output_encoding(p, self.output_compression)).collect();
if self.compression_level != DEFAULT_COMPRESSION_LEVEL
&& output_encodings.iter().all(|e| *e == OutputEncoding::Plain)
{
warn!(
"--compression-level {} has no effect: every output is plain text (name it \
*.gz or pass --output-compression bgzf for BGZF).",
self.compression_level
);
}
let cfg = PipelineConfig {
num_mates,
num_outputs: outputs.len(),
output_index,
output_encodings,
read_structures: self.read_structures.clone(),
discard_unsupported_segments: self.discard_unsupported_segments,
adapters,
use_pe_overlap: !self.no_overlap_detection && num_mates == 2,
overlap_min_length: self.overlap_min_length,
overlap_max_mismatch_rate: self.overlap_max_mismatch_rate,
overlap_diagnostic_length: self.overlap_diagnostic_length,
overlap_trust_max_chance: (self.overlap_trust_max_chance < 1.0)
.then_some(self.overlap_trust_max_chance),
overlap_adapter_library,
expected_insert_size: self.expected_insert_size,
insert_size_stats: self.insert_size_stats,
adapter_min_length: self.adapter_min_length,
adapter_mismatch_rate: self.adapter_mismatch_rate,
polyg_min_run,
polyx_min_run: self.trim_polyx,
quality_trim_3p: self.quality_trim_3p,
quality_trim_5p: self.quality_trim_5p,
filter_length: self.filter_length,
filter_max_ns: self.filter_max_ns,
filter_mean_qual: self.filter_mean_qual,
filter_low_qual: self.filter_low_qual,
};
let compression_level = CompressionLevel::new(
u8::try_from(self.compression_level).expect("compression level validated in 1..=12"),
)?;
let n_workers = self.threads;
let num_outputs = outputs.len();
let stdout_closed = AtomicBool::new(false);
let reports_requested = self.metrics.is_some() || self.json.is_some();
let mut agg = thread::scope(|s| -> Result<WorkerAggregate> {
let (batch_tx, batch_rx) = bounded::<WorkPacket>(n_workers * 2);
let (order_txs, order_rxs): (Vec<_>, Vec<_>) = (0..num_outputs)
.map(|_| bounded::<oneshot::Receiver<Result<Vec<u8>>>>(n_workers * 4))
.unzip();
let mut worker_handles = Vec::with_capacity(n_workers);
for _ in 0..n_workers {
let rx = batch_rx.clone();
let cfg_ref = &cfg;
worker_handles.push(s.spawn(move || worker_loop(rx, cfg_ref, compression_level)));
}
drop(batch_rx);
let mut writer_handles = Vec::with_capacity(num_outputs);
for (idx, order_rx) in order_rxs.into_iter().enumerate() {
let target = if outputs[idx].as_os_str() == "-" {
OutputTarget::Stdout
} else {
OutputTarget::File(outputs[idx].clone())
};
let encoding = cfg.output_encodings[idx];
let stdout_closed = &stdout_closed;
writer_handles.push(s.spawn(move || {
writer_loop(order_rx, target, encoding, stdout_closed, reports_requested)
}));
}
let first_batch_len = (first_batch.records.len() / num_mates) as u64;
let first_submit = submit_batch(first_batch, &batch_tx, &order_txs);
let mut records_read = first_batch_len;
let mut errors: Vec<anyhow::Error> = Vec::new();
match first_submit {
Err(e) => errors.push(e),
Ok(()) => loop {
if stdout_closed.load(Ordering::Relaxed) {
break;
}
let batch = match fill_batch_from_iters(
&mut iters,
batch_size,
num_mates,
interleaved_rule,
&mut split_name_check,
records_read,
) {
Ok(b) => b,
Err(e) => {
errors.push(e);
break;
}
};
if batch.records.is_empty() {
break; }
records_read += (batch.records.len() / num_mates) as u64;
if let Err(e) = submit_batch(batch, &batch_tx, &order_txs) {
errors.push(e);
break;
}
if records_read.is_multiple_of(LOG_EVERY) {
info!(
"[chelae trim] read {} {}",
fmt_count(records_read),
if num_mates == 1 { "reads" } else { "pairs" }
);
}
},
}
drop(batch_tx);
drop(order_txs);
let mut agg = WorkerAggregate::new(num_mates);
for h in worker_handles {
match h.join() {
Ok(Ok(partial)) => agg.merge(partial),
Ok(Err(e)) => errors.push(e),
Err(payload) => {
errors.push(anyhow!("worker thread panicked: {}", panic_message(&payload)));
}
}
}
for h in writer_handles {
match h.join() {
Ok(Ok(())) => {}
Ok(Err(e)) => errors.push(e),
Err(payload) => {
errors.push(anyhow!("writer thread panicked: {}", panic_message(&payload)));
}
}
}
if let Some(e) = select_most_specific_error(errors) {
return Err(e);
}
Ok(agg)
})?;
let mut metrics = agg.metrics;
flatten_mate_stats(&agg.mate_before, &agg.mate_after, &mut metrics);
if let Some(path) = &self.metrics {
DelimFile::default().write_tsv(path, std::iter::once(&metrics))?;
}
if let Some(path) = &self.json {
let insert_size = if self.insert_size_stats {
InsertSizeStats::from_aggregate(
std::mem::take(&mut agg.insert_histogram),
agg.insert_unknown,
)
} else {
None
};
let report =
FastpJsonReport::build(&metrics, &agg.mate_before, &agg.mate_after, insert_size);
write_json_report(path, &report)?;
}
self.emit_summary(&metrics, num_mates, start.elapsed());
Ok(())
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, ValueEnum)]
enum OutputCompression {
Auto,
Bgzf,
None,
}
impl std::fmt::Display for OutputCompression {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
self.to_possible_value().expect("no skipped variants").get_name().fmt(f)
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum OutputEncoding {
Bgzf,
Plain,
}
enum OutputTarget {
File(PathBuf),
Stdout,
}
#[derive(Debug, Clone, Copy)]
pub(crate) struct QualityTrim {
pub window: usize,
pub threshold: u8,
}
impl FromStr for QualityTrim {
type Err = String;
fn from_str(s: &str) -> Result<Self, Self::Err> {
let (w, q) = s
.split_once(':')
.ok_or_else(|| format!("quality trim must be WINDOW:QUAL (e.g. 4:20), got {s:?}"))?;
let window: usize = w.parse().map_err(|e| format!("invalid window {w:?}: {e}"))?;
let threshold: u8 = q.parse().map_err(|e| format!("invalid qual {q:?}: {e}"))?;
if window == 0 {
return Err("quality trim WINDOW must be >= 1".into());
}
Ok(QualityTrim { window, threshold })
}
}
#[derive(Debug, Clone, Copy)]
struct LengthFilter {
pub min: usize,
pub max: Option<usize>,
}
impl FromStr for LengthFilter {
type Err = String;
fn from_str(s: &str) -> Result<Self, Self::Err> {
if let Some((min_s, max_s)) = s.split_once(':') {
let min: usize = min_s.parse().map_err(|e| format!("invalid min {min_s:?}: {e}"))?;
let max: usize = max_s.parse().map_err(|e| format!("invalid max {max_s:?}: {e}"))?;
if max < min {
return Err(format!("length-filter max ({max}) < min ({min})"));
}
Ok(LengthFilter { min, max: Some(max) })
} else {
let min: usize = s.parse().map_err(|e| format!("invalid min {s:?}: {e}"))?;
Ok(LengthFilter { min, max: None })
}
}
}
#[derive(Debug, Clone, Copy)]
struct LowQualFilter {
pub threshold: u8,
pub max_fraction: f64,
}
impl FromStr for LowQualFilter {
type Err = String;
fn from_str(s: &str) -> Result<Self, Self::Err> {
let (q, f) = s
.split_once(':')
.ok_or_else(|| format!("low-qual filter must be Q:F (e.g. 15:0.4), got {s:?}"))?;
let threshold: u8 = q.parse().map_err(|e| format!("invalid Q {q:?}: {e}"))?;
let max_fraction: f64 = f.parse().map_err(|e| format!("invalid F {f:?}: {e}"))?;
if !(0.0..=1.0).contains(&max_fraction) {
return Err(format!("low-qual filter F must be in 0.0..=1.0, got {max_fraction}"));
}
Ok(LowQualFilter { threshold, max_fraction })
}
}
#[derive(Debug, Clone, Default, Serialize, Deserialize)]
struct TrimMetrics {
pub reads_in: u64,
pub reads_out: u64,
pub bases_in: u64,
pub bases_out: u64,
pub bases_trimmed_read_structure: u64,
pub reads_filtered_length: u64,
pub reads_filtered_quality: u64,
pub reads_filtered_n: u64,
pub reads_filtered_low_qual: u64,
pub bases_trimmed_adapter: u64,
pub bases_trimmed_quality: u64,
pub bases_trimmed_polyg: u64,
pub bases_trimmed_polyx: u64,
pub bases_filtered: u64,
pub reads_with_adapter_trimmed: u64,
pub q20_before_r1: u64,
pub q30_before_r1: u64,
pub q20_after_r1: u64,
pub q30_after_r1: u64,
pub total_bases_before_r1: u64,
pub total_bases_after_r1: u64,
pub q20_before_r2: u64,
pub q30_before_r2: u64,
pub q20_after_r2: u64,
pub q30_after_r2: u64,
pub total_bases_before_r2: u64,
pub total_bases_after_r2: u64,
}
impl TrimMetrics {
fn merge_totals(&mut self, other: &Self) {
self.reads_in += other.reads_in;
self.reads_out += other.reads_out;
self.bases_in += other.bases_in;
self.bases_out += other.bases_out;
self.bases_trimmed_read_structure += other.bases_trimmed_read_structure;
self.reads_filtered_length += other.reads_filtered_length;
self.reads_filtered_quality += other.reads_filtered_quality;
self.reads_filtered_n += other.reads_filtered_n;
self.reads_filtered_low_qual += other.reads_filtered_low_qual;
self.bases_trimmed_adapter += other.bases_trimmed_adapter;
self.bases_trimmed_quality += other.bases_trimmed_quality;
self.bases_trimmed_polyg += other.bases_trimmed_polyg;
self.bases_trimmed_polyx += other.bases_trimmed_polyx;
self.bases_filtered += other.bases_filtered;
self.reads_with_adapter_trimmed += other.reads_with_adapter_trimmed;
}
}
struct PipelineConfig {
num_mates: usize,
num_outputs: usize,
output_index: [usize; 2],
output_encodings: Vec<OutputEncoding>,
read_structures: Vec<ReadStructure>,
discard_unsupported_segments: bool,
adapters: AdapterSet,
use_pe_overlap: bool,
overlap_min_length: usize,
overlap_max_mismatch_rate: f64,
overlap_diagnostic_length: usize,
overlap_trust_max_chance: Option<f64>,
overlap_adapter_library: OverlapAdapterLibrary,
expected_insert_size: Option<usize>,
insert_size_stats: bool,
adapter_min_length: usize,
adapter_mismatch_rate: f64,
polyg_min_run: Option<usize>,
polyx_min_run: Option<usize>,
quality_trim_3p: Option<QualityTrim>,
quality_trim_5p: Option<QualityTrim>,
filter_length: LengthFilter,
filter_max_ns: Option<usize>,
filter_mean_qual: Option<u8>,
filter_low_qual: Option<LowQualFilter>,
}
struct Pipeline<'a> {
cfg: &'a PipelineConfig,
agg: WorkerAggregate,
overlap_stats: OverlapStats,
rs_seq_scratch: Vec<u8>,
rs_qual_scratch: Vec<u8>,
overlap_scratch: OverlapScratch,
umi_parts: Vec<Vec<u8>>,
serialize_bufs: Vec<Vec<u8>>,
}
impl<'a> Pipeline<'a> {
fn new(cfg: &'a PipelineConfig) -> Self {
Self {
cfg,
agg: WorkerAggregate::new(cfg.num_mates),
overlap_stats: OverlapStats::new(cfg.expected_insert_size),
rs_seq_scratch: Vec::new(),
rs_qual_scratch: Vec::new(),
overlap_scratch: OverlapScratch::default(),
umi_parts: Vec::new(),
serialize_bufs: (0..cfg.num_outputs)
.map(|_| Vec::with_capacity(bgzf::BGZF_BLOCK_SIZE * 2))
.collect(),
}
}
fn reset_batch_bufs(&mut self) {
for buf in &mut self.serialize_bufs {
buf.clear();
}
}
fn run(&mut self, records: &mut [OwnedRecord]) -> Result<()> {
let cfg = self.cfg;
let num_mates = cfg.num_mates;
let mut stage_bases = sum_seq_bases(records);
self.agg.metrics.bases_in += stage_bases;
self.agg.metrics.reads_in += 1;
for (i, rec) in records.iter().enumerate() {
let s = observe_stats(&rec.seq, &rec.qual);
self.agg.mate_before[i].absorb(&s);
}
if let Some(min_run) = cfg.polyg_min_run {
for rec in records.iter_mut() {
self.agg.metrics.bases_trimmed_polyg += trim_polyx_tail(rec, b'G', min_run);
}
}
let _ = stage_bases;
stage_bases = sum_seq_bases(records);
let mut pre_adapter_lens = [0usize; 2];
for (i, rec) in records.iter().enumerate() {
pre_adapter_lens[i] = rec.seq.len();
}
let overlap_result = if cfg.use_pe_overlap {
let center_shift = self.overlap_stats.center_shift(records[1].seq.len());
let result = detect_pe_overlap(
&records[0].seq,
&records[1].seq,
cfg.overlap_min_length,
cfg.overlap_max_mismatch_rate,
cfg.overlap_diagnostic_length,
&cfg.overlap_adapter_library,
center_shift,
cfg.insert_size_stats,
cfg.overlap_trust_max_chance,
&mut self.overlap_scratch,
);
self.overlap_stats.observe(result, cfg.insert_size_stats);
result
} else {
WalkResult { inferred_insert: None }
};
let overlap_fired = overlap_result.inferred_insert.is_some();
if let Some(insert_len) = overlap_result.inferred_insert {
if insert_len < records[0].seq.len() {
records[0].seq.truncate(insert_len);
records[0].qual.truncate(insert_len);
}
if insert_len < records[1].seq.len() {
records[1].seq.truncate(insert_len);
records[1].qual.truncate(insert_len);
}
}
let max_k = if cfg.use_pe_overlap { Some(cfg.overlap_min_length) } else { None };
if !cfg.adapters.is_empty() && !overlap_fired {
for (i, rec) in records.iter_mut().enumerate() {
let mate_adapters = cfg.adapters.for_mate(i);
if mate_adapters.is_empty() {
continue;
}
if let Some(pos) = find_best_adapter_match(
&rec.seq,
mate_adapters,
cfg.adapter_min_length,
cfg.adapter_mismatch_rate,
max_k,
) {
rec.seq.truncate(pos);
rec.qual.truncate(pos);
}
}
}
let after_adapter = sum_seq_bases(records);
self.agg.metrics.bases_trimmed_adapter += stage_bases - after_adapter;
for (i, rec) in records.iter().enumerate() {
if rec.seq.len() < pre_adapter_lens[i] {
self.agg.metrics.reads_with_adapter_trimmed += 1;
}
}
stage_bases = after_adapter;
if !cfg.read_structures.is_empty() {
self.umi_parts.clear();
let mut rs_too_short = false;
for (i, rec) in records.iter_mut().enumerate() {
match apply_read_structure(
&cfg.read_structures[i],
rec,
cfg.discard_unsupported_segments,
&mut self.umi_parts,
&mut self.rs_seq_scratch,
&mut self.rs_qual_scratch,
)? {
ApplyRsOutcome::Applied => {}
ApplyRsOutcome::TooShort => rs_too_short = true,
}
}
if rs_too_short {
self.agg.metrics.bases_filtered += sum_seq_bases(records);
self.agg.metrics.reads_filtered_length += 1;
return Ok(());
}
if !self.umi_parts.is_empty() {
let umi_suffix = join_umi(&self.umi_parts);
for rec in records.iter_mut() {
append_umi_to_head(&mut rec.head, &umi_suffix)?;
}
}
}
let after_rs = sum_seq_bases(records);
self.agg.metrics.bases_trimmed_read_structure += stage_bases - after_rs;
if let Some(min_run) = cfg.polyx_min_run {
for rec in records.iter_mut() {
let best =
b"ACT".iter().map(|&x| find_polyx_tail_len(&rec.seq, x)).max().unwrap_or(0);
if best >= min_run && best > 0 {
let new_len = rec.seq.len() - best;
self.agg.metrics.bases_trimmed_polyx += best as u64;
rec.seq.truncate(new_len);
rec.qual.truncate(new_len);
}
}
}
if let Some(qt) = cfg.quality_trim_5p {
for rec in records.iter_mut() {
self.agg.metrics.bases_trimmed_quality +=
trim_quality_sliding_5prime(rec, qt.window, qt.threshold);
}
}
if let Some(qt) = cfg.quality_trim_3p {
for rec in records.iter_mut() {
self.agg.metrics.bases_trimmed_quality +=
trim_quality_sliding_3prime(rec, qt.window, qt.threshold);
}
}
let mut post_stats = [BaseStats::default(); 2];
for (i, rec) in records.iter().enumerate() {
post_stats[i] = observe_stats(&rec.seq, &rec.qual);
}
match evaluate_filters(
records,
&post_stats[..num_mates],
cfg.filter_length,
cfg.filter_max_ns,
cfg.filter_mean_qual,
cfg.filter_low_qual,
) {
None => {
for (i, rec) in records.iter().enumerate() {
self.agg.metrics.bases_out += post_stats[i].total;
self.agg.mate_after[i].absorb(&post_stats[i]);
rec.write(&mut self.serialize_bufs[cfg.output_index[i]])
.map_err(|e| anyhow!("failed to serialize record: {e}"))?;
}
self.agg.metrics.reads_out += 1;
}
Some(reason) => {
self.agg.metrics.bases_filtered += sum_seq_bases(records);
match reason {
FilterReject::Length => self.agg.metrics.reads_filtered_length += 1,
FilterReject::NBases => self.agg.metrics.reads_filtered_n += 1,
FilterReject::Quality => self.agg.metrics.reads_filtered_quality += 1,
FilterReject::LowQual => self.agg.metrics.reads_filtered_low_qual += 1,
}
}
}
Ok(())
}
}
struct Batch {
records: Vec<OwnedRecord>,
}
struct WorkPacket {
batch: Batch,
result_txs: Vec<oneshot::Sender<Result<Vec<u8>>>>,
}
#[derive(Debug)]
struct WorkerAggregate {
metrics: TrimMetrics,
mate_before: Vec<MateStats>,
mate_after: Vec<MateStats>,
insert_histogram: Vec<u64>,
insert_unknown: u64,
}
impl WorkerAggregate {
fn new(num_mates: usize) -> Self {
Self {
metrics: TrimMetrics::default(),
mate_before: vec![MateStats::default(); num_mates],
mate_after: vec![MateStats::default(); num_mates],
insert_histogram: Vec::new(),
insert_unknown: 0,
}
}
fn merge(&mut self, other: WorkerAggregate) {
self.metrics.merge_totals(&other.metrics);
for (dst, src) in self.mate_before.iter_mut().zip(other.mate_before.iter()) {
dst.merge(src);
}
for (dst, src) in self.mate_after.iter_mut().zip(other.mate_after.iter()) {
dst.merge(src);
}
if self.insert_histogram.len() < other.insert_histogram.len() {
self.insert_histogram.resize(other.insert_histogram.len(), 0);
}
for (dst, src) in self.insert_histogram.iter_mut().zip(other.insert_histogram.iter()) {
*dst += src;
}
self.insert_unknown += other.insert_unknown;
}
}
#[derive(Debug, Clone, Copy, Default)]
struct MateStats {
reads: u64,
bases: u64,
q20_bases: u64,
q30_bases: u64,
gc_bases: u64,
}
impl MateStats {
fn absorb(&mut self, s: &BaseStats) {
self.reads += 1;
self.bases += s.total;
self.q20_bases += s.q20;
self.q30_bases += s.q30;
self.gc_bases += s.gc;
}
fn merge(&mut self, other: &Self) {
self.reads += other.reads;
self.bases += other.bases;
self.q20_bases += other.q20_bases;
self.q30_bases += other.q30_bases;
self.gc_bases += other.gc_bases;
}
}
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
struct BaseStats {
total: u64,
q20: u64,
q30: u64,
n_bases: u64,
gc: u64,
}
#[derive(Debug, Serialize)]
struct FastpJsonReport<'a> {
summary: SummarySection<'a>,
filtering_result: FilteringResultSection,
adapter_cutting: AdapterCuttingSection,
#[serde(skip_serializing_if = "Option::is_none")]
read1_before_filtering: Option<FilteringStats>,
#[serde(skip_serializing_if = "Option::is_none")]
read1_after_filtering: Option<FilteringStats>,
#[serde(skip_serializing_if = "Option::is_none")]
read2_before_filtering: Option<FilteringStats>,
#[serde(skip_serializing_if = "Option::is_none")]
read2_after_filtering: Option<FilteringStats>,
#[serde(skip_serializing_if = "Option::is_none")]
insert_size: Option<InsertSizeStats>,
command: String,
}
#[derive(Debug, Serialize)]
struct InsertSizeStats {
peak: usize,
unknown: u64,
histogram: Vec<u64>,
}
impl InsertSizeStats {
fn from_aggregate(histogram: Vec<u64>, unknown: u64) -> Option<Self> {
if histogram.iter().all(|&c| c == 0) && unknown == 0 {
return None;
}
let peak = histogram
.iter()
.enumerate()
.reduce(|best, cur| if cur.1 > best.1 { cur } else { best })
.map(|(i, _)| i)
.unwrap_or(0);
Some(Self { peak, unknown, histogram })
}
}
impl FastpJsonReport<'_> {
fn build(
metrics: &TrimMetrics,
before: &[MateStats],
after: &[MateStats],
insert_size: Option<InsertSizeStats>,
) -> Self {
let paired = before.len() == 2;
let sequencing = if paired { "paired end" } else { "single end" };
let mates = if paired { 2 } else { 1 };
let read1_before = before.first().map(FilteringStats::from_mate);
let read1_after = after.first().map(FilteringStats::from_mate);
let read2_before = before.get(1).map(FilteringStats::from_mate);
let read2_after = after.get(1).map(FilteringStats::from_mate);
let command = std::env::args().collect::<Vec<_>>().join(" ");
Self {
summary: SummarySection {
fastp_version: env!("CARGO_PKG_VERSION"),
sequencing,
before_filtering: SummaryFilteringStats::sum(before),
after_filtering: SummaryFilteringStats::sum(after),
},
filtering_result: FilteringResultSection {
passed_filter_reads: metrics.reads_out * mates,
low_quality_reads: (metrics.reads_filtered_quality
+ metrics.reads_filtered_low_qual)
* mates,
too_many_n_reads: metrics.reads_filtered_n * mates,
too_short_reads: metrics.reads_filtered_length * mates,
too_long_reads: 0, },
adapter_cutting: AdapterCuttingSection {
adapter_trimmed_reads: metrics.reads_with_adapter_trimmed,
adapter_trimmed_bases: metrics.bases_trimmed_adapter,
},
read1_before_filtering: read1_before,
read1_after_filtering: read1_after,
read2_before_filtering: read2_before,
read2_after_filtering: read2_after,
insert_size,
command,
}
}
}
#[derive(Debug, Serialize)]
struct SummarySection<'a> {
fastp_version: &'a str,
sequencing: &'a str,
before_filtering: SummaryFilteringStats,
after_filtering: SummaryFilteringStats,
}
#[derive(Debug, Serialize)]
struct SummaryFilteringStats {
total_reads: u64,
total_bases: u64,
q20_bases: u64,
q30_bases: u64,
q20_rate: f64,
q30_rate: f64,
read1_mean_length: u64,
read2_mean_length: u64,
gc_content: f64,
}
impl SummaryFilteringStats {
fn sum(mates: &[MateStats]) -> Self {
let mut agg = MateStats::default();
for m in mates {
agg.reads += m.reads;
agg.bases += m.bases;
agg.q20_bases += m.q20_bases;
agg.q30_bases += m.q30_bases;
agg.gc_bases += m.gc_bases;
}
let mean_len =
|m: Option<&MateStats>| m.and_then(|m| m.bases.checked_div(m.reads)).unwrap_or(0);
Self {
total_reads: agg.reads,
total_bases: agg.bases,
q20_bases: agg.q20_bases,
q30_bases: agg.q30_bases,
q20_rate: ratio(agg.q20_bases, agg.bases),
q30_rate: ratio(agg.q30_bases, agg.bases),
read1_mean_length: mean_len(mates.first()),
read2_mean_length: mean_len(mates.get(1)),
gc_content: ratio(agg.gc_bases, agg.bases),
}
}
}
#[derive(Debug, Serialize)]
struct FilteringStats {
total_reads: u64,
total_bases: u64,
q20_bases: u64,
q30_bases: u64,
q20_rate: f64,
q30_rate: f64,
}
impl FilteringStats {
fn from_mate(stats: &MateStats) -> Self {
Self {
total_reads: stats.reads,
total_bases: stats.bases,
q20_bases: stats.q20_bases,
q30_bases: stats.q30_bases,
q20_rate: ratio(stats.q20_bases, stats.bases),
q30_rate: ratio(stats.q30_bases, stats.bases),
}
}
}
#[derive(Debug, Serialize)]
struct FilteringResultSection {
passed_filter_reads: u64,
low_quality_reads: u64,
#[serde(rename = "too_many_N_reads")]
too_many_n_reads: u64,
too_short_reads: u64,
too_long_reads: u64,
}
#[derive(Debug, Serialize)]
struct AdapterCuttingSection {
adapter_trimmed_reads: u64,
adapter_trimmed_bases: u64,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum FilterReject {
Length,
NBases,
Quality,
LowQual,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum ApplyRsOutcome {
Applied,
TooShort,
}
#[derive(Debug, Default, Clone)]
struct AdapterSet {
r1: Vec<Adapter>,
r2: Vec<Adapter>,
}
impl AdapterSet {
fn is_empty(&self) -> bool {
self.r1.is_empty() && self.r2.is_empty()
}
fn for_mate(&self, i: usize) -> &[Adapter] {
if i == 0 { &self.r1 } else { &self.r2 }
}
}
#[derive(Debug, Clone)]
pub(crate) struct Adapter {
pub(crate) bytes: Vec<u8>,
pub(crate) pure_acgt: bool,
}
impl Adapter {
pub(crate) fn new(bytes: Vec<u8>) -> Self {
let pure_acgt = bytes
.iter()
.all(|&b| matches!(b, b'A' | b'C' | b'G' | b'T' | b'a' | b'c' | b'g' | b't'));
Self { bytes, pure_acgt }
}
}
#[derive(Debug, Clone)]
pub(crate) struct OverlapStats {
sum_insert: u64,
count_detect: u64,
pairs_since_update: u64,
expected_insert: Option<usize>,
histogram: Vec<u64>,
unknown: u64,
}
impl OverlapStats {
pub(crate) fn new(hint: Option<usize>) -> Self {
Self {
sum_insert: 0,
count_detect: 0,
pairs_since_update: 0,
expected_insert: hint,
histogram: Vec::new(),
unknown: 0,
}
}
pub(crate) fn center_shift(&self, r2_len: usize) -> isize {
match self.expected_insert {
Some(i) => (i as isize) - (r2_len as isize),
None => isize::MIN,
}
}
pub(crate) fn observe(&mut self, result: WalkResult, stats_on: bool) {
if let Some(insert) = result.inferred_insert {
self.sum_insert += insert as u64;
self.count_detect += 1;
if stats_on {
if self.histogram.len() <= insert {
self.histogram.resize(insert + 1, 0);
}
self.histogram[insert] += 1;
}
} else if stats_on {
self.unknown += 1;
}
self.pairs_since_update += 1;
self.maybe_update_expected_insert();
}
fn maybe_update_expected_insert(&mut self) {
if self.pairs_since_update < INSERT_STATS_UPDATE_INTERVAL {
return;
}
self.pairs_since_update = 0;
if self.count_detect < INSERT_STATS_MIN_DETECTIONS {
return;
}
let mean_insert = (self.sum_insert as f64) / (self.count_detect as f64);
let new_estimate = mean_insert.round() as usize;
match self.expected_insert {
Some(old) => {
let margin = (old / 20).max(1);
if new_estimate.abs_diff(old) >= margin {
self.expected_insert = Some(new_estimate);
}
}
None => self.expected_insert = Some(new_estimate),
}
}
}
#[derive(Debug, Clone, Copy)]
pub(crate) struct WalkResult {
pub(crate) inferred_insert: Option<usize>,
}
#[derive(Debug, Clone, Default)]
pub(crate) struct OverlapAdapterLibrary {
r1_prefixes: Vec<Vec<u8>>,
r2_prefixes: Vec<Vec<u8>>,
chance_table: Vec<f64>,
}
impl OverlapAdapterLibrary {
const TAIL_LENS: usize = ADAPTER_EVIDENCE_PROBE_LEN + 1;
const MISMATCH_COUNTS: usize = 2 * ADAPTER_EVIDENCE_PROBE_LEN + 1;
fn new(r1_prefixes: Vec<Vec<u8>>, r2_prefixes: Vec<Vec<u8>>) -> Self {
let (k1, k2) = (r1_prefixes.len(), r2_prefixes.len());
let pmfs1: Vec<Vec<f64>> = (0..Self::TAIL_LENS).map(|n| best_mismatch_pmf(n, k1)).collect();
let pmfs2: Vec<Vec<f64>> = (0..Self::TAIL_LENS).map(|n| best_mismatch_pmf(n, k2)).collect();
let mut chance_table = vec![0.0; Self::TAIL_LENS * Self::TAIL_LENS * Self::MISMATCH_COUNTS];
for (n1, pmf1) in pmfs1.iter().enumerate() {
for (n2, pmf2) in pmfs2.iter().enumerate() {
let base = (n1 * Self::TAIL_LENS + n2) * Self::MISMATCH_COUNTS;
for x in 0..Self::MISMATCH_COUNTS {
chance_table[base + x] = pmf1
.iter()
.enumerate()
.flat_map(|(i, p1)| {
pmf2.iter().take((x + 1).saturating_sub(i)).map(move |p2| p1 * p2)
})
.sum::<f64>()
.min(1.0);
}
}
}
Self { r1_prefixes, r2_prefixes, chance_table }
}
fn is_empty(&self) -> bool {
self.r1_prefixes.is_empty() && self.r2_prefixes.is_empty()
}
fn chance(&self, n1: usize, n2: usize, mismatches: usize) -> f64 {
if self.chance_table.is_empty() || n1 + n2 == 0 {
return 1.0;
}
let x = mismatches.min(Self::MISMATCH_COUNTS - 1);
self.chance_table[(n1 * Self::TAIL_LENS + n2) * Self::MISMATCH_COUNTS + x]
}
}
#[derive(Debug, Clone, Copy)]
struct AcceptedOverlap {
insert: usize,
probe_mismatches: usize,
probe_len: usize,
tail_mismatches: usize,
tail_compared: usize,
chance: f64,
}
impl AcceptedOverlap {
fn trustworthy(&self, max_chance: f64) -> bool {
if self.tail_compared > 0 { self.chance <= max_chance } else { self.probe_mismatches == 0 }
}
fn better_than(&self, other: &AcceptedOverlap, max_chance: f64) -> bool {
let trusted = self.trustworthy(max_chance);
if trusted != other.trustworthy(max_chance) {
return trusted;
}
let (m1, n1) =
(self.probe_mismatches + self.tail_mismatches, self.probe_len + self.tail_compared);
let (m2, n2) =
(other.probe_mismatches + other.tail_mismatches, other.probe_len + other.tail_compared);
match (m1 * n2).cmp(&(m2 * n1)) {
std::cmp::Ordering::Less => true,
std::cmp::Ordering::Greater => false,
std::cmp::Ordering::Equal => self.insert > other.insert,
}
}
}
enum ProbeOutcome {
Accept(AcceptedOverlap),
EvidenceFail,
ProbeFail,
}
#[derive(Debug, Default)]
pub(crate) struct OverlapScratch {
r2_rc: Vec<u8>,
screen: NegShiftScreen,
}
#[derive(Debug, Default)]
struct NegShiftScreen {
survivors: Vec<u64>,
state: ScreenState,
}
impl NegShiftScreen {
fn reset(&mut self) {
self.state = ScreenState::Pending;
}
fn ensure(
&mut self,
r1: &[u8],
r2_rc: &[u8],
max_abs: usize,
max_mm_rate: f64,
diagnostic_len: usize,
) -> bool {
if self.state == ScreenState::Pending {
self.state = self.screen(r1, r2_rc, max_abs, max_mm_rate, diagnostic_len);
}
self.state == ScreenState::Ready
}
fn screen(
&mut self,
r1: &[u8],
r2_rc: &[u8],
max_abs: usize,
max_mm_rate: f64,
diagnostic_len: usize,
) -> ScreenState {
let probe_len = |abs: usize| (r2_rc.len() - abs).min(r1.len()).min(diagnostic_len);
if probe_len(max_abs) < 16 {
return ScreenState::Unavailable;
}
self.survivors.clear();
self.survivors.resize(max_abs / 64 + 1, 0);
let prefix: [u8x16; 16] = std::array::from_fn(|j| u8x16::splat(r1[j] | 0x20));
for base in (0..=max_abs).step_by(16) {
let budget = (probe_len(base) as f64 * max_mm_rate).floor() as usize;
let mut mask = match r2_rc.get(base..base + 31) {
Some(window) if budget < 16 => {
screen_16_offsets(&prefix, window.try_into().expect("31-byte window"), budget)
}
_ => u16::MAX,
};
let lanes = max_abs - base + 1;
if lanes < 16 {
mask &= (1 << lanes) - 1;
}
self.survivors[base / 64] |= u64::from(mask) << (base % 64);
}
ScreenState::Ready
}
fn next_at_or_above(&self, from: usize) -> Option<usize> {
let mut word = from / 64;
let mut bits = self.survivors.get(word)? & (u64::MAX << (from % 64));
while bits == 0 {
word += 1;
bits = *self.survivors.get(word)?;
}
Some(word * 64 + bits.trailing_zeros() as usize)
}
fn prev_at_or_below(&self, from: usize) -> Option<usize> {
let mut word = from / 64;
let mut bits = self.survivors[word] & (u64::MAX >> (63 - from % 64));
while bits == 0 {
word = word.checked_sub(1)?;
bits = self.survivors[word];
}
Some(word * 64 + 63 - bits.leading_zeros() as usize)
}
}
#[derive(Debug, Default, Clone, Copy, PartialEq, Eq)]
enum ScreenState {
#[default]
Pending,
Ready,
Unavailable,
}
fn fill_batch_from_iters<I>(
iters: &mut [I],
batch_size: usize,
num_mates: usize,
interleaved_rule: Option<PairingRule>,
split_name_check: &mut SplitNameCheck,
seen_before: u64,
) -> Result<Batch>
where
I: Iterator<Item = Result<OwnedRecord>>,
{
let mut records: Vec<OwnedRecord> = Vec::with_capacity(batch_size * num_mates);
for slot_idx in 0..batch_size {
let record_idx = seen_before + slot_idx as u64 + 1;
let pulled = match interleaved_rule {
Some(rule) => match pull_pair_interleaved(&mut iters[0], rule, record_idx)? {
Some((r1, r2)) => {
records.push(r1);
records.push(r2);
true
}
None => false,
},
None => {
pull_per_file_slot(iters, num_mates, record_idx, split_name_check, &mut records)?
}
};
if !pulled {
break;
}
}
Ok(Batch { records })
}
fn pull_per_file_slot<I>(
iters: &mut [I],
num_mates: usize,
record_idx: u64,
split_name_check: &mut SplitNameCheck,
records: &mut Vec<OwnedRecord>,
) -> Result<bool>
where
I: Iterator<Item = Result<OwnedRecord>>,
{
let slot_start = records.len();
let mut eof_count = 0usize;
for iter in iters.iter_mut() {
match iter.next() {
Some(Ok(rec)) => records.push(rec),
Some(Err(e)) => return Err(e),
None => eof_count += 1,
}
}
if eof_count == num_mates {
return Ok(false);
}
let mates = &records[slot_start..];
anyhow::ensure!(
mates.len() == num_mates,
"FASTQ files are out of sync: {}/{} files produced a record at record {record_idx}",
mates.len(),
num_mates,
);
if num_mates == 2 {
split_name_check.check(&mates[0].head, &mates[1].head, record_idx)?;
}
Ok(true)
}
fn submit_batch(
batch: Batch,
batch_tx: &Sender<WorkPacket>,
order_txs: &[Sender<oneshot::Receiver<Result<Vec<u8>>>>],
) -> Result<()> {
let mut result_txs = Vec::with_capacity(order_txs.len());
for order_tx in order_txs {
let (tx, rx) = oneshot::channel::<Result<Vec<u8>>>();
order_tx.send(rx).map_err(|_| anyhow!("writer exited before receiving batch slot"))?;
result_txs.push(tx);
}
batch_tx
.send(WorkPacket { batch, result_txs })
.map_err(|_| anyhow!("workers exited before receiving batch"))?;
Ok(())
}
fn worker_loop(
batch_rx: Receiver<WorkPacket>,
cfg: &PipelineConfig,
compression_level: CompressionLevel,
) -> Result<WorkerAggregate> {
let output_encodings = &cfg.output_encodings;
let num_outputs = cfg.num_outputs;
let mut compressors: Vec<Compressor> =
(0..num_outputs).map(|_| Compressor::new(compression_level)).collect();
let mut pipeline = Pipeline::new(cfg);
while let Ok(packet) = batch_rx.recv() {
let WorkPacket { mut batch, result_txs } = packet;
debug_assert_eq!(result_txs.len(), num_outputs);
pipeline.reset_batch_bufs();
let processed: Result<()> = (|| {
debug_assert_eq!(batch.records.len() % cfg.num_mates, 0, "partial slot in batch");
for mates in batch.records.chunks_exact_mut(cfg.num_mates) {
pipeline.run(mates)?;
}
Ok(())
})();
let mut per_output = match processed {
Ok(()) => match encode_outputs(
&mut compressors,
&mut pipeline.serialize_bufs,
output_encodings,
) {
Ok(v) => v.into_iter().map(Ok).collect::<Vec<_>>(),
Err(e) => {
let msg = format!("{e}");
(0..num_outputs).map(|_| Err(anyhow!("{msg}"))).collect()
}
},
Err(e) => {
let msg = format!("{e}");
(0..num_outputs).map(|_| Err(anyhow!("{msg}"))).collect()
}
};
for tx in result_txs.into_iter() {
let payload = per_output.remove(0);
if tx.send(payload).is_err() {
return Err(anyhow!("writer dropped before worker could deliver batch output"));
}
}
}
pipeline.agg.insert_histogram = std::mem::take(&mut pipeline.overlap_stats.histogram);
pipeline.agg.insert_unknown = pipeline.overlap_stats.unknown;
Ok(pipeline.agg)
}
fn encode_outputs(
compressors: &mut [Compressor],
serialize_bufs: &mut [Vec<u8>],
encodings: &[OutputEncoding],
) -> Result<Vec<Vec<u8>>> {
let mut out: Vec<Vec<u8>> = Vec::with_capacity(serialize_bufs.len());
let mut block_buf: Vec<u8> = Vec::with_capacity(bgzf::BGZF_BLOCK_SIZE);
for (m, buf) in serialize_bufs.iter_mut().enumerate() {
match encodings[m] {
OutputEncoding::Bgzf => {
let mut compressed = Vec::with_capacity(buf.len().max(1024));
let mut offset = 0;
while offset < buf.len() {
let end = (offset + bgzf::BGZF_BLOCK_SIZE).min(buf.len());
block_buf.clear();
compressors[m]
.compress(&buf[offset..end], &mut block_buf)
.map_err(|e| anyhow!("BGZF compression failed: {e}"))?;
compressed.extend_from_slice(&block_buf);
offset = end;
}
out.push(compressed);
}
OutputEncoding::Plain => {
let capacity = buf.capacity();
out.push(std::mem::replace(buf, Vec::with_capacity(capacity)));
}
}
}
Ok(out)
}
fn writer_loop(
order_rx: Receiver<oneshot::Receiver<Result<Vec<u8>>>>,
target: OutputTarget,
encoding: OutputEncoding,
stdout_closed: &AtomicBool,
reports_requested: bool,
) -> Result<()> {
let is_stdout = matches!(target, OutputTarget::Stdout);
let note_stdout_closed = || {
stdout_closed.store(true, Ordering::Relaxed);
if reports_requested {
warn!(
"stdout closed by downstream reader; stopping output early. Counts in \
--metrics/--json may include reads chelae processed that never made it out \
before the pipe closed."
);
} else {
info!("stdout closed by downstream reader; stopping output early");
}
};
let mut writer: BufWriter<Box<dyn Write>> = match target {
OutputTarget::File(path) => {
let file = File::create(&path).map_err(|e| anyhow!("creating output {path:?}: {e}"))?;
BufWriter::with_capacity(256 * 1024, Box::new(file))
}
OutputTarget::Stdout => {
BufWriter::with_capacity(
256 * 1024,
Box::new(std::io::stdout().lock()) as Box<dyn Write>,
)
}
};
while let Ok(result_rx) = order_rx.recv() {
let bytes = result_rx
.recv()
.map_err(|_| anyhow!("worker dropped without delivering a batch"))??;
if is_stdout && stdout_closed.load(Ordering::Relaxed) {
continue; }
if let Err(e) = writer.write_all(&bytes) {
if is_stdout && e.kind() == std::io::ErrorKind::BrokenPipe {
note_stdout_closed();
continue;
}
return Err(e.into());
}
}
if is_stdout && stdout_closed.load(Ordering::Relaxed) {
return Ok(()); }
let broken_pipe = |e: &std::io::Error| is_stdout && e.kind() == std::io::ErrorKind::BrokenPipe;
if encoding == OutputEncoding::Bgzf {
let mut eof = Vec::with_capacity(28);
Compressor::append_eof(&mut eof);
if let Err(e) = writer.write_all(&eof) {
if !broken_pipe(&e) {
return Err(e.into());
}
note_stdout_closed();
return Ok(());
}
}
if let Err(e) = writer.flush() {
if !broken_pipe(&e) {
return Err(e.into());
}
note_stdout_closed();
}
Ok(())
}
const GZIP_OUTPUT_EXTENSIONS: [&str; 2] = ["gz", "bgz"];
fn resolve_output_encoding(path: &Path, mode: OutputCompression) -> OutputEncoding {
let has_gz_ext = path
.extension()
.and_then(|e| e.to_str())
.is_some_and(|e| GZIP_OUTPUT_EXTENSIONS.iter().any(|gz| e.eq_ignore_ascii_case(gz)));
match mode {
OutputCompression::Bgzf => OutputEncoding::Bgzf,
OutputCompression::None => {
if has_gz_ext {
warn!(
"--output-compression none with {path:?}: writing uncompressed data to a \
.gz/.bgz-suffixed path."
);
}
OutputEncoding::Plain
}
OutputCompression::Auto => {
if has_gz_ext {
OutputEncoding::Bgzf
} else {
OutputEncoding::Plain
}
}
}
}
fn sum_seq_bases(records: &[OwnedRecord]) -> u64 {
records.iter().map(|r| r.seq().len() as u64).sum()
}
fn flatten_mate_stats(before: &[MateStats], after: &[MateStats], m: &mut TrimMetrics) {
if let Some(b) = before.first() {
m.total_bases_before_r1 = b.bases;
m.q20_before_r1 = b.q20_bases;
m.q30_before_r1 = b.q30_bases;
}
if let Some(a) = after.first() {
m.total_bases_after_r1 = a.bases;
m.q20_after_r1 = a.q20_bases;
m.q30_after_r1 = a.q30_bases;
}
if let Some(b) = before.get(1) {
m.total_bases_before_r2 = b.bases;
m.q20_before_r2 = b.q20_bases;
m.q30_before_r2 = b.q30_bases;
}
if let Some(a) = after.get(1) {
m.total_bases_after_r2 = a.bases;
m.q20_after_r2 = a.q20_bases;
m.q30_after_r2 = a.q30_bases;
}
}
fn observe_stats(seq: &[u8], qual: &[u8]) -> BaseStats {
const PHRED33: u8 = 33;
debug_assert_eq!(seq.len(), qual.len());
let total = qual.len() as u64;
let mut q20 = 0u64;
let mut q30 = 0u64;
let mut n_bases = 0u64;
let mut gc = 0u64;
let q20_thr = u8x16::splat(PHRED33 + 20);
let q30_thr = u8x16::splat(PHRED33 + 30);
let case_mask = u8x16::splat(0x20);
let n_lc = u8x16::splat(b'n');
let g_lc = u8x16::splat(b'g');
let c_lc = u8x16::splat(b'c');
let mut qual_chunks = qual.chunks_exact(16);
let mut seq_chunks = seq.chunks_exact(16);
loop {
let mut q20_lanes = u8x16::splat(0);
let mut q30_lanes = u8x16::splat(0);
let mut n_lanes = u8x16::splat(0);
let mut gc_lanes = u8x16::splat(0);
let mut chunks_in_fold = 0usize;
for (qchunk, schunk) in qual_chunks.by_ref().zip(seq_chunks.by_ref()).take(255) {
let qv = u8x16::new(qchunk.try_into().unwrap());
let sv = u8x16::new(schunk.try_into().unwrap()) | case_mask;
q20_lanes -= qv.simd_ge(q20_thr);
q30_lanes -= qv.simd_ge(q30_thr);
n_lanes -= sv.simd_eq(n_lc);
gc_lanes -= sv.simd_eq(g_lc) | sv.simd_eq(c_lc);
chunks_in_fold += 1;
}
q20 += sum_u8_lanes(q20_lanes);
q30 += sum_u8_lanes(q30_lanes);
n_bases += sum_u8_lanes(n_lanes);
gc += sum_u8_lanes(gc_lanes);
if chunks_in_fold < 255 {
break;
}
}
for (&q, &s) in qual_chunks.remainder().iter().zip(seq_chunks.remainder().iter()) {
let phred = q.saturating_sub(PHRED33);
if phred >= 20 {
q20 += 1;
}
if phred >= 30 {
q30 += 1;
}
let s_lc = s | 0x20;
if s_lc == b'n' {
n_bases += 1;
}
if s_lc == b'g' || s_lc == b'c' {
gc += 1;
}
}
BaseStats { total, q20, q30, n_bases, gc }
}
fn sum_u8_lanes(lanes: u8x16) -> u64 {
u64::from(lanes.to_array().iter().map(|&lane| u16::from(lane)).sum::<u16>())
}
pub(crate) fn count_mismatches_ci_bounded(a: &[u8], b: &[u8], limit: usize) -> usize {
debug_assert_eq!(a.len(), b.len());
let case_mask = u8x16::splat(0x20);
let mut count = 0usize;
let chunks = a.len() / 16;
for i in 0..chunks {
let start = i * 16;
let av = u8x16::new(a[start..start + 16].try_into().unwrap()) | case_mask;
let bv = u8x16::new(b[start..start + 16].try_into().unwrap()) | case_mask;
let matches = av.simd_eq(bv).to_bitmask().count_ones() as usize;
count += 16 - matches;
if count > limit {
return count;
}
}
for j in (chunks * 16)..a.len() {
if !a[j].eq_ignore_ascii_case(&b[j]) {
count += 1;
if count > limit {
return count;
}
}
}
count
}
fn screen_16_offsets(prefix: &[u8x16; 16], window: &[u8; 31], budget: usize) -> u16 {
let case_mask = u8x16::splat(0x20);
let mut matches = u8x16::splat(0);
for (j, base) in prefix.iter().enumerate() {
let shifted = u8x16::new(window[j..j + 16].try_into().unwrap()) | case_mask;
matches -= base.simd_eq(shifted);
}
matches.simd_gt(u8x16::splat((15 - budget) as u8)).to_bitmask() as u16
}
fn reverse_complement_acgt_into(seq: &[u8], out: &mut Vec<u8>) {
const COMP_LUT: [u8; 16] = [
b'N', b'T', b'N', b'G', b'A', b'N', b'N', b'C', b'N', b'N', b'N', b'N', b'N', b'N', b'N',
b'N',
];
const REV_IDX: [u8; 16] = [15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0];
let n = seq.len();
out.clear();
out.resize(n, 0u8);
let lut = u8x16::new(COMP_LUT);
let rev_idx = u8x16::new(REV_IDX);
let nibble_mask = u8x16::splat(0x0F);
let case_mask = u8x16::splat(0x20);
let upper_mask = u8x16::splat(0xDF);
let chunks = n / 16;
for i in 0..chunks {
let start = i * 16;
let v = u8x16::new(seq[start..start + 16].try_into().unwrap());
let case_bits = v & case_mask; let upper = v & upper_mask;
let nibs = upper & nibble_mask;
let complemented = lut.shuffle(nibs);
let cased = complemented | case_bits;
let reversed = cased.shuffle(rev_idx);
let out_start = n - start - 16;
out[out_start..out_start + 16].copy_from_slice(reversed.as_array());
}
let tail_start = chunks * 16;
for j in tail_start..n {
let b = seq[j];
let case_bit = b & 0x20;
let upper = b & 0xDF;
let comp = COMP_LUT[(upper & 0x0F) as usize] | case_bit;
out[n - 1 - j] = comp;
}
}
pub(crate) fn find_polyx_tail_len(seq: &[u8], x: u8) -> usize {
let target = u8x16::splat(x | 0x20);
let case = u8x16::splat(0x20);
let mut count = 0usize;
let mut pos = seq.len();
while pos >= 16 {
let chunk = u8x16::new(seq[pos - 16..pos].try_into().unwrap()) | case;
let bits = chunk.simd_eq(target).to_bitmask() as u16;
let trailing = bits.leading_ones() as usize;
count += trailing;
if trailing < 16 {
return count;
}
pos -= 16;
}
for &b in seq[..pos].iter().rev() {
if (b | 0x20) == (x | 0x20) {
count += 1;
} else {
break;
}
}
count
}
fn trim_polyx_tail(rec: &mut OwnedRecord, x: u8, min_run: usize) -> u64 {
let tail = find_polyx_tail_len(&rec.seq, x);
if tail >= min_run && tail > 0 {
let new_len = rec.seq.len() - tail;
rec.seq.truncate(new_len);
rec.qual.truncate(new_len);
tail as u64
} else {
0
}
}
fn trim_quality_sliding_3prime(rec: &mut OwnedRecord, window: usize, threshold: u8) -> u64 {
const PHRED33: u8 = 33;
let qual = &rec.qual;
if qual.len() < window || window == 0 {
return 0;
}
let offset_total = (threshold as u16) + (PHRED33 as u16);
let max_s = qual.len() - window; let mut trim_pos = qual.len();
let can_simd = window <= 16 && offset_total <= 127 && qual.len() >= window + 16;
let mut scan_start: usize = max_s;
let mut done = false;
if can_simd {
let offset_vec = i8x16::splat(offset_total as i8);
let zero_vec = i8x16::splat(0);
let mut s_end = max_s; while s_end >= 15 {
let s_start = s_end - 15;
let mut sums = i8x16::splat(0);
for i in 0..window {
let bytes: [u8; 16] = qual[s_start + i..s_start + i + 16].try_into().unwrap();
let v = i8x16::new(bytemuck::cast(bytes)).saturating_sub(offset_vec);
sums = sums.saturating_add(v);
}
let bad_bits = sums.simd_lt(zero_vec).to_bitmask() as u16;
if bad_bits == 0xFFFF {
trim_pos = s_start;
if s_start == 0 {
done = true;
break;
}
s_end = s_start - 1;
} else {
let good_bits = !bad_bits;
let p = 15 - good_bits.leading_zeros() as usize;
if p < 15 {
trim_pos = s_start + p + 1;
}
done = true;
break;
}
}
scan_start = if done { 0 } else { s_end };
}
if !done {
if scan_start + window <= qual.len() {
let win = window as u32;
let threshold_total = u32::from(threshold) * win;
let mut s = scan_start;
let mut sum: u32 =
qual[s..s + window].iter().map(|&q| u32::from(q.saturating_sub(PHRED33))).sum();
if sum < threshold_total {
trim_pos = s;
while s > 0 {
sum += u32::from(qual[s - 1].saturating_sub(PHRED33));
sum -= u32::from(qual[s - 1 + window].saturating_sub(PHRED33));
s -= 1;
if sum < threshold_total {
trim_pos = s;
} else {
break;
}
}
}
}
}
let trimmed = (rec.seq.len() - trim_pos) as u64;
rec.seq.truncate(trim_pos);
rec.qual.truncate(trim_pos);
trimmed
}
fn trim_quality_sliding_5prime(rec: &mut OwnedRecord, window: usize, threshold: u8) -> u64 {
let Some(cut_at) = cut_right_quality_position(&rec.qual, window, threshold) else {
return 0;
};
let removed = (rec.seq.len() - cut_at) as u64;
rec.seq.truncate(cut_at);
rec.qual.truncate(cut_at);
removed
}
pub(crate) fn cut_right_quality_position(
qual: &[u8],
window: usize,
threshold: u8,
) -> Option<usize> {
const PHRED33: u8 = 33;
if qual.len() < window || window == 0 {
return None;
}
let win = window as u32;
let threshold_total = u32::from(threshold) * win;
let max_s = qual.len() - window;
let mut sum: u32 = qual[..window].iter().map(|&q| u32::from(q.saturating_sub(PHRED33))).sum();
if sum < threshold_total {
return Some(0);
}
for s in 1..=max_s {
sum -= u32::from(qual[s - 1].saturating_sub(PHRED33));
sum += u32::from(qual[s + window - 1].saturating_sub(PHRED33));
if sum < threshold_total {
return Some(s);
}
}
None
}
fn count_bases_below_q(qual: &[u8], threshold: u8) -> u64 {
const PHRED33: u8 = 33;
let cutoff = PHRED33.saturating_add(threshold);
let cutoff_vec = u8x32::splat(cutoff);
let mut count = 0u64;
let chunks = qual.chunks_exact(32);
let tail = chunks.remainder();
for chunk in chunks {
let v = u8x32::new(chunk.try_into().unwrap());
count += v.simd_lt(cutoff_vec).to_bitmask().count_ones() as u64;
}
for &q in tail {
if q < cutoff {
count += 1;
}
}
count
}
fn evaluate_filters(
records: &[OwnedRecord],
stats: &[BaseStats],
length: LengthFilter,
filter_max_ns: Option<usize>,
filter_mean_qual: Option<u8>,
filter_low_qual: Option<LowQualFilter>,
) -> Option<FilterReject> {
debug_assert_eq!(records.len(), stats.len());
for rec in records {
let len = rec.seq.len();
if len == 0 || len < length.min {
return Some(FilterReject::Length);
}
if let Some(max) = length.max
&& len > max
{
return Some(FilterReject::Length);
}
}
if let Some(n_max) = filter_max_ns {
for s in stats {
if s.n_bases > n_max as u64 {
return Some(FilterReject::NBases);
}
}
}
if let Some(min_q) = filter_mean_qual {
const PHRED33: u32 = 33;
for rec in records {
let raw: u32 = rec.qual.iter().map(|&q| u32::from(q)).sum();
let threshold = (u64::from(min_q) + u64::from(PHRED33)) * rec.qual.len() as u64;
if u64::from(raw) < threshold {
return Some(FilterReject::Quality);
}
}
}
if let Some(lq) = filter_low_qual {
for rec in records {
let below = count_bases_below_q(&rec.qual, lq.threshold);
if (below as f64) > lq.max_fraction * rec.qual.len() as f64 {
return Some(FilterReject::LowQual);
}
}
}
None
}
pub(crate) fn validate_adapter_bases(seq: &[u8]) -> Result<(), String> {
for (i, &b) in seq.iter().enumerate() {
if IUPAC_MASKS[b.to_ascii_uppercase() as usize] == 0 {
return Err(format!("invalid base {:?} at position {}", b as char, i));
}
}
Ok(())
}
#[inline]
fn base_matches_iupac(read_base: u8, adapter_base: u8) -> bool {
let r = IUPAC_MASKS[read_base.to_ascii_uppercase() as usize];
let a = IUPAC_MASKS[adapter_base.to_ascii_uppercase() as usize];
r != 0 && a != 0 && (r & a) != 0
}
pub(crate) fn find_adapter_3prime(
read: &[u8],
adapter: &Adapter,
min_length: usize,
max_mm_rate: f64,
max_k: Option<usize>,
) -> Option<usize> {
if read.len() < min_length || adapter.bytes.is_empty() {
return None;
}
let max_start = read.len() - min_length;
let max_start = max_k.map(|m| m.min(max_start)).unwrap_or(max_start);
let budget = |alignment_len: usize| (alignment_len as f64 * max_mm_rate).floor() as usize;
let matches_at = |k: usize| {
let alignment_len = adapter.bytes.len().min(read.len() - k);
if alignment_len < min_length {
return false;
}
let max_mm = budget(alignment_len);
let mismatches = if adapter.pure_acgt {
count_mismatches_ci_bounded(
&read[k..k + alignment_len],
&adapter.bytes[..alignment_len],
max_mm,
)
} else {
read[k..k + alignment_len]
.iter()
.zip(adapter.bytes[..alignment_len].iter())
.filter(|(r, a)| !base_matches_iupac(**r, **a))
.count()
};
mismatches <= max_mm
};
let mut first_unscreened = 0;
if adapter.pure_acgt && adapter.bytes.len() >= 16 && read.len() >= 16 {
let screen_end = max_start.min(read.len() - 16);
let prefix: [u8x16; 16] = std::array::from_fn(|j| u8x16::splat(adapter.bytes[j] | 0x20));
for base in (0..=screen_end).step_by(16) {
let window_budget = budget(adapter.bytes.len().min(read.len() - base));
let mut survivors = if window_budget < 16 {
let mut window = [0u8; 31];
let available = (read.len() - base).min(31);
window[..available].copy_from_slice(&read[base..base + available]);
screen_16_offsets(&prefix, &window, window_budget)
} else {
u16::MAX
};
let lanes = screen_end - base + 1;
if lanes < 16 {
survivors &= (1 << lanes) - 1;
}
while survivors != 0 {
let k = base + survivors.trailing_zeros() as usize;
if matches_at(k) {
return Some(k);
}
survivors &= survivors - 1;
}
}
first_unscreened = screen_end + 1;
}
(first_unscreened..=max_start).find(|&k| matches_at(k))
}
fn find_best_adapter_match(
read: &[u8],
adapters: &[Adapter],
min_length: usize,
max_mm_rate: f64,
max_k: Option<usize>,
) -> Option<usize> {
let mut best: Option<usize> = None;
for adapter in adapters {
if let Some(k) = find_adapter_3prime(read, adapter, min_length, max_mm_rate, max_k)
&& best.is_none_or(|b| k < b)
{
best = Some(k);
}
}
best
}
#[allow(clippy::too_many_arguments)]
pub(crate) fn detect_pe_overlap(
r1: &[u8],
r2: &[u8],
min_overlap: usize,
max_mm_rate: f64,
diagnostic_len: usize,
adapter_library: &OverlapAdapterLibrary,
center: isize,
stats_on: bool,
trust_max_chance: Option<f64>,
scratch: &mut OverlapScratch,
) -> WalkResult {
if r1.len() < min_overlap || r2.len() < min_overlap {
return WalkResult { inferred_insert: None };
}
reverse_complement_acgt_into(r2, &mut scratch.r2_rc);
scratch.screen.reset();
walk_overlap(
r1,
r2,
&scratch.r2_rc,
min_overlap,
max_mm_rate,
diagnostic_len,
adapter_library,
center,
stats_on,
trust_max_chance,
&mut scratch.screen,
)
}
#[inline]
fn try_shift_neg(
r1: &[u8],
r2: &[u8],
r2_rc: &[u8],
shift: isize,
max_mm_rate: f64,
diagnostic_len: usize,
adapter_library: &OverlapAdapterLibrary,
) -> ProbeOutcome {
let r1_len = r1.len();
let r2_len = r2_rc.len();
let abs_shift = shift.unsigned_abs();
let overlap_len = (r2_len - abs_shift).min(r1_len);
let probe_len = overlap_len.min(diagnostic_len);
let max_mm = (probe_len as f64 * max_mm_rate).floor() as usize;
let mismatches = count_mismatches_ci_bounded(
&r1[..probe_len],
&r2_rc[abs_shift..abs_shift + probe_len],
max_mm,
);
if mismatches > max_mm {
return ProbeOutcome::ProbeFail;
}
let insert = r2_len - abs_shift;
let mut accepted = AcceptedOverlap {
insert,
probe_mismatches: mismatches,
probe_len,
tail_mismatches: 0,
tail_compared: 0,
chance: 1.0,
};
if shift < 0 && !adapter_library.is_empty() {
let r1_post = if insert < r1_len { &r1[insert..] } else { &[] as &[u8] };
let r2_post = if insert < r2.len() { &r2[insert..] } else { &[] as &[u8] };
let r1_best = post_cut_best_match(r1_post, &adapter_library.r1_prefixes);
let r2_best = post_cut_best_match(r2_post, &adapter_library.r2_prefixes);
let ((mm1, n1), (mm2, n2)) = match (r1_best, r2_best) {
(Some(b1), Some(b2)) => {
if b1.0 + b2.0 > combined_evidence_budget(b1.1 + b2.1) {
return ProbeOutcome::EvidenceFail;
}
(b1, b2)
}
(b1, b2) => (b1.unwrap_or((0, 0)), b2.unwrap_or((0, 0))),
};
accepted.tail_mismatches = mm1 + mm2;
accepted.tail_compared = n1 + n2;
accepted.chance = adapter_library.chance(n1, n2, mm1 + mm2);
}
ProbeOutcome::Accept(accepted)
}
#[inline]
fn try_shift_pos(
r1: &[u8],
r2_rc: &[u8],
shift: isize,
max_mm_rate: f64,
diagnostic_len: usize,
) -> ProbeOutcome {
let r1_len = r1.len();
let r2_len = r2_rc.len();
let abs_shift = shift as usize;
let overlap_len = (r1_len - abs_shift).min(r2_len);
let probe_len = overlap_len.min(diagnostic_len);
let max_mm = (probe_len as f64 * max_mm_rate).floor() as usize;
let mismatches = count_mismatches_ci_bounded(
&r1[abs_shift..abs_shift + probe_len],
&r2_rc[..probe_len],
max_mm,
);
if mismatches > max_mm {
return ProbeOutcome::ProbeFail;
}
ProbeOutcome::Accept(AcceptedOverlap {
insert: r2_len + abs_shift,
probe_mismatches: mismatches,
probe_len,
tail_mismatches: 0,
tail_compared: 0,
chance: 1.0,
})
}
#[allow(clippy::too_many_arguments)]
fn walk_overlap(
r1: &[u8],
r2: &[u8],
r2_rc: &[u8],
min_overlap: usize,
max_mm_rate: f64,
diagnostic_len: usize,
adapter_library: &OverlapAdapterLibrary,
center: isize,
stats_on: bool,
trust_max_chance: Option<f64>,
screen: &mut NegShiftScreen,
) -> WalkResult {
let r2_len = r2_rc.len();
let lo = -((r2_len - min_overlap) as isize);
let hi = if stats_on { (r1.len() - min_overlap) as isize } else { 0 };
let first = if stats_on {
walk_overlap_full(
r1,
r2,
r2_rc,
lo,
hi,
center,
max_mm_rate,
diagnostic_len,
adapter_library,
)
} else {
walk_overlap_neg(
r1,
r2,
r2_rc,
lo,
center,
max_mm_rate,
diagnostic_len,
adapter_library,
screen,
)
};
let inferred_insert = first.map(|first| match trust_max_chance {
Some(max_chance) if !first.trustworthy(max_chance) => {
best_overlap(
r1,
r2,
r2_rc,
lo,
hi,
max_mm_rate,
diagnostic_len,
adapter_library,
max_chance,
screen,
)
.unwrap_or(first)
.insert
}
_ => first.insert,
});
WalkResult { inferred_insert }
}
#[allow(clippy::too_many_arguments)]
fn best_overlap(
r1: &[u8],
r2: &[u8],
r2_rc: &[u8],
lo: isize,
hi: isize,
max_mm_rate: f64,
diagnostic_len: usize,
adapter_library: &OverlapAdapterLibrary,
max_chance: f64,
screen: &mut NegShiftScreen,
) -> Option<AcceptedOverlap> {
let mut best: Option<AcceptedOverlap> = None;
let mut consider = |outcome: ProbeOutcome| {
if let ProbeOutcome::Accept(candidate) = outcome
&& best.is_none_or(|b| candidate.better_than(&b, max_chance))
{
best = Some(candidate);
}
};
let neg =
|shift| try_shift_neg(r1, r2, r2_rc, shift, max_mm_rate, diagnostic_len, adapter_library);
let max_abs = lo.unsigned_abs();
if screen.ensure(r1, r2_rc, max_abs, max_mm_rate, diagnostic_len) {
let mut abs = screen.prev_at_or_below(max_abs);
while let Some(a) = abs {
consider(neg(-(a as isize)));
abs = a.checked_sub(1).and_then(|below| screen.prev_at_or_below(below));
}
} else {
for shift in lo..=0 {
consider(neg(shift));
}
}
for shift in 1..=hi {
consider(try_shift_pos(r1, r2_rc, shift, max_mm_rate, diagnostic_len));
}
best
}
#[inline]
#[allow(clippy::too_many_arguments)]
fn walk_overlap_neg(
r1: &[u8],
r2: &[u8],
r2_rc: &[u8],
lo: isize,
center: isize,
max_mm_rate: f64,
diagnostic_len: usize,
adapter_library: &OverlapAdapterLibrary,
screen: &mut NegShiftScreen,
) -> Option<AcceptedOverlap> {
if lo > 0 {
return None;
}
let c = center.clamp(lo, 0);
macro_rules! visit_neg {
($shift:expr) => {{
let s: isize = $shift;
if let ProbeOutcome::Accept(accepted) =
try_shift_neg(r1, r2, r2_rc, s, max_mm_rate, diagnostic_len, adapter_library)
{
return Some(accepted);
}
}};
}
visit_neg!(c);
let mut k: isize = 1;
loop {
if k > NEAR_WALK_SHIFTS
&& screen.ensure(r1, r2_rc, lo.unsigned_abs(), max_mm_rate, diagnostic_len)
{
break;
}
let mut any = false;
let down = c - k;
if down >= lo {
visit_neg!(down);
any = true;
}
let up = c + k;
if up <= 0 {
visit_neg!(up);
any = true;
}
if !any {
return None;
}
k += 1;
}
let abs_c = c.unsigned_abs();
let k = k as usize;
let mut down = screen.next_at_or_above(abs_c + k);
let mut up = abs_c.checked_sub(k).and_then(|a| screen.prev_at_or_below(a));
loop {
let abs = match (down, up) {
(Some(d), Some(u)) if d - abs_c <= abs_c - u => {
down = screen.next_at_or_above(d + 1);
d
}
(Some(d), None) => {
down = screen.next_at_or_above(d + 1);
d
}
(_, Some(u)) => {
up = u.checked_sub(1).and_then(|a| screen.prev_at_or_below(a));
u
}
(None, None) => return None,
};
visit_neg!(-(abs as isize));
}
}
#[inline]
#[allow(clippy::too_many_arguments)]
fn walk_overlap_full(
r1: &[u8],
r2: &[u8],
r2_rc: &[u8],
lo: isize,
hi: isize,
center: isize,
max_mm_rate: f64,
diagnostic_len: usize,
adapter_library: &OverlapAdapterLibrary,
) -> Option<AcceptedOverlap> {
if lo > hi {
return None;
}
let c = center.clamp(lo, hi);
macro_rules! visit {
($shift:expr) => {{
let s: isize = $shift;
let outcome = if s <= 0 {
try_shift_neg(r1, r2, r2_rc, s, max_mm_rate, diagnostic_len, adapter_library)
} else {
try_shift_pos(r1, r2_rc, s, max_mm_rate, diagnostic_len)
};
if let ProbeOutcome::Accept(accepted) = outcome {
return Some(accepted);
}
}};
}
visit!(c);
let mut k: isize = 1;
loop {
let mut any = false;
let down = c - k;
if down >= lo {
visit!(down);
any = true;
}
let up = c + k;
if up <= hi {
visit!(up);
any = true;
}
if !any {
break;
}
k += 1;
}
None
}
fn post_cut_best_match(post_cut: &[u8], library: &[Vec<u8>]) -> Option<(usize, usize)> {
if post_cut.is_empty() {
return None;
}
let mut best: Option<(usize, usize)> = None;
for prefix in library {
let n = post_cut.len().min(prefix.len()).min(ADAPTER_EVIDENCE_PROBE_LEN);
if n == 0 {
continue;
}
let limit = best.map(|(m, _)| m).unwrap_or(n);
let mm = count_mismatches_ci_bounded(&post_cut[..n], &prefix[..n], limit);
match best {
Some((m, _)) if mm >= m => {}
_ => best = Some((mm, n)),
}
}
best
}
fn combined_evidence_budget(n_total: usize) -> usize {
(n_total * ADAPTER_EVIDENCE_MAX_MM + ADAPTER_EVIDENCE_PROBE_LEN)
/ (2 * ADAPTER_EVIDENCE_PROBE_LEN)
}
fn best_mismatch_pmf(n: usize, k: usize) -> Vec<f64> {
if n == 0 || k == 0 {
return vec![1.0];
}
let mut cdf = Vec::with_capacity(n + 1);
let mut acc = 0.0;
let mut choose = 1.0; for j in 0..=n {
if j > 0 {
choose = choose * (n - j + 1) as f64 / j as f64;
}
acc += choose * 0.75f64.powi(j as i32) * 0.25f64.powi((n - j) as i32);
cdf.push(acc.min(1.0));
}
let at_least = |j: usize| if j == 0 { 1.0 } else { (1.0 - cdf[j - 1]).powi(k as i32) };
(0..=n).map(|j| at_least(j) - if j < n { at_least(j + 1) } else { 0.0 }).collect()
}
fn build_adapter_set(
adapter_sequence: &[String],
adapter_fasta: &Option<PathBuf>,
kits: &[String],
num_mates: usize,
) -> Result<AdapterSet> {
let mut r1: Vec<Vec<u8>> = Vec::new();
let mut r2: Vec<Vec<u8>> = Vec::new();
if let Some(s) = adapter_sequence.first() {
r1.push(s.as_bytes().to_ascii_uppercase());
}
if let Some(s) = adapter_sequence.get(1) {
r2.push(s.as_bytes().to_ascii_uppercase());
}
for kit_name in kits {
let kits = expand_kit_name(kit_name).ok_or_else(|| anyhow!("Unknown kit {kit_name:?}"))?;
for kit in kits {
push_kit_adapters(kit, &mut r1, &mut r2);
}
}
if let Some(path) = adapter_fasta {
let seqs = load_adapter_fasta(path)?;
for seq in seqs {
let upper = seq.to_ascii_uppercase();
r1.push(upper.clone());
if num_mates >= 2 {
r2.push(upper);
}
}
}
dedupe_sort_by_len(&mut r1);
dedupe_sort_by_len(&mut r2);
Ok(AdapterSet {
r1: r1.into_iter().map(Adapter::new).collect(),
r2: r2.into_iter().map(Adapter::new).collect(),
})
}
fn load_adapter_fasta(path: &Path) -> Result<Vec<Vec<u8>>> {
Ok(load_adapter_fasta_with_names(path)?.into_iter().map(|(_, s)| s).collect())
}
pub(crate) fn load_adapter_fasta_with_names(path: &Path) -> Result<Vec<(String, Vec<u8>)>> {
let reader = Io::new(5, BUFFER_SIZE)
.new_reader(path)
.map_err(|e| anyhow!("Failed to open adapter FASTA {path:?}: {e}"))?;
let mut out: Vec<(String, Vec<u8>)> = Vec::new();
let mut current_name: Option<String> = None;
let mut current_seq: Vec<u8> = Vec::new();
let mut in_record = false;
let mut synthetic_counter: usize = 0;
let mut push = |name: Option<String>, seq: Vec<u8>, counter: &mut usize| {
let n = name.unwrap_or_else(|| {
*counter += 1;
format!("record_{counter}")
});
out.push((n, seq));
};
for line in reader.lines() {
let line = line.map_err(|e| anyhow!("Read error in {path:?}: {e}"))?;
let trimmed = line.trim_end_matches(&['\r', '\n'][..]);
if let Some(header) = trimmed.strip_prefix('>') {
if in_record {
push(current_name.take(), std::mem::take(&mut current_seq), &mut synthetic_counter);
}
let name = header.split_whitespace().next().unwrap_or("").to_string();
current_name = if name.is_empty() { None } else { Some(name) };
in_record = true;
} else {
for &b in trimmed.as_bytes() {
if !b.is_ascii_whitespace() {
current_seq.push(b);
in_record = true;
}
}
}
}
if in_record {
push(current_name, current_seq, &mut synthetic_counter);
}
for (i, (name, seq)) in out.iter().enumerate() {
if seq.is_empty() {
return Err(anyhow!(
"--adapter-fasta record {} ({name:?}): empty sequence — header without a body",
i + 1,
));
}
validate_adapter_bases(seq)
.map_err(|m| anyhow!("--adapter-fasta record {}: {m} (sequence {seq:?})", i + 1))?;
}
Ok(out)
}
fn build_overlap_adapter_library(
adapter_sequence: &[String],
adapter_fasta: &Option<PathBuf>,
user_adapters: &AdapterSet,
) -> Result<OverlapAdapterLibrary> {
fn push(dst: &mut Vec<Vec<u8>>, seq: &[u8]) {
let n = seq.len().min(ADAPTER_EVIDENCE_PROBE_LEN);
if n >= 8 {
dst.push(seq[..n].to_ascii_uppercase());
}
}
let mut r1_prefixes: Vec<Vec<u8>> = Vec::new();
let mut r2_prefixes: Vec<Vec<u8>> = Vec::new();
for kit in chelae_lib::adapter_db::ALL_KITS {
push(&mut r1_prefixes, kit.seq_r1);
if let Some(s2) = kit.seq_r2 {
push(&mut r2_prefixes, s2);
}
}
if let Some(s) = adapter_sequence.first() {
push(&mut r1_prefixes, s.as_bytes());
}
if let Some(s) = adapter_sequence.get(1) {
push(&mut r2_prefixes, s.as_bytes());
}
if let Some(path) = adapter_fasta {
let seqs = load_adapter_fasta(path)?;
for seq in &seqs {
push(&mut r1_prefixes, seq);
push(&mut r2_prefixes, seq);
}
}
for ad in &user_adapters.r1 {
push(&mut r1_prefixes, &ad.bytes);
}
for ad in &user_adapters.r2 {
push(&mut r2_prefixes, &ad.bytes);
}
r1_prefixes.sort();
r1_prefixes.dedup();
r2_prefixes.sort();
r2_prefixes.dedup();
Ok(OverlapAdapterLibrary::new(r1_prefixes, r2_prefixes))
}
fn push_kit_adapters(kit: &KitAdapter, r1: &mut Vec<Vec<u8>>, r2: &mut Vec<Vec<u8>>) {
r1.push(kit.seq_r1.to_vec());
if let Some(s2) = kit.seq_r2 {
r2.push(s2.to_vec());
}
}
fn dedupe_sort_by_len(v: &mut Vec<Vec<u8>>) {
v.sort();
v.dedup();
v.sort_by_key(|x| std::cmp::Reverse(x.len()));
}
fn apply_read_structure(
rs: &ReadStructure,
rec: &mut OwnedRecord,
discard_unsupported: bool,
umi_parts: &mut Vec<Vec<u8>>,
template_seq: &mut Vec<u8>,
template_qual: &mut Vec<u8>,
) -> Result<ApplyRsOutcome> {
template_seq.clear();
template_qual.clear();
let extracted = match rs.extract(&rec.seq, &rec.qual, SkipHandling::Exclude) {
Ok(iter) => iter,
Err(ReadStructureError::ReadTooShort { .. }) => return Ok(ApplyRsOutcome::TooShort),
Err(e) => {
return Err(anyhow!(
"Error applying read-structure {rs} to read {}: {e}",
String::from_utf8_lossy(&rec.head)
));
}
};
for (seg, seg_seq, seg_qual) in extracted {
match seg.kind {
SegmentType::Template => {
template_seq.extend_from_slice(seg_seq);
template_qual.extend_from_slice(seg_qual);
}
SegmentType::MolecularBarcode => umi_parts.push(seg_seq.to_vec()),
SegmentType::SampleBarcode | SegmentType::CellularBarcode => {
debug_assert!(discard_unsupported);
}
SegmentType::Skip => {
debug_assert!(false, "Skip segments should be elided by include_skips=false");
}
_ => {
log::warn!(
"Unknown read-structure segment kind {:?} treated as skip — \
update chelae to handle this kind explicitly.",
seg.kind
);
}
}
}
std::mem::swap(&mut rec.seq, template_seq);
std::mem::swap(&mut rec.qual, template_qual);
Ok(ApplyRsOutcome::Applied)
}
fn join_umi(parts: &[Vec<u8>]) -> Vec<u8> {
let total = parts.iter().map(|p| p.len()).sum::<usize>() + parts.len().saturating_sub(1);
let mut out = Vec::with_capacity(total);
for (i, part) in parts.iter().enumerate() {
if i > 0 {
out.push(UMI_JOIN);
}
out.extend_from_slice(part);
}
out
}
fn append_umi_to_head(head: &mut Vec<u8>, umi: &[u8]) -> Result<()> {
let space_idx = head.iter().position(|&b| b == b' ');
let name_end = space_idx.unwrap_or(head.len());
let name = &head[..name_end];
let colons = name.iter().filter(|&&b| b == UMI_ID_SEP).count();
if colons + 1 > MAX_READ_ID_FIELDS {
return Err(anyhow!(
"Cannot append UMI to read-id with more than {MAX_READ_ID_FIELDS} colon-delimited \
fields: {}",
String::from_utf8_lossy(name)
));
}
let joiner = if colons + 1 == MAX_READ_ID_FIELDS { UMI_JOIN } else { UMI_ID_SEP };
head.splice(name_end..name_end, std::iter::once(joiner).chain(umi.iter().copied()));
Ok(())
}
fn pct(num: u64, denom: u64) -> f64 {
ratio(num, denom) * 100.0
}
fn ratio(num: u64, denom: u64) -> f64 {
if denom == 0 { 0.0 } else { num as f64 / denom as f64 }
}
fn write_json_report(path: &Path, report: &FastpJsonReport) -> Result<()> {
let file =
File::create(path).map_err(|e| anyhow!("Failed to create JSON report at {path:?}: {e}"))?;
let writer = BufWriter::new(file);
serde_json::to_writer_pretty(writer, report)
.map_err(|e| anyhow!("Failed to serialize JSON report: {e}"))?;
Ok(())
}
fn panic_message(payload: &Box<dyn std::any::Any + Send>) -> String {
if let Some(s) = payload.downcast_ref::<&'static str>() {
(*s).to_string()
} else if let Some(s) = payload.downcast_ref::<String>() {
s.clone()
} else {
"(non-string panic payload)".to_string()
}
}
fn select_most_specific_error(mut errors: Vec<anyhow::Error>) -> Option<anyhow::Error> {
let is_symptom = |e: &anyhow::Error| {
let m = e.to_string();
m.contains("exited before") || m.contains("dropped")
};
if let Some(i) = errors.iter().position(|e| !is_symptom(e)) {
Some(errors.swap_remove(i))
} else if errors.is_empty() {
None
} else {
Some(errors.swap_remove(0))
}
}
#[cfg(test)]
fn reverse_complement_into(seq: &[u8], out: &mut Vec<u8>) {
out.clear();
out.reserve(seq.len());
for &b in seq.iter().rev() {
out.push(match b {
b'A' => b'T',
b'T' | b'U' => b'A',
b'C' => b'G',
b'G' => b'C',
b'a' => b't',
b't' | b'u' => b'a',
b'c' => b'g',
b'g' => b'c',
b'R' => b'Y',
b'Y' => b'R',
b'r' => b'y',
b'y' => b'r',
b'S' => b'S',
b'W' => b'W',
b's' => b's',
b'w' => b'w',
b'K' => b'M',
b'M' => b'K',
b'k' => b'm',
b'm' => b'k',
b'B' => b'V',
b'V' => b'B',
b'b' => b'v',
b'v' => b'b',
b'D' => b'H',
b'H' => b'D',
b'd' => b'h',
b'h' => b'd',
b'N' | b'n' => b,
other => other,
});
}
}
#[cfg(test)]
fn reverse_complement(seq: &[u8]) -> Vec<u8> {
let mut out = Vec::new();
reverse_complement_into(seq, &mut out);
out
}
#[cfg(test)]
mod tests {
use super::*;
use seq_io::fastq::OwnedRecord;
use seq_io::fastq::Reader as FastqReader;
use tempfile::TempDir;
fn fq_lines(prefix: &str, reads: &[&str]) -> Vec<String> {
reads
.iter()
.enumerate()
.flat_map(|(i, &seq)| {
vec![
format!("@{prefix}_{i}"),
seq.to_string(),
"+".to_string(),
"I".repeat(seq.len()),
]
})
.collect()
}
fn write_fastq(tmp: &TempDir, name: &str, lines: &[String]) -> PathBuf {
let path = tmp.path().join(format!("{name}.fq"));
Io::default().write_lines(&path, lines).unwrap();
path
}
fn read_fastq(path: &Path) -> Vec<OwnedRecord> {
let io = Io::default();
FastqReader::new(io.new_reader(path).unwrap())
.into_records()
.collect::<Result<Vec<_>, seq_io::fastq::Error>>()
.unwrap()
}
fn fq_record(name: &str, seq: &str) -> String {
format!("@{name}\n{seq}\n+\n{}\n", "I".repeat(seq.len()))
}
fn interleaved_fq_text(n: usize, r1_seq: &str, r2_seq: &str) -> String {
let mut s = String::new();
for i in 0..n {
s += &fq_record(&format!("pair{i}/1"), r1_seq);
s += &fq_record(&format!("pair{i}/2"), r2_seq);
}
s
}
fn write_bytes(tmp: &TempDir, name: &str, bytes: &[u8]) -> PathBuf {
let path = tmp.path().join(name);
std::fs::write(&path, bytes).unwrap();
path
}
use crate::commands::utils::test_gzip as gzip_bytes;
fn trim_cmd(inputs: Vec<PathBuf>, outputs: Vec<PathBuf>, metrics: Option<PathBuf>) -> Trim {
Trim {
inputs,
outputs,
output_compression: OutputCompression::Auto,
threads: 2,
compression_level: 1,
metrics,
read_structures: vec![],
discard_unsupported_segments: false,
adapter_sequence: vec![],
adapter_fasta: None,
kit: vec![],
no_overlap_detection: true,
overlap_min_length: 5,
overlap_max_mismatch_rate: 0.1,
overlap_diagnostic_length: usize::MAX,
overlap_trust_max_chance: 1e-4,
adapter_min_length: 5,
adapter_mismatch_rate: 0.1,
trim_polyg: 0, trim_polyx: None,
quality_trim_3p: None,
quality_trim_5p: None,
filter_length: LengthFilter { min: 0, max: None },
filter_max_ns: None,
filter_mean_qual: None,
filter_low_qual: None,
json: None,
batch_size: 1024,
expected_insert_size: None,
insert_size_stats: false,
}
}
fn owned_rec(head: &str, seq: &str, qual: &str) -> OwnedRecord {
assert_eq!(seq.len(), qual.len());
OwnedRecord {
head: head.as_bytes().to_vec(),
seq: seq.as_bytes().to_vec(),
qual: qual.as_bytes().to_vec(),
}
}
fn rs(s: &str) -> ReadStructure {
s.parse().unwrap()
}
fn apply_rs(
rs_spec: &ReadStructure,
rec: &mut OwnedRecord,
discard_unsupported: bool,
umi_parts: &mut Vec<Vec<u8>>,
) -> Result<ApplyRsOutcome> {
let mut s = Vec::new();
let mut q = Vec::new();
apply_read_structure(rs_spec, rec, discard_unsupported, umi_parts, &mut s, &mut q)
}
fn detect_overlap(r1: &[u8], r2: &[u8], min_overlap: usize, max_mm_rate: f64) -> Option<usize> {
let mut scratch = OverlapScratch::default();
let empty_library = OverlapAdapterLibrary::default();
detect_pe_overlap(
r1,
r2,
min_overlap,
max_mm_rate,
usize::MAX,
&empty_library,
isize::MIN,
false,
None,
&mut scratch,
)
.inferred_insert
}
fn synth_pair(
molecule: &[u8],
adapter_r1: &[u8],
adapter_r2: &[u8],
r1_len: usize,
r2_len: usize,
) -> (Vec<u8>, Vec<u8>) {
let i = molecule.len();
let mut r1 = Vec::with_capacity(r1_len);
if i >= r1_len {
r1.extend_from_slice(&molecule[..r1_len]);
} else {
r1.extend_from_slice(molecule);
r1.extend_from_slice(&adapter_r1[..r1_len - i]);
}
let bottom = rc_bytes(molecule);
let mut r2 = Vec::with_capacity(r2_len);
if i >= r2_len {
r2.extend_from_slice(&bottom[..r2_len]);
} else {
r2.extend_from_slice(&bottom);
r2.extend_from_slice(&adapter_r2[..r2_len - i]);
}
(r1, r2)
}
fn make_template(len: usize, seed: u64) -> Vec<u8> {
const BASES: &[u8; 4] = b"ACGT";
let mut out = Vec::with_capacity(len);
let mut x = seed.wrapping_add(1);
for _ in 0..len {
x = x.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
out.push(BASES[((x >> 33) & 3) as usize]);
}
out
}
fn rc_bytes(seq: &[u8]) -> Vec<u8> {
seq.iter()
.rev()
.map(|&b| match b {
b'A' => b'T',
b'T' => b'A',
b'C' => b'G',
b'G' => b'C',
b'N' => b'N',
other => other,
})
.collect()
}
fn walk_overlap_test(
r1: &[u8],
r2: &[u8],
center: isize,
stats_on: bool,
adapter_library: &OverlapAdapterLibrary,
) -> Option<usize> {
let mut scratch = OverlapScratch::default();
detect_pe_overlap(
r1,
r2,
11,
0.05,
usize::MAX,
adapter_library,
center,
stats_on,
None,
&mut scratch,
)
.inferred_insert
}
fn eval_filters(
recs: &[OwnedRecord],
length: LengthFilter,
filter_max_ns: Option<usize>,
filter_mean_qual: Option<u8>,
) -> Option<FilterReject> {
let stats: Vec<BaseStats> = recs.iter().map(|r| observe_stats(&r.seq, &r.qual)).collect();
evaluate_filters(recs, &stats, length, filter_max_ns, filter_mean_qual, None)
}
#[test]
fn validation_rejects_missing_input() {
let tmp = TempDir::new().unwrap();
let cmd =
trim_cmd(vec![tmp.path().join("nope.fq")], vec![tmp.path().join("out.fq.gz")], None);
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("does not exist"), "{err}");
}
#[test]
fn sniffed_se_input_with_two_outputs_errors() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let cmd = trim_cmd(vec![r1], vec![o1.clone(), o2.clone()], None);
let err = cmd.execute().unwrap_err().to_string();
assert!(err.contains("single-end"), "{err}");
assert!(!o1.exists() && !o2.exists(), "no output should be created on a layout error");
}
#[test]
fn validation_rejects_output_overwriting_input() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let cmd = trim_cmd(vec![r1.clone()], vec![r1.clone()], None);
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("refusing to overwrite"), "{err}");
}
#[test]
fn validation_rejects_metrics_overwriting_input() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let cmd = trim_cmd(vec![r1.clone()], vec![tmp.path().join("out.fq.gz")], Some(r1));
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("refusing to overwrite"), "{err}");
}
#[test]
fn validation_rejects_json_overwriting_input() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1.clone()], vec![tmp.path().join("out.fq.gz")], None);
cmd.json = Some(r1);
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("refusing to overwrite"), "{err}");
}
#[test]
fn validation_rejects_more_than_two_inputs() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let cmd = trim_cmd(vec![r1.clone(), r1.clone(), r1], vec![tmp.path().join("o.fq")], None);
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("--inputs accepts at most 2 paths; got 3"), "{err}");
}
#[test]
fn validation_rejects_more_than_two_outputs() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let outs = (1..=3).map(|i| tmp.path().join(format!("o{i}.fq"))).collect();
let err = trim_cmd(vec![r1], outs, None).validate().unwrap_err().to_string();
assert!(err.contains("--outputs accepts at most 2 paths; got 3"), "{err}");
}
#[test]
fn validation_rejects_dotdot_aliased_duplicate_outputs() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
std::fs::create_dir(tmp.path().join("sub")).unwrap();
let outs = vec![tmp.path().join("sub/../dup.fq"), tmp.path().join("dup.fq")];
let err = trim_cmd(vec![r1], outs, None).validate().unwrap_err().to_string();
assert!(err.contains("must be distinct"), "{err}");
}
#[test]
fn validation_rejects_same_input_twice() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let outs = vec![tmp.path().join("o1.fq"), tmp.path().join("o2.fq")];
let err = trim_cmd(vec![r1.clone(), r1], outs, None).validate().unwrap_err().to_string();
assert!(err.contains("R1 and R2 must be different files"), "{err}");
}
#[test]
fn validation_rejects_duplicate_output_paths() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq.gz");
let cmd = trim_cmd(vec![r1], vec![out.clone(), out], None);
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("must be distinct"), "{err}");
}
#[test]
fn validation_rejects_output_and_metrics_same_path() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let shared = tmp.path().join("shared.out");
let cmd = trim_cmd(vec![r1], vec![shared.clone()], Some(shared));
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("must be distinct"), "{err}");
}
#[test]
fn validation_rejects_metrics_dash() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], Some(PathBuf::from("-")));
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("--metrics does not support '-'"), "{err}");
}
#[test]
fn validation_rejects_json_dash() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], None);
cmd.json = Some(PathBuf::from("-"));
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("--json does not support '-'"), "{err}");
}
#[test]
fn validation_rejects_output_parent_missing() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let cmd =
trim_cmd(vec![r1], vec![tmp.path().join("nonexistent_dir").join("out.fq.gz")], None);
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("parent directory"), "{err}");
}
#[test]
fn validation_rejects_too_few_threads() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], None);
cmd.threads = 0;
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("Threads must be at least 1"), "{err}");
}
#[test]
fn validation_rejects_bad_compression_level() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], None);
cmd.compression_level = 0;
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("Compression level"), "{err}");
}
#[test]
fn validation_aggregates_multiple_errors() {
let tmp = TempDir::new().unwrap();
let mut cmd =
trim_cmd(vec![tmp.path().join("missing.fq")], vec![tmp.path().join("out.fq.gz")], None);
cmd.threads = 0;
cmd.compression_level = 99;
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("does not exist"), "{err}");
assert!(err.contains("Threads must be at least 1"), "{err}");
assert!(err.contains("Compression level"), "{err}");
}
#[test]
fn execute_single_end_passes_through() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["ACGTACGTAC", "GGGGAAAACCCC", "TTTT"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("read", &reads));
let out = tmp.path().join("out.fq.gz");
let metrics_path = tmp.path().join("trim-metrics.txt");
let cmd = trim_cmd(vec![r1], vec![out.clone()], Some(metrics_path.clone()));
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), reads.len());
for (rec, &expected) in written.iter().zip(reads.iter()) {
assert_eq!(rec.seq.as_slice(), expected.as_bytes());
}
assert!(metrics_path.exists());
}
#[test]
fn execute_paired_end_passes_through() {
let tmp = TempDir::new().unwrap();
let r1_seqs: Vec<&str> = vec!["AAAAAAAA", "CCCCCCCC", "GGGGGGGG"];
let r2_seqs: Vec<&str> = vec!["TTTTTTTT", "TTTTTTTT", "TTTTTTTT"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &r1_seqs));
let r2 = write_fastq(&tmp, "r2", &fq_lines("r", &r2_seqs));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let cmd = trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None);
cmd.execute().unwrap();
let w1 = read_fastq(&o1);
let w2 = read_fastq(&o2);
assert_eq!(w1.len(), r1_seqs.len());
assert_eq!(w2.len(), r2_seqs.len());
for (rec, &seq) in w1.iter().zip(r1_seqs.iter()) {
assert_eq!(rec.seq.as_slice(), seq.as_bytes());
}
for (rec, &seq) in w2.iter().zip(r2_seqs.iter()) {
assert_eq!(rec.seq.as_slice(), seq.as_bytes());
}
}
#[test]
fn execute_empty_input_yields_empty_output() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &[]);
let out = tmp.path().join("out.fq.gz");
let metrics_path = tmp.path().join("trim-metrics.txt");
let cmd = trim_cmd(vec![r1], vec![out.clone()], Some(metrics_path.clone()));
cmd.execute().unwrap();
assert!(out.exists());
let written = read_fastq(&out);
assert_eq!(written.len(), 0);
assert!(metrics_path.exists());
}
#[test]
fn execute_errors_on_mismatched_pair_counts() {
let tmp = TempDir::new().unwrap();
let r1_seqs: Vec<&str> = vec!["ACGT"; 5];
let r2_seqs: Vec<&str> = vec!["ACGT"; 3];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &r1_seqs));
let r2 = write_fastq(&tmp, "r2", &fq_lines("r", &r2_seqs));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let cmd = trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None);
let err = cmd.execute().unwrap_err().to_string();
assert!(err.contains("out of sync"), "{err}");
assert!(!o1.exists(), "o1 should not exist after pair-mismatch abort");
assert!(!o2.exists(), "o2 should not exist after pair-mismatch abort");
}
#[test]
fn execute_writes_metrics_counts() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAA", "CCCCCC"]; let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let metrics_path = tmp.path().join("metrics.txt");
let cmd = trim_cmd(vec![r1], vec![out], Some(metrics_path.clone()));
cmd.execute().unwrap();
let contents = std::fs::read_to_string(&metrics_path).unwrap();
let lines: Vec<&str> = contents.lines().collect();
assert_eq!(lines.len(), 2);
assert!(lines[0].contains("reads_in"));
assert!(lines[0].contains("bases_in"));
let values: Vec<&str> = lines[1].split('\t').collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let idx = |name: &str| header.iter().position(|h| *h == name).unwrap();
assert_eq!(values[idx("reads_in")], "2");
assert_eq!(values[idx("reads_out")], "2");
assert_eq!(values[idx("bases_in")], "10");
assert_eq!(values[idx("bases_out")], "10");
}
#[test]
fn join_umi_single_part() {
assert_eq!(join_umi(&[b"AACC".to_vec()]), b"AACC");
}
#[test]
fn join_umi_two_parts() {
assert_eq!(join_umi(&[b"AACC".to_vec(), b"GGTT".to_vec()]), b"AACC-GGTT");
}
#[test]
fn join_umi_three_parts() {
assert_eq!(join_umi(&[b"AA".to_vec(), b"CC".to_vec(), b"GG".to_vec()]), b"AA-CC-GG");
}
#[test]
fn append_umi_to_short_head() {
let mut head = b"readname".to_vec();
append_umi_to_head(&mut head, b"AACCGG").unwrap();
assert_eq!(head, b"readname:AACCGG");
}
#[test]
fn append_umi_to_illumina_7_field_head() {
let mut head = b"INSTR:123:FLOWCELL:1:1101:1000:2000".to_vec();
append_umi_to_head(&mut head, b"AACCGG").unwrap();
assert_eq!(head, b"INSTR:123:FLOWCELL:1:1101:1000:2000:AACCGG");
}
#[test]
fn append_umi_extends_existing_field_8() {
let mut head = b"INSTR:123:FLOWCELL:1:1101:1000:2000:TTTT".to_vec();
append_umi_to_head(&mut head, b"AACCGG").unwrap();
assert_eq!(head, b"INSTR:123:FLOWCELL:1:1101:1000:2000:TTTT-AACCGG");
}
#[test]
fn append_umi_preserves_comment_after_space() {
let mut head = b"INSTR:123:FLOWCELL:1:1101:1000:2000 1:N:0:CTAG".to_vec();
append_umi_to_head(&mut head, b"AACCGG").unwrap();
assert_eq!(head, b"INSTR:123:FLOWCELL:1:1101:1000:2000:AACCGG 1:N:0:CTAG");
}
#[test]
fn append_umi_rejects_too_many_fields() {
let mut head = b"a:b:c:d:e:f:g:h:i".to_vec();
let err = append_umi_to_head(&mut head, b"AACCGG").unwrap_err().to_string();
assert!(err.contains("more than"), "{err}");
}
#[test]
fn apply_read_structure_template_only() {
let mut rec = owned_rec("read1", "ACGTACGT", "IIIIIIII");
let mut umis = vec![];
apply_rs(&rs("+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"ACGTACGT");
assert_eq!(rec.qual, b"IIIIIIII");
assert!(umis.is_empty());
}
#[test]
fn apply_read_structure_hard_trim_5s_plus_t() {
let mut rec = owned_rec("read1", "AAAAACCCC", "112233445");
let mut umis = vec![];
apply_rs(&rs("5S+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"CCCC");
assert_eq!(rec.qual, b"3445");
assert!(umis.is_empty());
}
#[test]
fn apply_read_structure_extracts_umi() {
let mut rec = owned_rec("read1", "AACCGGTTTTAA", "111222333444");
let mut umis = vec![];
apply_rs(&rs("4M+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"GGTTTTAA");
assert_eq!(rec.qual, b"22333444");
assert_eq!(umis, vec![b"AACC".to_vec()]);
}
#[test]
fn apply_read_structure_skip_segment() {
let mut rec = owned_rec("read1", "SSSTEMPL", "12345678");
let mut umis = vec![];
apply_rs(&rs("3S+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"TEMPL");
assert_eq!(rec.qual, b"45678");
assert!(umis.is_empty());
}
#[test]
fn apply_read_structure_too_short_returns_too_short() {
let mut rec = owned_rec("read1", "ACG", "III");
let mut umis = vec![];
let outcome = apply_rs(&rs("5M+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(outcome, ApplyRsOutcome::TooShort);
}
#[test]
fn apply_read_structure_multiple_m_segments() {
let mut rec = owned_rec("read1", "AAABBBCCCDDD", "123456789012");
let mut umis = vec![];
apply_rs(&rs("3M3S3M+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"DDD");
assert_eq!(rec.qual, b"012");
assert_eq!(umis, vec![b"AAA".to_vec(), b"CCC".to_vec()]);
}
#[test]
fn apply_read_structure_discards_b_when_permitted() {
let mut rec = owned_rec("read1", "BBBBTEMPL", "123456789");
let mut umis = vec![];
apply_rs(&rs("4B+T"), &mut rec, true, &mut umis).unwrap();
assert_eq!(rec.seq, b"TEMPL");
assert_eq!(rec.qual, b"56789");
}
#[test]
fn validation_rejects_b_segment() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], None);
cmd.read_structures = vec![rs("4B+T")];
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("sample barcode (B)"), "{err}");
assert!(err.contains("fqtk demux"), "{err}");
}
#[test]
fn validation_rejects_c_segment() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], None);
cmd.read_structures = vec![rs("4C+T")];
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("cellular barcode (C)"), "{err}");
}
#[test]
fn validation_allows_b_with_discard_flag() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGTACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("out.fq.gz")], None);
cmd.read_structures = vec![rs("4B+T")];
cmd.discard_unsupported_segments = true;
cmd.validate().unwrap();
}
#[test]
fn validation_rejects_wrong_read_structure_count() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(
vec![r1, r2],
vec![tmp.path().join("o1.fq.gz"), tmp.path().join("o2.fq.gz")],
None,
);
cmd.read_structures = vec![rs("+T")]; let err = cmd.execute().unwrap_err().to_string();
assert!(err.contains("must be 0 or equal to the number of mates"), "{err}");
}
#[test]
fn execute_se_hard_trim() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAAATEMPL", "BBBBBSHORT"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("read", &reads));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.read_structures = vec![rs("5S+T")];
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 2);
assert_eq!(written[0].seq.as_slice(), b"TEMPL");
assert_eq!(written[1].seq.as_slice(), b"SHORT");
}
#[test]
fn execute_se_umi_extraction_appends_to_head() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@A:1:B:1:1:1:1".to_string(),
"AACCGGTTTTAA".to_string(),
"+".to_string(),
"IIIIIIIIIIII".to_string(),
];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.read_structures = vec![rs("4M+T")];
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 1);
assert_eq!(written[0].head.as_slice(), b"A:1:B:1:1:1:1:AACC");
assert_eq!(written[0].seq.as_slice(), b"GGTTTTAA");
}
#[test]
fn execute_pe_umi_from_both_mates_is_joined() {
let tmp = TempDir::new().unwrap();
let r1_lines = vec![
"@A:1:B:1:1:1:1".to_string(),
"AAAAGGGGGG".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
];
let r2_lines = vec![
"@A:1:B:1:1:1:1".to_string(),
"TTTTCCCCCC".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
];
let r1 = write_fastq(&tmp, "r1", &r1_lines);
let r2 = write_fastq(&tmp, "r2", &r2_lines);
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let mut cmd = trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None);
cmd.read_structures = vec![rs("4M+T"), rs("4M+T")];
cmd.execute().unwrap();
let w1 = read_fastq(&o1);
let w2 = read_fastq(&o2);
assert_eq!(w1.len(), 1);
assert_eq!(w2.len(), 1);
assert_eq!(w1[0].head.as_slice(), b"A:1:B:1:1:1:1:AAAA-TTTT");
assert_eq!(w2[0].head.as_slice(), b"A:1:B:1:1:1:1:AAAA-TTTT");
assert_eq!(w1[0].seq.as_slice(), b"GGGGGG");
assert_eq!(w2[0].seq.as_slice(), b"CCCCCC");
}
#[test]
fn apply_read_structure_rejects_read_exactly_at_fixed_length() {
let mut rec = owned_rec("read1", "AAAA", "IIII");
let mut umis = vec![];
let outcome = apply_rs(&rs("4M+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(outcome, ApplyRsOutcome::TooShort);
}
#[test]
fn apply_read_structure_minimum_variable_length() {
let mut rec = owned_rec("read1", "AAAAT", "11112");
let mut umis = vec![];
apply_rs(&rs("4M+T"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"T");
assert_eq!(rec.qual, b"2");
assert_eq!(umis, vec![b"AAAA".to_vec()]);
}
#[test]
fn apply_read_structure_tail_skip() {
let mut rec = owned_rec("read1", "ACGTACGTACGTNNNNN", "11111111111111111");
let mut umis = vec![];
apply_rs(&rs("+T5S"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"ACGTACGTACGT");
assert_eq!(rec.qual, b"111111111111");
assert!(umis.is_empty());
}
#[test]
fn apply_read_structure_head_umi_and_tail_skip() {
let mut rec = owned_rec("read1", "UUUUUUACGTACGTNNNNN", "1111112222222233333");
let mut umis = vec![];
apply_rs(&rs("6M+T5S"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"ACGTACGT");
assert_eq!(rec.qual, b"22222222");
assert_eq!(umis, vec![b"UUUUUU".to_vec()]);
}
#[test]
fn apply_read_structure_middle_plus_requires_post_bases() {
let mut rec = owned_rec("read1", "ACGT", "IIII");
let mut umis = vec![];
let outcome = apply_rs(&rs("+T5S"), &mut rec, false, &mut umis).unwrap();
assert_eq!(outcome, ApplyRsOutcome::TooShort);
}
#[test]
fn apply_read_structure_head_skip_and_tail_skip() {
let mut rec = owned_rec("read1", "HHHHHACGTACGTACTTTTT", "11111BBBBBBBBBBCCCCC");
let mut umis = vec![];
apply_rs(&rs("5S+T5S"), &mut rec, false, &mut umis).unwrap();
assert_eq!(rec.seq, b"ACGTACGTAC");
assert_eq!(rec.qual, b"BBBBBBBBBB");
assert!(umis.is_empty());
}
#[test]
fn append_umi_to_head_extends_existing_field_8_idempotently() {
let mut head = b"A:1:B:1:1:1:1:AAAA".to_vec();
append_umi_to_head(&mut head, b"BBBB").unwrap();
assert_eq!(head, b"A:1:B:1:1:1:1:AAAA-BBBB");
}
#[test]
fn reverse_complement_basic() {
assert_eq!(reverse_complement(b"ACGT"), b"ACGT");
assert_eq!(reverse_complement(b"AAAA"), b"TTTT");
assert_eq!(reverse_complement(b"ACCGTTT"), b"AAACGGT");
}
#[test]
fn reverse_complement_preserves_n() {
assert_eq!(reverse_complement(b"ACNGT"), b"ACNGT");
}
#[test]
fn reverse_complement_empty() {
assert_eq!(reverse_complement(b""), b"");
}
#[test]
fn iupac_exact_match() {
assert!(base_matches_iupac(b'A', b'A'));
assert!(base_matches_iupac(b'C', b'C'));
assert!(!base_matches_iupac(b'A', b'C'));
}
#[test]
fn iupac_n_in_adapter_matches_any_read_base() {
assert!(base_matches_iupac(b'A', b'N'));
assert!(base_matches_iupac(b'C', b'N'));
assert!(base_matches_iupac(b'G', b'N'));
assert!(base_matches_iupac(b'T', b'N'));
}
#[test]
fn iupac_ambiguity_code_matches_compatible_reads() {
assert!(base_matches_iupac(b'A', b'R'));
assert!(base_matches_iupac(b'G', b'R'));
assert!(!base_matches_iupac(b'C', b'R'));
assert!(!base_matches_iupac(b'T', b'R'));
}
#[test]
fn iupac_case_insensitive() {
assert!(base_matches_iupac(b'a', b'A'));
assert!(base_matches_iupac(b'a', b'n'));
}
fn ad(bytes: &[u8]) -> Adapter {
Adapter::new(bytes.to_vec())
}
#[test]
fn adapter_3prime_full_adapter_at_end() {
let read = b"TTTTTTTTTTAGATCGGAAGAG";
let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, None), Some(10));
}
#[test]
fn adapter_3prime_partial_at_end() {
let read = b"TTTTTTTTTTAGAT"; let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 4, 0.1, None), Some(10));
}
#[test]
fn adapter_3prime_partial_below_min_overlap_no_match() {
let read = b"TTTTTTTTTTAGA"; let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, None), None);
}
#[test]
fn adapter_3prime_tolerates_one_mismatch() {
let read = b"TTTTTTTTTTTGATCGGAAGAG";
let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, None), Some(10));
}
#[test]
fn adapter_3prime_rejects_too_many_mismatches() {
let read = b"TTTTTTTTTTCCCCCCCC";
let adapter = ad(b"AGATCGGA");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, None), None);
}
#[test]
fn adapter_3prime_iupac_n_in_adapter() {
let read = b"TTTTTTTTTTAGATCGGAAGAG";
let adapter = ad(b"AGANCGGAAGAG"); assert!(!adapter.pure_acgt);
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.0, None), Some(10));
}
#[test]
fn adapter_3prime_pure_acgt_takes_simd_path() {
let adapter = ad(b"AGATCGGAAGAG");
assert!(adapter.pure_acgt);
}
#[test]
fn adapter_3prime_empty_adapter_is_none() {
let adapter = ad(b"");
assert_eq!(find_adapter_3prime(b"ACGT", &adapter, 5, 0.1, None), None);
}
#[test]
fn adapter_3prime_read_shorter_than_min_overlap() {
let adapter = ad(b"AGATCGG");
assert_eq!(find_adapter_3prime(b"ACGT", &adapter, 5, 0.1, None), None);
}
#[test]
fn adapter_3prime_adapter_longer_than_read() {
let read = b"AGATC";
let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, None), Some(0));
}
#[test]
fn adapter_3prime_read_entirely_adapter() {
let read = b"AGATCGGAAGAG";
let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.0, None), Some(0));
}
#[test]
fn adapter_3prime_max_k_blocks_late_matches() {
let read = b"TTTTTTTTTTTTTTTTTTTTAGATCGGAAGAG";
let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, None), Some(20));
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, Some(10)), None);
}
#[test]
fn adapter_3prime_max_k_admits_in_range_match() {
let read = b"TTTTTAGATCGGAAGAG";
let adapter = ad(b"AGATCGGAAGAG");
assert_eq!(find_adapter_3prime(read, &adapter, 5, 0.1, Some(10)), Some(5));
}
#[test]
fn adapter_3prime_finds_adapter_with_mismatches_in_its_first_16_bases() {
let mut read: Vec<u8> = b"GATTACA".iter().copied().cycle().take(40).collect();
read.extend_from_slice(b"TGAACGGAAGAGCTCACGTCTGAACTCCAGTCA");
read.extend(b"CCCCTTTT".iter().cycle().take(150 - read.len()));
let adapter = ad(b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA");
assert_eq!(find_adapter_3prime(&read, &adapter, 3, 0.125, None), Some(40));
}
fn unscreened_adapter_scan(
read: &[u8],
adapter: &Adapter,
min_length: usize,
max_mm_rate: f64,
max_k: Option<usize>,
) -> Option<usize> {
if read.len() < min_length || adapter.bytes.is_empty() {
return None;
}
let max_start = read.len() - min_length;
let max_start = max_k.map(|m| m.min(max_start)).unwrap_or(max_start);
(0..=max_start).find(|&k| {
let alignment_len = adapter.bytes.len().min(read.len() - k);
let max_mm = (alignment_len as f64 * max_mm_rate).floor() as usize;
let matches = |r: &u8, a: &u8| {
if adapter.pure_acgt {
r.eq_ignore_ascii_case(a)
} else {
base_matches_iupac(*r, *a)
}
};
alignment_len >= min_length
&& read[k..k + alignment_len]
.iter()
.zip(&adapter.bytes[..alignment_len])
.filter(|(r, a)| !matches(r, a))
.count()
<= max_mm
})
}
#[test]
fn adapter_3prime_matches_unscreened_scan_on_random_reads() {
let adapters = [
ad(b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"),
ad(b"CTGTCTCTTATACACATCT"),
ad(b"TGGAATTCTCGGGTGCCAAGG"),
ad(b"AGATCGGAAGAGC"),
ad(b"AGATCGGAAGAGCRCACG"),
];
let mut rng: u64 = 0xADA9_7E55_C0FF_EE11;
let mut next = |n: usize| {
rng ^= rng << 13;
rng ^= rng >> 7;
rng ^= rng << 17;
(rng % n as u64) as usize
};
for case in 0..20_000 {
let adapter = &adapters[next(adapters.len())];
let read_len = 1 + next(260);
let start = next(read_len + 1);
let mut read: Vec<u8> = (0..start).map(|_| b"ACGT"[next(4)]).collect();
read.extend(adapter.bytes.iter().chain(b"ACGTACGTTTTT").cycle().take(read_len - start));
let error_rate = [0, 2, 8, 20][next(4)];
for base in read.iter_mut() {
match next(100) {
r if r < error_rate => *base = b"ACGTN"[next(5)],
99 => base.make_ascii_lowercase(),
_ => {}
}
}
let min_length = 1 + next(20);
let max_mm_rate = [0.0, 0.1, 0.125, 0.2][next(4)];
let max_k = [None, Some(next(read_len + 1))][next(2)];
assert_eq!(
find_adapter_3prime(&read, adapter, min_length, max_mm_rate, max_k),
unscreened_adapter_scan(&read, adapter, min_length, max_mm_rate, max_k),
"case {case}: min_length {min_length}, rate {max_mm_rate}, max_k {max_k:?}\n \
adapter {}\n read {}",
String::from_utf8_lossy(&adapter.bytes),
String::from_utf8_lossy(&read),
);
}
}
#[test]
fn best_adapter_match_picks_earliest_trim() {
let read = b"INSERTAGAT"; let adapters = vec![ad(b"CCCC"), ad(b"AGAT"), ad(b"AGATCGG")];
assert_eq!(find_best_adapter_match(read, &adapters, 4, 0.1, None), Some(6));
}
#[test]
fn best_adapter_match_none_match() {
let adapters = vec![ad(b"TTTT")];
assert_eq!(find_best_adapter_match(b"AAAACCCC", &adapters, 4, 0.1, None), None);
}
#[test]
fn pe_overlap_short_insert() {
let r1 = b"TTTTAA";
let r2 = b"AAAACC";
assert_eq!(detect_overlap(r1, r2, 4, 0.0), Some(4));
}
#[test]
fn pe_overlap_no_overlap_when_insert_longer_than_read() {
let r1 = b"AAAAAAAA";
let r2 = b"CCCCCCCC";
assert_eq!(detect_overlap(r1, r2, 4, 0.0), None);
}
#[test]
fn pe_overlap_tolerates_mismatches() {
let r1 = b"AAAATCCC";
let r2 = b"AGTTTTGG"; assert_eq!(detect_overlap(r1, r2, 5, 0.25), Some(5));
assert_eq!(detect_overlap(r1, r2, 5, 0.1), None);
}
#[test]
fn pe_overlap_respects_min_overlap() {
let r1 = b"TTTTAA";
let r2 = b"AAAACC";
assert_eq!(detect_overlap(r1, r2, 5, 0.0), None);
}
#[test]
fn walk_overlap_full_overlap_symmetric() {
let molecule = make_template(32, 1);
let (r1, r2) = synth_pair(&molecule, b"", b"", 32, 32);
let lib = OverlapAdapterLibrary::default();
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, false, &lib), Some(32));
}
#[test]
fn walk_overlap_short_insert_with_adapter() {
let molecule = make_template(20, 2);
let adapter_r1 = make_template(16, 100);
let adapter_r2 = make_template(16, 200);
let (r1, r2) = synth_pair(&molecule, &adapter_r1, &adapter_r2, 32, 32);
let lib = OverlapAdapterLibrary::default();
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, false, &lib), Some(20));
}
#[test]
fn walk_overlap_inner_overlap_requires_stats_on() {
let molecule = make_template(50, 3);
let (r1, r2) = synth_pair(&molecule, b"", b"", 32, 32);
let lib = OverlapAdapterLibrary::default();
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, false, &lib), None);
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, true, &lib), Some(50));
}
#[test]
fn walk_overlap_asymmetric_reads_short_insert() {
let molecule = make_template(18, 4);
let adapter_r1 = make_template(14, 100);
let adapter_r2 = make_template(14, 200);
let (r1, r2) = synth_pair(&molecule, &adapter_r1, &adapter_r2, 30, 25);
assert_eq!(r1.len(), 30);
assert_eq!(r2.len(), 25);
let lib = OverlapAdapterLibrary::default();
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, false, &lib), Some(18));
}
#[test]
fn walk_overlap_asymmetric_reads_inner_overlap() {
let molecule = make_template(40, 5);
let (r1, r2) = synth_pair(&molecule, b"", b"", 30, 25);
let lib = OverlapAdapterLibrary::default();
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, false, &lib), None);
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, true, &lib), Some(40));
}
#[test]
fn walk_overlap_exhausts_when_no_overlap_exists() {
let molecule = make_template(80, 6);
let (r1, r2) = synth_pair(&molecule, b"", b"", 32, 32);
let lib = OverlapAdapterLibrary::default();
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, false, &lib), None);
assert_eq!(walk_overlap_test(&r1, &r2, isize::MIN, true, &lib), None);
}
fn default_overlap_library() -> OverlapAdapterLibrary {
let adapters = build_adapter_set(&[], &None, &[], 2).unwrap();
build_overlap_adapter_library(&[], &None, &adapters).unwrap()
}
fn walk_with_defaults(
r1: &[u8],
r2: &[u8],
center: isize,
trust: Option<f64>,
) -> Option<usize> {
walk_with_library(&default_overlap_library(), r1, r2, center, trust, false)
}
fn walk_with_library(
lib: &OverlapAdapterLibrary,
r1: &[u8],
r2: &[u8],
center: isize,
trust: Option<f64>,
stats_on: bool,
) -> Option<usize> {
detect_pe_overlap(
r1,
r2,
30,
0.10,
64,
lib,
center,
stats_on,
trust,
&mut OverlapScratch::default(),
)
.inferred_insert
}
#[test]
fn best_mismatch_pmf_sums_to_one() {
for (n, k) in [(1, 1), (5, 4), (16, 5)] {
let total: f64 = best_mismatch_pmf(n, k).iter().sum();
assert!((total - 1.0).abs() < 1e-12, "n={n} k={k}: {total}");
}
}
#[test]
fn chance_is_one_when_no_tail_bases_compared() {
assert_eq!(default_overlap_library().chance(0, 0, 0), 1.0);
}
#[test]
fn chance_matches_binomial_model() {
let lib = default_overlap_library();
let short = lib.chance(5, 5, 2);
let long = lib.chance(8, 8, 2);
assert!((3e-3..9e-3).contains(&short), "{short}");
assert!((2e-6..1e-5).contains(&long), "{long}");
}
#[test]
fn chance_grows_with_allowed_mismatches() {
let lib = default_overlap_library();
assert!(lib.chance(6, 6, 0) < lib.chance(6, 6, 1));
assert!(lib.chance(6, 6, 1) < lib.chance(6, 6, 2));
}
const SATELLITE_R1: &[u8] = b"GAATGGAATGGAATGGAACGGAACGGAACGGAACGGAACGGAAAGAAATGCAATGGAATGGAATGGAATGGAACGGAACGGAAAGGAATGGAATGGAATGGAA";
const SATELLITE_R2: &[u8] = b"TTCCATTCCATTCCATTCCTTTCCGTTCCGTTCCATTCCATTCCATTCCATTGCATTTCTTTCCGTTCCGTTCCGTTCCGTTCCGTTCCATTCCATTCCATTC";
#[test]
fn first_acceptable_overlap_on_satellite_depends_on_walk_start() {
assert_eq!(walk_with_defaults(SATELLITE_R1, SATELLITE_R2, isize::MIN, None), Some(98));
assert_eq!(walk_with_defaults(SATELLITE_R1, SATELLITE_R2, 0, None), Some(103));
}
#[test]
fn trusted_overlap_on_satellite_is_the_full_overlap_from_any_walk_start() {
for center in [isize::MIN, -5, 0] {
assert_eq!(
walk_with_defaults(SATELLITE_R1, SATELLITE_R2, center, Some(1e-4)),
Some(103),
"center {center}"
);
}
}
fn telomere_pair() -> (Vec<u8>, Vec<u8>) {
let template: Vec<u8> = b"CCCTAA".iter().cycle().take(93).copied().collect();
let mut r2 = rc_bytes(&template);
r2.extend_from_slice(b"CTGTCTCTTATACACATCTGACGCTGCCGACGAGTTCTGTCATGTGTAGATCTCGGTG");
(template, r2)
}
#[test]
fn trusted_overlap_on_telomere_is_the_true_insert_from_any_walk_start() {
let (r1, r2) = telomere_pair();
let one_unit_off = 99 - r2.len() as isize;
assert_eq!(walk_with_defaults(&r1, &r2, one_unit_off, None), Some(99));
for center in [isize::MIN, one_unit_off, 0] {
assert_eq!(
walk_with_defaults(&r1, &r2, center, Some(1e-4)),
Some(93),
"center {center}"
);
}
}
fn adapter_like_repeat_pair() -> (Vec<u8>, Vec<u8>) {
let template: Vec<u8> = b"TCCGATCT".iter().cycle().take(100).copied().collect();
let adapter = b"AGATCGGAAG";
let r1 = [&template[..], &adapter[..]].concat();
let mut r2_rc = [&reverse_complement(adapter)[..], &template[..]].concat();
for base in &mut r2_rc[66..72] {
*base = if *base == b'A' { b'C' } else { b'A' };
}
(r1, reverse_complement(&r2_rc))
}
#[test]
fn trustworthy_overlap_outranks_a_better_aligned_untrustworthy_one_from_any_walk_start() {
let (r1, r2) = adapter_like_repeat_pair();
assert_eq!(walk_with_defaults(&r1, &r2, -2, None), Some(108));
for center in [isize::MIN, -10, -2, 0] {
assert_eq!(
walk_with_defaults(&r1, &r2, center, Some(1e-4)),
Some(100),
"center {center}"
);
}
}
fn unscreened_walk(
lib: &OverlapAdapterLibrary,
r1: &[u8],
r2: &[u8],
center: isize,
trust: Option<f64>,
stats_on: bool,
) -> Option<usize> {
if r1.len() < 30 || r2.len() < 30 {
return None;
}
let mut r2_rc = Vec::new();
reverse_complement_acgt_into(r2, &mut r2_rc);
let lo = -((r2.len() - 30) as isize);
let hi = if stats_on { (r1.len() - 30) as isize } else { 0 };
let probe = |shift| {
let outcome = if shift <= 0 {
try_shift_neg(r1, r2, &r2_rc, shift, 0.10, 64, lib)
} else {
try_shift_pos(r1, &r2_rc, shift, 0.10, 64)
};
match outcome {
ProbeOutcome::Accept(accepted) => Some(accepted),
_ => None,
}
};
let c = center.clamp(lo, hi);
let mut order = vec![c];
for k in 1..=(r1.len() + r2.len()) as isize {
order.extend([c - k, c + k].into_iter().filter(|s| (lo..=hi).contains(s)));
}
let first = order.into_iter().find_map(&probe)?;
match trust {
Some(max_chance) if !first.trustworthy(max_chance) => {
let best = (lo..=hi).filter_map(&probe).reduce(|best, other| {
if other.better_than(&best, max_chance) { other } else { best }
});
Some(best.unwrap_or(first).insert)
}
_ => Some(first.insert),
}
}
#[test]
fn screened_walk_matches_unscreened_walk_on_random_pairs() {
const ADAPTER_R1: &[u8] = b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA";
const ADAPTER_R2: &[u8] = b"AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT";
let lib = default_overlap_library();
let mut rng: u64 = 0x5EED_0FF5_E7C0_FFEE;
let mut next = |n: usize| {
rng ^= rng << 13;
rng ^= rng >> 7;
rng ^= rng << 17;
(rng % n as u64) as usize
};
for case in 0..5_000 {
let read_len = 40 + next(220);
let insert = 5 + next(2 * read_len);
let unit: Vec<u8> = (0..1 + next(6)).map(|_| b"ACGT"[next(4)]).collect();
let fragment: Vec<u8> = if next(4) == 0 {
unit.iter().copied().cycle().take(insert).collect()
} else {
(0..insert).map(|_| b"ACGT"[next(4)]).collect()
};
let mut read = |template: Vec<u8>, adapter: &[u8]| -> Vec<u8> {
let mut seq: Vec<u8> =
template.into_iter().chain(adapter.iter().copied()).collect();
seq.resize(read_len.max(seq.len()), b'A');
seq.truncate(read_len - next(read_len / 4));
let error_rate = [0, 1, 5, 12][next(4)];
for base in seq.iter_mut() {
match next(100) {
r if r < error_rate => *base = b"ACGTN"[next(5)],
99 => base.make_ascii_lowercase(),
_ => {}
}
}
seq
};
let r1 = read(fragment.clone(), ADAPTER_R1);
let r2 = read(reverse_complement(&fragment), ADAPTER_R2);
let center = next(r2.len() + 20) as isize - r2.len() as isize;
for (trust, stats_on) in
[(None, false), (Some(1e-4), false), (None, true), (Some(1e-4), true)]
{
assert_eq!(
walk_with_library(&lib, &r1, &r2, center, trust, stats_on),
unscreened_walk(&lib, &r1, &r2, center, trust, stats_on),
"case {case}: center {center}, trust {trust:?}, stats {stats_on}\n r1 {}\n \
r2 {}",
String::from_utf8_lossy(&r1),
String::from_utf8_lossy(&r2),
);
}
}
}
#[test]
fn screened_walk_finds_an_overlap_far_from_center() {
let fragment = b"TTGACCGATAGCTTACGGATCCAGTTAGCAATGCCTGAAC";
let mut r1 = [&fragment[..], b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"].concat();
let mut r2 =
[&reverse_complement(fragment)[..], b"AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT"].concat();
r1.resize(150, b'A');
r2.resize(150, b'A');
assert_eq!(walk_with_defaults(&r1, &r2, 0, None), Some(40));
}
#[test]
fn screen_is_unavailable_when_a_probe_is_shorter_than_16_bases() {
let r1 = vec![b'A'; 100];
let r2_rc = vec![b'A'; 100];
let mut screen = NegShiftScreen::default();
assert!(!screen.ensure(&r1, &r2_rc, 88, 0.10, 64));
screen.reset();
assert!(screen.ensure(&r1, &r2_rc, 84, 0.10, 64));
}
#[test]
fn insert_size_stats_empty_aggregate_yields_none() {
assert!(InsertSizeStats::from_aggregate(Vec::new(), 0).is_none());
assert!(InsertSizeStats::from_aggregate(vec![0, 0, 0], 0).is_none());
}
#[test]
fn insert_size_stats_only_unknowns_still_emits() {
let stats = InsertSizeStats::from_aggregate(vec![0, 0, 0], 1234).unwrap();
assert_eq!(stats.unknown, 1234);
assert_eq!(stats.peak, 0); assert_eq!(stats.histogram, vec![0, 0, 0]);
}
#[test]
fn insert_size_stats_picks_max_count_as_peak() {
let mut hist = vec![0u64; 10];
hist[3] = 100;
hist[5] = 250;
hist[7] = 175;
let stats = InsertSizeStats::from_aggregate(hist, 42).unwrap();
assert_eq!(stats.peak, 5);
assert_eq!(stats.unknown, 42);
}
#[test]
fn build_adapter_set_from_kit() {
let set = build_adapter_set(&[], &None, &["truseq".to_string()], 2).unwrap();
assert!(!set.r1.is_empty());
assert!(!set.r2.is_empty());
assert!(set.r1.iter().any(|a| a.bytes == b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"));
assert!(set.r2.iter().any(|a| a.bytes == b"AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT"));
}
#[test]
fn build_adapter_set_explicit_sequences_lowercase_normalized() {
let set =
build_adapter_set(&["agat".to_string(), "cgat".to_string()], &None, &[], 2).unwrap();
let r1_bytes: Vec<Vec<u8>> = set.r1.iter().map(|a| a.bytes.clone()).collect();
let r2_bytes: Vec<Vec<u8>> = set.r2.iter().map(|a| a.bytes.clone()).collect();
assert_eq!(r1_bytes, vec![b"AGAT".to_vec()]);
assert_eq!(r2_bytes, vec![b"CGAT".to_vec()]);
}
#[test]
fn build_adapter_set_dedupes() {
let set = build_adapter_set(
&["AGATCGGAAGAGCACACGTCTGAACTCCAGTCA".to_string()],
&None,
&["truseq".to_string()],
1,
)
.unwrap();
let count =
set.r1.iter().filter(|a| a.bytes == b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA").count();
assert_eq!(count, 1);
}
#[test]
fn build_adapter_set_classifies_pure_acgt() {
let set = build_adapter_set(&["AGATCGG".to_string()], &None, &[], 1).unwrap();
assert!(set.r1.iter().all(|a| a.pure_acgt));
let set = build_adapter_set(&["AGANCGG".to_string()], &None, &[], 1).unwrap();
assert!(!set.r1[0].pure_acgt);
}
#[test]
fn build_adapter_set_is_empty_without_configuration() {
let set = build_adapter_set(&[], &None, &[], 2).unwrap();
assert!(set.is_empty());
}
#[test]
fn build_adapter_set_small_rna_leaves_r2_empty() {
let set = build_adapter_set(&[], &None, &["small-rna".to_string()], 2).unwrap();
assert_eq!(set.r1.len(), 1);
assert!(set.r2.is_empty());
}
#[test]
fn build_adapter_set_all_unions_every_kit() {
let set = build_adapter_set(&[], &None, &["all".to_string()], 2).unwrap();
assert_eq!(set.r1.len(), 5);
assert_eq!(set.r2.len(), 4);
}
#[test]
fn reverse_complement_iupac_codes() {
assert_eq!(reverse_complement(b"RMKY"), b"RMKY");
assert_eq!(reverse_complement(b"ARGC"), b"GCYT");
}
#[test]
fn validation_rejects_bad_adapter_mismatch_rate() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("o.fq.gz")], None);
cmd.adapter_mismatch_rate = 1.5;
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("adapter-mismatch-rate"), "{err}");
}
#[test]
fn validation_rejects_unknown_kit() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("o.fq.gz")], None);
cmd.kit = vec!["unicorn".to_string()];
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("not a recognized preset"), "{err}");
}
#[test]
fn validation_rejects_invalid_adapter_base() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let mut cmd = trim_cmd(vec![r1], vec![tmp.path().join("o.fq.gz")], None);
cmd.adapter_sequence = vec!["ACGZ".to_string()];
let err = cmd.validate().unwrap_err().to_string();
assert!(err.contains("invalid base"), "{err}");
}
#[test]
fn execute_se_trims_explicit_adapter() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAAAAAAAAAGATCGGAAGAG", "CCCCCCCCCCAGATCGGAAGAG"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.adapter_sequence = vec!["AGATCGGAAGAG".to_string()];
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 2);
assert_eq!(written[0].seq.as_slice(), b"AAAAAAAAAA");
assert_eq!(written[1].seq.as_slice(), b"CCCCCCCCCC");
}
#[test]
fn execute_se_trims_truseq_kit() {
let tmp = TempDir::new().unwrap();
let truseq = b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA";
let insert = b"TTTTTTTTTT";
let mut read = Vec::from(&insert[..]);
read.extend_from_slice(truseq);
let reads_str = String::from_utf8(read).unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &[reads_str.as_str()]));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.kit = vec!["truseq".to_string()];
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written[0].seq.as_slice(), insert);
}
#[test]
fn execute_pe_overlap_trims_both_mates() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("p", &["ACGTACGTAGAT"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("p", &["ACGTACGTAGAT"]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let mut cmd = trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None);
cmd.no_overlap_detection = false;
cmd.overlap_min_length = 5;
cmd.execute().unwrap();
let w1 = read_fastq(&o1);
let w2 = read_fastq(&o2);
assert_eq!(w1[0].seq.as_slice(), b"ACGTACGT");
assert_eq!(w2[0].seq.as_slice(), b"ACGTACGT");
}
#[test]
fn execute_pe_overlap_evidence_check_fires_for_custom_adapter() {
let tmp = TempDir::new().unwrap();
let insert_r1_half = b"AAAACCCCGGGGTTTTAAAA"; let insert_r2_half = b"TTTTAAAACCCCGGGGTTTT"; let adapter_r1 = b"CAGCAGATCTCGGTGG"; let adapter_r2 = b"GGAGATCAGCAGTCGC"; let r1_seq: Vec<u8> = [&insert_r1_half[..], &adapter_r1[..], &[b'N'; 4][..]].concat();
let r2_seq: Vec<u8> = [&insert_r2_half[..], &adapter_r2[..], &[b'N'; 4][..]].concat();
let r1_path =
write_fastq(&tmp, "r1", &fq_lines("p", &[std::str::from_utf8(&r1_seq).unwrap()]));
let r2_path =
write_fastq(&tmp, "r2", &fq_lines("p", &[std::str::from_utf8(&r2_seq).unwrap()]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let mut cmd = trim_cmd(vec![r1_path, r2_path], vec![o1.clone(), o2.clone()], None);
cmd.no_overlap_detection = false;
cmd.overlap_min_length = 10;
cmd.adapter_sequence = vec![
std::str::from_utf8(adapter_r1).unwrap().to_string(),
std::str::from_utf8(adapter_r2).unwrap().to_string(),
];
cmd.execute().unwrap();
let w1 = read_fastq(&o1);
let w2 = read_fastq(&o2);
assert_eq!(w1.len(), 1);
assert_eq!(w2.len(), 1);
assert!(w1[0].seq.len() < 40, "R1 should be trimmed, got {}", w1[0].seq.len());
assert!(w2[0].seq.len() < 40, "R2 should be trimmed, got {}", w2[0].seq.len());
}
#[test]
fn execute_adapter_trim_tracks_bases_metric() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAAAAAAAAAGATCGGAAGAG"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let metrics_path = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1], vec![out], Some(metrics_path.clone()));
cmd.adapter_sequence = vec!["AGATCGGAAGAG".to_string()];
cmd.execute().unwrap();
let contents = std::fs::read_to_string(&metrics_path).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let idx = |name: &str| header.iter().position(|h| *h == name).unwrap();
assert_eq!(values[idx("bases_in")], "22");
assert_eq!(values[idx("bases_out")], "10");
assert_eq!(values[idx("bases_trimmed_adapter")], "12");
}
#[test]
fn execute_fasta_adapter_trimming() {
let tmp = TempDir::new().unwrap();
let fasta_path = tmp.path().join("adapters.fa");
let fasta_lines = vec![
">adapter1".to_string(),
"AGATCGGAAGAG".to_string(),
">adapter2".to_string(),
"TGTCTCTTATAC".to_string(),
];
Io::default().write_lines(&fasta_path, &fasta_lines).unwrap();
let reads: Vec<&str> = vec!["AAAAAAAAAAAGATCGGAAGAG", "CCCCCCCCCCTGTCTCTTATAC"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.adapter_fasta = Some(fasta_path);
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written[0].seq.as_slice(), b"AAAAAAAAAA");
assert_eq!(written[1].seq.as_slice(), b"CCCCCCCCCC");
}
#[test]
fn execute_hard_trim_tracks_bases() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAAATEMPL"]; let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let metrics_path = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1], vec![out], Some(metrics_path.clone()));
cmd.read_structures = vec![rs("5S+T")];
cmd.execute().unwrap();
let contents = std::fs::read_to_string(&metrics_path).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let idx = |name: &str| header.iter().position(|h| *h == name).unwrap();
assert_eq!(values[idx("bases_in")], "10");
assert_eq!(values[idx("bases_out")], "5");
assert_eq!(values[idx("bases_trimmed_read_structure")], "5");
}
#[test]
fn quality_trim_parse_ok() {
let qt: QualityTrim = "4:20".parse().unwrap();
assert_eq!(qt.window, 4);
assert_eq!(qt.threshold, 20);
}
#[test]
fn quality_trim_parse_rejects_missing_colon() {
let err: String = "4".parse::<QualityTrim>().unwrap_err();
assert!(err.contains("WINDOW:QUAL"), "{err}");
}
#[test]
fn quality_trim_parse_rejects_zero_window() {
assert!("0:20".parse::<QualityTrim>().is_err());
}
#[test]
fn length_filter_parse_min_only() {
let lf: LengthFilter = "15".parse().unwrap();
assert_eq!(lf.min, 15);
assert!(lf.max.is_none());
}
#[test]
fn length_filter_parse_min_max() {
let lf: LengthFilter = "15:150".parse().unwrap();
assert_eq!(lf.min, 15);
assert_eq!(lf.max, Some(150));
}
#[test]
fn length_filter_parse_rejects_max_less_than_min() {
let err: String = "50:10".parse::<LengthFilter>().unwrap_err();
assert!(err.contains("max"), "{err}");
}
#[test]
fn polyx_tail_all_g() {
assert_eq!(find_polyx_tail_len(b"ACGTGGGGGGGG", b'G'), 8);
}
#[test]
fn polyx_tail_no_g() {
assert_eq!(find_polyx_tail_len(b"ACGTACGTACGT", b'G'), 0);
}
#[test]
fn polyx_tail_strict_stops_at_first_non_g() {
assert_eq!(find_polyx_tail_len(b"ACGTGGGGTGGGGGGG", b'G'), 7);
}
#[test]
fn polyx_tail_single_trailing_base() {
assert_eq!(find_polyx_tail_len(b"AAAAAAG", b'G'), 1);
}
#[test]
fn polyx_tail_empty() {
assert_eq!(find_polyx_tail_len(b"", b'G'), 0);
}
#[test]
fn trim_polyg_tail_applied() {
let mut rec = owned_rec("r", "ACGTGGGGGGGGGGGGGGGGG", "IIIIIIIIIIIIIIIIIIIII");
let trimmed = trim_polyx_tail(&mut rec, b'G', 10);
assert_eq!(trimmed, 17);
assert_eq!(rec.seq, b"ACGT");
}
#[test]
fn trim_polyg_tail_below_min_run_untouched() {
let mut rec = owned_rec("r", "ACGTGGGG", "IIIIIIII");
let trimmed = trim_polyx_tail(&mut rec, b'G', 10);
assert_eq!(trimmed, 0);
assert_eq!(rec.seq, b"ACGTGGGG");
}
#[test]
fn quality_sliding_trims_tail() {
let mut rec = owned_rec("r", "AAAAAAAAAACCCCCCCCCC", "IIIIIIIIII&&&&&&&&&&");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 11);
assert_eq!(rec.seq, b"AAAAAAAAA");
}
#[test]
fn quality_sliding_no_trim_for_high_quality_read() {
let mut rec = owned_rec("r", "AAAAAAAA", "IIIIIIII");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 0);
}
#[test]
fn quality_sliding_shorter_than_window_is_noop() {
let mut rec = owned_rec("r", "AAA", "III");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 0);
}
#[test]
fn quality_sliding_window_equals_read_length() {
let mut rec = owned_rec("r", "AAAA", "IIII");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 0);
}
#[test]
fn quality_sliding_preserves_good_bases_past_mid_dip() {
let mut rec =
owned_rec("r", "AAAAAAAAAACCCCAAAAAAAAAAAAAAAA", "IIIIIIIIII&&&&IIIIIIIIIIIIIIII");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 0, "mid-dip must not trim trailing good bases");
assert_eq!(rec.seq.len(), 30);
}
#[test]
fn quality_sliding_trims_only_trailing_bad_with_mixed_interior() {
let mut rec = owned_rec("r", "AAAAAACCCAAAAAAACCCC", "IIIIII&&&IIIIIII&&&&");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 5);
assert_eq!(rec.seq.len(), 15);
}
#[test]
fn quality_sliding_window_8_trims_tail() {
let seq = "A".repeat(50);
let qual = format!("{}{}", "I".repeat(40), "&".repeat(10));
let mut rec = owned_rec("r", &seq, &qual);
let trimmed = trim_quality_sliding_3prime(&mut rec, 8, 20);
assert_eq!(trimmed, 13);
assert_eq!(rec.seq.len(), 37);
}
#[test]
fn quality_sliding_all_bad_drops_everything() {
let mut rec = owned_rec("r", "AAAAAAAAAAAAAAAAAAAA", "&&&&&&&&&&&&&&&&&&&&");
let trimmed = trim_quality_sliding_3prime(&mut rec, 4, 20);
assert_eq!(trimmed, 20);
assert_eq!(rec.seq.len(), 0);
}
#[test]
fn quality_sliding_scalar_matches_simd_on_random_data() {
fn reference(qual: &[u8], window: usize, threshold: u8) -> usize {
let mut trim_pos = qual.len();
if qual.len() < window {
return trim_pos;
}
let max_s = qual.len() - window;
let threshold_total = (threshold as u32) * (window as u32);
let mut s = max_s as i64;
while s >= 0 {
let start = s as usize;
let sum: u32 = qual[start..start + window]
.iter()
.map(|&q| u32::from(q.saturating_sub(33)))
.sum();
if sum < threshold_total {
trim_pos = start;
s -= 1;
} else {
break;
}
}
trim_pos
}
let mut rng: u64 = 0xC0FFEE1234CAFE99;
for &window in &[1usize, 4, 8, 16] {
for _ in 0..100 {
rng ^= rng << 13;
rng ^= rng >> 7;
rng ^= rng << 17;
let len = window + 20 + (rng as usize % 25);
let mut qual: Vec<u8> = Vec::with_capacity(len);
for _ in 0..len {
rng ^= rng << 13;
rng ^= rng >> 7;
rng ^= rng << 17;
qual.push(33 + (rng as u8 % 41));
}
let seq = vec![b'A'; len];
let mut rec = owned_rec(
"r",
std::str::from_utf8(&seq).unwrap(),
std::str::from_utf8(&qual).unwrap(),
);
let expected_trim_pos = reference(&qual, window, 20);
let _trimmed = trim_quality_sliding_3prime(&mut rec, window, 20);
assert_eq!(
rec.seq.len(),
expected_trim_pos,
"trim_pos mismatch on window={window} len={len} qual={qual:?}"
);
}
}
}
fn observe_stats_scalar(seq: &[u8], qual: &[u8]) -> BaseStats {
const PHRED33: u8 = 33;
assert_eq!(seq.len(), qual.len());
let mut s = BaseStats { total: qual.len() as u64, ..Default::default() };
for (&q, &b) in qual.iter().zip(seq.iter()) {
let phred = q.saturating_sub(PHRED33);
if phred >= 20 {
s.q20 += 1;
}
if phred >= 30 {
s.q30 += 1;
}
let b_lc = b | 0x20;
if b_lc == b'n' {
s.n_bases += 1;
}
if b_lc == b'g' || b_lc == b'c' {
s.gc += 1;
}
}
s
}
#[test]
fn observe_stats_empty_is_default() {
assert_eq!(observe_stats(b"", b""), BaseStats::default());
}
#[test]
fn observe_stats_single_byte() {
assert_eq!(
observe_stats(b"A", b"5"),
BaseStats { total: 1, q20: 1, q30: 0, n_bases: 0, gc: 0 }
);
assert_eq!(
observe_stats(b"N", b"?"),
BaseStats { total: 1, q20: 1, q30: 1, n_bases: 1, gc: 0 }
);
assert_eq!(
observe_stats(b"n", b"!"),
BaseStats { total: 1, q20: 0, q30: 0, n_bases: 1, gc: 0 }
);
assert_eq!(
observe_stats(b"G", b"I"),
BaseStats { total: 1, q20: 1, q30: 1, n_bases: 0, gc: 1 }
);
assert_eq!(
observe_stats(b"c", b"!"),
BaseStats { total: 1, q20: 0, q30: 0, n_bases: 0, gc: 1 }
);
}
#[test]
fn observe_stats_exact_chunk_boundary() {
let seq = b"ACGTACGTNNNNACGT";
let qual = b"IIIIIIIIIIIIIIII";
let s = observe_stats(seq, qual);
assert_eq!(s.total, 16);
assert_eq!(s.q20, 16);
assert_eq!(s.q30, 16);
assert_eq!(s.n_bases, 4);
}
#[test]
fn observe_stats_one_past_chunk() {
let seq = b"ACGTACGTNNNNACGTN";
let qual = b"IIIIIIIIIIIIIIII!"; let s = observe_stats(seq, qual);
assert_eq!(s.total, 17);
assert_eq!(s.q20, 16);
assert_eq!(s.q30, 16);
assert_eq!(s.n_bases, 5);
}
#[test]
fn observe_stats_sub_chunk_length() {
let seq = b"NNACG";
let qual = b"!!III";
let s = observe_stats(seq, qual);
assert_eq!(s, BaseStats { total: 5, q20: 3, q30: 3, n_bases: 2, gc: 2 });
}
#[test]
fn observe_stats_pe150_length_matches_scalar() {
let seq: Vec<u8> =
(0..150).map(|i| if i % 13 == 0 { b'N' } else { b"ACGT"[i % 4] }).collect();
let qual: Vec<u8> = (0..150).map(|i| 33 + ((i * 7 + 3) % 41) as u8).collect();
assert_eq!(observe_stats(&seq, &qual), observe_stats_scalar(&seq, &qual));
}
#[test]
fn observe_stats_random_lengths_match_scalar() {
let mut state: u64 = 0x1234_5678_9abc_def0;
let mut rng = || {
state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
state
};
for len in 0..=200 {
let seq: Vec<u8> = (0..len).map(|_| b"ACGTNn"[(rng() % 6) as usize]).collect();
let qual: Vec<u8> = (0..len).map(|_| 33 + (rng() % 45) as u8).collect();
assert_eq!(
observe_stats(&seq, &qual),
observe_stats_scalar(&seq, &qual),
"mismatch at len={len}"
);
}
}
#[test]
fn observe_stats_counts_reads_longer_than_255_chunks() {
for len in [255 * 16 - 1, 255 * 16, 255 * 16 + 1, 2 * 255 * 16 + 7] {
let seq = vec![b'G'; len];
let qual = vec![b'I'; len];
let stats = observe_stats(&seq, &qual);
assert_eq!(stats, observe_stats_scalar(&seq, &qual), "len={len}");
assert_eq!(stats.gc, len as u64, "len={len}");
assert_eq!(stats.q30, len as u64, "len={len}");
}
}
fn ci_mismatches_scalar(a: &[u8], b: &[u8], limit: usize) -> usize {
assert_eq!(a.len(), b.len());
let mut count = 0;
for (&x, &y) in a.iter().zip(b.iter()) {
if !x.eq_ignore_ascii_case(&y) {
count += 1;
if count > limit {
return count;
}
}
}
count
}
#[test]
fn ci_bounded_identical_returns_zero() {
let a = b"ACGTACGTACGTACGT";
assert_eq!(count_mismatches_ci_bounded(a, a, 0), 0);
assert_eq!(count_mismatches_ci_bounded(a, a, 10), 0);
}
#[test]
fn ci_bounded_case_insensitive() {
let a = b"ACGTACGTACGTACGT";
let b = b"acgtacgtacgtacgt";
assert_eq!(count_mismatches_ci_bounded(a, b, 0), 0);
}
#[test]
fn ci_bounded_all_mismatch_stops_early() {
let a = b"AAAAAAAAAAAAAAAA";
let b = b"TTTTTTTTTTTTTTTT";
let mm = count_mismatches_ci_bounded(a, b, 3);
assert!(mm > 3, "got {mm}");
}
#[test]
fn ci_bounded_empty_is_zero() {
assert_eq!(count_mismatches_ci_bounded(b"", b"", 0), 0);
}
#[test]
fn ci_bounded_sub_chunk_tail() {
let a = b"AAGCT";
let b = b"ATGCT"; assert_eq!(count_mismatches_ci_bounded(a, b, 5), 1);
let mm = count_mismatches_ci_bounded(a, b, 0);
assert!(mm > 0);
}
#[test]
fn ci_bounded_boundary_lengths_match_scalar() {
let base_a = b"ACGTACGTACGTACGTNNNNAAAAACCCCTGG";
let base_b = b"ACGTACGTacgTACGTnnnnAAAATCCCCTGG"; for len in 0..=base_a.len() {
let a = &base_a[..len];
let b = &base_b[..len];
for limit in [0, 1, 3, 5, 16, 64] {
assert_eq!(
count_mismatches_ci_bounded(a, b, limit),
ci_mismatches_scalar(a, b, limit),
"len={len} limit={limit}"
);
}
}
}
fn rc_acgt(seq: &[u8]) -> Vec<u8> {
let mut out = Vec::new();
reverse_complement_acgt_into(seq, &mut out);
out
}
#[test]
fn rc_acgt_empty() {
assert_eq!(rc_acgt(b""), b"");
}
#[test]
fn rc_acgt_palindrome() {
assert_eq!(rc_acgt(b"ACGT"), b"ACGT");
}
#[test]
fn rc_acgt_simple_uppercase() {
assert_eq!(rc_acgt(b"AAAA"), b"TTTT");
assert_eq!(rc_acgt(b"CCCC"), b"GGGG");
assert_eq!(rc_acgt(b"ACCGTTT"), b"AAACGGT");
}
#[test]
fn rc_acgt_preserves_n() {
assert_eq!(rc_acgt(b"ACNGT"), b"ACNGT");
}
#[test]
fn rc_acgt_preserves_lowercase() {
assert_eq!(rc_acgt(b"aaaa"), b"tttt");
assert_eq!(rc_acgt(b"aCgT"), b"AcGt");
assert_eq!(rc_acgt(b"nACG"), b"CGTn");
}
#[test]
fn rc_acgt_exact_chunk_boundary() {
let seq = b"AAAACCCCGGGGTTTT";
assert_eq!(rc_acgt(seq), b"AAAACCCCGGGGTTTT");
}
#[test]
fn rc_acgt_one_past_chunk() {
let seq = b"AAAACCCCGGGGTTTTN";
assert_eq!(rc_acgt(seq), b"NAAAACCCCGGGGTTTT");
}
#[test]
fn rc_acgt_thirty_three_bytes() {
let seq = b"AAAACCCCGGGGTTTTACGTACGTACGTACGTN";
assert_eq!(rc_acgt(seq), b"NACGTACGTACGTACGTAAAACCCCGGGGTTTT");
}
#[test]
fn rc_acgt_iupac_lossy_as_documented() {
assert_eq!(rc_acgt(b"R"), b"N");
}
#[test]
fn rc_acgt_matches_scalar_on_acgt_n_inputs() {
let inputs: &[&[u8]] = &[
b"",
b"A",
b"AC",
b"ACG",
b"ACGT",
b"ACGTN",
b"ACGTACGTAC",
b"ACGTACGTACGTACGT",
b"ACGTACGTACGTACGTA",
b"NNNNAAAACCCCGGGGTTTTacgt",
];
for input in inputs {
assert_eq!(rc_acgt(input), reverse_complement(input), "input = {:?}", input);
}
}
#[test]
fn filters_length_reject_min() {
let recs = vec![owned_rec("r", "AAAA", "IIII")];
let res = eval_filters(&recs, LengthFilter { min: 10, max: None }, None, None);
assert_eq!(res, Some(FilterReject::Length));
}
#[test]
fn filters_length_reject_max() {
let recs = vec![owned_rec("r", "AAAAAAAAAA", "IIIIIIIIII")];
let res = eval_filters(&recs, LengthFilter { min: 0, max: Some(5) }, None, None);
assert_eq!(res, Some(FilterReject::Length));
}
#[test]
fn filters_n_base_rejects_over_limit() {
let recs = vec![owned_rec("r", "ANNNNN", "IIIIII")];
let res = eval_filters(&recs, LengthFilter { min: 0, max: None }, Some(2), None);
assert_eq!(res, Some(FilterReject::NBases));
}
#[test]
fn filters_mean_quality_rejects_low() {
let recs = vec![owned_rec("r", "AAAA", "!!!!")];
let res = eval_filters(&recs, LengthFilter { min: 0, max: None }, None, Some(20));
assert_eq!(res, Some(FilterReject::Quality));
}
#[test]
fn filters_none_passes() {
let recs = vec![owned_rec("r", "AAAAAAAAAA", "IIIIIIIIII")];
let res = eval_filters(&recs, LengthFilter { min: 5, max: Some(20) }, Some(5), Some(20));
assert!(res.is_none());
}
#[test]
fn filters_pe_any_mate_fails_rejects() {
let recs = vec![owned_rec("r", "AAAAAAAAAA", "IIIIIIIIII"), owned_rec("r", "AAA", "III")];
let res = eval_filters(&recs, LengthFilter { min: 5, max: None }, None, None);
assert_eq!(res, Some(FilterReject::Length));
}
#[test]
fn execute_polyg_trims_tail_when_enabled() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@NB501234:1:FLOWCELL:1:1101:1000:2000".to_string(),
"ACGTGGGGGGGGGGGGGGGGG".to_string(), "+".to_string(),
"I".repeat(21),
];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let cmd = trim_cmd(vec![r1], vec![out.clone()], None);
let mut cmd = cmd;
cmd.trim_polyg = 10;
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written[0].seq.as_slice(), b"ACGT");
}
#[test]
fn execute_read_structure_head_and_tail_skip_pe() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("p", &["AAACGTACGTTTT", "AAAGATCGATCTT"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("p", &["AAACGTACGTTTT", "AAAGATCGATCTT"]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let metrics = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], Some(metrics.clone()));
cmd.read_structures = vec![rs("3S+T3S"), rs("3S+T3S")];
cmd.execute().unwrap();
let w1 = read_fastq(&o1);
let w2 = read_fastq(&o2);
assert_eq!(w1.len(), 2);
assert_eq!(w2.len(), 2);
assert_eq!(w1[0].seq.as_slice(), b"CGTACGT");
assert_eq!(w1[1].seq.as_slice(), b"GATCGAT");
assert_eq!(w2[0].seq.as_slice(), b"CGTACGT");
assert_eq!(w2[1].seq.as_slice(), b"GATCGAT");
let contents = std::fs::read_to_string(&metrics).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let get = |name: &str| -> u64 {
let idx = header.iter().position(|h| *h == name).unwrap();
values[idx].parse().unwrap()
};
assert_eq!(get("reads_in"), 2);
assert_eq!(get("reads_out"), 2);
assert_eq!(get("reads_filtered_length"), 0);
assert_eq!(get("bases_in"), 52);
assert_eq!(get("bases_out"), 28);
assert_eq!(get("bases_trimmed_read_structure"), 24);
}
#[test]
fn execute_read_structure_too_short_drops_read_as_length_filter() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["ACG", "ACGTACGTAC"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let metrics = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], Some(metrics.clone()));
cmd.read_structures = vec![rs("5M+T")];
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 1, "expected only the 10-base read to survive");
assert_eq!(written[0].seq.as_slice(), b"CGTAC");
let contents = std::fs::read_to_string(&metrics).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let idx = |name: &str| header.iter().position(|h| *h == name).unwrap();
assert_eq!(values[idx("reads_in")], "2");
assert_eq!(values[idx("reads_out")], "1");
assert_eq!(values[idx("reads_filtered_length")], "1");
}
#[test]
fn execute_length_filter_drops_short_reads() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["ACGTACGTACGT", "AAA", "CCCCCCCCCC"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let metrics = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], Some(metrics.clone()));
cmd.filter_length = LengthFilter { min: 5, max: None };
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 2);
assert_eq!(written[0].seq.as_slice(), b"ACGTACGTACGT");
assert_eq!(written[1].seq.as_slice(), b"CCCCCCCCCC");
let contents = std::fs::read_to_string(&metrics).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let idx = |name: &str| header.iter().position(|h| *h == name).unwrap();
assert_eq!(values[idx("reads_in")], "3");
assert_eq!(values[idx("reads_out")], "2");
assert_eq!(values[idx("reads_filtered_length")], "1");
}
#[test]
fn execute_pe_one_mate_short_drops_both() {
let tmp = TempDir::new().unwrap();
let r1_lines = vec![
"@pair1".to_string(),
"AAAAAAAAAA".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
"@pair2".to_string(),
"CCCCCCCCCC".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
];
let r2_lines = vec![
"@pair1".to_string(),
"TTT".to_string(),
"+".to_string(),
"III".to_string(),
"@pair2".to_string(),
"GGGGGGGGGG".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
];
let r1 = write_fastq(&tmp, "r1", &r1_lines);
let r2 = write_fastq(&tmp, "r2", &r2_lines);
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let mut cmd = trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None);
cmd.filter_length = LengthFilter { min: 5, max: None };
cmd.execute().unwrap();
let w1 = read_fastq(&o1);
let w2 = read_fastq(&o2);
assert_eq!(w1.len(), 1);
assert_eq!(w2.len(), 1);
assert_eq!(w1[0].seq.as_slice(), b"CCCCCCCCCC");
assert_eq!(w2[0].seq.as_slice(), b"GGGGGGGGGG");
}
#[test]
fn execute_quality_trim_sliding_window() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@read1".to_string(),
"AAAAAAAAAACCCCCCCCCC".to_string(),
"+".to_string(),
"IIIIIIIIII&&&&&&&&&&".to_string(), ];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.quality_trim_3p = Some(QualityTrim { window: 4, threshold: 20 });
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written[0].seq.as_slice(), b"AAAAAAAAA");
}
#[test]
fn execute_mean_qual_filter_drops_low_quality_reads() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@good".to_string(),
"AAAAAAAAAA".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
"@bad".to_string(),
"CCCCCCCCCC".to_string(),
"+".to_string(),
"!!!!!!!!!!".to_string(), ];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.filter_mean_qual = Some(20);
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 1);
assert_eq!(written[0].seq.as_slice(), b"AAAAAAAAAA");
}
#[test]
fn execute_base_accounting_identity_holds() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@A:1:B:1:1:1:1".to_string(),
"AAAAAAAAAAACGTACGTAGATCGGAAGAG".to_string(),
"+".to_string(),
"I".repeat(30),
"@A:1:B:1:1:1:2".to_string(),
"CCCCCCCC".to_string(),
"+".to_string(),
"IIIIIIII".to_string(),
];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let metrics = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], Some(metrics.clone()));
cmd.read_structures = vec![rs("5S+T")];
cmd.adapter_sequence = vec!["AGATCGGAAGAG".to_string()];
cmd.filter_length = LengthFilter { min: 5, max: None };
cmd.execute().unwrap();
let contents = std::fs::read_to_string(&metrics).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let get = |name: &str| -> u64 {
let idx = header.iter().position(|h| *h == name).unwrap();
values[idx].parse().unwrap()
};
let bases_in = get("bases_in");
let bases_out = get("bases_out");
let sum_trimmed = get("bases_trimmed_read_structure")
+ get("bases_trimmed_adapter")
+ get("bases_trimmed_polyg")
+ get("bases_trimmed_polyx")
+ get("bases_trimmed_quality");
let filtered = get("bases_filtered");
assert_eq!(bases_in, bases_out + sum_trimmed + filtered);
}
#[test]
fn execute_polyx_picks_longest_single_base_tail() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["NNNNAAATTT"]));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.trim_polyx = Some(3);
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written[0].seq.as_slice(), b"NNNNAAA");
}
#[test]
fn execute_n_base_filter() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["ACGTACGT", "NNNNNNNN"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.filter_max_ns = Some(2);
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 1);
assert_eq!(written[0].seq.as_slice(), b"ACGTACGT");
}
#[test]
fn execute_writes_json_report_with_fastp_schema_keys() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@NB501234:1:FC:1:1:1:1".to_string(),
"AAAAAAAAAA".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
"@NB501234:1:FC:1:1:1:2".to_string(),
"CCCCCCCCCC".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(),
];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let json: serde_json::Value =
serde_json::from_str(&std::fs::read_to_string(&json_path).unwrap()).unwrap();
assert!(json.get("summary").is_some());
assert!(json.get("filtering_result").is_some());
assert!(json.get("adapter_cutting").is_some());
assert!(json.get("read1_before_filtering").is_some());
assert!(json.get("read1_after_filtering").is_some());
assert!(json.get("read2_before_filtering").is_none());
assert!(json.get("read2_after_filtering").is_none());
let summary = &json["summary"];
assert_eq!(summary["sequencing"], "single end");
assert_eq!(summary["before_filtering"]["total_reads"], 2);
assert_eq!(summary["before_filtering"]["total_bases"], 20);
assert_eq!(summary["before_filtering"]["q20_rate"], 1.0);
assert_eq!(summary["before_filtering"]["q30_rate"], 1.0);
}
#[test]
fn execute_pe_json_includes_read2() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["AAAAAAAA"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("r", &["TTTTTTTT"]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1, r2], vec![o1, o2], None);
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let json: serde_json::Value =
serde_json::from_str(&std::fs::read_to_string(&json_path).unwrap()).unwrap();
assert_eq!(json["summary"]["sequencing"], "paired end");
assert_eq!(json["read1_before_filtering"]["total_bases"], 8);
assert_eq!(json["read2_before_filtering"]["total_bases"], 8);
assert_eq!(json["summary"]["before_filtering"]["total_bases"], 16);
}
#[test]
fn execute_zero_reads_produces_finite_json() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &[]);
let out = tmp.path().join("out.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let raw = std::fs::read_to_string(&json_path).unwrap();
assert!(!raw.contains("NaN"));
let json: serde_json::Value = serde_json::from_str(&raw).unwrap();
assert_eq!(json["summary"]["before_filtering"]["total_reads"], 0);
assert_eq!(json["summary"]["before_filtering"]["q20_rate"], 0.0);
}
#[test]
fn execute_adapter_trimming_populates_adapter_cutting() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAAAAAAAAAGATCGGAAGAG", "CCCCCCCCCC"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.adapter_sequence = vec!["AGATCGGAAGAG".to_string()];
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let json: serde_json::Value =
serde_json::from_str(&std::fs::read_to_string(&json_path).unwrap()).unwrap();
assert_eq!(json["adapter_cutting"]["adapter_trimmed_reads"], 1);
assert_eq!(json["adapter_cutting"]["adapter_trimmed_bases"], 12);
}
#[test]
fn execute_json_rates_are_fractions_not_percentages() {
let tmp = TempDir::new().unwrap();
let lines =
vec!["@read".to_string(), "AAAA".to_string(), "+".to_string(), "I5I5".to_string()];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("out.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let json: serde_json::Value =
serde_json::from_str(&std::fs::read_to_string(&json_path).unwrap()).unwrap();
let q20 = json["summary"]["before_filtering"]["q20_rate"].as_f64().unwrap();
let q30 = json["summary"]["before_filtering"]["q30_rate"].as_f64().unwrap();
assert!((0.0..=1.0).contains(&q20), "q20_rate {q20} out of range");
assert!((0.0..=1.0).contains(&q30), "q30_rate {q30} out of range");
assert!((q20 - 1.0).abs() < 1e-9, "q20_rate should be 1.0, got {q20}");
assert!((q30 - 0.5).abs() < 1e-9, "q30_rate should be 0.5, got {q30}");
}
#[test]
fn execute_json_has_fastp_version_and_valid_json() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let raw = std::fs::read_to_string(&json_path).unwrap();
assert!(!raw.contains("NaN"));
assert!(!raw.contains("Infinity"));
let json: serde_json::Value = serde_json::from_str(&raw).unwrap();
let version = json["summary"]["fastp_version"].as_str().unwrap();
assert!(!version.is_empty(), "fastp_version should be populated");
}
#[test]
fn execute_pe_adapter_trimmed_counts_per_mate() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["AAAAAAAAAAAGATCGGAAGAG"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("r", &["CCCCCCCCCCAGATCGGAAGAG"]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1, r2], vec![o1, o2], None);
cmd.adapter_sequence = vec!["AGATCGGAAGAG".to_string(), "AGATCGGAAGAG".to_string()];
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let json: serde_json::Value =
serde_json::from_str(&std::fs::read_to_string(&json_path).unwrap()).unwrap();
assert_eq!(json["adapter_cutting"]["adapter_trimmed_reads"], 2);
}
#[test]
fn execute_filter_populates_filtering_result() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["ACGTACGTACGT", "AAA", "CCCCCCCCCC"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let json_path = tmp.path().join("trim.json");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.filter_length = LengthFilter { min: 5, max: None };
cmd.json = Some(json_path.clone());
cmd.execute().unwrap();
let json: serde_json::Value =
serde_json::from_str(&std::fs::read_to_string(&json_path).unwrap()).unwrap();
assert_eq!(json["filtering_result"]["passed_filter_reads"], 2);
assert_eq!(json["filtering_result"]["too_short_reads"], 1);
assert_eq!(json["filtering_result"]["low_quality_reads"], 0);
assert!(
json["filtering_result"].get("too_many_N_reads").is_some(),
"too_many_N_reads (uppercase N) is required for MultiQC fastp-module compat"
);
assert!(
json["filtering_result"].get("too_many_n_reads").is_none(),
"lowercase too_many_n_reads would be silently dropped by MultiQC"
);
}
#[test]
fn flatten_mate_stats_paired() {
let before = vec![
MateStats { reads: 10, bases: 100, q20_bases: 90, q30_bases: 80, gc_bases: 40 },
MateStats { reads: 10, bases: 100, q20_bases: 85, q30_bases: 70, gc_bases: 42 },
];
let after = vec![
MateStats { reads: 9, bases: 90, q20_bases: 88, q30_bases: 80, gc_bases: 36 },
MateStats { reads: 9, bases: 90, q20_bases: 82, q30_bases: 70, gc_bases: 38 },
];
let mut m = TrimMetrics::default();
flatten_mate_stats(&before, &after, &mut m);
assert_eq!(m.q20_before_r1, 90);
assert_eq!(m.q30_before_r1, 80);
assert_eq!(m.total_bases_before_r1, 100);
assert_eq!(m.q20_after_r1, 88);
assert_eq!(m.q30_after_r1, 80);
assert_eq!(m.q20_before_r2, 85);
assert_eq!(m.q30_before_r2, 70);
assert_eq!(m.q20_after_r2, 82);
assert_eq!(m.total_bases_after_r2, 90);
}
#[test]
fn flatten_mate_stats_single_end_leaves_r2_zero() {
let before =
vec![MateStats { reads: 5, bases: 50, q20_bases: 50, q30_bases: 45, gc_bases: 22 }];
let after =
vec![MateStats { reads: 4, bases: 40, q20_bases: 40, q30_bases: 35, gc_bases: 18 }];
let mut m = TrimMetrics::default();
flatten_mate_stats(&before, &after, &mut m);
assert_eq!(m.q20_before_r1, 50);
assert_eq!(m.q20_before_r2, 0);
assert_eq!(m.total_bases_before_r2, 0);
}
#[test]
fn ratio_zero_denom_returns_zero() {
assert_eq!(ratio(5, 0), 0.0);
}
#[test]
fn pct_zero_denom_returns_zero() {
assert_eq!(pct(5, 0), 0.0);
}
#[test]
fn low_qual_filter_parses_q_f() {
let f: LowQualFilter = "15:0.4".parse().unwrap();
assert_eq!(f.threshold, 15);
assert!((f.max_fraction - 0.4).abs() < 1e-9);
}
#[test]
fn low_qual_filter_rejects_fraction_out_of_range() {
assert!("15:1.5".parse::<LowQualFilter>().unwrap_err().contains("0.0..=1.0"));
assert!("15:-0.1".parse::<LowQualFilter>().unwrap_err().contains("0.0..=1.0"));
}
#[test]
fn low_qual_filter_rejects_missing_colon() {
assert!("15".parse::<LowQualFilter>().unwrap_err().contains("Q:F"));
}
#[test]
fn count_bases_below_q_simple() {
assert_eq!(count_bases_below_q(b"!!!!", 1), 4);
assert_eq!(count_bases_below_q(b"IIII", 40), 0);
assert_eq!(count_bases_below_q(b"IIII", 41), 4);
assert_eq!(count_bases_below_q(b"5555", 20), 0);
assert_eq!(count_bases_below_q(b"5555", 21), 4);
}
#[test]
fn count_bases_below_q_mixed_crosses_simd_boundary() {
let mut qual = vec![b'!'; 20];
qual.extend(std::iter::repeat_n(b'I', 20));
assert_eq!(count_bases_below_q(&qual, 10), 20);
}
#[test]
fn quality_trim_5p_truncates_at_first_bad_window() {
let mut rec =
OwnedRecord { head: vec![], seq: b"ACGTACGTAC".to_vec(), qual: b"IIIIII!!!!".to_vec() };
let trimmed = trim_quality_sliding_5prime(&mut rec, 4, 20);
assert_eq!(trimmed, 5);
assert_eq!(rec.seq.as_slice(), b"ACGTA");
assert_eq!(rec.qual.as_slice(), b"IIIII");
}
#[test]
fn quality_trim_5p_empties_all_bad_read() {
let mut rec = OwnedRecord { head: vec![], seq: b"ACGT".to_vec(), qual: b"!!!!".to_vec() };
let trimmed = trim_quality_sliding_5prime(&mut rec, 2, 20);
assert_eq!(trimmed, 4);
assert!(rec.seq.is_empty());
assert!(rec.qual.is_empty());
}
#[test]
fn quality_trim_5p_noop_when_all_good() {
let mut rec =
OwnedRecord { head: vec![], seq: b"ACGTACGT".to_vec(), qual: b"IIIIIIII".to_vec() };
let trimmed = trim_quality_sliding_5prime(&mut rec, 4, 20);
assert_eq!(trimmed, 0);
assert_eq!(rec.seq.as_slice(), b"ACGTACGT");
assert_eq!(rec.qual.as_slice(), b"IIIIIIII");
}
#[test]
fn quality_trim_5p_noop_when_shorter_than_window() {
let mut rec = OwnedRecord { head: vec![], seq: b"ACG".to_vec(), qual: b"!!!".to_vec() };
let trimmed = trim_quality_sliding_5prime(&mut rec, 4, 20);
assert_eq!(trimmed, 0);
assert_eq!(rec.seq.as_slice(), b"ACG");
}
#[test]
fn execute_filter_low_qual_drops_reads_with_too_many_lowq_bases() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@good".to_string(),
"AAAAAAAAAA".to_string(),
"+".to_string(),
"IIIIIIIIII".to_string(), "@bad".to_string(),
"CCCCCCCCCC".to_string(),
"+".to_string(),
"!!!!!!!!!!".to_string(), ];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("o.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.filter_low_qual = Some(LowQualFilter { threshold: 15, max_fraction: 0.4 });
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 1);
assert_eq!(written[0].seq.as_slice(), b"AAAAAAAAAA");
}
#[test]
fn polyg_zero_disables_trim() {
let tmp = TempDir::new().unwrap();
let lines =
vec!["@r".to_string(), "ACGTGGGGGGGGGGGG".to_string(), "+".to_string(), "I".repeat(16)];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("o.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.trim_polyg = 0;
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written[0].seq.as_slice(), b"ACGTGGGGGGGGGGGG");
}
#[test]
fn panic_message_downcasts_static_str() {
let payload = std::panic::catch_unwind(|| panic!("static str panic")).unwrap_err();
assert_eq!(panic_message(&payload), "static str panic");
}
#[test]
fn panic_message_downcasts_string() {
let payload = std::panic::catch_unwind(|| panic!("string panic: {}", 42)).unwrap_err();
assert_eq!(panic_message(&payload), "string panic: 42");
}
#[test]
fn panic_message_falls_back_on_unknown_payload() {
#[derive(Debug)]
struct Weird;
let payload = std::panic::catch_unwind(|| std::panic::panic_any(Weird)).unwrap_err();
assert_eq!(panic_message(&payload), "(non-string panic payload)");
}
#[test]
fn select_most_specific_error_empty_returns_none() {
assert!(select_most_specific_error(vec![]).is_none());
}
#[test]
fn select_most_specific_error_prefers_non_symptom() {
let errors = vec![
anyhow!("writer exited before receiving batch slot"),
anyhow!("disk full: ENOSPC"),
];
let e = select_most_specific_error(errors).unwrap();
assert!(e.to_string().contains("disk full"));
}
#[test]
fn select_most_specific_error_falls_back_to_first_when_all_symptoms() {
let errors = vec![
anyhow!("writer exited before receiving batch slot"),
anyhow!("worker dropped before handoff"),
];
let e = select_most_specific_error(errors).unwrap();
assert!(e.to_string().contains("writer exited"));
}
fn drive_expected_insert(stats: &mut OverlapStats, mean_l: usize) {
stats.sum_insert = 0;
stats.count_detect = 0;
stats.pairs_since_update = 0;
for _ in 0..INSERT_STATS_MIN_DETECTIONS {
stats.sum_insert += mean_l as u64;
stats.count_detect += 1;
}
stats.pairs_since_update = INSERT_STATS_UPDATE_INTERVAL;
stats.maybe_update_expected_insert();
}
#[test]
fn overlap_stats_expected_insert_does_not_flap_within_margin() {
let mut stats = OverlapStats::new(None);
drive_expected_insert(&mut stats, 300);
assert_eq!(stats.expected_insert, Some(300));
for &mean in &[305usize, 295, 310, 290, 313, 287] {
let before = stats.expected_insert;
drive_expected_insert(&mut stats, mean);
assert_eq!(stats.expected_insert, before, "expected_insert flapped at mean={mean}");
}
drive_expected_insert(&mut stats, 400);
assert_eq!(stats.expected_insert, Some(400));
}
#[test]
fn overlap_stats_no_update_below_min_detections() {
let mut stats = OverlapStats::new(None);
for _ in 0..(INSERT_STATS_MIN_DETECTIONS - 1) {
stats.sum_insert += 50;
stats.count_detect += 1;
}
stats.pairs_since_update = INSERT_STATS_UPDATE_INTERVAL;
stats.maybe_update_expected_insert();
assert_eq!(stats.expected_insert, None);
}
#[test]
fn expected_insert_size_hint_takes_effect_via_center_shift() {
let stats = OverlapStats::new(Some(250));
assert_eq!(stats.expected_insert, Some(250));
assert_eq!(stats.center_shift(150), 100);
assert_eq!(stats.center_shift(125), 125);
let no_hint = OverlapStats::new(None);
assert_eq!(no_hint.expected_insert, None);
assert_eq!(no_hint.center_shift(150), isize::MIN);
}
#[test]
fn center_shift_uses_r2_length_for_asymmetric_pairs() {
let stats = OverlapStats::new(Some(300));
assert_eq!(stats.center_shift(100), 200);
assert_eq!(stats.center_shift(200), 100);
}
#[test]
fn overlap_stats_margin_scales_with_insert_estimate() {
let mut stats = OverlapStats::new(None);
drive_expected_insert(&mut stats, 1000);
assert_eq!(stats.expected_insert, Some(1000));
drive_expected_insert(&mut stats, 1049);
assert_eq!(stats.expected_insert, Some(1000), "moved at +49 (inside margin)");
drive_expected_insert(&mut stats, 1050);
assert_eq!(stats.expected_insert, Some(1050), "did not move at +50 (at margin)");
}
#[test]
fn overlap_stats_learned_mean_can_replace_hint() {
let mut stats = OverlapStats::new(Some(300));
assert_eq!(stats.expected_insert, Some(300));
drive_expected_insert(&mut stats, 350);
assert_eq!(stats.expected_insert, Some(350));
}
#[test]
fn execute_errors_when_r1_and_r2_have_different_record_counts() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("p", &["ACGT", "ACGT", "ACGT"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("p", &["ACGT", "ACGT"]));
let o1 = tmp.path().join("o1.fq.gz");
let o2 = tmp.path().join("o2.fq.gz");
let cmd = trim_cmd(vec![r1, r2], vec![o1, o2], None);
let err = cmd.execute().unwrap_err().to_string();
assert!(err.contains("out of sync"), "expected out-of-sync error, got: {err}");
}
#[test]
fn execute_filter_drops_all_reads_accounting_holds() {
let tmp = TempDir::new().unwrap();
let lines = vec![
"@a".to_string(),
"AAAA".to_string(),
"+".to_string(),
"!!!!".to_string(),
"@b".to_string(),
"CCCC".to_string(),
"+".to_string(),
"!!!!".to_string(),
];
let r1 = write_fastq(&tmp, "r1", &lines);
let out = tmp.path().join("o.fq.gz");
let metrics = tmp.path().join("m.txt");
let mut cmd = trim_cmd(vec![r1], vec![out], Some(metrics.clone()));
cmd.filter_mean_qual = Some(20);
cmd.execute().unwrap();
let contents = std::fs::read_to_string(&metrics).unwrap();
let lines: Vec<&str> = contents.lines().collect();
let header: Vec<&str> = lines[0].split('\t').collect();
let values: Vec<&str> = lines[1].split('\t').collect();
let get = |name: &str| -> u64 {
let idx = header.iter().position(|h| *h == name).unwrap();
values[idx].parse().unwrap()
};
assert_eq!(get("reads_in"), 2);
assert_eq!(get("reads_out"), 0);
assert_eq!(get("reads_filtered_quality"), 2);
assert_eq!(get("bases_out"), 0);
let sum_trimmed = get("bases_trimmed_read_structure")
+ get("bases_trimmed_adapter")
+ get("bases_trimmed_polyg")
+ get("bases_trimmed_polyx")
+ get("bases_trimmed_quality");
assert_eq!(get("bases_in"), get("bases_out") + sum_trimmed + get("bases_filtered"));
}
#[test]
fn execute_iupac_adapter_trims_via_scalar_fallback() {
let tmp = TempDir::new().unwrap();
let reads: Vec<&str> = vec!["AAAAAAAAAAAGATCGGAAGAG"];
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &reads));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.adapter_sequence = vec!["AGATCGGNAGAG".to_string()];
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 1);
assert_eq!(written[0].seq.as_slice(), b"AAAAAAAAAA");
}
#[test]
fn interleaved_input_detected_matches_split_input_output() {
let tmp = TempDir::new().unwrap();
let text = interleaved_fq_text(4, "ACGTACGTAC", "TGCATGCATG");
let interleaved = write_bytes(&tmp, "interleaved.fq", text.as_bytes());
let out1a = tmp.path().join("split_out1.fq");
let out2a = tmp.path().join("split_out2.fq");
trim_cmd(vec![interleaved], vec![out1a.clone(), out2a.clone()], None).execute().unwrap();
let mut r1_text = String::new();
let mut r2_text = String::new();
for i in 0..4 {
r1_text += &fq_record(&format!("pair{i}/1"), "ACGTACGTAC");
r2_text += &fq_record(&format!("pair{i}/2"), "TGCATGCATG");
}
let r1 = write_bytes(&tmp, "r1.fq", r1_text.as_bytes());
let r2 = write_bytes(&tmp, "r2.fq", r2_text.as_bytes());
let out1b = tmp.path().join("real_split_out1.fq");
let out2b = tmp.path().join("real_split_out2.fq");
trim_cmd(vec![r1, r2], vec![out1b.clone(), out2b.clone()], None).execute().unwrap();
assert_eq!(read_fastq(&out1a), read_fastq(&out1b));
assert_eq!(read_fastq(&out2a), read_fastq(&out2b));
}
#[test]
fn interleaved_input_to_interleaved_output_roundtrips() {
let tmp = TempDir::new().unwrap();
let text = interleaved_fq_text(3, "AAAACCCCGG", "TTTTGGGGCC");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
trim_cmd(vec![interleaved], vec![out.clone()], None).execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 6);
for i in 0..3 {
assert_eq!(written[2 * i].head, format!("pair{i}/1").as_bytes());
assert_eq!(written[2 * i + 1].head, format!("pair{i}/2").as_bytes());
}
}
#[test]
fn interleaved_input_spanning_multiple_batches_preserves_pairing() {
let tmp = TempDir::new().unwrap();
let text = interleaved_fq_text(10, "AAAACCCCGG", "TTTTGGGGCC");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
let mut cmd = trim_cmd(vec![interleaved], vec![out.clone()], None);
cmd.batch_size = 2;
cmd.execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 20);
for i in 0..10 {
assert_eq!(written[2 * i].head, format!("pair{i}/1").as_bytes());
assert_eq!(written[2 * i + 1].head, format!("pair{i}/2").as_bytes());
}
}
#[test]
fn split_input_to_interleaved_output_orders_r1_then_r2_per_pair() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("pair", &["AAAA", "CCCC"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("pair", &["GGGG", "TTTT"]));
let out = tmp.path().join("interleaved_out.fq");
trim_cmd(vec![r1, r2], vec![out.clone()], None).execute().unwrap();
let written = read_fastq(&out);
assert_eq!(written.len(), 4);
assert_eq!(written[0].seq.as_slice(), b"AAAA");
assert_eq!(written[1].seq.as_slice(), b"GGGG");
assert_eq!(written[2].seq.as_slice(), b"CCCC");
assert_eq!(written[3].seq.as_slice(), b"TTTT");
}
#[test]
fn interleaved_input_to_split_outputs() {
let tmp = TempDir::new().unwrap();
let text = interleaved_fq_text(3, "AAAAAAAAAA", "CCCCCCCCCC");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out1 = tmp.path().join("out1.fq");
let out2 = tmp.path().join("out2.fq");
trim_cmd(vec![interleaved], vec![out1.clone(), out2.clone()], None).execute().unwrap();
let r1 = read_fastq(&out1);
let r2 = read_fastq(&out2);
assert_eq!(r1.len(), 3);
assert_eq!(r2.len(), 3);
assert!(r1.iter().all(|r| r.seq == b"AAAAAAAAAA"));
assert!(r2.iter().all(|r| r.seq == b"CCCCCCCCCC"));
}
#[test]
fn single_input_with_distinct_names_treated_as_single_end() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "se", &fq_lines("distinct", &["ACGT", "TGCA", "AAAA"]));
let out = tmp.path().join("out.fq");
trim_cmd(vec![r1], vec![out.clone()], None).execute().unwrap();
assert_eq!(read_fastq(&out).len(), 3);
}
#[test]
fn slash_suffixed_mate_names_detected_as_interleaved() {
let tmp = TempDir::new().unwrap();
let text = interleaved_fq_text(2, "AAAAAAAAAA", "CCCCCCCCCC");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out1 = tmp.path().join("out1.fq");
let out2 = tmp.path().join("out2.fq");
trim_cmd(vec![interleaved], vec![out1.clone(), out2.clone()], None).execute().unwrap();
assert_eq!(read_fastq(&out1).len(), 2);
assert_eq!(read_fastq(&out2).len(), 2);
}
#[test]
fn interleaved_input_name_mismatch_mid_stream_errors() {
let tmp = TempDir::new().unwrap();
let mut text = String::new();
text += &fq_record("pair0/1", "ACGTACGTAC");
text += &fq_record("pair0/2", "TGCATGCATG");
text += &fq_record("pair1/1", "ACGTACGTAC");
text += &fq_record("pair1/2", "TGCATGCATG");
text += &fq_record("pair2/1", "ACGTACGTAC");
text += &fq_record("totally_different", "TGCATGCATG");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
let err = trim_cmd(vec![interleaved], vec![out], None).execute().unwrap_err().to_string();
assert!(err.contains("out of sync"), "{err}");
assert!(err.contains("pair 3"), "expected the error to mention the pair index: {err}");
assert!(
err.contains("file records 5/6"),
"expected the error to mention the file record indices: {err}"
);
}
#[test]
fn interleaved_input_odd_record_count_errors() {
let tmp = TempDir::new().unwrap();
let mut text = String::new();
text += &fq_record("pair0/1", "ACGTACGTAC");
text += &fq_record("pair0/2", "TGCATGCATG");
text += &fq_record("pair1/1", "ACGTACGTAC"); let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
let err = trim_cmd(vec![interleaved], vec![out], None).execute().unwrap_err().to_string();
assert!(err.contains("odd record count"), "{err}");
}
#[test]
fn interleaved_input_reversed_pair_mid_stream_errors() {
let tmp = TempDir::new().unwrap();
let mut text = String::new();
text += &fq_record("pair0/1", "ACGTACGTAC");
text += &fq_record("pair0/2", "TGCATGCATG");
text += &fq_record("pair1/1", "ACGTACGTAC");
text += &fq_record("pair1/2", "TGCATGCATG");
text += &fq_record("pair2/2", "TGCATGCATG"); text += &fq_record("pair2/1", "ACGTACGTAC");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
let err = trim_cmd(vec![interleaved], vec![out], None).execute().unwrap_err().to_string();
assert!(err.contains("out of sync"), "{err}");
}
#[test]
fn split_pe_unrecognized_naming_pairs_by_position() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("foo", &["AAAA", "CCCC"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("bar", &["GGGG", "TTTT"]));
let o1 = tmp.path().join("o1.fq");
let o2 = tmp.path().join("o2.fq");
trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None).execute().unwrap();
assert_eq!(read_fastq(&o1).len(), 2);
assert_eq!(read_fastq(&o2).len(), 2);
}
#[test]
fn split_pe_unrecognized_naming_still_errors_on_record_count_mismatch() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("foo", &["ACGT", "ACGT", "ACGT"]));
let r2 = write_fastq(&tmp, "r2", &fq_lines("bar", &["ACGT", "ACGT"]));
let o1 = tmp.path().join("o1.fq");
let o2 = tmp.path().join("o2.fq");
let err = trim_cmd(vec![r1, r2], vec![o1, o2], None).execute().unwrap_err().to_string();
assert!(err.contains("out of sync"), "{err}");
}
#[test]
fn interleaved_input_in_mate_2_first_order_errors() {
let tmp = TempDir::new().unwrap();
let mut text = String::new();
text += &fq_record("pair0/2", "ACGTACGTAC");
text += &fq_record("pair0/1", "TGCATGCATG");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
let err = trim_cmd(vec![interleaved], vec![out.clone()], None).execute().unwrap_err();
assert!(err.to_string().contains("mate-2/mate-1 order"), "{err}");
assert!(!out.exists());
}
#[test]
fn split_pe_swapped_inputs_error() {
let tmp = TempDir::new().unwrap();
let r1 = write_bytes(&tmp, "r1.fq", fq_record("pair0/1", "ACGT").as_bytes());
let r2 = write_bytes(&tmp, "r2.fq", fq_record("pair0/2", "ACGT").as_bytes());
let o1 = tmp.path().join("o1.fq");
let o2 = tmp.path().join("o2.fq");
let err = trim_cmd(vec![r2, r1], vec![o1, o2], None).execute().unwrap_err().to_string();
assert!(err.contains("wrong order"), "{err}");
}
#[test]
fn split_pe_mid_stream_name_mismatch_errors() {
let tmp = TempDir::new().unwrap();
let mut r1_text = String::new();
r1_text += &fq_record("pair0/1", "ACGT");
r1_text += &fq_record("pair1/1", "ACGT");
let mut r2_text = String::new();
r2_text += &fq_record("pair0/2", "ACGT");
r2_text += &fq_record("mismatched", "ACGT");
let r1 = write_bytes(&tmp, "r1.fq", r1_text.as_bytes());
let r2 = write_bytes(&tmp, "r2.fq", r2_text.as_bytes());
let o1 = tmp.path().join("o1.fq");
let o2 = tmp.path().join("o2.fq");
let err = trim_cmd(vec![r1, r2], vec![o1, o2], None).execute().unwrap_err().to_string();
assert!(err.contains("R1/R2 read names do not correspond"), "{err}");
}
#[test]
fn split_pe_sra_default_defline_names_accepted() {
let tmp = TempDir::new().unwrap();
let mut r1_text = String::new();
let mut r2_text = String::new();
for i in 1..=3 {
r1_text += &fq_record(&format!("SRR390728.{i} {i} length=10"), "ACGTACGTAC");
r2_text += &fq_record(&format!("SRR390728.{i} {i} length=10"), "TGCATGCATG");
}
let r1 = write_bytes(&tmp, "r1.fq", r1_text.as_bytes());
let r2 = write_bytes(&tmp, "r2.fq", r2_text.as_bytes());
let o1 = tmp.path().join("o1.fq");
let o2 = tmp.path().join("o2.fq");
trim_cmd(vec![r1, r2], vec![o1.clone(), o2.clone()], None).execute().unwrap();
assert_eq!(read_fastq(&o1).len(), 3);
assert_eq!(read_fastq(&o2).len(), 3);
}
#[test]
fn interleaved_sra_default_defline_names_detected_as_pe() {
let tmp = TempDir::new().unwrap();
let mut text = String::new();
for i in 1..=3 {
text += &fq_record(&format!("SRR390728.{i} {i} length=10"), "ACGTACGTAC");
text += &fq_record(&format!("SRR390728.{i} {i} length=10"), "TGCATGCATG");
}
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let o1 = tmp.path().join("o1.fq");
let o2 = tmp.path().join("o2.fq");
trim_cmd(vec![interleaved], vec![o1.clone(), o2.clone()], None).execute().unwrap();
let (mate1, mate2) = (read_fastq(&o1), read_fastq(&o2));
assert_eq!(mate1.len(), 3);
assert_eq!(mate2.len(), 3);
assert!(mate1.iter().all(|r| r.seq == b"ACGTACGTAC"));
assert!(mate2.iter().all(|r| r.seq == b"TGCATGCATG"));
}
#[test]
fn empty_single_input_treated_as_single_end() {
let tmp = TempDir::new().unwrap();
let empty = write_bytes(&tmp, "empty.fq", b"");
let out = tmp.path().join("out.fq");
trim_cmd(vec![empty], vec![out.clone()], None).execute().unwrap();
assert!(read_fastq(&out).is_empty());
}
#[test]
fn empty_single_input_with_two_outputs_yields_two_valid_empty_outputs() {
let tmp = TempDir::new().unwrap();
let empty = write_bytes(&tmp, "empty.fq", b"");
let out1 = tmp.path().join("out1.fq");
let out2 = tmp.path().join("out2.fq");
trim_cmd(vec![empty], vec![out1.clone(), out2.clone()], None).execute().unwrap();
assert!(out1.exists());
assert!(out2.exists());
assert!(read_fastq(&out1).is_empty());
assert!(read_fastq(&out2).is_empty());
}
#[test]
fn empty_single_input_with_two_read_structures_succeeds() {
let tmp = TempDir::new().unwrap();
let empty = write_bytes(&tmp, "empty.fq", b"");
let out = tmp.path().join("out.fq");
let mut cmd = trim_cmd(vec![empty], vec![out.clone()], None);
cmd.read_structures = vec![rs("+T"), rs("+T")];
cmd.execute().unwrap();
assert!(read_fastq(&out).is_empty());
}
#[test]
fn plain_extension_output_is_uncompressed() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq"); trim_cmd(vec![r1], vec![out.clone()], None).execute().unwrap();
let bytes = std::fs::read(&out).unwrap();
assert!(!bytes.starts_with(&[0x1f, 0x8b]), "expected plain text, got gzip-magic bytes");
assert!(bytes.starts_with(b"@r_0"));
}
#[test]
fn gz_extension_output_is_bgzf() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq.gz");
trim_cmd(vec![r1], vec![out.clone()], None).execute().unwrap();
let bytes = std::fs::read(&out).unwrap();
assert!(bytes.starts_with(&[0x1f, 0x8b]), "expected gzip/BGZF magic bytes");
}
#[test]
fn bgz_extension_output_is_bgzf() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq.bgz");
trim_cmd(vec![r1], vec![out.clone()], None).execute().unwrap();
let bytes = std::fs::read(&out).unwrap();
assert!(bytes.starts_with(&[0x1f, 0x8b]), "expected gzip/BGZF magic bytes for .bgz");
}
#[test]
fn uppercase_gz_extension_output_is_bgzf() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq.GZ");
trim_cmd(vec![r1], vec![out.clone()], None).execute().unwrap();
let bytes = std::fs::read(&out).unwrap();
assert!(bytes.starts_with(&[0x1f, 0x8b]), "expected gzip/BGZF magic bytes for .GZ");
}
#[test]
fn output_compression_override_forces_bgzf() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq"); let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.output_compression = OutputCompression::Bgzf;
cmd.execute().unwrap();
let bytes = std::fs::read(&out).unwrap();
assert!(bytes.starts_with(&[0x1f, 0x8b]), "expected gzip/BGZF magic bytes despite no .gz");
}
#[test]
fn output_compression_override_forces_none() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "r1", &fq_lines("r", &["ACGT"]));
let out = tmp.path().join("out.fq.gz");
let mut cmd = trim_cmd(vec![r1], vec![out.clone()], None);
cmd.output_compression = OutputCompression::None;
cmd.execute().unwrap();
let bytes = std::fs::read(&out).unwrap();
assert!(!bytes.starts_with(&[0x1f, 0x8b]), "expected plain text despite .gz name");
assert!(bytes.starts_with(b"@r_0"));
}
#[test]
fn gzip_content_with_plain_name_is_read() {
let tmp = TempDir::new().unwrap();
let text = fq_lines("r", &["ACGT", "TTTT"]).join("\n") + "\n";
let gz_bytes = gzip_bytes(text.as_bytes());
let r1 = write_bytes(&tmp, "input.txt", &gz_bytes);
let out = tmp.path().join("out.fq");
trim_cmd(vec![r1], vec![out.clone()], None).execute().unwrap();
assert_eq!(read_fastq(&out).len(), 2);
}
#[test]
fn interleaved_input_with_two_read_structures_succeeds() {
let tmp = TempDir::new().unwrap();
let text = interleaved_fq_text(2, "AAAACCCCTT", "GGGGTTTTAA");
let interleaved = write_bytes(&tmp, "in.fq", text.as_bytes());
let out = tmp.path().join("out.fq");
let mut cmd = trim_cmd(vec![interleaved], vec![out], None);
cmd.read_structures = vec![rs("+T"), rs("+T")];
cmd.execute().unwrap();
}
#[test]
fn sniffed_se_input_with_two_read_structures_errors() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "se", &fq_lines("distinct", &["ACGT"]));
let out = tmp.path().join("out.fq");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.read_structures = vec![rs("+T"), rs("+T")];
let err = cmd.execute().unwrap_err().to_string();
assert!(err.contains("number of mates"), "{err}");
}
#[test]
fn sniffed_se_input_with_two_adapter_sequences_errors() {
let tmp = TempDir::new().unwrap();
let r1 = write_fastq(&tmp, "se", &fq_lines("distinct", &["ACGTACGTACGT"]));
let out = tmp.path().join("out.fq");
let mut cmd = trim_cmd(vec![r1], vec![out], None);
cmd.adapter_sequence = vec!["AGATCGGAAG".to_string(), "AGATCGGAAG".to_string()];
let err = cmd.execute().unwrap_err().to_string();
assert!(err.contains("mate(s)"), "{err}");
}
}