2023 - figures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures\n",
    "\n",
    "Peaks / lengths / differences / jaccard etc."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from IPython.display import display\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import re\n",
    "from itertools import product\n",
    "\n",
    "sns.set_theme(style=\"whitegrid\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def bedl(file):\n",
    "    try:\n",
    "        tf = pd.read_csv(file, sep='\\t', header=None)\n",
    "        return tf[2] - tf[1]\n",
    "    except:\n",
    "        return np.zeros(0)  # Empty file\n",
    "\n",
    "\n",
    "def lines(file):\n",
    "    try:\n",
    "        tf = pd.read_csv(file, sep='\\t', header=None)\n",
    "        return len(tf)\n",
    "    except:\n",
    "        return 0  # Empty file\n",
    "\n",
    "\n",
    "def d(a, b):\n",
    "    return a / b if b != 0 else 0\n",
    "\n",
    "\n",
    "def last_col(file):\n",
    "    try:\n",
    "        cols = len(pd.read_csv(file, sep='\\t', nrows=1, header=None).columns)\n",
    "        return pd.read_csv(file, sep='\\t', header=None, usecols=[cols - 1])[cols - 1]\n",
    "    except:\n",
    "        return np.zeros(0)  # Empty file\n",
    "\n",
    "\n",
    "def sorted_file(file):\n",
    "    ts = tempfile.mktemp()\n",
    "    ! cat {file} | sort -k1,1 -k2,2n -k3,3n > {ts}\n",
    "    return ts"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": "# Load peaks\n",
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "GSE26320_PATH = os.path.expanduser('~/data/2023_GSE26320')\n",
    "GSE26320_CELLS = ['GM12878', 'HMEC', 'HSMM', 'K562', 'NHEK', 'NHLF', 'H1', 'Huvec', 'HepG2']\n",
    "MODIFICATIONS = ['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']\n",
    "GSE26320_REPS = ['rep1', 'rep2']\n",
    "\n",
    "IMMUNE_PATH = os.path.expanduser('~/data/2023_Immune')\n",
    "IMMUNE_CELLS = ['CD4ABT', 'TCellBB', 'BCell', 'TCell', 'Monocyte', 'PBMC', 'NK', 'CD34', 'CD4', ]  ## Longest first\n",
    "IMMUNE_REPS = ['rep1', 'rep2', 'rep3', '']\n",
    "\n",
    "ABF_PATH = os.path.expanduser('~/data/2018_chipseq_y20o20')\n",
    "ABF_MODIFICATIONS = ['k4me3', 'k27ac', 'k4me1', 'k27me3', 'k36me3']\n",
    "ABF_REPLICATES = [f'OD{i}' for i in range(1, 30)] + [f'YD{i}' for i in range(1, 30)]\n",
    "\n",
    "! mkdir -p {GSE26320_PATH}/pics\n",
    "! mkdir -p {IMMUNE_PATH}/pics\n",
    "! mkdir -p {ABF_PATH}/pics"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "TOOLS = ['MACS2', 'MACS2 broad', 'SPAN', 'SICER']\n",
    "palette = plt.get_cmap('tab10')\n",
    "TOOLS_PALETTE = {t: palette(i) for i, t in enumerate(TOOLS)}"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def load_peaks(path, suffix, modifications, cells, replicates):\n",
    "    df = pd.DataFrame(columns=['file', 'modification', 'cell', 'replicate', 'peaks'],\n",
    "                      dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith(suffix):\n",
    "            continue\n",
    "        mod = next((m for m in modifications if f'{m}_' in f), None)\n",
    "        cell = '' if '' in cells else next((c for c in cells if f'{c}_' in f), None)\n",
    "        rep = next((r for r in replicates if f'_{r}_' in f or f'_{r}-' in f or f.startswith(f'{r}_')), None)\n",
    "        if rep is None and '' in replicates:\n",
    "            rep = ''\n",
    "        if mod and cell is not None and rep is not None:\n",
    "            peaks_path = os.path.join(path, f)\n",
    "            ps = lines(peaks_path)\n",
    "            df.loc[len(df)] = (peaks_path, mod, cell, rep, ps)\n",
    "    return df"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "def load_peaks_path(path, modifications, cells, replicates):\n",
    "    df_macs2 = load_peaks(os.path.join(path, 'macs2'), '.narrowPeak', modifications, cells, replicates)\n",
    "    df_macs2['tool'] = 'MACS2'\n",
    "    print('MACS2', len(df_macs2))\n",
    "\n",
    "    df_macs2broad = load_peaks(os.path.join(path, 'macs2'), '.broadPeak', modifications, cells, replicates)\n",
    "    df_macs2broad['tool'] = 'MACS2 broad'\n",
    "    print('MACS2 broad', len(df_macs2broad))\n",
    "\n",
    "    df_sicer = load_peaks(os.path.join(path, 'sicer'), 'summary-FDR0.01', modifications, cells, replicates)\n",
    "    df_sicer['tool'] = 'SICER'\n",
    "    print('SICER', len(df_sicer))\n",
    "\n",
    "    df_span = load_peaks(os.path.join(path, 'span'), '.peak', modifications, cells, replicates)\n",
    "    df_span['tool'] = 'SPAN'\n",
    "    print('SPAN', len(df_span))\n",
    "\n",
    "    return pd.concat([df_macs2, df_macs2broad, df_sicer, df_span]).reset_index(drop=True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_encode_peaks = load_peaks_path(GSE26320_PATH, MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS)\n",
    "df_encode_peaks['dataset'] = 'ENCODE'\n",
    "df_encode_peaks.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_immune_peaks = load_peaks_path(IMMUNE_PATH, MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS)\n",
    "df_immune_peaks['dataset'] = 'Roadmap'\n",
    "df_immune_peaks.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def update_abf_modifications(df):\n",
    "    df.loc[df['modification'] == 'k4me3', 'modification'] = 'H3K4me3'\n",
    "    df.loc[df['modification'] == 'k27ac', 'modification'] = 'H3K27ac'\n",
    "    df.loc[df['modification'] == 'k4me1', 'modification'] = 'H3K4me1'\n",
    "    df.loc[df['modification'] == 'k27me3', 'modification'] = 'H3K27me3'\n",
    "    df.loc[df['modification'] == 'k36me3', 'modification'] = 'H3K36me3'"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_abf_peaks_all = load_peaks_path(ABF_PATH, ABF_MODIFICATIONS, [''], ABF_REPLICATES)\n",
    "update_abf_modifications(df_abf_peaks_all)\n",
    "df_abf_peaks_all['cell'] = ['OD' if 'OD' in r else 'YD' for r in df_abf_peaks_all['replicate']]\n",
    "df_abf_peaks_all['dataset'] = 'ABF'\n",
    "# Drop outliers\n",
    "ABF_OUTLIERS = {'k27ac': ['YD1', 'YD6'],\n",
    "                'k27me3': ['OD9', 'YD1', 'YD2', 'YD3', 'YD4', 'YD9', 'YD10', 'YD11'],\n",
    "                'k36me3': ['OD3', 'OD6', 'OD12', 'OD18', 'OD20', 'YD1', 'YD3', 'YD4', 'YD5'],\n",
    "                'k4me3': ['OD6', 'OD7', 'OD14', 'YD2', 'YD3', 'YD10', 'YD14']}\n",
    "t_outliers = [any(f'{o}_h3{m}'.lower() in f.lower() or f'{o}_{m}'.lower() in f.lower()\n",
    "                  for m, ol in ABF_OUTLIERS.items() for o in ol)\n",
    "              for f in df_abf_peaks_all['file']]\n",
    "print('Total outliers', sum(t_outliers))\n",
    "df_abf_peaks = df_abf_peaks_all[[not to for to in t_outliers]]\n",
    "\n",
    "df_abf_peaks.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": "df_peaks_full = pd.concat([df_encode_peaks, df_immune_peaks, df_abf_peaks]).reset_index(drop=True)",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "sns.set_style('whitegrid')\n",
    "plt.figure(figsize=(16, 4))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.barplot(data=df_peaks_full[df_peaks_full['modification'] == m], x='tool', y='peaks', hue='dataset',\n",
    "                capsize=.2, errwidth=2, edgecolor=\"black\",\n",
    "                # palette=TOOLS_PALETTE,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=df_peaks_full[df_peaks_full['modification'] == m], x='tool', y='peaks', hue='dataset',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.5, legend=False,\n",
    "                  order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('peaks')\n",
    "    if i < 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    # ax.set_ylim(0, 1.9e7)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(GSE26320_PATH_HG38 + f'/pics/library_mapped_detailed_all.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Outliers"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# print('Outliers')\n",
    "# for m in MODIFICATIONS:\n",
    "#     print(m)\n",
    "#     display(df_peaks_full[(df_peaks_full['modification'] == m)].sort_values(by=['peaks'], ascending=False).head(10))\n",
    "#     display(df_peaks_full[(df_peaks_full['modification'] == m)].sort_values(by=['peaks'], ascending=True).head(10))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# sns.set_style('whitegrid')\n",
    "# import math\n",
    "# for m, ds in product(MODIFICATIONS, ['ENCODE', 'Roadmap']):\n",
    "#     print(m, ds)\n",
    "#     t = df_peaks_full[(df_peaks_full['dataset'] == ds) & (df_peaks_full['modification'] == m)]\n",
    "#     experiments = [(c, r) for _, (c, r) in t[['cell', 'replicate']].drop_duplicates().iterrows()]\n",
    "#     ncols = int(math.floor(math.sqrt(len(experiments)) + 4))\n",
    "#     nrows = int(math.ceil(len(experiments) / ncols))\n",
    "#     plt.figure(figsize=(3 * ncols, 3 * nrows))\n",
    "#     axs = [plt.subplot(nrows, ncols, i + 1) for i in range(len(experiments))]\n",
    "#     for i, (c, r) in enumerate(experiments):\n",
    "#         ax = axs[i]\n",
    "#         ax.title.set_text(f'{ds} {m} {c} {r}')\n",
    "#         g_results = sns.barplot(data=t[(t['cell'] == c) & (t['replicate']==r)],\n",
    "#                             x='tool', y='peaks',\n",
    "#                             ax=ax,\n",
    "#                             capsize=.2, errwidth=2,\n",
    "#                             order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER']\n",
    "#                             )\n",
    "#         ax.xaxis.set_tick_params(rotation=90)\n",
    "#         ax.set_ylabel('Number')\n",
    "#         ax.set_ylim(0, t['peaks'].max() * 1.1)\n",
    "#     plt.tight_layout()\n",
    "#     plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Remove outliers\n",
    "OUTLIERS = [\n",
    "    ('H3K4me3', 'NK', ''),\n",
    "    ('H3K27ac', 'CD34', 'rep2'),\n",
    "    ('H3K4me1', 'NK', ''),\n",
    "    ('H3K27me3', 'BCell', ''),\n",
    "    ('H3K27me3', 'TCell', ''),\n",
    "]\n",
    "df_peaks = df_peaks_full.loc[[(m, c, r) not in OUTLIERS\n",
    "                              for _, (m, c, r) in\n",
    "                              df_peaks_full[['modification', 'cell', 'replicate']].iterrows()]].copy()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "sns.set_style('whitegrid')\n",
    "plt.figure(figsize=(14, 3.5))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    dfm = df_peaks[df_peaks['modification'] == m]\n",
    "    sns.barplot(data=dfm, x='tool', y='peaks', hue='dataset',\n",
    "                capsize=.2, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                # palette=TOOLS_PALETTE,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=dfm, x='tool', y='peaks', hue='dataset',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.5,\n",
    "                  order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  legend=False,\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('peaks')\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    # ax.set_ylim(0, 1.9e7)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(GSE26320_PATH_HG38 + f'/pics/library_mapped_detailed_all.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_peaks[(df_peaks['modification'] == 'H3K27ac') &\n",
    "         (df_peaks['dataset'] == 'Roadmap') &\n",
    "         (df_peaks['tool'] == 'MACS2 broad')].sort_values(by=['peaks'], ascending=False)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(14, 3.5))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.barplot(data=df_peaks[df_peaks['modification'] == m], x='tool', y='peaks', hue='dataset',\n",
    "                capsize=.2, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                # palette=TOOLS_PALETTE,\n",
    "                order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=df_peaks[df_peaks['modification'] == m], x='tool', y='peaks', hue='dataset',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.5, legend=False,\n",
    "                  order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                  # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  ax=ax)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    ax.title.set_text(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('peaks')\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    # ax.set_ylim(0, 1.9e7)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(GSE26320_PATH_HG38 + f'/pics/library_mapped_detailed_all.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Lengths"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "print('Load lengths')\n",
    "ts = []\n",
    "for dataset, mod, cell, rep, file, tool in tqdm(zip(\n",
    "        df_peaks['dataset'], df_peaks['modification'], df_peaks['cell'], df_peaks['replicate'], df_peaks['file'],\n",
    "        df_peaks['tool']\n",
    ")):\n",
    "    lengths = bedl(file)\n",
    "    t = pd.DataFrame(dict(length=lengths))\n",
    "    t = t.sample(min(len(t), 10_000)).copy()\n",
    "    t['dataset'] = dataset\n",
    "    t['modification'] = mod\n",
    "    t['cell'] = cell\n",
    "    t['replicate'] = rep\n",
    "    t['tool'] = tool\n",
    "    t['footprint'] = sum(lengths)\n",
    "    ts.append(t)\n",
    "df_lens = pd.concat(ts).reset_index(drop=True)\n",
    "del ts, t\n",
    "df_lens.sample(10)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for ds in ['ENCODE', 'Roadmap']:\n",
    "    print(ds)\n",
    "    plt.figure(figsize=(18, 3))\n",
    "    axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "    for i, m in enumerate(MODIFICATIONS):\n",
    "        ax = axs[i]\n",
    "        g_results = sns.histplot(data=df_lens[(df_lens['dataset'] == ds) & (df_lens['modification'] == m)],\n",
    "                                 x='length', hue='tool', ax=ax,\n",
    "                                 stat='density', common_bins=False, common_norm=False,\n",
    "                                 bins=50, kde=True, log_scale=True, alpha=0.2,\n",
    "                                 palette=TOOLS_PALETTE,\n",
    "                                 hue_order=TOOLS)\n",
    "        g_results.set(xscale='log')\n",
    "        g_results.set_ylim(0, 3)\n",
    "        g_results.set_xlim(1e2, 7e4)\n",
    "        ax.title.set_text(m)\n",
    "        if i < len(axs) - 1:\n",
    "            ax.legend().set_visible(False)\n",
    "        else:\n",
    "            sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for m in ['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']:\n",
    "    print(m)\n",
    "    plt.figure(figsize=(12, 3))\n",
    "    axs = [plt.subplot(1, 3, i + 1) for i in range(3)]\n",
    "    for i, tool in enumerate(['MACS2', 'MACS2 broad', 'SICER']):\n",
    "        ax = axs[i]\n",
    "        g_results = sns.histplot(data=df_lens[(df_lens['modification'] == m) & (df_lens['tool'] == tool)],\n",
    "                                 x='length', hue='dataset', ax=ax,\n",
    "                                 stat='density', common_bins=False, common_norm=False,\n",
    "                                 bins=30, kde=True, log_scale=True, alpha=0.5,\n",
    "                                 # palette=TOOLS_PALETTE,\n",
    "                                 )\n",
    "        g_results.set(xscale='log')\n",
    "        g_results.set_ylim(0, 2.5)\n",
    "        g_results.set_xlim(200, 5e4)\n",
    "        ax.title.set_text(f'{tool}')\n",
    "        if i < len(axs) - 1:\n",
    "            ax.legend().set_visible(False)\n",
    "        else:\n",
    "            sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "# plt.savefig(f'{GSE26320_PATH_HG38}/pics/H3K27ac_wiggle_mix_lengths.pdf', bbox_inches='tight', dpi=300)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(14, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    g_results = sns.histplot(data=df_lens[(df_lens['modification'] == m) & (df_lens['tool'] == 'SPAN')],\n",
    "                             x='length', hue='dataset', ax=ax,\n",
    "                             stat='density', common_bins=False, common_norm=False,\n",
    "                             bins=30, kde=True, log_scale=True, alpha=0.5,\n",
    "                             # palette=TOOLS_PALETTE,\n",
    "                             )\n",
    "    g_results.set(xscale='log')\n",
    "    g_results.set_ylim(0, 1.5)\n",
    "    g_results.set_xlim(200, 7e4)\n",
    "    ax.title.set_text(m)\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# plt.savefig(f'{GSE26320_PATH_HG38}/pics/H3K27ac_wiggle_mix_lengths.pdf', bbox_inches='tight', dpi=300)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_theme(style=\"whitegrid\")\n",
    "plt.figure(figsize=(14, 3.5))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=df_lens[df_lens['modification'] == m], x='tool', y='length', hue='dataset',\n",
    "                showfliers=False,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    # ax.set(yscale='log')\n",
    "\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    if m not in ['H3K36me3']:\n",
    "        ax.set_ylim(bottom=10, top=1e5)\n",
    "    else:\n",
    "        ax.set_ylim(bottom=10, top=1e5)\n",
    "    ax.set(yscale='log')\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_theme(style=\"whitegrid\")\n",
    "plt.figure(figsize=(14, 3.5))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=df_lens[df_lens['modification'] == m], x='tool', y='length', hue='dataset',\n",
    "                showfliers=False,\n",
    "                order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF', ],\n",
    "                ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    # ax.set(yscale='log')\n",
    "\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    if m not in ['H3K36me3']:\n",
    "        ax.set_ylim(bottom=20, top=1e5)\n",
    "    else:\n",
    "        ax.set_ylim(bottom=20, top=1e5)\n",
    "    ax.set(yscale='log')\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# df_lens['cell_rep'] = df_lens['cell'] + ' ' + df_lens['replicate']\n",
    "#\n",
    "# for m in MODIFICATIONS:\n",
    "#     print(m)\n",
    "#     plt.figure(figsize=(24, 12))\n",
    "#     t = list(df_lens[df_lens['modification'] == m].groupby(['dataset', 'cell_rep']))\n",
    "#     plots = len(t)\n",
    "#     ncols = int(math.ceil(math.sqrt(plots)))\n",
    "#     nrows = int(math.ceil(plots / ncols))\n",
    "#     axs = [plt.subplot(nrows, ncols, i + 1) for i in range(plots)]\n",
    "#\n",
    "#     for i, ((ds, cellrep), dt) in enumerate(t):\n",
    "#         ax = axs[i]\n",
    "#         sns.histplot(data=dt,\n",
    "#                      x='length', hue='tool',\n",
    "#                      stat='density', common_bins=False, common_norm=False,\n",
    "#                      bins=30, kde=True, log_scale=True, alpha=0.5,\n",
    "#                      palette=TOOLS_PALETTE,\n",
    "#                      ax=ax)\n",
    "#\n",
    "#         ax.set_title(f'{ds} {cellrep}')\n",
    "#         ax.set(xscale='log')\n",
    "#         ax.set_ylim(0, 2.5)\n",
    "#         ax.set_xlim(200, 7e4)\n",
    "#         if (i + 1) == plots:\n",
    "#             sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "#         else:\n",
    "#             ax.legend().set_visible(False)\n",
    "#         if i % ncols != 0:\n",
    "#             ax.set_ylabel(None)\n",
    "#         if i < ncols:\n",
    "#             # ax.set_xticks([], minor=False)\n",
    "#             # ax.set_xticklabels([])\n",
    "#             ax.set_xlabel(None)\n",
    "#     plt.tight_layout()\n",
    "#     plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": "# Average length",
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "from collections import Counter\n",
    "\n",
    "ts = []\n",
    "for (ds, m, c, r, tool), dft in df_lens[\n",
    "    ['dataset', 'modification', 'cell', 'replicate', 'tool', 'length', 'footprint']\n",
    "].groupby(['dataset', 'modification', 'cell', 'replicate', 'tool']):\n",
    "    dft.sort_values(by=['length'], inplace=True)\n",
    "    ts.append(dft.iloc[int(len(dft) * 0.1): int(len(dft) * 0.9)])\n",
    "    # ts.append(dft)\n",
    "ts = pd.concat(ts).reset_index(drop=True)\n",
    "df_lens_avg = ts.groupby(['dataset', 'modification', 'cell', 'replicate', 'tool']).agg('mean').reset_index()\n",
    "df_lens_avg['ds_c_r'] = df_lens_avg['dataset'] + ' ' + df_lens_avg['cell'] + ' ' + df_lens_avg['replicate']\n",
    "df_lens_avg.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=df_lens_avg[df_lens_avg['modification'] == m], x='tool', y='length', hue='dataset',\n",
    "                showfliers=False,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=df_lens_avg[df_lens_avg['modification'] == m], x='tool', y='length', hue='dataset',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.8, legend=False,\n",
    "                  order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    if m not in ['H3K27me3', 'H3K36me3']:\n",
    "        ax.set_ylim(bottom=200, top=8e3)\n",
    "    else:\n",
    "        ax.set_ylim(bottom=30, top=6e4)\n",
    "    ax.set(yscale='log')\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_lens_avg[(df_lens_avg['modification'] == 'H3K27me3') &\n",
    "            (df_lens_avg['dataset'] == 'Roadmap') &\n",
    "            (df_lens_avg['tool'] == 'SPAN')].sort_values(by=['length'], ascending=True)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_theme(style=\"whitegrid\")\n",
    "plt.figure(figsize=(14, 3.5))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.barplot(data=df_lens_avg[df_lens_avg['modification'] == m], x='tool', y='length', hue='dataset',\n",
    "                # order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                capsize=.2, err_kws={'linewidth': 2}, edgecolor=\"black\",\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=df_lens_avg[df_lens_avg['modification'] == m], x='tool', y='length', hue='dataset',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.5, legend=False,\n",
    "                  order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                  # order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    # ax.set(yscale='log')\n",
    "\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    ax.set_ylim(bottom=10, top=1.2e4)\n",
    "    ax.set(yscale='log')\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=df_lens_avg[df_lens_avg['modification'] == m], x='tool', y='footprint', hue='dataset',\n",
    "                showfliers=False,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    sns.stripplot(data=df_lens_avg[df_lens_avg['modification'] == m], x='tool', y='footprint', hue='dataset',\n",
    "                  dodge=True, size=2, palette='dark:black', alpha=0.8, legend=False,\n",
    "                  order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                  hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                  ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    # ax.set(yscale='log')\n",
    "\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    # if m not in ['H3K27me3', 'H3K36me3']:\n",
    "    #     ax.set_ylim(1, 4e3)\n",
    "    # else:\n",
    "    #     ax.set_ylim(1, 1.7e4)\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Distance between peaks"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Load distances')\n",
    "ts = []\n",
    "for dataset, mod, cell, rep, file, tool in tqdm(list(zip(\n",
    "        df_peaks['dataset'], df_peaks['modification'], df_peaks['cell'], df_peaks['replicate'], df_peaks['file'],\n",
    "        df_peaks['tool']\n",
    "))):\n",
    "    try:\n",
    "        tf = pd.read_csv(file, sep='\\t', header=None).sort_values(by=[0, 1, 2])\n",
    "        tss = []\n",
    "        for chr, tfc in tf.groupby(0):\n",
    "            starts = tfc[1]\n",
    "            ends = tfc[2]\n",
    "            ends = np.roll(ends, 1)\n",
    "            ds = starts - ends\n",
    "            t = pd.DataFrame(dict(distance=ds))\n",
    "            # Ignore first and last\n",
    "            tss.append(t.iloc[1:min(1_000, len(t) - 1)])\n",
    "        t = pd.concat(tss).reset_index(drop=True)\n",
    "        t['dataset'] = dataset\n",
    "        t['modification'] = mod\n",
    "        t['cell'] = cell\n",
    "        t['replicate'] = rep\n",
    "        t['tool'] = tool\n",
    "        ts.append(t)\n",
    "    except:\n",
    "        pass\n",
    "df_distances = pd.concat(ts).reset_index(drop=True)\n",
    "del ts, tss, t\n",
    "df_distances.sample(10)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=df_distances[df_distances['modification'] == m],\n",
    "                x='tool', y='distance', hue='dataset',\n",
    "                showfliers=False,\n",
    "                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                # hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                hue_order=['ENCODE', 'Roadmap', 'ABF'],\n",
    "                ax=ax)\n",
    "    ax.title.set_text(m)\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    # if m not in ['H3K27me3', 'H3K36me3']:\n",
    "    #     ax.set_ylim(bottom=200, top=8e3)\n",
    "    # else:\n",
    "    #     ax.set_ylim(bottom=30, top=6e4)\n",
    "    ax.set(yscale='log')\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Average distance"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Keep % of peaks, since outliers may significantly shift results\n",
    "ts = []\n",
    "for (ds, m, c, r, tool), dft in df_distances[\n",
    "    ['dataset', 'modification', 'cell', 'replicate', 'tool', 'distance']\n",
    "].groupby(['dataset', 'modification', 'cell', 'replicate', 'tool']):\n",
    "    dft.sort_values(by=['distance'], inplace=True)\n",
    "    ts.append(dft.iloc[int(len(dft) * 0.2): int(len(dft) * 0.8)])\n",
    "    # ts.append(dft)\n",
    "ts = pd.concat(ts).reset_index(drop=True)\n",
    "df_distances_avg = ts.groupby(['dataset', 'modification', 'cell', 'replicate', 'tool']).agg('mean').reset_index()\n",
    "df_distances_avg['ds_c_r'] = df_distances_avg['dataset'] + ' ' + df_distances_avg['cell'] + ' ' + df_distances_avg[\n",
    "    'replicate']\n",
    "df_distances_avg.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.set_theme(style=\"white\")\n",
    "df_lens['mindex'] = 0\n",
    "for i, m in enumerate(reversed(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3'])):\n",
    "    df_distances_avg.loc[df_distances_avg['modification'] == m, 'mindex'] = i\n",
    "plt.figure(figsize=(6, 3))\n",
    "# t = df_lens[df_lens['dataset'].isin(['ENCODE', 'Roadmap', 'ABF'])]\n",
    "t = df_distances_avg[df_distances_avg['dataset'].isin(['ENCODE', 'Roadmap', ])]\n",
    "g_results = sns.lineplot(data=t.sort_values(by=['mindex', 'tool', 'distance']),\n",
    "                         x='distance', y='mindex', hue='tool', orient='y', estimator='median',\n",
    "                         alpha=0.8,\n",
    "                         hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                         palette=TOOLS_PALETTE,\n",
    "                         errorbar=None,\n",
    "                         # errorbar=('pi', 25),\n",
    "                         marker='o', markersize=5\n",
    "                         )\n",
    "g_results.axes.set_yticks(range(5), minor=False)\n",
    "g_results.axes.set_yticklabels(reversed(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']))\n",
    "g_results.axes.set_ylabel('Target')\n",
    "g_results.axes.set_xlabel('Distance')\n",
    "g_results.set(xscale='log')\n",
    "# g_results.set_xlim(200, 1.0e4)\n",
    "g_results.axes.title.set_text(f'Peaks distance')\n",
    "g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "# plt.savefig(f'{ENCODE_REPROCESSED_PATH}/pics/lengths_all.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()\n",
    "sns.set_theme(style=\"whitegrid\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Avg length vs distance"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "sns.set_theme(style=\"white\")\n",
    "t = pd.merge(df_lens_avg[['dataset', 'modification', 'cell', 'replicate', 'tool', 'length']],\n",
    "             df_distances_avg[['dataset', 'modification', 'cell', 'replicate', 'tool', 'distance']],\n",
    "             on=['dataset', 'modification', 'cell', 'replicate', 'tool'])\n",
    "\n",
    "plt.figure(figsize=(15, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    # sns.kdeplot(t[(t['modification'] == m)],\n",
    "    #             palette=TOOLS_PALETTE,\n",
    "    #             hue_order=TOOLS,\n",
    "    #             x='length', y='distance', hue='tool',\n",
    "    #             alpha=0.3,\n",
    "    #             thresh=0.1,\n",
    "    #             fill=True,\n",
    "    #             legend=False,\n",
    "    #             ax = ax)\n",
    "    sns.scatterplot(t[(t['modification'] == m)],\n",
    "                    palette=TOOLS_PALETTE,\n",
    "                    hue_order=TOOLS,\n",
    "                    style='dataset',\n",
    "                    x='distance', y='length', hue='tool', alpha=0.8,\n",
    "                    ax=ax)\n",
    "    ax.set_title(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    ax.set(xscale='log')\n",
    "    ax.set(yscale='log')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "sns.set_theme(style=\"whitegrid\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "sns.set_theme(style=\"white\")\n",
    "t = pd.merge(df_lens_avg[['dataset', 'modification', 'cell', 'replicate', 'tool', 'length']],\n",
    "             df_distances_avg[['dataset', 'modification', 'cell', 'replicate', 'tool', 'distance']],\n",
    "             on=['dataset', 'modification', 'cell', 'replicate', 'tool'])\n",
    "\n",
    "plt.figure(figsize=(15, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    # sns.kdeplot(t[(t['modification'] == m)],\n",
    "    #             palette=TOOLS_PALETTE,\n",
    "    #             hue_order=TOOLS,\n",
    "    #             x='length', y='distance', hue='tool',\n",
    "    #             alpha=0.3,\n",
    "    #             thresh=0.1,\n",
    "    #             fill=True,\n",
    "    #             legend=False,\n",
    "    #             ax = ax)\n",
    "    sns.scatterplot(t[(t['modification'] == m) & (t['tool'] != 'SPAN')],\n",
    "                    palette=TOOLS_PALETTE,\n",
    "                    hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                    style='dataset',\n",
    "                    x='distance', y='length',  hue='tool', alpha=0.8,\n",
    "                    ax=ax)\n",
    "    ax.set_title(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    ax.set(xscale='log')\n",
    "    ax.set(yscale='log')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "sns.set_theme(style=\"whitegrid\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "# ATAC-seq, CTCF info\n",
    "\n"
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def simple_load_peaks(path, suffix):\n",
    "    df = pd.DataFrame(columns=['file', 'peaks'], dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith(suffix):\n",
    "            continue\n",
    "        file = os.path.join(path, f)\n",
    "        ps = lines(file)\n",
    "        df.loc[len(df)] = (file, ps)\n",
    "    return df\n",
    "\n",
    "\n",
    "def simple_load_peaks_path(path):\n",
    "    df_macs2 = simple_load_peaks(os.path.join(path, 'macs2'), '.narrowPeak')\n",
    "    df_macs2['tool'] = 'MACS2'\n",
    "    print('MACS2', len(df_macs2))\n",
    "\n",
    "    df_macs2_broad = simple_load_peaks(os.path.join(path, 'macs2'), '.broadPeak')\n",
    "    df_macs2_broad['tool'] = 'MACS2 broad'\n",
    "    print('MACS2 broad', len(df_macs2_broad))\n",
    "\n",
    "    df_sicer = simple_load_peaks(os.path.join(path, 'sicer'), '.scoreisland')\n",
    "    df_sicer['tool'] = 'SICER'\n",
    "    print('SICER', len(df_sicer))\n",
    "\n",
    "    df_span = simple_load_peaks(os.path.join(path, 'span'), '.peak')\n",
    "    df_span['tool'] = 'SPAN'\n",
    "    print('SPAN', len(df_span))\n",
    "\n",
    "    return pd.concat([df_macs2, df_macs2_broad, df_sicer, df_span]).reset_index(drop=True)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ATAC_PATH = os.path.expanduser('~/data/2024_Immgen')\n",
    "dfa = simple_load_peaks_path(ATAC_PATH)\n",
    "dfa"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for tool, file in tqdm(zip(dfa['tool'], dfa['file'])):\n",
    "    lengths = bedl(file)\n",
    "    t = pd.DataFrame(dict(length=lengths))\n",
    "    t.sort_values(by=['length'], inplace=True)\n",
    "    t = t.iloc[: int(len(t) * 0.9)]  # ATAC-seq may include large domains, skip them\n",
    "    t = t.sample(min(len(t), 10_000)).copy()\n",
    "    t['tool'] = tool\n",
    "    ts.append(t)\n",
    "dfa_lens = pd.concat(ts).reset_index(drop=True)\n",
    "dfa_lens['modification'] = 'ATAC'\n",
    "del ts\n",
    "display(dfa_lens.sample(10))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "CTCF_PATH = os.path.expanduser('~/data/2024_TFs')\n",
    "df_ctcf = simple_load_peaks_path(CTCF_PATH)\n",
    "display(df_ctcf)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "for tool, file in tqdm(zip(df_ctcf['tool'], df_ctcf['file'])):\n",
    "    lengths = bedl(file)\n",
    "    t = pd.DataFrame(dict(length=lengths))\n",
    "    t = t.sample(min(len(t), 10_000)).copy()\n",
    "    t['tool'] = tool\n",
    "    ts.append(t)\n",
    "dfa_ctcf_lens = pd.concat(ts).reset_index(drop=True)\n",
    "dfa_ctcf_lens['modification'] = 'CTCF'\n",
    "del ts\n",
    "display(dfa_ctcf_lens.sample(3))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Keep % of peaks, since outliers may significantly shift results\n",
    "tracks_per_dataset = 20\n",
    "ds_counter = Counter()\n",
    "ts = []\n",
    "for (ds, m, tool), dft in df_lens_avg.groupby(\n",
    "        ['dataset', 'modification', 'tool']):\n",
    "    if ds_counter[(ds, m, tool)] >= tracks_per_dataset:\n",
    "        break\n",
    "    ds_counter[(ds, m, tool)] += 1\n",
    "    ts.append(dft)\n",
    "df_lens_balanced_avg = pd.concat(ts).reset_index()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "t = pd.concat([\n",
    "    df_lens_balanced_avg,\n",
    "    dfa_ctcf_lens.groupby(['modification', 'tool']).agg('mean').reset_index(),\n",
    "    dfa_lens.groupby(['modification', 'tool']).agg('mean').reset_index()\n",
    "]).reset_index(drop=True)\n",
    "\n",
    "sns.set_theme(style=\"white\")\n",
    "t['mindex'] = 0\n",
    "for i, m in enumerate(reversed(['CTCF', 'ATAC', 'H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3'])):\n",
    "    t.loc[t['modification'] == m, 'mindex'] = i\n",
    "plt.figure(figsize=(6, 3))\n",
    "g_results = sns.lineplot(data=t.sort_values(by=['mindex', 'tool', 'length']),\n",
    "                         x='length', y='mindex', hue='tool', orient='y', estimator='median',\n",
    "                         alpha=0.8,\n",
    "                         hue_order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                         # hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "                         palette=TOOLS_PALETTE,\n",
    "                         errorbar=None,\n",
    "                         marker='o', markersize=5\n",
    "                         )\n",
    "g_results.axes.set_yticks(range(7), minor=False)\n",
    "g_results.axes.set_yticklabels(reversed(['CTCF', 'ATAC', 'H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']))\n",
    "g_results.axes.set_ylabel('Target')\n",
    "g_results.axes.set_xlabel('Length')\n",
    "g_results.set(xscale='log')\n",
    "g_results.set_xlim(100, 1.0e4)\n",
    "g_results.axes.title.set_text(f'')\n",
    "g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "# plt.savefig(f'{ENCODE_REPROCESSED_PATH}/pics/lengths_all.pdf', bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# BAMs"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def load_bams(path, modifications, cells, replicates, control='Input'):\n",
    "    df_bams = pd.DataFrame(columns=['bam_file', 'modification', 'cell', 'replicate', 'reads'], dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if not f.endswith('.bam'):\n",
    "            continue\n",
    "        cell = next((c for c in cells if c in f), None)\n",
    "        rep = next((r for r in replicates if r in f), None)\n",
    "        mod = next((m for m in modifications if m in f), None)\n",
    "        if mod is None and control is not None and control in f:\n",
    "            mod = 'Control'\n",
    "        if mod and cell is not None and rep is not None:\n",
    "            file = os.path.join(path, f)\n",
    "            reads = ! samtools view -F 0x04 -c {file}\n",
    "            reads = int(reads[0])\n",
    "            df_bams.loc[len(df_bams)] = (file, mod, cell, rep, reads)\n",
    "    return df_bams"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_encode_bams = load_bams(GSE26320_PATH + '/bams', MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS, 'Input')\n",
    "df_encode_bams['dataset'] = 'ENCODE'\n",
    "df_immune_bams = load_bams(IMMUNE_PATH + '/bams', MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS, 'Control')\n",
    "df_immune_bams['dataset'] = 'Roadmap'\n",
    "df_abf_bams = load_bams(ABF_PATH + '/bams', ABF_MODIFICATIONS, [''], ABF_REPLICATES, 'input')\n",
    "update_abf_modifications(df_abf_bams)\n",
    "df_abf_bams['dataset'] = 'ABF'\n",
    "reads = ! samtools view -F 0x04 -c {ABF_PATH}/bams/input.bam\n",
    "reads = int(reads[0])\n",
    "df_abf_bams.loc[len(df_abf_bams)] = (f'{ABF_PATH}/bams/input.bam', 'Control', '', '', reads, 'ABF')\n",
    "\n",
    "df_bams = pd.concat([df_encode_bams, df_immune_bams, df_abf_bams]).reset_index(drop=True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "reads",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.barplot(data=df_bams, hue='modification', x='dataset', y='reads',\n",
    "            capsize=.05, errwidth=2, edgecolor=\"black\",\n",
    "            hue_order=MODIFICATIONS + ['Control'],\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_bams, hue='modification', x='dataset', y='reads',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.5, legend=False,\n",
    "              hue_order=MODIFICATIONS + ['Control'],\n",
    "              ax=ax)\n",
    "ax.set_title('Library size')\n",
    "# Put a legend to the right of the current axis\n",
    "ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_peaks_lens = pd.merge(df_peaks[['dataset', 'modification', 'cell', 'replicate', 'tool', 'peaks']],\n",
    "                         df_bams[['dataset', 'modification', 'cell', 'replicate', 'reads']],\n",
    "                         on=['dataset', 'modification', 'cell', 'replicate'], how='inner')\n",
    "df_peaks_lens = pd.merge(df_peaks_lens,\n",
    "                         df_lens_avg[['dataset', 'modification', 'cell', 'replicate', 'tool', 'length', 'footprint']],\n",
    "                         on=['dataset', 'modification', 'cell', 'replicate', 'tool'], how='inner')\n",
    "df_peaks_lens.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for x, y in [('peaks', 'length'), ('peaks', 'reads')]:\n",
    "    print(f'\\n\\n{x} vs {y}')\n",
    "    plt.figure(figsize=(15, 3))\n",
    "    axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "    for i, m in enumerate(['H3K4me3', 'H3K27ac', 'H3K4me1', 'H3K27me3', 'H3K36me3']):\n",
    "        ax = axs[i]\n",
    "        sns.kdeplot(\n",
    "            df_peaks_lens[df_peaks_lens['modification'] == m], x=x, y=y, hue='tool',\n",
    "            alpha=0.8,\n",
    "            style='dataset',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=TOOLS,\n",
    "            thresh=0.05,\n",
    "            fill=True,\n",
    "            ax=ax)\n",
    "        ax.title.set_text(m)\n",
    "        if i < len(axs) - 1:\n",
    "            ax.legend().set_visible(False)\n",
    "        else:\n",
    "            sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "        if x != 'reads':\n",
    "            ax.xaxis.set_tick_params(rotation=90)\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": "# BigWigs stats",
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "CHROM_SIZES = {\n",
    "    c: s for _, (c, s) in pd.read_csv(os.path.join(GSE26320_PATH, 'hg38.chrom.sizes'),\n",
    "                                      sep='\\t', names=['chr', 'size']).iterrows() if '_' not in c\n",
    "}\n",
    "\n",
    "\n",
    "def load_bws(path, modification, cells, replicates, control='Input'):\n",
    "    df_bws = pd.DataFrame(columns=['file', 'modification', 'cell', 'replicate'], dtype=object)\n",
    "    for f in tqdm(os.listdir(path)):\n",
    "        if '.bw' not in f:\n",
    "            continue\n",
    "        mod = next((m for m in modification if m in f), None)\n",
    "        if mod is None and control is not None and control in f:\n",
    "            mod = 'Control'\n",
    "        cell = next((c for c in cells if c in f), None)\n",
    "        rep = next((r for r in replicates if r in f), None)\n",
    "        if mod and cell and rep is not None:\n",
    "            df_bws.loc[len(df_bws)] = ((os.path.join(path, f)), mod, cell, rep)\n",
    "    return df_bws"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "encode_bws = load_bws(GSE26320_PATH + '/bw', MODIFICATIONS, GSE26320_CELLS, GSE26320_REPS, 'Input')\n",
    "encode_bws['dataset'] = 'ENCODE'\n",
    "display(encode_bws.sample(3))\n",
    "\n",
    "immune_bws = load_bws(IMMUNE_PATH + '/bw', MODIFICATIONS, IMMUNE_CELLS, IMMUNE_REPS, 'Control')\n",
    "immune_bws['dataset'] = 'Roadmap'\n",
    "display(immune_bws.sample(3))\n",
    "all_bws = pd.concat([encode_bws, immune_bws]).reset_index(drop=False)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "import pyBigWig\n",
    "\n",
    "total_coverages = {}\n",
    "ts = []\n",
    "for _, (ds, m, c, r) in tqdm(list(all_bws[['dataset', 'modification', 'cell', 'replicate']].iterrows())):\n",
    "    print(ds, m, c, r)\n",
    "    t = all_bws[(all_bws['modification'] == m) &\n",
    "                (all_bws['cell'] == c) &\n",
    "                (all_bws['replicate'] == r) &\n",
    "                (all_bws['dataset'] == ds)]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    bw_path = t['file'].values[0]\n",
    "    # print(bw_path)\n",
    "\n",
    "    try:\n",
    "        with pyBigWig.open(bw_path) as bw:\n",
    "            total_coverage = sum(\n",
    "                bw.stats(chr, exact=True, type='sum')[0] for chr in CHROM_SIZES.keys() if '_' not in chr)\n",
    "            print('Total coverage', total_coverage)\n",
    "            # Multiplier to align BAM coverage with BigWig estimation\n",
    "            total_coverage *= 1.1e-2\n",
    "            total_coverages[(ds, m, c, r)] = total_coverage\n",
    "    except Exception:\n",
    "        pass\n",
    "        # print('ERROR')\n",
    "\n",
    "total_coverages_df = pd.DataFrame([(ds, m, c, r, n) for (ds, m, c, r), n in total_coverages.items()],\n",
    "                                  columns=['dataset', 'modification', 'cell', 'replicate', 'reads'])\n",
    "total_coverages_df.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# Prepare BAM files\n",
    "WD = os.path.expanduser('~/data/2023_wiggle')\n",
    "\n",
    "BLACKLIST_PATH = WD + '/hg38-blacklist.v2.bed'\n",
    "blacklist_df = pd.read_csv(BLACKLIST_PATH, sep='\\t', header=None)\n",
    "blacklist_df.rename({0: 'chr', 1: 'start', 2: 'end'}, axis=1, inplace=True)\n",
    "blacklist_df"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# def collect_bw_stats(bw_path, control_bw_path, chrom_sizes, blacklist, region_len, top_regions, work_regions, resolution):\n",
    "#     print('Prepare chromosome indexes')\n",
    "#     chr_indx = {}\n",
    "#     indx = 0\n",
    "#     for chr, chr_size in chrom_sizes.items():\n",
    "#         if '_' in chr:\n",
    "#             continue\n",
    "#         bins = int(math.floor(chr_size / region_len))\n",
    "#         chr_indx[chr] = (indx, indx + bins)\n",
    "#         indx += bins\n",
    "#     chromosomes = list(chr_indx.keys())\n",
    "# \n",
    "#     print('Compute coverage in regions')\n",
    "#     region_coverages = np.zeros(indx)\n",
    "#     with pyBigWig.open(bw_path) as bw:\n",
    "#         for chr, (istart, iend) in chr_indx.items():\n",
    "#             region_coverages[istart: iend] = \\\n",
    "#                 np.array([x or 0 for x in bw.stats(chr, 1, chrom_sizes[chr], nBins=iend - istart, exact=True, type='sum')])\n",
    "#     \n",
    "#     if blacklist is not None:\n",
    "#         print('Processing blacklisted regions')\n",
    "#         blacklist_marked = 0\n",
    "#         for chr, start, end in zip(blacklist['chr'], blacklist['start'], blacklist['end']):\n",
    "#             if chr not in chr_indx:\n",
    "#                 continue\n",
    "#             offset = chr_indx[chr][0]\n",
    "#             for x in range(int(math.floor(start / region_len)), int(math.ceil(end / region_len)) + 1):\n",
    "#                 region_coverages[min(x - offset, len(region_coverages) - 1)] = 0\n",
    "#                 blacklist_marked += 1\n",
    "#         print(f'Marked {blacklist_marked} / {len(region_coverages)} bins as blacklist')\n",
    "# \n",
    "#     if top_regions > len(region_coverages):\n",
    "#         raise Exception(f'Too many top regions {top_regions} > {len(region_coverages)}')\n",
    "#     \n",
    "#     print('Split by chromosomes')\n",
    "#     chr_xs = [[] for _ in range(len(chr_indx))]  # Create new array for each chromosome index\n",
    "#     ichr = 0\n",
    "#     for x in np.sort(np.argpartition(region_coverages, -top_regions)[-top_regions:]):\n",
    "#         while chr_indx[chromosomes[ichr]][1] < x:\n",
    "#             ichr += 1\n",
    "#         chr_xs[ichr].append(x)\n",
    "# \n",
    "#     print('Collect final regions')\n",
    "#     peaks = []\n",
    "#     for ichr, xs in enumerate(chr_xs):\n",
    "#         chr = chromosomes[ichr]\n",
    "#         offset = chr_indx[chr][0]\n",
    "#         for i, x in enumerate(xs):\n",
    "#             peaks.append((chr, (x - offset) * region_len, (x + 1 - offset) * region_len))\n",
    "#     df = pd.DataFrame(peaks, columns=['chromosome', 'start', 'end'])\n",
    "#     if len(df) > work_regions:\n",
    "#         print(f'Pick {work_regions} / {top_regions} uniform regions for computation speedup')\n",
    "#         step = int(math.ceil(len(df) / float(work_regions)))\n",
    "#         df = df.iloc[range(0, len(df), step)].copy()\n",
    "#     \n",
    "#     print('Collecting bigwig stats')\n",
    "#     sums, means, stds, control_sums, control_means, control_stds = [], [], [], [], [], []\n",
    "#     with pyBigWig.open(bw_path) as bw:\n",
    "#         for _, (chr, start, end) in df.iterrows():\n",
    "#             stats = np.array(bw.stats(\n",
    "#                 chr, start, end, nBins=int(math.ceil((end - start) / resolution)), exact=True, type='sum'))\n",
    "#             sums.append(np.sum(stats))\n",
    "#             means.append(np.mean(stats))\n",
    "#             stds.append(np.std(stats))\n",
    "#     df['sum'] = sums\n",
    "#     df['mean'] = means\n",
    "#     df['std'] = stds\n",
    "#     if control_bw_path is not None:\n",
    "#         with pyBigWig.open(control_bw_path) as cbw:\n",
    "#                 for _, (chr, start, end) in df[['chromosome', 'start', 'end']].iterrows():\n",
    "#                     control_stats = np.array(cbw.stats(\n",
    "#                         chr, start, end, nBins=int(math.ceil((end - start) / resolution)), exact=True, type='sum'))\n",
    "#                     control_sums.append(np.sum(control_stats))\n",
    "#                     control_means.append(np.mean(control_stats))\n",
    "#                     control_stds.append(np.std(control_stats))\n",
    "#         df['control_sum'] = control_sums\n",
    "#         df['control_mean'] = control_means\n",
    "#         df['control_std'] = control_stds\n",
    "#     return df"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from queue import PriorityQueue\n",
    "\n",
    "\n",
    "def collect_bw_stats(bw_path, control_bw_path, chrom_sizes, blacklist, region_len, top_regions, work_regions,\n",
    "                     resolution):\n",
    "    regions = sum(\n",
    "        int(math.floor(chr_size / region_len)) if '_' not in chr else 0 for chr, chr_size in chrom_sizes.items())\n",
    "    if top_regions > regions:\n",
    "        raise Exception(f'Too many top regions {top_regions} > {len(regions)}')\n",
    "\n",
    "    print('Compute coverage in regions')\n",
    "    region_coverages = PriorityQueue(maxsize=top_regions)\n",
    "    blacklist_marked = 0\n",
    "    with pyBigWig.open(bw_path) as bw:\n",
    "        for chr, chr_size in chrom_sizes.items():\n",
    "            if '_' in chr:\n",
    "                continue\n",
    "            bins = int(math.floor(chr_size / region_len))\n",
    "            blacklisted_bins = np.zeros(bins)\n",
    "            if blacklist is not None:\n",
    "                for _, (start, end) in blacklist[blacklist['chr'] == chr][['start', 'end']].iterrows():\n",
    "                    for i in range(int(math.floor(start / region_len)), int(math.ceil(end / region_len)) + 1):\n",
    "                        blacklisted_bins[min(i, bins - 1)] = 1\n",
    "                        blacklist_marked += 1\n",
    "            for i, cov in enumerate(bw.stats(chr, 0, chr_size, nBins=bins, exact=True, type='sum')):\n",
    "                if cov is None:\n",
    "                    continue  # VERY strange\n",
    "                if blacklisted_bins[i] == 0:\n",
    "                    if region_coverages.full():\n",
    "                        region_coverages.get()  # Removes the element with lowest coverage \n",
    "                    region_coverages.put((cov, chr, i))\n",
    "\n",
    "    if blacklist is not None:\n",
    "        print(f'Marked {blacklist_marked} / {regions} blacklisted regions')\n",
    "\n",
    "    peaks = []\n",
    "    while not region_coverages.empty():\n",
    "        cov, chr, i = region_coverages.get()\n",
    "        peaks.append((chr, i * region_len, (i + 1) * region_len, cov))\n",
    "\n",
    "    df = pd.DataFrame(peaks, columns=['chromosome', 'start', 'end', 'cov'])\n",
    "    df.sort_values(by=['cov'], ascending=False, inplace=True)\n",
    "\n",
    "    if len(df) > work_regions:\n",
    "        print(f'Pick {work_regions} / {top_regions} uniform regions for computation speedup')\n",
    "        step = int(math.ceil(len(df) / float(work_regions)))\n",
    "        df = df.iloc[range(0, len(df), step)].copy()\n",
    "\n",
    "    print('Collecting bigwig stats')\n",
    "    sums, means, stds, control_sums, control_means, control_stds = [], [], [], [], [], []\n",
    "    with pyBigWig.open(bw_path) as bw:\n",
    "        for _, (chr, start, end) in df[['chromosome', 'start', 'end']].iterrows():\n",
    "            stats = np.array(bw.stats(\n",
    "                chr, start, end, nBins=int(math.ceil((end - start) / resolution)), exact=True, type='sum'))\n",
    "            sums.append(np.sum(stats))\n",
    "            means.append(np.mean(stats))\n",
    "            stds.append(np.std(stats))\n",
    "    df['sum'] = sums\n",
    "    df['mean'] = means\n",
    "    df['std'] = stds\n",
    "    if control_bw_path is not None:\n",
    "        with pyBigWig.open(control_bw_path) as cbw:\n",
    "            for _, (chr, start, end) in df[['chromosome', 'start', 'end']].iterrows():\n",
    "                control_stats = np.array(cbw.stats(\n",
    "                    chr, start, end, nBins=int(math.ceil((end - start) / resolution)), exact=True, type='sum'))\n",
    "                control_sums.append(np.sum(control_stats))\n",
    "                control_means.append(np.mean(control_stats))\n",
    "                control_stds.append(np.std(control_stats))\n",
    "        df['control_sum'] = control_sums\n",
    "        df['control_mean'] = control_means\n",
    "        df['control_std'] = control_stds\n",
    "    return df"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "import math\n",
    "\n",
    "REGION_LEN = 10_000\n",
    "TOP_REGIONS = 10_000\n",
    "WORK_REGIONS = 200\n",
    "RESOLUTION = 100\n",
    "\n",
    "# Compute on top chromosomes for speed\n",
    "chrom_sizes_part = {f'chr{i + 1}': CHROM_SIZES[f'chr{i + 1}'] for i in range(3)}\n",
    "\n",
    "ts = []\n",
    "for _, (bw_path, ds, m, c, r) in all_bws[['file', 'dataset', 'modification', 'cell', 'replicate']].iterrows():\n",
    "    # if not (ds == 'Roadmap' and c in ['TCell', 'BCell', 'Monocyte']):\n",
    "    #     continue\n",
    "    if m == 'Control':\n",
    "        continue\n",
    "    cr = r\n",
    "    ct = all_bws[(all_bws['modification'] == 'Control') &\n",
    "                 (all_bws['cell'] == c) &\n",
    "                 (all_bws['replicate'] == cr) &\n",
    "                 (all_bws['dataset'] == ds)]\n",
    "    # Control can be different for Roadmap\n",
    "    if len(ct) == 0:\n",
    "        cr = ''\n",
    "        ct = all_bws[(all_bws['modification'] == 'Control') &\n",
    "                     (all_bws['cell'] == c) &\n",
    "                     (all_bws['replicate'] == cr) &\n",
    "                     (all_bws['dataset'] == ds)]\n",
    "    if len(ct) == 0:\n",
    "        continue\n",
    "    print(ds, m, c, r, bw_path)\n",
    "    cbw_path = ct['file'].values[0]\n",
    "    print(bw_path, cbw_path)\n",
    "    df = collect_bw_stats(bw_path, cbw_path, chrom_sizes_part, blacklist_df, REGION_LEN, TOP_REGIONS, WORK_REGIONS,\n",
    "                          RESOLUTION)\n",
    "    df['dataset'] = ds\n",
    "    df['modification'] = m\n",
    "    df['cell'] = c\n",
    "    df['replicate'] = r\n",
    "    df['control_replicate'] = cr\n",
    "    df['rpkm'] = df['sum'] / ((REGION_LEN / 1e3) * ((total_coverages[(ds, m, c, r)]) / 1e6))\n",
    "    df['control_reads'] = total_coverages[(ds, 'Control', c, cr)]\n",
    "    df['control_rpkm'] = df['control_sum'] / ((REGION_LEN / 1e3) * ((total_coverages[(ds, 'Control', c, cr)]) / 1e6))\n",
    "    df['signal_to_control'] = [d(s, c) for s, c in zip(df['rpkm'], df['control_rpkm'])]\n",
    "    ts.append(df)\n",
    "df_bw_stats = pd.concat(ts).reset_index(drop=True)\n",
    "df_bw_stats.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_bw_stats['std_mean'] = [d(std, mean) for std, mean in zip(df_bw_stats['std'], df_bw_stats['mean'])]\n",
    "df_bw_stats['control_std_mean'] = \\\n",
    "    [d(std, mean) for std, mean in zip(df_bw_stats['control_std'], df_bw_stats['control_mean'])]"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": "# Signal-to-noise, signal-to-control",
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "print('Assign groups by coverage percentiles individually in each track')\n",
    "ts = []\n",
    "for (ds, m, c, r), group in df_bw_stats.groupby(['dataset', 'modification', 'cell', 'replicate']):\n",
    "    group = group.sort_values(by=['mean'], ascending=False)\n",
    "    group['percentile'] = [10 * int(10.0 * i / len(group)) for i in range(len(group))]\n",
    "    ts.append(group)\n",
    "df_bw_stats_perc = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "df_bw_stats_perc.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Compute Signal-to-noise ratio')\n",
    "t = df_bw_stats_perc[\n",
    "    ['dataset', 'modification', 'cell', 'replicate', 'percentile', 'mean']\n",
    "].groupby(['dataset', 'modification', 'cell', 'replicate', 'percentile']).agg('mean').reset_index()\n",
    "ts = []\n",
    "for (ds, m, c, r), dt in t.groupby(['dataset', 'modification', 'cell', 'replicate']):\n",
    "    snr = dt[dt['percentile'] <= 20]['mean'].mean() / dt[dt['percentile'] >= 80]['mean'].mean()\n",
    "    ts.append((ds, m, c, r, snr))\n",
    "df_snr = pd.DataFrame(ts, columns=['dataset', 'modification', 'cell', 'replicate', 'signal_to_noise'])\n",
    "df_snr.sample(3)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Pick only top covered regions, where most of the peaks are located')\n",
    "df_bw_stats_signal_aggr = df_bw_stats_perc[df_bw_stats_perc['percentile'] <= 20][\n",
    "    ['dataset', 'modification', 'cell', 'replicate', 'sum', 'mean', 'std', 'rpkm', 'std_mean', 'signal_to_control',\n",
    "     'control_reads', 'control_sum', 'control_mean', 'control_std', 'control_rpkm', 'control_std_mean', ]\n",
    "].groupby(['dataset', 'modification', 'cell', 'replicate']).agg('mean').reset_index()\n",
    "df_bw_stats_signal_aggr.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_bw_stats_signal_snr_aggr = pd.merge(df_snr, df_bw_stats_signal_aggr,\n",
    "                                       on=['dataset', 'modification', 'cell', 'replicate'], how='inner')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(3, 4))\n",
    "ax = plt.axes()\n",
    "t = df_bw_stats_signal_snr_aggr\n",
    "sns.boxplot(data=t, x='modification', y='std_mean', hue='dataset',\n",
    "            order=MODIFICATIONS,\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "ax.title.set_text('Std/mean')\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_bw_stats_signal_snr_aggr[(df_bw_stats_signal_snr_aggr['modification'] == 'H3K4me3') &\n",
    "                            (df_bw_stats_signal_snr_aggr['cell'] == 'CD4')].sort_values(by=['std_mean'],\n",
    "                                                                                        ascending=False).head(50)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Compare vs peaks"
  },
  {
   "cell_type": "code",
   "source": [
    "df_aggr = pd.merge(df_bw_stats_signal_snr_aggr,\n",
    "                   df_peaks_lens,\n",
    "                   on=['dataset', 'modification', 'cell', 'replicate'], how='inner')\n",
    "df_aggr['reads_min'] = [min(s, c) for s, c in zip(df_aggr['reads'], df_aggr['control_reads'])]\n",
    "df_aggr"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_aggr[['dataset', 'modification', 'cell', 'replicate']].drop_duplicates().groupby(['dataset', 'modification']).size()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for m, cells in [\n",
    "    ('H3K4me1', ['H1']),\n",
    "    ('H3K27me3', ['CD4', 'Monocyte']),\n",
    "    ('H3K36me3', ['CD4']),\n",
    "    ('H3K27me3', ['HepG2']),\n",
    "]:\n",
    "    print(m, cells)\n",
    "    display(df_aggr[(df_aggr['modification'] == m) &\n",
    "                    (df_aggr['cell'].isin(cells))][['cell', 'replicate', 'reads', 'rpkm', 'std_mean',\n",
    "                                                    'control_reads', 'control_rpkm', 'control_std_mean',\n",
    "                                                    'reads_min',\n",
    "                                                    'signal_to_noise', 'signal_to_control', ]].drop_duplicates())\n",
    "    display(df_peaks_lens[(df_peaks_lens['modification'] == m) &\n",
    "                          (df_peaks_lens['cell'].isin(cells))][['cell', 'replicate', 'tool', 'peaks', 'length']])"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "features = ['reads', 'mean', 'rpkm', 'std_mean',\n",
    "            'control_reads', 'control_mean', 'control_rpkm', 'control_std_mean',\n",
    "            'reads_min',\n",
    "            'signal_to_noise', 'signal_to_control', ]\n",
    "ts = []\n",
    "ts.append(ts.append(\n",
    "    df_aggr[df_aggr['tool'] == 'SPAN'][['dataset', 'modification', 'cell', 'replicate'] + features].copy().reset_index(\n",
    "        drop=True)))\n",
    "for tool in TOOLS:\n",
    "    t = df_aggr[df_aggr['tool'] == tool][['peaks', 'length']].copy().reset_index(drop=True)\n",
    "    t.rename({'peaks': f'{tool} peaks', 'length': f'{tool} length'}, axis=1, inplace=True)\n",
    "ts.append(t)\n",
    "tt = pd.concat(ts, axis=1).reset_index(drop=True)\n",
    "display(tt.head(10))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "features = ['reads', 'mean', 'rpkm', 'std_mean',\n",
    "            'control_reads', 'control_mean', 'control_rpkm', 'control_std_mean',\n",
    "            'reads_min',\n",
    "            'signal_to_noise', 'signal_to_control', ]\n",
    "ts = []\n",
    "ts.append(ts.append(\n",
    "    df_aggr[df_aggr['tool'] == 'SPAN'][['dataset', 'modification', 'cell', 'replicate'] + features].copy().reset_index(\n",
    "        drop=True)))\n",
    "for tool in TOOLS:\n",
    "    t = df_aggr[df_aggr['tool'] == tool][['peaks', 'length']].copy().reset_index(drop=True)\n",
    "    t.rename({'peaks': f'{tool} peaks', 'length': f'{tool} length'}, axis=1, inplace=True)\n",
    "    ts.append(t)\n",
    "tt = pd.concat(ts, axis=1).reset_index(drop=True)\n",
    "display(tt.head(10))\n",
    "features = list(tt.columns)[4:]\n",
    "del ts\n",
    "ncols = int(math.ceil(math.sqrt(len(features) + 3)))\n",
    "nrows = int(math.ceil(len(features) / ncols))\n",
    "plt.figure(figsize=(2.5 * ncols, 2 * nrows))\n",
    "axs = [plt.subplot(nrows, ncols, i + 1) for i in range(len(features))]\n",
    "for i, y in enumerate(features):\n",
    "    ax = axs[i]\n",
    "    sns.boxplot(data=tt, x='modification', y=y,\n",
    "                hue='dataset',\n",
    "                order=MODIFICATIONS,\n",
    "                showfliers=False,\n",
    "                ax=ax)\n",
    "    # sns.stripplot(data=tt, x='modification', y=y, hue='dataset',\n",
    "    #               order=MODIFICATIONS,\n",
    "    #               dodge=True, size=1.5, alpha=0.5, palette='dark:black', legend=False,\n",
    "    #               ax=ax)\n",
    "    ax.title.set_text(y)\n",
    "    ax.set_ylabel(None)\n",
    "    if i < (nrows - 1) * ncols:\n",
    "        ax.set_xlabel(None)\n",
    "        ax.set_xticks([], minor=False)\n",
    "        ax.set_xticklabels([])\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "    if i % ncols == ncols - 1:\n",
    "        sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    else:\n",
    "        ax.legend().set_visible(False)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "sns.pairplot(data=tt[['dataset'] + features], hue='dataset', diag_kind='hist')\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "corrs = tt[features].corr(method='spearman')\n",
    "sns.clustermap(corrs,\n",
    "               col_cluster=True, row_cluster=True,\n",
    "               figsize=(6, 6),\n",
    "               cmap=plt.cm.seismic,\n",
    "               vmin=-1, vmax=1)\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "t = np.zeros(shape=(len(features) - 8, 4 * 5 * 3))\n",
    "\n",
    "for i, (tool, m) in enumerate(product(TOOLS, MODIFICATIONS)):\n",
    "    ttm = tt[tt['modification'] == m]\n",
    "    t[:, i * 3] = ttm.corr(method='spearman').loc[features[:-8]][f'{tool} peaks'].values\n",
    "    t[:, i * 3 + 1] = ttm[ttm['dataset'] == 'ENCODE'][features].corr(method='spearman').loc[features[:-8]][\n",
    "        f'{tool} peaks'].values\n",
    "    t[:, i * 3 + 2] = ttm[ttm['dataset'] == 'Roadmap'][features].corr(method='spearman').loc[features[:-8]][\n",
    "        f'{tool} peaks'].values\n",
    "t = pd.DataFrame(t,\n",
    "                 columns=[f'{s} {m} {tool}' for tool, m, s in\n",
    "                          product(TOOLS, MODIFICATIONS, ['Full', 'ENCODE', 'Roadmap'])],\n",
    "                 index=features[:-8])\n",
    "sns.clustermap(t,\n",
    "               col_cluster=False, row_cluster=True,\n",
    "               figsize=(20, 6),\n",
    "               cmap=plt.cm.seismic,\n",
    "               vmin=-1, vmax=1)\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Predict number of peaks"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "for tool in TOOLS:\n",
    "    plt.figure(figsize=(15, 3.2))\n",
    "    axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "    for i, m in enumerate(MODIFICATIONS):\n",
    "        ax = axs[i]\n",
    "        tt = df_aggr[(df_aggr['modification'] == m) & (df_aggr['tool'] == tool)].copy()\n",
    "        tt.fillna(0, inplace=True)\n",
    "        features_predict = [f for f in features if f not in ['peaks', 'length', 'footprint']]\n",
    "        # features_predict = [f for f in features if 'control' in features or f == 'reads_min']\n",
    "        x_all = tt[features_predict]\n",
    "        y_all = tt['peaks']\n",
    "        X, X_val, y, y_val = train_test_split(x_all, y_all)\n",
    "        Rf = RandomForestRegressor(n_estimators=10, n_jobs=-1)\n",
    "        Rf.fit(X, y)\n",
    "        print(f'{tool} {m} R2: train {Rf.score(X, y):.2} validate: {Rf.score(X_val, y_val):.2}')\n",
    "        ax.title.set_text(f'{tool} {m} R2 {Rf.score(x_all, y_all):.2}')\n",
    "        # plot the model\n",
    "        y_pred = Rf.predict(x_all)\n",
    "        sns.scatterplot(x=y_all, y=y_pred, ax=ax)\n",
    "        ax.set_xlabel('real')\n",
    "        if i > 0:\n",
    "            ax.set_ylabel(None)\n",
    "        else:\n",
    "            ax.set_ylabel('predicted')\n",
    "        # generate and graph y = x line\n",
    "        maxxy = max(np.max(y_pred), np.max(y_all))\n",
    "        minxy = min(np.min(y_pred), np.min(y_all))\n",
    "        sns.lineplot(x=[minxy, maxxy], y=[minxy, maxxy], color='black', ax=ax, alpha=0.5, linestyle='dotted')\n",
    "        ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# ATAC-seq, TFs BigWigs"
  },
  {
   "cell_type": "code",
   "source": [
    "REGION_LEN = 10_000\n",
    "TOP_REGIONS = 10_000\n",
    "WORK_REGIONS = 200\n",
    "RESOLUTION = 100\n",
    "\n",
    "CHROM_SIZES_MM10 = {\n",
    "    c: s for _, (c, s) in pd.read_csv(os.path.join(ATAC_PATH, 'mm10.chrom.sizes'),\n",
    "                                      sep='\\t', names=['chr', 'size']).iterrows() if '_' not in c\n",
    "}\n"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "atac_bws = load_bws(ATAC_PATH + '/bw', ['ATAC'], ['Monocyte'], ['']).head(5)\n",
    "atac_bws['dataset'] = 'ATAC-seq'\n",
    "atac_bws['replicate'] = atac_bws.index  # Override replicate\n",
    "atac_bws"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Compute on top chromosomes for speed\n",
    "chrom_sizes_part_mm10 = {f'chr{i + 1}': CHROM_SIZES_MM10[f'chr{i + 1}'] for i in range(3)}\n",
    "\n",
    "ts = []\n",
    "for _, (bw_path, ds, m, c, r) in atac_bws[['file', 'dataset', 'modification', 'cell', 'replicate']].iterrows():\n",
    "    print(ds, m, c, r, bw_path)\n",
    "    # continue\n",
    "    df = collect_bw_stats(\n",
    "        bw_path, None, chrom_sizes_part_mm10, None,\n",
    "        REGION_LEN, TOP_REGIONS, WORK_REGIONS, RESOLUTION\n",
    "    )\n",
    "    df['dataset'] = ds\n",
    "    df['modification'] = m\n",
    "    df['cell'] = c\n",
    "    df['replicate'] = r\n",
    "    ts.append(df)\n",
    "df_atac_bw_stats = pd.concat(ts).reset_index(drop=True)\n",
    "df_atac_bw_stats.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "CTCF_PATH = os.path.expanduser('~/data/2024_TFs')\n",
    "ctcf_bws = load_bws(CTCF_PATH + '/bw', ['CTCF'], ['H1'], ['']).head(5)\n",
    "ctcf_bws['dataset'] = 'CTCF'\n",
    "ctcf_bws['replicate'] = ctcf_bws.index  # Override replicate\n",
    "ctcf_bws"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Compute on top chromosomes for speed\n",
    "chrom_sizes_part = {f'chr{i + 1}': CHROM_SIZES[f'chr{i + 1}'] for i in range(3)}\n",
    "\n",
    "ts = []\n",
    "for _, (bw_path, ds, m, c, r) in ctcf_bws[['file', 'dataset', 'modification', 'cell', 'replicate']].iterrows():\n",
    "    print(ds, m, c, r, bw_path)\n",
    "    # continue\n",
    "    df = collect_bw_stats(\n",
    "        bw_path, None, chrom_sizes_part, None,\n",
    "        REGION_LEN, TOP_REGIONS, WORK_REGIONS, RESOLUTION\n",
    "    )\n",
    "    df['dataset'] = ds\n",
    "    df['modification'] = m\n",
    "    df['cell'] = c\n",
    "    df['replicate'] = r\n",
    "    ts.append(df)\n",
    "df_ctcf_bw_stats = pd.concat(ts).reset_index(drop=True)\n",
    "df_ctcf_bw_stats.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "df_atac_tf_bw_stats = pd.concat([df_atac_bw_stats, df_ctcf_bw_stats]).reset_index(drop=True)",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_atac_tf_bw_stats['std_mean'] = [\n",
    "    d(std, mean) for std, mean in zip(df_atac_tf_bw_stats['std'], df_atac_tf_bw_stats['mean'])\n",
    "]"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Assign groups by coverage percentiles individually in each track')\n",
    "ts = []\n",
    "for (ds, m, c, r), group in df_atac_tf_bw_stats.groupby(['dataset', 'modification', 'cell', 'replicate']):\n",
    "    group = group.sort_values(by=['mean'], ascending=False)\n",
    "    group['percentile'] = [10 * int(10.0 * i / len(group)) for i in range(len(group))]\n",
    "    ts.append(group)\n",
    "df_atac_tf_bw_stats_perc = pd.concat(ts).reset_index(drop=True)\n",
    "del ts\n",
    "df_atac_tf_bw_stats_perc.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Compute Signal-to-noise ratio')\n",
    "t = df_atac_tf_bw_stats_perc[\n",
    "    ['dataset', 'modification', 'cell', 'replicate', 'percentile', 'mean']\n",
    "].groupby(['dataset', 'modification', 'cell', 'replicate', 'percentile']).agg('mean').reset_index()\n",
    "ts = []\n",
    "for (ds, m, c, r), dt in t.groupby(['dataset', 'modification', 'cell', 'replicate']):\n",
    "    snr = dt[dt['percentile'] <= 20]['mean'].mean() / dt[dt['percentile'] >= 80]['mean'].mean()\n",
    "    ts.append((ds, m, c, r, snr))\n",
    "df_atac_tf_snr = pd.DataFrame(ts, columns=['dataset', 'modification', 'cell', 'replicate', 'signal_to_noise'])\n",
    "df_atac_tf_snr.sample(3)\n"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('Pick only top covered regions, where most of the peaks are located')\n",
    "df_atac_tf_bw_stats_signal_aggr = df_atac_tf_bw_stats_perc[df_atac_tf_bw_stats_perc['percentile'] <= 20][\n",
    "    ['dataset', 'modification', 'cell', 'replicate', 'sum', 'mean', 'std', 'std_mean', ]\n",
    "].groupby(['dataset', 'modification', 'cell', 'replicate']).agg('mean').reset_index()\n",
    "df_atac_tf_bw_stats_signal_aggr.sample(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_atac_tf_bw_stats_signal_snr_aggr = pd.merge(df_atac_tf_snr, df_atac_tf_bw_stats_signal_aggr,\n",
    "                                               on=['dataset', 'modification', 'cell', 'replicate'], how='inner')\n",
    "df_atac_tf_bw_stats_signal_snr_aggr.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(2, 3))\n",
    "ax = plt.axes()\n",
    "t = pd.concat([df_atac_tf_bw_stats_signal_snr_aggr, df_bw_stats_signal_snr_aggr]).reset_index(drop=True)\n",
    "# t = df_atac_tf_bw_stats_signal_snr_aggr\n",
    "sns.boxplot(data=t, x='modification', y='std_mean',\n",
    "            order=['CTCF', 'ATAC'] + MODIFICATIONS,\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "ax.title.set_text('Std/mean')\n",
    "ax.xaxis.set_tick_params(rotation=90)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "STOP",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Difference"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_diffs(diff_path, df):\n",
    "    tools = df['tool'].unique()\n",
    "    rows = []\n",
    "    for (m, c, r), tmcr in tqdm(list(df.groupby(['modification', 'cell', 'replicate']))):\n",
    "        print(m, c, r)\n",
    "        t = tmcr[tmcr['tool'] == 'SPAN']\n",
    "        if len(t) == 0:\n",
    "            continue\n",
    "        span_file = sorted_file(t['file'].values[0])\n",
    "        # Processing single tools information\n",
    "        for tool in tools:\n",
    "            if tool == 'SPAN':\n",
    "                continue\n",
    "            t = tmcr[tmcr['tool'] == tool]\n",
    "            if len(t) == 0:\n",
    "                continue\n",
    "            peaks_file = sorted_file(t['file'].values[0])\n",
    "            for name, args in [\n",
    "                (f'SPAN - {tool}', f' -a {span_file} -b {peaks_file} '),\n",
    "                (f'{tool} - SPAN', f' -b {span_file} -a {peaks_file} ')]:\n",
    "                diff_file = f'{diff_path}/{name.replace(\" \", \"_\")}_{m}_{c}_{r}.bed'\n",
    "                ! bedtools intersect {args} -wa -v > {diff_file}\n",
    "                peaks = lines(diff_file)\n",
    "                rows.append((m, c, r, name, peaks))\n",
    "    return pd.DataFrame(rows, columns=['modification', 'cell', 'replicate', 'name', 'peaks'])"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def plot_diff(df_peaks, df_diff):\n",
    "    t = df_peaks[['modification', 'cell', 'replicate', 'tool', 'peaks']].copy()\n",
    "    t.rename({'tool': 'name'}, axis=1, inplace=True)\n",
    "    t = pd.concat([t, df_diff]).reset_index(drop=True)\n",
    "    t = t.loc[\n",
    "        [(m, c, r) not in OUTLIERS for _, (m, c, r) in t[['modification', 'cell', 'replicate']].iterrows()]].copy()\n",
    "\n",
    "    plt.figure(figsize=(18, 4))\n",
    "    axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "\n",
    "    for i, m in enumerate(MODIFICATIONS):\n",
    "        ax = axs[i]\n",
    "        tm = t[t['modification'] == m]\n",
    "        print(m, len(t))\n",
    "        ax.set_title(m)\n",
    "        sns.barplot(data=tm, x='name', y='peaks',\n",
    "                    color='grey',\n",
    "                    capsize=.2, err_kws={'linewidth': 2},\n",
    "                    order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER',\n",
    "                           'SPAN - MACS2', 'SPAN - MACS2 broad', 'SPAN - SICER',\n",
    "                           'MACS2 - SPAN', 'MACS2 broad - SPAN', 'SICER - SPAN',\n",
    "                           ],\n",
    "                    ax=ax)\n",
    "        sns.stripplot(data=tm, x='name', y='peaks',\n",
    "                      size=2, alpha=0.8, color='black',\n",
    "                      order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER',\n",
    "                             'SPAN - MACS2', 'SPAN - MACS2 broad', 'SPAN - SICER',\n",
    "                             'MACS2 - SPAN', 'MACS2 broad - SPAN', 'SICER - SPAN',\n",
    "                             ],\n",
    "                      legend=False,\n",
    "                      ax=ax)\n",
    "        ax.xaxis.set_tick_params(rotation=90)\n",
    "        ax.set_ylabel('Peaks number')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## ENCODE"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "dp = f'{GSE26320_PATH}/diff'\n",
    "! mkdir -p {dp}\n",
    "df_encode_diff = compute_diffs(dp, df_encode_peaks)\n",
    "df_encode_diff.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": "plot_diff(df_encode_peaks, df_encode_diff)",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Roadmap"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "dp = f'{IMMUNE_PATH}/diff'\n",
    "! mkdir -p {dp}\n",
    "df_immune_diff = compute_diffs(dp, df_immune_peaks)\n",
    "df_immune_diff.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": "plot_diff(df_immune_peaks, df_immune_diff)",
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_immune_diff[(df_immune_diff['modification'] == 'H3K27ac') &\n",
    "               (df_immune_diff['name'] == 'MACS2 - SPAN')].sort_values(by=['peaks'], ascending=False).head(10)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## ABF"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "dp = f'{ABF_PATH}/diff'\n",
    "! mkdir -p {dp}\n",
    "\n",
    "df_abf_diff = compute_diffs(dp, df_abf_peaks)\n",
    "df_abf_diff.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "plot_diff(df_abf_peaks, df_abf_diff)\n",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Jaccard & Overlap"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "import pyranges as pr\n",
    "\n",
    "\n",
    "def compute_reps_overlap(df_peaks):\n",
    "    reps_overlap = pd.DataFrame(columns=['modification', 'cell', 'tool', 'rep1', 'rep2',\n",
    "                                         'peaks1', 'peaks1_len', 'peaks2', 'peaks2_len',\n",
    "                                         'peaks1_overlap', 'peaks2_overlap', 'peaks_overlap_len'], dtype=object)\n",
    "\n",
    "    for (m, c, tool), tm in tqdm(df_peaks.groupby(['modification', 'cell', 'tool'])):\n",
    "        reps = list(sorted(set(tm['replicate'])))\n",
    "        print(m, c, tool, len(reps))\n",
    "        files = [None] * len(reps)\n",
    "        peaks = [-1] * len(reps)\n",
    "        peaks_lens = [-1] * len(reps)\n",
    "        # print('Load')\n",
    "        for i in range(len(reps)):\n",
    "            t1 = tm[(tm['tool'] == tool) & (tm['replicate'] == reps[i])]\n",
    "            if len(t1) == 0:\n",
    "                continue\n",
    "            file = t1['file'].values[0]\n",
    "            peaks[i] = lines(file)\n",
    "            if peaks[i] == 0:\n",
    "                peaks_lens[i] = 0\n",
    "                continue\n",
    "            files[i] = pr.read_bed(sorted_file(file))\n",
    "            peaks_lens[i] = files[i].lengths().sum()\n",
    "            # print('Intersect')\n",
    "        for i, j in product(range(len(reps)), range(len(reps))):\n",
    "            if i >= j:\n",
    "                continue\n",
    "            rep1, rep2 = reps[i], reps[j]\n",
    "            file1, file2 = files[i], files[j]\n",
    "            peaks1, peaks2 = peaks[i], peaks[j]\n",
    "            peaks1_len, peaks2_len = peaks_lens[i], peaks_lens[j]\n",
    "            if peaks1 <= 0 or peaks2 <= 0:\n",
    "                reps_overlap.loc[len(reps_overlap)] = \\\n",
    "                    (m, c, tool, rep1, rep2, peaks1, peaks1_len, peaks2, peaks2_len,\n",
    "                     0, 0, 0)\n",
    "                continue\n",
    "            overlap1 = len(file1.overlap(file2))\n",
    "            overlap2 = len(file2.overlap(file1))\n",
    "            overlap_len = file1.intersect(file2).lengths().sum()\n",
    "            # assert file2.intersect(file1).lengths().sum() == overlap_len\n",
    "            reps_overlap.loc[len(reps_overlap)] = \\\n",
    "                (m, c, tool, rep1, rep2, peaks1, peaks1_len, peaks2, peaks2_len,\n",
    "                 overlap1, overlap2, overlap_len)\n",
    "    reps_overlap['jaccard'] = [\n",
    "        d(lo, l1 + l2 - lo)\n",
    "        for l1, l2, lo in zip(reps_overlap['peaks1_len'], reps_overlap['peaks2_len'], reps_overlap['peaks_overlap_len'])\n",
    "    ]\n",
    "    return reps_overlap"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Jaccard ENCODE "
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "t = df_encode_peaks.loc[[(m, c, r) not in OUTLIERS\n",
    "                         for _, (m, c, r) in df_encode_peaks[['modification', 'cell', 'replicate']].iterrows()]]\n",
    "# t = df_encode_peaks\n",
    "df_encode_reps = compute_reps_overlap(t)\n",
    "df_encode_reps.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=df_encode_reps, x='modification', hue='tool', y='jaccard',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_encode_reps, x='modification', hue='tool', y='jaccard',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.8,\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "ax.set_title('Jaccard between replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(14, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    g_results = sns.lineplot(data=df_encode_reps[df_encode_reps['modification'] == m],\n",
    "                             x='cell', y='jaccard', hue='tool',\n",
    "                             hue_order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                             palette=TOOLS_PALETTE,\n",
    "                             errorbar='se', marker='o', markersize=8,\n",
    "                             ax=ax)\n",
    "    g_results.axes.title.set_text(f'{m}')\n",
    "    ax.set_title(m)\n",
    "    if i != 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "t = pd.concat([df_encode_reps[['modification', 'cell', 'tool']],\n",
    "               df_encode_reps[['modification', 'cell', 'tool']]]).reset_index(drop=True)\n",
    "overlaps = [d(o, p) for o, p in zip(df_encode_reps['peaks1_overlap'], df_encode_reps['peaks1'])] + \\\n",
    "           [d(o, p) for o, p in zip(df_encode_reps['peaks2_overlap'], df_encode_reps['peaks2'])]\n",
    "t['overlap'] = overlaps\n",
    "print('ENCODE')\n",
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            showfliers=False\n",
    "            )\n",
    "sns.stripplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.8,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "\n",
    "ax.set_title('Overlap between replicates')\n",
    "ax.set_ylabel('Overlap')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "df_encode_reps['peaks_diff'] = np.fabs((df_encode_reps['peaks1'] - df_encode_reps['peaks2']) /\n",
    "                                       (df_encode_reps['peaks1'] + df_encode_reps['peaks2']) * 200)\n",
    "plt.figure(figsize=(14, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    g_results = sns.scatterplot(data=df_encode_reps[df_encode_reps['modification'] == m],\n",
    "                                x='peaks_diff', y='jaccard', hue='tool',\n",
    "                                hue_order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                                palette=TOOLS_PALETTE,\n",
    "                                ax=ax)\n",
    "    g_results.axes.title.set_text(f'{m}')\n",
    "    ax.set_title(m)\n",
    "    if i != 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "df_encode_reps['reps_peaks_diff'] = np.fabs((df_encode_reps['peaks1'] - df_encode_reps['peaks2']) /\n",
    "                                            (df_encode_reps['peaks1'] + df_encode_reps['peaks2']))\n",
    "\n",
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=df_encode_reps, x='modification', hue='tool', y='reps_peaks_diff',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_encode_reps, x='modification', hue='tool', y='reps_peaks_diff',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.8,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'], legend=False,\n",
    "              ax=ax)\n",
    "ax.set_title('Peaks difference between replicates')\n",
    "ax.set_ylabel('Peaks difference')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "rows = []\n",
    "for _, (modification, tool, jaccard, p1, p2, l1, l2) in df_encode_reps[\n",
    "    ['modification', 'tool', 'jaccard', 'peaks1', 'peaks2', 'peaks1_len', 'peaks2_len']].iterrows():\n",
    "    rows.append((modification, tool, jaccard, d(l1, p1)))\n",
    "    rows.append((modification, tool, jaccard, d(l2, p2)))\n",
    "df_jaccard_avg_len_encode = pd.DataFrame(rows, columns=['modification', 'tool', 'jaccard', 'average_len'])\n",
    "del rows"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(15, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.scatterplot(\n",
    "        df_jaccard_avg_len_encode[df_jaccard_avg_len_encode['modification'] == m],\n",
    "        palette=TOOLS_PALETTE,\n",
    "        hue_order=TOOLS,\n",
    "        x='jaccard', y='average_len', hue='tool', alpha=0.8,\n",
    "        ax=ax,\n",
    "        # fill=True\n",
    "    )\n",
    "    ax.set_title(m)\n",
    "    ax.set_xlim(-0.1, 1.1)\n",
    "    ax.set_title(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('average length')\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for m in ['H3K27ac', 'H3K4me1']:\n",
    "    print(m)\n",
    "    plt.figure(figsize=(5, 3))\n",
    "    ax = plt.axes()\n",
    "    sns.scatterplot(\n",
    "        df_jaccard_avg_len_encode[df_jaccard_avg_len_encode['modification'] == m],\n",
    "        palette=TOOLS_PALETTE,\n",
    "        hue_order=TOOLS,\n",
    "        x='jaccard', y='average_len', hue='tool', alpha=0.8,\n",
    "        ax=ax,\n",
    "        # fill=True\n",
    "    )\n",
    "    ax.set_title(m)\n",
    "    ax.set_xlim(-0.1, 1.1)\n",
    "    ax.set_title(m)\n",
    "    ax.set_ylabel('average length')\n",
    "    # Put a legend to the right of the current axis\n",
    "    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Jaccard Roadmap"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "t = df_immune_peaks.loc[[(m, c, r) not in OUTLIERS\n",
    "                         for _, (m, c, r) in df_immune_peaks[['modification', 'cell', 'replicate']].iterrows()]]\n",
    "# t = df_immune_peaks\n",
    "df_immune_reps = compute_reps_overlap(t)\n",
    "df_immune_reps.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Immune')\n",
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=df_immune_reps, x='modification', hue='tool', y='jaccard',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_immune_reps, x='modification', hue='tool', y='jaccard',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.8,\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              ax=ax)\n",
    "ax.set_title('Jaccard between replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(12, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    g_results = sns.lineplot(data=df_immune_reps[df_immune_reps['modification'] == m],\n",
    "                             x='cell', y='jaccard', hue='tool',\n",
    "                             hue_order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                             palette=TOOLS_PALETTE,\n",
    "                             errorbar='se', marker='o', markersize=8,\n",
    "                             ax=ax)\n",
    "    g_results.axes.title.set_text(f'{m}')\n",
    "    ax.set_title(m)\n",
    "    if i != 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    ax.xaxis.set_tick_params(rotation=90)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Immune')\n",
    "t = pd.concat([df_immune_reps[['modification', 'cell', 'tool']],\n",
    "               df_immune_reps[['modification', 'cell', 'tool']]]).reset_index(drop=True)\n",
    "overlaps = [d(o, p) for o, p in zip(df_immune_reps['peaks1_overlap'], df_immune_reps['peaks1'])] + \\\n",
    "           [d(o, p) for o, p in zip(df_immune_reps['peaks2_overlap'], df_immune_reps['peaks2'])]\n",
    "t['overlap'] = overlaps\n",
    "plt.figure(figsize=(8, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            ax=ax)\n",
    "sns.stripplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.8,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              ax=ax)\n",
    "\n",
    "ax.set_title('Overlap between replicates')\n",
    "ax.set_ylabel('Overlap')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "df_immune_reps['peaks_diff'] = np.fabs((df_immune_reps['peaks1'] - df_immune_reps['peaks2']) /\n",
    "                                       (df_immune_reps['peaks1'] + df_immune_reps['peaks2']) * 200)\n",
    "\n",
    "plt.figure(figsize=(14, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    print(m)\n",
    "    ax = axs[i]\n",
    "    g_results = sns.scatterplot(data=df_immune_reps[df_immune_reps['modification'] == m],\n",
    "                                x='peaks_diff', y='jaccard', hue='tool',\n",
    "                                hue_order=['MACS2', 'MACS2 broad', 'SICER', 'SPAN'],\n",
    "                                palette=TOOLS_PALETTE,\n",
    "                                ax=ax)\n",
    "    g_results.axes.title.set_text(f'{m}')\n",
    "    ax.set_title(m)\n",
    "    if i != 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "sns.set_style('whitegrid')\n",
    "print('Immune')\n",
    "df_immune_reps['peaks_diff'] = np.fabs((df_immune_reps['peaks1'] - df_immune_reps['peaks2']) /\n",
    "                                       (df_immune_reps['peaks1'] + df_immune_reps['peaks2']))\n",
    "\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=df_immune_reps, x='modification', hue='tool', y='peaks_diff',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_immune_reps, x='modification', hue='tool', y='peaks_diff',\n",
    "              dodge=True, size=3, palette='dark:black', alpha=0.8,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "ax.set_title('Peaks difference between replicates')\n",
    "ax.set_ylabel('Peaks difference')\n",
    "ax.set_xlabel('Modification')\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "rows = []\n",
    "for _, (modification, tool, jaccard, p1, p2, l1, l2) in df_immune_reps[\n",
    "    ['modification', 'tool', 'jaccard', 'peaks1', 'peaks2', 'peaks1_len', 'peaks2_len']].iterrows():\n",
    "    rows.append((modification, tool, jaccard, d(l1, p1)))\n",
    "    rows.append((modification, tool, jaccard, d(l2, p2)))\n",
    "df_jaccard_avg_len_immune = pd.DataFrame(rows, columns=['modification', 'tool', 'jaccard', 'average_len'])\n",
    "del rows"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(15, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.scatterplot(\n",
    "        df_jaccard_avg_len_immune[df_jaccard_avg_len_immune['modification'] == m],\n",
    "        palette=TOOLS_PALETTE,\n",
    "        hue_order=TOOLS,\n",
    "        x='jaccard', y='average_len', hue='tool', alpha=0.8,\n",
    "        ax=ax,\n",
    "    )\n",
    "    ax.set_title(m)\n",
    "    ax.set_xlim(-0.1, 1.1)\n",
    "    ax.set_title(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('average length')\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        #     Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Jaccard ABF"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_abf_reps = compute_reps_overlap(df_abf_peaks)\n",
    "df_abf_reps.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ABF')\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "# sns.boxplot(data=df_abf_reps, x='modification', hue='tool', y='jaccard',\n",
    "#             palette=TOOLS_PALETTE,\n",
    "#             order=MODIFICATIONS,\n",
    "#             hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "#             showfliers=False,\n",
    "#             ax=ax)\n",
    "sns.barplot(data=df_abf_reps, x='modification', hue='tool', y='jaccard',\n",
    "            order=MODIFICATIONS,\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_abf_reps, x='modification', hue='tool', y='jaccard',\n",
    "              dodge=True, size=1, palette='dark:black', alpha=0.2,\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "ax.set_title('Jaccard between ULI replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ABF')\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "# sns.boxplot(data=df_abf_reps, x='modification', hue='tool', y='jaccard',\n",
    "#             palette=TOOLS_PALETTE,\n",
    "#             order=MODIFICATIONS,\n",
    "#             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "#             showfliers=False,\n",
    "#             ax=ax)\n",
    "sns.barplot(data=df_abf_reps, x='modification', hue='tool', y='jaccard',\n",
    "            order=MODIFICATIONS,\n",
    "            palette=TOOLS_PALETTE,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "sns.stripplot(data=df_abf_reps, x='modification', hue='tool', y='jaccard',\n",
    "              dodge=True, size=1, palette='dark:black', alpha=0.2,\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "ax.set_title('Jaccard between ULI replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## All"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "sns.set_style(\"whitegrid\")\n",
    "t = pd.concat([df_encode_reps, df_immune_reps]).reset_index(drop=True)\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "# sns.boxplot(data=t, x='modification', hue='tool', y='jaccard',\n",
    "#             palette=TOOLS_PALETTE,\n",
    "#             order=MODIFICATIONS,\n",
    "#             hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "#             showfliers=False,\n",
    "#             ax=ax)\n",
    "sns.barplot(data=t, x='modification', hue='tool', y='jaccard',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "sns.stripplot(data=t, x='modification', hue='tool', y='jaccard',\n",
    "              dodge=True, size=2, palette='dark:black',\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False, alpha=0.4,\n",
    "              ax=ax)\n",
    "ax.set_title('Jaccard between replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "sns.set_style(\"whitegrid\")\n",
    "# t = pd.concat([df_encode_reps, df_immune_reps, df_abf_reps]).reset_index(drop=True)\n",
    "t = pd.concat([df_encode_reps, df_immune_reps]).reset_index(drop=True)\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "# sns.boxplot(data=t, x='modification', hue='tool', y='jaccard',\n",
    "#             palette=TOOLS_PALETTE,\n",
    "#             order=MODIFICATIONS,\n",
    "#             hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "#             showfliers=False,\n",
    "#             ax=ax)\n",
    "sns.barplot(data=t, x='modification', hue='tool', y='jaccard',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "            capsize=.2, err_kws={'linewidth': 2},\n",
    "            ax=ax)\n",
    "sns.stripplot(data=t, x='modification', hue='tool', y='jaccard',\n",
    "              dodge=True, size=2, palette='dark:black', alpha=0.4,\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "ax.set_title('Jaccard between replicates')\n",
    "ax.set_ylabel('Jaccard')\n",
    "ax.set_xlabel('Modification')\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "ts = []\n",
    "# for df in [df_encode_reps, df_immune_reps, df_abf_reps]:\n",
    "for df in [df_encode_reps, df_immune_reps]:\n",
    "    t = pd.concat([df[['modification', 'cell', 'tool']],\n",
    "                   df[['modification', 'cell', 'tool']]]).reset_index(drop=True)\n",
    "    overlaps = [d(o, p) for o, p in zip(df['peaks1_overlap'], df['peaks1'])] + \\\n",
    "               [d(o, p) for o, p in zip(df['peaks2_overlap'], df['peaks2'])]\n",
    "    t['overlap'] = overlaps\n",
    "    ts.append(t)\n",
    "t = pd.concat(ts).reset_index(drop=True)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('All')\n",
    "sns.set_style(\"whitegrid\")\n",
    "\n",
    "plt.figure(figsize=(6, 3))\n",
    "ax = plt.axes()\n",
    "sns.boxplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "            palette=TOOLS_PALETTE,\n",
    "            order=MODIFICATIONS,\n",
    "            hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "            showfliers=False,\n",
    "            ax=ax)\n",
    "sns.stripplot(data=t, x='modification', hue='tool', y='overlap',\n",
    "              dodge=True, size=0.2, palette='dark:black', alpha=0.8,\n",
    "              order=MODIFICATIONS,\n",
    "              hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "              legend=False,\n",
    "              ax=ax)\n",
    "\n",
    "ax.set_title('Overlap between replicates')\n",
    "ax.set_ylabel('Overlap')\n",
    "ax.set_xlabel('Modification')\n",
    "ax.xaxis.set_tick_params(rotation=45)\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print('ENCODE + Roadmap')\n",
    "sns.set_style(\"white\")\n",
    "# t = df_jaccard_avg_len_encode\n",
    "df_jaccard_avg_len_encode['dataset'] = 'ENCODE'\n",
    "df_jaccard_avg_len_immune['dataset'] = 'Roadmap'\n",
    "t = pd.concat([df_jaccard_avg_len_encode, df_jaccard_avg_len_immune]).reset_index(drop=True)\n",
    "plt.figure(figsize=(16, 3))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    sns.kdeplot(\n",
    "        t[t['modification'] == m],\n",
    "        palette=TOOLS_PALETTE,\n",
    "        hue_order=TOOLS,\n",
    "        x='jaccard', y='average_len', hue='tool', alpha=0.3,\n",
    "        common_norm=False, common_grid=False,\n",
    "        thresh=.05,\n",
    "        fill=True,\n",
    "        ax=ax,\n",
    "    )\n",
    "    sns.scatterplot(\n",
    "        t[t['modification'] == m],\n",
    "        palette=TOOLS_PALETTE,\n",
    "        hue_order=TOOLS,\n",
    "        style='dataset',\n",
    "        x='jaccard', y='average_len', hue='tool', alpha=0.8,\n",
    "        ax=ax,\n",
    "    )\n",
    "    # ax.set_title(m)\n",
    "    ax.set_xlim(-0.1, 1.1)\n",
    "    ax.set_title(m)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel(None)\n",
    "    else:\n",
    "        ax.set_ylabel('average length')\n",
    "    if i < len(axs) - 1:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Difference coverage"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import pyBigWig\n",
    "import tempfile\n",
    "\n",
    "tf = tempfile.mktemp()\n",
    "\n",
    "coverage_data = []\n",
    "\n",
    "all_bws_nc = all_bws[all_bws['modification'] != 'Control']\n",
    "all_bws_c = all_bws[all_bws['modification'] == 'Control']\n",
    "for _, (ds, m, c, r, bw_file) in tqdm(\n",
    "        list(all_bws_nc[['dataset', 'modification', 'cell', 'replicate', 'file']].iterrows())):\n",
    "    print(ds, m, c, r)\n",
    "    t = all_bws_c[(all_bws_c['dataset'] == ds) &\n",
    "                  (all_bws_c['cell'] == c) &\n",
    "                  (all_bws_c['replicate'] == r)]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    t = all_bws_c[(all_bws_c['dataset'] == ds) &\n",
    "                  (all_bws_c['cell'] == c)]\n",
    "    if len(t) == 0:\n",
    "        continue\n",
    "    bw_control_file = t['file'].values[0]\n",
    "    for tool in set(df_peaks['tool']):\n",
    "        t = df_peaks[(df_peaks['dataset'] == ds) &\n",
    "                     (df_peaks['modification'] == m) &\n",
    "                     (df_peaks['cell'] == c) &\n",
    "                     (df_peaks['replicate'] == r) &\n",
    "                     (df_peaks['tool'] == tool)]\n",
    "        if len(t) == 0:\n",
    "            continue\n",
    "        peaks_file = sorted_file(t['file'].values[0])\n",
    "        if lines(peaks_file) > 0:\n",
    "            t = pd.read_csv(peaks_file, sep='\\t', header=None)\n",
    "            # Speedup computations\n",
    "            t = t.sample(min(1_000, len(t)))\n",
    "            positions = [f'{chr}:{start}-{end}' for chr, start, end in zip(t[0], t[1], t[2]) if chr in CHROM_SIZES]\n",
    "            lengths = [end - start for chr, start, end in zip(t[0], t[1], t[2]) if chr in CHROM_SIZES]\n",
    "\n",
    "            with pyBigWig.open(bw_file) as bw:\n",
    "                total_coverage = sum(\n",
    "                    bw.stats(chr, type='sum', exact=True)[0] for chr in t[0].unique() if chr in CHROM_SIZES)\n",
    "                peaks_coverage = [bw.stats(chr, start, end, type='sum', exact=True)[0]\n",
    "                                  for chr, start, end in zip(t[0], t[1], t[2]) if chr in CHROM_SIZES]\n",
    "            for pos, length, peak_coverage in zip(positions, lengths, peaks_coverage):\n",
    "                coverage_data.append((ds, m, c, r, tool, pos, length, 'treatment', peak_coverage, total_coverage))\n",
    "\n",
    "            with pyBigWig.open(bw_control_file) as bwc:\n",
    "                total_coverage = sum(bwc.stats(chr, type='sum', exact=True)[0] for chr in t[0].unique())\n",
    "                peaks_coverage = [bwc.stats(chr, start, end, type='sum', exact=True)[0]\n",
    "                                  for chr, start, end in zip(t[0], t[1], t[2])]\n",
    "            for pos, length, peak_coverage in zip(positions, lengths, peaks_coverage):\n",
    "                coverage_data.append((ds, m, c, r, tool, pos, length, 'control', peak_coverage, total_coverage))\n",
    "        if tool == 'SPAN':\n",
    "            span_file = peaks_file\n",
    "\n",
    "    # Processing single tools information\n",
    "    for tool in TOOLS:\n",
    "        if tool == 'SPAN':\n",
    "            continue\n",
    "        t = df_peaks[(df_peaks['dataset'] == ds) &\n",
    "                     (df_peaks['modification'] == m) &\n",
    "                     (df_peaks['cell'] == c) &\n",
    "                     (df_peaks['replicate'] == r) &\n",
    "                     (df_peaks['tool'] == tool)]\n",
    "        if len(t) == 0:\n",
    "            continue\n",
    "        peaks_file = sorted_file(t['file'].values[0])\n",
    "        for name, args in [\n",
    "            (f'SPAN - {tool}', f' -a {span_file} -b {peaks_file} '),\n",
    "            (f'{tool} - SPAN', f' -b {span_file} -a {peaks_file} ')\n",
    "        ]:\n",
    "            ! bedtools intersect {args} -wa -v > {tf}\n",
    "            if lines(tf) > 0:\n",
    "                t = pd.read_csv(tf, sep='\\t', header=None)\n",
    "                # Speedup computations\n",
    "                t = t.sample(min(1_000, len(t)))\n",
    "                positions = [f'{chr}:{start}-{end}' for chr, start, end in zip(t[0], t[1], t[2]) if chr in CHROM_SIZES]\n",
    "                lengths = [end - start for chr, start, end in zip(t[0], t[1], t[2]) if chr in CHROM_SIZES]\n",
    "\n",
    "                with pyBigWig.open(bw_file) as bw:\n",
    "                    total_coverage = sum(\n",
    "                        bw.stats(chr, type='sum', exact=True)[0] for chr in t[0].unique() if chr in CHROM_SIZES)\n",
    "                    peaks_coverage = [bw.stats(chr, start, end, type='sum', exact=True)[0]\n",
    "                                      for chr, start, end in zip(t[0], t[1], t[2]) if chr in CHROM_SIZES]\n",
    "                for pos, length, peak_coverage in zip(positions, lengths, peaks_coverage):\n",
    "                    coverage_data.append((ds, m, c, r, name, pos, length, 'treatment', peak_coverage, total_coverage))\n",
    "\n",
    "                with pyBigWig.open(bw_control_file) as bwc:\n",
    "                    total_coverage = sum(bwc.stats(chr, type='sum', exact=True)[0] for chr in t[0].unique())\n",
    "                    peaks_coverage = [bwc.stats(chr, start, end, type='sum', exact=True)[0]\n",
    "                                      for chr, start, end in zip(t[0], t[1], t[2])]\n",
    "                for pos, length, peak_coverage in zip(positions, lengths, peaks_coverage):\n",
    "                    coverage_data.append((ds, m, c, r, name, pos, length, 'control', peak_coverage, total_coverage))\n",
    "\n",
    "df_coverage = pd.DataFrame(\n",
    "    coverage_data,\n",
    "    columns=['dataset', 'modification', 'cell', 'replicate', 'name', 'position', 'length', 'type', 'coverage',\n",
    "             'total_coverage']\n",
    ")\n",
    "del coverage_data\n",
    "df_coverage.sample(5)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "df_coverage['rpm'] = df_coverage['coverage'] / (df_coverage['total_coverage'] / 1e6)\n",
    "df_coverage['rpkm'] = df_coverage['rpm'] / (df_coverage['length'] / 1e3)\n",
    "\n",
    "for ds in ['ENCODE', 'Roadmap']:\n",
    "    for m in MODIFICATIONS:\n",
    "        plt.figure(figsize=(4.5, 4))\n",
    "        g_results = sns.boxplot(data=df_coverage[(df_coverage['dataset'] == ds) &\n",
    "                                                 (df_coverage['modification'] == m)], x='name', y='rpkm', hue='type',\n",
    "                                order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER',\n",
    "                                       'SPAN - MACS2', 'SPAN - MACS2 broad', 'SPAN - SICER',\n",
    "                                       'MACS2 - SPAN', 'MACS2 broad - SPAN', 'SICER - SPAN',\n",
    "                                       ],\n",
    "                                showfliers=False)\n",
    "        ax = g_results.axes\n",
    "        ax.xaxis.set_tick_params(rotation=90)\n",
    "        ax.set_ylabel('RPKM')\n",
    "        # ax.set_ylim(-0.8, 15)\n",
    "        # Put a legend to the right of the current axis\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(0.62, 0.82))\n",
    "        plt.title(f'Peaks coverage {ds} {m}')\n",
    "        plt.tight_layout()\n",
    "        # plt.savefig(f'{PATH}/pics/H3K4me1_coverage.pdf', bbox_inches='tight', dpi=300)\n",
    "        plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# Signal density"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "df_coverage",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "for ds in ['ENCODE', 'Roadmap']:\n",
    "    print(ds)\n",
    "    t = df_coverage[(df_coverage['dataset'] == ds) &\n",
    "                    # (df_coverage['name'].isin(['MACS2', 'MACS2 broad', 'SICER', 'SPAN'])\n",
    "                    (df_coverage['name'].isin(['SPAN']))\n",
    "                    ]\n",
    "    ts = []\n",
    "    for m, dfm in t.groupby(['modification']):\n",
    "        ts.append(dfm.sample(min(len(dfm), 5000)))\n",
    "    t = pd.concat(ts).reset_index(drop=True)\n",
    "    del ts\n",
    "    plt.figure(figsize=(6, 4))\n",
    "    ax = plt.axes()\n",
    "    t['rpkm'].clip(lower=0.2, inplace=True)\n",
    "    t['length'].clip(upper=3e5, inplace=True)\n",
    "    # ticks = [100, 1000, 10_000, 100_000]\n",
    "    # sns.kdeplot(\n",
    "    #     t, x='length', y='rpkm',\n",
    "    #     hue='modification',\n",
    "    #     alpha=0.8,\n",
    "    #     hue_order=MODIFICATIONS,\n",
    "    #     thresh=0.05,\n",
    "    # fill=True,\n",
    "    # ax=ax)\n",
    "    # \n",
    "    sns.scatterplot(data=t, x='length', y='rpkm', hue='modification', alpha=0.3,\n",
    "                    size=0.1,\n",
    "                    hue_order=MODIFICATIONS,\n",
    "                    ax=ax)\n",
    "    ax.set(xscale='log')\n",
    "    ax.set(yscale='log')\n",
    "    # ax.set(xticks=ticks)\n",
    "    # ax.set(xticklabels=ticks)\n",
    "    # Put a legend to the right of the current axis\n",
    "    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.tight_layout()\n",
    "    # plt.savefig(f'{GSE26320_PATH}/pics/peaks_length_vs_rpkm.pdf', bbox_inches='tight', dpi=300)\n",
    "    plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# IDR"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "def compute_reps_idr(df_peaks, modifications, cells):\n",
    "    reps_idr = pd.DataFrame(columns=['modification', 'cell', 'tool', 'rep1', 'rep2', 'peaks1', 'peaks2',\n",
    "                                     'idr_loaded_peaks1', 'idr_loaded_peaks2', 'idr_peaks', 'idr_loaded_peaks',\n",
    "                                     'idr_percentage'\n",
    "                                     ], dtype=object)\n",
    "\n",
    "    tf = tempfile.mktemp()\n",
    "    tf1 = tempfile.mktemp()\n",
    "    tf2 = tempfile.mktemp()\n",
    "\n",
    "    tools = list(sorted(set(df_peaks['tool'])))\n",
    "    for m, c in tqdm(list(product(modifications, cells))):\n",
    "        print(m, c)\n",
    "        tm = df_peaks[(df_peaks['modification'] == m) & (df_peaks['cell'] == c)]\n",
    "        reps = list(sorted(set(tm['replicate'])))\n",
    "        for tool in tools:\n",
    "            for i in range(len(reps)):\n",
    "                for j in range(i + 1, len(reps)):\n",
    "                    rep1, rep2 = reps[i], reps[j]\n",
    "                    t1 = tm[(tm['tool'] == tool) & (tm['replicate'] == rep1)]\n",
    "                    t2 = tm[(tm['tool'] == tool) & (tm['replicate'] == rep2)]\n",
    "                    if len(t1) == 0 or len(t2) == 0:\n",
    "                        continue\n",
    "                    try:\n",
    "                        file1 = t1['file'].values[0]\n",
    "                        file2 = t2['file'].values[0]\n",
    "                        print(tool, file1, file2)\n",
    "                        peaks1 = t1['peaks'].values[0]\n",
    "                        peaks2 = t2['peaks'].values[0]\n",
    "                        # Process SICER peaks, add missing 8th column, as -log10 qvalue\n",
    "                        if tool == 'SICER':\n",
    "                            sf1 = pd.read_csv(file1, sep='\\t', header=None)\n",
    "                            sf1[8] = -np.log10(sf1[7])\n",
    "                            sf1.to_csv(tf1, sep='\\t', index=None, header=False)\n",
    "                            file1 = tf1\n",
    "                            sf2 = pd.read_csv(file2, sep='\\t', header=None)\n",
    "                            sf2[8] = -np.log10(sf2[7])\n",
    "                            sf2.to_csv(tf2, sep='\\t', index=None, header=False)\n",
    "                            file2 = tf2\n",
    "                        ! idr --input-file-type bed --rank 8 --samples {file1} {file2} 2>&1 | tee {tf}\n",
    "                        with open(tf) as f:\n",
    "                            idr_out = f.read()\n",
    "                        # print('IDR OUT', tf, idr_out)\n",
    "                        lp1, lp2 = re.findall('Number of reported peaks - (\\\\d+)/(\\\\d+)', idr_out)[0]\n",
    "                        lp1, lp2 = int(lp1), int(lp2)\n",
    "                        ip, ilp, iperc = \\\n",
    "                            re.findall('Number of peaks passing IDR cutoff of 0.05 - (\\\\d+)/(\\\\d+) \\\\(([\\\\d\\\\.]+)',\n",
    "                                       idr_out)[0]\n",
    "                        ip, ilp, iperc = int(ip), int(ilp), float(iperc)\n",
    "                        print(m, c, tool, rep1, rep2, peaks1, peaks2, lp1, lp2, ip, ilp, iperc)\n",
    "                        reps_idr.loc[len(reps_idr)] = \\\n",
    "                            (m, c, tool, rep1, rep2, peaks1, peaks2, lp1, lp2, ip, ilp, iperc)\n",
    "                    except Exception as e:\n",
    "                        # raise e\n",
    "                        # print(e)\n",
    "                        pass  # Ignore errors for now\n",
    "    return reps_idr"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## IDR ENCODE"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# t = df_encode_peaks.loc[[(m, c, r) not in OUTLIERS for _, (m, c, r) in df_encode_peaks[['modification', 'cell', 'replicate']].iterrows()]]\n",
    "t = df_encode_peaks\n",
    "df_encode_idr_reps = compute_reps_idr(t, MODIFICATIONS, GSE26320_CELLS)\n",
    "df_encode_idr_reps.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('ENCODE')\n",
    "plt.figure(figsize=(8, 3))\n",
    "g_result = sns.boxplot(data=df_encode_idr_reps, x='modification', hue='tool', y='idr_percentage',\n",
    "                       palette=TOOLS_PALETTE,\n",
    "                       hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                       )\n",
    "ax = g_result.axes\n",
    "ax.set_title('IDR between replicates')\n",
    "ax.set_ylabel('IDR %')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "g_result = sns.boxplot(data=df_encode_idr_reps, x='modification', hue='tool', y='idr_peaks',\n",
    "                       palette=TOOLS_PALETTE,\n",
    "                       hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                       )\n",
    "ax = g_result.axes\n",
    "ax.set_title('IDR between replicates')\n",
    "ax.set_ylabel('IDR peaks')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "g_result = sns.boxplot(data=df_encode_idr_reps, x='modification', hue='tool', y='idr_loaded_peaks',\n",
    "                       palette=TOOLS_PALETTE,\n",
    "                       hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                       )\n",
    "ax = g_result.axes\n",
    "ax.set_title('IDR between replicates')\n",
    "ax.set_ylabel('IDR loaded peaks')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for m in MODIFICATIONS:\n",
    "    plt.figure(figsize=(6, 2))\n",
    "    g_result = sns.boxplot(data=df_encode_idr_reps[df_encode_idr_reps['modification'] == m],\n",
    "                           x='idr_peaks', y='tool',\n",
    "                           palette=TOOLS_PALETTE,\n",
    "                           order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                           )\n",
    "    ax = g_result.axes\n",
    "    ax.set_title(f'IDR between {m} replicates')\n",
    "    ax.set_ylabel('Tool')\n",
    "    ax.set_xlabel('IDR peaks')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(10, 2))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    t = df_encode_idr_reps[df_encode_idr_reps['modification'] == m]\n",
    "    g_results = sns.scatterplot(data=t,\n",
    "                                x='idr_peaks', y='idr_percentage',\n",
    "                                hue='tool',\n",
    "                                palette=TOOLS_PALETTE,\n",
    "                                alpha=0.8,\n",
    "                                ax=ax,\n",
    "                                )\n",
    "    for x, y, cell in zip(t['idr_peaks'], t['idr_percentage'], t['cell']):\n",
    "        ax.text(x, y, cell, fontsize=4)\n",
    "    ax.title.set_text(m)\n",
    "    if i != 0:\n",
    "        ax.set_ylabel(None)\n",
    "    if i != 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## IDR Roadmap"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# t = df_immune_peaks.loc[[(m, c, r) not in OUTLIERS for _, (m, c, r) in df_immune_peaks[['modification', 'cell', 'replicate']].iterrows()]]\n",
    "t = df_immune_peaks\n",
    "df_immune_idr_reps = compute_reps_idr(t, MODIFICATIONS, IMMUNE_CELLS)\n",
    "df_immune_idr_reps.sample(5)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Immune')\n",
    "plt.figure(figsize=(8, 3))\n",
    "g_result = sns.boxplot(data=df_immune_idr_reps, x='modification', hue='tool', y='idr_percentage',\n",
    "                       palette=TOOLS_PALETTE,\n",
    "                       hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                       )\n",
    "ax = g_result.axes\n",
    "ax.set_title('IDR between replicates')\n",
    "ax.set_ylabel('IDR %')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "print('Immune')\n",
    "plt.figure(figsize=(8, 3))\n",
    "g_result = sns.boxplot(data=df_immune_idr_reps, x='modification', hue='tool', y='idr_peaks',\n",
    "                       palette=TOOLS_PALETTE,\n",
    "                       hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                       )\n",
    "ax = g_result.axes\n",
    "ax.set_title('IDR between replicates')\n",
    "ax.set_ylabel('IDR peaks')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "g_result = sns.boxplot(data=df_immune_idr_reps, x='modification', hue='tool', y='idr_loaded_peaks',\n",
    "                       palette=TOOLS_PALETTE,\n",
    "                       hue_order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                       )\n",
    "ax = g_result.axes\n",
    "ax.set_title('IDR between replicates')\n",
    "ax.set_ylabel('IDR loaded peaks')\n",
    "ax.set_xlabel('Modification')\n",
    "sns.move_legend(ax, loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for m in MODIFICATIONS:\n",
    "    plt.figure(figsize=(6, 2))\n",
    "    g_result = sns.boxplot(data=df_immune_idr_reps[df_immune_idr_reps['modification'] == m],\n",
    "                           x='idr_peaks', y='tool',\n",
    "                           palette=TOOLS_PALETTE,\n",
    "                           order=['MACS2', 'MACS2 broad', 'SPAN', 'SICER'],\n",
    "                           )\n",
    "    ax = g_result.axes\n",
    "    ax.set_title(f'IDR between {m} replicates')\n",
    "    ax.set_ylabel('Tool')\n",
    "    ax.set_xlabel('IDR peaks')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "plt.figure(figsize=(10, 2))\n",
    "axs = [plt.subplot(1, 5, i + 1) for i in range(5)]\n",
    "for i, m in enumerate(MODIFICATIONS):\n",
    "    ax = axs[i]\n",
    "    t = df_immune_idr_reps[df_immune_idr_reps['modification'] == m]\n",
    "    g_results = sns.scatterplot(data=t,\n",
    "                                x='idr_peaks', y='idr_percentage',\n",
    "                                hue='tool',\n",
    "                                palette=TOOLS_PALETTE,\n",
    "                                alpha=0.8,\n",
    "                                ax=ax,\n",
    "                                )\n",
    "    for x, y, cell in zip(t['idr_peaks'], t['idr_percentage'], t['cell']):\n",
    "        ax.text(x, y, cell, fontsize=4)\n",
    "    ax.title.set_text(m)\n",
    "    if i != 0:\n",
    "        ax.set_ylabel(None)\n",
    "    if i != 4:\n",
    "        ax.legend().set_visible(False)\n",
    "    else:\n",
    "        # Put a legend to the right of the current axis\n",
    "        g_results.axes.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}