treeevolve

Simulation of species and gene trees.

species

Simulation of dated species trees.

species_tree_n(n: int, model: str = 'yule', innovation: bool = False, planted: bool = True, remove_extinct: bool = False, birth_rate: float | None = None, death_rate: float | None = None, episodes: list[EpisodeType] | None = None, contraction_probability: float = 0.0, contraction_proportion: float = 0.0, contraction_bias: str | None = None, bias_strength: float = 1.0, **kwargs: object) -> Tree

Simulates a species tree S with n leaves.

The tree is sampled under the Yule, BDP, or EBDP model conditioned on the number n of surviving species [1, 2, 3].

Parameters:

Name	Type	Description	Default
n	int	Number of leaves in the resulting tree that correspond to extant species.	required
model	str	Simulation model to be applied, the default is 'yule', see [1]. Other available models are birth-death process ('BDP', [2]), and episodic birth-death process ('EBDP', [3]).	'yule'
innovation	bool	If True, use the innovation model, see [4], to sample a lineage for the next speciation event. Only available for the Yule model. The default is False, in which case the lineage is chosen uniformly at random among the currently existing lineages. The innovation model is not available for the BDP and EBDP models.	False
planted	bool	Add a planted root that has the canonical root as its single neighbor.	True
remove_extinct	bool	Remove all branches that lead to extinctions, only relevant for some models.	False
birth_rate	float | None	The birth rate for models such as 'yule' and 'BDP'. The default is None, in which case the birth rate is set to 1.0 unless episodes are specified.	None
death_rate	float | None	The death rate for models such as 'BDP'. The default is None, in which case the death rate is set to 0.0 unless episodes are specified.	None
episodes	list[EpisodeType] | None	The episodes for the model 'EBDP'.	None
contraction_probability	float	Probability that an inner edge is contracted. The default is 0.0, in which case the tree is binary. Only one of this parameter and 'contraction_proportion' may be non-zero.	0.0
contraction_proportion	float	The proportion of inner edges to be contracted. The default is 0.0, in which case the tree is binary. Only one of this parameter and 'contraction_probability' may be non-zero.	0.0
contraction_bias	str | None	Specifies whether shorter edges, i.e., with a smaller difference t of the time stamps, have a higher probability to be contracted. Only relevant if 'contraction_proportion' > 0.0. The default is None, in which case all edges have the same probability to be contracted. The options 'inverse' and 'exponential' mean that an edge is sampled weighted by 1/t or e^(-t), respectively.	None
bias_strength	float	Intensity factor for preferring shorter edges to be contracted. The default is 1.0.	1.0

Raises:

Type	Description
ValueError	If unknown or invalid parameter values are passed.
ValueError	If the innovation model is selected for a model that does not support it.

Returns:

Type	Description
Tree	The simulated species tree.

References

1. G. U. Yule. A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F. R. S. In: Phil. Trans. R. Soc. Lond. B, 1924, 213, 21-87. doi:10.1098/rstb.1925.0002.
2. D. G. Kendall. On the Generalized "Birth-and-Death" Process. In: Ann. Math. Statist. 1948, 19, 1-15. doi:10.1214/aoms/1177730285.
3. T. Stadler. Simulating trees with a fixed number of extant species. In: Syst. Biol. 2011, 60, 676-684. doi:10.1093/sysbio/syr029.
4. S. Keller-Schmidt, K. Klemm. A model of macroevolution as a branching process based on innovations. In: Adv. Complex Syst.,2012, 15, 1250043. doi:10.1142/S0219525912500439.

Source code in asymmetree/treeevolve/species.py

def species_tree_n(
    n: int,
    model: str = "yule",
    innovation: bool = False,
    planted: bool = True,
    remove_extinct: bool = False,
    birth_rate: float | None = None,
    death_rate: float | None = None,
    episodes: list[EpisodeType] | None = None,
    contraction_probability: float = 0.0,
    contraction_proportion: float = 0.0,
    contraction_bias: str | None = None,
    bias_strength: float = 1.0,
    **kwargs: object,
) -> Tree:
    """Simulates a species tree S with n leaves.

    The tree is sampled under the Yule, BDP, or EBDP model conditioned on the number n of surviving
    species [1, 2, 3].

    Args:
        n: Number of leaves in the resulting tree that correspond to extant species.
        model: Simulation model to be applied, the default is 'yule', see [1]. Other available
            models are birth-death process ('BDP', [2]), and episodic birth-death process ('EBDP',
            [3]).
        innovation: If True, use the innovation model, see [4], to sample a lineage for the next
            speciation event. Only available for the Yule model. The default is False, in which case
            the lineage is chosen uniformly at random among the currently existing lineages.
            The innovation model is not available for the BDP and EBDP models.
        planted: Add a planted root that has the canonical root as its single neighbor.
        remove_extinct: Remove all branches that lead to extinctions, only relevant for some models.
        birth_rate: The birth rate for models such as 'yule' and 'BDP'. The default is None, in
            which case the birth rate is set to 1.0 unless episodes are specified.
        death_rate: The death rate for models such as 'BDP'. The default is None, in which case the
            death rate is set to 0.0 unless episodes are specified.
        episodes: The episodes for the model 'EBDP'.
        contraction_probability: Probability that an inner edge is contracted. The default is 0.0,
            in which case the tree is binary. Only one of this parameter and
            'contraction_proportion' may be non-zero.
        contraction_proportion: The proportion of inner edges to be contracted. The default is 0.0,
            in which case the tree is binary. Only one of this parameter and
            'contraction_probability' may be non-zero.
        contraction_bias: Specifies whether shorter edges, i.e., with a smaller difference t of the
            time stamps, have a higher probability to be contracted. Only relevant if
            'contraction_proportion' > 0.0. The default is None, in which case all edges have the
            same probability to be contracted. The options 'inverse' and 'exponential' mean that an
            edge is sampled weighted by 1/t or e^(-t), respectively.
        bias_strength: Intensity factor for preferring shorter edges to be contracted. The default
            is 1.0.

    Raises:
        ValueError: If unknown or invalid parameter values are passed.
        ValueError: If the innovation model is selected for a model that does not support it.

    Returns:
        The simulated species tree.

    References:
        1. G. U. Yule. A mathematical theory of evolution, based on the conclusions of Dr. J.
           C. Willis, F. R. S. In: Phil. Trans. R. Soc. Lond. B, 1924, 213, 21-87.
           doi:10.1098/rstb.1925.0002.
        2. D. G. Kendall. On the Generalized "Birth-and-Death" Process. In: Ann. Math. Statist.
           1948, 19, 1-15. doi:10.1214/aoms/1177730285.
        3. T. Stadler. Simulating trees with a fixed number of extant species. In: Syst. Biol.
           2011, 60, 676-684. doi:10.1093/sysbio/syr029.
        4. S. Keller-Schmidt, K. Klemm. A model of macroevolution as a branching process based
           on innovations. In: Adv. Complex Syst.,2012, 15, 1250043. doi:10.1142/S0219525912500439.
    """
    # parameter checking
    if not isinstance(n, int) or n < 0:
        raise ValueError("n must be an int >=0")
    elif n == 0:
        return Tree(None)

    if not isinstance(model, str):
        raise ValueError("model must be of type 'str'")

    # choice of the main simulation algorithm
    if model.lower() == "yule":
        tree = _yule_n(n, birth_rate, innovation)
    elif model.upper() == "BDP":
        if innovation:
            raise ValueError("innovation model is not available for BDP")

        tree = _BDP_n(n, birth_rate, death_rate)
    elif model.upper() == "EBDP":
        if innovation:
            raise ValueError("innovation model is not available for EBDP")

        tree = _EBDP_n(n, birth_rate, death_rate, episodes)
    else:
        raise ValueError(f"model '{model}' is not available")

    # remove extinct branches for models that include losses
    if remove_extinct and model.upper() in ("BDP", "EBDP"):
        delete_losses_and_contract(tree, inplace=True)

    # remove planted edge for models that are planted by construction
    if not planted and model.upper() in ("YULE", "BDP", "EBDP"):
        remove_planted_root(tree, inplace=True)

    # make tree non_binary by random contraction of edges
    nonbinary(
        tree,
        contraction_probability=contraction_probability,
        contraction_proportion=contraction_proportion,
        contraction_bias=contraction_bias,
        bias_strength=bias_strength,
        inplace=True,
    )

    # assign the distance attribute to all vertices
    distance_from_timing(tree)

    return tree

species_tree_age(age: float, model: str = 'yule', innovation: bool = False, birth_rate: float | None = None, death_rate: float | None = None, episodes: list[EpisodeType] | None = None, contraction_probability: float = 0.0, contraction_proportion: float = 0.0, contraction_bias: str | None = None, bias_strength: float = 1.0, **kwargs: object) -> Tree

Simulates a (planted) species tree S of the specified age.

Parameters:

Name	Type	Description	Default
age	float	Simulation time, i.e., the time span from the root of the tree to the leaves that correspond to extant species.	required
model	str	Simulation model to be applied, the default is 'yule', see [1]. Other available models are birth-death process ('BDP', [2]), and episodic birth-death process ('EBDP', [3]).	'yule'
innovation	bool	If True, use the innovation model, see [4], to sample a lineage for the next speciation event. The default is False, in which case the lineage is chosen uniformly at random among the currently existing lineages.	False
birth_rate	float | None	The birth rate for models such as 'yule' and 'BDP'. The default is None, in which case the birth rate is set to 1.0 unless episodes are specified.	None
death_rate	float | None	The death rate for models such as 'BDP'. The default is None, in which case the death rate is set to 0.0 unless episodes are specified.	None
episodes	list[EpisodeType] | None	The episodes for the model 'EBDP'.	None
contraction_probability	float	Probability that an inner edge is contracted. The default is 0.0, in which case the tree is binary. Only one of this parameter and 'contraction_proportion' may be non-zero.	0.0
contraction_proportion	float	The proportion of inner edges to be contracted. The default is 0.0, in which case the tree is binary. Only one of this parameter and 'contraction_probability' may be non-zero.	0.0
contraction_bias	str | None	Specifies whether shorter edges, i.e., with a smaller difference t of the time stamps, have a higher probability to be contracted. Only relevant if 'contraction_proportion' > 0.0. The default is None, in which case all edges have the same probability to be contracted. The options 'inverse' and 'exponential' mean that an edge is sampled weighted by 1/(a * t) or e^(-a * t), respectively, where a is a user-defined factor.	None
bias_strength	float	Intensity factor for preferring shorter edges to be contracted. The default is 1.0.	1.0

Raises:

Type	Description
ValueError	If unknown or invalid parameter values are passed.

Returns:

Type	Description
Tree	The simulated species tree.

References

1. G. U. Yule. A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F. R. S. In: Phil. Trans. R. Soc. Lond. B, 1924, 213, 21-87. doi:10.1098/rstb.1925.0002.
2. D. G. Kendall. On the Generalized "Birth-and-Death" Process. In: Ann. Math. Statist. 1948, 19, 1-15. doi:10.1214/aoms/1177730285.
3. T. Stadler. Simulating trees with a fixed number of extant species. In: Syst. Biol. 2011, 60, 676-684. doi:10.1093/sysbio/syr029.
4. S. Keller-Schmidt, K. Klemm. A model of macroevolution as a branching process based on innovations. In: Adv. Complex Syst.,2012, 15, 1250043. doi:10.1142/S0219525912500439.

Source code in asymmetree/treeevolve/species.py

