Python_Course_MeInBio_MATPLOTLIB-SEABORN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"IMAGES_PYTHON_COURSE_2018/LOGO_MEINBIO_TRANSPARENT.png\" width=\"200\" height=\"200\" />\n",
    "# <center>Python course 2018\n",
    "### <center>MeInBio Training Group \n",
    "#### <center> by Florian Heyl & Francesco Ferrari "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VISUALIZATION and PLOTTING\n",
    "\n",
    "Let's say it honestly. At the end of the day, what we really care about, is a cool plot to show to our boss (well, yes, data should also make sense, but hey ...)\n",
    "\n",
    "So, let's delve into the secret art of making plots with python. \n",
    "\n",
    "We will use two libraries for this task, MATPLOTLIB and SEABORN. The first is the most widely used library for visualizations, and it offers plenty of functions and endless flexibility. The second is built upon the first and offers easy functions that can give you really nice and elegant plots with very little coding effort. \n",
    "\n",
    "In this notebook we will learn:\n",
    "\n",
    "- how to make scatterplots and histograms\n",
    "- how to make barplots\n",
    "- how to make boxplots\n",
    "- how to make heatmaps\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### load all required modules\n",
    "\n",
    "import sys\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import random\n",
    "sys.path.append(os.path.abspath(\"./Book_and_other_resources/\"))\n",
    "import python_course_functions as pcf\n",
    "%matplotlib inline\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SCATTER PLOTS\n",
    "\n",
    "https://matplotlib.org/api/_as_gen/matplotlib.pyplot.scatter.html\n",
    "\n",
    "When you want to make a plot (or anything else), you should always refer to the official library documentation. There, you will find what kind os data structure is accepted by the function, and what are all the other parameters that you can tweek to adjust the appearence of your graph. Click on the link above, and scroll through the documentation of plt.scatter.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### easy scatterplot\n",
    "\n",
    "x=[1,2,3,4,5,6]\n",
    "y=[3.4,10,4.6,2,3,6.9]\n",
    "\n",
    "#here we set the paraments for how our image will look like\n",
    "fig = plt.figure(figsize=(4, 4), dpi= 150, facecolor='white', edgecolor='r')\n",
    "# dpi: dots per inch; # try to change facecolor to \"pink\", or \"blue\" and see what happens\n",
    "\n",
    "plt.scatter(x,y, label = \"# of reads per peak\") # here we plot our x and y coordinates \n",
    "                                                # and give them a label (for the legend)\n",
    "plt.title(\"My (first) scatterplot\\n\") # here we set the title\n",
    "plt.xlabel(\"x variable\") # here we choose the x-axis label\n",
    "plt.ylabel(\"y variable\") # here we choose the y-axis label\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig(\"./my_scatterplot.png\") # here we save the figure\n",
    "\n",
    "# if you check in your current folder, you will find a file with your plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise n° 1 (easy peasy)\n",
    "\n",
    "Bring the appropriate changes to the code below to achieve the following result:\n",
    "\n",
    "- the shape of the dots should be a triagle pointing upward;\n",
    "- the color of the triagles should be red;\n",
    "- the size of the triangles should be much bigger than the current dot size;\n",
    "- insert a legend\n",
    "\n",
    "Try to understand how to do that by reading the documentation  \n",
    "(Beware! Francesco and Florian are not being sloppy here. We want you to understand how to find your way through the documentation, so that you can learn to use any kind of plotting function!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### modify the code below\n",
    "\n",
    "x=[1,2,3,4,5,6]\n",
    "y=[3.4,10,4.6,2,3,6.9]\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(4, 4), dpi= 150, facecolor='w', edgecolor='r')\n",
    "\n",
    "plt.scatter(x,y) \n",
    "plt.title(\"My (first) scatterplot\\n\") \n",
    "plt.xlabel(\"x variable\") \n",
    "plt.ylabel(\"y variable\") \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:  \n",
    "<img src=\"IMAGES_PYTHON_COURSE_2018/my_scatterplot_mod.png\" width=\"600\" height=\"600\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HISTOGRAMS\n",
    "\n",
    "https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hist.html\n",
    "\n",
    "A histogram is a particularly useful tool that allows you to visualize the distribution of your data. Please, click on the link above and quickly browse through plt.hist documentation.\n",
    "\n",
    "Here under is how to make a very simple histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### easy histogram\n",
    "\n",
    "fig = plt.figure(figsize=(6, 4), dpi= 150, facecolor='w', edgecolor='r')\n",
    "\n",
    "# here we generate a random gaussian distribution with mean = 0 and sd = 2\n",
    "x = np.random.normal(loc = 0, scale = 2, size = 50000) \n",
    "\n",
    "# here we generate a random gaussian distribution with mean = 0 and sd = 1\n",
    "x_1 = np.random.normal(loc = 0, scale = 1, size = 50000)\n",
    "\n",
    "# here we generate a random gaussian distribution with mean = 0 and sd = 0.5\n",
    "x_2 = np.random.normal(loc = 0, scale = 0.5, size = 50000)\n",
    "\n",
    "# here we make the plots\n",
    "plt.hist(x, bins=100, label = r'Gauss. $\\mu$ = 0; $\\sigma$ = 2')\n",
    "plt.hist(x_1, bins=100, label = r'Gauss. $\\mu$ = 0; $\\sigma$ = 1')\n",
    "plt.hist(x_2, bins=100, label = r'Gauss. $\\mu$ = 0; $\\sigma$ = 0.5')\n",
    "\n",
    "plt.title(\"SOME RANDOM GAUSSIANS\\n\") # we set the title\n",
    "plt.xlabel(\"Value\") # we set the x-label\n",
    "plt.ylabel(\"Count\") # we set the y-label\n",
    "plt.legend() # we set the legend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise n° 2 (medium)\n",
    "\n",
    "We provide you with 4 lists containing the expression level of 10.000 genes from a RNA-Seq experiment. Gene expression level is expressed as log10(RPKM + 0.000001 pseudocount)\n",
    "\n",
    "Your task is  \n",
    "1) to plot the distribution of the normalized read counts and  \n",
    "2) to produce a scatter plot of gene expression level of ctr_1 VS ctr_2 and ctr_1 VS tr_1.\n",
    "\n",
    "If you are unsure, please refer to the documentation. If you get stuck, ask to the instructors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# expression level of ~ 10000 genes expressed in log10(RPKM + 0.000001 pseudocount) \n",
    "ctr_1 = np.log10(np.array(pcf.extract_column(\"Book_and_other_resources/Gene_Expression_table.tsv\", \"RPKM_Control_1\"))+0.000001)\n",
    "ctr_2 = np.log10(np.array(pcf.extract_column(\"Book_and_other_resources/Gene_Expression_table.tsv\", \"RPKM_Control_2\"))+0.000001)\n",
    "tr_1 = np.log10(np.array(pcf.extract_column(\"Book_and_other_resources/Gene_Expression_table.tsv\", \"RPKM_Treated_1\"))+0.000001)\n",
    "tr_2 = np.log10(np.array(pcf.extract_column(\"Book_and_other_resources/Gene_Expression_table.tsv\", \"RPKM_Treated_2\"))+0.000001)\n",
    "\n",
    "\n",
    "# 1) insert the appropriate code to plot the distribution of the read counts of ctr_1\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected Output:\n",
    "\n",
    "<img src=\"IMAGES_PYTHON_COURSE_2018/histogram_ctr_1.png\" width=\"600\" height=\"600\" />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2) insert the appropriate code to produce a scatterplot of the normalized read counts \n",
    "# of ctr_1 VS tr_1 AND of ctr_1 VS ctr_1\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected Output:\n",
    "<img src=\"./IMAGES_PYTHON_COURSE_2018/scatter_RNASEQ.png\" width=\"600\" height=\"600\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BARPLOTS\n",
    "\n",
    "https://matplotlib.org/api/_as_gen/matplotlib.pyplot.bar.html\n",
    "\n",
    "Very often we need to represent our data using a barplot, which represent the mean and the standard deviation (or the 95% confidence interval) of our data. \n",
    "\n",
    "In the example below, we want to plot the solid tumor size measured in mice treated with an immunotherapic medication vs ctr (data are invented).\n",
    "\n",
    "Have a look at the code below, which illustrates how you can produce nice, customizable barplots with matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### barplot with error bar\n",
    "\n",
    "# some made-up data representing tumor size in cm of three treatment groups (N=6)\n",
    "A_ctr = np.array([3, 2.3, 2.45, 2.89, 3.2, 2.4])\n",
    "A_tr_1 = np.array([1.9, 1.2, 1.4, 0.56, 1.1, 1])\n",
    "A_tr_2 = np.array([2.6, 1.34, 3.23, 2.9, 2.8, 2])\n",
    "\n",
    "# calculate the mean\n",
    "mean_ctr = np.mean(A_ctr)\n",
    "mean_tr_1 = np.mean(A_tr_1)\n",
    "mean_tr_2 = np.mean(A_tr_2)\n",
    "\n",
    "#calculate the standard deviation\n",
    "std_ctr = np.std(A_ctr)\n",
    "std_tr_1 = np.std(A_tr_1)\n",
    "std_tr_2 = np.std(A_tr_2)\n",
    "\n",
    "# Create Arrays for the plot\n",
    "tr_groups = ['Control', 'Treat 1', 'Treat 2']\n",
    "x_pos = np.arange(len(tr_groups))\n",
    "gr_means = [mean_ctr, mean_tr_1, mean_tr_2]\n",
    "error = [std_ctr, std_tr_1, std_tr_2]\n",
    "\n",
    "# Build the plot\n",
    "fig, ax = plt.subplots(figsize=(6, 4), dpi= 110)\n",
    "\n",
    "ax.bar(x_pos, gr_means, yerr=error, \n",
    "       align='center', alpha=0.5, ecolor='black', capsize=5, color=[\"blue\", \"red\",\"violet\"])\n",
    "ax.set_ylabel('Tumor size (cm)')\n",
    "ax.set_xticks(x_pos)\n",
    "ax.set_xticklabels(tr_groups)\n",
    "ax.set_title('Tumor size in three different treatment groups')\n",
    "ax.yaxis.grid(True)\n",
    "\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise n°3 (medium)\n",
    "\n",
    "You want to test the efficacy of a novel anti-tumor immunotherapic drug on your mouse system. You have been measuring how many days your mice survive after cancer onset in two treatment groups, a placebo and your drug of interest. You also want to see if the treatment is more or less active depending on how soon you start injecting the mice after cancer onset. So you further divide your mice into two treatment regimes, one with injections of drug/pacebo starting from 24h post cancer onset, and one with injections of drug/placebo starting from day 6 after cancer onset.\n",
    "\n",
    "Using the data provied, represent them using barplots. \n",
    "\n",
    "Hint: \n",
    "https://matplotlib.org/gallery/units/bar_unit_demo.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### measurements (days)\n",
    "random.seed(2)\n",
    "\n",
    "# here we generate our fake data\n",
    "placebo_24h = np.random.normal(11, 0.3, 10)\n",
    "drug_24h = np.random.normal(20,0.7, 10)\n",
    "\n",
    "placebo_6days = np.random.normal(9, 1, 10)\n",
    "drug_6days = np.random.normal(10, 1.9, 10)\n",
    "\n",
    "### insert your code in the space below\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected Outcome:  \n",
    "\n",
    "<img src=\"./IMAGES_PYTHON_COURSE_2018/placeboVSdrug.png\" width=\"600\" height=\"600\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BOXPLOT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://matplotlib.org/api/_as_gen/matplotlib.pyplot.boxplot.html\n",
    "\n",
    "A Boxplot is another usefull tool that allows you to visualize and compare the distribution of your data. \n",
    "\n",
    "To illustrate the code for building boxplots, we come back to the expression dataset and we try to visualize the distribution of our normalized read counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Visualize the distribution of normalized read counts using boxplots\n",
    "\n",
    "all_RPKM = [ctr_1, ctr_2, tr_1, tr_2]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 4), dpi= 130)\n",
    "\n",
    "boxplots = ax.boxplot(all_RPKM, patch_artist=True)\n",
    "ax.set_ylabel('Log10(RPKM + pseudocount)')\n",
    "ax.set_xticklabels(['CTR_1','CTR_2','TR_1','TR_2'])\n",
    "ax.set_title('Normalized Read Counts Distribution\\n')\n",
    "ax.yaxis.grid(True)\n",
    "\n",
    "for patch, color in zip(boxplots['boxes'], ['royalblue','dodgerblue','indianred','lightcoral']):\n",
    "    patch.set_facecolor(color)\n",
    "\n",
    "# or alternatively, change colors of the boxes one by one \n",
    "#boxplots['boxes'][0].set_facecolor('royalblue')  \n",
    "#boxplots['boxes'][1].set_facecolor('dodgerblue')  \n",
    "#boxplots['boxes'][2].set_facecolor('indianred')  \n",
    "#boxplots['boxes'][3].set_facecolor('lightcoral')  \n",
    "\n",
    "for line in boxplots['medians']:\n",
    "    line.set_color('k')\n",
    "\n",
    "# or alternatively, change colors of the median line one by one     \n",
    "#boxplots['medians'][0].set_color('k')\n",
    "#boxplots['medians'][1].set_color('k')\n",
    "#boxplots['medians'][2].set_color('k')\n",
    "#boxplots['medians'][3].set_color('k')\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "#print(boxplots)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise n°4 (medium)\n",
    "\n",
    "Now you want to visualize the distribution of your expression data based on the chromosomes in which genes map. \n",
    "\n",
    "When you are working with expression data, you usually import the normalized count table in the form of a dataframe using another module called **pandas** (unfortunately, we don't have time to introduce you to pandas, but we strongly reccomend you to read the documentation). Pandas dataframes offer quite a number of methods that allow for quick visualizations. \n",
    "\n",
    "Using the dataframe that we prepared for you, produce a boxplot of the normalized read counts distribution of Control_1 for each chromosome. Read the documentation of the following function:\n",
    "\n",
    "https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.boxplot.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### import data, add pseudocaunt and Log-transform \n",
    "all_RPKM_df = pd.read_csv(\"Book_and_other_resources/Gene_Expression_table.tsv\", sep=\"\\t\", index_col=0)\n",
    "all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]] = all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]] + 0.000001\n",
    "all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]] = np.log10(all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]])\n",
    "\n",
    "# insert your code in the space below (you can solve this exercise in 1 line)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "\n",
    "<img src=\"./IMAGES_PYTHON_COURSE_2018/BOXPLOT_CHR.png\" width=\"600\" height=\"600\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HEATMAPS\n",
    "\n",
    "https://seaborn.pydata.org/generated/seaborn.clustermap.html  \n",
    "https://seaborn.pydata.org/generated/seaborn.heatmap.html\n",
    "\n",
    "Heatmaps are fancy. They allow you to visualize numbers using a customizable color-code. \n",
    "\n",
    "Here, we show you how to make elegant heatmaps using some usefull seaborn functions (sns.clustermap and sns.heatmap). These functions not only allow you to plot heatmaps, but they perform clustering as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### create a clustered heatmap (hierarchical clustering) of our expression data\n",
    "\n",
    "### import data, add pseudocaunt and Log-transform \n",
    "all_RPKM_df = pd.read_csv(\"Book_and_other_resources/Gene_Expression_table.tsv\", sep=\"\\t\", index_col=0)\n",
    "all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]] = all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]] + 0.000001\n",
    "all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]] = np.log10(all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]])\n",
    "\n",
    "# select only relevant columns \n",
    "all_RPKM_df_1 = all_RPKM_df[[\"RPKM_Control_1\",\"RPKM_Control_2\", \"RPKM_Treated_1\",\"RPKM_Treated_2\"]]\n",
    "\n",
    "# plot a cluster heatmap of the data\n",
    "sns.set(font_scale = 2) # set labels size\n",
    "g = sns.clustermap(all_RPKM_df_1, yticklabels=False, figsize=(15, 15)) # do clustemap and set figure size\n",
    "ax = g.ax_heatmap  # get axis of figure\n",
    "ax.set_ylabel(\"GENES\")  # set y label \n",
    "\n",
    "\n",
    "# Does the result of the clustering make sense?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a normal heatmap (without clustering) \n",
    "\n",
    "# import data\n",
    "dat = pd.read_csv('Book_and_other_resources/heatmap_data.csv', sep=';', index_col=0)\n",
    "\n",
    "# select only interesting columns and clean dataframe\n",
    "dat1 = dat[list(dat)[:11]]\n",
    "dat1 = dat1.drop('Unnamed: 11', axis=1)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(18,9))\n",
    "ax = sns.heatmap(dat1, cmap='RdBu_r', vmin=0, vmax=200)\n",
    "\n",
    "ax.xaxis.set_ticks_position('top')\n",
    "plt.xticks(rotation=45, ha='left')\n",
    "\n",
    "#plt.savefig('heatmap.png', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:anaconda3]",
   "language": "python",
   "name": "conda-env-anaconda3-py"
  },
  "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
}