seqevolve

Simulation of nucleotide and amino acid sequences along a gene tree.

evolver

Evolve sequences along a tree using models for substitutions, indels, and rate heterogeneity.

References

1. Z. Yang. Computational molecular evolution. Oxford series in ecology and evolution. Oxford University Press, 2006. ISBN 978-0-19-856699-1 978-0-19-856702-8.

Evolver

Class for evolving sequences along a tree.

Source code in asymmetree/seqevolve/evolver.py

class Evolver:
    """Class for evolving sequences along a tree."""

    def __init__(
        self,
        subst_model: SubstModel,
        indel_model: IndelModel | None = None,
        het_model: HetModel | None = None,
        gillespie: str | bool = "auto",
        **kwargs,
    ) -> None:
        """Constructor for the Evolver class.

        Args:
            subst_model: Substitution model.
            indel_model: Model for insertions and deletions. If not provided, no indels are
                simulated.
            het_model: Model for substitution rate heterogeneity. If not provided, no rate
                heterogeneity is simulated.
            gillespie: If True, the Gillespie algorithm is used instead of constructing the exchange
                probability matrix P via matrix diagonalization. The Gillespie algorithm may be
                faster if the rate heterogeneity is large. The default is 'auto', in which case the
                exchange probability matrix is used unless the heterogeneity model has more than one
                class or is sitewise, in which case the Gillespie algorithm is used.
        """
        self.subst_model = subst_model

        self.indel_model = indel_model
        self.het_model = het_model

        if gillespie == "auto":
            if het_model and (het_model.classes > 1 or het_model.sitewise):
                self._gillespie = True
            else:
                self._gillespie = False
        else:
            self._gillespie = bool(gillespie)

    def evolve_along_tree(
        self,
        tree: Tree,
        start_length: int = 200,
        start_seq: str | None = None,
    ) -> dict[TreeNode, EvoSeq]:
        """Evolve sequence along a tree.

        If `start_seq` is not supplied, then the sequence at the root is constructed randomly using
        the stationary probabilities of the specified substitution model.

        Args:
            tree: The tree along which sequences are simulated.
            start_length: The length of the sequence at the root. The default is 200. This parameter
                is ignored if `start_seq` is provided.
            start_seq: The start sequence for the root of the tree. The default is None, in which
                case a sequence is constructed at random.

        Returns:
            A dictionary containing the TreeNode instances in the tree as keys and the simulated
                sequences as values (instances of type `EvoSeq`).
        """
        self.tree = tree
        self.site_counter = 0
        self.sequences = {}

        if start_seq is None:
            root_seq = self._random_sequence(start_length)
        else:
            root_seq = self._initialize_with_sequence(start_seq)

        if self.het_model:
            self.het_model.assign(root_seq)

        for v in tree.preorder():
            if v.parent is None:
                self.sequences[v] = root_seq
            else:
                self.sequences[v] = self._evolve(self.sequences[v.parent], v.dist)

        return self.sequences

    def get_sequences(self, include_inner: bool = True) -> dict[TreeNode, str]:
        """Return the simulated sequence as strings.

        Args:
            include_inner: If True, the sequences of the inner nodes of the tree are also included.

        Returns:
            A dictionary with the tree nodes as keys and the str sequence (without gaps) as values.
        """
        if include_inner:
            return {v: self.subst_model.to_sequence(seq) for v, seq in self.sequences.items()}

        return {
            v: self.subst_model.to_sequence(seq) for v, seq in self.sequences.items() if v.is_leaf()
        }

    def write_sequences(self, file_path: str | Path, include_inner: bool = True) -> None:
        """Write the simulated sequences into a file in fasta format.

        Args:
            file_path: The path and filename to the file into which the sequences shall be written.
            include_inner: If True, the sequences of the inner nodes of the tree are also written
                to the file.
        """
        write_fasta(file_path, self.get_sequences(include_inner=include_inner))

    def true_alignment(
        self,
        include_inner: bool = True,
        write_to: str | Path | None = None,
        alignment_format: str = "phylip",
    ) -> dict[TreeNode, str]:
        """Construct the 'true' alignment of the simulated sequences.

        The true alignment is constructed using the information about the history of substitutions
        and indels that is available for the simulated sequences.

        Args:
            include_inner: If True, the alignment also contains the sequences of the inner nodes of
                the tree.
            write_to: The path and filename to the file into which the alignment shall be written.
                The default is None, in which case nothing is written to a file.
            alignment_format: Format of the alignment that is written to a file. The default is
                'phylip'.

        Returns:
            A dictionary containing the TreeNode instances in the tree as keys and the str sequences
                as values that include the necessary gaps.
        """
        alg_builder = AlignmentBuilder(
            self.tree, self.sequences, self.subst_model.alphabet, include_inner=include_inner
        )

        alignment = alg_builder.build()

        if write_to is not None:
            write_alignment(write_to, alignment, alignment_format=alignment_format)

        return alignment

    def _random_positions(self, n: int) -> np.ndarray:
        """Draw n random positions using the stationary probabilities of the substitution model.

        Args:
            n: The number of random positions to draw.

        Returns:
            An array of length n containing the drawn positions as indices corresponding to the
                characters in the alphabet of the substitution model.
        """
        return np.random.choice(len(self.subst_model.alphabet), n, p=self.subst_model.freqs)

    def _random_sequence(self, n: int) -> EvoSeq:
        """Generate a random sequence using the stationary probabilities of the substitution model.

        Args:
            n: The length of the random sequence to generate.

        Returns:
            An instance of type `EvoSeq` containing the generated sequence.
        """
        seq = EvoSeq()

        for x in self._random_positions(n):
            seq.append(x, State.ROOT, self.site_counter)
            self.site_counter += 1

        return seq

    def _initialize_with_sequence(self, sequence: str) -> EvoSeq:
        """Initialize the sequence at the root of the tree with a given sequence.

        Args:
            sequence: The sequence to initialize the root with. This should be a string containing
                characters from the alphabet of the substitution model.

        Returns:
            An instance of type `EvoSeq` containing the initialized sequence.
        """
        seq = EvoSeq()

        for x in self.subst_model.to_indices(sequence):
            seq.append(x, State.ROOT, self.site_counter)
            self.site_counter += 1

        return seq

    def _evolve(self, parent_seq: EvoSeq, distance: float) -> EvoSeq:
        """Evolve a sequence along a branch of the tree.

        Args:
            parent_seq: The sequence at the parent node of the branch along which the sequence is
                evolved.
            distance: The evolutionary distance along the branch.

        Returns:
            An instance of type `EvoSeq` containing the evolved sequence.
        """
        child_seq = parent_seq.clone()

        if self.indel_model:
            self._generate_indels(child_seq, distance)

            if self.het_model:
                self.het_model.assign(child_seq)

        if not self._gillespie:
            # use transition prob. matrix
            self._substitute(child_seq, distance)
        else:
            # use Gillespie algorithm
            self._substitute_gillespie(child_seq, distance)

        return child_seq

    # ----------------------------------------------------------------------------------------------
    #                       Substitution by transition probability matrix
    # ----------------------------------------------------------------------------------------------

    def _substitute(self, sequence: EvoSeq, t: float) -> None:
        """Substitute characters in a sequence using the transition probability matrix.

        Args:
            sequence: The sequence to be evolved.
            t: The evolutionary time.
        """
        # (cumulative) transition probability matrices
        P = {}

        r = np.random.random(len(sequence))

        for i, site in enumerate(sequence):
            if site.status == State.INHERITED:
                # mutate according to matrix of the corresponding class

                d = site.rate_factor * t

                if d not in P:
                    if d > 0.0:
                        P[d] = np.cumsum(self.subst_model.transition_prob_matrix(d), axis=1)
                    else:
                        continue

                site._value = np.argmax(P[d][site._value, :] > r[i])

            else:
                # choose random character
                site._value = np.argmax(self.subst_model.freqs_cumulative > r[i])

    # ----------------------------------------------------------------------------------------------
    #                            Substitution by Gillespie algorithm
    # ----------------------------------------------------------------------------------------------

    def _substitute_gillespie(self, sequence: EvoSeq, t: float) -> None:
        """Substitute characters in a sequence using the Gillespie algorithm.

        Args:
            sequence: The sequence to be evolved.
            t: The evolutionary time.
        """
        for site in sequence:
            # choose random characters for insertions
            if site.status == State.INSERTION:
                site._value = np.argmax(self.subst_model.freqs_cumulative > np.random.random())

            # Gillespie algorithm for inherited positions
            else:
                current_time = 0.0
                while current_time < t:
                    # use negative value on the diagonal of the rate matrix Q
                    rate = -site.rate_factor * self.subst_model.Q[site._value, site._value]

                    current_time += np.random.exponential(1 / rate) if rate > 0.0 else float("inf")

                    if current_time >= t:
                        break

                    r_mutation = -np.random.random() * self.subst_model.Q[site._value, site._value]
                    current_sum = 0.0

                    for i in range(len(self.subst_model.alphabet)):
                        # skip the entry on diagonal
                        if i == site._value:
                            continue

                        current_sum += self.subst_model.Q[site._value, i]

                        if current_sum > r_mutation:
                            mutation = i
                            break

                    site._value = mutation

    # ----------------------------------------------------------------------------------------------
    #                                          Indels
    # ----------------------------------------------------------------------------------------------

    def _generate_indels(self, sequence: EvoSeq, t: float) -> None:
        """Generates indels using the Gillespie algorithm.

        Args:
            sequence: The sequence to be evolved.
            t: The evolutionary time.
        """
        current_time = 0.0

        while current_time < t:
            ins_rate, del_rate = self.indel_model.get_rates(len(sequence))
            total_rate = ins_rate + del_rate

            current_time += (
                np.random.exponential(1 / total_rate) if total_rate > 0.0 else float("inf")
            )

            if current_time >= t:
                break

            if np.random.random() < ins_rate / total_rate:
                self._insertion(sequence)
            else:
                self._deletion(sequence)

    def _insertion(self, sequence: EvoSeq) -> None:
        """Insert characters into a sequence.

        Args:
            sequence: The sequence to be evolved.
        """
        d = self.indel_model.draw_length()
        if d == 0:
            return

        pos = np.random.randint(-1, high=len(sequence))

        # initialize insertion before the first item
        if pos == -1:
            current_site = sequence.append_left(None, State.INSERTION, self.site_counter)
            self.site_counter += 1
            d -= 1

        # go to the element at position pos
        else:
            current_site = sequence.node_at(pos)

        for _ in range(d):
            current_site = sequence.insert_right_of(
                current_site, None, State.INSERTION, self.site_counter
            )

            self.site_counter += 1

    def _deletion(self, sequence: EvoSeq) -> None:
        """Delete characters from a sequence.

        Args:
            sequence: The sequence to be evolved.
        """
        d = self.indel_model.draw_length()
        if d == 0:
            return

        pos = np.random.randint(-d + 1, high=len(sequence))

        # deletion begins before or at the start of the sequence
        if pos <= 0:
            d += pos
            pos = 0

        sequence.remove_range(pos, d)