def species_tree_age(
    age: float,
    model: str = "yule",
    innovation: bool = False,
    birth_rate: float | None = None,
    death_rate: float | None = None,
    episodes: list[EpisodeType] | None = None,
    contraction_probability: float = 0.0,
    contraction_proportion: float = 0.0,
    contraction_bias: str | None = None,
    bias_strength: float = 1.0,
    **kwargs: object,
) -> Tree:
    """Simulates a (planted) species tree S of the specified age.

    Args:
        age: Simulation time, i.e., the time span from the root of the tree to the leaves that
            correspond to extant species.
        model: Simulation model to be applied, the default is 'yule', see [1]. Other available
            models are birth-death process ('BDP', [2]), and episodic birth-death process ('EBDP',
            [3]).
        innovation: If True, use the innovation model, see [4], to sample a lineage for the next
            speciation event. The default is False, in which case the lineage is chosen uniformly at
            random among the currently existing lineages.
        birth_rate: The birth rate for models such as 'yule' and 'BDP'. The default is None, in
            which case the birth rate is set to 1.0 unless episodes are specified.
        death_rate: The death rate for models such as 'BDP'. The default is None, in which case the
            death rate is set to 0.0 unless episodes are specified.
        episodes: The episodes for the model 'EBDP'.
        contraction_probability: Probability that an inner edge is contracted. The default is 0.0,
            in which case the tree is binary. Only one of this parameter and
            'contraction_proportion' may be non-zero.
        contraction_proportion: The proportion of inner edges to be contracted. The default is 0.0,
            in which case the tree is binary. Only one of this parameter and
            'contraction_probability' may be non-zero.
        contraction_bias: Specifies whether shorter edges, i.e., with a smaller difference t of the
            time stamps, have a higher probability to be contracted. Only relevant if
            'contraction_proportion' > 0.0. The default is None, in which case all edges have the
            same probability to be contracted. The options 'inverse' and 'exponential' mean that an
            edge is sampled weighted by 1/(a * t) or e^(-a * t), respectively, where a is a
            user-defined factor.
        bias_strength: Intensity factor for preferring shorter edges to be contracted. The default
            is 1.0.

    Raises:
        ValueError: If unknown or invalid parameter values are passed.

    Returns:
        The simulated species tree.

    References:
        1. G. U. Yule. A mathematical theory of evolution, based on the conclusions of Dr. J.
           C. Willis, F. R. S. In: Phil. Trans. R. Soc. Lond. B, 1924, 213, 21-87.
           doi:10.1098/rstb.1925.0002.
        2. D. G. Kendall. On the Generalized "Birth-and-Death" Process. In: Ann. Math. Statist.
           1948, 19, 1-15. doi:10.1214/aoms/1177730285.
        3. T. Stadler. Simulating trees with a fixed number of extant species. In: Syst. Biol.
           2011, 60, 676-684. doi:10.1093/sysbio/syr029.
        4. S. Keller-Schmidt, K. Klemm. A model of macroevolution as a branching process based
           on innovations. In: Adv. Complex Syst.,2012, 15, 1250043. doi:10.1142/S0219525912500439.
    """
    # parameter checking
    if not isinstance(age, (float, int)) or age <= 0.0:
        raise ValueError("age must be a number >0")
    elif isinstance(age, int):
        age = float(age)

    if not isinstance(model, str):
        raise ValueError("model must be of type 'str'")

    # main simulation algorithm
    if model.lower() == "yule":
        tree = _yule_age(age, birth_rate, innovation)
    elif model.upper() == "BDP":
        tree = _BDP_age(age, innovation, birth_rate, death_rate)
    elif model.upper() == "EBDP":
        tree = _EBDP_age(age, innovation, birth_rate, death_rate, episodes)
    else:
        raise ValueError(f"model '{model}' is not available")

    # make tree non_binary by random contraction of edges
    nonbinary(
        tree,
        contraction_probability=contraction_probability,
        contraction_proportion=contraction_proportion,
        contraction_bias=contraction_bias,
        bias_strength=bias_strength,
        inplace=True,
    )

    # assign the distance attribute to all vertices
    distance_from_timing(tree)

    return tree

species_tree_n_age(n: int, age: float, model: str = 'yule', innovation: bool = False, birth_rate: float | None = None, death_rate: float | None = None, contraction_probability: float = 0.0, contraction_proportion: float = 0.0, contraction_bias: str | None = None, bias_strength: float = 1.0, **kwargs: object) -> Tree

Simulate a (planted) species tree S with n leaves and of the specified age.

The tree is sampled under the Yule or BDP model conditioned on the number n of surviving species and with a given age [1]. The resulting tree does not contain loss leaves even if the specified death rate is > 0 as a consequence of the tree sampling method.

Parameters:

Name	Type	Description	Default
n	int	Number of leaves in the resulting tree that correspond to extant species.	required
age	float	Simulation time, i.e., the time span from the root of the tree to the leaves that correspond to extant species.	required
model	str	Simulation model to be applied, the default is 'yule', see [2]. Alternatively the birth-death process ('BDP', [3]) is available.	'yule'
innovation	bool	If True, use the innovation model, see [4], to sample a lineage for the next speciation event. The default is False, in which case the lineage is chosen uniformly at random among the currently existing lineages.	False
birth_rate	float | None	The birth rate for models such as 'yule' and 'BDP'. The default is None, in which case the birth rate is set to 1.0 unless episodes are specified.	None
death_rate	float | None	The death rate for models such as 'BDP'. The default is None, in which case the death rate is set to 0.0 unless episodes are specified.	None
contraction_probability	float	Probability that an inner edge is contracted. The default is 0.0, in which case the tree is binary. Only one of this parameter and 'contraction_proportion' may be non-zero.	0.0
contraction_proportion	float	The proportion of inner edges to be contracted. The default is 0.0, in which case the tree is binary. Only one of this parameter and 'contraction_probability' may be non-zero.	0.0
contraction_bias	str | None	Specifies whether shorter edges, i.e., with a smaller difference t of the time stamps, have a higher probability to be contracted. Only relevant if 'contraction_proportion' > 0.0. The default is None, in which case all edges have the same probability to be contracted. The options 'inverse' and 'exponential' mean that an edge is sampled weighted by 1/(a * t) or e^(-a * t), respectively, where a is a user-defined factor.	None
bias_strength	float	Intensity factor for preferring shorter edges to be contracted. The default is 1.0.	1.0

Raises:

Type	Description
ValueError	If unknown or invalid parameter values are passed.

Returns:

Type	Description
Tree	The simulated species tree.

References

1. G. U. Yule. A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F. R. S. In: Phil. Trans. R. Soc. Lond. B, 1924, 213, 21-87. doi:10.1098/rstb.1925.0002.
2. D. G. Kendall. On the Generalized "Birth-and-Death" Process. In: Ann. Math. Statist. 1948, 19, 1-15. doi:10.1214/aoms/1177730285.
3. T. Stadler. Simulating trees with a fixed number of extant species. In: Syst. Biol. 2011, 60, 676-684. doi:10.1093/sysbio/syr029.
4. S. Keller-Schmidt, K. Klemm. A model of macroevolution as a branching process based on innovations. In: Adv. Complex Syst.,2012, 15, 1250043. doi:10.1142/S0219525912500439.

Source code in asymmetree/treeevolve/species.py

def species_tree_n_age(
    n: int,
    age: float,
    model: str = "yule",
    innovation: bool = False,
    birth_rate: float | None = None,
    death_rate: float | None = None,
    contraction_probability: float = 0.0,
    contraction_proportion: float = 0.0,
    contraction_bias: str | None = None,
    bias_strength: float = 1.0,
    **kwargs: object,
) -> Tree:
    """Simulate a (planted) species tree S with n leaves and of the specified age.

    The tree is sampled under the Yule or BDP model conditioned on the number n of surviving species
    and with a given age [1]. The resulting tree does not contain loss leaves even if the specified
    death rate is > 0 as a consequence of the tree sampling method.

    Args:
        n: Number of leaves in the resulting tree that correspond to extant species.
        age: Simulation time, i.e., the time span from the root of the tree to the leaves that
            correspond to extant species.
        model: Simulation model to be applied, the default is 'yule', see [2]. Alternatively the
            birth-death process ('BDP', [3]) is available.
        innovation: If True, use the innovation model, see [4], to sample a lineage for the next
            speciation event. The default is False, in which case the lineage is chosen uniformly at
            random among the currently existing lineages.
        birth_rate: The birth rate for models such as 'yule' and 'BDP'. The default is None, in
            which case the birth rate is set to 1.0 unless episodes are specified.
        death_rate: The death rate for models such as 'BDP'. The default is None, in which case the
            death rate is set to 0.0 unless episodes are specified.
        contraction_probability: Probability that an inner edge is contracted. The default is 0.0,
            in which case the tree is binary. Only one of this parameter and
            'contraction_proportion' may be non-zero.
        contraction_proportion: The proportion of inner edges to be contracted. The default is 0.0,
            in which case the tree is binary. Only one of this parameter and
            'contraction_probability' may be non-zero.
        contraction_bias: Specifies whether shorter edges, i.e., with a smaller difference t of the
            time stamps, have a higher probability to be contracted. Only relevant if
            'contraction_proportion' > 0.0. The default is None, in which case all edges have the
            same probability to be contracted. The options 'inverse' and 'exponential' mean that an
            edge is sampled weighted by 1/(a * t) or e^(-a * t), respectively, where a is a
            user-defined factor.
        bias_strength: Intensity factor for preferring shorter edges to be contracted. The default
            is 1.0.

    Raises:
        ValueError: If unknown or invalid parameter values are passed.

    Returns:
        The simulated species tree.

    References:
        1. G. U. Yule. A mathematical theory of evolution, based on the conclusions of Dr. J.
           C. Willis, F. R. S. In: Phil. Trans. R. Soc. Lond. B, 1924, 213, 21-87.
           doi:10.1098/rstb.1925.0002.
        2. D. G. Kendall. On the Generalized "Birth-and-Death" Process. In: Ann. Math. Statist.
           1948, 19, 1-15. doi:10.1214/aoms/1177730285.
        3. T. Stadler. Simulating trees with a fixed number of extant species. In: Syst. Biol.
           2011, 60, 676-684. doi:10.1093/sysbio/syr029.
        4. S. Keller-Schmidt, K. Klemm. A model of macroevolution as a branching process based
           on innovations. In: Adv. Complex Syst.,2012, 15, 1250043. doi:10.1142/S0219525912500439.
    """
    # parameter checking
    if not isinstance(age, (float, int)) or age <= 0.0:
        raise ValueError("age must be a number >0")
    elif isinstance(age, int):
        age = float(age)

    if not isinstance(n, int) or n < 0:
        raise ValueError("n must be an int >=0")
    elif n == 0:
        return Tree(None)

    if not isinstance(model, str):
        raise ValueError("model must be of type 'str'")

    # main simulation algorithm
    if model.lower() == "yule":
        tree = _yule_n_age(n, age, birth_rate, innovation)
    elif model.upper() == "BDP":
        tree = _BDP_n_age(n, age, birth_rate, death_rate, innovation)
    else:
        raise ValueError(f"model '{model}' is not available")

    # make tree non_binary by random contraction of edges
    nonbinary(
        tree,
        contraction_probability=contraction_probability,
        contraction_proportion=contraction_proportion,
        contraction_bias=contraction_bias,
        bias_strength=bias_strength,
        inplace=True,
    )

    # assign the distance attribute to all vertices
    distance_from_timing(tree)

    return tree

