{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using `ncbi.datasets` library to download and parse genome annotation data in GFF3 format\n", "\n", "The objective of this notebook is to use the `ncbi.datasets` python library to demonstrate how to download and parse genome annotation data in gff3 format from NCBI Datasets.\n", "\n", "There are two major types of genome data available from NCBI Datasets:\n", "\n", "1. Genome datasets, which include genome, transcript, and protein sequences, genome annotation data in gff3 format, an assembly data report (genome assembly and annotation metadata) and a sequence report (a list of sequences that comprise the genome assembly). \n", "\n", "2. Genome summaries, which are sets of key metadata that describe the genome datasets \n", "\n", "As an example, we will download gff3 files for five _Lactobacillus_ genome assemblies, then parse those gff3 files to get information about crispr gene order in the corresponding genome assemblies.\n", "\n", "The notebook can be broken down into the following four major tasks:\n", "1. Query NCBI for _Lactobacillus_ assemblies and download the genome summaries for all available genome assemblies\n", "2. Parse the genome summaries to make a list of genome assemblies annotated in 2020\n", "3. Download sequence and annotation (gff3 and protein fasta) for five of those assemblies\n", "4. Parse the gff3 files to look at crispr gene order in those assemblies" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "After importing the various python modules we will need, create the assembly and download API instances." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from collections import defaultdict, Counter\n", "from datetime import datetime\n", "import json\n", "import os\n", "from textwrap import dedent\n", "import zipfile\n", "\n", "import gffutils\n", "import matplotlib.pyplot as plt\n", "import ncbi.datasets\n", "import pandas as pd\n", "from pyfaidx import Fasta\n", "\n", "plt.style.use('ggplot')\n", "\n", "genome_api = ncbi.datasets.GenomeApi(ncbi.datasets.ApiClient())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Group genome assemblies based on annotation date\n", "Let's look at RefSeq *lactobacillus* genome assemblies, and see how we can group them based on annotation date. Annotation date, when available, is described in the genome summary.\n", "\n", "This section will include the following two major tasks:\n", "1. Query NCBI for Lactobacillus assemblies and download the genome summaries for all available genome assemblies\n", "2. Make a list of genome assemblies annotated in 2020\n", "\n", "Using the `ncbi.datasets` library, let's first get the count of available RefSeq *lactobacillus* assemblies. By setting `limit='none'`, we tell the API to return a genome summary with only the count of genome assemblies for the specified NCBI Taxonomy ID." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of RefSeq lactobacillus genome assemblies: 777\n" ] } ], "source": [ "taxid = 1578 ## This is the NCBI Taxonomy ID for lactobacillus\n", "genome_summary = genome_api.assembly_descriptors_by_taxon(\n", " taxon=taxid,\n", " limit='none',\n", " filters_refseq_only=True)\n", "\n", "print(f\"Number of RefSeq lactobacillus genome assemblies: {genome_summary.total_count}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's download the genome summaries for these genome assemblies. To get the metadata for all of the genomes in scope, set `limit='all'`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Genome summaries for 777 genome assemblies downloaded and saved to genome_summary.\n", "As an example, here's the first genome summary for the first genome assembly in the list:\n", " {'assembly': {'annotation_metadata': {'file': [{'estimated_size': '111943',\n", " 'type': 'GENOME_GFF'},\n", " {'estimated_size': '1187873',\n", " 'type': 'GENOME_GBFF'},\n", " {'estimated_size': '301643',\n", " 'type': 'PROT_FASTA'},\n", " {'estimated_size': '126907',\n", " 'type': 'GENOME_GTF'}],\n", " 'name': 'From INSDC submitter',\n", " 'release_date': 'Oct 07, 2019',\n", " 'release_number': None,\n", " 'report_url': None,\n", " 'source': 'Korea University'},\n", " 'assembly_accession': 'GCF_008831485.1',\n", " 'assembly_category': 'representative genome',\n", " 'assembly_level': 'Complete Genome',\n", " 'bioproject_lineages': [{'bioprojects': [{'accession': 'PRJNA566216',\n", " 'parent_accession': None,\n", " 'title': 'Lactobacillus '\n", " 'acetotolerans '\n", " 'Genome '\n", " 'sequencing '\n", " 'and '\n", " 'assembly'}]}],\n", " 'chromosomes': ['ANONYMOUS'],\n", " 'contig_n50': 1683905,\n", " 'display_name': 'ASM883148v1',\n", " 'estimated_size': '2223244',\n", " 'org': {'assembly_count': None,\n", " 'assembly_counts': {'node': 10, 'subtree': 14},\n", " 'blast_node': None,\n", " 'breed': None,\n", " 'children': None,\n", " 'common_name': None,\n", " 'counts': None,\n", " 'cultivar': None,\n", " 'ecotype': None,\n", " 'icon': None,\n", " 'isolate': None,\n", " 'key': '1600',\n", " 'max_ord': None,\n", " 'merged': None,\n", " 'merged_tax_ids': None,\n", " 'min_ord': None,\n", " 'parent_tax_id': '1578',\n", " 'rank': 'SPECIES',\n", " 'sci_name': 'Lactobacillus acetotolerans',\n", " 'search_text': None,\n", " 'sex': None,\n", " 'strain': 'LA749',\n", " 'tax_id': '1600',\n", " 'title': 'Lactobacillus acetotolerans'},\n", " 'seq_length': '1683905',\n", " 'submission_date': '2019-10-07'},\n", " 'messages': None}\n" ] } ], "source": [ "## download the genome summaries for all RefSeq lactobacillus genome assemblies\n", "genome_summary = genome_api.assembly_descriptors_by_taxon(\n", " taxon=taxid,\n", " limit='all',\n", " filters_refseq_only=True)\n", "print(f\"Genome summaries for {genome_summary.total_count} genome assemblies downloaded and saved to genome_summary.\")\n", "print(f\"As an example, here's the first genome summary for the first genome assembly in the list:\\n\", genome_summary.assemblies[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we'll plot the assemblies by the year they were annotated." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAF0CAYAAAAggv9WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqKUlEQVR4nO3dfXCUhYHH8e8my4uQEggLaIBKI1AVBbRBUaog7s3dWEWOU+ZO9A5BKVJlUKpSXzh7vhC1AaoFbfGljt7Vq57El6G1E1NwNOM1CJkiHgj1FVFCsiEmQhCSvT8cU1MSIC/78sTvZ8YZs0mefDcE88vjk91QPB6PI0mSJH3DZaQ6QJIkSUoHDmNJkiQJh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCYBwqgO+snPnzoQcNxKJUFlZmZBjJ0LQeiF4zUHrBZuTIWi9YHMyBK0XbE6GoPVC8JoT2Zubm9vq6zxjLEmSJOEwliRJkgCHsSRJkgSk0TXGfysej1NfX09jYyOhUKjdx9m1axf79+/vxLLE2rVrF1988QU9e/bs0P2WJElS26TtMK6vr6dbt26Ewx1LDIfDZGZmdlJV4oXD4aYfCo455phU50iSJH1jpO2lFI2NjR0exUEVDodpbGxMdYYkSdI3StoO42/6ZQTf9PsvSZKUbGk7jCVJkqRkCsy1Cg1XT2nf+7Vye+aqF9of00lWrVrF5Zdf7rXEkiRJacAzxin0yCOPsG/fvlRnSJIkCYfxET3zzDNEo1Gi0SjXXXcdH330EZdeeinRaJTp06fz8ccfA7BgwQJeeumlpvcbMWIEAKWlpVxyySVcffXVnHvuuVx77bXE43EeffRRdu3axaWXXsoll1ySkvsmSZKkvwrMpRSpsHXrVn7+85/zwgsvkJOTQ3V1NQsWLODSSy9l+vTpPP3009x+++089thjhz3OW2+9RUlJCcceeywXX3wxZWVlzJ49m1/96lc888wz5OTkJOkeSZIkqTWeMT6M119/nQsvvLBpuPbr148333yTf/zHfwTgn/7pn/jTn/50xOOMHTuW3NxcMjIyGDVqFB999FFCuyVJktR2DuNO8vXHHm5sbOTAgQNNr+vevXvTv2dmZnLw4MGk90mSJOnwvJTiMCZMmMDs2bOZM2dO06UU+fn5PP/881xyySU899xznHnmmQAMGTKETZs2MWXKFP7whz80G8atycrKoq6uzkspJElSm7XlEbt2teG46fDIXakSmGHc3j+kcDjc7jO03/3ud5k/fz6XXHIJGRkZnHLKKdx1111cf/31PPzww+Tk5LBs2TIAZsyYwZVXXkk0GuW8886jV69eRzz+jBkzmDFjBoMGDeLZZ59tV6MkSZI6Rygej8dTHQGwc+fOZi/v3bv3qMblkXRkGKfCV72ddf+TIRKJUFlZmeqMoxa0XrA5GYLWCzYnQ9B6weZkSJfe9j7Hw5GkwxnjRH6Oc3NzW32d1xhLkiRJHMWlFCtXrmTDhg1kZ2dTWFjYdPvvfvc7Xn75ZTIyMjj99NO5/PLLAVi9ejUlJSVkZGRw5ZVXMnbs2ITFS5IkSZ3liMN40qRJ/MM//AMrVqxouu2tt95i/fr13H///XTr1o2amhoAduzYQWlpKUuXLqW6upo777yTn//852RktP3EdJpc4ZEy3/T7L0mSlGxHXKwnn3wyWVlZzW77wx/+wMUXX0y3bt0AyM7OBqCsrIyzzz6bbt26MXDgQI499li2b9/evrCMjEBdG9yZDh482K4fJiRJktR+7XpUik8++YQtW7bw9NNP061bN6644gqGDx9OLBZreipkgJycHGKxWIvHKC4upri4GICCggIikUiz18fjcWKxWIfHcWNjY6DOvjY2NtKtWzcGDRpEKBRKdc5RCYfDh/z5pbOg9YLNyRC0XrA5GYLWCzYnQ7r0tuUh2NoiHe5bqj7H7RrGjY2N1NXVcffdd/OXv/yFZcuW8Ytf/KJNx4hGo0Sj0aaXW/vNw8zMzPYkNkmX3xw9Wl/1VlVVpTrlqAX1cxwkNide0HrB5mQIWi/YnAxB622rdLhvgXpUipycHM444wxCoRDDhw8nIyOD2tpacnJymg26WCzmk1dIkiQpENo1jMeNG8fmzZuBLx9/+ODBg3zrW98iPz+f0tJSDhw4QEVFBZ988gnDhw/v1GBJkiQpEY54KcXy5ct5++23qa2tZe7cuUyfPp3JkyezcuVKFi5cSDgc5kc/+hGhUIihQ4dy1llnccMNN5CRkcHs2bP9JTJJkiQFwhGH8YIFC1q8ff78+S3ePm3aNKZNm9ahKEmSJCnZPJ0rSZIk4TCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEkAhFMdIEmSpK6v4eopR/22u9p47MxVL7TxPVrmGWNJkiQJh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCTiKYbxy5UquuuoqFi5ceMjrXnzxRaZPn85nn30GQDwe57HHHuO6667jxz/+Me+++27nF0uSJEkJcMRhPGnSJG655ZZDbq+srOTPf/4zkUik6baNGzfy6aef8sADDzBnzhweeeSRzq2VJEmSEuSIw/jkk08mKyvrkNufeOIJZsyYQSgUarpt/fr1nHvuuYRCIUaOHMnnn39OdXV15xZLkiRJCdCuJ/goKysjJyeHYcOGNbs9Fos1O4Pcv39/YrEY/fr1O+QYxcXFFBcXA1BQUNDs/TpTOBxO2LETIWi9ELzmoPWCzckQtF6wORmC1gs2J0O69Lb1STCOVqLuW6J6ofOa2zyM9+/fz+rVq7nttts69IGj0SjRaLTp5crKyg4drzWRSCRhx06EoPVC8JqD1gs2J0PQesHmZAhaL9icDEHrbasg3re2NOfm5rb6ujYP4127dlFRUcGNN94IQFVVFTfffDNLliwhJyenWVhVVRU5OTlt/RCSJElS0rV5GH/7299u9kt1P/rRj1iyZAl9+vQhPz+f3//+90yYMIFt27bRq1evFi+jkCRJktLNEYfx8uXLefvtt6mtrWXu3LlMnz6dyZMnt/i2p512Ghs2bGD+/Pl0796defPmdXqwJEmSlAhHHMYLFiw47OtXrFjR9O+hUIirrrqqw1GSJElSsvnMd5IkSRIOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkSAOFUB0iSpK6l4eopbXr7XW1428xVL7QtRmoDzxhLkiRJOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRJwFI9jvHLlSjZs2EB2djaFhYUAPPnkk7z55puEw2EGDRrEvHnz6N27NwCrV6+mpKSEjIwMrrzySsaOHZvQOyBJkiR1hiOeMZ40aRK33HJLs9tGjx5NYWEhP/vZzzjuuONYvXo1ADt27KC0tJSlS5dy66238uijj9LY2JiYckmSJKkTHXEYn3zyyWRlZTW7bcyYMWRmZgIwcuRIYrEYAGVlZZx99tl069aNgQMHcuyxx7J9+/YEZEuSJEmdq8NPCV1SUsLZZ58NQCwWY8SIEU2vy8nJaRrNf6u4uJji4mIACgoKiEQiHU1pUTgcTtixEyFovRC85qD1gs3JELResDkZgtYL6dHclqd4bqtU3zdIj88xJO7znKj7FoSviw4N4+eee47MzEzOOeecNr9vNBolGo02vVxZWdmRlFZFIpGEHTsRgtYLwWsOWi/YnAxB6wWbkyFovRDM5rZIh/vm5zj9tKU5Nze31de1+1Ep1q5dy5tvvsn8+fMJhULAl2eIq6qqmt4mFouRk5PT3g8hSZIkJU27hnF5eTnPP/88N998Mz169Gi6PT8/n9LSUg4cOEBFRQWffPIJw4cP77RYSZIkKVGOeCnF8uXLefvtt6mtrWXu3LlMnz6d1atXc/DgQe68804ARowYwZw5cxg6dChnnXUWN9xwAxkZGcyePZuMDB8qWZIkSenviMN4wYIFh9w2efLkVt9+2rRpTJs2rUNRkiRJUrJ5OleSJEnCYSxJkiQBDmNJkiQJcBhLkiRJgMNYkiRJAhzGkiRJEtDBp4SWJEnqChqunnLUb7urDcfNXPVC22OUMp4xliRJknAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEQPhIb7By5Uo2bNhAdnY2hYWFANTV1bFs2TJ2797NgAEDuP7668nKyiIej/P444+zceNGevTowbx588jLy0v4nZAkSZI66ohnjCdNmsQtt9zS7LaioiJOPfVUHnjgAU499VSKiooA2LhxI59++ikPPPAAc+bM4ZFHHklItCRJktTZjjiMTz75ZLKysprdVlZWxsSJEwGYOHEiZWVlAKxfv55zzz2XUCjEyJEj+fzzz6murk5AtiRJktS5jngpRUtqamro168fAH379qWmpgaAWCxGJBJperv+/fsTi8Wa3vbriouLKS4uBqCgoKDZ+3WmcDicsGMnQtB6IXjNQesFm5MhaL1gczIErRfSo3lXAo+dqPuWqOZE/lkErTkIXxftGsZfFwqFCIVCbX6/aDRKNBptermysrKjKS2KRCIJO3YiBK0XgtcctF6wORmC1gs2J0PQeiGYzW0RtPsWtF7o+s25ubmtvq5dj0qRnZ3ddIlEdXU1ffr0ASAnJ6dZWFVVFTk5Oe35EJIkSVJStWsY5+fns27dOgDWrVvHuHHjmm5/9dVXicfjvPPOO/Tq1avFyygkSZKkdHPESymWL1/O22+/TW1tLXPnzmX69OlMnTqVZcuWUVJS0vRwbQCnnXYaGzZsYP78+XTv3p158+Yl/A5IkiRJneGIw3jBggUt3r548eJDbguFQlx11VUdjpIkSZKSzWe+kyRJknAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIACHfknV966SVKSkoIhUIMHTqUefPmsWfPHpYvX05tbS15eXlcd911hMMd+jCSJElSwrX7jHEsFuN3v/sdBQUFFBYW0tjYSGlpKU899RQ/+MEPePDBB+nduzclJSWd2StJkiQlRIcupWhsbOSLL76goaGBL774gr59+7J582bGjx8PwKRJkygrK+uUUEmSJCmR2n2NQ05ODhdddBHXXHMN3bt3Z8yYMeTl5dGrVy8yMzOb3iYWi7X4/sXFxRQXFwNQUFBAJBJpb8phhcPhhB07EYLWC8FrDlov2JwMQesFm5MhaL2QHs27EnjsRN23RDUn8s8iaM1B+Lpo9zCuq6ujrKyMFStW0KtXL5YuXUp5eflRv380GiUajTa9XFlZ2d6Uw4pEIgk7diIErReC1xy0XrA5GYLWCzYnQ9B6IZjNbRG0+xa0Xuj6zbm5ua2+rt3DeNOmTQwcOJA+ffoAcOaZZ7J161b27t1LQ0MDmZmZxGIxcnJy2vshJEmSpKRp9zXGkUiEbdu2sX//fuLxOJs2bWLIkCGMGjWKN954A4C1a9eSn5/fabGSJElSorT7jPGIESMYP348N998M5mZmQwbNoxoNMrpp5/O8uXLefrpp/nOd77D5MmTO7NXkiRJSogOPcDw9OnTmT59erPbBg0axJIlSzoUJUmSJCWbz3wnSZIk4TCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBEA41QGSJCVTw9VTjvptd7XhuJmrXmh7jKS04hljSZIkCYexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSgA4+893nn3/Oww8/zEcffUQoFOKaa64hNzeXZcuWsXv3bgYMGMD1119PVlZWZ/VKkiRJCdGhYfz4448zduxYFi5cyMGDB9m/fz+rV6/m1FNPZerUqRQVFVFUVMTll1/eWb2SJElSQrT7Uoq9e/fyf//3f0yePBmAcDhM7969KSsrY+LEiQBMnDiRsrKyzimVJEmSEqjdZ4wrKiro06cPK1eu5IMPPiAvL4+ZM2dSU1NDv379AOjbty81NTUtvn9xcTHFxcUAFBQUEIlE2ptyWOFwOGHHToSg9ULwmoPWCzYnQ9B6web22pWg46b6fn2lK3+OIXGf5yB+XQStOQhfF+0exg0NDbz33nvMmjWLESNG8Pjjj1NUVNTsbUKhEKFQqMX3j0ajRKPRppcrKyvbm3JYkUgkYcdOhKD1QvCag9YLNidD0HrB5nSTLverK3+OIX0+z0craL3Q9Ztzc3NbfV27L6Xo378//fv3Z8SIEQCMHz+e9957j+zsbKqrqwGorq6mT58+7f0QkiRJUtK0exj37duX/v37s3PnTgA2bdrEkCFDyM/PZ926dQCsW7eOcePGdU6pJEmSlEAdelSKWbNm8cADD3Dw4EEGDhzIvHnziMfjLFu2jJKSkqaHa5MkSZLSXYeG8bBhwygoKDjk9sWLF3fksJIk6Wsarp5y1G/bll9wylz1QttjpC7MZ76TJEmScBhLkiRJgMNYkiRJAhzGkiRJEuAwliRJkgCHsSRJkgQ4jCVJkiTAYSxJkiQBDmNJkiQJcBhLkiRJgMNYkiRJAhzGkiRJEuAwliRJkgCHsSRJkgQ4jCVJkiTAYSxJkiQBDmNJkiQJcBhLkiRJgMNYkiRJAhzGkiRJEuAwliRJkgCHsSRJkgQ4jCVJkiTAYSxJkiQBDmNJkiQJcBhLkiRJgMNYkiRJAhzGkiRJEuAwliRJkgAId/QAjY2NLFq0iJycHBYtWkRFRQXLly+ntraWvLw8rrvuOsLhDn8YSZIkKaE6fMZ4zZo1DB48uOnlp556ih/84Ac8+OCD9O7dm5KSko5+CEmSJCnhOjSMq6qq2LBhA+effz4A8XiczZs3M378eAAmTZpEWVlZxyslSZKkBOvQNQ6//vWvufzyy9m3bx8AtbW19OrVi8zMTABycnKIxWItvm9xcTHFxcUAFBQUEIlEOpLSqnA4nLBjJ0LQeiF4zUHrBZuTIWi9YHN77UrQcRN5v4LWnKheCF6zXxd/FYSvi3YP4zfffJPs7Gzy8vLYvHlzm98/Go0SjUabXq6srGxvymFFIpGEHTsRgtYLwWsOWi/YnAxB6wWb000Q75fNiRe0Xuj6zbm5ua2+rt3DeOvWraxfv56NGzfyxRdfsG/fPn7961+zd+9eGhoayMzMJBaLkZOT094PIUmSJCVNu4fxZZddxmWXXQbA5s2befHFF5k/fz5Lly7ljTfeYMKECaxdu5b8/PxOi5UkSZISpdMfx3jGjBm89NJLXHfdddTV1TF58uTO/hCSJElSp+uUBxgeNWoUo0aNAmDQoEEsWbKkMw4rSZIkJY3PfCdJkiThMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIAh7EkSZIEOIwlSZIkwGEsSZIkAQ5jSZIkCXAYS5IkSYDDWJIkSQIcxpIkSRLgMJYkSZIACKc6QJL0Vw1XTznqt93VhuNmrnqh7TGS9A3jGWNJkiSJDpwxrqysZMWKFezZs4dQKEQ0GuWCCy6grq6OZcuWsXv3bgYMGMD1119PVlZWZzZLkiRJna7dwzgzM5MrrriCvLw89u3bx6JFixg9ejRr167l1FNPZerUqRQVFVFUVMTll1/emc2SJElSp2v3pRT9+vUjLy8PgGOOOYbBgwcTi8UoKytj4sSJAEycOJGysrLOKZUkSZISqFOuMa6oqOC9995j+PDh1NTU0K9fPwD69u1LTU1NZ3wISZIkKaE6/KgU9fX1FBYWMnPmTHr16tXsdaFQiFAo1OL7FRcXU1xcDEBBQQGRSKSjKS0Kh8MJO3YiBK0XgtecLr27/vHso3/bNhx30OrStsckQLp8no9WuvS25c+6LdLhvkF6fJ6D+DkOWnOieiF4zX5d/FUQvi46NIwPHjxIYWEh55xzDmeeeSYA2dnZVFdX069fP6qrq+nTp0+L7xuNRolGo00vV1ZWdiSlVZFIJGHHToSg9ULwmoPW21bpct+C9nkOWm9bpct968qf5yDeL5sTL2i90PWbc3NzW31duy+liMfjPPzwwwwePJgLL7yw6fb8/HzWrVsHwLp16xg3blx7P4QkSZKUNO0+Y7x161ZeffVVvv3tb3PjjTcC8C//8i9MnTqVZcuWUVJS0vRwbZIkSVK6a/cwPvHEE/ntb3/b4usWL17c7iBJkiQpFXzmO0mSJAmHsSRJkgR0wsO1SfrmaLh6ylG/bVselidz1Qttj1FaaMvXBPh1ISm9ecZYkiRJwmEsSZIkAQ5jSZIkCfAaY3URibr2FbzOUZKkbwrPGEuSJEk4jCVJkiTAYSxJkiQBDmNJkiQJcBhLkiRJgMNYkiRJAhzGkiRJEuAwliRJkgCf4ENSF+YTv0iS2sIzxpIkSRIOY0mSJAlwGEuSJEmAw1iSJEkCHMaSJEkS4DCWJEmSAIexJEmSBDiMJUmSJMBhLEmSJAEOY0mSJAlwGEuSJEkAhFMd8E3QcPWUo37bXW04buaqF9oeI0mSpBZ5xliSJEkigWeMy8vLefzxx2lsbOT8889n6tSpifpQSgDPckuSpG+ahJwxbmxs5NFHH+WWW25h2bJlvP766+zYsSMRH0qSJEnqFAk5Y7x9+3aOPfZYBg0aBMDZZ59NWVkZQ4YM6fCx23ImEzybKUmSpKMTisfj8c4+6BtvvEF5eTlz584F4NVXX2Xbtm3Mnj276W2Ki4spLi4GoKCgoLMTJEmSpDZJ2S/fRaNRCgoKEj6KFy1alNDjd7ag9ULwmoPWCzYnQ9B6weZkCFov2JwMQeuF4DWnqjchwzgnJ4eqqqqml6uqqsjJyUnEh5IkSZI6RUKG8QknnMAnn3xCRUUFBw8epLS0lPz8/ER8KEmSJKlTJOSX7zIzM5k1axZ33303jY2NnHfeeQwdOjQRH+qIotFoSj5uewWtF4LXHLResDkZgtYLNidD0HrB5mQIWi8ErzlVvQn55TtJkiQpaHzmO0mSJAmHsSRJkgQ4jCVJkiTAYSxJkiQBXXwY/9d//VeqEyRJkhQQCXm4tlR47LHHDrnt1Vdfpb6+HoBZs2YlO+mIKisr6dOnD927dycej7N27Vree+89hgwZwvnnn09mZmaqEw+xfv16Ro8eTffu3VOdctTq6+spLy+nsrKSjIwMcnNzGT16NBkZ6ftz4ccff0xZWRmxWAz48klz8vPzGTJkSIrL2uaPf/wj5513XqozWvTxxx8Ti8UYMWIEPXv2bLq9vLycsWPHpi7sMLZv3w7A8OHD2bFjB+Xl5eTm5nL66aenuOzo/OIXv+Daa69NdcZR27JlC9u3b2fo0KGMGTMm1Tkt2rZtG4MHD6ZXr1588cUXFBUV8e677zJkyBCmTZtGr169Up3YzJo1azjjjDOIRCKpTjlqBw8e5PXXX6dfv36MHj2a1157ja1btzJ48GCi0SjhcHpOqV27dvG///u/VFVVkZGRwXHHHcf3v//9tPua+Eq6fN/rMg/Xds0113DSSScxZswYvrpLTz75JFdccQUAkyZNSmFdyxYuXMg999xDjx49eOqpp9i1axfjxo3jrbfeAmDevHkpLjzUjBkz6NmzJ2PHjmXChAmMHTs2rQdmaWkpL774IscffzybN29m5MiRxONxPvzwQ+bPn8+3v/3tVCceoqioiNdff50JEyY0PWNkLBZrum3q1KmpDWyDa665hoceeijVGYdYs2YNL7/8MoMHD+aDDz5g5syZjBs3DoCbb76Ze++9N8WFh3rmmWcoLy+noaGB0aNHs23bNkaNGsWmTZsYM2YM06ZNS3ViM3/7OYzH42zevJlTTjkF+PLznG5+8pOfsGTJEgCKi4t5+eWXOeOMM/jzn//M9773vbT8u3fDDTdw//33k5mZyS9/+Ut69OjB+PHj2bRpEx988AE//vGPU53YzL/927/Rs2dPBg0axIQJEzjrrLPo06dPqrMO64EHHqChoYH9+/fTu3dv6uvrOfPMM9m0aRPxeDwtf9hbs2YNGzZs4KSTTmLjxo0MGzaM3r1786c//YmrrrqKUaNGpTqxmXT6vpeeP+a0w9KlS/nv//5vysvLueKKK8jJyeHZZ59Ny0H8lcbGRnr06AHApk2bWLJkCRkZGZx77rnceOONKa5r2eDBg1m8eDFvvPEGL730Eg899BDjxo3j+9//PieffHKq8w7x3HPPcffdd9OjRw8+++wzHnzwQW699VY++OADfvWrX3HXXXelOvEQf/zjHyksLDzkLMSFF17IDTfckHbfnFv7xhuPx6mpqUlyzdF55ZVXuPfee+nZsycVFRUsXbqU3bt3c8EFF5Cu5wreeOMN7r//fg4cOMCcOXN46KGH6NWrF1OmTOGWW25Ju2Eci8UYPHgw559/PqFQiHg8zrvvvstFF12U6rRWNTQ0NP37K6+8wu23306fPn246KKLuPXWW9Pu7x58+ffsq/+7+O677zb9QHLiiSem5feRQYMGUVBQwKZNmygtLeW3v/0teXl5TJgwgTPPPJNjjjkm1YmH+PDDD/nZz35GQ0MDc+fO5Ze//CUZGRmcc845afk5hi+/fu+//34yMjK48MILWbJkCXfccQd/93d/x3333cd9992X6sRm0un7XpcZxscccwwzZ87k3Xff5cEHH+S0005L229wX4lEIrz11luccsopDBgwgKqqKgYMGEBtbW2q01oVCoXIysoiGo0SjUbZs2cPpaWl/Od//iexWCztzg7G4/Gmyz569uzZNNSOP/549u3bl8q0VoVCIaqrqxkwYECz26urqwmFQimqal1NTQ233norvXv3bnZ7PB7n9ttvT1HV4cXj8abLJwYOHMgdd9xBYWEhu3fvTtv/bmRmZpKRkUGPHj0YNGhQ0/8O7d69e1p+XSxZsoQ1a9bw3HPPccUVVzBs2DC6d++elj9AfyUej1NXV0c8HicejzedyezZs2daXtoGMHTo0KZLlo4//nj+8pe/cMIJJ7Bz5860/F/8oVCIjIwMxowZw5gxYzh48CDl5eW89tprPPnkkzz66KOpTjxEPB7n4MGD1NfXs3//fvbu3UtWVhYHDhxo9sNUumloaCAjI4MDBw40XVYaiUTSsjmdvu+l39+aDsrLy2Px4sW8/PLLfPe73011zmH98Ic/ZMWKFTzzzDMcc8wx3HTTTQwbNozPP/+cf/3Xf011Xov+djT07duXCy64gAsuuIDdu3enqKp1p512Gvfccw8nnXQS5eXljB8/HqDpm186mjlzJv/xH//BcccdR//+/YEvr0f/9NNPmT17dorrDnX66adTX1/PsGHDDnlduo6g7Oxs3n///abmnj17smjRIh566CE+/PDD1Ma1IhwOs3//fnr06EFBQUHT7Xv37k3Ly5m+OlN11lln8cQTT5CdnZ2W35C/bu/evSxatIh4PN70jbpfv37U19en7X8v5s6dy+OPP85zzz3Ht771LW677Tb69+9P//79+eEPf5jqvEP87ecxHA6Tn59Pfn4++/fvT1HV4Z133nksWLCAxsZG/vmf/5mlS5cycOBAtm3bxtlnn53qvBadf/75/OQnP2H48OFs2bKFiy++GIDPPvuMrKysFNcdKp2+73WZa4y/smfPnmYXbvft2ze1QUdhx44dfPLJJzQ0NNC/f39OOOGEtPxGB7B58+a0uzbpSDZs2MCOHTsYNmwYo0ePBr68jKWhoYFu3bqluK5ljY2NbN++vdnX8vDhw9P26yJoqqqqyMzMbPG/D1u2bOHEE09MftQRHDhwoMWv188++4w9e/ak5fXyX7dhwwa2bNnCZZddluqUNtu/fz81NTUMHDgw1Smt2rt3LxUVFTQ2Nqb1976dO3eSm5ub6ow2+/p/iz///HM2bdpEJBJh+PDhKS5r3UcffcTHH3/M0KFDGTx4cKpzjihdvu91mWH8/vvvs2rVKvbu3dt04XZVVRW9e/dm9uzZ5OXlpbiwdUEc80FrDlpva+rr65s9gkK6C1ov2JwMQesFm5MhaL1gc1fUZYbxjTfeyJw5cxgxYkSz29955x1WrVrF/fffn6Ky1gVxzAetOWi9R5Kuj/LQmqD1gs3JELResDkZgtYLNneWr34hPhaLMXbsWGbMmNF0ycfXHy0mGbrMNcb79+8/ZBQDjBw5sumi83SzYsWKVsf8Qw89lJZjPmjNQesFeOmll1q8PR6Pp+XXctB6weZkCFov2JwMQesFm5PhkUce4dJLL2XEiBG88sorLF68mJtuuoljjz026b+b0GUuWBw7dixLliyhtLSUrVu3snXrVkpLS1myZEnaPlh/EMd80JqD1gvwm9/8hrq6Ovbt29fsn3T9BaCg9YLNyRC0XrA5GYLWCzYnQ319PWPHjqV3795MmTKFWbNmcc899/DOO+/4qBTtNWvWLDZu3HjIs6b8/d//fdo+K9RXY37ixIlNv4VZVVXFunXr0nbMB605aL0A3/nOdzjjjDNavMyjpKQkBUWHF7ResDkZgtYLNidD0HrB5mTZu3dv08NQnnLKKSxcuJDCwkLq6uqS2tFlrjEOqpbGfH5+ftqOeQhec9B6d+7cSVZWVovPBrVnz560+8XBoPWCzckQtF6wORmC1gs2J8Nrr73GwIEDGTlyZLPbKysrefbZZ5k7d27SWrrMMN67dy+rV69m/fr17Nmzh1AoRHZ2Nvn5+UydOvWQJx+QJEmSvq7LDOO7776bUaNGMWnSpKafhPbs2cPatWt56623uO2221Ib2IIgjvmgNQetF/7aXFZWRk1NTdo3B60XbE6GoPWCzckQtF6wORnSqjfeRcyfP79dr0ulu+66K7569ep4dXV1023V1dXx1atXx++8887UhR1G0JqD1huPB685aL3xuM3JELTeeNzmZAhabzxuczKkU2+XeVSKAQMG8Pzzz7Nnz56m2/bs2UNRURGRSCR1YYdRUVHB1KlTm13r07dvX6ZOnZqWT68MwWsOWi8ErzlovWBzMgStF2xOhqD1gs3JkE69XWYYL1iwgNraWu644w6uvPJKrrzySn76059SV1fH9ddfn+q8FgVxzAetOWi9ELzmoPWCzckQtF6wORmC1gs2J0M69XaZa4wBPv74Y6qqqhg5cmSzpzssLy9Py4fmqquro6ioiPXr11NTUwN8+RPS9773PaZOndr0rC/pJGjNQeuF4DUHrRdsToag9YLNyRC0XrA5GdKpt8sM4zVr1vDyyy8zePBgPvjgA2bOnMm4ceMAuPnmm7n33ntTXNiyoI15CF5z0HoheM1B6wWbkyFovWBzMgStF2xOhnTp7TKXUrzyyivce++93HTTTfz7v/87//M//8OaNWsA0vJZXuDLMX/ffffx+9//noULF1JWVtb0ut/85jcpLGtd0JqD1gvBaw5aL9icDEHrBZuTIWi9YHMypFNvl3nmu3g83vQTxsCBA7njjjsoLCxk9+7daTuMvxrzPXv2pKKigqVLl7J7924uuOACmztJ0HoheM1B6wWbkyFovWBzMgStF2xOhnTq7TJnjLOzs3n//febXu7ZsyeLFi2itraWDz/8MHVhh9HSmN+4cSNPPPFEWn7hQvCag9YLwWsOWi/YnAxB6wWbkyFovWBzMqRTb5cZxtdee+0hT3GYmZnJtddey09/+tPURB1BEMd80JqD1gvBaw5aL9icDEHrBZuTIWi9YHMypFNvl/nluyCqqqoiMzOzxecs37JlCyeeeGLyo44gaM1B64XgNQetF2xOhqD1gs3JELResDkZ0qnXYSxJkiTRhS6lkCRJkjrCYSxJkiThMJYkSZIAh7EkSZIEwP8DzbWzrVRVvbsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def annotation_year(annotation_metadata):\n", " return datetime.strptime(annotation_metadata.release_date, '%b %d, %Y').year\n", "\n", "annots_by_year = Counter()\n", "no_annot_assms = []\n", "for d in map(lambda d: d.assembly, genome_summary.assemblies):\n", " if not d.annotation_metadata:\n", " no_annot_assms.append(d.assembly_accession)\n", " continue\n", " annots_by_year[annotation_year(d.annotation_metadata)] += 1\n", "\n", "\n", "if len(no_annot_assms) > 0:\n", " print(dedent(f'''\n", " WARNING!\n", " Some assemblies do not have annotation. \n", " Most likely, this is because of an indexing delay. Skipping {len(no_annot_assms)} assemblies.\n", "'''))\n", "\n", "df = pd.DataFrame.from_dict(annots_by_year, orient='index', columns=['count']).sort_index()\n", "df.plot(kind='bar', y='count', figsize=(12,6))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can make a list of the genome assemblies that were annotated in 2020 and count the number of genomes in the list." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "145 genome assemblies were annotated in 2020.\n", "Here are the first five genome assemblies in the list: ['GCF_014648715.1', 'GCF_012848655.1', 'GCF_013867605.1', 'GCF_013867555.1', 'GCF_013342945.1']\n" ] } ], "source": [ "year_of_interest = 2020\n", "genome_accessions = []\n", "for d in map(lambda d: d.assembly, genome_summary.assemblies):\n", " if d.annotation_metadata:\n", " if annotation_year(d.annotation_metadata) == year_of_interest:\n", " genome_accessions.append(d.assembly_accession)\n", " \n", "print(f'{len(genome_accessions)} genome assemblies were annotated in {year_of_interest}.')\n", "print(f\"Here are the first five genome assemblies in the list: {genome_accessions[0:5]}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could download genome datasets for all of the returned assemblies, based on the accessions stored in `genome_accessions`. If there are quite a few of them, a `download_api.download_assembly_package_post()` request would be required (see docs here)[https://github.com/ncbi/datasets/blob/master/client_docs/python/docs/DownloadApi.md#download_assembly_package_post]. \n", "\n", "For demonstration purposes, we'll just download genome datasets for five genome assemblies through a GET request." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Download genome datasets for the first 5 accessions: ['GCF_009857225.1', 'GCF_014830125.1', 'GCF_013394955.1', 'GCF_013370935.1', 'GCF_013249185.1']\n", "Download saved to ncbi_genomes.zip\n", "Archive: ncbi_genomes.zip\n", " Length Date Time Name\n", "--------- ---------- ----- ----\n", " 661 11-19-2020 12:03 README.md\n", " 5811 11-19-2020 12:03 ncbi_dataset/data/assembly_data_report.jsonl\n", " 1668120 11-19-2020 12:03 ncbi_dataset/data/GCF_009857225.1/GCF_009857225.1_ASM985722v1_genomic.fna\n", " 898917 11-19-2020 12:03 ncbi_dataset/data/GCF_009857225.1/genomic.gff\n", " 597383 11-19-2020 12:03 ncbi_dataset/data/GCF_009857225.1/protein.faa\n", " 1838908 11-19-2020 12:03 ncbi_dataset/data/GCF_013249185.1/GCF_013249185.1_ASM1324918v1_genomic.fna\n", " 1062407 11-19-2020 12:03 ncbi_dataset/data/GCF_013249185.1/genomic.gff\n", " 594208 11-19-2020 12:03 ncbi_dataset/data/GCF_013249185.1/protein.faa\n", " 1786869 11-19-2020 12:03 ncbi_dataset/data/GCF_013370935.1/GCF_013370935.1_ASM1337093v1_genomic.fna\n", " 1011101 11-19-2020 12:03 ncbi_dataset/data/GCF_013370935.1/genomic.gff\n", " 590951 11-19-2020 12:03 ncbi_dataset/data/GCF_013370935.1/protein.faa\n", " 1779401 11-19-2020 12:03 ncbi_dataset/data/GCF_013394955.1/GCF_013394955.1_ASM1339495v1_genomic.fna\n", " 1010905 11-19-2020 12:03 ncbi_dataset/data/GCF_013394955.1/genomic.gff\n", " 589700 11-19-2020 12:03 ncbi_dataset/data/GCF_013394955.1/protein.faa\n", " 2064899 11-19-2020 12:03 ncbi_dataset/data/GCF_014830125.1/GCF_014830125.1_ASM1483012v1_genomic.fna\n", " 1167106 11-19-2020 12:03 ncbi_dataset/data/GCF_014830125.1/genomic.gff\n", " 694896 11-19-2020 12:03 ncbi_dataset/data/GCF_014830125.1/protein.faa\n", " 11623 11-19-2020 12:03 ncbi_dataset/data/GCF_009857225.1/sequence_report.jsonl\n", " 19774 11-19-2020 12:03 ncbi_dataset/data/GCF_013249185.1/sequence_report.jsonl\n", " 7913 11-19-2020 12:03 ncbi_dataset/data/GCF_013370935.1/sequence_report.jsonl\n", " 7909 11-19-2020 12:03 ncbi_dataset/data/GCF_013394955.1/sequence_report.jsonl\n", " 34330 11-19-2020 12:03 ncbi_dataset/data/GCF_014830125.1/sequence_report.jsonl\n", " 2536 11-19-2020 12:03 ncbi_dataset/data/dataset_catalog.json\n", "--------- -------\n", " 17446328 23 files\n", "CPU times: user 54.5 ms, sys: 44 ms, total: 98.5 ms\n", "Wall time: 1.07 s\n" ] } ], "source": [ "%%time\n", "## Download data\n", "# Let's use accessions #21-25 in the list because these all have at least one crispr gene\n", "print(f'Download genome datasets for 5 accessions: {genome_accessions[20:25]}')\n", "dl_response = genome_api.download_assembly_package(\n", " genome_accessions[20:25],\n", " exclude_sequence=False,\n", " include_annotation_type=['GENOME_GFF', 'PROT_FASTA'],\n", " _preload_content=False )\n", "\n", "## write to a zip file \n", "zipfn = 'ncbi_genomes.zip'\n", "with open(zipfn, 'wb') as f:\n", " f.write(dl_response.data)\n", "print(f'Download saved to {zipfn}')\n", "!unzip -l {zipfn}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we're going to look for crispr genes in these genome datasets, and extract the corresponding crispr gene and protein sequences. \n", "We're going to break this down into a few steps. \n", "First, let's make a list of the data files that are available for each genome assembly. We can use the `dataset_catalog.json` file to tell us which data files are available for each genome assembly." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Catalog found with metadata for 5 assemblies\n", "Assemblies:\n", "GCF_009857225.1, GCF_013249185.1, GCF_013370935.1, GCF_013394955.1, GCF_014830125.1\n", "Files for type: GFF3\n", "defaultdict(,\n", " { 'GCF_009857225.1': [ 'ncbi_dataset/data/GCF_009857225.1/genomic.gff'],\n", " 'GCF_013249185.1': [ 'ncbi_dataset/data/GCF_013249185.1/genomic.gff'],\n", " 'GCF_013370935.1': [ 'ncbi_dataset/data/GCF_013370935.1/genomic.gff'],\n", " 'GCF_013394955.1': [ 'ncbi_dataset/data/GCF_013394955.1/genomic.gff'],\n", " 'GCF_014830125.1': [ 'ncbi_dataset/data/GCF_014830125.1/genomic.gff']})\n", "Files for type: FASTA\n", "defaultdict(, {})\n", "Files for type: GENOMIC_NUCLEOTIDE_FASTA\n", "defaultdict(,\n", " { 'GCF_009857225.1': [ 'ncbi_dataset/data/GCF_009857225.1/GCF_009857225.1_ASM985722v1_genomic.fna'],\n", " 'GCF_013249185.1': [ 'ncbi_dataset/data/GCF_013249185.1/GCF_013249185.1_ASM1324918v1_genomic.fna'],\n", " 'GCF_013370935.1': [ 'ncbi_dataset/data/GCF_013370935.1/GCF_013370935.1_ASM1337093v1_genomic.fna'],\n", " 'GCF_013394955.1': [ 'ncbi_dataset/data/GCF_013394955.1/GCF_013394955.1_ASM1339495v1_genomic.fna'],\n", " 'GCF_014830125.1': [ 'ncbi_dataset/data/GCF_014830125.1/GCF_014830125.1_ASM1483012v1_genomic.fna']})\n", "Files for type: PROTEIN_FASTA\n", "defaultdict(,\n", " { 'GCF_009857225.1': [ 'ncbi_dataset/data/GCF_009857225.1/protein.faa'],\n", " 'GCF_013249185.1': [ 'ncbi_dataset/data/GCF_013249185.1/protein.faa'],\n", " 'GCF_013370935.1': [ 'ncbi_dataset/data/GCF_013370935.1/protein.faa'],\n", " 'GCF_013394955.1': [ 'ncbi_dataset/data/GCF_013394955.1/protein.faa'],\n", " 'GCF_014830125.1': [ 'ncbi_dataset/data/GCF_014830125.1/protein.faa']})\n" ] } ], "source": [ "## function to make list of files \n", "from ncbi.datasets.v1alpha1 import catalog_pb2\n", "import pprint\n", "pp = pprint.PrettyPrinter(indent=4)\n", "\n", "def retrieve_data_catalog(zip_file):\n", " with zipfile.ZipFile(zip_file, 'r') as zip:\n", " data_catalog = json.loads(zip.read('ncbi_dataset/data/dataset_catalog.json'))\n", " print(f\"Catalog found with metadata for {len(data_catalog['assemblies'])-1} assemblies\")\n", " return data_catalog\n", "\n", "def get_assemblies(data_catalog):\n", " return [x['accession'] for x in data_catalog['assemblies'] if 'accession' in x]\n", "\n", "# Temporary hack to support GENOMIC_NUCLEOTIDE_FASTA & PROTEIN_FASTA \n", "# which will be present in the next release\n", "def get_file_list(data_catalog, desired_filetype):\n", " desired_filetype = desired_filetype.upper()\n", " if desired_filetype not in catalog_pb2.File.FileType.keys():\n", " raise Exception(f'Filetype {desired_filetype} is invalid.')\n", " \n", " files = defaultdict(list)\n", " for asm in data_catalog['assemblies']:\n", " if 'accession' in asm:\n", " acc = asm['accession'] \n", " for f in asm['files']:\n", " filepath = os.path.join('ncbi_dataset', 'data', f['filePath'])\n", " if f['fileType'] == 'FASTA' and desired_filetype in ('GENOMIC_NUCLEOTIDE_FASTA') and filepath.endswith('fna'):\n", " files[acc].append(filepath)\n", " continue\n", " if f['fileType'] == 'GFF3':\n", " if desired_filetype in ('PROTEIN_FASTA') and filepath.endswith('faa'):\n", " files[acc].append(filepath)\n", " if desired_filetype in ('GFF3') and filepath.endswith('gff'):\n", " files[acc].append(filepath)\n", " continue\n", " if f['fileType'] == desired_filetype:\n", " files[acc].append(filepath)\n", " \n", " return files\n", "\n", "data_catalog = retrieve_data_catalog(zipfn)\n", "print(f'Assemblies:')\n", "print(', '.join(get_assemblies(data_catalog)))\n", "for file_type in ['GFF3', 'FASTA', 'GENOMIC_NUCLEOTIDE_FASTA', 'PROTEIN_FASTA']:\n", " file_list = get_file_list(data_catalog, file_type)\n", " print(f'Files for type: {file_type}')\n", " pp.pprint(file_list)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're going to create temporary files to store the genome annotation (gff3 format) and protein sequence data. We will put the files in a temporary directory named tempfiles." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "## setting up files and directories\n", "\n", "## temporary files; will be deleted at the end\n", "temp_dir = 'tempfiles'\n", "temp_gff = temp_dir + '/temp.gff'\n", "temp_fa = temp_dir + '/temp.fa'\n", "\n", "!mkdir -p {temp_dir}\n", "\n", "## final output files \n", "genes_fn = 'crispr_genes.fna'\n", "prots_fn = 'crispr_proteins.faa'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need a few functions to process the data. One of these functions will create a gff3 database in memory using the gffutils python package. Another combines fasta files, a third function will extract gene information by gene symbol from the gff3 database. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "def create_gff3_db(files_by_assembly, temp_file, zfh):\n", " '''\n", " create gff3 db in memory, per assembly\n", " okay for bacterial assemblies but use caution\n", " when parsing large assemblies like human\n", " '''\n", "\n", " db = {}\n", " for assembly_accession, files in files_by_assembly.items():\n", " with open(temp_file, 'wb') as f:\n", " for file in files:\n", " print(f'\\tWrite {file} to {temp_file}')\n", " f.write(zfh.read(file))\n", " db[assembly_accession] = gffutils.create_db(\n", " temp_file,\n", " dbfn = ':memory:',\n", " force=True,\n", " keep_order=True,\n", " merge_strategy='merge',\n", " sort_attribute_values=True)\n", "\n", "\n", " return db\n", "\n", "\n", "def combine_fasta(files_by_assembly, temp_file, zfh):\n", " '''\n", " Combine fasta for *all* FASTA files, nt & protein\n", " '''\n", " with open(temp_file, 'wb') as f:\n", " for assembly, files in files_by_assembly.items():\n", " for file in files:\n", " print(f'\\tAppending {file} to {temp_file}')\n", " f.write(zfh.read(file))\n", "\n", " print(f'Create FASTA object for {temp_file}.')\n", " return Fasta(temp_file, read_long_names=False, duplicate_action='first')\n", " \n", "\n", "def extract_genes(gff3_db, assemblies, desired_genes):\n", " crispr_order = defaultdict(list)\n", " crispr_genes = {}\n", "\n", " for assembly_accession in assemblies:\n", " if assembly_accession in assemblies:\n", " for gene in gff3_db[assembly_accession].features_of_type('gene'):\n", " gene_name = gene.attributes.get('Name', None)[0]\n", " if gene_name[:4] not in desired_genes:\n", " continue\n", " gene_range = (gene.start, gene.end)\n", " prot_acc = None\n", " if gene.attributes['gene_biotype'][0] == 'protein_coding':\n", " cds = list(gff3_db[assembly_accession].children(gene, featuretype='CDS'))\n", " prot_acc = cds[0].attributes.get('protein_id', None)[0]\n", "\n", " crispr_genes[gene_name] = ([gene.chrom, gene.strand, gene_range, prot_acc])\n", " crispr_order[assembly_accession].append(gene_name)\n", " return crispr_genes, crispr_order\n", "\n", "def write_fasta(fh, defline, content):\n", " fh.write('>' + defline + '\\n')\n", " fh.write(content + '\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use the above functions to process the datasets zip file. First, create the gff3 database, then extract crispr gene information from the gff3 data, then for each gene, store the crispr gene data and write the fasta sequences." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Create GFF3 Database\n", "\tWrite ncbi_dataset/data/GCF_009857225.1/genomic.gff to tempfiles/temp.gff\n", "\tWrite ncbi_dataset/data/GCF_013249185.1/genomic.gff to tempfiles/temp.gff\n", "\tWrite ncbi_dataset/data/GCF_013370935.1/genomic.gff to tempfiles/temp.gff\n", "\tWrite ncbi_dataset/data/GCF_013394955.1/genomic.gff to tempfiles/temp.gff\n", "\tWrite ncbi_dataset/data/GCF_014830125.1/genomic.gff to tempfiles/temp.gff\n", "Extract crispr genes (12 genes)\n", "Write genes to crispr_genes.fna\n", "Create genomes object\n", "\tAppending ncbi_dataset/data/GCF_009857225.1/GCF_009857225.1_ASM985722v1_genomic.fna to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_013249185.1/GCF_013249185.1_ASM1324918v1_genomic.fna to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_013370935.1/GCF_013370935.1_ASM1337093v1_genomic.fna to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_013394955.1/GCF_013394955.1_ASM1339495v1_genomic.fna to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_014830125.1/GCF_014830125.1_ASM1483012v1_genomic.fna to tempfiles/temp.fa\n", "Create FASTA object for tempfiles/temp.fa.\n", "Create proteome object\n", "\tAppending ncbi_dataset/data/GCF_009857225.1/protein.faa to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_013249185.1/protein.faa to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_013370935.1/protein.faa to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_013394955.1/protein.faa to tempfiles/temp.fa\n", "\tAppending ncbi_dataset/data/GCF_014830125.1/protein.faa to tempfiles/temp.fa\n", "Create FASTA object for tempfiles/temp.fa.\n", "CPU times: user 4.83 s, sys: 59.4 ms, total: 4.89 s\n", "Wall time: 4.96 s\n" ] } ], "source": [ "%%time\n", "crispr_genes = set(['cas3', 'cse1', 'cse2', 'cas7', 'cas5', 'cas6', 'cas4', 'cas1', 'cas2', 'cas5', 'cas8', 'cas9', 'csn2'])\n", "\n", "\n", "## create empty files to add data to\n", "open(genes_fn, 'w').close()\n", "open(prots_fn, 'w').close()\n", "\n", "with zipfile.ZipFile(zipfn, 'r') as zfh:\n", "\n", " print('Create GFF3 Database')\n", " gff3_db = create_gff3_db(get_file_list(data_catalog, 'GFF3'), temp_gff, zfh)\n", " \n", " print(f'Extract crispr genes ({len(crispr_genes)} genes)')\n", " crispr_genes, crispr_order = extract_genes(gff3_db, get_assemblies(data_catalog), crispr_genes)\n", " print(f'Write genes to {genes_fn}')\n", "\n", " print('Create genomes object')\n", " genomes = combine_fasta(get_file_list(data_catalog, 'GENOMIC_NUCLEOTIDE_FASTA'), temp_fa, zfh)\n", " \n", " # Write all genes to genes_fn\n", " with open(genes_fn, 'w') as fh:\n", " for gene_name, gene_info in crispr_genes.items():\n", " chrom, strand, gene_range, prot_acc = gene_info\n", " reverse_complement = True if strand == '-' else False\n", " gene_fasta = genomes.get_seq(chrom, gene_range[0], gene_range[1], rc=reverse_complement)\n", " write_fasta(fh, f'{gene_fasta.fancy_name}|{gene_name}', gene_fasta.seq)\n", "\n", " print('Create proteome object')\n", " proteome = combine_fasta(get_file_list(data_catalog, 'PROTEIN_FASTA'), temp_fa, zfh)\n", " \n", " # Write protein FASTA for each CRISPR gene\n", " with open(prots_fn, 'w') as fh:\n", " for gene_name, gene_info in crispr_genes.items():\n", " chrom, strand, gene_range, prot_acc = gene_info\n", " if prot_acc is not None:\n", " prot_fasta = proteome[prot_acc][:]\n", " write_fasta(fh, f'{prot_fasta.name}|{gene_name}', prot_fasta.seq)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "crispr_genes.fna and crispr_proteins.faa contain the genomic and protein FASTA sequences, respectively, for the crispr genes we identified in all of the five genome assemblies." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GCF_009857225.1 ['csn2']\n", "GCF_013249185.1 ['cas2e', 'cas1e', 'cas6e', 'cas5e', 'cas7e', 'cas3']\n", "GCF_013370935.1 ['cas9', 'cas1', 'cas2', 'csn2']\n", "GCF_013394955.1 ['csn2', 'cas2', 'cas1', 'cas9']\n", "GCF_014830125.1 ['cas3', 'cas5c', 'cas8c', 'cas7c', 'cas4', 'cas1c', 'cas2']\n" ] } ], "source": [ "## order of crispr genes in various assemblies\n", "for k,v in crispr_order.items():\n", " print(k, v)" ] } ], "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.8.1" } }, "nbformat": 4, "nbformat_minor": 4 }