02b_Adding_reaction_annotations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import cobra\n",
    "from cobra import Model, Reaction, Metabolite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Abbreviation</th>\n",
       "      <th>Name</th>\n",
       "      <th>Reaction</th>\n",
       "      <th>Reactions with MetNames</th>\n",
       "      <th>GPR</th>\n",
       "      <th>Genes</th>\n",
       "      <th>Protein</th>\n",
       "      <th>Subsystems</th>\n",
       "      <th>Reversible</th>\n",
       "      <th>Lower bound</th>\n",
       "      <th>Upper bound</th>\n",
       "      <th>Objective</th>\n",
       "      <th>Confidence Score</th>\n",
       "      <th>EC. Number</th>\n",
       "      <th>Notes</th>\n",
       "      <th>References</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rxn03253</td>\n",
       "      <td>acyl-CoA dehydrogenase (decanoyl-CoA)</td>\n",
       "      <td>cpd00015[c] + cpd03128[c] -&gt; cpd00982[c] + cpd...</td>\n",
       "      <td>FAD + Decanoyl-CoA -&gt; FADH2 + (2E)-Decenoyl-CoA</td>\n",
       "      <td>((PA14_52900) or (PA14_06600) or (PA14_31580))</td>\n",
       "      <td>PA14_52900 PA14_06600 PA14_31580</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Fatty acid metabolism</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>*JAB 06/23/14 - switched out NAD for FAD based...</td>\n",
       "      <td>Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>rxn02720</td>\n",
       "      <td>acyl-CoA dehydrogenase (dodecanoyl-CoA)</td>\n",
       "      <td>cpd00015[c] + cpd01260[c] -&gt; cpd00982[c] + cpd...</td>\n",
       "      <td>FAD + Lauroyl-CoA -&gt; FADH2 + (2E)-Dodecenoyl-CoA</td>\n",
       "      <td>((PA14_52900) or (PA14_06600) or (PA14_31580))</td>\n",
       "      <td>PA14_52900 PA14_06600 PA14_31580</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Fatty acid metabolism</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>*JAB 06/23/14 - switched out NAD for FAD based...</td>\n",
       "      <td>Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Abbreviation                                     Name  \\\n",
       "0     rxn03253    acyl-CoA dehydrogenase (decanoyl-CoA)   \n",
       "1     rxn02720  acyl-CoA dehydrogenase (dodecanoyl-CoA)   \n",
       "\n",
       "                                            Reaction  \\\n",
       "0  cpd00015[c] + cpd03128[c] -> cpd00982[c] + cpd...   \n",
       "1  cpd00015[c] + cpd01260[c] -> cpd00982[c] + cpd...   \n",
       "\n",
       "                            Reactions with MetNames  \\\n",
       "0   FAD + Decanoyl-CoA -> FADH2 + (2E)-Decenoyl-CoA   \n",
       "1  FAD + Lauroyl-CoA -> FADH2 + (2E)-Dodecenoyl-CoA   \n",
       "\n",
       "                                              GPR  \\\n",
       "0  ((PA14_52900) or (PA14_06600) or (PA14_31580))   \n",
       "1  ((PA14_52900) or (PA14_06600) or (PA14_31580))   \n",
       "\n",
       "                              Genes Protein             Subsystems  \\\n",
       "0  PA14_52900 PA14_06600 PA14_31580     NaN  Fatty acid metabolism   \n",
       "1  PA14_52900 PA14_06600 PA14_31580     NaN  Fatty acid metabolism   \n",
       "\n",
       "   Reversible  Lower bound  Upper bound  Objective  Confidence Score  \\\n",
       "0           0          0.0       1000.0          0               NaN   \n",
       "1           0          0.0       1000.0          0               NaN   \n",
       "\n",
       "  EC. Number                                              Notes  \\\n",
       "0        NaN  *JAB 06/23/14 - switched out NAD for FAD based...   \n",
       "1        NaN  *JAB 06/23/14 - switched out NAD for FAD based...   \n",
       "\n",
       "                                          References  \n",
       "0  Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...  \n",
       "1  Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reactions = pd.read_excel('reactions_extended.xlsx')\n",
    "reactions.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Abbreviation</th>\n",
       "      <th>Name</th>\n",
       "      <th>Reaction</th>\n",
       "      <th>GPR</th>\n",
       "      <th>Subsystems</th>\n",
       "      <th>Reversible</th>\n",
       "      <th>Lower bound</th>\n",
       "      <th>Upper bound</th>\n",
       "      <th>Objective</th>\n",
       "      <th>EC. Number</th>\n",
       "      <th>Notes</th>\n",
       "      <th>References</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rxn03253</td>\n",
       "      <td>acyl-CoA dehydrogenase (decanoyl-CoA)</td>\n",
       "      <td>cpd00015[c] + cpd03128[c] -&gt; cpd00982[c] + cpd...</td>\n",
       "      <td>((PA14_52900) or (PA14_06600) or (PA14_31580))</td>\n",
       "      <td>Fatty acid metabolism</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>*JAB 06/23/14 - switched out NAD for FAD based...</td>\n",
       "      <td>Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>rxn02720</td>\n",
       "      <td>acyl-CoA dehydrogenase (dodecanoyl-CoA)</td>\n",
       "      <td>cpd00015[c] + cpd01260[c] -&gt; cpd00982[c] + cpd...</td>\n",
       "      <td>((PA14_52900) or (PA14_06600) or (PA14_31580))</td>\n",
       "      <td>Fatty acid metabolism</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>*JAB 06/23/14 - switched out NAD for FAD based...</td>\n",
       "      <td>Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Abbreviation                                     Name  \\\n",
       "0     rxn03253    acyl-CoA dehydrogenase (decanoyl-CoA)   \n",
       "1     rxn02720  acyl-CoA dehydrogenase (dodecanoyl-CoA)   \n",
       "\n",
       "                                            Reaction  \\\n",
       "0  cpd00015[c] + cpd03128[c] -> cpd00982[c] + cpd...   \n",
       "1  cpd00015[c] + cpd01260[c] -> cpd00982[c] + cpd...   \n",
       "\n",
       "                                              GPR             Subsystems  \\\n",
       "0  ((PA14_52900) or (PA14_06600) or (PA14_31580))  Fatty acid metabolism   \n",
       "1  ((PA14_52900) or (PA14_06600) or (PA14_31580))  Fatty acid metabolism   \n",
       "\n",
       "   Reversible  Lower bound  Upper bound  Objective EC. Number  \\\n",
       "0           0          0.0       1000.0          0        NaN   \n",
       "1           0          0.0       1000.0          0        NaN   \n",
       "\n",
       "                                               Notes  \\\n",
       "0  *JAB 06/23/14 - switched out NAD for FAD based...   \n",
       "1  *JAB 06/23/14 - switched out NAD for FAD based...   \n",
       "\n",
       "                                          References  \n",
       "0  Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...  \n",
       "1  Guzik, M. W., Narancic, T., Ilic-Tomic, T., Vo...  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reactions = reactions.drop(['Reactions with MetNames', 'Genes', 'Protein', 'Confidence Score'], axis=1)\n",
    "reactions.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use libSBML for adding reaction annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from libsbml import * \n",
    "reader = SBMLReader()\n",
    "document = reader.readSBMLFromFile('created_PA14.xml')\n",
    "model = document.getModel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Species: 1295\n",
      "Reactions:  1518\n",
      "Compartments:  2\n"
     ]
    }
   ],
   "source": [
    "print('Species:', model.getNumSpecies())\n",
    "print('Reactions: ', model.getNumReactions())\n",
    "print('Compartments: ', model.getNumCompartments())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "reaction_list = model.getListOfReactions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in the BiGG dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "792\n"
     ]
    }
   ],
   "source": [
    "bigg_dict = {}\n",
    "with open('final_bigg_dict.txt') as f:\n",
    "    for line in f:\n",
    "        (key, val) = line.split('\\t')\n",
    "        key = key.replace(':', '')\n",
    "        #val = val.replace('[', '')\n",
    "        #val = val.replace(']', '')\n",
    "        #val = val.replace(\"'\", '')\n",
    "        val = val.replace('\\n', '')\n",
    "        #val = [val]\n",
    "        #print(key, val)\n",
    "        bigg_dict[key] = val\n",
    "print(len(bigg_dict))\n",
    "#print(bigg_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in the KEGG dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "958\n"
     ]
    }
   ],
   "source": [
    "kegg_dict = {}\n",
    "with open('final_kegg_dict.txt') as f:\n",
    "    for line in f:\n",
    "        (key, val) = line.split('\\t')\n",
    "        key = key.replace(':', '')\n",
    "        val = val.replace('[', '')\n",
    "        val = val.replace(']', '')\n",
    "        val = val.replace(\"'\", '')\n",
    "        val = val.replace('\\n', '')\n",
    "        val = [val]\n",
    "        #print(key, val)\n",
    "        kegg_dict[key] = val\n",
    "print(len(kegg_dict))\n",
    "#print(kegg_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rxn03253\n",
      "R00010\n",
      "R06103\n",
      "['\"KEGG:R00010\"', ' \"KEGG:R06103\"']\n"
     ]
    }
   ],
   "source": [
    "'''Testing'''\n",
    "react_test = reaction_list[0].getId()[2:]\n",
    "print(react_test)\n",
    "val_test = kegg_dict['rxn00007'][0]\n",
    "val_test = val_test.split(',')\n",
    "for val in val_test:\n",
    "    val = val.strip()\n",
    "    val = val[6:-1]\n",
    "    print(val)\n",
    "#val_test = val_test[6:-1]\n",
    "#val_test = val_test.replace('\"', '')\n",
    "#val_test = val_test.replace('KEGG:', '')\n",
    "print(val_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding the KEGG Reaction Identifier as Annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "664\n",
      "664\n"
     ]
    }
   ],
   "source": [
    "success = 0\n",
    "bigg_ids = 0\n",
    "for react in reaction_list:\n",
    "    react_id = react.getId()[2:]     #get the id of the reaction, without the Prefix 'R_'\n",
    "    if react_id in kegg_dict.keys():\n",
    "        val = kegg_dict[react_id][0]\n",
    "        #val = val[6:-1]\n",
    "        if not val: \n",
    "            continue\n",
    "        else: \n",
    "            if ',' in val:\n",
    "                #more_than_one += 1\n",
    "                val = val.split(',')\n",
    "                for val_part in val:\n",
    "                    val_part = val_part.strip()\n",
    "                    val_part = val_part[6:-1]\n",
    "                    \n",
    "                    bigg_ids += 1\n",
    "                    \n",
    "                    cv = CVTerm()\n",
    "                    cv.setQualifierType(BIOLOGICAL_QUALIFIER)\n",
    "                    cv.setBiologicalQualifierType(BQB_IS)\n",
    "                    cv.addResource('http://identifiers.org/kegg.reaction/'+val_part)\n",
    "                    \n",
    "                    status = react.addCVTerm(cv)\n",
    "                    if status == LIBSBML_OPERATION_SUCCESS:\n",
    "                        success += 1\n",
    "                    else:\n",
    "                        print(react, val)\n",
    "\n",
    "            else: \n",
    "                bigg_ids += 1\n",
    "                val = val[6:-1]\n",
    "                cv = CVTerm()\n",
    "                cv.setQualifierType(BIOLOGICAL_QUALIFIER)\n",
    "                cv.setBiologicalQualifierType(BQB_IS)\n",
    "                cv.addResource('http://identifiers.org/kegg.reaction/'+val)\n",
    "\n",
    "                status = react.addCVTerm(cv)\n",
    "                if status == LIBSBML_OPERATION_SUCCESS:\n",
    "                    success += 1\n",
    "                else:\n",
    "                    print(react, val)\n",
    "\n",
    "print(bigg_ids)           \n",
    "print(success)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newdocument = model.getSBMLDocument()\n",
    "writeSBMLToFile(newdocument,'created_PA14.xml') # 1 means success, 0 means failure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding the BiGG Reaction Identifiers as Annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "document = reader.readSBMLFromFile('created_PA14.xml')\n",
    "model = document.getModel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "reaction_list = model.getListOfReactions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "792\n",
      "792\n"
     ]
    }
   ],
   "source": [
    "success = 0\n",
    "bigg_ids = 0\n",
    "for react in reaction_list:\n",
    "    react_id = react.getId()[2:]     #get the id of the reaction, without the Prefix 'R_'\n",
    "    if react_id in bigg_dict.keys():\n",
    "        #val = bigg_dict[react_id][0] #With the original dictionary, not the manually edited\n",
    "        val = bigg_dict[react_id]\n",
    "        #val = val[6:-1]\n",
    "        if not val: \n",
    "            continue\n",
    "        else: \n",
    "            if ',' in val:\n",
    "                val = val.split(',')\n",
    "                for val_part in val:\n",
    "                    val_part = val_part.strip()\n",
    "                    val_part = val_part[6:-1]\n",
    "                    \n",
    "                    bigg_ids += 1\n",
    "                    \n",
    "                    cv = CVTerm()\n",
    "                    cv.setQualifierType(BIOLOGICAL_QUALIFIER)\n",
    "                    cv.setBiologicalQualifierType(BQB_IS)\n",
    "                    cv.addResource('http://identifiers.org/bigg.reaction/'+val_part)\n",
    "                    \n",
    "                    status = react.addCVTerm(cv)\n",
    "                    if status == LIBSBML_OPERATION_SUCCESS:\n",
    "                        success += 1\n",
    "                    else:\n",
    "                        print(react, val)\n",
    "            else: \n",
    "                bigg_ids += 1\n",
    "                #val = val[6:-1] Necessary for the original bigg_dictionary, not the manually edited one\n",
    "                cv = CVTerm()\n",
    "                cv.setQualifierType(BIOLOGICAL_QUALIFIER)\n",
    "                cv.setBiologicalQualifierType(BQB_IS)\n",
    "                cv.addResource('http://identifiers.org/bigg.reaction/'+val)\n",
    "\n",
    "                status = react.addCVTerm(cv)\n",
    "                if status == LIBSBML_OPERATION_SUCCESS:\n",
    "                    success += 1\n",
    "                else:\n",
    "                    print(react, val)\n",
    "\n",
    "print(bigg_ids)           \n",
    "print(success)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newdocument = model.getSBMLDocument()\n",
    "writeSBMLToFile(newdocument,'created_PA14.xml') # 1 means success, 0 means failure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Looking for missing identifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "PA14_model = cobra.io.read_sbml_model('created_PA14.xml')\n",
    "#PA14_model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "792\n",
      "726\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "found = 0\n",
    "not_found = 0\n",
    "more_than_one = 0\n",
    "more_than_one_dict = {}\n",
    "not_found_list = []\n",
    "\n",
    "for react in PA14_model.reactions:\n",
    "    if react.id in bigg_dict.keys():\n",
    "        found += 1\n",
    "        if ',' in bigg_dict[react.id][0]:\n",
    "            more_than_one += 1\n",
    "            more_than_one_dict[react.id] = bigg_dict[react.id]\n",
    "    else: \n",
    "        not_found += 1\n",
    "        not_found_list.append(react.id)\n",
    "        \n",
    "        \n",
    "        \n",
    "print(found)\n",
    "print(not_found)\n",
    "print(more_than_one)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['rxn01387', 'rxn13812', 'rJB00270', 'rxn00779', 'rxn13726', 'rJB00264', 'rJB00276', 'rxn05354', 'rxn05351', 'rxn05436']\n"
     ]
    }
   ],
   "source": [
    "print(not_found_list[:10])"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}