Generate-Insulation-tracks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2f87134e-1293-4063-86f5-0e5e83d066f3",
   "metadata": {},
   "source": [
    "# Let's run insulation on all of the samples here - just in case ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f26519ff-e8c1-4307-866a-36c99cf5df76",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# 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",
    "\n",
    "# https://superfastpython.com/numpy-number-blas-threads/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f47b0626-f3d2-4d2e-b2c9-e6fcefe836b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import standard python libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import os, subprocess\n",
    "import seaborn as sns\n",
    "import multiprocess as mp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2be7179-fd58-44c3-8a2b-16e538203834",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import python package for working with cooler files and tools for analysis\n",
    "import cooler\n",
    "import cooltools.lib.plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86c1fe7f-22c4-4f54-9e8f-db43ca887064",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "# from saddle import saddleplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25c45aff-f437-4f6c-8609-18f97ea534b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# download test data\n",
    "# this file is 145 Mb, and may take a few seconds to download\n",
    "import cooltools\n",
    "from cooltools import insulation\n",
    "\n",
    "import bioframe\n",
    "from matplotlib.colors import LogNorm\n",
    "from helper_func import saddleplot\n",
    "from data_catalog import bws, bws_vlim, telo_dict, telo_reps_dict\n",
    "from helper_func import get_stack, show_stacks\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99d4774e-0092-4ff1-8deb-0cc4bfd13a14",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpire import WorkerPool\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4ad01fa-60b8-40f9-9252-f788a950546f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "from tqdm.notebook import trange, tqdm\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8566723-4cd5-4d8d-b66a-07d523fb72ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define genomic view that will be used to call dots and pre-compute expected\n",
    "\n",
    "# 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": "502a3ba5-c2fb-45ae-84f5-b4ed392ed8d9",
   "metadata": {},
   "source": [
    "# We'll just do 10kb binsize and 50/100 diamond-size and also 25 - diamond: 100; 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18b1e4e3-9788-4052-8723-5e77c32546f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# binsize = 10_000\n",
    "# windows = [5*binsize, 10*binsize]\n",
    "\n",
    "binsize10 = 10_000\n",
    "binsize25 = 25_000\n",
    "\n",
    "bins_diamonds_dict = {\n",
    "    binsize10: [5*binsize10, 10*binsize10],\n",
    "    binsize25: [4*binsize25, 6*binsize25],\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fc983d5-1ef2-43c9-9b5f-a92ddde54cec",
   "metadata": {},
   "source": [
    "## run the analysis on pooled replicates ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8be31637-d4b9-4913-adb7-44c760b17f36",
   "metadata": {},
   "outputs": [],
   "source": [
    "telo_ins = {}\n",
    "for binsize, windows in bins_diamonds_dict.items():\n",
    "    print(f\"working on {binsize} ...\")\n",
    "    # Pooled replicate coolers - the main working horse -ones ...\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",
    "\n",
    "    def _job(packed_data, sample):\n",
    "        # unpack data\n",
    "        clr_dict, wins, the_view = packed_data\n",
    "        ins_kwargs = dict(chunksize=20000000, nproc=12)\n",
    "        from cooltools import insulation\n",
    "        import warnings\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "            warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
    "            _clr = clr_dict[sample]\n",
    "            _ins = insulation(\n",
    "                _clr,\n",
    "                window_bp=wins,\n",
    "                view_df=the_view,\n",
    "                **ins_kwargs\n",
    "            )\n",
    "        return (sample, _ins)\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=(telo_clrs, windows, hg38_arms),\n",
    "        start_method=\"forkserver\",\n",
    "        use_dill=True,\n",
    "    ) as wpool:\n",
    "        results = wpool.map(_job, list(telo_clrs.keys()), progress_bar=True)\n",
    "\n",
    "    # sort out the results ...\n",
    "    telo_ins[binsize] = {sample: _ins for sample, _ins in results}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b11ca0a4-796c-46e4-a3c2-dccae49ce91d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# telo_ins[25_000][\"mMito\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf9c1d92-4a4b-4de9-b0cd-8b6709ea4c21",
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir ins_bedgraph"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f0c4b88-3a7f-45b7-ae94-c40d020230d6",
   "metadata": {},
   "source": [
    "### store bedGraphs ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdde662c-3c22-45f3-9e6e-020951b4e42a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for binsize in [10_000, 25_000,]:\n",
    "    # Write the insulation track as a bigwig:\n",
    "    for k, _fname in telo_dict.items():\n",
    "        # derive output name\n",
    "        out_fname = f\"{k}.{binsize//1_000}kb.bed\"\n",
    "        telo_ins[binsize][k].to_csv(\n",
    "            f\"ins_bedgraph/{out_fname}\",\n",
    "            sep=\"\\t\",\n",
    "            index=False,\n",
    "        )\n",
    "        print(f\"generated {out_fname} ...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67fa4535-6574-4b38-92d0-a95705913f03",
   "metadata": {},
   "source": [
    "## run analysis on replicates - separately ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d46c8e4-4132-4e5d-a43a-538b9f09e698",
   "metadata": {},
   "outputs": [],
   "source": [
    "telo_reps_ins = {}\n",
    "for binsize, windows in bins_diamonds_dict.items():\n",
    "    print(f\"working on {binsize} ...\")\n",
    "    # Pooled replicate coolers - the main working horse -ones ...\n",
    "    # cooler files that we'll work on :\n",
    "    telo_reps_clrs = { _k: cooler.Cooler(f\"{_path}::/resolutions/{binsize}\") for _k, _path in telo_reps_dict.items() }\n",
    "\n",
    "    def _job(packed_data, sample):\n",
    "        # unpack data\n",
    "        clr_dict, wins, the_view = packed_data\n",
    "        ins_kwargs = dict(chunksize=20000000, nproc=12)\n",
    "        from cooltools import insulation\n",
    "        import warnings\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "            warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
    "            _clr = clr_dict[sample]\n",
    "            _ins = insulation(\n",
    "                _clr,\n",
    "                window_bp=wins,\n",
    "                view_df=the_view,\n",
    "                **ins_kwargs\n",
    "            )\n",
    "        return (sample, _ins)\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=(telo_reps_clrs, windows, hg38_arms),\n",
    "        start_method=\"forkserver\",\n",
    "        use_dill=True,\n",
    "    ) as wpool:\n",
    "        results = wpool.map(_job, list(telo_reps_clrs.keys()), progress_bar=True)\n",
    "\n",
    "    # sort out the results ...\n",
    "    telo_reps_ins[binsize] = {sample: _ins for sample, _ins in results}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f6bcdf9-3fea-4af5-8795-cc53576ecfdf",
   "metadata": {},
   "source": [
    "### store R1/R2 bedGraphs ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8661b22c-d960-4db8-9a44-b132c2f8c48b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ! mkdir ins_r1r2_bedgraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97068a72-4854-4fce-99f0-8effb1e6d6d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "for binsize in [10_000, 25_000,]:\n",
    "    # Write the insulation track as a bigwig:\n",
    "    for k, _fname in telo_reps_dict.items():\n",
    "        # derive output name\n",
    "        out_fname = f\"{k}.{binsize//1_000}kb.bed\"\n",
    "        telo_reps_ins[binsize][k].to_csv(\n",
    "            f\"ins_r1r2_bedgraph/{out_fname}\",\n",
    "            sep=\"\\t\",\n",
    "            index=False,\n",
    "        )\n",
    "        print(f\"generated {out_fname} ...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db27eda4-033b-4e7b-b020-e5a946954834",
   "metadata": {},
   "source": [
    "# TODO - save insulation tracks in different ways and work on them a bit ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ff53075-7f3f-40e4-bb3d-44a3ad1394b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls *_bedgraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8411884-1466-4068-b69a-7da8c32a9515",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls *_bedraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7493aae-15ca-4ced-84f8-f8ff9cbe3286",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b78dc0b-3618-4540-b88a-1902d1d31fe6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faae1300-0bf0-4a78-a48b-f46751c10de4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_samples_m = [\n",
    "    \"mMito\",\n",
    "    \"mTelo\",\n",
    "    \"mCyto\",\n",
    "    \"m5hR1R2\",\n",
    "    \"m10hR1R2\",\n",
    "]\n",
    "sub_samples_p = [\n",
    "    \"pMito\",\n",
    "    \"pTelo\",\n",
    "    \"pCyto\",\n",
    "    \"p5hR1R2\",\n",
    "    \"p10hR1R2\",\n",
    "]\n",
    "\n",
    "fig, axs = plt.subplots(\n",
    "    nrows=len(sub_samples_m),\n",
    "    ncols=2,\n",
    "    figsize=(16,10),\n",
    "    sharex=True,\n",
    "    sharey=True,\n",
    ")\n",
    "\n",
    "# a contiguous region ...\n",
    "the_region = hg38_arms.set_index(\"name\").loc[\"chr6_q\"]\n",
    "the_region = \"chr6:80,000,000-130,000,000\"\n",
    "the_region = \"chr6:80,000,000-130,000,000\"\n",
    "_start = 120_000_000\n",
    "_width = 55_000_000\n",
    "the_region = (\"chr6\",_start,_start+_width)\n",
    "\n",
    "value_col = \"log2_insulation_score_100000\"\n",
    "\n",
    "for sample_m, sample_p, (i, axs) in zip(sub_samples_m, sub_samples_p, enumerate(axs)):\n",
    "    axm,axp = axs\n",
    "    mins = telo_ins[25_000][sample_m]\n",
    "    pins = telo_ins[25_000][sample_p]\n",
    "    # select for a region\n",
    "    mins = bioframe.select(mins, the_region)\n",
    "    pins = bioframe.select(pins, the_region)\n",
    "    axm.plot([0,len(mins)],[0,0],'k',lw=0.25)\n",
    "    axm.plot( np.arange(0, len(mins)), mins[value_col], label=value_col,linewidth=0.5, color=\"gray\")\n",
    "\n",
    "    axp.plot([0,len(pins)],[0,0],'k',lw=0.25)\n",
    "    axp.plot( np.arange(0, len(pins)), pins[value_col], label=value_col,linewidth=0.5, color=\"gray\")\n",
    "\n",
    "    axm.set_xlim(0, len(mins))\n",
    "    axp.set_xlim(0, len(pins))\n",
    "\n",
    "    axm.set_ylim(-1.1,1.1)\n",
    "    axp.set_ylim(-1.1,1.1)\n",
    "\n",
    "    axm.set_xticks([])\n",
    "    axm.set_yticks([])\n",
    "    axp.set_xticks([])\n",
    "    axp.set_yticks([])\n",
    "    if i == 0:\n",
    "        axm.set_title(\"m\")\n",
    "        axp.set_title(\"p\")\n",
    "    axm.set_ylabel(sample_m.lstrip(\"m\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c41a066a-96c6-463e-af43-d7d02043a2a3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0185d153-a819-4c46-b974-acc77b763b81",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a32a30f-3d0c-4db7-a911-6358c47b8584",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_samples_m = [\n",
    "    \"mMito_R1\",\n",
    "    \"mMito_R2\",\n",
    "    \"mTelo\",\n",
    "    \"mCyto\",\n",
    "    \"m5h_R1\",\n",
    "    \"m5h_R2\",\n",
    "    \"m10h_R1\",\n",
    "    \"m10h_R2\",\n",
    "]\n",
    "sub_samples_p = [\n",
    "    \"pMito_R1\",\n",
    "    \"pMito_R2\",\n",
    "    \"pTelo\",\n",
    "    \"pCyto\",\n",
    "    \"p5h_R1\",\n",
    "    \"p5h_R2\",\n",
    "    \"p10h_R1\",\n",
    "    \"p10h_R2\",\n",
    "]\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(\n",
    "    nrows=len(sub_samples_m),\n",
    "    ncols=2,\n",
    "    figsize=(16,10),\n",
    "    sharex=True,\n",
    "    sharey=True,\n",
    ")\n",
    "\n",
    "# a contiguous region ...\n",
    "the_region = hg38_arms.set_index(\"name\").loc[\"chr6_q\"]\n",
    "the_region = \"chr6:80,000,000-130,000,000\"\n",
    "the_region = \"chr6:80,000,000-130,000,000\"\n",
    "_start = 120_000_000\n",
    "_width = 55_000_000\n",
    "the_region = (\"chr6\",_start,_start+_width)\n",
    "\n",
    "value_col = \"log2_insulation_score_100000\"\n",
    "\n",
    "for sample_m, sample_p, (i, axs) in zip(sub_samples_m, sub_samples_p, enumerate(axs)):\n",
    "    axm,axp = axs\n",
    "    mins = telo_reps_ins[25_000][sample_m]\n",
    "    pins = telo_reps_ins[25_000][sample_p]\n",
    "    # select for a region\n",
    "    mins = bioframe.select(mins, the_region)\n",
    "    pins = bioframe.select(pins, the_region)\n",
    "    axm.plot([0,len(mins)],[0,0],'k',lw=0.25)\n",
    "    axm.plot( np.arange(0, len(mins)), mins[value_col], label=value_col,linewidth=0.5, color=\"gray\")\n",
    "\n",
    "    axp.plot([0,len(pins)],[0,0],'k',lw=0.25)\n",
    "    axp.plot( np.arange(0, len(pins)), pins[value_col], label=value_col,linewidth=0.5, color=\"gray\")\n",
    "\n",
    "    axm.set_xlim(0, len(mins))\n",
    "    axp.set_xlim(0, len(pins))\n",
    "\n",
    "    axm.set_ylim(-1.1,1.1)\n",
    "    axp.set_ylim(-1.1,1.1)\n",
    "\n",
    "    axm.set_xticks([])\n",
    "    axm.set_yticks([])\n",
    "    axp.set_xticks([])\n",
    "    axp.set_yticks([])\n",
    "    if i == 0:\n",
    "        axm.set_title(\"m\")\n",
    "        axp.set_title(\"p\")\n",
    "    axm.set_ylabel(sample_m.lstrip(\"m\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b82c455d-ef2b-4b92-943d-6a0b492f8034",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eac7e67b-471b-4121-85bd-2668bf3b3741",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec25b71f-7f2d-4d52-9656-4f37d6119a99",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cae927b3-8483-486b-b83c-de7f9b75753b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddd = {}\n",
    "value_col = \"log2_insulation_score_100000\"\n",
    "\n",
    "for sample in sub_samples_m+sub_samples_p:\n",
    "    ins = telo_ins[sample]\n",
    "\n",
    "    ddd[sample] = bioframe.select(ins, the_region)[value_col]\n",
    "    # ddd[sample] = ins[value_col]\n",
    "\n",
    "sns.heatmap(pd.DataFrame(ddd).corr(method='pearson'), annot=True, cmap=\"Reds\",vmin=0,vmax=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f4860a4-d7c5-4512-a697-9e80691ec00d",
   "metadata": {},
   "outputs": [],
   "source": [
    "telo_ins[\"pCyto\"][\"log2_insulation_score_100000\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55248356-b1af-48ff-8db5-d62b2397c4c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import decomposition\n",
    "# unused but required import for doing 3d projections with matplotlib < 3.2\n",
    "import mpl_toolkits.mplot3d  # noqa: F401\n",
    "\n",
    "_key = \"log2_insulation_score_100000\"\n",
    "\n",
    "_common_mask = pd.concat(\n",
    "    [_df[_key].rename(_sample) for _sample, _df in telo_ins.items()],\n",
    "    axis=1,\n",
    ").isna().any(axis=1).to_numpy()\n",
    "\n",
    "_pca_samples = [\n",
    "    # 'mTelo',\n",
    "    # 'mCyto',\n",
    "    'm5hR1R2',\n",
    "    'm10hR1R2',\n",
    "    # 'pTelo',\n",
    "    # 'pCyto',\n",
    "    'p5hR1R2',\n",
    "    'p10hR1R2',\n",
    "    # 'mp10hR1R2',\n",
    "    'N93m5',\n",
    "    'N93m10',\n",
    "    'N93p5',\n",
    "    'N93p10',\n",
    "    # 'N93mp10',\n",
    "]\n",
    "_rest_samples = [_sample for _sample in telo_ins if (_sample not in _pca_samples)]\n",
    "\n",
    "\n",
    "X = pd.concat(\n",
    "    [telo_ins[_sample][_key].rename(_sample) for _sample in _pca_samples],\n",
    "    axis=1,\n",
    ").to_numpy()[~_common_mask].T\n",
    "\n",
    "\n",
    "\n",
    "X_rest = pd.concat(\n",
    "    [telo_ins[_sample][_key].rename(_sample) for _sample in _rest_samples],\n",
    "    axis=1,\n",
    ").to_numpy()[~_common_mask].T\n",
    "\n",
    "\n",
    "\n",
    "pca = decomposition.PCA(n_components=5)\n",
    "print(\"running PCA ...\")\n",
    "pca.fit(X)\n",
    "X_trans = pca.transform(X)\n",
    "X_rest_trans = pca.transform(X_rest)\n",
    "\n",
    "print(pca.explained_variance_ratio_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68b92aa6-a346-40ee-9f4c-31ffa4483944",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "_x, _y, _z = X_trans[:,0], X_trans[:,1], X_trans[:,2]\n",
    "plt.scatter(_x, _y, s=50, color=\"red\")\n",
    "ax = plt.gca()\n",
    "for i, txt in enumerate(_pca_samples):\n",
    "    ax.annotate(txt, (_x[i], _y[i]))\n",
    "\n",
    "_x, _y, _z = X_rest_trans[:,0], X_rest_trans[:,1], X_rest_trans[:,2]\n",
    "plt.scatter(_x, _y, s=50, color=\"blue\")\n",
    "for i, txt in enumerate(_rest_samples):\n",
    "    ax.annotate(txt, (_x[i], _y[i]))\n",
    "\n",
    "\n",
    "\n",
    "_1, _2, _3 = pca.explained_variance_ratio_[:3]\n",
    "ax.set_xlabel(f\"pc1 {_1:.2f}\")\n",
    "ax.set_ylabel(f\"pc2 {_2:.2f}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8af21047-3216-4df2-84e7-bc7d81e7f69d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# # Write the insulation track as a bigwig:\n",
    "# for k, _fname in telo_dict.items():\n",
    "#     # # derive output name based on cooler's name\n",
    "#     # kkk = clr_fnames[k]\n",
    "#     _cname = _fname.split(\"/\")[-1]\n",
    "#     out_fname = \".\".join( _cname.split(\".\")[:-1] )\n",
    "#     # apparently insulation sometimes reports the same bin\n",
    "#     df = telo_ins[k].drop_duplicates(subset=[\"chrom\",\"start\",\"end\"])\n",
    "#     # let's save only 50_000 bp window derived insulation ...\n",
    "#     w = windows[0]\n",
    "#     # store in bigwig ...\n",
    "#     bioframe.to_bigwig(\n",
    "#         df,\n",
    "#         hg38_chromsizes,\n",
    "#         f\"{out_fname}.b{binsize}.insul.w{w}.bw\",\n",
    "#         value_field=f\"log2_insulation_score_{w}\",\n",
    "#         path_to_binary=\"./bedGraphToBigWig\"\n",
    "#     )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0c0e440-910f-40b4-988e-07089694196a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.filters import threshold_li, threshold_otsu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6332256e-5be8-4ac8-b91e-44712027135b",
   "metadata": {},
   "outputs": [],
   "source": [
    "histkwargs = dict(\n",
    "    bins=10**np.linspace(-4,1,200),\n",
    "    histtype='step',\n",
    "    lw=2,\n",
    ")\n",
    "\n",
    "f, axs = plt.subplots(len(windows), 1, sharex=True, figsize=(6,5))\n",
    "thresholds_li = {}\n",
    "thresholds_otsu = {}\n",
    "for i, (w, ax) in enumerate(zip(windows, axs)):\n",
    "    _sample = \"m5hR1R2\"\n",
    "    _ins = telo_ins[_sample]\n",
    "    ax.hist(\n",
    "        _ins[f'boundary_strength_{w}'],\n",
    "        label=_sample,\n",
    "        **histkwargs\n",
    "    )\n",
    "    _sample = \"p5hR1R2\"\n",
    "    _ins = telo_ins[_sample]\n",
    "    ax.hist(\n",
    "        _ins[f'boundary_strength_{w}'],\n",
    "        label=_sample,\n",
    "        **histkwargs\n",
    "    )\n",
    "    ax.text(0.01, 0.9,\n",
    "             f'Window {w//1000}kb',\n",
    "             ha='left',\n",
    "             va='top',\n",
    "             transform=ax.transAxes)\n",
    "    ax.set(\n",
    "        xscale='log',\n",
    "        ylabel='# boundaries'\n",
    "    )\n",
    "    ax.legend(frameon=False)\n",
    "\n",
    "axs[-1].set(xlabel='Boundary strength')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48f99412-bdcc-4820-868e-7bc39dcb4e8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "histkwargs = dict(\n",
    "    bins=10**np.linspace(-4,1,200),\n",
    "    histtype='step',\n",
    "    lw=2,\n",
    ")\n",
    "\n",
    "_sample = \"m5hR1R2\"\n",
    "_ins = telo_ins[_sample]\n",
    "f, axs = plt.subplots(len(windows), 1, sharex=True, figsize=(6,4))\n",
    "f.suptitle(_sample)\n",
    "for i, (w, ax) in enumerate(zip(windows, axs)):\n",
    "    _ins_values = _ins[f'boundary_strength_{w}']\n",
    "    ax.hist(\n",
    "        _ins_values,\n",
    "        label=_sample,\n",
    "        **histkwargs\n",
    "    )\n",
    "    _li = threshold_li(_ins_values.dropna().values)\n",
    "    _otsu = threshold_otsu(_ins_values.dropna().values)\n",
    "    n_boundaries_li = (_ins_values.dropna()>=_li).sum()\n",
    "    n_boundaries_otsu = (_ins_values.dropna()>=_otsu).sum()\n",
    "    ax.axvline(_li, c='green')\n",
    "    ax.axvline(_otsu, c='magenta')\n",
    "    ax.text(0.01, 0.9,\n",
    "             f'Window {w//1000}kb',\n",
    "             ha='left',\n",
    "             va='top',\n",
    "             transform=ax.transAxes)\n",
    "    ax.text(0.01, 0.7,\n",
    "            f'{n_boundaries_otsu} boundaries (Otsu)',\n",
    "            c='magenta',\n",
    "            ha='left',\n",
    "            va='top',\n",
    "            transform=ax.transAxes)\n",
    "    ax.text(0.01, 0.5,\n",
    "            f'{n_boundaries_li} boundaries (Li)',\n",
    "            c='green',\n",
    "            ha='left',\n",
    "            va='top',\n",
    "            transform=ax.transAxes)\n",
    "    ax.set(\n",
    "        xscale='log',\n",
    "        ylabel='# boundaries'\n",
    "    )\n",
    "axs[-1].set(xlabel='Boundary strength')\n",
    "\n",
    "\n",
    "_sample = \"p5hR1R2\"\n",
    "_ins = telo_ins[_sample]\n",
    "f, axs = plt.subplots(len(windows), 1, sharex=True, figsize=(6,4))\n",
    "f.suptitle(_sample)\n",
    "for i, (w, ax) in enumerate(zip(windows, axs)):\n",
    "    _ins_values = _ins[f'boundary_strength_{w}']\n",
    "    ax.hist(\n",
    "        _ins_values,\n",
    "        label=_sample,\n",
    "        **histkwargs\n",
    "    )\n",
    "    _li = threshold_li(_ins_values.dropna().values)\n",
    "    _otsu = threshold_otsu(_ins_values.dropna().values)\n",
    "    n_boundaries_li = (_ins_values.dropna()>=_li).sum()\n",
    "    n_boundaries_otsu = (_ins_values.dropna()>=_otsu).sum()\n",
    "    ax.axvline(_li, c='green')\n",
    "    ax.axvline(_otsu, c='magenta')\n",
    "    ax.text(0.01, 0.9,\n",
    "             f'Window {w//1000}kb',\n",
    "             ha='left',\n",
    "             va='top',\n",
    "             transform=ax.transAxes)\n",
    "    ax.text(0.01, 0.7,\n",
    "            f'{n_boundaries_otsu} boundaries (Otsu)',\n",
    "            c='magenta',\n",
    "            ha='left',\n",
    "            va='top',\n",
    "            transform=ax.transAxes)\n",
    "    ax.text(0.01, 0.5,\n",
    "            f'{n_boundaries_li} boundaries (Li)',\n",
    "            c='green',\n",
    "            ha='left',\n",
    "            va='top',\n",
    "            transform=ax.transAxes)\n",
    "    ax.set(\n",
    "        xscale='log',\n",
    "        ylabel='# boundaries'\n",
    "    )\n",
    "axs[-1].set(xlabel='Boundary strength')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da00334f-e85c-4dc0-b823-65c540ca2023",
   "metadata": {},
   "source": [
    "# Let's save some of those boundaries ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f66a259e-d66d-4ede-951d-ab836d44da36",
   "metadata": {},
   "outputs": [],
   "source": [
    "_sample = \"m5hR1R2\"\n",
    "_ins = telo_ins[_sample]\n",
    "_ins_values = _ins[f'boundary_strength_50000']\n",
    "_li = threshold_li(_ins_values.dropna().values)\n",
    "_otsu = threshold_otsu(_ins_values.dropna().values)\n",
    "# pick one threshold and extract the boundaries ...\n",
    "_df = (_ins[_ins_values >= _otsu][[\"chrom\",\"start\",\"end\"]]).reset_index(drop=True)\n",
    "_df.to_csv(f\"{_sample}_otsu_ins_boundaries_w50kb.bed\",index=False,sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caaa07c3-bb28-4540-8295-7af090125380",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls *.bed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a22235e7-6476-43b6-bd92-c1e5f2cb5e54",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4339aa2e-9b45-4940-b6ad-aef76107b996",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8679d300-4359-440e-ad5b-e5662b4fafe7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "016e8256-89dc-4d3c-a7b6-91257702154d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84f3b8f6-d0fc-4815-90e1-96f5c0a8f863",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b533522-9ef8-4a76-b7a2-65c933f282ad",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d27a0f0-5774-406d-a1e2-6e44bb930988",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f78386d9-5944-4f33-a0b7-4bed4203f97c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cbd454c-12fa-4694-bbee-43f99be9f1ee",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c683a35a-8f9c-4ec8-b578-705c6dc64cf2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f04fef4-e43d-49b2-956b-97c6cd3de861",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c3190c1-4dfe-4b4c-844e-3e200d61af02",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f597c11-15c6-48bb-b12e-f2fd9a10c17b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4608c471-853d-4923-8412-a2bb81beda1e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57794820-d0c9-4fbd-8ee2-b602ee5a049e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd29d81a-60a9-412d-8192-861b2cfeb45c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd257ac4-60f0-47c8-8180-5801fb7623dd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7014cf61-a147-4eb1-87ef-aed49677c152",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13481040-5945-473c-95ce-30a9a02a9aeb",
   "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
}