What_we_see_in_Hi-C_4Jun2023.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "97dced82",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "Welcome to the workshop \"What can we see in Hi-C?\" by Leonid Mirny, Aleksandra Galitsyna and Emily Navarrete. \n",
    "Your assistants for this workshop are Aleksandra (Alex) Galitsyna and Emily Navarrete. \n",
    "For any questions, feel free to reach out to galitsyn at mit dot edu\n",
    "\n",
    "\n",
    "## Workshop outline\n",
    "\n",
    "\n",
    "1. [Qualitative analysis](Part-I.-Qualitative-analysis-of-Hi-C-maps)\n",
    "\n",
    "    - [Databases for Hi-C data](#1.-Databases-for-Hi-C-data)\n",
    "    \n",
    "    - [Interactive browsers for Hi-C](#2.-Interactive-browsers-for-Hi-C)\n",
    "\n",
    "2. [Quantitative analysis](#Part-II.-Quantitative-analysis-of-Hi-C-maps)\n",
    "\n",
    "    - [Load and navigate Hi-C in Python](#1.-Load-and-navigate-Hi-C-in-Python)\n",
    "    \n",
    "    - [Insulation](#2.-Insulation)\n",
    "    \n",
    "    - [Compartments](#3.-Compartments)\n",
    "    \n",
    "    - [Scaling plots](#4.-Scaling-plots)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34a11a42",
   "metadata": {},
   "source": [
    "## Part I. Qualitative analysis of Hi-C maps\n",
    "\n",
    "[go top](#Workshop-outline)\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 withinteractive browser resgen, loading the data from 4DN data portal and ENCODE. \n",
    "\n",
    "\n",
    "### 1. Databases for Hi-C 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",
    "[go top](#Workshop-outline)\n",
    "\n",
    "Let's now explore Hi-C 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 **HFF6c and h1-ESC cells probed by Micro-C** (high-resolution version of Hi-C). \n",
    "\n",
    "\n",
    "#### 2.1. Access resgen server.\n",
    "\n",
    "\n",
    "- Go to resgen: https://resgen.io/\n",
    "\n",
    "- Register (sign up) there with your e-mail.\n",
    "\n",
    "- Create your own test porject or join the project: https://resgen.io/mirnylab/workshop2023/\n",
    "\n",
    "- Enjoy the creation!\n",
    "\n",
    "\n",
    "#### 2.2. Create and save your view.\n",
    "\n",
    "- Open your test project.\n",
    "\n",
    "- **Load h1-ESC Micro-C dataset to the central panel of your view**. \n",
    "\n",
    "- What is the reference genome of this dataset? Hint: it's hg38, displyed in the bottom right corner. \n",
    "\n",
    "- Look up for hg38 reference genome coordinates in the left upper corner search box of resgen. Add them to the upper panel of the view. \n",
    "\n",
    "- Add gene annotation for the hg38 in the same manner. \n",
    "\n",
    "- Search for CTCF motifs available for hg38 genome. Add them to the upper panel.\n",
    "\n",
    "- Add CTCF ChIP-Seq coverage track for h1-ESC (bigwig format) from the available datasets on the left. \n",
    "\n",
    "- Learn how to zoom in, zoom out, and move around the view. Spot interesting TADs, compartments, dots in the dataset.\n",
    "\n",
    "- Now, it's time to **add HFF6c cells for the comparison**! Multiply the view (copy symbol on the right of the view panel). And load Micro-C and CTCF ChIP-Seq coverage track for HFF6C instead of h1ESC of the right panel. You may close the tracks that you don't need anymore. \n",
    "\n",
    "▸ **Practice task.**\n",
    "Select the region you like (it might be a gene or the genome locus, something with an interesting difference or effect). Save the view (lower right panel).\n",
    "\n",
    "\n",
    "⏸ *5-10 minutes: take your time and prepare the view you like.*\n",
    "\n",
    "#### 2.3. Add custom data from ENCODE (optional).\n",
    "\n",
    "ENCODE provides a vastful resource of ChIP-Seq and other epigenetics data. Let's try to add your favorime histone modification or transctiption factor binding to the view. \n",
    "\n",
    "- Go to [ENCODE project webpage](https://www.encodeproject.org/).\n",
    "\n",
    "- Navigate to \"Functional genomics\" page with available data table. \n",
    "\n",
    "- Select any ChIP-Seq for either transcription factor (TF) or histone modification) for either H1 or HFF cells (or two files for both). *Note that the genome assembly should be the same as our Micro-C (hg38)*\n",
    "\n",
    "- Open the link to the selected ChIP-Seq dataset(s) and navigate to \"File details\". \n",
    "You need the link to the **bigWig file \"fold change over control\" or \"signal p-value\"**. Copy the link. \n",
    "\n",
    "- Go back to resgen. Press \"Add dataset\" and then \"URL\". \n",
    "Paste your copied URL from ENCODE and change the file name so that it contains: (i) the name of the cell line, (ii) the name of the target factor for ChIP-Seq, (iii) ENCODE ID of the file or total dataset. \n",
    "\n",
    "- Add your new dataset to resgen view. Note that you may add the same target factor for both h1-ESC and HFF, and you may do interesting comparison between the two. \n",
    "\n",
    "▸ **Practice task.** Do you see any association of the selected factor with the chromatin architecture? \n",
    "Find the region you like and save the view with it. \n",
    "\n",
    "Save the view (lower right panel) and put your name to the title of it.\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",
    "[go top](#Workshop-outline)\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",
    "![Open2C](https://avatars.githubusercontent.com/u/70977326?s=400&u=8c4db6ec3adde1b507f751f166b0b3465c4fd65a&v=4 \"Open2C\")\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 and preprints: 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 HFF6c dataset. You may repeat the same analysis for H1-ESC and compare the results, if you want.\n",
    "- We selected only two chromosomes (chr2 and chr17) to make the computations faster. "
   ]
  },
  {
   "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/GenomeArchitectureAndFunction2023/blob/main/What_we_see_in_Hi-C_4Jun2023.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 cooler==0.9.2 cooltools==0.5.4 bioframe==0.4.1 pysam\n",
    "# !pip install aiohttp\n",
    "# !pip install --pre higlass-python==1.0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6559f3ef-17f9-47bd-99ca-1aefab6c30d2",
   "metadata": {},
   "source": [
    "#### Plan B. Open one of MIT servers, launch terminal, git-clone and open the notebook\n",
    "(dependencies are pre-installed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6b699a-4994-4608-9942-642914b5c6ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !git clone https://github.com/mirnylab/GenomeArchitectureAndFunction2023.git"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35f72da",
   "metadata": {},
   "source": [
    "### 1. Load and navigate Hi-C in Python\n",
    "[go top](#Workshop-outline)\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.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. \n",
    "you need HFF and hESC Micro-C:"
   ]
  },
  {
   "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 HFF cells: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc5abee9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# download test data\n",
    "# this file is 145 Mb, and may take a few seconds to download\n",
    "\n",
    "data_dir = './data/'\n",
    "cool_file = cooltools.download_data(\"HFF_MicroC\", cache=True, data_dir=data_dir)    \n",
    "print(cool_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eba170ab",
   "metadata": {},
   "source": [
    "Repeat the same for h1-ESC cells: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67d490e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "cool_file_hesc = cooltools.download_data(\"hESC_MicroC\", cache=True, data_dir=data_dir)\n",
    "print(cool_file_hesc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0ddd04f",
   "metadata": {},
   "source": [
    "The file we just downloaded, test.mcool, contains Micro-C data from HFF cells for two chromosomes in a [multi-resolution mcool format](https://cooler.readthedocs.io/en/latest/schema.html?highlight=mcool#multi-resolution)."
   ]
  },
  {
   "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(f'./data/test.mcool')"
   ]
  },
  {
   "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'./data/test.mcool::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 hESC cooler: \n",
    "clr_hesc = cooler.Cooler(f'./data/test_hESC.mcool::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 subregions\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('chr17')\n",
    "\n",
    "ax = axs[1]\n",
    "im = ax.matshow(\n",
    "    data_chr17,\n",
    "    vmax=5e-2, \n",
    "    extent=(0, clr.chromsizes['chr17'], clr.chromsizes['chr17'], 0)\n",
    ");\n",
    "plt.colorbar(im, ax=ax ,fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title('chr17', 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 = ('chr17', 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'chr17:{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('chr17')\n",
    "\n",
    "ax = axs[1]\n",
    "im = ax.matshow(\n",
    "    data_chr17,\n",
    "    norm=norm,\n",
    "    extent=(0, clr.chromsizes['chr17'], clr.chromsizes['chr17'], 0)\n",
    ");\n",
    "plt.colorbar(im, ax=ax ,fraction=0.046, pad=0.04, label='norm counts');\n",
    "ax.set_title('chr17', 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 = ('chr17', 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'./data/test.mcool::resolutions/10000')\n",
    "\n",
    "region = 'chr17:30,000,000-35,000,000'\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 HFF and h1ESC cells in Python\n",
    "\n",
    "▸ **Practice task.** \n",
    "\n",
    "Select the region you like (might be another one or the region you explored in resgen, if that happened to be in chromosome 2 or 17). \n",
    "\n",
    "Visualize that region with the color scheme and scale that you find instructuve and interesting. \n",
    "\n",
    "Repeat the same for both HFF and h1ESC 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": [
    "# Your code goes here:\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74944c26-7dd2-47ff-892c-04bb68328444",
   "metadata": {},
   "source": [
    "#### 1.6. HiGlass for Python\n",
    "\n",
    "[go top](#Workshop-outline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8684bedd-b885-4e87-849b-7b922d6d12cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import higlass as hg\n",
    "\n",
    "# Configure remote data sources (tilesets)\n",
    "tileset1 = hg.remote(\n",
    "    uid=\"R-Xq4_8GT-uUwKTt8kNuCA\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"Krietenstein et al. (2021) Micro-C for HFF1\",\n",
    ")\n",
    "\n",
    "tileset2 = hg.remote(\n",
    "    uid=\"O0q8If8qQwWqCv7egiBw6A\",\n",
    "    server=\"https://resgen.io/api/v1\",\n",
    "    name=\"Krietenstein et al. (2021) Micro-C for H1-ESC\",\n",
    ")\n",
    "\n",
    "# Create a HeatmapTrack for each tileset\n",
    "track1 = tileset1.track(\"heatmap\")\n",
    "track2 = tileset2.track(\"heatmap\")\n",
    "\n",
    "# Create two independent Views, one for each heatmap\n",
    "view1 = hg.view(track1, width=6)\n",
    "view2 = hg.view(track2, width=6)\n",
    "\n",
    "# Lock zoom & location for each View\n",
    "view_lock = hg.lock(view1, view2)\n",
    "\n",
    "# Concatenate Views side-by-side, and apply synchronization lock\n",
    "(view1 | view2).locks(view_lock)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1222348b",
   "metadata": {},
   "source": [
    "### 2. Insulation\n",
    "\n",
    "[go top](#Workshop-outline)\n"
   ]
  },
  {
   "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 Hi-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 = 10000  \n",
    "clr = cooler.Cooler(f'./data/test.mcool::resolutions/{resolution}')\n",
    "windows = [3*resolution, 5*resolution, 10*resolution, 25*resolution]\n",
    "insulation_table = insulation(clr, windows, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb0ba2af",
   "metadata": {},
   "outputs": [],
   "source": [
    "### repeat the same for h1-esc cells\n",
    "# clr_hesc = cooler.Cooler(f'./data/test_hESC.mcool::resolutions/{resolution}')\n",
    "# windows = [3*resolution, 5*resolution, 10*resolution, 25*resolution]\n",
    "# insulation_table_hesc = insulation(clr, windows, verbose=True)"
   ]
  },
  {
   "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:1005]"
   ]
  },
  {
   "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 = ('chr2', start, end)\n",
    "norm = LogNorm(vmax=0.1, vmin=0.001)\n",
    "data = clr.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[['start', 'end']].mean(axis=1), \n",
    "            insul_region['log2_insulation_score_'+str(windows[0])],\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[['start', 'end']].mean(axis=1), \n",
    "                insul_region[f'log2_insulation_score_{res}'], \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": "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_{windows[0]}'], label=f'Window {windows[0]} bp')\n",
    "\n",
    "# Define weak and strong boundaries:\n",
    "boundaries = insul_region[~np.isnan(insul_region[f'boundary_strength_{windows[0]}'])]\n",
    "weak_boundaries = boundaries[~boundaries[f'is_boundary_{windows[0]}']]\n",
    "strong_boundaries = boundaries[boundaries[f'is_boundary_{windows[0]}']]\n",
    "\n",
    "# Plot the boundaries:\n",
    "ins_ax.scatter(weak_boundaries[['start', 'end']].mean(axis=1), \n",
    "            weak_boundaries[f'log2_insulation_score_{windows[0]}'], label='Weak boundaries')\n",
    "ins_ax.scatter(strong_boundaries[['start', 'end']].mean(axis=1), \n",
    "            strong_boundaries[f'log2_insulation_score_{windows[0]}'], 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. Compare insulating boundaries for HFF and h1ESC cells in Python\n",
    "\n",
    "▸ **Practice task.** \n",
    "\n",
    "Select the region you like (might be another one or the region you explored in resgen, if that happened to be in chromosome 2 or 17), or the same one as above.\n",
    "\n",
    "Visualize that region with insulation and boundaries. \n",
    "\n",
    "Repeat the same for both HFF and h1ESC cells. Do you see the difference? \n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "922dd10f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here:\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d7f4ca",
   "metadata": {},
   "source": [
    "### 3. Compartments\n",
    "    \n",
    "[go top](#Workshop-outline)\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 per-chromosome compartmentalization\n",
    "\n",
    "We first load the Hi-C data at 100 kbp resolution. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef651481-bc4f-45e8-b86b-b7a44aa57670",
   "metadata": {},
   "outputs": [],
   "source": [
    "clr = cooler.Cooler('./data/test.mcool::resolutions/100000')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dfd8097",
   "metadata": {},
   "source": [
    "Download human 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('./hg38.fa'):\n",
    "    ## note downloading a ~1Gb file can take a minute\n",
    "    subprocess.call('wget https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.fa.gz', shell=True,\n",
    "                    stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT)\n",
    "    subprocess.call('gunzip hg38.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.bins()[:]\n",
    "hg38_genome = bioframe.load_fasta('./hg38.fa');\n",
    "## note the next command may require installing pysam\n",
    "gc_cov = bioframe.frac_gc(bins[['chrom', 'start', 'end']], hg38_genome)\n",
    "gc_cov.to_csv('hg38_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.chromnames,\n",
    "                        'start': 0,\n",
    "                        'end': clr.chromsizes.values,\n",
    "                        'name': clr.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( \n",
    "                        clr, \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.matrix()[:], \n",
    "    norm=norm,  \n",
    "    cmap='fall'\n",
    "); \n",
    "plt.axis([0,500,500,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('chr2:0-50Mb')\n",
    "ax.xaxis.set_visible(False)\n",
    "\n",
    "ax1 = divider.append_axes(\"top\", size=\"20%\", pad=0.25, sharex=ax)\n",
    "weights = clr.bins()[:]['weight'].values\n",
    "ax1.plot([0,500], [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([0, 500],[i,i],'k',lw=0.5)\n",
    "    ax.plot([i,i],[0, 500],'k',lw=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9990ccd6",
   "metadata": {},
   "source": [
    "#### 3.2. Compare compartments for HFF and h1ESC cells in Python\n",
    "\n",
    "▸ **Practice task.** \n",
    "\n",
    "Repeat the compartmental calling for h1ESC cells. Do you see the difference? \n",
    "\n",
    "⏸ 5-10 minutes: take your time and ask any questions. \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",
    "[go top](#Workshop-outline)\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).\n",
    "\n",
    "These datasets are available through cooltools interface, and we will use just three chromosomes for fast download and computations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2de2e64-c500-4c32-bb3b-217527e7172f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# download test data\n",
    "# this file is 145 Mb, and may take a few seconds to download\n",
    "\n",
    "data_dir = './data/'\n",
    "cooltools.download_data(\"mESC_dRad21_UT\", cache=True, data_dir=data_dir);\n",
    "cooltools.download_data(\"mESC_dRad21_IAA\", cache=True, data_dir=data_dir);\n",
    "cooltools.download_data(\"mESC_dWapl_IAA\", cache=True, data_dir=data_dir);"
   ]
  },
  {
   "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('./data/dRAD21_UT.mm10.mapq_30.mcool::/resolutions/'+str(resolution))\n",
    "clr_rad21 = cooler.Cooler('./data/dRAD21_IAA.mm10.mapq_30.mcool::/resolutions/'+str(resolution))\n",
    "clr_wapl = cooler.Cooler('./data/dWAPL_IAA.mm10.mapq_30.mcool::/resolutions/'+str(resolution))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c31cd23",
   "metadata": {},
   "source": [
    "Load mouse genome:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f47bdd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Use bioframe to fetch the genomic features from the UCSC.\n",
    "mm10_chromsizes = bioframe.fetch_chromsizes('mm10')\n",
    "mm10_cens = bioframe.fetch_centromeres('mm10')\n",
    "# create a view with chromosome arms using chromosome sizes and definition of centromeres\n",
    "mm10_arms = bioframe.make_chromarms(mm10_chromsizes, mm10_cens)\n",
    "\n",
    "# select only those chromosomes available in cooler\n",
    "mm10_arms = mm10_arms[mm10_arms.chrom.isin(clr.chromnames)].reset_index(drop=True)\n",
    "display(mm10_arms.head())"
   ]
  },
  {
   "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_arms,\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_arms['name']:\n",
    "    ax.loglog(\n",
    "        cvd['s_bp'].loc[cvd['region1']==region],\n",
    "        cvd['balanced.avg'].loc[cvd['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": "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 == contacts-vs-distance\n",
    "cvd_smooth_agg = cooltools.expected_cis(\n",
    "    clr=clr,\n",
    "    view_df=mm10_arms,\n",
    "    smooth=True,\n",
    "    aggregate_smoothed=True,\n",
    "    nproc=1 #if you do not have multiple cores available, set to 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_arms['name']:\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_arms,\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_arms,\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)\n",
    "\n",
    "\n",
    "# f.savefig('Rad21-wapl-degron.pdf')\n",
    "# f.savefig('Rad21-wapl-degron.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00958921",
   "metadata": {},
   "source": [
    "▸ **Practice task.** (Optional) \n",
    "\n",
    "Calculate contact-versus-distnace for HFF and h1-ESC cells and compare the results to the mouse data. \n",
    "Note that it will require you to change the genome to hg38 and cooler files!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85968dff-d59d-4910-88c0-a9185e6fdd1c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "test",
   "language": "python",
   "name": "test"
  },
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}