class15_reveal_python3_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 15:  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 > 14.  That should give us 13 genes to work with (to keep running times reasonable for an in-class notebook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the Python modules that we will need for this exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4473, 414)\n"
     ]
    }
   ],
   "source": [
    "gene_matrix_for_network_df = \n",
    "gene_matrix_for_network = \n"
   ]
  },
  {
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_keep = \n",
    "matrix_filt = \n",
    "\n"
   ]
  },
  {
   "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`.  Use `numpy.tile` and `numpy.mean` and `transpose`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(13, 414)\n"
     ]
    }
   ],
   "source": [
    "gene_matrix_binarized = \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your matrix by printing the first four columns of the first four rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "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": 6,
   "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] = \n",
    "    \n",
    "    # make a \"M x n\" dataframe from the transpose of the matrix binary_vecs\n",
    "    binary_df = \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 = \n",
    "    \n",
    "    # divide the matrix of counts by M, to get a probability matrix\n",
    "    probvec = \n",
    "    \n",
    "    # compute the shannon entropy using the formula\n",
    "    hvec = \n",
    "    return "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This test case should produce the value 3.938:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.471945238516729\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": 8,
   "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": 10,
   "metadata": {},
   "outputs": [],
   "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",
    "results = [[index, value] for index, value in enumerate(regulators)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, [1]],\n",
       " [1, [0]],\n",
       " [2, [8]],\n",
       " [3, [1, 2, 4, 6]],\n",
       " [4, [0, 3, 5, 11]],\n",
       " [5, [1, 2, 6, 10]],\n",
       " [6, [7]],\n",
       " [7, [6]],\n",
       " [8, [2]],\n",
       " [9, [7, 8, 10]],\n",
       " [10, [11]],\n",
       " [11, [10]],\n",
       " [12, [1, 6, 9]]]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results"
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}