__init__(subst_model: SubstModel, indel_model: IndelModel | None = None, het_model: HetModel | None = None, gillespie: str | bool = 'auto', **kwargs) -> None

Constructor for the Evolver class.

Parameters:

Name	Type	Description	Default
subst_model	SubstModel	Substitution model.	required
indel_model	IndelModel | None	Model for insertions and deletions. If not provided, no indels are simulated.	None
het_model	HetModel | None	Model for substitution rate heterogeneity. If not provided, no rate heterogeneity is simulated.	None
gillespie	str | bool	If True, the Gillespie algorithm is used instead of constructing the exchange probability matrix P via matrix diagonalization. The Gillespie algorithm may be faster if the rate heterogeneity is large. The default is 'auto', in which case the exchange probability matrix is used unless the heterogeneity model has more than one class or is sitewise, in which case the Gillespie algorithm is used.	'auto'

Source code in asymmetree/seqevolve/evolver.py

def __init__(
    self,
    subst_model: SubstModel,
    indel_model: IndelModel | None = None,
    het_model: HetModel | None = None,
    gillespie: str | bool = "auto",
    **kwargs,
) -> None:
    """Constructor for the Evolver class.

    Args:
        subst_model: Substitution model.
        indel_model: Model for insertions and deletions. If not provided, no indels are
            simulated.
        het_model: Model for substitution rate heterogeneity. If not provided, no rate
            heterogeneity is simulated.
        gillespie: If True, the Gillespie algorithm is used instead of constructing the exchange
            probability matrix P via matrix diagonalization. The Gillespie algorithm may be
            faster if the rate heterogeneity is large. The default is 'auto', in which case the
            exchange probability matrix is used unless the heterogeneity model has more than one
            class or is sitewise, in which case the Gillespie algorithm is used.
    """
    self.subst_model = subst_model

    self.indel_model = indel_model
    self.het_model = het_model

    if gillespie == "auto":
        if het_model and (het_model.classes > 1 or het_model.sitewise):
            self._gillespie = True
        else:
            self._gillespie = False
    else:
        self._gillespie = bool(gillespie)

evolve_along_tree(tree: Tree, start_length: int = 200, start_seq: str | None = None) -> dict[TreeNode, EvoSeq]

Evolve sequence along a tree.

If start_seq is not supplied, then the sequence at the root is constructed randomly using the stationary probabilities of the specified substitution model.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree along which sequences are simulated.	required
start_length	int	The length of the sequence at the root. The default is 200. This parameter is ignored if start_seq is provided.	200
start_seq	str | None	The start sequence for the root of the tree. The default is None, in which case a sequence is constructed at random.	None

Returns:

Type	Description
dict[TreeNode, EvoSeq]	A dictionary containing the TreeNode instances in the tree as keys and the simulated sequences as values (instances of type EvoSeq).

Source code in asymmetree/seqevolve/evolver.py

def evolve_along_tree(
    self,
    tree: Tree,
    start_length: int = 200,
    start_seq: str | None = None,
) -> dict[TreeNode, EvoSeq]:
    """Evolve sequence along a tree.

    If `start_seq` is not supplied, then the sequence at the root is constructed randomly using
    the stationary probabilities of the specified substitution model.

    Args:
        tree: The tree along which sequences are simulated.
        start_length: The length of the sequence at the root. The default is 200. This parameter
            is ignored if `start_seq` is provided.
        start_seq: The start sequence for the root of the tree. The default is None, in which
            case a sequence is constructed at random.

    Returns:
        A dictionary containing the TreeNode instances in the tree as keys and the simulated
            sequences as values (instances of type `EvoSeq`).
    """
    self.tree = tree
    self.site_counter = 0
    self.sequences = {}

    if start_seq is None:
        root_seq = self._random_sequence(start_length)
    else:
        root_seq = self._initialize_with_sequence(start_seq)

    if self.het_model:
        self.het_model.assign(root_seq)

    for v in tree.preorder():
        if v.parent is None:
            self.sequences[v] = root_seq
        else:
            self.sequences[v] = self._evolve(self.sequences[v.parent], v.dist)

    return self.sequences

get_sequences(include_inner: bool = True) -> dict[TreeNode, str]

Return the simulated sequence as strings.

Parameters:

Name	Type	Description	Default
include_inner	bool	If True, the sequences of the inner nodes of the tree are also included.	True

Returns:

Type	Description
dict[TreeNode, str]	A dictionary with the tree nodes as keys and the str sequence (without gaps) as values.

Source code in asymmetree/seqevolve/evolver.py

def get_sequences(self, include_inner: bool = True) -> dict[TreeNode, str]:
    """Return the simulated sequence as strings.

    Args:
        include_inner: If True, the sequences of the inner nodes of the tree are also included.

    Returns:
        A dictionary with the tree nodes as keys and the str sequence (without gaps) as values.
    """
    if include_inner:
        return {v: self.subst_model.to_sequence(seq) for v, seq in self.sequences.items()}

    return {
        v: self.subst_model.to_sequence(seq) for v, seq in self.sequences.items() if v.is_leaf()
    }

write_sequences(file_path: str | Path, include_inner: bool = True) -> None

Write the simulated sequences into a file in fasta format.

Parameters:

Name	Type	Description	Default
file_path	str | Path	The path and filename to the file into which the sequences shall be written.	required
include_inner	bool	If True, the sequences of the inner nodes of the tree are also written to the file.	True

Source code in asymmetree/seqevolve/evolver.py

def write_sequences(self, file_path: str | Path, include_inner: bool = True) -> None:
    """Write the simulated sequences into a file in fasta format.

    Args:
        file_path: The path and filename to the file into which the sequences shall be written.
        include_inner: If True, the sequences of the inner nodes of the tree are also written
            to the file.
    """
    write_fasta(file_path, self.get_sequences(include_inner=include_inner))

true_alignment(include_inner: bool = True, write_to: str | Path | None = None, alignment_format: str = 'phylip') -> dict[TreeNode, str]

Construct the 'true' alignment of the simulated sequences.

The true alignment is constructed using the information about the history of substitutions and indels that is available for the simulated sequences.

Parameters:

Name	Type	Description	Default
include_inner	bool	If True, the alignment also contains the sequences of the inner nodes of the tree.	True
write_to	str | Path | None	The path and filename to the file into which the alignment shall be written. The default is None, in which case nothing is written to a file.	None
alignment_format	str	Format of the alignment that is written to a file. The default is 'phylip'.	'phylip'

Returns:

Type	Description
dict[TreeNode, str]	A dictionary containing the TreeNode instances in the tree as keys and the str sequences as values that include the necessary gaps.

Source code in asymmetree/seqevolve/evolver.py

def true_alignment(
    self,
    include_inner: bool = True,
    write_to: str | Path | None = None,
    alignment_format: str = "phylip",
) -> dict[TreeNode, str]:
    """Construct the 'true' alignment of the simulated sequences.

    The true alignment is constructed using the information about the history of substitutions
    and indels that is available for the simulated sequences.

    Args:
        include_inner: If True, the alignment also contains the sequences of the inner nodes of
            the tree.
        write_to: The path and filename to the file into which the alignment shall be written.
            The default is None, in which case nothing is written to a file.
        alignment_format: Format of the alignment that is written to a file. The default is
            'phylip'.

    Returns:
        A dictionary containing the TreeNode instances in the tree as keys and the str sequences
            as values that include the necessary gaps.
    """
    alg_builder = AlignmentBuilder(
        self.tree, self.sequences, self.subst_model.alphabet, include_inner=include_inner
    )

    alignment = alg_builder.build()

    if write_to is not None:
        write_alignment(write_to, alignment, alignment_format=alignment_format)

    return alignment

evolving_sequence

Data structure for evolving sequences.

The EvoSeq class represents a sequence that evolves along a tree. It is implemented as a doubly linked list, which allows for efficient insertions and deletions of sites, which are necessary to simulate indels. Each site in the sequence is represented by an EvoSeqNode, which contains information about the state of the site (whether it is inherited from the parent sequence, inserted, or part of the root sequence), its site ID, and its parent element in the case of inherited sites.

