Neighbor-Joining Anchors

Anchor trees are neighbor-joining trees built from the gene alignments, with corrected pairwise distances (JC69 for DNA, Kimura 1983 for protein).

nj_tree

nj_tree(alignment: 'Alignment', taxon_table: TaxonTable) -> HifukuTree

Compute a neighbor-joining anchor tree from a compressed alignment.

The NJ algorithm runs on all of the alignment's local taxa. The returned HifukuTree is restricted to the taxa currently in taxon_table, so the caller controls the effective namespace by pre-populating taxon_table.

Negative branch lengths produced by the NJ algorithm (a known artifact on non-tree-metric distance matrices) are clamped to zero.

Parameters:

Name	Type	Description	Default
alignment	'Alignment'	Compressed alignment (DNA or AA) produced by load_alignment.	required
taxon_table	TaxonTable	Shared global TaxonTable. Leaf taxa from the NJ tree that are in this table are kept; the rest are pruned. Any kept taxa not already in the table are registered via require.	required

Returns:

Type	Description
HifukuTree	NJ tree restricted to taxon_table, with global_leaf_indices set.

Notes

Distance model (per seq_type):

• 'DNA': Jukes-Cantor 1969 correction, d = -(3/4) ln(1-(4/3)p). Jukes & Cantor (1969).
• 'AA': Kimura 1983 correction, d = -ln(1-p-0.2p^2). Kimura (1983).

Saturated pairs are capped at a large finite distance; see _SATURATED_DIST.

Neighbor-joining: saitou1987nj.

Source code in src/hifuku/nj.py

def nj_tree(alignment: "Alignment", taxon_table: TaxonTable) -> HifukuTree:
    """
    Compute a neighbor-joining anchor tree from a compressed alignment.

    The NJ algorithm runs on all of the alignment's local taxa. The returned
    HifukuTree is restricted to the taxa currently in *taxon_table*, so the
    caller controls the effective namespace by pre-populating *taxon_table*.

    Negative branch lengths produced by the NJ algorithm (a known artifact
    on non-tree-metric distance matrices) are clamped to zero.

    Parameters
    ----------
    alignment :
        Compressed alignment (DNA or AA) produced by ``load_alignment``.
    taxon_table :
        Shared global TaxonTable.  Leaf taxa from the NJ tree that are in
        this table are kept; the rest are pruned.  Any kept taxa not already
        in the table are registered via ``require``.

    Returns
    -------
    HifukuTree
        NJ tree restricted to *taxon_table*, with ``global_leaf_indices`` set.

    Notes
    -----
    Distance model (per seq_type):

    - ``'DNA'``: Jukes-Cantor 1969 correction, d = -(3/4) ln(1-(4/3)p).
      Jukes & Cantor (1969).
    - ``'AA'``: Kimura 1983 correction, d = -ln(1-p-0.2p^2).
      Kimura (1983).

    Saturated pairs are capped at a large finite distance; see _SATURATED_DIST.

    Neighbor-joining: saitou1987nj.
    """
    labels = alignment.index.labels
    if len(labels) < 3:
        raise ValueError(
            f"Alignment has only {len(labels)} taxa; NJ requires at least 3."
        )

    # Build the corrected pairwise distance matrix.
    p_dist = _p_distances(alignment)
    if alignment.seq_type == "DNA":
        d_corr = _jc69_correct(p_dist)
    else:
        d_corr = _kimura83_correct(p_dist)

    # Run neighbor joining and parse the newick into a temporary local
    # TaxonTable, so that pruned taxa never reach the global namespace.
    newick = _neighbor_join(list(labels), d_corr)
    local_tt = TaxonTable()
    tree = read_tree(newick, local_tt, from_string=True)

    # Clamp negative branch lengths (NJ artifact on non-tree-metric input).
    neg_mask = tree.edge_len < 0.0
    if np.any(neg_mask):
        tree.edge_len[neg_mask] = 0.0

    # Restrict to the labels currently in the caller's taxon_table (the
    # shared global namespace).  Using a local_tt above avoids accidentally
    # registering pruned taxa into the global namespace.
    ns_set = set(taxon_table.labels)
    keep = [lbl for lbl in tree.leaf_labels if lbl in ns_set]
    if len(keep) < 3:
        raise ValueError(
            f"Fewer than 3 NJ tree taxa are in the global namespace "
            f"({len(keep)} kept from {tree.n_leaves} in the tree)."
        )
    result = tree if len(keep) == tree.n_leaves else tree.restrict(keep)

    # Register kept taxa into the global namespace and set global_leaf_indices.
    for lbl in result.leaf_labels:
        taxon_table.require(lbl)
    glidx = np.full(result.n_nodes, -1, dtype=np.int32)
    for i in range(result.n_leaves):
        glidx[i] = taxon_table.index_of(result.names[i])
    result.global_leaf_indices = glidx

    return result