From a586056571208cdda5bd18de8495ae86641390c4 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Thu, 10 Sep 2026 12:52:07 +0200 Subject: [PATCH] feat: add bqc sniff index for reusable Salmon indexes Persistent Salmon 2.x transcriptome index from a FASTA (plain or gzipped), written beside the FASTA as .index or to --output-directory/--prefix. Same sniff-strand feature as sniff strand. --- README.md | 2 +- assets/changelog/changelog.md | 12 ++ assets/docs/usage.md | 46 ++++-- bqc/Cargo.lock | 2 +- bqc/Cargo.toml | 11 +- bqc/README.md | 21 ++- bqc/src/cli.rs | 261 +++++++++++++++++++++++++++++++++- bqc/src/sniff.rs | 17 +-- bqc/src/sniff/strand.rs | 43 ++++-- bqc/tests/strand.rs | 147 ++++++++++++++++++- 10 files changed, 508 insertions(+), 54 deletions(-) diff --git a/README.md b/README.md index f2c6cbd..4cdaa16 100644 --- a/README.md +++ b/README.md @@ -83,7 +83,7 @@ > [!IMPORTANT] > - πŸ“„ **CBQ-native**: reads CBQ directly. > - ⚑~ 4-8x faster than the field on typical data (fastp, cutadapt, atropos, trimmomatic), with πŸ“‰ ~90% less memory. -> - 🧬 **Reference-free sniff adapters** (curated 234-known adapters) and πŸ‘ƒ **sniff strandeness** for RNA against a Salmon index (given or produced on-the-fly through `--transcriptome `). +> - 🧬 **Reference-free sniff adapters** (curated 234-known adapters) and πŸ‘ƒ **sniff strandeness** for RNA against a Salmon index (`bqc sniff index --sequence tx.fa`, or on-the-fly through `--transcriptome `). > - πŸ”„ **Everything in one pass**: 3' adapter removal (incl. `--allow-indels`), quality/positional/homopolymer trim, polyA/G tails, overlap error correction from the mate, internal-adapter split, UMI extraction, filtering β€” configurable. > - πŸ”’ Deterministic + safe: single static binary diff --git a/assets/changelog/changelog.md b/assets/changelog/changelog.md index 4bf10ad..8299ed8 100644 --- a/assets/changelog/changelog.md +++ b/assets/changelog/changelog.md @@ -71,6 +71,18 @@ All notable changes to this project are documented in this file. +## 0.0.4 β€” 2026-09-10 + +### Added + +- **`sniff index`** β€” builds a reusable Salmon 2.x transcriptome index from + a FASTA (plain or gzipped) for `sniff strand --index`. Default destination + is `.index` beside the FASTA; `--output-directory` and `--prefix` + override location and name. `--force` replaces an existing directory; + without it an existing dest is refused before the build. Same optional + `sniff-strand` feature as `sniff strand`. `--transcriptome` on `sniff + strand` remains the throwaway one-shot. + ## 0.0.3 β€” 2026-08-13 ### Added diff --git a/assets/docs/usage.md b/assets/docs/usage.md index 5459a55..ee2d624 100644 --- a/assets/docs/usage.md +++ b/assets/docs/usage.md @@ -86,6 +86,7 @@ bqc workflow reads.cbq -o out.cbq --config illumina.toml -T 8 bqc umi reads.cbq -o out.cbq --umi-location read1 --umi-length 8 bqc dedup reads.cbq -o unique.cbq -T 8 bqc sniff adapters reads.cbq # inspect, never modify +bqc sniff index --sequence transcripts.fa # reusable Salmon index bqc sniff strand reads.cbq --index salmon-index # RNA-seq orientation ``` @@ -98,8 +99,8 @@ dataset-level QC reports; `bqc` itself only reads and writes CBQ. cargo install bqc ``` -`bqc sniff strand` is optional and off by default β€” its mapping engine -pulls in a large dependency tree and needs a C compiler: +`bqc sniff strand` and `bqc sniff index` are optional and off by default β€” +the Salmon mapper pulls in a large dependency tree and needs a C compiler: ```bash cargo install bqc --features sniff-strand @@ -605,10 +606,9 @@ bqc sniff strand reads.cbq --transcriptome transcripts.fa # index built on the `--index` and `--transcriptome` are mutually exclusive; exactly one is required. **`--transcriptome`** takes a transcriptome FASTA instead of a ready index and builds a Salmon index from it on the fly in a temporary directory, -discarded when the run ends. There is no persistent index cache β€” the index is -rebuilt on every run β€” so reach for `--index` when the same reference is used -repeatedly. The file must exist and be a FASTA; a missing path is a -configuration error. +discarded when the run ends. For a reusable index, run `bqc sniff index` +once and pass the directory to `--index`. The file must exist and be a FASTA; +a missing path is a configuration error. Two related results are reported: @@ -666,10 +666,39 @@ The report also records index provenance β€” reference count, k-mer length, decoy state and content hashes β€” so a result can be matched back to the index that produced it. +### 5.11 `bqc sniff index` + +Builds a reusable Salmon 2.x transcriptome index from a FASTA (plain or +gzipped). Same feature as `sniff strand`. The output is a directory, suitable +as `bqc sniff strand --index`. + +```bash +bqc sniff index --sequence transcripts.fa +bqc sniff index --sequence transcripts.fa.gz --output-directory /idx --prefix hg38 +idx=$(bqc sniff index --sequence transcripts.fa) +bqc sniff strand reads.cbq --index "$idx" +``` + +```text +--sequence transcriptome FASTA, optionally gzipped +--output-directory parent directory of the index [FASTA's directory] +--prefix index directory name [.index] +-T, --threads worker threads; 0 uses every core +--force replace an existing index directory +``` + +Default name: strip a trailing `.gz`, then `.fa` / `.fasta` / `.fna` +(any case), and append `.index`. `--prefix` is that last path component as +given β€” it does not add `.index`. It is a name, not a path; use +`--output-directory` to choose the location. The destination path is printed +to stdout. Without `--force`, an existing directory is refused before the +index is built. + ## 6. Sniff reports -Both sniff commands accept `--format text` (default), `json` or `tsv`, and -`-o` writes the report atomically (stdout otherwise). +`sniff adapters` and `sniff strand` accept `--format text` (default), `json` +or `tsv`, and `-o` writes the report atomically (stdout otherwise). `sniff +index` prints the index directory path to stdout and has no report. * **JSON** is the stable pipeline interface. It carries a `schema_version` and a common envelope: `tool`, `input` (path, records, schema flags), @@ -818,6 +847,7 @@ cardinality. bqc sniff adapters sample.cbq --require-confident \ --emit-config sample.adapters.toml --format json -o sample.adapters.json # 2. Is the library stranded? (needs a Salmon index) +bqc sniff index --sequence transcripts.fa --prefix reference.salmon bqc sniff strand sample.cbq --index reference.salmon \ --require-confident --format json -o sample.strand.json # 3. Clean the data in one pass using the detected adapters. diff --git a/bqc/Cargo.lock b/bqc/Cargo.lock index 63465cb..0f0bc32 100644 --- a/bqc/Cargo.lock +++ b/bqc/Cargo.lock @@ -309,7 +309,7 @@ dependencies = [ [[package]] name = "bqc" -version = "0.0.3" +version = "0.0.4" dependencies = [ "binseq", "bitflags 2.13.1", diff --git a/bqc/Cargo.toml b/bqc/Cargo.toml index 5012a57..f1d1ebe 100644 --- a/bqc/Cargo.toml +++ b/bqc/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "bqc" -version = "0.0.3" +version = "0.0.4" edition = "2024" description = "CBQ-native all-in-one quality control tool" license = "MIT" @@ -13,10 +13,11 @@ rust-version = "1.91" [features] default = [] -# `bqc sniff strand`. Opt-in: the Salmon mapper pulls in roughly 330 crates -# and needs a C compiler, which is a large cost to impose on every build for one -# subcommand. This gates optional compilation, not a choice of backend β€” there -# is one implementation, and no subprocess or FASTQ fallback exists. +# `bqc sniff strand` and `bqc sniff index`. Opt-in: the Salmon mapper pulls in +# roughly 330 crates and needs a C compiler, which is a large cost to impose on +# every build for one pair of subcommands. This gates optional compilation, not +# a choice of backend β€” there is one implementation, and no subprocess or FASTQ +# fallback exists. sniff-strand = [ "dep:piscem-rs", "dep:salmon-core", diff --git a/bqc/README.md b/bqc/README.md index 596436f..f7a510a 100644 --- a/bqc/README.md +++ b/bqc/README.md @@ -14,6 +14,7 @@ bqc filter reads.cbq -o out.cbq --min-length 30 --failed rejected.cbq bqc correct pairs.cbq -o out.cbq --correction-log corrections.tsv bqc workflow reads.cbq -o out.cbq --config illumina.toml -T 8 bqc sniff adapters reads.cbq # inspect, never modify +bqc sniff index --sequence transcripts.fa # reusable Salmon index bqc sniff strand reads.cbq --index salmon-index # RNA-seq orientation ``` @@ -31,6 +32,7 @@ bqc sniff strand reads.cbq --index salmon-index # RNA-seq orientation | Split reads at internal adapters into separate records | `bqc segment` | | All of the above fused into one pass | `bqc workflow` | | Non-destructive inspection: which adapters are in this file? | `bqc sniff adapters` | +| Reusable Salmon transcriptome index from a FASTA | `bqc sniff index` | | RNA-seq library strandedness, against a transcriptome | `bqc sniff strand` | `bqtools qc` describes a dataset; `bqc filter` decides whether an individual @@ -303,9 +305,11 @@ three_prime = "GGATCCTAAGCC" ## Sniffing -`bqc sniff` inspects a file without changing it. It never trims, filters, -reorders or rewrites anything; the input is opened read-only and is -byte-identical afterwards. +`bqc sniff adapters` and `bqc sniff strand` inspect a file without changing +it. They never trim, filter, reorder or rewrite anything; the input is opened +read-only and is byte-identical afterwards. `bqc sniff index` builds a reusable +Salmon transcriptome index from a FASTA so `sniff strand --index` does not +rebuild one on every run. ```bash bqc sniff adapters reads.cbq # human-readable summary @@ -430,6 +434,7 @@ bqc workflow sample.cbq --config sample.toml -o sample.clean.cbq ### Strandedness ```bash +bqc sniff index --sequence transcripts.fa # reusable index bqc sniff strand reads.cbq --index salmon-index --format json -o strand.json ``` @@ -477,10 +482,12 @@ Pair orientation is reported as measured. An outward or matching library is surfaced with a warning rather than collapsed into inward, because that is exactly the finding worth seeing. -**The index must be a Salmon 2.x index.** An index built by salmon 1.x uses the -pufferfish format and cannot be read; `bqc` says so and gives the rebuild -command. Its reference count, k, decoy state and content hashes are recorded with -the result, so a report can be matched back to the index that produced it. +**The index must be a Salmon 2.x index.** `bqc sniff index --sequence tx.fa` +builds one. An index built by salmon 1.x uses the pufferfish format and cannot +be read; `bqc` says so and gives the rebuild command. `--transcriptome` on +`sniff strand` still builds a throwaway index for a one-shot run. Its reference +count, k, decoy state and content hashes are recorded with the result, so a +report can be matched back to the index that produced it. Sampling stops at a batch boundary once `--target-informative` (50000) observations have accumulated, so a high mapping rate does not cost a full pass β€” diff --git a/bqc/src/cli.rs b/bqc/src/cli.rs index 694be9e..cd2127f 100644 --- a/bqc/src/cli.rs +++ b/bqc/src/cli.rs @@ -10,6 +10,8 @@ use std::collections::HashSet; use std::ops::Range; use std::path::{Path, PathBuf}; +#[cfg(feature = "sniff-strand")] +use std::sync::atomic::{AtomicU64, Ordering}; #[cfg(feature = "sniff-strand")] use clap::ArgGroup; @@ -61,12 +63,13 @@ pub enum Command { Sniff(Box), } -/// Non-destructive inspection. +/// Non-destructive inspection, plus Salmon index prep for `sniff strand`. /// /// These are subcommands rather than flags on one command because their inputs -/// differ: adapter discovery is reference-free and works on any library, while +/// differ: adapter discovery is reference-free and works on any library, /// strandedness needs a transcriptome (or its index) and only means anything -/// for RNA. +/// for RNA, and `index` builds that transcriptome index without inspecting a +/// CBQ. #[derive(Debug, Args)] pub struct SniffCommand { #[command(subcommand)] @@ -77,11 +80,38 @@ pub struct SniffCommand { pub enum SniffKind { /// Infer which adapter sequences contaminate the reads. Adapters(Box), + /// Build a reusable Salmon transcriptome index from a FASTA. + #[cfg(feature = "sniff-strand")] + Index(Box), /// Infer RNA-seq library strandedness against a Salmon transcriptome index. #[cfg(feature = "sniff-strand")] Strand(Box), } +#[cfg(feature = "sniff-strand")] +#[derive(Debug, Args)] +pub struct SniffIndexCommand { + /// Transcriptome FASTA, optionally gzipped. + #[arg(long, value_name = "FASTA")] + pub sequence: PathBuf, + + /// Directory the index is written into. Defaults to the FASTA's directory. + #[arg(long, value_name = "PATH")] + pub output_directory: Option, + + /// Index directory name. Defaults to .index. + #[arg(long, value_name = "NAME")] + pub prefix: Option, + + /// Worker threads; 0 uses every available core. + #[arg(short = 'T', long, value_name = "INT")] + pub threads: Option, + + /// Replace an existing index directory. + #[arg(long)] + pub force: bool, +} + #[cfg(feature = "sniff-strand")] #[derive(Debug, Args)] #[command(group( @@ -1200,10 +1230,191 @@ fn run_sniff(command: &SniffCommand) -> Result { match &command.kind { SniffKind::Adapters(command) => run_sniff_adapters(command), #[cfg(feature = "sniff-strand")] + SniffKind::Index(command) => run_sniff_index(command), + #[cfg(feature = "sniff-strand")] SniffKind::Strand(command) => run_sniff_strand(command), } } +#[cfg(feature = "sniff-strand")] +fn run_sniff_index(command: &SniffIndexCommand) -> Result { + let dest = index_destination( + &command.sequence, + command.output_directory.as_deref(), + command.prefix.as_deref(), + )?; + if !command.sequence.is_file() { + return Err(Error::config(format!( + "{} is not a file; give a transcriptome FASTA", + command.sequence.display() + ))); + } + if dest.exists() && !command.force { + return Err(Error::OutputExists(dest)); + } + let tmp = reserve_index_tmp(&dest)?; + let mut scratch = ScratchDir { + path: tmp, + keep: false, + }; + eprintln!("building Salmon index β†’ {}", dest.display()); + crate::sniff::strand::build_index( + &command.sequence, + &scratch.path, + resolve_threads(command.threads), + )?; + publish_index_dir(&scratch.path, &dest, command.force)?; + scratch.keep = true; + println!("{}", dest.display()); + Ok(Outcome::Success) +} + +/// Last path component of the index directory. +#[cfg(feature = "sniff-strand")] +fn index_destination( + sequence: &Path, + output_directory: Option<&Path>, + prefix: Option<&str>, +) -> Result { + let name = match prefix { + Some(prefix) => validate_index_prefix(prefix)?.to_string(), + None => fasta_index_name(sequence)?, + }; + let parent = match output_directory { + Some(directory) => directory.to_path_buf(), + None => match sequence.parent() { + Some(parent) if !parent.as_os_str().is_empty() => parent.to_path_buf(), + _ => PathBuf::from("."), + }, + }; + Ok(parent.join(name)) +} + +#[cfg(feature = "sniff-strand")] +fn fasta_index_name(sequence: &Path) -> Result { + let name = sequence + .file_name() + .and_then(|name| name.to_str()) + .ok_or_else(|| { + Error::config(format!( + "--sequence {} does not have a usable file name", + sequence.display() + )) + })?; + let stem = fasta_stem(name); + if stem.is_empty() { + return Err(Error::config(format!( + "--sequence {} has no file stem to name the index; pass --prefix", + sequence.display() + ))); + } + Ok(format!("{stem}.index")) +} + +#[cfg(feature = "sniff-strand")] +fn fasta_stem(name: &str) -> &str { + let name = strip_suffix_ci(name, ".gz").unwrap_or(name); + strip_suffix_ci(name, ".fasta") + .or_else(|| strip_suffix_ci(name, ".fna")) + .or_else(|| strip_suffix_ci(name, ".fa")) + .unwrap_or(name) +} + +#[cfg(feature = "sniff-strand")] +fn strip_suffix_ci<'a>(name: &'a str, suffix: &str) -> Option<&'a str> { + let end = name.len().checked_sub(suffix.len())?; + name.get(end..) + .filter(|tail| tail.eq_ignore_ascii_case(suffix)) + .map(|_| &name[..end]) +} + +#[cfg(feature = "sniff-strand")] +fn validate_index_prefix(prefix: &str) -> Result<&str> { + if prefix.is_empty() + || prefix == "." + || prefix == ".." + || prefix.contains('/') + || prefix.contains('\\') + { + return Err(Error::config( + "--prefix is a directory name, not a path; use --output-directory to choose the location", + )); + } + Ok(prefix) +} + +#[cfg(feature = "sniff-strand")] +static INDEX_TMP: AtomicU64 = AtomicU64::new(0); + +#[cfg(feature = "sniff-strand")] +fn reserve_index_tmp(dest: &Path) -> Result { + let name = dest + .file_name() + .ok_or_else(|| Error::config(format!("{} is not a valid output path", dest.display())))?; + let parent = match dest.parent() { + Some(parent) if !parent.as_os_str().is_empty() => parent, + _ => Path::new("."), + }; + std::fs::create_dir_all(parent).map_err(|error| Error::write(parent, error))?; + for _ in 0..128 { + let sequence = INDEX_TMP.fetch_add(1, Ordering::Relaxed); + let tmp = parent.join(format!( + ".{}.bqc-tmp-{}-{sequence}", + name.to_string_lossy(), + std::process::id() + )); + match std::fs::create_dir(&tmp) { + Ok(()) => return Ok(tmp), + Err(error) if error.kind() == std::io::ErrorKind::AlreadyExists => {} + Err(error) => return Err(Error::write(&tmp, error)), + } + } + Err(Error::write( + dest, + std::io::Error::new( + std::io::ErrorKind::AlreadyExists, + "could not reserve a unique temporary index directory", + ), + )) +} + +#[cfg(feature = "sniff-strand")] +fn publish_index_dir(tmp: &Path, dest: &Path, force: bool) -> Result<()> { + if dest.exists() { + if !force { + return Err(Error::OutputExists(dest.to_path_buf())); + } + let remove = if dest.is_dir() { + std::fs::remove_dir_all(dest) + } else { + std::fs::remove_file(dest) + }; + remove.map_err(|error| Error::write(dest, error))?; + } + match std::fs::rename(tmp, dest) { + Ok(()) => Ok(()), + Err(error) if error.kind() == std::io::ErrorKind::AlreadyExists => { + Err(Error::OutputExists(dest.to_path_buf())) + } + Err(error) => Err(Error::write(dest, error)), + } +} + +#[cfg(feature = "sniff-strand")] +struct ScratchDir { + path: PathBuf, + keep: bool, +} + +#[cfg(feature = "sniff-strand")] +impl Drop for ScratchDir { + fn drop(&mut self) { + if !self.keep { + let _ = std::fs::remove_dir_all(&self.path); + } + } +} + #[cfg(feature = "sniff-strand")] fn run_sniff_strand(command: &SniffStrandCommand) -> Result { use crate::sniff::strand; @@ -1768,4 +1979,48 @@ mod tests { assert_eq!(options.min_overlap, Some(10)); assert!(options.r1.is_none()); } + + #[cfg(feature = "sniff-strand")] + #[test] + fn fasta_stems_strip_known_suffixes() { + assert_eq!(fasta_stem("tx.fa"), "tx"); + assert_eq!(fasta_stem("tx.fa.gz"), "tx"); + assert_eq!(fasta_stem("tx.fasta"), "tx"); + assert_eq!(fasta_stem("tx.FASTA.GZ"), "tx"); + assert_eq!(fasta_stem("tx.fna.gz"), "tx"); + assert_eq!(fasta_stem("weird.txt"), "weird.txt"); + assert_eq!(fasta_stem(".fa"), ""); + } + + #[cfg(feature = "sniff-strand")] + #[test] + fn index_destination_composes_directory_and_prefix() { + assert_eq!( + index_destination(Path::new("tx.fa"), None, None).unwrap(), + PathBuf::from("./tx.index") + ); + assert_eq!( + index_destination(Path::new("/ref/tx.fa.gz"), None, None).unwrap(), + PathBuf::from("/ref/tx.index") + ); + assert_eq!( + index_destination(Path::new("tx.fa"), Some(Path::new("/idx")), None).unwrap(), + PathBuf::from("/idx/tx.index") + ); + assert_eq!( + index_destination(Path::new("tx.fa"), None, Some("hg38")).unwrap(), + PathBuf::from("./hg38") + ); + assert_eq!( + index_destination(Path::new("tx.fa"), Some(Path::new("/idx")), Some("hg38")).unwrap(), + PathBuf::from("/idx/hg38") + ); + let err = index_destination(Path::new("tx.fa"), None, Some("a/b")).unwrap_err(); + assert!( + format!("{err}").contains("--prefix is a directory name"), + "{err}" + ); + let err = index_destination(Path::new(".fa"), None, None).unwrap_err(); + assert!(format!("{err}").contains("no file stem"), "{err}"); + } } diff --git a/bqc/src/sniff.rs b/bqc/src/sniff.rs index 02b6246..136dcd8 100644 --- a/bqc/src/sniff.rs +++ b/bqc/src/sniff.rs @@ -1,21 +1,22 @@ // Copyright (c) 2026 The Hiller Lab at the Senckenberg Gesellschaft fΓΌr Naturforschung // Distributed under the terms of the GNU General Public License, Version 3.0. -//! `bqc sniff` β€” non-destructive inspection of a CBQ file. +//! `bqc sniff` β€” non-destructive inspection of a CBQ file, plus index prep. //! -//! Sniffing never trims, filters, reorders or rewrites anything. It samples a -//! deterministic subset of the input, infers a property from it, and writes a -//! report. The input file is opened read-only and is byte-identical afterwards. -//! -//! Two subcommands share the sampling, execution and report conventions and -//! nothing else, because their evidence has nothing in common: +//! `adapters` and `strand` never trim, filter, reorder or rewrite anything. +//! They sample a deterministic subset of the input, infer a property from it, +//! and write a report. The input file is opened read-only and is byte-identical +//! afterwards. //! //! ```text //! sniff adapters reference-free, sequence tails -> candidate adapters +//! sniff index transcriptome FASTA -> reusable Salmon index //! sniff strand needs a transcriptome, mappings -> library orientation //! ``` //! -//! Both may return an inconclusive answer, and both say so rather than guessing. +//! `index` does not inspect a CBQ: it exists so `strand` can reuse a reference +//! instead of rebuilding one on every run. `adapters` and `strand` may return +//! an inconclusive answer, and both say so rather than guessing. //! `--require-confident` turns an inconclusive answer into a distinct exit code //! so a pipeline can branch on it without parsing text. diff --git a/bqc/src/sniff/strand.rs b/bqc/src/sniff/strand.rs index 6a9be35..84a56b7 100644 --- a/bqc/src/sniff/strand.rs +++ b/bqc/src/sniff/strand.rs @@ -379,28 +379,25 @@ fn load_index(path: &std::path::Path) -> Result { }) } -/// Builds a Salmon index from a transcriptome FASTA into a fresh temp dir and -/// returns that dir's path, ready for [`load_index`]. The index stays alive for -/// the process and OS temp cleanup removes it afterwards. -// : no persistent index cache β€” rebuilt every run; cache on -// transcriptome path+mtime when build time starts to matter. -pub fn build_temp_index( +/// Builds a Salmon index from a transcriptome FASTA into `output`. +/// +/// Gzipped FASTA is accepted: `salmon_index` reads through needletail. +/// Persistent indexes are `bqc sniff index`; this is the shared builder. +pub fn build_index( transcriptome: &std::path::Path, + output: &std::path::Path, threads: usize, -) -> Result { +) -> Result<()> { if !transcriptome.is_file() { return Err(Error::config(format!( - "--transcriptome {} is not a file; give a transcriptome FASTA", + "{} is not a file; give a transcriptome FASTA", transcriptome.display() ))); } - let output = std::env::temp_dir().join(format!( - "bqc-strand-{}-{:x}", - std::process::id(), - now_unique() - )); - let mut options = - salmon_index::IndexBuildOptions::new(vec![transcriptome.to_path_buf()], output.clone()); + let mut options = salmon_index::IndexBuildOptions::new( + vec![transcriptome.to_path_buf()], + output.to_path_buf(), + ); options.threads = threads; salmon_index::build(&options).map_err(|error| { Error::config(format!( @@ -408,6 +405,22 @@ pub fn build_temp_index( transcriptome.display() )) })?; + Ok(()) +} + +/// Builds a Salmon index into a fresh temp dir for one `sniff strand` run. +/// Reuse [`build_index`] via `bqc sniff index` when the same reference is +/// needed more than once. +pub fn build_temp_index( + transcriptome: &std::path::Path, + threads: usize, +) -> Result { + let output = std::env::temp_dir().join(format!( + "bqc-strand-{}-{:x}", + std::process::id(), + now_unique() + )); + build_index(transcriptome, &output, threads)?; Ok(output) } diff --git a/bqc/tests/strand.rs b/bqc/tests/strand.rs index b2c7122..f2a3b86 100644 --- a/bqc/tests/strand.rs +++ b/bqc/tests/strand.rs @@ -47,15 +47,15 @@ fn revcomp(sequence: &[u8]) -> Vec { .collect() } -/// A small transcriptome and an index built from it. +/// A small transcriptome and, optionally, an index built from it. struct Reference { dir: TempDir, transcripts: Vec>, } impl Reference { - /// Builds a Salmon index in-process. No `salmon` binary is involved. - fn build() -> Self { + /// Writes a transcriptome FASTA. No index. + fn with_fasta() -> Self { use std::fmt::Write as _; let dir = tempfile::tempdir().expect("tempdir"); let transcripts: Vec> = (0..8).map(|i| transcript(0x7A01 + i * 2, 2000)).collect(); @@ -68,12 +68,17 @@ impl Reference { text.push('\n'); } std::fs::write(&fasta, text).expect("write transcriptome"); + Self { dir, transcripts } + } - let index_dir = dir.path().join("index"); - let mut options = salmon_index::IndexBuildOptions::new(vec![fasta], index_dir); + /// Builds a Salmon index in-process. No `salmon` binary is involved. + fn build() -> Self { + let this = Self::with_fasta(); + let mut options = + salmon_index::IndexBuildOptions::new(vec![this.transcriptome()], this.index()); options.threads = 2; salmon_index::build(&options).expect("build a salmon index"); - Self { dir, transcripts } + this } fn index(&self) -> PathBuf { @@ -398,3 +403,133 @@ fn a_transcriptome_builds_the_index_on_the_fly() { ]) .expect("sniff strand runs with a transcriptome"); } + +#[test] +fn sniff_index_writes_a_reusable_index() { + let reference = Reference::with_fasta(); + let records = reference.pairs(2000, true, PAIRED); + let input = reference.path("reads.cbq"); + write_cbq(&input, PAIRED, &records, BLOCK); + + bqc(&[ + "sniff", + "index", + "--sequence", + reference.transcriptome().to_str().unwrap(), + "-T", + "2", + ]) + .expect("sniff index runs"); + + let index = reference.path("transcripts.index"); + assert!(index.is_dir(), "index directory was not published"); + assert!( + !leftover_tmp(&reference), + "temporary index directory left behind" + ); + + let out = reference.path("strand.json"); + bqc(&[ + "sniff", + "strand", + input.to_str().unwrap(), + "--index", + index.to_str().unwrap(), + "--min-informative", + "100", + "--format", + "json", + "-o", + out.to_str().unwrap(), + ]) + .expect("sniff strand runs against the published index"); + + let result = &json(&out)["result"]; + assert_eq!(result["salmon_library_type"], "ISR", "{result}"); + assert_eq!(result["strandedness"], "reverse", "{result}"); +} + +#[test] +fn sniff_index_respects_output_directory_and_prefix() { + let reference = Reference::with_fasta(); + let outdir = reference.path("indexes"); + bqc(&[ + "sniff", + "index", + "--sequence", + reference.transcriptome().to_str().unwrap(), + "--output-directory", + outdir.to_str().unwrap(), + "--prefix", + "hg38", + "-T", + "2", + ]) + .expect("sniff index runs"); + assert!(outdir.join("hg38").is_dir()); + assert!(!reference.path("transcripts.index").exists()); +} + +#[test] +fn sniff_index_refuses_to_overwrite_without_force() { + let reference = Reference::with_fasta(); + let dest = reference.path("transcripts.index"); + std::fs::create_dir(&dest).expect("mkdir"); + std::fs::write(dest.join("marker"), b"keep").expect("marker"); + + let error = bqc(&[ + "sniff", + "index", + "--sequence", + reference.transcriptome().to_str().unwrap(), + ]) + .expect_err("existing index is refused"); + assert!( + format!("{error}").contains("refusing to overwrite"), + "{error}" + ); + assert_eq!( + std::fs::read(dest.join("marker")).expect("marker survived"), + b"keep" + ); + + bqc(&[ + "sniff", + "index", + "--sequence", + reference.transcriptome().to_str().unwrap(), + "--force", + "-T", + "2", + ]) + .expect("sniff index --force replaces the directory"); + assert!(dest.is_dir()); + assert!(!dest.join("marker").exists(), "old index contents survived"); +} + +#[test] +fn sniff_index_missing_fasta_is_a_config_error() { + let dir = tempfile::tempdir().expect("tempdir"); + let missing = dir.path().join("no-such.fa"); + let error = bqc(&["sniff", "index", "--sequence", missing.to_str().unwrap()]) + .expect_err("missing FASTA is refused"); + let text = format!("{error}"); + assert!(text.contains("is not a file"), "{error}"); +} + +#[test] +fn sniff_index_rejects_a_prefix_that_is_a_path() { + let error = bqc(&["sniff", "index", "--sequence", "tx.fa", "--prefix", "a/b"]) + .expect_err("prefix is not a path"); + assert!( + format!("{error}").contains("--prefix is a directory name"), + "{error}" + ); +} + +fn leftover_tmp(reference: &Reference) -> bool { + std::fs::read_dir(reference.dir.path()) + .expect("list reference dir") + .filter_map(Result::ok) + .any(|entry| entry.file_name().to_string_lossy().contains(".bqc-tmp-")) +}