Smith_2017_Figure7E.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I need to install qiime on the cluster, then I go:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CentOS Linux release 7.2.1511 (Core) \r\n",
      "NAME=\"CentOS Linux\"\r\n",
      "VERSION=\"7 (Core)\"\r\n",
      "ID=\"centos\"\r\n",
      "ID_LIKE=\"rhel fedora\"\r\n",
      "VERSION_ID=\"7\"\r\n",
      "PRETTY_NAME=\"CentOS Linux 7 (Core)\"\r\n",
      "ANSI_COLOR=\"0;31\"\r\n",
      "CPE_NAME=\"cpe:/o:centos:centos:7\"\r\n",
      "HOME_URL=\"https://www.centos.org/\"\r\n",
      "BUG_REPORT_URL=\"https://bugs.centos.org/\"\r\n",
      "\r\n",
      "CENTOS_MANTISBT_PROJECT=\"CentOS-7\"\r\n",
      "CENTOS_MANTISBT_PROJECT_VERSION=\"7\"\r\n",
      "REDHAT_SUPPORT_PRODUCT=\"centos\"\r\n",
      "REDHAT_SUPPORT_PRODUCT_VERSION=\"7\"\r\n",
      "\r\n",
      "CentOS Linux release 7.2.1511 (Core) \r\n",
      "CentOS Linux release 7.2.1511 (Core) \r\n"
     ]
    }
   ],
   "source": [
    "!cat /etc/*-release "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To install QIIME, I am following instructions from [here](https://gist.github.com/alanorth/100dd95f809223ae0100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!module load Python\n",
    "#!module load R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting qiime\n",
      "  Using cached qiime-1.9.1.tar.gz\n",
      "Requirement already satisfied (use --upgrade to upgrade): numpy>=1.9.0 in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from qiime)\n",
      "Collecting scipy>=0.14.0 (from qiime)\n",
      "  Using cached scipy-0.17.0-cp27-cp27m-manylinux1_x86_64.whl\n",
      "Collecting cogent==1.5.3 (from qiime)\n",
      "  Using cached cogent-1.5.3.tgz\n",
      "Collecting natsort<4.0.0 (from qiime)\n",
      "  Using cached natsort-3.5.6-py2.py3-none-any.whl\n",
      "Requirement already satisfied (use --upgrade to upgrade): matplotlib!=1.4.2,>=1.1.0 in /beegfs/common/software/Python/2.7.9/lib/python2.7/site-packages (from qiime)\n",
      "Collecting pynast==1.2.2 (from qiime)\n",
      "  Using cached pynast-1.2.2.tar.gz\n",
      "Collecting qcli<0.2.0,>=0.1.1 (from qiime)\n",
      "  Using cached qcli-0.1.1.tar.gz\n",
      "Collecting gdata (from qiime)\n",
      "  Using cached gdata-2.0.18.tar.gz\n",
      "Collecting biom-format<2.2.0,>=2.1.4 (from qiime)\n",
      "  Downloading biom-format-2.1.5.tar.gz (307kB)\n",
      "\u001b[K    100% |################################| 317kB 1.7MB/s \n",
      "\u001b[?25hCollecting emperor<1.0.0,>=0.9.51 (from qiime)\n",
      "  Using cached emperor-0.9.51.tar.gz\n",
      "Collecting scikit-bio<0.3.0,>=0.2.3 (from qiime)\n",
      "  Using cached scikit-bio-0.2.3.tar.gz\n",
      "Collecting burrito-fillings<0.2.0,>=0.1.1 (from qiime)\n",
      "  Downloading burrito-fillings-0.1.1.tar.gz (214kB)\n",
      "\u001b[K    100% |################################| 215kB 1.7MB/s \n",
      "\u001b[?25hRequirement already satisfied (use --upgrade to upgrade): pandas>=0.13.1 in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from qiime)\n",
      "Collecting burrito<1.0.0,>=0.9.1 (from qiime)\n",
      "  Using cached burrito-0.9.1.tar.gz\n",
      "Collecting qiime-default-reference<0.2.0,>=0.1.2 (from qiime)\n",
      "  Using cached qiime-default-reference-0.1.3.tar.gz\n",
      "Requirement already satisfied (use --upgrade to upgrade): six>=1.4 in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): python-dateutil in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): pytz in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): pyparsing>=1.5.6 in /beegfs/common/software/Python/2.7.9/lib/python2.7/site-packages (from matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): nose>=0.11.1 in /beegfs/common/software/Python/2.7.9/lib/python2.7/site-packages (from matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): mock in /beegfs/common/software/Python/2.7.9/lib/python2.7/site-packages (from matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Collecting click (from biom-format<2.2.0,>=2.1.4->qiime)\n",
      "  Using cached click-6.6.tar.gz\n",
      "Collecting future>=0.14.3 (from biom-format<2.2.0,>=2.1.4->qiime)\n",
      "  Using cached future-0.15.2.tar.gz\n",
      "Collecting pyqi (from biom-format<2.2.0,>=2.1.4->qiime)\n",
      "  Using cached pyqi-0.3.2.tar.gz\n",
      "Requirement already satisfied (use --upgrade to upgrade): IPython in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): funcsigs in /beegfs/common/software/Python/2.7.9/lib/python2.7/site-packages (from mock->matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): pbr>=0.11 in /beegfs/common/software/Python/2.7.9/lib/python2.7/site-packages (from mock->matplotlib!=1.4.2,>=1.1.0->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): decorator in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): simplegeneric>0.8 in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): pexpect in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): traitlets in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): pickleshare in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): ptyprocess>=0.5 in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from pexpect->IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): ipython-genutils in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from traitlets->IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Requirement already satisfied (use --upgrade to upgrade): path.py>=6.2 in /beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages (from pickleshare->IPython->scikit-bio<0.3.0,>=0.2.3->qiime)\n",
      "Installing collected packages: scipy, cogent, natsort, pynast, qcli, gdata, click, future, pyqi, biom-format, scikit-bio, emperor, burrito, burrito-fillings, qiime-default-reference, qiime\n",
      "  Running setup.py install for cogent ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \bdone\n",
      "\u001b[?25h  Running setup.py install for pynast ... \u001b[?25l-\b \b\\\b \bdone\n",
      "\u001b[?25h  Running setup.py install for qcli ... \u001b[?25l-\b \b\\\b \bdone\n",
      "\u001b[?25h  Running setup.py install for gdata ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \bdone\n",
      "\u001b[?25h  Running setup.py install for click ... \u001b[?25l-\b \b\\\b \b|\b \bdone\n",
      "\u001b[?25h  Running setup.py install for future ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \bdone\n",
      "\u001b[?25h  Running setup.py install for pyqi ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \bdone\n",
      "\u001b[?25h  Running setup.py install for biom-format ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \bdone\n",
      "\u001b[?25h  Running setup.py install for scikit-bio ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \bdone\n",
      "\u001b[?25h  Running setup.py install for emperor ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \bdone\n",
      "\u001b[?25h  Running setup.py install for burrito ... \u001b[?25l-\b \b\\\b \bdone\n",
      "\u001b[?25h  Running setup.py install for burrito-fillings ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \bdone\n",
      "\u001b[?25h  Running setup.py install for qiime-default-reference ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \bdone\n",
      "\u001b[?25h  Running setup.py install for qiime ... \u001b[?25l-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \b|\b \b/\b \b-\b \b\\\b \bdone\n",
      "\u001b[?25hSuccessfully installed biom-format-2.1.5 burrito-0.9.1 burrito-fillings-0.1.1 click-6.6 cogent-1.5.3 emperor-0.9.51 future-0.15.2 gdata-2.0.18 natsort-3.5.6 pynast-1.2.2 pyqi-0.3.2 qcli-0.1.1 qiime-1.9.1 qiime-default-reference-0.1.3 scikit-bio-0.2.3 scipy-0.17.0\n"
     ]
    }
   ],
   "source": [
    "!pip install --user qiime -I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "System information\n",
      "==================\n",
      "         Platform:\tlinux2\n",
      "   Python version:\t2.7.9 (default, Apr 10 2015, 15:04:25)  [GCC 4.8.2 20140120 (Red Hat 4.8.2-16)]\n",
      "Python executable:\t/software/Python/2.7.9/bin/python\n",
      "\n",
      "QIIME default reference information\n",
      "===================================\n",
      "For details on what files are used as QIIME's default references, see here:\n",
      " https://github.com/biocore/qiime-default-reference/releases/tag/0.1.3\n",
      "\n",
      "Dependency versions\n",
      "===================\n",
      "          QIIME library version:\t1.9.1\n",
      "           QIIME script version:\t1.9.1\n",
      "qiime-default-reference version:\t0.1.3\n",
      "                  NumPy version:\t1.10.4\n",
      "                  SciPy version:\t0.17.0\n",
      "                 pandas version:\t0.18.0\n",
      "             matplotlib version:\t1.4.3\n",
      "            biom-format version:\t2.1.5\n",
      "                   h5py version:\tNot installed.\n",
      "                   qcli version:\t0.1.1\n",
      "                   pyqi version:\t0.3.2\n",
      "             scikit-bio version:\t0.2.3\n",
      "                 PyNAST version:\t1.2.2\n",
      "                Emperor version:\t0.9.51\n",
      "                burrito version:\t0.9.1\n",
      "       burrito-fillings version:\t0.1.1\n",
      "              sortmerna version:\tSortMeRNA version 2.0, 29/11/2014\n",
      "              sumaclust version:\tSUMACLUST Version 1.0.00\n",
      "                  swarm version:\tSwarm 1.2.19 [Apr 20 2016 13:38:51]\n",
      "                          gdata:\tInstalled.\n",
      "\n",
      "QIIME config values\n",
      "===================\n",
      "For definitions of these settings and to learn how to configure QIIME, see here:\n",
      " http://qiime.org/install/qiime_config.html\n",
      " http://qiime.org/tutorials/parallel_qiime.html\n",
      "\n",
      "                     blastmat_dir:\tNone\n",
      "      pick_otus_reference_seqs_fp:\t/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta\n",
      "                         sc_queue:\tall.q\n",
      "      topiaryexplorer_project_dir:\tNone\n",
      "     pynast_template_alignment_fp:\t/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set_aligned/85_otus.pynast.fasta\n",
      "                  cluster_jobs_fp:\tstart_parallel_jobs.py\n",
      "pynast_template_alignment_blastdb:\tNone\n",
      "assign_taxonomy_reference_seqs_fp:\t/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta\n",
      "                     torque_queue:\tfriendlyq\n",
      "                    jobs_to_start:\t1\n",
      "                       slurm_time:\tNone\n",
      "            denoiser_min_per_core:\t50\n",
      "assign_taxonomy_id_to_taxonomy_fp:\t/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/taxonomy/97_otu_taxonomy.txt\n",
      "                         temp_dir:\t/tmp/\n",
      "                     slurm_memory:\tNone\n",
      "                      slurm_queue:\tNone\n",
      "                      blastall_fp:\tblastall\n",
      "                 seconds_to_sleep:\t1\n",
      "\n",
      "QIIME base install test results\n",
      "===============================\n",
      ".........\n",
      "----------------------------------------------------------------------\n",
      "Ran 9 tests in 0.145s\n",
      "\n",
      "OK\n"
     ]
    }
   ],
   "source": [
    "!print_qiime_config.py -t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goal: To match IDs to OTUs - this will enable me to select OTU differentially enriched based on the LEfSE analysis,  \n",
    "and return transcripts, based on the reference correlation performed here: /Volumes/group_dv/personal/DValenzano/papers/uBiome0/data/regr.csv  \n",
    "First, I take the OTU genus reference table, normalise it and save it as a csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#!export PATH=/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/qiime/:$PATH\n",
    "#!export PATH=/export/apps/qiime/1.8.0/bin:$PATH\n",
    "#!export PATH=$PATH:/software/R/3.2.4/bin/\n",
    "#!export PATH=$PATH:/software/R/3.2.4/bin\n",
    "#!echo $PATH\n",
    "#ls ~/.bashrc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To install all the required R packages, follow instructions from http://qiime.org/install/install.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#!module load R\n",
    "#!which R\n",
    "get_ipython().system(u'normalize_table.py -i /beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6.biom -a CSS -o /beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.biom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This file I just normalized has the OTUs **at the genus level** for all the individuals for which we also have the rnaseq data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "get_ipython().system(u'biom convert -i /beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.biom -o /beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.txt --to-tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "get_ipython().system(u'biom convert -i /beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.biom -o /beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.txt --to-tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "get_ipython().system(u'biom convert -i /beegfs/group_dv/home/DValenzano/uBiome/data/16wk_otu_table.biom -o /beegfs/group_dv/home/DValenzano/uBiome/data/16wk_otu_table.txt --to-tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t = open('/beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.txt', 'rU').read()\n",
    "t = t.replace('\\t',',')\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.csv', 'w')\n",
    "z.write(t)\n",
    "z.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#t = open('/beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.csv', 'rU').read()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now I need to generate a new correlation between fpkm and OTU, individual by individual  \n",
    "The OTU Table is [this](/beegfs/group_dv/home/DValenzano/uBiome/data/Description_otu_table_sorted_L6n.csv), which now needs to be transposed and sorted  \n",
    "I am following the same guidelines as in [here](https://github.com/dvalenzano/Interactive-Sessions/blob/master/11-Mar-2016.py)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t1 = t.replace('_','.').replace('YI', 'WT')\n",
    "t2 = t1.split('\\n')[1:]\n",
    "t3 = [ i.split(',') for i in t2 ]\n",
    "t4 = [t3[0]]+[i for i in t3[1:] if sum(map(float, (i[1:]))) != 0.0] #this removes 0-only OTUs\n",
    "tt = zip(*t4)\n",
    "tT = [ list(i) for i in tt ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "u = open('/beegfs/group_dv/home/DValenzano/uBiome/data/rnaseq_fpkm_genes.csv', 'rU').read()\n",
    "u2 = u.split('\\n')\n",
    "u3 = [ i.split(',') for i in u2]\n",
    "u4 = [u3[0]]+[i for i in u3[1:] if sum(map(float, (i[1:]))) != 0.0] #this removes 0-only transcripts\n",
    "ut = zip(*u4)\n",
    "uT = [ list(i) for i in ut ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/ipykernel/__main__.py:1: DeprecationWarning: the sets module is deprecated\n",
      "  if __name__ == '__main__':\n"
     ]
    }
   ],
   "source": [
    "from sets import Set\n",
    "uh = [i[0] for i in uT[1:]]\n",
    "th = [i[0] for i in tT[1:]]\n",
    "h = list(Set(uh) & Set(th))\n",
    "tT1 = [tT[0]] + [i for i in tT if i[0] in h]\n",
    "uT1 = [uT[0]] + [i for i in uT if i[0] in h]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "uTs = [uT1[0]] + sorted(uT1[1:], key = lambda x: (x[0]))\n",
    "tTs = [tT1[0]] + sorted(tT1[1:], key = lambda x: (x[0]))\n",
    "\n",
    "#[i[0] for i in uTs][1:] == [i[0] for i in tTs][1:] # this tests whether the two matrices have headers in same order\n",
    "\n",
    "uj = ','.join([','.join(i)+'\\n' for i in uTs]).replace('\\n,','\\n')[:-1]\n",
    "tj = ','.join([','.join(i)+'\\n' for i in tTs]).replace('\\n,','\\n')[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zu = open('/beegfs/group_dv/home/DValenzano/uBiome/data/rnaseq_fpkm_genes_sorted.csv', 'w')\n",
    "zu.write(uj)\n",
    "zu.close()\n",
    "\n",
    "zt = open('/beegfs/group_dv/home/DValenzano/uBiome/data/rnaseq_merged_otu_table_sorted.csv', 'w')\n",
    "zt.write(tj)\n",
    "zt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[i.split(',')[0] for i in tj.split('\\n')[1:-1]] == [i.split(',')[0] for i in uj.split('\\n')[1:-1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "There is a difference from the previous analysis (March 17-18 on README) since 3 samples were removed due to low depth  \n",
    "Good - Now continuing from [here](https://github.com/dvalenzano/Interactive-Sessions/blob/master/15-Mar-2016.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats.stats import pearsonr\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "from scipy import stats, integrate\n",
    "import matplotlib.pyplot as plt\n",
    "#import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rna = uj\n",
    "otu = tj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rnal = [i.split(',') for i in rna.split('\\n')] #that's a \"l\" as in \"L\"\n",
    "otul = [i.split(',') for i in otu.split('\\n')]\n",
    "\n",
    "rnat = [ list(i) for i in zip(*rnal)]\n",
    "otut = [list(i) for i in zip(*otul) ]\n",
    "\n",
    "rnat_z = [ map(float, i[1:]) for i in rnat[1:] if i.count('0') < 7]\n",
    "otut_z = [ map(float, i[1:]) for i in otut[1:] if i.count('0.0') < 7]\n",
    "\n",
    "\n",
    "# In[10]:\n",
    "\n",
    "rna_key = [i[0] for i in rnat[1:] if i.count('0') < 7]\n",
    "otu_key = [ i[0] for i in otut[1:] if i.count('0') < 7]\n",
    "\n",
    "rna_value = rnat_z\n",
    "otu_value = otut_z\n",
    "\n",
    "rna_d = dict(zip(rna_key, rna_value))\n",
    "otu_d = dict(zip(otu_key, otu_value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def loop_one(m1,otu): # this returns all the correlation coefficients between m1 and a given otu (otu)\n",
    "    ls = []\n",
    "    for i in m1:\n",
    "        ls.append(stats.linregress(i, otu))    #this function returns slope, intercept, r_value, p_value, std_err\n",
    "    return ','.join(map(str, ls))+'\\n'\n",
    "\n",
    "def loop_two(m1, m2):\n",
    "    ls = []\n",
    "    for j in m2:\n",
    "        ls.append(loop_one(m1, j))\n",
    "    return ','.join(ls).replace('\\n,','\\n')\n",
    "\n",
    "fpkmvsotu = loop_two(rnat_z, otut_z)\n",
    "fpkmvsotu1 = fpkmvsotu.replace('LinregressResult', '').replace('slope=','').replace('intercept=', '').replace('rvalue=', '').replace('pvalue=', '').replace('stderr=', '') \n",
    "\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/corr.csv','w')\n",
    "z.write(fpkmvsotu1)\n",
    "z.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class reg(object):\n",
    "    \n",
    "    \"\"\"\n",
    "    This class uses as input a linear regression file obtained by the loop_two function and returns \n",
    "    regression tables and plots, given specific parameteres, such as R-squared value, p-value, OTU or \n",
    "    transcript\n",
    "    \"\"\"\n",
    "            \n",
    "#     version = 0.1\n",
    "    \n",
    "#     def __init__(self, inp):\n",
    "#         self.inp = inp \n",
    "#         self.inp2 = [ i.split('LinregressResult')[1:] for i in self.inp.split('\\n')[:-1]] \n",
    "#         def loop1(inp, n):\n",
    "#             ls = []\n",
    "#             for i in inp:\n",
    "#                 ls.append(i.split(',')[n].split('=')[1])\n",
    "#             return ','.join(ls).replace(')', '')\n",
    "#         def loop2(inp, n):\n",
    "#             ls = []\n",
    "#             for i in inp:\n",
    "#                 ls.append(loop1(i, n))\n",
    "#             return ls\n",
    "#         self.slope =  [map(float, i.split(',')) for i in loop2(self.inp2, 0)]\n",
    "#         self.intercept = [map(float, i.split(',')) for i in loop2(self.inp2, 1)]\n",
    "#         self.rvalue = [map(float, i.split(',')) for i in loop2(self.inp2, 2)]\n",
    "#         self.rsquare = [ map(lambda x: x**2, i) for i in self.rvalue]\n",
    "#         self.pvalue = [map(float, i.split(',')) for i in loop2(self.inp2, 3)]\n",
    "#         self.stderr0 = [map(float, i.split(',')) for i in loop2(self.inp2, 4)  ]  \n",
    "\n",
    "    version = 0.2\n",
    "    \n",
    "    def __init__(self, inp):\n",
    "        self.inp = inp \n",
    "        self.inp2 = [ i.split('(') for i in self.inp.split('\\n')[:-1]] \n",
    "    \n",
    "        def loop1(inp, n):\n",
    "            ls = []\n",
    "            for i in inp[1:]:\n",
    "                ls.append(i.split(',')[n])\n",
    "            return ','.join(ls).replace(')', '')\n",
    "        def loop2(inp, n):\n",
    "            ls = []\n",
    "            for i in inp:\n",
    "                ls.append(loop1(i, n))\n",
    "            return ls\n",
    "       \n",
    "        self.slope =  [map(float, i.split(',')) for i in loop2(self.inp2, 0)]\n",
    "        self.intercept = [map(float, i.split(',')) for i in loop2(self.inp2, 1)]\n",
    "        self.rvalue = [map(float, i.split(',')) for i in loop2(self.inp2, 2)]\n",
    "        self.rsquare = [ map(lambda x: x**2, i) for i in self.rvalue]\n",
    "        self.pvalue = [map(float, i.split(',')) for i in loop2(self.inp2, 3)]\n",
    "        self.stderr = [map(float, i.split(',')) for i in loop2(self.inp2, 4)  ]   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def filt(inp, thr, sign): #filter p-values, r-values and so forth based on a chosen threshold\n",
    "    ar = np.array(inp)\n",
    "    ls = []\n",
    "    if sign == '>':\n",
    "        ls = zip(*np.where(ar > thr))\n",
    "    elif sign == '<':\n",
    "        ls = zip(*np.where(ar < thr))\n",
    "    else: \n",
    "        print \"Error: check your input file\"\n",
    "    return ls\n",
    "\n",
    "\n",
    "def filt_(inp, thr, sign): #filter p-values, r-values and so forth based on a chosen threshold\n",
    "    ar = np.array(inp)\n",
    "    li = []\n",
    "    lp = []\n",
    "    if sign == '>':\n",
    "        lp = np.where(ar > thr)\n",
    "    elif sign == '<':\n",
    "        lp = np.where(ar < thr)\n",
    "    else: \n",
    "        print \"Error: check your input file\"\n",
    "    p1_l = [ list(i) for i in lp]\n",
    "    la = list(ar[p1_l])\n",
    "    p2 = p1_l+[la]\n",
    "    p3 = [ list(i) for i in zip(*p2)]\n",
    "    p4 = sorted(p3, key=lambda x: x[2])\n",
    "    return ','.join([ ','.join(map(str, i))+'\\n' for i in p4]).replace('\\n,','\\n')[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#fpkmvsotu = open('/beegfs/group_dv/home/DValenzano/uBiome/data/corr.csv', 'rU').read()\n",
    "#fpkmvsotu1 = fpkmvsotu.replace('LinregressResult', '').replace('slope=','').replace('intercept=', '').replace('rvalue=', '').replace('pvalue=', '').replace('stderr=', '') \n",
    "\n",
    "rf = reg(fpkmvsotu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rpval = filt_(rf.pvalue, 0.001, '<')\n",
    "P2 = [otu_key[int(i.split(',')[0])]+','+rna_key[int(i.split(',')[1])] +','+i.split(',')[2] for i in rpval.split('\\n') ]\n",
    "P3 = ','.join([ i+'\\n' for i in P2]).replace('\\n,','\\n')[:-1]\n",
    "len(P3.split('\\n'))\n",
    "\n",
    "zpvalue = 'OTU,transcr,value\\n'+','.join([i+'\\n'for i in P3.split('\\n')[:100]]).replace('\\n,','\\n')\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/TOP_otu-tr_pvalue.csv', 'w')\n",
    "z.write(zpvalue)\n",
    "z.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Below the instructions to plot "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Rf = reg(fpkmvsotu1)\n",
    "Rpval = filt_(Rf.pvalue, 1.0, '<') #all the pvals \n",
    "Rrsq = filt_(Rf.rsquare, 0.0, '>') #all the rsquare "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This won't work from the server - for now\n",
    "# Xp = []\n",
    "# Yp = []\n",
    "# Otu = []\n",
    "# Rna = []\n",
    "# for i in range(12):\n",
    "#     Otu.append(Rpval.split('\\n')[i].split(',')[0])\n",
    "#     Rna.append(Rpval.split('\\n')[i].split(',')[1])\n",
    "#     Xp.append(otu_d[otu_key[int(Rpval.split('\\n')[i].split(',')[0])]])\n",
    "#     Yp.append(rna_d[rna_key[int(Rpval.split('\\n')[i].split(',')[1])]])\n",
    "\n",
    "# rng = zip(range(10), range(331,340))\n",
    "\n",
    "# fig = plt.figure()\n",
    "# fig.suptitle('Top 9 regressions based on p value ', fontsize=14, fontweight='bold')\n",
    "\n",
    "# for i,j in rng:\n",
    "#     ax = fig.add_subplot(j)\n",
    "#     fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "#     #ax.set_title('axes title')\n",
    "#     ax.set_xlabel(Otu[i])\n",
    "#     ax.set_ylabel(Rna[i])\n",
    "#     plt.scatter(Xp[i],Yp[i])\n",
    "    \n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This won't work from the server - for now\n",
    "# Xr = []\n",
    "# Yr = []\n",
    "# Otu = []\n",
    "# Rna = []\n",
    "# for i in range((len(Rrsq.split('\\n')) - 12),len(Rrsq.split('\\n'))):\n",
    "#     Otu.append(Rrsq.split('\\n')[i].split(',')[0])\n",
    "#     Rna.append(Rrsq.split('\\n')[i].split(',')[1])\n",
    "#     Xr.append(otu_d[otu_key[int(Rrsq.split('\\n')[i].split(',')[0])]])\n",
    "#     Yr.append(rna_d[rna_key[int(Rrsq.split('\\n')[i].split(',')[1])]])\n",
    "    \n",
    "# fig = plt.figure()\n",
    "# fig.suptitle('Top 9 regressions based on R2 ', fontsize=14, fontweight='bold')\n",
    "\n",
    "# for i,j in rng:\n",
    "#     ax = fig.add_subplot(j)\n",
    "#     fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "#     #ax.set_title('axes title')\n",
    "#     ax.set_xlabel(Otu[i])\n",
    "#     ax.set_ylabel(Rna[i])\n",
    "#     plt.scatter(Xr[i],Yr[i])\n",
    "    \n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Below what I ran on Apr 26 2016 **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fpkmvsotu0 = open('/beegfs/group_dv/home/DValenzano/uBiome/data/corr.csv', 'rU').read()\n",
    "\n",
    "uj = open('/beegfs/group_dv/home/DValenzano/uBiome/data/rnaseq_fpkm_genes_sorted.csv', 'rU').read()\n",
    "\n",
    "tj = open('/beegfs/group_dv/home/DValenzano/uBiome/data/rnaseq_merged_otu_table_sorted.csv','rU').read()\n",
    "\n",
    "rna = uj\n",
    "otu = tj\n",
    "\n",
    "rnal = [i.split(',') for i in rna.split('\\n')] #that's a \"l\" as in \"L\"\n",
    "otul = [i.split(',') for i in otu.split('\\n')]\n",
    "\n",
    "rnat = [ list(i) for i in zip(*rnal)]\n",
    "otut = [list(i) for i in zip(*otul) ]\n",
    "\n",
    "rnat_z = [ map(float, i[1:]) for i in rnat[1:] if i.count('0') < 7]\n",
    "otut_z = [ map(float, i[1:]) for i in otut[1:] if i.count('0.0') < 7]\n",
    "\n",
    "\n",
    "# In[10]:\n",
    "\n",
    "rna_key = [i[0] for i in rnat[1:]  if i.count('0') < 7]\n",
    "otu_key = [ i[0] for i in otut[1:] if i.count('0.0') < 7]\n",
    "\n",
    "rna_value = rnat_z\n",
    "otu_value = otut_z\n",
    "\n",
    "rna_d = dict(zip(rna_key, rna_value))\n",
    "otu_d = dict(zip(otu_key, otu_value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fpkmvsotu1 = fpkmvsotu0.replace('LinregressResult', '').replace('slope=','').replace('intercept=', '').replace('rvalue=', '').replace('pvalue=', '').replace('stderr=', '') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Rf = reg(fpkmvsotu1)\n",
    "Rpval001 = filt_(Rf.pvalue, 0.001, '<') #all the pvals < 0.00001\n",
    "Rrsq06 = filt_(Rf.rsquare, 0.6, '>') #all the rsquare > 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Rf = reg(fpkmvsotu1)\n",
    "Rpval = filt_(Rf.pvalue, 1.0, '<') #all the pvals \n",
    "Rrsq = filt_(Rf.rsquare, 0.0, '>') #all the rsquare "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "65"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prova = Rf.rsquare\n",
    "len(prova)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/software/Python/2.7.9/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    }
   ],
   "source": [
    "Xp = []\n",
    "Yp = []\n",
    "Otu = []\n",
    "Rna = []\n",
    "for i in range(12):\n",
    "    Otu.append(Rpval.split('\\n')[i].split(',')[0])\n",
    "    Rna.append(Rpval.split('\\n')[i].split(',')[1])\n",
    "    Xp.append(otu_d[otu_key[int(Rpval.split('\\n')[i].split(',')[0])]])\n",
    "    Yp.append(rna_d[rna_key[int(Rpval.split('\\n')[i].split(',')[1])]])\n",
    "    \n",
    "rng = zip(range(10), range(331,340))\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Top 9 regressions based on p value ', fontsize=14, fontweight='bold')\n",
    "\n",
    "# for i,j in rng:\n",
    "#     ax = fig.add_subplot(j)\n",
    "#     fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "#     #ax.set_title('axes title')\n",
    "#     ax.set_xlabel(Otu[i])\n",
    "#     ax.set_ylabel(Rna[i])\n",
    "#     plt.scatter(Xr[i],Yr[i])\n",
    "    \n",
    "# plt.show()\n",
    "\n",
    "for i,j in rng:\n",
    "    ax = fig.add_subplot(j)\n",
    "    fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "    #ax.set_title('axes title')\n",
    "    ax.set_xlabel(Otu[i])\n",
    "    ax.set_ylabel(Rna[i])\n",
    "    plt.scatter(Xp[i],Yp[i])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Xr = []\n",
    "Yr = []\n",
    "Otu = []\n",
    "Rna = []\n",
    "for i in range((len(Rrsq.split('\\n')) - 12),len(Rrsq.split('\\n'))):\n",
    "    Otu.append(Rrsq.split('\\n')[i].split(',')[0])\n",
    "    Rna.append(Rrsq.split('\\n')[i].split(',')[1])\n",
    "    Xr.append(otu_d[otu_key[int(Rrsq.split('\\n')[i].split(',')[0])]])\n",
    "    Yr.append(rna_d[rna_key[int(Rrsq.split('\\n')[i].split(',')[1])]])\n",
    "    \n",
    "fig = plt.figure()\n",
    "fig.suptitle('Top 9 regressions based on R2 ', fontsize=14, fontweight='bold')\n",
    "\n",
    "for i,j in rng:\n",
    "    ax = fig.add_subplot(j)\n",
    "    fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "    #ax.set_title('axes title')\n",
    "    ax.set_xlabel(Otu[i])\n",
    "    ax.set_ylabel(Rna[i])\n",
    "    plt.scatter(Xr[i],Yr[i])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pval = 'OTU,transcript,pvalue\\n'+','.join([ otu_key[int(i.split(',')[0])]+','+rna_key[int(i.split(',')[1])]+','+i.split(',')[2]+'\\n' for i in Rpval.split('\\n')]).replace('\\n,','\\n')\n",
    "rsq = 'OTU,transcript,rsquare\\n'+','.join([ otu_key[int(i.split(',')[0])]+','+rna_key[int(i.split(',')[1])]+','+i.split(',')[2]+'\\n' for i in Rrsq.split('\\n')]).replace('\\n,','\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/all_p.csv', 'w')\n",
    "z.write(pval)\n",
    "z.close()\n",
    "\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/all_r2.csv', 'w')\n",
    "z.write(rsq)\n",
    "z.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#The order of the individuals on each plot is:\n",
    "oder = [i[0] for i in uTs][1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**27-Apr-2016 - Goal: to plot correlations with individuals from various groups with different colors**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Groups = [i[0][:-2] for i in uTs][1:] # this establishes the groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Xr = []\n",
    "Yr = []\n",
    "Otu = []\n",
    "Rna = []\n",
    "for i in range((len(Rrsq.split('\\n')) - 12),len(Rrsq.split('\\n'))):\n",
    "    Otu.append(Rrsq.split('\\n')[i].split(',')[0])\n",
    "    Rna.append(Rrsq.split('\\n')[i].split(',')[1])\n",
    "    Xr.append(otu_d[otu_key[int(Rrsq.split('\\n')[i].split(',')[0])]])\n",
    "    Yr.append(rna_d[rna_key[int(Rrsq.split('\\n')[i].split(',')[1])]])\n",
    "    \n",
    "fig = plt.figure()\n",
    "fig.suptitle('Top 9 regressions based on R2 ', fontsize=14, fontweight='bold')\n",
    "\n",
    "for i,j in rng:\n",
    "    ax = fig.add_subplot(j)\n",
    "    fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "    #ax.set_title('axes title')\n",
    "    ax.set_xlabel(Otu[i])\n",
    "    ax.set_ylabel(Rna[i])\n",
    "    plt.scatter(Xr[i],Yr[i])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Below the group coloring works for a single plot**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfR = pd.DataFrame(dict(x=Xr[1], y=Yr[1], label=Groups))\n",
    "groupsR = df.groupby('label')\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.margins(0.05) # Optional, just adds 5% padding to the autoscaling\n",
    "for name, group in groupsR:\n",
    "    ax.plot(group.x, group.y, marker='o', linestyle='', ms=12, label=name)\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I am trying to achieve the same for the subplots with the 9 top correlations - it works, although the legend  \n",
    "does not look too nice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "fig.suptitle('Top 9 regressions based on R2 ', fontsize=14, fontweight='bold')\n",
    "\n",
    "for i,j in rng:\n",
    "    ax = fig.add_subplot(j)\n",
    "    fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "    ax.set_title('axes title')\n",
    "    ax.set_xlabel(Otu[i])\n",
    "    ax.set_ylabel(Rna[i])\n",
    "    dfR = pd.DataFrame(dict(x=Xr[i], y=Yr[i], label=Groups))\n",
    "    groupsR = dfR.groupby('label')\n",
    "    for name, group in groupsR:\n",
    "        ax.plot(group.x, group.y, marker='o', linestyle='', ms=6, label=name)\n",
    "ax.legend()    \n",
    "# plt.legend( loc = 'lower center', bbox_to_anchor = (0,-0.1,1,1),\n",
    "#             bbox_transform = plt.gcf().transFigure )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The decision is that the top plots will be displayed in b&w, while the single ones can be displayed separating group by group  \n",
    "\n",
    "Now I will work on the calculation of the **q values** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pval = open('/beegfs/group_dv/home/DValenzano/uBiome/data/all_p.csv', 'rU').read()\n",
    "rsq = open('/beegfs/group_dv/home/DValenzano/uBiome/data/all_r2.csv', 'rU').read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pvals = np.array([ float(i.split(',')[-1]) for i in pval.split('\\n')[1:-1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from statsmodels.sandbox.stats.multicomp import multipletests #from here: http://statsmodels.sourceforge.net/devel/generated/statsmodels.sandbox.stats.multicomp.multipletests.html#statsmodels.sandbox.stats.multicomp.multipletests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "prova = multipletests(pvals, alpha=0.5, method='fdr_tsbky', is_sorted=True, returnsorted=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import rpy2.robjects as robjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589793"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "robjects.r['pi'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below the density plot for the pvalue distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(pvals, bins = 50, histtype='stepfilled', color= 'grey')\n",
    "plt.title(\"p values density histogram\")\n",
    "plt.xlabel(\"Value\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r2vals = np.array([ float(i.split(',')[-1]) for i in rsq.split('\\n')[1:-1]])\n",
    "plt.hist(r2vals, bins = 50, histtype='stepfilled', color= 'grey')\n",
    "plt.title(\"p values density histogram\")\n",
    "plt.xlabel(\"Value\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I am now going to log2-transform all the fpkm values and then recalculate the correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import math\n",
    "#[ math.log(i)/math.log(2) for i in rna_d[rna_key[int(Rrsq.split('\\n')[-2].split(',')[1])]]]\n",
    "#rna_d[rna_key[int(Rrsq.split('\\n')[-2].split(',')[1])]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rnat_z2 = [ [ math.log(j+1)/math.log(2) for j in i ] for i in rnat_z ]\n",
    "\n",
    "rna_value2 = rnat_z2\n",
    "otu_value = otut_z\n",
    "\n",
    "rna_d2 = dict(zip(rna_key, rna_value2))\n",
    "otu_d2 = dict(zip(otu_key, otu_value))\n",
    "\n",
    "def loop_one(m1,otu): # this returns all the correlation coefficients between m1 and a given otu (otu)\n",
    "    ls = []\n",
    "    for i in m1:\n",
    "        ls.append(stats.linregress(i, otu))    #this function returns slope, intercept, r_value, p_value, std_err\n",
    "    return ','.join(map(str, ls))+'\\n'\n",
    "\n",
    "def loop_two(m1, m2):\n",
    "    ls = []\n",
    "    for j in m2:\n",
    "        ls.append(loop_one(m1, j))\n",
    "    return ','.join(ls).replace('\\n,','\\n')\n",
    "\n",
    "fpkmvsotu2 = loop_two(rnat_z2, otut_z)\n",
    "fpkmvsotu2_1 = fpkmvsotu2.replace('LinregressResult', '').replace('slope=','').replace('intercept=', '').replace('rvalue=', '').replace('pvalue=', '').replace('stderr=', '') \n",
    "\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/corr_log2.csv','w')\n",
    "z.write(fpkmvsotu2_1)\n",
    "z.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Rf2 = reg(fpkmvsotu2_1)\n",
    "Rpval2 = filt_(Rf2.pvalue, 1.0, '<') #all the pvals \n",
    "Rrsq2 = filt_(Rf2.rsquare, 0.0, '>') #all the rsquare "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Xp2 = []\n",
    "Yp2 = []\n",
    "Otu2 = []\n",
    "Rna2 = []\n",
    "for i in range(12):\n",
    "    Otu2.append(Rpval2.split('\\n')[i].split(',')[0])\n",
    "    Rna2.append(Rpval2.split('\\n')[i].split(',')[1])\n",
    "    Xp2.append(otu_d2[otu_key[int(Rpval2.split('\\n')[i].split(',')[0])]])\n",
    "    Yp2.append(rna_d2[rna_key[int(Rpval2.split('\\n')[i].split(',')[1])]])\n",
    "    \n",
    "rng = zip(range(10), range(331,340))\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Top 9 regressions based on p value ', fontsize=14, fontweight='bold')\n",
    "\n",
    "for i,j in rng:\n",
    "    ax = fig.add_subplot(j)\n",
    "    fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "    #ax.set_title('axes title')\n",
    "    ax.set_xlabel(Otu2[i])\n",
    "    ax.set_ylabel(Rna2[i])\n",
    "    plt.scatter(Xp2[i],Yp2[i])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Xp2 = []\n",
    "Yp2 = []\n",
    "Otu2 = []\n",
    "Rna2 = []\n",
    "for i in range((len(Rrsq2.split('\\n')) - 12),len(Rrsq2.split('\\n'))):\n",
    "    Otu2.append(Rrsq2.split('\\n')[i].split(',')[0])\n",
    "    Rna2.append(Rrsq2.split('\\n')[i].split(',')[1])\n",
    "    Xp2.append(otu_d2[otu_key[int(Rrsq2.split('\\n')[i].split(',')[0])]])\n",
    "    Yp2.append(rna_d2[rna_key[int(Rrsq2.split('\\n')[i].split(',')[1])]])\n",
    "    \n",
    "rng = zip(range(10), range(331,340))\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Top 9 regressions based on R square value ', fontsize=14, fontweight='bold')\n",
    "\n",
    "for i,j in rng:\n",
    "    ax = fig.add_subplot(j)\n",
    "    fig.subplots_adjust(top=0.85, hspace=0.8, wspace=0.8)\n",
    "    #ax.set_title('axes title')\n",
    "    ax.set_xlabel(Otu2[i])\n",
    "    ax.set_ylabel(Rna2[i])\n",
    "    plt.scatter(Xp2[i],Yp2[i])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pval2 = 'OTU,transcript,pvalue\\n'+','.join([ otu_key[int(i.split(',')[0])]+','+rna_key[int(i.split(',')[1])]+','+i.split(',')[2]+'\\n' for i in Rpval2.split('\\n')]).replace('\\n,','\\n')\n",
    "rsq2 = 'OTU,transcript,rsquare\\n'+','.join([ otu_key[int(i.split(',')[0])]+','+rna_key[int(i.split(',')[1])]+','+i.split(',')[2]+'\\n' for i in Rrsq2.split('\\n')]).replace('\\n,','\\n')\n",
    "\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/all_p2.csv', 'w')\n",
    "z.write(pval2)\n",
    "z.close()\n",
    "\n",
    "z = open('/beegfs/group_dv/home/DValenzano/uBiome/data/all_r2_2.csv', 'w')\n",
    "z.write(rsq2)\n",
    "z.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pvals2 = np.array([ float(i.split(',')[-1]) for i in pval2.split('\\n')[1:-1]])\n",
    "plt.hist(pvals2, bins = 50, histtype='stepfilled', color= 'grey')\n",
    "plt.title(\"p values density histogram\")\n",
    "plt.xlabel(\"Value\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r2vals2 = np.array([ float(i.split(',')[-1]) for i in rsq2.split('\\n')[1:-1]])\n",
    "plt.hist(r2vals, bins = 50, histtype='stepfilled', color= 'grey')\n",
    "plt.title(\"p values density histogram\")\n",
    "plt.xlabel(\"Value\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top_p2 = Rpval2.split('\\n')[:100]\n",
    "top_r2 = Rrsq2.split('\\n')[-100:][::-1] #this reverses the list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p2_k = [','.join(i.split(',')[:2]) for i in top_p2]\n",
    "p2_v = [ i.split(',')[2] for i in top_p2]\n",
    "p2_d = dict(zip(p2_k, p2_v))\n",
    "\n",
    "r2_k = [','.join(i.split(',')[:2]) for i in top_r2]\n",
    "r2_v = [ i.split(',')[2] for i in top_r2]\n",
    "r2_d = dict(zip(r2_k, r2_v))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/beegfs/group_dv/home/DValenzano/.Python/2.7/lib/python2.7/site-packages/ipykernel/__main__.py:1: DeprecationWarning: the sets module is deprecated\n",
      "  if __name__ == '__main__':\n"
     ]
    }
   ],
   "source": [
    "from sets import Set\n",
    "top = Set(r2_k) & Set(p2_k)"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}