signal-vs-EV1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# %matplotlib inline\n",
    "\n",
    "# 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",
    "\n",
    "# import IPython.display as display\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "import matplotlib as mpl\n",
    "import seaborn as sns\n",
    "\n",
    "import multiprocess as mp\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.linalg import toeplitz\n",
    "# import cooltools\n",
    "import cooler\n",
    "import os\n",
    "# import re\n",
    "import time\n",
    "from functools import partial\n",
    "from cooltools.lib.numutils import LazyToeplitz\n",
    "\n",
    "# supress divide by zero warnings and log(0) ...\n",
    "np.seterr(divide='ignore',invalid='ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### here are some custom functions that we are using downstream - just for clarity ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_chunks(interval,chsize):\n",
    "    \"\"\"\n",
    "    split interval (input as tuple) into\n",
    "    several intervals of size chsize,\n",
    "    accounting for boundary condition\n",
    "    \"\"\"\n",
    "    start,end = interval\n",
    "    int_size = end - start\n",
    "    for i in range(start, end, chsize):\n",
    "        if i + chsize < end:\n",
    "            yield (i,i+chsize)\n",
    "        else:\n",
    "            yield (i,end)\n",
    "\n",
    "\n",
    "def get_extrapolated_array(vals, factor, skip_last=None):\n",
    "    \"\"\"\n",
    "    values - factor - how many times to do it...\n",
    "    basically turning an array [1,2,3,4] into:\n",
    "    [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4] when factor=4\n",
    "    \"\"\"\n",
    "    _vals = np.asarray(vals)\n",
    "    res = np.zeros(factor*len(_vals))\n",
    "    for shift in range(factor):\n",
    "        res[shift::factor] = _vals\n",
    "    # sometimes one'd need to skip several last elements:\n",
    "    if skip_last is None:\n",
    "        return res\n",
    "    else:\n",
    "        return res[:-skip_last]\n",
    "\n",
    "def get_proper_interpolated_compartments(clr, ycmp, chrom):\n",
    "    \"\"\"\n",
    "    this is just a wrapper around 'get_extrapolated_array'\n",
    "    just expand, like 100kb into 1kb array ...\n",
    "    \"\"\"\n",
    "    bin_interval = clr.extent(chrom)\n",
    "    indices = range(*bin_interval)\n",
    "    if len(indices) < factor*len(ycmp):\n",
    "        skipping = (factor*len(ycmp) - len(indices))\n",
    "        return get_extrapolated_array(ycmp, factor, skipping)\n",
    "    else:\n",
    "        return get_extrapolated_array(ycmp, factor)\n",
    "    \n",
    "    \n",
    "def fillcolor_compartment_style(ycmp, ax, bin_range=None, lw=0.5, level=.0, color_less=\"blue\", color_more=\"red\"):\n",
    "    \"\"\"\n",
    "    fill colors between - compartment style ...\n",
    "    like, yo like, listen, like, you know, like uhmmmmm,\n",
    "    color like A-like into like red-like, and B-like in\n",
    "    bluish like you know\n",
    "    you got the jigitist you know like right yo ?!\n",
    "    \"\"\"\n",
    "    if bin_range is None:\n",
    "        y = ycmp\n",
    "        y_size = len(y)\n",
    "        x = np.arange(y_size)\n",
    "        y0 = np.zeros_like(x) + level\n",
    "    else:\n",
    "        y = ycmp[slice(*bin_range)]\n",
    "        y_size = len(y)\n",
    "        x = np.arange(*bin_range)\n",
    "        y0 = np.zeros_like(x) + level\n",
    "        \n",
    "    ax.plot(x,y,color=\"black\",linewidth=lw, alpha=0.5,linestyle=\"solid\",markersize=0)\n",
    "    ax.fill_between(x,y0,y,where=(y<=level),color = color_less)\n",
    "    ax.fill_between(x,y0,y,where=(y>level),color = color_more)\n",
    "    \n",
    "    return 0\n",
    "\n",
    "\n",
    "def get_corr_slice(bin_range, region, clr, chr_expect, comps):\n",
    "    \"\"\"\n",
    "    given an interval of bin indices st_slice - first computes\n",
    "    O/E for each of these bins and correlates it with compartment\n",
    "    signal.\n",
    "    \"\"\"\n",
    "    O = clr.matrix()[ slice(*bin_range), slice(*clr.extent(region)) ]\n",
    "    if isinstance(chr_expect, LazyToeplitz):\n",
    "        local_chrom = region.split(\":\")[0]\n",
    "        chr_offset = clr.offset(local_chrom)\n",
    "        chr_bin_range = [i-chr_offset for i in bin_range]\n",
    "        chr_region = [i-chr_offset for i in clr.extent(region)]\n",
    "        E = chr_expect[slice(*chr_bin_range),slice(*chr_region)]\n",
    "    else:\n",
    "        # E = get_exp_slice_toeplitz(clr, expect, st_slice, chromosome)\n",
    "        raise NotImplemented(\"not implemented yet ...\")\n",
    "    # we get a bunch of obs/exp ...\n",
    "    oes = pd.DataFrame((O/E).T,columns=range(*bin_range))\n",
    "    return oes.corrwith(pd.Series(comps))\n",
    "\n",
    "\n",
    "def get_corr_slice_flank(bin_range, chrom, clr, chr_expect, comps, mask):\n",
    "    \"\"\"\n",
    "    given an interval of bin indices st_slice - first computes\n",
    "    O/E for each of these bins and correlates it with compartment\n",
    "    signal.\n",
    "    \"\"\"\n",
    "    O = clr.matrix()[ slice(*bin_range), slice(*clr.extent(chrom)) ]\n",
    "    chr_offset = clr.offset(chrom)\n",
    "    chr_bin_range = [i-chr_offset for i in bin_range]\n",
    "    chr_extent = [i-chr_offset for i in clr.extent(chrom)]\n",
    "    # check mask and expected - in the lazytoeplitz form:\n",
    "    if isinstance(chr_expect, LazyToeplitz) and isinstance(mask,LazyToeplitz):\n",
    "        E = chr_expect[slice(*chr_bin_range),slice(*chr_extent)]\n",
    "        # mask for expected ...\n",
    "        M = mask[slice(*chr_bin_range),slice(*chr_extent)]\n",
    "        # apply mask:\n",
    "        E[~M] = np.nan\n",
    "    elif mask is None:\n",
    "        E = chr_expect[slice(*chr_bin_range),slice(*chr_extent)]\n",
    "    else:\n",
    "        # E = get_exp_slice_toeplitz(clr, expect, st_slice, chromosome)\n",
    "        raise NotImplemented(\"not implemented yet ...\")\n",
    "    # we get a bunch of obs/exp ...\n",
    "    oes = pd.DataFrame((O/E).T,columns=range(*bin_range))\n",
    "    return oes.corrwith(pd.Series(comps))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is how we are going to run expected tomorrow ...\n",
    "# cooltools compute-expected -p 8 -o exp.out -t cis --balance --ignore-diags 2 input.cool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### pick some data to deal with - i like clone 16 2%FA, like, yo!\n",
    "\n",
    "fabulous, yo\n",
    "\n",
    "we're like picking yo a compartment track at a given resolution like\n",
    "and we're also like yo picking a matrix/cooler at the same resolution or higher ...\n",
    "\n",
    "if it's higher yo - we just like extrapolate compartment track using `factor` yo !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "binsize = 100_000\n",
    "res_human = f\"{int(binsize/1000)}kb\"\n",
    "\n",
    "comp_path = \"../MNP/analysis/compartment\"\n",
    "comp_name = f\"poolMNP-bettercis.{res_human}.eigs.cis.vecs.txt\"\n",
    "\n",
    "\n",
    "factor = 1\n",
    "binsize = 100_000\n",
    "res_human = f\"{int(binsize/1000)}kb\"\n",
    "\n",
    "# comp_path = \"MNP/analysis/compartment\"\n",
    "# comp_name = f\"poolMNP-bettercis.{res_human}.eigs.cis.vecs.txt\"\n",
    "\n",
    "cool_path = \"../MNP/data/20200303_MNP_triplicates/results/coolers_library\"\n",
    "cool_name = f\"MNP-DT40-1-3-16-2p-R1-T1__galGal5.galGal5.mapq_30.1000.mcool::/resolutions/{binsize}\"\n",
    "\n",
    "exp_path = f\"../MNP/analysis/expected/{res_human}\"\n",
    "exp_name = f\"clone16-2p.{res_human}.expected.cis.tsv\"\n",
    "\n",
    "\n",
    "clr = cooler.Cooler(os.path.join(cool_path,cool_name))\n",
    "exp = pd.read_csv(os.path.join(exp_path,exp_name),sep='\\t')\n",
    "cmp = pd.read_csv(os.path.join(comp_path,comp_name),sep='\\t')\n",
    "\n",
    "# MNP-DT40-WT1-R1-T1__galGal5.galGal5.mapq_30.1000.mcool\n",
    "# MNP-DT40-WT2-R1-T1__galGal5.galGal5.mapq_30.1000.mcool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get chromosomal compartments and interpolated version of those ...\n",
    "chrom = \"chr1\"\n",
    "bin_interval = clr.extent(chrom)\n",
    "cmp_vals = cmp[cmp[\"chrom\"]==chrom]['E1'].values\n",
    "interpolated_cmp = get_proper_interpolated_compartments(clr, cmp_vals, chrom)\n",
    "\n",
    "# lazy expected here - just in case ...\n",
    "lazy_exp = LazyToeplitz(exp[exp['chrom']==chrom][\"balanced.avg\"].values)\n",
    "# #\n",
    "# flank_bins = 350\n",
    "# flank_mask = np.zeros(len(interpolated_cmp),dtype=np.bool)\n",
    "# flank_mask[:flank_bins+1] = True\n",
    "# mask = LazyToeplitz(flank_mask)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### that's it yo!\n",
    "\n",
    "some parameters set are irrelevant for later, as we did use ipywidgets yo\n",
    "which is cool like yo\n",
    "\n",
    "so you can change parameters on the fly..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fb021b7ea8904b3ca7c1d52d60321d78",
       "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": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "423f1b784c1d4c38a390af3b7c2bf757",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=-0.8, description='vmin', max=0.1, min=-1.0), FloatSlider(value=2.0, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize compartments:\n",
    "fig0 = plt.figure(figsize=(9,4),constrained_layout=True)\n",
    "\n",
    "# nice gridspec source:\n",
    "# http://www.sc.eso.org/~bdias/pycoffee/codes/20160407/gridspec_demo.html\n",
    "spec0 = gridspec.GridSpec(ncols=2, nrows=2, height_ratios=[2,.7], width_ratios=[1,1],figure=fig0)\n",
    "# # Also make sure the margins and spacing are apropriate\n",
    "spec0.update(left=0.05, right=0.95, bottom=0.08, top=0.93, wspace=0.02, hspace=0.03)\n",
    "# # BUT: this is irrelevant for the saved image, if using bbox_inches='tight'in savefig !\n",
    "ax10 = fig0.add_subplot(spec0[0,0])\n",
    "ax20 = fig0.add_subplot(spec0[1,:])\n",
    "ax30 = fig0.add_subplot(spec0[0,1])\n",
    "# show some compartments and reconstructed checkerboard ...\n",
    "reconstructed_heatmap = np.outer(cmp_vals,cmp_vals)\n",
    "\n",
    "# get some O/E ...\n",
    "oo = clr.matrix().fetch(chrom)\n",
    "ee = toeplitz(exp[exp[\"chrom\"]==chrom][\"balanced.avg\"].values)\n",
    "oe_mat = np.log(oo/ee)\n",
    "\n",
    "@widgets.interact(vmin=(-1,0.1),vmax=(1,3,.1),level=(-1,1,.1))\n",
    "def update(vmin=-0.8,vmax=2,level=0):\n",
    "    ax10.clear()\n",
    "    ax20.clear()\n",
    "    ax30.clear()\n",
    "    # refresh plotting\n",
    "    ax3m = ax30.imshow(reconstructed_heatmap,cmap=\"YlOrBr\",vmin=vmin,vmax=vmax,aspect='equal')\n",
    "    ax1m = ax10.imshow(oe_mat,cmap=\"YlOrBr\",vmin=vmin,vmax=vmax)\n",
    "    fillcolor_compartment_style(cmp_vals, ax20, level=level);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### let's look at some individual loci ...\n",
    "\n",
    "individual loci to come here right after I'me done with some useless typing\n",
    "\n",
    "alright!\n",
    "\n",
    "individual locus interacting with the entire chromosome in cis or just with a large region from it ...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "######\n",
    "# THIS IS USED TO BE WHERE WE SET A BUNCH OF PARAMETERS FOR PLOTTING\n",
    "# NOW IT'S ALL DONE USING IPYWIDGETS \n",
    "######\n",
    "\n",
    "\n",
    "# chr_offset = clr.offset(chrom)\n",
    "# chr_extent = [i-chr_offset for i in clr.extent(chrom)]\n",
    "\n",
    "# # small locus region:\n",
    "# start, end = 54_900_000, 55_000_000\n",
    "# locus = f\"{chrom}:{start}-{end}\"\n",
    "# locus_extent = [i-chr_offset for i in clr.extent(locus)]\n",
    "\n",
    "# # larger region:\n",
    "# flank = flank_bins * binsize\n",
    "# print(flank)\n",
    "# start = start - flank if start - flank > 0 else 0\n",
    "# end = end + flank if end + flank < clr.chromsizes[chrom] else clr.chromsizes[chrom]\n",
    "# # start, end = 55_900_000, 89_500_000\n",
    "# region = f\"{chrom}:{start}-{end}\"\n",
    "# region_extent = [i-chr_offset for i in clr.extent(region)]\n",
    "\n",
    "# # locus relative to the region\n",
    "# region_offset = clr.offset(region)\n",
    "# locus_region_extent = [i - region_offset for i in clr.extent(locus)]\n",
    "\n",
    "# # genomic corrdinates - i.e. x:\n",
    "# x = np.arange(*region_extent)\n",
    "\n",
    "# print(f\"locus {locus} interacting with region {region}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d144495aca80434b92f64c3e77f16977",
       "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": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "486853dd85634abfa384ed9692279727",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='flank_bins', max=1000, min=10, step=10), IntSlider(valu…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig = plt.figure(figsize=(8,4),constrained_layout=True)\n",
    "spec = gridspec.GridSpec(ncols=1, nrows=5, height_ratios=[1,1,1,1,2],figure=fig)\n",
    "# # Also make sure the margins and spacing are apropriate\n",
    "spec.update(left=0.05, right=0.95, bottom=0.08, top=0.93, wspace=0.02, hspace=0.03)\n",
    "# # BUT: this is irrelevant for the saved image, if using bbox_inches='tight'in savefig !\n",
    "ax1 = fig.add_subplot(spec[0])\n",
    "ax2 = fig.add_subplot(spec[1])\n",
    "ax3 = fig.add_subplot(spec[2])\n",
    "ax4 = fig.add_subplot(spec[3])\n",
    "ax5 = fig.add_subplot(spec[4])\n",
    "\n",
    "@widgets.interact(flank_bins=(10,1000,10),end_add_bins = (1,100,1),chrom=list(clr.chromnames),start=(0,clr.chromsizes[chrom],binsize))\n",
    "def update(flank_bins=10,end_add_bins=1,start=54_900_000,chrom=\"chr1\"):\n",
    "    #\n",
    "    chr_offset = clr.offset(chrom)\n",
    "    chr_extent = [i-chr_offset for i in clr.extent(chrom)]\n",
    "\n",
    "    # small locus region:\n",
    "    start = start\n",
    "    end = start+end_add_bins*binsize\n",
    "    locus = f\"{chrom}:{start}-{end}\"\n",
    "    locus_extent = [i-chr_offset for i in clr.extent(locus)]\n",
    "\n",
    "    # larger region:\n",
    "    flank = flank_bins * binsize\n",
    "    #     print(flank)\n",
    "    start = start - flank if start - flank > 0 else 0\n",
    "    end = end + flank if end + flank < clr.chromsizes[chrom] else clr.chromsizes[chrom]\n",
    "    # start, end = 55_900_000, 89_500_000\n",
    "    region = f\"{chrom}:{start}-{end}\"\n",
    "    region_extent = [i-chr_offset for i in clr.extent(region)]\n",
    "\n",
    "    # locus relative to the region\n",
    "    region_offset = clr.offset(region)\n",
    "    locus_region_extent = [i - region_offset for i in clr.extent(locus)]\n",
    "\n",
    "    # genomic corrdinates - i.e. x:\n",
    "    x = np.arange(*region_extent)\n",
    "    \n",
    "    d = clr.matrix().fetch(locus,region)\n",
    "    eee = lazy_exp[slice(*locus_extent),slice(*region_extent)]\n",
    "\n",
    "    #\n",
    "    ax1.clear()\n",
    "    ax2.clear()\n",
    "    ax3.clear()\n",
    "    ax4.clear()\n",
    "    ax5.clear()\n",
    "    #\n",
    "    vmin,vmax = -9,-5\n",
    "    # extent: left, right, bottom, top\n",
    "    extent = [*region_extent,*locus_extent[::-1]]\n",
    "    ax1.imshow(np.log(d),cmap=\"YlOrBr\",vmin=vmin,vmax = vmax,extent = extent, origin=\"upper\")\n",
    "    ax2.imshow(np.log(eee),cmap=\"YlOrBr\",vmin=vmin,vmax = vmax,extent = extent)\n",
    "    ax3.imshow(np.log(d/eee),cmap=\"coolwarm\",vmin=-2,vmax = 2,extent = extent)\n",
    "    ax4.imshow(d/eee,cmap=\"coolwarm\",vmin=0,vmax = 2,extent = extent)\n",
    "    #\n",
    "    fillcolor_compartment_style(interpolated_cmp, ax5, region_extent);\n",
    "    ax5.axvspan(*locus_extent,alpha=0.3,label=region)\n",
    "    ax5.legend(loc=\"best\")\n",
    "    ax5.set_xlim(region_extent)\n",
    "    fig.suptitle(f\"{locus} interacting with {region}, log(o), log(e), log(o/e), o/e, EV1, bins {binsize} bp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3fa1c1a1a841478fbe950ee3c5d39821",
       "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": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f768a7152544b0aa9a1e2d8a83b88a4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=11, description='flank_left_bins', max=1000, min=3, step=10), IntSlider(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig2 = plt.figure(figsize=(10,4),constrained_layout=True)\n",
    "spec2 = gridspec.GridSpec(ncols=3, nrows=3, height_ratios=[1,1,1],width_ratios=[2,1,0.05],figure=fig2)\n",
    "# # Also make sure the margins and spacing are apropriate\n",
    "spec2.update(left=0.05, right=0.95, bottom=0.08, top=0.93, wspace=0.02, hspace=0.03)\n",
    "# # BUT: this is irrelevant for the saved image, if using bbox_inches='tight'in savefig !\n",
    "ax12 = fig2.add_subplot(spec2[0,0])\n",
    "ax22 = fig2.add_subplot(spec2[1,0])\n",
    "ax32 = fig2.add_subplot(spec2[2,0])\n",
    "ax42 = fig2.add_subplot(spec2[:,1])\n",
    "ax2_cmap = fig2.add_subplot(spec2[:,2])\n",
    "\n",
    "@widgets.interact(flank_left_bins=(3,1000,10),flank_right_bins=(3,1000,10),chrom=list(clr.chromnames),locus=(0,clr.chromsizes[chrom],binsize),scale=['log','linear'])\n",
    "def update_oe(flank_left_bins=11,flank_right_bins=11,locus=44_900_000,chrom=\"chr1\",scale='linear'):\n",
    "    ax12.clear()\n",
    "    ax22.clear()\n",
    "    ax32.clear()\n",
    "    ax42.clear()\n",
    "    ax2_cmap.clear()\n",
    "\n",
    "    # depends only on the chromosome\n",
    "    chr_offset = clr.offset(chrom)\n",
    "    \n",
    "    # get a single bp level focus:\n",
    "    locus_start = clr.offset((chrom,locus,None)) #locus_extent[0] + chr_offset\n",
    "    # locus name:\n",
    "    locus_name = f\"{chrom}:{locus}\"\n",
    "\n",
    "    # larger region:\n",
    "    flank_left = flank_left_bins * binsize\n",
    "    flank_right = flank_right_bins * binsize\n",
    "    start = locus - flank_left if locus - flank_left > 0 else 0\n",
    "    end = locus + flank_right if locus + flank_right < clr.chromsizes[chrom] else clr.chromsizes[chrom]\n",
    "    region = f\"{chrom}:{start}-{end}\"\n",
    "    region_extent = [i-chr_offset for i in clr.extent(region)]\n",
    "    # genomic corrdinates - i.e. x:\n",
    "    x = np.arange(*region_extent)\n",
    "\n",
    "    # extracting a 1D observed ...\n",
    "    d, = clr.matrix()[locus_start, slice(*clr.extent(region))]\n",
    "    eee, = lazy_exp[locus_start - chr_offset, slice(*region_extent)]\n",
    "    if scale == \"linear\":\n",
    "        oe = d/eee\n",
    "        oe_label = \"o/e\"\n",
    "        base_level = [1,1]\n",
    "    elif scale == \"log\":\n",
    "        oe = np.log(d/eee)\n",
    "        oe_label = \"log(o/e)\"\n",
    "        base_level = [0,0]\n",
    "    else:\n",
    "        print(\"whaaaaat?!\")\n",
    "    \n",
    "    ax12.plot(x,np.log(d),'ro-',label=\"log(o)\",alpha=0.5,markersize=3)\n",
    "    ax12.plot(x,np.log(eee),'b-',label=\"log(e)\")\n",
    "    ax12.set_xlim(region_extent)\n",
    "    ax12.legend(loc='best',frameon=False)\n",
    "\n",
    "    ax22.plot(x,oe,'go-',label=oe_label,alpha=0.5,markersize=3)\n",
    "    ax22.plot(region_extent,base_level,'k-',label=\"o=e\")\n",
    "    ax22.set_xlim(region_extent)\n",
    "    ax22.legend(loc='best',frameon=False)\n",
    "\n",
    "    fillcolor_compartment_style(interpolated_cmp, ax32, region_extent);\n",
    "    ax32.axvline(locus_start - chr_offset, label=locus)\n",
    "    ax32.set_xlim(region_extent)\n",
    "    ax32.legend(loc=\"best\",frameon=False)\n",
    "\n",
    "    inter_cmp_inrange = interpolated_cmp[slice(*region_extent)]\n",
    "    some_corrs = pd.DataFrame({\"x\":inter_cmp_inrange,\"y\":oe}).corr()\n",
    "    r1 = some_corrs[\"x\"][\"y\"]\n",
    "\n",
    "    sss = np.abs(locus_start - chr_offset - x)*binsize/1_000_000\n",
    "    scat = ax42.scatter(inter_cmp_inrange,oe,c=sss,alpha=0.7,label=f\"r={r1}\",cmap=\"cividis\")\n",
    "    ax42.set_ylabel(oe_label)\n",
    "    ax42.legend(loc=\"best\")\n",
    "    fig2.colorbar(scat,cax=ax2_cmap,label=f\"genomic distance from {locus_name}, mb\")\n",
    "    # this one below - was causing all y-axes to shrink tremendously:\n",
    "    # ax2_cmap.set_ylabel(f\"genomic distance from locus {locus}, megabases\")\n",
    "    \n",
    "    fig2.suptitle(f\"{locus_name} interacting with {region}: obs, exp, EV1, o/e and EV1 correlations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### insertion site 1; chr1:72,565,520\n",
    "#### insertion site 3; chr1: 72,536,061, Gallus_gallus-5.0/galGal5 Assembly Dec 2015\n",
    "\n",
    "\n",
    "#### Insertion site 2; chr1: 72,548,590 bp Gallus_gallus-5.0/galGal5 Assembly Dec 2015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/venevs/miniconda3/envs/cooler-env/lib/python3.6/site-packages/ipykernel_launcher.py:8: 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": "6677efc92da44f19a43b72bb26994347",
       "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": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "521cc7f9fa8e443ea84ca54e2126f1a3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=90, description='locus', max=1963), Output()), _dom_classes=('widget-int…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eee, = lazy_exp[locus_start - chr_offset, slice(*region_extent)]\n",
    "\n",
    "\n",
    "region1 = \"chr1:5000000-12500000\"\n",
    "region2 = \"chr1:65000000-74626886\"\n",
    "\n",
    "\n",
    "ff = plt.figure(figsize=(9,3),constrained_layout=True)\n",
    "ss = gridspec.GridSpec(ncols=2, nrows=1,figure=ff)\n",
    "# # Also make sure the margins and spacing are apropriate\n",
    "ss.update(left=0.05, right=0.95, bottom=0.08, top=0.93, wspace=0.02, hspace=0.03)\n",
    "# # BUT: this is irrelevant for the saved image, if using bbox_inches='tight'in savefig !\n",
    "axxl = ff.add_subplot(ss[0,0])\n",
    "axxr = ff.add_subplot(ss[0,1],sharey=axxl)\n",
    "\n",
    "@widgets.interact(locus=clr.extent(chrom))\n",
    "def ppp(locus=90):\n",
    "    d1, = clr.matrix()[locus, slice(*clr.extent(region1))]\n",
    "    d2, = clr.matrix()[locus, slice(*clr.extent(region2))]\n",
    "    x1 = range(*clr.extent(region1))\n",
    "    x2 = range(*clr.extent(region2))\n",
    "\n",
    "    axxl.clear()\n",
    "    axxl.plot(x1,np.log(d1),'ro-',alpha=0.4)\n",
    "    axxl.set_title(x1[np.nanargmax(d1)])\n",
    "    axxl.axvline(x1[np.nanargmax(d1)])\n",
    "#     axxl.set_xticklabels(range(*clr.extent(region1)))\n",
    "    axxr.clear()\n",
    "    axxr.plot(x2,np.log(d2),'ro-',alpha=0.4)\n",
    "    axxr.axvline(x2[np.nanargmax(d2)])\n",
    "    axxr.set_title(x2[np.nanargmax(d2)])\n",
    "#     axxr.set_xticklabels(range(*clr.extent(region2)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### there is an insertion at bin 90 - ~chr1:9_000_000 and deletion of that same piece from bin 722 ~chr1:72_200_000\n",
    "\n",
    "### which is just a tiny bit upstream of the insertion sites 1,2 and 3 that are located at bins 725: ~chr1:72_5XX_XXX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### now we correlate it on a massive scale ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tackling 4 chunks of 600 bins using 16 workers ...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37d34c440f8647b8b7f6175988d3d63f",
       "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": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5fe44d2b53184c8f9fca059ca6604f85",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=300, continuous_update=False, description='flank_bins', max=10000, min=1…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# chsize must be small enough to fit in the memory ...\n",
    "chunk_size = 600\n",
    "nproc = 16\n",
    "work_chunks = list(get_chunks(bin_interval,chunk_size))\n",
    "print(f\"tackling {len(work_chunks)} chunks of {chunk_size} bins using {nproc} workers ...\")\n",
    "\n",
    "\n",
    "buffer = {}\n",
    "#############\n",
    "\n",
    "figC = plt.figure(figsize=(10,4),constrained_layout=True)\n",
    "specC = gridspec.GridSpec(ncols=2, nrows=3, height_ratios=[1,1,1],width_ratios=[4,1],figure=figC)\n",
    "# # Also make sure the margins and spacing are apropriate\n",
    "specC.update(left=0.05, right=0.95, bottom=0.08, top=0.93, wspace=0.02, hspace=0.03)\n",
    "# # BUT: this is irrelevant for the saved image, if using bbox_inches='tight'in savefig !\n",
    "ax1C = figC.add_subplot(specC[0,0])\n",
    "ax2C = figC.add_subplot(specC[1,0])\n",
    "ax3C = figC.add_subplot(specC[2,0])\n",
    "\n",
    "# ax4C = figC.add_subplot(specC[0,1])\n",
    "# ax5C = figC.add_subplot(specC[1,1])\n",
    "# # ax_cmap = fig.add_subplot(spec[:,2])\n",
    "\n",
    "chrom_offset = clr.offset(chrom)\n",
    "region_extent = [ i-chrom_offset for i in clr.extent(chrom) ]\n",
    "\n",
    "x = np.arange(len(interpolated_cmp))\n",
    "y0 = np.zeros_like(interpolated_cmp)\n",
    "\n",
    "\n",
    "\n",
    "@widgets.interact(\n",
    "    flank_bins=widgets.IntSlider(\n",
    "        value=300,min=100,max=10000,step=100,continuous_update=False),\n",
    "    ylim = widgets.FloatRangeSlider(\n",
    "        value=[-0.7, 0.8],min=-1,max=1.2,step=0.01,continuous_update=False),\n",
    "    rerun = False,\n",
    "    lw = (0,1,0.1),\n",
    "    range_zoom=widgets.IntRangeSlider(\n",
    "        value=region_extent,min=0,max=region_extent[1],step=100,continuous_update=False)\n",
    ")\n",
    "def update(flank_bins, ylim, rerun, range_zoom, lw=0.1):\n",
    "    ax1C.clear()\n",
    "    ax2C.clear()\n",
    "    ax3C.clear()\n",
    "    \n",
    "    if rerun:\n",
    "        # re-define mask on the fly ...\n",
    "        flank_mask = np.zeros(len(interpolated_cmp),dtype=np.bool)\n",
    "        flank_mask[:flank_bins+1] = True\n",
    "        mask = LazyToeplitz(flank_mask)\n",
    "\n",
    "        p = mp.Pool(processes=nproc)\n",
    "\n",
    "        t1 = time.time()\n",
    "        # get corrs flanked ...\n",
    "        job = partial(\n",
    "                get_corr_slice_flank,\n",
    "                chrom=chrom,\n",
    "                clr=clr,\n",
    "                chr_expect=lazy_exp,\n",
    "                comps=interpolated_cmp,\n",
    "                mask = mask\n",
    "               )\n",
    "        corrs = p.map(job, work_chunks)#, chunksize=chunks)\n",
    "        t2 = time.time()\n",
    "        # print(f\"done corrs in {t2-t1} sec\")\n",
    "        display(f\"done corrs in {t2-t1} sec\")\n",
    "\n",
    "\n",
    "        t1 = time.time()\n",
    "        # get corrs gw ...\n",
    "        job = partial(\n",
    "                get_corr_slice,\n",
    "                region=chrom,\n",
    "                clr=clr,\n",
    "                chr_expect=lazy_exp,\n",
    "                comps=interpolated_cmp,\n",
    "               )\n",
    "        corrs2 = p.map(job, work_chunks)#, chunksize=chunks)\n",
    "        t2 = time.time()\n",
    "        display(f\"done corrs2 in {t2-t1} sec\")\n",
    "\n",
    "        p.close()\n",
    "        p.terminate()\n",
    "\n",
    "        corrs = pd.concat(corrs)\n",
    "        delta=pd.concat(corrs2)\n",
    "        \n",
    "        buffer[\"corrs\"] = corrs\n",
    "        buffer[\"delta\"] = delta\n",
    "\n",
    "    # ax1\n",
    "    fillcolor_compartment_style(interpolated_cmp, ax1C, range_zoom, lw);\n",
    "    #ax1C.legend(loc=\"best\")\n",
    "    ax1C.set_xlim(range_zoom)\n",
    "    ax1C.set_ylabel(\"EV1\")\n",
    "\n",
    "    # ax2\n",
    "    fillcolor_compartment_style(buffer[\"corrs\"], ax2C, range_zoom, lw);\n",
    "    #ax2C.legend(loc=\"best\")\n",
    "    ax2C.set_xlim(range_zoom)\n",
    "    ax2C.set_ylabel(\"pearson-R-flank\")\n",
    "    ax2C.set_ylim(*ylim)\n",
    "\n",
    "    # ax3\n",
    "    fillcolor_compartment_style(buffer[\"delta\"], ax3C, range_zoom, lw);\n",
    "    #ax3C.legend(loc=\"best\")\n",
    "    ax3C.set_xlim(range_zoom)\n",
    "    ax3C.set_ylabel(\"pearson-R-gw\")\n",
    "    ax3C.set_ylim(*ylim)\n",
    "\n",
    "    # ax4.scatter(corrs,interpolated_cmp,alpha=0.2)\n",
    "    # ax4.axvline(0,color='red')\n",
    "    # ax4.axhline(0,color='red')\n",
    "    # ax5.scatter(delta,interpolated_cmp,alpha=0.2)\n",
    "    # ax5.axvline(0,color='red')\n",
    "    # ax5.axhline(0,color='red')\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "# t1 = time.time()\n",
    "# #get deltas ...\n",
    "# job = partial(\n",
    "#         get_delta_slice,\n",
    "#         chromosome=region,\n",
    "#         clr=clr,\n",
    "#         chr_expect=lazy_exp,\n",
    "#         comps=inter_cmp_inrange \n",
    "#        )\n",
    "# delta = p.map(job, work_chunks)\n",
    "# t2 = time.time()\n",
    "# print(f\"done deltas in {t2-t1} sec\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## legacy and records ...\n",
    "\n",
    "### figure out flanking mask here - just for records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8756, 8898)\n",
      "[False False False False False False]\n",
      "24\n",
      "17\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1669f76cff604114abdb052585144a46",
       "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": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "chrom = \"chr20\"\n",
    "local_chrom = chrom\n",
    "flank_bins = 3\n",
    "\n",
    "len_ll = len(exp[exp['chrom']==chrom][\"balanced.avg\"].values)\n",
    "M = np.zeros(len_ll,dtype=np.bool)\n",
    "M[:flank_bins+1] = True\n",
    "llll = LazyToeplitz(M)\n",
    "\n",
    "# plt.imshow(toeplitz([True,False,False,False],[True,True,True,True,True,True,True,False,False,False,False,False,False,False,False,False]))\n",
    "print(clr.extent(\"chr20\"))\n",
    "\n",
    "\n",
    "# bin_range=[8756,8756+10]\n",
    "bin_range=[8756+20,8756+26]\n",
    "# bin_range=[8898-26,clr.extent(\"chr20\")[1]]\n",
    "# bin_range = clr.extent(\"chr20\")\n",
    "chrom = \"chr20\"\n",
    "local_chrom = chrom\n",
    "chr_expect=lazy_exp\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "O = clr.matrix()[ slice(*bin_range), slice(*clr.extent(chrom)) ]\n",
    "chr_offset = clr.offset(local_chrom)\n",
    "chr_bin_range = [i-chr_offset for i in bin_range]\n",
    "chr_extent = [i-chr_offset for i in clr.extent(chrom)]\n",
    "E = chr_expect[slice(*chr_bin_range),slice(*chr_extent)]\n",
    "# mask for expected ...\n",
    "mask_rows = (np.arange(*chr_bin_range)<=flank_bins)\n",
    "start = chr_bin_range[0]-flank_bins if chr_bin_range[0]>flank_bins else 0\n",
    "end = chr_bin_range[0]+flank_bins+1 if chr_bin_range[0]+flank_bins<chr_extent[1] else None\n",
    "print(mask_rows)\n",
    "print(end)\n",
    "print(start)\n",
    "mask_cols = np.zeros(chr_extent[1],dtype=np.int)\n",
    "mask_cols[start:end] = 1\n",
    "\n",
    "\n",
    "M1 = toeplitz(mask_rows,mask_cols)\n",
    "M1 = np.array(M1,dtype=np.bool)\n",
    "\n",
    "# OO = O.copy()\n",
    "\n",
    "# OO[~M1] = np.nan\n",
    "\n",
    "\n",
    "M2 = llll[slice(*chr_bin_range),slice(*chr_extent)]\n",
    "\n",
    "for i in range(*chr_bin_range):\n",
    "    j = i\n",
    "    i -= chr_bin_range[0]\n",
    "    M1[i,j] = False\n",
    "    M2[i,j] = False\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "# # we get a bunch of obs/exp ...\n",
    "f,(ax1,ax2,ax3) = plt.subplots(nrows=3)\n",
    "# oes = pd.DataFrame((O/E).T,columns=range(*bin_range))\n",
    "f.set_size_inches(48,5)\n",
    "\n",
    "ax1.imshow(np.log(O/E))\n",
    "ax2.imshow(M1)\n",
    "ax3.imshow(M2)\n",
    "\n",
    "(M1==M2).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "## some legacy functions from previously used \"locus_exp_module\", just for records:\n",
    "\n",
    "## from locus_exp_module import *\n",
    "\n",
    "\n",
    "# def get_locus_obs(clr, locus_bins_id, region):\n",
    "#     \"\"\"\n",
    "#     given a bin_id of a locus\n",
    "#     yield a 1D array of interaction frequency\n",
    "#     with the region ...\n",
    "#     \"\"\"\n",
    "#     # unpack chrom start stop\n",
    "#     reg_start,reg_end = clr.extent(region)\n",
    "#     # it should be 1 row  - so we'll get it to avoid extracting matrix:\n",
    "#     obs_bin, = clr.matrix()[locus_bins_id, reg_start:reg_end]\n",
    "#     return obs_bin\n",
    "\n",
    "\n",
    "# def get_locus_exp(clr,clr_exp,locus_bins_id, chrom):\n",
    "#     \"\"\"\n",
    "#     given a bin_id of a locus\n",
    "#     yield a 1D array of \"expected\" interaction frequency\n",
    "#     with the region ...\n",
    "#     seems correct based on LazyToeplitz, but rather inefficient\n",
    "#     \"\"\"\n",
    "#     # extract only the relevant part of expected:\n",
    "#     chrom_exp = clr_exp[clr_exp['chrom']==chrom][[\"diag\",\"balanced.avg\"]].set_index(\"diag\")[\"balanced.avg\"]\n",
    "\n",
    "#     # our clocus exp to be filled ...\n",
    "#     locus_exp = np.zeros(len(chrom_exp))\n",
    "\n",
    "#     # our locus id will be aligned with diag=0\n",
    "#     # upstream of that we need to fill with reversed\n",
    "#     # distance decay of size: \"locus_bin - chrom_start\"\n",
    "#     # and downstream of that we'll fill with the\n",
    "#     # distance decay of size: \"chrom_end - locus_bin\":\n",
    "#     chrom_start,chrom_end = clr.extent(chrom)\n",
    "\n",
    "#     downstream = locus_bins_id - chrom_start\n",
    "#     upstream = chrom_end - locus_bins_id\n",
    "    \n",
    "#     # filling in:\n",
    "#     # reversed downstream\n",
    "#     locus_exp[0:downstream+1] = chrom_exp.loc[0:downstream][::-1].values\n",
    "#     # upstream:\n",
    "#     locus_exp[downstream+1:] = chrom_exp.loc[1:upstream-1].values\n",
    "\n",
    "#     # beware of numpy indexing and pandas .loc indexing differences\n",
    "#     # numpy.array([1,2,3,4])[:3] -> [1,2,3]\n",
    "#     # pd.Series([1,2,3,4],index=[0,1,2,3]).loc[:3] = [1,2,3,4]\n",
    "#     return locus_exp\n",
    "\n",
    "\n",
    "\n",
    "# def get_local_score(arr,idx,radius=100):\n",
    "#     \"\"\"\n",
    "#     get some local averge score around idx in arr\n",
    "#     kinda useless\n",
    "#     \"\"\"\n",
    "#     if idx<radius:\n",
    "#         return np.nanmean(arr[:idx+radius+1])\n",
    "#     elif idx > len(arr)-radius:\n",
    "#         return np.nanmean(arr[idx-radius:])\n",
    "#     else:\n",
    "#         return np.nanmean(arr[idx-radius:idx+radius+1])\n",
    "\n",
    "# def get_corr(st, chromosome, aclr, expect, comps):\n",
    "#     O, = aclr.matrix()[ st, slice(*aclr.extent(chromosome)) ]\n",
    "#     E = get_locus_exp(aclr, expect, st, chromosome)\n",
    "#     a_corr = pd.DataFrame({\"x\":O/E, \"y\":comps}).dropna().corr(method=\"pearson\")['x']['y']\n",
    "#     return float(a_corr)\n",
    "\n",
    "\n",
    "# def get_delta(st, chromosome, aclr, expect, comps):\n",
    "#     O, = aclr.matrix()[ st, slice(*aclr.extent(chromosome)) ]\n",
    "#     E = get_locus_exp(aclr, expect, st, chromosome)\n",
    "#     return float(get_delta_score(O/E,comps))\n",
    "\n",
    "\n",
    "# def get_delta_score(OE, cmp_status):\n",
    "#     \"\"\"get some local averge score around idx in arr\"\"\"\n",
    "#     return np.nanmean(OE[cmp_status>0]) - np.nanmean(OE[cmp_status<0])\n",
    "\n",
    "\n",
    "# def get_delta_slice(st_slice, chromosome, clr, expect, comps, level=.0):\n",
    "#     O = clr.matrix()[ slice(*st_slice), slice(*clr.extent(chromosome)) ]\n",
    "#     E = get_exp_slice_toeplitz(clr, expect, st_slice, chromosome)\n",
    "#     return np.nanmean((O/E)[:,comps>level],axis=1) - np.nanmean((O/E)[:,comps<=level],axis=1)\n",
    "\n",
    "# def get_exp_slice(clr, clr_exp, st_slice, chrom):\n",
    "#     \"\"\"\n",
    "#     given a slice of bin_id-s for several loci\n",
    "#     yield a 2D array of \"expected\" interaction frequency\n",
    "#     within the region ...\n",
    "#     seems correct based on LazyToeplitz, but probably inefficient\n",
    "#     because of the \"for\"-loop\n",
    "#     \"\"\"\n",
    "#     # extract only the relevant part of expected:\n",
    "#     chrom_exp = clr_exp[clr_exp['chrom']==chrom][[\"diag\",\"balanced.avg\"]].set_index(\"diag\")[\"balanced.avg\"]\n",
    "\n",
    "#     # our clocus exp to be filled ...\n",
    "#     start,end = st_slice\n",
    "#     shape = (end-start, len(chrom_exp))\n",
    "#     locus_exp = np.zeros(shape)\n",
    "\n",
    "#     # our locus id will be aligned with diag=0\n",
    "#     # upstream of that we need to fill with reversed\n",
    "#     # distance decay of size: \"locus_bin - chrom_start\"\n",
    "#     # and downstream of that we'll fill with the\n",
    "#     # distance decay of size: \"chrom_end - locus_bin\":\n",
    "#     chrom_start, chrom_end = clr.extent(chrom)\n",
    "    \n",
    "#     for locus_bins_id in range(start, end):\n",
    "#         downstream = locus_bins_id - chrom_start\n",
    "#         upstream = chrom_end - locus_bins_id\n",
    "        \n",
    "#         rel_id = locus_bins_id - start\n",
    "\n",
    "#         # filling in:\n",
    "#         # reversed downstream\n",
    "#         locus_exp[rel_id, 0:downstream+1] = chrom_exp.loc[0:downstream][::-1].values\n",
    "#         # upstream:\n",
    "#         locus_exp[rel_id, downstream+1:] = chrom_exp.loc[1:upstream-1].values\n",
    "#     return locus_exp\n",
    "\n",
    "\n",
    "# def get_exp_slice_toeplitz(clr, clr_exp, st_slice, chrom):\n",
    "#     \"\"\"\n",
    "#     given a slice of bin_id-s for several loci\n",
    "#     yield a 2D array of \"expected\" interaction frequency\n",
    "#     within the region ...\n",
    "#     seems correct based on LazyToeplitz, more efficient\n",
    "#     because toeplitz function usage ...\n",
    "#     \"\"\"\n",
    "#     # extract only the relevant part of expected:\n",
    "#     chrom_exp = clr_exp[clr_exp['chrom']==chrom][[\"diag\",\"balanced.avg\"]].set_index(\"diag\")[\"balanced.avg\"]\n",
    "\n",
    "#     # our clocus exp to be filled ...\n",
    "#     start,end = st_slice\n",
    "#     shape = (end-start, len(chrom_exp))\n",
    "#     cols_exp = np.zeros(len(chrom_exp))\n",
    "\n",
    "#     chrom_start, chrom_end = clr.extent(chrom)\n",
    "    \n",
    "#     locus_bins_id = start\n",
    "#     downstream = locus_bins_id - chrom_start\n",
    "#     upstream = chrom_end - locus_bins_id\n",
    "#     # filling in:\n",
    "#     # reversed downstream\n",
    "#     cols_exp[0:downstream+1] = chrom_exp.loc[0:downstream][::-1].values\n",
    "#     # upstream:\n",
    "#     cols_exp[downstream+1:] = chrom_exp.loc[1:upstream-1].values\n",
    "    \n",
    "#     rows_exp = chrom_exp.loc[downstream:downstream+(end-start)-1].values\n",
    "#     return toeplitz(rows_exp, cols_exp)\n",
    "\n",
    "\n",
    "# def get_corr_slice_old(st_slice, chromosome, clr, expect, comps):\n",
    "#     \"\"\"\n",
    "#     given an interval of bin indices st_slice - first computes\n",
    "#     O/E for each of these bins and correlates it with compartment\n",
    "#     signal.\n",
    "#     \"\"\"\n",
    "#     O = clr.matrix()[ slice(*st_slice), slice(*clr.extent(chromosome)) ]\n",
    "#     E = get_exp_slice_toeplitz(clr, expect, st_slice, chromosome)\n",
    "#     # we get a bunch of obs/exp ...\n",
    "#     oes = pd.DataFrame((O/E).T,columns=range(*st_slice))\n",
    "#     return oes.corrwith(pd.Series(comps))\n"
   ]
  }
 ],
 "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": 2
}