co-occurrence_of_enzyme_domains_(preparing_data).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping to Brenda to Uniprot\n",
    "First I grabbed the \"brenda_download.txt\" file from the Brenda website (need to accept their user license)\n",
    "\n",
    "First, I used awk to extact the Uniprot accessions from this file and map them to Brenda's e.c. numbers. The regular expressions was based on that descibed here: https://www.uniprot.org/help/accession_numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "###extract uniprot accessions out of brenda_download.txt\n",
    "awk -vOFS=\"\\t\" '$0==\"///\"{current_ec=\"unknown\"};$1==\"ID\"{current_ec=$2};patsplit($0,uniprot_accessions,/[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}/){for(i in uniprot_accessions)print(uniprot_accessions[i],current_ec)};' brenda_download.txt > brenda_uniprot_to_ec.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I took the Uniprot accessions from the first column of \"brenda_uniprot_to_ec.tsv\", and used the \"retrieve ID\" function on the Uniprot website, made sure the Pfam column was enabled, and exported it into \"brenda_uniprot.tsv\".\n",
    "Then we merge the data in that file with \"brenda_uniprot_to_ec.tsv\" as shown below to generate \"uniprot_to_brenda_merged.tsv\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "brenda_from_uniprot_df=pandas.read_table(\"brenda_uniprot.tsv\",dtype=\"str\")\n",
    "brenda_from_uniprot_df.drop(columns=brenda_from_uniprot_df.columns[0],inplace=True)\n",
    "brenda_from_uniprot_df\n",
    "brenda_to_ec_df=pandas.read_table(\"brenda_uniprot_to_ec.tsv\",dtype=\"str\",names=[\"uniprot_accession\",\"brenda_ec_number\"])\n",
    "brenda_to_ec_deduplicated=brenda_to_ec_df.groupby(\"uniprot_accession\")['brenda_ec_number'].apply(lambda x:\";\".join(frozenset(x)))\n",
    "brenda_uniprot_merged=pandas.merge(brenda_from_uniprot_df,brenda_to_ec_deduplicated,left_on=[\"Entry\"],right_on=[\"uniprot_accession\"])\n",
    "brenda_uniprot_merged.set_index(\"Entry\",inplace=True)\n",
    "brenda_uniprot_merged.to_csv(\"uniprot_to_brenda_merged.tsv\",sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping MIBiG to NCBI protein accessions, and those in turn to Uniprot\n",
    "Downloaded the MIBiG json files from the mibig website and extracted them into a folder. These json files do not give uniprot accessions, but only contain NCBI nucleotide accessions, and sometimes start and stop locations. So, we need to map these NCBI nucleotide accessions to NCBI protein accessions (Using NCBI's efetch tool), and then map NCBI protein accessions to Uniprot (Using Uniprot's ID mapping service)\n",
    "\n",
    "You ought to fill in your own e-mail address and api-key for NCBI:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "email_for_ncbi=\"your@email.com\"\n",
    "ncbi_api_key=\"get_this_from_your_NCBI_account\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Run this cell only once\n",
    "\n",
    "#The NCBI efetch service will sometimes time out (maybe when it's too busy?) which stops the loop in the next cell.\n",
    "#Using this cache, we do not need to request data we have already requested before.\n",
    "cache={}#map urls to response text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#If this cell fails, you may want to try running it again and see if it gets further along.\n",
    "\n",
    "###this cell queries NCBI for all the protein IDs in Mibig\n",
    "\n",
    "import requests\n",
    "import pandas\n",
    "import glob\n",
    "import json\n",
    "import re\n",
    "import sys\n",
    "\n",
    "session = requests.Session()\n",
    "session.headers.update({'User-Agent': 'programmatic access by '+email_for_ncbi,\n",
    "                       'From':email_for_ncbi,\n",
    "                       })\n",
    "ncbi_efetch_base_url = \"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?api_key=\"+ncbi_api_key+\"&email=\"+email_for_ncbi+\"&db=nuccore&rettype=ft&conwithfeat=on\"\n",
    "#the &conwithfeat=on part is not documented by NCBI, but it's required to make certain large records download that otherwise would just be empty \n",
    "\n",
    "mibig_file_list=glob.glob(\"mibig_json_2.0/*.json\")\n",
    "\n",
    "results_df=pandas.DataFrame(columns=[\"ncbi_protein_accession\",\"mibig_accession\",\"protein_id_type\",\"mibig_taxid\"])\n",
    "results_df.set_index([\"ncbi_protein_accession\",\"mibig_accession\"],inplace=True)\n",
    "\n",
    "gb_regex=re.compile(\"protein_id.+?(gb|dbj|ref|gp|emb|tpe|tpg|tpd)\\|(.{7,13}?)\\..{1,2}\\|\")\n",
    "for filename in mibig_file_list:\n",
    "#     if(int(filename[-10:-5])<500):   #useful if you need to start somewhere in the middle\n",
    "#          continue\n",
    "    with open(filename) as f:\n",
    "        mibig_file_dict=json.load(f)\n",
    "        print(\"retrieving ncbi protein IDs for {}\".format(mibig_file_dict[\"cluster\"][\"mibig_accession\"]))\n",
    "        mibig_locus_dict=mibig_file_dict[\"cluster\"][\"loci\"]\n",
    "        mibig_accession=mibig_file_dict[\"cluster\"][\"mibig_accession\"]\n",
    "        if \"start_coord\" in mibig_locus_dict:\n",
    "            request_url=ncbi_efetch_base_url+\"&id={}&seq_start={}&seq_stop={}\".format(mibig_locus_dict[\"accession\"],\n",
    "                                                                                      mibig_locus_dict[\"start_coord\"],\n",
    "                                                                                      mibig_locus_dict[\"end_coord\"]\n",
    "                                                                                     )\n",
    "        else:\n",
    "            request_url=ncbi_efetch_base_url+\"&id=\"+mibig_locus_dict[\"accession\"]\n",
    "            \n",
    "        print(\"requesting \"+request_url)\n",
    "        if request_url in cache:\n",
    "            print(\"Grabbing from cache\")\n",
    "            response_text=cache[request_url]        \n",
    "        else:\n",
    "            response_text=None\n",
    "            r=session.get(request_url)\n",
    "            if(r.ok):\n",
    "                print(\"response in \"+str(r.elapsed.total_seconds())+\" s\")\n",
    "                if(r.text):\n",
    "                    response_text=r.text\n",
    "                    cache[request_url]=response_text\n",
    "                else:\n",
    "                    print(\"no results\")\n",
    "            else:\n",
    "                print(\"error: \"+str(r.status_code))\n",
    "            time.sleep(5)\n",
    "        if(response_text):\n",
    "            total_proteins=0\n",
    "            for regex_match in gb_regex.finditer(response_text):\n",
    "                total_proteins+=1\n",
    "                protein_id_type=regex_match.group(1)\n",
    "                protein_accession=regex_match.group(2)\n",
    "                if (protein_accession,mibig_accession) in results_df.index:\n",
    "                    print(\"duplicates!!\")\n",
    "                else:  \n",
    "                    results_df.loc[(protein_accession,mibig_accession),\n",
    "                                   [\"protein_id_type\",\"mibig_taxid\"]\n",
    "                                   ]=[protein_id_type,mibig_file_dict[\"cluster\"][\"ncbi_tax_id\"]]\n",
    "                                  \n",
    "            print(\"{}: found {} protein accessions\".format(mibig_accession,total_proteins))\n",
    "            print(\"table now has {} entries\".format(len(results_df.index)))\n",
    "\n",
    "results_df.to_csv(\"ncbi_accessions_to_mibig.tsv\",sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we will use the file that this created, \"ncbi_accessions_to_mibig.tsv\", later on, and will use the dataframe in-memory in this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generate 2 files for NCBI protein accessions. One for refseq accessions and one for genbank/EMB/DDBJ accessions\n",
    "\n",
    "only_ref_df=results_df[results_df[\"protein_id_type\"]==\"ref\"].reset_index().set_index(\"ncbi_protein_accession\")\n",
    "print(\"length of submit_to_uniprot_refseq_protein.txt: {} entries\".format(len(only_ref_df)))\n",
    "only_ref_df.to_csv(\"submit_to_uniprot_refseq_protein.txt\",columns=[],header=False)\n",
    "\n",
    "not_ref_df=results_df[results_df[\"protein_id_type\"]!=\"ref\"].reset_index().set_index(\"ncbi_protein_accession\")\n",
    "print(\"length of submit_to_uniprot_embl_genbank_cds.txt: {} entries\".format(len(not_ref_df)))\n",
    "not_ref_df.to_csv(\"submit_to_uniprot_embl_genbank_cds.txt\",columns=[],header=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We submit these two files to the UniProt ID mapping service:\n",
    "\n",
    "- Use the Refseq protein input option for \"submit_to_uniprot_refseq_protein.txt\". The results of this were saved in \"mibig_uniprot_refseq_proteins.tsv\"\n",
    "- And use the EMBL/Genbank CDS input option for \"submit_to_uniprot_embl_genbank_cds.txt\". For this file there are too many accession numbers to all submit at once, so I submitted in three batches and concatenated the results (saved in \"mibig_uniprot_embl_genbank_cds_proteins_partA.tsv\", etc., as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "#combine the three separate requests to uniprot\n",
    "cat mibig_uniprot_embl_genbank_cds_proteins_partA.tsv > mibig_uniprot_embl_genbank_cds_proteins_combined.tsv\n",
    "tail -n+2 mibig_uniprot_embl_genbank_cds_proteins_partB.tsv >> mibig_uniprot_embl_genbank_cds_proteins_combined.tsv\n",
    "tail -n+2 mibig_uniprot_embl_genbank_cds_proteins_partC.tsv >> mibig_uniprot_embl_genbank_cds_proteins_combined.tsv\n",
    "tail -n+2 mibig_uniprot_embl_genbank_cds_proteins_partD.tsv >> mibig_uniprot_embl_genbank_cds_proteins_combined.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly we combine the metadata from the two different Uniprot ID mapping types, and merge that data with the MIBiG to NCBI protein accession mapping, to finally give a mapping from MIBiG BGC# to Uniprot accession in the \"uniprot_to_mibig_merged_on_protein_id.tsv\" file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Table matching protein ids to mibig has {} entries\".format(len(results_df)))\n",
    "\n",
    "uniprot_refseq_df=pandas.read_table(\"mibig_uniprot_refseq_proteins.tsv\",dtype=\"str\")\n",
    "uniprot_refseq_df.rename(columns={uniprot_refseq_df.columns[1]:\"protein_id\"},inplace=True)\n",
    "uniprot_refseq_df.set_index([\"Entry\",\"protein_id\"],inplace=True)\n",
    "print(\"Refseq table has {} entries\".format(len(uniprot_refseq_df)))\n",
    "\n",
    "uniprot_embl_genbank_df=pandas.read_table(\"mibig_uniprot_embl_genbank_cds_proteins_combined.tsv\",dtype=\"str\")\n",
    "uniprot_embl_genbank_df.rename(columns={uniprot_embl_genbank_df.columns[0]:\"protein_id\"},inplace=True)\n",
    "uniprot_embl_genbank_df.drop(columns=uniprot_embl_genbank_df.columns[2],inplace=True)\n",
    "uniprot_embl_genbank_df.set_index([\"Entry\",\"protein_id\"],inplace=True)\n",
    "uniprot_embl_genbank_df.drop_duplicates(inplace=True)\n",
    "print(\"EMBL/Genbank CDS table has {} entries\".format(len(uniprot_embl_genbank_df)))\n",
    "\n",
    "uniprot_df=pandas.concat([uniprot_refseq_df,uniprot_embl_genbank_df],verify_integrity=True)\n",
    "print(\"Combined uniprot table has {} entries\".format(len(uniprot_df)))\n",
    "\n",
    "#Merging on only protein_id. Useful for quickly checking if a uniprot ID is in mibig\n",
    "merged_df_on_protein_id=pandas.merge(uniprot_df.reset_index(),results_df.reset_index(),left_on=[\"protein_id\"],right_on=[\"ncbi_protein_accession\"])\n",
    "merged_df_on_protein_id.set_index(\"Entry\",inplace=True)\n",
    "merged_df_on_protein_id.to_csv(\"uniprot_to_mibig_merged_on_protein_id.tsv\",sep=\"\\t\")\n",
    "print(\"Merging on protein ID only leaves {} entries\".format(len(merged_df_on_protein_id)))\n",
    "\n",
    "\n",
    "#Merging on both protein_id and taxid.\n",
    "#Unfortunately we do lose some entries where NONE of the uniprot results match the taxid in mibig\n",
    "merged_df_on_protein_id_and_taxid=pandas.merge(uniprot_df.reset_index(),results_df.reset_index(),left_on=[\"protein_id\",\"Taxonomic lineage IDs\"],right_on=[\"ncbi_protein_accession\",\"mibig_taxid\"])\n",
    "merged_df_on_protein_id_and_taxid.set_index(\"Entry\",inplace=True)\n",
    "merged_df_on_protein_id_and_taxid.to_csv(\"uniprot_to_mibig_merged_on_protein_id_and_taxid.tsv\",sep=\"\\t\")\n",
    "print(\"Merging on protein ID and taxid leaves {} entries\".format(len(merged_df_on_protein_id_and_taxid)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uniprot entries curated to have catalytic activity\n",
    "Search uniprot for 'annotation:(type:\"catalytic activity\" evidence:manual)', made sure the Pfam column was turned on, and exported it as \"Uniprot_curated_catalytic_activity.tsv\".\n",
    "\n",
    "Uniprot also has an older index of proteins involved in metabolic pathways, called pathway.txt: https://www.uniprot.org/docs/pathway\n",
    "I extracted the uniprot accessions out of that file using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "egrep -o '(......)' pathway.txt | awk '{print substr($0,2,6)}'>pathway_accessions.txt\n",
    "#This is not ideal and will miss a couple of entries; in the future I should use the same regular expression as in the \"Mapping to Brenda to Uniprot\" section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "then I queried the uniprot KB using those accessions in 3 batches of less than 50,000, and concatenated the files to generate \"uniprot_pathway_complete.tsv\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-process all pfam annotations in Uniprot\n",
    "For this we actually use the protein2ipr.dat.gz file from the Interpro FTP server, which is conveniently sorted by uniprot accession, speeding up downstream processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##parsing protein2ipr.dat.gz file. uniprot IDs are already sorted alphabetically in this file\n",
    "#the output files require about 2 Gb\n",
    "\n",
    "protein2ipr_filename=\"D:/domain_analysis/protein2ipr.dat.gz\"\n",
    "output_folder=\"D:/domain_analysis/pfam_to_uniprot_files/\"\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "import gzip\n",
    "import collections\n",
    "\n",
    "def write_progress_and_clear_lists(list_dict):\n",
    "    print(\"writing progress...\")\n",
    "    for k,v in list_dict.items():\n",
    "        with gzip.open(output_folder+k+\".gz\", 'ta') as output_file: #make this gzip.open for gzipped files\n",
    "            for l in v:\n",
    "                output_file.write(l+\"\\n\")\n",
    "            v.clear()\n",
    "    print(\"...done writing\")\n",
    "\n",
    "total_lines=872493862 # figured this out using \"zcat protein2ipr.dat.gz | wc -l\"\n",
    "with gzip.open(protein2ipr_filename, \"tr\") as input_file:\n",
    "    domains_to_uniprot_list_dict=collections.defaultdict(list)\n",
    "    counter=0\n",
    "    for line in input_file:\n",
    "        fields=line.split(\"\\t\")\n",
    "        domain_accession=fields[3]\n",
    "        if domain_accession[0:2]==\"PF\":\n",
    "            domains_to_uniprot_list_dict[domain_accession].append(fields[0])\n",
    "        if(counter%5000000==0 and counter>0):\n",
    "            print(\"On line #{} ({:.2f}%) now\".format(counter,100.0*counter/total_lines))\n",
    "            if(counter%25000000==0):\n",
    "                write_progress_and_clear_lists(domains_to_uniprot_list_dict)\n",
    "        counter+=1\n",
    "write_progress_and_clear_lists(domains_to_uniprot_list_dict)"
   ]
  }
 ],
 "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
}