Finding Core Fragment of a Molecule Series

Problem

You want to find the core fragment in a molecule series i.e. find the largest common substructure of the molecules in a series. See example in Figure 1.

../_images/core-example.png — **Figure 1: Example maximum common subgraph of a molecule set**

Ingredients

OEChem TK - cheminformatics toolkit
GraphSim TK - fingerprinting toolkit

Difficulty Level

🌶️ 🌶️

Download

Download code

get_common_core.py and supporting file test.ism

See also Usage subsection.

Source Code

get_common_core#!/usr/bin/env python3
# (C) 2023 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.

"""Identify the "best" i.e. largest common fragment of a set of molecules."""

import argparse
import os
import pathlib
import sys

import rich.console
from openeye import oechem, oegraphsim
from rich_argparse import HelpPreviewAction, RichHelpFormatter

__SCRIPT_NAME__ = pathlib.Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Extract common largest core of molecules."
__SCRIPT_TOOLKITS__ = ["oechem", "oegraphsim"]
__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]",
    )
    parser.add_argument("--help-image", action=HelpPreviewAction)

    # input options
    input_group = parser.add_argument_group("Input options")
    input_group.add_argument(
        "-in",
        "--input",
        type=str,
        required=True,
        metavar="MOL-FILE",
        help="Input molecule file.",
    )
    core_group = parser.add_argument_group("Core identifying options")
    core_group.add_argument(
        "--min-bonds",
        type=int,
        required=False,
        default=5,
        choices=range(1, 10),
        help="Minimum number of bonds in common core",
    )
    core_group.add_argument(
        "--max-bonds",
        type=int,
        required=False,
        default=18,
        choices=range(4, 18),
        help="Maximum number of bonds in common core",
    )
    return parser.parse_args()


def main() -> int:
    """Identify common fragments."""
    args = parse_options()

    # check input file
    ifs = oechem.oemolistream()
    if not ifs.open(args.input):
        oechem.OEThrow.Fatal("Cannot open input file!")

    # read molecules
    mol_list: list[oechem.OEMolBase] = []
    for mol in ifs.GetOEGraphMols():
        oechem.OESuppressHydrogens(mol)
        mol_list.append(oechem.OEGraphMol(mol))

    console = rich.console.Console()
    console.print(f"Number of molecules = {len(mol_list)}")

    if (
        best_core := get_largest_common_core(mol_list, args.min_bonds, args.max_bonds)
    ) is not None:
        console.print(
            f"Largest core fragment of the molecules = {oechem.OEMolToSmiles(best_core)}"
        )
    else:
        console.print("No core fragment is identified!")
    return os.EX_OK


def get_reference_molecule(mol_list: list[oechem.OEMolBase]) -> oechem.OEMolBase | None:
    """Return a reference "smallest" molecule."""
    ref_mol = None
    for mol in mol_list:
        if ref_mol is None or mol.NumAtoms() < ref_mol.NumAtoms():
            ref_mol = mol
    return ref_mol


def get_fragments(
    mol: oechem.OEMolBase, min_bonds: int, max_bonds: int
) -> list[oechem.OEMolBase]:
    """Return a list of all possible fragments of a molecule."""
    frags: list[oechem.OEMolBase] = []
    fp_type = oegraphsim.OEGetFPType(
        f"Tree,ver=2.0.0,size=4096,bonds={min_bonds}-{max_bonds},atype=AtmNum,btype=Order"
    )

    unique_fragments = True
    for ab_set in oegraphsim.OEGetFPCoverage(mol, fp_type, unique_fragments):
        frag_atom_pred = oechem.OEIsAtomMember(ab_set.GetAtoms())

        frag = oechem.OEGraphMol()
        adjust_hcount = True
        oechem.OESubsetMol(frag, mol, frag_atom_pred, adjust_hcount)
        oechem.OEFindRingAtomsAndBonds(frag)
        frags.append(oechem.OEGraphMol(frag))

    return frags


def get_core_fragments(
    mol_list: list[oechem.OEMolBase],
    frags: list[oechem.OEMolBase],
    atom_expr: int = oechem.OEExprOpts_DefaultAtoms,
    bond_expr: int = oechem.OEExprOpts_DefaultBonds,
) -> list[oechem.OEMolBase]:
    """Identify fragments that are present in all molecules."""
    core_frags: list[oechem.OEMolBase] = []
    for frag in frags:
        ss = oechem.OESubSearch(frag, atom_expr, bond_expr)
        if not ss.IsValid():
            continue

        valid_core = True
        for mol in mol_list:
            valid_core = ss.SingleMatch(mol)
            if not valid_core:
                break

        if valid_core:
            core_frags.append(frag)

    return core_frags