State

Bases: Enum

State of the sites in evolving sequences.

Source code in asymmetree/seqevolve/evolving_sequence.py

class State(Enum):
    """State of the sites in evolving sequences."""

    ROOT = auto()
    INSERTION = auto()
    INHERITED = auto()

EvoSeqNode

Bases: DLListNode

Site in an evolving sequence.

Source code in asymmetree/seqevolve/evolving_sequence.py

class EvoSeqNode(DLListNode):
    """Site in an evolving sequence."""

    __slots__ = ("status", "site_id", "parent_el", "rate_class", "rate_factor")

    def __init__(
        self,
        value: int,
        status: State,
        site_id: int,
        prev_node: EvoSeqNode | None = None,
        next_node: EvoSeqNode | None = None,
        parent_el: EvoSeqNode | None = None,
        rate_class: int = 0,
        rate_factor: float = 1.0,
    ) -> None:
        """Constructor for the EvoSeqNode class.

        Args:
            value: The character at the site represented as an integer (index in the alphabet) or
                None for gaps.
            status: The state of the site (whether it is inherited from the parent sequence,
                inserted, or part of the root sequence).
            site_id: The site ID, which is a unique identifier for the site across the entire tree.
            prev_node: The previous site in the sequence.
            next_node: The next site in the sequence.
        """
        super().__init__(value, prev_node=prev_node, next_node=next_node)

        self.status = status
        self.site_id = site_id
        self.parent_el = parent_el
        self.rate_class = rate_class
        self.rate_factor = rate_factor

__init__(value: int, status: State, site_id: int, prev_node: EvoSeqNode | None = None, next_node: EvoSeqNode | None = None, parent_el: EvoSeqNode | None = None, rate_class: int = 0, rate_factor: float = 1.0) -> None

Constructor for the EvoSeqNode class.

Parameters:

Name	Type	Description	Default
value	int	The character at the site represented as an integer (index in the alphabet) or None for gaps.	required
status	State	The state of the site (whether it is inherited from the parent sequence, inserted, or part of the root sequence).	required
site_id	int	The site ID, which is a unique identifier for the site across the entire tree.	required
prev_node	EvoSeqNode | None	The previous site in the sequence.	None
next_node	EvoSeqNode | None	The next site in the sequence.	None

Source code in asymmetree/seqevolve/evolving_sequence.py

def __init__(
    self,
    value: int,
    status: State,
    site_id: int,
    prev_node: EvoSeqNode | None = None,
    next_node: EvoSeqNode | None = None,
    parent_el: EvoSeqNode | None = None,
    rate_class: int = 0,
    rate_factor: float = 1.0,
) -> None:
    """Constructor for the EvoSeqNode class.

    Args:
        value: The character at the site represented as an integer (index in the alphabet) or
            None for gaps.
        status: The state of the site (whether it is inherited from the parent sequence,
            inserted, or part of the root sequence).
        site_id: The site ID, which is a unique identifier for the site across the entire tree.
        prev_node: The previous site in the sequence.
        next_node: The next site in the sequence.
    """
    super().__init__(value, prev_node=prev_node, next_node=next_node)

    self.status = status
    self.site_id = site_id
    self.parent_el = parent_el
    self.rate_class = rate_class
    self.rate_factor = rate_factor

EvoSeq

Bases: DLList

Evolving sequence.

The EvoSeq class represents a sequence that evolves along a tree. It is implemented as a doubly linked list, which allows for efficient insertions and deletions of sites, which are necessary to simulate indels.

Source code in asymmetree/seqevolve/evolving_sequence.py

class EvoSeq(DLList):
    """Evolving sequence.

    The EvoSeq class represents a sequence that evolves along a tree. It is implemented as a doubly
    linked list, which allows for efficient insertions and deletions of sites, which are necessary
    to simulate indels.
    """

    def __init__(self) -> None:
        """Constructor for the EvoSeq class."""
        super().__init__()

    def __iter__(self):
        """Iterator for the sites in the evolving sequence."""
        return EvoSeqIterator(self)

    def __next__(self):
        """Raise an error if next is called on the EvoSeq object itself."""
        raise TypeError("EvoSeq object is not an iterator. Use iter(EvoSeq) to get an iterator.")

    def append(
        self,
        value: int | None,
        status: State,
        site_id: int,
        parent_el: EvoSeqNode | None = None,
        rate_class: int = 0,
        rate_factor: float = 1.0,
    ) -> EvoSeqNode:
        """Append a new site to the end of the evolving sequence.

        Args:
            value: The character at the site represented as an integer (index in the alphabet) or
                None for gaps.
            status: The state of the site (whether it is inherited from the parent sequence,
                inserted, or part of the root sequence).
            site_id: The site ID, which is a unique identifier for the site across the entire tree.
            parent_el: The parent element in the case of inherited sites.
            rate_class: The rate class of the site.
            rate_factor: The rate factor of the site.

        Returns:
            The newly appended `EvoSeqNode`.
        """
        new_end = EvoSeqNode(
            value,
            status,
            site_id,
            prev_node=self._last,
            parent_el=parent_el,
            rate_class=rate_class,
            rate_factor=rate_factor,
        )
        if self._last:
            self._last._next = new_end
        self._last = new_end
        if not self._first:
            self._first = new_end
        self._count += 1

        return new_end

    def append_left(
        self,
        value: int | None,
        status: State,
        site_id: int,
        parent_el: EvoSeqNode | None = None,
        rate_class: int = 0,
        rate_factor: float = 1.0,
    ) -> EvoSeqNode:
        """Append a new site to the beginning of the evolving sequence.

        Args:
            value: The character at the site represented as an integer (index in the alphabet) or
                None for gaps.
            status: The state of the site (whether it is inherited from the parent sequence,
                inserted, or part of the root sequence).
            site_id: The site ID, which is a unique identifier for the site across the entire tree.
            parent_el: The parent element in the case of inherited sites.
            rate_class: The rate class of the site.
            rate_factor: The rate factor of the site.

        Returns:
            The newly appended `EvoSeqNode` (at the beginning of the sequence).
        """
        new_start = EvoSeqNode(
            value,
            status,
            site_id,
            next_node=self._first,
            parent_el=parent_el,
            rate_class=rate_class,
            rate_factor=rate_factor,
        )
        if self._first:
            self._first._prev = new_start
        self._first = new_start
        if not self._last:
            self._last = new_start
        self._count += 1

        return new_start

    def insert_right_of(
        self,
        element: EvoSeqNode,
        value: int | None,
        status: State,
        site_id: int,
        parent_el: EvoSeqNode | None = None,
        rate_class: int = 0,
        rate_factor: float = 1.0,
    ) -> EvoSeqNode:
        """Insert a new site to the right of the given element in the evolving sequence.

        Args:
            element: The element to the left of the new site.
            value: The character at the site represented as an integer (index in the alphabet) or
                None for gaps.
            status: The state of the site (whether it is inherited from the parent sequence,
                inserted, or part of the root sequence).
            site_id: The site ID, which is a unique identifier for the site across the entire tree.
            parent_el: The parent element in the case of inherited sites.
            rate_class: The rate class of the site.
            rate_factor: The rate factor of the site.

        Returns:
            The newly inserted `EvoSeqNode`.
        """
        if element is self._last:
            new_element = self.append(
                value,
                status,
                site_id,
                parent_el=parent_el,
                rate_class=rate_class,
                rate_factor=rate_factor,
            )

        else:
            new_element = EvoSeqNode(
                value,
                status,
                site_id,
                prev_node=element,
                next_node=element._next,
                parent_el=parent_el,
                rate_class=rate_class,
                rate_factor=rate_factor,
            )
            new_element._next._prev = new_element
            element._next = new_element
            self._count += 1

        return new_element

    def clone(self) -> EvoSeq:
        """Clone the evolving sequence to create a child sequence from a parent sequence.

        All sites in the child sequence are marked as inherited, and they reference the
        corresponding sites in the parent sequence.

        Returns:
            A new `EvoSeq` that is a clone of the current sequence.
        """
        child_seq = EvoSeq()

        parent_site = self._first
        while parent_site:
            child_seq.append(
                parent_site._value,
                State.INHERITED,
                parent_site.site_id,
                parent_el=parent_site,
                rate_class=parent_site.rate_class,
                rate_factor=parent_site.rate_factor,
            )
            parent_site = parent_site._next

        return child_seq

    def element_pairs(self) -> Iterator[tuple[EvoSeqNode, EvoSeqNode]]:
        """Generator for subsequent element pairs.

        Yields:
            Tuples of subsequent elements in the evolving sequence.
        """
        if self._count >= 2:
            first = self._first
            second = self._first._next

            while second:
                yield (first, second)
                first = second
                second = second._next

    def count_status(self, status: State) -> int:
        """Count the number of sites with a given status.

        Args:
            status: The status to count.

        Returns:
            The number of sites with the given status.
        """
        counter = 0

        site = self._first
        while site:
            if site.status == status:
                counter += 1
            site = site._next

        return counter

__init__() -> None

Constructor for the EvoSeq class.

Source code in asymmetree/seqevolve/evolving_sequence.py

def __init__(self) -> None:
    """Constructor for the EvoSeq class."""
    super().__init__()

__iter__()

Iterator for the sites in the evolving sequence.

Source code in asymmetree/seqevolve/evolving_sequence.py

def __iter__(self):
    """Iterator for the sites in the evolving sequence."""
    return EvoSeqIterator(self)

__next__()

