diff --git a/pkg/obialign/locatepattern.go b/pkg/obialign/locatepattern.go index d5170e6..304cfbe 100644 --- a/pkg/obialign/locatepattern.go +++ b/pkg/obialign/locatepattern.go @@ -28,10 +28,18 @@ func buffIndex(i, j, width int) int { // // The function returns the start and end positions of the best // match, as well as the number of errors in the best match. +// +// When the sequence is too short relative to the pattern for the +// backtracking to reconstruct a valid alignment (e.g. the pattern +// is longer than the sequence, or the match sits too close to a +// sequence boundary), no reliable position can be computed. In that +// case the function returns the sentinel (-1, -1, -1) instead of a +// guessed, potentially wrong, position: callers must treat this as +// "no match" rather than use the returned coordinates. func LocatePattern(id string, pattern, sequence []byte) (int, int, int) { - if len(pattern) >= len(sequence) { - log.Panicf("Sequence %s:Pattern %s must be shorter than sequence %s", id, pattern, sequence) + if len(sequence) == 0 { + log.Panicf("Sequence %s:Pattern %s must not be empty", id, pattern) } // Pattern spreads over the columns @@ -158,5 +166,18 @@ func LocatePattern(id string, pattern, sequence []byte) (int, int, int) { // obilog.Warnf("from : %d to: %d error: %d match: %v", // i, end+1, -buffer[buffIndex(len(sequence)-1, len(pattern)-1, width)], // string(sequence[i:(end+1)])) + + if i < 0 || end == -1 { + // i < 0: the backtracking ran off the start of the sequence + // without fully consuming the pattern. + // end == -1: the backtracking loop never ran at all (e.g. a + // single-base pattern, jmax == 0), so no alignment boundary + // was ever established. + // Either way, no valid alignment exists for this (pattern, + // sequence) pair: signal it explicitly instead of returning + // an out-of-bounds or uncomputed position. + return -1, -1, -1 + } + return i, end + 1, -buffer[buffIndex(len(sequence)-1, len(pattern)-1, width)] } diff --git a/pkg/obialign/locatepattern_test.go b/pkg/obialign/locatepattern_test.go new file mode 100644 index 0000000..70fcbfe --- /dev/null +++ b/pkg/obialign/locatepattern_test.go @@ -0,0 +1,123 @@ +package obialign + +import ( + "math/rand" + "testing" +) + +func TestLocatePatternNormal(t *testing.T) { + // Pattern fully and exactly present in the middle of a longer sequence. + start, end, nerr := LocatePattern("id", []byte("ACGT"), []byte("TTTTACGTTTTT")) + if start != 4 || end != 8 || nerr != 0 { + t.Errorf("got start=%d end=%d nerr=%d, want start=4 end=8 nerr=0", start, end, nerr) + } +} + +func TestLocatePatternOneMismatch(t *testing.T) { + start, end, nerr := LocatePattern("id", []byte("ACGT"), []byte("TTTTACTTTTTT")) + if nerr != 1 { + t.Errorf("got nerr=%d, want 1 (start=%d end=%d)", nerr, start, end) + } +} + +// The real-world case that used to panic: pattern longer than the sequence +// fragment extracted for indel relocation. +func TestLocatePatternPatternLongerThanSequence(t *testing.T) { + start, end, nerr := LocatePattern("id", + []byte("GGGCAATCCTGAGCCAAATC"), + []byte("tcctgagccaaatcacgtt")) + + if start != -1 || end != -1 || nerr != -1 { + t.Errorf("got start=%d end=%d nerr=%d, want the (-1,-1,-1) sentinel", start, end, nerr) + } +} + +func TestLocatePatternSequenceLengthOne(t *testing.T) { + start, end, nerr := LocatePattern("id", []byte("AB"), []byte("A")) + + if start < 0 || end < 0 || nerr < 0 { + t.Fatalf("got start=%d end=%d nerr=%d, expected a valid (non-sentinel) result", start, end, nerr) + } + if start != 0 || end != 1 || nerr != 1 { + t.Errorf("got start=%d end=%d nerr=%d, want start=0 end=1 nerr=1", start, end, nerr) + } +} + +// A pattern much longer than the sequence can still yield a mathematically +// valid (in-bounds) alignment: the extra pattern length is absorbed as gaps, +// driving the error count high enough that the caller's maxerr threshold +// rejects it. The function itself must still return consistent bounds. +func TestLocatePatternPatternMuchLongerThanSequence(t *testing.T) { + start, end, nerr := LocatePattern("id", []byte("ACGTACGTACGTACGTACGT"), []byte("ACG")) + + isSentinel := start == -1 && end == -1 && nerr == -1 + isValid := start >= 0 && end > start && end <= 3 && nerr >= 0 + + if !isSentinel && !isValid { + t.Errorf("got start=%d end=%d nerr=%d, want either the sentinel or consistent in-bounds values", start, end, nerr) + } +} + +func TestLocatePatternNeverReturnsOutOfBounds(t *testing.T) { + patterns := []string{"A", "AC", "ACG", "ACGT", "ACGTA", "ACGTAC", "ACGTACG", "ACGTACGT"} + sequences := []string{"A", "AC", "ACG", "ACGT", "ACGTA", "ACGTAC", "ACGTACG", "ACGTACGT"} + + for _, p := range patterns { + for _, s := range sequences { + start, end, nerr := LocatePattern("id", []byte(p), []byte(s)) + + if start == -1 && end == -1 && nerr == -1 { + // Explicit "no reliable match" sentinel: always acceptable. + continue + } + + if start < 0 || end < 0 || start >= end || end > len(s) || nerr < 0 { + t.Errorf("pattern=%q sequence=%q -> start=%d end=%d nerr=%d is out of bounds / inconsistent", + p, s, start, end, nerr) + } + } + } +} + +// Randomized property test over a wide range of pattern/sequence length +// combinations, including pattern >= sequence, to make sure the function +// never panics and never returns anything but the sentinel or fully +// consistent, in-bounds coordinates. +func TestLocatePatternRandomizedNeverInvalid(t *testing.T) { + const bases = "ACGT" + rng := rand.New(rand.NewSource(42)) + + randSeq := func(n int) []byte { + b := make([]byte, n) + for i := range b { + b[i] = bases[rng.Intn(len(bases))] + } + return b + } + + for trial := 0; trial < 5000; trial++ { + patLen := 1 + rng.Intn(15) + seqLen := 1 + rng.Intn(15) + + pattern := randSeq(patLen) + sequence := randSeq(seqLen) + + func() { + defer func() { + if r := recover(); r != nil { + t.Fatalf("panic for pattern=%q sequence=%q: %v", pattern, sequence, r) + } + }() + + start, end, nerr := LocatePattern("id", pattern, sequence) + + isSentinel := start == -1 && end == -1 && nerr == -1 + isValid := start >= 0 && end > start && end <= len(sequence) && nerr >= 0 + + if !isSentinel && !isValid { + t.Errorf("pattern=%q sequence=%q -> start=%d end=%d nerr=%d is neither the sentinel nor consistent", + pattern, sequence, start, end, nerr) + } + }() + } +} diff --git a/pkg/obiapat/pattern.go b/pkg/obiapat/pattern.go index 25a3c7b..b58473f 100755 --- a/pkg/obiapat/pattern.go +++ b/pkg/obiapat/pattern.go @@ -373,17 +373,7 @@ func (pattern ApatPattern) BestMatch(sequence ApatSequence, begin, length int) ( cpattern := (*[1 << 30]byte)(unsafe.Pointer(pattern.pointer.pointer.cpat)) frg := sequence.pointer.reference.Sequence()[start:end] - - if len(frg) <= int(pattern.pointer.pointer.patlen) { - // The sequence is too short (e.g. the match is near one of its - // ends) to extract a fragment longer than the pattern, which - // obialign.LocatePattern requires. Keep the original match - // instead of refining it. - start = best[0] - end = best[1] - log.Debugln("Fragment too short for indel relocation, keeping original match", start, end, nerr) - return - } + fragStart := start log.Debugln( string(frg), @@ -395,11 +385,20 @@ func (pattern ApatPattern) BestMatch(sequence ApatSequence, begin, length int) ( (*cpattern)[0:int(pattern.pointer.pointer.patlen)], frg) - // olderr := m[2] + if from < 0 { + // obialign.LocatePattern could not reconstruct a reliable + // alignment (e.g. the fragment is too short relative to the + // pattern). Reporting a guessed position would risk placing + // the primer boundary incorrectly, so treat it as no match + // at all rather than falling back to an unrefined position. + matched = false + log.Debugln("No reliable indel relocation, discarding match", sequence.pointer.reference.Id()) + return + } nerr = score - start = start + from - end = start + to + start = fragStart + from + end = fragStart + to log.Debugf("BestMatch on %s : score=%d [%d..%d]", sequence.pointer.reference.Id(), score, start, nerr) return } @@ -478,6 +477,7 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int) for _, m := range res { // Recompute the start and end position of the match // when the pattern allows for indels + valid := true if m[2] > 0 && pattern.pointer.pointer.hasIndel { // obilog.Warnf("Locating indel on sequence %s[%s]", sequence.pointer.reference.Id(), pattern.String()) start := m[0] - m[2]*2 @@ -491,17 +491,20 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int) cpattern := (*[1 << 30]byte)(unsafe.Pointer(pattern.pointer.pointer.cpat)) frg := sequence.pointer.reference.Sequence()[start:end] - // obialign.LocatePattern requires the fragment to be strictly - // longer than the pattern. When the match sits near one of the - // sequence ends, the fragment can be clamped to the sequence - // boundaries and end up too short; in that case keep the - // original (unrefined) match instead of crashing. - if len(frg) > int(pattern.pointer.pointer.patlen) { - pb, pe, score := obialign.LocatePattern( - sequence.pointer.reference.Id(), - (*cpattern)[0:int(pattern.pointer.pointer.patlen)], - frg) + pb, pe, score := obialign.LocatePattern( + sequence.pointer.reference.Id(), + (*cpattern)[0:int(pattern.pointer.pointer.patlen)], + frg) + if pb < 0 { + // obialign.LocatePattern could not reconstruct a + // reliable alignment (e.g. the match sits too close + // to a sequence end for the fragment to be usable). + // Reporting a guessed position risks placing the + // primer boundary incorrectly, so drop the match + // entirely instead of keeping an unrefined guess. + valid = false + } else { // olderr := m[2] m[2] = score m[0] = start + pb @@ -512,7 +515,7 @@ func (pattern ApatPattern) AllMatches(sequence ApatSequence, begin, length int) } } - if int(pattern.pointer.pointer.maxerr) >= m[2] { + if valid && int(pattern.pointer.pointer.maxerr) >= m[2] { res[j] = m j++ }