01-data_extr_property_calc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: Pandarallel will run on 8 workers.\n",
      "INFO: Pandarallel will use Memory file system to transfer data between the main process and workers.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from matplotlib.colors import ListedColormap\n",
    "from faerun import Faerun\n",
    "from rdkit import Chem \n",
    "from rdkit.Chem import Lipinski, Descriptors, rdMolDescriptors, AllChem, PandasTools, rdchem\n",
    "from pandarallel import pandarallel\n",
    "import numpy as np\n",
    "pandarallel.initialize(progress_bar=False)\n",
    "import tmap as tm\n",
    "from map4 import MAP4Calculator\n",
    "import os\n",
    "import seaborn as sns\n",
    "import joblib\n",
    "import matplotlib.pyplot as plt\n",
    "from rdkit.Chem.AtomPairs import Pairs\n",
    "folder = \"/data/coconut/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "bacteria_names = [\"bacillus\", \"bacta\", \"bacteria\"]\n",
    "fungi_names = [\"aspergillus\"]\n",
    "def simple_tax(textTaxa):\n",
    "    Taxa_clean = []\n",
    "    textTaxa = textTaxa.split(\",\")\n",
    "    for tax in textTaxa:\n",
    "        tax = tax.strip().replace(\"[\", \"\").replace(\"]\", \"\").replace(\" \",\"\").replace(\"\\t\",\"\")\n",
    "        #print(tax)\n",
    "        Taxa_clean.append(tax)\n",
    "    \n",
    "    if \"Homosapiens\" in Taxa_clean:\n",
    "        #print(\"found Homo sapiens\")\n",
    "        simple_tax = \"Homo_sapiens\"\n",
    "    elif \"animal\" in Taxa_clean or \"animals\" in Taxa_clean:\n",
    "        simple_tax = \"animal\"\n",
    "    elif \"bacteria\" in Taxa_clean:\n",
    "        simple_tax = \"bacteria\"\n",
    "    elif \"fungi\" in Taxa_clean:\n",
    "        simple_tax = \"fungi\"\n",
    "    elif \"plants\" in Taxa_clean or \"plant\" in Taxa_clean:\n",
    "        simple_tax = \"plants\"\n",
    "    elif \"marine\" in Taxa_clean:\n",
    "        simple_tax = \"marine\"\n",
    "    else:\n",
    "        simple_tax = \"other\"\n",
    "        for tax in Taxa_clean:\n",
    "            for bacteria_name in bacteria_names:\n",
    "                if bacteria_name in tax.lower():\n",
    "                    simple_tax = \"bacteria\"\n",
    "                    break\n",
    "        for tax in Taxa_clean:\n",
    "            for fungi_name in fungi_names:\n",
    "                if fungi_name in tax.lower():\n",
    "                    simple_tax = \"fungi\"\n",
    "                    break\n",
    "    return simple_tax \n",
    "    \n",
    "def clean_source(citationDOI):\n",
    "    if len(citationDOI) >=10:\n",
    "        return 1\n",
    "    elif len(citationDOI) >=3:\n",
    "        return -1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def canonical_smiles(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    if mol:\n",
    "        smiles = Chem.MolToSmiles(mol)\n",
    "        return smiles\n",
    "    else:\n",
    "        return np.nan\n",
    "    \n",
    "MAP4 = MAP4Calculator(dimensions=1024)\n",
    "def calc_map4(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    smiles = Chem.MolToSmiles(mol, isomericSmiles=False)\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    map4 = MAP4.calculate(mol)\n",
    "    return np.array(map4)\n",
    "\n",
    "def calc_hba(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    hba = Lipinski.NumHAcceptors(mol)\n",
    "    return hba\n",
    "\n",
    "def calc_mw(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    mw = Descriptors.ExactMolWt(mol)\n",
    "    return mw\n",
    "\n",
    "def calc_hbd(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    hbd = Lipinski.NumHDonors(mol)\n",
    "    return hbd\n",
    "\n",
    "def calc_fcsp3(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    fcsp3 = Lipinski.FractionCSP3(mol)\n",
    "    return fcsp3\n",
    "\n",
    "def calc_logp(smiles):\n",
    "    mol = Chem.MolFromSmiles(smiles)\n",
    "    logp = Descriptors.MolLogP(mol)\n",
    "    return logp\n",
    "\n",
    "def is_lipinski(row):\n",
    "    rules = 0\n",
    "    if row.MW >= 500:\n",
    "        rules += 1\n",
    "    if row.HBA > 10:\n",
    "        rules += 1\n",
    "    if row.HBD > 5:\n",
    "        rules += 1\n",
    "    if row.aLogP > 5:\n",
    "        rules += 1\n",
    "    if rules > 1 :\n",
    "        return 0\n",
    "    else:\n",
    "        return 1\n",
    "\n",
    "def create_label(row):\n",
    "    smiles = row.SMILES\n",
    "    ID = row.coconut_id.strip()\n",
    "    link = f'https://coconut.naturalproducts.net/compound/coconut_id/{ID}'\n",
    "    label = f'{smiles}__<a target=\"_blank\" href={link}>{ID}</a>'\n",
    "    #label = smiles+\"__\"+NPAID+\"__\"+link\n",
    "    return label\n",
    "\n",
    "def has_aminoacid(smiles):\n",
    "    m = Chem.MolFromSmiles(smiles)\n",
    "    # generic dipeptide (does not hit Pro)\n",
    "    if m.HasSubstructMatch(Chem.MolFromSmarts('[NX3,NX4+][CH1,CH2][CX3](=[OX1])[NX3,NX4+][CH1,CH2][CX3](=[OX1])[O,N]')):\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "def has_sugar(smiles):\n",
    "    m = Chem.MolFromSmiles(smiles)\n",
    "    # acetal moiety\n",
    "    if m.HasSubstructMatch(Chem.MolFromSmarts('[CR][OR][CHR]([OR0,NR0])[CR]')):\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "def calc_ecfp4(smiles):\n",
    "    m = Chem.MolFromSmiles(smiles)\n",
    "    return AllChem.GetMorganFingerprintAsBitVect(m,2,nBits=1024)\n",
    "\n",
    "def calc_ap(smiles, l=1024):\n",
    "    m = Chem.MolFromSmiles(smiles)\n",
    "    sparse = Pairs.GetAtomPairFingerprintAsBitVect(m).GetOnBits()\n",
    "    ap = np.zeros(l)\n",
    "    for i in sparse:\n",
    "        ap[i%l] = 1\n",
    "    return ap\n",
    "\n",
    "\n",
    "origins = [\"plants\", \"fungi\", \"bacteria\", \"animal\", \"Homo_sapiens\"]\n",
    "def origin(simple_tax):\n",
    "    if simple_tax in origins:\n",
    "        return int(origins.index(simple_tax))\n",
    "    else:\n",
    "        return 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DB extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(folder+\"coconut_futher_analysis.pkl\"):\n",
    "    sdfFile = os.path.join(\"data/COCONUT_DB.sdf\")\n",
    "    coconut_ = PandasTools.LoadSDF(sdfFile)\n",
    "    print(len(coconut_))\n",
    "\n",
    "    coconut_ = coconut_.query(\"textTaxa != '[notax]'\").copy()\n",
    "    print(len(coconut_))\n",
    "\n",
    "    coconut_[\"canonical_smiles\"] = coconut_.SMILES.parallel_map(canonical_smiles)\n",
    "    coconut_ = coconut_.dropna().copy()\n",
    "    print(len(coconut_))\n",
    "\n",
    "    coconut_[\"has_source\"] = coconut_.citationDOI.map(clean_source)\n",
    "    #coconut = coconut.query(\"has_source == 1\").copy()\n",
    "    print(len(coconut_))\n",
    "    pd.to_pickle(coconut_, folder + \"coconut_futher_analysis.pkl\")\n",
    "else:\n",
    "    coconut_ = pd.read_pickle(folder + \"coconut_futher_analysis.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RDKit ERROR: [12:45:03] Explicit valence for atom # 34 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:03] Explicit valence for atom # 65 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:04] Explicit valence for atom # 25 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:04] Explicit valence for atom # 37 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:04] Explicit valence for atom # 36 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 16 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 8 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 8 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 44 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 33 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 39 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:05] Explicit valence for atom # 18 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:06] Explicit valence for atom # 35 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:06] Explicit valence for atom # 20 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:07] Explicit valence for atom # 32 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:07] Explicit valence for atom # 35 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:07] Explicit valence for atom # 11 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:07] Explicit valence for atom # 33 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:07] Explicit valence for atom # 26 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:07] Explicit valence for atom # 3 O, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:08] Explicit valence for atom # 22 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:08] Explicit valence for atom # 35 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:09] Explicit valence for atom # 36 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:09] Explicit valence for atom # 16 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:09] Explicit valence for atom # 18 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:09] Explicit valence for atom # 26 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:09] Explicit valence for atom # 0 O, 3, is greater than permitted\n",
      "RDKit ERROR: [12:45:09] Explicit valence for atom # 8 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:10] Explicit valence for atom # 68 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:10] Explicit valence for atom # 8 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:10] Explicit valence for atom # 27 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:10] Explicit valence for atom # 8 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:11] Explicit valence for atom # 8 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:11] Explicit valence for atom # 60 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:11] Explicit valence for atom # 16 O, 3, is greater than permitted\n",
      "RDKit ERROR: [12:45:11] Explicit valence for atom # 20 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:12] Explicit valence for atom # 27 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:12] Explicit valence for atom # 17 N, 5, is greater than permitted\n",
      "RDKit ERROR: [12:45:12] Explicit valence for atom # 8 N, 5, is greater than permitted\n"
     ]
    }
   ],
   "source": [
    "coconut_[\"canonical_smiles\"] = coconut_.SMILES.parallel_map(canonical_smiles)\n",
    "coconut_ = coconut_.dropna().copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "coconut_extra = coconut_.query(\"has_source == 0 or has_source == -1\").copy()\n",
    "coconut = coconut_.query(\"has_source == 1\").copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "400837\n",
    "135091\n",
    "135091\n",
    "135091"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "67730\n"
     ]
    }
   ],
   "source": [
    "print(len(coconut))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "135052\n"
     ]
    }
   ],
   "source": [
    "print(len(coconut_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "67322\n"
     ]
    }
   ],
   "source": [
    "print(len(coconut_extra))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Properties calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(folder+\"coconut_futher_analysis_prop.pkl\"):\n",
    "    coconut_extra[\"simple_tax\"] = coconut_extra.textTaxa.map(simple_tax)\n",
    "    coconut_extra[\"MAP4\"] = coconut_extra.canonical_smiles.parallel_map(calc_map4)\n",
    "    coconut_extra[\"fcsp3\"] = coconut_extra.canonical_smiles.parallel_map(calc_fcsp3)\n",
    "    coconut_extra[\"HBA\"] = coconut_extra.canonical_smiles.parallel_map(calc_hba)\n",
    "    coconut_extra[\"HBD\"] = coconut_extra.canonical_smiles.parallel_map(calc_hbd)\n",
    "    coconut_extra[\"aLogP\"] = coconut_extra.canonical_smiles.parallel_map(calc_logp)\n",
    "    coconut_extra[\"MW\"] = coconut_extra.canonical_smiles.parallel_map(calc_mw)\n",
    "    coconut_extra[\"isLipinski\"] = coconut_extra.apply(is_lipinski, axis=1)\n",
    "    coconut_extra[\"isPeptide\"] = coconut_extra.canonical_smiles.parallel_map(has_aminoacid)\n",
    "    coconut_extra[\"hasSugar\"] = coconut_extra.canonical_smiles.parallel_map(has_sugar)\n",
    "    coconut_extra[\"TMAPlabel\"] = coconut_extra.apply(create_label, axis=1)\n",
    "    coconut_extra[\"origin\"] = coconut_extra.simple_tax.map(origin)\n",
    "    coconut_extra[\"ecfp4\"] = coconut_extra.canonical_smiles.parallel_map(calc_ecfp4)\n",
    "    coconut_extra[\"ap\"] = coconut_extra.canonical_smiles.parallel_map(calc_ap)\n",
    "    pd.to_pickle(coconut_extra, folder + \"coconut_futher_analysis_prop.pkl\")\n",
    "else:\n",
    "    coconut_extra = pd.read_pickle(folder + \"coconut_futher_analysis_prop.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "124         plants\n",
       "2254        plants\n",
       "2307        plants\n",
       "2770      bacteria\n",
       "3325      bacteria\n",
       "            ...   \n",
       "401608      plants\n",
       "401609    bacteria\n",
       "401617      plants\n",
       "401621      plants\n",
       "401623      plants\n",
       "Name: simple_tax, Length: 67322, dtype: object"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coconut_extra.simple_tax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(folder+\"coconut_prop.pkl\"):    \n",
    "    coconut[\"simple_tax\"] = coconut.textTaxa.map(simple_tax)\n",
    "    coconut[\"MAP4\"] = coconut.canonical_smiles.parallel_map(calc_map4)\n",
    "    coconut[\"fcsp3\"] = coconut.canonical_smiles.parallel_map(calc_fcsp3)\n",
    "    coconut[\"HBA\"] = coconut.canonical_smiles.parallel_map(calc_hba)\n",
    "    coconut[\"HBD\"] = coconut.canonical_smiles.parallel_map(calc_hbd)\n",
    "    coconut[\"aLogP\"] = coconut.canonical_smiles.parallel_map(calc_logp)\n",
    "    coconut[\"MW\"] = coconut.canonical_smiles.parallel_map(calc_mw)\n",
    "    coconut[\"isLipinski\"] = coconut.apply(is_lipinski, axis=1)\n",
    "    coconut[\"isPeptide\"] = coconut.canonical_smiles.parallel_map(has_aminoacid)\n",
    "    coconut[\"hasSugar\"] = coconut.canonical_smiles.parallel_map(has_sugar)\n",
    "    coconut[\"TMAPlabel\"] = coconut.apply(create_label, axis=1)\n",
    "    coconut[\"origin\"] = coconut.simple_tax.map(origin)\n",
    "    coconut[\"ecfp4\"] = coconut.canonical_smiles.parallel_map(calc_ecfp4)\n",
    "    coconut[\"ap\"] = coconut.canonical_smiles.parallel_map(calc_ap)\n",
    "    pd.to_pickle(coconut, folder + \"coconut_prop.pkl\")\n",
    "else:\n",
    "    coconut = pd.read_pickle(folder + \"coconut_prop.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['coconut_id', 'inchi', 'inchikey', 'SMILES', 'sugar_free_smiles',\n",
       "       'molecular_formula', 'molecular_weight', 'citationDOI', 'textTaxa',\n",
       "       'name', 'synonyms', 'NPL_score', 'number_of_carbons',\n",
       "       'number_of_nitrogens', 'number_of_oxygens', 'number_of_rings',\n",
       "       'total_atom_number', 'bond_count', 'found_in_databases',\n",
       "       'murko_framework', 'alogp', 'apol', 'topoPSA', 'ID', 'ROMol',\n",
       "       'has_source', 'canonical_smiles', 'simple_tax', 'MAP4', 'fcsp3', 'HBA',\n",
       "       'HBD', 'aLogP', 'MW', 'isLipinski', 'isPeptide', 'hasSugar',\n",
       "       'TMAPlabel', 'origin', 'has_lipo_chain', 'has_ring', 'has_sterol',\n",
       "       'n_rings_nosugar', 'mw_nolipidtail_nosugar', 'smarts', 'is_terminal',\n",
       "       'ecfp4', 'ap'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coconut.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}