<jupyter_start><jupyter_text>Compound clustering Aim of this notebookCompound filtering: Cluster set of compounds in order to extract a diverse set of compounds for downstream tasks, e.g. docking.<jupyter_code>%load_ext autoreload
%autoreload 2
from pathlib import Path

import numpy as np
import pandas as pd
from rdkit import Chem, DataStructs
from rdkit.Chem import PandasTools, rdFingerprintGenerator, Draw
from rdkit.ML.Cluster import Butina
DATA_PATH = Path('.') / '..' / '..' / 'data'
MOLECULES_PATH = DATA_PATH / 'focused_library_similar_proteins' / 'focused_library.csv'<jupyter_output><empty_output><jupyter_text>Load molecule library<jupyter_code>dataset = pd.read_csv(MOLECULES_PATH)
dataset.shape
PandasTools.AddMoleculeColumnToFrame(dataset, smilesCol='canonical_smiles')<jupyter_output><empty_output><jupyter_text>Generate fingerprints<jupyter_code>rdkit_gen = rdFingerprintGenerator.GetRDKitFPGenerator(maxPath=5)
dataset['fingerprint'] = dataset.apply(lambda x: rdkit_gen.GetFingerprint(x.ROMol), axis=1)<jupyter_output><empty_output><jupyter_text>Cluster molecules (Butina algorithm) Cluster molecules with different distance cutoffs - how many clusters come out?<jupyter_code>def cluster_molecules(fingerprints, cutoff=0.6):
    """
    Cluster fingerprints using the Butina algorithm.
    
    Parameters
    ----------
    fingerprints : list of rdkit.DataStructs.cDataStructs.ExplicitBitVect
        List of fingerprints.
    cutoff : float
        Distance cutoff Butina clustering.
        
    Returns
    -------
    list of tuple of int
        List of clusters, whereby each cluster is described by its cluster member IDs.
    """
    
    # Calculate Tanimoto distance matrix
    distance_matrix = _tanimoto_distance_matrix(fingerprints)
    
    # Now cluster the data with the implemented Butina algorithm:
    clusters = Butina.ClusterData(
        distance_matrix,
        len(fingerprints),
        cutoff,
        isDistData=True
    )
    
    # Sort clusters by size
    clusters = sorted(clusters, key=len, reverse=True)
    
    return clusters

def _tanimoto_distance_matrix(fingerprints):
    """
    Calculate distance matrix for list of fingerprints.
    
    Parameters
    ----------
    fingerprints : list of rdkit.DataStructs.cDataStructs.ExplicitBitVect
        List of fingerprints.
        
    Returns
    -------
    list of floats
        Distance matrix (a triangular similarity matrix in the form of a list)
    """
    
    fingerprints = list(fingerprints)
    distance_matrix = []
    
    for i in range(1,len(fingerprints)):
        similarities = DataStructs.BulkTanimotoSimilarity(fingerprints[i], fingerprints[:i])
        
        # Since we need a distance matrix, calculate 1-x for every element in similarity matrix
        distance_matrix.extend([1-x for x in similarities])
    
    return distance_matrix
# Perform clustering with different distance cutoffs

cluster_by_cutoff = {}

