<jupyter_start><jupyter_text>Copyright 2020 NVIDIA Corporation SPDX-License-Identifier: Apache-2.0<jupyter_code>from dask.distributed import Client

client = Client()
# client = Client("tcp://127.0.0.1:36005")
client
import cudf, cuml
from cuml import KMeans, UMAP
from dask import delayed
import dask.dataframe as dd
import dask.bag as db
import pandas as pd
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem
from chembl_webresource_client.new_client import new_client
from chembl_webresource_client.utils import utils
from PIL import Image
import numpy as np
import os, sys, wget, gzip
import hashlib
from ftplib import FTP
import matplotlib.pyplot as plt
import progressbar<jupyter_output><empty_output><jupyter_text>ChEMBL SearchIllustrating the efficient lazy search and retrieval API for ChEMBL<jupyter_code>molecule = new_client.molecule

# filter for approved drugs, filter for molfile present
approved_drugs = molecule.filter(max_phase=4).filter(molecule_structures__isnull=False)

# filter for all small molecules, filter for molfile present
small_molecules = molecule.filter(molecule_type='Small molecule').filter(molecule_structures__isnull=False)

print(len(approved_drugs))
print(len(small_molecules))<jupyter_output>3322
1914279<jupyter_text>RDKit MolFromMolBlockWe've filtered above to ensure that molecule structure is present.  From the molfile, we can get an RDKit mol.<jupyter_code>mol = approved_drugs[0]['molecule_structures']['molfile']
m = Chem.MolFromMolBlock(mol)
m<jupyter_output><empty_output><jupyter_text>Fingerprint to ArrayThe 4 in ECFP4 corresponds to the diameter of the atom environments considered, while the Morgan fingerprints take a radius parameter. So a Morgan fingerprint with radius=2 is roughly equivalent to ECFP4 and FCFP4.<jupyter_code># Get a ECFP4/Morgan fingerprint
mol = approved_drugs[0]['molecule_structures']['molfile']
m = Chem.MolFromMolBlock(mol)
fp = AllChem.GetMorganFingerprintAsBitVect(m, radius=2)
ar = np.array(fp)
ar<jupyter_output><empty_output><jupyter_text>List of BitVectorsLet's just get a list of a few fingerprints so we can test putting into a DataFrame<jupyter_code>fingerprints = []
max_mols = 20

for chembl_mol in approved_drugs[:max_mols]:
    # Get a ECFP4/Morgan fingerprint
    mol = chembl_mol['molecule_structures']['molfile']
    m = Chem.MolFromMolBlock(mol)
    fp = AllChem.GetMorganFingerprintAsBitVect(m, radius=2)
    ar = np.array(fp)
    fingerprints.append(ar)

fingerprints = np.asarray(fingerprints, dtype=np.float32)
len(fingerprints)<jupyter_output><empty_output><jupyter_text>cuDF GPU DataFrame of FingerprintsThis is a test of a function to convert a list of fingerprints to a GPU DataFrame<jupyter_code>def np2cudf(df):
    # convert numpy array to cuDF dataframe
    df = pd.DataFrame({'fea%d'%i:df[:,i] for i in range(df.shape[1])})
    pdf = cudf.DataFrame()
    for c,column in enumerate(df):
        pdf[str(c)] = df[column]
    return pdf

gdf = np2cudf(fingerprints)
gdf.head()<jupyter_output><empty_output><jupyter_text>Clustering with K-Means<jupyter_code>n_clusters = 3

kmeans_float = KMeans(n_clusters=n_clusters)
kmeans_float.fit(gdf)
print("labels:")
print(kmeans_float.labels_)
print("cluster_centers:")
print(kmeans_float.cluster_centers_)<jupyter_output>labels:
0     0
1     0
2     2
3     0
4     0
5     0
6     0
7     2
8     2
9     0
10    0
11    0
12    0
13    0
14    0
15    0
16    1
17    0
18    0
19    0
dtype: int32
cluster_centers:
   0      1     2     3     4     5     6     7     8     9     ...    2038  \
0   0.0  0.125   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  ...  0.0625   
1   0.0  0.000   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  ...  0.0000   
2   0.0  0.000   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  ...  0.0000   

   2039  2040  2041  2042  2043  2044      2045  2046  2047  
0   0.0   0.0   0.0   0.0   0.0   0.0  0.000000   0.0   0.0  
1   0.0   0.0   0.0   0.0   1.0   0.0  0.000000   0.0   0.0  
2   0.0   0.0   0.0   0.0   0.0   0.0  0.333333   0.0   0.0  

[3 rows x 2048 columns]<jupyter_text>Dask Computation of FingerprintsThe first time this is run, the ChEMBL molecules will be downloaded so it will take a while.  The data is automatically cached.  The next time it is run, this should only take about a minute.<jupyter_code>def GetMorganFP(chembl_mol, radius=2):
    molfile = chembl_mol['molecule_structures']['molfile']
    m = Chem.MolFromMolBlock(molfile)
    fp = AllChem.GetMorganFingerprintAsBitVect(m, radius=radius)
    ar = np.array(fp)
    return ar


def ToNpArray(fingerprints):
    fingerprints = np.asarray(fingerprints, dtype=np.float32)
    return fingerprints
