--- a +++ b/.ipynb_checkpoints/scGEM_SCOT_alignment-checkpoint.ipynb @@ -0,0 +1,197 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Notebook for running SCOT on scGEM Data\n", + "**Note:** This version of the notebook runs a new setting for SCOT, where we use correlation as a metric for building kNN graphs and use connectivity information from this graph in intra-domain similarity matrices fed into the optimal transport algorithm. \n", + "\n", + "SCOT software has been updated on 20 September 2020. It now outputs error statements for convergence issues at low epsilon values. When it runs into numerical instabilities in convergence, it outputs None, None instead of X_new, y_new. If you run into such an error, please try using a larger epsilon value for the entropic regularization. \n", + "\n", + "If you have any questions, e-mail: ritambhara@brown.edu, pinar_demetci@brown.edu, rebecca_santorella@brown.edu " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "import src.utils as ut\n", + "import src.evals as evals\n", + "from src.scot import *" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Dimensions of input datasets are: X= (177, 27) y= (177, 34)\n" + ] + } + ], + "source": [ + "X=np.genfromtxt(\"data/scGEM_methylation.csv\", delimiter=\",\")\n", + "y=np.genfromtxt(\"data/scGEM_expression.csv\", delimiter=\",\")\n", + "print(\"Dimensions of input datasets are: \", \"X= \", X.shape, \" y= \", y.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "It. |Err \n", + "-------------------\n", + " 0|4.960354e-03|\n", + " 10|3.044009e-04|\n", + " 20|1.238559e-04|\n", + " 30|3.991449e-05|\n", + " 40|1.207972e-05|\n", + " 50|3.565486e-06|\n", + " 60|1.043156e-06|\n", + " 70|3.043384e-07|\n", + " 80|8.871365e-08|\n", + " 90|2.585305e-08|\n", + " 100|7.533550e-09|\n", + " 110|2.195218e-09|\n", + " 120|6.396653e-10|\n" + ] + } + ], + "source": [ + "X=ut.zscore_standardize(X)\n", + "y=ut.zscore_standardize(y)\n", + "X_new,y_new= scot(X, y, k=35, e=5e-3, mode=\"connectivity\", metric=\"correlation\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Evaluating Alignment" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Average FOSCTTM score for this alignment is: 0.216454802259887\n" + ] + } + ], + "source": [ + "fracs=evals.calc_domainAveraged_FOSCTTM(X_new, y_new)\n", + "print(\"Average FOSCTTM score for this alignment is: \", np.mean(fracs))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hU1dbA4d+idwRDr5EWUyCEIiodkQgqoDQbYLsWxPJdG3rFhlcpXtu1oSIqCIhyQQUFEUUUQRKKFAFpShDpSugkWd8fezIkIRkCZEqS9T7PPDOnzDlrDmHWnH32WVtUFWOMMYVXkWAHYIwxJrgsERhjTCFnicAYYwo5SwTGGFPIWSIwxphCrliwAzhdYWFhWr9+/WCHYYwx+UpiYuJuVa2S3bJ8lwjq169PQkJCsMMwxph8RUR+y2mZNQ0ZY0whZ4nAGGMKOUsExhhTyOW7awTZOX78OElJSRw5ciTYoRhTKJQqVYratWtTvHjxYIdi8kCBSARJSUmUL1+e+vXrIyLBDseYAk1V2bNnD0lJSYSHhwc7HJMHCkTT0JEjRzj33HMtCRgTACLCueeea2fgBUiBSASAJQFjAsj+vxUsBSYRGGOMOTOWCPLIM888Q1RUFE2bNiU2NpbFixcD7kL2ww8/TKNGjYiOjqZ169Z88cUXAPz9998MHDiQBg0a0KBBAwYOHMjff//NypUriY2NJTY2lsqVKxMeHk5sbCyXXHJJruOpX78+u3fvBuCiiy7K+w98lpYvX86sWbOyXfbtt99SsWJF7zHI+LnHjh1LREQEERERtG7dmu+//9677PPPP6d58+Y0a9aMyMhI3nzzTe+y999/n+joaKKiooiMjGTMmDEMGTKE2NhYIiMjKV26tHd/IkJsbCwNGzbMFMfChQvp2LEjdevWJeM4Hr169aJcuXJ+OErGBIiq5qtHixYtNKs1a9acNC+QFi5cqG3atNEjR46oququXbt027Ztqqr60EMP6cCBA73L/vzzT50yZYqqql599dX6+OOPe7czfPhw7dOnT6ZtDxo0SKdOnXraMdWrV0937dp1Jh8nIN59910dMmRItsu++eYb7dGjx0nzP/vsM42Li/N+rsTERK1Tp45u375djx07pjVq1NCtW7eqquqRI0d07dq1qqo6a9Ysbd68ufff5PDhwzp27Fjvdjdv3qxRUVG5iqNDhw4aExOjCxYsUFXVffv2aevWrbVs2bKnewjyvWD/vytUkpNVb7xR9dtvz3gTQILm8L1qZwR5YPv27YSFhVGyZEkAwsLCqFmzJocOHeKtt97ilVde8S6rVq0a/fr1Y8OGDSQmJvLYY495tzN8+HASEhLYuHFjrvfdq1cvWrRoQVRUFGPHjs12nfRfq2lpadx5551ERUVx+eWX0717dz7++GPAnUE8/vjjxMXFERMTw9q1awF44oknGDRoEJdeein169dn2rRpPPjgg8TExBAfH8/x48cBSExMpEOHDrRo0YJu3bqxfft2ADp27MhDDz1E69atady4MQsWLODYsWMMHz6cKVOmEBsby5QpU3L1WUeOHMno0aMJCwsDIC4ujkGDBvHqq6+SnJxMSkoK5557LgAlS5akSZMmADz77LOMGTOGmjVrAq7r46233prrY5zVgAEDmDx5MgDTpk3jqquuOuNtGZMrO3fCu+/C5s1+2bwlgjxw6aWXsnXrVho3bsydd97J/PnzAdiwYQN169alQoUKJ71nzZo1xMbGUrRoUe+8okWLEhsby+rVq3O973HjxpGYmEhCQgIvv/wye/bsyXHdadOmsWXLFlauXMnbb7/Njz/+mGl5WFgYS5cu5Y477mDMmDHe+Rs3bmTmzJnMmDGD66+/nk6dOrFy5UpKly7NzJkzOX78OEOHDuXjjz8mMTGRm266iUcffdT7/pSUFH766SdefPFFnnzySUqUKMFTTz1F//79Wb58Of379z8p1gULFnibZJ555hkAVq9eTYsWLTKt17JlS1avXk3lypW58sorqVevHtdccw0TJ04kLS0NgFWrVp30vrPRpUsXvvvuO1JTU5k8eXK28RuTp/budc+eHzp5rUDcRxBs5cqVIzExkQULFvDNN9/Qv39/nnvuOeLi4nJ8j6pm2/Mip/k5efnll/nf//4HwNatW/n111+9v4qz+v777+nbty9FihShevXqdOrUKdPy9F+2LVq0YNq0ad75l112GcWLFycmJobU1FTi4+MBiImJYcuWLaxbt45Vq1bRtWtXAFJTU6lRo0a2292yZUuuPle7du34/PPPT7lexuP19ttvs3LlSubOncuYMWP46quvGD9+fK72dzqKFi1K27ZtmTJlCocPH8aq4Rq/S08ElSv7ZfOWCPJI0aJF6dixIx07diQmJob33nuPfv368fvvv5OcnEz58uUzrR8VFcWyZctIS0ujSBF3YpaWlsaKFSs4//zzc7XPb7/9lrlz5/Ljjz9SpkwZOnbs6LNvt2a4wJmd9OarokWLkpKSctL8IkWKULx4ce8Xb5EiRUhJSUFViYqKOukM41TbPV2RkZEkJibSuXNn77ylS5cSGRnpnY6JiSEmJoYbbriB8PBwxo8fT1RU1EnvO1sDBgygd+/ePPHEE3m2TWNy5OdEYE1DeWDdunX8+uuv3unly5dTr149ypQpw80338zdd9/NsWPHAHc9YcKECTRs2JDmzZszYsQI7/tGjBhBXFwcDRs2zNV+//77bypVqkSZMmVYu3YtixYt8rl+27Zt+eSTT0hLS2PHjh18++23p/9hs9GkSRN27drlTQTHjx8/ZfNW+fLlSU5OPq39PPjggzz00EPe5q/ly5czfvx47rzzTg4cOJDp86T/GwAMGzaMBx98kD///BOAo0eP8vLLL5/WvrNq164dw4YN45prrjmr7RiTK0ePQrly+fOMQETigZeAosDbqvpcluUVgQlAXU8sY1T1XX/G5A8HDhxg6NCh/PXXXxQrVoyGDRt6L9yOGDGCf/3rX0RGRlKqVCnKli3LU089BcA777zD0KFDadiwIarKhRdeyDvvvJPr/cbHx/PGG2/QtGlTmjRpQps2bXyuf/XVV/P1118THR1N48aNueCCC6hYseKZf3CPEiVK8PHHH3P33Xfz999/k5KSwr333ktUVFSO7+nUqRPPPfccsbGxDBs2LFft7FdeeSXbtm3joosuQkQoX748EyZMoEaNGiQnJzNq1Chuu+02SpcuTdmyZb3NQt27d2fHjh1ccskl3qakm2666aw+s4hw//33n9U2jMm1QYPcw0/kVM0FZ7xhkaLAeqArkAQsAa5R1TUZ1nkEqKiqD4lIFWAdUF1Vj+W03ZYtW2rWgWl++eWXXDenFHYHDhygXLly7Nmzh9atW/PDDz9QvXr1YIdl8iH7f5e/iEiiqrbMbpk/zwhaAxtUdZMniMlAT2BNhnUUKC+u0bkcsBc480Zkc0qXX345f/31F8eOHeOxxx6zJGBMfjBmDGzbBi+84JfN+zMR1AK2ZphOAi7Iss5/gU+BP4DyQH9VTcu6IRH5B/APgLp16556zx07njyvXz+48044dAi6dz95+eDB7rF7N/Tpk3lZHrWlh4K8ui5gjAmgb78Fz705/uDPi8XZ9YHM2g7VDVgO1ARigf+KyEmd7lV1rKq2VNWWVapkO/ayMcYUXHv3+u1CMeC/EhPAhcDsDNPDgGFZ1pkJtMswPQ9o7Wu7oVhiorDxV/mK8ePHa8OGDbVhw4Y6fvz4bNd5/vnn9fzzz9eYmBjt3LmzbtmyxbusW7duWrFixZPKQqSlpekjjzyijRo10oiICH3ppZcyLf/pp5+0SJEimUp51KtXT6Ojo7VZs2aa3d/c6NGjFfAeh2PHjunAgQM1OjpaIyIi9N///rd33YSEBI2OjtYGDRro0KFDNS0tzednWbZsmbZp00YjIyM1JiZGJ0+e7N3WoEGDtH79+tqsWTNt1qyZLlu2TFVV9+7dq7169dKYmBht1aqVrly50vue//znPxoZGalRUVE6YMAAPXz4sHc/F1xwgfczLl68WFVVJ0yY4N1+s2bNVES8+8nI/t8FUJMmqv36ndUm8FFiwp+JoBiwCQgHSgArgKgs67wOPOF5XQ3YBoT52m5hTQQpKSnBDsHLH4lgz549Gh4ernv27NG9e/dqeHi47t2796T15s2bpwcPHlRV1ddee037ZfjPMXfuXP30009PSgTjxo3TG264QVNTU1VVdceOHd5lKSkp2qlTJ73ssstOSgQ5fcbff/9dL730Uq1bt653nYkTJ2r//v1VVfXgwYNar1493bx5s6qqtmrVShcuXKhpaWkaHx+vs2bN8vlZ1q1bp+vXr1dV1W3btmn16tV13759qppz7an7779fn3jiCVVV/eWXX7Rz586qqpqUlKT169fXQ4cOqapq37599d1331VV1a5du3pjmTlzpnbo0OGk7f78888aHh6e7XEoDP/vQkZYmOodd5zVJnwlAr81DalqCnAXMBv4BfhIVVeLyO0icrtntaeBi0RkJfA18JCq7vZXTIFyxx130LJlS6Kionj88ccB+OKLL+jXr593nW+//ZYrrrgCgDlz5nDhhRcSFxdH3759OXDgAODq/zz11FO0bduWqVOn8tZbb9GqVSuaNWvG1VdfzaFDhwBXAqJNmza0atWK4cOHZ6qEOXr0aFq1akXTpk29sWT0+uuv8+CDD3qnx48fz9ChQ4FT1zHasmUL0dHR3ukxY8Z4b7DauHEj8fHxtGjRgnbt2nlrF+Vk9uzZdO3alcqVK1OpUiW6du3Kl19+edJ6nTp1okyZMgC0adOGpKQk77IuXbqcdONe+mccPny498a9qlWrepe98sorXH311Znmncp9993HqFGjMt0BLiIcPHiQlJQUDh8+TIkSJahQoQLbt29n//79XHjhhYgIAwcOZPr06T4/S+PGjWnUqBEANWvWpGrVquzatctnTGvWrKFLly4AREREsGXLFnbs2AHgjSklJYVDhw55ay6JCPv37wfcPSnp8zOaNGmS3SsRbKpQtSrk5vrome8j+BVFT+eRH84I9uzZo6ru12aHDh10xYoVevz4ca1Tp44eOHBAVVVvv/12/eCDD3TXrl3arl077/znnntOn3zySVV1v0pHjhzp3e7u3bu9rx999FF9+eWXVVW1R48e+uGHH6qq6uuvv+6thDl79my99dZbNS0tTVNTU7VHjx46f/78TLHu3LlTGzRo4J2Oj4/3VtZM/xyHDh3SqKgo7/7Tfy1nrdo5evRobzXVzp07e3/VLlq0SDt16qSqqjNmzNDHHnvspGM2evRoffrpp73TTz31lI4ePdrncR4yZEim96hmXzG0cuXKOmLECG3RooXGx8d740pKStL27dtrSkrKSb+069evr82bN9e4uDh98803vfNnzJihd999d6bjoOqahvr3769hYWFapkwZ73uWLFmiXbp08b7/u+++y7ayanafRVV18eLFGhER4T2bGTRokDZu3FhjYmL03nvv9Va1HTZsmN53333e9xQtWlQTEhJUVfXFF1/UsmXLalhYmF577bXeba9Zs0br1KmjtWvX1po1a2ZqZkt33nnnZWpmyijU/t8Z37Dqo4H10UcfERcXR/PmzVm9ejVr1qyhWLFixMfH89lnn5GSksLMmTPp2bMnixYtYs2aNVx88cXExsby3nvv8dtvv3m3lfFGq1WrVtGuXTtiYmKYOHGi9+7dH3/8kb59+wJw7bXXetefM2cOc+bMoXnz5sTFxbF27dpMd0ADVKlShfPOO49FixaxZ88e1q1bx8UXXwy4OkbNmjWjTZs23jpGuXHgwAEWLlxI3759iY2N5bbbbvNWI73yyiu9N9Rl5P5OM/NVc2nChAkkJCTwwAMPnDKeo0ePUqpUKRISErj11lu9N5Pde++9jBw5MlPhv3Q//PADS5cu5YsvvuDVV1/lu+++49ChQzzzzDPZxv/TTz9RtGhR/vjjDzZv3szzzz/Ppk2bcvW5cvos27dv54YbbuDdd9/1ns08++yzrF27liVLlrB3715GjhwJwMMPP8y+ffuIjY3llVdeoXnz5hQrVox9+/YxY8YMNm/ezB9//MHBgweZMGEC4M6UXnjhBbZu3coLL7zAzTffnGn/ixcvpkyZMpnO+kzBZLWG8tjmzZsZM2YMS5YsoVKlSgwePNhb/6d///68+uqrVK5cmVatWlG+fHlUla5duzJp0qRst1e2bFnv68GDBzN9+nSaNWvG+PHjT9kVVFUZNmwYt912m8/1+vfvz0cffURERAS9e/dGRHJVx6hYsWLeCp+Ad3laWhrnnHMOy5cv97nfjGrXrp3p8yQlJdExu27AwNy5c3nmmWeYP3++t47RqbZ99dVXA9C7d29uvPFGABISEhgwYAAAu3fvZtasWRQrVoxevXp5m0mqVq1K7969+emnn6hUqRKbN2+mWbNm3hjj4uL46aef+PDDD4mPj6d48eJUrVqViy++mISEBNq1a5ep+SopKSlTE0xOn2X//v306NGDESNGZLpjPL2YX8mSJbnxxhu9VWIrVKjAu++6m/JVlfDwcMLDw5k9ezbh4eGk97a76qqrWLhwIddffz3vvfceL730EgB9+/bllltuyXTcJk+ebM1CoWD1arjrLhg1Clq18ssu7Iwgj+3fv5+yZctSsWJFduzY4R2NDFxt/qVLl/LWW295f+m3adOGH374gQ0bNgBw6NAh1q9fn+22k5OTqVGjBsePH2fixIne+W3atOGTTz4B8NbJB+jWrRvjxo3zXnPYtm0bO3fuPGm7V111FdOnT2fSpEneuHJTx6hatWrs3LmTPXv2cPToUW+10AoVKhAeHs7UqVMB98W0YsUKn8etW7duzJkzh3379rFv3z7mzJlDt27dTlpv2bJl3HbbbXz66ae5btfv1asX8+bNA2D+/Pk0btwYcEl7y5YtbNmyhT59+vDaa6/Rq1cvDh486K2DdPDgQebMmUN0dDQxMTHs3LnT+57atWuzdOlSqlevTt26dZk3bx6qysGDB1m0aBERERHUqFGD8uXLs2jRIlSV999/n549e/r8LMeOHaN3794MHDjQe6aXLv3MSlWZPn2699d6+k2C4Kqwtm/fngoVKlC3bl0WLVrEoUOHUFW+/vpr793ANWvW9JZMnzdvnve6BLhkPnXqVG+iNEG0bZu7j+BYjgUXzl5ObUah+sgP1wgGDRqkERER2r17d+3du7e3l4aqawsuW7ast7eIqurXX3+tLVu21JiYGI2JidEZM2ao6sk9V1577TWtX7++dujQQe+66y4dNGiQqqquX79eW7dura1atdInnnhCa9as6X3Piy++qNHR0RodHa1t2rTRDRs2ZBtzjx49MvUOOXLkiMbHx2tMTIz26dNHO3TooN98881Jcb300kvaoEEDveSSS3TQoEHeawSbNm3Sbt26adOmTfX888/3XvfI6RqBquo777yjDRo00AYNGui4ceO88x977DHvMenSpYtWrVrV27Xxiiuu8K7Xtm1bDQsL01KlSmmtWrX0yy+/VFU3ilj37t29x2D58uUn7TvjNYKNGzdq06ZNtWnTphoZGakjRozINt6MxyE5OVn79OmjkZGRev755+uoUaO86y1ZskSjoqL0vPPO0yFDhni7j+b0WT744AMtVqxYpi6c6d03O3XqpNHR0RoVFaXXXXedJicnq6obJa9hw4bapEkT7d27d6YeV8OHD9cmTZpoVFSUXn/99d7rCgsWLNC4uDht2rSptm7d2ntNQdVda7nggguy/dzpQu3/XYE1aZIqqJ7l8cbHNQK/1RryF6s1dLJDhw5RunRpRITJkyczadIkZsyYEeywTAFX2P/fBcxrr8GQIfDnn1Ct2hlvJli1hkyAJCYmctddd6GqnHPOOYwbNy7YIRlj8kr6WASVKvltF5YICoB27dqdsg3eGJNPVaoEF1wAJUr4bRcF5mJxfmviMiY/s/9vATRkCJxi0KmzVSASQalSpdizZ4/9cRoTAKrKnj17KFWqVLBDMXmkQDQN1a5dm6SkpFPehm+MyRulSpWidu3awQ6jcOjcGdq3Bz+Oj10gEkHx4sUJDw8PdhjGGJO3UlLg+++hdWu/7qZANA0ZY0yB9PvvcPw4ZLjZzx8sERhjTKjyVBywRGCMMYVVeqHHhg39uhtLBMYYE6qqVYMePcBTbNBfLBEYY0yo6tMHPv8cfJRkzwuWCIwxJlSlpgZkN5YIjDEmFKWmQsWK4Bl8yJ8sERhjTCjauhUOHoRzz/X7riwRGGNMKErvMeTnrqNgicAYY0LTpk3uuUEDv+/KEoExxoSiffvcc1iY33dlicAYY0JR8+auBHXJkn7fVYEoOmeMMQVOt27uEQB2RmCMMaHo0CFXcC4ALBEYY0woGjwYmjYNyK4sERhjTCg6cADKlw/IriwRGGNMKEpOhnLlArIrSwTGGBOKDhywRGCMMYVacnLAmoas+6gxxoSiIUOgbt2A7MoSgTHGhKL77gvYrqxpyBhjQo0qJCXB4cMB2Z0lAmOMCTUHDkCdOvDqqwHZnSUCY4wJNQcOuGe7j8AYYwqp9ERg3UeNMaaQSk52z3ZGYIwxhVRBOiMQkXgRWSciG0Tk4RzW6Sgiy0VktYjM92c8xhiTL9SvD2PGQOPGAdmd3+4jEJGiwKtAVyAJWCIin6rqmgzrnAO8BsSr6u8iUtVf8RhjTL5Rty78858B250/zwhaAxtUdZOqHgMmAz2zrHMtME1VfwdQ1Z1+jMcYY/KHPXtg/XpITQ3I7vyZCGoBWzNMJ3nmZdQYqCQi34pIoogMzG5DIvIPEUkQkYRdu3b5KVxjjAkRH3wATZrA/v0B2Z0/E4FkM0+zTBcDWgA9gG7AYyJyUqOYqo5V1Zaq2rJKlSp5H6kxxoSSAF8szvEagYhU9vVGVd17im0nAXUyTNcG/shmnd2qehA4KCLfAc2A9afYtjHGFFzJyW7Q+uLFA7I7XxeLd+O+qFM80xl/4Stw3im2vQRoJCLhwDZgAO6aQEYzgP+KSDGgBHAB8ELuQjfGmAIqgGMRgO9E8ArQEfgBmAR8r6pZm3ZypKopInIXMBsoCoxT1dUicrtn+Ruq+ouIfAn8DKQBb6vqqjP7KMYYU0AEcJhKAPH13S4igksG1+B6Ac0BXlfVzQGJLhstW7bUhISEYO3eGGP8b/58+PNP6N8/zzYpIomq2jK7ZT7vI/CcAXwjIstwTTtPA78Cb+VZdMYYYzLr0CGgu/N1sbgsrt9/f6AKMA2IU9WtOb3HGGNMHli+3F0jaNgwILvzdUawE/frfxKwAXeBuJWItAJQ1Wn+D88YYwqhG25w5SU++SQgu/OVCD7yPEd4Hhkp7gzBGGNMXktODpleQ5/Zr35jjAmCAHcf9XVn8b8CFoUxxpgTQigRGGOMCbTjx+Ho0YDeR+CraShCRH7OZr7gepY29VNMxhhTeIm4i8QRWS/N+o+vRLAZuCJQgRhjjAGKFYOrrgrsLn0sO6qqvwUsEmOMMfDKK67raLduAdulr2sEzQIWhTHGGJg2De65B957L6C79ZUIsrs+YIwxJq+pwksvudpCF1wA77wT0N37ahrKdaVRY4wxZygtDbp2hXnzoGdPdzZQunRAQ/B1RtBURPZn80gWkcCMn2aMMQVRaiq89pp7XaQIxMXBf//rmoYqVgx4OL7OCFaqavOARWKMMYXFRx/BkCFwySXuwvDo0UENx24oM8aYQEpJgSeegJiYgFUXPRVfZwRTAxaFMcYUFhMnwvr1rhmoSGj8FvcVRaX0YSUzEpH7RGSkH2MyxpiC6fhxePJJd02gV69gR+PlKxH0AMZmM/8lzzJjjDGnY+NGlwyeftqVkggRPruPqmpaNjPTPGMZG2OMOR0REbBhA5QoEexIMvF1RnBIRBplnemZd9h/IRljTAGxbRusXAlffw033wz790PJkiF1NgC+E8Fw4AsRGSwiMZ7HjcBMzzJjjDHZWbgQOneG2rWhaVPXTXTqVDcWcQjKsWlIVb8QkV7AA8BQz+xVwNWqujIQwRljTL5UrhxUqwYjRrjmoJIloVMnKFs22JFly9c1AlR1FTAofVpEKgF/+TsoY4zJ15o2hUmTgh1FruXYNCQiw0UkwvO6pIjMAzYCO0TkkkAFaIwx+cqzz8KKFcGO4rT4ukbQH1jneT3Is24VoAPwbz/HZYwx+c/PP8Mjj8BXXwU7ktPiKxEcU9X0CqTdgEmqmqqqv3CKJiVjjCl0VOGBB9z1gZtuCnY0p8VXIjgqItEiUgXoBMzJsKyMf8Myxph8ZsIEmDPHNQ1VrhzsaE6Lr1/29wAf45qDXlDVzQAi0h1YFoDYjDEmf9i1C+69Fy68EO64I9jRnDZf3UcXAxHZzJ8FzPJnUMYYk6+EhcGAATB0KBQtGuxoTpvP0neepqH3RCRBRJZ4XscEKjhjjAlZx4/Dww+7+kEi8Oqr7p6BfCjHMwIR6QmMAZ4FngcEaAFME5H7VXVGYEI0xpggWb0atmyBHp46m9Onw759sHUrjB8Pmze7G8fuuy+YUZ41X9cIngK6quqWDPNWeO4nmOF5GGNMwXLgAEyZAm+/DYsWQZs2JxLB8OGudhC4EhIvvwyXXx68WPOIr0RQPEsSAEBVt4hIcf+FZIwxQfLVV3DNNbBnD5x/PvznP3DVVSeWf/mlaxIqW9ZdFyggfCWC4yJSV1V/zzhTROoBKf4NyxhjgqBIEahXD2bMgIsuOrlKaM2awYnLz3wlgseBuSLybyARUKAV8DDwUABiM8aYwEhJgWLFoEsXWLIkZIaQDJQcP62qTgf6Ap2B8cD7uBvL+nmWGWNM/peaCu3bw6hRbrqQJQHw3WuomKquAAYGMB5jjAmsZ56BH3+Ee+4JdiRB4yv1/ZT+QkReCUAsxhgTWJMnw+OPww03QL9+wY4maHwlgoxXSS4+k42LSLyIrBORDSLysI/1WolIqoj0OZP9GGPMaZsyBQYPds1Cb70VcsNHBpKvRKA+lp2SiBQFXgUuAyKBa0QkMof1RgKzz2Z/xhhzWlJSoFUrmDbNjSBWiPnqNRQhIj/jzgwaeF7jmVZVbQYBvGQAABYCSURBVHqKbbcGNqjqJgARmQz0BNZkWW8o8AmuR5IxxgTGdde5ewYK4cXhrHwlgvPPctu1gK0ZppOACzKuICK1gN64nkk5JgIR+QfwD4C6deueZVjGmEKtTx9XJfSf/7Qk4OGr+uhvZ7nt7BrcsjY3vQg8pKqp4qN9TlXHAmMBWrZseVZNVsaYQiglBRISYO1a+OQTVzbCePlzpLEkoE6G6drAH1nWaQlM9iSBMKC7iKTYfQrGmDz1+utw993udY0acOedwY0nxPjzvGgJ0EhEwkWkBDAA+DTjCqoarqr1VbU+bhCcOy0JGGPy3J9/QuvWMHeuKyRXxgZZzMhvZwSqmiIid+F6AxUFxqnqahG53bP8DX/t2xhjMnnmGXj6absmkANfdxavxEcX0lz0Gsp2NLOcEoCqDj7V9owx5rT99Recc44lAR98nRGkF9ke4nn+wPN8HXDIbxEZY8zZ2rcPPvrIDSDz9ttw/fUwZkywowpZp+w1JCIXq2rGO4sfFpEfcAPXGGNMaHn5ZXjoIThyxJ0FlCjhqoqaHOXmXKmsiLRNnxCRi4Cy/gvJGGPOQoMGbtSwxERXWfTwYbjssmBHFdJyc7H4ZmCciFTEXTP4G7jJr1EZY0xu7N/vCsd9+CH8/TcsW+aGlUwfWtLkyikTgaomAs1EpAIgqvq3/8Myxhgf1qxxbf5TpsChQxAVBU2aBDuqfOuUTUMiUk1E3gGmqOrfIhIpIjcHIDZjjHF27ID+/d0FYIA//oCpU129oMWL3YDyn3wS3Bjzsdw0DY0H3gUe9UyvB6YA7/gpJmOMgVWr4Ouv4dgxePFFN6D8gAFuWefOsH07lCsX3BgLiNwkgjBV/UhEhoH3RrFUP8dljCmMDh+G0qXd6yefhI8/dq8bNYIvvoCmntuXihSxJJCHctNr6KCInIvn5jIRaYO7YGyMMXlnwQLX42fdOjc9ahQkJbkbwn755UQSMHkuN2cE/4erEdTAc/9AFdyg9sYYc3bmz4fnnnNnAt9/7xJBeiXi8PDgxlaI5CYRrAY6AE1wpaXX4d9idcaYgmzlSqhdGypWhHvvhd273Zf+TTe5nkAVKgQ7wkInN4ngR1WNwyUEAERkKRDnt6iMMQVTYqIbI/iHHyA2FqZPh0qV7Ms/yHwVnauOG2WstIg058RAMxUAq+FqjMmdvXthzhzX/PPooxAWBlWrumX16gU3NgP4PiPoBgzGDSjzPCcSwX7gEf+GZYzJt9atg3HjXL//uDh44w2XAMBVAZ0zB2rWDG6MJhNfRefeE5EPgGtUdWIAYzLG5FcjR8LDD0OxYnDFFW7erbfCJZdAlSrubKB8+eDGaE7i86KvqqYBtwUoFmNMfjZliksC/fq58s9tPbUqq1Rxo4OFh1sSCFG5uVj8lYjcj7ub+GD6TFXd67eojDH5y+LF7q7ftm3h/fehZMlgR2ROQ266gd6EG5zmOyDR80jwZ1DGmHzg2DH44ANQhVq13P0AM2ZYEsiHclN91O7qMMaccOwYvPcevP66K/vcoAFcdJEbDMbkS6dMBCJSHLgDaO+Z9S3wpqoe92NcxphQlJQEffq4pqDoaNcMdNFFwY7KnKXcXCN4HSgOvOaZvsEz7xZ/BWWMCRG7d7tf/e3buyafG2+E1atdOeg+fU6UgzD5Wm4SQStVbZZhep6IrPBXQMaYIEtLg7lz3aDv06fD8eNuPICqVV2voFq1ICIi2FGaPJSbRJAqIg1UdSOAiJwHWBlqYwqixES4+mr47Tc491wYMgSuvNKVgQAbBL6Ayk0ieAD4RkQ24e4urgfc6NeojDGBdeQIlCoFDRu6tv9Ro6BnT+sBVEjkptfQ1yLSiBPVR9eq6lG/R2aM8a+0NPjmG3jrLfj5Z1i0yFUE/fzzYEdmAizH+whEpJWn8ByeL/5Y4ClgtIhUDlB8xpi8lpbmhn5s1MiVfpg9Gy67zJWFMIWSrxvK3gSOAYhIe+A54H3c6GRj/R+aMcYvZs+G++6DOnVg4kQ3EPzzz0MZKypcWPn6CVA0QxmJ/sBYVf0E+ERElvs/NGOMX1x2GcybBx07WvdPA/g+IygqIumJogswL8MyO4c0Jr9Ztgy++8697tTJkoDx8vWFPgmYLyK7gcPAAgARaYgNXm9M/vPPf8KaNbBli+shZIyHr/EInhGRr4EawBxVVc+iIsDQQARnjDkD27a5m74Abr7ZlYNIS4NffnEXiS0JmCx8NvGo6qJs5q33XzjGmLMyY4YbD2DmTNcjqG5d+NtzAt++Pdxmw4uYk1lbvzEFxfr1cO21blD49EJwjz8e3JhMvpCb8QiMMfnB88+7JqDp060rqDktlgiMKQj27nWDxFx3HdSoEexoTD5jicCYgmDJEvc81PpxmNNn1wiMyc+OHoX9+6FbN/jzT6hQIdgRmXzIEoEx+UlqKkyb5gaMWb/eNQddfjmMH29JwJwxvyYCEYkHXgKKAm+r6nNZll8HpA90egC4Q1Vt0BtjcvKf/8CDD7rXxYtDr15www3Bjcnke35LBCJSFHgV6AokAUtE5FNVXZNhtc1AB1XdJyKX4YrZXeCvmIzJ91q3hnvugWHDoGxZKFcu2BGZAsCfZwStgQ2quglARCYDPQFvIlDVhRnWXwTU9mM8xuRfqq42UIcO7mFMHvJnr6FawNYM00meeTm5GfgiuwUi8g8RSRCRhF27duVhiMbkE/fc42oFeSu9GJN3/JkIsittmO1fsYh0wiWCh7JbrqpjVbWlqrasUqVKHoZoTAhJTXU3hIErC/H662784K5d4ZVX3HyrGGr8wJ+JIAmok2G6NvBH1pVEpCnwNtBTVff4MR5jQsfx4zBnDnz0kSsE17y5GyFsrGfMp0WL4M47Yfly2LXLDSI/alRwYzYFlj+vESwBGolIOLANGABcm3EFEakLTANusGJ2plA4fBhKl4bkZLjySncfAEBcHDz2GLRq5aYvuQQSE12CsLMA42d+SwSqmiIidwGzcd1Hx6nqahG53bP8DWA4cC7wmrg/9hRVbemvmIwJqBUrYO1aVxF0+nQ4csT1+lm0CCpXhu+/dzWBSpeG8PDM7y1a1CUHYwLAr/cRqOosYFaWeW9keH0LcIs/YzAmoPbvP3Fj15tvunb+c85xff2rVoVzz3XXAYoUgZb2m8eEBruz2Ji8sm0btGkD999/oq//bbdB48buV78xIcoSgTFnY+ZMV/IBYOFC+OsvNyg8QJ067mFMiLNEYIwvf/3lLtqmpsLXX8Pkya5r5+efQ9u2sGmT6/0D7lf/1KnQrFlwYzbmNFkiMCartDQ4cMC19W/a5HrwgLuA2707nHceVK/u5g0daqWfTb5nicCYgwdh2TL3q//772HcOFfGYdw413Nn/nx3cbdRI6hWLdjRGpPnLBGYgi81FebOhdmz3WtwX+wvvOBeX3qpa99P17kz9OhxYrp9+8DFakwQWCIwBd+jj8LIkVCqlHuAa+ZJTwSPP+5u7Cpb1vXnz9qn35gCzhKBKRhU3Y1an38Ox47BypUwfDhcdBEMGgQtWrg7eUuWPPm9l14a+HiNCSGWCEz+t2OHu4i7dKn7pV+ypLt5a+dOt/z8893DGJMtG7ze5D9paa7L5r/+5co2hIVBZKS7k3ffPnfxd/NmN3qXMeaU7IzA5C/ffQeDB7svehF44AHX7v/BB8GOzJh8yxKBCR3JyfDxx/D77yfmNW4M11zjXj/wgLvA26CBu7GrV6/s2/yNMafFEoEJvKNHM3/ZN2zontu0gTVrMq97+eUnEsHEiXDVVfD22ycKuxljzpolAhNY+/bBxRfDL7+cmJeW5pp5nnsOqlRxpZqzq8G/bZvV5jfGDywRmMA5dsz9ot+40Q29WKlS5uVXXOH7/ZYEjPELSwTG/3bscKUZUlPdL/533oHrrw92VMYYD0sE5szs2gUpKSemK1Rwd+bOmQNPP32ilMOxY66Oz6pVri//lCn2y96YEGP3EZjT89tvbmStqlWhZs0Tj/fec8vr1HEJolw596hcGR555EQzkCUBY0KOnRGYnB04cGJw9SJF3Jd5zZru+bnn3BCM6dq2dc/nnw8//hj4WI0xZ8wSgTlZWhqMGAFPPuleg/ulv2yZG3P3q6+CG58xJk9ZIjAn+/JLV5Gzf3/X1TPdnj0uERhjChRLBIXZ/v3uS//QIdels0gRdxbQvbsblrFTJ2vTN6YQsERQmKSkuOqcIvDEEzB6tEsC4Ob16XNi3c6dgxKiMSbwLBEUVNu2uSEWwVXj/Ogj9yt/40Y38EqVKnDtta6AW61arodPWFhQQzbGBIclgoJG1f26X74crrvuxPx69eD++0904xwyJDjxGWNCjiWC/GzxYle8bcECV5Bt715XwuHjj93g6+vWufWKFIHzznPPxhiThSWC/CQlxX3pd+rkpl980ZVjLlECeveGiAh3d+/Ro66pp3Hj4MZrjMkXLBHkFz//DH37wvr1kJgIcXEwapQbmL1WrZMLuBljTC5ZW0F+MGGCq9WfnAxTp0JMjJtfpw5ER1sSMMacFTsjCBU7d7o7d9NddBGULw8PPui6ebZv7wq2Va8evBiNMQWSnREE2p9/uiadRYvc9Lp1bhSuWrUgPv7EY/Nmt7x5c/i//4O5cy0JGGP8ws4I/CW9Gye4Add37XJNPJ995ko0T53qlu3d60bruuceNzBLiRJufoMG7vmaa04M1WiMMX5giSAvpaTArFlu4JWwMPcM7gt+/35Xuvmf/4SbboImTdyyCy90N3kZY0yQWCLIC9dcA99848o1JCe7Jpw77jixfOZM14e/VSsoXjx4cRpjTDYsEeSFwYOhZEnXh//SS13Rtoxf+Om1+o0xJgRZIjhT69a5m7tuuQW6dXMPY4zJhywR5Nby5e4u3nRTp7pmoKuucsMxGmNMPlW4E4Gq679frZqbXrMGjh93JRo+/RSmT4cbboCHHnIXf196yb0H3Jf/p59aEjDG5HuFMxHs2wcffghvv+0u7H7xhZvfrRskJbnXRYpAly5Qt66brl0bDh8OTrzGGONHfk0EIhIPvAQUBd5W1eeyLBfP8u7AIWCwqi71W0AjRsCkSbBpExw54ur1DB58Yvlbb7kvexFo2dJ9+RtjTAHnt0QgIkWBV4GuQBKwREQ+VdU1GVa7DGjkeVwAvO559o/q1SEyEi65BAYNcokgo/h4v+3aGGNClT/PCFoDG1R1E4CITAZ6AhkTQU/gfVVVYJGInCMiNVR1u18iuuUW9zDGGOPlz1pDtYCtGaaTPPNOdx1E5B8ikiAiCbt27crzQI0xpjDzZyKQbObpGayDqo5V1Zaq2rJKlSp5EpwxxhjHn4kgCaiTYbo28McZrGOMMcaP/JkIlgCNRCRcREoAA4BPs6zzKTBQnDbA3367PmCMMSZbfrtYrKopInIXMBvXfXScqq4Wkds9y98AZuG6jm7AdR+90V/xGGOMyZ5f7yNQ1Vm4L/uM897I8FqBIf6MwRhjjG82QpkxxhRylgiMMaaQE9WTemuGNBHZBfx2hm8PA3bnYTj+lp/itVj9Iz/FCvkr3sIWaz1Vzbb/fb5LBGdDRBJUtWWw48it/BSvxeof+SlWyF/xWqwnWNOQMcYUcpYIjDGmkCtsiWBssAM4TfkpXovVP/JTrJC/4rVYPQrVNQJjjDEnK2xnBMYYY7KwRGCMMYVcoUkEIhIvIutEZIOIPBzseDISkToi8o2I/CIiq0XkHs/8J0Rkm4gs9zy6BztWABHZIiIrPTEleOZVFpGvRORXz3OlYMcJICJNMhy/5SKyX0TuDZVjKyLjRGSniKzKMC/HYykiwzx/w+tEpFsIxDpaRNaKyM8i8j8ROcczv76IHM5wfN/IecsBizXHf/NgHlcf8U7JEOsWEVnumZ/3x1ZVC/wDV/RuI3AeUAJYAUQGO64M8dUA4jyvywPrgUjgCeD+YMeXTbxbgLAs80YBD3tePwyMDHacOfwd/AnUC5VjC7QH4oBVpzqWnr+JFUBJINzzN100yLFeChTzvB6ZIdb6GdcLkeOa7b95sI9rTvFmWf48MNxfx7awnBF4h81U1WNA+rCZIUFVt6vqUs/rZOAXshmpLcT1BN7zvH4P6BXEWHLSBdioqmd6Z3qeU9XvgL1ZZud0LHsCk1X1qKpuxlXtbR2QQMk+VlWdo6opnslFuDFFgi6H45qToB5X8B2viAjQD5jkr/0XlkSQqyExQ4GI1AeaA4s9s+7ynHaPC5XmFtwocnNEJFFE/uGZV009Y0l4nqsGLbqcDSDzf6ZQPLaQ87EM9b/jm4AvMkyHi8gyEZkvIu2CFVQW2f2bh/pxbQfsUNVfM8zL02NbWBJBrobEDDYRKQd8AtyrqvuB14EGQCywHXd6GAouVtU44DJgiIi0D3ZApyJucKQrgameWaF6bH0J2b9jEXkUSAEmemZtB+qqanPg/4APRaRCsOLzyOnfPGSPq8c1ZP4Bk+fHtrAkgpAfElNEiuOSwERVnQagqjtUNVVV04C3CPDpak5U9Q/P807gf7i4dohIDQDP887gRZity4ClqroDQvfYeuR0LEPy71hEBgGXA9eppxHb08yyx/M6Edfu3jh4Ufr8Nw/J4wogIsWAq4Ap6fP8cWwLSyLIzbCZQeNpA3wH+EVV/5Nhfo0Mq/UGVmV9b6CJSFkRKZ/+GnexcBXueA7yrDYImBGcCHOU6VdVKB7bDHI6lp8CA0SkpIiEA42An4IQn5eIxAMPAVeq6qEM86uISFHP6/NwsW4KTpTemHL6Nw+545rBJcBaVU1Kn+GXYxvIK+PBfOCGxFyPy56PBjueLLG1xZ2K/gws9zy6Ax8AKz3zPwVqhECs5+F6WKwAVqcfS+Bc4GvgV89z5WDHmiHmMsAeoGKGeSFxbHHJaTtwHPfL9GZfxxJ41PM3vA64LARi3YBrX0//u33Ds+7Vnr+PFcBS4IoQiDXHf/NgHtec4vXMHw/cnmXdPD+2VmLCGGMKucLSNGSMMSYHlgiMMaaQs0RgjDGFnCUCY4wp5CwRGGNMIWeJwJhTEJHqIjJZRDaKyBoRmSUiOd7AIyIHPM/1M1aTNCZUWSIwxgfPzX7/A75V1QaqGgk8AlQLbmTG5B1LBMb41gk4rqremu+qulxVF4jIAyKyxFPE7ElfGxGRKBH5yVM//mcRaeT3yI3JJUsExvgWDSRmnSkil+Ju7W+NK2LW4hTF924HXlLVWKAl7u5RY0JCsWAHYEw+dannscwzXQ6XGL7LYf0fgUdFpDYwTTOXFDYmqOyMwBjfVgMtspkvwLOqGut5NFTVd3LaiKp+iCuDfRiYLSKd/ROuMafPEoExvs0DSorIrekzRKQVsB+4yTOGBCJSS0RyHIzHUyVyk6q+jCt41tS/YRuTe9Y0ZIwPqqoi0ht4UUQeBo7gxmy+F/gL+NF1LOIAcD05j8PQH7heRI7jxk1+ys+hG5NrVn3UGGMKOWsaMsaYQs4SgTHGFHKWCIwxppCzRGCMMYWcJQJjjCnkLBEYY0whZ4nAGGMKuf8H/PuiQwzs09YAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "legend_label=\"SCOT alignment FOSCTTM \\n average value: \"+str(np.mean(fracs))\n", + "plt.plot(np.arange(len(fracs)), np.sort(fracs), \"r--\", label=legend_label)\n", + "plt.legend()\n", + "plt.xlabel(\"Cells\")\n", + "plt.ylabel(\"Sorted FOSCTTM\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualizing Alignment" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from sklearn.decomposition import PCA\n", + "\n", + "pca=PCA(n_components=2)\n", + "Xy_pca=pca.fit_transform(np.concatenate((X_new, y_new), axis=0))\n", + "X_pca=Xy_pca[0: 177,]\n", + "y_pca=Xy_pca[177:,]\n", + "\n", + "### Read in cell type information:\n", + "# Xlabels=np.genfromtxt(\"data/scGEM_typeExpression.txt\")\n", + "# ylabels=np.genfromtxt(\"data/scGEM_typeMethylation.txt\")\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2)\n", + "ax1.scatter(X_pca[:,0], X_pca[:,1], c=\"k\", s=15, label=\"Gene Expression\")\n", + "ax1.scatter(y_pca[:,0], y_pca[:,1], c=\"r\", s=15, label=\"DNA Methylation\")\n", + "ax1.legend()\n", + "ax1.set_title(\"Colored based on domains\")\n", + "\n", + "ax2.scatter(X_pca[:,0], X_pca[:,1], cmap=\"rainbow\")# , c=Xlabels, s=15)\n", + "ax2.scatter(y_pca[:,0], y_pca[:,1], cmap=\"rainbow\")# , c=ylabels, s=15)\n", + "ax2.set_title(\"Colored based on cell type\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}