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..f7318be6 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_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 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, rg_id); self.writer.write_alignment_record(&self.header, &record)?; } @@ -124,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)?; - self.writer.write_alignment_record(&self.header, &record)?; - Ok(()) - } - /// Write batch of buffered records (for parallel processing) /// /// # Arguments @@ -185,10 +174,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 +197,8 @@ impl SamWriter { // Quality scores *record.quality_scores_mut() = QualityScores::from(read_qual.to_vec()); + maybe_insert_rg_tag(&mut record, rg_id); + Ok(record) } @@ -243,10 +236,12 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + 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() { - let record = transcript_to_record( + let mut record = transcript_to_record( transcript, read_name, read_seq, @@ -257,6 +252,7 @@ impl SamWriter { hit_index + 1, // 1-based &attrs, )?; + maybe_insert_rg_tag(&mut record, rg_id); records.push(record); } @@ -291,7 +287,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 +303,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + 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); @@ -317,7 +315,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 +331,11 @@ impl SamWriter { combined_score, &attrs, )?; + maybe_insert_rg_tag(&mut rec1, 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 +351,7 @@ impl SamWriter { combined_score, &attrs, )?; + maybe_insert_rg_tag(&mut rec2, rg_id); records.push(rec2); } @@ -389,6 +389,8 @@ impl SamWriter { if params.out_sam_strand_field != "intronMotif" { attrs.remove("XS"); } + 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]; @@ -491,6 +493,7 @@ impl SamWriter { ); data.insert(Tag::new(b'M', b'D'), Value::String(BString::from(md))); } + maybe_insert_rg_tag(&mut mapped_rec, rg_id); // --- Build unmapped mate record --- let mut unmapped_rec = RecordBuf::default(); @@ -530,6 +533,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, rg_id); // Order: mate1 first, mate2 second if mate1_is_mapped { @@ -550,8 +554,11 @@ impl SamWriter { mate1_qual: &[u8], mate2_seq: &[u8], mate2_qual: &[u8], + params: &Parameters, ) -> Result, Error> { let mut records = Vec::with_capacity(2); + let rg_id_owned = params.primary_rg_id()?; + let rg_id = rg_id_owned.as_deref(); // Mate1 record let mut rec1 = RecordBuf::default(); @@ -567,6 +574,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, rg_id); records.push(rec1); // Mate2 record @@ -583,6 +591,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, rg_id); records.push(rec2); Ok(records) @@ -590,7 +599,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 +619,54 @@ 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 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))); + } +} + /// Convert Transcript to SAM record #[allow(clippy::too_many_arguments)] fn transcript_to_record( @@ -1130,28 +1181,106 @@ mod tests { } #[test] - fn test_sam_writer_creation() { + fn test_build_sam_header_with_rg() { let genome = make_test_genome(); - let params = Parameters::parse_from(vec!["ruSTAR", "--readFilesIn", "test.fq"]); - - let tmpfile = NamedTempFile::new().unwrap(); - let writer = SamWriter::create(tmpfile.path(), &genome, ¶ms); - assert!(writer.is_ok()); + 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_write_unmapped() { + 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"]); + 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 read_seq = vec![0, 1, 2, 3]; // ACGT - let read_qual = vec![30, 30, 30, 30]; + 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(); + let params = Parameters::parse_from(vec!["ruSTAR", "--readFilesIn", "test.fq"]); - let result = writer.write_unmapped("read1", &read_seq, &read_qual); - assert!(result.is_ok()); + let tmpfile = NamedTempFile::new().unwrap(); + let writer = SamWriter::create(tmpfile.path(), &genome, ¶ms); + assert!(writer.is_ok()); } #[test] @@ -1207,6 +1336,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 +1344,7 @@ mod tests { &mate1_qual, &mate2_seq, &mate2_qual, + ¶ms, ) .unwrap(); diff --git a/src/lib.rs b/src/lib.rs index 227b706a..78a1628c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -517,6 +517,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_id_owned = params.primary_rg_id()?; // Buffer for BySJout mode: accumulate all results before filtering let mut bysj_buffer: Vec = Vec::new(); @@ -576,6 +577,7 @@ fn align_reads_single_end( &read.name, &clipped_seq, &clipped_qual, + rg_id_owned.as_deref(), )?; buffer.push(record); } @@ -643,6 +645,7 @@ fn align_reads_single_end( &read.name, &clipped_seq, &clipped_qual, + rg_id_owned.as_deref(), )?; buffer.push(record); } @@ -893,6 +896,7 @@ fn align_reads_paired_end( &m1_qual, &m2_seq, &m2_qual, + params, )?; for record in records { buffer.push(record); @@ -1030,6 +1034,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 10e9d5da..88a99331 100644 --- a/src/params.rs +++ b/src/params.rs @@ -306,6 +306,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, @@ -606,8 +612,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) @@ -624,6 +633,83 @@ 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 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()); + } + 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 + ))); + } + Ok(block.join("\t")) + }) + .collect() + } + + /// 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> { + 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) } } @@ -718,6 +804,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(()) } @@ -1003,6 +1101,98 @@ mod tests { } #[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")); + } + fn run_rng_seed_override() { let p = parse(&["--readFilesIn", "r.fq", "--runRNGseed", "42"]); assert_eq!(p.run_rng_seed, 42);