From 056124d400c8decbe84d780e680e1a7688cefe10 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 13 Mar 2026 16:22:17 +0100
Subject: [PATCH 01/25] added function for pairwise distance compared to
 reference

---
 msaexplorer/explore.py | 193 +++++++++++++++++++++++++++--------------
 1 file changed, 129 insertions(+), 64 deletions(-)

diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 32e6d0a..2d9a176 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -173,9 +173,92 @@ def _determine_aln_type(alignment) -> str:
         else:
             return 'AA'
 
+    @staticmethod
+    def _create_distance_calculation_function_mapping() -> Dict[str, Callable[[str, str, int], float]]:
+        """
+        create a mapping of distance types to distance calculation functions
+        :return: dictionary of mappings
+        """
+
+        def ghd(seq1: str, seq2: str, aln_length: int) -> float:
+            """
+            global hamming distance - defined as percentage of total number of matches
+            """
+            return sum(c1 == c2 for c1, c2 in zip(seq1, seq2)) / aln_length * 100
+
+        def lhd(seq1: str, seq2: str, aln_length: int) -> float:
+            """
+            local hamming distance - defined as total number of matches excluding terminal gaps
+            """
+            # Trim gaps from both sides
+            i, j = 0, aln_length - 1
+            while i < aln_length and (seq1[i] == '-' or seq2[i] == '-'):
+                i += 1
+            while j >= 0 and (seq1[j] == '-' or seq2[j] == '-'):
+                j -= 1
+            if i > j:
+                return 0.0
+
+            seq1_, seq2_ = seq1[i:j + 1], seq2[i:j + 1]
+            matches = sum(c1 == c2 for c1, c2 in zip(seq1_, seq2_))
+            length = j - i + 1
+
+            return (matches / length) * 100 if length > 0 else 0.0
+
+        def ged(seq1: str, seq2: str, aln_length: int = None) -> float:
+            """
+            gap excluded distance - defined as percentage of total number of matches excluding all gaps
+            """
+
+            matches, mismatches = 0, 0
+
+            for c1, c2 in zip(seq1, seq2):
+                if c1 != '-' and c2 != '-':
+                    if c1 == c2:
+                        matches += 1
+                    else:
+                        mismatches += 1
+            return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
+
+        def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
+            """
+            gap compressed distance - defined as percentage of total number of matches with sequential gap mismatches
+            counting as a single mismatch
+            """
+            matches = 0
+            mismatches = 0
+            in_gap = False
+
+            for char1, char2 in zip(seq1, seq2):
+                if char1 == '-' and char2 == '-':  # Shared gap: do nothing
+                    continue
+                elif char1 == '-' or char2 == '-':  # Gap in only one sequence
+                    if not in_gap:  # Start of a new gap stretch
+                        mismatches += 1
+                        in_gap = True
+                else:  # No gaps
+                    in_gap = False
+                    if char1 == char2:  # Matching characters
+                        matches += 1
+                    else:  # Mismatched characters
+                        mismatches += 1
+
+            return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
+
+
+        # Map distance type to corresponding function
+        distance_functions: Dict[str, Callable[[str, str, int], float]] = {
+            'ghd': ghd,
+            'lhd': lhd,
+            'ged': ged,
+            'gcd': gcd
+        }
+
+        return distance_functions
+
     # Properties with setters
     @property
-    def reference_id(self):
+    def reference_id(self) -> str:
         return self._reference_id
 
     @reference_id.setter
@@ -1052,71 +1135,11 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
         **4) gcd (gap compressed distance)**: All consecutive gaps are compressed to one mismatch.
         \ndistance = matches / gap_compressed_alignment_length * 100
 
+        :param distance_type: type of distance computation technique
         :return: array with pairwise distances.
         """
 
-        def hamming_distance(seq1: str, seq2: str) -> int:
-            return sum(c1 == c2 for c1, c2 in zip(seq1, seq2))
-
-        def ghd(seq1: str, seq2: str) -> float:
-            return hamming_distance(seq1, seq2) / self.length * 100
-
-        def lhd(seq1, seq2):
-            # Trim gaps from both sides
-            i, j = 0, self.length - 1
-            while i < self.length and (seq1[i] == '-' or seq2[i] == '-'):
-                i += 1
-            while j >= 0 and (seq1[j] == '-' or seq2[j] == '-'):
-                j -= 1
-            if i > j:
-                return 0.0
-
-            seq1_, seq2_ = seq1[i:j + 1], seq2[i:j + 1]
-            matches = sum(c1 == c2 for c1, c2 in zip(seq1_, seq2_))
-            length = j - i + 1
-            return (matches / length) * 100 if length > 0 else 0.0
-
-        def ged(seq1: str, seq2: str) -> float:
-
-            matches, mismatches = 0, 0
-
-            for c1, c2 in zip(seq1, seq2):
-                if c1 != '-' and c2 != '-':
-                    if c1 == c2:
-                        matches += 1
-                    else:
-                        mismatches += 1
-            return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
-
-        def gcd(seq1: str, seq2: str) -> float:
-            matches = 0
-            mismatches = 0
-            in_gap = False
-
-            for char1, char2 in zip(seq1, seq2):
-                if char1 == '-' and char2 == '-':  # Shared gap: do nothing
-                    continue
-                elif char1 == '-' or char2 == '-':  # Gap in only one sequence
-                    if not in_gap:  # Start of a new gap stretch
-                        mismatches += 1
-                        in_gap = True
-                else:  # No gaps
-                    in_gap = False
-                    if char1 == char2:  # Matching characters
-                        matches += 1
-                    else:  # Mismatched characters
-                        mismatches += 1
-
-            return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
-
-
-        # Map distance type to corresponding function
-        distance_functions: Dict[str, Callable[[str, str], float]] = {
-            'ghd': ghd,
-            'lhd': lhd,
-            'ged': ged,
-            'gcd': gcd
-        }
+        distance_functions = self._create_distance_calculation_function_mapping()
 
         if distance_type not in distance_functions:
             raise ValueError(f"Invalid distance type '{distance_type}'. Choose from {list(distance_functions.keys())}.")
@@ -1132,12 +1155,54 @@ def gcd(seq1: str, seq2: str) -> float:
             seq1 = sequences[i]
             for j in range(i, n):
                 seq2 = sequences[j]
-                dist = distance_func(seq1, seq2)
+                dist = distance_func(seq1, seq2, self.length)
                 distance_matrix[i, j] = dist
                 distance_matrix[j, i] = dist
 
         return distance_matrix
 
+    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> tuple[str, list, ndarray]:
+        """
+        Calculate pairwise identities between reference and all sequences in the alignment. Same computation as calc_pairwise_identity_matrix but compared to a single sequence. Supported distance computation methods.
+
+        **1) ghd (global hamming distance)**: At each alignment position, check if characters match:
+        \ndistance = matches / alignment_length * 100
+
+        **2) lhd (local hamming distance)**: Restrict the alignment to the region in both sequences that do not start and end with gaps:
+        \ndistance = matches / min(5'3' ungapped seq1, 5'3' ungapped seq2) * 100
+
+        **3) ged (gap excluded distance)**: All gaps are excluded from the alignment
+        \ndistance = matches / (matches + mismatches) * 100
+
+        **4) gcd (gap compressed distance)**: All consecutive gaps are compressed to one mismatch.
+        \ndistance = matches / gap_compressed_alignment_length * 100
+
+        :param distance_type: type of distance computation technique
+        :return: tuple with reference id, sequence ids and pairwise distances.
+        """
+
+        distance_functions = self._create_distance_calculation_function_mapping()
+
+        if distance_type not in distance_functions:
+            raise ValueError(f"Invalid distance type '{distance_type}'. Choose from {list(distance_functions.keys())}.")
+
+        distance_func = distance_functions[distance_type]
+        aln = self.alignment
+        ref_id = self.reference_id
+
+        ref_seq = aln[ref_id] if ref_id is not None else self.get_consensus()
+        distances = []
+        distance_names = []
+
+        for seq_id in aln:
+            if seq_id == ref_id:
+                continue
+            distance_names.append(seq_id)
+            distances.append(distance_func(ref_seq, aln[seq_id], self.length))
+
+        return ref_id if ref_id is not None else 'consensus', distance_names, np.array(distances)
+
+
     def get_snps(self, include_ambig:bool=False) -> dict:
         """
         Calculate snps similar to snp-sites (output is comparable):

From ea61e3d9fc5fbae09edacebe810fdcc5f4ed3329 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 13 Mar 2026 16:50:20 +0100
Subject: [PATCH 02/25] planned simplot

---
 msaexplorer/draw.py              | 23 ++++++++++++++++++-
 msaexplorer/explore.py           | 32 ++++++++++++++++++++++----
 tests/test_stats_calculations.py | 39 +++++++++++++++++++++++++++++++-
 3 files changed, 88 insertions(+), 6 deletions(-)

diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 35ff57f..2a8233c 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -1198,4 +1198,25 @@ def consensus_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, threshold
     else:
         ax.set_yticks([])
 
-    return ax
\ No newline at end of file
+    return ax
+
+def simplot(aln: explore.MSA | str, ref: str | None, colors: str | list,  window_size: int = 50,
+            ax: plt.Axes | None = None, show_x_label: bool = False) -> plt.Axes:
+    """
+    Calculate binned pairwise distances (similarity) between sequences and a reference sequence
+    over a sliding window in the zoomed region of the alignment. Multiple distance
+    calculation options are supported. Each sequence is plotted as a line displaying the similarity
+    to the reference. The reference sequence can be either set to None (here it will be calculated to a
+    consensus) or to a specific reference id.
+    a majority consensus.
+
+    :param aln: alignment MSA class or path
+    :param ref: reference sequence id or None
+    :param colors: color for each sequence. can be a single color or a list of colors
+    :param window_size: window size for sliding window
+    :param ax: matplotlib axes
+    :param show_x_label: whether to show the x-axis label
+    :return:
+    """
+
+    pass
\ No newline at end of file
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 2d9a176..b183eac 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -14,6 +14,7 @@
 import math
 import collections
 import re
+from dataclasses import dataclass
 from typing import Callable, Dict
 
 # installed
@@ -37,6 +38,26 @@ def _get_line_iterator(source):
     else:
         return io.StringIO(source)
 
+
+@dataclass(frozen=True)
+class PairwiseDistanceToReferenceResult:
+    """
+    Result container for pairwise identity values between a reference/consensus
+    sequence and each sequence in the alignment.
+
+    The object remains iterable so it can be unpacked as a tuple:
+    ``reference_label, sequence_ids, distances = result``.
+    """
+
+    reference_label: str
+    sequence_ids: list[str]
+    distances: ndarray
+
+    def __iter__(self):
+        yield self.reference_label
+        yield self.sequence_ids
+        yield self.distances
+
 class MSA:
     """
     An alignment class that allows computation of several stats. Supported inputs are file paths to alignments in "fasta",
@@ -1161,7 +1182,7 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
 
         return distance_matrix
 
-    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> tuple[str, list, ndarray]:
+    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> PairwiseDistanceToReferenceResult:
         """
         Calculate pairwise identities between reference and all sequences in the alignment. Same computation as calc_pairwise_identity_matrix but compared to a single sequence. Supported distance computation methods.
 
