From 6a494bd91f213c574a4e59b9d323d84e4bed6b2a Mon Sep 17 00:00:00 2001 From: Rui Xu Date: Thu, 8 Oct 2026 00:46:19 +0800 Subject: [PATCH 1/6] update seqspec-guide docs --- docs/seqspec-guide.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/seqspec-guide.md b/docs/seqspec-guide.md index 8daaefa..87e86ab 100644 --- a/docs/seqspec-guide.md +++ b/docs/seqspec-guide.md @@ -402,6 +402,7 @@ The processing path is determined by the `strand` field and read lengths. 2. **All barcode regions must have whitelists.** Unlike short reads, long-read barcode extraction requires every barcode region to have a non-empty `onlist`. 3. **`min_len` / `max_len` accuracy matters.** They determine the end-window sampling size and spacer lengths in composite alignment. 4. **Target region placement.** The target (`cdna` or `gdna`) separates 5' and 3' processing. Ensure it is correctly positioned in the region order. +5. **Describe the read, not the original library design.** For long reads `precellar` treats the library structure as the read structure directly, so a region that is absent from the FASTQ must not be listed in `library_spec`. Adapters that upstream tooling already trimmed are the common case. ## 9. Multi-Modal Assays @@ -669,6 +670,7 @@ Use this checklist to verify a seqspec YAML before running `precellar`. Each ite - [ ] Total fixed-sequence length ≥ 12 bp in each end-region collection that contains a barcode. (**error** — `BarcodeExtractor::new()`) - [ ] Fixed regions have their actual nucleotide sequence in the `sequence` field, written in 5'→3' forward orientation. (**not checked** — wrong sequences or orientation silently break composite alignment) - [ ] `min_len` and `max_len` are accurate — they determine end-window size and composite alignment spacer lengths. (**not checked**) +- [ ] No region describes sequence that upstream trimming already removed from the FASTQ. Trimmed adapters declared as `fixed` become anchors no read can match. (**not checked** — silently shifts composite alignment) ### RNA Modality From e806be8ec33d6e7c681af260c0b035e456b64ec2 Mon Sep 17 00:00:00 2001 From: Rui Xu Date: Thu, 8 Oct 2026 18:36:17 +0800 Subject: [PATCH 2/6] small fix --- precellar-core/src/transcriptome/align.rs | 114 +++++++++++++++++++++- 1 file changed, 112 insertions(+), 2 deletions(-) diff --git a/precellar-core/src/transcriptome/align.rs b/precellar-core/src/transcriptome/align.rs index 4c83f73..7b5d22a 100644 --- a/precellar-core/src/transcriptome/align.rs +++ b/precellar-core/src/transcriptome/align.rs @@ -670,8 +670,14 @@ impl TxAligner { Some(TxAlignment::Intergenic) } } else { - // Now we know there's exactly one gene, we choose the first non-discordant alignment - let (read1, read2, alignment) = aln.into_iter().next().unwrap(); + // Now we know there's exactly one gene, we choose the first non-discordant + // alignment. `gene_id` is set by any record pair that produced an + // alignment, which need not be the first pair, so the first pair may + // carry none. + let (read1, read2, alignment) = aln + .into_iter() + .find(|(_, _, alignment)| !alignment.is_empty()) + .unwrap(); Some(TxAlignment::PeAligned { read1, read2, @@ -770,10 +776,114 @@ impl TxAligner { #[cfg(test)] mod tests { + use crate::align::MultiMap; use crate::transcriptome::Exons; + use noodles_core::Position; + use noodles_sam::alignment::record::cigar::Op; + use noodles_sam::alignment::record::data::field::tag::Tag; + use noodles_sam::alignment::record::Flags; + use noodles_sam::alignment::record_buf::data::field::value::Value; + use noodles_sam::header::record::value::{map::ReferenceSequence, Map}; + use std::num::NonZeroUsize; use super::*; + #[test] + fn test_align_pe() { + let header = sam::Header::builder() + .add_reference_sequence( + "chr1", + Map::::new(NonZeroUsize::new(10_000).unwrap()), + ) + .build(); + + let transcript = Transcript { + id: "tx1".to_string(), + chrom: "chr1".to_string(), + start: 100, + end: 600, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { + start: 100, + end: 200, + }, + Exon { + start: 300, + end: 400, + }, + Exon { + start: 500, + end: 600, + }, + ]), + }; + + let record = |start: u64, cigar: Vec<(Kind, usize)>, is_read2: bool| { + let mut rec = RecordBuf::default(); + *rec.reference_sequence_id_mut() = Some(0); + *rec.alignment_start_mut() = Position::new(start as usize + 1); + *rec.flags_mut() = if is_read2 { + Flags::SEGMENTED | Flags::REVERSE_COMPLEMENTED + } else { + Flags::SEGMENTED + }; + *rec.cigar_mut() = cigar + .into_iter() + .map(|(kind, len)| Op::new(kind, len)) + .collect(); + rec.data_mut() + .insert(Tag::CELL_BARCODE_ID, Value::String("ACGT".into())); + rec + }; + + let discordant = |is_read2: bool| { + record( + 120, + vec![(Kind::Match, 60), (Kind::Skip, 140), (Kind::Match, 60)], + is_read2, + ) + }; + let intronic = |is_read2: bool| record(220, vec![(Kind::Match, 30)], is_read2); + + let read1 = MultiMap::new(discordant(false), Some(vec![intronic(false)])); + let read2 = MultiMap::new(discordant(true), Some(vec![intronic(true)])); + + let aligner = TxAligner::new( + vec![transcript], + header, + Some(ChemistryStrandedness::Forward), + ); + let result = aligner + ._align_pe(&read1, &read2, ChemistryStrandedness::Forward) + .expect("the secondary pair aligns, so the read is not dropped"); + + let TxAlignment::PeAligned { + read1: reported1, + read2: reported2, + alignment, + gene_id, + .. + } = &result + else { + panic!("expected PeAligned, got {result:?}"); + }; + + assert_eq!(gene_id, "gene1"); + assert_eq!( + alignment.as_slice(), + &[("tx1".to_string(), AlignType::Intronic)] + ); + assert_eq!(reported1.start(), 220); + assert_eq!(reported2.start(), 220); + + assert!(!result.is_exonic_only()); + assert!(result.is_intronic_only()); + assert!(!result.is_spanning()); + } + #[test] fn test_align() { let transcript = Transcript { From 92f7ec7b7667bf529ec3f090c21c48f09ff35545 Mon Sep 17 00:00:00 2001 From: Rui Xu Date: Sun, 11 Oct 2026 02:10:23 +0800 Subject: [PATCH 3/6] add ambiguous layer --- precellar-core/src/qc.rs | 2 + precellar-core/src/transcriptome/align.rs | 14 +++ precellar-core/src/transcriptome/annotate.rs | 110 ++++++++++++++++-- precellar-core/src/transcriptome/mod.rs | 2 +- .../src/transcriptome/quantification.rs | 60 +++++++--- 5 files changed, 167 insertions(+), 21 deletions(-) diff --git a/precellar-core/src/qc.rs b/precellar-core/src/qc.rs index a78b230..7f65622 100644 --- a/precellar-core/src/qc.rs +++ b/precellar-core/src/qc.rs @@ -415,6 +415,7 @@ pub struct QcGeneQuant { pub(crate) num_unique_umi: u64, pub(crate) num_spliced: u64, pub(crate) num_unspliced: u64, + pub(crate) num_ambiguous: u64, } impl QcGeneQuant { @@ -462,6 +463,7 @@ impl From for Value { "frac_duplicates": 1.0 - qc.num_unique_umi as f64 / num_transcriptomic as f64, "frac_spliced": qc.num_spliced as f64 / qc.num_unique_umi as f64, "frac_unspliced": qc.num_unspliced as f64 / qc.num_unique_umi as f64, + "frac_ambiguous": qc.num_ambiguous as f64 / qc.num_unique_umi as f64, }, }) } diff --git a/precellar-core/src/transcriptome/align.rs b/precellar-core/src/transcriptome/align.rs index 7b5d22a..c325d07 100644 --- a/precellar-core/src/transcriptome/align.rs +++ b/precellar-core/src/transcriptome/align.rs @@ -141,6 +141,20 @@ impl SplicedRecord { })) } + /// A record with the given segments, for tests that exercise logic keyed on + /// the transcript alignment list rather than on read coordinates. + #[cfg(test)] + pub(crate) fn for_test(chrom: &str, segments: Vec<(u64, u64)>) -> Self { + Self { + chrom: chrom.to_string(), + segments: segments + .into_iter() + .map(|(start, end)| Segment { start, end }) + .collect(), + strand: bed_utils::bed::Strand::Forward, + } + } + /// The leftmost position of the alignment. fn start(&self) -> u64 { self.segments.first().map_or(0, |segment| segment.start) diff --git a/precellar-core/src/transcriptome/annotate.rs b/precellar-core/src/transcriptome/annotate.rs index 7ccf1ee..60b233d 100644 --- a/precellar-core/src/transcriptome/annotate.rs +++ b/precellar-core/src/transcriptome/annotate.rs @@ -15,6 +15,14 @@ pub enum Annotation { Ambiguous, } +/// Per-gene UMI counts of one cell, split by splicing status. +#[derive(Debug, Default)] +pub struct SplicedCounts { + pub spliced: Vec<(GeneIndex, u32)>, + pub unspliced: Vec<(GeneIndex, u32)>, + pub ambiguous: Vec<(GeneIndex, u32)>, +} + /// A group of transcript alignments corresponding to a single UMI. struct AlignmentGroup<'a>(Vec<&'a TxAlignment>); @@ -156,12 +164,16 @@ impl<'a> GeneCounter<'a> { .collect() } - pub fn to_spliced_counts(&self) -> (Vec<(GeneIndex, u32)>, Vec<(GeneIndex, u32)>) { - let mut spliced_counts = Vec::new(); - let mut unspliced_counts = Vec::new(); + /// Splits the UMI counts of this cell by splicing status. + /// + /// The three categories are mutually exclusive and together account for + /// every UMI, so their sums equal [`GeneCounter::to_counts`]. + pub fn to_spliced_counts(&self) -> SplicedCounts { + let mut counts = SplicedCounts::default(); self.0.iter().for_each(|(gene_idx, umi_group)| { let mut spliced: u32 = 0; let mut unspliced: u32 = 0; + let mut ambiguous: u32 = 0; umi_group .iter() .for_each(|reads| match reads.to_annotation() { @@ -171,16 +183,21 @@ impl<'a> GeneCounter<'a> { Annotation::Unspliced => { unspliced = unspliced.checked_add(1).unwrap(); } - Annotation::Ambiguous => {} + Annotation::Ambiguous => { + ambiguous = ambiguous.checked_add(1).unwrap(); + } }); if spliced > 0 { - spliced_counts.push((*gene_idx, spliced)); + counts.spliced.push((*gene_idx, spliced)); } if unspliced > 0 { - unspliced_counts.push((*gene_idx, unspliced)); + counts.unspliced.push((*gene_idx, unspliced)); + } + if ambiguous > 0 { + counts.ambiguous.push((*gene_idx, ambiguous)); } }); - (spliced_counts, unspliced_counts) + counts } pub fn to_fragments(&self) -> impl Iterator + '_ { @@ -189,3 +206,82 @@ impl<'a> GeneCounter<'a> { .flat_map(|(_, umi_group)| umi_group.iter().flat_map(|reads| reads.to_fragments())) } } + +#[cfg(test)] +mod tests { + use super::*; + use crate::transcriptome::align::{AlignType, SplicedRecord}; + + /// One read of a molecule, given as the per-transcript alignment types it + /// is compatible with. Coordinates are irrelevant here: `to_annotation` + /// reads only the `(transcript_id, AlignType)` list. + fn read(umi: &str, alignment: &[(&str, AlignType)]) -> TxAlignment { + TxAlignment::SeAligned { + read: SplicedRecord::for_test("chr1", vec![(0, 10)]), + alignment: alignment + .iter() + .map(|(tx, ty)| (tx.to_string(), ty.clone())) + .collect(), + gene_id: "gene1".to_string(), + barcode: "ACGT".to_string(), + umi: Some(umi.to_string()), + } + } + + /// The three categories partition the UMIs of a gene, so the layers sum to + /// the total count. + #[test] + fn splicing_status_partitions_the_umis_of_a_gene() { + let genes = IndexMap::from([("gene1".to_string(), "GeneName1".to_string())]); + + let alignments = vec![ + // Spliced: every transcript is exonic (rule 2). + read("AAAAAAAA", &[("tx1", AlignType::Exonic)]), + // Unspliced: a spanning read decides the molecule (rule 1). + read("CCCCCCCC", &[("tx1", AlignType::Spanning)]), + // Ambiguous: tx1 says exonic, tx2 says intronic, and no read spans, + // so neither rule 2 nor rule 3 holds. + read( + "GGGGGGGG", + &[("tx1", AlignType::Exonic), ("tx2", AlignType::Intronic)], + ), + ]; + + let counter = GeneCounter::new(&alignments, &genes); + let counts = counter.to_spliced_counts(); + + assert_eq!(counts.spliced, vec![(0, 1)]); + assert_eq!(counts.unspliced, vec![(0, 1)]); + assert_eq!(counts.ambiguous, vec![(0, 1)]); + + // The invariant that makes the three layers comparable: nothing is + // dropped, so `X == spliced + unspliced + ambiguous`. + let total: u32 = counter.to_counts().iter().map(|(_, c)| c).sum(); + let split: u32 = counts + .spliced + .iter() + .chain(&counts.unspliced) + .chain(&counts.ambiguous) + .map(|(_, c)| c) + .sum(); + assert_eq!(total, 3); + assert_eq!(total, split); + } + + /// A molecule with no spanning read is still unspliced when every + /// compatible transcript has intronic evidence. + #[test] + fn intronic_evidence_on_all_transcripts_is_unspliced() { + let genes = IndexMap::from([("gene1".to_string(), "GeneName1".to_string())]); + let alignments = vec![read( + "AAAAAAAA", + &[("tx1", AlignType::Intronic), ("tx2", AlignType::Intronic)], + )]; + + let counts = GeneCounter::new(&alignments, &genes).to_spliced_counts(); + + assert_eq!(counts.unspliced, vec![(0, 1)]); + assert!(counts.spliced.is_empty()); + assert!(counts.ambiguous.is_empty()); + } +} diff --git a/precellar-core/src/transcriptome/mod.rs b/precellar-core/src/transcriptome/mod.rs index e71954e..859b7c1 100644 --- a/precellar-core/src/transcriptome/mod.rs +++ b/precellar-core/src/transcriptome/mod.rs @@ -11,7 +11,7 @@ pub mod gtf; mod quantification; pub use align::{TxAlignResult, TxAligner, TxAlignment}; -pub use annotate::GeneCounter; +pub use annotate::{GeneCounter, SplicedCounts}; pub use quantification::Quantifier; use anyhow::{bail, ensure, Result}; diff --git a/precellar-core/src/transcriptome/quantification.rs b/precellar-core/src/transcriptome/quantification.rs index dc9c87d..d27cd11 100644 --- a/precellar-core/src/transcriptome/quantification.rs +++ b/precellar-core/src/transcriptome/quantification.rs @@ -67,6 +67,8 @@ impl Quantifier { ExternalChunkBuilder::new(tempfile::tempfile()?, 3)?; let mut unspliced_cache: ExternalChunkBuilder>> = ExternalChunkBuilder::new(tempfile::tempfile()?, 3)?; + let mut ambiguous_cache: ExternalChunkBuilder>> = + ExternalChunkBuilder::new(tempfile::tempfile()?, 3)?; let mut coverage_cache: ExternalChunkBuilder>> = ExternalChunkBuilder::new(tempfile::tempfile()?, 3)?; @@ -116,11 +118,12 @@ impl Quantifier { .chunk_by(|(chunk_id, _)| *chunk_id); aln_chunks.into_iter().for_each(|(_, chunk)| { - let ((t, s), (uns, cov)) = self.count_one_cell(&genome_index, chunk.map(|x| x.1)); - umi_cache.add(t).unwrap(); - spliced_cache.add(s).unwrap(); - unspliced_cache.add(uns).unwrap(); - coverage_cache.add(cov).unwrap(); + let counts = self.count_one_cell(&genome_index, chunk.map(|x| x.1)); + umi_cache.add(counts.umi).unwrap(); + spliced_cache.add(counts.spliced).unwrap(); + unspliced_cache.add(counts.unspliced).unwrap(); + ambiguous_cache.add(counts.ambiguous).unwrap(); + coverage_cache.add(counts.coverage).unwrap(); }); bar.abandon(); @@ -156,6 +159,16 @@ impl Quantifier { }); adata.layers().add_iter("unspliced", unspliced_count)?; + let ambiguous_count = ambiguous_cache.finish()?.map(|x| { + let x = x.unwrap(); + qc.num_ambiguous += x + .iter() + .flat_map(|x| x.iter().map(|(_, c)| *c as u64)) + .sum::(); + make_csr(x, num_cols) + }); + adata.layers().add_iter("ambiguous", ambiguous_count)?; + let coverage = coverage_cache.finish()?.map(|x| { let x = x.unwrap(); make_csr(x, genome_index.len()) @@ -175,30 +188,51 @@ impl Quantifier { Ok(qc) } + /// Counts a batch of cells, one entry per cell in each returned vector. fn count_one_cell( &self, index: &GenomeBaseIndex, alignments: impl Iterator>, - ) -> ( - (Vec>, Vec>), - (Vec>, Vec>), - ) { + ) -> CellCounts { let alignments = alignments.collect::>(); let chunk_size = (alignments.len() / self.num_threads).max(1); - alignments + let per_cell: Vec<_> = alignments .par_chunks(chunk_size) .flat_map_iter(|chunk| { chunk.into_iter().map(|aln| { let umi_group = GeneCounter::new(aln, &self.genes); let coverage = index.count_fragments(umi_group.to_fragments()); - let (spliced, unspliced) = umi_group.to_spliced_counts(); - ((umi_group.to_counts(), spliced), (unspliced, coverage)) + ( + umi_group.to_counts(), + umi_group.to_spliced_counts(), + coverage, + ) }) }) - .unzip() + .collect(); + + let mut counts = CellCounts::default(); + for (umi, spliced, coverage) in per_cell { + counts.umi.push(umi); + counts.spliced.push(spliced.spliced); + counts.unspliced.push(spliced.unspliced); + counts.ambiguous.push(spliced.ambiguous); + counts.coverage.push(coverage); + } + counts } } +/// Counts for a batch of cells, each vector holding one entry per cell. +#[derive(Default)] +struct CellCounts { + umi: Vec>, + spliced: Vec>, + unspliced: Vec>, + ambiguous: Vec>, + coverage: Vec>, +} + fn sort_alignments( alignments: I, chunk_size: usize, From ce358607f56ed25c89ff60f6d72ae89759a9b875 Mon Sep 17 00:00:00 2001 From: Rui Xu Date: Mon, 12 Oct 2026 13:32:59 +0800 Subject: [PATCH 4/6] add terminal exon exemption and micro-intron tolerance to transcript alignment --- precellar-core/src/transcriptome/align.rs | 378 ++++++++++++++++++++-- precellar-core/src/transcriptome/mod.rs | 2 +- src/align.rs | 7 +- src/sinks.rs | 31 +- 4 files changed, 394 insertions(+), 24 deletions(-) diff --git a/precellar-core/src/transcriptome/align.rs b/precellar-core/src/transcriptome/align.rs index c325d07..a680eb6 100644 --- a/precellar-core/src/transcriptome/align.rs +++ b/precellar-core/src/transcriptome/align.rs @@ -49,6 +49,46 @@ impl AlignType { } } +/// Tolerances applied when matching read structure against a transcript model. +#[derive(Debug, Copy, Clone, PartialEq, Eq)] +pub struct TxAlignOpts { + /// Annotated introns at or below this length are treated as alignment + /// artifacts rather than retained introns: aligners routinely fail to + /// report short skips, so a read covering one looks like intron retention + /// when it is not. Zero-length introns (abutting exons) are always + /// exempted, since they contain no intronic sequence at all. IsoQuant + /// uses 50. + pub micro_intron_length: u64, + /// Tolerance for matching observed junctions against annotated exon + /// boundaries. Long-read indel noise routinely shifts a junction by a + /// few bases, which must not be read as a structural mismatch. Only the + /// inner edges of the alignment are affected; transcript ends are + /// governed by the unconditional exemption above. IsoQuant uses 4 for + /// PacBio and 6 for ONT. Not yet implemented; reserved. + pub junction_tolerance: u64, +} + +impl Default for TxAlignOpts { + fn default() -> Self { + Self { + micro_intron_length: 50, + junction_tolerance: 0, + } + } +} + +impl TxAlignOpts { + /// Options that disable every tolerance, restoring exact structural + /// matching. Used to reproduce historical classifications. Zero-length + /// introns remain exempted, as they hold no intronic sequence. + pub fn exact() -> Self { + Self { + micro_intron_length: 0, + junction_tolerance: 0, + } + } +} + #[derive(Eq, PartialEq, Debug, Clone, Encode, Decode)] pub struct TxAlignResult { pub gene_id: String, @@ -170,6 +210,7 @@ impl SplicedRecord { &self, transcript: &Transcript, strandedness: ChemistryStrandedness, + opts: &TxAlignOpts, ) -> TxAlignResult { let is_antisense = match strandedness { ChemistryStrandedness::Forward => transcript.strand != self.strand, @@ -186,7 +227,7 @@ impl SplicedRecord { gene_id: transcript.gene_id.clone(), transcript_id: transcript.id.clone(), strand, - align_type: self.align_helper(transcript), + align_type: self.align_helper(transcript, opts), } } @@ -224,12 +265,20 @@ impl SplicedRecord { /// Exonic: XXXXX----XXXX--------XXXXX /// Intronic: XX /// Spanning: XXXXXXXXXXXX-------XXXXX - fn align_helper(&self, transcript: &Transcript) -> AlignType { + /// + /// Outer transcript edges are exempt from the containment checks: TSS/polyA + /// annotations (also long-read end trimming) are unreliable, so a read overhanging + /// a terminal exon says nothing about splicing. Inner edges (junctions) are still + /// checked exactly. The exemption is keyed on the transcript's terminal exons, not + /// the first overlapped exon, so that a read starting in intron 2 still counts as + /// intron retention rather than an overhang. + fn align_helper(&self, transcript: &Transcript, opts: &TxAlignOpts) -> AlignType { if self.end() <= transcript.start || self.start() >= transcript.end { panic!("Read does not overlap transcript"); } - let exons = find_exons(transcript.exons(), self.start(), self.end() - 1); + let tx_exons = transcript.exons(); + let exons = find_exons(tx_exons, self.start(), self.end() - 1); if exons.is_empty() { // No exon overlaps the read, but the read overlaps the transcript. return AlignType::Intronic; @@ -238,10 +287,19 @@ impl SplicedRecord { let segments = &self.segments; let n_ex = exons.len(); let n_segment = segments.len(); + + // The left outer edge is exempt only at the transcript's start, and the + // right outer edge only at the transcript's end. Exon boundaries are + // unique within a transcript (exons are non-overlapping and in order), + // so coordinate equality identifies the terminal exons. + let at_tx_start = exons[0].start == tx_exons[0].start; + let at_tx_end = exons[n_ex - 1].end == tx_exons[tx_exons.len() - 1].end; if n_ex == n_segment { if n_ex == 1 { - if segments[0].start >= exons[0].start && segments[0].end <= exons[0].end { - // Fully contained within the exon. + if (at_tx_start || segments[0].start >= exons[0].start) + && (at_tx_end || segments[0].end <= exons[0].end) + { + // Fully contained within the exon (terminal edges exempt). AlignType::Exonic } else { AlignType::Spanning @@ -256,8 +314,8 @@ impl SplicedRecord { && segments[n_segment - 1].start == exons[n_ex - 1].start { // first and last segments also match. - if segments[0].start >= exons[0].start - && segments[n_segment - 1].end <= exons[n_ex - 1].end + if (at_tx_start || segments[0].start >= exons[0].start) + && (at_tx_end || segments[n_segment - 1].end <= exons[n_ex - 1].end) { AlignType::Exonic } else { @@ -295,10 +353,34 @@ impl SplicedRecord { && right_idx == left_idx + 1 }); - if is_compatible { - AlignType::Spanning - } else { + if !is_compatible { AlignType::Discordant + } else { + // Every observed junction lands on adjacent exon boundaries, so + // the read differs from a plain exonic match only by the introns + // bridged within single segments. An intron short enough that + // aligners routinely fail to report it as a skip is an + // alignment artifact, not retained-intron evidence; if the read + // bridges no genuine intron, it is exonic. Zero-length introns + // (abutting exons) never count: they hold no intronic sequence. + let only_artifacts = exons + .windows(2) + .filter(|pair| { + let (left, right) = (&pair[0], &pair[1]); + // The intron is bridged by a single segment that aligns + // straight through it without reporting a junction. A + // segment that only reaches partway into the intron is + // genuine retained-intron sequence, not a missed skip. + segments + .iter() + .any(|segment| segment.start <= left.end && segment.end >= right.start) + }) + .all(|pair| pair[1].start - pair[0].end <= opts.micro_intron_length); + if only_artifacts { + AlignType::Exonic + } else { + AlignType::Spanning + } } } else { // There are more segments than exons. Cannot be concordant. @@ -453,6 +535,7 @@ pub struct TxAligner { transcripts: GIntervalMap, strandedness: Option, header: sam::Header, + opts: TxAlignOpts, } impl TxAligner { @@ -470,9 +553,16 @@ impl TxAligner { transcripts, header, strandedness, + opts: TxAlignOpts::default(), } } + /// Applies custom tolerance options to the annotator. + pub fn with_opts(mut self, opts: TxAlignOpts) -> Self { + self.opts = opts; + self + } + pub fn transcripts(&self) -> impl Iterator { self.transcripts.iter().map(|(_, t)| t) } @@ -711,7 +801,7 @@ impl TxAligner { let region = GenomicRange::new(&rec.chrom, rec.start() as u64, rec.end() as u64); self.transcripts .find(®ion) - .map(|(_, transcript)| rec.align_transcript(transcript, strandedness)) + .map(|(_, transcript)| rec.align_transcript(transcript, strandedness, &self.opts)) .collect() } @@ -933,33 +1023,283 @@ mod tests { strand: bed_utils::bed::Strand::Forward, }; - let mut test_fun = |segments: Vec<(u64, u64)>, expected: AlignType| { + let mut test_fun = |segments: Vec<(u64, u64)>, opts: &TxAlignOpts, expected: AlignType| { record.segments = segments .into_iter() .map(|(s, e)| Segment { start: s, end: e }) .collect(); - assert_eq!(record.align_helper(&transcript), expected); + assert_eq!(record.align_helper(&transcript, opts), expected); }; + let exact = TxAlignOpts::exact(); + let default = TxAlignOpts::default(); // Incompatible - test_fun(vec![(12, 18), (32, 38)], AlignType::Discordant); - test_fun(vec![(21, 29), (32, 38)], AlignType::Discordant); + test_fun(vec![(12, 18), (32, 38)], &exact, AlignType::Discordant); + test_fun(vec![(12, 18), (32, 38)], &default, AlignType::Discordant); + test_fun(vec![(21, 29), (32, 38)], &exact, AlignType::Discordant); + test_fun(vec![(21, 29), (32, 38)], &default, AlignType::Discordant); + test_fun( + vec![(5, 40), (50, 60), (500, 700), (1400, 1501)], + &exact, + AlignType::Discordant, + ); test_fun( vec![(5, 40), (50, 60), (500, 700), (1400, 1501)], + &default, AlignType::Discordant, ); // Exonic - test_fun(vec![(12, 20), (30, 38)], AlignType::Exonic); + test_fun(vec![(12, 20), (30, 38)], &exact, AlignType::Exonic); + test_fun(vec![(12, 20), (30, 38)], &default, AlignType::Exonic); // Intronic - test_fun(vec![(6, 8)], AlignType::Intronic); + test_fun(vec![(6, 8)], &exact, AlignType::Intronic); + test_fun(vec![(6, 8)], &default, AlignType::Intronic); + + // Reaching left of a middle exon is genuine intron retention. + test_fun(vec![(29, 38)], &exact, AlignType::Spanning); + test_fun(vec![(29, 38)], &default, AlignType::Spanning); - // Spanning - test_fun(vec![(29, 38)], AlignType::Spanning); + // The first segment bridges a 10bp intron (exons 10-20 and 30-40): + // an aligner artifact under the default micro-intron exemption, + // retained-intron evidence under exact matching. test_fun( vec![(5, 40), (50, 60), (500, 700), (1500, 1501)], + &exact, + AlignType::Spanning, + ); + test_fun( + vec![(5, 40), (50, 60), (500, 700), (1500, 1501)], + &default, + AlignType::Exonic, + ); + } + + /// The outer edge of the transcript's first and last exon is TSS/polyA + /// territory, not splicing evidence, and is exempt unconditionally. The + /// same edges on interior exons face introns and stay checked. + #[test] + fn test_terminal_exemption() { + let transcript = Transcript { + id: "tx1".to_string(), + chrom: "chr1".to_string(), + start: 100, + end: 1100, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { + start: 100, + end: 200, + }, + Exon { + start: 500, + end: 600, + }, + Exon { + start: 1000, + end: 1100, + }, + ]), + }; + + let mut record = SplicedRecord { + chrom: "chr1".to_string(), + segments: vec![], + strand: bed_utils::bed::Strand::Forward, + }; + let mut test_fun = |segments: Vec<(u64, u64)>, expected: AlignType| { + record.segments = segments + .into_iter() + .map(|(s, e)| Segment { start: s, end: e }) + .collect(); + // The exemption is unconditional: it holds under exact matching too. + assert_eq!( + record.align_helper(&transcript, &TxAlignOpts::exact()), + expected + ); + assert_eq!( + record.align_helper(&transcript, &TxAlignOpts::default()), + expected + ); + }; + + // Reaching past the transcript start / end. + test_fun(vec![(80, 150)], AlignType::Exonic); + test_fun(vec![(1050, 1150)], AlignType::Exonic); + test_fun(vec![(80, 200), (500, 600), (1000, 1150)], AlignType::Exonic); + + // Extending into intron 1 from inside the first exon: the right edge + // of the first exon is a junction, checked exactly. + test_fun(vec![(120, 300)], AlignType::Spanning); + // Read starts in intron 1 and reaches left of exon 2: exon 2 is not + // the transcript's first exon, so this is genuine intron retention. + test_fun(vec![(450, 590)], AlignType::Spanning); + // Multi-segment read reaching left of a middle exon. + test_fun(vec![(450, 600), (1000, 1100)], AlignType::Spanning); + + // A single-exon transcript: both edges of its only exon are terminal. + let single = Transcript { + id: "tx2".to_string(), + chrom: "chr1".to_string(), + start: 100, + end: 200, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![Exon { + start: 100, + end: 200, + }]), + }; + record.segments = vec![Segment { + start: 80, + end: 220, + }]; + assert_eq!( + record.align_helper(&single, &TxAlignOpts::exact()), + AlignType::Exonic + ); + assert_eq!( + record.align_helper(&single, &TxAlignOpts::default()), + AlignType::Exonic + ); + } + + /// An intron short enough that aligners routinely miss the skip is an + /// alignment artifact, not retained-intron evidence. A segment that only + /// reaches partway into an intron holds genuine intronic sequence and is + /// never exempt. + #[test] + fn test_micro_intron() { + // 10bp intron. + let micro = Transcript { + id: "tx1".to_string(), + chrom: "chr1".to_string(), + start: 10, + end: 40, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { start: 10, end: 20 }, + Exon { start: 30, end: 40 }, + ]), + }; + // 40bp and 40bp introns: each within the default threshold, but their + // sum (80bp) exceeds it, so the length check must be per-gap. + let micro_chain = Transcript { + id: "tx2".to_string(), + chrom: "chr1".to_string(), + start: 10, + end: 120, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { start: 10, end: 20 }, + Exon { start: 60, end: 70 }, + Exon { + start: 110, + end: 120, + }, + ]), + }; + // 50bp intron, exactly at the default threshold. + let boundary = Transcript { + id: "tx3".to_string(), + chrom: "chr1".to_string(), + start: 10, + end: 80, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { start: 10, end: 20 }, + Exon { start: 70, end: 80 }, + ]), + }; + // 200bp intron: genuine retention evidence. + let real = Transcript { + id: "tx4".to_string(), + chrom: "chr1".to_string(), + start: 10, + end: 230, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { start: 10, end: 20 }, + Exon { + start: 220, + end: 230, + }, + ]), + }; + // Zero-length intron (abutting exons): holds no intronic sequence. + let abutting = Transcript { + id: "tx5".to_string(), + chrom: "chr1".to_string(), + start: 10, + end: 30, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { start: 10, end: 20 }, + Exon { start: 20, end: 30 }, + ]), + }; + + let mut record = SplicedRecord { + chrom: "chr1".to_string(), + segments: vec![], + strand: bed_utils::bed::Strand::Forward, + }; + let mut test_fun = |transcript: &Transcript, + segments: Vec<(u64, u64)>, + opts: &TxAlignOpts, + expected: AlignType| { + record.segments = segments + .into_iter() + .map(|(s, e)| Segment { start: s, end: e }) + .collect(); + assert_eq!(record.align_helper(transcript, opts), expected); + }; + let exact = TxAlignOpts::exact(); + let default = TxAlignOpts::default(); + + // Single segment through a 10bp intron: artifact at default, retained + // at exact. + test_fun(µ, vec![(12, 38)], &default, AlignType::Exonic); + test_fun(µ, vec![(12, 38)], &exact, AlignType::Spanning); + // Multi-segment read whose first segment bridges a 40bp micro-intron; + // the observed junction between exons 2 and 3 is not bridging evidence. + test_fun( + µ_chain, + vec![(12, 70), (110, 115)], + &default, + AlignType::Exonic, + ); + test_fun( + µ_chain, + vec![(12, 70), (110, 115)], + &exact, AlignType::Spanning, ); + // 50bp intron sits exactly at the default threshold (`<=`). + test_fun(&boundary, vec![(12, 75)], &default, AlignType::Exonic); + test_fun(&boundary, vec![(12, 75)], &exact, AlignType::Spanning); + // 200bp intron: never exempt. + test_fun(&real, vec![(12, 225)], &default, AlignType::Spanning); + test_fun(&real, vec![(12, 225)], &exact, AlignType::Spanning); + // Zero-length intron: exempt even under exact matching. + test_fun(&abutting, vec![(12, 28)], &default, AlignType::Exonic); + test_fun(&abutting, vec![(12, 28)], &exact, AlignType::Exonic); + // Reaching partway into a 10bp intron is genuine intronic sequence. + test_fun(µ, vec![(12, 25)], &default, AlignType::Spanning); + test_fun(µ, vec![(12, 25)], &exact, AlignType::Spanning); } } diff --git a/precellar-core/src/transcriptome/mod.rs b/precellar-core/src/transcriptome/mod.rs index 859b7c1..a685131 100644 --- a/precellar-core/src/transcriptome/mod.rs +++ b/precellar-core/src/transcriptome/mod.rs @@ -10,7 +10,7 @@ mod annotate; pub mod gtf; mod quantification; -pub use align::{TxAlignResult, TxAligner, TxAlignment}; +pub use align::{TxAlignOpts, TxAlignResult, TxAligner, TxAlignment}; pub use annotate::{GeneCounter, SplicedCounts}; pub use quantification::Quantifier; diff --git a/src/align.rs b/src/align.rs index 7d07bcc..1a9da64 100644 --- a/src/align.rs +++ b/src/align.rs @@ -475,11 +475,14 @@ impl AlignmentJobState { let mut aligner = AlignerRef::try_from(self.aligner.bind(py).clone())?; let header = aligner.header(); let transcript_annotator = if sink.needs_transcriptome() { + let opts = sink.tx_align_opts(); match sink.annotation_source() { - Some(source) => Some(source.load(header.clone(), strandedness)?), + Some(source) => Some(source.load(header.clone(), strandedness)?.with_opts(opts)), // No explicit source: STAR carries its own annotation, other // aligners have none and fail when the sink asks for it. - None => aligner.transcript_annotator(strandedness)?, + None => aligner + .transcript_annotator(strandedness)? + .map(|annotator| annotator.with_opts(opts)), } } else { None diff --git a/src/sinks.rs b/src/sinks.rs index 3ed19b5..c4741ce 100644 --- a/src/sinks.rs +++ b/src/sinks.rs @@ -6,7 +6,7 @@ use precellar::{ align::MultiMapR, fragment::{IntoFragOpts, IntoFragments}, qc::QcFragment, - transcriptome::Quantifier, + transcriptome::{Quantifier, TxAlignOpts}, }; use pyo3::prelude::*; use pyo3::types::PyDict; @@ -40,6 +40,12 @@ pub(crate) trait AlignmentSink { None } + /// Tolerance options for transcript structural matching, applied to the + /// transcript annotator before quantification. + fn tx_align_opts(&self) -> TxAlignOpts { + TxAlignOpts::default() + } + fn strandedness( &self, _assays: &[seqspec::Assay], @@ -295,6 +301,7 @@ pub(crate) struct GeneQuantificationSink { output: PathBuf, strandedness: Option, annotation_source: Option, + micro_intron_length: Option, } #[pymethods] @@ -319,13 +326,20 @@ impl GeneQuantificationSink { /// zstd-compressed. Works with any aligner. Its chromosome names must /// match the aligner index; mismatches fail before reads are processed. /// Mutually exclusive with `transcriptome_index`. + /// micro_intron_length : int, optional + /// Annotated introns at or below this length are treated as alignment + /// artifacts rather than retained introns: aligners routinely fail to + /// report short skips, so a read covering one looks like intron + /// retention when it is not. Set 0 to disable the exemption and + /// reproduce the historical strict behavior. Defaults to 50. #[new] - #[pyo3(signature = (output, *, strandedness=None, transcriptome_index=None, gtf=None))] + #[pyo3(signature = (output, *, strandedness=None, transcriptome_index=None, gtf=None, micro_intron_length=None))] fn new( output: PathBuf, strandedness: Option, transcriptome_index: Option, gtf: Option, + micro_intron_length: Option, ) -> Result { // Fail at construction time rather than after the alignment setup has // already run. @@ -344,6 +358,7 @@ impl GeneQuantificationSink { output, strandedness, annotation_source, + micro_intron_length, }) } } @@ -357,6 +372,18 @@ impl AlignmentSink for GeneQuantificationSink { self.annotation_source.clone() } + fn tx_align_opts(&self) -> TxAlignOpts { + // Delegate the default to `TxAlignOpts::default()` so that the Python + // binding never pins its own copy of the core default value. + match self.micro_intron_length { + Some(v) => TxAlignOpts { + micro_intron_length: v, + ..TxAlignOpts::default() + }, + None => TxAlignOpts::default(), + } + } + fn strandedness( &self, assays: &[seqspec::Assay], From aadea8c064cfefd02c98b666151596d7b25d53a8 Mon Sep 17 00:00:00 2001 From: Rui Xu Date: Mon, 12 Oct 2026 20:08:01 +0800 Subject: [PATCH 5/6] add junction tolerance --- precellar-core/src/transcriptome/align.rs | 164 ++++++++++++++++++++-- src/sinks.rs | 31 ++-- 2 files changed, 176 insertions(+), 19 deletions(-) diff --git a/precellar-core/src/transcriptome/align.rs b/precellar-core/src/transcriptome/align.rs index a680eb6..09cb2cf 100644 --- a/precellar-core/src/transcriptome/align.rs +++ b/precellar-core/src/transcriptome/align.rs @@ -64,7 +64,7 @@ pub struct TxAlignOpts { /// few bases, which must not be read as a structural mismatch. Only the /// inner edges of the alignment are affected; transcript ends are /// governed by the unconditional exemption above. IsoQuant uses 4 for - /// PacBio and 6 for ONT. Not yet implemented; reserved. + /// PacBio and 6 for ONT; `0` requires exact equality. pub junction_tolerance: u64, } @@ -87,6 +87,12 @@ impl TxAlignOpts { junction_tolerance: 0, } } + + /// Whether an observed junction coordinate matches an annotated exon + /// boundary within tolerance. + fn junction_matches(&self, observed: u64, annotated: u64) -> bool { + observed.abs_diff(annotated) <= self.junction_tolerance + } } #[derive(Eq, PartialEq, Debug, Clone, Encode, Decode)] @@ -269,9 +275,9 @@ impl SplicedRecord { /// Outer transcript edges are exempt from the containment checks: TSS/polyA /// annotations (also long-read end trimming) are unreliable, so a read overhanging /// a terminal exon says nothing about splicing. Inner edges (junctions) are still - /// checked exactly. The exemption is keyed on the transcript's terminal exons, not - /// the first overlapped exon, so that a read starting in intron 2 still counts as - /// intron retention rather than an overhang. + /// checked, within `opts.junction_tolerance`. The exemption is keyed on the + /// transcript's terminal exons, not the first overlapped exon, so that a read + /// starting in intron 2 still counts as intron retention rather than an overhang. fn align_helper(&self, transcript: &Transcript, opts: &TxAlignOpts) -> AlignType { if self.end() <= transcript.start || self.start() >= transcript.end { panic!("Read does not overlap transcript"); @@ -307,11 +313,14 @@ impl SplicedRecord { } else if segments[1..n_segment - 1] .iter() .zip(exons[1..n_ex - 1].iter()) - .all(|(s, e)| s.start == e.start && s.end == e.end) + .all(|(s, e)| { + opts.junction_matches(s.start, e.start) && opts.junction_matches(s.end, e.end) + }) { - // All internal segments (except the 1st and last) match exactly. - if segments[0].end == exons[0].end - && segments[n_segment - 1].start == exons[n_ex - 1].start + // All internal segments (except the 1st and last) match within + // tolerance. + if opts.junction_matches(segments[0].end, exons[0].end) + && opts.junction_matches(segments[n_segment - 1].start, exons[n_ex - 1].start) { // first and last segments also match. if (at_tx_start || segments[0].start >= exons[0].start) @@ -348,8 +357,8 @@ impl SplicedRecord { return false; }; - left_segment.end == left_exon.end - && right_segment.start == right_exon.start + opts.junction_matches(left_segment.end, left_exon.end) + && opts.junction_matches(right_segment.start, right_exon.start) && right_idx == left_idx + 1 }); @@ -1302,4 +1311,139 @@ mod tests { test_fun(µ, vec![(12, 25)], &default, AlignType::Spanning); test_fun(µ, vec![(12, 25)], &exact, AlignType::Spanning); } + + /// Splice sites carry a GT-AG motif, so a junction that misses the + /// annotated boundary is evidence against the isoform -- but long-read + /// indel noise shifts junctions by a few bases, and treating that as a + /// structural mismatch loses the read entirely (`Discordant` is filtered + /// out before counting). The tolerance applies only to inner edges. + #[test] + fn test_junction_tolerance() { + let transcript = Transcript { + id: "tx1".to_string(), + chrom: "chr1".to_string(), + start: 100, + end: 1100, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { + start: 100, + end: 200, + }, + Exon { + start: 500, + end: 600, + }, + Exon { + start: 1000, + end: 1100, + }, + ]), + }; + + let mut record = SplicedRecord { + chrom: "chr1".to_string(), + segments: vec![], + strand: bed_utils::bed::Strand::Forward, + }; + let mut test_fun = |segments: Vec<(u64, u64)>, tolerance: u64, expected: AlignType| { + record.segments = segments + .into_iter() + .map(|(s, e)| Segment { start: s, end: e }) + .collect(); + let opts = TxAlignOpts { + junction_tolerance: tolerance, + ..TxAlignOpts::exact() + }; + assert_eq!( + record.align_helper(&transcript, &opts), + expected, + "segments={:?} tolerance={}", + record.segments, + tolerance + ); + }; + + // An exact three-exon match is unaffected by tolerance. + for tol in [0, 4, 6] { + test_fun( + vec![(120, 200), (500, 600), (1000, 1080)], + tol, + AlignType::Exonic, + ); + } + + // Inner edge of the first segment off by 3bp: recovered at 4 and 6. + test_fun( + vec![(120, 197), (500, 600), (1000, 1080)], + 0, + AlignType::Discordant, + ); + test_fun( + vec![(120, 197), (500, 600), (1000, 1080)], + 4, + AlignType::Exonic, + ); + + // Off by 5bp: still discordant at 4, recovered at 6 (PacBio vs ONT). + test_fun( + vec![(120, 195), (500, 600), (1000, 1080)], + 4, + AlignType::Discordant, + ); + test_fun( + vec![(120, 195), (500, 600), (1000, 1080)], + 6, + AlignType::Exonic, + ); + + // Off by 7bp exceeds every supported tolerance. + for tol in [0, 4, 6] { + test_fun( + vec![(120, 193), (500, 600), (1000, 1080)], + tol, + AlignType::Discordant, + ); + } + + // Both edges of the internal segment off by 3bp. + test_fun( + vec![(120, 200), (503, 597), (1000, 1080)], + 0, + AlignType::Discordant, + ); + test_fun( + vec![(120, 200), (503, 597), (1000, 1080)], + 4, + AlignType::Exonic, + ); + + // Inner edge of the last segment off by 3bp. + test_fun( + vec![(120, 200), (500, 600), (1003, 1080)], + 0, + AlignType::Discordant, + ); + test_fun( + vec![(120, 200), (500, 600), (1003, 1080)], + 4, + AlignType::Exonic, + ); + + // Tolerance must not rescue a junction that no boundary explains: 300 + // is 100bp past the end of exon 1. + for tol in [0, 4, 6] { + test_fun(vec![(120, 300), (500, 600)], tol, AlignType::Discordant); + } + + // Outer edges are governed by the terminal exemption, not tolerance: + // overhanging the transcript start by 20bp stays exonic at tolerance 0. + test_fun( + vec![(80, 200), (500, 600), (1000, 1120)], + 0, + AlignType::Exonic, + ); + } } diff --git a/src/sinks.rs b/src/sinks.rs index c4741ce..2c16ad0 100644 --- a/src/sinks.rs +++ b/src/sinks.rs @@ -302,6 +302,7 @@ pub(crate) struct GeneQuantificationSink { strandedness: Option, annotation_source: Option, micro_intron_length: Option, + junction_tolerance: Option, } #[pymethods] @@ -332,14 +333,23 @@ impl GeneQuantificationSink { /// report short skips, so a read covering one looks like intron /// retention when it is not. Set 0 to disable the exemption and /// reproduce the historical strict behavior. Defaults to 50. + /// junction_tolerance : int, optional + /// Permitted deviation, in bases, between an observed splice junction + /// and the annotated exon boundary. Long-read indel noise routinely + /// shifts a junction by a few bases; without tolerance such a read is + /// discordant with every isoform and drops out of the counts + /// entirely. Set by sequencing chemistry: 4 for PacBio, 6 for ONT. + /// Only inner edges are affected -- transcript ends are always exempt. + /// Defaults to 0 (exact matching), which suits short reads. #[new] - #[pyo3(signature = (output, *, strandedness=None, transcriptome_index=None, gtf=None, micro_intron_length=None))] + #[pyo3(signature = (output, *, strandedness=None, transcriptome_index=None, gtf=None, micro_intron_length=None, junction_tolerance=None))] fn new( output: PathBuf, strandedness: Option, transcriptome_index: Option, gtf: Option, micro_intron_length: Option, + junction_tolerance: Option, ) -> Result { // Fail at construction time rather than after the alignment setup has // already run. @@ -359,6 +369,7 @@ impl GeneQuantificationSink { strandedness, annotation_source, micro_intron_length, + junction_tolerance, }) } } @@ -373,14 +384,16 @@ impl AlignmentSink for GeneQuantificationSink { } fn tx_align_opts(&self) -> TxAlignOpts { - // Delegate the default to `TxAlignOpts::default()` so that the Python - // binding never pins its own copy of the core default value. - match self.micro_intron_length { - Some(v) => TxAlignOpts { - micro_intron_length: v, - ..TxAlignOpts::default() - }, - None => TxAlignOpts::default(), + // Start from the core default so the Python binding never pins its own + // copy of a default value; override only what the caller supplied. + let defaults = TxAlignOpts::default(); + TxAlignOpts { + micro_intron_length: self + .micro_intron_length + .unwrap_or(defaults.micro_intron_length), + junction_tolerance: self + .junction_tolerance + .unwrap_or(defaults.junction_tolerance), } } From 5a378324aae889b6ec01a781167a17a8cf688d1b Mon Sep 17 00:00:00 2001 From: Rui Xu Date: Tue, 13 Oct 2026 16:08:39 +0800 Subject: [PATCH 6/6] minor fix --- precellar-core/src/transcriptome/align.rs | 60 ++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/precellar-core/src/transcriptome/align.rs b/precellar-core/src/transcriptome/align.rs index 09cb2cf..01385b5 100644 --- a/precellar-core/src/transcriptome/align.rs +++ b/precellar-core/src/transcriptome/align.rs @@ -385,7 +385,14 @@ impl SplicedRecord { .any(|segment| segment.start <= left.end && segment.end >= right.start) }) .all(|pair| pair[1].start - pair[0].end <= opts.micro_intron_length); - if only_artifacts { + // Bridged introns being artifacts is not enough: the outer + // edges must also be contained, exactly as in the branch above. + // A segment running past the last overlapped exon without + // reaching the next one holds retained-intron sequence that no + // gap test above can see, since that intron is never bridged. + let ends_contained = (at_tx_start || segments[0].start >= exons[0].start) + && (at_tx_end || segments[n_segment - 1].end <= exons[n_ex - 1].end); + if only_artifacts && ends_contained { AlignType::Exonic } else { AlignType::Spanning @@ -1446,4 +1453,55 @@ mod tests { AlignType::Exonic, ); } + + /// The micro-intron exemption must not swallow a segment that runs past + /// the last exon it overlaps without reaching the next one. That overhang + /// is retained-intron sequence, and the per-gap length test cannot see it: + /// the intron it sits in is never bridged, so it is never a candidate gap. + #[test] + fn test_micro_intron_requires_contained_ends() { + let transcript = Transcript { + id: "tx1".to_string(), + chrom: "chr1".to_string(), + start: 10, + end: 1010, + strand: bed_utils::bed::Strand::Forward, + gene_id: "gene1".to_string(), + gene_name: "GeneName1".to_string(), + exons: Exons(vec![ + Exon { start: 10, end: 20 }, + Exon { start: 30, end: 40 }, + Exon { + start: 1000, + end: 1010, + }, + ]), + }; + + let mut record = SplicedRecord { + chrom: "chr1".to_string(), + segments: vec![], + strand: bed_utils::bed::Strand::Forward, + }; + let mut test_fun = |segments: Vec<(u64, u64)>, expected: AlignType| { + record.segments = segments + .into_iter() + .map(|(s, e)| Segment { start: s, end: e }) + .collect(); + assert_eq!( + record.align_helper(&transcript, &TxAlignOpts::default()), + expected, + "segments={:?}", + record.segments + ); + }; + + // Bridges the 10bp intron and stops inside exon 2: a missed skip. + test_fun(vec![(12, 38)], AlignType::Exonic); + // Same bridge, but then runs 460bp into the intron before exon 3. + test_fun(vec![(12, 500)], AlignType::Spanning); + // Reaching exon 3 without a junction bridges the 960bp intron too, + // which the gap test already rejects. + test_fun(vec![(12, 1005)], AlignType::Spanning); + } }