From 0232413cda140f88d9e8e1d88fdf5e82484e45a5 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 11 Aug 2026 22:20:22 +0200 Subject: [PATCH 1/3] feat(solo): STARsolo per-read barcode and gene SAM tags (Phase 14.7) `--outSAMattributes` now accepts `CR CY UR UY CB UB gx gn sM sS sQ sF` alongside the existing `GX GN`, emitted in BAM output only as STAR does (`ReadAlign_outputTranscriptSAM.cpp` drops them from the SAM text path). Read-time tags (raw barcode/UMI with qualities, the `sM` cbMatch code, the whole barcode read in `sS`/`sQ`) are written on every record of a read, mapped or unmapped, in the single-end and paired-end solo loops. `CB`/`UB` follow STAR's readInfo route: count records carry their input read index, UMI collapsing records what each read was counted as, and the buffered sorted-BAM records are rewritten before the sort, with `-` for reads that were not counted. Validation enforces STAR's rules: sorted BAM output, a gene-level first `--soloFeatures` entry, and no `UB` under `CB_samTagOut`. `gx`/`gn` list every gene of their own alignment and `sF` carries `(overlap type, genes for the read)`, from `align_genes`, which keeps STAR's per-alignment gene sets and `ovType` priority. A mate pair counts as one alignment, so both mates carry the pair's genes. Also in this change: - `--soloType CB_samTagOut`: whitelist correction into `CB` with no gene model, no UMI collapsing and no `Solo.out`. - `--readFilesIn` accepts STAR's three-file solo layout (`cDNA_read1 cDNA_read2 barcode_read`), so paired-end cDNA works with a separate barcode read; a third file without `--soloType` is refused. - `--soloBarcodeReadLength` is honoured: the default requires the barcode read to be exactly CB+UMI long and treats any other length as a fatal input error, while `0` disables the check and pads a short read with `N`/`H`. - `--soloUMIdedup 1MM_Directional`/`1MM_Directional_UMItools` count distinct corrected UMIs off STAR's absorb chain (`umiArrayCorrect_Directional`) rather than counting unabsorbed UMIs. Co-Authored-By: Claude Opus 5 (1M context) --- src/io/bam.rs | 11 + src/io/sam.rs | 229 +++++++++++++++++ src/lib.rs | 253 +++++++++++++++++-- src/params/mod.rs | 306 +++++++++++++++++++++- src/params/sam.rs | 60 ++++- src/quant/mod.rs | 4 + src/solo/count.rs | 238 +++++++++++++++-- src/solo/gene.rs | 182 +++++++++++++ src/solo/mod.rs | 576 ++++++++++++++++++++++++++++++++++++++---- src/solo/whitelist.rs | 36 +++ 10 files changed, 1792 insertions(+), 103 deletions(-) diff --git a/src/io/bam.rs b/src/io/bam.rs index b1e6717f..92d86941 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -164,6 +164,12 @@ impl SortedBamWriter { Ok(()) } + /// The buffered records, for post-processing before the sort (STARsolo's + /// `CB`/`UB` tags, which are only known after the counting pass). + pub fn records_mut(&mut self) -> &mut Vec { + &mut self.records + } + /// Estimate memory used by buffered records (rough: 400 bytes/record for 150bp reads). fn estimated_ram(&self) -> u64 { self.records.len() as u64 * 400 @@ -424,6 +430,11 @@ impl SortedBamStdoutWriter { Ok(()) } + /// The buffered records, for STARsolo's post-counting `CB`/`UB` pass. + pub fn records_mut(&mut self) -> &mut Vec { + &mut self.records + } + pub fn finish(&mut self) -> Result<(), Error> { if self.limit_bam_sort_ram > 0 { let est = self.records.len() as u64 * 400; diff --git a/src/io/sam.rs b/src/io/sam.rs index 189c4bb5..87215706 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -1112,6 +1112,179 @@ pub fn add_gene_tags(records: &mut [RecordBuf], gx: &str, gn: &str, attrs: SamAt } } +/// The per-read STARsolo barcode tag values of one read +/// (`ReadAlign_alignBAM.cpp:389-473`). +/// +/// `CB` is filled here only when the barcode is corrected as the read is +/// processed (`--soloType CB_samTagOut`); in a counting run it, and `UB`: are +/// added when the sorted BAM is written, once UMI collapsing has run. +pub struct SoloBarcodeTagValues<'a> { + /// `CR`: raw cell barcode. + pub cb_seq: &'a str, + /// `CY`: raw cell-barcode quality. + pub cb_qual: &'a str, + /// `UR`: raw UMI. + pub umi_seq: &'a str, + /// `UY`: raw UMI quality. + pub umi_qual: &'a str, + /// `CB`: corrected cell barcode, when known at read time. + pub cb_corrected: Option<&'a str>, + /// `sM`: STAR's `cbMatch` code. + pub cb_match: i32, + /// `sS`: whole barcode read. + pub barcode_seq: &'a str, + /// `sQ`: whole barcode-read quality. + pub barcode_qual: &'a str, +} + +/// Add the STARsolo per-read barcode tags to every record of one read (mapped +/// alignments and the unmapped record alike, as in STAR). `attrs` should already +/// be narrowed to the tags this run emits, see `Parameters::solo_sam_tags`. +pub fn add_solo_barcode_tags( + records: &mut [RecordBuf], + values: &SoloBarcodeTagValues<'_>, + attrs: SamAttributes, +) { + if !attrs.intersects(SamAttributes::SOLO_TAGS) { + return; + } + let str_tags: [(SamAttributes, [u8; 2], &str); 5] = [ + (SamAttributes::CR, *b"CR", values.cb_seq), + (SamAttributes::CY, *b"CY", values.cb_qual), + (SamAttributes::UR, *b"UR", values.umi_seq), + (SamAttributes::UY, *b"UY", values.umi_qual), + (SamAttributes::SS, *b"sS", values.barcode_seq), + ]; + for rec in records.iter_mut() { + // An empty value means the read carries no such sequence at all (a + // barcode read too short to hold a CB+UMI); STAR has no empty tags, so + // the tag is left off rather than written blank. + for (flag, tag, value) in str_tags { + if attrs.contains(flag) && !value.is_empty() { + rec.data_mut().insert( + Tag::new(tag[0], tag[1]), + Value::String(BString::from(value)), + ); + } + } + // sQ is a quality string, written verbatim like sS. + if attrs.contains(SamAttributes::SQ) && !values.barcode_qual.is_empty() { + rec.data_mut().insert( + Tag::new(b's', b'Q'), + Value::String(BString::from(values.barcode_qual)), + ); + } + if attrs.contains(SamAttributes::SM) { + rec.data_mut() + .insert(Tag::new(b's', b'M'), Value::Int32(values.cb_match)); + } + if let (true, Some(cb)) = (attrs.contains(SamAttributes::CB), values.cb_corrected) { + rec.data_mut() + .insert(Tag::new(b'C', b'B'), Value::String(BString::from(cb))); + } + } +} + +/// Add the per-alignment STARsolo gene tags (`gx`, `gn`, `sF`) to a read's +/// records. +/// +/// `tags` is one entry per alignment; `records_per_align` is how many records +/// each alignment produced (1 single-end, 2 for a mate pair), so both mates of a +/// pair carry the pair's genes. Records past the end of `tags` (there are none +/// in practice) are left alone. +pub fn add_align_gene_tags( + records: &mut [RecordBuf], + tags: &[crate::solo::AlignGeneTag], + records_per_align: usize, + attrs: SamAttributes, +) { + use noodles::sam::alignment::record_buf::data::field::value::Array; + + if !attrs.intersects(SamAttributes::GXM | SamAttributes::GNM | SamAttributes::SF) + || records_per_align == 0 + { + return; + } + for (i, rec) in records.iter_mut().enumerate() { + let Some(tag) = tags.get(i / records_per_align) else { + continue; + }; + if attrs.contains(SamAttributes::GXM) { + rec.data_mut().insert( + Tag::new(b'g', b'x'), + Value::String(BString::from(tag.gx.as_str())), + ); + } + if attrs.contains(SamAttributes::GNM) { + rec.data_mut().insert( + Tag::new(b'g', b'n'), + Value::String(BString::from(tag.gn.as_str())), + ); + } + if attrs.contains(SamAttributes::SF) { + rec.data_mut().insert( + Tag::new(b's', b'F'), + Value::Array(Array::Int32(tag.sf.to_vec())), + ); + } + } +} + +/// Private aux tag carrying the input read index on a buffered record until the +/// sorted BAM is written, where it becomes `CB`/`UB`. STAR does the same thing +/// by encoding `iReadAll` in the record's trailing bytes +/// (`SoloFeature_addBAMtags.cpp:8`); a local `z`-namespace tag survives sorting +/// without a parallel array. +const SOLO_READ_INDEX_TAG: [u8; 2] = *b"zR"; + +/// Stamp the input read index on each of a read's records, for the CB/UB pass. +pub fn add_solo_read_index(records: &mut [RecordBuf], read_index: u32) { + for rec in records.iter_mut() { + rec.data_mut().insert( + Tag::new(SOLO_READ_INDEX_TAG[0], SOLO_READ_INDEX_TAG[1]), + Value::UInt32(read_index), + ); + } +} + +/// Replace the private read-index tag with `CB`/`UB`, read out of STAR's +/// readInfo once UMI collapsing has run. +/// +/// Both tags are written whenever either was requested, and both fall back to +/// `"-"`, exactly as `SoloFeature::addBAMtags` does: a read that was not counted +/// (no whitelist cell, no valid UMI, no gene) has no cell or molecule to name. +pub fn apply_solo_read_info( + records: &mut [RecordBuf], + read_info: &[crate::solo::ReadInfo], + whitelist: &crate::solo::CbWhitelist, + umi_len: usize, +) { + let tag = Tag::new(SOLO_READ_INDEX_TAG[0], SOLO_READ_INDEX_TAG[1]); + for rec in records.iter_mut() { + let Some(Value::UInt32(read_index)) = rec.data().get(&tag).cloned() else { + continue; + }; + rec.data_mut().remove(&tag); + let info = read_info + .get(read_index as usize) + .copied() + .unwrap_or_default(); + let cb = (info.cb != u32::MAX) + .then(|| whitelist.barcode_string(info.cb)) + .flatten() + .unwrap_or_else(|| "-".to_string()); + let ub = if info.umi == u64::MAX { + "-".to_string() + } else { + crate::solo::whitelist::unpack_barcode(info.umi, umi_len) + }; + rec.data_mut() + .insert(Tag::new(b'C', b'B'), Value::String(BString::from(cb))); + rec.data_mut() + .insert(Tag::new(b'U', b'B'), Value::String(BString::from(ub))); + } +} + /// Apply `--outSAMflagOR` / `--outSAMflagAND` to a mapped record's FLAG: /// `(FLAG & flagAND) | flagOR`. Matches STAR/STAR-rs, which apply this only to /// mapped-mate records; unmapped and transcriptome-BAM records are untouched. @@ -1752,6 +1925,62 @@ mod tests { } } + /// The solo barcode tags land on every record of the read, and only the + /// requested ones are written. Values with no sequence behind them (a + /// barcode read too short to hold a CB+UMI) are left off entirely. + #[test] + fn solo_barcode_tags_are_added_per_requested_attribute() { + let values = SoloBarcodeTagValues { + cb_seq: "ACGTACGTACGTACGT", + cb_qual: "IIIIIIIIIIIIIIII", + umi_seq: "ACGTACGTAC", + umi_qual: "JJJJJJJJJJ", + cb_corrected: Some("ACGTACGTACGTACGA"), + cb_match: 1, + barcode_seq: "", + barcode_qual: "", + }; + let mut records = vec![RecordBuf::default(), RecordBuf::default()]; + let attrs = SamAttributes::CR + | SamAttributes::CY + | SamAttributes::UR + | SamAttributes::SM + | SamAttributes::SS + | SamAttributes::CB; + add_solo_barcode_tags(&mut records, &values, attrs); + + for rec in &records { + let data = rec.data(); + let get = |t: [u8; 2]| data.get(&Tag::new(t[0], t[1])).cloned(); + assert_eq!( + get(*b"CR"), + Some(Value::String(BString::from("ACGTACGTACGTACGT"))) + ); + assert_eq!( + get(*b"CY"), + Some(Value::String(BString::from("IIIIIIIIIIIIIIII"))) + ); + assert_eq!( + get(*b"UR"), + Some(Value::String(BString::from("ACGTACGTAC"))) + ); + assert_eq!(get(*b"sM"), Some(Value::Int32(1))); + assert_eq!( + get(*b"CB"), + Some(Value::String(BString::from("ACGTACGTACGTACGA"))) + ); + // Not requested, and no value: absent. + assert_eq!(get(*b"UY"), None); + assert_eq!(get(*b"sS"), None); + assert_eq!(get(*b"sQ"), None); + } + + // No solo attributes requested → untouched records. + let mut untouched = vec![RecordBuf::default()]; + add_solo_barcode_tags(&mut untouched, &values, SamAttributes::STANDARD); + assert!(untouched[0].data().is_empty()); + } + #[test] fn test_build_sam_header() { let genome = make_test_genome(); diff --git a/src/lib.rs b/src/lib.rs index 5086fcbb..d3477e27 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -154,6 +154,16 @@ trait AlignmentWriter: Send { fn finish(&mut self) -> Result<(), error::Error> { Ok(()) } + /// Records still held in memory before `finish` sorts and writes them. + /// + /// Only the coordinate-sorted writers have any: STARsolo's `CB`/`UB` tags + /// are filled here, after the counting pass, which is exactly why STAR + /// restricts them to sorted BAM output. + fn buffered_records_mut( + &mut self, + ) -> Option<&mut Vec> { + None + } } /// Drops quality strings on the way to the real writer (`--outSAMmode NoQS`). @@ -179,6 +189,12 @@ impl AlignmentWriter for NoQsWriter { fn finish(&mut self) -> Result<(), error::Error> { self.0.finish() } + + fn buffered_records_mut( + &mut self, + ) -> Option<&mut Vec> { + self.0.buffered_records_mut() + } } /// Null writer that discards all output (for two-pass mode pass 1) @@ -224,6 +240,11 @@ impl AlignmentWriter for crate::io::bam::SortedBamWriter { fn finish(&mut self) -> Result<(), error::Error> { self.finish() } + fn buffered_records_mut( + &mut self, + ) -> Option<&mut Vec> { + Some(self.records_mut()) + } } impl AlignmentWriter for crate::io::sam::SamStdoutWriter { @@ -257,6 +278,11 @@ impl AlignmentWriter for crate::io::bam::SortedBamStdoutWriter { fn finish(&mut self) -> Result<(), error::Error> { self.finish() } + fn buffered_records_mut( + &mut self, + ) -> Option<&mut Vec> { + Some(self.records_mut()) + } } fn align_reads(params: &Parameters) -> anyhow::Result<()> { @@ -780,14 +806,14 @@ fn run_single_pass( // The dedicated solo loop reads the barcode read in lockstep, quantifies // per cell, and otherwise emits the cDNA alignments like the SE path. if let Some(sctx) = solo_ctx { - // `--soloBarcodeMate 1` (5' 10x): barcode on mate 1, both mates aligned as - // a pair. Otherwise the standard SE-solo path (barcode on a separate read). - if params.solo_barcode_on_mate1() { + // Paired-end cDNA: either `--soloBarcodeMate 1` (5' 10x, barcode on mate + // 1) or a third `--readFilesIn` barcode file. Otherwise the single-end + // solo path (one cDNA read + a barcode read). + if params.solo_paired_cdna() { align_reads_solo_pe(params, index, writer.as_mut(), &stats, &sj_stats, sctx)?; } else { align_reads_solo(params, index, writer.as_mut(), &stats, &sj_stats, sctx)?; } - writer.finish()?; if let Some(ref mut w) = tr_writer { w.finish()?; } @@ -798,7 +824,32 @@ fn run_single_pass( } // Per-cell count matrices (raw + filtered), Summary.csv, and the SJ // feature matrix — written here where sj_stats is available. - write_solo_output(sctx, params, &stats, &sj_stats, index)?; + // + // With `CB`/`UB` requested this has to precede the alignment output: + // both tags come out of the readInfo that UMI collapsing fills, which is + // why STAR only offers them for the (post-counting) sorted BAM. + sctx.reserve_read_info(stats.total_reads() as usize); + // `CB_samTagOut` produces no Solo.out matrices at all (`Solo.cpp:13`). + if params.solo_type != params::SoloType::CbSamTagOut { + write_solo_output(sctx, params, &stats, &sj_stats, index)?; + } + if sctx.read_info_enabled() + && let Some(records) = writer.buffered_records_mut() + { + let info = sctx + .read_info + .as_ref() + .expect("readInfo enabled") + .lock() + .unwrap(); + crate::io::sam::apply_solo_read_info( + records, + &info, + &sctx.whitelist, + params.solo_umi_len as usize, + ); + } + writer.finish()?; stats.print_summary(); return Ok(stats); } @@ -2078,6 +2129,15 @@ fn align_reads_solo( // — a large saving for solo runs that only need the count matrix. let emit_sam = params.emits_alignments(); let output_unmapped = emit_sam && params.out_sam_unmapped != params::OutSamUnmapped::None; + // STARsolo barcode tags requested for this run (BAM output only, as in STAR). + let solo_tags = if emit_sam { + params.solo_sam_tags() + } else { + crate::params::SamAttributes::empty() + }; + let track_read_info = emit_sam && solo_ctx.read_info_enabled(); + // `--soloType CB_samTagOut`: barcodes become SAM tags, nothing is counted. + let tag_out_only = params.solo_type == params::SoloType::CbSamTagOut; // Shared, 'static parameters for the per-batch aligner tasks spawned below. let params_arc = Arc::new(params.clone()); @@ -2260,13 +2320,35 @@ fn align_reads_solo( Vec::new() }; - // Solo quantification (CB match + UMI check + gene assignment). - let outcome = solo.process_read( - &transcripts, - transcripts.len(), - sread.barcode.as_ref(), - &junctions, - ); + // Solo quantification (CB match + UMI check + gene + // assignment). `CB_samTagOut` skips all of it: the + // barcode is corrected for the CB tag and nothing is + // counted. + let (mut outcome, cb_corrected) = if tag_out_only { + let mut outcome = crate::solo::SoloReadOutcome::default(); + let corrected = sread.barcode.as_ref().map(|bc| { + let (tags, corrected) = solo.tag_barcode(bc); + outcome.barcode = Some(tags); + corrected + }); + (outcome, corrected) + } else { + ( + solo.process_read( + &transcripts, + transcripts.len(), + sread.barcode.as_ref(), + &junctions, + ), + None, + ) + }; + // CB/UB: tie this read's count records to its input + // index so collapsing can fill STAR's readInfo. + let read_index = (base + read_idx as u64) as u32; + if track_read_info { + outcome.set_read_index(read_index); + } // Build SAM records for the cDNA alignment (same as SE path). // Skipped entirely under `--outSAMtype None` (count-only). @@ -2299,7 +2381,7 @@ fn align_reads_solo( n_for_mapq, )?; // STARsolo GX/GN gene tags (Gene-feature assignment). - if params.out_sam_attributes.intersects( + if solo_tags.intersects( crate::params::SamAttributes::GX | crate::params::SamAttributes::GN, ) { @@ -2308,7 +2390,21 @@ fn align_reads_solo( &mut records, gx, gn, - params.out_sam_attributes, + solo_tags, + ); + } + // Per-alignment gx/gn + the sF feature status. + if solo_tags.intersects( + crate::params::SamAttributes::GXM + | crate::params::SamAttributes::GNM + | crate::params::SamAttributes::SF, + ) { + let per_align = solo.align_gene_tags(&transcripts); + crate::io::sam::add_align_gene_tags( + &mut records, + &per_align, + 1, + solo_tags, ); } for record in records { @@ -2317,6 +2413,29 @@ fn align_reads_solo( } } + // STARsolo per-read barcode tags (CR/CY/UR/UY/sM/sS/sQ, + // plus CB for a CB_samTagOut run) go on every record of + // the read, mapped or unmapped, as in STAR. + if !solo_tags.is_empty() { + let values = crate::solo::SoloTagStrings::build( + sread.barcode.as_ref(), + sread.barcode_read.as_ref(), + outcome.barcode, + cb_corrected, + ); + crate::io::sam::add_solo_barcode_tags( + &mut buffer.records, + &values.as_values(), + solo_tags, + ); + } + if track_read_info { + crate::io::sam::add_solo_read_index( + &mut buffer.records, + read_index, + ); + } + Ok(SoloReadProduct { sam_records: buffer, per_feature: outcome.per_feature, @@ -2393,6 +2512,14 @@ fn align_reads_solo_pe( let max_multimaps = params.out_filter_multimap_nmax as usize; let emit_sam = params.emits_alignments(); let output_unmapped = emit_sam && params.out_sam_unmapped != params::OutSamUnmapped::None; + let solo_tags = if emit_sam { + params.solo_sam_tags() + } else { + crate::params::SamAttributes::empty() + }; + let track_read_info = emit_sam && solo_ctx.read_info_enabled(); + let tag_out_only = params.solo_type == params::SoloType::CbSamTagOut; + let keep_barcode_read = params.solo_keeps_barcode_read(); let params_arc = Arc::new(params.clone()); struct SoloReadProduct { @@ -2600,26 +2727,50 @@ fn align_reads_solo_pe( // Solo quantification: union both mates (strand from mate 1) // for a both-mapped pair; fall back to the mapped mate for // half-mapped. - let outcome = if !both_mapped.is_empty() { + // `CB_samTagOut` corrects the barcode for the CB tag and + // counts nothing, exactly as in the single-end loop. + let (mut outcome, cb_corrected) = if tag_out_only { + let mut outcome = crate::solo::SoloReadOutcome::default(); + let corrected = pread.barcode.as_ref().map(|bc| { + let (tags, corrected) = solo.tag_barcode(bc); + outcome.barcode = Some(tags); + corrected + }); + (outcome, corrected) + } else if !both_mapped.is_empty() { let pairs: Vec<_> = both_mapped .iter() .map(|pa| (&pa.mate1_transcript, &pa.mate2_transcript)) .collect(); - solo.process_read_pe(&pairs, pread.barcode.as_ref(), &junctions) + ( + solo.process_read_pe( + &pairs, + pread.barcode.as_ref(), + &junctions, + ), + None, + ) } else if let Some(PairedAlignmentResult::HalfMapped { mapped_transcript, .. }) = results.first() { - solo.process_read( - std::slice::from_ref(mapped_transcript), - 1, - pread.barcode.as_ref(), - &junctions, + ( + solo.process_read( + std::slice::from_ref(mapped_transcript), + 1, + pread.barcode.as_ref(), + &junctions, + ), + None, ) } else { - solo.process_read(&[], 0, pread.barcode.as_ref(), &[]) + (solo.process_read(&[], 0, pread.barcode.as_ref(), &[]), None) }; + let read_index = (base + pair_idx as u64) as u32; + if track_read_info { + outcome.set_read_index(read_index); + } // SAM records (skipped under `--outSAMtype None`). if !emit_sam { @@ -2672,7 +2823,7 @@ fn align_reads_solo_pe( .collect(); // Soft-clip the fixed per-mate clips against the original // mate reads (matching SE/PE non-solo). Inert at default 10x. - let records = SamWriter::build_paired_records( + let mut records = SamWriter::build_paired_records( &out_read_name, &pread.mate1.sequence, &pread.mate1.quality, @@ -2687,11 +2838,67 @@ fn align_reads_solo_pe( params, n_for_mapq, )?; + // STARsolo gene tags: a pair is one alignment, so + // both mates carry the pair's genes. + if solo_tags + .intersects(crate::params::SamAttributes::SOLO_GENE_TAGS) + { + let pairs: Vec<_> = both_mapped + .iter() + .map(|pa| (&pa.mate1_transcript, &pa.mate2_transcript)) + .collect(); + if solo_tags.intersects( + crate::params::SamAttributes::GX + | crate::params::SamAttributes::GN, + ) { + let (gx, gn) = solo.gene_tags_pe(&pairs); + crate::io::sam::add_gene_tags( + &mut records, + &gx, + &gn, + solo_tags, + ); + } + let per_align = solo.align_gene_tags_pe(&pairs); + crate::io::sam::add_align_gene_tags( + &mut records, + &per_align, + 2, + solo_tags, + ); + } for record in records { buffer.push(record); } } + // STARsolo per-read barcode tags. sS/sQ come from the + // separate barcode read when there is one, else from + // the unclipped mate 1 (`--soloBarcodeMate 1`). + if !solo_tags.is_empty() { + let barcode_read = pread + .barcode_read + .as_ref() + .or_else(|| keep_barcode_read.then_some(&pread.mate1)); + let values = crate::solo::SoloTagStrings::build( + pread.barcode.as_ref(), + barcode_read, + outcome.barcode, + cb_corrected, + ); + crate::io::sam::add_solo_barcode_tags( + &mut buffer.records, + &values.as_values(), + solo_tags, + ); + } + if track_read_info { + crate::io::sam::add_solo_read_index( + &mut buffer.records, + read_index, + ); + } + Ok(SoloReadProduct { sam_records: buffer, per_feature: outcome.per_feature, diff --git a/src/params/mod.rs b/src/params/mod.rs index 3b507310..43f26702 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -526,8 +526,9 @@ pub struct Parameters { pub genome_transform_vcf: Option, // ── Read files ────────────────────────────────────────────────────── - /// Input read file(s); second file is mate 2 for paired-end - #[arg(long = "readFilesIn", num_args = 1..=2)] + /// Input read file(s): mate 1, then mate 2 for paired-end. A solo run adds + /// the barcode read as the last file (`cDNA_read [cDNA_read2] barcode_read`). + #[arg(long = "readFilesIn", num_args = 1..=3)] pub read_files_in: Vec, /// Command to decompress input files (e.g. "zcat" for .gz) @@ -1342,6 +1343,64 @@ impl Parameters { !matches!(self.out_std, OutStd::None) || self.out_sam_type.format != OutSamFormat::None } + /// Whether the primary alignment output is BAM (file or stdout). + /// + /// STAR emits the STARsolo barcode/gene tags (`CR CY UR UY CB UB GX GN sM sS + /// sQ`) in BAM records only, `ReadAlign_outputTranscriptSAM.cpp` drops them + /// from the SAM text path. + pub fn bam_output(&self) -> bool { + // `--outStd` replaces the file output, so it decides the format: an + // `--outStd SAM` run writes SAM text whatever `--outSAMtype` says. + match self.out_std { + OutStd::Sam => false, + OutStd::BamUnsorted | OutStd::BamSortedByCoordinate => true, + OutStd::None => self.out_sam_type.format == OutSamFormat::Bam, + } + } + + /// Whether a coordinate-sorted BAM is produced (file or stdout), STAR's + /// `outBAMcoord`. `CB`/`UB` can only be written there, since both are known + /// only after the solo counting pass. + pub fn bam_sorted_output(&self) -> bool { + match self.out_std { + OutStd::Sam | OutStd::BamUnsorted => false, + OutStd::BamSortedByCoordinate => true, + OutStd::None => { + self.out_sam_type.format == OutSamFormat::Bam + && self.out_sam_type.sort_order == Some(OutSamSortOrder::SortedByCoordinate) + } + } + } + + /// The STARsolo tags this run actually emits: those requested via + /// `--outSAMattributes` that reach a BAM record. Empty for SAM text output, + /// which STAR never decorates with them. + pub fn solo_sam_tags(&self) -> SamAttributes { + if self.bam_output() { + self.out_sam_attributes & SamAttributes::SOLO_TAGS + } else { + SamAttributes::empty() + } + } + + /// Whether the run has to track STAR's readInfo, the per-read (cell, + /// corrected UMI) pair that UMI collapsing fills and the sorted-BAM writer + /// turns into `CB`/`UB`. `CB_samTagOut` corrects the barcode inline instead, + /// and does no collapsing at all. + pub fn solo_read_info_needed(&self) -> bool { + self.solo_type != SoloType::CbSamTagOut + && self + .solo_sam_tags() + .intersects(SamAttributes::CB | SamAttributes::UB) + } + + /// Whether the whole barcode read has to be kept around (the `sS`/`sQ` tags + /// are the only consumers). + pub fn solo_keeps_barcode_read(&self) -> bool { + self.solo_sam_tags() + .intersects(SamAttributes::SS | SamAttributes::SQ) + } + /// Whether `--chimOutType` includes `Junctions` (write Chimeric.out.junction). pub fn chim_out_junctions(&self) -> bool { self.chim_out_type.iter().any(|s| s == "Junctions") @@ -1569,6 +1628,16 @@ impl Parameters { )); } + // A third read file only means anything to a solo run, where it is the + // barcode read (`ParametersSolo.cpp:124-132`). + if params.read_files_in.len() > 2 && !params.solo_enabled() { + return Err(command.error( + ErrorKind::InvalidValue, + "--readFilesIn takes at most two files (mate 1, mate 2); a third file is \ + the barcode read of a --soloType run", + )); + } + // --genomeTransformType: Haploid and Diploid are implemented. Both require // a VCF, and are incompatible with a GTF (STAR itself doesn't combine // genomeTransform with sjdb annotation at genomeGenerate). @@ -1846,16 +1915,21 @@ impl Parameters { "--soloType SmartSeq requires --readFilesManifest (a TSV of read1read2cellID per cell)", )); } - // CB_UMI_Simple needs exactly two read files: cDNA + barcode read. + // Barcode-read chemistries take `cDNA_read [cDNA_read2] barcode_read`: + // two files single-end, three when the cDNA is paired-end. + // `--soloBarcodeMate 1` is the exception (barcode on mate 1, two cDNA + // files, no barcode file), handled below. if matches!( params.solo_type, SoloType::CbUmiSimple | SoloType::CbUmiComplex | SoloType::CbSamTagOut - ) && params.read_files_in.len() != 2 + ) && params.solo_barcode_mate == 0 + && !matches!(params.read_files_in.len(), 2 | 3) { return Err(command.error( ErrorKind::InvalidValue, format!( - "--soloType {} requires exactly two --readFilesIn files (cDNA read then barcode read); got {}", + "--soloType {} requires two --readFilesIn files (cDNA read then barcode read), \ + or three for paired-end cDNA (cDNA read 1, cDNA read 2, barcode read); got {}", params.solo_type, params.read_files_in.len() ), @@ -1872,6 +1946,15 @@ impl Parameters { "--soloBarcodeMate 1 is only supported with --soloType CB_UMI_Simple", )); } + if params.read_files_in.len() != 2 { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--soloBarcodeMate 1 requires exactly two --readFilesIn cDNA mate files; got {}", + params.read_files_in.len() + ), + )); + } } other => { return Err(command.error( @@ -1895,6 +1978,67 @@ impl Parameters { )); } } + // CB / UB SAM tags (ParametersSolo.cpp:403-435). `CB_samTagOut` + // corrects the barcode as the read is processed, so it needs neither + // a sorted BAM nor a gene feature, but it has no UMI collapsing, + // hence no UB. + let cb_ub = params + .out_sam_attributes + .intersects(SamAttributes::CB | SamAttributes::UB); + if params.solo_type == SoloType::CbSamTagOut { + if params.out_sam_attributes.contains(SamAttributes::UB) { + return Err(command.error( + ErrorKind::InvalidValue, + "UB attribute (corrected UMI) in --outSAMattributes cannot be used with \ + --soloType CB_samTagOut; use UR (uncorrected UMI) instead", + )); + } + } else if cb_ub { + if !params.bam_sorted_output() { + return Err(command.error( + ErrorKind::InvalidValue, + "CB and/or UB attributes in --outSAMattributes can only be output in the \ + sorted BAM file; re-run with --outSAMtype BAM SortedByCoordinate", + )); + } + // STAR fills readInfo from the FIRST feature on the --soloFeatures + // list, which therefore has to be a gene-level one. + let first_feature = params.solo_features.first().map(String::as_str); + if !matches!( + first_feature, + Some("Gene" | "GeneFull" | "GeneFull_Ex50pAS" | "GeneFull_ExonOverIntron") + ) { + return Err(command.error( + ErrorKind::InvalidValue, + "CB and/or UB attributes in --outSAMattributes require the first \ + --soloFeatures entry to be Gene, GeneFull, GeneFull_Ex50pAS, or \ + GeneFull_ExonOverIntron", + )); + } + } + // `CB_samTagOut` only corrects the barcode against the whitelist, so + // the posterior-based multi-match modes have nothing to resolve + // (`ParametersSolo.cpp:678`), and there is no counting to spread + // multi-gene reads over (`ParametersSolo.cpp:483`). + if params.solo_type == SoloType::CbSamTagOut { + if !matches!(params.solo_cb_match_wl_type.as_str(), "Exact" | "1MM") { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "--soloCBmatchWLtype {} does not work with --soloType CB_samTagOut; \ + use Exact or 1MM", + params.solo_cb_match_wl_type + ), + )); + } + if params.solo_multi_mappers.iter().any(|m| m != "Unique") { + return Err(command.error( + ErrorKind::InvalidValue, + "multimapping options do not work for --soloType CB_samTagOut; \ + use --soloMultiMappers Unique", + )); + } + } // soloMultiMappers values. for m in ¶ms.solo_multi_mappers { if !matches!( @@ -1911,10 +2055,13 @@ impl Parameters { } // Gene-level features need a gene model (SJ does not — junctions come // from the alignments). - let needs_gtf = params - .solo_features - .iter() - .any(|f| f == "Gene" || f == "GeneFull" || f == "Velocyto"); + // `CB_samTagOut` counts nothing (`Solo.cpp:13`: no SoloFeature is + // even constructed), so it needs no gene model. + let needs_gtf = params.solo_type != SoloType::CbSamTagOut + && params + .solo_features + .iter() + .any(|f| f == "Gene" || f == "GeneFull" || f == "Velocyto"); if needs_gtf && params.sjdb_gtf_file.is_none() { return Err(command.error( ErrorKind::MissingRequiredArgument, @@ -2091,11 +2238,13 @@ impl Parameters { self.read_files_in.first() } - /// Path to the barcode (CB+UMI) read file — the SECOND `--readFilesIn` - /// file when solo is enabled. `None` if absent. + /// Path to the barcode (CB+UMI) read file: the LAST `--readFilesIn` file + /// when solo is enabled (`cDNA_read [cDNA_read2] barcode_read`). `None` if + /// absent, or when the barcode sits on mate 1 instead. pub fn barcode_read_file(&self) -> Option<&PathBuf> { - if self.solo_enabled() { - self.read_files_in.get(1) + if self.solo_enabled() && !self.solo_barcode_on_mate1() { + self.read_files_in + .get(self.read_files_in.len().checked_sub(1)?) } else { None } @@ -2107,6 +2256,13 @@ impl Parameters { self.solo_enabled() && self.solo_barcode_mate == 1 } + /// True when the solo run aligns a cDNA mate PAIR: either the barcode is on + /// mate 1 (`--soloBarcodeMate 1`, two files) or a third `--readFilesIn` file + /// carries the barcode read (STAR's `cDNA_read1 cDNA_read2 barcode_read`). + pub fn solo_paired_cdna(&self) -> bool { + self.solo_barcode_on_mate1() || (self.solo_enabled() && self.read_files_in.len() == 3) + } + /// The two cDNA mate files (mate 1, mate 2) for a `--soloBarcodeMate 1` run. pub fn solo_cdna_mate_files(&self) -> Option<(&PathBuf, &PathBuf)> { match (self.read_files_in.first(), self.read_files_in.get(1)) { @@ -2933,6 +3089,130 @@ mod tests { ); } + /// The STARsolo barcode tags parse into their own flags + /// (`Parameters_samAttributes.cpp:111-165`). + #[test] + fn solo_sam_attributes_parse() { + let p = try_parse(&[ + "--readFilesIn", + "r.fq", + "--outSAMattributes", + "NH", + "HI", + "CR", + "CY", + "UR", + "UY", + "sM", + "sS", + "sQ", + "GX", + "GN", + ]) + .unwrap(); + let a = p.out_sam_attributes; + for flag in [ + SamAttributes::CR, + SamAttributes::CY, + SamAttributes::UR, + SamAttributes::UY, + SamAttributes::SM, + SamAttributes::SS, + SamAttributes::SQ, + SamAttributes::GX, + SamAttributes::GN, + ] { + assert!(a.contains(flag), "missing {flag:?}"); + } + // Not in the preset sets. + assert!(!SamAttributes::ALL.intersects(SamAttributes::SOLO_TAGS)); + assert!(try_parse(&["--readFilesIn", "r.fq", "--outSAMattributes", "Cb"]).is_err()); + } + + /// `CB`/`UB` are filled when the sorted BAM is written, so STAR refuses them + /// with any other output type, and refuses `UB` outright for `CB_samTagOut` + /// (no UMI collapsing there). `ParametersSolo.cpp:403-435`. + #[test] + fn cb_ub_attributes_require_a_sorted_bam_and_a_gene_feature() { + let solo = [ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + "CB_UMI_Simple", + "--sjdbGTFfile", + "genes.gtf", + "--soloCBwhitelist", + "wl.txt", + ]; + let with = |extra: &[&str]| { + let mut v = solo.to_vec(); + v.extend_from_slice(extra); + try_parse(&v) + }; + + let sorted = ["--outSAMtype", "BAM", "SortedByCoordinate"]; + let mut ok = solo.to_vec(); + ok.extend_from_slice(&sorted); + ok.extend_from_slice(&["--outSAMattributes", "NH", "CB", "UB"]); + assert!(try_parse(&ok).is_ok()); + + // SAM text and unsorted BAM: refused. + assert!(with(&["--outSAMattributes", "NH", "CB"]).is_err()); + assert!( + with(&[ + "--outSAMtype", + "BAM", + "Unsorted", + "--outSAMattributes", + "NH", + "UB" + ]) + .is_err() + ); + + // `--outStd SAM` replaces the sorted BAM with SAM text, which STAR + // never tags: refused, even though --outSAMtype still says sorted BAM. + let mut std_sam = solo.to_vec(); + std_sam.extend_from_slice(&sorted); + std_sam.extend_from_slice(&["--outStd", "SAM", "--outSAMattributes", "NH", "CB"]); + assert!(try_parse(&std_sam).is_err()); + + // Sorted BAM but a non-gene first feature: refused. + let mut sj_first = solo.to_vec(); + sj_first.extend_from_slice(&sorted); + sj_first.extend_from_slice(&[ + "--soloFeatures", + "SJ", + "Gene", + "--outSAMattributes", + "NH", + "CB", + ]); + assert!(try_parse(&sj_first).is_err()); + + // CB_samTagOut corrects the barcode inline: CB needs no sorted BAM, and + // UB does not exist at all. (Its only allowed match types are Exact and + // 1MM, so the 1MM_multi default has to be overridden.) + let tag_out = [ + "--readFilesIn", + "cdna.fq", + "bc.fq", + "--soloType", + "CB_samTagOut", + "--soloCBwhitelist", + "wl.txt", + "--soloCBmatchWLtype", + "1MM", + ]; + let mut cb_only = tag_out.to_vec(); + cb_only.extend_from_slice(&["--outSAMattributes", "NH", "CB"]); + assert!(try_parse(&cb_only).is_ok()); + let mut with_ub = tag_out.to_vec(); + with_ub.extend_from_slice(&["--outSAMattributes", "NH", "UB"]); + assert!(try_parse(&with_ub).is_err()); + } + #[test] fn out_sam_order_accepts_star_values_rejects_others() { assert!(try_parse(&["--readFilesIn", "r.fq", "--outSAMorder", "Paired"]).is_ok()); diff --git a/src/params/sam.rs b/src/params/sam.rs index 0f363f61..dcff28d5 100644 --- a/src/params/sam.rs +++ b/src/params/sam.rs @@ -10,7 +10,7 @@ bitflags::bitflags! { /// Each bit corresponds to one tag the writer may emit. `STANDARD` and /// `ALL` are convenience aliases matching STAR's preset names. #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] - pub struct SamAttributes: u16 { + pub struct SamAttributes: u32 { const NH = 1 << 0; const HI = 1 << 1; const AS = 1 << 2; @@ -34,6 +34,49 @@ bitflags::bitflags! { const VA = 1 << 13; const VG = 1 << 14; + // ---- STARsolo barcode tags (BAM output only, like STAR) ---- + /// `CR:Z`: raw (uncorrected) cell barcode sequence. + const CR = 1 << 15; + /// `CY:Z`: quality string of the raw cell barcode. + const CY = 1 << 16; + /// `UR:Z`: raw (uncorrected) UMI sequence. + const UR = 1 << 17; + /// `UY:Z`: quality string of the raw UMI. + const UY = 1 << 18; + /// `CB:Z`: whitelist-corrected cell barcode. Filled at sorting time from + /// the solo read info (except `--soloType CB_samTagOut`, which corrects + /// the barcode as the read is processed). + const CB = 1 << 19; + /// `UB:Z`: collapsed (corrected) UMI. Only known after UMI collapsing, + /// so it is added when the sorted BAM is written. + const UB = 1 << 20; + /// `sM:i`: STAR's `cbMatch` code (its barcode/UMI assessment). + const SM = 1 << 21; + /// `sS:Z`: full barcode-read sequence (CB + UMI + any adapter). + const SS = 1 << 22; + /// `sQ:Z`: full barcode-read quality string. + const SQ = 1 << 23; + /// `gx:Z`: gene ids of THIS alignment, `;`-joined (multi-gene allowed, + /// unlike the read-level unique-gene `GX`). + const GXM = 1 << 24; + /// `gn:Z`: gene names of this alignment, `;`-joined. + const GNM = 1 << 25; + /// `sF:B:i`: `(overlap type, number of genes)` for the read. + const SF = 1 << 26; + + /// Every STARsolo barcode/gene tag. STAR emits these in BAM output only. + const SOLO_TAGS = + Self::CR.bits() | Self::CY.bits() | Self::UR.bits() | Self::UY.bits() + | Self::CB.bits() | Self::UB.bits() + | Self::SM.bits() | Self::SS.bits() | Self::SQ.bits() + | Self::GX.bits() | Self::GN.bits() + | Self::GXM.bits() | Self::GNM.bits() | Self::SF.bits(); + + /// The tags derived from the gene model rather than the barcode read. + const SOLO_GENE_TAGS = + Self::GX.bits() | Self::GN.bits() + | Self::GXM.bits() | Self::GNM.bits() | Self::SF.bits(); + // STAR `Standard` = NH HI AS nM (the mismatch count nM, NOT edit-distance NM). const STANDARD = Self::NH.bits() | Self::HI.bits() | Self::AS.bits() @@ -69,6 +112,18 @@ impl FromStr for SamAttributes { "RG" => Self::RG, "GX" => Self::GX, "GN" => Self::GN, + "CR" => Self::CR, + "CY" => Self::CY, + "UR" => Self::UR, + "UY" => Self::UY, + "CB" => Self::CB, + "UB" => Self::UB, + "sM" => Self::SM, + "sS" => Self::SS, + "sQ" => Self::SQ, + "gx" => Self::GXM, + "gn" => Self::GNM, + "sF" => Self::SF, "vW" => Self::VW, "vA" => Self::VA, "vG" => Self::VG, @@ -117,7 +172,8 @@ impl clap::Args for SamAttributes { .default_values(["Standard"]) .help( "SAM optional tags: Standard, All, None, or any combination of \ - NH HI AS NM nM MD jM jI XS RG vW vA vG.", + NH HI AS NM nM MD jM jI XS RG vW vA vG, plus the STARsolo tags \ + CR CY UR UY CB UB GX GN gx gn sM sS sQ sF (BAM output only).", ), ) } diff --git a/src/quant/mod.rs b/src/quant/mod.rs index efd0694f..6df6f46e 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -25,6 +25,10 @@ use crate::junction::gtf::GtfRecord; // --------------------------------------------------------------------------- /// Per-gene annotation built from GTF exon records. +/// +/// `Default` is the empty annotation, no genes, no intervals, used by runs +/// that need a solo context without a gene model (`--soloType CB_samTagOut`). +#[derive(Default)] pub struct GeneAnnotation { /// gene_id strings in GTF-file order (index = gene_idx). pub gene_ids: Vec, diff --git a/src/solo/count.rs b/src/solo/count.rs index 4ddb4143..c87054d1 100644 --- a/src/solo/count.rs +++ b/src/solo/count.rs @@ -246,33 +246,140 @@ fn connected_components(umis: &HashMap, umi_len: usize) -> u64 { } /// 1MM_Directional: a lower-count UMI within Hamming-1 of a hub whose count -/// satisfies `count_hub >= 2*count_leaf + dir_count_add` is absorbed; the -/// molecule count is the number of surviving (non-absorbed) UMIs. +/// satisfies `count_hub >= 2*count_leaf + dir_count_add` is absorbed into that +/// hub's own (already corrected) UMI; the molecule count is the number of +/// distinct surviving UMIs, STAR's `umiArrayCorrect_Directional`, which counts +/// `umiC.size()` over the corrected values. fn directional(umis: &HashMap, umi_len: usize, dir_count_add: i64) -> u64 { - // Sort by count desc, then by UMI value for determinism. - let mut items: Vec<(u64, u32)> = umis.iter().map(|(&u, &c)| (u, c)).collect(); - items.sort_by(|a, b| b.1.cmp(&a.1).then(a.0.cmp(&b.0))); - let n = items.len(); - let mut absorbed = vec![false; n]; + let corrections = directional_correction_map(umis, umi_len, dir_count_add); + umis.keys() + .map(|u| corrections.get(u).copied().unwrap_or(*u)) + .collect::>() + .len() as u64 +} + +// --------------------------------------------------------------------------- +// UMI correction maps (STAR's `umiCorrected`, for the readInfo / UB SAM tag) +// --------------------------------------------------------------------------- + +/// Swap the low and high halves of a packed UMI, as STAR's `umiSwapHalves` does +/// before its second 1MM scan. The half-swapped value is also the order STAR's +/// graph collapse walks the UMIs in, and so decides ties between equal-count +/// representatives. +fn swap_halves(umi: u64, umi_len: usize) -> u64 { + let half_bits = umi_len; // umi_len bases → 2*umi_len bits, half = umi_len bits + let mask_low = (1u64 << half_bits) - 1; + let high = umi >> half_bits; + ((umi & mask_low) << half_bits) | high +} + +/// `raw UMI -> corrected UMI` for one `(cell, gene)`, matching what STAR records +/// in `umiCorrected` for the active dedup method +/// (`SoloFeature_collapseUMIall.cpp`, `SoloFeature_collapseUMI_Graph.cpp`). +/// Only UMIs that actually change are listed; `Exact`/`NoDedup` correct nothing. +#[allow(clippy::implicit_hasher)] // always called with the default hasher +pub fn umi_correction_map( + umis: &HashMap, + method: UmiDedup, + umi_len: usize, +) -> HashMap { + let mut map = match method { + UmiDedup::Exact | UmiDedup::NoDedup => HashMap::default(), + UmiDedup::OneMmCr => cellranger_1mm_map(umis, umi_len), + UmiDedup::OneMmAll => graph_correction_map(umis, umi_len), + UmiDedup::OneMmDirectional => directional_correction_map(umis, umi_len, 0), + UmiDedup::OneMmDirectionalUmiTools => directional_correction_map(umis, umi_len, -1), + }; + map.retain(|raw, corrected| raw != corrected); + map +} + +/// 1MM_All: every UMI of a connected component is corrected to the component's +/// highest-count UMI, ties going to the one STAR meets first, the smallest +/// half-swapped value (`umiArrayCorrect_Graph`'s `umiBest` scan). +fn graph_correction_map(umis: &HashMap, umi_len: usize) -> HashMap { + let keys: Vec = umis.keys().copied().collect(); + let n = keys.len(); + let mut map = HashMap::default(); + if n <= 1 { + return map; + } + let mut parent: Vec = (0..n).collect(); + fn find(parent: &mut [usize], mut x: usize) -> usize { + while parent[x] != x { + parent[x] = parent[parent[x]]; + x = parent[x]; + } + x + } for i in 0..n { - if absorbed[i] { - continue; + for j in (i + 1)..n { + if hamming1(keys[i], keys[j], umi_len) { + let (ri, rj) = (find(&mut parent, i), find(&mut parent, j)); + if ri != rj { + parent[ri] = rj; + } + } } - let hub_count = i64::from(items[i].1); - for j in 0..n { - if i == j || absorbed[j] { - continue; + } + // Best (count, then smallest swapped value) UMI per component. + let mut best: HashMap = HashMap::default(); + for (i, &umi) in keys.iter().enumerate() { + let root = find(&mut parent, i); + let count = umis[&umi]; + let swapped = swap_halves(umi, umi_len); + match best.get(&root) { + Some(&(bc, bs, _)) if bc > count || (bc == count && bs <= swapped) => {} + _ => { + best.insert(root, (count, swapped, umi)); } - let leaf_count = i64::from(items[j].1); - if leaf_count <= hub_count - && hub_count >= 2 * leaf_count + dir_count_add - && hamming1(items[i].0, items[j].0, umi_len) + } + } + // A component of one is never coloured by STAR, so it is never recorded. + let mut comp_size: HashMap = HashMap::default(); + for i in 0..n { + let root = find(&mut parent, i); + *comp_size.entry(root).or_insert(0) += 1; + } + for (i, &umi) in keys.iter().enumerate() { + let root = find(&mut parent, i); + if comp_size[&root] > 1 + && let Some(&(_, _, rep)) = best.get(&root) + { + map.insert(umi, rep); + } + } + map +} + +/// 1MM_Directional: each UMI, scanned in descending count order, is corrected to +/// the *corrected* value of the first earlier (higher-count) UMI within one +/// mismatch whose count satisfies `hub >= 2*leaf + dir_count_add`, STAR's +/// `umiArrayCorrect_Directional`, chain and all. +fn directional_correction_map( + umis: &HashMap, + umi_len: usize, + dir_count_add: i64, +) -> HashMap { + let mut items: Vec<(u64, u32)> = umis.iter().map(|(&u, &c)| (u, c)).collect(); + items.sort_by(|a, b| b.1.cmp(&a.1).then(a.0.cmp(&b.0))); + let mut corrected: Vec = items.iter().map(|&(u, _)| u).collect(); + for iu in 1..items.len() { + for iuu in 0..iu { + if i64::from(items[iuu].1) >= 2 * i64::from(items[iu].1) + dir_count_add + && hamming1(items[iu].0, items[iuu].0, umi_len) { - absorbed[j] = true; + corrected[iu] = corrected[iuu]; + break; } } } - (n - absorbed.iter().filter(|&&a| a).count()) as u64 + items + .iter() + .map(|&(u, _)| u) + .zip(corrected) + .filter(|(raw, corr)| raw != corr) + .collect() } // --------------------------------------------------------------------------- @@ -363,6 +470,7 @@ fn build_matrix_body( pseudocount: f64, dir: &Path, n_features: usize, + record_read_info: bool, ) -> Result<(tempfile::NamedTempFile, MatrixStats), Error> { let mut body_tmp = tempfile::Builder::new() .prefix(".matrix_body") @@ -385,6 +493,7 @@ fn build_matrix_body( cb, umi: m.umi, gene: m.gene, + read_index: m.read_index, }); } } @@ -414,6 +523,9 @@ fn build_matrix_body( body: Vec, stat: Option, genes: Vec, + /// `(read index, readInfo entry)` for this cell's reads, when the + /// `CB`/`UB` SAM tags asked for STAR's readInfo. + read_info: Vec<(u32, crate::solo::ReadInfo)>, } let cell_outs: Vec = bounds .par_iter() @@ -457,6 +569,35 @@ fn build_matrix_body( }; cell_entries.sort_unstable_by_key(|&(g, _)| g); + // STAR's readInfo: every read counted for this cell records the + // cell it landed in and the UMI it collapsed onto + // (`SoloFeature_collapseUMIall.cpp:250-265`). + let mut read_info: Vec<(u32, crate::solo::ReadInfo)> = Vec::new(); + if record_read_info { + let mut per_gene: HashMap> = HashMap::default(); + for (&umi, genes) in &umi_genes { + for (&gene, &rc) in genes { + *per_gene.entry(gene).or_default().entry(umi).or_insert(0) += rc; + } + } + let corrections: HashMap> = per_gene + .iter() + .map(|(&gene, umis)| (gene, umi_correction_map(umis, method, umi_len))) + .collect(); + read_info.reserve(j - i); + for r in &records[i..j] { + if r.read_index == crate::solo::NO_READ_INDEX { + continue; + } + let umi = corrections + .get(&r.gene) + .and_then(|m| m.get(&r.umi)) + .copied() + .unwrap_or(r.umi); + read_info.push((r.read_index, crate::solo::ReadInfo { cb, umi })); + } + } + let n_reads = (j - i) as u64; let n_genes = cell_entries.len() as u32; let mut n_umis = 0u64; @@ -477,11 +618,15 @@ fn build_matrix_body( body: cbody, stat, genes, + read_info, } }) .collect(); // Sequential merge: byte order preserved (CB-ascending, gene-ascending). + let mut info_guard = record_read_info + .then(|| ctx.read_info.as_ref().map(|m| m.lock().unwrap())) + .flatten(); for co in cell_outs { body.write_all(&co.body).map_err(|e| Error::io(e, dir))?; nnz += co.genes.len(); @@ -491,7 +636,15 @@ fn build_matrix_body( if let Some(s) = co.stat { cell_stats.push(s); } + if let Some(info) = info_guard.as_mut() { + for (read_index, entry) in co.read_info { + if let Some(slot) = info.get_mut(read_index as usize) { + *slot = entry; + } + } + } } + drop(info_guard); body.flush().map_err(|e| Error::io(e, dir))?; } @@ -1369,7 +1522,7 @@ pub fn write_gene_matrix( let multi_methods = MultiMethod::parse_list(¶ms.solo_multi_mappers); // One {prefix}{soloOutFileNames[0]}/{raw,filtered}/ per feature. - for (feature, recorder) in ctx.features.iter().zip(&ctx.recorders) { + for (fi, (feature, recorder)) in ctx.features.iter().zip(&ctx.recorders).enumerate() { let feature_dir = params.output_path(&format!("{solo_dir}{}/", feature.dir_name())); let raw_dir = feature_dir.join("raw"); std::fs::create_dir_all(&raw_dir).map_err(|e| Error::io(e, &raw_dir))?; @@ -1385,6 +1538,9 @@ pub fn write_gene_matrix( pseudocount, &raw_dir, n_genes, + // STAR fills readInfo from one feature only: the first on the + // --soloFeatures list (`ParametersSolo.cpp:423-434`). + fi == 0 && ctx.read_info_enabled(), )?; write_features( &raw_dir.join(&features_name), @@ -2473,6 +2629,48 @@ mod tests { assert_eq!(dedup_count(&c, UmiDedup::OneMmDirectional, 4), 2); } + /// The correction map behind the `UB` SAM tag: 1MM_All sends every UMI of a + /// connected component to the component's highest-count member, and leaves + /// UMIs that collapse with nothing untouched. + #[test] + fn graph_correction_maps_each_component_to_its_top_umi() { + // AAAA(5)–AAAC(1)–AACC(2) is one component (AAAC bridges); TTTT alone. + let c = counts(&[("AAAA", 5), ("AAAC", 1), ("AACC", 2), ("TTTT", 3)]); + let map = umi_correction_map(&c, UmiDedup::OneMmAll, 4); + assert_eq!(map.get(&umi("AAAC")), Some(&umi("AAAA"))); + assert_eq!(map.get(&umi("AACC")), Some(&umi("AAAA"))); + // The representative and the isolated UMI are not corrections. + assert_eq!(map.get(&umi("AAAA")), None); + assert_eq!(map.get(&umi("TTTT")), None); + // The molecule count agrees with the map: 2 components. + assert_eq!(dedup_count(&c, UmiDedup::OneMmAll, 4), 2); + } + + /// Directional follows the chain: a leaf takes the *corrected* UMI of the + /// hub that absorbed it, so a two-step chain lands on the head. + #[test] + fn directional_correction_follows_the_chain() { + // AAAA(9) ← AAAC(4) ← AAAG(1): AAAC absorbs into AAAA (9 >= 2*4), + // AAAG into AAAC's corrected value (4 >= 2*1) → all three are one. + let c = counts(&[("AAAA", 9), ("AAAC", 4), ("AAAG", 1)]); + let map = umi_correction_map(&c, UmiDedup::OneMmDirectional, 4); + assert_eq!(map.get(&umi("AAAC")), Some(&umi("AAAA"))); + assert_eq!(map.get(&umi("AAAG")), Some(&umi("AAAA"))); + assert_eq!(dedup_count(&c, UmiDedup::OneMmDirectional, 4), 1); + } + + /// Exact and NoDedup correct nothing, so `UB` is the read's own UMI. + #[test] + fn exact_dedup_corrects_nothing() { + let c = counts(&[("AAAA", 5), ("AAAC", 1)]); + assert!(umi_correction_map(&c, UmiDedup::Exact, 4).is_empty()); + assert!(umi_correction_map(&c, UmiDedup::NoDedup, 4).is_empty()); + // 1MM_CR does correct, and only lists UMIs that actually change. + let cr = umi_correction_map(&c, UmiDedup::OneMmCr, 4); + assert_eq!(cr.get(&umi("AAAC")), Some(&umi("AAAA"))); + assert_eq!(cr.get(&umi("AAAA")), None); + } + #[test] fn cellranger_1mm_collapses_neighbor() { // AAAA (5) and AAAC (1) are 1MM → low-count corrected to high-count → diff --git a/src/solo/gene.rs b/src/solo/gene.rs index 5d4e4a60..f89cd93d 100644 --- a/src/solo/gene.rs +++ b/src/solo/gene.rs @@ -304,6 +304,131 @@ pub fn classify_read( }) } +/// STAR's `ReadAnnotFeature::overlapTypes` (`ReadAnnotations.h:13`), the first +/// value of the `sF` SAM tag. The read takes the lowest-numbered type any of its +/// alignments reaches. +/// +/// `ExonicSense50p` / `Exonic50pAntisense` are STAR's `GeneFull_Ex50pAS` types +/// and never arise here, since that feature is not implemented. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[repr(i32)] +pub enum OverlapType { + None = 0, + ExonicSense = 1, + ExonicAntisense = 2, + ExonicSense50p = 3, + Exonic50pAntisense = 4, + IntronicSense = 5, + IntronicAntisense = 6, + Intergenic = 7, +} + +impl OverlapType { + /// The sense-strand types, for which STAR reports `sF:B:i,-1,-1` on an + /// alignment that carries no gene of its own + /// (`ReadAlign_alignBAM.cpp:441-447`). + pub fn is_sense(self) -> bool { + matches!( + self, + Self::ExonicSense | Self::ExonicSense50p | Self::IntronicSense + ) + } +} + +/// The genes each individual alignment of a read belongs to, for the `gx`/`gn` +/// SAM tags (STAR's `readAnnot.annotFeatures[f].fAlign`), plus the read-level +/// gene set size and overlap type behind `sF`. +#[derive(Debug, Clone, Default)] +pub struct AlignGenes { + /// One sorted gene-index list per alignment, parallel to the transcripts. + pub per_align: Vec>, + /// Union over alignments (STAR's `fSet`). + pub gene_set: Vec, + /// STAR's `ovType`. + pub ov_type: i32, +} + +/// Per-alignment gene assignment for `gx`/`gn`/`sF`, on the `Gene` (exon, +/// transcript-concordant) or `GeneFull` (gene body) basis. +/// +/// Mirrors the per-transcript half of [`classify_read`]: the same concordance +/// and strand rules, kept per alignment instead of unioned over the read. +pub fn align_genes( + transcripts: &[Transcript], + gene_ann: &GeneAnnotation, + strand: SoloStrand, + feature: SoloFeature, +) -> AlignGenes { + let want_exon = feature == SoloFeature::Gene; + let mut out = AlignGenes { + per_align: Vec::with_capacity(transcripts.len()), + ..Default::default() + }; + if transcripts.is_empty() { + out.ov_type = OverlapType::None as i32; + return out; + } + let (mut exon_sense, mut exon_anti) = (false, false); + let (mut body_sense, mut body_anti) = (false, false); + let mut raw: Vec = Vec::new(); + + for tr in transcripts { + let mut genes: Vec = Vec::new(); + // Exon overlap decides the exonic/intronic half of `ovType` even for the + // GeneFull basis, which is what makes an intronic GeneFull read report + // `intronic` rather than `exonic`. + gene_ann.overlapping_genes_into(tr, &mut raw); + for &g in &raw { + let concordant = tr + .exons + .iter() + .all(|b| gene_ann.block_is_exonic(g, b.genome_start, b.genome_end)); + if !concordant { + continue; + } + if strand_keeps(strand, gene_ann.gene_is_reverse[g], tr.is_reverse) { + exon_sense = true; + if want_exon { + genes.push(g as u32); + } + } else { + exon_anti = true; + } + } + gene_ann.overlapping_genes_full_into(tr, &mut raw); + for &g in &raw { + if strand_keeps(strand, gene_ann.gene_is_reverse[g], tr.is_reverse) { + body_sense = true; + if !want_exon { + genes.push(g as u32); + } + } else { + body_anti = true; + } + } + genes.sort_unstable(); + genes.dedup(); + out.gene_set.extend_from_slice(&genes); + out.per_align.push(genes); + } + out.gene_set.sort_unstable(); + out.gene_set.dedup(); + + // Lower types win, as in STAR's `otFinal` scan. + out.ov_type = if exon_sense { + OverlapType::ExonicSense + } else if exon_anti { + OverlapType::ExonicAntisense + } else if body_sense { + OverlapType::IntronicSense + } else if body_anti { + OverlapType::IntronicAntisense + } else { + OverlapType::Intergenic + } as i32; + out +} + /// Assign a single-end (cDNA) read to a gene from its alignment set, using the /// `Gene` (exonic) or `GeneFull` (gene-body, intron-inclusive) overlap basis. /// Thin wrapper over [`classify_read`] for the single-feature case (and tests). @@ -490,6 +615,63 @@ mod tests { } } + /// `gx`/`gn` are per alignment, so two alignments of one read report their + /// own genes; `sF` reports the read-level overlap type and gene count, and + /// falls back to `(-1, -1)` on an alignment with no gene of its own. + #[test] + fn align_genes_are_per_alignment_with_a_read_level_overlap_type() { + // Ga (+) exons [100,200); Gb (+) exons [400,500). Intergenic elsewhere. + let g = genome(); + let exons = vec![gtf_exon(101, 200, '+', "Ga"), gtf_exon(401, 500, '+', "Gb")]; + let ann = GeneAnnotation::from_gtf_exons(&exons, &g); + + let hits = vec![read_at(110, 150, false), read_at(410, 450, false)]; + let out = align_genes(&hits, &ann, SoloStrand::Forward, SoloFeature::Gene); + assert_eq!(out.per_align, vec![vec![0], vec![1]]); + assert_eq!(out.gene_set, vec![0, 1]); + assert_eq!(out.ov_type, OverlapType::ExonicSense as i32); + + // Second alignment intergenic: it keeps no gene, but the read is still + // exonic-sense from the first. + let mixed = vec![read_at(110, 150, false), read_at(900, 950, false)]; + let out = align_genes(&mixed, &ann, SoloStrand::Forward, SoloFeature::Gene); + assert_eq!(out.per_align, vec![vec![0], Vec::::new()]); + assert_eq!(out.ov_type, OverlapType::ExonicSense as i32); + assert!(OverlapType::ExonicSense.is_sense()); + + // Antisense read over the same exon: exonicAS, no gene kept. + let anti = vec![read_at(110, 150, true)]; + let out = align_genes(&anti, &ann, SoloStrand::Forward, SoloFeature::Gene); + assert_eq!(out.per_align, vec![Vec::::new()]); + assert_eq!(out.ov_type, OverlapType::ExonicAntisense as i32); + assert!(!OverlapType::ExonicAntisense.is_sense()); + + // No gene anywhere near: intergenic. + let far = vec![read_at(900, 950, false)]; + let out = align_genes(&far, &ann, SoloStrand::Forward, SoloFeature::Gene); + assert_eq!(out.ov_type, OverlapType::Intergenic as i32); + assert!(out.gene_set.is_empty()); + } + + /// A read inside an intron is intronic-sense, and carries the gene only on + /// the `GeneFull` basis. + #[test] + fn intronic_reads_report_the_intronic_overlap_type() { + // One gene, two exons: body [100,500), intron [200,400). + let g = genome(); + let exons = vec![gtf_exon(101, 200, '+', "Ga"), gtf_exon(401, 500, '+', "Ga")]; + let ann = GeneAnnotation::from_gtf_exons(&exons, &g); + let intronic = vec![read_at(250, 300, false)]; + + let full = align_genes(&intronic, &ann, SoloStrand::Forward, SoloFeature::GeneFull); + assert_eq!(full.per_align, vec![vec![0]]); + assert_eq!(full.ov_type, OverlapType::IntronicSense as i32); + + let gene = align_genes(&intronic, &ann, SoloStrand::Forward, SoloFeature::Gene); + assert_eq!(gene.per_align, vec![Vec::::new()]); + assert_eq!(gene.ov_type, OverlapType::IntronicSense as i32); + } + #[test] fn classify_read_regions_and_antisense() { // Ga (+): exons [100,200) and [400,500) → body [100,500), intron [200,400). diff --git a/src/solo/mod.rs b/src/solo/mod.rs index 22d63e02..a13daf8d 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -217,6 +217,126 @@ fn decode_seq(encoded: &[u8]) -> String { encoded.iter().map(|&b| decode_base(b) as char).collect() } +/// Decoded per-read values behind the STARsolo barcode SAM tags, built once per +/// read and borrowed by [`crate::io::sam::add_solo_barcode_tags`]. +#[derive(Debug, Default)] +pub struct SoloTagStrings { + cb_seq: String, + cb_qual: String, + umi_seq: String, + umi_qual: String, + barcode_seq: String, + barcode_qual: String, + cb_corrected: Option, + cb_match: i32, +} + +impl SoloTagStrings { + /// Decode one read's barcode into tag values. + /// + /// `barcode` is `None` when the barcode read was too short to carry a + /// CB+UMI: STAR pads such a read with `N`s, which always scores `cbMatch=-2` + /// (Ns in the barcode), so that is what `sM` reports, but there is no + /// barcode sequence to put in `CR`/`CY`/`UR`/`UY`, and those are left off. + pub fn build( + barcode: Option<&CellBarcode>, + barcode_read: Option<&EncodedRead>, + tags: Option, + cb_corrected: Option, + ) -> Self { + let mut out = Self { + cb_match: tags.map_or(-2, |t| t.cb_match), + cb_corrected, + ..Default::default() + }; + if let Some(bc) = barcode { + out.cb_seq = bc.cb_string(); + out.cb_qual = String::from_utf8_lossy(&bc.cb_qual).into_owned(); + out.umi_seq = bc.umi_string(); + out.umi_qual = String::from_utf8_lossy(&bc.umi_qual).into_owned(); + } + if let Some(read) = barcode_read { + out.barcode_seq = decode_seq(&read.sequence); + out.barcode_qual = String::from_utf8_lossy(&read.quality).into_owned(); + } + out + } + + /// Borrowed view for the SAM writer. + pub fn as_values(&self) -> crate::io::sam::SoloBarcodeTagValues<'_> { + crate::io::sam::SoloBarcodeTagValues { + cb_seq: &self.cb_seq, + cb_qual: &self.cb_qual, + umi_seq: &self.umi_seq, + umi_qual: &self.umi_qual, + cb_corrected: self.cb_corrected.as_deref(), + cb_match: self.cb_match, + barcode_seq: &self.barcode_seq, + barcode_qual: &self.barcode_qual, + } + } +} + +/// STAR's `--soloBarcodeReadLength` handling for a separate barcode read +/// (`SoloReadBarcode_getCBandUMI.cpp:228-241`). +/// +/// The default (`1`) means the barcode read must be exactly CB+UMI long, and any +/// other length is a fatal input error. `0` turns the check off, and a read +/// shorter than CB+UMI is then padded with `N` (quality `H`) so it still parses +/// (and, having Ns, is never counted). +#[derive(Debug, Clone, Copy, Default)] +pub struct BarcodeReadLength { + expected: Option, + cbumi_len: usize, +} + +impl BarcodeReadLength { + pub fn from_params(params: &Parameters) -> Self { + // The barcode inside a cDNA mate has the mate's length, so STAR turns the + // check off there (`ParametersSolo.cpp:143`), as it does for the + // variable-geometry Complex chemistry. + if params.solo_barcode_on_mate1() || params.solo_type == SoloType::CbUmiComplex { + return Self::default(); + } + let cbumi_len = params.solo_cb_len as usize + params.solo_umi_len as usize; + let expected = match params.solo_barcode_read_length { + 0 => None, + 1 => Some(cbumi_len), + n if n > 0 => Some(n as usize), + _ => None, + }; + Self { + expected, + cbumi_len, + } + } + + /// Check the barcode read's length, padding it up to CB+UMI when checking is + /// off and it falls short. + pub fn apply(&self, read: &mut EncodedRead) -> Result<(), Error> { + if let Some(expected) = self.expected { + if read.sequence.len() != expected { + return Err(Error::from(std::io::Error::new( + std::io::ErrorKind::InvalidData, + format!( + "solo: barcode read '{}' is {} bases, not the expected {}. \ + If the CB+UMI length is not the barcode read length, set \ + --soloBarcodeReadLength to that length, or 0 to skip the check.", + read.name, + read.sequence.len(), + expected + ), + ))); + } + } else if read.sequence.len() < self.cbumi_len { + // 4 = N in the genome encoding; 'H' is STAR's filler quality. + read.sequence.resize(self.cbumi_len, 4); + read.quality.resize(self.cbumi_len, b'H'); + } + Ok(()) + } +} + /// Reads cDNA reads and their paired barcode reads in lockstep from two FASTQ /// files. The cDNA read flows into the normal alignment path; the barcode read /// is parsed into a [`CellBarcode`] (or `None` when too short). @@ -224,6 +344,8 @@ pub struct SoloReadReader { cdna: FastqReader, barcode: FastqReader, layout: SoloBarcodeLayout, + keep_barcode_read: bool, + barcode_len: BarcodeReadLength, } /// One cDNA read paired with its (optional) extracted barcode. @@ -231,6 +353,8 @@ pub struct SoloRead { pub cdna: EncodedRead, /// `None` when the barcode read was too short to extract CB+UMI. pub barcode: Option, + /// The whole barcode read, kept only for the `sS`/`sQ` SAM tags. + pub barcode_read: Option, } impl SoloReadReader { @@ -240,11 +364,15 @@ impl SoloReadReader { barcode_path: &Path, layout: SoloBarcodeLayout, decompress_cmd: Option<&str>, + keep_barcode_read: bool, + barcode_len: BarcodeReadLength, ) -> Result { Ok(Self { cdna: FastqReader::open(cdna_path, decompress_cmd)?, barcode: FastqReader::open(barcode_path, decompress_cmd)?, layout, + keep_barcode_read, + barcode_len, }) } @@ -254,9 +382,15 @@ impl SoloReadReader { let cdna_opt = self.cdna.next_encoded()?; let barcode_opt = self.barcode.next_encoded()?; match (cdna_opt, barcode_opt) { - (Some(cdna), Some(bc)) => { + (Some(cdna), Some(mut bc)) => { + self.barcode_len.apply(&mut bc)?; let barcode = self.layout.extract(&bc); - Ok(Some(SoloRead { cdna, barcode })) + let barcode_read = self.keep_barcode_read.then_some(bc); + Ok(Some(SoloRead { + cdna, + barcode, + barcode_read, + })) } (None, None) => Ok(None), (Some(_), None) => Err(Error::from(std::io::Error::new( @@ -289,7 +423,7 @@ impl SoloReadReader { pub fn open_reader(params: &Parameters) -> Result { debug_assert!(matches!( params.solo_type, - SoloType::CbUmiSimple | SoloType::CbUmiComplex + SoloType::CbUmiSimple | SoloType::CbUmiComplex | SoloType::CbSamTagOut )); let cdna = params.cdna_read_file().ok_or_else(|| { Error::from(std::io::Error::new( @@ -304,7 +438,14 @@ pub fn open_reader(params: &Parameters) -> Result { )) })?; let layout = SoloBarcodeLayout::from_params(params); - SoloReadReader::open(cdna, barcode, layout, params.read_files_command.as_deref()) + SoloReadReader::open( + cdna, + barcode, + layout, + params.read_files_command.as_deref(), + params.solo_keeps_barcode_read(), + BarcodeReadLength::from_params(params), + ) } /// One paired-end solo read for `--soloBarcodeMate 1` (5' 10x): both mates carry @@ -312,41 +453,73 @@ pub fn open_reader(params: &Parameters) -> Result { pub struct SoloPairedRead { pub mate1: EncodedRead, pub mate2: EncodedRead, - /// `None` when mate 1 was too short to extract CB+UMI. + /// `None` when the barcode read (or mate 1) was too short to extract CB+UMI. pub barcode: Option, + /// The whole barcode read, kept only for the `sS`/`sQ` SAM tags. `None` for + /// a `--soloBarcodeMate 1` run, where mate 1 itself is the barcode read. + pub barcode_read: Option, } -/// Reads the two cDNA mate files in lockstep for a `--soloBarcodeMate 1` run, -/// extracting the barcode from the start of mate 1. +/// Reads a solo cDNA mate pair in lockstep, taking the barcode either from the +/// start of mate 1 (`--soloBarcodeMate 1`) or from a third barcode-read file +/// (`--readFilesIn cDNA_read1 cDNA_read2 barcode_read`). pub struct SoloPairedReader { mate1: FastqReader, mate2: FastqReader, + /// The separate barcode read, when the run has one. + barcode: Option, layout: SoloBarcodeLayout, + keep_barcode_read: bool, + barcode_len: BarcodeReadLength, } impl SoloPairedReader { pub fn open( mate1_path: &Path, mate2_path: &Path, + barcode_path: Option<&Path>, layout: SoloBarcodeLayout, decompress_cmd: Option<&str>, + keep_barcode_read: bool, + barcode_len: BarcodeReadLength, ) -> Result { Ok(Self { mate1: FastqReader::open(mate1_path, decompress_cmd)?, mate2: FastqReader::open(mate2_path, decompress_cmd)?, + barcode: barcode_path + .map(|p| FastqReader::open(p, decompress_cmd)) + .transpose()?, layout, + keep_barcode_read, + barcode_len, }) } - /// Fetch the next (mate1, mate2) pair with the barcode extracted from mate 1. + /// Fetch the next (mate1, mate2) pair with its barcode. pub fn next_read(&mut self) -> Result, Error> { match (self.mate1.next_encoded()?, self.mate2.next_encoded()?) { (Some(mate1), Some(mate2)) => { - let barcode = self.layout.extract(&mate1); + let (barcode, barcode_read) = match &mut self.barcode { + // Separate barcode read: must stay in lockstep with the mates. + Some(reader) => { + let mut bc = reader.next_encoded()?.ok_or_else(|| { + Error::from(std::io::Error::new( + std::io::ErrorKind::UnexpectedEof, + "solo: barcode read file has fewer reads than the cDNA mate files", + )) + })?; + self.barcode_len.apply(&mut bc)?; + let extracted = self.layout.extract(&bc); + (extracted, self.keep_barcode_read.then_some(bc)) + } + // `--soloBarcodeMate 1`: the barcode is a prefix of mate 1. + None => (self.layout.extract(&mate1), None), + }; Ok(Some(SoloPairedRead { mate1, mate2, barcode, + barcode_read, })) } (None, None) => Ok(None), @@ -376,11 +549,19 @@ pub fn open_paired_reader(params: &Parameters) -> Result, pub umi: u64, pub gene: u32, + /// Input read index, as in [`SoloCountRecord`]. + pub read_index: u32, } /// A read that mapped to multiple genes (gene-ambiguous). Distributed across its @@ -580,6 +770,29 @@ pub struct SoloContext { /// since building the transcriptome costs a GTF pass. pub transcriptome: Option, pub transcript3p: Option>, + /// STAR's `readInfo` (`ParametersSolo.cpp:418-435`): the cell and corrected + /// UMI each read ended up counted under, filled by UMI collapsing and read + /// back when the sorted BAM is written to fill `CB`/`UB`. `None` unless one + /// of those tags was requested. + pub read_info: Option>>, +} + +/// One `readInfo` entry: what a read was counted as, once collapsing has run. +/// Both fields keep STAR's "undefined" sentinel (all ones), which surfaces as +/// `CB:Z:-` / `UB:Z:-`. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct ReadInfo { + pub cb: u32, + pub umi: u64, +} + +impl Default for ReadInfo { + fn default() -> Self { + Self { + cb: u32::MAX, + umi: u64::MAX, + } + } } /// Per-region read tallies for the `Summary.csv` mapping funnel (uniquely-mapped @@ -602,6 +815,44 @@ pub struct SoloReadOutcome { pub sj: Vec, /// Velocyto record for this read (resolved CB, gene-assigned), if enabled. pub velocyto: Option, + /// Barcode facts for the per-read SAM tags. `None` when the barcode read was + /// too short to extract a CB+UMI at all. + pub barcode: Option, +} + +impl SoloReadOutcome { + /// Stamp the input read index onto every count record this read produced, so + /// UMI collapsing can fill STAR's readInfo (the `CB`/`UB` SAM tags). + pub fn set_read_index(&mut self, read_index: u32) { + for fo in &mut self.per_feature { + if let Some(r) = &mut fo.record { + r.read_index = read_index; + } + if let Some(m) = &mut fo.multi { + m.read_index = read_index; + } + } + } +} + +/// The `gx`/`gn`/`sF` tag values of one alignment of a read. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct AlignGeneTag { + pub gx: String, + pub gn: String, + pub sf: [i32; 2], +} + +/// What the barcode of one read resolved to, for the STARsolo SAM tags. +/// +/// `cb_match` is STAR's `cbMatch` code (the `sM` tag); `cb_index` is the +/// whitelist entry the barcode corrected to, which also seeds the `CB` tag of a +/// `CB_samTagOut` run and the readInfo entry that fills `CB`/`UB` in the sorted +/// BAM. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct SoloBarcodeTags { + pub cb_match: i32, + pub cb_index: Option, } /// The record(s) one read produces for a single feature. @@ -651,27 +902,35 @@ impl SoloContext { }; // Gene model from the GTF (validated to be present for Gene/GeneFull). - let gtf_path = params.sjdb_gtf_file.as_ref().ok_or_else(|| { - Error::from(std::io::Error::new( - std::io::ErrorKind::InvalidInput, - "STARsolo Gene feature requires --sjdbGTFfile", - )) - })?; - let exons = crate::junction::gtf::parse_gtf_configured( - gtf_path, - ¶ms.sjdb_gtf_feature_exon, - ¶ms.sjdb_gtf_chr_prefix, - )?; - let gene_ann = GeneAnnotation::from_gtf_exons_configured( - &exons, - genome, - ¶ms.sjdb_gtf_tag_exon_parent_gene, - ); - log::info!( - "STARsolo: {} genes loaded from {}", - gene_ann.n_genes(), - gtf_path.display() - ); + // `CB_samTagOut` quantifies nothing (`Solo.cpp:13` builds no + // SoloFeature), so it runs without a gene model. + let tag_out_only = params.solo_type == SoloType::CbSamTagOut; + let (exons, gene_ann) = if tag_out_only { + (Vec::new(), GeneAnnotation::default()) + } else { + let gtf_path = params.sjdb_gtf_file.as_ref().ok_or_else(|| { + Error::from(std::io::Error::new( + std::io::ErrorKind::InvalidInput, + "STARsolo Gene feature requires --sjdbGTFfile", + )) + })?; + let exons = crate::junction::gtf::parse_gtf_configured( + gtf_path, + ¶ms.sjdb_gtf_feature_exon, + ¶ms.sjdb_gtf_chr_prefix, + )?; + let gene_ann = GeneAnnotation::from_gtf_exons_configured( + &exons, + genome, + ¶ms.sjdb_gtf_tag_exon_parent_gene, + ); + log::info!( + "STARsolo: {} genes loaded from {}", + gene_ann.n_genes(), + gtf_path.display() + ); + (exons, gene_ann) + }; let strand: SoloStrand = params.solo_strand.parse().map_err(|e: String| { Error::from(std::io::Error::new(std::io::ErrorKind::InvalidInput, e)) @@ -679,15 +938,19 @@ impl SoloContext { // Quantified gene features (Gene, GeneFull). Validation guarantees these // parse; default to Gene if somehow empty. - let features: Vec = params - .solo_features - .iter() - .filter_map(|f| f.parse().ok()) - .collect(); - let features = if features.is_empty() { - vec![SoloFeature::Gene] + let features: Vec = if tag_out_only { + Vec::new() } else { - features + let parsed: Vec = params + .solo_features + .iter() + .filter_map(|f| f.parse().ok()) + .collect(); + if parsed.is_empty() { + vec![SoloFeature::Gene] + } else { + parsed + } }; let recorders = features.iter().map(|_| SoloRecorder::new()).collect(); let feature_reads = features.iter().map(|_| AtomicU64::new(0)).collect(); @@ -734,9 +997,55 @@ impl SoloContext { .transpose()?, transcript3p: transcript3p .then(|| Mutex::new(crate::solo::transcript3p::Transcript3pAcc::new())), + // Sized once the read count is known (`reserve_read_info`). + read_info: params + .solo_read_info_needed() + .then(|| Mutex::new(Vec::new())), }) } + /// Whether this run has to track STAR's readInfo (the `CB`/`UB` SAM tags). + pub fn read_info_enabled(&self) -> bool { + self.read_info.is_some() + } + + /// `--soloType CB_samTagOut`: match the barcode to the whitelist and stop + /// there. Returns the `sM` facts plus the corrected barcode for the `CB` + /// tag, the whitelist entry for an Exact/1MM hit, the raw barcode when + /// there is no whitelist to correct against, and `"-"` otherwise + /// (`SoloReadBarcode_getCBandUMI.cpp:311-328`). + pub fn tag_barcode(&self, bc: &CellBarcode) -> (SoloBarcodeTags, String) { + let cb_match = self + .whitelist + .match_cb(&bc.cb_seq, &bc.cb_qual, self.match_type); + self.stats.record_cb(&cb_match); + let corrected = match cb_match.resolved_index() { + Some(idx) => self + .whitelist + .barcode_string(idx) + // No whitelist: STAR passes the barcode through uncorrected. + .unwrap_or_else(|| bc.cb_string()), + None => "-".to_string(), + }; + ( + SoloBarcodeTags { + cb_match: cb_match.star_code(), + cb_index: cb_match.resolved_index(), + }, + corrected, + ) + } + + /// Size the readInfo array for `n_reads` input reads, once the alignment + /// pass has counted them. No-op unless `CB`/`UB` were requested. + pub fn reserve_read_info(&self, n_reads: usize) { + if let Some(info) = &self.read_info { + let mut info = info.lock().unwrap(); + info.clear(); + info.resize(n_reads, ReadInfo::default()); + } + } + /// Process one solo read: match the cell barcode, validate the UMI, assign /// a gene, and (on success) produce a count record. Stats are recorded /// here; the returned records are appended to the recorder by the caller. @@ -744,7 +1053,7 @@ impl SoloContext { /// `(gene_id, gene_name)` when uniquely assigned, else `("-", "-")` /// (STARsolo convention). Drives `--outSAMattributes GX GN`. pub fn gene_tags<'a>(&'a self, transcripts: &[Transcript]) -> (&'a str, &'a str) { - match assign_gene_se(transcripts, &self.gene_ann, self.strand, SoloFeature::Gene) { + match assign_gene_se(transcripts, &self.gene_ann, self.strand, self.tag_feature()) { GeneAssignment::Gene(g) => ( self.gene_ann.gene_ids[g as usize].as_str(), self.gene_ann.gene_names[g as usize].as_str(), @@ -753,6 +1062,106 @@ impl SoloContext { } } + /// `GX`/`GN` for a paired-end solo read: the pair's single gene, or `("-", + /// "-")` when it has none or several. Owned strings, since the pair's + /// effective transcripts are built here. + pub fn gene_tags_pe(&self, pairs: &[(&Transcript, &Transcript)]) -> (String, String) { + let mut eff: Vec = Vec::with_capacity(pairs.len() * 2); + for (m1, m2) in pairs { + let mut m2c = (*m2).clone(); + m2c.is_reverse = m1.is_reverse; + eff.push((*m1).clone()); + eff.push(m2c); + } + let (gx, gn) = self.gene_tags(&eff); + (gx.to_string(), gn.to_string()) + } + + /// The feature the gene SAM tags are computed from: STAR's `samAttrFeature`, + /// which is the first entry of `--soloFeatures` + /// (`ParametersSolo.cpp:423`). + fn tag_feature(&self) -> SoloFeature { + self.features.first().copied().unwrap_or(SoloFeature::Gene) + } + + /// Per-alignment `gx`/`gn`/`sF` values for one read. + /// + /// `gx`/`gn` list every gene of *that* alignment (`;`-joined, `"-"` when it + /// has none), unlike `GX`/`GN`, which name the read's single gene. `sF` + /// carries `(overlap type, genes for the read)`, or `(-1, -1)` when the read + /// overlaps a sense-strand feature but this particular alignment has no gene + /// (`ReadAlign_alignBAM.cpp:441-473`). + pub fn align_gene_tags(&self, transcripts: &[Transcript]) -> Vec { + let genes = crate::solo::gene::align_genes( + transcripts, + &self.gene_ann, + self.strand, + self.tag_feature(), + ); + self.gene_tag_values(&genes, 1) + } + + /// The same for a paired-end solo read: an alignment is a mate pair, so the + /// two mates' genes are merged and both records carry the pair's tags. + pub fn align_gene_tags_pe(&self, pairs: &[(&Transcript, &Transcript)]) -> Vec { + // Both mates evaluated against the pair's (mate 1's) strand, as in + // `process_read_pe`. + let mut eff: Vec = Vec::with_capacity(pairs.len() * 2); + for (m1, m2) in pairs { + let mut m2c = (*m2).clone(); + m2c.is_reverse = m1.is_reverse; + eff.push((*m1).clone()); + eff.push(m2c); + } + let genes = + crate::solo::gene::align_genes(&eff, &self.gene_ann, self.strand, self.tag_feature()); + self.gene_tag_values(&genes, 2) + } + + /// Render `AlignGenes` into per-alignment tag strings, merging every + /// `mates_per_align` consecutive entries into one alignment. + fn gene_tag_values( + &self, + genes: &crate::solo::gene::AlignGenes, + mates_per_align: usize, + ) -> Vec { + let sense = matches!( + genes.ov_type, + x if x == crate::solo::gene::OverlapType::ExonicSense as i32 + || x == crate::solo::gene::OverlapType::ExonicSense50p as i32 + || x == crate::solo::gene::OverlapType::IntronicSense as i32 + ); + let n_genes = genes.gene_set.len() as i32; + genes + .per_align + .chunks(mates_per_align.max(1)) + .map(|mates| { + let mut gs: Vec = mates.concat(); + gs.sort_unstable(); + gs.dedup(); + let join = |names: &[String]| { + if gs.is_empty() { + "-".to_string() + } else { + gs.iter() + .map(|&g| names[g as usize].as_str()) + .collect::>() + .join(";") + } + }; + AlignGeneTag { + gx: join(&self.gene_ann.gene_ids), + gn: join(&self.gene_ann.gene_names), + sf: if sense && gs.is_empty() { + [-1, -1] + } else { + [genes.ov_type, n_genes] + }, + } + }) + .collect() + } + pub fn process_read( &self, cdna_transcripts: &[Transcript], @@ -810,6 +1219,17 @@ impl SoloContext { .match_cb(&bc.cb_seq, &bc.cb_qual, self.match_type); self.stats.record_cb(&cb_match); + // Per-read SAM-tag facts, recorded before any early return so that a + // rejected barcode still tags its alignments. STAR reports a UMI + // rejection in place of the CB code (`getCBandUMI.cpp:304`). + let umi_status = check_umi(&bc.umi_seq); + out.barcode = Some(SoloBarcodeTags { + cb_match: umi_status + .star_code() + .unwrap_or_else(|| cb_match.star_code()), + cb_index: cb_match.resolved_index(), + }); + let cb_resolved: Option = match &cb_match { CbMatch::Exact(idx) | CbMatch::Corrected(idx) => Some(*idx), CbMatch::Multi(_) => None, // deferred to collation @@ -817,7 +1237,7 @@ impl SoloContext { }; // UMI validity. - let umi = match check_umi(&bc.umi_seq) { + let umi = match umi_status { UmiCheck::Ok(packed) => { self.stats.record_umi(&UmiCheck::Ok(packed)); packed @@ -897,12 +1317,20 @@ impl SoloContext { // valid-barcode reads (STARsolo "Reads Mapped to "). self.feature_reads[fi].fetch_add(1, Ordering::Relaxed); match (cb_resolved, &cb_match) { - (Some(cb), _) => fo.record = Some(SoloCountRecord { cb, umi, gene }), + (Some(cb), _) => { + fo.record = Some(SoloCountRecord { + cb, + umi, + gene, + read_index: NO_READ_INDEX, + }); + } (None, CbMatch::Multi(cands)) => { fo.multi = Some(SoloMultiRecord { candidates: cands.clone(), umi, gene, + read_index: NO_READ_INDEX, }); } (None, _) => unreachable!("non-multi unresolved CB returned early"), @@ -1154,6 +1582,48 @@ mod tests { assert!(bc.umi_has_n()); } + /// STAR checks the barcode read's length by default and refuses anything + /// else; with the check off, a short read is padded with `N` (which then + /// scores as an N-containing barcode rather than being silently dropped). + #[test] + fn barcode_read_length_is_checked_then_padded() { + use crate::io::fastq::encode_base; + + let read = |seq: &str| EncodedRead { + name: "r1".to_string(), + sequence: seq.bytes().map(encode_base).collect(), + quality: vec![b'I'; seq.len()], + }; + // Default: exactly CB+UMI = 26 bases. + let checked = BarcodeReadLength { + expected: Some(26), + cbumi_len: 26, + }; + assert!(checked.apply(&mut read("A".repeat(26).as_str())).is_ok()); + let err = checked + .apply(&mut read("A".repeat(20).as_str())) + .unwrap_err() + .to_string(); + assert!(err.contains("20 bases, not the expected 26"), "{err}"); + // Longer is refused too. + assert!(checked.apply(&mut read("A".repeat(30).as_str())).is_err()); + + // --soloBarcodeReadLength 0: pad up to CB+UMI with N (encoded 4) / 'H'. + let unchecked = BarcodeReadLength { + expected: None, + cbumi_len: 26, + }; + let mut short = read("ACGT"); + unchecked.apply(&mut short).unwrap(); + assert_eq!(short.sequence.len(), 26); + assert!(short.sequence[4..].iter().all(|&b| b == 4)); + assert!(short.quality[4..].iter().all(|&q| q == b'H')); + // A read at or past the length is left alone. + let mut long = read("A".repeat(30).as_str()); + unchecked.apply(&mut long).unwrap(); + assert_eq!(long.sequence.len(), 30); + } + #[test] fn reader_pairs_cdna_and_barcode() { use std::io::Write; @@ -1177,7 +1647,15 @@ mod tests { .unwrap(); bc.flush().unwrap(); - let mut reader = SoloReadReader::open(cdna.path(), bc.path(), v2_layout(), None).unwrap(); + let mut reader = SoloReadReader::open( + cdna.path(), + bc.path(), + v2_layout(), + None, + false, + BarcodeReadLength::default(), + ) + .unwrap(); let batch = reader.read_batch(10).unwrap(); assert_eq!(batch.len(), 2); assert_eq!(batch[0].cdna.name, "r1"); @@ -1209,7 +1687,15 @@ mod tests { .unwrap(); bc.flush().unwrap(); - let mut reader = SoloReadReader::open(cdna.path(), bc.path(), v2_layout(), None).unwrap(); + let mut reader = SoloReadReader::open( + cdna.path(), + bc.path(), + v2_layout(), + None, + false, + BarcodeReadLength::default(), + ) + .unwrap(); assert!(reader.read_batch(10).is_err()); } } diff --git a/src/solo/whitelist.rs b/src/solo/whitelist.rs index 1d882ae8..adab4e77 100644 --- a/src/solo/whitelist.rs +++ b/src/solo/whitelist.rs @@ -175,6 +175,42 @@ pub enum CbMatch { MultMatchRejected, } +impl CbMatch { + /// STAR's `cbMatch` code for this outcome, as reported in the `sM` SAM tag + /// (`SoloReadBarcode_getCBandUMI.cpp:9-90`). Multi-match carries the number + /// of whitelist candidates. + pub fn star_code(&self) -> i32 { + match self { + Self::Exact(_) => 0, + Self::Corrected(_) => 1, + Self::Multi(cands) => cands.len() as i32, + Self::NoMatch => -1, + Self::NinCb => -2, + Self::MultMatchRejected => -3, + } + } + + /// Whitelist index when the barcode resolved to a single cell. + pub fn resolved_index(&self) -> Option { + match self { + Self::Exact(i) | Self::Corrected(i) => Some(*i), + _ => None, + } + } +} + +impl UmiCheck { + /// STAR's `umiCheck` code, which overwrites `cbMatch` (and so the `sM` tag) + /// when the UMI is rejected. A valid UMI leaves the CB code in place. + pub fn star_code(&self) -> Option { + match self { + Self::Ok(_) => None, + Self::NinUmi => Some(-23), + Self::Homopolymer => Some(-24), + } + } +} + // --------------------------------------------------------------------------- // UMI validity (matches STAR umiCheck=-23 / -24) // --------------------------------------------------------------------------- From 4479fd82e5e7daf48a6e1649f1d6f5008e513e65 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 11 Aug 2026 22:20:32 +0200 Subject: [PATCH 2/3] test(solo): cover the barcode SAM tags, CB_samTagOut and 3-file solo input Three integration tests on the synthetic genome: - `test_starsolo_barcode_sam_tags`: every tag on a counting run, including the collapsed `UB` (both UMI clouds resolve to the higher-count UMI) and the absence of the internal read-index tag from the output. - `test_solo_cb_sam_tag_out`: exact, 1MM-corrected and unmatched barcodes give `CB` / `sM` of `(cb, 0)`, `(cb, 1)` and `("-", -1)`, with no `UB` and no `Solo.out`. - `test_solo_paired_cdna_with_separate_barcode_read`: paired-end cDNA with a third barcode-read file tags both mates of every pair. Unit tests cover the per-alignment gene sets and overlap types, the UMI correction maps behind `UB`, the barcode-read length check and padding, and the new `--outSAMattributes` tokens and their validation. Co-Authored-By: Claude Opus 5 (1M context) --- tests/alignment_features.rs | 371 ++++++++++++++++++++++++++++++++++++ 1 file changed, 371 insertions(+) diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index 86691d31..97f075a1 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2204,3 +2204,374 @@ fn test_read_name_separator_cuts_the_qname_and_is_configurable() { ); } } + +// --------------------------------------------------------------------------- +// Test, STARsolo (Phase 14.7): CB/UB/CR/CY/UR/UY/GX/GN/sM/sS/sQ SAM tags +// --------------------------------------------------------------------------- + +/// Read every record's optional tags out of a BAM as `(tag, string value)`. +/// Integer tags are rendered decimal, so one helper covers `sM` too. +fn bam_tag_maps(path: &Path) -> Vec> { + use noodles::sam::alignment::record::data::field::Value; + let mut reader = bam::io::Reader::new(fs::File::open(path).unwrap()); + let _header = reader.read_header().expect("BAM header readable"); + let mut out = Vec::new(); + for rec in reader.records() { + let rec = rec.expect("valid BAM record"); + let mut map = std::collections::HashMap::new(); + for field in rec.data().iter() { + let (tag, value) = field.expect("valid aux field"); + let key = String::from_utf8(tag.as_ref().to_vec()).unwrap(); + let rendered = match value { + Value::String(s) => String::from_utf8_lossy(s.as_ref()).into_owned(), + Value::Int8(v) => v.to_string(), + Value::UInt8(v) => v.to_string(), + Value::Int16(v) => v.to_string(), + Value::UInt16(v) => v.to_string(), + Value::Int32(v) => v.to_string(), + Value::UInt32(v) => v.to_string(), + // `B:i` arrays (sF) render as their comma-joined values. + Value::Array(array) => { + use noodles::sam::alignment::record::data::field::value::Array as A; + match array { + A::Int32(values) => values + .iter() + .map(|v| v.expect("valid array element").to_string()) + .collect::>() + .join(","), + other => format!("{other:?}"), + } + } + other => format!("{other:?}"), + }; + map.insert(key, rendered); + } + out.push(map); + } + out +} + +/// A counting solo run with the barcode tags requested writes them on every +/// alignment: the raw barcode/UMI and their qualities, the whole barcode read, +/// the match code, the gene, and, from the post-counting readInfo, the +/// corrected cell barcode and collapsed UMI. +#[test] +fn test_starsolo_barcode_sam_tags() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let gtf = write_gtf(&tmpdir); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", Some(>f)); + + let cdna_path = tmpdir.path().join("cdna.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + + let cb = "AAAACCCCGGGGTTTT"; + // Two UMIs one mismatch apart, 3 reads vs 1 read: 1MM_All collapses them + // into a single molecule, whose UMI is the higher-count one. + let umi_hi = "ACGTACGTAC"; + let umi_lo = "ACGTACGTAG"; + let n_reads = 4usize; + { + let mut cf = fs::File::create(&cdna_path).unwrap(); + let mut bf = fs::File::create(&barcode_path).unwrap(); + let exon1 = &genome[10000..10050]; + for i in 0..n_reads { + writeln!(cf, "@read{i}").unwrap(); + cf.write_all(exon1).unwrap(); + writeln!(cf, "\n+\n{}", "I".repeat(50)).unwrap(); + + let umi = if i < 3 { umi_hi } else { umi_lo }; + writeln!(bf, "@read{i}").unwrap(); + writeln!(bf, "{cb}{umi}").unwrap(); + writeln!(bf, "+\n{}", "I".repeat(26)).unwrap(); + } + } + fs::write(&wl_path, format!("{cb}\nCCCCGGGGTTTTAAAA\n")).unwrap(); + + let output_dir = tmpdir.path().join("out_tags"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_UMI_Simple", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloFeatures", + "Gene", + "--sjdbGTFfile", + gtf.to_str().unwrap(), + "--outSAMtype", + "BAM", + "SortedByCoordinate", + "--outSAMattributes", + "NH", + "HI", + "AS", + "nM", + "CR", + "CY", + "UR", + "UY", + "CB", + "UB", + "GX", + "GN", + "gx", + "gn", + "sM", + "sS", + "sQ", + "sF", + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let bam_path = output_dir.join("Aligned.sortedByCoord.out.bam"); + let tags = bam_tag_maps(&bam_path); + assert_eq!(tags.len(), n_reads, "expected one record per read"); + + let mut collapsed = 0usize; + for t in &tags { + assert_eq!(t.get("CR").map(String::as_str), Some(cb)); + assert_eq!( + t.get("CY").map(String::as_str), + Some("I".repeat(16).as_str()) + ); + assert_eq!(t.get("CB").map(String::as_str), Some(cb), "corrected CB"); + assert_eq!(t.get("GX").map(String::as_str), Some("G1")); + assert_eq!(t.get("GN").map(String::as_str), Some("G1")); + // Per-alignment gene lists: one gene here, so the same value as GX/GN. + assert_eq!(t.get("gx").map(String::as_str), Some("G1")); + assert_eq!(t.get("gn").map(String::as_str), Some("G1")); + // sF = (overlap type, genes for the read) = (exonic sense, 1). + assert_eq!(t.get("sF").map(String::as_str), Some("1,1")); + // Exact whitelist match → STAR's cbMatch code 0. + assert_eq!(t.get("sM").map(String::as_str), Some("0")); + // sS/sQ carry the whole barcode read. + assert_eq!(t.get("sS").unwrap().len(), 26); + assert_eq!( + t.get("sQ").map(String::as_str), + Some("I".repeat(26).as_str()) + ); + + let ur = t.get("UR").expect("UR tag"); + assert!(ur == umi_hi || ur == umi_lo, "unexpected raw UMI {ur}"); + assert_eq!( + t.get("UY").map(String::as_str), + Some("I".repeat(10).as_str()) + ); + // Every read collapses onto the 3-read UMI. + assert_eq!( + t.get("UB").map(String::as_str), + Some(umi_hi), + "collapsed UB" + ); + if ur == umi_lo { + collapsed += 1; + } + // The private read-index tag never reaches the output. + assert!(!t.contains_key("zR"), "internal zR tag leaked into the BAM"); + } + assert_eq!( + collapsed, 1, + "expected the single low-count UMI to be corrected" + ); +} + +/// `--soloType CB_samTagOut` corrects the barcode as the read is processed and +/// counts nothing: CB comes out without a sorted BAM or a gene model, and no +/// `Solo.out` directory is written. +#[test] +fn test_solo_cb_sam_tag_out() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + let cdna_path = tmpdir.path().join("cdna.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + + let cb = "AAAACCCCGGGGTTTT"; + // One exact barcode, one a single mismatch away (corrected under 1MM), one + // unrelated (no match → CB:Z:-). + let observed = [cb, "AAAACCCCGGGGTTTA", "TTTTGGGGCCCCAAAA"]; + { + let mut cf = fs::File::create(&cdna_path).unwrap(); + let mut bf = fs::File::create(&barcode_path).unwrap(); + let exon1 = &genome[10000..10050]; + for (i, bc) in observed.iter().enumerate() { + writeln!(cf, "@read{i}").unwrap(); + cf.write_all(exon1).unwrap(); + writeln!(cf, "\n+\n{}", "I".repeat(50)).unwrap(); + writeln!(bf, "@read{i}").unwrap(); + writeln!(bf, "{bc}ACGTACGTAC").unwrap(); + writeln!(bf, "+\n{}", "I".repeat(26)).unwrap(); + } + } + fs::write(&wl_path, format!("{cb}\nCCCCGGGGTTTTAAAA\n")).unwrap(); + + let output_dir = tmpdir.path().join("out_tagout"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + cdna_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_samTagOut", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloCBmatchWLtype", + "1MM", + "--outSAMtype", + "BAM", + "Unsorted", + "--outSAMattributes", + "NH", + "HI", + "CR", + "UR", + "CB", + "sM", + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let tags = bam_tag_maps(&output_dir.join("Aligned.out.bam")); + assert_eq!(tags.len(), observed.len()); + let by_cr: std::collections::HashMap<&str, &std::collections::HashMap> = + tags.iter().map(|t| (t["CR"].as_str(), t)).collect(); + + // Exact match: cbMatch 0, corrected to itself. + assert_eq!(by_cr[cb]["CB"], cb); + assert_eq!(by_cr[cb]["sM"], "0"); + // One mismatch: corrected to the whitelist barcode, cbMatch 1. + assert_eq!(by_cr["AAAACCCCGGGGTTTA"]["CB"], cb); + assert_eq!(by_cr["AAAACCCCGGGGTTTA"]["sM"], "1"); + // No match within one edit: "-" and cbMatch -1. + assert_eq!(by_cr["TTTTGGGGCCCCAAAA"]["CB"], "-"); + assert_eq!(by_cr["TTTTGGGGCCCCAAAA"]["sM"], "-1"); + // No UB (there is no UMI collapsing at all) and no count matrices. + assert!(tags.iter().all(|t| !t.contains_key("UB"))); + assert!( + !output_dir.join("Solo.out").exists(), + "CB_samTagOut must not write Solo.out" + ); +} + +/// Paired-end cDNA with a separate barcode read +/// (`--readFilesIn cDNA_read1 cDNA_read2 barcode_read`): both mates align as a +/// pair and both carry the barcode tags. Run under `CB_samTagOut`, which STAR +/// documents for exactly this three-file layout. +#[test] +fn test_solo_paired_cdna_with_separate_barcode_read() { + let tmpdir = TempDir::new().unwrap(); + let genome = build_genome(); + let fasta = write_fasta(&tmpdir, &genome); + let genome_dir = tmpdir.path().join("genome"); + build_index(&fasta, &genome_dir, "7", None); + + let mate1_path = tmpdir.path().join("mate1.fq"); + let mate2_path = tmpdir.path().join("mate2.fq"); + let barcode_path = tmpdir.path().join("barcode.fq"); + let wl_path = tmpdir.path().join("whitelist.txt"); + let cb = "AAAACCCCGGGGTTTT"; + let n_pairs = 4usize; + { + let mut f1 = fs::File::create(&mate1_path).unwrap(); + let mut f2 = fs::File::create(&mate2_path).unwrap(); + let mut fb = fs::File::create(&barcode_path).unwrap(); + for i in 0..n_pairs { + // FR pair: mate 1 forward at p, mate 2 the reverse complement of the + // fragment's right end. + let p = 500 + i * 200; + let seq1 = &genome[p..p + 50]; + let seq2 = rc(&genome[p + 150..p + 200]); + + writeln!(f1, "@pair{i}").unwrap(); + f1.write_all(seq1).unwrap(); + writeln!(f1, "\n+\n{}", "I".repeat(50)).unwrap(); + writeln!(f2, "@pair{i}").unwrap(); + f2.write_all(&seq2).unwrap(); + writeln!(f2, "\n+\n{}", "I".repeat(50)).unwrap(); + writeln!(fb, "@pair{i}").unwrap(); + writeln!(fb, "{cb}ACGTACGTAC").unwrap(); + writeln!(fb, "+\n{}", "I".repeat(26)).unwrap(); + } + } + fs::write(&wl_path, format!("{cb}\nCCCCGGGGTTTTAAAA\n")).unwrap(); + + let output_dir = tmpdir.path().join("out_pe_solo"); + fs::create_dir_all(&output_dir).unwrap(); + let prefix = format!("{}/", output_dir.display()); + + cargo_bin_cmd!("rustar-aligner") + .args([ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + mate1_path.to_str().unwrap(), + mate2_path.to_str().unwrap(), + barcode_path.to_str().unwrap(), + "--soloType", + "CB_samTagOut", + "--soloCBwhitelist", + wl_path.to_str().unwrap(), + "--soloCBmatchWLtype", + "1MM", + "--outSAMtype", + "BAM", + "Unsorted", + "--outSAMattributes", + "NH", + "HI", + "CR", + "UR", + "CB", + "sM", + "sS", + "--outFileNamePrefix", + &prefix, + ]) + .assert() + .success(); + + let tags = bam_tag_maps(&output_dir.join("Aligned.out.bam")); + // Two records (one per mate) for every pair. + assert_eq!(tags.len(), n_pairs * 2, "expected both mates per pair"); + for t in &tags { + assert_eq!(t.get("CR").map(String::as_str), Some(cb)); + assert_eq!(t.get("CB").map(String::as_str), Some(cb)); + assert_eq!(t.get("UR").map(String::as_str), Some("ACGTACGTAC")); + assert_eq!(t.get("sM").map(String::as_str), Some("0")); + // sS is the whole barcode read, not either cDNA mate. + assert_eq!(t.get("sS").map(|s| s.len()), Some(26)); + } + assert!(!output_dir.join("Solo.out").exists()); +} From c9e649c02f2737e7a98850ce68b8b25bb1fe1fdc Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 11 Aug 2026 22:20:39 +0200 Subject: [PATCH 3/3] docs(solo): record Phase 14.7 and close out Phase 14 Marks sub-phase 14.7 complete and Phase 14 (STARsolo) with it, describes the tag implementation in ROADMAP.md and CLAUDE.md, adds the CHANGELOG entries, and documents the new `--outSAMattributes` tokens plus the three-file solo `--readFilesIn` layout on the docs site. Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 40 +++++++++++++++++++ CLAUDE.md | 2 +- ROADMAP.md | 10 +++-- docs-old/phase14_starsolo.md | 2 +- .../content/docs/reference/cli-parameters.md | 4 +- 5 files changed, 51 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 14b4346b..781ab0c8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,46 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Features +- **STARsolo per-read barcode SAM tags (Phase 14.7).** `--outSAMattributes` + now accepts `CR CY UR UY CB UB gx gn sM sS sQ sF` next to the existing + `GX GN`, and emits them in BAM output only, as STAR does. + + - `CR`/`CY`/`UR`/`UY` (raw barcode and UMI with their qualities), `sM` + (STAR's `cbMatch` assessment code) and `sS`/`sQ` (the whole barcode read) + are written as each read is processed, on mapped and unmapped records + alike, for both the single-end and the `--soloBarcodeMate 1` solo path. + - `CB` (corrected barcode) and `UB` (collapsed UMI) come from STAR's + readInfo: UMI collapsing records what each read was counted as, and the + buffered records are rewritten before the sort. Reads that were not + counted get `-`, as in STAR. Both tags therefore require + `--outSAMtype BAM SortedByCoordinate` and a gene-level first + `--soloFeatures` entry, which is now validated. + - `--soloType CB_samTagOut` is implemented: whitelist correction into `CB` + with no gene model, no UMI collapsing and no `Solo.out` output. As in + STAR, it rejects `UB` and accepts only `Exact`/`1MM` for + `--soloCBmatchWLtype`. + - `gx`/`gn` name every gene of the alignment they sit on (`;`-joined), and + `sF` reports `(overlap type, genes for the read)`, falling back to + `(-1, -1)` on a sense-strand read whose alignment has no gene of its own. + All the gene tags now reach the paired-end solo path too, where a mate pair + counts as one alignment. + - `--soloUMIdedup 1MM_Directional`/`1MM_Directional_UMItools` now count + distinct corrected UMIs off STAR's absorb chain + (`umiArrayCorrect_Directional`) rather than counting unabsorbed UMIs; the + two agree except where a chain is longer than one step. + +- **Paired-end cDNA with a separate barcode read.** `--readFilesIn` accepts + STAR's three-file solo layout (`cDNA_read1 cDNA_read2 barcode_read`), so + paired-end cDNA now works with a separate barcode read for every barcode + chemistry, `CB_samTagOut` included. A third file without `--soloType` is + refused rather than ignored. + +- **`--soloBarcodeReadLength` is honoured.** As in STAR, the default expects the + barcode read to be exactly CB+UMI long and treats any other length as a fatal + input error naming the read; `0` turns the check off and pads a short read + with `N` (quality `H`), which then scores as an N-containing barcode instead + of being dropped silently. + - **CLI and output parity: SAM/SJ/read-input knobs and the STAR limit surface** — 30 further STAR 2.7.11b parameters. (`--outSAMorder` came from #145.) diff --git a/CLAUDE.md b/CLAUDE.md index 149b5a3d..47a8089c 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -32,7 +32,7 @@ Always run `cargo clippy --all-targets`, `cargo fmt --check`, and `cargo test` b ## Current Status -**396 tests passing, 0 clippy warnings.** SE: 8613/8926 compare_sam.py (96.5%; note: lower due to seeded-RNG tie-break PR diverging from STAR's mt19937), **99.815% faithfulness (tie-adjusted)** (8611/8627 non-tie reads exact), 299 tie-breaking diffs excluded. 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **0 STAR-only / 0 rustar-aligner-only SE reads**. PE: **8390 both-mapped** (STAR: 8390), **0 half-mapped**, 0 MAPQ inflations / 0 deflations, **99.883% PE exact faithfulness (tie-adjusted)** (16284/16306, 475 tie-breaking diffs excluded), **0 proper-pair diffs**, **0 NH diffs**. Phase 17.A: `scoreSeedBest` pre-extension. Phase 17.B: per-mate seeding. Phase 17.C: STAR-faithful SCORE-GATE + mappedFilter. Phase 17.D: combined-span penalty fix + dedup ordering. Phase 17.8: `--quantMode GeneCounts`. Phase E fix (2026-04-21): mate_id-aware diagonal dedup. Phase E2 (2026-04-22): STAR-faithful combined-read seeding. Phase E3 (2026-04-22): combined-threshold half-mapped fallback. Phase E4 (2026-04-22): PE-CHECK2 unconditional. Phase E5 (2026-04-23): split_combined_wt n_mismatch propagation. Phase E6 (2026-04-24): tie-adjusted faithfulness metric in assess_faithfulness.py. Phase F1: --runRNGseed + seeded primary tie-break (PR #5). Phase F2: --outSAMattrRGline (PR #6). Phase F3: --quantMode TranscriptomeSAM (PR #7). Phase F4: SJDB insertion into Genome+SA at genomeGenerate (PR #8). Phase G1 (2026-04-29): junction_shifts fix in split_combined_wt (rDNA cross-copy false-splice filter). Phase G2 (2026-04-29): MAX_RECURSION 10k→100k + sa_pos_to_forward overflow fix (ERR12389696.7118031 NH=3→9). Phase 17.2 (2026-04-29): coordinate-sorted BAM output (`--outSAMtype BAM SortedByCoordinate` → `Aligned.sortedByCoord.out.bam`). Phase 17.4 (2026-04-29): `--outReadsUnmapped Fastx` → `Unmapped.out.mate1` / `Unmapped.out.mate2`; writes unmapped + TooManyLoci reads; PE writes both mates for fully-unmapped and half-mapped pairs. Phase 17.6 (2026-05-01): `--outStd SAM/BAM_Unsorted/BAM_SortedByCoordinate` — routes primary alignment output to stdout via `Box` trait dispatch; `SamStdoutWriter`, `BamStdoutWriter`, `SortedBamStdoutWriter` in sam.rs/bam.rs; verified with samtools pipe (967 records). Phase G3 (2026-05-01): SA tie-breaking fix — `compare_suffixes` tie-breaker changed from `pos_b.cmp(&pos_a)` to `packed_a.cmp(&packed_b)` (ascending by packed SA value with strand bit); rustar-aligner SA is now **byte-for-byte identical** to STAR's SA for the yeast genome (10,862 → 0 entry diffs). diff AS: 6→4 cases (4 remaining are rustar-aligner improvements: .844151 VIII 0mm vs STAR VII 6mm, .4972950 spliced vs unspliced mate2). Phase 17.3 (2026-05-01): PE chimeric detection — `detect_inter_mate_chimeric` in `chimeric/detect.rs`; intra-mate multi-cluster chimeric via cluster splitting + mate2 read_pos adjustment; inter-mate chimeric for discordant pairs (diff chr, same strand, or >1Mb); `align_paired_read` returns 4-tuple including `Vec`; no benchmark regression (8390 both-mapped, 0 half-mapped). Phase 17.11 (2026-05-01): `--chimOutType WithinBAM` — chimeric alignments written as supplementary records (FLAG 0x800) in primary BAM; donor record has full SEQ + SA tag; acceptor has FLAG 0x800 + SA tag + empty SEQ; `build_within_bam_records` in `chimeric/output.rs`; `chim_out_junctions()` / `chim_out_within_bam()` helpers in params.rs; supports mixed `--chimOutType Junctions WithinBAM`. Phase 17.7 (2026-05-01): GTF tag parameters — `--sjdbGTFchrPrefix`, `--sjdbGTFfeatureExon`, `--sjdbGTFtagExonParentTranscript`, `--sjdbGTFtagExonParentGene`; `_configured` variants in `junction/gtf.rs`, `quant/mod.rs`, `quant/transcriptome.rs`, `junction/mod.rs`; all 4 production paths thread params; backward-compat wrappers preserve zero test disruption. Phase 17.9 (2026-05-01): `--outBAMcompression` (BGZF level -1–9, default 1; -1/0=NONE, 1-8=flate2 levels, ≥9=BEST) + `--limitBAMsortRAM` (bytes, 0=unlimited; aborts sort if ~400 bytes/record estimate exceeds limit); `bgzf_compression()` + `make_bgzf_writer()` helpers in `io/bam.rs`; threaded through all 4 BAM writers (unsorted file, sorted file, unsorted stdout, sorted stdout). PE chimericDetectionOld (2026-05-01): per-mate `detect_chimeric_old` called on `all_m1_transcripts` / `all_m2_transcripts` pools after `filter_paired_transcripts` in `read_align.rs`. Phase 17.12 (2026-05-01): BySJout disk buffering — `BySJReadMeta` struct + `NamedTempFile` SAM temp file replaces `Vec`; `create_bysj_writer` / `bysj_write_records` / `bysj_read_n_records` helpers in `io/sam.rs`; `tempfile` moved to `[dependencies]`. Phase 17.13 (2026-05-01): 8 integration tests in `tests/alignment_features.rs` — synthetic 20kb genome with planted GT-AG intron; tests cover BAM output, PE alignment, spliced reads, BySJout, GeneCounts, unmapped output, two-pass mode. Phase 12.2 (2026-05-04): SE chimeric Tier 1b soft-clip re-mapping — `detect_from_soft_clips` in `chimeric/detect.rs` re-seeds the primary alignment's soft-clipped bases when `detect_chimeric_old` finds no partner; `adjust_read_positions` helper shifts sub-seq coords into full-read space for right clips; called as Step 3c in `read_align.rs`. Phase 17.10 (2026-05-04): Chimeric Tier 3 — `detect_from_chimeric_residuals` in `chimeric/detect.rs` re-seeds outer uncovered read regions (before donor / after acceptor) of each found chimeric pair; enables 3-way gene-fusion detection; called as Step 3d in `read_align.rs`. See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs-old/](docs-old/) for per-phase development notes. The published Astro Starlight docs site is in [docs/](docs/). +**642 tests passing (599 lib + 43 integration), 0 clippy warnings.** SE: 8613/8926 compare_sam.py (96.5%; note: lower due to seeded-RNG tie-break PR diverging from STAR's mt19937), **99.815% faithfulness (tie-adjusted)** (8611/8627 non-tie reads exact), 299 tie-breaking diffs excluded. 1 CIGAR-only disagree (ERR12389696.13573895, insertion placement, seed-level tie). **0 STAR-only / 0 rustar-aligner-only SE reads**. PE: **8390 both-mapped** (STAR: 8390), **0 half-mapped**, 0 MAPQ inflations / 0 deflations, **99.883% PE exact faithfulness (tie-adjusted)** (16284/16306, 475 tie-breaking diffs excluded), **0 proper-pair diffs**, **0 NH diffs**. Phase 17.A: `scoreSeedBest` pre-extension. Phase 17.B: per-mate seeding. Phase 17.C: STAR-faithful SCORE-GATE + mappedFilter. Phase 17.D: combined-span penalty fix + dedup ordering. Phase 17.8: `--quantMode GeneCounts`. Phase E fix (2026-04-21): mate_id-aware diagonal dedup. Phase E2 (2026-04-22): STAR-faithful combined-read seeding. Phase E3 (2026-04-22): combined-threshold half-mapped fallback. Phase E4 (2026-04-22): PE-CHECK2 unconditional. Phase E5 (2026-04-23): split_combined_wt n_mismatch propagation. Phase E6 (2026-04-24): tie-adjusted faithfulness metric in assess_faithfulness.py. Phase F1: --runRNGseed + seeded primary tie-break (PR #5). Phase F2: --outSAMattrRGline (PR #6). Phase F3: --quantMode TranscriptomeSAM (PR #7). Phase F4: SJDB insertion into Genome+SA at genomeGenerate (PR #8). Phase G1 (2026-04-29): junction_shifts fix in split_combined_wt (rDNA cross-copy false-splice filter). Phase G2 (2026-04-29): MAX_RECURSION 10k→100k + sa_pos_to_forward overflow fix (ERR12389696.7118031 NH=3→9). Phase 17.2 (2026-04-29): coordinate-sorted BAM output (`--outSAMtype BAM SortedByCoordinate` → `Aligned.sortedByCoord.out.bam`). Phase 17.4 (2026-04-29): `--outReadsUnmapped Fastx` → `Unmapped.out.mate1` / `Unmapped.out.mate2`; writes unmapped + TooManyLoci reads; PE writes both mates for fully-unmapped and half-mapped pairs. Phase 17.6 (2026-05-01): `--outStd SAM/BAM_Unsorted/BAM_SortedByCoordinate` — routes primary alignment output to stdout via `Box` trait dispatch; `SamStdoutWriter`, `BamStdoutWriter`, `SortedBamStdoutWriter` in sam.rs/bam.rs; verified with samtools pipe (967 records). Phase G3 (2026-05-01): SA tie-breaking fix — `compare_suffixes` tie-breaker changed from `pos_b.cmp(&pos_a)` to `packed_a.cmp(&packed_b)` (ascending by packed SA value with strand bit); rustar-aligner SA is now **byte-for-byte identical** to STAR's SA for the yeast genome (10,862 → 0 entry diffs). diff AS: 6→4 cases (4 remaining are rustar-aligner improvements: .844151 VIII 0mm vs STAR VII 6mm, .4972950 spliced vs unspliced mate2). Phase 17.3 (2026-05-01): PE chimeric detection — `detect_inter_mate_chimeric` in `chimeric/detect.rs`; intra-mate multi-cluster chimeric via cluster splitting + mate2 read_pos adjustment; inter-mate chimeric for discordant pairs (diff chr, same strand, or >1Mb); `align_paired_read` returns 4-tuple including `Vec`; no benchmark regression (8390 both-mapped, 0 half-mapped). Phase 17.11 (2026-05-01): `--chimOutType WithinBAM` — chimeric alignments written as supplementary records (FLAG 0x800) in primary BAM; donor record has full SEQ + SA tag; acceptor has FLAG 0x800 + SA tag + empty SEQ; `build_within_bam_records` in `chimeric/output.rs`; `chim_out_junctions()` / `chim_out_within_bam()` helpers in params.rs; supports mixed `--chimOutType Junctions WithinBAM`. Phase 17.7 (2026-05-01): GTF tag parameters — `--sjdbGTFchrPrefix`, `--sjdbGTFfeatureExon`, `--sjdbGTFtagExonParentTranscript`, `--sjdbGTFtagExonParentGene`; `_configured` variants in `junction/gtf.rs`, `quant/mod.rs`, `quant/transcriptome.rs`, `junction/mod.rs`; all 4 production paths thread params; backward-compat wrappers preserve zero test disruption. Phase 17.9 (2026-05-01): `--outBAMcompression` (BGZF level -1–9, default 1; -1/0=NONE, 1-8=flate2 levels, ≥9=BEST) + `--limitBAMsortRAM` (bytes, 0=unlimited; aborts sort if ~400 bytes/record estimate exceeds limit); `bgzf_compression()` + `make_bgzf_writer()` helpers in `io/bam.rs`; threaded through all 4 BAM writers (unsorted file, sorted file, unsorted stdout, sorted stdout). PE chimericDetectionOld (2026-05-01): per-mate `detect_chimeric_old` called on `all_m1_transcripts` / `all_m2_transcripts` pools after `filter_paired_transcripts` in `read_align.rs`. Phase 17.12 (2026-05-01): BySJout disk buffering — `BySJReadMeta` struct + `NamedTempFile` SAM temp file replaces `Vec`; `create_bysj_writer` / `bysj_write_records` / `bysj_read_n_records` helpers in `io/sam.rs`; `tempfile` moved to `[dependencies]`. Phase 17.13 (2026-05-01): 8 integration tests in `tests/alignment_features.rs` — synthetic 20kb genome with planted GT-AG intron; tests cover BAM output, PE alignment, spliced reads, BySJout, GeneCounts, unmapped output, two-pass mode. Phase 12.2 (2026-05-04): SE chimeric Tier 1b soft-clip re-mapping — `detect_from_soft_clips` in `chimeric/detect.rs` re-seeds the primary alignment's soft-clipped bases when `detect_chimeric_old` finds no partner; `adjust_read_positions` helper shifts sub-seq coords into full-read space for right clips; called as Step 3c in `read_align.rs`. Phase 17.10 (2026-05-04): Chimeric Tier 3 — `detect_from_chimeric_residuals` in `chimeric/detect.rs` re-seeds outer uncovered read regions (before donor / after acceptor) of each found chimeric pair; enables 3-way gene-fusion detection; called as Step 3d in `read_align.rs`. Phase 14.7 (2026-08-11): STARsolo per-read barcode SAM tags. `--outSAMattributes` gains `CR CY UR UY CB UB sM sS sQ` (BAM output only, as in STAR); read-time tags are added in both solo loops via `SoloTagStrings` + `add_solo_barcode_tags`, while `CB`/`UB` follow STAR's readInfo path (`SoloCountRecord.read_index` → `SoloContext.read_info` filled by UMI collapsing → `apply_solo_read_info` rewrites the buffered sorted-BAM records, `-` for uncounted reads), hence the sorted-BAM + gene-feature validation. `umi_correction_map` in `solo/count.rs` supplies the per-method UMI corrections; `1MM_Directional` now counts distinct corrected UMIs off STAR's absorb chain. `--soloType CB_samTagOut` implemented (whitelist correction only, no gene model, no `Solo.out`). Also: per-alignment `gx`/`gn` + `sF` via `align_genes` in `solo/gene.rs` (STAR's `fAlign` sets + `ovType` priority), gene tags on the PE solo path (a pair is one alignment), three-file solo `--readFilesIn` (`cDNA_read1 cDNA_read2 barcode_read`) so PE cDNA works with a separate barcode read, and `--soloBarcodeReadLength` enforced STAR-style (`BarcodeReadLength`: exact-length check by default, N-padding when set to 0). See [ROADMAP.md](ROADMAP.md) for detailed phase tracking and [docs-old/](docs-old/) for per-phase development notes. The published Astro Starlight docs site is in [docs/](docs/). ## Source Layout diff --git a/ROADMAP.md b/ROADMAP.md index d419f612..c1ac18ce 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -308,11 +308,11 @@ See [docs-old/phase17_features.md](docs-old/phase17_features.md) for sub-phase t --- -## Phase 14: STARsolo (Single-Cell) — SUBSTANTIALLY COMPLETE +## Phase 14: STARsolo (Single-Cell) (COMPLETE) **Prerequisite met**: position agreement >99% (SE 99.815% tie-adj, PE 99.883%). Phase unblocked 2026-06-10. -**Status**: A working, STARsolo-faithful single-cell pipeline — Gene count matrix **byte-identical to STARsolo's**, with GeneFull/SJ/Velocyto features, CB_UMI_Simple/Complex/SmartSeq chemistries, EmptyDrops_CR + CellRanger2.2 cell calling, multi-mapper resolution, and `Summary.csv`. Validated in a native three-way benchmark against STARsolo and CellRanger (see below). Remaining: per-record `CB`/`UB`/`GX`/`GN` SAM tags + `CB_samTagOut` (14.7). +**Status**: A STARsolo-faithful single-cell pipeline, Gene count matrix **byte-identical to STARsolo's**, with GeneFull/SJ/Velocyto features, CB_UMI_Simple/Complex/SmartSeq chemistries, EmptyDrops_CR + CellRanger2.2 cell calling, multi-mapper resolution, `Summary.csv`, and the per-read barcode SAM tags. Validated in a native three-way benchmark against STARsolo and CellRanger (see below). Single-cell quantification layered around the existing aligner: the cDNA read aligns through the normal SE path; a paired **barcode read** (R1 = cell barcode + UMI) is parsed, corrected against a whitelist, assigned to a gene, UMI-deduplicated, and emitted as a sparse per-cell count matrix. Target: faithful port of STARsolo (all features). See [docs-old/phase14_starsolo.md](docs-old/phase14_starsolo.md) for the full design and sub-phase tracking. @@ -325,7 +325,7 @@ Single-cell quantification layered around the existing aligner: the cDNA read al | 14.CR | CellRanger 4/5-matching flags (`1MM_CR`, `MultiGeneUMI_CR`, `1MM_multi_Nbase_pseudocounts`, `CellRanger4` clip) | ✅ Complete | | 14.5 | `Summary.csv` (STARsolo-faithful; CellRanger funnel split out) | ✅ Complete | | 14.6 | Cell filtering (`--soloCellFilter`: CellRanger2.2, TopCells, EmptyDrops_CR MC rescue) | ✅ Complete | -| 14.7 | `CB`/`UB`/`GX`/`GN` SAM tags + `CB_samTagOut` | ⬜ Planned | +| 14.7 | `CB`/`UB`/`GX`/`GN` SAM tags + `CB_samTagOut` | ✅ Complete | | 14.8 | More features: GeneFull, SJ, Velocyto (spliced/unspliced/ambiguous) | ✅ Complete | | 14.9 | Multi-gene resolution (`--soloMultiMappers`: Uniform/PropUnique/EM/Rescue) | ✅ Complete | | 14.10 | Other chemistries: CB_UMI_Complex, SmartSeq (SE + PE fragment counts) | ✅ Complete | @@ -347,4 +347,8 @@ Single-cell quantification layered around the existing aligner: the cDNA read al **Phase 14.5–14.11 + performance** (2026-07): completed the feature-parity set — `Summary.csv` (STARsolo-faithful, CellRanger funnel split to its own file), `--soloCellFilter` CellRanger2.2/TopCells/**EmptyDrops_CR** (Monte-Carlo ambient rescue in the `filtered/` writer), `--soloFeatures` **GeneFull/SJ/Velocyto** (spliced/unspliced/ambiguous per Sullivan 2025), `--soloMultiMappers` Uniform/PropUnique/EM/Rescue, chemistries **CB_UMI_Complex** (multi-segment) and **SmartSeq** (plate-based, SE + PE fragment counts), and a rustar-vs-STARsolo SJ + multi-mapper diff harness. Performance: pipelined solo FASTQ decode, parallelized matrix build + EmptyDrops MC, libdeflate/zlib-rs for matrix gzip + BGZF, and an **O(log n + k) segment-tree gene-overlap query** (replacing STAR's linear scan — the #1 solo hotspot, ~14% wall reduction). Sparse suffix array (`--genomeSAsparseD`, byte-identical to STAR's D=2) for a 31% smaller index. 516 tests, 0 clippy warnings. +**Phase 14.7, barcode SAM tags + `CB_samTagOut`** (2026-08-11): the per-read STARsolo tags. `--outSAMattributes` gains `CR CY UR UY CB UB sM sS sQ` alongside the existing `GX GN`, and all of them are emitted in BAM output only, as STAR does (`ReadAlign_outputTranscriptSAM.cpp` drops them from the SAM text path). The read-time tags (`CR`/`CY`/`UR`/`UY`, the `sM` cbMatch code, the whole barcode read in `sS`/`sQ`) go on every record of a read, mapped or unmapped, in both the SE and the `--soloBarcodeMate 1` PE solo loop. `CB`/`UB` follow STAR's readInfo route: count records carry their input read index (`SoloCountRecord.read_index`, free in the struct's padding), UMI collapsing fills `SoloContext.read_info` with each read's cell and corrected UMI, and the buffered sorted-BAM records, tagged with a private `zR` index tag, mirroring STAR's `iReadAll` encoding, are rewritten into `CB`/`UB` before the sort, with `"-"` for reads that were not counted. That ordering is why STAR restricts both tags to `--outSAMtype BAM SortedByCoordinate`, which validation now enforces along with the gene-feature requirement and STAR's refusal of `UB` under `CB_samTagOut`. UMI correction maps were added for every dedup method (`umi_correction_map`: graph components map to their highest-count UMI, directional follows STAR's absorb chain, `1MM_CR` reuses the existing map); `directional` now derives its molecule count from that chain, as STAR does. `--soloType CB_samTagOut` runs as its own mode: whitelist correction only, no gene model, no counting, no `Solo.out`. + +Completing the set: the per-alignment gene tags `gx`/`gn` (every gene of *that* alignment, `;`-joined) and the feature-status tag `sF` (`(overlap type, genes for the read)`, `(-1,-1)` on a sense-strand read whose alignment carries no gene) come from `align_genes` in `solo/gene.rs`, which keeps STAR's per-alignment `fAlign` sets and its `ovType` priority (exonic > exonicAS > intronic > intronicAS > intergenic; the two `GeneFull_Ex50pAS` types cannot arise since that feature is not implemented). The gene tags now also reach the paired-end solo path, where a mate pair is one alignment and both mates carry the pair's genes. `--readFilesIn` accepts STAR's three-file solo layout (`cDNA_read1 cDNA_read2 barcode_read`), so paired-end cDNA works with a separate barcode read for every barcode chemistry, `CB_samTagOut` included; a third file without `--soloType` is refused. `--soloBarcodeReadLength` is now honoured as STAR does (`BarcodeReadLength`): the default expects a barcode read exactly CB+UMI long and treats any other length as a fatal input error, while `0` turns the check off and pads a short read with `N`/`H`, which retires the one divergence this phase had opened. 599 lib + 27 integration tests, 0 clippy warnings. + **Native three-way benchmark** (2026-07, `test/aws/`): fresh single-instance EC2 comparison on a real 10x dataset (`5k_Mouse_PBMCs_5p_gem-x_GEX`, 5′ GEM-X, GRCm39-2024-A), all native x86_64, 10 threads, NVMe, page cache dropped, no BAM. Wall / peak RSS / cells: **STARsolo 2.7.11b** 87 s / 28.3 GB / 4,061; **rustar-aligner** 121 s / 25.7 GB / 3,689 (→ ~105 s with the segment-tree query merged after this run); **rustar `--genomeSAsparseD 2`** 119 s / **17.7 GB** / 3,692; **CellRanger 10.0.0** 347 s / 13.1 GB / 3,858. `Gene/raw` matrix byte-identical to STARsolo's. Supersedes the earlier Docker-emulation numbers above (those were penalized by Rosetta/virtiofs). Remaining gap to STARsolo is small and output-identical; rustar owns the memory frontier via sparse SA. diff --git a/docs-old/phase14_starsolo.md b/docs-old/phase14_starsolo.md index 230190b1..807d4ce2 100644 --- a/docs-old/phase14_starsolo.md +++ b/docs-old/phase14_starsolo.md @@ -54,7 +54,7 @@ two files but is a *single-end alignment* run. | 14.4 | UMI dedup + raw `matrix.mtx` (**MVP complete**) | ✅ Complete | | 14.5 | `Summary.csv` / `Barcodes.stats` / `Features.stats` | ⬜ Planned | | 14.6 | Cell filtering (`filtered/` matrix) | ⬜ Planned | -| 14.7 | `CB`/`UB`/`GX`/`GN` SAM tags + `CB_samTagOut` | ⬜ Planned | +| 14.7 | `CB`/`UB`/`GX`/`GN` SAM tags + `CB_samTagOut` | ✅ Complete | | 14.8 | More features: GeneFull, SJ, Velocyto | ⬜ Planned | | 14.9 | Multi-gene resolution (`--soloMultiMappers`) | ⬜ Planned | | 14.10 | Other chemistries: CB_UMI_Complex, SmartSeq | ⬜ Planned | diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index f749039b..e4974044 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -29,7 +29,7 @@ Run `rustar-aligner --help` for the full machine-generated listing. | Parameter | Default | Description | |-----------|---------|-------------| -| `--readFilesIn` | — | Input FASTQ file(s); second file is mate 2 for paired-end (required for `alignReads`). | +| `--readFilesIn` | — | Input FASTQ file(s); second file is mate 2 for paired-end (required for `alignReads`). A `--soloType` run appends the barcode read as the last file: `cDNA_read barcode_read`, or `cDNA_read1 cDNA_read2 barcode_read` for paired-end cDNA. | | `--readFilesCommand` | — | Decompression command, e.g. `zcat` for `.gz`. | | `--readMapNumber` | `-1` | Number of reads to map (`-1` = all). | | `--clip5pNbases` | `0` | Bases to clip from the 5' end of each mate. | @@ -50,7 +50,7 @@ Run `rustar-aligner --help` for the full machine-generated listing. | Parameter | Default | Description | |-----------|---------|-------------| | `--outSAMstrandField` | `None` | `None` or `intronMotif` (sets XS tag from junction motifs). | -| `--outSAMattributes` | `Standard` | Tags to include: `Standard`, `All`, `None`, or an explicit list (e.g. `NH HI AS NM nM MD`). | +| `--outSAMattributes` | `Standard` | Tags to include: `Standard`, `All`, `None`, or an explicit list (e.g. `NH HI AS NM nM MD`). The STARsolo tags `CR CY UR UY CB UB GX GN gx gn sM sS sQ sF` are also accepted, and (as in STAR) written to BAM output only; `CB`/`UB` additionally need `--outSAMtype BAM SortedByCoordinate`. | | `--outSAMattrRGline` | `-` | Read group line(s). Multiple blocks separated by a literal `,`. | | `--outSAMunmapped` | `None` | Unmapped reads in SAM: `None`, `Within`, or `Within KeepPairs`. | | `--outSAMmapqUnique` | `255` | MAPQ value for uniquely-mapping reads. |