From 3ba83628be82fc08d143165685fb23b770b57e81 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 26 Aug 2026 23:34:35 +0200 Subject: [PATCH 1/2] feat(cli): six more STAR parameters, each with behaviour rather than acceptance Closes six of the twenty-five STAR 2.7.11b names still listed in NOT_YET_ACCEPTED, chosen as the ones no open theme PR already owns. - --parametersFiles: STAR-format parameter files (name value..., # and // comments). Files are expanded before clap sees the arguments, and a flag the user also passes on the command line is dropped from the file side, so the command line wins even for multi-value parameters where clap would append. An unknown name, an empty value, an unreadable file, or a nested parametersFiles is fatal and names the file and line, as in STAR. - --versionGenome: genomeParameters.txt is checked before the index is read. Versions compare by component, so 2.7.10a sorts after 2.7.4a rather than before it as a string comparison would have it. - --sysShell, and --readFilesCommand now runs through a shell. Spawning the whole command string as one program name meant "gunzip -c", the example in STAR's own documentation, failed as a missing program. Where no POSIX shell is guaranteed the command is split into words instead. - --outFilterMismatchNoverReadLmax: mismatch-to-read-length ratio, SE and PE. - --alignTranscriptsPerReadNmax: caps alignments per read before the filters. - --alignSoftClipAtReferenceEnds: No prohibits clipping past a chromosome end. Machine-checked STAR parameter coverage: 184/203, with the floor asserted by a new test that measures the figure from clap rather than restating it. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 8 + .../content/docs/reference/cli-parameters.md | 8 +- src/align/read_align.rs | 69 ++++ src/index/io.rs | 125 +++++++ src/io/fastq.rs | 164 ++++++++- src/lib.rs | 8 +- src/params/mod.rs | 348 +++++++++++++++++- src/solo/mod.rs | 8 +- tests/cli_surface.rs | 278 ++++++++++++++ tests/parameter_surface.rs | 37 +- 10 files changed, 1028 insertions(+), 25 deletions(-) create mode 100644 tests/cli_surface.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346b..e7d9ec67 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,14 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **Six more STAR 2.7.11b parameters, all with behaviour rather than + acceptance**: `--parametersFiles` (STAR-format parameter files, command + line wins, unknown name or empty value is fatal), `--versionGenome` (an + index older than the requested version is refused instead of misread), + `--sysShell` plus a `--readFilesCommand` that now runs through a shell so + multi-word commands such as `gunzip -c` work, `--outFilterMismatchNoverReadLmax`, + `--alignTranscriptsPerReadNmax`, and `--alignSoftClipAtReferenceEnds`. + Machine-checked STAR parameter coverage rises to 184/203. - **CLI and output parity: SAM/SJ/read-input knobs and the STAR limit surface** — 30 further STAR 2.7.11b parameters. (`--outSAMorder` came from #145.) diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index f749039b..0824f7aa 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -24,13 +24,16 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--genomeSAindexNbases` | `14` | Length of the SA pre-indexing string (log2). Lower for small genomes. | | `--genomeChrBinNbits` | `18` | Log2 of chromosome bin size. | | `--genomeSAsparseD` | `1` | SA sparsity (higher = less RAM, slower mapping). | +| `--versionGenome` | `2.7.4a` | Earliest genome index version accepted; an older index is refused rather than misread. | ## Read input | Parameter | Default | Description | |-----------|---------|-------------| | `--readFilesIn` | — | Input FASTQ file(s); second file is mate 2 for paired-end (required for `alignReads`). | -| `--readFilesCommand` | — | Decompression command, e.g. `zcat` for `.gz`. | +| `--readFilesCommand` | — | Decompression command, e.g. `zcat` or `gunzip -c` for `.gz`. Run through a shell, so multi-word commands work. | +| `--sysShell` | `-` | Shell that runs `--readFilesCommand`. `-` uses `/bin/sh` on Unix; on Windows the command is split into words and run directly. | +| `--parametersFiles` | `-` | STAR-format parameter file(s): `parameterName value...` per line, `#` and `//` comments. Command-line values win. `-` means none. | | `--readMapNumber` | `-1` | Number of reads to map (`-1` = all). | | `--clip5pNbases` | `0` | Bases to clip from the 5' end of each mate. | | `--clip3pNbases` | `0` | Bases to clip from the 3' end of each mate. | @@ -65,6 +68,7 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--outFilterMultimapScoreRange` | `1` | Score range for keeping multi-mappers within best score. | | `--outFilterMismatchNmax` | `10` | Max mismatches per pair. | | `--outFilterMismatchNoverLmax` | `0.3` | Max ratio of mismatches to mapped length. | +| `--outFilterMismatchNoverReadLmax` | `1.0` | Max ratio of mismatches to read length. | | `--outFilterScoreMin` | `0` | Min absolute alignment score. | | `--outFilterScoreMinOverLread` | `0.66` | Min alignment score normalized to read length. | | `--outFilterMatchNmin` | `0` | Min absolute matched bases. | @@ -100,6 +104,8 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--alignSJoverhangMin` | `5` | Min overhang for novel splice junctions. | | `--alignSJDBoverhangMin` | `3` | Min overhang for annotated junctions. | | `--alignSJstitchMismatchNmax` | `0 -1 0 0` | Max mismatches for SJ stitching `[noncan, GC/AG, AT/AC, noncan]`. | +| `--alignTranscriptsPerReadNmax` | `10000` | Max alignments kept per read before the score filters. | +| `--alignSoftClipAtReferenceEnds` | `Yes` | `No` prohibits soft-clipping past a chromosome end (Cufflinks compatibility). | ## Scoring penalties diff --git a/src/align/read_align.rs b/src/align/read_align.rs index ec66a166..4ac4b2a7 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -132,6 +132,27 @@ pub enum PairedAlignmentResult { }, } +/// Whether `t`'s soft-clipped ends would extend past the boundaries of its +/// chromosome: whether the clipped bases have nowhere to sit on the reference. +/// This is what `--alignSoftClipAtReferenceEnds No` prohibits. +fn clips_past_reference_end(t: &Transcript, index: &GenomeIndex, read_len: usize) -> bool { + let Some(first) = t.exons.first() else { + return false; + }; + let Some(last) = t.exons.last() else { + return false; + }; + let chr_start = index.genome.chr_start[t.chr_idx]; + let chr_end = chr_start + index.genome.chr_length[t.chr_idx]; + + // Bases clipped before the first aligned base, and after the last one. + let left_clip = first.read_start as u64; + let right_clip = (read_len - last.read_end) as u64; + + first.genome_start < chr_start.saturating_add(left_clip) && left_clip > 0 + || last.genome_end.saturating_add(right_clip) > chr_end && right_clip > 0 +} + /// Align a read to the genome. /// /// # Algorithm @@ -392,6 +413,19 @@ pub fn align_read( ); } + // Hard cap on how many alignments one read may carry into the filters + // (STAR's alignTranscriptsPerReadNmax). Transcripts are already ordered + // best-score-first, so the cap keeps the best ones. + if transcripts.len() > params.align_transcripts_per_read_nmax { + log::debug!( + "Read {}: {} alignments capped to alignTranscriptsPerReadNmax={}", + read_name, + transcripts.len(), + params.align_transcripts_per_read_nmax + ); + transcripts.truncate(params.align_transcripts_per_read_nmax); + } + // Score-range filter: keep only alignments within outFilterMultimapScoreRange of the best. // (STAR's multMapSelect step — must run before quality filters.) if !transcripts.is_empty() { @@ -410,7 +444,21 @@ pub fn align_read( let pre_filter_count = transcripts.len(); let mut filter_reasons = std::collections::HashMap::new(); + let prohibit_ref_end_clip = params + .align_soft_clip_at_reference_ends + .eq_ignore_ascii_case("No"); + transcripts.retain(|t| { + // A soft clip that would hang past the start or the end of the + // chromosome (STAR's alignSoftClipAtReferenceEnds No, needed for + // Cufflinks-compatible output). + if prohibit_ref_end_clip && clips_past_reference_end(t, index, read_seq.len()) { + *filter_reasons + .entry("soft_clip_at_reference_end") + .or_insert(0) += 1; + return false; + } + // Absolute score threshold if t.score < params.out_filter_score_min { *filter_reasons.entry("score_min").or_insert(0) += 1; @@ -437,6 +485,25 @@ pub fn align_read( return false; } + // Mismatch count over the *read* length (STAR's + // outFilterMismatchNoverReadLmax). STAR takes + // outFilterMismatchNoverLmax over the *mapped* length instead; the + // check above still divides by the read length, so the two behave + // identically here until issue #238 is fixed. + let read_mismatch_rate = t.n_mismatch as f64 / read_length; + if read_mismatch_rate > params.out_filter_mismatch_nover_read_lmax { + *filter_reasons.entry("mismatch_rate_read").or_insert(0) += 1; + log::debug!( + "Filtered {}: {:.1}% read-length mismatch rate > {:.1}% max ({}/{} bases)", + read_name, + read_mismatch_rate * 100.0, + params.out_filter_mismatch_nover_read_lmax * 100.0, + t.n_mismatch, + read_length + ); + return false; + } + // Relative mismatch count (mismatches / read_length) let mismatch_rate = t.n_mismatch as f64 / read_length; if mismatch_rate > params.out_filter_mismatch_nover_lmax { @@ -1475,6 +1542,8 @@ fn filter_paired_transcripts(paired_alns: &mut Vec, params: &Pa if combined_nm > params.out_filter_mismatch_nmax || (combined_nm as f64) > params.out_filter_mismatch_nover_lmax * (mate1_len + mate2_len) + || (combined_nm as f64) + > params.out_filter_mismatch_nover_read_lmax * (mate1_len + mate2_len) { paired_alns.clear(); return; diff --git a/src/index/io.rs b/src/index/io.rs index ce6de640..4ffdce69 100644 --- a/src/index/io.rs +++ b/src/index/io.rs @@ -21,6 +21,11 @@ impl GenomeIndex { pub fn load(genome_dir: &Path, params: &Parameters) -> Result { log::info!("Loading genome from {}...", genome_dir.display()); + // Refuse an index older than --versionGenome before reading a byte of + // it: an older layout read with today's reader is silently wrong, not + // loudly wrong. + check_genome_version(genome_dir, ¶ms.version_genome)?; + // Load Genome file let genome = load_genome(genome_dir, params)?; log::info!( @@ -150,6 +155,74 @@ impl GenomeIndex { } } +/// Parse a STAR genome version string (`2.7.4a`) into a comparable tuple: +/// the numeric components, then the trailing letter suffix (`a` -> 1). +/// +/// STAR compares these versions as strings, which orders `2.7.10a` before +/// `2.7.4a`; comparing components avoids that trap. +fn parse_genome_version(v: &str) -> Option<(Vec, u32)> { + let v = v.trim(); + if v.is_empty() { + return None; + } + let digits_end = v.rfind(|c: char| c.is_ascii_digit()).map_or(0, |i| i + 1); + let (numeric, suffix) = v.split_at(digits_end); + let mut parts = Vec::new(); + for p in numeric.split('.') { + parts.push(p.parse::().ok()?); + } + let suffix_rank = suffix.bytes().next().map_or(0, |b| { + u32::from(b.to_ascii_lowercase().saturating_sub(b'a')) + 1 + }); + Some((parts, suffix_rank)) +} + +/// Read `versionGenome` from `genomeParameters.txt`. +fn read_genome_version(genome_dir: &Path) -> Result, Error> { + let path = genome_dir.join("genomeParameters.txt"); + let contents = match std::fs::read_to_string(&path) { + Ok(s) => s, + Err(e) if e.kind() == std::io::ErrorKind::NotFound => return Ok(None), + Err(e) => return Err(Error::io(e, &path)), + }; + for line in contents.lines() { + if let Some(rest) = line.strip_prefix("versionGenome") + && let Some(value) = rest.split_whitespace().next() + { + return Ok(Some(value.to_string())); + } + } + Ok(None) +} + +/// Fail when the index on disk is older than `required` (`--versionGenome`). +/// +/// An index with no recorded version, or a version neither side can parse, is +/// accepted with a warning: refusing it would break directories that load +/// correctly today, and the version line is advisory metadata, not a checksum. +pub(crate) fn check_genome_version(genome_dir: &Path, required: &str) -> Result<(), Error> { + let Some(found) = read_genome_version(genome_dir)? else { + log::warn!( + "{} has no versionGenome line; skipping the version check", + genome_dir.join("genomeParameters.txt").display() + ); + return Ok(()); + }; + let (Some(found_v), Some(required_v)) = + (parse_genome_version(&found), parse_genome_version(required)) + else { + log::warn!("could not compare genome version '{found}' against '{required}'"); + return Ok(()); + }; + if found_v < required_v { + return Err(Error::Index(format!( + "genome index in {} has versionGenome {found}, older than the required {required}. Regenerate the index with --runMode genomeGenerate, or lower --versionGenome if you know the layout is compatible", + genome_dir.display() + ))); + } + Ok(()) +} + /// Read `genomeFileSizes\t ` from genomeParameters.txt /// and return the first field (total genome byte count, including Gsj if /// sjdb was baked in). Returns `Ok(None)` if the file or line is absent, @@ -403,4 +476,56 @@ mod tests { assert_eq!(loaded_index.suffix_array.get(i), index.suffix_array.get(i)); } } + + // ── --versionGenome ────────────────────────────────────────────────── + + fn write_genome_params(dir: &std::path::Path, version_line: &str) { + std::fs::write(dir.join("genomeParameters.txt"), version_line).unwrap(); + } + + #[test] + fn version_genome_rejects_an_older_index() { + let dir = tempfile::tempdir().unwrap(); + write_genome_params(dir.path(), "versionGenome\t2.7.1a\n"); + let err = check_genome_version(dir.path(), "2.7.4a").unwrap_err(); + let msg = err.to_string(); + assert!( + msg.contains("2.7.1a"), + "error should quote the found version: {msg}" + ); + assert!( + msg.contains("2.7.4a"), + "error should quote the required version: {msg}" + ); + } + + #[test] + fn version_genome_accepts_equal_and_newer_indices() { + let dir = tempfile::tempdir().unwrap(); + write_genome_params(dir.path(), "versionGenome\t2.7.4a\n"); + assert!(check_genome_version(dir.path(), "2.7.4a").is_ok()); + + write_genome_params(dir.path(), "versionGenome\t2.7.10b\n"); + assert!(check_genome_version(dir.path(), "2.7.4a").is_ok()); + } + + #[test] + fn version_genome_orders_by_component_not_lexically() { + // A string comparison puts "2.7.10a" before "2.7.4a"; component + // comparison must not. + let older = parse_genome_version("2.7.4a").unwrap(); + let newer = parse_genome_version("2.7.10a").unwrap(); + assert!(newer > older, "2.7.10a must sort after 2.7.4a"); + assert!(parse_genome_version("2.7.4b").unwrap() > older); + } + + #[test] + fn version_genome_missing_line_is_accepted_with_a_warning() { + let dir = tempfile::tempdir().unwrap(); + write_genome_params(dir.path(), "genomeType\tFull\n"); + assert!(check_genome_version(dir.path(), "2.7.4a").is_ok()); + // No genomeParameters.txt at all is also accepted. + std::fs::remove_file(dir.path().join("genomeParameters.txt")).unwrap(); + assert!(check_genome_version(dir.path(), "2.7.4a").is_ok()); + } } diff --git a/src/io/fastq.rs b/src/io/fastq.rs index 0c795de0..bb0b040d 100644 --- a/src/io/fastq.rs +++ b/src/io/fastq.rs @@ -87,6 +87,75 @@ pub struct FastqReader { name_separators: Vec, } +/// A `--readFilesCommand`, plus the shell it is run through (`--sysShell`). +/// +/// STAR runs the command as `sysShell -c " "`, so a multi-word +/// command such as `gunzip -c` works. Spawning the whole string as one program +/// name does not: it fails with "no such file or directory" on the very +/// commands STAR's own documentation gives as examples. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct ReadCommand { + /// The command line as written by the user, without the input file. + pub command: String, + /// Shell to run it through. `None` means run it directly, splitting the + /// command into words: the fallback where no POSIX shell is guaranteed. + pub shell: Option, +} + +impl ReadCommand { + /// Build from `--readFilesCommand` and `--sysShell`. + /// + /// `sys_shell` is STAR's `-` sentinel for "the platform default", which is + /// `/bin/sh` on Unix and no shell on Windows. + pub fn new(command: &str, sys_shell: &str) -> Self { + let shell = if sys_shell != "-" && !sys_shell.is_empty() { + Some(sys_shell.to_string()) + } else if cfg!(unix) { + Some("/bin/sh".to_string()) + } else { + None + }; + ReadCommand { + command: command.to_string(), + shell, + } + } + + /// The `Command` that streams `path` on stdout. + fn build(&self, path: &Path) -> Result { + if let Some(shell) = &self.shell { + let quoted = shlex::try_quote(&path.to_string_lossy()) + .map_err(|_| { + Error::Parameter(format!( + "input path contains a NUL byte and cannot be passed to \ + --readFilesCommand: {}", + path.display() + )) + })? + .into_owned(); + let mut c = Command::new(shell); + c.arg("-c").arg(format!("{} {quoted}", self.command)); + return Ok(c); + } + + // No shell: split the command line into words ourselves so that a + // multi-word command still runs. + let mut words = shlex::split(&self.command).ok_or_else(|| { + Error::Parameter(format!( + "could not parse --readFilesCommand '{}' into words", + self.command + )) + })?; + if words.is_empty() { + return Err(Error::Parameter("--readFilesCommand is empty".to_string())); + } + let program = words.remove(0); + let mut c = Command::new(program); + c.args(words).arg(path); + Ok(c) + } +} + impl FastqReader { /// Open a FASTQ file (plain or gzip compressed) /// @@ -96,7 +165,7 @@ impl FastqReader { /// /// # Returns /// A FastqReader that iterates over encoded reads - pub fn open(path: &Path, decompress_cmd: Option<&str>) -> Result { + pub fn open(path: &Path, decompress_cmd: Option<&ReadCommand>) -> Result { let reader: Box = if let Some(cmd) = decompress_cmd { // Use external decompression command Self::open_with_command(path, cmd)? @@ -153,9 +222,9 @@ impl FastqReader { } /// Open FASTQ file using external decompression command - fn open_with_command(path: &Path, cmd: &str) -> Result, Error> { - let mut child = Command::new(cmd) - .arg(path) + fn open_with_command(path: &Path, cmd: &ReadCommand) -> Result, Error> { + let mut child = cmd + .build(path)? .stdout(Stdio::piped()) .spawn() .map_err(|e| Error::io(e, path))?; @@ -250,7 +319,11 @@ impl PairedFastqReader { /// /// # Returns /// A PairedFastqReader that iterates over paired reads with name validation - pub fn open(path1: &Path, path2: &Path, decompress_cmd: Option<&str>) -> Result { + pub fn open( + path1: &Path, + path2: &Path, + decompress_cmd: Option<&ReadCommand>, + ) -> Result { let reader1 = FastqReader::open(path1, decompress_cmd)?; let reader2 = FastqReader::open(path2, decompress_cmd)?; @@ -778,4 +851,85 @@ mod tests { let batch3 = reader.read_paired_batch(3).unwrap(); assert_eq!(batch3.len(), 0); } + + // ── --readFilesCommand / --sysShell ────────────────────────────────── + + fn write_gzipped_fastq(dir: &std::path::Path) -> std::path::PathBuf { + use flate2::Compression; + use flate2::write::GzEncoder; + let path = dir.join("reads.fq.gz"); + let f = std::fs::File::create(&path).unwrap(); + let mut enc = GzEncoder::new(f, Compression::fast()); + enc.write_all(b"@r1\nACGTACGTAC\n+\nIIIIIIIIII\n").unwrap(); + enc.finish().unwrap(); + path + } + + #[cfg(unix)] + #[test] + fn read_files_command_accepts_a_multi_word_command() { + // `gunzip -c` is the example STAR's own documentation gives. Spawning + // the whole string as one program name fails; running it through a + // shell is what STAR does. + let dir = tempfile::tempdir().unwrap(); + let path = write_gzipped_fastq(dir.path()); + let cmd = ReadCommand::new("gunzip -c", "-"); + let mut reader = FastqReader::open(&path, Some(&cmd)).unwrap(); + let record = reader.next_encoded().unwrap().expect("one record"); + assert_eq!(record.name, "r1"); + assert_eq!(record.sequence.len(), 10); + } + + #[cfg(unix)] + #[test] + fn read_files_command_runs_through_the_shell_sys_shell_selects() { + let dir = tempfile::tempdir().unwrap(); + let path = write_gzipped_fastq(dir.path()); + let cmd = ReadCommand::new("gunzip -c", "/bin/sh"); + assert_eq!(cmd.shell.as_deref(), Some("/bin/sh")); + let mut reader = FastqReader::open(&path, Some(&cmd)).unwrap(); + assert!(reader.next_encoded().unwrap().is_some()); + } + + #[cfg(unix)] + #[test] + fn read_files_command_handles_a_path_with_spaces() { + let dir = tempfile::tempdir().unwrap(); + let sub = dir.path().join("a directory with spaces"); + std::fs::create_dir(&sub).unwrap(); + let path = write_gzipped_fastq(&sub); + let cmd = ReadCommand::new("gunzip -c", "-"); + let mut reader = FastqReader::open(&path, Some(&cmd)).unwrap(); + assert!( + reader.next_encoded().unwrap().is_some(), + "the file path must be quoted when handed to the shell" + ); + } + + #[test] + fn read_files_command_without_a_shell_splits_the_command_itself() { + // The Windows path: no POSIX shell is guaranteed, so the command line + // is split into program plus arguments rather than handed to a shell. + let cmd = ReadCommand { + command: "gunzip -c".to_string(), + shell: None, + }; + let built = cmd.build(std::path::Path::new("reads.fq.gz")).unwrap(); + assert_eq!(built.get_program(), "gunzip"); + let args: Vec<_> = built + .get_args() + .map(|a| a.to_string_lossy().into_owned()) + .collect(); + assert_eq!(args, vec!["-c".to_string(), "reads.fq.gz".to_string()]); + } + + #[test] + fn sys_shell_default_is_platform_dependent() { + let cmd = ReadCommand::new("zcat", "-"); + if cfg!(unix) { + assert_eq!(cmd.shell.as_deref(), Some("/bin/sh")); + } else { + assert_eq!(cmd.shell, None); + } + } } diff --git a/src/lib.rs b/src/lib.rs index 5086fcbb..08677ae0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -525,7 +525,8 @@ fn run_smartseq( counts.add(ci, g); } }; - let cmd = params.read_files_command.as_deref(); + let read_command = params.read_command(); + let cmd = read_command.as_ref(); for (ci, cell) in cells.iter().enumerate() { match &cell.read2 { @@ -1414,8 +1415,7 @@ fn align_reads_single_end( let read_file = ¶ms.read_files_in[0]; info!("Reading single-end from {}", read_file.display()); - let reader = - FastqReader::open(read_file, params.read_files_command.as_deref())?.with_params(params); + let reader = FastqReader::open(read_file, params.read_command().as_ref())?.with_params(params); // Create chimeric output writer if enabled let chimeric_writer = if params.chim_segment_min > 0 && params.chim_out_junctions() { @@ -2753,7 +2753,7 @@ fn align_reads_paired_end( let reader = PairedFastqReader::open( ¶ms.read_files_in[0], ¶ms.read_files_in[1], - params.read_files_command.as_deref(), + params.read_command().as_ref(), )? .with_params(params); diff --git a/src/params/mod.rs b/src/params/mod.rs index 3b507310..58717d12 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -530,6 +530,24 @@ pub struct Parameters { #[arg(long = "readFilesIn", num_args = 1..=2)] pub read_files_in: Vec, + /// STAR-format parameter file(s): lines of `parameterName value...`, + /// `#` and `//` comments, one parameter per line. `-` means none. Values + /// given on the command line override values read from a file, and + /// `parametersFiles` itself may only appear on the command line. + #[arg(long = "parametersFiles", num_args = 1.., default_value = "-")] + pub parameters_files: Vec, + + /// Shell used to run `--readFilesCommand`; `-` selects the platform + /// default (`/bin/sh` on Unix). STAR's `sysShell`. + #[arg(long = "sysShell", default_value = "-")] + pub sys_shell: String, + + /// Earliest genome index version this run accepts. A genome directory + /// written by an older STAR / rustar-aligner is rejected rather than read + /// with the wrong layout. + #[arg(long = "versionGenome", default_value = "2.7.4a")] + pub version_genome: String, + /// Command to decompress input files (e.g. "zcat" for .gz) #[arg(long = "readFilesCommand")] pub read_files_command: Option, @@ -855,6 +873,12 @@ pub struct Parameters { #[arg(long = "outFilterMismatchNoverLmax", default_value_t = 0.3)] pub out_filter_mismatch_nover_lmax: f64, + /// Max ratio of mismatches to *read* length (STAR's + /// `outFilterMismatchNoverReadLmax`). Applied to the full read length, + /// where `outFilterMismatchNoverLmax` is applied to the mapped length. + #[arg(long = "outFilterMismatchNoverReadLmax", default_value_t = 1.0)] + pub out_filter_mismatch_nover_read_lmax: f64, + /// Min alignment score (absolute) #[arg(long = "outFilterScoreMin", default_value_t = 0)] pub out_filter_score_min: i32, @@ -931,6 +955,18 @@ pub struct Parameters { #[arg(long = "alignEndsType", default_value = "Local")] pub align_ends_type: String, + /// Max number of alignments kept for one read before the score filters + /// run (STAR's `alignTranscriptsPerReadNmax`). A read that produces more + /// keeps only the highest-scoring `N`. + #[arg(long = "alignTranscriptsPerReadNmax", default_value_t = 10000)] + pub align_transcripts_per_read_nmax: usize, + + /// `Yes` (default) allows an alignment to soft-clip past the end of a + /// chromosome; `No` prohibits it, which is what Cufflinks-compatible + /// output needs. + #[arg(long = "alignSoftClipAtReferenceEnds", default_value = "Yes")] + pub align_soft_clip_at_reference_ends: String, + /// Min overlap (bases) between mates required to trigger merge-and-realign; 0 = off #[arg(long = "peOverlapNbasesMin", default_value_t = 0)] pub pe_overlap_nbases_min: u64, @@ -1499,6 +1535,166 @@ impl Parameters { } /// Parse and validate parameter combinations. + /// Long-flag names clap accepts, walked recursively so flattened argument + /// groups are included. Shared with the parameter-file loader so a file + /// can only set parameters this build actually has. + fn recognised_long_flags(command: &clap::Command) -> std::collections::BTreeSet { + fn walk(cmd: &clap::Command, out: &mut std::collections::BTreeSet) { + for arg in cmd.get_arguments() { + if let Some(long) = arg.get_long() { + out.insert(long.to_string()); + } + for alias in arg.get_all_aliases().unwrap_or_default() { + out.insert(alias.to_string()); + } + } + for sub in cmd.get_subcommands() { + walk(sub, out); + } + } + let mut out = std::collections::BTreeSet::new(); + walk(command, &mut out); + out + } + + /// Parse one STAR parameter file into `--name value...` arguments. + /// + /// STAR's `Parameters::scanOneLine`: an empty line is skipped, a first + /// token starting with `#` or `//` is a comment, an unrecognised name is + /// fatal, and a name with no value is fatal. `parametersFiles` inside a + /// file is refused, as in STAR, where it is command-line only. + fn parse_parameters_file( + path: &str, + known: &std::collections::BTreeSet, + command: &clap::Command, + ) -> Result, clap::Error> { + use clap::error::ErrorKind; + + let text = std::fs::read_to_string(path).map_err(|e| { + let mut cmd = command.clone(); + cmd.error( + ErrorKind::Io, + format!("could not read --parametersFiles '{path}': {e}"), + ) + })?; + + let mut out: Vec = Vec::new(); + for (lineno, line) in text.lines().enumerate() { + let mut words = line.split_whitespace(); + let Some(name) = words.next() else { + continue; // blank line + }; + if name.starts_with('#') || name.starts_with("//") { + continue; // comment + } + if name == "parametersFiles" { + let mut cmd = command.clone(); + return Err(cmd.error( + ErrorKind::InvalidValue, + format!( + "parametersFiles cannot be set inside a parameter file ({path}, line {}); it is a command-line-only parameter", + lineno + 1 + ), + )); + } + if !known.contains(name) { + let mut cmd = command.clone(); + return Err(cmd.error( + ErrorKind::UnknownArgument, + format!( + "unrecognised parameter name '{name}' in --parametersFiles {path} (line {})", + lineno + 1 + ), + )); + } + let values: Vec<&str> = words.collect(); + if values.is_empty() { + let mut cmd = command.clone(); + return Err(cmd.error( + ErrorKind::InvalidValue, + format!( + "empty value for parameter '{name}' in --parametersFiles {path} (line {})", + lineno + 1 + ), + )); + } + out.push(std::ffi::OsString::from(format!("--{name}"))); + for v in values { + out.push(std::ffi::OsString::from(v)); + } + } + Ok(out) + } + + /// Expand `--parametersFiles` into leading arguments. + /// + /// Returns `args` unchanged when no file is requested. Otherwise the file + /// arguments are placed before the user's own arguments, and any flag the + /// user also passed on the command line is dropped from the file side, so + /// the command line wins even for the multi-value parameters where clap + /// would otherwise append rather than replace. + fn expand_parameters_files( + args: &[std::ffi::OsString], + command: &clap::Command, + ) -> Result, clap::Error> { + let lossy: Vec = args + .iter() + .map(|a| a.to_string_lossy().into_owned()) + .collect(); + + // Collect the values of every --parametersFiles occurrence. + let mut files: Vec = Vec::new(); + let mut i = 0; + while i < lossy.len() { + if lossy[i] == "--parametersFiles" { + let mut j = i + 1; + while j < lossy.len() && !lossy[j].starts_with("--") { + files.push(lossy[j].clone()); + j += 1; + } + i = j; + } else { + i += 1; + } + } + files.retain(|f| f != "-"); + if files.is_empty() { + return Ok(args.to_vec()); + } + + let known = Self::recognised_long_flags(command); + let mut file_args: Vec = Vec::new(); + for f in &files { + file_args.extend(Self::parse_parameters_file(f, &known, command)?); + } + + // Flags the user set directly; the file must not also set them. + let user_flags: std::collections::BTreeSet<&str> = + lossy.iter().filter_map(|a| a.strip_prefix("--")).collect(); + let mut filtered: Vec = Vec::new(); + let mut skipping = false; + for a in file_args { + let text = a.to_string_lossy().into_owned(); + if let Some(flag) = text.strip_prefix("--") { + skipping = user_flags.contains(flag); + if skipping { + continue; + } + filtered.push(a); + } else if !skipping { + filtered.push(a); + } + } + + let mut out: Vec = Vec::with_capacity(filtered.len() + args.len()); + if let Some(program) = args.first() { + out.push(program.clone()); + } + out.extend(filtered); + out.extend(args.iter().skip(1).cloned()); + Ok(out) + } + pub fn try_parse() -> Result { Self::try_parse_from(std::env::args_os()) } @@ -1512,7 +1708,11 @@ impl Parameters { let args: Vec<_> = args.into_iter().map(Into::into).collect(); let mut command = ::command(); - let matches = command.clone().get_matches_from(args.iter()); + // `--parametersFiles` is expanded before clap sees the arguments: + // file values become leading arguments, so a value repeated on the + // command line wins (STAR reads files first, then the command line). + let expanded = Self::expand_parameters_files(&args, &command)?; + let matches = command.clone().get_matches_from(expanded.iter()); let mut params = ::from_arg_matches(&matches)?; params.command_line = { @@ -1669,6 +1869,35 @@ impl Parameters { } } + // alignSoftClipAtReferenceEnds is Yes/No; anything else would be + // silently read as Yes. + if !["Yes", "No"] + .iter() + .any(|v| v.eq_ignore_ascii_case(¶ms.align_soft_clip_at_reference_ends)) + { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unsupported --alignSoftClipAtReferenceEnds '{}'; expected Yes or No", + params.align_soft_clip_at_reference_ends + ), + )); + } + + if params.align_transcripts_per_read_nmax == 0 { + return Err(command.error( + ErrorKind::InvalidValue, + "--alignTranscriptsPerReadNmax must be > 0", + )); + } + + if params.out_filter_mismatch_nover_read_lmax < 0.0 { + return Err(command.error( + ErrorKind::InvalidValue, + "--outFilterMismatchNoverReadLmax must be >= 0", + )); + } + // quantMode GeneCounts requires a GTF file if params.quant_gene_counts() && params.sjdb_gtf_file.is_none() { return Err(command.error( @@ -2069,6 +2298,17 @@ impl Parameters { Ok(params) } + /// The `--readFilesCommand`, bound to the shell `--sysShell` selects. + /// + /// `None` when no command was given, in which case compression is + /// detected from the file extension. + pub fn read_command(&self) -> Option { + self.read_files_command + .as_deref() + .filter(|c| !c.trim().is_empty() && *c != "-") + .map(|c| crate::io::fastq::ReadCommand::new(c, &self.sys_shell)) + } + /// Returns true if `--quantMode GeneCounts` was requested. pub fn quant_gene_counts(&self) -> bool { self.quant_mode.iter().any(|m| m == "GeneCounts") @@ -2185,6 +2425,112 @@ mod tests { Parameters::try_parse_from(&full) } + // ── --parametersFiles ──────────────────────────────────────────────── + + fn write_param_file(dir: &std::path::Path, name: &str, body: &str) -> String { + let path = dir.join(name); + std::fs::write(&path, body).unwrap(); + path.to_string_lossy().into_owned() + } + + #[test] + fn parameters_file_sets_values_and_skips_comments() { + let dir = tempfile::tempdir().unwrap(); + let f = write_param_file( + dir.path(), + "star.params", + "# a comment\n\ + // another comment\n\ + \n\ + outFilterMismatchNmax 3\n\ + runThreadN 4\n\ + outSAMattributes NH HI AS nM\n", + ); + let p = try_parse(&["--readFilesIn", "r.fq", "--parametersFiles", &f]).unwrap(); + assert_eq!(p.out_filter_mismatch_nmax, 3); + assert_eq!(p.run_thread_n, NonZeroUsize::new(4).unwrap()); + } + + #[test] + fn parameters_file_is_overridden_by_the_command_line() { + let dir = tempfile::tempdir().unwrap(); + let f = write_param_file(dir.path(), "star.params", "outFilterMismatchNmax 3\n"); + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--parametersFiles", + &f, + "--outFilterMismatchNmax", + "7", + ]) + .unwrap(); + assert_eq!(p.out_filter_mismatch_nmax, 7); + } + + #[test] + fn parameters_file_multi_value_is_replaced_not_appended() { + // clap appends repeated occurrences of a multi-value argument, so the + // loader has to drop the file's copy rather than rely on ordering. + let dir = tempfile::tempdir().unwrap(); + let f = write_param_file(dir.path(), "star.params", "readFilesIn a.fq b.fq\n"); + let p = try_parse(&["--parametersFiles", &f, "--readFilesIn", "c.fq"]).unwrap(); + assert_eq!(p.read_files_in, vec![PathBuf::from("c.fq")]); + } + + #[test] + fn parameters_file_dash_means_none() { + let p = try_parse(&["--readFilesIn", "r.fq", "--parametersFiles", "-"]).unwrap(); + assert_eq!(p.parameters_files, vec!["-".to_string()]); + } + + #[test] + fn parameters_file_rejects_unknown_parameter_name() { + let dir = tempfile::tempdir().unwrap(); + let f = write_param_file(dir.path(), "star.params", "notAParameter 1\n"); + let err = try_parse(&["--readFilesIn", "r.fq", "--parametersFiles", &f]).unwrap_err(); + assert!( + err.to_string().contains("notAParameter"), + "error should name the offending parameter: {err}" + ); + } + + #[test] + fn parameters_file_rejects_empty_value() { + let dir = tempfile::tempdir().unwrap(); + let f = write_param_file(dir.path(), "star.params", "outFilterMismatchNmax\n"); + let err = try_parse(&["--readFilesIn", "r.fq", "--parametersFiles", &f]).unwrap_err(); + assert!( + err.to_string().contains("empty value"), + "error should say the value is empty: {err}" + ); + } + + #[test] + fn parameters_file_rejects_nested_parameters_files() { + let dir = tempfile::tempdir().unwrap(); + let f = write_param_file(dir.path(), "star.params", "parametersFiles other.params\n"); + let err = try_parse(&["--readFilesIn", "r.fq", "--parametersFiles", &f]).unwrap_err(); + assert!( + err.to_string().contains("command-line-only"), + "error should explain the restriction: {err}" + ); + } + + #[test] + fn parameters_file_missing_is_an_error_naming_the_path() { + let err = try_parse(&[ + "--readFilesIn", + "r.fq", + "--parametersFiles", + "/nonexistent/star.params", + ]) + .unwrap_err(); + assert!( + err.to_string().contains("/nonexistent/star.params"), + "error should name the unreadable file: {err}" + ); + } + #[test] fn defaults() { let p = try_parse(&["--readFilesIn", "reads.fq"]).unwrap(); diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 22d63e02..4d8f659b 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -239,7 +239,7 @@ impl SoloReadReader { cdna_path: &Path, barcode_path: &Path, layout: SoloBarcodeLayout, - decompress_cmd: Option<&str>, + decompress_cmd: Option<&crate::io::fastq::ReadCommand>, ) -> Result { Ok(Self { cdna: FastqReader::open(cdna_path, decompress_cmd)?, @@ -304,7 +304,7 @@ pub fn open_reader(params: &Parameters) -> Result { )) })?; let layout = SoloBarcodeLayout::from_params(params); - SoloReadReader::open(cdna, barcode, layout, params.read_files_command.as_deref()) + SoloReadReader::open(cdna, barcode, layout, params.read_command().as_ref()) } /// One paired-end solo read for `--soloBarcodeMate 1` (5' 10x): both mates carry @@ -329,7 +329,7 @@ impl SoloPairedReader { mate1_path: &Path, mate2_path: &Path, layout: SoloBarcodeLayout, - decompress_cmd: Option<&str>, + decompress_cmd: Option<&crate::io::fastq::ReadCommand>, ) -> Result { Ok(Self { mate1: FastqReader::open(mate1_path, decompress_cmd)?, @@ -380,7 +380,7 @@ pub fn open_paired_reader(params: &Parameters) -> Result Vec { + let bases: [u8; 4] = *b"ACGT"; + let mut state = seed; + let mut seq = Vec::with_capacity(length); + for _ in 0..length { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + seq.push(bases[((state >> 16) & 3) as usize]); + } + seq +} + +fn build_genome() -> Vec { + let mut genome = lcg_seq(88888, 20_000); + let block: Vec = genome[REPEAT_SRC..REPEAT_SRC + READ_LEN].to_vec(); + for &dst in &REPEAT_COPIES { + genome[dst..dst + READ_LEN].copy_from_slice(&block); + } + genome +} + +fn write_fasta(dir: &Path, genome: &[u8]) -> PathBuf { + let path = dir.join("genome.fa"); + let mut f = fs::File::create(&path).unwrap(); + writeln!(f, ">chr1").unwrap(); + f.write_all(genome).unwrap(); + writeln!(f).unwrap(); + path +} + +fn build_index(fasta: &Path, genome_dir: &Path) { + fs::create_dir_all(genome_dir).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "genomeGenerate", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--genomeFastaFiles", + fasta.to_str().unwrap(), + "--genomeSAindexNbases", + "7", + ]) + .assert() + .success(); +} + +fn write_fastq(path: &Path, reads: &[(String, Vec)]) { + let mut f = fs::File::create(path).unwrap(); + for (name, seq) in reads { + writeln!(f, "@{name}").unwrap(); + f.write_all(seq).unwrap(); + writeln!(f).unwrap(); + writeln!(f, "+").unwrap(); + writeln!(f, "{}", "I".repeat(seq.len())).unwrap(); + } +} + +/// Run the aligner and return the SAM records (header lines dropped). +fn align(genome_dir: &Path, fastq: &Path, prefix: &str, extra: &[&str]) -> Vec { + let mut cmd = cargo_bin_cmd!("rustar-aligner"); + cmd.args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq.to_str().unwrap(), + "--outFileNamePrefix", + prefix, + ]); + cmd.args(extra); + cmd.assert().success(); + + fs::read_to_string(format!("{prefix}Aligned.out.sam")) + .unwrap() + .lines() + .filter(|l| !l.starts_with('@')) + .map(str::to_string) + .collect() +} + +fn mismatched_read(genome: &[u8], n_mismatches: usize) -> Vec { + let mut seq = genome[UNIQUE_START..UNIQUE_START + READ_LEN].to_vec(); + // Spread the substitutions out so they cannot all be soft-clipped off one + // end instead of being counted as mismatches. + for i in 0..n_mismatches { + let pos = 15 + i * 20; + seq[pos] = match seq[pos] { + b'A' => b'C', + b'C' => b'G', + b'G' => b'T', + _ => b'A', + }; + } + seq +} + +// ── outFilterMismatchNoverReadLmax ────────────────────────────────────────── + +#[test] +fn mismatch_nover_read_lmax_filters_by_read_length_ratio() { + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = build_genome(); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + + let fq = root.join("mm.fq"); + write_fastq(&fq, &[("mm3".to_string(), mismatched_read(&genome, 3))]); + + // Control: the read maps with the default ratio of 1.0. + let permissive = align( + &genome_dir, + &fq, + &format!("{}/permissive_", root.display()), + &["--outFilterMismatchNoverReadLmax", "1.0"], + ); + assert_eq!(permissive.len(), 1, "read should map by default"); + let control_flag: u32 = permissive[0].split('\t').nth(1).unwrap().parse().unwrap(); + assert_eq!( + control_flag & 0x4, + 0, + "control read must be mapped: {}", + permissive[0] + ); + + // 0.01 of 100 bases allows one mismatch; the read carries three. + let strict = align( + &genome_dir, + &fq, + &format!("{}/strict_", root.display()), + &["--outFilterMismatchNoverReadLmax", "0.01"], + ); + let mapped = strict + .iter() + .filter(|r| { + let flag: u32 = r.split('\t').nth(1).unwrap().parse().unwrap(); + flag & 0x4 == 0 + }) + .count(); + assert_eq!( + mapped, 0, + "a 3-mismatch read must fail a 0.01 read-length ratio: {strict:?}" + ); +} + +// ── alignTranscriptsPerReadNmax ───────────────────────────────────────────── + +#[test] +fn transcripts_per_read_nmax_caps_the_alignments_kept() { + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = build_genome(); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + + let fq = root.join("multi.fq"); + let repeat = genome[REPEAT_SRC..REPEAT_SRC + READ_LEN].to_vec(); + write_fastq(&fq, &[("multi".to_string(), repeat)]); + + // Control: all five copies are reported when nothing caps them. + let uncapped = align( + &genome_dir, + &fq, + &format!("{}/uncapped_", root.display()), + &["--outFilterMultimapNmax", "20"], + ); + assert!( + uncapped.len() >= 3, + "the planted repeat should multimap, got {} record(s)", + uncapped.len() + ); + + let capped = align( + &genome_dir, + &fq, + &format!("{}/capped_", root.display()), + &[ + "--outFilterMultimapNmax", + "20", + "--alignTranscriptsPerReadNmax", + "2", + ], + ); + assert!( + capped.len() < uncapped.len(), + "the cap must reduce the alignments kept: {} capped vs {} uncapped", + capped.len(), + uncapped.len() + ); + assert!( + capped.len() <= 2, + "at most 2 alignments may survive --alignTranscriptsPerReadNmax 2, got {}", + capped.len() + ); +} + +// ── alignSoftClipAtReferenceEnds ──────────────────────────────────────────── + +#[test] +fn soft_clip_at_reference_ends_no_prohibits_clipping_past_the_chromosome() { + let dir = TempDir::new().unwrap(); + let root = dir.path(); + let genome = build_genome(); + let genome_dir = root.join("genome"); + build_index(&write_fasta(root, &genome), &genome_dir); + + // A read whose last 20 bases fall past the end of chr1: 80 bases of the + // genome tail, then 20 bases that exist nowhere, so the aligner has to + // soft-clip them past the reference end. + let mut seq = genome[genome.len() - 80..].to_vec(); + seq.extend_from_slice(b"ACACACACACGTGTGTGTGT"); + let fq = root.join("edge.fq"); + write_fastq(&fq, &[("edge".to_string(), seq)]); + + let allowed = align( + &genome_dir, + &fq, + &format!("{}/allowed_", root.display()), + &["--alignSoftClipAtReferenceEnds", "Yes"], + ); + let mapped_allowed = allowed + .iter() + .filter(|r| { + let flag: u32 = r.split('\t').nth(1).unwrap().parse().unwrap(); + flag & 0x4 == 0 + }) + .count(); + assert_eq!( + mapped_allowed, 1, + "control: the read maps with a soft clip past the end: {allowed:?}" + ); + assert!( + allowed[0].split('\t').nth(5).unwrap().contains('S'), + "control alignment should carry a soft clip: {}", + allowed[0] + ); + + let prohibited = align( + &genome_dir, + &fq, + &format!("{}/prohibited_", root.display()), + &["--alignSoftClipAtReferenceEnds", "No"], + ); + let mapped_prohibited = prohibited + .iter() + .filter(|r| { + let flag: u32 = r.split('\t').nth(1).unwrap().parse().unwrap(); + flag & 0x4 == 0 + }) + .count(); + assert_eq!( + mapped_prohibited, 0, + "--alignSoftClipAtReferenceEnds No must reject it: {prohibited:?}" + ); +} diff --git a/tests/parameter_surface.rs b/tests/parameter_surface.rs index df3dd940..dd2a6bc4 100644 --- a/tests/parameter_surface.rs +++ b/tests/parameter_surface.rs @@ -84,20 +84,10 @@ const ACCEPTED_BUT_INERT: &[(&str, &str)] = &[ /// Adding a name here must always be a deliberate act. Removing one is what /// progress looks like. const NOT_YET_ACCEPTED: &[&str] = &[ - // STAR meta-parameters, none of them accepted here. clap rejects them, so - // a user who passes one is told rather than quietly ignored, which is the - // behaviour these three need most: silently dropping `--parametersFiles` - // would discard every parameter in that file. - "parametersFiles", - "sysShell", - "versionGenome", // Aligner core (annotated-junction stitching, alignEndsType, in-recursion // length penalty). "alignEndsProtrude", "alignInsertionFlush", - "alignSoftClipAtReferenceEnds", - "alignTranscriptsPerReadNmax", - "outFilterMismatchNoverReadLmax", "seedNoneLociPerWindow", "seedSplitMin", // Long reads. @@ -192,6 +182,33 @@ fn star_parameter_surface_is_fully_accounted_for() { ); } +/// The coverage figure itself, measured from clap rather than asserted by +/// hand: how many of STAR 2.7.11b's parameter names this CLI accepts. The +/// denominator is whatever `star_2.7.11b_params.txt` holds, not a number +/// written here. +/// +/// The floor rises as the port advances. It exists so that a regression that +/// silently drops a parameter fails here rather than in a user's pipeline. +#[test] +fn star_parameter_coverage_meets_the_floor() { + const FLOOR: usize = 180; + + let star = star_parameter_names(); + let ours = recognised_flags(); + let accepted = star.iter().filter(|n| ours.contains(n.as_str())).count(); + + println!( + "STAR parameter coverage: {accepted}/{} accepted, {} still missing", + star.len(), + star.len() - accepted + ); + assert!( + accepted >= FLOOR, + "STAR parameter coverage fell to {accepted}/{}; the floor is {FLOOR}", + star.len() + ); +} + #[test] fn inert_parameters_are_actually_accepted() { // A name in ACCEPTED_BUT_INERT that the CLI does not accept would be a From 60e92e050580b38a3b53d719d4c49a777c4251b0 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 27 Aug 2026 00:00:59 +0200 Subject: [PATCH 2/2] refactor(cli): narrow this PR to the parameters no other PR owns The aligner-core PR implements outFilterMismatchNoverReadLmax, alignTranscriptsPerReadNmax and alignSoftClipAtReferenceEnds as part of its own theme, and it is the older and larger change, so this PR yields them rather than racing it. What remains here is the parameter-file loader, the genome version check, and the shell-run readFilesCommand fix. Machine-checked coverage for this branch alone: 181/203. --- CHANGELOG.md | 12 +- .../content/docs/reference/cli-parameters.md | 3 - src/align/read_align.rs | 69 ----- src/params/mod.rs | 47 --- tests/cli_surface.rs | 278 ------------------ tests/parameter_surface.rs | 8 +- 6 files changed, 12 insertions(+), 405 deletions(-) delete mode 100644 tests/cli_surface.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index e7d9ec67..effbe06c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,14 +21,14 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features -- **Six more STAR 2.7.11b parameters, all with behaviour rather than +- **Three more STAR 2.7.11b parameters, all with behaviour rather than acceptance**: `--parametersFiles` (STAR-format parameter files, command line wins, unknown name or empty value is fatal), `--versionGenome` (an - index older than the requested version is refused instead of misread), - `--sysShell` plus a `--readFilesCommand` that now runs through a shell so - multi-word commands such as `gunzip -c` work, `--outFilterMismatchNoverReadLmax`, - `--alignTranscriptsPerReadNmax`, and `--alignSoftClipAtReferenceEnds`. - Machine-checked STAR parameter coverage rises to 184/203. + index older than the requested version is refused instead of misread), and + `--sysShell`, alongside a `--readFilesCommand` fix: it now runs through a + shell, so multi-word commands such as `gunzip -c` work instead of failing + as a missing program. Machine-checked STAR parameter coverage rises to + 181/203. - **CLI and output parity: SAM/SJ/read-input knobs and the STAR limit surface** — 30 further STAR 2.7.11b parameters. (`--outSAMorder` came from #145.) diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index 0824f7aa..d97b61bc 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -68,7 +68,6 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--outFilterMultimapScoreRange` | `1` | Score range for keeping multi-mappers within best score. | | `--outFilterMismatchNmax` | `10` | Max mismatches per pair. | | `--outFilterMismatchNoverLmax` | `0.3` | Max ratio of mismatches to mapped length. | -| `--outFilterMismatchNoverReadLmax` | `1.0` | Max ratio of mismatches to read length. | | `--outFilterScoreMin` | `0` | Min absolute alignment score. | | `--outFilterScoreMinOverLread` | `0.66` | Min alignment score normalized to read length. | | `--outFilterMatchNmin` | `0` | Min absolute matched bases. | @@ -104,8 +103,6 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--alignSJoverhangMin` | `5` | Min overhang for novel splice junctions. | | `--alignSJDBoverhangMin` | `3` | Min overhang for annotated junctions. | | `--alignSJstitchMismatchNmax` | `0 -1 0 0` | Max mismatches for SJ stitching `[noncan, GC/AG, AT/AC, noncan]`. | -| `--alignTranscriptsPerReadNmax` | `10000` | Max alignments kept per read before the score filters. | -| `--alignSoftClipAtReferenceEnds` | `Yes` | `No` prohibits soft-clipping past a chromosome end (Cufflinks compatibility). | ## Scoring penalties diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 4ac4b2a7..ec66a166 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -132,27 +132,6 @@ pub enum PairedAlignmentResult { }, } -/// Whether `t`'s soft-clipped ends would extend past the boundaries of its -/// chromosome: whether the clipped bases have nowhere to sit on the reference. -/// This is what `--alignSoftClipAtReferenceEnds No` prohibits. -fn clips_past_reference_end(t: &Transcript, index: &GenomeIndex, read_len: usize) -> bool { - let Some(first) = t.exons.first() else { - return false; - }; - let Some(last) = t.exons.last() else { - return false; - }; - let chr_start = index.genome.chr_start[t.chr_idx]; - let chr_end = chr_start + index.genome.chr_length[t.chr_idx]; - - // Bases clipped before the first aligned base, and after the last one. - let left_clip = first.read_start as u64; - let right_clip = (read_len - last.read_end) as u64; - - first.genome_start < chr_start.saturating_add(left_clip) && left_clip > 0 - || last.genome_end.saturating_add(right_clip) > chr_end && right_clip > 0 -} - /// Align a read to the genome. /// /// # Algorithm @@ -413,19 +392,6 @@ pub fn align_read( ); } - // Hard cap on how many alignments one read may carry into the filters - // (STAR's alignTranscriptsPerReadNmax). Transcripts are already ordered - // best-score-first, so the cap keeps the best ones. - if transcripts.len() > params.align_transcripts_per_read_nmax { - log::debug!( - "Read {}: {} alignments capped to alignTranscriptsPerReadNmax={}", - read_name, - transcripts.len(), - params.align_transcripts_per_read_nmax - ); - transcripts.truncate(params.align_transcripts_per_read_nmax); - } - // Score-range filter: keep only alignments within outFilterMultimapScoreRange of the best. // (STAR's multMapSelect step — must run before quality filters.) if !transcripts.is_empty() { @@ -444,21 +410,7 @@ pub fn align_read( let pre_filter_count = transcripts.len(); let mut filter_reasons = std::collections::HashMap::new(); - let prohibit_ref_end_clip = params - .align_soft_clip_at_reference_ends - .eq_ignore_ascii_case("No"); - transcripts.retain(|t| { - // A soft clip that would hang past the start or the end of the - // chromosome (STAR's alignSoftClipAtReferenceEnds No, needed for - // Cufflinks-compatible output). - if prohibit_ref_end_clip && clips_past_reference_end(t, index, read_seq.len()) { - *filter_reasons - .entry("soft_clip_at_reference_end") - .or_insert(0) += 1; - return false; - } - // Absolute score threshold if t.score < params.out_filter_score_min { *filter_reasons.entry("score_min").or_insert(0) += 1; @@ -485,25 +437,6 @@ pub fn align_read( return false; } - // Mismatch count over the *read* length (STAR's - // outFilterMismatchNoverReadLmax). STAR takes - // outFilterMismatchNoverLmax over the *mapped* length instead; the - // check above still divides by the read length, so the two behave - // identically here until issue #238 is fixed. - let read_mismatch_rate = t.n_mismatch as f64 / read_length; - if read_mismatch_rate > params.out_filter_mismatch_nover_read_lmax { - *filter_reasons.entry("mismatch_rate_read").or_insert(0) += 1; - log::debug!( - "Filtered {}: {:.1}% read-length mismatch rate > {:.1}% max ({}/{} bases)", - read_name, - read_mismatch_rate * 100.0, - params.out_filter_mismatch_nover_read_lmax * 100.0, - t.n_mismatch, - read_length - ); - return false; - } - // Relative mismatch count (mismatches / read_length) let mismatch_rate = t.n_mismatch as f64 / read_length; if mismatch_rate > params.out_filter_mismatch_nover_lmax { @@ -1542,8 +1475,6 @@ fn filter_paired_transcripts(paired_alns: &mut Vec, params: &Pa if combined_nm > params.out_filter_mismatch_nmax || (combined_nm as f64) > params.out_filter_mismatch_nover_lmax * (mate1_len + mate2_len) - || (combined_nm as f64) - > params.out_filter_mismatch_nover_read_lmax * (mate1_len + mate2_len) { paired_alns.clear(); return; diff --git a/src/params/mod.rs b/src/params/mod.rs index 58717d12..5b5fba10 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -873,12 +873,6 @@ pub struct Parameters { #[arg(long = "outFilterMismatchNoverLmax", default_value_t = 0.3)] pub out_filter_mismatch_nover_lmax: f64, - /// Max ratio of mismatches to *read* length (STAR's - /// `outFilterMismatchNoverReadLmax`). Applied to the full read length, - /// where `outFilterMismatchNoverLmax` is applied to the mapped length. - #[arg(long = "outFilterMismatchNoverReadLmax", default_value_t = 1.0)] - pub out_filter_mismatch_nover_read_lmax: f64, - /// Min alignment score (absolute) #[arg(long = "outFilterScoreMin", default_value_t = 0)] pub out_filter_score_min: i32, @@ -955,18 +949,6 @@ pub struct Parameters { #[arg(long = "alignEndsType", default_value = "Local")] pub align_ends_type: String, - /// Max number of alignments kept for one read before the score filters - /// run (STAR's `alignTranscriptsPerReadNmax`). A read that produces more - /// keeps only the highest-scoring `N`. - #[arg(long = "alignTranscriptsPerReadNmax", default_value_t = 10000)] - pub align_transcripts_per_read_nmax: usize, - - /// `Yes` (default) allows an alignment to soft-clip past the end of a - /// chromosome; `No` prohibits it, which is what Cufflinks-compatible - /// output needs. - #[arg(long = "alignSoftClipAtReferenceEnds", default_value = "Yes")] - pub align_soft_clip_at_reference_ends: String, - /// Min overlap (bases) between mates required to trigger merge-and-realign; 0 = off #[arg(long = "peOverlapNbasesMin", default_value_t = 0)] pub pe_overlap_nbases_min: u64, @@ -1869,35 +1851,6 @@ impl Parameters { } } - // alignSoftClipAtReferenceEnds is Yes/No; anything else would be - // silently read as Yes. - if !["Yes", "No"] - .iter() - .any(|v| v.eq_ignore_ascii_case(¶ms.align_soft_clip_at_reference_ends)) - { - return Err(command.error( - ErrorKind::InvalidValue, - format!( - "unsupported --alignSoftClipAtReferenceEnds '{}'; expected Yes or No", - params.align_soft_clip_at_reference_ends - ), - )); - } - - if params.align_transcripts_per_read_nmax == 0 { - return Err(command.error( - ErrorKind::InvalidValue, - "--alignTranscriptsPerReadNmax must be > 0", - )); - } - - if params.out_filter_mismatch_nover_read_lmax < 0.0 { - return Err(command.error( - ErrorKind::InvalidValue, - "--outFilterMismatchNoverReadLmax must be >= 0", - )); - } - // quantMode GeneCounts requires a GTF file if params.quant_gene_counts() && params.sjdb_gtf_file.is_none() { return Err(command.error( diff --git a/tests/cli_surface.rs b/tests/cli_surface.rs deleted file mode 100644 index e0e115b8..00000000 --- a/tests/cli_surface.rs +++ /dev/null @@ -1,278 +0,0 @@ -//! End-to-end checks for the STAR parameters closed in this change, so that -//! "accepted by the CLI" and "actually changes the output" stay different -//! claims. -//! -//! Genome: 20 kb of LCG(88888) background, the same generator the other -//! integration tests use, with one 100 bp segment planted at four extra -//! positions so a read from it multimaps. - -use assert_cmd::cargo::cargo_bin_cmd; -use std::fs; -use std::io::Write; -use std::path::{Path, PathBuf}; -use tempfile::TempDir; - -const READ_LEN: usize = 100; -/// Where the repeated 100 bp block is copied to. Four copies plus the -/// original gives a read from it five alignments. -const REPEAT_SRC: usize = 1_000; -const REPEAT_COPIES: [usize; 4] = [4_000, 8_000, 12_000, 16_000]; -/// A unique region used for the mismatch tests. -const UNIQUE_START: usize = 6_000; - -fn lcg_seq(seed: u32, length: usize) -> Vec { - let bases: [u8; 4] = *b"ACGT"; - let mut state = seed; - let mut seq = Vec::with_capacity(length); - for _ in 0..length { - state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); - seq.push(bases[((state >> 16) & 3) as usize]); - } - seq -} - -fn build_genome() -> Vec { - let mut genome = lcg_seq(88888, 20_000); - let block: Vec = genome[REPEAT_SRC..REPEAT_SRC + READ_LEN].to_vec(); - for &dst in &REPEAT_COPIES { - genome[dst..dst + READ_LEN].copy_from_slice(&block); - } - genome -} - -fn write_fasta(dir: &Path, genome: &[u8]) -> PathBuf { - let path = dir.join("genome.fa"); - let mut f = fs::File::create(&path).unwrap(); - writeln!(f, ">chr1").unwrap(); - f.write_all(genome).unwrap(); - writeln!(f).unwrap(); - path -} - -fn build_index(fasta: &Path, genome_dir: &Path) { - fs::create_dir_all(genome_dir).unwrap(); - cargo_bin_cmd!("rustar-aligner") - .args([ - "--runMode", - "genomeGenerate", - "--genomeDir", - genome_dir.to_str().unwrap(), - "--genomeFastaFiles", - fasta.to_str().unwrap(), - "--genomeSAindexNbases", - "7", - ]) - .assert() - .success(); -} - -fn write_fastq(path: &Path, reads: &[(String, Vec)]) { - let mut f = fs::File::create(path).unwrap(); - for (name, seq) in reads { - writeln!(f, "@{name}").unwrap(); - f.write_all(seq).unwrap(); - writeln!(f).unwrap(); - writeln!(f, "+").unwrap(); - writeln!(f, "{}", "I".repeat(seq.len())).unwrap(); - } -} - -/// Run the aligner and return the SAM records (header lines dropped). -fn align(genome_dir: &Path, fastq: &Path, prefix: &str, extra: &[&str]) -> Vec { - let mut cmd = cargo_bin_cmd!("rustar-aligner"); - cmd.args([ - "--runMode", - "alignReads", - "--genomeDir", - genome_dir.to_str().unwrap(), - "--readFilesIn", - fastq.to_str().unwrap(), - "--outFileNamePrefix", - prefix, - ]); - cmd.args(extra); - cmd.assert().success(); - - fs::read_to_string(format!("{prefix}Aligned.out.sam")) - .unwrap() - .lines() - .filter(|l| !l.starts_with('@')) - .map(str::to_string) - .collect() -} - -fn mismatched_read(genome: &[u8], n_mismatches: usize) -> Vec { - let mut seq = genome[UNIQUE_START..UNIQUE_START + READ_LEN].to_vec(); - // Spread the substitutions out so they cannot all be soft-clipped off one - // end instead of being counted as mismatches. - for i in 0..n_mismatches { - let pos = 15 + i * 20; - seq[pos] = match seq[pos] { - b'A' => b'C', - b'C' => b'G', - b'G' => b'T', - _ => b'A', - }; - } - seq -} - -// ── outFilterMismatchNoverReadLmax ────────────────────────────────────────── - -#[test] -fn mismatch_nover_read_lmax_filters_by_read_length_ratio() { - let dir = TempDir::new().unwrap(); - let root = dir.path(); - let genome = build_genome(); - let genome_dir = root.join("genome"); - build_index(&write_fasta(root, &genome), &genome_dir); - - let fq = root.join("mm.fq"); - write_fastq(&fq, &[("mm3".to_string(), mismatched_read(&genome, 3))]); - - // Control: the read maps with the default ratio of 1.0. - let permissive = align( - &genome_dir, - &fq, - &format!("{}/permissive_", root.display()), - &["--outFilterMismatchNoverReadLmax", "1.0"], - ); - assert_eq!(permissive.len(), 1, "read should map by default"); - let control_flag: u32 = permissive[0].split('\t').nth(1).unwrap().parse().unwrap(); - assert_eq!( - control_flag & 0x4, - 0, - "control read must be mapped: {}", - permissive[0] - ); - - // 0.01 of 100 bases allows one mismatch; the read carries three. - let strict = align( - &genome_dir, - &fq, - &format!("{}/strict_", root.display()), - &["--outFilterMismatchNoverReadLmax", "0.01"], - ); - let mapped = strict - .iter() - .filter(|r| { - let flag: u32 = r.split('\t').nth(1).unwrap().parse().unwrap(); - flag & 0x4 == 0 - }) - .count(); - assert_eq!( - mapped, 0, - "a 3-mismatch read must fail a 0.01 read-length ratio: {strict:?}" - ); -} - -// ── alignTranscriptsPerReadNmax ───────────────────────────────────────────── - -#[test] -fn transcripts_per_read_nmax_caps_the_alignments_kept() { - let dir = TempDir::new().unwrap(); - let root = dir.path(); - let genome = build_genome(); - let genome_dir = root.join("genome"); - build_index(&write_fasta(root, &genome), &genome_dir); - - let fq = root.join("multi.fq"); - let repeat = genome[REPEAT_SRC..REPEAT_SRC + READ_LEN].to_vec(); - write_fastq(&fq, &[("multi".to_string(), repeat)]); - - // Control: all five copies are reported when nothing caps them. - let uncapped = align( - &genome_dir, - &fq, - &format!("{}/uncapped_", root.display()), - &["--outFilterMultimapNmax", "20"], - ); - assert!( - uncapped.len() >= 3, - "the planted repeat should multimap, got {} record(s)", - uncapped.len() - ); - - let capped = align( - &genome_dir, - &fq, - &format!("{}/capped_", root.display()), - &[ - "--outFilterMultimapNmax", - "20", - "--alignTranscriptsPerReadNmax", - "2", - ], - ); - assert!( - capped.len() < uncapped.len(), - "the cap must reduce the alignments kept: {} capped vs {} uncapped", - capped.len(), - uncapped.len() - ); - assert!( - capped.len() <= 2, - "at most 2 alignments may survive --alignTranscriptsPerReadNmax 2, got {}", - capped.len() - ); -} - -// ── alignSoftClipAtReferenceEnds ──────────────────────────────────────────── - -#[test] -fn soft_clip_at_reference_ends_no_prohibits_clipping_past_the_chromosome() { - let dir = TempDir::new().unwrap(); - let root = dir.path(); - let genome = build_genome(); - let genome_dir = root.join("genome"); - build_index(&write_fasta(root, &genome), &genome_dir); - - // A read whose last 20 bases fall past the end of chr1: 80 bases of the - // genome tail, then 20 bases that exist nowhere, so the aligner has to - // soft-clip them past the reference end. - let mut seq = genome[genome.len() - 80..].to_vec(); - seq.extend_from_slice(b"ACACACACACGTGTGTGTGT"); - let fq = root.join("edge.fq"); - write_fastq(&fq, &[("edge".to_string(), seq)]); - - let allowed = align( - &genome_dir, - &fq, - &format!("{}/allowed_", root.display()), - &["--alignSoftClipAtReferenceEnds", "Yes"], - ); - let mapped_allowed = allowed - .iter() - .filter(|r| { - let flag: u32 = r.split('\t').nth(1).unwrap().parse().unwrap(); - flag & 0x4 == 0 - }) - .count(); - assert_eq!( - mapped_allowed, 1, - "control: the read maps with a soft clip past the end: {allowed:?}" - ); - assert!( - allowed[0].split('\t').nth(5).unwrap().contains('S'), - "control alignment should carry a soft clip: {}", - allowed[0] - ); - - let prohibited = align( - &genome_dir, - &fq, - &format!("{}/prohibited_", root.display()), - &["--alignSoftClipAtReferenceEnds", "No"], - ); - let mapped_prohibited = prohibited - .iter() - .filter(|r| { - let flag: u32 = r.split('\t').nth(1).unwrap().parse().unwrap(); - flag & 0x4 == 0 - }) - .count(); - assert_eq!( - mapped_prohibited, 0, - "--alignSoftClipAtReferenceEnds No must reject it: {prohibited:?}" - ); -} diff --git a/tests/parameter_surface.rs b/tests/parameter_surface.rs index dd2a6bc4..bd707ba4 100644 --- a/tests/parameter_surface.rs +++ b/tests/parameter_surface.rs @@ -84,12 +84,16 @@ const ACCEPTED_BUT_INERT: &[(&str, &str)] = &[ /// Adding a name here must always be a deliberate act. Removing one is what /// progress looks like. const NOT_YET_ACCEPTED: &[&str] = &[ - // Aligner core (annotated-junction stitching, alignEndsType, in-recursion - // length penalty). + // Aligner core; implemented by the aligner-core PR, not here. "alignEndsProtrude", "alignInsertionFlush", + "alignSoftClipAtReferenceEnds", + "alignTranscriptsPerReadNmax", + "outFilterMismatchNoverReadLmax", "seedNoneLociPerWindow", "seedSplitMin", + // Aligner core (annotated-junction stitching, alignEndsType, in-recursion + // length penalty). // Long reads. "winReadCoverageBasesMin", // Chimeric multimapping.