# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/filizolalab/CADD2021/blob/main/Tutorial_Example_BSR3101_Davide.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="pVJlMNH_X67c"
# # (Example Tutorial for BRS3101) Molecule Featurization and Datasets
#
# **What will we learn**
#
# In this tutorial we introduce some key concepts regarding molecular datasets for applicaiton of machine learning to drug design and molecular modeling.
# 1. Different featurization techniques (description of molecular structures for ML)
# 2. Datasets, iterators, train/validation/test splitting
#
# **What will we use**
#
# We will use `deepchem` and `rdKit` within the [Google Colaboratory](https://colab.research.google.com/) framework. 
# * Google Colaboratory is a cloud-based platform to run python code without the need to setup your local machine with a python environment. It allows you to write and execute Python in your browser, with zero configuration required, Free access to GPUs, Easy sharing. The code you develop and run in colab will run in any python environment.
#
# * [`deepchem`](https://deepchem.io/) is a pyhton library that implements tools to create high quality, open source tools for drug discovery, materials science, quantum chemistry, and biology. 
#
# * [`rdKit`](https://www.rdkit.org/) is a library of chemoinformatics 

# + [markdown] id="A8dJDVuLblLN"
# ## Install `deepchem` and `rdkit`

# + colab={"base_uri": "https://localhost:8080/"} id="BV21K0DGYGF3" outputId="02de5e5a-d649-4e6d-def7-9c31a73e1ad4"
# !pip install --pre deepchem
# !pip install --pre rdkit-pypi 
# !pip install --pre PubChemPy

# + colab={"base_uri": "https://localhost:8080/", "height": 35} id="Y5djSAYnYNjl" outputId="cfc4b32c-bff9-41ee-a290-fe468d3f5ca1"
import deepchem as dc
dc.__version__

# + [markdown] id="89dF5Fy_dNkr"
# ## Featurizers

# + [markdown] id="0OiU2hAtiBiM"
# Molecules can be represented in different ways.
# * fingerprints (e.g. Circular Fingerprints, )
# * molecular properties 
# * Graph convolution models [Duvenaud et al.](https://arxiv.org/abs/1509.09292)
# * Tokenizers
#
# [list of featurizers](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6917215/)
#

# + [markdown] id="CLAKI0TfWhy1"
# ### Fingerprints

# + id="yEoQQgrldPeD"
from rdkit import Chem
smiles = ['C1=CC=CC=C1']

# Example 1: (size = 2048, radius = 4)
featurizer = dc.feat.CircularFingerprint(size=2048, radius=4)
features = featurizer.featurize(smiles)


# + colab={"base_uri": "https://localhost:8080/"} id="Isinq0e4dTO3" outputId="9c389c8d-69e9-48a2-ee9f-10eb8d8663cb"
features.shape

# + colab={"base_uri": "https://localhost:8080/"} id="MBxHSehuWLUK" outputId="451e912b-4b8d-4fac-a49f-2a6c1835e33f"
for i,fp in enumerate(features[0]): 
  if fp==1: print(i)

# + [markdown] id="MS3PO5nIWll7"
# ### Molecular properties

# + id="B03cYIWHWniy"
smiles = ['CC(=O)OC1=CC=CC=C1C(=O)O']
featurizer = dc.feat.RDKitDescriptors()
features = featurizer.featurize(smiles)

# + colab={"base_uri": "https://localhost:8080/"} id="ZJQG5hQ6XAz0" outputId="b971bf23-1a88-4788-9004-36f091476e92"
features

# + [markdown] id="dU3au4NykMaV"
# ## Datasets
# https://deepchem.readthedocs.io/en/latest/api_reference/moleculenet.html
#

# + [markdown] id="fOeieHOObfqg"
# ## Load an example dataset (Delaney dataset of molecular solubilities)
#
#