def get_largest_common_core(
    mol_list: list[oechem.OEMolBase],
    min_bonds: int,
    max_bonds: int,
    atom_expr: int = oechem.OEExprOpts_DefaultAtoms,
    bond_expr: int = oechem.OEExprOpts_DefaultBonds,
) -> oechem.OEMolBase | None:
    """Enumerate all the fragments to find a largest common core."""
    ref_mol = get_reference_molecule(mol_list)
    if ref_mol is None:
        oechem.OEThrow.Fatal("No reference molecule identified.")
        return None

    frags = get_fragments(ref_mol, min_bonds, max_bonds)
    if len(frags) == 0:
        oechem.OEThrow.Fatal(
            f"No fragment is enumerated with bonds {min_bonds}-{max_bonds}!"
        )

    common_frags = get_core_fragments(mol_list, frags, atom_expr, bond_expr)
    if len(common_frags) == 0:
        oechem.OEThrow.Error("No common fragment is found!")

    largest_core = None
    for frag in common_frags:
        if largest_core is None or get_fragment_score(
            largest_core
        ) < get_fragment_score(frag):
            largest_core = frag
    return largest_core


def get_fragment_score(frag: oechem.OEMolBase) -> float:
    """Score the fragment (larger is better)."""
    score = 0.0
    score += 2.0 * oechem.OECount(frag, oechem.OEAtomIsInRing())
    score += 1.0 * oechem.OECount(frag, oechem.OENotAtom(oechem.OEAtomIsInRing()))
    return score


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 obvious choice to solve this problem is to use a maximum common substructure (MCS) search algorithm.

Definition:: The maximum common subgraph is the largest possible common subgraph of two graphs, i.e. it can not be extended to another common subgraph by the addition of vertices or edges.

The OEChem TK provides the OEMCSSearch class that solves the MCS problem for two molecules. See example in Figure 2.

../_images/mcss.png — **Figure 2: Example maximum common subgraph of two molecules**

MCS is commonly used in cheminformatics to assess molecule similarity. One advantage of MCS over the fingerprinting methods is that it provides an atom - atom mapping between two molecules that is essential to solve this problem.

However, MCS identification belongs to the class of NP-complete problems for which no algorithm with polynomial-time complexity is known for the general case. The identification of all subgraphs containing k nodes that are common to two graphs containing n and m nodes, respectively, requires

\(\frac{m!n!}{(m-k)!(n-k)!k!}\)

atom-by-atom comparisons. The identification of the largest common subgraph is achieved by carrying out this set of comparisons with \(\forall k : 1 < k < min(m, n)\) until it is not possible to identify a larger common substructure. Even though chemical structures do not have high connectivity, which would lead to exponential behavior, it still can be a very expensive operation for large molecules.

The other problem is that the MCS search algorithm implemented in the OEMCSSearch class only works for two molecules. Identifying a MCS for a molecule set would require performing \(n x n\) comparisons. Combining the returned \(n x n\) MCS atom mappings is not a trivial problem either. See example in Figure 3.

../_images/core-mcss-ABC.png — **Figure 3: Example maximum common subgraph of multiple molecules**

However the problem we want to solve is quite specific that does not require to perform \(n x n\) MCS searches. If we can assume that a common core does exist in the molecule series than the following process can identify this core fragment:

Exhaustively enumerate all fragments in an arbitrary molecule of the series. Choosing a small molecule is preferred since it will have fewer fragments.
Perform a substructure search with each fragment to verify whether it is a common fragment, i.e. a substructure of each molecule in the series.
Chose the largest common fragment as the core.

The first step of finding a core of a molecule series is to select an arbitrary reference molecule. The get_reference_molecule function simply selects the molecule with the smallest number of atoms.

def get_reference_molecule(mol_list: list[oechem.OEMolBase]) -> oechem.OEMolBase | None:
    """Return a reference "smallest" molecule."""
    ref_mol = None
    for mol in mol_list:
        if ref_mol is None or mol.NumAtoms() < ref_mol.NumAtoms():
            ref_mol = mol
    return ref_mol

