From d6af096ebf0ebbc58c53b1e11b652c2c9c5fd194 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 29 Sep 2026 22:19:38 +0200 Subject: [PATCH] fix(chimeric): well-formed single-end WithinBAM records, STAR hard-clip rules Soft-clip (Tier 1b) and residual (Tier 3) segments kept sub-sequence CIGARs, so a donor could carry 40M with a 100 bp SEQ and crash BAM output. CIGARs are now lifted to full-read orientation and padded, SE segments are ordered by read position, the representative follows STAR's chimRepresent, and the supplementary record is hard-clipped on its junction side (bamHardClip, SoftClip keeps soft clips). Also fixes Tier 1b clip side for reverse primaries and strand-aware segment ordering in Tiers 1b, 2 and 3. Closes #279 Co-Authored-By: Claude Opus 5.5 --- src/chimeric/detect.rs | 138 +++++++++++---- src/chimeric/output.rs | 327 ++++++++++++++++++++++++++++++------ src/lib.rs | 34 +++- src/params/mod.rs | 6 + tests/alignment_features.rs | 124 ++++++++++++++ 5 files changed, 537 insertions(+), 92 deletions(-) diff --git a/src/chimeric/detect.rs b/src/chimeric/detect.rs index a32af77a..9840b47e 100644 --- a/src/chimeric/detect.rs +++ b/src/chimeric/detect.rs @@ -43,7 +43,16 @@ impl<'a> ChimericDetector<'a> { } let read_len = read_seq.len(); - let [left_clip, right_clip] = transcript.count_soft_clips(); + // Soft clips in read (5'->3') orientation: a reverse-strand CIGAR is laid out on the + // reverse-complemented read, so its left clip is the read's 3' end. + let [left_clip, right_clip] = { + let [l, r] = transcript.count_soft_clips(); + if transcript.is_reverse { + [r, l] + } else { + [l, r] + } + }; let score_min = params.chim_score_min; let score_drop_max = params.chim_score_drop_max; let non_gtag_penalty = params.chim_score_junction_non_gtag; @@ -98,16 +107,12 @@ impl<'a> ChimericDetector<'a> { continue; } - // Shift sub-seq read coords into full-read space for right clips - let clip_tr = if is_right { - adjust_read_positions(clip_tr_raw, clip_start) - } else { - clip_tr_raw - }; + // Lift sub-seq read coords and CIGAR into full-read space + let clip_tr = lift_clip_transcript(clip_tr_raw, clip_start, clip_len, read_len); - // Determine donor / acceptor by read order - let primary_rs = transcript.exons[0].read_start; - let clip_rs = clip_tr.exons[0].read_start; + // Determine donor / acceptor by read order (5'->3' of the original read) + let primary_rs = ro_coords(transcript, read_len).0; + let clip_rs = ro_coords(&clip_tr, read_len).0; let (tr_donor, tr_acceptor): (&Transcript, &Transcript) = if primary_rs <= clip_rs { (transcript, &clip_tr) } else { @@ -221,17 +226,19 @@ impl<'a> ChimericDetector<'a> { let intron_max = params.align_intron_max as u64; let overhang_min = params.chim_junction_overhang_min as usize; - // Outer boundaries of the existing chimeric pair in read space - let left_covered = chim.donor.read_start.min(chim.acceptor.read_start); - let right_covered = chim.donor.read_end.max(chim.acceptor.read_end); + // Outer boundaries of the existing chimeric pair in read (5'->3') space + let (donor_ro_start, donor_ro_end) = segment_ro_span(&chim.donor, read_len); + let (acceptor_ro_start, acceptor_ro_end) = segment_ro_span(&chim.acceptor, read_len); + let left_covered = donor_ro_start.min(acceptor_ro_start); + let right_covered = donor_ro_end.max(acceptor_ro_end).min(read_len); // Which segment is at the left / right boundary - let left_partner = if chim.donor.read_start <= chim.acceptor.read_start { + let left_partner = if donor_ro_start <= acceptor_ro_start { &chim.donor } else { &chim.acceptor }; - let right_partner = if chim.donor.read_end >= chim.acceptor.read_end { + let right_partner = if donor_ro_end >= acceptor_ro_end { &chim.donor } else { &chim.acceptor @@ -246,7 +253,7 @@ impl<'a> ChimericDetector<'a> { let mut results = Vec::new(); for (clip_start, clip_end, partner_seg) in candidates { - let clip_len = clip_end - clip_start; + let clip_len = clip_end.saturating_sub(clip_start); if clip_len < min_seg { continue; } @@ -284,19 +291,16 @@ impl<'a> ChimericDetector<'a> { continue; } - // Shift sub-seq read coords into full-read space - let clip_tr = if clip_start > 0 { - adjust_read_positions(clip_tr_raw, clip_start) - } else { - clip_tr_raw - }; + // Lift sub-seq read coords and CIGAR into full-read space + let clip_tr = lift_clip_transcript(clip_tr_raw, clip_start, clip_len, read_len); let new_seg = transcript_to_segment(&clip_tr) .map_err(|e| Error::Chimeric(format!("tier3 segment: {e}")))?; // Donor / acceptor ordered by read position let (donor_seg, acceptor_seg): (&ChimericSegment, &ChimericSegment) = - if new_seg.read_start <= partner_seg.read_start { + if segment_ro_span(&new_seg, read_len).0 <= segment_ro_span(partner_seg, read_len).0 + { (&new_seg, partner_seg) } else { (partner_seg, &new_seg) @@ -441,8 +445,9 @@ impl<'a> ChimericDetector<'a> { return Ok(None); } - // Determine donor/acceptor based on read position - let (donor_t, acceptor_t) = if t1.exons[0].read_start < t2.exons[0].read_start { + // Determine donor/acceptor based on read position (5'->3' of the original read) + let read_len = read_seq.len(); + let (donor_t, acceptor_t) = if ro_coords(t1, read_len).0 < ro_coords(t2, read_len).0 { (t1, t2) } else { (t2, t1) @@ -611,18 +616,62 @@ pub fn detect_inter_mate_chimeric( Some(chim) } -/// Shift all exon read_start/read_end values in a transcript by `offset`. +/// Lift a transcript stitched against a sub-sequence of the read into full-read space. /// -/// Used when a transcript was stitched against a sub-slice of the read (e.g. a right soft-clip -/// at position `offset`) so that its read coordinates become relative to the full read. -fn adjust_read_positions(mut tr: Transcript, offset: usize) -> Transcript { +/// `clip_start`/`clip_len` locate the sub-sequence on the original (forward) read. Read +/// coordinates of a transcript follow its CIGAR orientation (reverse-strand transcripts are +/// laid out on the reverse-complemented read), so the sub-sequence starts at `clip_start` for a +/// forward transcript and at `read_len - clip_start - clip_len` for a reverse one. The CIGAR is +/// padded with soft clips so that it spans the whole read, like transcripts from the main path. +fn lift_clip_transcript( + mut tr: Transcript, + clip_start: usize, + clip_len: usize, + read_len: usize, +) -> Transcript { + use noodles::sam::alignment::record::cigar::op::{Kind, Op}; + + let after = read_len.saturating_sub(clip_start + clip_len); + let (left_pad, right_pad) = if tr.is_reverse { + (after, clip_start) + } else { + (clip_start, after) + }; for exon in &mut tr.exons { - exon.read_start += offset; - exon.read_end += offset; + exon.read_start += left_pad; + exon.read_end += left_pad; + } + if left_pad > 0 { + match tr.cigar.first_mut() { + Some(op) if op.kind() == Kind::SoftClip => { + *op = Op::new(Kind::SoftClip, op.len() + left_pad); + } + _ => tr.cigar.insert(0, Op::new(Kind::SoftClip, left_pad)), + } + } + if right_pad > 0 { + match tr.cigar.last_mut() { + Some(op) if op.kind() == Kind::SoftClip => { + *op = Op::new(Kind::SoftClip, op.len() + right_pad); + } + _ => tr.cigar.push(Op::new(Kind::SoftClip, right_pad)), + } } tr } +/// Read-orientation (5'->3' of the original read) half-open span `[start, end)` of a segment. +fn segment_ro_span(seg: &ChimericSegment, read_len: usize) -> (usize, usize) { + if seg.is_reverse { + ( + read_len.saturating_sub(seg.read_end), + read_len.saturating_sub(seg.read_start), + ) + } else { + (seg.read_start, seg.read_end) + } +} + /// Compute read-orientation (5'→3' of original read) start/end for a transcript. /// /// STAR uses "ro" coords so that clipping amounts are always measured from the 5' end of the @@ -1489,4 +1538,31 @@ mod tests { assert_eq!(with.len(), 1, "diffMates should waive the inter-mate gap"); assert_ne!(with[0].donor.chr_idx, with[0].acceptor.chr_idx); } + + /// Sub-sequence transcripts (soft-clip / residual re-seeding) are lifted into + /// full-read CIGAR orientation with a CIGAR spanning the whole read (#279). + #[test] + fn test_lift_clip_transcript_forward_and_reverse() { + use crate::align::transcript::cigar_to_string; + + // Forward: sub-sequence read[60..100) of a 100 bp read. + let fwd = lift_clip_transcript(make_transcript(0, 10, 50, false), 60, 40, 100); + assert_eq!(fwd.exons[0].read_start, 60); + assert_eq!(fwd.exons[0].read_end, 100); + assert_eq!(cigar_to_string(&fwd.cigar), "60S40M"); + assert_eq!(ro_coords(&fwd, 100), (60, 99)); + + // Reverse: same sub-sequence; on the reverse-complemented read it is the + // first 40 bases, so the padding goes to the right. + let rev = lift_clip_transcript(make_transcript(0, 10, 50, true), 60, 40, 100); + assert_eq!(rev.exons[0].read_start, 0); + assert_eq!(rev.exons[0].read_end, 40); + assert_eq!(cigar_to_string(&rev.cigar), "40M60S"); + assert_eq!(ro_coords(&rev, 100), (60, 99)); + + // Left sub-sequence read[0..30) on the reverse strand. + let rev_left = lift_clip_transcript(make_transcript(0, 10, 40, true), 0, 30, 100); + assert_eq!(cigar_to_string(&rev_left.cigar), "70S30M"); + assert_eq!(ro_coords(&rev_left, 100), (0, 29)); + } } diff --git a/src/chimeric/output.rs b/src/chimeric/output.rs index 627c5a8a..fdfd5b9b 100644 --- a/src/chimeric/output.rs +++ b/src/chimeric/output.rs @@ -1,11 +1,12 @@ // Chimeric.out.junction file writer and WithinBAM record builder +use crate::align::transcript::cigar_to_string; use crate::chimeric::segment::{ChimericAlignment, ChimericSegment}; 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::{MappingQuality, cigar}; use noodles::sam::alignment::record_buf::data::field::Value; use noodles::sam::alignment::record_buf::{QualityScores, RecordBuf, Sequence}; use std::fs::File; @@ -127,56 +128,167 @@ impl ChimericJunctionWriter { } } -/// Build two SAM records for `--chimOutType WithinBAM`. +/// Build the 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. +/// Mirrors STAR's `ChimericAlign::chimericBAMoutput` + `ReadAlign::alignBAM`. +/// Returns two records: `[donor, acceptor]` for paired-end, and for +/// single-end the two segments in read order (STAR's `trChim[0]`, +/// `trChim[1]`): +/// +/// - Every CIGAR accounts for the whole read (the part outside the segment is +/// clipped), so CIGAR query length always matches SEQ. +/// - Single-end (`chimType==3`): the segment with the higher score is the +/// representative record (donor only if its score is strictly higher); the +/// other one is supplementary (FLAG 0x800). With `hard_clip` (STAR's +/// default `HardClip`) the supplementary record is hard-clipped on its +/// chimeric-junction side (alignType -11/-12) and its SEQ drops those +/// bases; with `SoftClip` (alignType -13) it keeps soft clips and the full +/// SEQ. +/// - Paired-end: donor is representative, acceptor is supplementary with an +/// empty SEQ (unchanged legacy behavior). +/// - Each of the two records carries an SA tag describing the other one +/// (`chr,pos,strand,CIGAR,MAPQ,NM;`), using that record's final CIGAR. pub fn build_within_bam_records( alignment: &ChimericAlignment, genome: &Genome, mapq: u8, + single_end: bool, + hard_clip: bool, ) -> Result, Error> { - let donor = &alignment.donor; - let acceptor = &alignment.acceptor; - - 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, - false, - &acceptor_sa, - )?; - let acceptor_record = build_segment_record( - &alignment.read_name, - &alignment.read_seq, - acceptor, - genome, - mapq, - true, - &donor_sa, - )?; - - Ok(vec![donor_record, acceptor_record]) + use cigar::op::Kind; + + let read_len = alignment.read_seq.len(); + let mut segs = [&alignment.donor, &alignment.acceptor]; + let mut cigars = segs.map(|seg| full_length_cigar(seg, read_len)); + + if single_end { + // STAR orders trChim by read position (roStart); order by the 5' clip of + // the full-length CIGAR so this holds whichever detection tier built the pair. + let ro_start = |seg: &ChimericSegment, ops: &[cigar::Op]| { + let clip = |op: Option<&cigar::Op>| { + op.filter(|op| op.kind() == Kind::SoftClip) + .map_or(0, |op| op.len()) + }; + if seg.is_reverse { + clip(ops.last()) + } else { + clip(ops.first()) + } + }; + if ro_start(segs[0], &cigars[0]) > ro_start(segs[1], &cigars[1]) { + segs.swap(0, 1); + cigars.swap(0, 1); + } + } + + // STAR SE: chimRepresent = (trChim[0].maxScore > trChim[1].maxScore) ? 0 : 1 + let represent = usize::from(single_end && segs[0].score <= segs[1].score); + // Bases hard-clipped from the start/end of SEQ (in CIGAR orientation). + let mut seq_trim = [(0usize, 0usize); 2]; + + if single_end && hard_clip { + let isuppl = 1 - represent; + let seg = segs[isuppl]; + // STAR: alignType = (itr%2 == Str) ? -12 (hard clip right) : -11 (hard clip left) + let hard_right = (isuppl == 1) == seg.is_reverse; + let ops = &mut cigars[isuppl]; + if hard_right { + if let Some(last) = ops.last_mut() + && last.kind() == Kind::SoftClip + { + seq_trim[isuppl].1 = last.len(); + *last = cigar::Op::new(Kind::HardClip, last.len()); + } + } else if let Some(first) = ops.first_mut() + && first.kind() == Kind::SoftClip + { + seq_trim[isuppl].0 = first.len(); + *first = cigar::Op::new(Kind::HardClip, first.len()); + } + } + + let sa = [0, 1].map(|i| format_sa_entry(segs[i], &cigars[i], genome, mapq)); + + let mut records = Vec::with_capacity(2); + for i in 0..2 { + let is_supplementary = i != represent; + let with_seq = single_end || !is_supplementary; + records.push(build_segment_record( + &alignment.read_name, + &alignment.read_seq, + segs[i], + &cigars[i], + if with_seq { Some(seq_trim[i]) } else { None }, + genome, + mapq, + is_supplementary, + &sa[1 - i], + )?); + } + + Ok(records) +} + +/// Segment CIGAR padded with soft clips so that it covers the whole read. +/// +/// Segments built from full-read transcripts already do. Segments re-seeded +/// from a sub-sequence (soft-clip / residual re-mapping) may not; their +/// `read_start` is the left clip in CIGAR orientation. +fn full_length_cigar(seg: &ChimericSegment, read_len: usize) -> Vec { + use cigar::op::Kind; + + let consumes_read = |k: Kind| { + matches!( + k, + Kind::Match + | Kind::Insertion + | Kind::SoftClip + | Kind::SequenceMatch + | Kind::SequenceMismatch + ) + }; + let query_len: usize = seg + .cigar + .iter() + .filter(|op| consumes_read(op.kind())) + .map(|op| op.len()) + .sum(); + if query_len == read_len { + return seg.cigar.clone(); + } + + // Strip existing clips, then re-pad to the full read length. + let core: Vec = seg + .cigar + .iter() + .copied() + .filter(|op| !matches!(op.kind(), Kind::SoftClip | Kind::HardClip)) + .collect(); + let core_len: usize = core + .iter() + .filter(|op| consumes_read(op.kind())) + .map(|op| op.len()) + .sum(); + let left = seg.read_start.min(read_len.saturating_sub(core_len)); + let right = read_len.saturating_sub(left + core_len); + + let mut ops = Vec::with_capacity(core.len() + 2); + if left > 0 { + ops.push(cigar::Op::new(Kind::SoftClip, left)); + } + ops.extend(core); + if right > 0 { + ops.push(cigar::Op::new(Kind::SoftClip, right)); + } + ops } /// 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 +fn format_sa_entry(seg: &ChimericSegment, ops: &[cigar::Op], genome: &Genome, mapq: u8) -> String { + let chr = &genome.chr_name[seg.chr_idx]; + let pos = seg.genome_start - genome.chr_start[seg.chr_idx] + 1; // 1-based per-chr let strand = if seg.is_reverse { '-' } else { '+' }; - let cigar = seg.cigar_string(); + let cigar = cigar_to_string(ops); format!( "{},{},{},{},{},{};", chr, pos, strand, cigar, mapq, seg.n_mismatch @@ -184,10 +296,16 @@ fn format_sa_entry( } /// Build one SAM record for a chimeric segment. +/// +/// `seq_trim` is `None` for an empty SEQ (`*`), or the number of bases to drop +/// from the start / end of the CIGAR-oriented read (hard clips). +#[allow(clippy::too_many_arguments)] fn build_segment_record( read_name: &str, read_seq: &[u8], seg: &ChimericSegment, + ops: &[cigar::Op], + seq_trim: Option<(usize, usize)>, genome: &Genome, mapq: u8, is_supplementary: bool, @@ -219,21 +337,20 @@ fn build_segment_record( *record.mapping_quality_mut() = MappingQuality::new(mapq); - *record.cigar_mut() = seg.cigar.iter().copied().collect(); + *record.cigar_mut() = ops.iter().copied().collect(); - // Primary record carries the full read sequence; supplementary uses * (empty). - if !is_supplementary { - if seg.is_reverse { - let seq_bytes: Vec = read_seq + if let Some((trim_left, trim_right)) = seq_trim { + let oriented: Vec = if seg.is_reverse { + read_seq .iter() .rev() .map(|&b| decode_base(complement_base(b))) - .collect(); - *record.sequence_mut() = Sequence::from(seq_bytes); + .collect() } else { - let seq_bytes: Vec = read_seq.iter().map(|&b| decode_base(b)).collect(); - *record.sequence_mut() = Sequence::from(seq_bytes); - } + read_seq.iter().map(|&b| decode_base(b)).collect() + }; + let end = oriented.len().saturating_sub(trim_right).max(trim_left); + *record.sequence_mut() = Sequence::from(oriented[trim_left..end].to_vec()); // Leave QUAL empty (not available for chimeric segments) *record.quality_scores_mut() = QualityScores::default(); } @@ -527,7 +644,7 @@ mod tests { "READ_001".to_string(), ); let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); + let records = build_within_bam_records(&alignment, &genome, 255, false, true).unwrap(); assert_eq!(records.len(), 2); } @@ -567,7 +684,7 @@ mod tests { "READ_001".to_string(), ); let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); + let records = build_within_bam_records(&alignment, &genome, 255, false, true).unwrap(); let donor_flags = records[0].flags(); let acceptor_flags = records[1].flags(); @@ -618,7 +735,7 @@ mod tests { "READ_001".to_string(), ); let genome = make_genome_2chr(); - let records = build_within_bam_records(&alignment, &genome, 255).unwrap(); + let records = build_within_bam_records(&alignment, &genome, 255, false, true).unwrap(); // Donor record's SA tag should point to acceptor let sa_tag = Tag::new(b'S', b'A'); @@ -676,7 +793,7 @@ mod tests { 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(); + let records = build_within_bam_records(&alignment, &genome, 255, false, true).unwrap(); // Donor has sequence, acceptor has empty sequence (*) assert!( @@ -688,4 +805,104 @@ mod tests { "supplementary record must have empty SEQ" ); } + + fn cigar_str(rec: &RecordBuf) -> String { + cigar_to_string(rec.cigar().as_ref()) + } + + fn sa_str(rec: &RecordBuf) -> String { + match rec.data().get(b"SA") { + Some(Value::String(s)) => s.to_string(), + other => panic!("missing SA tag: {other:?}"), + } + } + + /// Issue #279: SE segments whose CIGAR only covers the segment (e.g. from + /// soft-clip re-seeding) must be padded to the read length, and the + /// supplementary one hard-clipped on its junction side, like STAR. + fn se_alignment(acceptor_reverse: bool) -> ChimericAlignment { + 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), Op::new(Kind::SoftClip, 37)], + score: 63, + n_mismatch: 0, + }; + // Partial CIGAR (37M): only the segment, as built from a sub-sequence. + let acceptor = ChimericSegment { + chr_idx: 1, + genome_start: 600, + genome_end: 637, + is_reverse: acceptor_reverse, + read_start: if acceptor_reverse { 0 } else { 63 }, + read_end: if acceptor_reverse { 37 } else { 100 }, + cigar: vec![Op::new(Kind::Match, 37)], + score: 37, + n_mismatch: 1, + }; + ChimericAlignment::new( + donor, + acceptor, + 0, + 0, + 0, + (0..100u8).map(|i| i % 4).collect(), + "READ_SE".to_string(), + ) + } + + #[test] + fn test_within_bam_se_hard_clip_forward() { + let genome = make_genome_2chr(); + let records = + build_within_bam_records(&se_alignment(false), &genome, 255, true, true).unwrap(); + // Representative = higher-scoring donor, full SEQ, soft clips. + assert!(!records[0].flags().is_supplementary()); + assert_eq!(cigar_str(&records[0]), "63M37S"); + assert_eq!(records[0].sequence().len(), 100); + // Supplementary acceptor: 5' side (junction) hard-clipped, SEQ trimmed. + assert!(records[1].flags().is_supplementary()); + assert_eq!(cigar_str(&records[1]), "63H37M"); + assert_eq!(records[1].sequence().len(), 37); + assert_eq!(sa_str(&records[0]), "chr22,89,+,63H37M,255,1;"); + assert_eq!(sa_str(&records[1]), "chr9,101,+,63M37S,255,0;"); + } + + #[test] + fn test_within_bam_se_hard_clip_reverse_and_soft_clip() { + let genome = make_genome_2chr(); + // Reverse acceptor: CIGAR laid out on the reverse-complemented read, + // so the junction side is on the right. + let records = + build_within_bam_records(&se_alignment(true), &genome, 255, true, true).unwrap(); + assert_eq!(cigar_str(&records[1]), "37M63H"); + assert_eq!(records[1].sequence().len(), 37); + + // SoftClip mode (STAR alignType -13): soft clips and full SEQ. + let records = + build_within_bam_records(&se_alignment(true), &genome, 255, true, false).unwrap(); + assert!(records[1].flags().is_supplementary()); + assert_eq!(cigar_str(&records[1]), "37M63S"); + assert_eq!(records[1].sequence().len(), 100); + } + + #[test] + fn test_within_bam_se_representative_is_higher_score() { + let genome = make_genome_2chr(); + let mut aln = se_alignment(false); + aln.acceptor.score = 70; + let records = build_within_bam_records(&aln, &genome, 255, true, true).unwrap(); + // STAR: chimRepresent = trChim[0].maxScore > trChim[1].maxScore ? 0 : 1 + assert!(records[0].flags().is_supplementary()); + assert!(!records[1].flags().is_supplementary()); + // Donor (5' segment, forward) hard-clipped on its 3' (right) side. + assert_eq!(cigar_str(&records[0]), "63M37H"); + assert_eq!(records[0].sequence().len(), 63); + assert_eq!(cigar_str(&records[1]), "63S37M"); + } } diff --git a/src/lib.rs b/src/lib.rs index 7b60e9b9..7f4655e0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1639,7 +1639,13 @@ 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)?; + let supp = build_within_bam_records( + chim_aln, + &index.genome, + 255, + true, + params.chim_out_bam_hard_clip(), + )?; writer.write_batch(&supp)?; } } @@ -1705,8 +1711,13 @@ 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)?; + let supp = build_within_bam_records( + chim_aln, + &index.genome, + 255, + true, + params.chim_out_bam_hard_clip(), + )?; writer.write_batch(&supp)?; } } @@ -2962,7 +2973,13 @@ 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)?; + let supp = build_within_bam_records( + chim_aln, + &index.genome, + 255, + false, + params.chim_out_bam_hard_clip(), + )?; writer.write_batch(&supp)?; } } @@ -3022,8 +3039,13 @@ 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)?; + let supp = build_within_bam_records( + chim_aln, + &index.genome, + 255, + false, + params.chim_out_bam_hard_clip(), + )?; writer.write_batch(&supp)?; } } diff --git a/src/params/mod.rs b/src/params/mod.rs index 3b507310..dd89b386 100644 --- a/src/params/mod.rs +++ b/src/params/mod.rs @@ -1352,6 +1352,12 @@ impl Parameters { self.chim_out_type.iter().any(|s| s == "WithinBAM") } + /// Whether WithinBAM supplementary records are hard-clipped (STAR's + /// `pCh.out.bamHardClip`): `HardClip` is the default, `SoftClip` turns it off. + pub fn chim_out_bam_hard_clip(&self) -> bool { + !self.chim_out_type.iter().any(|s| s == "SoftClip") + } + /// 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] != "-" diff --git a/tests/alignment_features.rs b/tests/alignment_features.rs index dc94067f..4d950f20 100644 --- a/tests/alignment_features.rs +++ b/tests/alignment_features.rs @@ -2209,3 +2209,127 @@ fn test_read_name_separator_cuts_the_qname_and_is_configurable() { ); } } + +// --------------------------------------------------------------------------- +// Issue #279: single-end `--chimOutType WithinBAM` must write well-formed +// records (CIGAR query length == SEQ length). As in STAR +// (ChimericAlign_chimericBAMoutput.cpp + ReadAlign_alignBAM.cpp), the +// non-representative segment is a supplementary record (0x800) hard-clipped +// on its junction side, and both records carry an SA tag. +// --------------------------------------------------------------------------- + +#[test] +fn test_se_chim_within_bam_records_well_formed() { + use noodles::sam::alignment::record::cigar::op::Kind; + + 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); + + // 100 bp chimeric reads: 60 bp forward from one locus fused to the reverse + // complement of 40 bp from a distant locus (strand switch => chimeric), + // in both orders. + let read_len = 100usize; + let fastq_path = tmpdir.path().join("chim.fq"); + { + let mut f = fs::File::create(&fastq_path).unwrap(); + for i in 0..20usize { + let a = 1000 + i * 150; + let b = 14000 + i * 150; + let mut first = genome[a..a + 60].to_vec(); + first.extend_from_slice(&rc(&genome[b..b + 40])); + let mut second = rc(&genome[b..b + 40]); + second.extend_from_slice(&genome[a..a + 60]); + for (tag, s) in [("a", &first), ("b", &second)] { + writeln!(f, "@chim{i}{tag}").unwrap(); + f.write_all(s).unwrap(); + writeln!(f, "\n+\n{}", "I".repeat(s.len())).unwrap(); + } + } + } + + let run = |name: &str, extra: &[&str]| -> Vec { + let out = tmpdir.path().join(name); + fs::create_dir_all(&out).unwrap(); + let prefix = format!("{}/", out.display()); + let mut args: Vec = [ + "--runMode", + "alignReads", + "--genomeDir", + genome_dir.to_str().unwrap(), + "--readFilesIn", + fastq_path.to_str().unwrap(), + "--outSAMtype", + "BAM", + "Unsorted", + "--chimSegmentMin", + "20", + "--outFileNamePrefix", + &prefix, + "--chimOutType", + "WithinBAM", + ] + .iter() + .map(|s| (*s).to_string()) + .collect(); + args.extend(extra.iter().map(|s| (*s).to_string())); + cargo_bin_cmd!("rustar-aligner") + .args(&args) + .assert() + .success(); + let mut reader = bam::io::Reader::new(fs::File::open(out.join("Aligned.out.bam")).unwrap()); + let header = reader.read_header().expect("BAM header readable"); + reader + .record_bufs(&header) + .map(|r| r.expect("valid BAM record")) + .collect() + }; + + let query_len = |rec: &noodles::sam::alignment::RecordBuf, kinds: &[Kind]| -> usize { + rec.cigar() + .as_ref() + .iter() + .filter(|op| kinds.contains(&op.kind())) + .map(|op| op.len()) + .sum() + }; + let consumes = [ + Kind::Match, + Kind::Insertion, + Kind::SoftClip, + Kind::SequenceMatch, + Kind::SequenceMismatch, + ]; + + for (name, extra, hard) in [ + ("out_hard", &[][..], true), + ("out_soft", &["SoftClip"][..], false), + ] { + let records = run(name, extra); + let mut n_suppl = 0usize; + for rec in &records { + let qlen = query_len(rec, &consumes); + let hclip = query_len(rec, &[Kind::HardClip]); + assert_eq!( + qlen, + rec.sequence().len(), + "{name}: CIGAR query length != SEQ length for {:?}", + rec.name() + ); + assert_eq!(qlen + hclip, read_len, "{name}: {:?}", rec.name()); + if rec.flags().is_supplementary() { + n_suppl += 1; + assert_eq!(hclip > 0, hard, "{name}: supplementary clip type"); + assert!(rec.data().get(b"SA").is_some(), "suppl must carry SA"); + } else { + assert_eq!(hclip, 0, "{name}: only supplementary is hard-clipped"); + } + } + assert!( + n_suppl > 0, + "{name}: expected supplementary chimeric records" + ); + } +}