# + [markdown] id="wamatJyyhc5x"
# The Delaney (ESOL) dataset a regression dataset containing structures and water solubility data for 1128 compounds. The dataset is widely used to validate machine learning models on estimating solubility directly from molecular structures (as encoded in SMILES strings).
#
# The raw data csv file contains columns below:
# * “Compound ID” - Name of the compound
# * “smiles” - SMILES representation of the molecular structure
# * “measured log solubility in mols per litre” - Log-scale water solubility of the compound, used as label
#
# ### the Dataset structure
# Now let's consider the contents of the Dataset. Every Dataset stores a list of samples. Very roughly speaking, a sample is a single data point. In this case, each sample is a molecule. In other datasets a sample might correspond to an experimental assay, a cell line, an image, or many other things. For every sample the dataset stores the following information.
#
# * The features, referred to as `X`. This is the input that should be fed into a model to represent the sample.
# * The labels, referred to as `y`. This is the desired output from the model. During training, it tries to make the model's output for each sample as close as possible to y.
# * The weights, referred to as `w`. This can be used to indicate that some data values are more important than others. In later tutorials we will see examples of how this is useful.
# * An `ID`, which is a unique identifier for the sample. This can be anything as long as it is unique. Sometimes it is just an integer index, but in this dataset the ID is a SMILES string describing the molecule.
#
#
# The final piece of information listed in the output is `task_names`. Some datasets contain multiple pieces of information for each sample. For example, if a sample represents a molecule, the dataset might record the results of several different experiments on that molecule. This dataset has only a single task: "measured log solubility in mols per litre". 
#

# + id="rYyXsl85YQut"
tasks, datasets, transformers = dc.molnet.load_delaney(featurizer='GraphConv')

# + id="RckBYgdXYcRC"
train_dataset, valid_dataset, test_dataset = datasets

# + colab={"base_uri": "https://localhost:8080/"} id="tVtNbU7PYwJ8" outputId="9a0b72f2-d6ca-4e97-a2a4-a8b00fe60537"
train_dataset

# + [markdown] id="AUD5gMnJ062I"
# ## Prepare and featurize your own dataset
#
#
# Here we will read a file containing smiles and labels and prepare a dataset for machine learning.
#
# * read from local (see [here](https://colab.research.google.com/notebooks/io.ipynb) for colab file IO)
#

# + colab={"resources": {"http://localhost:8080/nbextensions/google.colab/files.js": {"data": "<KEY>", "ok": true, "headers": [["content-type", "application/javascript"]], "status": 200, "status_text": "OK"}}, "base_uri": "https://localhost:8080/", "height": 74} id="55FmnVlW1AJD" outputId="601df1fd-5f5e-4e7a-d8fa-73f9c2f56244"
from google.colab import files
import csv
uploaded = files.upload()

# + id="IJ4t4QxT1Yd2"

lines = [line.decode("utf-8") for line in uploaded['test1.csv'].splitlines()]
reader = csv.reader(lines)
molecules = list(reader)


# + id="FFxnWLI1140i"
IDs = [molecule[0] for molecule in molecules]
smiles = [molecule[1] for molecule in molecules]
label = [molecule[2] for molecule in molecules]

features = featurizer.featurize(smiles)
dataset = dc.data.NumpyDataset(X=features, y=label, ids=IDs)


# + id="BHDkDXZt_UA6"
splitter = dc.splits.RandomSplitter()

train_dataset, valid_dataset, test_dataset = splitter.train_valid_test_split(dataset)

# + id="P0ggduKsCH4f"
print(f'total {len(dataset)}')

print(f'train {len(train_dataset)}')
print(f'valid {len(valid_dataset)}')
print(f'test {len(test_dataset)}')



# + id="tRQlZgukEm1O"
train_dataset.X.shape

# + [markdown] id="zHo_thGghuSB"
# ## Train a deep-learning model to predict 

# + id="4ky5fHosDhQi"
model = dc.models.MultitaskClassifier(n_tasks=1, n_features=2048, layer_sizes=[1000])

# + id="OtZEzMxbEZin"

model.fit(train_dataset, nb_epoch=10)
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)



