<jupyter_start><jupyter_text>Calculating Atomic Contributions for Molecules Based on Graph Convolutional  QSAR Model. A tutorial.  Why do the DeepChem Tutorial?  You will learn   to find atoms of molecules important for the activity, and visualize them. We will learn on two tasks: classification and regression.  Contributions are also known as "attributions" , "coloration" etc. in literature. Atomic contributions are  defined as the difference in activity of the whole molecule and fragment remaining after atom removal. This is one of model interpretation methods [1], analogous to Similarity maps [2] in QSAR domain, or occlusion methods in other fields (image classification etc). <NAME>, <NAME>. Institute of Molecular and Translational Medicine, Palacky University, Olomouc, Czech Republic.<jupyter_code>import pandas as pd
import deepchem as dc
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw, PyMol, rdFMCS
from rdkit.Chem.Draw import IPythonConsole
from rdkit import rdBase
from deepchem import metrics
from IPython.display import Image, display
from rdkit.Chem.Draw import SimilarityMaps
import tensorflow as tf<jupyter_output><empty_output><jupyter_text>1. First, let's build a classification QSAR model for Blood-brain barrier permeability BBB permeability is the ability of compounds to  enter  central nervous system. Here we use dataset of relatively small  compounds which are transported by diffusion without any  carriers.  The property is defined as log10 (concentration in brain / concentration in blood). Compounds with positive  value (and 0) are labeled active, others inactive.  After modelling we will identify atoms favorable and unfavorable for diffusion.<jupyter_code># 1. Create dataset
DATASET_FILE ='atomic_contributions_tutorial_data/logBB.sdf'
# create RDKit mol objects, we will need them later
mols = [m for m in Chem.SDMolSupplier(DATASET_FILE) if m is not None ]
len(mols)
loader = dc.data.SDFLoader(tasks=["logBB_class"], 

                           featurizer=dc.feat.ConvMolFeaturizer(),
                           sanitize = True)

dataset = loader.create_dataset(DATASET_FILE, shard_size=2000)


dataset.ids.shape
# 2. Build a model
np.random.seed(2020)
tf.random.set_seed(2020)
m = dc.models.GraphConvModel(1, mode="classification", nb_epoch=10,batch_normalize=False, restore=False, batch_size = 100)
m.fit(dataset)

# 3. Test set performance
TEST_DATASET_FILE = 'atomic_contributions_tutorial_data/logBB_test_.sdf'
loader = dc.data.SDFLoader(tasks=["p_np"], sanitize= True,
                           featurizer=dc.feat.ConvMolFeaturizer())

test_dataset = loader.create_dataset(TEST_DATASET_FILE, shard_size=2000)


pred =  m.predict(test_dataset)
pred = np.argmax(np.squeeze(pred),axis= 1)
ba = metrics.balanced_accuracy_score(y_true=test_dataset.y, y_pred=pred)
# balanced accuracy is high enough to procced to model interpretation - i.e. atomic contributions
ba<jupyter_output><empty_output><jupyter_text>2. Now, let's prepare dataset of fragments based on training set (any other unseen data set of interest can be used). These fragments will be used to evaluate contributions of atoms. For each molecule we will generate a list of ConvMol objects: featurizations for single atom-depleted versions of the  molecule (iterate over all heavy atoms in the molecule) All we need is to load again the training set,  but when featurizing it,  set per_atom_fragmentation=True<jupyter_code>loader = dc.data.SDFLoader(tasks=[],# dont need task (moreover, passing the task can lead to inconsitencies in data shapes)
                        featurizer=dc.feat.ConvMolFeaturizer(per_atom_fragmentation=True),
                        sanitize=True) 

frag_dataset = loader.create_dataset(DATASET_FILE, shard_size=5000)
# dataset has the same number of  rows  than original training set, 
# but each row is a fragment list of corresponding molecule. 
# Molecules go in the same order as   original input file. 
# IMPORTANT: fragments order depends on input format: if sdf, then fragment order is the same, as atom order in corresponding mol blocks
# if SMILES (i.e. csv with molecules represented as SMILES), then order is given by RDKit CanonicalRankAtoms

