fix: validate inputs and fix coordinate offsets in pattern matching

Introduce explicit input validation and a `(-1, -1, -1)` sentinel return value to signal invalid or unreliable matches. Add pre-call length checks and early returns to prevent panics during indel relocation. Fix an index offset bug in coordinate calculation by preserving the original fragment position before applying alignment deltas. Update match filtering to enforce reliability checks only on successfully relocated alignments. Comprehensive tests validate edge cases, boundary constraints, and error thresholds.
2026-08-24 05:41:19 +00:00 · 2026-08-19 17:25:44 +02:00
parent 438893d910
commit 620646d412
3 changed files with 174 additions and 27 deletions
@@ -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)]
 }
@@ -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)
+			}
+		}()
+	}
+}
@@ -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)

+			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++
 		}