nonbinary(tree: Tree, contraction_probability: float = 0.0, contraction_proportion: float = 0.0, contraction_bias: str | None = None, bias_strength: float = 1.0, inplace: bool = False) -> Tree

Introduce multifurcation into a tree by contraction of inner edges.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree whose edges shall be contracted.	required
contraction_probability	float	Probability that an inner edge is contracted; results in non-binary tree; default is 0.0. Only one of this parameter and 'contraction_proportion' may be non-zero.	0.0
contraction_proportion	float	The proportion of inner edges to be contracted; results in non-binary tree; default is 0.0. Only one of this parameter and 'contraction_probability' may be non-zero.	0.0
contraction_bias	str | None	Specifies whether shorter edges, i.e., with a smaller difference t of the time stamps, have a higher probability to be contracted. Only relevant if 'contraction_proportion' > 0.0. The default is None, in which case all edges have the same probability to be contracted. The options 'inverse' and 'exponential' mean that an edge is sampled weighted by 1/(a * t) or e^(-a * t), respectively, where a is a user-defined factor.	None
bias_strength	float	Intensity factor for preferring shorter edges to be contracted.	1.0
inplace	bool	If True, the edges are contracted in the original tree instance; otherwise a copy of the tree in which edges are contracted is returned.	False

Returns:

Type	Description
Tree	The tree whose edges have been contracted according to the parameters.

Raises:

Type	Description
ValueError	If both contraction_probability and contraction_proportion are non-zero or do not lie in the interval [0, 1].

Source code in asymmetree/treeevolve/species.py

def nonbinary(
    tree: Tree,
    contraction_probability: float = 0.0,
    contraction_proportion: float = 0.0,
    contraction_bias: str | None = None,
    bias_strength: float = 1.0,
    inplace: bool = False,
) -> Tree:
    """Introduce multifurcation into a tree by contraction of inner edges.

    Args:
        tree: The tree whose edges shall be contracted.
        contraction_probability: Probability that an inner edge is contracted; results in non-binary
            tree; default is 0.0. Only one of this parameter and 'contraction_proportion' may be
            non-zero.
        contraction_proportion: The proportion of inner edges to be contracted; results in
            non-binary tree; default is 0.0. Only one of this parameter and
            'contraction_probability' may be non-zero.
        contraction_bias: Specifies whether shorter edges, i.e., with a smaller difference t of the
            time stamps, have a higher probability to be contracted. Only relevant if
            'contraction_proportion' > 0.0. The default is None, in which case all edges have the
            same probability to be contracted. The options 'inverse' and 'exponential' mean that an
            edge is sampled weighted by 1/(a * t) or e^(-a * t), respectively, where a is a
            user-defined factor.
        bias_strength: Intensity factor for preferring shorter edges to be contracted.
        inplace: If True, the edges are contracted in the original tree instance; otherwise a copy
            of the tree in which edges are contracted is returned.

    Returns:
        The tree whose edges have been contracted according to the parameters.

    Raises:
        ValueError: If both contraction_probability and contraction_proportion are non-zero or do
            not lie in the interval [0, 1].
    """
    if contraction_probability != 0.0 and contraction_proportion != 0.0:
        raise ValueError(
            "only one of parameters 'contraction_probability' and 'contraction_proportion' may be "
            "non-zero"
        )

    if contraction_probability < 0.0 or contraction_probability > 1.0:
        raise ValueError("contraction probability must be in [0.0, 1.0]")

    if contraction_proportion < 0.0 or contraction_proportion > 1.0:
        raise ValueError("contraction proportion must be in [0.0, 1.0]")

    if not isinstance(bias_strength, (int, float)) or bias_strength <= 0.0:
        raise ValueError("factor for contraction bias must be > 0")

    if not inplace:
        tree = tree.copy()

    edges = None
    if contraction_probability > 0.0:
        edges = _select_edges_by_probability(
            tree, contraction_probability, exclude_planted_edge=True
        )
    elif contraction_proportion > 0.0:
        edges = _select_edges_by_proportion(
            tree, contraction_proportion, contraction_bias, bias_strength, exclude_planted_edge=True
        )

    if edges:
        tree.contract(edges)
        distance_from_timing(tree)

    return tree

rescale(tree: Tree, height: float, inplace: bool = True) -> Tree

Rescale the time stamps of a tree.

The time stamps of all vertices are multiplied by a scaling factor such that the height of the tree is equal to the specified height. The height of the tree is defined as the time stamp of the root vertex.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree whose time stamps shall be rescaled.	required
height	float	The new height of the tree.	required
inplace	bool	If True, the time stamps are rescaled in the original tree instance; otherwise a copy of the tree with rescaled time stamps is returned.	True

Returns:

Type	Description
Tree	The tree with rescaled time stamps.

Source code in asymmetree/treeevolve/species.py

def rescale(tree: Tree, height: float, inplace: bool = True) -> Tree:
    """Rescale the time stamps of a tree.

    The time stamps of all vertices are multiplied by a scaling factor such that the height of the
    tree is equal to the specified height. The height of the tree is defined as the time stamp of
    the root vertex.

    Args:
        tree: The tree whose time stamps shall be rescaled.
        height: The new height of the tree.
        inplace: If True, the time stamps are rescaled in the original tree instance; otherwise a
            copy of the tree with rescaled time stamps is returned.

    Returns:
        The tree with rescaled time stamps.
    """
    if not inplace:
        tree = tree.copy()

    old_height = tree.root.tstamp

    # not available for trees that only consist of a root
    if old_height <= 0.0:
        raise RuntimeError(f"cannot rescale tree of height {old_height}")

    scaling_factor = height / old_height

    for v in tree.preorder():
        v.tstamp *= scaling_factor

    return tree

assign_losses(tree: Tree, proportion: float) -> None

Randomly assigns a specified proportion of leaves as losses.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree whose leaves shall be assigned as losses.	required
proportion	float	The proportion of leaves to be assigned as losses.	required

Source code in asymmetree/treeevolve/species.py

def assign_losses(tree: Tree, proportion: float) -> None:
    """Randomly assigns a specified proportion of leaves as losses.

    Args:
        tree: The tree whose leaves shall be assigned as losses.
        proportion: The proportion of leaves to be assigned as losses.
    """
    for leaf in tree.random_leaves(proportion):
        leaf.event = "L"

genes

Simulation of dated gene trees.

This module contains the class GeneTreeSimulator for simulating dated gene trees along a given dated species tree. The simulation is based on a birth-death process with speciations as additional branching events, and allows for the simulation of gene duplications, losses, horizontal gene transfers, and gene conversions. A pruning function is provided to remove all branches that lead to losses only, and to suppress superfluous vertices.

GeneTreeSimulator

Class for simulating dated gene trees along a given dated species tree.

Source code in asymmetree/treeevolve/genes.py

