From e8e26c34eae7b340f6ffb077707360088b2e86f2 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 8 Oct 2026 20:56:15 +0200 Subject: [PATCH 1/3] fix(sam): write STAR's @HD VN:1.4 header version STAR hard-codes `@HD VN:1.4` (samHeaders.cpp:94); rustar wrote noodles' default 1.6, so the first header line differed from STAR's in every SAM and BAM. A VN given in --outSAMheaderHD now sets the version instead of being stored as an extra field. Co-Authored-By: Claude Opus 5.5 --- src/io/sam.rs | 39 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/src/io/sam.rs b/src/io/sam.rs index 13ebc9f7..9c9f88f6 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -955,10 +955,14 @@ where // @HD line. `--outSAMheaderHD` replaces it wholesale, given as the // tab-separated fields STAR expects (`@HD VN:1.4 SO:coordinate`). + // STAR writes `VN:1.4` (`samHeaders.cpp:94`), not the latest spec version. + let star_version = sam::header::record::value::map::header::Version::new(1, 4); if params.out_sam_header_hd.is_empty() { - builder = builder.set_header(Map::default()); + builder = builder.set_header(Map::::new( + star_version, + )); } else { - let mut hd = Map::::default(); + let mut hd = Map::::new(star_version); for field in ¶ms.out_sam_header_hd { let field = field.trim(); // STAR takes the leading `@HD` as part of the value list; ignore it. @@ -975,6 +979,12 @@ where "--outSAMheaderHD tag '{tag}' is not two characters" ))); } + if tag == "VN" { + *hd.version_mut() = value + .parse() + .map_err(|e| Error::Parameter(format!("invalid @HD version '{value}': {e}")))?; + continue; + } 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}")))?; @@ -1769,6 +1779,31 @@ mod tests { assert_eq!(header.reference_sequences().len(), 1); } + #[test] + fn test_build_sam_header_hd_version_matches_star() { + use sam::header::record::value::map::header::Version; + let genome = make_test_genome(); + let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "test.fq"]); + let header = build_sam_header(&genome, ¶ms).unwrap(); + // STAR's default @HD line is `@HD VN:1.4` (samHeaders.cpp:94). + assert_eq!(header.header().unwrap().version(), Version::new(1, 4)); + + // A VN given in --outSAMheaderHD is the version, not an extra field. + let params = Parameters::parse_from([ + "rustar-aligner", + "--readFilesIn", + "test.fq", + "--outSAMheaderHD", + "@HD", + "VN:1.6", + "SO:unsorted", + ]); + let header = build_sam_header(&genome, ¶ms).unwrap(); + let hd = header.header().unwrap(); + assert_eq!(hd.version(), Version::new(1, 6)); + assert_eq!(hd.other_fields().len(), 1, "only SO besides VN"); + } + #[test] fn test_build_sam_header_pg_line_populated() { let genome = make_test_genome(); From 55a0cdae7c347afda1bc166b6b62e6e5c03648d9 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 8 Oct 2026 21:15:51 +0200 Subject: [PATCH 2/3] fix(sam): append SO:coordinate after the user's @HD fields; accept --outSAMheaderHD - STAR builds the sorted header as samHeaderHD + "\tSO:coordinate" (samHeaders.cpp:100), so a user @HD keeps its field order and SO comes last. '-' is STAR's default for --outSAMheaderHD and was rejected as a malformed field. Co-Authored-By: Claude Opus 5.5 --- src/io/bam.rs | 34 ++++++++++++++++++++++++++++++---- src/io/sam.rs | 20 +++++++++++++++++++- 2 files changed, 49 insertions(+), 5 deletions(-) diff --git a/src/io/bam.rs b/src/io/bam.rs index cad89eb0..a164965f 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -589,16 +589,18 @@ fn render_sam_text_lenient(header: &sam::Header, sort_order: Option<&str>) -> Ve if let Some(hd) = header.header() { buf.extend_from_slice(b"@HD\tVN:"); buf.extend_from_slice(hd.version().to_string().as_bytes()); - if let Some(so) = sort_order { - buf.extend_from_slice(b"\tSO:"); - buf.extend_from_slice(so.as_bytes()); - } for (tag, value) in hd.other_fields() { buf.push(b'\t'); buf.extend_from_slice(tag.as_ref()); buf.push(b':'); buf.extend_from_slice(value); } + // Appended after the user's fields, even when they carry their own SO, + // as STAR's `samHeaderSortedCoord` does (`samHeaders.cpp:100`). + if let Some(so) = sort_order { + buf.extend_from_slice(b"\tSO:"); + buf.extend_from_slice(so.as_bytes()); + } buf.push(b'\n'); } @@ -1146,4 +1148,28 @@ mod tests { "estimate must account for SEQ and QUAL" ); } + + #[test] + fn sorted_header_appends_so_after_the_user_hd_fields() { + let genome = create_test_genome(); + let render = |extra: &[&str]| { + let mut args = vec!["rustar-aligner", "--readFilesIn", "test.fq"]; + args.extend_from_slice(extra); + let params = Parameters::parse_from(args); + let header = crate::io::sam::build_sam_header(&genome, ¶ms).unwrap(); + let text = render_sam_text_lenient(&header, Some("coordinate")); + String::from_utf8(text) + .unwrap() + .lines() + .next() + .unwrap() + .to_string() + }; + // STAR: samHeaderHD + "\tSO:coordinate" (samHeaders.cpp:100). + assert_eq!(render(&[]), "@HD\tVN:1.4\tSO:coordinate"); + assert_eq!( + render(&["--outSAMheaderHD", "@HD", "VN:1.6", "GO:query"]), + "@HD\tVN:1.6\tGO:query\tSO:coordinate" + ); + } } diff --git a/src/io/sam.rs b/src/io/sam.rs index 9c9f88f6..5f819c74 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -957,7 +957,8 @@ where // tab-separated fields STAR expects (`@HD VN:1.4 SO:coordinate`). // STAR writes `VN:1.4` (`samHeaders.cpp:94`), not the latest spec version. let star_version = sam::header::record::value::map::header::Version::new(1, 4); - if params.out_sam_header_hd.is_empty() { + // `-` is STAR's default value for --outSAMheaderHD, meaning "no override". + if params.out_sam_header_hd.is_empty() || params.out_sam_header_hd == ["-"] { builder = builder.set_header(Map::::new( star_version, )); @@ -1804,6 +1805,23 @@ mod tests { assert_eq!(hd.other_fields().len(), 1, "only SO besides VN"); } + #[test] + fn test_build_sam_header_hd_dash_is_the_default() { + use sam::header::record::value::map::header::Version; + let genome = make_test_genome(); + let params = Parameters::parse_from([ + "rustar-aligner", + "--readFilesIn", + "test.fq", + "--outSAMheaderHD", + "-", + ]); + let header = build_sam_header(&genome, ¶ms).unwrap(); + let hd = header.header().unwrap(); + assert_eq!(hd.version(), Version::new(1, 4)); + assert_eq!(hd.other_fields().len(), 0); + } + #[test] fn test_build_sam_header_pg_line_populated() { let genome = make_test_genome(); From c855c8a189cfbf7d4b771f6eea28a98bec051bac Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Thu, 8 Oct 2026 21:50:09 +0200 Subject: [PATCH 3/3] Match STAR's SAM/BAM header text: @CO command line, verbatim comment file, user @HD, transcriptome BAM STAR's header layout (samHeaders.cpp): @HD as the user wrote it (VN:1.4 by default, no VN added to a user @HD), @SQ, user @PG, program @PG, comment-file lines verbatim (blank lines skipped), @RG, then @CO user command line. The transcriptome BAM carries only @SQ and @RG. The text is rendered once from the parameters and shared by the SAM and BAM writers. Co-Authored-By: Claude Opus 5.5 --- src/io/bam.rs | 178 +++++++++++++++++++------------------------ src/io/sam.rs | 110 +++++++++++++++++++++++--- src/params/mod.rs | 20 +++++ tests/determinism.rs | 4 +- 4 files changed, 199 insertions(+), 113 deletions(-) diff --git a/src/io/bam.rs b/src/io/bam.rs index a164965f..60604089 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -1,6 +1,7 @@ /// BAM output writer with noodles (streaming, unsorted) use crate::error::Error; use crate::genome::Genome; +use crate::io::sam::HeaderKind; use crate::params::Parameters; use crate::quant::transcriptome::TranscriptomeIndex; use byteorder::{LittleEndian, WriteBytesExt}; @@ -167,6 +168,7 @@ struct CoordinateSorter { buffered_bytes: u64, runs: Vec, header: sam::Header, + header_text: Vec, compression: i32, ram_limit: u64, temp_dir: std::path::PathBuf, @@ -174,8 +176,9 @@ struct CoordinateSorter { } impl CoordinateSorter { - fn new(header: sam::Header, params: &Parameters) -> Self { + fn new(header: sam::Header, header_text: Vec, params: &Parameters) -> Self { Self { + header_text, records: Vec::new(), buffered_bytes: 0, runs: Vec::new(), @@ -257,7 +260,7 @@ impl CoordinateSorter { self.records.sort_by_key(sort_key); let mut bgzf = make_bgzf_writer(out, self.compression); - write_bam_header_lenient(&mut bgzf, &self.header, Some("coordinate"))?; + write_bam_header_lenient(&mut bgzf, &self.header, &self.header_text)?; let mut writer = bam::io::Writer::from(bgzf); let written = if self.runs.is_empty() { @@ -429,22 +432,21 @@ impl BamWriter { fn with_header( output_path: &Path, header: sam::Header, + header_text: &[u8], compression: i32, ) -> Result { let buf_writer = BufWriter::new(File::create(output_path)?); let mut bgzf = make_bgzf_writer(buf_writer, compression); - write_bam_header_lenient(&mut bgzf, &header, None)?; + write_bam_header_lenient(&mut bgzf, &header, header_text)?; let writer = bam::io::Writer::from(bgzf); Ok(Self { writer, header }) } /// Create a new BAM writer with header from genome index. pub fn create(output_path: &Path, genome: &Genome, params: &Parameters) -> Result { - Self::with_header( - output_path, - crate::io::sam::build_sam_header(genome, params)?, - params.out_bam_compression, - ) + let header = crate::io::sam::build_sam_header(genome, params)?; + let text = crate::io::sam::render_header_text(&header, params, HeaderKind::Main)?; + Self::with_header(output_path, header, &text, params.out_bam_compression) } /// Create a BAM writer whose @SQ header lists the transcripts (not @@ -460,11 +462,9 @@ impl BamWriter { .iter() .zip(tr_idx.tr_length.iter()) .map(|(id, len)| (id.as_str(), *len as usize)); - Self::with_header( - output_path, - crate::io::sam::build_sam_header_from_refs(refs, params)?, - params.out_bam_compression, - ) + let header = crate::io::sam::build_sam_header_from_refs(refs, params)?; + let text = crate::io::sam::render_header_text(&header, params, HeaderKind::Transcriptome)?; + Self::with_header(output_path, header, &text, params.out_bam_compression) } /// Write batch of buffered records (for parallel processing) @@ -495,8 +495,9 @@ impl SortedBamWriter { params: &Parameters, ) -> Result { let header = crate::io::sam::build_sam_header(genome, params)?; + let text = crate::io::sam::render_header_text(&header, params, HeaderKind::MainSorted)?; Ok(Self { - sorter: CoordinateSorter::new(header, params), + sorter: CoordinateSorter::new(header, text, params), output_path: output_path.to_path_buf(), }) } @@ -530,22 +531,18 @@ impl SortedBamWriter { /// only constraint there is "no interior nul", which is enforced upstream /// via the usual UTF-8 input. /// -/// `sort_order`: if `Some("coordinate")`, injects `SO:coordinate` into the -/// @HD line. Pass `None` for unsorted output. +/// `text`: the SAM header text, from `io::sam::render_header_text`. fn write_bam_header_lenient( writer: &mut W, header: &sam::Header, - sort_order: Option<&str>, + text: &[u8], ) -> Result<(), Error> { const MAGIC: &[u8; 4] = b"BAM\x01"; writer.write_all(MAGIC)?; - // Build the SAM text block byte-for-byte identical to - // `sam::io::Writer::write_header` minus the name validator: - // `@HD`, `@SQ` (one per reference), `@RG` (if any), `@PG` (if any), - // `@CO` (if any), each line terminated by `\n`. - let text = render_sam_text_lenient(header, sort_order); + // The SAM text block is rendered by `io::sam::render_header_text` in STAR's + // line order, without noodles' name validator. let l_text = i32::try_from(text.len()).map_err(|_| { Error::Index(format!( "BAM SAM-text header exceeds i32::MAX bytes: {} bytes", @@ -553,7 +550,7 @@ fn write_bam_header_lenient( )) })?; writer.write_i32::(l_text)?; - writer.write_all(&text)?; + writer.write_all(text)?; // Binary reference list: n_ref then (l_name, name\0, l_ref) per ref. let refs = header.reference_sequences(); @@ -580,79 +577,6 @@ fn write_bam_header_lenient( Ok(()) } -fn render_sam_text_lenient(header: &sam::Header, sort_order: Option<&str>) -> Vec { - let mut buf: Vec = Vec::new(); - - // @HD line. noodles' Map
in this version doesn't expose - // sort_order/group_order via dedicated accessors; we serialize through - // the generic `other_fields` map (which includes SO/GO when set). - if let Some(hd) = header.header() { - buf.extend_from_slice(b"@HD\tVN:"); - buf.extend_from_slice(hd.version().to_string().as_bytes()); - for (tag, value) in hd.other_fields() { - buf.push(b'\t'); - buf.extend_from_slice(tag.as_ref()); - buf.push(b':'); - buf.extend_from_slice(value); - } - // Appended after the user's fields, even when they carry their own SO, - // as STAR's `samHeaderSortedCoord` does (`samHeaders.cpp:100`). - if let Some(so) = sort_order { - buf.extend_from_slice(b"\tSO:"); - buf.extend_from_slice(so.as_bytes()); - } - buf.push(b'\n'); - } - - // @SQ lines. Use the raw bytes so forbidden characters pass through. - for (name, rs) in header.reference_sequences() { - buf.extend_from_slice(b"@SQ\tSN:"); - buf.extend_from_slice(name); - buf.extend_from_slice(b"\tLN:"); - buf.extend_from_slice(usize::from(rs.length()).to_string().as_bytes()); - // Other optional @SQ fields (AH, AN, AS, DS, M5, SP, TP, UR) — - // rustar-aligner doesn't set any today, so skip. - buf.push(b'\n'); - } - - // @RG lines. - for (id, rg) in header.read_groups() { - buf.extend_from_slice(b"@RG\tID:"); - buf.extend_from_slice(id); - for (tag, value) in rg.other_fields() { - buf.push(b'\t'); - buf.extend_from_slice(tag.as_ref()); - buf.push(b':'); - buf.extend_from_slice(value); - } - buf.push(b'\n'); - } - - // @PG lines — noodles' map doesn't guarantee insertion order; for - // rustar-aligner we emit a single @PG with id "rustar-aligner". If more are added - // later, pipe them in here. - for (id, pg) in header.programs().as_ref() { - buf.extend_from_slice(b"@PG\tID:"); - buf.extend_from_slice(id); - for (tag, value) in pg.other_fields() { - buf.push(b'\t'); - buf.extend_from_slice(tag.as_ref()); - buf.push(b':'); - buf.extend_from_slice(value); - } - buf.push(b'\n'); - } - - // @CO lines (comments). - for comment in header.comments() { - buf.extend_from_slice(b"@CO\t"); - buf.extend_from_slice(comment); - buf.push(b'\n'); - } - - buf -} - /// Streaming unsorted BAM writer that writes to stdout. pub struct BamStdoutWriter { writer: bam::io::Writer>>, @@ -662,11 +586,12 @@ pub struct BamStdoutWriter { impl BamStdoutWriter { pub fn create(genome: &crate::genome::Genome, params: &Parameters) -> Result { let header = crate::io::sam::build_sam_header(genome, params)?; + let header_text = crate::io::sam::render_header_text(&header, params, HeaderKind::Main)?; let mut bgzf = make_bgzf_writer( BufWriter::new(std::io::stdout()), params.out_bam_compression, ); - write_bam_header_lenient(&mut bgzf, &header, None)?; + write_bam_header_lenient(&mut bgzf, &header, &header_text)?; let writer = bam::io::Writer::from(bgzf); Ok(Self { writer, header }) } @@ -693,8 +618,9 @@ pub struct SortedBamStdoutWriter { impl SortedBamStdoutWriter { pub fn create(genome: &crate::genome::Genome, params: &Parameters) -> Result { let header = crate::io::sam::build_sam_header(genome, params)?; + let text = crate::io::sam::render_header_text(&header, params, HeaderKind::MainSorted)?; Ok(Self { - sorter: CoordinateSorter::new(header, params), + sorter: CoordinateSorter::new(header, text, params), }) } @@ -1157,7 +1083,8 @@ mod tests { args.extend_from_slice(extra); let params = Parameters::parse_from(args); let header = crate::io::sam::build_sam_header(&genome, ¶ms).unwrap(); - let text = render_sam_text_lenient(&header, Some("coordinate")); + let text = crate::io::sam::render_header_text(&header, ¶ms, HeaderKind::MainSorted) + .unwrap(); String::from_utf8(text) .unwrap() .lines() @@ -1172,4 +1099,57 @@ mod tests { "@HD\tVN:1.6\tGO:query\tSO:coordinate" ); } + + #[test] + fn header_text_follows_stars_layout() { + let genome = create_test_genome(); + let dir = std::env::temp_dir().join(format!("hdr_text_{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let com = dir.join("com.txt"); + std::fs::write(&com, "@CO\tfirst\n\n \n@CO\tsecond\n").unwrap(); + let com = com.to_str().unwrap().to_string(); + let params = Parameters::parse_from([ + "rustar-aligner", + "--readFilesIn", + "test.fq", + "--outSAMattrRGline", + "ID:x", + "SM:y z", + "--outSAMheaderPG", + "@PG", + "ID:foo", + "--outSAMheaderCommentFile", + &com, + "--outSAMheaderHD", + "@HD", + "SO:foo", + ]); + let header = crate::io::sam::build_sam_header(&genome, ¶ms).unwrap(); + let text = |k| { + String::from_utf8(crate::io::sam::render_header_text(&header, ¶ms, k).unwrap()) + .unwrap() + }; + let main = text(HeaderKind::Main); + let tags: Vec<&str> = main.lines().map(|l| &l[..3]).collect(); + // User @HD verbatim (no VN), @SQ, user @PG, @PG, comments, @RG, @CO. + assert_eq!(main.lines().next().unwrap(), "@HD\tSO:foo"); + let pos = |t: &str| tags.iter().position(|x| *x == t).unwrap(); + assert!(pos("@SQ") < pos("@PG") && pos("@PG") < pos("@RG") && pos("@RG") < tags.len() - 1); + assert_eq!( + main.lines().nth(tags.len() - 1).unwrap()[..21].to_string(), + "@CO\tuser command line" + ); + assert!(main.contains("@PG\tID:foo\n")); + // Comment file lines are copied as they are, blank lines skipped. + assert!(main.contains("\n@CO\tfirst\n@CO\tsecond\n@RG")); + assert!(!main.contains("@CO\t@CO")); + // The transcriptome BAM header has no @HD, @PG or @CO. + let tr = text(HeaderKind::Transcriptome); + assert!( + tr.lines() + .all(|l| l.starts_with("@SQ") || l.starts_with("@RG")) + ); + assert!(tr.contains("@RG\tID:x")); + std::fs::remove_dir_all(&dir).ok(); + } } diff --git a/src/io/sam.rs b/src/io/sam.rs index 5f819c74..e1265fc4 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -24,7 +24,7 @@ use noodles::sam::header::record::value::{ use std::collections::HashSet; use std::fmt::Write as FmtWrite; use std::fs::File; -use std::io::BufWriter; +use std::io::{BufWriter, Write as _}; use std::num::NonZeroUsize; use std::path::Path; @@ -98,7 +98,9 @@ impl SamWriter { let header = build_sam_header(genome, params)?; let mut writer = sam::io::Writer::new(buf_writer); - writer.write_header(&header)?; + writer + .get_mut() + .write_all(&render_header_text(&header, params, HeaderKind::Main)?)?; Ok(Self { writer, header }) } @@ -875,7 +877,9 @@ impl SamStdoutWriter { pub fn create(genome: &Genome, params: &Parameters) -> Result { let header = build_sam_header(genome, params)?; let mut writer = sam::io::Writer::new(BufWriter::new(std::io::stdout())); - writer.write_header(&header)?; + writer + .get_mut() + .write_all(&render_header_text(&header, params, HeaderKind::Main)?)?; Ok(Self { writer, header }) } @@ -1083,19 +1087,101 @@ where 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(¶ms.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()) +} + +/// Which SAM header text to render. +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +pub enum HeaderKind { + /// `Aligned.out.sam` / unsorted BAM. + Main, + /// Coordinate-sorted BAM: `SO:coordinate` appended to the `@HD` line. + MainSorted, + /// `Aligned.toTranscriptome.out.bam`: `@SQ` (transcripts) and `@RG` only. + Transcriptome, +} + +/// Render the SAM header text in STAR's line order (`samHeaders.cpp`): +/// `@HD`, `@SQ`, user `@PG` (`--outSAMheaderPG`), the program `@PG`, the +/// `--outSAMheaderCommentFile` lines verbatim, `@RG`, then +/// `@CO user command line: ...`. The `@HD` line is the user's +/// `--outSAMheaderHD` text verbatim (no `VN` is added), else `@HD\tVN:1.4`. +/// The transcriptome BAM carries only `@SQ` and `@RG` lines. +pub fn render_header_text( + header: &sam::Header, + params: &Parameters, + kind: HeaderKind, +) -> Result, Error> { + let mut buf: Vec = Vec::new(); + + if kind != HeaderKind::Transcriptome { + if params.out_sam_header_hd.is_empty() || params.out_sam_header_hd[0] == "-" { + buf.extend_from_slice(b"@HD\tVN:1.4"); + } else { + buf.extend_from_slice(params.out_sam_header_hd.join("\t").as_bytes()); + } + if kind == HeaderKind::MainSorted { + buf.extend_from_slice(b"\tSO:coordinate"); + } + buf.push(b'\n'); + } + + // Raw bytes so reference names with characters noodles forbids pass through. + for (name, rs) in header.reference_sequences() { + buf.extend_from_slice(b"@SQ\tSN:"); + buf.extend_from_slice(name); + buf.extend_from_slice(b"\tLN:"); + buf.extend_from_slice(usize::from(rs.length()).to_string().as_bytes()); + buf.push(b'\n'); + } + + if kind != HeaderKind::Transcriptome { + if !params.out_sam_header_pg.is_empty() && params.out_sam_header_pg[0] != "-" { + buf.extend_from_slice(params.out_sam_header_pg.join("\t").as_bytes()); + buf.push(b'\n'); + } + if let Some(pg) = header.programs().as_ref().get(&b"rustar-aligner"[..]) { + buf.extend_from_slice(b"@PG\tID:rustar-aligner"); + for (tag, value) in pg.other_fields() { + buf.push(b'\t'); + buf.extend_from_slice(tag.as_ref()); + buf.push(b':'); + buf.extend_from_slice(value); + } + buf.push(b'\n'); + } + if params.out_sam_header_comment_file != "-" { + let path = std::path::Path::new(¶ms.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() { + // STAR copies non-blank lines as they are (each should start with @CO). + if !line.trim().is_empty() { + buf.extend_from_slice(line.as_bytes()); + buf.push(b'\n'); + } } } } - Ok(builder.build()) + for (id, rg) in header.read_groups() { + buf.extend_from_slice(b"@RG\tID:"); + buf.extend_from_slice(id); + for (tag, value) in rg.other_fields() { + buf.push(b'\t'); + buf.extend_from_slice(tag.as_ref()); + buf.push(b':'); + buf.extend_from_slice(value); + } + buf.push(b'\n'); + } + + if kind != HeaderKind::Transcriptome { + buf.extend_from_slice(b"@CO\tuser command line: "); + buf.extend_from_slice(params.star_command_line.as_deref().unwrap_or("").as_bytes()); + buf.push(b'\n'); + } + + Ok(buf) } /// Insert `RG:Z:` on the record when an ID is set. `Parameters::try_parse_from` diff --git a/src/params/mod.rs b/src/params/mod.rs index a3539e08..067d8807 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1346,6 +1346,12 @@ pub struct Parameters { /// Full command line as invoked, embedded in the BAM `@PG` `CL:` field. #[arg(skip)] pub command_line: Option, + + /// The command line as STAR records it (`Parameters.cpp`, `commandLine`): + /// the arguments joined by single spaces, a value argument containing + /// white space wrapped in double quotes. Written to the `@CO` header line. + #[arg(skip)] + pub star_command_line: Option, } impl Parameters { @@ -1546,6 +1552,20 @@ impl Parameters { let matches = command.clone().get_matches_from(args.iter()); let mut params = ::from_arg_matches(&matches)?; + params.star_command_line = Some( + args.iter() + .enumerate() + .map(|(i, a)| { + let a = a.to_string_lossy(); + if i > 0 && !a.starts_with("--") && a.contains([' ', '\t']) { + format!("\"{a}\"") + } else { + a.into_owned() + } + }) + .collect::>() + .join(" "), + ); params.command_line = { let args: Vec<_> = args.iter().map(|s| s.to_string_lossy()).collect(); shlex::try_join(args.iter().map(AsRef::as_ref)).ok() diff --git a/tests/determinism.rs b/tests/determinism.rs index a7787481..55917ab9 100644 --- a/tests/determinism.rs +++ b/tests/determinism.rs @@ -132,12 +132,12 @@ fn align( .iter() .map(|name| { let bytes = fs::read(format!("{prefix}{name}")).unwrap_or_default(); - // Drop the @PG header line: it records the command line, which + // Drop the @PG and @CO header lines: they record the command line, which // differs by the thread count itself. let filtered: Vec = if *name == "Aligned.out.sam" { String::from_utf8_lossy(&bytes) .lines() - .filter(|l| !l.starts_with("@PG")) + .filter(|l| !l.starts_with("@PG") && !l.starts_with("@CO")) .collect::>() .join("\n") .into_bytes()