Raise an error if next is called on the EvoSeq object itself.

Source code in asymmetree/seqevolve/evolving_sequence.py

def __next__(self):
    """Raise an error if next is called on the EvoSeq object itself."""
    raise TypeError("EvoSeq object is not an iterator. Use iter(EvoSeq) to get an iterator.")

append(value: int | None, status: State, site_id: int, parent_el: EvoSeqNode | None = None, rate_class: int = 0, rate_factor: float = 1.0) -> EvoSeqNode

Append a new site to the end of the evolving sequence.

Parameters:

Name	Type	Description	Default
value	int | None	The character at the site represented as an integer (index in the alphabet) or None for gaps.	required
status	State	The state of the site (whether it is inherited from the parent sequence, inserted, or part of the root sequence).	required
site_id	int	The site ID, which is a unique identifier for the site across the entire tree.	required
parent_el	EvoSeqNode | None	The parent element in the case of inherited sites.	None
rate_class	int	The rate class of the site.	0
rate_factor	float	The rate factor of the site.	1.0

Returns:

Type	Description
EvoSeqNode	The newly appended EvoSeqNode.

Source code in asymmetree/seqevolve/evolving_sequence.py

def append(
    self,
    value: int | None,
    status: State,
    site_id: int,
    parent_el: EvoSeqNode | None = None,
    rate_class: int = 0,
    rate_factor: float = 1.0,
) -> EvoSeqNode:
    """Append a new site to the end of the evolving sequence.

    Args:
        value: The character at the site represented as an integer (index in the alphabet) or
            None for gaps.
        status: The state of the site (whether it is inherited from the parent sequence,
            inserted, or part of the root sequence).
        site_id: The site ID, which is a unique identifier for the site across the entire tree.
        parent_el: The parent element in the case of inherited sites.
        rate_class: The rate class of the site.
        rate_factor: The rate factor of the site.

    Returns:
        The newly appended `EvoSeqNode`.
    """
    new_end = EvoSeqNode(
        value,
        status,
        site_id,
        prev_node=self._last,
        parent_el=parent_el,
        rate_class=rate_class,
        rate_factor=rate_factor,
    )
    if self._last:
        self._last._next = new_end
    self._last = new_end
    if not self._first:
        self._first = new_end
    self._count += 1

    return new_end

append_left(value: int | None, status: State, site_id: int, parent_el: EvoSeqNode | None = None, rate_class: int = 0, rate_factor: float = 1.0) -> EvoSeqNode

Append a new site to the beginning of the evolving sequence.

Parameters:

Name	Type	Description	Default
value	int | None	The character at the site represented as an integer (index in the alphabet) or None for gaps.	required
status	State	The state of the site (whether it is inherited from the parent sequence, inserted, or part of the root sequence).	required
site_id	int	The site ID, which is a unique identifier for the site across the entire tree.	required
parent_el	EvoSeqNode | None	The parent element in the case of inherited sites.	None
rate_class	int	The rate class of the site.	0
rate_factor	float	The rate factor of the site.	1.0

Returns:

Type	Description
EvoSeqNode	The newly appended EvoSeqNode (at the beginning of the sequence).

Source code in asymmetree/seqevolve/evolving_sequence.py

def append_left(
    self,
    value: int | None,
    status: State,
    site_id: int,
    parent_el: EvoSeqNode | None = None,
    rate_class: int = 0,
    rate_factor: float = 1.0,
) -> EvoSeqNode:
    """Append a new site to the beginning of the evolving sequence.

    Args:
        value: The character at the site represented as an integer (index in the alphabet) or
            None for gaps.
        status: The state of the site (whether it is inherited from the parent sequence,
            inserted, or part of the root sequence).
        site_id: The site ID, which is a unique identifier for the site across the entire tree.
        parent_el: The parent element in the case of inherited sites.
        rate_class: The rate class of the site.
        rate_factor: The rate factor of the site.

    Returns:
        The newly appended `EvoSeqNode` (at the beginning of the sequence).
    """
    new_start = EvoSeqNode(
        value,
        status,
        site_id,
        next_node=self._first,
        parent_el=parent_el,
        rate_class=rate_class,
        rate_factor=rate_factor,
    )
    if self._first:
        self._first._prev = new_start
    self._first = new_start
    if not self._last:
        self._last = new_start
    self._count += 1

    return new_start

insert_right_of(element: EvoSeqNode, value: int | None, status: State, site_id: int, parent_el: EvoSeqNode | None = None, rate_class: int = 0, rate_factor: float = 1.0) -> EvoSeqNode

Insert a new site to the right of the given element in the evolving sequence.

Parameters:

Name	Type	Description	Default
element	EvoSeqNode	The element to the left of the new site.	required
value	int | None	The character at the site represented as an integer (index in the alphabet) or None for gaps.	required
status	State	The state of the site (whether it is inherited from the parent sequence, inserted, or part of the root sequence).	required
site_id	int	The site ID, which is a unique identifier for the site across the entire tree.	required
parent_el	EvoSeqNode | None	The parent element in the case of inherited sites.	None
rate_class	int	The rate class of the site.	0
rate_factor	float	The rate factor of the site.	1.0

Returns:

Type	Description
EvoSeqNode	The newly inserted EvoSeqNode.

Source code in asymmetree/seqevolve/evolving_sequence.py

def insert_right_of(
    self,
    element: EvoSeqNode,
    value: int | None,
    status: State,
    site_id: int,
    parent_el: EvoSeqNode | None = None,
    rate_class: int = 0,
    rate_factor: float = 1.0,
) -> EvoSeqNode:
    """Insert a new site to the right of the given element in the evolving sequence.

    Args:
        element: The element to the left of the new site.
        value: The character at the site represented as an integer (index in the alphabet) or
            None for gaps.
        status: The state of the site (whether it is inherited from the parent sequence,
            inserted, or part of the root sequence).
        site_id: The site ID, which is a unique identifier for the site across the entire tree.
        parent_el: The parent element in the case of inherited sites.
        rate_class: The rate class of the site.
        rate_factor: The rate factor of the site.

    Returns:
        The newly inserted `EvoSeqNode`.
    """
    if element is self._last:
        new_element = self.append(
            value,
            status,
            site_id,
            parent_el=parent_el,
            rate_class=rate_class,
            rate_factor=rate_factor,
        )

    else:
        new_element = EvoSeqNode(
            value,
            status,
            site_id,
            prev_node=element,
            next_node=element._next,
            parent_el=parent_el,
            rate_class=rate_class,
            rate_factor=rate_factor,
        )
        new_element._next._prev = new_element
        element._next = new_element
        self._count += 1

    return new_element

clone() -> EvoSeq

Clone the evolving sequence to create a child sequence from a parent sequence.

All sites in the child sequence are marked as inherited, and they reference the corresponding sites in the parent sequence.

Returns:

Type	Description
EvoSeq	A new EvoSeq that is a clone of the current sequence.

Source code in asymmetree/seqevolve/evolving_sequence.py

def clone(self) -> EvoSeq:
    """Clone the evolving sequence to create a child sequence from a parent sequence.

    All sites in the child sequence are marked as inherited, and they reference the
    corresponding sites in the parent sequence.

    Returns:
        A new `EvoSeq` that is a clone of the current sequence.
    """
    child_seq = EvoSeq()

    parent_site = self._first
    while parent_site:
        child_seq.append(
            parent_site._value,
            State.INHERITED,
            parent_site.site_id,
            parent_el=parent_site,
            rate_class=parent_site.rate_class,
            rate_factor=parent_site.rate_factor,
        )
        parent_site = parent_site._next

    return child_seq

element_pairs() -> Iterator[tuple[EvoSeqNode, EvoSeqNode]]

Generator for subsequent element pairs.

Yields:

Type	Description
tuple[EvoSeqNode, EvoSeqNode]	Tuples of subsequent elements in the evolving sequence.

Source code in asymmetree/seqevolve/evolving_sequence.py

def element_pairs(self) -> Iterator[tuple[EvoSeqNode, EvoSeqNode]]:
    """Generator for subsequent element pairs.

    Yields:
        Tuples of subsequent elements in the evolving sequence.
    """
    if self._count >= 2:
        first = self._first
        second = self._first._next

        while second:
            yield (first, second)
            first = second
            second = second._next

count_status(status: State) -> int

Count the number of sites with a given status.

Parameters:

Name	Type	Description	Default
status	State	The status to count.	required

Returns:

Type	Description
int	The number of sites with the given status.

Source code in asymmetree/seqevolve/evolving_sequence.py

def count_status(self, status: State) -> int:
    """Count the number of sites with a given status.

    Args:
        status: The status to count.

    Returns:
        The number of sites with the given status.
    """
    counter = 0

    site = self._first
    while site:
        if site.status == status:
            counter += 1
        site = site._next

    return counter

EvoSeqIterator

Bases: DLListIterator

Iterator class for evolving sequence.

Source code in asymmetree/seqevolve/evolving_sequence.py

class EvoSeqIterator(DLListIterator):
    """Iterator class for evolving sequence."""

    def __init__(self, evoseq: EvoSeq) -> None:
        """Constructor for the EvoSeqIterator class.

        Args:
            evoseq: The evolving sequence to iterate over.
        """
        super().__init__(evoseq)

    def __next__(self) -> EvoSeqNode:
        """Next site element in the evolving sequence.

        Overrides the super class methods and returns the element instead of the value only.

        Returns:
            The next `EvoSeqNode` in the evolving sequence.
        """
        if self._current:
            x = self._current
            self._current = self._current._next
            return x
        else:
            raise StopIteration

