allele_specific_reference.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "require([\"codemirror/keymap/sublime\", \"notebook/js/cell\", \"base/js/namespace\"],\n",
       "    function(sublime_keymap, cell, IPython) {\n",
       "        cell.Cell.options_default.cm_config.keyMap = 'sublime';\n",
       "    \n",
       "        var cells = IPython.notebook.get_cells();\n",
       "        for(var cl=0; cl< cells.length ; cl++){\n",
       "            cells[cl].code_mirror.setOption('keyMap', 'sublime');\n",
       "        }\n",
       "    }\n",
       ");\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "require([\"codemirror/keymap/sublime\", \"notebook/js/cell\", \"base/js/namespace\"],\n",
    "    function(sublime_keymap, cell, IPython) {\n",
    "        cell.Cell.options_default.cm_config.keyMap = 'sublime';\n",
    "    \n",
    "        var cells = IPython.notebook.get_cells();\n",
    "        for(var cl=0; cl< cells.length ; cl++){\n",
    "            cells[cl].code_mirror.setOption('keyMap', 'sublime');\n",
    "        }\n",
    "    }\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# i had to install this beauty https://github.com/matplotlib/ipympl\n",
    "# to make following to work ...\n",
    "%matplotlib widget\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# %matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "import os\n",
    "import re\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import cooler\n",
    "\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bioframe\n",
    "from bioframe.io import resources\n",
    "\n",
    "# cooltools used dekkerlab-fork, develop-branch, commit-af508bba9c6577b8271b716ba27c330bead15981\n",
    "\n",
    "from cooler.tools import split, partition\n",
    "from functools import partial\n",
    "from cooltools import expected\n",
    "\n",
    "from cooltools import numutils\n",
    "from cooltools.eigdecomp import cooler_cis_eig\n",
    "\n",
    "from helper_module import fillcolor_compartment_style, to_uscs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "# READ SOME COMPARTMENTS IN ...\n",
    "# we have only 100kb externally called compartments ...\n",
    "binsize = 100_000\n",
    "res_human = f\"{int(binsize/1000)}kb\"\n",
    "\n",
    "comp_path = \"/home/venevs/DOTS_TESTING/alv-chicken/MNP/analysis/compartment\"\n",
    "comp_name = f\"poolMNP-bettercis.{res_human}.eigs.cis.vecs.txt\"\n",
    "\n",
    "\n",
    "cmp_vals50kb = None\n",
    "\n",
    "cmp = pd.read_csv(os.path.join(comp_path,comp_name),sep='\\t')\n",
    "\n",
    "chrom=\"chr1\"\n",
    "# get chromosomal compartments and interpolated version of those ...\n",
    "cmp_vals = cmp[cmp[\"chrom\"]==chrom]['E1'].values\n",
    "# # interpolated_cmp = get_proper_interpolated_compartments(clr, cmp_vals, chrom)\n",
    "# cmp_vals = cmp_vals[::10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ls ../newref-al1\n",
    "# all-allele-spe-1percent-highcis.galGal5.mapq_30.1000.mcool\n",
    "# all-allele-spe-2percent-highcis.galGal5.mapq_30.1000.mcool\n",
    "# clone16highcis.galGal5.mapq_30.1000.mcool\n",
    "# clone17highcis.galGal5.mapq_30.1000.mcool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "binsize = 50_000\n",
    "\n",
    "get_cpath = lambda p,c: os.path.join(p,c+f\"::/resolutions/{binsize}\")\n",
    "\n",
    "c1 = cooler.Cooler(get_cpath(\"../newref-al1/\",\"clone17highcis.galGal5.mapq_30.1000.mcool\"))\n",
    "c2 = cooler.Cooler(get_cpath(\"../newref-al2/\",\"clone17highcis.galGal5.mapq_30.1000.mcool\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a bunch of decisions/assignments, based on the binsize of the cooler:\n",
    "\n",
    "if binsize == 50_000:\n",
    "    cmp_vals_tmp = np.zeros(2*len(cmp_vals))\n",
    "    cmp_vals_tmp[0::2] = cmp_vals\n",
    "    cmp_vals_tmp[1::2] = cmp_vals\n",
    "    cmp_vals = cmp_vals_tmp\n",
    "    bad_bins =  np.r_[np.arange(180,189),np.arange(650,655),np.arange(1440,1446)]\n",
    "elif binsize == 100_000:\n",
    "    bad_bins =  np.r_[np.arange(90,95),np.arange(325,330)]\n",
    "    cmp_vals = cmp_vals\n",
    "elif binsize == 250_000:\n",
    "    bad_bins =  np.array([9*4,9*4+1,32*4,32*4+1])\n",
    "    if cmp_vals50kb is None:\n",
    "        cmp_vals = cmp_vals[::3] # probably wrong ...\n",
    "    else:\n",
    "        cmp_vals = cmp_vals50kb[::5].copy()\n",
    "elif binsize == 500_000:\n",
    "    bad_bins =  np.array([9*2,9*2+1,32*2,32*2+1,72*2,72*2+1])\n",
    "    cmp_vals = cmp_vals[::5]\n",
    "elif binsize == 1_000_000:\n",
    "    bad_bins =  np.array([9,32,72])\n",
    "    cmp_vals = cmp_vals[::10]\n",
    "else:\n",
    "    print(\"uh oh - dunno !\")\n",
    "#  \n",
    "# insertion site, bin-id:    \n",
    "#\n",
    "\n",
    "if binsize == 50_000:\n",
    "    ins_row = 1451 # obvious only in raw data, balaned has al1 and al2 separated !\n",
    "elif binsize == 100_000:\n",
    "    ins_row = 725 # 290*2.5\n",
    "elif binsize == 250_000:\n",
    "    ins_row = 290 # ok\n",
    "elif binsize == 500_000:\n",
    "    ins_row = 145\n",
    "elif binsize == 1_000_000:\n",
    "    ins_row = 72\n",
    "else:\n",
    "    print(\"uh oh - dunno where is the insertion!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "parm = (\"chr1\",0,74626886)\n",
    "qarm = (\"chr1\",75126886,196202544)\n",
    "centro = (\"chr1\",74626886,75126886)\n",
    "balance = True\n",
    "\n",
    "# bad_bins = None\n",
    "ignore_diags = 2\n",
    "theregion=parm\n",
    "theregion_name=\"parm\"\n",
    "\n",
    "bad_bin_region = (\"chr1\",9000000,9500000)\n",
    "\n",
    "weight_name = \"weight\"\n",
    "weight1 = weight_name + \"1\"\n",
    "weight2 = weight_name + \"2\"\n",
    "transforms = {\"balanced\": lambda p: p[\"count\"] * p[weight1] * p[weight2]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20431 20431\n",
      "83083621 83083180\n",
      "498961285 498965657\n"
     ]
    }
   ],
   "source": [
    "print(c1.info[\"nbins\"],c2.info[\"nbins\"])\n",
    "print(c1.info[\"nnz\"],c2.info[\"nnz\"])\n",
    "print(c1.info[\"sum\"],c2.info[\"sum\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert existing eigen track as a phasing track:\n",
    "bins = c1.bins()[:].copy()\n",
    "# no support for bad bins in eigdecomp ...\n",
    "cmp_phase = np.zeros(len(bins))\n",
    "cmp_phase[:len(cmp_vals)] = cmp_vals\n",
    "bins.insert(4,\"cmp\",cmp_phase)\n",
    "\n",
    "eigvals1, eigvec_table1 = cooler_cis_eig(\n",
    "    clr=c1,\n",
    "    bins=bins,\n",
    "    regions= [parm, qarm],\n",
    "    n_eigs=2,\n",
    "    bad_bins = bad_bins,\n",
    "    phasing_track_col=\"cmp\",\n",
    "    clip_percentile=99.9,\n",
    "    sort_metric=None,\n",
    ")\n",
    "\n",
    "\n",
    "\n",
    "# insert existing eigen track as a phasing track:\n",
    "bins = c2.bins()[:].copy()\n",
    "# no support for bad bins in eigdecomp ...\n",
    "cmp_phase = np.zeros(len(bins))\n",
    "cmp_phase[:len(cmp_vals)] = cmp_vals\n",
    "bins.insert(4,\"cmp\",cmp_phase)\n",
    "\n",
    "eigvals2, eigvec_table2 = cooler_cis_eig(\n",
    "    clr=c2,\n",
    "    bins=bins,\n",
    "    regions= [parm, qarm],\n",
    "    n_eigs=2,\n",
    "    bad_bins = bad_bins,\n",
    "    phasing_track_col=\"cmp\",\n",
    "    clip_percentile=99.9,\n",
    "    sort_metric=None,\n",
    ")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/venevs/miniconda3/envs/newexp/lib/python3.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "98d9de5cdd3b4727b3288d2dd1e71ba4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20431\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fc8337bdb7bc4e8698a4524179bbae7a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntRangeSlider(value=(0, 3925), continuous_update=False, description='pq_range_zoom', ma…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot it, just in case:\n",
    "f = plt.figure(figsize=(10,5))\n",
    "ax1 = f.add_subplot(3,1,1)\n",
    "ax2 = f.add_subplot(3,1,2)\n",
    "ax3 = f.add_subplot(3,1,3)\n",
    "print(len(eigvec_table1))\n",
    "\n",
    "\n",
    "\n",
    "pq_range_zoom=widgets.IntRangeSlider(\n",
    "        value=c1.extent(chrom),\n",
    "        min=0,\n",
    "        max=c1.extent(chrom)[1],\n",
    "        step=1,\n",
    "        continuous_update=False)\n",
    "@widgets.interact(pq_range_zoom=pq_range_zoom)\n",
    "def update(pq_range_zoom):\n",
    "    ax1.clear()\n",
    "    ax2.clear()\n",
    "    ax3.clear()\n",
    "    \n",
    "    ax1.set_xlim(*pq_range_zoom)\n",
    "    ax2.set_xlim(*pq_range_zoom)\n",
    "    ax3.set_xlim(*pq_range_zoom)\n",
    "\n",
    "    fillcolor_compartment_style(eigvec_table1[\"E1\"].values, ax1)\n",
    "    #     # fillcolor_compartment_style(cmp_vals, ax)\n",
    "    #     for bb in bad_bins:\n",
    "    #         ax.axvline(bb)\n",
    "    fillcolor_compartment_style(eigvec_table2[\"E1\"].values, ax2)\n",
    "#     fillcolor_compartment_style(np.log(np.abs((eigvec_table1[\"E1\"]-eigvec_table2[\"E1\"]).values)), ax3)\n",
    "    fillcolor_compartment_style((eigvec_table1[\"E1\"]-eigvec_table2[\"E1\"]).values, ax3)\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### seeming indifference between EV1's of the same dataset mapped to 2 different genomes\n",
    "\n",
    "`EV1(al1)-EV1(al2)` appear `=0` upstream of the insertion site, but i think now it's just a relative thing - it isn't `0` exactly - visible on `log`-scale ...\n",
    "\n",
    "the difference is indeed very small, but that's I guess fine and expected given the fact only small portion of the heatmap if affected by the change ...\n",
    "\n",
    "see the diagram of affected areas below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/venevs/miniconda3/envs/newexp/lib/python3.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "520761f9031849bcae15af97f64a3c5f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fdfce9180d0>"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot it, just in case:\n",
    "f1 = plt.figure(figsize=(3,3))\n",
    "ggg, = f1.add_gridspec(1,1)\n",
    "ax = f1.add_subplot(ggg)\n",
    "\n",
    "mat = np.zeros((100,100))\n",
    "\n",
    "# highlight diagonal for visual guidance:\n",
    "mat[np.diag_indices(mat.shape[0])] +=1\n",
    "\n",
    "# highlight \"affected\" areas:\n",
    "mat[:,-10:] += 1\n",
    "mat[-10:,:] += 1\n",
    "ax.imshow(mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # redefine cmp_vals:\n",
    "# cmp_vals = eigvec_table1[eigvec_table1[\"chrom\"]==chrom][\"E1\"].values\n",
    "# if binsize == 50_000:\n",
    "#    cmp_vals50kb = cmp_vals.copy()"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}