annotate_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the required packages \n",
    "\n",
    "import cobra # package for constraint-based metabolic modelling\n",
    "# See here for documentation: https://cobrapy.readthedocs.io/en/latest/\n",
    "\n",
    "import requests # for sending HTTP requests; accessing the uniprot website\n",
    "import xmltodict # to work with the xml requiry from uniprot \n",
    "\n",
    "import pandas as pd # for working with dataframes\n",
    "import numpy as np # for general maths operations\n",
    "\n",
    "import multiprocessing as mp # for running long tasks in parallel\n",
    "import time # for checking runtimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The model has 13096 reactions and 8400 metabolites and 3628 genes.\n"
     ]
    }
   ],
   "source": [
    "# Load the metabolic model\n",
    "model = cobra.io.read_sbml_model(\"./Human-GEM.xml\")\n",
    "# The model was obtained from Robinson et al. 2020\n",
    "# DOI: 10.1126/scisignal.aaz1482 \n",
    "\n",
    "# Print model details\n",
    "rxns = model.reactions\n",
    "mets = model.metabolites\n",
    "gens = model.genes\n",
    "print(\"The model has {} reactions and {} metabolites and {} genes.\".format(len(rxns),len(mets),len(gens)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def string_reaction(r):\n",
    "    \"\"\"Returns a chemical reaction as an easily readable string\"\"\"\n",
    "    met_dict = {}\n",
    "    met_dict_name = {}\n",
    "    for m in r.metabolites:\n",
    "        coeff = r.get_coefficient(m.id)\n",
    "        met_dict_name[m.name] = coeff\n",
    "    substrates_names = dict(filter(lambda x: x[1] < 0.0, met_dict_name.items())) \n",
    "    products_names = dict(filter(lambda x: x[1] > 0.0, met_dict_name.items())) \n",
    "    subs = [' '.join([str(substrates_names[key]),str(key)]) for key in substrates_names.keys()]\n",
    "    subs = ' + '.join(subs)\n",
    "    prods = [' '.join([str(products_names[key]),str(key)]) for key in products_names.keys()]\n",
    "    prods = ' + '.join(prods)\n",
    "    if r.lower_bound < 0.0:\n",
    "        if r.upper_bound == 0.0:\n",
    "            reaction_string = \"{} <-- {}\".format(subs,prods)\n",
    "        else:\n",
    "            reaction_string = \"{} <--> {}\".format(subs,prods)\n",
    "    else: \n",
    "        reaction_string = \"{} --> {}\".format(subs,prods)\n",
    "    return(reaction_string)\n",
    "\n",
    "def convert_ENSG_to_uniprot(ENSG_ID):\n",
    "    \"\"\"Convert ENSEMBL Id to Uniprot Id\"\"\"\n",
    "    server = \"http://www.uniprot.org/uniprot/?query={}&format=xml\".format(ENSG_ID)\n",
    "    r = requests.get(server, headers={\"Content-Type\": \"text/xml\"})\n",
    "    if not r.ok:\n",
    "        r.raise_for_status()\n",
    "        sys.exit()\n",
    "    response = r.text\n",
    "    if len(response) > 0:\n",
    "        # if a response was found for the search item\n",
    "        info = xmltodict.parse(response)\n",
    "        try:\n",
    "            uni_id = info['uniprot']['entry']['accession']\n",
    "            uni_name = info['uniprot']['entry']['name']\n",
    "            if isinstance(uni_id, list):\n",
    "                uni_id = uni_id[0]\n",
    "            return uni_id,uni_name\n",
    "        except TypeError:\n",
    "            uni_id = info['uniprot']['entry'][0]['accession']\n",
    "            if isinstance(uni_id, list):\n",
    "                uni_id = uni_id[0]\n",
    "            uni_name = info['uniprot']['entry'][0]['name']\n",
    "            return uni_id,uni_name\n",
    "    else:\n",
    "        # if no result was found for the search item \n",
    "        return \"NA\",\"NA\"\n",
    "    \n",
    "def get_protein_info(r):\n",
    "    \"\"\"Get protein information from Uniprot for that reaction.\"\"\"\n",
    "    try: # Check for EC number in model annotation\n",
    "        ec = r.annotation['ec-code']\n",
    "    except KeyError:# If no EC number is found, annotate with NA\n",
    "        ec = 'NA'\n",
    "    # Update the reaction rule with Uniprot IDs\n",
    "    rxn_rule = r.gene_reaction_rule\n",
    "    u_names = [] # Enzyme name on uniprot \n",
    "    for g in r.genes:\n",
    "        # One reaction may be catalyzed by different isoforms\n",
    "        u_id, u_name = convert_ENSG_to_uniprot(g.id)\n",
    "        u_names.append(u_name)\n",
    "        rxn_rule = rxn_rule.replace(g.id,u_id)\n",
    "    unique_u_name = list(set(u_names))\n",
    "    return [rxn_rule, unique_u_name, ec]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Runtime: 10250.82 second.\n"
     ]
    }
   ],
   "source": [
    "tic = time.time()\n",
    "\n",
    "# Get protein information for all model reactions using parallel processing\n",
    "pool = mp.Pool(int(mp.cpu_count()/2)) # Adjust this to wanted # of cores\n",
    "results = pool.map(get_protein_info,[r for r in model.reactions])\n",
    "pool.close()\n",
    "\n",
    "toc = time.time()\n",
    "print(\"Runtime: {} second.\".format(round(toc-tic,2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gather all of the reaction/protein info into a dataframe \n",
    "reaction_ids = [r.id for r in model.reactions]\n",
    "reaction_str = [string_reaction(r) for r in model.reactions]\n",
    "reaction_rules = [r[0] for r in results]\n",
    "uniprot_names = [r[1] for r in results]\n",
    "reaction_ecs = [r[2] for r in results]\n",
    "\n",
    "# Save the dataframe as a csv file \n",
    "data = {\"ReactionID\":reaction_ids, \"Reaction\":reaction_str, \"UniprotID\":uniprot_names,\n",
    "        \"EC\":reaction_ecs,\"ReactionRule\":reaction_rules}\n",
    "df = pd.DataFrame(data)\n",
    "df.to_csv('model_proteins.csv',index=False,header=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Annotate the model and save the above results \n",
    "# Can now use the annotated model to access protein information directly\n",
    "for i,r in enumerate(model.reactions):\n",
    "    r.name = reaction_str[i]\n",
    "    r.annotation['ec-code'] = reaction_ecs[i]\n",
    "    r.annotation['uniprot.name'] = uniprot_names[i]\n",
    "    r.notes['protein_reaction_rule'] = reaction_rules[i]\n",
    "    \n",
    "# Save the model to file\n",
    "cobra.io.write_sbml_model(model, \"Human-GEM-annotated.xml\")\n",
    "# Use the update model for further analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMR_3905\n",
      "-1.0 ethanol + -1.0 NAD+ --> 1.0 acetaldehyde + 1.0 H+ + 1.0 NADH\n",
      "{'sbo': 'SBO:0000176', 'ec-code': ['1.1.1.1', '1.1.1.71'], 'kegg.reaction': 'R00754', 'bigg.reaction': 'ALCD2x', 'metanetx.reaction': 'MNXR95725', 'uniprot.name': ['HOT_HUMAN', 'ADH1G_HUMAN', 'ADHX_HUMAN', 'PTGR3_HUMAN', 'ADH1B_HUMAN', 'ADH1A_HUMAN', 'ADH6_HUMAN', 'ADH7_HUMAN', 'ADH4_HUMAN']}\n",
      "{'Confidence Level': '0', 'AUTHORS': 'PMID:10868354;PMID:12491384;PMID:12818203;PMID:14674758;PMID:15289102;PMID:15299346;PMID:15327949;PMID:15682493;PMID:15713978', 'protein_reaction_rule': 'Q8IWW8 or P28332 or Q8N4Q0 or P07327 or P40394 or P00325 or P11766 or P08319 or P00326'}\n",
      " \n",
      "HMR_3907\n",
      "-1.0 ethanol + -1.0 NADP+ --> 1.0 acetaldehyde + 1.0 H+ + 1.0 NADPH\n",
      "{'sbo': 'SBO:0000176', 'ec-code': '1.1.1.2', 'kegg.reaction': 'R00746', 'bigg.reaction': 'ALCD2y', 'metanetx.reaction': 'MNXR95726', 'uniprot.name': 'AK1A1_HUMAN'}\n",
      "{'Confidence Level': '0', 'protein_reaction_rule': 'P14550'}\n",
      " \n",
      "HMR_4097\n",
      "-1.0 acetate + -1.0 ATP + -1.0 CoA --> 1.0 acetyl-CoA + 1.0 AMP + 1.0 PPi\n",
      "{'sbo': 'SBO:0000176', 'ec-code': '6.2.1.1', 'kegg.reaction': 'R00235', 'bigg.reaction': 'ACS', 'metanetx.reaction': 'MNXR95413', 'uniprot.name': 'ACSA_HUMAN'}\n",
      "{'Confidence Level': '0', 'AUTHORS': 'PMID:10843999;PMID:11150295;PMID:2884217;PMID:4737256;PMID:5995', 'protein_reaction_rule': 'Q9NR19'}\n",
      " \n"
     ]
    }
   ],
   "source": [
    "# CHECK: reload the model and view added info\n",
    "model = cobra.io.read_sbml_model(\"./Human-GEM-annotated.xml\")\n",
    "for r in model.reactions[0:3]:\n",
    "    print(r.id)\n",
    "    print(r.name)\n",
    "    print(r.annotation)\n",
    "    print(r.notes)\n",
    "    print(\" \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8 (with cplex)",
   "language": "python",
   "name": "python3.8"
  },
  "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}