__init__(evoseq: EvoSeq) -> None

Constructor for the EvoSeqIterator class.

Parameters:

Name	Type	Description	Default
evoseq	EvoSeq	The evolving sequence to iterate over.	required

Source code in asymmetree/seqevolve/evolving_sequence.py

def __init__(self, evoseq: EvoSeq) -> None:
    """Constructor for the EvoSeqIterator class.

    Args:
        evoseq: The evolving sequence to iterate over.
    """
    super().__init__(evoseq)

__next__() -> EvoSeqNode

Next site element in the evolving sequence.

Overrides the super class methods and returns the element instead of the value only.

Returns:

Type	Description
EvoSeqNode	The next EvoSeqNode in the evolving sequence.

Source code in asymmetree/seqevolve/evolving_sequence.py

def __next__(self) -> EvoSeqNode:
    """Next site element in the evolving sequence.

    Overrides the super class methods and returns the element instead of the value only.

    Returns:
        The next `EvoSeqNode` in the evolving sequence.
    """
    if self._current:
        x = self._current
        self._current = self._current._next
        return x
    else:
        raise StopIteration

substitution

Substitution models for nucleotide and amino acid evolution.

The module contains the class SubstModel that can be used to represent a substitution model for nucleotide or amino acid evolution. The class provides methods to compute the transition probability matrix for a given time and to convert sequences between their str representation and a list of indices.

For amino acid substitution models, some empirical models are available that can be loaded via the EMPIRICAL_MODELS object. The values in the matrices and the equilibrium frequencies are taken from [2].

References

1. Z. Yang. Computational molecular evolution. Oxford series in ecology and evolution. Oxford University Press, 2006. ISBN 978-0-19-856699-1 978-0-19-856702-8.
2. http://giphy.pasteur.fr/empirical-models-of-amino-acid-substitution/

SubstModel

Substitution model for nucleotide or amino acid evolution.

Source code in asymmetree/seqevolve/substitution.py

class SubstModel:
    """Substitution model for nucleotide or amino acid evolution."""

    nuc_models = {"JC69", "K80", "GTR", "CUSTOM"}
    aa_models = {"JC69", "CUSTOM"}

    nucleotides = "ACGT"
    amino_acids = "ARNDCQEGHILKMFPSTWYV"

    def __init__(self, model_type: str, model_name: str, **kwargs) -> None:
        """Constructor for the SubstModel class.

        Some models require additional parameters, e.g. 'K80' requires the parameter 'kappa' and
        'GTR' requires 'abcdef' and 'f', see documentation for details. For a custom model, the
        keyword parameter 'filename' has to be specified with a path to a custom model in paml
        format. All of these parameters are passed as keyword arguments.

        Args:
            model_type: Available options are 'n' for nucleotide and 'a' for amino acid sequences.
            model_name: Available are e.g. 'JC69' (nuc./a.a.), 'K80' (nuc.), 'GTR' (nuc.), 'WAG'
                (a.a.), 'JTT' (a.a.), 'BLOSUM62' (a.a.), 'LG' (a.a.), and 'DAYHOFF' (a.a.). For the
                option 'CUSTOM', the keyword parameter 'filename' has to be specified with a path to
                a custom model in paml format.

        Raises:
            ValueError: If an invalid model was specified.
        """
        self.model_type = model_type.lower()
        self.model_name = model_name.upper()

        self._params = kwargs

        if (
            self.model_type in ("n", "nuc", "nucleotide")
            and self.model_name in SubstModel.nuc_models
        ):
            self.model_type = "n"
            self.alphabet = SubstModel.nucleotides

        elif self.model_type in ("a", "aa", "amino", "aminoacid", "protein") and (
            self.model_name in SubstModel.aa_models or self.model_name in EMPIRICAL_MODELS
        ):
            self.model_type = "a"
            self.alphabet = SubstModel.amino_acids

        else:
            raise ValueError(f"model '{model_type}', '{model_name}' is not available")

        self.alphabet_dict = {item: index for index, item in enumerate(self.alphabet)}

        self._load_exchangeability_and_freqs()
        self._build_rate_matrix()

    def _load_exchangeability_and_freqs(self):
        """Load the exchangeability matrix S and the stationary frequencies π.

        The exchangeability matrix S and the stationary frequencies π are loaded according to the
        specified model. For custom models, they are loaded from a provided PAML file.
        """
        # a custom model (via a paml file) was specified
        if self.model_name == "CUSTOM":
            if self._params is None or "filename" not in self._params:
                raise ValueError("custom model requires the parameter 'filename'")

            S, freqs = parse_paml(self._params["filename"], model_type=self.model_type)
            self.S, self.freqs = np.asarray(S), np.asarray(freqs)

        # non-empirical nucleotide models
        elif self.model_type == "n":
            if self.model_name == "JC69":
                self.S, self.freqs = _JC69_nuc()

            elif self.model_name == "K80":
                if self._params is None or "kappa" not in self._params:
                    raise ValueError("model 'K80' requires the parameter 'kappa'")

                self.S, self.freqs = _K80_nuc(self._params["kappa"])

            elif self.model_name == "GTR":
                if self._params is None or "abcdef" not in self._params or "f" not in self._params:
                    raise ValueError("model 'GTR' requires the parameters 'abcdef' and 'f'")

                self.S, self.freqs = _GTR_nuc(self._params["abcdef"], self._params["f"])

        # non-empirical and empirical amino acid models
        elif self.model_type == "a":
            if self.model_name == "JC69":
                self.S, self.freqs = _JC69_aa()
            elif self.model_name in EMPIRICAL_MODELS:
                self.S, self.freqs = EMPIRICAL_MODELS[self.model_name]

        # make sure stationary frequencies sum to 1
        self.freqs /= np.sum(self.freqs)

        # compute cumulative frequencies for evolver
        self.freqs_cumulative = np.cumsum(self.freqs)

    def _build_rate_matrix(self):
        """Compute the rate matrix Q.

        The rate matrix Q is computed from the exchangeability matrix S and the stationary
        frequencies π according to the formula Q[i, j] = S[i, j] * π[j] for i != j and
        Q[i, i] = -sum(Q[i, :]). The matrix is then normalized so that the expected number of
        substitutions per unit time is 1.
        """
        # rate matrix from exchangeability matrix and stationary frequencies
        Q = self.S @ np.diag(self.freqs)

        # compute diagonals, rows sum to 0
        for i in range(Q.shape[0]):
            Q[i, i] = -(np.sum(Q[i, :]) - Q[i, i])

        # normalize matrix
        k = -np.dot(self.freqs, np.diagonal(Q))
        Q /= k

        self.Q = Q

    def eigensystem(self):
        """Compute the eigensystem of the rate matrix Q."""
        if not hasattr(self, "eigenvals"):
            self.eigenvals, self.U, self.U_inv = diagonalize(self.Q, self.freqs)

        return self.eigenvals, self.U, self.U_inv

    def transition_prob_matrix(self, t: float) -> np.ndarray:
        """Calculate the transition probability matrix P(t).

        The transition probability matrix P(t) is computed from the eigensystem of the rate matrix Q
        according to the formula P(t) = U x e^(Λ * t) x U^(-1) where U is the matrix of eigenvectors
        and Λ is the diagonal matrix of eigenvalues.

        Args:
            t: The time / evolutionary distance for which to compute the transition probability
                matrix.

        Returns:
            The transition probability matrix P(t) as a `numpy.ndarray`.
        """
        # ensure that eigensystem has been computed
        self.eigensystem()

        # first multiplication element-wise, since corresponding matrix only has non-zero entries on
        # the diagonal

        return (self.U * np.exp(self.eigenvals * t)) @ self.U_inv

    def to_indices(self, sequence: str) -> list[int]:
        """Convert a sequence into a list of indices.

        The letters in the sequence are converted to their corresponding indices in the alphabet of
        the model.

        Args:
            sequence: The sequence to be converted.

        Returns:
            The list of indices corresponding to the letters in the sequence.

        Raises:
            ValueError: If the sequence contains a letter that is not in the alphabet of the model.
        """
        try:
            result = [self.alphabet_dict[letter] for letter in sequence]
        except KeyError:
            raise ValueError("invalid sequence for the specified model")

        return result

    def to_sequence(self, evoseq: EvoSeq) -> str:
        """Construct the str representation of a sequence.

        Args:
            evoseq: The sequence to be converted as an `EvoSeq` object.

        Returns:
            The nucleotide or amino acid sequence depending on the model.
        """
        return "".join(self.alphabet[x._value] for x in evoseq)

__init__(model_type: str, model_name: str, **kwargs) -> None

Constructor for the SubstModel class.

Some models require additional parameters, e.g. 'K80' requires the parameter 'kappa' and 'GTR' requires 'abcdef' and 'f', see documentation for details. For a custom model, the keyword parameter 'filename' has to be specified with a path to a custom model in paml format. All of these parameters are passed as keyword arguments.

Parameters:

Name	Type	Description	Default
model_type	str	Available options are 'n' for nucleotide and 'a' for amino acid sequences.	required
model_name	str	Available are e.g. 'JC69' (nuc./a.a.), 'K80' (nuc.), 'GTR' (nuc.), 'WAG' (a.a.), 'JTT' (a.a.), 'BLOSUM62' (a.a.), 'LG' (a.a.), and 'DAYHOFF' (a.a.). For the option 'CUSTOM', the keyword parameter 'filename' has to be specified with a path to a custom model in paml format.	required

Raises:

Type	Description
ValueError	If an invalid model was specified.

Source code in asymmetree/seqevolve/substitution.py