class GeneTreeSimulator:
    """Class for simulating dated gene trees along a given dated species tree."""

    def __init__(self, S: Tree) -> None:
        """Initialize the gene tree simulator with a given dated species tree.

        Args:
            S: A dated species tree with a non-empty root.

        Raises:
            TypeError: If S is not a non-empty tree of type 'Tree'.
        """
        if not isinstance(S, Tree) or not S.root:
            raise TypeError("'S' must be a non-empty tree of type 'Tree'")

        self.S = S  # species tree
        self.lca_S = LCA(S)  # LCA data structure for the species tree

        # list of speciations sorted by time stamp
        self.sorted_speciations = sorted_nodes(S)

        # sort edges (u,v) of the species tree by tstamp of u
        self.sorted_edges = sorted_edges(S)

        self._analyze_species_tree()

    def simulate(
        self,
        dupl_rate: float = 0.0,
        loss_rate: float = 0.0,
        hgt_rate: float = 0.0,
        gc_rate: float = 0.0,
        dupl_polytomy: float = 0.0,
        prohibit_extinction: str | bool = "per_species",
        replace_prob: float = 0.0,
        additive_transfer_distance_bias: str | None = None,
        replacing_transfer_distance_bias: str | None = None,
        transfer_distance_bias: str | None = None,
        transfer_distance_bias_strength: float = 1.0,
        gc_distance_bias: str | None = None,
        gc_distance_bias_strength: float = 1.0,
        **kwargs,
    ) -> Tree:
        """Simulate a gene tree along the specified species tree.

        Args:
            dupl_rate: Duplication rate, default is 0.0.
            loss_rate: Loss rate, default is 0.0.
            hgt_rate: Horizontal gene transfer rate, default is 0.0.
            gc_rate: Gene conversion rate, default is 0.0.
            dupl_polytomy: Allows non-binary duplication events by specifying the lambda parameter
                for a poisson distribution (copy number = drawn number + 2); default is 0.0.
            prohibit_extinction: Avoid the extinction of all members in any species ('per_species'),
                of the complete gene family ('per_family'), or no constraints (False); default is
                'per_species'.
            replace_prob: Enables replacing HGT events, probability by which one random homolog in
                the receiving branch of the receiving branch gets lost immediately after the
                transfer. The default is 0.0, in which case all HGT events are additive.
            additive_transfer_distance_bias: Specifies whether closer related species have a higher
                probability to be the recipient species in an additive HGT event. The default is
                None, in which case the recipient species is chosen at random among the co-existing
                species. The options 'inverse' and 'exponential' mean that a species branch is
                sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time
                between the last common ancestor of the two species branches and the time of the
                event, see [1], and a is a user-defined factor.
            replacing_transfer_distance_bias: Specifies whether closer related gene branches have a
                higher probability to be replaced in a replacing HGT event. The default is None, in
                which case the replaced gene is chosen at random among the co-existing gene
                branches. The options 'inverse' and 'exponential' mean that a gene branch is sampled
                weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time between
                the last common ancestor of the two gene branches and the time of the event, see
                [1], and a is a user-defined factor.
            transfer_distance_bias: Set a common bias mode for additive and replacing HGT, see
                description of parameters 'additive_transfer_distance_bias' and
                'replacing_transfer_distance_bias'. This parameter is used if the bias mode shall be
                the same for both types of HGT and the respective parameters are not set, and
                overrides the default bias mode for replacing HGT.
            transfer_distance_bias_strength: Intensity of the transfer distance bias (factor a) for
                additive and replacing HGT, default is 1.0.
            gc_distance_bias: Specifies whether closer related gene branches have a higher
                probability to be replaced in a gene conversion event. The default is None, in which
                case the replaced gene is chosen at random among the paralogs in the respective
                species lineage. The options 'inverse' and 'exponential' mean that a paralog is
                sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time
                between the last common ancestor of the two gene branches and the time of the event,
                see [1], and a is a user-defined factor.
            gc_distance_bias_strength: Intensity of the distance bias (factor a) for gene
                conversion, default is 1.0.

        Returns:
            The simulated gene tree.

        References:
            1. S. Kundu, M. S. Bansal. SaGePhy: an improved phylogenetic simulation framework
               for gene and subgene evolution. In: Bioinformatics, 35(18), 2019, 3496-3498.
               doi:10.1093/bioinformatics/btz081.
        """
        self.rate_sum = 0.0
        for rate in (dupl_rate, loss_rate, hgt_rate, gc_rate):
            self.rate_sum += rate
            if rate < 0.0:
                raise ValueError(f"negative rate {rate}")

        self.d = dupl_rate
        self.l = loss_rate
        self.h = hgt_rate
        self.gc = gc_rate

        if prohibit_extinction not in (False, "per_family", "per_species"):
            raise ValueError(f"unknown mode prohibit_extinction attribute: {prohibit_extinction}")

        self._prohibit_extinction = prohibit_extinction

        self._dupl_polytomy = dupl_polytomy

        self._replace_prob = replace_prob

        for m in (
            additive_transfer_distance_bias,
            replacing_transfer_distance_bias,
            transfer_distance_bias,
            gc_distance_bias,
        ):
            if m not in (None, "inverse", "exponential"):
                raise ValueError(f"unknown mode for transfer distance bias: {m}")

        # mode for transfer distance bias for additive and replacing HGT
        self._additive_transfer_distance_bias = additive_transfer_distance_bias
        self._replacing_transfer_distance_bias = replacing_transfer_distance_bias
        if transfer_distance_bias:
            if not additive_transfer_distance_bias:
                self._additive_transfer_distance_bias = transfer_distance_bias
            if not replacing_transfer_distance_bias:
                self._replacing_transfer_distance_bias = transfer_distance_bias

        # intensity factor for transfer distance bias
        if (
            not isinstance(transfer_distance_bias_strength, (int, float))
            or transfer_distance_bias_strength <= 0.0
        ):
            raise ValueError("factor for transfer distance bias must be > 0")
        self._transfer_distance_bias_strength = transfer_distance_bias_strength

        self._gc_distance_bias = gc_distance_bias
        if (
            not isinstance(gc_distance_bias_strength, (int, float))
            or gc_distance_bias_strength <= 0.0
        ):
            raise ValueError("factor for g.c. distance bias must be > 0")
        self._gc_distance_bias_strength = gc_distance_bias_strength

        return self._run()

    def _analyze_species_tree(self) -> None:
        """Determine the number of surviving species lineages in the subtree of each node."""
        self.S_subtree_survivors = {}

        for u in self.S.postorder():
            if not u.children:
                if u.event == "L":
                    self.S_subtree_survivors[u] = 0
                else:
                    self.S_subtree_survivors[u] = 1
            else:
                self.S_subtree_survivors[u] = sum(self.S_subtree_survivors[v] for v in u.children)

        if not self.S_subtree_survivors[self.S.root]:
            warnings.warn("species tree has no non-loss leaves", category=warnings.UserWarning)

    def _new_branch(self, parent: TreeNode, S_edge: TreeNode, transferred: int) -> _Branch:
        """Create a new gene branch as child of the given parent node.

        Args:
            parent: The parent node in the gene tree.
            S_edge: The species tree edge into which the new branch is embedded (represented by its
                lower node).
            transferred: 1 if the new branch is a transferred branch in an HGT event, 0 otherwise.

        Returns:
            The newly created gene branch.
        """
        new_branch = _Branch(self.id_counter, len(self.branches), parent, S_edge, transferred)
        self.branches.append(new_branch)

        self.ES_to_b[S_edge].append(new_branch)
        self.id_counter += 1

        if self.S_subtree_survivors[S_edge]:
            self.surv_non_loss_lineages.add(new_branch)

        return new_branch

    def _remove_branch(self, branch: _Branch) -> None:
        """Remove a gene branch from the tree.

        Args:
            branch: The gene branch to be removed.
        """
        # efficient removal from list
        if branch.list_id != len(self.branches) - 1:
            branch2 = self.branches[-1]
            self.branches[branch.list_id] = branch2
            branch2.list_id = branch.list_id
        self.branches.pop()

        self.ES_to_b[branch.S_edge].remove(branch)
        self.surv_non_loss_lineages.discard(branch)

    def _initiatialize_tree(self) -> Tree:
        """Initialize the gene tree with a root node and branches for the first speciations."""
        root = TreeNode(
            label=0,
            event="S" if len(self.S.root.children) > 1 else None,
            reconc=self.S.root.label,
            tstamp=self.S.root.tstamp,
        )

        T = Tree(root)
        self.id_counter += 1
        self.spec_queue.popleft()

        for S_edge in self.S.root.children:
            self._new_branch(T.root, S_edge, 0)

        return T

    def _get_branch_and_type(self) -> tuple[_Branch, str]:
        """Randomly sample a branch and an event type according to the rates.

        Returns:
            A tuple containing the sampled branch and the event type ('D', 'L', 'H', or 'GC').
        """
        total_rate = len(self.branches) * self.rate_sum
        r = np.random.uniform(high=total_rate)
        i = int(len(self.branches) * r / total_rate)
        r2 = r - i * self.rate_sum

        if r2 <= self.d:
            event_type = "D"
        elif r2 <= self.d + self.l:
            event_type = "L"
        elif r2 <= self.d + self.l + self.h:
            event_type = "H"
        else:
            event_type = "GC"

        return self.branches[i], event_type

    def _run(self) -> Tree:
        """Run the gene tree simulation.

        The simulation proceeds by iteratively sampling the next event (speciation, duplication,
        loss, HGT, or gene conversion) and updating the gene tree accordingly until all speciation
        events have been processed.

        Returns:
            The simulated gene tree.
        """
        # queue for speciation events
        self.spec_queue = deque(self.sorted_speciations)
        self.id_counter = 0

        # keep track of surving branches that are in species branches with at least 1 surviving
        # species leaf
        self.surv_non_loss_lineages = set()

        self.branches = []

        # maps species tree branches to existing gene branches
        self.ES_to_b = {S_edge: [] for _, S_edge in self.sorted_edges}

        self.T = self._initiatialize_tree()
        t = self.T.root.tstamp

        while self.spec_queue:
            total_rate = len(self.branches) * self.rate_sum
            event_tstamp = t - np.random.exponential(1 / total_rate) if total_rate > 0.0 else -1
            next_spec_tstamp = self.spec_queue[0].tstamp

            # speciation
            if event_tstamp <= next_spec_tstamp:
                self._speciation()
                t = next_spec_tstamp

            else:
                branch, event_type = self._get_branch_and_type()

                # duplication
                if event_type == "D":
                    self._duplication(event_tstamp, branch)

                # loss
                elif event_type == "L":
                    self._loss(event_tstamp, branch)

                # HGT
                elif event_type == "H":
                    self._hgt(event_tstamp, branch)

                # gene conversion
                elif event_type == "GC":
                    self._gene_conversion(event_tstamp, branch)

                t = event_tstamp

        # add the 'dist' attribute to the nodes
        distance_from_timing(self.T)

        return self.T

    def _speciation(self) -> None:
        """Execute the next speciation event in the queue.

        This method also handles loss and leaf nodes of the species tree, as they are represented by
        speciation events with one or zero children, respectively, in the species tree.
        """
        S_edge = self.spec_queue.popleft()

        # copy since we modify this list
        for branch in self.ES_to_b[S_edge].copy():
            spec_node = TreeNode(
                label=branch.label,
                event="S",
                reconc=S_edge.label,
                tstamp=S_edge.tstamp,
                transferred=branch.transferred,
            )
            branch.parent.add_child(spec_node)
            self._remove_branch(branch)

            for S_w in S_edge.children:
                self._new_branch(spec_node, S_w, 0)

            # loss leaves if it was a species extinction event
            if (not S_edge.children) and S_edge.label == "L":
                spec_node.event = "L"

    def _duplication(self, event_tstamp: float, branch: _Branch) -> None:
        """Execute a duplication event on the given branch at the given time stamp.

        This method also allows for non-binary duplication events if the parameter 'dupl_polytomy'
        is set to a positive value.

        Args:
            event_tstamp: The time stamp of the duplication event.
            branch: The gene branch on which the duplication event occurs.
        """
        S_edge = branch.S_edge

        dupl_node = TreeNode(
            label=branch.label,
            event="D",
            reconc=(S_edge.parent.label, S_edge.label),
            tstamp=event_tstamp,
            transferred=branch.transferred,
        )
        branch.parent.add_child(dupl_node)
        self._remove_branch(branch)

        copy_number = (
            2 if self._dupl_polytomy <= 0.0 else 2 + np.random.poisson(lam=self._dupl_polytomy)
        )

        for i in range(copy_number):
            self._new_branch(dupl_node, S_edge, 0)

    def _loss(self, event_tstamp: float, branch: _Branch) -> None:
        """Execute a loss event on the given branch at the given time stamp.

        This method respects the extinction prohibition settings.

        Args:
            event_tstamp: The time stamp of the loss event.
            branch: The gene branch on which the loss event occurs.
        """
        # not executing the loss event if extinction is prohibited is equivalent to temporarily
        # setting the loss rate to zero in the respective branches
        if (
            self._prohibit_extinction == "per_family"
            and len(self.surv_non_loss_lineages) == 1
            and next(iter(self.surv_non_loss_lineages)) is branch
        ):
            return

        if self._prohibit_extinction == "per_species" and len(self.ES_to_b[branch.S_edge]) <= 1:
            return

        S_edge = branch.S_edge

        loss_node = TreeNode(
            label=branch.label,
            event="L",
            reconc=(S_edge.parent.label, S_edge.label),
            tstamp=event_tstamp,
            transferred=branch.transferred,
        )
        branch.parent.add_child(loss_node)
        self._remove_branch(branch)

    def _coexisting_species_edges(self, tstamp: float, exclude_edge: _Branch = None) -> list:
        """Return list of edges for the given timestamp.

        Args:
            tstamp: The time stamp for which to find coexisting species edges.
            exclude_edge: An optional species edge to exclude from the returned list (e.g., the edge
                into which a branch is currently embedded, to avoid self-transfers in HGT events).

        Returns:
            A list of species tree edges that are coexisting at the given time stamp and do not
                match the exclude_edge.
        """
        valid_species = []

        for S_u, S_v in self.sorted_edges:
            if S_u.tstamp <= tstamp:
                break
            elif S_v.tstamp < tstamp and S_v != exclude_edge:
                valid_species.append(S_v)

        return valid_species

    def _sample_recipient(
        self,
        event_tstamp: float,
        branch: _Branch,
    ) -> tuple[TreeNode | None, _Branch | None]:
        """Sample a recipient species edge for an HGT event on the given branch at the given time.

        This method samples a recipient species edge according to the settings for transfer distance
        bias and the replace_prob parameter, and also returns the replaced gene branch in case of a
        replacing HGT event.

        Args:
            event_tstamp: The time stamp of the HGT event.
            branch: The gene branch on which the HGT event occurs.

        Returns:
            A tuple containing the recipient species edge and the replaced gene branch (if any).
        """
        S_edge = branch.S_edge
        trans_edge, replaced_gene_branch = None, None

        valid_species = self._coexisting_species_edges(event_tstamp, exclude_edge=S_edge)

        if not valid_species:
            return None, None

        # ---- ADDITIVE HGT EVENT ----
        if np.random.random() >= self._replace_prob:
            # ---- no transfer distance bias ---
            if not self._additive_transfer_distance_bias:
                trans_edge = random.choice(valid_species)
            # ---- transfer distance bias ---
            else:
                distances = [(self.lca_S(S_edge, e).tstamp - event_tstamp) for e in valid_species]
                a = self._transfer_distance_bias_strength
                if self._additive_transfer_distance_bias == "inverse":
                    weights = 1 / (a * np.asarray(distances))
                elif self._additive_transfer_distance_bias == "exponential":
                    weights = np.exp(-a * np.asarray(distances))

                trans_edge = random.choices(valid_species, weights=weights)[0]

        # ---- REPLACING HGT EVENT ----
        else:
            valid_genes = [b for e in valid_species for b in self.ES_to_b[e]]

            if not valid_genes:
                return None, None

            # ---- no transfer distance bias ---
            if not self._replacing_transfer_distance_bias:
                replaced_gene_branch = random.choice(valid_genes)
            # ---- transfer distance bias ---
            else:
                lca_T = LCA(self.T)
                distances = [
                    (lca_T(branch.parent, b.parent).tstamp - event_tstamp) for b in valid_genes
                ]
                a = self._transfer_distance_bias_strength
                if self._replacing_transfer_distance_bias == "inverse":
                    weights = 1 / (a * np.asarray(distances))
                elif self._replacing_transfer_distance_bias == "exponential":
                    weights = np.exp(-a * np.asarray(distances))

                replaced_gene_branch = random.choices(valid_genes, weights=weights)[0]

            trans_edge = replaced_gene_branch.S_edge

        return trans_edge, replaced_gene_branch

    def _hgt(self, event_tstamp: float, branch: _Branch) -> None:
        """Execute a horizontal gene transfer event on the given branch at the given time stamp.

        Args:
            event_tstamp: The time stamp of the HGT event.
            branch: The gene branch on which the HGT event occurs.
        """
        S_edge = branch.S_edge

        trans_edge, replaced_gene_branch = self._sample_recipient(event_tstamp, branch)

        if trans_edge:
            hgt_node = TreeNode(
                label=branch.label,
                event="H",
                reconc=(S_edge.parent.label, S_edge.label),
                tstamp=event_tstamp,
                transferred=branch.transferred,
            )
            branch.parent.add_child(hgt_node)
            self._remove_branch(branch)

            # original branch
            self._new_branch(hgt_node, S_edge, 0)

            # receiving branch
            self._new_branch(hgt_node, trans_edge, 1)

        # replacing HGT leads to loss in the recipient species
        if replaced_gene_branch:
            self._loss(event_tstamp, replaced_gene_branch)
            # save replaced gene information in HGT node
            hgt_node.replaced_gene = replaced_gene_branch.label

    def _sample_replaced_gene_gc(self, event_tstamp: float, branch: _Branch) -> _Branch | None:
        """Sample a gene branch to be replaced in a gene conversion event.

        Args:
            event_tstamp: The time stamp of the gene conversion event.
            branch: The gene branch on which the gene conversion event occurs.

        Returns:
            The gene branch to be replaced, or None if no replacement is possible.
        """
        candidates = [b for b in self.ES_to_b[branch.S_edge] if b is not branch]

        # no gene conversion possible if there is currently no other branch in the species lineage
        if not candidates:
            return None

        # no distance bias, sample randomly among the candidates
        if not self._gc_distance_bias:
            return random.choice(candidates)

        # distance bias, sample weighted by distance to the common ancestor in the gene tree

        lca_T = LCA(self.T)
        distances = [(lca_T(branch.parent, b.parent).tstamp - event_tstamp) for b in candidates]
        a = self._gc_distance_bias_strength

        if self._gc_distance_bias == "inverse":
            weights = 1 / (a * np.asarray(distances))
        elif self._gc_distance_bias == "exponential":
            weights = np.exp(-a * np.asarray(distances))
        else:
            raise ValueError(f"unknown mode for g.c. distance bias: {self._gc_distance_bias}")

        return random.choices(candidates, weights=weights)[0]

    def _gene_conversion(self, event_tstamp: float, branch: _Branch) -> None:
        """Execute a gene conversion event on the given branch at the given time stamp.

        Args:
            event_tstamp: The time stamp of the gene conversion event.
            branch: The gene branch on which the gene conversion event occurs.
        """
        S_edge = branch.S_edge

        replaced_gene_branch = self._sample_replaced_gene_gc(event_tstamp, branch)

        if replaced_gene_branch:
            gc_node = TreeNode(
                label=branch.label,
                event="GC",
                reconc=(S_edge.parent.label, S_edge.label),
                tstamp=event_tstamp,
                transferred=branch.transferred,
            )
            branch.parent.add_child(gc_node)
            self._remove_branch(branch)

            self._new_branch(gc_node, S_edge, 0)
            self._new_branch(gc_node, S_edge, 0)

            # gene conversion leads to loss of a homologous gene
            self._loss(event_tstamp, replaced_gene_branch)
            # save replaced gene information in GC node
            gc_node.replaced_gene = replaced_gene_branch.label

