cli parity #147 with fixes

CHANGELOG.md

@@ -21,6 +21,30 @@ Sections commonly used: Features, Bug fixes, Other changes.
 
 ### Features
 
+- **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.)
+
+  - Implemented: `--outSAMmode` (`Full`/`NoQS`/`None`), `--outSJtype None`,
+    `--outSJfilterReads Unique`, `--outSAMheaderHD`, `--outSAMheaderPG`,
+    `--outSAMheaderCommentFile`, `--readFilesPrefix`, `--readNameSeparator`,
+    `--readQualityScoreBase`, `--outQSconversionAdd`.
+  - Refused loudly rather than accepted and ignored:
+    `--outSAMfilter KeepOnlyAddedReferences` / `KeepAllAddedReferences` and
+    `--readFilesType SAM`, both of which need machinery this aligner does not
+    have.
+  - Accepted and documented as output-neutral: the `--limit*` family,
+    `--outTmpDir`, `--outTmpKeep`, `--runDirPerm`, `--genomeFileSizes`,
+    `--outBAMsorting*`, `--readMatesLengthsIn`.
+
+  `tests/parameter_surface.rs` now checks every STAR 2.7.11b parameter name
+  against the CLI, so the coverage figure is machine-checked and surface drift
+  fails a test.
+
+### Bug fixes
+
+- Read names are cut at `--readNameSeparator` (default `/`), as STAR does. A
+  read named `foo/1` was previously emitted as `foo/1` where STAR emits `foo`.
+
 - **STARsolo single-cell quantification (`--soloType`)** — the 10x
   Chromium / plate-based count-matrix pipeline, ported from STAR and
   verified against real STARsolo (#90).

src/io/fastq.rs

@@ -72,6 +72,19 @@ pub struct PairedRead {
 /// FASTQ reader that handles decompression and base encoding
 pub struct FastqReader {
     inner: fastq::io::Reader<Box<dyn BufRead + Send>>,
+    /// Signed shift applied to every input quality byte so the rest of the
+    /// pipeline always sees Phred+33.
+    ///
+    /// `--readQualityScoreBase 64` contributes `-31`, converting Solexa/Illumina
+    /// 1.3 input. `--outQSconversionAdd` contributes its own value, which STAR
+    /// documents as an output conversion; applying it here rather than at the
+    /// writer means it also reaches anything else that reads qualities. That is
+    /// only observable when both it and STARsolo are in use, which STAR itself
+    /// does not support either.
+    qual_shift: i32,
+    /// Characters that terminate a read name (`--readNameSeparator`). Empty
+    /// means keep the whole name.
+    name_separators: Vec<u8>,
 }
 
 impl FastqReader {
@@ -115,9 +128,30 @@ impl FastqReader {
 
         Ok(Self {
             inner: fastq_reader,
+            qual_shift: 0,
+            name_separators: vec![b'/'],
         })
     }
 
+    /// Apply the read-input knobs: `--readQualityScoreBase`,
+    /// `--outQSconversionAdd` and `--readNameSeparator`.
+    #[must_use]
+    pub fn with_params(mut self, params: &crate::params::Parameters) -> Self {
+        let base = if params.read_quality_score_base == 0 {
+            33
+        } else {
+            params.read_quality_score_base
+        };
+        self.qual_shift = (33 - base) + params.out_qs_conversion_add;
+        self.name_separators = params
+            .read_name_separator
+            .iter()
+            .filter(|s| s.as_str() != "-")
+            .filter_map(|s| s.as_bytes().first().copied())
+            .collect();
+        self
+    }
+
     /// Open FASTQ file using external decompression command
     fn open_with_command(path: &Path, cmd: &str) -> Result<Box<dyn BufRead + Send>, Error> {
         let mut child = Command::new(cmd)
@@ -150,7 +184,25 @@ impl FastqReader {
 
                 let sequence = record.sequence().iter().map(|&b| encode_base(b)).collect();
 
-                let quality = record.quality_scores().to_vec();
+                let name = match self
+                    .name_separators
+                    .iter()
+                    .filter_map(|&sep| name.as_bytes().iter().position(|&b| b == sep))
+                    .min()
+                {
+                    Some(cut) => name[..cut].to_string(),
+                    None => name,
+                };
+
+                let quality = if self.qual_shift == 0 {
+                    record.quality_scores().to_vec()
+                } else {
+                    record
+                        .quality_scores()
+                        .iter()
+                        .map(|&b| (b as i32 + self.qual_shift).clamp(33, 126) as u8)
+                        .collect()
+                };
 
                 Ok(Some(EncodedRead {
                     name,
@@ -205,6 +257,15 @@ impl PairedFastqReader {
         Ok(Self { reader1, reader2 })
     }
 
+    /// Apply the read-input knobs to both mates. See
+    /// [`FastqReader::with_params`].
+    #[must_use]
+    pub fn with_params(mut self, params: &crate::params::Parameters) -> Self {
+        self.reader1 = self.reader1.with_params(params);
+        self.reader2 = self.reader2.with_params(params);
+        self
+    }
+
     /// Get next paired read with name validation
     ///
     /// # Returns
@@ -565,9 +626,11 @@ mod tests {
 
         let pair1 = reader.next_paired().unwrap().unwrap();
         assert_eq!(pair1.name, "read1");
-        assert_eq!(pair1.mate1.name, "read1/1");
+        // Read names are cut at `/`, STAR's default --readNameSeparator, so
+        // both mates report the shared name rather than the /1 and /2 forms.
+        assert_eq!(pair1.mate1.name, "read1");
         assert_eq!(pair1.mate1.sequence, vec![0, 1, 2, 3]); // ACGT
-        assert_eq!(pair1.mate2.name, "read1/2");
+        assert_eq!(pair1.mate2.name, "read1");
         assert_eq!(pair1.mate2.sequence, vec![2, 2, 1, 1]); // GGCC
 
         let pair2 = reader.next_paired().unwrap().unwrap();

src/io/sam.rs

@@ -950,8 +950,35 @@ where
 {
     let mut builder = sam::Header::builder();
 
-    // @HD line (default version and unsorted)
-    builder = builder.set_header(Map::default());
+    // @HD line. `--outSAMheaderHD` replaces it wholesale, given as the
+    // tab-separated fields STAR expects (`@HD VN:1.4 SO:coordinate`).
+    if params.out_sam_header_hd.is_empty() {
+        builder = builder.set_header(Map::default());
+    } else {
+        let mut hd = Map::<sam::header::record::value::map::Header>::default();
+        for field in &params.out_sam_header_hd {
+            let field = field.trim();
+            // STAR takes the leading `@HD` as part of the value list; ignore it.
+            if field.is_empty() || field == "@HD" {
+                continue;
+            }
+            let Some((tag, value)) = field.split_once(':') else {
+                return Err(Error::Parameter(format!(
+                    "--outSAMheaderHD field '{field}' is not TAG:value"
+                )));
+            };
+            if tag.len() != 2 {
+                return Err(Error::Parameter(format!(
+                    "--outSAMheaderHD tag '{tag}' is not two characters"
+                )));
+            }
+            let tag_bytes: [u8; 2] = tag.as_bytes()[..2].try_into().unwrap();
+            let other_tag: HeaderOtherTag<_> = HeaderOtherTag::try_from(tag_bytes)
+                .map_err(|e| Error::Parameter(format!("invalid @HD tag '{tag}': {e}")))?;
+            hd.other_fields_mut().insert(other_tag, value.into());
+        }
+        builder = builder.set_header(hd);
+    }
 
     // @SQ lines for each reference
     for (name, length) in refs {
@@ -1009,6 +1036,51 @@ where
         .insert(program_tag::COMMAND_LINE, BString::from(cl));
     builder = builder.add_program("rustar-aligner", pg);
 
+    // `--outSAMheaderPG` adds one further @PG line, given the same way.
+    if !params.out_sam_header_pg.is_empty() {
+        let mut extra = Map::<Program>::default();
+        let mut id: Option<String> = None;
+        for field in &params.out_sam_header_pg {
+            let field = field.trim();
+            if field.is_empty() || field == "@PG" {
+                continue;
+            }
+            let Some((tag, value)) = field.split_once(':') else {
+                return Err(Error::Parameter(format!(
+                    "--outSAMheaderPG field '{field}' is not TAG:value"
+                )));
+            };
+            if tag == "ID" {
+                id = Some(value.to_string());
+                continue;
+            }
+            if tag.len() != 2 {
+                return Err(Error::Parameter(format!(
+                    "--outSAMheaderPG tag '{tag}' is not two characters"
+                )));
+            }
+            let tag_bytes: [u8; 2] = tag.as_bytes()[..2].try_into().unwrap();
+            let other_tag: HeaderOtherTag<_> = HeaderOtherTag::try_from(tag_bytes)
+                .map_err(|e| Error::Parameter(format!("invalid @PG tag '{tag}': {e}")))?;
+            extra.other_fields_mut().insert(other_tag, value.into());
+        }
+        let id =
+            id.ok_or_else(|| Error::Parameter("--outSAMheaderPG needs an ID: field".to_string()))?;
+        builder = builder.add_program(id, extra);
+    }
+
+    // `--outSAMheaderCommentFile` contributes one @CO line per line of the
+    // named file. `-` (the default) means no comments.
+    if params.out_sam_header_comment_file != "-" {
+        let path = std::path::Path::new(&params.out_sam_header_comment_file);
+        let contents = std::fs::read_to_string(path).map_err(|e| Error::io(e, path))?;
+        for line in contents.lines() {
+            if !line.is_empty() {
+                builder = builder.add_comment(line);
+            }
+        }
+    }
+
     Ok(builder.build())
 }
 

src/junction/sj_output.rs

@@ -48,6 +48,12 @@ pub(crate) struct SjCounts {
 pub struct SpliceJunctionStats {
     /// Thread-safe map for parallel accumulation
     junctions: DashMap<SjKey, SjCounts>,
+    /// `--outSJfilterReads Unique`: drop junctions from multi-mapping reads at
+    /// record time rather than discounting them at output time. STAR gates the
+    /// recording itself (`ReadAlign::recordSJ`: the whole read is skipped
+    /// unless `outSJfilterReads=="All" || nTrO==1`), so a multimapper
+    /// contributes nothing at all — not its counts, and not its overhang.
+    unique_reads_only: bool,
 }
 
 impl Clone for SpliceJunctionStats {
@@ -66,7 +72,10 @@ impl Clone for SpliceJunctionStats {
                 },
             );
         }
-        Self { junctions: new_map }
+        Self {
+            junctions: new_map,
+            unique_reads_only: self.unique_reads_only,
+        }
     }
 }
 
@@ -75,6 +84,15 @@ impl SpliceJunctionStats {
     pub fn new() -> Self {
         Self {
             junctions: DashMap::new(),
+            unique_reads_only: false,
+        }
+    }
+
+    /// [`Self::new`] honouring `--outSJfilterReads`.
+    pub fn with_params(params: &Parameters) -> Self {
+        Self {
+            junctions: DashMap::new(),
+            unique_reads_only: params.out_sj_filter_reads == "Unique",
         }
     }
 
@@ -101,6 +119,10 @@ impl SpliceJunctionStats {
         overhang: u32,
         annotated: bool,
     ) {
+        if self.unique_reads_only && !is_unique {
+            return; // STAR records nothing for this read under `Unique`.
+        }
+
         let key = SjKey {
             chr_idx,
             intron_start: start,
@@ -132,11 +154,12 @@ impl SpliceJunctionStats {
             .map(|entry| {
                 let key = entry.key().clone();
                 let counts = entry.value();
+                let multi = counts.multi_count.load(Ordering::Relaxed);
                 (
                     key,
                     counts.annotated,
                     counts.unique_count.load(Ordering::Relaxed),
-                    counts.multi_count.load(Ordering::Relaxed),
+                    multi,
                     counts.max_overhang.load(Ordering::Relaxed),
                 )
             })
@@ -281,11 +304,12 @@ impl SpliceJunctionStats {
             .map(|entry| {
                 let key = entry.key().clone();
                 let counts = entry.value();
+                let multi = counts.multi_count.load(Ordering::Relaxed);
                 (
                     key,
                     counts.annotated,
                     counts.unique_count.load(Ordering::Relaxed),
-                    counts.multi_count.load(Ordering::Relaxed),
+                    multi,
                     counts.max_overhang.load(Ordering::Relaxed),
                 )
             })