Manual#

physquirrel reconstructs level-1 phylogenetic networks from quartet data using the Squirrel algorithm. This guide covers the full workflow: preparing input data, running the algorithm, using convenience shortcuts, and working with the output network. Note that physquirrel is built on top of PhyloZoo; see the PhyloZoo documentation for details on working with the core datatypes (SemiDirectedPhyNetwork, Split, Quartet, DistanceMatrix, and MSA). The main datatypes from phylozoo used in physquirrel are re-exported in the top-level namespace for convenience.

The standard import convention used throughout is:

import physquirrel as psq

Input data#

physquirrel works with three input types, depending on how far upstream you want to start:

Input

Type

Entry point

Quarnets (4-leaf subnetworks)

SqQuartetProfileSet

squirrel()

Pairwise distance matrix

DistanceMatrix

squirrel_from_distances()

Multiple-sequence alignment

MSA

squirrel_from_msa()

From a set of quarnets#

A quarnet is a 4-leaf subnetwork that can be either a split (tree-like) or a 4-cycle (network-like). In physquirrel, quarnets are represented by SqQuartetProfile objects, which are then collected into a SqQuartetProfileSet. These objects build on top of phylozoo’s QuartetProfile: the standard object to represent quartets.

For full control, construct an SqQuartetProfileSet from individual SqQuartetProfile objects:

import physquirrel as psq

# A split quarnet: AB|CD
q = psq.Quartet(psq.Split({'A', 'B'}, {'C', 'D'}))
profile = psq.SqQuartetProfile([q])

# A cycle quarnet: two resolved quartets sharing four taxa
q1 = psq.Quartet(psq.Split({'A', 'B'}, {'C', 'D'}))
q2 = psq.Quartet(psq.Split({'A', 'D'}, {'B', 'C'}))
cycle_profile = psq.SqQuartetProfile([q1, q2], reticulation_leaf='A')

profileset = psq.SqQuartetProfileSet(profiles=[profile, cycle_profile])

Each SqQuartetProfile holds exactly 1 or 2 resolved quartets on the same four taxa. A 2-quartet profile encodes a 4-cycle; an optional reticulation_leaf identifies which taxon is below the hybrid node. Quarnets can also be weighted by passing a profile_weight to the constructor; this is used in cycle resolution and scoring to downweight low-confidence profiles.

From a distance matrix#

The delta heuristic infers a profile set from pairwise distances, represented by a DistanceMatrix. For each 4-taxon subset it computes a score \(\delta \in [0,1]\); below the threshold lam the quartet is classified as a split, otherwise as a cycle.

import numpy as np
import physquirrel as psq

matrix = np.array([
    [0.0, 0.1, 0.9, 0.9],
    [0.1, 0.0, 0.9, 0.9],
    [0.9, 0.9, 0.0, 0.1],
    [0.9, 0.9, 0.1, 0.0],
])
dm = psq.DistanceMatrix(matrix, labels=['A', 'B', 'C', 'D'])

profileset = psq.delta_heuristic(dm, lam=0.3, weight=True)

The weight parameter controls whether profile confidence (distance to the threshold) is stored as profile-level weights in the returned SqQuartetProfileSet.

From an MSA#

When starting from aligned sequences, these first need to be converted into a distance matrix. This can be done through phylozoo’s hamming_distances() function, which computes pairwise Hamming distances between sequences in the MSA. The resulting distance matrix can then be passed to the delta heuristic as described above. Alternatively, physquirrel provides a convenience function that handles this internally: it computes Hamming distances and applies the delta heuristic in one call.

Use the lam parameter to control the split/cycle classification threshold, and the optional weight parameter to store profile confidence as weights in the returned SqQuartetProfileSet.

MSAs are represented by MSA. Both FASTA and NEXUS formats are supported. Load from a file with load(), or construct directly from a dictionary of sequences:

import physquirrel as psq

# From a FASTA file (.fasta, .fa, .fas)
msa = psq.MSA.load("alignment.fasta")

# From a NEXUS file (.nexus, .nex, .nxs) — specify format explicitly
msa = psq.MSA.load("alignment.nex", format="nexus")

# From a dictionary of sequences
msa = psq.MSA({"A": "ACGT", "B": "ACGT", "C": "TGCA", "D": "TGCA"})

profileset = psq.delta_heuristic_from_msa(msa, lam=0.3)

Running Squirrel#

The main entry point is the function squirrel(). It accepts an SqQuartetProfileSet (see previous section)and returns a SemiDirectedPhyNetwork:

import physquirrel as psq

network = psq.squirrel(profileset)

print(f"eNewick string:  {network.to_string()}")

The pipeline runs internally:

  1. T* tree (tstar_tree()) — B*-compatible split system → initial tree.

  2. Adapted quartet joining (adapted_quartet_joining()) — refines the tree using quartet support.

  3. Unresolve tree (unresolve_tree()) — iteratively contracts low-support splits, yielding a sequence of contracted trees from fully resolved to star.

  4. Cycle resolution (resolve_cycles()) — inserts reticulation cycles at internal vertices of each contracted tree using a TSP-based circular ordering.

  5. Scoring (sqprofileset_similarity()) — each candidate network is scored by profile similarity (consistency with the input profile set); the best is returned.

Rooting at an outgroup#

Pass outgroup to return a DirectedPhyNetwork rooted at the edge incident to the specified leaf:

import physquirrel as psq

rooted = psq.squirrel(profileset, outgroup='A')
assert isinstance(rooted, psq.DirectedPhyNetwork)
assert rooted.root_node is not None

