COMP90014__Assignment_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COMP90014 Assignment 2\n",
    "### Semester 2, 2019\n",
    "\n",
    "**Task 6 has been ammended - please use this version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This assignment should be completed by each student individually. Make sure you read this entire document, and ask for help if anything is not clear. Any changes or clarification to this document will be announced via the LMS.\n",
    "\n",
    "Please make sure you are aware of the University's rules on academic honesty and plagiarism, which are very strict: https://academichonesty.unimelb.edu.au/ \n",
    "\n",
    "Make sure you **do not** copy any code either from other students or from the internet. This is considered plagiarism. It is generally a good idea to avoid looking at any solutions as you may find it surprisingly difficult to generate your own solution to the problem once you have seen somebody else's.\n",
    "\n",
    "Your completed notebook file containing all your answers will be turned in via LMS. No other files or formats will be accepted - only upload the completed `.ipynb` file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overview\n",
    "To complete the assignment you will need to finish the tasks in this notebook. There are multiple tasks that are connected in a logical order.\n",
    "\n",
    "The tasks are a combination of writing your own implementations of algorithms we've discussed in lectures, writing your own code to use library implementations of these algorithms and interpreting the results in short answer format. Each short answer question has a word limit that will be strictly enforced! **Please note that for this assignment, you may be awarded a mark of zero for a question if you go over the word limit.**\n",
    "\n",
    "In some case, we have provided test input and test output that you can use to try out your solutions. These tests are just samples and are not exhaustive - they may warn you if you've made a mistake, but they are not guaranteed to. It's up to you to decide whether your code is correct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marking\n",
    "\n",
    "Cells that must be completed to receive marks are clearly labeled. There are 19 graded cells, some of which are code cells, in which you must complete the code to solve a problem, and some of which are markdown cells, in which you must write your answers to short-answer questions. \n",
    "\n",
    "In this assignment, every graded cell is worth 2 marks. In addition to the graded cells, up to 7 marks will be given for code style, readability, efficiency and comments. \n",
    "\n",
    "The total marks for the assignment add up to 45, and it will be worth 15% of your overall subject grade.\n",
    "\n",
    "Please make sure that you do not edit the \"GRADED CELL\" comments in either the code or the markdown cells, as this will disrupt the marking system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background and data \n",
    "\n",
    "WGCNA stands for weighted gene co-expression network analysis. It is a data analysis technique used for studying biological networks based on pairwise correlations of gene expression data. WGCNA is good at identifying clusters of genes that may be co-regulated, and therefore may have shared biological function.\n",
    "\n",
    "For this assignment, you will primarily be using the [FlyAtlas](http://flyatlas.org) dataset. For this assignment, instead of using the probe-wise dataset, we will be using the expression value for each gene.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 0 - Setup "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import altair"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import scipy\n",
    "import re\n",
    "from io import StringIO\n",
    "from sklearn.decomposition import PCA\n",
    "from scipy.cluster.hierarchy import linkage, dendrogram, fcluster\n",
    "from scipy.spatial.distance import squareform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_expression = pd.read_csv('flyatlas_subset.csv.gz', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_expression.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data frame has 18952 rows (measurements) and 136 columns (samples) so it is certainly high dimensional.These 136 columns represent 4 replicates each from 34 different tissue types."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will set numpy and pandas to display numbers to just two decimal places in this notebook - this won't affect the actual numbers, just their display, and you can change it if you prefer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=2)\n",
    "pd.options.display.precision=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The actual stored numbers have not changed\n",
    "raw_expression.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_expression = np.log(raw_expression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code snippet removes the replicate name from each sample, so we can use these labels as categories for plotting later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tissues_list = [re.match('(.+?)(( biological)? rep\\d+)', c).group(1)\n",
    "                     for c in raw_expression.columns]\n",
    "tissues = pd.Series(tissues_list, index = raw_expression.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(tissues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1 - Building a correlation matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1(a)**\n",
    "\n",
    "The [FlyAtlas](http://flyatlas.org) dataset contains four biological replicates for each tissue. Combine the biological replicates by calculating the mean expression value for each gene in each tissue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def average_by_tissue(expression, tissues):\n",
    "    '''\n",
    "    Given a DataFrame of gene expression data, \n",
    "    and a list, array or Series of tissues corresponding to the columns of the dataframe,\n",
    "    average over the expression values in each gene for each tissue type and\n",
    "    return the resulting dataframe. \n",
    "    The columns of the new dataframe should correspond to the provided tissues.\n",
    "    '''\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return\n",
    "\n",
    "```\n",
    "      A    B\n",
    "0   4.5  2.5\n",
    "1  10.0  7.0\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df = pd.DataFrame([[5,4,3,2],[10,10,6,8]])\n",
    "print(average_by_tissue(test_df, ['A','A','B','B']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate expression for each tissue in the flyatlas data\n",
    "tissue_expression = average_by_tissue(log_expression, tissues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tissue_expression.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1(b)**\n",
    "\n",
    "WGCNA starts by building a pairwise correlation matrix of genes. Using the matrix you just created, produce an *unsigned* correlation matrix where each cell contains the absolute value of the correlation coefficients.\n",
    "\n",
    "You can calculate the Pearson correlation values yourself, or look up a numpy or scipy function to do so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def calculate_unsigned_correlation(expression):\n",
    "    ''' \n",
    "    Produce the unsigned correlation matrix for a table of gene expression values.\n",
    "    Assume that the columns of the expression matrix are samples and the rows are\n",
    "    genes, and return an array of arrays giving the Pearson correlation between each pair of genes,\n",
    "    in the same order as the rows of the expression table.\n",
    "    '''\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 1.  ,  0.95,  0.96,  0.44,  0.3 ,  0.15],\n",
    "       [ 0.95,  1.  ,  1.  ,  0.71,  0.59,  0.46],\n",
    "       [ 0.96,  1.  ,  1.  ,  0.67,  0.54,  0.41],\n",
    "       [ 0.44,  0.71,  0.67,  1.  ,  0.99,  0.95],\n",
    "       [ 0.3 ,  0.59,  0.54,  0.99,  1.  ,  0.99],\n",
    "       [ 0.15,  0.46,  0.41,  0.95,  0.99,  1.  ]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df = pd.DataFrame([[ 3.8,  2.7,  4.5],\n",
    "                       [ 4.3,  3.4,  6.2],\n",
    "                       [ 5.3,  4.3,  7. ],\n",
    "                       [ 4.6,  6. ,  7.7],\n",
    "                       [ 5.2,  7.3,  8.8],\n",
    "                       [ 6.2,  8.5,  9.4]], \n",
    "                         columns=['Tissue1', 'Tissue2', 'Tissue3'],\n",
    "                         index=['GeneA', 'GeneB', 'GeneC', 'GeneD', 'GeneE', 'GeneF'])\n",
    "calculate_unsigned_correlation(test_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 1.  ,  0.95,  0.3 ,  0.15],\n",
    "       [ 0.95,  1.  ,  0.59,  0.46],\n",
    "       [ 0.3 ,  0.59,  1.  ,  0.99],\n",
    "       [ 0.15,  0.46,  0.99,  1.  ]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df = pd.DataFrame([[ 3.8,  2.7,  4.5],\n",
    "                       [ 4.3,  3.4,  6.2],\n",
    "                       [ 5.2,  7.3,  8.8],\n",
    "                       [ 6.2,  8.5,  9.4]], \n",
    "                         columns=['Tissue1', 'Tissue2', 'Tissue3'],\n",
    "                         index=['GeneA', 'GeneB', 'GeneC', 'GeneD'])\n",
    "calculate_unsigned_correlation(test_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the correlation matrix for the flyatlas data\n",
    "unsigned_correlation = calculate_unsigned_correlation(tissue_expression)\n",
    "_ = plt.hist(unsigned_correlation.flatten(), bins=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1(c)**\n",
    "\n",
    "Why are we using an unsigned correlation matrix instead of a signed correlation matrix? (max 50 words)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*# ~~ GRADED CELL - your answer here --*\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2 - Building an adjacency matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the correlation matrix to create a network, we will transform it into an adjacency matrix. You will create two types of adjacency matrix, a binary adjacency matrix and a weighted adjacency matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2(a)**\n",
    "\n",
    "To create the binary adjacency matrix, transform the correlation matrix such that every correlation greater than or equal to a given threshold value is considered adjacent (represented by a 1 in the matrix), and every correlation below that value is considered not adjacent (represented by a 0). Set the diagonal of the adjacency matrix to 0, so that we don't consider a node to be adjacent to itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def calculate_binary_adjacencies(correlation, threshold):\n",
    "    '''\n",
    "    Given a correlation matrix between genes of shape (N,N),\n",
    "    return the corresponding binary adjacency matrix of shape (N,N),\n",
    "    where correlation values are above the given threshold.\n",
    "    '''\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 0.,  1.,  0.,  0.],\n",
    "       [ 1.,  0.,  1.,  0.],\n",
    "       [ 0.,  1.,  0.,  1.],\n",
    "       [ 0.,  0.,  1.,  0.]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_corr = np.array([[ 1.  ,  0.95,  0.3 ,  0.15],\n",
    "       [ 0.95,  1.  ,  0.59,  0.46],\n",
    "       [ 0.3 ,  0.59,  1.  ,  0.99],\n",
    "       [ 0.15,  0.46,  0.99,  1.  ]])\n",
    "calculate_binary_adjacencies(test_corr, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 0.,  1.,  0.,  0.],\n",
    "       [ 1.,  0.,  0.,  0.],\n",
    "       [ 0.,  0.,  0.,  1.],\n",
    "       [ 0.,  0.,  1.,  0.]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_corr = np.array([[ 1.  ,  0.95,  0.3 ,  0.15],\n",
    "       [ 0.95,  1.  ,  0.59,  0.46],\n",
    "       [ 0.3 ,  0.59,  1.  ,  0.99],\n",
    "       [ 0.15,  0.46,  0.99,  1.  ]])\n",
    "calculate_binary_adjacencies(test_corr, 0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the binary adjacency matrix for the flyatlas data\n",
    "adjacency_binary = calculate_binary_adjacencies(unsigned_correlation, 0.85)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2(b)**\n",
    "\n",
    "Calculate the connectivity of the adjacency matrix by dividing the total number of edges by the number of possible edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def calculate_connectivity(adjacency):\n",
    "    '''\n",
    "    Calculate the number of edges that exist in a given binary adjacency matrix,\n",
    "    divided by the total number of possible edges between all nodes.\n",
    "    '''\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 0.5\n",
    "calculate_connectivity(np.array([[ 0.,  1.,  0.,  0.],\n",
    "                                   [ 1.,  0.,  1.,  0.],\n",
    "                                   [ 0.,  1.,  0.,  1.],\n",
    "                                   [ 0.,  0.,  1.,  0.]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return 0.33\n",
    "calculate_connectivity(np.array([[ 0.,  1.,  0.,  0.],\n",
    "                                   [ 1.,  0.,  0.,  0.],\n",
    "                                   [ 0.,  0.,  0.,  1.],\n",
    "                                   [ 0.,  0.,  1.,  0.]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calculate_connectivity(adjacency_binary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2(c)**\n",
    "\n",
    "The weighted adjacency matrix can be created by raising the correlation matrix to some power. Write a function that raises the correlation matrix to some power, `beta`, and sets the diagonal to `0`. For the rest of the assignment we will use `beta = 4` but your function should accept any integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def calculate_weighted_adjacencies(correlation, beta):\n",
    "    '''\n",
    "    Given a correlation matrix between genes of shape (N,N),\n",
    "    return the corresponding binary adjacency matrix of shape (N,N),\n",
    "    where we use a power-law soft threshold with parameter beta.\n",
    "    '''\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 0.  ,  0.9 ,  0.09,  0.02],\n",
    "       [ 0.9 ,  0.  ,  0.35,  0.21],\n",
    "       [ 0.09,  0.35,  0.  ,  0.98],\n",
    "       [ 0.02,  0.21,  0.98,  0.  ]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_corr = np.array([[ 1.  ,  0.95,  0.3 ,  0.15],\n",
    "       [ 0.95,  1.  ,  0.59,  0.46],\n",
    "       [ 0.3 ,  0.59,  1.  ,  0.99],\n",
    "       [ 0.15,  0.46,  0.99,  1.  ]])\n",
    "calculate_weighted_adjacencies(test_corr, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 0.  ,  0.86,  0.03,  0.  ],\n",
    "       [ 0.86,  0.  ,  0.21,  0.1 ],\n",
    "       [ 0.03,  0.21,  0.  ,  0.97],\n",
    "       [ 0.  ,  0.1 ,  0.97,  0.  ]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_corr = np.array([[ 1.  ,  0.95,  0.3 ,  0.15],\n",
    "       [ 0.95,  1.  ,  0.59,  0.46],\n",
    "       [ 0.3 ,  0.59,  1.  ,  0.99],\n",
    "       [ 0.15,  0.46,  0.99,  1.  ]])\n",
    "calculate_weighted_adjacencies(test_corr, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the weighted adjacency matrix for the flyatlas data\n",
    "adjacency_weighted = calculate_weighted_adjacencies(unsigned_correlation, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2(d)**\n",
    "\n",
    "How do you expect the network connectivity would change if the threshold for the binary adjacency matrix is increased or decreased? (max 50 words)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*# ~~ GRADED CELL - your answer here --*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3 - Defining modules with hierarchical clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will implement a distance or dissimilarity function between genes that makes use of the gene network as specified by the adjacency matrix, and these distances to carry out hierarchical clustering. We'll use scipy's hierarchical clustering functions `linkage()` and `fcluster()`, as they provide us with an easy way to draw the dendrogram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distance or dissimilarity function we'll implement is based on that given in lectures:\n",
    "\n",
    "If $i=j$, \n",
    "\n",
    "$$d_{ij} = 0$$\n",
    "\n",
    "otherwise\n",
    "\n",
    "$$d_{ij} = 1 - \\frac{l_{ij} + a_{ij}}{min(k_i,k_j) + 1 - a_{ij}}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$l_{ij} = \\Sigma_u a_{iu} a_{uj}$$\n",
    "\n",
    "$$k_i = \\Sigma_u a_{iu}$$\n",
    "\n",
    "and $a_{ij}$ refers to the $i$,$j$th element of the adjacency matrix.\n",
    "\n",
    "Note we have set $d_{ij}$ to $0$ if $i=j$ as this is the distance from a node to itself.\n",
    "\n",
    "In the functions below we'll refer to $k_i$ as the vertex connectivity of node $i$ (i.e. gene $i$), and $l_{ij}$ as the neighbour connectivity between $i$ and $j$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 3(a)**\n",
    "\n",
    "In the distance metric above, $l_{ij}$ sums over every node in the graph, but only a subset of nodes contribute. Based on the equations when applied to the binary adjacency matrix, explain which subset of nodes contribute to $d_{ij}$ for a given $i$ and $j$. (max 50 words)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*# ~~ GRADED CELL - your answer here --*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 3(b)**\n",
    "\n",
    "Complete the `vertex_connectivity()` and `neighbour_connectivity()` functions below. The `distance_matrix()` function has been provided, and should work correctly once the other two functions are complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def vertex_connectivity(adjacency):\n",
    "    '''\n",
    "    Given an adjacency matrix of shape (N,N), calculate the\n",
    "    vertex connectivity k_i of every node i and return these as an array of\n",
    "    shape (N,1).\n",
    "    '''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def neighbour_connectivity(adjacency):\n",
    "    '''\n",
    "    Given the adjacency matrix of shape (N,N), calculate the sum\n",
    "    of path weights from every i to every j via a single neighbouring node,\n",
    "    i.e. Sum a_iu a_uj over all nodes u.\n",
    "    Return these path weights from i to j in a matrix of shape (N,N)\n",
    "    '''\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_matrix(adjacency):\n",
    "    '''\n",
    "    Given the adjacency matrix of shape (N,N), calculate the distance \n",
    "    between every i and j based on the dissamilarity formula provided in lectures,\n",
    "    and return these distances in a matrix of shape (N,N).\n",
    "    '''\n",
    "    l_ij = neighbour_connectivity(adjacency)\n",
    "    k_i = vertex_connectivity(adjacency)[:,np.newaxis]\n",
    "    print(k_i)\n",
    "    k_j = vertex_connectivity(adjacency)[np.newaxis,:]\n",
    "    print(k_j)\n",
    "    print(np.minimum(k_i,k_j))\n",
    "    d_ij = 1 - ((l_ij + adjacency) / (np.minimum(k_i,k_j) + 1 - adjacency))\n",
    "    np.fill_diagonal(d_ij, 0)\n",
    "    return d_ij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([ 0.89,  1.17,  1.21,  1.07])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vertex_connectivity(np.array([[ 0.  ,  0.86,  0.03,  0.  ],\n",
    "                           [ 0.86,  0.  ,  0.21,  0.1 ],\n",
    "                           [ 0.03,  0.21,  0.  ,  0.97],\n",
    "                           [ 0.  ,  0.1 ,  0.97,  0.  ]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 0.74,  0.01,  0.18,  0.12],\n",
    "       [ 0.01,  0.79,  0.12,  0.2 ],\n",
    "       [ 0.18,  0.12,  0.99,  0.02],\n",
    "       [ 0.12,  0.2 ,  0.02,  0.95]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbour_connectivity(np.array([[ 0.  ,  0.86,  0.03,  0.  ],\n",
    "                           [ 0.86,  0.  ,  0.21,  0.1 ],\n",
    "                           [ 0.03,  0.21,  0.  ,  0.97],\n",
    "                           [ 0.  ,  0.1 ,  0.97,  0.  ]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below test case should return (if displayed to a precision of two decimal places)\n",
    "\n",
    "```\n",
    "array([[ 0.  ,  0.16,  0.89,  0.94],\n",
    "       [ 0.16,  0.  ,  0.83,  0.85],\n",
    "       [ 0.89,  0.83,  0.  ,  0.1 ],\n",
    "       [ 0.94,  0.85,  0.1 ,  0.  ]])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_matrix(np.array([[ 0.  ,  0.86,  0.03,  0.  ],\n",
    "                       [ 0.86,  0.  ,  0.21,  0.1 ],\n",
    "                       [ 0.03,  0.21,  0.  ,  0.97],\n",
    "                       [ 0.  ,  0.1 ,  0.97,  0.  ]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the distance matrix for flyatlas data\n",
    "distances = distance_matrix(adjacency_weighted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now carry out hierarchical clustering. Here, scipy's `linkage()` function performs agglomerative clustering based on the provided distances. \n",
    "\n",
    "The most important function below is `linkage()`. Other functions you see are `squareform()`, a utility function to transform the distance matrix into the format scipy requires, and `transpose()`, which is here used to compensate for small floating-point errors and ensure the distance matrix is exactly symmetric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = linkage(squareform((distances+distances.transpose())/2), 'ward')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy provides a function to draw the dendrogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(25, 10))\n",
    "_ = dendrogram(Z, no_labels=True)\n",
    "plt.axhline(y = 4.5, color = 'r', linestyle = '--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can extract the desired number of flat clusters with `fcluster()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = fcluster(Z, 6, criterion='maxclust')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a numpy array assigning each gene to one of 6 clusters, labelled 1 to 6:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 3(c)**\n",
    "\n",
    "Write a function to convert this flat list of labels into a set of modules in the form of lists of gene names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def module_lists(genes, labels):\n",
    "    '''\n",
    "    Given an array or series of gene names and an array or series of cluster labels,\n",
    "    return a list of lists where each list represents the genes in a cluster.\n",
    "    '''\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return [['GeneA', 'GeneC'], ['GeneB'], ['GeneD']]\n",
    "module_lists(['GeneA', 'GeneB', 'GeneC', 'GeneD'], [1,2,1,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get our flyatlas module lists\n",
    "modules_hierarchical = module_lists(tissue_expression.index, labels)\n",
    "\n",
    "# The sizes of our modules\n",
    "[len(module) for module in modules_hierarchical]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The modules themselves\n",
    "modules_hierarchical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4 - Clustering with k-medoids "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this task, we'll use a different clustering method and see what modules we get.\n",
    "\n",
    "We don't have a Euclidean distance space, but we have a pairwise matrix of distances, which allows us to implement k-medoids. Recall that in k-medoids, one of our data points acts as the centroid of each cluster - in this case there will be a centroid gene for each cluster. The algorithm we will implement for k-medoids is:\n",
    "\n",
    "1. Initialise the centroids randomly.\n",
    "2. Assign each gene to the closest centroid (in this case, using our network-based distances).\n",
    "3. Choose the most-central gene in each cluster to be the new centroid. This is the gene which minimises the sum of squared distances within the cluster.\n",
    "4. Repeat from step (2) until the algorithm converges and the centroids no longer change.\n",
    "\n",
    "The functions `initialise_centroids()` and `calculate_new_centroids()` are provided for you. You need to implement `assign_points()` and complete the function `kmedoids()` itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialise_centroids(N, k):\n",
    "    \"\"\"\n",
    "    Select k centroid indices randomly given that there are N data points.\n",
    "    We are only selecting the indices, so we don't need the actual data.\n",
    "    \"\"\"\n",
    "    centroid_indices = np.random.choice(list(range(N)), size=k, replace=False)\n",
    "    return centroid_indices "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_centroid_index(within_cluster_distances):\n",
    "    \"\"\"\n",
    "    Take in distance array of size (C,C) where C is the size of the cluster.\n",
    "    Return a centroid index (of type int) which is the most-central data point.\n",
    "    \"\"\"\n",
    "    # We'll minimise the sum of square distances\n",
    "    # Calculate this quantity for every point, and pick the best as the centroid\n",
    "    sse = np.sum((within_cluster_distances**2),axis=1)\n",
    "    return np.argmin(sse)\n",
    "\n",
    "def calculate_new_centroids(distances, assignments, k):\n",
    "    \"\"\"\n",
    "    Take distances of shape (N,N) and cluster assignments of shape (N) and\n",
    "    return centroid indices array of shape (k).\n",
    "    \"\"\"\n",
    "    assert np.max(assignments) < k\n",
    "    centroid_list = []\n",
    "    for c in range(k):\n",
    "        in_cluster = assignments==c\n",
    "        cluster_distances = distances[in_cluster,:][:,in_cluster]\n",
    "        centroid_index_in_cluster = calculate_centroid_index(cluster_distances)\n",
    "        centroid_index = np.array(list(range(len(distances))))[in_cluster][centroid_index_in_cluster]\n",
    "        centroid_list.append(centroid_index)\n",
    "    return np.array(centroid_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def assign_points(centroid_indices, distances):\n",
    "    \"\"\"\n",
    "    Assign each point to its closest centroid.\n",
    "    Take in an array of centroid indices of length k, and an array representing\n",
    "    the distance matrix, of shape (N,N).\n",
    "    Return a 1D array of length N representing cluster assignments.\n",
    "    Each value in the returned array should be a number from 0 to k-1,\n",
    "    indicating which cluster (centroid) this data point has been assigned to.\n",
    "    \"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~\n",
    "\n",
    "def kmedoids(distances, k):\n",
    "    \"\"\"\n",
    "    Implement k-medoids clustering on a given set of points, by taking in\n",
    "    a pre-computed distance matrix of size (N,N), and a number of clusters k.\n",
    "    Returns a tuple of (centroid_indices, cluster_assignments)\n",
    "    where \n",
    "    centroid_indices is a list of indices specifying which data points are now\n",
    "     centroids, and\n",
    "    cluster_assignments is a 1D array of length N, where the values of the\n",
    "     array are numbers from 0 to k-1 and represent cluster assignments.\n",
    "    \"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return array([0, 0, 1, 1])\n",
    "\n",
    "test_dist = np.array([[ 0.  ,  0.16,  0.89,  0.94],\n",
    "       [ 0.16,  0.  ,  0.83,  0.85],\n",
    "       [ 0.89,  0.83,  0.  ,  0.1 ],\n",
    "       [ 0.94,  0.85,  0.1 ,  0.  ]])\n",
    "assign_points([0,2], test_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should return array([0, 1, 1, 1])\n",
    "\n",
    "test_dist = np.array([[ 0.  ,  0.16,  0.89,  0.94],\n",
    "       [ 0.16,  0.  ,  0.83,  0.85],\n",
    "       [ 0.89,  0.83,  0.  ,  0.1 ],\n",
    "       [ 0.94,  0.85,  0.1 ,  0.  ]])\n",
    "assign_points([0,1], test_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# k-medoids is stochastic, so you're not guaranteed to get the correct result every \n",
    "# time for this test. However the most likely outcome is\n",
    "# (array([2, 0]), array([1, 1, 0, 0]))\n",
    "# or equivalently\n",
    "# (array([0, 2]), array([0, 0, 1, 1]))\n",
    "\n",
    "test_dist = np.array([[ 0.  ,  0.16,  0.89,  0.94],\n",
    "       [ 0.16,  0.  ,  0.83,  0.85],\n",
    "       [ 0.89,  0.83,  0.  ,  0.1 ],\n",
    "       [ 0.94,  0.85,  0.1 ,  0.  ]])\n",
    "kmedoids(test_dist, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can cluster our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cluster our flyatlas data\n",
    "centroids, kmedoids_labels = kmedoids(distances, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# genes that are the centroids\n",
    "tissue_expression.index[centroids]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the module_lists function you defined earlier\n",
    "modules_kmedoids = module_lists(tissue_expression.index, kmedoids_labels)\n",
    "\n",
    "# Module sizes\n",
    "[len(m) for m in modules_kmedoids]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The modules themselves\n",
    "modules_kmedoids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5 - Dimension Reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this task we will be performing Priciple Components Analysis to determine which gene in the first principle component has the highest contribution to the varience."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 5(a)**\n",
    "\n",
    "Perform a Principle Componanets Analysis on the log_expression matrix with the correct number of components and print the explained variance by component list (like we saw in the week 9 tutorial)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 5(b)**\n",
    "\n",
    "Print the gene that contributes most to the first eigen vector (the first principle component) of the PCA. The word limit on the second cell is 50 words."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - complete this cell ~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 6 - Centrality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this task we will be computing differing measures of centrality in the gene found above in Task 5. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 6 (a)**\n",
    "\n",
    "Degree Centrality\n",
    "\n",
    "Using the adjacency_binary matrix (with a threshold value of 0.85) claculate the degree centrality of the gene found in Task 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - your code here --"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 6 (b)**\n",
    "\n",
    "Betweeness Centrality\n",
    "\n",
    "Using the adjacency_binary matrix generated above(with a threshold value of 0.85), turn the matrix into a graph object (using networkx) and then claculate the betweeness centrality of the gene found in Task 5 using the [networkx function for betweeness centrality](https://networkx.github.io/documentation/networkx-2.2/reference/algorithms/generate/networkx.algorithms.centrality.betweenness_centrality.html#networkx.algorithms.centrality.betweenness_centrality)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - your code here --"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 6 (c)**\n",
    "\n",
    "Closeness Centrality\n",
    "\n",
    "Using the [dijkstra_path](https://networkx.github.io/documentation/networkx-1.10/reference/generated/networkx.algorithms.shortest_paths.weighted.dijkstra_path.html) function from networkx claculate the closeness centrality of the gene found in Task 5 using the following formula:\n",
    "\n",
    "\n",
    "$$C(u) = \\frac{n - 1}{\\sum_{v=1}^{n-1} d(v, u)}$$\n",
    "\n",
    "where d(v, u) is the shortest-path distance between v and u, and n is the number of nodes that can reach u."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~ GRADED CELL - your code here --"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 6 (d)**\n",
    "\n",
    "What does each measurement say about the gene's centrality? Is it relatively central? \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*# ~~ GRADED CELL - your answer here --*\n"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}