What_we_see_in_Hi-C_8Jul2024.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "97dced82",
   "metadata": {},
   "source": [
    "# What we see in Hi-C: Introduction to chromatin biology and sequencing assays, Hi-C data browsing and feature extraction\n",
    "\n",
    "Welcome to the workshop \"What can we see in Hi-C?\" by Aleksandra Galitsyna.\n",
    "Your assistants for this workshop are Sofya Gaydukova and Joseph Finkelberg.\n",
    "\n",
    "This is a two-part intro to Hi-C data analysis:\n",
    "\n",
    "- Part I covers qualitative Hi-C data analysis in interactive browsers (higlass, resgen)\n",
    "- Part II shows how to use flexibility of Python programming for whole-genome visualization, compartments, domains and scaling analysis of Hi-C with [Open2C tools](https://github.com/open2c).\n",
    "\n",
    "Disclaimer: although we formally refer to whole-genome chromosome capture data as Hi-C, we will mostly work with Micro-C here.\n",
    "\n",
    "Author: Aleksandra Galitsyna, based on the materials from [Open2C](https://github.com/open2c) and [4DN](https://www.4dnucleome.org/software/) communities. \n",
    "\n",
    "For any questions, feel free to reach out to galitsyn at mit dot edu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4521d59-0d8d-45ad-8a3b-22301b74c52a",
   "metadata": {},
   "source": [
    "### 0. Open this notebook on a machine and install the requirements\n",
    "\n",
    "You have two options to run this notebook: \n",
    "\n",
    "#### Plan A. Open in Google Colab and run the following lines:\n",
    "\n",
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/mirnylab/GenomeArchitectureAndFunction2024/blob/main/What_we_see_in_Hi-C_8Jul2024.ipynb\">\n",
    "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
    "</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca854ba8-63c2-4eb8-8318-fea6077d93ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install ipykernel cooler==0.10.2 bioframe==0.7.2 higlass-python==1.1.0 pysam\n",
    "# !pip install git+https://github.com/open2c/cooltools.git # Note that this will ensure you can use automatic data download API\n",
    "# !pip install matplotlib==3.8.* seaborn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6559f3ef-17f9-47bd-99ca-1aefab6c30d2",
   "metadata": {},
   "source": [
    "#### Plan B. Open some Linux-based server, launch terminal, install dependencies, git-clone and open the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29ad5a23-d441-4156-ae64-a3c35ee97294",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%bash \n",
    "# conda create -n hic_open2c_workshop python=3.11\n",
    "# conda activate hic_open2c_workshop\n",
    "# pip install ipykernel cooler==0.10.2 bioframe==0.7.2 higlass-python==1.1.0 pysam\n",
    "# pip install git+https://github.com/open2c/cooltools.git\n",
    "# pip install matplotlib==3.8.* seaborn\n",
    "# python -m ipykernel install --user --name hic_open2c_workshop --display-name \"hic_open2c_workshop\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6b699a-4994-4608-9942-642914b5c6ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !git clone https://github.com/mirnylab/GenomeArchitectureAndFunction2024.git"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3ed7344-97ea-4c66-b66e-70620e72e3da",
   "metadata": {},
   "source": [
    "Select the 'hic_open2c_workshop' enviroment in your notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32772246-7ed0-4780-bcad-9d694d7f8949",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "Please, note that the workshop is designed to run in the Linux-based system with conda and pip access. \n",
    "If you are using MacOS or Windows, you may install the dependencies in the way preferable for these systems. \n",
    "However, we cannot guarantee that all elements of this notebook will work in these systems. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e39d180e-8452-40e7-83ea-2b9ed2308c00",
   "metadata": {},
   "source": [
    "#### Check that imports run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bb157c3-2a33-423c-96a6-d4d12c709fe0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import higlass as hg\n",
    "import cooltools\n",
    "import cooler\n",
    "import bioframe\n",
    "import matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8190c25-fb14-420c-be86-5a0473552fc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Matplotlib: \", matplotlib.__version__)\n",
    "print(\"HiGlass Python: \", hg.__version__)\n",
    "print(\"Cooler: \", cooler.__version__)\n",
    "print(\"Bioframe: \", bioframe.__version__)\n",
    "print(\"Cooltools: \", cooltools.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9d69743-f806-4d33-aaba-0e66b0959f49",
   "metadata": {},
   "source": [
    "No error during running these cells means that the installation and preparation of the environment was successful.\n",
    "Warnings are also fine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a4b6b4e-a93b-4d40-a1bf-8a63318dc62f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's ignore warnings for now:\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34a11a42",
   "metadata": {},
   "source": [
    "## Part I. Qualitative analysis of Hi-C maps\n",
    "\n",
    "Qualitative Hi-C data anlaysis is the first step for generation of hypothesis, quality checks, etc. It involves eye-balling, loading and comparing various datasets from genomics and epigenetics. Here, we will learn how to do that with interactive browsers. \n",
    "\n",
    "\n",
    "### 1. Databases for Hi-C-type data\n",
    "\n",
    "#### 1.1. Find public Hi-C dataset.\n",
    "\n",
    "Hi-C data can be found at various public databases, including broad-range sequencing storages (GEO, SRA, ENA, ArrayExpress, DDBJ), [ENCODE](https://www.encodeproject.org/matrix/?type=Experiment&control_type!=*&status=released&perturbed=false&assay_title=intact+Hi-C&assay_title=in+situ+Hi-C&assay_title=dilution+Hi-C) and [4DN Data Portal](https://data.4dnucleome.org/). \n",
    "\n",
    "We will practice with 4DN. \n",
    "[4DN](https://data.4dnucleome.org/) deposits raw sequencing data, processed data in widely accepted formats (.hic and .mcool) and annotation tracks (such as insulation and compartments).\n",
    "\n",
    "\n",
    "▸ **Practice task. Get familiar with 4DN database.**\n",
    "\n",
    "- Go to the 4DN data portal website: https://data.4dnucleome.org/\n",
    "- Select Micro-C for h1-ESC cells (\"Micro-C\" -> \"Browse\" -> Line \"4DNES21D8SP8\")\n",
    "- Explore the data (press \"Explore data button\"m you will be redirected to HiGlass view with Micro-C for h1-ESC cells).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21312b2a-842e-40ab-b1c2-aa3acb7d29ca",
   "metadata": {},
   "source": [
    "### 2. Interactive browsers for Hi-C\n",
    "\n",
    "Let's now explore Hi-C-like data and learn how to create your own views with genomics and epigenetics datasets. \n",
    "In this part, you will explore and compare chromatin organization of **mESC (mouse embryonic stem cells) cells probed by Micro-C** (high-resolution version of Hi-C). \n",
    "\n",
    "We will explore Micro-C data from untreated cells and from Rad21, Wapl adn CTCF rapid degratation in mESC. \n",
    "You can look up the details on the expertiment in the original publication: [Hsieh et al. 2022](https://www.nature.com/articles/s41588-022-01223-8).\n",
    "\n",
    "#### 2.1. Explore Hi-C and epigenetics data with higlass and resgen (follow the lecture).\n",
    "\n",
    "For your reference, we will cover navigating the datasets in:\n",
    "\n",
    "- [higlass.io](https://higlass.io/), interactive browser for 3D and other genomics data\n",
    "- [resgen.io](https://resgen.io/), server version of higlass, with option to deposit your own data, organize it into projects and store the views\n",
    "\n",
    "** Disclaimer: it is highly recommended not to access resgen.io at the time of the workshop, as it might have an effect of \"spam attack\" on the server. **\n",
    "\n",
    "#### 2.2. Explore Hi-C data in your notebook with higlass-python.\n",
    "\n",
    "Next, let's explore Hi-C data of mESC in our notebooks with [higlass-python](https://github.com/higlass/higlass-python), a Python wrapper for higlass. \n",
    "\n",
    "You may notice how similar are the generated python widgets to the higlass.io interface.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1732705d-4a16-46e3-99ab-1623ac6b15c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import higlass as hg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61708a91-b86f-40d4-a164-c74a00dcd639",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define the data tilesets:\n",
    "mm10_genes = hg.remote(\n",
    "    uid=\"QDutvmyiSrec5nX4pA5WGQ\",\n",
    "    server=\"//higlass.io/api/v1\",\n",
    ")\n",
    "\n",
    "tileset = hg.remote(\n",
    "    uid=\"WmNzNu-XTieALMScNwrCZA\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"Hsieh et al. (2021) mESC Micro-C\",\n",
    ")\n",
    "\n",
    "### Define the view:\n",
    "view = hg.view(\n",
    "    mm10_genes.track(\"gene-annotations\", height=100).opts(\n",
    "        minHeight = 24,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    tileset.track(\"heatmap\", height=500).opts(\n",
    "        minHeight = 500,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    genomePositionSearchBox = hg.GenomePositionSearchBox(\n",
    "        autocompleteServer=\"//higlass.io/api/v1\",\n",
    "        autocompleteId=\"QDutvmyiSrec5nX4pA5WGQ\",\n",
    "        chromInfoId=\"mm10\",\n",
    "        chromInfoServer=\"//higlass.io/api/v1\",\n",
    "        visible=True)\n",
    ")\n",
    "\n",
    "# Help to get access to track sources from higlass.io data sources:\n",
    "list_of_track_source_servers = [\n",
    "    \"//higlass.io/api/v1\",\n",
    "    \"https://resgen.io/api/v1\",\n",
    "    \"https://resgen.io/api/v1/gt/paper-data\"\n",
    "  ]\n",
    "\n",
    "# Show the view:\n",
    "view.viewconf(trackSourceServers = list_of_track_source_servers, exportViewUrl = \"/api/v1/viewconfs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7554bc81-349e-4c50-940f-483c4c1c5a52",
   "metadata": {},
   "source": [
    "▸ **Practice task.** Zoom in and out Micro-C map. Change the color bar. \n",
    "\n",
    "Wait for several seconds if the view is still updating.\n",
    "\n",
    "What chromatin structures can you find? Find inter-chromosomal compartments, intra-chromosomal compartments, then domains, dots and stripes. \n",
    "\n",
    "⏸ *5-10 mintes, feel free to ask questions.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14bd13ee-de47-4047-b7ec-da3da18bf124",
   "metadata": {},
   "source": [
    "#### 2.3. Compare Micro-C data between mESC and rapid Rad21 degradation.\n",
    "\n",
    "Let's compare two views, for untreated mESC (left) and IAA-treated Rad21-degron mESCs:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab4a1d86-b901-499e-bd14-10126c8366ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "mm10_genes = hg.remote(\n",
    "    uid=\"QDutvmyiSrec5nX4pA5WGQ\",\n",
    "    server=\"//higlass.io/api/v1\",\n",
    ")\n",
    "\n",
    "tileset_heatmap_1 = hg.remote(\n",
    "    uid=\"WmNzNu-XTieALMScNwrCZA\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"Hsieh et al. (2021) mESC Micro-C\",\n",
    ")\n",
    "\n",
    "tileset_heatmap_2 = hg.remote(\n",
    "    uid=\"HJBk7MibTpaAo2ZAP7lG_A\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC Micro-C with Rad21 IAA\",\n",
    ")\n",
    "\n",
    "view1 = hg.view(\n",
    "    mm10_genes.track(\"gene-annotations\", height=100).opts(\n",
    "        minHeight = 24,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    tileset_heatmap_1.track(\"heatmap\", height=500).opts(\n",
    "        minHeight = 500,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    genomePositionSearchBox = hg.GenomePositionSearchBox(\n",
    "        autocompleteServer=\"//higlass.io/api/v1\",\n",
    "        autocompleteId=\"QDutvmyiSrec5nX4pA5WGQ\",\n",
    "        chromInfoId=\"mm10\",\n",
    "        chromInfoServer=\"//higlass.io/api/v1\",\n",
    "        visible=True), \n",
    "    width=6\n",
    ")\n",
    "\n",
    "view2 = hg.view(\n",
    "    mm10_genes.track(\"gene-annotations\", height=100).opts(\n",
    "        minHeight = 24,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    tileset_heatmap_2.track(\"heatmap\", height=500).opts(\n",
    "        minHeight = 500,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    width=6\n",
    ")\n",
    "\n",
    "#In order to get access to track sources from higlass.io data sources\n",
    "list_of_track_source_servers = [\n",
    "    \"//higlass.io/api/v1\",\n",
    "    \"https://resgen.io/api/v1\",\n",
    "    \"https://resgen.io/api/v1/gt/paper-data\"\n",
    "  ]\n",
    "\n",
    "view_lock = hg.lock(view1, view2)\n",
    "viewconf = (view1 | view2).locks(view_lock, inplace=True)\n",
    "viewconf.trackSourceServers = list_of_track_source_servers\n",
    "viewconf.exportViewUrl = \"/api/v1/viewconfs\"\n",
    "\n",
    "viewconf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9d05665-8ab0-46e5-b426-e2985ffef4c2",
   "metadata": {},
   "source": [
    "▸ **Practice task.** What changes between Micro-C maps do you see? \n",
    "\n",
    "Do domains, dots and stripes change? And compartments? \n",
    "\n",
    "⏸ *5-10 mintes, feel free to ask questions.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "664b1de4-2dac-4f35-8d61-c52151d2ba44",
   "metadata": {},
   "source": [
    "#### 2.4. Compare Micro-C datasets and associate the changes with epigenetics\n",
    "\n",
    "Let's add more information to the view: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffce0f59-b81a-45a4-9434-3bffceff44f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "mm10_genes = hg.remote(\n",
    "    uid=\"QDutvmyiSrec5nX4pA5WGQ\",\n",
    "    server=\"//higlass.io/api/v1\",\n",
    ")\n",
    "\n",
    "tileset_heatmap_1 = hg.remote(\n",
    "    uid=\"WmNzNu-XTieALMScNwrCZA\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"Hsieh et al. (2021) mESC Micro-C\",\n",
    ")\n",
    "\n",
    "tileset_heatmap_2 = hg.remote(\n",
    "    uid=\"HJBk7MibTpaAo2ZAP7lG_A\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC Micro-C with Rad21 IAA\",\n",
    ")\n",
    "\n",
    "tileset_ctcf = hg.remote(\n",
    "    uid=\"H65ZUeoUSXudbspjx4LJhA\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC CTCF ChIP-Seq coverage\",\n",
    ")\n",
    "\n",
    "tileset_H3K4me1 = hg.remote(\n",
    "    uid=\"R-8frGVgSSejDb3HTdHcvw\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC H3K4me1 ChIP-Seq coverage\",\n",
    ")\n",
    "\n",
    "tileset_H3K9me3 = hg.remote(\n",
    "    uid=\"G4twI42bS-GIWIymaLreBA\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC H3K9me3 ChIP-Seq coverage\",\n",
    ")\n",
    "\n",
    "tileset_H3K27ac = hg.remote(\n",
    "    uid=\"WXUOiZkYRTiUtBIAP9Cy5Q\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC H3K27ac ChIP-Seq coverage\",\n",
    ")\n",
    "\n",
    "tileset_H3K4me3 = hg.remote(\n",
    "    uid=\"CmtC5DBbTEKhicuuL1nNOw\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC H3K4me3 ChIP-Seq coverage\",\n",
    ")\n",
    "\n",
    "tileset_H3K27me3 = hg.remote(\n",
    "    uid=\"PMar87awRW2MaLLsTvimcw\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"mESC H3K27me3 ChIP-Seq coverage\",\n",
    ")\n",
    "\n",
    "# Epilogos, or chromatin states are more complicated to include: \n",
    "from typing import ClassVar\n",
    "from typing_extensions import Literal\n",
    "\n",
    "class MultivecTrack(hg.PluginTrack):\n",
    "    type: Literal[\n",
    "        \"horizontal-stacked-bar\",\n",
    "    ]\n",
    "    plugin_url: ClassVar[str] = (\n",
    "        \"https://unpkg.com/higlass-multivec/dist/higlass-multivec.min.js\"\n",
    "    )\n",
    "\n",
    "track_epilogos = MultivecTrack(\n",
    "    **{\n",
    "        \"type\": \"horizontal-stacked-bar\",\n",
    "        \"tilesetUid\": \"WmHBFc78TmOgTCIePEwawQ\",\n",
    "        \"server\": \"https://resgen.io/api/v1\",\n",
    "        \"height\": 100,\n",
    "        \"minHeight\": 100,\n",
    "    }\n",
    ")\n",
    "\n",
    "view1 = hg.view(\n",
    "    mm10_genes.track(\"gene-annotations\", height=100).opts(\n",
    "        minHeight = 24,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    (track_epilogos, \"top\"),\n",
    "    tileset_heatmap_1.track(\"heatmap\", height=500).opts(\n",
    "        minHeight = 500,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    (tileset_ctcf.track(\"horizontal-bar\", width=30).opts(\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "        barFillColor=\"black\",\n",
    "    ), \"left\"),\n",
    "    (tileset_H3K4me1.track(\"horizontal-bar\", width=30).opts(\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "        barFillColor=\"yellow\",\n",
    "    ), \"left\"),\n",
    "    (tileset_H3K4me3.track(\"horizontal-bar\", width=30).opts(\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "        barFillColor=\"red\",\n",
    "    ), \"left\"),\n",
    "    (tileset_H3K27ac.track(\"horizontal-bar\", width=30).opts(\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "        barFillColor=\"orange\",\n",
    "    ), \"left\"),\n",
    "    (tileset_H3K9me3.track(\"horizontal-bar\", width=30).opts(\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "        barFillColor=\"gray\",\n",
    "    ), \"left\"),\n",
    "    (tileset_H3K27me3.track(\"horizontal-bar\", width=30).opts(\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "        barFillColor=\"blue\",\n",
    "    ), \"left\"),\n",
    "    genomePositionSearchBox = hg.GenomePositionSearchBox(\n",
    "        autocompleteServer=\"//higlass.io/api/v1\",\n",
    "        autocompleteId=\"QDutvmyiSrec5nX4pA5WGQ\",\n",
    "        chromInfoId=\"mm10\",\n",
    "        chromInfoServer=\"//higlass.io/api/v1\",\n",
    "        visible=True), \n",
    "    width=6\n",
    ")\n",
    "\n",
    "view2 = hg.view(\n",
    "    mm10_genes.track(\"gene-annotations\", height=100).opts(\n",
    "        minHeight = 24,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    (track_epilogos, \"top\"),\n",
    "    tileset_heatmap_2.track(\"heatmap\", height=500).opts(\n",
    "        minHeight = 500,\n",
    "        showMousePosition=True,\n",
    "        showTooltip=True,\n",
    "    ),\n",
    "    width=6\n",
    ")\n",
    "\n",
    "#In order to get access to track sources from higlass.io data sources\n",
    "list_of_track_source_servers = [\n",
    "    \"//higlass.io/api/v1\",\n",
    "    \"https://resgen.io/api/v1\",\n",
    "    \"https://resgen.io/api/v1/gt/paper-data\"\n",
    "  ]\n",
    "\n",
    "view_lock = hg.lock(view1, view2)\n",
    "viewconf = (view1 | view2).locks(view_lock, inplace=True)\n",
    "viewconf.trackSourceServers = list_of_track_source_servers\n",
    "viewconf.exportViewUrl = \"/api/v1/viewconfs\"\n",
    "\n",
    "viewconf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b8d2a2b-293e-4a11-a361-ec288f82815d",
   "metadata": {},
   "source": [
    "! Warning: \n",
    "Run the cell again if it produces an error! \n",
    "\n",
    "! Warning: Navigation might be slower due to larger amounts of data you view at the same time! "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7964335c-a08f-4f5d-a01c-6972cb8c513d",
   "metadata": {},
   "source": [
    "▸ **Practice task.** Navigate to the gene-scale view. For example, search for one of the genes CSNK1A1L, HMCN1, GPATCH2, or LPIN1, and zom out to +- 1-5 Mb around that gene. \n",
    "\n",
    "What histone modifications are enriched around these genes? \n",
    "\n",
    "What structures are associated with these modifications? \n",
    "\n",
    "Are these structures affected by Rad21 degradation? \n",
    "\n",
    "**Hints and further directions:**\n",
    "\n",
    "Find H3K27ac-enriched regions. What happens to them? \n",
    "\n",
    "Find H3K27me3-enriched regions. What happens to them? \n",
    "\n",
    "⏸ *10-15 mintes break to rest, ask questions and complete this open-end task.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b6f752e",
   "metadata": {},
   "source": [
    "## Part II. Quantitative analysis of Hi-C maps\n",
    "\n",
    "Computationally, we think of Hi-C map as a 2-dimensional matrix with numbers. Each row and column corresponds to the bin in the genome, and the number represents the probability of interactions between corresponding two bins. Let's learn how to work with it in Python. \n",
    "\n",
    "This part of workshop is based on the efforts of [Open Chromosome Community (Open2C)](https://open2c.github.io/) of developers who create the tools for simpler and better analysis of Hi-C. Check out our software: https://github.com/open2c and Twitter where we frequently post tutorials, preprints and papers: https://twitter.com/Open2C_team \n",
    "\n",
    "\n",
    "\n",
    "For more detailed examples of Hi-C data analysis, see Open2C examples: https://github.com/open2c/open2c_examples\n",
    "\n",
    "\n",
    "**Important notes:**\n",
    "- In this practical task, you will learn how to work with the same mESC Micro-C datasets.\n",
    "- We selected only two chromosomes (chr1 and chr12) from this rich data to make the computations faster. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35f72da",
   "metadata": {},
   "source": [
    "### 1. Load and navigate Hi-C in Python\n",
    "\n",
    "Import python libraries that we will need in almost all insulation/compartments/scalings analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3de4400",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Standard python libraries:\n",
    "\n",
    "# Data analysis/manipulation:\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# visualization:\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['font.size'] = 8\n",
    "from matplotlib import colors\n",
    "import mpl_toolkits\n",
    "import seaborn as sns\n",
    "\n",
    "# Operation system and subprocesses control:\n",
    "import os\n",
    "import subprocess\n",
    "\n",
    "###  Libraries for Hi-C and epigenetics data analysis\n",
    "import cooler\n",
    "import cooltools\n",
    "import bioframe"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a45f9b21",
   "metadata": {},
   "source": [
    "Now, let's list the datasets avaliable to you. We provide them as part of cooltools standard data. You need mESC Micro-C data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d96fdb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cooltools.print_available_datasets()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31d4cbb7",
   "metadata": {},
   "source": [
    "Download test data for mESC cells, or set up the oaths manually to the downloaded folder from Google Drive: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc5abee9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# download test data\n",
    "# this file is ~200 Mb, and may take a few seconds to download\n",
    "\n",
    "data_dir = './data/'\n",
    "\n",
    "# Online version: \n",
    "cool_file_ut = cooltools.download_data(\"mESC_dRAD21_UT\", cache=True, data_dir=data_dir)    \n",
    "\n",
    "# Offline version:\n",
    "# Or: cool_file_ut = \"./data/mESC_dRAD21_UT.chr1-chr12-only.mm10.mapq_30.1000.mcool\"\n",
    "print(cool_file_ut)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba170ab",
   "metadata": {},
   "source": [
    "Repeat the same for Rad21 and other degradations: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67d490e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "cool_file_rad21 = cooltools.download_data(\"mESC_dRAD21_IAA\", cache=True, data_dir=data_dir)    \n",
    "cool_file_ctcf = cooltools.download_data(\"mESC_dCTCF_IAA\", cache=True, data_dir=data_dir)    \n",
    "cool_file_wapl = cooltools.download_data(\"mESC_dWAPL_IAA\", cache=True, data_dir=data_dir)    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "583176e2-2291-4aa4-b515-2e9831b2eca6",
   "metadata": {},
   "source": [
    "Let's learn more about `cool` data and its content:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a310890c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to print which resolutions are stored in the mcool, use list_coolers\n",
    "cooler.fileops.list_coolers(cool_file_ut)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa06f506",
   "metadata": {},
   "source": [
    "Load cooler file to Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3835eaf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "### to load a cooler with a specific resolution use the following syntax:\n",
    "clr = cooler.Cooler(f'{cool_file_ut}::resolutions/1000000')\n",
    "\n",
    "### to print chromosomes and binsize for this cooler\n",
    "print(f'chromosomes: {clr.chromnames}, binsize: {clr.binsize}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b7dbaf9-d826-4808-a02f-24817bc0f3cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "clr.bins()[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7278b350-fcb8-406f-82e4-b11d3755c8f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "clr.pixels()[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ab9ed2",
   "metadata": {},
   "outputs": [],
   "source": [
    "### load Rad21 IAA cooler: \n",
    "clr_rad21 = cooler.Cooler(f'{cool_file_rad21}::resolutions/1000000')\n",
    "\n",
    "### to print chromosomes and binsize for this cooler\n",
    "print(f'chromosomes: {clr.chromnames}, binsize: {clr.binsize}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f843df0d",
   "metadata": {},
   "source": [
    "#### 1.1. View multiple chromosomes\n",
    "\n",
    "Create the list of available chromosomes for visualization:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f3b6fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "### to make a list of chromosome start/ends in bins:\n",
    "chromstarts = []\n",
    "for i in clr.chromnames:\n",
    "    print(f'{i} : {clr.extent(i)}')\n",
    "    chromstarts.append(clr.extent(i)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78805d07",
   "metadata": {},
   "source": [
    "Load the whole matrix into the memory: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb5e4391",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = clr.matrix(balance=True)[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae472935",
   "metadata": {},
   "source": [
    "Visualize both available chromosomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f089d467",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an image:\n",
    "f, ax = plt.subplots(figsize=(7,6))\n",
    "\n",
    "# Display cooler matrix: \n",
    "im = ax.matshow(data, vmax=1e-2); \n",
    "\n",
    "plt.colorbar(im, fraction=0.046, pad=0.04, label='normalized counts');\n",
    "\n",
    "ax.set(xticks=chromstarts, \n",
    "       xticklabels=clr.chromnames,\n",
    "       xlabel='position, chrom#', \n",
    "       ylabel='position, bin#')\n",
    "ax.xaxis.set_label_position('top')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "757354f0-5b3b-453c-b85d-7890992a489d",
   "metadata": {},
   "source": [
    "#### 1.2. View regions of chromosomes\n",
    "\n",
    "Below, we fetch and plot an individual chromosome (left) and a region of a chromosome (right) using `clr.fetch()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ead5be8-21d3-4987-ad08-bbb6c606cf87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to plot ticks in terms of megabases we use the EngFormatter \n",
    "# https://matplotlib.org/gallery/api/engineering_formatter.html\n",
    "from matplotlib.ticker import EngFormatter\n",
    "bp_formatter = EngFormatter('b')\n",
    "\n",
    "def format_ticks(ax, x=True, y=True, rotate=True):\n",
    "    if y:\n",
    "        ax.yaxis.set_major_formatter(bp_formatter)\n",
    "    if x:\n",
    "        ax.xaxis.set_major_formatter(bp_formatter)\n",
    "        ax.xaxis.tick_bottom()\n",
    "    if rotate:\n",
    "        ax.tick_params(axis='x',rotation=45)\n",
    "\n",
    "# Create the figure: \n",
    "f, axs = plt.subplots(\n",
    "    figsize=(14,4),\n",
    "    ncols=3)\n",
    "\n",
    "# Left figure: all data\n",
    "data = clr.matrix(balance=True)[:]\n",
    "\n",
    "ax = axs[0]\n",
    "im = ax.matshow(data, vmax=1e-2); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_xticks(chromstarts)\n",
    "ax.set_xticklabels(clr.chromnames)\n",
    "ax.set_yticks(chromstarts)\n",
    "ax.set_yticklabels(clr.chromnames)\n",
    "ax.xaxis.tick_bottom()\n",
    "ax.set_title('All data')\n",
    "\n",
    "# Center figure: chr17\n",
    "data_chr17 = clr.matrix(balance=True).fetch('chr1')\n",
    "\n",
    "ax = axs[1]\n",
    "im = ax.matshow(\n",
    "    data_chr17,\n",
    "    vmax=5e-2, \n",
    "    extent=(0, clr.chromsizes['chr1'], clr.chromsizes['chr1'], 0)\n",
    ");\n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title('chr1', y=1.08)\n",
    "ax.set_ylabel('position, Mb')\n",
    "format_ticks(ax)\n",
    "\n",
    "# Right figure: region of chr17\n",
    "start, end = 30_000_000, 60_000_000\n",
    "region = ('chr1', start, end)\n",
    "data_region = clr.matrix(balance=True).fetch(region)\n",
    "\n",
    "ax = axs[2]\n",
    "im = ax.matshow(\n",
    "    data_region,\n",
    "    vmax=5e-2, \n",
    "    extent=(start, end, end, start)\n",
    "); \n",
    "ax.set_title(f'chr1:{start:,}-{end:,}', y=1.08)\n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "format_ticks(ax)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d8407ef-8e29-4c30-b371-f9ca2dc541b6",
   "metadata": {},
   "source": [
    "#### 1.3. Logarithmic color scale\n",
    "\n",
    "Hi-C/Micro-C data has a high dynamic range, we often plot the data in log-scale. This enables simultaneous visualization of features near and far from the diagonal in a consistent colorscale. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3e8e993-3d27-4761-a787-490d97a92f15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot heatmaps at megabase resolution with 3 levels of zoom in log-scale with a consistent colormap#\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "# Create the figure: \n",
    "f, axs = plt.subplots(\n",
    "    figsize=(14,4),\n",
    "    ncols=3)\n",
    "bp_formatter = EngFormatter('b')\n",
    "\n",
    "# Normalized scale:\n",
    "norm = LogNorm(vmax=1e-1)\n",
    "\n",
    "# Left figure: all data\n",
    "data = clr.matrix(balance=True)[:]\n",
    "\n",
    "ax = axs[0]\n",
    "im = ax.matshow(\n",
    "    data, \n",
    "    norm=norm,\n",
    ") \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_xticks(chromstarts)\n",
    "ax.set_xticklabels(clr.chromnames)\n",
    "ax.set_yticks(chromstarts)\n",
    "ax.set_yticklabels(clr.chromnames)\n",
    "ax.xaxis.tick_bottom()\n",
    "ax.set_title('All data')\n",
    "\n",
    "# Center figure: chr17\n",
    "data_chr17 = clr.matrix(balance=True).fetch('chr1')\n",
    "\n",
    "ax = axs[1]\n",
    "im = ax.matshow(\n",
    "    data_chr17,\n",
    "    norm=norm,\n",
    "    extent=(0, clr.chromsizes['chr1'], clr.chromsizes['chr1'], 0)\n",
    ");\n",
    "plt.colorbar(im, ax=ax ,fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title('chr1', y=1.08)\n",
    "ax.set(ylabel='position, Mb', xlabel='position, Mb')\n",
    "format_ticks(ax)\n",
    "\n",
    "\n",
    "# Right figure: region of chr17\n",
    "start, end = 30_000_000, 60_000_000\n",
    "region = ('chr1', start, end)\n",
    "data_region = clr.matrix(balance=True).fetch(region)\n",
    "\n",
    "ax = axs[2]\n",
    "im = ax.matshow(\n",
    "    data_region,\n",
    "    norm=norm,\n",
    "    extent=(start, end, end, start)\n",
    "); \n",
    "ax.set_title(f'chr17:{start:,}-{end:,}', y=1.08)\n",
    "plt.colorbar(im, ax=ax ,fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set(xlabel='position, Mb')\n",
    "format_ticks(ax)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b88af6b-dddf-4150-9318-98a69eff4bdd",
   "metadata": {},
   "source": [
    "#### 1.4. Colormaps\n",
    "\n",
    "In different papers people prefer to use various color schemes that represent the data better. \n",
    "\n",
    "`cooltools.lib.plotting` registers a set of colormaps that are useful for visualizing C data.\n",
    "In particular, the `fall` colormap (inspired by [colorbrewer](https://colorbrewer2.org/#type=sequential&scheme=YlOrRd&n=9)) offers a high dynamic range, linear, option for visualizing Hi-C matrices. This often displays features more clearly than red colormaps.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0474792-442d-4c4f-9f29-fc031384fa5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot the corrected data in fall heatmap and compare to the white-red colormap ###\n",
    "### thanks for the alternative colormap naming to https://twitter.com/HiC_memes/status/1286326919122825221/photo/1###\n",
    "import cooltools.lib.plotting\n",
    "\n",
    "vmax = 5e-2\n",
    "norm = LogNorm(vmax=vmax)\n",
    "fruitpunch = sns.blend_palette(['white', 'red'], as_cmap=True)\n",
    "\n",
    "f, axs = plt.subplots(\n",
    "    figsize=(13, 10),\n",
    "    nrows=2, \n",
    "    ncols=2,\n",
    "    sharex=True, sharey=True)\n",
    "\n",
    "ax = axs[0, 0]\n",
    "ax.set_title('Pumpkin Spice')\n",
    "im = ax.matshow(clr.matrix(balance=True)[:], vmax=vmax, cmap='fall'); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts (linear)');\n",
    "plt.xticks(chromstarts,clr.chromnames);\n",
    "\n",
    "ax = axs[0, 1]\n",
    "ax.set_title('Fruit Punch')\n",
    "im3 = ax.matshow(clr.matrix(balance=True)[:], vmax=vmax, cmap=fruitpunch); \n",
    "plt.colorbar(im3, ax=ax, fraction=0.046, pad=0.04, label='norm counts (linear)');\n",
    "plt.xticks(chromstarts,clr.chromnames);\n",
    "\n",
    "ax = axs[1, 0]\n",
    "im = ax.matshow(clr.matrix(balance=True)[:], norm=norm, cmap='fall'); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts (log)');\n",
    "plt.xticks(chromstarts,clr.chromnames);\n",
    "\n",
    "ax = axs[1, 1]\n",
    "im3 = ax.matshow(clr.matrix(balance=True)[:], norm=norm, cmap=fruitpunch); \n",
    "plt.colorbar(im3, ax=ax, fraction=0.046, pad=0.04, label='norm counts (log)');\n",
    "plt.xticks(chromstarts,clr.chromnames);\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08a68106-df9d-4767-ba3d-71ae0bb2dd1f",
   "metadata": {},
   "source": [
    "The utility of fall colormaps becomes more noticeable at higher resolutions and higher degrees of zoom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be3bc963-1128-41bf-90c7-883c175b978e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cooltools.lib.plotting\n",
    "clr_10kb = cooler.Cooler(f'{cool_file_ut}::resolutions/10000')\n",
    "\n",
    "region = 'chr12:30,000,000-35,000,000'\n",
    "chrom, start, end = bioframe.parse_region_string(region)\n",
    "extents = (start, end, end, start)\n",
    "\n",
    "vmax=5e-2\n",
    "norm = LogNorm(vmax=vmax)\n",
    "\n",
    "f, axs = plt.subplots(\n",
    "    figsize=(13, 10),\n",
    "    nrows=2, \n",
    "    ncols=2,\n",
    "    sharex=True, \n",
    "    sharey=True\n",
    ")\n",
    "\n",
    "ax = axs[0, 0]\n",
    "im = ax.matshow(\n",
    "    clr_10kb.matrix(balance=True).fetch(region),\n",
    "    cmap='fall',\n",
    "    vmax=vmax,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "\n",
    "ax = axs[0, 1]\n",
    "im2 = ax.matshow(\n",
    "    clr_10kb.matrix(balance=True).fetch(region),\n",
    "    cmap=fruitpunch, \n",
    "    vmax=vmax,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im2, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "\n",
    "ax = axs[1, 0]\n",
    "im = ax.matshow(\n",
    "    clr_10kb.matrix(balance=True).fetch(region),\n",
    "    cmap='fall',\n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "\n",
    "ax = axs[1, 1]\n",
    "im2 = ax.matshow(\n",
    "    clr_10kb.matrix(balance=True).fetch(region),\n",
    "    cmap=fruitpunch, \n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im2, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "\n",
    "for ax in axs.ravel():\n",
    "    format_ticks(ax, rotate=False)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "811dbe5a",
   "metadata": {},
   "source": [
    "#### 1.5. Compare Micro-C for untreated and IAA-treated Rad21 degron cells\n",
    "\n",
    "▸ **Practice task.** \n",
    "\n",
    "Select the region you like (might be the region you explored in resgen). Make sure you select ~ 1-5 Mb size region from chr1 or chr12.\n",
    "\n",
    "Visualize that region with the color scheme and scale that you find instructuve and interesting. \n",
    "\n",
    "Repeat the same for both mESC and IAA-treated Rad21 degron cells. What is the major difference between them? \n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db120749",
   "metadata": {},
   "outputs": [],
   "source": [
    "clr_10kb = cooler.Cooler(f'{cool_file_ut}::resolutions/10000')\n",
    "clr_10kb_rad21 = cooler.Cooler(f'{cool_file_rad21}::resolutions/10000')\n",
    "\n",
    "region = ... #### you may select any other region from chr1 or chr 12\n",
    "\n",
    "chrom, start, end = bioframe.parse_region_string(region)\n",
    "extents = (start, end, end, start)\n",
    "\n",
    "vmax=5e-2\n",
    "norm = LogNorm(vmax=vmax)\n",
    "\n",
    "f, axs = plt.subplots(\n",
    "    figsize=(13, 10),\n",
    "    nrows=1, \n",
    "    ncols=2,\n",
    "    sharex=True, \n",
    "    sharey=True\n",
    ")\n",
    "\n",
    "ax = axs[0]\n",
    "im = ax.matshow(\n",
    "    clr_10kb.matrix(balance=True).fetch(region),\n",
    "    cmap='fall',\n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title(\"Untreated\")\n",
    "\n",
    "ax = axs[1]\n",
    "im2 = ax.matshow(\n",
    "    clr_10kb_rad21.matrix(balance=True).fetch(region),\n",
    "    cmap='fall', \n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im2, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title(\"dRad21 IAA\")\n",
    "\n",
    "for ax in axs.ravel():\n",
    "    format_ticks(ax, rotate=False)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51e03e6f-4c63-4b30-884b-e4f3a080bbc9",
   "metadata": {},
   "source": [
    "▸ **Practice task.**\n",
    "\n",
    "Load Wapl and CTCF degraded datasets. Show the same region for them.\n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "704d59a1-1948-44dd-839c-3bde6c7250fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "clr_10kb_wapl = cooler.Cooler(f'{cool_file_wapl}::resolutions/10000')\n",
    "clr_10kb_ctcf = cooler.Cooler(f'{cool_file_ctcf}::resolutions/10000')\n",
    "\n",
    "region = ... #### you may select any other region\n",
    "\n",
    "chrom, start, end = bioframe.parse_region_string(region)\n",
    "extents = (start, end, end, start)\n",
    "\n",
    "vmax=5e-2\n",
    "norm = LogNorm(vmax=vmax)\n",
    "\n",
    "f, axs = plt.subplots(\n",
    "    figsize=(17, 10),\n",
    "    nrows=1, \n",
    "    ncols=4,\n",
    "    sharex=True, \n",
    "    sharey=True\n",
    ")\n",
    "\n",
    "ax = axs[0]\n",
    "im = ax.matshow(\n",
    "    clr_10kb.matrix(balance=True).fetch(region),\n",
    "    cmap='fall',\n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title(\"Untreated\")\n",
    "\n",
    "ax = axs[1]\n",
    "im2 = ax.matshow(\n",
    "    clr_10kb_rad21.matrix(balance=True).fetch(region),\n",
    "    cmap='fall', \n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im2, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title(\"dRad21 IAA\")\n",
    "\n",
    "ax = axs[2]\n",
    "im3 = ax.matshow(\n",
    "    clr_10kb_wapl.matrix(balance=True).fetch(region),\n",
    "    cmap='fall',\n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im3, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title(\"dWAPL IAA\")\n",
    "\n",
    "ax = axs[3]\n",
    "im4 = ax.matshow(\n",
    "    clr_10kb_ctcf.matrix(balance=True).fetch(region),\n",
    "    cmap='fall', \n",
    "    norm=norm,\n",
    "    extent=extents\n",
    "); \n",
    "plt.colorbar(im4, ax=ax, fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title(\"dCTCF IAA\")\n",
    "\n",
    "for ax in axs.ravel():\n",
    "    format_ticks(ax, rotate=False)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1222348b",
   "metadata": {},
   "source": [
    "### 2. Insulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a797f980-55fd-477a-830f-8197e74edea3",
   "metadata": {},
   "source": [
    "Insulation is a simple concept, yet a powerful way to look at Hi-C/Micro-C data. Insulation is one aspect of locus-specific contact frequency at small genomic distances, and reflects the segmentation of the genome into domains.\n",
    "\n",
    "Insulation can be computed with multiple methods. One of the most common methods involves using a diamond-window score to generate an ***insulation profile***. To compute this profile, slide a diamond-shaped window along the genome, with one of the corners on the main diagonal of the matrix, and sum up the contacts within the window for each position.\n",
    "\n",
    "Insulation profiles reveal that certain locations have lower scores, reflecting lowered contact frequencies between upstream and downstream loci. These positions are often referred to as ***boundaries***, and are also obtained with multiple methods. Here we illustrate one thresholding method for determining boundaries from an insulation profile.\n",
    "\n",
    "<img src=\"https://www.biorxiv.org/content/biorxiv/early/2022/11/01/2022.10.31.514564/F4.large.jpg\"\n",
    "     alt=\"Insulation and boundaries\"\n",
    "     style=\"float: left; margin-right: 10px;\" />\n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85a4ceaa-84f9-4dd8-8e41-9b017cf9df90",
   "metadata": {},
   "source": [
    "#### 2.1. Calculate genome-wide contact insulation\n",
    "Here we load the Micro-C data at 10 kbp resolution and calculate insulation score with 4 different window sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07f27cd1-c1d9-44b3-83eb-4e3c12bc5012",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cooltools.lib.plotting\n",
    "from cooltools import insulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd21ade2-1b7a-49cf-bd7a-663c1c7fc8c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 10_000  \n",
    "clr_10k = cooler.Cooler(f'{cool_file_ut}::resolutions/{resolution}')\n",
    "windows = [3*resolution, 5*resolution, 10*resolution, 25*resolution]\n",
    "insulation_table = insulation(clr_10k, windows, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae34475-e9be-484f-8f80-da79490a94cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "windows # These are insulation windows that we will run analysis with"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed44a3e6-08f5-448e-b425-727eda0805f1",
   "metadata": {},
   "source": [
    "`insulation_table` is a dataframe where rows correspond to genomic bins of the cooler.\n",
    "\n",
    "The columns of this insulation dataframe report the insulation score, the number of valid (non-nan) pixels, whether the given bin is valid, the boundary prominence (strength) and whether locus is called as a boundary after thresholding, for each of the window sizes provided to the function.\n",
    "\n",
    "Below we print the information returned for any window size, as well as the specific information for the largest window used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e6bdb41-6233-43ee-81a1-db3c620ea0f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "first_window_summary = insulation_table.columns[[ str(windows[-1]) in i for i in insulation_table.columns]]\n",
    "\n",
    "insulation_table[['chrom','start','end','region','is_bad_bin']+list(first_window_summary)].iloc[1000:1007]\n",
    "# Note how insulation boundary looks like in this dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c6c2d67",
   "metadata": {},
   "source": [
    "Let's visualize insulation. For that, we will rotate the contact matrix by 45 degrees, and we need some custom (yet helpful) code for that: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "997cb37a-d5ed-42a3-a2ef-f21c40e3a218",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions to help with plotting\n",
    "def pcolormesh_45deg(ax, matrix_c, start=0, resolution=1, *args, **kwargs):\n",
    "    start_pos_vector = [start+resolution*i for i in range(len(matrix_c)+1)]\n",
    "    import itertools\n",
    "    n = matrix_c.shape[0]\n",
    "    t = np.array([[1, 0.5], [-1, 0.5]])\n",
    "    matrix_a = np.dot(np.array([(i[1], i[0])\n",
    "                                for i in itertools.product(start_pos_vector[::-1],\n",
    "                                                           start_pos_vector)]), t)\n",
    "    x = matrix_a[:, 1].reshape(n + 1, n + 1)\n",
    "    y = matrix_a[:, 0].reshape(n + 1, n + 1)\n",
    "    im = ax.pcolormesh(x, y, np.flipud(matrix_c), *args, **kwargs)\n",
    "    im.set_rasterized(True)\n",
    "    return im\n",
    "\n",
    "from matplotlib.ticker import EngFormatter\n",
    "bp_formatter = EngFormatter('b')\n",
    "def format_ticks(ax, x=True, y=True, rotate=True):\n",
    "    if y:\n",
    "        ax.yaxis.set_major_formatter(bp_formatter)\n",
    "    if x:\n",
    "        ax.xaxis.set_major_formatter(bp_formatter)\n",
    "        ax.xaxis.tick_bottom()\n",
    "    if rotate:\n",
    "        ax.tick_params(axis='x',rotation=45)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c52ce345-6eed-4191-b291-34ba544ae10f",
   "metadata": {},
   "source": [
    "Let's see what the insulation track at the highest resolution looks like, next to a rotated Hi-C matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8d944a3-ff0f-465e-8dad-5b606187d53a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import additional plotting utilities and set up some parameters: \n",
    "from matplotlib.colors import LogNorm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "plt.rcParams['font.size'] = 8\n",
    "\n",
    "# Load the data:\n",
    "start = 10_500_000\n",
    "end = start+ 90*windows[0]\n",
    "region = ('chr1', start, end)\n",
    "norm = LogNorm(vmax=0.1, vmin=0.001)\n",
    "data = clr_10k.matrix(balance=True).fetch(region)\n",
    "\n",
    "# Create the figure: \n",
    "f, ax = plt.subplots(figsize=(18, 6))\n",
    "im = pcolormesh_45deg(ax, data, start=region[1], resolution=resolution, norm=norm, cmap='fall')\n",
    "ax.set_aspect(0.5)\n",
    "ax.set_ylim(0, 10*windows[0])\n",
    "format_ticks(ax, rotate=False)\n",
    "ax.xaxis.set_visible(False)\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"1%\", pad=0.1, aspect=6)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# Select insulation for a region:\n",
    "insul_region = bioframe.select(insulation_table, region)\n",
    "\n",
    "# Plot insulation: \n",
    "ins_ax = divider.append_axes(\"bottom\", size=\"50%\", pad=0., sharex=ax)\n",
    "ins_ax.set_prop_cycle(plt.cycler(\"color\", plt.cm.plasma(np.linspace(0,1,5))))\n",
    "ins_ax.plot(insul_region.loc[:, ['start', 'end']].mean(axis=1).values, \n",
    "            insul_region.loc[:, 'log2_insulation_score_'+str(windows[0])].values,\n",
    "            label=f'Window {windows[0]} bp')\n",
    "\n",
    "ins_ax.legend(bbox_to_anchor=(0., -1), loc='lower left', ncol=4);\n",
    "\n",
    "format_ticks(ins_ax, y=False, rotate=False)\n",
    "ax.set_xlim(region[1], region[2])\n",
    "\n",
    "# Add multiple resolutions:\n",
    "for res in windows[1:]:\n",
    "    ins_ax.plot(insul_region.loc[:, ['start', 'end']].mean(axis=1).values, \n",
    "                insul_region.loc[:, f'log2_insulation_score_{res}'].values, \n",
    "                label=f'Window {res} bp')\n",
    "ins_ax.legend(bbox_to_anchor=(0., -1), loc='lower left', ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b38afe61-3301-40e7-be99-47d2cc76141c",
   "metadata": {},
   "source": [
    "This result highlights how much the result is dependent on window size: smaller windows are sensitive to local structure, whereas large windows capture regions that insulate at larger scales."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42bb265e-4182-4cfa-a695-6ce57586ce8d",
   "metadata": {},
   "source": [
    "#### 2.2. Boundary calling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d87d7bb-eb3b-454b-bf59-3b45cd498428",
   "metadata": {},
   "source": [
    "The insulation table also has annotations for valleys of the insulation score, which correspond to highly insulating regions, such as TAD boundaries. All potential boundaries have an assigned `boundary_strength_` column. Additionally, this strength is thresholded to find regions that insulate particularly strongly, and this is recorded in the `is_boundary_` columns."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd75c4b9-df1d-45f9-9fd2-77a4a799723f",
   "metadata": {},
   "source": [
    "Let's repeat the previous plot and show where we found the boundaries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b0c77cc-98d3-4fce-8534-ee8d97be8e6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_window_size = 100_000 # More smooth and representative window size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61c2f470-9bba-44af-82c1-fbce3e5fb882",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure and plot Micro-C on it:\n",
    "f, ax = plt.subplots(figsize=(20, 10))\n",
    "im = pcolormesh_45deg(ax, data, start=region[1], resolution=resolution, norm=norm, cmap='fall')\n",
    "ax.set_aspect(0.5)\n",
    "ax.set_ylim(0, 10*windows[0])\n",
    "format_ticks(ax, rotate=False)\n",
    "ax.xaxis.set_visible(False)\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"1%\", pad=0.1, aspect=6)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# Select insulation table for a region:\n",
    "insul_region = bioframe.select(insulation_table, region)\n",
    "\n",
    "# Plot insulation: \n",
    "ins_ax = divider.append_axes(\"bottom\", size=\"50%\", pad=0., sharex=ax)\n",
    "\n",
    "ins_ax.plot(insul_region[['start', 'end']].mean(axis=1),\n",
    "            insul_region[f'log2_insulation_score_{selected_window_size}'], label=f'Window {selected_window_size} bp')\n",
    "\n",
    "# Define weak and strong boundaries:\n",
    "boundaries = insul_region[~np.isnan(insul_region[f'boundary_strength_{selected_window_size}'])]\n",
    "weak_boundaries = boundaries[~boundaries[f'is_boundary_{selected_window_size}']]\n",
    "strong_boundaries = boundaries[boundaries[f'is_boundary_{selected_window_size}']]\n",
    "\n",
    "# Plot the boundaries:\n",
    "ins_ax.scatter(weak_boundaries[['start', 'end']].mean(axis=1), \n",
    "            weak_boundaries[f'log2_insulation_score_{selected_window_size}'], label='Weak boundaries')\n",
    "ins_ax.scatter(strong_boundaries[['start', 'end']].mean(axis=1), \n",
    "            strong_boundaries[f'log2_insulation_score_{selected_window_size}'], label='Strong boundaries')\n",
    "\n",
    "ins_ax.legend(bbox_to_anchor=(0., -1), loc='lower left', ncol=4);\n",
    "\n",
    "format_ticks(ins_ax, y=False, rotate=False)\n",
    "ax.set_xlim(region[1], region[2])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4abafdcf",
   "metadata": {},
   "source": [
    "#### 2.3. Pileup plots at insulation boundaries\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d1d3698-200f-4a3f-b156-2b7c84a45ae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundaries = insulation_table[~np.isnan(insulation_table[f'boundary_strength_{selected_window_size}'])]\n",
    "strong_boundaries = boundaries[boundaries[f'is_boundary_{selected_window_size}']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d69e8a8f-4d23-460c-ba9e-38080054e6ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(strong_boundaries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b265612a-ab25-4268-a1b2-a6e1744bcdcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "flank = 300_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "922dd10f",
   "metadata": {},
   "outputs": [],
   "source": [
    "stack = cooltools.pileup(clr_10k, strong_boundaries, flank=flank, nproc=1)\n",
    "\n",
    "# Aggregate. Note that some pixels might be converted to NaNs after IC, thus we aggregate by nanmean: \n",
    "mtx = np.nanmean(stack, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc0e9352-0b1b-4905-af3f-b5d1e3e2ee4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "norm = mpl.colors.LogNorm(vmin=1e-3, vmax=0.02)\n",
    "plt.imshow(\n",
    "    mtx,\n",
    "    norm=norm,\n",
    "    cmap='fall',\n",
    "    interpolation='none')\n",
    "\n",
    "plt.colorbar(label = 'log10 mean ICed Hi-C')\n",
    "ticks_pixels = np.linspace(0, flank*2//resolution,5)\n",
    "ticks_kbp = ((ticks_pixels-ticks_pixels[-1]/2)*resolution//1000).astype(int)\n",
    "plt.xticks(ticks_pixels, ticks_kbp)\n",
    "plt.yticks(ticks_pixels, ticks_kbp)\n",
    "plt.xlabel('relative position, kbp')\n",
    "plt.ylabel('relative position, kbp')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4074a764-8af8-4032-8e0b-7fde57234f7a",
   "metadata": {},
   "source": [
    "▸ **Practice task.** \n",
    "\n",
    "Create the same pileups for Rad21 degradation. What are the changes?\n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c68563a-4c24-4ba9-9b84-ab2e5008afd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 10_000  \n",
    "clr_10k_rad21 = cooler.Cooler(...) # Open cooler file with dRad21 degron (IAA-treated) data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fffe9aa3-23f7-4a81-8880-3c554951e849",
   "metadata": {},
   "outputs": [],
   "source": [
    "stack = cooltools.pileup(clr_10k_rad21, strong_boundaries, flank=flank, nproc=1)\n",
    "mtx = np.nanmean(stack, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c453f52c-eb32-4212-96a9-38cd8f62dfa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "norm = mpl.colors.LogNorm(vmin=1e-3, vmax=0.02)\n",
    "plt.imshow(\n",
    "    mtx,\n",
    "    norm=norm,\n",
    "    cmap='fall',\n",
    "    interpolation='none')\n",
    "\n",
    "plt.colorbar(label = 'log10 mean ICed Hi-C')\n",
    "ticks_pixels = np.linspace(0, flank*2//resolution,5)\n",
    "ticks_kbp = ((ticks_pixels-ticks_pixels[-1]/2)*resolution//1000).astype(int)\n",
    "plt.xticks(ticks_pixels, ticks_kbp)\n",
    "plt.yticks(ticks_pixels, ticks_kbp)\n",
    "plt.xlabel('relative position, kbp')\n",
    "plt.ylabel('relative position, kbp')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d7f4ca",
   "metadata": {},
   "source": [
    "### 3. Compartments\n",
    "\n",
    "Compartments is another prominent feature of Hi-C/Micro-C. It involves long-range interactions of genomic loci and can be seen as a checkerboard pattern in maps. \n",
    "\n",
    "For compartmental calling, we perform eigendecomposition of Hi-C matrices.\n",
    "\n",
    "Since the orientation of eigenvectors is determined up to a sign, the convention for Hi-C data anaylsis is to orient eigenvectors to be positively correlated with a binned profile of GC content as a 'phasing track'. \n",
    "\n",
    "In humans and mice, GC content is useful for phasing because it typically has a strong correlation at the 100kb-1Mb bin level with the eigenvector. In other organisms, other phasing tracks have been used to orient\n",
    "eigenvectors from Hi-C data. \n",
    "\n",
    "\n",
    "\n",
    "<img src=\"https://www.biorxiv.org/content/biorxiv/early/2022/11/01/2022.10.31.514564/F2.large.jpg\"\n",
    "     alt=\"Compartments and saddles\"\n",
    "     style=\"float: left; margin-right: 10px;\" />\n",
    "     \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c1c9ccf-aac6-484e-b885-26eea183abeb",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### 3.1. Calculating compartmentalization in cis\n",
    "\n",
    "We first load the Hi-C data at 100 kb resolution. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef651481-bc4f-45e8-b86b-b7a44aa57670",
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 100_000\n",
    "clr_100k = cooler.Cooler(f\"{cool_file_ut}::resolutions/{resolution}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dfd8097",
   "metadata": {},
   "source": [
    "Download mouse genome and create the phasing track out of GC-content: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86080e9b-f887-4201-9f92-3989ed3d39ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "## fasta sequence is required for calculating binned profile of GC conent\n",
    "if not os.path.isfile('./mm10.fa'):\n",
    "    ## note downloading a ~1Gb file can take a minute\n",
    "    subprocess.call('wget https://hgdownload.cse.ucsc.edu/goldenpath/mm10/bigZips/mm10.fa.gz', shell=True,\n",
    "                    stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT)\n",
    "    subprocess.call('gunzip mm10.fa.gz', shell=True,\n",
    "                    stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "494178c0-c479-4fbf-80ef-d8aa76ac455b",
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = clr_100k.bins()[:]\n",
    "mm10_genome = bioframe.load_fasta('./mm10.fa')\n",
    "\n",
    "## note the next command may require installing pysam\n",
    "gc_cov = bioframe.frac_gc(bins[['chrom', 'start', 'end']], mm10_genome)\n",
    "gc_cov.to_csv('mm10_gc_cov_100kb.tsv',index=False,sep='\\t')\n",
    "display(gc_cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab666409-96ca-458b-a984-d5ab21529c6f",
   "metadata": {},
   "source": [
    "Cooltools also allows a view to be passed for eigendecomposition to limit to a certain set of regions. The following code creates the simplest view, of the two chromosomes in this cooler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2580405-2f1a-4433-a533-35fc2be91081",
   "metadata": {},
   "outputs": [],
   "source": [
    "view_df = pd.DataFrame({'chrom': clr_100k.chromnames,\n",
    "                        'start': 0,\n",
    "                        'end': clr_100k.chromsizes.values,\n",
    "                        'name': clr_100k.chromnames}\n",
    "                      )\n",
    "display(view_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d90e12d-39cc-4a5a-bbef-0940870f8deb",
   "metadata": {},
   "source": [
    "To capture the pattern of compartmentalization within-chromosomes, in cis, cooltools `eigs_cis` first removes\n",
    "the dependence of contact frequency by distance, and then performs eigenedecompostion. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d444c23d-58d6-4084-a4f1-411f76785e09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# obtain first 3 eigenvectors\n",
    "cis_eigs = cooltools.eigs_cis(clr_100k, \n",
    "                        gc_cov,\n",
    "                        view_df=view_df, \n",
    "                        n_eigs=3,\n",
    "                        )\n",
    "\n",
    "# cis_eigs[0] returns eigenvalues, here we focus on eigenvectors\n",
    "eigenvector_track = cis_eigs[1][['chrom', 'start', 'end', 'E1']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "989fe312-9f1c-499f-9fed-ef39f9c913fe",
   "metadata": {},
   "source": [
    "Plotting the first eigenvector next to the Hi-C map allows us to see how this captures the plaid pattern. \n",
    "\n",
    "To better visualize this relationship, we overlay the map with a binary segmentation of the eigenvector. Eigenvectors can be segmented by a variety of methods. The simplest segmentation, shown here, is to simply binarize eigenvectors, and term everything above zero the \"A-compartment\" and everything below 0 the \"B-compartment\".\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3125d232-fa49-4ad8-9484-3ac4804cf816",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "f, ax = plt.subplots(\n",
    "    figsize=(15, 10),\n",
    ")\n",
    "\n",
    "norm = LogNorm(vmax=0.1)\n",
    "\n",
    "im = ax.matshow(\n",
    "    clr_100k.matrix()[:], \n",
    "    norm=norm,  \n",
    "    cmap='fall'\n",
    "); \n",
    "plt.axis([0, 700, 700, 0])\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.01)\n",
    "plt.colorbar(im, cax=cax, label='corrected frequencies');\n",
    "ax.set_ylabel('chr1:0-100Mb')\n",
    "ax.xaxis.set_visible(False)\n",
    "\n",
    "ax1 = divider.append_axes(\"top\", size=\"20%\", pad=0.25, sharex=ax)\n",
    "weights = clr_100k.bins()[:]['weight'].values\n",
    "ax1.plot([0, 700], [0, 0], 'k',lw=0.25)\n",
    "ax1.plot( eigenvector_track['E1'].values, label='E1')\n",
    "\n",
    "ax1.set_ylabel('E1')\n",
    "ax1.set_xticks([]);\n",
    "\n",
    "\n",
    "for i in np.where(np.diff( (cis_eigs[1]['E1']>0).astype(int)))[0]:\n",
    "    ax.plot([i, 700], [i, i], 'k', lw=0.5)\n",
    "    ax.plot([i, i], [0, i], 'k', lw=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9990ccd6",
   "metadata": {},
   "source": [
    "#### 3.2. Compare compartments for untreated and Rad21 degradation\n",
    "\n",
    "▸ **Practice task.** \n",
    "\n",
    "Repeat the compartmental calling for Rad21 degradation. Do you see the difference? \n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions. \n",
    "\n",
    "Hint: you don't need to re-calculate the GC content for that. Just run cis_eig on the cool file with Rad21 degradation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eef5e774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here:\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b53e3e18",
   "metadata": {},
   "source": [
    "### 4. Scaling plots\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39ac9e6c-8603-4ed7-98e0-b0e573bf527b",
   "metadata": {},
   "source": [
    "In Hi-C maps, contact frequency decreases very strongly with **genomic separation** (also referred to as **genomic distance**). In the Hi-C field, this decay is often interchangeably referred to as the: \n",
    "\n",
    "* **expected** because one \"expects\" a certain average contact frequency at a given genomic separation\n",
    "* **scaling** which is borrowed from the polymer physics literature\n",
    "* **P(s) curve** contact *probability*, *P*, as a function of genomic *separation*, *s*. \n",
    "\n",
    "The rate of decay of contacts with genomic separation reflects the polymeric nature of chromosomes and can tell us about the global folding patterns of the genome. \n",
    "\n",
    "This decay has been observed to vary through the cell cycle, across cell types, and after degredation of structural maintenance of chromosomes complexes (SMCs) in both interphase and mitosis. \n",
    "\n",
    "In this chapter, we will: \n",
    "* calculate the P(s) of a given cooler\n",
    "* plot the P(s) curve\n",
    "* smooth the P(s) curve with logarithmic binning\n",
    "* plot the derivative of P(s)\n",
    "\n",
    "We will analyse the data with Rad21 and Wapl degratation in mouse from [Hsieh et al. 2022](https://www.nature.com/articles/s41588-022-01223-8), as above.\n",
    "\n",
    "Again, these datasets are available through cooltools interface, and we will use just two chromosomes for fast download and computations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69782420",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load a Hi-C map at a 1kb resolution from a cooler file.\n",
    "resolution = 1000 # note this might be slightly slow on a laptop\n",
    "                  # and could be lowered to 10kb for increased speed\n",
    "clr = cooler.Cooler(f'{cool_file_ut}::/resolutions/{resolution}')\n",
    "clr_rad21 = cooler.Cooler(f'{cool_file_rad21}::/resolutions/{resolution}')\n",
    "clr_wapl = cooler.Cooler(f'{cool_file_wapl}::/resolutions/{resolution}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c31cd23",
   "metadata": {},
   "source": [
    "Load mouse genome:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f47bdd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### If you have internet access:\n",
    "# Use bioframe to fetch the genomic features from the UCSC.\n",
    "mm10_chromsizes = bioframe.fetch_chromsizes('mm10')\n",
    "mm10_chromsizes = bioframe.make_viewframe(mm10_chromsizes)\n",
    "\n",
    "# select only those chromosomes available in cooler\n",
    "mm10_chromsizes = mm10_chromsizes[mm10_chromsizes.chrom.isin(clr.chromnames)].reset_index(drop=True)\n",
    "display(mm10_chromsizes.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e345057-f7b0-4aa9-88e1-c9f5ea1f4339",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Alternative if you're running locally offline:\n",
    "mm10_chromsizes = bioframe.make_viewframe( clr.chromsizes )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dabc1ae5-90d4-4568-969e-1a11ac7b0cd6",
   "metadata": {},
   "source": [
    "#### 4.1. Calculate the P(s) curves"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a7175ee-ed9e-40d2-ad02-fc7d91f78e4b",
   "metadata": {},
   "source": [
    "To calculate the average contact frequency as a function of genomic separation, we use the fact that each diagonal of a Hi-C map records contacts between loci separated by the same genomic distance. For example, the 3rd diagonal of our matrix contains contacts between loci separated by 3-4kb (note that diagonals are 0-indexed). Thus, we calculate the average contact frequency, *P(s)*, at a given genomic distance, *s*, as the average value of all pixels of the corresponding diagonal. This operation is performed by `cooltools.expected_cis`.\n",
    "\n",
    "Note that for human genome, we would calculate the *P(s)* separately for each chromosomal **arm**, by providing  `hg38_arms` as a `view_df`. This way we will ignore contacts accross the centromere, which is generally a good idea, since such contacts have a slightly different decay versus genomic separation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06b20fb0-9338-4e0a-b669-33af20143221",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cvd == contacts-vs-distance\n",
    "cvd = cooltools.expected_cis(\n",
    "    clr=clr,\n",
    "    view_df=mm10_chromsizes,\n",
    "    smooth=False,\n",
    "    aggregate_smoothed=False,\n",
    "    nproc=1 #if you do not have multiple cores available, set to 1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a395c81-9d5b-43e5-99ea-0e30bae33359",
   "metadata": {},
   "source": [
    "This function calculates average contact frequency for raw and normalized interactions ( `count.avg` and `balanced.avg`) for each diagonal and each regions in the `hg38_arms` of a Hi-C map. It aslo keeps the sum of raw and normalized interaction counts (`count.sum` and `balanced.sum`) as well as the number of valid (i.e. non-masked) pixels at each diagonal, `n_valid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3559058a-22d7-4c04-b886-059dcfb95ce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(cvd.head(4))\n",
    "display(cvd.tail(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28dbfc7e-ad4f-412a-97b0-dd1f531a4383",
   "metadata": {},
   "source": [
    "Note that the data from the first couple of diagonals are masked. This is done intentionally, since interactions at these diagonals (very short-ranged) are contaminated by non-informative Hi-C byproducts - dangling ends and self-circles. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "102ba897",
   "metadata": {},
   "source": [
    "Make sure also to convert distance to basepairs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d205e67f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cvd['s_bp'] = cvd['dist'] * resolution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b8bc2ec-cd04-4682-a2e1-711cfffbd975",
   "metadata": {},
   "source": [
    "#### 4.2. Plot the P(s) curve"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "509df096-aefe-4f3d-b1aa-4e406d16edcc",
   "metadata": {},
   "source": [
    "Time to plot *P(s)* !\n",
    "\n",
    "The first challenge  is that Hi-C has a very wide dynamic range. Hi-C probes genomic separations ranging from 100s to 100,000,000s of basepairs and contact frequencies also tend to span many orders of magnitude.\n",
    "\n",
    "Plotting such data in the linear scale would reveal only a part of the whole picture. Instead, we typically switch to double logarithmic (aka log-log) plots, where the x and y coordinates vary by orders of magnitude."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a61545c-c7a5-4745-b47d-c9d976b4d6e2",
   "metadata": {},
   "source": [
    "With the flags used above, `expected_cis()` does not smooth or aggregate across regions. This can lead to noisy P(s) curves for each region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0321bf87-df0b-4fd5-b919-ae2bf1a98c6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1,1)\n",
    "\n",
    "for region in mm10_chromsizes['name']:\n",
    "    ax.loglog(\n",
    "        cvd['s_bp'].loc[cvd['region1']==region],\n",
    "        cvd['balanced.avg'].loc[cvd['region1']==region],\n",
    "        alpha=0.7,\n",
    "        linewidth=1,\n",
    "    )\n",
    "    ax.set(\n",
    "        xlabel='separation, bp', \n",
    "        ylabel='IC contact frequency')\n",
    "    ax.set_aspect(1.0)\n",
    "    ax.grid(lw=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd345eaf-4782-4ea7-81c1-89327ed78a85",
   "metadata": {},
   "source": [
    "#### 4.3. Smoothing & aggregating P(s) curves"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ef00585-4176-4438-b4fd-9f7d827f7f1b",
   "metadata": {},
   "source": [
    "We notices how each chromosome gets too noisy signal at its ends. We can improve that by performing aggregation and smoothing of interactions.\n",
    "\n",
    "We can pass flags to `expected_cis()` that return smoothed and aggregated columns for futher analysis (which are on by default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9baa2053",
   "metadata": {},
   "outputs": [],
   "source": [
    "cvd_smooth_agg = cooltools.expected_cis(\n",
    "    clr=clr,\n",
    "    view_df=mm10_chromsizes,\n",
    "    smooth=True,\n",
    "    aggregate_smoothed=True,\n",
    "    nproc=1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "637798b1-5570-4b68-910d-81d8b37fafde",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(cvd_smooth_agg.head(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20473d6e-b51d-4b95-8404-d6e97d357ba9",
   "metadata": {},
   "source": [
    "The `balanced.avg.smoothed.agg` is averaged across regions, and shows the same exact curve for each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b7b55f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "cvd_smooth_agg['s_bp'] = cvd_smooth_agg['dist']* resolution\n",
    "cvd_smooth_agg['balanced.avg.smoothed.agg'].loc[cvd_smooth_agg['dist'] < 2] = np.nan\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3eb7391-7357-4fc3-8f0f-8a8a973a4719",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1,1)\n",
    "\n",
    "for region in mm10_chromsizes['chrom']:\n",
    "    ax.loglog(\n",
    "        cvd_smooth_agg['s_bp'].loc[cvd_smooth_agg['region1']==region],\n",
    "        cvd_smooth_agg['balanced.avg.smoothed.agg'].loc[cvd_smooth_agg['region1']==region],\n",
    "    )\n",
    "    ax.set(\n",
    "        xlabel='separation, bp', \n",
    "        ylabel='IC contact frequency')\n",
    "    ax.set_aspect(1.0)\n",
    "    ax.grid(lw=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da125b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat the same for rad21 and wapl degrons: \n",
    "\n",
    "## Your code goes here: \n",
    "\n",
    "cvd_smooth_agg_rad21 = cooltools.expected_cis(\n",
    "    clr=clr_rad21,\n",
    "    view_df=mm10_chromsizes,\n",
    "    smooth=True,\n",
    "    aggregate_smoothed=True,\n",
    "    nproc=1 # change to larger number of cores if they are available\n",
    ")\n",
    "\n",
    "cvd_smooth_agg_wapl = cooltools.expected_cis(\n",
    "    clr=clr_wapl,\n",
    "    view_df=mm10_chromsizes,\n",
    "    smooth=True,\n",
    "    aggregate_smoothed=True,\n",
    "    nproc=1 # change to larger number of cores if they are available\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3892bb9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Add resolutions: \n",
    "cvd_smooth_agg_rad21['s_bp'] = cvd_smooth_agg_rad21['dist'] * resolution\n",
    "cvd_smooth_agg_rad21.loc[cvd_smooth_agg_rad21['dist'] < 2, 'balanced.avg.smoothed.agg'] = np.nan\n",
    "cvd_smooth_agg_wapl['s_bp'] = cvd_smooth_agg_wapl['dist'] * resolution\n",
    "cvd_smooth_agg_wapl.loc[cvd_smooth_agg_wapl['dist'] < 2, 'balanced.avg.smoothed.agg'] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8d5133e",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your code goes here: \n",
    "display(cvd_smooth_agg_rad21.head(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cabc5519-55f3-44f4-99a4-1b62649fa54e",
   "metadata": {},
   "source": [
    "#### 4.4. Plot the smoothed P(s) curve and its derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f4978b8-85ed-4444-b986-3ae27e3b5339",
   "metadata": {},
   "source": [
    "Logbin-smoothing of P(s) reduces the \"fanning\" at longer s and enables us to plot the derivative of the P(s) curve in the log-log space. This derivative is extremely informative, as it can be compared to predictions from various polymer models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15525348-ea94-4ddc-9b11-60ce1952ac82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just take a single value for each genomic separation\n",
    "cvd_merged = cvd_smooth_agg.drop_duplicates(subset=['dist'])[['s_bp', 'balanced.avg.smoothed.agg']]\n",
    "\n",
    "# Calculate derivative in log-log space\n",
    "der = np.gradient(np.log(cvd_merged['balanced.avg.smoothed.agg']),\n",
    "                  np.log(cvd_merged['s_bp']))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daa734bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(\n",
    "    figsize=(6.5,13),\n",
    "    nrows=2, \n",
    "    gridspec_kw={'height_ratios':[6,2]}, \n",
    "    sharex=True)\n",
    "\n",
    "\n",
    "ax = axs[0]\n",
    "ax.loglog(\n",
    "    cvd_merged['s_bp'],\n",
    "    cvd_merged['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"Untreated\",\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    ylabel='IC contact frequency',\n",
    "    xlim=(1e3,1e8)\n",
    ")\n",
    "ax.set_aspect(1.0)\n",
    "ax.grid(lw=0.5)\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "ax = axs[1]\n",
    "ax.semilogx(\n",
    "    cvd_merged['s_bp'],\n",
    "    der,\n",
    "    alpha=0.5\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    xlabel='separation, bp', \n",
    "    ylabel='slope')\n",
    "\n",
    "ax.grid(lw=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578591e0",
   "metadata": {},
   "source": [
    "▸ **Practice task.** \n",
    "\n",
    "Calculate contact-versus-distance and derivatives for Rad21 and Wapl degrons.  \n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b50a7d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat the same for rad21 and wapl:\n",
    "\n",
    "# Repeat for degrons: \n",
    "cvd_merged_rad21 = cvd_smooth_agg_rad21.drop_duplicates(subset=['dist'])[['s_bp', 'balanced.avg.smoothed.agg']]\n",
    "der_rad21 = np.gradient(\n",
    "    np.log(cvd_merged_rad21['balanced.avg.smoothed.agg']),\n",
    "    np.log(cvd_merged_rad21['s_bp']))\n",
    "\n",
    "cvd_merged_wapl = cvd_smooth_agg_wapl.drop_duplicates(subset=['dist'])[['s_bp', 'balanced.avg.smoothed.agg']]\n",
    "der_wapl = np.gradient(\n",
    "    np.log(cvd_merged_wapl['balanced.avg.smoothed.agg']),\n",
    "    np.log(cvd_merged_wapl['s_bp']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d27d5a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(\n",
    "    figsize=(6.5,13),\n",
    "    nrows=2, \n",
    "    gridspec_kw={'height_ratios':[6,2]}, \n",
    "    sharex=True)\n",
    "\n",
    "\n",
    "ax = axs[0]\n",
    "ax.loglog(\n",
    "    cvd_merged['s_bp'],\n",
    "    cvd_merged['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"Untreated\",\n",
    ")\n",
    "ax.loglog(\n",
    "    cvd_merged_rad21['s_bp'],\n",
    "    cvd_merged_rad21['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"Rad21 Degron\",\n",
    ")\n",
    "ax.loglog(\n",
    "    cvd_merged_wapl['s_bp'],\n",
    "    cvd_merged_wapl['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"WAPL Degron\",\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    ylabel='IC contact frequency',\n",
    "    xlim=(1e3,1e8)\n",
    ")\n",
    "ax.set_aspect(1.0)\n",
    "ax.grid(lw=0.5)\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "ax = axs[1]\n",
    "ax.semilogx(\n",
    "    cvd_merged['s_bp'],\n",
    "    der,\n",
    "    alpha=0.5\n",
    ")\n",
    "ax.semilogx(\n",
    "    cvd_merged_rad21['s_bp'],\n",
    "    der_rad21,\n",
    "    alpha=0.5\n",
    ")\n",
    "ax.semilogx(\n",
    "    cvd_merged_wapl['s_bp'],\n",
    "    der_wapl,\n",
    "    alpha=0.5\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    xlabel='separation, bp', \n",
    "    ylabel='slope')\n",
    "\n",
    "ax.grid(lw=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2574fae8-b4d3-45e6-9fbd-96ed449ac354",
   "metadata": {},
   "source": [
    "▸ **Practice task (optional).** \n",
    "\n",
    "Compare contact-versus-distance and derivatives for Rad21, Wapl and CTCF degrons. \n",
    "\n",
    "Does the average loop size changes when CTCF is degraded? \n",
    "\n",
    "Hint: Read dCTCF-IAA data and add to the visualization.\n",
    "\n",
    "⏸ 10-15 minutes: take your time and ask any questions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abdd5ac3-df75-4405-8590-e3030733be3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read cooler for dCTCF_IAA\n",
    "clr_ctcf = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76841138-5a8c-43b9-b7ad-9f508e81f03f",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Generated table of expected:\n",
    "cvd_smooth_agg_ctcf = cooltools.expected_cis(\n",
    "    clr=clr_ctcf,\n",
    "    view_df=mm10_chromsizes,\n",
    "    smooth=True,\n",
    "    aggregate_smoothed=True,\n",
    "    nproc=1 # change to larger number of cores if they are available\n",
    ")\n",
    "\n",
    "cvd_smooth_agg_ctcf['s_bp'] = cvd_smooth_agg_ctcf['dist'] * resolution\n",
    "cvd_smooth_agg_ctcf.loc[cvd_smooth_agg_ctcf['dist'] < 2, 'balanced.avg.smoothed.agg'] = np.nan\n",
    "\n",
    "cvd_merged_ctcf = cvd_smooth_agg_ctcf.drop_duplicates(subset=['dist'])[['s_bp', 'balanced.avg.smoothed.agg']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbcc3afa-7456-43a3-a2bd-5b56860e9e5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Calculate derivative:\n",
    "der_ctcf = np.gradient(\n",
    "    np.log(cvd_merged_ctcf['balanced.avg.smoothed.agg']),\n",
    "    np.log(cvd_merged_ctcf['s_bp']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17bc6764-95b9-4330-9551-45c3ca985d19",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot:\n",
    "f, axs = plt.subplots(\n",
    "    figsize=(6.5,13),\n",
    "    nrows=2, \n",
    "    gridspec_kw={'height_ratios':[6,2]}, \n",
    "    sharex=True)\n",
    "\n",
    "\n",
    "ax = axs[0]\n",
    "ax.loglog(\n",
    "    cvd_merged['s_bp'],\n",
    "    cvd_merged['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"Untreated\",\n",
    ")\n",
    "ax.loglog(\n",
    "    cvd_merged_rad21['s_bp'],\n",
    "    cvd_merged_rad21['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"Rad21 Degron\",\n",
    ")\n",
    "ax.loglog(\n",
    "    cvd_merged_wapl['s_bp'],\n",
    "    cvd_merged_wapl['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"WAPL Degron\",\n",
    ")\n",
    "ax.loglog(\n",
    "    cvd_merged_ctcf['s_bp'],\n",
    "    cvd_merged_ctcf['balanced.avg.smoothed.agg'],\n",
    "    '-',\n",
    "    markersize=5,\n",
    "    label=\"CTCF Degron\",\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    ylabel='IC contact frequency',\n",
    "    xlim=(1e3,1e8)\n",
    ")\n",
    "ax.set_aspect(1.0)\n",
    "ax.grid(lw=0.5)\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "ax = axs[1]\n",
    "ax.semilogx(\n",
    "    cvd_merged['s_bp'],\n",
    "    der,\n",
    "    alpha=0.5\n",
    ")\n",
    "ax.semilogx(\n",
    "    cvd_merged_rad21['s_bp'],\n",
    "    der_rad21,\n",
    "    alpha=0.5\n",
    ")\n",
    "ax.semilogx(\n",
    "    cvd_merged_wapl['s_bp'],\n",
    "    der_wapl,\n",
    "    alpha=0.5\n",
    ")\n",
    "ax.semilogx(\n",
    "    cvd_merged_ctcf['s_bp'],\n",
    "    der_ctcf,\n",
    "    alpha=0.5\n",
    ")\n",
    "\n",
    "ax.set(\n",
    "    xlabel='separation, bp', \n",
    "    ylabel='slope')\n",
    "\n",
    "ax.grid(lw=0.5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hic_open2c_workshop",
   "language": "python",
   "name": "hic_open2c_workshop"
  },
  "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}