# + colab={"base_uri": "https://localhost:8080/"} id="N-fp9sfTEdTd" outputId="5caa9f2b-004c-4f3d-c420-20bcd5eaba74"
metric

# + colab={"base_uri": "https://localhost:8080/"} id="NxkJswnCFOK7" outputId="94951e77-16f7-4ba2-ad41-704b193ba375"
print('training set score:', model.evaluate(train_dataset, metric))
print('test set score:', model.evaluate(test_dataset, metric))

# + [markdown] id="1yvUGIWmGrq8"
# Let us plot the results
#
# [how to plot in Colab](https://colab.research.google.com/github/jakevdp/PythonDataScienceHandbook/blob/master/notebooks/04.02-Simple-Scatter-Plots.ipynb#scrollTo=RmV1JUqGGH0X)
#
#

# + id="oZA1UVNYFWRt"
# %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np

xy = model.predict(test_dataset).reshape(32,2)

# + colab={"base_uri": "https://localhost:8080/", "height": 262} id="843tcd50Fn6E" outputId="1220707b-a356-48d6-8879-7b8e6c9186b4"
plt.plot(xy[:,0], xy[:,1], 'o', color='black');


# + colab={"base_uri": "https://localhost:8080/"} id="6IxTfC6UGP6e" outputId="83e6b188-2dba-4f89-eef9-98e42423fff1"
xy

# + [markdown] id="wToNIL5emi9W"
# # Cluster ligands by similarity in rdKit
#
# [examples of basic `rdkit` use](https://www.rdkit.org/docs/GettingStartedInPython.html)
#

# + id="BFUNxDY_GcyS"
from rdkit import Chem
from rdkit.Chem import AllChem


# + id="TafBgmgjmqID"
# Define clustering setup
def cluster_fingerprints(fps,cutoff=0.2):
    from rdkit import DataStructs
    from rdkit.ML.Cluster import Butina

    # first generate the distance matrix:
    dists = []
    nfps = len(fps)
    for i in range(1,nfps):
        sims = DataStructs.BulkTanimotoSimilarity(fps[i],fps[:i])
        dists.extend([1-x for x in sims])

    # now cluster the data:
    cs = Butina.ClusterData(dists,nfps,cutoff,isDistData=True)
    return cs



# + id="dfoN19_hoIi4"
featurizer = dc.feat.CircularFingerprint(size=2048, radius=4)
tasks, datasets, transformers = dc.molnet.load_delaney(featurizer=featurizer)

# + colab={"base_uri": "https://localhost:8080/"} id="ae0CL8tgowao" outputId="5f34dff3-022a-4eff-e496-3f47a8428624"
datasets[0].X

# + id="F2hRms69nUrS"
ms = [Chem.MolFromSmiles(x) for x in smiles]

# + id="xLPKgimToOFS"
fps = [AllChem.GetMorganFingerprintAsBitVect(x,4,2048) for x in ms]

# + id="PyLfx7OHpFSW"
clusters = cluster_fingerprints(fps,cutoff=0.4)

# + colab={"base_uri": "https://localhost:8080/"} id="FuokegX7pFx4" outputId="0ff99c3e-d3f3-4174-e516-a771d5e5e005"
#show one of the clusters
print(clusters[0])



# + colab={"base_uri": "https://localhost:8080/", "height": 1000} id="7aAgnQ6Npd3R" outputId="67180713-9b1b-4673-e32f-45fac1c44246"
#now display structures from one of the clusters
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole


mols=[ms[i] for i in clusters[2]]
Draw.MolsToGridImage(mols)

# + id="MigWx3ICpkJP"
from rdkit.Chem import Descriptors

# + colab={"base_uri": "https://localhost:8080/"} id="HE7GyXEVtL7Q" outputId="090e89d8-c613-4054-a904-22ab90657f71"

[Descriptors.MolLogP(mol) for mol in ms]

# + id="lMv5YQ1btOKj"


StarcoderdataJupyterScriptsDedupFiltered