notebooks/10xResults.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "                                          10xResults v0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "with open(\"config.json\") as json_file:\n",
    "    parameter = json.load(json_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pickle\n",
    "\n",
    "\n",
    "with open(parameter[\"SAVE_DIR\"]+\"TCC_matrix.dat\", 'rb') as f:\n",
    "    T=pickle.load(f)\n",
    "with open(parameter[\"SAVE_DIR\"]+\"pwise_dist_L1.dat\", 'rb') as f:\n",
    "    D_l1=pickle.load(f)\n",
    "with open(parameter[\"SAVE_DIR\"]+\"nonzero_ec.dat\", 'rb') as f:\n",
    "    nonzero_ec=pickle.load(f)\n",
    "\n",
    "ecfile_dir = parameter[\"kallisto\"][\"TCC_output\"]+'matrix.ec'\n",
    "eclist=np.loadtxt(ecfile_dir,dtype=str)    \n",
    "\n",
    "from sklearn.preprocessing import normalize\n",
    "\n",
    "TCC=T.T\n",
    "T_norm = normalize(T, norm='l1', axis=0) \n",
    "T_normT = T_norm.transpose()\n",
    "    \n",
    "NUM_OF_CELLS=np.shape(T)[1]\n",
    "print(\"NUM_OF_CELLS =\", NUM_OF_CELLS\n)",
    "print(\"NUM_OF_nonzero_EC =\", np.shape(T)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "EC_dict = {}\n",
    "for i in range(np.shape(eclist)[0]):\n",
    "    EC_dict[i] = [int(x) for x in eclist[i,1].split(',')]\n",
    "    \n",
    "union=set()\n",
    "for i in nonzero_ec:\n",
    "    new = [tx for tx in EC_dict[i] if tx not in union] # filter out previously seen transcripts\n",
    "    union.update(new) \n",
    "union_list=list(union) #union of all transctipt ids seen in nonzero eq.classes\n",
    "NUM_OF_TX_inTCC = len(union)\n",
    "print(\"NUM_OF_Transcripts =\", NUM_OF_TX_inTCC #number of distinct transcripts in nonzero eq. classes )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#sort eq. classes based on size\n",
    "size_of_ec = [len(EC_dict[i]) for i in nonzero_ec]\n",
    "ec_idx = [i[0] for i in sorted(enumerate(size_of_ec), key=lambda x:x[1])]\n",
    "index_ec = np.array(ec_idx)\n",
    "\n",
    "ec_sort_map={}\n",
    "nonzero_ec_srt=[]; #init\n",
    "for i in range(len(nonzero_ec)):\n",
    "    nonzero_ec_srt+=[nonzero_ec[index_ec[i]]]\n",
    "    ec_sort_map[nonzero_ec[index_ec[i]]] = i\n",
    "nonzero_ec_srt = np.array(nonzero_ec_srt)\n",
    "\n",
    "ec_size_sort=np.array(size_of_ec)[index_ec]\n",
    "sumi = np.array(T.sum(axis=1))\n",
    "sumi_sorted=sumi[index_ec]\n",
    "total_num_of_umis = int(sumi_sorted.sum())\n",
    "total_num_of_umis_per_cell = np.array(T.sum(axis=0))[0,:]\n",
    "\n",
    "print(\"Total number of UMIs =\", total_num_of_umis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(sorted(total_num_of_umis_per_cell)[::-1], 'b-',linewidth=2.0)\n",
    "ax1.set_title('UMI counts per cell')\n",
    "ax1.set_xlabel('cells (sorted by UMI counts)')\n",
    "ax1.set_ylabel('UMI counts')\n",
    "ax1.set_yscale(\"log\", nonposy='clip')\n",
    "# ax1.set_xscale(\"log\", nonposy='clip')\n",
    "ax1.grid(True)\n",
    "ax1.grid(True,'minor')\n",
    "\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(sorted(sumi.reshape(np.shape(sumi)[0]))[::-1], 'r-',linewidth=2.0)\n",
    "ax1.set_title('UMI counts per eq. class')\n",
    "ax1.set_xlabel('ECs (sorted by UMI counts)')\n",
    "ax1.set_ylabel('UMI counts')\n",
    "ax1.set_yscale(\"log\", nonposy='clip')\n",
    "# ax1.set_xscale(\"log\", nonposy='clip')\n",
    "ax1.grid(True)\n",
    "ax1.grid(True,'minor')\n",
    "\n",
    "\n",
    "cell_nonzeros=np.array(((T_norm!=0)).sum(axis=0))[0]\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(total_num_of_umis_per_cell,cell_nonzeros, '.g',linewidth=2.0)\n",
    "ax1.set_title('UMI counts vs nonzero ECs')\n",
    "ax1.set_xlabel('total num of umis per cell')\n",
    "ax1.set_ylabel('total num of nonzero ecs per cell')\n",
    "ax1.set_yscale(\"log\", nonposy='clip')\n",
    "ax1.set_xscale(\"log\", nonposy='clip')\n",
    "ax1.grid(True)\n",
    "ax1.grid(True,'minor')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## TCC MEAN-VARIANCE \n",
    "def meanvar_plot(TCC_,alph=0.05):\n",
    "\n",
    "    TCC_var=np.var(TCC_.todense(),axis=0)\n",
    "    TCC_mean=np.mean(TCC_.todense(),axis=0)\n",
    "    TCC_mean=np.array(TCC_mean)[0]\n",
    "    TCC_var=np.array(TCC_var)[0]\n",
    "\n",
    "    fig = plt.figure()\n",
    "    N=TCC_.sum()\n",
    "    C=TCC_.shape[0]\n",
    "    ax = plt.gca()\n",
    "\n",
    "    ax.plot(TCC_mean ,TCC_var,'.', c='blue', alpha=alph, markeredgecolor='none')\n",
    "\n",
    "    xlims=[0.0001,10*TCC_mean.max()]\n",
    "    ax.set_xlim(xlims)\n",
    "    ax.set_ylim([0.0001,10*TCC_var.max()])\n",
    "    ax.set_yscale('log')\n",
    "    ax.set_xscale('log')\n",
    "    ax.plot(xlims, [(C-1)*(xlims[0])**2, (C-1)*(xlims[1])**2], color='g', linestyle='-', linewidth=2)\n",
    "    ax.plot(xlims, [(xlims[0]), (xlims[1])], color='k', linestyle='--', linewidth=1)\n",
    "    ax.set_title(\"TCC Mean-Variance [\"+str(TCC_.shape[1])+\" TCCs in \"+str(C)+\" Cells]\")\n",
    "    ax.set_xlabel(\"mean(TCC)\")\n",
    "    ax.set_ylabel(\"var(TCC)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "meanvar_plot(TCC,alph=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn import cluster,manifold\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def AffinityProp(D,pref,damp):\n",
    "    aff= cluster.AffinityPropagation(affinity='precomputed',\n",
    "                                     preference=pref,damping=damp, verbose=True)\n",
    "    labels=aff.fit_predict(D)\n",
    "    return labels\n",
    "\n",
    "def spectral(k,D):\n",
    "    spectral = cluster.SpectralClustering(n_clusters=k,affinity='precomputed')\n",
    "    spectral.fit(D)\n",
    "    labels = spectral.labels_\n",
    "    return labels\n",
    "\n",
    "def tSNE_pairwise(D):\n",
    "    tsne = manifold.TSNE(n_components=2, random_state=213, metric='precomputed', n_iter=2000, verbose=1);\n",
    "    X_tsne = tsne.fit_transform(D);\n",
    "    return X_tsne"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### t-SNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_tsne=tSNE_pairwise(D_l1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot function with (nonrandom) colors corresponding to labels (star a specific label)\n",
    "def stain_plot(X,labels,stain,title,nc=2,ax_lim=0,marksize=46):\n",
    "    \n",
    "    unique_labels = np.unique(labels)\n",
    "    N = len(unique_labels)\n",
    "    max_value = 16581375 #255**3\n",
    "    interval = int(max_value / N)\n",
    "    colors = [hex(I)[2:].zfill(6) for I in range(0, max_value, interval)]    \n",
    "    color= [(int(i[:2], 16)/float(255), int(i[2:4], 16)/float(255), \n",
    "             int(i[4:], 16)/float(255)) for i in colors]\n",
    "    i=0;\n",
    "    plt.figure(figsize=(15,10))\n",
    "    for label in unique_labels:\n",
    "        ind = np.squeeze(labels == label)\n",
    "        if label in stain: \n",
    "            plt.scatter(X[ind,0],X[ind,1],c='red',s=146,edgecolor='black',\n",
    "                        lw = 0.5, alpha=1,marker='*',label=label)\n",
    "        else:\n",
    "            plt.scatter(X[ind,0],X[ind,1],c=color[i],s=marksize,edgecolor='lightgray',\n",
    "                        lw = 0.5,label=label)        \n",
    "        i+=1   \n",
    "    plt.title(title)\n",
    "    plt.gray()\n",
    "    plt.legend(loc='upper right',bbox_to_anchor=(1.18, 1.01),ncol=nc)\n",
    "    if ax_lim>0:\n",
    "        plt.xlim([-ax_lim,ax_lim])  \n",
    "        plt.ylim([-ax_lim,ax_lim])  \n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot function with colors(heatmap) corresponding to the total number of umis per cell  (log10) \n",
    "def umi_counts_plot(X,total_num_of_umis_per_cell,title,ax_lim=0,,marksize=26):\n",
    "    \n",
    "    plt.figure(figsize=(15,10))\n",
    "    plt.scatter(X_tsne[:,0],X_tsne[:,1],c=np.log10(total_num_of_umis_per_cell),s=marksize,edgecolor='black',lw=0.25,cmap='OrRd')\n",
    "    plt.colorbar()\n",
    "    plt.title(title)\n",
    "    if ax_lim>0:\n",
    "        plt.xlim([-ax_lim,ax_lim])  \n",
    "        plt.ylim([-ax_lim,ax_lim])  \n",
    "    plt.axis('off')\n",
    "\n",
    "def color_plot(X,colorvec,title,ax_lim=0,marksize=26):\n",
    "    \n",
    "    plt.figure(figsize=(10,6))\n",
    "    plt.scatter(X[:,0],X[:,1],c=np.log10(1+colorvec),s=marksize,edgecolor='black',lw=0.25,cmap='OrRd')\n",
    "    plt.colorbar()\n",
    "    plt.title(title)\n",
    "    if ax_lim>0:\n",
    "        plt.xlim([-ax_lim,ax_lim])  \n",
    "        plt.ylim([-ax_lim,ax_lim])  \n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Spectral Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_of_clusters=2\n",
    "similarity_mat=D_l1.max()-D_l1\n",
    "labels_spectral = spectral(num_of_clusters,similarity_mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stain_plot(X_tsne,labels_spectral,[],\"TCC -- tSNE, spectral clustering\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Affinity Propagation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pref = -np.median(D_l1)*np.ones(NUM_OF_CELLS)\n",
    "labels_aff=AffinityProp(-D_l1,pref,0.5)\n",
    "np.unique(labels_aff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stain_plot(X_tsne,labels_aff,[],\"TCC -- tSNE, affinity propagation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=2)\n",
    "X_pca = pca.fit_transform(T_normT.todense())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stain_plot(X_pca,labels_aff,[],\"TCC -- PCA, affinity propagation\")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}