frag_dataset.X.shape

frag_dataset.ids.shape# ids indicate name of the molecule which  fragments ocome from
# transform dataset to flatten fragments lists to make it suitable for prediction 
tr = dc.trans.FlatteningTransformer(frag_dataset)
frag_dataset = tr.transform(frag_dataset)
# dataset has now number of  rows equal total number of fragments in all molecules (i.e. flat)
frag_dataset.X.shape

# the ids are expanded so that each fragment has an id of parent molecule
frag_dataset.ids.shape<jupyter_output><empty_output><jupyter_text>3. Now we will predict the activity for  molecules and for their fragments. Then, for each fragment,  we'll find the activity difference - that is atomic contribution of removed atom, corresponding to that fragment. Note: here, in classification context,  we use probabilistic output of the model as the activity. So, contribution is probability difference, i.e. "how much a given atom increases/decreases the probability of the molecule to be active"<jupyter_code># whole  molecules
pred = np.squeeze(m.predict(dataset))[:, 1] # probabilitiy of class 1
pred = pd.DataFrame(pred, index=dataset.ids, columns=["Molecule"])  # turn to dataframe for convinience
# fragments
pred_frags =np.squeeze(m.predict(frag_dataset))[:, 1]
pred_frags = pd.DataFrame(pred_frags, index=frag_dataset.ids, columns=["Fragment"])<jupyter_output><empty_output><jupyter_text>4.  Let's finally find atomic contributions:<jupyter_code># merge 2 dataframes by molecule names
df = pd.merge(pred_frags, pred, right_index=True, left_index=True)
# find contribs
df['Contrib'] = df["Molecule"] - df["Fragment"]
df<jupyter_output><empty_output><jupyter_text>5. Visualize contributions with RDKit Similarity Maps<jupyter_code>def vis_contribs(mols, df, smi_or_sdf = "sdf"): 
    # input format of file, which was used to create dataset  determines the order of atoms, 
    # so we take it into account for correct mapping!
    maps = []
    for mol  in mols:
        wt = {}
        if smi_or_sdf == "smi":
            for n,atom in enumerate(Chem.rdmolfiles.CanonicalRankAtoms(mol)):
                wt[atom] = df.loc[mol.GetProp("_Name"),"Contrib"][n]
        if smi_or_sdf == "sdf":        
            for n,atom in enumerate(range(mol.GetNumHeavyAtoms())):
                wt[atom] = df.loc[Chem.MolToSmiles(mol),"Contrib"][n]
        maps.append(SimilarityMaps.GetSimilarityMapFromWeights(mol,wt))
    return maps<jupyter_output><empty_output><jupyter_text>Let's look at some pictures:<jupyter_code>np.random.seed(2000)
maps = vis_contribs(np.random.choice(np.array(mols),20), df)<jupyter_output><empty_output><jupyter_text>We can see, that aromatics or aliphatics  have positive impact on membrane permeability, while polar or charged heteroatoms  have negative influence. This is generally consistent with literature data.  6. Let's take  a regression task, aquatic toxicity (towards water organism T. pyriformis). Toxicity is defined as log10 (IGC50) (concentration that inhibits colony growth by 50 %). Toxicophores for T. pyriformis will be identified by atomic contributions. All the above steps are the same: load data, featurize, build a model, create dataset of fragments, find contributions and visualize them Note: this time as it is regression, contributions will be in activity units, not probability<jupyter_code>DATASET_FILE ='atomic_contributions_tutorial_data/Tetrahymena_pyriformis_Work_set_OCHEM.sdf'
# create RDKit mol objects, we will need them later
mols = [m for m in Chem.SDMolSupplier(DATASET_FILE) if m is not None ]
len(mols)
loader = dc.data.SDFLoader(tasks=["IGC50"], 
                           featurizer=dc.feat.ConvMolFeaturizer(), sanitize = True)

