class19_correlnet_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CS446/546 - Class Session 19 - Correlation network\n",
    "\n",
    "In this class session we are going to analyze gene expression data from a human bladder cancer cohort, using python. We will load a data matrix of expression measurements of 4,473 genes in 414 different bladder cancer samples. These genes have been selected because they are differentially expressed between normal bladder and bladder cancer (thus more likely to have a function in bladder cancer specifically), but the columns in the data matrix are restricted to bladder cancer samples (not normal bladder) because we want to obtain a network representing variation across cancers. The measurements in the matrix have already been normalized to account for inter-sample heterogeneity and then log2 transformed. Our job is to compute Pearson correlation coefficients between all pairs of genes, obtain Fisher-transformed *z*-scores for all pairs of genes, test each pair of genes for significance of the *z* score, adjust for multiple hypothesis testing, filter to eliminate any pair for which *R* < 0.75 or *P*adj > 0.01, load the graph into an `igraph.Graph` object, and plot the degree distribution on log-log scale. We will then answer two questions: (1) does the network look to be scale-free? and (2) what is it's best-fit scaling exponent?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start by importing all of the modules that we will need for this notebook. Note the difference in language-design philosophy between R (which requires loading one package for this analysis) and python (where we have to load seven modules). Python keeps its core minimal, whereas R has a lot of statistical and plotting functions in the base language (or in packages that are loaded by default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import scipy.stats\n",
    "import matplotlib\n",
    "import pylab\n",
    "import numpy\n",
    "import statsmodels.sandbox.stats.multicomp\n",
    "import igraph\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `pandas.read_csv`, load the tab-deliminted text file of gene expression measurements (rows correspond to genes, columns correspond to bladder tumor samples), into a data frame `gene_matrix_for_network_df`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_matrix_for_network_df = pandas.read_csv(\"shared/bladder_cancer_genes_tcga.txt\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `pandas.DataFrame.as_matrix` method to make a matrix `gene_matrix_for_network`. Print out the dimensions of the matrix, by accessing its `shape` variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4473, 414)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_matrix_for_network = gene_matrix_for_network_df.as_matrix()\n",
    "gene_matrix_for_network.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `del` to delete the data frame, since we no longer need it (save memory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "del gene_matrix_for_network_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the online help for the `numpy.corrcoef` function, using `help(numpy.corrcoef)`.  When you pass a single argument `x` which is a 2D \"array\" (i.e., a matrix), by default does `corrcoef` compute coefficients for pairs of rows, or pairs of columns?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function corrcoef in module numpy.lib.function_base:\n",
      "\n",
      "corrcoef(x, y=None, rowvar=True, bias=<class 'numpy._globals._NoValue'>, ddof=<class 'numpy._globals._NoValue'>)\n",
      "    Return Pearson product-moment correlation coefficients.\n",
      "    \n",
      "    Please refer to the documentation for `cov` for more detail.  The\n",
      "    relationship between the correlation coefficient matrix, `R`, and the\n",
      "    covariance matrix, `C`, is\n",
      "    \n",
      "    .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} * C_{jj} } }\n",
      "    \n",
      "    The values of `R` are between -1 and 1, inclusive.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    x : array_like\n",
      "        A 1-D or 2-D array containing multiple variables and observations.\n",
      "        Each row of `x` represents a variable, and each column a single\n",
      "        observation of all those variables. Also see `rowvar` below.\n",
      "    y : array_like, optional\n",
      "        An additional set of variables and observations. `y` has the same\n",
      "        shape as `x`.\n",
      "    rowvar : bool, optional\n",
      "        If `rowvar` is True (default), then each row represents a\n",
      "        variable, with observations in the columns. Otherwise, the relationship\n",
      "        is transposed: each column represents a variable, while the rows\n",
      "        contain observations.\n",
      "    bias : _NoValue, optional\n",
      "        Has no effect, do not use.\n",
      "    \n",
      "        .. deprecated:: 1.10.0\n",
      "    ddof : _NoValue, optional\n",
      "        Has no effect, do not use.\n",
      "    \n",
      "        .. deprecated:: 1.10.0\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    R : ndarray\n",
      "        The correlation coefficient matrix of the variables.\n",
      "    \n",
      "    See Also\n",
      "    --------\n",
      "    cov : Covariance matrix\n",
      "    \n",
      "    Notes\n",
      "    -----\n",
      "    Due to floating point rounding the resulting array may not be Hermitian,\n",
      "    the diagonal elements may not be 1, and the elements may not satisfy the\n",
      "    inequality abs(a) <= 1. The real and imaginary parts are clipped to the\n",
      "    interval [-1,  1] in an attempt to improve on that situation but is not\n",
      "    much help in the complex case.\n",
      "    \n",
      "    This function accepts but discards arguments `bias` and `ddof`.  This is\n",
      "    for backwards compatibility with previous versions of this function.  These\n",
      "    arguments had no effect on the return values of the function and can be\n",
      "    safely ignored in this and previous versions of numpy.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(numpy.corrcoef)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the 4,473 x 4,473 matrix of gene-gene Pearson correlation coefficients, using `numpy.corrcoef` (this function treats each row as a variable, so you don't have to do any transposing of the matrix, unlike the situation in R)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_matrix_for_network_cor = numpy.corrcoef(gene_matrix_for_network)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the online help for `numpy.fill_diagonal`.  Does it return the modified matrix or modify the matrix argument in place?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function fill_diagonal in module numpy.lib.index_tricks:\n",
      "\n",
      "fill_diagonal(a, val, wrap=False)\n",
      "    Fill the main diagonal of the given array of any dimensionality.\n",
      "    \n",
      "    For an array `a` with ``a.ndim >= 2``, the diagonal is the list of\n",
      "    locations with indices ``a[i, ..., i]`` all identical. This function\n",
      "    modifies the input array in-place, it does not return a value.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    a : array, at least 2-D.\n",
      "      Array whose diagonal is to be filled, it gets modified in-place.\n",
      "    \n",
      "    val : scalar\n",
      "      Value to be written on the diagonal, its type must be compatible with\n",
      "      that of the array a.\n",
      "    \n",
      "    wrap : bool\n",
      "      For tall matrices in NumPy version up to 1.6.2, the\n",
      "      diagonal \"wrapped\" after N columns. You can have this behavior\n",
      "      with this option. This affects only tall matrices.\n",
      "    \n",
      "    See also\n",
      "    --------\n",
      "    diag_indices, diag_indices_from\n",
      "    \n",
      "    Notes\n",
      "    -----\n",
      "    .. versionadded:: 1.4.0\n",
      "    \n",
      "    This functionality can be obtained via `diag_indices`, but internally\n",
      "    this version uses a much faster implementation that never constructs the\n",
      "    indices and uses simple slicing.\n",
      "    \n",
      "    Examples\n",
      "    --------\n",
      "    >>> a = np.zeros((3, 3), int)\n",
      "    >>> np.fill_diagonal(a, 5)\n",
      "    >>> a\n",
      "    array([[5, 0, 0],\n",
      "           [0, 5, 0],\n",
      "           [0, 0, 5]])\n",
      "    \n",
      "    The same function can operate on a 4-D array:\n",
      "    \n",
      "    >>> a = np.zeros((3, 3, 3, 3), int)\n",
      "    >>> np.fill_diagonal(a, 4)\n",
      "    \n",
      "    We only show a few blocks for clarity:\n",
      "    \n",
      "    >>> a[0, 0]\n",
      "    array([[4, 0, 0],\n",
      "           [0, 0, 0],\n",
      "           [0, 0, 0]])\n",
      "    >>> a[1, 1]\n",
      "    array([[0, 0, 0],\n",
      "           [0, 4, 0],\n",
      "           [0, 0, 0]])\n",
      "    >>> a[2, 2]\n",
      "    array([[0, 0, 0],\n",
      "           [0, 0, 0],\n",
      "           [0, 0, 4]])\n",
      "    \n",
      "    The wrap option affects only tall matrices:\n",
      "    \n",
      "    >>> # tall matrices no wrap\n",
      "    >>> a = np.zeros((5, 3),int)\n",
      "    >>> fill_diagonal(a, 4)\n",
      "    >>> a\n",
      "    array([[4, 0, 0],\n",
      "           [0, 4, 0],\n",
      "           [0, 0, 4],\n",
      "           [0, 0, 0],\n",
      "           [0, 0, 0]])\n",
      "    \n",
      "    >>> # tall matrices wrap\n",
      "    >>> a = np.zeros((5, 3),int)\n",
      "    >>> fill_diagonal(a, 4, wrap=True)\n",
      "    >>> a\n",
      "    array([[4, 0, 0],\n",
      "           [0, 4, 0],\n",
      "           [0, 0, 4],\n",
      "           [0, 0, 0],\n",
      "           [4, 0, 0]])\n",
      "    \n",
      "    >>> # wide matrices\n",
      "    >>> a = np.zeros((3, 5),int)\n",
      "    >>> fill_diagonal(a, 4, wrap=True)\n",
      "    >>> a\n",
      "    array([[4, 0, 0, 0, 0],\n",
      "           [0, 4, 0, 0, 0],\n",
      "           [0, 0, 4, 0, 0]])\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(numpy.fill_diagonal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the diagonal elements of the matrix to zero, using `numpy.fill_diagonal`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "numpy.fill_diagonal(gene_matrix_for_network_cor, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the online help for `numpy.multiply`.  Does it do element-wise multiplication or matrix multiplication?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on ufunc object:\n",
      "\n",
      "multiply = class ufunc(builtins.object)\n",
      " |  Functions that operate element by element on whole arrays.\n",
      " |  \n",
      " |  To see the documentation for a specific ufunc, use `info`.  For\n",
      " |  example, ``np.info(np.sin)``.  Because ufuncs are written in C\n",
      " |  (for speed) and linked into Python with NumPy's ufunc facility,\n",
      " |  Python's help() function finds this page whenever help() is called\n",
      " |  on a ufunc.\n",
      " |  \n",
      " |  A detailed explanation of ufuncs can be found in the docs for :ref:`ufuncs`.\n",
      " |  \n",
      " |  Calling ufuncs:\n",
      " |  ===============\n",
      " |  \n",
      " |  op(*x[, out], where=True, **kwargs)\n",
      " |  Apply `op` to the arguments `*x` elementwise, broadcasting the arguments.\n",
      " |  \n",
      " |  The broadcasting rules are:\n",
      " |  \n",
      " |  * Dimensions of length 1 may be prepended to either array.\n",
      " |  * Arrays may be repeated along dimensions of length 1.\n",
      " |  \n",
      " |  Parameters\n",
      " |  ----------\n",
      " |  *x : array_like\n",
      " |      Input arrays.\n",
      " |  out : ndarray, None, or tuple of ndarray and None, optional\n",
      " |      Alternate array object(s) in which to put the result; if provided, it\n",
      " |      must have a shape that the inputs broadcast to. A tuple of arrays\n",
      " |      (possible only as a keyword argument) must have length equal to the\n",
      " |      number of outputs; use `None` for outputs to be allocated by the ufunc.\n",
      " |  where : array_like, optional\n",
      " |      Values of True indicate to calculate the ufunc at that position, values\n",
      " |      of False indicate to leave the value in the output alone.\n",
      " |  **kwargs\n",
      " |      For other keyword-only arguments, see the :ref:`ufunc docs <ufuncs.kwargs>`.\n",
      " |  \n",
      " |  Returns\n",
      " |  -------\n",
      " |  r : ndarray or tuple of ndarray\n",
      " |      `r` will have the shape that the arrays in `x` broadcast to; if `out` is\n",
      " |      provided, `r` will be equal to `out`. If the function has more than one\n",
      " |      output, then the result will be a tuple of arrays.\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __call__(self, /, *args, **kwargs)\n",
      " |      Call self as a function.\n",
      " |  \n",
      " |  __repr__(self, /)\n",
      " |      Return repr(self).\n",
      " |  \n",
      " |  __str__(self, /)\n",
      " |      Return str(self).\n",
      " |  \n",
      " |  accumulate(...)\n",
      " |      accumulate(array, axis=0, dtype=None, out=None, keepdims=None)\n",
      " |      \n",
      " |      Accumulate the result of applying the operator to all elements.\n",
      " |      \n",
      " |      For a one-dimensional array, accumulate produces results equivalent to::\n",
      " |      \n",
      " |        r = np.empty(len(A))\n",
      " |        t = op.identity        # op = the ufunc being applied to A's  elements\n",
      " |        for i in range(len(A)):\n",
      " |            t = op(t, A[i])\n",
      " |            r[i] = t\n",
      " |        return r\n",
      " |      \n",
      " |      For example, add.accumulate() is equivalent to np.cumsum().\n",
      " |      \n",
      " |      For a multi-dimensional array, accumulate is applied along only one\n",
      " |      axis (axis zero by default; see Examples below) so repeated use is\n",
      " |      necessary if one wants to accumulate over multiple axes.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      array : array_like\n",
      " |          The array to act on.\n",
      " |      axis : int, optional\n",
      " |          The axis along which to apply the accumulation; default is zero.\n",
      " |      dtype : data-type code, optional\n",
      " |          The data-type used to represent the intermediate results. Defaults\n",
      " |          to the data-type of the output array if such is provided, or the\n",
      " |          the data-type of the input array if no output array is provided.\n",
      " |      out : ndarray, None, or tuple of ndarray and None, optional\n",
      " |          A location into which the result is stored. If not provided or `None`,\n",
      " |          a freshly-allocated array is returned. For consistency with\n",
      " |          :ref:`ufunc.__call__`, if given as a keyword, this may be wrapped in a\n",
      " |          1-element tuple.\n",
      " |      \n",
      " |          .. versionchanged:: 1.13.0\n",
      " |             Tuples are allowed for keyword argument.\n",
      " |      keepdims : bool\n",
      " |          Has no effect. Deprecated, and will be removed in future.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      r : ndarray\n",
      " |          The accumulated values. If `out` was supplied, `r` is a reference to\n",
      " |          `out`.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      1-D array examples:\n",
      " |      \n",
      " |      >>> np.add.accumulate([2, 3, 5])\n",
      " |      array([ 2,  5, 10])\n",
      " |      >>> np.multiply.accumulate([2, 3, 5])\n",
      " |      array([ 2,  6, 30])\n",
      " |      \n",
      " |      2-D array examples:\n",
      " |      \n",
      " |      >>> I = np.eye(2)\n",
      " |      >>> I\n",
      " |      array([[ 1.,  0.],\n",
      " |             [ 0.,  1.]])\n",
      " |      \n",
      " |      Accumulate along axis 0 (rows), down columns:\n",
      " |      \n",
      " |      >>> np.add.accumulate(I, 0)\n",
      " |      array([[ 1.,  0.],\n",
      " |             [ 1.,  1.]])\n",
      " |      >>> np.add.accumulate(I) # no axis specified = axis zero\n",
      " |      array([[ 1.,  0.],\n",
      " |             [ 1.,  1.]])\n",
      " |      \n",
      " |      Accumulate along axis 1 (columns), through rows:\n",
      " |      \n",
      " |      >>> np.add.accumulate(I, 1)\n",
      " |      array([[ 1.,  1.],\n",
      " |             [ 0.,  1.]])\n",
      " |  \n",
      " |  at(...)\n",
      " |      at(a, indices, b=None)\n",
      " |      \n",
      " |      Performs unbuffered in place operation on operand 'a' for elements\n",
      " |      specified by 'indices'. For addition ufunc, this method is equivalent to\n",
      " |      `a[indices] += b`, except that results are accumulated for elements that\n",
      " |      are indexed more than once. For example, `a[[0,0]] += 1` will only\n",
      " |      increment the first element once because of buffering, whereas\n",
      " |      `add.at(a, [0,0], 1)` will increment the first element twice.\n",
      " |      \n",
      " |      .. versionadded:: 1.8.0\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      a : array_like\n",
      " |          The array to perform in place operation on.\n",
      " |      indices : array_like or tuple\n",
      " |          Array like index object or slice object for indexing into first\n",
      " |          operand. If first operand has multiple dimensions, indices can be a\n",
      " |          tuple of array like index objects or slice objects.\n",
      " |      b : array_like\n",
      " |          Second operand for ufuncs requiring two operands. Operand must be\n",
      " |          broadcastable over first operand after indexing or slicing.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      Set items 0 and 1 to their negative values:\n",
      " |      \n",
      " |      >>> a = np.array([1, 2, 3, 4])\n",
      " |      >>> np.negative.at(a, [0, 1])\n",
      " |      >>> print(a)\n",
      " |      array([-1, -2, 3, 4])\n",
      " |      \n",
      " |      ::\n",
      " |      \n",
      " |      Increment items 0 and 1, and increment item 2 twice:\n",
      " |      \n",
      " |      >>> a = np.array([1, 2, 3, 4])\n",
      " |      >>> np.add.at(a, [0, 1, 2, 2], 1)\n",
      " |      >>> print(a)\n",
      " |      array([2, 3, 5, 4])\n",
      " |      \n",
      " |      ::\n",
      " |      \n",
      " |      Add items 0 and 1 in first array to second array,\n",
      " |      and store results in first array:\n",
      " |      \n",
      " |      >>> a = np.array([1, 2, 3, 4])\n",
      " |      >>> b = np.array([1, 2])\n",
      " |      >>> np.add.at(a, [0, 1], b)\n",
      " |      >>> print(a)\n",
      " |      array([2, 4, 3, 4])\n",
      " |  \n",
      " |  outer(...)\n",
      " |      outer(A, B, **kwargs)\n",
      " |      \n",
      " |      Apply the ufunc `op` to all pairs (a, b) with a in `A` and b in `B`.\n",
      " |      \n",
      " |      Let ``M = A.ndim``, ``N = B.ndim``. Then the result, `C`, of\n",
      " |      ``op.outer(A, B)`` is an array of dimension M + N such that:\n",
      " |      \n",
      " |      .. math:: C[i_0, ..., i_{M-1}, j_0, ..., j_{N-1}] =\n",
      " |         op(A[i_0, ..., i_{M-1}], B[j_0, ..., j_{N-1}])\n",
      " |      \n",
      " |      For `A` and `B` one-dimensional, this is equivalent to::\n",
      " |      \n",
      " |        r = empty(len(A),len(B))\n",
      " |        for i in range(len(A)):\n",
      " |            for j in range(len(B)):\n",
      " |                r[i,j] = op(A[i], B[j]) # op = ufunc in question\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      A : array_like\n",
      " |          First array\n",
      " |      B : array_like\n",
      " |          Second array\n",
      " |      kwargs : any\n",
      " |          Arguments to pass on to the ufunc. Typically `dtype` or `out`.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      r : ndarray\n",
      " |          Output array\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      numpy.outer\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.multiply.outer([1, 2, 3], [4, 5, 6])\n",
      " |      array([[ 4,  5,  6],\n",
      " |             [ 8, 10, 12],\n",
      " |             [12, 15, 18]])\n",
      " |      \n",
      " |      A multi-dimensional example:\n",
      " |      \n",
      " |      >>> A = np.array([[1, 2, 3], [4, 5, 6]])\n",
      " |      >>> A.shape\n",
      " |      (2, 3)\n",
      " |      >>> B = np.array([[1, 2, 3, 4]])\n",
      " |      >>> B.shape\n",
      " |      (1, 4)\n",
      " |      >>> C = np.multiply.outer(A, B)\n",
      " |      >>> C.shape; C\n",
      " |      (2, 3, 1, 4)\n",
      " |      array([[[[ 1,  2,  3,  4]],\n",
      " |              [[ 2,  4,  6,  8]],\n",
      " |              [[ 3,  6,  9, 12]]],\n",
      " |             [[[ 4,  8, 12, 16]],\n",
      " |              [[ 5, 10, 15, 20]],\n",
      " |              [[ 6, 12, 18, 24]]]])\n",
      " |  \n",
      " |  reduce(...)\n",
      " |      reduce(a, axis=0, dtype=None, out=None, keepdims=False)\n",
      " |      \n",
      " |      Reduces `a`'s dimension by one, by applying ufunc along one axis.\n",
      " |      \n",
      " |      Let :math:`a.shape = (N_0, ..., N_i, ..., N_{M-1})`.  Then\n",
      " |      :math:`ufunc.reduce(a, axis=i)[k_0, ..,k_{i-1}, k_{i+1}, .., k_{M-1}]` =\n",
      " |      the result of iterating `j` over :math:`range(N_i)`, cumulatively applying\n",
      " |      ufunc to each :math:`a[k_0, ..,k_{i-1}, j, k_{i+1}, .., k_{M-1}]`.\n",
      " |      For a one-dimensional array, reduce produces results equivalent to:\n",
      " |      ::\n",
      " |      \n",
      " |       r = op.identity # op = ufunc\n",
      " |       for i in range(len(A)):\n",
      " |         r = op(r, A[i])\n",
      " |       return r\n",
      " |      \n",
      " |      For example, add.reduce() is equivalent to sum().\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      a : array_like\n",
      " |          The array to act on.\n",
      " |      axis : None or int or tuple of ints, optional\n",
      " |          Axis or axes along which a reduction is performed.\n",
      " |          The default (`axis` = 0) is perform a reduction over the first\n",
      " |          dimension of the input array. `axis` may be negative, in\n",
      " |          which case it counts from the last to the first axis.\n",
      " |      \n",
      " |          .. versionadded:: 1.7.0\n",
      " |      \n",
      " |          If this is `None`, a reduction is performed over all the axes.\n",
      " |          If this is a tuple of ints, a reduction is performed on multiple\n",
      " |          axes, instead of a single axis or all the axes as before.\n",
      " |      \n",
      " |          For operations which are either not commutative or not associative,\n",
      " |          doing a reduction over multiple axes is not well-defined. The\n",
      " |          ufuncs do not currently raise an exception in this case, but will\n",
      " |          likely do so in the future.\n",
      " |      dtype : data-type code, optional\n",
      " |          The type used to represent the intermediate results. Defaults\n",
      " |          to the data-type of the output array if this is provided, or\n",
      " |          the data-type of the input array if no output array is provided.\n",
      " |      out : ndarray, None, or tuple of ndarray and None, optional\n",
      " |          A location into which the result is stored. If not provided or `None`,\n",
      " |          a freshly-allocated array is returned. For consistency with\n",
      " |          :ref:`ufunc.__call__`, if given as a keyword, this may be wrapped in a\n",
      " |          1-element tuple.\n",
      " |      \n",
      " |          .. versionchanged:: 1.13.0\n",
      " |             Tuples are allowed for keyword argument.\n",
      " |      keepdims : bool, optional\n",
      " |          If this is set to True, the axes which are reduced are left\n",
      " |          in the result as dimensions with size one. With this option,\n",
      " |          the result will broadcast correctly against the original `arr`.\n",
      " |      \n",
      " |          .. versionadded:: 1.7.0\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      r : ndarray\n",
      " |          The reduced array. If `out` was supplied, `r` is a reference to it.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.multiply.reduce([2,3,5])\n",
      " |      30\n",
      " |      \n",
      " |      A multi-dimensional array example:\n",
      " |      \n",
      " |      >>> X = np.arange(8).reshape((2,2,2))\n",
      " |      >>> X\n",
      " |      array([[[0, 1],\n",
      " |              [2, 3]],\n",
      " |             [[4, 5],\n",
      " |              [6, 7]]])\n",
      " |      >>> np.add.reduce(X, 0)\n",
      " |      array([[ 4,  6],\n",
      " |             [ 8, 10]])\n",
      " |      >>> np.add.reduce(X) # confirm: default axis value is 0\n",
      " |      array([[ 4,  6],\n",
      " |             [ 8, 10]])\n",
      " |      >>> np.add.reduce(X, 1)\n",
      " |      array([[ 2,  4],\n",
      " |             [10, 12]])\n",
      " |      >>> np.add.reduce(X, 2)\n",
      " |      array([[ 1,  5],\n",
      " |             [ 9, 13]])\n",
      " |  \n",
      " |  reduceat(...)\n",
      " |      reduceat(a, indices, axis=0, dtype=None, out=None)\n",
      " |      \n",
      " |      Performs a (local) reduce with specified slices over a single axis.\n",
      " |      \n",
      " |      For i in ``range(len(indices))``, `reduceat` computes\n",
      " |      ``ufunc.reduce(a[indices[i]:indices[i+1]])``, which becomes the i-th\n",
      " |      generalized \"row\" parallel to `axis` in the final result (i.e., in a\n",
      " |      2-D array, for example, if `axis = 0`, it becomes the i-th row, but if\n",
      " |      `axis = 1`, it becomes the i-th column).  There are three exceptions to this:\n",
      " |      \n",
      " |      * when ``i = len(indices) - 1`` (so for the last index),\n",
      " |        ``indices[i+1] = a.shape[axis]``.\n",
      " |      * if ``indices[i] >= indices[i + 1]``, the i-th generalized \"row\" is\n",
      " |        simply ``a[indices[i]]``.\n",
      " |      * if ``indices[i] >= len(a)`` or ``indices[i] < 0``, an error is raised.\n",
      " |      \n",
      " |      The shape of the output depends on the size of `indices`, and may be\n",
      " |      larger than `a` (this happens if ``len(indices) > a.shape[axis]``).\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      a : array_like\n",
      " |          The array to act on.\n",
      " |      indices : array_like\n",
      " |          Paired indices, comma separated (not colon), specifying slices to\n",
      " |          reduce.\n",
      " |      axis : int, optional\n",
      " |          The axis along which to apply the reduceat.\n",
      " |      dtype : data-type code, optional\n",
      " |          The type used to represent the intermediate results. Defaults\n",
      " |          to the data type of the output array if this is provided, or\n",
      " |          the data type of the input array if no output array is provided.\n",
      " |      out : ndarray, None, or tuple of ndarray and None, optional\n",
      " |          A location into which the result is stored. If not provided or `None`,\n",
      " |          a freshly-allocated array is returned. For consistency with\n",
      " |          :ref:`ufunc.__call__`, if given as a keyword, this may be wrapped in a\n",
      " |          1-element tuple.\n",
      " |      \n",
      " |          .. versionchanged:: 1.13.0\n",
      " |             Tuples are allowed for keyword argument.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      r : ndarray\n",
      " |          The reduced values. If `out` was supplied, `r` is a reference to\n",
      " |          `out`.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      A descriptive example:\n",
      " |      \n",
      " |      If `a` is 1-D, the function `ufunc.accumulate(a)` is the same as\n",
      " |      ``ufunc.reduceat(a, indices)[::2]`` where `indices` is\n",
      " |      ``range(len(array) - 1)`` with a zero placed\n",
      " |      in every other element:\n",
      " |      ``indices = zeros(2 * len(a) - 1)``, ``indices[1::2] = range(1, len(a))``.\n",
      " |      \n",
      " |      Don't be fooled by this attribute's name: `reduceat(a)` is not\n",
      " |      necessarily smaller than `a`.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      To take the running sum of four successive values:\n",
      " |      \n",
      " |      >>> np.add.reduceat(np.arange(8),[0,4, 1,5, 2,6, 3,7])[::2]\n",
      " |      array([ 6, 10, 14, 18])\n",
      " |      \n",
      " |      A 2-D example:\n",
      " |      \n",
      " |      >>> x = np.linspace(0, 15, 16).reshape(4,4)\n",
      " |      >>> x\n",
      " |      array([[  0.,   1.,   2.,   3.],\n",
      " |             [  4.,   5.,   6.,   7.],\n",
      " |             [  8.,   9.,  10.,  11.],\n",
      " |             [ 12.,  13.,  14.,  15.]])\n",
      " |      \n",
      " |      ::\n",
      " |      \n",
      " |       # reduce such that the result has the following five rows:\n",
      " |       # [row1 + row2 + row3]\n",
      " |       # [row4]\n",
      " |       # [row2]\n",
      " |       # [row3]\n",
      " |       # [row1 + row2 + row3 + row4]\n",
      " |      \n",
      " |      >>> np.add.reduceat(x, [0, 3, 1, 2, 0])\n",
      " |      array([[ 12.,  15.,  18.,  21.],\n",
      " |             [ 12.,  13.,  14.,  15.],\n",
      " |             [  4.,   5.,   6.,   7.],\n",
      " |             [  8.,   9.,  10.,  11.],\n",
      " |             [ 24.,  28.,  32.,  36.]])\n",
      " |      \n",
      " |      ::\n",
      " |      \n",
      " |       # reduce such that result has the following two columns:\n",
      " |       # [col1 * col2 * col3, col4]\n",
      " |      \n",
      " |      >>> np.multiply.reduceat(x, [0, 3], 1)\n",
      " |      array([[    0.,     3.],\n",
      " |             [  120.,     7.],\n",
      " |             [  720.,    11.],\n",
      " |             [ 2184.,    15.]])\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  identity\n",
      " |      The identity value.\n",
      " |      \n",
      " |      Data attribute containing the identity element for the ufunc, if it has one.\n",
      " |      If it does not, the attribute value is None.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.add.identity\n",
      " |      0\n",
      " |      >>> np.multiply.identity\n",
      " |      1\n",
      " |      >>> np.power.identity\n",
      " |      1\n",
      " |      >>> print(np.exp.identity)\n",
      " |      None\n",
      " |  \n",
      " |  nargs\n",
      " |      The number of arguments.\n",
      " |      \n",
      " |      Data attribute containing the number of arguments the ufunc takes, including\n",
      " |      optional ones.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      Typically this value will be one more than what you might expect because all\n",
      " |      ufuncs take  the optional \"out\" argument.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.add.nargs\n",
      " |      3\n",
      " |      >>> np.multiply.nargs\n",
      " |      3\n",
      " |      >>> np.power.nargs\n",
      " |      3\n",
      " |      >>> np.exp.nargs\n",
      " |      2\n",
      " |  \n",
      " |  nin\n",
      " |      The number of inputs.\n",
      " |      \n",
      " |      Data attribute containing the number of arguments the ufunc treats as input.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.add.nin\n",
      " |      2\n",
      " |      >>> np.multiply.nin\n",
      " |      2\n",
      " |      >>> np.power.nin\n",
      " |      2\n",
      " |      >>> np.exp.nin\n",
      " |      1\n",
      " |  \n",
      " |  nout\n",
      " |      The number of outputs.\n",
      " |      \n",
      " |      Data attribute containing the number of arguments the ufunc treats as output.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      Since all ufuncs can take output arguments, this will always be (at least) 1.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.add.nout\n",
      " |      1\n",
      " |      >>> np.multiply.nout\n",
      " |      1\n",
      " |      >>> np.power.nout\n",
      " |      1\n",
      " |      >>> np.exp.nout\n",
      " |      1\n",
      " |  \n",
      " |  ntypes\n",
      " |      The number of types.\n",
      " |      \n",
      " |      The number of numerical NumPy types - of which there are 18 total - on which\n",
      " |      the ufunc can operate.\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      numpy.ufunc.types\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.add.ntypes\n",
      " |      18\n",
      " |      >>> np.multiply.ntypes\n",
      " |      18\n",
      " |      >>> np.power.ntypes\n",
      " |      17\n",
      " |      >>> np.exp.ntypes\n",
      " |      7\n",
      " |      >>> np.remainder.ntypes\n",
      " |      14\n",
      " |  \n",
      " |  signature\n",
      " |  \n",
      " |  types\n",
      " |      Returns a list with types grouped input->output.\n",
      " |      \n",
      " |      Data attribute listing the data-type \"Domain-Range\" groupings the ufunc can\n",
      " |      deliver. The data-types are given using the character codes.\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      numpy.ufunc.ntypes\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      >>> np.add.types\n",
      " |      ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',\n",
      " |      'LL->L', 'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D',\n",
      " |      'GG->G', 'OO->O']\n",
      " |      \n",
      " |      >>> np.multiply.types\n",
      " |      ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',\n",
      " |      'LL->L', 'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D',\n",
      " |      'GG->G', 'OO->O']\n",
      " |      \n",
      " |      >>> np.power.types\n",
      " |      ['bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L',\n",
      " |      'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D', 'GG->G',\n",
      " |      'OO->O']\n",
      " |      \n",
      " |      >>> np.exp.types\n",
      " |      ['f->f', 'd->d', 'g->g', 'F->F', 'D->D', 'G->G', 'O->O']\n",
      " |      \n",
      " |      >>> np.remainder.types\n",
      " |      ['bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L',\n",
      " |      'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'OO->O']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(numpy.multiply)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the online help for `numpy.tri`.  Does it modify a matrix argument in-place or return a matrix?  What is in the matrix that it returns?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function tri in module numpy.lib.twodim_base:\n",
      "\n",
      "tri(N, M=None, k=0, dtype=<class 'float'>)\n",
      "    An array with ones at and below the given diagonal and zeros elsewhere.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    N : int\n",
      "        Number of rows in the array.\n",
      "    M : int, optional\n",
      "        Number of columns in the array.\n",
      "        By default, `M` is taken equal to `N`.\n",
      "    k : int, optional\n",
      "        The sub-diagonal at and below which the array is filled.\n",
      "        `k` = 0 is the main diagonal, while `k` < 0 is below it,\n",
      "        and `k` > 0 is above.  The default is 0.\n",
      "    dtype : dtype, optional\n",
      "        Data type of the returned array.  The default is float.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    tri : ndarray of shape (N, M)\n",
      "        Array with its lower triangle filled with ones and zero elsewhere;\n",
      "        in other words ``T[i,j] == 1`` for ``i <= j + k``, 0 otherwise.\n",
      "    \n",
      "    Examples\n",
      "    --------\n",
      "    >>> np.tri(3, 5, 2, dtype=int)\n",
      "    array([[1, 1, 1, 0, 0],\n",
      "           [1, 1, 1, 1, 0],\n",
      "           [1, 1, 1, 1, 1]])\n",
      "    \n",
      "    >>> np.tri(3, 5, -1)\n",
      "    array([[ 0.,  0.,  0.,  0.,  0.],\n",
      "           [ 1.,  0.,  0.,  0.,  0.],\n",
      "           [ 1.,  1.,  0.,  0.,  0.]])\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(numpy.tri)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the upper-triangle of the matrix to zero, using `numpy.multiply` and `numpy.tri`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_matrix_for_network_cor = numpy.multiply(gene_matrix_for_network_cor, numpy.tri(*gene_matrix_for_network_cor.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `numpy.where`, get a tuple of two numpy.arrays containing the row/col indices of the entries of the matrix for which *R* >= 0.75. Use array indexing to obtain the *R* values for these matrix entries, as a numpy array `cor_coeff_values_above_thresh`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "inds_correl_above_thresh = numpy.where(gene_matrix_for_network_cor >= 0.75)\n",
    "cor_coeff_values_above_thresh = gene_matrix_for_network_cor[inds_correl_above_thresh]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refer to Eq. (13.5) in the assigned readding for today's class (p9 of the PDF).  Obtain a numpy array of the correlation coefficients that exceeded 0.75, and Fisher-transform the correlation coefficient values to get a vector `z_scores` of *z* scores. Each of these *z* scores will correspond to an **edge** in the network, unless the absolute *z* score is too small such that we can't exclude the null hypothesis that the corresponding two genes' expression values are indepdenent (we will perform that check in the next step)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_scores = 0.5*numpy.log((1 + cor_coeff_values_above_thresh)/\n",
    "                         (1 - cor_coeff_values_above_thresh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete the correlation matrix object in order to save memory (we won't need it from here on out)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "del gene_matrix_for_network_cor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that under the null hypothesis that two genes are independent, then sqrt(M-3)z for the pair of genes is an independent sample from the normal distribution with zero mean and unit variance, where M is the number of samples used to compute the Pearson correlation coefficient (i.e., M = 414). For each entry in `z_scores` compute a *P* value as the area under two tails of the normal distribution N(x),  where the two tails are x < -sqrt(M-3)z and x > sqrt(M-3)z. (You'll know you are doing it right if z=0 means you get a P value of 1).  You will want to use the functions `numpy.abs` and `scipy.stats.norm.cdf`, as well as the `math.sqrt` function (in order to compute the square root)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = gene_matrix_for_network.shape[1]\n",
    "P_values = 2*scipy.stats.norm.cdf(-numpy.abs(z_scores)*math.sqrt(M-3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adjust the P values for multiple hypothesis testing, using the `statsmodels.sandbox.stats.multicomp.multipletests` function wth `method=\"fdr_bh\"`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_values_adj = statsmodels.sandbox.stats.multicomp.multipletests(P_values, method=\"fdr_bh\")[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verify that we don't need to drop any entries due to the adjusted P value not being small enough (use `numpy.where` and `len`); this should produce zero since we have M=414 samples per gene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(numpy.where(P_values_adj >= 0.01)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the online help for the function `zip`.  What does it do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class zip in module builtins:\n",
      "\n",
      "class zip(object)\n",
      " |  zip(iter1 [,iter2 [...]]) --> zip object\n",
      " |  \n",
      " |  Return a zip object whose .__next__() method returns a tuple where\n",
      " |  the i-th element comes from the i-th iterable argument.  The .__next__()\n",
      " |  method continues until the shortest iterable in the argument sequence\n",
      " |  is exhausted and then it raises StopIteration.\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __getattribute__(self, name, /)\n",
      " |      Return getattr(self, name).\n",
      " |  \n",
      " |  __iter__(self, /)\n",
      " |      Implement iter(self).\n",
      " |  \n",
      " |  __new__(*args, **kwargs) from builtins.type\n",
      " |      Create and return a new object.  See help(type) for accurate signature.\n",
      " |  \n",
      " |  __next__(self, /)\n",
      " |      Implement next(self).\n",
      " |  \n",
      " |  __reduce__(...)\n",
      " |      Return state information for pickling.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(zip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to pass our tuple of numpy arrays containing row and column indices to `Graph.TupleList`; however, `Graph.TupleList` accepts a tuple list, not a tuple of numpy arrays. So we need to make a tuple list, using `zip`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_col_inds_tuple_list = zip(inds_correl_above_thresh[0], inds_correl_above_thresh[1])\n",
    "\n",
    "## [note this can be done more elegantly using the unary \"*\" operator:\n",
    "##    row_col_inds_tuple_list = zip(*inds_correl_above_thresh)  \n",
    "##  see how we only need to type the variable name once, if we use the unary \"*\" ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make an undirected graph from the row/column indices of the (upper-triangle) gene pairs whose  correlations were above our threshold, using `igraph.Graph.TupleList`. Print a summary of the network, as a sanity check, using the `igraph.Graph.summary` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH UN-- 1394 9916 -- \\n+ attr: name (v)'"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "final_network = igraph.Graph.TupleList(row_col_inds_tuple_list)\n",
    "final_network.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the degree distribution on log-log scale; does it appear to be scale-free?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEQCAYAAACqduMIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGG1JREFUeJzt3X+MHPV5x/HPc8tC11ZzB6n/qA9DSB0ddWJRkysgOf1B\nftRGiYNjSIkTIjW1sIhE1EqRVbuNkqhNhSMrUkGiSd0E0SipgRDrioOjayQToUaQcs6FOoQ4WFRg\nnyvhBM5t8Ck+n5/+sbvnvfV3d2d2d3Zmdt8vCYmdnZn9LozuszPP94e5uwAAqDeUdgMAANlEQAAA\ngggIAEAQAQEACCIgAABBBAQAIIiAAAAEERAAgKDMBISZ/a6ZfcXMHjOzT6bdHgAYdIkGhJk9aGav\nmtlP6rZvNLOjZnbMzHZKkru/4O53S/pTSeuTbBcAoLWk7yAekrSxdoOZFSQ9IOkWSWskbTWzNZX3\nPijpCUkHE24XAKCFRAPC3Z+S9Frd5hskHXP3l9z9rKSHJd1a2f9xd79F0seSbBcAoLVLUvjMUUnH\na16fkHSjmf2xpC2SLlOTOwgz2y5puyQtX778nddee21yLQWAPnT48OFfuPuKVvulERBB7v59Sd+P\nsN9eSXslaXx83KemppJtGAD0GTN7Ocp+afRimpG0qub1lZVtAIAMSSMgnpX0NjO7xswulfQRSY/H\nOYGZbTKzvadPn06kgQCA5Lu57pP0tKQxMzthZtvc/ZykeyRNSnpB0qPu/nyc87r7AXffPjw83P1G\nAwAkJVyDcPetDbYfFF1ZASDTMjOSOg4eMQFA8nIZEDxiAoDk5TIgAADJIyAAAEG5DAhqEACQvFwG\nBDUIAEheLgMCAJA8AgIAEJTLgKAGAQDJy2VAdFqDmJie0frdh3TNzie0fvchTUwzVyAA1MvMdN+9\nMjE9o137j2hufkGSNDM7p137j0iSNq8bTbNpAJApubyD6MSeyaOL4VA1N7+gPZNHU2oRAGTTwAXE\nydm5WNsBYFDlMiA6KVKvHCnF2g4AgyqXAdFJkXrHhjGVioUl20rFgnZsGOtW8wCgLwxckbpaiN4z\neVQnZ+e0cqSkHRvGKFADQJ2BCwipHBIEAgA0l8tHTACA5BEQAIAgAgIAEJTLgGAuJgBIXi4DgvUg\nACB5uQwIAEDyCAgAQBABAQAIIiAAAEEEBAAgiIAAAATlMiAYBwEAyctlQDAOAgCSl8uAAAAkj4AA\nAAQREACAIAICABBEQAAAgggIAEAQAQEACCIgAABBBAQAIIiAAAAE5TIgmIsJAJKXy4BgLiYASN4l\naTcAYRPTM9ozeVQnZ+e0cqSkHRvGtHndaNrNAjBACIgMmpie0a79RzQ3vyBJmpmd0679RySJkADQ\nM7l8xNTv9kweXQyHqrn5Be2ZPJpSiwAMIgIig07OzsXaDgBJICAyaOVIKdZ2AEgCAZFBOzaMqVQs\nLNlWKha0Y8NYSi0CMIgoUmdQtRBNLyYAaSIgMmrzulECAUCqeMQEAAgiIAAAQQQEACCIgAAABBEQ\nAIAgAgIAEJSpbq5mtlnS+yW9SdLX3P3fU24SAAysxO8gzOxBM3vVzH5St32jmR01s2NmtlOS3H3C\n3e+SdLekO5JuGwCgsV48YnpI0sbaDWZWkPSApFskrZG01czW1Ozymcr7AICUJB4Q7v6UpNfqNt8g\n6Zi7v+TuZyU9LOlWK/uipO+6+49C5zOz7WY2ZWZTp06dSrbxADDA0ipSj0o6XvP6RGXbpyS9V9Lt\nZnZ36EB33+vu4+4+vmLFiuRbCgADKlNFane/X9L9abcDAJDeHcSMpFU1r6+sbIvEzDaZ2d7Tp093\nvWEAgLK0AuJZSW8zs2vM7FJJH5H0eNSD3f2Au28fHh5OrIH9bGJ6Rut3H9I1O5/Q+t2HNDEdOZsB\nDJBedHPdJ+lpSWNmdsLMtrn7OUn3SJqU9IKkR939+aTbgnI47Np/RDOzc3JJM7Nz2rX/CCEB4CKJ\n1yDcfWuD7QclHUz687HUnsmjmptfWLJtbn5BeyaPsv4EgCVyOdUGNYj2nZydi7UdwODKZUBQg2jf\nypFSrO0ABlcuAwLt27FhTKViYcm2UrGgHRvGUmoRgKzK1DgIJK9aZ9gzeVQnZ+e0cqSkHRvGqD8A\nuEguA8LMNknatHr16rSbkkub140SCABayuUjJmoQAJC8XAYEACB5BAQAIIgaBNoyMT1DoRvoc7m8\ng6AGkS6m6wAGQy4DAulqNl0HgP5BQCA2pusABgMBgdiYrgMYDLkMCCbrSxfTdQCDIZcBQZE6XZvX\njereLWs1OlKSSRodKeneLWvpxQT0mVx2c0X6mK4D6H+5vIMAACSPgAAABPGICT3D6GsgXwgI9ER1\n9HV1gF119LUkQgLIqFw+YqKba/4w+hrIn1wGBN1c84fR10D+5DIgkD+Mvgbyh4BATzD6GsgfitTo\niWohml5MQH4QEOiZOKOv6RILpI+AQObQJRbIBmoQyBy6xALZECsgzOxyM3u7mb3VzFILF8ZB9De6\nxALZ0PKPvJkNm9lfm9kRSc9I+idJj0p62cy+ZWY3J93IeoyD6G90iQWyIcpdwGOSjkv6A3cfc/d3\nufu4u6+StFvSrWa2LdFWYqDQJRbIhpZFand/X5P3Dks63NUWYeBF7RJLTycgWZF7MZnZ37r7Z2te\nFyR93d0/lkjLMNBadYmlpxOQvDiF5lVmtkuSzOwySfslvZhIq4AW6OkEJC9OQPy5pLWVkDgg6Ul3\n/3wirQJaoKcTkLyWj5jM7Pqal/ep3IvpB5KeMrPr3f1HSTUOaGTlSEkzgTCo7enUrEZB/QJoLUoN\n4kt1r1+XtKay3SW9u9uNAlrZsWFsSQ1CWtrTqVmNQhL1CyCCKL2Yej7OAWilVU+nVjWKRu8REMAF\nUR4x3Snpm+7uDd7/HUm/7e7/0e3GAc006+nUTo2C+gWwVJRHTG+W9GMzq455OCXpNyStlvRHkn4h\naWdiLQwws02SNq1evbqXH4scaVWjaFW/ABChF5O73yfpekn7JK2Q9J7K6xlJH3f329y9p91dmWoD\nrTQbjc1IbSCaSAPl3H1B0vcq/wCZF2U0di96MdFbCnlmDUoLF3Yw+2yTt93d/667TYpufHzcp6am\n0vp4oKn6nlRS+U7l3i1rCQmkyswOu/t4q/2i3EG8Edi2XNI2lesTqQUEkKZWdwfNelIREMiDKN1c\nF8dBmNlvSvoLSZ+Q9LAuHiMBDIQoc0Ex2ht5F2mqDTO7wsy+IOm/VA6V6939r9z91URbB2RUlLmg\nWNcCeRdlwaA9kp6V9H+S1rr759399cRbBmRYlLsDeksh76LUID4t6deSPiPpb8ysut1ULlK/KaG2\nAZlRX28YWVbU62fmL9qv9u4g6roWQFZFqUGktvY0kAWhekNxyFQsmOYXLvQCDN0dtFrXAsiyyAsG\nAYMqVG+YP+8aKRW1/LJLErs7YAwF0kZAAC00qjecnpvXjz/3J4l8JivmIQsICKBG6Fd7lLUnmh3f\nzh90xlAgC6gvABXVX+0zs3NyXfjVfvO1KyL1Rmp0/MT0TOy2MIYCWUBAABWNfrU/+bNTunfLWo2O\nlGSSRkdKwekyurlONmMokAU8YgIqmv1qj9IbKcqv/vpHUDdfu0JP/uyUTs7OabhUlJk0e2Zew6Xi\nRb2kTOW7kvW7Dy0+uor6SCupgjeF9P5GQAAVcWoN7RwfKjx/45lXFvebnZtf8u/FIdPllfEWpvL6\nvtXjdu0/oqmXX9O3D8+0LGQnVfCmkN7/MvOIyczeamZfM7PH0m4LBlOnI59bHR96BNXM/HnXsksv\n0ehISfVzLs/NL+gbz7wSfKT1l4/8WOt3H1qsfTR69PXpR5+LVB+ZmJ7R+t2HdM3OJyKdt51Hasim\nRO8gzOxBSR+Q9Kq7v6Nm+0ZJ90kqSPqqu+9295ckbSMgkJZORz63Or6dAnO7RenaX/ONzrHg3vIX\nf7O7BArp/a/lehAdndzsDyX9StLXqwFhZgVJP5f0PkknVJ7naau7/7Ty/mPufnuU87MeBPJk/e5D\nwUdQzYw2WSK1W8ePjpT0g53vDr7XqM3NztvsfMiGqOtBJPqIyd2fkvRa3eYbJB1z95fc/azK04bf\nGvWcZrbdzKbMbOrUqVNdbC2QrNAjqGaaLZEa1cnZuZbHN/rFPzE90zBYZmbn9Mavz6lYsCXbmYyw\nv6RRgxiVdLzm9QlJo2b2ZjP7iqR1Zrar0cHuvtfdx919fMWKFUm3FeiazetGL+oue+dNVy2+HikV\ndfmy4kVdaWuPi2vlSGnx+IJZw33qVR8tNTM7Ny+5gm1Gf8hMLyZ3/6Wku9NuB5Ckdifvqx7XaBnT\n2945uqRHU3V79dd89TNDx4Z+8UctqFcL6dOfTWbKEaQrjYCYkbSq5vWVlW2RmdkmSZtWr17dzXYB\nmdesED5+9RVNC+xxivBxCs0UpftXokVqSTKzt0j6Tk2R+hKVi9TvUTkYnpX0UXd/Pu65KVIDzbU7\nkC1OQZ2idP5kokhtZvskPS1pzMxOmNk2dz8n6R5Jk5JekPRoO+EAoLlO5oYKFbara2DUoijd3xJ9\nxOTuWxtsPyjpYJKfDQy6TmaEbfQ4KrSNonT/ykyROg5qEEBrnQ5ka1RQJxAGRy4Dwt0PSDowPj5+\nV9ptAbKmWndoVF2snRuq3buBKMdW95mZnVPBTAvuGq2boJC7kGxLvEidJIrUwFKhbrC1SsWC7t2y\nVlK4y2uUcQyNutrWHtuqHaE2ERK9k4kiNYDeajZ+oXYgWycT7UU5Ns7EhEzwl125fMREDQIIa1Rf\nMGlJV9RO6hNRjo07NoKxFNmUyzsIdz/g7tuHh4fTbgqQKVFXomu035DZRdN6t/MZcVe+Y6W8bMpl\nQAAIi7qmRaMJ/BbcW46ZiPIZcSYYZCxFdhEQQB8JTQgYKgDX7xeayK9RbSDKZ9RPMFg9f/0EhUzw\nl2257MVUU4O468UXX0y7OUDuXbPziWC3WJP037vf3+vmIGF93YuJGgTQXVFrFxgsuezFBKC7dmwY\nC45bOHP2nCamZ5o+Amo0IK7ZALiJ6Rl9/vHny2tKSFpWHNJlxYJmz8wzeC5DCAgAi3+Ma/9oS9Lr\nZ+abrltdPyBuofLIunbt6tAI6x3fek7z5y881Dozf15n5s+3PBa9lctHTAC6b/O6US2/7OLfjM0G\nsjUbENfouD2TR5eEQ5xj0Vu5DAgz22Rme0+fPp12U4C+EncAXasBbqH3ow6KY/Bc+nIZEBSpgWTE\nLVa3KmKH3o9a+KZAnr5cBgSAZEQdaNds/1bH7dgwpuLQxeMuon4meociNYBFcdatrt8/ai+mUEGc\nXkzZlMuBclVM9w10Ls7aDr1YwyH0WVLvV7Lr5XfutagD5QgIYIC1u7ZDUms4hD6rOGSSSfMLF/5W\nJb2GRC+/cxr6eiQ1gO5od22HpLqhhj5r/rwvCYckP79ZOwax620uA4JurkB3dLK2QxLdUOOcM8lu\nsL38zlmWy4CgmyvQHZ2s7ZBEN9Q450yyGyxzU5XlMiAAdEe7azsk1Q019FnFIVOxsLRbbNLdYHv5\nnbOMbq7AAIvSrTVu19ck2tOrz2/Vjn4oUMdBLyYAGDBRezFxBwGg6+KMZWg1XXjtuYZLRZlpcUDd\nzdeu0Hee+5/FAXeXLyvqc5ve3vCzorZ10O4UGuEOAkBXxRnLcNs7R/XtwzPBGWFbvd/IkEmFIYs0\nbqLfxzs0wjgIAKmIM5Zh3w+PN50uvNn7jZx3RR43wXiH5nIZEIyDALIrzliBhRZPMFq9H0ecqccH\nbbxDI7kMCMZBANkVZ6xAwZrP6trq/TjiTD0+aOMdGsllQADIrjhjGbbeuKrpdOHN3m9kyBR53ATj\nHZqjFxOAroo7lmH86iua9mKqvp9ELybGOzRHLyYAGDD0YgIAdISAAAAEERAAgCCK1AC6ptNpK5pN\nuyFdWPvaJFWrp8uK5d+5Z+bPSyr3Yjrv0kipqDd+Pa/KZklSwcrHnXfJVO7tdLYyqK5a4J56+TV9\n85lXFs+//NKC/v5Daxen/ahdS7v280YrRfMnf3aqbwreFKkBdEWn01aEjq8KTdXRS4Uh09YbVumR\n/zyu+fPR25DVaTsoUgPoqU6nrQgdXxWaqqOXFs679v0wXjhI+Z+2I5cBwVQbQPZ0Om1F1qe3aHfa\nj6x/r2ZyGRBMtQFkT6fTVmR9eot2p/3I+vdqJpcBASB7Op22InR8VWiqjl4qDJm23riqXAuJIe/T\ndhAQALpi87pR3btlrUZHSjKVe/XEKdDWHi9d+MU+OlLSng9fpz23X7f4Xu2f6WXFocWeTFK5V5FU\n7sVUrPsLV7AL75ukS2tC5/JlRf3DHb+nO2+6asn5l19a0Jc+fJ2+sHmt9nz4Oo2UikvOWT3f6EhJ\nd950VdvfP4voxQQAA4ZeTACAjhAQAIAgAgIAEERAAACCCAgAQBABAQAIIiAAAEEEBAAgiIAAAAQR\nEACAIAICABBEQAAAgjKzJrWZLZf0j5LOSvq+u38z5SYB6JJO16qOcs7qetCh9azrP6t+benqutK1\na13XrlH9jWdeWXJ8db+RUlFnzy0srofdSHFIarbLkEkfvfEqSdK+Hx7XgrsKVp5i/Aub1y75vlG+\nX7ckOpurmT0o6QOSXnX3d9Rs3yjpPkkFSV91991m9nFJs+5+wMwecfc7Wp2f2VyB7Ot0reqo52yk\n/rMmpme041vPxV4+NC133nSVxq++ouH3bee/ZVZmc31I0sbaDWZWkPSApFskrZG01czWSLpS0vHK\nbq3/rwPIhU7Xqo56zkbqP2vP5NHchINUvqNo9n2TXPc60YBw96ckvVa3+QZJx9z9JXc/K+lhSbdK\nOqFySDRtl5ltN7MpM5s6depUEs0G0EWdrlXdjWNr98/bGtEL7i3bnNR3SqNIPaoLdwpSORhGJe2X\ndJuZfVnSgUYHu/tedx939/EVK1Yk21IAHet0repuHFu7f97WiC6YtWxzUt8pM72Y3P0Nd/+Eu3+S\nAjXQPzpdqzrqORup/6wdG8Ziry2dpq03rmr6fZNc9zqNXkwzklbVvL6ysi0yM9skadPq1au72S4A\nCagWT7vZiyl0zqi9mKr/nrdeTNXv2ze9mCTJzN4i6TvVXkxmdomkn0t6j8rB8Kykj7r783HPTS8m\nAIgvE72YzGyfpKcljZnZCTPb5u7nJN0jaVLSC5IebSccAADJSvQRk7tvbbD9oKSDSX42AKAzmSlS\nx2Fmm8xs7+nTp9NuCgD0rVwGhLsfcPftw8PDaTcFAPpWLgMCAJC8XAYEj5gAIHm5DAgeMQFA8nIZ\nEACA5BEQAIAgAgIAEJTLgKBIDQDJy2VAUKQGgOTlMiAAAMkjIAAAQQQEACAojQWDOlZdMEjS/5rZ\ni5XNw5IaVa0bvfdbkn7R/RZ2VbPvlYVzt3OOOMdE2bfVPlwb6Z0/z9dHs+Oyfn20+s5XRzqLu/fF\nP5L2xn1P0lTa7e7ke2Xh3O2cI84xUfZttQ/XBtdHO++3uG4yfX1069rop0dMB9p8L+uSbHs3zt3O\nOeIcE2XfVvtwbaR3/jxfHwN/bSS+5GiWmdmUR1h2D4OHawPNDMr10U93EO3Ym3YDkFlcG2hmIK6P\ngb6DAAA0Nuh3EACABggIAEAQAQEACCIgKsxsuZn9i5n9s5l9LO32IFvM7K1m9jUzeyzttiBbzGxz\n5e/GI2b2J2m3p5v6OiDM7EEze9XMflK3faOZHTWzY2a2s7J5i6TH3P0uSR/seWPRc3GuD3d/yd23\npdNS9FrMa2Oi8nfjbkl3pNHepPR1QEh6SNLG2g1mVpD0gKRbJK2RtNXM1ki6UtLxym4LPWwj0vOQ\nol8fGCwPKf618ZnK+32jrwPC3Z+S9Frd5hskHav8Ijwr6WFJt0o6oXJISH3+3wVlMa8PDJA414aV\nfVHSd939R71ua5IG8Q/hqC7cKUjlYBiVtF/SbWb2ZeV7iD06E7w+zOzNZvYVSevMbFc6TUPKGv3t\n+JSk90q63czuTqNhScnlbK5JcPc3JH0i7XYgm9z9lyo/YwaWcPf7Jd2fdjuSMIh3EDOSVtW8vrKy\nDZC4PtDYwF0bgxgQz0p6m5ldY2aXSvqIpMdTbhOyg+sDjQzctdHXAWFm+yQ9LWnMzE6Y2TZ3Pyfp\nHkmTkl6Q9Ki7P59mO5EOrg80wrVRxmR9AICgvr6DAAC0j4AAAAQREACAIAICABBEQAAAgggIAEAQ\nAQF0kZm9pX6KaCCvCAgAQBABASSksgrdtJn9ftptAdrBbK5AAsxsTOX1Av7M3Z9Luz1AOwgIoPtW\nSPo3SVvc/adpNwZoF4+YgO47LekVSe9KuyFAJ7iDALrvrKQPSZo0s1+5+7+m3SCgHQQEkAB3f8PM\nPiDpe5WQ6Ot1A9CfmO4bABBEDQIAEERAAACCCAgAQBABAQAIIiAAAEEEBAAgiIAAAAQREACAoP8H\nbowHrDaBsvwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc22c86a390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "degree_dist = final_network.degree_distribution()\n",
    "xs, ys = zip(*[(left, count) for left, _, count in degree_dist.bins()])\n",
    "matplotlib.pyplot.scatter(xs, ys, marker=\"o\")\n",
    "ax = matplotlib.pyplot.gca()\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xscale(\"log\")\n",
    "matplotlib.pyplot.ylim((0.5,1000))\n",
    "pylab.xlabel(\"k\")\n",
    "pylab.ylabel(\"N(k)\")\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `igraph.statistics.power_law_fit` function to estimate the scaling exponent *alpha* of the degree distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.7453577551868795"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "igraph.statistics.power_law_fit(final_network.degree()).alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## extra challenge:\n",
    "If you got this far, see if you can scatter plot the relationship between R (as the independent variable) and -log10(P) value (as the dependent variable).  When the effect size variable (e.g., *R*) can range from negative to positive, this plot is sometimes called a \"volcano plot\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGg1JREFUeJzt3X2wXVWZ5/HvkxAgIDYBAgUJIQyiNgiGnhTgMD0yUjRI\niyBt89Jgo1LSNYXTQ2NT8tYDdovQQxuwukdnYGAERV5U+so02FRawS6dgBW4gQhMJLwoHJGgEHmL\nCOGZP86+eHJz7nm59+zz+v1Unbr7rL3PyTo7yfndtdfaa0VmIknSZLN6XQFJUn8yICRJdRkQkqS6\nDAhJUl0GhCSpLgNCklSXASFJqsuAkCTVZUBIkuraotcVmImddtopFy9e3OtqSNJAuffee3+RmfOb\nHTfQAbF48WJWrlzZ62pI0kCJiJ+0cpyXmCRJdRkQkqS6DAhJUl0GhCSpLgNCklTXQI9ikqRBNjZe\n4bI71vCz9RvYbfu5nH3EOzj2gAW9rtabDAhJ6oGx8Qrn3rKaDa9tBKCyfgPn3rIaoG9CorRLTBGx\ndUT8MCLuj4gHI+IzRfmeEXFPRKyNiJsiYsuifKvi+dpi/+Ky6iZJvXbZHWveDIcJG17byGV3rOlR\njTZXZh/Eq8D7MvPdwBLgyIg4GPhb4PLMfBvwPHBacfxpwPNF+eXFcZI0lH62fkNb5b1QWkBk1UvF\n0znFI4H3Ad8oyq8Fji22jymeU+w/LCKirPpJUi/ttv3ctsp7odRRTBExOyJWAeuA5cCjwPrMfL04\n5Clg4mLbAuBJgGL/r4Ad67zn6RGxMiJWPvvss2VWX5JKMTZe4eVXX9+sfO6c2Zx9xDt6UKP6Sg2I\nzNyYmUuAhcCBwDs78J5XZubSzFw6f37TuaYkqa+MjVc4++v3s37Da5uUz9tmDpcct1/fdFBDl+6D\nyMz1wJ3Ae4DtI2Ji9NRCoFJsV4DdAYr9vwP8shv1k6RuuejWB3ntjdysPLN/Ri9NKHMU0/yI2L7Y\nngscDjxMNSg+XBx2KvCtYvvW4jnF/u9m5uZnUZIG2OSWQ7PyXirzPohdgWsjYjbVILo5M/8pIh4C\nboyIzwLjwNXF8VcDX4mItcBzwIkl1k2S1ERpAZGZDwAH1Cl/jGp/xOTyXwN/XFZ9JKkfzNtmDs+/\nsnlrYd42c3pQm8a8k1qSumBiWo164TBndnDh0fv2oFaNGRCSVLLJ02oABNUbwxb04RxMEwwISSpZ\nvWk1JsLhB+e8rzeVaoHTfUtSyQZhWo16DAhJKtkgTKtRjwEhSSU7+4h3MHfO7E3K+m1ajXrsg5Ck\nkk10QPfz4kD1GBCSVJJ+XzGuGQNCkkowCCvGNWMfhCSV4LxbHuj7FeOaMSAkqcMOX3YXr7z2Rt19\n/T60tZYBIUkddMHYah5Z9/KU+/t9aGstA0KSOuiGe55suL/fh7bWMiAkqYM2NljGJmJwOqjBgJCk\njhkbrzTcf/JBi7pUk84wICSpAyaGtU5l75235bPH7tfFGs2cASFJHVBvxtYJpxy8iOVnHdrdCnWA\nASFJHTDV8NWAgWs5TDAgJKkDBnXG1kYMCEnqgEGdsbUR52KSpA4Y1BlbGzEgJKlDjj1gwUAHwmQG\nhCRN06BP592MASFJ03DB2Gquv/unTNw3PYjTeTdjJ7UktWlsvLJJOEwYtOm8mzEgJKlNl92xZrNw\nmDBI03k3Y0BIUpsahcAg3/cwmQEhSW0YG68wK6LuvmCwpvNuxk5qSWrR5I7pWgGcfPCioemgBgNC\nklpywdhqvnr3T+vumx3B549/91CFA3iJSZKaGhuvTBkOAG9kDl04gAEhSU2dd8sDDfcPU8d0LQNC\nkhoYG6/wymtvTLl/2DqmaxkQktTAp7/ZuPUwbB3TtQwISZrC4cvu4tXXp249zJ0za2AXA2qFASFJ\ndYyNV3hk3csNj7nkuP27VJveKC0gImL3iLgzIh6KiAcj4r8U5RdFRCUiVhWPo2pec25ErI2INRFx\nRFl1k6Rmms2pdMoQX1qaUOZ9EK8Dn8rM+yJiO+DeiFhe7Ls8M/+u9uCI2Ac4EdgX2A34l4h4e2bW\nXwVckkoyNl6h0mROpWG+tDShtIDIzKeBp4vtFyPiYaBR3B4D3JiZrwKPR8Ra4EBgRVl1lKTJJu6W\nbuSQvXboUm16qyt9EBGxGDgAuKco+mREPBAR10TEvKJsAfBkzcueonGgSFJHTdwQN9VMrQB777wt\n13/iPV2rUy+VHhAR8Rbgm8CZmfkC8CVgL2AJ1RbG59t8v9MjYmVErHz22Wc7Xl9Jo6vZDXFXnLCE\n5Wcd2p3K9IFSAyIi5lANh+sz8xaAzHwmMzdm5hvAVVQvIwFUgN1rXr6wKNtEZl6ZmUszc+n8+fPL\nrL6kEdLshrgF288d+k7pycocxRTA1cDDmbmspnzXmsM+BPyo2L4VODEitoqIPYG9gR+WVT9JqtVs\n1NKw3i3dSJmjmA4BPgKsjohVRdl5wEkRsQRI4AngzwAy88GIuBl4iOoIqDMcwSSpWxotArTNnFkj\n13qAckcxfZ/qNCWT3d7gNRcDF5dVJ0mqZ2IRoI1Zv3v6c0N+Q9xUXA9C0kgbtUWA2mFASBpZY+OV\nKcNhWBcBaodzMUkaWZfdsWbKex6GdRGgdtiCkDRyGl1WmjCsiwC1w4CQNFIarS09YZgXAWqHl5gk\njZRWwmGUO6Zr2YKQNDJOvqr53J+Xn7DEcCjYgpA0Mn7w6HMN94/idBqNGBCSRsLY+GZTu23GfodN\nGRCSht7YeIVzb1nd8JhRWCGuXfZBSBpqJ1+1oumlpb133nYkVohrlwEhaWgdvuwuHln3csNjDtlr\nh5FZAKhdXmKSNJTGxitNw2HB9nMNhwYMCElD6aJbH2y4f+6c2XZKN2FASBo6J1+1gvUbXmt4zCXH\n7WendBMGhKSh0kqn9CF77WA4tMCAkDQ0Wh2xZL9DawwISUOhlXCYt80clp91aHcqNAQMCEkDr5Vw\nALjw6H27UJvhYUBIGmithoP9Du0zICQNrLHxSsvhYL9D+wwISQPrrJtWNT3GcJg+A0LSQDro4uW8\n0eQYw2FmDAhJA+fkq1bwzIu/aXiM4TBzBoSkgXLB2Oqm/Q5bzg7DoQOczVXSwGh1xNJ/+/C7u1Cb\n4ddSQETEzsAhwG7ABuBHwMrMbHYJUJI6opWpu6F6p7TDWTujYUBExH8EzgF2AMaBdcDWwLHAXhHx\nDeDzmflC2RWVNLouGFvdUjjsst2W3indQc1aEEcBn8jMn07eERFbAB8ADge+WULdJAmAr9692VfQ\nZuyU7ryGAZGZZzfY9zow1vEaSVKNk69a0fSYbebMMhxK0HAUU0QcFBH3R8RLEbEiIvbpVsUkjbax\n8QqHXPrdljqlP3fc/l2o0ehpNsz1vwN/CewILAMuL71Gkkbe2HiFc29ZTWX9hqbHnnLwIjulS9Ks\nD2JWZi4vtr8eEeeWXSFJo21svMJf3LyKzObHXnHCEsOhRM0CYvuIOG6q55l5SznVkjSKLhhb3VKH\nNDg7azc0C4jvAUdP8TwBA0JSR4yNV9oKBzuly9dsFNPHulURSaPtsjvWND3GS0rd1WwU0ykRMeUx\nEbFXRPz7zldL0ii5YKx5h/S8beYYDl3W7BLTjsB4RNwL3As8S/VO6rcB7wV+QfVO681ExO7AdcAu\nVC9HXZmZX4iIHYCbgMXAE8Dxmfl8RATwBao3570CfDQz75vRp5PU18bGK5x3ywO88lrzWXtcLrT7\nGrYgMvMLwO8BNwDzgcOK5xXgI5n5R5n5yBQvfx34VGbuAxwMnFHcR3EO8J3M3Bv4Dr8NmPcDexeP\n04EvzeSDSepvF4yt5sybVrUUDnZI90bTyfoycyOwvHi0LDOfBp4utl+MiIeBBcAxwKHFYdcCdwGf\nLsqvy8wE7o6I7SNi1+J9JA2RVjukI+Dkgxbx2WP360KtNFmzyfq2AE6jOjnfRHxXgG8BV2fma638\nIRGxGDgAuAfYpeZL/+dUL0FRvP+TNS97qigzIKQhc/4/rm56zOwIHr3kqC7URlNp1oL4CrAe+AzV\nL2yAhcCpwFeBE5r9ARHxFqqT+Z2ZmS9UuxqqMjMjooXbYTZ5v9OpXoJi0aJF7bxUUo+NjVe47I41\nvPybjU2PPemg3btQIzXSLCD+bWa+fVLZU1QvAf242ZtHxByq4XB9zU11z0xcOoqIXalOIQ7Vlknt\nv4iFRdkmMvNK4EqApUuXthUuknqn1cV+oDp9hpeVeq/ZXEzPRcQf1w51jYhZEXEC8HyjFxajkq4G\nHs7MZTW7bqXaAqH4+a2a8j+NqoOBX9n/IA2HdsLhihOWGA59olkL4kTgb4EvRsREIGwP3Fnsa+QQ\n4CPA6ohYVZSdB1wK3BwRpwE/AY4v9t1OdYjrWqrDXL1JTxoC7YSDo5X6S7M7qZ+g6GeIiB2Lsl+2\n8saZ+X0gpth9WJ3jEzijlfeWNBgOung5z7z4m6bHzY7gpIN2t+XQZ1pakxo2D4aIOLxmpldJ2kSr\n4TBvmzmM/9c/6EKN1K5mfRCNXN2xWkgaGmPjFfb5q2+3FA7gHdL9rNl9ELdOtYvqNByS9Kax8Qpn\nf+N+XtvY2gBDF/vpb80uMf0+cArw0qTyAA4spUaSBtLYeIUzb1rV/MCCQ1n7X7OAuBt4JTO/N3lH\nRDSfm1fSSDh82V08su7llo4N4HKn7R4IzUYxvb/Bvv/Q+epIGjTthMMu223JPecfXnKN1Cktj2KS\npMlaHakEXlIaRC0FRES8SHVNh1q/AlZSndL7sU5XTFL/Ghuv8Bc3rdrsS2EqrgQ3mFptQVxBdQ6m\nr1G9hHgisBdwH3ANv52+W9KQu2BsdctrRwPsvfO2hsOAavU+iA9m5v/MzBcz84ViwrwjMvMmYF6J\n9ZPUR6YTDsvPOrS8CqlUrQbEKxFxfDFR36yIOB74dbHPGVWlEdDqIj8TTjl4keEw4Fq9xHQy1fWi\nv1g8XwGcEhFzgU+WUTFJ/eGCsdVcf/dPW/5NcFbAsuPtcxgGLQVE0Ql99BS7v9+56kjqJ+0MYQWH\nsQ6bli4xRcTCiPjHiFhXPL4ZEQvLrpyk3mk3HE45eJHhMGRa7YP431QX9NmtePyfokzSkBkbr/C7\nf/XttsLhkL128B6HIdRqH8T8zKwNhC9HxJllVEhS77Tb3wDVcLj+E+8prU7qnVYD4pcRcQpwQ/H8\nJKClhYMkDYZ2h7Buu+VsLv7QfnZGD7FWA+LjwN8Dl1Md1vp/cUlQaWjsf+E/88KrG1s+3vsbRkOr\no5h+Anyw5LpI6rJ2p8wA51QaJc0WDPp7GtwIl5l/3vEaSeqKdkYpzQr4k4MMhlHTrAWxsk7ZrsDT\nJdRFUhe0u7DPNnNm8dDfTDnzv4ZYs/Ugrp1cFhH3ZebvlVclSWV55/m38+sWlwOF6jj4zx23f3kV\nUl9r9T6IWtHxWkgq1dh4hcXn3NZWOAAsc5rukTadBYOu6ngtJJWm3UtKAG/dajYPfObIkmqkQdF2\nQGTmF5sfJakftDtdBnjjm37LJUelITSdYACHsGpTBoQ0RNq9G3pCAI9f+oedr5AGmgEhDYnF59w2\nrdd5SUlTMSCkAWerQWUxIKQBZqtBZTIgpAE03VYDwBO2GtQiA0IaINO5p2GCl5TULgNCGhAzaTU4\nfFXTYUBIfe7kq1bwg0efm9ZrXbdBM2FASH3sbefexuvtTZ/0JlsNmikDQupD7a7wVstWgzrFgJD6\nyEz6GcBWgzqrtICIiGuADwDrMvNdRdlFwCeAZ4vDzsvM24t95wKnARuBP8/MO8qqm9SPZtJqMBhU\nhjJbEF8G/gG4blL55Zn5d7UFEbEPcCKwL7Ab8C8R8fbMnN7/FmmAzCQYHLqqMpUWEJn5rxGxuMXD\njwFuzMxXgccjYi1wILCipOpJPdfu6m6TeTe0ytaLPohPRsSfUl3v+lOZ+TywALi75pinirLNRMTp\nwOkAixYtKrmqUufNZNgqwNazg/938VEdrJFUX7cD4kvA3wBZ/Pw88PF23iAzrwSuBFi6dOn0f/2S\numymHdAGg7qtqwGRmc9MbEfEVcA/FU8rwO41hy4syqShMN1J9SZc4drQ6oGuBkRE7JqZTxdPPwT8\nqNi+FfhaRCyj2km9N/DDbtZNKsNMg8FOaPVSmcNcbwAOBXaKiKeAC4FDI2IJ1UtMTwB/BpCZD0bE\nzcBDwOvAGY5g0iCbaTB4s5v6QWQO7mX8pUuX5sqVK3tdDelNM5kaY4LTcatsEXFvZi5tdpx3Uksd\nsOc5tzHTX7Uctqp+Y0BIMzDTS0kAu2y3Jfecf3gHaiN1lgEhTcNMh6yCw1bV/wwIqQ2duJQE9jNo\nMBgQUgs6cSkJnFRPg8WAkBroxKgkMBg0mAwIqY5OtRjsgNYgMyCkGp0KBm900zAwICRsMUj1GBAa\nWYcvu4tH1r3ckfdyziQNIwNCI6dTrQUwGDTcDAiNjE4GA3gvg4afAaGh18lg2CJg7SUGg0aDAaGh\nNNNlPSez81mjyIDQUOn0ZSRvcNMoMyA0FOxfkDrPgNDA6uQw1QkGg/RbBoQGTqdbC2/dajYPfObI\njr6nNAwMCA2MTgcD2GKQGjEg1NfKCAVbDFJrDAj1JVsLUu8ZEOobZYQCGAzSdBkQ6qn9L/xnXnh1\nY8ff1/WepZkzINQTthak/mdAqGvKCgUwGKQyGBAq1TvPv51fb+zAos51GApSuQwIlcLWgjT4DAh1\nTJmh4KR5UvcZEJqRMkMBbC1IvWRAqG2GgjQaDAi1ZM9zbqOcrubfMhik/mJAaEqdXpWtHkNB6l8G\nhDZhKEiaYEColIV36jEYpMFiQIyoMm9gq2UoSIPLgBghZY8+mmAoSMOhtICIiGuADwDrMvNdRdkO\nwE3AYuAJ4PjMfD4iAvgCcBTwCvDRzLyvrLqNEkNB0nSV2YL4MvAPwHU1ZecA38nMSyPinOL5p4H3\nA3sXj4OALxU/NQ2GgqROKC0gMvNfI2LxpOJjgEOL7WuBu6gGxDHAdZmZwN0RsX1E7JqZT5dVv2Fj\nKEjqtG73QexS86X/c2CXYnsB8GTNcU8VZQbEFLrVyQyGgjSqetZJnZkZEW1/w0XE6cDpAIsWLep4\nvfpZt1oJYChI6n5APDNx6SgidgXWFeUVYPea4xYWZZvJzCuBKwGWLl3anV+he8hQkNQr3Q6IW4FT\ngUuLn9+qKf9kRNxItXP6V6Pa/9DNQABDQdLUyhzmegPVDumdIuIp4EKqwXBzRJwG/AQ4vjj8dqpD\nXNdSHeb6sbLq1Y8MBUn9qMxRTCdNseuwOscmcEZZdek33Q4EMBQktc87qbvEUJA0aAyIkhgIkgad\nAdEhvQgEMBQklceAmKZeBQIYCpK6w4Bo0dvOvY3Xe3TXhYEgqRcMiBbYnyBpFBkQTTgJnqRRZUD0\niIEgqd8ZEF1kKEgaJAZEiQwESYPMgOggA0HSMDEgmnji0j+csqPaQJA0zAyIFhgEkkbRrF5XQJLU\nnwwISVJdBoQkqS4DQpJUlwEhSaorqqt9DqaIeJbq2tbdshPwiy7+eYPAc1Kf52VznpP6enFe9sjM\n+c0OGuiA6LaIWJmZS3tdj37iOanP87I5z0l9/XxevMQkSarLgJAk1WVAtOfKXlegD3lO6vO8bM5z\nUl/fnhf7ICRJddmCkCTVZUAAEXFkRKyJiLURcU6d/ZdHxKri8eOIWF+zb2PNvlu7W/NytXBeFkXE\nnRExHhEPRMRRNfvOLV63JiKO6G7NyzPdcxIRiyNiQ82/lf/R/dqXp4XzskdEfKc4J3dFxMKafadG\nxCPF49Tu1rw8Mzwn/fG9kpkj/QBmA48C/wbYErgf2KfB8f8ZuKbm+Uu9/gy9Oi9Ur53+p2J7H+CJ\nmu37ga2APYv3md3rz9Tjc7IY+FGvP0MPz8vXgVOL7fcBXym2dwAeK37OK7bn9foz9fKcFM/74nvF\nFgQcCKzNzMcy8zfAjcAxDY4/CbihKzXrrVbOSwJvLbZ/B/hZsX0McGNmvpqZjwNri/cbdDM5J8Os\nlfOyD/DdYvvOmv1HAMsz87nMfB5YDhzZhTqXbSbnpG8YELAAeLLm+VNF2WYiYg+qvxF/t6Z464hY\nGRF3R8Sx5VWz61o5LxcBp0TEU8DtVFtXrb52EM3knADsWVx6+l5E/H6pNe2uVs7L/cBxxfaHgO0i\nYscWXzuIZnJOoE++VwyI9pwIfCMzN9aU7ZHVuyD/BLgiIvbqTdV64iTgy5m5EDgK+EpEjPq/qanO\nydPAosw8ADgL+FpEvLXB+wybvwTeGxHjwHuBCrCx8UuGXqNz0hffK6P+nxmqfym71zxfWJTVcyKT\nLi9lZqX4+RhwF3BA56vYE62cl9OAmwEycwWwNdV5Zdo5p4Nk2uekuNz2y6L8XqrXp99eeo27o+l5\nycyfZeZxRUCeX5Stb+W1A2om56R/vld63QnS6wfVZVcfo3rpaKIzad86x70TeILi3pGibB6wVbG9\nE/AIDTq4B+nRynkBvg18tNj+XarX2wPYl007qR9jODqpZ3JO5k+cA6odlxVgh15/pi6el52AWcX2\nxcBfF9s7AI8X/5fmFdsDf15meE765nul5yeyHx5ULwX8mOpvdecXZX8NfLDmmIuASye97t8Bq4u/\n/NXAab3+LN08L1Q72X5QfP5VwB/UvPb84nVrgPf3+rP0+pwAfwQ8WJTdBxzd68/S5fPy4eKL7sfA\n/5r4Aiz2fZzqQIa1wMd6/Vl6fU766XvFO6klSXXZByFJqsuAkCTVZUBIkuoyICRJdRkQkqS6tuh1\nBaRhEhEbqQ5N3ILqmP6PZHHzkzRobEFInbUhM5dk5ruA54Azel0haboMCKk8KxiOiec0ogwIqQQR\nMRs4DBiqRaQ0WgwIqbPmRsQq4OfALlTXN5AGkgEhddaGzFwC7EF1kj77IDSwnItJ6qCIeCkz31Js\nHwCMAXtl5uu9rZnUPlsQUkkycxx4gOoiQtLAsQUhSarLFoQkqS4DQpJUlwEhSarLgJAk1WVASJLq\nMiAkSXUZEJKkugwISVJd/x8yYTkuXMYhhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc22c7717b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "inds_use = numpy.where(P_values_adj > 0)\n",
    "matplotlib.pyplot.scatter(cor_coeff_values_above_thresh[inds_use], -numpy.log10(P_values_adj[inds_use]))\n",
    "pylab.xlabel(\"R\")\n",
    "pylab.ylabel(\"-log10(P)\")\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## extra-extra challenge\n",
    "For each of the gene pairs for which R>0.75, see if you can compute the t-test P value for each correlation coefficient (don't bother adjusting for false discovery rate control).  Compare to the (un-adjusted) P values that you got using the Fisher transformation, using a scatter plot.  How do they compare?  Which test has better statistical power, for this case where M = 414? (If you are wondering, *general* advice is to use Fisher if M>=10; for very small numbers of samples, use the Student t test)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHeFJREFUeJzt3XucXWV97/HPN5MBJlwc0iAHcukgjfhCuTqFULRNRURS\nNUCtgHgE6yGnHrQNXoPwMngsBUVRrKfYWBGoEcIljigqRg4pFkloYAIDKBLllgEJHiUgRC7hd/5Y\nzyabyZ6ZvVdm7ev3/XrtV9Y8e+29n1ls8s1zWc+jiMDMzGykSY2ugJmZNScHhJmZVeSAMDOzihwQ\nZmZWkQPCzMwqckCYmVlFDggzM6vIAWFmZhU5IMzMrKLJja7Atpg2bVr09fU1uhpmZi3ltttu+01E\n7DbeeS0dEH19faxZs6bR1TAzaymSHqzmPHcxmZlZRQ4IMzOryAFhZmYVOSDMzKwiB4SZmVXU0rOY\nzMw6zcDgMOdffy+PPLGJPXt7+NhR+3DMQdML+SwHhJlZixgYHOaM5UNsen4zAMNPbOKM5UMAhYSE\nu5jMzFrE+dff+1I4lGx6fjPnX39vIZ/ngDAzaxGPPLGppvJt5S4mM7MmM9o4w569PQxXCIM9e3sK\nqYcDwsysSQwMDvPp797N7555/qWy8nGGjx21z8vGIAB6urv42FH7FFIfdzGZmTWB0gB0eTiUlMYZ\njjloOucetx/Te3sQML23h3OP28+zmMzM2lmlAehypXGGYw6aXlggjOQWhJlZExhvoLmocYaxuAVh\nZlYnY93kNtoANBQ7zjAWtyDMzOqgNMYw/MQmgi2DzwODw0A2AN3T3bXV63p7ugsdZxiLWxBmZnUw\n1k1u5eMK9VpGoxoOCDOzOqjmJrd6DkBXw11MZmZ1MNogcyMGn6vlgDAzq4NKYwyNGnyulruYzMy2\nQbXLbzfjGMN4CgsISTOBy4DdgQCWRMSFkqYCy4A+4AHgXRHxO0kCLgTmAc8Ap0TE7UXVz8xsW9W6\n/HazjTGMp8gupheAj0TEvsAc4DRJ+wKLgBsiYjZwQ/oZ4GhgdnosAC4qsG5mZrmdNTDEXouuY+Gy\ntXVdfrveCmtBRMSjwKPp+ClJPwOmA/OBuem0S4GVwCdS+WUREcAqSb2S9kjvY2bWcAODw5z57SGe\nfm70JTGguOW3660ug9SS+oCDgNXA7mV/6f+arAsKsvB4uOxl61OZmVnDlbqTxgsHaO6ZSbUoPCAk\n7QRcAyyMiCfLn0uthajx/RZIWiNpzeOPPz6BNTUzG914i+mVNPvMpFoUOotJUjdZOCyNiOWp+LFS\n15GkPYANqXwYmFn28hmp7GUiYgmwBKC/v7+mcDEzq9bI2UmjrZNUrktq2LIYRSisBZFmJX0d+FlE\nXFD21LXAyen4ZOA7ZeXvVWYOsNHjD2bWCJXWTdI4r+maJL7wrgPaJhyg2BbE4cB/B4YkrU1lnwTO\nA66U9H7gQeBd6bnvk01xXUc2zfV9BdbNzGxUlbqTAhCV+8R33K6Lc45tn5ZDSZGzmP4TRg3dIyqc\nH8BpRdXHzKxao81CCrJd3FrlRrdt5TupzcxGGG3MYXpvDzcvelMDatQYXovJzGyEVlw3qQhuQZiZ\njdCK6yYVwQFhZm2v2gX1yrXauklFcECYWdsaGBzm09+9m9898/xLZeMtqGdbeAzCzNpS6V6G8nAo\naacF9YrkgDCztjTe0hjtsqBekdzFZGZtodalMdplQb0iOSDMrKWNNs4w2l3P0JlTVvNwQJhZyzpr\nYIilqx6qGASjLY3R29PN2e94rQeoq+CAMLOWMzA4zNnX3s0Tm7YegC7XaUtjTDQHhJm1lLFaDSN1\n2tIYE80BYWZN76yBIS5f/TCbo/otYDzOsO0cEGbW1E762i3c/Mvf1vQajzNMDAeEmTWtswaGagoH\nASfNmcU/HrNfcZXqIA4IM2tKA4PDLF31UNXn7zqlm8Vvd6thIjkgzKxplN/sNkmqaiDarYbiOCDM\nrKFKoTDy5rZqBqSne+pqoRwQZtYwA4PDfOzqO3h+cxYG1c5ROnzvqSw99bDiKmaAA8LMGmRgcJgP\nX7mWF6ufuerupDorbDVXSRdL2iDprrKyAyWtkrRW0hpJh6RySfqypHWS7pR0cFH1MrPGKy3FXU04\ndEmIrDvpi8cf6HCooyJbEJcAXwEuKyv7HPDpiPiBpHnp57nA0cDs9DgUuCj9aWZtaLyluEt6urs4\n97j9PMbQIIUFRETcJKlvZDGwSzp+BfBIOp4PXBYRAayS1Ctpj4h4tKj6mVl9VNrus5q9GDwA3Xj1\nHoNYCFwv6fNk3Vt/lsqnAw+Xnbc+lTkgzFrUWNt99k7prrjTG0D3JHH+3xzgYGgC9d5R7gPA6REx\nEzgd+HqtbyBpQRq/WPP4449PeAXNbNuNt91nRNZ9NFJvT7fDoYnUOyBOBpan46uAQ9LxMDCz7LwZ\nqWwrEbEkIvojon+33XYrrKJmlt94YwwbNz3Pucftx/TenpcGoL90/IGsXfwWh0MTqXcX0yPAXwAr\ngTcB96Xya4EPSrqCbHB6o8cfzFpHtfszlOzZ28MxB013GDS5wgJC0uVkM5SmSVoPLAZOBS6UNBn4\nA7Agnf59YB6wDngGeF9R9TKziXXWwBDfrGHNJC/D3TqKnMV04ihPvb7CuQGcVlRdzGzi1RoM4GW4\nW43vpDazmtUaDp6y2pocEGZWtfKF9arlbT9blwPCzKpSmrpazR3QJd1d8nhDC3NAmNmo8uwFXeIN\nfFqfA8LMtjIwOMyZ3x7i6eeqby2UvMerrbYNB4SZvUyeriSAHbfr4pxjvbBeO3FAmBmQbwDaq622\nNweEmXHWwBBLVz1U9Y5ugpdWZnU4tC8HhFmHO+lrt3DzL39b9flTuidxz2eOLrBG1iwcEGYdKO/s\npEmCfzpu/4JqZc3GAWHWYWptMXRJbI7w3dAdyAFh1iHyjDN88fgDHQgdzAFh1gFqbTUAnDRnlsOh\nwzkgzNrYWQNDfGv1Q7xY443Qh+891Te7mQPCrB0NDA7z8avv4LnNtSVDl8SJh850OBjggDBrOwOD\nwyxctram1xy+91SWnnpYQTWyVuWAMGsjtYaDl8ewsTggzFpcrftBl7jVYONxQJi1qIHBYT65/E6e\nef7Fml7X0z2Jc4/b360GG5cDwqwFDQwO85Gr7mBzjdOT3GqwWhQWEJIuBt4GbIiI15WVfwg4DdgM\nXBcRH0/lZwDvT+V/HxHXF1U3s1aVd9oqOBysdkW2IC4BvgJcViqQ9JfAfOCAiHhW0itT+b7ACcBr\ngT2BH0t6dUTUvluJWZs6a2CIb656qObXTRK8+1Bv4mO1KywgIuImSX0jij8AnBcRz6ZzNqTy+cAV\nqfx+SeuAQ4BbiqqfWavIs08DQG9PN2e/w1t+Wn71HoN4NfBGSecAfwA+GhH/BUwHVpWdtz6VmXW0\nvLu7uTvJJkK9A2IyMBWYA/wpcKWkV9XyBpIWAAsAZs2aNeEVNGsGeVsNIltDyd1JNhHqHRDrgeUR\nEcCtkl4EpgHDwMyy82aksq1ExBJgCUB/f3+OoTqz5pV3nGHXKd0sfru7k2xijRsQkrYvjRmMVVal\nAeAvgRslvRrYDvgNcC3wLUkXkA1SzwZuzfH+Zi0rTzj4ngYrUjUtiFuAg6soexlJlwNzgWmS1gOL\ngYuBiyXdBTwHnJxaE3dLuhK4B3gBOM0zmKzTXL764arO6+nu4tzjvDyGFW/UgJD038gGinskHUTW\nvQmwCzBlvDeOiBNHeeo9o5x/DnDOeO9r1k5q3frTu7pZPY3VgjgKOIVsPOALbAmIJ4FPFlsts/aW\nZ5mMXad0c/OiNxVYK7OXGzUgIuJS4FJJfx0R19SxTmZtK++d0F2TxOK3v7aYSpmNYtwxCIeD2cTw\nDCVrNV6sz6xOvrW6unDokvjlufMKro3Z+Oo9zdWso9Q6CA1w4qEzxz/JrA4Km+Zq1snydie9x3dB\nWxMpbJqrWScaGBzm41ffwXObaxuF7u4S57/zAI8zWFOpdprrBWXlT+FprmZbOfKCldy34emaX+dB\naGtWnuZqNgFqDQcBe/qmN2ty1YxBfE/Su4G+8vMj4n8XVSmzVpEtx30nm2q44c3jDNYqqgmI7wAb\ngdsAz1wyS/IMRB++91SHg7WMagJiRkS8tfCamLWIk752Czf/8rc1vaZ7Epz/Nwe6O8laSjUB8VNJ\n+0XEUOG1MWtiA4PDfPSqO3ihxnUy3KVkraqagHgDcIqk+8m6mAREROxfaM3MmkieVkOX4AvvcqvB\nWlc1AXF04bUwa1J572twq8HaQTWL9T0o6Q3A7Ij4hqTdgJ2Kr5pZ4+RZjhtgx+26OOdYb+Zj7aGa\ntZgWA/3APsA3gG7gm8DhxVbNrDG8TIZZppoupmOBg4DbASLiEUk7F1orswbJczf0JOCC4z3WYO2n\nmoB4LiJCUgBI2rHgOpnVXd4ZSofvPZWlpx5WUK3MGmtSFedcKelfgV5JpwI/Br5WbLXM6mdgcJiP\nXZ1v+qrDwdpZNYPUn5d0JNle1PsAn4qIFeO9TtLFwNuADRHxuhHPfQT4PLBbRPxGkoALgXnAM8Ap\nEXF7zb+NWQ3yTF0Ftxqsc1S1o1wKhHFDYYRLgK8Al5UXSpoJvAUoHwU8GpidHocCF6U/zQpx6Dkr\neOyp52p+nQeirZOMtR/EU8Cobe6I2GWsN46ImyT1VXjqi8DHydZ4KpkPXBYRAayS1Ctpj4h4dKzP\nMMvjyAtW1hQOvuHNOtVYy33vDCDpM8CjwL+T3UV9ErBHng+TNB8Yjog7sl6ll0wHHi77eX0qc0DY\nhMmz8ursV+7Iig/PLa5SZk2smi6md0TEAWU/XyTpDuBTtXyQpClkGw29pZbXVXifBcACgFmzZm3L\nW1mHyHPTW29PN2e/w5v4WGerJiCelnQScAVZl9OJQO3bZsHewF5AqfUwA7hd0iHAMFC+U/uMVLaV\niFgCLAHo7++vbdqJdRTf8Ga2baoJiHeTzTC6kCwgbk5lNUmrwb6y9LOkB4D+NIvpWuCDkq4gG5ze\n6PEH2xZ5B6G9X4PZFtVMc32AbBC5JpIuB+YC0yStBxZHxNdHOf37ZFNc15FNc31frZ9nBll30unL\n1o4+u2IUHog221o1azF9gwqzmSLib8d6XUScOM7zfWXHAZw2Xl3MxpLnvoZJgncf6i4ls0qq2pO6\n7HgHsrWZHimmOmb5vObM7/MHL8ltNqGq6WK6pvzn1HX0n4XVyKwG+y/+IU8+u7mm10zpnsQ/Hbe/\nu5PMxlHVndQjzKZssNmsEQYGh1m4bG3Nr3Orwax61YxBjLyj+tfAJwqrkdk48obDl7wkt1lNquli\n8t4P1hQGBof5yJVrqXGogd133o7VZx5ZTKXM2lg1LYgbIuKI8crMipRnrAG88qrZthhrsb4dgClk\n9zHsSrYOE8AuZOskmRUu793QAu4/768mvkJmHWSsFsT/BBYCewK3sSUgniRbxtusUH2Lrsv1Oi+w\nZzYxxlrN9ULgQkkfioh/rmOdrMPlHYSeLFh3rlsNZhNlrC6mPwUeLoWDpPcCfw08CJwdEbVvxWU2\njrytBk9fNZt4Y3Ux/SvwZgBJfw6cB3wIOJBsNdV3Fl476xh5xxoAHvBYg1khxgqIrrJWwvHAknRX\n9TWSam//m40i7wylHbrEz8+ZV0CNzAzGCQhJkyPiBeAI0iY9VbzOrGp5upQEfNE3vZkVbqy/6C8H\n/kPSb4BNwE8AJP0JsLEOdbM2lrdLyTOUzOpnrFlM50i6gWz/6R+lJbkBJpGNRZjl4oFos9YwZldR\nRKyqUPaL4qpj7SzPfg3gsQazRvFYgtWFb3ozaz0OCCvUkRes5L4NT9f8Oq+hZNZ4DggrTJ5d3sDL\ncps1CweETbi8rYZdtu/izk+/tYAamVkek4p6Y0kXS9og6a6ysvMl/VzSnZK+Lam37LkzJK2TdK+k\no4qqlxWrb9F1ucLhPXNmORzMmkxhAQFcAoz8P34F8LqI2B/4BXAGgKR9gROA16bX/IukrgLrZhNs\nYHA490D0A+f9laevmjWhwrqYIuImSX0jyn5U9uMqtqznNB+4IiKeBe6XtA44BLilqPrZxMk71uCd\n3syaWyPHIP4WWJaOp5MFRsl6vClRS9iWVoOZNbeGBISkM4EXgKU5XruAtC7UrFmzJrhmVq1Dz1nB\nY089V/Pr3Gowax11DwhJpwBvA44oW75jGJhZdtqMVLaViFhCttw4/f39tfdr2DbL22rw9FWz1lLX\ngJD0VuDjwF9ExDNlT10LfEvSBWRbnM4Gbq1n3Wx8eZfKAHcpmbWiwgJC0uXAXGCapPXAYrJZS9sD\nKyQBrIqIv4uIuyVdCdxD1vV0WkTUvkGAFSZvq8F3RJu1Lm3p5Wk9/f39sWbNmkZXo63l3R8a3Gow\na1aSbouI/vHO853UNqq8A9GTBevOdTiYtToHhFXkgWgzc0DYy+Td6Q3cpWTWbhwQ9pK8rQZ3KZm1\nJweEeSDazCpyQHS4vOsoeWlus/bngOhgXkfJzMbigOhAf3LGdbyQ8/YXh4NZ5yhyPwhrQn2L8oXD\n7jtv53Aw6zBuQXQQdymZWS0cEB3AXUpmlocDos251WBmeTkg2tT+i3/Ik8/mWxDX4WBm4IBoS3lb\nDd7tzczKOSDajLuUzGyiOCDaxF6LriPPOLSA+x0OZlaBA6INuNVgZkXwjXItzuFgZkVxC6JF5Q0G\ncDiYWXUcEC3IrQYzq4fCupgkXSxpg6S7ysqmSloh6b70566pXJK+LGmdpDslHVxUvVqdw8HM6qXI\nMYhLgJEbBiwCboiI2cAN6WeAo4HZ6bEAuKjAerWkvkXX5QoH4XAws3wK62KKiJsk9Y0ong/MTceX\nAiuBT6TyyyIigFWSeiXtERGPFlW/VuJWg5k1Qr1nMe1e9pf+r4Hd0/F04OGy89anso7ncDCzRmnY\nIHVEhKSa7+2StICsG4pZs2ZNeL2aRd5gOHzvqSw99bAJro2ZdaJ6B8Rjpa4jSXsAG1L5MDCz7LwZ\nqWwrEbEEWALQ39+fcxHr5uZWg5k1g3p3MV0LnJyOTwa+U1b+3jSbaQ6wsRPHH84aGHI4mFnTKKwF\nIelysgHpaZLWA4uB84ArJb0feBB4Vzr9+8A8YB3wDPC+ourVrPIGw5eOP5BjDvJwjZlNvCJnMZ04\nylNHVDg3gNOKqkuzc6vBzJqR12JqMIeDmTUrL7XRIA4GM2t2bkE0gMPBzFqBA6LOHA5m1ircxVQn\nDgYzazVuQdSBw8HMWpEDomAOBzNrVe5iKojXUjKzVueAKIBbDWbWDtzFNMEcDmbWLtyCmCB7LbqO\nvEvLOhzMrBk5ICaAWw1m1o7cxbSNHA5m1q4cENvA4WBm7cxdTDk4GMysE7gFUSOHg5l1CgdEDRwO\nZtZJ3MVUhYHBYRYuW1vz6xwMZtbKHBDjcKvBzDqVu5jG4HAws07WkICQdLqkuyXdJelySTtI2kvS\naknrJC2TtF0j6lbicDCzTlf3LiZJ04G/B/aNiE2SrgROAOYBX4yIKyR9FXg/cFG96+dgMDPLNKqL\naTLQI2kyMAV4FHgTcHV6/lLgmHpXyuFgZrZF3QMiIoaBzwMPkQXDRuA24ImIeCGdth6YXs96ORzM\nzF6uEV1MuwLzgb2AJ4CrgLfW8PoFwAKAWbNmTUid8oSDg8HM2l0jprm+Gbg/Ih4HkLQcOBzolTQ5\ntSJmAMOVXhwRS4AlAP39/XlX2AbcajAzG0sjxiAeAuZImiJJwBHAPcCNwDvTOScD3ymyEg4HM7Ox\nNWIMYjXZYPTtwFCqwxLgE8CHJa0D/gj4elF1cDiYmY2vIXdSR8RiYPGI4l8BhxT92R5vMDOrTkfd\nSe1wMDOrXkcFRK0cDmbWyRwQo3A4mFmn82quIzgYzMwybkGUcTiYmW3RUQExVgA4HMzMXq7jupgc\nBGZm1emoFoSZmVXPAWFmZhU5IMzMrCIHhJmZVeSAMDOzihSxTVsqNJSkx4EHG12PZBrwm0ZXosn4\nmlTm67I1X5PKiroufxwRu413UksHRDORtCYi+htdj2bia1KZr8vWfE0qa/R1cReTmZlV5IAwM7OK\nHBATZ0mjK9CEfE0q83XZmq9JZQ29Lh6DMDOzityCMDOzihwQVZJ0saQNku4qK5sqaYWk+9Kfu6Zy\nSfqypHWS7pR0cONqXpxRrsnZkoYlrU2PeWXPnZGuyb2SjmpMrYslaaakGyXdI+luSf+Qyjv9uzLa\ndenY74ukHSTdKumOdE0+ncr3krQ6/e7LJG2XyrdPP69Lz/cVXsmI8KOKB/DnwMHAXWVlnwMWpeNF\nwGfT8TzgB4CAOcDqRte/jtfkbOCjFc7dF7gD2B7YC/gl0NXo36GAa7IHcHA63hn4RfrdO/27Mtp1\n6djvS/pvvlM67gZWp+/AlcAJqfyrwAfS8f8CvpqOTwCWFV1HtyCqFBE3Ab8dUTwfuDQdXwocU1Z+\nWWRWAb2S9qhPTetnlGsymvnAFRHxbETcD6wDDimscg0SEY9GxO3p+CngZ8B0/F0Z7bqMpu2/L+m/\n+e/Tj93pEcCbgKtT+cjvSuk7dDVwhCQVWUcHxLbZPSIeTce/BnZPx9OBh8vOW8/Y/zO0mw+m7pKL\nS10pdOA1SV0AB5H9y9DflWTEdYEO/r5I6pK0FtgArCBrKT0RES+kU8p/75euSXp+I/BHRdbPATFB\nImv3eUoYXATsDRwIPAp8obHVaQxJOwHXAAsj4sny5zr5u1LhunT09yUiNkfEgcAMshbSaxpcpZdx\nQGybx0rdAenPDal8GJhZdt6MVNb2IuKx9KV/EfgaW7oFOuaaSOom+0twaUQsT8Ud/12pdF38fclE\nxBPAjcBhZN2Mpd0+y3/vl65Jev4VwP8rsl4OiG1zLXByOj4Z+E5Z+XvTDJU5wMay7oW2NqL//Fig\nNMPpWuCENBNjL2A2cGu961e01Cf8deBnEXFB2VMd/V0Z7bp08vdF0m6SetNxD3Ak2djMjcA702kj\nvyul79A7gf+bWqPFafRIfqs8gMvJmsDPk/ULvp+s/+8G4D7gx8DU2DI74f+Q9ScOAf2Nrn8dr8m/\np9/5zvSF3qPs/DPTNbkXOLrR9S/omryBrPvoTmBteszzd2XU69Kx3xdgf2Aw/e53AZ9K5a8iC8N1\nwFXA9ql8h/TzuvT8q4quo++kNjOzitzFZGZmFTkgzMysIgeEmZlV5IAwM7OKHBBmZlaRA8JakqTN\nZSuArpXUJ6lf0pfHeM1cSd+bwDrMlfRnE/V+Y3zOa9LvOChp76I/z6xk8vinmDWlTZEtUVDuAWBN\nUR8oaXJsWSMHYC7we+CnVZy7LY4Bro6If6zm5HRTmiK7O9ksN7cgrG2UtxAk/UVZ62JQ0s7ptJ0k\nXS3p55KWllbDlPR6Sf8h6TZJ15cti7FS0pckrQH+oeyz+oC/A05Pn/FGSZdI+qqk1cDnJB0i6Zb0\n+T+VtE967SmSlkv6obL9IT6XyrvSe9wlaUjS6cr2R1gIfEDSjem8D6dz7pK0sFQfZfsmXEZ209VM\nSb+XdL6yvQZ+nOqzUtKvJL2j4P8c1g4afTehH37keQCb2XJH7rdT2Vzge+n4u8Dh6XgnstbyXLIV\nMGeQ/ePoFrI7fLvJWgG7pfOPBy5OxyuBfxmlDmdTtpcBcAnwPdK+BcAuwOR0/GbgmnR8CvArsrV0\ndgAeJFtj5/XAirL36x35OemcIWDH9HvdTbYyah/wIjCn7PVBugMZ+Dbwo/S7HgCsbfR/Qz+a/+Eu\nJmtVlbqYyt0MXCBpKbA8ItanxsKtEbEeIC2z3Ac8AbwOWJHO6SJbQqRkWQ31uioiNqfjVwCXSppN\n9pd1d9l5N0TExlSPe4A/JvvL/lWS/hm4juwv9JHeQBaIT6fXLgfeSLZMxYOR7SlR8hzww3Q8BDwb\nEc9LGkq/t9mY3MVkbSkizgP+B9AD3CyptIzys2WnbSZrWQi4OyIOTI/9IuItZec9XcNHl5/7GeDG\niHgd8Hay1kLJVvWIiN+R/et+JVn31b/V8LmV6vl8RJTW0nmx9JmRjU34H4c2LgeEtSVJe0fEUER8\nFvgvxl5n/15gN0mHpdd2S3ptFR/zFNn2maN5BVuWaj6lijpPAyZFxDXAWWTbuY70E+AYSVMk7Ui2\nAupPqqirWc0cENauFqZB3DvJVpv9wWgnRsRzZMsnf1bSHWTjGtVMX/0ucGxpkLrC858DzpU0SHX/\nYp8OrExdX98EzqhQ19vJxjpuJduR7d8iYrCK9zarmVdzNTOzityCMDOzihwQZmZWkQPCzMwqckCY\nmVlFDggzM6vIAWFmZhU5IMzMrCIHhJmZVfT/AeLDt3euiWgYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc22c418dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ts = numpy.divide(cor_coeff_values_above_thresh * math.sqrt(M - 2), numpy.sqrt(1 - cor_coeff_values_above_thresh**2))\n",
    "P_values_studentT = 2*scipy.stats.t.cdf(-ts, M-2)\n",
    "inds_use = numpy.where(numpy.logical_and(P_values > 0, P_values_studentT > 0))\n",
    "matplotlib.pyplot.scatter(-numpy.log10(P_values[inds_use]),\n",
    "                          -numpy.log10(P_values_studentT[inds_use]))\n",
    "pylab.xlabel(\"Fisher transform\")\n",
    "pylab.ylabel(\"Student t\")\n",
    "pylab.show()"
   ]
  }
 ],
 "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
}