FigureSupp2-MCD-detection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "100a1b4f-b85a-48d3-9374-77358d8fd8da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # set the number of threads for many common libraries\n",
    "# from os import environ\n",
    "# N_THREADS = '1'\n",
    "# environ['OMP_NUM_THREADS'] = N_THREADS\n",
    "# environ['OPENBLAS_NUM_THREADS'] = N_THREADS\n",
    "# environ['MKL_NUM_THREADS'] = N_THREADS\n",
    "# environ['VECLIB_MAXIMUM_THREADS'] = N_THREADS\n",
    "# environ['NUMEXPR_NUM_THREADS'] = N_THREADS\n",
    "# # https://superfastpython.com/numpy-number-blas-threads/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e3e7973",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from itertools import chain\n",
    "\n",
    "# Hi-C utilities imports:\n",
    "import cooler\n",
    "import bioframe\n",
    "import cooltools\n",
    "from cooltools.lib.numutils import fill_diag\n",
    "\n",
    "# Visualization imports:\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "import matplotlib.patches as patches\n",
    "from matplotlib.ticker import EngFormatter\n",
    "\n",
    "from itertools import cycle\n",
    "\n",
    "# from ipywidgets import interact, fixed\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d743c5",
   "metadata": {},
   "source": [
    "### Import modified \"guts\" of the dotfinder submodule from `helper_func` file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e1c6ae5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from helper_func import draw_kernel\n",
    "\n",
    "# turns out still need some of the dotfinder guts in here\n",
    "from cooltools.api.dotfinder import bp_to_bins, generate_tiles_diag_band\n",
    "from cooltools.lib.numutils import LazyToeplitz\n",
    "from cooltools.lib.common import assign_regions\n",
    "\n",
    "from datashader.mpl_ext import dsshow, alpha_colormap\n",
    "import datashader as ds\n",
    "import datashader.transfer_functions as tf\n",
    "from functools import partial\n",
    "from data_catalog import bws, bws_vlim, telo_dict\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6727d9cd",
   "metadata": {},
   "source": [
    "### pick a dataset and binsize to work on ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48b49590",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# # ! pip install --upgrade --no-cache --no-deps --ignore-install cooler\n",
    "# # ls /home/dekkerlab/dots-test\n",
    "# # import higlass as hg\n",
    "# import jscatter\n",
    "import scipy\n",
    "import logging\n",
    "import multiprocess as mp\n",
    "# import mpire for nested multi-processing\n",
    "from mpire import WorkerPool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "474819ac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ! ls /data/proj_sync/data_ranger/finalcoolers_mega/ranGAP1-aux-G1s-MEGA.hg38.mapq_30.1000.mcool\n",
    "\n",
    "# 10 kb is a resolution at which one can clearly see \"dots\":\n",
    "binsize = 10_000\n",
    "# cooler files that we'll work on :\n",
    "telo_clrs = { _k: cooler.Cooler(f\"{_path}::/resolutions/{binsize}\") for _k, _path in telo_dict.items() }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76cb8d1b",
   "metadata": {},
   "source": [
    "### pick a region to work on ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "372737ae-d867-427f-8517-6cd7719b7359",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Use bioframe to fetch the genomic features from the UCSC.\n",
    "hg38_chromsizes = bioframe.fetch_chromsizes('hg38')\n",
    "hg38_cens = bioframe.fetch_centromeres('hg38')\n",
    "hg38_arms_full = bioframe.make_chromarms(hg38_chromsizes, hg38_cens)\n",
    "# # remove \"bad\" chromosomes and near-empty arms ...\n",
    "# excluded_arms = [\"chr13_p\", \"chr14_p\", \"chr15_p\", \"chr21_p\", \"chr22_p\", \"chrM_p\", \"chrY_p\", \"chrY_q\", \"chrX_p\", \"chrX_q\"]\n",
    "# hg38_arms = hg38_arms_full[~hg38_arms_full[\"name\"].isin(excluded_arms)].reset_index(drop=True)\n",
    "\n",
    "# can do 1 chromosome (or arm) as well ..\n",
    "included_arms = [\"chr1_q\", \"chr2_p\", \"chr4_q\", \"chr6_q\"]\n",
    "included_arms = hg38_arms_full[\"name\"].to_list()[:44] # all autosomal ones ...\n",
    "hg38_arms = hg38_arms_full[hg38_arms_full[\"name\"].isin(included_arms)].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a14ab6eb",
   "metadata": {},
   "source": [
    "### pre-calculate expected for the cooler ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09fe2528",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def _job(packed_data, sample):\n",
    "    # packed data -> exp_kwargs and a dict with coolers for each sample\n",
    "    exp_kwargs, clr_dict = packed_data\n",
    "    _clr = clr_dict[sample]\n",
    "    # in order to use spawn/forkserver we have to import for worker\n",
    "    from cooltools import expected_cis\n",
    "    _exp = expected_cis( _clr, **exp_kwargs)\n",
    "    return (sample, _exp)\n",
    "\n",
    "# define expected parameters in the form of kwargs-dict:\n",
    "exp_kwargs = dict(\n",
    "    view_df=hg38_arms,\n",
    "    intra_only=False,\n",
    "    nproc=12\n",
    ")\n",
    "\n",
    "# have to use daemon=False, because _job is multiprocessing-based already ...\n",
    "with WorkerPool(\n",
    "    n_jobs=8,\n",
    "    daemon=False,\n",
    "    shared_objects=( exp_kwargs, telo_clrs ),\n",
    "    start_method=\"forkserver\",  # little faster than spawn, fork is the fastest\n",
    "    use_dill=True,\n",
    ") as wpool:\n",
    "    results = wpool.map(_job, telo_clrs, progress_bar=True)\n",
    "\n",
    "# sort out the results ...\n",
    "telo_exps_cis = {sample: _exp for sample, _exp in results}\n",
    "# # old way of doing it\n",
    "# telo_exps_cis = {k: cooltools.expected_cis( _clr, **exp_kwargs) for k, _clr in telo_clrs.items()}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "760aba61",
   "metadata": {},
   "source": [
    "# generate custom kernels - similar to original from hiccups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b27dc025-4ae0-4ab8-8615-7c75bf70bfde",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# define stripy kernels for small compartment detection ...\n",
    "def get_stripy_kernels_new(halfwidth):\n",
    "    \"\"\"\n",
    "    halfwidth : int\n",
    "        half width of the kernel, kernel size must be odd number in both dimensions\n",
    "\n",
    "    returns :\n",
    "    dictionaty with kernels\n",
    "    \"\"\"\n",
    "    # kernel width defined - odd dimensions ...\n",
    "    kwidth = (2*halfwidth + 1)\n",
    "    # define str|ipe width\n",
    "    stripe_width = kwidth // 3\n",
    "\n",
    "    # create a grid of coordinates from -h to +h, to define round kernels\n",
    "    x, y = np.meshgrid(\n",
    "        np.linspace(-halfwidth, halfwidth, kwidth),\n",
    "        np.linspace(-halfwidth, halfwidth, kwidth),\n",
    "    )\n",
    "\n",
    "    # define horizontal and vertical stripes\n",
    "    maskv = ((x < stripe_width - halfwidth) | (x > halfwidth - stripe_width))\n",
    "    maskv = maskv & ((y >= stripe_width - halfwidth) & (y <= halfwidth - stripe_width))\n",
    "    maskvmid = ~maskv & ((y >= stripe_width - halfwidth) & (y <= halfwidth - stripe_width))\n",
    "    maskh = ((y < stripe_width - halfwidth) | (y > halfwidth - stripe_width))\n",
    "    maskh = maskh & ((x >= stripe_width - halfwidth) & (x <= halfwidth - stripe_width))\n",
    "    maskhmid = ~maskh & ((x >= stripe_width - halfwidth) & (x <= halfwidth - stripe_width))\n",
    "\n",
    "    # new kernels with more round donut and lowleft masks:\n",
    "    return {\n",
    "        f'mid': maskvmid,\n",
    "        f'v{halfwidth}': maskv,\n",
    "        f'h{halfwidth}': maskh,\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab8e923b-d766-4ff8-9dba-2f65f064266d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# draw_kernel??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17dbc877-fc9c-42c4-b93a-7ae40a2e2294",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# define stripy kernels of different sizes\n",
    "k4 = get_stripy_kernels_new(halfwidth=3)\n",
    "k7 = get_stripy_kernels_new(halfwidth=7)\n",
    "kl = get_stripy_kernels_new(halfwidth=10)\n",
    "\n",
    "\n",
    "# plot rounded kernels\n",
    "fig, axs = plt.subplots(ncols=len(k4), nrows=len([k4, k7, kl]), figsize=(len(k4)*2.5, len([k4, k7, kl])*2.5), squeeze=False)\n",
    "for ax_row, ks in zip(axs, [k4, k7, kl]):\n",
    "    for ax, (ktype, kernel) in zip(ax_row, ks.items()):\n",
    "        imk = draw_kernel(kernel, ax, kernel_name=ktype,cmap=\"plasma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b0e701a-1852-49c0-92ff-d1cd7cbc066a",
   "metadata": {},
   "source": [
    "## Work on a particular clr/exp pair - mostly the 5hr sample ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75443b94-5c4b-4d35-8d69-e99539e822cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "clr = telo_clrs[\"p5hR1R2\"]\n",
    "exp = telo_exps_cis[\"p5hR1R2\"]\n",
    "exp_indexed = exp.set_index([\"region1\", \"region2\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd15c333",
   "metadata": {},
   "source": [
    "# Working on a Figure here - just load the calls:\n",
    " - enriched pixels\n",
    " - clustered pixels\n",
    " - existing anchors ...\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "888c5d03-abcf-441a-b485-b9e77251dc3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# id_anchor_fnames = {\n",
    "#     \"mega_2X_enrichment\": \"ID_anchors/mega_2X_enrichment.fourth_mega.max_size.bed\",\n",
    "#     \"5hr_2X_enrichment_old\": \"ID_anchors/5hr_2X_enrichment.second_bulk.max_size.bed\",\n",
    "#     \"5hr_2X_enrichment\": \"ID_anchors/5hr_2X_enrichment.pixel_derived.bed\",\n",
    "#     \"5hr_2X_enrichment_nosing\": \"ID_anchors/5hr_2X_enrichment.pixel_derived.no_singletons.bed\",\n",
    "#     \"5hr_notinCyto_2X_enrichment_signal\": \"ID_anchors/p5notin_pCyto_anchors_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"5hr_2X_enrichment_signal\": \"ID_anchors/5hr_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"10hr_2X_enrichment_signal\": \"ID_anchors/10hrs_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"N93p5_2X_enrichment_signal\": \"ID_anchors/N93p5_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"pCyto_2X_enrichment_signal\": \"ID_anchors/pCyto_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"mCyto_2X_enrichment_signal\": \"ID_anchors/mCyto_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"mega_3X_enrichment\": \"ID_anchors/mega_3X_enrichment.fifth_mega3x.max_size.bed\",\n",
    "#     \"MEGA_2X_enrichment\": \"ID_anchors/MEGAp5_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"MEGA_weaker_2X_enrichment\": \"ID_anchors/MEGA_plus_weak_anchors_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"MEGAN93_2X_enrichment\": \"ID_anchors/MEGAN93p5_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"MEGAminus_2X_enrichment\": \"ID_anchors/MEGA_minus_ctrl_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "#     \"cyto_2x_enrichment\": \"ID_anchors/cyto_2x_enrichment.third_mCyto.max_size.bed\",\n",
    "# }\n",
    "\n",
    "# id_anchors_dict = {}\n",
    "# for id_name, fname in id_anchor_fnames.items():\n",
    "#     id_anchors_dict[id_name] = pd.read_csv(fname, sep=\"\\t\")\n",
    "#     # ...\n",
    "#     print(f\"loaded {len(id_anchors_dict[id_name]):5d} ID anchors {id_name:>20} in BED format ...\")\n",
    "\n",
    "\n",
    "# # _anchors5 = id_anchors_dict[\"MEGA_2X_enrichment\"]\n",
    "# _anchors5 = id_anchors_dict[\"5hr_2X_enrichment_signal\"]\n",
    "# _anchorsCyto = id_anchors_dict[\"pCyto_2X_enrichment_signal\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d94db032-68ad-4a5d-82b6-42ab1b27d88b",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls native_comps_10kb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f819e13b-5606-4d31-b5f9-1993b2e87a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls enriched_pixels_10kb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48111554",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# loading enriched pixels ...\n",
    "enrich_fnames = {\n",
    "    \"5hr_2X_enrichment\": \"enriched_pixels_10kb/5hr_2X_aug24.binpe\",\n",
    "}\n",
    "# let's load them all into a dictionary ...\n",
    "enriched_dict = {}\n",
    "for id_name, fname in enrich_fnames.items():\n",
    "    enriched_dict[id_name] = pd.read_csv(fname, sep=\"\\t\")\n",
    "    # ...\n",
    "    print(f\"loaded {len(enriched_dict[id_name]):5d} enriched pixels {id_name:>20} in BEDPE format ...\")\n",
    "\n",
    "\n",
    "# clustered pixels loading\n",
    "clust_fnames = {\n",
    "    \"5hr_2X_enrichment\": \"clustered_pixels_10kb/5hr_2X_aug24.binpe\",\n",
    "}\n",
    "clustered_dict = {}\n",
    "for clst_name, fname in clust_fnames.items():\n",
    "    clustered_dict[clst_name] = pd.read_csv(fname, sep=\"\\t\")\n",
    "    # ...\n",
    "    print(f\"loaded {len(clustered_dict[clst_name]):5d} clustered ID interactions {clst_name:>20} in BEDPE format ...\")\n",
    "\n",
    "\n",
    "# clustered pixels loading\n",
    "filt_fnames = {\n",
    "    \"5hr_2X_enrichment\": \"enriched_pixels_10kb/5hr_2X_second.no_singletons.binpe\",\n",
    "}\n",
    "filt_dict = {}\n",
    "for filt_name, fname in filt_fnames.items():\n",
    "    filt_dict[filt_name] = pd.read_csv(fname, sep=\"\\t\")\n",
    "    # ...\n",
    "    print(f\"loaded {len(filt_dict[filt_name]):5d} filtered (no singletons) ID interactions {filt_name:>20} in BEDPE format ...\")\n",
    "\n",
    "# loading anchors here ...\n",
    "print(\"...\")\n",
    "id_anchor_fnames = {\n",
    "    \"mega_2X_enrichment\": \"ID_anchors/mega_2X_enrichment.fourth_mega.max_size.bed\",\n",
    "    \"5hr_2X_enrichment_old\": \"ID_anchors/5hr_2X_enrichment.second_bulk.max_size.bed\",\n",
    "    \"5hr_2X_enrichment\": \"ID_anchors/5hr_2X_enrichment.pixel_derived.bed\",\n",
    "    \"5hr_2X_enrichment_nosing\": \"ID_anchors/5hr_2X_enrichment.pixel_derived.no_singletons.bed\",\n",
    "    \"5hr_notinCyto_2X_enrichment_signal\": \"ID_anchors/p5notin_pCyto_anchors_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"5hr_2X_enrichment_signal\": \"ID_anchors/5hr_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"10hr_2X_enrichment_signal\": \"ID_anchors/10hrs_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"N93p5_2X_enrichment_signal\": \"ID_anchors/N93p5_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"pCyto_2X_enrichment_signal\": \"ID_anchors/pCyto_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"mCyto_2X_enrichment_signal\": \"ID_anchors/mCyto_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"mega_3X_enrichment\": \"ID_anchors/mega_3X_enrichment.fifth_mega3x.max_size.bed\",\n",
    "    \"MEGA_2X_enrichment\": \"ID_anchors/MEGAp5_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"MEGA_weaker_2X_enrichment\": \"ID_anchors/MEGA_plus_weak_anchors_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"MEGAN93_2X_enrichment\": \"ID_anchors/MEGAN93p5_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"MEGAminus_2X_enrichment\": \"ID_anchors/MEGA_minus_ctrl_2X_enrichment.pixel_derived.signal_peaks.bed\",\n",
    "    \"cyto_2x_enrichment\": \"ID_anchors/cyto_2x_enrichment.third_mCyto.max_size.bed\",\n",
    "}\n",
    "\n",
    "id_anchors_dict = {}\n",
    "for id_name, fname in id_anchor_fnames.items():\n",
    "    id_anchors_dict[id_name] = pd.read_csv(fname, sep=\"\\t\")\n",
    "    # ...\n",
    "    print(f\"loaded {len(id_anchors_dict[id_name]):5d} ID anchors {id_name:>20} in BED format ...\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68a9a230-7021-4dbb-9a8d-da02276ff4cb",
   "metadata": {},
   "source": [
    "# Pick one condition to explore `5hr_2X_enrichment` ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d35f3901-5af8-4c44-988b-5f8800df735b",
   "metadata": {},
   "outputs": [],
   "source": [
    "_condition = \"5hr_2X_enrichment\"\n",
    "\n",
    "enriched_pixels = enriched_dict[_condition]\n",
    "clustered_pixels = clustered_dict[_condition]\n",
    "filtered_pixels = filt_dict[_condition]\n",
    "\n",
    "final_mcds = id_anchors_dict[f\"{_condition}_signal\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7deda64a-b64e-49a1-b232-42a6c1404bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "_region1 = ('chr6', 129_000_000, 129_000_000+1_500_000)\n",
    "_region2 = ('chr6', 129_000_000, 129_000_000+3_700_000)\n",
    "\n",
    "_region1_zoom = ('chr6', 129_550_000, 129_550_000 + 280_000)\n",
    "_region2_zoom = ('chr6', 130_290_000, 130_290_000 + 310_000)\n",
    "\n",
    "region1_name = bioframe.select(hg38_arms, _region1).iat[0,-1]\n",
    "region2_name = bioframe.select(hg38_arms, _region2).iat[0,-1]\n",
    "assert region1_name == region2_name\n",
    "region_name = region2_name\n",
    "\n",
    "tile_span_i = clr.extent(_region1)\n",
    "tile_span_j = clr.extent(_region2)\n",
    "_the_tile = (region_name, tile_span_i, tile_span_j )\n",
    "\n",
    "tile_start_ij = (tile_span_i[0], tile_span_j[0])\n",
    "lazy_exp = LazyToeplitz(\n",
    "    exp_indexed.loc[region_name, region_name][\"balanced.avg\"].to_numpy()\n",
    ")\n",
    "# RAW observed matrix slice:\n",
    "observed = clr.matrix()[slice(*tile_span_i), slice(*tile_span_j)]\n",
    "expected = lazy_exp[slice(*tile_span_i), slice(*tile_span_j)]\n",
    "\n",
    "# let's figure out slices' coordinates ....\n",
    "_bins_i = clr.bins()[slice(*tile_span_i)]\n",
    "_bins_j = clr.bins()[slice(*tile_span_j)]\n",
    "_chrom_i, _start_i, _end_i = _bins_i.iloc[0][\"chrom\"], _bins_i.iloc[0][\"start\"], _bins_i.iloc[-1][\"end\"]\n",
    "_chrom_j, _start_j, _end_j = _bins_j.iloc[0][\"chrom\"], _bins_j.iloc[0][\"start\"], _bins_j.iloc[-1][\"end\"]\n",
    "\n",
    "\n",
    "gOE = scipy.ndimage.gaussian_filter(\n",
    "    (observed/expected),\n",
    "    sigma=0.4,\n",
    "    order=0,\n",
    "    mode='reflect',\n",
    "    cval=0.0,\n",
    "    # radius=3,\n",
    "    truncate=1.0,\n",
    ")\n",
    "\n",
    "region1_name_zoom = bioframe.select(hg38_arms, _region1_zoom).iat[0,-1]\n",
    "region2_name_zoom = bioframe.select(hg38_arms, _region2_zoom).iat[0,-1]\n",
    "assert region1_name_zoom == region2_name_zoom\n",
    "region_name_zoom = region2_name_zoom\n",
    "\n",
    "tile_span_zoom_i = clr.extent(_region1_zoom)\n",
    "tile_span_zoom_j = clr.extent(_region2_zoom)\n",
    "_the_tile_zoom = (region_name_zoom, tile_span_zoom_i, tile_span_zoom_j )\n",
    "\n",
    "tile_start_zoom_ij = (tile_span_zoom_i[0], tile_span_zoom_j[0])\n",
    "lazy_exp = LazyToeplitz(\n",
    "    exp_indexed.loc[region_name_zoom, region_name_zoom][\"balanced.avg\"].to_numpy()\n",
    ")\n",
    "# RAW observed matrix slice:\n",
    "observed_zoom = clr.matrix()[slice(*tile_span_zoom_i), slice(*tile_span_zoom_j)]\n",
    "expected_zoom = lazy_exp[slice(*tile_span_zoom_i), slice(*tile_span_zoom_j)]\n",
    "\n",
    "# let's figure out slices' coordinates ....\n",
    "_bins_zoom_i = clr.bins()[slice(*tile_span_zoom_i)]\n",
    "_bins_zoom_j = clr.bins()[slice(*tile_span_zoom_j)]\n",
    "_chrom_zoom_i, _start_zoom_i, _end_zoom_i = _bins_zoom_i.iloc[0][\"chrom\"], _bins_zoom_i.iloc[0][\"start\"], _bins_zoom_i.iloc[-1][\"end\"]\n",
    "_chrom_zoom_j, _start_zoom_j, _end_zoom_j = _bins_zoom_j.iloc[0][\"chrom\"], _bins_zoom_j.iloc[0][\"start\"], _bins_zoom_j.iloc[-1][\"end\"]\n",
    "\n",
    "gOE_zoom = observed_zoom/expected_zoom\n",
    "\n",
    "# select MCDs from the regions here ...\n",
    "_anchors = final_mcds[[\"chrom\",\"peak_start\",\"peak_end\",\"cluster\"]].rename(columns={\"peak_start\":\"start\",\"peak_end\":\"end\"})\n",
    "_anchors_reg = bioframe.select(_anchors, _region2).reset_index(drop=True)\n",
    "_anchors_reg[\"cluster\"] = _anchors_reg[\"cluster\"] - _anchors_reg[\"cluster\"].min()\n",
    "_anchors_reg = calculate_valencies(\n",
    "    _anchors_reg,   # must be output of bedpe_to_anchors, which in turn is a clustering inside\n",
    "    clustered_pixels,\n",
    "    cluster_colname = \"cluster\",\n",
    "    valency_colname = \"valency\",\n",
    "    bed_cols = [\"chrom\", \"start\", \"end\"],\n",
    "    bedpe_cols1 = [\"chrom1\", \"start1\", \"end1\"],\n",
    "    bedpe_cols2 = [\"chrom2\", \"start2\", \"end2\"],\n",
    ")\n",
    "\n",
    "_bedpe_region = bioframe.pair_by_distance(\n",
    "    _anchors_reg,\n",
    "    min_sep=0,\n",
    "    max_sep=(_region2[2] - _region2[1])+100_000_000,\n",
    "    suffixes=(\"1\",\"2\"),\n",
    "    keep_order=True,\n",
    ")\n",
    "_bedpe_region[\"bin1_id\"] = _bedpe_region[[\"chrom1\",\"start1\",\"end1\"]].apply(clr.offset,axis=1,result_type=\"expand\")\n",
    "_bedpe_region[\"bin1_width\"] = _bedpe_region[[\"chrom1\",\"start1\",\"end1\"]].apply(clr.extent,axis=1,result_type=\"expand\").apply(np.diff,axis=1,result_type=\"expand\")[0]\n",
    "\n",
    "_bedpe_region[\"bin2_id\"] = _bedpe_region[[\"chrom2\",\"start2\",\"end2\"]].apply(clr.offset,axis=1,result_type=\"expand\")\n",
    "_bedpe_region[\"bin2_width\"] = _bedpe_region[[\"chrom2\",\"start2\",\"end2\"]].apply(clr.extent,axis=1,result_type=\"expand\").apply(np.diff,axis=1,result_type=\"expand\")[0]\n",
    "\n",
    "\n",
    "def get_bin_coverage(df):\n",
    "    \"\"\"\n",
    "    df with bin1_id and bin2_id columns\n",
    "    coverage for every bin that is out there ...\n",
    "    \"\"\"\n",
    "    # simply count \"valencies\" of enriched pixels i-s and j-s - and sum tham up togeher ...\n",
    "    b1cov = df.groupby(\"bin1_id\").size()\n",
    "    b2cov = df.groupby(\"bin2_id\").size()\n",
    "    b1cov.index.name = \"bin\"\n",
    "    b2cov.index.name = \"bin\"\n",
    "    return b1cov.add(b2cov, fill_value=0)\n",
    "\n",
    "# redefine kernel drawing for some tweaks ...\n",
    "def draw_kernel(kernel, axis=None, kernel_name=\"default\", cmap='viridis'):\n",
    "    if axis is None:\n",
    "        f, axis = plt.subplots()\n",
    "    # kernel:\n",
    "    imk = axis.imshow(\n",
    "                    kernel[::-1,::-1],  # flip it, as in convolution\n",
    "                    alpha=0.85,\n",
    "                    cmap=cmap,\n",
    "                    interpolation='nearest')\n",
    "    # draw a square around the target pixel:\n",
    "    x0 = kernel.shape[0] // 2 - 0.5\n",
    "    y0 = kernel.shape[1] // 2 - 0.5\n",
    "    rect = patches.Rectangle((x0, y0), 1, 1, lw=1, ec='r', fc='r')\n",
    "    axis.add_patch(rect)\n",
    "\n",
    "    # clean axis:\n",
    "    axis.set_xticks([])\n",
    "    axis.set_yticks([])\n",
    "    axis.set_xticklabels('',visible=False)\n",
    "    axis.set_yticklabels('',visible=False)\n",
    "    axis.set_title(f\"{kernel_name}\", fontsize=12)\n",
    "    # add a checkerboard to highlight pixels:\n",
    "    checkerboard = np.add.outer(range(kernel.shape[0]),\n",
    "                                range(kernel.shape[1])) % 2\n",
    "    # show it:\n",
    "    axis.imshow(checkerboard,\n",
    "            cmap='gray',\n",
    "            interpolation='nearest',\n",
    "            alpha=0.3)\n",
    "\n",
    "    return imk\n",
    "\n",
    "\n",
    "# more drawing functions ...\n",
    "def rectangles_around_dots(dots_bins_df, the_tile, loc=\"upper\", lw=1, ec=\"cyan\", fc=\"none\"):\n",
    "    rectangle_kwargs = dict(lw=lw, ec=ec, fc=fc)\n",
    "    # parse the tile\n",
    "    _, tspan1, tspan2 = the_tile\n",
    "    # select only visible pixels :\n",
    "    _the_dots = dots_bins_df \\\n",
    "        .query(\"\"\"(@tspan1[0] < bin1_id < @tspan1[1]) & \\\n",
    "                  (@tspan2[0] < bin2_id < @tspan2[1]) \"\"\") \\\n",
    "        .eval(\"\"\"b1 = bin1_id - @tspan1[0] - 0.5\n",
    "                 b2 = bin2_id - @tspan2[0] - 0.5 \"\"\")\n",
    "    print(f\"{len(_the_dots)} pixels are visible out of {len(dots_bins_df)} ...\")\n",
    "    # iterate over visible pixels...\n",
    "    for b1, b2 in _the_dots[[\"b1\", \"b2\"]].itertuples(index=False):\n",
    "        w1 = w2 = 1\n",
    "        if loc == \"upper\":\n",
    "            yield patches.Rectangle((b2, b1), w2, w1, **rectangle_kwargs)\n",
    "        elif loc == \"lower\":\n",
    "            yield patches.Rectangle((b1, b2), w1, w2, **rectangle_kwargs)\n",
    "        else:\n",
    "            raise ValueError(\"loc has to be uppper or lower\")\n",
    "\n",
    "# in a specific region, and exposing importnat plotting parameters\n",
    "def rectangles_around_dots_ww(dots_bins_df, the_tile, loc=\"upper\", lw=1, ec=\"cyan\", fc=\"none\", halo=30_000, ext_width=0):\n",
    "    rectangle_kwargs = dict(lw=lw, ec=ec, fc=fc)\n",
    "    # parse the tile\n",
    "    _, tspan1, tspan2 = the_tile\n",
    "    # select only visible \"boxes\" :\n",
    "    _the_dots = dots_bins_df \\\n",
    "        .query(\"\"\"(@tspan1[0] - @halo < bin1_id < @tspan1[1] + @halo) & \\\n",
    "                  (@tspan2[0] - @halo < bin2_id < @tspan2[1] + @halo) \"\"\") \\\n",
    "        .eval(\"\"\"\n",
    "                b1 = bin1_id - @tspan1[0] - @ext_width - 0.5\n",
    "                b2 = bin2_id - @tspan2[0] - @ext_width - 0.5\n",
    "                bin1_width = bin1_width + 2*@ext_width\n",
    "                bin2_width = bin2_width + 2*@ext_width\n",
    "            \"\"\")\n",
    "    print(f\"{len(_the_dots)} pixels are visible out of {len(dots_bins_df)} ...\")\n",
    "    for b1, b2, w1, w2 in _the_dots[[\"b1\", \"b2\", \"bin1_width\", \"bin2_width\"]].itertuples(index=False):\n",
    "        if loc == \"upper\":\n",
    "            yield patches.Rectangle((b2, b1), w2+1, w1+1, **rectangle_kwargs)\n",
    "        elif loc == \"lower\":\n",
    "            yield patches.Rectangle((b1, b2), w1+1, w2+1, **rectangle_kwargs)\n",
    "        else:\n",
    "            raise ValueError(\"loc has to be uppper or lower\")\n",
    "\n",
    "from scipy.signal import find_peaks\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.colors import Normalize, TwoSlopeNorm\n",
    "\n",
    "# https://stackoverflow.com/questions/48625475/python-shifted-logarithmic-colorbar-white-color-offset-to-center\n",
    "class MidPointLogNorm(LogNorm):\n",
    "    def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):\n",
    "        LogNorm.__init__(self,vmin=vmin, vmax=vmax, clip=clip)\n",
    "        self.midpoint=midpoint\n",
    "\n",
    "    def __call__(self, value, clip=None):\n",
    "        # I'm ignoring masked values and all kinds of edge cases to make a\n",
    "        # simple example...\n",
    "        vmin, midpoint, vmax = self.vmin, self.midpoint, self.vmax\n",
    "        x, y = [np.log(vmin), np.log(midpoint), np.log(vmax)], [0, 0.5, 1]\n",
    "        return np.ma.masked_array(np.interp(np.log(value), x, y))\n",
    "\n",
    "    def inverse(self, value):\n",
    "        if not self.scaled():\n",
    "            raise ValueError(\"Not invertible until scaled\")\n",
    "        # t_vmin, t_midpoint, t_vmax = np.log(self.vmin), np.log(self.midpoint), np.log(self.vmax)\n",
    "        vmin, midpoint, vmax = self.vmin, self.midpoint, self.vmax\n",
    "\n",
    "        x, y = [0, 0.5, 1], [np.log(vmin), np.log(midpoint), np.log(vmax)]\n",
    "        # # return np.ma.masked_array(np.interp(np.log(value), x, y))\n",
    "        # if np.iterable(value):\n",
    "        #     val = np.ma.asarray(value)\n",
    "        #     return np.ma.power(val, 1. / gamma) * (vmax - vmin) + vmin\n",
    "        # else:\n",
    "        # return pow(value, 1. / gamma) * (vmax - vmin) + vmin\n",
    "        return np.exp(np.interp(value, x, y))\n",
    "\n",
    "def rectangles_around_dots_pileup(matrix, dots_bins_df, the_tile, halo=30_000, half_width=10):\n",
    "    # parse the tile\n",
    "    _, tspan1, tspan2 = the_tile\n",
    "    # select only visible \"boxes\" :\n",
    "    _the_dots = dots_bins_df \\\n",
    "        .query(\"\"\"(@tspan1[0] - @halo < bin1_id < @tspan1[1] + @halo) & \\\n",
    "                  (@tspan2[0] - @halo < bin2_id < @tspan2[1] + @halo) \"\"\") \\\n",
    "        .eval(\"\"\"\n",
    "                b1 = bin1_id + bin1_width/2 - @tspan1[0] - @half_width\n",
    "                b2 = bin2_id + bin2_width/2 - @tspan2[0] - @half_width\n",
    "                bin1_width = 2*@half_width + 1\n",
    "                bin2_width = 2*@half_width + 1\n",
    "            \"\"\")\n",
    "    _acc = []\n",
    "    for b1, b2, w1, w2 in _the_dots[[\"b1\", \"b2\", \"bin1_width\", \"bin2_width\"]].itertuples(index=False):\n",
    "        # print(b1,(b1+w1),b2,(b2+w2))\n",
    "        _mat = matrix[int(b1):int(b1+w1),int(b2):int(b2+w2)]\n",
    "        # print(_mat)\n",
    "        # print(_mat.shape)\n",
    "        _acc.append(_mat)\n",
    "    return np.asarray(_acc)\n",
    "\n",
    "\n",
    "# some common parameters, kwargs and things of that nature ...\n",
    "mycmap = plt.cm.coolwarm\n",
    "mycmap.set_over(_over_color)\n",
    "mycmap.set_under(_under_color)\n",
    "\n",
    "imshow_kwargs = dict(\n",
    "    # norm=LogNorm(vmin=1/5, vmax=5),\n",
    "    norm=TwoSlopeNorm(vcenter=1, vmin=0.0025, vmax=2.5),\n",
    "    # norm=MidPointLogNorm(vmin=1/10, vmax=3, midpoint=1),\n",
    "    # cmap=\"RdBu_r\",\n",
    "    cmap=mycmap,\n",
    "    interpolation=\"antialiased\",\n",
    "    interpolation_stage=\"rgba\",\n",
    ")\n",
    "_over_color = \"crimson\"\n",
    "_under_color = \"steelblue\"\n",
    "\n",
    "suptitle_kwargs = dict(\n",
    "    x=0.05,\n",
    "    y=.975,\n",
    "    horizontalalignment='left',\n",
    "    verticalalignment='top',\n",
    "    fontsize = 12,\n",
    ")\n",
    "\n",
    "#101820FF\n",
    "_pixel_boxes_kwargs = dict(loc=\"upper\", lw=1, ec=\"black\", fc=\"black\")\n",
    "# _pixel_boxes_kwargs = dict(loc=\"upper\", lw=1, ec=\"#EDFF00FF\", fc=\"#EDFF00FF\")\n",
    "_big_boxes_kwargs = dict(loc=\"upper\", lw=1.5, ec=\"#EDFF00FF\", fc=\"none\", halo=0, ext_width=4)\n",
    "\n",
    "\n",
    "\n",
    "# Run pileups themselves for the little stackup cartoon overhere ...\n",
    "_pstack = rectangles_around_dots_pileup(\n",
    "    gOE,\n",
    "    _bedpe_region,\n",
    "    _the_tile,\n",
    "    halo=0,\n",
    "    half_width=10,\n",
    ")\n",
    "\n",
    "_code_color=\"tab:blue\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "559718fb-1e11-4e5a-aa63-7b035f8af228",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with a square Figure.\n",
    "fig = plt.figure(figsize=(7.25, 11), layout=\"none\", facecolor=\"none\")\n",
    "# make like 6 subfigure on top of each other ...\n",
    "subfigs = fig.subfigures(6, 3, hspace=0.025, wspace=0.05, width_ratios=[1,1.6,0.7])\n",
    "\n",
    "#################################################\n",
    "f1 = subfigs[0,0]\n",
    "f1.suptitle(\"$\\\\bf{a.}$ \", **suptitle_kwargs)\n",
    "ax = f1.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "ax.text(\n",
    "    .0,\n",
    "    .5,\n",
    "    \"\"\"\n",
    "observed/expected\n",
    "contact map used\n",
    "for quantification\n",
    "of pixel enrichments\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "\n",
    "# nice megabase formatiing ...\n",
    "_mb = lambda coord: f\"{coord/1_000_000:.1f}\" if coord%1_000_000 else f\"{coord//1_000_000}\"\n",
    "\n",
    "f1 = subfigs[0, 1]\n",
    "ax = f1.add_axes([0,0,1,1])\n",
    "hm = ax.imshow(gOE, **imshow_kwargs)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "_region1_title = f\"{_region1[0]}:{_mb(_region1[1])}—{_mb(_region1[2])}MB\"\n",
    "_region2_title = f\"{_region2[0]}:{_mb(_region2[1])}—{_mb(_region2[2])}MB\"\n",
    "ax.set_ylabel(_region1_title, fontsize=9, labelpad=1)\n",
    "ax.set_title(_region2_title, fontsize=9, pad=2)\n",
    "####################################################\n",
    "\n",
    "\n",
    "####################################################\n",
    "# second step ...\n",
    "f2 = subfigs[1,0]\n",
    "f2.suptitle(\"$\\\\bf{b.}$ \", **suptitle_kwargs)\n",
    "# f2.set_facecolor(\"blue\")\n",
    "k_fig = {\"M\": k7[\"mid\"], \"V\": k7[\"h7\"], \"H\": k7[\"v7\"]}\n",
    "# ax = f2.add_subplot(1, 1, 1)\n",
    "spec = f2.add_gridspec(\n",
    "    ncols=len(k_fig),\n",
    "    nrows=2,\n",
    "    height_ratios=[1,0.7],\n",
    "    wspace=0.07,\n",
    "    hspace=0.07,\n",
    ")\n",
    "# plot rounded kernels\n",
    "for ii, (ktype, kernel) in enumerate(k_fig.items()):\n",
    "    ax = f2.add_subplot(spec[0, ii])\n",
    "    imk = draw_kernel(kernel, ax, kernel_name=ktype, cmap=\"plasma\")\n",
    "# just below - write some stuff ...\n",
    "ax = f2.add_subplot(spec[1, :])\n",
    "ax.set_axis_off()\n",
    "ax.text(\n",
    "    .0,\n",
    "    .7,\n",
    "    \"\"\"\n",
    "use simple thresholding\n",
    "to detect enriched pixels\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ");\n",
    "ax.text(\n",
    "    .0,\n",
    "    .1,\n",
    "    \" M > 2*(V | H)\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\",\n",
    "    fontname='Monospace',\n",
    "    color=_code_color,\n",
    "    fontsize=12\n",
    ");\n",
    "# ...\n",
    "f2 = subfigs[1, 1]\n",
    "ax = f2.add_axes([0,0,1,1])\n",
    "hm = ax.imshow(gOE, **imshow_kwargs)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "for box in rectangles_around_dots(\n",
    "    enriched_pixels,\n",
    "    _the_tile,\n",
    "    **_pixel_boxes_kwargs,\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "# ...\n",
    "f2 = subfigs[1,2]\n",
    "ax = f2.add_axes([0,0,0.9,1])\n",
    "hm = ax.imshow(gOE_zoom, **imshow_kwargs)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "for box in rectangles_around_dots(\n",
    "    enriched_pixels,\n",
    "    _the_tile_zoom,\n",
    "    **dict(loc=\"upper\", lw=1, ec=\"black\", fc=\"none\"),\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "####################################################\n",
    "\n",
    "\n",
    "####################################################\n",
    "# third step ...\n",
    "f3 = subfigs[2,0]\n",
    "f3.suptitle(\"$\\\\bf{c.}$ \", **suptitle_kwargs)\n",
    "ax = f3.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "ax.text(\n",
    "    .0,\n",
    "    .7,\n",
    "    \"\"\"\n",
    "filter 'singletons' and\n",
    "noisy clusters using\n",
    "density based clustering\n",
    "of enriched pixels:\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "ax.text(\n",
    "    .0,\n",
    "    .25,\n",
    "    \"\"\"\n",
    "sklearn.cluster.OPTICS(\n",
    " min_samples = 5,\n",
    " max_eps = 33_000,\n",
    ").fit(enriched_pixels)\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\",\n",
    "    fontname='Monospace',\n",
    "    fontsize=9,\n",
    "    color=_code_color,\n",
    ");\n",
    "\n",
    "f3 = subfigs[2, 1]\n",
    "ax = f3.add_axes([0,0,1,1])\n",
    "hm = ax.imshow(gOE, **imshow_kwargs)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "# draw boxes around clustered pixels ...\n",
    "for box in rectangles_around_dots_ww(\n",
    "    clustered_pixels,\n",
    "    _the_tile,\n",
    "    **_big_boxes_kwargs,\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "# draw enriched pixels themselves again ...\n",
    "for box in rectangles_around_dots(\n",
    "    filtered_pixels,\n",
    "    _the_tile,\n",
    "    **_pixel_boxes_kwargs,\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "\n",
    "f3 = subfigs[2,2]\n",
    "ax = f3.add_axes([0,0,0.9,1])\n",
    "hm = ax.imshow(gOE_zoom, **imshow_kwargs)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "for box in rectangles_around_dots(\n",
    "    filtered_pixels,\n",
    "    _the_tile_zoom,\n",
    "    **dict(loc=\"upper\", lw=1, ec=\"black\", fc=\"none\"),\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "# draw boxes around clustered pixels ...\n",
    "for box in rectangles_around_dots_ww(\n",
    "    clustered_pixels,\n",
    "    _the_tile_zoom,\n",
    "    **dict(loc=\"upper\", lw=3, ec=\"#EDFF00FF\", fc=\"none\", halo=0, ext_width=2),\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "#################################################\n",
    "\n",
    "\n",
    "#################################################\n",
    "# fifth step ...\n",
    "f4 = subfigs[3,0]\n",
    "f4.suptitle(\"$\\\\bf{d.}$ \", **suptitle_kwargs)\n",
    "ax = f4.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "ax.text(\n",
    "    .0,\n",
    "    .5,\n",
    "    \"\"\"\n",
    "use 'coverage' of\n",
    "filtered pixels to\n",
    "detect MCD-anchors,\n",
    "as microcompartments\n",
    "keep contributing\n",
    "to the same anchors\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "\n",
    "f4 = subfigs[3, 1]\n",
    "ax = f4.add_axes([0,0.1,1,0.7])\n",
    "# calculate pixels coverage ...\n",
    "rich_pix_vals = get_bin_coverage(filtered_pixels)\n",
    "_from, _to = clr.extent(_region2)\n",
    "rich_pix_vals_plot = pd.Series(index=np.arange(_from,_to)).add(rich_pix_vals.loc[_from: _to], fill_value=0)\n",
    "# ax.plot(rich_pix_vals_plot.fillna(0), marker=\",\", markersize=2, lw=0.75, label=\"all enriched pixels\", color=\"black\")\n",
    "ax.plot(rich_pix_vals_plot.fillna(0), marker=\".\", markersize=1, lw=0.5, color=\"black\")\n",
    "ax.set_xlim(_from, _to)\n",
    "# coverage track params\n",
    "cov_ylim = (-3, 77)\n",
    "cov_yticks = [0,70]\n",
    "ax.set_ylim(cov_ylim)\n",
    "ax.set_yticks(cov_yticks)\n",
    "ax.set_xticks([])\n",
    "ax.set_xlabel(_region2_title, fontsize=9, labelpad=2)\n",
    "####################################################\n",
    "\n",
    "\n",
    "#################################################\n",
    "f5 = subfigs[4,0]\n",
    "f5.suptitle(\"$\\\\bf{e.}$ \", **suptitle_kwargs)\n",
    "ax = f5.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "ax.text(\n",
    "    .0,\n",
    "    .65,\n",
    "    \"\"\"\n",
    "detect coverage peaks,\n",
    "i.e. MCD-anchors, using\n",
    "simple 1D peak detection:\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "ax.text(\n",
    "    .0,\n",
    "    .2,\n",
    "    \"\"\"\n",
    "scipy.signal.find_peaks(\n",
    " coverage,\n",
    " height=7,\n",
    " distance=5,\n",
    ")\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\",\n",
    "    fontname='Monospace',\n",
    "    fontsize=9,\n",
    "    color=_code_color,\n",
    ");\n",
    "\n",
    "f5 = subfigs[4, 1]\n",
    "ax = f5.add_axes([0,0.1,1,0.7])\n",
    "# # take empty bins and fill non-zero coverage with rich_clust_pix_vals ...\n",
    "# # use value threshold to define the \"floor\" of the pixel coverage ...\n",
    "_bins = clr.bins()[:]\n",
    "_bins.index.name = \"bin\"\n",
    "_bins[\"cov\"] = rich_pix_vals\n",
    "_bins[\"cov\"] = _bins[\"cov\"].fillna(0)\n",
    "_arr = _bins[\"cov\"].to_numpy()\n",
    "\n",
    "# detect praks on the coverage track ...\n",
    "_val_thresh = 7\n",
    "_distance_thresh = 5\n",
    "_peaks, _props = find_peaks(\n",
    "    _arr,\n",
    "    height=_val_thresh,\n",
    "    prominence=(None,None),\n",
    "    distance=_distance_thresh,\n",
    ")\n",
    "# extract left/right boundaries of every peak ...\n",
    "_lefts = _props[\"left_bases\"]\n",
    "_rights = _props[\"right_bases\"]\n",
    "_arr_clipped = np.clip(_arr, _val_thresh, None)\n",
    "\n",
    "# ax.plot(_arr_clipped[_from: _to], color=\"black\", lw=0.75, marker=None)\n",
    "ax.plot(_arr[_from: _to], color=\"black\", lw=0.75, marker=None)\n",
    "ax.axhline(_val_thresh, lw=0.5, color=\"dimgray\", linestyle=\"--\")\n",
    "ax.plot(_peaks - _from, _arr_clipped[_peaks], marker=\".\", lw=0, markersize=7, color=\"orange\")\n",
    "ax.plot(_lefts - _from, _arr_clipped[_lefts], marker=8, lw=0, markersize=7, color=\"dimgray\" )\n",
    "ax.plot(_rights - _from, _arr_clipped[_rights], marker=9, lw=0, markersize=7, color=\"dimgray\" )\n",
    "# ax.plot(_lefts - _from, _arr_clipped[_lefts], marker=\"|\", lw=0, markersize=7, color=\"tab:blue\" )\n",
    "# ax.plot(_rights - _from, _arr_clipped[_rights], marker=\"|\", lw=0, markersize=7, color=\"tab:red\" )\n",
    "# ax.plot(_lefts - _from, _arr_clipped[_lefts], marker=\"|\", lw=0, markersize=12, color=\"dimgray\" )\n",
    "# ax.plot(_rights - _from, _arr_clipped[_rights], marker=\"|\", lw=0, markersize=12, color=\"dimgray\" )\n",
    "ax.set_xlim(0, _to-_from )\n",
    "ax.set_ylim(cov_ylim)\n",
    "# ax.set_yticks(cov_yticks)\n",
    "ax.set_yticks([0,_val_thresh,70])\n",
    "ax.set_xticks([])\n",
    "ax.set_xlabel(_region2_title, fontsize=9, labelpad=2)\n",
    "\n",
    "f5 = subfigs[4, 2]\n",
    "ax = f5.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "\n",
    "x_foot = 0.03\n",
    "x_summit = 0.1\n",
    "y_foot = 0.37\n",
    "y_summit = 0.63\n",
    "_radius = 0.05\n",
    "_dist = 0.17\n",
    "x_left = x_foot\n",
    "# ax.add_patch(patches.Arrow(x_foot+_radius, y_foot, _dist-2*_radius, 0, color=\"grey\", width=0.1))\n",
    "ax.scatter(x=[x_foot+0.05], y=[y_foot], s=80, marker=8, linewidths=1.5, color=\"dimgray\")\n",
    "ax.scatter(x=[x_foot+_dist-0.05], y=[y_foot], s=80, marker=9, linewidths=1.5, color=\"dimgray\")\n",
    "ax.text(\n",
    "    x=x_left+0.27,\n",
    "    y=y_foot,\n",
    "    s=\"footprint\",\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "ax.scatter(x=[x_summit], y=[y_summit], s=80, marker=\".\", linewidths=1.5, color=\"orange\")\n",
    "ax.text(\n",
    "    x=x_left+0.27,\n",
    "    y=y_summit,\n",
    "    s=\"summit\",\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "ax.set_aspect(1)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_ylim(0,1)\n",
    "####################################################\n",
    "\n",
    "\n",
    "####################################################\n",
    "# third step ...\n",
    "f6 = subfigs[5,0]\n",
    "f6.suptitle(\"$\\\\bf{f.}$ \", **suptitle_kwargs)\n",
    "ax = f6.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "ax.text(\n",
    "    .0,\n",
    "    .5,\n",
    "    \"\"\"\n",
    "explore all-by-all grid\n",
    "of detected anchors,\n",
    "i.e. microcompartment\n",
    "domains.\n",
    "\n",
    "use their 'summits' for\n",
    "centering stackups and\n",
    "pileups.\n",
    "    \"\"\",\n",
    "    transform=ax.transAxes,\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\"\n",
    ")\n",
    "\n",
    "f6 = subfigs[5, 1]\n",
    "ax = f6.add_axes([0,0,1,1])\n",
    "hm = ax.imshow(gOE, **imshow_kwargs)\n",
    "# bounds yields (x0, y0, width, height)\n",
    "x0, y0, dx, dy = ax.get_position().bounds\n",
    "# get figure bounds ...\n",
    "_, _, _h, _w = f6.figbbox.bounds\n",
    "awidth = 0.06\n",
    "aoffset = 0.02\n",
    "\n",
    "ax_x = f6.add_axes([x0, y0-awidth-aoffset, dx, awidth])\n",
    "ax_y = f6.add_axes([x0-(awidth+aoffset)*(_h/_w), y0, awidth*(_h/_w), dy])\n",
    "\n",
    "_reg2w = np.diff(clr.extent(_region2)).item()\n",
    "_reg1w = np.diff(clr.extent(_region1)).item()\n",
    "ax_x.set_xlim(0, _reg2w)\n",
    "ax_y.set_ylim(0, _reg1w)\n",
    "\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax_x.set_xticks([])\n",
    "ax_x.set_yticks([])\n",
    "ax_y.set_xticks([])\n",
    "ax_y.set_yticks([])\n",
    "\n",
    "# draw boxes around clustered pixels ...\n",
    "_big_boxes_kwargs[\"ec\"] = \"black\"\n",
    "_big_boxes_kwargs[\"lw\"] = 1\n",
    "for box in rectangles_around_dots_ww(\n",
    "    _bedpe_region,\n",
    "    _the_tile,\n",
    "    **_big_boxes_kwargs,\n",
    "):\n",
    "    ax.add_patch(box)\n",
    "    x0, y0, dx, dy = box.get_bbox().bounds\n",
    "    ax_x.axvspan(x0, x0+dx, facecolor='black', alpha=0.7)\n",
    "    ax_x.axvspan(y0, y0+dy, facecolor='black', alpha=0.7)\n",
    "    ax_y.axhspan(y0, y0+dy, facecolor='black', alpha=0.7)\n",
    "\n",
    "ax_y.invert_yaxis()\n",
    "\n",
    "f6 = subfigs[5,2]\n",
    "stack = ttt\n",
    "ooe_norm = imshow_kwargs[\"norm\"]\n",
    "# create a 21 x 21 vertex mesh\n",
    "*_, dimens = stack.shape\n",
    "X, Y = np.meshgrid(np.linspace(0,1,dimens), np.linspace(0,1,dimens))\n",
    "# create vertices for a rotated mesh (3D rotation matrix)\n",
    "Z = np.zeros_like(X)\n",
    "# show the 3D rotated projection\n",
    "ax = f6.add_subplot(1,1,1,projection='3d')\n",
    "num_to_show = 6\n",
    "for i, _matrix in enumerate( stack ):\n",
    "    if i > num_to_show:\n",
    "        break\n",
    "    is_top_element = i!=num_to_show\n",
    "    # get the heatmap for the region\n",
    "    ax.plot_wireframe(X, Y, Z+i/num_to_show, linewidths=4, rstride=1, cstride=1,color=\"black\")\n",
    "    ax.plot_surface(\n",
    "        X,\n",
    "        Y,\n",
    "        Z+i/num_to_show,\n",
    "        rstride=1,\n",
    "        cstride=1,\n",
    "        facecolors=mycmap(ooe_norm(_matrix)),\n",
    "        shade=is_top_element,\n",
    "    )\n",
    "ax.set_axis_off()\n",
    "ax.set_position([0,0,1,1])\n",
    "\n",
    "fig.savefig(\"FigureSupp2-mcd_detection.pdf\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "756b627c-a4b1-45d8-8bea-81e9b684d1c4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45100d78-4b63-4626-a8c6-809cdf283dd3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6c4ac0f-6038-4190-92b6-088011b9731a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b44f2cb-2836-461f-94c6-2ca0653b0527",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c2ac1e3-67d6-4852-85ba-a034f1e70973",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8772912-38ea-4073-903b-29404ddb559f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3357bc2-4fc3-4275-ba19-fea2e05a70ff",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac6b09d4-6edb-426b-89cd-c95e2adb998f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc0f8b39-be34-44d6-bc87-feb625acc581",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb9f01b7-1361-4edb-be53-acc4edbb0760",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c84855b5-19fc-4604-bb64-9a60fd67e784",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d19c950-68d9-421c-987d-1bc208cd906a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45673341-6ff1-4483-9341-47b4cda76352",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cd33766f-9922-4314-95be-3b849feccbaa",
   "metadata": {},
   "source": [
    "# Legacy stuff abandoned ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e06e118-2f8b-4300-aa50-cad109df8078",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rectangles_around_dots(dots_bins_df, the_tile, loc=\"upper\", lw=1, ec=\"cyan\", fc=\"none\"):\n",
    "    rectangle_kwargs = dict(lw=lw, ec=ec, fc=fc)\n",
    "    # parse the tile\n",
    "    _, tspan1, tspan2 = the_tile\n",
    "    # select only visible pixels :\n",
    "    _the_dots = dots_bins_df \\\n",
    "        .query(\"\"\"(@tspan1[0] < bin1_id < @tspan1[1]) & \\\n",
    "                  (@tspan2[0] < bin2_id < @tspan2[1]) \"\"\") \\\n",
    "        .eval(\"\"\"b1 = bin1_id - @tspan1[0]\n",
    "                 b2 = bin2_id - @tspan2[0] \"\"\")\n",
    "    print(f\"{len(_the_dots)} pixels are visible out of {len(dots_bins_df)} ...\")\n",
    "    # iterate over visible pixels...\n",
    "    for b1, b2 in _the_dots[[\"b1\", \"b2\"]].itertuples(index=False):\n",
    "        w1 = w2 = 1\n",
    "        if loc == \"upper\":\n",
    "            yield patches.Rectangle((b2, b1), w2, w1, **rectangle_kwargs)\n",
    "        elif loc == \"lower\":\n",
    "            yield patches.Rectangle((b1, b2), w1, w2, **rectangle_kwargs)\n",
    "        else:\n",
    "            raise ValueError(\"loc has to be uppper or lower\")\n",
    "\n",
    "# in a specific region, and exposing importnat plotting parameters\n",
    "def rectangles_around_dots_ww(dots_bins_df, the_tile, loc=\"upper\", lw=1, ec=\"cyan\", fc=\"none\", halo=30_000, ext_width=0):\n",
    "    rectangle_kwargs = dict(lw=lw, ec=ec, fc=fc)\n",
    "    # parse the tile\n",
    "    _, tspan1, tspan2 = the_tile\n",
    "    # select only visible \"boxes\" :\n",
    "    _the_dots = dots_bins_df \\\n",
    "        .query(\"\"\"(@tspan1[0] - @halo < bin1_id < @tspan1[1] + @halo) & \\\n",
    "                  (@tspan2[0] - @halo < bin2_id < @tspan2[1] + @halo) \"\"\") \\\n",
    "        .eval(\"\"\"\n",
    "                b1 = bin1_id - @tspan1[0] - @ext_width\n",
    "                b2 = bin2_id - @tspan2[0] - @ext_width\n",
    "                bin1_width = bin1_width + @ext_width\n",
    "                bin2_width = bin2_width + @ext_width\n",
    "            \"\"\")\n",
    "    print(f\"{len(_the_dots)} pixels are visible out of {len(dots_bins_df)} ...\")\n",
    "    for b1, b2, w1, w2 in _the_dots[[\"b1\", \"b2\", \"bin1_width\", \"bin2_width\"]].itertuples(index=False):\n",
    "        if loc == \"upper\":\n",
    "            yield patches.Rectangle((b2, b1), w2+1, w1+1, **rectangle_kwargs)\n",
    "        elif loc == \"lower\":\n",
    "            yield patches.Rectangle((b1, b2), w1+1, w2+1, **rectangle_kwargs)\n",
    "        else:\n",
    "            raise ValueError(\"loc has to be uppper or lower\")\n",
    "\n",
    "def bedpe_to_anchors(\n",
    "    bedpe_df,\n",
    "    view_df,  # for sorting !\n",
    "    cols1 = [\"chrom1\", \"start1\", \"end1\"],\n",
    "    cols2 = [\"chrom2\", \"start2\", \"end2\"],\n",
    "    mode=\"cluster\"\n",
    "):\n",
    "    \"\"\"\n",
    "    turning bedpe interactions to a bed of anchors - the simple way\n",
    "\n",
    "    mode - allow for several way to merge upstream and downstream\n",
    "    anchors. cluster, max_size, max_valency, median\n",
    "    \"\"\"\n",
    "    _cols = [\"chrom\", \"start\", \"end\"]\n",
    "    _cluster_cols = [\"chrom\", \"cluster_start\", \"cluster_end\"]\n",
    "    # concat left and right anchors ...\n",
    "    _bed = pd.concat(\n",
    "        [\n",
    "            bedpe_df[cols1].rename(columns={c1:c for c1,c in zip(cols1, _cols)}),\n",
    "            bedpe_df[cols2].rename(columns={c2:c for c2,c in zip(cols2, _cols)}),\n",
    "        ],\n",
    "        ignore_index=True,\n",
    "    )\n",
    "    # clustering anchors - define clusters of overlaping anchors ...\n",
    "    _anchors = bioframe.cluster(\n",
    "        bioframe.sort_bedframe(_bed, view_df=view_df),\n",
    "        # min_dist=None,\n",
    "        min_dist=binsize+1,\n",
    "        return_input=True,\n",
    "    ).reset_index(drop=True)\n",
    "    if mode == \"cluster\":\n",
    "        # simply return resulting clusters - i.e. total footprint of clustered anchors ...\n",
    "        _anchors = _anchors.drop_duplicates(subset=_cluster_cols).reset_index(drop=True)\n",
    "        # calculate size just in case\n",
    "        _anchors[\"size\"] = _anchors[_cluster_cols[2]] - _anchors[_cluster_cols[1]]\n",
    "        # return _anchors with coordinates rename as needed !\n",
    "        return _anchors.drop(columns=[\"start\",\"end\"]).rename(columns={\"cluster_start\":\"start\", \"cluster_end\":\"end\"})\n",
    "    elif mode == \"max_size\":\n",
    "        # return the largest anchor per cluster\n",
    "        # size of anchors, not clusters !!!\n",
    "        _anchors[\"size\"] = _anchors[_cols[2]] - _anchors[_cols[1]]\n",
    "        _largest_anchor_idx = _anchors.groupby(\"cluster\")[\"size\"].idxmax()\n",
    "        _anchors = _anchors.loc[_largest_anchor_idx]\n",
    "        # return _anchors - i.e. the largest anchor per cluster of overlaping anchors\n",
    "        return _anchors.drop(columns=[\"cluster_start\",\"cluster_end\"]).reset_index(drop=True)\n",
    "    elif mode == \"median\":\n",
    "        # return the of start and end coords per cluster of overlaping anchors ...\n",
    "        _anchors = _anchors.groupby(\"cluster\").agg({\"chrom\":\"first\", \"start\":\"median\", \"end\":\"median\"})\n",
    "        _anchors = _anchors.reset_index().astype({\"start\":int, \"end\":int})\n",
    "        _anchors[\"size\"] = _anchors[_cols[2]] - _anchors[_cols[1]]\n",
    "        # return _anchors - i.e. the largest anchor per cluster of overlaping anchors\n",
    "        return _anchors.reset_index(drop=True)\n",
    "    else:\n",
    "        raise ValueError()\n",
    "\n",
    "\n",
    "def calculate_valencies(\n",
    "    bed_df,  # must be output of bedpe_to_anchors, which in turn is a clustering inside\n",
    "    bedpe_df,\n",
    "    cluster_colname = \"cluster\",\n",
    "    valency_colname = \"valency\",\n",
    "    bed_cols = [\"chrom\", \"cluster_start\", \"cluster_end\"],\n",
    "    bedpe_cols1 = [\"chrom1\", \"start1\", \"end1\"],\n",
    "    bedpe_cols2 = [\"chrom2\", \"start2\", \"end2\"],\n",
    "):\n",
    "    \"\"\"\n",
    "    calculate valencies of a given anchors, given the bedpe ...\n",
    "    \"\"\"\n",
    "\n",
    "    if cluster_colname not in bed_df.columns:\n",
    "        raise ValueError(\"bed_df does not seem to be the result of bedpe_to_anchors/clustering ...\")\n",
    "\n",
    "    # overlap combined anchors with the left anchors to see how many \"dots\" we overlap ...\n",
    "    anchors_left = bioframe.overlap(\n",
    "        bed_df,\n",
    "        bedpe_df,\n",
    "        how='left',\n",
    "        cols1=bed_cols,\n",
    "        cols2=bedpe_cols1,\n",
    "    ).dropna( subset=[f\"{c}_\" for c in bedpe_cols1] )\n",
    "\n",
    "    # overlap combined anchors with the right anchors to see how many \"dots\" we overlap ...\n",
    "    anchors_right = bioframe.overlap(\n",
    "        bed_df,\n",
    "        bedpe_df,\n",
    "        how='left',\n",
    "        cols1=bed_cols,\n",
    "        cols2=bedpe_cols2,\n",
    "    ).dropna( subset=[f\"{c}_\" for c in bedpe_cols2] )\n",
    "\n",
    "    _num_clusters = len(bed_df)\n",
    "    # sanity check here ... - make sure we cover all of the cluster that are available ...\n",
    "    assert ( bed_df[cluster_colname].sort_values() == np.arange(_num_clusters) ).all()\n",
    "    # ...\n",
    "    _empty_clust_series = pd.Series(\n",
    "        data=np.zeros(_num_clusters),\n",
    "        index=pd.Index(data=np.arange(_num_clusters), name=cluster_colname),\n",
    "        name=\"count\"\n",
    "    )\n",
    "\n",
    "    # calculate valencies ...\n",
    "    _valencies = (_empty_clust_series + anchors_left[cluster_colname].value_counts()).fillna(0) \\\n",
    "                + (_empty_clust_series + anchors_right[cluster_colname].value_counts()).fillna(0)\n",
    "\n",
    "    # assign valencies back to anchors bed_df - carefully !\n",
    "    # _valencies are indexed using cluster_id - s\n",
    "    bed_df_clust_indexed = bed_df.set_index(cluster_colname)\n",
    "    bed_df_clust_indexed[valency_colname] = _valencies.astype(int)\n",
    "    #\n",
    "    return bed_df_clust_indexed.reset_index()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "790efa86-94e6-464a-817a-cd2449f1e490",
   "metadata": {},
   "source": [
    "# Let's draw an all-by-all grid of IDs ...\n",
    "... to demonstrate what's being called and what's not ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5eba5c9-7502-411b-bcf5-629ed68f5afd",
   "metadata": {},
   "outputs": [],
   "source": [
    "_anchors5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b9e9c5e-5381-4c3f-8eed-7fda6442d26b",
   "metadata": {},
   "source": [
    "### turn anchors into all-by-grid that fits in the region ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1332f365-40ce-4e59-b694-a07eb5f2aad2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "_region1 = ('chr6', 129_000_000, 129_000_000+1_500_000)\n",
    "_region2 = ('chr6', 129_000_000, 129_000_000+2_700_000)\n",
    "\n",
    "\n",
    "_anchors = _anchors5[[\"chrom\",\"peak_start\",\"peak_end\",\"cluster\"]].rename(columns={\"peak_start\":\"start\",\"peak_end\":\"end\"})\n",
    "_anchors_reg = bioframe.select(_anchors, _region2).reset_index(drop=True)\n",
    "_anchors_reg[\"cluster\"] = _anchors_reg[\"cluster\"] - _anchors_reg[\"cluster\"].min()\n",
    "_anchors_reg = calculate_valencies(\n",
    "    _anchors_reg,   # must be output of bedpe_to_anchors, which in turn is a clustering inside\n",
    "    clustered_pixels,\n",
    "    cluster_colname = \"cluster\",\n",
    "    valency_colname = \"valency\",\n",
    "    bed_cols = [\"chrom\", \"start\", \"end\"],\n",
    "    bedpe_cols1 = [\"chrom1\", \"start1\", \"end1\"],\n",
    "    bedpe_cols2 = [\"chrom2\", \"start2\", \"end2\"],\n",
    ")\n",
    "\n",
    "# # print(_anchors_reg)\n",
    "# _anchors_reg = _anchors_reg.query(\"valency > 1\").reset_index(drop=True)\n",
    "# # _anchors_reg = _anchors_reg.sample(n=10, weights=\"valency\").sort_index()\n",
    "\n",
    "_bedpe_region = bioframe.pair_by_distance(\n",
    "    _anchors_reg,\n",
    "    min_sep=0,\n",
    "    # max_sep=(_end - _start)+100_000,\n",
    "    max_sep=(_end - _start)+100_000_000,\n",
    "    suffixes=(\"1\",\"2\"),\n",
    "    keep_order=True,\n",
    ")\n",
    "\n",
    "# now find those all-by-allers that were actually called in the screening procedure ...\n",
    "\n",
    "_overlap1 = bioframe.overlap(\n",
    "    _bedpe_region,\n",
    "    clustered_pixels,\n",
    "    how='left',\n",
    "    cols1=[\"chrom1\",\"start1\",\"end1\"],\n",
    "    cols2=[\"chrom1\",\"start1\",\"end1\"],\n",
    "    suffixes=('', '_'),\n",
    ").dropna().set_index([\"cluster1\",\"cluster2\"])\n",
    "\n",
    "\n",
    "_overlap2 = bioframe.overlap(\n",
    "    _bedpe_region,\n",
    "    clustered_pixels,\n",
    "    how='left',\n",
    "    cols1=[\"chrom2\",\"start2\",\"end2\"],\n",
    "    cols2=[\"chrom2\",\"start2\",\"end2\"],\n",
    "    suffixes=('', '_'),\n",
    ").dropna().set_index([\"cluster1\",\"cluster2\"])\n",
    "\n",
    "\n",
    "_2d_mask = []\n",
    "for _c1, _c2 in _bedpe_region[[\"cluster1\",\"cluster2\"]].itertuples(index=False):\n",
    "    # print(_c1, _c2)\n",
    "    try:\n",
    "        _labs1 = _overlap1.loc[(_c1,_c2), \"labels_\"].to_list()\n",
    "        _labs2 = _overlap2.loc[(_c1,_c2), \"labels_\"].to_list()\n",
    "        _call_overlap = bool(set(_labs1) & set(_labs2))\n",
    "        _2d_mask.append(_call_overlap)\n",
    "    except KeyError:\n",
    "        _2d_mask.append(False)\n",
    "\n",
    "print(_2d_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f31c6f-ad80-48a6-88c7-405ec25ec1c4",
   "metadata": {},
   "source": [
    "# annotate bedpe region with bins and bin widths ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6db6e0a-67bf-46b0-9a57-cefe0d773694",
   "metadata": {},
   "outputs": [],
   "source": [
    "_bedpe_region[\"bin1_id\"] = _bedpe_region[[\"chrom1\",\"start1\",\"end1\"]].apply(clr.offset,axis=1,result_type=\"expand\")\n",
    "_bedpe_region[\"bin1_width\"] = _bedpe_region[[\"chrom1\",\"start1\",\"end1\"]].apply(clr.extent,axis=1,result_type=\"expand\").apply(np.diff,axis=1,result_type=\"expand\")[0]\n",
    "\n",
    "\n",
    "_bedpe_region[\"bin2_id\"] = _bedpe_region[[\"chrom2\",\"start2\",\"end2\"]].apply(clr.offset,axis=1,result_type=\"expand\")\n",
    "_bedpe_region[\"bin2_width\"] = _bedpe_region[[\"chrom2\",\"start2\",\"end2\"]].apply(clr.extent,axis=1,result_type=\"expand\").apply(np.diff,axis=1,result_type=\"expand\")[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a98f11e0-bde0-4061-a417-4c8f7cc8eda9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.colors import Normalize, TwoSlopeNorm\n",
    "\n",
    "# # https://stackoverflow.com/questions/48625475/python-shifted-logarithmic-colorbar-white-color-offset-to-center\n",
    "# class MidPointLogNorm(LogNorm):\n",
    "#     def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):\n",
    "#         LogNorm.__init__(self,vmin=vmin, vmax=vmax, clip=clip)\n",
    "#         self.midpoint=midpoint\n",
    "#     def __call__(self, value, clip=None):\n",
    "#         # I'm ignoring masked values and all kinds of edge cases to make a\n",
    "#         # simple example...\n",
    "#         x, y = [np.log(self.vmin), np.log(self.midpoint), np.log(self.vmax)], [0, 0.5, 1]\n",
    "#         return np.ma.masked_array(np.interp(np.log(value), x, y))\n",
    "\n",
    "\n",
    "\n",
    "# https://stackoverflow.com/questions/48625475/python-shifted-logarithmic-colorbar-white-color-offset-to-center\n",
    "class MidPointLogNorm(LogNorm):\n",
    "    def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):\n",
    "        LogNorm.__init__(self,vmin=vmin, vmax=vmax, clip=clip)\n",
    "        self.midpoint=midpoint\n",
    "\n",
    "    def __call__(self, value, clip=None):\n",
    "        # I'm ignoring masked values and all kinds of edge cases to make a\n",
    "        # simple example...\n",
    "        vmin, midpoint, vmax = self.vmin, self.midpoint, self.vmax\n",
    "        x, y = [np.log(vmin), np.log(midpoint), np.log(vmax)], [0, 0.5, 1]\n",
    "        return np.ma.masked_array(np.interp(np.log(value), x, y))\n",
    "\n",
    "    def inverse(self, value):\n",
    "        if not self.scaled():\n",
    "            raise ValueError(\"Not invertible until scaled\")\n",
    "        # t_vmin, t_midpoint, t_vmax = np.log(self.vmin), np.log(self.midpoint), np.log(self.vmax)\n",
    "        vmin, midpoint, vmax = self.vmin, self.midpoint, self.vmax\n",
    "\n",
    "        x, y = [0, 0.5, 1], [np.log(vmin), np.log(midpoint), np.log(vmax)]\n",
    "        # # return np.ma.masked_array(np.interp(np.log(value), x, y))\n",
    "        # if np.iterable(value):\n",
    "        #     val = np.ma.asarray(value)\n",
    "        #     return np.ma.power(val, 1. / gamma) * (vmax - vmin) + vmin\n",
    "        # else:\n",
    "        # return pow(value, 1. / gamma) * (vmax - vmin) + vmin\n",
    "        return np.exp(np.interp(value, x, y))\n",
    "\n",
    "\n",
    "\n",
    "imshow_kwargs = dict(\n",
    "        # norm=LogNorm(vmin=1/5, vmax=5),\n",
    "        norm=TwoSlopeNorm(vcenter=1, vmin=-0.1, vmax=3),\n",
    "        # norm=MidPointLogNorm(vmin=1/10, vmax=3, midpoint=1),\n",
    "        # cmap=\"RdBu_r\",\n",
    "        cmap=\"coolwarm\",\n",
    "        # interpolation=\"nearest\",\n",
    "        interpolation=\"none\",\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fab6985-c39e-4e91-b71f-9c879c24c220",
   "metadata": {},
   "outputs": [],
   "source": [
    "region1_name = bioframe.select(hg38_arms, _region1).iat[0,-1]\n",
    "region2_name = bioframe.select(hg38_arms, _region2).iat[0,-1]\n",
    "assert region1_name == region2_name\n",
    "region_name = region2_name\n",
    "#\n",
    "tile_span_i = clr.extent(_region1)\n",
    "tile_span_j = clr.extent(_region2)\n",
    "_the_tile = (region_name, tile_span_i, tile_span_j )\n",
    "\n",
    "tile_start_ij = (tile_span_i[0], tile_span_j[0])\n",
    "lazy_exp = LazyToeplitz(\n",
    "    exp_indexed.loc[region_name, region_name][\"balanced.avg\"].to_numpy()\n",
    ")\n",
    "# RAW observed matrix slice:\n",
    "observed = clr.matrix()[slice(*tile_span_i), slice(*tile_span_j)]\n",
    "expected = lazy_exp[slice(*tile_span_i), slice(*tile_span_j)]\n",
    "\n",
    "# let's figure out slices' coordinates ....\n",
    "_bins_i = clr.bins()[slice(*tile_span_i)]\n",
    "_bins_j = clr.bins()[slice(*tile_span_j)]\n",
    "_chrom_i, _start_i, _end_i = _bins_i.iloc[0][\"chrom\"], _bins_i.iloc[0][\"start\"], _bins_i.iloc[-1][\"end\"]\n",
    "_chrom_j, _start_j, _end_j = _bins_j.iloc[0][\"chrom\"], _bins_j.iloc[0][\"start\"], _bins_j.iloc[-1][\"end\"]\n",
    "\n",
    "\n",
    "# Start with a square Figure.\n",
    "fig = plt.figure(figsize=(18, 8))\n",
    "\n",
    "outer_grid = fig.add_gridspec(4, 1, wspace=0, hspace=0)\n",
    "\n",
    "ax = fig.add_subplot(outer_grid[0])\n",
    "axb = fig.add_subplot(outer_grid[1])\n",
    "axc = fig.add_subplot(outer_grid[2])\n",
    "axd = fig.add_subplot(outer_grid[3])\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "ax_x = divider.append_axes(\"bottom\", size=0.12, pad=0.00, sharex=ax)\n",
    "ax_y = divider.append_axes(\"left\", size=0.12, pad=0.00, sharey=ax)\n",
    "# create new axes on the right and on the top of the current axes.\n",
    "divider = make_axes_locatable(axb)\n",
    "axb_x = divider.append_axes(\"bottom\", size=0.12, pad=0.00, sharex=axb)\n",
    "axb_y = divider.append_axes(\"left\", size=0.12, pad=0.00, sharey=axb)\n",
    "# create new axes on the right and on the top of the current axes.\n",
    "divider = make_axes_locatable(axc)\n",
    "axc_x = divider.append_axes(\"bottom\", size=0.12, pad=0.00, sharex=axc)\n",
    "axc_y = divider.append_axes(\"left\", size=0.12, pad=0.00, sharey=axc)\n",
    "# create new axes on the right and on the top of the current axes.\n",
    "divider = make_axes_locatable(axd)\n",
    "axd_x = divider.append_axes(\"bottom\", size=0.12, pad=0.00, sharex=axd)\n",
    "axd_y = divider.append_axes(\"left\", size=0.12, pad=0.00, sharey=axd)\n",
    "\n",
    "\n",
    "gOE = scipy.ndimage.gaussian_filter(\n",
    "    (observed/expected),\n",
    "    sigma=0.4,\n",
    "    order=0,\n",
    "    mode='reflect',\n",
    "    cval=0.0,\n",
    "    # radius=3,\n",
    "    truncate=1.0,\n",
    ")\n",
    "\n",
    "\n",
    "hm_axes = [ax, axb, axc, axd]\n",
    "x_axes = [ax_x, axb_x, axc_x, axd_x]\n",
    "y_axes = [ax_y, axb_y, axc_y, axd_y]\n",
    "_boxes_to_draw = [clustered_pixels_old, enriched_pixels, clustered_pixels, _bedpe_region]\n",
    "_boxes_type = [\"box\", \"pixel\", \"box\", \"box\"]\n",
    "_alphas = [0.7, 0.2, 0.7, 0.7]\n",
    "\n",
    "for _ax, _axx, _axy, _bt, _boxes, _alpha in zip(\n",
    "    hm_axes,\n",
    "    x_axes,\n",
    "    y_axes,\n",
    "    _boxes_type,\n",
    "    _boxes_to_draw,\n",
    "    _alphas,\n",
    "):\n",
    "    #\n",
    "    _ax.imshow(\n",
    "        # (observed/expected),\n",
    "        gOE,\n",
    "        **imshow_kwargs,\n",
    "    )\n",
    "    # ...\n",
    "    if _bt == \"pixel\":\n",
    "        _draw_boxes = rectangles_around_dots\n",
    "        _boxes_kwargs = dict(loc=\"upper\", lw=1, ec=\"green\", fc=\"green\")\n",
    "    else:\n",
    "        _draw_boxes = rectangles_around_dots_ww\n",
    "        _boxes_kwargs = dict(loc=\"upper\", lw=1.5, ec=\"springgreen\", fc=\"none\", halo=0, ext_width=3)\n",
    "    # ...\n",
    "    for box in _draw_boxes(\n",
    "        _boxes,\n",
    "        _the_tile,\n",
    "        **_boxes_kwargs,\n",
    "    ):\n",
    "        _ax.add_patch(box)\n",
    "        x0, y0, dx, dy = box.get_bbox().bounds\n",
    "        _axx.axvspan(x0, x0+dx, facecolor='black', alpha=_alpha)\n",
    "        _axx.axvspan(y0, y0+dy, facecolor='black', alpha=_alpha)\n",
    "        _axy.axhspan(y0, y0+dy, facecolor='black', alpha=_alpha)\n",
    "        # _axy.axhspan(x0, x0+dx, facecolor='black', alpha=_alpha)\n",
    "        ########################################################\n",
    "    _ax.set_xticks([])\n",
    "    _ax.set_yticks([])\n",
    "    _axx.set_xticks([])\n",
    "    _axx.set_yticks([])\n",
    "    _axy.set_xticks([])\n",
    "    _axy.set_yticks([])\n",
    "\n",
    "\n",
    "plt.savefig(\"grid_csh_coolwarm_600.pdf\",dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eb965c6-f2a8-446e-b9fa-b92cb7dd7113",
   "metadata": {},
   "outputs": [],
   "source": [
    "_region = (_chrom , _start, _end)\n",
    "region_name = bioframe.select(hg38_arms, _region).iat[0,-1]\n",
    "tile_span_i = clr.extent(_region)\n",
    "tile_span_j = clr.extent(_region)\n",
    "\n",
    "_the_tile = (region_name,tile_span_i, tile_span_j )\n",
    "\n",
    "tile_start_ij = (tile_span_i[0], tile_span_j[0])\n",
    "lazy_exp = LazyToeplitz(\n",
    "    exp_indexed.loc[region_name, region_name][\"balanced.avg\"].to_numpy()\n",
    ")\n",
    "# RAW observed matrix slice:\n",
    "observed = clr.matrix()[slice(*tile_span_i), slice(*tile_span_j)]\n",
    "expected = lazy_exp[slice(*tile_span_i), slice(*tile_span_j)]\n",
    "\n",
    "# let's figure out slices' coordinates ....\n",
    "_bins_i = clr.bins()[slice(*tile_span_i)]\n",
    "_bins_j = clr.bins()[slice(*tile_span_j)]\n",
    "_chrom_i, _start_i, _end_i = _bins_i.iloc[0][\"chrom\"], _bins_i.iloc[0][\"start\"], _bins_i.iloc[-1][\"end\"]\n",
    "_chrom_j, _start_j, _end_j = _bins_j.iloc[0][\"chrom\"], _bins_j.iloc[0][\"start\"], _bins_j.iloc[-1][\"end\"]\n",
    "\n",
    "\n",
    "f, (axleft, axright) = plt.subplots(nrows=1,ncols=2,figsize=(24,11),sharex=True,sharey=True)\n",
    "# f.suptitle(f\"tile # {_tile_id} {(_chrom_i, _start_i, _end_i)} {(_chrom_j, _start_j, _end_j)}\",y=0.9)\n",
    "\n",
    "# print(f\"tile # {_tile_id} {(_chrom_i, _start_i, _end_i)} {(_chrom_j, _start_j, _end_j)}\")\n",
    "\n",
    "axleft.imshow(\n",
    "    observed,\n",
    "    cmap=\"YlOrBr\",\n",
    "    interpolation=\"none\",\n",
    "    norm=LogNorm(0.0001, 0.01)\n",
    ")\n",
    "axright.imshow(\n",
    "    (observed/expected),\n",
    "    cmap=\"coolwarm\",\n",
    "    interpolation=\"none\",\n",
    "    norm=LogNorm(0.25, 4)\n",
    ")\n",
    "\n",
    "\n",
    "for box in rectangles_around_dots_ww(\n",
    "    _bedpe_region,\n",
    "    _the_tile,\n",
    "    loc=\"upper\",\n",
    "    lw=2,\n",
    "    ec=\"darkgreen\",\n",
    "    fc=\"none\",\n",
    "    halo=100\n",
    "):\n",
    "    axleft.add_patch(box)\n",
    "for box in rectangles_around_dots_ww(\n",
    "    yyy,\n",
    "    _the_tile,\n",
    "    loc=\"upper\",\n",
    "    lw=2,\n",
    "    ec=\"crimson\",\n",
    "    fc=\"none\",\n",
    "    halo=100\n",
    "):\n",
    "    axleft.add_patch(box)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f053a07-c034-429f-9dd7-640311b7b05d",
   "metadata": {},
   "outputs": [],
   "source": [
    "_nnn = len(_bedpe_region[\"cluster1\"].unique())+1\n",
    "\n",
    "assert (_nnn*_nnn - _nnn)/2 == len(_bedpe_region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "614a15f5-ca6e-4ef8-9217-df7bf3668443",
   "metadata": {},
   "outputs": [],
   "source": [
    "_flank = 75_000 # Length of flank to one side from the boundary, in basepairs\n",
    "# create the stack of snips:\n",
    "_region_stack = cooltools.pileup(\n",
    "    clr,\n",
    "    _bedpe_region,\n",
    "    view_df=hg38_arms,\n",
    "    expected_df=exp,\n",
    "    flank=_flank,\n",
    "    nproc=12\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13f02613-1a49-40ee-bbff-5d0b2e952c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "_number = len(_bedpe_region[\"cluster1\"].unique())+1\n",
    "_width_unit = 2\n",
    "_height_unit = 2\n",
    "nrows, ncols = _number, _number\n",
    "f,axs = plt.subplots(\n",
    "    nrows=nrows,\n",
    "    ncols=ncols,\n",
    "    figsize = (nrows*_height_unit, ncols*_width_unit),\n",
    "    sharex=True,\n",
    "    sharey=True,\n",
    ")\n",
    "\n",
    "imshow_kwargs = dict(\n",
    "    norm=LogNorm(vmin=1/4,vmax=4),\n",
    "    # norm = mpl.colors.Normalize(vmin=1-0.9999,vmax=1+0.9999),\n",
    "    interpolation=\"none\",\n",
    "    # extent=[-_flank//1000, _flank//1000, -_flank//1000, _flank//1000],\n",
    "    cmap='RdBu_r',\n",
    ")\n",
    "\n",
    "_counter = 0\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        _ax = axs[i,j]\n",
    "        if j == i:\n",
    "            _ax.set_axis_off()\n",
    "            # _ax.set_xticks([])\n",
    "            # _ax.set_yticks([])\n",
    "        elif j > i:\n",
    "            _ax.set_xticks([])\n",
    "            _ax.set_yticks([])\n",
    "            _ax.imshow(_region_stack[_counter], **imshow_kwargs)\n",
    "            if _2d_mask[_counter]:\n",
    "                _ax.spines[\"left\"].set_color(\"limegreen\")\n",
    "                _ax.spines[\"right\"].set_color(\"limegreen\")\n",
    "                _ax.spines[\"top\"].set_color(\"limegreen\")\n",
    "                _ax.spines[\"bottom\"].set_color(\"limegreen\")\n",
    "                _ax.spines[\"left\"].set_linewidth(5)\n",
    "                _ax.spines[\"right\"].set_linewidth(5)\n",
    "                _ax.spines[\"top\"].set_linewidth(5)\n",
    "                _ax.spines[\"bottom\"].set_linewidth(5)\n",
    "            _counter += 1\n",
    "            #\n",
    "        else:\n",
    "            _ax.set_axis_off()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0865e4b0-d245-4df9-99b4-4843def77329",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12351d4a-19ff-4728-9c41-86211262dfc6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "328c5b7d-2e74-463a-a6f3-754598feabd8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}