Enumerating Fragment Combinations

Problem

You want to enumerate all the adjacent fragment combinations returned by the fragmentation methods of OEMedChem TK.

../_images/enumfrags.png

Figure 1: Example molecule with fragment highlighting

For example, in Figure 1 you want to return fragment combinations: A+B, B+C, …, A+B+C, …, A+B+C+D, … etc. However, you do not want disconnected fragment combinations such as D.A+B. See examples in Table 1.

Table 1. Example of fragment combinations

valid - all fragments adjacent

invalid - disconnected fragments

../_images/enumfrags-valid.png ../_images/enumfrags-invalid.png

Ingredients

Difficulty Level

🌶️ 🌶️

Download

Download code

enumfrags.py

See also Usage subsection.

Source Code

enumfrags
#!/usr/bin/env python3
# (C) 2026 Cadence Design Systems, Inc. (Cadence)
# All rights reserved.
# TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
# provided to current licensees or subscribers of Cadence products or
# SaaS offerings (each a "Customer").
# Customer is hereby permitted to use, copy, and modify the Sample Code,
# subject to these terms. Cadence claims no rights to Customer's
# modifications. Modification of Sample Code is at Customer's sole and
# exclusive risk. Sample Code may require Customer to have a then
# current license or subscription to the applicable Cadence offering.
# THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
# NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
# PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
# liable for any damages or liability in connection with the Sample Code
# or its use.

"""Enumerate connected fragment combinations."""

import argparse
import enum
import os
import sys
from collections.abc import Callable, Iterator
from itertools import combinations
from pathlib import Path

from openeye import oechem, oemedchem
from rich.console import Console
from rich_argparse import HelpPreviewAction, RichHelpFormatter

__SCRIPT_NAME__ = Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Enumerate connected fragment combinations."
__SCRIPT_TOOLKITS__ = ["oechem", "oemedchem"]
__SCRIPT_CATEGORIES__ = ["cheminformatics"]


def parse_options() -> argparse.Namespace:
    """Set up command line options."""
    parser = argparse.ArgumentParser(
        add_help=True,
        formatter_class=RichHelpFormatter,
        description="[yellow]" + __SCRIPT_DESC__ + "[/yellow]",
    )

    io_group = parser.add_argument_group("Input/output options")
    io_group.add_argument(
        "--in",
        dest="input",
        type=str,
        required=True,
        metavar="MOL-FILE",
        help="input molecule file",
    )
    io_group.add_argument(
        "--out",
        dest="output",
        type=str,
        required=True,
        metavar="MOL-FILE",
        help="output molecule file",
    )

    frag_group = parser.add_argument_group("Fragmentation options")
    frag_group.add_argument(
        "--frag-type",
        type=FragmentationType,
        default=FragmentationType.FunctionalGroup,
        choices=list(FragmentationType),
        help="fragmentation type (default: %(default)s)",
    )

    parser.add_argument("--help-image", action=HelpPreviewAction)
    parser.add_argument(
        "--save-console-svg",
        default=False,
        action="store_true",
        help=f"run command and capture console output in {__SCRIPT_NAME__}.svg file",
    )

    return parser.parse_args()


def main() -> int:
    """Enumerate connected fragment combinations."""
    args = parse_options()
    console = Console(record=args.save_console_svg)

    ifs = oechem.oemolistream()
    if not ifs.open(args.input):
        oechem.OEThrow.Fatal(f"Cannot open input file: {args.input}")

    ofs = oechem.oemolostream()
    if not ofs.open(args.output):
        oechem.OEThrow.Fatal(f"Cannot open output file: {args.output}")

    mol = oechem.OEGraphMol()
    if not oechem.OEReadMolecule(ifs, mol):
        oechem.OEThrow.Fatal(f"Cannot read molecule from: {args.input}")
    console.print(
        f"Input molecule: {oechem.OEMolToSmiles(mol)}", highlight=False, markup=False
    )

    frag_func = _get_fragmentation_function(args.frag_type)

    frags = list(frag_func(mol))
    console.print(f"{len(frags)} fragments generated")

    frag_combs = get_fragment_combinations(mol, frags)
    console.print(f"{len(frag_combs)} fragment combinations generated")

    for frag in frag_combs:
        oechem.OEWriteMolecule(ofs, frag)

    if args.save_console_svg:
        console.save_svg(f"{__SCRIPT_NAME__}.svg", title="output")

    return os.EX_OK


