<jupyter_start><jupyter_text>Prepare Data SetsPrepare the necessary reference data set.1. Natural products from ChEMBL (ChEMBL NPs)2. Approved and experimental drugs from DrugBank (DrugBank)3. Advanced Screening Collection (Enamine)<jupyter_code>%reload_ext autoreload
%autoreload 2
# def warn(*args, **kwargs):
# warn    pass  # to silence scikit-learn warnings

import warnings
warnings.filterwarnings('ignore')
# warnings.warn = warn

# Global Imports
# from collections import Counter
# import glob
from pathlib import Path
import sys

import pandas as pd
import numpy as np
# from pandarallel import pandarallel
# pandarallel.initialize()

# from scipy.stats import median_absolute_deviation as mad

from rdkit import DataStructs
from rdkit.Chem import AllChem as Chem
# from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole

import matplotlib.pyplot as plt
import seaborn as sns

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

# Project-local Imports
PROJECT_DIR = list(Path("..").absolute().parents)[1]
sys.path.append(str(PROJECT_DIR))
import plt_style
import utils as u
from utils import lp

Path("../tmp").mkdir(exist_ok=True)<jupyter_output><empty_output><jupyter_text>ChEMBL NPTwo versions of ChEMBL26 were downloaded from the FTP server (http://ftp.ebi.ac.uk/pub/databases/chembl/ChEMBLdb/releases/chembl_26/) on 25-Mar-2020:* `chembl_26_postgresql.tar.gz`* `chembl_26.sdf.gz`The PostgreSQL database was installed locally and all entries which had ("J. Nat. Prod" OR "J Nat Prod") as the journal source were extracted with the following SQL statement:```sql = """select distinct(md.chembl_id) as Chembl_Idfrom docs d,     molecule_dictionary md,     compound_records crwhere md.molregno = cr.molregno  and cr.doc_id = d.doc_id  and d.journal in ('J. Nat. Prod.', 'J Nat Prod')"""```This list was saved as `chembl_26_natprot_ids_only.tsv` in the `Input Data` folder of this project and was used for filtering the ChEMBL_Ids from the SD file in the following procedure.* 40525 cpds.* 35026 cpds. after deglycosylation, standardization and deduplication  LoadingLoad the result of the SQL query, the ChEMBL IDs that have "J. Nat. Prod." or "J Nat Prod" listed as journal:<jupyter_code>chembl_np_ids = set(pd.read_csv("../Input Data/chembl_26_natprot_ids_only.tsv", sep="\t")["Chembl_Id"].values)
lp(chembl_np_ids)<jupyter_output>Length                                            :      40547<jupyter_text>Load the downloaded ChEMBL 26 SD file (gzipped) and keep only the natural product IDs:<jupyter_code># Takes 10min
chembl_np = u.read_sdf("/home/pahl/comas/notebooks/sdf/cpd_libs_2003/chembl_26.sdf.gz", merge_prop="chembl_id", merge_list=chembl_np_ids)
lp(chembl_np, "ChEMBL NP")
chembl_np.head()<jupyter_output>Shape ChEMBL NP                                   :      40525 /    3  [ Smiles, Name, chembl_id ]<jupyter_text>Initial StandardizationStandardization implemented in MolVS:* standardization* remove stereo information* choose largest fragment<jupyter_code># Takes 3min
chembl_np = u.apply_to_smiles(chembl_np, "Smiles", {"Smiles_Std": lambda x: u.standardize_mol(x, remove_stereo=True)})
chembl_np.drop("Smiles", axis=1, inplace=True)
chembl_np = chembl_np.rename(columns={"Smiles_Std": "Smiles"})
lp(chembl_np, "ChEMBL NP")<jupyter_output>Shape ChEMBL NP                                   :      40525 /    3  [ Name, chembl_id, Smiles ]  ( NAN values in 1 col(s) )<jupyter_text>Filters (MedChem, Duplicates)The data set is filtered for Isotopes, MedChemAtoms, MinHeavyAtoms and Duplicates.  The deduplication is performed here as well to remove as many compounds as possible at early stages of the processing.  The filter for MaxHeavyAtoms is deferred until after the deglycosylation.<jupyter_code># Takes 2min
chembl_np = u.filter_mols(chembl_np, "Smiles", ["Isotopes", "MedChemAtoms", "MinHeavyAtoms", "Duplicates"])
lp(chembl_np, "ChEMBL NP (MedChem filters)")

chembl_np.head(10)
chembl_np.to_csv("../tmp/chembl_np.tsv", sep="\t", index=False)<jupyter_output><empty_output><jupyter_text>DeglycosylationDeglycosylation is performed using the KNIME Sugar Remover node (CDK). KNIME v 4.2.3.  Here, the result of the deglycosylation is loaded again.<jupyter_code>chembl_np = pd.read_csv("../tmp/chembl_np_deglyco.tsv", sep="\t")
lp(chembl_np, "ChEMBL NP (Deglyco)")<jupyter_output>Shape ChEMBL NP (Deglyco)                         :      34479 /    3  [ Name, chembl_id, Smiles ]<jupyter_text>StandardizationCleanup after deglycosylation:* standardization* choose largest fragment<jupyter_code># Takes 1:20min
chembl_np_std = u.apply_to_smiles(chembl_np, "Smiles", {"Smiles_Std": u.standardize_mol})
chembl_np_std.drop("Smiles", axis=1, inplace=True)
chembl_np_std = chembl_np_std.rename(columns={"Smiles_Std": "Smiles"})
lp(chembl_np_std, "ChEMBL NP (Deglyco)")
chembl_np_std = u.filter_mols(chembl_np_std, "Smiles", ["MinHeavyAtoms", "MaxHeavyAtoms", "Duplicates"])
lp(chembl_np_std, "ChEMBL NP (Filters)")<jupyter_output>Applying filters (3)...<jupyter_text>Canonical Tautomer<jupyter_code># Takes 28min
chembl_np_std_tauto = u.apply_to_smiles(chembl_np_std, "Smiles", {"Smiles_Std": lambda x: u.standardize_mol(x, canonicalize_tautomer=True)}, parallel=True)
chembl_np_std_tauto.drop("Smiles", axis=1, inplace=True)
chembl_np_std_tauto = chembl_np_std_tauto.rename(columns={"Smiles_Std": "Smiles"})<jupyter_output><empty_output><jupyter_text>Filter DuplicatesRemove duplicates that might have been formed by canonicalization.<jupyter_code>chembl_np_std_tauto = u.filter_mols(chembl_np_std_tauto, "Smiles", "Duplicates")
lp(chembl_np_std_tauto, "ChEMBL NP (Filter duplicates)")

chembl_np_std_tauto.head()<jupyter_output>Applying filters (1)...<jupyter_text>Write Data Set<jupyter_code>chembl_np_std_tauto.to_csv("../tmp/chembl_np_std.tsv", sep="\t", index=False)<jupyter_output><empty_output><jupyter_text>Create Subset<jupyter_code>chembl_np = pd.read_csv("../tmp/chembl_np_std.tsv", sep="\t")
chembl_np_subset = chembl_np.sample(n=5000, random_state=0xc0ffee)
chembl_np_subset.to_csv("../tmp/chembl_np_std_subset.tsv", sep="\t", index=False)
lp(chembl_np_subset, "ChEMBL NP subset")<jupyter_output>Shape ChEMBL NP subset                            :       5000 /    3  [ Name, chembl_id, Smiles ]<jupyter_text>DrugBankApproved and experimental drugs, downloaded from DrugBank on 25-Mar-2020, v 5.1.5. (https://www.drugbank.com/)* 8479 cpds.* 7472 cpds. after standardization and deduplication  LoadingLoad both SD files of approved and experimental drugs.<jupyter_code>drugbank = u.read_sdf(["/home/pahl/comas/notebooks/sdf/drugbank/drugbank_5.1.5_appr.sdf", "/home/pahl/comas/notebooks/sdf/drugbank/drugbank_5.1.5_exp.sdf"])
drugbank =drugbank[["DATABASE_ID", "Smiles"]]
lp(drugbank, "DrugBank")

drugbank.head()<jupyter_output>RDKit ERROR: [21:16:13] ERROR: Explicit valence for atom # 6 K, 2, is greater than permitted
RDKit ERROR: [21:17:24] Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [21:17:24] ERROR: Could not sanitize molecule ending on line 19217
RDKit ERROR: [21:17:24] ERROR: Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [21:17:25] Explicit valence for atom # 10 N, 4, is greater than permitted
RDKit ERROR: [21:17:25] ERROR: Could not sanitize molecule ending on line 318021
RDKit WARNING: [21:17:26] MRV_IMPLICIT_H SGroup on atom without aromatic bonds, 6, ignored.
RDKit WARNING: [21:17:26] MRV_IMPLICIT_H SGroup on atom without aromatic bonds, 11, ignored.
RDKit WARNING: [21:17:26] MRV_IMPLICIT_H SGroup on atom without aromatic bonds, 10, ignored.
RDKit WARNING: [21:17:26] MRV_IMPLICIT_H SGroup on atom without aromatic bonds, 14, ignored.
RDKit ERROR: [21:17:25] ERROR: Explicit valence for atom # 10 N, 4, is greater than permitted
RDKit ERROR: [2[...]<jupyter_text>StandardizationStandardization implemented in MolVS:* standardization* remove stereo information* choose largest fragment* canonicalize tautomer<jupyter_code>drugbank = u.apply_to_smiles(drugbank, "Smiles", {"Smiles_Std": lambda x: u.standardize_mol(x, remove_stereo=True, canonicalize_tautomer=True)})
drugbank.drop("Smiles", axis=1, inplace=True)
drugbank = drugbank.rename(columns={"Smiles_Std": "Smiles"})
lp(drugbank, "DrugBank")<jupyter_output><empty_output><jupyter_text>Filters (Isotopes, MedChem, NumHeavyAtoms, Duplicates)The data set is filtered for Isotopes, MedChemAtoms, MinHeavyAtoms, MaxHeavyAtoms and Duplicates.<jupyter_code>drugbank = u.filter_mols(drugbank, "Smiles", ["Isotopes", "MedChemAtoms", "MinHeavyAtoms", "MaxHeavyAtoms", "Duplicates"])
lp(drugbank, "DrugBank (MedChem filters)")

drugbank.head(10)<jupyter_output>Applying filters (5)...<jupyter_text>Write Data Set<jupyter_code>drugbank.to_csv("../tmp/drugbank_std.tsv", sep="\t", index=False)
lp(drugbank, "DrugBank (MedChem filters)")

drugbank.head()<jupyter_output>Shape DrugBank (MedChem filters)                  :       7472 /    2  [ DATABASE_ID, Smiles ]<jupyter_text>Create Subset<jupyter_code>drugbank = pd.read_csv("../tmp/drugbank_std.tsv", sep="\t")
drugbank_subset = drugbank.sample(n=5000, random_state=0xc0ffee)
drugbank_subset.to_csv("../tmp/drugbank_std_subset.tsv", sep="\t", index=False)
lp(drugbank_subset, "Drugbank subset")<jupyter_output>Shape Drugbank subset                             :       5000 /    2  [ DATABASE_ID, Smiles ]<jupyter_text>Enamine Advanced Screening Collection* Enamine Advanced Screening Collection (527411 compounds)* Downloaded from [Enamine](https://enamine.net/hit-finding/compound-collections/screening-collection/advanced-collection) on 07-Dec-2020SD file extracted and gzippedRandom selection of 50k, then standardized and filtered.  LoadingLoad the downloaded gzipped SD file into a DataFrame<jupyter_code>enamine = u.read_sdf(["/home/pahl/comas/notebooks/sdf/enamine/enamine_adv.gz"])
enamine = enamine[["idnumber", "Smiles"]]
lp(enamine, "Enamine")

enamine.head()<jupyter_output>RDKit ERROR: [17:22:53] Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:22:53] ERROR: Could not sanitize molecule ending on line 24061869
RDKit ERROR: [17:22:53] ERROR: Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:22:53] Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:22:53] ERROR: Could not sanitize molecule ending on line 24061948
RDKit ERROR: [17:22:53] ERROR: Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:23:13] Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:23:13] ERROR: Could not sanitize molecule ending on line 28077193
RDKit ERROR: [17:23:13] ERROR: Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:23:35] Explicit valence for atom # 7 N, 4, is greater than permitted
RDKit ERROR: [17:23:35] ERROR: Could not sanitize molecule ending on line 32668058
RDKit ERROR: [17:23:35] ERROR: Explicit valen[...]<jupyter_text>Create SubsetA random subset of 50k cpds. is selected.<jupyter_code>enamine_subset = enamine.sample(n=50000, random_state=0xc0ffee)
enamine_subset.to_csv("../tmp/enamine_subset.tsv", sep="\t", index=False)
# enamine_subset = pd.read_csv("../tmp/enamine_subset.tsv", sep="\t")<jupyter_output><empty_output><jupyter_text>StandardizationStandardization implemented in MolVS:* standardization* remove stereo information* choose largest fragment* canonicalize tautomer<jupyter_code>enamine_subset = u.apply_to_smiles(enamine_subset, "Smiles", {"Smiles_Std": lambda x: u.standardize_mol(x, remove_stereo=True, canonicalize_tautomer=True)}, parallel=True)
enamine_subset.drop("Smiles", axis=1, inplace=True)
enamine_subset = enamine_subset.rename(columns={"Smiles_Std": "Smiles"})
lp(enamine_subset, "Enamine (Std)")<jupyter_output>INFO: Pandarallel will run on 6 workers.
INFO: Pandarallel will use Memory file system to transfer data between the main process and workers.<jupyter_text>FiltersThe data set is filtered for Isotopes, MedChemAtoms, MinHeavyAtoms, MaxHeavyAtoms and Duplicates.<jupyter_code>enamine_subset = u.filter_mols(enamine_subset, "Smiles", ["Isotopes", "MedChemAtoms", "MinHeavyAtoms", "MaxHeavyAtoms", "Duplicates"])
lp(enamine_subset, "Enamine (MedChem filters)")