dataset = loader.create_dataset(DATASET_FILE, shard_size=5000)
np.histogram(dataset.y)
np.random.seed(2020)
tf.random.set_seed(2020)
m = dc.models.GraphConvModel(1, mode="regression", nb_epoch=40,batch_normalize=False, restore=False)
m.fit(dataset)
TEST_DATASET_FILE = 'atomic_contributions_tutorial_data/Tetrahymena_pyriformis_Test_set_OCHEM.sdf'
loader = dc.data.SDFLoader(tasks=["IGC50"], sanitize= True,
        
                           featurizer=dc.feat.ConvMolFeaturizer())

test_dataset = loader.create_dataset(TEST_DATASET_FILE, shard_size=2000)

pred = m.predict(test_dataset)
mse=metrics.mean_squared_error(y_true=test_dataset.y, y_pred=pred)
r2=metrics.r2_score(y_true=test_dataset.y, y_pred=pred)
mse
r2 # r2 is not very high, but acceptable
# load again the training set, 
# but when featurizing it,  set per_atom_fragmentation=True

loader = dc.data.SDFLoader(tasks=[], # dont need any task
                           sanitize = True,
                           featurizer=dc.feat.ConvMolFeaturizer(per_atom_fragmentation=True)) 

frag_dataset = loader.create_dataset(DATASET_FILE, shard_size=5000)
tr = dc.trans.FlatteningTransformer(frag_dataset) # flatten dataset and add ids to each fragment
frag_dataset = tr.transform(frag_dataset)
# whole molecules
pred = m.predict(dataset)
pred = pd.DataFrame(pred, index=dataset.ids, columns=["Molecule"])  # turn to dataframe for convinience
# fragments
pred_frags = m.predict(frag_dataset)
pred_frags = pd.DataFrame(pred_frags, index=frag_dataset.ids, columns=["Fragment"])  # turn to dataframe for convinience
# merge 2 dataframes by molecule names
df = pd.merge(pred_frags, pred, right_index=True, left_index=True)
# find contribs
df['Contrib'] = df["Molecule"] - df["Fragment"]
# Lets take some molecules with moderate activity (not extremely active/inactive)
maps = vis_contribs([mol for mol in mols if float(mol.GetProp("IGC50"))>3 and float(mol.GetProp("IGC50"))<4], df)<jupyter_output><empty_output><jupyter_text>We can see, that known toxicophores are in green, namely, nitro-aromatics, halo-aromatics, long alkyl chains, aldehyde; while carboxylic group, alcohol, amino are detoxyfying, as is consistent with literature [3]  Appendix  In this tutorial we operated on sdf files. However,  If we use csv with SMILES as input , the order of atoms in the dataframe DOES NOT correspond to original atom order. If we  want to recover the original atom order for each molecule (to have it in our main dataframe), we need to use RDKit's  Chem.rdmolfiles.CanonicalRankAtoms. Here are some utils to do this. We can add a column with atom ids (as in input molecules) and use the resulting dataframe for analysis with any other  software, outside "python-rdkit-deepchem" environment<jupyter_code>def get_mapping(mols, mol_names): 
    """perform mapping:
    atom number original <-> atom number(position) 
    after ranking (both 1-based)"""
    # mols - RDKit mols
    # names  - any seq of strings
    # return list of nested lists: [[molecule, [atom , atom, ..], [...]]
    assert(len(mols)==len(mol_names))
    mapping = []
    for m,n in zip(mols, mol_names):
        atom_ids = [i+1 for i in list(
            Chem.rdmolfiles.CanonicalRankAtoms(m))]
        mapping.append([n, atom_ids
                   ])
    
    return mapping
    
def append_atomid_col(df, mapping):
    # add column with CORRECT atom number(position)
    for i in mapping:
        df.loc[i[0],"AtomID"] = i[1]
    return df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08279419478929882, 'fraction_numerical': 0.008742787200559538, 'mean_word_length': 4.643315244203256, 'pattern_counts': {'":': 2, '<': 41, '<?xml version=': 0, '>': 41, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '634967', 'license': 'MIT', 'n_tokens_mistral': 3621, 'n_tokens_neox': 3409, 'n_words': 1165, 'path': 'examples/tutorials/28_Atomic_contribtuions_tutorial.ipynb', 'repo_name': 'maclandrol/deepchem', 'size': '11458'}