Pharmacophore Utilities

Pharmacophore Utilities#

Modules for pharmacophore feature extraction and manipulation.

Pharmacophore container and priority labels#

get_pharmacophores() returns a Pharmacophore dataclass. For backwards compatibility it still unpacks and indexes as the original (types, positions, vectors) 3-tuple.

When built with return_atom_ids=True (or when priority_atoms is provided), each pharmacophore retains the set of contributing atom indices on .atom_ids. This enables priority_labels to compute 0/1 priority labels lazily against any number of priority-atom sets without re-extracting pharmacophores. When priority_atoms is passed directly to get_pharmacophores(), the labels are computed immediately and stored on .labels.

Labeling semantics (priority_pharm_labels()):

  • Point features (Acceptor, Donor, Halogen, Cation, Anion, ZnBinder): label 1 if any contributing atom appears in priority_atoms.

  • Aromatic and ring-derived Hydrophobe: require at least min_ring_priority_atoms heavy atoms (default 3) from a shared ring to also be in priority_atoms. This prevents a single priority atom on an aromatic ring from flagging the whole ring feature. Use min_ring_priority_atoms=1 for looser labeling.

The same options are available on get_pharmacophore(). Priority pharmacophore array indices (not atom indices) can be passed to ConditionalEval via priority_pharm_indices for subset Tversky scoring after full-set pharmacophore alignment.

Pharmacophore Extraction#

class shepherd_score.pharm_utils.pharmacophore.Pharmacophore(types, positions, vectors, mol=None, atom_ids=None, labels=None)[source]#

Bases: object

Container for the pharmacophores extracted from a molecule.

For backwards compatibility with the original 3-tuple return of get_pharmacophores(), instances unpack and index as (types, positions, vectors):

X, P, V = get_pharmacophores(mol)   # still works

Named attributes (.types, .positions, .vectors) are also available. When built with return_atom_ids=True, the per-pharmacophore atom-id sets are retained on .atom_ids and priority labels can be computed lazily — for any number of priority-atom sets or thresholds — without re-extracting pharmacophores.

Parameters:
types#

Pharmacophore type indices (P_TYPES order), shape (N,).

Type:

np.ndarray

positions#

Anchor positions, shape (N, 3).

Type:

np.ndarray

vectors#

Relative unit vectors, shape (N, 3).

Type:

np.ndarray

mol#

Source molecule, needed for ring-aware priority_labels().

Type:

rdkit.Chem.Mol or None

atom_ids#

Per-pharmacophore atom-id sets aligned with types; None unless the container was built with return_atom_ids=True.

Type:

list of set or None

labels#

Priority labels aligned with types, populated when get_pharmacophores is called with priority_atoms; None otherwise.

Type:

np.ndarray or None

priority_labels(priority_atoms, min_ring_priority_atoms=3)[source]#

Compute a 0/1 priority label per pharmacophore against priority_atoms.

Requires the container to have been built with return_atom_ids=True. See priority_pharm_labels() for the labeling semantics.

Parameters:

  • priority_atoms (iterable of int) – Atom indices considered “priority”.

  • min_ring_priority_atoms (int, optional) – Minimum heavy ring atoms in priority_atoms before an aromatic or aromatic-derived hydrophobe is labeled 1. Default is 3.

Return type:

np.ndarray, shape (N,), dtype int64

Generate pharmacophores from a RDKit conformer.

Parts of code adapted from Francois Berenger / Tsuda Lab and RDKit.

References:

shepherd_score.pharm_utils.pharmacophore.pattern_of_smarts(s)[source]#
shepherd_score.pharm_utils.pharmacophore.find_hydrophobes(mol, cluster_hydrophobic=True, return_atom_ids=False)[source]#

Find hydrophobes and cluster them.