enamine_subset.head(10)
enamine_subset.to_csv("../tmp/enamine_subset_std.tsv", sep="\t", index=False)
enamine_subset[["idnumber"]].to_csv("../Input Data/enamine_subset_std_ids.tsv", sep="\t", index=False)
lp(enamine_subset, "Enamine (filters)")<jupyter_output>Shape Enamine (filters)                           :      49968 /    2  [ idnumber, Smiles ]<jupyter_text>Internal CompoundsApproved and experimental drugs, downloaded from DrugBank on 25-Mar-2020, v 5.1.5.* 244 individual cpds.* standardization, removal of stereo information and MedChem filters: `internal_cpds_mc.tsv`    * this data set is used for the chemical class-similarity calculations* deduplication, based on InChIKeys: `internal_cpds_std.tsv` (168 cpds.)    * analogous workflow to the reference data sets, without deglycosylation    * this data set is used for the NP Likeness calculations and the PCA analysis  LoadingLoad Load list of 244 compounds with their Smiles.<jupyter_code>cpds = pd.read_csv("../Input Data/internal_cpds.tsv", sep="\t")
lp(cpds, "Internal compounds")

cpds.head()<jupyter_output>Shape Internal compounds                          :        244 /    3  [ Compound_Id, CpdClass, Smiles ]<jupyter_text>StandardizationStandardization implemented in MolVS:* standardization* remove stereo information* choose largest fragment* canonicalize tautomer<jupyter_code>cpds = u.apply_to_smiles(cpds, "Smiles", {"Smiles_Std": lambda x: u.standardize_mol(x, remove_stereo=True, canonicalize_tautomer=True)})
cpds.drop("Smiles", axis=1, inplace=True)
cpds = cpds.rename(columns={"Smiles_Std": "Smiles"})
lp(cpds, "Internal compounds")<jupyter_output><empty_output><jupyter_text>Filters (MedChem)The data set is filtered for Isotopes, MedChemAtoms, MinHeavyAtoms, MaxHeavyAtoms.<jupyter_code>cpds = u.filter_mols(cpds, "Smiles", ["Isotopes", "MedChemAtoms", "MinHeavyAtoms", "MaxHeavyAtoms"])
lp(cpds, "Internal compounds (MedChem filters)")