fingerprints = []
max_mols = 100

for x in approved_drugs[:max_mols]:
    y = delayed(GetMorganFP)(x)
    fingerprints.append(y)
fingerprints = delayed(ToNpArray)(fingerprints)

results = fingerprints.compute()
print(results.shape)<jupyter_output>(100, 2048)<jupyter_text>Saving SMILES StringsThis is inefficient.  The [ftp site](ftp://ftp.ebi.ac.uk/pub/databases/chembl/ChEMBLdb/releases/chembl_26) is a better way to download ChEMBL data.  There, 2048 radius-2 Morgan FPS, SMILES, InChI etc can be downloaded directly.<jupyter_code>I_really_want_to_do_this = False

if I_really_want_to_do_this:
    
    fn_aprvd_mols = 'aprvd_molecule_smiles.txt'
    fn_small_mols = 'small_molecule_smiles.txt'
    with open(fn_aprvd_mols,'w') as fp:
        for compound in approved_drugs:
            fp.write(compound['molecule_structures']['canonical_smiles']+'\n')
    with open(fn_small_mols,'w') as fp:
        for compound in small_molecules:
            fp.write(compound['molecule_structures']['canonical_smiles']+'\n')<jupyter_output><empty_output><jupyter_text>Download SMILES from FTP<jupyter_code>chemreps_url = 'ftp://ftp.ebi.ac.uk/pub/databases/chembl/ChEMBLdb/releases/chembl_26/chembl_26_chemreps.txt.gz'
chemreps_local_path = '/data/chembl_26_chemreps.txt.gz'
chemreps_sha256 = '0585b113959592453c2e1bb6f63f2fc9d5dd34be8f96a3a3b3f80e78d5dbe1bd'
chemreps_exists_and_is_good = False

while not chemreps_exists_and_is_good:
    if os.path.exists(chemreps_local_path):
        with open(chemreps_local_path, 'rb') as file:
            local_sha256 = hashlib.sha256(file.read()).hexdigest()
        if chemreps_sha256==local_sha256:
            chemreps_exists_and_is_good = True
            print('chembl chemreps file found locally, SHA256 matches')
    if not chemreps_exists_and_is_good:
        print('downloading chembl chemreps file...')
        wget.download(chemreps_url, chemreps_local_path)<jupyter_output>chembl chemreps file found locally, SHA256 matches<jupyter_text>Morgan fingerprints from saved SMILES strings<jupyter_code>def MorganFromSmiles(smiles, radius=2, nBits=512):
    m = Chem.MolFromSmiles(smiles)
    fp = AllChem.GetMorganFingerprintAsBitVect(m, radius=radius, nBits=nBits)
    ar = np.array(fp)
    return ar


def ToNpArray(fingerprints):
    fingerprints = np.asarray(fingerprints, dtype=np.float32)
    return fingerprints
smiles_list = []
chemblID_list = []
count=1
max=10000

with gzip.open(chemreps_local_path, 'rb') as fp:
    fp.__next__()
    for i,line in enumerate(fp):
        fields = line.split()
        chemblID_list.append(fields[0].decode("utf-8"))
        smiles_list.append(fields[1].decode("utf-8"))
        count+=1
        if count>max: 
            break
%%time
# create dask bag
results = db.from_sequence(smiles_list).map(MorganFromSmiles).compute()
m = Chem.MolFromSmiles(smiles_list[1])
%%time
# take list of len n_mols of np.arrays of len nBits and stack
np_array_fingerprints = np.stack(results).astype(np.float32)

# take np.array shape (n_mols, nBits) for GPU DataFrame
gdf = np2cudf(np_array_fingerprints)

n_clusters = 7

kmeans_float = KMeans(n_clusters=n_clusters)
kmeans_float.fit(gdf)
gdf.head()
print(smiles_list[0])
umap = UMAP()
Xt = umap.fit_transform(gdf)
X = Xt.as_matrix()
# plt.scatter(X[:,0], X[:,1], c=kmeans_float.labels_, s=.2)
# Add data back to gdf for visualization
gdf.add_column('x', Xt[0].to_array())
gdf.add_column('y', Xt[1].to_array())

gdf.add_column('cluster', kmeans_float.labels_.to_array())

# Add a id column, tobe used for point selection in interactive UI
gdf.head(10)
import chemvisualize
import importlib
importlib.reload(chemvisualize)

v = chemvisualize.ChemVisualization(gdf.copy(), n_clusters, chemblID_list.copy())

v.start('0.0.0.0')<jupyter_output>Dash is running on http://0.0.0.0:5000/

Dash is running on http://0.0.0.0:5000/

 * Serving Flask app "chemvisualize" (lazy loading)
 * Environment: production
[31m   WARNING: This is a development server. Do not use it in a production deployment.[0m
[2m   Use a production WSGI server instead.[0m
 * Debug mode: off
{'fraction_non_alphanumeric': 0.1007309273640932, 'fraction_numerical': 0.05150753768844221, 'mean_word_length': 4.136070381231671, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 4, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '544972', 'license': 'Apache-2.0', 'n_tokens_mistral': 3523, 'n_tokens_neox': 3233, 'n_words': 768, 'path': 'demo.ipynb', 'repo_name': 'ohadmo/cheminformatics', 'size': '8756'}