# ---
# 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="Se3M6z2uQIGU" colab_type="text"
# # Setup
#

# + [markdown] id="DIMxNnjfQU2G" colab_type="text"
# Google’s updated their Colab instances to use TensorFlow 2.X by default. The last stable version of the Signaturizers is still on TensorFlow 1.X.

# + id="Z605c7pvP3tx" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="905b68b9-3f58-4677-c21d-05fb5af8c2db"
# %tensorflow_version 1.x

# + [markdown] id="1PA55FW9M4x1" colab_type="text"
# ## Install the _signaturizer_ package

# + id="xpnKeewTM9yT" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 732} outputId="91f45a7e-3237-4c4f-f1d0-09abf5711728"
# !pip install signaturizer

# + [markdown] id="B4OHQNcBMvp_" colab_type="text"
# ## Install RDKit

# + id="ivgaJ2f5M1x7" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 1000} outputId="a2b35479-ce11-47ae-b13c-e8f39ace895d"
# !wget -c https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
# !chmod +x Miniconda3-latest-Linux-x86_64.sh
# !bash ./Miniconda3-latest-Linux-x86_64.sh -b -f -p /usr/local
# !conda install -q -y -c conda-forge rdkit

# + [markdown] id="vZmnpVrAjWSS" colab_type="text"
# Add RDKit to system PATH

# + id="0QBnPcl_Q9i5" colab_type="code" colab={}
import sys
sys.path.append('/usr/local/lib/python3.7/site-packages/')

# + [markdown] id="T5MDJ-6diK9f" colab_type="text"
# # Basic usage

# + [markdown] id="PH17HwDMiGSJ" colab_type="text"
# ## Generating Bioactivity Signatures

# + [markdown] id="DtKnPHu5NILM" colab_type="text"
# Let's start loading the signaturizer for B1 space which generates the Mode of Action bioactivity descriptor. The MoA is represented by the exact molecular target/s or the physiological outcomes (if the target has not yet been elucidated or is subject to debate). 
#
# The MoA of "Aspirin", for example, would be "inhibitor of COX-1" while its broader mode of action would be "blocks prostaglandin synthesis".

# + id="aSXE4FRsOthC" colab_type="code" colab={}
from signaturizer import Signaturizer
sign = Signaturizer('B1')

# + [markdown] id="57CPy_YBNS_Z" colab_type="text"
# To generate signatures for a compound we need it's structure represented in [SMILES format](https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system).
#
# Let's take for example the smiles for 3 common pharmacy drugs:
# - Aspirin: CC(=O)OC1=CC=CC=C1C(=O)O
# - Lovastatin: CCC(C)C(=O)OC1CC(C=C2C1C(C(C=C2)C)CCC3CC(CC(=O)O3)O)C
# - Simvastatin: CC(C)C1=C(C(=C(N1CCC(CC(CC(=O)O)O)O)C2=CC=C(C=C2)F)C3=CC=CC=C3)C(=O)NC4=CC=CC=C4

# + id="Io5wNV5LuWSV" colab_type="code" colab={}
aspirin = 'CC(=O)OC1=CC=CC=C1C(=O)O'
lovastatin ='CCC(C)C(=O)OC1CC(C=C2C1C(C(C=C2)C)CCC3CC(CC(=O)O3)O)C'
simvastatin = 'CC(C)C1=C(C(=C(N1CCC(CC(CC(=O)O)O)O)C2=CC=C(C=C2)F)C3=CC=CC=C3)C(=O)NC4=CC=CC=C4'

# + id="MS9SfFfYNavH" colab_type="code" colab={}
smiles = [aspirin, lovastatin, simvastatin]
results = sign.predict(smiles)

# + id="Wp0nC6j6iWHb" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="0b12ac46-52c0-478f-e16b-407734cbd157"
print(results.signature.shape)

# + id="7z9lfmdqPvNp" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 1000} outputId="1de1b27c-ea59-4bb6-e718-48e753d04b05"
print(results.signature)

# + [markdown] id="elO_JGQLWL3O" colab_type="text"
# We can measure the distance between signatures to highlight that the two statins have indeed a similar MoA, which is different from the aspirin.

# + id="SWoQuLE8WMEM" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 69} outputId="2bb10d02-6a34-407d-9d17-f6cd189a7250"
from scipy.spatial.distance import euclidean
print(euclidean(results.signature[0], results.signature[1]))
print(euclidean(results.signature[0], results.signature[2]))
print(euclidean(results.signature[1], results.signature[2]))

# + [markdown] id="wj6b_Xl4imiK" colab_type="text"
# ## Generating multiple Signatures in one batch

# + [markdown] id="ov7mmmaYkgiu" colab_type="text"
# We might want to describe out molecules by multiple point of view using different descriptors. The signaturizer package allow easy and modular bioactivity space cambinations.
#
# Signaturizers are loaded as TF hub modules so they are cached locally, in this way model download (takes time) is happening only once

