Repository: marcelm/cutadapt

Issue:
Penalize indels more
Indels are enabled by default, but cost the same as mismatches. This often leads to unexpected results, typically longer sequences being trimmed than expected, see #597.

We could penalize indels more, but should also re-consider whether optimizing the number of matches is the best thing to do.

Discussion:
luchaoqi: Any updates on this PR? Love the way it penalizes different cases:
```
        self._match_score = 1
        self._mismatch_score = -1
        self._insertion_score = -2
        self._deletion_score = -2
```
I found CI fails and has not merged successfully. Is this because cutadapt has dropped support for python 3.6?
marcelm: The algorithm isn’t quite finished. See the To Do above: It’s not guaranteed that the leftmost match is found if there are multiple occurrences of the adapter.

Please have some patience. I will continue work on this when I have time.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 986a674..0a14734 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -6,6 +6,11 @@ Changelog
 development version
 -------------------
 
+* :issue:`604`, :pr:`608`: The `alignment algorithm was tweaked <algorithm-indel-scores>`_
+  to penalize indels more and to more accurately pick the leftmost adapter
+  occurrence if there are multiple. This will normally affect very few
+  reads, but should generally lead to fewer surprising results in cases
+  where it matters.
 * :issue:`607`: Print an error when an output file was specified
   multiple times (for example, for ``--untrimmed-output`` and
   ``--too-short-output``). Sending output from different filters to
diff --git a/doc/algorithms.rst b/doc/algorithms.rst
index d835367..dc85107 100644
--- a/doc/algorithms.rst
+++ b/doc/algorithms.rst
@@ -37,7 +37,7 @@ has the disadvantage that they are not at all intuitive: What does a total score
 of *x* mean? Is that good or bad? How should a threshold be chosen in order to
 avoid finding alignments with too many errors?
 
-For Cutadapt, the adapter alignment algorithm uses *unit costs* instead.
+For Cutadapt, the adapter alignment algorithm primarily uses *unit costs* instead.
 This means that mismatches, insertions and deletions are counted as one error, which
 is easier to understand and allows to specify a single parameter for the
 algorithm (the maximum error rate) in order to describe how many errors are
@@ -60,7 +60,8 @@ Conceptually, the procedure is as follows:
 2. Keep only those alignments that do not exceed the specified maximum error
    rate.
 3. Then, keep only those alignments that have a maximal number of matches
-   (that is, there is no alignment with more matches).
+   (that is, there is no alignment with more matches). (Note: This has been
+   changed, see the section below for an update.)
 4. If there are multiple alignments with the same number of matches, then keep
    only those that have the smallest error rate.
 5. If there are still multiple candidates left, choose the alignment that starts
@@ -70,6 +71,89 @@ In Step 1, the different adapter types are taken into account: Only those
 overlaps that are actually allowed by the adapter type are actually considered.
 
 
+.. _algorithm-indel-scores:
+
+Alignment algorithm changes in Cutadapt 4
+=========================================
+
+The above algorithm has been tweaked slightly in Cutadapt 4.
+The main problem was that the idea of maximizing the number of matches
+(criterion 3 in the section above) sometimes leads to unintuitive results.
+
+For example, the previous algorithm would prefer an alignment such as this one::
+
+    CCAGTCCTTTCCTGAGAGT   Read
+    ||||||||   ||
+    CCAGTCCT---CT         5' adapter
+
+This alignment was considered to be the best one because it contains 10 matches,
+which is the maximum possible.
+The three consecutive deletions are ignored when making that decision.
+To the user, the unexpected result is visible because the read would end up as
+``GAGAGT`` after trimming.
+
+With the tuned algorithm, the alignment is more sensible::
+
+    CCAGTCCTTTCCTGAGAGT   Read
+    ||||||||X|
+    CCAGTCCTCT            5' adapter
+
+The trimmed read is now ``CCTGAGAGT``, which is what one would likely expect.
+
+The alignment algorithm in Cutadapt can perhaps now be described as
+a *hybrid* algorithm that uses both edit distance and score:
+
+- Edit distance is used to fill out the dynamic programming matrix.
+  Conceptually, this can be seen as computing the edit distance for all
+  possible overlaps between the read and the adapter.
+  We need to use the edit distance as optimization criterion at this
+  stage because we want to be able to let the user provide a maximum
+  error rate (``-e``).
+  Also, using edit distance (that is, unit costs) allows using some
+  optimizations while filling in the matrix (Ukkonen’s trick).
+- A second matrix with scores is filled in simultaneously.
+  The value in a cell is the score of the edit-distance-based alignment,
+  the score is not used as optimization criterion.
+- Finally, the score is used to decide which of the overlaps between
+  read and adapter is the best one.
+  (This means looking into the last row and column of the score matrix.)
+
+The score function is currently: match: +1, mismatch: -1, indel: -2
+
+A second change in the alignment algorithm is relevant if there are
+multiple adapter occurrences in a read (such as adapter dimers).
+With the new algorithm, leftmost (earlier) adapter occurrences are now more
+reliably preferred even if a later match has fewer errors.
+
+Here are two examples from the SRR452441 dataset (R1 only),
+trimmed with the standard Illumina adapter.
+The top row shows the alignment as found by the previous algorithm,
+the middle row shows the sequencing read,
+and the last row shows the alignment as found by the updated algorithm. ::
+
+    @SRR452441.2151945
+                                             AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC  Previous alignment
+                                             ||||||||||||||||||||||||||||||||||
+    -GATCGGAAGAGCACACGTCTGAACTCCAGTCACGCACACGAGATCGGAAGAGCACACGTCTGAACTCCAGTCACGCACACGAATCTCGTATGCCGTCTTCT
+    X|||||||||||||||||||||||||||||||||
+    AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC  New alignment
+
+Previously the read was trimmed to the first 40 bases,
+now the earlier, nearly full-length occurrence is taken into account,
+and the read is empty after trimming. ::
+
+    @SRR452441.2157038
+                                    AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC  Previous alignment
+                                    ||||||||||||||||||||||||||||||||||
+    -GATCGGAAGAGCACACGTCTGAACTCCAGTCAGATCGGAAGAGCACACGTCTGAACTCCAGTCACGCACACGAATCTCGTATGCCGTCTTCTGCTTGAAAA
+    X||||||||||||||||||||||||||||||||X
+    AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC  New alignment
+
+Only very few reads should be affected by the above changes (in SRR452441, which has 2.2 million
+reads, only four reads were trimmed differently).
+In those cases where it matters, however, there should now be fewer surprises.
+
+
 .. _quality-trimming-algorithm:
 
 Quality trimming algorithm
