From 6740954094071e2713edb9eeb07e68158aeaab26 Mon Sep 17 00:00:00 2001 From: Phil Ewels Date: Fri, 17 Apr 2026 14:21:59 +0200 Subject: [PATCH 1/7] Add --outSAMattrRGline flag and emit @RG header + RG:Z tags Parses STAR's comma-separated multi-RG syntax (first field must be ID: on each block, `,` separates blocks, tokens are tab-joined into the SAM @RG body). Replicates a single RG line across N input files like STAR does, and errors on count mismatch or RG-in-outSAMattributes without a line. Auto-appends RG to sam_attribute_set() when a line is set (matches Parameters_samAttributes.cpp:201). SAM/BAM headers now emit one @RG per unique ID, and every record (mapped, unmapped, half-mapped, both- mates) carries RG:Z: when attrs contain RG. Unblocks nf-core/rnaseq, which invokes STAR with `--outSAMattrRGline 'ID:X' 'SM:X'` and previously errored on the unknown flag. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/io/bam.rs | 9 ++- src/io/sam.rs | 182 ++++++++++++++++++++++++++++++++++++++++++++--- src/lib.rs | 7 ++ src/params.rs | 192 +++++++++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 378 insertions(+), 12 deletions(-) diff --git a/src/io/bam.rs b/src/io/bam.rs index aa70cdc9..abfb2f41 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -125,9 +125,10 @@ mod tests { let read_qual = vec![30, 30, 30, 30]; // Build record using SAM builder - let record = - crate::io::sam::SamWriter::build_unmapped_record(read_name, &read_seq, &read_qual) - .unwrap(); + let record = crate::io::sam::SamWriter::build_unmapped_record( + read_name, &read_seq, &read_qual, None, + ) + .unwrap(); let result = writer.write_batch(&[record]); assert!(result.is_ok(), "Writing unmapped read should succeed"); @@ -204,12 +205,14 @@ mod tests { "read1", &[0, 1, 2, 3], &[30, 30, 30, 30], + None, ) .unwrap(), crate::io::sam::SamWriter::build_unmapped_record( "read2", &[0, 1, 2, 3], &[30, 30, 30, 30], + None, ) .unwrap(), ]; diff --git a/src/io/sam.rs b/src/io/sam.rs index 35ca5ac2..551cd945 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -15,7 +15,10 @@ use noodles::sam::alignment::record::data::field::Tag; use noodles::sam::alignment::record_buf::data::field::Value; use noodles::sam::alignment::record_buf::data::field::value::Array; use noodles::sam::alignment::record_buf::{QualityScores, RecordBuf, Sequence}; -use noodles::sam::header::record::value::{Map, map::Program}; +use noodles::sam::header::record::value::{ + Map, + map::{Program, ReadGroup, tag::Other as HeaderOtherTag}, +}; use std::collections::HashSet; use std::fmt::Write as FmtWrite; use std::fs::File; @@ -104,9 +107,11 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + let rg_ids = params.rg_ids()?; + let rg_id = rg_ids.first().map(String::as_str); for (hit_index, transcript) in transcripts.iter().take(max_output).enumerate() { - let record = transcript_to_record( + let mut record = transcript_to_record( transcript, read_name, read_seq, @@ -117,6 +122,7 @@ impl SamWriter { hit_index + 1, // 1-based &attrs, )?; + maybe_insert_rg_tag(&mut record, &attrs, rg_id); self.writer.write_alignment_record(&self.header, &record)?; } @@ -136,7 +142,7 @@ impl SamWriter { read_seq: &[u8], read_qual: &[u8], ) -> Result<(), Error> { - let record = Self::build_unmapped_record(read_name, read_seq, read_qual)?; + let record = Self::build_unmapped_record(read_name, read_seq, read_qual, None)?; self.writer.write_alignment_record(&self.header, &record)?; Ok(()) } @@ -185,10 +191,12 @@ impl SamWriter { /// * `read_name` - Read identifier /// * `read_seq` - Read sequence (encoded) /// * `read_qual` - Quality scores + /// * `rg_id` - Optional read group ID to emit as `RG:Z:` tag pub fn build_unmapped_record( read_name: &str, read_seq: &[u8], read_qual: &[u8], + rg_id: Option<&str>, ) -> Result { let mut record = RecordBuf::default(); @@ -206,6 +214,12 @@ impl SamWriter { // Quality scores *record.quality_scores_mut() = QualityScores::from(read_qual.to_vec()); + if let Some(id) = rg_id { + record + .data_mut() + .insert(Tag::READ_GROUP, Value::String(BString::from(id))); + } + Ok(record) } @@ -243,10 +257,12 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + let rg_ids = params.rg_ids()?; + let rg_id = rg_ids.first().map(String::as_str); let mut records = Vec::with_capacity(max_output); for (hit_index, transcript) in transcripts.iter().take(max_output).enumerate() { - let record = transcript_to_record( + let mut record = transcript_to_record( transcript, read_name, read_seq, @@ -257,6 +273,7 @@ impl SamWriter { hit_index + 1, // 1-based &attrs, )?; + maybe_insert_rg_tag(&mut record, &attrs, rg_id); records.push(record); } @@ -291,7 +308,7 @@ impl SamWriter { if paired_alignments.is_empty() { // Both mates unmapped return Self::build_paired_unmapped_records( - read_name, mate1_seq, mate1_qual, mate2_seq, mate2_qual, + read_name, mate1_seq, mate1_qual, mate2_seq, mate2_qual, params, ); } @@ -307,6 +324,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + let rg_ids = params.rg_ids()?; + let rg_id = rg_ids.first().map(String::as_str); let mut records = Vec::with_capacity(max_output * 2); @@ -317,7 +336,7 @@ impl SamWriter { let combined_score = paired_aln.combined_wt_score; // Create record for mate1 (this=mate1, mate=mate2) - let rec1 = build_paired_mate_record( + let mut rec1 = build_paired_mate_record( read_name, mate1_seq, mate1_qual, @@ -333,10 +352,11 @@ impl SamWriter { combined_score, &attrs, )?; + maybe_insert_rg_tag(&mut rec1, &attrs, rg_id); records.push(rec1); // Create record for mate2 (this=mate2, mate=mate1) - let rec2 = build_paired_mate_record( + let mut rec2 = build_paired_mate_record( read_name, mate2_seq, mate2_qual, @@ -352,6 +372,7 @@ impl SamWriter { combined_score, &attrs, )?; + maybe_insert_rg_tag(&mut rec2, &attrs, rg_id); records.push(rec2); } @@ -389,6 +410,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + let rg_ids = params.rg_ids()?; + let rg_id = rg_ids.first().map(String::as_str); // Compute mapped mate's per-chr position for co-location let chr_start = genome.chr_start[mapped_transcript.chr_idx]; @@ -491,6 +514,7 @@ impl SamWriter { ); data.insert(Tag::new(b'M', b'D'), Value::String(BString::from(md))); } + maybe_insert_rg_tag(&mut mapped_rec, &attrs, rg_id); // --- Build unmapped mate record --- let mut unmapped_rec = RecordBuf::default(); @@ -530,6 +554,7 @@ impl SamWriter { let unmapped_seq_bytes: Vec = unmapped_seq.iter().map(|&b| decode_base(b)).collect(); *unmapped_rec.sequence_mut() = Sequence::from(unmapped_seq_bytes); *unmapped_rec.quality_scores_mut() = QualityScores::from(unmapped_qual.to_vec()); + maybe_insert_rg_tag(&mut unmapped_rec, &attrs, rg_id); // Order: mate1 first, mate2 second if mate1_is_mapped { @@ -550,8 +575,12 @@ impl SamWriter { mate1_qual: &[u8], mate2_seq: &[u8], mate2_qual: &[u8], + params: &Parameters, ) -> Result, Error> { let mut records = Vec::with_capacity(2); + let attrs = params.sam_attribute_set(); + let rg_ids = params.rg_ids()?; + let rg_id = rg_ids.first().map(String::as_str); // Mate1 record let mut rec1 = RecordBuf::default(); @@ -567,6 +596,7 @@ impl SamWriter { let seq1_bytes: Vec = mate1_seq.iter().map(|&b| decode_base(b)).collect(); *rec1.sequence_mut() = Sequence::from(seq1_bytes); *rec1.quality_scores_mut() = QualityScores::from(mate1_qual.to_vec()); + maybe_insert_rg_tag(&mut rec1, &attrs, rg_id); records.push(rec1); // Mate2 record @@ -583,6 +613,7 @@ impl SamWriter { let seq2_bytes: Vec = mate2_seq.iter().map(|&b| decode_base(b)).collect(); *rec2.sequence_mut() = Sequence::from(seq2_bytes); *rec2.quality_scores_mut() = QualityScores::from(mate2_qual.to_vec()); + maybe_insert_rg_tag(&mut rec2, &attrs, rg_id); records.push(rec2); Ok(records) @@ -590,7 +621,7 @@ impl SamWriter { } /// Build paired SAM header from genome -pub fn build_sam_header(genome: &Genome, _params: &Parameters) -> Result { +pub fn build_sam_header(genome: &Genome, params: &Parameters) -> Result { let mut builder = sam::Header::builder(); // @HD line (default version and unsorted) @@ -610,12 +641,52 @@ pub fn build_sam_header(genome: &Genome, _params: &Parameters) -> Result = HashSet::new(); + for line in &rg_lines { + let mut fields = line.split('\t'); + let id = fields + .next() + .and_then(|f| f.strip_prefix("ID:")) + .ok_or_else(|| Error::Parameter(format!("malformed RG line '{}'", line)))?; + if !seen_ids.insert(id.to_string()) { + continue; + } + let mut map = Map::::default(); + for field in fields { + if field.len() < 3 || &field[2..3] != ":" { + return Err(Error::Parameter(format!( + "RG field '{}' is not TAG:value", + field + ))); + } + let tag_bytes: [u8; 2] = field.as_bytes()[..2].try_into().unwrap(); + let other_tag: HeaderOtherTag<_> = + HeaderOtherTag::try_from(tag_bytes).map_err(|e| { + Error::Parameter(format!("invalid RG tag '{}': {}", &field[..2], e)) + })?; + map.other_fields_mut().insert(other_tag, field[3..].into()); + } + builder = builder.add_read_group(id, map); + } + // @PG line builder = builder.add_program("ruSTAR", Map::::default()); Ok(builder.build()) } +/// Insert `RG:Z:` on the record when RG is in `attrs` and an ID is set. +fn maybe_insert_rg_tag(record: &mut RecordBuf, attrs: &HashSet, rg_id: Option<&str>) { + if let (true, Some(id)) = (attrs.contains("RG"), rg_id) { + record + .data_mut() + .insert(Tag::READ_GROUP, Value::String(BString::from(id))); + } +} + /// Convert Transcript to SAM record #[allow(clippy::too_many_arguments)] fn transcript_to_record( @@ -1129,6 +1200,99 @@ mod tests { assert_eq!(header.reference_sequences().len(), 1); } + #[test] + fn test_build_sam_header_with_rg() { + let genome = make_test_genome(); + let params = Parameters::parse_from(vec![ + "ruSTAR", + "--readFilesIn", + "test.fq", + "--outSAMattrRGline", + "ID:rg0", + "SM:sample0", + "LB:lib0", + ]); + let header = build_sam_header(&genome, ¶ms).unwrap(); + let rgs = header.read_groups(); + assert_eq!(rgs.len(), 1); + assert!(rgs.contains_key(&b"rg0"[..])); + let map = rgs.get(&b"rg0"[..]).unwrap(); + // SM and LB should be present as other_fields + let sm_tag = HeaderOtherTag::<_>::try_from([b'S', b'M']).unwrap(); + let lb_tag = HeaderOtherTag::<_>::try_from([b'L', b'B']).unwrap(); + let sm: &[u8] = map.other_fields().get(&sm_tag).unwrap().as_ref(); + let lb: &[u8] = map.other_fields().get(&lb_tag).unwrap().as_ref(); + assert_eq!(sm, b"sample0"); + assert_eq!(lb, b"lib0"); + } + + #[test] + fn test_sam_output_includes_rg_header_and_tag() { + use crate::align::transcript::Exon; + use std::io::Read; + + let genome = make_test_genome(); + let params = Parameters::parse_from(vec![ + "ruSTAR", + "--readFilesIn", + "test.fq", + "--outSAMattrRGline", + "ID:rg0", + "SM:sample0", + ]); + + let tmpfile = NamedTempFile::new().unwrap(); + let mut writer = SamWriter::create(tmpfile.path(), &genome, ¶ms).unwrap(); + + let transcript = Transcript { + chr_idx: 0, + genome_start: 0, + genome_end: 4, + is_reverse: false, + exons: vec![Exon { + genome_start: 0, + genome_end: 4, + read_start: 0, + read_end: 4, + }], + cigar: vec![CigarOp::Match(4)], + score: 100, + n_mismatch: 0, + n_gap: 0, + n_junction: 0, + junction_motifs: vec![], + junction_annotated: vec![], + read_seq: vec![0, 1, 2, 3], + }; + + writer + .write_alignment( + "read1", + &[0, 1, 2, 3], + &[30, 30, 30, 30], + &[transcript], + &genome, + ¶ms, + 1, + ) + .unwrap(); + drop(writer); + + let mut contents = String::new(); + std::fs::File::open(tmpfile.path()) + .unwrap() + .read_to_string(&mut contents) + .unwrap(); + assert!( + contents.contains("@RG\tID:rg0\tSM:sample0"), + "missing @RG header; got:\n{contents}" + ); + assert!( + contents.contains("RG:Z:rg0"), + "missing RG:Z tag on record; got:\n{contents}" + ); + } + #[test] fn test_sam_writer_creation() { let genome = make_test_genome(); @@ -1207,6 +1371,7 @@ mod tests { let mate1_qual = vec![30, 30, 30, 30]; let mate2_seq = vec![3, 2, 1, 0]; // TGCA let mate2_qual = vec![30, 30, 30, 30]; + let params = Parameters::parse_from(vec!["ruSTAR", "--readFilesIn", "t.fq"]); let records = SamWriter::build_paired_unmapped_records( "read1", @@ -1214,6 +1379,7 @@ mod tests { &mate1_qual, &mate2_seq, &mate2_qual, + ¶ms, ) .unwrap(); diff --git a/src/lib.rs b/src/lib.rs index 6c930dd8..e67f8c92 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -472,6 +472,8 @@ fn align_reads_single_end( let max_multimaps = params.out_filter_multimap_nmax as usize; let output_unmapped = params.out_sam_unmapped != params::OutSamUnmapped::None; let by_sjout = params.out_filter_type == OutFilterType::BySJout; + let rg_ids = params.rg_ids()?; + let rg_id_owned = rg_ids.first().cloned(); // Buffer for BySJout mode: accumulate all results before filtering let mut bysj_buffer: Vec = Vec::new(); @@ -531,6 +533,7 @@ fn align_reads_single_end( &read.name, &clipped_seq, &clipped_qual, + rg_id_owned.as_deref(), )?; buffer.push(record); } @@ -597,6 +600,7 @@ fn align_reads_single_end( &read.name, &clipped_seq, &clipped_qual, + rg_id_owned.as_deref(), )?; buffer.push(record); } @@ -777,6 +781,7 @@ fn align_reads_paired_end( let max_multimaps = params.out_filter_multimap_nmax as usize; let output_unmapped = params.out_sam_unmapped != params::OutSamUnmapped::None; let by_sjout = params.out_filter_type == OutFilterType::BySJout; + let _ = params.rg_ids()?; // Validate early; per-call uses params inside builders. // Buffer for BySJout mode let mut bysj_buffer: Vec = Vec::new(); @@ -847,6 +852,7 @@ fn align_reads_paired_end( &m1_qual, &m2_seq, &m2_qual, + params, )?; for record in records { buffer.push(record); @@ -984,6 +990,7 @@ fn align_reads_paired_end( &m1_qual, &m2_seq, &m2_qual, + params, )?; for record in records { buffer.push(record); diff --git a/src/params.rs b/src/params.rs index f9f94553..2a9db479 100644 --- a/src/params.rs +++ b/src/params.rs @@ -302,6 +302,12 @@ pub struct Parameters { #[arg(long = "outSAMattributes", num_args = 1.., default_values_t = vec!["Standard".to_string()])] pub out_sam_attributes: Vec, + /// Read group line(s) for the `@RG` SAM header. Space-separated fields; + /// a bare `,` separates multiple RG blocks. Each block must start with `ID:`. + /// Default `-` means no `@RG` line (matches STAR). + #[arg(long = "outSAMattrRGline", num_args = 1.., default_values_t = vec!["-".to_string()])] + pub out_sam_attr_rg_line: Vec, + /// Unmapped reads in SAM output: None or Within #[arg(long = "outSAMunmapped", default_value = "None")] pub out_sam_unmapped: OutSamUnmapped, @@ -602,8 +608,11 @@ impl Parameters { /// - `"All"` → {NH, HI, AS, NM, nM, MD, jM, jI, XS} /// - `"None"` → {} (empty) /// - Explicit list (e.g. ["NH", "AS"]) → collected as-is + /// + /// `RG` is auto-appended when `--outSAMattrRGline` is set (STAR behavior, + /// `Parameters_samAttributes.cpp:201`). pub fn sam_attribute_set(&self) -> HashSet { - match self + let mut attrs: HashSet = match self .out_sam_attributes .iter() .map(String::as_str) @@ -620,6 +629,82 @@ impl Parameters { .collect(), ["None"] => HashSet::new(), tags => tags.iter().map(|s| s.to_string()).collect(), + }; + if self.rg_line_set() { + attrs.insert("RG".to_string()); + } + attrs + } + + /// True if the user provided a non-default `--outSAMattrRGline`. + pub fn rg_line_set(&self) -> bool { + !self.out_sam_attr_rg_line.is_empty() && self.out_sam_attr_rg_line[0] != "-" + } + + /// Parse `--outSAMattrRGline` into one tab-joined body per `@RG` block. + /// + /// Mirrors `Parameters_readFilesInit.cpp:65-82`: tokens are joined by tabs + /// until a bare `,` token starts a new block. Each block's first token must + /// begin with `ID:`. + pub fn parsed_rg_lines(&self) -> Result, crate::error::Error> { + if !self.rg_line_set() { + return Ok(Vec::new()); + } + let mut lines: Vec = Vec::new(); + let mut i = 0; + while i < self.out_sam_attr_rg_line.len() { + let tok = &self.out_sam_attr_rg_line[i]; + if i > 0 && tok == "," { + i += 1; + if i >= self.out_sam_attr_rg_line.len() { + break; + } + } + let first = &self.out_sam_attr_rg_line[i]; + if !first.starts_with("ID:") { + return Err(crate::error::Error::Parameter(format!( + "--outSAMattrRGline: first field of each RG line must start with 'ID:', got '{}'", + first + ))); + } + let mut body = first.clone(); + i += 1; + while i < self.out_sam_attr_rg_line.len() && self.out_sam_attr_rg_line[i] != "," { + body.push('\t'); + body.push_str(&self.out_sam_attr_rg_line[i]); + i += 1; + } + lines.push(body); + } + Ok(lines) + } + + /// Per-file read group ID, replicated from a single RG line if needed. + /// Returns empty vec when no RG line is set. + pub fn rg_ids(&self) -> Result, crate::error::Error> { + let lines = self.parsed_rg_lines()?; + if lines.is_empty() { + return Ok(Vec::new()); + } + let ids: Vec = lines + .iter() + .map(|body| { + let first = body.split('\t').next().unwrap_or(""); + first.trim_start_matches("ID:").to_string() + }) + .collect(); + let n_files = self.read_files_in.len().max(1); + if ids.len() > 1 && ids.len() != n_files { + return Err(crate::error::Error::Parameter(format!( + "--outSAMattrRGline: {} RG entries does not match --readFilesIn count {} (must be 1 or N)", + ids.len(), + n_files + ))); + } + if ids.len() == 1 { + Ok(vec![ids[0].clone(); n_files]) + } else { + Ok(ids) } } @@ -714,6 +799,18 @@ impl Parameters { )); } + // Read group: `RG` in outSAMattributes without an RG line is a fatal + // error (STAR: Parameters_samAttributes.cpp:206). STAR's "All" preset + // does NOT include RG, so only match a literal "RG" token here. Also + // parse the RG line to validate its ID: prefix and per-file RG count. + let user_wants_rg_attr = self.out_sam_attributes.iter().any(|a| a == "RG"); + if !self.rg_line_set() && user_wants_rg_attr { + return Err(crate::error::Error::Parameter( + "--outSAMattributes contains RG tag, but --outSAMattrRGline is not set".into(), + )); + } + self.rg_ids()?; + Ok(()) } @@ -997,6 +1094,99 @@ mod tests { assert_eq!(p.win_bin_window_dist(), 81_920); // 2^14 * 5 } + #[test] + fn rg_line_default_unset() { + let p = parse(&["--readFilesIn", "r.fq"]); + assert!(!p.rg_line_set()); + assert_eq!(p.parsed_rg_lines().unwrap(), Vec::::new()); + assert_eq!(p.rg_ids().unwrap(), Vec::::new()); + assert!(!p.sam_attribute_set().contains("RG")); + } + + #[test] + fn rg_line_single() { + let p = parse(&[ + "--readFilesIn", + "r.fq", + "--outSAMattrRGline", + "ID:foo", + "SM:bar", + "LB:lib1", + ]); + assert!(p.rg_line_set()); + assert_eq!( + p.parsed_rg_lines().unwrap(), + vec!["ID:foo\tSM:bar\tLB:lib1".to_string()] + ); + assert_eq!(p.rg_ids().unwrap(), vec!["foo".to_string()]); + assert!(p.sam_attribute_set().contains("RG")); + } + + #[test] + fn rg_line_multi() { + let p = parse(&[ + "--readFilesIn", + "r1.fq", + "r2.fq", + "--outSAMattrRGline", + "ID:a", + "SM:a", + ",", + "ID:b", + "LB:x", + ]); + let lines = p.parsed_rg_lines().unwrap(); + assert_eq!( + lines, + vec!["ID:a\tSM:a".to_string(), "ID:b\tLB:x".to_string()] + ); + assert_eq!(p.rg_ids().unwrap(), vec!["a".to_string(), "b".to_string()]); + } + + #[test] + fn rg_line_single_replicates_for_multi_file() { + let p = parse(&[ + "--readFilesIn", + "r1.fq", + "r2.fq", + "--outSAMattrRGline", + "ID:foo", + ]); + assert_eq!( + p.rg_ids().unwrap(), + vec!["foo".to_string(), "foo".to_string()] + ); + } + + #[test] + fn rg_line_missing_id_prefix_errors() { + let p = parse(&["--readFilesIn", "r.fq", "--outSAMattrRGline", "SM:oops"]); + let err = p.parsed_rg_lines().unwrap_err(); + assert!(err.to_string().contains("ID:")); + } + + #[test] + fn rg_line_count_mismatch_errors() { + // 1 input file, 2 RG entries — mismatch (ids.len()>1 && != n_files). + let p = parse(&[ + "--readFilesIn", + "r1.fq", + "--outSAMattrRGline", + "ID:a", + ",", + "ID:b", + ]); + let err = p.rg_ids().unwrap_err(); + assert!(err.to_string().contains("does not match")); + } + + #[test] + fn validate_rg_attr_without_line_errors() { + let p = parse(&["--readFilesIn", "r.fq", "--outSAMattributes", "NH", "RG"]); + let err = p.validate().unwrap_err(); + assert!(err.to_string().contains("RG")); + } + #[test] fn sj_stitch_mismatch() { let p = parse(&[ From 716a8162fc78e9c62cca2b825aa96fcfa6264a87 Mon Sep 17 00:00:00 2001 From: Phil Ewels Date: Fri, 17 Apr 2026 15:22:12 +0200 Subject: [PATCH 2/7] chore(simplify): drop dead rg_ids() call and redundant attrs gate - Remove the `let _ = params.rg_ids()?;` early-validation line in `align_reads_paired_end`. `validate()` already calls `rg_ids()` at startup, so the second call was a no-op. - Simplify `maybe_insert_rg_tag` to `(record, rg_id)` by dropping the `attrs: &HashSet` parameter and its `attrs.contains("RG")` gate. `sam_attribute_set()` auto-inserts "RG" whenever an RG line is configured, so `rg_id.is_some()` already implies the tag is wanted. - Replace the inline `if let Some(id) = rg_id { ... }` block in `build_unmapped_record` with a call to the simplified helper, so every writer-side RG insertion goes through the same path. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/io/sam.rs | 31 ++++++++++++++----------------- src/lib.rs | 1 - 2 files changed, 14 insertions(+), 18 deletions(-) diff --git a/src/io/sam.rs b/src/io/sam.rs index 551cd945..f6c3b51c 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -122,7 +122,7 @@ impl SamWriter { hit_index + 1, // 1-based &attrs, )?; - maybe_insert_rg_tag(&mut record, &attrs, rg_id); + maybe_insert_rg_tag(&mut record, rg_id); self.writer.write_alignment_record(&self.header, &record)?; } @@ -214,11 +214,7 @@ impl SamWriter { // Quality scores *record.quality_scores_mut() = QualityScores::from(read_qual.to_vec()); - if let Some(id) = rg_id { - record - .data_mut() - .insert(Tag::READ_GROUP, Value::String(BString::from(id))); - } + maybe_insert_rg_tag(&mut record, rg_id); Ok(record) } @@ -273,7 +269,7 @@ impl SamWriter { hit_index + 1, // 1-based &attrs, )?; - maybe_insert_rg_tag(&mut record, &attrs, rg_id); + maybe_insert_rg_tag(&mut record, rg_id); records.push(record); } @@ -352,7 +348,7 @@ impl SamWriter { combined_score, &attrs, )?; - maybe_insert_rg_tag(&mut rec1, &attrs, rg_id); + maybe_insert_rg_tag(&mut rec1, rg_id); records.push(rec1); // Create record for mate2 (this=mate2, mate=mate1) @@ -372,7 +368,7 @@ impl SamWriter { combined_score, &attrs, )?; - maybe_insert_rg_tag(&mut rec2, &attrs, rg_id); + maybe_insert_rg_tag(&mut rec2, rg_id); records.push(rec2); } @@ -514,7 +510,7 @@ impl SamWriter { ); data.insert(Tag::new(b'M', b'D'), Value::String(BString::from(md))); } - maybe_insert_rg_tag(&mut mapped_rec, &attrs, rg_id); + maybe_insert_rg_tag(&mut mapped_rec, rg_id); // --- Build unmapped mate record --- let mut unmapped_rec = RecordBuf::default(); @@ -554,7 +550,7 @@ impl SamWriter { let unmapped_seq_bytes: Vec = unmapped_seq.iter().map(|&b| decode_base(b)).collect(); *unmapped_rec.sequence_mut() = Sequence::from(unmapped_seq_bytes); *unmapped_rec.quality_scores_mut() = QualityScores::from(unmapped_qual.to_vec()); - maybe_insert_rg_tag(&mut unmapped_rec, &attrs, rg_id); + maybe_insert_rg_tag(&mut unmapped_rec, rg_id); // Order: mate1 first, mate2 second if mate1_is_mapped { @@ -578,7 +574,6 @@ impl SamWriter { params: &Parameters, ) -> Result, Error> { let mut records = Vec::with_capacity(2); - let attrs = params.sam_attribute_set(); let rg_ids = params.rg_ids()?; let rg_id = rg_ids.first().map(String::as_str); @@ -596,7 +591,7 @@ impl SamWriter { let seq1_bytes: Vec = mate1_seq.iter().map(|&b| decode_base(b)).collect(); *rec1.sequence_mut() = Sequence::from(seq1_bytes); *rec1.quality_scores_mut() = QualityScores::from(mate1_qual.to_vec()); - maybe_insert_rg_tag(&mut rec1, &attrs, rg_id); + maybe_insert_rg_tag(&mut rec1, rg_id); records.push(rec1); // Mate2 record @@ -613,7 +608,7 @@ impl SamWriter { let seq2_bytes: Vec = mate2_seq.iter().map(|&b| decode_base(b)).collect(); *rec2.sequence_mut() = Sequence::from(seq2_bytes); *rec2.quality_scores_mut() = QualityScores::from(mate2_qual.to_vec()); - maybe_insert_rg_tag(&mut rec2, &attrs, rg_id); + maybe_insert_rg_tag(&mut rec2, rg_id); records.push(rec2); Ok(records) @@ -678,9 +673,11 @@ pub fn build_sam_header(genome: &Genome, params: &Parameters) -> Result` on the record when RG is in `attrs` and an ID is set. -fn maybe_insert_rg_tag(record: &mut RecordBuf, attrs: &HashSet, rg_id: Option<&str>) { - if let (true, Some(id)) = (attrs.contains("RG"), rg_id) { +/// Insert `RG:Z:` on the record when an ID is set. `sam_attribute_set()` +/// auto-adds `RG` to the attribute set whenever an RG line is configured, so +/// `rg_id.is_some()` implies the attribute is wanted — no extra gate needed. +fn maybe_insert_rg_tag(record: &mut RecordBuf, rg_id: Option<&str>) { + if let Some(id) = rg_id { record .data_mut() .insert(Tag::READ_GROUP, Value::String(BString::from(id))); diff --git a/src/lib.rs b/src/lib.rs index e67f8c92..f2a72cf7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -781,7 +781,6 @@ fn align_reads_paired_end( let max_multimaps = params.out_filter_multimap_nmax as usize; let output_unmapped = params.out_sam_unmapped != params::OutSamUnmapped::None; let by_sjout = params.out_filter_type == OutFilterType::BySJout; - let _ = params.rg_ids()?; // Validate early; per-call uses params inside builders. // Buffer for BySJout mode let mut bysj_buffer: Vec = Vec::new(); From f23e00db381170f2013aa43ff0a951e0cb46fb5a Mon Sep 17 00:00:00 2001 From: Phil Ewels Date: Fri, 17 Apr 2026 15:23:42 +0200 Subject: [PATCH 3/7] chore(simplify): introduce primary_rg_id helper for writer callsites Add `Parameters::primary_rg_id() -> Result, Error>` that returns the first RG line's `ID:` value directly, and use it to replace the `let rg_ids = params.rg_ids()?; let rg_id = rg_ids.first()...` two-step pattern at every SAM/BAM writer callsite (5 in sam.rs, 1 in lib.rs). `rg_ids()` is retained because `validate()` still uses it for the 1-vs-N count check against `--readFilesIn`. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/io/sam.rs | 20 ++++++++++---------- src/lib.rs | 3 +-- src/params.rs | 11 +++++++++++ 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/io/sam.rs b/src/io/sam.rs index f6c3b51c..0c812beb 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -107,8 +107,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } - let rg_ids = params.rg_ids()?; - let rg_id = rg_ids.first().map(String::as_str); + let rg_id_owned = params.primary_rg_id()?; + let rg_id = rg_id_owned.as_deref(); for (hit_index, transcript) in transcripts.iter().take(max_output).enumerate() { let mut record = transcript_to_record( @@ -253,8 +253,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } - let rg_ids = params.rg_ids()?; - let rg_id = rg_ids.first().map(String::as_str); + let rg_id_owned = params.primary_rg_id()?; + let rg_id = rg_id_owned.as_deref(); let mut records = Vec::with_capacity(max_output); for (hit_index, transcript) in transcripts.iter().take(max_output).enumerate() { @@ -320,8 +320,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } - let rg_ids = params.rg_ids()?; - let rg_id = rg_ids.first().map(String::as_str); + let rg_id_owned = params.primary_rg_id()?; + let rg_id = rg_id_owned.as_deref(); let mut records = Vec::with_capacity(max_output * 2); @@ -406,8 +406,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } - let rg_ids = params.rg_ids()?; - let rg_id = rg_ids.first().map(String::as_str); + let rg_id_owned = params.primary_rg_id()?; + let rg_id = rg_id_owned.as_deref(); // Compute mapped mate's per-chr position for co-location let chr_start = genome.chr_start[mapped_transcript.chr_idx]; @@ -574,8 +574,8 @@ impl SamWriter { params: &Parameters, ) -> Result, Error> { let mut records = Vec::with_capacity(2); - let rg_ids = params.rg_ids()?; - let rg_id = rg_ids.first().map(String::as_str); + let rg_id_owned = params.primary_rg_id()?; + let rg_id = rg_id_owned.as_deref(); // Mate1 record let mut rec1 = RecordBuf::default(); diff --git a/src/lib.rs b/src/lib.rs index f2a72cf7..31f13909 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -472,8 +472,7 @@ fn align_reads_single_end( let max_multimaps = params.out_filter_multimap_nmax as usize; let output_unmapped = params.out_sam_unmapped != params::OutSamUnmapped::None; let by_sjout = params.out_filter_type == OutFilterType::BySJout; - let rg_ids = params.rg_ids()?; - let rg_id_owned = rg_ids.first().cloned(); + let rg_id_owned = params.primary_rg_id()?; // Buffer for BySJout mode: accumulate all results before filtering let mut bysj_buffer: Vec = Vec::new(); diff --git a/src/params.rs b/src/params.rs index 2a9db479..e6e83596 100644 --- a/src/params.rs +++ b/src/params.rs @@ -679,6 +679,17 @@ impl Parameters { Ok(lines) } + /// Read group ID emitted on SAM records (the first RG line's `ID:` value). + /// Returns `None` when no RG line is configured. + pub fn primary_rg_id(&self) -> Result, crate::error::Error> { + Ok(self.parsed_rg_lines()?.first().and_then(|body| { + body.split('\t') + .next()? + .strip_prefix("ID:") + .map(str::to_owned) + })) + } + /// Per-file read group ID, replicated from a single RG line if needed. /// Returns empty vec when no RG line is set. pub fn rg_ids(&self) -> Result, crate::error::Error> { From 845398f86f16e154c98c91a26a90e130c21324aa Mon Sep 17 00:00:00 2001 From: Phil Ewels Date: Fri, 17 Apr 2026 15:24:25 +0200 Subject: [PATCH 4/7] chore(simplify): rewrite parsed_rg_lines with Vec::split Replace the manual `while i < len` index-walk with a `out_sam_attr_rg_line.split(|tok| tok == ",")` that maps each block to its tab-joined body and propagates any `ID:`-prefix error via `collect::>()`. Same public signature and return value as before, but with half the local state. Behavior note: an empty block (adjacent `,`s or a trailing `,`) now returns a `Parameter` error instead of being silently skipped. This matches STAR's `Parameters_readFilesInit.cpp` which requires every RG block to begin with `ID:`. Existing tests still pass. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/params.rs | 48 +++++++++++++++++++----------------------------- 1 file changed, 19 insertions(+), 29 deletions(-) diff --git a/src/params.rs b/src/params.rs index e6e83596..307c76ce 100644 --- a/src/params.rs +++ b/src/params.rs @@ -643,40 +643,30 @@ impl Parameters { /// Parse `--outSAMattrRGline` into one tab-joined body per `@RG` block. /// - /// Mirrors `Parameters_readFilesInit.cpp:65-82`: tokens are joined by tabs - /// until a bare `,` token starts a new block. Each block's first token must - /// begin with `ID:`. + /// Mirrors `Parameters_readFilesInit.cpp:65-82`: tokens are split on bare + /// `,` separators, and each resulting block's first token must begin with + /// `ID:`. An empty block (adjacent commas or a trailing comma) is an error. pub fn parsed_rg_lines(&self) -> Result, crate::error::Error> { if !self.rg_line_set() { return Ok(Vec::new()); } - let mut lines: Vec = Vec::new(); - let mut i = 0; - while i < self.out_sam_attr_rg_line.len() { - let tok = &self.out_sam_attr_rg_line[i]; - if i > 0 && tok == "," { - i += 1; - if i >= self.out_sam_attr_rg_line.len() { - break; + self.out_sam_attr_rg_line + .split(|tok| tok == ",") + .map(|block| { + let first = block.first().ok_or_else(|| { + crate::error::Error::Parameter( + "--outSAMattrRGline: empty RG block".into(), + ) + })?; + if !first.starts_with("ID:") { + return Err(crate::error::Error::Parameter(format!( + "--outSAMattrRGline: first field of each RG line must start with 'ID:', got '{}'", + first + ))); } - } - let first = &self.out_sam_attr_rg_line[i]; - if !first.starts_with("ID:") { - return Err(crate::error::Error::Parameter(format!( - "--outSAMattrRGline: first field of each RG line must start with 'ID:', got '{}'", - first - ))); - } - let mut body = first.clone(); - i += 1; - while i < self.out_sam_attr_rg_line.len() && self.out_sam_attr_rg_line[i] != "," { - body.push('\t'); - body.push_str(&self.out_sam_attr_rg_line[i]); - i += 1; - } - lines.push(body); - } - Ok(lines) + Ok(block.join("\t")) + }) + .collect() } /// Read group ID emitted on SAM records (the first RG line's `ID:` value). From 471518ae59f250badde265b8a9f445fa90945fc9 Mon Sep 17 00:00:00 2001 From: Phil Ewels Date: Fri, 17 Apr 2026 15:24:59 +0200 Subject: [PATCH 5/7] chore: remove dead SamWriter::write_unmapped MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `SamWriter::write_unmapped` had no production callers — the SE/PE write paths build unmapped records via `build_unmapped_record` and batch them through `write_batch`. The only caller was the method's own self-test. Delete both. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/io/sam.rs | 32 -------------------------------- 1 file changed, 32 deletions(-) diff --git a/src/io/sam.rs b/src/io/sam.rs index 0c812beb..f7318be6 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -130,23 +130,6 @@ impl SamWriter { Ok(()) } - /// Write unmapped read - /// - /// # Arguments - /// * `read_name` - Read identifier - /// * `read_seq` - Read sequence (encoded) - /// * `read_qual` - Quality scores - pub fn write_unmapped( - &mut self, - read_name: &str, - read_seq: &[u8], - read_qual: &[u8], - ) -> Result<(), Error> { - let record = Self::build_unmapped_record(read_name, read_seq, read_qual, None)?; - self.writer.write_alignment_record(&self.header, &record)?; - Ok(()) - } - /// Write batch of buffered records (for parallel processing) /// /// # Arguments @@ -1300,21 +1283,6 @@ mod tests { assert!(writer.is_ok()); } - #[test] - fn test_write_unmapped() { - let genome = make_test_genome(); - let params = Parameters::parse_from(vec!["ruSTAR", "--readFilesIn", "test.fq"]); - - let tmpfile = NamedTempFile::new().unwrap(); - let mut writer = SamWriter::create(tmpfile.path(), &genome, ¶ms).unwrap(); - - let read_seq = vec![0, 1, 2, 3]; // ACGT - let read_qual = vec![30, 30, 30, 30]; - - let result = writer.write_unmapped("read1", &read_seq, &read_qual); - assert!(result.is_ok()); - } - #[test] fn test_transcript_to_record() { let genome = make_test_genome(); From 8dc42620b16c15b10e1f8d74beaf659353b9c026 Mon Sep 17 00:00:00 2001 From: Phil Ewels Date: Fri, 17 Apr 2026 15:47:40 +0200 Subject: [PATCH 6/7] chore: apply cargo fmt Fixes formatting drift inherited from main. No semantic changes. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/align/read_align.rs | 45 +++++++++++++----------------- src/align/stitch.rs | 3 -- src/lib.rs | 62 +++++++++++++++++++++++++++++++++++------ src/quant/mod.rs | 33 ++++++++++++---------- 4 files changed, 91 insertions(+), 52 deletions(-) diff --git a/src/align/read_align.rs b/src/align/read_align.rs index c53daef6..524f5f0c 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -1,9 +1,7 @@ /// Read alignment driver function use crate::align::score::{AlignmentScorer, SpliceMotif}; use crate::align::seed::Seed; -use crate::align::stitch::{ - cluster_seeds, stitch_seeds_with_jdb_debug, -}; +use crate::align::stitch::{cluster_seeds, stitch_seeds_with_jdb_debug}; use crate::align::transcript::Transcript; use crate::error::Error; use crate::index::GenomeIndex; @@ -621,7 +619,10 @@ pub fn align_paired_read( // This correctly represents mate2 on - strand for FR pairs without explicit RC handling. let debug_name: &str = if debug_pe { read_name } else { "" }; let mut mate1_transcripts: Vec = Vec::new(); - for cluster in mate1_clusters.iter().take(params.align_windows_per_read_nmax) { + for cluster in mate1_clusters + .iter() + .take(params.align_windows_per_read_nmax) + { let ts = stitch_seeds_with_jdb_debug( cluster, mate1_seq, @@ -635,7 +636,10 @@ pub fn align_paired_read( } let mut mate2_transcripts: Vec = Vec::new(); - for cluster in mate2_clusters.iter().take(params.align_windows_per_read_nmax) { + for cluster in mate2_clusters + .iter() + .take(params.align_windows_per_read_nmax) + { let ts = stitch_seeds_with_jdb_debug( cluster, mate2_seq, @@ -656,14 +660,9 @@ pub fn align_paired_read( let mut joint_pairs: Vec = Vec::new(); for t1 in &mate1_transcripts { for t2 in &mate2_transcripts { - if let Some(pair) = try_pair_transcripts( - t1, - t2, - len1, - len2, - params, - combined_score_threshold, - ) { + if let Some(pair) = + try_pair_transcripts(t1, t2, len1, len2, params, combined_score_threshold) + { joint_pairs.push(pair); } } @@ -840,12 +839,10 @@ pub fn align_paired_read( // Half-mapped fallback: report the best-scoring single-mate transcript. // Apply per-mate quality threshold (outFilterScoreMinOverLread * (len - 1)). - let thresh1 = - ((params.out_filter_score_min_over_lread * (len1 as f64 - 1.0)) as i32) - .max(params.out_filter_score_min); - let thresh2 = - ((params.out_filter_score_min_over_lread * (len2 as f64 - 1.0)) as i32) - .max(params.out_filter_score_min); + let thresh1 = ((params.out_filter_score_min_over_lread * (len1 as f64 - 1.0)) as i32) + .max(params.out_filter_score_min); + let thresh2 = ((params.out_filter_score_min_over_lread * (len2 as f64 - 1.0)) as i32) + .max(params.out_filter_score_min); let best_m1 = mate1_transcripts .into_iter() @@ -1041,9 +1038,7 @@ fn filter_paired_transcripts(paired_alns: &mut Vec, params: &Pa // the read is unmapped. // // Step 1: find the best pair and check quality thresholds on it. - let best_pa = paired_alns - .iter() - .max_by_key(|pa| pa.combined_wt_score); + let best_pa = paired_alns.iter().max_by_key(|pa| pa.combined_wt_score); if let Some(best) = best_pa { let mate1_len = (best.mate1_region.1 - best.mate1_region.0) as f64; let mate2_len = (best.mate2_region.1 - best.mate2_region.0) as f64; @@ -1059,8 +1054,7 @@ fn filter_paired_transcripts(paired_alns: &mut Vec, params: &Pa // per-mate spans are small (penalty ≈ −2 each → sum penalty −4 vs combined penalty −2). let combined_score = best.combined_wt_score; if combined_score < params.out_filter_score_min - || combined_score - < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32 + || combined_score < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32 { paired_alns.clear(); return; @@ -1074,8 +1068,7 @@ fn filter_paired_transcripts(paired_alns: &mut Vec, params: &Pa // which directly mirrors STAR's joint transcript nMatch without extension inflation. let combined_match = best.combined_n_match; if combined_match < params.out_filter_match_nmin - || combined_match - < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32 + || combined_match < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32 { paired_alns.clear(); return; diff --git a/src/align/stitch.rs b/src/align/stitch.rs index e2ceeffc..77ac2ece 100644 --- a/src/align/stitch.rs +++ b/src/align/stitch.rs @@ -2407,7 +2407,6 @@ pub(crate) fn stitch_seeds_with_jdb_debug( Ok(transcripts) } - /// Shared core: preprocessing + recursive stitcher, returns working transcripts + context. #[allow(clippy::too_many_arguments)] fn stitch_seeds_core( @@ -2750,8 +2749,6 @@ fn stitch_seeds_core( )) } - - #[cfg(test)] mod tests { use super::*; diff --git a/src/lib.rs b/src/lib.rs index 31f13909..78a1628c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -144,7 +144,10 @@ fn align_reads(params: &Parameters) -> anyhow::Result<()> { let quant_ctx: Option> = if params.quant_gene_counts() { let gtf_path = params.sjdb_gtf_file.as_ref().unwrap(); - info!("quantMode GeneCounts: building gene annotation from {}", gtf_path.display()); + info!( + "quantMode GeneCounts: building gene annotation from {}", + gtf_path.display() + ); let ctx = crate::quant::QuantContext::build(gtf_path, &index.genome)?; Some(std::sync::Arc::new(ctx)) } else { @@ -223,8 +226,22 @@ fn run_single_pass( // Route to single-end or paired-end mode match params.read_files_in.len() { - 1 => align_reads_single_end(params, index, &mut writer, &stats, &sj_stats, quant.as_ref()), - 2 => align_reads_paired_end(params, index, &mut writer, &stats, &sj_stats, quant.as_ref()), + 1 => align_reads_single_end( + params, + index, + &mut writer, + &stats, + &sj_stats, + quant.as_ref(), + ), + 2 => align_reads_paired_end( + params, + index, + &mut writer, + &stats, + &sj_stats, + quant.as_ref(), + ), n => anyhow::bail!("Invalid number of read files: {} (expected 1 or 2)", n), }?; } @@ -241,8 +258,22 @@ fn run_single_pass( // Route to single-end or paired-end mode (same functions as SAM, generic!) match params.read_files_in.len() { - 1 => align_reads_single_end(params, index, &mut writer, &stats, &sj_stats, quant.as_ref()), - 2 => align_reads_paired_end(params, index, &mut writer, &stats, &sj_stats, quant.as_ref()), + 1 => align_reads_single_end( + params, + index, + &mut writer, + &stats, + &sj_stats, + quant.as_ref(), + ), + 2 => align_reads_paired_end( + params, + index, + &mut writer, + &stats, + &sj_stats, + quant.as_ref(), + ), n => anyhow::bail!("Invalid number of read files: {} (expected 1 or 2)", n), }?; @@ -345,8 +376,22 @@ fn run_pass1( // Align reads (single-end or paired-end); no quant counting in pass 1 match params.read_files_in.len() { - 1 => align_reads_single_end(¶ms_pass1, index, &mut null_writer, &stats, &sj_stats, None)?, - 2 => align_reads_paired_end(¶ms_pass1, index, &mut null_writer, &stats, &sj_stats, None)?, + 1 => align_reads_single_end( + ¶ms_pass1, + index, + &mut null_writer, + &stats, + &sj_stats, + None, + )?, + 2 => align_reads_paired_end( + ¶ms_pass1, + index, + &mut null_writer, + &stats, + &sj_stats, + None, + )?, n => anyhow::bail!("Invalid number of read files: {} (expected 1 or 2)", n), } @@ -574,7 +619,8 @@ fn align_reads_single_end( // Gene-level quantification (lock-free atomic counts) if let Some(ref q) = quant { - q.counts.count_se_read(&transcripts, n_for_mapq, &q.gene_ann); + q.counts + .count_se_read(&transcripts, n_for_mapq, &q.gene_ann); } // Record junction statistics diff --git a/src/quant/mod.rs b/src/quant/mod.rs index 41c2d7a5..f9884f96 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -234,13 +234,8 @@ impl GeneCounts { } /// Write `ReadsPerGene.out.tab` in STAR's format. - pub fn write_output( - &self, - path: &Path, - gene_ann: &GeneAnnotation, - ) -> Result<(), Error> { - let mut file = std::fs::File::create(path) - .map_err(|e| Error::io(e, path))?; + pub fn write_output(&self, path: &Path, gene_ann: &GeneAnnotation) -> Result<(), Error> { + let mut file = std::fs::File::create(path).map_err(|e| Error::io(e, path))?; macro_rules! wl { ($($arg:tt)*) => { @@ -350,10 +345,7 @@ impl QuantContext { let exons = crate::junction::gtf::parse_gtf(gtf_path)?; let gene_ann = GeneAnnotation::from_gtf_exons(&exons, genome); let n = gene_ann.n_genes(); - log::info!( - "quantMode GeneCounts: {} genes loaded from GTF", - n - ); + log::info!("quantMode GeneCounts: {} genes loaded from GTF", n); let counts = GeneCounts::new(n); Ok(QuantContext { gene_ann, counts }) } @@ -381,7 +373,13 @@ mod tests { } } - fn make_gtf_exon(seqname: &str, start: u64, end: u64, strand: char, gene_id: &str) -> GtfRecord { + fn make_gtf_exon( + seqname: &str, + start: u64, + end: u64, + strand: char, + gene_id: &str, + ) -> GtfRecord { let mut attrs = std::collections::HashMap::new(); attrs.insert("gene_id".to_string(), gene_id.to_string()); attrs.insert("transcript_id".to_string(), "T1".to_string()); @@ -401,7 +399,12 @@ mod tests { genome_start: gs, genome_end: ge, is_reverse, - exons: vec![Exon { genome_start: gs, genome_end: ge, read_start: 0, read_end: (ge - gs) as usize }], + exons: vec![Exon { + genome_start: gs, + genome_end: ge, + read_start: 0, + read_end: (ge - gs) as usize, + }], cigar: vec![], score: 100, n_mismatch: 0, @@ -417,7 +420,7 @@ mod tests { fn test_gene_annotation_basic() { let genome = make_genome(); let exons = vec![ - make_gtf_exon("chr1", 101, 200, '+', "G1"), // 0-based: [100, 200) + make_gtf_exon("chr1", 101, 200, '+', "G1"), // 0-based: [100, 200) make_gtf_exon("chr1", 301, 400, '+', "G1"), make_gtf_exon("chr1", 501, 600, '-', "G2"), ]; @@ -426,7 +429,7 @@ mod tests { assert_eq!(ann.gene_ids[0], "G1"); assert_eq!(ann.gene_ids[1], "G2"); assert!(!ann.gene_is_reverse[0]); // G1 is + - assert!(ann.gene_is_reverse[1]); // G2 is - + assert!(ann.gene_is_reverse[1]); // G2 is - assert_eq!(ann.chr_exons[0].len(), 3); } From e3060461623c68a6bbcbdbe6de83c9fe38122978 Mon Sep 17 00:00:00 2001 From: James Ferguson Date: Mon, 20 Apr 2026 13:30:19 +1000 Subject: [PATCH 7/7] Fix merge conflict squashing } braket --- src/params.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/params.rs b/src/params.rs index b89d4e50..88a99331 100644 --- a/src/params.rs +++ b/src/params.rs @@ -1191,6 +1191,8 @@ mod tests { let p = parse(&["--readFilesIn", "r.fq", "--outSAMattributes", "NH", "RG"]); let err = p.validate().unwrap_err(); assert!(err.to_string().contains("RG")); + } + fn run_rng_seed_override() { let p = parse(&["--readFilesIn", "r.fq", "--runRNGseed", "42"]); assert_eq!(p.run_rng_seed, 42);