def __init__(self, model_type: str, model_name: str, **kwargs) -> None:
    """Constructor for the SubstModel class.

    Some models require additional parameters, e.g. 'K80' requires the parameter 'kappa' and
    'GTR' requires 'abcdef' and 'f', see documentation for details. For a custom model, the
    keyword parameter 'filename' has to be specified with a path to a custom model in paml
    format. All of these parameters are passed as keyword arguments.

    Args:
        model_type: Available options are 'n' for nucleotide and 'a' for amino acid sequences.
        model_name: Available are e.g. 'JC69' (nuc./a.a.), 'K80' (nuc.), 'GTR' (nuc.), 'WAG'
            (a.a.), 'JTT' (a.a.), 'BLOSUM62' (a.a.), 'LG' (a.a.), and 'DAYHOFF' (a.a.). For the
            option 'CUSTOM', the keyword parameter 'filename' has to be specified with a path to
            a custom model in paml format.

    Raises:
        ValueError: If an invalid model was specified.
    """
    self.model_type = model_type.lower()
    self.model_name = model_name.upper()

    self._params = kwargs

    if (
        self.model_type in ("n", "nuc", "nucleotide")
        and self.model_name in SubstModel.nuc_models
    ):
        self.model_type = "n"
        self.alphabet = SubstModel.nucleotides

    elif self.model_type in ("a", "aa", "amino", "aminoacid", "protein") and (
        self.model_name in SubstModel.aa_models or self.model_name in EMPIRICAL_MODELS
    ):
        self.model_type = "a"
        self.alphabet = SubstModel.amino_acids

    else:
        raise ValueError(f"model '{model_type}', '{model_name}' is not available")

    self.alphabet_dict = {item: index for index, item in enumerate(self.alphabet)}

    self._load_exchangeability_and_freqs()
    self._build_rate_matrix()

eigensystem()

Compute the eigensystem of the rate matrix Q.

Source code in asymmetree/seqevolve/substitution.py

def eigensystem(self):
    """Compute the eigensystem of the rate matrix Q."""
    if not hasattr(self, "eigenvals"):
        self.eigenvals, self.U, self.U_inv = diagonalize(self.Q, self.freqs)

    return self.eigenvals, self.U, self.U_inv

transition_prob_matrix(t: float) -> np.ndarray

Calculate the transition probability matrix P(t).

The transition probability matrix P(t) is computed from the eigensystem of the rate matrix Q according to the formula P(t) = U x e^(Λ * t) x U^(-1) where U is the matrix of eigenvectors and Λ is the diagonal matrix of eigenvalues.

Parameters:

Name	Type	Description	Default
t	float	The time / evolutionary distance for which to compute the transition probability matrix.	required

Returns:

Type	Description
ndarray	The transition probability matrix P(t) as a numpy.ndarray.

Source code in asymmetree/seqevolve/substitution.py

def transition_prob_matrix(self, t: float) -> np.ndarray:
    """Calculate the transition probability matrix P(t).

    The transition probability matrix P(t) is computed from the eigensystem of the rate matrix Q
    according to the formula P(t) = U x e^(Λ * t) x U^(-1) where U is the matrix of eigenvectors
    and Λ is the diagonal matrix of eigenvalues.

    Args:
        t: The time / evolutionary distance for which to compute the transition probability
            matrix.

    Returns:
        The transition probability matrix P(t) as a `numpy.ndarray`.
    """
    # ensure that eigensystem has been computed
    self.eigensystem()

    # first multiplication element-wise, since corresponding matrix only has non-zero entries on
    # the diagonal

    return (self.U * np.exp(self.eigenvals * t)) @ self.U_inv

to_indices(sequence: str) -> list[int]

Convert a sequence into a list of indices.

The letters in the sequence are converted to their corresponding indices in the alphabet of the model.

Parameters:

Name	Type	Description	Default
sequence	str	The sequence to be converted.	required

Returns:

Type	Description
list[int]	The list of indices corresponding to the letters in the sequence.

Raises:

Type	Description
ValueError	If the sequence contains a letter that is not in the alphabet of the model.

Source code in asymmetree/seqevolve/substitution.py

def to_indices(self, sequence: str) -> list[int]:
    """Convert a sequence into a list of indices.

    The letters in the sequence are converted to their corresponding indices in the alphabet of
    the model.

    Args:
        sequence: The sequence to be converted.

    Returns:
        The list of indices corresponding to the letters in the sequence.

    Raises:
        ValueError: If the sequence contains a letter that is not in the alphabet of the model.
    """
    try:
        result = [self.alphabet_dict[letter] for letter in sequence]
    except KeyError:
        raise ValueError("invalid sequence for the specified model")

    return result

to_sequence(evoseq: EvoSeq) -> str

Construct the str representation of a sequence.

Parameters:

Name	Type	Description	Default
evoseq	EvoSeq	The sequence to be converted as an EvoSeq object.	required

Returns:

Type	Description
str	The nucleotide or amino acid sequence depending on the model.

Source code in asymmetree/seqevolve/substitution.py

def to_sequence(self, evoseq: EvoSeq) -> str:
    """Construct the str representation of a sequence.

    Args:
        evoseq: The sequence to be converted as an `EvoSeq` object.

    Returns:
        The nucleotide or amino acid sequence depending on the model.
    """
    return "".join(self.alphabet[x._value] for x in evoseq)

