<jupyter_start><jupyter_text>R-group AnalysisIn this notebook, we will use an R-group decomposition to perform exploratory data analysis on a dataset from the ChEMBL database. We will define a common scaffold for a set of molecules and use the function [rdRGroupDecomposition.RGroupDecompose](https://www.rdkit.org/docs/source/rdkit.Chem.rdRGroupDecomposition.html) from the RDKit to create a table of substituents attached at different positions on the scaffold. To perform this analysis, we will carry out the following steps.- Read the input data from a SMILES file- Cluster the input data to identify similar molecules- View the largest cluster and identify the common scaffold- Perform the R-group decomposition- View the R-groups and their frequency  SetupInstall the necessary Python packages<jupyter_code>!pip install useful_rdkit_utils mols2grid<jupyter_output>Requirement already satisfied: useful_rdkit_utils in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (0.1.5.post0.dev1)
Requirement already satisfied: mols2grid in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (0.2.2)
Requirement already satisfied: Pillow in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (from useful_rdkit_utils) (8.3.1)
Requirement already satisfied: pandas in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (from useful_rdkit_utils) (1.3.2)
Requirement already satisfied: pystow in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (from useful_rdkit_utils) (0.2.8)
Requirement already satisfied: jupyter in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (from useful_rdkit_utils) (1.0.0)
Requirement already satisfied: py3DMol in /opt/anaconda3/envs/rdkit_2021_08/lib/python3.9/site-packages (from useful_rdkit_utils) (1.7.0)
Requirement already satisfied: pytest in /opt/anaconda3/envs/rdki[...]<jupyter_text>Import the necessary Python packages<jupyter_code>from rdkit import Chem
from rdkit.Chem import rdFMCS
from rdkit.Chem import rdRGroupDecomposition
import pandas as pd
import mols2grid
import useful_rdkit_utils as uru
from ipywidgets import interact<jupyter_output><empty_output><jupyter_text>Read the Input DataWe begin by reading a set of molecules from a csv file.  The csv file has the SMILES, Name, and pIC50 for the molecules.  We'd like to find molecules with a common scaffold, then use that scaffold and decompose the molecules into R-groups around the common scaffold.<jupyter_code>df = pd.read_csv("https://raw.githubusercontent.com/PatWalters/practical_cheminformatics_tutorials/main/data/CHEMBL1075104.smi",names=["SMILES","Name","pIC50"])<jupyter_output><empty_output><jupyter_text>Cluster the Input Data to Identify Similar MoleculesWe'll begin by adding an RDKit molecule and a fingerprint to each molecule.<jupyter_code>df['mol'] = df.SMILES.apply(Chem.MolFromSmiles)
df['fp'] = df.mol.apply(uru.mol2morgan_fp)<jupyter_output><empty_output><jupyter_text>In order to find molecules with a common scaffold, we'll cluster the molecules using the [Butina clustering](https://www.rdkit.org/docs/source/rdkit.ML.Cluster.Butina.html) method in the RDKit.<jupyter_code>df['cluster'] = uru.taylor_butina_clustering(df.fp.values)<jupyter_output><empty_output><jupyter_text>Let's look at the sizes of the clusters.<jupyter_code>df.cluster.value_counts()<jupyter_output><empty_output><jupyter_text>View the Largest Cluster and Identify the Common ScaffoldLet's look at cluster 0, the largest cluster.<jupyter_code>mols2grid.display(df.query("cluster == 0"),mol_col="mol",subset=["img"])<jupyter_output><empty_output><jupyter_text>At first glance, it looks like the structure below is the scaffold.  Try searching the structures with that scaffold. - Click on the arrow next to the magnifying glass in the grid above and select **SMARTS**. - Paste the SMILES from the box below into the **Search** box.<jupyter_code>Chem.MolFromSmiles("[*]C(=O)c1ccc(N[*])c(O[*])c1")<jupyter_output><empty_output><jupyter_text>We can sketch the scaffold with R-groups in place and save the core with R-groups to a molfile. I sketched this with [Marvin Sketch](https://marvinjs-demo.chemaxon.com/latest/examples/demo.html) but you should be able to use any chemical sketcher.<jupyter_code>mb = """
  Mrv2004 04242220562D          

 13 13  0  0  0  0            999 V2000
   -1.0632    2.0028    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.7777    1.5903    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.7777    0.7652    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.0632    0.3527    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3488    0.7652    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3488    1.5903    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.3321    1.9834    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
   -2.4922    2.0028    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.4922    2.8278    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.3057    0.2630    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.2067    1.5903    0.0000 R#  0  0  0  0  0  0  0  0  0  0  0  0
    1.0466    1.5709    0.0000 R#  0  0  0  0  0  0  0  0  0  0  0  0
    0.3057   -0.5620    0.0000 R#  0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0  0  0  0
  2  3  1  0  0  0  0
  3  4  2  0  0  0  0
  4  5  1  0  0  0  0
  5  6  2  0  0  0  0
  1  6  1  0  0  0  0
  6  7  1  0  0  0  0
  2  8  1  0  0  0  0
  8  9  2  0  0  0  0
  5 10  1  0  0  0  0
  8 11  1  0  0  0  0
  7 12  1  0  0  0  0
 10 13  1  0  0  0  0
M  RGP  3  11   1  12   2  13   3
M  END
"""<jupyter_output><empty_output><jupyter_text>We can then convert the molfile to an RDKit molecule.<jupyter_code>core = Chem.MolFromMolBlock(mb)<jupyter_output><empty_output><jupyter_text>Here's the core molecule.<jupyter_code>core<jupyter_output><empty_output><jupyter_text>Now let's create a new Pandas dataframe from the molecules in cluster 0.<jupyter_code>df_0 = df.query("cluster == 0").copy()<jupyter_output><empty_output><jupyter_text>We'll add an index column to keep track of things<jupyter_code>df_0['index'] = range(0,len(df_0))<jupyter_output><empty_output><jupyter_text>Perform the R-group DecompositionAs mentioned above, we're using the function [rdRGroupDecomposition.RGroupDecompose](https://www.rdkit.org/docs/source/rdkit.Chem.rdRGroupDecomposition.html) from the RDKit. Note that RDKit returns two values from this function.  - rgd - a dictionary containing the results of the R-group decomposition. This dictionary has keys containing the core with a key "Core", and the R-groups in keys named "R1", "R2", etc. Each dictionary key links to a list of cores or R-groups corresponding input molecules that matched the core (didn't fail).- failed - a list containing the indices of molecules that did not match the core.<jupyter_code>rgd,failed = rdRGroupDecomposition.RGroupDecompose([core],df_0.mol.values,asRows=False)<jupyter_output>RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches
RDKit DEBUG: [11:45:20] No core matches
[11:45:20] No core matches<jupyter_text>Let's look at the first core, all the rest should be the same.<jupyter_code>rgd['Core'][0]<jupyter_output><empty_output><jupyter_text>Note this core has 4 R-groups, while our input core has 3.  Scroll up and see for yourself.  This is very cool!  The RDKit noticed another point of substitution and added an R-group for us.  Here are the indices of the molecules that failed.<jupyter_code>failed<jupyter_output><empty_output><jupyter_text>Let's look at the molecules that failed and try to figure out why. These results make sense.  - Molecule 4 has a **Br** in place of the **O** in the core- Molecule 5 has a **methyl** in place of the oxygen in the core- Molecule 6 does not have an **O** connected to the core- etcIn some cases, we may want to go back and modify our core definition to include these molecules.  In this case, we'll just skip them.<jupyter_code>mols2grid.display(df_0.iloc[failed],mol_col="mol",subset=["img","index"])<jupyter_output><empty_output><jupyter_text>Now we want to delete the rows corresponding to molecules that failed the R-group decomposition from **df_O**. Since we created **df_O** from **df**, it inherited the original dataframe's index.  We'll call **reset_index** to create a new index that simply numbers the rows from 0 to len(df_0).<jupyter_code>df_0.reset_index(inplace=True)<jupyter_output><empty_output><jupyter_text>Delete the failed rows from **df_0**.<jupyter_code>df_0 = df_0.drop(failed,axis=0)<jupyter_output><empty_output><jupyter_text>Get the names of the R-groups<jupyter_code>r_groups = [x for x in rgd.keys() if x != "Core"]
r_groups<jupyter_output><empty_output><jupyter_text>Add the SMILES for the R-groups to **df_0**<jupyter_code>for r in sorted(r_groups):
    df_0[r] = rgd[r]
    df_0[r] = df_0[r].apply(Chem.MolToSmiles)<jupyter_output><empty_output><jupyter_text>View the R-groups and Their FrequencyCreate an interactive viewer for the R-groups that shows how many times each R-group was used. Use the menu below to toggle between R1, R2, R3, and R4.<jupyter_code>@interact(rg=r_groups)
def display_rgroups(rg):
    vc_df = uru.value_counts_df(df_0,rg)
    return mols2grid.display(vc_df,smiles_col=rg,subset=["img","count"])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09221146085552866, 'fraction_numerical': 0.07677562550443906, 'mean_word_length': 4.031979695431472, 'pattern_counts': {'":': 0, '<': 87, '<?xml version=': 0, '>': 87, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, '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': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '431037', 'license': 'MIT', 'n_tokens_mistral': 3975, 'n_tokens_neox': 3610, 'n_words': 1199, 'path': 'sar_analysis/R_group_analysis.ipynb', 'repo_name': 'PatWalters/practical_cheminformatics_tutorials', 'size': '9912'}