Test aligner

.github/workflows/unittest.yaml

@@ -28,6 +28,6 @@ jobs:
       - name: Install uv
         uses: astral-sh/setup-uv@v7
       - name: Install the project
-        run: uv sync --extra unittest --dev
+        run: uv sync --dev
       - name: Run tests
         run: uv run -m test

pyproject.toml

@@ -32,5 +32,4 @@ repository = "https://github.com/BioinfoUNIBA/REDItools3"
 issues = "https://github.com/BioinfoUNIBA/REDItools3/issues"
 
 [project.optional-dependencies]
-unittest = ["biopython"]
 wps = ["flake8", "wemake-python-styleguide"]

test/aligner.py

@@ -0,0 +1,102 @@
+class Aligner:
+    def __init__(self, match=1, mismatch=1, gap=1):
+        self.match = 1
+        self.mismatch = 1
+        self.gap = 1
+
+    def align(self, ref_seq, query_seq):
+        matrix = NWMatrix(
+            ref_seq,
+            query_seq,
+            self.match,
+            self.mismatch,
+            self.gap,
+        )
+        matrix.run_dp()
+        return self.trace_matrix(
+            ref_seq,
+            query_seq,
+            matrix.trace_matrix,
+        )
+
+    def trace_matrix(self, ref_seq, qry_seq, trace_mat):
+        row_idx = len(ref_seq)
+        col_idx = len(qry_seq)
+
+        ref_align = []
+        qry_align = []
+        while row_idx > 0 or col_idx > 0:
+            trace_val = trace_mat[row_idx][col_idx]
+
+            if trace_val in [2, 5, 6, 9]:
+                row_idx -= 1
+                col_idx -= 1
+                ref_align.insert(0, ref_seq[row_idx])
+                qry_align.insert(0, qry_seq[col_idx])
+            elif trace_val in [3, 7]:
+                row_idx -= 1
+                ref_align.insert(0, ref_seq[row_idx])
+                qry_align.insert(0, '-')
+            elif trace_val == 4:
+                col_idx -= 1
+                ref_align.insert(0, '-')
+                qry_align.insert(0, qry_seq[col_idx])
+
+        return (
+            ''.join(ref_align),
+            ''.join(qry_align),
+        )
+
+class NWMatrix:
+    def __init__(self, ref_seq, query_seq, match, mismatch, gap):
+        self.ref_seq = ref_seq
+        self.qry_seq = query_seq
+        self.gap = gap
+        self.match = match
+        self.mismatch = mismatch
+
+        self.nw_matrix = self.init_nw_matrix(ref_seq, query_seq)
+        self.trace_matrix = self.init_trace_matrix(ref_seq, query_seq)
+
+    def assess_cell(self, col_idx, ref_base, row_idx, query_base):
+        align_val = self.match if ref_base == query_base else -self.mismatch
+        t_list = [
+            self.nw_matrix[row_idx][col_idx] + align_val,
+            self.nw_matrix[row_idx][col_idx + 1] - self.gap,
+            self.nw_matrix[row_idx + 1][col_idx] - self.gap,
+        ]
+        t_max = max(t_list)                
+        self.nw_matrix[row_idx + 1][col_idx + 1] = t_max
+        self.trace_matrix[row_idx + 1][col_idx + 1] += sum((
+            idx + 2 for idx, tv in enumerate(t_list) if tv == t_max
+        ))
+
+    def run_dp(self):
+        for col_idx, ref_base in enumerate(self.qry_seq):
+            for row_idx, query_base in enumerate(self.ref_seq):
+                self.assess_cell(col_idx, ref_base, row_idx, query_base)
+
+    def init_nw_matrix(self, ref_seq, query_seq):
+        nw_mat = [
+            [-self.gap * (row_idx + 1)] + \
+                [0 for _ in range(len(query_seq))]
+            for row_idx in range(len(ref_seq))
+        ]
+        nw_mat.insert(
+            0, [
+                -self.gap * _
+                for _ in range(len(query_seq) + 1)
+            ],
+        )
+        return nw_mat
+
+    def init_trace_matrix(self, ref_seq, query_seq):
+        trace_mat = [
+            [3] + [0 for _ in range(len(query_seq))]
+            for _ in range(len(ref_seq))
+        ]
+        trace_mat.insert(
+            0,
+            [4 for _ in range(len(query_seq) + 1)],
+        )
+        return trace_mat

test/sam_gen.py

@@ -3,8 +3,8 @@
 import re
 from dataclasses import InitVar, dataclass
 from tempfile import NamedTemporaryFile
+from test.aligner import Aligner
 
-from Bio.Align import PairwiseAligner
 from pysam import samtools
 
 
@@ -53,10 +53,6 @@ class Sequence:
     read_n = 0
     flag_reverse_strand = 16
     phred_default = 30
-    aligner = PairwiseAligner(
-        mismatch_score=-1,
-        query_internal_open_gap_score=-1,
-    )
 
     def __post_init__(self, phred, qname):
         if phred is None:
@@ -89,11 +85,11 @@ def tlen(self, ref_seq):
     def cigar_str(self, ref_seq):
         if self._cigar_str is not None:
             return self._cigar_str
-        alignment = Sequence.aligner.align(
+        alignment = Aligner().align(
             ref_seq[self.start:self.start + len(self)],
             str(self),
-        )[0]
-        cigar_iter = self.assemble_cigar_list(alignment[0], alignment[1])
+        )
+        cigar_iter = self.assemble_cigar_list(*alignment)
         cigar_pieces = [f'{length}{op}' for length, op in cigar_iter]
         self._cigar_str = ''.join(cigar_pieces)
         return self._cigar_str