CLASSIC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Perform a k-fold cross-validation of linear SVC vs CADD**\n",
    "\n",
    "Configuration\n",
    "- n_cores: set to the number of CPU cores available on your machine; will spawn this number of parallel processes\n",
    "- k: number of folds to use in the cross-validation\n",
    "- confidence: (out of 1; NOT %) for the interval of the mean difference between linear SVC and CADD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_cores = 32\n",
    "k = 1000\n",
    "confidence = 0.999"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** No need to edit below here ** Just run each cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import csv, math, numpy as np, matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.metrics import roc_curve, roc_auc_score, auc, matthews_corrcoef\n",
    "from sklearn import cross_validation\n",
    "from multiprocessing import Pool\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the tab-separated values, and partition them according to the class.\n",
    "\n",
    "0 = non-pathogenic (ESP)\n",
    "\n",
    "1 = pathogenic (ClinVar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = [[], []]\n",
    "\n",
    "src = open(\"./all-scores.tsv\")\n",
    "for line in csv.reader(src, dialect=\"excel-tab\"):\n",
    "    y = int(line[0]);\n",
    "    map(float, line);\n",
    "    X[y].append(line[1:]);\n",
    "\n",
    "for i in range(2):\n",
    "    X[i] = np.array(X[i], dtype='float')\n",
    "\n",
    "target = np.array([0]*len(X[0]) + [1]*len(X[1])) # vector of 0s then 1s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate k folds with sklearn. This is done for both classes to ensure a balanced input. Training / test groups from each fold are paired (non-pathogenic / pathogenic) to be passed to a parallel pool for cross-validation.\n",
    "\n",
    "The ROC AUC is calculated using a linear SVC for combination of CADD, RVIS, and PROVEAN. This is compared to the CADD-only AUC for the same test-set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# select a subset of features to test independent contributions; 0=CADD, 1=RVIS, 2=PROVEAN\n",
    "featureCols = [0,1,2]\n",
    "\n",
    "# folds is a (2,1) list of 1 of the folds generated by cross_validation.KFolds on both X[0] and X[1]\n",
    "# these are paired by class to maintain a balance from each\n",
    "# foldType corresponds to 0=train and 1=test\n",
    "def getFoldFeatures(folds, foldType):\n",
    "    features = X[0][folds[0][foldType]];\n",
    "    features = np.append(features, X[1][folds[1][foldType]], axis=0)\n",
    "    return features;\n",
    "\n",
    "# as with getFoldFeatures, but returns a vector similar to targets, but shorter\n",
    "def getFoldTargets(folds, foldType):\n",
    "    return [0]*len(folds[0][foldType]) + [1]*len(folds[1][foldType]);\n",
    "\n",
    "# function to be mapped over all paired folds\n",
    "def trainAndValidate(folds):\n",
    "    X_train = getFoldFeatures(folds, 0);\n",
    "    y_train = getFoldTargets(folds, 0);\n",
    "    assert(len((X_train)==len(y_train)))\n",
    "    \n",
    "    X_test = getFoldFeatures(folds, 1);\n",
    "    y_test = getFoldTargets(folds, 1);\n",
    "    assert(len((X_test)==len(y_test)))\n",
    "\n",
    "    svc = SVC(kernel='linear').fit(X_train[:,featureCols], y_train);\n",
    "    decision = svc.decision_function(X_test[:,featureCols])\n",
    "    \n",
    "    # X_test[:,0] corresponds to CADD scores alone\n",
    "    return ({'y': y_test, 'decision': list(decision)}, roc_auc_score(y_test, decision), roc_auc_score(y_test, X_test[:,0]))\n",
    "    \n",
    "\n",
    "if __name__ == '__main__':\n",
    "    p = Pool(n_cores)\n",
    "\n",
    "    kf = list(); # folds for each class in X; shape=(2,1)\n",
    "    for i in range(2):\n",
    "        kf.append(list(cross_validation.KFold(n=len(X[i]), n_folds=k, shuffle=True, random_state=314159))); # constant random_state allows for reproduction of results; 314159 (pi * 10^5) is a nothing-up-my-sleeve number https://en.wikipedia.org/wiki/Nothing_up_my_sleeve_number\n",
    "    \n",
    "    pairedFolds = list(); # list of tuples containing fold for each class\n",
    "    for f in range(k):\n",
    "        pairedFolds.append((kf[0][f], kf[1][f]));\n",
    "    \n",
    "    results = p.map(trainAndValidate, pairedFolds);\n",
    "    pairedAUC = [];\n",
    "    for r in results:\n",
    "        pairedAUC.append([r[1], r[2]])\n",
    "        \n",
    "    pairedAUC = np.reshape(np.array(pairedAUC, dtype='float'), (k, 2));\n",
    "    p.close()\n",
    "    \n",
    "    CADD = pairedAUC[:,1]\n",
    "    CLASSIC = pairedAUC[:,0]\n",
    "    improvements = CLASSIC - CADD;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate statistics of a paired t-test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean AUC CADD:    0.92787\n",
      "mean AUC CLASSIC: 0.95185\n",
      "n = 1000\n",
      "mean diff: 0.02399\n",
      "t-statistic: 17.05394\n",
      "df: 999\n",
      "p-value: 1.9717e-57 (two-sided)\n",
      "99.9% CI: 0.01934 0.02863\n"
     ]
    }
   ],
   "source": [
    "print('mean AUC CADD:    %0.5f' % CADD.mean())\n",
    "print('mean AUC CLASSIC: %0.5f' % CLASSIC.mean())\n",
    "print('n = %u' % len(improvements))\n",
    "print('mean diff: %0.5f' % (CLASSIC.mean() - CADD.mean()))\n",
    "\n",
    "t = stats.ttest_rel(CLASSIC, CADD)\n",
    "df = len(improvements) - 1;\n",
    "sem = stats.sem(improvements)\n",
    "\n",
    "print('t-statistic: %0.5f' % t[0])\n",
    "print('df: %u' % df)\n",
    "print('p-value: %0.5g (two-sided)' % t[1])\n",
    "\n",
    "critical = sem * stats.t._ppf((1+confidence)/2, df)\n",
    "print('%0.1f%% CI: %0.5f %0.5f' % (confidence*100, improvements.mean() - critical, improvements.mean() + critical))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the Matthews Correlation Coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MCC: 0.78\n"
     ]
    }
   ],
   "source": [
    "y_decision = []\n",
    "decision = []\n",
    "\n",
    "for r in results:\n",
    "    y_decision += r[0]['y']\n",
    "    decision += r[0]['decision']\n",
    "\n",
    "assert(len(y_decision)==len(decision))\n",
    "    \n",
    "y = []\n",
    "for d in decision:\n",
    "    y.append(1 if d>0 else 0)\n",
    "    \n",
    "print('MCC: %0.2f' % matthews_corrcoef(y_decision, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Produce ROC curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rocCLASSIC = roc_curve(y_decision, decision)\n",
    "rocCADD = roc_curve(target, list(X[0][:,0]) + list(X[1][:,0]))\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "plt.plot(rocCLASSIC[0], rocCLASSIC[1], 'k', label='CLASSIC (%0.1f%%)' % (auc(rocCLASSIC[0], rocCLASSIC[1])*100))\n",
    "plt.plot(rocCADD[0], rocCADD[1], 'k--', label='CADD (%0.1f%%)' % (auc(rocCADD[0], rocCADD[1])*100))\n",
    "\n",
    "plt.plot([0, 1], [0, 1], 'k:')\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.0])\n",
    "\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.legend(loc=\"lower right\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a master model from all folds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "svc = SVC(kernel='linear').fit(np.append(X[0], X[1], axis=0), target);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.55141463 -0.44124492 -0.2671904 ]]\n",
      "[-1.94571548]\n"
     ]
    }
   ],
   "source": [
    "print(svc.coef_) # CADD / RVIS / PROVEAN\n",
    "print(svc.intercept_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Output some meta data for inclusion in the manuscript"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 14435; 88.37% of all CADD-tested sites\n",
      "neutral n=6968\n",
      "pathogenic n=7467\n",
      "0.5172843782473155\n",
      "3.4.0 (default, Jun 19 2015, 14:20:21) \n",
      "[GCC 4.8.2]\n",
      "0.16.1\n"
     ]
    }
   ],
   "source": [
    "print('n = %u; %0.2f%% of all CADD-tested sites' % (len(target), len(target)/16334*100))\n",
    "print('neutral n=%u' % len(X[0]))\n",
    "print('pathogenic n=%u' % len(X[1]))\n",
    "print(max(len(X[0]), len(X[1]))/len(target))\n",
    "\n",
    "import sys, sklearn\n",
    "print(sys.version)\n",
    "print(sklearn.__version__)"
   ]
  }
 ],
 "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.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}