datalab/genomics/Getting started with the Genomics API.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Copyright 2015 Google Inc. All rights reserved. -->\n",
    "\n",
    "<!-- Licensed under the Apache License, Version 2.0 (the \"License\"); -->\n",
    "<!-- you may not use this file except in compliance with the License. -->\n",
    "<!-- You may obtain a copy of the License at -->\n",
    "\n",
    "<!--     http://www.apache.org/licenses/LICENSE-2.0 -->\n",
    "\n",
    "<!-- Unless required by applicable law or agreed to in writing, software -->\n",
    "<!-- distributed under the License is distributed on an \"AS IS\" BASIS, -->\n",
    "<!-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -->\n",
    "<!-- See the License for the specific language governing permissions and -->\n",
    "<!-- limitations under the License. -->\n",
    "\n",
    "# Getting started with the Google Genomics API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we'll cover how to make authenticated requests to the [Google Genomics API](https://cloud.google.com/genomics/reference/rest/).\n",
    "\n",
    "----\n",
    "\n",
    "NOTE:\n",
    "\n",
    "* If you're new to notebooks, or want to check out additional samples, check out the full [list](../) of general notebooks.\n",
    "* For additional Genomics samples, check out the full [list](./) of Genomics notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install Python libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be using the [Google Python API client](https://github.com/google/google-api-python-client) for interacting with Genomics API. We can install this library, or any other 3rd-party Python libraries from the [Python Package Index (PyPI)](https://pypi.python.org/pypi) using the `pip` package manager.\n",
    "\n",
    "There are [50+ Google APIs](http://api-python-client-doc.appspot.com/) that you can work against with the Google Python API Client, but we'll focus on the Genomics API in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already up-to-date: google-api-python-client in /usr/local/lib/python2.7/dist-packages\n",
      "Cleaning up...\n"
     ]
    }
   ],
   "source": [
    "!pip install --upgrade google-api-python-client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create an Authenticated Client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we construct a Python object that we can use it to make requests. \n",
    "\n",
    "The following snippet shows how we can authenticate using the service account on the Datalab host.  For more detail about authentication from Python, see [Using OAuth 2.0 for Server to Server Applications](https://developers.google.com/api-client-library/python/auth/service-accounts)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<httplib2.Http at 0x7f6433997a90>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<httplib2.Http at 0x7f6433997a90>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from httplib2 import Http\n",
    "from oauth2client.client import GoogleCredentials\n",
    "credentials = GoogleCredentials.get_application_default()\n",
    "http = Http()\n",
    "credentials.authorize(http)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then we create a client for the Genomics API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from apiclient.discovery import build\n",
    "genomics = build('genomics', 'v1', http=http)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Send a request to the Genomics API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a Python client for the Genomics API, we can access a variety of different resources. For details about each available resource, see the python client [API docs here](https://google-api-client-libraries.appspot.com/documentation/genomics/v1/python/latest/index.html).\n",
    "\n",
    "Using our `genomics` client, we'll demonstrate fetching a Dataset resource by ID (the [1000 Genomes dataset](http://googlegenomics.readthedocs.org/en/latest/use_cases/discover_public_data/1000_genomes.html) in this case).\n",
    "\n",
    "First, we need to construct a request object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "request = genomics.datasets().get(datasetId='10473108253681171589')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll send this request to the Genomics API by calling the `request.execute()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "response = request.execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will need enable the Genomics API for your project if you have not done so previously.  Click on [this link](https://console.developers.google.com/flows/enableapi?apiid=genomics) to enable the API in your project."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The response object returned is simply a Python dictionary. Let's take a look at the properties returned in the response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "projectId => genomics-public-data\n",
      "id => 10473108253681171589\n",
      "createTime => 1970-01-01T00:00:00.000Z\n",
      "name => 1000 Genomes\n"
     ]
    }
   ],
   "source": [
    "for entry in response.items():\n",
    "    print \"%s => %s\" % entry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Success! We can see the name of the specified Dataset and a few other pieces of metadata.\n",
    "\n",
    "Accessing other Genomics API resources will follow this same set of steps. The full [list of available resources within the API is here](https://google-api-client-libraries.appspot.com/documentation/genomics/v1/python/latest/index.html). Each resource has details about the different verbs that can be applied (e.g., [Dataset methods](https://google-api-client-libraries.appspot.com/documentation/genomics/v1/python/latest/genomics_v1.datasets.html))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Access Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this portion of the notebook, we implement [this same example](https://github.com/googlegenomics/getting-started-with-the-api/tree/master/python) implemented as a python script.  First let's define a few constants to use within the examples that follow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dataset_id = '10473108253681171589' # This is the 1000 Genomes dataset ID\n",
    "sample = 'NA12872'\n",
    "reference_name = '22'\n",
    "reference_position = 51003835"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get read bases for a sample at specific a position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First find the read group set ID for the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "request = genomics.readgroupsets().search(\n",
    "  body={'datasetIds': [dataset_id], 'name': sample},\n",
    "  fields='readGroupSets(id)')\n",
    "read_group_sets = request.execute().get('readGroupSets', [])\n",
    "if len(read_group_sets) != 1:\n",
    "  raise Exception('Searching for %s didn\\'t return '\n",
    "                  'the right number of read group sets' % sample)\n",
    "\n",
    "read_group_set_id = read_group_sets[0]['id']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have the read group set ID, lookup the reads at the position in which we are interested."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "request = genomics.reads().search(\n",
    "  body={'readGroupSetIds': [read_group_set_id],\n",
    "        'referenceName': reference_name,\n",
    "        'start': reference_position,\n",
    "        'end': reference_position + 1,\n",
    "        'pageSize': 1024},\n",
    "  fields='alignments(alignment,alignedSequence)')\n",
    "reads = request.execute().get('alignments', [])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we print out the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NA12872 bases on 22 at 51003835 are\n",
      "C: 1\n",
      "G: 13\n"
     ]
    }
   ],
   "source": [
    "# Note: This is simplistic - the cigar should be considered for real code\n",
    "bases = [read['alignedSequence'][\n",
    "           reference_position - int(read['alignment']['position']['position'])]\n",
    "         for read in reads]\n",
    "\n",
    "print '%s bases on %s at %d are' % (sample, reference_name, reference_position)\n",
    "\n",
    "from collections import Counter\n",
    "for base, count in Counter(bases).items():\n",
    "  print '%s: %s' % (base, count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get variants for a sample at specific a position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First find the call set ID for the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "request = genomics.callsets().search(\n",
    "  body={'variantSetIds': [dataset_id], 'name': sample},\n",
    "  fields='callSets(id)')\n",
    "resp = request.execute()\n",
    "call_sets = resp.get('callSets', [])\n",
    "if len(call_sets) != 1:\n",
    "  raise Exception('Searching for %s didn\\'t return '\n",
    "                  'the right number of call sets' % sample)\n",
    "\n",
    "call_set_id = call_sets[0]['id']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have the call set ID, lookup the variants that overlap the position in which we are interested."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "request = genomics.variants().search(\n",
    "  body={'callSetIds': [call_set_id],\n",
    "        'referenceName': reference_name,\n",
    "        'start': reference_position,\n",
    "        'end': reference_position + 1},\n",
    "  fields='variants(names,referenceBases,alternateBases,calls(genotype))')\n",
    "variant = request.execute().get('variants', [])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we print out the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the called genotype is G,G for rs131767\n"
     ]
    }
   ],
   "source": [
    "variant_name = variant['names'][0]\n",
    "genotype = [variant['referenceBases'] if g == 0\n",
    "            else variant['alternateBases'][g - 1]\n",
    "            for g in variant['calls'][0]['genotype']]\n",
    "\n",
    "print 'the called genotype is %s for %s' % (','.join(genotype), variant_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}