# + id="7wYXLKpMid9u" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="d3ae470f-092e-47cf-9bb0-0b799a24d4c9"
sign = Signaturizer(['B1', 'E4', 'A5', 'D3'])
results = sign.predict(smiles)
print(results.signature.shape)

# + [markdown] id="n9lWWYyMizlF" colab_type="text"
# ## Generating Global Bioactivity Signatures

# + [markdown] id="8Koa5I3PN09_" colab_type="text"
# A _global_ signature consists of stacking all the bioactivity spaces described by the Chemical Checker, this is convinient when we want a rich _general purpose_ descriptor.

# + id="NsLm7iKPisYF" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="04bd771c-9f06-4886-9237-0d0926e24a28"
sign = Signaturizer('GLOBAL')
results = sign.predict(smiles)
print(results.signature.shape)

# + [markdown] id="NCbAJPnYjuv6" colab_type="text"
# ## Saving results to h5

# + [markdown] id="uEsM2iHbk7iH" colab_type="text"
# When generating signatures for a lot of molecules it is better to save them to disk and avoid consuming all system memory.

# + id="NQtJidPIi38B" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="1ab97787-51b6-4d6b-aae0-3b7fd9aa5d09"
# !rm destination.h5
results = sign.predict(smiles, 'destination.h5')

# + [markdown] id="AC6vY-CmxlyI" colab_type="text"
# We can read the results anytime using the [h5py package](https://www.h5py.org/).

# + id="kUi2c5TajZVB" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="38b146e0-e381-46a7-f6ea-53b3ed6e762e"
import h5py
with h5py.File('destination.h5', 'r') as fh:
  signatures = fh['signature'][:]
print(signatures.shape)

# + [markdown] id="OJ6BIyrXlYx6" colab_type="text"
# ## Applicability predictions

# + [markdown] id="_y6i4vmIlfkJ" colab_type="text"
# The `applicability` score helps identifying molecules for which the signature prediction is extremely accurate

# + id="uSKHZhj_j5AG" colab_type="code" colab={}
results = sign.predict(smiles, applicability=True)

# + [markdown] id="J3WvU9Y8nWBG" colab_type="text"
# Each bioactivity signature comes with an applicability score, since we are using the global signature we will have 25 applicability score for each molecule.

# + id="2xlYd59enHWq" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="daee9527-dcc5-47dd-a149-54d6f75895b5"
print(results.applicability.shape)

# + [markdown] id="kVXJ1RDUoAf8" colab_type="text"
# NB: In `results.dataset` we store the bioactivity dataset code.

# + id="zBrA7oVmn6cQ" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 52} outputId="df053f6d-252a-4335-e1ec-44f4e2b3c297"
print(results.dataset)

# + [markdown] id="KvJa2Migp0tc" colab_type="text"
# We can highligth most reliable predictions (applicability is tuned as measure of correlation, i.e. how much we think the signature is related to real values)

# + id="Wr4S21_foqoQ" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 72} outputId="05df2037-f7ac-454e-a50f-a3e68dd6b31a"
# %matplotlib inline
import matplotlib
import seaborn as sns
import matplotlib.pyplot as plt

# + id="iLvF9WaToznw" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 88} outputId="324c3565-e7c3-484a-dbf1-c4f8ce8aa8eb"
plt.imshow(results.applicability)
plt.gca().set_xticks(range(len(results.dataset)));
plt.gca().set_xticklabels(results.dataset);

# + [markdown] id="lQXDRC09p2dk" colab_type="text"
# # Usage example

# + [markdown] id="oeM4qUz11LlG" colab_type="text"
# ## Data source
# In this notebook we will analize a representative subset of the [FooDB](https://foodb.ca/) compound library exploiting the Chemical Checker signatures generated _via_ our Signaturizers. Let's first download and inspect the  data:

# + id="KCUPgpB1p7Pb" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 208} outputId="1e44cd48-e5b9-424c-a808-85586f36534e"
# !wget http://gitlabsbnb.irbbarcelona.org/packages/signaturizer/raw/master/notebook/foodb_repr.csv

# + id="3zZwIfVs1Dfv" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 424} outputId="db2cc453-5113-46ea-9298-661c0a079808"
import pandas as pd
df = pd.read_csv('./foodb_repr.csv')
df

# + [markdown] id="Kw7UMnSq1WNy" colab_type="text"
# We have different column in this file, the most important for our task is the **SMILES** columns, which is the input required to generate signatures.

# + id="FmA6lA7j1SIL" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 121} outputId="50a5c22b-8b9e-4a54-9837-2911680911f8"
df.nunique()

# + [markdown] id="6kQEAen-1cGO" colab_type="text"
# Comparing the total number of rows (47954) to the number of unique entries (6325) we realize that molecules are repeated, and indeed different food can contain repeated compounds. For this reason we'll use a non redundant set of SMILES.

# + id="QiXculT_1Ysx" colab_type="code" colab={}
smiles = df.SMILES.unique()

