fix(features): correct positional modification indexing, and optional terminal composition

deeplc/_features.py

@@ -31,6 +31,7 @@
 def encode_peptidoform(
     peptidoform: Peptidoform | str,
     add_ccs_features: bool = False,
+    add_terminal_composition: bool = False,
     padding_length: int = 60,
     positions: set[int] | None = None,
     positions_pos: set[int] | None = None,
@@ -48,6 +49,12 @@ def encode_peptidoform(
         The peptidoform to encode, either as a Peptidoform object or a string.
     add_ccs_features
         Whether to include CCS features. Default is False.
+    add_terminal_composition
+        Whether to append the N- and C-terminal group compositions to
+        ``matrix_global``, adding 12 values. Default is False, which keeps
+        ``matrix_global`` at its existing width so models trained against it
+        remain valid. Without this, a terminal modification and a side-chain
+        modification on the same residue are indistinguishable.
     padding_length
         The maximum length of the sequence after padding. Default is 60.
     positions
@@ -125,6 +132,10 @@ def encode_peptidoform(
     matrix_sum = _compute_rolling_sum(std_matrix.T, n=2)[:, ::2].T
 
     matrix_global = np.concatenate([matrix_all, pos_matrix.flatten()])
+    if add_terminal_composition:
+        matrix_global = np.concatenate(
+            [matrix_global, _terminal_composition(peptidoform, dict_index).flatten()]
+        )
 
     return {
         "matrix": std_matrix,
@@ -206,6 +217,55 @@ def _fill_pos_matrix(
     return pos_mat
 
 
+def _positional_rows(i: int, seq_len: int, positions: set[int]) -> list[int]:
+    """
+    Rows of the positional matrix that sequence position ``i`` occupies.
+
+    Two corrections over indexing ``pos_mat`` by ``i`` directly:
+
+    * ``_fill_pos_matrix`` lays rows out as ``sorted(positions)``, so row 0 is
+      ``min(positions)``. Raw indexing puts an N-terminal delta in the row that
+      means "fourth residue from the C-terminus".
+    * In a short peptide one residue can occupy both a positive and a negative
+      row, for example index 3 of a 7-mer is also -4. ``_fill_pos_matrix``
+      writes the base residue to both, so a delta must reach both as well.
+    """
+    offset = min(positions)
+    rows = []
+    if i in positions:
+        rows.append(i - offset)
+    if (i - seq_len) in positions:
+        rows.append((i - seq_len) - offset)
+    return rows
+
+
+def _terminal_composition(
+    peptidoform: Peptidoform,
+    dict_index: dict[str, int],
+) -> np.ndarray:
+    """
+    Composition of the N- and C-terminal groups, as two stacked atom vectors.
+
+    Terminal groups are already folded into ``matrix`` and the positional block,
+    but there they are indistinguishable from a modification on the side chain
+    of the first or last residue. ``[Acetyl]-PEPTIDEK`` and ``P[Acetyl]EPTIDEK``
+    otherwise produce identical features, and they do not elute alike.
+    """
+    out = np.zeros((2, len(dict_index)), dtype=np.float16)
+    for row, key in ((0, "n_term"), (1, "c_term")):
+        for tag in peptidoform.properties.get(key) or []:
+            try:
+                composition = tag.composition
+            except Exception:
+                warnings.warn(f"No composition for terminal modification {tag}", stacklevel=2)
+                continue
+            for atom, change in composition.items():
+                index = dict_index.get(atom, dict_index.get(sub(r"\[.*?\]", "", atom)))
+                if index is not None:
+                    out[row, index] += change
+    return out
+
+
 def _apply_composition_to_matrices(
     mat: np.ndarray,
     pos_mat: np.ndarray,
@@ -216,23 +276,26 @@ def _apply_composition_to_matrices(
     dict_index_pos: dict[str, int],
     positions: set[int],
 ) -> None:
-    """Apply a composition delta to the standard and positional matrices."""
+    """
+    Apply a composition delta to the standard and positional matrices.
+
+    Positional rows come from :func:`_positional_rows`, which applies the same
+    offset and the same both-ends handling that :func:`_fill_pos_matrix` uses
+    for base residue compositions.
+    """
+    rows = _positional_rows(i, seq_len, positions)
     for atom_comp, change in composition.items():
         try:
             mat[i, dict_index[atom_comp]] += change
-            if i in positions:
-                pos_mat[i, dict_index_pos[atom_comp]] += change
-            elif (i - seq_len) in positions:
-                pos_mat[i - seq_len, dict_index_pos[atom_comp]] += change
+            for row in rows:
+                pos_mat[row, dict_index_pos[atom_comp]] += change
         except KeyError:
             try:
                 warnings.warn(f"Replacing pattern for atom: {atom_comp}", stacklevel=2)
                 atom_comp_clean = sub(r"\[.*?\]", "", atom_comp)
                 mat[i, dict_index[atom_comp_clean]] += change
-                if i in positions:
-                    pos_mat[i, dict_index_pos[atom_comp_clean]] += change
-                elif (i - seq_len) in positions:
-                    pos_mat[i - seq_len, dict_index_pos[atom_comp_clean]] += change
+                for row in rows:
+                    pos_mat[row, dict_index_pos[atom_comp_clean]] += change
             except KeyError:
                 warnings.warn(f"Ignoring atom {atom_comp} at pos {i}", stacklevel=2)
                 continue

tests/test_features.py

@@ -262,3 +262,43 @@ def test_single_residue(self):
     def test_two_residues_no_crash(self):
         result = encode_peptidoform("AC")
         assert result["matrix_global"].shape == (_GLOBAL_BASE_LEN,)
+
+
+def test_positional_delta_uses_same_row_as_base_residue():
+    """A modification delta must land on the row of the residue it modifies.
+
+    Regression test: ``_fill_pos_matrix`` offsets rows by ``min(positions)`` so
+    that row 0 means position -4, while modification deltas were written with
+    raw indexing, putting an N-terminal delta four residues from the other end.
+    """
+    from deeplc._features import DEFAULT_POSITIONS, encode_peptidoform
+
+    order = sorted(DEFAULT_POSITIONS)
+    plain = encode_peptidoform("PEPTIDEK")["matrix_global"][7:55].reshape(8, 6)
+    n_term = encode_peptidoform("[Acetyl]-PEPTIDEK")["matrix_global"][7:55].reshape(8, 6)
+    c_term = encode_peptidoform("PEPTIDEK-[Amidated]")["matrix_global"][7:55].reshape(8, 6)
+
+    changed_n = [order[r] for r in range(8) if (n_term[r] - plain[r]).any()]
+    changed_c = [order[r] for r in range(8) if (c_term[r] - plain[r]).any()]
+
+    assert changed_n == [0], f"N-terminal delta landed at {changed_n}, expected position 0"
+    assert changed_c == [-1], f"C-terminal delta landed at {changed_c}, expected position -1"
+
+
+def test_terminal_composition_is_opt_in_and_separates_terminal_from_side_chain():
+    """``[Acetyl]-PEPTIDEK`` and ``P[Acetyl]EPTIDEK`` are chemically different."""
+    import numpy as np
+
+    from deeplc._features import encode_peptidoform
+
+    default_width = encode_peptidoform("PEPTIDEK")["matrix_global"].shape[0]
+    assert default_width == 55, "default global width must not change"
+
+    terminal = encode_peptidoform("[Acetyl]-PEPTIDEK", add_terminal_composition=True)
+    side_chain = encode_peptidoform("P[Acetyl]EPTIDEK", add_terminal_composition=True)
+
+    assert terminal["matrix_global"].shape[0] == 67
+    assert not np.allclose(terminal["matrix_global"], side_chain["matrix_global"])
+    # the acetyl composition C2H2O appears in the N-terminal block only when terminal
+    assert terminal["matrix_global"][55:61].tolist() == [2, 2, 0, 1, 0, 0]
+    assert side_chain["matrix_global"][55:61].tolist() == [0, 0, 0, 0, 0, 0]