Parameters:

  • mol (rdkit Mol object with a conformer.)

  • cluster_hydrophobic (bool (default=True) to cluster hydrophobic atoms if they fall within 2A.)

  • return_atom_ids (bool (default=False)) – When True, returns a list of (center_tuple, aromatic_atom_ids_set) pairs instead of plain center tuples. aromatic_atom_ids_set is the union of atom indices from aromatic-ring matches that ended up in the cluster; it is an empty set for clusters composed entirely of non-aromatic matches.

Returns:

  • list of tuples containing coordinates for the locations for each hydrophobe,

  • or (when return_atom_ids=True) a list of (center, set[int]) pairs.

shepherd_score.pharm_utils.pharmacophore.get_pharmacophores_dict(mol, multi_vector=True, exclude=[], check_access=False, scale=1.0, return_atom_ids=False)[source]#

Get the positions of pharmacophore anchors and their associated unit vectors.

Returns a dictionary. Adapted from rdkit.Chem.Features.ShowFeats.ShowMolFeats.

Parameters:

  • mol (rdkit.Chem.Mol) – RDKit Mol object with a conformer.

  • multi_vector (bool, optional) – Whether to represent pharmacophores with multiple vectors. Default is True.

  • exclude (list, optional) – List of atom indices to not include as a HBD. Default is [].

  • check_access (bool, optional) – Check if HBD/HBA are accessible to the molecular surface. Default is False.

  • scale (float, optional) – Length of the vector in Angstroms. Default is 1.0.

  • return_atom_ids (bool, optional) – When True, each family sub-dict also contains an 'A' key holding a list of atom-id sets (one set per emitted pharmacophore, aligned with 'P'). For hydrophobes, only aromatic-ring-derived clusters carry non-empty sets; aliphatic clusters have empty sets. Default is False.

Returns:

Dictionary with format {'FeatureName': {'P': [(anchor coord), ...], 'V': [(rel. vec), ...]}}. When return_atom_ids=True, each entry also has 'A': [set_of_int, ...] aligned with 'P'.

Return type:

dict

shepherd_score.pharm_utils.pharmacophore.priority_pharm_labels(mol, atom_ids_per_pharm, pharm_types, priority_atoms, min_ring_priority_atoms=3)[source]#

Compute a 0/1 priority label for each pharmacophore.

Parameters:

  • mol (rdkit.Chem.Mol) – Molecule used to resolve ring membership for aromatic/hydrophobe labels.

  • atom_ids_per_pharm (list of sets) – One set of atom indices per pharmacophore, aligned with the X/P/V arrays.

  • pharm_types (np.ndarray) – Pharmacophore type indices aligned with atom_ids_per_pharm (same order as P_TYPES).

  • priority_atoms (iterable of int) – Atom indices considered “priority”.

  • min_ring_priority_atoms (int, optional) – Minimum number of heavy ring atoms that must also be in priority_atoms before an aromatic or aromatic-derived hydrophobe is labeled 1. Use 1 to treat any single priority atom in the ring as sufficient. Default is 3.

Returns:

1 where the pharmacophore is priority, else 0. Non-ring pharmacophores use simple atom-id intersection. Aromatic and aromatic-derived hydrophobe pharmacophores additionally require at least min_ring_priority_atoms heavy atoms from a shared ring to appear in priority_atoms.

Return type:

np.ndarray, shape (N,), dtype int64

shepherd_score.pharm_utils.pharmacophore.get_pharmacophores(mol, multi_vector=True, exclude=[], check_access=False, scale=1.0, return_atom_ids=False, priority_atoms=None, min_ring_priority_atoms=3)[source]#

Get the identity, anchor positions, and relative unit vectors for each pharmacophore.

Pharmacophore ordering for indexing: (‘Acceptor’, ‘Donor’, ‘Aromatic’, ‘Hydrophobe’, ‘Cation’, ‘Anion’, ‘ZnBinder’)

Notes

The check_access parameter is currently based on whether interaction points sampled from a sphere’s surface with a radius of 1.8A from the acceptor/donor atom falls outside the solvent accessible surface defined by the vdW radius + 0.8A of the neighboring atoms. This works for buried acceptors/donors, but may be prone to false positives. For example, CN(C)C would have its sole HBA rejected. Other approaches such as buried volume should be considered in the future.

