--- a +++ b/.ipynb_checkpoints/SNAREseq_replicate-checkpoint.ipynb @@ -0,0 +1,190 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Notebook for running SCOT on SNARE-seq Cell Mixture Data\n", + "Access to the raw dataset: Gene Expression Omnibus accession no GSE126074 \n", + "SNARE-seq data in `/data` folder containes the version with dimensionality reduction techniques applied from the original SNARE-seq paper (https://www.nature.com/articles/s41587-019-0290-0) \n", + "SCOT software has been updated on 20 September 2020. It now outputs error statements for convergence issues. 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", + "If you have any questions, e-mail: ritambhara@brown.edu, pinar_demetci@brown.edu, rebecca_santorella@brown.edu " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Import source code:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import src.utils as ut\n", + "import src.evals as evals\n", + "from src.scot import *" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Read in the data:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Dimensions of input datasets are: X= (1047, 10) y= (1047, 19)\n" + ] + } + ], + "source": [ + "X=np.exp(np.load(\"data/scrna_feat.npy\"))\n", + "y=np.load(\"data/scatac_feat.npy\")\n", + "print(\"Dimensions of input datasets are: \", \"X= \", X.shape, \" y= \", y.shape)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Perform normalization (optional):" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "X=ut.unit_normalize(X)\n", + "y=ut.unit_normalize(y)\n", + "\n", + "## If you'd like to apply z-score normalization instead:\n", + "# X=ut.zscore_standardize()\n", + "# y=ut.zscore_standardize()\n", + "# Note that zscore_standardize doesn't yield as good results on this dataset and MMD-MA and UnionCom comparisons \n", + "# also used unit (l-2) normalization" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Set hyperparameters of the algorithm:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "It. |Err \n", + "-------------------\n", + " 0|1.425627e-03|\n", + " 10|1.977665e-04|\n", + " 20|7.060937e-05|\n", + " 30|8.834972e-06|\n", + " 40|7.413147e-07|\n", + " 50|6.468501e-08|\n", + " 60|5.819180e-09|\n", + " 70|5.293865e-10|\n" + ] + } + ], + "source": [ + "# Set hyperparameters of the algorithm:\n", + "k=24\n", + "e=0.0038 \n", + "# Other values to try for very similar alignment results:\n", + "# k=25 with e=0.0018, 0.00182, 0.00185, 0.002, or k=30 with \n", + "# Combinations from a range of k=20 to k=30 and e=0.0015 to e= 0.0040 (and aounrd 0.01 for k=30 setting) seem to \n", + "# yield the best results on this dataset, so if you'd like to perform hyperparameter tuning, \n", + "# you can set a grid between these values.\n", + "X_new,y_new= scot(X, y, k, e)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Average FOSCTTM score for this alignment is: 0.19851720567368114\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": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1sUlEQVR4nO3dd3hUdfb48fchNOlVhNBCMwkEQggYUBZQmqyiLILgKij81K8IousidtTVteHaUUEQV1lAEYRFVFQEUVFIKCKhhbISQEpEOpKQ8/vjMxlCSIYAmcwkc17PM09yy9x7LpfMmfupoqoYY4wJXSUCHYAxxpjAskRgjDEhzhKBMcaEOEsExhgT4iwRGGNMiCsZ6ADOVo0aNbRhw4aBDsMYY4qUpKSkvapaM7dtRS4RNGzYkMTExECHYYwxRYqI/C+vbVY0ZIwxIc4SgTHGhDhLBMYYE+KKXB1BbtLT00lNTeXYsWOBDsUYk4uyZctSt25dSpUqFehQTC6KRSJITU2lYsWKNGzYEBEJdDjGmGxUlbS0NFJTU4mIiAh0OCYXxaJo6NixY1SvXt2SgDFBSESoXr26PbEHsWKRCABLAsYEMfv7DG7FJhEYY4w5N5YICshTTz1F8+bNadmyJbGxsfz444+Aq8i+//77adq0KXFxcbRv355PP/0UgP379zNo0CCaNGlC48aNGTRoEPv372f16tXExsYSGxtLtWrViIiIIDY2lq5du+Y7noYNG7J3714AOnToUPAXfJ5WrlzJvHnzct22cOFCKleu7P03yH7d48ePJzIyksjISNq1a8e3337r3TZ37lxat25Nq1atiI6O5q233vJu+/e//02LFi2IiYmhdevWjB07ljvvvJPY2Fiio6O54IILvOcTEWJjY2nSpMkpcXz//fd07tyZ+vXrk30ej2uvvZYKFSr44V/JGI89e+Cyy2DuXP8cX1WL1KtNmzaaU3Jy8mnrCtP333+vCQkJeuzYMVVV3bNnj27fvl1VVUePHq2DBg3ybvv11191+vTpqqrat29fHTNmjPc4jz76qF533XWnHHvw4MH64YcfnnVMDRo00D179pzL5RSKd955R++8885ct3399df65z//+bT1//3vfzUuLs57XUlJSVqvXj3duXOnHj9+XGvXrq3btm1TVdVjx47punXrVFV13rx52rp1a+89OXbsmI4fP9573C1btmjz5s3zFUenTp00JiZGFy9erKqq+/bt03bt2mn58uXP9p8g5AT677RIO3xY9YMPVFNSzvkQQKLm8blqTwQFYOfOndSoUYMyZcoAUKNGDerUqcORI0eYMGECr776qndbrVq16N+/PykpKSQlJfHII494j/Poo4+SmJjIpk2b8n3ua6+9ljZt2tC8eXPGjx+f6z5Z31YzMzMZNmwYkZGRdOvWjV69ejFjxgzAPUGMGTOGuLg4YmJiWLduHQCPPfYYgwcPpmPHjjRo0ICZM2dy3333ERMTQ8+ePUlPTwcgKSmJTp060aZNG3r06MHOnTsB6Ny5M6NHj6Zdu3Y0a9aMxYsXc/z4cR599FGmT59ObGws06dPz9e1Pvvsszz//PPUqFEDgLi4OAYPHszrr7/OwYMHycjIoHr16gCUKVOGiy++GICnn36asWPHUqdOHe+2W2+9Nd//xjkNGDCAadOmATBz5kz+8pe/nPOxjMmXcuWgXz9o3Ngvh7dEUAC6d+/Otm3baNasGcOGDWPRokUApKSkUL9+fSpVqnTae5KTk4mNjSUsLMy7LiwsjNjYWNasWZPvc0+aNImkpCQSExN55ZVXSEtLy3PfmTNnsnXrVpKTk3nvvfdYsmTJKdtr1KjB8uXLueOOOxg7dqx3/aZNm1iwYAFz5szhxhtvpEuXLqxevZoLLriATz75hPT0dEaMGMGMGTNISkpiyJAhPPTQQ973Z2RksHTpUl566SUef/xxSpcuzRNPPMH111/PypUruf7660+LdfHixd4imaeeegqANWvW0KZNm1P2i4+PZ82aNVSrVo3evXvToEEDBg4cyJQpU8jMzATg559/Pu195+OKK67gm2++4cSJE0ybNi3X+I0pUAcPwhdfwO7dfjl8sehHEGgVKlQgKSmJxYsX8/XXX3P99dfzzDPPEBcX5/dzv/LKK8yaNQuAbdu2sXHjRu+34py+/fZb+vXrR4kSJbjooovo0qXLKduzvtm2adOGmTNnetdfeeWVlCpVipiYGE6cOEHPnj0BiImJYevWraxfv56ff/6Zbt26AXDixAlq166d63G3bt2ar+vq2LEjc8+yPPTtt99m9erVfPnll4wdO5YvvviCyZMnn9Ux8iMsLIzLLruMadOmcfToUWw0XON3mzZB9+4wcyb06VPgh7dEUEDCwsLo3LkznTt3JiYmhnfffZf+/fvzyy+/cODAgdOeCqKjo1m5ciWZmZmUKOEezDIzM1m5ciXR0dH5OufChQv58ssvWbJkCeXKlaNz587n1VY7q/gqLCyMjIyM09aXKFGCUqVKeZsClihRgoyMDFSV5s2bn/aEcabjnq3o6GiSkpK4/PLLveuSkpJo3ry5dzkmJoaYmBhuuukmIiIimDx5Ms2bNz/tfedrwIAB9OnTh8cee6zAjmlMnjxFsJQu7ZfDW9FQAVi/fj0bN270Lq9cuZIGDRpQrlw5hg4dysiRIzl+/DgAe/bs4cMPP6RJkya0bt2aJ5980vu+J598kri4OJo0aZKv8+7fv5+qVatSrlw51q1bxw8//OBz/0svvZSPPvqIzMxMdu3axcKFC8/+YnNx8cUXs2fPHm8iSE9PP2PxVsWKFTl48OBZnee+++5j9OjR3uKvlStXMnnyZIYNG8ahQ4dOuZ6sewDwwAMPMGrUKH799VcAjh8/zttvv31W586pY8eOPPDAAwwcOPC8jmNMvng+P/DTEB32RFAADh06xIgRI/j9998pWbIkTZo08VbcPvnkkzz88MNER0dTtmxZypcvzxNPPAHAxIkTGTFiBI09FUDt27dn4sSJ+T5vz549efPNN4mKiuLiiy8mISHB5/59+/blq6++Ijo6mnr16hEXF0flypXP8apPKl26NDNmzOCuu+5i//79ZGRkcPfdd5/yTT2nLl268MwzzxAbG8sDDzyQr3L23r17s337djp06ICIULFiRd5//31q167NwYMHee6557j99tu54IILKF++vLdYqFevXuzatYuuXbuiqogIQ4YMOa9rFhH+/ve/n9cxjMm3rCcCPyUC0WztoYuC+Ph4zTkxzdq1a4mKigpQREXLoUOHqFChAmlpabRr147vvvuOiy66KNBhmRBgf6fnYf586NEDvv0WLr30nA4hIkmqGp/bNnsiCDFXXXUVv//+O8ePH+eRRx6xJGBMURAXB/PmQT7rD89W8UwEnTufvq5/fxg2DI4cgV69Tt9+883utXcvXHfdqdsKqCw9GBRUvYAxphDVqAFXXum3w1tlsTHGBLtffoFZs1x/An/Iq8txsL6CcYiJUOOv4SsmT56sTZo00SZNmujkyZNz3eeDDz7Q6OhoFRFdtmyZd/0ff/yhN998s7Zo0UJbtmypX3/9tXfbf/7zH23RooXGxMRojx49vLGPGTNG69Spo61atdJWrVrpJ598oqqqe/fu1c6dO2v58uVPGQbjwIED3n1btWql1atX15EjR6qq6gsvvKBRUVEaExOjl19+uW7duvWM19WpUydt1qyZ93i7du3ybps+fbpGRUVpdHS0Dhw4UFVVFyxYcMr5y5Qpo7NmzVJV1SFDhmjLli01JiZG+/btqwcPHvR5LFXVUaNGaXR0tEZGRuqIESM0MzNTVVUffPBBrVu37mnDZixatEhbt26tYWFhuQ57sn//fg0PD89z6BD7Oz0P772nCqobNpzzIfAxxETAP9jP9hWqiSAjIyPQIXj5IxGkpaVpRESEpqWl6W+//aYRERH622+/nbZfcnKyrlu3Tjt16nRKInjttdf05ptvVlXVXbt2aVxcnJ44cULT09O1Zs2a3nhHjRrlHd9pzJgx+vzzz592jkOHDunixYv1jTfeyPNDTVU1Li5OFy1apKruQ/rw4cOqqjpu3Djt37//Ga8r5zVk2bBhg8bGxnr3y54gsv97Va1a1XvO/fv3e7fdc889+vTTT/s81nfffacdOnTQjIwMzcjI0ISEBG/yXLJkie7YseO0RLBlyxZdtWqV3nTTTbkmgrvuuksHDhxoicAfJk1yH9dbtpzzIXwlAisa8oM77riD+Ph4mjdvzpgxYwD47LPP6Nevn3efhQsXctVVVwEwf/582rdvT1xcHP369ePQoUOAG/9n9OjRxMXF8eGHHzJhwgTatm1Lq1at6Nu3L0eOHAHcEBAJCQnExMTw8MMPnzIS5vPPP0/btm1p2bKlN5bs3nzzTUaNGuVdnjx5MsOHDwfOPI7R1q1badGihXd57Nix3g5WmzZtomfPnrRp04aOHTt6xy7Ky+eff063bt2oVq0aVatWpVu3bnz22Wen7ZfVVDan5ORkb4exCy+8kCpVqpCYmOj9j3748GFUlQMHDnjHHMpL+fLlueyyyyhbtmye+2zYsIHdu3fTsWNHwDWHLVeuHAAJCQmkpqae1XVlN2HCBO68806qVq3qvZ6cZsyYwZVXXuk9Z1aHRVXl6NGj3k5/eR1LRDh27BjHjx/njz/+ID09nVq1annjz94zPEvDhg1p2bKltwNkdklJSezatYvu3bv7vDZzjg4fdj8997ugWSLwg6eeeorExER++uknFi1axE8//UTXrl358ccfOey5odOnT2fAgAHs3buXJ598ki+//JLly5cTHx/Pv/71L++xqlevzvLlyxkwYAB/+ctfWLZsGatWrSIqKsrb52DkyJGMHDmS1atXU7duXe9758+fz8aNG1m6dCkrV64kKSmJb7755pRY+/bt6x2iIntccHbjGOV022238eqrr5KUlMTYsWMZNmwYAHPmzOHRRx89bf/t27dTr14973LdunXZvn17vs/XqlUr5syZQ0ZGBlu2bCEpKYlt27ZRqlQp3njjDWJiYqhTpw7JyckMHTrU+77XXnuNli1bMmTIEPbt25fv82WNMZTbhCsTJ07kSk/F3pmu65ZbbiE2NpZ//OMf7hEdl2Q2bNjApZdeSkJCQq6JY9q0aad1Zrvlllu46KKLWLduHSNGjPB5rPbt29OlSxdq165N7dq16dGjxzk37czMzOTee+89ZXwqU8Cy6gYqVvTL4S0R+MEHH3xAXFwcrVu3Zs2aNSQnJ1OyZEl69uzJf//7XzIyMvjkk0+45ppr+OGHH0hOTubSSy8lNjaWd999l//973/eY2XvaPXzzz/TsWNHYmJimDJlirf37pIlS7xPGzfccIN3//nz5zN//nxat25NXFwc69atO6UHNEDNmjVp1KgRP/zwA2lpaaxbt45LPe2UX3nlFVq1akVCQoJ3HKP8OHToEN9//z39+vUjNjaW22+/3Tsaae/evb0d6grSkCFDqFu3LvHx8dx999106NCBsLAw0tPTeeONN1ixYgU7duygZcuWPP3004B7ctu0aRMrV66kdu3a3Hvvvfk+X24fxADvv/8+iYmJpzxl5WXKlCmsXr2axYsXs3jxYt577z3ADdK3ceNGFi5cyNSpU7n11lv5/fffve/buXMnq1evpkePHqcc75133mHHjh1ERUV5R3TN61gpKSmsXbuW1NRUtm/fzoIFC1i8eHG+rz+7cePG0atXr1O+hJgCdvAghIWBj6fU81E8m48G0JYtWxg7dizLli2jatWq3Hzzzd7xfwYMGMBrr71GtWrViI+Pp2LFiqgq3bp1Y+rUqbker3z58t7fb775Zj7++GNatWrF5MmTz9gUVFV54IEHuP32233uN2DAAD744AMiIyPp06cPIpKvcYxKlizpHeET8G7PzMykSpUqrFy50ud5swsPDz/lelJTU+mcWzPgPJQsWZIXX3zRu9yhQweaNWvmjSGr93b//v155plnALxFIQC33nqrt6juTFatWkVGRsZpI5p++eWXPPXUUyxatMg7vpKv6woPDwfccBs33HADS5cuZdCgQdStW5dLLrmEUqVKERERQbNmzdi4cSNt27YF3BeNPn36UCqXXqZhYWEMGDCA5557jltuuSXPYy1cuJCEhARvMeKVV17JkiVLvEVdZ2PJkiUsXryYcePGcejQIY4fP06FChW8/86mANxxh2v27qcpP+2JoIAdOHCA8uXLU7lyZXbt2uWdjQygU6dOLF++nAkTJniLXxISEvjuu+9ISUkB4PDhw2zYsCHXYx88eJDatWuTnp7OlClTvOsTEhL46KOPALzj5AP06NGDSZMmeesctm/fzu5chrHt06cPs2fPZurUqd648jOOUa1atdi9ezdpaWn88ccf3tFCK1WqREREBB9++CHgEtKqVat8/rv16NGD+fPns2/fPvbt28f8+fNP+8bry5EjR7zFbl988QUlS5YkOjqa8PBwkpOT2bNnj3dbVhFI1lMKwKxZs06p7/Bl6tSppz0NrFixgttvv505c+acUqaf13VlZGR4Z5BLT09n7ty53vNfe+213uSxd+9eNmzYQKNGjfI8v6p6//+oKnPmzCEyMtLnserXr8+iRYvIyMggPT2dRYsWnXPR0JQpU/jll1/YunUrY8eOZdCgQZYEClq9em6GMn/JqxY5WF9FodXQ4MGDtWnTpnr55Zdrnz599J133vFuu/POO7V8+fLe1h6qql999ZXGx8drTEyMxsTE6OzZs1X19NY548aN04YNG2rbtm11+PDhOnjwYFV1LUPatWunMTExOmrUKK1Tp473PS+99JK2aNFCW7RooQkJCZqSxwxHf/7znzUiIsK7fOzYMe3Zs6dGRkbqNddco506dfK2Kske18svv6yNGjXSjh076uDBg70tcjZv3qw9evTQli1balRUlD7++OOqqjp79mx95JFHco1h4sSJ2rhxY23cuLFOmjTJu37o0KHe1jUzZ87U8PBwLV26tF544YXavXt3VXUtWpo1a6aRkZF6xRVXnNJ884033tDIyEiNiYnRq666Svfu3auqqjfeeKO3WenVV1+tO3bs8L6nQYMGWrVqVS1fvryGh4frmjVrvNsiIiJ07dq1p8R+xRVX6IUXXuht2nn11Vf7vK5Dhw5pXFycxsTEaHR0tN51113elmGZmZl6zz33aFRUlLZo0UKnTp3qPdaWLVu0Tp06euLECe+6EydOaIcOHbRFixbavHlzveGGG7ytiPI6VkZGht52220aGRmpUVFRes8993iPN2rUKA0PD1cR0fDwcO89Xbp0qYaHh2u5cuW0WrVqGh0dfdo99DXzXLD9nRYpn3+uOnfueR0CH62G/DrWkIj0BF4GwoC3VfWZHNvrA+8CVTz73K+quU9k62FjDZ3uyJEjXHDBBYgI06ZNY+rUqcyePTvQYRlzilD/Oz0vPXrA/v1whhGGfQnIWEMiEga8DnQDUoFlIjJHVZOz7fYw8IGqviEi0cA8oKG/YiqukpKSGD58OKpKlSpVmDRpUqBDMsYUpIMH/dZiCPxbWdwOSFHVzQAiMg24BsieCBTImrGlMrDDj/EUWx07djxjGbwxpgg7eBCyNW4oaP6sLA4HtmVbTvWsy+4x4EYRScU9DYzI7UAicpuIJIpIYlalX07+LOIyxpwf+/s8T35+Igh0q6GBwGRVrQv0At4TkdNiUtXxqhqvqvE1a9Y87SBly5YlLS3N/rMZE4RUlbS0NJ89tc0ZFOGioe1AvWzLdT3rshsK9ARQ1SUiUhaoAZzextGHunXrkpqaSl5PC8aYwCpbtqx1ODsf338P2foUFTR/JoJlQFMRicAlgAHADTn2+QW4ApgsIlFAWeCsP82zOssYY0yxlMv4WgXJb0VDqpoBDAc+B9biWgetEZEnRKS3Z7d7gVtFZBUwFbhZrXzHGGNO2rABnn8ecukMWlD8OsSEp0/AvBzrHs32ezJwbhNwGmNMKFi+HO67D666CnIZibYgBLqy2BhjjC8HDriflSr53u88WCIwxphgtn+/+1m5st9OYYnAGGOC2f79bghqP7YaskRgjDHBbPt2uOgivw1BDTYfgTHGBLe334ZsQ6b7gz0RGGNMMAsLAz93xrNEYIwxwerYMRg+3PUs9iNLBMYYE6zWrIHXX4dNm/x6GksExhgTrBYvdj+7dvXraSwRGGNMsMrqTJbLqMsFyRKBMcYEq4MH4YILoKR/G3haIjDGmGAl4voQ+JklAmOMCVbPPQebN/v9NJYIjDEmxFkiMMaYYKQK3bvD9Ol+P5UlAmOMCUYHDsAXX8C2bX4/lSUCY4wJRlkzktWq5fdTWSIwxphgtGOH++mnWcmys0RgjDHBaMYMKFHC7xPXgyUCY4wJTrVrw7Bh0LCh309l8xEYY0wwevDBQjuVPREYY0ywOXgQjh8vtNNZIjDGmGBz331+n4wmO0sExhgTbHbuLJQxhrJYIjDGmGCzc6erLC4klgiMMSbYWCIwxpgQpuo6kxViIrDmo8YYE2xeeAGaNSu001kiMMaYYCICI0cW6imtaMgYY4LJ8uXw7beFeso8nwhEpJqvN6rqbwUfjjHGhLDDh6F/f6heHX78sdBO66toaC+QCmR4liXbNgUa+SsoY4wJSe+/D5s2wauvFuppfSWCV4AuwHfAVOBbVdVCicoYY0LN0aPwf/8HDRpAz56Feuo86whU9W4gFvgQuAlYISLPiUhE4YRmjDEhZO9e97N3b1dhXIh8Vhar8zVwH/AmcAvQtTACM8aYkFKyJFxzDfTpU+inzjMRiEh5EblBRGYD84AKQBtVnZDfg4tITxFZLyIpInJ/Hvv0F5FkEVkjIv856yswxpjioHZt+Phj6NKl0E/tq45gN7ARmOb5qUC8iMQDqOpMXwcWkTDgdaAbrtJ5mYjMUdXkbPs0BR4ALlXVfSLi/znZjDHGnMJX0dAHwArgYuAq4Opsr6vycex2QIqqblbV47iEck2OfW4FXlfVfQCquvvswjfGmGLiyy/d/MTLlxf6qX09Efz3TN/6zyAc2JZtORW4JMc+zQBE5DsgDHhMVT/LeSARuQ24DaB+/frnEZIxxgSp/fthzx5XV1DIfD0RPFwI5y8JNAU6AwOBCSJSJedOqjpeVeNVNb5mzZqFEJYxxhSyTZvcz0qVCv3U/hxiYjtQL9tyXc+67FKBOaqarqpbgA24xGCMMaFjxQp46CG47DLXj6CQ+XoGiRSRn3JZL7iWpS3PcOxlQFNPv4PtwADghhz7fIx7EnhHRGrgioo25ydwY4wpNh56CCpUcK2GCrkPAfhOBFtwFcPnRFUzRGQ48Dmu/H+Sqq4RkSeARFWd49nWXUSSgRPAKFVNO9dzGmNMkfT++24ymurVA3J6yWvUCBFZrqpxhRzPGcXHx2tiYmKgwzDGmCJFRJJUNT63bb7qCFr5KR5jjDFZ5s6Fp56CEycCFoKvRJBb/YAxxpiCNHUqjBsHYWEBC8FXIrCRRo0xxp8yM2HePLj88oCG4auyuKWIHMhlfVarocJv7GqMMcXJtm3w+++u2WgA+UoEq1W1daFFYowxoWbSJPez5Zla4/uXzVlsjDGBUrEi9OgB7dsHNAxfieDDQovCGGNC0d//Dp+dNrxaofOVCKqJyO05V4rI7SLyjB9jMsaY4u/AgYA2Gc3OVyLoAozPZf0E8jcMtTHGmNyouikpmzYNimTgKxGUyW2yelXNxLUcMsYYcy6WLIFFi+DuuwPafyCLr0Rw1DOD2Ck86476LyRjjCnmvv/e/bwh5zicgeGr+eijwKci8iSQ5FkXj5ta8m4/x2WMMcVTejrMmgUREVCjRqCjAXwkAlX9VESuBUYBIzyrfwb6qurqQojNGGOKn/fec08Er78e6Ei88hx9NNedRaoCv+dWd1BYbPRRY0yRd/Cg60NQiM5p9FEReVREIj2/lxGRBcAmYJeIdPVPqMYYEwIKOQmcia/K4uuB9Z7fB3v2rQl0Av7p57iMMab4WbUKLrnEzUQWRHxVFh/PVgTUA5iqqieAtSLi633GGGNyOnwYunRxU1HWrx/oaE7h64ngDxFpISI1cZ3L5mfbVs6/YRljTDEzdy7s2wfTp0NccE3+6Oub/UhgBq446EVV3QIgIr2AFYUQmzHGFA/798PIkVCpUsCHnM6Nr+ajPwKRuayfB8zzZ1DGGFOsHD8OR47AzJlQtmygozmNz7J+EWmB60fQ3LNqDTDW+hEYY8xZqFkTli+HJk0CHUmufDUfvQaYBSwChnhei4CZnm3GGGPO5OOPITExaJMA+H4ieALopqpbs637ydOfYLbnZYwxJi/79sEtt8DVV8O//x3oaPLkq9VQyRxJAADPulL+CsgYY4qNt95ycxL/7W+BjsQnX4kgQ0ROa+wqIg2ADP+FZIwxxcTixdC8OcTGBjoSn3wVDY0BvhSRf3Lq6KP3A6P9HZgxxhR569ZBu3aBjuKMfDUf/VhEtgD3cnL00TVAf1VdVRjBGWNMkaUKbdtCdHSgIzmjPBOBiJT0fOAPKsR4jDGmeBBxFcVF4InAVx3B0qxfROTVQojFGGOKh9at4a67oEcPqFo10NGcka9EkH1e4kv9HYgxxhQbW7cGOoKz4quyOGCTzxhjTJGkCi+/7JqMBsk0lPnhKxFEishPuCeDxp7f8Syrqrb0e3TGGFOUjBjhpqBs1AiuvTbQ0eSbr0QQVWhRGGNMcdC7N5QvD089BSWLzrQtvpqP/q8wAzHGmCLp0CFYsQI6doTu3d2riPFVWXzeRKSniKwXkRQRud/Hfn1FREUk14mVjTEmKKWkQGQk/OlPMG1aoKM5Z35LBCISBrwOXAlEAwNF5LSeFSJSETcJzo/+isUYYwpcSgpcein89hvMmgX9+wc6onPmzyeCdkCKqm5W1ePANCC34av/ATwLHPNjLMYYU3CefRaaNYO0NDfZzLXXQgm/FrD4la+exavx0YQ0H62GwoFt2ZZTgUtynCMOqKeqn4jIKB+x3AbcBlA/yCZ9NsaEoOHDIT0dBg+GevUCHc1581WtfZXn552en+95fv61IE4sIiWAfwE3n2lfVR0PjAeIj4+3/g3GmMBYuND1FG7VCh5+ONDRFJgzthoSkW6q2jrbpvtFZDluFFJftgPZU2Vdz7osFYEWwEIRAbgImCMivVU1Mf+XYIwxfnbsGHz0kXsCaNzYjSoqcub3FRH5KdQSEbk020KHfL5vGdBURCJEpDQwAJiTtVFV96tqDVVtqKoNgR8ASwLGmOCiCtddBzfeCOXKwbvvFqskAGeYvN5jKDBJRCp7ln/HzV/sk6pmiMhw4HMgDJikqmtE5AkgUVXn+D6CMcYEgY8/hk8+gSefhHvvhbJlAx1RgRPV/BW5ZyUCVd3v14jOID4+XhMT7aHBGFNINm2C66+HH34oUr2FcxKRJFXNta/WGYt4RKSWiEwEpqnqfhGJFpGhBR6lMcYEiz/+cMNE7Nzp6gQWLCjSSeBM8lPWPxlXvFPHs7wBuNtP8RhjTOB99plrFbTKMxljpUqBjcfP8pMIaqjqB0AmuLJ/4IRfozLGmEA5fBjuuAOqVYMrrgh0NIUiP886h0WkOp7OZSKSAAS0nsAYY/zm1ltdkdALL0CpUoGOplDkJxH8Ddfss7GIfAfUBPr5NSpjjAmEpUth6lQ3btDf/hboaApNfhLBGqATcDFuUpr1+HnUUmOMKXR//AHx8TB7tptzOITk5wN9iapmqOoaVf1ZVdOBJf4OzBhjCsWePfD//h8MHeoGjuvdu1iMH3Q2fA06dxFu4LgLRKQ1JyezrwSUK4TYjDHG/0aPhnfeKVZjB50tX0VDPXADwtUFXuBkIjgAPOjfsIwxphCsXu2GjBgxAv7xj0BHEzC+Bp17V0TeAwaq6pRCjMkYY/zv4EH3NCACY8YEOpqA8llHoKqZwD2FFIsxxhSe0qVd/cDYsVC9eqCjCaj8tBr6UkT+DkwHDmetVNXf/BaVMcb4w759blrJ3r2hRg1YsqRYDx2RX/lpNXQ9bnKab4Akz8tGfTPGFB2qbnL5evVc66Bnn3XrLQkA+XgiUNWIwgjEGGP8ZtgwePNN6NABHnkEunULdERB5YyJQERKAXcAf/KsWgi85elPYIwxwemPP6BMGTh+3CWBG26AiROL5XwC5ys/RUNvAG2AcZ5XG886Y4wJLqqQmAh9+sDll7t1pUvDgw/C669bEshDfgrI2qpqq2zLC0Rklb8CMsaYc/Lrr3DllbBypfvwv+celxhE3NwCJk/5SQQnRKSxqm4CEJFG2DDUxphgkp7uksDatfDaa/DXv0KVKoGOqsjITyIYBXwtIptxvYsbALf4NSpjjMmPzEw3PlCpUhAVBQ895CaaN2clP62GvhKRprjRRwHWq+of/g3LGGN82LPHTST/6acwbhz06wdTprhiIHPW8qwsFpG2noHn8HzwxwL/AJ4XkWqFE54xxmSzYAFccglceCG89x707AnR0W6bJYFz5qvV0FvAcQAR+RPwDPBv3Oxk4/0fmjHG4GYLS/e0Vv/nPyElxVX+Ll3qkkHz5oGNrxjwVTQUlm0YieuB8ar6EfCRiKz0e2TGmNB15Aj88AO8/z7MmAEff+yag06Y4MYFKuaTyRc2n4lAREp6Jqu/Argtn+8zxphz1707fP01ZGRA+fLQqxfUr++2RdhAB/7g6wN9KrBIRPYCR4HFACLSBJu83hhTUFRh0SK47DI39k+pUq4iuHVrNxRENauS9Ddf8xE8JSJfAbWB+aqqnk0lgBGFEZwxphjLyIDvvoMbb4TUVHjgAVcH8MkngY4s5Pgs4lHVH3JZt8F/4Rhjiq3jx903/hIlXNn/sGFucpgyZeCFF2DkyEBHGLLyM9aQMcacm19+cRO/tGnjPvC3bXPr//c/VxcwcSJs2AB/+xuEhQU21hBmlb7GmIK3e7cb7vk//4FDh1wl74gRrvIX4L77XF2ACQqWCIwxBSctDSpWdENAZ2RA166uyKdTp1M7fFkSCCqWCIwx52bfPlfZm5LiXsuWude337oJYCZODHSEJp8sERhjzkwV9u6Fb75xQzpERcH69XD11W57xYrQsiWMGXOyzb8pMiwRGGNyl5EBr7wC8+fD4sWuty/A8OHw6qsQE+N6/zZq5CaCt7F+iixLBMYY55dfYNYsOHoU7r/fNfV87TU3ycuQIa7CNyHBtQACV/F7ySWBjdkUCEsExoSyI0fciJ5Tp7oWPgAdO8Lo0e4b/pIlUKtWYGM0fufXfgQi0lNE1otIiojcn8v2v4lIsoj8JCJfiUgDf8ZjTMhLT3dTOf76q1t+7TVXzj91qpvaccMGVw+QVcxjSSAk+O2JQETCgNeBbkAqsExE5qhqcrbdVgDxqnpERO4AnsONdGqMKQjp6TBzJnz2GWzeDMuXu3b9H30Ef/mLm9IxLg7atoXKlQMdrQkQfxYNtQNSVHUzgIhMA64BvIlAVb/Otv8PwI1+jMeY4ksVNm503+ZXr4YGDVxv3cOHYcAAN3RzdDTcdBO0b3+ybD883L1MSPNnIggHtmVbTgV81SwNBT7NbYOI3IZnGOz61jTNhLrjx+G33+Cii9zygAGwcCHs2uWWS5Rwo3eCG7f/u+/cB78N4WDyEBSVxSJyIxAPdMptu6qOxzMrWnx8vOa2jzHF1q+/wk8/uRm5FixwHbZq1YJVq9wQzW3buuTQrZvryduo0ckP/RIlXOcuY3zwZyLYDtTLtlzXs+4UItIVeAjo5Jkb2ZjQ9ccfbkjmlBRXhFOpEkye7IZoBtd2/667IDYWMjPdunvvPfkEYMw58GciWAY0FZEIXAIYANyQfQcRaY2bG7mnqu72YyzGBJeDB92Hfo0a7oN/0CD3zX/TJvftHuCxx1xP3X79oF07137fKnSNH/gtEahqhogMBz4HwoBJqrpGRJ4AElV1DvA8UAH4UFxztV9Utbe/YjImYLZvd8U7n3/uWuykpp7soXvhha51T3Q0/PnPbjL2evXcDF0AjRu7lzF+IicnHisa4uPjNTExMdBhGOPbrl2uvX7nzu6bf8uWbgx+cOPwd+zotl12WQCDNKFERJJUNT63bUFRWWxMkbd5M7z1FnzxBezYcbIFz+bNbmiGhQvdU0DTptZJywQdSwTGnAtVNy5PZKQr0lm7Fv71L/dNv21b9+HfocPJJp4NG7qXMUHIEoEx+fXNNzBnjvvQX7zYVfg+8gg88QR06eKGabbKXFMEWSIwJktiohtj/9df3Yf6oUNw7BhMmOC23323K/ePiYGePd1r0CC3rVy5QEVtzHmzRGBCh6qrsP3xR1dmv2kTJCe74ZdLlIBx4+Cdd9y+pUpBhQpw8cUn3z9limvhU716QMI3xl8sEZjiJy3NfbPftctV0F53HdSuDU8+CY8+6vapVMmV73fp4ppulikDjz/uxuGvUcP12M0pKqpwr8OYQmKJwBRd+/a5b/aRke5D+oMPXAes9evdt/8sMTEuEVx3nfs2f8kl0KqVm3glu3r1MCYUWSIwwU3V9bQtU8ZNovLXv8KWLe6bflqa2+ell1wi6NbNDbk8cKBruXPhhW7+3KpV3X5RUfat3phcWCIwweXIEdf2fsECN6jaf//reth+/jmULeuSQHi4mzKxcWNo0eJkp6yqVWHatMDGb0wRZInABE5GBmzb5j74L7/czYo1bBi8+67bXquW64A1cKBbLlHCtdoxxhQoSwSmcCUmurL8n3925ftHj7r1a9e6sv4774QePVwRTmxsICM1JmRYIjAFKzMTdu92A6x9+61rqrltGzz/vBtQbdcuePll1/N26FBXadus2clZstq2dS9jTKGxRGDO37p1rljn4ovdKJtZs8iVKOFa7ERFufJ9cJ2wjh5124wxQcESgTl7s2e7cv116+Drr91cuQMHwn/+AzVrwuuvu4rcyy6D8uVPfa9Nl2hM0LFEYE6n6opwvv3WVc4uWuQ+wBcudNufftoV+VSt6lrv3HGHK/YB981/2LBARW6MOQeWCIyTkuJGxyxZEq691g2uBq7Ip00buOKKk/t+9BFccIFLBG5CIWNMEWaJIJSdOAHz57tB1ebMcRW80dHw7LOuOWdUFPzpTyfL97NkVewaY4oFSwSh4MQJ1zmrdGnXAWv9evdBv3u3a8tfs6ab/Dxr7PzISPcyxoQESwTF2dq1rqnm1Klw4ID7sB871jXd7N7dffDHxcE117gkYYwJSZYIioudO91kKf37u+XRo+G559zv114LfftC165uuXTpk8MtG2NCniWCokoV1qyBr75yH+qrVrn1V13lJklp2NANuXzjjW6YBmOMyYMlgqLit98gKQkaNHA9cZcudU03wbXZf+QR6N37ZMXuHXcELlZjTJFiiSDYnDjh2uxnZLj2+Nu2uSGXU1LctImjRrkin9q14dVXoVcvaNQo0FEbY4owSwSBpgorVsCHH7qhl2vWhLlzXXv+xES3T5Mm0KmTK+tv3dqtq18fhg8PWNjGmOLDEkEgTZzopk/cutUtx8W5qROzLF8ekLCMMaHFRv4qDJmZ8N13MGKEK9/PGlO/RAnXgWviRNe2PynJNfE0xphCZE8E/rR7tyvPnz7dlfOXLeumUyxTxm2/5Rb3MsaYALJEUBBU3Qf96tWwZIkr57/rLihVCt56Czp3hmeeca16KlYMdLTGGHMKSwTnIj3dfcgDtG/vKnUzMtxyWJhr2QNuULZdu1y7fmOMCVKWCM7Gjh3w4otuIvUlS9xY+716uQreOnVcZW90NFSpcvI9lgSMMUHOEsGZHD3qJlxZtAimTHEVv/36nXwCeOSRwMZnjDHnyVoN5SZreOYsI0bAF1+4Gbe++MJNvl65cuDiM8aYAmRPBNllZrqxe557Dn75xc3CVaWKG8Wzfn2bhMUYUyzZE0GWFSvcROvdu7vfr78eKlVy2xo0sCRgjCm2QveJYO9e19qnZ0+33Lcv7N8Pb74JN998sq2/McYUc359IhCRniKyXkRSROT+XLaXEZHpnu0/ikhDf8YDuEHcHn/cFfXccMPJ9c8+60b0vP12SwLGmJDit0QgImHA68CVQDQwUESic+w2FNinqk2AF4Fn/RUPJ064D/voaHjsMejRAz799OT2fv3ccM7GGBNi/PlE0A5IUdXNqnocmAZck2Ofa4B3Pb/PAK4Q8VNh/OrV8OCDbkL2NWtg1iy45BK/nMoYY4oSf9YRhAPbsi2nAjk/eb37qGqGiOwHqgN7s+8kIrcBtwHUr1//3KKJjXW9fKtVc4O9GWOMAYpIqyFVHa+q8aoaX7NmzXM/UI0algSMMSYHf34qbgfqZVuu61mX6z4iUhKoDKT5MSZjjDE5+DMRLAOaikiEiJQGBgBzcuwzBxjs+f06YIGqqh9jMsYYk4Pf6gg8Zf7Dgc+BMGCSqq4RkSeARFWdA0wE3hORFOA3XLIwxhhTiPzaoUxV5wHzcqx7NNvvx4B+/ozBGGOMb1ZzaowxIc4SgTHGhDhLBMYYE+IsERhjTIiTotZaU0T2AP87x7fXIEev5WIqVK4TQuda7TqLl0BcZwNVzbVHbpFLBOdDRBJVNT7QcfhbqFwnhM612nUWL8F2nVY0ZIwxIc4SgTHGhLhQSwTjAx1AIQmV64TQuVa7zuIlqK4zpOoIjDHGnC7UngiMMcbkYInAGGNCXMgkAhHpKSLrRSRFRO4PdDznQ0TqicjXIpIsImtEZKRnfTUR+UJENnp+VvWsFxF5xXPtP4lIXGCv4OyISJiIrBCRuZ7lCBH50XM90z3DnCMiZTzLKZ7tDQMa+FkQkSoiMkNE1onIWhFpXxzvp4jc4/k/+7OITBWRssXlforIJBHZLSI/Z1t31vdQRAZ79t8oIoNzO1dBC4lEICJhwOvAlUA0MFBEogMb1XnJAO5V1WggAbjTcz33A1+palPgK88yuOtu6nndBrxR+CGfl5HA2mzLzwIvqmoTYB8w1LN+KLDPs/5Fz35FxcvAZ6oaCbTCXW+xup8iEg7cBcSragvc8PQDKD73czLQM8e6s7qHIlINGIOb1rcdMCYrefiVqhb7F9Ae+Dzb8gPAA4GOqwCvbzbQDVgP1Pasqw2s9/z+FjAw2/7e/YL9hZvZ7ivgcmAuILgemSVz3lvc3BftPb+X9Owngb6GfFxjZWBLzliL2/3k5Bzl1Tz3Zy7QozjdT6Ah8PO53kNgIPBWtvWn7OevV0g8EXDyP2CWVM+6Is/zuNwa+BGopao7PZt+BWp5fi/K1/8ScB+Q6VmuDvyuqhme5ezX4r1Oz/b9nv2DXQSwB3jHUwT2toiUp5jdT1XdDowFfgF24u5PEsXvfmZ3tvcwIPc2VBJBsSQiFYCPgLtV9UD2beq+ThTptsEichWwW1WTAh2Ln5UE4oA3VLU1cJiTRQhAsbmfVYFrcImvDlCe04tSiq1gvoehkgi2A/WyLdf1rCuyRKQULglMUdWZntW7RKS2Z3ttYLdnfVG9/kuB3iKyFZiGKx56GagiIlmz62W/Fu91erZXBtIKM+BzlAqkquqPnuUZuMRQ3O5nV2CLqu5R1XRgJu4eF7f7md3Z3sOA3NtQSQTLgKae1gmlcRVUcwIc0zkTEcHN97xWVf+VbdMcIKuVwWBc3UHW+kGelgoJwP5sj6tBS1UfUNW6qtoQd88WqOpfga+B6zy75bzOrOu/zrN/UH4Dy05VfwW2icjFnlVXAMkUs/uJKxJKEJFynv/DWddZrO5nDmd7Dz8HuotIVc8TVHfPOv8KdOVKIVbi9AI2AJuAhwIdz3ley2W4R8yfgJWeVy9c+elXwEbgS6CaZ3/BtZraBKzGtdoI+HWc5TV3BuZ6fm8ELAVSgA+BMp71ZT3LKZ7tjQId91lcXyyQ6LmnHwNVi+P9BB4H1gE/A+8BZYrL/QSm4uo+0nFPeUPP5R4CQzzXnALcUhix2xATxhgT4kKlaMgYY0weLBEYY0yIs0RgjDEhzhKBMcaEOEsExhgT4iwRGHMGInKRiEwTkU0ikiQi80SkmY/9D3l+Nsw+EqUxwarkmXcxJnR5Oj7NAt5V1QGeda1wY8ZsCGRsxhQUeyIwxrcuQLqqvpm1QlVXqepiERklIss848k/7usgItJcRJaKyErP/k39Hrkx+WSJwBjfWuBGyDyFiHTHjSXfDtcruI2I/MnHcf4PeFlVY4F4XM9TY4KCFQ0Zc266e14rPMsVcInhmzz2XwI8JCJ1gZmqutH/IRqTP/ZEYIxva4A2uawX4GlVjfW8mqjqxLwOoqr/AXoDR4F5InK5f8I15uxZIjDGtwVAGRG5LWuFiLQEDgBDPHNCICLhInJhXgcRkUbAZlV9BTcCZUv/hm1M/lnRkDE+qKqKSB/gJREZDRwDtgJ3A78DS1zDIg4BN3JyvPmc+gM3iUg6bqaqf/o1cGPOgo0+aowxIc6KhowxJsRZIjDGmBBnicAYY0KcJQJjjAlxlgiMMSbEWSIwxpgQZ4nAGGNC3P8H3ZUiafpFMlsAAAAASUVORK5CYII=\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()" + ] + } + ], + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}