cpds.head(10)
cpds.to_csv("../tmp/internal_cpds_mc.tsv", sep="\t", index=False)
lp(cpds, "Internal compounds (MedChem filters)")

cpds.head()<jupyter_output>Shape Internal compounds (MedChem filters)        :        244 /    3  [ Compound_Id, CpdClass, Smiles ]<jupyter_text>Filters (Duplicates)The data set is filtered for Duplicates (based on InChIKeys (the stereo centers have been removed in the previous standardization step)).  This is done for consistency with the standardization and filtering of the reference sets.  This data set is used for the NP Likeness calculations and the PCA analysis.<jupyter_code>cpds = u.filter_mols(cpds, "Smiles", ["Duplicates"])
lp(cpds, "Internal compounds (Duplicates)")

cpds.head(10)
cpds.to_csv("../tmp/internal_cpds_std.tsv", sep="\t", index=False)
lp(cpds, "Internal compounds (MedChem filters)")

cpds.head()
!cp ../tmp/internal_cpds_mc.tsv "../Input Data/internal_cpds_mc.tsv"
!cp ../tmp/internal_cpds_std.tsv "../Input Data/internal_cpds_std.tsv"<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11412436798213933, 'fraction_numerical': 0.027381968612515597, 'mean_word_length': 5.888285843509724, 'pattern_counts': {'":': 12, '<': 81, '<?xml version=': 0, '>': 81, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '99245', 'license': 'MIT', 'n_tokens_mistral': 5704, 'n_tokens_neox': 5224, 'n_words': 1334, 'path': '1. Data Prep/prepare_data_sets.ipynb', 'repo_name': 'jose-manuel/code_manuscript', 'size': '15229'}