ncbi_lumberjack_og.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "ncbi_lumberjack_2.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/biothomme/Limonene/blob/master/ncbi_lumberjack_og.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oTz6BHxh85rC",
        "colab_type": "text"
      },
      "source": [
        "#NCBI Lumberjack OG\n",
        "###Uppsala, Dec 2019\n",
        "In the NCBI database genbank it can be very tiring to get your fasta sequences... If genbank was a jungle, this script would be lumberjack, enabling a great overview of the jungle, to get your data - fast and fancy.\n",
        "\n",
        "Before you start, you should ensure, to aim fasta sequences for <b>one specific</b> (but also multiple) <b>genes</b> of a <b>monophyletic taxon</b>. The taxonomic level to retrieve entries can be set to genus or species.\n",
        "\n",
        "Like OG Simpson may be more than Homer S. <b>NCBI Lumberjack OG</b> outcompetes the 'old' NCBI Lumberjack 2: It is possible to retrieve sequences for multiple outgroups.\n",
        "\n",
        "---\n",
        "---\n",
        "\n",
        "So once again: This script enables the datamining of the NCBI database for genes in monophyletic groups. Follow the blocks straight forward to gain your result. You can mine x gene entries per Genus/Species of your monophylum and download a folder containing <b>all sequences</b> in a <b>fasta-file</b> together with 2 informative files:\n",
        "- <b>\\*stat.csv:</b> includes all taxIDs (only on level genus/species) and Accession IDs for your gene within the monophylum.\n",
        "- <b>\\*data.csv:</b> corresponds to your fasta file and includes very important data of the sequences\n",
        "<br> \n",
        "<br> \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cWUhz7j9au8U",
        "colab_type": "text"
      },
      "source": [
        "#New to Colab?\n",
        "Colab is a project of Google, which allows to share Notebooks with interactive code-commands (e.g. in Python). Best is to log in with a Google account and if necessary press the button on the upper left > Open in Playground <.\n",
        "\n",
        "Afterwards please follow all paragraphs and execute the blocks with squared brackets on the upper left corner, which look like:\n",
        "\n",
        "# [ <font color=\"white\">˘</font> ]\n",
        "\n",
        " This can be performed by clicking on the field between the brackets. If the block is a form, please fill in the fields first and execute to confirm your input.\n",
        "\n",
        "<b>To run the script successful please follow each sneaky block, read carefully and only skip fields if declared as <font color=\"grey\">OPTIONAL</font> and share <font color=\"grey\">GREY</font> font color!</b>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fxSVqGIgsHtc",
        "colab_type": "text"
      },
      "source": [
        "#0 Usage\n",
        "This script works with the help of two powerful instutitions: biopython (https://biopython.org/wiki/Documentation) and through a wonderful backdoor to NCBI databases: the Entrez E-utilities service (https://www.ncbi.nlm.nih.gov/books/NBK25497/). If you use the sequence data for scientific purpose (and not only for fun ;-) ), please cite those sources properly (especially biopython requires a citation of Cock et al., 2009!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EmS3wSJxmgXn",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "#1 `ncbi_miner`"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dbx6paTIvMLC",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "\n",
        "At first the biopython (Cock et al., 2009) package needs to be installed. \n",
        "\n",
        "**Execute the block!**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9jhrXDkepk1q",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip install biopython"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "94XwSgMTfD9N"
      },
      "source": [
        "##1.1 Mandatory arguments for search\n",
        "The following field is the heart of the lumberjack. Nearly all functions are encoded within.\n",
        "\n",
        "**Execute the block!**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Gg8cPOoXm-o-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title {display-mode: 'form'}\n",
        "#@markdown The heart of Lumberjack is the function `ncbi_miner`. \n",
        "#@markdown It feeds through the NCBI jungle like the hyperactive caterpillar of a leaf miner. Do not forget to execute!\n",
        "\"\"\"\n",
        "Created on Tue Nov 26 09:29:22 2019\n",
        "\n",
        "@author: Thomsn\n",
        "\"\"\"\n",
        "\n",
        "__author__ = 'thomas m. huber'\n",
        "__email__ = ['thomas.huber@evobio.eu']\n",
        "\n",
        "def ncbi_miner(taxon_query,\n",
        "               gene_name,\n",
        "               gene_length,\n",
        "               my_mail,\n",
        "               outgroups = [],\n",
        "               coding_sequence = True,\n",
        "               exclude_query = None,\n",
        "               my_API_key = None,\n",
        "               tree_resolution = 'genus',\n",
        "               length_tolerance_percent = 20,\n",
        "               upper_tolerance = 2,\n",
        "               search_limit = 10000,\n",
        "               entries_per_genus = 1,\n",
        "               entries_per_tax = 1,\n",
        "               taxonlist_path = None,\n",
        "               random_mining = False,\n",
        "               strict_search = True):\n",
        "    import os\n",
        "    import pandas as pd\n",
        "    from Bio import SeqIO\n",
        "    from Bio import Entrez\n",
        "    from datetime import datetime\n",
        "    from urllib.error import HTTPError\n",
        "    if random_mining:\n",
        "        from random import shuffle\n",
        "    start_time = datetime.now()\n",
        "\n",
        "\n",
        "    newpath = f'{datetime.now().strftime(\"%Y_%m_%d\")}_{taxon_query}_{gene_name}'\n",
        "    if not os.path.exists(newpath):\n",
        "        os.makedirs(newpath)\n",
        "        print(f'Step 0:\\n >> Directory {newpath} was created.')\n",
        "\n",
        "\n",
        "    gl_lowend = gene_length * (100 - length_tolerance_percent)/100\n",
        "    gl_highend = gene_length * (100 + upper_tolerance * length_tolerance_percent)/100\n",
        "    \n",
        "    tree_resolution = str(tree_resolution).lower()\n",
        "\n",
        "    tax_columns = ['taxID', 'genus', 'epithet', 'entry_UIDs', 'count', 'outgroup']\n",
        "\n",
        "    if tree_resolution == 'genus':\n",
        "        taxon_level = 'Genus'\n",
        "    elif tree_resolution == 'species':\n",
        "        taxon_level = 'species'\n",
        "    else:\n",
        "        print(f'Mining data with tree resolution on the taxonomic level \\'{tree_resolution}\\' \\\n",
        "is not possible with this script. Try \\'genus\\' (default) or \\'species\\'.')\n",
        "\n",
        "    if taxonlist_path:\n",
        "        try:\n",
        "            old_taxon_list = pd.read_csv(taxonlist_path)\n",
        "        except FileNotFoundError:\n",
        "            print(f'It was not possible to find input file {taxonlist_path}, please check \\\n",
        "the path and restart \\'ncbi_miner\\'.')\n",
        "            return\n",
        "        else:\n",
        "            if all([(any(old_taxon_list.keys() == i )) for i in tax_columns]):\n",
        "                print(f'The csv-file {taxonlist_path} was loaded. Step 1 will be \\\n",
        "skipped. Follow up steps will use those taxa to search for sequences.')\n",
        "\n",
        "                Entrez.email = my_mail\n",
        "                Entrez.api_key = my_API_key\n",
        "                all_taxaIDs = list(old_taxon_list['taxID'])\n",
        "        \n",
        "                progress_criterion = len(all_taxaIDs) // 20\n",
        "                percentage_factor = 20/100\n",
        "        \n",
        "                print(f'Step 2:\\n >> Esearch for gene entries of all species in the taxon {taxon_query}. \\\n",
        "This may take time, so keep the internet connection, chill down, drink a tea. \\n The progress is ...')\n",
        "                taxon_list = pd.DataFrame(columns = tax_columns)\n",
        "                for i, taxon in enumerate(all_taxaIDs):\n",
        "                    k = i + 1\n",
        "                    if (k % progress_criterion) == 1:\n",
        "                        print(f' - {int(i // (progress_criterion*percentage_factor))} % -')\n",
        "                    else:\n",
        "                        pass\n",
        "\n",
        "                    if coding_sequence:\n",
        "                        query = f'{gene_name}[Gene Name] txid{str(taxon)}'\n",
        "                    else:\n",
        "                        query = f'{gene_name} txid{str(taxon)}'\n",
        "                    if exclude_query:\n",
        "                        exclude_string = f' NOT {exclude_query}'\n",
        "                        query = f'{query}{exclude_query}'\n",
        "                    try:\n",
        "                        with Entrez.esearch(db='nuccore', term=query, retmax=search_limit, sort='date released') as handle:\n",
        "                            gene_record = Entrez.read(handle)\n",
        "                    except HTTPError:\n",
        "                            print('New Error, but RUN Forrest RUN!')\n",
        "                    else:\n",
        "                            count = gene_record['Count']\n",
        "                            if count == '0':\n",
        "                                pass\n",
        "                            else:\n",
        "                                taxon_list.loc[taxon] = [taxon,\n",
        "                                                     old_taxon_list['genus'].iloc[i],\n",
        "                                                     old_taxon_list['epithet'].iloc[i],\n",
        "                                                     gene_record['IdList'],\n",
        "                                                     int(gene_record['Count']),\n",
        "                                                     old_taxon_list['outgroup'].iloc[i]]\n",
        "                taxon_list.to_csv(f'{newpath}/{taxon_query}_{gene_name}_stat.csv')\n",
        "                og_list = taxon_list.iloc[:,:][taxon_list.outgroup]\n",
        "                print(og_list)\n",
        "                taxon_list = taxon_list.drop(og_list.index)\n",
        "                print(taxon_list)\n",
        "                print(f' >> Done - entry database successfully established. Summary was saved as \\\n",
        "{taxon_query}_{gene_name}_stat.csv')\n",
        "\n",
        "            ### OLD ###\n",
        "\n",
        "            else:\n",
        "                print(f'The input file {taxonlist_path} does not fit with the conditions (header).\\\n",
        "Please change and restart \\'ncbi_miner\\'.')\n",
        "                return\n",
        "\n",
        "    else:\n",
        "        #### > STEP 1 < ####\n",
        "    \n",
        "        Entrez.email = my_mail\n",
        "        Entrez.api_key = my_API_key\n",
        "        print(f'{taxon_query}[orgn]')\n",
        "        try:\n",
        "            with Entrez.esearch(db=\"taxonomy\", term=f'{taxon_query}[orgn]', retmax=search_limit) as handle:\n",
        "                record = Entrez.read(handle)\n",
        "        except HTTPError:\n",
        "            return str('Database error, try later...')\n",
        "        else:\n",
        "            all_taxaIDs = record['IdList']\n",
        "              \n",
        "                #try:\n",
        "                 #       with Entrez.esearch(db=\"taxonomy\", term=f'{taxon_query}[orgn]', retmax=search_limit) as handle:\n",
        "                  #          record = Entrez.read(handle)\n",
        "                #except HTTPError:\n",
        "                 #       return str('Database error, try later...')\n",
        "                #else:\n",
        "                 #       outgroup_taxaID = pd.Series([record['IdList']])\n",
        "\n",
        "            print('Step 1:\\n >> Done - taxon database successfully established.')\n",
        "        \n",
        "    \n",
        "        #### > STEP 2 < ####\n",
        "    \n",
        "        Entrez.email = my_mail\n",
        "        Entrez.api_key = my_API_key\n",
        "        \n",
        "        progress_criterion = len(all_taxaIDs) // 20\n",
        "        percentage_factor = 20/100\n",
        "        print(f'Step 2:\\n >> Esearch for gene entries of all species in the taxon {taxon_query}. \\\n",
        "This may take time, so keep the internet connection, chill down, drink a tea. \\n The progress is ...')\n",
        "        taxon_list = pd.DataFrame(columns = tax_columns)\n",
        "        for i, taxon in enumerate(all_taxaIDs, 1):\n",
        "            if (i % progress_criterion) == 1:\n",
        "                print(f' - {int(i // (progress_criterion*percentage_factor))} % -')\n",
        "\n",
        "            try:\n",
        "                with Entrez.esummary(db=\"taxonomy\", id=taxon, retmax=search_limit) as handle:\n",
        "                    record = Entrez.read(handle)[0]\n",
        "            except IndexError:\n",
        "                print('New Error. Nemas Problemas!')\n",
        "                break\n",
        "            except HTTPError:\n",
        "                print('New Error, but bro, stay seated, we skip it and keep searching...')\n",
        "            else:\n",
        "                if record['Rank'] == 'species' and record['Genus'] != '':\n",
        "                    try:\n",
        "                        if coding_sequence:\n",
        "                                query = f'{gene_name}[Gene Name]+txid{str(taxon)}'\n",
        "                        else:\n",
        "                                query = f'{gene_name}+txid{str(taxon)}'\n",
        "                        if exclude_query:\n",
        "                                exclude_string = f'+NOT+{exclude_query}'\n",
        "                                query = f'{query}{exclude_query}'\n",
        "                        with Entrez.esearch(db='nuccore', term=query, retmax=search_limit, sort='pub+date') as handle:\n",
        "                            gene_record = Entrez.read(handle)\n",
        "                    except HTTPError:\n",
        "                        print('New Error! But nothing big to worry about.')\n",
        "                    else:\n",
        "                        count = gene_record['Count']\n",
        "                        if count == '0':\n",
        "                            pass\n",
        "                        else:\n",
        "                            taxon_list.loc[taxon] = [taxon,\n",
        "                                                 record['Genus'],\n",
        "                                                 record['Species'],\n",
        "                                                 gene_record['IdList'],\n",
        "                                                 int(gene_record['Count']),\n",
        "                                                 'False']\n",
        "                else:\n",
        "                    pass\n",
        "        og_list = pd.DataFrame(columns = tax_columns)\n",
        "        if len(outgroups) > 0:\n",
        "            for og in outgroups:\n",
        "                taxon = og\n",
        "                if coding_sequence:\n",
        "                    query = f'{gene_name}[Gene Name]+{og}[Organism]'\n",
        "                else:\n",
        "                    query = f'{gene_name}+{og}[Organism]'\n",
        "                if exclude_query:\n",
        "                    exclude_string = f'+NOT+{exclude_query}'\n",
        "                    query = f'{query}{exclude_query}'\n",
        "                try:\n",
        "                    with Entrez.esearch(db='nuccore', term=query, retmax=search_limit, sort='pub+date') as handle:\n",
        "                          gene_record = Entrez.read(handle)\n",
        "                except HTTPError:\n",
        "                    print('New Error! But nothing big to worry about.')\n",
        "                else:\n",
        "                    count = gene_record['Count']\n",
        "                    if count == '0':\n",
        "                        print('There are no entries of the gene for outgroup {og}.')\n",
        "                        pass\n",
        "                    else:\n",
        "                        og_list.loc[len(og_list)] = [taxon,\n",
        "                                            f'{taxon}_gen',\n",
        "                                            f'{taxon}_sp',\n",
        "                                            gene_record['IdList'],\n",
        "                                            int(gene_record['Count']),\n",
        "                                            'True']\n",
        "\n",
        "        else:\n",
        "            print('Warning: You have no outgroups!')\n",
        "        #for taxon in enumerate(all_taxaIDs, 1):\n",
        "            taxon_list.to_csv(f'{newpath}/{taxon_query}_{gene_name}_stat_ingroup.csv')\n",
        "        print(f' >> Done - entry database successfully established. Summary (without outgroups) was saved as \\\n",
        "{taxon_query}_{gene_name}_stat_ingroup.csv')\n",
        "        \n",
        "    \n",
        "    #### > STEP 3 < ####\n",
        "\n",
        "    print(f'Step 3:\\n >> Efetch for each species of the given genera with the most entries \\\n",
        "for the given gene. This may take some time as well, I would answer some mails in the meantime ;-) \\n The progress is ...')\n",
        "    data_list = pd.DataFrame(columns = ['taxID',\n",
        "                                        'accession',\n",
        "                                        'length',\n",
        "                                        'date',\n",
        "                                        'organism',\n",
        "                                        'reference',\n",
        "                                        'gene_information',\n",
        "                                        'sampling_locality',\n",
        "                                        'outgroup',\n",
        "                                        'genus',\n",
        "                                        'epithet'])\n",
        "    ### NEW ###\n",
        "    full_taxon_list = pd.concat([taxon_list, og_list])\n",
        "    if taxon_level == 'Genus':\n",
        "        genera = full_taxon_list['genus'].unique()\n",
        "    else:\n",
        "        genera = full_taxon_list['taxID'].unique()\n",
        "    ### OLD ###\n",
        "\n",
        "    progress_criterion = -( -len(genera) // 10)\n",
        "    percentage_factor = 10/100\n",
        "    \n",
        "    \n",
        "    for j, genus_ID in enumerate(genera):\n",
        "        if (j % progress_criterion) == 1:\n",
        "            print(f' - {int(j // (progress_criterion*percentage_factor))} % -')\n",
        "        else:\n",
        "            pass\n",
        "        if taxon_level == 'Genus':\n",
        "          all_species = full_taxon_list[full_taxon_list['genus'] == genus_ID]\n",
        "        else:\n",
        "          all_species = full_taxon_list[full_taxon_list['taxID'] == genus_ID]\n",
        "        all_species = all_species.sort_values(by=['count'], ascending=False)\n",
        "        if len(list(all_species['entry_UIDs'])) < entries_per_genus:\n",
        "            epg = len(list(all_species['entry_UIDs']))\n",
        "        else:\n",
        "            epg = entries_per_genus\n",
        "        for entry in list(range(epg)):\n",
        "        ### NEW ###\n",
        "            entry_list = list(all_species['entry_UIDs'])[entry]\n",
        "            if random_mining:\n",
        "                shuffle(entry_list)\n",
        "        ### OLD ###\n",
        "            Entrez.email = my_mail\n",
        "            Entrez.api_key = my_API_key\n",
        "            entry_counter = 0\n",
        "            for i, acc_id in enumerate(entry_list):\n",
        "                try:\n",
        "                    with Entrez.efetch(db=\"nuccore\", id=acc_id, retmode='text', rettype=\"gb\") as handle:\n",
        "                        record = SeqIO.read(handle, \"genbank\")\n",
        "                except HTTPError:\n",
        "                    print('New Error, but chill down, everything is soft')\n",
        "                else:\n",
        "                    if gl_lowend < len(record) < gl_highend:\n",
        "                        features = record.features\n",
        "                        filterframe = [x.type == 'gene' for x in features]\n",
        "                        if any(filterframe) == True:\n",
        "                            keyword = 'gene'\n",
        "                        else:\n",
        "                            filterframe = [x.type != 'source' for x in features]\n",
        "                            keyword = 'product'\n",
        "                        gene_features = [x for i, x in enumerate(features) if filterframe[i]==True]\n",
        "                        try:\n",
        "                            gene_info = [(x.qualifiers.get(keyword))[0] for x in gene_features]\n",
        "                        except TypeError:\n",
        "                            gene_info = ['']\n",
        "                        if strict_search:\n",
        "                            if (len(gene_info) == 1 and gene_name.lower() in gene_info[0].lower()):\n",
        "                                sample_locality = features[0].qualifiers.get('country')\n",
        "                                try:\n",
        "                                    referenza = record.annotations['references'][0]\n",
        "                                except TypeError:\n",
        "                                    referenza = ['']\n",
        "                                selection = full_taxon_list[full_taxon_list['taxID'] == genus_ID]\n",
        "                                data_list.loc[f'{j+1}_{record.id}'] = [genus_ID,\n",
        "                                                    record.id,\n",
        "                                                    len(record),\n",
        "                                                    record.annotations['date'],\n",
        "                                                    record.annotations['organism'],\n",
        "                                                    referenza,\n",
        "                                                    gene_info,\n",
        "                                                    sample_locality,\n",
        "                                                    selection['outgroup'].iat[0],\n",
        "                                                    selection['genus'].iat[0],\n",
        "                                                    selection['epithet'].iat[0]]\n",
        "                                with open(f'{newpath}/{taxon_query}_{gene_name}.fasta', 'a') as finalfasta:\n",
        "                                    rawseq = str(record.seq)\n",
        "                                    if taxon_level == 'species':\n",
        "                                        fasta_head = str(record.annotations['organism'])\n",
        "                                        fasta_head = fasta_head.replace(' ', '_').replace('.', '').replace('-', '_')\n",
        "                                    else:\n",
        "                                        fasta_head = genus_ID\n",
        "                                    if entries_per_tax == 1 and entries_per_genus == 1:\n",
        "                                        finalfasta.write(f'>{fasta_head}\\n')\n",
        "                                    else:\n",
        "                                        finalfasta.write(f'>{fasta_head}_{record.id}\\n')\n",
        "                                    finalfasta.write(f'{rawseq}\\n\\n')\n",
        "                        ### NEW ###\n",
        "                                entry_counter += 1\n",
        "                                if entry_counter == entries_per_tax:\n",
        "                                    break\n",
        "                                else:\n",
        "                                    pass\n",
        "                            else: \n",
        "                                pass\n",
        "                        else:\n",
        "                            sample_locality = features[0].qualifiers.get('country')\n",
        "                            try:\n",
        "                                referenza = record.annotations['references'][0]\n",
        "                            except TypeError:\n",
        "                                referenza = ['']\n",
        "                            selection = full_taxon_list[full_taxon_list['taxID'] == genus_ID]\n",
        "                            data_list.loc[len(data_list)] = [genus_ID,\n",
        "                                                record.id,\n",
        "                                                len(record),\n",
        "                                                record.annotations['date'],\n",
        "                                                record.annotations['organism'],\n",
        "                                                referenza,\n",
        "                                                gene_info,\n",
        "                                                sample_locality,\n",
        "                                                selection['outgroup'].iat[0],\n",
        "                                                selection['genus'].iat[0],\n",
        "                                                selection['epithet'].iat[0]]\n",
        "                            with open(f'{newpath}/{taxon_query}_{gene_name}.fasta', 'a') as finalfasta:\n",
        "                                rawseq = str(record.seq)\n",
        "                                if taxon_level == 'species':\n",
        "                                    fasta_head = str(record.annotations['organism'])\n",
        "                                    fasta_head = fasta_head.replace(' ', '_').replace('.', '').replace('-', '_')\n",
        "                                else:\n",
        "                                    fasta_head = genus_ID\n",
        "                                if entries_per_tax == 1 and entries_per_genus == 1:\n",
        "                                    finalfasta.write(f'>{fasta_head}\\n')\n",
        "                                else:\n",
        "                                    finalfasta.write(f'>{fasta_head}_{record.id}\\n')\n",
        "                                finalfasta.write(f'{rawseq}\\n\\n')\n",
        "                    ### NEW ###\n",
        "                            entry_counter += 1\n",
        "                            if entry_counter == entries_per_tax:\n",
        "                                break\n",
        "                            else:\n",
        "                                pass\n",
        "                ### OLD ###\n",
        "                    else:\n",
        "                        pass\n",
        "    print('Updating the outgroup data...')\n",
        "    out_selection = data_list[data_list['outgroup'] == 'True']\n",
        "    for i, acc in enumerate(out_selection['accession']):\n",
        "        org = out_selection.loc[out_selection['accession'] == acc, ['organism']]\n",
        "        org = ''.join(str(org.values[0])).replace(' ','_').replace(\"'\",'')\n",
        "        que = f'{org}'\n",
        "        print(que)\n",
        "        try:\n",
        "            with Entrez.esearch(db=\"taxonomy\", term=que, retmax=search_limit) as handle:\n",
        "                record = Entrez.read(handle)\n",
        "        except HTTPError:\n",
        "            print(f'Error, did not find {org}.')\n",
        "        else:\n",
        "            print(record)\n",
        "            txid = record['IdList'][0]\n",
        "            print(txid)\n",
        "            try:\n",
        "                with Entrez.esummary(db=\"taxonomy\", id=txid, retmax=search_limit) as handle:\n",
        "                    spef_record = Entrez.read(handle)[0]\n",
        "            except IndexError:\n",
        "                print('New Error. Nemas Problemas!')\n",
        "                break\n",
        "            except HTTPError:\n",
        "                print('Errore furore, no problemo spaghetto!')\n",
        "            else:\n",
        "                try:\n",
        "                    if coding_sequence:\n",
        "                            query = f'{gene_name}[Gene Name]+txid{str(txid)}'\n",
        "                    else:\n",
        "                            query = f'{gene_name}+txid{str(txid)}'\n",
        "                    if exclude_query:\n",
        "                            exclude_string = f'+NOT+{exclude_query}'\n",
        "                            query = f'{query}{exclude_query}'\n",
        "                    with Entrez.esearch(db='nuccore', term=query, retmax=search_limit, sort='pub+date') as handle:\n",
        "                          gene_record = Entrez.read(handle)\n",
        "                except HTTPError:\n",
        "                    print('New Error! But nothing big to worry about.')\n",
        "                else:\n",
        "                    count = gene_record['Count']\n",
        "                    if count == '0':\n",
        "                        print('Ulala, there was an error, check the outgroups in the stat/data-csv.files afterwards.')\n",
        "                        pass\n",
        "                    else:\n",
        "                        og_gen = spef_record['Genus']\n",
        "                        og_sp = spef_record['Species']\n",
        "                        data_list.at[data_list['accession'] == acc, ['taxID']] = txid\n",
        "                        data_list.at[data_list['accession'] == acc, ['genus']] = og_gen\n",
        "                        data_list.at[data_list['accession'] == acc, ['epithet']] = og_sp\n",
        "                        taxon_list.loc[txid] = [txid,\n",
        "                                              og_gen,\n",
        "                                              og_sp,\n",
        "                                              gene_record['IdList'],\n",
        "                                              int(gene_record['Count']),\n",
        "                                              'True']\n",
        "\n",
        "    data_list.to_csv(f'{newpath}/{taxon_query}_{gene_name}_data.csv')\n",
        "    taxon_list.to_csv(f'{newpath}/{taxon_query}_{gene_name}_stat.csv')\n",
        "    print(f' >> Done - most recent fasta sequences were collected and successfully concatenated. \\\n",
        "It was saved in the file {taxon_query}_{gene_name}.fasta and is proove for nexus conversion. \\\n",
        "Summary of used sequences was saved as {taxon_query}_{gene_name}_data.csv. In addition the \\\n",
        "outgroup was added to the stat-file and saved as: {taxon_query}_{gene_name}_stat.csv.')\n",
        "    stop_time = datetime.now()\n",
        "    process_time = stop_time - start_time\n",
        "    print(f' >> NCBI mining finished. It took {process_time.seconds // 60} min, \\\n",
        "{process_time.seconds % 60} sec. The files are stored in the directory {newpath}.')\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lpyNsbnYmkfB",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "#2 Parameters\n",
        "NCBI Lumberjack uses 4 mandatory parameters <b>(see 2.1)</b> and several optional parameters <b>(see 2.2)</b>. You need to execute all blocks of both paragraphs!\n",
        "\n",
        "<font color=\"grey\"><b>OPTIONAL:</b> If you have already performed a NCBI Deforestation with the Lumberjack you can upload a \\*data.csv or \\*stat.csv file for your analysis in <b>paragraph 2.3<b>.</font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eJSMMe8jmqox",
        "colab_type": "text"
      },
      "source": [
        "##2.1 Mandatory arguments for search\n",
        "Please enter your monophylum as query and run the block."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-xKON2cEpWP9",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "#@markdown Name of monophyletic taxon:\n",
        "taxon_query = 'Bufo' #@param {type:\"string\"}\n",
        "#@markdown Name of gene:\n",
        "gene_name = 'coi' #@param {type:\"string\"}\n",
        "#@markdown Estimate of gene length (check at NCBI, pubmed, ...):\n",
        "gene_length = 600 #@param {type:\"number\"}\n",
        "#@markdown Enter your mail adress (mandatory for NCBI search):\n",
        "my_mail = 'antilope.booty@evo.com' #@param {type:\"string\"}\n",
        "#@markdown <br>**Execute the block!**"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rpcmQWallgIg",
        "colab_type": "text"
      },
      "source": [
        "##2.2 Optional arguments for search\n",
        "The following parameters can be adjusted. Please run the block afterwards."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z--W6bb0zgLQ",
        "colab_type": "text"
      },
      "source": [
        "###2.2.1 Do you want an outgroup?\n",
        "Please type taxon names (lowest level: genus) in quotation marks, within square brackets; e.g. ['Soldanella', 'Primula', 'Gentiana'].\n",
        "If not used, leave empty square brackets ([]) or type just `None`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "91ePz-Rbo_F4",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "outgroups = ['Rana'] #@param {type:\"raw\"}\n",
        "#@markdown **Execute the block!**"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "b5KNW_-27fjm"
      },
      "source": [
        "###2.2.2 Is your gene_query for a protein coding sequence (i.e. gene)?\n",
        "This is a very important question. For instance 'COI' encodes the first subunit of cytochrome oxidase. RNA coding sequences like '28S' ar here <b>no coding sequence</b>, so please untick the box for queries like this. Otherwise the fasta-file will be empty :-P\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CJ35I-yx7Prv",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "coding_sequence = True #@param {type:\"boolean\"}\n",
        "#@markdown **Execute the block!**"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SWiWkzOI0eTf",
        "colab_type": "text"
      },
      "source": [
        "###2.2.3 Do you want to search faster?\n",
        "\n",
        "Enter your NCBI API key (increases search pace by factor ~ 3). Read more here : https://ncbiinsights.ncbi.nlm.nih.gov/2017/11/02/new-api-keys-for-the-e-utilities/\n",
        "\n",
        "(!) You will need to set it within quotation marks (e.g. 'AZNE192930N8D...NDJE9D0').\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8NSSu0YZzM7I",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "my_API_key =  None#@param {type:\"raw\"}\n",
        "if my_API_key:\n",
        "    my_API_key = my_API_key.replace(' ', '')\n",
        "#@markdown **Execute the block!**<br>\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8eaPMCXa10yJ",
        "colab_type": "text"
      },
      "source": [
        "###2.2.4 Do you want entries per genera or species of your monophylum?\n",
        "Which resolution should your tree have (default: genus)? i.e. taxonomic level of external branches."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KQoryYH41a3w",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "tree_resolution = 'species' #@param ['genus', 'species']\n",
        "#@markdown **Execute the block!**\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bifY9z3R21Kf",
        "colab_type": "text"
      },
      "source": [
        "###2.2.5 Do you want to only retrieve proved entries for your gene?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rj-YkyRe2PeV",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "#@markdown This is the most important filter and recommended, especially for large monophyla! Inactivate if your \\*data.csv file contains much less taxa than the \\*stat.csv file.\n",
        "strict_search = True #@param {type:\"boolean\"}\n",
        "#@markdown **Execute the block!**\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "gBHU7fAg5aDa"
      },
      "source": [
        "###2.2.6 How many entries do you want to retrieve?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "cellView": "form",
        "id": "XDL5AMSB5PIf",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "#@markdown How many species per genus do you want to retrieve at max? The species within a genus will be ranked after count of entries for the given gene. This only works if `tree_resolution` is genus.\n",
        "entries_per_genus = 1 #@param {type:\"slider\", min:1, max:20, step:1}\n",
        "#@markdown How many entries do you want to retrieve per species/taxon at max?\n",
        "entries_per_tax = 1 #@param {type:\"slider\", min:1, max:20, step:1}\n",
        "#@markdown **Execute the block!**\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jv69QRL4574G",
        "colab_type": "text"
      },
      "source": [
        "### 2.2.7 Other parameters"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7k6D-CJuAu-E",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "#@markdown Usually the results will be sorted by increasing age, so the most recent entries should be retrieved. Do you want to change it to random order?\n",
        "random_mining = False #@param {type:\"boolean\"}\n",
        "#@markdown Filter: Tolerance (in percent) of estimated gene length to examine the retrieved sequences setting lower limit (default: 50 % tolerance). This filter is weak and not needed if you use `strict_search`. But you can play with it.\n",
        "length_tolerance_percent = 100 #@param {type:\"slider\", min:0, max:100, step:1}\n",
        "#@markdown Filteradjustment: Scaling the upper limit for tolerance of `length_tolerance_percent`.\n",
        "upper_tolerance = 100 #@param {type:\"number\"}\n",
        "#@markdown Enter the limit of retrieved results per search (default: 10000), super unimportant - do not change.\n",
        "search_limit = 10000 #@param {type:\"slider\", min:1, max:10000, step:1}\n",
        "\n",
        "taxonlist_path = None\n",
        "#@markdown **Execute the block!**\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5IF2_ttS6RiR",
        "colab_type": "text"
      },
      "source": [
        "Check if you have executed all blocks properly after entering the input. \n",
        "\n",
        "If you do not want to upload a csv file as query, you should jump now to paragraph 3 and start your datamining. Otherwise continue with 2.3 (**OPTIONAL**)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T9o6IdNeFFYc",
        "colab_type": "text"
      },
      "source": [
        "##<font color=\"grey\">2.3 Upload of accession list as csv</font>\n",
        "<font color=\"grey\"><b>OPTIONAL:</b> Do you already have a set of selected accession IDs as a result of a previous use of lumberjack (i.e. deforestation) and you want to use it for a new search (e.g. different genes)? Then you should upload the \\*stat.csv (= unfiltered) or \\*data.csv (= filtered) file here:</font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m7ZRW1-oFEHo",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title <font color=\"grey\">--> Upload here! <-- </font>{display-mode: 'form'}\n",
        "#@markdown <font color=\"grey\">First run the block and then select the file, you want to use. It will be uploaded. Once you done that, you need to rerun the paragraphs 2.1 and 2.2 to perform a new neutral analysis without this file as query. </font>\n",
        "\n",
        "from google.colab import files\n",
        "uploaded = files.upload() \n",
        "\n",
        "for fn in uploaded.keys():\n",
        "  taxonlist_path = fn\n",
        "  print(f'Upload successful, {fn} will be used in the following ncbi_miner session')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Cptmqa2FlQTL",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "#3 Let's go!\n",
        "It is getting serious. The following block is the query command: it will start your datamining. Once more, check if you have the right imput arguments und start it. This will take a lot of time. Sometimes more than one hour for large monophyla. Important is to keep the connection to web, but luckily the analysis runs on a fancy cloud outside your window, so your computer will not be bothered. The datamining runs in 3 steps and will talk to you:\n",
        "\n",
        "\n",
        "1.   Getting a taxon database within the monophylum\n",
        "2.   Getting a gene sequences database for your taxa\n",
        "3.   Collecting filtered sequences for your taxa\n",
        "\n",
        "If you want to read more about it, check out the paragraph beyond the execution block.\n",
        "\n",
        "Warning: Sometimes there are problems in the connection, check your input, wait some time and it will be fine...\n",
        "\n",
        "**And now: execute the block!**\n",
        "\n",
        "<font color=\"grey\">P.S. The sometimes weird percentages during the analysis are coded like this on purpose, to keep you entertained... ;-)</font>\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2mGkjHEioV05",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "#@markdown On the left you can see how the function is called. Run it, to conduct the deforestation of the NCBI jungle...\n",
        "ncbi_miner(taxon_query = taxon_query,\n",
        "               gene_name = gene_name,\n",
        "               gene_length = gene_length,\n",
        "               my_mail = my_mail,\n",
        "               outgroups = outgroups,\n",
        "               coding_sequence = coding_sequence,\n",
        "               my_API_key = my_API_key, \n",
        "               tree_resolution = tree_resolution,\n",
        "               length_tolerance_percent = length_tolerance_percent,\n",
        "               upper_tolerance = upper_tolerance,\n",
        "               search_limit = search_limit,\n",
        "               entries_per_genus = entries_per_genus,\n",
        "               entries_per_tax = entries_per_tax,\n",
        "               taxonlist_path = taxonlist_path,\n",
        "               random_mining = random_mining,\n",
        "               strict_search = strict_search)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KrE40XWxrl9g",
        "colab_type": "text"
      },
      "source": [
        "After the Lumberjack terminated download your search of perform a ReSearch (optional)...\n",
        "\n",
        "###Background: mechanics of the Lumberjack OG:\n",
        "In Step 1 your `taxon_query` (+ outgroup, but for the latter it is more complicated) is searched against the NCBI taxonomy database (esearch) and all taxa within the monophylum will be collected. \n",
        "\n",
        "Afterwards Step 2 begins with single NCBI taxonomy database requests (esummary) and all taxa are filtered for species or genus level (parameter: `tree_resolution`). This wonderful preselection will be run for the `gene_name` against the NCBI nucleotide database (i.e. genbank; esearch) and all entries are collected (sorted from latest to oldest or random; parameter: `random_mining`) and all together is saved as the \\*stat.csv file. \n",
        "\n",
        "In the following Step 3, a maximum specified number of \n",
        "\n",
        "*   species per genus (parameter: `entries_per_genus`; works only if `tree_resolution` is genus)\n",
        "*   entries per species (parameter: `entries_per_tax`)\n",
        "\n",
        "will be retrieved from the NCBI nucleotide database (efetch). Therefore some filters (length, ...) are used.\n",
        "\n",
        "In the end the \\*data.csv file will be produced, all sequences concatenated in a single fasta-file and all 3 files are ready to download in a ZIP folder (see 5.)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n0DLOibEjy3w",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "#<font color=\"grey\">4 ReSearch with different gene but same accessions</font>\n",
        "<font color=\"grey\"><b>OPTIONAL:</b> This block uses the data which was previously produced on the server to conduct a new search for a different gene. It is like you search for all taxa which had an entry in the preceding (first!) search you did with the last taxon, you searched for.\n",
        "\n",
        "Note: If you want to upload a file from your computer, change `gene_name`and `gene_length` in 2.1, upload the file in 2.3, run all blocks in the paragraphs 2.1 - 2.3 and finally `ncbi_miner` in 3.</font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "UKLpTN7WwKfz"
      },
      "source": [
        "##<font color=\"grey\">4.1 Mandatory arguments for ReSearch</font>\n",
        "<font color=\"grey\">Please enter your new gene name as query type in the date of your preceding search and run the next two blocks.</font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u_Zg6ODdjyCU",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "#@markdown <font color=\"grey\">What is the name of your new gene?</font>\n",
        "new_gene_name = '28s' #@param {type:\"string\"}\n",
        "#@markdown <font color=\"grey\">Estimate the length of your new gene:</font>\n",
        "new_gene_length = 800 #@param {type:\"number\"}\n",
        "#@markdown <font color=\"grey\">When did you perform the search for accessions for the old gene?</font>\n",
        "date_of_old_search = \"2020-01-06\"#@param {type:\"date\"}\n",
        "date_underscore = date_of_old_search.replace('-', '_')\n",
        "coding_sequence = False #@param {type:\"boolean\"}\n",
        "\n",
        "old_file = f'{date_underscore}_{taxon_query}_{gene_name}/{taxon_query}_{gene_name}_stat.csv'"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "fAuQW_S5ww34"
      },
      "source": [
        "##<font color=\"grey\">4.2 Run the ReSearch</font>\n",
        "<font color=\"grey\">Check if you adjusted the parameters in 4.1 right and execute the next block. \n",
        "\n",
        "Note: It will not overwrite the files produced by the search before!</font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ff0X95nNmvmz",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#{display-mode: 'form'}\n",
        "#@markdown <font color=\"grey\"> Run the search with new gene - old accessions (same parameters).</font>\n",
        "ncbi_miner(taxon_query = taxon_query,\n",
        "               gene_name = new_gene_name,\n",
        "               gene_length = new_gene_length,\n",
        "               my_mail = my_mail,\n",
        "               coding_sequence = coding_sequence,\n",
        "               my_API_key = my_API_key,\n",
        "               tree_resolution = tree_resolution,\n",
        "               length_tolerance_percent = length_tolerance_percent,\n",
        "               upper_tolerance = upper_tolerance,\n",
        "               search_limit = search_limit,\n",
        "               entries_per_genus = entries_per_genus,\n",
        "               entries_per_tax = entries_per_tax,\n",
        "               taxonlist_path = old_file,\n",
        "               random_mining = random_mining)\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IIqrAl76mPsE",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "#5 Download\n",
        "After producing your data, you can download it in a zip file containing fasta file (sequences) and the corresponding data (\\*data.csv) as well as information about the gene entries of your monophylum in genbank (\\*stat.csv)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nonflVdjs_Bk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#@title 5.1 Zip file for (Re)Search\\* {display-mode: 'form'}\n",
        "from datetime import datetime\n",
        "import shutil\n",
        "from google.colab import files\n",
        "from os import listdir\n",
        "#@markdown When did you perform the (re)search you want to download?\n",
        "date_of_search = \"2020-01-07\"#@param {type:\"date\"}\n",
        "date_uscore = date_of_search.replace('-', '_')\n",
        "#@markdown Did You perform a search (paragraph 3) or  ReSearch (paragraph 4)?\n",
        "option = 'search' #@param ['search', 'research']\n",
        "\n",
        "if option == 'search':\n",
        "  zip_gene = gene_name\n",
        "else:\n",
        "  zip_gene = new_gene_name\n",
        "directory = f'{date_uscore}_{taxon_query}_{zip_gene}'\n",
        "\n",
        "folder = shutil.make_archive(directory, 'zip', directory)\n",
        "\n",
        "\n",
        "files.download(folder)\n",
        "\n",
        "\n",
        "#@markdown <br><br>The data you produced will be stored into a ZIP file and downloaded. <br> <br> \\*) Search and ReSearch store their files in different folders. Those would need to be downloaded subsequently.\n",
        "#@markdown <br><br><b>Execute this block for download!</b>\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G_BsxlXzuqvE",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "---\n",
        "#References\n",
        "<i>Cock PA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, Friedberg I, Hamelryck T, Kauff F, Wilczynski B and de Hoon MJL (2009):</i> Biopython - freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25, 1422-1423"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xkQoOeLJtEQB",
        "colab_type": "text"
      },
      "source": [
        "---\n",
        "<font color=\"grey\">Thank you for using the NCBI lumberjack, have fun with the sequences!<br>\n",
        "`Questions: thomas.huber{ett}evobio.eu`</font>"
      ]
    }
  ]
}