A multiflow global optimization approach for plasmid binning in short read assembly graph: PlasBin-HMF

Install

[!IMPORTANT] PlasBin-HMF relies on the commercial Gurobi MILP solver. If you are an academic, you can get a free license, see https://www.gurobi.com/academics#licenses.

In a virtual environment:

pip install plasbin-hmf

Or if you want to install it via uv:

# As a dependency of your project:
uv add plasbin-hmf
# As an uv tool:
uv tool install plasbin-hmf

Usage

In your virtual environment:

plasbin-hmf --help

You can also use uv tools:

uvx plasbin-hmf --help

To get the plasmid bins:

uvx plasbin-hmf run --help

Options

Algorithm parameters can be given via CLI or you can use --config config.yaml:

bins:
    circular:
        max_with_seeds: <int> | auto    # default: auto
        max_optional_seeds: <int>       # default: 2
    partially_circular:
        max_with_seeds: <int> | auto    # default: auto
        max_optional_seeds: <int>       # default: 0
constraints:
    min_bin_size: <int>                             # default: 1000
    min_coverage: <(0, inf)>                        # default: 0.3
    plasmidness_coefficient: <[0, 1]>               # default: 0
    max_coverage_ratio: <(1, inf]> | null           # default: 20
    partially_circular_s_t_constraints: <bool>      # default: false
input:
    input_format: PlasBin-HMF | PlasBin-flow    # default: PlasBin-HMF
    gfa_origin: Unicycler  | SKESA              # default: Unicycler
    default_plasmidness: <[-1, 1]>              # default: 0
    contig_length_threshold: <int>              # default: 100
output:
    keep_optimization_traces: <bool>    # default: false
    plasbin_flow_output: <bool>         # default: false

Gurobi can also be configured, via CLI of via --gurobi-config gurobi_cfg.yaml:

mip_gap: <float>            # default: 0.0
time_limit: <float> | null  # default: null
threads: <int> | all        # default: all

Input

PlasBin-HMF takes as input:

• An assembly graph (GFA v1 file produced by Unicycler or SKESA, see --gfa-origin option), can be compressed.
• A file listing the seed contigs
• A file associating to each contig a plasmidness score

Seed contigs (e.g. seeds.tsv, tabulation-separated values):

Column ID	Type	Description
contig_id	str	Contig ID (corresponds to the ID of the segment line in the GFA)

Plasmidness scores (e.g. plasmidness.tsv, tabulation-separated values):

Column ID	Type	Description
contig_id	str	Contig ID (corresponds to the ID of the segment line in the GFA)
plasmidness_score	float	Plasmidness score

[!TIP] If you use PlasBin-flow (another plasmid binning tool), you can use the same seed and plasmidness files via the option --input-format PlasBin-flow. Otherwise (--input-format PlasBin-HMF), PlasBin-HMF expects the following files.

Output

Suppose --output-dir/-o output_dir is provided.

PlasBin-HMF keeps your configuration for future reproducibility:

📂 output_dir
├── ℹ️ configs.yaml
└── ℹ️ gurobi_cfg.yaml

Solution

📂 output_dir
├── ℹ️ solution_metadata.yaml
├── 📄 bins.tsv
└── 📂 subgraphs
    ├── 📄 link_arcs.tsv
    ├── 📄 source_arcs.tsv
    └── 📄 sink_arcs.tsv

[!TIP] You can also get the extracted bins under the PlasBin-flow format with the option --plasbin-flow-output. In that case, 📂 output_dir will also contain 📄 plasbin_flow_bins.tsv file in the PlasBin-flow format. The Flow column corresponds to the normalization factor (which is the minimum positive arc flow, and not the total flow which can be 0 for circular bins)

Predicted bins

The predicted bins are stored in the 📄 bins.tsv file (tabulation separated columns):

Column ID	Type	Description
bin_id	str	Bin ID
contig_id	str	Contig ID in the bin (corresponds to the ID of the segment line in the GFA)
contig_norm_cov	int	Normalized coverage of the contig

[!NOTE] A contig can belong to several bins. In that case there are multiple rows with the same contig_id value but different bin_id.

Solution networks

PlasBin-HMF uses the assembly graph topology to produce the bins. Folder 📂 subgraphs contains the arcs defining the subnetworks associated with the bins.

For each subnetwork (bin), the set of arcs is the union of:

• a set of link-arcs (arcs between two oriented contigs);
• a set of source-arcs (arcs from the source to the oriented contigs);
• a set of sink-arcs (arcs from the oriented contigs to the sink).

File 📄 link_arcs.tsv contains the solution link arcs for all the subnetworks. A link is an arc that corresponds to a link in the assembly graph (between two oriented contigs).

Column ID	Type	Description
bin_id	str	Bin ID
pred_id	str	Predecessor oriented contig ID
pred_or	{'+', '-'}	Predecessor oriented contig orientation
succ_id	str	Successor oriented contig ID
succ_or	{'+', '-'}	Successor oriented contig orientation
flow	float	Flow on the arc (always strictly positive)

File 📄 source_arcs.tsv contains the solution source arcs for all the subnetworks.

Column ID	Type	Description
bin_id	str	Bin ID
succ_id	str	Successor oriented contig ID
succ_or	{'+', '-'}	Successor oriented contig orientation
flow	float	Flow on the arc (can be 0)

File 📄 sink_arcs.tsv contains the solution sink arcs for all the subnetworks

Column ID	Type	Description
bin_id	str	Bin ID
pred_id	str	Predecessor oriented contig ID
pred_or	{'+', '-'}	Predecessor oriented contig orientation
flow	float	Flow on the arc (can be 0)

Solution metadata

File ℹ️ solution_metadata.yaml is a YAML stream. Each document (separated by --- at the top) describes one predicted bin:

---
bin_id: <str> # e.g. Bin_1
property:
    topology: circular|partially_circular
    seed_constraint: required|optional
stats:
    active_contigs: <int>               # Number of active (unoriented) contig
    active_contig_vertices: <int>       # Number of active oriented contig
    active_source_arcs: <int>           # Number of active source arcs
    active_sink_arcs: <int>             # Number of active link arcs
    active_link_arcs: <int>             # Number of active sink arcs
    cumulative_contig_length: <int>     # Do not account for the multiplicity of a contig
    plasmidness_score: <float>          # A good bin has a plasmidness score greater than 0
    minimum_positive_arc_flow: <float>  # This can be used as a coverage normalization factor
    maximum_contig_coverage: <float>    # Maximum contig coverage in the bin
origin:                 # only if --keep-optimization-traces
    instance_id: <str>  # e.g. instance_1
    multiflow_id: <str> # e.g. multiflow_2
    flow_id: <str>      # e.g. flow_1

If there is no solution

In case PlasBin-HMF was not able to find any bin under the parametrized constraints, the program produces:

📂 output_dir
└── ℹ️ no_solution.yaml

Where no_solution.yaml contains:

status: no_solution

Optimization traces

If you provide --keep-optimization-traces, the optimization traces will be kept in the output directory.

📂 output_dir
└── 📂 optimization_traces
    ├── ℹ️ traces_metadata.yaml
    ├── 🗄️ instance_1.tar.xz # or .zip
    ├── 🗄️ instance_2.tar.xz # or .zip
    └── ...

[!NOTE] The archive format depends on the user choice via --compression-format or on the first supported compression backend available on the system (by default auto).

Compressed instance results

Directory 📂 optimization_traces keeps trace for each instance of the steps to find the bins corresponding to that instance. An instance is defined by:

• an ID (prefixed with instance_ and the same ID given in the instance_id key in solution_metadata.yaml);
• a network such that removing the source and the sink vertices still gives a connected subgraph;
• a topology (circular or partially circular);
• a seed constraint (required or optional);

[!TIP] The instance IDs reflect the tree of instances: instance_2 is associated with the second connected component, instance_2_0 is associated with the first component after solving the parent instance instance_2.

Full network
├── instance_1
│   ├── instance_1_1
│   │   └── ...
│   └── ...
└── ...

One compressed file corresponds to one instance. Decompressing it gives one directory (for example instance_1.<archive>):

📂 instance_1  # Instance ID prefixed with `instance_`
├── ℹ️ metadata.yaml
├── 📂 multiflow_1  # The 1-multiflow attempt was succesful
│   ├── 📃 gurobi.log
│   ├── 📄 gurobi.sol
│   └── ℹ️ solution_metadata.yaml
└── 📂 multiflow_2  # The 2-multiflow attempt failed
    ├── 📃 gurobi.log
    └── ℹ️ infeasible.yaml

File ℹ️ metadata.yaml describes the instance:

instance_id: <str>
multiflow_size_upper_bound: <int>
required_bin_property:
    topology: circular|partially_circular
    seed_constraint: required|optional
input_stats:
    number_of_contigs: <int>
    number_of_links: <int>
    number_of_seeds: <int>
    null_plasmidness: <bool>  # True if every contig has null plasmidness. The multiflow score function depends on this boolean.
parent_instance: <str> | null # Instance ID, if any. Otherwise null

For each instance, there is at least one optimization attempt. An attempt multiflow_n means an attempt to find the best (n)-multiflow with the instance properties.

In a multiflow attempt (multiflow_1 means a multiflow of size 1), file ℹ️ solution_metadata.yaml contains:

multiflow_id: <str>
multiflow_size: <int>
status: <str>  # Gurobi status, see https://docs.gurobi.com/projects/optimizer/en/current/reference/numericcodes/statuscodes.html
multiflow_score: <float>
flows:
    flow_ID:  # Flow ID e.g. flow_1
        active_contigs: <int>               # Number of active (unoriented) contig
        active_contig_vertices: <int>       # Number of active oriented contig
        active_source_arcs: <int>           # Number of active source arcs
        active_sink_arcs: <int>             # Number of active link arcs
        active_link_arcs: <int>             # Number of active sink arcs
        total_flow: <float>                 # Should be 0 in the circular context
        plasmidness_score: <float>          # A good bin has a plasmidness score greater than 0
        cumulative_contig_length: <int>     # Do not account for the multiplicity of a contig
        min_positive_arc_flow: <float>      # Corresponds to the minimum positive arc flow
        max_positive_arc_flow: <float>      # Corresponds to the maximum positive arc flow
        min_positive_contig_inflow: <float> # Corresponds to the minimum positive contig inflow (forward + reverse orientations)
        max_positive_contig_inflow: <float> # Corresponds to the maximum positive contig inflow (forward + reverse orientations)

In case PlasBin-HMF was not able to find any bin under the parametrized constraints, the program produces file ℹ️ infeasible.yaml:

multiflow_id: <str>
multiflow_size: <int>
status: <str>  # Gurobi status, see https://docs.gurobi.com/projects/optimizer/en/current/reference/numericcodes/statuscodes.html

Traces metadata

File ℹ️ traces_metadata.yaml is a YAML stream of all the instance_id/metadata.yaml, see instance metadata.

Each document (separated by --- at the top) describes one optimization instance.

---
instance_id: <str>
multiflow_size_upper_bound: <int>
required_bin_property:
    topology: circular|partially_circular
    seed_constraint: required|optional
input_stats:
    number_of_contigs: <int>
    number_of_links: <int>
    number_of_seeds: <int>
parent_instance: <str> | null # Instance ID, if any. Otherwise null
# ---
# ... Another instance