LDA.ipynb

{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  },
  "name": "",
  "signature": "sha256:51c1a517adf6cf924f23bc24b134ead7ec06cfe132e01fc957a555442a8c4976"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Final Project - CS 281"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ankit Gupta"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Import Statements"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pandas as pd\n",
      "import numpy as np\n",
      "from scipy import stats as scistats\n",
      "import matplotlib.pyplot as plt\n",
      "import matplotlib.mlab as mlab\n",
      "from IPython.display import Image\n",
      "import lda\n",
      "from scipy.sparse import csr_matrix\n",
      "\n",
      "%matplotlib inline\n",
      "from CTM import CTM\n",
      "from CTMParallel import CTMParallel\n",
      "from sklearn.cross_validation import train_test_split"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pickle\n",
      "import sys"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import csv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Returns the min, max, and average accuracy across the folds\n",
      "def KFoldCrossValidation(X, Y, num_folds, print_outputs=True):\n",
      "    accuracies = []\n",
      "    for train, test in KFold(X.shape[0], num_folds, shuffle=True):\n",
      "        X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test]\n",
      "        logistic_classifier = LogisticRegression()\n",
      "        logistic_classifier.fit(X_train, Y_train)\n",
      "        accuracies.append(accuracy_score(logistic_classifier.predict(X_test), Y_test))\n",
      "    if print_outputs:\n",
      "        print \"Max Accuracy: \", max(accuracies), \"\\nMin Accuracy: \", min(accuracies), \"\\nAverage Accuracy: \", sum(accuracies)/len(accuracies)\n",
      "    return max(accuracies), min(accuracies), sum(accuracies)/len(accuracies)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 38
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "This part performs LDA on the Rat and Neuron Data, and saves the resultant models to pickle objects."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pandas as pd\n",
      "import numpy as np\n",
      "from scipy import stats as scistats\n",
      "import matplotlib.pyplot as plt\n",
      "import matplotlib.mlab as mlab\n",
      "from IPython.display import Image\n",
      "import lda\n",
      "from scipy.sparse import csr_matrix\n",
      "\n",
      "%matplotlib inline\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = pd.read_csv('counts_rat.csv').set_index('GeneID')\n",
      "genes_wanted = (data > 1).sum(axis=1) > 5\n",
      "genes_wanted = data.var(axis=1).divide(data.mean(axis=1)).sort(inplace=False, ascending=False)[:1000].index\n",
      "counts_newsetup = data.ix[genes_wanted, :].T\n",
      "classes = np.array(map(lambda x: x.split(\"_\")[0], counts_newsetup.index))\n",
      "vocab = counts_newsetup.columns\n",
      "counts = counts_newsetup.values\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "topic_options = [15, 25, 35, 45]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for topics in topic_options:\n",
      "    print \"Working on topic\", topics\n",
      "    sys.stdout.flush()\n",
      "    model = lda.LDA(n_topics=topics, n_iter=10)\n",
      "    model.fit(counts)\n",
      "    print \"    Saving pickled model\"\n",
      "    sys.stdout.flush()\n",
      "    pickle.dump(vocab, open(\"pickled_objects/lda_vocab_rat_\" + str(topics) + \"topics.p\", \"wb\"))\n",
      "    pickle.dump(counts, open(\"pickled_objects/lda_counts_rat_\" + str(topics) + \"topics.p\", \"wb\"))\n",
      "    pickle.dump(classes, open(\"pickled_objects/lda_classes_rat_\" + str(topics) + \"topics.p\", \"wb\"))\n",
      "    pickle.dump(model, open(\"pickled_objects/lda_model_rat_\" + str(topics) + \"topics.p\", \"wb\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 15\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 25\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 35\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 45\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "features =  pd.read_table(\"GSE60361_C1-3005-Expression.txt\", sep='\\t', usecols=range(3006)).set_index('cell_id').T\n",
      "lines = []\n",
      "with open(\"expression_mRNA_17-Aug-2014.txt\", 'r') as f:\n",
      "    reader = csv.reader(f, delimiter='\\t')\n",
      "    for line in reader:\n",
      "        lines.append(line)\n",
      "metadata = pd.DataFrame(lines)\n",
      "metadata = metadata.set_index(0)\n",
      "metadata.columns = metadata.ix['cell_id', :]\n",
      "classification = metadata.ix['group #', :]\n",
      "\n",
      "# This contains the classification for each gene into 9 major classes of genes.\n",
      "\n",
      "genes_of_interest = (-features.sum()).sort(inplace=False)[20:].index\n",
      "remove_low_expression_genes = features.sum() > 100\n",
      "genes_wanted = remove_low_expression_genes[remove_low_expression_genes].index.intersection(genes_of_interest)\n",
      "genes_wanted = features.T.var(axis=1).sort(inplace=False, ascending=False)[:1000].index.intersection(genes_wanted)\n",
      "print \"Number of genes being investigated\", len(genes_wanted)\n",
      "features = features.ix[:, genes_wanted]\n",
      "classification = classification[features.index]\n",
      "vocab = features.columns\n",
      "X_train, X_test, y_train, y_test, samples_train, samples_test = train_test_split(features, classification, features.index, train_size=500, test_size=100, random_state=42)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Number of genes being investigated 980\n"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "neuron_topic_options = [5, 15, 25, 35, 45, 55]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for topics in neuron_topic_options:\n",
      "    print \"Working on topic\", topics\n",
      "    sys.stdout.flush()\n",
      "    model = lda.LDA(n_topics=topics, n_iter=20)\n",
      "    model.fit(X_train.values)\n",
      "    print \"    Saving pickled model\"\n",
      "    sys.stdout.flush()\n",
      "    pickle.dump(vocab, open(\"pickled_objects/lda_vocab_neuron_\" + str(topics) + \"topics.p\", \"wb\"))\n",
      "    pickle.dump(X_train.values, open(\"pickled_objects/lda_counts_neuron_\" + str(topics) + \"topics.p\", \"wb\"))\n",
      "    pickle.dump(y_train.values, open(\"pickled_objects/lda_classes_neuron_\" + str(topics) + \"topics.p\", \"wb\"))\n",
      "    pickle.dump(model, open(\"pickled_objects/lda_model_neuron_\" + str(topics) + \"topics.p\", \"wb\"))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 5\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 15\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 25\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 35\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 45\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Working on topic 55\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "    Saving pickled model\n"
       ]
      }
     ],
     "prompt_number": 58
    }
   ],
   "metadata": {}
  }
 ]
}