def is_adjacent_atom_bond_sets(
    frag_a: oechem.OEAtomBondSet, frag_b: oechem.OEAtomBondSet
) -> bool:
    """Check if two atom/bond sets share an adjacent bond."""
    for atom_a in frag_a.GetAtoms():
        for atom_b in frag_b.GetAtoms():
            if atom_a.GetBond(atom_b) is not None:
                return True
    return False


def is_adjacent_atom_bond_set_combination(
    frag_list: list[oechem.OEAtomBondSet],
) -> bool:
    """Check if a list of fragments forms a single connected component."""
    parts = [0] * len(frag_list)
    num_parts = 0

    for idx, frag in enumerate(frag_list):
        if parts[idx] != 0:
            continue

        num_parts += 1
        parts[idx] = num_parts
        traverse_fragments(frag, frag_list, parts, num_parts)

    return num_parts == 1


def traverse_fragments(
    act_frag: oechem.OEAtomBondSet,
    frag_list: list[oechem.OEAtomBondSet],
    parts: list[int],
    num_parts: int,
) -> None:
    """Recursively traverse adjacent fragments to assign connected components."""
    for idx, frag in enumerate(frag_list):
        if parts[idx] != 0:
            continue

        if not is_adjacent_atom_bond_sets(act_frag, frag):
            continue

        parts[idx] = num_parts
        traverse_fragments(frag, frag_list, parts, num_parts)


def combine_and_connect_atom_bond_sets(
    frag_list: tuple[oechem.OEAtomBondSet, ...],
) -> oechem.OEAtomBondSet:
    """Combine fragment atom/bond sets and add connecting bonds."""
    combined = oechem.OEAtomBondSet()
    for frag in frag_list:
        for atom in frag.GetAtoms():
            combined.AddAtom(atom)
        for bond in frag.GetBonds():
            combined.AddBond(bond)

    for atom_a in combined.GetAtoms():
        for atom_b in combined.GetAtoms():
            if atom_a.GetIdx() < atom_b.GetIdx():
                continue

            bond = atom_a.GetBond(atom_b)
            if bond is None:
                continue
            if combined.HasBond(bond):
                continue

            combined.AddBond(bond)

    return combined


def get_fragment_atom_bond_set_combinations(
    frag_list: list[oechem.OEAtomBondSet],
) -> list[oechem.OEAtomBondSet]:
    """Generate all valid adjacent fragment combinations."""
    frag_combs: list[oechem.OEAtomBondSet] = []

    num_frags = len(frag_list)
    for n in range(2, num_frags):
        for frag_comb in combinations(frag_list, n):
            if is_adjacent_atom_bond_set_combination(list(frag_comb)):
                frag = combine_and_connect_atom_bond_sets(frag_comb)
                frag_combs.append(frag)

    return frag_combs


def get_fragment_combinations(
    mol: oechem.OEMolBase, frag_list: list[oechem.OEAtomBondSet]
) -> list[oechem.OEGraphMol]:
    """Generate all connected fragment combination molecules from a fragmented molecule."""
    fragments: list[oechem.OEGraphMol] = []
    frag_combs = get_fragment_atom_bond_set_combinations(frag_list)

    for f in frag_combs:
        frag_atom_pred = oechem.OEIsAtomMember(f.GetAtoms())
        frag_bond_pred = oechem.OEIsBondMember(f.GetBonds())

        fragment = oechem.OEGraphMol()
        adjust_h_count = True
        oechem.OESubsetMol(
            fragment, mol, frag_atom_pred, frag_bond_pred, adjust_h_count
        )
        fragments.append(fragment)

    return fragments


class FragmentationType(enum.Enum):
    """Molecule fragmentation type."""

    FunctionalGroup = "func-group"
    RingChain = "ring-chain"
    RingLinkerSideChain = "ring-linker-sidechain"

    def __str__(self) -> str:
        """Convert to string representation."""
        return self.value