@@ -104,4 +188,3 @@ The numbers in parentheses are not computed (because 8 is greater than zero),
 but shown here for completeness. The position of the minimum (-25) is used as
 the trimming position. Therefore, the read is trimmed to the first four bases,
 which have quality values 42, 40, 26, 27.
-
diff --git a/src/cutadapt/_align.pyx b/src/cutadapt/_align.pyx
index c47282e..0269ea4 100644
--- a/src/cutadapt/_align.pyx
+++ b/src/cutadapt/_align.pyx
@@ -1,7 +1,11 @@
 # cython: profile=False, emit_code_comments=False, language_level=3
 from cpython.mem cimport PyMem_Malloc, PyMem_Free, PyMem_Realloc
 
-DEF MATCH_SCORE = 16384
+DEF MATCH_SCORE = +1
+DEF MISMATCH_SCORE = -1
+DEF INSERTION_SCORE = -2
+DEF DELETION_SCORE = -2
+
 
 # structure for a DP matrix entry
 ctypedef struct _Entry:
@@ -13,7 +17,6 @@ ctypedef struct _Entry:
 ctypedef struct _Match:
     int origin
     int cost
-    int matches
     int score
     int ref_stop
     int query_stop
@@ -118,15 +121,13 @@ class DPMatrix:
         return '\n'.join(rows)
 
 
-cdef int score_to_matches(int score) nogil:
-    return (score + MATCH_SCORE - 1) // MATCH_SCORE
-
-
 cdef class Aligner:
     """
     Find a full or partial occurrence of a query string in a reference string
     allowing errors (mismatches, insertions, deletions).
 
+    This is a hybrid alignment algorithm that uses both costs and scores.
+
     By default, unit costs are used, meaning that mismatches, insertions and
     deletions are counted as one error (edit distance).
 
@@ -140,7 +141,7 @@ cdef class Aligner:
 
     Conceptually, the algorithm used here instead tests all possible overlaps
     between the two sequences and chooses the overlap which maximizes the
-    number of matches in the overlapping part, while the error rate must not
+    score in the overlapping part, while the error rate must not
     go above a threshold.
 
     TODO working here
@@ -176,16 +177,16 @@ cdef class Aligner:
     An optimal alignment fulfills all of these criteria:
 
     - its error_rate is at most max_error_rate
-    - Among those alignments with error_rate <= max_error_rate, the alignment contains
-      a maximal number of matches (there is no alignment with more matches).
-    - If there are multiple alignments with the same no. of matches, then one that
+    - Among those alignments with error_rate <= max_error_rate, the alignment has
+      highest score
+    - If there are multiple alignments with the same score, then one that
       has minimal no. of errors is chosen.
     - If there are still multiple candidates, choose the alignment that starts at the
       leftmost position within the read.
 
     The alignment itself is not returned, only the tuple
-    (ref_start, ref_stop, query_start, query_stop, matches, errors).
-    matches is the number of matches and errors is the number of
+    (ref_start, ref_stop, query_start, query_stop, score, errors).
+    score is the total score and errors is the number of
     errors in the alignment.
 
     It is always the case that at least one of query_start and reference_start is zero.
@@ -216,6 +217,10 @@ cdef class Aligner:
         bytes _reference  # internal, bytes version of reference (possibly translated to a non-ASCII representation)
         readonly int effective_length
         int* n_counts  # n_counts[i] == number of N characters in reference[:i]
+        int _match_score
+        int _mismatch_score
+        int _insertion_score
+        int _deletion_score
 
     def __cinit__(
         self,
@@ -243,6 +248,11 @@ cdef class Aligner:
         self._insertion_cost = indel_cost
         self._deletion_cost = indel_cost
 
+        self._match_score = MATCH_SCORE
+        self._mismatch_score = MISMATCH_SCORE
+        self._insertion_score = INSERTION_SCORE
+        self._deletion_score = DELETION_SCORE
+
     def _compute_flags(self):
         cdef int flags = 0
         if self.start_in_reference:
@@ -311,15 +321,14 @@ cdef class Aligner:
 
     def locate(self, str query):
         """
