Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/align/read_align.rs
Original file line number Diff line number Diff line change
Expand Up @@ -311,7 +311,7 @@ pub fn align_read(
junction_db,
params.align_transcripts_per_window_nmax,
debug_name,
)?;
);
if debug_read {
eprintln!(
"[DEBUG {}] Cluster[{}]: {} transcripts from DP",
Expand Down Expand Up @@ -814,7 +814,7 @@ pub fn align_paired_read(
params.align_transcripts_per_window_nmax,
params.align_mates_gap_max.into(),
debug_name,
)?;
);

for wt in &wts {
let split_result =
Expand Down
24 changes: 12 additions & 12 deletions src/align/score.rs
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,7 @@ impl AlignmentScorer {
genome_pos
};
let motif = self.detect_splice_motif(donor, len, genome);
let score = self.score_splice_junction(&motif);
let score = self.score_splice_junction(motif);
(
score,
GapType::SpliceJunction {
Expand Down Expand Up @@ -235,7 +235,7 @@ impl AlignmentScorer {
rc_donor
};
let motif = self.detect_splice_motif(donor, del_len, genome);
let score = self.score_splice_junction(&motif);
let score = self.score_splice_junction(motif);
(
score,
GapType::SpliceJunction {
Expand Down Expand Up @@ -389,7 +389,7 @@ impl AlignmentScorer {
donor_sa
};
let motif = self.detect_splice_motif(donor_fwd, del as u32, genome);
let motif_score = self.score_splice_junction(&motif);
let motif_score = self.score_splice_junction(motif);
let score2 = score1 + motif_score;

if score2 > max_score2 {
Expand Down Expand Up @@ -493,7 +493,7 @@ impl AlignmentScorer {
donor_sa
};
best_motif = self.detect_splice_motif(donor_fwd, del as u32, genome);
best_motif_score = self.score_splice_junction(&best_motif);
best_motif_score = self.score_splice_junction(best_motif);
}
}

Expand All @@ -518,7 +518,7 @@ impl AlignmentScorer {
}

/// Score a splice junction based on motif
pub(crate) fn score_splice_junction(&self, motif: &SpliceMotif) -> i32 {
pub(crate) fn score_splice_junction(&self, motif: SpliceMotif) -> i32 {
match motif {
SpliceMotif::GtAg | SpliceMotif::CtAc => self.score_gap,
SpliceMotif::GcAg | SpliceMotif::CtGc => self.score_gap_gcag,
Expand Down Expand Up @@ -695,7 +695,7 @@ mod tests {
let motif = scorer.detect_splice_motif(2, 12, &genome);
assert_eq!(motif, SpliceMotif::GtAg);

let score = scorer.score_splice_junction(&motif);
let score = scorer.score_splice_junction(motif);
assert_eq!(score, 0); // Canonical
}

Expand Down Expand Up @@ -738,7 +738,7 @@ mod tests {
let motif = scorer.detect_splice_motif(2, 12, &genome);
assert_eq!(motif, SpliceMotif::GcAg);

let score = scorer.score_splice_junction(&motif);
let score = scorer.score_splice_junction(motif);
assert_eq!(score, -4);
}

Expand Down Expand Up @@ -781,7 +781,7 @@ mod tests {
let motif = scorer.detect_splice_motif(2, 12, &genome);
assert_eq!(motif, SpliceMotif::AtAc);

let score = scorer.score_splice_junction(&motif);
let score = scorer.score_splice_junction(motif);
assert_eq!(score, -8);
}

Expand Down Expand Up @@ -822,7 +822,7 @@ mod tests {
let motif = scorer.detect_splice_motif(2, 12, &genome);
assert_eq!(motif, SpliceMotif::NonCanonical);

let score = scorer.score_splice_junction(&motif);
let score = scorer.score_splice_junction(motif);
assert_eq!(score, -8);
}

Expand Down Expand Up @@ -1017,7 +1017,7 @@ mod tests {
let motif = scorer.detect_splice_motif(2, 12, &genome_ctac);
assert_eq!(motif, SpliceMotif::CtAc);
// Should score same as canonical GT-AG
assert_eq!(scorer.score_splice_junction(&motif), 0);
assert_eq!(scorer.score_splice_junction(motif), 0);

// CT-GC motif: (1,3,2,1) — reverse complement of GC-AG
let seq_ctgc = vec![
Expand All @@ -1030,7 +1030,7 @@ mod tests {
let genome_ctgc = make_test_genome(&seq_ctgc);
let motif = scorer.detect_splice_motif(2, 12, &genome_ctgc);
assert_eq!(motif, SpliceMotif::CtGc);
assert_eq!(scorer.score_splice_junction(&motif), -4);
assert_eq!(scorer.score_splice_junction(motif), -4);

// GT-AT motif: (2,3,0,3) — reverse complement of AT-AC
let seq_gtat = vec![
Expand All @@ -1043,7 +1043,7 @@ mod tests {
let genome_gtat = make_test_genome(&seq_gtat);
let motif = scorer.detect_splice_motif(2, 12, &genome_gtat);
assert_eq!(motif, SpliceMotif::GtAt);
assert_eq!(scorer.score_splice_junction(&motif), -8);
assert_eq!(scorer.score_splice_junction(motif), -8);
}

#[test]
Expand Down
51 changes: 24 additions & 27 deletions src/align/seed.rs
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ impl Seed {
false,
debug_name,
&mut seeds,
)?;
);

// Cap check between directions (STAR: seedPerReadNmax applies across both)
if seeds.len() >= params.seed_per_read_nmax {
Expand All @@ -87,7 +87,7 @@ impl Seed {
true,
debug_name,
&mut seeds,
)?;
);

// STAR's storeAligns dedup: same rStart + same Length → skip duplicate.
// Multiple istart chains can find the same (read_pos, length, direction) seed.
Expand Down Expand Up @@ -217,7 +217,7 @@ fn search_direction_sparse(
is_rc: bool,
debug_name: &str,
seeds: &mut Vec<Seed>,
) -> Result<(), Error> {
) {
let read_len = read_seq.len();

// STAR (ReadAlign_mapOneRead.cpp lines 41-42):
Expand Down Expand Up @@ -262,7 +262,7 @@ fn search_direction_sparse(
}

let result =
find_seed_at_position(read_seq, pos, index, min_seed_length, false, params)?;
find_seed_at_position(read_seq, pos, index, min_seed_length, false, params);

if !debug_name.is_empty() {
let dir = if is_rc { "RC" } else { "FWD" };
Expand Down Expand Up @@ -292,16 +292,14 @@ fn search_direction_sparse(
seeds.push(seed);

if seeds.len() >= params.seed_per_read_nmax {
return Ok(());
return;
}
}

pos += result.advance; // Always advance by MMP length (matches STAR)
// Remaining-length check at loop top: stop when < seedMapMin bases remain
}
}

Ok(())
}

/// Find a seed starting at a specific position in the read.
Expand All @@ -317,23 +315,23 @@ fn find_seed_at_position(
min_seed_length: usize,
is_reverse: bool,
params: &Parameters,
) -> Result<MmpResult, Error> {
) -> MmpResult {
if read_pos >= read_seq.len() {
return Ok(MmpResult {
return MmpResult {
seed: None,
advance: 1,
});
};
}

// Extract k-mer for SAindex lookup
let sa_nbases = index.sa_index.nbases as usize;
let remaining = read_seq.len() - read_pos;

if remaining < min_seed_length {
return Ok(MmpResult {
return MmpResult {
seed: None,
advance: 1,
});
};
}

// Build k-mer for SAindex lookup, stopping at first N base
Expand All @@ -351,10 +349,10 @@ fn find_seed_at_position(
}

if actual_len == 0 {
return Ok(MmpResult {
return MmpResult {
seed: None,
advance: 1,
}); // First base is N
}; // First base is N
}

// Hierarchical SAindex lookup (STAR's maxMappableLength2strands approach).
Expand All @@ -367,17 +365,17 @@ fn find_seed_at_position(
.hierarchical_lookup(kmer_idx, actual_len as u32, n_sa);

let Some((sa_start, sa_end, matched_level, bounds_tight)) = result else {
return Ok(MmpResult {
return MmpResult {
seed: None,
advance: 1,
});
};
};

if sa_start >= sa_end {
return Ok(MmpResult {
return MmpResult {
seed: None,
advance: 1,
});
};
}

// STAR short-circuit (maxMappableLength2strands.cpp):
Expand Down Expand Up @@ -407,14 +405,14 @@ fn find_seed_at_position(
// Uses narrowed range (accurate loci count, not overestimated k-mer range)
let n_loci = narrowed_end - narrowed_start;
if n_loci > params.seed_multimap_nmax {
return Ok(MmpResult {
return MmpResult {
seed: None,
advance,
});
};
}

if match_length >= min_seed_length {
Ok(MmpResult {
MmpResult {
seed: Some(Seed {
read_pos,
length: match_length,
Expand All @@ -425,12 +423,12 @@ fn find_seed_at_position(
mate_id: 2, // Single-end default
}),
advance,
})
}
} else {
Ok(MmpResult {
MmpResult {
seed: None,
advance,
})
}
}
}

Expand Down Expand Up @@ -1023,16 +1021,15 @@ mod tests {
// Count how many seeds dense would produce (every position that has a match)
let mut dense_count = 0;
for read_pos in 0..read.len() {
let result = find_seed_at_position(&read, read_pos, &index, 4, false, &params).unwrap();
let result = find_seed_at_position(&read, read_pos, &index, 4, false, &params);
if result.seed.is_some() {
dense_count += 1;
}
}
// Also count R→L dense seeds
let rc_read = reverse_complement_read(&read);
for rc_pos in 0..rc_read.len() {
let result =
find_seed_at_position(&rc_read, rc_pos, &index, 4, false, &params).unwrap();
let result = find_seed_at_position(&rc_read, rc_pos, &index, 4, false, &params);
if result.seed.is_some() {
dense_count += 1;
}
Expand Down
21 changes: 10 additions & 11 deletions src/align/stitch.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,6 @@
use crate::align::score::AlignmentScorer;
use crate::align::seed::Seed;
use crate::align::transcript::{CigarOpExt as _, Transcript};
use crate::error::Error;
use crate::index::GenomeIndex;
use noodles::sam::alignment::record::cigar;

Expand Down Expand Up @@ -1500,7 +1499,7 @@ pub fn stitch_seeds(
read_seq: &[u8],
index: &GenomeIndex,
scorer: &AlignmentScorer,
) -> Result<Vec<Transcript>, Error> {
) -> Vec<Transcript> {
stitch_seeds_with_jdb(cluster, read_seq, index, scorer, None, 1)
}

Expand All @@ -1519,7 +1518,7 @@ pub fn stitch_seeds_with_jdb(
scorer: &AlignmentScorer,
junction_db: Option<&crate::junction::SpliceJunctionDb>,
max_transcripts_per_window: usize,
) -> Result<Vec<Transcript>, Error> {
) -> Vec<Transcript> {
stitch_seeds_with_jdb_debug(
cluster,
read_seq,
Expand Down Expand Up @@ -2469,7 +2468,7 @@ pub(crate) fn stitch_seeds_with_jdb_debug(
junction_db: Option<&crate::junction::SpliceJunctionDb>,
max_transcripts_per_window: usize,
debug_read_name: &str,
) -> Result<Vec<Transcript>, Error> {
) -> Vec<Transcript> {
let (working_transcripts, stitch_cluster, stitch_is_reverse, stitch_read) = stitch_seeds_core(
cluster,
read_seq,
Expand All @@ -2479,7 +2478,7 @@ pub(crate) fn stitch_seeds_with_jdb_debug(
max_transcripts_per_window,
0,
debug_read_name,
)?;
);

// Finalize working transcripts → Transcript (filtering by overhang+repeat check)
let mut transcripts: Vec<Transcript> = Vec::with_capacity(working_transcripts.len());
Expand Down Expand Up @@ -2583,7 +2582,7 @@ pub(crate) fn stitch_seeds_with_jdb_debug(
});
transcripts.truncate(max_transcripts_per_window);

Ok(transcripts)
transcripts
}

/// Shared core: preprocessing + recursive stitcher, returns working transcripts + context.
Expand All @@ -2597,7 +2596,7 @@ pub(crate) fn stitch_seeds_core(
max_transcripts_per_window: usize,
align_mates_gap_max: u64,
debug_read_name: &str,
) -> Result<(Vec<WorkingTranscript>, SeedCluster, bool, Vec<u8>), Error> {
) -> (Vec<WorkingTranscript>, SeedCluster, bool, Vec<u8>) {
let debug = !debug_read_name.is_empty();

// Include ALL seeds (anchor and non-anchor) in the stitcher.
Expand Down Expand Up @@ -2730,12 +2729,12 @@ pub(crate) fn stitch_seeds_core(
}

if wa_entries.is_empty() {
return Ok((
return (
Vec::new(),
stitch_cluster,
stitch_is_reverse,
stitch_read_owned,
));
);
}

if debug {
Expand Down Expand Up @@ -2924,12 +2923,12 @@ pub(crate) fn stitch_seeds_core(
);
}

Ok((
(
working_transcripts,
stitch_cluster,
stitch_is_reverse,
stitch_read_owned,
))
)
}

#[cfg(test)]
Expand Down
Loading
Loading