<jupyter_start><jupyter_text>T009 · Ligand-based pharmacophoresAuthors:- <NAME>, CADD seminar, 2017, Charité/FU Berlin- <NAME>, CADD seminar, 2018, Charité/FU Berlin- <NAME>, [Volkamer lab](https://volkamerlab.org/), Charité- <NAME>, [Volkamer lab](https://volkamerlab.org/), Charité __Talktorial T009__: This talktorial is part of the TeachOpenCADD pipeline described in the [first TeachOpenCADD paper](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0351-x), comprising of talktorials T001-T010. **Note**: Please run this notebook cell by cell. Running all cells in one is possible also, however, part of the nglview 3D representations might be missing.  Aim of this talktorialIn this talktorial, we use known EGFR ligands, which were selected and aligned in the previous talktorial, to identify donor, acceptor, and hydrophobic pharmacophoric features for each ligand. Those features are then clustered to define an ensemble pharmacophore, which represents the properties of the set of known EGFR ligands and can be used to search for novel EGFR ligands via virtual screening.  Learning goals  Contents in *Theory** Pharmacophore modeling  * Structure- and ligand-based pharmacophore modeling* Virtual screening with pharmacophores* Clustering: k-means  Contents in *Practical** Get pre-aligned ligands from previous talktorial* Show ligands with NGLView* Extract pharmacophore features* Show the pharmacophore features of all ligands  * Hydrogen bond donors  * Hydrogen bond acceptors  * Hydrophobic contacts* Collect coordinates of features per feature type* Generate ensemble pharmacophores  * Set static parameters for k-means clustering  * Set static parameters for cluster selection  * Define k-means clustering and cluster selection functions  * Cluster features  * Select relevant clusters  * Get selected cluster coordinates* Show clusters  * Hydrogen bond donors  * Hydrogen bond acceptors  * Hydrophobic contacts* Show ensemble pharmacophore  References* IUPAC pharmacophore definition ([Pure & Appl. Chem (1998), 70, 1129-43](https://www.degruyter.com/view/journals/pac/70/5/article-p1129.xml))* 3D pharmacophores in LigandScout ([J. Chem. Inf. Model. (2005), 45, 160-9](http://pubs.acs.org/doi/pdf/10.1021/ci049885e))* Book chapter: Pharmacophore Perception and Applications ([Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), **1**, 259-82](https://onlinelibrary.wiley.com/doi/10.1002/9783527806539.ch6f))* Book chapter: Structure-Based Virtual Screening ([Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), **1**, 313-31](https://onlinelibrary.wiley.com/doi/10.1002/9783527806539.ch6h)).* <NAME> and the origin of the pharmacophore concept ([Internet Electron. J. Mol. Des. (2007), 6, 271-9](http://biochempress.com/Files/iejmd_2007_6_0271.pdf))* <NAME>'s demonstration of pharmacophore modeling with RDKit ([RDKit UGM 2016 on GitHub](https://github.com/rdkit/UGM_2016/blob/master/Notebooks/Stiefl_RDKitPh4FullPublication.ipynb))   Theory  PharmacophoresIn computer-aided drug design, the description of drug-target interactions with pharmacophores is a well-established method. The term pharmacophore was defined in 1998 by a IUPAC working party:"A pharmacophore is the ensemble of steric and electronic features that is necessary to ensure the optimal supramolecular interactions with a specific biological target structure and to trigger (or to block) its biological response." ([Pure & Appl. Chem. (1998), 70, 1129-43](https://www.degruyter.com/view/journals/pac/70/5/article-p1129.xml))In other words, a *pharmacophore* consists of several *pharmacophoric features*, which describe important steric and physico-chemical properties of a ligand observed to bind a target under investigation. Such *physico-chemical properties* (also called feature types) can be hydrogen bond donors/acceptors, hydrophobic/aromatic interactions, or positively/negatively charged groups, and the *steric properties* are defined by the 3D arrangement of these features.  Structure- and ligand-based pharmacophore modelingIn pharmacophore modeling, two main approaches are used, depending on the biological question and available data sources, i.e. structure- and ligand-based pharmacophore modeling.*Structure-based pharmacophore models* are derived from protein-ligand complexes. Features are defined by observed interactions between the protein and ligand, ensuring that only those ligand moieties are used for virtual screening that have already been shown to be involved in ligand binding. However, structures of protein-ligand complexes are not available for all targets. In this case, either complex structures can be generated by modeling the ligand into the target binding site, e.g. via molecular docking, or pharmacophore modeling methods can be invoked that only use the target binding site to detect potential protein-ligand interaction sites.*Ligand-based pharmacophore models* are based on a set of ligands known to bind the target under investigation. The common chemical features of these ligands build the pharmacophore model. This method is used for targets with multiple known ligands and in case of missing protein-ligand complex structures. In this talktorial, we will use ligand-based pharmacophore modeling using a set of known EGFR ligands.For more information on pharmacophore modeling, we recommend [Pharmacophore Perception and Applications: Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), **1**, 259-82](https://onlinelibrary.wiley.com/doi/10.1002/9783527806539.ch6f) and [J. Chem. Inf. Model. (2005), 45, 160-9](http://pubs.acs.org/doi/pdf/10.1021/ci049885e).![Structure-based pharmacophore representing protein-ligand interactions (figure by <NAME>)](images/target+ligand+pharma_3x_feather.png)*Figure 1*: Structure-based pharmacophore representing protein-ligand interactions. Figure created by <NAME>.  Virtual screening with pharmacophoresAs described earlier in **Talktorial 4**, virtual screening (VS) describes the screening of a query, e.g. a pharmacophore model (T009) or a query compound (T004), against a large library of compounds, in order to identify the small molecules in the library that are most likely to bind a target under investigation (represented by the query). In pharmacophore-based virtual screening, the compound library is matched compound-by-compound into a pharmacophore model and ranked by the best matching results ([Structure-Based Virtual Screening: Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), **1**, 313-31](https://onlinelibrary.wiley.com/doi/10.1002/9783527806539.ch6h)).  Clustering: k-meansIn this talktorial we will generate an ensemble pharmacophore by clustering the feature points of several ligand-based pharmacophores. The clustering algorithm uses the k-means clustering, which aims to cluster a data set into k clusters:1. k different centroids are selected and each point of the data set is assigned to its closest centroids.2. New centroids are calculated based on the current clusters and each point of the data set is newly assigned to its closest centroids.3. This procedure is repeated until the centroids are stable.More info about [k-means in Wikipedia](https://de.wikipedia.org/wiki/K-Means-Algorithmus).  Practical<jupyter_code>import collections
from collections import Counter
from pathlib import Path
import operator
import time
import math

import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets, cluster

from rdkit import RDConfig, Chem, Geometry, DistanceGeometry
from rdkit.Chem import (
    ChemicalFeatures,
    rdDistGeom,
    Draw,
    rdMolTransforms,
    AllChem,
)
from rdkit.Chem.Draw import DrawingOptions
from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib
from rdkit.Numerics import rdAlignment
import nglview as nv<jupyter_output><empty_output><jupyter_text>Set the Path to the current notebook.<jupyter_code>HERE = Path(_dh[-1])
DATA = HERE / "data"<jupyter_output><empty_output><jupyter_text>Get pre-aligned ligands from previous talktorialWe retrieve all ligands that were aligned in the previous talktorial.First, we get the file paths to all ligand PDB files.<jupyter_code>mol_files = sorted(Path(HERE / "../T008_query_pdb/data/").glob("*_lig.pdb"))
print([f.name for f in mol_files])
# NBVAL_CHECK_OUTPUT
pdb_ids = [Path(file).stem.split("_")[0] for file in mol_files]
pdb_ids<jupyter_output><empty_output><jupyter_text>Second, we read all ligands from these PDB files using RDKit.<jupyter_code># NBVAL_CHECK_OUTPUT
molecules = []
for mol_file in mol_files:
    molecule = Chem.MolFromPDBFile(str(mol_file), removeHs=False)
    if molecule is None:
        print(mol_file, "Could not be read")
    else:
        Chem.SanitizeMol(molecule)
        print(Chem.MolToSmiles(molecule))
        molecules.append(molecule)
print(f"Number of molecules: {len(molecules)}")
Draw.MolsToGridImage(
    molecules,
    molsPerRow=4,
    legends=[f"From PDB ID: {pdb_id}" for pdb_id in pdb_ids],
)<jupyter_output><empty_output><jupyter_text>Upon inspection of the structures we notice, that we encountered a problem here: When loading ligands from a PDB file, RDKit does not assign bonds in structures like aromatic rings correctly. We use the RDKit function `AssignBondOrdersFromTemplate`, which assigns bonds to a molecule based on a reference molecule, e.g. in our case based on the SMILES pattern of the molecule. Check for further information: [RDKit discussion on *Aromaticity of non-protein molecules in PDB not detected*](https://github.com/rdkit/rdkit/issues/1031) and [RDKit documentation on `AssignBondOrdersFromTemplate`](http://rdkit.org/docs/source/rdkit.Chem.AllChem.htmlAssignBondOrdersFromTemplate).<jupyter_code># Load SMILES for PDB ligand structures
ligands = pd.read_csv(HERE / "../T008_query_pdb/data/PDB_top_ligands.csv")

# Get SMILES in the same order as in pdb_ids
ligand_smiles = [
    ligands[ligands["@structureId"] == pdb_id]["smiles"].values[0] for pdb_id in pdb_ids
]

# Generate RDKit Mol object from SMILES
reference_molecules = [Chem.MolFromSmiles(smiles) for smiles in ligand_smiles]

# Assign bond orders to molecules based on SMILES patterns of reference molecules
molecules = [
    AllChem.AssignBondOrdersFromTemplate(reference_molecule, molecule)
    for reference_molecule, molecule in zip(reference_molecules, molecules)
]
Draw.MolsToGridImage(
    molecules,
    molsPerRow=4,
    legends=[f"From PDB ID: {pdb_id}" for pdb_id in pdb_ids],
)<jupyter_output><empty_output><jupyter_text>Now the correct bonds are included.We can also have a look at the molecules in 2D. In order to keep the original coordinates, we make a copy for the 2D representation.<jupyter_code>molecules_2d = []
for molecule in molecules:
    tmp = Chem.Mol(molecule)
    AllChem.Compute2DCoords(tmp)
    molecules_2d.append(tmp)
Draw.MolsToGridImage(
    molecules_2d,
    molsPerRow=4,
    legends=[f"From PDB ID: {pdb_id}" for pdb_id in pdb_ids],
)<jupyter_output><empty_output><jupyter_text>Visualize with `nglview` In the next step we visualize the ligand molecules with nglview using the following function.The molecules will be drawn on top of each other in a single view. This way we are able to visually compare strucural differences between the ligands.<jupyter_code>def show_ligands(molecules):
    """Generate a view of the ligand molecules.

    Parameters
    -----------
    molecules: list of rdkit.Chem.rdchem.Mol

    Returns
    ----------
    nglview.widget.NGLWidget
    """
    view = nv.NGLWidget()
    for molecule in molecules:
        component = view.add_component(molecule)
        time.sleep(0.1)
        component.clear()
        component.add_ball_and_stick(multipleBond=True)
    return view
view = show_ligands(molecules)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Extract pharmacophore featuresAs described above, the aim of this talktorial is to generate a ligand-based ensemble pharmacophore from a set of ligands. First, we need to extract pharmacophore features per ligand.Therefore, we load a feature factory. We use the default feature definitions.See also the [RDKit documentation on chemical features and pharmacophores](https://rdkit.readthedocs.io/en/latest/GettingStartedInPython.htmlchemical-features-and-pharmacophores).<jupyter_code>feature_factory = AllChem.BuildFeatureFactory(str(Path(RDConfig.RDDataDir) / "BaseFeatures.fdef"))<jupyter_output><empty_output><jupyter_text>We take a look at the pharmacophore features that are implemented in RDKit:<jupyter_code>list(feature_factory.GetFeatureDefs().keys())<jupyter_output><empty_output><jupyter_text>As an example, we extract all features from the first molecule.<jupyter_code>example_molecule = molecules[0]
example_molecule
# NBVAL_CHECK_OUTPUT
features = feature_factory.GetFeaturesForMol(example_molecule)
print(f"Number of features found: {len(features)}")<jupyter_output>Number of features found: 17<jupyter_text>The type of a feature is called family in RDKit and can be retrieved with `GetFamily()`.<jupyter_code># NBVAL_CHECK_OUTPUT
features[0].GetFamily()<jupyter_output><empty_output><jupyter_text>First we get the frequency of feature types for our example molecule.<jupyter_code># NBVAL_CHECK_OUTPUT
feature_frequency = collections.Counter([feature.GetFamily() for feature in features])
feature_frequency<jupyter_output><empty_output><jupyter_text>Now apply the functions shown above to all molecules in our ligand set. We calculate the feature type frequency per molecule, then we display it in a DataFrame.<jupyter_code>molecule_feature_frequencies = []
for molecule in molecules:
    features = [feature.GetFamily() for feature in feature_factory.GetFeaturesForMol(molecule)]
    feature_frequency = collections.Counter(features)
    molecule_feature_frequencies.append(feature_frequency)

feature_frequencies_df = (
    pd.DataFrame(
        molecule_feature_frequencies,
        index=[f"Mol{i}" for i, _ in enumerate(molecules, 1)],
    )
    .fillna(0)
    .astype(int)
)
feature_frequencies_df.transpose()<jupyter_output><empty_output><jupyter_text>Furtheron, we concentrate in this talktorial only on the following feature types: hydrogen bond acceptors (acceptors), hydrogen bond donors (donors), and hydrophobic contacts (hydrophobics).We retrieve the feature RDKit objects per feature type and per molecule.<jupyter_code>acceptors = []
donors = []
hydrophobics = []

for molecule in molecules:
    acceptors.append(feature_factory.GetFeaturesForMol(molecule, includeOnly="Acceptor"))
    donors.append(feature_factory.GetFeaturesForMol(molecule, includeOnly="Donor"))
    hydrophobics.append(feature_factory.GetFeaturesForMol(molecule, includeOnly="Hydrophobe"))

features = {
    "donors": donors,
    "acceptors": acceptors,
    "hydrophobics": hydrophobics,
}<jupyter_output><empty_output><jupyter_text>Show the pharmacophore features of all ligandsPharmacophore feature types usually are displayed in defined colors, e.g. usually hydrogen bond donors, hydrogen bond acceptors, and hydrophobic contacts are colored green, red, and yellow, respectively.<jupyter_code>feature_colors = {
    "donors": (0, 0.9, 0),  # Green
    "acceptors": (0.9, 0, 0),  # Red
    "hydrophobics": (1, 0.9, 0),  # Yellow
}
def visualize_features(
    molecules,
    features,
    feature_type="features",
    color="yellow",
    sphere_radius=0.5,
):
    """Generate a view of the molecules highlighting the specified feature type.

    Parameters
    -----------
    molecules: list of rdkit.Chem.rdchem.Mol
        molecules to be visualized
    features: list of tuples of rdkit.Chem.rdMolChemicalFeatures.MolChemicalFeature
        extracted features from molecules chosen to be highlighted
    feature_type: string, optional
        name of the feature to be highlighted
    color: string, optional
        color used to display the highlighted features
    sphere_radius: float, optional
        display size of the highlighted features

    Returns
    ----------
    nglview.widget.NGLWidget
    """
    print(f"Number of {feature_type} in all ligands: {sum([len(i) for i in features])}")
    view = show_ligands(molecules)
    for i, feature_set in enumerate(features, 1):
        for feature in feature_set:
            loc = list(feature.GetPos())
            label = f"{feature_type}_{i}"
            view.shape.add_sphere(loc, color, sphere_radius, label)
    return view<jupyter_output><empty_output><jupyter_text>We use this function to visualize the features for the feature types under consideration.  Hydrogen bond donors<jupyter_code>feature_type = "donors"
view = visualize_features(
    molecules,
    features[feature_type],
    feature_type,
    feature_colors[feature_type],
)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Hydrogen bond acceptors<jupyter_code>feature_type = "acceptors"
view = visualize_features(
    molecules,
    features[feature_type],
    feature_type,
    feature_colors[feature_type],
)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Hydrophobic contacts<jupyter_code>feature_type = "hydrophobics"
view = visualize_features(
    molecules,
    features[feature_type],
    feature_type,
    feature_colors[feature_type],
)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Collect coordinates of features per feature typeWe want to cluster features per feature type. To do this we collect all coordinates of features and store them per feature type.<jupyter_code>features_coord = {
    "donors": [list(item.GetPos()) for sublist in features["donors"] for item in sublist],
    "acceptors": [list(item.GetPos()) for sublist in features["acceptors"] for item in sublist],
    "hydrophobics": [
        list(item.GetPos()) for sublist in features["hydrophobics"] for item in sublist
    ],
}<jupyter_output><empty_output><jupyter_text>Now, we have the positions of e.g. all acceptor features:<jupyter_code>features_coord["acceptors"]<jupyter_output><empty_output><jupyter_text>Generate ensemble pharmacophoresIn order to generate ensemble pharmacophores, we use k-means clustering to cluster features per feature type.  Set static parameters for k-means clustering`kq`: With this paramter, we determine the number of clusters `k` per feature type depending on the number of feature points, i.e. per feature type:`k` = number_of_features / `kq``kq` should be selected so that `k` (feature clusters) is for all clusters at least 1 and not larger than 4-5 clustersWe set our 'kq' used to determine k in k-means as follows<jupyter_code>kq = 7<jupyter_output><empty_output><jupyter_text>Set static parameters for cluster selection`min_cluster_size`: We only want to retain clusters that potentially contain features from most molecules in our ligand ensemble. Therefore, we set this variable to 75% of the number of molecules in our ligand ensemble.`top_cluster_number`: With this parameter, we select only the largest clusters.<jupyter_code>min_cluster_size = int(len(molecules) * 0.75)
top_cluster_number = 4<jupyter_output><empty_output><jupyter_text>Define k-means clustering and cluster selection functionsWe define a function that calculates the centers of clusters, which are derived from k-means clustering.<jupyter_code>def clustering(feature_coord, kq):
    """
    Compute the k-means clustering of input feature coordinates.

    Parameters
    -----------
    feature_coord: list
        feature coordinates extracted from the features of a set of molecules
    kq: int
        Quotient used when calculating number of features k

    Returns
    ----------
    sklearn.cluster._kmeans.KMeans
    """
    # Define parameter k as feature number divided by "k quotient"
    k = math.ceil(len(feature_coord) / kq)
    # Tailor-made adaption of k for hydrophobics in for the example in this talktorial
    k = 2 if k == 1 else k
    print(f"Clustering: \nVariable k in k-means: {k} of {len(feature_coord)} points\n")

    # Initialize k-means and compute clustering
    k_means = cluster.KMeans(n_clusters=k)
    k_means.fit(feature_coord)

    return k_means<jupyter_output><empty_output><jupyter_text>We define a function that sorts the clusters by size and outputs a list of indices of the largest clusters.<jupyter_code>def get_clusters(k_means, min_cluster_size, top_cluster_number):
    """
    Select clusters from a k-means clustering and returns the indices of the selected clusters.
    The selection is based on cluster size. The largest clusters are chosen.

    Parameters
    -----------
    k_means: sklearn.cluster._kmeans.KMeans
        k-means clustering of feature coordinates
    min_cluster_size: int
        minimum size needed for a cluster to be viable for selection
    top_cluster_number: int
        number of clusters to be selected

    Returns
    ----------
    list of int
        indices of selected clusters
    """

    # Sort clusters by size and only show largest
    feature_labels = k_means.labels_

    feature_labels_count = Counter(feature_labels)

    feature_labels_count = sorted(
        feature_labels_count.items(), key=operator.itemgetter(1), reverse=True
    )
    print(f"Clusters sorted by size: \n{feature_labels_count}\n")

    # Get number of the largest clusters, which are larger then the threshold (selected clusters)
    cluster_indices_sel = []

    for cluster_index, cluster_size in feature_labels_count:
        if cluster_size >= min_cluster_size and top_cluster_number > 0:
            cluster_indices_sel.append(cluster_index)
            top_cluster_number -= 1

    print(f"Cluster indices of selected clusters: \n{cluster_indices_sel}\n")

    return cluster_indices_sel<jupyter_output><empty_output><jupyter_text>Cluster featuresFor each feature type, we perform the k-means clustering with our defined `clustering` function.<jupyter_code># NBVAL_CHECK_OUTPUT
k_means = {
    "donors": clustering(features_coord["donors"], kq),
    "acceptors": clustering(features_coord["acceptors"], kq),
    "hydrophobics": clustering(features_coord["hydrophobics"], kq),
}<jupyter_output>Clustering: 
Variable k in k-means: 2 of 11 points

Clustering: 
Variable k in k-means: 4 of 23 points

Clustering: 
Variable k in k-means: 2 of 7 points<jupyter_text>Select relevant clustersFor each feature type, we select relevant clusters with our defined `get_clusters` function.<jupyter_code>print("Hydrogen bond donors\n")
cluster_indices_sel_donors = get_clusters(k_means["donors"], min_cluster_size, top_cluster_number)
print("Hydrogen bond acceptors\n")
cluster_indices_sel_acceptors = get_clusters(
    k_means["acceptors"], min_cluster_size, top_cluster_number
)
print("Hydrophobic contacts\n")
cluster_indices_sel_hydrophobic = get_clusters(
    k_means["hydrophobics"], min_cluster_size, top_cluster_number
)
cluster_indices_sel = {
    "donors": cluster_indices_sel_donors,
    "acceptors": cluster_indices_sel_acceptors,
    "hydrophobics": cluster_indices_sel_hydrophobic,
}<jupyter_output><empty_output><jupyter_text>Get selected cluster coordinates<jupyter_code>def get_selected_cluster_center_coords(k_means, cluster_indices_sel):
    """Retrieve cluster center coordinates for selected clusters.
    The selected clusters are defined by their indices.

    Parameters
    -----------
    k_means: sklearn.cluster._kmeans.KMeans
        k-means clustering of feature coordinates
    cluster_indices_sel: list of int
        indices of the chosen clusters

    Returns
    ----------
    two-dimensional list of floats
        cluster center coordinates
    """
    cluster_centers = k_means.cluster_centers_

    # Cast to list and then to pandas Series (for element selection by indices)
    cluster_centers = pd.Series(cluster_centers.tolist())

    # Select cluster centers by indices of selected clusters
    cluster_centers_sel = cluster_centers[cluster_indices_sel]

    return list(cluster_centers_sel)
cluster_centers_sel = {
    "donors": get_selected_cluster_center_coords(k_means["donors"], cluster_indices_sel["donors"]),
    "acceptors": get_selected_cluster_center_coords(
        k_means["acceptors"], cluster_indices_sel["acceptors"]
    ),
    "hydrophobics": get_selected_cluster_center_coords(
        k_means["hydrophobics"], cluster_indices_sel["hydrophobics"]
    ),
}
cluster_centers_sel["acceptors"]<jupyter_output><empty_output><jupyter_text>Show clustersPer feature type, we visualize cluster centers alongside with all molecules and all feature points.<jupyter_code>def visualize_clusters(
    molecules,
    features,
    clusters,
    feature_type="features",
    color="yellow",
    feature_radius=0.5,
    cluster_radius=1,
):
    """Visualize the cluster centers.

    Parameters
    -----------
     molecules: list of rdkit.Chem.rdchem.Mol
        molecules to be visualized
    features: list of tuples of rdkit.Chem.rdMolChemicalFeatures.MolChemicalFeature
        extracted features from molecules chosen to be highlighted
    clusters: Two-dimensional list of floats
        cluster center coordinates
    feature_type: string, optional
        feature type being visualized
    color: string, optional
        color used for rendering the features and cluster centers
    feature_radius: float, optional
        rendering size of the features
    cluster_radius: float, optional
        rendering size of the cluster centers

    Returns
    ----------
    nglview.widget.NGLWidget
    """
    view = visualize_features(
        molecules,
        features,
        feature_type,
        color=color,
        sphere_radius=feature_radius,
    )
    for i, center in enumerate(clusters, 1):
        view.shape.add_sphere(list(center), color, cluster_radius, f"cluster_{feature_type}_{i}")
    return view<jupyter_output><empty_output><jupyter_text>Hydrogen bond donors<jupyter_code>feature_type = "donors"
view = visualize_clusters(
    molecules,
    features[feature_type],
    cluster_centers_sel[feature_type],
    feature_type,
    feature_colors[feature_type],
)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Hydrogen bond acceptor<jupyter_code>feature_type = "acceptors"
view = visualize_clusters(
    molecules,
    features[feature_type],
    cluster_centers_sel[feature_type],
    feature_type,
    feature_colors[feature_type],
)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Hydrophobic contacts<jupyter_code>feature_type = "hydrophobics"
view = visualize_clusters(
    molecules,
    features[feature_type],
    cluster_centers_sel[feature_type],
    feature_type,
    feature_colors[feature_type],
)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output><jupyter_text>Show ensemble pharmacophoreIn this last step, we combine the clustered pharmacophoric features (i.e. hydrogen bond donors and acceptors as well as hydrophobic contacts), to one ensemble pharmacophore, representing the pharmacophoric properties of the four selected ligands.<jupyter_code>view = show_ligands(molecules)
# Load clusters
for feature_type in cluster_indices_sel.keys():
    centers = cluster_centers_sel[feature_type]
    for i, loc in enumerate(centers):
        sphere_radius = 1
        feature_color = feature_colors[feature_type]
        label = f"{feature_type}_c{i}"
        view.shape.add_sphere(loc, feature_color, sphere_radius, label)
view
view.render_image(trim=True, factor=2, transparent=True);
view._display_image()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08344948357814232, 'fraction_numerical': 0.013223258639791288, 'mean_word_length': 4.953617021276596, 'pattern_counts': {'":': 18, '<': 135, '<?xml version=': 0, '>': 137, 'https://': 13, 'lorem ipsum': 0, 'www.': 2, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '525587', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 8767, 'n_tokens_neox': 7906, 'n_words': 2589, 'path': 'teachopencadd/talktorials/T009_compound_ensemble_pharmacophores/talktorial.ipynb', 'repo_name': 'PabloSoto1995/teachopencadd', 'size': '28049'}