# + [markdown] id="lV-O6j6W1icV" colab_type="text"
# ## Signatures
# Signatures are numerical representation of molecules. Having molecules represented as vector is needed if we want to do any quantitative analysis of our library (e.g. comparison, clustering, visualization ...). The advantage of using a signaturizer with respect to standard chemical representation (e.g. morgan fingerprint) is that we can leverage on the wealth of public bioactivity data that goes beyond chemical information.
#
# In the Chemical Checker we defined 25 bioactivity spaces which capture multiple aspect of a molecule. The stacked version of the signatures for these 25 spaces is what we call the _global_ signature, i.e. a general purpose entry point for our analisys.

# + id="7p9R91uG1er7" colab_type="code" colab={}
# load the global bioactivity signaturizer
from signaturizer import Signaturizer
sign = Signaturizer('GLOBAL', verbose=False)

# + [markdown] id="TwkcsgWp1rH0" colab_type="text"
# For the sake of time we can just consider a random subsampling of these molecules.

# + id="N1b5vRbd1poj" colab_type="code" colab={}
import numpy as np
np.random.seed(42)
subsample = np.random.choice(smiles, 3000, replace=False)
result = sign.predict(subsample)

# + id="-tW3NoYK1t9r" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="fc13947c-0f26-40a6-9b06-134eddb8a073"
result.signature.shape

# + [markdown] id="vb8_CN1C11bB" colab_type="text"
# ## Visualization
# Only imagination limits us in the way we can visualize the results generated.

# + id="ibYq88w91v4W" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 72} outputId="d138d732-0d78-4c0c-efdf-46f0913127ef"
# %matplotlib inline
import matplotlib
import seaborn as sns
from matplotlib import gridspec, colors, colorbar, patches, lines
import matplotlib.pyplot as plt

# + [markdown] id="zJRH0Q6f17cU" colab_type="text"
# We can plot the overall distribution of values for all generated signatures, and appreciate their bimodal distribution.

# + id="KHPBRfiq14rI" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 286} outputId="bd6ab4f1-ae29-4edf-e649-5c5d5f16afed"
sns.distplot(result.signature.ravel())

# + [markdown] id="VM6cuYsV2Brp" colab_type="text"
# Or we can take a peek at the raw matrix.

# + id="1KUGTeLS19-k" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 173} outputId="b3e54a2a-3bbe-456f-860d-52326d0a12b4"
plt.imshow(result.signature)

# + [markdown] id="MdAJjHDY2KxJ" colab_type="text"
# The signatures we generated have 3200 dimensions but looking at the raw data does not inform us on how similar/different are our molecules.
#
# In Cheminformatics we often need a graphical representation of the chemical space covered by a set of molecules.  In this representation, similar molecules will be close and molecules that are different will be far apart.
#
# We can use scikit-learn to create a lower bi-dimensional representation which will enable such visualization.

# + id="VhPVjLAa2EP6" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 416} outputId="12a35524-8d53-4552-8f1b-c1b33b843cf1"
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
pca = PCA(n_components=128)
pca_red = pca.fit_transform(result.signature)
print(pca.explained_variance_ratio_)
print(np.sum(pca.explained_variance_ratio_))
projection = TSNE(n_components=2).fit_transform(pca_red)

# + id="PubHSuwAEaqk" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="bbe6f8db-152f-4ed2-d670-1c9f6e7ae2d2"
print(projection.shape)

# + id="BhOw47IB2WQ6" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 285} outputId="98ebf5a3-86aa-4640-eb04-a8e73a5768f6"
plt.scatter(projection[:,0],projection[:,1], color='grey', alpha=.4)

# + [markdown] id="YdFc6ir02c6k" colab_type="text"
# Now we can appreciate the diversity of our molecules and also that they are 'organized' with some form of structure. We can annotate molecules with the food of origin to highlight this structure.
#
# For example let's color all molecules coming from animal origin (i.e. meat):

# + id="rmGJXIqY2Y9n" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 285} outputId="a55554ac-514b-4377-c3d1-dbc88e24106b"
cat_smiles = df[df.food_group == 'Animal foods'].SMILES.unique()
cat_mask = np.isin(subsample, cat_smiles)
plt.scatter(projection[:,0],projection[:,1], color='lightgrey', alpha=.2)
plt.scatter(projection[:,0][cat_mask],projection[:,1][cat_mask], color='red', alpha=0.2)

# + [markdown] id="nvUr_fOL2jVB" colab_type="text"
# We clearly see that molecules from this food source are localized in two clusers of our bioactivity space.
#
# If we go for compounds of vegetable orgin we can see, as expectable, that the two set of molecules do not overlap.

# + id="LrDzqFA02e7I" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 282} outputId="c4f3027b-6dbf-4f6c-d374-704272e53798"
cat_smiles = df[df.food_group == 'Vegetables'].SMILES.unique()
cat_mask = np.isin(subsample, cat_smiles)
plt.scatter(projection[:,0],projection[:,1], color='lightgrey', alpha=.2)
plt.scatter(projection[:,0][cat_mask],projection[:,1][cat_mask], color='green', alpha=0.2)

# + id="O63aNP4D2mkt" colab_type="code" colab={}


StarcoderdataJupyterScriptsDedupFiltered