__init__(S: Tree) -> None

Initialize the gene tree simulator with a given dated species tree.

Parameters:

Name	Type	Description	Default
S	Tree	A dated species tree with a non-empty root.	required

Raises:

Type	Description
TypeError	If S is not a non-empty tree of type 'Tree'.

Source code in asymmetree/treeevolve/genes.py

def __init__(self, S: Tree) -> None:
    """Initialize the gene tree simulator with a given dated species tree.

    Args:
        S: A dated species tree with a non-empty root.

    Raises:
        TypeError: If S is not a non-empty tree of type 'Tree'.
    """
    if not isinstance(S, Tree) or not S.root:
        raise TypeError("'S' must be a non-empty tree of type 'Tree'")

    self.S = S  # species tree
    self.lca_S = LCA(S)  # LCA data structure for the species tree

    # list of speciations sorted by time stamp
    self.sorted_speciations = sorted_nodes(S)

    # sort edges (u,v) of the species tree by tstamp of u
    self.sorted_edges = sorted_edges(S)

    self._analyze_species_tree()

simulate(dupl_rate: float = 0.0, loss_rate: float = 0.0, hgt_rate: float = 0.0, gc_rate: float = 0.0, dupl_polytomy: float = 0.0, prohibit_extinction: str | bool = 'per_species', replace_prob: float = 0.0, additive_transfer_distance_bias: str | None = None, replacing_transfer_distance_bias: str | None = None, transfer_distance_bias: str | None = None, transfer_distance_bias_strength: float = 1.0, gc_distance_bias: str | None = None, gc_distance_bias_strength: float = 1.0, **kwargs) -> Tree

Simulate a gene tree along the specified species tree.

Parameters:

Name	Type	Description	Default
dupl_rate	float	Duplication rate, default is 0.0.	0.0
loss_rate	float	Loss rate, default is 0.0.	0.0
hgt_rate	float	Horizontal gene transfer rate, default is 0.0.	0.0
gc_rate	float	Gene conversion rate, default is 0.0.	0.0
dupl_polytomy	float	Allows non-binary duplication events by specifying the lambda parameter for a poisson distribution (copy number = drawn number + 2); default is 0.0.	0.0
prohibit_extinction	str | bool	Avoid the extinction of all members in any species ('per_species'), of the complete gene family ('per_family'), or no constraints (False); default is 'per_species'.	'per_species'
replace_prob	float	Enables replacing HGT events, probability by which one random homolog in the receiving branch of the receiving branch gets lost immediately after the transfer. The default is 0.0, in which case all HGT events are additive.	0.0
additive_transfer_distance_bias	str | None	Specifies whether closer related species have a higher probability to be the recipient species in an additive HGT event. The default is None, in which case the recipient species is chosen at random among the co-existing species. The options 'inverse' and 'exponential' mean that a species branch is sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time between the last common ancestor of the two species branches and the time of the event, see [1], and a is a user-defined factor.	None
replacing_transfer_distance_bias	str | None	Specifies whether closer related gene branches have a higher probability to be replaced in a replacing HGT event. The default is None, in which case the replaced gene is chosen at random among the co-existing gene branches. The options 'inverse' and 'exponential' mean that a gene branch is sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time between the last common ancestor of the two gene branches and the time of the event, see [1], and a is a user-defined factor.	None
transfer_distance_bias	str | None	Set a common bias mode for additive and replacing HGT, see description of parameters 'additive_transfer_distance_bias' and 'replacing_transfer_distance_bias'. This parameter is used if the bias mode shall be the same for both types of HGT and the respective parameters are not set, and overrides the default bias mode for replacing HGT.	None
transfer_distance_bias_strength	float	Intensity of the transfer distance bias (factor a) for additive and replacing HGT, default is 1.0.	1.0
gc_distance_bias	str | None	Specifies whether closer related gene branches have a higher probability to be replaced in a gene conversion event. The default is None, in which case the replaced gene is chosen at random among the paralogs in the respective species lineage. The options 'inverse' and 'exponential' mean that a paralog is sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time between the last common ancestor of the two gene branches and the time of the event, see [1], and a is a user-defined factor.	None
gc_distance_bias_strength	float	Intensity of the distance bias (factor a) for gene conversion, default is 1.0.	1.0

Returns:

Type	Description
Tree	The simulated gene tree.

References

1. S. Kundu, M. S. Bansal. SaGePhy: an improved phylogenetic simulation framework for gene and subgene evolution. In: Bioinformatics, 35(18), 2019, 3496-3498. doi:10.1093/bioinformatics/btz081.

Source code in asymmetree/treeevolve/genes.py

