From 5381739aa512b0724d5bc201f28e9af5939d150c Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Tue, 3 Nov 2026 00:23:33 +0100 Subject: [PATCH] Half-mapped pairs go through the same selection and filters as pairs A pair with one aligned mate kept one location and skipped the multimapper filters. Pairs and single-mate alignments now share one list, score range, primary, mappedFilter and too-many-loci check, and one SAM builder. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 + src/align/read_align.rs | 388 ++++++++++++++++------------------ src/io/sam.rs | 241 ++++++++++----------- src/lib.rs | 253 +++++++--------------- src/quant/mod.rs | 30 +-- src/quant/splice_status.rs | 8 +- src/solo/mod.rs | 20 +- src/wasp/mod.rs | 14 +- tests/half_mapped_multimap.rs | 277 ++++++++++++++++++++++++ 9 files changed, 685 insertions(+), 552 deletions(-) create mode 100644 tests/half_mapped_multimap.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index b18a76e..51a0432 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,12 @@ Sections commonly used: Features, Bug fixes, Other changes. ### Bug fixes +- A pair with one aligned mate is no longer reported as one unique location. Single-mate + alignments now compete with pairs for the score range, the primary and `NH`, and pass the + same `mappedFilter` and `--outFilterMultimapNmax` checks: a mate in a 3-copy repeat gets + `NH:3 MAPQ:1` on all 3 copies. The unmapped mate's record is written only under + `--outSAMunmapped Within`, after the read's alignments; `Within KeepPairs` writes one after + each single-mate alignment. `NH` of a pair also counts alignments past `--outSAMmultNmax`. - The residual differences in DIVERGENCE.md §5 are fixed. On the parity suite's yeast (10k and 100k, SE and PE, chimeric, `--outSAMunmapped Within`), human chr21 (SE, PE) and fusion runs, alignments, `SJ.out.tab` diff --git a/src/align/read_align.rs b/src/align/read_align.rs index 35e450a..601b9e5 100644 --- a/src/align/read_align.rs +++ b/src/align/read_align.rs @@ -239,6 +239,62 @@ pub enum PairedAlignmentResult { }, } +impl PairedAlignmentResult { + /// The score alignments of one read compete on: the pair's, or the single mate's. + pub fn score(&self) -> i32 { + match self { + Self::BothMapped(pa) => pa.combined_wt_score, + Self::HalfMapped { + mapped_transcript, .. + } => mapped_transcript.score, + } + } + + /// Rank in the output order (see [`mult_map_order`]). + pub fn output_rank(&self) -> u32 { + match self { + Self::BothMapped(pa) => pa.output_rank, + Self::HalfMapped { + mapped_transcript, .. + } => mapped_transcript.output_rank, + } + } + + fn output_rank_mut(&mut self) -> &mut u32 { + match self { + Self::BothMapped(pa) => &mut pa.output_rank, + Self::HalfMapped { + mapped_transcript, .. + } => &mut mapped_transcript.output_rank, + } + } + + fn g_length(&self) -> u64 { + match self { + Self::BothMapped(pa) => pa.star_g_length(), + Self::HalfMapped { + mapped_transcript: t, + .. + } => t.genome_end - t.genome_start, + } + } + + /// The transcript of mate 1 and of mate 2, where that mate is aligned. + fn mates(&self) -> [Option<&Transcript>; 2] { + match self { + Self::BothMapped(pa) => [Some(&pa.mate1_transcript), Some(&pa.mate2_transcript)], + Self::HalfMapped { + mapped_transcript, + mate1_is_mapped: true, + } => [Some(mapped_transcript), None], + Self::HalfMapped { + mapped_transcript, + mate1_is_mapped: false, + } => [None, Some(mapped_transcript)], + } + } +} + /// Align a read to the genome. /// /// # Algorithm @@ -697,8 +753,8 @@ type PairedAlignResult = ( /// 3. Tag seeds with mate_id (0=mate1, 1=mate2); adjust read_pos to combined-read coords /// 4. Cluster and stitch combined seeds → WorkingTranscripts spanning both mates /// 5. Split each WT by mate_id → finalize each half → pair -/// 6. Decision tree: dedup → score-range → TooManyLoci → quality filter -/// 7. Half-mapped fallback from single-mate WTs +/// 6. Decision tree over pairs and single mates alike: dedup → score-range → +/// quality filter → TooManyLoci /// /// # Arguments /// * `mate1_seq` - First mate sequence (encoded) @@ -782,9 +838,9 @@ pub fn align_paired_read( let combined_score_threshold = (params.out_filter_score_min_over_lread * (len1 + len2) as f64) as i32; - let mut joint_pairs: Vec = Vec::new(); - let mut single_mate1_transcripts: Vec = Vec::new(); - let mut single_mate2_transcripts: Vec = Vec::new(); + // Pairs and single-mate alignments in one list, in window order: they compete + // for the score range, the primary and NH alike. + let mut alignments: Vec = Vec::new(); // STAR's `trAll` for chimeric detection: per window, the combined-read // transcripts — whole pairs and single mates alike. let mut chim_windows: Vec> = Vec::new(); @@ -956,10 +1012,10 @@ pub fn align_paired_read( )); } if let Some(pair) = pair { - joint_pairs.push(pair); + alignments.push(PairedAlignmentResult::BothMapped(Box::new(pair))); } } else { - // Single-mate WT: save for half-mapped fallback + // Single-mate WT let all_m1 = wt.exons.iter().all(|e| e.mate_id == 0); let all_m2 = wt.exons.iter().all(|e| e.mate_id == 1); if all_m1 { @@ -1000,7 +1056,10 @@ pub fn align_paired_read( combined_len, )); } - single_mate1_transcripts.push(t); + alignments.push(PairedAlignmentResult::HalfMapped { + mapped_transcript: t, + mate1_is_mapped: true, + }); } } else if all_m2 { let (read_slice, orig_rev) = if stitch_is_reverse { @@ -1040,7 +1099,10 @@ pub fn align_paired_read( combined_len, )); } - single_mate2_transcripts.push(t); + alignments.push(PairedAlignmentResult::HalfMapped { + mapped_transcript: t, + mate1_is_mapped: false, + }); } } } @@ -1072,7 +1134,7 @@ pub fn align_paired_read( // rescoring from scratch in the PE frame (`crate::align::pe_overlap`). STAR's semantics // (mirrored here, per `star_pe_overlap.rs` in the sister project STAR-rs) are an // unconditional overwrite: if the merge succeeds and at least one transcript converts, - // the separate-mate `joint_pairs` computed above is replaced outright — no additional + // the separate-mate `alignments` computed above are replaced outright — no additional // score comparison, and the multimap-range filter below is not skipped, it just runs once // on the replaced set instead of twice. if params.pe_overlap_nbases_min > 0 { @@ -1097,7 +1159,7 @@ pub fn align_paired_read( merge_params.out_filter_score_min_over_lread = 0.0; let (merged_transcripts, _merged_chim, _n_mapq, _unmapped, _best) = align_read_inner(&merge.merged, read_name, index, &merge_params, false)?; - let mut converted: Vec = Vec::new(); + let mut converted: Vec = Vec::new(); for t in &merged_transcripts { let Some((m1, m2)) = crate::align::pe_overlap::convert_merged_transcript_to_pe( t, @@ -1118,56 +1180,56 @@ pub fn align_paired_read( + scorer.genomic_length_penalty(combined_span); let is_proper_pair = check_proper_pair(&m1, &m2, params); let insert_size = calculate_insert_size(&m1, &m2); - converted.push(PairedAlignment { - output_rank: 0, - mate1_transcript: m1, - mate2_transcript: m2, - mate1_region: (0, len1), - mate2_region: (0, len2), - is_proper_pair, - insert_size, - combined_wt_score, - }); + converted.push(PairedAlignmentResult::BothMapped(Box::new( + PairedAlignment { + output_rank: 0, + mate1_transcript: m1, + mate2_transcript: m2, + mate1_region: (0, len1), + mate2_region: (0, len2), + is_proper_pair, + insert_size, + combined_wt_score, + }, + ))); } if !converted.is_empty() { pe_overlap_replaced = true; - joint_pairs = converted; + alignments = converted; } } } - // --- Decision tree: dedup, score-filter, quality-filter, then half-mapped fallback --- + // --- Decision tree: dedup, score-filter, quality-filter, too many loci --- - // Step 1: position dedup — remove exact (chr, mate1_pos, mate2_pos, strand, CIGAR) duplicates. + // Step 1: position dedup — remove exact (chr, mate positions, strands, CIGARs) duplicates. // Run dedup BEFORE score-range filter so the backup pool is already deduplicated. // (STAR's ordering is multMapSelect → dedup, but dedup before multMapSelect is equivalent // since removing exact duplicates doesn't change the best score.) // Order-preserving exact dedup (keep first in discovery order). We must not // sort by position here: STAR's primary tie-break falls back to the - // earliest-discovered pair, so discovery order has to survive to the + // earliest-discovered alignment, so discovery order has to survive to the // primary-selection sort below. { // Keyed on the mate positions, with the CIGARs of one key compared - // directly: no CIGAR string is built per pair. - let mut firsts: rustc_hash::FxHashMap<(usize, u64, bool, u64, bool), Vec> = + // directly: no CIGAR string is built per alignment. + type MatePos = Option<(usize, u64, bool)>; + let mut firsts: rustc_hash::FxHashMap<[MatePos; 2], Vec> = rustc_hash::FxHashMap::default(); - let keep: Vec = joint_pairs + let keep: Vec = alignments .iter() .enumerate() - .map(|(i, p)| { - let (m1, m2) = (&p.mate1_transcript, &p.mate2_transcript); + .map(|(i, a)| { + let mates = a.mates(); let same = firsts - .entry(( - m1.chr_idx, - m1.genome_start, - m1.is_reverse, - m2.genome_start, - m2.is_reverse, - )) + .entry(mates.map(|m| m.map(|t| (t.chr_idx, t.genome_start, t.is_reverse)))) .or_default(); let dup = same.iter().any(|&j| { - let q = &joint_pairs[j]; - q.mate1_transcript.cigar == m1.cigar && q.mate2_transcript.cigar == m2.cigar + alignments[j] + .mates() + .iter() + .zip(&mates) + .all(|(q, m)| q.map(|t| &t.cigar) == m.map(|t| &t.cigar)) }); if !dup { same.push(i); @@ -1176,91 +1238,68 @@ pub fn align_paired_read( }) .collect(); let mut keep = keep.into_iter(); - joint_pairs.retain(|_| keep.next().unwrap_or(true)); + alignments.retain(|_| keep.next().unwrap_or(true)); } // Step 2: score-range filter (STAR's multMapSelect). - if !joint_pairs.is_empty() { - let best_score = joint_pairs - .iter() - .map(|pa| pa.combined_wt_score) - .max() - .unwrap_or(0); + if let Some(best_score) = alignments.iter().map(PairedAlignmentResult::score).max() { let score_threshold = best_score - params.out_filter_multimap_score_range; - joint_pairs.retain(|pa| pa.combined_wt_score >= score_threshold); + alignments.retain(|a| a.score() >= score_threshold); } if params.out_multimapper_order == MultimapperOrder::Random { // Under Random, the tied top-score prefix is shuffled with a per-read seed // (deterministic per read, so thread-count invariant). - joint_pairs.sort_by(|a, b| { - b.combined_wt_score.cmp(&a.combined_wt_score).then_with(|| { - ( - a.mate1_transcript.chr_idx, - a.mate1_transcript.genome_start, - a.mate1_transcript.is_reverse, - ) - .cmp(&( - b.mate1_transcript.chr_idx, - b.mate1_transcript.genome_start, - b.mate1_transcript.is_reverse, - )) - }) + alignments.sort_by_key(|a| { + let t = a.mates().into_iter().flatten().next(); + ( + std::cmp::Reverse(a.score()), + t.map(|t| (t.chr_idx, t.genome_start, t.is_reverse)), + ) }); shuffle_tied_prefix( - &mut joint_pairs, - |pa| pa.combined_wt_score, + &mut alignments, + PairedAlignmentResult::score, per_read_seed(params.run_rng_seed, read_name), ); - for (rank, pair) in joint_pairs.iter_mut().enumerate() { - pair.output_rank = rank as u32; + for (rank, alignment) in alignments.iter_mut().enumerate() { + *alignment.output_rank_mut() = rank as u32; } } else { mult_map_order( - &mut joint_pairs, - |pa| pa.combined_wt_score, - PairedAlignment::star_g_length, - |pa| &mut pa.output_rank, + &mut alignments, + PairedAlignmentResult::score, + PairedAlignmentResult::g_length, + PairedAlignmentResult::output_rank_mut, params.out_sam_mult_nmax >= 0, ); } // `--peOverlapNbasesMin` may have replaced the pairs: fold the final ones into // trBest too (the best pair is STAR's trBest whenever a pair exists). - for p in joint_pairs.iter().filter(|_| pe_overlap_replaced) { - let score = p.combined_wt_score; + for a in alignments.iter().filter(|_| pe_overlap_replaced) { + let score = a.score(); if best_tr.is_none_or(|(b, _)| score > b.score) { - best_tr = Some(( - BestTr { - score, - n_mm: p.mate1_transcript.n_mismatch + p.mate2_transcript.n_mismatch, - }, - 0, - )); + let n_mm = a.mates().iter().flatten().map(|t| t.n_mismatch).sum(); + best_tr = Some((BestTr { score, n_mm }, 0)); } } let best_tr = best_tr.map_or_else(BestTr::default, |(b, _)| b); - // Step 4: quality filter (mappedFilter). The best pair is STAR's `trBest` - // whenever a pair exists (a single mate cannot clear the combined-length - // gates), so a failure here unmaps the read with that reason rather than - // falling through to the half-mapped fallback. - if let Some(reason) = filter_paired_transcripts( - &mut joint_pairs, - mate1_seq, - mate2_seq, - &index.genome, - params, - ) { + // Step 4: quality filter (mappedFilter). It looks at the best alignment only, + // and a failure unmaps the read with that reason. + if let Some(reason) = + filter_paired_transcripts(&mut alignments, mate1_seq, mate2_seq, &index.genome, params) + { return Ok((Vec::new(), pe_chimeric, 0, Some(reason), best_tr)); } // Step 5: too-many-loci — STAR checks `multi` AFTER mappedFilter, only when the - // best pair passes the score/match/mismatch gates (ReadAlign_mappedFilter.cpp). - // Pairs that survive quality but exceed outFilterMultimapNmax → too many loci; - // a read whose best pair fails quality was already emptied above (→ unmapped). - if joint_pairs.len() > params.out_filter_multimap_nmax as usize { - let n_loci = joint_pairs.len(); + // best alignment passes the score/match/mismatch gates (ReadAlign_mappedFilter.cpp). + // Alignments that survive quality but exceed outFilterMultimapNmax → too many loci; + // a read whose best alignment fails quality was already emptied above (→ unmapped). + let n_loci = alignments.len(); + if n_loci > params.out_filter_multimap_nmax as usize { return Ok(( Vec::new(), pe_chimeric, @@ -1270,92 +1309,22 @@ pub fn align_paired_read( )); } - if !joint_pairs.is_empty() { - let pe_mapq_n = joint_pairs.len().max(1); - let results = joint_pairs - .into_iter() - .map(|pa| PairedAlignmentResult::BothMapped(Box::new(pa))) - .collect(); - return Ok((results, pe_chimeric, pe_mapq_n, None, best_tr)); - } - - // Half-mapped fallback: report the best-scoring single-mate transcript. - // STAR applies the quality filter to the COMBINED read (Lread-1 = len1+len2), so we - // use the same threshold for each mate here. - let single_mate_threshold = combined_score_threshold.max(params.out_filter_score_min); - - let best_m1 = single_mate1_transcripts - .into_iter() - .filter(|t| t.score >= single_mate_threshold) - .max_by_key(|t| t.score); - let best_m2 = single_mate2_transcripts - .into_iter() - .filter(|t| t.score >= single_mate_threshold) - .max_by_key(|t| t.score); - - match (best_m1, best_m2) { - (Some(t1), None) => Ok(( - vec![PairedAlignmentResult::HalfMapped { - mapped_transcript: t1, - mate1_is_mapped: true, - }], - pe_chimeric, - 1, - None, - best_tr, - )), - (None, Some(t2)) => Ok(( - vec![PairedAlignmentResult::HalfMapped { - mapped_transcript: t2, - mate1_is_mapped: false, - }], - pe_chimeric, - 1, - None, - best_tr, - )), - (Some(t1), Some(t2)) => { - // Both have single-mate alignments but couldn't form a valid pair. - // Report the higher-scoring mate as half-mapped. - if t1.score >= t2.score { - Ok(( - vec![PairedAlignmentResult::HalfMapped { - mapped_transcript: t1, - mate1_is_mapped: true, - }], - pe_chimeric, - 1, - None, - best_tr, - )) - } else { - Ok(( - vec![PairedAlignmentResult::HalfMapped { - mapped_transcript: t2, - mate1_is_mapped: false, - }], - pe_chimeric, - 1, - None, - best_tr, - )) - } - } + if alignments.is_empty() { // STAR splits this case in two (`ReadAlign_mappedFilter.cpp`): a read // with no good window at all is `unmappedOther` ("other"), and only a // read that *had* a window whose best transcript failed the score or // length thresholds is `unmappedShort` ("too short"). Reporting every // unmapped pair as too short empties the "other" bucket, which is what // #48 measured on paired-end data. - (None, None) => { - let reason = if clusters.is_empty() || !any_transcript { - UnmappedReason::Other - } else { - UnmappedReason::TooShort - }; - Ok((Vec::new(), pe_chimeric, 0, Some(reason), best_tr)) - } + let reason = if clusters.is_empty() || !any_transcript { + UnmappedReason::Other + } else { + UnmappedReason::TooShort + }; + return Ok((Vec::new(), pe_chimeric, 0, Some(reason), best_tr)); } + + Ok((alignments, pe_chimeric, n_loci, None, best_tr)) } /// Attempt to pair two per-mate transcripts into a PairedAlignment. @@ -1477,10 +1446,10 @@ pub(crate) fn calculate_insert_size(mate1_trans: &Transcript, mate2_trans: &Tran } } -/// Filter paired transcripts by quality thresholds. +/// Filter a pair's alignments by quality thresholds. /// STAR's mappedFilter applies ALL quality checks to trBest (the highest-scoring transcript). fn filter_paired_transcripts( - paired_alns: &mut Vec, + alignments: &mut Vec, mate1_seq: &[u8], mate2_seq: &[u8], genome: &crate::genome::Genome, @@ -1491,47 +1460,44 @@ fn filter_paired_transcripts( // all transcripts in the score window are included (they affect NH/MAPQ). If trBest fails, // the read is unmapped. // - // Returns why the read is unmapped when the best pair fails, in STAR's + // Returns why the read is unmapped when the best alignment fails, in STAR's // order: too short (score or matched bases), then too many mismatches. - // The pairs arrive sorted by score (and shuffled among ties under + // The alignments arrive sorted by score (and shuffled among ties under // `--outMultimapperOrder Random`), so the first is the primary, as in SE. - // `max_by_key` would return the *last* of tied pairs. - let best_pa = paired_alns.first(); - if let Some(best) = best_pa { - let mate1_len = (best.mate1_region.1 - best.mate1_region.0) as f64; - let mate2_len = (best.mate2_region.1 - best.mate2_region.0) as f64; - let combined_lread_m1 = mate1_len + mate2_len; - let combined_nm = best.mate1_transcript.n_mismatch + best.mate2_transcript.n_mismatch; - let combined_score = best.combined_wt_score; - - if combined_score < params.out_filter_score_min - || combined_score < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32 - { - paired_alns.clear(); - return Some(UnmappedReason::TooShort); - } + // `max_by_key` would return the *last* of tied alignments. + let best = alignments.first()?; + // The combined read's length, also when only one mate aligned. + let combined_lread_m1 = (mate1_seq.len() + mate2_seq.len()) as f64; + let mates = best.mates(); + let combined_score = best.score(); + + if combined_score < params.out_filter_score_min + || combined_score < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32 + { + alignments.clear(); + return Some(UnmappedReason::TooShort); + } - // STAR's `nMatch` (read == genome) and `rLength` (aligned bases) of the - // combined transcript: the mismatch ratio is over `rLength`, not the - // read length (`ReadAlign_mappedFilter.cpp`). - let combined_match = star_n_match(&best.mate1_transcript, mate1_seq, genome) - + star_n_match(&best.mate2_transcript, mate2_seq, genome); - let combined_r_length = - star_r_length(&best.mate1_transcript) + star_r_length(&best.mate2_transcript); - if combined_match < params.out_filter_match_nmin - || combined_match < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32 - { - paired_alns.clear(); - return Some(UnmappedReason::TooShort); - } + // STAR's `nMatch` (read == genome) and `rLength` (aligned bases) of the + // combined transcript: the mismatch ratio is over `rLength`, not the + // read length (`ReadAlign_mappedFilter.cpp`). + let combined_match = mates[0].map_or(0, |t| star_n_match(t, mate1_seq, genome)) + + mates[1].map_or(0, |t| star_n_match(t, mate2_seq, genome)); + if combined_match < params.out_filter_match_nmin + || combined_match < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32 + { + alignments.clear(); + return Some(UnmappedReason::TooShort); + } - if combined_nm > params.out_filter_mismatch_nmax - || f64::from(combined_nm) / f64::from(combined_r_length.max(1)) - > params.out_filter_mismatch_nover_lmax - { - paired_alns.clear(); - return Some(UnmappedReason::TooManyMismatches); - } + let combined_nm: u32 = mates.iter().flatten().map(|t| t.n_mismatch).sum(); + let combined_r_length: u32 = mates.iter().flatten().map(|t| star_r_length(t)).sum(); + if combined_nm > params.out_filter_mismatch_nmax + || f64::from(combined_nm) / f64::from(combined_r_length.max(1)) + > params.out_filter_mismatch_nover_lmax + { + alignments.clear(); + return Some(UnmappedReason::TooManyMismatches); } // Too many loci is the caller's to decide, after these gates and with the // loci count intact; clearing here made it unreachable and turned those diff --git a/src/io/sam.rs b/src/io/sam.rs index 898bf33..5ca2b1d 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -1,12 +1,12 @@ //! SAM/BAM output writer with noodles -use crate::align::read_align::PairedAlignment; +use crate::align::read_align::PairedAlignmentResult; use crate::align::transcript::{Transcript, cigar_to_string}; use crate::error::Error; use crate::genome::Genome; use crate::io::fastq::{complement_base, decode_base}; use crate::junction::encode_motif; use crate::mapq::calculate_mapq; -use crate::params::{Parameters, SamAttributes}; +use crate::params::{OutSamUnmapped, Parameters, SamAttributes}; use crate::stats::{BestTr, UnmappedReason}; use bstr::BString; use noodles::sam; @@ -399,7 +399,7 @@ impl SamWriter { /// * `mate1_qual` - First mate quality scores /// * `mate2_seq` - Second mate sequence (encoded) /// * `mate2_qual` - Second mate quality scores - /// * `paired_alignments` - Paired alignments + /// * `paired_alignments` - The read's alignments: pairs and single mates, primary first /// * `genome` - Genome index /// * `params` - Parameters #[allow(clippy::too_many_arguments)] @@ -414,7 +414,7 @@ impl SamWriter { m1_clip3p: usize, m2_clip5p: usize, m2_clip3p: usize, - paired_alignments: &[PairedAlignment], + paired_alignments: &[PairedAlignmentResult], genome: &Genome, params: &Parameters, n_for_mapq: usize, @@ -448,16 +448,62 @@ impl SamWriter { let mut records = Vec::with_capacity(max_output * 2); let best_score = paired_alignments .iter() - .map(|p| p.combined_wt_score) + .map(PairedAlignmentResult::score) .max() .unwrap_or(i32::MIN); + // Observed from STAR: `Within` writes one unmapped-mate record, after all the + // read's alignments, when the primary has one mate; `Within KeepPairs` writes + // one after every single-mate alignment. + let mut primary_unmapped_mate = None; - for (pair_idx, (is_primary, paired_aln)) in - output_order(paired_alignments, |p| p.output_rank) + for (pair_idx, (is_primary, alignment)) in + output_order(paired_alignments, PairedAlignmentResult::output_rank) .take(max_output) .enumerate() { let hit_index = pair_idx + 1; // 1-based + let paired_aln = match alignment { + PairedAlignmentResult::BothMapped(pa) => pa, + PairedAlignmentResult::HalfMapped { + mapped_transcript, + mate1_is_mapped, + } => { + let (mut mapped, mut unmapped) = Self::build_half_mapped_records( + read_name, + mate1_seq, + mate1_qual, + mate2_seq, + mate2_qual, + m1_clip5p, + m1_clip3p, + m2_clip5p, + m2_clip3p, + mapped_transcript, + *mate1_is_mapped, + genome, + params, + mapq, + n_alignments, + hit_index, + )?; + // Observed from STAR: --outSAMflagOR/AND leave an unmapped mate alone. + set_secondary(&mut mapped, !is_primary); + apply_sam_flag_or_and(&mut mapped, params); + set_secondary(&mut unmapped, !is_primary); + for rec in [&mut mapped, &mut unmapped] { + apply_primary_flag(rec, mapped_transcript.score, best_score, params); + } + records.push(mapped); + match params.out_sam_unmapped { + OutSamUnmapped::WithinKeepPairs => records.push(unmapped), + OutSamUnmapped::Within if is_primary => { + primary_unmapped_mate = Some(unmapped); + } + _ => {} + } + continue; + } + }; // STAR reports the pre-split combined WT score (with length penalty) as AS. // This is stored as combined_wt_score, matching STAR's primaryScore. let combined_score = paired_aln.combined_wt_score; @@ -482,7 +528,7 @@ impl SamWriter { paired_aln.insert_size, params, ), - max_output, // NH = number of reported alignments + n_alignments, // NH counts them all, reported or not (observed from STAR) hit_index, params.out_sam_attr_ih_start, combined_score, @@ -513,7 +559,7 @@ impl SamWriter { -paired_aln.insert_size, // Negative for mate2 (mode 1 default) params, ), - max_output, // NH = number of reported alignments + n_alignments, // NH counts them all, reported or not (observed from STAR) hit_index, params.out_sam_attr_ih_start, combined_score, @@ -538,13 +584,16 @@ impl SamWriter { records.push(rec2); } } + records.extend(primary_unmapped_mate); Ok(records) } - /// Build SAM records for a half-mapped pair (one mate mapped, one unmapped). + /// Build the SAM records of a single-mate alignment of a pair. /// - /// Returns 2 records: mate1 first, mate2 second (regardless of which is mapped). + /// Returns the mapped mate's record and the unmapped mate's. The caller sets + /// the flags that depend on the read's other alignments, and decides whether + /// the unmapped mate's record is written. /// /// **Mapped mate:** Normal alignment with FLAG 0x8 (mate unmapped). /// RNEXT `*`, PNEXT 0, as STAR. @@ -552,8 +601,7 @@ impl SamWriter { /// **Unmapped mate:** FLAG 0x4, uT:A:4, MAPQ 0, RNAME `*`, POS 0, RNEXT/PNEXT /// at the mapped mate. SEQ/QUAL in forward orientation (no RC). #[allow(clippy::too_many_arguments)] - #[allow(clippy::too_many_arguments)] - pub fn build_half_mapped_records( + fn build_half_mapped_records( read_name: &str, mate1_seq: &[u8], mate1_qual: &[u8], @@ -567,13 +615,10 @@ impl SamWriter { mate1_is_mapped: bool, genome: &Genome, params: &Parameters, - n_for_mapq: usize, - ) -> Result, Error> { - let mut records = Vec::with_capacity(2); - - let n_alignments = 1usize; - let effective_n = n_alignments.max(n_for_mapq); - let mapq = calculate_mapq(effective_n, params.out_sam_mapq_unique); + mapq: u8, + n_alignments: usize, + hit_index: usize, + ) -> Result<(RecordBuf, RecordBuf), Error> { let attrs = params.out_sam_attributes; let rg_id_owned = params.primary_rg_id()?; let rg_id = rg_id_owned.as_deref(); @@ -654,7 +699,7 @@ impl SamWriter { if attrs.contains(SamAttributes::HI) { data.insert( Tag::HIT_INDEX, - Value::from(params.out_sam_attr_ih_start as i32), + Value::from((hit_index - 1) as i32 + params.out_sam_attr_ih_start as i32), ); } if attrs.contains(SamAttributes::AS) { @@ -701,7 +746,6 @@ impl SamWriter { data.insert(Tag::new(b'M', b'D'), Value::String(BString::from(md))); } maybe_insert_rg_tag(&mut mapped_rec, rg_id); - apply_sam_flag_or_and(&mut mapped_rec, params); // Soft-clip the mapped mate's clipped bases (STAR convention), matching the // both-mapped path. No-op when the mapped mate has no clip. if mapped_clip5p > 0 || mapped_clip3p > 0 { @@ -763,16 +807,7 @@ impl SamWriter { }, ); - // Order: mate1 first, mate2 second - if mate1_is_mapped { - records.push(mapped_rec); - records.push(unmapped_rec); - } else { - records.push(unmapped_rec); - records.push(mapped_rec); - } - - Ok(records) + Ok((mapped_rec, unmapped_rec)) } /// Build transcriptome-space SAM records for `--quantMode TranscriptomeSAM`. @@ -5134,13 +5169,50 @@ mod tests { ); } + /// The records of a one-mate alignment under `--outSAMunmapped Within`. + fn half_mapped_records( + mapped_transcript: &Transcript, + mate1_is_mapped: bool, + genome: &Genome, + ) -> Vec { + let params = Parameters::parse_from([ + "rustar-aligner", + "--readFilesIn", + "r1.fq", + "r2.fq", + "--outSAMunmapped", + "Within", + ]); + let mate1_seq = [0, 1, 2, 3]; + let mate2_seq = [3, 2, 1, 0]; + let qual = [30, 30, 30, 30]; + SamWriter::build_paired_records( + "read1", + &mate1_seq, + &qual, + &mate2_seq, + &qual, + 0, + 0, + 0, + 0, + &[PairedAlignmentResult::HalfMapped { + mapped_transcript: mapped_transcript.clone(), + mate1_is_mapped, + }], + genome, + ¶ms, + 1, + ) + .unwrap() + } + #[test] fn test_build_half_mapped_flags() { use crate::align::transcript::Exon; use cigar::op::{Kind, Op}; let genome = make_test_genome(); - let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r1.fq", "r2.fq"]); let mapped_transcript = Transcript { output_rank: 0, @@ -5165,29 +5237,8 @@ mod tests { junction_strand: vec![], }; - let mate1_seq = vec![0, 1, 2, 3]; // ACGT - let mate1_qual = vec![30, 30, 30, 30]; - let mate2_seq = vec![3, 2, 1, 0]; // TGCA - let mate2_qual = vec![30, 30, 30, 30]; - // mate1 is mapped - let records = SamWriter::build_half_mapped_records( - "read1", - &mate1_seq, - &mate1_qual, - &mate2_seq, - &mate2_qual, - 0, - 0, - 0, - 0, - &mapped_transcript, - true, // mate1_is_mapped - &genome, - ¶ms, - 1, - ) - .unwrap(); + let records = half_mapped_records(&mapped_transcript, true, &genome); assert_eq!(records.len(), 2); @@ -5214,7 +5265,6 @@ mod tests { use cigar::op::{Kind, Op}; let genome = make_test_genome(); - let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r1.fq", "r2.fq"]); let mapped_transcript = Transcript { output_rank: 0, @@ -5239,28 +5289,7 @@ mod tests { junction_strand: vec![], }; - let mate1_seq = vec![0, 1, 2, 3]; - let mate1_qual = vec![30, 30, 30, 30]; - let mate2_seq = vec![3, 2, 1, 0]; - let mate2_qual = vec![30, 30, 30, 30]; - - let records = SamWriter::build_half_mapped_records( - "read1", - &mate1_seq, - &mate1_qual, - &mate2_seq, - &mate2_qual, - 0, - 0, - 0, - 0, - &mapped_transcript, - true, - &genome, - ¶ms, - 1, - ) - .unwrap(); + let records = half_mapped_records(&mapped_transcript, true, &genome); // mapped_pos = genome_start(2) - chr_start(0) + 1 = 3 let expected_pos = 3usize; @@ -5304,7 +5333,6 @@ mod tests { use cigar::op::{Kind, Op}; let genome = make_test_genome(); - let params = Parameters::parse_from(["rustar-aligner", "--readFilesIn", "r1.fq", "r2.fq"]); let mapped_transcript = Transcript { output_rank: 0, @@ -5329,60 +5357,23 @@ mod tests { junction_strand: vec![], }; - let mate1_seq = vec![0, 1, 2, 3]; - let mate1_qual = vec![30, 30, 30, 30]; - let mate2_seq = vec![3, 2, 1, 0]; - let mate2_qual = vec![30, 30, 30, 30]; - - // When mate1 is mapped: mate1 comes first, mate2 second - let records_m1 = SamWriter::build_half_mapped_records( - "read1", - &mate1_seq, - &mate1_qual, - &mate2_seq, - &mate2_qual, - 0, - 0, - 0, - 0, - &mapped_transcript, - true, - &genome, - ¶ms, - 1, - ) - .unwrap(); + // The mapped mate comes first, whichever it is (observed from STAR). + let records_m1 = half_mapped_records(&mapped_transcript, true, &genome); assert!(records_m1[0].flags().is_first_segment()); // First record = mate1 assert!(records_m1[1].flags().is_last_segment()); // Second record = mate2 - // When mate2 is mapped: mate1 still comes first, mate2 second - let records_m2 = SamWriter::build_half_mapped_records( - "read1", - &mate1_seq, - &mate1_qual, - &mate2_seq, - &mate2_qual, - 0, - 0, - 0, - 0, - &mapped_transcript, - false, - &genome, - ¶ms, - 1, - ) - .unwrap(); - assert!(records_m2[0].flags().is_first_segment()); // First record = mate1 (unmapped) - assert!(records_m2[1].flags().is_last_segment()); // Second record = mate2 (mapped) - assert!(records_m2[0].flags().is_unmapped()); // mate1 is unmapped - assert!(!records_m2[1].flags().is_unmapped()); // mate2 is mapped + let records_m2 = half_mapped_records(&mapped_transcript, false, &genome); + assert!(records_m2[0].flags().is_last_segment()); // First record = mate2 (mapped) + assert!(records_m2[1].flags().is_first_segment()); // Second record = mate1 (unmapped) + assert!(!records_m2[0].flags().is_unmapped()); + assert!(records_m2[1].flags().is_unmapped()); } /// STAR writes the mate on the forward strand first (the mate of the /// transcript's first exon), so mate2 leads when mate1 is reverse. #[test] fn test_build_paired_records_forward_mate_first() { + use crate::align::read_align::PairedAlignment; use crate::align::transcript::Exon; use cigar::op::{Kind, Op}; @@ -5433,7 +5424,7 @@ mod tests { 0, 0, 0, - &[p], + &[PairedAlignmentResult::BothMapped(Box::new(p))], &genome, ¶ms, 1, diff --git a/src/lib.rs b/src/lib.rs index 1d0c993..f483293 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1318,6 +1318,28 @@ struct Stage2 { read_idx: Vec, } +/// Count a mapped pair: by its number of alignments, as half-mapped when the primary +/// has one mate, and with its transcripts' statistics when it is unique. +fn record_pair_alignment_stats( + stats: &crate::stats::AlignmentStats, + results: &[crate::align::read_align::PairedAlignmentResult], + max_multimaps: usize, +) { + use crate::align::read_align::PairedAlignmentResult; + stats.record_alignment(results.len(), max_multimaps); + if matches!( + results.first(), + Some(PairedAlignmentResult::HalfMapped { .. }) + ) { + stats.record_half_mapped(); + } + if results.len() == 1 { + for t in pair_result_transcripts(results) { + stats.record_transcript_stats(t); + } + } +} + /// Every transcript (one per mate) of the alignments of a pair, in order. fn pair_result_transcripts( results: &[crate::align::read_align::PairedAlignmentResult], @@ -2918,7 +2940,7 @@ fn align_reads_solo_pe( sj_stats: &std::sync::Arc, solo_ctx: &std::sync::Arc, ) -> anyhow::Result<()> { - use crate::align::read_align::{PairedAlignment, PairedAlignmentResult, align_paired_read}; + use crate::align::read_align::{PairedAlignmentResult, align_paired_read}; use crate::io::fastq::clip_read; use crate::io::sam::{BufferedSamRecords, SamWriter}; use crate::solo::{SoloCountRecord, SoloMultiRecord}; @@ -3137,49 +3159,19 @@ fn align_reads_solo_pe( stats.record_unmapped_reason( unmapped_reason.unwrap_or(crate::stats::UnmappedReason::Other), ); - } else if has_half_mapped { - stats.record_alignment(1, max_multimaps); - stats.record_half_mapped(); - if let Some(PairedAlignmentResult::HalfMapped { - mapped_transcript, - .. - }) = results.first() - { - stats.record_transcript_stats(mapped_transcript); - } } else { - let n = both_mapped.len(); - stats.record_alignment(n, max_multimaps); - if n == 1 { - stats.record_transcript_stats(&both_mapped[0].mate1_transcript); - stats.record_transcript_stats(&both_mapped[0].mate2_transcript); - } + record_pair_alignment_stats(&stats, &results, max_multimaps); } - let is_unique = - both_mapped.len() == 1 || (has_half_mapped && results.len() == 1); + let is_unique = results.len() == 1; // Per-read junction dedup (fix A): a junction crossed by // both mates / several loci counts once. - let mut read_trs: Vec<&crate::align::transcript::Transcript> = - Vec::new(); - for result in &results { - match result { - PairedAlignmentResult::BothMapped(pair) => { - read_trs.push(&pair.mate1_transcript); - read_trs.push(&pair.mate2_transcript); - } - PairedAlignmentResult::HalfMapped { - mapped_transcript, .. - } => { - read_trs.push(mapped_transcript); - } - } - } + let read_trs = pair_result_transcripts(&results); record_read_junctions(read_trs, &index, &sj_stats, is_unique); // SJ feature: junctions crossed by the unique pair (both mates). let junctions: Vec<(u64, u64)> = - if solo.sj_enabled && both_mapped.len() == 1 { + if solo.sj_enabled && both_mapped.len() == 1 && is_unique { let pair = both_mapped[0]; let mut js = Vec::new(); for tr in [&pair.mate1_transcript, &pair.mate2_transcript] { @@ -3197,8 +3189,8 @@ fn align_reads_solo_pe( }; // Solo quantification: union both mates (strand from mate 1) - // for a both-mapped pair; fall back to the mapped mate for - // half-mapped. + // for a both-mapped pair; the mapped mate for a single-mate + // alignment. // `CB_samTagOut` corrects the barcode for the CB tag and // counts nothing, exactly as in the single-end loop. let (mut outcome, cb_corrected) = if tag_out_only { @@ -3214,28 +3206,10 @@ fn align_reads_solo_pe( corrected }); (outcome, corrected) - } else if !both_mapped.is_empty() { - let pairs: Vec<_> = both_mapped - .iter() - .map(|pa| (&pa.mate1_transcript, &pa.mate2_transcript)) - .collect(); + } else if !results.is_empty() { ( solo.process_read_pe( - &pairs, - pread.barcode.as_ref().map_err(|m| *m), - &junctions, - ), - None, - ) - } else if let Some(PairedAlignmentResult::HalfMapped { - mapped_transcript, - .. - }) = results.first() - { - ( - solo.process_read( - std::slice::from_ref(mapped_transcript), - 1, + &results, pread.barcode.as_ref().map_err(|m| *m), &junctions, ), @@ -3276,37 +3250,7 @@ fn align_reads_solo_pe( buffer.push(record); } } - } else if has_half_mapped { - if let Some(PairedAlignmentResult::HalfMapped { - mapped_transcript, - mate1_is_mapped, - }) = results.first() - { - let records = SamWriter::build_half_mapped_records( - &out_read_name, - &pread.mate1.sequence, - &pread.mate1.quality, - &pread.mate2.sequence, - &pread.mate2.quality, - clip5p_m1, - clip3p_m1, - clip5p_m2, - clip3p_m2, - mapped_transcript, - *mate1_is_mapped, - &index.genome, - params, - n_for_mapq, - )?; - for record in records { - buffer.push(record); - } - } - } else if both_mapped.len() <= max_multimaps { - let paired_alns: Vec = both_mapped - .iter() - .map(|pa| PairedAlignment::clone(pa)) - .collect(); + } else { // Soft-clip the fixed per-mate clips against the original // mate reads (matching SE/PE non-solo). Inert at default 10x. let mut records = SamWriter::build_paired_records( @@ -3319,7 +3263,7 @@ fn align_reads_solo_pe( clip3p_m1, clip5p_m2, clip3p_m2, - &paired_alns, + &results, &index.genome, params, n_for_mapq, @@ -3327,8 +3271,11 @@ fn align_reads_solo_pe( // STARsolo gene tags: a pair is one alignment, so // both mates carry the pair's tags; pairs in output // order, as `build_paired_records` writes them. - if solo_tags - .intersects(crate::params::SamAttributes::SOLO_GENE_TAGS) + // shortcut: a read with a single-mate alignment gets no + // gene tags, upgrade when these tags take such a read. + if !has_half_mapped + && solo_tags + .intersects(crate::params::SamAttributes::SOLO_GENE_TAGS) { let alns: Vec<_> = both_mapped .iter() @@ -3624,7 +3571,7 @@ fn align_reads_paired_end( stats.record_alignment(0, max_multimaps); stats.record_unmapped_reason(crate::stats::UnmappedReason::Other); if let Some(ref q) = quant { - q.count_pe_read(&[], true, false); + q.count_pe_read(&[], 0); } if output_unmapped { // Full original mates for unmapped pairs (STAR convention). @@ -3755,8 +3702,7 @@ fn align_reads_paired_end( // towards the novel junction list, and a read with an // unannotated junction is held (STAR `outFilterBySJout`). if by_sjout && stage2.is_none() { - let is_unique = - both_mapped.len() == 1 || (has_half_mapped && results.len() == 1); + let is_unique = results.len() == 1; record_read_junctions( pair_result_transcripts(&results), &index, @@ -3785,25 +3731,8 @@ fn align_reads_paired_end( stats.record_unmapped_reason( unmapped_reason.unwrap_or(crate::stats::UnmappedReason::Other), ); - } else if has_half_mapped { - // Half-mapped: count as mapped for the mapped mate - stats.record_alignment(1, max_multimaps); - stats.record_half_mapped(); - // Record transcript stats from the mapped mate only - if let Some(PairedAlignmentResult::HalfMapped { - mapped_transcript, .. - }) = results.first() - { - stats.record_transcript_stats(mapped_transcript); - } } else { - // Both-mapped pairs - let n = both_mapped.len(); - stats.record_alignment(n, max_multimaps); - if n == 1 { - stats.record_transcript_stats(&both_mapped[0].mate1_transcript); - stats.record_transcript_stats(&both_mapped[0].mate2_transcript); - } + record_pair_alignment_stats(&stats, &results, max_multimaps); } // Chimeric stats @@ -3816,12 +3745,11 @@ fn align_reads_paired_end( // Dereference Box to get &PairedAlignment slice. let bm_deref: Vec<&crate::align::read_align::PairedAlignment> = both_mapped.iter().map(AsRef::as_ref).collect(); - q.count_pe_read(&bm_deref, results.is_empty(), has_half_mapped); + q.count_pe_read(&bm_deref, results.len()); } // Record junction statistics - let is_unique = - both_mapped.len() == 1 || (has_half_mapped && results.len() == 1); + let is_unique = results.len() == 1; // Per-read junction dedup (fix A): a junction crossed by // both mates / several loci counts once. let read_trs = pair_result_transcripts(&results); @@ -3830,7 +3758,7 @@ fn align_reads_paired_end( // --outWigType bedGraph: both-mapped pairs' reported alignments, // both mates per pair (mate1 sense, mate2 flipped). signal_n_tr is // the pair NH; each mate record is counted in the writer thread. - let signal_n_tr = both_mapped.len(); + let signal_n_tr = results.len(); let signal_contrib: Vec<(crate::align::transcript::Transcript, bool)> = if params.out_wig_bedgraph() && !both_mapped.is_empty() @@ -3871,40 +3799,7 @@ fn align_reads_paired_end( buffer.push(record); } } - } else if has_half_mapped { - // Half-mapped pair - if let Some(PairedAlignmentResult::HalfMapped { - mapped_transcript, - mate1_is_mapped, - }) = results.first() - { - let records = SamWriter::build_half_mapped_records( - &out_read_name, - &paired_read.mate1.sequence, - &paired_read.mate1.quality, - &paired_read.mate2.sequence, - &paired_read.mate2.quality, - m1_clip5p, - m1_clip3p, - m2_clip5p, - m2_clip3p, - mapped_transcript, - *mate1_is_mapped, - &index.genome, - params, - n_for_mapq, - )?; - for record in records { - buffer.push(record); - } - } - } else if both_mapped.len() <= max_multimaps { - // Both-mapped pairs (within multimap limit) - // Extract PairedAlignments for the existing build_paired_records - let paired_alns: Vec = both_mapped - .iter() - .map(|pa| PairedAlignment::clone(pa)) - .collect(); + } else { let mut records = SamWriter::build_paired_records( &out_read_name, &paired_read.mate1.sequence, @@ -3915,38 +3810,45 @@ fn align_reads_paired_end( m1_clip3p, m2_clip5p, m2_clip3p, - &paired_alns, + &results, &index.genome, params, n_for_mapq, )?; - // WASP allele-specific filtering: stamp vW/vA/vG by - // re-mapping the allele-swapped pair (STAR waspMap). - if let Some(ctx) = &*wasp_ctx { - crate::wasp::annotate_records_pe( - &mut records, - &paired_alns, - &m1_seq, - &m2_seq, - &paired_read.name, - &index, - ctx, - params.out_sam_attributes, - )?; - } - // --outSAMsplicingStatus Yes: splicing status tag. - if let Some(tx) = quant.as_ref().and_then(|q| q.splice_tag.as_ref()) { - crate::quant::splice_status::tag_records_pe( - &mut records, - &paired_alns, - tx, - ); + // shortcut: a read with a single-mate alignment gets neither of + // these tags, upgrade when they take such a read. + if !has_half_mapped { + let paired_alns: Vec = both_mapped + .iter() + .map(|pa| PairedAlignment::clone(pa)) + .collect(); + // WASP allele-specific filtering: stamp vW/vA/vG by + // re-mapping the allele-swapped pair (STAR waspMap). + if let Some(ctx) = &*wasp_ctx { + crate::wasp::annotate_records_pe( + &mut records, + &paired_alns, + &m1_seq, + &m2_seq, + &paired_read.name, + &index, + ctx, + params.out_sam_attributes, + )?; + } + // --outSAMsplicingStatus Yes: splicing status tag. + if let Some(tx) = quant.as_ref().and_then(|q| q.splice_tag.as_ref()) { + crate::quant::splice_status::tag_records_pe( + &mut records, + &paired_alns, + tx, + ); + } } for record in records { buffer.push(record); } } - // else: too many loci, skip output // Transcriptome SAM projection (both-mapped pairs only) let transcriptome_records: TrRecords = if let Some(ref tidx) = tr_local { @@ -3966,10 +3868,11 @@ fn align_reads_paired_end( TrRecords::default() }; // Collect unmapped mates for --outReadsUnmapped Fastx. - // Write both mates if: pair is fully unmapped OR half-mapped. + // Write both mates if: pair is fully unmapped OR the primary has one mate. // STAR writes both mates of half-mapped pairs to the unmapped files. let (unmapped_mate1, unmapped_mate2) = if write_unmapped_fastq { - let pair_unmapped = results.is_empty() || has_half_mapped; + let pair_unmapped = + !matches!(results.first(), Some(PairedAlignmentResult::BothMapped(_))); if pair_unmapped { ( vec![(fastx_name1.clone(), m1_seq.clone(), m1_qual.clone())], diff --git a/src/quant/mod.rs b/src/quant/mod.rs index f3c39f6..2b96353 100644 --- a/src/quant/mod.rs +++ b/src/quant/mod.rs @@ -470,24 +470,19 @@ impl GeneCounts { self.count_unique_hit(&transcripts[0], transcripts[0].is_reverse, gene_ann); } - /// Count a paired-end fragment (both mates uniquely mapped). + /// Count a paired-end fragment. `n_alignments` counts single-mate alignments too. /// - /// - `unmapped` or empty with no half-mapped → N_unmapped - /// - `both_mapped.len() > 1` → N_multimapping - /// - `both_mapped.len() == 1` → unique; union of exon overlaps from both mates, + /// - `n_alignments > 1` → N_multimapping + /// - no pair among them → N_unmapped + /// - one pair → unique; union of exon overlaps from both mates, /// strand determined by mate1 pub fn count_pe_read( &self, both_mapped: &[&PairedAlignment], - unmapped: bool, - half_mapped: bool, + n_alignments: usize, gene_ann: &GeneAnnotation, ) { - if unmapped || (both_mapped.is_empty() && half_mapped) { - self.n_unmapped.fetch_add(1, Ordering::Relaxed); - return; - } - if both_mapped.len() > 1 { + if n_alignments > 1 { self.n_multimapping.fetch_add(1, Ordering::Relaxed); return; } @@ -740,23 +735,18 @@ impl QuantContext { } /// Count a read pair in every enabled quantification. - pub fn count_pe_read( - &self, - both_mapped: &[&PairedAlignment], - unmapped: bool, - half_mapped: bool, - ) { + pub fn count_pe_read(&self, both_mapped: &[&PairedAlignment], n_alignments: usize) { if let Some(g) = &self.gene { if let Some(c) = &g.counts { - c.count_pe_read(both_mapped, unmapped, half_mapped, &g.gene_ann); + c.count_pe_read(both_mapped, n_alignments, &g.gene_ann); } - if let Some(c) = &g.coverage { + if let Some(c) = g.coverage.as_ref().filter(|_| n_alignments == 1) { c.coverage.count_pe_read(both_mapped, &g.gene_ann); } } if let Some(v) = &self.splicing { v.counts - .count_pe_read(both_mapped, unmapped, &v.transcriptome); + .count_pe_read(both_mapped, n_alignments, &v.transcriptome); } } diff --git a/src/quant/splice_status.rs b/src/quant/splice_status.rs index 4c3afe1..8dfa5ff 100644 --- a/src/quant/splice_status.rs +++ b/src/quant/splice_status.rs @@ -594,13 +594,13 @@ impl SplicingCounts { pub fn count_pe_read( &self, both_mapped: &[&PairedAlignment], - unmapped: bool, + n_alignments: usize, idx: &TranscriptomeIndex, ) { - if unmapped || both_mapped.is_empty() { - self.n_unmapped.fetch_add(1, Ordering::Relaxed); - } else if both_mapped.len() > 1 { + if n_alignments > 1 { self.n_multimapping.fetch_add(1, Ordering::Relaxed); + } else if both_mapped.is_empty() { + self.n_unmapped.fetch_add(1, Ordering::Relaxed); } else { self.record(&classify_read(&AlignBlocks::from_pair(both_mapped[0]), idx)); } diff --git a/src/solo/mod.rs b/src/solo/mod.rs index d12a3d2..4923bf5 100644 --- a/src/solo/mod.rs +++ b/src/solo/mod.rs @@ -2233,21 +2233,29 @@ impl SoloContext { /// strand — matching STAR's PE quantification (cf. [`crate::quant`]'s /// `count_pe_read`). Each aligned pair counts as one locus. /// - /// `pairs` are the `(mate1, mate2)` transcripts of every BothMapped pair (>1 - /// pair ⇒ multimapper). `junctions` are the pair's crossed junctions. + /// `results` are the read's alignments, pairs and single mates (>1 ⇒ + /// multimapper). `junctions` are the pair's crossed junctions. pub fn process_read_pe( &self, - pairs: &[(&Transcript, &Transcript)], + results: &[crate::align::read_align::PairedAlignmentResult], barcode: Result<&CellBarcode, BarcodeMiss>, junctions: &[(u64, u64)], ) -> SoloReadOutcome { // A mate pair is one alignment in STAR (one `Transcript`, both mates' // blocks), so the Gene feature asks one transcript to contain the pair. - let alns: Vec = pairs + use crate::align::read_align::PairedAlignmentResult; + let alns: Vec = results .iter() - .map(|(m1, m2)| pair_alignment(m1, m2)) + .map(|r| match r { + PairedAlignmentResult::BothMapped(pa) => { + pair_alignment(&pa.mate1_transcript, &pa.mate2_transcript) + } + PairedAlignmentResult::HalfMapped { + mapped_transcript, .. + } => mapped_transcript.clone(), + }) .collect(); - self.process_read(&alns, pairs.len(), barcode, junctions) + self.process_read(&alns, results.len(), barcode, junctions) } } diff --git a/src/wasp/mod.rs b/src/wasp/mod.rs index cf0e241..73f3561 100644 --- a/src/wasp/mod.rs +++ b/src/wasp/mod.rs @@ -492,20 +492,12 @@ pub fn wasp_type_pe( let (results, _chim, _n_for_mapq, _reason, _best) = align_paired_read(&remap1, &remap2, read_name, index, remap_params)?; - let both: Vec<&PairedAlignment> = results - .iter() - .filter_map(|r| match r { - PairedAlignmentResult::BothMapped(pa) => Some(pa.as_ref()), - PairedAlignmentResult::HalfMapped { .. } => None, - }) - .collect(); - if both.is_empty() { + let Some(PairedAlignmentResult::BothMapped(pa)) = results.first() else { return Ok(4); // re-map unmapped / half-mapped / too-many-loci - } - if both.len() > 1 { + }; + if results.len() > 1 { return Ok(5); } - let pa = both[0]; if !same_exons(&pa.mate1_transcript, tr1) || !same_exons(&pa.mate2_transcript, tr2) { return Ok(6); } diff --git a/tests/half_mapped_multimap.rs b/tests/half_mapped_multimap.rs new file mode 100644 index 0000000..563409d --- /dev/null +++ b/tests/half_mapped_multimap.rs @@ -0,0 +1,277 @@ +//! A pair with one aligned mate goes through the same selection and filters as a +//! pair with two. Expected values are STAR 2.7.11b's observed output on this +//! fixture (same genome and reads, `--genomeSAindexNbases 7`). + +use assert_cmd::cargo::cargo_bin_cmd; +use std::fs; +use std::path::Path; +use tempfile::TempDir; + +const COPIES: [usize; 3] = [5000, 10000, 15000]; +const COPY_LEN: usize = 150; + +fn lcg_seq(seed: u32, length: usize) -> Vec { + let mut state = seed; + (0..length) + .map(|_| { + state = state.wrapping_mul(1_103_515_245).wrapping_add(12345); + b"ACGT"[((state >> 16) & 3) as usize] + }) + .collect() +} + +fn reverse_complement(seq: &[u8]) -> Vec { + seq.iter() + .rev() + .map(|b| match b { + b'A' => b'T', + b'C' => b'G', + b'G' => b'C', + b'T' => b'A', + _ => b'N', + }) + .collect() +} + +fn genome() -> Vec { + let mut g = lcg_seq(11, 20_900); + let copy = lcg_seq(23, COPY_LEN); + for start in COPIES { + g[start..start + COPY_LEN].copy_from_slice(©); + } + g +} + +/// A run of one pair, its `Aligned.out.sam` split into fields. +struct Run { + records: Vec>, + log: String, + genes: Option, +} + +fn run(mate1: &[u8], mate2: &[u8], options: &[&str]) -> Run { + let dir = TempDir::new().unwrap(); + let write_read = |name: &str, seq: &[u8]| { + let path = dir.path().join(name); + let qual = "I".repeat(seq.len()); + fs::write( + &path, + format!("@r\n{}\n+\n{qual}\n", String::from_utf8_lossy(seq)), + ) + .unwrap(); + path + }; + let r1 = write_read("r1.fq", mate1); + let r2 = write_read("r2.fq", mate2); + let fasta = dir.path().join("g.fa"); + fs::write( + &fasta, + format!(">c1\n{}\n", String::from_utf8_lossy(&genome())), + ) + .unwrap(); + let gtf = dir.path().join("g.gtf"); + fs::write( + >f, + "c1\tt\texon\t4001\t16000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\";\n", + ) + .unwrap(); + let index = dir.path().join("i"); + fs::create_dir_all(&index).unwrap(); + cargo_bin_cmd!("rustar-aligner") + .args(["--runMode", "genomeGenerate", "--genomeSAindexNbases", "7"]) + .args(["--genomeDir", index.to_str().unwrap()]) + .args(["--genomeFastaFiles", fasta.to_str().unwrap()]) + .args(["--sjdbGTFfile", gtf.to_str().unwrap()]) + .assert() + .success(); + let prefix = format!("{}/", dir.path().display()); + cargo_bin_cmd!("rustar-aligner") + .args(["--genomeDir", index.to_str().unwrap()]) + .args(["--readFilesIn", r1.to_str().unwrap(), r2.to_str().unwrap()]) + .args(["--outFileNamePrefix", &prefix, "--outSAMtype", "SAM"]) + .args(options) + .assert() + .success(); + let read = |name: &str| fs::read_to_string(Path::new(&prefix).join(name)); + Run { + records: read("Aligned.out.sam") + .unwrap() + .lines() + .filter(|l| !l.starts_with('@')) + .map(|l| l.split('\t').map(String::from).collect()) + .collect(), + log: read("Log.final.out").unwrap(), + genes: read("ReadsPerGene.out.tab").ok(), + } +} + +impl Run { + /// FLAG, POS, MAPQ, PNEXT of every record, in file order. + fn core(&self) -> Vec<(u32, u32, u32, u32)> { + self.records + .iter() + .map(|r| { + let n = |i: usize| r[i].parse().unwrap(); + (n(1), n(3), n(4), n(7)) + }) + .collect() + } + + fn tag(&self, record: usize, tag: &str) -> &str { + let found = self.records[record] + .iter() + .find(|f| f.starts_with(tag)) + .unwrap(); + &found[tag.len() + 3..] + } + + fn count(&self, label: &str) -> u64 { + self.log + .lines() + .find_map(|l| { + let (name, value) = l.split_once('|')?; + (name.trim() == label).then(|| value.trim().parse().unwrap()) + }) + .unwrap() + } +} + +fn copy() -> Vec { + lcg_seq(23, COPY_LEN) +} + +fn no_mate() -> Vec { + vec![b'N'; 15] +} + +/// STAR on this fixture: the copies in the order 5001, 15001, 10001, MAPQ 1. +const ON_ALL_COPIES: [(u32, u32, u32, u32); 3] = + [(73, 5001, 1, 0), (329, 15001, 1, 0), (329, 10001, 1, 0)]; + +#[test] +fn mate_in_a_repeat_is_a_multimapper_on_every_copy() { + let run = run(©(), &no_mate(), &[]); + assert_eq!(run.core(), ON_ALL_COPIES); + for (i, hi) in ["1", "2", "3"].iter().enumerate() { + assert_eq!(run.tag(i, "NH"), "3"); + assert_eq!(run.tag(i, "HI"), *hi); + assert_eq!(run.tag(i, "AS"), "148"); + } + assert_eq!(run.count("Uniquely mapped reads number"), 0); + assert_eq!(run.count("Number of reads mapped to multiple loci"), 1); +} + +#[test] +fn unmapped_mate_record_only_on_request() { + let within = run(©(), &no_mate(), &["--outSAMunmapped", "Within"]); + let mut expected = ON_ALL_COPIES.to_vec(); + expected.push((133, 0, 0, 5001)); + assert_eq!(within.core(), expected); + assert_eq!(within.tag(3, "uT"), "4"); + + let keep = run( + ©(), + &no_mate(), + &["--outSAMunmapped", "Within", "KeepPairs"], + ); + assert_eq!( + keep.core(), + [ + (73, 5001, 1, 0), + (133, 0, 0, 5001), + (329, 15001, 1, 0), + (389, 0, 0, 15001), + (329, 10001, 1, 0), + (389, 0, 0, 10001), + ] + ); +} + +#[test] +fn single_mates_are_held_to_the_loci_and_mismatch_limits() { + let loci = run(©(), &no_mate(), &["--outFilterMultimapNmax", "2"]); + assert!(loci.records.is_empty()); + assert_eq!(loci.count("Number of reads mapped to too many loci"), 1); + + let mut mutated = copy(); + for p in (40..120).step_by(30) { + mutated[p] = if mutated[p] == b'A' { b'C' } else { b'A' }; + } + let mismatches = run(&mutated, &no_mate(), &["--outFilterMismatchNmax", "1"]); + assert!(mismatches.records.is_empty()); + // The mismatch limit also stops the extensions, so STAR calls this too short. + assert_eq!(mismatches.count("Number of reads unmapped: too short"), 1); +} + +#[test] +fn pair_and_single_mates_share_the_score_range() { + let g = genome(); + let near = reverse_complement(&g[5250..5270]); + let run = run(©(), &near, &["--outFilterMultimapScoreRange", "30"]); + assert_eq!( + run.core(), + [ + (99, 5001, 1, 5251), + (147, 5251, 1, 5001), + (329, 15001, 1, 0), + (329, 10001, 1, 0), + ] + ); + assert_eq!(run.tag(0, "AS"), "168"); + assert_eq!(run.tag(2, "AS"), "148"); + assert_eq!(run.tag(3, "HI"), "3"); +} + +#[test] +fn nh_counts_alignments_past_the_output_limit() { + let g = genome(); + let near = reverse_complement(&g[5250..5270]); + let run = run( + ©(), + &near, + &[ + "--outFilterMultimapScoreRange", + "30", + "--outSAMmultNmax", + "2", + ], + ); + assert_eq!( + run.core(), + [ + (99, 5001, 1, 5251), + (147, 5251, 1, 5001), + (329, 15001, 1, 0), + ] + ); + for i in 0..3 { + assert_eq!(run.tag(i, "NH"), "3"); + } +} + +#[test] +fn mates_that_cannot_pair_are_two_alignments() { + let g = genome(); + let run = run( + &g[1000..1100], + &g[2000..2100], + &[ + "--outFilterScoreMinOverLread", + "0.3", + "--outFilterMatchNminOverLread", + "0.3", + ], + ); + assert_eq!(run.core(), [(73, 1001, 3, 0), (393, 2001, 3, 0)]); + assert_eq!(run.tag(1, "NH"), "2"); +} + +/// Not compared with STAR: its `--quantMode` fails on this machine. By the rule +/// that an alignment count above one is a multimapper. +#[test] +fn repeat_mate_counts_as_multimapping_in_gene_counts() { + let run = run(©(), &no_mate(), &["--quantMode", "GeneCounts"]); + let genes = run.genes.unwrap(); + assert!(genes.contains("N_multimapping\t1\t1\t1"), "{genes}"); + assert!(genes.contains("N_unmapped\t0\t0\t0"), "{genes}"); +}