Parameters:

  • mol (rdkit.Chem.Mol) – RDKit Mol object with conformer.

  • multi_vector (bool, optional) – Whether to represent pharmacophores with multiple vectors. Default is True.

  • exclude (list, optional) – List of hydrogen indices to not include as a HBD. Default is [].

  • check_access (bool, optional) – Check if HBD/HBA are accessible to the molecular surface. Default is False.

  • scale (float, optional) – Length of a pharmacophore vector in Angstroms. Default is 1.0.

  • return_atom_ids (bool, optional) – When True, the returned Pharmacophore retains the per-pharmacophore atom-id sets on .atom_ids, enabling lazy priority labeling via Pharmacophore.priority_labels(). Implied by priority_atoms. Default is False.

  • priority_atoms (iterable of int, optional) – Atom indices considered “priority”. When provided, priority labels are computed in this single call and stored on .labels of the returned container (atom ids are retained automatically). See Pharmacophore.priority_labels() for the labeling semantics. Default is None.

  • min_ring_priority_atoms (int, optional) – Only used when priority_atoms is provided. Minimum heavy ring atoms in priority_atoms before an aromatic or aromatic-derived hydrophobe is labeled 1. Set to 1 to label any pharmacophore whose ring shares a single priority atom. Default is 3.

Returns:

Container that unpacks as the original (X, P, V) 3-tuple for backwards compatibility, where:

  • X (.types): pharmacophore type indices (P_TYPES order), shape (N,).

  • P (.positions): anchor positions, shape (N, 3).

  • V (.vectors): relative unit vectors, shape (N, 3); adding P and V gives the extended point.

When return_atom_ids=True (or priority_atoms is given), .atom_ids holds per-pharmacophore atom-id sets and .priority_labels(priority_atoms, ...) computes 0/1 priority labels. When priority_atoms is given, .labels is also populated in this call.

Return type:

Pharmacophore

Pharmacophore Vectors#

Generate the vector features for pharmacophores from a rdkit conformer.

Adapted from rdkit:

rdkit/rdkit rdkit/rdkit

Changed to return anchor position and relative unit vector.

shepherd_score.pharm_utils.pharmvec.GetAromaticFeatVects(conf, featAtoms, featLoc, return_both=False, scale=1.0)[source]#

Compute the direction vector for an aromatic feature

Changed: only return one vector, process later for visualization and scoring

Parameters:

  • conf (a conformer)

  • featAtoms (list of atom IDs that make up the feature)

  • featLoc (location of the aromatic feature specified as point3d)

  • return_both (bool for whether to return both vectors or just one.)

  • scale (the size of the direction vector)

Returns:

list of anchor position(s) as rdkit Point3D list of relative unit vector(s) as rdkit Point3D

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetDonorFeatVects(conf, featAtoms, scale=1.0, exclude=[])[source]#

Get vectors for hydrogen bond donors in the direction of the hydrogens.

Parameters:

  • conf (rdkit Mol object with a conformer.)

  • featAtoms (list containing rdkit Atom object of atom attributed as a donor.)

  • scale (float (default = 1.) length of direction vector.)

  • exclude (list of atom indices that should not be included as a donatable H.)

Returns:

list of anchor position(s) as rdkit Point3D or None list of relative unit vector(s) as rdkit Point3D or None list of neighboring hydrogens or None

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptorFeatVects(conf, featAtoms, scale=1.0)[source]#

Get the anchor positions and relative unit vectors of an acceptor atom.

Assumes HBA’s are only O and N as defined by smarts_features.fdef. If HBA is not one of those, then it assumes the atom has one lone pair.

Parameters:

  • conf (Chem.Mol) – RDKit Mol object with a conformer.

  • featAtoms (list) – List containing RDKit Atom object of atom attributed as an acceptor.

  • scale (float, optional) – Length of direction vector. Default is 1.0.

Returns:

(list of anchor position(s) as RDKit Point3D or [None], list of relative unit vector(s) as RDKit Point3D or [None])

Return type:

tuple

shepherd_score.pharm_utils.pharmvec.GetHalogenFeatVects(conf, featAtoms, scale=1.0)[source]#

Get the anchor positions and relative unit vectors of a halogen atom. Assumes only one connection.

Parameters:

  • conf (rdkit Mol object with a conformer.)

  • featAtoms (list containing rdkit Atom object of atom attributed as an acceptor.)

  • scale (float (default = 1.) length of direction vector.)

Returns:

list of anchor position(s) as rdkit Point3D or [None] list of relative unit vector(s) as rdkit Point3D or [None]

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetDonor1FeatVects_single(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Donor of type 1. Made to generate a single vector representation.

This is a donor with one heavy atom. It is not clear where we should we should be putting the direction vector for this. It should probably be a cone. In this case we will just use the direction vector from the donor atom to the heavy atom.

Changed: conditioning based on the number of hydrogens 1. If 1 hydrogen, vector should point in the direction of the hydrogen. 2. If 2 hydrogens, vector should point in a bisecting direction of the two hydrogens. 3. If 3 hydrogens, point in the direction of the bond.

Parameters:

  • conformer (conf - rdkit Mol object with)

  • feature (featAtoms - list of atoms that are part of the)

  • 1.0) (scale - float for length of the direction vector (default =)

Returns:

anchor position as rdkit Point3D or None relative unit vector(s) as rdkit Point3D or None list of hydrogen rdkit Atom objects

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetDonor2FeatVects_single(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Donor of type 2. Made to generate a single vector representation.

This is a donor with two heavy atoms as neighbors. The atom may are may not have hydrogen on it. Here are the situations with the neighbors that will be considered here 1. two heavy atoms and two hydrogens: we will assume a sp3 arrangement here 2. two heavy atoms and one hydrogen: this can either be sp2 or sp3 3. two heavy atoms and no hydrogens

Changed: conditioning based on the number of hydrogens 1. For case 1, point in the direction bisecting the two hydrogens. 2. For case 2, point in the direction of the hydrogen. 3. For case 3, no changes.

Parameters:

  • conf (rdkit Mol object with conformer)

  • featAtoms (list of atoms that are part of the feature)

  • scale (float for length of the direction vector (default = 1.0))

Returns:

anchor position as rdkit Point3D or None relative unit vector(s) as rdkit Point3D or None list of hydrogen rdkit Atom objects

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetDonor3FeatVects_single(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Donor of type 3. Made to generate a single vector representation.

This is a donor with three heavy atoms as neighbors. We will assume a tetrahedral arrangement of these neighbors. So the direction we are seeking is the last fourth arm of the sp3 arrangement

Changed: Return anchor and relative unit vector tuple

Parameters:

  • conf (rdkit Mol object with conformer)

  • featAtoms (list of atoms that are part of the feature)

  • scale (float for length of the direction vector (default = 1.0))

Returns:

anchor position as rdkit Point3D or None relative unit vector(s) as rdkit Point3D or None list of hydrogen rdkit Atom objects

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptor1FeatVects_single(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Acceptor of type 1 (single vector representation).

This is an acceptor with one heavy atom neighbor. There are two possibilities:

  • The bond to the heavy atom is a single bond (e.g. CO): We use the inversion of this bond direction and mark it as a ‘cone’.

  • The bond to the heavy atom is a double bond (e.g. C=O): We have two possible directions except in some special cases (e.g. SO2) where we use bond direction.

Notes

Modified to condition on the number of hydrogens with methanamine fix:

  • Case 1: If one hydrogen, vector points in the opposite direction of the bisection of the acute angle formed by the heavy-acceptor-hydrogen. If two hydrogens, assume sp3 and project in that lone-pair direction. If not tetrahedral, return None.

  • Case 2: Return the bisecting vector of the two lone-pairs.

Parameters:

  • conf (Chem.Mol) – RDKit Mol object with conformer.

  • featAtoms (list) – List of atoms that are part of the feature.

  • scale (float, optional) – Length of the direction vector. Default is 1.0.

Returns:

(anchor position as RDKit Point3D or None, relative unit vector(s) as RDKit Point3D or None)

Return type:

tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptor2FeatVects_single(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Acceptor of type 2. Made to generate a single vector representation.

This is the acceptor with two adjacent heavy atoms. We will special case a few things here. If the acceptor atom is an oxygen we will assume a sp3 hybridization the acceptor directions (two of them) reflect that configurations. Otherwise the direction vector in plane with the neighboring heavy atoms

Changed: Only generate one vector Rather than generating 2 vectors for sp3 oxygen, just keep the average vector.

Parameters:

  • conf (rdkit Mol object with conformer)

  • featAtoms (list of atoms that are part of the feature)

  • scale (float for length of the direction vector (default = 1.0))

Returns:

anchor position as rdkit Point3D or None relative unit vector(s) as rdkit Point3D or None

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptor3FeatVects_single(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Donor of type 3. Made to generate a single vector representation.

This is a donor with three heavy atoms as neighbors. We will assume a tetrahedral arrangement of these neighbors. So the direction we are seeking is the last fourth arm of the sp3 arrangement

Changed: to return anchor and relative unit vector tuple

Parameters:

  • conf (rdkit Mol object with conformer)

  • featAtoms (list of atoms that are part of the feature)

  • scale (float for length of the direction vector (default = 1.0))

Returns:

anchor position as rdkit Point3D or None relative unit vector(s) as rdkit Point3D or None

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptor1FeatVects(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Acceptor of type 1 (multi-vector representation).

This is an acceptor with one heavy atom neighbor. There are two possibilities:

  • The bond to the heavy atom is a single bond (e.g. CO): We use the inversion of this bond direction and mark it as a ‘cone’.

  • The bond to the heavy atom is a double bond (e.g. C=O): We have two possible directions except in some special cases (e.g. SO2) where we use bond direction.

Notes

Modified to change return format, with fixes for methanamine and two vectors for hydroxyls.

Parameters:

  • conf (Chem.Mol) – RDKit Mol object with a conformer.

  • featAtoms (list) – List containing RDKit Atom object of atom attributed as an acceptor.

  • scale (float, optional) – Length of direction vector. Default is 1.0.

Returns:

(list of anchor position(s) as RDKit Point3D or None, list of relative unit vector(s) as RDKit Point3D or None)

Return type:

tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptor2FeatVects(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Acceptor of type 2. Made to generate a single vector representation.

This is the acceptor with two adjacent heavy atoms. We will special case a few things here. If the acceptor atom is an oxygen we will assume a sp3 hybridization the acceptor directions (two of them) reflect that configurations. Otherwise the direction vector in plane with the neighboring heavy atoms

Changed: return format

Parameters:

  • conf (rdkit Mol object with a conformer.)

  • featAtoms (list containing rdkit Atom object of atom attributed as an acceptor.)

  • scale (float (default = 1.) length of direction vector.)

Returns:

list of anchor position(s) as rdkit Point3D or None list of relative unit vector(s) as rdkit Point3D or None

Return type:

Tuple

shepherd_score.pharm_utils.pharmvec.GetAcceptor3FeatVects(conf, featAtoms, scale=1.0)[source]#

Get the direction vectors for Donor of type 3. Made to generate a single vector representation.

This is a donor with three heavy atoms as neighbors. We will assume a tetrahedral arrangement of these neighbors. So the direction we are seeking is the last fourth arm of the sp3 arrangement

Changed: return format

Parameters:

  • conf (rdkit Mol object with a conformer.)

  • featAtoms (list containing rdkit Atom object of atom attributed as an acceptor.)

  • scale (float (default = 1.) length of direction vector.)

Returns:

list of anchor position(s) as rdkit Point3D or None list of relative unit vector(s) as rdkit Point3D or None

Return type:

Tuple