def simulate(
    self,
    dupl_rate: float = 0.0,
    loss_rate: float = 0.0,
    hgt_rate: float = 0.0,
    gc_rate: float = 0.0,
    dupl_polytomy: float = 0.0,
    prohibit_extinction: str | bool = "per_species",
    replace_prob: float = 0.0,
    additive_transfer_distance_bias: str | None = None,
    replacing_transfer_distance_bias: str | None = None,
    transfer_distance_bias: str | None = None,
    transfer_distance_bias_strength: float = 1.0,
    gc_distance_bias: str | None = None,
    gc_distance_bias_strength: float = 1.0,
    **kwargs,
) -> Tree:
    """Simulate a gene tree along the specified species tree.

    Args:
        dupl_rate: Duplication rate, default is 0.0.
        loss_rate: Loss rate, default is 0.0.
        hgt_rate: Horizontal gene transfer rate, default is 0.0.
        gc_rate: Gene conversion rate, default is 0.0.
        dupl_polytomy: Allows non-binary duplication events by specifying the lambda parameter
            for a poisson distribution (copy number = drawn number + 2); default is 0.0.
        prohibit_extinction: Avoid the extinction of all members in any species ('per_species'),
            of the complete gene family ('per_family'), or no constraints (False); default is
            'per_species'.
        replace_prob: Enables replacing HGT events, probability by which one random homolog in
            the receiving branch of the receiving branch gets lost immediately after the
            transfer. The default is 0.0, in which case all HGT events are additive.
        additive_transfer_distance_bias: Specifies whether closer related species have a higher
            probability to be the recipient species in an additive HGT event. The default is
            None, in which case the recipient species is chosen at random among the co-existing
            species. The options 'inverse' and 'exponential' mean that a species branch is
            sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time
            between the last common ancestor of the two species branches and the time of the
            event, see [1], and a is a user-defined factor.
        replacing_transfer_distance_bias: Specifies whether closer related gene branches have a
            higher probability to be replaced in a replacing HGT event. The default is None, in
            which case the replaced gene is chosen at random among the co-existing gene
            branches. The options 'inverse' and 'exponential' mean that a gene branch is sampled
            weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time between
            the last common ancestor of the two gene branches and the time of the event, see
            [1], and a is a user-defined factor.
        transfer_distance_bias: Set a common bias mode for additive and replacing HGT, see
            description of parameters 'additive_transfer_distance_bias' and
            'replacing_transfer_distance_bias'. This parameter is used if the bias mode shall be
            the same for both types of HGT and the respective parameters are not set, and
            overrides the default bias mode for replacing HGT.
        transfer_distance_bias_strength: Intensity of the transfer distance bias (factor a) for
            additive and replacing HGT, default is 1.0.
        gc_distance_bias: Specifies whether closer related gene branches have a higher
            probability to be replaced in a gene conversion event. The default is None, in which
            case the replaced gene is chosen at random among the paralogs in the respective
            species lineage. The options 'inverse' and 'exponential' mean that a paralog is
            sampled weighted by 1/(a * t) or e^(-(a * t)), resp., where t is the elapsed time
            between the last common ancestor of the two gene branches and the time of the event,
            see [1], and a is a user-defined factor.
        gc_distance_bias_strength: Intensity of the distance bias (factor a) for gene
            conversion, default is 1.0.

    Returns:
        The simulated gene tree.

    References:
        1. S. Kundu, M. S. Bansal. SaGePhy: an improved phylogenetic simulation framework
           for gene and subgene evolution. In: Bioinformatics, 35(18), 2019, 3496-3498.
           doi:10.1093/bioinformatics/btz081.
    """
    self.rate_sum = 0.0
    for rate in (dupl_rate, loss_rate, hgt_rate, gc_rate):
        self.rate_sum += rate
        if rate < 0.0:
            raise ValueError(f"negative rate {rate}")

    self.d = dupl_rate
    self.l = loss_rate
    self.h = hgt_rate
    self.gc = gc_rate

    if prohibit_extinction not in (False, "per_family", "per_species"):
        raise ValueError(f"unknown mode prohibit_extinction attribute: {prohibit_extinction}")

    self._prohibit_extinction = prohibit_extinction

    self._dupl_polytomy = dupl_polytomy

    self._replace_prob = replace_prob

    for m in (
        additive_transfer_distance_bias,
        replacing_transfer_distance_bias,
        transfer_distance_bias,
        gc_distance_bias,
    ):
        if m not in (None, "inverse", "exponential"):
            raise ValueError(f"unknown mode for transfer distance bias: {m}")

    # mode for transfer distance bias for additive and replacing HGT
    self._additive_transfer_distance_bias = additive_transfer_distance_bias
    self._replacing_transfer_distance_bias = replacing_transfer_distance_bias
    if transfer_distance_bias:
        if not additive_transfer_distance_bias:
            self._additive_transfer_distance_bias = transfer_distance_bias
        if not replacing_transfer_distance_bias:
            self._replacing_transfer_distance_bias = transfer_distance_bias

    # intensity factor for transfer distance bias
    if (
        not isinstance(transfer_distance_bias_strength, (int, float))
        or transfer_distance_bias_strength <= 0.0
    ):
        raise ValueError("factor for transfer distance bias must be > 0")
    self._transfer_distance_bias_strength = transfer_distance_bias_strength

    self._gc_distance_bias = gc_distance_bias
    if (
        not isinstance(gc_distance_bias_strength, (int, float))
        or gc_distance_bias_strength <= 0.0
    ):
        raise ValueError("factor for g.c. distance bias must be > 0")
    self._gc_distance_bias_strength = gc_distance_bias_strength

    return self._run()

dated_gene_tree(S: Tree, **kwargs: object) -> Tree

Simulate a gene tree along the specified species tree.

Parameters:

Name	Type	Description	Default
S	Tree	A dated species tree of type 'Tree' with a non-empty root.	required
**kwargs	object	Keyword arguments passed to :meth:GeneTreeSimulator.simulate. Accepted parameters are dupl_rate, loss_rate, hgt_rate, gc_rate, dupl_polytomy, prohibit_extinction, replace_prob, additive_transfer_distance_bias, replacing_transfer_distance_bias, transfer_distance_bias, transfer_distance_bias_strength, gc_distance_bias, and gc_distance_bias_strength. See :meth:GeneTreeSimulator.simulate for the full parameter descriptions.	{}

Returns:

Type	Description
Tree	The simulated gene tree.

Source code in asymmetree/treeevolve/genes.py

def dated_gene_tree(S: Tree, **kwargs: object) -> Tree:
    """Simulate a gene tree along the specified species tree.

    Args:
        S: A dated species tree of type 'Tree' with a non-empty root.
        **kwargs: Keyword arguments passed to :meth:`GeneTreeSimulator.simulate`. Accepted
            parameters are ``dupl_rate``, ``loss_rate``, ``hgt_rate``, ``gc_rate``,
            ``dupl_polytomy``, ``prohibit_extinction``, ``replace_prob``,
            ``additive_transfer_distance_bias``, ``replacing_transfer_distance_bias``,
            ``transfer_distance_bias``, ``transfer_distance_bias_strength``,
            ``gc_distance_bias``, and ``gc_distance_bias_strength``. See
            :meth:`GeneTreeSimulator.simulate` for the full parameter descriptions.

    Returns:
        The simulated gene tree.
    """
    gene_tree_simulator = GeneTreeSimulator(S)

    return gene_tree_simulator.simulate(**kwargs)

prune_losses(tree: Tree) -> Tree

Prune all loss branches.

Returns a copy of the tree with all branches that lead to losses only removed and superfluous vertices suppressed. Additionally, if the root of the tree has only a single child, then this 'planted edge' is also removed.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree to be pruned.	required

Returns:

Type	Description
Tree	A pruned version of the tree.

Source code in asymmetree/treeevolve/genes.py

def prune_losses(tree: Tree) -> Tree:
    """Prune all loss branches.

    Returns a copy of the tree with all branches that lead to losses only removed and superfluous
    vertices suppressed. Additionally, if the root of the tree has only a single child, then this
    'planted edge' is also removed.

    Args:
        tree: The tree to be pruned.

    Returns:
        A pruned version of the tree.
    """
    pruned_tree = delete_losses_and_contract(tree, inplace=False)

    remove_planted_root(pruned_tree, inplace=True)

    return pruned_tree

observable_tree(tree: Tree) -> Tree

Prune all loss branches (deprecated).

This method is deprecated. Use prune_losses() instead.

Parameters:

Name	Type	Description	Default
tree	Tree	The tree to be pruned.	required

Returns:

Type	Description
Tree	A pruned version of the tree.

Source code in asymmetree/treeevolve/genes.py

def observable_tree(tree: Tree) -> Tree:
    """Prune all loss branches (deprecated).

    This method is deprecated. Use prune_losses() instead.

    Args:
        tree: The tree to be pruned.

    Returns:
        A pruned version of the tree.
    """
    warn("This method is deprecated. Use prune_losses() instead.", DeprecationWarning, stacklevel=2)

    return prune_losses(tree)

rate_heterogeneity

Evolution rate heterogeneity.

This module contains functions to introduce evolution rate asymmetries and autocorrelation. The main function is 'rate_heterogeneity()', which modifies the edge lengths of a gene tree according to three sources of heterogeneity: (i) gene-family-specific heterogeneity, (ii) species-specific heterogeneity, and (iii) paralog-specific heterogeneity. The function 'gene_trees()' can be used to simulate gene trees with rate heterogeneity along a species tree.

gene_trees(S: Tree, n: int = 1, dupl_rate: float | tuple = 0.0, loss_rate: float | tuple = 0.0, hgt_rate: float | tuple = 0.0, base_rate: float | tuple = 1.0, **kwargs: object) -> list[Tree]

Simulates dated gene trees with non-ultrametric edge lengths along a species tree.

Parameters:

Name	Type	Description	Default
S	Tree	The species tree along which the gene trees are simulated.	required
n	int	Number of gene trees to be simulated, default is 1, in which case a tree is returned, otherwise a list is returned.	1
dupl_rate	float | tuple	The (distribution for the) duplication rate, see documentation for available option. The default is constant 0.0.	0.0
loss_rate	float | tuple	The (distribution for the) loss rate, see documentation for available option. The default is constant 0.0.	0.0
hgt_rate	float | tuple	The (distribution for the) HGT rate, see documentation for available option. The default is constant 0.0.	0.0
base_rate	float | tuple	The (distribution for the) evolution rate at the roots of the gene trees, see documentation for available options. The default is constant 1.0.	1.0
**kwargs	object	Keyword arguments passed to :meth:GeneTreeSimulator.simulate and :func:rate_heterogeneity. Accepted parameters are dupl_rate, loss_rate, hgt_rate, gc_rate, dupl_polytomy, prohibit_extinction, replace_prob, additive_transfer_distance_bias, replacing_transfer_distance_bias, transfer_distance_bias, transfer_distance_bias_strength, gc_distance_bias, gc_distance_bias_strength, autocorr_variance, rate_increase, and CSN_weights. See :meth:GeneTreeSimulator.simulate and :func:rate_heterogeneity documentations for the full parameter descriptions.	{}

Returns:

Type	Description
list[Tree]	A list of simulated gene tree.

Source code in asymmetree/treeevolve/rate_heterogeneity.py

def gene_trees(
    S: Tree,
    n: int = 1,
    dupl_rate: float | tuple = 0.0,
    loss_rate: float | tuple = 0.0,
    hgt_rate: float | tuple = 0.0,
    base_rate: float | tuple = 1.0,
    **kwargs: object,
) -> list[Tree]:
    """Simulates dated gene trees with non-ultrametric edge lengths along a species tree.

    Args:
        S: The species tree along which the gene trees are simulated.
        n: Number of gene trees to be simulated, default is 1, in which case a tree is returned,
            otherwise a list is returned.
        dupl_rate: The (distribution for the) duplication rate, see documentation for available
            option. The default is constant 0.0.
        loss_rate: The (distribution for the) loss rate, see documentation for available option.
            The default is constant 0.0.
        hgt_rate: The (distribution for the) HGT rate, see documentation for available option. The
            default is constant 0.0.
        base_rate: The (distribution for the) evolution rate at the roots of the gene trees, see
            documentation for available options. The default is constant 1.0.
        **kwargs: Keyword arguments passed to :meth:`GeneTreeSimulator.simulate` and
            :func:`rate_heterogeneity`. Accepted parameters are ``dupl_rate``, ``loss_rate``,
            ``hgt_rate``, ``gc_rate``, ``dupl_polytomy``, ``prohibit_extinction``, ``replace_prob``,
            ``additive_transfer_distance_bias``, ``replacing_transfer_distance_bias``,
            ``transfer_distance_bias``, ``transfer_distance_bias_strength``,
            ``gc_distance_bias``, ``gc_distance_bias_strength``, ``autocorr_variance``,
            ``rate_increase``, and ``CSN_weights``. See :meth:`GeneTreeSimulator.simulate` and
             :func:`rate_heterogeneity` documentations for the full parameter descriptions.

    Returns:
        A list of simulated gene tree.
    """
    gene_trees = []
    simulator = GeneTreeSimulator(S)

    dupl_rate_sampler = Sampler(dupl_rate)
    loss_rate_sampler = Sampler(loss_rate)
    hgt_rate_sampler = Sampler(hgt_rate)
    base_rate_sampler = Sampler(base_rate)

    # autocorrelation between genes of the same or related species
    autocorr_variance = kwargs.pop("autocorr_variance", 0.0)
    _, autocorr_factors = autocorrelation_factors(S, autocorr_variance)

    # main simulation and imbalancing
    for _ in range(n):
        TGT = simulator.simulate(
            dupl_rate=dupl_rate_sampler(),
            loss_rate=loss_rate_sampler(),
            hgt_rate=hgt_rate_sampler(),
            **kwargs,
        )
        rate_heterogeneity(
            TGT, S, base_rate=base_rate_sampler(), autocorr_factors=autocorr_factors, **kwargs
        )
        gene_trees.append(TGT)

    return gene_trees