Convenience shortcuts#

For common workflows, physquirrel provides one-call convenience functions that compose delta_heuristic() and squirrel() internally.

From a distance matrix:

import physquirrel as psq

network = psq.squirrel_from_distances(dm, lam=0.3, outgroup='A')

From an MSA:

import physquirrel as psq

network = psq.squirrel_from_msa(msa, lam=0.3)

All keyword arguments accepted by squirrel() (e.g. rho, tsp_threshold) can be forwarded through both convenience functions.

Parallelization#

The squirrel() algorithm can be parallelized by handling the candidate trees in parallel during cycle resolution and scoring. This is controlled by the optional parallel argument, which accepts a ParallelConfig object from phylozoo. By default, the algorithm runs sequentially.

Threading (shared memory, lower overhead, good for 2–4 cores):

from phylozoo.utils.parallel import ParallelConfig, ParallelBackend
import physquirrel as psq

# Use threading with 4 worker threads
parallel_config = ParallelConfig(backend=ParallelBackend.THREADING, n_jobs=4)
network = psq.squirrel(profileset, parallel=parallel_config)

Multiprocessing (true parallelism, higher overhead, good for 4+ cores):

from phylozoo.utils.parallel import ParallelConfig, ParallelBackend
import physquirrel as psq

# Use multiprocessing with 8 worker processes
parallel_config = ParallelConfig(backend=ParallelBackend.MULTIPROCESSING, n_jobs=8)
network = psq.squirrel(profileset, parallel=parallel_config)

Sequential (default, no parallelization):

import physquirrel as psq

# Omit the parallel argument or pass None
network = psq.squirrel(profileset)
# or explicitly:
network = psq.squirrel(profileset, parallel=None)

Notes:

  • Threading is recommended for most cases due to lower overhead; it works well for 2–4 cores.

  • Multiprocessing enables true parallelism but incurs higher overhead (serialization, process spawning); use it for large datasets or 4+ cores.

  • The number of jobs (n_jobs) should typically match your CPU core count. Use n_jobs=-1 to automatically detect the number of available cores.

I/O and plotting#

Loading and saving quartet profile sets#

SqQuartetProfileSet supports two file formats. The class uses the IOMixin interface from phylozoo, which infers the format from the file extension by default. Alternatively, the format can be specified explicitly when loading or saving. The four main PhyloZoo functions for I/O are: save, to_string, load, and from_string.

psq format (.psq, default)#

The psq format is a compact JSON file that stores all information losslessly: the full taxon set, per-profile quartet weights, reticulation leaf assignments, and profile-level confidence weights. It is the recommended format for archiving and exchanging profile sets between physquirrel runs.

import physquirrel as psq

# Infer a profile set from a distance matrix
profileset = psq.delta_heuristic(dm, lam=0.3)

# Save to file — format inferred from .psq extension
profileset.save("profiles.psq")

# Load back
loaded = psq.SqQuartetProfileSet.load("profiles.psq")

A profiles.psq file looks like:

{
  "taxa": ["A", "B", "C", "D"],
  "profiles": [
    {
      "taxa": ["A", "B", "C", "D"],
      "quartets": [
        {"quartet": {"type": "resolved",
                     "split": {"set1": ["A", "B"], "set2": ["C", "D"]}},
         "weight": 1.0}
      ],
      "reticulation_leaf": null,
      "profile_weight": 0.85
    }
  ]
}

profile_list format (.txt)#

The profile_list format is a plain-text file with one profile per line. Each line begins with a type tag followed by four taxon names and an optional weight:

# physquirrel profile list
# Lines starting with '#' are comments and are ignored.
SQ A B C D 0.85    # split quarnet: split AB|CD, confidence weight 0.85
4C A B C E 0.60    # 4-cycle: circular order A,B,C,E; A is the reticulation leaf

  • SQ a b c d [weight] — a split quarnet with split \(ab|cd\).

  • 4C a b c d [weight] — a 4-cycle with circular order \(a,b,c,d\) where \(a\) is the reticulation leaf (below the hybrid node).

The weight field is optional and defaults to 1.0. This format is compatible with the quarnet file format used in earlier versions of physquirrel.

import physquirrel as psq

# Save as profile_list
profileset.save("profiles.txt", format="profile_list")

# Or get the string directly
txt = profileset.to_string(format="profile_list")
print(txt)

# Load back — format must be specified explicitly for .txt
loaded = psq.SqQuartetProfileSet.load("profiles.txt", format="profile_list")
loaded2 = psq.from_profile_list(txt)

Saving networks#

SemiDirectedPhyNetwork and DirectedPhyNetwork are provided by phylozoo and support the eNewick format (default) via the standard IOMixin interface:

import physquirrel as psq

# Serialize to an eNewick string
enewick_str = network.to_string()              # eNewick by default

# Write directly to file (.nwk, .newick, .enewick, .enw)
network.save("network.nwk")

See the PhyloZoo documentation for the full list of supported formats and options.

Visualization#

physquirrel networks can be plotted directly using phylozoo’s unified plot() function, which dispatches by network type and uses matplotlib:

from phylozoo.viz import plot

plot(network)                          # semi-directed network, twopi layout
plot(network, layout="pz-radial")      # radial layout (trees only)
plot(rooted_network)                   # :class:`~phylozoo.core.network.dnetwork.base.DirectedPhyNetwork`, DAG layout

The viz extra dependency must be installed (pip install physquirrel[viz] or pip install physquirrel[graphviz]); see the Installation Guide for details.