Skip to content

Commit e09748b

Browse files
pinin4fjordsclaude
andcommitted
fix(stitch): use sjdb_score in place of motif_score on annotated junctions
STAR's Transcript_alignScore.cpp adds either sjdbScore (annotated) or motifScore (unannotated) to the alignment score, not both. rustar was adding both on annotated junctions, over-crediting the score by motif_score per junction. Currently invisible because #27 keeps is_annotated false in practice; becomes the dominant remaining AS divergence on annotated splices as soon as #27 lands. Fixes #50 Co-Authored-By: Claude <noreply@anthropic.com>
1 parent 70be24d commit e09748b

1 file changed

Lines changed: 37 additions & 1 deletion

File tree

src/align/stitch.rs

Lines changed: 37 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1313,9 +1313,10 @@ fn stitch_align_to_transcript(
13131313
|| db.is_annotated(cluster.chr_idx, donor_fwd, acceptor_fwd, 2)
13141314
});
13151315

1316-
d_score += motif_score;
13171316
if is_annotated {
13181317
d_score += scorer.sjdb_score;
1318+
} else {
1319+
d_score += motif_score;
13191320
}
13201321

13211322
new_wt.n_junction += 1;
@@ -3430,4 +3431,39 @@ mod tests {
34303431
bin_start_edge = bin_start_edge.saturating_sub(win_flank_nbins);
34313432
assert_eq!(bin_start_edge, 0, "Flanking should saturate at 0");
34323433
}
3434+
3435+
#[test]
3436+
fn test_junction_score_annotated_uses_sjdb_not_motif() {
3437+
use crate::align::score::AlignmentScorer;
3438+
3439+
let scorer = AlignmentScorer::from_params_minimal();
3440+
3441+
for motif_score in [0_i32, scorer.score_gap_gcag, scorer.score_gap_atac] {
3442+
let baseline: i32 = 100;
3443+
3444+
let mut d_score_annot = baseline;
3445+
let is_annotated = true;
3446+
if is_annotated {
3447+
d_score_annot += scorer.sjdb_score;
3448+
} else {
3449+
d_score_annot += motif_score;
3450+
}
3451+
assert_eq!(d_score_annot, baseline + scorer.sjdb_score);
3452+
3453+
let mut d_score_novel = baseline;
3454+
let is_annotated = false;
3455+
if is_annotated {
3456+
d_score_novel += scorer.sjdb_score;
3457+
} else {
3458+
d_score_novel += motif_score;
3459+
}
3460+
assert_eq!(d_score_novel, baseline + motif_score);
3461+
}
3462+
3463+
let baseline: i32 = 100;
3464+
let motif_score = scorer.score_gap_gcag;
3465+
let additive_buggy = baseline + motif_score + scorer.sjdb_score;
3466+
let replacement_correct = baseline + scorer.sjdb_score;
3467+
assert_ne!(additive_buggy, replacement_correct);
3468+
}
34333469
}

0 commit comments

Comments
 (0)