for cutoff in np.arange(0, 1, 0.1):
    
    cutoff = round(cutoff, 1)
    
    clusters = cluster_molecules(dataset.fingerprint, cutoff)
    cluster_by_cutoff[cutoff] = clusters<jupyter_output><empty_output><jupyter_text>Get clustering stats (per distance cutoff) Use these statistics to decide on distance cutoff.<jupyter_code>cluster_stats = {
    cutoff: [
        cutoff,  
        sum([len(cluster) > 100 for cluster in clusters]),
        sum([len(cluster) > 25 for cluster in clusters]),
        sum([len(cluster) > 10 for cluster in clusters]),
        sum([len(cluster) > 5 for cluster in clusters]),
        sum([len(cluster) == 1 for cluster in clusters]),
        len(clusters)
    ] for cutoff, clusters in cluster_by_cutoff.items()
}
cluster_stats = pd.DataFrame(
    cluster_stats.values(), 
    columns=['Distance cutoff', '# Clusters > 100', '# Clusters > 25', '# Clusters > 10', '# Clusters > 5', '# Clusters = 1', '# Clusters']
)
cluster_stats
cluster_stats.set_index('Distance cutoff').transpose().plot(figsize=(9,6))<jupyter_output><empty_output><jupyter_text>Draw clusters This function allows you to draw each cluster of a clustering based on a selected distance cutoff.<jupyter_code>def draw_clusters(mols, clusters, max_mols=1000):
    """
    Draw molecules sorted by descending cluster size.
    
    Parameters
    ----------
    mols : list of rdkit.Chem.rdchem.Mol
        Molecules.
    clusters : list of tuple of int
        List of clusters, whereby each cluster is described by its cluster member IDs.
    max_mols : int
        Number of molecules shown per cluster ID (default 1000).
    
    Returns
    -------
    dict of PIL.PngImagePlugin.PngImageFile
        Dictionary of molecules for each cluster (legend = "cluster ID | molecule ID").
    """
    
    img_dict = {}
    
    for cluster_id, cluster in enumerate(clusters):
        
        clustered_mols = []
        legends = []
        
        clustered_mols.extend([mols[mol_id] for mol_id in cluster])
        legends.extend([f'{cluster_id+1} | {mol_id}' for mol_id in cluster])
        
        img = Draw.MolsToGridImage(clustered_mols, legends=legends,  molsPerRow=6, maxMols=max_mols)
        img_dict[cluster_id+1] = img
        
    return img_dict
img_dict = draw_clusters(dataset.ROMol, cluster_by_cutoff[0.6], max_mols=50)
# Show cluster 1 (largest cluster) - by default "only" the first 1000 molecules are shown per cluster
img_dict[1]<jupyter_output><empty_output><jupyter_text>Get cluster centroids First molecule in each cluster is cluster centroid - thus, extract first molecule per cluster!<jupyter_code>cluster_centroids_by_cutoff = {cutoff: [cluster[0] for cluster in clusters] for cutoff, clusters in cluster_by_cutoff.items()} 
def draw_cluster_centroids(mols, cluster_centroids):
    """
    Draw fragments sorted by descending cluster size.
    
    Parameters
    ----------
    mols : list of rdkit.Chem.rdchem.Mol
        List of top x fragments.
    cluster_centroids : list of int
        List of cluster centroids (their molecule IDs).
    
    Returns
    -------
    PIL.PngImagePlugin.PngImageFile
        Image of cluster centroids (legend = "cluster ID | molecule ID").
    """
    
    cluster_centroids_mols = []
    legends = []
    
    for cluster_id, molecule_id in enumerate(cluster_centroids):
        
        cluster_centroids_mols.append(mols[molecule_id])
        legends.append(f'{cluster_id+1} | {molecule_id}')
        
    img = Draw.MolsToGridImage(cluster_centroids_mols, legends=legends,  molsPerRow=6)
    
    return img
draw_cluster_centroids(dataset.ROMol, cluster_centroids_by_cutoff[0.6])<jupyter_output><empty_output><jupyter_text>Save diverse compounds set to file<jupyter_code>dataset.iloc[cluster_centroids_by_cutoff[0.6]].to_csv(MOLECULES_PATH.parent / 'focused_library_diverse.csv')
print(f'Number of moleules in diverse set: {len(cluster_centroids_by_cutoff[0.6])}')<jupyter_output>Number of moleules in divers set: 109
{'fraction_non_alphanumeric': 0.0937905200979686, 'fraction_numerical': 0.009508716323296355, 'mean_word_length': 3.4932038834951458, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '523445', 'license': 'MIT', 'n_tokens_mistral': 2233, 'n_tokens_neox': 2087, 'n_words': 550, 'path': 'code/compound_clustering/compound_clustering.ipynb', 'repo_name': 'volkamerlab/covid19-SBapproach', 'size': '6941'}