@@ -1178,7 +1199,7 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> tuple[
         \ndistance = matches / gap_compressed_alignment_length * 100
 
         :param distance_type: type of distance computation technique
-        :return: tuple with reference id, sequence ids and pairwise distances.
+        :return: dataclass with reference label, sequence ids and pairwise distances.
         """
 
         distance_functions = self._create_distance_calculation_function_mapping()
@@ -1200,8 +1221,11 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> tuple[
             distance_names.append(seq_id)
             distances.append(distance_func(ref_seq, aln[seq_id], self.length))
 
-        return ref_id if ref_id is not None else 'consensus', distance_names, np.array(distances)
-
+        return PairwiseDistanceToReferenceResult(
+            reference_label=ref_id if ref_id is not None else 'consensus',
+            sequence_ids=distance_names,
+            distances=np.array(distances)
+        )
 
     def get_snps(self, include_ambig:bool=False) -> dict:
         """
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index 97cc3f4..ce16aeb 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -2,7 +2,7 @@
 
 import pytest
 from conftest import create_alignment
-from msaexplorer.explore import MSA
+from msaexplorer.explore import MSA, PairwiseDistanceToReferenceResult
 import numpy as np
 
 
@@ -295,6 +295,43 @@ def test_invalid_distance_type_raises(self):
             msa.calc_pairwise_identity_matrix(distance_type="invalid")
 
 
+class TestCalcPairwiseDistanceToReference:
+    """Tests for calc_pairwise_distance_to_reference."""
+
+    @pytest.mark.parametrize(
+        "sequences",
+        [
+            {"ref": "AAAA", "q1": "AAAT", "q2": "AATT"},
+            {"q1": "AAAT", "ref": "AAAA", "q2": "AATT"},
+            {"q1": "AAAT", "q2": "AATT", "ref": "AAAA"},
+        ],
+    )
+    def test_returns_dataclass_for_all_reference_positions(self, sequences):
+        msa = MSA(create_alignment(sequences), reference_id="ref")
+
+        result = msa.calc_pairwise_distance_to_reference(distance_type="ghd")
+
+        assert isinstance(result, PairwiseDistanceToReferenceResult)
+        assert result.reference_label == "ref"
+        assert result.sequence_ids == ["q1", "q2"]
+        assert np.allclose(result.distances, np.array([75.0, 50.0]))
+
+    def test_uses_consensus_when_reference_id_is_not_set(self):
+        msa = MSA(create_alignment({"q1": "AAAA", "q2": "AAAT", "q3": "AATT"}))
+
+        result = msa.calc_pairwise_distance_to_reference(distance_type="ghd")
+
+        assert result.reference_label == "consensus"
+        assert result.sequence_ids == ["q1", "q2", "q3"]
+        assert np.allclose(result.distances, np.array([75.0, 100.0, 75.0]))
+
+    def test_invalid_distance_type_raises(self):
+        msa = MSA(create_alignment({"s1": "ACGT", "s2": "ACGT"}))
+
+        with pytest.raises(ValueError, match="Invalid distance type"):
+            msa.calc_pairwise_distance_to_reference(distance_type="invalid")
+
+
 class TestCalcPositionMatrix:
     """Tests for calc_position_matrix."""
 

From 178e56f692374f1def53584cd7b413679fd372ec Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 15 Mar 2026 13:30:52 +0100
Subject: [PATCH 03/25] added simplot and new pairwise distance calculations

---
 README.md                        |  26 +++---
 msaexplorer/draw.py              | 128 ++++++++++++++++++++++----
 msaexplorer/explore.py           | 152 ++++++++++++++++++++++++-------
 tests/test_stats_calculations.py |   8 +-
 4 files changed, 246 insertions(+), 68 deletions(-)

diff --git a/README.md b/README.md
index a494637..3c3aa90 100644
--- a/README.md
+++ b/README.md
@@ -87,15 +87,15 @@ options:
   --version    show program's version number and exit  
 ```
 
-- :white_check_mark: The app runs solely in your browser. No need to install anything, just might take a few seconds to load.
-- :white_check_mark: Use the app offline (after loading it).
-- :white_check_mark: Analyse alignments on your smartphone or tablet.
-- :white_check_mark: Download alignment statistics (e.g. entropy, SNPs, coverage, consensus, ORFs and more).
-- :white_check_mark: Annotate the alignment by additionally reading in gb, gff or bed files.
-- :white_check_mark: Flexibility to customize plots and colors.
-- :white_check_mark: Easily export the plot as pdf.
-- :white_check_mark: Generate plots of the whole alignment as well as just parts of it.
-- :white_check_mark: Publication ready figures with just a few clicks.
+- The app runs solely in your browser. No need to install anything, just might take a few seconds to load.
+- Use the app offline (after loading it).
+- Analyse alignments on your smartphone or tablet.
+- Download alignment statistics (e.g. entropy, SNPs, coverage, consensus, ORFs and more).
+- Annotate the alignment by additionally reading in gb, gff or bed files.
+- Flexibility to customize plots and colors.
+- Easily export the plot as pdf.
+- Generate plots of the whole alignment as well as just parts of it.
+- Publication ready figures with just a few clicks.
 
 | ![](readme_assets/upload_tab.png) | ![](readme_assets/plot_tab.png) | ![](readme_assets/plot2_tab.png) | ![](readme_assets/analysis_tab.png) |
 |-----------------------------------|---------------------------------|----------------------------------|-------------------------------------|
@@ -112,11 +112,9 @@ shinylive export ./ site/  # you should now have a new 'site' folder with the ap
 ```
 
 ## Features of MSAexplorer as a python package ([full documentation](https://jonas-fuchs.github.io/MSAexplorer/docs/msaexplorer.html))
-- :white_check_mark: Access MSAexplorer as a python package
-- :white_check_mark: Seamlessly integrates with Biopython.
-- :white_check_mark: Maximum flexibility for the plotting and analysis features while retaining minimal syntax.
-- :white_check_mark: Integrates seamlessly with matplotlib.
-- :white_check_mark: Minimal requirements.
+- Access MSAexplorer as a python package
+- Seamlessly integrates with Biopython and matplotlib.
+- Maximum flexibility for the plotting and analysis features while retaining minimal syntax.
 
 ```python
 ### Minimal analysis example ###
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 2a8233c..52dcb5c 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -10,21 +10,18 @@
 ## Functions
 
 """
-import pathlib
 # built-in
 from itertools import chain
 from typing import Callable, Dict
 from copy import deepcopy
 import os
 
-import matplotlib
-from numpy import ndarray
-
 # MSAexplorer
 from msaexplorer import explore, config
 
 # libs
 import numpy as np
+from numpy import ndarray
 import matplotlib.pyplot as plt
 import matplotlib.patches as patches
 from matplotlib.cm import ScalarMappable
@@ -1200,23 +1197,122 @@ def consensus_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, threshold
 
     return ax
 
-def simplot(aln: explore.MSA | str, ref: str | None, colors: str | list,  window_size: int = 50,
-            ax: plt.Axes | None = None, show_x_label: bool = False) -> plt.Axes:
+def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None, colors: str | list | None = None,
+            window_size: int = 200, step_size: int = 20, distance_calculation: str = 'ghd', line_width: int | float = 0.5,
+            show_legend: bool = False, bbox_to_anchor: tuple[float|int, float|int] | list= (1, 1), show_x_label: bool = False) -> plt.Axes:
     """
     Calculate binned pairwise distances (similarity) between sequences and a reference sequence
-    over a sliding window in the zoomed region of the alignment. Multiple distance
-    calculation options are supported. Each sequence is plotted as a line displaying the similarity
-    to the reference. The reference sequence can be either set to None (here it will be calculated to a
-    consensus) or to a specific reference id.
-    a majority consensus.
+    over a stepwise sliding window in the zoomed region of the alignment. This is inspired by simplot that helps to identify
+    recombination sites. For proper recombination analyses use simplot or simplot++ (https://github.com/Stephane-S/Simplot_PlusPlus).
+    Each sequence is plotted as a line displaying the similarity to the reference. The reference sequence can be either
+    set to None (compared to consensus) or to a specific reference id.
 
     :param aln: alignment MSA class or path
-    :param ref: reference sequence id or None
-    :param colors: color for each sequence. can be a single color or a list of colors
-    :param window_size: window size for sliding window
     :param ax: matplotlib axes
+    :param ref: reference sequence id or None. For None all computations are compared to a majority consensus
+    :param colors: color for each sequence. can be a single named color or a list of named colors or a plt.colormap or None (auto coloring)
+    :param window_size: window size for sliding window
+    :param step_size: step size for sliding window
+    :param distance_calculation: distance calculation method. Supported: ghd (global hamming distance), ged (gap excluded distance) and for nt: jc69(Jukes-Cantor 1969) and k2p (Kimura 2-Parameter / K80). For more information see: explore.MSA.calc_pairwise_identity_matrix()
+    :param line_width: width of the plotted lines
+    :param show_legend: whether to show the legend
+    :param bbox_to_anchor: bounding box coordinates for the legend - see: https://matplotlib.org/stable/api/legend_api.html
     :param show_x_label: whether to show the x-axis label
-    :return:
+
+    :return: matplotlib axes
     """
 
-    pass
\ No newline at end of file
+    # validate inputs
+    aln, ax = _validate_input_parameters(aln=aln, ax=ax)
+    if not isinstance(window_size, int) or window_size <= 0:
+        raise ValueError('window_size has to be a positive integer')
+    if window_size > aln.length:
+        raise ValueError('window_size can not be larger than the (zoomed) alignment length')
+    if ref is not None and ref not in aln.alignment:
+        raise ValueError(f'Reference {ref} not in alignment')
+    if distance_calculation not in ['ghd', 'ged', 'jc69', 'k2p']:
+        raise ValueError(f'Distance calculation method {distance_calculation} not supported. Supported: ghd, ged, jc69, k2p')
+    if distance_calculation in ['jc69', 'k2p'] and aln.aln_type == 'AA':
+        raise ValueError(f'Distance calculation method {distance_calculation} only supported for nucleotide alignments')
+
+    # get the sequence ids to plot
+    sequence_ids = [key for key in aln.alignment.keys() if key != ref]
+
+    # validate colors
+    if colors is not None:
+        # named color or colormap
+        if isinstance(colors, str):
+            # first try potential colormap
+            try:
+                cmap = plt.get_cmap(colors)
+                cmap_colors = cmap(np.linspace(0, 1, len(sequence_ids)))
+                color_map = {seq_id: color for seq_id, color in zip(sequence_ids, cmap_colors)}
+            # single color
+            except ValueError:
+                _validate_color(colors)
+                color_map = {seq_id: colors for seq_id in sequence_ids}
+        # list of colors
+        elif isinstance(colors, list):
+            if len(colors) != len(sequence_ids):
+                raise ValueError('colors list length has to match the number of plotted sequences')
+            for color in colors:
+                _validate_color(color)
+            color_map = {seq_id: color for seq_id, color in zip(sequence_ids, colors)}
+        else:
+            raise ValueError('colors has to be either a single named color, a list of named colors, a plt colormap or None (auto)')
+    else:
+        color_map = None
+
+    # define region and window
+    region_start = aln.zoom[0] if aln.zoom is not None else 0
+    half_window_size = window_size // 2
+
+    # copy alignment and set reference
+    aln_tmp = deepcopy(aln)
+    aln_tmp.reference_id = ref
+    # reset zoom (later the alignment dictionary will be replaced for a slice of the alignment)
+    aln_tmp._zoom = None
+
+    # remember x positions and distances
+    x_positions = []
+    distance_traces = {seq_id: [] for seq_id in sequence_ids}
+
+    for point_to_plot in range(0, aln.length + 1, step_size):
+        # define windows
+        left_side = point_to_plot - half_window_size if point_to_plot - half_window_size > 0 else 0
+        right_side = point_to_plot + half_window_size if point_to_plot + half_window_size < aln.length else aln.length
+        # slice alignment and replace the original alignment
+        aln_tmp._alignment = {
+            seq_id: seq[left_side:right_side]
+            for seq_id, seq in aln.alignment.items()
+        }
+        window_result = aln_tmp.calc_pairwise_distance_to_reference(distance_type=distance_calculation)
+        value_map = {seq_id: value for seq_id, value in zip(window_result.sequence_ids, window_result.distances)}
+
+        x_positions.append(region_start + point_to_plot)
+        for seq_id in sequence_ids:
+            distance_traces[seq_id].append(value_map[seq_id])
+
+    for seq_id in sequence_ids:
+        if color_map is not None:
+            ax.plot(x_positions, distance_traces[seq_id],
+                    color=color_map[seq_id],
+                    linewidth=line_width, label=seq_id)
+        else:
+            ax.plot(x_positions, distance_traces[seq_id],
+                    linewidth=line_width, label=seq_id)
+    # Format axis
+    ax.set_ylabel('similarity (%)')
+    ax.set_ylim(0, 102)
+    _format_x_axis(aln=aln, ax=ax, show_x_label=show_x_label, show_left=True)
+    # add legend for all sequences on top
+    if show_legend:
+        leg1 = ax.legend(frameon=False, loc='lower right', bbox_to_anchor=bbox_to_anchor, ncols=3)
+        ax.add_artist(leg1)
+    # add legend for query
+    ref_label = ref if ref is not None else 'consensus'
+    ref_handle = plt.Line2D([0], [0], linewidth=0, label=f'query sequence: {ref_label}')
+    leg2 = ax.legend(handles=[ref_handle], frameon=False, bbox_to_anchor=(1, 0), loc='lower right')
+    ax.add_artist(leg2)
+
+    return ax
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index b183eac..0a1a425 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -39,25 +39,6 @@ def _get_line_iterator(source):
         return io.StringIO(source)
 
 
-@dataclass(frozen=True)
-class PairwiseDistanceToReferenceResult:
-    """
-    Result container for pairwise identity values between a reference/consensus
-    sequence and each sequence in the alignment.
-
-    The object remains iterable so it can be unpacked as a tuple:
-    ``reference_label, sequence_ids, distances = result``.
-    """
-
-    reference_label: str
-    sequence_ids: list[str]
-    distances: ndarray
-
-    def __iter__(self):
-        yield self.reference_label
-        yield self.sequence_ids
-        yield self.distances
-
 class MSA:
     """
     An alignment class that allows computation of several stats. Supported inputs are file paths to alignments in "fasta",
@@ -77,6 +58,25 @@ def __init__(self, alignment_string: str | MultipleSeqAlignment, reference_id: s
         self._zoom = self._validate_zoom(zoom_range, self._alignment)
         self._aln_type = self._determine_aln_type(self._alignment)
 
+    @dataclass(frozen=True)
+    class PairwiseDistanceResult:
+        """
+        Result container for pairwise identity values between a reference/consensus
+        sequence and each sequence in the alignment.
+
+        The object remains iterable so it can be unpacked as a tuple:
+        ``reference_label, sequence_ids, distances = result``.
+        """
+
+        reference_id: str
+        sequence_ids: list[str]
+        distances: ndarray
+
+        def __iter__(self):
+            yield self.reference_id
+            yield self.sequence_ids
+            yield self.distances
+
     # Static methods
     @staticmethod
     def _read_alignment(source: str | MultipleSeqAlignment) -> dict:
@@ -231,15 +231,15 @@ def ged(seq1: str, seq2: str, aln_length: int = None) -> float:
             gap excluded distance - defined as percentage of total number of matches excluding all gaps
             """
 
-            matches, mismatches = 0, 0
+            diff, total = 0, 0
 
             for c1, c2 in zip(seq1, seq2):
                 if c1 != '-' and c2 != '-':
-                    if c1 == c2:
-                        matches += 1
-                    else:
-                        mismatches += 1
-            return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
+                    total += 1
+                    if c1 != c2:
+                        diff += 1
+
+            return (1 - diff / total) * 100 if total > 0 else 0
 
         def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
             """
@@ -266,13 +266,75 @@ def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
 
             return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
 
+        def jc69(seq1: str, seq2: str, aln_length: int = None) -> float:
+            """
+            Jukes-Cantor 1969 (JC69) corrected identity.
+            Gaps are excluded. The proportion of differing sites (p-distance) is corrected
+            for multiple hits: d = -(3/4) * ln(1 - (4/3) * p).
+            Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
+            Returns 0 when p >= 0.75 (formula undefined / sequence saturated).
+            """
+            diff, total = 0, 0
+            for c1, c2 in zip(seq1, seq2):
+                if c1 != '-' and c2 != '-':
+                    total += 1
+                    if c1 != c2:
+                        diff += 1
+            if total == 0:
+                return 0.0
+            p = diff / total
+            if p == 0.0:
+                return 100.0
+            correction = 1.0 - (4.0 / 3.0) * p
+            if correction <= 0.0:  # saturated – formula undefined
+                return 0.0
+            d = -(3.0 / 4.0) * math.log(correction)
+            return max(0.0, (1.0 - d) * 100.0)
+
+        def k2p(seq1: str, seq2: str, aln_length: int = None) -> float:
+            """
+            Kimura 2-Parameter (K2P / K80) corrected identity.
+            Gaps are excluded. Transitions and transversions (Tv) are weighted separately:
+            d = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q)
+            where P = Ti / total and Q = Tv / total.
+            Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
+            Returns 0 when the logarithm arguments become non-positive (saturated).
+            """
+            transitions = [{'A', 'G'}, {'C', 'T'}]
+
+            ts, tv, total = 0, 0, 0
+            for c1, c2 in zip(seq1, seq2):
+                if c1 != '-' and c2 != '-':
+                    total += 1
+                    if c1 != c2:
+                        if {c1, c2} in transitions:
+                            ts += 1
+                        else:
+                            tv += 1
+            if total == 0:
+                return 0.0
+            if ts == 0 and tv == 0:
+                return 100.0
+            P = ts / total   # transition proportion
+            Q = tv / total   # transversion proportion
+            term1 = 1.0 - 2.0 * P - Q
+            term2 = 1.0 - 2.0 * Q
+            # saturated – formula undefined
+            if term1 <= 0.0 or term2 <= 0.0:
+                return 0.0
+            # calculate distance
+            d = -0.5 * math.log(term1) - 0.25 * math.log(term2)
+
+            return max(0.0, (1 - d) * 100.0)
 
         # Map distance type to corresponding function
         distance_functions: Dict[str, Callable[[str, str, int], float]] = {
             'ghd': ghd,
             'lhd': lhd,
             'ged': ged,
-            'gcd': gcd
+            'gcd': gcd,
+            'jc69': jc69,
+            'k2p': k2p,
         }
 
         return distance_functions
@@ -831,7 +893,7 @@ def calc_reverse_complement_alignment(self) -> dict | TypeError:
 
         return reverse_complement_dict
 
-    def calc_numerical_alignment(self, encode_mask:bool=False, encode_ambiguities:bool=False):
+    def calc_numerical_alignment(self, encode_mask:bool=False, encode_ambiguities:bool=False) -> ndarray:
         """
         Transforms the alignment to numerical values. Ambiguities are encoded as -3, mask as -2 and the
         remaining chars with the idx + 1 of config.CHAR_COLORS[self.aln_type]['standard'].
@@ -861,7 +923,7 @@ def calc_numerical_alignment(self, encode_mask:bool=False, encode_ambiguities:bo
 
         return numerical_matrix
 
-    def calc_identity_alignment(self, encode_mismatches:bool=True, encode_mask:bool=False, encode_gaps:bool=True, encode_ambiguities:bool=False, encode_each_mismatch_char:bool=False) -> np.ndarray:
+    def calc_identity_alignment(self, encode_mismatches:bool=True, encode_mask:bool=False, encode_gaps:bool=True, encode_ambiguities:bool=False, encode_each_mismatch_char:bool=False) -> ndarray:
         """
         Converts alignment to identity array (identical=0) compared to majority consensus or reference:\n
 
@@ -922,7 +984,7 @@ def calc_identity_alignment(self, encode_mismatches:bool=True, encode_mask:bool=
 
         return identity_matrix
 
-    def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=True) -> np.ndarray:
+    def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=True) -> ndarray:
         """
         Calculate the similarity score between the alignment and the reference sequence, with normalization to highlight
         differences. The similarity scores are scaled to the range [0, 1] based on the substitution matrix values for the
@@ -995,7 +1057,7 @@ def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=Tr
 
         return similarity_array
 
-    def calc_position_matrix(self, matrix_type:str='PWM') -> np.ndarray | ValueError:
+    def calc_position_matrix(self, matrix_type:str='PWM') -> ndarray | ValueError:
         """
         Calculates a position matrix of the specified type for the given alignment. The function
         supports generating matrices of types Position Frequency Matrix (PFM), Position Probability
@@ -1156,6 +1218,14 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
         **4) gcd (gap compressed distance)**: All consecutive gaps are compressed to one mismatch.
         \ndistance = matches / gap_compressed_alignment_length * 100
 
+        **5) jc69 (Jukes-Cantor 1969)**: Gaps excluded. Applies the JC69 substitution model to correct
+        the p-distance for multiple hits (assumes equal base frequencies and substitution rates).
+        \ncorrected_identity = (1 - d_JC69) * 100,  where d = -(3/4) * ln(1 - (4/3) * p)
+
+        **6) k2p (Kimura 2-Parameter / K80)**: Gaps excluded. Distinguishes transitions (Ti) and
+        transversions (Tv). Returns (1 - d_K2P) * 100 as corrected percent identity.
+        \nd = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q),  P = Ti/total, Q = Tv/total
+
         :param distance_type: type of distance computation technique
         :return: array with pairwise distances.
         """
@@ -1165,8 +1235,12 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
         if distance_type not in distance_functions:
             raise ValueError(f"Invalid distance type '{distance_type}'. Choose from {list(distance_functions.keys())}.")
 
-        # Compute pairwise distances
         aln = self.alignment
+
+        if self.aln_type == 'AA' and distance_type in ['jc69', 'k2p']:
+            raise ValueError(f"JC69 and K2P are not supported for {self.aln_type} alignment.")
+
+        # Compute pairwise distances
         distance_func = distance_functions[distance_type]
         distance_matrix = np.zeros((len(aln), len(aln)))
 
@@ -1182,7 +1256,7 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
 
         return distance_matrix
 
-    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> PairwiseDistanceToReferenceResult:
+    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> PairwiseDistanceResult:
         """
         Calculate pairwise identities between reference and all sequences in the alignment. Same computation as calc_pairwise_identity_matrix but compared to a single sequence. Supported distance computation methods.
 
@@ -1198,6 +1272,12 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
         **4) gcd (gap compressed distance)**: All consecutive gaps are compressed to one mismatch.
         \ndistance = matches / gap_compressed_alignment_length * 100
 
+        **5) jc69 (Jukes-Cantor 1969)**: Gaps excluded. JC69 substitution-model corrected identity.
+        \ncorrected_identity = (1 - d_JC69) * 100
+
+        **6) k2p (Kimura 2-Parameter / K80)**: Gaps excluded. Distinguishes transitions and transversions.
+        \ncorrected_identity = (1 - d_K2P) * 100
+
         :param distance_type: type of distance computation technique
         :return: dataclass with reference label, sequence ids and pairwise distances.
         """
@@ -1207,6 +1287,10 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
         if distance_type not in distance_functions:
             raise ValueError(f"Invalid distance type '{distance_type}'. Choose from {list(distance_functions.keys())}.")
 
+        if self.aln_type == 'AA' and distance_type in ['jc69', 'k2p']:
+            raise ValueError(f"JC69 and K2P are not supported for {self.aln_type} alignment.")
+
+        # Compute pairwise distances
         distance_func = distance_functions[distance_type]
         aln = self.alignment
         ref_id = self.reference_id
@@ -1221,8 +1305,8 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
             distance_names.append(seq_id)
             distances.append(distance_func(ref_seq, aln[seq_id], self.length))
 
-        return PairwiseDistanceToReferenceResult(
-            reference_label=ref_id if ref_id is not None else 'consensus',
+        return self.PairwiseDistanceResult(
+            reference_id=ref_id if ref_id is not None else 'consensus',
             sequence_ids=distance_names,
             distances=np.array(distances)
         )
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index ce16aeb..8ffd603 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -2,7 +2,7 @@
 
 import pytest
 from conftest import create_alignment
-from msaexplorer.explore import MSA, PairwiseDistanceToReferenceResult
+from msaexplorer.explore import MSA
 import numpy as np
 
 
@@ -311,8 +311,8 @@ def test_returns_dataclass_for_all_reference_positions(self, sequences):
 
         result = msa.calc_pairwise_distance_to_reference(distance_type="ghd")
 
-        assert isinstance(result, PairwiseDistanceToReferenceResult)
-        assert result.reference_label == "ref"
+        assert isinstance(result, MSA.PairwiseDistanceResult)
+        assert result.reference_id == "ref"
         assert result.sequence_ids == ["q1", "q2"]
         assert np.allclose(result.distances, np.array([75.0, 50.0]))
 
@@ -321,7 +321,7 @@ def test_uses_consensus_when_reference_id_is_not_set(self):
 
         result = msa.calc_pairwise_distance_to_reference(distance_type="ghd")
 
-        assert result.reference_label == "consensus"
+        assert result.reference_id == "consensus"
         assert result.sequence_ids == ["q1", "q2", "q3"]
         assert np.allclose(result.distances, np.array([75.0, 100.0, 75.0]))
 

From 76a3a83300316b42505d210c19208a57692dd3fb Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 15 Mar 2026 13:41:43 +0100
Subject: [PATCH 04/25] made refernce displaying id optional

---
 msaexplorer/draw.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 52dcb5c..c1733b0 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -1199,7 +1199,8 @@ def consensus_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, threshold
 
 def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None, colors: str | list | None = None,
             window_size: int = 200, step_size: int = 20, distance_calculation: str = 'ghd', line_width: int | float = 0.5,
-            show_legend: bool = False, bbox_to_anchor: tuple[float|int, float|int] | list= (1, 1), show_x_label: bool = False) -> plt.Axes:
+            show_legend: bool = False, show_reference: bool = True, bbox_to_anchor: tuple[float|int, float|int] | list= (1, 1),
+            show_x_label: bool = False) -> plt.Axes:
     """
     Calculate binned pairwise distances (similarity) between sequences and a reference sequence
     over a stepwise sliding window in the zoomed region of the alignment. This is inspired by simplot that helps to identify
@@ -1216,6 +1217,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
     :param distance_calculation: distance calculation method. Supported: ghd (global hamming distance), ged (gap excluded distance) and for nt: jc69(Jukes-Cantor 1969) and k2p (Kimura 2-Parameter / K80). For more information see: explore.MSA.calc_pairwise_identity_matrix()
     :param line_width: width of the plotted lines
     :param show_legend: whether to show the legend
+    :param show_reference: whether to show the reference id in the right lower corner
     :param bbox_to_anchor: bounding box coordinates for the legend - see: https://matplotlib.org/stable/api/legend_api.html
     :param show_x_label: whether to show the x-axis label
 
@@ -1310,9 +1312,10 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
         leg1 = ax.legend(frameon=False, loc='lower right', bbox_to_anchor=bbox_to_anchor, ncols=3)
         ax.add_artist(leg1)
     # add legend for query
-    ref_label = ref if ref is not None else 'consensus'
-    ref_handle = plt.Line2D([0], [0], linewidth=0, label=f'query sequence: {ref_label}')
-    leg2 = ax.legend(handles=[ref_handle], frameon=False, bbox_to_anchor=(1, 0), loc='lower right')
-    ax.add_artist(leg2)
+    if show_reference:
+        ref_label = ref if ref is not None else 'consensus'
+        ref_handle = plt.Line2D([0], [0], linewidth=0, label=f'reference: {ref_label}')
+        leg2 = ax.legend(handles=[ref_handle], frameon=False, bbox_to_anchor=(1, 0), loc='lower right')
+        ax.add_artist(leg2)
 
     return ax

From cb34f4a3187dab6b460c6b07c5d46f89f39f015e Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 15 Mar 2026 14:40:00 +0100
Subject: [PATCH 05/25] moved stuff

---
 msaexplorer/_helpers.py          | 191 +++++++++++++++++++++++++++
 msaexplorer/_msa_data_classes.py |  28 ++++
 msaexplorer/draw.py              |  23 +---
 msaexplorer/explore.py           | 220 ++++---------------------------
 msaexplorer/export.py            |  13 +-
 tests/test_stats_calculations.py |   5 +-
 6 files changed, 251 insertions(+), 229 deletions(-)
 create mode 100644 msaexplorer/_helpers.py
 create mode 100644 msaexplorer/_msa_data_classes.py

diff --git a/msaexplorer/_helpers.py b/msaexplorer/_helpers.py
new file mode 100644
index 0000000..4473c78
--- /dev/null
+++ b/msaexplorer/_helpers.py
@@ -0,0 +1,191 @@
+"""
+This contains helper functions, not intended to be used outside of this package.
+"""
+
+import os, io, math
+from typing import Callable, Dict
+from matplotlib.colors import is_color_like
+
+# explore helpers
+def _get_line_iterator(source):
+    """
+    allow reading in both raw string or paths
+    """
+    if isinstance(source, str) and os.path.exists(source):
+        return open(source, 'r')
+    else:
+        return io.StringIO(source)
+
+
+def _create_distance_calculation_function_mapping() -> Dict[str, Callable[[str, str, int], float]]:
+    """
+    create a mapping of distance types to distance calculation functions
+    :return: dictionary of mappings
+    """
+
+    def ghd(seq1: str, seq2: str, aln_length: int) -> float:
+        """
+        global hamming distance - defined as percentage of total number of matches
+        """
+        return sum(c1 == c2 for c1, c2 in zip(seq1, seq2)) / aln_length * 100
+
+    def lhd(seq1: str, seq2: str, aln_length: int) -> float:
+        """
+        local hamming distance - defined as total number of matches excluding terminal gaps
+        """
+        # Trim gaps from both sides
+        i, j = 0, aln_length - 1
+        while i < aln_length and (seq1[i] == '-' or seq2[i] == '-'):
+            i += 1
+        while j >= 0 and (seq1[j] == '-' or seq2[j] == '-'):
+            j -= 1
+        if i > j:
+            return 0.0
+
+        seq1_, seq2_ = seq1[i:j + 1], seq2[i:j + 1]
+        matches = sum(c1 == c2 for c1, c2 in zip(seq1_, seq2_))
+        length = j - i + 1
+
+        return (matches / length) * 100 if length > 0 else 0.0
+
+    def ged(seq1: str, seq2: str, aln_length: int = None) -> float:
+        """
+        gap excluded distance - defined as percentage of total number of matches excluding all gaps
+        """
+
+        diff, total = 0, 0
+
+        for c1, c2 in zip(seq1, seq2):
+            if c1 != '-' and c2 != '-':
+                total += 1
+                if c1 != c2:
+                    diff += 1
+
+        return (1 - diff / total) * 100 if total > 0 else 0
+
+    def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
+        """
+        gap compressed distance - defined as percentage of total number of matches with sequential gap mismatches
+        counting as a single mismatch
+        """
+        matches = 0
+        mismatches = 0
+        in_gap = False
+
+        for char1, char2 in zip(seq1, seq2):
+            if char1 == '-' and char2 == '-':  # Shared gap: do nothing
+                continue
+            elif char1 == '-' or char2 == '-':  # Gap in only one sequence
+                if not in_gap:  # Start of a new gap stretch
+                    mismatches += 1
+                    in_gap = True
+            else:  # No gaps
+                in_gap = False
+                if char1 == char2:  # Matching characters
+                    matches += 1
+                else:  # Mismatched characters
+                    mismatches += 1
+
+        return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
+
+    def jc69(seq1: str, seq2: str, aln_length: int = None) -> float:
+        """
+        Jukes-Cantor 1969 (JC69) corrected identity.
+        Gaps are excluded. The proportion of differing sites (p-distance) is corrected
+        for multiple hits: d = -(3/4) * ln(1 - (4/3) * p).
+        Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
+        Returns 0 when p >= 0.75 (formula undefined / sequence saturated).
+        """
+        diff, total = 0, 0
+        for c1, c2 in zip(seq1, seq2):
+            if c1 != '-' and c2 != '-':
+                total += 1
+                if c1 != c2:
+                    diff += 1
+        if total == 0:
+            return 0.0
+        p = diff / total
+        if p == 0.0:
+            return 100.0
+        correction = 1.0 - (4.0 / 3.0) * p
+        if correction <= 0.0:  # saturated – formula undefined
+            return 0.0
+        d = -(3.0 / 4.0) * math.log(correction)
+        return max(0.0, (1.0 - d) * 100.0)
+
+    def k2p(seq1: str, seq2: str, aln_length: int = None) -> float:
+        """
+        Kimura 2-Parameter (K2P / K80) corrected identity.
+        Gaps are excluded. Transitions and transversions (Tv) are weighted separately:
+        d = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q)
+        where P = Ti / total and Q = Tv / total.
+        Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
+        Returns 0 when the logarithm arguments become non-positive (saturated).
+        """
+        transitions = [{'A', 'G'}, {'C', 'T'}]
+
+        ts, tv, total = 0, 0, 0
+        for c1, c2 in zip(seq1, seq2):
+            if c1 != '-' and c2 != '-':
+                total += 1
+                if c1 != c2:
+                    if {c1, c2} in transitions:
+                        ts += 1
+                    else:
+                        tv += 1
+        if total == 0:
+            return 0.0
+        if ts == 0 and tv == 0:
+            return 100.0
+        P = ts / total   # transition proportion
+        Q = tv / total   # transversion proportion
+        term1 = 1.0 - 2.0 * P - Q
+        term2 = 1.0 - 2.0 * Q
+        # saturated – formula undefined
+        if term1 <= 0.0 or term2 <= 0.0:
+            return 0.0
+        # calculate distance
+        d = -0.5 * math.log(term1) - 0.25 * math.log(term2)
+
+        return max(0.0, (1 - d) * 100.0)
+
+    # Map distance type to corresponding function
+    distance_functions: Dict[str, Callable[[str, str, int], float]] = {
+        'ghd': ghd,
+        'lhd': lhd,
+        'ged': ged,
+        'gcd': gcd,
+        'jc69': jc69,
+        'k2p': k2p,
+    }
+
+    return distance_functions
+
+# export helpers
+def _check_and_create_path(path: str):
+    """
+    Check and create path if it doesn't exist.
+    :param path: string to file
+    """
+    if path is not None:
+        output_dir = os.path.dirname(path)
+        if output_dir and not os.path.exists(output_dir):
+            os.makedirs(output_dir)
+
+# draw helpers
+def _validate_color(c):
+    """
+    validate color and raise error
+    """
+    if not is_color_like(c):
+        raise ValueError(f'{c} is not a color')
+
+def _get_contrast_text_color(rgba_color):
+    """
+    compute the brightness of a color
+    """
+    r, g, b, a = rgba_color
+    brightness = (r * 299 + g * 587 + b * 114) / 1000
+
+    return 'white' if brightness < 0.5 else 'black'
+
diff --git a/msaexplorer/_msa_data_classes.py b/msaexplorer/_msa_data_classes.py
new file mode 100644
index 0000000..a91563d
--- /dev/null
+++ b/msaexplorer/_msa_data_classes.py
@@ -0,0 +1,28 @@
+"""
+this contains the dataclasses used to store the data for the msa explorer. these are not meant to be used outside of this package.
+"""
+
+# build-in
+from dataclasses import dataclass
+
+# libs
+from numpy import ndarray
+
+@dataclass(frozen=True)
+class PairwiseDistanceResult:
+    """
+    Result container for pairwise identity values between a reference/consensus
+    sequence and each sequence in the alignment.
+
+    The object remains iterable so it can be unpacked as a tuple:
+    ``reference_label, sequence_ids, distances = result``.
+    """
+
+    reference_id: str
+    sequence_ids: list[str]
+    distances: ndarray
+
+    def __iter__(self):
+        yield self.reference_id
+        yield self.sequence_ids
+        yield self.distances
\ No newline at end of file
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index c1733b0..58e9171 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -18,6 +18,7 @@
 
 # MSAexplorer
 from msaexplorer import explore, config
+from msaexplorer._helpers import _validate_color, _get_contrast_text_color
 
 # libs
 import numpy as np
@@ -25,7 +26,7 @@
 import matplotlib.pyplot as plt
 import matplotlib.patches as patches
 from matplotlib.cm import ScalarMappable
-from matplotlib.colors import is_color_like, Normalize, to_rgba, LinearSegmentedColormap
+from matplotlib.colors import Normalize, to_rgba, LinearSegmentedColormap
 from matplotlib.collections import PatchCollection, PolyCollection
 from matplotlib.text import TextPath
 from matplotlib.patches import PathPatch
@@ -72,14 +73,6 @@ def _validate_input_parameters(aln: explore.MSA | str, ax: plt.Axes, annotation:
         return aln, ax, annotation
 
 
-def _validate_color(c):
-    """
-    validate color and raise error
-    """
-    if not is_color_like(c):
-        raise ValueError(f'{c} is not a color')
-
-
 def _validate_color_scheme(scheme: str | None, aln: explore.MSA):
     """
     validates colorscheme
@@ -208,16 +201,6 @@ def _create_legend(color_scheme: str, aln_colors: dict, aln: explore.MSA, detect
     )
 
 
-def _get_contrast_text_color(rgba_color):
-        """
-        compute the brightness of a color
-        """
-        r, g, b, a = rgba_color
-        brightness = (r * 299 + g * 587 + b * 114) / 1000
-
-        return 'white' if brightness < 0.5 else 'black'
-
-
 def _create_alignment(aln: explore.MSA, ax: plt.Axes, matrix: ndarray, aln_colors: dict | ScalarMappable, fancy_gaps: bool, create_identity_patch: bool,
                       show_gaps: bool, show_different_sequence: bool, show_sequence_all: bool, reference_color: str | None, values_to_plot: list,
                       identical_value: int | float = 0):
@@ -1214,7 +1197,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
     :param colors: color for each sequence. can be a single named color or a list of named colors or a plt.colormap or None (auto coloring)
     :param window_size: window size for sliding window
     :param step_size: step size for sliding window
-    :param distance_calculation: distance calculation method. Supported: ghd (global hamming distance), ged (gap excluded distance) and for nt: jc69(Jukes-Cantor 1969) and k2p (Kimura 2-Parameter / K80). For more information see: explore.MSA.calc_pairwise_identity_matrix()
+    :param distance_calculation: distance calculation method. Supported: ghd (global hamming distance), ged (gap excluded distance) and for nt additionally: jc69(Jukes-Cantor 1969) and k2p (Kimura 2-Parameter / K80). For more information see: explore.MSA.calc_pairwise_identity_matrix()
     :param line_width: width of the plotted lines
     :param show_legend: whether to show the legend
     :param show_reference: whether to show the reference id in the right lower corner
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 0a1a425..580143a 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -9,12 +9,9 @@
 """
 
 # built-in
-import os
-import io
 import math
 import collections
 import re
-from dataclasses import dataclass
 from typing import Callable, Dict
 
 # installed
@@ -27,16 +24,8 @@
 
 # msaexplorer
 from msaexplorer import config
-
-
-def _get_line_iterator(source):
-    """
-    allow reading in both raw string or paths
-    """
-    if isinstance(source, str) and os.path.exists(source):
-        return open(source, 'r')
-    else:
-        return io.StringIO(source)
+from msaexplorer._msa_data_classes import PairwiseDistanceResult
+from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping
 
 
 class MSA:
@@ -58,25 +47,6 @@ def __init__(self, alignment_string: str | MultipleSeqAlignment, reference_id: s
         self._zoom = self._validate_zoom(zoom_range, self._alignment)
         self._aln_type = self._determine_aln_type(self._alignment)
 
-    @dataclass(frozen=True)
-    class PairwiseDistanceResult:
-        """
-        Result container for pairwise identity values between a reference/consensus
-        sequence and each sequence in the alignment.
-
-        The object remains iterable so it can be unpacked as a tuple:
-        ``reference_label, sequence_ids, distances = result``.
-        """
-
-        reference_id: str
-        sequence_ids: list[str]
-        distances: ndarray
-
-        def __iter__(self):
-            yield self.reference_id
-            yield self.sequence_ids
-            yield self.distances
-
     # Static methods
     @staticmethod
     def _read_alignment(source: str | MultipleSeqAlignment) -> dict:
@@ -194,151 +164,6 @@ def _determine_aln_type(alignment) -> str:
         else:
             return 'AA'
 
-    @staticmethod
-    def _create_distance_calculation_function_mapping() -> Dict[str, Callable[[str, str, int], float]]:
-        """
-        create a mapping of distance types to distance calculation functions
-        :return: dictionary of mappings
-        """
-
-        def ghd(seq1: str, seq2: str, aln_length: int) -> float:
-            """
-            global hamming distance - defined as percentage of total number of matches
-            """
-            return sum(c1 == c2 for c1, c2 in zip(seq1, seq2)) / aln_length * 100
-
-        def lhd(seq1: str, seq2: str, aln_length: int) -> float:
-            """
-            local hamming distance - defined as total number of matches excluding terminal gaps
-            """
-            # Trim gaps from both sides
-            i, j = 0, aln_length - 1
-            while i < aln_length and (seq1[i] == '-' or seq2[i] == '-'):
-                i += 1
-            while j >= 0 and (seq1[j] == '-' or seq2[j] == '-'):
-                j -= 1
-            if i > j:
-                return 0.0
-
-            seq1_, seq2_ = seq1[i:j + 1], seq2[i:j + 1]
-            matches = sum(c1 == c2 for c1, c2 in zip(seq1_, seq2_))
-            length = j - i + 1
-
-            return (matches / length) * 100 if length > 0 else 0.0
-
-        def ged(seq1: str, seq2: str, aln_length: int = None) -> float:
-            """
-            gap excluded distance - defined as percentage of total number of matches excluding all gaps
-            """
-
-            diff, total = 0, 0
-
-            for c1, c2 in zip(seq1, seq2):
-                if c1 != '-' and c2 != '-':
-                    total += 1
-                    if c1 != c2:
-                        diff += 1
-
-            return (1 - diff / total) * 100 if total > 0 else 0
-
-        def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
-            """
-            gap compressed distance - defined as percentage of total number of matches with sequential gap mismatches
-            counting as a single mismatch
-            """
-            matches = 0
-            mismatches = 0
-            in_gap = False
-
-            for char1, char2 in zip(seq1, seq2):
-                if char1 == '-' and char2 == '-':  # Shared gap: do nothing
-                    continue
-                elif char1 == '-' or char2 == '-':  # Gap in only one sequence
-                    if not in_gap:  # Start of a new gap stretch
-                        mismatches += 1
-                        in_gap = True
-                else:  # No gaps
-                    in_gap = False
-                    if char1 == char2:  # Matching characters
-                        matches += 1
-                    else:  # Mismatched characters
-                        mismatches += 1
-
-            return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
-
-        def jc69(seq1: str, seq2: str, aln_length: int = None) -> float:
-            """
-            Jukes-Cantor 1969 (JC69) corrected identity.
-            Gaps are excluded. The proportion of differing sites (p-distance) is corrected
-            for multiple hits: d = -(3/4) * ln(1 - (4/3) * p).
-            Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
-            Returns 0 when p >= 0.75 (formula undefined / sequence saturated).
-            """
-            diff, total = 0, 0
-            for c1, c2 in zip(seq1, seq2):
-                if c1 != '-' and c2 != '-':
-                    total += 1
-                    if c1 != c2:
-                        diff += 1
-            if total == 0:
-                return 0.0
-            p = diff / total
-            if p == 0.0:
-                return 100.0
-            correction = 1.0 - (4.0 / 3.0) * p
-            if correction <= 0.0:  # saturated – formula undefined
-                return 0.0
-            d = -(3.0 / 4.0) * math.log(correction)
-            return max(0.0, (1.0 - d) * 100.0)
-
-        def k2p(seq1: str, seq2: str, aln_length: int = None) -> float:
-            """
-            Kimura 2-Parameter (K2P / K80) corrected identity.
-            Gaps are excluded. Transitions and transversions (Tv) are weighted separately:
-            d = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q)
-            where P = Ti / total and Q = Tv / total.
-            Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
-            Returns 0 when the logarithm arguments become non-positive (saturated).
-            """
-            transitions = [{'A', 'G'}, {'C', 'T'}]
-
-            ts, tv, total = 0, 0, 0
-            for c1, c2 in zip(seq1, seq2):
-                if c1 != '-' and c2 != '-':
-                    total += 1
-                    if c1 != c2:
-                        if {c1, c2} in transitions:
-                            ts += 1
-                        else:
-                            tv += 1
-            if total == 0:
-                return 0.0
-            if ts == 0 and tv == 0:
-                return 100.0
-            P = ts / total   # transition proportion
-            Q = tv / total   # transversion proportion
-            term1 = 1.0 - 2.0 * P - Q
-            term2 = 1.0 - 2.0 * Q
-            # saturated – formula undefined
-            if term1 <= 0.0 or term2 <= 0.0:
-                return 0.0
-            # calculate distance
-            d = -0.5 * math.log(term1) - 0.25 * math.log(term2)
-
-            return max(0.0, (1 - d) * 100.0)
-
-        # Map distance type to corresponding function
-        distance_functions: Dict[str, Callable[[str, str, int], float]] = {
-            'ghd': ghd,
-            'lhd': lhd,
-            'ged': ged,
-            'gcd': gcd,
-            'jc69': jc69,
-            'k2p': k2p,
-        }
-
-        return distance_functions
-
     # Properties with setters
     @property
     def reference_id(self) -> str:
@@ -1061,7 +886,7 @@ def calc_position_matrix(self, matrix_type:str='PWM') -> ndarray | ValueError:
         """
         Calculates a position matrix of the specified type for the given alignment. The function
         supports generating matrices of types Position Frequency Matrix (PFM), Position Probability
-        Matrix (PPM), Position Weight Matrix (PWM), and cummulative Information Content (IC). It validates
+        Matrix (PPM), Position Weight Matrix (PWM), and cumulative Information Content (IC). It validates
         the provided matrix type and includes pseudo-count adjustments to ensure robust calculations.
 
         :param matrix_type: Type of position matrix to calculate. Accepted values are 'PFM', 'PPM',
@@ -1102,8 +927,6 @@ def calc_position_matrix(self, matrix_type:str='PWM') -> ndarray | ValueError:
         if matrix_type == 'IC':
             return ic
 
-        return None
-
     def calc_percent_recovery(self) -> dict:
         """
         Recovery per sequence either compared to the majority consensus seq
@@ -1204,7 +1027,7 @@ def calc_character_frequencies(self) -> dict:
 
     def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
         """
-        Calculate pairwise identities for an alignment. As there are different definitions of sequence identity, there are different options implemented:
+        Calculate pairwise identities for an alignment. Different options are implemented:
 
         **1) ghd (global hamming distance)**: At each alignment position, check if characters match:
         \ndistance = matches / alignment_length * 100
@@ -1213,24 +1036,26 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
         \ndistance = matches / min(5'3' ungapped seq1, 5'3' ungapped seq2) * 100
 
         **3) ged (gap excluded distance)**: All gaps are excluded from the alignment
-        \ndistance = matches / (matches + mismatches) * 100
+        \ndistance = (1 - mismatches / total) * 100
 
         **4) gcd (gap compressed distance)**: All consecutive gaps are compressed to one mismatch.
         \ndistance = matches / gap_compressed_alignment_length * 100
 
+        RNA/DNA only:
+
         **5) jc69 (Jukes-Cantor 1969)**: Gaps excluded. Applies the JC69 substitution model to correct
         the p-distance for multiple hits (assumes equal base frequencies and substitution rates).
         \ncorrected_identity = (1 - d_JC69) * 100,  where d = -(3/4) * ln(1 - (4/3) * p)
 
-        **6) k2p (Kimura 2-Parameter / K80)**: Gaps excluded. Distinguishes transitions (Ti) and
+        **6) k2p (Kimura 2-Parameter / K80)**: Gaps excluded. Distinguishes transitions (Ts) and
         transversions (Tv). Returns (1 - d_K2P) * 100 as corrected percent identity.
-        \nd = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q),  P = Ti/total, Q = Tv/total
+        \nd = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q),  P = Ts/total, Q = Tv/total
 
-        :param distance_type: type of distance computation technique
+        :param distance_type: type of distance computation: ghd, lhd, ged, gcd and nucleotide only: jc69 and k2p
         :return: array with pairwise distances.
         """
 
-        distance_functions = self._create_distance_calculation_function_mapping()
+        distance_functions = _create_distance_calculation_function_mapping()
 
         if distance_type not in distance_functions:
             raise ValueError(f"Invalid distance type '{distance_type}'. Choose from {list(distance_functions.keys())}.")
@@ -1258,7 +1083,7 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
 
     def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> PairwiseDistanceResult:
         """
-        Calculate pairwise identities between reference and all sequences in the alignment. Same computation as calc_pairwise_identity_matrix but compared to a single sequence. Supported distance computation methods.
+        Calculate pairwise identities between reference and all sequences in the alignment. Same computation as calc_pairwise_identity_matrix but compared to a single sequence. Different options are implemented:
 
         **1) ghd (global hamming distance)**: At each alignment position, check if characters match:
         \ndistance = matches / alignment_length * 100
@@ -1267,22 +1092,27 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
         \ndistance = matches / min(5'3' ungapped seq1, 5'3' ungapped seq2) * 100
 
         **3) ged (gap excluded distance)**: All gaps are excluded from the alignment
-        \ndistance = matches / (matches + mismatches) * 100
+        \ndistance = (1 - mismatches / total) * 100
 
         **4) gcd (gap compressed distance)**: All consecutive gaps are compressed to one mismatch.
         \ndistance = matches / gap_compressed_alignment_length * 100
 
-        **5) jc69 (Jukes-Cantor 1969)**: Gaps excluded. JC69 substitution-model corrected identity.
-        \ncorrected_identity = (1 - d_JC69) * 100
+        RNA/DNA only:
 
-        **6) k2p (Kimura 2-Parameter / K80)**: Gaps excluded. Distinguishes transitions and transversions.
-        \ncorrected_identity = (1 - d_K2P) * 100
+        **5) jc69 (Jukes-Cantor 1969)**: Gaps excluded. Applies the JC69 substitution model to correct
+        the p-distance for multiple hits (assumes equal base frequencies and substitution rates).
+        \ncorrected_identity = (1 - d_JC69) * 100,  where d = -(3/4) * ln(1 - (4/3) * p)
 
-        :param distance_type: type of distance computation technique
+        **6) k2p (Kimura 2-Parameter / K80)**: Gaps excluded. Distinguishes transitions (Ts) and
+        transversions (Tv). Returns (1 - d_K2P) * 100 as corrected percent identity.
+        \nd = -(1/2) * ln(1 - 2P - Q) - (1/4) * ln(1 - 2Q),  P = Ts/total, Q = Tv/total
+
+        :param distance_type: type of distance computation: ghd, lhd, ged, gcd and nucleotide only: jc69 and k2p
+        :return: array with pairwise distances.
         :return: dataclass with reference label, sequence ids and pairwise distances.
         """
 
-        distance_functions = self._create_distance_calculation_function_mapping()
+        distance_functions = _create_distance_calculation_function_mapping()
 
         if distance_type not in distance_functions:
             raise ValueError(f"Invalid distance type '{distance_type}'. Choose from {list(distance_functions.keys())}.")
@@ -1305,7 +1135,7 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
             distance_names.append(seq_id)
             distances.append(distance_func(ref_seq, aln[seq_id], self.length))
 
-        return self.PairwiseDistanceResult(
+        return PairwiseDistanceResult(
             reference_id=ref_id if ref_id is not None else 'consensus',
             sequence_ids=distance_names,
             distances=np.array(distances)
diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index 5b2051f..7234737 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -6,20 +6,9 @@
 ## Functions:
 """
 
-import os
 from numpy import ndarray
 from msaexplorer import config
-
-
-def _check_and_create_path(path: str):
-    """
-    Check and create path if it doesn't exist.
-    :param path: string to file
-    """
-    if path is not None:
-        output_dir = os.path.dirname(path)
-        if output_dir and not os.path.exists(output_dir):
-            os.makedirs(output_dir)
+from msaexplorer._helpers import _check_and_create_path
 
 
 def snps(snp_dict: dict, format_type: str = 'vcf', path: str | None = None) -> str | None | ValueError:
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index 8ffd603..442188c 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -1,9 +1,10 @@
 """Tests for alignment statistics calculation methods in ``MSA``."""
 
 import pytest
+import numpy as np
 from conftest import create_alignment
 from msaexplorer.explore import MSA
-import numpy as np
+from msaexplorer._msa_data_classes import PairwiseDistanceResult
 
 
 class TestCalcEntropy:
@@ -311,7 +312,7 @@ def test_returns_dataclass_for_all_reference_positions(self, sequences):
 
         result = msa.calc_pairwise_distance_to_reference(distance_type="ghd")
 
-        assert isinstance(result, MSA.PairwiseDistanceResult)
+        assert isinstance(result, PairwiseDistanceResult)
         assert result.reference_id == "ref"
         assert result.sequence_ids == ["q1", "q2"]
         assert np.allclose(result.distances, np.array([75.0, 50.0]))

From ab3cce3f1991f4ccddcebec3dfb66714cea18e65 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Tue, 17 Mar 2026 15:25:36 +0100
Subject: [PATCH 06/25] finalized integration into app

---
 app_src/shiny_plots.py                        | 19 +++++++++-
 app_src/shiny_server.py                       | 31 +++++++++++-----
 ...{_msa_data_classes.py => _data_classes.py} |  0
 msaexplorer/_helpers.py                       | 11 +++---
 msaexplorer/draw.py                           |  2 +-
 msaexplorer/explore.py                        | 36 +++++++++----------
 tests/test_stats_calculations.py              |  2 +-
 7 files changed, 66 insertions(+), 35 deletions(-)
 rename msaexplorer/{_msa_data_classes.py => _data_classes.py} (100%)

diff --git a/app_src/shiny_plots.py b/app_src/shiny_plots.py
index 3f5db7a..f67d4a0 100644
--- a/app_src/shiny_plots.py
+++ b/app_src/shiny_plots.py
@@ -44,7 +44,7 @@ def create_msa_plot(aln, ann, inputs, fig_size=None) -> plt.Figure | None:
     plot_functions = []
 
     # First plot
-    if inputs['stat_type'] not in ['Off', 'sequence logo']:
+    if inputs['stat_type'] not in ['Off', 'sequence logo', 'simplot-like similarity']:
         height_ratios.append(inputs['plot_1_size'])
         plot_functions.append(
             lambda ax: draw.stat_plot(
@@ -66,6 +66,23 @@ def create_msa_plot(aln, ann, inputs, fig_size=None) -> plt.Figure | None:
                 color_scheme = inputs['logo_coloring']
             )
         )
+    elif inputs['stat_type'] == 'simplot-like similarity':
+        height_ratios.append(inputs['plot_1_size'])
+        # define window size and step size automatically
+        if aln.length > 100:
+            window_size = 50
+        else:
+            window_size = 2
+        window_step = int(max([window_size/10, 1]))
+        plot_functions.append(
+            lambda ax: draw.simplot(
+                aln, ax=ax, ref=aln.reference_id,
+                window_size=window_size,
+                step_size=window_step,
+                distance_calculation='ged' if aln.aln_type == 'AA' else 'k2p',
+                colors=None
+            )
+        )
 
     # Second plot
     if inputs['alignment_type'] != 'Off':
diff --git a/app_src/shiny_server.py b/app_src/shiny_server.py
index 71ae5b7..0633a49 100644
--- a/app_src/shiny_server.py
+++ b/app_src/shiny_server.py
@@ -333,7 +333,7 @@ def finalize_loaded_alignment(aln, annotation_file):
 
         if aln.aln_type == 'AA':
             ui.update_selectize('stat_type',
-                                choices=['Off', 'sequence logo', 'entropy', 'coverage', 'identity', 'similarity'],
+                                choices=['Off', 'sequence logo', 'entropy', 'coverage', 'identity', 'similarity', 'simplot-like similarity'],
                                 selected='Off')
             ui.update_selectize('download_type',
                                 choices=['alignment','SNPs', 'consensus', 'character frequencies', '% recovery', 'entropy',
@@ -345,7 +345,7 @@ def finalize_loaded_alignment(aln, annotation_file):
         else:
             ui.update_selectize('stat_type',
                                 choices=['Off', 'sequence logo', 'gc', 'entropy', 'coverage', 'identity', 'similarity',
-                                         'ts tv score', 'gap frequency'], selected='Off')
+                                         'ts tv score', 'gap frequency', 'simplot-like similarity'], selected='Off')
             ui.update_selectize('download_type', choices=['alignment', 'SNPs', 'consensus', 'character frequencies', '% recovery',
                                                           'reverse complement alignment', 'conserved orfs', 'gc',
                                                           'entropy', 'gap frequency', 'coverage', 'mean identity', 'mean similarity',
@@ -914,22 +914,31 @@ def update_additional_options_left():
         """
         Update UI for the left plot in the analysis tab
         """
+
+        aln = reactive.alignment.get()
+        if aln is None:
+            return None
+
         # ensure that it is switched back
         if input.analysis_plot_type_left() == 'Off':
             ui.remove_ui(selector="div:has(> #additional_analysis_options_left)")
             ui.remove_ui(selector="div:has(> #additional_analysis_options_left-label)")
             ui.remove_ui(selector="div:has(> #analysis_info_left)")
         if input.analysis_plot_type_left() == 'Pairwise identity':
+            choices = {
+                'ghd': 'global hamming distance',
+                'lhd': 'local hamming distance',
+                'ged': 'gap excluded distance',
+                'gcd': 'gap compressed distance'
+            }
+            if aln.aln_type != 'AA':
+                choices['jc69'] = 'Jukes-Cantor 1969 distance'
+                choices['k2p'] = 'Kimura 2-Parameter (K2P / K80) distance'
             ui.insert_ui(
                 ui.input_selectize(
                     'additional_analysis_options_left',
                     label='Options left',
-                    choices={
-                        'ghd': 'global hamming distance',
-                        'lhd': 'local hamming distance',
-                        'ged': 'gap excluded distance',
-                        'gcd': 'gap compressed distance'
-                    },
+                    choices=choices,
                     selected='ghd'
                 ),
                 selector='#analysis_plot_type_right-label',
@@ -958,7 +967,11 @@ def analysis_info_left():
         elif selected_option == 'ged':
             return 'INFO ged (gap excluded distance):\n\nAll gaps are excluded from the \nalignment\n\ndistance = matches / (matches + mismatches) * 100'
         elif selected_option == 'gcd':
-            return 'INFO gcd (gap compressed distance):\n\nAll consecutive gaps arecompressed to\none mismatch.\n\ndistance = matches / gap_compressed_alignment_length * 100'
+            return 'INFO gcd (gap compressed distance):\n\nAll consecutive gaps are compressed to\none mismatch.\n\ndistance = matches / gap_compressed_alignment_length * 100'
+        elif selected_option == 'jc69':
+            return 'INFO jc69 (Jukes-Cantor 1969 distance):\n\nDistance calculation based ungapped sequences assuming equal substitution rates.'
+        elif selected_option == 'k2p':
+            return 'INFO k2p (Kimura 2-Parameter (K2P / K80) distance):\n\nDistance calculation based ungapped sequences assuming unequal substitution rates\n depending on transitions or transversions.'
         else:
             return None
 
diff --git a/msaexplorer/_msa_data_classes.py b/msaexplorer/_data_classes.py
similarity index 100%
rename from msaexplorer/_msa_data_classes.py
rename to msaexplorer/_data_classes.py
diff --git a/msaexplorer/_helpers.py b/msaexplorer/_helpers.py
index 4473c78..acc66d6 100644
--- a/msaexplorer/_helpers.py
+++ b/msaexplorer/_helpers.py
@@ -3,6 +3,7 @@
 """
 
 import os, io, math
+import numpy as np
 from typing import Callable, Dict
 from matplotlib.colors import is_color_like
 
@@ -46,7 +47,7 @@ def lhd(seq1: str, seq2: str, aln_length: int) -> float:
         matches = sum(c1 == c2 for c1, c2 in zip(seq1_, seq2_))
         length = j - i + 1
 
-        return (matches / length) * 100 if length > 0 else 0.0
+        return (matches / length) * 100 if length > 0 else np.nan
 
     def ged(seq1: str, seq2: str, aln_length: int = None) -> float:
         """
@@ -61,7 +62,7 @@ def ged(seq1: str, seq2: str, aln_length: int = None) -> float:
                 if c1 != c2:
                     diff += 1
 
-        return (1 - diff / total) * 100 if total > 0 else 0
+        return (1 - diff / total) * 100 if total > 0 else np.nan
 
     def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
         """
@@ -86,7 +87,7 @@ def gcd(seq1: str, seq2: str, aln_length: int = None) -> float:
                 else:  # Mismatched characters
                     mismatches += 1
 
-        return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else 0
+        return matches / (matches + mismatches) * 100 if (matches + mismatches) > 0 else np.nan
 
     def jc69(seq1: str, seq2: str, aln_length: int = None) -> float:
         """
@@ -103,7 +104,7 @@ def jc69(seq1: str, seq2: str, aln_length: int = None) -> float:
                 if c1 != c2:
                     diff += 1
         if total == 0:
-            return 0.0
+            return np.nan
         p = diff / total
         if p == 0.0:
             return 100.0
@@ -134,7 +135,7 @@ def k2p(seq1: str, seq2: str, aln_length: int = None) -> float:
                     else:
                         tv += 1
         if total == 0:
-            return 0.0
+            return np.nan
         if ts == 0 and tv == 0:
             return 100.0
         P = ts / total   # transition proportion
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 58e9171..d4a15da 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -1181,7 +1181,7 @@ def consensus_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, threshold
     return ax
 
 def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None, colors: str | list | None = None,
-            window_size: int = 200, step_size: int = 20, distance_calculation: str = 'ghd', line_width: int | float = 0.5,
+            window_size: int = 200, step_size: int = 20, distance_calculation: str = 'ged', line_width: int | float = 1,
             show_legend: bool = False, show_reference: bool = True, bbox_to_anchor: tuple[float|int, float|int] | list= (1, 1),
             show_x_label: bool = False) -> plt.Axes:
     """
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 580143a..f42935a 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -24,10 +24,10 @@
 
 # msaexplorer
 from msaexplorer import config
-from msaexplorer._msa_data_classes import PairwiseDistanceResult
+from msaexplorer._data_classes import PairwiseDistanceResult
 from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping
 
-
+#TODO: Move outputs to dataclasses
 class MSA:
     """
     An alignment class that allows computation of several stats. Supported inputs are file paths to alignments in "fasta",
@@ -882,7 +882,7 @@ def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=Tr
 
         return similarity_array
 
-    def calc_position_matrix(self, matrix_type:str='PWM') -> ndarray | ValueError:
+    def calc_position_matrix(self, matrix_type:str='PWM') -> None | ndarray | ValueError:
         """
         Calculates a position matrix of the specified type for the given alignment. The function
         supports generating matrices of types Position Frequency Matrix (PFM), Position Probability
@@ -1223,10 +1223,10 @@ def calc_transition_transversion_score(self) -> list:
 
 class Annotation:
     """
-    An annotation class that allows to read in gff, gb or bed files and adjust its locations to that of the MSA.
+    An annotation class that allows to read in gff, gb, or bed files and adjust its locations to that of the MSA.
     """
 
-    def __init__(self, aln: MSA, annotation_path: str):
+    def __init__(self, aln: MSA, annotation: str | GenBankIterator):
         """
         The annotation class. Lets you parse multiple standard formats
         which might be used for annotating an alignment. The main purpose
@@ -1236,11 +1236,11 @@ def __init__(self, aln: MSA, annotation_path: str):
         and the MSA have to partly match.
 
         :param aln: MSA class
-        :param annotation_path: path to annotation file (gb, bed, gff) or raw string
+        :param annotation: path to file (gb, bed, gff) or raw string or GenBankIterator from biopython
 
         """
 
-        self.ann_type, self._seq_id, self.locus, self.features  = self._parse_annotation(annotation_path, aln)  # read annotation
+        self.ann_type, self._seq_id, self.locus, self.features  = self._parse_annotation(annotation, aln)  # read annotation
         self._gapped_seq = self._MSA_validation_and_seq_extraction(aln, self._seq_id)  # extract gapped sequence
         self._position_map = self._build_position_map()  # build a position map
         self._map_to_alignment()  # adapt feature locations
@@ -1259,15 +1259,15 @@ def _MSA_validation_and_seq_extraction(aln: MSA, seq_id: str) -> str:
             return aln._alignment[seq_id]
 
     @staticmethod
-    def _parse_annotation(annotation_path: str, aln: MSA) -> tuple[str, str, str, Dict]:
+    def _parse_annotation(annotation: str | GenBankIterator, aln: MSA) -> tuple[str, str, str, Dict]:
 
-        def detect_annotation_type(handle) -> str:
+        def detect_annotation_type(handle: str | GenBankIterator) -> str:
             """
             Detect the type of annotation file (GenBank, GFF, or BED) based
             on the first relevant line (excluding empty and #). Also recognizes
             Bio.SeqIO iterators as GenBank format.
 
-            :param file_path: Path to the annotation file or Bio.SeqIO iterator for genbank records read with biopython.
+            :param handle: Path to the annotation file or Bio.SeqIO iterator for genbank records read with biopython.
             :return: The detected file type ('gb', 'gff', or 'bed').
 
             :raises ValueError: If the file type cannot be determined.
@@ -1300,22 +1300,22 @@ def detect_annotation_type(handle) -> str:
             raise ValueError(
                 "File type could not be determined. Ensure the file follows a recognized format (GenBank, GFF, or BED).")
 
-        def parse_gb(file_path) -> dict:
+        def parse_gb(file: str | GenBankIterator) -> dict:
             """
             Parse a GenBank file into the same dictionary structure used by the annotation pipeline.
 
-            :param file_path: path to genbank file, raw string, or Bio.SeqIO iterator
+            :param file: path to genbank file, raw string, or Bio.SeqIO iterator
             :return: nested dictionary
             """
             records = {}
 
             # Check if input is a GenBankIterator
-            if isinstance(file_path, GenBankIterator):
+            if isinstance(file, GenBankIterator):
                 # Direct GenBankIterator input
-                seq_records = list(file_path)
+                seq_records = list(file)
             else:
                 # File path or string input
-                with _get_line_iterator(file_path) as handle:
+                with _get_line_iterator(file) as handle:
                     seq_records = list(SeqIO.parse(handle, "genbank"))
 
             for seq_record in seq_records:
@@ -1459,12 +1459,12 @@ def parse_bed(file_path) -> dict:
         }
         # determine the annotation content -> should be standard formatted
         try:
-            annotation_type = detect_annotation_type(annotation_path)
+            annotation_type = detect_annotation_type(annotation)
         except ValueError as err:
             raise err
 
         # read in the annotation
-        annotations = parse_functions[annotation_type](annotation_path)
+        annotations = parse_functions[annotation_type](annotation)
 
         # sanity check whether one of the annotation ids and alignment ids match
         annotation_found = False
@@ -1480,7 +1480,7 @@ def parse_bed(file_path) -> dict:
                     break
 
         if not annotation_found:
-            raise ValueError(f'the annotations of {annotation_path} do not match any ids in the MSA')
+            raise ValueError(f'the annotations of {annotation} do not match any ids in the MSA')
 
         # return only the annotation that has been found, the respective type and the seq_id to map to
         return annotation_type, aln_id, annotations[annotation]['locus'], annotations[annotation]['features']
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index 442188c..eebea2e 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -4,7 +4,7 @@
 import numpy as np
 from conftest import create_alignment
 from msaexplorer.explore import MSA
-from msaexplorer._msa_data_classes import PairwiseDistanceResult
+from msaexplorer._data_classes import PairwiseDistanceResult
 
 
 class TestCalcEntropy:

From 45bce680133710cfb862d1bdf24a82b073af9345 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Wed, 18 Mar 2026 11:45:10 +0100
Subject: [PATCH 07/25] implemented data class for stats

---
 app_src/shiny_plots.py           |  2 +-
 app_src/shiny_server.py          |  9 ++--
 msaexplorer/_data_classes.py     | 69 +++++++++++++++++++++++++++----
 msaexplorer/draw.py              |  9 ++--
 msaexplorer/explore.py           | 52 +++++++++++++++++-------
 msaexplorer/export.py            | 17 ++++++--
 tests/test_stats_calculations.py | 70 ++++++++++++++++----------------
 7 files changed, 161 insertions(+), 67 deletions(-)

diff --git a/app_src/shiny_plots.py b/app_src/shiny_plots.py
index f67d4a0..585df33 100644
--- a/app_src/shiny_plots.py
+++ b/app_src/shiny_plots.py
@@ -199,7 +199,7 @@ def create_analysis_custom_heatmap(aln, inputs):
     else:
         figure_size = int(inputs['dimensions']['width'] * 0.7)
 
-    matrix = aln.calc_pairwise_identity_matrix(inputs['additional_analysis_options_left'])
+    matrix = aln.calc_pairwise_identity_matrix(inputs['additional_analysis_options_left']).distances
     labels = [x.split(' ')[0] for x in list(aln.alignment.keys())]
 
     # generate hover text
diff --git a/app_src/shiny_server.py b/app_src/shiny_server.py
index 0633a49..3dc644d 100644
--- a/app_src/shiny_server.py
+++ b/app_src/shiny_server.py
@@ -34,6 +34,7 @@
 
 # msaexplorer
 from msaexplorer import explore, config, export, draw
+from msaexplorer._data_classes import AlignmentStats
 
 
 def server(input, output, session):
@@ -715,7 +716,7 @@ def _consensus_option():
 
         def _stat_option():
             # create function mapping
-            stat_functions: Dict[str, Callable[[], list | ndarray]] = {
+            stat_functions: Dict[str, Callable[[], AlignmentStats | ndarray]] = {
                 'gc': aln.calc_gc,
                 'entropy': aln.calc_entropy,
                 'coverage': aln.calc_coverage,
@@ -735,8 +736,10 @@ def _stat_option():
                     break
             # use correct function
             data = stat_functions[stat_type]()
-            # calculate the mean (identical to draw module of msaexplorer)
-            if stat_type in ['mean identity', 'mean similarity']:
+            if isinstance(data, AlignmentStats):
+                data = data.values
+            # calculate the mean for identity or similarity (identical to draw module of msaexplorer)
+            else:
                 # for the mean nan values get handled as the lowest possible number in the matrix
                 data = np.nan_to_num(data, True, -1 if stat_type == 'identity' else 0)
                 data = np.mean(data, axis=0)
diff --git a/msaexplorer/_data_classes.py b/msaexplorer/_data_classes.py
index a91563d..b64881c 100644
--- a/msaexplorer/_data_classes.py
+++ b/msaexplorer/_data_classes.py
@@ -8,21 +8,76 @@
 # libs
 from numpy import ndarray
 
+
 @dataclass(frozen=True)
-class PairwiseDistanceResult:
+class AlignmentStats:
+    """
+    Generic result container for position-based statistics.
     """
-    Result container for pairwise identity values between a reference/consensus
-    sequence and each sequence in the alignment.
 
-    The object remains iterable so it can be unpacked as a tuple:
-    ``reference_label, sequence_ids, distances = result``.
+    stat_name: str
+    positions: ndarray
+    values: ndarray
+    aln_type: str
+    reference_id: str | None
+
+    def __post_init__(self):
+        if self.positions.shape != self.values.shape:
+            raise ValueError("positions and values must have the same shape")
+
+    # dunder methods
+    def __len__(self) -> int:
+        return len(self.values)
+
+    def __getitem__(self, index: int) -> float:
+        return self.values[index]
+
+    def __iter__(self):
+        yield self.stat_name
+        yield self.positions
+        yield self.values
+
+    def __contains__(self, item: float) -> bool:
+        return item in self.positions
+
+    # normal mehods
+    def as_array(self) -> ndarray:
+        return self.values
+
+    def as_list(self) -> list:
+        return self.values.tolist()
+
+
+@dataclass(frozen=True)
+class PairwiseDistance:
+    """
+    Result container for Pairwise distances. Array can either be a 2D (compared to reference) or 3D array.
     """
 
-    reference_id: str
+    reference_id: str | None
     sequence_ids: list[str]
     distances: ndarray
 
+    # dunder methods
     def __iter__(self):
         yield self.reference_id
         yield self.sequence_ids
-        yield self.distances
\ No newline at end of file
+        yield self.distances
+
+    def __len__(self) -> int:
+        return len(self.sequence_ids)
+
+    def __getitem__(self, index: int | str) -> float:
+        """
+        Different ways to access the distance matrix.
+        - pd[0] -> Distance first sequence to all other sequences
+        - pd['seq_name'] -> Distance of a specific sequence to all other sequences
+        """
+        if isinstance(index, str):
+            idx = self.sequence_ids.index(index)
+            return self.distances[idx]
+        return self.distances[index]
+
+    def __contains__(self, item: str) -> bool:
+        """Seq ID present"""
+        return item in self.sequence_ids
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index d4a15da..80d3741 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -18,6 +18,7 @@
 
 # MSAexplorer
 from msaexplorer import explore, config
+from msaexplorer._data_classes import AlignmentStats
 from msaexplorer._helpers import _validate_color, _get_contrast_text_color
 
 # libs
@@ -675,7 +676,7 @@ def stat_plot(aln: explore.MSA | str, stat_type: str, ax: plt.Axes | None = None
     aln, ax = _validate_input_parameters(aln, ax)
 
     # define possible functions to calc here
-    stat_functions: Dict[str, Callable[[], list | ndarray]] = {
+    stat_functions: Dict[str, Callable[[], AlignmentStats | ndarray]] = {
         'gc': aln.calc_gc,
         'entropy': aln.calc_entropy,
         'coverage': aln.calc_coverage,
@@ -694,7 +695,9 @@ def stat_plot(aln: explore.MSA | str, stat_type: str, ax: plt.Axes | None = None
 
     # generate input data
     array = stat_functions[stat_type]()
-
+    if isinstance(array, AlignmentStats):
+        array = array.values
+    # define possible spans for values
     if stat_type == 'identity':
         min_value, max_value = -1, 0
     elif stat_type == 'ts tv score':
@@ -1193,7 +1196,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
 
     :param aln: alignment MSA class or path
     :param ax: matplotlib axes
-    :param ref: reference sequence id or None. For None all computations are compared to a majority consensus
+    :param ref: reference sequence id or None. For 'None' all computations are compared to a majority consensus
     :param colors: color for each sequence. can be a single named color or a list of named colors or a plt.colormap or None (auto coloring)
     :param window_size: window size for sliding window
     :param step_size: step size for sliding window
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index f42935a..8a4bf9d 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -24,7 +24,7 @@
 
 # msaexplorer
 from msaexplorer import config
-from msaexplorer._data_classes import PairwiseDistanceResult
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats
 from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping
 
 #TODO: Move outputs to dataclasses
@@ -215,6 +215,24 @@ def alignment(self) -> dict:
             return self._alignment
 
     # functions for different alignment stats
+    def _create_position_stat_result(self, stat_name: str, values: list | ndarray) -> AlignmentStats:
+        """
+        Build a shared dataclass for position-wise statistics.
+        """
+        values_array = np.asarray(values, dtype=float)
+        if self.zoom is None:
+            positions = np.arange(self.length, dtype=int)
+        else:
+            positions = np.arange(self.zoom[0], self.zoom[0] + self.length, dtype=int)
+
+        return AlignmentStats(
+            stat_name=stat_name,
+            positions=positions,
+            values=values_array,
+            aln_type=self.aln_type,
+            reference_id=self.reference_id,
+        )
+
     def get_reference_coords(self) -> tuple[int, int]:
         """
         Determine the start and end coordinates of the reference sequence
@@ -575,7 +593,7 @@ def calc_length_stats(self) -> dict:
                 'max length': int(np.max(seq_lengths))
                 }
 
-    def calc_entropy(self) -> list:
+    def calc_entropy(self) -> AlignmentStats:
         """
         Calculate the normalized shannon's entropy for every position in an alignment:
 
@@ -642,9 +660,9 @@ def shannons_entropy(character_list: list, states: int, aln_type: str) -> float:
                 pos.append(aln[record][nuc_pos])
             entropys.append(shannons_entropy(pos, states, self.aln_type))
 
-        return entropys
+        return self._create_position_stat_result('entropy', entropys)
 
-    def calc_gc(self) -> list | TypeError:
+    def calc_gc(self) -> AlignmentStats | TypeError:
         """
         Determine the GC content for every position in an nt alignment.
         :return: GC content for every position.
@@ -674,9 +692,9 @@ def calc_gc(self) -> list | TypeError:
                 sum([nucleotides.count(x) * to_count[x] for x in to_count]) / len(nucleotides)
             )
 
-        return gc
+        return self._create_position_stat_result('gc', gc)
 
-    def calc_coverage(self) -> list:
+    def calc_coverage(self) -> AlignmentStats:
         """
         Determine the coverage of every position in an alignment.
         This is defined as:
@@ -692,15 +710,15 @@ def calc_coverage(self) -> list:
                 pos = pos + aln[record][nuc_pos]
             coverage.append(1 - pos.count('-') / len(pos))
 
-        return coverage
+        return self._create_position_stat_result('coverage', coverage)
 
-    def calc_gap_frequency(self) -> list:
+    def calc_gap_frequency(self) -> AlignmentStats:
         """
         Determine the gap frequency for every position in an alignment. This is the inverted coverage.
         """
         coverage = self.calc_coverage()
 
-        return [1 - x for x in coverage]
+        return self._create_position_stat_result('gap frequency', 1 - coverage.values)
 
     def calc_reverse_complement_alignment(self) -> dict | TypeError:
         """
@@ -1025,7 +1043,7 @@ def calc_character_frequencies(self) -> dict:
 
         return freqs
 
-    def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
+    def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> PairwiseDistance:
         """
         Calculate pairwise identities for an alignment. Different options are implemented:
 
@@ -1079,9 +1097,13 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> ndarray:
                 distance_matrix[i, j] = dist
                 distance_matrix[j, i] = dist
 
-        return distance_matrix
+        return PairwiseDistance(
+            reference_id=None,
+            sequence_ids=list(aln.keys()),
+            distances=distance_matrix
+        )
 
-    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> PairwiseDistanceResult:
+    def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> PairwiseDistance:
         """
         Calculate pairwise identities between reference and all sequences in the alignment. Same computation as calc_pairwise_identity_matrix but compared to a single sequence. Different options are implemented:
 
@@ -1135,7 +1157,7 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
             distance_names.append(seq_id)
             distances.append(distance_func(ref_seq, aln[seq_id], self.length))
 
-        return PairwiseDistanceResult(
+        return PairwiseDistance(
             reference_id=ref_id if ref_id is not None else 'consensus',
             sequence_ids=distance_names,
             distances=np.array(distances)
@@ -1194,7 +1216,7 @@ def get_snps(self, include_ambig:bool=False) -> dict:
 
         return snp_dict
 
-    def calc_transition_transversion_score(self) -> list:
+    def calc_transition_transversion_score(self) -> AlignmentStats:
         """
         Based on the snp positions, calculates a transition/transversions score.
         A positive score means higher ratio of transitions and negative score means
@@ -1218,7 +1240,7 @@ def calc_transition_transversion_score(self) -> list:
                 else:
                     score[pos] -= snps['POS'][pos]['ALT'][alt]['AF']
 
-        return score
+        return self._create_position_stat_result('ts tv score', score)
 
 
 class Annotation:
diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index 7234737..c85282f 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -6,8 +6,10 @@
 ## Functions:
 """
 
+import numpy as np
 from numpy import ndarray
 from msaexplorer import config
+from msaexplorer._data_classes import AlignmentStats
 from msaexplorer._helpers import _check_and_create_path
 
 
@@ -147,11 +149,11 @@ def _validate_sequence(seq: str):
         return fasta_formated_sequence
 
 
-def stats(stat_data: list | ndarray, seperator: str = '\t', path: str | None = None) -> str | None:
+def stats(stat_data: AlignmentStats | list | ndarray, seperator: str = '\t', path: str | None = None) -> str | None:
     """
     Export a list of stats per nucleotide to tabular or csv format.
 
-    :param stat_data: list of stat values
+    :param stat_data: position statistic dataclass or list/array of values
     :param seperator: seperator for values and index
     :param path: path to save the file
     :return: tabular/csv formatted string
@@ -161,8 +163,15 @@ def stats(stat_data: list | ndarray, seperator: str = '\t', path: str | None = N
 
     lines = [f'position{seperator}value']
 
-    for idx, stat_val in enumerate(stat_data):
-        lines.append(f'{idx}{seperator}{stat_val}')
+    if isinstance(stat_data, AlignmentStats):
+        positions = stat_data.positions
+        values = stat_data.values
+    else:
+        values = stat_data
+        positions = np.arange(len(values), dtype=int)
+
+    for position, stat_val in zip(positions, values):
+        lines.append(f'{position}{seperator}{stat_val}')
 
     if path is not None:
         with open(path, 'w') as out_file:
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index eebea2e..beac784 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -4,18 +4,18 @@
 import numpy as np
 from conftest import create_alignment
 from msaexplorer.explore import MSA
-from msaexplorer._data_classes import PairwiseDistanceResult
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats
 
 
 class TestCalcEntropy:
     """Tests for calc_entropy."""
 
-    def test_returns_list_of_correct_length(self):
-        """Entropy list has length equal to alignment length."""
+    def test_returns_dataclass_of_correct_length(self):
+        """Entropy output has length equal to alignment length."""
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT"}))
         entropy = msa.calc_entropy()
 
-        assert isinstance(entropy, list)
+        assert isinstance(entropy, AlignmentStats)
         assert len(entropy) == 8
 
     def test_identical_sequences_is_zero(self):
@@ -23,46 +23,46 @@ def test_identical_sequences_is_zero(self):
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT"}))
         entropy = msa.calc_entropy()
 
-        assert all(e == 0 for e in entropy)
+        assert np.all(entropy.values == 0)
 
     def test_perfectly_mixed_position(self):
         """Position with equal frequencies has maximum entropy."""
         msa = MSA(create_alignment({"s1": "AAAA", "s2": "CCCC", "s3": "GGGG", "s4": "TTTT"}))
         entropy = msa.calc_entropy()
 
-        assert all(e == 1 for e in entropy)
+        assert np.all(entropy.values == 1)
 
     def test_normalized_between_zero_and_one(self):
         """Entropy values are normalized to [0, 1]."""
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT", "s3": "GCGTACGT"}))
         entropy = msa.calc_entropy()
 
-        assert all(0 <= e <= 1 for e in entropy)
+        assert np.all((entropy.values >= 0) & (entropy.values <= 1))
 
     def test_single_gap_position(self):
         """Gaps are handled correctly (ignored in entropy calculation)."""
         msa = MSA(create_alignment({"s1": "A-GT", "s2": "ACGT", "s3": "ACGT"}))
         entropy = msa.calc_entropy()
 
-        assert entropy[1] == 0
+        assert entropy.values[1] == 0
 
     def test_with_ambiguous_nucleotides(self):
         """Ambiguous nucleotides are handled in entropy calculation."""
         msa = MSA(create_alignment({"s1": "ARGT", "s2": "ACGT"}))
         entropy = msa.calc_entropy()
 
-        assert entropy[1] == 0.75
+        assert entropy.values[1] == 0.75
 
 
 class TestCalcGC:
     """Tests for calc_gc."""
 
-    def test_returns_list_of_correct_length(self):
-        """GC list has length equal to alignment length."""
+    def test_returns_dataclass_of_correct_length(self):
+        """GC output has length equal to alignment length."""
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT"}))
         gc = msa.calc_gc()
 
-        assert isinstance(gc, list)
+        assert isinstance(gc, AlignmentStats)
         assert len(gc) == 8
 
     def test_mixed_bases_partial_gc(self):
@@ -70,10 +70,10 @@ def test_mixed_bases_partial_gc(self):
         msa = MSA(create_alignment({"s1": "AAGT", "s2": "ACGT"}))
         gc = msa.calc_gc()
 
-        assert gc[0] == 0.0
-        assert gc[1] == 0.5
-        assert gc[2] == 1.0
-        assert gc[3] == 0.0
+        assert gc.values[0] == 0.0
+        assert gc.values[1] == 0.5
+        assert gc.values[2] == 1.0
+        assert gc.values[3] == 0.0
 
     def test_raises_error_for_aa_alignment(self):
         """GC calculation raises TypeError for amino acid alignments."""
@@ -86,18 +86,18 @@ def test_with_ambiguous_nucleotides(self):
         msa = MSA(create_alignment({"s1": "SWGT", "s2": "ACGT"}))
         gc = msa.calc_gc()
 
-        assert all(g == 0.5 for g in gc[0:1])
+        assert np.all(gc.values[0:1] == 0.5)
 
 
 class TestCalcCoverage:
     """Tests for calc_coverage."""
 
-    def test_returns_list_of_correct_length(self):
-        """Coverage list has length equal to alignment length."""
+    def test_returns_dataclass_of_correct_length(self):
+        """Coverage output has length equal to alignment length."""
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT"}))
         coverage = msa.calc_coverage()
 
-        assert isinstance(coverage, list)
+        assert isinstance(coverage, AlignmentStats)
         assert len(coverage) == 8
 
     def test_with_gaps(self):
@@ -105,20 +105,20 @@ def test_with_gaps(self):
         msa = MSA(create_alignment({"s1": "AT--", "s2": "C---", "s3": "G---", "s4": "GT--"}))
         coverage = msa.calc_coverage()
 
-        assert coverage[0] == 1.0
-        assert coverage[1] == 0.5
-        assert all(c == 0.0 for c in coverage[2:])
+        assert coverage.values[0] == 1.0
+        assert coverage.values[1] == 0.5
+        assert np.all(coverage.values[2:] == 0.0)
 
 
 class TestCalcTransitionTransversionScore:
     """Tests for calc_transition_transversion_score."""
 
-    def test_returns_list_of_correct_length(self):
-        """TS/TV score list has length equal to alignment length."""
+    def test_returns_dataclass_of_correct_length(self):
+        """TS/TV score output has length equal to alignment length."""
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT"}))
         score = msa.calc_transition_transversion_score()
 
-        assert isinstance(score, list)
+        assert isinstance(score, AlignmentStats)
         assert len(score) == 8
 
     def test_identical_sequences_is_zero(self):
@@ -126,23 +126,23 @@ def test_identical_sequences_is_zero(self):
         msa = MSA(create_alignment({"s1": "ACGTACGT", "s2": "ACGTACGT"}))
         score = msa.calc_transition_transversion_score()
 
-        assert all(s == 0 for s in score)
+        assert np.all(score.values == 0)
 
     def test_transition_is_positive(self):
         """Transition substitutions (A<->G, C<->T) are positive."""
         msa = MSA(create_alignment({"s1": "AGAG", "s2": "AAAA"}))
         score = msa.calc_transition_transversion_score()
 
-        assert score[1] == 0.5
-        assert score[3] == 0.5
+        assert score.values[1] == 0.5
+        assert score.values[3] == 0.5
 
     def test_transversion_is_negative(self):
         """Transversion substitutions (A<->C, A<->T, G<->C, G<->T) are negative."""
         msa = MSA(create_alignment({"s1": "ACAC", "s2": "AAAA"}))
         score = msa.calc_transition_transversion_score()
 
-        assert score[1] == -0.5
-        assert score[3] == -0.5
+        assert score.values[1] == -0.5
+        assert score.values[3] == -0.5
 
 
     def test_raises_error_for_aa_alignment(self):
@@ -156,7 +156,7 @@ def test_rna_alignment(self):
         msa = MSA(create_alignment({"s1": "AGAGAG", "s2": "AUAAAA"}))
         score = msa.calc_transition_transversion_score()
 
-        assert score[1] == -0.5
+        assert score.values[1] == -0.5
 
 
 class TestGetSnps:
@@ -282,8 +282,10 @@ class TestCalcPairwiseIdentityMatrix:
     )
     def test_distance_types(self, distance_type, expected_offdiag):
         msa = MSA(create_alignment({"s1": "-A-CGT-", "s2": "TAACG--"}))
-        matrix = msa.calc_pairwise_identity_matrix(distance_type=distance_type)
+        pairwise_distance = msa.calc_pairwise_identity_matrix(distance_type=distance_type)
+        assert isinstance(pairwise_distance, PairwiseDistance)
 
+        matrix = pairwise_distance.distances
         assert matrix.shape == (2, 2)
         assert matrix[0, 0] == 100.0
         assert matrix[1, 1] == 100.0
@@ -312,7 +314,7 @@ def test_returns_dataclass_for_all_reference_positions(self, sequences):
 
         result = msa.calc_pairwise_distance_to_reference(distance_type="ghd")
 
-        assert isinstance(result, PairwiseDistanceResult)
+        assert isinstance(result, PairwiseDistance)
         assert result.reference_id == "ref"
         assert result.sequence_ids == ["q1", "q2"]
         assert np.allclose(result.distances, np.array([75.0, 50.0]))

From 0647056ce43d0b218112352ccbb833615e78163a Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Wed, 18 Mar 2026 12:41:17 +0100
Subject: [PATCH 08/25] implemented smaller helpers to alignment class

---
 msaexplorer/explore.py | 67 ++++++++++++++++++++++++++----------------
 1 file changed, 42 insertions(+), 25 deletions(-)

diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 8a4bf9d..7eece3a 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -35,9 +35,10 @@ class MSA:
     compatibility with Biopython.
     """
 
+    # dunder methods
     def __init__(self, alignment_string: str | MultipleSeqAlignment, reference_id: str = None, zoom_range: tuple | int = None):
         """
-        Initialise an Alignment object.
+        Initialize an Alignment object.
         :param alignment_string: Path to alignment file or raw alignment string
         :param reference_id: reference id
         :param zoom_range: start and stop positions to zoom into the alignment
@@ -47,6 +48,16 @@ def __init__(self, alignment_string: str | MultipleSeqAlignment, reference_id: s
         self._zoom = self._validate_zoom(zoom_range, self._alignment)
         self._aln_type = self._determine_aln_type(self._alignment)
 
+    def __len__(self) -> int:
+        return len(next(iter(self.alignment.values())))
+
+    def __getitem__(self, key: str) -> str:
+        return self.alignment[key]
+
+    def __contains__(self, item: str) -> bool:
+        """Seq ID present"""
+        return item in self.sequence_ids
+
     # Static methods
     @staticmethod
     def _read_alignment(source: str | MultipleSeqAlignment) -> dict:
@@ -196,6 +207,10 @@ def aln_type(self) -> str:
         """
         return self._aln_type
 
+    @property
+    def sequence_ids(self) -> list:
+        return list(self.alignment.keys())
+
     # On the fly properties without setters
     @property
     def length(self) -> int:
@@ -233,6 +248,14 @@ def _create_position_stat_result(self, stat_name: str, values: list | ndarray) -
             reference_id=self.reference_id,
         )
 
+    def _to_array(self) -> ndarray:
+        """convert alignment to numpy array"""
+        return np.array([list(self.alignment[seq_id]) for seq_id in self.sequence_ids])
+
+    def _get_reference_seq(self) -> str:
+        """get the sequence of the reference sequence or majority consensus"""
+        return self.alignment[self.reference_id] if self.reference_id is not None else self.get_consensus()
+
     def get_reference_coords(self) -> tuple[int, int]:
         """
         Determine the start and end coordinates of the reference sequence
@@ -385,7 +408,7 @@ def get_conserved_orfs(self, min_length: int = 100, identity_cutoff: float | Non
             - all ungapped seqs[start:stop] must have at least min_length
             - no ungapped seq can have a Stop in between Start Stop
 
-        Conservation is measured by number of positions with identical characters divided by
+        Conservation is measured by the number of positions with identical characters divided by
         orf slice of the alignment.
 
         **Algorithm overview:**
@@ -413,7 +436,8 @@ def determine_conserved_start_stops(alignment: dict, alignment_length: int) -> t
             stops = config.STOP_CODONS[self.aln_type]
 
             list_of_starts, list_of_stops = [], []
-            ref = alignment[list(alignment.keys())[0]]
+            # define one sequence (first) as reference (it does not matter which one)
+            ref = alignment[self.sequence_ids[0]]
             for nt_position in range(alignment_length):
                 if ref[nt_position:nt_position + 3] in starts:
                     conserved_start = True
@@ -706,7 +730,7 @@ def calc_coverage(self) -> AlignmentStats:
 
         for nuc_pos in range(self.length):
             pos = str()
-            for record in aln.keys():
+            for record in self.sequence_ids:
                 pos = pos + aln[record][nuc_pos]
             coverage.append(1 - pos.count('-') / len(pos))
 
@@ -746,8 +770,7 @@ def calc_numerical_alignment(self, encode_mask:bool=False, encode_ambiguities:bo
         :returns matrix
         """
 
-        aln = self.alignment
-        sequences = np.array([list(aln[seq_id]) for seq_id in list(aln.keys())])
+        sequences = self._to_array()
         # ini matrix
         numerical_matrix = np.full(sequences.shape, np.nan, dtype=float)
         # first encode mask
@@ -778,11 +801,10 @@ def calc_identity_alignment(self, encode_mismatches:bool=True, encode_mask:bool=
         :return: identity alignment
         """
 
-        aln = self.alignment
-        ref = aln[self.reference_id] if self.reference_id is not None else self.get_consensus()
+        ref = self._get_reference_seq()
 
         # convert alignment to array
-        sequences = np.array([list(aln[seq_id]) for seq_id in list(aln.keys())])
+        sequences = self._to_array()
         reference = np.array(list(ref))
         # ini matrix
         identity_matrix = np.full(sequences.shape, 0, dtype=float)
@@ -861,8 +883,8 @@ def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=Tr
             If the specified substitution matrix is not available for the given alignment type.
         """
 
-        aln = self.alignment
-        ref = aln[self.reference_id] if self.reference_id is not None else self.get_consensus()
+        ref = self._get_reference_seq()
+
         if matrix_type is None:
             if self.aln_type == 'AA':
                 matrix_type = 'BLOSUM65'
@@ -878,7 +900,7 @@ def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=Tr
 
         # set dtype and convert alignment to a NumPy array for vectorized processing
         dtype = np.dtype(float, metadata={'matrix': matrix_type})
-        sequences = np.array([list(aln[seq_id]) for seq_id in list(aln.keys())])
+        sequences = self._to_array()
         reference = np.array(list(ref))
         valid_chars = list(subs_matrix.keys())
         similarity_array = np.full(sequences.shape, np.nan, dtype=dtype)
@@ -916,11 +938,10 @@ def calc_position_matrix(self, matrix_type:str='PWM') -> None | ndarray | ValueE
         """
 
         # ini
-        aln = self.alignment
         if matrix_type not in ['PFM', 'PPM', 'IC', 'PWM']:
             raise ValueError('Matrix_type must be PFM, PPM, IC or PWM.')
         possible_chars = list(config.CHAR_COLORS[self.aln_type]['standard'].keys())
-        sequences = np.array([list(aln[seq_id]) for seq_id in list(aln.keys())])
+        sequences = self._to_array()
 
         # calc position frequency matrix
         pfm = np.array([np.sum(sequences == char, 0) for char in possible_chars])
@@ -931,7 +952,7 @@ def calc_position_matrix(self, matrix_type:str='PWM') -> None | ndarray | ValueE
         pseudo_count = 0.0001  # to avoid 0 values
         pfm = pfm + pseudo_count
         ppm_non_char_excluded = pfm/np.sum(pfm, axis=0)  # use this for pwm/ic calculation
-        ppm = pfm/len(aln.keys())  # calculate the frequency based on row number
+        ppm = pfm/len(self.sequence_ids)  # calculate the frequency based on row number
         if matrix_type == 'PPM':
             return ppm
 
@@ -959,11 +980,7 @@ def calc_percent_recovery(self) -> dict:
         """
 
         aln = self.alignment
-
-        if self.reference_id is not None:
-            ref = aln[self.reference_id]
-        else:
-            ref = self.get_consensus()  # majority consensus
+        ref = self._get_reference_seq()
 
         if not any(char != '-' for char in ref):
             raise ValueError("Reference sequence is entirely gapped, cannot calculate recovery.")
@@ -977,7 +994,7 @@ def calc_percent_recovery(self) -> dict:
         cumulative_length = len(non_gap_positions)
 
         # Calculate recovery
-        for seq_id in aln:
+        for seq_id in self.sequence_ids:
             if seq_id == self.reference_id:
                 continue
             seq = aln[seq_id]
@@ -1099,7 +1116,7 @@ def calc_pairwise_identity_matrix(self, distance_type:str='ghd') -> PairwiseDist
 
         return PairwiseDistance(
             reference_id=None,
-            sequence_ids=list(aln.keys()),
+            sequence_ids=self.sequence_ids,
             distances=distance_matrix
         )
 
@@ -1174,8 +1191,8 @@ def get_snps(self, include_ambig:bool=False) -> dict:
         :return: dictionary containing snp positions and their variants including their frequency.
         """
         aln = self.alignment
-        ref = aln[self.reference_id] if self.reference_id is not None else self.get_consensus()
-        aln = {x: aln[x] for x in aln.keys() if x != self.reference_id}
+        ref = self._get_reference_seq()
+        aln = {x: aln[x] for x in self.sequence_ids if x != self.reference_id}
         seq_ids = list(aln.keys())
         snp_dict = {'#CHROM': self.reference_id if self.reference_id is not None else 'consensus', 'POS': {}}
 
@@ -1185,7 +1202,7 @@ def get_snps(self, include_ambig:bool=False) -> dict:
                 if reference_char in config.AMBIG_CHARS[self.aln_type] and reference_char != '-':
                     continue
             alt_chars, snps = [], []
-            for i, seq_id in enumerate(aln.keys()):
+            for i, seq_id in enumerate(seq_ids):
                 alt_chars.append(aln[seq_id][pos])
                 if reference_char != aln[seq_id][pos]:
                     snps.append(i)

From 2c33483397857e08176864e80a91f154e3becbae Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Thu, 19 Mar 2026 01:31:03 +0100
Subject: [PATCH 09/25] finalized dunder methods intergration

---
 app_src/shiny_plots.py          |  4 +-
 app_src/shiny_server.py         | 18 ++++----
 msaexplorer/__init__.py         |  2 +-
 msaexplorer/_helpers.py         | 58 ++++++++++++++++++++++++++
 msaexplorer/draw.py             | 39 +++++++++---------
 msaexplorer/explore.py          | 73 ++++++---------------------------
 msaexplorer/export.py           |  2 +-
 tests/README.md                 |  2 +-
 tests/test_alignment_parsing.py |  9 ++--
 9 files changed, 108 insertions(+), 99 deletions(-)

diff --git a/app_src/shiny_plots.py b/app_src/shiny_plots.py
index 585df33..066306f 100644
--- a/app_src/shiny_plots.py
+++ b/app_src/shiny_plots.py
@@ -13,7 +13,7 @@ def set_aln(aln, inputs):
     # set the reference sequence
     if 'reference' in inputs:
         if inputs['reference'] == 'first':
-            aln.reference_id = list(aln.alignment.keys())[0]
+            aln.reference_id = next(iter(aln))
         elif inputs['reference'] == 'consensus':
             aln.reference_id = None
         else:
@@ -200,7 +200,7 @@ def create_analysis_custom_heatmap(aln, inputs):
         figure_size = int(inputs['dimensions']['width'] * 0.7)
 
     matrix = aln.calc_pairwise_identity_matrix(inputs['additional_analysis_options_left']).distances
-    labels = [x.split(' ')[0] for x in list(aln.alignment.keys())]
+    labels = [x.split(' ')[0] for x in list(aln)]
 
     # generate hover text
     hover_text = [
diff --git a/app_src/shiny_server.py b/app_src/shiny_server.py
index 3dc644d..a806239 100644
--- a/app_src/shiny_server.py
+++ b/app_src/shiny_server.py
@@ -128,7 +128,7 @@ def prepare_inputs():
                 complete_size += input.plot_1_size()
             if inputs['annotation'] != 'Off':
                 complete_size += input.plot_3_size()
-            relative_msa_height = input.plot_2_size() * increase_height / complete_size * window_height / len(aln.alignment)
+            relative_msa_height = input.plot_2_size() * increase_height / complete_size * window_height / len(aln)
             try:
                 relative_msa_width = window_width / (inputs['zoom_range'][1] - inputs['zoom_range'][0])
             except ZeroDivisionError:
@@ -292,10 +292,10 @@ def finalize_loaded_alignment(aln, annotation_file):
         All the necessary steps to load an alignment into the app are done here.
         """
 
-        aln.reference_id = list(aln.alignment.keys())[0]
+        aln.reference_id = next(iter(aln))
         reactive.alignment.set(aln)
 
-        alignment_length = len(next(iter(aln.alignment.values()))) - 1
+        alignment_length = aln.length - 1
         ui.update_slider('zoom_range', max=alignment_length - 1, value=(0, alignment_length - 1))
 
         ui.remove_ui(selector="#orf_column")
@@ -315,7 +315,7 @@ def finalize_loaded_alignment(aln, annotation_file):
             )
 
         for id in ['reference', 'reference_2']:
-            ui.update_selectize(id=id, choices=['first', 'consensus'] + list(aln.alignment.keys()), selected='first')
+            ui.update_selectize(id=id, choices=['first', 'consensus'] + list(aln), selected='first')
 
         ui.update_selectize(
             id='matrix',
@@ -323,7 +323,7 @@ def finalize_loaded_alignment(aln, annotation_file):
             selected='BLOSUM65' if aln.aln_type == 'AA' else 'TRANS',
         )
 
-        aln_len, seq_threshold = len(aln.alignment.keys()), 5
+        aln_len, seq_threshold = len(aln), 5
         for ratio in config.STANDARD_HEIGHT_RATIOS.keys():
             if aln_len >= ratio:
                 seq_threshold = ratio
@@ -634,7 +634,7 @@ def update_download_options():
                 where='beforeBegin'
             ) if aln is None else ui.insert_ui(
                 ui.input_selectize(
-                    id='reference_2', label='Reference', choices=['first', 'consensus'] + list(aln.alignment.keys()), selected='first'
+                    id='reference_2', label='Reference', choices=['first', 'consensus'] + list(aln), selected='first'
                 ),
                 selector='#download_format-label',
                 where='beforeBegin'
@@ -689,7 +689,7 @@ def download_stats():
         # helper functions
         def _snp_option():
             if input.reference_2() == 'first':
-                aln.reference_id = list(aln.alignment.keys())[0]
+                aln.reference_id = next(iter(aln))
             elif input.reference_2() == 'consensus':
                 aln.reference_id = None
             else:
@@ -776,7 +776,7 @@ def _char_freq_option():
 
         def _percent_recovery_option():
             if input.reference_2() == 'first':
-                aln.reference_id = list(aln.alignment.keys())[0]
+                aln.reference_id = next(iter(aln))
             elif input.reference_2() == 'consensus':
                 aln.reference_id = None
             else:
@@ -879,7 +879,7 @@ def number_of_seq():
         if aln is None:
             return None
 
-        return len(aln.alignment)
+        return len(aln)
 
     @render.ui
     def aln_len():
diff --git a/msaexplorer/__init__.py b/msaexplorer/__init__.py
index a53d13d..165af04 100644
--- a/msaexplorer/__init__.py
+++ b/msaexplorer/__init__.py
@@ -296,7 +296,7 @@
 annotation = explore.Annotation(msa, 'annotation.gff3')
 
 # Set parameters
-msa.reference_id = list(msa.alignment.keys())[0]
+msa.reference_id = next(iter(msa))
 msa.zoom = (0, 2000)
 
 # Compute statistics
diff --git a/msaexplorer/_helpers.py b/msaexplorer/_helpers.py
index acc66d6..ac70f15 100644
--- a/msaexplorer/_helpers.py
+++ b/msaexplorer/_helpers.py
@@ -4,8 +4,12 @@
 
 import os, io, math
 import numpy as np
+from msaexplorer import config
 from typing import Callable, Dict
 from matplotlib.colors import is_color_like
+from Bio.Align import MultipleSeqAlignment
+from Bio import AlignIO
+
 
 # explore helpers
 def _get_line_iterator(source):
@@ -18,6 +22,60 @@ def _get_line_iterator(source):
         return io.StringIO(source)
 
 
+def _read_alignment(source: str | MultipleSeqAlignment) -> dict:
+    """
+    Parse MSA alignment using Biopython with automatic format detection.
+    :param source: Path to alignment file, raw alignment string, or Bio.Align.MultipleSeqAlignment object
+    :possible_chars: list of possible characters in the alignment
+    :return: dictionary with ids as keys and sequences as values
+    """
+    # Handle Bio.Align.MultipleSeqAlignment objects
+    if isinstance(source, MultipleSeqAlignment):
+        aln_dict = {record.id: str(record.seq).upper() for record in source}
+    else:
+        # Try multiple formats in order of likelihood
+        formats_to_try = ["fasta", "clustal", "phylip", "stockholm", "nexus"]
+
+        aln_dict = None
+        for fmt in formats_to_try:
+            try:
+                with _get_line_iterator(source) as handle:
+                    alignment = AlignIO.read(handle, fmt)
+                    aln_dict = {record.id: str(record.seq).upper() for record in alignment}
+                    break
+            except Exception:
+                continue
+
+        if aln_dict is None:
+            # If no format worked, raise an error
+            raise ValueError(f"Alignment file could not be parsed. Supported formats: {', '.join(formats_to_try)}")
+
+    # Validate alignment
+    if not aln_dict:
+        raise ValueError(f"Alignment does not contain any sequences.")
+
+    if len(aln_dict) < 2:
+        raise ValueError("Alignment must contain more than one sequence.")
+
+    # Check for non-allowed characters
+    for sequence_id, seq in aln_dict.items():
+        invalid_chars = set(seq) - set(config.POSSIBLE_CHARS)
+        if invalid_chars:
+            raise ValueError(
+                f"{sequence_id} contains invalid characters: {', '.join(invalid_chars)}. Allowed chars are: {config.POSSIBLE_CHARS}."
+            )
+
+    # Validate all sequences have same length
+    first_seq_len = len(next(iter(aln_dict.values())))
+    for sequence_id, seq in aln_dict.items():
+        if len(seq) != first_seq_len:
+            raise ValueError(
+                f"All alignment sequences must have the same length. Sequence '{sequence_id}' has length {len(seq)}, expected {first_seq_len}."
+            )
+
+    return aln_dict
+
+
 def _create_distance_calculation_function_mapping() -> Dict[str, Callable[[str, str, int], float]]:
     """
     create a mapping of distance types to distance calculation functions
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 80d3741..e731a6f 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -132,20 +132,20 @@ def _seq_names(aln: explore.MSA, ax: plt.Axes, custom_seq_names: tuple, show_seq
         show_seq_names = True
         if not isinstance(custom_seq_names, tuple):
             raise ValueError('configure your custom names list: custom_names=(name1, name2...)')
-        if len(custom_seq_names) != len(aln.alignment.keys()):
+        if len(custom_seq_names) != len(aln):
             raise ValueError('length of sequences not equal to number of custom names')
     if show_seq_names:
         ax.yaxis.set_ticks_position('none')
-        ax.set_yticks(np.arange(len(aln.alignment)))
+        ax.set_yticks(np.arange(len(aln)))
         if custom_seq_names:
             names = custom_seq_names[::-1]
         else:
-            names = [x.split(' ')[0] for x in list(aln.alignment.keys())[::-1]]
+            names = [x.split(' ')[0] for x in aln.sequence_ids[::-1]]
         if include_consensus:
             names = names + ['consensus']
-            y_ticks = np.arange(len(aln.alignment) + 1)
+            y_ticks = np.arange(len(aln) + 1)
         else:
-            y_ticks = np.arange(len(aln.alignment))
+            y_ticks = np.arange(len(aln))
         ax.set_yticks(y_ticks)
         ax.set_yticklabels(names)
     else:
@@ -318,7 +318,7 @@ def _plot_sequence_text(aln: explore.MSA, seq_name: str, ref_name: str | None, a
         else:
             different_cols = [False] * aln.length
 
-        for idx, (character, value) in enumerate(zip(aln.alignment[seq_name], values)):
+        for idx, (character, value) in enumerate(zip(aln[seq_name], values)):
             if value != value_to_skip and character != '-' or seq_name == ref_name and character != '-' or character == '-' and not show_gaps or always_text and character != '-':
                 if seq_name == ref_name and ref_color is not None:
                     text_color = _get_contrast_text_color(to_rgba(ref_color))
@@ -343,9 +343,9 @@ def _plot_sequence_text(aln: explore.MSA, seq_name: str, ref_name: str | None, a
     polygons, polygon_colors, patch_list = [], [], []
     # determine zoom
     zoom = (0, aln.length) if aln.zoom is None else aln.zoom
-    for i, seq_name in enumerate(aln.alignment):
+    for i, seq_name in enumerate(aln):
         # define initial y position
-        y_position = len(aln.alignment) - i - 1.4
+        y_position = len(aln) - i - 1.4
         # now plot relevant stuff for the current row
         row = matrix[i]
         # plot a line below everything for fancy gaps
@@ -376,7 +376,7 @@ def _plot_sequence_text(aln: explore.MSA, seq_name: str, ref_name: str | None, a
         # add sequence text
         if show_different_sequence or show_sequence_all:
             _plot_sequence_text(
-                aln=aln, seq_name=list(aln.alignment.keys())[i], ref_name=aln.reference_id, always_text=show_sequence_all,
+                aln=aln, seq_name=seq_name, ref_name=aln.reference_id, always_text=show_sequence_all,
                 values=matrix[i], matrix=matrix, ax=ax, zoom=zoom, y_position=y_position, value_to_skip=identical_value, ref_color=reference_color,
                 show_gaps=show_gaps, aln_colors=aln_colors
             )
@@ -451,9 +451,9 @@ def alignment(aln: explore.MSA | str, ax: plt.Axes | None = None, show_sequence_
         consensus_plot(aln=aln, ax=ax, show_x_label=show_x_label, show_name=False, show_sequence=show_sequence_all,
             color_scheme='standard', basic_color=basic_color, mask_color=mask_color, ambiguity_color=ambiguity_color
         )
-        ax.set_ylim(-0.5, len(aln.alignment) + 1)
+        ax.set_ylim(-0.5, len(aln) + 1)
     else:
-        ax.set_ylim(-0.5, len(aln.alignment))
+        ax.set_ylim(-0.5, len(aln))
 
     _seq_names(aln=aln, ax=ax, custom_seq_names=custom_seq_names, show_seq_names=show_seq_names, include_consensus=show_consensus)
     # configure axis
@@ -537,9 +537,9 @@ def identity_alignment(aln: explore.MSA | str, ax: plt.Axes | None = None, show_
                        color_scheme='standard', basic_color=basic_color, mask_color=mask_color,
                        ambiguity_color=ambiguity_color
                        )
-        ax.set_ylim(-0.5, len(aln.alignment) + 1)
+        ax.set_ylim(-0.5, len(aln) + 1)
     else:
-        ax.set_ylim(-0.5, len(aln.alignment))
+        ax.set_ylim(-0.5, len(aln))
 
     _seq_names(aln=aln, ax=ax, custom_seq_names=custom_seq_names, show_seq_names=show_seq_names,
                include_consensus=show_consensus)
@@ -619,9 +619,9 @@ def similarity_alignment(aln: explore.MSA | str, ax: plt.Axes | None = None, mat
     if show_consensus:
         consensus_plot(aln=aln, ax=ax, show_x_label=show_x_label, show_name=False, show_sequence=any([show_sequence_all, show_similarity_sequence]),
                        color_scheme='standard', basic_color=basic_color)
-        ax.set_ylim(-0.5, len(aln.alignment) + 1)
+        ax.set_ylim(-0.5, len(aln) + 1)
     else:
-        ax.set_ylim(-0.5, len(aln.alignment))
+        ax.set_ylim(-0.5, len(aln))
 
     _seq_names(aln=aln, ax=ax, custom_seq_names=custom_seq_names, show_seq_names=show_seq_names,
                include_consensus=show_consensus)
@@ -851,7 +851,6 @@ def variant_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, lollisize:
     return ax
 
 
-
 def _plot_annotation(annotation_dict: dict, ax: plt.Axes, direction_marker_size: int | None, color: str | ScalarMappable):
     """
     Plot annotation rectangles
@@ -1130,7 +1129,7 @@ def consensus_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, threshold
 
     zoom_offset = aln.zoom[0] if aln.zoom is not None else 0
 
-    y_position = len(aln.alignment) - 0.4
+    y_position = len(aln) - 0.4
 
     for pos, char in enumerate(consensus):
         x = pos + zoom_offset
@@ -1216,7 +1215,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
         raise ValueError('window_size has to be a positive integer')
     if window_size > aln.length:
         raise ValueError('window_size can not be larger than the (zoomed) alignment length')
-    if ref is not None and ref not in aln.alignment:
+    if ref is not None and ref not in aln:
         raise ValueError(f'Reference {ref} not in alignment')
     if distance_calculation not in ['ghd', 'ged', 'jc69', 'k2p']:
         raise ValueError(f'Distance calculation method {distance_calculation} not supported. Supported: ghd, ged, jc69, k2p')
@@ -1224,7 +1223,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
         raise ValueError(f'Distance calculation method {distance_calculation} only supported for nucleotide alignments')
 
     # get the sequence ids to plot
-    sequence_ids = [key for key in aln.alignment.keys() if key != ref]
+    sequence_ids = [seq_id for seq_id in aln if seq_id != ref]
 
     # validate colors
     if colors is not None:
@@ -1272,7 +1271,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
         # slice alignment and replace the original alignment
         aln_tmp._alignment = {
             seq_id: seq[left_side:right_side]
-            for seq_id, seq in aln.alignment.items()
+            for seq_id, seq in aln.items()
         }
         window_result = aln_tmp.calc_pairwise_distance_to_reference(distance_type=distance_calculation)
         value_map = {seq_id: value for seq_id, value in zip(window_result.sequence_ids, window_result.distances)}
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 7eece3a..1aecbe4 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -25,7 +25,7 @@
 # msaexplorer
 from msaexplorer import config
 from msaexplorer._data_classes import PairwiseDistance, AlignmentStats
-from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping
+from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping, _read_alignment
 
 #TODO: Move outputs to dataclasses
 class MSA:
@@ -43,75 +43,24 @@ def __init__(self, alignment_string: str | MultipleSeqAlignment, reference_id: s
         :param reference_id: reference id
         :param zoom_range: start and stop positions to zoom into the alignment
         """
-        self._alignment = self._read_alignment(alignment_string)
+        self._alignment = _read_alignment(alignment_string)
         self._reference_id = self._validate_ref(reference_id, self._alignment)
         self._zoom = self._validate_zoom(zoom_range, self._alignment)
         self._aln_type = self._determine_aln_type(self._alignment)
 
     def __len__(self) -> int:
-        return len(next(iter(self.alignment.values())))
+        return len(self.sequence_ids)
 
     def __getitem__(self, key: str) -> str:
         return self.alignment[key]
 
     def __contains__(self, item: str) -> bool:
-        """Seq ID present"""
         return item in self.sequence_ids
 
-    # Static methods
-    @staticmethod
-    def _read_alignment(source: str | MultipleSeqAlignment) -> dict:
-        """
-        Parse MSA alignment using Biopython with automatic format detection.
-        :param source: Path to alignment file, raw alignment string, or Bio.Align.MultipleSeqAlignment object
-        :return: dictionary with ids as keys and sequences as values
-        """
-        # Handle Bio.Align.MultipleSeqAlignment objects
-        if isinstance(source, MultipleSeqAlignment):
-            aln_dict = {record.id: str(record.seq).upper() for record in source}
-        else:
-            # Try multiple formats in order of likelihood
-            formats_to_try = ["fasta", "clustal", "phylip", "stockholm", "nexus"]
-
-            aln_dict = None
-            for fmt in formats_to_try:
-                try:
-                    with _get_line_iterator(source) as handle:
-                        alignment = AlignIO.read(handle, fmt)
-                        aln_dict = {record.id: str(record.seq).upper() for record in alignment}
-                        break
-                except Exception:
-                    continue
-
-            if aln_dict is None:
-                # If no format worked, raise an error
-                raise ValueError(f"Alignment file could not be parsed. Supported formats: {', '.join(formats_to_try)}")
-
-        # Validate alignment
-        if not aln_dict:
-            raise ValueError(f"Alignment does not contain any sequences.")
-
-        if len(aln_dict) < 2:
-            raise ValueError("Alignment must contain more than one sequence.")
-
-        # Check for non-allowed characters
-        for sequence_id, seq in aln_dict.items():
-            invalid_chars = set(seq) - set(config.POSSIBLE_CHARS)
-            if invalid_chars:
-                raise ValueError(
-                    f"{sequence_id} contains invalid characters: {', '.join(invalid_chars)}. Allowed chars are: {config.POSSIBLE_CHARS}"
-                )
-
-        # Validate all sequences have same length
-        first_seq_len = len(next(iter(aln_dict.values())))
-        for sequence_id, seq in aln_dict.items():
-            if len(seq) != first_seq_len:
-                raise ValueError(
-                    f"All alignment sequences must have the same length. Sequence '{sequence_id}' has length {len(seq)}, expected {first_seq_len}."
-                )
-
-        return aln_dict
+    def __iter__(self):
+        return iter(self.alignment)
 
+    # Static methods
     @staticmethod
     def _validate_ref(reference: str | None, alignment: dict) -> str | None | ValueError:
         """
@@ -229,7 +178,6 @@ def alignment(self) -> dict:
         else:
             return self._alignment
 
-    # functions for different alignment stats
     def _create_position_stat_result(self, stat_name: str, values: list | ndarray) -> AlignmentStats:
         """
         Build a shared dataclass for position-wise statistics.
@@ -249,20 +197,23 @@ def _create_position_stat_result(self, stat_name: str, values: list | ndarray) -
         )
 
     def _to_array(self) -> ndarray:
-        """convert alignment to numpy array"""
+        """convert alignment to a numpy array"""
         return np.array([list(self.alignment[seq_id]) for seq_id in self.sequence_ids])
 
     def _get_reference_seq(self) -> str:
         """get the sequence of the reference sequence or majority consensus"""
         return self.alignment[self.reference_id] if self.reference_id is not None else self.get_consensus()
 
+    def items(self):
+        return self.alignment.items()
+
     def get_reference_coords(self) -> tuple[int, int]:
         """
         Determine the start and end coordinates of the reference sequence
         defined as the first/last nucleotide in the reference sequence
         (excluding N and gaps).
 
-        :return: start, end
+        :return: Start, End
         """
         start, end = 0, self.length
 
@@ -1508,7 +1459,7 @@ def parse_bed(file_path) -> dict:
         # sanity check whether one of the annotation ids and alignment ids match
         annotation_found = False
         for annotation in annotations.keys():
-            for aln_id in aln.alignment.keys():
+            for aln_id in aln:
                 aln_id_sanitized = aln_id.split(' ')[0]
                 # check in both directions
                 if aln_id_sanitized in annotation:
diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index c85282f..c203321 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -12,7 +12,7 @@
 from msaexplorer._data_classes import AlignmentStats
 from msaexplorer._helpers import _check_and_create_path
 
-
+#TODO Include tests for export functions
 def snps(snp_dict: dict, format_type: str = 'vcf', path: str | None = None) -> str | None | ValueError:
     """
     Export a SNP dictionary to a VCF or tabular format. Importantly, the input dictionary has to be in the standard
diff --git a/tests/README.md b/tests/README.md
index d4578bd..38099c5 100644
--- a/tests/README.md
+++ b/tests/README.md
@@ -1,6 +1,6 @@
 # Tests
 
-This test suite validates the `msaexplorer.explore` module.
+This test suite validates the `msaexplorer` module.
 
 Run with:
 
diff --git a/tests/test_alignment_parsing.py b/tests/test_alignment_parsing.py
index d6c0838..3e26360 100644
--- a/tests/test_alignment_parsing.py
+++ b/tests/test_alignment_parsing.py
@@ -5,6 +5,7 @@
 from Bio.Align import MultipleSeqAlignment
 
 from msaexplorer.explore import MSA
+from msaexplorer._helpers import _read_alignment
 
 DATA_DIR = pathlib.Path(__file__).parent / "data"
 
@@ -19,7 +20,7 @@
 def test_read_alignment_parses_multiple_formats_from_file(format_name: str):
     """Test that MSA can read different alignment formats from a file."""
     alignment_file = DATA_DIR / f"alignment.{format_name}"
-    parsed_alignment = MSA._read_alignment(str(alignment_file))
+    parsed_alignment = _read_alignment(str(alignment_file))
 
     assert parsed_alignment == EXPECTED_ALIGNMENT
 
@@ -29,7 +30,7 @@ def test_read_alignment_parses_multiple_formats_from_string(format_name: str):
     """Test that MSA can read different alignment formats from a raw string."""
     alignment_file = DATA_DIR / f"alignment.{format_name}"
     alignment_content = alignment_file.read_text(encoding="utf-8")
-    parsed_alignment = MSA._read_alignment(alignment_content)
+    parsed_alignment = _read_alignment(alignment_content)
 
     assert parsed_alignment == EXPECTED_ALIGNMENT
 
@@ -44,7 +45,7 @@ def test_read_alignment_accepts_bio_alignment_object():
     bio_alignment = MultipleSeqAlignment(records)
 
     # Pass it directly to MSA._read_alignment
-    parsed_alignment = MSA._read_alignment(bio_alignment)
+    parsed_alignment = _read_alignment(bio_alignment)
 
     assert parsed_alignment == EXPECTED_ALIGNMENT
 
@@ -69,4 +70,4 @@ def test_msa_initializes_with_bio_alignment_object():
 def test_read_alignment_raises_for_unparseable_content() -> None:
     """Test that MSA raises ValueError when given unparseable content."""
     with pytest.raises(ValueError, match="could not be parsed"):
-        MSA._read_alignment("this is not an alignment")
+        _read_alignment("this is not an alignment")

From 0c2ef3c41046f5c598da937fb1249646dec222c5 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Thu, 19 Mar 2026 10:43:12 +0100
Subject: [PATCH 10/25] added tests for output functions

---
 msaexplorer/export.py |   2 +-
 tests/test_export.py  | 214 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 215 insertions(+), 1 deletion(-)
 create mode 100644 tests/test_export.py

diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index c203321..c85282f 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -12,7 +12,7 @@
 from msaexplorer._data_classes import AlignmentStats
 from msaexplorer._helpers import _check_and_create_path
 
-#TODO Include tests for export functions
+
 def snps(snp_dict: dict, format_type: str = 'vcf', path: str | None = None) -> str | None | ValueError:
     """
     Export a SNP dictionary to a VCF or tabular format. Importantly, the input dictionary has to be in the standard
diff --git a/tests/test_export.py b/tests/test_export.py
new file mode 100644
index 0000000..ab18877
--- /dev/null
+++ b/tests/test_export.py
@@ -0,0 +1,214 @@
+"""Tests for export helpers in ``msaexplorer.export``."""
+
+import numpy as np
+import pytest
+
+from msaexplorer import export
+from msaexplorer._data_classes import AlignmentStats
+
+
+@pytest.fixture
+def snp_dict():
+    """Returns a SNP dictionary for testing."""
+    return {
+        "#CHROM": "ref",
+        "POS": {
+            3: {
+                "ref": "C",
+                "ALT": {
+                    "T": {"AF": 1.0, "SEQ_ID": ["s1"]},
+                },
+            },
+            1: {
+                "ref": "A",
+                "ALT": {
+                    "G": {"AF": 0.5, "SEQ_ID": ["s1", "s2"]},
+                    "-": {"AF": 0.5, "SEQ_ID": ["s3"]},
+                },
+            },
+        },
+    }
+
+
+class TestSnpsExport:
+    def test_vcf_output_is_correct_and_sorted(self, snp_dict):
+        """Test that the VCF has correct vcf format."""
+        result = export.snps(snp_dict, format_type="vcf")
+
+        assert result == "\n".join(
+            [
+                "##fileformat=VCFv4.2",
+                "##source=MSAexplorer",
+                "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
+                "ref\t2\t.\tA\tG,-\t.\t.\tAF=0.5,0.5;SEQ_ID=s1|s2,s3",
+                "ref\t4\t.\tC\tT\t.\t.\tAF=1.0;SEQ_ID=s1",
+            ]
+        )
+
+    def test_tabular_output_is_correct(self, snp_dict):
+        """Test that the VCF has correct tabular format."""
+        result = export.snps(snp_dict, format_type="tabular")
+
+        assert result == "\n".join(
+            [
+                "CHROM\tPOS\tREF\tALT\tAF\tSEQ_ID",
+                "ref\t2\tA\tG\t0.5\ts1,s2",
+                "ref\t2\tA\t-\t0.5\ts3",
+                "ref\t4\tC\tT\t1.0\ts1",
+            ]
+        )
+
+    def test_invalid_input_raises_value_error(self):
+        """Test that invalid input raises a ValueError."""
+        with pytest.raises(ValueError, match="must be a dictionary"):
+            export.snps([], format_type="vcf")
+
+        with pytest.raises(ValueError, match="Missing required key"):
+            export.snps({"POS": {}}, format_type="vcf")
+
+        with pytest.raises(ValueError, match="Invalid format_type"):
+            export.snps({"#CHROM": "ref", "POS": {}}, format_type="csv")
+
+    def test_file_export_creates_expected_extension(self, snp_dict, tmp_path):
+        """Test that the file extension is correct when written."""
+        path_without_ext = tmp_path / "nested" / "snps_output"
+
+        result = export.snps(snp_dict, format_type="vcf", path=str(path_without_ext))
+
+        assert result is None
+        written = path_without_ext.with_suffix(path_without_ext.suffix + ".vcf")
+        assert written.exists()
+        assert "##fileformat=VCFv4.2" in written.read_text()
+
+
+class TestFastaExport:
+    def test_single_sequence_to_string(self):
+        """Test that a single sequence is exported correctly."""
+        result = export.fasta("ACGT", header="s1")
+
+        assert result == ">s1\nACGT"
+
+    def test_dictionary_sequence_to_string(self):
+        """Test that a dictionary of sequences is exported correctly."""
+        result = export.fasta({"s1": "ACGT", "s2": "A-GT"})
+
+        assert result == ">s1\nACGT\n>s2\nA-GT"
+
+    def test_invalid_sequence_raises(self):
+        """Test that invalid sequences raise a ValueError."""
+        with pytest.raises(ValueError, match="invalid characters"):
+            export.fasta("ACGTZ", header="s1")
+
+    def test_file_export_writes_content(self, tmp_path):
+        """Test that the file is written correctly."""
+        out = tmp_path / "subdir" / "seq.fasta"
+
+        result = export.fasta("ACGT", header="s1", path=str(out))
+
+        assert result is None
+        assert out.exists()
+        assert out.read_text() == ">s1\nACGT"
+
+
+class TestStatsExport:
+    def test_alignment_stats_dataclass_output(self):
+        """Test that the AlignmentStats dataclass is exported correctly."""
+        stat_data = AlignmentStats(
+            stat_name="entropy",
+            positions=np.array([2, 4, 6]),
+            values=np.array([0.1, 0.2, 0.3]),
+            aln_type="dna",
+            reference_id=None,
+        )
+
+        result = export.stats(stat_data, seperator=",")
+
+        assert result == "\n".join([
+            "position,value",
+            "2,0.1",
+            "4,0.2",
+            "6,0.3",
+        ])
+
+    def test_plain_array_output_uses_zero_based_positions(self):
+        """Test that plain arrays are exported correctly."""
+        result = export.stats([10, 20, 30])
+
+        assert result == "\n".join([
+            "position\tvalue",
+            "0\t10",
+            "1\t20",
+            "2\t30",
+        ])
+
+
+class TestOrfExport:
+    def test_orf_output_is_correct(self):
+        """Test that the ORF output is correct."""
+        orf_dict = {
+            "orf_1": {
+                "location": [(10, 50)],
+                "frame": 1,
+                "strand": "+",
+                "conservation": 87.125,
+                "internal": False,
+            }
+        }
+
+        result = export.orf(orf_dict, chrom="chr1")
+
+        assert result == "chr1\t10\t50\torf_1\t87.12\t+"
+
+    def test_orf_invalid_input_raises(self):
+        """Test that invalid input raises a ValueError."""
+        with pytest.raises(ValueError, match="empty"):
+            export.orf({}, chrom="chr1")
+
+
+class TestCharacterFreqExport:
+    def test_character_freq_output_is_correct(self):
+        """Test that the character frequency output is correct."""
+        data = {
+            "total": {},
+            "seq1": {
+                "A": {"counts": 2, "% of alignment": 50.0, "% of non-gapped": 66.67},
+                "C": {"counts": 1, "% of alignment": 25.0, "% of non-gapped": 33.33},
+                "-": {"counts": 1, "% of alignment": 25.0, "% of non-gapped": 0.0},
+            },
+        }
+
+        result = export.character_freq(data, seperator=",")
+
+        assert result == "\n".join([
+            "sequence,char,counts,% of non-gapped",
+            "seq1,A,2,66.67",
+            "seq1,C,1,33.33",
+        ])
+
+    def test_character_freq_invalid_char_raises(self):
+        """Test that invalid characters raise a ValueError."""
+        bad = {"seq1": {"?": {"counts": 1, "% of alignment": 10.0, "% of non-gapped": 10.0}}}
+
+        with pytest.raises(ValueError, match="invalid"):
+            export.character_freq(bad)
+
+
+class TestPercentRecoveryExport:
+    def test_percent_recovery_output_is_correct(self):
+        """Test that the percent recovery output is correct."""
+        rec = {"seq1": 75.0, "seq2": 90.5}
+
+        result = export.percent_recovery(rec)
+
+        assert result == "\n".join([
+            "sequence\t% recovery",
+            "seq1\t75.0",
+            "seq2\t90.5",
+        ])
+
+    def test_percent_recovery_value_must_be_float(self):
+        """Test that the percent recovery value must be a float."""
+        with pytest.raises(ValueError, match="invalid"):
+            export.percent_recovery({"seq1": 75})
+
+

From 56dccaed5b846f45848939bedc3455e3a55b7751 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Thu, 19 Mar 2026 17:21:56 +0100
Subject: [PATCH 11/25] added orf class

---
 msaexplorer/_data_classes.py | 117 ++++++++++++++++++++++++++++------
 msaexplorer/draw.py          |  19 +++++-
 msaexplorer/explore.py       |  80 ++++++++++++-----------
 msaexplorer/export.py        |  24 ++++---
 tests/test_export.py         |  35 +++++------
 tests/test_orf_detection.py  | 119 ++++++++++++++++++-----------------
 6 files changed, 247 insertions(+), 147 deletions(-)

diff --git a/msaexplorer/_data_classes.py b/msaexplorer/_data_classes.py
index b64881c..847ad9b 100644
--- a/msaexplorer/_data_classes.py
+++ b/msaexplorer/_data_classes.py
@@ -3,7 +3,7 @@
 """
 
 # build-in
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 
 # libs
 from numpy import ndarray
@@ -18,8 +18,6 @@ class AlignmentStats:
     stat_name: str
     positions: ndarray
     values: ndarray
-    aln_type: str
-    reference_id: str | None
 
     def __post_init__(self):
         if self.positions.shape != self.values.shape:
@@ -32,21 +30,9 @@ def __len__(self) -> int:
     def __getitem__(self, index: int) -> float:
         return self.values[index]
 
-    def __iter__(self):
-        yield self.stat_name
-        yield self.positions
-        yield self.values
-
     def __contains__(self, item: float) -> bool:
         return item in self.positions
 
-    # normal mehods
-    def as_array(self) -> ndarray:
-        return self.values
-
-    def as_list(self) -> list:
-        return self.values.tolist()
-
 
 @dataclass(frozen=True)
 class PairwiseDistance:
@@ -59,11 +45,6 @@ class PairwiseDistance:
     distances: ndarray
 
     # dunder methods
-    def __iter__(self):
-        yield self.reference_id
-        yield self.sequence_ids
-        yield self.distances
-
     def __len__(self) -> int:
         return len(self.sequence_ids)
 
@@ -81,3 +62,99 @@ def __getitem__(self, index: int | str) -> float:
     def __contains__(self, item: str) -> bool:
         """Seq ID present"""
         return item in self.sequence_ids
+
+
+@dataclass(frozen=True)
+class OpenReadingFrame:
+    """
+    Represents a single conserved ORF detected across an alignment.
+
+    Attributes:
+        orf_id:       Unique identifier, e.g. ``'ORF_0'``.
+        location:     Main ORF boundaries as a tuple of ``(start, stop)`` pairs
+                      (0-based, half-open).  Typically a single pair, but may
+                      carry additional coordinates for split ORFs.
+        frame:        Reading frame (0, 1, or 2).
+        strand:       ``'+'`` for forward, ``'-'`` for reverse complement.
+        conservation: Percentage of fully identical alignment columns inside the ORF.
+        internal:     Tuple of ``(start, stop)`` pairs for nested (internal) ORFs
+                      that share the same stop codon.
+    """
+
+    orf_id: str
+    location: tuple[tuple[int, int], ...]
+    frame: int
+    strand: str
+    conservation: float
+    internal: tuple[tuple[int, int], ...] = field(default_factory=tuple)
+
+    def __post_init__(self):
+        if self.strand not in ('+', '-'):
+            raise ValueError(f"strand must be '+' or '-', got {self.strand!r}")
+        if not (0 <= self.frame <= 2):
+            raise ValueError(f"frame must be 0, 1, or 2, got {self.frame!r}")
+
+    def __len__(self) -> int:
+        """Length of the main ORF in alignment columns."""
+        return self.location[0][1] - self.location[0][0]
+
+    def __contains__(self, position: int) -> bool:
+        """True if *position* (0-based) falls inside the main ORF."""
+        start, stop = self.location[0]
+        return start <= position < stop
+
+
+@dataclass(frozen=True)
+class OrfContainer:
+    """
+    Ordered collection of `OpenReadingFrame` objects returned by
+    `MSA.get_conserved_orfs` or`MSA.get_non_overlapping_conserved_orfs`.
+
+    The class intentionally mimics a *read-only dict* interface
+    """
+
+    orfs: tuple[OpenReadingFrame, ...] = field(default_factory=tuple)
+
+    # dict-like interface
+    def keys(self) -> list[str]:
+        """Return ORF identifiers in insertion order."""
+        return [orf.orf_id for orf in self.orfs]
+
+    def values(self) -> list[OpenReadingFrame]:
+        """Return :class:`OpenReadingFrame` objects in insertion order."""
+        return list(self.orfs)
+
+    def items(self) -> list[tuple[str, OpenReadingFrame]]:
+        """Return ``(orf_id, OpenReadingFrame)`` pairs in insertion order."""
+        return [(orf.orf_id, orf) for orf in self.orfs]
+
+    # dunder methods
+    def __len__(self) -> int:
+        return len(self.orfs)
+
+    def __bool__(self) -> bool:
+        return len(self.orfs) > 0
+
+    def __iter__(self):
+        """Iterate over ORF identifiers (mirrors ``dict.__iter__``)."""
+        return iter(orf.orf_id for orf in self.orfs)
+
+    def __getitem__(self, key: str | int) -> OpenReadingFrame:
+        """
+        Access an ORF by identifier string or integer index.
+
+        Examples::
+            orfs['ORF_0']   # by identifier
+            orfs[0]         # by index
+        """
+        if isinstance(key, int):
+            return self.orfs[key]
+        for orf in self.orfs:
+            if orf.orf_id == key:
+                return orf
+        raise KeyError(key)
+
+    def __contains__(self, orf_id: str) -> bool:
+        """True if an ORF with *orf_id* is present in the collection."""
+        return any(orf.orf_id == orf_id for orf in self.orfs)
+
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index e731a6f..82e486b 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -945,12 +945,25 @@ def orf_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, min_length: int
     aln_temp = deepcopy(aln)
     aln_temp.zoom = None
     if non_overlapping_orfs:
-        annotation_dict = aln_temp.get_non_overlapping_conserved_orfs(min_length=min_length)
+        orf_collection = aln_temp.get_non_overlapping_conserved_orfs(min_length=min_length)
     else:
-        annotation_dict = aln_temp.get_conserved_orfs(min_length=min_length)
+        orf_collection = aln_temp.get_conserved_orfs(min_length=min_length)
+
+    # Normalize ORF dataclasses to the annotation dict shape consumed by plotting helpers.
+    annotation_dict = {}
+    for orf_id, orf_data in orf_collection.items():
+        annotation_dict[orf_id] = {
+            'location': orf_data.location,
+            'strand': orf_data.strand,
+            'conservation': orf_data.conservation,
+        }
     # filter dict for zoom
     if aln.zoom is not None:
-        annotation_dict = {key:val for key, val in annotation_dict.items() if max(val['location'][0][0], aln.zoom[0]) <= min(val['location'][0][1], aln.zoom[1])}
+        annotation_dict = {
+            key: val
+            for key, val in annotation_dict.items()
+            if max(val['location'][0][0], aln.zoom[0]) <= min(val['location'][0][1], aln.zoom[1])
+        }
     # add track for plotting
     _add_track_positions(annotation_dict)
     # plot
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 1aecbe4..ad5f689 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -24,7 +24,7 @@
 
 # msaexplorer
 from msaexplorer import config
-from msaexplorer._data_classes import PairwiseDistance, AlignmentStats
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, OpenReadingFrame, OrfContainer
 from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping, _read_alignment
 
 #TODO: Move outputs to dataclasses
@@ -192,8 +192,6 @@ def _create_position_stat_result(self, stat_name: str, values: list | ndarray) -
             stat_name=stat_name,
             positions=positions,
             values=values_array,
-            aln_type=self.aln_type,
-            reference_id=self.reference_id,
         )
 
     def _to_array(self) -> ndarray:
@@ -204,9 +202,6 @@ def _get_reference_seq(self) -> str:
         """get the sequence of the reference sequence or majority consensus"""
         return self.alignment[self.reference_id] if self.reference_id is not None else self.get_consensus()
 
-    def items(self):
-        return self.alignment.items()
-
     def get_reference_coords(self) -> tuple[int, int]:
         """
         Determine the start and end coordinates of the reference sequence
@@ -350,7 +345,7 @@ def get_ambiguous_char(nucleotides: list) -> str:
 
         return consensus
 
-    def get_conserved_orfs(self, min_length: int = 100, identity_cutoff: float | None = None) -> dict:
+    def get_conserved_orfs(self, min_length: int = 100, identity_cutoff: float | None = None) -> OrfContainer:
         """
         **conserved ORF definition:**
             - conserved starts and stops
@@ -459,7 +454,8 @@ def calculate_identity(identity_matrix: ndarray, aln_slice:list) -> float:
         aln_len = self.length
 
         orf_counter = 0
-        orf_dict = {}
+        # use mutable dicts during construction and convert to dataclass at the end for immutability
+        temp_orfs: list[dict] = []
 
         for aln, direction in zip(alignments, ['+', '-']):
             # check for starts and stops in the first seq and then check if these are present in all seqs
@@ -483,28 +479,42 @@ def calculate_identity(identity_matrix: ndarray, aln_slice:list) -> float:
                     # if no stop codon between start and stop --> write to dictionary
                     if not additional_stops(ungapped_sliced_seqs):
                         if direction == '+':
-                            positions = [start, next_stop + 3]
+                            positions = (start, next_stop + 3)
                         else:
-                            positions = [aln_len - next_stop - 3, aln_len - start]
+                            positions = (aln_len - next_stop - 3, aln_len - start)
                         if last_stop != next_stop:
                             last_stop = next_stop
-                            conservation = calculate_identity(identities, positions)
+                            conservation = calculate_identity(identities, list(positions))
                             if identity_cutoff is not None and conservation < identity_cutoff:
                                 continue
-                            orf_dict[f'ORF_{orf_counter}'] = {'location': [positions],
-                                                              'frame': frame,
-                                                              'strand': direction,
-                                                              'conservation': conservation,
-                                                              'internal': []
-                                                              }
+                            temp_orfs.append({
+                                'orf_id': f'ORF_{orf_counter}',
+                                'location': [positions],
+                                'frame': frame,
+                                'strand': direction,
+                                'conservation': conservation,
+                                'internal': [],
+                            })
                             orf_counter += 1
                         else:
-                            if orf_dict:
-                                orf_dict[f'ORF_{orf_counter - 1}']['internal'].append(positions)
-
-        return orf_dict
+                            if temp_orfs:
+                                temp_orfs[-1]['internal'].append(positions)
+
+        # convert mutable intermediate dicts to frozen dataclasses
+        orf_list = [
+            OpenReadingFrame(
+                orf_id=t['orf_id'],
+                location=tuple(t['location']),
+                frame=t['frame'],
+                strand=t['strand'],
+                conservation=t['conservation'],
+                internal=tuple(t['internal']),
+            )
+            for t in temp_orfs
+        ]
+        return OrfContainer(orfs=tuple(orf_list))
 
-    def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cutoff:float = None) -> dict:
+    def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cutoff:float = None) -> OrfContainer:
         """
         First calculates all ORFs and then searches from 5'
         all non-overlapping orfs in the fw strand and from the
@@ -521,37 +531,33 @@ def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cut
 
             frame:    3      2           1      2       1
 
-        :return: dictionary with non-overlapping orfs
+        :return: OrfContainer with non-overlapping orfs
         """
-        orf_dict = self.get_conserved_orfs(min_length, identity_cutoff)
+        all_orfs = self.get_conserved_orfs(min_length, identity_cutoff)
 
         fw_orfs, rw_orfs = [], []
 
-        for orf in orf_dict:
-            if orf_dict[orf]['strand'] == '+':
-                fw_orfs.append((orf, orf_dict[orf]['location'][0]))
+        for orf in all_orfs:
+            orf_obj = all_orfs[orf]
+            if orf_obj.strand == '+':
+                fw_orfs.append((orf, orf_obj.location[0]))
             else:
-                rw_orfs.append((orf, orf_dict[orf]['location'][0]))
+                rw_orfs.append((orf, orf_obj.location[0]))
 
         fw_orfs.sort(key=lambda x: x[1][0])  # sort by start pos
         rw_orfs.sort(key=lambda x: x[1][1], reverse=True)  # sort by stop pos
-        non_overlapping_orfs = []
+        non_overlapping_ids = []
         for orf_list, strand in zip([fw_orfs, rw_orfs], ['+', '-']):
             previous_stop = -1 if strand == '+' else self.length + 1
             for orf in orf_list:
                 if strand == '+' and orf[1][0] >= previous_stop:
-                    non_overlapping_orfs.append(orf[0])
+                    non_overlapping_ids.append(orf[0])
                     previous_stop = orf[1][1]
                 elif strand == '-' and orf[1][1] <= previous_stop:
-                    non_overlapping_orfs.append(orf[0])
+                    non_overlapping_ids.append(orf[0])
                     previous_stop = orf[1][0]
 
-        non_overlap_dict = {}
-        for orf in orf_dict:
-            if orf in non_overlapping_orfs:
-                non_overlap_dict[orf] = orf_dict[orf]
-
-        return non_overlap_dict
+        return OrfContainer(orfs=tuple(all_orfs[orf_id] for orf_id in all_orfs if orf_id in non_overlapping_ids))
 
     def calc_length_stats(self) -> dict:
         """
diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index c85282f..13a2045 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -9,7 +9,7 @@
 import numpy as np
 from numpy import ndarray
 from msaexplorer import config
-from msaexplorer._data_classes import AlignmentStats
+from msaexplorer._data_classes import AlignmentStats, OrfContainer
 from msaexplorer._helpers import _check_and_create_path
 
 
@@ -180,33 +180,31 @@ def stats(stat_data: AlignmentStats | list | ndarray, seperator: str = '\t', pat
         return '\n'.join(lines)
 
 
-def orf(orf_dict: dict, chrom: str, path: str | None = None) -> str | ValueError:
+def orf(orf_dict: OrfContainer, chrom: str, path: str | None = None) -> str | ValueError | None:
     """
-    Exports the ORF dictionary to a .bed file.
+    Exports the ORF collection to a .bed file.
 
-    :param orf_dict: Dictionary containing ORF information.
+    :param orf_dict: OrfContainer instance
     :param chrom: CHROM identifier for bed format.
     :param path: Path to the output .bed file.
-    :param : Reference name
     """
-    if not orf_dict:
-        raise ValueError("The ORF dictionary is empty. Nothing to export.")
-    else:
-        if list(orf_dict[list(orf_dict.keys())[0]].keys()) != ['location', 'frame', 'strand', 'conservation', 'internal']:
-            raise ValueError("The ORF dictionary has not the right format.")
+    if not isinstance(orf_dict, OrfContainer):
+        raise ValueError('The ORF collection must be an instance of msaexplorer._data_classes.OrfContainer.')
 
     _check_and_create_path(path)
 
     lines = []
 
     for orf_id, orf_data in orf_dict.items():
-        lines.append(
-            f"{chrom}\t{orf_data['location'][0][0]}\t{orf_data['location'][0][1]}\t{orf_id}\t{orf_data['conservation']:.2f}\t{orf_data['strand']}"
-        )
+        loc = orf_data.location[0]
+        conservation = orf_data.conservation
+        strand = orf_data.strand
+        lines.append(f"{chrom}\t{loc[0]}\t{loc[1]}\t{orf_id}\t{conservation:.2f}\t{strand}")
 
     if path is not None:
         with open(path, 'w') as out_file:
             out_file.write('\n'.join(lines))
+        return None
     else:
         return '\n'.join(lines)
 
diff --git a/tests/test_export.py b/tests/test_export.py
index ab18877..602b987 100644
--- a/tests/test_export.py
+++ b/tests/test_export.py
@@ -4,7 +4,7 @@
 import pytest
 
 from msaexplorer import export
-from msaexplorer._data_classes import AlignmentStats
+from msaexplorer._data_classes import AlignmentStats, OpenReadingFrame, OrfContainer
 
 
 @pytest.fixture
@@ -117,8 +117,6 @@ def test_alignment_stats_dataclass_output(self):
             stat_name="entropy",
             positions=np.array([2, 4, 6]),
             values=np.array([0.1, 0.2, 0.3]),
-            aln_type="dna",
-            reference_id=None,
         )
 
         result = export.stats(stat_data, seperator=",")
@@ -143,25 +141,26 @@ def test_plain_array_output_uses_zero_based_positions(self):
 
 
 class TestOrfExport:
-    def test_orf_output_is_correct(self):
-        """Test that the ORF output is correct."""
-        orf_dict = {
-            "orf_1": {
-                "location": [(10, 50)],
-                "frame": 1,
-                "strand": "+",
-                "conservation": 87.125,
-                "internal": False,
-            }
-        }
-
-        result = export.orf(orf_dict, chrom="chr1")
+    @pytest.fixture
+    def orf_collection(self):
+        return OrfContainer(orfs=(
+            OpenReadingFrame(
+                orf_id='orf_1',
+                location=((10, 50),),
+                frame=1,
+                strand='+',
+                conservation=87.125,
+            ),
+        ))
+
+    def test_orf_output_is_correct(self, orf_collection):
+        result = export.orf(orf_collection, chrom="chr1")
 
         assert result == "chr1\t10\t50\torf_1\t87.12\t+"
 
     def test_orf_invalid_input_raises(self):
-        """Test that invalid input raises a ValueError."""
-        with pytest.raises(ValueError, match="empty"):
+        with pytest.raises(ValueError, match="instance"):
+            export.orf(OrfContainer(), chrom="chr1")
             export.orf({}, chrom="chr1")
 
 
diff --git a/tests/test_orf_detection.py b/tests/test_orf_detection.py
index 9ca5e0b..b999f32 100644
--- a/tests/test_orf_detection.py
+++ b/tests/test_orf_detection.py
@@ -4,6 +4,7 @@
 
 import pytest
 from msaexplorer.explore import MSA
+from msaexplorer._data_classes import OpenReadingFrame, OrfContainer
 from conftest import create_alignment
 
 
@@ -22,10 +23,10 @@ def test_basic_orf_detection(self):
         orfs = aln.get_conserved_orfs(min_length=9)
 
         assert len(orfs) > 0
-        orf_key = list(orfs.keys())[0]
-        assert orfs[orf_key]['strand'] == '+'
-        assert orfs[orf_key]['location'][0][0] == 0  # Start at position 0
-        assert orfs[orf_key]['location'][0][1] == 12  # End at position 9
+        orf_obj = orfs[0]
+        assert orf_obj.strand == '+'
+        assert orf_obj.location[0][0] == 0  # Start at position 0
+        assert orf_obj.location[0][1] == 12  # End at position 12
 
     def test_with_gaps(self):
         """Test ORF detection handles gaps correctly."""
@@ -40,7 +41,7 @@ def test_with_gaps(self):
         assert len(orfs) == 1
 
     def test_no_orfs(self):
-        """Test that empty dict is returned when no ORFs found."""
+        """Test that empty collection is returned when no ORFs found."""
         alignment_dict = {
             'seq1': 'AAAAAAAAAAA',
             'seq2': 'AAAAAAAAAAA',
@@ -65,7 +66,6 @@ def test_with_internal_stop(self):
 
     def test_multiple_frames(self):
         """Test ORF detection across different reading frames."""
-        # One ORF in frame 1
         alignment_dict = {
             'seq1': 'CATGAAATAAGATGCCCTAG',
             'seq2': 'CATGAAATAAGATGCCCTAG',
@@ -75,12 +75,11 @@ def test_multiple_frames(self):
         orfs = aln.get_conserved_orfs(min_length=9)
 
         # Should detect ORFs in frame 1
-        frames = [orf['frame'] for orf in orfs.values()]
+        frames = [orf.frame for orf in orfs.values()]
         assert 1 in frames
 
     def test_reverse_strand(self):
         """Test ORF detection on reverse complement strand."""
-        # Should find ORF in both strands
         alignment_dict = {
             'seq1': 'ATGTTATTTCATTAA',
             'seq2': 'ATGTTATTTCATTAA',
@@ -88,15 +87,15 @@ def test_reverse_strand(self):
         }
         aln = MSA(create_alignment(alignment_dict))
         orfs = aln.get_conserved_orfs(min_length=9)
-        strands = [orf['strand'] for orf in orfs.values()]
+        strands = [orf.strand for orf in orfs.values()]
 
         assert strands == ['+', '-']
 
     def test_non_conserved_start(self):
         """Test that non-conserved start codons are rejected."""
         alignment_dict = {
-            'seq1': 'ATGAAATAA',  # Has ATG
-            'seq2': 'ATGAAATAA',  # Has ATG
+            'seq1': 'ATGAAATAA',
+            'seq2': 'ATGAAATAA',
             'seq3': 'TTGAAATAA'   # Has TTG (not a start codon)
         }
         aln = MSA(create_alignment(alignment_dict))
@@ -107,8 +106,8 @@ def test_non_conserved_start(self):
     def test_non_conserved_stop(self):
         """Test that non-conserved stop codons are rejected."""
         alignment_dict = {
-            'seq1': 'ATGAAATAA',  # Has TAA stop
-            'seq2': 'ATGAAATAA',  # Has TAA stop
+            'seq1': 'ATGAAATAA',
+            'seq2': 'ATGAAATAA',
             'seq3': 'ATGAAACAA'   # Has CAA (not a stop)
         }
         aln = MSA(create_alignment(alignment_dict))
@@ -144,8 +143,7 @@ def test_internal_orfs_detected(self):
         aln = MSA(create_alignment(alignment_dict))
         orfs = aln.get_conserved_orfs(min_length=9)
 
-        assert len(orfs['ORF_0']['internal']) == 1
-
+        assert len(orfs['ORF_0'].internal) == 1
 
     def test_rna_alignment(self):
         """Test ORF detection works with RNA alignments (U instead of T)."""
@@ -197,16 +195,14 @@ def test_invalid_identity_cutoff(self):
         }
         aln = MSA(create_alignment(alignment_dict))
 
-        # Negative cutoff
         with pytest.raises(ValueError, match='conservation cutoff must be between 0 and 100'):
             aln.get_conserved_orfs(min_length=9, identity_cutoff=-1.0)
 
-        # > 100 cutoff
         with pytest.raises(ValueError, match='conservation cutoff must be between 0 and 100'):
             aln.get_conserved_orfs(min_length=9, identity_cutoff=101.0)
 
     def test_return_structure(self):
-        """Test that returned ORF dictionary has correct structure."""
+        """Test that the return value is an OrfContainer instance with correct structure."""
         alignment_dict = {
             'seq1': 'ATGAAATAA',
             'seq2': 'ATGAAATAA',
@@ -215,24 +211,56 @@ def test_return_structure(self):
         aln = MSA(create_alignment(alignment_dict))
         orfs = aln.get_conserved_orfs(min_length=9)
 
-        assert isinstance(orfs, dict)
+        assert isinstance(orfs, OrfContainer)
+
+        orf_obj = orfs[0]
+        assert isinstance(orf_obj, OpenReadingFrame)
+        assert isinstance(orf_obj.location, tuple)
+        assert isinstance(orf_obj.frame, int)
+        assert orf_obj.strand in ['+', '-']
+        assert isinstance(orf_obj.conservation, float)
+        assert isinstance(orf_obj.internal, tuple)
+
+    def test_getitem_by_name(self):
+        """OrfContainer supports access by orf_id string."""
+        alignment_dict = {
+            'seq1': 'ATGAAATAA',
+            'seq2': 'ATGAAATAA',
+            'seq3': 'ATGAAATAA'
+        }
+        aln = MSA(create_alignment(alignment_dict))
+        orfs = aln.get_conserved_orfs(min_length=9)
+
+        assert orfs['ORF_0'] is orfs[0]
+
+    def test_contains(self):
+        """'ORF_0' should be contained in result, 'ORF_99' should not."""
+        alignment_dict = {
+            'seq1': 'ATGAAATAA',
+            'seq2': 'ATGAAATAA',
+            'seq3': 'ATGAAATAA'
+        }
+        aln = MSA(create_alignment(alignment_dict))
+        orfs = aln.get_conserved_orfs(min_length=9)
 
-        # first orf
-        orf = orfs[list(orfs.keys())[0]]
+        assert 'ORF_0' in orfs
+        assert 'ORF_99' not in orfs
 
-        # Check required keys
-        assert 'location' in orf
-        assert 'frame' in orf
-        assert 'strand' in orf
-        assert 'conservation' in orf
-        assert 'internal' in orf
+    def test_orf_len_and_contains_position(self):
+        """OpenReadingFrame.__len__ and __contains__ work correctly."""
+        alignment_dict = {
+            'seq1': 'ATGAAATAA',
+            'seq2': 'ATGAAATAA',
+            'seq3': 'ATGAAATAA'
+        }
+        aln = MSA(create_alignment(alignment_dict))
+        orfs = aln.get_conserved_orfs(min_length=9)
+        orf_obj = orfs[0]
 
-        # Check types
-        assert isinstance(orf['location'], list)
-        assert isinstance(orf['frame'], int)
-        assert orf['strand'] in ['+', '-']
-        assert isinstance(orf['conservation'], float)
-        assert isinstance(orf['internal'], list)
+        assert len(orf_obj) == 9        # 0..9
+        assert 0 in orf_obj             # start is inside
+        assert 8 in orf_obj             # last position inside
+        assert 9 not in orf_obj         # stop is exclusive
 
 
 class TestGetNonOverlappingConservedOrfs:
@@ -240,7 +268,6 @@ class TestGetNonOverlappingConservedOrfs:
 
     def test_basic(self):
         """Test basic non-overlapping ORF selection."""
-        # Two non-overlapping ORFs (first and second frame) but one internal in the first
         alignment_dict = {
             'seq1': 'ATGAAAATGAAATAACCCTATGGGGTAG',
             'seq2': 'ATGAAAATGAAATAACCCTATGGGGTAG',
@@ -254,7 +281,7 @@ def test_basic(self):
             assert len(orfs) == 2
 
     def test_adjecent_orfs(self):
-        """Test if ORFs that are directly adjacent (touching boundaries) are retained as non-overlapping for both frames."""
+        """Test if ORFs that are directly adjacent are retained as non-overlapping."""
         alignment_dict = {
             'seq1': 'ATGAAAATGTAAATGAAAATGTAA',
             'seq2': 'ATGAAAATGTAAATGAAAATGTAA',
@@ -266,25 +293,6 @@ def test_adjecent_orfs(self):
             orfs = aln.get_non_overlapping_conserved_orfs(min_length=9)
             assert len(orfs) == 2
 
-    def test_preserves_structure(self):
-        """Test that non-overlapping ORFs preserve the same data structure."""
-        alignment_dict = {
-            'seq1': 'ATGAAATAA',
-            'seq2': 'ATGAAATAA',
-            'seq3': 'ATGAAATAA'
-        }
-        aln = MSA(create_alignment(alignment_dict))
-        orfs = aln.get_non_overlapping_conserved_orfs(min_length=9)
-        # first orf
-        orf = orfs[list(orfs.keys())[0]]
-
-        # Check all required keys are present
-        assert 'location' in orf
-        assert 'frame' in orf
-        assert 'strand' in orf
-        assert 'conservation' in orf
-        assert 'internal' in orf
-
     def test_with_identity_cutoff(self):
         """Test non-overlapping ORFs with identity cutoff."""
         alignment_dict = {
@@ -295,6 +303,5 @@ def test_with_identity_cutoff(self):
         aln = MSA(create_alignment(alignment_dict))
         orfs = aln.get_non_overlapping_conserved_orfs(min_length=9, identity_cutoff=90.0)
 
-        # All returned ORFs should meet identity cutoff
         for orf in orfs.values():
-            assert orf['conservation'] >= 90.0
+            assert orf.conservation >= 90.0

From 8b9dc151ea2cddc4ab02fc8e31a908e1565e7611 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 20 Mar 2026 01:54:18 +0100
Subject: [PATCH 12/25] added Variant class

---
 msaexplorer/_data_classes.py     |  35 ++++++++-
 msaexplorer/draw.py              |  51 +++++++------
 msaexplorer/explore.py           |  37 ++++++----
 msaexplorer/export.py            | 122 ++++++++++++-------------------
 tests/test_export.py             |  62 +++++++---------
 tests/test_orf_detection.py      |   4 +-
 tests/test_stats_calculations.py |  37 +++++-----
 7 files changed, 174 insertions(+), 174 deletions(-)

diff --git a/msaexplorer/_data_classes.py b/msaexplorer/_data_classes.py
index 847ad9b..8d8f4cd 100644
--- a/msaexplorer/_data_classes.py
+++ b/msaexplorer/_data_classes.py
@@ -2,13 +2,44 @@
 this contains the dataclasses used to store the data for the msa explorer. these are not meant to be used outside of this package.
 """
 
-# build-in
+# built-in
 from dataclasses import dataclass, field
 
 # libs
 from numpy import ndarray
 
 
+@dataclass(frozen=True)
+class SingleNucleotidePolymorphism:
+    """
+    SNP data for one alignment position.
+    """
+
+    ref: str
+    alt: dict[str, tuple[float, tuple[str, ...]]] = field(default_factory=dict)
+
+
+@dataclass(frozen=True)
+class VariantCollection:
+    """Container for SNPs"""
+
+    chrom: str
+    positions: dict[int, SingleNucleotidePolymorphism] = field(default_factory=dict)
+
+    def __post_init__(self):
+        """sort the positions by position (key)"""
+        object.__setattr__(self, 'positions', dict(sorted(self.positions.items())))
+
+    def __len__(self) -> int:
+        return len(self.positions)
+
+    def __iter__(self):
+        return iter(self.positions)
+
+    def __contains__(self, position: int) -> bool:
+        return position in self.positions
+
+
 @dataclass(frozen=True)
 class AlignmentStats:
     """
@@ -105,7 +136,7 @@ def __contains__(self, position: int) -> bool:
 
 
 @dataclass(frozen=True)
-class OrfContainer:
+class OrfCollection:
     """
     Ordered collection of `OpenReadingFrame` objects returned by
     `MSA.get_conserved_orfs` or`MSA.get_non_overlapping_conserved_orfs`.
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 82e486b..49899f9 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -778,32 +778,31 @@ def variant_plot(aln: explore.MSA | str, ax: plt.Axes | None = None, lollisize:
     # define where to plot (each ref type gets a separate line)
     ref_y_positions, y_pos, detected_var = {}, 0, set()
 
-    # iterate over snp dict
-    for pos in snps['POS']:
-        for identifier in snps['POS'][pos]:
-            # fill in y pos dict
-            if identifier == 'ref':
-                if snps['POS'][pos]['ref'] not in ref_y_positions:
-                    ref_y_positions[snps['POS'][pos]['ref']] = y_pos
-                    y_pos += 1.1
-                    continue
-            # plot
-            if identifier == 'ALT':
-                for alt in snps['POS'][pos]['ALT']:
-                    ax.vlines(x=pos + aln.zoom[0] if aln.zoom is not None else pos,
-                              ymin=ref_y_positions[snps['POS'][pos]['ref']],
-                              ymax=ref_y_positions[snps['POS'][pos]['ref']] + snps['POS'][pos]['ALT'][alt]['AF'],
-                              color=colors[alt],
-                              zorder=100,
-                              linewidth=lollisize[0]
-                              )
-                    ax.plot(pos + aln.zoom[0] if aln.zoom is not None else pos,
-                            ref_y_positions[snps['POS'][pos]['ref']] + snps['POS'][pos]['ALT'][alt]['AF'],
-                            color=colors[alt],
-                            marker='o',
-                            markersize=lollisize[1]
-                            )
-                    detected_var.add(alt)
+    # iterate over SNPs
+    for pos, snp_pos in snps.positions.items():
+        if snp_pos.ref not in ref_y_positions:
+            ref_y_positions[snp_pos.ref] = y_pos
+            y_pos += 1.1
+
+        for alt, (af, _) in snp_pos.alt.items():
+            x_pos = pos + aln.zoom[0] if aln.zoom is not None else pos
+            y_ref = ref_y_positions[snp_pos.ref]
+            ax.vlines(
+                x=x_pos,
+                ymin=y_ref,
+                ymax=y_ref + af,
+                color=colors[alt],
+                zorder=100,
+                linewidth=lollisize[0]
+            )
+            ax.plot(
+                x_pos,
+                y_ref + af,
+                color=colors[alt],
+                marker='o',
+                markersize=lollisize[1]
+            )
+            detected_var.add(alt)
 
     # plot hlines
     for y_char in ref_y_positions:
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index ad5f689..aaf4fe0 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -24,7 +24,7 @@
 
 # msaexplorer
 from msaexplorer import config
-from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, OpenReadingFrame, OrfContainer
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, OpenReadingFrame, OrfCollection, SingleNucleotidePolymorphism, VariantCollection
 from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping, _read_alignment
 
 #TODO: Move outputs to dataclasses
@@ -345,7 +345,7 @@ def get_ambiguous_char(nucleotides: list) -> str:
 
         return consensus
 
-    def get_conserved_orfs(self, min_length: int = 100, identity_cutoff: float | None = None) -> OrfContainer:
+    def get_conserved_orfs(self, min_length: int = 100, identity_cutoff: float | None = None) -> OrfCollection:
         """
         **conserved ORF definition:**
             - conserved starts and stops
@@ -512,9 +512,9 @@ def calculate_identity(identity_matrix: ndarray, aln_slice:list) -> float:
             )
             for t in temp_orfs
         ]
-        return OrfContainer(orfs=tuple(orf_list))
+        return OrfCollection(orfs=tuple(orf_list))
 
-    def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cutoff:float = None) -> OrfContainer:
+    def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cutoff:float = None) -> OrfCollection:
         """
         First calculates all ORFs and then searches from 5'
         all non-overlapping orfs in the fw strand and from the
@@ -557,7 +557,7 @@ def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cut
                     non_overlapping_ids.append(orf[0])
                     previous_stop = orf[1][0]
 
-        return OrfContainer(orfs=tuple(all_orfs[orf_id] for orf_id in all_orfs if orf_id in non_overlapping_ids))
+        return OrfCollection(orfs=tuple(all_orfs[orf_id] for orf_id in all_orfs if orf_id in non_overlapping_ids))
 
     def calc_length_stats(self) -> dict:
         """
@@ -1137,7 +1137,7 @@ def calc_pairwise_distance_to_reference(self, distance_type:str='ghd') -> Pairwi
             distances=np.array(distances)
         )
 
-    def get_snps(self, include_ambig:bool=False) -> dict:
+    def get_snps(self, include_ambig:bool=False) -> VariantCollection:
         """
         Calculate snps similar to snp-sites (output is comparable):
         https://github.com/sanger-pathogens/snp-sites
@@ -1145,7 +1145,7 @@ def get_snps(self, include_ambig:bool=False) -> dict:
         The SNPs are compared to a majority consensus sequence or to a reference if it has been set.
 
         :param include_ambig: Include ambiguous snps (default: False)
-        :return: dictionary containing snp positions and their variants including their frequency.
+        :return: dataclass containing SNP positions and their variants including frequency.
         """
         aln = self.alignment
         ref = self._get_reference_seq()
@@ -1172,7 +1172,7 @@ def get_snps(self, include_ambig:bool=False) -> dict:
                 snps = [x for x in snps if alt_chars[x] not in config.AMBIG_CHARS[self.aln_type]]
                 if not snps:
                     continue
-            if pos not in snp_dict:
+            if pos not in snp_dict['POS']:
                 snp_dict['POS'][pos] = {'ref': reference_char, 'ALT': {}}
             for snp in snps:
                 if alt_chars[snp] not in snp_dict['POS'][pos]['ALT']:
@@ -1188,7 +1188,15 @@ def get_snps(self, include_ambig:bool=False) -> dict:
                 for alt in snp_dict['POS'][pos]['ALT']:
                     snp_dict['POS'][pos]['ALT'][alt]['AF'] /= len(aln)
 
-        return snp_dict
+        snp_positions = {}
+        for pos, pos_info in snp_dict['POS'].items():
+            alt = {
+                allele: (details['AF'], tuple(details['SEQ_ID']))
+                for allele, details in pos_info['ALT'].items()
+            }
+            snp_positions[pos] = SingleNucleotidePolymorphism(ref=pos_info['ref'], alt=alt)
+
+        return VariantCollection(chrom=snp_dict['#CHROM'], positions=snp_positions)
 
     def calc_transition_transversion_score(self) -> AlignmentStats:
         """
@@ -1205,14 +1213,13 @@ def calc_transition_transversion_score(self) -> AlignmentStats:
         snps = self.get_snps()
         score = [0]*self.length
 
-        for pos in snps['POS']:
-            t_score_temp = 0
-            for alt in snps['POS'][pos]['ALT']:
+        for pos, snp in snps.positions.items():
+            for alt, (af, _) in snp.alt.items():
                 # check the type of substitution
-                if snps['POS'][pos]['ref'] + alt in ['AG', 'GA', 'CT', 'TC', 'CU', 'UC']:
-                    score[pos] += snps['POS'][pos]['ALT'][alt]['AF']
+                if snp.ref + alt in ['AG', 'GA', 'CT', 'TC', 'CU', 'UC']:
+                    score[pos] += af
                 else:
-                    score[pos] -= snps['POS'][pos]['ALT'][alt]['AF']
+                    score[pos] -= af
 
         return self._create_position_stat_result('ts tv score', score)
 
diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index 13a2045..646ab4f 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -9,111 +9,79 @@
 import numpy as np
 from numpy import ndarray
 from msaexplorer import config
-from msaexplorer._data_classes import AlignmentStats, OrfContainer
+from msaexplorer._data_classes import AlignmentStats, OrfCollection, VariantCollection
 from msaexplorer._helpers import _check_and_create_path
 
 
-def snps(snp_dict: dict, format_type: str = 'vcf', path: str | None = None) -> str | None | ValueError:
+def snps(snp_data: VariantCollection, format_type: str = 'vcf', path: str | None = None) -> str | None | ValueError:
     """
-    Export a SNP dictionary to a VCF or tabular format. Importantly, the input dictionary has to be in the standard
-    format that MSAexplorer produces.
+    Export SNP data from a VariantCollection to VCF or tabular format.
 
-    :param snp_dict: Dictionary containing SNP positions and variant information.
+    :param snp_data: VariantCollection containing SNP positions and variant information.
     :param format_type: Format type ('vcf' or 'tabular'). Default is 'vcf'.
     :param path: Path to output VCF or tabular format. (optional)
     :return: A string containing the SNP data in the requested format.
-    :raises ValueError: if the input dictionary is missing required keys or format_type is invalid.
+    :raises ValueError: if the input type is invalid or format_type is invalid.
     """
 
     def _validate():
-        if not isinstance(snp_dict, dict):
-            raise ValueError('Input SNP data must be a dictionary.')
-        for key in ['#CHROM', 'POS']:
-            if key not in snp_dict:
-                raise ValueError(f"Missing required key '{key}' in SNP data.")
-        if not isinstance(snp_dict['POS'], dict):
-            raise ValueError('Expected the \'POS\' key to contain a dictionary of positions.')
+        if not isinstance(snp_data, VariantCollection):
+            raise ValueError('Input SNP data must be a VariantCollection dataclass.')
         if format_type not in ['vcf', 'tabular']:
             raise ValueError('Invalid format_type.')
         _check_and_create_path(path)
 
-    def _vcf_format(snp_dict: dict) -> list:
-        """
-        Produce  vcf formatted SNP data.
-        :param snp_dict: dictionary containing SNP positions and variant information.
-        :return: list of lines to write
-        """
-        output_lines = []
-        # VCF header
-        output_lines.append('##fileformat=VCFv4.2')
-        output_lines.append('##source=MSAexplorer')
-        output_lines.append('#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO')
-        # process each SNP position in sorted order
-        for pos in sorted(snp_dict['POS'].keys()):
-            pos_info = snp_dict['POS'][pos]
-            ref = pos_info.get('ref', '.')
-            alt_dict = pos_info.get('ALT', {})
-            # Create comma-separated list of alternative alleles
-            alt_alleles = ",".join(alt_dict.keys()) if alt_dict else "."
-            # Prepare INFO field: include allele frequencies and sequence IDs
-            afs = []
-            seq_ids = []
-            for alt, details in alt_dict.items():
-                af = details.get('AF', 0)
-                afs.append(str(af))
-                seq_ids.append("|".join(details.get('SEQ_ID', [])))
+    def _vcf_format(data: VariantCollection) -> list[str]:
+        """Produce VCF formatted SNP data."""
+        output_lines = [
+            '##fileformat=VCFv4.2',
+            '##source=MSAexplorer',
+            '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO',
+        ]
+
+        for pos in sorted(data.positions.keys()):
+            pos_info = data.positions[pos]
+            alt_dict = pos_info.alt
+            alt_alleles = ','.join(alt_dict.keys()) if alt_dict else '.'
+
+            afs = [str(af) for af, _seq_ids in alt_dict.values()]
+            seq_ids = ['|'.join(seq_ids) for _af, seq_ids in alt_dict.values()]
             info_fields = []
             if afs:
-                info_fields.append("AF=" + ",".join(afs))
+                info_fields.append('AF=' + ','.join(afs))
             if seq_ids:
-                info_fields.append("SEQ_ID=" + ",".join(seq_ids))
-            info = ";".join(info_fields) if info_fields else "."
+                info_fields.append('SEQ_ID=' + ','.join(seq_ids))
+            info = ';'.join(info_fields) if info_fields else '.'
 
-            # VCF is 1-indexed; we assume pos is 0-indexed and add 1
-            line = f"{snp_dict['#CHROM']}\t{pos + 1}\t.\t{ref}\t{alt_alleles}\t.\t.\t{info}"
-            output_lines.append(line)
+            output_lines.append(
+                f'{data.chrom}\t{pos + 1}\t.\t{pos_info.ref}\t{alt_alleles}\t.\t.\t{info}'
+            )
 
         return output_lines
 
-    def _tabular_format(snp_dict: dict) -> list:
-        """
-        Produce  tabular formatted SNP data.
-
-        :param snp_dict: dictionary containing SNP positions and variant information.
-        :return: list of lines to write
-        """
-        output_lines = []
-        # Create a header for the tabular output
-        output_lines.append('CHROM\tPOS\tREF\tALT\tAF\tSEQ_ID')
-
-        # Process each SNP position and each alternative allele
-        for pos in sorted(snp_dict['POS'].keys()):
-            pos_info = snp_dict['POS'][pos]
-            ref = pos_info.get('ref', '.')
-            alt_dict = pos_info.get('ALT', {})
-            for alt, details in alt_dict.items():
-                af = details.get('AF', 0)
-                seq_id = ",".join(details.get('SEQ_ID', []))
-                output_lines.append(f"{snp_dict['#CHROM']}\t{pos + 1}\t{ref}\t{alt}\t{af}\t{seq_id}")
+    def _tabular_format(data: VariantCollection) -> list[str]:
+        """Produce tabular formatted SNP data."""
+        output_lines = ['CHROM\tPOS\tREF\tALT\tAF\tSEQ_ID']
+
+        for pos in sorted(data.positions.keys()):
+            pos_info = data.positions[pos]
+            for alt, (af, seq_ids) in pos_info.alt.items():
+                output_lines.append(
+                    f'{data.chrom}\t{pos + 1}\t{pos_info.ref}\t{alt}\t{af}\t{",".join(seq_ids)}'
+                )
 
         return output_lines
 
-    # validate correct input format
     _validate()
+    lines = _vcf_format(snp_data) if format_type == 'vcf' else _tabular_format(snp_data)
 
-    # generate line data
-    if format_type == 'vcf':
-        lines = _vcf_format(snp_dict)
-    else:
-        lines = _tabular_format(snp_dict)
-
-    # export to file or return plain text
     if path is not None:
-        out_path = f"{path}.{format_type}"
+        out_path = f'{path}.{format_type}'
         with open(out_path, 'w') as out_file:
             out_file.write('\n'.join(lines))
-    else:
-        return '\n'.join(lines)
+        return None
+
+    return '\n'.join(lines)
 
 
 def fasta(sequence: str | dict, header: str | None = None, path: str | None = None) -> str | None:
@@ -180,7 +148,7 @@ def stats(stat_data: AlignmentStats | list | ndarray, seperator: str = '\t', pat
         return '\n'.join(lines)
 
 
-def orf(orf_dict: OrfContainer, chrom: str, path: str | None = None) -> str | ValueError | None:
+def orf(orf_dict: OrfCollection, chrom: str, path: str | None = None) -> str | ValueError | None:
     """
     Exports the ORF collection to a .bed file.
 
@@ -188,7 +156,7 @@ def orf(orf_dict: OrfContainer, chrom: str, path: str | None = None) -> str | Va
     :param chrom: CHROM identifier for bed format.
     :param path: Path to the output .bed file.
     """
-    if not isinstance(orf_dict, OrfContainer):
+    if not isinstance(orf_dict, OrfCollection):
         raise ValueError('The ORF collection must be an instance of msaexplorer._data_classes.OrfContainer.')
 
     _check_and_create_path(path)
diff --git a/tests/test_export.py b/tests/test_export.py
index 602b987..854f9dd 100644
--- a/tests/test_export.py
+++ b/tests/test_export.py
@@ -4,36 +4,34 @@
 import pytest
 
 from msaexplorer import export
-from msaexplorer._data_classes import AlignmentStats, OpenReadingFrame, OrfContainer
+from msaexplorer._data_classes import AlignmentStats, OpenReadingFrame, OrfCollection, SingleNucleotidePolymorphism, VariantCollection
 
 
 @pytest.fixture
-def snp_dict():
-    """Returns a SNP dictionary for testing."""
-    return {
-        "#CHROM": "ref",
-        "POS": {
-            3: {
-                "ref": "C",
-                "ALT": {
-                    "T": {"AF": 1.0, "SEQ_ID": ["s1"]},
-                },
-            },
-            1: {
-                "ref": "A",
-                "ALT": {
-                    "G": {"AF": 0.5, "SEQ_ID": ["s1", "s2"]},
-                    "-": {"AF": 0.5, "SEQ_ID": ["s3"]},
+def snp_result():
+    """Returns a SNP dataclass object for testing."""
+    return VariantCollection(
+        chrom="ref",
+        positions={
+            3: SingleNucleotidePolymorphism(
+                ref="C",
+                alt={"T": (1.0, ("s1",))},
+            ),
+            1: SingleNucleotidePolymorphism(
+                ref="A",
+                alt={
+                    "G": (0.5, ("s1", "s2")),
+                    "-": (0.5, ("s3",)),
                 },
-            },
+            ),
         },
-    }
+    )
 
 
 class TestSnpsExport:
-    def test_vcf_output_is_correct_and_sorted(self, snp_dict):
+    def test_vcf_output_is_correct_and_sorted(self, snp_result):
         """Test that the VCF has correct vcf format."""
-        result = export.snps(snp_dict, format_type="vcf")
+        result = export.snps(snp_result, format_type="vcf")
 
         assert result == "\n".join(
             [
@@ -45,9 +43,9 @@ def test_vcf_output_is_correct_and_sorted(self, snp_dict):
             ]
         )
 
-    def test_tabular_output_is_correct(self, snp_dict):
+    def test_tabular_output_is_correct(self, snp_result):
         """Test that the VCF has correct tabular format."""
-        result = export.snps(snp_dict, format_type="tabular")
+        result = export.snps(snp_result, format_type="tabular")
 
         assert result == "\n".join(
             [
@@ -60,20 +58,18 @@ def test_tabular_output_is_correct(self, snp_dict):
 
     def test_invalid_input_raises_value_error(self):
         """Test that invalid input raises a ValueError."""
-        with pytest.raises(ValueError, match="must be a dictionary"):
+        with pytest.raises(ValueError, match="must be a VariantCollection"):
             export.snps([], format_type="vcf")
 
-        with pytest.raises(ValueError, match="Missing required key"):
-            export.snps({"POS": {}}, format_type="vcf")
-
+    def test_invalid_format_raises_value_error(self, snp_result):
         with pytest.raises(ValueError, match="Invalid format_type"):
-            export.snps({"#CHROM": "ref", "POS": {}}, format_type="csv")
+            export.snps(snp_result, format_type="csv")
 
-    def test_file_export_creates_expected_extension(self, snp_dict, tmp_path):
+    def test_file_export_creates_expected_extension(self, snp_result, tmp_path):
         """Test that the file extension is correct when written."""
         path_without_ext = tmp_path / "nested" / "snps_output"
 
-        result = export.snps(snp_dict, format_type="vcf", path=str(path_without_ext))
+        result = export.snps(snp_result, format_type="vcf", path=str(path_without_ext))
 
         assert result is None
         written = path_without_ext.with_suffix(path_without_ext.suffix + ".vcf")
@@ -143,7 +139,7 @@ def test_plain_array_output_uses_zero_based_positions(self):
 class TestOrfExport:
     @pytest.fixture
     def orf_collection(self):
-        return OrfContainer(orfs=(
+        return OrfCollection(orfs=(
             OpenReadingFrame(
                 orf_id='orf_1',
                 location=((10, 50),),
@@ -160,7 +156,7 @@ def test_orf_output_is_correct(self, orf_collection):
 
     def test_orf_invalid_input_raises(self):
         with pytest.raises(ValueError, match="instance"):
-            export.orf(OrfContainer(), chrom="chr1")
+            export.orf(OrfCollection(), chrom="chr1")
             export.orf({}, chrom="chr1")
 
 
@@ -209,5 +205,3 @@ def test_percent_recovery_value_must_be_float(self):
         """Test that the percent recovery value must be a float."""
         with pytest.raises(ValueError, match="invalid"):
             export.percent_recovery({"seq1": 75})
-
-
diff --git a/tests/test_orf_detection.py b/tests/test_orf_detection.py
index b999f32..6d87ac4 100644
--- a/tests/test_orf_detection.py
+++ b/tests/test_orf_detection.py
@@ -4,7 +4,7 @@
 
 import pytest
 from msaexplorer.explore import MSA
-from msaexplorer._data_classes import OpenReadingFrame, OrfContainer
+from msaexplorer._data_classes import OpenReadingFrame, OrfCollection
 from conftest import create_alignment
 
 
@@ -211,7 +211,7 @@ def test_return_structure(self):
         aln = MSA(create_alignment(alignment_dict))
         orfs = aln.get_conserved_orfs(min_length=9)
 
-        assert isinstance(orfs, OrfContainer)
+        assert isinstance(orfs, OrfCollection)
 
         orf_obj = orfs[0]
         assert isinstance(orf_obj, OpenReadingFrame)
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index beac784..75eaea4 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -4,7 +4,7 @@
 import numpy as np
 from conftest import create_alignment
 from msaexplorer.explore import MSA
-from msaexplorer._data_classes import PairwiseDistance, AlignmentStats
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, VariantCollection
 
 
 class TestCalcEntropy:
@@ -167,8 +167,9 @@ def test_no_variants(self):
         msa = MSA(create_alignment({"s1": "ACGT", "s2": "ACGT"}))
         snps = msa.get_snps()
 
-        assert snps["#CHROM"] == "consensus"
-        assert snps["POS"] == {}
+        assert isinstance(snps, VariantCollection)
+        assert snps.chrom == "consensus"
+        assert snps.positions == {}
 
     def test_exclude_ambiguous(self):
         """Test if no variants are present for include_ambig=False"""
@@ -182,10 +183,10 @@ def test_exclude_ambiguous(self):
         )
         snps = msa.get_snps(include_ambig=False)
 
-        assert snps["POS"] == {}
+        assert snps.positions == {}
 
     def test_include_ambiguous(self):
-        """Test if variants are present in correct frequency for include_ambig=False"""
+        """Test if variants are present in correct frequency for include_ambig=True"""
         msa = MSA(
             create_alignment({
                 "ref": "AAAA",
@@ -196,12 +197,12 @@ def test_include_ambiguous(self):
         )
 
         snps = msa.get_snps(include_ambig=True)
-        alts = snps["POS"][1]["ALT"]
+        alts = snps.positions[1].alt
 
-        assert 1 in snps["POS"]
+        assert 1 in snps
         assert set(alts.keys()) == {"C", "N"}
-        assert alts["C"]["AF"] == 0.5
-        assert alts["N"]["AF"] == 0.5
+        assert alts["C"][0] == 0.5
+        assert alts["N"][0] == 0.5
 
     def test_alt_gaps_exclude_ambiguous(self):
         """Test that no gaps are present for include_ambig=False"""
@@ -216,7 +217,7 @@ def test_alt_gaps_exclude_ambiguous(self):
 
         snps = msa.get_snps(include_ambig=False)
 
-        assert snps["POS"] == {}
+        assert snps.positions == {}
 
     def test_alt_gaps_include_ambiguous(self):
         """Test that there are gaps present for include_ambig=True"""
@@ -230,11 +231,11 @@ def test_alt_gaps_include_ambiguous(self):
         )
 
         snps = msa.get_snps(include_ambig=True)
-        alts = snps["POS"][1]["ALT"]
+        alts = snps.positions[1].alt
 
-        assert 1 in snps["POS"]
+        assert 1 in snps
         assert set(alts.keys()) == {'-', 'C'}
-        assert alts["-"]["AF"] == 0.5
+        assert alts["-"][0] == 0.5
 
     def test_gaps_in_reference(self):
         """Test that gaps are always included if in the reference"""
@@ -249,8 +250,8 @@ def test_gaps_in_reference(self):
 
         snps = msa.get_snps(include_ambig=False)
 
-        assert snps["POS"][1]['ref'] == '-'
-        assert snps["POS"][2]['ref'] == '-'
+        assert snps.positions[1].ref == '-'
+        assert snps.positions[2].ref == '-'
 
     def test_correct_sequence_identifiers(self):
         msa = MSA(
@@ -263,9 +264,9 @@ def test_correct_sequence_identifiers(self):
         )
         snps = msa.get_snps()
 
-        assert snps["#CHROM"] == "ref"
-        assert set(snps["POS"][1]["ALT"]["C"]["SEQ_ID"]) == {"q1"}
-        assert set(snps["POS"][1]["ALT"]["G"]["SEQ_ID"]) == {"q2"}
+        assert snps.chrom == "ref"
+        assert set(snps.positions[1].alt["C"][1]) == {"q1"}
+        assert set(snps.positions[1].alt["G"][1]) == {"q2"}
 
 
 class TestCalcPairwiseIdentityMatrix:

From 56445a3e5fde24eda023f557458fde9c7b9a5a0c Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 20 Mar 2026 12:08:35 +0100
Subject: [PATCH 13/25] added project specific copilot instructions

---
 .github/copilot-instructions.md | 123 ++++++++++++++++++++++++++++++++
 1 file changed, 123 insertions(+)
 create mode 100644 .github/copilot-instructions.md

diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md
new file mode 100644
index 0000000..d1e0af6
--- /dev/null
+++ b/.github/copilot-instructions.md
@@ -0,0 +1,123 @@
+# Copilot Instructions for MSAexplorer
+
+This file defines repository-specific standards for AI-assisted edits.
+Use these rules for all code, tests, and documentation changes.
+
+## 1) Core Principles
+
+- Keep changes small, readable, and reviewable.
+- Prefer simple and explicit code over clever abstractions.
+- Avoid overengineering: introduce new classes/indirection only when they clearly improve maintainability.
+- Keep behavior stable unless the task explicitly requests a breaking change.
+- Preserve existing architecture and naming patterns.
+
+## 2) Repository Structure and Interaction
+
+This repository contains two Python packages with different roles:
+
+- `msaexplorer/` (core library)
+  - Analysis and data model (`explore.py`, `_data_classes.py`)
+  - Plotting (`draw.py`)
+  - Export helpers (`export.py`)
+  - CLI entrypoint (`cli.py`)
+- `app_src/` (Shiny app frontend)
+  - UI/server glue for interactive exploration (`shiny_user_interface.py`, `shiny_server.py`, `shiny_plots.py`)
+
+How they interact:
+
+- `app_src` should consume stable APIs from `msaexplorer`.
+- `msaexplorer` must remain independently usable without the app.
+- Feature work starts in `msaexplorer`; `app_src` is adapted afterwards if needed.
+- Keep app-only concerns out of core modules.
+
+## 3) Code Cleanliness Standards
+
+- Use single quotes for normal strings: `'text'`.
+- Use double quotes only when required (e.g., quote escaping readability) or for docstrings.
+- Keep functions focused on one task; split long functions when logic becomes hard to scan.
+- Prefer descriptive names over abbreviations.
+- Remove dead code, unused imports, and outdated comments.
+- Add comments only for non-obvious reasoning, not for trivial operations.
+- Keep public APIs typed (input and return types) where practical.
+- Reuse existing helpers/dataclasses before adding new structures.
+
+## 4) Naming Conventions
+
+- Functions/methods: `snake_case`
+- Variables: `snake_case`
+- Classes/dataclasses: `PascalCase`
+- Constants: `UPPER_SNAKE_CASE`
+- Internal helpers: prefix with `_` (module-private intent)
+- Boolean names should read naturally (`include_ambig`, `show_legend`, `is_valid`)
+
+## 5) Docstring Standards
+
+Use concise, practical docstrings.
+
+- Use triple double quotes for docstrings.
+- First line: short summary sentence in imperative/present style.
+- Describe parameters, return values, and raised exceptions when relevant.
+- Keep docstrings short and precise; avoid tutorial-length blocks in function docstrings.
+- Document every public class and public function.
+- For internal/private helpers, add docstrings when behavior is not obvious.
+
+Recommended structure:
+
+```python
+"""
+Calculate pairwise distances against the current reference.
+
+:param distance_type: Distance metric key.
+:return: Pairwise distances container.
+:raises ValueError: If the metric is unsupported.
+"""
+```
+
+## 6) Documentation Guidelines (pdoc-oriented)
+
+- API docs are generated primarily from docstrings (pdoc).
+- Keep module/class/function docstrings accurate and synchronized with behavior.
+- Prefer small executable examples over long narrative snippets.
+- Any example added to docs (including package `__init__.py` examples) must be tested locally before integration.
+- Do not document planned behavior as if already implemented.
+- When refactoring signatures or return types, update related docstrings in the same change.
+
+## 7) Testing Requirements
+
+Every code change should be validated before merge.
+
+Minimum expectations:
+
+- Run targeted tests for edited modules.
+- Add/adjust tests when behavior, output format, or API contracts change.
+- Ensure export/plot/stat outputs are validated with deterministic assertions.
+
+Typical commands:
+
+```bash
+pytest -q
+```
+
+For local source verification during development:
+
+```bash
+PYTHONPATH="/absolute/path/to/MSAexplorer" pytest -q
+```
+
+Guidance:
+
+- Prefer explicit regression tests over broad smoke-only checks.
+- If a bug is fixed, include a test that fails before and passes after.
+- Keep fixtures minimal and representative.
+
+## 8) Change Discipline
+
+Before submitting changes:
+
+- Confirm naming and style consistency.
+- Confirm docstrings match implementation.
+- Confirm tests pass.
+- Confirm no unrelated edits were introduced.
+
+If uncertain, prefer the simpler design.
+

From 98b76f2690b47773b6788bb28b6064dc97956888 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 20 Mar 2026 12:10:26 +0100
Subject: [PATCH 14/25] updated version

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 88eba53..542934c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = 'msaexplorer'
-version = '0.4.6'
+version = '0.5'
 description = 'A tool for exploring multiple sequence alignments.'
 keywords = ['msa, multiple sequence alignments, plotting, shiny app, alignment analysis']
 dependencies = ['numpy>=2.0', 'matplotlib>=3.8', 'biopython>=1.81']

From f374400fda1806eb8993b6bebf9a3dae74b1d0ca Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 20 Mar 2026 12:40:09 +0100
Subject: [PATCH 15/25] added test suite for drawing and populating axis for
 the different plots

---
 msaexplorer/draw.py      |   2 +-
 tests/test_draw_plots.py | 191 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 192 insertions(+), 1 deletion(-)
 create mode 100644 tests/test_draw_plots.py

diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 49899f9..2ed3f3f 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -1283,7 +1283,7 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
         # slice alignment and replace the original alignment
         aln_tmp._alignment = {
             seq_id: seq[left_side:right_side]
-            for seq_id, seq in aln.items()
+            for seq_id, seq in aln.alignment.items()
         }
         window_result = aln_tmp.calc_pairwise_distance_to_reference(distance_type=distance_calculation)
         value_map = {seq_id: value for seq_id, value in zip(window_result.sequence_ids, window_result.distances)}
diff --git a/tests/test_draw_plots.py b/tests/test_draw_plots.py
new file mode 100644
index 0000000..568506c
--- /dev/null
+++ b/tests/test_draw_plots.py
@@ -0,0 +1,191 @@
+"""Tests for the main plotting helpers in ``msaexplorer.draw``."""
+
+from pathlib import Path
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import pytest
+from matplotlib.axes import Axes
+
+from conftest import create_alignment
+from msaexplorer import draw
+from msaexplorer.explore import MSA
+
+
+@pytest.fixture
+def ax() -> Axes:
+    """Create a fresh matplotlib axis for each test."""
+    fig, axis = plt.subplots()
+    yield axis
+    plt.close(fig)
+
+
+@pytest.fixture
+def dna_msa() -> MSA:
+    """Create a small nucleotide alignment with mismatches, a gap and a mask character."""
+    return MSA(
+        create_alignment(
+            {
+                'ref': 'ATGAAATTTTAA',
+                'seq1': 'ATGAAATTTTAA',
+                'seq2': 'ATGACATTTTAA',
+                'seq3': 'ATGAAAT-TNAA',
+            }
+        ),
+        reference_id='ref',
+    )
+
+
+@pytest.fixture
+def aa_msa() -> MSA:
+    """Create a small amino-acid alignment for AA plotting tests."""
+    return MSA(
+        create_alignment(
+            {
+                'ref_aa': 'MKTWQALV',
+                'seq1_aa': 'MKTWQALV',
+                'seq2_aa': 'MKAWQALV',
+                'seq3_aa': 'MKTWQ-LV',
+            }
+        ),
+        reference_id='ref_aa',
+    )
+
+
+@pytest.fixture
+def orf_msa() -> MSA:
+    """Create a small nucleotide alignment with one conserved ORF."""
+    return MSA(
+        create_alignment(
+            {
+                'ref': 'ATGAAATTTTAA',
+                'seq1': 'ATGAAATTTTAA',
+                'seq2': 'ATGAAGTTTTAA',
+            }
+        ),
+        reference_id='ref',
+    )
+
+
+@pytest.fixture
+def bed_annotation_file(tmp_path: Path) -> Path:
+    """Create a BED annotation file that matches the test alignment ids."""
+    annotation_file = tmp_path / 'regions.bed'
+    annotation_file.write_text(
+        'ref\t0\t4\n'
+        'ref\t6\t10\n',
+        encoding='utf-8',
+    )
+    return annotation_file
+
+
+def _assert_axis_is_populated(ax: Axes) -> None:
+    """Assert that plotting added visible artists to the axis."""
+    assert ax.collections or ax.patches or ax.lines or ax.texts
+
+
+def test_alignment_returns_axis_and_adds_collections(dna_msa: MSA, ax: Axes) -> None:
+    """Test that ``alignment`` plots onto the provided axis."""
+    returned_ax = draw.alignment(dna_msa, ax=ax, show_sequence_all=False)
+
+    assert returned_ax is ax
+    assert len(ax.collections) >= 1
+    _assert_axis_is_populated(ax)
+
+
+def test_identity_alignment_returns_axis_and_adds_collections(dna_msa: MSA, ax: Axes) -> None:
+    """Test that ``identity_alignment`` plots onto the provided axis."""
+    returned_ax = draw.identity_alignment(dna_msa, ax=ax, show_identity_sequence=True)
+
+    assert returned_ax is ax
+    assert len(ax.collections) >= 1
+    _assert_axis_is_populated(ax)
+
+
+def test_similarity_alignment_returns_axis_and_adds_collections(aa_msa: MSA, ax: Axes) -> None:
+    """Test that ``similarity_alignment`` plots onto the provided axis."""
+    returned_ax = draw.similarity_alignment(aa_msa, ax=ax, show_similarity_sequence=True)
+
+    assert returned_ax is ax
+    assert len(ax.collections) >= 1
+    _assert_axis_is_populated(ax)
+
+
+@pytest.mark.parametrize('stat_type', ['entropy', 'identity'])
+def test_stat_plot_returns_axis_for_main_stat_types(dna_msa: MSA, ax: Axes, stat_type: str) -> None:
+    """Test that ``stat_plot`` renders both scalar and matrix-based statistics."""
+    returned_ax = draw.stat_plot(dna_msa, stat_type=stat_type, ax=ax, rolling_average=1)
+
+    assert returned_ax is ax
+    assert len(ax.collections) >= 1
+    _assert_axis_is_populated(ax)
+
+
+def test_variant_plot_returns_axis_and_draws_variant_markers(dna_msa: MSA, ax: Axes) -> None:
+    """Test that ``variant_plot`` draws lollipop markers for detected variants."""
+    returned_ax = draw.variant_plot(dna_msa, ax=ax, show_legend=True)
+
+    assert returned_ax is ax
+    assert len(ax.lines) >= 1
+    assert ax.get_ylabel() == 'reference'
+    _assert_axis_is_populated(ax)
+
+
+def test_annotation_plot_returns_axis_and_draws_feature_patches(dna_msa: MSA, ax: Axes, bed_annotation_file: Path) -> None:
+    """Test that ``annotation_plot`` draws annotation boxes from a BED file."""
+    returned_ax = draw.annotation_plot(dna_msa, str(bed_annotation_file), feature_to_plot='ignored', ax=ax)
+
+    assert returned_ax is ax
+    assert len(ax.patches) >= 1
+    assert 'bed regions' in ax.get_title(loc='left')
+    _assert_axis_is_populated(ax)
+
+
+def test_orf_plot_returns_axis_and_draws_orf_patches(orf_msa: MSA, ax: Axes) -> None:
+    """Test that ``orf_plot`` draws at least one conserved ORF."""
+    returned_ax = draw.orf_plot(orf_msa, ax=ax, min_length=9, non_overlapping_orfs=True)
+
+    assert returned_ax is ax
+    assert len(ax.patches) >= 1
+    assert ax.get_title(loc='left') == 'conserved orfs'
+    _assert_axis_is_populated(ax)
+
+
+@pytest.mark.parametrize('plot_type', ['stacked', 'logo'])
+def test_sequence_logo_returns_axis_for_both_plot_modes(dna_msa: MSA, ax: Axes, plot_type: str) -> None:
+    """Test that ``sequence_logo`` supports both stacked and logo rendering."""
+    returned_ax = draw.sequence_logo(dna_msa, ax=ax, plot_type=plot_type)
+
+    assert returned_ax is ax
+    _assert_axis_is_populated(ax)
+
+
+def test_consensus_plot_returns_axis_and_draws_collection_and_text(dna_msa: MSA, ax: Axes) -> None:
+    """Test that ``consensus_plot`` draws a consensus row on the provided axis."""
+    returned_ax = draw.consensus_plot(dna_msa, ax=ax, show_sequence=True)
+
+    assert returned_ax is ax
+    assert len(ax.collections) >= 1
+    assert len(ax.texts) == dna_msa.length
+    _assert_axis_is_populated(ax)
+
+
+def test_simplot_returns_axis_and_draws_one_line_per_non_reference_sequence(dna_msa: MSA, ax: Axes) -> None:
+    """Test that ``simplot`` draws one trace per non-reference sequence."""
+    returned_ax = draw.simplot(
+        dna_msa,
+        ref='ref',
+        ax=ax,
+        window_size=4,
+        step_size=2,
+        distance_calculation='k2p',
+        show_legend=True,
+        show_reference=True,
+    )
+
+    assert returned_ax is ax
+    assert len(ax.lines) == len(dna_msa) - 1
+    assert ax.get_ylabel() == 'similarity (%)'
+    _assert_axis_is_populated(ax)
+

From d5b60b18e9907e4c94a60e0f8087c73edeb44dcd Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 20 Mar 2026 13:15:39 +0100
Subject: [PATCH 16/25] simplified variant collection building

---
 app_src/shiny_server.py      |  4 ++--
 msaexplorer/_data_classes.py |  8 +++----
 msaexplorer/explore.py       | 46 ++++++++++++++++--------------------
 tests/test_draw_plots.py     |  3 ++-
 4 files changed, 29 insertions(+), 32 deletions(-)

diff --git a/app_src/shiny_server.py b/app_src/shiny_server.py
index a806239..a5fb8c3 100644
--- a/app_src/shiny_server.py
+++ b/app_src/shiny_server.py
@@ -871,7 +871,7 @@ def zoom_range_analysis():
 
         aln = set_aln(aln, prepare_minimal_inputs())
 
-        return f'{aln.zoom[0]} - {aln.zoom[1]}'
+        return f'{aln.zoom[0]} - {aln.zoom[1]-1}'
 
     @render.ui
     def number_of_seq():
@@ -909,7 +909,7 @@ def snps():
 
         aln = set_aln(aln, prepare_minimal_inputs(ref=True))
 
-        return len(aln.get_snps()['POS'])
+        return len(aln.get_snps())
 
     @reactive.Effect
     @reactive.event(input.analysis_plot_type_left)
diff --git a/msaexplorer/_data_classes.py b/msaexplorer/_data_classes.py
index 8d8f4cd..5f6903e 100644
--- a/msaexplorer/_data_classes.py
+++ b/msaexplorer/_data_classes.py
@@ -26,10 +26,6 @@ class VariantCollection:
     chrom: str
     positions: dict[int, SingleNucleotidePolymorphism] = field(default_factory=dict)
 
-    def __post_init__(self):
-        """sort the positions by position (key)"""
-        object.__setattr__(self, 'positions', dict(sorted(self.positions.items())))
-
     def __len__(self) -> int:
         return len(self.positions)
 
@@ -38,6 +34,10 @@ def __iter__(self):
 
     def __contains__(self, position: int) -> bool:
         return position in self.positions
+    
+    def __getitem__(self, position: int) -> SingleNucleotidePolymorphism:
+        """Access a SNP by position."""
+        return self.positions[position]
 
 
 @dataclass(frozen=True)
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index aaf4fe0..3c2759c 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -1151,7 +1151,9 @@ def get_snps(self, include_ambig:bool=False) -> VariantCollection:
         ref = self._get_reference_seq()
         aln = {x: aln[x] for x in self.sequence_ids if x != self.reference_id}
         seq_ids = list(aln.keys())
-        snp_dict = {'#CHROM': self.reference_id if self.reference_id is not None else 'consensus', 'POS': {}}
+        chrom = self.reference_id if self.reference_id is not None else 'consensus'
+        snp_positions = {}
+        aln_size = len(aln)
 
         for pos in range(self.length):
             reference_char = ref[pos]
@@ -1160,11 +1162,13 @@ def get_snps(self, include_ambig:bool=False) -> VariantCollection:
                     continue
             alt_chars, snps = [], []
             for i, seq_id in enumerate(seq_ids):
-                alt_chars.append(aln[seq_id][pos])
-                if reference_char != aln[seq_id][pos]:
+                alt_char = aln[seq_id][pos]
+                alt_chars.append(alt_char)
+                if reference_char != alt_char:
                     snps.append(i)
             if not snps:
                 continue
+            # Filter out ambiguous snps if not included
             if include_ambig:
                 if all(alt_chars[x] in config.AMBIG_CHARS[self.aln_type] for x in snps):
                     continue
@@ -1172,31 +1176,23 @@ def get_snps(self, include_ambig:bool=False) -> VariantCollection:
                 snps = [x for x in snps if alt_chars[x] not in config.AMBIG_CHARS[self.aln_type]]
                 if not snps:
                     continue
-            if pos not in snp_dict['POS']:
-                snp_dict['POS'][pos] = {'ref': reference_char, 'ALT': {}}
-            for snp in snps:
-                if alt_chars[snp] not in snp_dict['POS'][pos]['ALT']:
-                    snp_dict['POS'][pos]['ALT'][alt_chars[snp]] = {
-                        'AF': 1,
-                        'SEQ_ID': [seq_ids[snp]]
-                    }
-                else:
-                    snp_dict['POS'][pos]['ALT'][alt_chars[snp]]['AF'] += 1
-                    snp_dict['POS'][pos]['ALT'][alt_chars[snp]]['SEQ_ID'].append(seq_ids[snp])
-            # calculate AF
-            if pos in snp_dict['POS']:
-                for alt in snp_dict['POS'][pos]['ALT']:
-                    snp_dict['POS'][pos]['ALT'][alt]['AF'] /= len(aln)
-
-        snp_positions = {}
-        for pos, pos_info in snp_dict['POS'].items():
+            # Build allele dict with counts
+            alt_dict = {}
+            for snp_idx in snps:
+                alt_char = alt_chars[snp_idx]
+                if alt_char not in alt_dict:
+                    alt_dict[alt_char] = {'count': 0, 'seq_ids': []}
+                alt_dict[alt_char]['count'] += 1
+                alt_dict[alt_char]['seq_ids'].append(seq_ids[snp_idx])
+            
+            # Convert to final format: alt_char -> (frequency, seq_ids_tuple)
             alt = {
-                allele: (details['AF'], tuple(details['SEQ_ID']))
-                for allele, details in pos_info['ALT'].items()
+                alt_char: (data['count'] / aln_size, tuple(data['seq_ids']))
+                for alt_char, data in alt_dict.items()
             }
-            snp_positions[pos] = SingleNucleotidePolymorphism(ref=pos_info['ref'], alt=alt)
+            snp_positions[pos] = SingleNucleotidePolymorphism(ref=reference_char, alt=alt)
 
-        return VariantCollection(chrom=snp_dict['#CHROM'], positions=snp_positions)
+        return VariantCollection(chrom=chrom, positions=snp_positions)
 
     def calc_transition_transversion_score(self) -> AlignmentStats:
         """
diff --git a/tests/test_draw_plots.py b/tests/test_draw_plots.py
index 568506c..4d1c78d 100644
--- a/tests/test_draw_plots.py
+++ b/tests/test_draw_plots.py
@@ -1,6 +1,7 @@
 """Tests for the main plotting helpers in ``msaexplorer.draw``."""
 
 from pathlib import Path
+from typing import Any, Generator
 
 import matplotlib
 matplotlib.use('Agg')
@@ -14,7 +15,7 @@
 
 
 @pytest.fixture
-def ax() -> Axes:
+def ax() -> Generator[Axes, Any, None]:
     """Create a fresh matplotlib axis for each test."""
     fig, axis = plt.subplots()
     yield axis

From df57b7e3cad8f32a2339843e3a5a54bef5d26e8c Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Fri, 20 Mar 2026 15:40:46 +0100
Subject: [PATCH 17/25] introduced dataclass for length stats

---
 msaexplorer/_data_classes.py     | 11 +++++++++++
 msaexplorer/explore.py           | 21 +++++++++++----------
 tests/test_stats_calculations.py | 23 ++++++++++++++++++++++-
 3 files changed, 44 insertions(+), 11 deletions(-)

diff --git a/msaexplorer/_data_classes.py b/msaexplorer/_data_classes.py
index 5f6903e..153f238 100644
--- a/msaexplorer/_data_classes.py
+++ b/msaexplorer/_data_classes.py
@@ -65,6 +65,17 @@ def __contains__(self, item: float) -> bool:
         return item in self.positions
 
 
+@dataclass(frozen=True)
+class LengthStats:
+    """Summary statistics for ungapped sequence lengths in an alignment."""
+
+    n_sequences: int
+    mean_length: float
+    std_length: float
+    min_length: int
+    max_length: int
+
+
 @dataclass(frozen=True)
 class PairwiseDistance:
     """
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 3c2759c..854a574 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -24,7 +24,7 @@
 
 # msaexplorer
 from msaexplorer import config
-from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, OpenReadingFrame, OrfCollection, SingleNucleotidePolymorphism, VariantCollection
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, LengthStats, OpenReadingFrame, OrfCollection, SingleNucleotidePolymorphism, VariantCollection
 from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping, _read_alignment
 
 #TODO: Move outputs to dataclasses
@@ -559,20 +559,21 @@ def get_non_overlapping_conserved_orfs(self, min_length: int = 100, identity_cut
 
         return OrfCollection(orfs=tuple(all_orfs[orf_id] for orf_id in all_orfs if orf_id in non_overlapping_ids))
 
-    def calc_length_stats(self) -> dict:
+    def calc_length_stats(self) -> LengthStats:
         """
         Determine the stats for the length of the ungapped seqs in the alignment.
-        :return: dictionary with length stats
+        :return: dataclass with length stats
         """
 
         seq_lengths = [len(self.alignment[x].replace('-', '')) for x in self.alignment]
 
-        return {'number of seq': len(self.alignment),
-                'mean length': float(np.mean(seq_lengths)),
-                'std length': float(np.std(seq_lengths)),
-                'min length': int(np.min(seq_lengths)),
-                'max length': int(np.max(seq_lengths))
-                }
+        return LengthStats(
+            n_sequences=len(self.alignment),
+            mean_length=float(np.mean(seq_lengths)),
+            std_length=float(np.std(seq_lengths)),
+            min_length=int(np.min(seq_lengths)),
+            max_length=int(np.max(seq_lengths)),
+        )
 
     def calc_entropy(self) -> AlignmentStats:
         """
@@ -923,7 +924,7 @@ def calc_position_matrix(self, matrix_type:str='PWM') -> None | ndarray | ValueE
         if matrix_type == 'IC':
             return ic
 
-    def calc_percent_recovery(self) -> dict:
+    def calc_percent_recovery(self) -> dict[str, float]:
         """
         Recovery per sequence either compared to the majority consensus seq
         or the reference seq.\n
diff --git a/tests/test_stats_calculations.py b/tests/test_stats_calculations.py
index 75eaea4..0124866 100644
--- a/tests/test_stats_calculations.py
+++ b/tests/test_stats_calculations.py
@@ -4,7 +4,28 @@
 import numpy as np
 from conftest import create_alignment
 from msaexplorer.explore import MSA
-from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, VariantCollection
+from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, LengthStats, VariantCollection
+
+class TestCalcLengthStats:
+    """Tests for calc_length_stats."""
+
+    def test_returns_length_stats_dataclass(self):
+        msa = MSA(create_alignment({'s1': 'A-CG', 's2': 'AT-G', 's3': 'ATCG'}))
+
+        result = msa.calc_length_stats()
+
+        assert isinstance(result, LengthStats)
+
+    def test_length_stats_values_are_correct(self):
+        msa = MSA(create_alignment({'s1': 'A-CG', 's2': 'AT-G', 's3': 'ATCG'}))
+
+        result = msa.calc_length_stats()
+
+        assert result.n_sequences == 3
+        assert result.mean_length == pytest.approx(3.3333333333)
+        assert result.std_length == pytest.approx(0.4714045208)
+        assert result.min_length == 3
+        assert result.max_length == 4
 
 
 class TestCalcEntropy:

From 14feabd974c51e0920b1a75e6378e5c8670aa33d Mon Sep 17 00:00:00 2001
From: Jonas Fuchs <78491186+jonas-fuchs@users.noreply.github.com>
Date: Sat, 21 Mar 2026 14:41:28 +0100
Subject: [PATCH 18/25] Update tests/test_export.py

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 tests/test_export.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/test_export.py b/tests/test_export.py
index 854f9dd..fd36610 100644
--- a/tests/test_export.py
+++ b/tests/test_export.py
@@ -156,7 +156,6 @@ def test_orf_output_is_correct(self, orf_collection):
 
     def test_orf_invalid_input_raises(self):
         with pytest.raises(ValueError, match="instance"):
-            export.orf(OrfCollection(), chrom="chr1")
             export.orf({}, chrom="chr1")
 
 

From a73ee9697163710d2fc408de90955e3bbdf07454 Mon Sep 17 00:00:00 2001
From: Jonas Fuchs <78491186+jonas-fuchs@users.noreply.github.com>
Date: Sat, 21 Mar 2026 14:46:13 +0100
Subject: [PATCH 19/25] Update msaexplorer/_data_classes.py

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 msaexplorer/_data_classes.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/msaexplorer/_data_classes.py b/msaexplorer/_data_classes.py
index 153f238..ef3a0a1 100644
--- a/msaexplorer/_data_classes.py
+++ b/msaexplorer/_data_classes.py
@@ -61,8 +61,8 @@ def __len__(self) -> int:
     def __getitem__(self, index: int) -> float:
         return self.values[index]
 
-    def __contains__(self, item: float) -> bool:
-        return item in self.positions
+    def __contains__(self, position: int) -> bool:
+        return position in self.positions
 
 
 @dataclass(frozen=True)

From b59da2490bddbd725abbfa17c1a6c5d3dfd2e9d5 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sat, 21 Mar 2026 17:21:36 +0100
Subject: [PATCH 20/25] finalized copilot review

---
 app_src/shiny_server.py | 11 +++--------
 msaexplorer/_helpers.py |  2 +-
 msaexplorer/draw.py     |  4 ++++
 msaexplorer/export.py   | 15 +++++++++------
 4 files changed, 17 insertions(+), 15 deletions(-)

diff --git a/app_src/shiny_server.py b/app_src/shiny_server.py
index a5fb8c3..9243a1f 100644
--- a/app_src/shiny_server.py
+++ b/app_src/shiny_server.py
@@ -295,8 +295,7 @@ def finalize_loaded_alignment(aln, annotation_file):
         aln.reference_id = next(iter(aln))
         reactive.alignment.set(aln)
 
-        alignment_length = aln.length - 1
-        ui.update_slider('zoom_range', max=alignment_length - 1, value=(0, alignment_length - 1))
+        ui.update_slider('zoom_range', max=aln.length - 1, value=(0, aln.length - 1))
 
         ui.remove_ui(selector="#orf_column")
         if aln.aln_type != 'AA':
@@ -685,7 +684,6 @@ def download_stats():
         """
         Download various files in standard format
         """
-
         # helper functions
         def _snp_option():
             if input.reference_2() == 'first':
@@ -730,10 +728,7 @@ def _stat_option():
                 raise ValueError('Rolling_average must be between 1 and length of alignment.')
             # define seperator
             seperator = '\t' if input.download_format() == 'tabular' else ','
-            # define which stat type to exprt
-            for stat_type in ['entropy', 'mean similarity', 'coverage', 'mean identity', 'ts tv score', 'gc']:
-                if stat_type == input.download_type():
-                    break
+            stat_type = str(input.download_type())
             # use correct function
             data = stat_functions[stat_type]()
             if isinstance(data, AlignmentStats):
@@ -741,7 +736,7 @@ def _stat_option():
             # calculate the mean for identity or similarity (identical to draw module of msaexplorer)
             else:
                 # for the mean nan values get handled as the lowest possible number in the matrix
-                data = np.nan_to_num(data, True, -1 if stat_type == 'identity' else 0)
+                data = np.nan_to_num(data, True, -1 if stat_type == 'mean identity' else 0)
                 data = np.mean(data, axis=0)
             # apply rolling average
             data = draw._moving_average(data, input.download_type_options_1(), None, aln.length)[0]
diff --git a/msaexplorer/_helpers.py b/msaexplorer/_helpers.py
index ac70f15..429803d 100644
--- a/msaexplorer/_helpers.py
+++ b/msaexplorer/_helpers.py
@@ -181,7 +181,7 @@ def k2p(seq1: str, seq2: str, aln_length: int = None) -> float:
         Returns (1 - d) * 100 as a corrected percent identity (100 = identical).
         Returns 0 when the logarithm arguments become non-positive (saturated).
         """
-        transitions = [{'A', 'G'}, {'C', 'T'}]
+        transitions = [{'A', 'G'}, {'C', 'T'}, {'C', 'U'}]
 
         ts, tv, total = 0, 0, 0
         for c1, c2 in zip(seq1, seq2):
diff --git a/msaexplorer/draw.py b/msaexplorer/draw.py
index 2ed3f3f..9d049c0 100644
--- a/msaexplorer/draw.py
+++ b/msaexplorer/draw.py
@@ -1227,6 +1227,10 @@ def simplot(aln: explore.MSA | str, ref: str | None, ax: plt.Axes | None = None,
         raise ValueError('window_size has to be a positive integer')
     if window_size > aln.length:
         raise ValueError('window_size can not be larger than the (zoomed) alignment length')
+    if not isinstance(step_size, int) or step_size <= 0:
+        raise ValueError('step_size has to be a positive integer')
+    if step_size > aln.length:
+        raise ValueError('step_size can not be larger than the (zoomed) alignment length')
     if ref is not None and ref not in aln:
         raise ValueError(f'Reference {ref} not in alignment')
     if distance_calculation not in ['ghd', 'ged', 'jc69', 'k2p']:
diff --git a/msaexplorer/export.py b/msaexplorer/export.py
index 646ab4f..a2a085f 100644
--- a/msaexplorer/export.py
+++ b/msaexplorer/export.py
@@ -13,7 +13,7 @@
 from msaexplorer._helpers import _check_and_create_path
 
 
-def snps(snp_data: VariantCollection, format_type: str = 'vcf', path: str | None = None) -> str | None | ValueError:
+def snps(snp_data: VariantCollection, format_type: str = 'vcf', path: str | None = None) -> str | None:
     """
     Export SNP data from a VariantCollection to VCF or tabular format.
 
@@ -21,7 +21,7 @@ def snps(snp_data: VariantCollection, format_type: str = 'vcf', path: str | None
     :param format_type: Format type ('vcf' or 'tabular'). Default is 'vcf'.
     :param path: Path to output VCF or tabular format. (optional)
     :return: A string containing the SNP data in the requested format.
-    :raises ValueError: if the input type is invalid or format_type is invalid.
+    :raises ValueError: If the input type is invalid or format_type is invalid.
     """
 
     def _validate():
@@ -148,7 +148,7 @@ def stats(stat_data: AlignmentStats | list | ndarray, seperator: str = '\t', pat
         return '\n'.join(lines)
 
 
-def orf(orf_dict: OrfCollection, chrom: str, path: str | None = None) -> str | ValueError | None:
+def orf(orfs: OrfCollection, chrom: str, path: str | None = None) -> str | ValueError | None:
     """
     Exports the ORF collection to a .bed file.
 
@@ -156,14 +156,17 @@ def orf(orf_dict: OrfCollection, chrom: str, path: str | None = None) -> str | V
     :param chrom: CHROM identifier for bed format.
     :param path: Path to the output .bed file.
     """
-    if not isinstance(orf_dict, OrfCollection):
-        raise ValueError('The ORF collection must be an instance of msaexplorer._data_classes.OrfContainer.')
+    if not isinstance(orfs, OrfCollection):
+        raise ValueError('The ORF collection must be an instance of msaexplorer._data_classes.OrfCollection.')
+
+    if not orfs:
+        raise ValueError('The ORF collection is empty.')
 
     _check_and_create_path(path)
 
     lines = []
 
-    for orf_id, orf_data in orf_dict.items():
+    for orf_id, orf_data in orfs.items():
         loc = orf_data.location[0]
         conservation = orf_data.conservation
         strand = orf_data.strand

From efa4efd705d3552739dfd95c7e50471d022175f4 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sat, 21 Mar 2026 19:22:05 +0100
Subject: [PATCH 21/25] fixed zoom bug

---
 app_src/shiny_server.py         |  5 +++++
 msaexplorer/explore.py          | 18 +++++++++++-------
 tests/test_alignment_parsing.py | 18 ++++++++++++++++++
 3 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/app_src/shiny_server.py b/app_src/shiny_server.py
index 9243a1f..d9b88da 100644
--- a/app_src/shiny_server.py
+++ b/app_src/shiny_server.py
@@ -824,6 +824,11 @@ def _percent_recovery_option():
                 str(error),
                 style="color: red; font-weight: bold;"
             ), duration=10)
+            # create dummy download
+            with tempfile.NamedTemporaryFile(prefix='download_error_', suffix='.txt', delete=False) as tmpfile:
+                tmpfile.write(b'')
+                tmpfile.flush()
+                return tmpfile.name
 
     @output
     @render.download
diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 854a574..729aacc 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -90,18 +90,22 @@ def _validate_zoom(zoom: tuple | int, original_aln: dict) -> ValueError | tuple
             # check if only over value is provided -> stop is alignment length
             if isinstance(zoom, int):
                 if 0 <= zoom < aln_length:
-                    return zoom, aln_length - 1
+                    return zoom, aln_length
                 else:
                     raise ValueError('Zoom start must be within the alignment length range.')
             # check if more than 2 values are provided
             if len(zoom) != 2:
                 raise ValueError('Zoom position have to be (zoom_start, zoom_end)')
-            # validate zoom start/stop
-            for position in zoom:
-                if type(position) != int:
-                    raise ValueError('Zoom positions have to be integers.')
-                if position not in range(0, aln_length):
-                    raise ValueError('Zoom position out of range')
+            start, end = zoom
+            # validate zoom start/stop for Python slicing semantics [start:end)
+            if type(start) is not int or type(end) is not int:
+                raise ValueError('Zoom positions have to be integers.')
+            if not (0 <= start < aln_length):
+                raise ValueError('Zoom position out of range')
+            if not (0 < end <= aln_length):
+                raise ValueError('Zoom position out of range')
+            if start >= end:
+                raise ValueError('Zoom position have to be (zoom_start, zoom_end) with zoom_start < zoom_end')
 
         return zoom
 
diff --git a/tests/test_alignment_parsing.py b/tests/test_alignment_parsing.py
index 3e26360..afddfaf 100644
--- a/tests/test_alignment_parsing.py
+++ b/tests/test_alignment_parsing.py
@@ -71,3 +71,21 @@ def test_read_alignment_raises_for_unparseable_content() -> None:
     """Test that MSA raises ValueError when given unparseable content."""
     with pytest.raises(ValueError, match="could not be parsed"):
         _read_alignment("this is not an alignment")
+
+
+def test_zoom_allows_exclusive_end_equal_to_alignment_length() -> None:
+    """Accept zoom ranges that use Python slicing semantics [start:end)."""
+    msa = MSA(str(DATA_DIR / 'alignment.fasta'))
+
+    msa.zoom = (0, msa.length)
+
+    assert msa.length == len(EXPECTED_ALIGNMENT['seq1'])
+
+
+def test_zoom_rejects_end_beyond_alignment_length() -> None:
+    """Reject zoom ranges whose end exceeds the alignment length."""
+    msa = MSA(str(DATA_DIR / 'alignment.fasta'))
+
+    with pytest.raises(ValueError, match='Zoom position out of range'):
+        msa.zoom = (0, msa.length + 1)
+

From 7355c5f33c5960e0234eb2ecacda32ea435bcb49 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 22 Mar 2026 11:19:26 +0100
Subject: [PATCH 22/25] removed todo

---
 msaexplorer/explore.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 729aacc..2007449 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -27,7 +27,7 @@
 from msaexplorer._data_classes import PairwiseDistance, AlignmentStats, LengthStats, OpenReadingFrame, OrfCollection, SingleNucleotidePolymorphism, VariantCollection
 from msaexplorer._helpers import _get_line_iterator, _create_distance_calculation_function_mapping, _read_alignment
 
-#TODO: Move outputs to dataclasses
+
 class MSA:
     """
     An alignment class that allows computation of several stats. Supported inputs are file paths to alignments in "fasta",

From 6aa68bc6486e5a7781343ef17c03967250337f88 Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 22 Mar 2026 11:28:31 +0100
Subject: [PATCH 23/25] smaller cleanup

---
 msaexplorer/explore.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/msaexplorer/explore.py b/msaexplorer/explore.py
index 2007449..c01530d 100644
--- a/msaexplorer/explore.py
+++ b/msaexplorer/explore.py
@@ -17,7 +17,6 @@
 # installed
 import numpy as np
 from numpy import ndarray
-from Bio import AlignIO
 from Bio import SeqIO
 from Bio.Align import MultipleSeqAlignment
 from Bio.SeqIO.InsdcIO import GenBankIterator
@@ -884,7 +883,7 @@ def calc_similarity_alignment(self, matrix_type:str|None=None, normalize:bool=Tr
 
         return similarity_array
 
-    def calc_position_matrix(self, matrix_type:str='PWM') -> None | ndarray | ValueError:
+    def calc_position_matrix(self, matrix_type:str='PWM') -> None | ndarray:
         """
         Calculates a position matrix of the specified type for the given alignment. The function
         supports generating matrices of types Position Frequency Matrix (PFM), Position Probability

From 4238d6b1321fef6b904b530c6c3199a7d0f786ec Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 22 Mar 2026 15:52:32 +0100
Subject: [PATCH 24/25] updated examples in init.py

---
 msaexplorer/__init__.py | 411 +++++++++-------------------------------
 1 file changed, 94 insertions(+), 317 deletions(-)

diff --git a/msaexplorer/__init__.py b/msaexplorer/__init__.py
index 165af04..fcc6316 100644
--- a/msaexplorer/__init__.py
+++ b/msaexplorer/__init__.py
@@ -1,382 +1,159 @@
 r"""
-# What is MSAexplorer?
+# MSAexplorer
 
-MSAexplorer is a comprehensive Python package for analyzing and visualizing multiple sequence alignments (MSAs).
-It combines powerful statistical analysis, publication-quality plotting, and flexible data export in a simple,
-dependency-light package. Perfect for both interactive use and integration into bioinformatics pipelines.
+MSAexplorer is a lightweight toolkit to **explore**, **plot**, and **export** multiple sequence alignments.
+The API is organized into three modules:
 
-## Key Features
+- `explore`: parse alignments/annotations and compute statistics.
+- `draw`: create matplotlib plots directly from `MSA` objects or file paths.
+- `export`: serialize computed results (SNPs, stats, FASTA, ORFs, recovery tables).
 
-- **Multiple input formats**: FASTA, CLUSTAL, PHYLIP, STOCKHOLM, NEXUS, or direct Biopython objects
-- **Annotation support**: GenBank, GFF3, BED formats with automatic coordinate mapping
-- **Rich statistics**: Entropy, GC content, coverage, pairwise identity, ORF detection, SNP analysis, and more
-- **Publication-ready plots**: Identity matrices, statistical plots, similarity heatmaps, and custom visualizations
-- **Flexible export**: Generate reports, export alignments in multiple formats, process with external tools
-- **Interactive web app**: Shiny-based interface with no installation required
-- **Biopython integration**: Direct support for Bio.SeqIO and Bio.AlignIO objects for seamless workflows
+All examples below are tested with files from `example_alignments/`.
 
-## Installation
+## Quick Start (`explore`)
 
-### Via pip (recommended)
-```bash
-pip install msaexplorer
-# or with optional sequence processing tools:
-pip install msaexplorer[process]  # adds pyfamsa and pytrimal support
-```
-
-### From source
-```bash
-git clone https://github.com/jonas-fuchs/MSAexplorer
-cd MSAexplorer
-pip install .  # or pip install .[process]
-```
-
-## Quick Start
-
-### As a Web Application
-
-Launch the interactive shiny app:
-```bash
-msaexplorer --run
-```
-
-Or use the web version (no installation): [GitHub Pages](https://jonas-fuchs.github.io/MSAexplorer/app/)
-
-Export as a static site:
-```bash
-pip install shinylive
-shinylive export ./ site/
-```
-
-### As a Python Package
-
-#### Basic Analysis
 ```python
+from pathlib import Path
 from msaexplorer import explore
 
-# Load alignment (supports FASTA, CLUSTAL, PHYLIP, STOCKHOLM, NEXUS)
-msa = explore.MSA('alignment.fasta')
-
-# Or from a Biopython object
-from Bio import AlignIO
-bio_alignment = AlignIO.read('alignment.fasta', 'fasta')
-msa = explore.MSA(bio_alignment)
+base = Path('example_alignments')
+aln = explore.MSA(str(base / 'DNA.fasta'))
 
-# Get basic statistics
-print(f"Sequences: {len(msa.alignment)}")
-print(f"Length: {msa.length} bp")
-print(f"Type: {msa.aln_type}")  # DNA, RNA, or AA
+# set reference and zoom window (start inclusive, end exclusive)
+aln.reference_id = aln.sequence_ids[0]
+aln.zoom = (0, 300)
 
-# Compute statistics
-entropy = msa.calc_entropy()
-gc_content = msa.calc_gc()
-coverage = msa.calc_coverage()
-pairwise_identity = msa.calc_pairwise_identity_matrix()
+print(aln.aln_type)        # 'DNA'
+print(len(aln), aln.length)
+print(aln.get_reference_coords())
 ```
 
-#### With Annotations
+### Biopython Interoperability
+
 ```python
+from pathlib import Path
+from Bio import AlignIO, SeqIO
 from msaexplorer import explore
 
-# Load alignment
-msa = explore.MSA('alignment.fasta')
+base = Path('example_alignments')
 
-# Load annotation (GenBank, GFF3, or BED format)
-annotation = explore.Annotation(msa, 'annotation.gb')
+# alignment from Bio.Align.MultipleSeqAlignment
+bio_aln = AlignIO.read(str(base / 'DNA.fasta'), 'fasta')
+aln = explore.MSA(bio_aln, reference_id=bio_aln[0].id, zoom_range=(0, 200))
 
-# Or from a Biopython GenBank iterator
-from Bio import SeqIO
-gb_iterator = SeqIO.parse('annotation.gb', 'genbank')
-annotation = explore.Annotation(msa, gb_iterator)
-
-# Access annotation features
-print(annotation.features.keys())
-print(f"Annotation type: {annotation.ann_type}")
+# annotation from Bio.SeqIO GenBank iterator
+gb_iter = SeqIO.parse(str(base / 'DNA_RNA.gb'), 'genbank')
+ann = explore.Annotation(aln, gb_iter)
+print(ann.ann_type, list(ann.features.keys())[:3])
 ```
 
-#### Advanced Analysis
+### Working with Downstream Dataclasses
+
 ```python
+from pathlib import Path
 from msaexplorer import explore
 
-msa = explore.MSA('alignment.fasta')
+aln = explore.MSA(str(Path('example_alignments') / 'DNA.fasta'), zoom_range=(0, 300))
+aln.reference_id = aln.sequence_ids[0]
 
-# Set reference and zoom range
-msa.reference_id = 'seq1'
-msa.zoom = (0, 1000)
+entropy = aln.calc_entropy()                              # AlignmentStats
+length_stats = aln.calc_length_stats()                    # LengthStats
+dist_to_ref = aln.calc_pairwise_distance_to_reference()   # PairwiseDistance
+variants = aln.get_snps()                                 # VariantCollection
+orfs = aln.get_non_overlapping_conserved_orfs(90)         # OrfCollection
 
-# Statistical analyses
-length_stats = msa.calc_length_stats()
-snps = msa.get_snps(include_ambig=False)
-consensus = msa.get_consensus(threshold=0.7, use_ambig_nt=True)
+print(entropy.stat_name, entropy.positions[:3], entropy.values[:3])
+print(length_stats.mean_length, length_stats.std_length)
+print(dist_to_ref.reference_id, dist_to_ref.sequence_ids[:2], dist_to_ref.distances[:2])
+print(variants.chrom, len(variants))
+print(orfs.keys()[:3])
+```
 
-# For nucleotide alignments
-identity_matrix = msa.calc_identity_alignment()
-similarity_matrix = msa.calc_similarity_alignment()
-position_matrix = msa.calc_position_matrix(matrix_type='PWM')
+## Plotting (`draw`)
 
-# For DNA/RNA alignments
-reverse_complement = msa.calc_reverse_complement_alignment()
-conserved_orfs = msa.get_conserved_orfs(min_length=100)
-ts_tv_score = msa.calc_transition_transversion_score()
+All plotting functions return a matplotlib `Axes`.
+You can either:
 
-# Recovery statistics
-recovery = msa.calc_percent_recovery()
-char_frequencies = msa.calc_character_frequencies()
-```
+1. pass an existing axis (best for multi-panel figures), or
+2. pass only an alignment/path for a one-liner plot.
 
-### Plotting and Visualization
+### One-Liner Examples (path input)
 
-#### Basic Identity Plot
 ```python
+from msaexplorer import draw
 import matplotlib.pyplot as plt
-from msaexplorer import explore, draw
 
-# Load and prepare alignment
-aln = explore.MSA('alignment.fasta', reference_id='seq1')
-
-# Create identity visualization
-fig, ax = plt.subplots(figsize=(14, 8))
-draw.identity_alignment(
-    aln,
-    ax,
-    show_gaps=False,
-    show_mask=True,
-    show_mismatches=True,
-    color_scheme='purine_pyrimidine',
-    show_seq_names=True,
-    show_legend=True
-)
-plt.tight_layout()
+draw.identity_alignment('example_alignments/DNA.fasta')
 plt.show()
-```
-
-#### Statistical Plots
-```python
-import matplotlib.pyplot as plt
-from msaexplorer import explore, draw
-
-aln = explore.MSA('alignment.fasta')
-
-fig, axes = plt.subplots(nrows=4, figsize=(14, 10), sharex=True)
-
-# Entropy plot
-draw.stat_plot(aln, axes[0], stat_type='entropy', rolling_average=5)
-axes[0].set_ylabel('Entropy')
-
-# GC content plot
-draw.stat_plot(aln, axes[1], stat_type='gc', rolling_average=5)
-axes[1].set_ylabel('GC Content')
-
-# Coverage plot
-draw.stat_plot(aln, axes[2], stat_type='coverage', rolling_average=5)
-axes[2].set_ylabel('Coverage')
-
-# SNP frequency
-snps = aln.get_snps()
-draw.snp_plot(aln, axes[3], snps=snps)
-axes[3].set_ylabel('SNP Count')
-
-plt.tight_layout()
+draw.stat_plot('example_alignments/DNA.fasta', stat_type='entropy', rolling_average=5)
 plt.show()
 ```
 
-#### Similarity Heatmaps
-```python
-import matplotlib.pyplot as plt
-from msaexplorer import explore, draw
-
-aln = explore.MSA('alignment.fasta', reference_id='seq1')
+### Multi-Panel Figure with Main Plot Types
 
-fig, ax = plt.subplots(figsize=(10, 8))
-draw.similarity_alignment(
-    aln,
-    ax,
-    matrix_type='similarity'  # or 'identity', 'numerical'
-)
-plt.tight_layout()
-plt.show()
-```
-
-#### Comparison Plots
 ```python
+# import necessary packages
+from pathlib import Path
 import matplotlib.pyplot as plt
 from msaexplorer import explore, draw
 
-aln1 = explore.MSA('alignment1.fasta', reference_id='seq1')
-aln2 = explore.MSA('alignment2.fasta', reference_id='seq1')
-
-fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 6))
+# Example 1
+# load DNA alignment
+base = Path('example_alignments')
+aln = explore.MSA(str(base / 'DNA.fasta'), zoom_range=(0, 100))
 
-draw.identity_alignment(aln1, ax1)
-draw.identity_alignment(aln2, ax2)
-
-ax1.set_title('Alignment 1')
-ax2.set_title('Alignment 2')
+# ini figure
+fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(14, 15), height_ratios=[1, 1, 10])
+# plot
+draw.stat_plot(aln, ax=ax[0], stat_type='coverage', rolling_average=1, show_title=True)
+draw.sequence_logo(aln, ax=ax[1], plot_type='logo')
+draw.identity_alignment(aln, ax=ax[2], show_consensus=True, show_seq_names=True, fancy_gaps=True, show_legend=True, color_scheme='standard', show_identity_sequence=True)
 
 plt.tight_layout()
 plt.show()
-```
-
-### Data Export and Processing
-
-```python
-from msaexplorer import explore, export
-
-msa = explore.MSA('alignment.fasta')
-
-# Export to different formats
-export.to_fasta(msa, 'output.fasta')
-export.to_phylip(msa, 'output.phy')
-export.to_nexus(msa, 'output.nex')
-
-# Export statistics
-stats = msa.calc_length_stats()
-# Process statistics as needed
-
-# Get alignment data
-alignment_dict = msa.alignment  # dict[seq_id: sequence]
-consensus = msa.get_consensus()
-identity_matrix = msa.calc_pairwise_identity_matrix()
-```
-
-## Input Format Support
-
-### Alignments
-Automatically detects format or specify explicitly:
-- **FASTA** - Most common, widely supported
-- **CLUSTAL** - From CLUSTAL/MUSCLE alignments
-- **PHYLIP** - PHYLIP sequential format
-- **STOCKHOLM** - Pfam/HMMER format
-- **NEXUS** - PAUP/MrBayes format
-- **Bio.Align.MultipleSeqAlignment** - Direct Biopython objects
-
-### Annotations
-Supports standard bioinformatics formats:
-- **GenBank** (.gb, .gbk) - Feature-rich annotation format
-- **GFF3** - General Feature Format v3
-- **BED** - Browser Extensible Data format
-- **Bio.SeqIO.GenBankIterator** - Direct Biopython iterators
-
-## Configuration and Customization
-
-### Set Reference and Zoom
-```python
-from msaexplorer import explore
-
-msa = explore.MSA('alignment.fasta')
-
-# Set reference sequence for identity calculations
-msa.reference_id = 'my_reference'
-
-# Focus on a region
-msa.zoom = (100, 500)  # or just start position
-msa.zoom = 100  # equals (100, alignment_end)
 
-# Reset zoom
-msa.zoom = None
-```
-
-### Color Schemes
-Available color schemes for plotting:
-- `purine_pyrimidine` - Distinguishes chemical properties
-- `nucleotide` - Standard ATCG coloring
-- `clustalx` - ClustalX color scheme
-- `taylor` - Taylor amino acid coloring
-- And more...
-
-## Workflow Examples
-
-### Pipeline: Analyze and Visualize
-```python
-from msaexplorer import explore, draw
-import matplotlib.pyplot as plt
-
-# Load data
-msa = explore.MSA('alignment.fasta')
-annotation = explore.Annotation(msa, 'annotation.gff3')
-
-# Set parameters
-msa.reference_id = next(iter(msa))
-msa.zoom = (0, 2000)
-
-# Compute statistics
-stats = {
-    'entropy': msa.calc_entropy(),
-    'gc': msa.calc_gc(),
-    'coverage': msa.calc_coverage(),
-    'snps': msa.get_snps(),
-    'identity': msa.calc_pairwise_identity_matrix()
-}
-
-# Create comprehensive figure
-fig = plt.figure(figsize=(16, 12))
-gs = fig.add_gridspec(3, 2, hspace=0.3, wspace=0.3)
+# Example 2
+# load AA alignment
+aln = explore.MSA(str(base / 'AS.fasta'))
+# ini figure
+fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(15, 5), height_ratios=[1, 2])
 
-# Plot 1: Entropy
-ax1 = fig.add_subplot(gs[0, :])
-draw.stat_plot(msa, ax1, stat_type='entropy')
-ax1.set_ylabel('Entropy')
-
-# Plot 2: Identity alignment
-ax2 = fig.add_subplot(gs[1:, :])
-draw.identity_alignment(msa, ax2, show_seq_names=True)
+draw.stat_plot(aln, ax=ax[0], stat_type='entropy', rolling_average=5, show_title=True)
+draw.identity_alignment(aln, ax=ax[1], show_consensus=True, show_seq_names=True, fancy_gaps=True, show_legend=True, color_scheme='hydrophobicity')
 
 plt.tight_layout()
 plt.show()
 ```
 
-### Pipeline: Comparative Genomics
-```python
-from msaexplorer import explore
-
-# Load multiple alignments
-alignments = {
-    'gene_a': explore.MSA('gene_a.fasta'),
-    'gene_b': explore.MSA('gene_b.fasta'),
-    'gene_c': explore.MSA('gene_c.fasta'),
-}
-
-# Compute comparative statistics
-results = {}
-for gene_name, msa in alignments.items():
-    results[gene_name] = {
-        'length_stats': msa.calc_length_stats(),
-        'entropy': msa.calc_entropy(),
-        'snps': msa.get_snps(),
-        'pairwise_identity': msa.calc_pairwise_identity_matrix()
-    }
-
-# Use results for downstream analysis
-for gene_name, stats in results.items():
-    print(f"{gene_name}:")
-    print(f"  Mean length: {stats['length_stats']['mean length']:.0f} bp")
-    print(f"  SNP count: {len(stats['snps']['POS'])}")
-```
-
-## Documentation
+## Export (`export`)
 
-For detailed API documentation, use Python's built-in help:
 ```python
-from msaexplorer import explore, draw, export
+from pathlib import Path
+from msaexplorer import explore, export
 
-help(explore.MSA)
-help(explore.Annotation)
-help(draw.identity_alignment)
-help(export.to_fasta)
-```
+aln = explore.MSA(str(Path('example_alignments') / 'DNA.fasta'), zoom_range=(0, 300))
+aln.reference_id = aln.sequence_ids[0]
 
-Or view the full documentation at the [GitHub repository](https://github.com/jonas-fuchs/MSAexplorer).
+variants = aln.get_snps()
+entropy = aln.calc_entropy()
+consensus = aln.get_consensus()
 
-## Citation
+# return as strings
+vcf_text = export.snps(variants, format_type='vcf')
+stats_text = export.stats(entropy)
+fasta_text = export.fasta(consensus, header='consensus_dna')
 
-If you use MSAexplorer in your research, please cite:
-```
-[Citation information to be added]
+# or write to disk (path argument)
+# export.snps(variants, format_type='tabular', path='results/snps')
+# export.stats(entropy, path='results/entropy.tsv')
 ```
 
-## License
+## Notes
+
+- `zoom` can be reset with `aln.zoom = None`.
+- If no `reference_id` is set, methods that need a reference use the consensus.
+- For API-level details, inspect module/class docstrings in `explore`, `draw`, and `export`.
 
-MIT License - See LICENSE file for details.
 """
 
 from importlib.metadata import version, PackageNotFoundError

From b03f0a254ee832bb4bc30c48e38cd6b1eed3bbaf Mon Sep 17 00:00:00 2001
From: jonas-fuchs <jonas.michael.fuchs@googlemail.com>
Date: Sun, 22 Mar 2026 17:02:03 +0100
Subject: [PATCH 25/25] fixed logo bug

---
 .github/workflows/deploy_page.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/deploy_page.yaml b/.github/workflows/deploy_page.yaml
index 3f3004b..a3cb434 100644
--- a/.github/workflows/deploy_page.yaml
+++ b/.github/workflows/deploy_page.yaml
@@ -50,7 +50,7 @@ jobs:
       - name: Generate API Docs
         run: |
           pdoc ./msaexplorer \
-            --logo ../logo.svg \
+            --logo https://raw.githubusercontent.com/jonas-fuchs/MSAexplorer/master/app_src/www/img/logo.svg \
             -o docs/
 
       - name: Prepare Deployment Directory