python_notebooks/1kg_rna_analysis_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GA4GH RNA Quantification API Data Analysis Example\n",
    "This example uses the `rna_quantification_service` to retrieve data for a couple of common RNASeq analyses.  First, a simple spearman correlation between 2 samples will be calculated.  This will be followed by a differential expression calculation using DESeq2 and an expression-based clustering.  For the DESeq2 analysis and clustering this notebook will roughly follow the example RNASeq workflow shown at http://www.bioconductor.org/help/workflows/rnaseqGene/."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize client\n",
    "In this step we create a client object which will be used to communicate with the server.  This example uses the GA4GH 1000 Genomes server.  It will also require that `scipy` and `rpy2` be installed to do the example calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy.stats\n",
    "from ga4gh.client import client\n",
    "import rpy2.robjects as robjects\n",
    "import rpy2.robjects.packages as rpackages\n",
    "local_client = client.HttpClient(\"http://1kgenomes.ga4gh.org\")\n",
    "#local_client = client.HttpClient(\"http://localhost:8000\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search RNA Quantification Sets\n",
    "Begin by retrieving a list of RNA quantification sets in the dataset.  RNA quantification sets are a way to associate a group of related RNA quantifications.  Note that the `dataset_id` is obtained as described in the `1kg_metadata_service` notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dataset_id = local_client.search_datasets().next().id\n",
    "rna_quantification_sets = []\n",
    "for rna_quant_set in local_client.search_rna_quantification_sets(dataset_id=dataset_id):\n",
    "    rna_quantification_sets.append(rna_quant_set.id)\n",
    "    print(\"id: {}\".format(rna_quant_set.id))\n",
    "    print(\"name: {}\\n\".format(rna_quant_set.name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search RNA Quantifications\n",
    "We can now list all of the RNA quantifications in the RNA quantification set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rna_set_id = rna_quantification_sets[0]\n",
    "rna_quantification_ids = []\n",
    "counter = 0\n",
    "max_num_quantifications = 5\n",
    "for rna_quant in local_client.search_rna_quantifications(rna_quantification_set_id=rna_set_id):\n",
    "    if counter < max_num_quantifications:\n",
    "        # In order to make later examples run faster we will store the IDs of the first\n",
    "        # few quantifications returned by the server.\n",
    "        rna_quantification_ids.append(rna_quant.id)\n",
    "    counter += 1\n",
    "    print(\"({}): {}\\n\".format(rna_quant.id, rna_quant.name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select RNA Quantifications\n",
    "The RNA Quantification Set contains several RNA Quantifications.  Before we move on let's verify that the RNA Quantifications to be compared were processed with the same annotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "feature_sets = set()\n",
    "for rna_quant_id in rna_quantification_ids:\n",
    "    for feature_set_id in local_client.get_rna_quantification(\n",
    "            rna_quantification_id=rna_quant_id).feature_set_ids:\n",
    "        feature_sets.add(feature_set_id)\n",
    "print(\"If == 1 we don't have to cull from the list --> {}\".format(len(feature_sets)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search Expression Levels\n",
    "The feature level expression data for each RNA quantification is reported as a set of Expression Levels.  Here is an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def getUnits(unitType):\n",
    "    units = [\"\", \"FPKM\", \"TPM\"]\n",
    "    return units[unitType]\n",
    "\n",
    "counter = 0\n",
    "expression_levels = []\n",
    "for expression in local_client.search_expression_levels(\n",
    "        rna_quantification_id=rna_quantification_ids[0]):\n",
    "    if counter > 5:\n",
    "        break\n",
    "    counter += 1\n",
    "    expression_levels.append(expression)\n",
    "    print(\"Expression Level: {}\".format(expression.name))\n",
    "    print(\" id: {}\".format(expression.id))\n",
    "    print(\" expression: {} {}\".format(expression.expression, getUnits(expression.units)))\n",
    "    print(\" read_count: {}\".format(expression.raw_read_count))\n",
    "    print(\" confidence_interval: {} - {}\\n\".format(\n",
    "            expression.conf_interval_low, expression.conf_interval_high))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Finding Features in the other RNA Quantifications\n",
    "We can also easily examine the expression level of these features in the other RNA Quantifications by specifying the name in the `expressionLevelSearch`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "names = [expression_levels[1].name]\n",
    "for rna_quantification_id in rna_quantification_ids[1:]:\n",
    "    for expression in local_client.search_expression_levels(\n",
    "            rna_quantification_id=rna_quantification_id, names=names):\n",
    "        print(\"RNA Quantification: {}\".format(rna_quantification_id))\n",
    "        print(\"Expression Level: {}\".format(expression.name))\n",
    "        print(\" id: {}\".format(expression.id))\n",
    "        print(\" expression: {} {}\\n\".format(expression.expression, getUnits(expression.units)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Spearman Correlation\n",
    "Now that we can retrieve expression levels for the same feature for all the RNA Quantifications we can calculate the correlation between the Quantifications.  For this we will compute the Spearman Correlation using the `scipy` package.  Further, as this is just an example, we will just calculate the correlation using a small list of features instead of the entire quantification sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def build_expression_dict(rna_quantification_id, max_features=50):\n",
    "    counter = 0\n",
    "    expression_dict = {}\n",
    "    for expression in local_client.search_expression_levels(\n",
    "            rna_quantification_id=rna_quantification_id):\n",
    "        if counter > max_features:\n",
    "            break\n",
    "        counter += 1\n",
    "        expression_dict[expression.name] = expression.expression\n",
    "    return expression_dict\n",
    "\n",
    "\n",
    "expressions_dict_1 = build_expression_dict(rna_quantification_ids[0])\n",
    "featureNames = set(expressions_dict_1.keys())\n",
    "expressions_dict_2 = build_expression_dict(rna_quantification_ids[1])\n",
    "featureNames = featureNames.intersection(set(expressions_dict_2.keys()))\n",
    "sample_1 = []\n",
    "sample_2 = []\n",
    "featureNameList = list(featureNames) # preserve feature order\n",
    "print(\"Comparing {} features\".format(len(featureNameList)))\n",
    "for feature_name in featureNameList:\n",
    "    sample_1.append(expressions_dict_1[feature_name])\n",
    "    sample_2.append(expressions_dict_2[feature_name])\n",
    "\n",
    "scipy.stats.spearmanr(sample_1, sample_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Differential Expression\n",
    "Another common analysis is to look at the differential expression of features between datasets.  This is accomplished by using DESeq2 which is an R package.  This first step will use the `biocLite.R` script to install Bioconductor is needed or update it if already installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "utils = rpackages.importr('utils')\n",
    "utils.chooseBioCmirror(ind=1)\n",
    "robjects.r.source(\"https://bioconductor.org/biocLite.R\")\n",
    "robjects.r.biocLite(robjects.r.c(\"RColorBrewer\", \"pheatmap\", \"DESeq2\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the count matrix\n",
    "DESeq2 operates on a count matrix in which each row contains a feature and each column corresponds to an RnaQuantification.  The values in the matrix are the counts for each feature for each Quantification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def build_expression_dict(rna_quantification_id, max_features=50):\n",
    "    \"\"\"\n",
    "        We are going to rewrite this to return count as well as the quantification\n",
    "        name so that we can build the required matrix.  Also, zero count features\n",
    "        are going to be filtered out.\n",
    "    \"\"\"\n",
    "    counter = 0\n",
    "    expression_dict = {}\n",
    "    quantification = local_client.get_rna_quantification(\n",
    "        rna_quantification_id=rna_quantification_id)\n",
    "    for expression in local_client.search_expression_levels(\n",
    "            rna_quantification_id=rna_quantification_id):\n",
    "        if counter > max_features:\n",
    "            break\n",
    "        counter += 1\n",
    "        if expression.raw_read_count > 0:\n",
    "            expression_dict[expression.name] = expression.raw_read_count\n",
    "    return quantification.name, expression_dict\n",
    "\n",
    "\n",
    "rna_quant_names = []\n",
    "rna_quant_values = []\n",
    "conditions = []\n",
    "for i in range(5):\n",
    "    quant_name, expressions_dict = build_expression_dict(rna_quantification_ids[i])\n",
    "    rna_quant_names.append(quant_name)\n",
    "    rna_quant_values.append(expressions_dict)\n",
    "    conditions.append(\"sample_{}\".format(i+1))\n",
    "    if i == 0:\n",
    "        featureNames = set(expressions_dict.keys())\n",
    "    else:\n",
    "        featureNames = featureNames.intersection(set(expressions_dict.keys()))\n",
    "\n",
    "def get_count_vector(featureList, quantification):\n",
    "    \"\"\"\n",
    "        Extracts counts from the quantification dictionary using featureList as keys.\n",
    "        Returns an robjects.IntVector of counts.\n",
    "    \"\"\"\n",
    "    count_list = [quantification[feature_name] for feature_name in featureList]\n",
    "    return robjects.IntVector(count_list)\n",
    "\n",
    "\n",
    "featureNameList = list(featureNames) # preserve feature order\n",
    "countData = robjects.r.cbind(get_count_vector(featureNameList, rna_quant_values[0]),\n",
    "                             get_count_vector(featureNameList, rna_quant_values[1]))\n",
    "for quantification in rna_quant_values[2:]:\n",
    "    countData = robjects.r.cbind(countData, get_count_vector(featureNameList, quantification))\n",
    "countData.rownames = robjects.StrVector(featureNameList)\n",
    "countData.colnames = robjects.StrVector(rna_quant_names)\n",
    "coldata = robjects.r.cbind(robjects.StrVector(rna_quant_names),\n",
    "                           robjects.StrVector(conditions))\n",
    "coldata.colnames = robjects.StrVector([\"library\", \"sample\"])\n",
    "print(\"matrix created\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ready for DESeq\n",
    "Now that the count matrix and metadata tables are done the next step is to create the `DESeqDataSet` and then run the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "robjects.r.library(\"DESeq2\")\n",
    "design = robjects.Formula(\"~ library\")\n",
    "ddsMat = robjects.r.DESeqDataSetFromMatrix(countData=countData, colData=coldata,\n",
    "                                           design=design)\n",
    "deseq_result = robjects.r.DESeq(ddsMat)\n",
    "# robjects.r.results shows the last result\n",
    "print(robjects.r.results(deseq_result))\n",
    "# can also look any of the results by specifying contrast:\n",
    "print(robjects.r.results(deseq_result, contrast=robjects.r.c(\"library\", rna_quant_names[0],\n",
    "                                                             rna_quant_names[1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Clustering\n",
    "It is often interesting to see how samples cluster with respect to gene expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "robjects.r.X11()\n",
    "robjects.r.library(\"pheatmap\")\n",
    "rld = robjects.r.rlog(deseq_result, blind=False)\n",
    "print(robjects.r.head(robjects.r.assay(rld), 3))\n",
    "robjects.r.pheatmap(robjects.r.assay(rld))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}