rate_heterogeneity(T: Tree, S: Tree, base_rate: float = 1.0, autocorr_factors: dict | None = None, autocorr_variance: float = 0.0, rate_increase: float | tuple = ('gamma', 0.5, 2.2), CSN_weights: tuple = (1, 1, 1), inplace: bool = True, **kwargs) -> Tree

Introduces evolution rate heterogeneity into a gene tree.

The function applies rate multiplier for three sources of evolution rate heterogeneity: (i) gene-family-specific heterogeneity is modeled through the base_rate parameter; (ii) for species-specific heterogeneity an autocorrelated relaxed molecular clock model as in [1] is used; and (iii) for paralog-specific heterogeneity it is chosen between the three modes conservation, subfunctionalization, and neofunctionalization for the rates of the descendant lineages of a duplication event. The assigned rates are used to modify the length ('dist' attribute) of the edges of the (originally ultrametric) dated gene tree.

Parameters:

Name	Type	Description	Default
T	Tree	The gene tree.	required
S	Tree	The species tree.	required
base_rate	float	Mean of substitution rate for conserved genes.	1.0
autocorr_factors	dict | None	A dictionary containing autocorrelation rate factors for the edges of S (key = v.label for edge (v.parent, v); value = the rate as a float). The default is None, in which case autocorrelation factors are generated if 'autocorr_variance' > 0.0, or no such modification is applied. See [1] for theoretical background.	None
autocorr_variance	float	Autocorrelation variance factor for a lognormal distribution, only considered if autocorrelation_rates are not supplied. See [1] for theoretical background.	0.0
rate_increase	float | tuple	Distribution of the (relative) rate increase (w.r.t. the base rate) for divergent genes, i.e. to a factor 1 + x. The default is a Gamma distribution with shape 0.5 and scale 2.2, which was fitted to the data in [2].	('gamma', 0.5, 2.2)
CSN_weights	tuple	Weights for choice between conservation, subfunctionalization and neofunctionalization. The default is (1, 1, 1), i.e., all three modes are equally likely to be chosen at each duplication event.	(1, 1, 1)
inplace	bool	If False, copy the tree before imbalancing. The default is True.	True

Returns:

Type	Description
Tree	The original instance of the gene tree (inplace=True) or a copy of the gene tree (inplace=False) with modified dist attributes of the nodes.

References

1. H. Kishino, J. L. Thorne, and W. J. Bruno. Performance of a Divergence Time Estimation Method under a Probabilistic Model of Rate Evolution. In: Molecular Biology and Evolution, 18(3):352-361, March 2001. doi: 10.1093/oxfordjournals.molbev.a003811.
2. K. P. Byrne and K. H. Wolfe. Consistent Patterns of Rate Asymmetry and Gene Loss Indicate Widespread Neofunctionalization of Yeast Genes After Whole-Genome Duplication. In: Genetics, 175(3):1341-1350, March 2007. doi: 10.1534/genetics.106.066951.

Source code in asymmetree/treeevolve/rate_heterogeneity.py

def rate_heterogeneity(
    T: Tree,
    S: Tree,
    base_rate: float = 1.0,
    autocorr_factors: dict | None = None,
    autocorr_variance: float = 0.0,
    rate_increase: float | tuple = ("gamma", 0.5, 2.2),
    CSN_weights: tuple = (1, 1, 1),
    inplace: bool = True,
    **kwargs,
) -> Tree:
    """Introduces evolution rate heterogeneity into a gene tree.

    The function applies rate multiplier for three sources of evolution rate heterogeneity: (i)
    gene-family-specific heterogeneity is modeled through the `base_rate` parameter; (ii) for
    species-specific heterogeneity an autocorrelated relaxed molecular clock model as in [1] is
    used; and (iii) for paralog-specific heterogeneity it is chosen between the three modes
    conservation, subfunctionalization, and neofunctionalization for the rates of the descendant
    lineages of a duplication event. The assigned rates are used to modify the length ('dist'
    attribute) of the edges of the (originally ultrametric) dated gene tree.

    Args:
        T: The gene tree.
        S: The species tree.
        base_rate: Mean of substitution rate for conserved genes.
        autocorr_factors: A dictionary containing autocorrelation rate factors for the edges of S
            (key = v.label for edge (v.parent, v); value = the rate as a float). The default is
            None, in which case autocorrelation factors are generated if 'autocorr_variance' > 0.0,
            or no such modification is applied. See [1] for theoretical background.
        autocorr_variance: Autocorrelation variance factor for a lognormal distribution, only
            considered if `autocorrelation_rates` are not supplied. See [1] for theoretical
            background.
        rate_increase: Distribution of the (relative) rate increase (w.r.t. the base rate) for
            divergent genes, i.e. to a factor 1 + x. The default is a Gamma distribution with shape
            0.5 and scale 2.2, which was fitted to the data in [2].
        CSN_weights: Weights for choice between conservation, subfunctionalization and
            neofunctionalization. The default is (1, 1, 1), i.e., all three modes are equally likely
            to be chosen at each duplication event.
        inplace: If False, copy the tree before imbalancing. The default is True.

    Returns:
        The original instance of the gene tree (`inplace=True`) or a copy of the gene tree
            (`inplace=False`) with modified `dist` attributes of the nodes.

    References:
        1. H. Kishino, J. L. Thorne, and W. J. Bruno. Performance of a Divergence Time
           Estimation Method under a Probabilistic Model of Rate Evolution. In: Molecular Biology
           and Evolution, 18(3):352-361, March 2001. doi: 10.1093/oxfordjournals.molbev.a003811.
        2. K. P. Byrne and K. H. Wolfe. Consistent Patterns of Rate Asymmetry and Gene Loss
           Indicate Widespread Neofunctionalization of Yeast Genes After Whole-Genome Duplication.
           In: Genetics, 175(3):1341-1350, March 2007. doi: 10.1534/genetics.106.066951.
    """
    if not inplace:
        T = T.copy()

    # paralog-specific rate heterogeneity (neo- and subfunctionalization)
    _divergent_rates(
        T, S, Sampler(rate_increase, shift=1.0), np.asarray(CSN_weights) / sum(CSN_weights)
    )

    # species-specific rate heterogeneity (autocorrelation)
    if not autocorr_factors and autocorr_variance > 0.0:
        _, autocorr_factors = autocorrelation_factors(S, autocorr_variance)
    if autocorr_factors:
        for v in T.preorder():
            if v.parent:
                edge_ID = v.reconc[1] if isinstance(v.reconc, (tuple, list)) else v.reconc
                v.dist *= autocorr_factors[edge_ID]

    # gene-family-specific rate heterogeneity (base rate)
    for v in T.preorder():
        v.dist *= base_rate

    return T

assign_rates(T: Tree, S: Tree, **kwargs: object) -> Tree

Deprecated alias for rate_heterogeneity().

Parameters:

Name	Type	Description	Default
T	Tree	The gene tree.	required
S	Tree	The species tree.	required
**kwargs	object	Keyword arguments passed to :func:rate_heterogeneity. Accepted parameters are base_rate, autocorr_factors, autocorr_variance, rate_increase, and CSN_weights. See :func:rate_heterogeneity documentation for the full parameter descriptions.	{}

Returns:

Type	Description
Tree	The original instance of the gene tree with modified dist attributes of the nodes.

Source code in asymmetree/treeevolve/rate_heterogeneity.py

def assign_rates(T: Tree, S: Tree, **kwargs: object) -> Tree:
    """Deprecated alias for rate_heterogeneity().

    Args:
        T: The gene tree.
        S: The species tree.
        **kwargs: Keyword arguments passed to :func:`rate_heterogeneity`. Accepted parameters are
            ``base_rate``, ``autocorr_factors``, ``autocorr_variance``, ``rate_increase``, and
            ``CSN_weights``. See :func:`rate_heterogeneity` documentation for the full parameter
            descriptions.

    Returns:
        The original instance of the gene tree with modified `dist` attributes of the nodes.
    """
    warn(
        "This method is deprecated. Use rate_heterogeneity() instead.",
        DeprecationWarning,
        stacklevel=2,
    )
    return rate_heterogeneity(T, S, **kwargs)

autocorrelation_factors(tree: Tree, variance: float) -> tuple[dict, dict]

Lognormal model to assign rate factors to a species tree.

The parameter variance is a hyperparameter for a log-normal distribution from which offspring rates are drawn. The overall variance of this distribution is variance * divergence time. The rates are first computed for the nodes, the rates of the edges are assigned afterwards as the arithmetic mean of the rates of the two incident nodes.

Parameters:

Name	Type	Description	Default
tree	Tree	The species tree.	required
variance	float	The hyperparameter for a log-normal distribution from which offspring rates are drawn at each node.	required

Returns:

Type	Description
dict	A dict mapping the labels of the nodes to their assigned rates.
dict	A dict mapping the labels of the edges (v.parent, v) to the assigned rates of the edges.

References

1. H. Kishino, J. L. Thorne, and W. J. Bruno. Performance of a Divergence Time Estimation Method under a Probabilistic Model of Rate Evolution. In: Molecular Biology and Evolution, 18(3):352-361, March 2001. doi: 10.1093/oxfordjournals.molbev.a003811.

Source code in asymmetree/treeevolve/rate_heterogeneity.py

def autocorrelation_factors(
    tree: Tree,
    variance: float,
) -> tuple[dict, dict]:
    """Lognormal model to assign rate factors to a species tree.

    The parameter `variance` is a hyperparameter for a log-normal distribution from which offspring
    rates are drawn. The overall variance of this distribution is `variance` * divergence time.
    The rates are first computed for the nodes, the rates of the edges are assigned afterwards as
    the arithmetic mean of the rates of the two incident nodes.

    Args:
        tree: The species tree.
        variance: The hyperparameter for a log-normal distribution from which offspring rates are
            drawn at each node.

    Returns:
        A `dict` mapping the labels of the nodes to their assigned rates.
        A `dict` mapping the labels of the edges `(v.parent, v)` to the assigned rates of the edges.

    References:
        1. H. Kishino, J. L. Thorne, and W. J. Bruno. Performance of a Divergence Time
           Estimation Method under a Probabilistic Model of Rate Evolution. In: Molecular Biology
           and Evolution, 18(3):352-361, March 2001. doi: 10.1093/oxfordjournals.molbev.a003811.
    """
    node_rates = {}  # maps node v --> rate of v
    edge_rates = {}  # maps v of edge (u,v) --> rate of (u,v)

    for v in tree.preorder():
        if not v.parent:
            # assign factor 1.0 to root (= expected value for all other nodes
            # and edges)
            node_rates[v.label] = 1.0
            edge_rates[v.label] = 1.0
        else:
            var = variance * v.dist
            # ensure that exp. value is equal to parent's rate
            mu = np.log(node_rates[v.parent.label]) - var / 2

            node_rates[v.label] = np.exp(np.random.normal(mu, np.sqrt(var)))

            # edge rate as arithmetic mean of u and v
            edge_rates[v.label] = (node_rates[v.parent.label] + node_rates[v.label]) / 2

    return node_rates, edge_rates

distance_noise

Simulation of noisy distance matrices.

Distances derived from (real-life) gene or protein sequences are burdened with noise. Such data can either be modeled by simulating sequences, or by disturbing the distances specified by a given tree directly. This module implements functions for the latter alternative.

noisy_matrix(orig_matrix: np.ndarray, sd: float, metric_repair: str = 'reject') -> np.ndarray

Disturb a distance matrix (which must be a metric) with random noise.

Parameters:

