tutorials/esmecata_method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EsMeCaTa method tutorial\n",
    "\n",
    "The idea of this tutorial is to present the method behind EsMeCaTa. So it is not a tutorial on how to use it (it will not show the command) but more on the functions behind the workflow and what is done.\n",
    "\n",
    "As an example, EsMeCaTa will be apply to the [Buchnera example file](https://github.com/AuReMe/esmecata/blob/master/test/buchnera_workflow.tsv) containing the following taxonomic affiliation:\n",
    "\n",
    "`'cellular organisms;Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Erwiniaceae;Buchnera;Buchnera aphidicola'`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EsMeCaTa proteomes\n",
    "\n",
    "### From taxonomic affiliations to taxon ID\n",
    "\n",
    "First EsMeCaTa iterates through the taxonomic affiliations to find the taxon ID (in the NCBI Taxonomy database using ete3) for each taxon name. These taxon IDs will be used to search for proteomes in [UniProt Proteomes database](https://www.uniprot.org/proteomes/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{2: 'Bacteria',\n",
       " 32199: 'Buchnera',\n",
       " 46073: 'Buchnera',\n",
       " 9: 'Buchnera aphidicola',\n",
       " 131567: 'cellular organisms',\n",
       " 91347: 'Enterobacterales',\n",
       " 1903409: 'Erwiniaceae',\n",
       " 1236: 'Gammaproteobacteria',\n",
       " 1224: 'Proteobacteria'}"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.proteomes import taxonomic_affiliation_to_taxon_id\n",
    "from ete3 import NCBITaxa\n",
    "ncbi = NCBITaxa()\n",
    "\n",
    "taxonomic_affiliation = 'cellular organisms;Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Erwiniaceae;Buchnera;Buchnera aphidicola'\n",
    "\n",
    "tax_ids_to_names, taxon_ids = taxonomic_affiliation_to_taxon_id('Buchnera', taxonomic_affiliation, ncbi)\n",
    "\n",
    "tax_ids_to_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each taxon in the taxonomic affiliations, we have a taxon ID. If no corresponding taxon IDs were found for a taxon name, it will shown 'not found'.\n",
    "\n",
    "But there is a risk that some taxon names are associated with multiple taxon IDs so we have to disambiguate the taxon IDs found."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Disambiguate taxon IDs found\n",
    "\n",
    "A taxon name can be associated with multiple taxon IDs.\n",
    "\n",
    "For example, we will use the taxonomic affiliation `Bacteria;Gammaproteobacteria;Yersinia`. The taxon Yersinia is associated with 2 taxon IDs (one mantid (taxId: 444888) and one bacterial (taxId: 629)). EsMeCaTa will not able to differentiate them.\n",
    "\n",
    "EsMeCaTa will compare all the taxons of the taxonomic affiliation (here: 2 (Bacteria) and 1236 (Gammaproteobacteria)) to the lineage associated with the two taxIDs (for bacteria Yersinia: `[1, 131567, 2, 1224, 1236, 91347, 1903411, 629]` and for the mantid one: `[1, 131567, 2759, 33154, 33208, 6072, 33213, 33317, 1206794, 88770, 6656, 197563, 197562, 6960, 50557, 85512, 7496, 33340, 33341, 6970, 7504, 7505, 267071, 444888]`).\n",
    "In this example, there is 2 matches for the bacteria one (2 and 1236) and 0 for the mantid one. So EsMeCaTa will select the taxon ID associated with the bacteria (629).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before Disambiguation:\n",
      " {'Yersinia': OrderedDict([('Bacteria', [2]), ('Gammaproteobacteria', [1236]), ('Yersinia', [629, 444888])])}\n",
      "After Disambiguation:\n",
      " {'Yersinia': OrderedDict([('Bacteria', [2]), ('Gammaproteobacteria', [1236]), ('Yersinia', [629])])}\n"
     ]
    }
   ],
   "source": [
    "from esmecata.core.proteomes import disambiguate_taxon, taxonomic_affiliation_to_taxon_id\n",
    "\n",
    "yersinia_taxonomic_affiliation = 'Bacteria;Gammaproteobacteria;Yersinia'\n",
    "\n",
    "yersinia_tax_ids_to_names, yersinia_taxon_ids = taxonomic_affiliation_to_taxon_id('Yersinia', yersinia_taxonomic_affiliation, ncbi)\n",
    "\n",
    "yersinia_json_taxonomic_affiliations = {'Yersinia': yersinia_taxon_ids}\n",
    "print('Before Disambiguation:\\n', yersinia_json_taxonomic_affiliations)\n",
    "yersinia_json_taxonomic_affiliations = disambiguate_taxon(yersinia_json_taxonomic_affiliations, ncbi)\n",
    "print('After Disambiguation:\\n', yersinia_json_taxonomic_affiliations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, Yersinia is now associated with only one taxonID.\n",
    "\n",
    "With the taxonomic affiliation of Buchnera, this will not change anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Buchnera': OrderedDict([('cellular organisms', [131567]),\n",
       "              ('Bacteria', [2]),\n",
       "              ('Proteobacteria', [1224]),\n",
       "              ('Gammaproteobacteria', [1236]),\n",
       "              ('Enterobacterales', [91347]),\n",
       "              ('Erwiniaceae', [1903409]),\n",
       "              ('Buchnera', [32199]),\n",
       "              ('Buchnera aphidicola', [9])])}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.proteomes import disambiguate_taxon\n",
    "\n",
    "json_taxonomic_affiliations = {'Buchnera': taxon_ids}\n",
    "json_taxonomic_affiliations = disambiguate_taxon(json_taxonomic_affiliations, ncbi)\n",
    "json_taxonomic_affiliations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Limit taxonomic rank used\n",
    "\n",
    "Then if a rank_limit has been given, EsMeCaTa will shrink the taxonomic affiliations to keep only the taxon equal or below the specified taxonomic rank. By default, it is not applied so there will be no changed.\n",
    "\n",
    "But as an example, we will try with a rank_limit on the `family`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Buchnera': OrderedDict([('Erwiniaceae', [1903409]),\n",
       "              ('Buchnera', [32199]),\n",
       "              ('Buchnera aphidicola', [9])])}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.proteomes import filter_rank_limit\n",
    "\n",
    "rank_limit = 'family'\n",
    "\n",
    "json_taxonomic_affiliations_limited = filter_rank_limit(json_taxonomic_affiliations, ncbi, rank_limit)\n",
    "\n",
    "json_taxonomic_affiliations_limited"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see all the ranks above the family `Erwiniaceae` (all the ranks superior to `family`) have been removed from the taxonomic affliations.\n",
    "\n",
    "This option is useful to avoid using high taxonomic ranks (containing too many organisms to make a valid estimation of the metabolism). For example, if EsMeCaTa finds only taxon IDs at the `phylum` rank, using this option with a rank_limit of `order` will remove this taxonomic affiliation from the dataset as there is no taxonomic rank equal or below the taxonomic rank `order`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search for proteomes in UniProt\n",
    "Then we use the taxon IDs and the taxonomic affiliations to search for proteomes in UniProt.\n",
    "\n",
    "EsMeCaTa searches for proteomes with a score of [BUSCO score](https://www.uniprot.org/help/assessing_proteomes) higher or equal to 80%. Furthermore, proteomes described by UniProt as [redundant](https://www.uniprot.org/help/proteome_redundancy_faq) or [excluded](https://www.uniprot.org/help/proteome_exclusion_reasons) will not be selected.\n",
    "\n",
    "Also EsMeCaTa will first search for [reference](https://www.uniprot.org/help/reference_proteome) proteomes. If no reference proteome is available (or not enough reference proteoems are available), it will then search in non-reference proteome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Proteomes for Buchnera aphidicola:\n",
      "Reference proteomes (with Busco score superior to 80):  UP000000601,UP000001806\n",
      "All proteomes (with Busco score superior to 80):  UP000000416,UP000001269,UP000006139,UP000006904,UP000019087,UP000066321,UP000093070,UP000244884,UP000271533,UP000298585,UP000298654,UP000298660,UP000298759,UP000325981,UP000502958,UP000503183,UP001163441,UP000502374,UP001163094,UP001163440\n",
      "Used proteomes:  22\n"
     ]
    }
   ],
   "source": [
    "from esmecata.core.proteomes import rest_query_proteomes\n",
    "\n",
    "busco_percentage_keep = 80\n",
    "all_proteomes = None\n",
    "\n",
    "observation_name = 'Buchnera'\n",
    "tax_name = 'Buchnera aphidicola'\n",
    "tax_id = json_taxonomic_affiliations[observation_name][tax_name][0]\n",
    "proteomes, organism_ids, data_proteomes = rest_query_proteomes(observation_name, tax_id, tax_name, busco_percentage_keep, all_proteomes)\n",
    "reference_proteomes = [proteome_data[0] for proteome_data in data_proteomes if proteome_data[4] is True and proteome_data[1] >= 80]\n",
    "non_reference_proteomes  = [proteome_data[0] for proteome_data in data_proteomes if proteome_data[4] is False and proteome_data[1] >= 80]\n",
    "if tax_name in ['Buchnera aphidicola']:\n",
    "    print('Proteomes for Buchnera aphidicola:')\n",
    "    print('Reference proteomes (with Busco score superior to 80): ', ','.join(reference_proteomes))\n",
    "    print('All proteomes (with Busco score superior to 80): ', ','.join(non_reference_proteomes))\n",
    "    print('Used proteomes: ', len(reference_proteomes) + len(non_reference_proteomes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we retrieve 2 reference proteomes with `Buchnera aphidicola`. The minimal number of proteomes required by EsMeCaTa is 5.\n",
    "So EsMeCaTa will use reference (2) and non-reference proteomes (30) for the next step.\n",
    "\n",
    "If the number of proteomes is superior to a limit (by default 100 proteomes), a subsampling procedure will be performed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subsampling procedures\n",
    "\n",
    "We will used an artifical example to show how the subsampling works.\n",
    "\n",
    "In this artificial example, we have searched for the genus Escherichia and have found 20 proteomes.\n",
    "Among these 20 proteomes:\n",
    "- 10 are associated with *Escherichia coli* (taxon ID 562).\n",
    "- 5 to *Escherichia albertii* (taxon ID 208962).\n",
    "- 4 to *Escherichia fergusonii* (taxon ID 564).\n",
    "- 1 to *Escherichia alba* (taxon ID 2562891).\n",
    "\n",
    "We have put a limit of proteomes at 10 (by default it is 100).\n",
    "\n",
    "EsMeCaTa will subsample a number close to 10 but it will try to keep the taxonomic diversity observed.\n",
    "So it will compute the percentage of contribution of each organism to the number of proteomes associated with *Escherichia*. This will give the following percentages:\n",
    "\n",
    "- *Escherichia coli*: 10 proteomes / 20 = 50%.\n",
    "- *Escherichia albertii*: 5 proteomes / 20 = 25%.\n",
    "- *Escherichia fergusonii*: 4 proteomes / 20 = 20%.\n",
    "- *Escherichia alba*: 1 proteome / 20 = 5%.\n",
    "\n",
    "But there is another rule: even if an organism contributes to less than 1 percent, it must be present in the proteomes returned by EsMeCaTa. This rule and the round applied to the number of proteomes lead to variation in the number of proteomes returned (which can be different compared to the input limit_maximal_number_proteomes).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "208962 3\n",
      "2562891 1\n",
      "562 5\n",
      "564 2\n"
     ]
    }
   ],
   "source": [
    "from esmecata.core.proteomes import subsampling_proteomes\n",
    "from collections import Counter\n",
    "\n",
    "organism_ids = {'2562891': ['UP000477739'], '208962': ['UP000292187', 'UP000193938', 'UP000650750', 'UP000407502', 'UP000003042'],\n",
    "                '562': ['UP000000625', 'UP000000558', 'UP000464341', 'UP000219757', 'UP000092491',\n",
    "                        'UP000567387', 'UP000234906', 'UP000016096', 'UP000017268', 'UP000017618'],\n",
    "                '564': ['UP000510927', 'UP000392711', 'UP000000745', 'UP000033773']}\n",
    "revert_organism_ids ={}\n",
    "for org_id, proteomes in organism_ids.items():\n",
    "    revert_organism_ids.update({proteome_id: org_id for proteome_id in proteomes})\n",
    "limit_maximal_number_proteomes = 10\n",
    "\n",
    "selected_proteomes = subsampling_proteomes(organism_ids, limit_maximal_number_proteomes, ncbi)\n",
    "selected_organisms = [revert_organism_ids[proteome] for proteome in selected_proteomes]\n",
    "\n",
    "for org_id in Counter(selected_organisms):\n",
    "    print(org_id, Counter(selected_organisms)[org_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the subsampling selects 11 proteomes among the 20:\n",
    "\n",
    "- 5 proteomes from *Escherichia coli* (taxon ID 562).\n",
    "- 3 proteomes from *Escherichia albertii* (taxon ID 208962).\n",
    "- 2 proteomes from *Escherichia fergusonii* (taxon ID 564).\n",
    "- 1 proteome from *Escherichia alba* (taxon ID 2562891)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download proteomes\n",
    "\n",
    "After these steps, EsMeCaTa will select the lowest taxonomic rank in the taxonomic affiliations associated with at least 5 proteomes.\n",
    "\n",
    "In our example:\n",
    "\n",
    "`'cellular organisms;Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Erwiniaceae;Buchnera;Buchnera aphidicola'`\n",
    "\n",
    "The lowest taxonomic rank with 5 proteomes is the species `Buchnera aphidicola` and its 32 proteomes `'['UP000000601', 'UP000001806', 'UP000000416', 'UP000001269', 'UP000006139', 'UP000006904', 'UP000019087', 'UP000066321', 'UP000502374', 'UP000503183', 'UP000093070', 'UP000244884', 'UP000271533', 'UP000298564', 'UP000298575', 'UP000298585',\n",
    " 'UP000298594', 'UP000298613', 'UP000298654', 'UP000298660', 'UP000298673', 'UP000298688', 'UP000298716',\n",
    " 'UP000298738', 'UP000298745', 'UP000298759', 'UP000298773', 'UP000298791', 'UP000298793', 'UP000325981',\n",
    " 'UP000326914', 'UP000502958']`.\n",
    "\n",
    "\n",
    "Then the selected proteomes are downloaded."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EsMeCaTa clustering\n",
    "\n",
    "### Protein clusters with MMseqs2\n",
    "\n",
    "Using these proteomes, EsMeCaTa will use [MMseqs2](https://github.com/soedinglab/MMseqs2) to cluster the proteins.\n",
    "\n",
    "First, mmseqs creates a sequence database from the proteomes with the command `createdb`.\n",
    "\n",
    "Then it performs the clustering using `cluster` with the following default options `--min-seq-id 0.3 -c 0.8`. These options search for a minimal sequence identity of 30 % and a coverage of 80%. The idea of these options is to be close to similar threshold for homology between sequences.\n",
    "\n",
    "From the result of the clustering `result2profile` will create profiles. Using these profiles mmseqs `profile2consensus` will create consensus sequences. These consensus are put into a fasta file with `convert2fasta`.\n",
    "\n",
    "A fasta file containing the representative sequences is also created with `convert2fasta`.\n",
    "\n",
    "The results of the clustering are stored into a tsv file with `createtsv`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tsv file containing the results is as follow:\n",
    "\n",
    "|         |         |         |         |\n",
    "|---------|---------|---------|---------|\n",
    "|  Protein_1 | Protein_1  | Protein_1a  | Protein_1b  |\n",
    "|  Protein_2 | Protein_2  | Protein_2a  | Protein_2b  |\n",
    "|  Protein_3 | Protein_3  |   |   |\n",
    "|  Protein_4 | Protein_4  | Protein_4a  |  |\n",
    "\n",
    "One row corresponds to one protein cluster created by mmseqs.\n",
    "\n",
    "The first column corresponds to the name of the cluster (which is the ID of the representative protein of the cluster). Then the following column contains the protein of this cluster (the first ID is the ID of the representative protein of the cluster, that is why it occurs two times in a row).\n",
    "\n",
    "Using this file, EsMeCaTa wil create the following dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_clusters = {'Protein_1a': ['Protein_1a', 'Protein_1b', 'Protein_1c'],\n",
    "                    'Protein_2a': ['Protein_2a', 'Protein_2b', 'Protein_2c'],\n",
    "                   'Protein_3a': ['Protein_3a'],\n",
    "                   'Protein_4a': ['Protein_4a', 'Protein_4b']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cluster filtering according to proteome representation\n",
    "\n",
    "These clusters will be filtered by EsMecaTa according to the representation of each proteomes associated with a taxon in the cluster.\n",
    "\n",
    "This is performed in the `esmecata.clustering.filter_protein_cluster` function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Protein_1a {'Proteome_B', 'Proteome_C', 'Proteome_A'}\n",
      "Protein_2a {'Proteome_B', 'Proteome_C', 'Proteome_A'}\n",
      "Protein_3a {'Proteome_A'}\n",
      "Protein_4a {'Proteome_B', 'Proteome_A'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'Protein_1a': 1.0,\n",
       " 'Protein_2a': 1.0,\n",
       " 'Protein_3a': 0.3333333333333333,\n",
       " 'Protein_4a': 0.6666666666666666}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "organism_prots = {'Protein_1a': 'Proteome_A', 'Protein_2a': 'Proteome_A',\n",
    "                'Protein_3a': 'Proteome_A', 'Protein_4a': 'Proteome_A',\n",
    "                 'Protein_1b': 'Proteome_B', 'Protein_2b': 'Proteome_B',\n",
    "                 'Protein_4b': 'Proteome_B', 'Protein_1c': 'Proteome_C',\n",
    "                  'Protein_2c': 'Proteome_C'}\n",
    "number_proteomes = 3\n",
    "\n",
    "computed_threshold_cluster = {}\n",
    "for rep_protein in protein_clusters:\n",
    "    proteomes_cluster = set([organism_prots[prot] for prot in protein_clusters[rep_protein]])\n",
    "    print(rep_protein, proteomes_cluster)\n",
    "    computed_threshold_cluster[rep_protein] = len(proteomes_cluster) / number_proteomes\n",
    "\n",
    "computed_threshold_cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cluster `Protein_1a` contains 3 proteins (`Protein_1a` from the proteome `Proteome_A`, `Protein_1b` from `Proteome_B` and protein `Protein_1c` from proteome `Proteome_C`). So it contains all the proteome given as input to MMseqs2.\n",
    "\n",
    "Whereas the cluster `Protein_3a` contains only 1 protein (`Protein_3a`) from 1 proteome (`Proteome_A`). So it contains only 33% of the proteome given to MMseqs2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Protein_1a', 'Protein_2a', 'Protein_4a']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reference_threshold = 0.5\n",
    "kept_clusters = []\n",
    "for rep_protein in protein_clusters:\n",
    "    if computed_threshold_cluster[rep_protein] >= reference_threshold:\n",
    "        kept_clusters.append(rep_protein)\n",
    "kept_clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then a threshold is used to filter each clustering. By default it is 0.5 to select protein clusters containing protein from most of the proteomes.\n",
    "\n",
    "So in the previous example, the cluster `Protein_3a` will be removed as it contains only a protein from one proteome.\n",
    "\n",
    "If the threshold is put to 0, it will select all the clusters associated with the proteomes of the taxon.\n",
    "\n",
    "The protein clusters kept are then written into different files and will be used during the next step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EsMeCaTa annotation\n",
    "\n",
    "There are two methods to retrieve annotation.\n",
    "A first one using UniProt annotations and a method to propagate the annotion in protein clusters.\n",
    "And a second one using eggnog-mapper to annotate consensus sequences associated with protein clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Retrieve annotation using eggnog-mapper\n",
    "\n",
    "By default, EsMeCaTa will use eggnog-mapper. With the protein clusters kept during the previous step, EsMeCaTa will use eggnog-mapper to annotate the consensus protein sequences associated with the protein clusters. \n",
    "\n",
    "It will return:\n",
    "\n",
    "- protein name\n",
    "- list of GO Terms\n",
    "- list of Enzyme Commissions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Retrieve annotation from UniProt (less accurate)\n",
    "\n",
    "Using the protein cluster kept during the previous step, EsMeCaTa will query UniProt to get the annotations of the proteins inside each cluster.\n",
    "\n",
    "It will return:\n",
    "\n",
    "- protein name\n",
    "- if the protein is reviewed or not\n",
    "- list of GO Terms\n",
    "- list of Enzyme Commissions\n",
    "- list of InterPro domains\n",
    "- list of Rhea IDs\n",
    "- gene name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Q89AL1': ['Ribosomal large subunit pseudouridine synthase B',\n",
       "  'UniProtKB reviewed (Swiss-Prot)',\n",
       "  ['GO:0000455', 'GO:0160139', 'GO:0003723'],\n",
       "  ['5.4.99.22'],\n",
       "  ['IPR042092',\n",
       "   'IPR018496',\n",
       "   'IPR020103',\n",
       "   'IPR002942',\n",
       "   'IPR006145',\n",
       "   'IPR050343',\n",
       "   'IPR036986',\n",
       "   'IPR000748',\n",
       "   'IPR020094'],\n",
       "  ['RHEA:42520'],\n",
       "  'rluB'],\n",
       " 'P57369': ['Ribosomal large subunit pseudouridine synthase B',\n",
       "  'UniProtKB reviewed (Swiss-Prot)',\n",
       "  ['GO:0001522', 'GO:0006364', 'GO:0160139', 'GO:0003723'],\n",
       "  ['5.4.99.22'],\n",
       "  ['IPR042092',\n",
       "   'IPR018496',\n",
       "   'IPR020103',\n",
       "   'IPR006145',\n",
       "   'IPR050343',\n",
       "   'IPR036986',\n",
       "   'IPR000748',\n",
       "   'IPR020094'],\n",
       "  ['RHEA:42520'],\n",
       "  'rluB']}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.annotation import query_uniprot_annotation_rest\n",
    "\n",
    "protein_queries = ['Q89AL1', 'P57369']\n",
    "\n",
    "output_dict = {}\n",
    "output_dict = query_uniprot_annotation_rest(protein_queries, output_dict)\n",
    "output_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Propagation of annotation in the cluster\n",
    "\n",
    "Using annotation present in each proteins of the cluster, EsMeCaTa will propagate the annotation according to the propagation `-p` option. This will impact only the GO Terms and EC. For the protein function and the gene name, EsMeCaTa will select the protein function/gene name with the most occurrences in the protein cluster.\n",
    "\n",
    "By default, `-p` is set to 1 meaning an annotation is kept if it is found in all the proteins of the cluster.\n",
    "\n",
    "For example this is the case for the previous cell code. `Q89AL1` and `P57369` have both the same GO Terms and EC so the cluster will be associated to these annotations.\n",
    "\n",
    "We will use another example to explain the other option. The following protein cluster contains 3 proteins with different annotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dict = {'prot_1': ['function_1', True, ['GO:0031522', 'GO:0004765'], ['7.4.2.8'], ['IPR027417'], [], 'gene_1'],\n",
    "                'prot_2': ['function_2', True, ['GO:0031522', 'GO:0005737'], ['7.4.2.8', '2.7.1.71'], ['IPR027417'], [], 'gene_1'],\n",
    "                 'prot_3': ['function_1', True, ['GO:0031522', 'GO:0005737'], ['7.4.2.8'], [], [], 'gene_3']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using the default propagation `-p 1` we keep annotations occurring in all proteins:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'prot_1': ['function_1', ['GO:0031522'], ['7.4.2.8'], 'gene_1']}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.annotation import propagate_annotation_in_cluster\n",
    "\n",
    "propagate_annotation = 1\n",
    "reference_proteins = {'prot_1': ['prot_1', 'prot_2', 'prot_3']}\n",
    "uniref_output_dict = None\n",
    "    \n",
    "protein_annotations = propagate_annotation_in_cluster(output_dict, reference_proteins, propagate_annotation, uniref_output_dict)\n",
    "\n",
    "protein_annotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using the propagation `-p 0.66` we keep annotations if they occur at least in 66 % proteins of the cluster (here 2 proteisn among the 3).\n",
    "\n",
    "In this way, we add the GO Term `GO:0005737`, which is present in prot_1 and prot_2 but not in prot_3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'prot_1': ['function_1', ['GO:0005737', 'GO:0031522'], ['7.4.2.8'], 'gene_1']}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.annotation import propagate_annotation_in_cluster\n",
    "\n",
    "propagate_annotation = 0.66\n",
    "reference_proteins = {'prot_1': ['prot_1', 'prot_2', 'prot_3']}\n",
    "uniref_output_dict = None\n",
    "\n",
    "protein_annotations = propagate_annotation_in_cluster(output_dict, reference_proteins, propagate_annotation, uniref_output_dict)\n",
    "\n",
    "protein_annotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the propagation `-p 0` will keep all the annotations present in the cluster.\n",
    "\n",
    "This adds the GO Term `GO:0004765` (only present in prot_1) and the EC `2.7.1.71` (only present in prot_2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'prot_1': ['function_1',\n",
       "  ['GO:0004765', 'GO:0005737', 'GO:0031522'],\n",
       "  ['2.7.1.71', '7.4.2.8'],\n",
       "  'gene_1']}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from esmecata.core.annotation import propagate_annotation_in_cluster\n",
    "\n",
    "propagate_annotation = 0\n",
    "reference_proteins = {'prot_1': ['prot_1', 'prot_2', 'prot_3']}\n",
    "uniref_output_dict = None\n",
    "    \n",
    "protein_annotations = propagate_annotation_in_cluster(output_dict, reference_proteins, propagate_annotation, uniref_output_dict)\n",
    "\n",
    "protein_annotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then these annotations are written into a tsv file and pathologic file (for Pathway Tools)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}