diagonalize(Q: np.ndarray, freqs: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Diagonalize a matrix.

Diagonalize the rate matrix Q using the stationary frequencies freqs. If Q is symmetric, it is diagonalized directly. If Q is not symmetric but the model is time-reversible, it is first transformed to a symmetric matrix using the stationary frequencies and then diagonalized.

Parameters:

Name	Type	Description	Default
Q	ndarray	Rate matrix to be diagonalized.	required
freqs	ndarray	The stationary frequencies π corresponding to the rate matrix Q.	required

Returns:

Type	Description
ndarray	A 1-dimensional array containing the eigenvalues.
ndarray	A 2-dimensional array containing the eigenvectors as columns. The i-th column corresponds to the i-th eigenvalue.
ndarray	A 2-dimensional array containing the inverse of the eigenvector matrix.

Source code in asymmetree/seqevolve/substitution.py

def diagonalize(Q: np.ndarray, freqs: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Diagonalize a matrix.

    Diagonalize the rate matrix Q using the stationary frequencies freqs. If Q is symmetric, it is
    diagonalized directly. If Q is not symmetric but the model is time-reversible, it is first
    transformed to a symmetric matrix using the stationary frequencies and then diagonalized.

    Args:
        Q: Rate matrix to be diagonalized.
        freqs: The stationary frequencies π corresponding to the rate matrix Q.

    Returns:
        A 1-dimensional array containing the eigenvalues.
        A 2-dimensional array containing the eigenvectors as columns. The i-th column corresponds to
            the i-th eigenvalue.
        A 2-dimensional array containing the inverse of the eigenvector matrix.
    """
    # matrix is already symmetric
    if np.allclose(Q, Q.T):
        eigenvals, U = linalg.eigh(Q)
        U_inv = linalg.inv(U)

    # matrix is not symmetric but model is time-reversible
    else:
        Phi = np.diag(np.sqrt(freqs))
        Phi_inv = linalg.inv(Phi)

        B = Phi @ Q @ Phi_inv

        eigenvals, R = linalg.eigh(B)

        U = Phi_inv @ R
        U_inv = linalg.inv(R) @ Phi

    return eigenvals, U, U_inv

heterogeneity

Heterogeneity models for sequence evolution.

Not all sites in a sequence evolve at the same rate. This module provides a class for modeling rate heterogeneity across sites, which is a common feature of sequence evolution. The class supports rate heterogeneity based on a Gamma distribution ('+Gamma'-model) as well as invariant sites ('+I'-model), see [1].

References

1. Z. Yang. Computational molecular evolution. Oxford series in ecology and evolution. Oxford University Press, 2006. ISBN 978-0-19-856699-1 978-0-19-856702-8.

HetModel

Heterogeneity model.

Supports rate heterogeneity based on a Gamma distribution ('+Gamma'-model) as well as invariant sites ('+I'-model), see [1].

References: 1. Z. Yang. Computational molecular evolution. Oxford series in ecology and evolution. Oxford University Press, 2006. ISBN 978-0-19-856699-1 978-0-19-856702-8.

Source code in asymmetree/seqevolve/heterogeneity.py

class HetModel:
    """Heterogeneity model.

    Supports rate heterogeneity based on a Gamma distribution ('+Gamma'-model) as well as invariant
    sites ('+I'-model), see [1].

    References:
    1. Z. Yang. Computational molecular evolution. Oxford series in ecology and evolution.
       Oxford University Press, 2006. ISBN 978-0-19-856699-1 978-0-19-856702-8.
    """

    def __init__(
        self,
        alpha: float,
        classes: int = 5,
        sitewise: bool = False,
        invariant: float = 0.0,
    ) -> None:
        """Constructor for the HetModel class.

        Args:
            alpha: Parameter of the Gamma distribution (with mean 1) from which the rate factors are
                drawn.
            classes: Number of classes such that sites in the same class share the same rate factor.
            sitewise: If True, ignore the 'classes' attribute and treat each site as its own class.
            invariant: Expected proportion of invariant sites, i.e., sites where no mutations
                happen.

        Raises:
            ValueError: If alpha is not a float > 0.0, if classes is not an int > 0, or if
                invariant is not a float in [0.0, 1.0].
        """
        if not isinstance(alpha, float) or alpha <= 0.0:
            raise ValueError("heterogeneity parameter alpha must be a float > 0.0")

        self._alpha = alpha

        if not isinstance(classes, int) or classes <= 0:
            raise ValueError("number of classes must be an int >0")

        self.classes = classes
        self.sitewise = sitewise

        if not isinstance(invariant, float) or invariant < 0.0 or invariant > 1.0:
            raise ValueError("proportion of invariant sites must be in [0.0, 1.0]")

        self._invariant = invariant

        if not self.sitewise:
            self._initialize_classes()

    def assign(self, sequence: EvoSeq, exclude_inherited: bool = True) -> None:
        """Assign rate classes and rate factors to the sites of a sequence.

        Args:
            sequence: The sequence to which the rate classes and rate factors should be assigned.
            exclude_inherited: If True, do not assign new classes and rates to inherited sites,
                i.e., only if the sequence corresponds to the root of a tree along which is
                simulated or if new sites are added as a result of an insertion event.
        """
        n = len(sequence)
        if exclude_inherited:
            n -= sequence.count_status(State.INHERITED)

        rate_classes, rate_factors = self._draw(n)

        pos = 0
        for site in sequence:
            if not exclude_inherited or site.status != State.INHERITED:
                site.rate_class = (
                    rate_classes[pos] if rate_classes[pos] is not None else site.site_id
                )
                site.rate_factor = rate_factors[pos]
                pos += 1

    def _initialize_classes(self) -> None:
        """Initialize the rate factors for the classes."""
        if self.classes == 1:
            self._class_rates = [1.0]
        else:
            self._class_rates = np.random.gamma(
                self._alpha, scale=1 / self._alpha, size=self.classes
            )

    def _draw(self, n) -> tuple[list[int | str | None], np.ndarray]:
        """Draw rate classes and rate factors for n sites.

        Args:
            n: The number of sites for which to draw the rate classes and rate factors.

        Returns:
            A list of rate classes.
            An array of rate factors for the `n` sites.
        """
        # mode 1: sitewise heterogeneity
        if self.sitewise:
            rate_classes = [None for _ in range(n)]
            rate_factors = np.random.gamma(self._alpha, scale=1 / self._alpha, size=n)

        # mode 2: one or multiple  classes
        else:
            rate_classes = np.random.randint(self.classes, size=n).tolist()
            rate_factors = np.asarray([self._class_rates[c] for c in rate_classes])

        if self._invariant:
            drawn_variable = np.random.random(n) > self._invariant

            for i in range(n):
                if not drawn_variable[i]:
                    rate_classes[i] = "invariant"

            rate_factors *= drawn_variable

        return rate_classes, rate_factors

__init__(alpha: float, classes: int = 5, sitewise: bool = False, invariant: float = 0.0) -> None

Constructor for the HetModel class.

Parameters:

Name	Type	Description	Default
alpha	float	Parameter of the Gamma distribution (with mean 1) from which the rate factors are drawn.	required
classes	int	Number of classes such that sites in the same class share the same rate factor.	5
sitewise	bool	If True, ignore the 'classes' attribute and treat each site as its own class.	False
invariant	float	Expected proportion of invariant sites, i.e., sites where no mutations happen.	0.0

Raises:

Type	Description
ValueError	If alpha is not a float > 0.0, if classes is not an int > 0, or if invariant is not a float in [0.0, 1.0].

Source code in asymmetree/seqevolve/heterogeneity.py

def __init__(
    self,
    alpha: float,
    classes: int = 5,
    sitewise: bool = False,
    invariant: float = 0.0,
) -> None:
    """Constructor for the HetModel class.

    Args:
        alpha: Parameter of the Gamma distribution (with mean 1) from which the rate factors are
            drawn.
        classes: Number of classes such that sites in the same class share the same rate factor.
        sitewise: If True, ignore the 'classes' attribute and treat each site as its own class.
        invariant: Expected proportion of invariant sites, i.e., sites where no mutations
            happen.

    Raises:
        ValueError: If alpha is not a float > 0.0, if classes is not an int > 0, or if
            invariant is not a float in [0.0, 1.0].
    """
    if not isinstance(alpha, float) or alpha <= 0.0:
        raise ValueError("heterogeneity parameter alpha must be a float > 0.0")

    self._alpha = alpha

    if not isinstance(classes, int) or classes <= 0:
        raise ValueError("number of classes must be an int >0")

    self.classes = classes
    self.sitewise = sitewise

    if not isinstance(invariant, float) or invariant < 0.0 or invariant > 1.0:
        raise ValueError("proportion of invariant sites must be in [0.0, 1.0]")

    self._invariant = invariant

    if not self.sitewise:
        self._initialize_classes()

assign(sequence: EvoSeq, exclude_inherited: bool = True) -> None

Assign rate classes and rate factors to the sites of a sequence.

Parameters:

Name	Type	Description	Default
sequence	EvoSeq	The sequence to which the rate classes and rate factors should be assigned.	required
exclude_inherited	bool	If True, do not assign new classes and rates to inherited sites, i.e., only if the sequence corresponds to the root of a tree along which is simulated or if new sites are added as a result of an insertion event.	True

Source code in asymmetree/seqevolve/heterogeneity.py

def assign(self, sequence: EvoSeq, exclude_inherited: bool = True) -> None:
    """Assign rate classes and rate factors to the sites of a sequence.

    Args:
        sequence: The sequence to which the rate classes and rate factors should be assigned.
        exclude_inherited: If True, do not assign new classes and rates to inherited sites,
            i.e., only if the sequence corresponds to the root of a tree along which is
            simulated or if new sites are added as a result of an insertion event.
    """
    n = len(sequence)
    if exclude_inherited:
        n -= sequence.count_status(State.INHERITED)

    rate_classes, rate_factors = self._draw(n)

    pos = 0
    for site in sequence:
        if not exclude_inherited or site.status != State.INHERITED:
            site.rate_class = (
                rate_classes[pos] if rate_classes[pos] is not None else site.site_id
            )
            site.rate_factor = rate_factors[pos]
            pos += 1

indels

Indel models for sequence evolution.

Indel events (insertions and deletions) are an important aspect of sequence evolution. This module provides a class for modeling indel events, which can be used to simulate sequence evolution along a tree. The class supports different distributions for the length of the indels.

References

1. Z. Yang. Computational molecular evolution. Oxford series in ecology and evolution. Oxford University Press, 2006. ISBN 978-0-19-856699-1 978-0-19-856702-8.
2. R. A. Cartwright. DNA assembly with gaps (Dawg): Simulating sequence evolution. In: Bioinformatics, 21(Suppl 3):iii31-iii38, November 2005. doi:10.1093/bioinformatics/bti1200.

IndelModel

Indel model

Source code in asymmetree/seqevolve/indels.py

class IndelModel:
    """Indel model"""

    def __init__(
        self,
        insertion_rate: float,
        deletion_rate: float,
        length_distr: tuple = ("zipf", 1.821),  # (Chang and Benner 2004)
        min_length: int = 1,
        max_length: int | None = None,
        **kwargs,
    ) -> None:
        """Constructor for the IndelModel class.

        Args:
            insertion_rate: Insertion rate w.r.t. units of evolutionary distance.
            deletion_rate: Deletion rate w.r.t. units of evolutionary distance.
            length_distr: Distribution of the length of the indels, see documentation for available
                options. The default is a Zipf distribution as observed for empirical data in [1].
            min_length: The minimal length of an indel, default is 1.
            max_length: The maximal length of an indel. The default is None in which case the
                distribution is not truncated (on the upper side).

        Raises:
            ValueError: If insertion_rate or deletion_rate is negative.

        References:
            1. M. S. Chang and S. A. Benner. Empirical Analysis of Protein Insertions and
               Deletions Determining Parameters for the Correct Placement of Gaps in Protein
               Sequence Alignments. In: Journal of Molecular Biology, 341(2):617-631, August 2004.
               doi:10.1016/j.jmb.2004.05.045.
        """
        if insertion_rate < 0.0 or deletion_rate < 0.0:
            raise ValueError("insertion and deletion rates must be non-negative")

        self._ins_rate = insertion_rate  # insertion rate per site
        self._del_rate = deletion_rate  # deletion rate per site

        if isinstance(length_distr, Sampler):
            self.sampler = length_distr
        else:
            self.sampler = Sampler(
                length_distr, minimum=min_length, maximum=max_length, discrete=True
            )

    def get_rates(self, seq_length: int) -> tuple[float, float]:
        """Return the current insertion and deletion rate.

        Computes the current insertion and deletion rate according to the model and the length of
        the sequence.

        Args:
            seq_length: The length of the sequence.

        Returns:
            The total insertion rate for the given sequence length, see 1.
            The total deletion rate for the given sequence length, see 1.

        References:
            1. R. A. Cartwright. DNA assembly with gaps (Dawg): Simulating sequence evolution.
               In: Bioinformatics, 21(Suppl 3):iii31-iii38, November 2005.
               doi:10.1093/bioinformatics/bti1200.
        """
        # expected value may be infinite for zipf distribution, in which case the minimal length is
        # used, see Sampler class for details
        return (
            (seq_length + 1) * self._ins_rate,
            (seq_length + self.sampler._exp_val - 1) * self._del_rate,
        )

    def draw_length(self) -> int:
        """Draw the length for an indel.

        Returns:
            The drawn length for the indel.
        """
        return self.sampler()

__init__(insertion_rate: float, deletion_rate: float, length_distr: tuple = ('zipf', 1.821), min_length: int = 1, max_length: int | None = None, **kwargs) -> None

Constructor for the IndelModel class.

Parameters:

Name	Type	Description	Default
insertion_rate	float	Insertion rate w.r.t. units of evolutionary distance.	required
deletion_rate	float	Deletion rate w.r.t. units of evolutionary distance.	required
length_distr	tuple	Distribution of the length of the indels, see documentation for available options. The default is a Zipf distribution as observed for empirical data in [1].	('zipf', 1.821)
min_length	int	The minimal length of an indel, default is 1.	1
max_length	int | None	The maximal length of an indel. The default is None in which case the distribution is not truncated (on the upper side).	None

Raises:

Type	Description
ValueError	If insertion_rate or deletion_rate is negative.

References

1. M. S. Chang and S. A. Benner. Empirical Analysis of Protein Insertions and Deletions Determining Parameters for the Correct Placement of Gaps in Protein Sequence Alignments. In: Journal of Molecular Biology, 341(2):617-631, August 2004. doi:10.1016/j.jmb.2004.05.045.

Source code in asymmetree/seqevolve/indels.py

def __init__(
    self,
    insertion_rate: float,
    deletion_rate: float,
    length_distr: tuple = ("zipf", 1.821),  # (Chang and Benner 2004)
    min_length: int = 1,
    max_length: int | None = None,
    **kwargs,
) -> None:
    """Constructor for the IndelModel class.

    Args:
        insertion_rate: Insertion rate w.r.t. units of evolutionary distance.
        deletion_rate: Deletion rate w.r.t. units of evolutionary distance.
        length_distr: Distribution of the length of the indels, see documentation for available
            options. The default is a Zipf distribution as observed for empirical data in [1].
        min_length: The minimal length of an indel, default is 1.
        max_length: The maximal length of an indel. The default is None in which case the
            distribution is not truncated (on the upper side).

    Raises:
        ValueError: If insertion_rate or deletion_rate is negative.

    References:
        1. M. S. Chang and S. A. Benner. Empirical Analysis of Protein Insertions and
           Deletions Determining Parameters for the Correct Placement of Gaps in Protein
           Sequence Alignments. In: Journal of Molecular Biology, 341(2):617-631, August 2004.
           doi:10.1016/j.jmb.2004.05.045.
    """
    if insertion_rate < 0.0 or deletion_rate < 0.0:
        raise ValueError("insertion and deletion rates must be non-negative")

    self._ins_rate = insertion_rate  # insertion rate per site
    self._del_rate = deletion_rate  # deletion rate per site

    if isinstance(length_distr, Sampler):
        self.sampler = length_distr
    else:
        self.sampler = Sampler(
            length_distr, minimum=min_length, maximum=max_length, discrete=True
        )

get_rates(seq_length: int) -> tuple[float, float]

Return the current insertion and deletion rate.

Computes the current insertion and deletion rate according to the model and the length of the sequence.

Parameters:

Name	Type	Description	Default
seq_length	int	The length of the sequence.	required

Returns:

Type	Description
float	The total insertion rate for the given sequence length, see 1.
float	The total deletion rate for the given sequence length, see 1.

References

1. R. A. Cartwright. DNA assembly with gaps (Dawg): Simulating sequence evolution. In: Bioinformatics, 21(Suppl 3):iii31-iii38, November 2005. doi:10.1093/bioinformatics/bti1200.

Source code in asymmetree/seqevolve/indels.py

def get_rates(self, seq_length: int) -> tuple[float, float]:
    """Return the current insertion and deletion rate.

    Computes the current insertion and deletion rate according to the model and the length of
    the sequence.

    Args:
        seq_length: The length of the sequence.

    Returns:
        The total insertion rate for the given sequence length, see 1.
        The total deletion rate for the given sequence length, see 1.

    References:
        1. R. A. Cartwright. DNA assembly with gaps (Dawg): Simulating sequence evolution.
           In: Bioinformatics, 21(Suppl 3):iii31-iii38, November 2005.
           doi:10.1093/bioinformatics/bti1200.
    """
    # expected value may be infinite for zipf distribution, in which case the minimal length is
    # used, see Sampler class for details
    return (
        (seq_length + 1) * self._ins_rate,
        (seq_length + self.sampler._exp_val - 1) * self._del_rate,
    )

draw_length() -> int

Draw the length for an indel.

Returns:

Type	Description
int	The drawn length for the indel.

Source code in asymmetree/seqevolve/indels.py

def draw_length(self) -> int:
    """Draw the length for an indel.

    Returns:
        The drawn length for the indel.
    """
    return self.sampler()

alignment

Module for constructing the true alignment of simulated sequences.

In the simulated sequences, the true history of substitutions and indels is known. This module uses this information to construct the true alignment of the sequences, which can for instance be used for evaluating the performance of alignment methods.

AlignmentBuilder

Construction of the true alignment of simulated sequences.

Source code in asymmetree/seqevolve/alignment.py

class AlignmentBuilder:
    """Construction of the true alignment of simulated sequences."""

    def __init__(
        self,
        tree: Tree,
        sequence_dict: dict[TreeNode, EvoSeq],
        alphabet: str,
        include_inner: bool = True,
    ) -> None:
        """Constructor for the AlignmentBuilder class.

        Args:
            tree: The tree along which the sequences where simulated.
            sequence_dict: A dict containing the TreeNode instances in the tree as keys and the
                simulated sequences as values.
            alphabet: The alphabet of the substitution model that was used.
            include_inner: If True, the alignment also contains the sequences of the inner nodes of
                the tree.
        """
        self.tree = tree
        self.sequence_dict = sequence_dict
        self.include_inner = include_inner

        # add gap character to the alphabet if not already present
        self.alphabet = alphabet if alphabet[-1] == "-" else alphabet + "-"

    def build(self) -> dict[TreeNode, str]:
        """Build the true alignment.

        Returns:
            A `dict` containing the `TreeNode` instances in the tree as keys and the str sequences
                as values that include the necessary gaps.
        """
        self._get_preorder()
        self._sort_sites()
        self._alignment_matrix()

        return self._sequences()

    def _get_preorder(self) -> None:
        """Get the preorder of the tree."""
        self.nodes = [v for v in self.tree.preorder() if self.include_inner or not v.children]
        self.node_index = {item: index for index, item in enumerate(self.nodes)}

    def _sort_sites(self) -> None:
        """Sort the sites of the sequences in a way that they can be correctly aligned."""
        G = nx.DiGraph()

        for node in self.nodes:
            seq = self.sequence_dict[node]

            # handle case |seq| = 1 correctly
            if len(seq) == 1:
                G.add_node(seq[0].site_id)

            # add nodes and edges to auxiliary graph
            for first, second in seq.element_pairs():
                G.add_edge(first.site_id, second.site_id)

        self.sites = list(nx.topological_sort(G))

    def _alignment_matrix(self) -> None:
        """Construct the alignment matrix."""
        self.alignment = np.zeros((len(self.nodes), len(self.sites)), dtype=np.int8)
        positions = [self.sequence_dict[v]._first for v in self.nodes]

        for j in range(len(self.sites)):
            for i in range(len(positions)):
                if positions[i] and self.sites[j] == positions[i].site_id:
                    self.alignment[i, j] = positions[i]._value
                    positions[i] = positions[i]._next
                else:
                    self.alignment[i, j] = -1

    def _sequences(self) -> dict[TreeNode, str]:
        """Convert the alignment matrix to a dictionary of sequences.

        Returns:
            A `dict` containing the `TreeNode` instances in the tree as keys and the str sequences
                as values that include the necessary gaps.
        """
        result = {}

        for i in range(len(self.nodes)):
            sequence = "".join(self.alphabet[x] for x in self.alignment[i, :])
            result[self.nodes[i]] = sequence

        return result

__init__(tree: Tree, sequence_dict: dict[TreeNode, EvoSeq], alphabet: str, include_inner: bool = True) -> None

Constructor for the AlignmentBuilder class.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree along which the sequences where simulated.	required
sequence_dict	dict[TreeNode, EvoSeq]	A dict containing the TreeNode instances in the tree as keys and the simulated sequences as values.	required
alphabet	str	The alphabet of the substitution model that was used.	required
include_inner	bool	If True, the alignment also contains the sequences of the inner nodes of the tree.	True

Source code in asymmetree/seqevolve/alignment.py

def __init__(
    self,
    tree: Tree,
    sequence_dict: dict[TreeNode, EvoSeq],
    alphabet: str,
    include_inner: bool = True,
) -> None:
    """Constructor for the AlignmentBuilder class.

    Args:
        tree: The tree along which the sequences where simulated.
        sequence_dict: A dict containing the TreeNode instances in the tree as keys and the
            simulated sequences as values.
        alphabet: The alphabet of the substitution model that was used.
        include_inner: If True, the alignment also contains the sequences of the inner nodes of
            the tree.
    """
    self.tree = tree
    self.sequence_dict = sequence_dict
    self.include_inner = include_inner

    # add gap character to the alphabet if not already present
    self.alphabet = alphabet if alphabet[-1] == "-" else alphabet + "-"

build() -> dict[TreeNode, str]

Build the true alignment.

Returns:

Type	Description
dict[TreeNode, str]	A dict containing the TreeNode instances in the tree as keys and the str sequences as values that include the necessary gaps.

Source code in asymmetree/seqevolve/alignment.py

def build(self) -> dict[TreeNode, str]:
    """Build the true alignment.

    Returns:
        A `dict` containing the `TreeNode` instances in the tree as keys and the str sequences
            as values that include the necessary gaps.
    """
    self._get_preorder()
    self._sort_sites()
    self._alignment_matrix()

    return self._sequences()