-        locate(query) -> (refstart, refstop, querystart, querystop, matches, errors)
+        locate(query) -> (refstart, refstop, querystart, querystop, score, errors)
 
         Find the query within the reference associated with this aligner. The
         intervals (querystart, querystop) and (refstart, refstop) give the
         location of the match.
 
         That is, the substrings query[querystart:querystop] and
-        self.reference[refstart:refstop] were found to align best to each other,
-        with the given number of matches and the given number of errors.
+        self.reference[refstart:refstop] were found to align best to each other.
 
         The alignment itself is not returned.
         """
@@ -409,7 +418,6 @@ cdef class Aligner:
         best.cost = m + n + 1
         best.origin = 0
         best.score = 0
-        best.matches = 0
 
         # Ukkonen's trick: index of the last cell that is at most k
         cdef int last = min(m, k + 1)
@@ -420,11 +428,12 @@ cdef class Aligner:
             int cost_diag
             int cost_deletion
             int cost_insertion
-            int origin, cost, matches, score
+            int origin, cost, score
             int length
             int ref_start
             int cur_effective_length
             int last_filled_i = 0
+            int best_length
             bint characters_equal
             bint is_acceptable
             # We keep only a single column of the DP matrix in memory.
@@ -454,14 +463,8 @@ cdef class Aligner:
                         cost = diag_entry.cost
                         origin = diag_entry.origin
                         # Among the optimal alignments whose edit distance is within the
-                        # maximum allowed error rate, we prefer those that have the
-                        # maximum number of matches. And among those, we want to
-                        # prefer those with mismatches over those with indels.
-                        # To achieve this, scores are as follows
-                        # - match: MATCH_SCORE
-                        # - mismatch: 0
-                        # - indel: -1
-                        score = diag_entry.score + MATCH_SCORE
+                        # maximum allowed error rate, we prefer those with maximal score.
+                        score = diag_entry.score + self._match_score
                     else:
                         # Characters do not match.
                         cost_diag = diag_entry.cost + 1
@@ -472,19 +475,19 @@ cdef class Aligner:
                             # MISMATCH
                             cost = cost_diag
                             origin = diag_entry.origin
-                            score = diag_entry.score
+                            score = diag_entry.score + self._mismatch_score
                         elif cost_insertion <= cost_deletion:
                             # INSERTION
                             cost = cost_insertion
                             origin = column[i-1].origin
                             # penalize insertions slightly
-                            score = column[i-1].score - 1
+                            score = column[i-1].score + self._insertion_score
                         else:
                             # DELETION
                             cost = cost_deletion
                             origin = column[i].origin
                             # penalize deletions slightly
-                            score = column[i].score - 1
+                            score = column[i].score + self._deletion_score
 
                     # Remember the current cell for next iteration
                     diag_entry = column[i]
@@ -506,7 +509,10 @@ cdef class Aligner:
                 elif stop_in_query:
                     # Found a match. If requested, find best match in last row.
                     # length of the aligned part of the reference
-                    length = m + min(column[m].origin, 0)
+                    cost = column[m].cost
+                    score = column[m].score
+                    origin = column[m].origin
+                    length = m + min(origin, 0)
                     cur_effective_length = length
                     if self.wildcard_ref:
                         if length < m:
@@ -515,44 +521,38 @@ cdef class Aligner:
                             cur_effective_length = length - self.n_counts[length]
                         else:
                             cur_effective_length = self.effective_length
-                    cost = column[m].cost
-                    score = column[m].score
-                    matches = score_to_matches(score)
-                    origin = column[m].origin
                     is_acceptable = (
                         length >= self._min_overlap
                         and cost <= cur_effective_length * max_error_rate
                     )
+                    best_length = m + min(best.origin, 0)
+
+                    # Update if
+                    # - this is the first occurrence
+                    # - or this occurrence is longer
+                    # - or if this occurrence overlaps the previous best one and has a higher score
                     if is_acceptable and (
-                        # no best match recorded so far
-                        (best.cost == m + n + 1)
-                        # same start position, use score to judge whether this is better
-                        or (origin == best.origin and score > best.score)
-                        # different start position, only matches count
-                        or (matches > best.matches or (matches == best.matches and cost < best.cost))
+                        (best.cost == m + n + 1)  # No best match recorded so far, this is the first one
+                        or (origin <= best.origin + m // 2 and score > best.score)  # This match overlaps the previous best one sufficiently (and has higher score)
+                        or (length > best_length and score > best.score)  # Length is greater than best length so far
                     ):
-                        # The case "matches == best.matches and cost < best.cost" applies
-                        # when the query contains the reference twice, and where the right
-                        # occurrence contains fewer errors
                         best.score = score
-                        best.matches = matches
                         best.cost = cost
                         best.origin = origin
                         best.ref_stop = m
                         best.query_stop = j
-                        if cost == 0 and matches == m:
+                        if cost == 0 and origin >= 0:
                             # exact match, stop early
                             break
                 # column finished
 
         if max_n == n:
             first_i = 0 if self.stop_in_reference else m
-            # search in last column # TODO last?
-            for i in range(first_i, last_filled_i + 1):
+            # search in last column
+            for i in reversed(range(first_i, last_filled_i + 1)):
                 length = i + min(column[i].origin, 0)
                 cost = column[i].cost
                 score = column[i].score
-                matches = score_to_matches(score)
                 if self.wildcard_ref:
                     if length < m:
                         # Recompute effective length so that it only takes into
@@ -571,15 +571,13 @@ cdef class Aligner:
                     length >= self._min_overlap
                     and cost <= cur_effective_length * max_error_rate
                 )
+                best_length = best.ref_stop + min(best.origin, 0)
+
                 if is_acceptable and (
-                    # no best match recorded so far
                     (best.cost == m + n + 1)
-                    # same start position, use score to judge whether this is better
-                    or (origin == best.origin and score > best.score)
-                    # different start position, only matches count
-                    or (matches > best.matches or (matches == best.matches and cost < best.cost))
+                    or (origin <= best.origin + m // 2 and score > best.score)
+                    or (length > best_length and score > best.score)
                 ):
-                    best.matches = matches
                     best.score = score
                     best.cost = cost
                     best.origin = column[i].origin
@@ -599,7 +597,7 @@ cdef class Aligner:
             ref_start = -best.origin
             query_start = 0
 
-        return (ref_start, best.ref_stop, query_start, best.query_stop, best.matches, best.cost)
+        return (ref_start, best.ref_stop, query_start, best.query_stop, best.score, best.cost)
 
     def __dealloc__(self):
         PyMem_Free(self.column)
@@ -670,17 +668,18 @@ cdef class PrefixComparer:
         This is used to find an anchored 5' adapter (type 'FRONT') in the 'no indels' mode.
         This is very simple as only the number of errors needs to be counted.
 
-        This function returns a tuple compatible with what Aligner.locate outputs.
+        This function returns a tuple compatible with what Aligner.locate returns.
         """
         cdef:
             bytes query_bytes = query.encode('ascii')
             char* r_ptr = self.reference
             char* q_ptr
