IntroToML_Eng_HighSchool.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Intro to ML, High school level - Practical session ###\n",
    "# Alice Schoenauer Sebag, Feb. 2018\n",
    "# Altschuler&Wu lab, UCSF\n",
    "\n",
    "## First of all, let us download the data. The original dataset [1] is publicly available \n",
    "# (see below), and it was pre-processed using the method described in [2]. \n",
    "# References: \n",
    "# [1] [Caie et al., Molecular Cancer Therapeutics, 2010], [Ljosa et al., Nature Methods, 2012]\n",
    "#      https://data.broadinstitute.org/bbbc/BBBC021/\n",
    "# [2] [Kang et al., Nature Biotech., 2016]\n",
    "#      https://www.nature.com/articles/nbt.3419\n",
    "\n",
    "import os, pandas\n",
    "import numpy as np\n",
    "if not os.path.isfile('IntroToML_TD.tar.gz'):\n",
    "    print('Be patient while the data downloads')\n",
    "    os.system('wget -c http://members.cbio.mines-paristech.fr/~aschoenauer/IntroToML_TD.tar.gz')\n",
    "    os.system('tar -xvzf IntroToML_TD.tar.gz')\n",
    "    print('Done!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## We have downloaded 2 files: the X (features) as well as the Y vector (drug classes, and some additional information).\n",
    "# Let us have a look at what is inside these files.\n",
    "Y = pandas.read_csv('IntroToML_TD/Caie_metadata.csv')\n",
    "\n",
    "# This will show how the file looks like: it contains what we want to predict, the DrugClass, and some additional info.\n",
    "Y.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Let's see how many datapoints we have:\n",
    "print('Nb of datapoints:',Y.shape[0])\n",
    "\n",
    "# **Question** Can you tell me how many drug classes are in the data?\n",
    "print('Nb of drug classes?') #your answer :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Let's do the same with the feature matrix\n",
    "X = pandas.read_csv('IntroToML_TD/Caie_features.csv')\n",
    "X.head()\n",
    "\n",
    "# **Question** How many features do we have?\n",
    "print('Nb of features?')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Now we would like to plot the data! In a nutshell, we're going to transform it first, \n",
    "# for the plot to be more informative.\n",
    "\n",
    "# ** Depending on your math level, read or skip this! **\n",
    "# Rather than have the axes be the features, we will have the axes be combinations of features, \n",
    "# the combinations along which the data varies most. This method is called PCA, for Principal \n",
    "# Component Analysis. Never forget: before applying PCA, we need to standardize the data because \n",
    "# PCA relies on Euclidean distances between datapoints.\n",
    "\n",
    "# To learn more about PCA, https://en.wikipedia.org/wiki/Principal_component_analysis\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "import matplotlib.pyplot as p\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "pca = PCA(n_components=3)\n",
    "# Data standardization\n",
    "nX = (X-np.mean(X,0))/np.std(X, 0)\n",
    "# Let us find the best projection and use it to transform the data.\n",
    "npX = pca.fit_transform(nX)\n",
    "\n",
    "# Now we can plot it.\n",
    "f=p.figure(figsize=(18,18))\n",
    "ax = f.gca(projection='3d')\n",
    "ax.scatter(npX[:,0], npX[:,1], npX[:,2])\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is all very well, but we would be interested to see how the drug classes are distributed in this space, not \n",
    "# only the datapoints. => Let us add colors\n",
    "\n",
    "# Picking some colors\n",
    "color_names = ['red', 'blue', 'orange', 'magenta', 'green','black', '0.5', '#bcbddc','yellow','cyan', '#a6cee3','#dd1c77', '#1f78b4', '#b2df8a', '#33a02c', '#fb9a99', '#e31a1c', '#fdbf6f', '#ff7f00', '#cab2d6', '#6a3d9a']\n",
    "\n",
    "# These are the distinct drug classes present in the data.\n",
    "distinct_drug_classes = sorted(list(Y.DrugClass.unique()))\n",
    "\n",
    "# Picking the right color for each datapoint.\n",
    "colors = []\n",
    "drug_classes = Y.DrugClass.values\n",
    "for drug_class in drug_classes:\n",
    "    index = distinct_drug_classes.index(drug_class)\n",
    "    colors.append(color_names[index])\n",
    "    \n",
    "# Let us plot again\n",
    "f=p.figure(figsize=(18,18))\n",
    "ax = f.gca(projection='3d')\n",
    "ax.scatter(npX[:,0], npX[:,1], npX[:,2], color=colors)\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# 3D plots can sometimes be misleading because of the perspective. Furthermore, what do you think of \n",
    "# plots without a legend? => Let us plot in 2D, add a legend and label the axes.\n",
    "\n",
    "f, subplots=p.subplots(3,1, figsize=(15,22))\n",
    "subplots[0].scatter(npX[:,0], npX[:,1], color=colors, alpha=0.5)\n",
    "subplots[0].set_xlabel('Principal component 1')\n",
    "subplots[0].set_ylabel('Principal component 2')\n",
    "subplots[1].scatter(npX[:,0], npX[:,2], color=colors, alpha=0.5)\n",
    "subplots[1].set_xlabel('Principal component 1')\n",
    "subplots[1].set_ylabel('Principal component 3')\n",
    "subplots[2].scatter(npX[:,1], npX[:,2], color=colors, alpha=0.5)\n",
    "subplots[2].set_xlabel('Principal component 2')\n",
    "subplots[2].set_ylabel('Principal component 3')\n",
    "# Adding the legend\n",
    "for i, drug_class in enumerate(distinct_drug_classes):\n",
    "    subplots[2].scatter(0,0, color=color_names[i], label=drug_class)\n",
    "subplots[2].legend()\n",
    "\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# **Question** If you look into what is in Y, you will see that you also have the information whether a datapoint is \n",
    "# a positive control, a negative control, or an actual experiment. Can you tell the difference between the three?\n",
    "# Better, can you redo the above plots, but with the colors corresponding to positive control, negative control,\n",
    "# or other?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Finally, the actual classification! For this, we will retrict ourselves to the subset of the data\n",
    "# whose drug class is actually known, and which show activity at the dose used.\n",
    "\n",
    "wh = np.where((Y.DrugClass=='Unknown')|(Y.Active==False))[0]\n",
    "# Deleting the useless lines in the data\n",
    "active_Y = Y.drop(wh, axis=0)\n",
    "active_X = np.delete(X.values, wh, 0)\n",
    "\n",
    "print('Nb of datapoints retained: ', active_X.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Now we are going to choose linear Support Vector Machines as our hypothesis space, ie, type of model.\n",
    "from sklearn.svm import SVC\n",
    "model = SVC(kernel='linear')\n",
    "\n",
    "# Let us fit the model to the data\n",
    "model.fit(active_X, active_Y.DrugClass)\n",
    "predicted_Y = model.predict(active_X)\n",
    "\n",
    "# Let us look at the error! We will compute the accuracy of the model and its confusion matrix. \n",
    "# The entry (i,j) of the confusion matrix will tell how many datapoints with true label i were predicted to be j.\n",
    "\n",
    "from sklearn.metrics import confusion_matrix, recall_score\n",
    "score = recall_score(active_Y.DrugClass, predicted_Y, average='micro')\n",
    "print(\"Model accuracy: \", score)\n",
    "conf_mat = confusion_matrix(active_Y.DrugClass, predicted_Y)\n",
    "print(conf_mat)\n",
    "\n",
    "# Pretty good! But... wait... what did we say about overfitting? Is this the right way to estimate the error\n",
    "# of our model on new data??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# **Question** Split the data into a training set and a test set as we discussed, and get a more accurate estimate\n",
    "# of the model's accuracy on data it has never seen."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}