<jupyter_start><jupyter_text>Talktorial 9 Ligand-based pharmacophores Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin <NAME>, <NAME> and <NAME> **Note**: Please run this notebook cell by cell. Running all cells in one is possible also, however, a few PyMol images might not turn out as intended.  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 Theory* Pharmacophore modeling  * Structure- and ligand-based pharmacophore modeling* Virtual screening with pharmacophores* Clustering: k means Practical* Get pre-aligned ligands from previous talktorial* Start PyMoL* Show ligands with PyMol* 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://iupac.org/publications/pac/70/5/1129/))* 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))* PyMol integration with RDKit ([rdkit.Chem.PyMol documentation](http://rdkit.org/docs/source/rdkit.Chem.PyMol.html))* <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://iupac.org/publications/pac/70/5/1129/))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)). Figure 1: Structure-based pharmacophore representing protein-ligand interactions (figure by <NAME>). Virtual screening with pharmacophoresAs described earlier in **talktorial 4**, virtual screening (VS) describes the screening of a query (e.g. here in **talktorial 9** a pharmacophore model or in **talktorial 4** a query compound) against a large library of compounds, in order to identify those 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 used is 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.([K means wikipedia](https://de.wikipedia.org/wiki/K-Means-Algorithmus))  Practical<jupyter_code>import os, glob

# RDKit
from rdkit import RDConfig, Chem, Geometry, DistanceGeometry
from rdkit.Chem import ChemicalFeatures, rdDistGeom, Draw, rdMolTransforms, AllChem
from rdkit.Chem.Draw import IPythonConsole, DrawingOptions
from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib
from rdkit.Numerics import rdAlignment
IPythonConsole.ipython_useSVG=True

# PyMOL related
from rdkit.Chem import PyMol
from pymol import *
import time # Needed for waiting a second
from PIL import Image # For export the image to disk

import collections
import pandas as pd
import math

from sklearn import datasets, cluster
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from collections import Counter # For handling the labels
import operator<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 = []
for file in glob.glob("../data/T8/*_lig.pdb"):
    mol_files.append(file)
mol_files
pdb_ids = [i.split("/")[-1].split("_")[0] for i in mol_files]
pdb_ids<jupyter_output><empty_output><jupyter_text>Second, we read all ligands from these PDB files using RDKit.<jupyter_code>mols = []
for mol_file in mol_files:
    mol = Chem.MolFromPDBFile(mol_file, removeHs=False)
    if mol is None:
        print(mol_file, 'could not be read')
    else:
        Chem.SanitizeMol(mol)
        print(Chem.MolToSmiles(mol))
        mols.append(mol)
rangeMols = range(1, len(mols)+1)
print('Number of molecules: ', len(mols))
Draw.MolsToGridImage(mols, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])<jupyter_output><empty_output><jupyter_text>We encounter a problem here: When loading ligands from a PDB file, RDKit does not assign e.g. aromatic rings to the ligand. 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
ligs = pd.read_csv("../data/T8/PDB_top_ligands.csv", sep="\t")

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

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

# Assign bond orders to molecules (mols) based on SMILES patterns (refmols)
mols = [AllChem.AssignBondOrdersFromTemplate(refmol, mol) for refmol, mol in zip(refmols, mols)]
Draw.MolsToGridImage(mols, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])<jupyter_output><empty_output><jupyter_text>We can also have a look at the molecules in 2D (we copy the molecules for this example to keep the original coordinates).<jupyter_code>mols_2D = [] 
for mol in mols:
    tmp=Chem.Mol(mol)
    AllChem.Compute2DCoords(tmp)
    mols_2D.append(tmp)