-            int i, matches = 0
+            int i
             int n = len(query_bytes)
             int length = min(self.m, n)
             bint compare_ascii = False
-            int errors
+            int errors = 0
+            int score
 
         if self.wildcard_query:
             query_bytes = query_bytes.translate(IUPAC_TABLE)
@@ -693,17 +692,17 @@ cdef class PrefixComparer:
 
         if compare_ascii:
             for i in range(length):
-                if r_ptr[i] == q_ptr[i]:
-                    matches += 1
+                if r_ptr[i] != q_ptr[i]:
+                    errors += 1
         else:
             for i in range(length):
-                if (r_ptr[i] & q_ptr[i]) != 0:
-                    matches += 1
+                if (r_ptr[i] & q_ptr[i]) == 0:
+                    errors += 1
 
-        errors = length - matches
         if errors > self.max_k or length < self.min_overlap:
             return None
-        return (0, length, 0, length, matches, length - matches)
+        score = (length - errors) * MATCH_SCORE + errors * MISMATCH_SCORE
+        return (0, length, 0, length, score, errors)
 
 
 cdef class SuffixComparer(PrefixComparer):
@@ -723,5 +722,5 @@ cdef class SuffixComparer(PrefixComparer):
         result = super().locate(query[::-1])
         if result is None:
             return None
-        _, length, _, _, matches, errors = result
-        return (self.m - length, self.m, n - length, n, matches, errors)
+        _, length, _, _, score, errors = result
+        return (self.m - length, self.m, n - length, n, score, errors)
diff --git a/src/cutadapt/adapters.py b/src/cutadapt/adapters.py
index 5012e69..9faadb0 100644
--- a/src/cutadapt/adapters.py
+++ b/src/cutadapt/adapters.py
@@ -308,7 +308,7 @@ class SingleMatch(Match, ABC):
         "astop",
         "rstart",
         "rstop",
-        "matches",
+        "score",
         "errors",
         "adapter",
         "sequence",
@@ -322,7 +322,7 @@ class SingleMatch(Match, ABC):
         astop: int,
         rstart: int,
         rstop: int,
-        matches: int,
+        score: int,
         errors: int,
         adapter: "SingleAdapter",
         sequence: str,
@@ -331,7 +331,7 @@ class SingleMatch(Match, ABC):
         self.astop: int = astop
         self.rstart: int = rstart
         self.rstop: int = rstop
-        self.matches: int = matches
+        self.score: int = score
         self.errors: int = errors
         self.adapter: SingleAdapter = adapter
         self.sequence = sequence
@@ -341,8 +341,10 @@ class SingleMatch(Match, ABC):
         self.length: int = astop - astart
 
     def __repr__(self):
-        return "SingleMatch(astart={}, astop={}, rstart={}, rstop={}, matches={}, errors={})".format(
-            self.astart, self.astop, self.rstart, self.rstop, self.matches, self.errors
+        return (
+            f"{self.__class__.__name__}(astart={self.astart}, astop={self.astop}, "
+            f"rstart={self.rstart}, rstop={self.rstop}, "
+            f"score={self.score}, errors={self.errors})"
         )
 
     def __eq__(self, other) -> bool:
@@ -352,7 +354,7 @@ class SingleMatch(Match, ABC):
             and self.astop == other.astop
             and self.rstart == other.rstart
             and self.rstop == other.rstop
-            and self.matches == other.matches
+            and self.score == other.score
             and self.errors == other.errors
             and self.adapter is other.adapter
             and self.sequence == other.sequence
@@ -410,11 +412,6 @@ class SingleMatch(Match, ABC):
 class RemoveBeforeMatch(SingleMatch):
     """A match that removes sequence before the match"""
 
-    def __repr__(self):
-        return "RemoveBeforeMatch(astart={}, astop={}, rstart={}, rstop={}, matches={}, errors={})".format(
-            self.astart, self.astop, self.rstart, self.rstop, self.matches, self.errors
-        )
-
     def rest(self) -> str:
         """
         Return the part of the read before this match if this is a
@@ -448,11 +445,6 @@ class RemoveBeforeMatch(SingleMatch):
 class RemoveAfterMatch(SingleMatch):
     """A match that removes sequence after the match"""
 
-    def __repr__(self):
-        return "RemoveAfterMatch(astart={}, astop={}, rstart={}, rstop={}, matches={}, errors={})".format(
-            self.astart, self.astop, self.rstart, self.rstop, self.matches, self.errors
-        )
-
     def rest(self) -> str:
         """
         Return the part of the read before this match if this is a
@@ -910,14 +902,14 @@ class LinkedMatch(Match):
         )
 
     @property
-    def matches(self):
+    def score(self):
         """Number of matching bases"""
-        m = 0
+        s = 0
         if self.front_match is not None:
-            m += self.front_match.matches
+            s += self.front_match.score
         if self.back_match is not None:
-            m += self.back_match.matches
-        return m
+            s += self.back_match.score
+        return s
 
     @property
     def errors(self):
@@ -1070,13 +1062,12 @@ class MultipleAdapters(Matchable):
             if match is None:
                 continue
 
-            # the no. of matches determines which adapter fits best
+            # the score determines which adapter fits best
             if (
                 best_match is None
-                or match.matches > best_match.matches
+                or match.score > best_match.score
                 or (
-                    match.matches == best_match.matches
-                    and match.errors < best_match.errors
+                    match.score == best_match.score and match.errors < best_match.errors
                 )
             ):
                 best_match = match
@@ -1258,13 +1249,13 @@ class IndexedPrefixAdapters(IndexedAdapters):
             raise ValueError("Only 5' anchored adapters are allowed")
         return super()._accept(adapter)
 
-    def _make_match(self, adapter, length, matches, errors, sequence):
+    def _make_match(self, adapter, length, score, errors, sequence):
         return RemoveBeforeMatch(
             astart=0,
             astop=len(adapter.sequence),
             rstart=0,
             rstop=length,
-            matches=matches,
+            score=score,
             errors=errors,
             adapter=adapter,
             sequence=sequence,
@@ -1285,13 +1276,13 @@ class IndexedSuffixAdapters(IndexedAdapters):
             raise ValueError("Only anchored 3' adapters are allowed")
         return super()._accept(adapter)
 
-    def _make_match(self, adapter, length, matches, errors, sequence):
+    def _make_match(self, adapter, length, score, errors, sequence):
         return RemoveAfterMatch(
             astart=0,
             astop=len(adapter.sequence),
             rstart=len(sequence) - length,
             rstop=len(sequence),
-            matches=matches,
+            score=score,
             errors=errors,
             adapter=adapter,
             sequence=sequence,
diff --git a/src/cutadapt/align.py b/src/cutadapt/align.py
index 44110f1..c8129c8 100644
--- a/src/cutadapt/align.py
+++ b/src/cutadapt/align.py
@@ -102,7 +102,8 @@ def hamming_environment(s: str, k: int) -> Iterator[Tuple[str, int, int]]:
 
 def naive_edit_environment(s: str, k: int) -> Iterator[str]:
     """
-    Apply all possible edits up to edit distance k to string s.
+    Apply all possible edits up to edit distance k to string s
+    and yield the resulting strings.
     A string may be returned more than once.
     """
     yield s
@@ -126,8 +127,8 @@ def edit_environment(s: str, k: int) -> Iterator[Tuple[str, int, int]]:
     Find all strings t for which the edit distance between s and t is at most k,
     assuming the alphabet is A, C, G, T.
 
-    Yield tuples (t, e, m), where e is the edit distance between s and t and
-    m is the number of matches in the optimal alignment.
+    Yield tuples (t, e, score), where e is the edit distance between s and t and
+    score is the score of the optimal alignment.
     """
     rate = k / len(s) if s else 0
     aligner = Aligner(s, max_error_rate=rate, flags=0, min_overlap=len(s))
@@ -137,8 +138,8 @@ def edit_environment(s: str, k: int) -> Iterator[Tuple[str, int, int]]:
             continue
         seen.add(t)
         result = aligner.locate(t)
-        matches, errors = result[-2:]  # type: ignore
-        yield t, errors, matches
+        score, errors = result[-2:]  # type: ignore
+        yield t, errors, score
 
 
 def slow_edit_environment(s: str, k: int) -> Iterator[Tuple[str, int, int]]:
@@ -148,6 +149,8 @@ def slow_edit_environment(s: str, k: int) -> Iterator[Tuple[str, int, int]]:
 
     Yield tuples (t, e, m), where e is the edit distance between s and t and
     m is the number of matches in the optimal alignment.
+
+    This is slow and only used in testing.
     """
     n = len(s)
     alphabet = "TGCA"
diff --git a/src/cutadapt/modifiers.py b/src/cutadapt/modifiers.py
index ac418bf..d177a36 100644
--- a/src/cutadapt/modifiers.py
+++ b/src/cutadapt/modifiers.py
@@ -302,9 +302,9 @@ class ReverseComplementer(SingleEndModifier):
             reverse_read
         )
 
-        forward_match_count = sum(m.matches for m in forward_matches)
-        reverse_match_count = sum(m.matches for m in reverse_matches)
-        use_reverse_complement = reverse_match_count > forward_match_count
+        forward_score = sum(m.score for m in forward_matches)
+        reverse_score = sum(m.score for m in reverse_matches)
+        use_reverse_complement = reverse_score > forward_score
 
         if use_reverse_complement:
             self.reverse_complemented += 1
@@ -404,7 +404,7 @@ class PairedAdapterCutter(PairedEndModifier):
         self, sequence1: str, sequence2: str
     ) -> Optional[Tuple[Match, Match]]:
         best = None