The get_fragments function is then called with the reference molecule to enumerate all of its fragments by using functions of GraphSim TK. When generating fingerprints, a molecule graph is exhaustively traversed, enumerating various fragments (using the path, circular or tree method), and then hashing them into a fixed-length bit-vector. The OEGetFPCoverage function of GraphSim TK provides access to these fragments by returning an iterator over OEAtomBondSet objects, each of which stores the atoms and bonds of a specific fragment. See example in Figure 4 that shows how tree fragments are enumerated up to a given length.

../_images/TreeEnumeration.png — **Figure 4: Example of enumerating tree fragments with various lengths**

The get_fragments function first generates a tree fingerprint type from a string representation using the OEGetFPType function. Then it loops over the OEAtomBondSet objects returned by the OEGetFPCoverage function and extracts the corresponding molecule fragment using the OESubsetMol function.

def get_fragments(
    mol: oechem.OEMolBase, min_bonds: int, max_bonds: int
) -> list[oechem.OEMolBase]:
    """Return a list of all possible fragments of a molecule."""
    frags: list[oechem.OEMolBase] = []
    fp_type = oegraphsim.OEGetFPType(
        f"Tree,ver=2.0.0,size=4096,bonds={min_bonds}-{max_bonds},atype=AtmNum,btype=Order"
    )

    unique_fragments = True
    for ab_set in oegraphsim.OEGetFPCoverage(mol, fp_type, unique_fragments):
        frag_atom_pred = oechem.OEIsAtomMember(ab_set.GetAtoms())

        frag = oechem.OEGraphMol()
        adjust_hcount = True
        oechem.OESubsetMol(frag, mol, frag_atom_pred, adjust_hcount)
        oechem.OEFindRingAtomsAndBonds(frag)
        frags.append(oechem.OEGraphMol(frag))

    return frags

Note

When using the OEGetFPCoverage function only fragment size parameter of the fingerprint type is considered. The other parameters such as atom and bond typing are utilized only when hashing the enumerated fragment into a fixed-length bit-vector.

Discussion

Warning

The above algorithm is unable to identify disconnected MCS of multiple molecules, neither does the algorithm that is implemented in the OEMCSSearch class for two molecules. See example in Figure 5.
The above algorithm is also not suitable to cluster molecules with more than one core fragments. This is a more complicate problem that can be tackled with performing n x n MCS searches.

../_images/disconnected-mcs.png — **Figure 5: Example of disconnected MCS of two molecules**

When performing substructure searches in the get_largest_common_core function the following two options are used: (see Example (A) in Table 1)

The OEExprOpts_DefaultAtoms option means that two atoms are considered to be equivalent if they have the same atomic number, aromaticity, and formal charge.
The OEExprOpts_DefaultBonds option means that two bonds can be mapped to each other if they have the same bond order and aromaticity.

By modifying the atom and bond expression options, very diverse pattern matching can be performed. Example(B) - Example (D) in Table 1 show examples where the discrimination capability of the OEExprOpts_DefaultAtoms and OEExprOpts_DefaultBonds options are decreased by using various modifiers. This results in identifying larger and larger “common” core fragments. For example, using the OEExprOpts_EqAromatic modifier, atoms in any aromatic ring systems are considered equivalent. As a result, a larger core fragment is identified since the pyridine and pyrimidine rings are considered equivalent.

**Table 1. Examples of the effects of using various atom and bond expressions**

Example (A) – DefaultAtoms and DefaultBonds

Example (B) – DefaultAtoms \| EqAromatic and DefaultBonds

Example (C) – DefaultAtoms \| EqHalogen and DefaultBonds

Example (D) – DefaultAtoms \| EqAromatic \| EqCAliphaticONS and DefaultBonds \| EqSingleDouble

Usage

See Download section to download the script and test data.

> get_common_core --help

The following example shows the output of the get_common_core script for the example molecule file get_common_core.ism:

> get_common_core --input get_common_core.ism

Number of molecules = 459
Largest core fragment of the molecules = Cc1cc2cc(ccc2nc1N)O

Finding Core Fragment of a Molecule Series

Problem

Ingredients

Difficulty Level

Download

Source Code

Solution

Discussion

Usage

See also in OEChem TK manual

See also in GraphSim TK manual

See also