Name	Type	Description	Default
orig_matrix	ndarray	The 2-dimensional distance matrix.	required
sd	float	Disturbance parameter. Serves as standard deviation of a normal distrubition with mean 1.0.	required
metric_repair	str	Strategy to ensure that the resulting matrix is still a metric and, in particular, satisfies the triangle inequality. The default is 'reject', see [1]. Other available options are 'general' and 'DOMR', see [2].	'reject'

Returns:

Type	Description
ndarray	The disturbed distance matrix.

Raises:

Type	Description
ValueError	If the argument 'metric_repair' is not one of 'reject', 'general', or 'DOMR'.

References

1. P. F. Stadler, M. Geiß, D. Schaller, A. López Sánchez, M. González Laffitte, D. Valdivia, M. Hellmuth, M. Hernández Rosales. From pairs of most similar sequences to phylogenetic best matches. In: Alg. Mol. Biol., 2020, 15, 5. doi: 10.1186/s13015-020-00165-2.
2. A. C. Gilbert, L. Jain. If it ain't broke, don't fix it: Sparse metric repair. In: 2017 55th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 612-619, Monticello, IL, USA, October 2017. IEEE. ISBN 978-1-5386-3266-6. doi: 10.1109/ALLERTON.2017.8262793.

Source code in asymmetree/treeevolve/distance_noise.py

def noisy_matrix(orig_matrix: np.ndarray, sd: float, metric_repair: str = "reject") -> np.ndarray:
    """Disturb a distance matrix (which must be a metric) with random noise.

    Args:
        orig_matrix: The 2-dimensional distance matrix.
        sd: Disturbance parameter. Serves as standard deviation of a normal distrubition with mean
            1.0.
        metric_repair: Strategy to ensure that the resulting matrix is still a metric and, in
            particular, satisfies the triangle inequality. The default is 'reject', see [1]. Other
            available options are 'general' and 'DOMR', see [2].

    Returns:
        The disturbed distance matrix.

    Raises:
        ValueError: If the argument 'metric_repair' is not one of 'reject', 'general', or 'DOMR'.

    References:
        1. P. F. Stadler, M. Geiß, D. Schaller, A. López Sánchez, M. González Laffitte, D.
           Valdivia, M. Hellmuth, M. Hernández Rosales. From pairs of most similar sequences to
           phylogenetic best matches. In: Alg. Mol. Biol., 2020, 15, 5.
           doi: 10.1186/s13015-020-00165-2.
        2. A. C. Gilbert, L. Jain. If it ain't broke, don't fix it: Sparse metric repair.
           In: 2017 55th Annual Allerton Conference on Communication, Control, and Computing
           (Allerton), pages 612-619, Monticello, IL, USA, October 2017. IEEE.
           ISBN 978-1-5386-3266-6. doi: 10.1109/ALLERTON.2017.8262793.
    """
    if metric_repair == "general":
        return _noise_general_metric_repair(orig_matrix, sd)
    elif metric_repair == "DOMR":
        return _noise_metric_repair_DOMR(orig_matrix, sd)
    elif metric_repair == "reject":
        return _noise_reject_method(orig_matrix, sd)

    raise ValueError(f"illegal argument '{metric_repair}'")

convex_linear_comb(D1: np.ndarray, D2: np.ndarray, alpha: float = 0.05, first_only: bool = False) -> np.ndarray | tuple[np.ndarray, np.ndarray]

Convex linear combination of distance matrices.

Returns the convex linear combinations of two distance matrices D1' = (1-alpha) * D1 + alpha * D2 D2' = (1-alpha) * D2 + alpha * D1

The motivation for this function is to simulate non-random noise from a wrong topology, which can be used to test the robustness of tree reconstruction methods.

Parameters:

Name	Type	Description	Default
D1	ndarray	The first distance matrix.	required
D2	ndarray	The second distance matrix.	required
alpha	float	Disturbance parameter, the default is 0.05.	0.05
first_only	bool	If True, only return the disturbed first matrix.	False

Returns:

Type	Description
ndarray | tuple[ndarray, ndarray]	The disturbed distance matrix or matrices as a convex linear combination of the input matrices. If first_only is False, a tuple of the two disturbed matrices is returned. If first_only is True, only the disturbed first matrix is returned.

References

1. P. F. Stadler, M. Geiß, D. Schaller, A. López Sánchez, M. González Laffitte, D. Valdivia, M. Hellmuth, M. Hernández Rosales. From pairs of most similar sequences to phylogenetic best matches. In: Alg. Mol. Biol., 2020, 15, 5. doi: 10.1186/s13015-020-00165-2.

Source code in asymmetree/treeevolve/distance_noise.py

def convex_linear_comb(
    D1: np.ndarray,
    D2: np.ndarray,
    alpha: float = 0.05,
    first_only: bool = False,
) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
    """Convex linear combination of distance matrices.

    Returns the convex linear combinations of two distance matrices
        D1' = (1-alpha) * D1 + alpha * D2
        D2' = (1-alpha) * D2 + alpha * D1

    The motivation for this function is to simulate non-random noise from a wrong topology, which
    can be used to test the robustness of tree reconstruction methods.

    Args:
        D1: The first distance matrix.
        D2: The second distance matrix.
        alpha: Disturbance parameter, the default is 0.05.
        first_only: If True, only return the disturbed first matrix.

    Returns:
        The disturbed distance matrix or matrices as a convex linear combination of the input
            matrices. If `first_only` is False, a tuple of the two disturbed matrices is returned.
            If `first_only` is True, only the disturbed first matrix is returned.

    References:
        1. P. F. Stadler, M. Geiß, D. Schaller, A. López Sánchez, M. González Laffitte, D.
           Valdivia, M. Hellmuth, M. Hernández Rosales. From pairs of most similar sequences to
           phylogenetic best matches. In: Alg. Mol. Biol., 2020, 15, 5.
           doi: 10.1186/s13015-020-00165-2.
    """
    if not first_only:
        if D1.shape == D2.shape:
            return ((1 - alpha) * D1 + alpha * D2, (1 - alpha) * D2 + alpha * D1)

        if D1.shape > D2.shape:
            D1, D2 = D2, D1  # now D1 is smaller
            swap = True
        else:
            swap = False

        N1, N2 = D1.shape[0], D2.shape[0]
        D1_alpha = (1 - alpha) * D1 + alpha * D2[:N1, :N1]
        D2_alpha = (1 - alpha) * D2

        for i in range(N2 // N1):
            for j in range(N2 // N1):
                D2_alpha[i * N1 : i * N1 + N1, j * N1 : j * N1 + N1] += alpha * D1
                if i == (N2 // N1) - 1 and i == j:
                    D2_alpha[i * N1 + N1 : N2, j * N1 + N1 : N2] += (
                        alpha * D1[: (N2 % N1), : (N2 % N1)]
                    )
                if i == (N2 // N1) - 1:
                    D2_alpha[i * N1 + N1 : N2, j * N1 : j * N1 + N1] += alpha * D1[: (N2 % N1), :]
                if j == (N2 // N1) - 1:
                    D2_alpha[i * N1 : i * N1 + N1, j * N1 + N1 : N2] += alpha * D1[:, : (N2 % N1)]

        if swap:
            return D2_alpha, D1_alpha
        else:
            return D1_alpha, D2_alpha

    else:
        if D1.shape == D2.shape:
            return (1 - alpha) * D1 + alpha * D2

        N1, N2 = D1.shape[0], D2.shape[0]
        if N1 < N2:
            return (1 - alpha) * D1 + alpha * D2[:N1, :N1]
        else:
            D1_alpha = (1 - alpha) * D1
            for i in range(N1 // N2):
                for j in range(N1 // N2):
                    D1_alpha[i * N2 : i * N2 + N2, j * N2 : j * N2 + N2] += alpha * D2
                    if i == (N1 // N2) - 1 and i == j:
                        D1_alpha[i * N2 + N2 : N1, j * N2 + N2 : N1] += (
                            alpha * D2[: (N1 % N2), : (N1 % N2)]
                        )
                    if i == (N1 // N2) - 1:
                        D1_alpha[i * N2 + N2 : N1, j * N2 : j * N2 + N2] += (
                            alpha * D2[: (N1 % N2), :]
                        )
                    if j == (N1 // N2) - 1:
                        D1_alpha[i * N2 : i * N2 + N2, j * N2 + N2 : N1] += (
                            alpha * D2[:, : (N1 % N2)]
                        )
            return D1_alpha

wrong_topology_matrix(PGT: Tree) -> np.ndarray

Return a wrong topology matrix by rearranging the edges of a binary tree.

Parameters:

Name	Type	Description	Default
PGT	Tree	The input binary tree, which is typically a pruned gene tree (PGT).	required

Returns:

Type	Description
ndarray	The distance matrix of the wrong topology.

Source code in asymmetree/treeevolve/distance_noise.py

def wrong_topology_matrix(PGT: Tree) -> np.ndarray:
    """Return a wrong topology matrix by rearranging the edges of a binary tree.

    Args:
        PGT: The input binary tree, which is typically a pruned gene tree (PGT).

    Returns:
        The distance matrix of the wrong topology.
    """
    # do not include root
    distances = [v.dist for v in PGT.preorder()][1:]

    # pruned gene tree (PGT) is binary and not planted, hence |E| should be even
    if len(distances) % 2 != 0:
        raise ValueError("The input tree is not binary!")

    random.shuffle(distances)

    random_tree = Tree(TreeNode(label=0, dist=0.0))
    id_counter = 1
    current_leaves = [random_tree.root]

    while distances:
        v = current_leaves.pop(random.randint(0, len(current_leaves) - 1))
        dist1, dist2 = distances.pop(), distances.pop()
        new_child1 = TreeNode(label=id_counter, dist=dist1)
        new_child2 = TreeNode(label=id_counter + 1, dist=dist2)
        v.add_child(new_child1)
        v.add_child(new_child2)
        current_leaves.extend(v.children)
        id_counter += 2

    # implicit random bijection to original tree by random leaf order
    random_leaves = [leaf for leaf in random_tree.leaves()]
    random.shuffle(random_leaves)

    _, D = distance_matrix(random_tree, leaf_order=random_leaves)

    return D

check_metric(matrix: np.ndarray, verbose: bool = False) -> bool

Check whether a given matrix is a metric.

Parameters:

Name	Type	Description	Default
matrix	ndarray	The distance matrix to check.	required
verbose	bool	If True, print the reason for failure if the matrix is not a metric.	False

Returns:

Type	Description
bool	True if the matrix is a metric, False otherwise.

Source code in asymmetree/treeevolve/distance_noise.py

def check_metric(matrix: np.ndarray, verbose: bool = False) -> bool:
    """Check whether a given matrix is a metric.

    Args:
        matrix: The distance matrix to check.
        verbose: If True, print the reason for failure if the matrix is not a metric.

    Returns:
        True if the matrix is a metric, False otherwise.
    """
    N = matrix.shape[0]

    for i in range(N):
        if matrix[i, i] != 0.0:
            if verbose:
                print("Not all diagonal elements zero!")
            return False

    for i in range(N - 1):
        for j in range(i + 1, N):
            if matrix[i, j] != matrix[j, i]:
                if verbose:
                    print("Not symmetrical for", i, j)
                return False

            if abs(matrix[i, j] - np.min(matrix[i, :] + matrix[:, j])) > 1e-8:
                if verbose:
                    print(
                        "Violation of triangle inequality!",
                        i,
                        j,
                        np.min(matrix[i, :] + matrix[:, j]),
                        matrix[i, j],
                    )
                return False

    return True