def _get_fragmentation_function(
    frag_type: FragmentationType,
) -> Callable[[oechem.OEMolBase], Iterator[oechem.OEAtomBondSet]]:
    """Return the fragmentation function for the given type."""
    match frag_type:
        case FragmentationType.RingChain:
            return oemedchem.OEGetRingChainFragments
        case FragmentationType.RingLinkerSideChain:
            return oemedchem.OEGetRingLinkerSideChainFragments
    return oemedchem.OEGetFuncGroupFragments


setattr(main, "__SCRIPT_NAME__", __SCRIPT_NAME__)
setattr(main, "__SCRIPT_DESC__", __SCRIPT_DESC__)
setattr(main, "__SCRIPT_TOOLKITS__", __SCRIPT_TOOLKITS__)
setattr(main, "__SCRIPT_CATEGORIES__", __SCRIPT_CATEGORIES__)

if __name__ == "__main__":
    sys.exit(main())

Solution

The OEMedChem TK currently provides three ways to partition a molecule into fragments (see examples in Table 2):

Table 2. Example of the fragmentation methods of OEMedChem TK

OEGetFuncGroupFragments

OEGetRingChainFragments

OEGetRingLinkerSideChainFragments

../_images/frags2img-01.svg ../_images/frags2img-02.svg ../_images/frags2img-03.svg

The get_fragment_combinations function takes a molecule and a list of OEAtomBondSet objects. This list is generated by using one of the fragmentation methods shown in Table 2. Each OEAtomBondSet object stores the atoms and bonds of a fragment generated for a given molecule. These parameters are passed to the get_fragment_atom_bond_set_combinations function, which returns the list of adjacent fragment combinations. Each fragment combination is then passed to the OESubsetMol function using an OEIsAtomMember and an OEIsBondMember predicate. The subset molecule generated by OESubsetMol is then appended to the returned molecule list.

def get_fragment_combinations(
    mol: oechem.OEMolBase, frag_list: list[oechem.OEAtomBondSet]
) -> list[oechem.OEGraphMol]:
    """Generate all connected fragment combination molecules from a fragmented molecule."""
    fragments: list[oechem.OEGraphMol] = []
    frag_combs = get_fragment_atom_bond_set_combinations(frag_list)

    for f in frag_combs:
        frag_atom_pred = oechem.OEIsAtomMember(f.GetAtoms())
        frag_bond_pred = oechem.OEIsBondMember(f.GetBonds())

        fragment = oechem.OEGraphMol()
        adjust_h_count = True
        oechem.OESubsetMol(
            fragment, mol, frag_atom_pred, frag_bond_pred, adjust_h_count
        )
        fragments.append(fragment)

    return fragments

The get_fragment_atom_bond_set_combinations function below generates all fragment combinations in the range of [2 - (num_frags-1)] using the itertools.combinations() function. For example, if five fragments are generated for a molecule, the following number of fragment combinations are enumerated:

\(\dbinom{5}{2} + \dbinom{5}{3} + \dbinom{5}{4} = \frac{5!}{2!3!} + \frac{5!}{3!2!} + \frac{5!}{4!1!} = 10 + 10 + 5 = 25\)

However, only adjacent fragment combinations are kept, i.e., no disconnected molecule fragment is generated (see Table 1). If the fragments in a given combination are all adjacent, i.e., the is_adjacent_atom_bond_set_combination function returns True, then the atoms and bonds of these fragments are combined by calling the combine_and_connect_atom_bond_sets function, and the combined fragment is added to the list returned by get_fragment_atom_bond_set_combinations.

def get_fragment_atom_bond_set_combinations(
    frag_list: list[oechem.OEAtomBondSet],
) -> list[oechem.OEAtomBondSet]:
    """Generate all valid adjacent fragment combinations."""
    frag_combs: list[oechem.OEAtomBondSet] = []

    num_frags = len(frag_list)
    for n in range(2, num_frags):
        for frag_comb in combinations(frag_list, n):
            if is_adjacent_atom_bond_set_combination(list(frag_comb)):
                frag = combine_and_connect_atom_bond_sets(frag_comb)
                frag_combs.append(frag)

    return frag_combs

