#! /usr/bin/env python # -*- coding: utf-8 -*- ############################################################################## ## DendroPy Phylogenetic Computing Library. ## ## Copyright 2010-2015 Jeet Sukumaran and Mark T. Holder. ## All rights reserved. ## ## See "LICENSE.rst" for terms and conditions of usage. ## ## If you use this work or any portion thereof in published work, ## please cite it as: ## ## Sukumaran, J. and M. T. Holder. 2010. DendroPy: a Python library ## for phylogenetic computing. Bioinformatics 26: 1569-1571. ## ############################################################################## """ Taxon-to-taxon phylogenetic distances. """ import math import collections import csv from dendropy.calculate import statistics from dendropy.utility import GLOBAL_RNG from dendropy.utility import container from dendropy.utility import error import dendropy class PhylogeneticDistanceMatrix(object): """ Calculates and maintains patristic distance information of taxa on a tree. """ @classmethod def from_tree(cls, tree, *args, **kwargs): """ Creates and returns a |PhylogeneticDistanceMatrix| based on the given tree. Note that this creates a "snapshot" of the current state of the tree. Subsequent changes to the tree will not be reflected in |PhylogeneticDistanceMatrix| instances previously created. Also note that syntactically you may prefer to use:: pdm = tree.phylogenetic_distance_matrix() instead of:: pdm = PhylogeneticDistanceMatrix.from_tree(tree) Parameters ---------- tree : a |Tree| instance The |Tree| from which to get the phylogenetic distances. Returns ------- pdm : A |PhylogeneticDistanceMatrix| instance Examples -------- :: import dendropy tree = dendropy.Tree.get(path="tree.nex", schema="nexus") pdm1 = dendropy.PhylogeneticDistanceMatrix.from_tree(tree) # following is equivalent to above and probably preferred: pdm2 = tree.phylogenetic_distance_matrix() """ pdm = cls(*args, **kwargs) pdm.compile_from_tree(tree=tree) return pdm @classmethod def from_csv(cls, src, taxon_namespace=None, is_allow_new_taxa=None, is_first_row_column_names=True, is_first_column_row_names=True, default_data_type=float, label_transform_fn=None, **csv_reader_kwargs ): """ Instantiates a new PhylogeneticDistanceMatrix instance with data from an external source. Parameters ---------- src : file or file-like Source of data. This is a token delimited-file (e.g., a comma-delimited or tab-delimited file) providing a table which lists taxon labels in both rows and columns. The cells of the table are numeric (typically real) values that indicate the distance between the taxa of the current row and column. Note that *only* the upper right section of the table is considered. The diagonals values are typically zeroes and, in either case, ignored along with the lower diagonal. Despite being ignored by the PhylogeneticDistanceMatrix object, the values are parsed by the underlying reader and thus have to be valid numerical values. taxon_namespace : |TaxonNamespace| instance The taxon namespace with which to manage taxa. If this has not already been pre-populated with the taxon names, then ``is_allow_new_taxa`` should be set to |True|. is_allow_new_taxa : bool If |False|: we do *not* expect to encounter any new taxa in the data file, and it is an error if we do. If |True|: we do expect to encounter new taxa in the data file. The default value of this depends on the value passed to ``taxon_namespace``. If ``taxon_namespace`` is ``None`` or an empty |TaxonNamespace| instance, then unless explicitly set to |False|, ``is_allow_new_taxa`` will default to |True|: allowing of creation of new taxa corresponding to labels found in the data source. On the other hand, if ``taxon_namespace`` is not None and its value is a |TaxonNamespace| instance with at least one taxon, unless explicitly set to |True|, ``is_allow_new_taxa`` will default to |False|, and it will be an error if taxon labels are found in the data source that do not correspond (exactly) to |Taxon| objects defined in the taxon namespace. This is to err on the side of caution, to avoid (or rather, highlight) problems due to incorrect or mismatching labels between the data source and the current taxon namespace. is_first_row_column_names : bool By default |True|: assumes that first row lists the taxon names. Set to |False| if there is no header row. is_first_column_row_names : bool By default |True|: assumes that first column lists the taxon names. Set to |False| if there is now row name column. label_transform_fn : function object If not None, this should be a function object that takes a string as an argument and returns another string. This function will be applied to row and column labels before they are matched to taxon labels in the |TaxonNamespace| instance given by ``taxon_namespace``. \*\*csv_reader_kwargs : keyword arguments This arguments will be passed to the underlying CSV reader. The most important one is probably 'delimiter'. Returns ------- pdm : A |PhylogeneticDistanceMatrix| instance Examples -------- :: import dendropy pdm1 = dendropy.PhylogeneticDistanceMatrix.from_csv( src=open("data.csv"), delimiter=",") pdm2 = dendropy.PhylogeneticDistanceMatrix.from_csv( src=open("data.tsv"), delimiter="\t") """ if taxon_namespace is None: taxon_namespace = dendropy.TaxonNamespace() if len(taxon_namespace) == 0 and is_allow_new_taxa is None: is_allow_new_taxa = True old_taxon_namespace_mutability = taxon_namespace.is_mutable taxon_namespace.is_mutable = is_allow_new_taxa data_table = container.DataTable.from_csv( src, is_first_row_column_names=is_first_row_column_names, is_first_column_row_names=is_first_column_row_names, default_data_type=default_data_type, label_transform_fn=label_transform_fn, **csv_reader_kwargs ) if not is_first_row_column_names and not is_first_column_row_names: distances = {} seen_row_labels = set() seen_row_taxa_labels = set() for i1, t1_label in enumerate(data_table.row_name_iter()): if len(taxon_namespace) <= i1: t1 = taxon_namespace.require_taxon(label=t1_label) else: t1 = taxon_namespace[i1] seen_row_labels.add(t1_label) assert t1.label not in seen_row_taxa_labels seen_row_taxa_labels.add(t1.label) distances[t1] = {} for i2, t2_label in enumerate(data_table.column_name_iter()): if t2_label in seen_row_labels: continue if len(taxon_namespace) <= i2: t2 = taxon_namespace.require_taxon(label=t2_label) else: t2 = taxon_namespace[i2] distances[t1][t2] = data_table[t1_label, t2_label] else: distances = {} seen_taxa = set() taxa = [] if is_first_column_row_names: name_iter = data_table.row_name_iter() else: name_iter = data_table.column_name_iter() for label in name_iter: t1 = taxon_namespace.require_taxon(label=label) assert t1 not in seen_taxa seen_taxa.add(t1) taxa.append(t1) seen_row_taxa = set() for i1, t1 in enumerate(taxa): assert t1 in seen_taxa seen_row_taxa.add(t1) distances[t1] = {} for i2, t2 in enumerate(taxa): if t2 in seen_row_taxa: continue distances[t1][t2] = data_table[i1, i2] # else: # raise NotImplementedError() taxon_namespace.is_mutable = old_taxon_namespace_mutability pdm = cls() pdm.compile_from_dict( distances=distances, taxon_namespace=taxon_namespace) return pdm def __init__(self, is_store_path_edges=False): self.clear() self.is_store_path_edges = is_store_path_edges def clear(self): self.taxon_namespace = None self._mapped_taxa = set() self._all_distinct_mapped_taxa_pairs = set() self._tree_length = None self._num_edges = None self._taxon_phylogenetic_distances = {} self._taxon_phylogenetic_path_steps = {} self._taxon_phylogenetic_path_edges = {} self._mrca = {} def compile_from_tree(self, tree): """ Calculates the distances. Note that the path length (in number of steps) between taxa that span the root will be off by one if the tree is unrooted. """ self.clear() self.taxon_namespace = tree.taxon_namespace # for i1, t1 in enumerate(self.taxon_namespace): # self._taxon_phylogenetic_distances[t1] = {} # self._taxon_phylogenetic_path_steps[t1] = {} # self._mrca[t1] = {} self._tree_length = 0.0 self._num_edges = 0 if self.is_store_path_edges: default_pedges = [] else: default_pedges = None for node in tree.postorder_node_iter(): try: self._tree_length += node.edge.length except TypeError: # None for edge length pass self._num_edges += 1 children = node.child_nodes() if len(children) == 0: node.desc_paths = {node : (0,0, default_pedges)} else: node.desc_paths = {} for cidx1, c1 in enumerate(children): for desc1, (desc1_plen, desc1_psteps, desc1_pedges) in c1.desc_paths.items(): if c1.edge_length is None: c1_edge_length = 0.0 else: c1_edge_length = c1.edge.length if self.is_store_path_edges: pedges = list(desc1_pedges + [c1.edge]) else: pedges = default_pedges node.desc_paths[desc1] = (desc1_plen + c1_edge_length, desc1_psteps + 1, pedges) assert desc1.taxon is not None if desc1.taxon not in self._taxon_phylogenetic_distances: self._mapped_taxa.add(desc1.taxon) self._taxon_phylogenetic_distances[desc1.taxon] = {} self._taxon_phylogenetic_distances[desc1.taxon][desc1.taxon] = 0.0 self._taxon_phylogenetic_path_steps[desc1.taxon] = {} self._taxon_phylogenetic_path_steps[desc1.taxon][desc1.taxon] = 0 if self.is_store_path_edges: self._taxon_phylogenetic_path_edges[desc1.taxon] = {} self._taxon_phylogenetic_path_edges[desc1.taxon][desc1.taxon] = [] self._mrca[desc1.taxon] = {desc1.taxon: desc1} for c2 in children[cidx1+1:]: for desc2, (desc2_plen, desc2_psteps, desc2_pedges) in c2.desc_paths.items(): self._mapped_taxa.add(desc2.taxon) self._mrca[desc1.taxon][desc2.taxon] = c1.parent_node # self._all_distinct_mapped_taxa_pairs.add( tuple([desc1.taxon, desc2.taxon]) ) self._all_distinct_mapped_taxa_pairs.add( frozenset([desc1.taxon, desc2.taxon]) ) if c2.edge_length is None: c2_edge_length = 0.0 else: c2_edge_length = c2.edge.length pat_dist = node.desc_paths[desc1][0] + desc2_plen + c2_edge_length self._taxon_phylogenetic_distances[desc1.taxon][desc2.taxon] = pat_dist path_steps = node.desc_paths[desc1][1] + desc2_psteps + 1 self._taxon_phylogenetic_path_steps[desc1.taxon][desc2.taxon] = path_steps if self.is_store_path_edges: pedges = tuple(node.desc_paths[desc1][2] + [c2.edge] + desc2_pedges[::-1]) self._taxon_phylogenetic_path_edges[desc1.taxon][desc2.taxon] = pedges del(c1.desc_paths) self._mirror_lookups() # assert self._tree_length == tree.length() def compile_from_dict(self, distances, taxon_namespace): self.clear() self.taxon_namespace = taxon_namespace for t1 in distances: self._mapped_taxa.add(t1) self._taxon_phylogenetic_distances[t1] = {} for t2 in distances[t1]: self._taxon_phylogenetic_distances[t1][t2] = distances[t1][t2] self._mirror_lookups() def _mirror_lookups(self): for ddata in ( self._taxon_phylogenetic_distances, self._taxon_phylogenetic_path_steps, self._mrca, ): for taxon1 in ddata: for taxon2 in ddata[taxon1]: # assert taxon1 is not taxon2 if taxon2 not in ddata: ddata[taxon2] = {} ddata[taxon2][taxon1] = ddata[taxon1][taxon2] for taxon1 in self._taxon_phylogenetic_path_edges: for taxon2 in self._taxon_phylogenetic_path_edges: if taxon2 not in self._taxon_phylogenetic_path_edges: self._taxon_phylogenetic_path_edges[taxon2][taxon1] = {} if taxon1 not in self._taxon_phylogenetic_path_edges[taxon2]: self._taxon_phylogenetic_path_edges[taxon2][taxon1] = tuple(reversed(self._taxon_phylogenetic_path_edges[taxon1][taxon2])) def __eq__(self, o): if self.taxon_namespace is not o.taxon_namespace: return False return (True and (self._mapped_taxa == o._mapped_taxa) and (self._all_distinct_mapped_taxa_pairs == o._all_distinct_mapped_taxa_pairs) and (self._taxon_phylogenetic_distances == o._taxon_phylogenetic_distances) and (self._taxon_phylogenetic_path_steps == o._taxon_phylogenetic_path_steps) and (self._taxon_phylogenetic_path_edges == o._taxon_phylogenetic_path_edges) and (self._mrca == o._mrca) and (self._tree_length == o._tree_length) and (self._num_edges == o._num_edges) ) def __hash__(self): return id(self) def __call__(self, taxon1, taxon2): return self.patristic_distance(taxon1, taxon2) def __copy__(self): return self.clone() def __iter__(self): for taxon in self._taxon_phylogenetic_distances: yield taxon def clone(self): o = self.__class__() o.taxon_namespace = self.taxon_namespace o._mapped_taxa = set(self._mapped_taxa) o._all_distinct_mapped_taxa_pairs = set(self._all_distinct_mapped_taxa_pairs) o._tree_length = self._tree_length o._num_edges = self._num_edges for src, dest in ( (self._taxon_phylogenetic_distances, o._taxon_phylogenetic_distances,), (self._taxon_phylogenetic_path_steps, o._taxon_phylogenetic_path_steps,), (self._taxon_phylogenetic_path_edges, o._taxon_phylogenetic_path_edges,), (self._mrca, o._mrca,), ): for t1 in src: dest[t1] = {} for t2 in src[t1]: dest[t1][t2] = src[t1][t2] return o def mrca(self, taxon1, taxon2): """ Returns MRCA of two taxon objects. """ return self._mrca[taxon1][taxon2] def distance(self, taxon1, taxon2, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Returns distance between taxon1 and taxon2. """ if is_weighted_edge_distances: return self.patristic_distance(taxon1, taxon2, is_normalize_by_tree_size=is_normalize_by_tree_size) else: return self.path_edge_count(taxon1, taxon2, is_normalize_by_tree_size=is_normalize_by_tree_size) def patristic_distance(self, taxon1, taxon2, is_normalize_by_tree_size=False): """ Returns patristic distance between two taxon objects. """ if taxon1 is taxon2: return 0.0 d = self._taxon_phylogenetic_distances[taxon1][taxon2] if is_normalize_by_tree_size: return d / self._tree_length else: return d def path_edge_count(self, taxon1, taxon2, is_normalize_by_tree_size=False): """ Returns the number of edges between two taxon objects. """ if taxon1 is taxon2: return 0 d = self._taxon_phylogenetic_path_steps[taxon1][taxon2] if is_normalize_by_tree_size: return float(d) / self._num_edges else: return d def path_edges(self, taxon1, taxon2): """ Returns the edges between two taxon objects. """ return self._taxon_phylogenetic_path_edges[taxon1][taxon2] def distances(self, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Returns list of patristic distances. """ dmatrix, normalization_factor = self._get_distance_matrix_and_normalization_factor( is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size, ) results = [] for t1, t2 in self._all_distinct_mapped_taxa_pairs: results.append(dmatrix[t1][t2]/normalization_factor) return results def max_pairwise_distance_taxa(self, is_weighted_edge_distances=True): if is_weighted_edge_distances: dists = self._taxon_phylogenetic_distances else: dists = self._taxon_phylogenetic_path_steps max_dist = None max_dist_taxa = None for t1, t2 in self._all_distinct_mapped_taxa_pairs: pat_dist = dists[t1][t2] if max_dist is None or pat_dist > max_dist: max_dist = pat_dist max_dist_taxa = (t1, t2) return max_dist_taxa def sum_of_distances(self, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Returns sum of patristic distances on tree. """ return sum(self.distances(is_weighted_edge_distances=is_weighted_edge_distances,is_normalize_by_tree_size=is_normalize_by_tree_size)) def taxon_iter(self, filter_fn=None): """ Iterates over taxa in matrix. Note that this could be a subset of the taxa in the associated taxon namespace. """ for t1 in self._mapped_taxa: if not filter_fn or filter_fn(t1): yield t1 def distinct_taxon_pair_iter(self, filter_fn=None): """ Iterates over all distinct pairs of taxa in matrix. """ for t1, t2 in self._all_distinct_mapped_taxa_pairs: if not filter_fn or (filter_fn(t1) and filter_fn(t2)): yield t1, t2 def mean_pairwise_distance(self, filter_fn=None, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Calculates the phylogenetic ecology statistic "MPD"[1,2] for the tree (only considering taxa for which ``filter_fn`` returns True when applied if ``filter_fn`` is specified). The mean pairwise distance (mpd) is given by: .. math:: mpd = \\frac{ \\sum_{i}^{n} \\sum_{j}^{n} \\delta_{i,j} }{n \\choose 2}, where :math:`i \\neq j`, :math:`\\delta_{i,j}` is the phylogenetic distance between species :math:`i` and :math:`j`, and :math:`n` is the number of species in the sample. Parameters ---------- filter_fn : function object or None If |None|, then all leaves will be considered. Otherwise should be a function object that takes a Taxon instance as an argument and returns |True| if it is to be included in the calculation or |False| otherwise. In trees sampled from multiple communites, ``filter_fn`` can be used to restrict the calculation to only one community based on some criteria. is_weighted_edge_distances : bool If |True| then the edge-weighted distance, i.e., considering edge lengths, is returned. Otherwise the the path steps or the number of edges rather then the sum of is_weighted_edge_distances edges, connecting two taxa is considered. is_normalize_by_tree_size : bool If |True| then the results are normalized by the total tree length or steps/edges (depending on whether edge-weighted or unweighted distances are used, respectively). Otherwise, raw distances are used. Returns ------- mpd : float The Mean Pairwise Distance (MPD) statistic for the daata. Examples -------- :: import dendropy tree = dendropy.Tree.get(path="data.nex", schema="nexus") pdm = dendropy.PhylogeneticDistanceMatrix(tree) # consider all tips mpd1 = pdm.mean_pairwise_distance() # only tips within the same community, based on the node annotation # "community" mpds_by_community = {} for community_label in ("1", "2", "3",): filter_fn = lambda x: x.annotations["community"] == community_label mpd = pdm.mean_pairwise_distance(filter_fn=filter_fn) mpds_by_community[community_label] = mpd References ---------- [1] Webb, C.O. 2000. Exploring the phylogenetic structure of ecological communities: An example for rainforest trees. The American Naturalist 156: 145-155. [2] Swenson, N.G. Functional and Phylogenetic Ecology in R. """ comparison_regime = self.distinct_taxon_pair_iter(filter_fn=filter_fn) return self._calculate_mean_pairwise_distance( comparison_regime=comparison_regime, is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size, ) def mean_nearest_taxon_distance(self, filter_fn=None, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Calculates the phylogenetic ecology statistic "MNTD"[1,2] for the tree (only considering taxa for which ``filter_fn`` returns True when applied if ``filter_fn`` is specified). The mean nearest taxon distance (mntd) is given by: .. math:: mntd = \\frac{ \\sum_{i}^{n} min(\\delta_{i,j}) }{n}, where :math:`i \\neq j`, :math:`\\delta_{i,j}` is the phylogenetic distance between species :math:`i` and :math:`j`, and :math:`n` is the number of species in the sample. Parameters ---------- filter_fn : function object or None If |None|, then all leaves will be considered. Otherwise should be a function object that takes a Taxon instance as an argument and returns |True| if it is to be included in the calculation or |False| otherwise. In trees sampled from multiple communites, ``filter_fn`` can be used to restrict the calculation to only one community based on some criteria. is_weighted_edge_distances : bool If |True| then the edge-weighted distance, i.e., considering edge lengths, is returned. Otherwise the the path steps or the number of edges rather then the sum of is_weighted_edge_distances edges, connecting two taxa is considered. is_normalize_by_tree_size : bool If |True| then the results are normalized by the total tree length or steps/edges (depending on whether edge-weighted or unweighted distances are used, respectively). Otherwise, raw distances are used. Returns ------- mntd : float The Mean Nearest Taxon Distance (MNTD) statistic for the daata. Examples -------- :: import dendropy tree = dendropy.Tree.get(path="data.nex", schema="nexus") pdm = dendropy.PhylogeneticDistanceMatrix(tree) # consider all tips mntd = pdm.mean_nearest_taxon_distance() # only tips within the same community, based on the node annotation # "community" mntds_by_community = {} for community_label in ("1", "2", "3",): filter_fn = lambda x: x.annotations["community"] == community_label mntd = pdm.mean_pairwise_distance(filter_fn=filter_fn) mntds_by_community[community_label] = mntd References ---------- [1] Webb, C.O. 2000. Exploring the phylogenetic structure of ecological communities: An example for rainforest trees. The American Naturalist 156: 145-155. [2] Swenson, N.G. Functional and Phylogenetic Ecology in R. """ comparison_regime = self._get_taxon_to_all_other_taxa_comparisons(filter_fn=filter_fn) return self._calculate_mean_nearest_taxon_distance( comparison_regime=comparison_regime, is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size,) def standardized_effect_size_mean_pairwise_distance(self, assemblage_memberships, num_randomization_replicates=1000, is_weighted_edge_distances=True, is_normalize_by_tree_size=False, is_skip_single_taxon_assemblages=False, null_model_type="taxa.label", rng=None): """ Returns the standardized effect size value for the MPD statistic under a null model under various community compositions. The S.E.S. is given by: .. math:: SES(statistic) = \\frac{observed - mean(model_{null})}{sd(model_{null})} This removes any bias associated with the decrease in variance in the MPD statistic value as species richness increases to the point where communities become saturated. Equivalent to -1 times the Nearest Relative Index (NRI) when using phylogenetic distances. In contrast to the function calculating the non-standardized effect size version of this statistic, which uses filter function to specify the subset of taxa to be considerd, here a collection of (multiple) sets of taxa constituting a community is specified. This to allow calculation of the null model statistic across all community sets for each randomization replicate. Parameters ---------- assemblage_memberships : iterable of iterable of |Taxon| objects A collection of collections, e.g. a list of sets, with the elements of each set being |Taxon| instances. Each set specifies the composition of a community. The standardized effect size of this statistic will be calculated for each community as specified by a set of |Taxon| instances. num_randomization_replicates : int Number of randomization replicates. is_weighted_edge_distances: bool If ``True`` then edge lengths will be considered for distances. Otherwise, just the number of edges. Returns ------- r : list of results A list of results, with each result corresponding to a community set given in ``assemblage_memberships``. Each result consists of a named tuple with the following elements: - obs : the observed value of the statistic - null_model_mean : the mean value of the statistic under the null model - null_model_sd : the standard deviation of the statistic under the null model - z : the standardized effect value of the statistic (= SES as defined in [1] above) - p : the p-value of the observed value of the Examples -------- :: import dendropy tree = dendropy.Tree.get_from_path( src="data/community.tree.newick", schema="newick", rooting="force-rooted") pdm = tree.phylogenetic_distance_matrix() assemblage_membership_definitions = pdm.assemblage_membership_definitions_from_csv("data/comm1.csv") results = pdm.standardized_effect_size_mean_pairwise_distance(assemblage_memberships=assemblage_membership_definitions.values()) print(results) """ if assemblage_memberships is None: assemblage_memberships = [ set(self._mapped_taxa) ] comparison_regimes = [] for idx, assemblage_membership in enumerate(assemblage_memberships): if len(assemblage_membership) == 1: if is_skip_single_taxon_assemblages: continue else: raise error.SingleTaxonAssemblageException("{}: {}".format(idx, assemblage_membership)) filter_fn = lambda taxon: taxon in assemblage_membership comparison_regime = list(self.distinct_taxon_pair_iter(filter_fn=filter_fn)) comparison_regimes.append(comparison_regime) results = self._calculate_standardized_effect_size( statisticf_name="_calculate_mean_pairwise_distance", is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size, comparison_regimes=comparison_regimes, null_model_type=null_model_type, num_randomization_replicates=num_randomization_replicates, rng=rng) return results def standardized_effect_size_mean_nearest_taxon_distance(self, assemblage_memberships, num_randomization_replicates=1000, is_weighted_edge_distances=True, is_normalize_by_tree_size=False, is_skip_single_taxon_assemblages=False, null_model_type="taxa.label", rng=None): """ Returns the standardized effect size value for the MNTD statistic under a null model under various community compositions. The S.E.S. is given by: .. math:: SES(statistic) = \\frac{observed - mean(model_{null})}{sd(model_{null})} This removes any bias associated with the decrease in variance in the MPD statistic value as species richness increases to the point where communities become saturated. Equivalent to -1 times the Nearest Taxon Index when using phylogenetic distances. In contrast to the function calculating the non-standardized effect size version of this statistic, which uses filter function to specify the subset of taxa to be considerd, here a collection of (multiple) sets of taxa constituting a community is specified. This to allow calculation of the null model statistic across all community sets for each randomization replicate. Parameters ---------- assemblage_memberships : iterable of iterable of |Taxon| objects A collection of collections, e.g. a list of sets, with the elements of each set being |Taxon| instances. Each set specifies the composition of a community. The standardized effect size of this statistic will be calculated for each community as specified by a set of |Taxon| instances. num_randomization_replicates : int Number of randomization replicates. is_weighted_edge_distances: bool If ``True`` then edge lengths will be considered for distances. Otherwise, just the number of edges. Returns ------- r : list of results A list of results, with each result corresponding to a community set given in ``assemblage_memberships``. Each result consists of a named tuple with the following elements: - obs : the observed value of the statistic - null_model_mean : the mean value of the statistic under the null model - null_model_sd : the standard deviation of the statistic under the null model - z : the standardized effect value of the statistic (= SES as defined in [1] above) - p : the p-value of the observed value of the statistic under the null model. Examples -------- :: import dendropy tree = dendropy.Tree.get_from_path( src="data/community.tree.newick", schema="newick", rooting="force-rooted") pdm = dendropy.PhylogeneticDistanceMatrix.from_tree(tree) assemblage_memberships = pdm.assemblage_membership_definitions_from_csv("data/comm1.csv") results = pdm.standardized_effect_size_mean_nearest_taxon_distance(assemblage_memberships=assemblage_memberships) print(results) """ if assemblage_memberships is None: assemblage_memberships = [ set(self._mapped_taxa) ] comparison_regimes = [] for idx, assemblage_membership in enumerate(assemblage_memberships): if len(assemblage_membership) == 1: if is_skip_single_taxon_assemblages: continue else: raise error.SingleTaxonAssemblageException("{}: {}".format(idx, assemblage_membership)) filter_fn = lambda taxon: taxon in assemblage_membership comparison_regime = self._get_taxon_to_all_other_taxa_comparisons(filter_fn=filter_fn) comparison_regimes.append(comparison_regime) results = self._calculate_standardized_effect_size( statisticf_name="_calculate_mean_nearest_taxon_distance", comparison_regimes=comparison_regimes, is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size, null_model_type=null_model_type, num_randomization_replicates=num_randomization_replicates, rng=rng) return results def shuffle_taxa(self, is_shuffle_phylogenetic_distances=True, is_shuffle_phylogenetic_path_steps=True, is_shuffle_mrca=True, rng=None): """ Randomly shuffles taxa in-situ. """ if rng is None: rng = GLOBAL_RNG reordered_taxa = list(self._mapped_taxa) rng.shuffle(reordered_taxa) current_to_shuffled_taxon_map = dict(zip(self._mapped_taxa, reordered_taxa)) to_shuffle = [] if is_shuffle_phylogenetic_distances: to_shuffle.append("_taxon_phylogenetic_distances") if is_shuffle_phylogenetic_path_steps: to_shuffle.append("_taxon_phylogenetic_path_steps") if is_shuffle_mrca: to_shuffle.append("_mrca") for attr_name in to_shuffle: src = getattr(self, attr_name) dest = {} ## 5m8.076s # for t1, t2 in self._all_distinct_mapped_taxa_pairs: # x1 = current_to_shuffled_taxon_map[t1] # x2 = current_to_shuffled_taxon_map[t2] # d = src[t1][t2] # try: # dest[x1][x2] = d # except KeyError: # dest[x1] = {x2: d} # if t1 in src[t1]: # dest[x1][x1] = src[t1][t1] # try: # dest[x2][x1] = d # except KeyError: # dest[x2] = {x1: d} # if t2 in src[t2]: # dest[x2][x2] = src[t2][t2] # setattr(self, attr_name, dest) # 4m48.025s for t1 in src: x1 = current_to_shuffled_taxon_map[t1] dest[x1] = {} for t2 in src[t1]: x2 = current_to_shuffled_taxon_map[t2] dest[x1][x2] = src[t1][t2] setattr(self, attr_name, dest) return current_to_shuffled_taxon_map def nj_tree(self, is_weighted_edge_distances=True, tree_factory=None, ): """ Returns an Neighbor-Joining (NJ) tree based on the distances in the matrix. Calculates and returns a tree under the Neighbor-Joining algorithm of Saitou and Nei (1987) for the data in the matrix. Parameters ---------- is_weighted_edge_distances: bool If ``True`` then edge lengths will be considered for distances. Otherwise, just the number of edges. Returns ------- t : |Tree| A |Tree| instance corresponding to the Neighbor-Joining (NJ) tree for this data. Examples -------- :: import dendropy # Read data from a CSV file into a PhylogeneticDistanceMatrix # object with open("distance_matrix.csv") as src: pdm = dendropy.PhylogeneticDistanceMatrix.from_csv( src, is_first_row_column_names=True, is_first_column_row_names=True, is_allow_new_taxa=True, delimiter=",", ) # Calculate the tree nj_tree = pdm.nj_tree() # Print it print(nj_tree.as_string("nexus")) References ---------- Saitou, N. and Nei, M. (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution, 4: 406-425. """ if is_weighted_edge_distances: original_dmatrix = self._taxon_phylogenetic_distances else: original_dmatrix = self._taxon_phylogenetic_path_steps if tree_factory is None: tree_factory = dendropy.Tree tree = tree_factory(taxon_namespace=self.taxon_namespace) tree.is_rooted = False # initialize node pool node_pool = [] for t1 in self._mapped_taxa: nd = tree.node_factory() nd.taxon = t1 nd._nj_distances = {} node_pool.append(nd) # initialize factor n = len(self._mapped_taxa) # cache calculations for nd1 in node_pool: nd1._nj_xsub = 0.0 for nd2 in node_pool: if nd1 is nd2: continue d = original_dmatrix[nd1.taxon][nd2.taxon] nd1._nj_distances[nd2] = d nd1._nj_xsub += d while n > 1: # calculate the Q-matrix min_q = None nodes_to_join = None for idx1, nd1 in enumerate(node_pool[:-1]): for idx2, nd2 in enumerate(node_pool[idx1+1:]): v1 = (n - 2) * nd1._nj_distances[nd2] qvalue = v1 - nd1._nj_xsub - nd2._nj_xsub if min_q is None or qvalue < min_q: min_q = qvalue nodes_to_join = (nd1, nd2) # create the new node new_node = tree.node_factory() # attach it to the tree for node_to_join in nodes_to_join: new_node.add_child(node_to_join) node_pool.remove(node_to_join) # calculate the distances for the new node new_node._nj_distances = {} new_node._nj_xsub = 0.0 for node in node_pool: # actual node-to-node distances v1 = 0.0 for node_to_join in nodes_to_join: v1 += node._nj_distances[node_to_join] v3 = nodes_to_join[0]._nj_distances[nodes_to_join[1]] dist = 0.5 * (v1 - v3) new_node._nj_distances[node] = dist node._nj_distances[new_node] = dist # Adjust/recalculate the values needed for the Q-matrix # calculations new_node._nj_xsub += dist node._nj_xsub += dist for node_to_join in nodes_to_join: node._nj_xsub -= node_to_join._nj_distances[node] # calculate the branch lengths if n > 2: v1 = 0.5 * nodes_to_join[0]._nj_distances[nodes_to_join[1]] v4 = 1.0/(2*(n-2)) * (nodes_to_join[0]._nj_xsub - nodes_to_join[1]._nj_xsub) delta_f = v1 + v4 delta_g = nodes_to_join[0]._nj_distances[nodes_to_join[1]] - delta_f nodes_to_join[0].edge.length = delta_f nodes_to_join[1].edge.length = delta_g else: d = nodes_to_join[0]._nj_distances[nodes_to_join[1]] nodes_to_join[0].edge.length = d / 2 nodes_to_join[1].edge.length = d / 2 # clean up for node_to_join in nodes_to_join: del node_to_join._nj_distances del node_to_join._nj_xsub # add the new node to the pool of nodes node_pool.append(new_node) # adjust count n -= 1 tree.seed_node = node_pool[0] del tree.seed_node._nj_distances del tree.seed_node._nj_xsub return tree def upgma_tree(self, is_weighted_edge_distances=True, tree_factory=None, ): """ Returns an Unweighted Pair Group Method with Arithmetic Mean (UPGMA) tree based on the distances in the matrix. Parameters ---------- is_weighted_edge_distances: bool If ``True`` then edge lengths will be considered for distances. Otherwise, just the number of edges. Returns ------- t : |Tree| A |Tree| instance corresponding to the UPGMA tree for this data. Examples -------- :: import dendropy # Read data from a CSV file into a PhylogeneticDistanceMatrix # object with open("distance_matrix.csv") as src: pdm = dendropy.PhylogeneticDistanceMatrix.from_csv( src, is_first_row_column_names=True, is_first_column_row_names=True, is_allow_new_taxa=True, delimiter=",", ) # Calculate the tree upgma_tree = pdm.upgma_tree() # Print it print(upgma_tree.as_string("nexus")) """ if is_weighted_edge_distances: original_dmatrix = self._taxon_phylogenetic_distances else: original_dmatrix = self._taxon_phylogenetic_path_steps if tree_factory is None: tree_factory = dendropy.Tree tree = tree_factory(taxon_namespace=self.taxon_namespace) tree.is_rooted = True node_pool = [] for t1 in self._mapped_taxa: nd = tree.node_factory() nd.taxon = t1 nd._upgma_cluster = set([nd]) nd._upgma_distance_from_tip = 0.0 nd._upgma_distances = {} node_pool.append(nd) for idx1, nd1 in enumerate(node_pool[:-1]): for idx2, nd2 in enumerate(node_pool[idx1+1:]): d = original_dmatrix[nd1.taxon][nd2.taxon] nd1._upgma_distances[nd2] = d nd2._upgma_distances[nd1] = d while len(node_pool) > 1: min_distance = None nodes_to_join = None for idx1, nd1 in enumerate(node_pool[:-1]): for idx2, nd2 in enumerate(node_pool[idx1+1:]): d = nd1._upgma_distances[nd2] if min_distance is None or d < min_distance: nodes_to_join = (nd1, nd2) min_distance = d new_node = tree.node_factory() new_node._upgma_cluster = set() new_node._upgma_distances = {} elen = min_distance / 2.0 for node_to_join in nodes_to_join: new_node.add_child(node_to_join) new_node._upgma_cluster.update(node_to_join._upgma_cluster) node_to_join.edge.length = elen - node_to_join._upgma_distance_from_tip node_pool.remove(node_to_join) new_node._upgma_distance_from_tip = nodes_to_join[0].edge.length + nodes_to_join[0]._upgma_distance_from_tip for idx1, nd1 in enumerate(node_pool): d1 = 0.0 count = 0.0 for node_to_join in nodes_to_join: d2 = node_to_join._upgma_distances[nd1] xc = len(node_to_join._upgma_cluster) d1 += (d2 * xc) count += xc d = d1 / count nd1._upgma_distances[new_node] = d new_node._upgma_distances[nd1] = d for node_to_join in nodes_to_join: del node_to_join._upgma_cluster del node_to_join._upgma_distance_from_tip del node_to_join._upgma_distances node_pool.append(new_node) tree.seed_node = node_pool[0] del tree.seed_node._upgma_cluster del tree.seed_node._upgma_distance_from_tip del tree.seed_node._upgma_distances return tree def as_data_table(self, is_weighted_edge_distances=True): """ Returns this as a table. """ if is_weighted_edge_distances: df = self.patristic_distance else: df = self.path_edge_count dt = container.DataTable() for t1 in self._mapped_taxa: dt.add_row(row_name=t1.label) dt.add_column(column_name=t1.label) for t1 in self._mapped_taxa: for t2 in self._mapped_taxa: dt[t1.label, t2.label] = df(t1, t2) return dt def write_csv(self, out, is_first_row_column_names=True, is_first_column_row_names=True, is_weighted_edge_distances=True, is_normalize_by_tree_size=True, label_transform_fn=None, **csv_writer_kwargs ): if isinstance(out, str): dest = open(out, "w") else: dest = out if label_transform_fn is None: label_transform_fn = lambda x: x dmatrix, normalization_factor = self._get_distance_matrix_and_normalization_factor( is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size, ) if "delimiter" not in csv_writer_kwargs: csv_writer_kwargs["delimiter"] = "," writer = csv.writer(dest, csv_writer_kwargs) if is_first_row_column_names: row = [] if is_first_column_row_names: row.append("") for taxon in self._mapped_taxa: row.append(label_transform_fn(taxon.label)) writer.writerow(row) # dest.write(delimiter.join(row)) # dest.write("\n") for taxon1 in self._mapped_taxa: row = [] if is_first_column_row_names: row.append(label_transform_fn(taxon1.label)) for taxon2 in self._mapped_taxa: d = dmatrix[taxon1][taxon2] / normalization_factor row.append("{}".format(d)) writer.writerow(row) # dest.write(delimiter.join(row)) # dest.write("\n") def assemblage_membership_definitions_from_csv( self, src, default_data_type=float, **csv_reader_kwargs): """ Convenience method to return list of community sets from a delimited file that lists taxon (labels) in columns and community presence/absences or abundances in rows. """ if isinstance(src, str): with open(src) as srcf: data_table = container.DataTable.from_csv( src, default_data_type=default_data_type, **csv_reader_kwargs ) else: data_table = container.DataTable.from_csv( src, default_data_type=default_data_type, **csv_reader_kwargs ) mapped_taxon_labels = set([taxon.label for taxon in self.taxon_iter()]) for column_name in data_table.column_name_iter(): assert column_name in mapped_taxon_labels assemblage_memberships = collections.OrderedDict() for row_name in data_table.row_name_iter(): assemblage_membership = set() for taxon in self.taxon_iter(): if data_table[row_name, taxon.label] > 0: assemblage_membership.add(taxon) assemblage_memberships[row_name] = assemblage_membership return assemblage_memberships def _get_taxon_to_all_other_taxa_comparisons(self, filter_fn=None): permutations = collections.defaultdict(list) for taxon1 in self._mapped_taxa: # permutations[taxon1] = [] if filter_fn and not filter_fn(taxon1): continue for taxon2 in self._mapped_taxa: if taxon1 is taxon2: continue if filter_fn and not filter_fn(taxon2): continue permutations[taxon1].append(taxon2) return permutations def _get_distance_matrix_and_normalization_factor(self, is_weighted_edge_distances, is_normalize_by_tree_size): if is_weighted_edge_distances: dmatrix = self._taxon_phylogenetic_distances if is_normalize_by_tree_size: normalization_factor = self._tree_length else: normalization_factor = 1.0 else: dmatrix = self._taxon_phylogenetic_path_steps if is_normalize_by_tree_size: normalization_factor = float(self._num_edges) else: normalization_factor = 1.0 return dmatrix, normalization_factor def _calculate_mean_pairwise_distance(self, comparison_regime, is_weighted_edge_distances, is_normalize_by_tree_size): dmatrix, normalization_factor = self._get_distance_matrix_and_normalization_factor( is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size,) distances = [] for taxon1, taxon2 in comparison_regime: distances.append(dmatrix[taxon1][taxon2]) if distances: return (sum(distances) / normalization_factor) / (len(distances) * 1.0) else: raise error.NullAssemblageException("No taxa in assemblage") def _calculate_mean_nearest_taxon_distance(self, comparison_regime, is_weighted_edge_distances, is_normalize_by_tree_size): dmatrix, normalization_factor = self._get_distance_matrix_and_normalization_factor( is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size,) distances = [] for taxon1 in comparison_regime: # subdistances = [dmatrix[taxon1][taxon2] for taxon2 in comparison_regime[taxon1]] # distances.append(min(subdistances)) min_distance = dmatrix[taxon1][comparison_regime[taxon1][0]] for taxon2 in comparison_regime[taxon1][1:]: d = dmatrix[taxon1][taxon2] if d < min_distance: min_distance = d distances.append(min_distance) if distances: return (sum(distances) / normalization_factor) / (len(distances) * 1.0) else: raise error.NullAssemblageException("No taxa in assemblage") def _calculate_standardized_effect_size(self, statisticf_name, comparison_regimes, is_weighted_edge_distances, is_normalize_by_tree_size, null_model_type="taxa.label", num_randomization_replicates=1000, rng=None): result_type = collections.namedtuple("PhylogeneticCommunityStandardizedEffectSizeStatisticCalculationResult", ["obs", "null_model_mean", "null_model_sd", "z", "rank", "p",]) statisticf_kwargs={ "is_weighted_edge_distances": is_weighted_edge_distances, "is_normalize_by_tree_size": is_normalize_by_tree_size } observed_stat_values = {} null_model_stat_values = {} null_model_matrix = self.clone() assert null_model_matrix == self if is_weighted_edge_distances: is_shuffle_phylogenetic_distances = True is_shuffle_phylogenetic_path_steps = False else: is_shuffle_phylogenetic_distances = False is_shuffle_phylogenetic_path_steps = True for rep_idx in range(num_randomization_replicates): null_model_matrix.shuffle_taxa( is_shuffle_phylogenetic_distances=is_shuffle_phylogenetic_distances, is_shuffle_phylogenetic_path_steps=is_shuffle_phylogenetic_distances, is_shuffle_mrca=False, rng=rng) for comparison_regime_idx, comparison_regime in enumerate(comparison_regimes): statisticf_kwargs["comparison_regime"] = comparison_regime if rep_idx == 0: observed_stat_values[comparison_regime_idx] = getattr(self, statisticf_name)(**statisticf_kwargs) null_model_stat_values[comparison_regime_idx] = [] stat_value = getattr(null_model_matrix, statisticf_name)(**statisticf_kwargs) null_model_stat_values[comparison_regime_idx].append(stat_value) results = [] for comparison_regime_idx, comparison_regime in enumerate(comparison_regimes): obs_value = observed_stat_values[comparison_regime_idx] stat_values = null_model_stat_values[comparison_regime_idx] null_model_mean, null_model_var = statistics.mean_and_sample_variance(stat_values) rank = statistics.rank( value_to_be_ranked=obs_value, value_providing_rank=stat_values) if null_model_var > 0: null_model_sd = math.sqrt(null_model_var) z = (obs_value - null_model_mean) / null_model_sd else: null_model_sd = 0.0 z = None p = float(rank) / len(stat_values) result = result_type( obs=obs_value, null_model_mean=null_model_mean, null_model_sd=null_model_sd, z=z, rank=rank, p=p) results.append(result) return results class NodeDistanceMatrix(object): @classmethod def from_tree(cls, tree): ndm = cls() ndm.compile_from_tree(tree=tree) return ndm def __init__(self): self.clear() def clear(self): self._tree_length = None self._num_edges = None self._node_phylogenetic_distances = {} self._node_phylogenetic_path_steps = {} self._mrca = {} def compile_from_tree(self, tree): self.clear() self._tree_length = 0.0 self._num_edges = 0 for node1 in tree.postorder_node_iter(): try: self._tree_length += node1.edge.length except TypeError: # None for edge length pass self._num_edges += 1 if node1 not in self._node_phylogenetic_distances: self._node_phylogenetic_distances[node1] = {node1: 0.0} self._node_phylogenetic_path_steps[node1] = {node1: 0} self._mrca[node1] = {node1: node1} children = node1.child_nodes() for ch_idx, ch1 in enumerate(children): ch1_elen = ch1.edge.length if ch1.edge.length is not None else 0.0 for ch1_subtree_node in list(self._node_phylogenetic_distances[ch1].keys()): if ch1_subtree_node not in self._node_phylogenetic_distances[node1]: d = self._node_phylogenetic_distances[ch1][ch1_subtree_node] + ch1_elen d2 = self._node_phylogenetic_path_steps[ch1][ch1_subtree_node] + 1 self._node_phylogenetic_distances[node1][ch1_subtree_node] = d self._node_phylogenetic_distances[ch1_subtree_node][node1] = d self._node_phylogenetic_path_steps[node1][ch1_subtree_node] = d2 self._node_phylogenetic_path_steps[ch1_subtree_node][node1] = d2 self._node_phylogenetic_distances[node1][ch1] = ch1_elen self._node_phylogenetic_distances[ch1][node1] = ch1_elen self._node_phylogenetic_path_steps[node1][ch1] = 1 self._node_phylogenetic_path_steps[ch1][node1] = 1 for ch2 in children[ch_idx+1:]: self._mrca[ch1][ch2] = node1 self._mrca[ch2][ch1] = node1 ch2_elen = ch2.edge.length if ch2.edge.length is not None else 0.0 d = ch1_elen + ch2_elen self._node_phylogenetic_distances[ch1][ch2] = d self._node_phylogenetic_distances[ch2][ch1] = d self._node_phylogenetic_path_steps[ch1][ch2] = 2 self._node_phylogenetic_path_steps[ch2][ch1] = 2 # Below is ugly, ugly, ugly. Basic idea is to link nodes of each # the subtrees of each of the child nodes of node1. Assumes # that any pairwise comparison of nodes descending from node1 # (as given by nodes in a pairwise comparison with node1) not # already made have their MRCA at node1. for snd1 in self._node_phylogenetic_distances[node1]: for snd2 in self._node_phylogenetic_distances[node1]: if snd1 is snd2: continue if snd1 not in self._node_phylogenetic_distances: self._node_phylogenetic_distances[snd1] = {} self._node_phylogenetic_path_steps[snd1] = {} if snd2 not in self._node_phylogenetic_distances: self._node_phylogenetic_distances[snd2] = {} self._node_phylogenetic_path_steps[snd2] = {} if snd2 not in self._node_phylogenetic_distances[snd1]: self._node_phylogenetic_distances[snd1][snd2] = self._node_phylogenetic_distances[node1][snd1] + self._node_phylogenetic_distances[node1][snd2] self._node_phylogenetic_path_steps[snd1][snd2] = self._node_phylogenetic_path_steps[node1][snd1] + self._node_phylogenetic_path_steps[node1][snd2] if snd1 not in self._node_phylogenetic_distances[snd2]: self._node_phylogenetic_distances[snd2][snd1] = self._node_phylogenetic_distances[node1][snd1] + self._node_phylogenetic_distances[node1][snd2] self._node_phylogenetic_path_steps[snd2][snd1] = self._node_phylogenetic_path_steps[node1][snd1] + self._node_phylogenetic_path_steps[node1][snd2] if snd1 not in self._mrca: self._mrca[snd1] = {} if snd2 not in self._mrca: self._mrca[snd2] = {} if snd2 not in self._mrca[snd1]: self._mrca[snd1][snd2] = node1 self._mrca[snd2][snd1] = node1 def __eq__(self, o): if self.node_namespace is not o.node_namespace: return False return (True and (self._node_phylogenetic_distances == o._node_phylogenetic_distances) and (self._node_phylogenetic_path_steps == o._node_phylogenetic_path_steps) and (self._mrca == o._mrca) and (self._tree_length == o._tree_length) and (self._num_edges == o._num_edges) ) def __iter__(self): for node in self._node_phylogenetic_distances: yield node def __hash__(self): return id(self) def __call__(self, node1, node2): return self.patristic_distance(node1, node2) def __copy__(self): return self.clone() def clone(self): o = self.__class__() o._tree_length = self._tree_length o._num_edges = self._num_edges for src, dest in ( (self._node_phylogenetic_distances, o._node_phylogenetic_distances,), (self._node_phylogenetic_path_steps, o._node_phylogenetic_path_steps,), (self._mrca, o._mrca,), ): for t1 in src: dest[t1] = {} for t2 in src[t1]: dest[t1][t2] = src[t1][t2] return o def mrca(self, node1, node2): """ Returns MRCA of two node objects. """ return self._mrca[node1][node2] def distance(self, node1, node2, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Returns distance between node1 and node2. """ if is_weighted_edge_distances: return self.patristic_distance(node1, node2, is_normalize_by_tree_size=is_normalize_by_tree_size) else: return self.path_edge_count(node1, node2, is_normalize_by_tree_size=is_normalize_by_tree_size) def patristic_distance(self, node1, node2, is_normalize_by_tree_size=False): """ Returns patristic distance between two node objects. """ if node1 is node2: return 0.0 d = self._node_phylogenetic_distances[node1][node2] if is_normalize_by_tree_size: return d / self._tree_length else: return d def path_edge_count(self, node1, node2, is_normalize_by_tree_size=False): """ Returns the number of edges between two node objects. """ if node1 is node2: return 0 d = self._node_phylogenetic_path_steps[node1][node2] if is_normalize_by_tree_size: return float(d) / self._num_edges else: return d def distances(self, is_weighted_edge_distances=True, is_normalize_by_tree_size=False): """ Returns list of patristic distances. """ dmatrix, normalization_factor = self._get_distance_matrix_and_normalization_factor( is_weighted_edge_distances=is_weighted_edge_distances, is_normalize_by_tree_size=is_normalize_by_tree_size, ) results = [] nodes = list(dmatrix.keys()) for node_idx1, node1 in enumerate(nodes[:-1]): for node_idx2, node2 in enumerate(nodes[node_idx1+1:]): results.append(dmatrix[node1][node2]/normalization_factor) return results def _get_distance_matrix_and_normalization_factor(self, is_weighted_edge_distances, is_normalize_by_tree_size): if is_weighted_edge_distances: dmatrix = self._node_phylogenetic_distances if is_normalize_by_tree_size: normalization_factor = self._tree_length else: normalization_factor = 1.0 else: dmatrix = self._node_phylogenetic_path_steps if is_normalize_by_tree_size: normalization_factor = float(self._num_edges) else: normalization_factor = 1.0 return dmatrix, normalization_factor