-        best_matches = None
+        best_score = None
         best_errors = None
         for adapter1, adapter2 in self._adapter_pairs:
             match1 = adapter1.match_to(sequence1)
@@ -413,15 +413,15 @@ class PairedAdapterCutter(PairedEndModifier):
             match2 = adapter2.match_to(sequence2)
             if match2 is None:
                 continue
-            total_matches = match1.matches + match2.matches
+            total_score = match1.score + match2.score
             total_errors = match1.errors + match2.errors
             if (
                 best is None
-                or total_matches > best_matches
-                or (total_matches == best_matches and total_errors < best_errors)
+                or total_score > best_score
+                or (total_score == best_score and total_errors < best_errors)
             ):
                 best = match1, match2
-                best_matches = total_matches
+                best_score = total_score
                 best_errors = total_errors
         return best
 


Test patch:
diff --git a/tests/test_adapters.py b/tests/test_adapters.py
index 5af9f86..6696f69 100644
--- a/tests/test_adapters.py
+++ b/tests/test_adapters.py
@@ -56,7 +56,7 @@ def test_issue_52():
         astop=17,
         rstart=5,
         rstop=21,
-        matches=15,
+        score=15,
         errors=2,
         adapter=adapter,
         sequence=sequence,
@@ -76,16 +76,21 @@ def test_issue_52():
 
 
 def test_issue_80():
-    # This issue turned out to not be an actual issue with the alignment
-    # algorithm. The following alignment is found because it has more matches
-    # than the 'obvious' one:
+    # This issue was at the time not considered to be an actual issue with the alignment
+    # algorithm. The following alignment with three errors was found because it had more
+    # matches than the 'obvious' one:
     #
     # TCGTATGCCGTCTTC
     # =========X==XX=
     # TCGTATGCCCTC--C
     #
-    # This is correct, albeit a little surprising, since an alignment without
-    # indels would have only two errors.
+    # The alignment algorithm has since been changed so that not the number of matches
+    # is relevant, but a score that penalizes indels. Now, the resulting alignment
+    # should be this one (with only two errors):
+    #
+    # TCGTATGCCGTCTTC
+    # =========X==X
+    # TCGTATGCCCTCC
 
     adapter = BackAdapter(
         sequence="TCGTATGCCGTCTTC",
@@ -95,12 +100,11 @@ def test_issue_80():
         adapter_wildcards=False,
     )
     result = adapter.match_to("TCGTATGCCCTCC")
-    assert result.errors == 3, result
+    assert result.errors == 2, result
     assert result.astart == 0, result
-    assert result.astop == 15, result
+    assert result.astop == 13, result
 
 
-@pytest.mark.xfail(strict=True)
 def test_back_adapter_indel_and_exact_occurrence():
     adapter = BackAdapter(
         sequence="GATCGGAAGA",
@@ -112,12 +116,12 @@ def test_back_adapter_indel_and_exact_occurrence():
     # GATCGTGAAGAGATCGGAAGA
     # GATCG-GAAGA
     #            GATCGGAAGA
-    assert match.errors == 0
-    assert match.matches == 10
     assert match.astart == 0
     assert match.astop == 10
     assert match.rstart == 0
-    assert match.rstop == 10
+    assert match.rstop == 11
+    assert match.errors == 1
+    assert match.score == 8
 
 
 def test_back_adapter_indel_and_mismatch_occurrence():
@@ -130,12 +134,12 @@ def test_back_adapter_indel_and_mismatch_occurrence():
     # CTGGATCGGA-GAGCCGTAGATCGGGAGAGGC
     #    ||||||| ||      ||||||X|||
     #    GATCGGAAGA      GATCGGAAGA
-    assert match.errors == 1
-    assert match.matches == 9
     assert match.astart == 0
     assert match.astop == 10
     assert match.rstart == 3
     assert match.rstop == 12
+    assert match.score == 7
+    assert match.errors == 1
 
 
 def test_str():
@@ -153,13 +157,16 @@ def test_prefix_with_indels_one_mismatch():
         adapter_wildcards=False,
         indels=True,
     )
+    # GCACATCGGAA
+    # |||||||X
+    # GCACATCT
     result = a.match_to("GCACATCGGAA")
-    assert result.errors == 1
-    assert result.matches == 7
     assert result.astart == 0
     assert result.astop == 8
     assert result.rstart == 0
     assert result.rstop == 8