Draw.MolsToGridImage(mols_2D, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])<jupyter_output><empty_output><jupyter_text>Start PyMolWe start PyMol within the terminal.<jupyter_code># Open PyMol in shell
os.popen('pymol -R')<jupyter_output><empty_output><jupyter_text>Note: If no separate PyMol window opens at this stage, something may be wrong with your PyMol installation and needs to be fixed first! You can also try to manually start PyMol by typing `pymol -R` in your shell. We need to wait until PyMol is launched completely. Then we can link PyMol to the Jupyter notebook via RDKit: `objPMV = PyMol.MolViewer()`. Check general functionalities of PyMol integration with RDKit on ([rdkit.Chem.PyMol documentation](http://rdkit.org/docs/source/rdkit.Chem.PyMol.html)).<jupyter_code># Error handling: wait until PyMol is loaded

nrTry = 0  # Number of current attempt
ttw = 10  # Time to wait in seconds

while nrTry < ttw:  # Try until PyMol is loaded and the object can be saved
    nrTry += 1
    try:  
        objPMV = PyMol.MolViewer()  # Save the PyMol object
        break  # Stop the loop when PyMol is loaded        
    except ConnectionRefusedError:  # Exception handling if PyMol is not loaded yet
        time.sleep(1) # Wait...
    
if nrTry == ttw: # After ttw trys: print error message
    print("Error: PyMol did not start correctly.\n" +
          "Try again and/or check if PyMol is installed completely.")<jupyter_output><empty_output><jupyter_text>This is the RDKit PyMol object, which we will be using in the following in order to control PyMol from within the Jupyter notebook:<jupyter_code>objPMV<jupyter_output><empty_output><jupyter_text>The two most important commands of the RDKit PyMol integration are: * `objPMV.ShowMol(object, name, showOnly)` to load the object in PyMol* `objPMV.GetPNG(h=height)` to show the figure in the jupyter notebook  Show ligands with PyMolWe show all ligands (pre-aligned in previous talktorial) with PyMol. We load each molecule individually into PyMol and set a unique PyMol name (m1, m2, ...).<jupyter_code>for mol, i in zip(mols, rangeMols):
    objPMV.ShowMol(mol, name='m%d'%i, showOnly=False)
    i += 1<jupyter_output><empty_output><jupyter_text>All ligands should be visible now in the PyMol window. We display a picture of the PyMol window in this talktorial via retrieving a *png* picture.<jupyter_code>objPMV.GetPNG(h=300)<jupyter_output><empty_output><jupyter_text>Next, we define a function from the steps described above. We add some styling PyMol commands. You can pass PyMol commands from RDKit to PyMol using:`objPMV.server.do("any_pymol_command")`<jupyter_code>def visualize_ligands(objPMV, molecules):
    '''
    This function shows all input molecules in PyMol (within the Jupyter notebook).
    '''
    
    # Initialize PyMol in order to remove all previous objects
    objPMV.server.do("reinitialize")
    
    # Load ligands
    rangeMols = range(1, len(molecules)+1)
    for mol, i in zip(molecules, rangeMols):
        objPMV.ShowMol(mol, name='mol_%d'%i, showOnly=False)
        toStickCmd='cmd.show("sticks","mol_'+str(i)+'")'
        objPMV.server.do(toStickCmd)
        i += 1
    
    # Turn camera
    objPMV.server.do("turn x, -40")
    
    # Set background to white
    objPMV.server.do("bg_color white")
    # Zoom in on ligands
    objPMV.server.do("zoom")
    # Turn on ray tracing for better image quality
    objPMV.server.do("ray 1800, 1000")
    
    # Export as PNG file
    outputPNG = objPMV.GetPNG(w=1800, h=1000)
    outputPNG.save("../data/T9/ligands.png", )
    
    # Display in Jupyter notebook
    return objPMV.GetPNG(h=300)
visualize_ligands(objPMV, mols)<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 (with the default feature definitions).See also [rdkit docu on chemical features and pharmacophores](https://rdkit.readthedocs.io/en/latest/GettingStartedInPython.htmlchemical-features-and-pharmacophores).<jupyter_code>ffact = AllChem.BuildFeatureFactory(os.path.join(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(ffact.GetFeatureDefs().keys())<jupyter_output><empty_output><jupyter_text>As an example, we get all feature for an example molecule.<jupyter_code>m1 = mols[0]
feats = ffact.GetFeaturesForMol(m1)
print('Number of features found:',len(feats))<jupyter_output>Number of features found: 14<jupyter_text>The type (in RDKit called family) of a feature can be retrieved with `GetFamily()`.<jupyter_code>feats[0].GetFamily()<jupyter_output><empty_output><jupyter_text>We get the frequency of features types for our example molecule.<jupyter_code>feats_freq = collections.Counter([x.GetFamily() for x in feats])
feats_freq<jupyter_output><empty_output><jupyter_text>We apply the functions shown above to all molecules in our ligand set. We display the frequency of feature types per molecule as DataFrame.<jupyter_code># Get feature type frequency per molecule
mols_feats_freq = []
for i in mols:
    feats = [x.GetFamily() for x in ffact.GetFeaturesForMol(i)]
    feats_freq = collections.Counter(feats)
    mols_feats_freq.append(feats_freq)

# Show data as DataFrame
p = pd.DataFrame(mols_feats_freq, index=["m"+str(i) for i in range(1, len(mols)+1)]).fillna(0).astype(int)
p.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 i in mols:
    acceptors.append(ffact.GetFeaturesForMol(i, includeOnly='Acceptor'))
    donors.append(ffact.GetFeaturesForMol(i, includeOnly='Donor'))
    hydrophobics.append(ffact.GetFeaturesForMol(i, 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<jupyter_output><empty_output><jupyter_text>RDKit's PyMol integration allows us to draw spheres (representing pharmacophore features) in PyMoL using the following command:`objPMV.server.sphere(loc, sphereRad, colors[i], label, 1)`<jupyter_code>def visualize_features(objPMV, molecules, feature_type, features, feature_color):
    '''
    This function displays all input molecules and all input features as spheres in PyMOL. 
    A png picture from the PyMOL window is loaded into the Jupyter Notebook and saved as file to disc.
    At the end, the PyMOL session is cleaned from all objects.
    '''
    
    # Initialize PyMol in order to remove all previous objects
    objPMV.server.do("reinitialize")
    
    print("Number of " + feature_type + " in all ligands: " + str(sum([len(i) for i in features])))
    
    # Load ligands
    rangeMols = range(1, len(molecules)+1)
    for mol, i in zip(molecules, rangeMols):
        objPMV.ShowMol(mol, name='mol_%d'%i, showOnly=False)
        toStickCmd='cmd.show("sticks","mol_'+str(i)+'")'
        objPMV.server.do(toStickCmd)
        i += 1
    
    # Load features
    for i in range(len(features)):
        
        for feature in features[i]:
            loc = list(feature.GetPos())
            sphere_radius = 0.5
            label = feature_type + '_%d'%(i+1)
            objPMV.server.sphere(loc, sphere_radius, feature_color, label, 1) # show the sphere (pharmacophore feature)
    
    # Turn camera
    objPMV.server.do("turn x, -40")
    
    # Set background to white
    objPMV.server.do("bg_color white")
    # Zoom in on ligands
    objPMV.server.do("zoom")
    # Turn on ray tracing for better image quality
    objPMV.server.do("ray 1800, 1000")
    
    # Export as PNG file
    outputPNG = objPMV.GetPNG(w=1800, h=1000)
    outputPNG.save("../data/T9/ligands_features_"+feature_type+".png", )
    
    # Display in Jupyter notebook
    return objPMV.GetPNG(h=300)<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"
visualize_features(objPMV, mols, feature_type, features[feature_type], feature_colors[feature_type])<jupyter_output>Number of donors in all ligands: 11<jupyter_text>Hydrogen bond acceptors<jupyter_code>feature_type = "acceptors"
visualize_features(objPMV, mols, feature_type, features[feature_type], feature_colors[feature_type])<jupyter_output>Number of acceptors in all ligands: 23<jupyter_text>Hydrophobic contacts<jupyter_code>feature_type = "hydrophobics"
visualize_features(objPMV, mols, feature_type, features[feature_type], feature_colors[feature_type])<jupyter_output>Number of hydrophobics in all ligands: 7<jupyter_text>Collect coordinates of features per feature typeSince we want to cluster features (per feature type), we now collect all coordinates of features (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`<jupyter_code># k quotient (kq) used to determine k in k-means: k = number of feature points / kq
# kq should be selected so that k (feature clusters) is for all clusters at least 1 and not larger than 4-5 clusters
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 cluster.<jupyter_code># Threshold for clustering: number = percentage of threshold value
min_cluster_size = int(len(mols) * 0.75)

# Show only top features
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, kd):
    '''
    This functions computes the k-means clustering of input feature coordinates.
    '''
    
    # Define parameter k as feature number divided by "k quotient"
    k = math.ceil(len(feature_coord) / kq)
    k = 2 if k == 1 else k  # Tailor-made adaption of k for hydrophobics in for the example in this talktorial
    print('Clustering: \nVariable k in k-means: %d of %d points\n'%(k, len(feature_coord)))
    
    # Initialize of k-means
    k_means = cluster.KMeans(n_clusters=k)
    
    # Compute the k-means clustering
    k_means.fit(feature_coord)
    
    # Return the clusters
    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):
    '''
    This function retrieves information on a input k-means clustering:
    * gets cluster label for each feature
    * counts cluster sizes and sort cluster indices by cluster size
    * selects clusters based on size
    * returns selected cluster indices
    '''
    
    # Sort clusters by size and only show largest
    feature_labels = k_means.labels_
    print('Cluster labels for all features: \n%s\n'% feature_labels)

    feature_labels_count = Counter(feature_labels)
    print('Cluster label counter: \n%s\n'% feature_labels_count)

    feature_labels_count = sorted(feature_labels_count.items(), 
                                  key=operator.itemgetter(1), 
                                  reverse=True)
    print('Sorted cluster label counters: \n%s\n'% feature_labels_count)

    # 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:  # feature_labels_count = list of (cluster_index, cluster_size)
        if cluster_size >= min_cluster_size and top_cluster_number > 0:
            cluster_indices_sel.append(cluster_index)
            top_cluster_number -= 1
            
    print('Cluster indices of selected clusters: \n%s\n'% cluster_indices_sel)
    
    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>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_don = get_clusters(k_means["donors"], min_cluster_size, top_cluster_number)
print("Hydrogen bond acceptors\n")
cluster_indices_sel_acc = get_clusters(k_means["acceptors"], min_cluster_size, top_cluster_number)
print("Hydrophobic contacts\n")
cluster_indices_sel_h = get_clusters(k_means["hydrophobics"], min_cluster_size, top_cluster_number)
cluster_indices_sel = {"donors": cluster_indices_sel_don, 
                       "acceptors": cluster_indices_sel_acc, 
                       "hydrophobics": cluster_indices_sel_h}<jupyter_output><empty_output><jupyter_text>Get selected cluster coordinates<jupyter_code>def get_selected_cluster_center_coords(k_means, cluster_indices_sel, feature_type):
    '''
    This function retrieves cluster center coordinates for selected clusters (by their indices).
    '''
    
    # Get cluster centers for a certain feature type
    cluster_centers = k_means[feature_type].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[feature_type]]
    
    # Cast to list and return
    return list(cluster_centers_sel)
cluster_centers_sel = {"donors": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "donors"),
                       "acceptors": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "acceptors"),
                       "hydrophobics": get_selected_cluster_center_coords(k_means, 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(objPMV, molecules, feature_type, features, cluster_centers_sel, feature_color):
    '''
    This function displays 
    * all input molecules, 
    * all input features as spheres, and 
    * the resulting cluster centers in PyMoL. 
    A png picture from the PyMoL window is loaded into the Jupyter notebook and saved as file to disc.
    '''
    
    # Initialize PyMol in order to remove all previous objects
    objPMV.server.do("reinitialize")
    
    print("Number of " + feature_type + " in all ligands: " + str(sum([len(i) for i in features])))
    
    # Load ligands
    rangeMols = range(1, len(molecules)+1)
    for mol, i in zip(molecules, rangeMols):
        objPMV.ShowMol(mol, name='mol_%d'%i, showOnly=False)
        toStickCmd='cmd.show("sticks","mol_'+str(i)+'")'
        objPMV.server.do(toStickCmd)
        i += 1
    
    # Load features
    for i in range(len(features)):
        
        for feature in features[i]:
            loc = list(feature.GetPos())
            sphere_radius = 0.5
            label = feature_type + '_%d'%(i+1)
            objPMV.server.sphere(loc, sphere_radius, feature_color, label, 1) # show the sphere (pharmacophore feature)
    
    # Load clusters
    for i in range(len(cluster_centers_sel)):
        loc = cluster_centers_sel[i]
        sphere_radius = 1
        label = feature_type + '_c%d'%(i+1)
        objPMV.server.sphere(loc, sphere_radius, feature_color, label, 1)
    
    # Turn camera
    objPMV.server.do("turn x, -40")
    
    # Set PyMol styling
    objPMV.server.do("bg_color white")
    objPMV.server.do("zoom")
    objPMV.server.do("ray 1800, 1000")
    
    # Export as PNG file
    outputPNG = objPMV.GetPNG(w=1800, h=1000)
    outputPNG.save("../data/T9/ligands_features_clusters_"+feature_type+".png")
    
    # Display in Jupyter notebook
    return objPMV.GetPNG(h=300)<jupyter_output><empty_output><jupyter_text>Hydrogen bond donors<jupyter_code>feature_type = "donors"
visualize_clusters(objPMV, mols, feature_type, features[feature_type], 
                   cluster_centers_sel[feature_type], 
                   feature_colors[feature_type])<jupyter_output>Number of donors in all ligands: 11<jupyter_text>Hydrogen bond acceptor<jupyter_code>feature_type = "acceptors"
visualize_clusters(objPMV, mols, feature_type, features[feature_type], 
                   cluster_centers_sel[feature_type],  
                   feature_colors[feature_type])<jupyter_output>Number of acceptors in all ligands: 23<jupyter_text>Hydrophobic contacts<jupyter_code>feature_type = "hydrophobics"
visualize_clusters(objPMV, mols, feature_type, features[feature_type], 
                   cluster_centers_sel[feature_type], 
                   feature_colors[feature_type])<jupyter_output>Number of hydrophobics in all ligands: 7<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># Initialize PyMol in order to remove all previous objects
objPMV.server.do("reinitialize")

# Load ligands
rangeMols = range(1, len(mols)+1)
for mol, i in zip(mols, rangeMols):
    objPMV.ShowMol(mol, name='mol_%d'%i, showOnly=False)
    toStickCmd='cmd.show("sticks","mol_'+str(i)+'")'
    objPMV.server.do(toStickCmd)
    i += 1

# Load clusters
for feature_type in cluster_indices_sel.keys():
    centers = cluster_centers_sel[feature_type]
    for i in range(len(centers)):
        loc = centers[i]
        sphere_radius = 1
        feature_color = feature_colors[feature_type]
        label = feature_type + '_c%d'%(i+1)
        objPMV.server.sphere(loc, sphere_radius, feature_color, label, 1)

# Turn camera
objPMV.server.do("turn x, -40")

# Set PyMol styling
objPMV.server.do("bg_color white")
objPMV.server.do("zoom")
objPMV.server.do("ray 1800, 1000")

# Export as PNG file
outputPNG = objPMV.GetPNG(w=1800, h=1000)
outputPNG.save("../data/T9/ligands_ensemble_ph4.png")

# Display in Jupyter notebook
objPMV.GetPNG(h=300)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup