diff --git a/CHANGELOG.md b/CHANGELOG.md index d66f6b9a..8965277e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -65,6 +65,30 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- **`--chimOutType WithinBAM` ported from STAR** (`chimericBAMoutput`, + `alignBAM`). The chimera now replaces the read's normal alignment, as in + STAR: one segment as a normal record and the other as a supplementary + (`0x800`), hard-clipped on the junction side unless `SoftClip`, with `SA` + tags built from the other record's final CIGAR. The representative segment + is chosen as STAR chooses it; paired segments covering both mates are + written as pairs, and mates on either side of the junction as two normal + records. Records carry STAR's flags, MAPQ, mate fields, qualities and + `NH`/`HI`/`AS`/`nM`/`NM`. Previously the normal alignment was kept and two + incomplete records were appended, so no chimeric read matched STAR's + output; on yeast 10k every chimeric read now does, apart from known ties. + As in STAR, such reads are not counted as mapped in `Log.final.out` and + feed neither `SJ.out.tab` nor gene counts. + - `WithinBAM` requires `--outSAMtype BAM` and adds `NM` to the output + attributes, as STAR does. + - `--chimOutType` accepts up to three values (e.g. `Junctions WithinBAM + HardClip`) and rejects unknown ones. `SeparateSAMold` is refused rather + than silently ignored; it is not implemented. + +- **Paired-end `nM` is the pair's mismatch count on both mates**, as STAR + writes it (`trOut.nMM`, like `AS`). We wrote each mate's own count, which + differed from STAR on 42% of yeast PE pairs; no comparison script checked + tags. + - **`--quantMode TranscriptomeSAM` no longer aborts on reads soft-clipped at both ends of a single match block** (e.g. `1S97M2S`). With the default `--quantTranscriptomeSAMoutput BanSingleEnd_BanIndels_ExtendSoftclip`, the diff --git a/docs/src/content/docs/guides/chimeric.md b/docs/src/content/docs/guides/chimeric.md index a30512a4..f99694df 100644 --- a/docs/src/content/docs/guides/chimeric.md +++ b/docs/src/content/docs/guides/chimeric.md @@ -5,7 +5,7 @@ description: Detect reads spanning two distant genomic locations — fusion cand A chimeric alignment is a read whose two halves map to different genomic locations — different chromosomes, the same chromosome with an unrealistically large gap, or the same chromosome but on opposite strands. These are the candidate evidence for gene fusions, large-scale structural variants, and circular RNA back-splicing. -rustar-aligner implements STAR's chimeric detection pipeline with four tiers, all of which run automatically when chimeric detection is enabled. +rustar-aligner uses STAR's two chimeric detectors, ported as written; which one runs depends on `--chimMultimapNmax`. ## Enabling chimeric detection @@ -25,14 +25,12 @@ Higher values (e.g. 20) produce fewer, more confident calls; lower values (e.g. ## What gets detected -The 4-tier pipeline runs in this order, stopping as soon as a chimeric pair is found: +Both detectors work on the read's window transcripts, the alignments stitching produced before any filtering. For paired-end reads these are whole-fragment transcripts (`mate1 | RC(mate2)`), so one side of a chimera can be a complete mate pair. -1. **Tier 1 — transcript-pair search.** Searches the read's existing transcript pool for two segments that together cover most of the read but map to incompatible locations. -2. **Tier 2 — multi-cluster.** When the seed pool produces multiple distinct alignment clusters, evaluates each pair as a candidate chimeric. -3. **Tier 1b — soft-clip re-mapping.** Takes the soft-clipped tail of the primary alignment and re-seeds it against the genome. Recovers chimeric pairs where the original aligner only kept one half. -4. **Tier 3 — residual outer re-seeding.** For reads where Tier 1/1b/2 has already found a chimeric pair, re-seeds the *remaining* uncovered regions of the read (before the donor / after the acceptor). Enables 3-way detection — gene fusions involving three loci, e.g. a complex rearrangement where two breakpoints are present in a single read. +- **Default (`--chimMultimapNmax 0`):** STAR's `chimericDetectionOld`. The best linear alignment is one segment, and the detector looks for the single best partner that covers the rest of the read. A chimera is reported only if it is clearly better than the runner-up (`--chimScoreSeparation`), and at most one per read. +- **`--chimMultimapNmax N > 0`:** STAR's `chimericDetectionMult`. Every pair of window transcripts is tried and re-scored after the junction is placed, and all chimeras within `--chimMultimapScoreRange` of the best are reported, up to `N` (more than `N` reports none). The junction file gains six columns: number of chimeras, maximum possible score, best non-chimeric score, this chimera's score, best chimeric score, and whether the mates were merged. -For paired-end data, additional inter-mate detection runs: if mate 1 maps confidently to one location and mate 2 to another that's incompatible with a normal proper pair (different chromosomes, same strand, or >1 Mb gap), the pair is reported as chimeric. +When the junction falls inside a read, its position is chosen by STAR's scan over the overlap, preferring GT/AG or CT/AC motifs; when the two mates of a pair lie on either side of it, the junction type is `-1`. ## Output formats @@ -49,10 +47,13 @@ Writes a `Chimeric.out.junction` file with one row per chimeric junction ### `WithinBAM` ```bash ---chimOutType WithinBAM +--chimOutType WithinBAM [HardClip|SoftClip] +--outSAMtype BAM Unsorted # or SortedByCoordinate; WithinBAM needs BAM output ``` -Embeds the chimeric segments as supplementary alignment records (FLAG `0x800`) in the primary BAM, with `SA` tags linking the donor and acceptor halves. This is the format expected by tools that process chimeric BAM directly (e.g. for fusion calling on already-sorted BAMs). +Writes the chimera into the main BAM **in place of the read's normal alignment**, as STAR does: one segment as a normal record and the other as a supplementary record (FLAG `0x800`), with `SA` tags pointing at each other. The supplementary segment is hard-clipped by default (`HardClip`) or soft-clipped (`SoftClip`). For paired-end reads a segment that covers both mates is written as a normal pair; when the mates lie on either side of the junction, both are written as normal records with no supplementary. `NM` is added to the output attributes, as STAR does. This is the input [Arriba](https://github.com/suhrig/arriba) recommends. + +Because these reads are written as chimeras, STAR does not count them as uniquely or multi-mapped in `Log.final.out`, and they do not contribute to `SJ.out.tab` or gene counts; rustar-aligner does the same. ### Mixed output @@ -60,7 +61,7 @@ Embeds the chimeric segments as supplementary alignment records (FLAG `0x800`) i --chimOutType Junctions WithinBAM ``` -Writes both the junction file and the supplementary BAM records. Useful when downstream tools have different format requirements. +Writes both the junction file and the chimeric BAM records. ## Tuning parameters @@ -87,11 +88,11 @@ The `Chimeric.out.junction` file has 14 tab-separated columns (STAR-compatible): | 4 | `chr_acceptorB` | Acceptor (right segment) chromosome | | 5 | `brkpt_acceptorB` | Acceptor breakpoint position | | 6 | `strand_acceptorB` | Acceptor strand | -| 7 | `junction_type` | -1 (encompassing PE) / 0 (non-canonical) / 1 (GT/AG) / 2 (CT/AC) | +| 7 | `junction_type` | -1 (mates on either side of the junction) / 0 (non-canonical) / 1 (GT/AG) / 2 (CT/AC) | | 8 | `repeat_left_lenA` | Length of repeat to the left | | 9 | `repeat_right_lenB` | Length of repeat to the right | | 10 | `read_name` | Source read name | -| 11 | `start_alnA` | Donor start position on the read | -| 12 | `cigar_alnA` | Donor CIGAR | -| 13 | `start_alnB` | Acceptor start position on the read | -| 14 | `cigar_alnB` | Acceptor CIGAR | +| 11 | `start_alnA` | Donor segment's first aligned genomic position | +| 12 | `cigar_alnA` | Donor CIGAR; a paired segment carries the second mate after a `p` operation (the genomic gap between mates, negative when they overlap) | +| 13 | `start_alnB` | Acceptor segment's first aligned genomic position | +| 14 | `cigar_alnB` | Acceptor CIGAR, same convention | diff --git a/docs/src/content/docs/reference/cli-parameters.md b/docs/src/content/docs/reference/cli-parameters.md index bbdd0f48..ddfd10e4 100644 --- a/docs/src/content/docs/reference/cli-parameters.md +++ b/docs/src/content/docs/reference/cli-parameters.md @@ -173,7 +173,7 @@ Run `rustar-aligner --help` for the full machine-generated listing. | `--chimSegmentReadGapMax` | `0` | Max read-space gap between chimeric segments. | | `--chimJunctionOverhangMin` | `20` | Min overhang at chimeric junction. | | `--chimScoreJunctionNonGTAG` | `-1` | Score penalty for non-GT/AG chimeric junctions. | -| `--chimOutType` | `Junctions` | `Junctions`, `WithinBAM`, or both space-separated. | +| `--chimOutType` | `Junctions` | `Junctions` and/or `WithinBAM`, optionally followed by `HardClip` (default) or `SoftClip` for the supplementary segment. `WithinBAM` requires `--outSAMtype BAM`. | ## Debug diff --git a/docs/src/content/docs/reference/output-files.md b/docs/src/content/docs/reference/output-files.md index 0d62a67b..85a5324d 100644 --- a/docs/src/content/docs/reference/output-files.md +++ b/docs/src/content/docs/reference/output-files.md @@ -71,7 +71,7 @@ Written when `--chimSegmentMin > 0`. Format depends on `--chimOutType`. ### Chimeric records in the primary BAM -When `--chimOutType` includes `WithinBAM`, the chimeric segments are embedded as supplementary alignment records (FLAG `0x800`) in the main BAM output, with `SA` tags linking the donor and acceptor halves. Tools like Arriba and STAR-Fusion know how to read either format. +When `--chimOutType` includes `WithinBAM`, a chimeric read is written into the main BAM in place of its normal alignment, as STAR does: one segment as a normal record and the other as a supplementary record (FLAG `0x800`, hard-clipped unless `SoftClip`), with `SA` tags linking the two. Paired reads whose mates lie on either side of the junction are written as two normal records. See the [chimeric guide](/rustar-aligner/guides/chimeric/). ## Quantification diff --git a/src/chimeric/bam.rs b/src/chimeric/bam.rs new file mode 100644 index 00000000..4f92a5cb --- /dev/null +++ b/src/chimeric/bam.rs @@ -0,0 +1,780 @@ +//! `--chimOutType WithinBAM`: chimeric records in the main BAM. +//! +//! A port of `ChimericAlign::chimericBAMoutput` and the mapped branch of +//! `ReadAlign::alignBAM`, which it calls once per chimeric segment. The segments +//! are STAR's whole transcripts ([`ChimBam`]), so a segment that covers both +//! mates of a pair produces two records, exactly as STAR's does. +//! +//! STAR writes these in place of the read's normal alignments +//! (`ReadAlign_oneRead.cpp:99`); the caller is responsible for that. + +use crate::chimeric::detect::{ChimExon, WinTr}; +use crate::chimeric::{ChimBam, ChimericAlignment}; +use crate::error::Error; +use crate::genome::Genome; +use crate::io::fastq::{complement_base, decode_base}; +use crate::io::sam::{apply_sam_flag_or_and, fastq_qual_to_phred, maybe_insert_rg_tag}; +use crate::params::{Parameters, SamAttributes}; +use bstr::BString; +use noodles::sam; +use noodles::sam::alignment::record::MappingQuality; +use noodles::sam::alignment::record::cigar::{Op, op::Kind}; +use noodles::sam::alignment::record::data::field::Tag; +use noodles::sam::alignment::record_buf::data::field::Value; +use noodles::sam::alignment::record_buf::data::field::value::Array; +use noodles::sam::alignment::record_buf::{QualityScores, RecordBuf, Sequence}; + +/// STAR's `canonSJ` for the gap between mates. +const SJ_MATE_GAP: i32 = -3; +/// STAR's `canonSJ` for a deletion. +const SJ_DELETION: i32 = -1; + +/// The read a chimera came from, as sequenced. +pub struct ChimReadInput<'a> { + pub name: &'a str, + /// Mates as sequenced (encoded bases), before clipping. `[1]` is empty for + /// single-end. + pub seq: [&'a [u8]; 2], + /// FASTQ qualities (Phred+33) for the same bases. + pub qual: [&'a [u8]; 2], + /// Bases clipped before alignment, per mate: `[5', 3']`. + pub clip: [[usize; 2]; 2], +} + +/// STAR's `alignType` for a chimeric record. +#[derive(Clone, Copy, PartialEq, Eq)] +enum AlignType { + /// `-10`: written as a normal alignment. + Normal, + /// `-11`: supplementary, hard-clipped on the left (junction on the left). + HardLeft, + /// `-12`: supplementary, hard-clipped on the right. + HardRight, + /// `-13`: supplementary, soft-clipped. + Soft, +} + +/// Mate information passed to `alignBAM` for a single-mate segment. +#[derive(Clone, Copy)] +struct MateRef { + chr: usize, + /// Absolute genome position of the mate's first base. + start: u64, + /// The mate aligns to the reverse strand. + reverse: bool, +} + +/// Every chimeric BAM record for one read, in STAR's order. +pub fn build_chimeric_bam_records( + chims: &[ChimericAlignment], + read: &ChimReadInput, + genome: &Genome, + params: &Parameters, +) -> Result, Error> { + let mut out = Vec::new(); + for chim in chims { + if let Some(b) = &chim.bam { + out.extend(chimeric_bam_output(b, read, genome, params)?); + } + } + Ok(out) +} + +fn first_frag(t: &WinTr) -> u8 { + t.exons[0].frag +} + +fn last_frag(t: &WinTr) -> u8 { + t.exons[t.exons.len() - 1].frag +} + +/// `ChimericAlign::chimericBAMoutput` (`ChimericAlign_chimericBAMoutput.cpp`). +fn chimeric_bam_output( + b: &ChimBam, + read: &ChimReadInput, + genome: &Genome, + params: &Parameters, +) -> Result, Error> { + let tr = &b.tr; + // Which segment is written as the representative (normal) alignment. + // chimType 1: one segment covers both mates; 2: the mates bracket the + // junction (both written as normal); 3: both segments in one mate (SE). + let (represent, chim_type): (Option, u8) = if first_frag(&tr[0]) != last_frag(&tr[0]) { + (Some(0), 1) + } else if first_frag(&tr[1]) != last_frag(&tr[1]) { + (Some(1), 1) + } else if first_frag(&tr[0]) != first_frag(&tr[1]) { + (None, 2) + } else { + (Some(usize::from(b.max_score[0] <= b.max_score[1])), 3) + }; + let hard_clip = params.chim_out_bam_hard_clip(); + + let mut recs: Vec = Vec::new(); + let (mut i_repr, mut i_suppl): (Option, Option) = (None, None); + for itr in 0..2 { + let other = &tr[1 - itr]; + let (mate, align_type); + if chim_type == 2 { + mate = Some(MateRef { + chr: other.chr, + start: other.exons[0].g, + reverse: other.str_ != other.exons[0].frag, + }); + align_type = AlignType::Normal; + } else if represent == Some(itr) { + mate = None; + align_type = AlignType::Normal; + // For a two-mate representative, the chimerically split mate is + // its second record. + i_repr = Some(recs.len() + usize::from(first_frag(&tr[itr]) != first_frag(other))); + } else { + align_type = if !hard_clip { + AlignType::Soft + } else if itr % 2 == tr[itr].str_ as usize { + AlignType::HardRight + } else { + AlignType::HardLeft + }; + i_suppl = Some(recs.len()); + mate = if chim_type == 1 { + // Mate info for the supplementary: the representative's mate + // that is not this segment's. + let r = &tr[represent.expect("chimType 1 has a representative")]; + let mut iex = 0; + while iex < r.exons.len() - 1 && r.exons[iex].frag == first_frag(&tr[itr]) { + iex += 1; + } + Some(MateRef { + chr: r.chr, + start: r.exons[iex].g, + reverse: r.str_ != r.exons[iex].frag, + }) + } else { + None + }; + } + recs.extend(align_bam( + &tr[itr], + b.max_score[itr], + b.n_mm[itr], + b, + mate, + align_type, + read, + genome, + params, + )?); + } + + // SA tags: the representative and the supplementary point at each other, + // using the other record's final CIGAR, MAPQ and NM. + if let (Some(r), Some(s)) = (i_repr, i_suppl) { + let sa_for_repr = sa_tag(&recs[s], genome); + let sa_for_suppl = sa_tag(&recs[r], genome); + recs[r] + .data_mut() + .insert(Tag::new(b'S', b'A'), Value::String(sa_for_repr)); + recs[s] + .data_mut() + .insert(Tag::new(b'S', b'A'), Value::String(sa_for_suppl)); + } + Ok(recs) +} + +/// `chr,pos,strand,CIGAR,MAPQ,NM;` from a finished record, as STAR builds it +/// with `bam_cigarString` and the record's `NM` (always present under +/// WithinBAM, which adds the attribute). +fn sa_tag(rec: &RecordBuf, genome: &Genome) -> BString { + use std::fmt::Write; + let chr = rec + .reference_sequence_id() + .map_or("*", |i| genome.chr_name[i].as_str()); + let pos = rec.alignment_start().map_or(0, usize::from); + let strand = if rec.flags().is_reverse_complemented() { + '-' + } else { + '+' + }; + let mut cigar = String::new(); + for op in rec.cigar().as_ref() { + let c = match op.kind() { + Kind::Match => 'M', + Kind::Insertion => 'I', + Kind::Deletion => 'D', + Kind::Skip => 'N', + Kind::SoftClip => 'S', + Kind::HardClip => 'H', + Kind::Pad => 'P', + Kind::SequenceMatch => '=', + Kind::SequenceMismatch => 'X', + }; + let _ = write!(cigar, "{}{c}", op.len()); + } + let mapq = rec.mapping_quality().map_or(255, u8::from); + let nm = match rec.data().get(&Tag::EDIT_DISTANCE) { + Some(v) => v.as_int().unwrap_or(0), + None => 0, + }; + format!("{chr},{pos},{strand},{cigar},{mapq},{nm};").into() +} + +/// The combined read in a transcript's strand frame (`Read1[0]` or `Read1[2]`), +/// over the clipped mates the exons index. +fn strand_frame(read: &ChimReadInput, b: &ChimBam, str_: u8) -> Vec { + let clipped = |m: usize| -> &[u8] { + let s = read.seq[m]; + let [c5, c3] = read.clip[m]; + if s.len() < c5 + c3 { + &[] + } else { + &s[c5..s.len() - c3] + } + }; + let mut fwd: Vec = clipped(0).to_vec(); + if b.paired { + fwd.push(crate::align::stitch::PE_SPACER_BASE); + fwd.extend(clipped(1).iter().rev().map(|&x| complement_base(x))); + } + if str_ == 0 { + fwd + } else { + fwd.iter().rev().map(|&x| complement_base(x)).collect() + } +} + +/// `ReadAlign::alignBAM` for a mapped chimeric transcript: one record per mate +/// the transcript covers. +#[allow(clippy::too_many_arguments)] +fn align_bam( + t: &WinTr, + max_score: i32, + n_mm: u32, + b: &ChimBam, + mate: Option, + align_type: AlignType, + read: &ChimReadInput, + genome: &Genome, + params: &Parameters, +) -> Result, Error> { + let n_ex = t.exons.len(); + let (mut n_mates, mut i_ex_mate) = (1usize, n_ex - 1); + for i in 0..n_ex - 1 { + if t.canon_sj[i] == SJ_MATE_GAP { + n_mates = 2; + i_ex_mate = i; + break; + } + } + let rl = b.read_length; + let rl_orig = [read.seq[0].len(), read.seq[1].len()]; + let left_mate = if b.paired { t.str_ as usize } else { 0 }; + let chr_start = genome.chr_start[t.chr]; + let frame = strand_frame(read, b, t.str_); + let attrs = params.out_sam_attributes; + let rg_owned = params.primary_rg_id()?; + + let mapq = match b.chim_n { + n if n >= 5 => 0, + n if n >= 3 => 1, + 2 => 3, + _ => params.out_sam_mapq_unique, + }; + + let mut out = Vec::with_capacity(n_mates); + for imate in 0..n_mates { + let (i_ex1, i_ex2) = if imate == 0 { + (0, i_ex_mate) + } else { + (i_ex_mate + 1, n_ex - 1) + }; + let e1 = t.exons[i_ex1]; + let e2 = t.exons[i_ex2]; + let mate_idx = e1.frag as usize; + let str_ = u16::from(t.str_); + + // FLAG + let mut flag: u16 = 0; + if b.paired { + flag = 0x1; + if i_ex_mate == n_ex - 1 { + if mate.is_none() { + flag |= 0x8; + } + } else { + flag |= 0x2; + } + } + if align_type != AlignType::Normal { + flag |= 0x800; + } + if !b.is_best { + flag |= 0x100; + } + if mate_idx == 0 { + flag |= str_ * 0x10; + if n_mates == 2 { + flag |= (1 - str_) * 0x20; + } + } else { + flag |= (1 - str_) * 0x10; + if n_mates == 2 { + flag |= str_ * 0x20; + } + } + if b.paired { + flag |= if mate_idx == 0 { 0x40 } else { 0x80 }; + if n_mates == 1 && mate.is_some_and(|m| m.reverse) { + flag |= 0x20; + } + } + + // Clipped bases on this record's left, in genome orientation. + let trim_l = match (t.str_, mate_idx) { + (0, 0) => read.clip[0][0], + (0, _) => read.clip[1][1], + (_, 0) => read.clip[0][1], + (_, _) => read.clip[1][0], + }; + + // CIGAR + let mut cigar: Vec = Vec::new(); + let trim_l1 = trim_l + e1.r + - if e1.r < rl[left_mate] { + 0 + } else { + rl[left_mate] + 1 + }; + if trim_l1 > 0 { + let k = if align_type == AlignType::HardLeft { + Kind::HardClip + } else { + Kind::SoftClip + }; + cigar.push(Op::new(k, trim_l1)); + } + let mut sj_motif: Vec = Vec::new(); + let mut sj_intron: Vec = Vec::new(); + for ii in i_ex1..=i_ex2 { + let x = t.exons[ii]; + if ii > i_ex1 { + let p = t.exons[ii - 1]; + let gap_g = x.g - (p.g + p.l as u64); + let gap_r = x.r - p.r - p.l; + if gap_r > 0 { + cigar.push(Op::new(Kind::Insertion, gap_r)); + } + if t.canon_sj[ii - 1] >= 0 || t.sj_annot[ii - 1] { + cigar.push(Op::new(Kind::Skip, gap_g as usize)); + sj_motif + .push((t.canon_sj[ii - 1] + if t.sj_annot[ii - 1] { 20 } else { 0 }) as i8); + sj_intron.push((p.g + p.l as u64 + 1 - chr_start) as i32); + sj_intron.push((x.g - chr_start) as i32); + } else if gap_g > 0 { + cigar.push(Op::new(Kind::Deletion, gap_g as usize)); + } + } + if x.l > 0 { + cigar.push(Op::new(Kind::Match, x.l)); + } + } + if sj_motif.is_empty() { + sj_motif.push(-1); + sj_intron.push(-1); + } + let trim_r1 = (if e1.r < rl[left_mate] { + rl_orig[left_mate] + } else { + rl[left_mate] + 1 + rl_orig[mate_idx] + }) as i64 + - (e2.r + e2.l) as i64 + - trim_l as i64; + if trim_r1 > 0 { + let k = if align_type == AlignType::HardRight { + Kind::HardClip + } else { + Kind::SoftClip + }; + cigar.push(Op::new(k, trim_r1 as usize)); + } + + let mut rec = RecordBuf::default(); + rec.name_mut().replace(read.name.into()); + *rec.flags_mut() = sam::alignment::record::Flags::from(flag); + *rec.reference_sequence_id_mut() = Some(t.chr); + let pos = (e1.g - chr_start + 1) as usize; + *rec.alignment_start_mut() = Some( + pos.try_into() + .map_err(|e| Error::Alignment(format!("invalid chimeric position {pos}: {e}")))?, + ); + *rec.mapping_quality_mut() = MappingQuality::new(mapq); + *rec.cigar_mut() = cigar.into_iter().collect(); + + // RNEXT / PNEXT / TLEN + if n_mates > 1 { + *rec.mate_reference_sequence_id_mut() = Some(t.chr); + let other = t.exons[if imate == 0 { i_ex_mate + 1 } else { 0 }]; + let mpos = (other.g - chr_start + 1) as usize; + *rec.mate_alignment_start_mut() = mpos.try_into().ok(); + let last = t.exons[n_ex - 1]; + let tlen = if params.out_sam_tlen == 2 { + let mate1_end = t.exons[i_ex_mate].g + t.exons[i_ex_mate].l as u64; + let end = (last.g + last.l as u64).max(mate1_end); + let start = t.exons[0].g.min(t.exons[i_ex_mate + 1].g); + let left_most = usize::from(t.exons[0].g > t.exons[i_ex_mate + 1].g); + let v = (end - start) as i32; + if imate == left_most { v } else { -v } + } else { + let v = (last.g + last.l as u64 - t.exons[0].g) as i32; + if imate == 0 { v } else { -v } + }; + *rec.template_length_mut() = tlen; + } else if let Some(m) = mate { + *rec.mate_reference_sequence_id_mut() = Some(m.chr); + let mpos = (m.start - genome.chr_start[m.chr] + 1) as usize; + *rec.mate_alignment_start_mut() = mpos.try_into().ok(); + } + + // SEQ / QUAL: the mate as sequenced, reverse-complemented when this + // record is on the other strand, then trimmed where hard-clipped. + let seq = read.seq[mate_idx]; + let qual = read.qual[mate_idx]; + let same_strand = mate_idx == t.str_ as usize; + let mut seq_out: Vec = if same_strand { + seq.iter().map(|&x| decode_base(x)).collect() + } else { + seq.iter() + .rev() + .map(|&x| decode_base(complement_base(x))) + .collect() + }; + let mut qual_out = fastq_qual_to_phred(qual); + if !same_strand { + qual_out.reverse(); + } + match align_type { + AlignType::HardLeft => { + let k = trim_l1.min(seq_out.len()); + seq_out.drain(..k); + if qual_out.len() >= k { + qual_out.drain(..k); + } + } + AlignType::HardRight => { + let k = (trim_r1.max(0) as usize).min(seq_out.len()); + seq_out.truncate(seq_out.len() - k); + if qual_out.len() >= k { + qual_out.truncate(qual_out.len() - k); + } + } + _ => {} + } + *rec.sequence_mut() = Sequence::from(seq_out); + *rec.quality_scores_mut() = QualityScores::from(qual_out); + + // Attributes, in the order the main writer uses. + let (nm, md) = nm_md(t, i_ex1, i_ex2, &frame, genome); + let data = rec.data_mut(); + if attrs.contains(SamAttributes::NH) { + data.insert(Tag::ALIGNMENT_HIT_COUNT, Value::from(b.chim_n as i32)); + } + if attrs.contains(SamAttributes::HI) { + data.insert( + Tag::HIT_INDEX, + Value::from(b.i_tr as i32 + params.out_sam_attr_ih_start as i32), + ); + } + if attrs.contains(SamAttributes::AS) { + data.insert(Tag::ALIGNMENT_SCORE, Value::from(max_score)); + } + if attrs.contains(SamAttributes::NMM) { + data.insert(Tag::new(b'n', b'M'), Value::from(n_mm as i32)); + } + if attrs.contains(SamAttributes::NM) { + data.insert(Tag::EDIT_DISTANCE, Value::from(nm as i32)); + } + if attrs.contains(SamAttributes::XS) { + let (p, m) = (t.intron_motifs[1], t.intron_motifs[2]); + if p > 0 && m == 0 { + data.insert(Tag::new(b'X', b'S'), Value::Character(b'+')); + } else if m > 0 && p == 0 { + data.insert(Tag::new(b'X', b'S'), Value::Character(b'-')); + } + } + if attrs.contains(SamAttributes::JM) { + data.insert(Tag::new(b'j', b'M'), Value::Array(Array::Int8(sj_motif))); + } + if attrs.contains(SamAttributes::JI) { + data.insert(Tag::new(b'j', b'I'), Value::Array(Array::Int32(sj_intron))); + } + if attrs.contains(SamAttributes::MD) { + data.insert(Tag::new(b'M', b'D'), Value::String(md.into())); + } + maybe_insert_rg_tag(&mut rec, rg_owned.as_deref()); + apply_sam_flag_or_and(&mut rec, params); + out.push(rec); + } + Ok(out) +} + +/// `ReadAlign::samAttrNM_MD`: edit distance and MD string over one mate's +/// exons. Unlike `alignScore`, an `N` on either side counts as a mismatch. +fn nm_md(t: &WinTr, i_ex1: usize, i_ex2: usize, frame: &[u8], genome: &Genome) -> (u32, String) { + use std::fmt::Write; + let nt = |b: u8| -> char { + match b { + 0 => 'A', + 1 => 'C', + 2 => 'G', + 3 => 'T', + _ => 'N', + } + }; + let mut md = String::new(); + let (mut match_n, mut n_mm, mut n_i, mut n_d) = (0u32, 0u32, 0u32, 0u32); + for iex in i_ex1..=i_ex2 { + let x: ChimExon = t.exons[iex]; + for ii in 0..x.l { + let r1 = frame.get(x.r + ii).copied().unwrap_or(4); + let g1 = genome.get_base(x.g + ii as u64).unwrap_or(4); + if r1 != g1 || r1 == 4 || g1 == 4 { + n_mm += 1; + let _ = write!(md, "{match_n}{}", nt(g1)); + match_n = 0; + } else { + match_n += 1; + } + } + if iex < i_ex2 { + let y = t.exons[iex + 1]; + if t.canon_sj[iex] < 0 { + n_d += (y.g - (x.g + x.l as u64)) as u32; + } + n_i += (y.r - x.r - x.l) as u32; + if t.canon_sj[iex] == SJ_DELETION { + let _ = write!(md, "{match_n}^"); + for g in x.g + x.l as u64..y.g { + md.push(nt(genome.get_base(g).unwrap_or(4))); + } + match_n = 0; + } + } + } + let _ = write!(md, "{match_n}"); + (n_mm + n_i + n_d, md) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::chimeric::detect::{ChimExon, WinTr}; + + const CHR_PAD: u64 = 1024; + + /// Two 1000-base chromosomes of pseudo-random sequence. + fn genome() -> Genome { + let n_genome = CHR_PAD * 2; + let mut seq = vec![5u8; 2 * n_genome as usize]; + let mut x: u64 = 0x9E37_79B9_7F4A_7C15; + for chr in 0..2u64 { + for i in 0..1000u64 { + x = x.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1); + seq[(chr * CHR_PAD + i) as usize] = ((x >> 33) % 4) as u8; + } + } + Genome { + transform_blocks: None, + sequence: seq.into(), + n_genome, + n_genome_real: n_genome, + n_chr_real: 2, + chr_name: vec!["chr0".into(), "chr1".into()], + chr_length: vec![1000, 1000], + chr_start: vec![0, CHR_PAD, n_genome], + } + } + + fn slice(g: &Genome, start: u64, len: usize) -> Vec { + (0..len as u64) + .map(|i| g.get_base(start + i).unwrap()) + .collect() + } + + fn params() -> Parameters { + Parameters::parse_from([ + "rustar-aligner", + "--readFilesIn", + "r.fq", + "--chimSegmentMin", + "12", + "--chimOutType", + "WithinBAM", + "--outSAMtype", + "BAM", + "Unsorted", + ]) + } + + fn tr(chr: usize, str_: u8, g: u64, r: usize, l: usize, frag: u8) -> WinTr { + WinTr { + chr, + str_, + exons: vec![ChimExon { g, r, l, frag }], + canon_sj: Vec::new(), + sj_annot: Vec::new(), + max_score: 0, + intron_motifs: [0; 3], + r_length: l, + ro_start: r, + } + } + + fn tag(rec: &RecordBuf, t: [u8; 2]) -> Option<&Value> { + rec.data().get(&Tag::new(t[0], t[1])) + } + + fn cigar_string(rec: &RecordBuf) -> String { + let mut s = String::new(); + for op in rec.cigar().as_ref() { + let c = match op.kind() { + Kind::Match => 'M', + Kind::SoftClip => 'S', + Kind::HardClip => 'H', + _ => '?', + }; + let _ = std::fmt::Write::write_fmt(&mut s, format_args!("{}{c}", op.len())); + } + s + } + + /// A single-end chimera, half chr0 and half chr1: the higher-scoring + /// segment is written as the normal record, the other as a hard-clipped + /// supplementary with only its own bases in SEQ, and each SA tag carries + /// the other record's final CIGAR. + #[test] + fn single_end_chimera_is_representative_plus_hard_clipped_supplementary() { + let g = genome(); + let mut read = slice(&g, 100, 60); + read.extend(slice(&g, CHR_PAD + 500, 60)); + let qual = vec![b'I'; 120]; + let b = ChimBam { + tr: [tr(0, 0, 100, 0, 60, 0), tr(1, 0, CHR_PAD + 500, 60, 60, 0)], + max_score: [57, 59], + n_mm: [0, 0], + i_tr: 0, + chim_n: 1, + is_best: true, + read_length: [120, 0], + paired: false, + }; + let input = ChimReadInput { + name: "r", + seq: [&read, &[]], + qual: [&qual, &[]], + clip: [[0, 0], [0, 0]], + }; + let recs = chimeric_bam_output(&b, &input, &g, ¶ms()).unwrap(); + assert_eq!(recs.len(), 2); + // Segment 0 scores lower, so it is the supplementary; it leads in the + // read and is forward, so its right end (the junction) is hard-clipped. + let (suppl, repr) = (&recs[0], &recs[1]); + assert_eq!(u16::from(suppl.flags()), 0x800); + assert_eq!(cigar_string(suppl), "60M60H"); + assert_eq!(suppl.sequence().len(), 60); + assert_eq!(u16::from(repr.flags()), 0); + assert_eq!(cigar_string(repr), "60S60M"); + assert_eq!(repr.sequence().len(), 120); + // WithinBAM adds NM; SA points at the other record's final CIGAR. + assert_eq!(tag(repr, *b"NM"), Some(&Value::from(0))); + assert_eq!( + tag(repr, *b"SA"), + Some(&Value::String("chr0,101,+,60M60H,255,0;".into())) + ); + assert_eq!( + tag(suppl, *b"SA"), + Some(&Value::String("chr1,501,+,60S60M,255,0;".into())) + ); + assert_eq!(tag(repr, *b"AS"), Some(&Value::from(59))); + } + + /// Mates on two chromosomes bracket the junction (chimType 2): both are + /// written as normal paired records pointing at each other, with no SA. + #[test] + fn bracketing_mates_are_two_normal_records_without_sa() { + let g = genome(); + let m1 = slice(&g, 100, 60); + let m2_rc = slice(&g, CHR_PAD + 500, 60); + let m2: Vec = m2_rc.iter().rev().map(|&x| 3 - x).collect(); + let q = vec![b'I'; 60]; + let b = ChimBam { + tr: [tr(0, 0, 100, 0, 60, 0), tr(1, 0, CHR_PAD + 500, 61, 60, 1)], + max_score: [59, 59], + n_mm: [0, 0], + i_tr: 0, + chim_n: 1, + is_best: true, + read_length: [60, 60], + paired: true, + }; + let input = ChimReadInput { + name: "r", + seq: [&m1, &m2], + qual: [&q, &q], + clip: [[0, 0], [0, 0]], + }; + let recs = chimeric_bam_output(&b, &input, &g, ¶ms()).unwrap(); + assert_eq!(recs.len(), 2); + // Mate1 forward, mate reverse: 0x1|0x20|0x40. Mate2 reverse: 0x1|0x10|0x80. + assert_eq!(u16::from(recs[0].flags()), 0x1 | 0x20 | 0x40); + assert_eq!(u16::from(recs[1].flags()), 0x1 | 0x10 | 0x80); + assert_eq!(recs[0].mate_reference_sequence_id(), Some(1)); + assert_eq!(recs[0].mate_alignment_start().map(usize::from), Some(501)); + assert_eq!(recs[1].mate_reference_sequence_id(), Some(0)); + assert_eq!(recs[1].mate_alignment_start().map(usize::from), Some(101)); + assert!(recs.iter().all(|r| tag(r, *b"SA").is_none())); + assert_eq!(cigar_string(&recs[0]), "60M"); + assert_eq!(cigar_string(&recs[1]), "60M"); + } + + /// `SoftClip` keeps the supplementary's clipped bases in SEQ. + #[test] + fn soft_clip_mode_keeps_the_whole_read() { + let g = genome(); + let mut read = slice(&g, 100, 60); + read.extend(slice(&g, CHR_PAD + 500, 60)); + let qual = vec![b'I'; 120]; + let b = ChimBam { + tr: [tr(0, 0, 100, 0, 60, 0), tr(1, 0, CHR_PAD + 500, 60, 60, 0)], + max_score: [57, 59], + n_mm: [0, 0], + i_tr: 0, + chim_n: 1, + is_best: true, + read_length: [120, 0], + paired: false, + }; + let input = ChimReadInput { + name: "r", + seq: [&read, &[]], + qual: [&qual, &[]], + clip: [[0, 0], [0, 0]], + }; + let p = Parameters::parse_from([ + "rustar-aligner", + "--readFilesIn", + "r.fq", + "--chimSegmentMin", + "12", + "--chimOutType", + "WithinBAM", + "SoftClip", + "--outSAMtype", + "BAM", + "Unsorted", + ]); + let recs = chimeric_bam_output(&b, &input, &g, &p).unwrap(); + assert_eq!(u16::from(recs[0].flags()), 0x800); + assert_eq!(cigar_string(&recs[0]), "60M60S"); + assert_eq!(recs[0].sequence().len(), 120); + } +} diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index 14eb9a23..f3f64c52 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -412,8 +412,20 @@ fn comp(b: u8) -> u8 { /// `Transcript::alignScore` (`Transcript_alignScore.cpp`): re-score a /// transcript from its exons. fn align_score(tr: &WinTr, read: &ChimRead, genome: &Genome, scorer: &AlignmentScorer) -> i32 { + align_score_mm(tr, read, genome, scorer).0 +} + +/// `alignScore` with the mismatch count it recomputes alongside (`nMM`), which +/// is what a chimeric BAM record reports as `nM`. +fn align_score_mm( + tr: &WinTr, + read: &ChimRead, + genome: &Genome, + scorer: &AlignmentScorer, +) -> (i32, u32) { let r: &[u8] = if tr.str_ == 0 { read.fwd } else { &read.rev }; let mut score = 0i32; + let mut n_mm = 0u32; for e in &tr.exons { for ii in 0..e.l { let r1 = r.get(e.r + ii).copied().unwrap_or(5); @@ -423,6 +435,7 @@ fn align_score(tr: &WinTr, read: &ChimRead, genome: &Genome, scorer: &AlignmentS score += 1; } else { score -= 1; + n_mm += 1; } } } @@ -451,7 +464,7 @@ fn align_score(tr: &WinTr, read: &ChimRead, genome: &Genome, scorer: &AlignmentS let span = tr.g_length().max(1) as f64; score += (span.log2() * scorer.score_genomic_length_log2_scale - 0.5).ceil() as i32; } - score + (score, n_mm) } /// `Transcript::generateCigarP` / `ReadAlign::outputTranscriptCIGARp`: the @@ -1153,12 +1166,16 @@ fn mate_segment(tr: &WinTr, junction_exon: usize, read: &ChimRead, score: i32) - } } +#[allow(clippy::too_many_arguments)] fn to_alignment( p: &Placed, score: i32, read: &ChimRead, genome: &Genome, scorer: &AlignmentScorer, + i_tr: usize, + chim_n: usize, + is_best: bool, ) -> ChimericAlignment { let e0 = if p.t0.str_ == 1 { 0 @@ -1195,9 +1212,41 @@ fn to_alignment( acceptor_start: p.t1.exons[0].g, acceptor_cigar: cigar_p(&p.t1, read), }); + // What `chimericBAMoutput` needs. Both detectors re-score the placed + // transcripts with `alignScore` before writing BAM, which is where `AS` and + // `nM` come from (`chimericDetectionOldOutput.cpp:12-13`, + // `ChimericAlign_chimericStitching.cpp:180`). + let (s0, mm0) = align_score_mm(&p.t0, read, genome, scorer); + let (s1, mm1) = align_score_mm(&p.t1, read, genome, scorer); + chim.bam = Some(Box::new(ChimBam { + tr: [p.t0.clone(), p.t1.clone()], + max_score: [s0, s1], + n_mm: [mm0, mm1], + i_tr, + chim_n, + is_best, + read_length: read.read_length, + paired: read.paired, + })); chim } +/// A chimera as STAR's `chimericBAMoutput` receives it: the two placed +/// transcripts in read order, their re-scored `maxScore` / `nMM`, and the +/// run-level `iTr` / `chimN` / `isBestChimAlign`. +#[derive(Debug, Clone)] +pub struct ChimBam { + pub tr: [WinTr; 2], + pub max_score: [i32; 2], + pub n_mm: [u32; 2], + pub i_tr: usize, + pub chim_n: usize, + pub is_best: bool, + /// `readLength[0..2]` of the (clipped) mates the exons index. + pub read_length: [usize; 2], + pub paired: bool, +} + /// `ReadAlign::chimericDetection`: run whichever STAR detector the parameters /// select over a read's window transcripts. pub fn chimeric_detection( @@ -1219,7 +1268,8 @@ pub fn chimeric_detection( // single-best bookkeeping use the re-scored segments. let score = align_score(&p.t0, read, genome, scorer) + align_score(&p.t1, read, genome, scorer); - return vec![to_alignment(&p, score, read, genome, scorer)]; + // `chimericBAMoutput(&trChim[0], &trChim[1], this, 0, 1, true, P)`. + return vec![to_alignment(&p, score, read, genome, scorer, 0, 1, true)]; } let Some(ib) = w.best() else { return Vec::new(); @@ -1233,21 +1283,32 @@ pub fn chimeric_detection( if best == 0 { return Vec::new(); } - let kept: Vec<&MultChim> = chims.iter().filter(|c| c.score >= min_score).collect(); + let kept: Vec<(usize, &MultChim)> = chims + .iter() + .enumerate() + .filter(|(_, c)| c.score >= min_score) + .collect(); if kept.len() > params.chim_multimap_nmax { return Vec::new(); } let chim_n = kept.len(); + // `bestChimAlign`: the chimera that last raised `chimScoreBest`, i.e. the + // first one to reach the final best score. + let best_idx = chims.iter().position(|c| c.score == best); kept.into_iter() - .map(|c| { - to_alignment(&c.p, c.score, read, genome, scorer).with_multimap(MultimapInfo { - chim_n, - max_possible_score: frag_len, - max_non_chim_score: max_non_chim, - chim_score: c.score, - best_chim_score: best, - pe_merged: false, - }) + .enumerate() + .map(|(i_tr, (idx, c))| { + let is_best = best_idx == Some(idx); + to_alignment(&c.p, c.score, read, genome, scorer, i_tr, chim_n, is_best).with_multimap( + MultimapInfo { + chim_n, + max_possible_score: frag_len, + max_non_chim_score: max_non_chim, + chim_score: c.score, + best_chim_score: best, + pe_merged: false, + }, + ) }) .collect() } diff --git a/src/chimeric/mod.rs b/src/chimeric/mod.rs index 74c404e0..ef4fb4ca 100644 --- a/src/chimeric/mod.rs +++ b/src/chimeric/mod.rs @@ -9,12 +9,14 @@ // `chimericDetectionMult` under `--chimMultimapNmax`, both over the read's // window transcripts — combined two-mate transcripts for a pair. +mod bam; mod detect; mod output; mod segment; -pub use detect::{ChimRead, WinTr, chimeric_detection}; -pub use output::{ChimericJunctionWriter, build_within_bam_records}; +pub use bam::{ChimReadInput, build_chimeric_bam_records}; +pub use detect::{ChimBam, ChimRead, WinTr, chimeric_detection}; +pub use output::ChimericJunctionWriter; pub use segment::{ChimericAlignment, ChimericSegment, ExonSpan, JunctionLine, MultimapInfo}; #[cfg(test)] diff --git a/src/chimeric/output.rs b/src/chimeric/output.rs index 9610ef06..36b769d4 100644 --- a/src/chimeric/output.rs +++ b/src/chimeric/output.rs @@ -1,13 +1,7 @@ -// Chimeric.out.junction file writer and WithinBAM record builder +// Chimeric.out.junction file writer (WithinBAM records are built in `bam.rs`) -use crate::chimeric::segment::{ChimericAlignment, ChimericSegment}; +use crate::chimeric::segment::ChimericAlignment; use crate::error::Error; -use crate::genome::Genome; -use bstr::BString; -use noodles::sam; -use noodles::sam::alignment::record::MappingQuality; -use noodles::sam::alignment::record_buf::data::field::Value; -use noodles::sam::alignment::record_buf::{QualityScores, RecordBuf, Sequence}; use std::fs::File; use std::io::{BufWriter, Write}; use std::path::PathBuf; @@ -188,199 +182,6 @@ impl ChimericJunctionWriter { } } -/// Build two SAM records for `--chimOutType WithinBAM`. -/// -/// Returns `[donor_record, acceptor_record]`: -/// - Donor: normal FLAGS; full read sequence; SA tag pointing to acceptor. -/// - Acceptor: FLAG 0x0800 (supplementary); empty SEQ/QUAL; SA tag pointing to donor. -pub fn build_within_bam_records( - alignment: &ChimericAlignment, - genome: &Genome, - mapq: u8, -) -> Result, Error> { - let donor = &alignment.donor; - let acceptor = &alignment.acceptor; - - // Which end of the supplementary segment's CIGAR covers the OTHER segment's - // bases. STAR hard-clips exactly that end (`chimericBAMoutput.cpp:55` - // picks -11 for a left junction, -12 for a right one) and drops those bases - // from SEQ (`ReadAlign_alignBAM.cpp:503-511`), which is the ordinary SAM - // rule: hard-clipped bases are absent from SEQ, soft-clipped ones present. - // - // In original-read order the other segment is 5' of this one when this one - // trails. The CIGAR is in reference orientation, so a reverse segment has - // that end on the opposite side. - let read_len = alignment.read_seq.len(); - let ro = |seg: &ChimericSegment| { - if seg.is_reverse { - read_len.saturating_sub(seg.read_end) - } else { - seg.read_start - } - }; - let acceptor_trails = ro(acceptor) > ro(donor); - let hard_leading = acceptor_trails != acceptor.is_reverse; - - let donor_sa = format_sa_entry(donor, &genome.chr_name, &genome.chr_start, mapq); - let acceptor_sa = format_sa_entry(acceptor, &genome.chr_name, &genome.chr_start, mapq); - - let donor_record = build_segment_record( - &alignment.read_name, - &alignment.read_seq, - donor, - genome, - mapq, - None, - &acceptor_sa, - )?; - let acceptor_record = build_segment_record( - &alignment.read_name, - &alignment.read_seq, - acceptor, - genome, - mapq, - Some(hard_leading), - &donor_sa, - )?; - - Ok(vec![donor_record, acceptor_record]) -} - -/// Format one SA tag entry: `chr,pos,strand,CIGAR,mapQ,NM;` -fn format_sa_entry( - seg: &ChimericSegment, - chr_names: &[String], - chr_starts: &[u64], - mapq: u8, -) -> String { - let chr = &chr_names[seg.chr_idx]; - let chr_start = chr_starts[seg.chr_idx]; - let pos = seg.genome_start - chr_start + 1; // 1-based per-chr - let strand = if seg.is_reverse { '-' } else { '+' }; - let cigar = seg.cigar_string(); - format!( - "{},{},{},{},{},{};", - chr, pos, strand, cigar, mapq, seg.n_mismatch - ) -} - -/// Build one SAM record for a chimeric segment. -fn build_segment_record( - read_name: &str, - read_seq: &[u8], - seg: &ChimericSegment, - genome: &Genome, - mapq: u8, - // `None` for the representative record; `Some(hard_leading)` for the - // supplementary one, naming which CIGAR end to hard-clip. - supplementary: Option, - sa_tag: &str, -) -> Result { - use crate::io::fastq::{complement_base, decode_base}; - use noodles::sam::alignment::record::data::field::Tag; - - let mut record = RecordBuf::default(); - record.name_mut().replace(read_name.into()); - - let mut flags = sam::alignment::record::Flags::empty(); - if seg.is_reverse { - flags |= sam::alignment::record::Flags::REVERSE_COMPLEMENTED; - } - if supplementary.is_some() { - flags |= sam::alignment::record::Flags::SUPPLEMENTARY; - } - *record.flags_mut() = flags; - - *record.reference_sequence_id_mut() = Some(seg.chr_idx); - - let chr_start = genome.chr_start[seg.chr_idx]; - let pos = (seg.genome_start - chr_start + 1) as usize; - *record.alignment_start_mut() = Some( - pos.try_into() - .map_err(|e| Error::Chimeric(format!("invalid chimeric position {pos}: {e}")))?, - ); - - *record.mapping_quality_mut() = MappingQuality::new(mapq); - - // SEQ in reference orientation, matching the CIGAR. - let mut seq_bytes: Vec = if seg.is_reverse { - read_seq - .iter() - .rev() - .map(|&b| decode_base(complement_base(b))) - .collect() - } else { - read_seq.iter().map(|&b| decode_base(b)).collect() - }; - - let mut cigar_ops = seg.cigar.clone(); - - // Pad the CIGAR out to the full read with soft clips where the segment does - // not already carry them. STAR always emits these (`trimL1`/`trimR1`, - // `ReadAlign_alignBAM.cpp:235,269`), but our segments sometimes arrive as a - // bare match block, and a record whose CIGAR claims fewer query bases than - // SEQ holds cannot be written. Read coordinates are in original-read order, - // so a reverse segment's leading clip is the one past its 3' end. - { - use noodles::sam::alignment::record::cigar::{Op, op::Kind}; - let q: usize = cigar_ops - .iter() - .filter(|op| op.kind().consumes_read()) - .map(|op| op.len()) - .sum(); - if q < read_seq.len() { - let before = seg.read_start; - let after = read_seq.len().saturating_sub(seg.read_end); - let (lead, trail) = if seg.is_reverse { - (after, before) - } else { - (before, after) - }; - let lead_present = matches!(cigar_ops.first(), Some(o) if o.kind() == Kind::SoftClip); - let trail_present = matches!(cigar_ops.last(), Some(o) if o.kind() == Kind::SoftClip); - if trail > 0 && !trail_present { - cigar_ops.push(Op::new(Kind::SoftClip, trail)); - } - if lead > 0 && !lead_present { - cigar_ops.insert(0, Op::new(Kind::SoftClip, lead)); - } - } - } - - // Supplementary: hard-clip the end covering the other segment and drop - // those bases from SEQ. Leaving them soft-clipped with an empty SEQ, as - // this did before, is not representable -- soft clips consume the query, so - // a 55S65M CIGAR asserts a 120-base SEQ and the record fails to write. - if let Some(hard_leading) = supplementary { - use noodles::sam::alignment::record::cigar::{Op, op::Kind}; - let idx = if hard_leading { 0 } else { cigar_ops.len() - 1 }; - if let Some(op) = cigar_ops.get(idx).copied() - && op.kind() == Kind::SoftClip - { - let n = op.len(); - cigar_ops[idx] = Op::new(Kind::HardClip, n); - if hard_leading { - seq_bytes.drain(..n.min(seq_bytes.len())); - } else { - let keep = seq_bytes.len().saturating_sub(n); - seq_bytes.truncate(keep); - } - } - } - - *record.cigar_mut() = cigar_ops.iter().copied().collect(); - *record.sequence_mut() = Sequence::from(seq_bytes); - // QUAL is not carried for chimeric segments. - *record.quality_scores_mut() = QualityScores::default(); - - let data = record.data_mut(); - data.insert(Tag::new(b'S', b'A'), Value::String(BString::from(sa_tag))); - data.insert(Tag::new(b'N', b'M'), Value::from(seg.n_mismatch as i32)); - data.insert(Tag::ALIGNMENT_SCORE, Value::from(seg.score)); - - Ok(record) -} - #[cfg(test)] mod tests { use super::*; @@ -766,407 +567,4 @@ mod tests { assert_eq!(fields[10], "501"); // donor start assert_eq!(fields[12], "1501"); // acceptor start } - - // --- build_within_bam_records tests --- - - fn make_genome_2chr() -> crate::genome::Genome { - use crate::genome::Genome; - Genome { - transform_blocks: None, - sequence: vec![0u8; 2048].into(), - n_genome: 1024, - n_genome_real: 1024, - n_chr_real: 2, - chr_name: vec!["chr9".to_string(), "chr22".to_string()], - chr_length: vec![512, 512], - chr_start: vec![0, 512, 1024], - } - } - - #[test] - fn test_within_bam_returns_two_records() { - use cigar::op::{Kind, Op}; - let donor = ChimericSegment { - chr_idx: 0, - genome_start: 100, - genome_end: 163, - is_reverse: false, - read_start: 0, - read_end: 63, - // Real segments span the whole read: the part belonging to the - // other segment is soft-clipped, not absent. - cigar: vec![Op::new(Kind::Match, 63)], - score: 63, - n_mismatch: 0, - first_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - last_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - }; - let acceptor = ChimericSegment { - chr_idx: 1, - genome_start: 600, - genome_end: 637, - is_reverse: false, - read_start: 63, - read_end: 100, - cigar: vec![Op::new(Kind::SoftClip, 63), Op::new(Kind::Match, 37)], - score: 37, - n_mismatch: 1, - first_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - last_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - }; - let alignment = ChimericAlignment::new( - donor, - acceptor, - 0, - 0, - 0, - vec![0u8; 100], - "READ_001".to_string(), - ); - let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); - - assert_eq!(records.len(), 2); - } - - #[test] - fn test_within_bam_donor_not_supplementary() { - use cigar::op::{Kind, Op}; - let donor = ChimericSegment { - chr_idx: 0, - genome_start: 100, - genome_end: 163, - is_reverse: false, - read_start: 0, - read_end: 63, - cigar: vec![Op::new(Kind::Match, 63)], - score: 63, - n_mismatch: 0, - first_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - last_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - }; - let acceptor = ChimericSegment { - chr_idx: 1, - genome_start: 600, - genome_end: 637, - is_reverse: false, - read_start: 63, - read_end: 100, - cigar: vec![Op::new(Kind::Match, 37)], - score: 37, - n_mismatch: 1, - first_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - last_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - }; - let alignment = ChimericAlignment::new( - donor, - acceptor, - 0, - 0, - 0, - vec![0u8; 100], - "READ_001".to_string(), - ); - let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); - - let donor_flags = records[0].flags(); - let acceptor_flags = records[1].flags(); - - assert!( - !donor_flags.is_supplementary(), - "donor must not be supplementary" - ); - assert!( - acceptor_flags.is_supplementary(), - "acceptor must be supplementary (0x800)" - ); - } - - #[test] - fn test_within_bam_sa_tag_format() { - use cigar::op::{Kind, Op}; - use noodles::sam::alignment::record::data::field::Tag; - let donor = ChimericSegment { - chr_idx: 0, - genome_start: 100, - genome_end: 163, - is_reverse: false, - read_start: 0, - read_end: 63, - cigar: vec![Op::new(Kind::Match, 63)], - score: 63, - n_mismatch: 2, - first_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - last_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - }; - let acceptor = ChimericSegment { - chr_idx: 1, - genome_start: 600, - genome_end: 637, - is_reverse: true, - read_start: 63, - read_end: 100, - cigar: vec![Op::new(Kind::Match, 37)], - score: 37, - n_mismatch: 1, - first_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - last_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - }; - let alignment = ChimericAlignment::new( - donor, - acceptor, - 0, - 0, - 0, - vec![0u8; 100], - "READ_001".to_string(), - ); - let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); - - // Donor record's SA tag should point to acceptor - let sa_tag = Tag::new(b'S', b'A'); - let donor_sa = records[0].data().get(&sa_tag).unwrap(); - let donor_sa_str = format!("{donor_sa:?}"); - // SA tag: chr22,89,-,37M,255,1; (pos = 600-512+1=89, strand=-, nm=1) - assert!( - donor_sa_str.contains("chr22"), - "SA tag must name acceptor chr" - ); - assert!( - donor_sa_str.contains("89"), - "SA tag must have per-chr position" - ); - assert!( - donor_sa_str.contains('-'), - "SA tag must reflect reverse strand" - ); - - // Acceptor record's SA tag should point to donor - let acceptor_sa = records[1].data().get(&sa_tag).unwrap(); - let acceptor_sa_str = format!("{acceptor_sa:?}"); - assert!( - acceptor_sa_str.contains("chr9"), - "SA tag must name donor chr" - ); - } - - /// A segment whose CIGAR is a bare match block must still produce a record - /// whose CIGAR spans the read. - /// - /// Detection sometimes hands us a segment with no clip ops at all (a bare - /// `59M` on a 120-base read). STAR always emits the surrounding soft clips - /// (`trimL1`/`trimR1`, `ReadAlign_alignBAM.cpp:235,269`); without them the - /// record claims fewer query bases than SEQ carries and BAM writing fails - /// with "read length-sequence length mismatch". - #[test] - fn test_within_bam_pads_a_bare_cigar_to_span_the_read() { - use cigar::op::{Kind, Op}; - let bare = |chr_idx: usize, gs: u64, rs: usize, re: usize, rev: bool| ChimericSegment { - chr_idx, - genome_start: gs, - genome_end: gs + (re - rs) as u64, - is_reverse: rev, - read_start: rs, - read_end: re, - // No clips, deliberately. - cigar: vec![Op::new(Kind::Match, re - rs)], - score: (re - rs) as i32, - n_mismatch: 0, - first_exon: ExonSpan { - genome_start: gs, - genome_end: gs + (re - rs) as u64, - read_start: rs, - read_end: re, - }, - last_exon: ExonSpan { - genome_start: gs, - genome_end: gs + (re - rs) as u64, - read_start: rs, - read_end: re, - }, - }; - // Reverse acceptor, mirroring the fixture read that first tripped this. - let alignment = ChimericAlignment::new( - bare(0, 100, 0, 61, false), - bare(1, 600, 61, 120, true), - 0, - 0, - 0, - vec![0u8; 120], - "READ_BARE".to_string(), - ); - let records = build_within_bam_records(&alignment, &make_genome_2chr(), 255).unwrap(); - - for rec in &records { - let q: usize = rec - .cigar() - .as_ref() - .iter() - .filter(|op| op.kind().consumes_read()) - .map(|op| op.len()) - .sum(); - assert_eq!( - q, - rec.sequence().len(), - "CIGAR query length must match SEQ: {}", - crate::align::transcript::cigar_to_string(rec.cigar().as_ref()) - ); - } - } - - #[test] - fn test_within_bam_supplementary_is_hard_clipped_with_trimmed_seq() { - use cigar::op::{Kind, Op}; - let donor = ChimericSegment { - chr_idx: 0, - genome_start: 100, - genome_end: 163, - is_reverse: false, - read_start: 0, - read_end: 63, - // Real segments span the whole read: the part belonging to the - // other segment is soft-clipped, not absent. - cigar: vec![Op::new(Kind::Match, 63), Op::new(Kind::SoftClip, 37)], - score: 63, - n_mismatch: 0, - first_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - last_exon: ExonSpan { - genome_start: 100, - genome_end: 163, - read_start: 0, - read_end: 63, - }, - }; - let acceptor = ChimericSegment { - chr_idx: 1, - genome_start: 600, - genome_end: 637, - is_reverse: false, - read_start: 63, - read_end: 100, - cigar: vec![Op::new(Kind::SoftClip, 63), Op::new(Kind::Match, 37)], - score: 37, - n_mismatch: 0, - first_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - last_exon: ExonSpan { - genome_start: 600, - genome_end: 637, - read_start: 63, - read_end: 100, - }, - }; - let read_seq = vec![0u8; 100]; // 100 A bases - let alignment = - ChimericAlignment::new(donor, acceptor, 0, 0, 0, read_seq, "READ_001".to_string()); - let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); - - // The representative record carries the whole read against a CIGAR that - // spans it. - assert_eq!(records[0].sequence().len(), 100, "donor keeps the full SEQ"); - assert_eq!( - crate::align::transcript::cigar_to_string(records[0].cigar().as_ref()), - "63M37S" - ); - - // The supplementary hard-clips the end covering the donor and drops - // those bases from SEQ, as STAR does (chimericBAMoutput.cpp:55, - // alignBAM.cpp:503-511). Leaving them soft-clipped with an empty SEQ is - // not representable: soft clips consume the query, so the record would - // claim 100 bases and carry none. - assert_eq!( - crate::align::transcript::cigar_to_string(records[1].cigar().as_ref()), - "63H37M", - "supplementary hard-clips the other segment's bases" - ); - assert_eq!( - records[1].sequence().len(), - 37, - "SEQ excludes hard-clipped bases" - ); - - // The invariant the SAM writer enforces: CIGAR query length == SEQ length. - for rec in &records { - let q: usize = rec - .cigar() - .as_ref() - .iter() - .filter(|op| op.kind().consumes_read()) - .map(|op| op.len()) - .sum(); - assert_eq!(q, rec.sequence().len(), "CIGAR query length must match SEQ"); - } - } } diff --git a/src/chimeric/segment.rs b/src/chimeric/segment.rs index 0c6c7dea..e9f85e85 100644 --- a/src/chimeric/segment.rs +++ b/src/chimeric/segment.rs @@ -143,6 +143,8 @@ pub struct ChimericAlignment { pub multimap: Option, /// The junction line, set by the STAR detector; see [`JunctionLine`]. pub junction_line: Option, + /// What `--chimOutType WithinBAM` writes from; see [`crate::chimeric::ChimBam`]. + pub bam: Option>, } impl ChimericAlignment { @@ -168,6 +170,7 @@ impl ChimericAlignment { read_name, multimap: None, junction_line: None, + bam: None, } } diff --git a/src/io/sam.rs b/src/io/sam.rs index 13ebc9f7..288b968d 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -1090,7 +1090,7 @@ where /// Insert `RG:Z:` on the record when an ID is set. `Parameters::try_parse_from` /// auto-ORs `SamAttributes::RG` into `out_sam_attributes` whenever an RG line is /// configured, so `rg_id.is_some()` implies the attribute is wanted. -fn maybe_insert_rg_tag(record: &mut RecordBuf, rg_id: Option<&str>) { +pub(crate) fn maybe_insert_rg_tag(record: &mut RecordBuf, rg_id: Option<&str>) { if let Some(id) = rg_id { record .data_mut() @@ -1118,7 +1118,7 @@ pub fn add_gene_tags(records: &mut [RecordBuf], gx: &str, gn: &str, attrs: SamAt /// 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. -fn apply_sam_flag_or_and(record: &mut RecordBuf, params: &Parameters) { +pub(crate) fn apply_sam_flag_or_and(record: &mut RecordBuf, params: &Parameters) { let or_bits = params.out_sam_flag_or as u16; let and_bits = params.out_sam_flag_and.min(u16::MAX as u32) as u16; let bits = u16::from(record.flags()); @@ -1164,7 +1164,7 @@ fn paired_tlen( /// the BAM binary QUAL field, per SAM spec §4.2.3. /// /// `saturating_sub` clamps malformed bytes < 33 to 0 rather than underflowing. -fn fastq_qual_to_phred(qual: &[u8]) -> Vec { +pub(crate) fn fastq_qual_to_phred(qual: &[u8]) -> Vec { qual.iter().map(|&b| b.saturating_sub(33)).collect() } @@ -1687,10 +1687,12 @@ fn build_paired_mate_record( data.insert(Tag::ALIGNMENT_SCORE, Value::from(combined_score)); } // nM (mismatch count) before NM (edit distance), matching STAR's tag order. + // nM is the whole pair's `trOut.nMM`, written on both mates (as AS is); NM + // below is per mate (`samAttrNM_MD` over this mate's exons). if attrs.contains(SamAttributes::NMM) { data.insert( Tag::new(b'n', b'M'), - Value::from(transcript.n_mismatch as i32), + Value::from((transcript.n_mismatch + mate_transcript.n_mismatch) as i32), ); } if attrs.contains(SamAttributes::NM) { @@ -2425,7 +2427,7 @@ mod tests { // Check TLEN assert_eq!(rec.template_length(), 250); - // AS is the combined score (STAR behavior); nM is per-mate mismatches + // AS and nM are both the whole pair's (STAR's trOut.maxScore / trOut.nMM) let data = rec.data(); assert_eq!( data.get(&Tag::ALIGNMENT_SCORE), @@ -2434,8 +2436,8 @@ mod tests { ); assert_eq!( data.get(&Tag::new(b'n', b'M')), - Some(&Value::from(0_i32)), - "nM should be 0 (no mismatches in this mate)" + Some(&Value::from(1_i32)), + "nM should be the pair's mismatches (0 + 1), on both mates" ); } @@ -3930,11 +3932,11 @@ mod tests { Some(&Value::from(180_i32)), "Mate1 AS should be combined score (100+80=180)" ); - // NM attribute maps to nM tag (mismatches only) + // nM is the pair's mismatch count (STAR's trOut.nMM), on both mates assert_eq!( rec1.data().get(&Tag::new(b'n', b'M')), - Some(&Value::from(0_i32)), - "Mate1 nM should be 0" + Some(&Value::from(2_i32)), + "Mate1 nM should be the pair's 0 + 2" ); // Mate2 also gets combined AS @@ -3967,7 +3969,7 @@ mod tests { assert_eq!( rec2.data().get(&Tag::new(b'n', b'M')), Some(&Value::from(2_i32)), - "Mate2 nM should be 2 (mismatches only, not edit distance)" + "Mate2 nM should be the pair's 0 + 2 (mismatches only, not edit distance)" ); } diff --git a/src/lib.rs b/src/lib.rs index 1b0024f3..1e2d3b8d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1669,13 +1669,6 @@ fn align_reads_single_end( )?; } } - if params.chim_out_within_bam() { - use crate::chimeric::build_within_bam_records; - for chim_aln in &batch.chimeric_alns { - let supp = build_within_bam_records(chim_aln, &index.genome, 255)?; - writer.write_batch(&supp)?; - } - } // Write unmapped FASTQ records if let Some(ref mut uw) = unmapped_writer { @@ -1736,14 +1729,6 @@ fn align_reads_single_end( )?; } } - if params.chim_out_within_bam() { - use crate::chimeric::build_within_bam_records; - for chim_aln in &meta.chimeric_alns { - let supp = - build_within_bam_records(chim_aln, &index.genome, 255)?; - writer.write_batch(&supp)?; - } - } } else { // Skip these records in the temp file (advance reader) crate::io::sam::bysj_read_n_records( @@ -1919,6 +1904,47 @@ fn align_reads_single_end( } } + // --chimOutType WithinBAM: STAR writes the chimera in place of + // the read's alignments and returns before outputAlignments + // (`ReadAlign_oneRead.cpp:99`), so the read gets none of the + // mapped-read accounting: no mapped stats, junctions, counts, + // transcriptome or unmapped output. Unmapped reasons were + // already counted by mappedFilter, which runs first. + if params.chim_out_within_bam() && !chimeric_alns.is_empty() { + if transcripts.is_empty() { + stats.record_alignment(n_for_mapq, max_multimaps); + if let Some(reason) = unmapped_reason { + stats.record_unmapped_reason(reason); + } + } else { + stats.record_input_read(); + } + let input = crate::chimeric::ChimReadInput { + name: &out_read_name, + seq: [&read.sequence, &[]], + qual: [&read.quality, &[]], + clip: [[clip5p, clip3p], [0, 0]], + }; + for record in crate::chimeric::build_chimeric_bam_records( + &chimeric_alns, + &input, + &index.genome, + params, + )? { + buffer.push(record); + } + return Ok(AlignmentBatchResults { + sam_records: buffer, + chimeric_alns, + primary_junction_keys: Vec::new(), + transcriptome_records: Vec::new(), + unmapped_mate1: Vec::new(), + unmapped_mate2: Vec::new(), + signal_contrib: Vec::new(), + signal_n_tr: 0, + }); + } + // Record stats (atomic, lock-free) // For too-many-loci, n_for_mapq carries the true loci count // while transcripts is empty @@ -3010,13 +3036,6 @@ fn align_reads_paired_end( )?; } } - if params.chim_out_within_bam() { - use crate::chimeric::build_within_bam_records; - for chim_aln in &batch.chimeric_alns { - let supp = build_within_bam_records(chim_aln, &index.genome, 255)?; - writer.write_batch(&supp)?; - } - } if let Some(ref mut uw1) = unmapped_writer1 { for (name, seq, qual) in &batch.unmapped_mate1 { uw1.write_record(name, seq, qual)?; @@ -3080,14 +3099,6 @@ fn align_reads_paired_end( )?; } } - if params.chim_out_within_bam() { - use crate::chimeric::build_within_bam_records; - for chim_aln in &meta.chimeric_alns { - let supp = - build_within_bam_records(chim_aln, &index.genome, 255)?; - writer.write_batch(&supp)?; - } - } } else { crate::io::sam::bysj_read_n_records( &mut reader, @@ -3290,6 +3301,45 @@ fn align_reads_paired_end( let (results, pe_chimeric, n_for_mapq, unmapped_reason) = align_paired_read(&m1_seq, &m2_seq, &paired_read.name, &index, params)?; + // --chimOutType WithinBAM: the chimera replaces the read's + // alignments and STAR returns before outputAlignments + // (`ReadAlign_oneRead.cpp:99`); see the single-end path. + if params.chim_out_within_bam() && !pe_chimeric.is_empty() { + stats.record_chimeric(); + if results.is_empty() { + stats.record_alignment(n_for_mapq, max_multimaps); + stats.record_unmapped_reason( + unmapped_reason.unwrap_or(crate::stats::UnmappedReason::Other), + ); + } else { + stats.record_input_read(); + } + let input = crate::chimeric::ChimReadInput { + name: &out_read_name, + seq: [&paired_read.mate1.sequence, &paired_read.mate2.sequence], + qual: [&paired_read.mate1.quality, &paired_read.mate2.quality], + clip: [[m1_clip5p, m1_clip3p], [m2_clip5p, m2_clip3p]], + }; + for record in crate::chimeric::build_chimeric_bam_records( + &pe_chimeric, + &input, + &index.genome, + params, + )? { + buffer.push(record); + } + return Ok(AlignmentBatchResults { + sam_records: buffer, + chimeric_alns: pe_chimeric, + primary_junction_keys: Vec::new(), + transcriptome_records: Vec::new(), + unmapped_mate1: Vec::new(), + unmapped_mate2: Vec::new(), + signal_contrib: Vec::new(), + signal_n_tr: 0, + }); + } + // Classify the result for stats and SAM output let has_half_mapped = results .iter() diff --git a/src/params/mod.rs b/src/params/mod.rs index a3539e08..330e7736 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1212,7 +1212,7 @@ pub struct Parameters { pub chim_score_junction_non_gtag: i32, /// Chimeric output type - #[arg(long = "chimOutType", num_args = 1..=2, default_values_t = vec!["Junctions".to_string()])] + #[arg(long = "chimOutType", num_args = 1.., default_values_t = vec!["Junctions".to_string()])] pub chim_out_type: Vec, // ── STARsolo (single-cell) ────────────────────────────────────────── @@ -1382,6 +1382,21 @@ impl Parameters { self.chim_out_type.iter().any(|s| s == "WithinBAM") } + /// `--chimOutType ... HardClip|SoftClip`: whether the supplementary chimeric + /// segment is hard-clipped. STAR reads the values in order, so the last of + /// the two wins; hard clipping is the default (`ParametersChimeric_initialize.cpp:10-30`). + pub fn chim_out_bam_hard_clip(&self) -> bool { + self.chim_out_type + .iter() + .rev() + .find_map(|s| match s.as_str() { + "HardClip" => Some(true), + "SoftClip" => Some(false), + _ => None, + }) + .unwrap_or(true) + } + /// True if the user provided a non-default `--outSAMattrRGline`. pub fn rg_line_set(&self) -> bool { !self.out_sam_attr_rg_line.is_empty() && self.out_sam_attr_rg_line[0] != "-" @@ -1700,6 +1715,44 @@ impl Parameters { } } + // --chimOutType (`ParametersChimeric_initialize.cpp`): unknown values are + // fatal. SeparateSAMold is a STAR value this build does not write, so it + // is refused rather than silently ignored. + if params.chim_segment_min > 0 { + for t in ¶ms.chim_out_type { + match t.as_str() { + "Junctions" | "WithinBAM" | "HardClip" | "SoftClip" => {} + "SeparateSAMold" => { + return Err(command.error( + ErrorKind::InvalidValue, + "--chimOutType SeparateSAMold is not supported; use Junctions and/or WithinBAM", + )); + } + other => { + return Err(command.error( + ErrorKind::InvalidValue, + format!( + "unknown --chimOutType value '{other}'; allowed: Junctions, \ + WithinBAM, HardClip, SoftClip" + ), + )); + } + } + } + if params.chim_out_within_bam() { + // WithinBAM needs BAM output (`:76-81`), and adds the NM attribute + // that the SA tag reads (`:99-102`). + if params.out_sam_type.format != OutSamFormat::Bam { + return Err(command.error( + ErrorKind::InvalidValue, + "--chimOutType WithinBAM requires BAM output: use --outSAMtype BAM Unsorted \ + or SortedByCoordinate", + )); + } + params.out_sam_attributes |= SamAttributes::NM; + } + } + // Unknown --quantMode values are fatal, as in STAR // (`Parameters.cpp:898-936`: "unrecognized option in --quantMode"); a // leading "-" means none. @@ -2627,6 +2680,46 @@ mod tests { ); } + #[test] + fn chim_out_type_follows_star() { + let base = ["--readFilesIn", "r.fq", "--chimSegmentMin", "12"]; + let with = |extra: &[&str]| { + let mut a = base.to_vec(); + a.extend_from_slice(extra); + try_parse(&a) + }; + // WithinBAM needs BAM output, and adds NM. + assert!(with(&["--chimOutType", "WithinBAM"]).is_err()); + let p = with(&[ + "--chimOutType", + "WithinBAM", + "--outSAMtype", + "BAM", + "Unsorted", + ]) + .unwrap(); + assert!(p.out_sam_attributes.contains(SamAttributes::NM)); + assert!(p.chim_out_bam_hard_clip()); + // Three values, and the last of HardClip / SoftClip wins. + let p = with(&[ + "--chimOutType", + "Junctions", + "WithinBAM", + "SoftClip", + "--outSAMtype", + "BAM", + "Unsorted", + ]) + .unwrap(); + assert!(p.chim_out_junctions() && p.chim_out_within_bam()); + assert!(!p.chim_out_bam_hard_clip()); + // Unknown values, and the unimplemented SeparateSAMold, are errors. + assert!(with(&["--chimOutType", "Junction"]).is_err()); + assert!(with(&["--chimOutType", "SeparateSAMold"]).is_err()); + // Without chimeric detection the value is not checked, as in STAR. + assert!(try_parse(&["--readFilesIn", "r.fq", "--chimOutType", "WithinBAM"]).is_ok()); + } + #[test] fn chimeric_params() { let p = try_parse(&[ @@ -2639,6 +2732,9 @@ mod tests { "--chimOutType", "WithinBAM", "SoftClip", + "--outSAMtype", + "BAM", + "Unsorted", ]) .unwrap(); assert_eq!(p.chim_segment_min, 20); diff --git a/src/stats.rs b/src/stats.rs index 636af48a..75fc9463 100644 --- a/src/stats.rs +++ b/src/stats.rs @@ -136,6 +136,15 @@ impl AlignmentStats { } } + /// Count an input read without placing it in any mapping category. + /// + /// For a read whose chimera is written into the BAM: STAR skips + /// `outputAlignments` for it (`ReadAlign_oneRead.cpp:99`), which is where + /// the uniquely/multi-mapped counts are kept, but still counts it as input. + pub fn record_input_read(&self) { + self.total_reads.fetch_add(1, Ordering::Relaxed); + } + /// Record the length of an input read (for average input read length) pub fn record_read_bases(&self, len: u64) { self.read_bases.fetch_add(len, Ordering::Relaxed);