+    assert result.score == 6  # 7 matches, 1 mismatch
+    assert result.errors == 1
 
 
 def test_prefix_with_indels_two_mismatches():
@@ -172,12 +179,15 @@ def test_prefix_with_indels_two_mismatches():
         indels=True,
     )
     result = a.match_to("GCACATCGGAA")
-    assert result.errors == 2
-    assert result.matches == 6
+    # GCACATCGGAA
+    # ||||||XX
+    # GCACATTT
     assert result.astart == 0
     assert result.astop == 8
     assert result.rstart == 0
     assert result.rstop == 8
+    assert result.score == 4
+    assert result.errors == 2
 
 
 def test_linked_adapter():
@@ -215,7 +225,7 @@ def test_info_record():
         astop=17,
         rstart=5,
         rstop=21,
-        matches=15,
+        score=15,
         errors=2,
         adapter=adapter,
         sequence=read.sequence,
@@ -329,7 +339,7 @@ def test_linked_matches_property():
         back_required=False,
         name="name",
     )
-    assert la.match_to("AAAATTTT").matches == 3
+    assert la.match_to("AAAATTTT").score == 3
 
 
 @pytest.mark.parametrize("adapter_class", [PrefixAdapter, SuffixAdapter])
@@ -445,9 +455,10 @@ def test_indexed_prefix_adapters_with_n_wildcard():
     sequence = "GGTCCAGA"
     ma = IndexedPrefixAdapters([PrefixAdapter(sequence, max_errors=1, indels=False)])
     for i in range(len(sequence)):
+        # N in the read should be counted as mismatch
         t = sequence[:i] + "N" + sequence[i + 1 :] + "TGCT"
         result = ma.match_to(t)
         assert isinstance(result, RemoveBeforeMatch)
         assert (result.rstart, result.rstop) == (0, 8)
         assert result.errors == 1