A fragment combination is considered adjacent only if all fragments are connected to each other within the combination. To determine this, the is_adjacent_atom_bond_set_combination function performs a depth-first search using the traverse_fragments function.

The num_parts variable keeps track of the number of connected components found in the fragment list. If all fragments are connected, num_parts will be one.

def is_adjacent_atom_bond_set_combination(
    frag_list: list[oechem.OEAtomBondSet],
) -> bool:
    """Check if a list of fragments forms a single connected component."""
    parts = [0] * len(frag_list)
    num_parts = 0

    for idx, frag in enumerate(frag_list):
        if parts[idx] != 0:
            continue

        num_parts += 1
        parts[idx] = num_parts
        traverse_fragments(frag, frag_list, parts, num_parts)

    return num_parts == 1

The traverse_fragments function is a recursive function that visits fragments connected to each other. The parts list is used not only to track whether a fragment has been visited, but also to record which connected component it belongs to.

def traverse_fragments(
    act_frag: oechem.OEAtomBondSet,
    frag_list: list[oechem.OEAtomBondSet],
    parts: list[int],
    num_parts: int,
) -> None:
    """Recursively traverse adjacent fragments to assign connected components."""
    for idx, frag in enumerate(frag_list):
        if parts[idx] != 0:
            continue

        if not is_adjacent_atom_bond_sets(act_frag, frag):
            continue

        parts[idx] = num_parts
        traverse_fragments(frag, frag_list, parts, num_parts)

The is_adjacent_atom_bond_sets function determines whether two fragments are connected. Two fragments are considered adjacent if an atom from each is connected by a bond.

def is_adjacent_atom_bond_sets(
    frag_a: oechem.OEAtomBondSet, frag_b: oechem.OEAtomBondSet
) -> bool:
    """Check if two atom/bond sets share an adjacent bond."""
    for atom_a in frag_a.GetAtoms():
        for atom_b in frag_b.GetAtoms():
            if atom_a.GetBond(atom_b) is not None:
                return True
    return False

If is_adjacent_atom_bond_set_combination returns True, i.e., the fragments are all connected, then combine_and_connect_atom_bond_sets can be used to combine the adjacent fragments by merging their atoms and bonds and adding the bonds that connect them.

def combine_and_connect_atom_bond_sets(
    frag_list: tuple[oechem.OEAtomBondSet, ...],
) -> oechem.OEAtomBondSet:
    """Combine fragment atom/bond sets and add connecting bonds."""
    combined = oechem.OEAtomBondSet()
    for frag in frag_list:
        for atom in frag.GetAtoms():
            combined.AddAtom(atom)
        for bond in frag.GetBonds():
            combined.AddBond(bond)

    for atom_a in combined.GetAtoms():
        for atom_b in combined.GetAtoms():
            if atom_a.GetIdx() < atom_b.GetIdx():
                continue

            bond = atom_a.GetBond(atom_b)
            if bond is None:
                continue
            if combined.HasBond(bond):
                continue

            combined.AddBond(bond)

    return combined

Download code

enumfrags.py

The following example shows the output fragment combinations of the molecule depicted in Figure 1.

Usage

See Download section to download the script.

> enumfrags --help
../_images/enumfrags-help.svg
> enumfrags --in molecule.ism --out .ism

Running the above command will generate the following output:

Input molecule: CC(C)NCC(COc1ccc(c(c1)Cc2ccccc2)CC(=O)N)O
5 fragments generated
11 fragment combinations generated
c1ccc(cc1)Cc2cccc(c2)O
c1cc(ccc1CC(=O)N)O
c1ccc(cc1)OCCO
CC(C)NCC(C)O
c1ccc(cc1)Cc2cc(ccc2CC(=O)N)O
c1ccc(cc1)Cc2cccc(c2)OCCO
c1cc(ccc1CC(=O)N)OCCO
CC(C)NCC(COc1ccccc1)O
c1ccc(cc1)Cc2cc(ccc2CC(=O)N)OCCO
CC(C)NCC(COc1cccc(c1)Cc2ccccc2)O
CC(C)NCC(COc1ccc(cc1)CC(=O)N)O

See also in Python documentation

  • itertools.combinations()

See also in OEChem TK manual

API

See also in OEMedChem TK manual

Theory

API

See also