CallingCompartments_and_Insulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "require([\"codemirror/keymap/sublime\", \"notebook/js/cell\", \"base/js/namespace\"],\n",
       "    function(sublime_keymap, cell, IPython) {\n",
       "        cell.Cell.options_default.cm_config.keyMap = 'sublime';\n",
       "    \n",
       "        var cells = IPython.notebook.get_cells();\n",
       "        for(var cl=0; cl< cells.length ; cl++){\n",
       "            cells[cl].code_mirror.setOption('keyMap', 'sublime');\n",
       "        }\n",
       "    }\n",
       ");\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "require([\"codemirror/keymap/sublime\", \"notebook/js/cell\", \"base/js/namespace\"],\n",
    "    function(sublime_keymap, cell, IPython) {\n",
    "        cell.Cell.options_default.cm_config.keyMap = 'sublime';\n",
    "    \n",
    "        var cells = IPython.notebook.get_cells();\n",
    "        for(var cl=0; cl< cells.length ; cl++){\n",
    "            cells[cl].code_mirror.setOption('keyMap', 'sublime');\n",
    "        }\n",
    "    }\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width: 90% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# change the cell width\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width: 90% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import os\n",
    "import cooler\n",
    "import cooltools\n",
    "import numpy as np\n",
    "from cooltools.api import eigdecomp\n",
    "import bioframe\n",
    "from pathlib import Path\n",
    "import multiprocess as mp\n",
    "\n",
    "from copy import copy\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.print_figure_kwargs={'facecolor' : \"w\"}\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### some of the input cooler files that we used to combine/recombine/subset etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ins_samples = {\n",
    "    \"Hap1-WT-combined.hg19\" : \"Hap1-WT-combined.mcool\",\n",
    "    \"PolII-IAA.hg19\" : \"PolII-IAA.hg19.mapq_30.1000.mcool\",\n",
    "    \"PolII-NT.hg19\" : \"PolII-NT.hg19.mapq_30.1000.mcool\",\n",
    "    \"all_CTCF_combined.hg19\" : \"all_combined.mcool\",\n",
    "    \"noCTCF_IAA_combined.hg19\" : \"noCTCF_IAA_combined.mcool\",\n",
    "    \n",
    "    \"all_CTCF_combined.hg19\" : \"all_combined.mcool\",\n",
    "    \"noCTCF_IAA_combined.hg19\" : \"noCTCF_IAA_combined.mcool\",\n",
    "    \n",
    "    \"CkoCT442_NT_pool.hg19\" : \"CkoCT442-NT-pool.mcool\",\n",
    "    \"CkoCT442_IAA_pool.hg19\" : \"CkoCT442-IAA-pool.mcool\",\n",
    "    # RAD21 degron for HCT from Rao et al 2017 ...\n",
    "    \"RAD21-IAA.hg19\" : \"RAD21-AID-IAA-6H.hg19.mapq_30.1000.mcool\",\n",
    "    \"RAD21-NT.hg19\" : \"RAD21-AID-NT.hg19.mapq_30.1000.mcool\",\n",
    "\n",
    "    # WAPL KO for HAP1 from Haarhuis 2017 ...\n",
    "    \"WAPL-KO1.hg19\" : \"Haarthuis-WAPL-KO1.14.hg19.mapq_30.1000.mcool\",\n",
    "    \"WAPL-KO3.hg19\" : \"Haarthuis-WAPL-KO3.3.hg19.mapq_30.1000.mcool\",\n",
    "    \"WAPL-WT.hg19\" : \"Haarthuis-HAP1-WT.hg19.mapq_30.1000.mcool\",\n",
    "\n",
    "    \"CkoC442-IAA48H-R1-T1.hg19\" : \"CkoC442-IAA48H-R1-T1__hg19.hg19.mapq_30.1000.mcool\",\n",
    "    \"CkoC442-NT-R1-T1.hg19\" : \"./wt_like_coolers_to_combine/CkoC442-NT-R1-T1__hg19.hg19.mapq_30.1000.mcool\",\n",
    "    \"CkoCT442-AAVS1sg2-4-NT-R1-T1.hg19\" : \"CkoCT442-AAVS1sg2-4-NT-R1-T1__hg19.hg19.mapq_30.1000.mcool\",\n",
    "    \"CkoCT442-DDX55sg2-B-NT-R1-T1.hg19\" : \"CkoCT442-DDX55sg2-B-NT-R1-T1__hg19.hg19.mapq_30.1000.mcool\",\n",
    "    \"G1-NT.hg19\" : \"G1-NT.hg19.1kb.mcool\",\n",
    "    \"G1-IAA.hg19\" : \"G1-IAA.hg19.1kb.mcool\",\n",
    "    \n",
    "    # Mutants time !!!\n",
    "    # Knock-outs - genes were deleted ...\n",
    "    \"mutControl-NT.hg19\":\"CkoCT442-AAVS1-NT-pool.mcool\",\n",
    "    \"mutControl-IAA.hg19\":\"CkoCT442-AAVS1-IAA-pool.mcool\",\n",
    "    \"mutDDX55-NT.hg19\":\"DDX55-clones-NT.hg19.mapq_30.1000.mcool\",\n",
    "    \"mutDDX55-IAA.hg19\":\"DDX55-clones-IAA.hg19.mapq_30.1000.mcool\",\n",
    "    \"mutTAF5L-NT.hg19\":\"TAF5L-clones-NT.hg19.mapq_30.1000.mcool\",\n",
    "    \"mutTAF5L-IAA.hg19\":\"TAF5L-clones-IAA.hg19.mapq_30.1000.mcool\",\n",
    "    # Knock-downs - genes were silenced/deplted (siRNA) ...\n",
    "    \"siControl-NT.hg19\":\"siCTRL-NT.hg19.mapq_30.1000.mcool\",\n",
    "    \"siControl-IAA.hg19\":\"siCTRL-IAA.hg19.mapq_30.1000.mcool\",\n",
    "    \"siDDX55-NT.hg19\":\"siDDX55-NT.hg19.mapq_30.1000.mcool\",\n",
    "    \"siDDX55-IAA.hg19\":\"siDDX55-IAA.hg19.mapq_30.1000.mcool\",\n",
    "    \"siTAF5L-NT.hg19\":\"siTAF5L-NT.hg19.mapq_30.1000.mcool\",\n",
    "    \"siTAF5L-IAA.hg19\":\"siTAF5L-IAA.hg19.mapq_30.1000.mcool\",\n",
    "    \n",
    "    # PlaB (splicing inhibition) two replicates pooled together - known to affect R-loops (mostly destabilize)\n",
    "    \"CtrlPlaB-NT.hg19\" : \"NT-hg19-combined-90000000.mcool\",\n",
    "    \"CtrlPlaB-IAA.hg19\" : \"IAA-hg19-combined-90000000.mcool\",\n",
    "    \"PlaB-NT.hg19\" : \"NT-PlaB-hg19-combined-90000000.mcool\",\n",
    "    \"PlaB-IAA.hg19\" : \"IAA-PlaB-hg19-combined-90000000.mcool\",\n",
    "    \n",
    "    ## Pool HAP1 WT, NO TIR1 and TIR1\n",
    "    # this is a special subsampling case for the standard CTCF-degron(+aux)\n",
    "    # but normalized, such that the number of reads, matches the NO-TIR1\n",
    "    \"Ctrl500M-noTIR1.hg19\" : \"CkoC44-NO-TIR1-pool.mcool\",\n",
    "    \"Ctrl500M-wtHAP1.hg19\" : \"Hap1-WT-combined-500000000.mcool\",\n",
    "    \"Ctrl500M-CT442-NT.hg19\" : \"CkoCT442-NT-pool-500000000.mcool\",\n",
    "    \"Ctrl500M-CT442-IAA.hg19\" : \"CkoCT442-IAA-pool-500000000.mcool\",\n",
    "    \n",
    "    ## A group for comparing PolII degron samples - i.e. with matched read depth to Hap1WT\n",
    "    # \"Hap1-WT-combined.hg19\" : \"Hap1-WT-combined.mcool\",\n",
    "    \"PolII-800M-IAA.hg19\" : \"PTB2539-IAA-pool-800000000.mcool\",\n",
    "    \"PolII-800M-NT.hg19\" : \"PTB2539-NT-pool-800000000.mcool\",\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%bash\n",
    "# cname=\"./CkoCT442-NT-pool.mcool\"\n",
    "# for uri in $(cooler ls $cname);\n",
    "# do\n",
    "#     echo $uri\n",
    "#     cooler balance --check $uri\n",
    "# done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "already exist !\n",
      "Hap1-WT-combined.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "PolII-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "PolII-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "all_CTCF_combined.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "noCTCF_IAA_combined.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CkoCT442_NT_pool.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CkoCT442_IAA_pool.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "RAD21-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "RAD21-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "WAPL-KO1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "WAPL-KO3.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "WAPL-WT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CkoC442-IAA48H-R1-T1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CkoC442-NT-R1-T1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CkoCT442-AAVS1sg2-4-NT-R1-T1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CkoCT442-DDX55sg2-B-NT-R1-T1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "G1-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "G1-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "mutControl-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "mutControl-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "mutDDX55-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "mutDDX55-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "mutTAF5L-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "mutTAF5L-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "siControl-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "siControl-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "siDDX55-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "siDDX55-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "siTAF5L-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "siTAF5L-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CtrlPlaB-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "CtrlPlaB-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "PlaB-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "PlaB-IAA.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "Ctrl500M-noTIR1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "Ctrl500M-wtHAP1.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "Ctrl500M-CT442-NT.hg19.5kb.100kb.bw\n",
      "already exist !\n",
      "Ctrl500M-CT442-IAA.hg19.5kb.100kb.bw\n",
      "working on ...\n",
      "PolII-800M-IAA.hg19 PTB2539-IAA-pool-800000000.mcool\n",
      "working on ...\n",
      "PolII-800M-NT.hg19 PTB2539-NT-pool-800000000.mcool\n"
     ]
    }
   ],
   "source": [
    "binsize = 5_000\n",
    "binsize_human = f\"{int(binsize/1_000)}kb\"\n",
    "diamond = 100_000\n",
    "diamond_human = f\"{int(diamond/1_000)}kb\"\n",
    "\n",
    "hg19_chromsizes = bioframe.fetch_chromsizes(\"hg19\")\n",
    "\n",
    "# iterate over samples to calculate insulation on:\n",
    "for k,cname in ins_samples.items():\n",
    "    target_bw_file = Path(f\"{k}.{binsize_human}.{diamond_human}.bw\")\n",
    "    if target_bw_file.is_file():\n",
    "        print(\"already exist !\")\n",
    "        print(target_bw_file)\n",
    "        continue\n",
    "    else:\n",
    "        print(\"working on ...\")\n",
    "        print(k,cname)\n",
    "        ######################################\n",
    "        ! cooltools insulation \\\n",
    "            -o {k}.{binsize_human}.{diamond_human}.tsv \\\n",
    "            {cname}::/resolutions/{binsize} \\\n",
    "            {diamond}\n",
    "        ###################################### \n",
    "        df = bioframe.read_table(\n",
    "                f\"{k}.{binsize_human}.{diamond_human}.tsv\",\n",
    "                header=0\n",
    "        )\n",
    "        ######################################\n",
    "        bioframe.to_bigwig(\n",
    "            df,\n",
    "            hg19_chromsizes,\n",
    "            f\"{k}.{binsize_human}.{diamond_human}.bw\",\n",
    "            value_field = f\"log2_insulation_score_{diamond}\"\n",
    "        )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # checking sequencing depth of several samples ...\n",
    "# !cooler info all_combined.mcool::/resolutions/25000 | grep sum\n",
    "# !cooler info all_combined.mcool::/resolutions/25000 | grep nnz\n",
    "# !cooler info noCTCF_IAA_combined.mcool::/resolutions/25000 | grep sum\n",
    "# !cooler info noCTCF_IAA_combined.mcool::/resolutions/25000 | grep nnz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating compartment status per arm for several samples ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rename_arms(name):\n",
    "    \"\"\"\n",
    "    turn UCSC genomic region name\n",
    "    into a human readable name of chrom arm\n",
    "    \"\"\"\n",
    "    chrom,coords = name.split(\":\")\n",
    "    start,end = coords.split(\"-\")\n",
    "    if start == \"0\":\n",
    "        return chrom+\"p\"\n",
    "    else:\n",
    "        return chrom+\"q\"\n",
    "\n",
    "hg19_chromsizes = bioframe.fetch_chromsizes('hg19', as_bed=True)\n",
    "hg19_cens = bioframe.fetch_centromeres('hg19')\n",
    "hg19_arms = bioframe.split(hg19_chromsizes, hg19_cens, cols_points=['chrom', 'mid'])\n",
    "# # Select only chromosomes that are present in the cooler. \n",
    "# # This step is typically not required! we call it only because the test data are reduced. \n",
    "# hg19_chromsizes = hg19_chromsizes.set_index(\"chrom\").loc[hg19_chromsizes.chrom].reset_index() \n",
    "hg19_arms = hg19_arms.set_index(\"chrom\").loc[hg19_chromsizes.chrom.values[:-1]].reset_index()\n",
    "# call this to automaticly assign names to chromosomal arms:\n",
    "hg19_arms = bioframe.parse_regions(hg19_arms)\n",
    "hg19_arms[\"name\"] = hg19_arms[\"name\"].apply(rename_arms)\n",
    "hg19_arms.head()\n",
    "hg19_arms.to_csv(\"hg19_arms.bed\",sep=\"\\t\",header=False,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "already exist !\n",
      "frac_gene_coverage.25kb.bedGraph\n"
     ]
    }
   ],
   "source": [
    "# get some cooler to extract bin table - for reference track extraction ...\n",
    "binsize = 25_000\n",
    "human_binsize = f\"{int(binsize/1_000)}kb\"\n",
    "db = \"hg19\"\n",
    "# input samples here\n",
    "# make sure we do not regenerate the same reference track again ...\n",
    "target_ref_track = Path(f\"frac_gene_coverage.{human_binsize}.bedGraph\")\n",
    "if target_ref_track.is_file():\n",
    "    print(\"already exist !\")\n",
    "    print(target_ref_track)\n",
    "else:\n",
    "    print(f\"working on ... {human_binsize} gene coverage for {db} ...\")\n",
    "    ######################################\n",
    "    # extract gene coverage for a given bin-size ...\n",
    "    # we need SOME cooler to get a bin table with that binsize ...\n",
    "    s= 'CkoCT442_NT_pool.hg19'\n",
    "    fname = ins_samples[s]\n",
    "    clr = cooler.Cooler(f\"{fname}::/resolutions/{binsize}\")\n",
    "    # gene cov extraction !\n",
    "    gene_cov_df = bioframe.frac_gene_coverage(\n",
    "                clr.bins()[:],\n",
    "                db,\n",
    "                )\n",
    "    # save coverage track ...\n",
    "    gene_cov_df[[\"chrom\",\"start\",\"end\",\"coverage\"]] \\\n",
    "        .to_csv(\n",
    "            target_ref_track,\n",
    "            sep=\"\\t\",\n",
    "            header=False,\n",
    "            index=False\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": ""
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "already exist !\n",
      "Hap1-WT-combined.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "PolII-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "PolII-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "all_CTCF_combined.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "noCTCF_IAA_combined.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "CkoCT442_NT_pool.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "CkoCT442_IAA_pool.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "RAD21-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "RAD21-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "CkoC442-IAA48H-R1-T1.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "CkoC442-NT-R1-T1.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "CkoCT442-AAVS1sg2-4-NT-R1-T1.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "CkoCT442-DDX55sg2-B-NT-R1-T1.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "G1-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "G1-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "mutControl-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "mutControl-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "mutDDX55-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "mutDDX55-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "mutTAF5L-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "mutTAF5L-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "siControl-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "siControl-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "siDDX55-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "siDDX55-IAA.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "siTAF5L-NT.hg19.EV.25kb.cis.bw\n",
      "already exist !\n",
      "siTAF5L-IAA.hg19.EV.25kb.cis.bw\n",
      "working on ...\n",
      "CtrlPlaB-NT.hg19 NT-hg19-combined-90000000.mcool\n",
      "to generate CtrlPlaB-NT.hg19.EV.25kb.cis.bw ...\n",
      "Tokenization took: 43.51 ms\n",
      "Type conversion took: 26.82 ms\n",
      "Parser memory cleanup took: 0.03 ms\n",
      "INFO:numexpr.utils:NumExpr defaulting to 8 threads.\n",
      "bedGraphToBigWig /tmp/tmpbay6v7ky.bg /tmp/tmpyi4ji_my.chrom.sizes CtrlPlaB-NT.hg19.EV.25kb.cis.bw\n",
      "working on ...\n",
      "CtrlPlaB-IAA.hg19 IAA-hg19-combined-90000000.mcool\n",
      "to generate CtrlPlaB-IAA.hg19.EV.25kb.cis.bw ...\n",
      "Tokenization took: 42.28 ms\n",
      "Type conversion took: 23.66 ms\n",
      "Parser memory cleanup took: 0.02 ms\n",
      "INFO:numexpr.utils:NumExpr defaulting to 8 threads.\n",
      "bedGraphToBigWig /tmp/tmpipe5et7e.bg /tmp/tmp16bdx9zg.chrom.sizes CtrlPlaB-IAA.hg19.EV.25kb.cis.bw\n",
      "working on ...\n",
      "PlaB-NT.hg19 NT-PlaB-hg19-combined-90000000.mcool\n",
      "to generate PlaB-NT.hg19.EV.25kb.cis.bw ...\n",
      "Tokenization took: 30.12 ms\n",
      "Type conversion took: 19.99 ms\n",
      "Parser memory cleanup took: 0.02 ms\n",
      "INFO:numexpr.utils:NumExpr defaulting to 8 threads.\n",
      "bedGraphToBigWig /tmp/tmpevejyov8.bg /tmp/tmp0tetqtut.chrom.sizes PlaB-NT.hg19.EV.25kb.cis.bw\n",
      "working on ...\n",
      "PlaB-IAA.hg19 IAA-PlaB-hg19-combined-90000000.mcool\n",
      "to generate PlaB-IAA.hg19.EV.25kb.cis.bw ...\n",
      "Tokenization took: 34.36 ms\n",
      "Type conversion took: 22.80 ms\n",
      "Parser memory cleanup took: 0.01 ms\n",
      "INFO:numexpr.utils:NumExpr defaulting to 8 threads.\n",
      "bedGraphToBigWig /tmp/tmpu60h05va.bg /tmp/tmpq34b_m1x.chrom.sizes PlaB-IAA.hg19.EV.25kb.cis.bw\n"
     ]
    }
   ],
   "source": [
    "binsize = 25_000\n",
    "human_binsize = f\"{int(binsize/1_000)}kb\"\n",
    "# iterate over samples to calculate insulation on:\n",
    "for k,cname in ins_samples.items():\n",
    "    target_bw_file = Path(f\"{k}.EV.{human_binsize}.cis.bw\")\n",
    "    if target_bw_file.is_file():\n",
    "        print(\"already exist !\")\n",
    "        print(target_bw_file)\n",
    "        continue\n",
    "    else:\n",
    "        print(\"working on ...\")\n",
    "        print(k, cname)\n",
    "        print(f\"to generate {target_bw_file} ...\")\n",
    "        ######################################\n",
    "        ! cooltools call-compartments \\\n",
    "          --reference-track frac_gene_coverage.{human_binsize}.bedGraph \\\n",
    "          --regions hg19_arms.bed \\\n",
    "          --contact-type cis \\\n",
    "          --n-eigs 3 \\\n",
    "          -v \\\n",
    "          -o {k}.EV.{human_binsize} \\\n",
    "          --bigwig \\\n",
    "        {cname}::/resolutions/{binsize}\n",
    "\n",
    "# output of call-compartments ...\n",
    "# out_prefix + \".\" + contact_type + \".bw\"\n",
    "# out_prefix + \".\" + contact_type + \".lam.txt\"\n",
    "# out_prefix + \".\" + contact_type + \".vecs.tsv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# binsize = 10_000\n",
    "# human_binsize = f\"{int(binsize/1_000)}kb\"\n",
    "# # input samples here\n",
    "# ev_samples = ['CkoCT442_NT_pool.hg19','CkoCT442_IAA_pool.hg19']\n",
    "# clrs = []\n",
    "# for s in ev_samples:\n",
    "#     fname = ins_samples[s]\n",
    "#     clr_uri = f\"{fname}::/resolutions/{binsize}\"\n",
    "#     clrs.append(cooler.Cooler(clr_uri))\n",
    "\n",
    "# # let's calculate gene coverage by pulling mRNA data from USCS\n",
    "# # for our hg19 genome ...\n",
    "# # this would need to be repeated for every genome, every binsize, etc\n",
    "# # probably better to download mRNA data locally to avoid redownloading ...\n",
    "# gene_cov_df = bioframe.frac_gene_coverage(\n",
    "#                 clrs[0].bins()[:],\n",
    "#                 \"hg19\"\n",
    "#                 )\n",
    "# hg19_chromsizes = bioframe.fetch_chromsizes(\"hg19\")\n",
    "\n",
    "# #saving to bigwig - whatever ...\n",
    "# for s,_r in zip(ev_samples, ev_res):\n",
    "#     # upack res, as store it:\n",
    "#     eigvals, eigvecs = _r\n",
    "#     bw_fname = f\"{s}.EV.{human_binsize}.bw\"\n",
    "#     # save compartment calls as bigwig\n",
    "#     bioframe.to_bigwig(\n",
    "#         eigvecs,\n",
    "#         hg19_chromsizes,\n",
    "#         bw_fname,\n",
    "#         value_field=\"E1\"\n",
    "#     )\n",
    "\n",
    "# # running eigdecomp ...\n",
    "# nproc = 1\n",
    "# if nproc > 1:\n",
    "#     pool = mp.Pool(nproc)\n",
    "#     map_ = pool.map\n",
    "# else:\n",
    "#     map_ = map\n",
    "\n",
    "\n",
    "# ev_res = []\n",
    "# for _clr in clrs:\n",
    "#     _gene_cov_df = copy(gene_cov_df)\n",
    "#     _res = eigdecomp.cooler_cis_eig(\n",
    "#         clr=_clr,\n",
    "#         bins=_gene_cov_df,\n",
    "#         regions=hg19_arms,\n",
    "#         n_eigs=3,\n",
    "#         phasing_track_col='coverage',\n",
    "#         balance='weight',\n",
    "#         ignore_diags=2,\n",
    "#         clip_percentile=99.9,\n",
    "#         sort_metric=None,\n",
    "#         map = map_\n",
    "#     )\n",
    "#     ev_res.append(_res)\n",
    "    \n",
    "# # plotting some stuff ...\n",
    "# from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "# chr1q_region = tuple(hg19_arms.iloc[1])[:-1]\n",
    "\n",
    "# for i,_r in enumerate(ev_res):\n",
    "#     a,b = _r\n",
    "#     ev1 = bioframe.select(b, chr1q_region)[\"E1\"].values\n",
    "\n",
    "#     fig = plt.figure(figsize=(15,6),constrained_layout=True)\n",
    "#     gs = fig.add_gridspec(1,3,width_ratios=[1,0.05,1])\n",
    "#     axl = fig.add_subplot(gs[0,0])\n",
    "#     cax = fig.add_subplot(gs[0,1])\n",
    "#     axr = fig.add_subplot(gs[0,2])\n",
    "\n",
    "#     hm = clrs[i].matrix(balance=True).fetch(chr1q_region)\n",
    "#     valid_mask = np.nan_to_num(hm).sum(axis=0) == 0.\n",
    "#     hm = hm[~valid_mask][:,~valid_mask]\n",
    "\n",
    "#     ccc = axl.imshow(hm,cmap=\"RdYlBu_r\",interpolation=\"nearest\",norm=LogNorm())\n",
    "#     plt.colorbar(ccc,cax=cax)\n",
    "\n",
    "#     divider = make_axes_locatable(axl)\n",
    "#     axev = divider.append_axes(\"bottom\", size=\"15%\", pad=0.08)\n",
    "\n",
    "#     axev.plot(ev1[~valid_mask],'b-')\n",
    "#     axev.axhline(0,color=\"grey\")\n",
    "#     axev.set_xlim(0,np.sum(~valid_mask))\n",
    "\n",
    "#     # color histogram\n",
    "#     _1, bin_edges, _2 = axr.hist(np.nan_to_num(np.log(hm[hm>0])).flatten(),bins=100,log=True)\n",
    "#     axr.bar(bin_edges[0],np.sum(hm==0),width=-.3,align=\"edge\",color=\"red\",alpha=0.4,label=\"zero-pixels\")\n",
    "#     axr.set_xlabel(\"heatmap_signal, excluding 0-pixels\")\n",
    "#     axr.legend(frameon=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}