-        assert result.matches == 7
+        assert result.score == 6
diff --git a/tests/test_align.py b/tests/test_align.py
index 4a80bec..a62e705 100644
--- a/tests/test_align.py
+++ b/tests/test_align.py
@@ -1,3 +1,5 @@
+from typing import NamedTuple
+
 import pytest
 
 from cutadapt.align import (
@@ -16,6 +18,15 @@ from cutadapt.adapters import Where
 from utils import binomial
 
 
+class AlignmentResult(NamedTuple):
+    ref_start: int
+    ref_end: int
+    query_start: int
+    query_end: int
+    score: int
+    errors: int
+
+
 # convenience function (to avoid having to instantiate an Aligner manually)
 def locate(
     reference,
@@ -59,12 +70,78 @@ class TestAligner:
         result = aligner.locate("")
         assert (0, 0, 0, 0, 0, 0) == result
 
-
-def test_polya():
+    def test_indels_penalized(self):
+        # Alignment in older versions:
+        # CCAGTCCTTTCCTGAGAGT
+        # CCAGTCCT---CT
+        #
+        # Should now be:
+        # CCAGTCCTTTCCTGAGAGT
+        # CCAGTCCTCT
+        aligner = Aligner("CCAGTCCTCT", 0.3, flags=Where.PREFIX)
+        result = aligner.locate("CCAGTCCTTTCCTGAGAGT")
+        assert (0, 10, 0, 10, 9 - 1, 1) == result
+        # refstart, refstop, querystart, querystop, score, errors)
+
+        # Alignment:
+        # TCGATGC
+        # TCGATC
+        aligner = Aligner("TCGATC", 1.5 / 6, flags=Where.PREFIX)
+        result = aligner.locate("TCGATGC")
+        assert (0, 6, 0, 6, 4, 1) == result
+
+    def test_align_illumina(self):
+        aligner = Aligner("GCCGAACTTCTTAGACTGCCTTAAGGACGT", 0.1, flags=Where.BACK)
+        result = AlignmentResult(
+            *aligner.locate("CAAATCACCAGAAGGCGCCTAACTTCTTAGACTGCC")
+        )
+        #                 GCCGAACTTCTTAGACTGCCTTAAGGACGT (ref)
+        #                 |||X||||||||||||||||
+        # CAAATCACCAGAAGGCGCCTAACTTCTTAGACTGCC         (query)
+        assert result.ref_start == 0
+        assert result.ref_end == 20
+        assert result.query_start == 16
+        assert result.query_end == 36
+        assert result.score == 18
+        assert result.errors == 1
+
+
+def test_poly_t():
+    aligner = Aligner("TTTT", 0.25, flags=Where.BACK)
+    result = AlignmentResult(*aligner.locate("CCTTTT"))
+    assert result.ref_start == 0
+    assert result.ref_end == 4
+    assert result.query_start == 2
+    assert result.query_end == 6
+    assert result.score == 4
+    assert result.errors == 0
+
+
+def test_poly_t_partial_match():
+    aligner = Aligner("TTTTTT", 0.25, flags=Where.BACK)
+    result = AlignmentResult(*aligner.locate("CCTTTT"))
+    assert result.ref_start == 0
+    assert result.ref_end == 4
+    assert result.query_start == 2
+    assert result.query_end == 6
+    assert result.score == 4
+    assert result.errors == 0
+
+
+def test_poly_t_2():
+    aligner = Aligner("TTT", 1 / 3, flags=Where.BACK)
+    result = AlignmentResult(*aligner.locate("CCTTTT"))
+    assert result.ref_start == 0
+    assert result.ref_end == 3
+    assert result.query_start == 2
+    assert result.query_end == 5
+
+
+def test_poly_a():
     s = "AAAAAAAAAAAAAAAAA"
     t = "ACAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"
     result = locate(s, t, 0.0, Where.BACK.value)
-    # start_s, stop_s, start_t, stop_t, matches, cost = result
+    # start_s, stop_s, start_t, stop_t, score, cost = result
     assert result == (0, len(s), 4, 4 + len(s), len(s), 0)
 
 
@@ -111,7 +188,7 @@ def compare_suffixes(ref, query, wildcard_ref=False, wildcard_query=False):
 
 
 def test_compare_prefixes():
-    assert compare_prefixes("AAXAA", "AAAAATTTTTTTTT") == (0, 5, 0, 5, 4, 1)
+    assert compare_prefixes("AAXAA", "AAAAATTTTTTTTT") == (0, 5, 0, 5, 3, 1)
     assert compare_prefixes("AANAA", "AACAATTTTTTTTT", wildcard_ref=True) == (
         0,
         5,
@@ -128,7 +205,7 @@ def test_compare_prefixes():
         5,
         0,
     )
-    assert compare_prefixes("XAAAAA", "AAAAATTTTTTTTT") == (0, 6, 0, 6, 4, 2)
+    assert compare_prefixes("XAAAAA", "AAAAATTTTTTTTT") == (0, 6, 0, 6, 2, 2)
 
     a = WILDCARD_SEQUENCES[0]
     for s in WILDCARD_SEQUENCES:
@@ -157,7 +234,7 @@ def test_compare_prefixes():
             result = compare_prefixes(
                 "CCCXTTXATC", r, wildcard_ref=wildc_ref, wildcard_query=wildc_query
             )
-            assert result == (0, 10, 0, 10, 8, 2)
+            assert result == (0, 10, 0, 10, 6, 2)
 
 
 def test_n_wildcard_in_ref_matches_n_wildcard_in_query_prefix():
@@ -184,7 +261,7 @@ def test_n_wildcard_in_ref_matches_n_wildcard_in_query_back():
 
 
 def test_compare_suffixes():
-    assert compare_suffixes("AAXAA", "TTTTTTTAAAAA") == (0, 5, 7, 12, 4, 1)
+    assert compare_suffixes("AAXAA", "TTTTTTTAAAAA") == (0, 5, 7, 12, 3, 1)
     assert compare_suffixes("AANAA", "TTTTTTTAACAA", wildcard_ref=True) == (
         0,
         5,
@@ -201,7 +278,7 @@ def test_compare_suffixes():
         5,
         0,
     )
-    assert compare_suffixes("AAAAAX", "TTTTTTTAAAAA") == (0, 6, 6, 12, 4, 2)
+    assert compare_suffixes("AAAAAX", "TTTTTTTAAAAA") == (0, 6, 6, 12, 2, 2)
 
 
 @pytest.mark.parametrize("upper", (True, False))
@@ -403,13 +480,14 @@ def test_edit_environment(k, s, environment_func):
     aligner = Aligner(s, max_error_rate=error_rate, flags=0, min_overlap=len(s))
     for t, dist, m in result:
         result = aligner.locate(t)
-        start1, stop1, start2, stop2, matches, errors = result
+        start1, stop1, start2, stop2, score, errors = result
         assert errors == dist
-        assert m == matches
         assert start1 == 0
         assert stop1 == len(s)
         assert start2 == 0
         assert stop2 == len(t)
         assert edit_distance(s, t) == dist
-        assert m <= len(s), (s, t, dist)
-        assert m <= len(t), (s, t, dist)
+        if environment_func is edit_environment:
+            assert m == score
+            assert m <= len(s), (s, t, dist)
+            assert m <= len(t), (s, t, dist)
diff --git a/tests/test_modifiers.py b/tests/test_modifiers.py
index b801ebb..29c4bdd 100644
--- a/tests/test_modifiers.py
+++ b/tests/test_modifiers.py
@@ -275,7 +275,7 @@ class TestRenamer:
                 astop=4,
                 rstart=8,
                 rstop=12,
-                matches=3,
+                score=3,
                 errors=1,
                 adapter=adapter,
                 sequence=sequence,
@@ -305,7 +305,7 @@ class TestRenamer:
             astop=4,
             rstart=0,
             rstop=4,
-            matches=3,
+            score=3,
             errors=1,
             adapter=adapter1,
             sequence=sequence,
@@ -315,7 +315,7 @@ class TestRenamer:
             astop=4,
             rstart=4,
             rstop=8,
-            matches=3,
+            score=3,
             errors=1,
             adapter=adapter2,
             sequence=sequence[4:],
@@ -392,7 +392,7 @@ class TestPairedEndRenamer:
                 astop=4,
                 rstart=1,
                 rstop=3,
-                matches=1,
+                score=1,
                 errors=1,
                 adapter=FrontAdapter("AT"),
                 sequence="AACC",
@@ -406,7 +406,7 @@ class TestPairedEndRenamer:
                 astop=4,
                 rstart=1,
                 rstop=3,
-                matches=1,
+                score=1,
                 errors=1,
                 adapter=FrontAdapter("GA"),
                 sequence="GGTT",