diff --git a/.github/workflows/unittest.yaml b/.github/workflows/unittest.yaml index a32f53d..e3e2400 100644 --- a/.github/workflows/unittest.yaml +++ b/.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 diff --git a/pyproject.toml b/pyproject.toml index 8368496..fc553ca 100644 --- a/pyproject.toml +++ b/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"] diff --git a/test/aligner.py b/test/aligner.py new file mode 100644 index 0000000..3ad036b --- /dev/null +++ b/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 diff --git a/test/sam_gen.py b/test/sam_gen.py index c456766..cce9de2 100644 --- a/test/sam_gen.py +++ b/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