Using the workflow to automate remote database queries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this recipe, we're going to develop a small [`Workflow`](http://scikit-bio.org/docs/0.5.0/generated/skbio.workflow.Workflow.html#skbio.workflow.Workflow) to process sequence data sourced from [RNAcentral](http://rnacentral.org/). The focus of this recipe is to highlight how the `Workflow` object can decompose complex logic into smaller blocks that are easier to verify. The structure of the Notebook is to first develop functions that allow us to query RNAcentral and NCBI, where the former is our primary source of data and the latter is for obtaining lineage details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we're going to import the `ElementTree` to facilitate XML parsing in one of the functions. Second, we're going to import the `requests` library which makes querying remote RESTful interfaces very easy. More detail on `requests` can be found [here](http://docs.python-requests.org/en/latest/). Third, our example will query for RNA sequence, so we'll grab the appropriate scikit-bio object. And last, we're going to import the `Workflow` object, two decorators and a helper object that defines a notion of `not None`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import xml.etree.ElementTree as ET\n",
    "\n",
    "import requests\n",
    "\n",
    "from skbio import RNA\n",
    "from skbio.workflow import Workflow, method, requires, not_none"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining functions for querying remote databases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's put together a function that allows us to query RNAcentral. Details about its RESTful interface can be found [here](http://rnacentral.org/api/), and is an excellent example of a well thought out web interface. This function will allow us to perform arbitrary queries, and it'll return a fair bit of detail about each sequence found. In the future, we may add a standard mechanism for querying RNAcentral, but for right now, we'll just roll our own."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def query_rnacentral(query, max_records=100):\n",
    "    \"\"\"Query RNACentral and yield parsed records\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    query : str\n",
    "        The query to use against RNACentral\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    generator\n",
    "        Yields `dict` of the individual records.\n",
    "    \n",
    "    See Also\n",
    "    --------\n",
    "    http://rnacentral.org/api\n",
    "    \"\"\"\n",
    "    def _request_to_json(url, params=None, headers=None):\n",
    "        # Perform the request, and if valid, fetch the JSON data. We're defining this as a\n",
    "        # function as we'll be performing this type of operation a few times.\n",
    "        r = requests.get(url, params=params, headers=headers)\n",
    "        if r.status_code != 200:\n",
    "            raise IOError(\"We received the following status code: %d\" % r.status_code)\n",
    "\n",
    "        return r.json()\n",
    "\n",
    "    url = 'http://rnacentral.org/api/v1/rna'\n",
    "    \n",
    "    # Tell RNAcentral that we'd like JSON data back\n",
    "    headers = {'Accept': 'application/json'}\n",
    "    \n",
    "    # These are the specific parameters for the query. For example, if our query was foo,\n",
    "    # the resulting URL would be (feel free to copy/paste into your browser too):\n",
    "    # http://rnacentral.org/api/v1/rna/?query=foo&page_size=25&flat=true\n",
    "    \n",
    "    # More information about the specific parameters and their values can be found on the\n",
    "    # API description page: http://rnacentral.org/api\n",
    "    params = {'query': query,   # The specific query\n",
    "              'page_size': 25,  # The number of items to return per page. \n",
    "              'flat': 'True'}   # This expands out the detail about each record\n",
    "    \n",
    "    data = _request_to_json(url, params, headers)\n",
    "    \n",
    "    # The key `next` contains the URL of the next \"page\" of data\n",
    "    next_page = data['next']\n",
    "    \n",
    "    # The key `count` contains the total number of items that match our query\n",
    "    count_in_payload = data['count']\n",
    "    \n",
    "    # And, let's setup a counter to track now many results have been yielded\n",
    "    count_yielded = 0\n",
    "    \n",
    "    while count_in_payload > 0 and count_yielded < max_records:\n",
    "        for record in data['results']:\n",
    "            # Let's only care about records that appear to have sequence and length data\n",
    "            sequence = record.get('sequence')\n",
    "            length = record.get('length')\n",
    "            \n",
    "            if sequence is None or length is None:\n",
    "                continue\n",
    "                \n",
    "            # The \"flat\" form from RNAcentral has a list of external references associated with each sequence.\n",
    "            # This list might contain more than a single entry if the sequence has been deposited multiple times,\n",
    "            # and for the purposes of our program, we're going to consider each one of these as independent.\n",
    "            sequence = {'sequence': sequence, 'length': length}    \n",
    "            for xref in record['xrefs']['results']:\n",
    "                # Now, let's build our notion of a record and yield it so a consumer can operate on it\n",
    "                to_yield = sequence.copy()\n",
    "                to_yield.update(xref)\n",
    "                yield to_yield\n",
    "                count_yielded += 1\n",
    "        \n",
    "        # Finally, if there is another page of data, let's grab it, otherwise we'll terminate the loop\n",
    "        if next_page:\n",
    "            data = _request_to_json(next_page)\n",
    "        else:\n",
    "            break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what do these records that we're getting from RNAcentral look like? Let's explore! For our test, let's search for \"tRNA\", and to be polite, we'll only request a single record as that's all we need right now. We're also going to use a helper function, `pprint`, to format the resulting `dict` better for human consumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'accession': {'citations': 'http://rnacentral.org/api/v1/accession/DQ668905.1:1..363:rRNA/citations',\n",
      "               'description': 'uncultured bacterium partial 16S ribosomal '\n",
      "                              'RNA',\n",
      "               'expert_db_url': '',\n",
      "               'external_id': '',\n",
      "               'gene': '',\n",
      "               'id': 'DQ668905.1:1..363:rRNA',\n",
      "               'optional_id': '',\n",
      "               'organelle': '',\n",
      "               'product': '16S ribosomal RNA',\n",
      "               'rna_type': 'rRNA',\n",
      "               'source_url': 'http://www.ebi.ac.uk/ena/data/view/Non-coding:DQ668905.1:1..363:rRNA',\n",
      "               'species': 'environmental samples uncultured bacterium',\n",
      "               'url': 'http://rnacentral.org/api/v1/accession/DQ668905.1:1..363:rRNA/info'},\n",
      " 'database': 'ENA',\n",
      " 'first_seen': '2014-05-29',\n",
      " 'is_active': True,\n",
      " 'is_expert_db': False,\n",
      " 'last_seen': '2016-03-15',\n",
      " 'length': 363,\n",
      " 'sequence': 'GUAUGCAACUUACCUUUUACUAGAGAAUAGCCAAGAGAAAUUUUGAUUAAUGCUCUAUGUUCUUAUUUACUCGCAUGAGUAAAUAAGCAAAGCUCCGGCGGUAAAAGAUGGGCAUGCGUCCUAUUAGCUUGUAGGUGAGGUAAUGGCUCACCUAAGCUCCGAUAGGUAGGGGUCCUGAGAGGGAGAUCCCCCACACUGGUACUGAGACACGGACCAGACUUCUACGGAAGGCAGCAGUAAGGAAUAUUGGACAAUGGAGGCAACUCUGAUCCAGCCAUGCCGCGUGAAGGAAGACGGCCUUAUGGGUUGUAAACUUCUUUUAUACAGGAAGAAACCUUUCCACGUGUGGAAAGCUGACGGUAC',\n",
      " 'taxid': 77133}\n"
     ]
    }
   ],
   "source": [
    "from pprint import pprint\n",
    "test_query = query_rnacentral('tRNA', max_records=1)\n",
    "pprint(next(test_query))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, that's interesting. Even though we queried for \"tRNA\", we ended up getting back \"rRNA\"! Unfortunately, while the search interface of RNAcentral is quite powerful (and we're intentionally not fully leveraging it here), it is not perfect. So, we'll need to do some data scrubbing on our end if we want to only get specific records. More on this in a little. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, it'd be great to get the lineage associated with the records. Luckily, RNAcentral provides a `taxid`, which corresponds to the NCBI taxonomy. Unfortunately, they don't provide the full lineage. But, we can query NCBI to gather these details. NCBI provides [EUtils](http://www.ncbi.nlm.nih.gov/books/NBK25501/) that let users programmatically query their resources similarly to RNAcentral. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def query_ncbi_lineage(taxon_id):\n",
    "    \"\"\"Obtain the NCBI lineage for a taxon ID\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    taxon_id : int\n",
    "        The taxon ID of interest\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    list or None\n",
    "        Each taxon name or None if unable to retreive the taxon details\n",
    "    \"\"\"\n",
    "    url = \"http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi\"\n",
    "    \n",
    "    # Define our parameters to use in our query\n",
    "    params = {'db': 'taxonomy',  # We want to query the taxonomy database\n",
    "              'id': taxon_id}    # We're requesting detail on the taxon ID specifically\n",
    "    \n",
    "    # Make the request\n",
    "    r = requests.get(url, params=params)\n",
    "    \n",
    "    # Bail if we received a bad status\n",
    "    if r.status_code != 200:\n",
    "        return None\n",
    "    \n",
    "    # NCBI returns XML, so we need to parse the \"content\" of our request into a usable structure\n",
    "    tree = ET.fromstring(r.content)\n",
    "    \n",
    "    # We are only interested in the Lineage key within the tree. There is a fair bit of other data\n",
    "    # present within the structure, however.\n",
    "    lineage = next(tree.iter('Lineage'))\n",
    "    \n",
    "    # The lineage, if we did get it, is delimited by \"; \", so let's go ahead and parse that into \n",
    "    # a friendlier structure\n",
    "    if lineage is not None:\n",
    "        return [v.strip() for v in lineage.text.split(';')]\n",
    "    else:\n",
    "        return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Now we can query NCBI, let's do a quick test to see what lineage we get back that would be associated with the previous RNAcentral record."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['cellular organisms', 'Bacteria', 'Proteobacteria', 'Gammaproteobacteria', 'Enterobacteriales', 'Enterobacteriaceae', 'Salmonella', 'Salmonella enterica', 'Salmonella enterica subsp. enterica']\n"
     ]
    }
   ],
   "source": [
    "print(query_ncbi_lineage(57045))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building a database record processor without the Workflow object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the querying functions are out of the way, we're going to construct a method that parses and filters the records based on some runtime critera. Specifically, this function will:\n",
    "\n",
    "- Obtain the record accession\n",
    "- Filter any records that appear to be environmental\n",
    "- Force the sequence to use the RNA alphabet\n",
    "- Filter any sequences that are to short or two long\n",
    "- Filter based on GC content\n",
    "- Obtain the taxonomy\n",
    "- Determine the \"group\" the sequence is in based on the taxonomy\n",
    "\n",
    "Most of these things are going to be optional, as is typical in a sequence processor where under some conditions you may want to filter, and may want to specify the specific criteria to filter by.\n",
    "\n",
    "In addition, we also want to track our \"failed\" records, which is also common as they can be used to debug issues in the pipeline or get dumped into some other processing stream.\n",
    "\n",
    "And last, we want to track processing statistics so we know how many items were filtered based on what critera."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "\n",
    "failures = []\n",
    "def failed(item):\n",
    "    \"\"\"Store a failed record\n",
    "    \n",
    "    This function is decomposed in case we want to add other logic for dealing\n",
    "    with failures\n",
    "    \"\"\"\n",
    "    global failures\n",
    "    failures.append(item)\n",
    "    \n",
    "def record_parser(items, options):\n",
    "    \"\"\"Parse items based on our runtime options\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    items : Iterable\n",
    "        Iterable of `dict` that contains the record information\n",
    "    options : dict\n",
    "        The options that describe filtering criteria, etc\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    generator\n",
    "        Yields `dict`s of the parsed records\n",
    "    \"\"\"\n",
    "    stats = Counter()\n",
    "    for item in items:\n",
    "        # We're going to be updating item inplace, so let's make a copy to be safe\n",
    "        item = item.copy()\n",
    "        \n",
    "        # The id is delimited by a ':', where the first component is the INSDC accession\n",
    "        id_ = item['accession']['id'].split(':', 1)[0]\n",
    "        \n",
    "        # If the options indicate to, let's filter by \"environmental\" criteria. For instance,\n",
    "        # if the description contains the word \"environmental\" then we want to ignore the record.\n",
    "        # This is common when parsing out named isolates.\n",
    "        if options['filter-description']:\n",
    "            acc = item['accession']\n",
    "            tokens = {t.strip().lower() for t in acc['description'].split()}\n",
    "            \n",
    "            if options['filter-description'].intersection(tokens):\n",
    "                failed(item)\n",
    "                stats['filter-description'] += 1\n",
    "                continue\n",
    "        \n",
    "        # Force the record to use an RNA alphabet\n",
    "        item['sequence'] = item['sequence'].replace('T', 'U')\n",
    "        \n",
    "        # Filter out records that don't meet a minimum length requirement\n",
    "        if options['filter-length'] and item['length'] < options['minimum-length']:\n",
    "            failed(item)\n",
    "            stats['minimum-length'] += 1\n",
    "            continue\n",
    "            \n",
    "        # Filter out records that don't meet a maximum length requirement\n",
    "        if options['filter-length'] and item['length'] > options['maximum-length']:\n",
    "            failed(item)\n",
    "            stats['maximum-length'] += 1\n",
    "            continue\n",
    "        \n",
    "        # Compute GC if necessary\n",
    "        if options['compute-gc']:\n",
    "            gc = float(item['sequence'].count('G') + item['sequence'].count('C'))\n",
    "            gc /= item['length']\n",
    "        else:\n",
    "            gc = None\n",
    "           \n",
    "        # Filter by a minimum GC content\n",
    "        if options['compute-gc'] and options['minimum-gc'] and gc is not None:\n",
    "            if gc < options['minimum-gc']:\n",
    "                failed(item)\n",
    "                stats['minimum-gc'] += 1\n",
    "                continue\n",
    "        \n",
    "        # Filter by a maximum GC content\n",
    "        if options['compute-gc'] and options['maximum-gc'] and gc is not None:\n",
    "            if gc < options['maximum-gc']:\n",
    "                failed(item)\n",
    "                stats['maximum-gc'] += 1\n",
    "                continue\n",
    "        \n",
    "        # If we have a taxonomy ID, then let's grab the lineage\n",
    "        if item.get('taxid', None) is not None:\n",
    "            lineage = query_ncbi_lineage(item['taxid'])\n",
    "        else:\n",
    "            lineage = None\n",
    "            \n",
    "        # Determine the group based on the lineage\n",
    "        group = lineage[-1] if lineage is not None else 'unassigned'\n",
    "        \n",
    "        # Finally, let's save a little bit more state and yield for a subsequent consumer\n",
    "        item['gc'] = gc\n",
    "        item['group'] = group\n",
    "        item['lineage'] = lineage\n",
    "        item['id'] = id_\n",
    "        item['sequence'] = RNA(item['sequence'], metadata={'id': item['id']}, validate=True)\n",
    "\n",
    "        yield item"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that's a pretty complex function. There are a lot of conditionals and short circuiting the loop with continues. Though this particular example probably is not extremely difficult to verify its correctness, it is not difficult to imagine a more complex process (such as QC for amplicon data) that would be much more difficult to verify. Additionally, we can't retreive the `stats` defined without making the return more complex, or the use of global variables. Let's quickly play with the function though so we can see how it behaves. First, lets not define any options and examine the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We raised an exception!\n",
      "<class 'KeyError'>: ('filter-description',)\n"
     ]
    }
   ],
   "source": [
    "options = {}\n",
    "items = query_rnacentral('tRNA')\n",
    "records = record_parser(items, options)\n",
    "\n",
    "try:\n",
    "    pprint(next(records))\n",
    "except Exception as e:\n",
    "    print(\"We raised an exception!\")\n",
    "    print(\"%s: %s\" % (type(e), e.args))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so that's a bit annoying. We could update the function to issue `dict.get` calls, and increase the complexity of our conditional logic, but for how, lets just define usable options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'accession': {'citations': 'http://rnacentral.org/api/v1/accession/DQ668905.1:1..363:rRNA/citations',\n",
      "               'description': 'uncultured bacterium partial 16S ribosomal '\n",
      "                              'RNA',\n",
      "               'expert_db_url': '',\n",
      "               'external_id': '',\n",
      "               'gene': '',\n",
      "               'id': 'DQ668905.1:1..363:rRNA',\n",
      "               'optional_id': '',\n",
      "               'organelle': '',\n",
      "               'product': '16S ribosomal RNA',\n",
      "               'rna_type': 'rRNA',\n",
      "               'source_url': 'http://www.ebi.ac.uk/ena/data/view/Non-coding:DQ668905.1:1..363:rRNA',\n",
      "               'species': 'environmental samples uncultured bacterium',\n",
      "               'url': 'http://rnacentral.org/api/v1/accession/DQ668905.1:1..363:rRNA/info'},\n",
      " 'database': 'ENA',\n",
      " 'first_seen': '2014-05-29',\n",
      " 'gc': 0.465564738292011,\n",
      " 'group': 'environmental samples',\n",
      " 'id': 'DQ668905.1',\n",
      " 'is_active': True,\n",
      " 'is_expert_db': False,\n",
      " 'last_seen': '2016-03-15',\n",
      " 'length': 363,\n",
      " 'lineage': ['cellular organisms', 'Bacteria', 'environmental samples'],\n",
      " 'sequence': RNA\n",
      "---------------------------------------------------------------------\n",
      "Metadata:\n",
      "    'id': 'DQ668905.1'\n",
      "Stats:\n",
      "    length: 363\n",
      "    has gaps: False\n",
      "    has degenerates: False\n",
      "    has definites: True\n",
      "    GC-content: 46.56%\n",
      "---------------------------------------------------------------------\n",
      "0   GUAUGCAACU UACCUUUUAC UAGAGAAUAG CCAAGAGAAA UUUUGAUUAA UGCUCUAUGU\n",
      "60  UCUUAUUUAC UCGCAUGAGU AAAUAAGCAA AGCUCCGGCG GUAAAAGAUG GGCAUGCGUC\n",
      "...\n",
      "300 UAUGGGUUGU AAACUUCUUU UAUACAGGAA GAAACCUUUC CACGUGUGGA AAGCUGACGG\n",
      "360 UAC,\n",
      " 'taxid': 77133}\n"
     ]
    }
   ],
   "source": [
    "options = {'filter-description': set(),\n",
    "           'filter-length': False,\n",
    "           'minimum-length': None,\n",
    "           'maximum-length': None,\n",
    "           'compute-gc': True,\n",
    "           'minimum-gc': None,\n",
    "           'maximum-gc': None}\n",
    "           \n",
    "items = query_rnacentral('tRNA')\n",
    "records = record_parser(items, options)\n",
    "pprint(next(records))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building a better database record processor with the Workflow object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That worked, but can we do better? One large issue with the above function is that there are so many different paths through the function, that it is difficult to validate the individual functions. Additionally, the boilerplate logic for checking options and state is replicated numerous times. In a real world situation, additional features will likely get added in as well further increasing the complexity of the method, and the likelihood of bugs creeping in.\n",
    "\n",
    "The goal of the `Workflow` is to centralize the boilerplate logic, and allow the relatively small blocks of logic that rely on shared state (such as the GC filtering) to be decomposed. To do this, the `Workflow` relies on function decorators: `method` and `requires`. The structure of the object is such that each individual block of code (such as the minimum length filter) can be explicitly unit-tested easily.\n",
    "\n",
    "`method` tells the `Workflow` that this is a method to be executed for every item being operated on. It takes a `priority` parameter that indicates its relative order of execution in the workflow, where a higher value means it will be executed earlier. If a method sets `failed`, the item is considered to have \"failed\" and, by default, none of the other methods are executed (e.g., shortcircuiting like the `continue` above). \n",
    "\n",
    "The second decorator, `requires`, performs the option checking to verify that specific runtime options are set, have particular values (e.g., execute method foo if option X has value Y), and to check on state requirements. If a requirement is not met, the function is not evaluated.\n",
    "\n",
    "Last, the `Workflow` allows the developer to specify a callback to be executed on a successful processing, and one that can be executed on a failed processing.\n",
    "\n",
    "Let's take a look at what this looks like. First, we're going to define two helper functions that simply test for the presence of some piece of state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def has_gc(item):\n",
    "    \"\"\"Test if state has gc computed\"\"\"\n",
    "    return item.get('gc') is not None\n",
    "\n",
    "def has_taxid(item):\n",
    "    \"\"\"Test if state has a valid taxid\"\"\"\n",
    "    return item.get('taxid') is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we're going to define the actual workflow itself. The first method we'll define is a special one that is executed for every record called `initialize_state`. This method is responsible for resetting the notion of `state`, which _should_ be independent for every item. Following this, we'll define the actual methods that'll do the work and decorate on their various requirements. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class RecordParser(Workflow):\n",
    "    def initialize_state(self, item):\n",
    "        \"\"\"Initialize state based on the item being processed\"\"\"\n",
    "        self.state = item.copy()\n",
    "        self.state['gc'] = None\n",
    "        self.state['group'] = None\n",
    "        self.state['lineage'] = None\n",
    "        self.state['id'] = None\n",
    "         \n",
    "    @method(priority=99)\n",
    "    @requires(option='filter-description', values=not_none)  # not_none means any option value except None is valid\n",
    "    def filter_description(self):\n",
    "        \"\"\"Filter records depending on keywords present in the description\"\"\"\n",
    "        acc = self.state['accession']\n",
    "        tokens = {t.strip().lower() for t in acc['description'].split()}\n",
    "       \n",
    "        if self.options['filter-description'].intersection(tokens):\n",
    "            self.failed = True\n",
    "            self.stats['filter-description'] += 1\n",
    "    \n",
    "    @method(priority=95)  # This method has no requirements\n",
    "    def force_to_rna(self):\n",
    "        \"\"\"Convert any thymine to uracil\"\"\" \n",
    "        self.state['sequence'] = self.state['sequence'].replace('T', 'U')\n",
    "    \n",
    "    @method(priority=89)\n",
    "    @requires(option='filter-length', values=True)  # Require 'filter-length' to be True\n",
    "    @requires(option='minimum-length', values=not_none)  # Require that the 'minimum-length' is any value other than None\n",
    "    def minimum_length(self):\n",
    "        if self.state['length'] < self.options['minimum-length']:\n",
    "            self.failed = True\n",
    "            self.stats['minimum-length'] += 1\n",
    "\n",
    "    @method(priority=88)\n",
    "    @requires(option='filter-length', values=True)\n",
    "    @requires(option='maximum-length', values=not_none)  # Require that the 'maximum-length' is any value other than None\n",
    "    def maximum_length(self):\n",
    "        if self.state['length'] > self.options['maximum-length']:\n",
    "            self.failed = True\n",
    "            self.stats['maximum-length'] += 1\n",
    "    \n",
    "    @method(priority=80)\n",
    "    @requires(option='compute-gc', values=True)\n",
    "    def compute_gc(self):\n",
    "        gc = self.state['sequence'].count('G') + self.state['sequence'].count('C')\n",
    "        self.state['gc'] = float(gc) / self.state['length']\n",
    "    \n",
    "    @method(priority=79)\n",
    "    @requires(option='minimum-gc', values=not_none)\n",
    "    @requires(state=has_gc)  # Require that has_gc(self.state) evaluates to True \n",
    "    def minimum_gc(self):\n",
    "        if self.state['gc'] < self.options['minimum-gc']:\n",
    "            self.failed = True\n",
    "            self.stats['minimum-gc'] += 1\n",
    "    \n",
    "    @method(priority=78)\n",
    "    @requires(option='maximum-gc', values=not_none)\n",
    "    @requires(state=has_gc)  # Require that has_gc(self.state) evaluates to True\n",
    "    def maximum_gc(self):\n",
    "        if self.state['gc'] > self.options['maximum-gc']:\n",
    "            self.failed = True\n",
    "            self.stats['maximum-gc'] += 1\n",
    "    \n",
    "    @method(priority=3)\n",
    "    @requires(state=has_taxid)  # Require that has_taxid(self.state) evaluates to True\n",
    "    def fetch_lineage(self):\n",
    "        self.state['lineage'] = query_ncbi_lineage(self.state['taxid'])    \n",
    "    \n",
    "    @method(priority=2)\n",
    "    def assign_to_group(self):\n",
    "        self.state['group'] = self.state['lineage'][-1] if self.state['lineage'] is not None else 'unassigned'\n",
    "        \n",
    "    @method(priority=1)\n",
    "    def cast_to_rnasequence(self):\n",
    "        self.state['sequence'] = RNA(self.state['sequence'], metadata={'id': self.state['id']}, validate=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the workflow itself is complex, each individual block of logic is small making it much easier to test out the components. Integration testing is of course still a very good idea.\n",
    "\n",
    "Now, let's take the object out for a spin. First, we'll instantiate the workflow. We can predefine `state` on instantiation which is useful if it is possible to preallocate memory (instaed of always creating an object in `initialize_state`), but this isn't necessary but a point of optimization. Second, we'll specify the same options we defined in our singular function above. And third, we're going to define a `stats` object to track information about failures (note, `state` is required, `options` is optional, and anything in `kwargs` gets set as a member variable).\n",
    "\n",
    "To get a record out, we pass in our `items` iterator (created from the `query_rna_central` function defined earlier), and simply get the next item. On `__call__`, we can also pass in callbacks which we'll show shortly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'accession': {'citations': 'http://rnacentral.org/api/v1/accession/DQ668905.1:1..363:rfam/citations',\n",
      "               'description': 'uncultured bacterium Bacterial small subunit '\n",
      "                              'ribosomal RNA',\n",
      "               'expert_db_url': 'http://rfam.xfam.org/family/RF00177',\n",
      "               'external_id': 'RF00177',\n",
      "               'gene': '',\n",
      "               'id': 'DQ668905.1:1..363:rfam',\n",
      "               'optional_id': 'SSU_rRNA_bacteria',\n",
      "               'organelle': '',\n",
      "               'product': 'Bacterial small subunit ribosomal RNA',\n",
      "               'rna_type': 'rRNA',\n",
      "               'source_url': '',\n",
      "               'species': 'uncultured bacterium',\n",
      "               'url': 'http://rnacentral.org/api/v1/accession/DQ668905.1:1..363:rfam/info'},\n",
      " 'database': 'Rfam',\n",
      " 'first_seen': '2016-03-15',\n",
      " 'gc': 0.465564738292011,\n",
      " 'group': 'environmental samples',\n",
      " 'id': None,\n",
      " 'is_active': True,\n",
      " 'is_expert_db': False,\n",
      " 'last_seen': '2016-03-15',\n",
      " 'length': 363,\n",
      " 'lineage': ['cellular organisms', 'Bacteria', 'environmental samples'],\n",
      " 'sequence': RNA\n",
      "---------------------------------------------------------------------\n",
      "Metadata:\n",
      "    'id': None\n",
      "Stats:\n",
      "    length: 363\n",
      "    has gaps: False\n",
      "    has degenerates: False\n",
      "    has definites: True\n",
      "    GC-content: 46.56%\n",
      "---------------------------------------------------------------------\n",
      "0   GUAUGCAACU UACCUUUUAC UAGAGAAUAG CCAAGAGAAA UUUUGAUUAA UGCUCUAUGU\n",
      "60  UCUUAUUUAC UCGCAUGAGU AAAUAAGCAA AGCUCCGGCG GUAAAAGAUG GGCAUGCGUC\n",
      "...\n",
      "300 UAUGGGUUGU AAACUUCUUU UAUACAGGAA GAAACCUUUC CACGUGUGGA AAGCUGACGG\n",
      "360 UAC,\n",
      " 'taxid': 77133}\n",
      "Counter()\n"
     ]
    }
   ],
   "source": [
    "records = RecordParser(state=None, options=options, stats=Counter())\n",
    "pprint(next(records(items)))\n",
    "pprint(records.stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, things appear to be working as we'd expect. To finish this notebook off, lets do three more things. First, we'll construct some example callback functions. Second, we're going to change the options to trigger some failures. And last, we're going to enable debug mode to gather additional context.\n",
    "\n",
    "By default, no failure callback is defined as it is assumed that failures are ignored. In addition, there actually is always a success callback but it is by default just a passthrough. Each callback is provided `self`, or the entire workflow object. By default, success simply just returns `self.state`. These callbacks can be handy if you want to automatically stream results to a database, file, or other consumer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We had a failure with DQ668905.1:1..363:silva\n",
      "Counter({'filter-description': 1})\n"
     ]
    }
   ],
   "source": [
    "def fail_callback(wf):\n",
    "    \"\"\"Print out the ID associated with the failure\"\"\"\n",
    "    print(\"We had a failure with %s\" % wf.state['accession']['id'])\n",
    "    \n",
    "def success_callback(wf):\n",
    "    \"\"\"Print out the ID associated with the success\"\"\"\n",
    "    print(\"We successfully processed %s\" % wf.state['accession']['id'])\n",
    "\n",
    "# Let's filter out those pesky ribosomal RNAs\n",
    "options['filter-description'] = {'ribosomal', 'rrna'}\n",
    "\n",
    "records = RecordParser(state=None, options=options, stats=Counter())\n",
    "next(records(items, success_callback=success_callback, fail_callback=fail_callback))\n",
    "pprint(records.stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would be cumbersome to interrogate `stats` everytime something failed in order to figure out why it failed. The workflow assumes that you typically don't care about why something failed specifically (as it also adds a fair bit of overhead), but it is very useful to know when debugging. We can enable this tracking it by setting `debug=True`. To make this more interesting, let's trigger a failure deep in the workflow as the description filter happens first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We had a failure with KC615826.1:1..558:rRNA\n"
     ]
    }
   ],
   "source": [
    "options['filter-description'] = None\n",
    "options['compute-gc'] = True\n",
    "options['minimum-gc'] = 1.0\n",
    "\n",
    "records = RecordParser(state=None, options=options, stats=Counter(), debug=True)\n",
    "next(records(items, success_callback=success_callback, fail_callback=fail_callback))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so we understandably had a failure as the minumum GC was set at 100%. What type of debug detail do we have? First is the trace, which contains `tuple`s of `(\"method_name\", order_of_execution)`.  A low value for the `order_of_execution` means it was executed earlier relative to a higher number. The value is assigned based on the priority, so in the subsequent example, there isn't a method associated with `0` as the highest priority method (`filter_description`) was not executed due to the options we set.\n",
    "\n",
    "What we can see is that the highest number in the execution order is associated with `minimum_gc` indicating that that method failed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{('compute_gc', 4), ('minimum_gc', 5), ('force_to_rna', 1)}\n"
     ]
    }
   ],
   "source": [
    "pprint(records.debug_trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get runtime information (in seconds), where the key in the runtime `dict` corresponds to a `tuple` in the `debug_trace`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{('compute_gc', 4): 1.0728836059570312e-05,\n",
      " ('force_to_rna', 1): 3.0994415283203125e-06,\n",
      " ('minimum_gc', 5): 8.821487426757812e-06}\n"
     ]
    }
   ],
   "source": [
    "pprint(records.debug_runtime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, we have pre- and post-state tracking which lets us examine `state` prior to and following each method execution. These are also keyed by the `tuple`s in the `debug_trace`. For the purposes of example, let's take a look at the `compute_gc` method as that one modifies state, whereas the `minimum_gc` method, though triggering the failure, does not modify the `state` object. And, to minimize cognitive load, let's just look at what happened to the `gc` key. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GC prior to the execution of `compute_gc`\n",
      "None\n",
      "\n",
      "GC following the execution of `compute_gc`\n",
      "0.2616487455197133\n"
     ]
    }
   ],
   "source": [
    "print(\"GC prior to the execution of `compute_gc`\")\n",
    "print(records.debug_pre_state[('compute_gc', 4)]['gc'])\n",
    "\n",
    "print(\"\\nGC following the execution of `compute_gc`\")\n",
    "print(records.debug_post_state[('compute_gc', 4)]['gc'])"
   ]
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}