patient_level_pca.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Dimensionality reduction on patient staining scores\n",
    "Each patient may have up to 10 stain scores associated with their cancer. We can do dimensionality reduction and clustering to group these patients. This script explores its use in Reverse the Odds data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn import cluster, decomposition\n",
    "import patient_level_lib as pll # use this functionality to load the data from local mongoDB database\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Settings \n",
    "set which measure to use for dimensionality reduction (look at patient_level_lib.combine_cores_per_patient)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "measure = 'aggregateSQSCorrected'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data\n",
    "into a pandas dataframe with a single column for each stain type\n",
    "\n",
    "df dims: nPatients * nStains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Cannot infer number of levels from empty list",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-483cd1400955>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mdf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpll\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread_patient_data_from_mongodb\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mD:\\GitHub\\CRUK\\rto_image_analysis\\patient_level_lib.py\u001b[0m in \u001b[0;36mread_patient_data_from_mongodb\u001b[1;34m(mongoFilter, projection)\u001b[0m\n\u001b[0;32m    114\u001b[0m     \u001b[1;31m# exclude patient_key and set up multi-index of (variable, stain)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    115\u001b[0m     \u001b[0mtups\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mliteral_eval\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 116\u001b[1;33m     \u001b[0mcols\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMultiIndex\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfrom_tuples\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtups\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    117\u001b[0m     \u001b[1;31m# set up new df with patient_key as index. You have to put into np array first, otherwise\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    118\u001b[0m     \u001b[1;31m# the existing columns and indices mess things up and you get nans only\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Anaconda2\\lib\\site-packages\\pandas\\core\\index.pyc\u001b[0m in \u001b[0;36mfrom_tuples\u001b[1;34m(cls, tuples, sortorder, names)\u001b[0m\n\u001b[0;32m   4785\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtuples\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   4786\u001b[0m             \u001b[1;31m# I think this is right? Not quite sure...\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 4787\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Cannot infer number of levels from empty list'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   4788\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   4789\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtuples\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mIndex\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mTypeError\u001b[0m: Cannot infer number of levels from empty list"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    "df = pll.read_patient_data_from_mongodb()[measure]\n",
    "print df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "exclude the test mre11 stain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'function' object has no attribute 'drop'",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-6-3abb9e015a75>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdrop\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'test mre11'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m: 'function' object has no attribute 'drop'"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    "df.drop('test mre11', axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualise missing data\n",
    "This would be patients for whom we did not measure a particular stain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = np.isnan(df)\n",
    "plt.pcolor(r)\n",
    "plt.set_xticklabels(df.columns)\n",
    "total_n = df.index.length\n",
    "print('Total n patients: ', total_n)\n",
    "for iStain in df.columns:\n",
    "    nMissing = np.sum(isnan(df[iStain]))\n",
    "    print(iStain, 'missing values: ', nMissing, ', completeness: ', 100*(total_n - nMissing)/total_n, '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot correlations between stains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-10-1b1d841faee5>, line 7)",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-10-1b1d841faee5>\"\u001b[1;36m, line \u001b[1;32m7\u001b[0m\n\u001b[1;33m    pd.tools.plotting.scatter_matrix(df, alpha=0.2, diagonal='kde'))\u001b[0m\n\u001b[1;37m                                                                   ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    "r = df.corr()\n",
    "plt.pcolor(r)\n",
    "plt.colorbar()\n",
    "plt.set_yticklabels(df.columns)\n",
    "plt.set_xticklabels(df.columns)\n",
    "plt.set_title('correlation between stains across patients, calculated pairwise')\n",
    "pd.tools.plotting.scatter_matrix(df, alpha=0.2, diagonal='kde'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster patients into groups (unsupervised)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'module' object has no attribute 'cluster'",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-15-ab0833b8c601>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mskl\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcluster\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mKMeans\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m: 'module' object has no attribute 'cluster'"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    "km = cluster.KMeans(n_clusters=3)\n",
    "km.fit_predict(df)\n",
    "print(km.cluster_centers_)\n",
    "print(km.labels_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal component analysis\n",
    "Will generate weights for each of the principal components for each patient. For example, one component might be high MRE11, low Ki67; then some patients will strongly express that component, other will not. Can then regress this on outcome variables. The more efficient way of doing that is probably PLS, but with PCA we also get a peek at what the principal components themselves are. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pca = decomposition.PCA()\n",
    "pca.fit(df)\n",
    "print(pca.components_.shape())\n",
    "print(pca.components_)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}