class21_reveal_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 21:  joint entropy and the REVEAL algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the bladder cancer gene expression data to test out the REVEAL algorithm. First, we'll load the data and filter to include only genes for which the median log2 expression level is > 12 (as we did in class session 20).  That should give us 164 genes to work with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the Python modules that we will need for this exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import numpy\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data file `shared/bladder_cancer_genes_tcga.txt` into a `pandas.DataFrame`, convert it to a `numpy.ndarray` matrix, and print the matrix dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4473, 414)\n"
     ]
    }
   ],
   "source": [
    "gene_matrix_for_network_df = pandas.read_csv(\"shared/bladder_cancer_genes_tcga.txt\", sep=\"\\t\")\n",
    "gene_matrix_for_network = gene_matrix_for_network_df.as_matrix()\n",
    "print(gene_matrix_for_network.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter the matrix to include only rows for which the column-wise median is > 14; matrix should now be 13 x 414."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_keep = numpy.where(numpy.median(gene_matrix_for_network, axis=1) > 14)\n",
    "matrix_filt = gene_matrix_for_network[genes_keep, ][0]\n",
    "matrix_filt.shape\n",
    "N = matrix_filt.shape[0]\n",
    "M = matrix_filt.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binarize the gene expression matrix using the mean value as a breakpoint, turning it into a NxM matrix of booleans (`True`/`False`).  Call it `gene_matrix_binarized`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(13, 414)\n"
     ]
    }
   ],
   "source": [
    "gene_matrix_binarized = numpy.tile(numpy.mean(matrix_filt, axis=1),(M,1)).transpose() < matrix_filt\n",
    "print(gene_matrix_binarized.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your matrix by printing the first four columns of the first four rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True, False, False],\n",
       "       [False,  True, False, False],\n",
       "       [ True,  True,  True, False],\n",
       "       [False,  True, False, False]], dtype=bool)"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_matrix_binarized[0:4,0:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The core part of the REVEAL algorithm is a function that can compute the joint entropy of a collection of binary (TRUE/FALSE) vectors X1, X2, ..., Xn (where length(X1) = length(Xi) = M).\n",
    "Write a function `entropy_multiple_vecs` that takes as its input a nxM matrix (where n is the number of variables, i.e., genes, and M is the number of samples in which gene expression was measured). The function should use the log2 definition of the Shannon entropy. It should return the joint entropy H(X1, X2, ..., Xn) as a scalar numeric value. I have created a skeleton version of this function for you, in which you can fill in the code. I have also created some test code that you can use to test your function, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy_multiple_vecs(binary_vecs):\n",
    "    ## use shape to get the numbers of rows and columns as [n,M]\n",
    "    [n, M] = binary_vecs.shape\n",
    "    \n",
    "    # make a \"M x n\" dataframe from the transpose of the matrix binary_vecs\n",
    "    binary_df = pandas.DataFrame(binary_vecs.transpose())\n",
    "    \n",
    "    # use the groupby method to obtain a data frame of counts of unique occurrences of the 2^n possible logical states\n",
    "    binary_df_counts = binary_df.groupby(binary_df.columns.values.tolist()).size().values\n",
    "    \n",
    "    # divide the matrix of counts by M, to get a probability matrix\n",
    "    probvec = binary_df_counts/M\n",
    "    \n",
    "    # compute the shannon entropy using the formula\n",
    "    hvec = -probvec*numpy.log2(probvec)\n",
    "    return numpy.sum(hvec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This test case should produce the value 3.938:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.81175425925\n"
     ]
    }
   ],
   "source": [
    "print(entropy_multiple_vecs(gene_matrix_binarized[0:4,]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example implementation of the REVEAL algorithm:\n",
    "We'll go through stage 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ratio_thresh = 0.1\n",
    "genes_to_fit = list(range(0,N))\n",
    "stage = 0\n",
    "regulators = [None]*N\n",
    "entropies_for_stages = [None]*N\n",
    "max_stage = 4\n",
    "\n",
    "entropies_for_stages[0] = numpy.zeros(N)\n",
    "\n",
    "for i in range(0,N):\n",
    "    single_row_matrix = gene_matrix_binarized[i,:,None].transpose()\n",
    "    entropies_for_stages[0][i] = entropy_multiple_vecs(single_row_matrix)\n",
    "    \n",
    "genes_to_fit = set(range(0,N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1],\n",
       " [0],\n",
       " [8],\n",
       " None,\n",
       " None,\n",
       " None,\n",
       " [7],\n",
       " [6],\n",
       " [2],\n",
       " [7, 8, 10],\n",
       " [11],\n",
       " [10],\n",
       " [1, 6, 9]]"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for stage in range(1,max_stage + 1):\n",
    "    for gene in genes_to_fit.copy():\n",
    "        # we are trying to find regulators for gene \"gene\"\n",
    "        poss_regs = set(range(0,N)) - set([gene])\n",
    "        poss_regs_combs = [list(x) for x in itertools.combinations(poss_regs, stage)]\n",
    "        HGX = numpy.array([ entropy_multiple_vecs(gene_matrix_binarized[[gene] + poss_regs_comb,:]) for poss_regs_comb in poss_regs_combs ])\n",
    "        HX = numpy.array([ entropy_multiple_vecs(gene_matrix_binarized[poss_regs_comb,:]) for poss_regs_comb in poss_regs_combs ])\n",
    "        HG = entropies_for_stages[0][gene]\n",
    "        min_value = numpy.min(HGX - HX)\n",
    "        if HG - min_value >= ratio_thresh * HG:\n",
    "            regulators[gene]=poss_regs_combs[numpy.argmin(HGX - HX)]\n",
    "            genes_to_fit.remove(gene)\n",
    "regulators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}