results.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sqlite3\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os, sys\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import sqlalchemy\n",
    "from sqlalchemy.ext.automap import automap_base\n",
    "from sqlalchemy.orm import Session\n",
    "from sqlalchemy import create_engine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cd ~/projects/black_spruce/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pd.set_option('display.max_colwidth',10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# conn.close()\n",
    "conn = sqlite3.connect(\"black_spruce.sqlite\")\n",
    "conn.row_factory = sqlite3.Row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts = pd.read_sql('''\n",
    "select unigene_id, sum(assembled_reads) as total \n",
    "from unigene_sample \n",
    "group by unigene_id\n",
    "''', conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples = pd.read_sql(\"select distinct sample_id from unigene_sample\", conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "total_singletons = len(unigene_counts[unigene_counts.total==1])\n",
    "total_singletons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    1102"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def run_select(conn, sql, args):\n",
    "    rows = []\n",
    "    cols = None\n",
    "    for row in conn.execute(sql, args):\n",
    "        if not cols:\n",
    "            cols = list(row.keys())\n",
    "        rows.append(row)\n",
    "    return pd.DataFrame(rows, columns=cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "singletons_by_sample = {}\n",
    "sql = \"select assembled_reads from unigene_sample where unigene_id=? and sample_id=?\"\n",
    "for unigene_id in unigene_counts[unigene_counts.total==1].unigene_id:\n",
    "    for sample_id in samples.sample_id:\n",
    "        if not sample_id in singletons_by_sample:\n",
    "            singletons_by_sample[sample_id] = 0\n",
    "        for count in conn.execute(sql, [unigene_id, sample_id]):\n",
    "            singletons_by_sample[sample_id] += count[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "singletons_by_sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    {u'32C': 189, u'32N': 525, u'40C': 246, u'40N': 142}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## get read stats for samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "run_select(conn, '''\n",
    "select \n",
    "sample_id, \n",
    "sum(raw_reads),\n",
    "sum(phred_reads),\n",
    "sum(length_reads),\n",
    "sum(seqclean_reads)\n",
    "from sample \n",
    "group by sample_id\n",
    "''', [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "run_select(conn,'''\n",
    "select * \n",
    "from sample_stats \n",
    "where key = ?\n",
    "''', [\"seqclean_reads\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_unigene_counts_by_sample(conn):\n",
    "    sql = '''\n",
    "        select unigene.unigene_id, sample_id, mapped_reads, dge, dge_tissue   \n",
    "        from unigene \n",
    "        inner join unigene_sample on unigene.unigene_id = unigene_sample.unigene_id\n",
    "        '''\n",
    "    df = pd.read_sql_query(sql, conn)\n",
    "    return df\n",
    "unigene_counts = get_unigene_counts_by_sample(conn)\n",
    "unigene_counts_all = unigene_counts.copy(deep=True)\n",
    "unigene_tissue = unigene_counts[['unigene_id', 'dge_tissue']].drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples = unigene_counts.sample_id.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts  = unigene_counts [unigene_counts.dge == 1].pivot(\"unigene_id\", \"sample_id\")[\"mapped_reads\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts['expression'] = unigene_counts.apply(np.sum, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts = unigene_counts.sort('expression', ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_unigene_annotations(conn):\n",
    "    sql = \"select * from unigene_anno\"\n",
    "    return pd.read_sql_query(sql, conn)\n",
    "unigene_anno = get_unigene_annotations(conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ontologies = unigene_anno.ontology.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_ontology(row, ont, anno):\n",
    "    anno_df = anno[(anno.unigene_id == row.name) & (anno.ontology == ont)]\n",
    "    if len(anno_df):\n",
    "        return';'.join(anno_df.anno)\n",
    "    \n",
    "for ontology in ontologies:\n",
    "    unigene_counts[ontology] = unigene_counts.apply(get_ontology, args=(ontology, unigene_anno), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_tissue.index = unigene_tissue.unigene_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "unigene_counts_tissue = unigene_counts.join(unigene_tissue)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts_tissue[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_blast_desc_for_unigenes(conn):\n",
    "    sql = \"select * from unigene_blast\"\n",
    "    return pd.read_sql_query(sql, conn)\n",
    "unigene_blast = get_blast_desc_for_unigenes(conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_blast_desc(row):\n",
    "    blast_df = unigene_blast[unigene_blast.unigene_id==row.name]\n",
    "    if len(blast_df):\n",
    "        return ';'.join(blast_df.desc)\n",
    "\n",
    "unigene_counts_tissue['blastx'] = unigene_counts_tissue.apply(get_blast_desc, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_counts_tissue[0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##A Venn diagram showing the number of common and differentially expressed genes between cambium and needle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unigene_mapped_totals = {}\n",
    "for row in conn.execute('''select unigene.unigene_id, sample_id, mapped_reads, dge  \n",
    "from unigene_sample\n",
    "inner join unigene on \n",
    "unigene.unigene_id = unigene_sample.unigene_id\n",
    "'''):\n",
    "    if not row[0] in unigene_mapped_totals:\n",
    "        unigene_mapped_totals[row[0]] = {\"C\":0, \"N\":0}\n",
    "    tissue = row[1][-1]\n",
    "    dge = row[-1]\n",
    "    unigene_mapped_totals[row[0]][tissue] += row[2]\n",
    "    unigene_mapped_totals[row[0]][\"dge\"] = dge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "min_count = 2\n",
    "tissue_keys = [\"C\", \"N\"]\n",
    "tissue_totals = {\"C\":0, \"N\":0, \"B\":0, \"X\":0, \"C_DGE\":0, \"N_DGE\":0, \"B_DGE\":0, \"B_DGE_C\":0, \"B_DGE_N\":0}\n",
    "unigene_keys = sorted(list(unigene_mapped_totals.keys()))\n",
    "for k in unigene_keys:\n",
    "    v = unigene_mapped_totals[k]\n",
    "    print v\n",
    "    if v[\"C\"] and v[\"N\"] and (v[\"C\"] + v[\"N\"] > min_count):\n",
    "        key = \"B\"\n",
    "        dge_key = \"B_DGE_C\"\n",
    "        if v[\"N\"] > v[\"C\"]:\n",
    "            dge_key = \"B_DGE_N\"\n",
    "        if v['dge']:\n",
    "            tissue_totals[dge_key] += 1\n",
    "            tissue_totals[\"B_DGE\"] += 1\n",
    "    elif v[\"C\"] and not v[\"N\"] and v[\"C\"] > min_count:\n",
    "        key = \"C\"\n",
    "        if v['dge']:\n",
    "            tissue_totals[\"C_DGE\"] += 1\n",
    "    elif v[\"N\"] and not v[\"C\"] and v[\"N\"] > min_count:\n",
    "        key = \"N\"\n",
    "        if v['dge']:\n",
    "            tissue_totals[\"N_DGE\"] += 1\n",
    "    else:\n",
    "        key = \"X\"\n",
    "        if v['dge']:\n",
    "            print v\n",
    "    tissue_totals[key] += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tissue_totals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib_venn import venn2\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(1)\n",
    "plt.subplot(211)\n",
    "v = venn2(subsets = (tissue_totals[\"C\"], tissue_totals[\"N\"], tissue_totals[\"B\"]),\n",
    "      set_labels = (\"Cambium\", \"Needle\"))\n",
    "plt.subplot(212)\n",
    "v = venn2(subsets = (tissue_totals[\"C_DGE\"], tissue_totals[\"N_DGE\"], tissue_totals[\"B_DGE\"]),\n",
    "      set_labels = (\"Cambium\", \"Needle\"))\n",
    "\n",
    "plt.savefig(\"manuscript/figures/unigene_venn.svg\")\n",
    "plt.savefig(\"manuscript/figures/unigene_venn.pdf\")\n",
    "plt.savefig(\"manuscript/figures/unigene_venn.png\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Base=automap_base()\n",
    "engine = create_engine(\"sqlite:////Users/chris/projects/black_spruce/black_spruce.sqlite\")\n",
    "session = Session(engine)\n",
    "Base.prepare(engine, reflect=True)\n",
    "for t in Base.metadata.tables:\n",
    "    print \"%s = Base.classes.%s\" % (t.capitalize(), t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Unigene_sample = Base.classes.unigene_sample\n",
    "Unigene_anno = Base.classes.unigene_anno\n",
    "Sample = Base.classes.sample\n",
    "Sample_stats = Base.classes.sample_stats\n",
    "Unigene_go = Base.classes.unigene_go\n",
    "Unigene_blast = Base.classes.unigene_blast\n",
    "Go = Base.classes.go\n",
    "Unigene = Base.classes.unigene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "no_hits = []\n",
    "for x in session.query(Unigene).outerjoin(Unigene_blast):\n",
    "    if len(x.unigene_blast_collection) < 1:\n",
    "        no_hits.append(x.unigene_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print len(no_hits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": "0",
   "metadata": {
    "collapsed": false
   },
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}